基于叠前直接反演的页岩气储层脆性及裂缝参数预测方法

doi:10.11720/wtyht.2025.1506

基于叠前直接反演的页岩气储层脆性及裂缝参数预测方法

石学文¹^,², 王畅¹^,², 张洞君¹^,², 冯艳雯¹^,²

1.中国石油西南油气田分公司页岩气研究院, 四川成都 610051

2.页岩气开发与评价四川省重点实验室, 四川成都 610051

A method for predicting the brittleness and fracture parameters of shale gas reservoirs based on prestack direct inversion

SHI Xue-Wen¹^,², WANG Chang¹^,², ZHANG Dong-Jun¹^,², FENG Yan-Wen¹^,²

1. Shale Gas Research Institute, Southwest Oil & Gasfield Company,PetroChina, Chengdu 610051, China

2. Sichuan Key Laboratory of Shale Gas Evaluation and Exploitation, Chengdu 610051, China

第一作者: 石学文(1982-),男,高级工程师,2005年毕业于西南石油大学资源勘查专业,2017年毕业于西南石油大学地质工程专业,获硕士学位,现从事页岩气勘探开发工作。

责任编辑: 叶佩

收稿日期: 2025-01-22 修回日期: 2025-05-9

基金资助:

中国石油天然气股份有限公司科技项目(2023ZZ21YJ04)

Received: 2025-01-22 Revised: 2025-05-9

摘要

目前常规各向异性介质的反演容易受小角度入射、地层弱性质变化的假设条件限制,并且岩石弹性参数的预测多通过线性反演弹性参数进而间接拟合得到,使得各向异性地层的岩石物理参数反演结果的精度与可靠性较低。为此,本文推导了具有垂直对称轴的横向各向同性(transverse isotropy with vertical axis of symmetry,VTI)介质杨氏模量、泊松比以及裂缝参数的岩石物理模型,基于精确VTI反射系数方程,提出了一种L₁范数约束的贝叶斯各向异性非线性直接反演方法,用于对页岩气储层的脆性及裂缝参数进行直接反演预测。该方法在西南页岩气工区取得了良好的应用效果,为页岩气的储层表征提供了新的方法。

关键词： VTI; 非线性; 直接反演; 脆性参数; 裂缝参数

Abstract

The conventional inversion of anisotropic media is often constrained by assumptions of narrow-angle incidence and weak changes in stratigraphic properties.Moreover,the prediction of rock elastic parameters typically involves linear inversion and indirect fitting,leading to less accurate and reliable inversion results of petrophysical parameters for anisotropic formations.Hence,this study derived the petrophysical models of Young's modulus,Poisson's ratio,and fracture parameters for vertical transverse isotropy(VTI) media.Based on the precise VTI reflection coefficient equation,this study proposed a Bayesian anisotropic nonlinear direct inversion method,constrained by the L₁ norm,to predict the brittleness and fracture parameters of shale gas reservoirs through direct inversion.The proposed method yielded satisfactory application results in a study area of shale gas in Southwest China,offering a novel technique for characterizing shale gas reservoirs.

Keywords： vertical transverse isotropy(VTI); nonlinear; direct inversion; brittleness parameter; fracture parameter

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

本文引用格式

石学文, 王畅, 张洞君, 冯艳雯. 基于叠前直接反演的页岩气储层脆性及裂缝参数预测方法[J]. 物探与化探, 2025, 49(4): 826-837 doi:10.11720/wtyht.2025.1506

SHI Xue-Wen, WANG Chang, ZHANG Dong-Jun, FENG Yan-Wen. A method for predicting the brittleness and fracture parameters of shale gas reservoirs based on prestack direct inversion[J]. Geophysical and Geochemical Exploration, 2025, 49(4): 826-837 doi:10.11720/wtyht.2025.1506

0 引言

中国页岩储层分布范围广,含气储量丰富,随着近年来的非常规油气勘探开发技术的不断成熟和进步,页岩气的开发与利用逐渐成为了当今能源与科学研究的热点^[1-2]。

含油气储层预测的地球物理方法多采用叠前地震反演岩性或流体识别敏感参数的手段,通过挖掘地震振幅中蕴含的丰富信息,为储层表征及流体识别提供可靠的预测结果。Zong等^[3]推导了关于杨氏模量、泊松比的弹性阻抗方程,并提出了一种利用叠前地震数据获得杨氏模量和泊松比的弹性阻抗反演方法。该方法避免了杨氏模量和泊松比间接计算带来的累计误差,为页岩储层脆性及可压裂性的研究提供了新的思路。Grana^[4]开发了贝叶斯线性化岩石物理反演方法,实现了对地震属性中的岩石和流体性质的估算,反演结果在保持较高计算效率的同时,提高了参数不确定性评估能力,为利用岩石物理模型从地震属性反演岩石和流体性质提供了一种新的有效途径。de Figueiredo等^[5]在叠前岩石物性参数反演研究的基础上,将地质信息和物理约束进一步融入反演过程,提出了一种用于从地震数据中估算弹性和岩石物理性质的联合贝叶斯反演方法。此方法推导了空间独立条件分布的解析表达式,并将其与地统计模拟方法相结合,计算空间相关的后验分布,提高了储层性质估算的准确性,为油气勘探开发提供了更有价值的信息。考虑地层黏弹性质导致的地震波传播的衰减,Chen等^[6]基于逆品质因子表达近似反射系数和衰减弹性阻建立了一种两步弹性阻抗反演方法,提出了从不同入射角和频率的观测地震数据中估算岩石弹性性质以及纵横波品质逆因子的叠前反演方法,为从地震数据中估算弹性参数和衰减特性提供了一种有效的手段。

