基于方位各向异性反演的裂缝型储层预测及流体识别方法

doi:10.11720/wtyht.2025.0121

基于方位各向异性反演的裂缝型储层预测及流体识别方法

吴怡^,¹, 周长所¹, 徐国贤¹, 袁俊亮¹, 宋晓麟², 曾勇坚^,², 王群武², 张奎²

1.中海油研究总院责任有限公司,北京 100015

2.北京普瑞斯安能源科技有限公司,北京 100015

A method for fracture density prediction and fluid identification of fractured reservoirs based on azimuthal anisotropic inversion

WU Yi^,¹, ZHOU Chang-Suo¹, XU Guo-Xian¹, YUAN Jun-Liang¹, SONG Xiao-Lin², ZENG Yong-Jian^,², WANG Qun-Wu², ZHANG Kui²

1. Research Institute Co.,Ltd.,CNOOC,Beijing 100015,China

2. Beijing Precise Energy Technology Co.,Ltd.,Beijing 100015,China

通讯作者: 曾勇坚(1991-),男,江西临川人,高级工程师,主要从事叠前地震预测工作。Email:zengyj91@163.com

第一作者: 吴怡(1987-),男,高级工程师,2013年毕业于中国石油大学(北京),就职于中海油研究总院有限责任公司,主要从事海上石油天然气勘探与开发等工作。Email:wuyi11@cnooc.com.cn

责任编辑: 叶佩

收稿日期: 2025-04-15 修回日期: 2025-07-29

基金资助:

中国海洋石油集团有限公司“永乐8区古潜山钻井地质风险分析研究”项目(CCL2024RCPS0396PSN)

Received: 2025-04-15 Revised: 2025-07-29

摘要

基于方位地震数据的叠前各向异性地震反演作为裂缝型储层裂缝检测与流体识别的核心方法,为裂缝型储层的勘探开发提供了有效的指导作用。然而,随着勘探目标精细刻画要求的提升,常规约束下的各向异性反演方法的稳定性与可靠性受到挑战,且在发育高角度裂缝的中国东部裂缝型储层工区,缺乏有针对性的裂缝密度预测与流体识别的直接反演方法。为此,本文首先利用多井交会明确了用于流体识别的敏感因子;之后以孔隙弹性理论为指导,基于线性化HTI介质反射系数方程,推导得到了能够反映流体因子与裂缝密度随入射角及方位角变化的方位地震反射系数方程;最后,依托McMC优化算法,发展了一种L₁范数约束的贝叶斯各向异性两步地震反演方法,以通过地震反演的手段实现了高角度裂缝的储层流体识别及裂缝密度的准确预测。反射系数方程精度分析证明了所推导方程的可行性,通过合成记录测试与实际工区应用,验证了所提出方法的合理性与可靠性,为发育高角度裂缝的储层裂缝密度预测与流体识别提供了新的方案。

关键词： 裂缝型储层; 各向异性反演; 流体识别; 裂缝密度

Abstract

The pre-stack anisotropic seismic inversion based on azimuthal seismic data is a key method for fracture detection and fluid identification of fractured reservoirs,which provides effective guidance for exploration and development of fractured reservoirs.However,the anisotropic inversion method under conventional constraints encounters challenges in terms of stability and reliability,with higher requirements being placed for fine characterization of exploration targets.Moreover,there exists a lack of direct inversion methods targeting the fracture density prediction and fluid identification of the fractured reservoirs with high-angle fractures in East China.Therefore,this paper first identified the sensitive factors for fluid identification using multi-well crossplots.Then,based on the poroelasticity theory and the linearized reflection coefficient equation for horizontal transverse isotropy(HTI) media,an azimuthal seismic reflection coefficient equation was deduced,which can reflect the variations of fluid factors and fracture density with the angle of incidence and azimuth angle.Finally,based on the advanced Markov Chain Monte Carlo(MCMC) algorithm,a two-step Bayesian anisotropic seismic inversion method with the L₁ norm constraint was proposed,achieving the fluid identification and accurate prediction of fracture density for reservoirs with high-angle fractures.The deduced reflection coefficient equation proved to be feasible through accuracy analysis.In addition,the method proved to be rational and reliable through synthetic seismogram testing and application in actual survey areas,offering a novel solution for the fracture density prediction and fluid identification in such reservoirs.

Keywords： fractured reservoir; anisotropic inversion; fluid identification; fracture density

