CN115586573B - 一种致密砂岩储层的动态约束物性参数地震反演方法 - Google Patents

Abstract

本申请公开了一种致密砂岩储层的动态约束物性参数地震反演方法，包括：基于双重孔隙模型，设置基本岩石物理参数；基于Biot‑Rayleigh控制方程平面波解，构建地震岩石物理正演算子；利用测井观测数据，计算物性参数的混合高斯分布；利用地震道集数据及弹性参数，构建主目标函数及子目标函数；分别基于主目标函数及子目标函数，采用快速模拟退火算法更新物性参数；计算似然函数对后验权系数的偏导数，调节后验权系数；降低温度，重复步骤五至步骤六，直到达到终止温度，迭代停止，输出最终物性参数结果。

Description

一种致密砂岩储层的动态约束物性参数地震反演方法

技术领域

本发明属于非常规储层地震勘探技术领域，特别涉及致密砂岩储层的动态约束物性参数地震反演方法。

背景技术

地震反演通过地震数据定量预测岩石弹性参数，是地震数据定量处理的重要技术，在此基础上，物性参数地震反演通过结合地震反演技术与岩石物理模型，旨在从地震数据直接定量预测储层物性参数，是油气地震预测的重要途径。另一方面，致密砂岩储层作为一种重要的非常规油气资源，已占油气探明总量的近30％，是当前油气勘探增储上产的重点领域。然而，不同于常规砂岩，致密砂岩因沉积、粘土及围压等因素，其岩石内部通常呈现孔隙结构及流体分布的非均质性，其地震响应复杂，常规岩石物理模型难以适用，制约了常规物性参数地震反演方法的适用性。

常规物性参数地震反演方法通常采用Gassmann模型作为岩石物理模型，例如文献[1]-[5]公开的现有技术，但Gassmann模型基于均匀介质假设，难以合理表征致密砂岩内部的非均质性，势必对物性参数反演结果造成偏差。另一方面，Biot-Rayleigh理论模型确凿描述了非均匀介质弹性波的传播过程，是研究复杂含流体孔隙介质中地震波传播规律的重要理论[6]-[7]。相比传统Gassmann模型，Biot-Rayleigh模型能更合理地描述致密砂岩储层的地震响应特征。目前，采用Biot-Rayleigh模型的物性参数预测主要基于岩石物理图版方法[8]-[9]，但这些预测方法基于弹性阻抗参数，依赖前期叠前地震反演结果，易对物性参数结果造成累积误差。此外，基于Biot-Rayleigh模型的地震岩石物理正演算子高度非线性，需借助正则化约束方法以稳定求解反演问题，目前常规正则化方法采用固定正则参数或单一约束项，且主要针对弹性参数地震预测[10]-[11]，无法适用于物性参数地震反演问题。

综上，针对致密砂岩储层的地震反演问题，需采用更准确的岩石物理模型，并借助有效的约束反演方法，因此，开展致密砂岩储层的动态约束物性参数地震反演技术方法研究，对提高致密砂岩类非常规储层的物性参数定量预测精度具有重要意义。

上述引用非专利文献如下：

[1]李红兵,张佳佳,潘豪杰等,2021,基于弹性阻抗的孔隙结构与物性参数非线性同步反演[J].中国科学·地球科学,51,1166–1180.

[2]李志勇,张家树,蔡涵鹏等,2017,基于Hampel三截尾函数的储层弹性和物性参数同步反演[J].石油物探,56,261–272.

[3]Li,K.,Yin,X.,Zong,Z.,et al.,2022,Estimation of porosity,fluid bulkmodulus,and stiff-pore volume fraction using a multitrace Bayesian amplitude-variation-with-offset petrophysics inversion in multiporosityreservoirs[J].Geophysics,87,M25–M41.

[4]Grana,D.,Fjeldstad,T.,Omre,H.,2017.Bayesian Gaussian mixturelinear inversion for geophysical inverseproblems[J].Mathematical Geosciences,49,493–515.

