基于各向异性梯度的裂缝密度反演

秦海旭 吴国忱

（中国山东青岛266580中国石油大学（华东）地球科学与技术学院）

引言

随着油气勘探工作的发展，勘探的难度不断加大，油气勘探进入复杂构造和岩性油气藏勘探阶段，勘探开发的重点落在越来越复杂的非均质储层上.近年来碳酸盐储层研究受到越来越多的关注（熊翥，2005），其中裂缝是许多油气藏中流体或者气体的重要通道，因此了解裂缝储集层的地震响应，获取裂缝发育方向，识别裂缝流体的性质对于油气藏勘探和开发具有重要意义.利用纵波资料检测裂缝是目前一个重要课题（Rüger，1997；Chichinina et al，2006；Downton，Gray，2006）.研究裂缝的常规方法是将裂缝看成散射体，用散射理论研究裂缝介质的地震响应（刁顺等，1994a，b），但此方法的效率非常低.近年来将裂缝等效为高陡倾角裂缝介质，运用等效理论（Crampin，1985；Sayers，1988；Bakulin et al，2000a，b，c）将裂缝介质等效为各向异性介质，根据各向异性介质理论研究裂缝介质的响应有了长足的发展.实际地层中裂缝的形成受多种因素控制，其物理属性复杂，表现出强烈的各向异性.由于上覆地层的压力使得水平或低陡倾角裂缝存在较少，大多数为高陡倾角裂缝，因此本文采用裂缝介质等效理论将高陡倾角裂缝等效为横向各向同性（HTI）介质便于实际应用.

裂缝密度等裂缝特征参数是描述裂缝性储集层的重要参数，从地震资料中得到裂缝特征参数对油气田的勘探开发具有重要意义.通过横波分裂等手段已经能够提取裂缝介质的裂缝特征参数（李文林，李承楚，1994；何樵登，陶春辉，1995），但横波或多波采集处理花费较大，不便于实际推广.研究表明，纵波对各向异性介质的振幅响应也很强烈（Mallick，Frazer，1991；Lynn et al，1995，1996a，b），且利用纵波资料也能识别裂缝特征响应（Thomsen，1986；Tsvankin，Thomsen，1995；Lynn et al，1996a，b）.国内很多研究人员对裂缝密度的地震响应进行了研究（朱培民等，2001；魏建新，2002；李琼等，2006；田锋，2007），均证明基于纵波资料反演裂缝密度具有一定的可行性.

本文提出一种计算裂缝密度的新方法.运用等效理论将裂缝介质等效为各向异性介质，对于高陡倾角裂缝，将裂缝介质等效为HTI介质.研究HTI介质的AVO（amplitude versus offset）响应与裂缝参数的关系，根据Rüger（1997）推导的各向异性反射系数方程，得到裂缝介质的各向异性梯度与裂缝密度的关系，最终求出裂缝介质的裂缝密度.

1 裂缝密度等效各向异性参数

利用裂缝介质各向异性梯度求取裂缝介质的裂缝密度，需要利用等效理论将裂缝介质等效为各向异性介质.目前广泛存在的裂缝等效理论主要有Thomsen等效理论（Thomsen，1986）、Hudson等效理论（Hudson，Liu，1999）及线性滑动理论（Bakulin et al，2000a，b，c）等.Thomsen等效理论相当于低频模型，假设流体流动无限慢，对于一些特定裂缝的等效结果是比较准确的，但是该方法限制条件较多不利于实际研究.Hudson等效理论与Thomsen等效理论正好相反，相当于高频模型，其假设流体流动无限快，也不利于实际研究.线性滑动理论是将裂缝看作地层中的一个柔性面，该柔性面可以用裂缝的切向柔度和法向柔度线性表示，且裂缝的切向柔度与法向柔度符合线性滑动边界条件便于求解；对于高陡倾角裂缝运用等效理论将裂缝介质等效为HTI介质.含有裂缝介质的岩石总柔度等于裂缝柔度与围岩柔度的和，假设围岩为各向同性，即