PDF (10779KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

吴怡, 周长所, 徐国贤, 袁俊亮, 宋晓麟, 曾勇坚, 王群武, 张奎. 基于方位各向异性反演的裂缝型储层预测及流体识别方法[J]. 物探与化探, 2025, 49(5): 1173-1189 doi:10.11720/wtyht.2025.0121

WU Yi, ZHOU Chang-Suo, XU Guo-Xian, YUAN Jun-Liang, SONG Xiao-Lin, ZENG Yong-Jian, WANG Qun-Wu, ZHANG Kui. A method for fracture density prediction and fluid identification of fractured reservoirs based on azimuthal anisotropic inversion[J]. Geophysical and Geochemical Exploration, 2025, 49(5): 1173-1189 doi:10.11720/wtyht.2025.0121

0 引言

AVAZ(amplitude variation with incident angle and azimuth)理论是振幅随入射角变化(amplitude variation with offset,AVO)理论向方位地震领域的拓展,在考虑振幅随偏移距变化的同时还兼顾振幅随方位的变化规律,进而通过地震数据来分析储层特征^[1-2]。为充分挖掘五维数据中的地下地质信息,进而对复杂非常规油气储层的勘探开发提供数据参考,基于AVAZ理论的方位地震反演技术应运而生。通过构建方位各向异性地震数据与储层参数之间的关系,从而由方位地震数据获取裂缝型非常规储层弹性参数、物性参数以及裂缝参数等用于储层表征与裂缝检测的有用信息^[3-8]。

实际储层条件下,稀疏垂直裂缝体系引发的弱水平横向各向异性(horizontal transverse isotropy,HTI)特征,易导致方位振幅差异数据的信噪比衰减,形成方法应用的主要技术瓶颈。因此,如何提高方位地震数据各向异性反演的稳定性与抗噪性是目前非常规油气勘探需要解决的重要目标。基于方位振幅差异的各向异性两步地震反演方法的技术原理在于利用方位各向异性反射系数的线性近似特性,通过差异方位数据组的差分运算消除各向同性分量,从而有效降低反演体系参数维度,理论上能够显著改善解的稳定性^[9-13]。

地下介质流体识别是勘探地球物理学的热点领域,也是油气地球物理勘探的最终目标。近年来,利用地震资料的振幅、信息进行油气识别的方法得到了广泛的发展。泊松比作为表征储层流体信息的参数之一,常用于裂缝型储层的流体识别^[14-15]。然而,由于泊松比耦合了多种岩石骨架特征的信息,使得其难以普适于大多数裂缝型储层工区。因此如何利用方位地震数据,提高复杂各向异性介质的油气识别的准确度,成为了新难题。为此,也有学者将基于Gassmann理论获得的Russell流体因子f用于HTI介质中的流体识别中。Pan等^[16]基于Gassmann 各向异性孔隙弹性理论和线性滑动模型,推导了含基于Russell流体因子的HTI介质反射系数方程,并通过AVO反演得到了较好的油气识别效果。Feng等^[17]在Russell流体因子上进一步考虑了喷射流机制,实现了考虑喷射流作用的HTI介质流体识别。Zuo等^[18]在水平横向各向同性介质假设下,重新推导了包含剪切液相参数和裂缝参数的各向异性反射系数方程,有效降低了方程的维度和病态性。同时,引入了自适应独立粘性McMC(Markov chain Monte Carlo)策略并结合吉布斯算法,构建了基于贝叶斯框架的改进算法,提升了反演的精度和稳定性,该方法在中国西南地区页岩气勘探区裂缝预测和流体识别方面具有良好的应用价值。然而此类各向异性反演方法虽然提高了流体识别精度,但是由于其方程中大多利用了裂缝参数来表征裂缝发育程度,难以在获得可靠的流体识别结果的同时定量或直接地给出裂缝密度信息。为解决裂缝密度与流体同时反演的问题,严彬等^[19] 基于弱化干岩石骨架作用影响的固液缝解耦流体因子开展了TTI(titled transversely isotropy)介质下的流体因子与裂缝密度同时反演,为裂缝型储层流体识别提供了有效的解决方案。然而该方法所用到的弹性阻抗反演手段相比于AVO反演将首先获取地层的弹性信息,进而再估计地层的流体信息与裂缝信息,在裂缝型储层参数的提取过程中容易产生累积误差。

面对上述研究现状,本研究首先针对发育高角度裂缝的中国东部某裂缝型储层裂缝密度定量预测及流体识别的需要,基于多井交会,优选能够区分含气储层与含水层的流体因子f。在此基础上,考虑到工区的裂缝型储层发育高角度裂缝,基于HTI介质线性化PP波反射系数方程,分别利用多孔介质理论中的流体因子与弹性模量以及裂缝弱度参数与裂缝密度参数之间的定量关系,得到了同时含有裂缝密度与流体因子f的新反射系数方程。为降低多参数反演的不稳定性以及提高反演结果的可靠性,我们基于贝叶斯框架与零阶熵正则化约束建立了反演的目标函数,并采用方位作差的两步AVO反演策略进行裂缝密度以及流体因子的提取。合成记录测试证明了所提出方法的可行性,该方法在实际高角度裂缝发育的储层工区取得了合理、可靠的裂缝密度与流体识别结果。

1 理论与方法

1.1 流体识别因子交会分析及优选

实际裂缝型储层工区储层流体识别因子的识别精度对于评价储层含气性,获取含气储层的地下分布范围与展布状态起到了至关重要的作用。因此,有必要针对实际工区选取合理、有效的储层流体识别因子。该工区位于中国东部渤海海域渤中背斜区域,研究面向的目标储层岩性以砂、泥岩为主,储层段发育高角度裂缝,储层流体以天然气为主。图1展示了常用于实际裂缝型储层工区流体识别的各向同性背景弹性参数(流体因子f、泊松比以及纵波阻抗)与纵波速度的多井测井曲线交会图。其中,流体因子f由f=(ρ_b $V_{p}^{2}$ )- $γ_{d r y}^{2}$ (ρ_b $V_{s}^{2}$ )计算,其中的V_p、V_s 、ρ_b以及 $γ_{d r y}^{2}$ 分别表示岩石的纵波速度、横波速度、密度以及干岩石骨架的纵横波速度比的平方。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 实际裂缝型储层工区各向同性背景弹性参数与纵波速度测井曲线交会

Fig.1 Cross plot of isotropic background elastic parameters and longitudinal wave velocity logging curves of actual fractured reservoir working area

交会展示了不同参数用于区分气层与水层的能力。可以看出,泊松比与纵波阻抗交会难以通过其高值异常或者低值异常将气层与水层有效区分,而在图1a中的椭圆所圈出的区域内,流体因子f能够利用其低值异常对气层和水层进行有效地区分。这一结果的产生依赖于流体因子f中包含了饱和背景岩石的流体信息,相比于泊松比和纵波阻抗,其排除了大部分由干岩石骨架信息带来的干扰,使其具有较高的流体识别能力。基于上述流体识别因子交会分析结果,本研究将优选流体因子f进行实际裂缝型储层的流体识别。

1.2 方位各向异性反射系数方程

针对发育高角度裂缝的裂缝型储层流体识别及裂缝发育程度评估的需要,本研究将采用发育平行垂直裂缝的HTI介质来解释裂缝型储层的岩石物理与地震波传播规律。经附录A中的详细推导过程,可得含裂缝密度与流体识别因子的反射系数方程如下:

(1)$\begin{array}{c}R_{\mathrm{pp}}^{\mathrm{HTI}}(\theta)=\left[\left(1-\frac{\gamma_{\mathrm{dry}}^{2}}{\gamma_{\mathrm{sat}}^{2}}\right) \frac{\sec ^{2} \theta}{4}\right] \frac{\Delta f}{f}+\left[\frac{\gamma_{\mathrm{dry}}^{2}}{4 \gamma_{\mathrm{sat}}^{2}} \sec ^{2} \theta-\frac{2}{\gamma_{\mathrm{sat}}^{2}} \sin ^{2} \theta\right] \times \frac{\Delta \mu_{\mathrm{b}}}{\mu_{\mathrm{b}}}+\left[\frac{1}{2}-\frac{\sec ^{2} \theta}{4}\right] \frac{\Delta \rho_{\mathrm{b}}}{\rho_{b}}+ \\\left\{\begin{array}{l}-\frac{\sec ^{2} \theta}{4}\left[1-\left(1-\frac{\lambda_{\mathrm{b}}}{M_{\mathrm{b}}}\right)\left(\sin ^{2} \theta \cos ^{2} \varphi+\cos ^{2} \theta\right)\right]^{2} \frac{4}{3\left(1-1 / \gamma_{\mathrm{dry}}^{2}\right) / \gamma_{\mathrm{dry}}^{2}} \\+\frac{1}{2}\left(1-\frac{\lambda_{\mathrm{b}}}{M_{\mathrm{b}}}\right) \sin ^{2} \theta \sin ^{2} \varphi\left(1-\tan ^{2} \theta \cos ^{2} \varphi\right) \frac{16}{3\left(3-2 / \gamma_{\mathrm{dry}}^{2}\right)} e\end{array}\right\} \Delta e\end{array}$,

作为参与反演的核心方程,式(1)给出了HTI介质下的裂缝密度与流体识别因子与PP波反射系数之间的定量联系,为采用方位地震数据进行裂缝型储层的裂缝密度估算与含油气性评价提供了理论支撑。

反射系数的精度与入射角敏感性和方位角敏感性是该方程能否获得合理、可靠的预测结果的关键。为此,我们设计了表1与表2所示的两套双层裂缝型储层模型进行反射系数方程的精度分析。表1与表2分别建立了含气砂岩和含水砂岩与泥岩的反射界面,能够模拟大多数工区内的地层反射特征,对于模拟实际地震波的反射系数具有指导意义。图2与图3分别展示了基于裂缝型含气砂岩和含水砂岩与泥岩的双层模型的反射系数方程(式(1)所示的含流体因子的HTI方程和式(5)所示的含弹性模量的HTI方程)精度对比及误差分析结果。观察方程随入射角和方位角的变化关系可以看出,无论是含气砂岩还是含水砂岩与泥岩的分界面处的反射系数,都反映出了具有入射角敏感性,还具有较好的方位角敏感性,反射系数能够随入射角与方位角变化而变化,为利用方位各向异性地震数据进行储层裂缝密度预测与流体识别奠定了基础。观察本文所提出的新方程与常规含弹性模量线性方程的误差绝对值不难发现,误差绝对值随着入射角度的增大而增大,但在0°~30°的入射角使用范围内,误差仍然控制在较小的范围内,误差分析进一步验证了方程的精度,为反演结果的准确性提供了理论支持。

表1 双层裂缝型含气砂岩与泥岩模型参数

Table 1 Model parameters of double-layer fractured gas-sand and mudstone

	M_b/ GPa	μ_b/ GPa	ρ_b/ (g·cm^-3)	f/GPa	δ_N	δ_T	e
含气砂岩	36.25	13.28	2.60	6.60	0.16	0.08	0.03
泥岩	40.75	14.22	2.70	9.00	0.05	0.03	0.01

新窗口打开| 下载CSV

表2 双层裂缝型含水砂岩与泥岩模型参数

Table 2 Model parameters of double-layer fractured water-sand and mudstone

	M_b/ GPa	μ_b/ GPa	ρ_b/ (g·cm^-3)	f/GPa	δ_N	δ_T	e
含水砂岩	37.72	13.28	2.62	8.20	0.16	0.08	0.03
泥岩	40.75	14.22	2.70	9.00	0.05	0.03	0.01

新窗口打开| 下载CSV

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 基于双层裂缝型含气砂岩与泥岩模型的反射系数方程精度对比及误差分析

Fig.2 Precision comparison and error analysis of reflection coefficient equations based on double-layer fractured gas-bearing sand and mudstone model

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 基于双层裂缝型含水砂岩与泥岩模型的反射系数方程精度对比及误差分析

Fig.3 Precision comparison and error analysis of reflection coefficient equations based on double-layer fractured water-bearing sand and mudstone model

1.3 基于L₁范数约束的方位各向异性贝叶斯AVO反演

线性反射系数方程(式(1))包含两个重要的角度信息,方位角与入射角,通过观察不难发现,方位角信息仅存在于各向异性项中。因此,当同一入射角下、不同方位角反射系数作差,反射系数的差值将仅包括各向异性项。基于上述结论,本研究将基于方位地震数据的各向异性反演分为两个步骤:第一步,利用不同方位角度作差后的地震数据反演裂缝密度e一个独立参数;第二步,将反演得到的裂缝密度作为已知,利用角度地震数据反演各向同性背景流体因子f等3个独立参数。该两步反演策略能够同时利用方位地震数据丰富的入射角与方位角信息,且降低每一步反演过程中的独立参数,提高多参数反演的稳定性。此外,AVO反演相比于弹性阻抗反演规避了由地层弹性信息转换到储层参数信息之间的累计误差的产生,反演结果的准确性更好。

贝叶斯框架能够联合待反演参数先验分布p(m)作为约束,以地震观测数据(实际地震记录)与待反演参数正演合成地震记录的残差p(d|m)作为似然函数,来表征待反演参数的联合后验概率分布p(m|d):

(2)

p (m | d) \propto p (m) p (d | m) = \frac{p (m) p (d | m)}{\int p (m) p (d | m) d m}

。

假设似然函数满足高斯线性模型,先验分布满足高斯分布,则式(2)中的后验概率分布表达如下:

(3)

\begin{array}{c} p (m | d) \propto N (d; G m, \sum_{d}) \cdot N (m; μ_{m}, \sum_{m}) = \\ N (d; G m, \sum_{d} + m; μ_{m}, \sum_{m}) \end{array}

式中:N表示正态分布;G表示正演算子;d与m分别表示观测数据(地震数据)与待反演模型参数;μ_m代表先验模型的均值; $\sum$ _d和 $\sum$ _m分别代表观测数据与先验模型的协方差。在第一步反演过程中,考虑多个方位角度作差,贝叶斯框架下的有效最大后验概率分布p_step1建立如下:

(4)$\begin{array}{c}p_{\text {step1 }}=\sum_{k=1}^{K} \sum_{l=1}^{L}\left[\Delta \boldsymbol{d}\left(\theta_{k}, \varphi_{l}\right)-\boldsymbol{W}_{k, l} \cdot \Delta \boldsymbol{R}_{\mathrm{pp}}^{\mathrm{HTI}}\left(\theta_{k}, \varphi_{l},\right.\right. \\\left.\left.\boldsymbol{m}_{\text {ani }}\right)\right]^{\mathrm{T}} \sum_{\Delta d}^{-1}\left[\Delta \boldsymbol{d}\left(\theta_{k}, \varphi_{l}\right)-\boldsymbol{W}_{k, l} \cdot \Delta \boldsymbol{R}_{\mathrm{pp}}^{\mathrm{HTI}}\right. \\\left.\left(\theta_{k}, \varphi_{l}, \boldsymbol{m}_{\text {ani }}\right)\right]+\left[\boldsymbol{m}_{\text {ani }}-\boldsymbol{m}_{\text {ani }, 0}\right]^{\mathrm{T}} \sum_{\boldsymbol{m}_{\text {ani }}}^{-1} \\{\left[\boldsymbol{m}_{\text {ani }}-\boldsymbol{m}_{\text {ani }, 0}\right]}\end{array}$,

式中:k和l分别代表第k个入射角和第l个方位角度的差值;K和L分别代表入射角与方位角的角度度差值的总个数;Δd(θ_k,φ_l)、W_k_,_l和ΔR_pp(θ_k,φ_l,m_ani)分别表示作差后的观测地震数据,地震子波矩阵和以流体因子、剪切模量、裂缝密度为待反演参数的作差后的PP波反射系数;m_ani=[ $\begin{array}{c} Δ e_{1} & Δ e_{2} & \dots & Δ e_{N} \end{array}$ ]^T和m_ani,0[ $\begin{array}{c} Δ e_{1,0} & Δ e_{2,0} & \dots & Δ e_{N, 0} \end{array}$ ]^T分别代表待反演的裂缝密度e及先验模型; $\sum_{Δ d}^{- 1}$ 和 $\sum_{m_{a n i}}^{- 1}$ 分别代表作差后地震数据与各向异性待反演参数的协方差。由于贝叶斯框架的最大后验概率解基于L₂范数,强调了反演结果的平滑特征,不利于反演结果分辨率的提升以及流体分布的边界刻画,因此我们采用零阶熵正则化稳定算子^[25]S_ZE(m)作为稀疏性L₁范数约束来表征流体边界、提高反演的分辨率和稳定性。零阶熵正则化约束下的第一步反演目标函数表达如下:

(5)$\begin{array}{c}F_{\text {step1 }}=\sum_{k=1}^{K} \sum_{l=1}^{L}\left[\Delta \boldsymbol{d}\left(\theta_{k}, \varphi_{l}\right)-\boldsymbol{W}_{k, l} \cdot \Delta \boldsymbol{R}_{\mathrm{pp}}^{\mathrm{HTI}}\left(\theta_{k}, \varphi_{l},\right.\right. \\\left.\left.\boldsymbol{m}_{\text {ani }}\right)\right]^{\mathrm{T}} \sum_{\Delta \mathbf{d}}^{-1}\left[\Delta \boldsymbol{d}\left(\theta_{k}, \varphi_{l}\right)-\boldsymbol{W}_{k, l} \cdot \Delta \boldsymbol{R}_{\mathrm{pp}}^{\mathrm{HTI}}\left(\theta_{k}, \varphi_{l},\right.\right. \\\left.\left.\boldsymbol{m}_{\text {ani }}\right)\right]+\left[\boldsymbol{m}_{\text {ani }}-\boldsymbol{m}_{\text {ani }, 0}\right]^{\mathrm{T}} \sum_{\mathbf{m}_{\text {ani }}}^{-1}\left[\boldsymbol{m}_{\text {ani }}-\boldsymbol{m}_{\text {ani }, 0}\right]+ \\S_{\mathrm{ZE}}\left(\boldsymbol{m}_{\text {ani }}\right),\end{array}$

(6)$\begin{array}S_{\mathrm{ZE}}\left(\boldsymbol{m}_{\mathrm{ani}}\right)=\sum_{i=1}^{K+L} \sum_{n=1}^{N}-\frac{\left|m_{\mathrm{ani}, K+L, n}\right|+\boldsymbol{\iota}}{\sum_{n=1}^{N}\left(\left|m_{\mathrm{ani}, K+L, n}\right|+\boldsymbol{\iota}\right)}\\\ln \left[\frac{\left|m_{\mathrm{ani}, K+L, n}\right|+\boldsymbol{\iota}}{\sum_{n=1}^{N}\left(\left|m_{\mathrm{ani}, K+L, n}\right|+\boldsymbol{\iota}\right)}\right],\end{array}$

式中:ι表示一个小值,作为阻尼因子防止分母为零。n和N分别代表第n个待反演参数和待反演参数的总采样点个数。

与第一步反演相同,基于贝叶斯框架,第二步反演各向同性背景流体因子f、剪切模量μ_b以及密度ρ_b的目标函数建立如下:

(7)$\begin{array}{c}F_{\text {step2 }}=\sum_{k=1}^{K} \sum_{q=1}^{Q}\left[\boldsymbol{d}\left(\theta_{k}, \varphi_{q}\right)-\boldsymbol{W}_{k, q} \cdot \boldsymbol{R}_{\mathrm{pp}}^{\mathrm{HTI}}\left(\theta_{k}, \varphi_{q},\right.\right. \\\left.\left.\boldsymbol{m}_{\text {iso }}\right)\right]^{\mathrm{T}} \sum_{\mathbf{d}}^{-1}\left[\mathbf{d}\left(\theta_{k}, \varphi_{q}\right)-\boldsymbol{W}_{k, q} \cdot \boldsymbol{R}_{\mathrm{pp}}^{\mathrm{HTI}}\left(\theta_{k}, \varphi_{q},\right.\right. \\\left.\left.\boldsymbol{m}_{\text {iso }}\right)\right]+\left[\boldsymbol{m}_{\text {iso }}-\boldsymbol{m}_{\text {iso }, 0}\right]^{\mathrm{T}} \sum_{\mathbf{m}_{\text {iso }}}^{-1}\left[\boldsymbol{m}_{\text {iso }}-\boldsymbol{m}_{\text {iso }, 0}\right]+ \\S_{Z E}\left(\boldsymbol{m}_{\text {iso }}\right),\end{array}$

(8)$\begin{aligned}S_{Z E}\left(\boldsymbol{m}_{\text {iso }}\right) & =\sum_{i=1}^{K+Q} \sum_{n=1}^{N}-\frac{\left|m_{\text {ani }, K+Q, n}\right|+\boldsymbol{\iota}}{\sum_{n=1}^{N}\left(\left|m_{\text {ani }, K+Q, n}\right|+\boldsymbol{\iota}\right)} \\& \ln \left[\frac{\left|m_{\text {ani }, K+Q, n}\right|+\boldsymbol{\iota}}{\sum_{n=1}^{N}\left(\left|m_{\text {ani }, K+Q, n}\right|+\boldsymbol{\iota}\right)}\right]\end{aligned},$

式中:$\begin{array}\boldsymbol{m}_{\text {iso }}=\left[\begin{array}{lllllll}\frac{\Delta f_{1}}{f_{1}} & \frac{\Delta f_{2}}{f_{2}} & \cdots & \frac{\Delta f_{N}}{f_{N}} & \frac{\Delta \mu_{b 1}}{\mu_{b 1}} & \frac{\Delta \mu_{b 2}}{\mu_{b 2}} & \cdots \\\frac{\Delta \rho_{b 1}}{\rho_{b 1}} & \frac{\Delta \rho_{b 2}}{\rho_{b 2}} & \cdots & \frac{\Delta \rho_{b N}}{\rho_{b N}}\end{array}\right]^{\mathrm{T}}\end{array}$和m_iso,0= $[\begin{array}{l} \frac{Δ f_{1,0}}{f_{1,0}} & \frac{Δ f_{2,0}}{f_{2,0}} & \dots \end{array} {\begin{array}{l} \frac{Δ f_{N, 0}}{f_{N, 0}} & \frac{Δ μ_{b 1,0}}{μ_{b 1,0}} & \frac{Δ μ_{b 2,0}}{μ_{b 2,0}} & \dots & \frac{Δ μ_{b N, 0}}{μ_{b N, 0}} & \frac{Δ ρ_{b 1,0}}{ρ_{b 1,0}} & \frac{Δ ρ_{b 2,0}}{ρ_{b 2,0}} & \dots & \frac{Δ ρ_{b N, 0}}{ρ_{b N, 0}} \end{array}]}^{T}$ 分别代表各向同性背景待反演参数及先验模型。d(θ_k,φ_l)、W_k_,_l和R_pp(θ_k,φ_l,m_ani)分别表示观测地震数据,地震子波矩阵和以裂缝密度作为待反演参数的PP波反射系数。 $\sum_{d}^{- 1}$ 和 $\sum_{m_{i s o}}^{- 1}$ 分别代表地震数据与各向同性背景待反演参数的协方差。q和Q分别代表第q个方位角与方位角的总数。

由于L₁范数的存在,常规的梯度类算法难以对目标函数进行求解,因此,采用McMC随机优化算法^[26-27]反演岩石物性参数。McMC首先初始化模型,指定马尔可夫链的数量,设置最大迭代次数,并结合先前的似然信息以及初始模型与观测数据之间的概率关系。在每次迭代中,根据提案分布生成候选模型,并根据当前模型评估其接受概率。该算法通过运行多个马尔可夫链,探索参数空间,能够避免求解过程陷入局部极小,并对采样模型进行统计分析,即可得到最终的反演结果。

为测试本研究所提出的理论和算法的有效性,基于表1与表2所示的双层模型,开展基于L₁范数约束的方位各向异性贝叶斯AVO反演测试。图4与图5分别展示了使用表1与表2所示的双层裂缝型含气砂岩与泥岩模型以及双层裂缝型含水砂岩与泥岩模型的反演模型测试结果。可以看出,采用L₁范数约

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 双层裂缝型含气砂岩与泥岩模型反演测试

Fig.4 Inversion test of double-layer fractured gas-sand and mudstone models

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 双层裂缝型含水砂岩与泥岩模型反演测试

Fig.5 Inversion test of double-layer fractured water-sand and mudstone models

束的方位各向异性贝叶斯AVO反演的无噪声反演结果能够与模型值具有较高的吻合度,当合成地震记录噪声时,反演结果的准确性有所降低,但仍然能够与模型具有较好的吻合程度,能够对裂缝型储层的地层变化进行较为准确地描述,为反射系数方程的可靠性与反演算法的有效性提供了有力依据。

2 合成地震记录测试

我们选取了实际裂缝型储层工区某井段的实际测井曲线进行合成地震记录测试,以测试不同反演方法的差异及本研究所提出方法的合理性和在实际工区应用的可行性。

图6展示了用于合成地震记录测试的岩石物理参数模型。在图6中,每一幅子图的横坐标代表模型参数的值域范围,纵坐标代表时间,采样间隔为2 ms采样。以图6所示的岩石物理模型作为输入参数,以含流体因子的HTI反射系数方程(式(A-20))作为正演算子进行正演模拟,可以得到如图7所示的叠前噪声合成方位地震记录。在此基础上,对图7所示的合成地震记录添加信噪比为10的随机噪声,可得如图8所示的信噪比为10的合成地震记录。图7~8分别展示了6个方位角度下的合成地震记录,每个方位角度下的合成地震记录的入射角取值范围为0°~30°,在每一幅子图中,其横坐标代表入射角范围,纵坐标与图6相同,代表采样间隔为2 ms采样的时间范围。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 合成记录测试模型

Fig.6 Models for synthetic record test

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 无噪声合成地震记录

Fig.7 Synthetic seismogram with noise free

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 信噪比为10的合成地震记录

Fig.8 Synthetic seismogram with signal to noise ratio 10

以图7所示的无噪声合成地震记录作为观测数据输入,分别采用贝叶斯两步方位地震反演方法与常规同步方位地震反演方法(常规同步反演方法的目标函数建立及反演优化求解详见附录B),无噪声的反演结果分别如图9中的蓝色和绿色曲线所示。从无噪声的反演结果可以看出,无论是新方法还是常规方法,反演结果都能够与模型值具有较高的吻合程度,能够准确地描述流体因子与裂缝密度的层间变化,验证了所提出的推导方程的精度。由于采用了贝叶斯框架与零阶熵正则化稳定算子进一步融合了先验信息以及待反演参数之间的统计学关系,且新方法的两步反演策略降低了反演过程中独立参数的数量,因此,对比两种方法的反演结果能够看出,新方法相比于常规方法的反演结果与模型值的吻合度更高。

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9 无噪声合成地震记录方位各向异性反演结果

Fig.9 Azimuth anisotropy inversion results from synthetic seismic record with noise free

为进一步测试反演方法的鲁棒性与可靠性,我们采用图8所示的信噪比为10的合成地震记录作为观测数据输入,分别采用贝叶斯两步方位地震反演方法与常规同步方位地震反演方法,得到了分别如图10中的蓝色和绿色曲线所示的信噪比为10的反演结果。可以看出,两种反演方法在含噪声地震数据作为观测数据时的反演结果仍然具有较高的精度,反演结果的精度相比于无噪声的反演结果有所下降,但仍然能够与真实模型具有一定的吻合程度,进一步证明了所推导的反射系数方程的准确性。对比两种反演方法的反演结果,得益于零阶熵正则化约束的贝叶斯框架下的两步反演策略在提升反演结果稳定性与鲁棒性方面的优势,新方法能够更为有效地描述地层流体因子与裂缝密度的层间变化情况,其反演结果的抗噪性更强。为进一步量化反演结果与真实模型之间的差异,图11与图12分别计算了合成地震记录无噪声与信噪比为10的方位各向异性反演结果与模型的误差绝对值。可以看出,无噪声的反演结果误差绝对值较小,随着观测数据噪声的加入,误差绝对值增加,但新方法在无噪声和含噪声反演结果的误差绝对值小于常规方法,其在反演结果的准确性与稳定性方面更具优势。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10 信噪比为10合成地震记录方位各向异性反演结果

Fig.10 Azimuth anisotropy inversion results from synthetic seismic record with signal to noise ratio 10

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11 无噪声方位各向异性反演结果与模型误差绝对值

Fig.11 Absolute values of error with noise free azimuth anisotropy inversion results and models

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12 信噪比为10的方位各向异性反演结果与模型误差绝对值

Fig.12 Absolute values of error with signal to noise ratio of 10’s azimuth anisotropy inversion results and models

合成记录测试验证了所提出方程的精度准确性以及贝叶斯两步反演方法应用于理论模型的可靠性,为本研究所提出的反演方法应用于实际工区获取准确、可靠的裂缝密度预测和流体识别结果提供了依据。

3 实际工区应用

我们选取了中国东部某裂缝型储层工区实际资料对不同反演方法进行测试。

图13展示了实际裂缝型储层工区的叠前方位地震剖面,图中选取了两个方位角(45°和135°)与3个入射角(12°、23°和30°),足以对反问题的求解提供充足的观测数据。地震数据为精细处理后的纯波剖面,能够反映真实振幅随入射角和方位角的变化关系。图11中,横坐标代表CDP点,纵坐标代表时间,采样间隔为2 ms采样。以图13所示的地震剖面作为观测数据输入,分别采用新贝叶斯两步方位地震反演方法与常规同步方位地震反演方法,可以输出并得到如图14与图15所示的实际地震资料反演结果。图14与图15中,横坐标与纵坐标与图11相同,分别代表CDP点和时间(2 ms采样)。此外,为验证实际资料反演结果的准确性与可靠性,在图中的井A和井B位置处,我们展示了与反演剖面相对应弹性参数或裂缝参数盲井的变密度显示的结果。测井中的流体因子 $f$ 与裂缝密度e分别由f=(ρ_b $V_{p}^{2}$ )- $γ_{d r y}^{2}$ (ρ_b $V_{s}^{2}$ )和e= $\frac{3 ϕ_{f}}{4 π χ}$ 计算得到,其中ϕ_f和χ分别为裂缝孔隙度和裂缝纵横比,可由测井解释得到。可以看出,图14中新方法的反演结果与测井解释结果具有较高的吻合度,反演结果的地层连续性较好,能够有效描述不同弹性与裂缝发育程度地层的横向展布与纵向变化,图14a与图14d中的黑色椭圆所圈出的地层内,流体因子f反演结果存在低值异常,与此同时,裂缝密度反演结果存在高值异常,反演结果能够可靠、稳定地刻画此处发育高角度裂缝的含气储层,具有较为准确的裂缝密度预测与流体识别结果。对比新方法与常规方法的反演结果,可以看出,常规方法由于未采用贝叶斯框架与零阶熵正则化,且未能在反演过程中降低每一步反演过程中的独立参数的数量,相比于新方法在井A和井B处的反演结果,其准确性有所降低。此外,图14~15中的每幅子图中的黄色箭头所指示的位置处,常规方法的反演结果对含气裂缝型储层的展布刻画效果欠佳,特别是地层的连续性方面与流体因子与裂缝参数的异常值响应方面,其效果相比于新方法较差。实际地震资料的二维剖面测试验证了所推导的反射系数方程的精度以及本研究所提出的反演方法的准确性与可靠性,为新方法应用于大范围工区奠定了基础。

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13 实际裂缝型储层工区叠前方位地震剖面

Fig.13 Prestack azimuth seismic profiles in an actual fractured working area

图14

新窗口打开| 下载原图ZIP| 生成PPT

图14 基于新贝叶斯两步方位地震反演方法的反演结果

Fig.14 Inversion results based on the new Bayesian two step azimuth seismic inversion method

图15

新窗口打开| 下载原图ZIP| 生成PPT

图15 基于常规同步方位地震反演方法的反演结果

Fig.15 Inversion results based on the conventional simultaneous azimuth seismic inversion method

在二维剖面测试的基础上,图16展示了将本研究所提出的方位各向异性反演方法用于三维工区的沿层切片结果。可以看出,低流体因子f与高裂缝密度e异常的区域能够有效地刻画发育高角度裂缝的含气砂岩的横向展布,且由于零阶熵正则化稳定算子作为L₁范数约束强调了反演结果中地层边界的稀疏性,裂缝型储层的横向展布具有较为清晰的边界。观察图14中的两口具有油气解释结果的盲井(井A和井B)不难发现,反演结果在高含气井(井A)的位置处具有明显的低流体因子与高裂缝密度异常,而在含水井(井B)的位置处,则不出现低流体因子异常的假象,反演结果具有较好的可靠性与合理性。值得注意的是,流体因子的反演结果与裂缝密度的反演结果具有一定的相关性,低流体因子分布异常的区域裂缝发育程度往往较高。产生该现象的原因在于沿层切片所选的位置位于含气性较高的主力裂缝型储层层段位置处,储层流体主要以天然气填充为主,较少发育水层;另一方面,裂缝作为天然气赋存的有效介质,裂缝发育较强的位置含气性也较强。此外,沿层切片的高裂缝密度异常区域大于低流体因子异常的区域,说明部分裂缝型储层未能较好赋存天然气,符合实际地质情况,预测结果具有一定的合理性。

图16

新窗口打开| 下载原图ZIP| 生成PPT

图16 基于新贝叶斯两步方位地震反演方法的流体因子f (a)和裂缝密度 $e$ (b)反演结果沿层切片

Fig.16 (a) Interlayer slices of fluid factor f (a) and fracture density e (b) inversion results based on the new Bayesian two step azimuth seismic inversion method

实际资料的反演结果进一步证实了所提出方法在中国东部某实际裂缝型储层工区的适用性,为发育高角度裂缝的含气砂岩的裂缝密度预测及流体识别提供了合理、有效的新方法。

4 结论

本研究基于测井交会优选了流体因子f作为流体识别敏感参数,为有效区分裂缝型储层含气及含水储层提供了理论依据,为流体识别的准确性提供了帮助。以流体敏感参数优选结果为参考的基础上,推导了具有较好的入射角与方位角敏感性的反射系数方程,该方程在保持较高精度的同时能够充分利用方位地震道集中的振幅信息与各种角度信息,为储层流体敏感参数与裂缝参数预测提供扎实的理论基础。

本研究采用贝叶斯框架将测井信息引入目标函数的约束中,采用零阶熵正则化L₁范数约束稳定算子作为提升反演的准确性,发展了一种宽方位两步各向异性贝叶斯AVO反演方法,该方法相比于常规各向异性同步AVO反演,能够提供更为准确、可靠的裂缝密度预测与流体识别结果。

本研究所提出的反演方法在中国东部某裂缝型储层工区资料的应用上表现出了较好应用效果,能够有效描述发育高角度裂缝的含气储层的空间分布状态,含气裂缝型储层分布边界刻画清晰,反演结果与测井解释结果具有较高的吻合度,为储层裂缝密度直接预测及流体识别提供了一种新的技术手段。

附录A

方位各向异性反射系数方程推导过程如下:

当裂缝型地层各向异性较弱且裂缝满足旋转不变的假设时,Schoenberg线性滑动模型^[20]给出了HTI介质刚度矩阵C_HTI与各向同性背景刚度矩阵C_b与裂缝矩阵S_c的定量关系如下:

(A-1)

C_{H T I} = C_{b} - C_{b} S_{c} C_{b}

(A-2)

C_{b} = [\begin{array}{c} M_{b} & λ_{b} & λ_{b} & 0 & 0 & 0 \\ λ_{b} & M_{b} & λ_{b} & 0 & 0 & 0 \\ λ_{b} & λ_{b} & M_{b} & 0 & 0 & 0 \\ 0 & 0 & 0 & μ_{b} & 0 & 0 \\ 0 & 0 & 0 & 0 & μ_{b} & 0 \\ 0 & 0 & 0 & 0 & 0 & μ_{b} \end{array}]

(A-3)

S_{c} \approx [\begin{array}{c} δ_{N} / M_{b} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & δ_{T} / μ_{b} & 0 \\ 0 & 0 & 0 & 0 & 0 & δ_{T} / μ_{b} \end{array}]

式中:M_b、μ_b和λ_b分别代表各向同性背景介质的纵波模量、剪切模量和拉梅参数;δ_N与δ_T分别代表垂直裂缝的法向弱度与切向弱度,随裂缝发育程度的增加而增大。将式(A-2)和式(A-3)代入式(A-1),可得由弹性模量与裂缝参数表示的HTI介质刚度矩阵如下:

(A-4)

C_{H T I} = [\begin{array}{c} M_{b} (1 - δ_{N}) & λ_{b} (1 - δ_{N}) & λ_{b} (1 - δ_{N}) & 0 & 0 & 0 \\ λ_{b} (1 - δ_{N}) & M_{b} (1 - \frac{λ_{b}^{2}}{M_{b}^{2}} δ_{N}) & λ_{b} (1 - \frac{λ_{b}}{M_{b}} δ_{N}) & 0 & 0 & 0 \\ λ_{b} (1 - δ_{N}) & λ_{b} (1 - \frac{λ_{b}}{M_{b}} δ_{N}) & M_{b} (1 - \frac{λ_{b}^{2}}{M_{b}^{2}} δ_{N}) & 0 & 0 & 0 \\ 0 & 0 & 0 & μ_{b} & 0 & 0 \\ 0 & 0 & 0 & 0 & μ_{b} (1 - δ_{T}) & 0 \\ 0 & 0 & 0 & 0 & 0 & μ_{b} (1 - δ_{T}) \end{array}]

。

地震波在地层界面处的能量分配有助于理解地震波的传播机理,为利用地震反射波振幅随入射角与方位角的变化关系来预测地下岩石弹性、物性及裂缝发育程度提供了基础。基于一阶散射波稳相法^[21],可以推导得到HTI介质分界面处纵波入射、纵波反射的反射系数方程 $R_{p p}^{H T I}$ (θ)如下:

(A-5)$\begin{array}{l}R_{\mathrm{pp}}^{\mathrm{HTI}}(\theta)=R_{\mathrm{pp}}^{\mathrm{HTI}-\mathrm{iso}}(\theta)+R_{\mathrm{pp}}^{\mathrm{HTI}-\mathrm{ani}}(\theta)=a_{M}(\theta) \frac{\Delta M_{\mathrm{b}}}{M_{\mathrm{b}}}+ \\a_{\mu}(\theta) \frac{\Delta \mu_{\mathrm{b}}}{\mu_{\mathrm{b}}}+a_{\rho}(\theta) \frac{\Delta \rho_{\mathrm{b}}}{\rho_{\mathrm{b}}}+b_{\mathrm{N}}(\theta) \Delta \delta_{\mathrm{N}}+b_{\mathrm{T}}(\theta) \Delta \delta_{\mathrm{T}},\end{array}$

(A-6)$a_{M}(\theta)=\frac{1}{4 \cos ^{2} \theta}, $

(A-7)$a_{\mu}(\theta)=-\left(1-\frac{\lambda_{\mathrm{b}}}{M_{\mathrm{b}}}\right) \sin ^{2} \theta $,

(A-8)$a_{\rho}(\theta)=\frac{1}{2}-\frac{1}{4 \cos ^{2} \theta},$

(A-9)$\begin{array}{c}b_{\mathrm{N}}(\theta)=-\frac{\sec ^{2} \theta}{4}\left[1-\left(1-\frac{\lambda_{\mathrm{b}}}{M_{\mathrm{b}}}\right)\left(\sin ^{2} \theta \cos ^{2} \varphi+\right.\right. \\\left.\left.\cos ^{2} \theta\right)\right]^{2},\end{array}$

(A-10)$b_{\mathrm{T}}(\theta)=\frac{1}{2}\left(1-\frac{\lambda_{\mathrm{b}}}{M_{\mathrm{b}}}\right) \sin ^{2} \theta \sin ^{2} \varphi\left(1-\tan ^{2} \theta \cos ^{2} \varphi\right) \text {, }$

式中:θ为入射角,φ为方位角(观测测线与裂缝法向的夹角),Δ表示上下层介质弹性参数的差值。

Russell 等^[22]结合Gassmann方程^[23]给出的岩石背景项纵波速度与表征流体项f和岩石骨架项s之间的定量关系如下:

(A-11)$V_{\mathrm{p}}=\sqrt{\frac{s+f}{\rho_{\mathrm{b}}}}=\sqrt{\frac{K_{\text {dry }}+(4 / 3) \mu_{\mathrm{b}}+\alpha^{2} \Upsilon}{\rho_{\mathrm{b}}}},$

(A-12)$\alpha=1-\frac{K_{\mathrm{dry}}}{K_{\mathrm{m}}},$

(A-13)$\Upsilon=1 /\left(\frac{\alpha-\phi}{K_{\mathrm{m}}}+\frac{\phi}{K_{\mathrm{f}}}\right),$,

(A-14)$f=\left(\rho_{\mathrm{b}} V_{\mathrm{p}}^{2}\right)-\gamma_{\mathrm{dry}}^{2}\left(\rho_{\mathrm{b}} V_{\mathrm{s}}^{2}\right),$

(A-15)$s=\gamma_{\mathrm{dry}}^{2} \mu_{\mathrm{b}},$

式中:ϕ代表孔隙度,K_dry、K_m和K_f分别代表各向同性背景中的干岩石骨架、岩石基质和孔隙流体的体积模量。 $γ_{d r y}^{2}$ 为各向同性背景干岩石骨架的纵横波速度比,一般取2.233。基于式(A-11)~(A-15)所述的岩石物理关系,采用模型参数化方法进行参数之间的微分近似,反射系数方程的各向同性部分可以被改写为:

(A-16)$\begin{array}{r}R_{\mathrm{pp}}^{\mathrm{HTI}-\mathrm{iso}}(\theta)=\left[\left(1-\frac{\gamma_{\mathrm{dry}}^{2}}{\gamma_{\mathrm{sat}}^{2}}\right) \frac{\sec ^{2} \theta}{4}\right] \frac{\Delta f}{f}+\left[\frac{\gamma_{\mathrm{dry}}^{2}}{4 \gamma_{\mathrm{sat}}^{2}} \sec ^{2} \theta-\right. \\\left.\quad \frac{2}{\gamma_{\mathrm{sat}}^{2}} \sin ^{2} \theta\right] \times \frac{\Delta \mu_{\mathrm{b}}}{\mu_{\mathrm{b}}}+\left[\frac{1}{2}-\frac{\sec ^{2} \theta}{4}\right] \frac{\Delta \rho_{\mathrm{b}}}{\rho_{\mathrm{b}}},\end{array}$

式中: $γ_{s a t}^{2}$ 为各向同性背景饱和流体填充岩石的纵横波速度比。f作为与孔隙流体相关项,其可作为流体识别因子进行储层的流体识别。裂缝参数作为间接反映裂缝发育程度的参数,利用Bakulin 等^[24]给出的干燥或含气裂缝填充模型可以得到裂缝密度e与裂缝参数之间的定量关系如下:

(A-17)$\delta_{\mathrm{N}}=\frac{4}{3\left(1-1 / \gamma_{\mathrm{dry}}^{2}\right) / \gamma_{\mathrm{dry}}^{2}} e,$

(A-18)$\delta_{\mathrm{T}}=\frac{16}{3\left(3-2 / \gamma_{\mathrm{dry}}^{2}\right)} e,$

将式(A-17)、(A-18)代入反射系数方程的各向异性部分,同样采用模型参数化方法进行裂缝参数与裂缝密度之间的微分近似,可得:

(A-19)$R_{\mathrm{pp}}^{\mathrm{HTI}-\mathrm{ani}}(\theta)=\left\{\begin{array}{l}-\frac{\sec ^{2} \theta}{4}\left[1-\left(1-\frac{\lambda_{\mathrm{b}}}{M_{\mathrm{b}}}\right)\left(\sin ^{2} \theta \cos ^{2} \varphi+\cos ^{2} \theta\right)\right]^{2} \frac{4}{3\left(1-1 / \gamma_{\mathrm{dry}}^{2}\right) / \gamma_{\mathrm{dry}}^{2}} \\+\frac{1}{2}\left(1-\frac{\lambda_{\mathrm{b}}}{M_{\mathrm{b}}}\right) \sin ^{2} \theta \sin ^{2} \varphi\left(1-\tan ^{2} \theta \cos ^{2} \varphi\right) \frac{16}{3\left(3-2 / \gamma_{\mathrm{dry}}^{2}\right)} e\end{array}\right\} \Delta e,$

联立式(A-16)与式(A-19),可得含裂缝密度与流体识别因子的反射系数方程如下:

(A-20)$\begin{array}{l}R_{\mathrm{pp}}^{\mathrm{HTI}}(\theta)=\left[\left(1-\frac{\gamma_{\mathrm{dry}}^{2}}{\gamma_{\mathrm{sat}}^{2}}\right) \frac{\sec ^{2} \theta}{4}\right] \frac{\Delta f}{f}+\left[\frac{\gamma_{\mathrm{dry}}^{2}}{4 \gamma_{\mathrm{sat}}^{2}} \sec ^{2} \theta-\frac{2}{\gamma_{\mathrm{sat}}^{2}} \sin ^{2} \theta\right] \times \frac{\Delta \mu_{\mathrm{b}}}{\mu_{\mathrm{b}}}+\left[\frac{1}{2}-\frac{\sec ^{2} \theta}{4}\right] \frac{\Delta \rho_{\mathrm{b}}}{\rho_{\mathrm{b}}} \\+\left\{\begin{array}{l}-\frac{\sec ^{2} \theta}{4}\left[1-\left(1-\frac{\lambda_{\mathrm{b}}}{M_{\mathrm{b}}}\right)\left(\sin ^{2} \theta \cos ^{2} \varphi+\cos ^{2} \theta\right)\right]^{2} \frac{4}{3\left(1-1 / \gamma_{\mathrm{dry}}^{2}\right) / \gamma_{\mathrm{dry}}^{2}} \\+\frac{1}{2}\left(1-\frac{\lambda_{\mathrm{b}}}{M_{\mathrm{b}}}\right) \sin ^{2} \theta \sin ^{2} \varphi\left(1-\tan ^{2} \theta \cos ^{2} \varphi\right) \frac{16}{3\left(3-2 / \gamma_{\mathrm{dry}}^{2}\right)} e\end{array}\right\} \Delta e\end{array}$

附录B

常规各向异性AVO同步反演方法将利用振幅同时随入射角和方位角的变化直接反演包括各向同性背景项中的流体因子f、剪切模量μ_b、密度ρ_b和各向异性项中的裂缝密度e等4个待反演的裂缝型储层参数。基于贝叶斯框架,其反演的目标函数p_sim建立如下:

(B-1)$\begin{array}{c}p_{\text {sim }}=\sum_{k=1}^{K} \sum_{u=1}^{U}\left[\boldsymbol{d}\left(\theta_{k}, \varphi_{u}\right)-\boldsymbol{W}_{k, u} \cdot \boldsymbol{R}_{\mathrm{pp}}^{\mathrm{HT}}\left(\theta_{k}, \varphi_{u}, \boldsymbol{m}_{\text {all }}\right)\right]^{\mathrm{T}} \\\sum_{d}^{-1}\left[\boldsymbol{d}\left(\theta_{k}, \varphi_{u}\right)-\boldsymbol{W}_{k, u} \cdot \boldsymbol{R}_{\mathrm{pp}}^{\mathrm{HTI}}\left(\theta_{k}, \varphi_{u}, \boldsymbol{m}_{\text {all }}\right)\right]+ \\{\left[\boldsymbol{m}_{\text {all }}-\boldsymbol{m}_{\text {all }, 0}\right]^{\mathrm{T}} \sum_{m_{\text {all }}}^{-1}\left[\boldsymbol{m}_{\text {all }}-\boldsymbol{m}_{\text {all }, 0}\right],}\end{array}$

式中:分别表示观测地震数据,地震子波矩阵和以裂缝密度作为待反演参数的PP波反射系数。和分别代表地震数据与4个待反演参数(流体因子f、剪切模量μ_b、密度ρ_b以及裂缝密度e)的协方差。k和K分别代表第k个入射角与入射角的总数,u和U分别代表第u个方位角与方位角的总数。m_all与m_all,0分别代表4个待反演参数及其初始模型,表达式如下:

(B-2)$\begin{array} \boldsymbol{m}_{\mathrm{all}}=\left[\begin{array}{lllllllllllllll}\frac{\Delta f_{1}}{f_{1}} & \frac{\Delta f_{2}}{f_{2}} & \cdots & \frac{\Delta f_{U}}{f_{U}} & \frac{\Delta \mu_{\mathrm{b} 1}}{\mu_{\mathrm{b} 1}} & \frac{\Delta \mu_{\mathrm{b} 2}}{\mu_{\mathrm{b} 2}} & \cdots & \frac{\Delta \mu_{\mathrm{b} U}}{\mu_{\mathrm{b} U}} & \frac{\Delta \rho_{\mathrm{b} 1}}{\rho_{\mathrm{b} 1}} & \frac{\Delta \rho_{\mathrm{b} 2}}{\rho_{\mathrm{b} 2}} & \cdots & \frac{\Delta \rho_{\mathrm{b} U}}{\rho_{\mathrm{b} U}} & \frac{\Delta e_{1}}{e_{1}} & \\ \frac{\Delta e_{2}}{e_{2}} & \cdots\end{array} \frac{\Delta e_{U}}{e_{U}}\right]^{\mathrm{T}}\end{array},$

(B-3)$\begin{array}\boldsymbol{m}_{\mathrm{all}, 0}=\left[\begin{array}{lllllllllll}\frac{\Delta f_{1,0}}{f_{1,0}} & \frac{\Delta f_{2,0}}{f_{2,0}} & \cdots & \frac{\Delta f_{U, 0}}{f_{U, 0}} & \frac{\Delta \mu_{\mathrm{b} 1,0}}{\mu_{\mathrm{b} 1,0}}& \frac{\Delta \mu_{\mathrm{b} 2,0}}{\mu_{\mathrm{b} 2,0}} & \cdots & \frac{\Delta \mu_{\mathrm{b} U, 0}}{\mu_{\mathrm{b} U, 0}} & \frac{\Delta \rho_{\mathrm{b} 1,0}}{\rho_{\mathrm{b} 1,0}} & \frac{\Delta \rho_{\mathrm{b} 2,0}}{\rho_{\mathrm{b} 2,0}} & \cdots \end{array} \frac{\Delta \rho_{\mathrm{b} U, 0}}{\rho_{\mathrm{b} U, 0}} \frac{\Delta e_{1,0}}{e_{1,0}} \frac{\Delta e_{2,0}}{e_{2,0}}\right.\\\left.\ldots \frac{\Delta e_{U, 0}}{e_{U, 0}}\right]^{\mathrm{T}}\end{array}$

目标函数的求解同样采用McMC随机优化算法^[26-27]。首先初始化模型,并指定马尔可夫链的数量,设置最大迭代次数。在每次迭代中,根据提案分布生成候选模型,并根据当前模型评估其接受概率。该算法通过运行多个马尔可夫链,探索参数空间,并对采样模型进行统计分析,即可得到最终的4个参数的同时反演结果。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[19]

严彬, 张广智, 李林, 等.

裂缝诱导的TTI介质固液解耦反射系数方程及裂缝和流体参数反演方法

[J]. 地球物理学报, 2023, 66(10):4349-4369.