[5]Bosch,M.,Mukerji,T.,Gonzalez,E.F.,2010,Seismic inversion forreservoir properties combining statisticalrockphysics andgeostatistics:Areview.Geophysics,75,165–176.

[6]巴晶,Carcione,J.M.,曹宏等,2012,非饱和岩石中的纵波频散与衰减：双重孔隙介质波传播方程.地球物理学报,55,219–231.

[7]Ba J.,Xu,W.,Fu,L.-Y.,et al.,2017,Rock anelasticity due to patchy-saturation and fabric heterogeneity:A double double-porosity model ofwavepropagation.Journal ofGeophysical Research-Solid Earth,122,1949–1971.

[8]Pang,M.,Ba,J.,Fu,L.-Y.,et al.,2020,Estimation of microfractureporosity in deep carbonate reservoirsbasedon 3Drock-physics templates[J].Interpretation,8,43–52.

[9]Picotti,S.,Carcione,J.M.,Ba,J.,2018,Rock-physics templates forseismic Q[J].Geophysics,84,MR13–MR23.

[10]Haber,E.,Oldenburg,D.W.,2000,A GCV based method for nonlinearill-posed problems[J].ComputationalGeosciences,4,41–63.

[11]Sen,M.K.,Roy,I.G.,2003,Computation of differential seismogramsand iteration adaptive regularizationinprestackwaveforminversion[J].Geophysics,68,2026–2039.

发明内容

解决的技术问题：本申请主要是提出一种致密砂岩储层的动态约束物性参数地震反演方法，解决现有技术中存在的常规岩石物理模型难以适用，难以合理表征致密砂岩内部的非均质性，势必对物性参数反演结果造成偏差，预测方法基于弹性阻抗参数，依赖前期叠前地震反演结果，易对物性参数结果造成累积误差，无法适用于物性参数地震反演问题等技术问题，旨在提高致密砂岩类非常规储层的物性参数定量预测精度。

技术方案：

一种致密砂岩储层的动态约束物性参数地震反演方法，具体包括如下步骤：

步骤一，基于双重孔隙模型，设置致密砂岩储层基本岩石物理参数；

步骤二，基于Biot-Rayleigh控制方程平面波解，得到地震道集数据及弹性参数，构建地震岩石物理正演算子；

步骤三，利用测井观测数据，计算物性参数的混合高斯分布；

步骤四，利用地震道集数据及弹性参数，构建主目标函数及子目标函数；

步骤五，分别基于主目标函数及子目标函数，采用快速模拟退火算法更新物性参数；

步骤六，计算似然函数对后验权系数的偏导数，调节后验权系数；

步骤七，降低温度，重复步骤五至步骤六，直到达到终止温度，迭代停止，输出最终物性参数结果。

作为本发明的一种优选技术方案，步骤一具体包括：Biot-Rayleigh双重孔隙模型将储层岩石内部分为背景介质与嵌入体两种结构，通过嵌入体中流体的涨缩效应描述岩石因结构及流体的非均质性对波传播的影响；基于双重孔隙模型，针对致密砂岩，设置基本岩石物理参数，包括岩石基质的体积模量、岩石基质的剪切模量、岩石基质的密度、烃类的体积模量、卤水的体积模量、流体的密度、圆孔的孔隙纵横比、裂隙的孔隙纵横比、嵌入体的半径、背景相的渗透率、嵌入体的渗透率与流体的粘度。

作为本发明的一种优选技术方案：所述步骤二具体包括：

Biot-Rayleigh控制方程是描述双重孔隙模型中波传播的波动方程，其平面波方程为

其中，

其中，A、N、Q₁、Q₂、R₁、R₂为6组不同的Biot刚度系数，ρ₁₁、ρ₁₂和ρ₁₃为背景的密度参数，ρ₂₂和ρ₃₃为嵌入体的密度参数，φ₁与φ₂分别是背景和嵌入体的孔隙度，k为波数，ω为角频率，i为虚数符号，x为衰减因子，a和b为波动方程的系数；

