一种连通孔隙含量的地震岩石物理非线性反演方法

本发明属于非常规储层地震勘探，特别涉及一种连通孔隙含量的地震岩石物理非线性反演方法。

背景技术：

1、地震岩石物理反演结合了岩石物理理论与地震反演技术，旨在通过地震数据定量预测储层物性参数，是油气勘探及开发的重要技术手段。随着油气勘探的不断推进，常规储层日益枯竭，碳酸盐岩、致密砂岩等非常规储层已成为主要的勘探目标，然而此类储层普遍发育复杂孔隙结构，常规地震岩石物理反演方法未考虑储层孔隙结构的复杂性，难以定量预测与渗透性密切相关的孔隙连通性，限制了此类储层的地震勘探与开发。

2、常规地震岩石物理反演方法主要将储层物性参数(孔隙度、泥质含量、含水饱和度等)作为待预测变量，例如文献[1]-[4]公开的现有技术，但上述方法无法对储层孔隙结构信息进行定量预测与评价。研究发现，碳酸盐岩、致密砂岩等非常规储层通常发育复杂孔隙结构，其孔隙空间通常包含连通孔隙及孤立孔隙，相关孔隙参数对储层渗透性评价具有重要意义，对储层油气开发具有指导价值(详见文献[5]-[7])。对此，考虑孔隙结构的地震岩石物理反演技术日益受到关注，例如zhao等(2013)提出了一种碳酸盐岩孔隙类型的地震预测方法(详见文献[8])；li等(2022)提出了一种预测致密砂岩硬孔孔隙百分比的地震反演方法(详见文献[9])；李红兵等(2021)发展了一种预测复杂储层孔隙纵横比的地震反演方法(详见文献[10])。虽然上述反演技术能够定量预测孔隙类型或结构参数，但在岩石物理建模中未考虑孔隙连通情况，无法对孔隙连通性进行定量评估。此外，为定量预测储层连通孔隙含量，需建立联合连通孔隙及孤立孔隙的地震岩石物理模型，其正演算子非线性程度高，因此相应的反演问题多解性强、解析解通常难以获取，常规拟线性或确定性反演算法难以适用(详见文献[11]-[12])，需借助非线性优化算法求解反问题，并依赖可靠的先验信息约束反演过程，提高预测结果精度。

3、综上，针对发育复杂孔隙结构的碳酸盐岩、致密砂岩储层，需在岩石物理建模及反演过程中考虑孔隙连通情况，并发展有效的非线性反演算法提高反演结果精度，因此，开展连通孔隙含量的地震岩石物理非线性反演方法，对提高此类非常规储层的地震定量预测精度具有重要意义。

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

5、[1]fjeldstad，t.，avseth，p.，omre，h.，2021，a one-stepbayesian inversionframework for 3d reservoir characterization based on a gaussian mixturemodel-a norwegian sea demonstration[j].geophysics，86(2)，r221-r236.

6、[2]luo，c.，ba，j.，guo，q.，2023，probabilistic seismic petrophysicalinversion with statistical double-porosity biot-rayleigh model[j].geophysics，88(3)，m 157-m171.

7、[3]桂金咏，高建虎，雍学善，李胜军，2015，基于双相介质理论的储层参数反演方法[j].地球物理学报，9，3424-3438.

8、[4]李坤，印兴耀，宗兆云.2020.岩石物理驱动的相约束叠前地震概率化反演方法[j].中国科学：地球科学，50，124-146.

9、[5]邓继新，周浩，王欢等.2015.基于储层砂岩微观孔隙结构特征的弹性波频散响应分析[j].地球物理学报，58，3389-3400.

10、[6]贺锡雷，林凯，张祖豪等.2018.孔隙度和孔隙结构对储层特性影响的定量比较[j].石油物探，57(2)：179-185.

