0 引言

裂缝通常会引起地震反射异常^[1]，增加岩石有效孔隙度，提高渗透率^[2-3]。裂缝密度是表征储层裂缝的关键参数之一^[4]。为准确获取裂缝密度参数，学者们利用方位各向异性地震资料开展了大量研究^[5-9]。Bakulin等^[10]认为当裂缝分别充填油、水、气三种不同类型的流体时，Schoenberg线性滑动模型参数和Thomsen各向异性参数均与裂缝密度存在定量关系。随后，Zhang等^[11]、Downton等^[12]实现了裂缝岩石物理参数与地震反射系数之间的相互转换，为裂缝密度定量预测奠定了基础。Isabel^[13]首次在横向各向同性(Horizontally Transverse Isotropic，HTI)介质的裂缝密度反演过程中添加了奇异值分解(Singular Value Decomposition，SVD)方法，但该方法中各向异性对反射系数的贡献度并不高。于是，刘财等^[14]运用不同方位相交测线反射系数取差的办法改善了这一缺点。Li等^[4]提出了一种基于SVD的受光滑先验弹性参数约束的裂缝密度和裂缝方位反演方法。由于各向异性梯度可作为裂缝密度的指示因子^[15]，Al-marzoug等^[16]将多方位振幅随炮检距(Amplitude Versus Offset，AVO)的梯度拟合成椭圆，并使用拟合椭圆的主轴方向和椭圆率分别表征裂缝方位和裂缝密度的分布特征。秦海旭等^[17]推导了各向异性梯度与裂缝密度之间的关系，通过各向异性梯度间接反演获得裂缝密度。为充分利用炮检距向量片(Offset Vector Tile，OVT)域地震数据包含的叠前各向异性信息，魏欣伟等^[18]提出了一种叠前AVAZ（振幅随炮检距和方位角变化）的裂缝密度预测方法。针对弹性波阻抗^[19]，Martin等^[20]指出弱各向异性介质中的弹性阻抗为各向同性弹性阻抗背景部分与各向异性弹性阻抗扰动部分的乘积。陈超等^[21]将组稀疏引入到方位弹性阻抗反演目标函数中，利用反演获得的裂缝弱度参数计算裂缝密度。

直接反演裂缝密度的方法比间接换算更具优势，可以有效地避免累积误差，提高反演结果的精度。根据这一思想，本文基于具有水平对称轴的HTI介质推导出具有裂缝密度表示的反射系数公式，并开展了裂缝密度参数的直接反演。同时，为了进一步提高反演结果的横向分辨率和稳定性，在ADMM(Alternate Direction Multiplier Method，交替方向乘子法)^[22]中引入了反距离加权插值算法^[23]的思想，对横向差分算子进行修改以增强反演方法的空间约束能力，进而提高反演目标的横向连续性。最后，通过理论数据测试和实际资料应用验证了本文提出的裂缝密度预测方法(简称“本文方法”)的有效性。

1 基本原理 1.1 含裂缝密度的反射系数近似方程

在长波长假设条件下，线性滑动模型与Hudson薄硬币状裂隙模型具有相同的一阶表示形式，因此根据Hudson模型一阶等效弹性刚度矩阵可建立线性滑动模型参数与裂缝岩石相关参数之间的关系^[24-25]

$ {\delta }_{\mathrm{N}}=\frac{4e}{3g\left(1-g\right)\left[1+\frac{1}{\mathrm{\pi }\left(1-g\right)}\left(\frac{k\mathrm{\text{'}}+\frac{4}{3}\mu \mathrm{\text{'}}}{\mu \chi }\right)\right]} $

(1)

