E-mail Alert Rss
 

物探与化探, 2025, 49(5): 1141-1154 doi: 10.11720/wtyht.2025.0120

方法研究信息处理仪器研制

基于麻雀搜索算法的致密砂岩储层参数非线性反演方法

曹绍贺1, 黄中群1, 袁春艳1, 马百征1, 王群武2, 张奎2

1.中国石化华北油气分公司,河南 郑州 450006

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

Nonlinear inversion of tight sandstone reservoir parameters using the sparrow search algorithm

CAO Shao-He1, HUANG Zhong-Qun1, YUAN Chun-Yan1, MA Bai-Zheng1, WANG Qun-Wu2, ZHANG Kui2

1. North China Petroleum Bureau,SINOPEC,Zhengzhou 450006,China

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

第一作者: 曹绍贺(1984-),女,2009年毕业于中国石油大学(北京)地球探测与信息技术专业,获硕士学位;就职于中国石化华北油气分公司,主要从事三维地震储层及含气性预测研究工作。

责任编辑: 叶佩

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

基金资助: 中国石化集团华北石油局有限公司“彬长区块太昌目标区石盒子组地震资料目标处理与解释”项目(34550000-23-ZC0611-0064)

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

摘要

随着油气储层勘探的日益发展,常规线性反演方法已经无法满足目前致密砂岩储层的勘探精度需要。为此,本文以致密砂岩储层为勘探背景并采用包裹体模型计算了饱和岩石下的纵、横波速度及密度表达式,基于精确Zoeppritz建立非线性精确反射系数方程,以L1范数作为反演目标泛函,采用麻雀搜索算法(sparrow search algorithm,SSA),发展了一种高精度的储层参数AVO反演预测方法,相比于常规最小二乘法求解的L2范数目标函数的反演结果,该方法能够提高反演结果的精度与分辨率,为致密砂岩储层的岩石物理参数预测提供了更为合理、可靠的结果,为致密砂岩储层孔隙发育与含油气性评估提供了更为有效的解决方案。模型测试及鄂尔多斯盆地致密砂岩储层实际资料测试的结果验证了该方法的有效性和实用性。

关键词: 麻雀搜索算法; 非线性; AVO反演; 高精度; 含油气储层预测

Abstract

With the ongoing exploration of hydrocarbon reservoirs, conventional linear inversion methods for conventional hydrocarbon reservoirs fail to meet the accuracy requirements for tight sandstone reservoirs.In response,this study,focusing on tight sandstone reservoirs,calculated expressions for P-wave velocity,S-wave velocity,and density in saturated rocks using inclusion models.A nonlinear reflection coefficient equation was established based on the exact Zoeppritz equations.Using the L1-norm as the inversion objective function,a high-precision amplitude variation with offset(AVO) inversion method for reservoir parameter prediction was developed using the sparrow search algorithm(SSA).Compared to the conventional least-squares inversion results with the L2-norm as the objective function,the proposed method improves the accuracy and resolution of the inversion results.It provides more reasonable and reliable predictions of petrophysical parameters and offers an effective approach for evaluating pore development and hydrocarbon content in tight sandstone reservoirs.The validity and practicality of this method are verified through model tests and application to actual data from tight sandstone reservoirs in the Ordos Basin.

Keywords: sparrow search algorithm(SSA); nonlinear inversion; amplitude versus offset(AVO) inversion; high accuracy; prediction of hydrocarbon-bearing reservoirs

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

本文引用格式

曹绍贺, 黄中群, 袁春艳, 马百征, 王群武, 张奎. 基于麻雀搜索算法的致密砂岩储层参数非线性反演方法[J]. 物探与化探, 2025, 49(5): 1141-1154 doi:10.11720/wtyht.2025.0120

CAO Shao-He, HUANG Zhong-Qun, YUAN Chun-Yan, MA Bai-Zheng, WANG Qun-Wu, ZHANG Kui. Nonlinear inversion of tight sandstone reservoir parameters using the sparrow search algorithm[J]. Geophysical and Geochemical Exploration, 2025, 49(5): 1141-1154 doi:10.11720/wtyht.2025.0120

0 引言

利用地震反演的手段对储层岩石物理参数进行估计进而开展流体识别及储层表征是勘探地震学中的一个重要研究方向,能够为储层的含油气性评价提供明确的引导作用[1-5]

依托岩石物理理论研究,选择并构建相应的储层敏感参数岩石物理模型进而开展地球物理反演,对地下岩性以及流体类型进行判断,已经成为了目前储层描述以及含油气评价的一项重要手段。Grana等[6]利用贝叶斯框架建立了相应的储层物性参数目标函数,并提出了统计学理论下的岩石物理建模及储层物性参数反演方法,以及岩性流体聚类策略。宗兆云等[7]利用Gassmann方程推导了纵、横波模量相关的储层流体识别反射系数方程,并发展了一种叠前弹性阻抗反演储层流体敏感参数的新方法。随着地球物理理论的快速发展以及勘探工区的日益复杂,考虑砂岩致密特征的岩石物理模型及反演逐步成为了研究的焦点。卢明辉等[8]在对比和分析等效介质理论中获得了能够描述致密砂岩孔隙结构的岩石物理模型。基于Biot相洽理论,印兴耀等[9]发展了一种针对裂缝型致密砂岩的弹性参数估计方法,为致密砂岩弹性性质的分析开辟了新的途径。刘倩等[10]在考虑孔隙连通性的基础上建立了含不连通孔隙度影响下的岩石物理模型,并对其储层的弹性特征进行了研究。针对含气致密砂岩储层,Fan等[11]将统一岩石骨架模型与Sun模型相结合,得到了新的孔隙结构计算方法,并对致密砂岩孔隙结构参数进行了贝叶斯概率化反演。上述岩石物理理论驱动下的反演方法为致密砂岩储层表征提供了可靠的解决方案,但是,由于大多数反演算法为基于线性近似式的线性反演,当储层复杂性变高、纵横向变化剧烈,且勘探精度要求较高时,线性反演结果难以给出可靠的储层弹性或者物性参数预测结果。为此,寻找更为合适的反演方法以提高反演结果的准确性是目前需要解决的一项重点问题。

线性反演在保证一定的反演精度的前提下,可以为储层预测提供一种高效、快捷的指导方式[12-13]。非线性反演在极大提高反射系数方程的同时,能够进一步改善线性反演受到多种假设条件的影响,为复杂油藏的精细刻画及高精度识别提供了有效的解决方案。Yin等[14]给出了一种基于Zoeppritz方程的新型精确PP反射系数并将其应用于地震数据及叠前非线性AVO反演。该方法在非线性反问题的求解中使用了蒙特卡洛求解策略的逆算子估计算法,克服了线性方程所存在的局限性,并减少待反演参数的数量,为相对值的精确反演提供了理论支撑。Cheng等[15]推导了拉梅参数和密度相关三阶Zoeppritz非线性近似方程,该研究为反演过程引入了一种逆算子估计算法,结合L1范数的约束,利用非线性反演逆算子优化反演算法提高了拉梅参数反演稳定性和准确性。逆算子优化算法在求解L1范数提升反演分辨率与准确性方面具有较强的优势,但在求解过程中具有较低的速度,时间成本更高。Zhou等[16]在精确Zoeppritz方程的基础上,提出了一种新的Russell孔隙流体因子相关方程的非线性纵横波组合反演方法,利用梯度下降反演算法,其方法在较大程度上克服了流体因子估计精度的限制,更符合地球物理反演问题的本质,且具有较高的运算速度,时间成本较低。但梯度下降优化算法适用于L2范数的求解,难以提高反演结果的分辨率,同时,待反演参数梯度下降的过程中容易陷入局部极小值,影响反演的准确性与稳定性。非线性方程反演算法的蓬勃发展得益于反演优化算法的不断进步,因此,选择合理、稳定的优化算法开展非线性反演对致密砂岩储层进行预测,能够为提高致密砂岩储层反演结果的准确性与可靠性提供技术支撑。