众多研究结果表明,页岩储层通常发育有一组或者多组的水平或垂直裂缝,存在一定的各向异性特征^[7-10]。将处于简单构造带上的页岩储层进行一定的近似与简化,可将其看作具有垂直对称轴的横向各向同性(transverse isotropy with vertical axis of symmetry,VTI)介质,经过叠前各向异性反演表征裂缝发育状态的各向异性参数,可对地下页岩储层的含气性进行评估。Plessix等^[11]在系统研究垂直横向各向异性对反射系数影响的基础上,利用纵横波联合的多波多分量叠前地震数据开展了VTI介质的多参数反演方法,对提高VTI介质参数的反演精度从而更好地利用叠前反演技术进行油气勘探和岩性识别提供了重要的理论基础与实践意义。在多波探测研究发展的浪潮中,侯栋甲等^[12]基于贝叶斯理论提出了一种适用于VTI介质中的多波数据的叠前纵横波联合反演方法,旨在采用Ruger近似来表示各向异性介质中的反射系数作为正演方程,并使用最小二乘优化联合反演密度比、速度比以及各向异性参数差等5个岩石物理参数;在反演中其进一步引入了贝叶斯理论,假设未知参数服从高斯先验信息和改进的柯西分布,与仅使用纵波数据的常规反演相比,提高了反演的精度、抗噪能力和稳定性。在上述研究基础上,Lu等^[13]发展了一种各向异性振幅随偏移距变化(amplitude variation with offset,AVO)的反演方法,该方法能够同时反演弹性参数和Thomsen各向异性参数,联合使用了PP和PS地震数据,考虑了各向异性AVO 反演中的局部极小值的问题。解决策略采用分步反演,在两个阶段进行等向性和各向异性AVO反演,并使用第一阶段的结果来约束第二阶段的反演,此外,研究中改进了基于欧几里德距离相似性的停止条件来控制噪声情况下的迭代。该研究在识别储层中细层的复杂岩性和流体信息方面具有很好的可行性与有效性,在处理复杂地质条件下的储层特征识别方面具有重要意义。基于精确VTI方程,Lang等^[14]给出了一种可以克服部分现有线性近似局限性的用于横向各向同性(VTI)介质的PP 波反射系数方程,新反射系数方程与精确各向异性反射系数方程相比,通过相对值替换,减少了反演参数的数量,之后,使用马尔可夫链蒙特卡罗(Markov chain Monte Carlo,McMC)算法的贝叶斯反演框架,得到了5个反射率参数的反演结果。该方法进一步提高了反演的稳定性,为含油气页岩储层各向异性地震解释提供了一种改进方案。

上述研究内容中的反演策略大多基于叠前线性反演算法以及基于相对值反演的非线性反演算法,缺乏针对VTI介质的弹性、物性参数的绝对值直接反演方法的研究,使得反演结果容易受到小角度入射、地层变化较弱等假设条件的限制,难以对地层岩石物理性质的真实变化进行可靠的描述。因此,本文基于精确VTI非线性反射系数方程,采用Schoenberg线性滑动模型进行裂缝型VTI介质的岩石物理建模,从而推导了含杨氏模量、泊松比以及岩石裂缝发育参数的非线性反射系数方程;基于贝叶斯框架下的有效最大后验概率解并结合L₁范数约束建立反演目标函数,采用McMC非线性优化算法进行目标函数的求解,获得目标函数取得极小值时的待反演参数的数值解。该方法在模型测试与实际资料测试的应用中均表现出了合理、稳定的各向异性反演结果,为页岩储层脆性及裂缝发育评价提供了更可靠的手段。

1 理论与方法

裂缝型页岩储层通常发育有水平层状薄裂隙,在页岩储层中,裂缝法向弱度与切向弱度可作为评估裂缝发育程度的裂缝参数,而杨氏模量、泊松比则常作为反映地层脆性的敏感参数。本研究将针对裂缝型页岩储层建立含脆性及裂缝参数的VTI介质刚度矩阵,在此基础上推导VTI介质下的非线性反射系数方程,并开展页岩储层脆性及裂缝参数的叠前非线性直接反演。具体理论与方法阐明如下。

1.1 裂缝型页岩储层脆性及裂缝参数模型构建

当各向同性介质中发育一组定向排列的裂隙时,Bakulin等^[15]认为其有效刚度矩阵可以由背景柔度矩阵S_b与附加裂隙柔度矩阵S_c之和描述:

(1)$ \begin{aligned}\boldsymbol{C}=\left[\boldsymbol{S}_{\mathrm{b}}+\boldsymbol{S}_{\mathrm{c}}\right]^{-1}= & {\left[\boldsymbol{S}_{\mathrm{b}}\left(\boldsymbol{I}+\boldsymbol{S}_{\mathrm{c}} \boldsymbol{S}_{\mathrm{b}}^{-1}\right)\right]^{-1}=\boldsymbol{C}_{\mathrm{b}}\left[\boldsymbol{I}+\boldsymbol{S}_{\mathrm{c}} \boldsymbol{C}_{\mathrm{b}}\right]^{-1}=} \\& \boldsymbol{C}_{\mathrm{b}}\left[\boldsymbol{I}+\sum_{k=0}^{\infty}\left(-\boldsymbol{S}_{\mathrm{c}} \boldsymbol{C}_{\mathrm{b}}\right)^{k}\right],\end{aligned}$