求解公式(1)，基于Biot-Rayleigh控制方程平面波解，计算岩石的弹性参数，将该计算过程表示为：

m＝R(r,E)+e (3)

其中，m为弹性参数，弹性参数包括纵波速度、横波速度和密度，r为物性参数，物性参数包括孔隙度、含水饱和度与泥质含量，R为Biot-Rayleigh控制方程平面波解对应的纵波速度与横波速度，E为基本岩石物理参数，e为随机误差；

将计算所得的弹性参数带入Zoeppritz方程，基于褶积模型，构建地震岩石物理正演算子，其表达式为：

d＝G(R(r,E))+e (4)

其中，d为地震道集数据，G为Zoeppritz方程正演模型。

作为本发明的一种优选技术方案：步骤三具体包括：

基于混合高斯模型，根据岩性种类数，设置高斯分布总数；

利用测井观测的孔隙度、含水饱和度与泥质含量，采用期望最大算法，计算各高斯分布的期望、协方差和权系数，获取物性参数的混合高斯分布，其概率密度函数为：

其中，r为物性参数，Г_k表示第k个标准高斯分布，μ_r|k、Σ_r|k和β_k分别为Г_k的期望、协方差和权系数，C为高斯分布总数(即岩性种类数)。

作为本发明的一种优选技术方案：步骤四具体包括：

基于贝叶斯框架，假设公式(4)中随机误差符合高斯分布，利用地震道集数据及弹性参数，结合物性参数混合高斯分布，则物性参数r对于地震道集数据d与弹性参数m的后验概率为

其中，

其中，P_k(r|m,d)表示第k个高斯分布的后验概率，const_k表示归一化常数，d为地震道集数据，m为弹性参数，r为物性参数，R为Biot-Rayleigh控制方程平面波解对应的纵波速度与横波速度，G为Zoeppritz方程正演模型，σ_m和σ_d分别为弹性参数与地震数据噪声的标准差，μ_r|k、Σ_r|k和β_k分别为第k个高斯分布的期望、协方差和权系数，C为高斯分布总数，T为转置符号；

物性参数r的反演问题等价于最大后验概率问题，根据公式(6)构建约束反演目标函数O，其一般形式为：

其中，d为地震道集数据，m为弹性参数，r为物性参数，R为Biot-Rayleigh控制方程平面波解对应的纵波速度与横波速度，G为Zoeppritz方程正演模型，σ_m和σ_d分别为弹性参数与地震数据噪声的标准差，μ_r|k、Σ_r|k和β_k分别为第k个高斯分布的期望、协方差和权系数，C为高斯分布总数，T为转置符号；

针对致密砂岩储层，令C＝2，主目标函数O_o可表示为：

O_o(r)＝||d-G(R(r))||₂+κ₁||m-R(r)||₂+κ₂[Θ₂+β₁(Θ₁-Θ₂)] (9)

Θ_k＝(r-μ_r|k)^T(Σ_r|k)^-1(r-μ_r|k),k＝1,2 (10)

其中，κ₁和κ₂为后验权系数，Θ_k为先验约束项通式，Θ₁和Θ₂为两种不同的模型先验约束项；

构建子目标函数O_i(i＝1,2,3)，其表达式分别为

O₁＝||d-G(R(r))||₂ (11)

O₂＝κ₁||m-R(r)||₂ (12)

O₃＝κ₂[Θ₂+β₁(Θ₁-Θ₂)] (13)

其中，Θ₁和Θ₂为两种不同的模型先验约束项。

作为本发明的一种优选技术方案：步骤五具体包括：

分别基于主目标函数与子目标函数，采用快速模拟退火算法更新物性参数：

其中，T^o为初始温度，T⁽ⁿ⁺¹⁾为第n+1次迭代的温度，r_j ⁽ⁿ⁺¹⁾和r_j ⁽ⁿ⁾分别为第n+1次和第n次迭代的物性参数更新值，其下标(j＝0,1,2,3)表示分别基于主目标函数(j＝0)与子目标函数(j＝1,2,3)的物性参数更新值，Δr为物性参数的扰动范围，rand₁和rand₂表示范围为[0,1]的随机数；