11、[7]lai，j.，wang，g.，wang，z.，et al.，2018，a review on pore structurecharacterization in tight sandstones[j].earth-science reviews，177，436-457.

12、[8]zhao，l.，nasser，m.，han，d.h.，2013，quantitative geophysical pore-typecharacterization and its geological implication in carbonate reservoirs[j].geophysical prospecting，61，827-841.

13、[9]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 multiporosity reservoirs[j].geophysics，87，m25-m41.

14、[10]李红兵，张佳佳，潘豪杰等，2021，基于弹性阻抗的孔隙结构与物性参数非线性同步反演[j].中国科学：地球科学，51，1166-1180.

15、[11]liu，q.，dong，n.，ji，y.，et al.，2018，direct reservoir propertyestimation based on prestack seismic inversion[j].journal ofpetroleum scienceand engineering，171，1475-1486.

16、[12]张佳佳，印兴耀，张广智等.2020.基于线性化岩石物理反演的物性参数预测方法[j].石油勘探与开发，47(1)，57-64。

技术实现思路

1、解决的技术问题：本发明提出了一种连通孔隙含量的地震岩石物理非线性反演方法，解决现有技术中存在的碳酸盐岩及致密砂岩储层普遍发育复杂孔隙结构，常规地震岩石物理反演方法未考虑储层孔隙结构的复杂性；仅将储层物性参数(孔隙度、泥质含量、含水饱和度等)作为待预测变量，无法对储层孔隙结构信息进行定量预测与评价；在岩石物理建模中未考虑孔隙连通情况，无法对孔隙连通性进行定量预测与评估；联立连通及孤立孔隙的地震岩石物理模型非线性程度高，反演问题多解性强、解析解通常难以获取，常规拟线性或确定性反演算法难以适用等技术问题，旨在提高碳酸盐岩及致密砂岩储层等非常规储层的地震定量预测精度。

2、技术方案：

3、本发明公开了一种连通孔隙含量的地震岩石物理非线性反演方法，所述地震岩石物理非线性反演方法包括以下步骤：

4、步骤a，计算岩石骨架、孔隙流体及含流体岩石弹性参数，构建地震岩石物理模型；

5、步骤b，基于测井数据和地震岩石物理模型，获取井旁总孔隙度及连通孔隙含量；

6、步骤c，根据井旁总孔隙度及连通孔隙含量，构建目标参数的蒙特卡洛模拟样本集；

7、步骤d，基于目标参数的蒙特卡洛模拟样本集，估算先验概率密度函数，构建目标参数的后验概率密度函数；

8、步骤e，对目标参数进行迭代更新；其中，在每次迭代过程中，依据目标参数的概率密度函数，随机生成目标参数的更新模型，并计算模型接收概率，根据模型接收概率判别是否接收更新模型；

9、步骤f，重复步骤e，直至重复次数达到最大迭代次数，输出目标参数模型，计算连通孔隙含量模型的期望均值，作为最终反演结果。

10、作为本发明的一种优选技术方案，步骤a中，计算岩石骨架、孔隙流体及含流体岩石弹性参数，构建地震岩石物理模型的过程包括以下子步骤：

11、步骤a1，设置基本岩石物理参数，基本岩石物理参数包括矿物的体积模量、剪切模量、密度，流体的体积模量和密度；基于岩心x射线光谱分析，确定岩石矿物百分含量，采用v-r-h模型，计算岩石基质的体积模量、剪切模量及密度：

12、

13、

14、

15、式中，km、gm和ρm分别为岩石基质的体积模量、剪切模量和密度，kmi、gmi、ρmi和vmi分别为第i种矿物的体积模量、剪切模量、密度和该矿物的百分含量，n为矿物种类总数；

16、步骤a2，将岩石孔隙划分为连通孔隙与孤立孔隙，采用自洽模型计算岩石骨架的等效弹性模量：

17、(1-φ)(km-k*)p*m+φvcn(kcn-k*)p*cn+φvds(kds-k*)p*ds＝0          (4)；