在弱各向异性近似假设下,考虑裂隙密度较小并仅保留式(1)中泰勒展开中的一次项,可得:

(2)$ \boldsymbol{C} \approx \boldsymbol{C}_{\mathrm{b}}-\boldsymbol{C}_{\mathrm{b}} \boldsymbol{S}_{\mathrm{c}} \boldsymbol{C}_{\mathrm{b}}$

我们将采用Schoenberg线性滑动模型^[16-17]对含裂隙储层其进行描述。该模型假定裂隙等效为理想情况下的一个零厚度平面(无限应变),并且满足位移不连续(跨越裂隙面的不连续位移与裂隙面上作用的最小应力线性相关)的特征。对于裂缝型储层,在小裂隙弱度假设下,裂隙的柔度矩阵可线性表示为:

(3)

${S}_{c}\approx \left[\begin{array}{llllll}{\delta }_{N}/{C}_{11b}& 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& {\delta }_{V}/{C}_{44b}& 0\\ 0& 0& 0& 0& 0& {\delta }_{H}/{C}_{66b}\end{array}\right]$

可得出在各向同性背景下发育垂直裂缝的表达式:

(4)$\boldsymbol{C}_{\mathrm{HTI}}=\boldsymbol{C}_{\mathrm{iso}}-\boldsymbol{C}_{\mathrm{iso}} \boldsymbol{S}_{\mathrm{c}} \boldsymbol{C}_{\mathrm{iso}}。$

针对页岩储层的含水平裂缝特征,根据VTI介质与具有水平对称轴的横向各向同性(transverse isotropy with horizontal symmetry axis, HTI)介质刚度矩阵转换关系,结合各向同性背景刚度矩阵,可得出各向同性背景下发育水平裂隙的裂缝型页岩储层刚度矩阵C_VTI表达式:

(5)

${C}_{VTI}=\left[\begin{array}{llllll}{M}_{b}(1-{\chi }_{b}^{2}{\delta }_{N})& {\lambda }_{b}(1-{\chi }_{b}{\delta }_{N})& {\lambda }_{b}(1-{\delta }_{N})& 0& 0& 0\\ {\lambda }_{b}(1-{\chi }_{b}{\delta }_{N})& {M}_{b}(1-{\chi }_{b}^{2}{\delta }_{N})& {\lambda }_{b}(1-{\delta }_{N})& 0& 0& 0\\ {\lambda }_{b}(1-{\delta }_{N})& {\lambda }_{b}(1-{\delta }_{N})& {M}_{b}(1-{\delta }_{N})& 0& 0& 0\\ 0& 0& 0& {\mu }_{b}(1-{\delta }_{T})& 0& 0\\ 0& 0& 0& 0& {\mu }_{b}(1-{\delta }_{T})& 0\\ 0& 0& 0& 0& 0& {\mu }_{b}\end{array}\right]$

式中:χ_b=$\frac{{\lambda }_{b}}{{M}_{b}}$=$\frac{{\sigma }_{b}}{(1-{\sigma }_{b})}$;M_b为背景介质纵波模量;σ_b表示各向同性背景岩石的泊松比;λ_b与μ_b分别为背景介质的拉梅参数;δ_N与δ_T分别表示裂隙的法向与切向弱度的裂缝参数,用于表征由裂隙引起的裂缝型介质的柔度量,可由如下表达式估算:

(6)

${\delta }_{N}=\frac{4e}{3g(1-g)\left[1+\frac{1}{\pi (1-g)}\left(\frac{M\text{'}}{\mu \chi }\right)\right]}$

(7)

${\delta }_{T}=\frac{16e}{3(3-2g)\left[1+\frac{4}{\pi (3-2g)}\left(\frac{\mu \text{'}}{\mu \chi }\right)\right]}$

式(6)~(7)中,g=μ_b/(M_b);e为裂缝密度;χ为裂隙纵横比;M'和μ'分别代表裂缝的纵波模量和剪切模量。对于页岩气储层,δ_N与δ_T可以简化为:

(8)

${\delta }_{N}=\frac{4e}{3g(1-g)}$

(9)

${\delta }_{T}=\frac{16e}{3(3-2g)}$

式(5)~(9)推导了含裂缝参数发育有水平层状薄裂隙的页岩气储层的等效刚度矩阵。在C_VTI的各向同性背景项中,杨氏模量、泊松比将取代拉梅参数与纵波模量,其转换表达式如下:

(10)

${\lambda }_{b}=\frac{{E}_{b}{\sigma }_{b}}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}$

(11)

${\mu }_{b}=\frac{{E}_{b}}{2(1+{\sigma }_{b})}$

(12)

${M}_{b}=\frac{{E}_{b}(1-{\sigma }_{b})}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}$

将式(10)~(12)代入式(5),可得:

(13)

${C}_{VTI}=\left[\begin{array}{ll}{C}_{VTI}^{1}& 0\\ 0& {C}_{VTI}^{2}\end{array}\right]$

(14)

${C}_{VTI}^{1}=\left[\begin{array}{lll}\frac{{E}_{b}(1-{\sigma }_{b})}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}\left(1-\frac{{\sigma }_{b}^{2}}{(1-{\sigma }_{b}{)}^{2}}{\delta }_{N}\right)& \frac{{E}_{b}{\sigma }_{b}}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}\left(1-\frac{{\sigma }_{b}}{(1-{\sigma }_{b})}{\delta }_{N}\right)& \frac{{E}_{b}{\sigma }_{b}}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}\left(1-\frac{{\sigma }_{b}}{(1-{\sigma }_{b})}{\delta }_{N}\right)\\ \frac{{E}_{b}{\sigma }_{b}}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}\left(1-\frac{{\sigma }_{b}}{(1-{\sigma }_{b})}{\delta }_{N}\right)& \frac{{E}_{b}{\sigma }_{b}}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}\left(1-\frac{{\sigma }_{b}}{(1-{\sigma }_{b})}{\delta }_{N}\right)& \frac{{E}_{b}{\sigma }_{b}}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}(1-{\delta }_{N})\\ \frac{{E}_{b}{\sigma }_{b}}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}(1-{\delta }_{N})& \frac{{E}_{b}{\sigma }_{b}}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}(1-{\delta }_{N})& \frac{{E}_{b}{\sigma }_{b}}{(1+{\sigma }_{b})(1-2{\sigma }_{b})}(1-{\delta }_{N})\end{array}\right]$