考虑到上述所描述的致密砂岩岩石物理建模理论在直接反演中具有较多的参数以及较高的复杂性,不利于储层参数的稳定提取,且基于梯度类的非线性求解方法容易在优化过程中陷入局部极小值,收敛性较差。为此,本文对提高致密砂岩储层参数的叠前反演精度进行了深入研究,具体如下:首先利用岩石物理理论及Kuster-Toksöz(KT)模型[17-18]下的球形包含物模型建立纵横波速、度密度与Kf、密度以及孔隙度之间的定量关系,之后,将此关系代入精确Zoeppritz方程[19],建立含3参数的新非线性精确反射系数方程,采用L1范数建立反演目标函数,采用麻雀搜索算法[20]进行优化求解得到目标函数取得极小值时的Kf、密度以及孔隙度的数值解。该方法在精度测试和模型测试中表现出了较高的合理性和可靠性,并且在鄂尔多斯盆地的致密砂岩储层的实际应用中获得了高精度、高分辨率的勘探效果。

1 理论与方法

1.1 致密砂岩储层岩石物理模型构建

考虑到Gassmann流体替换方程[21]联系了连通孔隙低频近似下的饱和岩石骨架、基质及赋存流体之间的定量关系,作为致密砂岩岩石物理模型,由于其孔隙的不连通特征使得常规的Gassmann流体替换方程不再适用。因此Kuster-Toksöz(KT) 模型[17-18]为计算具有可变孔隙度的多孔岩石的有效体积模量Keff和有效剪切模量μeff的计算提供了一种新的解决方案:

$\left(K_{\mathrm{eff}}-K_{m}\right)\left(\frac{K_{m}+\frac{4}{3} \mu_{m}}{K_{\mathrm{eff}}+\frac{4}{3} \mu_{m}}\right)=\sum_{i=1}^{n_{i}} x_{i}\left(K_{i}-K_{m}\right) P_{i},$
$\left(\mu_{\mathrm{eff}}-\mu_{m}\right)\left(\frac{\mu_{m}+\zeta_{m}}{\mu_{\mathrm{eff}}+\zeta_{m}}\right)=\sum_{i=1}^{n_{i}} x_{i}\left(\mu_{i}-\mu_{m}\right) Q_{i},$
$\zeta_{m}=\frac{\mu_{m}}{6} \frac{9 K_{m}+8 \mu_{m}}{K_{m}+2 \mu_{m}},$

式中:PiQi分别表示不同类型包含物的几何因子;Kmμm分别表示岩石基质的体积模量和剪切模量;nx分别表示包含物的数量和体积分数;Kiμi则分别表示第i类包含物的体积模量和剪切模量。假设包含物形状为球体,且岩石孔隙由流体完全填充(Ki=Kf,μi=0),由此可得饱和致密砂岩的模量如下:

$K_{\mathrm{sat}}=\frac{4 K_{m} \mu_{m}+3 K_{f} K_{m}+4 K_{f} \mu_{m} \phi-4 K_{m} \mu_{m} \phi}{3 K_{m} \phi+4 \mu_{m}+3 K_{f}-3 K_{f} \phi},$
$\mu_{\mathrm{sat}}=\mu_{\mathrm{dry}}=\frac{\mu_{m}\left(9 K_{m}+8 \mu_{m}\right)(1-\phi)}{9 K_{m}+8 \mu_{m}+6\left(K_{m}-2 K_{m}\right) \phi},$

式中:Ksatμsat分别代表饱和流体填充孔隙下的致密砂岩体积模量和剪切模量,由式(4)和式(5)可计算饱和流体填充状态下的致密砂岩的纵横波速度表达式如下:

$V_{\mathrm{p}}=\sqrt{\frac{K_{\text {sat }}+\frac{4}{3} \mu_{\text {sat }}}{\rho}},$
$V_{\mathrm{s}}=\sqrt{\frac{\mu_{\mathrm{sat}}}{\rho}},$

式中:VpVsρ分别代表纵、横波速度和密度。上述岩石物理模型建立了孔隙流体体积模量以及孔隙度与岩石弹性参数(纵、横波速度)之间的定量联系,为后续采用非线性反射系数方程直接反演致密砂岩储层的物性以及含流体性质参数奠定了基础。

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

Zoeppritz方程[19]阐明了各向同性完全弹性假设下地震波能量在界面处的分布,并为反射系数随入射角的变化提供了精确的描述,Aki等[22]给出了纵波入射、纵波反射的反射系数表达式如下:

$R_{\mathrm{pp}}=\left(E^{-} F^{+}-G^{+} H^{-} p^{2}\right) / D$

其中,

D=E+F++G-H-p2,
$E^{+}=b \frac{\cos \alpha}{V_{\mathrm{p} 1}}+c \frac{\cos \alpha^{\prime}}{V_{\mathrm{p} 2}},$
$E^{-}=b \frac{\cos \alpha}{V_{\mathrm{p} 1}}-c \frac{\cos \alpha^{\prime}}{V_{\mathrm{p} 2}}$
$F^{+}=b \frac{\cos \beta}{V_{\mathrm{s} 1}}+c \frac{\cos \beta^{\prime}}{V_{\mathrm{s} 2}},$
$G^{+}=a+d \frac{\cos \alpha}{V_{\mathrm{p} 1}} \frac{\cos \beta^{\prime}}{V_{\mathrm{s} 2}},$
$G^{-}=a-d \frac{\cos \beta}{V_{\mathrm{p} 1}} \frac{\cos \beta^{\prime}}{V_{\mathrm{s} 2}},$
$H^{-}=a-d \frac{\cos \alpha^{\prime}}{V_{\mathrm{p} 2}} \frac{\cos \beta}{V_{\mathrm{s} 1}},$
$a=\rho_{2}\left(1-2 V_{\mathrm{s} 2}^{2} p^{2}\right)-\rho_{1}\left(1-2 V_{\mathrm{s} 1}^{2} p^{2}\right),$
$b=\rho_{2}\left(1-2 V_{\mathrm{s} 2}^{2} p^{2}\right)+2 \rho_{1} V_{\mathrm{s} 1}^{2} p^{2},$
$c=\rho_{2}\left(1-2 V_{\mathrm{s} 1}^{2} p^{2}\right)+2 \rho_{1} V_{\mathrm{s} 2}^{2} p^{2},$
$d=2 \rho_{2} V_{\mathrm{s} 2}^{2}-2 \rho_{1} V_{\mathrm{s} 1}^{2},$