18、(1-φ)(gm-g*)q*m+φvcn(gcn-g*)q*cn+φvds(gds-g*)q*ds＝0          (5)；

19、式中，k*和g*分别为岩石骨架的等效体积模量和剪切模量；φ为井旁总孔隙度；vcn和vds分别为连通孔隙含量和孤立孔隙含量，满足vcn+vds＝1；kcn和gcn分别为连通孔隙包裹体的体积模量和剪切模量，kds和gds分别为孤立孔隙包裹体的体积模量和剪切模量；p*cn和q*cn分别为连通孔隙包裹体的几何因子，p*ds和q*ds分别为孤立孔隙包裹体的几何因子，p*m和q*m分别为岩石基质的几何因子；

20、步骤a3，设置连通孔隙为空孔隙，孤立孔隙为含流体孔隙，则：

21、kcn＝gcn＝gds＝0                              (6)；

22、kds＝kfl                                  (7)；

23、式中，kfl为流体的体积模量；

24、步骤a4，基于含孤立孔隙岩石骨架，针对连通孔隙，采用biot-gassmann模型计算储层岩石的体积模量和剪切模量：

25、

26、grk＝g*  (9)；

27、式中，krk和grk分别为储层岩石的体积模量和剪切模量，ksk为含孤立孔隙岩石骨架的体积模量：

28、

29、步骤a5，储层岩石的纵波速度、横波速度和密度分别表示为：

30、

31、

32、ρrk＝(1-φ)ρm+φρfl  (13)；

33、式中，vp、vs和ρrk分别为储层岩石的纵波速度、横波速度和密度，ρfl为流体的密度。

34、作为本发明的一种优选技术方案，步骤b中，基于测井数据，获取井旁总孔隙度及连通孔隙含量的过程包括以下子步骤：

35、步骤b1，基于测井观测的纵波阻抗及横波阻抗数据，通过匹配模型预测数据与测井观测数据，反演得到井旁总孔隙度及连通孔隙含量，构建误差函数：

36、

37、其中，

38、

39、

40、式中，vj＝[φj，vcnj]为目标参数，φj为第j个测井采样点的井旁总孔隙度，vcnj为第j个测井采样点的连通孔隙含量；和分别为第j个测井采样点的观测纵波阻抗和横波阻抗，和分别为第j个测井采样点的预测纵波阻抗和横波阻抗，vp、vs和ρrk分别为储层岩石的纵波速度、横波速度和密度，krk和grk分别为储层岩石的体积模量和剪切模量，j＝1，2，...，ns，ns为测井采样点总数，κ为误差函数的权系数；

41、步骤b2，为平衡纵波阻抗和横波阻抗的量级差异，通过以下方式估算误差函数的权系数：

42、

43、步骤b3，采用快速模拟退火算法，估算使误差函数极小的目标参数，获取第j个测井采样点的井旁总孔隙度及连通孔隙含量vj＝[φj，vcnj]。

44、作为本发明的一种优选技术方案，步骤c中，根据井旁总孔隙度及连通孔隙含量，构建目标参数样本集的过程包括以下子步骤：

45、步骤c1，基于井旁总孔隙度及连通孔隙含量，构建目标参数的高斯混合模型，采用期望最大化算法，估算高斯混合模型参数，则目标参数的概率分布函数为：

46、

47、式中，v＝[φ，vcn]为目标参数，φ为井旁总孔隙度，vcn为连通孔隙含量，nk表示第k个高斯分布，μv|k、∑v|k和σk分别为第k个高斯分布的期望、协方差和权系数，nc为高斯分布总数；

48、步骤c2，根据目标参数的概率分布函数，采用蒙特卡洛随机模拟，获取增广目标参数，将增广目标参数与井旁总孔隙度及连通孔隙含量合并，构建目标参数的蒙特卡洛模拟样本集vmci*＝[φi*，vcni*]，i*＝1，2，3，...，nmc，其中nmc为蒙特卡洛样本总数。