(15)${C}_{VTI}^{2}=\left[\begin{array}{lll}\frac{{E}_{b}}{2(1+{\sigma }_{b})}(1-{\delta }_{T})& 0& 0\\ 0& \frac{{E}_{b}}{2(1+{\sigma }_{b})}(1-{\delta }_{T})& 0\\ 0& 0& \frac{{E}_{b}}{2(1+{\sigma }_{b})}\end{array}\right]$。

式中:E_b和σ_b分别为各向同性背景中的杨氏模量与泊松比。式(13)~(15)建立了含脆性及裂缝参数的缝型页岩储层刚度矩阵,该刚度矩阵在量化裂缝参数作用的同时,也能够为进一步反映储层的脆性及裂缝发育程度提供基础。

1.2 VTI非线性反射系数方程构建

Zoeppritz等^[18]对各向同性完全弹性假设下的地震波反射系数随入射角的变化给予了精确的描述。与Zoeppritz类似,Rüger^[19]总结了VTI介质下的精确反射系数方程,表述如下:

(16)$M R=b$,

式中:R代表反射系数序列,R=[${R}_{pp}^{VTI}$${R}_{ps}^{VTI}$${T}_{pp}^{VTI}$${T}_{ps}^{VTI}$]^T;b=[$\begin{array}{llll}-{m}_{11}& -{m}_{21}& {m}_{31}& {m}_{41}\end{array}$]^T。M中的各个分量m_ij的精确表达式展开如下:

(17)

${m}_{11}={l}_{\alpha }^{\left(1\right)}$

(18)

${m}_{12}={m}_{\beta }^{\left(1\right)}$

(19)

${m}_{13}=-{l}_{\alpha }^{\left(2\right)}$

(20)

${m}_{14}=-{m}_{\beta }^{\left(2\right)}$

(21)

${m}_{31}={m}_{\alpha }^{\left(1\right)}$

(22)

${m}_{32}=-{l}_{\beta }^{\left(1\right)}$

(23)

${m}_{33}={m}_{\alpha }^{\left(2\right)}$

(24)

${m}_{34}=-{l}_{\beta }^{\left(2\right)}$

(25)

${m}_{21}=p{l}_{\alpha }^{\left(1\right)}{a}_{13}^{\left(1\right)}+{q}_{\alpha }^{\left(1\right)}{m}_{\alpha }^{\left(1\right)}{a}_{33}^{\left(1\right)}$

(26)

${m}_{22}=p{m}_{\beta }^{\left(1\right)}{a}_{13}^{\left(1\right)}-{q}_{\beta }^{\left(1\right)}{l}_{\beta }^{\left(1\right)}{a}_{33}^{\left(1\right)}$

(27)

${m}_{23}=-(p{l}_{\alpha }^{\left(2\right)}{a}_{13}^{\left(2\right)}+{q}_{\alpha }^{\left(2\right)}{m}_{\alpha }^{\left(2\right)}{a}_{33}^{\left(2\right)})$

(28)

${m}_{24}=-(p{m}_{\beta }^{\left(2\right)}{a}_{13}^{\left(2\right)}-{q}_{\beta }^{\left(2\right)}{l}_{\beta }^{\left(2\right)}{a}_{33}^{\left(2\right)})$

(29)

${m}_{41}={a}_{55}^{\left(1\right)}({q}_{\alpha }^{\left(1\right)}{l}_{\alpha }^{\left(1\right)}+p{m}_{\alpha }^{\left(1\right)})$

(30)

${m}_{42}={a}_{55}^{\left(1\right)}({q}_{\beta }^{\left(1\right)}{m}_{\beta }^{\left(1\right)}-p{l}_{\beta }^{\left(1\right)})$

(31)

${m}_{43}={a}_{55}^{\left(2\right)}({q}_{\alpha }^{\left(2\right)}{l}_{\alpha }^{\left(2\right)}+p{m}_{\alpha }^{\left(2\right)})$

(32)

${m}_{44}={a}_{55}^{\left(2\right)}({q}_{\beta }^{\left(2\right)}{m}_{\beta }^{\left(2\right)}-p{l}_{\beta }^{\left(2\right)})$