式(8)~(19)中:Rpp代表反射系数;p为射线参数;αα'β以及β'分别表示纵波入射角或反射角、纵波透射角、横波反射角以及横波透射角;VpVsρ的下标1和2则分别表示所在的上层和下层介质。联立式(4)~(19),考虑基质模量Kmμm为常数,可得致密砂岩背景下的3个参数(Kfϕρ)非线性精确反射系数方程如下(反射系数方程的推导过程及具体表达见附录A):

$R_{\mathrm{pp}}=f\left(K_{f, 1}, K_{f, 2}, \phi_{1}, \phi_{2}, \rho_{1}, \rho_{2}\right),$

式中:f代表反射系数与流体体积模量Kf、孔隙度ϕ以及密度ρ之间的非线性映射关系,变量中的下标1和2分别表示反射界面所在处的上层与下层介质。反射系数方程建立的地层界面处的纵波反射系数Rpp与界面两侧的流体体积模量、孔隙度和密度之间的非线性映射,并写成了下标1和2的形式。值得注意的是,基于致密砂岩岩石物理模型推导的反射系数方程仍然包含3个独立的待反演参数,当待反演参数序列长度为N时,纵波反射系数的长度为N-1,其非线性表征了一个NN-1的非线性映射,当反射界面的采样点数加1,上一个采样点处的反射界面下层介质将变成新反射界面处的上层介质。新方程相比于线性近似的优点在于不受地层弱对比度变化以及小角度入射的使用限制,相比线性近似,能够更加准确地反映地层的真实地震响应随不同储层参数的变化。

为测试新推导的非线性反射系数方程用于非线性反演的可行性,设计了模型参数如表1所示的双层模型对其进行精度测试,模型假设各向同性完全弹性,且上下层介质分别设置为含气砂岩和泥岩,测试结果如图1所示。

表1   含气砂岩与泥岩双层模型

Table 1  The double-layer model contains gas-bearing sandstone and mudstone

Kf/PaKm/Paμm/PaVp/(km·s-1)Vs/(km·s-1)ρ/(kg·m-3)ϕ
含气砂岩1.0×1083.2×10103.4×10104.693.0929000.10
泥岩2.40×1092.0×10107×1093.491.6025000.05

新窗口打开| 下载CSV


图1

图1   反射系数方程精度测试结果

Fig.1   The precision test results of reflection coefficient equations


图1a中,黑色直线代表精确Zoeppritz,空心圆圈和黑色虚线分别代表表1所示的模型参数下,本文所推导得到的非线性精确方程和Aki-Richards线性近似的反射系数随入射角的变化关系。图1b则计算了图1a中的非线性精确方程和Aki-Richards线性近似与精确Zoeppritz的反射系数的相对误差(误差绝对值)。可以看出,由于新方程的推导源于精确Zoeppritz,未受到线性化近似对方程精度的影响,使得新推导的非线性反射系数方程的精度与精确Zoeppritz几乎一致,相比于Aki-Richards线性近似,其方程精度特别是在大角度入射情况下明显提高。该测试结果为新方程能够应用于叠前反演来获取致密砂岩储层相关参数并提高反演的分辨率的可行性提供了理论支撑。

1.3 麻雀搜索算法与非线性反演

利用褶积模型对地震数据进行正演模拟,能够生成合成地震记录。褶积模型描述如下:

$d_{t}(\theta)=w_{t}(\theta) * R_{\mathrm{pp}, t}(\theta), $

式中:θ为入射角;*表示褶积算子;dt(θ)、wt(θ)、Rpp,t(θ)分别表示合成地震记录、地震子波以及反射系数。将式(21)所示的褶积运算写成矩阵形式:

$S=W \cdot \boldsymbol{R}_{\mathrm{pp}}+\text { noise }, $

式中:Rpp为地层不同反射界面处的非线性反射系数;S为地震记录向量;W为子波矩阵;noise表示噪声项。基于褶积模型构建观测数据与正演合成数据的残差项,并结合正则化约束进一步提高算法的鲁棒性,降低优化求解过程中出现过拟合现象,同时为保障反演结果的分辨率,我们将采用强调求解结果稀疏性L1范数[23-25]建立反演的目标函数:

J=S-W·Rpp|+m-m0|

式中:||代表L1范数,mm0分别为3个待反演参数(Kfϕρ)以及其初始模型值向量,其详细表达形式如下:

Rpp=[Rpp,1(θ1),Rpp,2(θ1),,Rpp,N-1(θ1),Rpp,1(θ2),Rpp,2(θ2),,Rpp,N-1(θ2),Rpp,1(θ3),Rpp,2(θ3),,Rpp,N-1(θ3)]T
S=[S1(θ1),S2(θ1),,SN-1(θ1),S1(θ2),S2(θ2),,SN-1(θ2),S1(θ3),S2(θ3),,SN-1(θ3)]T
$\boldsymbol{W}=\left[\begin{array}{ccccccccccc}w_{1}\left(\theta_{1}\right) & 0 & \cdots & 0 & & & & & & & \\w_{2}\left(\theta_{1}\right) & w_{1}\left(\theta_{1}\right) & \cdots & 0 & & & & & & & \\\vdots & \vdots & \ddots & \vdots & & & & & & & \\w_{N / 2}\left(\theta_{1}\right) & w_{N / 2-1}\left(\theta_{1}\right) & \cdots & w_{1}\left(\theta_{1}\right) & & & & & & & \\0 & w_{N / 2}\left(\theta_{1}\right) & \cdots & w_{2}\left(\theta_{1}\right) & & & & & & & \\\vdots & \vdots & \ddots & \vdots & & & & & & & \\0 & 0 & \cdots & w_{N-1}\left(\theta_{1}\right) & & & & & & & \\& & & & w_{1}\left(\theta_{2}\right) & 0 & \cdots & 0 & & & \\& & & & w_{2}\left(\theta_{2}\right) & w_{1}\left(\theta_{2}\right) & \cdots & 0 & & & \\& & & & \vdots & \vdots & \ddots & \vdots & & & \\& & & & w_{N / 2}\left(\theta_{2}\right) & w_{N / 2-1}\left(\theta_{2}\right) & \cdots & w_{1}\left(\theta_{2}\right) & & & \\& & & & 0 & w_{N / 2}\left(\theta_{2}\right) & \cdots & w_{2}\left(\theta_{2}\right) & & & \\& & & & \vdots & \vdots & \ddots & \vdots & & & \\& & & & 0 & 0 & \cdots & w_{N-1}\left(\theta_{2}\right) & & & \\& & & & & & & & w_{1}\left(\theta_{3}\right) & 0 & \cdots & 0 \\& & & & & & & & w_{2}\left(\theta_{3}\right) & w_{1}\left(\theta_{3}\right) & \cdots & 0 \\& & & & & & & & \vdots & \vdots & \ddots & \vdots \\& & & & & & & & w_{N / 2}\left(\theta_{3}\right) & w_{N / 2-1}\left(\theta_{3}\right) & \cdots & w_{1}\left(\theta_{3}\right) \\& & & & & & & & 0 & w_{N / 2}\left(\theta_{3}\right) & \cdots & w_{2}\left(\theta_{3}\right) \\& & & & & & & & \vdots & \vdots & \ddots & \vdots \\& & & & & & & & 0 & 0 & \cdots & w_{N-1}\left(\theta_{3}\right)\end{array}\right]$
m=[Kf,1,Kf,2,,Kf,N,ϕ1,ϕ2,,ϕN,ρ1,ρ2,,ρN]T,
m0=[Kf0,1,Kf0,2,,Kf0,N,ϕ0,1,ϕ0,2,,ϕ0,N,ρ0,1,ρ0,2,,ρ0,N]T