对于更新后的物性参数，基于Metropolis准则判别是否接收更新的物性参数，其接收概率为：

其中，O_j表示主目标函数(j＝0)与子目标函数(j＝1,2,3)。

作为本发明的一种优选技术方案：步骤六具体包括：

构建主目标函数O_o的似然函数，并计算该构建似然函数对后验权系数的偏导数

其中，κ＝[κ₁,κ₂]为后验权系数，L为主目标函数O_o的似然函数，Ω表示解空间，r为物性参数，K_z、K_e和K_β为归一化常数，O_i(i＝1,2,3)为子目标函数；

鉴于公式(16)无法解析计算，利用物性参数更新值r_j(j＝0,1,2,3)，采用后验均值法估算似然函数对后验权系数的偏导数，具体计算方式为

其中，N表示总迭代次数，r_o ⁽ⁿ⁾、r₂ ⁽ⁿ⁾和r₃ ⁽ⁿ⁾分别为基于主目标函数O_o、子目标函数O₂和子目标函数O₃的数物性参数更新值；

采用高斯-牛顿法调节后验权系数：

其中，κ⁽ⁿ⁺¹⁾和κ⁽ⁿ⁾分别为第n+1次和第n次迭代的后验权系数，δ为确保迭代收敛的常系数，T为转置符号。

作为本发明的一种优选技术方案：步骤七具体包括：

降低温度的公式为：

其中，T^o为初始温度，T⁽ⁿ⁺¹⁾为第n+1次迭代的温度，n为当前迭代次数，λ为温度衰减因子，M为模型参数维度；

重复步骤五至步骤六，直到达到终止温度，迭代停止，输出基于主目标函数的物性参数更新值，

即为最终反演结果。

作为本发明的一种优选技术方案：所述终止温度为初始温度的千分之一。

有益效果：本申请所述致密砂岩储层的动态约束物性参数地震反演方法采用的以上技术方案与现有技术相比，具有以下技术效果：

1、用于从地震数据反演储层物性参数，可以提高致密砂岩类非常规储层的物性参数定量预测精度。

2、本发明采用Biot-Rayleigh双重孔隙模型，考虑了致密砂岩储层内结构及流体的非均质影响；采用了动态约束的反演方法，通过调节后验权系数，提高反演过程的稳定性和准确性，相较于传统的物性参数地震反演方法，本发明有效提高了反演结果的精确度。

附图说明：

图1为本申请一种致密砂岩储层的动态约束物性参数地震反演方法的流程示意图；

图2为本申请实施例中地震道集数据；

图3为本申请实施例中物性参数的混合高斯分布；其中a为孔隙度的混合高斯分布图，b为含水饱和度的混合高斯分布图，c为泥质含量的混合高斯分布图；

图4为本申请实施例中利用常规方法的物性参数反演结果；

图5为本申请实施例中利用本发明方法的物性参数反演结果。

具体实施方式

下面将结合本发明的附图，对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动条件下所获得的所有其它实施例，都属于本发明保护的范围。

实施例1

如图1所示，本实施例提供了一种致密砂岩储层的动态约束物性参数地震反演方法，具体包括以下步骤：

一、构建基于Biot-Rayleigh双重孔隙模型的地震岩石物理正演算子

基于双重孔隙模型，针对致密砂岩，设置基本岩石物理参数，实施例中参数设置如下：岩石基质的体积模量K_m＝56GPa，剪切模量μ_m＝43GPa，岩石基质的密度ρ_m＝2.65g/cm³，烃类的体积模量K_h＝0.012GPa，卤水的体积模量K_w＝2.2GPa，流体的密度ρ_f＝0.82g/cm³，圆孔的孔隙纵横比α_s＝0.13，裂隙的孔隙纵横比α_c＝0.001，嵌入体的半径R_o＝0.05m，渗透率