(33)

${l}_{\alpha }=\sqrt{\frac{{a}_{33}{q}_{\alpha }^{2}+{a}_{55}{p}^{2}-1}{{a}_{11}{p}^{2}+{a}_{55}{q}_{\alpha }^{2}-1+{a}_{33}{q}_{\alpha }^{2}+{a}_{55}{p}^{2}-1}}$

(34)

${m}_{\alpha }=\sqrt{\frac{{a}_{55}{q}_{\alpha }^{2}+{a}_{11}{p}^{2}-1}{{a}_{11}{p}^{2}+{a}_{55}{q}_{\alpha }^{2}-1+{a}_{33}{q}_{\alpha }^{2}+{a}_{55}{p}^{2}-1}}$

(35)

${l}_{\beta }=\sqrt{\frac{{a}_{55}{q}_{\beta }^{2}+{a}_{11}{p}^{2}-1}{{a}_{11}{p}^{2}+{a}_{55}{q}_{\beta }^{2}-1+{a}_{33}{q}_{\beta }^{2}+{a}_{55}{p}^{2}-1}}$

(36)

${m}_{\beta }=\sqrt{\frac{{a}_{33}{q}_{\beta }^{2}+{a}_{55}{p}^{2}-1}{{a}_{11}{p}^{2}+{a}_{55}{q}_{\beta }^{2}-1+{a}_{33}{q}_{\beta }^{2}+{a}_{55}{p}^{2}-1}}$

式中:上角标(1)和(2)分别表示反射界面的上、下地层;p=sinθ/V_p(θ),表示水平慢度;a_ij=C_VTI,_ij/ρ,表示刚度与密度的比值;ρ表示密度;q_α与q_β分别为垂直对称轴方向(VTI介质竖直方向)的纵波慢度和横波慢度。

为降低反演过程的复杂性并提高反演的稳定性,将对精确VTI方程进行适当简化。基于弱各向异性假设,相速度公式可以简化并表达如下:

(37)

${V}_{p}\left(\theta \right)=\sqrt{{a}_{33}}(1+\delta si{n}^{2}\theta co{s}^{2}\theta +\epsilon si{n}^{4}\theta )$

(38)

$\epsilon =\frac{{C}_{VTI,11}-{C}_{VTI,33}}{2{C}_{VTI,33}}$

(39)

$\delta =\frac{({C}_{VTI,13}+{C}_{VTI,55}{)}^{2}-({C}_{VTI,33}-{C}_{VTI,55}{)}^{2}}{2{C}_{VTI,33}({C}_{VTI,33}-{C}_{VTI,55})}$

式中:ε与δ为另一种表征裂隙发育状态的Thomsen各向异性参数,将式(13)~(15)代入式(16)~(39)所示的VTI精确反射系数方程中,可得:

(40)

$\begin{array}{l}{R}_{pp}^{VTI}\left(\theta \right)=f({E}_{b}^{\left(1\right)}　{E}_{b}^{\left(2\right)}　{\sigma }_{b}^{\left(1\right)}　{\sigma }_{b}^{\left(2\right)}　{\rho }^{\left(1\right)}　{\rho }^{\left(2\right)}\\ {\delta }_{N}^{\left(1\right)}　{\delta }_{N}^{\left(2\right)}　{\delta }_{T}^{\left(1\right)}　{\delta }_{T}^{\left(2\right)})。\end{array}$

式(40)即为推导得到的含缝型页岩储层脆性及裂缝参数的非线性精确反射系数方程,其描述了地层反射系数${R}_{pp}^{VTI}$与5个独立岩石物理参数(E_b、σ_b、ρ、δ_N以及δ_T)之间的非线性关系。其中${R}_{pp}^{VTI}$的获取可由克莱姆法则或者矩阵求逆来得到,众多学者都给出过详细的表达式,如Graebner^[20],此处不再赘述。新推导的反射系数方程可适用于中、大角度入射以及地层变化较强的实际情况,但其使用了Schoenberg线性滑动模型来表征水平裂缝发育程度,且结合了近似相速度公式对反射系数方程进行了简化,因此反射系数方程的应用将受到有水平裂缝发育的页岩工区,且受到弱各向异性的假设条件限制。

方程的精度分析将用于验证所推导的反射系数方程的准确性与合理性。以表1与表2所示的两种上下含气页岩介质的岩石物理参数为依据,可以得到反射系数方程的精度分析结果,分别如图1与图2所示。在图1与图2中,红色实线表示精确VTI方程,蓝色加号表示本文中式(40)所推导的简化后的含脆性以及裂缝参数的非线性VTI方程。可以看出,在不同模型中,反射系数方程在小、中入射角度下都能够与精确VTI方程高度吻合,在大角度入射甚至入射角大于30°的情况下,误差逐渐增大,但是仍能够与精确方程具有较高的吻合度。产生该结果的原因在于反射系数方程的推导过程基于精确VTI方程,并且保留了反射系数与岩石物理参数之间的非线性关系。同时,由于方程中所包含的岩石物理参数直接定义为上、下层介质,且并未表示成反射率的形式,因此,其不受地层岩石物理参数弱变化的限制,能够在地层变化较为剧烈的区域为叠前反演给予更真实的方程基础。以上分析结果证明了新方程在精度方面的可靠性,证实了本文方程推导的合理性与准确性,为本文所提出的非线性方程应用于页岩气储层并实现储层脆性参数与裂缝参数的准确预测奠定了基础。