式中，S是岩石总的柔度矩阵，Sf是裂缝的柔度矩阵，Sb是围岩的柔度矩阵.其中柔度矩阵为弹性矩阵C的倒数即S＝C－1.根据线性滑动理论有

根据S＝Sb＋Sf与S＝C－1得C＝Cb－Cf，即

式中，en为裂缝的法向柔度，et为裂缝的切向柔度，λ和μ为拉梅常数.式（4）是将裂缝介质的参数等效为本构坐标系下各向异性介质的弹性参数.

根据Rüger（1997）给出的HTI介质各向异性参数δ（v）和γ的表达式

将式（4）带入式（5）可得到

当en≤1且et≤1时，式（6）可简化为

式中，η＝μ／（λ＋2μ），ζ＝λ／（λ＋2μ），λ＝ρ（α2－2β2），μ＝ρβ2，α和β分别为纵波和横波速度.

式（7）中，求取裂缝介质的各向异性参数必须已知裂缝介质的切向柔度与法向柔度.根据Rüger（1997）的理论，裂缝充填物只影响裂缝的法向柔度，而不影响裂缝的切向柔度.而裂缝的切向柔度可以表示为

式中，Φc为总裂缝孔隙度，A0为断裂面上裂缝面积所占总面积的比例.对于非黏滞性流体有

式中，a为裂缝纵横比，η＝β2／α2.将式（9）代入式（8），得到最终裂缝切向柔度为

裂缝介质饱和度为0时，法向柔量egn的表达式为

裂缝介质饱和度为100%时，法向柔量eln的表达式为

式中：Fn为裂缝流体充填因子，Fn＝ρlα2l／（ρα2），其中ρl，ρ，αl和α分别为裂缝充填流体和背景各向同性介质的密度和纵波速度.因此，裂缝法向柔度在部分饱和裂缝储层中的表达式为

式中Sl是裂缝中的流体饱和度.式（13）成立的条件是裂缝中多相流体是不相容的.根据Reid和MacBeth（2000）可知，天然裂缝储层中多相流体是不相容的.当Sl＝0时，en＝egn；当Sl＝1时，则en＝eln.对于多组裂缝同时存在的情况，根据矢量叠加原理，将裂缝的切向柔度et与裂缝的法向柔度en进行矢量叠加，求出裂缝组合的切向柔度与法向柔度，然后将裂缝介质等效为各向异性介质.

将式（13）和式（10）带入式（7），可以求取裂缝介质的各向异性参数为

式中，δ和ε为各向异性参数，η＝β2／α2，ζ＝λ／（λ＋2μ）＝1－2η，λ＝ρ（α2－2β2），μ＝ρβ2，Sl为饱和度.

根据Mavko等（1998）结果，裂缝密度e与孔隙度Φc的关系为

根据式（15）得

将式（16）带入式（14），可以得到裂缝密度与各向异性参数γ和δ的关系式为

式（17）中各向异性参数与裂缝密度e、饱和度Sl及裂缝类型A0有关.假设该地区的饱和度和裂缝类型已知，则可以得到各向异性参数与裂缝密度的关系.式（17）可表示为

其中

式中A，B，C为与背景参数和裂缝类型有关的参数.由式（18）可知各向异性参数与裂缝密度成线性关系，如图1所示.

图1 各向异性参数（δ和γ）与裂缝密度（e）的关系示意图Fig.1 Relationship between the anisotropic parameters（δandγ）and fracture density（e）

2 裂缝密度与各向异性梯度的关系

运用等效理论将裂缝介质等效为各向异性介质.本文针对高角度裂缝，将裂缝介质等效为HTI介质.根据HTI介质的反射系数方程可以求出裂缝介质的近似反射系数随偏移距的变化.Rüger（1997）给出了HTI介质中PP波近似反射系数随方位角的变化公式为

式中，RPP（θ，φ）为反射系数，θ为入射角，φ为对称轴与x方向的夹角.当各向同性背景上发育一组裂缝时，PP波反射系数近似为

其中

根据式（19）、（20）和（21）可得