嵌入体的渗透率

流体的粘度η＝0.001Pa·s；

Biot-Rayleigh控制方程是描述双重孔隙模型中波传播的波动方程，其平面波方程为

其中，

其中，A、N、Q₁、Q₂、R₁、R₂为6组不同的Biot刚度系数，ρ₁₁、ρ₁₂和ρ₁₃为背景的密度参数，ρ₂₂和ρ₃₃为嵌入体的密度参数，φ₁与φ₂分别是背景和嵌入体的孔隙度，k为波数，ω为角频率，i为虚数符号，x为衰减因子，a和b为波动方程的系数；

求解公式(1)，得到纵波速度与横波速度，将上述计算过程表示为

m＝R(r,E)+e (3)

其中，R为Biot-Rayleigh控制方程平面波解对应的纵波速度与横波速度，m为弹性参数(包括纵波速度、横波速度与密度)，r为物性参数(包括孔隙度、含水饱和度与泥质含量)，E为基本岩石物理参数，e为随机误差；

将计算所得的弹性参数带入Zoeppritz方程，基于褶积模型，构建地震岩石物理正演算子

d＝G(R(r,E))+e (4)

其中，d为地震道集数据，G为Zoeppritz方程正演模型；

二、计算物性参数先验概率分布及构建目标函数

利用测井观测的孔隙度、含水饱和度与泥质含量，采用期望最大算法，基于混合高斯模型，计算各高斯分布的期望、协方差和权系数，获取物性参数的混合高斯分布，其概率密度函数为

其中，r为物性参数，Г_k表示第k个标准高斯分布，μ_r|k、Σ_r|k和β_k分别为Г_k的期望、协方差和权系数，C为高斯分布总数(即岩性种类数)；实施例中，假设致密砂岩主要含砂岩与泥岩两种岩性，设置C＝2；

基于贝叶斯框架，假设公式(4)中随机误差符合高斯分布，结合物性参数混合高斯分布，则物性参数r对于地震道集数据d与弹性参数m的后验概率为

其中，

其中，P_k(r|m,d)表示第k个高斯分布的后验概率，const_k表示归一化常数，d为地震道集数据，m为弹性参数，r为物性参数，R为Biot-Rayleigh控制方程平面波解对应的纵波速度与横波速度，G为Zoeppritz方程正演模型，σ_m和σ_d分别为弹性参数与地震数据噪声的标准差，μ_r|k、Σ_r|k和β_k分别为第k个高斯分布的期望、协方差和权系数，C为高斯分布总数，T为转置符号；

物性参数r的反演问题等价于最大后验概率问题，即求取r使后验概率P(r|m,d)最大，根据公式(6)构建约束反演目标函数O，其一般形式为

其中，σ_m和σ_d和分别为弹性参数与地震数据噪声的标准差；

鉴于致密砂岩主要由砂岩与泥岩两种岩性构成，令C＝2，主目标函数O_o可表示为

O(r)＝||d-G(R(r))||₂+κ₁||m-R(r)||₂+κ₂[Θ₂+β₁(Θ₁-Θ₂)] (9)

Θ_k＝(r-μ_r|k)^T(Σ_r|k)^-1(r-μ_r|k),k＝1,2 (10)

其中，κ₁和κ₂为后验权系数，Θ_k为先验约束项通式，Θ₁和Θ₂为两种不同的模型先验约束项；

为后续调节后验权系数，构建子目标函数O_i(i＝1,2,3)，其表达式分别为

O₁＝||d-G(R(r))||₂ (11)

O₂＝κ₁||m-R(r)||₂ (12)

O₃＝κ₂[Θ₂+β₁(Θ₁-Θ₂)] (13)

其中，Θ₁和Θ₂为两种不同的模型先验约束项。

三、调节后验权系数及更新物性参数

分别基于主目标函数与子目标函数，采用快速模拟退火算法对物性参数进行更新