表1 含气页岩双层模型1

Table 1 The double-layer model 1 of gas-bearing shales

	E/GPa	σ	V_p/ (m·s^-1)	V_s/ (m·s^-1)	ρ/ (kg·m^-3)	δ_N	δ_T
上层页岩	21.28	0.272	3268	1829	2500	0.248	0.090
下层页岩	33.27	0.160	3678	2339	2620	0.056	0.024

新窗口打开| 下载CSV

表2 含气页岩双层模型2

Table 2 The double-layer model 2 of gas-bearing shales

	E/GPa	σ	V_p/ (m·s^-1)	V_s/ (m·s^-1)	ρ/ (kg·m^-3)	δ_N	δ_T
上层页岩	17.47	0.256	2960	1695	2420	0.181	0.068
下层页岩	9.57	0.324	2430	1240	2350	0.104	0.032

新窗口打开| 下载CSV

图1

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

图1 基于含气页岩双层模型1的反射系数方程精度分析结果

Fig.1 The accuracy analysis result of reflection coefficient equations based on double-layer gas-bearing shales model 1

图2

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

图2 基于含气页岩双层模型2的反射系数方程精度分析结果

Fig.2 The accuracy analysis result of reflection coefficient equations based on double-layer gas-bearing shales model 2

1.3 页岩气储层脆性及裂缝参数叠前直接反演

先验信息的加入能够为VTI介质各向异性非线性直接反演的稳定性提供保障。因此,我们将采用贝叶斯框架对反演过程施加约束。贝叶斯公式表达式如下:

(41)

$P\left(m\right|d)=\frac{P\left(d\right|m\left)P\right(m)}{P\left(d\right)}\propto P(d\left|m\right)P\left(m\right)$

式中:$P\left(m\right)$代表先验函数;$P\left(d\right|m)$为似然函数;$P(m\left|d\right)$为后验概率分布函数;d与m分别为观测数据与待反演数据;$P\left(d\right)$为常数。在叠前地震反演过程过中假设似然函数服从高斯分布,先验函数服从具有长尾巴分布的柯西约束,则似然函数与先验分布可以表示为:

(42)

$\begin{array}{l}{P}_{Gauss}\left(d\right|m)=\frac{1}{(2\pi {\sigma }_{d}^{2}{)}^{N/2}}·\\ exp\left\{\frac{-[d-G{\left(m\right)]}^{T}[d-G\left(m\right)]}{2{\sigma }_{d}^{2}}\right\}\end{array}$

(43)

${P}_{Cauchy}\left(m\right)=\frac{1}{\left(\pi \right.{\sigma }_{m}{)}^{Z}}\prod _{i=1}^{Z}\frac{1}{1+{m}_{i}^{2}/{\sigma }_{m}^{2}}$

式中:${\sigma }_{d}^{2}$与${\sigma }_{m}^{2}$分别为噪声的方差;G为正演算子;N与Z分别为叠前地震观测数据与待反演参数序列的总采样点个数。根据式(39)~(40)的描述,即可得到有效后验概率密度函数的表达式:

(44)$\begin{array}{l}P\left(m\right|d)\propto [d-G{\left(m\right)]}^{T}[d-G(m\left)\right]+\\ 2{\sigma }_{d}^{2}\sum _{i=1}^{Z}ln(1+{m}_{i}^{2}/{\sigma }_{m}^{2})\end{array}$。

为进一步提高页岩气储层展布刻画的分辨率,我们将式(41)加入L₁范数稀疏性约束,稀疏性约束能够为薄地层的识别问题及地层边界的特征突出提供解决方案,有助于提高页岩气藏的识别可靠性和稳定性,令$G(m,{\theta }_{k})=W\left({\theta }_{k}\right)·{R}_{pp}^{VTI}\left({\theta }_{k}\right)$,可得叠前非线性直接反演的目标函数J(m)如下:

(45)

$\begin{array}{l}J\left(m\right)=\sum _{k=1}^{K}[d\left({\theta }_{k}\right)-G(m,{\theta }_{k}{\left)\right]}^{T}[d\left({\theta }_{k}\right)-G(m,{\theta }_{k})]+\\ 2{\sigma }_{d}^{2}\sum _{i=1}^{Z}ln(1+{m}_{i}^{2}/{\sigma }_{m}^{2})+|{m}_{i}-{m}_{i0}|\end{array}$

式中:| |表示L₁范数;θ_k表示第k个角度下的入射角;K为入射角的个数;W(θ_k)为不同入射角下的子波矩阵。式中的其他向量的表达式展示如下:

(46)

$d\left({\theta }_{k}\right)=\left[\begin{array}{llll}{d}_{1}\left({\theta }_{k}\right)& {d}_{2}\left({\theta }_{k}\right)& \dots & {d}_{K}\left({\theta }_{k}\right)\end{array}\right]$

^T,

(47)

${R}_{pp}^{VTI}\left({\theta }_{k}\right)={\left[\begin{array}{llll}{R}_{pp,1}^{VTI}\left({\theta }_{k}\right)& {R}_{pp,2}^{VTI}\left({\theta }_{k}\right)& \dots & {R}_{pp,Z}^{VTI}\left({\theta }_{k}\right)\end{array}\right]}^{T}$

(48)

$m={\left[\begin{array}{llll}{m}_{1}& {m}_{2}& \dots & {m}_{Y+1}\end{array}\right]}^{T}$

