基于纵横波同步联合的孔隙模量三参数AVO 反演方法

2020-02-09杜启振刘力辉符力耘冷雪梅刘子煊

岩性油气藏 2020年1期

丁燕，杜启振，刘力辉，符力耘，冷雪梅，刘子煊

（1.中国石油大学（华东）深层油气重点实验室，山东青岛 266580；2.青岛海洋科学与技术国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛 266580；3.北京诺克斯达石油科技有限公司，北京 100083；4.中国石化国际勘探开发有限公司，北京 100029）

0 引言

近几年，随着叠前地震反演技术和多波多分量技术的发展，储层流体识别方法的有效性和准确度不断提高。利用叠前振幅随偏移距或角度变化信息来进行流体识别是一种重要的方法［1-4］。印兴耀等［5-6］通过含有Gassmann 流体项的弹性阻抗方程直接反演来提高流体识别的准确性；Zong 等［7］建立了纵波模量和横波模量表示的孔隙弹性理论和反射系数近似方程，并发展了基于纵波模量和横波模量AVO 反演的流体识别方法；张世鑫［8］发展了基于固液解耦的地震流体识别方法，克服了常规流体识别方法受孔隙度影响产生的流体识别假象，然而，这些方法仅使用了反射纵波地震数据，在岩性和流体的共同影响下，利用纵波得到的反演结果往往存在精度低、反演不稳定、边界刻画不清晰等问题。

与纵波速度相比，横波速度受流体影响较小，因此在反演中同时利用纵波和转换横波地震数据，可以降低反演的多解性，提高储层预测精度［9-11］。Larsen 等［12］提出了纵波和转换波AVO 联合加权叠加反演方法；陈天胜等［13］基于一种方向加速度优化算法，利用纵波和转换横波反演得到了稳定的速度比值；Du 等［14］以弹性孔隙介质理论为基础，将流体因子引入Russell 反射系数公式，实现了多波联合AVO 反演，但常规流体因子由于受孔隙度和孔隙流体的耦合作用影响，流体识别存在的多解性问题仍未得到有效解决。

综合前人的研究结果，首先讨论Gassmann 流体因子和流体体积模量随孔隙度的变化关系，然后针对复杂储层中Gassmann 流体因子受孔隙流体和孔隙度耦合影响而存在的流体识别假象问题，推导基于流体体积模量、剪切模量和孔隙度的PP 波和PS 波Zoeppritz 线性近似方程，建立一种能够稳定获取等效流体体积模量、剪切模量和孔隙度的纵横波同步联合反演方法，以期解决单纵波反演精度低、反演不稳定、边界刻画不清晰的问题。

1 理论方法

1.1 孔隙流体体积模量的敏感度分析

Domenico［15］假设流体的体积模量是气体和液体体积模量的体积加权函数

式中：Kf为流体体积模量，GPa；Sw为含水饱和度，%；Kw为液体体积模量，GPa；Kg为气体体积模量，GPa。

Han 等［16］提出了Gassmann 流体项经验公式

式中：f为流体因子项为岩石骨架矿物与孔隙度的综合增益函数；Φ为孔隙度，%；Kdry为干岩石体积模量，GPa；Km为固体矿物基质体积模量，GPa；Kn为干岩石体积模量与基质体积模量的比值。

Nur 等［17］给出了临界孔隙度模型

式中：μdry为干岩石剪切模量，GPa；μm为固体基质剪切模量（GPa）；Φc为临界孔隙度，%。对增益函数G(Φ)进行简化，得到

假设固体骨架矿物相同，孔隙中包含水和气，则Gassmann 流体因子和流体体积模量随含水饱和度和孔隙度的变化趋势（图1）便可由式（1）与式（2）计算得到。由图1 可看出，Gassmann 流体因子受孔隙流体和孔隙度的综合影响，即与孔隙度和含水饱和度呈非线性关系［图1（a）］，而流体体积模量在一定孔隙度下随含水饱和度的变化趋势则为完全线性［图1（b）］。因此，如果将流体体积模量与孔隙度分别单独作为变量，即可消除常规流体因子受孔隙度和孔隙流体的耦合影响，解决流体识别多解性的问题。

