程晓艳，胡曦，罗颖，李亚丁，高培丞，谢伟，戴赟，佟恺林，张奎，王群武

（1.四川页岩气勘探开发责任有限公司，四川 成都 610051；2.川庆钻探工程有限公司地质勘探开发研究院，四川 成都 610051；3.北京普瑞斯安能源科技有限公司，北京 100085）

0 引言

裂缝作为页岩储层中重要的储集空间，是油气勘探开发中的一项重要的预测刻画对象。针对裂缝预测研究，吸引了大量学者的关注［1-4］。在页岩地层中发育的裂缝，形成条件较为复杂，发育程度各异，使得裂缝的预测和识别较为困难。利用地震数据进行裂缝预测主要分为叠前和叠后预测2 大类。叠前预测主要利用裂缝方位各向异性在方位地震数据中形成的差异特征来预测裂缝发育程度、分布及裂缝密度等，例如基于方位地震的AVAZ 裂缝密度反演方法［5-8］。叠后预测主要利用裂缝在三维叠后地震属性中的特征来对裂缝进行识别和预测，例如蚂蚁追踪、曲率属性法、相干体、分形维数法等［9-11］。

地震相干技术作为一种三维地震解释技术，目前常用来进行裂缝的识别。到目前为止，应用比较广泛的相干算法是基于地震道互相干的C1 算法［12］、地震多道相似性的C2 算法［13］、基于结构相似性的C3 算法［14］等。以上方法的原理都是基于相邻地震道的不连续性，这些算法可以有效突出地震道之间不连续性的特征，可以更有效地进行裂缝及断层识别。除此之外，Bakker等［15-16］还提出了通过梯度结构张量方法进行相干体计算的方法。相干体属性通过计算地震道之间的波形相似性，可以将地震同相轴连续性较差的区域标识出来。目前，相干体技术广泛应用于断层、裂缝等特殊地质体的识别［17-45］。

页岩裂缝的发育方向、发育程度和展布特征是由构造运动、应力场及岩石的力学性质等决定的。由于裂缝的存在，地震响应特征也变得复杂，使得地震波传播速度及振幅发生变化，而且会出现横波双折射现象。地震波在各向异性介质中传播时，会出现横波分裂现象，Crampin 等［17-20］提出了裂缝介质横波分裂理论。除此之外，在裂缝介质中，由于方位地震波在经过各向异性介质时，会产生方位差异，裂缝探测可以利用这种方位的变化特征来实现。Thomsen［21］推导了VTI 介质的弹性模量计算公式，并得到了表征各向异性程度的Thomsen 参数。Tsvankin 等［22］将上述理论推广应用到HTI 中。Hudson 模型［23］是著名的等效介质理论之一，描述了包含薄而扁平的椭圆形裂缝背景介质的整体各向异性特征。根据弱各向异性理论假设，Rüger［24］基于各向异性介质岩石物理理论，推导了考虑方位变化的反射系数方程。Gray 等［25］提出了基于HTI 介质的裂缝方位和裂缝密度的AVAZ 计算方法。Schoenberg 等［26］基于线性滑动理论，建立了弱各向异性假设下的裂缝参数模型，裂缝密度由这些参数来表征。Bachrach 等［27-28］通过重构裂缝储层弹性和各向异性参数之间的关系，实现裂缝参数的地震反演。

本文提出一种考虑相干属性约束的页岩裂缝地震预测方法，在考虑地震属性相干作用的基础上，利用叠前地震数据反演裂缝参数，减小了裂缝参数反演的多解性，提高了反演结果的准确性。首先，通过空间变换离心窗口提取地震相干属性，提高相干体属性在噪声处的异常低相干值，改善裂缝边界的计算效果。其次，概率化相干属性结果，形成裂缝发生概率分布体。最后，通过相干属性体的概率化模型作为裂缝参数的先验分布，结合贝叶斯概率化反演算法，开展不同方位角下的各向异性弹性阻抗反演，以实现裂缝参数的后验概率分布估计及反演预测，进而提取工区的裂缝岩石物理参数，为储层裂缝预测和描述提供可靠的支撑。

1 方法原理

1.1 基于离心窗提取相干属性及裂缝密度概率化