(49)${m}_{i}={\left[\begin{array}{lllll}{E}_{b,i}& {\sigma }_{b,i}& {\rho }_{i}& {\delta }_{N,i}& {\delta }_{T,i}\end{array}\right]}^{T}$。

采用马尔科夫链蒙特卡洛(McMC)非线性优化算法对式(45)所示的非线性目标函数进行求解,在目标函数取得极小值的同时,即可得到相应页岩气储层脆性以及裂缝参数的反演结果。基于McMC算法的非线性优化求解过程描述如表3所示。

表3 McMC算法非线性优化求解流程

Table 3 Solution flow of McMC algorithm nonlinear optimization

流程	描述
步骤1	初始化算法,设置马尔可夫链的数量及每条马尔可夫链的最大循环数
步骤2	选择初始模型作为首次迭代的模型,并以此作为每条马尔可夫链的初始状态
步骤3	在每条马尔可夫链中,由特定均值和方差的建议分布生成一个新的候选模型
步骤4	根据所建立的后验概率目标函数,计算新的候选模型的接受概率
步骤5	以一定概率接受新的候选模型
步骤6	返回固定状态下的值作为单个马尔可夫链预测的结果,并以此计算每条马尔科夫链的候选模型及接受情况
步骤7	取多此循环后稳定的多条马尔科夫链的平均值作为最终的反演结果

新窗口打开| 下载CSV

2 模型测试

新方程在精度方面的可靠性将使用该方程进行叠前非线性直接反演进行页岩气储层弹性以及物性参数预测提供了可能,模型测试则进一步用于评估本文所提出的方法的可行性。我们选取了某页岩气田工区的实际测井井段进行模型测试。

图3展示了反射系数方程中所包含的5个独立的待反演岩石物理参数。图3a~e分别代表杨氏模量、泊松比、密度、裂缝参数δ_N与δ_T的模型,图中的横坐标代表各参数的数值范围,纵坐标代表采样时间。采用图3所示的岩石物理模型参数作为输入,结合式(40)以及褶积模型,正演得到了叠前无噪声与信噪比为10的正演合成地震记录,分别如图4a、b所示。图5a~e展示了以图3所示的岩石物理模型为基础,无噪声合成地震记录下的相应页岩储层岩石物理参数的叠前非线性直接反演结果。可以看出,杨氏模量与泊松比岩石弹性参数反演结果与模型具有较好的吻合度,能准确、细致反映地层的弹性性质变化。密度以及裂缝参数的反演结果相比于上述两个参数的吻合度较差,但是反演结果仍保持了较高的准确性,能够比较准确地描述地层岩石的密度变化及裂缝发育状态。为评价反演方法的稳定性与抗噪性,图6a~e展示了图3所示岩石物理模型为依托的信噪比为10的合成地震记录的反演结果。不难看出,在合成地震记录含噪声的情况下,由于噪声对输入的地震数据进行了扰动,反演结果的精度相比于无噪声情况有所下降,但是与真实模型依旧可以保持较为满意的吻合度,特别是杨氏模量、泊松比、以及裂缝参数的反演结果,都仍然能够较为准确地描述地层的岩石物理参数的变化。由于本文采用贝叶斯公式设计了叠前非线性反演的目标函数,因此反演结果的不确定性也是评价反演结果好坏的一项重要指标。图5与图6中的黑色虚线展示了相应岩石物理参数反演结果的95%置信区间,该区间包含了目标函数后验概率解的绝大部分的分布范围,不难看出,合成地震无噪声和含噪声下与置信区间都能够较为均匀地分布在有效最大后验概率解(反演结果)附近,且95%置信区间基本能够保持在较小的范围内,验证了反演方法的可靠性与稳定性。图7展示了图5与图6所示的相应岩石物理参数反演结果与真实模型的误差绝对值,其中红色曲线与蓝色曲线分别表示无噪声反演结果和信噪比为10的反演结果与真实模型之间的误差绝对值。可以观察到,依赖于方程精度的优势反演算法的鲁棒性,无论是合成记录无噪还是含有噪声,反演结果的误差都能够控制在较小的范围内,从量化的角度证明了反演方法的有效性。

图3

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

图3 待反演岩石物理参数模型

Fig.3 The rock physical models to be inverted

图4

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

图4 正演合成地震记录

Fig.4 Forward synthetic seismogram

图5

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

图5 无噪声合成地震记录反演结果

Fig.5 The inversion results at the synthesized seismogram of noise free

图6

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

图6 合成地震记录信噪比为10的反演结果

Fig.6 The inversion results at the synthesized seismogram with signal to noise ratio(SNR) equals 10

图7

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

图7 反演结果与真实模型误差绝对值

Fig.7 The absolute values of the error between the inversion results and the real models

模型测试结果所展示的反演方法在理论上的合理性、稳定性以及可靠性,为新方法应用于实际页岩储层进行岩石脆性、含裂缝性发育预测提供了理论依据。

3 实际资料测试

在模型测试的基础上,笔者将本文所提出的方法应用于中国西南某页岩气田,开展实际地震资料的反演测试工作,以评价该方法在实际工区的适用性。

