0 引言

随着油气资源的不断开发，常规油气藏储量日益减少，能源的开发逐渐由浅层转向深层、由常规油气藏转向特殊油气藏，其中的裂缝型油气藏逐渐成为重点研究对象^[1]。裂缝型储层是指以裂缝为主要储集空间和渗流通道的储层。裂缝可能对储层中分散、孤立的孔隙起连通作用，增加了有效孔隙度。因此对裂缝的精确预测成为裂缝型油气藏勘探和开发的重点和难点^[1]。裂缝分布复杂、规律性差，裂缝型储层呈较强的非均质性，利用常规地震反演方法预测裂缝具有一定的局限性^[2]。

如今，利用地震技术预测裂缝的方法很多，其中横波勘探方法最有效，因为横波分裂的原因是各向异性，但采集成本太高难以在生产中实现^[3]。基于横波分裂的多分量转换波裂缝检测方法虽然降低了成本，但对地震资料的品质要求很高^[4]，且对转换波地震资料的地质解释还处于探索阶段，因此也难以投入实际生产。叠后地震属性分析较容易实现，成本相对较低，但裂缝储层的复杂性导致单一地震属性预测结果的多解性很强。目前应用最为广泛的是纵波各向异性裂缝预测方法，考察由裂缝引起的纵波动力学属性随着方位角的变化特征——纵波的方位各向异性。

自20世纪80年代以来，随着对各向异性研究的深入，应用纵波各向异性预测裂缝技术发展迅速。Crampin^[5]证实裂缝介质呈现各向异性。Thomsen^[6]给出了VTI介质对称平面内水平反射界面的纵、横波NMO速度表达式。Tsvankin^[7]将Thomsen各向异性参数推广到HTI介质中，并给出了任意各向异性强度下水平界面纯纵波叠加速度的精确表达式。Rüger^[8]基于弱各向异性的概念，推导了HTI介质纵波反射系数与裂缝参数之间的解析关系，结果表明，纵波AVO梯度在平行于裂缝走向和垂直于裂缝走向的两个主方向上存在较大差异，这是进行纵波AVAZ裂缝检测的理论基础。曲寿利等^[9]提出了波阻抗随方位变化(IPVA)裂缝检测新方法，得到了稳定的反演结果。Zhu等^[10-11]将Thomsen的思想用于衰减各向异性反演及体波衰减系数反演。

根据动力学属性可知，纵波方位各向异性裂缝检测方法分为基于速度或旅行时方位各向异性、基于振幅方位各向异性、基于弹性波阻抗方位各向异性和基于衰减方位各向异性的裂缝检测方法^[12]。基于速度方位各向异性的裂缝检测方法较稳定，但只能识别大套储层，识别薄储层时分辨率不高^[13]。基于振幅方位各向异性的裂缝检测方法对介质的各向异性程度敏感^[14]，分辨率高、易操作，但抗干扰能力差，对地震数据有较高的要求^[15]。相对而言，基于弹性波阻抗各向异性的裂缝检测方法具有更高的稳定性和分辨率，但对前期资料要求高，受子波提取的影响较大。

基于衰减方位各向异性的裂缝检测方法抗噪性好，可以较精确地反映地下波场变化，其关键点也是难点在于计算衰减系数，目前很难求取精确的Q值，因此基于衰减方位各向异性的裂缝检测方法暂时无法得到广泛应用。

孔丽云等^[16]对三层裂缝孔隙模型的正演模拟发现，当入射角大于30°时，目的层顶界面振幅和底界面旅行时呈较强的各向异性特征。因此，本文采用顶界面振幅和底界面旅行时(速度)联合反演的方法。对于每一个微小时窗，将该时窗底界速度反演得到的各向异性梯度作为该时窗顶界振幅反演的约束，移动时窗进而得到一个反演体。该方法克服了速度反演分辨率低和振幅反演稳定性差的问题，可得到更好的反演效果。

1 方法原理 1.1 基于振幅各向异性裂缝反演

Rüger^[8]提出了HTI介质中振幅随方位角变化的纵波反射系数公式