其中，T^o为初始温度，T⁽ⁿ⁺¹⁾为第n+1次迭代的温度，r_j ⁽ⁿ⁺¹⁾和r_j ⁽ⁿ⁾分别为第n+1次和第n次迭代的物性参数，其下标(j＝0,1,2,3)表示分别基于主目标函数(j＝0)与子目标函数(j＝1,2,3)的物性参数更新值，Δr为物性参数的扰动范围，rand₁和rand₂表示范围为[0,1]的随机数；实施例中，设置T^o＝0.85，N＝1000，Δr＝[0.002,0.09,0.03]；

对于更新后的物性参数，利用Metropolis准则判别是否接收，其接收概率为

其中，O_j表示主目标函数(j＝0)与子目标函数(j＝1,2,3)；

构建主目标函数O_o的似然函数，并计算该似然函数对后验权系数的偏导数

其中，κ＝[κ₁,κ₂]为后验权系数，L为主目标函数O_o的似然函数，r为物性参数，Ω表示解空间，K_z、K_e和K_β为归一化常数；

鉴于公式(16)无法解析计算，利用物性参数更新值r_k(k＝0,1,2,3)，采用后验均值法估算似然函数对后验权系数的偏导数，具体计算方式为

其中，N表示总迭代次数，r_o ⁽ⁿ⁾、r₂ ⁽ⁿ⁾和r₃ ⁽ⁿ⁾分别为基于主目标函数O_o、子目标函数O₂和子目标函数O₃的数物性参数更新值；

采用高斯-牛顿法调节后验权系数

其中，κ⁽ⁿ⁺¹⁾和κ⁽ⁿ⁾分别为第n+1次和第n次迭代的后验权系数，δ为确保迭代收敛的常系数；实施例中，设置δ＝0.02；

降低温度的公式为

其中，T^o为初始温度，T⁽ⁿ⁺¹⁾为第n+1次迭代的温度，n为当前迭代次数，λ为温度衰减因子，M为模型参数维度；实施例中，设置λ＝0.96，M＝3；

重复步骤五至步骤六，直到达到终止温度，终止温度为初始温度的千分之一，迭代停止，输出基于主目标函数的物性参数更新值，即为最终反演结果。

图2为实施例中的地震道集数据。图3(a～c)为实施例中物性参数的混合高斯分布，包括孔隙度、含水饱和度与泥质含量；由于不同岩性中物性参数存在统计差异，因此分布曲线呈现多峰分布特征(即多种高斯分布的叠加形态)。图4为实施例中利用常规方法(Gassmann模型)的物性参数反演结果。从图4中可见，孔隙度的反演结果(虚线)与观测数据(实线)有较大偏差，仅能反映观测数据的背景趋势；含水饱和度与泥质含量的反演结果(虚线)不稳定，与观测数据(实线)相比，难以揭示真实的物性参数特征。图5为实施例中利用本发明方法的物性参数反演结果。从图5中可见，孔隙度的反演结果(虚线)与观测数据(实线)基本一致；含水饱和度与泥质含量的反演结果(虚线)的稳定性有较明显的提高，与测井数据(实线)吻合较好。经计算所得，图5中孔隙度与含水饱和度的反演结果与观测数据的平均误差分别为0.049与0.372，而图4中的平均误差分别为0.087与0.661，由此可见，本发明方法可以有效提高物性参数反演结果的准确度。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (9)

1.一种致密砂岩储层的动态约束物性参数地震反演方法，其特征在于，具体包括如下步骤：

步骤一，基于双重孔隙模型，设置致密砂岩储层基本岩石物理参数；

步骤二，基于Biot-Rayleigh控制方程平面波解，得到地震道集数据及弹性参数，构建地震岩石物理正演算子；

步骤三，利用测井观测数据，计算物性参数的混合高斯分布；

步骤四，利用地震道集数据及弹性参数，构建主目标函数及子目标函数；

