一种孔隙参数与储层参数的地震岩石物理解析反演方法

1.本发明属于非常规储层地震勘探技术领域，特别涉及孔隙参数与储层参数的地震岩石物理解析反演方法。


背景技术：

2.地震岩石物理反演结合了地震反演技术与岩石物理模型，通过地震数据或属性定量预测储层参数，旨在为油气储层的预测、勘探及开发提供技术支撑。另一方面，随着油气勘探的不断推进，常规储层日益枯竭，致密砂岩、碳酸盐岩等非常规储层已称为增储上产的重点领域，但此类储层普遍发育复杂孔隙结构，严重制约了常规地震岩石物理反演方法的适用性。
3.常规地震岩石物理反演方法主要预测储层参数(孔隙度、含水饱和度等)，例如文献[1]-[4]公开的现有技术，但上述方法在岩石物理建模中通常将孔隙参数(孔隙纵横比)作为已知参数，并在反演过程中采用固定的孔隙纵横比。研究发现，岩石的孔隙纵横比对地震响应具有重要影响，尤其针对发育复杂孔隙结构的储层岩石[5]-[6]；此外，非常规储层由于地质构造及沉积环境等因素，储层岩石的孔隙结构异常复杂[7]-[8]，因此基于固定孔隙参数的地震岩石物理反演方法难以适用。目前，考虑孔隙参数或结构的地震岩石物理反演方法研究起步较晚，例如li等(2022)提出了一种预测硬孔孔隙百分比的地震反演方法[9]，但该方法仅定性考虑了孔隙结构的影响，无法直接预测孔隙参数，难以定量评估储层的孔隙结构；李宏兵等(2013)和teillet等(2021)提出了预测孔隙纵横比的地震反演方法[10]-[11]，虽然该方法能够定量预测孔隙参数，但无法同步预测储层参数，限制了储层地震预测的准确度；李红兵等(2021)发展了孔隙扁度与含水饱和度的同步反演方法[12]，但该反演方法依赖全局优化，求解过程效率低，并且计算结果无解析表达，存在不确定性。
[0004]
综上，针对具有复杂孔隙结构的非常规储层，需定量考虑孔隙结构的影响，因此，开展同步预测孔隙参数与储层参数的解析反演技术方法研究，对提高致密砂岩与碳酸盐岩类非常规储层参数定量预测精度具有重要意义。
[0005]
上述引用非专利文献如下：
[0006]
[1]yin,x.,sun,r.,wang,b.,zhang,g.,2014,simultaneous inversion of petrophysical parameters based on geostatistical a priori information[j].applied geophysics,11,311
–
320.
[0007]
[2]de figueiredo,l.p.,grana,d.,santos,m.,et al.,2017,bayesian seismic inversion based on rock-physics prior modeling for the joint estimation of acoustic impedance,porosity and lithofacies[j].journal of computational physics,336,128
–
142.
[0008]
[3]李志勇,张家树,蔡涵鹏等,2017,基于hampel三截尾函数的储层弹性和物性参数同步反演[j].石油物探,56(2),261
–
272.
[0009]
[4]张佳佳,印兴耀,张广智等,2020,基于线性化岩石物理反演的物性参数预测方
法[j].石油勘探与开发,47(1),57
–
64.
[0010]
[5]weger,r.j.,eberli,g.p.,baechle,g.t.,et al.,2009,quantification of pore structure and its effect on sonic velocity and permeability in carbonates[j].aapg bulletin,93,1297
–
1317.
[0011]
[6]lehocki,i.,avseth,p.,2021,from cradle to grave:how burial history controls the rock-physics properties of quartzose sandstones[j].geophysical prospecting,69,629
–
649.
[0012]
[7]李王鹏,刘忠群,胡宗全等,2021,四川盆地川西坳陷新场须家河组二段致密砂岩储层裂缝发育特征及主控因素[j].石油与天然气地质,42(4),884
–
897.
[0013]
[8]ren,w.,shi,w.,xie,x.,et al.,2020,clay mineral content,type,and their effects on pore throat structure and reservoir properties:insight from the permian tight sandstones in the hangjinqi area,north ordos basin,china[j].marine and petroleum geology,115,104281.
[0014]
[9]li,k.,yin,x.,zong,z.,et al.,2022,estimation of porosity,fluid bulk modulus,and stiff-pore volume fraction using a multitrace bayesian amplitude-variation-with-offset petrophysics inversion in multiporosity reservoirs[j].geophysics,87,m25
–
m41.
[0015]
[10]李宏兵,张佳佳,姚逢昌,2013,岩石的等效孔隙纵横比反演及其应用[j].地球物理学报,56(2),608
–
615.
[0016]
[11]teillet,t.,fournier,f.,zhao,l.,et al.,2021,geophysical pore type inversion in carbonate reservoir:integration of cores,well logs,and seismic data(yadana field,offshore myanmar)[j].geophysics,86,b149
–
b164.
[0017]
[12]李红兵,张佳佳,潘豪杰等,2021,基于弹性阻抗的孔隙结构与物性参数非线性同步反演[j].中国科学
·
地球科学,51,1166
–
1180.