式中，G（φ）为介质总梯度，GI为各向同性梯度，GA为各向异性梯度.将式（18）带入式（22）得

由式（23）可得到各向异性梯度与裂缝密度的关系.

设计两个饱和度不同的裂缝介质模型，其中模型1的饱和度为0，模型2的饱和度为1.两个模型的纵波速度均为6 000m／s，横波速度均为3 150m／s，裂缝引起的饱和度固定；裂缝长度为1m，宽度为0.001m，A0＝0.1，Fn＝0.02；裂缝密度为1—20条／m.各向异性梯度与裂缝密度的关系如图2所示.图2a中饱和度为1，相当于裂缝中充满液体，此时裂缝的各向异性梯度与裂缝密度呈线性关系，直线的斜率为正值；图2b中饱和度为0，相当于裂缝中充满气体，此时各向异性梯度与裂缝密度也呈线性关系，直线的斜率为负值.可以看出，当裂缝完全含有气体时，各向异性梯度为负值；而当裂缝完全含液体时，各向异性梯度为正值.在饱和度未知的情况下，各向异性梯度可以作为流体的指示器，不同模型的各向异性梯度与裂缝密度变化趋势不同.Reid和MacBeth（2000）指出，该特征对于纵横波速比在1.81—2.65之间成立.大多数岩石满足这种情况.各向异性梯度的斜率也可以作为流体指示器，针对上述模型若直线斜率是负值可以认为裂缝中充填的是气体，直线斜率是正值可以认为裂缝中充填的是液体，不同模型的各向异性梯度变化规律不同.

图2 饱和度为1（a）和0（b）时各向异性梯度（GA）与裂缝密度（e）的关系Fig.2 Relationship between the anisotropic gradient（GA）and fracture density（e）when the saturation is 1（a）and 0（b）

3 裂缝密度反演方法

根据上面推导可知，裂缝密度与各向异性梯度存在一一对应的关系.实际生产中可以根据各向异性梯度得到裂缝介质的裂缝密度，具体步骤为：

1）提取每个共中心点（CMP）处不同方位的反射系数，进而求取不同方位的梯度；

2）利用G（φ）＝GI＋GAsin2φ，分别计算总梯度与各向同性梯度和各向异性梯度的关系；

4）根据式（15）—（23）计算裂缝介质的裂缝密度，其中计算裂缝密度的公式为

式（24）即为用各向异性梯度计算裂缝密度的公式.式中：GA为各向异性梯度；A，B，C为与裂缝有关的参数.在反演裂缝密度时需要已知多方位的地震数据、背景速度及裂缝类型.式（24）中用到一个含油水饱和度的信息，在裂缝密度反演之前需要已知该地区含油气情况，若目标区主要含气则饱和度为0，若目标区主要含油水则饱和度为1.

4 模型试算

4.1 单点模型试算

设计一简单模型，上层为各向同性介质，下层为裂缝介质，裂缝长度为1m，裂缝宽度为0.001m，A0＝0.1，Fn＝0.02.第一层纵波速度为5 000m／s，横波速度为2 500m／s，密度为2 600kg／m3；第二层纵波速度为4 000m／s，横波速度为2 200m／s，密度为2 100kg／m3，裂缝密度为2条／m.裂缝介质的饱和度为0，即完全充填的是气体.模型精确的反射系数如图3所示.图3a为不含裂缝的反射系数，可以看出当介质不含裂缝时不同方位的反射系数都一样；图3b为含有裂缝的反射系数，可以看出当介质含有裂缝时不同方位的反射系数差别很大，AVO类型会发生变化，这为利用不同方位的地震数据进行裂缝特征参数反演提供了依据.

图3 不同方位角（φ）情况下的模型反射系数示意图（a）不含裂缝介质；（b）含有裂缝介质Fig.3 The sketch of reflection coefficient in the case of several azimuths（a）No fractured medium；（b）Fractured medium