图1 流体因子随孔隙度与含水饱和度的变化趋势Fig.1 Trend of fluid factor with porosity and water saturation

1.2 Kf−μ −Φ 纵横波孔隙模量三参数AVO 近似方程和反演

1.2.1Kf−μ−Φ纵横波孔隙模量三参数AVO 近似方程推导

Russell 等［18］建立了弹性孔隙介质理论和AVO技术之间的关系，推导了关于流体因子、剪切模量和密度的三参数AVO 近似式

式中：θ为PP 波的入射角与透射角的平均值，（°）；φ为PS 波反射角与透射角的平均值，（°）；μ为剪切模量，GPa；ρ为密度，g/cm3；γsat为饱和岩石条件下纵横波速度比均值；γdry为干岩石情况下的纵横波速度比。

Gassmann 流体项及剪切模量表达式为

式中：vP，vS分别为饱和岩石纵波、横波速度，m/s。由式（7）和式（8）的流体因子、剪切模量变化率、纵横波速度变化率与密度变化率的关系，可以得到

根据Gardner 等［19］提出的经验关系式，得到纵波速度与密度变化率之间的关系

结合式（9）—（11），可以推导得到

将式（12）分别代入式（5）和式（6），并可将其简化得到

根据Han 等［16］总结得到的流体体积模量Kf与Gassmann 流体项f之间的关系，联立式（2）和式（4）可以得到

将式（15）代入式（13）和式（14），可以得到

式（16）和式（17）即为推导的PP 波和PS 波三参数Zoeppritz 线性近似公式。

1.2.2 近似方程精度分析

为验证推导方程的精度［20］，分别设计了3 种砂岩模型（表1），分别用精确Zoeppritz 方程［21］，Russell方程［18］与式（16），式（17）计算得到地层界面处PP波、PS 波的反射系数，其中=2.333。当入射角不超过30°时，基于Kf−μ−Φ的纵波近似式与Russell近似式精度相近，而转换波近似方程的精度则略高于Russell 近似式精度［图2（a）—（b），图2（e）—（f）］。基于Kf−μ−Φ的纵波近似式与Russell 近似精度相近，而转换波近似方程精度在大于30° 时变差［图2（c）—（d）］。综合以上分析，新推导的孔隙模量三参数AVO 近似方程在入射角0°～30°内具有较高的精度，因此在最佳入射角0°～30°时利用本文提出的Kf−μ−Φ方程进行纵横波同步联合反演具有可行性。

表1 砂岩模型参数Table 1 Sandstone model parameters

图2 3 种模型的PP 波、PS 波反射系数对比Fig.2 Comparison of PP wave and PS wave reflection coefficient of three kinds of models

1.3 基于贝叶斯理论的叠前同步联合反演

对PP 波和PS 波进行同步联合反演，考虑M层n个入射角的情况，将式（16）和式（17）写成矩阵形式

根据地震褶积模型理论，增加子波矩阵至式（18），则叠前地震道集可以通过矩阵乘法表示为

式中：WP和WS分别为PP 波和PS 波的子波矩阵。

式（19）可以写成

式中：d为PP 波和PS 波地震记录矩阵；G为含有子波的系数矩阵；r为待反演的由等效体积模量、剪切模量和孔隙度组成的列向量系数矩阵，即

为提高反演的稳定性，本文采用基于贝叶斯理论［22］的PP 波和PS 波同步联合反演。贝叶斯反演的目的就是在给定测量数据d（带噪声n）的情况下估计模型参数r，把似然函数赋给噪音分布函数

将（24）式代入到边缘化公式，取对数并进一步求导得到目标函数为

其中

令式（25）等于0，就得到问题的贝叶斯解

式中：μ为数据保真项与模型约束项之间权重的调节因子。