考虑到式(23)所示的目标函数具有强非线性及涉及到L1范数的求解,使得常规梯度类型的反演算法难以获得可靠的数值解,为此,我们采用麻雀搜索算法(SSA)来求解此类非线性L1范数问题。相比于常规梯度下降算法,SSA优化求解算法作为一种启发式算法,其不受目标函数可微、可导及连续性的限制,可以有效避免搜索陷入局部极小值,同时,其能够优化L1范数目标函数,具有较高的鲁棒性和较好的收敛性,为得到分辨率更高、准确性更好的孔隙流体体积模量、孔隙度以及密度反演结果提供了算法基础。求解过程阐述如下:

令不同待反演参数序列md=[Kf,1,Kf,2,…,Kf,N,ϕ1,ϕ2,…,ϕN,ρ1,ρ2,…,ρN]=[md,1,md,2,…,md,l]为不同种群,目标函数Jd(md,1,md,2,…,md,l)=Jd(|S-WRpp|+|m-m0|)作为相应种群的适应度(目标函数越小代表适应度越强)。

首先,初始化麻雀种群,设定种群数量为h,将待反演参数种群矩阵m*中的各分量mi,j看作麻雀,建立多个含不同待反演参数md首次迭代值的位置矩阵m*:

m*=mT1mT2mTh=m1,1m1,2m1,lm2,1m2,2m2,lmh,1mh,2mh,l,

式中:l为待反演参数序列的长度,其长度l=3N。基于目标函数J,可建立初始麻雀种群的适应度矩阵如下:

J*=J1(m1,1m1,2m1,l)J2(m2,1m2,2m1,l)Jh(mh,1mh,2mh,l)

将初始化得到的麻雀种群按照适应度从大到小进行排序后,开始进行麻雀种群的更新,优化过程首先对发现者进行更新,公式如下:

$m_{i, j}^{t+1}=\left\{\begin{array}{cl}m_{i, j}^{t} \cdot \exp \left(\frac{-i}{a \cdot \text { iter }_{\max }}\right) & \text { if } R_{2}<S T, \\m_{i, j}^{t}+\Gamma \cdot \boldsymbol{L} & \text { if } \mathrm{R}_{2} \geqslant \mathrm{ST}\end{array}\right.$

式中:titermax分别代表当前迭代和最大迭代次数;a为0~1之间的随机数;Γ为基于高斯分布的随机数;L为长度为l的各元素数值为1的列向量;R2ST则分别为介于0~1以及0.5~1之间的“警告值”和“安全阈值”常数。之后,将对跟随着的位置进行更新,迭代公式如下:

$m_{i, j}^{\iota+1}=\left\{\begin{array}{cc}\Gamma \cdot \exp \left(\frac{m_{\text {worst }}-m_{i, j}^{t}}{i^{2}}\right) & \text { if } i>d / 2 \\m_{p}^{\iota+1}+\left|m_{i, j}^{\iota}-m_{p}^{\iota+1}\right| \cdot \boldsymbol{A}^{+} \cdot \boldsymbol{L} & \text { 其他。 }\end{array}\right.$

式中:mpmworst分别为发现者中的最佳位置以及全局的最差位置;A+=AT(AAT)-1,为每个元素被随机分配为1或-1的长度为l的列向量A的广义逆。最后,基于两类种群的更新,将进行随机危险预警,其迭代公式如下:

$m_{i, j}^{t+1}=\left\{\begin{array}{cl}m_{\text {best }}^{t}+b \cdot\left|m_{i, j}^{t}-m_{\text {best }}^{t}\right| & \text { if } J_{i}>J_{g} \\m_{i, j}^{t}+\mathrm{K} \cdot\left(\frac{\left|m_{i, j}^{t}-m_{\text {worst }}^{t}\right|}{\left(J_{i}-J_{w}\right)+\varepsilon}\right) & \text { 其他。 }\end{array}\right.$

式中:Ji,Jg以及Jw分别为当前、全局最优以及最差的适应度值;K为-1~1之间的随机数;b为满足均值为0方差为1的高斯分布的随机数;ε作为一个极小的常数,防止分母为零。利用式(23)所建立的目标函数,联合式(29)、(30),对式(31)~(33)进行多次迭代,直至满足输出条件(最佳适应度Jg收敛至接近零值或达到最大迭代次数),即可求解得到目标函数取得极小时的全局最佳适应度Jg所对应的待反演参数mTg=[Kf,1,Kf,2,…,Kf,N,ϕ1,ϕ2,…,ϕN,ρ1,ρ2,…,ρN]=[mg,1,mg,2,…,mg,l]的数值解。上述基于SSA算法进行的叠前非线性反演的流程图将由附录B中的图B-1所示。

2 模型测试

本文将采用图2所示的岩石物理模型,对所设计的反演方法进行模型测试,以此验证新方法所得的反演结果的合理性和有效性。

图2

图2   岩石物理模型

Fig.2   The rock physical models


图2展示了用于模型测试的Kf、孔隙度以及密度3个参数的岩石物理模型。图3则展示了模型测试所用的无噪和信噪比为5的正演合成地震记录。图4图5分别展示了以图3a图3b为观测数据的无噪和信噪比为5的地震记录下不同反演方法下的3个参数反演结果。在图4图5中,红色曲线代表真实模型,青色曲线代表初始模型,蓝色曲线和绿色曲线分别代表新非线性反演方法和常规梯度下降反演方法(梯度下降反演算法的目标函数及算法流程见附录C)下的反演结果。

图3

图3   正演合成的地震记录

Fig.3   The synthesized seismogram


图4

图4   无噪声合成地震记录下的反演结果

Fig.4   The inversion results of the synthesized seismogram with noise free


图5

图5   合成地震记录信噪比为5情况下的反演结果

Fig.5   The inversion results of the synthesized seismogram with signal-to-noise ratio of 5


不难看出,合成地震记录无噪情况下的反演结果与真实模型的吻合度较高,其中新非线性方法相比于常规方法,其结果与真实模型更为接近,此结果证明了所推导的非线性方程拥有较高的精度,能够用于合理有效地反演致密砂岩储层的岩石物理参数,且新方法由于采用了L1范数约束下的麻雀搜索算法进行求解,其结果相比于常规方法更接近真实值,进一步证实了反演算法对于进一步提升反演精度的有效性。通过观察合成地震记录信噪比为5的含噪声反演结果可以看出新提出的反演算法相比于常规方法的准确性更高,具有更好的抗噪性,证实了新方法在提升含噪声复杂影响下的反演结果可靠性较高。

为进一步量化反演结果与真实模型之间的吻合度,图6图7分别展示了合成地震记录无噪声反演结果与合成地震记录信噪比为5时的反演结果与真实模型之间的误差绝对值。通过观察能够发现,采用新方法的无噪声合成地震记录反演结果其误差绝对值基本小于常规方法,随着合成地震记录噪声的增加,反演结果的误差将增大,但由于新方法采用了L1范数约束以及具有高鲁棒性的SSA优化求解算法,其反演结果的绝对误差小于常规方法,反演稳定性更强。

图6

图6   无噪声合成地震记录下反演结果与真实模型的误差绝对值

Fig.6   The absolute error between the inversion results and the true model under noise-free synthetic seismogram


图7

图7   合成地震记录信噪比为5情况下反演结果与真实模型的误差绝对值

Fig.7   The absolute error between the inversion results and the true model under a synthetic seismogram with a signal-to-noise ratio of 5


模型测试结果通过对比不同反演优化方法,不仅验证了所推导的反射系数方程的精度,更证实了所采用的L1范数约束的麻雀搜索算法在提高反演准确性与稳定性方面的优势。为所提出反演方法应用于致密砂岩储层实际地震资料的非线性反演及岩石物理模型参数的准确提取奠定了理论基础。

3 实际资料测试

为测试所提出的非线性反演方法的实用性和可靠性,本文将选用中国中部鄂尔多斯盆地某致密砂岩勘探工区实际资料开展反演测试研究。

图8展示了用于实际地震资料测试的3个角度部分叠加的过井地震剖面,地震剖面经过精细处理,保证了其能够反映振幅随入射角变化的纯波剖面。图中的黑色竖线位置代表测井所在位置。图9a~c分别展示了L1范数约束下的以麻雀搜索算法优化求解得到的Kf、孔隙度和密度的非线性反演结果。图10a~c分别展示了采用常规梯度下降求解算法进行反演所得到的Kf、孔隙度以及密度的反演结果。在图9图10中,图中的测井所在位置分别展示了相应的滤波后的变密度显示的测井曲线。可以看出,3个参数的反演结果在两种方法下都能够在测井所在位置与测井结果有较高的吻合度。对比两种反演方法,新方法所得到的反演结果能在箭头所指示的位置处对地层参数的变化描述表现出更好的效果,其反演结果的连续性、准确性以及分辨率都优于常规方法,对于复杂情况下的薄储层能够给出更好的展布刻画效果,能够改善现有方法勘探精度方面的不足,从而提升含油气储层的识别精度。为量化不同方法反演结果之间的差异,图11展示了不同反演方法下实际地震资料井旁道反演结果与滤波后测井曲线的对比情况。不难看出,基于L1范数约束的SSA优化算法得到的非线性反演结果相比于常规梯度下降算法的反演结果与滤波后的测井曲线吻合度更高,反演结果在储层参数预测精度方面更具优势。表2统计了利用式(27)计算的实际资料井旁道反演结果与滤波后测井曲线的归一化均方根误差(normalized root mean square error,NRMSE):

$N R M S E=\frac{\sqrt{(\text { 反演值 }- \text { 滤波后井曲线值 })^{\mathrm{T}} \times(\text { 反演值 }- \text { 滤波后井曲线值 }) / N}}{\max (\text { 滤波后井曲线值 })-\min (\text { 滤波后井曲线值 })}$。

图8

图8   部分叠加的实际地震剖面

Fig.8   The partial stacked real seismic profiles


图9

图9   新非线性方法下的实际地震资料反演结果

Fig.9   The inversion results of real seismic data under the new nonlinear method


图10

图10   常规方法下的实际地震资料反演结果

Fig.10   The inversion results of real seismic data under the conventional method


图11

图11   不同反演方法下实际地震资料井旁道反演结果与滤波后测井曲线对比

Fig.11   Comparison of the well bypass inversion results of actual seismic data under different inversion methods and the filtered logging curves


表2   实际资料井旁道反演结果与滤波后测井曲线归一化均方根误差统计

Table 2  Statistics of the normalized root mean square error between the inversion results of the sidetrack of the actual data and the filtered logging curve

Kf/Paϕρ/(kg·m-3)
新方法0.08090.07180.0747
常规方法0.12990.08150.0885

新窗口打开| 下载CSV


可以看出,无论是新SSA反演方法还是常规梯度下降反演方法,其反演结果与滤波后的测井曲线都具有相对较小的归一化均方根误差值,验证了所建立的岩石物理模型的可靠性以及所推导的反射系数方程的准确性。对比新方法、常规方法与滤波后测井曲线的归一化均方根误差,新方法对孔隙流体体积模量Kf、孔隙度ϕ以及密度ρ这3个储层参数的预测结果相比于常规方法均具有更小的均方根误差,反演结果的精度更高,验证了本研究所提出的反演算方法在提升反演准确性方面的优势,为提升储层表征的合理性奠定了基础。

上述不同反演方法间的对比分析证实了所提方法在实际工区应用是有效的以及新方法在提升反演结果分辨率和准确性方面的优势,本文进一步将新方法应用于鄂尔多斯盆地三维致密砂岩工区,以测试其在三维工区的应用效果。图12图13分别给出了新方法在三维工区反演的Kf与孔隙度结果的沿层切片。在该三维工区中,有3口盲井(井A、井B以及井C),其中井A与井C在该层段测井解释含油,而井B在该层段不含油。通过观察切片结果的低Kf与高孔隙度异常位置,能够清晰反映含油砂岩储层的横向展布预测,预测结果与测井所在位置是否含油气解释结果能够吻合,为利用本文方法进行致密砂岩储层精确预测给予了新的方案。

图12

图12   新非线性方法下的Kf反演结果沿层切片

Fig.12   The interlayer slice of the real seismic data Kf inversion result under the new nonlinear method


图13

图13   新非线性方法下的孔隙度反演结果沿层切片

Fig.13   The interlayer slice of the real seismic data porosity inversion result under the new nonlinear method


实际地震资料的反演结果首先证明了所建立的岩石物理模型以及所推导的反射系数方程能够用于实际工区反演得到合理、可靠的反演结果。其次,通过不同反演方法的对比,证明了新方法能够进一步提高反演结果的精度,进而满足鄂尔多斯盆地致密砂岩储层实际工区勘探精度的需要。该方法在鄂尔多斯盆地致密砂岩储层实际二维及三维工区的应用中,利用所预测的储层参数得到了准确的致密砂岩相对高孔隙度含油气储层的分布,为该方法应用于鄂尔多斯盆地致密砂岩工区进行储层表征奠定了方法基础。

4 结论

本文针对致密砂岩储层背景下勘探精度与分辨率提升的需要,采用Kuster-Toksöz (KT)模型联系了饱和岩石的纵横波速度、密度以及储层岩石物理参数,基于精确Zoepprtiz,推导得到了非线性精确反射系数方程,方程相比于Aki-Richards线性近似实现了反射系数精度的提升,并联系了反射系数与储层孔隙流体体积模量与孔隙度,为储层物性参数的直接反演预测提供理论依据。

以褶积模型为依托,采用L1范数约束下的麻雀搜算算法作为求解方法,发展了一种高精度非线性AVO反演方法,该方法相比于常规梯度下降算法能够进一步提升致密砂岩储层孔隙流体体积模量、孔隙度和密度反演结果的精度和分辨率,为提高储层表征的可靠性提供了技术基础。

将所提出的方法应用在鄂尔多斯盆地致密砂岩工区,得到了合理、可靠的孔隙流体体积模量与孔隙度预测结果,并在复杂致密砂岩的相对高孔隙度含油储层预测与展布刻画中得到了可靠的预测结果,为提升复杂致密砂岩储层实际勘探的准确性和分辨率给予了新的方案。

附录A

含孔隙流体体积模量Kf、孔隙度ϕ以及密度ρ的非线性三参数反射系数方程的推导首先将建立纵波速度Vp、横波速度Vs与3个待反演参数之间的非线性联系,将式(4)、(5)代入式(6)、(7),可得:

Vp=4Kmμm+3KfKm+4Kfμmϕ-4Kmμmϕ3ρ(Kmϕ+4μm+3Kf-3Kfϕ)+4μm9Km+8μm)(1-ϕ)3ρ9Km+8μm+6(Km-2Km)ϕ),
Vs=μm9Km+8μm)(1-ϕ)9ρ(Km+8μm+6(Km-2Km)ϕ),