$ {\delta }_{\mathrm{T}}=\frac{16e}{3\left(3-2g\right)\left[1+\frac{4}{\mathrm{\pi }\left(3-2g\right)}\left(\frac{\mu \mathrm{\text{'}}}{\mu \chi }\right)\right]} $

(2)

式中：$ {\delta }_{\mathrm{N}} $、$ {\delta }_{\mathrm{T}} $分别为裂缝法向弱度、切向弱度；e为裂缝密度；$ g={{V}_{\mathrm{S}}}^{2}/{{V}_{\mathrm{P}}}^{2}= $$ \mu /(\lambda +2\mu ) $为横纵波速度比的平方，其中$ {V}_{\mathrm{P}} $、$ {V}_{\mathrm{S}} $分别为围岩的纵、横波速度，$ \mu $、$ \lambda $均为不含裂隙各向同性岩石的拉梅参数; $ \chi $为裂缝纵横比；$ k\mathrm{\text{'}} $、$ \mu \mathrm{\text{'}} $分别是裂缝充填物的体积模量、剪切模量。

从式(2)可以看出，裂缝弱度与横纵波速度比、裂缝密度、裂缝填充物有关。当裂缝为干裂缝或充填气体时，可以将裂缝充填物的剪切模量和体积模量均设置为0(即$ k\mathrm{\text{'}}=\mu \mathrm{\text{'}}=0 $)，此时

$ \left\{\begin{array}{c}{\delta }_{\mathrm{N}}=\frac{4e}{3g\left(1-g\right)}\\ {\delta }_{\mathrm{T}}=\frac{16e}{3\left(3-2g\right)}\end{array}\right. $

(3)

由于流体剪切模量为0并且体积模量较小，因此当裂缝中饱含流体时，对于纵横比小的裂缝$ \chi \ll 1 $，$ (k\mathrm{\text{'}}+4/3\mu \mathrm{\text{'}})/\mu \chi \gg 1 $。此时有

$ \left\{\begin{array}{l}{\delta }_{\mathrm{N}}=0\\ {\delta }_{\mathrm{T}}=\frac{16e}{3\left(3-2g\right)}\end{array}\right. $

(4)

可见法向裂缝弱度受到裂缝中流体的影响，而切向裂缝弱度不论裂缝充填流体与否，都保持不变。定义裂缝流体因子为

$ f=\frac{1}{\left[1-g+\left(\frac{k\mathrm{\text{'}}+\frac{4}{3}\mu \mathrm{\text{'}}}{\mu \chi \mathrm{\pi }}\right)\right]} $

(5)

由于水的体积模量比油气大，因此当裂缝中含油气量增加时，裂缝流体因子增大。此时，法向裂缝弱度可表示为裂缝密度与裂缝流体因子的乘积，即

$ \begin{array}{l}{\delta }_{\mathrm{N}}=\frac{4}{3g}\frac{e}{1-g+\frac{k\mathrm{\text{'}}+\frac{4}{3}\mu \mathrm{\text{'}}}{\mu \chi \mathrm{\pi }}}\\ \;\;\;\;=\frac{4}{3g}ef=\frac{4}{3g}F\end{array} $

(6)

式中$ F $为可采性因子，即地质甜点(裂缝流体因子$ f $)和工程甜点(裂缝密度$ e $)的乘积($ F=f\times e $)。含裂缝碳酸盐岩可以通过HTI介质进行等效。在弱各向异性假设条件下，该等效介质分界面PP波反射系数由两部分组成^[26-29]。

$ {R}_{\mathrm{P}\mathrm{P}}\left(\theta ;\varphi \right)={R}_{\mathrm{P}\mathrm{P}}^{\mathrm{i}\mathrm{s}\mathrm{o}}\left(\theta \right)+{R}_{\mathrm{P}\mathrm{P}}^{\mathrm{a}\mathrm{n}\mathrm{i}}\left(\theta ;\varphi \right) $

(7)

$ \begin{array}{l}{R}_{\mathrm{P}\mathrm{P}}^{\mathrm{i}\mathrm{s}\mathrm{o}}\left(\theta \right)=\frac{1}{2}\mathrm{s}\mathrm{e}{\mathrm{c}}^{2}\stackrel{-}{\theta }\frac{\mathrm{\Delta }{V}_{\mathrm{P}}}{\overline{{V}_{\mathrm{P}}}}-4g\mathrm{s}\mathrm{i}\mathrm{n}\theta \frac{\mathrm{\Delta }{V}_{\mathrm{S}}}{\overline{{V}_{\mathrm{S}}}}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}\left(1-4g\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \right)\frac{\mathrm{\Delta }\rho }{\stackrel{-}{\rho }}\end{array} $

(8)

$ \begin{array}{l}{R}_{\mathrm{P}\mathrm{P}}^{\mathrm{a}\mathrm{n}\mathrm{i}}\left(\theta ;\varphi \right)=-\left[g\left(1-2g\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\varphi \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \right.\left(1+\mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\theta \right)+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left.{g}^{2}\mathrm{c}\mathrm{o}{\mathrm{s}}^{4}\varphi \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\theta \right]\mathrm{\Delta }{\delta }_{\mathrm{N}}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;g\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\varphi \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \left(1-\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\varphi \mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\theta \right)\mathrm{\Delta }{\delta }_{\mathrm{T}}\left(9\right)\end{array} $

式中：$ {R}_{\mathrm{P}\mathrm{P}} $为纵波反射系数；$ \theta $为入射角；$ \varphi $为裂缝法向与测线方向的夹角；$ \rho $为密度；上标iso、ani分别表示各向同性部分、各向异性部分；$ \mathrm{\Delta } $表示界面上、下地层参数之间的变化量；顶横线“-”表示界面上、下地层参数的平均值。

当$ \theta < {30}^{\circ } $时，可以忽略与$ \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \mathrm{t}\mathrm{a}{\mathrm{n}}^{2}\theta $相关的项，可得

$ {R}_{\mathrm{P}\mathrm{P}}^{\mathrm{a}\mathrm{n}\mathrm{i}}\left(\theta ;\varphi \right)=-g\left(1-2g\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\varphi \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \mathrm{\Delta }{\delta }_{\mathrm{N}}+\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;g\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\varphi \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \mathrm{\Delta }{\delta }_{\mathrm{T}} $

(10)

将式(3)、式(4)代入到式(10)，可得

$ \begin{array}{l}{R}_{\mathrm{P}\mathrm{P}}^{\mathrm{a}\mathrm{n}\mathrm{i}}\left(\theta ;\varphi \right)=-\frac{4}{3}\left(1-2g\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\varphi \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \mathrm{\Delta }F+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{16g}{3\left(3-2g\right)}\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\varphi \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta \mathrm{\Delta }e\end{array} $

(11)

1.2 基于空间约束的ADMM裂缝密度AVAZ反演

利用两个方位($ {\varphi }_{1} $和$ {\varphi }_{2} $)的地震数据取差，消除各向同性背景部分后的纵波反射系数可以表示为

$ \mathrm{\Delta }{R}_{\mathrm{P}\mathrm{P}}(\theta ;{\varphi }_{1},{\varphi }_{2})={A}_{1}(\theta ;{\varphi }_{1},{\varphi }_{2})\mathrm{\Delta }F+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{A}_{2}(\theta ;{\varphi }_{1},{\varphi }_{2})\mathrm{\Delta }e $

(12)

其中

$ \left\{\begin{array}{l}{A}_{1}(\theta ;{\varphi }_{1},{\varphi }_{2})=-\frac{4}{3}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta (1-2g)(\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}{\varphi }_{2}-\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}{\varphi }_{1})\\ {A}_{2}(\theta ;{\varphi }_{1},{\varphi }_{2})=\frac{16g}{3\left(3-2g\right)}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\theta (\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}{\varphi }_{2}-\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}{\varphi }_{1})\end{array}\right. $

(13)

式中$ \mathrm{\Delta }{R}_{\mathrm{P}\mathrm{P}} $为各向异性扰动部分的纵波反射系数。

方位差数据可以表示为

$ \mathrm{\Delta }\boldsymbol{S}=\boldsymbol{G}\boldsymbol{L} $

(14)

其中

$ \left\{\begin{array}{l}\mathrm{\Delta }\boldsymbol{S}={\left[\begin{array}{c}\begin{array}{c}{\boldsymbol{S}}_{\mathrm{P}\mathrm{P}}^{}({\theta }_{1};{\varphi }_{2})-{\boldsymbol{S}}_{\mathrm{P}\mathrm{P}}^{}({\theta }_{1};{\varphi }_{1})\\ {\boldsymbol{S}}_{\mathrm{P}\mathrm{P}}^{}({\theta }_{2};{\varphi }_{2})-{\boldsymbol{S}}_{\mathrm{P}\mathrm{P}}^{}({\theta }_{2};{\varphi }_{1})\end{array}\\ \begin{array}{c}⋮\\ {\boldsymbol{S}}_{\mathrm{P}\mathrm{P}}^{}({\theta }_{m};{\varphi }_{2})-{\boldsymbol{S}}_{\mathrm{P}\mathrm{P}}^{}({\theta }_{m};{\varphi }_{1})\end{array}\end{array}\right]}_{mn\times t}\\ \boldsymbol{L}={\left[\begin{array}{c}\mathrm{\Delta }F\\ \mathrm{\Delta }e\end{array}\right]}_{2n\times t}\\ {\boldsymbol{S}}_{\mathrm{P}\mathrm{P}}^{}={\left[\begin{array}{c}{S}_{{}_{\mathrm{P}\mathrm{P}}}^{1}\\ ⋮\\ {S}_{{}_{\mathrm{P}\mathrm{P}}}^{n}\end{array}\right]}_{n\times t}\end{array}\right. $

(15)

式中：$ \mathrm{\Delta }\boldsymbol{S} $为方位差数据；$ \boldsymbol{L} $为代反演参数；$ \boldsymbol{G} $为正演算子；m为入射角个数；n为采样点个数；t为反演剖面道数。正演算子可以表示为

$ \boldsymbol{G}={\left[\begin{array}{cc}{\boldsymbol{W}}_{0}{\boldsymbol{B}}_{1}\left({\theta }_{1}\right)& {\boldsymbol{W}}_{0}{\boldsymbol{B}}_{2}\left({\theta }_{1}\right)\\ {\boldsymbol{W}}_{0}{\boldsymbol{B}}_{1}\left({\theta }_{2}\right)& {\boldsymbol{W}}_{0}{\boldsymbol{B}}_{2}\left({\theta }_{2}\right)\\ ⋮& ⋮\\ {\boldsymbol{W}}_{0}{\boldsymbol{B}}_{1}\left({\theta }_{m}\right)& {\boldsymbol{W}}_{0}{\boldsymbol{B}}_{2}\left({\theta }_{m}\right)\end{array}\right]}_{mn\times 2n} $

(16)

其中

$ \left\{\begin{array}{l}{\boldsymbol{W}}_{0}^{}={\left[\begin{array}{cccc}{w}_{1}& & & \\ ⋮& {w}_{1}& & \\ {w}_{k}& ⋮& \ddots & \\ & {w}_{k}& \ddots & {w}_{1}\\ & & \ddots & ⋮\\ & & & {w}_{k}\end{array}\right]}_{(n+k-1)\times n}\\ {\boldsymbol{B}}_{1}\left({\theta }_{i}\right)={\left[\begin{array}{ccc}{A}_{1}({\theta }_{i};{\varphi }_{1},{\varphi }_{2})& & \\ & \ddots & \\ & & {A}_{1}({\theta }_{i};{\varphi }_{1},{\varphi }_{2})\end{array}\right]}_{n\times n}\\ {\boldsymbol{B}}_{2}\left({\theta }_{i}\right)={\left[\begin{array}{ccc}{A}_{2}({\theta }_{i};{\varphi }_{1},{\varphi }_{2})& & \\ & \ddots & \\ & & {A}_{2}({\theta }_{i};{\varphi }_{1},{\varphi }_{2})\end{array}\right]}_{n\times n}\\ i=\mathrm{1,2},\cdots ,m\end{array}\right. $

(17)

式中$ {w}_{k} $为子波元素, k为子波的长度。为了简化表示，将$ \boldsymbol{B} $记为

$ \boldsymbol{B}={\left[\begin{array}{cc}{\boldsymbol{B}}_{1}\left({\theta }_{1}\right)& {\boldsymbol{B}}_{2}\left({\theta }_{1}\right)\\ {\boldsymbol{B}}_{1}\left({\theta }_{2}\right)& {\boldsymbol{B}}_{2}\left({\theta }_{2}\right)\\ ⋮& ⋮\\ {\boldsymbol{B}}_{1}\left({\theta }_{m}\right)& {\boldsymbol{B}}_{2}\left({\theta }_{m}\right)\end{array}\right]}_{mn\times 2n} $

(18)

并引入纵向差分矩阵$ {\boldsymbol{D}}_{y} $，则将$ \boldsymbol{L} $可表示为

$ \boldsymbol{L}={\boldsymbol{D}}_{y}\boldsymbol{M} $

(19)

其中

$ \left\{\begin{array}{l}{\boldsymbol{D}}_{y}={\left[\begin{array}{ccc}{\boldsymbol{D}}_{1}& & \\ & \ddots & \\ & & {\boldsymbol{D}}_{1}\end{array}\right]}_{mn\times 2n}\\ {\boldsymbol{D}}_{1}={\left[\begin{array}{ccccc}-1& 1& & & \\ & -1& 1& & \\ & & \ddots & \ddots & \\ & & & -1& 1\end{array}\right]}_{n\times n}\\ \boldsymbol{M}={\left[\begin{array}{cccc}{F}_{11}& {F}_{21}& \cdots & {F}_{t1}\\ & & ⋮& \\ {F}_{1n}& {F}_{1n}& \cdots & {F}_{tn}\\ {e}_{11}& {e}_{21}& \cdots & {e}_{t1}\\ & & ⋮& \\ {e}_{1n}& {e}_{2n}& \cdots & {e}_{tn}\end{array}\right]}_{2n\times t}\end{array}\right. $

(20)

式中：D为差分矩阵；M为待反演参数构成的矩阵。

为了提高多个参数反演的横向稳定性，ATV(Anisotropic Total Variation)算法引入横向差分$ {\boldsymbol{D}}_{x} $算子，实现多个地震道的同时反演

$ \mathrm{A}\mathrm{T}\mathrm{V}\left(\boldsymbol{M}\right)={‖\boldsymbol{M}{\boldsymbol{D}}_{x}‖}_{1}+{‖{\boldsymbol{D}}_{y}\boldsymbol{M}‖}_{1} $

(21)

式中$ {‖\cdot ‖}_{1} $为$ {\mathrm{L}}_{1} $范数。在进行多道反演时，仅使用来自两个相邻道的差分约束反演结果可能会出现较大的偏差。因此，可以通过将当前道与多个相邻道之间的差，利用反距离加权(Inverse Distance Weighted，IDW)算法^[23]提高反演结果的横向连续性。

反距离加权插值法基于相近、相似原理，认为两个物体距离越近，性质就越相似，因而以插值点与样本点间的距离为权重进行加权平均，离插值点越近的样本点赋予的权重越大。反距离加权插值为

$ \widehat{Z}\left({h}_{0}\right)=\sum\limits _{j=1}^{N}{\vartheta }_{j}Z\left({h}_{j}\right) $

(22)

式中：$ \widehat{Z}\left({h}_{0}\right) $为预测点$ {h}_{0} $处的采用反距离加权插值的预测值；N为预测点周围的样点数；$ {\vartheta }_{j} $为第j道地震数据的权重值，可以通过不同点之间的距离来计算，即

$ {\vartheta }_{j}=\frac{{d}_{j,0}^{-p}}{\sum\limits _{j=1}^{N}{d}_{j,0}^{-p}} $

(23)

式中：p为指数值，一般取2；$ {d}_{j,0} $为预测点$ {h}_{0} $与已知样点$ {h}_{j} $之间的距离。

根据式(22)，利用邻近两道地震数据和初始模型进行约束，可以得到反演目标函数为

$ \begin{array}{l}J\left(\boldsymbol{M}\right)=\underset{\boldsymbol{M}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left[\right||\boldsymbol{W}\boldsymbol{B}{\boldsymbol{D}}_{y}\boldsymbol{M}-\mathrm{\Delta }\boldsymbol{S}|{|}_{2}^{2}+\kappa \left|\right|\boldsymbol{M}-{\boldsymbol{M}}_{0}|{|}_{2}^{2}+\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \alpha \left(\left|\right|\boldsymbol{M}{\boldsymbol{D}}_{x}|{|}_{1}+|\left|{\boldsymbol{D}}_{y}\boldsymbol{M}\right|{|}_{1}\right)]\end{array} $

(24)

其中

$ {\boldsymbol{D}}_{x}={{\left[\begin{array}{cccc}-1& & & \\ {\vartheta }_{12}& -1& & \\ {\vartheta }_{13}& \lambda {\vartheta }_{23}& \ddots & \\ & {\vartheta }_{24}& \ddots & -1\\ & & \ddots & {\vartheta }_{\left(n-2\right)\left(n-1\right)}\\ & & & {\vartheta }_{\left(n-2\right)n}\end{array}\right]}_{t\times }}_{t} $

(25)

式中：$ \boldsymbol{W}={\left[\begin{array}{ccc}{\boldsymbol{W}}_{0}& & \\ & \ddots & \\ & & {\boldsymbol{W}}_{0}\end{array}\right]}_{mn\times mn} $为子波矩阵；B为构成反射系数的角度系数矩阵；$ {\boldsymbol{M}}_{0} $为低频模型；$ \kappa $为低频约束项的系数；$ \alpha $为稀疏约束项的系数。

本文采用ADMM^[30]对含多个未知参数的待解目标函数(式(24)）求解，将目标函数分解成多个单一参数的子问题。首先引入矩阵$ {\boldsymbol{M}}_{x} $、$ {\boldsymbol{M}}_{y} $分别代替$ \boldsymbol{M}{\boldsymbol{D}}_{x} $、$ {\boldsymbol{D}}_{y}\boldsymbol{M} $，可以将式(24)简化为

$ \begin{array}{l}J\left(\boldsymbol{M}\right)=\underset{\boldsymbol{M}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left[\right||\boldsymbol{W}\boldsymbol{B}{\boldsymbol{D}}_{y}\boldsymbol{M}-\mathrm{\Delta }\boldsymbol{S}|{|}_{2}^{2}+\kappa \left|\right|\boldsymbol{M}-{\boldsymbol{M}}_{0}|{|}_{2}^{2}+\\ \;\;\;\;\;\;\;\;\;\;\;\alpha \left(\left|\right|{\boldsymbol{M}}_{x}|{|}_{1}+|\left|{\boldsymbol{M}}_{y}\right|{|}_{1}\right)]\mathrm{s}.\mathrm{t}.{\boldsymbol{M}}_{x}=\boldsymbol{M}{\boldsymbol{D}}_{\mathrm{x}}{\boldsymbol{M}}_{y}={\boldsymbol{D}}_{y}\boldsymbol{M}\end{array} $

(26)

再引入对偶项$ {\boldsymbol{C}}_{x} $、$ {\boldsymbol{C}}_{y} $，得到无约束优化问题

$ \begin{array}{l}J\left(\boldsymbol{M}\right)=\underset{\boldsymbol{M}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left[\right||\boldsymbol{W}\boldsymbol{B}{\boldsymbol{D}}_{y}\boldsymbol{M}-\mathrm{\Delta }\boldsymbol{S}|{|}_{2}^{2}+\kappa \left|\right|\boldsymbol{M}-{\boldsymbol{M}}_{0}|{|}_{2}^{2}+\\\;\;\;\;\;\;\;\;\;\;\; \alpha \left(\left|\right|{\boldsymbol{M}}_{x}|{|}_{1}+|\left|{\boldsymbol{M}}_{y}\right|{|}_{1}\right)+\eta {‖{\boldsymbol{M}}_{x}-\boldsymbol{M}{\boldsymbol{D}}_{x}-{\boldsymbol{C}}_{x}‖}_{2}^{2}+\\ \;\;\;\;\;\;\;\;\;\;\;\eta {‖{\boldsymbol{M}}_{y}-{\boldsymbol{D}}_{y}\boldsymbol{M}-{\boldsymbol{C}}_{y}‖}_{2}^{2}\left]\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\right(27)\end{array} $

式中$ \eta $为差分的拉格朗日约束参数。

于是，将式(27)拆解成分别关于$ \boldsymbol{M} $、$ {\boldsymbol{M}}_{x} $、$ {\boldsymbol{M}}_{y} $、和$ {\boldsymbol{C}}_{x} $、$ {\boldsymbol{C}}_{y} $的子优化问题，其中$ \boldsymbol{M} $的子目标函数为

$ \begin{array}{l}{K}_{1}\left(\boldsymbol{M}\right)=\underset{\boldsymbol{M}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left(\right||\boldsymbol{W}\boldsymbol{B}{\boldsymbol{D}}_{y}\boldsymbol{M}-\mathrm{\Delta }\boldsymbol{S}|{|}_{2}^{2}+\kappa \left|\right|\boldsymbol{M}-{\boldsymbol{M}}_{0}|{|}_{2}^{2}+\\\;\;\;\;\;\;\;\;\;\;\;\;\; \eta {‖{\boldsymbol{M}}_{x}-\boldsymbol{M}{\boldsymbol{D}}_{x}-{\boldsymbol{C}}_{x}‖}_{2}^{2}+\eta {‖{\boldsymbol{M}}_{y}-{\boldsymbol{D}}_{y}\boldsymbol{M}-{\boldsymbol{C}}_{y}‖}_{2}^{2})\end{array} $

(28)

$ {\boldsymbol{M}}_{x} $、$ {\boldsymbol{M}}_{y} $的子目标函数分别为

$ {K}_{2}\left({\boldsymbol{M}}_{x}\right)=\underset{{\boldsymbol{M}}_{x}}{\mathrm{m}\mathrm{i}\mathrm{n}}(\alpha {‖{\boldsymbol{M}}_{x}‖}_{1}+\eta {‖{\boldsymbol{M}}_{x}-\boldsymbol{M}{\boldsymbol{D}}_{x}-{\boldsymbol{C}}_{x}‖}_{2}^{2}) $

(29)

$ {K}_{3}\left({\boldsymbol{M}}_{y}\right)=\underset{{\boldsymbol{M}}_{y}}{\mathrm{m}\mathrm{i}\mathrm{n}}(\alpha {‖{\boldsymbol{M}}_{y}‖}_{1}+\eta {‖{\boldsymbol{M}}_{y}-{\boldsymbol{D}}_{y}\boldsymbol{M}-{\boldsymbol{C}}_{y}‖}_{2}^{2}) $

(30)

$ {\boldsymbol{C}}_{x} $、$ {\boldsymbol{C}}_{y} $的子目标函数分别为

$ {K}_{4}\left({\boldsymbol{C}}_{x}\right)=\underset{{\boldsymbol{C}}_{x}}{\mathrm{m}\mathrm{i}\mathrm{n}}\eta {‖{\boldsymbol{M}}_{x}-\boldsymbol{M}{\boldsymbol{D}}_{x}-{\boldsymbol{C}}_{x}‖}_{2}^{2} $

(31)

$ {K}_{5}\left({\boldsymbol{C}}_{y}\right)=\underset{{\boldsymbol{C}}_{y}}{\mathrm{m}\mathrm{i}\mathrm{n}}\eta {‖{\boldsymbol{M}}_{y}-{\boldsymbol{D}}_{y}\boldsymbol{M}-{\boldsymbol{C}}_{y}‖}_{2}^{2} $

(32)

针对式(28)关于$ \boldsymbol{M} $的子目标函数，可以确定该目标函数为一个非凸优化问题，求$ \boldsymbol{M} $偏导并令其为0，可得

$ \begin{array}{l}\left[\right(\boldsymbol{W}\boldsymbol{B}{\boldsymbol{D}}_{y}{)}^{\mathrm{T}}\left(\boldsymbol{W}\boldsymbol{B}{\boldsymbol{D}}_{y}\right)+\kappa \boldsymbol{E}+\eta {\boldsymbol{D}}_{y}^{\mathrm{T}}{\boldsymbol{D}}_{y}]\boldsymbol{M}+\boldsymbol{M}·\eta {\boldsymbol{D}}_{x}^{\mathrm{T}}{\boldsymbol{D}}_{x}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; =(\boldsymbol{W}\boldsymbol{B}{\boldsymbol{D}}_{y}{)}^{\mathrm{T}}\mathrm{\Delta }\boldsymbol{S}+\kappa {\boldsymbol{M}}_{0}+\eta {\boldsymbol{D}}_{y}^{\mathrm{T}}{\boldsymbol{M}}_{y}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \eta {\boldsymbol{M}}_{x}{\boldsymbol{D}}_{x}^{\mathrm{T}}-\eta {\boldsymbol{D}}_{y}^{\mathrm{T}}{\boldsymbol{C}}_{y}-\eta {\boldsymbol{C}}_{x}{\boldsymbol{D}}_{x}^{\mathrm{T}}\left(33\right)\end{array} $

式中“·”表示标量乘法运算。

由于式(33)符合二维西尔维特方程的基本形式，因此可通过特征值分解技术显式求解。

针对式(29)和式(30)同时含有L₁范数和L₂范数的子目标函数，引入软阈值收缩算法^[23]求解，可依次求解得

$ {\boldsymbol{M}}_{x}^{(r+1)}=\mathrm{m}\mathrm{a}\mathrm{x}(\left|{\boldsymbol{M}}^{(r+1)}{\boldsymbol{D}}_{x}+{\boldsymbol{C}}_{x}^{\left(r\right)}\right|-\frac{\alpha }{\eta },0)\mathrm{ }·\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left[{\boldsymbol{M}}^{(r+1)}{\boldsymbol{D}}_{x}+{\boldsymbol{C}}_{x}^{\left(r\right)}\right] $

(34)

$ {\boldsymbol{M}}_{y}^{(r+1)}=\mathrm{m}\mathrm{a}\mathrm{x}\left[\left|{\boldsymbol{D}}_{y}{\boldsymbol{M}}^{(r+1)}+{\boldsymbol{C}}_{y}^{\left(r\right)}\right|-\frac{\alpha }{\eta },0\right]·\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left[{\boldsymbol{D}}_{y}{\boldsymbol{M}}^{(r+1)}+{\boldsymbol{C}}_{y}^{\left(r\right)}\right] $

(35)

式中: $ \mathrm{上}\mathrm{标}\left(r\right)\mathrm{表}\mathrm{示}\mathrm{迭}\mathrm{代}\mathrm{次}\mathrm{数};\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(Z\right) $为符号函数，即

$ \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(Z\right)=\left\{\begin{array}{cc}1& Z > 0\\ 0& Z=0\\ -1& Z < 0\end{array}\right. $

(36)

针对式(31)和式(32)为仅含有2范数的凸优化问题，采用凸优化算法进行求解，得

$ {\boldsymbol{C}}_{x}^{(r+1)}={\boldsymbol{C}}_{x}^{\left(r\right)}+\left[{\boldsymbol{M}}^{(r+1)}{\boldsymbol{D}}_{x}-{\boldsymbol{M}}_{x}^{(r+1)}\right] $

(37)

$ {\boldsymbol{C}}_{y}^{(r+1)}={\boldsymbol{C}}_{y}^{\left(r\right)}+\left[{\boldsymbol{D}}_{y}{\boldsymbol{M}}^{(r+1)}-{\boldsymbol{M}}_{y}^{(r+1)}\right] $

(38)

本文反演方法流程及框架如表 1所示。

2 模型测试

测试模型为Marmousi Ⅱ，选取其中纵向300个采样点、横向650道的数据用于建立可采性因子F和裂缝密度e的理论模型(图 1)，采样时间间隔为2ms。该模型地层产状相对平缓，在剖面中部、底部位置分别发育透镜体状(图 1中红色箭头)和薄层状油气储层(图 1中黑色虚线圈)。油气储层均表现出F和e高值异常，表明该理论模型建立的合理性，同时存在的高值突变异常有利于验证本文方法具有能够改善地层横向连续性、提高横向分辨率的优点。

2.1 无噪测试

首先，将理论模型进行高斯低通滤波处理，构建用于补充背景信息待反演参数F和e的低频模型(图 2)。其次，基于式(8)、式(11)计算出四组包含各向异性扰动部分的纵波反射系数($ {\theta }_{1}=10°\mathrm{、}{\varphi }_{2}=20°;{\varphi }_{1}=\mathrm{ }0°\mathrm{、}{\varphi }_{2}=90° $)，并分别与频率为30 Hz的雷克子波褶积，以获得无噪情况下不同入射角、不同方位的合成地震记录(图 3)；两组正交方位$ {\varphi }_{1} $、$ {\varphi }_{2} $分别为0°、90°，将无噪情况下相同入射角且不同方位的地震记录取差，获得两组无噪情况下的方位差地震记录(图 4)。最后，针对待反演参数F和e，分别采用单道反演方法(简称“单道”)、传统多道方法(简称“传统多道”)和本文方法展开无噪测试，得到无噪情况下F和e反演结果(图 5)。

由图 5可知，在无噪情况下分别采用三种方法(单道、传统多道、本文方法)反演获得的F和e无明显差异，无法直观评价反演效果。为了更加准确地反映三种方法之间的差异，体现本文方法的优势，选取剖面的第100道数据，对比三种方法获得的无噪情况下F和e单道反演结果(图 6)。由图可见，在同一地层中(图 6中黑色虚线圈)，传统多道的反演精度比单道更高，但本文方法的反演结果比传统多道更加准确。同时，在地层边界处(图 6中红色箭头处)存在F和e的突变，本文方法能更准确刻画地层边界，表明本文方法针对参数突变(异常)具有更高的稳定性。

为了更好地体现这种差异，分别利用三种方法反演获得的F和e与对应的理论模型做差，得到残差剖面(图 7)。由图可见，单道(图 7a)、传统多道(图 7b)和本文方法(图 7c)得到的F残差剖面整体上逐渐“干净”，即残差值依次减小，且位于透镜体状和薄层状的油气储层(红色箭头处)的F残差也依次减小，这说明本文方法的反演结果最准确，传统多道方法次之，单道方法效果最差。同时，“竖条状”^[31]反射也依次减少，说明相对于单道和传统多道，本文方法反演的F结果(图 5c)地层的横向连续性最好，传统多道方法(图 5b)次之，单道方法(图 5a)最差。同理，三种方法反演e结果的残差剖面(图 7d、7e、7f)也表明了采用本文方法结果更准确，地层横向连续性更好。

2.2 加噪测试

为测试本文方法的抗噪性，对图 3中的合成地震记录添加随机高斯白噪，分别得到四组信噪比(SNR)为5的不同入射角、不同方位的地震记录(图 8)，并获得两组相应的方位差地震记录(图 9)；再分别采用单道、传统多道和本文方法反演，最后得到SNR=5地震资料的F和e反演结果(图 10)。

由图 10可见，三种方法对透镜体状和薄层状的油气储层(红色箭头处)形态逐渐刻画得更好，且储层四周的伪层厚度逐渐变薄，说明本文方法受到油气储层F高值异常影响最小，反演效果最好，而传统多道方法次之，单道方法效果最差。同时，剖面的“竖条状”反射(黑色虚线圈内)得到削弱，表明横向连续性更好，由于本文方法通过反距离加权的原理对相邻道之间的关系考虑更充分，所以采用本文方法的反演获得的F剖面的“竖条状”反射明显减少，这表明进一步优化了地层横向连续性，提高了剖面的横向分辨率。

同理，三组SNR=5地震资料的e反演结果(图 10d~图 10f)也表现出和F相同的效果。因此，加噪测试结果表明：本文方法能表现出最好的准确性和横向连续性，对裂缝密度的预测更加准确、可靠。

同样，SNR=5地震资料的三种方法F和e的第100道单道反演结果如图 11所示。由图可见，采用本文方法反演获得的F和e更准确(黑色虚线圈)，传统多道方法次之，单道方法最差，这反映了传统多道方法的抗噪性比单道更好，而本文方法在传统多道方法的基础上又更准确，证明了本文方法可有效提高了F和e反演结果的准确性，同时表现出更强的抗噪性。

另外，从图 12可见，三种方法均受到噪声不同程度的影响，但本文方法获得的F残差剖面中的“竖条状”反射得到明显削弱，说明了本文方法获得的剖面横向连续性最好(红色虚线内)。尤其是在透镜体状和薄层状储层的边界及内部(红色箭头处)，残差均明显得到削弱，进一步说明了本文方法预测的F和e结果最准确。因此，通过加噪测试，相比于单道和传统多道方法，本文方法表现出更高的准确性和抗噪性，地层横向连续性更好，对裂缝密度的预测更加准确、可靠，有利于油气储层的有效预测。

2 3定量评价

引入信噪比和均方根误差(RMSE)^[32]对三种方法的反演结果进行定量评价。信噪比和均方根误差表达式分别为

$ \mathrm{S}\mathrm{N}\mathrm{R}=10\times \mathrm{l}\mathrm{g}\frac{\sum\limits_{i=1}^{{N}_{\mathrm{t}\mathrm{r}}}\sum\limits_{j=1}^{{N}_{\mathrm{s}}}({H}_{i,j}{-\overline{H})}^{2}}{\sum\limits_{i=1}^{{N}_{\mathrm{t}\mathrm{r}}}\sum\limits_{j=1}^{{N}_{\mathrm{s}}}({H}_{i,j}-{X}_{i,j}{)}^{2}} $

(39)

$ \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{\sum\limits_{i=1}^{{t}_{}}\sum\limits_{j=1}^{{n}_{}}({H}_{i,j}-{X}_{i,j}{)}^{2}}{tn}} $

(40)

式中：$ H $为理论参数值；$ \overline{H} $为理论参数的均值；X为反演得到的参数。

表 2为无噪情况下三种方法F和e反演结果的定量评价结果，表 3为SNR=5的地震资料三种方法F和e反演结果的定量评价结果。由表 2、表 3可知，本文方法反演获得的剖面信噪比均比单道和传统多道方法更高，且均方根误差最小。

为了更直观评价以上三种方法的反演精度，对不同方法的残差进行统计，获得不加噪和SNR=5时的残差统计直方图，分别如图 13和图 14所示。由图 13、图 14可见，在相同信噪比条件下，三种反演方法的反演精度由高到低依次为本文方法、传统多道、单道的反演方法；随着信噪比降低，三种方法的残差总量均增加，反演精度相比下降。

因此，从定量角度看，本文方法比单道和传统多道方法的反演精度更高，抗噪性更好。

3 实际资料应用

选取塔里木盆地YM区块实际资料（图 15）测试本文方法的应用效果。该区块碳酸盐岩埋藏深，非均质性强，裂缝较为发育。储层段发育大量高角度裂缝，具有明显的HTI介质特征。研究区OVT道集炮检距范围为565~8456 m，主要分布在3000~8000 m，方位角分布均匀，非常有利于本文方法进行各向异性反演。

采用方位差数据反演e和F仅需两个方位即可。综合各方位覆盖次数、最大炮检距以及各角度资料质量，将OVT叠前道集等分为2个方位，并在这2个方位下进行了入射角部分叠加，共得到不同方位、不同入射角的4个地震数据。数据的中心方位角分别为90°(75°~105°)、135°(120°~150°)，中心入射角分别为10°(6°~14°)、19°(15°~23°)，叠加数据如图 15所示。由图可见，振幅在不同方位和不同入射角下均存在较大的差异，适合用于各向异性反演。

分别采用本文方法与常规单道反演方法预测裂缝密度，结果如图 16所示。由图可见，两种反演方法均表明在过水平井轨迹的位置裂缝密度较大；本文方法与常规单道反演方法相比，剖面在横向上存在的“竖条状”反射(红色箭头所指)明显削弱，同时本文方法反演结果中横向上地层具有更好的连续性(图 16中红色圆圈内)，反演结果的信噪比相对有所提升。

为了进一步验证反演结果的准确性，将反演结果进一步换算成裂缝流体因子f（图 16c)。由图可见，水平井段f为高值，表明整体含气；大值与气测曲线大值(黑色虚线圆圈)吻合较好，但水平井段气测具有差异，后段含气测显示更高，与裂缝发育密度相关；水平井段后段裂缝发育密度更高，与气测全烃结果对应较好。因此，可采性因子的预测结果与水平井段气测全烃显示具有很好的一致性。这证实了本文方法具有较高的精度，通过直接反演获得的裂缝密度更加准确。

4 结论

本文提出的基于空间约束的碳酸盐岩裂缝密度直接反演的方法实现了裂缝密度的直接反演与反演算法的优化，通过模型测试和实际资料应用，得出如下结论：

(1) 引入反距离加权插值优化的横向差分算子充分利用了相邻多个地震道之间的相关性，提高了反演结果中地层的横向连续性和稳定性；

(2) 通过推导裂缝密度与包含各向异性纵波反射系数之间的关系实现了裂缝密度的直接反演，提高了裂缝预测的准确性和可靠性；

(3) 本文方法可以间接求取裂缝流体因子，能为油气开采提供指导。