反演式（26）存在一定的弱非线性，为得到稳定的求解结果，采用共轭梯度法进行求解。

2 应用实例

2.1 Marmousi2 模型试算

为验证本方法的正确性，本文利用Marmousi2油藏模型［23］进行测试。图3 为利用纵、横波速度，密度以及孔隙度等参数，计算得到模型等效流体体积模量、剪切模量和孔隙度（Kf−μ−Φ）剖面。在此模型参数基础上，通过式（16）和式（17）正演得到反射系数，并与30 Hz 雷克子波褶积生成纵波和转换波角度道集，同时加入信噪比为2 的高斯噪音。图4是CDP 300 处PP 波与PS 波的角道集，其入射角为0°～30°。

图3 等效流体体积模量（a）、剪切模量（b）和孔隙度（c）模型Fig.3 Models of equivalent fluid bulk modulus（a），shear modulus（b）and porosity（c）

图4 CDP 300 处PP 波角道集（a）与PS 波角道集（b）Fig.4 PP wave angle gathers（a）and PS wave angle gathers（b）at CDP 300

利用中心角度为5°，15°，25°部分的叠加数据，分别进行单纵波和纵横波同步联合反演。由图5中CDP 300 处单纵波、纵横波同步联合反演与实际模型单道对比结果可以得知，本文方法的反演结果误差相对较小，稳定性高于单纵波反演结果；图6为单纵波和纵横波同步联合反演结果对比，其中，等效流体体积模量剖面能清楚地反映含油气砂岩特征，剪切模量剖面反映了砂岩与页岩的岩性差别，孔隙度剖面则准确反映了储层的孔隙特征。综合对比结果，采用本文纵横波同步联合反演方法得到的结果较单纵波反演结果具有更高的精确度，且能准确刻画地质体储层特征和含气特征。

图5 CDP 300 处单纵波、纵横波同步联合反演结果与实际模型单道对比Fig.5 Single channel comparison of single P-wave and PP wave and SS wave with actual model at CDP 300

图6 单纵波和纵横波联合反演对比Fig.6 Comparison of single P-wave and PP-wave and SS wave joint inversion

2.2 实际数据应用

利用中国LJ 地区A 油田M 测线的纵波和转换横波数据进行应用效果测试。该区位于凹陷南缓斜坡，受燕山运动和喜马拉雅运动的影响，形成了向北倾没的大型鼻状构造，受构造控制，该区油气藏类型主要为断块油气藏，兼有岩性、泥岩裂缝等特殊油气藏［24］。本区纵波地震资料主频为35 Hz，转换波地震资料主频为15 Hz。由于实际纵波和转换横波记录的反射振幅相差较大，为得到更可靠的联合反演结果，需要对纵波振幅和转换横波振幅做匹配。以纵波振幅为标准，每一角度计算一个校正系数［式（27）］，为了不改变整个角道集的AVO 特征，将各个角度的校正系数求平均作为最终的校正系数，图7 为在时间域内匹配一致的纵波和转换波角道集。

图7 CDP 1000 点PP 波角道集（a）与PS 波角道集（b）Fig.7 PP wave angle gathers（a）and PS wave angle gathers（b）at CDP 1000

利用反演得到的M 测线目标层段局部放大的等效流体体积模量［图8（a）］、剪切模量［图8（b）］和孔隙度剖面［图8（c）］。图中黑色线柱表示生产井段内的油层，紫色线柱表示生产井段内的水层，经过对比发现，油层对应等效流体体积模量剖面低值［图8（a）中红色所示］，孔隙度剖面高值［图8（c）中红色所示］，而水层等效流体体积模量值高于油层［图8（a）中蓝色所示］，相对应的孔隙度剖面低值［图8（c）中蓝色所示］。油水层所在砂岩与泥页岩区在剪切模量剖面上有较明显的区分［图8（b）］。综上所述，反演结果中储层特征清晰，油、水层得到了较好的区分，且油层范围与已知生产井的信息较吻合。