步骤五，分别基于主目标函数及子目标函数，采用快速模拟退火算法更新物性参数；

步骤六，计算似然函数对后验权系数的偏导数，调节后验权系数；

步骤七，降低温度，重复步骤五至步骤六，直到达到终止温度，迭代停止，输出最终物性参数结果。

2.根据权利要求1所述的致密砂岩储层的动态约束物性参数地震反演方法，其特征在于，步骤一具体包括：

Biot-Rayleigh双重孔隙模型将储层岩石内部分为背景介质与嵌入体两种结构，通过嵌入体中流体的涨缩效应描述岩石因结构及流体的非均质性对波传播的影响；基于双重孔隙模型，针对致密砂岩，设置基本岩石物理参数，包括岩石基质的体积模量、岩石基质的剪切模量、岩石基质的密度、烃类的体积模量、卤水的体积模量、流体的密度、圆孔的孔隙纵横比、裂隙的孔隙纵横比、嵌入体的半径、背景相的渗透率、嵌入体的渗透率与流体的粘度。

3.根据权利要求2所述的致密砂岩储层的动态约束物性参数地震反演方法，其特征在于，所述步骤二具体包括：

Biot-Rayleigh控制方程是描述双重孔隙模型中波传播的波动方程，其平面波方程为

其中，

其中，A、N、Q₁、Q₂、R₁、R₂为6组不同的Biot刚度系数，ρ₁₁、ρ₁₂、和ρ₁₃为背景介质的密度参数，ρ₂₂和ρ₃₃为嵌入体的密度参数，φ₁与φ₂分别是背景介质和嵌入体的孔隙度，k为波数，

ω为角频率，i为虚数符号，x为衰减因子，a和b为波动方程的系数；

求解公式(1)，基于Biot-Rayleigh控制方程平面波解，计算岩石的弹性参数，将该计算过程表示为：

m＝R(r, E)+e (3)

其中，m为弹性参数，弹性参数包括纵波速度、横波速度和密度，r为物性参数，物性参数包括孔隙度、含水饱和度与泥质含量，R为Biot-Rayleigh控制方程平面波解对应的纵波速度与横波速度，E为基本岩石物理参数，e为随机误差；

将计算所得的弹性参数带入Zoeppritz方程，基于褶积模型，构建地震岩石物理正演算子，其表达式为：

d＝G(R(r, E))+e (4)

其中，d为地震道集数据，G为Zoeppritz方程正演模型。

4.根据权利要求3所述的致密砂岩储层的动态约束物性参数地震反演方法，其特征在于，步骤三具体包括：

基于混合高斯模型，根据岩性种类数，设置高斯分布总数；

利用测井观测的孔隙度、含水饱和度与泥质含量，采用期望最大算法，计算各高斯分布的期望、协方差和权系数，获取物性参数的混合高斯分布，其概率密度函数为：

其中，r为物性参数，Г_k表示第k个标准高斯分布，μ_r|k、Σ_r|k和β_k分别为Г_k的期望、协方差和权系数，C为高斯分布总数即岩性种类数。

5.根据权利要求4所述的致密砂岩储层的动态约束物性参数地震反演方法，其特征在于，步骤四具体包括：

基于贝叶斯框架，假设公式(4)中随机误差符合高斯分布，利用地震道集数据及弹性参数，结合物性参数混合高斯分布，则物性参数r对于地震道集数据d与弹性参数m的后验概率为

其中，

其中，P_k(r|m,d)表示第k个标准高斯分布的后验概率，const_k表示归一化常数，d为地震道集数据，m为弹性参数，r为物性参数，R为Biot-Rayleigh控制方程平面波解对应的纵波速度与横波速度，G为Zoeppritz方程正演模型，σ_m和σ_d分别为弹性参数与地震数据噪声的标准差，μ_r|k、Σ_r|k和β_k分别为第k个标准高斯分布的期望、协方差和权系数，C为高斯分布总数，T为转置符号；

物性参数r的反演问题等价于最大后验概率问题，根据公式(6)构建约束反演目标函数O，其一般形式为：