考虑反射系数方程中反射界面两侧的上、下层介质,将式(A-1)、(A-2)代入式(8)~(19),可得含Kfϕρ的非线性三参数反射系数方程如下:

Rpp=f(Kf,1,Kf,2,ϕ1,ϕ2,ρ1,ρ2)=(E-F+-G+H-p2)/D,
D=E+F++G-H-p2,
$E^{+}=b \frac{\cos \alpha}{V_{\mathrm{p} 1}}+c \frac{\cos \alpha^{\prime}}{V_{\mathrm{p} 2}},$
$E^{-}=b \frac{\cos \alpha}{V_{\mathrm{p} 1}}-c \frac{\cos \alpha^{\prime}}{V_{\mathrm{p} 2}}$
$F^{+}=b \frac{\cos \beta}{V_{\mathrm{s} 1}}+c \frac{\cos \beta^{\prime}}{V_{\mathrm{s} 2}},$
$G^{+}=a+d \frac{\cos \alpha}{V_{\mathrm{p} 1}} \frac{\cos \beta^{\prime}}{V_{\mathrm{s} 2}},$
$G^{-}=a-d \frac{\cos \beta}{V_{\mathrm{p} 1}} \frac{\cos \beta^{\prime}}{V_{\mathrm{s} 2}},$
$H^{-}=a-d \frac{\cos \alpha^{\prime}}{V_{\mathrm{p} 2}} \frac{\cos \beta}{V_{\mathrm{s} 1}},$
$a=\rho_{2}\left(1-2 V_{\mathrm{s} 2}^{2} p^{2}\right)-\rho_{1}\left(1-2 V_{\mathrm{s} 1}^{2} p^{2}\right),$
$b=\rho_{2}\left(1-2 V_{\mathrm{s} 2}^{2} p^{2}\right)+2 \rho_{1} V_{\mathrm{s} 1}^{2} p^{2},$
$c=\rho_{2}\left(1-2 V_{\mathrm{s} 1}^{2} p^{2}\right)+2 \rho_{1} V_{\mathrm{s} 2}^{2} p^{2},$
$d=2 \rho_{2} V_{\mathrm{s} 2}^{2}-2 \rho_{1} V_{\mathrm{s} 1}^{2},$
$V_{\mathrm{p} 1}=\sqrt{\frac{4 K_{m} \mu_{m}+3 K_{f 1} K_{m}+4 K_{f 1} \mu_{m} \phi_{1}-4 K_{m} \mu_{m} \phi_{1}}{3 \rho_{1}\left(K_{m} \phi_{1}+4 \mu_{m}+3 K_{f 1}-3 K_{f 1} \phi_{1}\right)}+\frac{4 \mu_{m}\left(9 K_{m}+8 \mu_{m}\right)\left(1-\phi_{1}\right)}{3 \rho_{1}\left[9 K_{m}+8 \mu_{m}+6\left(K_{m}-2 K_{m}\right) \phi_{1}\right]}},$
$V_{\mathrm{p} 2}=\sqrt{\frac{4 K_{m} \mu_{m}+3 K_{f 2} K_{m}+4 K_{f 2} \mu_{m} \phi_{2}-4 K_{m} \mu_{m} \phi_{2}}{3 \rho_{2}\left(K_{m} \phi_{2}+4 \mu_{m}+3 K_{f 2}-3 K_{f 2} \phi_{2}\right)}+\frac{4 \mu_{m}\left(9 K_{m}+8 \mu_{m}\right)\left(1-\phi_{2}\right)}{3 \rho_{2}\left[9 K_{m}+8 \mu_{m}+6\left(K_{m}-2 K_{m}\right) \phi_{2}\right]}},$
$V_{\mathrm{s} 1}=\sqrt{\frac{\mu_{m}\left(9 K_{m}+8 \mu_{m}\right)\left(1-\phi_{1}\right)}{9 \rho_{1}\left[K_{m}+8 \mu_{m}+6\left(K_{m}-2 K_{m}\right) \phi_{1}\right]}},$
$V_{\mathrm{s} 2}=\sqrt{\frac{\mu_{m}\left(9 K_{m}+8 \mu_{m}\right)\left(1-\phi_{2}\right)}{9 \rho_{2}\left[K_{m}+8 \mu_{m}+6\left(K_{m}-2 K_{m}\right) \phi_{2}\right]}},$

式中:Kf1ϕ1ρ1分别表示反射界面上层的孔隙流体体积模量、孔隙度和密度;Kf2ϕ2ρ2分别表示反射界面下层的孔隙流体体积模量、孔隙度和密度。

附录B

基于SSA算法的叠前非线性反演孔隙流体体积模量Kf、孔隙度ϕ以及密度ρ的算法流程图描述如图B-1

图B-1

图B-1   基于SSA算法的叠前非线性反演流程

Fig.B-1   The flow chart of pre-stack nonlinear inversion on the basies of SSA algorithm


附录C

依托式(22)所示的褶积模型,基于褶积模型构建观测数据与正演合成数据的残差项,并采用常规L2范数正则化约束进一步提高算法的抗噪性并降低优化求解过程中出现过拟合现象,建立基于梯度下降反演算法的目标函数建立如下:

JG=S-W·Rpp+m-m0,

式中:‖‖代表L2范数, mm0分别表示3个待反演参数(Kf,ϕ,ρ)以及其初始模型值向量,式(C-1)中的各项详细表达描述如下:

$\begin{aligned}\boldsymbol{R}_{\mathrm{pp}}= & {\left[R_{\mathrm{pp}, 1}\left(\theta_{1}\right), R_{\mathrm{pp}, 2}\left(\theta_{1}\right), \cdots, R_{\mathrm{pp}, N-1}\left(\theta_{1}\right), R_{\mathrm{pp}, 1}\left(\theta_{2}\right), R_{\mathrm{pp}, 2}\left(\theta_{2}\right), \cdots, R_{\mathrm{pp}, N-1}\left(\theta_{2}\right), R_{\mathrm{pp}, 1}\left(\theta_{3}\right), R_{\mathrm{pp}, 2}\left(\theta_{3}\right),\right.} \\ & \left.\cdots, R_{\mathrm{pp}, N-1}\left(\theta_{3}\right)\right]^{\mathrm{T}},\end{aligned}$
$\boldsymbol{S}=\left[\mathrm{S}_{1}\left(\theta_{1}\right), \mathrm{S}_{2}\left(\theta_{1}\right), \cdots, \mathrm{S}_{N-1}\left(\theta_{1}\right), \mathrm{S}_{1}\left(\theta_{2}\right), \mathrm{S}_{2}\left(\theta_{2}\right), \cdots, \mathrm{S}_{N-1}\left(\theta_{2}\right), \mathrm{S}_{1}\left(\theta_{3}\right), \mathrm{S}_{2}\left(\theta_{3}\right), \cdots, \mathrm{S}_{N-1}\left(\theta_{3}\right)\right]^{\mathrm{T}},$
$W=\left[\begin{array}{ccccccccccc}\\w_{1}\left(\theta_{1}\right) & 0 & \cdots & 0 & & & & & & & \\w_{2}\left(\theta_{1}\right) & w_{1}\left(\theta_{1}\right) & \cdots & 0 & & & & & & & \\\vdots & \vdots & \ddots & \vdots & & & & & & & \\w_{N / 2}\left(\theta_{1}\right) & w_{N / 2-1}\left(\theta_{1}\right) & \cdots & w_{1}\left(\theta_{1}\right) & & & & & & & \\0 & w_{N / 2}\left(\theta_{1}\right) & \cdots & w_{2}\left(\theta_{1}\right) & & & & & & & \\\vdots & \vdots & \ddots & \vdots & & & & & & & \\0 & 0 & \cdots & w_{N-1}\left(\theta_{1}\right) & & & & & & & \\& & & & w_{1}\left(\theta_{2}\right) & 0 & \cdots & 0 & & & \\& & & & w_{2}\left(\theta_{2}\right) & w_{1}\left(\theta_{2}\right) & \cdots & 0 & & & \\& & & & \vdots & \vdots & \ddots & \vdots & & & \\& & & & w_{N / 2}\left(\theta_{2}\right) & w_{N / 2-1}\left(\theta_{2}\right) & \cdots & w_{1}\left(\theta_{2}\right) & & & \\& & & & 0 & w_{N / 2}\left(\theta_{2}\right) & \cdots & w_{2}\left(\theta_{2}\right) & & & \\& & & & \vdots & \vdots & \ddots & \vdots & & & \\& & & & 0 & 0 & \cdots & w_{N-1}\left(\theta_{2}\right) & & & \\& & & & & & & w_{1}\left(\theta_{3}\right) & 0 & \cdots & 0\\& & & & & & & w_{2}\left(\theta_{3}\right) & w_{1}\left(\theta_{3}\right) & \cdots & 0 \\& & & & & & & \vdots & \vdots & \ddots & \vdots \\& & & & & & & w_{N / 2}\left(\theta_{3}\right) & w_{N / 2-1}\left(\theta_{3}\right) & \cdots & w_{1}\left(\theta_{3}\right) \\& & & & & & & 0 & w_{N / 2}\left(\theta_{3}\right) & \cdots & w_{2}\left(\theta_{3}\right) \\& & & & & & &\vdots & \vdots & \ddots & \vdots \\& & & & & & & 0 & 0 & \cdots & w_{N-1}\left(\theta_{3}\right)\end{array}\right]$
m=[Kf,1,Kf,2,,Kf,N,ϕ1,ϕ2,,ϕN,ρ1,ρ2,,ρN]T
m0=[Kf0,1,Kf0,2,,Kf0,N,ϕ0,1,ϕ0,2,,ϕ0,N,ρ0,1,ρ0,2,,ρ0,N]T

采用梯度下降算法[16,26]求解式(C-1)所建立的目标函数,进行多次迭代,直至满足输出条件(目标函数JG收敛至接近零值或达到最大迭代次数),即可求解得到目标函数取得极小时时对应的待反演参数m=[Kf,1,Kf,2,…,Kf,N,ϕ1,ϕ2,…,ϕN,ρ1,ρ2,…,ρN]T的数值解。梯度下降算法优化求解流程如图C-1

图C-1

图C-1   基于梯度下降算法的叠前非线性反演流程

Fig.C-1   The flow chart of pre-stack nonlinear inversion on the basies of gradient descent algorithm


参考文献

Bosch M, Carvajal C, Rodrigues J, et al.

Petrophysical seismic inversion conditioned to well-log data:Methods and application to a gas reservoir

[J]. Geophysics, 2009, 74(2):O1-O15.

[本文引用: 1]

Ran B.

Linearized orthorhombic AVAZ inversion:Theoretical and practical consideration

[C]// Denver:SEG Technical Program Expanded Abstracts 2014, Society of Exploration Geophysicists,2014:528-532.

Zong Z Y, Yin X Y, Wu G C.

Geofluid discrimination incorporating poroelasticity and seismic reflection inversion

[J]. Surveys in Geophysics, 2015, 36(5):659-681.

Pan X P, Zhang G Z.

Model parameterization and PP-wave amplitude versus angle and azimuth(AVAZ) direct inversion for fracture quasi-weaknesses in weakly anisotropic elastic media

[J]. Surveys in Geophysics, 2018, 39(5):937-964.

Zong Z Y, Zhang J L, Chen F B.

Seismic prediction for formation pressure considering diagenesis effect

[J]. IEEE Transactions on Geoscience and Remote Sensing, 2024,62:4508813.

[本文引用: 1]

Grana D, Della Rossa E.

Probabilistic petrophysical-properties estimation integrating statistical rock physics with seismic inversion

[J]. Geophysics, 2010, 75(3):O21-O37.

[本文引用: 1]

宗兆云, 印兴耀, 吴国忱.

基于叠前地震纵横波模量直接反演的流体检测方法

[J]. 地球物理学报, 2012, 55(1):284-292.

[本文引用: 1]

Zong Z Y, Yin X Y, Wu G C.

Fluid detection method based on direct inversion of pre-stack seismic P-S wave modulus

[J]. Chinese Journal of Geophysics, 2012, 55(1):284-292.

[本文引用: 1]

卢明辉, 巴晶, 晏信飞.

致密砂岩的等效介质理论研究

[C]// 中国地球物理学会, 中国石油学会,2011:725-729.

[本文引用: 1]

Lu M, Ba J, Yan X F.

Theoretical study on equivalent medium of tight sandstone

[C]// Chinese Geophysical Society, Chinese Petroleum Society,2011:725-729.

[本文引用: 1]

印兴耀, 刘欣欣, 曹丹平.

基于Biot相洽理论的致密砂岩弹性参数计算方法

[J]. 石油物探, 2013, 52(5):445-451,441.

DOI:10.3969/j.issn.1000-1441.2013.05.001      [本文引用: 1]

致密砂岩储层中速度孔隙度关系复杂,针对致密砂岩中孔隙/裂隙并存的特点,提出了一种致密砂岩储层弹性参数的计算方法。根据Biot相洽理论,引入岩石骨架泊松比,在Biot相洽条件下对其进行求取,用于岩石骨架弹性模量的计算;结合约束反演求取的裂隙孔隙度,在Biot理论框架下求取孔隙和裂隙并存时致密砂岩的弹性参数。该方法无需对除裂隙以外的孔隙进行描述,避免了孔隙参数选取不合理带来的计算误差。使用该方法计算实际致密砂岩储层的弹性参数,用于横波速度计算以及流体敏感参数交会分析,取得了较好的效果,证明了方法的有效性。

Yin X Y, Liu X X, Cao D P.

Elastic parameters calculation for tight sand reservoir based on Biotconsistent theory

[J]. Geophysical Prospecting for Petroleum, 2013, 52(5):445-451,441.

[本文引用: 1]

刘倩, 印兴耀, 李超.

含不连通孔隙的致密砂岩储层岩石弹性模量预测方法

[J]. 石油物探, 2015, 54(6):635-642.

DOI:10.3969/j.issn.1000-1441.2015.06.001      [本文引用: 1]