实际地震资料中含有噪声，为了验证本文方法的有效性，对充满气体的模型分别设计信噪比为10，4，2和1.每个信噪比考虑40次实现，相当于40个共中心点（CMP）.分别计算每次实现的不同方位的反射系数，然后褶积一个地震子波生成不同方位的地震数据.从合成的地震记录出发，计算不同方位的AVO梯度，从而获取该点的各向异性梯度，最后根据式（24）计算裂缝介质的裂缝密度.图4—7分别给出了信噪比为10，4，2和1的计算结果.其中，a图表示从地震记录提取的反射系数，可以看出，不同方位的反射系数差别很大，且信噪比越低，提取的反射系数与精确的反射系数差别越大；b图表示从反射系数中提取的裂缝介质各向异性梯度，并由此计算出的裂缝密度.可以看出，信噪比越高，提取的反射系数更集中，计算得到的裂缝密度值在平均值附近比较集中，计算的裂缝密度与实际较接近.信噪比等于1时，也能得到比较好的结果，但是对于其中的一次实现有可能远离实际值.表1给出了不同信噪比情况下计算的平均裂缝密度.可以看出，不同信噪比对于多次实现得到的平均值比较接近实际情况，偏差不超过5%.但是当信噪比等于1时，一次实现其偏差可能超过30%.信噪比小于1的资料无法使用.从图4—7中可以看出，由一次实现提取的反射系数计算得到的裂缝介质各向异性梯度，并由此计算的裂缝介质的裂缝密度不同，但是由多次实现计算的结果求取的平均裂缝密度却基本相同，误差小于5%.实际计算中可选用多条测线反演裂缝密度的结果取平均值作为实际反演的结果.

表1 由不同信噪比地震记录计算的平均裂缝密度Table 1 The average fracture density from the seismic data with different signal to noise ratios

图4 信噪比为10的反射系数（a）与计算得到的裂缝密度（b）Fig.4 Reflection coefficients（a）and fracture density（b）with signal-to-noise ratio of 10

图5 信噪比为4的反射系数（a）与计算得到的裂缝密度（b）Fig.5 Reflection coefficients（a）and fracture density（b）with signal-to-noise ratio of 4

图6 信噪比为2的反射系数（a）与计算得到的裂缝密度（b）Fig.6 Reflection coefficients（a）and fracture density（b）with signal-to-noise ratio of 2

图7 信噪比为1的反射系数（a）与计算得到的裂缝密度（b）Fig.7 Reflection coefficients（a）and fracture density（b）with signal-to-noise ratio of 1

4.2 单道模型试算

抽取Marmous模型中其中一道作为模型数据（图8）.图8中，如果流体指示因子D等于0则孔隙中充填液体，即饱和度为1；若D小于0则孔隙中充填气，即饱和度为0.模型中设置两个含气位置，如红线所在区域.利用雷克子波合成不同信噪比不同方位的地震记录，然后从该记录出发，求取不同方位的梯度，得到每个采样点处的各向异性梯度，根据该地区的含流体情况，最终反演得到的裂缝密度结果如图9所示.可以看出：不同信噪比反演裂缝密度的位置都正确，只是当信噪比较低时，反演得到的裂缝密度不准确；信噪比等于1时得到的裂缝密度反演结果与已知裂缝密度差别较大，因此实际反演中应采用信噪比大于1的地震数据.

图8 模型示意图（a）纵波速度；（b）横波速度；（c）密度；（d）流体指示因子；（e）裂缝密度Fig.8 The sketch of models（a）P wave velocity；（b）S wave velocity；（c）Density；（d）Fluid factor；（e）Fracture density

图9 不同信噪比的地震记录得到的裂缝密度反演结果Fig.9 The inversion results for fracture density from the seismic data with several SNRs

4.3 多道反演试算

抽取Marmous模型中一部分作为模型数据（图10），利用雷克子波合成不同信噪比不同方位的地震记录，然后从该记录出发，求取每个采样点不同方位的梯度，得到每个采样点处的各向异性梯度，根据该地区的含流体情况，最终反演得到的裂缝密度结果如图11所示.可以看出，不同信噪比的地震记录反演得到的裂缝密度位置都正确，当信噪比等于或小于1时，反演得到的裂缝密度值与已知的裂缝密度值差别较大.因此，实际反演中应采用信噪比大于1的地震数据.