49、作为本发明的一种优选技术方案，步骤d中，基于目标参数样本集，估算先验概率密度函数，构建后验概率密度函数的过程包括以下子步骤：

50、步骤d1，基于目标参数样本集，采用epanechnikov核密度估计法，计算目标参数的先验概率密度函数：

51、

52、其中，

53、

54、式中，p(v)为目标参数的先验概率密度函数，v＝[φ，vcn]为目标参数，φ为井旁总孔隙度，vcn为连通孔隙含量，nmc为蒙特卡洛样本总数，l为目标参数的维度，f(x)为epanechnikov核函数，hv为核函数的缩放因子；

55、步骤d2，基于贝叶斯框架，目标参数的后验概率密度函数等价于：

56、

57、式中，z为观测阻抗数据，包括观测纵波阻抗和横波阻抗，由叠前地震反演获取；p(v|z)和p(z|v)分别为目标参数的后验概率密度函数和似然函数；p(z)为观测阻抗数据的边际概率，是一个常数；

58、步骤d3，假设观测阻抗数据误差满足期望为零的高斯分布，则目标参数的似然函数表示为：

59、

60、式中，d为观测阻抗数据的维度，联合式(19)、式(21)和式(22)，目标参数的后验概率密度函数表示为：

61、

62、式中，zpdt为预测阻抗数据，包括预测纵波阻抗和预测横波阻抗，由步骤a所构建的岩石物理模型预测所得；cnor为归一化常数，∑z为观测阻抗数据误差的协方差矩阵。

63、作为本发明的一种优选技术方案，步骤e中，对目标参数进行迭代更新的过程包括以下子步骤：

64、依据目标参数的先验概率密度函数，随机生成目标参数的更新模型，基于metropolis准则计算模型接收概率，第k次迭代时模型接收概率表示为：

65、

66、其中，v*为目标参数的更新模型，v(k)为第k次迭代更新前的目标参数，pacp为模型接收概率，p(v*|z)和p(v(k)|z)分别为目标参数的更新模型的后验概率密度函数和第k次迭代更新前的目标参数的后验概率密度函数；z为观测阻抗数据，包括观测纵波阻抗和横波阻抗；

67、生成介于0至1之间的随机数，若随机数大于等于pacp则接收目标参数的更新模型，令第k+1次迭代更新前的目标参数v(k+1)＝v*，否则目标参数模型保持不变，令第k+1次迭代更新前的目标参数v(k+1)＝v(k)。

68、作为本发明的一种优选技术方案，步骤f中，最大迭代次数kmax为大于等于20000的正整数。

69、作为本发明的一种优选技术方案，步骤f中，输出目标参数模型，计算连通孔隙含量模型的期望均值，作为最终反演结果的过程包括以下子步骤：

70、选择最后m次迭代时的目标参数模型集合vest＝[vw，vw+1，vw+2，...，vw+m-1]，vw为第w次迭代更新前的目标参数，w＝kmax-m-1，依据目标参数模型集合vest估算连通孔隙含量模型的期望均值，作为最终反演结果。

71、有益效果：

72、第一，本发明的连通孔隙含量的地震岩石物理非线性反演方法，用于反演复杂孔隙储层的连通孔隙含量，反演结果可对储层渗透率进行评价，可以提高致密砂岩与碳酸盐岩类非常规储层的地震定量预测精度。

73、第二，本发明的连通孔隙含量的地震岩石物理非线性反演方法，构建了考虑连通孔隙与孤立孔隙的地震岩石物理模型，可以用于反演连通孔隙含量；发展了地震非线性反演方法，结合目标参数先验概率密度函数稳定反演过程，相较于传统的地震岩石物理反演方法，本发明提高了复杂孔隙储层的地震定量预测精度。