致密砂岩储层与常规储层相比具有低孔、低渗、连通性差等特性。针对致密砂岩储层中不连通孔隙的影响,提出了一种适用于致密砂岩储层的岩石物理模型建立方法,重点是利用Raymer公式引入不连通孔隙度的影响来修正基质模量的方法,阐述了依照连通孔隙度Gassmann方程在干岩石骨架的孔隙空间进行流体充填,分析讨论了不连通孔隙的引入对于岩石弹性模量及流体响应特征的影响。应用该方法对实验测量数据及某实际工区井资料进行试算发现,与常规有效介质模型预测结果相比,基于含不连通孔隙致密砂岩岩石物理模型估算的纵、横波速度值与测井结果吻合更好,验证了方法的合理性。

Liu Q, Yin X Y, Li C.

Prediction method of rock elastic modulus of tight sandstone reservoirs with disconnected pores

[J]. Geophysical Prospecting for Petroleum, 2015, 54(6):635-642.

[本文引用: 1]

Fan X G, Zhang G Z, Zhang J J.

Prediction method of pore structure parameters of tight sandstone

[C]// San Antonio:SEG Technical Program Expanded Abstracts 2019, Society of Exploration Geophysicists,2019:3678-3682.

[本文引用: 1]

Cooke D A, Schneider W A.

Generalized linear inversion of reflection seismic data

[J]. Geophysics, 1983, 48(6):665-676.

[本文引用: 1]

Du B Y, Yang W Y, Zhang J, et al.

Matrix-fluid decoupling-based joint PP-PS-wave seismic inversion for fluid identification

[J]. Geophysics, 2019, 84(3):R477-R487.

[本文引用: 1]

Yin X Y, Cheng G S, Zong Z Y.

Non-linear AVO inversion based on a novel exact PP reflection coefficient

[J]. Journal of Applied Geophysics, 2018,159:408-417.

[本文引用: 1]

Cheng G S, Yin X Y, Zong Z Y.

Third-order AVO inversion for lamé parameter based on inverse operator estimation algorithm

[J]. Journal of Petroleum Science and Engineering, 2018,164:117-126.

[本文引用: 1]

Zhou L, Liu X Y, Li J Y, et al.

Robust AVO inversion for the fluid factor and shear modulus

[J]. Geophysics, 2021, 86(4):R471-R483.

[本文引用: 2]

Kuster G T, Toksöz M N.

Velocity and attenuation of seismic waves in two-phase media:Part i.Theoretical formulations

[J]. Geophysics, 1974, 39(5):587-606.

[本文引用: 2]

Berryman J G.

Long-wavelength propagation in composite elastic media II.Ellipsoidal inclusions

[J]. Acoustical Society of America Journal, 1980, 68(6):1820-1831.

[本文引用: 2]

Zoeppritz K, Erdbebnenwellen V.

On the reflection and penetration of seismic waves through unstable layers

[J]. Gottinger Nachrichten, 1919, 1(5):66-84.

[本文引用: 2]

Xue J K, Shen B.

A novel swarm intelligence optimization approach:Sparrow search algorithm

[J]. Systems Science & Control Engineering, 2020, 8(1):22-34.

[本文引用: 1]

Gassmann F.

Elastic waves through a packing of spheres

[J]. Geophysics, 1951, 16(4):673-685.

[本文引用: 1]

Aki K, Richards P. Quantitative seismology:Theory and methods[M]. W.H. Freeman and Co.,1980.

[本文引用: 1]

印兴耀, 邓炜, 宗兆云.

基于逆算子估计的AVO反演方法研究

[J]. 地球物理学报, 2016, 59(4):1457-1468.

DOI:10.6038/cjg20160426      [本文引用: 1]

传统反演算法以优化算法为主,而基于逆算子估计的AVO反演算法则利用了直接求逆的思路.算法的关键在于寻找存在逆函数的子域,进而可以在子域内直接求逆,这种解决反问题的思路不同于一般的优化类算法所采用的直接搜索解的方式,具有更高的效率.AVO反演利用了振幅随着偏移距的变化特征,反演的精度受到地震资料质量的影响,通过加入L1范数约束以及合理的初始模型有助于提高反演的稳定性以及准确度.模型测算和实际应用表明,基于逆算子估计的AVO反演方法具有较高的精确程度和可靠性.

Yin X Y, Deng W, Zong Z Y.

AVO inversion based on inverse operator estimation

[J]. Chinese J. Geophys., 2016, 59(4):1457-1468.

[本文引用: 1]

邓炜, 印兴耀, 宗兆云, .

基于逆算子估计的高阶AVO非线性反演

[J]. 石油地球物理勘探, 2016, 51(5):955-964,837.

Deng W, Yin X Y, Zong Z Y, et al.

High-order nonlinear AVO inversion based on estimated inverse operator

[J]. Oil Geophysical Prospecting, 2016, 51(5):955-964,837.

Cheng G S, Yin X Y, Zong Z Y.

Nonlinear amplitude-variation-with-offset inversion for Lamé parameters using a direct inversion method

[J]. Interpretation, 2017, 5(3):SL57-SL67.

[本文引用: 1]

Luo C, Li X Y, Huang G T.

Pre-stack AVA inversion by using propagator matrix forward modeling

[J]. Pure and Applied Geophysics, 2019, 176(10):4445-4476.

DOI:10.1007/s00024-019-02157-9      [本文引用: 1]

Most existing amplitude variation with angle (AVA) inversions are based on the exact Zoeppritz equation or its approximations. These modeling methods, which are ray-tracing-based and describe P-wave primary reflections only, lead to exacting requirements for pre-processing of the input data. Current processing is inadequate to satisfy these demands, especially for removing the effects of transmission losses, P-wave multiples and various converted-wave modes. By using input data with processing errors, inversion results of primary-only methods are predictably not accurate enough. The propagator matrix (PM), like the reflectivity method, uses an analytical solution to the wave equation and considers full-wave propagation effects in horizontal or nearly horizontal multilayered earth models. The numerical examples verify that a PM can effectively estimate transmission losses, multi-reflections and the comprehensive responses of thin interbedded layers, and also has higher reflection sensitivities to P-wave and S-wave velocity and density, as compared with raytracing-based AVA modelling. A pre-stack AVA three-parameter inversion by using a PM as the forward engine is proposed. Following a Bayesian approach, the inversion is stabilized by including the correlation of P-wave velocity, S-wave velocity and density. For inversion accuracy, the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) optimization method is used to solve the augmented function, and the generalized cross-validation (GCV) criterion (Huang et al. in J Geophys Eng 14(1):100-112, 2017) is introduced to adaptively acquire the regularization parameter. Theoretical model inversion analysis shows that the proposed inversion can make use of transmission losses, P-wave multiples and converted wave modes, which not only cost-effectively simplifies the pre-processing, but also generates reasonable inverted results for multilayered conditions. The proposed inversion is then applied to a set of real data, and a comparison with Zoeppritz equation-based inversion demonstrates that PM inversion is clearly superior to Zoeppritz equation-based inversion in terms of stability and accuracy.

/

京ICP备05055290号-3
版权所有 © 2021《物探与化探》编辑部
通讯地址:北京市学院路29号航遥中心 邮编:100083
电话:010-62060192;62060193 E-mail:whtbjb@sina.com , whtbjb@163.com