图8展示了5个角度(7°、14°、19°、26°与33°)的部分叠加地震剖面,其中的盲井位于图中CDP为100的黑色竖线所指示的位置,用于对反演结果的准确定性进行评价。地震数据经过一系列处理,保证了地震振幅随入射角的真实变化关系。图9a~e分别展示了以图8所示的地震数据作为观测数据输入,采用本文设计的叠前非线性反演方法得到的杨氏模量、泊松比、密度、裂缝参数δ_N与δ_T的实际资料反演结果。在图9中,横坐标代表CDP采样点,纵坐标代表时间,每幅子图中的CDP为100位置处分别展示了相应岩石物理参数的滤波后的测井曲线。不难观察到,反演结果具有较好的分辨率,能够较好地描述地层的岩石物性、弹性纵向性质变化和横向展布状态,能够对地层的细节进行较为敏锐地刻画。从测井位置处的对比来看,表征岩石脆性的杨氏模量与泊松比反演结果能够与测井结果具有较高的吻合度,指示裂缝发育状态的裂缝参数δ_N与δ_T反演结果也能够较好地对地层的裂缝发育情况变化进行预测。密度的反演结果相比于其余参数较差,但也能与滤波后的测井曲线大致吻合。图9中椭圆框圈出的高裂缝发育程度的优质含气页岩储层范围内,反演结果展示了高杨氏模量、低泊松比以及高裂缝参数的响应,能够较好地指示页岩气的赋存,并且满足精度和分辨率的需要,为页岩储层的含气性评价奠定了基础。为进一步对比本文所提出方法与常规线性反演弹性参数进而拟合得到的页岩气储层脆性以及裂缝参数的结果的差异,图10a~d分别展示了实际地震资料叠前常规线性反演弹性参数间接拟合得到的杨氏模量、泊松比、裂缝参数δ_N与δ_T结果。可以看出,常规方法能够在一定程度上描述地层的岩石物性、弹性纵向性质变化和横向展布状态,其反演结果与测井曲线具有一定的吻合度。对比本文所提出方法与常规方法的结果,可以看出,由于新方法采用了L₁范数约束的贝叶斯框架反演目标函数,且采用McMC非线性反演算法进行了目标函数的求解,使得新方法在描述地层岩石脆性以及含裂缝性的横向展布与纵向分辨率方面相比于常规方法具有明显的优势,特别是图9中椭圆框圈出的裂缝发育程度较好的含气页岩储层的区域,新方法拥有更优秀的含气页岩指示结果以及更高分辨率的含裂缝性指示结果。与此同时,由于新方法使用了基于精确VTI方程的更高精度和更少假设条件的非线性反射系数方程,使其反演结果与测井曲线的吻合度相比于常规方法反演结果能够更贴近滤波后的测井曲线,具有更高的精度。在二维剖面反演论述的基础上,将本研究所提出的叠前非线性反演方法应用于三维工区,得到了如图10所示的实际地震资料叠前非线性反演脆性参数及裂缝参数结果的沿层切片。可以看出,高杨氏模量、低泊松比以及高裂缝参数异常的区域能够有效地刻画裂缝型含气页岩的横向展布,且裂缝型页岩储层的边界清晰。通过观察图10中的含天然气井与无天然气井位置处的岩石脆性参数与裂缝参数的反演结果,不难看出,反演结果能在含天然气井(A井、W井)的位置处指示裂缝型页岩储层的发育,且在无天然气井(S井)处不产生裂缝型页岩储层发育的假象,反演结果合理、可靠。

图8

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

图8 实际叠前过井地震资料剖面

Fig.8 The actual pre-stack cross seismic data profiles

图9

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

图9 实际地震资料叠前非线性反演(左)及常规线性反演(右)弹性参数间接拟合结果

Fig.9 The nonlinear inversion results(left) and indirect fitting results from conventional linear inversion(right) elastic parameters of actual pre-stack seismic profiles

图10

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

图10 实际地震资料叠前非线性反演脆性参数及裂缝参数结果沿层切片

Fig.10 Interlayer slices of nonlinear actual pre-stack seismic inversion results of brittleness and fractural parameters

实际资料的反演结果进一步证实了所提出方法在实际页岩气工区的适用性,利用该方法,能够为含气页岩的储层预测以及储层表征开辟新的方案。

4 结论

针对发育水平裂缝的含气页岩储层预测的需要,本文基于精确VTI反射系数方程,结合所建立的VTI介质岩石物理模型,提出了一种基于叠前反演的页岩气储层脆性及裂缝参数预测方法,该方法以贝叶斯框架为基础建立了L₁范数的非线性反演目标函数,采用McMC优化算法对页岩工区的地层杨氏模量、泊松比以及δ_N、δ_T进行了反演预测。新方法能够直接反演获取含水平裂缝的各向异性地层的岩石物理参数(杨氏模量、泊松比、密度、裂缝参数δ_N与δ_T),且有效规避线性反演以及反射率参数反演中的小角度入射、地层弱对比度的限制,模型测试与二维实际资料测试表明,相比于常规线性反演弹性参数进而拟合得到岩石物理参数,新方法提高了预测结果的分辨率与准确性。利用所提出方法在西南某页岩气田的三维反演结果的高杨氏模量、低泊松比、高裂缝参数异常能够有效预测裂缝型含气页岩储层发育的横向沿层展布状态,测井解释结果验证了预测结果准确性与可靠性,为裂缝型页岩气储层的准确预测开辟了新的方法。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

董大忠, 王玉满, 李新景,

等.

中国页岩气勘探开发新突破及发展前景思考

[J]. 天然气工业, 2016, 36(1):19-32.