相干体是一个检测波形相似性的属性，通过计算地震道之间的波形相似度来识别地震同相轴连续性较差的区域，有利于识别断裂和特殊地质体。本文将基于离心窗的扫描方式提取相干属性，在每个窗口进行相干计算时，采用基于本征结构和倾角扫描的第3 代相干算法，循环多个窗口之后，最终选取相干值最大窗口所得的计算结果。计算流程如图1 所示。

图1 相干属性提取流程Fig.1 Process of coherent attribute extraction

第3 代相干算法作为本文使用的方法，可以按窗扫描地震数据的倾角或方位来进行，是一种基于本征结构的模式。首先，选取扫描窗中的地震数据，然后，以扫描窗为中心点建立样点矢量，由样点矢量的列在倾角p 和q 的平面上，进行相关计算。协方差矩阵为

式中：Cij（i，j 分别为1～9 的整数）为协方差矩阵元素。

Cij为计算窗口内垂向2 个地震波形的互相关。

式中：t 为样点的时间，s；xi，yi分别为第i 个样点与中心点之间距离的x，y 轴的投影值；ui为第i 个样点对应的振幅；u 为地震振幅；＜u＞为所有样点的均值振幅；k 为上下样点数。

对协方差矩阵进行特征值和特征向量分解，可以写成：

式中：ν 为特征向量；λ 为特征向量对应的特征值；上标m 为特征值对应的特征向量序号。

通过协方差矩阵的分解，本征值及对应的本征向量可以由计算获取。矩阵的主要变化量可以通过第1本征值和对应的向量表征，而剩余成分则通过第2 本征值代表。理论上讲，数据体95%的信息量一般只需少数本征值和向量即可代表。由于第1 本征值和矩阵的迹就能用来表征数据的主要信息量，所以第3 代相干算法也采用类似的原则，则基于本征结构的相干性估计可定义为

式中：Ce为基于本征结构的相干属性；λ1为第1 特征值；Cjj为协方差矩阵的对角矩阵元素。

相干体属性与裂缝密切相关，通过计算地震道之间的波形相似度来识别地震同相轴连续性较差的区域，有利于识别断裂和特殊地质体。本文将利用得到的相干体数据对裂缝密度反演进行约束，一般情况下，相干值越小表明存在裂缝的可能性越大，那么可以得到整个三维体发生裂缝的概率模型PFra（f） 等价为

式中：f 为相干值；fmax为相干最大值。

1.2 页岩储层裂缝密度方位弹性阻抗反演方法

裂缝中流体不同，Thomsen 各向异性参数与线性滑动模型参数和裂缝密度的关系也是有区别的［29］。当气体填充于裂缝当中，线性滑动模型参数ΔN，ΔT的表达式分别为

式中：e 为裂缝密度；g 为裂缝岩石纵、横波速度比倒数的平方。

当裂缝中填充的是油、水时，ΔN，ΔT可以表示为

因此，裂缝密度为

利用式（10），当地震数据反演出裂缝岩石弹性参数和线性滑动模型参数时，就可以计算出裂缝密度。

标准化的垂直裂缝方位弹性阻抗近似式，作为页岩储层裂缝密度地震预测的桥梁，有效地表征了方位弹性阻抗EI 和裂缝参数的关系：

其中：a（θ）=sec2θ；b（θ）=-8gsin2θ；c（θ，ϕ）=-2·（gcos2ϕsin2θ）（1-2g）；d（θ，ϕ）=2gcos2ϕsin2θ。

式中：IP，IS分别为纵、横波波阻抗，kg/m3·m/s；IP0，IS0分别为纵、横波波阻抗均值，kg/m3·m/s；θ 为地震入射角，（°）；ϕ 为方位角，（°）。

对式（11）两端取对数得：

式中：EI0为弹性阻抗均值，kg/m3·m/s。

求解裂缝岩石物理参数，由式（12）可以看出，需要4 个方位以上的地震数据（见式（13）。下标1，2，3，4 表示4 个不同方位，符号意义与上相同）。

利用式（13）建立裂缝岩石参数与方位弹性阻抗的确定性关系。将其简写成矩阵形式为

式中：d 为对应式（13）等式左边的方位弹性阻抗；G 为系数矩阵；m 为待求解岩石物理参数。

贝叶斯核心思想主要是利用似然函数和先验分布来获取后验概率。根据贝叶斯理论，加上裂缝发生概率模型的约束，本文的待反演参数的后验概率分布P 可以写成：