图10 模型示意图（a）纵波速度；（b）横波速度；（c）密度；（d）裂缝密度Fig.10 The sketch of models（a）P wave velocity；（b）S wave velocity；（c）Density；（d）Fracture density

5 讨论与结论

图11 裂缝密度反演结果（a）信噪比为10；（b）信噪比为4；（c）信噪比为2；（d）信噪比为1Fig.11 The inversion results of fracture density（a）SNR＝10；（b）SNR＝4；（c）SNR＝2；（d）SNR＝1

本文提出一种利用各向异性梯度计算裂缝密度的新方法.利用裂缝等效理论将高陡倾角裂缝介质等效为HTI介质，推导出各向异性梯度与裂缝密度的关系式，并根据此关系从地震数据反演得到裂缝密度.模型试算表明，该方法能得到裂缝密度，反演结果与初始模型基本一致证明了该方法的正确性；地震记录添加一定的信噪比后也能反演得到裂缝密度，证明了该方法的稳定性，但当该地区的信噪比等于1时，反演得到的裂缝密度结果与实际结果差别较大，故实际反演裂缝密度时应采用信噪比大于1的地震数据.需要说明的是：本文结果仅仅针对一组高陡倾角裂缝介质，而实际地层中裂缝的存在形式是多种多样的，需进一步研究实际裂缝介质中的裂缝密度反演方法；本文仅仅从各向异性梯度反演裂缝密度，有必要进一步建立地震数据或反射系数与裂缝密度的直接关系，减少反演累计误差，为裂缝性储层反演与储层预测提供依据；本文仅仅针对理论模型进行裂缝密度反演，下一步有必要进行实际数据的裂缝密度反演方法的研究.

刁顺，杨慧珠，许云.1994a.裂缝及介质非均匀散射［J］.石油物探，33（4）：49-55.

Diao S，Yang H Z，Xu Y.1994a.Scattering due to inhomogeneity of media［J］.Geophysical Prospecting for Petroleum，33（4）：49-55（in Chinese）.

刁顺，杨慧珠，许云.1994b.TI介质裂纹散射研究［J］.石油地球物理勘探，29（5）：545-551.

Diao S，Yang H Z，Xu Y.1994b.Research on crack scattering in transversally isotropic medium［J］.Oil Geophysical Prospecting，29（5）：545-551（in Chinese）.

何樵登，陶春辉.1995.用遗传算法反演裂隙各向异性介质［J］.石油物探，34（3）：46-50.

He Q D，Tao C H.1995.Inversing fractured anisotropic media through a genetic algorithm［J］.Geophysical Prospecting for Petroleum，34（3）：46-50（in Chinese）.

李琼，李勇，李正文，吴朝容.2006.基于GA-BP理论的储层视裂缝密度地震非线性反演方法［J］.地球物理学进展，21（2）：465-471.

Li Q，Li Y，Li Z W，Wu C R.2006.A seismic nonlinear inversion for apparent fracture density of hydrocarbon reservoir based on GA-BP theory［J］.Progress in Geophysics，21（2）：465-471（in Chinese）.

李文林，李承楚.1994.一种拾取裂隙密度的新方法［J］.石油地球物理勘探，29（3）：286-293.

Li W L，Li C C.1994.A new method for determing fracture density［J］.Oil Geophysical Prospecting，29（3）：286-293（in Chinese）.

田锋.2007.裂缝体的弹性模量和裂缝密度［J］.地质学报，81（10）：1338-1344.

Tian F.2007.Elasticity modulus and fracture density in the rocks with crevices［J］.Acta Geologica Sinica，81（10）：1338-1344（in Chinese）.

魏建新.2002.不同裂缝密度的物理模型研究［J］.石油物探，41（4）：433-438.

Wei J X.2002.A physical model study of different cracks densities［J］.Geophysical Prospecting for Petroleum，41（4）：433-438（in Chinese）.