其中，d为地震道集数据，m为弹性参数，r为物性参数，R为Biot-Rayleigh控制方程平面波解对应的纵波速度与横波速度，G为Zoeppritz方程正演模型，σ_m和σ_d分别为弹性参数与地震数据噪声的标准差，μ_r|k、Σ_r|k和β_k分别为第k个标准高斯分布的期望、协方差和权系数，C为高斯分布总数，T为转置符号；

针对致密砂岩储层，令C＝2，主目标函数O_o可表示为：

O_o(r)＝||d-G(R(r))||₂+κ₁||m-R(r)||₂+κ₂[Θ₂+β₁(Θ₁-Θ₂)] (9)

Θ_k＝(r-μ_r|k)^T(Σ_r|k)^-1(r-μ_r|k),k＝1,2 (10)

其中，κ₁和κ₂为后验权系数，Θ_k为先验约束项通式，Θ₁和Θ₂为两种不同的模型先验约束项；构建子目标函数O₁、O₂和O₃，其表达式分别为

O₁＝||d-G(R(r))||₂ (11)

O₂＝κ₁||m-R(r)||₂ (12)

O₃＝κ₂[Θ₂+β₁(Θ₁-Θ₂)] (13)

其中，Θ₁和Θ₂为两种不同的模型先验约束项。

6.根据权利要求1所述的致密砂岩储层的动态约束物性参数地震反演方法，其特征在于，步骤五具体包括：

分别基于主目标函数与子目标函数，采用快速模拟退火算法更新物性参数：

其中，T^o为初始温度，T⁽ⁿ⁺¹⁾为第n+1次迭代的温度，r_j ⁽ⁿ⁺¹⁾和r_j ⁽ⁿ⁾分别为第n+1次和第n次迭代时主目标函数与子目标函数的物性参数更新值，Δr为物性参数的扰动范围，rand₁和rand₂表示范围为[0,1]的随机数；

对于更新后的物性参数，基于Metropolis准则判别是否接收更新的物性参数，其接收概率为：

其中，O_j表示主目标函数O_o与子目标函数O₁、O₂和O₃。

7.根据权利要求1所述的致密砂岩储层的动态约束物性参数地震反演方法，其特征在于，步骤六具体包括：

构建主目标函数O_o的似然函数，并计算该构建似然函数对后验权系数的偏导数

其中，κ＝[κ₁,κ₂]为后验权系数，L为主目标函数O_o的似然函数，Ω表示解空间，r为物性参数，K_z、K_e和K_β为归一化常数，O₁、O₂和O₃为子目标函数；

鉴于公式(16)无法解析计算，利用主目标函数与子目标函数的物性参数更新值r_j，采用后验均值法估算似然函数对后验权系数的偏导数，具体计算方式为

其中，N表示总迭代次数，r_o ⁽ⁿ⁾、r₂ ⁽ⁿ⁾和r₃ ⁽ⁿ⁾分别为基于主目标函数O_o、子目标函数O₂和子目标函数O₃的数物性参数更新值；

采用高斯-牛顿法调节后验权系数：

其中，κ⁽ⁿ⁺¹⁾和κ⁽ⁿ⁾分别为第n+1次和第n次迭代的后验权系数，δ为确保迭代收敛的常系数，T为转置符号。

8.根据权利要求1所述的致密砂岩储层的动态约束物性参数地震反演方法，其特征在于，步骤七具体包括：

降低温度的公式为：

其中，T^o为初始温度，T⁽ⁿ⁺¹⁾为第n+1次迭代的温度，n为当前迭代次数，λ为温度衰减因子，M为模型参数维度；

重复步骤五至步骤六，直到达到终止温度，迭代停止，输出基于主目标函数的物性参数更新值，即为最终反演结果。

9.根据权利要求8所述的致密砂岩储层的动态约束物性参数地震反演方法，其特征在于，所述终止温度为初始温度的千分之一。

Images