技术实现要素：

[0018]
解决的技术问题：本技术主要是提出一种孔隙参数与储层参数的地震岩石物理解析反演方法，解决现有技术中存在的致密砂岩、碳酸盐岩储层普遍发育复杂孔隙结构，严重制约了常规地震岩石物理反演方法的适用性，非常规储层由于地质构造及沉积环境等因素，储层岩石的孔隙结构异常复杂，因此基于固定孔隙参数的地震岩石物理反演方法难以适用；仅定性考虑了孔隙结构的影响，无法直接预测孔隙参数，难以定量评估储层的孔隙结构；无法同步预测储层参数，限制了储层地震预测的准确度；反演方法依赖全局优化，求解过程效率低，并且计算结果无解析表达，存在不确定性等技术问题，旨在提高致密砂岩与碳酸盐岩类非常规储层参数定量预测精度。
[0019]
技术方案：
[0020]
一种孔隙参数与储层参数的地震岩石物理解析反演方法，具体包括如下步骤：
[0021]
步骤一，采用地震道集数据，通过叠前地震反演，获取观测弹性参数模型；
[0022]
步骤二，设置基本岩石物理参数，基于测井数据，反演得到井旁孔隙纵横比；
[0023]
步骤三，计算模型参数的概率分布函数，构建模型参数的初始模型；
[0024]
步骤四，利用地震岩石物理正演算子模拟弹性参数，计算模拟弹性参数与观测弹性参数的误差；
[0025]
步骤五，基于初始模型，计算精确正演算子对模型参数的偏导矩阵，构建线性正演算子；
[0026]
步骤六，基于线性正演算子，解析计算模型参数的后验期望值；
[0027]
步骤七，更新模型参数的初始模型为后验期望值，重复步骤四至步骤六，直至模拟与观测弹性参数的误差降低到预设范围内，迭代停止，输出模型参数的后验期望值，即为孔隙参数与储层参数的反演结果。
[0028]
作为本发明的一种优选技术方案，步骤一具体包括：
[0029]
获取工区pp波地震道集数据，通过叠前地震反演，构建目标函数
[0030][0031]
其中，d为地震道集数据，m为弹性参数，为弹性参数的平均值，σm为弹性参数的协方差矩阵，g为zoeppritz地震正演模型，σd为地震数据噪声的标准差，t为转置符号；弹性参数m包括纵波速度、横波速度和密度；
[0032]
采用贝叶斯线性反演方法，得到满足式(1)目标函数的最优解，获取纵波速度、横波速度和密度参数模型。
[0033]
作为本发明的一种优选技术方案：所述步骤二具体包括：
[0034]
设置基本岩石物理参数，包括岩石基质的体积模量、岩石基质的剪切模量、岩石基质的密度、烃类的体积模量、烃类的密度、卤水的体积模量和卤水的密度；
[0035]
基于测井观测的纵波速度与横波速度，通过匹配模型预测与井旁观测数据，反演得到井旁孔隙纵横比，构建目标函数
[0036][0037]
其中，α为孔隙纵横比，σ1与σ2为误差项权系数，和分别为模型预测和井旁观测的纵波速度，和分别为模型预测和井旁观测的横波速度，t为转置符号；
[0038]
采用粒子群优化算法，得到满足式(2)目标函数的最优解，获取井旁孔隙纵横比参数。
[0039]
作为本发明的一种优选技术方案：步骤三具体包括：
[0040]
通过测井观测数据获得井旁孔隙度与井旁含水饱和度，结合反演得到的井旁孔隙纵横比，构建模型参数的混合高斯模型，其概率分布函数为；
[0041][0042]
其中，z为模型参数，nk表示第k个高斯分布，和λk分别为nk的期望、协方差和权系数，c为高斯分布个数；模型参数包括孔隙纵横比、孔隙度与含水饱和度；
[0043]
通过期望最大算法计算混合高斯模型，得到c个高斯分布的期望、协方差和权系
数，再通过测井插值，得到模型参数的初始模型。
[0044]
作为本发明的一种优选技术方案：步骤四具体包括：
[0045]
利用解耦微分等效介质模型计算岩石骨架的弹性模量
[0046]
kd(φ)＝km(1-φ)
p
ꢀꢀ
(4)
[0047]
μd(φ)＝μm(1-φ)qꢀꢀ
(5)
[0048]
其中，φ为孔隙度，kd和μd分别为岩石骨架的体积模量和剪切模量，km和μm分别为岩石基质的体积模量和剪切模量，p和q为孔隙结构系数；
[0049]
基于岩石骨架的弹性模量，利用gassmann方程计算含流体岩石的弹性模量
[0050][0051]
μs＝μdꢀꢀ
(7)
[0052]
其中，ks和μs分别为含流体岩石的体积模量和剪切模量，kf为混合流体的体积模量，采用wood模型计算
[0053][0054]
其中，kw和kh分别为卤水和烃类的体积模量，sw为含水饱和度；
[0055]
根据含流体岩石的弹性模量，得到模拟弹性参数，
[0056][0057][0058]
ρ＝φρf+(1-φ)ρmꢀꢀ
(11)
[0059]
其中，v
p
、vs和ρ分别为纵波速度、横波速度和密度，ρf和ρm分为流体和岩石基质的密度；计算模拟弹性参数与观测弹性参数的误差ε
[0060]
ε＝||m
est-m
obs
||2ꢀꢀ
(12)
[0061]
其中，m
est
为模拟弹性参数，m
obs
为观测弹性参数即由步骤一获取的观测弹性参数模型。
[0062]
作为本发明的一种优选技术方案：步骤五具体包括：
[0063]
基于模型参数的初始模型，计算精确地震岩石物理正演算子对模型参数的偏导矩阵
[0064]
[0065]
其中，r为精确正演算子，通过联立公式(4)
–
(11)获得，z为模型参数，zo为模型参数的初始模型，g为精确正演算子的偏导矩阵，模型参数包括孔隙纵横比、孔隙度与含水饱和度，其具体计算方法如下：
[0066]
纵波速度对孔隙度的偏导数为：
[0067][0068]
其中，
[0069][0070][0071]
纵波速度对含水饱和度的偏导数为：
[0072][0073]
其中，
[0074][0075]
纵波速度对孔隙纵横比的偏导数为：
[0076][0077]
其中，
[0078][0079][0080]
其中，
[0081][0082][0083]
[0084]
横波速度对孔隙度的偏导数为：
[0085][0086]
横波速度对含水饱和度的偏导数为：
[0087][0088]
横波速度对孔隙纵横比的偏导数为：
[0089][0090]
其中，
[0091][0092]
其中，
[0093][0094][0095][0096][0097][0098][0099][0100][0101]
[0102][0103]
密度对孔隙度的偏导数为：
[0104][0105]
密度对含水饱和度的偏导数为：
[0106][0107]
密度对孔隙纵横比的偏导数为：
[0108][0109]
在公式(14)
–
(41)中，v
p
、vs和ρ分别为纵波速度、横波速度和密度，kd和μd分别为岩石骨架的体积模量和剪切模量，ks和μs分别为含流体岩石的体积模量和剪切模量，km和μm分别为岩石基质的体积模量和剪切模量，kf为混合流体的体积模量，p和q为孔隙结构系数，ρf为流体的密度，ρm为岩石基质的密度，ρw为卤水的密度，ρh为烃类的密度，φ为孔隙度，α为孔隙纵横比，sw为含水饱和度，f1、f2、f3、f4、f5、f6、f7、f8与f9为9个不同的孔隙因子系数，a和b为孔隙扁度系数，r为基质模量系数；
[0110]
基于精确正演算子的偏导矩阵，构建线性正演算子
[0111]
m＝gz+r(zo)-gzo+e
ꢀꢀ
(42)
[0112]
其中，m为弹性参数，e为随机误差。
[0113]
作为本发明的一种优选技术方案：步骤六具体包括：
[0114]
利用公式(42)线性正演算子，基于贝叶斯线性反演理论，解析计算每个高斯分布中模型参数的期望值
[0115][0116]
其中，g为精确正演算子的偏导矩阵，m为弹性参数，σe为随机误差的协方差，和为第k个高斯分布的期望和协方差，为第k个高斯分布的模型参数期望值，t为转置符号；
[0117]
计算模型参数的后验期望值
[0118][0119]
其中，nk表示第k个高斯分布，c为高斯分布个数，λk为nk的权系数，μ
post
为模型参数的后验期望值。
[0120]
作为本发明的一种优选技术方案：步骤七预设范围为观测弹性参数值的百分之一。
[0121]
有益效果：本技术所述孔隙参数与储层参数的地震岩石物理解析反演方法采用的以上技术方案与现有技术相比，具有以下技术效果：
[0122]
1、用于同步解析反演孔隙参数及储层参数，可以提高致密砂岩与碳酸盐岩类非常规储层的储层参数定量预测精度。
[0123]
2、本发明推导了精确正演算子对孔隙纵横比、孔隙度与含水饱和度的偏导矩阵，可以同步反演孔隙参数与储层参数；选用了线性反演策略以解析计算反演结果，提高了反演时效性及稳定性，相较于传统的地震岩石物理反演方法，本发明提高了复杂孔隙储层的储层参数反演精度。
附图说明：
[0124]
图1为本技术一种孔隙参数与储层参数的地震岩石物理解析反演方法的流程示意图；
[0125]
图2为本技术实施例中地震道集数据与观测弹性参数模型；其中a为实施例中的地震道集数据图，b为通过叠前地震反演获取的观测弹性参数图；
[0126]
图3为本技术实施例中观测储层参数与井旁孔隙纵横比；其中a为实施例中的观测储层参数图，b为反演得到的井旁孔隙纵横比图；
[0127]
图4为本技术实施例中利用常规方法的储层参数反演结果；
[0128]
图5为本技术实施例中利用本技术方法的孔隙参数与储层参数反演结果。
具体实施方式
[0129]
下面将结合本发明的附图，对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动条件下所获得的所有其它实施例，都属于本发明保护的范围。
[0130]
实施例1
[0131]
如图1所示，本实施例提供了一种孔隙参数与储层参数的地震岩石物理解析反演方法，具体包括以下步骤：
[0132]
一、采用地震道集数据，通过叠前地震反演，获取观测弹性参数模型及构建孔隙和储层参数初始模型：
[0133]
获取工区pp波地震道集数据，通过叠前地震反演，其目标函数如下
[0134][0135]
其中，d为地震道集数据，m为弹性参数，g为zoeppritz地震正演模型，σd为地震数据噪声的标准差，t为转置符号；和σm分别为弹性参数的平均值与协方差矩阵，可通过井旁弹性参数估算所得；弹性参数m包括纵波速度、横波速度和密度；
[0136]
采用贝叶斯线性反演方法，得到满足式(1)目标函数的最优解，获取纵波速度、横波速度和密度参数模型，即观测弹性参数模型；
[0137]
设置基本岩石物理参数，包括岩石基质的体积模量km、岩石基质的剪切模量μm、岩
石基质的密度ρm，烃类的体积模量kh、烃类的密度ρh、卤水的体积模量kw、卤水的密度ρw；实施例中基本岩石物理参数设置如下：km＝50gpa，μm＝23gpa，kh＝0.012gpa，kw＝2.5gpa，ρm＝2.62g/cm3，ρh＝0.02g/cm3，ρw＝1.10g/cm3；
[0138]
基于测井观测的纵波速度与横波速度，反演得到井旁孔隙纵横比，其目标函数如下
[0139][0140]
其中，α为孔隙纵横比，σ1与σ2为误差项权系数，和分别为模型预测和井旁观测的纵波速度，和分别为模型预测和井旁观测的横波速度，t为转置符号；
[0141]
采用粒子群优化算法，实施例中设置σ1＝0.6，σ1＝0.4，得到满足式(2)目标函数的最优解，获取井旁孔隙纵横比参数；
[0142]
通过测井观测数据获得井旁孔隙度与井旁含水饱和度，结合反演所得的井旁孔隙纵横比，构建混合高斯模型，其概率分布函数为
[0143][0144]
其中，z为模型参数，模型参数包括孔隙纵横比、孔隙度与含水饱和度，nk表示第k个高斯分布，和λk分别为nk的期望、协方差和权系数，c为高斯分布个数；
[0145]
实施例中设置c＝2，通过期望最大算法计算混合高斯模型，得到模型参数的协方差与均值，通过测井插值，得到模型参数的初始模型。
[0146]
二、计算精确正演算子对模型参数的偏导矩阵：
[0147]
利用解耦微分等效介质模型计算岩石骨架的弹性模量
[0148]
kd(φ)＝km(1-φ)
p
ꢀꢀ
(4)
[0149]
μd(φ)＝μm(1-φ)qꢀꢀ
(5)
[0150]
其中，φ为孔隙度，kd和μd分别为岩石骨架的体积模量和剪切模量，km和μm分别为岩石基质的体积模量和剪切模量，p和q为孔隙结构系数；
[0151]
基于岩石骨架的弹性模量，利用gassmann方程计算含流体岩石的弹性模量
[0152][0153]
μs＝μ
d (7)
[0154]
其中，ks和μs分别为含流体岩石的体积模量和剪切模量，kf为混合流体的体积模量，采用wood模型计算
[0155][0156]
其中，kw和kh分别为卤水和烃类的体积模量，sw为含水饱和度；
[0157]
基于含流体岩石的弹性模量，得到模拟弹性参数
[0158][0159][0160]
ρ＝φρf+(1-φ)ρmꢀꢀ
(11)
[0161]
其中，v
p
、vs和ρ分别为纵波速度、横波速度和密度，ρf和ρm分为流体和岩石基质的密度；
[0162]
计算模拟弹性参数与观测弹性参数的误差ε
[0163]
ε＝||m
est-m
obs
||2ꢀꢀ
(12)
[0164]
其中，m
est
为模拟弹性参数，m
obs
为观测弹性参数，即由步骤一获取的观测弹性参数模型；
[0165]
基于模型参数的初始模型，计算精确地震岩石物理正演算子对模型参数的偏导矩阵
[0166][0167]
其中，r为精确正演算子，通过联立公式(4)
–
(11)获得，z为模型参数，模型参数包括孔隙纵横比、孔隙度与含水饱和度，zo为模型参数的初始模型，g为精确正演算子的偏导矩阵，其具体计算方法如下：
[0168]
计算纵波速度对孔隙度的偏导数：
[0169][0170]
其中，
[0171][0172][0173]
计算纵波速度对含水饱和度的偏导数：
[0174][0175]
其中，
[0176][0177]
计算纵波速度对孔隙纵横比的偏导数：
[0178][0179]
其中，
[0180][0181][0182]
其中，
[0183][0184][0185][0186]
计算横波速度对孔隙度的偏导数：
[0187][0188]
计算横波速度对含水饱和度的偏导数：
[0189][0190]
计算横波速度对孔隙纵横比的偏导数：
[0191][0192]
其中，
[0193][0194]
其中，
[0195][0196][0197][0198][0199][0200][0201][0202][0203][0204][0205]
计算密度对孔隙度的偏导数：
[0206][0207]
计算密度对含水饱和度的偏导数：
[0208][0209]
计算密度对孔隙纵横比的偏导数：
[0210][0211]
在公式(14)
–
(41)中，v
p
、vs和ρ分别为纵波速度、横波速度和密度，kd和μd分别为岩
石骨架的体积模量和剪切模量，ks和μs分别为含流体岩石的体积模量和剪切模量，km和μm分别为岩石基质的体积模量和剪切模量，kf为混合流体的体积模量，p和q为孔隙结构系数，ρf为流体的密度，ρm为岩石基质的密度，ρw为卤水的密度，ρh为烃类的密度，φ为孔隙度，α为孔隙纵横比，sw为含水饱和度，f1、f2、f3、f4、f5、f6、f7、f8与f9为9个不同的孔隙因子系数，a和b为孔隙扁度系数，r为基质模量系数。
[0212]
基于公式(14)-(41)的计算结果，形成精确正演算子的偏导矩阵g，进而构建线性正演算子
[0213]
m＝gz+r(zo)-gzo+e
ꢀꢀ
(42)
[0214]
其中，m为弹性参数，e为随机误差。
[0215]
三、解析计算孔隙参数与储层参数的反演结果：
[0216]
利用公式(42)线性正演算子，基于贝叶斯线性反演理论，解析计算每个高斯分布中模型参数的期望值
[0217][0218]
其中，g为精确正演算子的偏导矩阵，m为弹性参数，σe为随机误差的协方差，和为第k个高斯分布的期望和协方差，为第k个高斯分布的模型参数期望值；
[0219]
计算模型参数的后验期望值
[0220][0221]
其中，nk表示第k个高斯分布，c为高斯分布个数，λk为nk的权系数，μ
post
为模型参数的后验期望值；
[0222]
更新模型参数的初始模型为后验期望值，重复上述步骤，直至模拟与观测弹性参数的误差降低到观测弹性参数值的百分之一，迭代停止，输出模型参数的后验期望值，即为孔隙参数与储层参数的反演结果。
[0223]
图2(a)为本实施例中的地震道集数据，图2(b)为通过叠前地震反演获取的观测弹性参数，包括纵波速度、横波速度和密度。图3(a)为本实施例中的观测储层参数，包括孔隙度与汗水饱和，图3(b)为反演得到的井旁孔隙纵横比。图4为本实施例中利用常规方法(即固定孔隙纵横比)的储层参数反演结果。从图4中可见，孔隙度与含水饱和度的反演结果(虚线)与观测数据(实线)均有较大偏差。图5为本实施例中利用本发明方法的孔隙参数与储层参数反演结果。从图5中可见，孔隙度与含水饱和度的反演结果(虚线)与观测数据(实线)基本一致，此外孔隙纵横比的反演结果(虚线)也与井旁数据(实线)基本吻合。经计算所得，图5中孔隙度与含水饱和度的反演结果与观测数据的平均误差分别为0.056与0.258，而图4中的平均误差分别为0.092与0.517，由此可见，本发明方法可以有效提高储层参数反演结果的准确度。
[0224]
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。