熊翥.2005.西部碳酸盐岩储层识别和特征描述［J］.勘探地球物理进展，28（4）：229-233.

Xiong Z.2005.Identification and characterization of carbonate reservoirs in western China［J］.Progress in Exploration Geophysics，28（4）：229-233（in Chinese）.

朱培民，王家映，於文辉，朱光明.2001.用纵波AVO数据反演储层裂隙密度参数［J］.石油物探，40（2）：1-12.

Zhu P M，Wang J Y，Yu W H，Zhu G M.2011.Inverting reservoir fracture density using P-wave AVO data［J］.Geophysical Prospecting for Petroleum，40（2）：1-12（in Chinese）.

Bakulin A，Grechka V，Tsvankin L.2000a.Estimation of fracture parameters from reflection seismic data.Part 1：HTI model due to a single fracture set［J］.Geophysics，65（6）：1788-1802.

Bakulin A，Grechka V，Tsvankin L.2000b.Estimation of fracture parameters from reflection seismic data.Part 2：Fractured models with orthorhombic symmetry［J］.Geophysics，65（6）：1803-1817.

Bakulin A，Grechka V，Tsvankin L.2000c.Estimation of fracture parameters from reflection seismic data.Part 3：Fractured models with monoclinic symmetry［J］.Geophysics，65（6）：1818-1830.

Chichinina T，Sabinin V，Ronquillo-Jarillo G.2006.QVOA analysis：P-wave attenuation anisotropy for fracture characterization［J］.Geophysics，71（3）：C37-C48.

Crampin S.1985.Evidence for aligned cracks in the Earth’s crust［J］.First Break，3（3）：12-15.

Downton J，Gray D.2006.AVAZ parameter uncertainty estimation［C］∥Expanded Abstracts of76th Annual Internat SEG Mtg.Houston：New Orleans：Society of Exploration Gephysicist：234-238.

Hudson J A，Liu E.1999.Effective elastic properties of heavily faulted structures［J］.Geophysics，64（4）：479-485.

Kangan F，James B R.1996.The relationship of Thomsen anisotropic parameters to the crack density of a cracked medium［J］.CREWES Research Report，8（6）：1-7.

Lynn H B，Bates C R，Simon M K，Van Dok R.1995.The effects of azimuthal anisotropy in P-wave 3-D seismic［C］∥Expanded Abstracts of65th Annual Internat SEG Mtg.Houston：Society of Exploration Gephysicist：293-296.

Lynn H B，Simon M K，Bates C R.1996a.Correlation between P-wave AVOA and S-wave traveltime anisotropy in naturally fractured gas reservoir［J］.The Leading Edge，15（8）：931-935.

Lynn H B，Simon M K，Bates C R，Van Dok R.1996b.Azimuthal anisotropy in P-wave 3-D（multiazimuth）seismic data［J］.The Leading Edge，15（8）：923-928.

Mallick S，Frazer L N.1991.Reflection／transmission coefficients and azimuthal anisotropy in marine seismic studies［J］.Geophys J Int，105（1）：241-252.

Mavko G，Mukerji T，Dvorkin J.1998.The Rock Physics Handbook：Tools for Seismic Analysis of Porous Media［M］.Cambridge：Cambridge University Press：133-139.

Reid F，MacBeth C.2000.Tests of vector fidelity in permanently installed multi-component sensors［C］∥Expanded Abstracts of70th Annual Internat SEG Mtg.Calgary：Society of Exploration Gephysicist：1213-1216.

Rüger A.1997.P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry［J］.Geophysics，62（3）：713-722.

Sayers C M.1988.Stress-induced ultrasonic wave velocity anisotropy in fractured rock［J］.Ultrasonics，26（6）：311-317.

Thomsen L.1986.Weak elastic anisotropy［J］.Geophysics，51（10）：1954-1966.

Tsvankin I，Thomsen L.1995.Inversion of reflection traveltimes for transverse isotropy［J］.Geophysics，60（4）：1095-1107.