$ \begin{array}{l} {R_{\rm{P}}}\left( {i, \varphi } \right) = \frac{1}{2}\frac{{\Delta z}}{{\bar z}} + \frac{1}{2}\left\{ {\frac{{\Delta \alpha }}{{\bar \alpha }} - {{\left( {\frac{{2\bar \beta }}{{\bar \alpha }}} \right)}^2}\frac{{\Delta G}}{{\bar G}} + \left[ {\Delta \delta + {{\left( {\frac{{2\bar \beta }}{{\bar \alpha }}} \right)}^2}\Delta \gamma } \right]{{\cos }^2}\varphi } \right\}{\sin ^2}i + \\ \frac{1}{2}\left( {\frac{{\Delta \alpha }}{{\bar \alpha }} + \Delta \varepsilon {{\cos }^4}\varphi + \Delta \delta {{\sin }^2}\varphi {{\cos }^2}\varphi } \right){\sin ^2}i{\tan ^2}i \end{array} $

(1)

式中；i为入射角；φ为观测方位与裂缝对称轴的夹角；z为纵波阻抗；G为剪切模量；α为纵波速度；β为横波速度；Δ·表示上、下介质岩性参数之差；·表示上、下介质岩性参数的平均值；δ、γ、ε为HTI介质Thomsen参数。在弱各向异性的假设条件下，当i较小时，舍去有关i的三角函数高次项，式(1)简化为

$ {R_{\rm{P}}}\left( {i, \varphi } \right) = A + {B^{{\rm{iso}}}}{\sin ^2}i + {B^{{\rm{ani}}}}{\sin ^2}i{\cos ^2}\left( {\phi - {\phi ^{{\rm{sym}}}}} \right) $

(2)

式中：φ=ϕ-ϕ^sym, ϕ为观测方位, ϕ^sym为裂缝对称轴的方位角；A为常数；B^iso为各向同性梯度；B^ani为各向异性梯度。

Mallick等^[17]对式(2)进一步简化，得

$ R\left( \varphi \right) = a + b\cos 2\varphi $

(3)

式中：a为各向同性项；bcos2φ为各向异性项。在实际应用中，一般采用椭圆拟合法预测裂缝密度和方位，通常以椭圆长轴(a+b)的方向指示裂缝走向，以椭圆扁率(a-b)/(a+b)指示裂缝密度。在实际操作中一般将入射角取为固定值或采用分方位叠后道集进行反演，并没有同时利用多个入射角和多个方位角信息^[18]。

1.2 基于振幅和速度各向异性联合反演

首先，对式(2)进行线性化处理，得

$ \begin{array}{l} {R_{\rm{P}}}\left( {i, \phi } \right) = A + {A_1}{\sin ^2}i{\cos ^2}\phi + {A_2}{\sin ^2}i{\sin ^2}\phi + \\ {A_3}{\sin ^2}i\sin \phi \cos \phi \end{array} $

(4)

以适应多方位角、多入射角情况，实现了从叠后分方位资料到叠前分方位资料的应用。其中