由地震相干值转换成的裂缝发生概率模型作为约束项加入到反演目标函数中，则目标函数J 可以写成相干约束和模型约束的形式：

式中：η1，η2为反演权重系数；Ml为参数模型。

则裂缝参数的最终目标泛函m* 可以表示为

综上所述，建立了基于相干约束的页岩裂缝OVT地震预测方法研究及应用流程，具体流程见图2。

图2 页岩裂缝OVT 地震预测方法研究及应用流程Fig.2 Research and application process of OVT seismic prediction method for shale fractures

2 实际应用

中国西南某盆地页岩经历了多期的构造演化和长期的变质作用，其油气成藏条件、成藏过程、分布与富集规律均十分复杂。该盆地页岩油气藏的裂缝系统具有多尺度、组合形式复杂、非均质性强等特征，单一技术很难实现裂缝的准确识别。因此，提出基于相干约束的页岩裂缝OVT 地震预测方法对研究工区的页岩进行应用。把工区内0°～120°方位角地震数据等分为4个数据体（见图3）。图4 为利用图3 中4 个方位地震数据反演得到的不同方位弹性阻抗。由图3，4 看出，不同方位弹性阻抗之间有一定的差异，可能是由于裂缝存在引起的方位差异。

图3 不同方位角的地震数据Fig.3 Seismic data of different azimuths

图5 为地震相干属性二维剖面及其等效的裂缝分布概率模型。由图5 看出，在相干值较小的地方，裂缝发育的概率大。图6 为不考虑相干约束、只考虑低频约束的4 个参数的反演结果。

图5 地震相干属性二维剖面及其等效裂缝分布概率模型Fig.5 Two-dimensional section of seismic coherence attributes and its equivalent fracture distribution probability model

图6 不考虑相干约束的裂缝参数反演结果Fig.6 Inversion results of fracture parameters without considering coherence constraints

图7 为考虑相干约束的反演结果。对比图6 和图7 的参数反演结果可以发现，考虑相干约束后，反演结果在横向与纵向分辨率都有了明显提升，有效挖掘了地震的非平稳特征，使得反演结果更加合理。

图6a 和图7a 分别是常规方法反演的纵波阻抗和加相干约束的纵波阻抗，由于相干作用的约束，图7a的纵波阻抗预测结果融入了图5a 中相干值的部分细节，使得相干值低的地方，也就是在同相轴非连续的地方预测的纵波阻抗会有一定的变化。图7 中其他参数，包括横波阻抗和各向异性参数都有同样的效果。而常规方法中的参数预测结果缺乏这种效果，更多的是因为常规方法更加依赖低频模型的约束，导致预测的结果更趋于平滑，缺少了细节上的预测。而在地震同相轴非连续的地方，往往参数也会有所差异，正是利用了相干约束，图7 的预测结果才更符合预期。

图8 为不考虑相干约束和考虑相干约束的裂缝密度反演结果对比，井上投的是裂缝解释结果，红色部分为裂缝发育段。

图8 裂缝密度反演结果对比Fig.8 Comparison of fracture density inversion results

图8a 是基于常规方法得到的裂缝密度结果，红圈内发育的3 处裂缝均没有较好吻合。图8b 是考虑相干约束的裂缝密度反演结果，在红圈内底部，反演效果有了明显的提升，在红圈上部也有一定改善，与测井解释的裂缝发育处能够较好吻合。黄圈部分为考虑相干约束的裂缝密度反演结果（见图8b），分辨率得到了有效提升。

图9 为裂缝密度的预测结果切片。其中，图9a 为不考虑相干约束的常规裂缝预测结果，图9b 为考虑相干约束的裂缝密度预测结果。由图9a、图9b 对比结果来看，考虑相干约束的裂缝密度预测细节更丰富，分辨率更高，在井周的切片显示与解释结果对应。

3 结束语

本文建立了基于相干约束的地震数据裂缝参数反演方法，有效地提升了裂缝参数预测的准确性。通过相干约束，改善了常规裂缝参数反演对于初始模型的依赖性，提高了地震数据对裂缝参数在纵向及横向上的分辨率，使得裂缝参数的预测精度有了明显的提高。本方法也是利用属性约束对裂缝参数反演的一种尝试，使得反演结果相对于常规的方法更具合理性。