$ \left\{ {\begin{array}{*{20}{l}} {{A_1} = {B^{{\rm{iso}}}} + {B^{{\rm{ani}}}}{{\cos }^2}{\phi ^{{\rm{sym}}}}}\\ {{A_2} = {B^{{\rm{iso}}}} + {B^{{\rm{ani}}}}{{\sin }^2}{\phi ^{{\rm{sym}}}}}\\ {{A_3} = {B^{{\rm{ani}}}}\sin 2{\phi ^{{\rm{sym}}}}} \end{array}} \right. $

(5)

因此

$ \left\{ {\begin{array}{*{20}{l}} {{\phi ^{{\rm{sym}}}} = \frac{1}{2}{\rm{a}}\tan 2\left( {{A_3}, {A_1} - {A_2}} \right)}\\ {{B^{{\rm{ani}}}} = \frac{{{A_3}}}{{\sin \left( {2{\phi ^{{\rm{sym}}}}} \right)}}}\\ {{B^{{\rm{iso}}}} = \frac{1}{2}\left[ {\left( {{A_1} + {A_2}} \right) - {B^{{\rm{ani}}}}} \right]} \end{array}} \right. $

(6)

式(6)中atan2为四象限反正切函数，其函数值分布在(-π, π)，所在象限由A₃和(A₁-A₂)的正负决定。裂缝密度与各向异性梯度B^ani之间有良好的对应关系^[19]，因而可以用B^ani表示裂缝密度。

由式(4)可形成不同入射角、不同方位角组成的方程组，方程组矩阵形式为

$ \left( {\begin{array}{*{20}{c}} {R\left( {{i_1}, {\phi _1}} \right)}\\ {R\left( {{i_1}, {\phi _2}} \right)}\\ \vdots \\ {R\left( {{i_M}, {\phi _{N - 1}}} \right)}\\ {R\left( {{i_M}, {\phi _N}} \right)} \end{array}} \right){\rm{ = }}\left( {\begin{array}{*{20}{c}} 1&{{{\sin }^2}{i_1}{{\cos }^2}{\phi _1}}&{{{\sin }^2}{i_1}{{\sin }^2}{\phi _1}}&{{{\sin }^2}{i_1}\sin {\phi _1}\cos {\phi _1}}\\ 1&{{{\sin }^2}{i_1}{{\cos }^2}{\phi _2}}&{{{\sin }^2}{i_1}{{\sin }^2}{\phi _2}}&{{{\sin }^2}{i_1}\sin {\phi _2}\cos {\phi _2}}\\ \vdots & \vdots & \vdots & \vdots \\ 1&{{{\sin }^2}{i_M}{{\cos }^2}{\phi _{N - 1}}}&{{{\sin }^2}{i_M}{{\sin }^2}{\phi _{N - 1}}}&{{{\sin }^2}{i_M}\sin {\phi _{N - 1}}\cos {\phi _{N - 1}}}\\ 1&{{{\sin }^2}{i_M}{{\cos }^2}{\phi _N}}&{{{\sin }^2}{i_M}{{\sin }^2}{\phi _N}}&{{{\sin }^2}{i_M}\sin {\phi _N}\cos {\phi _N}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} A\\ {{A_1}}\\ {{A_2}}\\ {{A_3}} \end{array}} \right) $

(7)

式中R(i_m, ϕ_n)为第m(m=1, …, M)个入射角、第n(n=1, …, N)个方位角的地震道对应的振幅, M为入射角总个数, N为方位角总个数。式(4)可用于叠前方位道集裂缝定量预测，同时利用了方位和入射角信息，精度更高。式(4)~式(7)为联合反演的基础，以速度反演得到的各向异性梯度作为约束。Tsvankin^[7]提出了HTI介质相速度公式

$ {V_{{\rm{P0}}}}\left( {i, \phi } \right) = {V_{{\rm{iso}}}}\left[ {1 + \delta {{\sin }^2}i{{\cos }^2}\left( {\phi - {\phi ^{{\rm{sym}}}}} \right) + \left( {\varepsilon - \delta } \right){{\sin }^4}i{{\cos }^4}\left( {\phi - {\phi ^{{\rm{sym}}}}} \right)} \right] $

(8)

式中：V_P0为纵波相速度，由拉平不同方位角、不同入射角道集得到的各向异性时差求得；V_iso为纵波各向同性速度。根据式(8)，利用共轭梯度法可反演δ、ε。由式(2)可知，B^ani与δ、ε之间的关系为

$ {B^{{\rm{ani}}}} = \frac{1}{2}\left[ {\Delta \delta + 2{{\left( {\frac{{2\bar \beta }}{{\bar \alpha }}} \right)}^2}\Delta \gamma } \right] $

(9)

无法通过反演得到各向异性参数γ，在一般情况下ε和γ的单调性一致，同时增、减或为零，因此可假设二者间具有线性关系^[20]。根据测井数据拟合ε与γ的线性关系即可得到γ。已知δ和γ，根据式(9)即可由速度各向异性反演得到各向异性梯度B_v^ani，再构建振幅反演的初始目标函数

$ \min {\left( {\mathit{\boldsymbol{S}} - \mathit{\boldsymbol{WUC}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{S}} - \mathit{\boldsymbol{WUC}}} \right) $

(10)

将B_v^ani作为约束加入之后得到新目标函数

$ \min \left[ {{{\left( {\mathit{\boldsymbol{S}} - \mathit{\boldsymbol{WUC}}} \right)}^{\rm{T}}}\left( {\mathit{\boldsymbol{S}} - \mathit{\boldsymbol{WUC}}} \right) + \lambda {{\left( {{\mathit{\boldsymbol{B}}^{{\rm{ani}}}} - \mathit{\boldsymbol{B}}_v^{{\rm{ani}}}} \right)}^{\rm{T}}}\left( {{\mathit{\boldsymbol{B}}^{{\rm{ani}}}} - \mathit{\boldsymbol{B}}_v^{{\rm{ani}}}} \right)} \right] $

(11)

其中

$ \mathit{\boldsymbol{U}}{\rm{ = }}\left( {\begin{array}{*{20}{c}} 1&{{{\sin }^2}{i_1}{{\cos }^2}{\phi _1}}&{{{\sin }^2}{i_1}{{\sin }^2}{\phi _1}}&{\sin {\phi _1}\cos {\phi _1}{{\sin }^2}{i_1}}\\ 1&{{{\sin }^2}{i_1}{{\cos }^2}{\phi _2}}&{{{\sin }^2}{i_1}{{\sin }^2}{\phi _2}}&{\sin {\phi _2}\cos {\phi _2}{{\sin }^2}{i_1}}\\ \vdots & \vdots & \vdots & \vdots \\ 1&{{{\sin }^2}{i_M}{{\cos }^2}{\phi _{N - 1}}}&{{{\sin }^2}{i_M}{{\sin }^2}{\phi _{N - 1}}}&{\sin {\phi _{N - 1}}\cos {\phi _{N - 1}}{{\sin }^2}{i_M}}\\ 1&{{{\sin }^2}{i_M}{{\cos }^2}{\phi _N}}&{{{\sin }^2}{i_M}{{\sin }^2}{\phi _N}}&{\sin {\phi _N}\cos {\phi _N}{{\sin }^2}{i_M}} \end{array}} \right) $

(12)

式中：B^ani、B_v^ani分别为B^ani、B_v^ani的矩阵形式；C=(A A₁ A₂ A₃)^T；S为振幅矩阵；W为子波矩阵。求解式(11)可以得到C，进而根据式(6)求得B^ani和ϕ^sym。图 1为基于振幅和速度各向异性联合反演流程。

2 模型测试

对测试模型(图 2)正演得到叠前全方位角道集(图 3)，利用叠前全方位角道集进行振幅和速度各向异性联合反演，得到裂缝密度(图 4)和裂缝方位(图 5)。其中：模型顶界面的裂缝密度反演值为0.1372(真实值为0.1387)，相对误差为1.08%；裂缝走向方位角反演值约为0°(真实值为0°)。可见反演结果较准确。在叠前全方位角道集中加入信噪比为10dB的高斯白噪声后进行振幅和速度各向异性联合反演，得到的顶界面裂缝密度反演值为0.1411(图 6)，相对误差为1.73%；裂缝走向方位角反演结果相对杂乱(图 7)，但仍以0°为主，说明方法切实可行。

3 实际资料应用 3.1 工区概况

工区位于四川盆地中部遂宁市安居区与合川市潼南县的交界处。构造位置位于四川盆地乐山—龙女寺古隆起东段斜坡部位。乐山—龙女寺古隆起是一个长期继承性发育的巨型隆起，桐湾及加里东运动使二叠系以下地层遭受不同程度剥蚀，特别是对震旦系灯影组(Z₂dn)及寒武系—奥陶系( -O)的风化剥蚀，形成了灯4段(Z₂dn₄)、下寒武统龙王庙组(₁l)储层。钻探成果表明：乐山—龙女寺古隆起的震旦系—寒武系(Z-)储层发育，Z顶发育4个大型背斜圈闭，其中高石梯构造面积最大，达到365km²，具有较好的含油气性。

高石梯构造白云岩天然裂缝具有多种成因类型、多种尺度。根据裂缝的规模和当前的测量和预测手段，一般可以预测两类裂缝。一是大尺度裂缝，延伸距离长，可达几十米到几千米，断距大，地震剖面可以直接解释及识别；二是中、小尺度裂缝，延伸距离小，可达几厘米到几十米，断距小，岩心和成像测井显示为小尺度裂缝。工区内共收集到13口井的常规测井资料、11口井的成像测井资料，其中8口有横波测井数据。经过资料筛选，在Z₂dn₄中共处理、解释10口井数据。根据测井、地质综合评价将Z₂dn₄储层划分为裂缝—孔隙—孔洞复合型、孔隙—孔洞型、孔隙型三类，其中有效储层的主要储集空间类型为前二类。

3.2 实际资料应用效果

工区观测系统纵横比为0.8325，选取包含3个入射角(5°、15°、30°)、12个方位角的叠前宽方位地震资料。

图 8为不同方法得到裂缝密度剖面。由图可见：①叠前振幅各向异性反演结果(图 8b)较叠后振幅各向异性反演结果(图 8a)的横向连续性更好，反映了更多的层状信息，精度大幅提高(图 8a、图 8b白框区域)，但前者(图 8b)的井曲线与剖面吻合度不高。②叠前速度各向异性反演结果(图 8d)与井裂缝曲线吻合较好，可作为振幅反演的约束，但较叠前振幅各向异性反演结果(图 8b)的分辨率低。③叠前振幅和速度各向异性联合反演结果(图 8c)较叠前振幅各向异性反演结果(图 8b)明显提高了与井曲线的吻合度(图 8b、图 8c蓝框区域)，说明前者有效提高了反演精度，且较叠前速度各向异性反演结果(图 8d)的分辨率高。④井裂缝密度曲线表明，在₁q~Z₂dn₃的上半层(₁q~2240ms)，gs1井、gs9井、gs11井裂缝密度依次降低；在₁q~Z₂dn₃的下半层(2240ms~Z₂dn₃)，gs9井、gs1井、gs11井裂缝密度依次降低。因此gs11井处裂缝不发育，各向异性不显著，导致反演效果不理想。

图 9为叠前振幅各向异性反演、叠前振幅与速度各向异性联合反演得到的₁q~Z₂dn₃上半层裂缝密度切片。由图可见：叠前振幅各向异性反演结果(图 9a)和图 8b的井裂缝密度曲线分布规律不一致；叠前振幅与速度各向异性联合反演结果(图 9b)和图 8c的井裂缝密度曲线分布规律较一致。图 10为叠前振幅各向异性反演、叠前振幅与速度各向异性联合反演得到的₁q~Z₂dn₃下半层裂缝密度切片。由图可见：叠前振幅各向异性反演(图 10a)、叠前振幅与速度各向异性联合反演(图 10b)在井点位置的反演结果分别与图 8b、图 8c的井裂缝密度曲线分布规律一致，但就整体分布而言，图 10a、图 10b的裂缝分别集中于gs1井、gs9井附近；裂缝密度曲线(图 8)表明，gs9井区的裂缝较gs1井区更发育，因此图 10b的反演效果更精确。

图 11为叠前振幅各向异性反演、叠前振幅与速度各向异性联合反演得到的₁q~Z₂dn₃裂缝走向方位角切片，图 12为gs1井裂缝走向玫瑰图。由图可见，叠前振幅各向异性反演结果(图 11a)、叠前振幅与速度各向异性联合反演结果(图 11b)相近，且在gs1井处裂缝走向方位角接近90°(红色)，与gs1井裂缝走向玫瑰图(图 12)较吻合。

综上所述，与叠前、叠后振幅各向异性反演相比，叠前振幅与速度各向异性联合反演的精度最高。

4 结论

本文分析了现有基于方位各向异性的裂缝检测方法的不足，提出了利用叠前振幅与速度各向异性的联合反演方法，模型测试和实际资料应用效果表明：

(1) 与叠前、叠后振幅各向异性反演相比，文中方法的反演精度最高。

(2) 在gs1井区，叠前速度各向异性反演结果与井裂缝曲线吻合较好，可作为振幅反演的约束，但较振幅各向异性反演的分辨率低，二者联合反演改善了速度反演分辨率低的问题。