CN114185090B - 岩相和弹性参数同步反演方法、装置、电子设备及介质 - Google Patents
岩相和弹性参数同步反演方法、装置、电子设备及介质 Download PDFInfo
- Publication number
- CN114185090B CN114185090B CN202010968412.XA CN202010968412A CN114185090B CN 114185090 B CN114185090 B CN 114185090B CN 202010968412 A CN202010968412 A CN 202010968412A CN 114185090 B CN114185090 B CN 114185090B
- Authority
- CN
- China
- Prior art keywords
- lithofacies
- inversion
- elastic
- lithology
- elastic parameter
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 54
- 230000001360 synchronised effect Effects 0.000 title claims abstract description 20
- 230000006870 function Effects 0.000 claims abstract description 49
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 26
- 239000011159 matrix material Substances 0.000 claims description 57
- 238000005070 sampling Methods 0.000 claims description 57
- 239000013598 vector Substances 0.000 claims description 28
- 239000011435 rock Substances 0.000 claims description 25
- 238000009826 distribution Methods 0.000 claims description 20
- 208000035126 Facies Diseases 0.000 claims description 18
- 230000008878 coupling Effects 0.000 claims description 10
- 238000010168 coupling process Methods 0.000 claims description 10
- 238000005859 coupling reaction Methods 0.000 claims description 10
- 230000007704 transition Effects 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 7
- 238000004590 computer program Methods 0.000 claims description 4
- 238000005056 compaction Methods 0.000 abstract description 5
- 239000004576 sand Substances 0.000 description 6
- 239000012530 fluid Substances 0.000 description 5
- 230000008569 process Effects 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 238000010276 construction Methods 0.000 description 3
- 238000004062 sedimentation Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 230000008021 deposition Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000005137 deposition process Methods 0.000 description 1
- 235000021185 dessert Nutrition 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 238000002948 stochastic simulation Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
- G01V2210/6242—Elastic parameters, e.g. Young, Lamé or Poisson
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/66—Subsurface modeling
Abstract
本发明提出了一种岩相和弹性参数同步反演方法、装置、电子设备及介质。反演方法包括：初始化待反演工区的岩相；利用岩相约束的贝叶斯线性反演计算弹性参数，并计算弹性参数关于岩相的似然函数；利用修正的Viterbi算法计算岩相关于弹性参数的后验概率最大值，并更新岩相，反复迭代直至收敛。本发明基于测井数据建立不同岩相的弹性低频趋势，可有效避免后续岩相预测的弹性压实“陷阱”。通过利用岩相反演结果约束弹性参数反演结果的低频和弹性参数之间的相关性，提高弹性参数反演结果的精度，进而更新弹性参数关于岩相的似然函数，并在贝叶斯网络和马尔科夫随机场先验约束下，更新岩相反演结果，提高岩相反演结果的稳定性和空间连续性。
Description
技术领域
本发明属于油气探勘领域，涉及地震资料解释与反演技术，更具体地涉及一种基于马尔科夫模型的岩相和弹性参数同步反演方法、反演装置、电子设备及存储介质。
背景技术
基于地震数据的岩相预测一直是油气地球物理在勘探开发应用中的重点和难点，通过结合测井数据的岩相分类和油气解释结论，可以实现基于地震数据的岩性、物性以及含油气性的定量分类，因此也可以更加直观的给出勘探开发的有利“甜点区域”。
现阶段基于地震数据的岩相预测大致分为两类：一类是基于监督或者无监督学习的属性融合方法，该方法主要是依靠地震数据或地震属性的空间结构的相似性，并借助机器学习对其挖掘，实现岩相分类。这种分类方法并不需要借助明确的地球物理意义，与其说是岩相预测，不如说是大套的地震相预测，因此其垂向分辨率可能较低。例如，波形聚类、PCA属性融合等方法属于此类的岩相预测。第二类是基于叠前弹性参数反演结果的基础上，并结合测井上获取的岩石物理关系，利用统计的方法给出岩相关于弹性参数的似然关系，进而对岩相进行分类预测，由于该方法需要依靠岩石物理关系，因此具有相对明确的地球物理意义。不难看出，岩相预测精度取决于弹性参数反演结果的精度，而弹性参数反演的精度又取决于弹性参数低频趋势构建的精度。现阶段，弹性参数低频趋势的构建仅仅是借助于测井曲线在层位控制下进行横向插值，其精度受限于参与井的多少、插值方式以及沉积关系等因素，另外这种低频趋势的构建模式，也不能精确地表征岩性压实趋势。事实上，不同岩相的弹性参数低频趋势是不一样的，弹性参数低频趋势的精度又和岩相的比例有着密切的关系。因此，岩相反演问题和弹性参数反演问题是一个耦合的关系，而现阶段的岩相反演方法大都采用分步预测，即先反演弹性参数，再基于岩石物理关系，利用贝叶斯准则反演岩相，这种方法在一定程度会降低岩相预测的精度。
另外，在有地质意义的背景下，相邻采样点之间的岩相具有一定耦合性和相关性，例如在沉积过程中，t时刻沉积的岩相必然和t时刻之前沉积的岩相有关，而同一沉积时刻的相邻岩相之间也是存在明显的相关性，即垂向和横向都具有明显的相关性，但是现阶段工业界所用的基于贝叶斯准则的岩相反演方法是逐道逐个采样点计算的，这种反演方法从本质上忽略了岩相在空间的相关性和耦合性，当地震数据信噪比较低时，岩相反演的稳定性较差，产生明显的“破碎”现象，进而增加了储层预测和流体识别的不确定性。
Eidsvik等(2004)利用一阶马尔科夫随机场描述岩性和流体类别的先验分布，并结合Gassmann方程以及Shuey近似公式建立似然函数，通过贝叶斯推断整合了层位数据、测井数据以及AVO属性，建立了岩性和流体类别的后验分布，在此基础上，利用MCMC随机模拟抽取了满足后验分布的大量样本，最后基于最大后验期望准则来预测岩性和流体类别的分布，但是该方法仅局限于平面预测。Larsen(2006)利用HMM模型来建立1D情况下地震岩相的空间相关性，并结合贝叶斯原理整合叠前地震数据与岩石物理关系，建立了地震岩相条件于地震数据的近似后验概率，最后利用“前向-后向递归算法”从叠前地震数据中模拟地震岩相，然而该方法仅适用1维地震道预测；Ulvmoen(2010)在Larson(2006)和Buland(2008)工作的基础上，进一步引入马尔科夫随机场描述地震岩相在横向的连续性，在此基础上，通过对给出近似的岩相后验概率分布，利用前后-后向递归算法结合Gibbs抽样从叠前地震数据中直接反演岩相，但是该方法需要大量的随机模拟才能获取稳定的岩相体，其计算效率较低。
针对以上问题，本领域急需一种岩相和弹性参数同步反演方法，进一步提高反演的计算效率。
发明内容
针对常规弹性参数反演结果精度不高，以及常规岩相反演方法易受弹性压实趋势的影响，且未考虑岩相在空间的耦合性和连续性造成岩相反演结果稳定性差等原因。本发明提出一种岩相和弹性参数的同步反演方法，在整个过程中，弹性参数和岩相是交替更新的。通过利用岩相反演结果约束弹性参数反演结果的低频和弹性参数之间的相关性，提高弹性参数反演结果的精度，进而更新弹性参数关于岩相的似然函数，并在贝叶斯网络和马尔科夫随机场先验约束下，更新岩相反演结果，提高岩相反演结果的稳定性和空间连续性。
根据本发明的一个方面，提供一种基于马尔科夫模型的岩相和弹性参数同步反演方法，包括：
初始化待反演工区的岩相；
利用岩相约束的贝叶斯线性反演计算弹性参数，并计算弹性参数关于岩相的似然函数；
利用修正的Viterbi算法计算岩相关于弹性参数的后验概率最大值，并更新岩相，反复迭代直至收敛。
进一步地，所述初始化待反演工区的岩相包括：结合测井数据中的岩相曲线，统计岩相的先验分布和岩相概率转移矩阵P，将岩相先验分布中概率最大的岩相赋值到工区中每一道的每一个采样点中，作为初始化岩相数据体。
进一步地，利用岩相约束的贝叶斯线性反演计算弹性参数包括：
结合测井数据中的弹性参数，利用指数函数从测井数据中拟合所有岩相的纵波阻抗低频趋势kIT、纵横波速度比低频趋势kγT、密度低频趋势kρT，并计算所有岩相的弹性参数协方差矩阵kΣm；
输入M个角度叠加地震数据和M个角度子波，并构建角度叠加地震数据向量d和角度子波褶积矩阵；
利用岩相约束的贝叶斯线性反演计算本次迭代的弹性参数反演结果。
进一步地，所述弹性参数协方差矩阵kΣm为：
kΣm表示岩相标签为k对应的弹性参数的协方差矩阵，kIT表示岩相标签为k(k＝1，2，…K，K表示岩相的总个数)的纵波阻抗的低频趋势；kγT表示岩相标签为k的纵横波速度比的低频趋势；kρT表示岩相标签为k的密度的低频趋势；kσI、kσγ、kσρ、kcIγ、kcIρ和kcγρ分别表示岩相标签为k的纵波阻抗标准差、纵横波速度比标准差、密度标准差、纵波阻抗和纵横波速度比相关系数、纵波阻抗和密度相关系数以及纵横波速度比和密度的相关系数。
进一步地，岩相约束的贝叶斯线性反演获取的弹性参数反演结果为：
其中， πI为纵波阻抗自然对数反演结果，πγ为纵横波速度比自然对数反演结果，πρ为密度自然对数反演结果；/> πIT、πγT和πρT分别表示利用前一次迭代更新的岩相π确定的纵波阻抗自然对数低频模型、纵横波速度比自然对数低频模型和密度自然对数低频模型；πΣm表示利用前一次迭代更新的岩相π确定的参数协方差矩阵；G表示由子波褶积矩阵、AVO近似公式系数以及差分矩阵确定的AVO正演矩阵；/>表示噪声的方差，I表示单位矩阵。
进一步地，弹性参数关于岩相的似然函数采用三变量高斯分布描述：
其中，表示第t个采样点的弹性参数反演结果关于第t个采样点的岩相的似然函数，/>表示第t个采样点的弹性参数反演结果向量，/>表示第t个采样点弹性参数低频趋势向量，其中tπ＝1,2,L K。
进一步地，所述修正的Viterbi算法包括：
(1)建立三维工区内所有道的随机访问顺序；
(2)按照访问顺序选择一道，根据周边相邻四道的岩相信息以及当前道的似然函数信息，利用常规Viterbi算法计算当前道岩相关于弹性参数反演结果后验概率的最大值，再进行回溯更新当前道的岩相；
(3)重复(2)更新工区内所有道的岩相；
(4)重复(1)至(3)，进行多次迭代，直至工区的岩相比例稳定。
进一步地，当前道岩相关于弹性参数反演结果的后验概率为：
其中表示当前道的后验概率，xπ表示当前道的岩相向量，/>表示当前道的弹性参数反演结果向量，δxπ表示与当前道相邻的四道岩相向量；p(x,tπ|x,t+1π)表示当前道的第t个采样点的岩相x,tπ关于当前道的第t+1个采样点的岩相x,t+1π的条件概率，为概率转移矩阵P的第t+1行、第t列元素；V(x,tπ,δx,tπ)表示势函数，表征同一采样点的相邻岩相之间的连续性，其中δx,tπ表示与当前道相邻四道的第t个采样点的岩相，β表示耦合系数，I(x,tπ＝i,tπ)表示指示函数；表示当前道的第t个采样点的弹性参数反演结果/>关于当前道的第t个采样点的岩相的似然函数。
根据本发明的另一个方面，提供一种基于马尔科夫模型的岩相和弹性参数同步反演装置，包括：
初始化模块，初始化待反演工区的岩相；
计算模块，利用岩相约束的贝叶斯线性反演计算弹性参数，并计算弹性参数关于岩相的似然函数；
更新迭代模块，利用修正的Viterbi算法计算岩相关于弹性参数的后验概率最大值，并更新岩相，反复迭代直至收敛。
根据本发明的另一个方面，提供一种电子设备，包括：
存储器，存储有可执行指令；
处理器，所述处理器运行所述存储器中的所述可执行指令，以实现所述的基于马尔科夫模型的岩相和弹性参数同步反演方法。
根据本发明的另一个方面，提供一种计算机可读存储介质，该计算机可读存储介质存储有计算机程序，该计算机程序被处理器执行时实现所述的基于马尔科夫模型的岩相和弹性参数同步反演方法。
本发明扩充了常规弹性参数反演和岩相反演，基于弹性参数和岩相交替更新的策略，实现了岩相和弹性参数的同步反演。在此过程中，通过引入岩相约束的贝叶斯线性反演，提高了弹性参数的反演精度，通过引入贝叶斯网络和马尔科夫随机场先验约束以及修正Viterbi算法，提高了岩相反演的稳定性和空间连续性。
附图说明
通过结合附图对本公开示例性实施方式进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开示例性实施方式中，相同的参考标号通常代表相同部件。
图1一种基于马尔科夫模型的岩相和弹性参数同步反演方法，其特征在于，包括：
图2为根据本发明实施例的基于马尔科夫模型的岩相和弹性参数同步反演流程图。
图3为本发明实施例所用的测井数据。其中“*”代表真实的测井弹性参数，“---”代表基于不同岩相的弹性参数拟合的低频趋势，本申请中用不同的颜色命名了图中不同灰度的部分，具体地红色代表气砂，蓝色代表水砂，绿色代表致密砂，黑色代表泥岩。
图4为本发明实施例所使用的角度叠加地震数据，其中图4(a)为小角度叠加地震数据(反射角为11度)、图4(b)为中角度叠加地震数据(反射角为14度)、图4(c)为大角度叠加地震数据(反射角为17度)。
图5为利用本发明实施例方法获取的反演结果，其中图5(a)为岩相反演结果，井A为验证井，可以看出井旁道岩相反演结果与井A的岩相吻合度较高，4套含气砂岩均被预测出来，且岩相剖面的横向连续性非常好；图5(b)、图5(c)和图5(d)分别为纵波阻抗反演结果、纵横波速度比反演结果以及密度反演结果，这三种弹性参数井旁道反演结果与井A弹性参数吻合度较好，说明了弹性参数反演结果的精度较高。
具体实施方式
下面将参照附图更详细地描述本公开的优选实施方式。虽然附图中显示了本公开的优选实施方式，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。
本发明涉及一种岩相和弹性参数同步反演方法。弹性参数和岩相反演对地球物理储层预测和流体识别具有重要的意义，然而常规弹性参数反演方法是基于工区统一岩相这一假设，因此弹性参数反演结果精度不高；常规岩相反演方法未考虑弹性压实“陷阱”以及岩相在空间的耦合性，因此岩相反演结果稳定性差。
本发明提出一种岩相和弹性参数同步反演方法，在该方法中弹性参数与岩相交替迭代更新，并使用贝叶斯网络和马尔科夫随机场描述岩相的先验分布，通过利用修正的Viterbi算法替代传统的随机模拟，进一步提高反演的计算效率。
具体地，通过结合测井数据，拟合不同岩相的弹性参数低频趋势和建立不同岩相的弹性参数协方差矩阵；在此基础上，利用迭代的方式交替更新弹性参数和岩相，即利用前一次迭代的岩相数据体约束弹性参数反演结果的低频趋势和弹性参数之间的相关性，并利用更新的弹性参数反演结果构建新的似然函数，在贝叶斯网络和马尔科夫随机场先验约束下，进一步更新岩相，直至弹性参数和岩相稳定收敛。最终可以获取高精度的弹性参数反演结果和具有明确地质意义的岩相反演结果。
本发明提供一种基于马尔科夫模型的岩相和弹性参数同步反演方法，包括：
初始化待反演工区的岩相；
利用岩相约束的贝叶斯线性反演计算弹性参数，并计算弹性参数关于岩相的似然函数；
利用修正的Viterbi算法计算岩相关于弹性参数的后验概率最大值，并更新岩相，反复迭代直至收敛。
具体地，初始化待反演工区的岩相步骤可以包括：首先结合测井数据中的岩相曲线，统计岩相的先验分布和岩相概率转移矩阵P；接下来将岩相先验分布中概率最大的岩相赋值到工区中每一道的每一个采样点中，作为初始化岩相数据体。
进一步地，利用岩相约束的贝叶斯线性反演计算弹性参数包括：
结合测井数据中的弹性参数，利用指数函数从测井数据中拟合所有岩相的纵波阻抗低频趋势kIT、纵横波速度比低频趋势kγT、密度低频趋势kρT，并计算所有岩相的弹性参数协方差矩阵kΣm；
输入M个角度叠加地震数据和M个角度子波，并构建角度叠加地震数据向量d和角度子波褶积矩阵；
利用岩相约束的贝叶斯线性反演计算本次迭代的弹性参数反演结果。
进一步地，所述弹性参数协方差矩阵kΣm为：
kΣm表示岩相标签为k对应的弹性参数的协方差矩阵，kIT表示岩相标签为k(k＝1，2，…K，K表示岩相的总个数)的纵波阻抗的低频趋势；kγT表示岩相标签为k的纵横波速度比的低频趋势；kρT表示岩相标签为k的密度的低频趋势；kσI、kσγ、kσρ、kcIγ、kcIρ和kcγρ分别表示岩相标签为k的纵波阻抗标准差、纵横波速度比标准差、密度标准差、纵波阻抗和纵横波速度比相关系数、纵波阻抗和密度相关系数以及纵横波速度比和密度的相关系数。
具体地，根据前一次迭代的岩相数据体、计算获取的所有岩相的弹性参数低频趋势和协方差矩阵以及输入的角度叠加数据和角度子波褶积矩阵，利用岩相约束的贝叶斯线性反演计算本次迭代的弹性参数反演结果。
岩相约束的贝叶斯线性反演，是指利用给定的岩相数据体，确定贝叶斯线性AVO反演所需要的弹性参数低频趋势和模型参数的协方差矩阵。在整个技术流程中岩相和弹性参数是不断更新迭代，通过反演的岩相更新弹性参数的低频趋势和协方差矩阵，再由反演的弹性参数更新岩相，经过多次迭代直至收敛。
根据高斯分布的性质，可以得到弹性参数关于地震数据和岩相的后验概率期望，也就是岩相约束的贝叶斯线性反演获取的弹性参数反演结果：
其中 πI为纵波阻抗自然对数反演结果，πγ为纵横波速度比自然对数反演结果，πρ为密度自然对数反演结果；/> πIT、πγT和πρT分别表示利用前一次迭代更新的岩相π确定的纵波阻抗自然对数低频模型、纵横波速度比自然对数低频模型和密度自然对数低频模型，以πIT为例，假如第t个采样点的岩相标签为k，那么πIT的第t个元素为πIT(t)＝kIT(t)，πγT和πρT的表达方式类似；πΣm表示利用前一次迭代更新的岩相π确定的参数协方差矩阵，为3T行3T列稀疏矩阵(T为一道的采样点个数)，假如第t个采样点的岩相标签为k，那么矩阵元素/> πΣm(t,t+3T)＝kcIγ kσI kσγ、πΣm(t,t+6T)＝kcIρ kσI kσρ、πΣm(t+3T,t+6T)＝kcγρ kσγ kσρ、πΣm(t+3T,t)＝kcIγ kσI kσγ、πΣm(t+6T,t)＝kcIρ kσI kσρ和πΣm(t+6T,t+3T)＝kcγρ kσγ kσρ；G表示由子波褶积矩阵、AVO近似公式系数以及差分矩阵确定的AVO正演矩阵；/>表示噪声的方差，I表示单位矩阵。
进一步地，根据弹性参数反演结果、所有岩相的弹性参数低频趋势以及协方差矩阵，计算工区所有道、所有采样点弹性参数反演结果关于岩相的似然函数。弹性参数关于岩相的似然函数采用三变量高斯分布描述：
其中，表示第t个采样点的弹性参数反演结果关于第t个采样点的岩相的似然函数，/>表示第t个采样点的弹性参数反演结果向量，/>表示第t个采样点弹性参数低频趋势向量，其中tπ＝1,2,L K。
接下来，利用修正的Viterbi算法，计算岩相关于弹性参数反演结果的最大后验概率，进而获取本次迭代的岩相数据体。
为了使得岩相的反演结果更具有地质意义，使用贝叶斯网络描述岩相在垂向上的沉积关系，使用马尔科夫随机场描述岩相在横向上的连续性，即空间中某一点的岩相只跟它周边相邻的五个点的岩相相关，分别是该点同一道的下一点(符合沉积过程)、该点相邻的前后左右道中并与该点属于同一沉积时间的四个采样点(沿构造走向、且相邻的四个采样点)。
岩相在垂向上属于贝叶斯网络且在横向上属于马尔科夫随机场，这一性质决定无法直接使用传统Viterbi算法计算最大概率解，因此提出修正的Viterbi算法：
(1)建立三维工区内所有道的随机访问顺序；
(2)按照访问顺序选择一道，根据周边相邻四道的岩相信息以及当前道的似然函数信息，利用常规Viterbi算法更新当前道的岩相；
当前道岩相关于弹性参数反演结果的后验概率为：
其中表示当前道的后验概率，xπ表示当前道的岩相向量，/>表示当前道的弹性参数反演结果向量，δxπ表示与当前道相邻的四道岩相向量；p(x,tπ|x,t+1π)表示当前道的第t个采样点的岩相x,tπ关于当前道的第t+1个采样点的岩相x,t+1π的条件概率，为概率转移矩阵P的第t+1行、第t列元素；V(x,tπ,δx,tπ)表示势函数，表征同一采样点的相邻岩相之间的连续性，其中δx,tπ表示与当前道相邻四道的第t个采样点的岩相，β表示耦合系数，I(x,tπ＝i,tπ)表示指示函数，即如果当前道的第t个采样点的岩相与第i道(其中i属于当前道的邻域)的第t个采样点的相同，那么为1，否则为0，；/>表示当前道的第t个采样点的弹性参数反演结果/>关于当前道的第t个采样点的岩相的似然函数；当计算当前道岩相xπ时，其相邻道岩相δxπ是已知的，因此可以利用常规Viterbi算法计算后验概率/>的最大值，再进行回溯更新当前道的岩相；
(3)重复(2)更新工区内所有道的岩相；
(4)重复(1)至(3)，进行多次迭代，直至工区的岩相比例稳定。
重复进行反演计算弹性参数，计算弹性参数关于岩相的似然函数，更新岩相，对弹性参数与岩相进行迭代更新，直至两者的反演结果稳定收敛。
本发明基于测井数据建立不同岩相的弹性低频趋势，可有效避免后续岩相预测的弹性压实“陷阱”；在此基础上，对弹性参数和岩相进行交替更新。通过利用岩相反演结果约束弹性参数反演结果的低频和弹性参数之间的相关性，提高弹性参数反演结果的精度，进而更新弹性参数关于岩相的似然函数，并在贝叶斯网络和马尔科夫随机场先验约束下，更新岩相反演结果，提高岩相反演结果的稳定性和空间连续性。
为便于理解本发明实施例的方案及其效果，以下给出具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。
实施例1
下面结合图2-5和实例对本发明进行详细的描述。
本实施例是以角度部分叠加地震数据体和测井数据作为输入数据。前期经过地震数据采集、静校正、去噪、反褶积、速度分析、去除多次波等处理流程后，通过真振幅叠前时间偏移获取共反射点道集，将其转化为角道集，根据实际工区情况，对入射角范围进行划分，通过角度部分叠加获取角度部分叠加数据M个，其中M≥3。在反演前还需要进行井震标定，并提取角度子波。在本实施例方法中，岩相约束的贝叶斯线性反演是逐道计算的，弹性参数关于岩相的似然函数是逐道、逐点计算的，基于修正的Viterbi算法反演岩相则是按照随机顺序访问工区的每个CDP点，再进行整道计算的。
图2给出基于马尔科夫模型的岩相和弹性参数同步反演技术流程，如图所示主要包括：根据测井数据，统计岩相的先验分布和概率转移矩阵；根据测井数据，拟合所有岩相的弹性参数低频趋势，并计算所有岩相的弹性参数协方差矩阵；输入角度部分叠加地震数据和角度子波，并构建角度叠加地震数据向量和角度子波褶积矩阵；根据步骤一获取的岩相先验分布的最大值对应的岩相作为岩相数据体的初始值；根据前一次迭代的岩相数据体，确定弹性参数反演所需的低频模型和弹性参数协方差矩阵，利用岩相约束的贝叶斯线性反演计算本次迭代的弹性参数反演结果；根据弹性参数反演结果，计算弹性参数关于岩相的似然函数；利用修正的Viterbi算法计算岩相关于弹性参数反演结果的最大后验概率，进而获取本次迭代的岩相数据体；重复步骤五至步骤七，对弹性参数反演结果和岩相反演结果进行交替迭代更新，直至两者的反演结果稳定收敛。
具体地，本实施例的反演方法包括以下步骤：
第一步，结合测井数据中的岩相曲线，统计岩相的先验分布和岩相概率转移矩阵P；
第二步，结合测井数据中的弹性参数，利用指数函数从测井数据中拟合所有岩相的纵波阻抗低频趋势、纵横波速度比低频趋势、密度低频趋势，并计算所有岩相的弹性参数协方差矩阵；如图3所示，其中“*”代表真实的测井弹性参数，“---”代表基于不同岩相的弹性参数拟合的低频趋势，红色代表气砂，蓝色代表水砂，绿色代表致密砂，黑色代表泥岩。
第三步，输入M个角度叠加地震数据和M个角度子波，并构建角度叠加地震数据向量d和角度子波褶积矩阵；如图4所示，其中图4(a)为小角度叠加地震数据(反射角为11度)、图4(b)为中角度叠加地震数据(反射角为14度)、图4(c)为大角度叠加地震数据(反射角为17度)。
第四步，根据步骤一获取的岩相先验分布中概率最大的岩相赋值到工区中每一道的每一个采样点中，作为初始化岩相数据体；
第五步，根据前一次迭代的岩相数据体、步骤二获取的所有岩相的弹性参数低频趋势和协方差矩阵以及步骤三输入的角度叠加数据和角度子波褶积矩阵，利用岩相约束的贝叶斯线性反演计算本次迭代的弹性参数反演结果；
岩相约束的贝叶斯线性反演获取的弹性参数反演结果可以表示为：其中/> πI为纵波阻抗自然对数反演结果，πγ为纵横波速度比自然对数反演结果，πρ为密度自然对数反演结果； πIT、πγT和πρT分别表示利用前一次迭代更新的岩相π确定的纵波阻抗自然对数低频模型、纵横波速度比自然对数低频模型和密度自然对数低频模型，以πIT为例，假如第t个采样点的岩相标签为k，那么πIT的第t个元素为πIT(t)＝kIT(t)，πγT和πρT的表达方式类似；πΣm表示利用前一次迭代更新的岩相π确定的参数协方差矩阵，为3T行3T列稀疏矩阵(T为一道的采样点个数)，假如第t个采样点的岩相标签为k，那么矩阵元素/> πΣm(t,t+3T)＝kcIγ kσI kσγ、πΣm(t,t+6T)＝kcIρ kσI kσρ、πΣm(t+3T,t+6T)＝kcγρ kσγ kσρ、πΣm(t+3T,t)＝kcIγ kσI kσγ、πΣm(t+6T,t)＝kcIρ kσI kσρ和πΣm(t+6T,t+3T)＝kcγρ kσγ kσρ；G表示由子波褶积矩阵、AVO近似公式系数以及差分矩阵确定的AVO正演矩阵；/>表示噪声的方差，I表示单位矩阵；
步骤六，根据步骤四获取的弹性参数反演结果、步骤二获取的所有岩相的弹性参数低频趋势以及协方差矩阵，计算工区所有道、所有采样点弹性参数反演结果关于岩相的似然函数；
弹性参数关于岩相的似然函数可以表示为： 表示第t个采样点的弹性参数反演结果关于第t个采样点的岩相的似然函数，/>表示第t个采样点的弹性参数反演结果向量，/>表示第t个采样点弹性参数低频趋势向量，其中tπ＝1,2,L K；
步骤七，利用修正的Viterbi算法，计算岩相关于弹性参数反演结果的最大后验概率，进而获取本次迭代的岩相数据体。
修正的Viterbi算法具体实施过程：
(1)建立三维工区内所有道的随机访问顺序；
(2)按照访问顺序选择一道，根据周边相邻四道的岩相信息以及当前道的似然函数信息，利用常规Viterbi算法更新当前道的岩相；
当前道岩相关于弹性参数反演结果的后验概率为：其中/>表示当前道的后验概率，xπ表示当前道的岩相向量，/>表示当前道的弹性参数反演结果向量，δxπ表示与当前道相邻的四道岩相向量；p(x,tπ|x,t+1π)表示当前道的第t个采样点的岩相x,tπ关于当前道的第t+1个采样点的岩相x,t+1π的条件概率，为概率转移矩阵P的第t+1行、第t列元素；V(x,tπ,δx,tπ)表示势函数，表征同一采样点的相邻岩相之间的连续性，其中δx,tπ表示与当前道相邻四道的第t个采样点的岩相，/>β表示耦合系数，I(x,tπ＝i,tπ)表示指示函数，即如果当前道的第t个采样点的岩相与第i道(其中i属于当前道的邻域)的第t个采样点的相同，那么为1，否则为0，；/>表示当前道的第t个采样点的弹性参数反演结果/>关于当前道的第t个采样点的岩相的似然函数；当计算当前道岩相xπ时，其相邻道岩相δxπ是已知的，因此可以利用常规Viterbi算法计算后验概率的最大值，再进行回溯更新当前道的岩相；
(3)重复(2)更新工区内所有道的岩相；
(4)重复(1)至(3)，进行多次迭代，直至工区的岩相比例稳定；
步骤八，重复步骤五到步骤七，对弹性参数与岩相进行迭代更新，直至两者的反演结果稳定收敛。
如图5所示，其中图5(a)为岩相反演结果，井A为验证井，可以看出井旁道岩相反演结果与井A的岩相吻合度较高，4套含气砂岩均被预测出来，且岩相剖面的细节丰富、横向连续性高，具有明确的地质意义；图5(b)、图5(c)和图5(d)分别为纵波阻抗反演结果、纵横波速度比反演结果以及密度反演结果，这三种弹性参数井旁道反演结果与井A弹性参数吻合度也非常好，验证了弹性参数反演结果的精度较高。
实施例2
本实施例提供一种电子设备，包括：存储器，存储有可执行指令；处理器，处理器运行存储器中的可执行指令，以实现上述的基于马尔科夫模型的岩相和弹性参数同步反演方法。
根据本公开实施例的电子设备包括存储器和处理器。
该存储器用于存储非暂时性计算机可读指令。具体地，存储器可以包括一个或多个计算机程序产品，该计算机程序产品可以包括各种形式的计算机可读存储介质，例如易失性存储器和/或非易失性存储器。该易失性存储器例如可以包括随机存取存储器(RAM)和/或高速缓冲存储器(cache)等。该非易失性存储器例如可以包括只读存储器(ROM)、硬盘、闪存等。
该处理器可以是中央处理单元(CPU)或者具有数据处理能力和/或指令执行能力的其它形式的处理单元，并且可以控制电子设备中的其它组件以执行期望的功能。在本公开的一个实施例中，该处理器用于运行该存储器中存储的该计算机可读指令。
本领域技术人员应能理解，为了解决如何获得良好用户体验效果的技术问题，本实施例中也可以包括诸如通信总线、接口等公知的结构，这些公知的结构也应包含在本公开的保护范围之内。
有关本实施例的详细说明可以参考前述各实施例中的相应说明，在此不再赘述。
实施例4
本公开实施例提供一种计算机可读存储介质，该计算机可读存储介质存储有计算机程序，该计算机程序被处理器执行时实现所述的基于马尔科夫模型的岩相和弹性参数同步反演方法。
根据本公开实施例的计算机可读存储介质，其上存储有非暂时性计算机可读指令。当该非暂时性计算机可读指令由处理器运行时，执行前述的本公开各实施例方法的全部或部分步骤。
上述计算机可读存储介质包括但不限于：光存储介质(例如：CD－ROM和DVD)、磁光存储介质(例如：MO)、磁存储介质(例如：磁带或移动硬盘)、具有内置的可重写非易失性存储器的媒体(例如：存储卡)和具有内置ROM的媒体(例如：ROM盒)。
本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。
Claims (8)
1.一种基于马尔科夫模型的岩相和弹性参数同步反演方法，其特征在于，包括：
初始化待反演工区的岩相；
利用岩相约束的贝叶斯线性反演计算弹性参数，并计算弹性参数关于岩相的似然函数；
利用修正的Viterbi算法计算岩相关于弹性参数的后验概率最大值，并更新岩相，反复迭代直至收敛；
其中，所述修正的Viterbi算法包括：
(1)建立三维工区内所有道的随机访问顺序；
(2)按照访问顺序选择一道，根据周边相邻四道的岩相信息以及当前道的似然函数信息，利用常规Viterbi算法计算当前道岩相关于弹性参数反演结果后验概率的最大值，再进行回溯更新当前道的岩相；
(3)重复(2)更新工区内所有道的岩相；
(4)重复(1)至(3)，进行多次迭代，直至工区的岩相比例稳定；
其中，所述当前道岩相关于弹性参数反演结果的后验概率为：
其中表示当前道的后验概率，xπ表示当前道的岩相向量，/>表示当前道的弹性参数反演结果向量，δxπ表示与当前道相邻的四道岩相向量；p(x,tπ|x,t+1π)表示当前道的第t个采样点的岩相x,tπ关于当前道的第t+1个采样点的岩相x,t+1π的条件概率，为概率转移矩阵P的第t+1行、第t列元素；V(x,tπ,δx,tπ)表示势函数，表征同一采样点的相邻岩相之间的连续性，其中δx,tπ表示与当前道相邻四道的第t个采样点的岩相，β表示耦合系数，I(x,tπ＝i,tπ)表示指示函数；表示当前道的第t个采样点的弹性参数反演结果/>关于当前道的第t个采样点的岩相的似然函数。
2.根据权利要求1所述的基于马尔科夫模型的岩相和弹性参数同步反演方法，其特征在于，所述初始化待反演工区的岩相包括：结合测井数据中的岩相曲线，统计岩相的先验分布和岩相概率转移矩阵P，将岩相先验分布中概率最大的岩相赋值到工区中每一道的每一个采样点中，作为初始化岩相数据体。
3.根据权利要求1所述的基于马尔科夫模型的岩相和弹性参数同步反演方法，其特征在于，利用岩相约束的贝叶斯线性反演计算弹性参数包括：
结合测井数据中的弹性参数，利用指数函数从测井数据中拟合所有岩相的纵波阻抗低频趋势kIT、纵横波速度比低频趋势kγT、密度低频趋势kρT，并计算所有岩相的弹性参数协方差矩阵kΣm；
输入M个角度叠加地震数据和M个角度子波，并构建角度叠加地震数据向量d和角度子波褶积矩阵；
利用岩相约束的贝叶斯线性反演计算本次迭代的弹性参数反演结果。
4.根据权利要求3所述的基于马尔科夫模型的岩相和弹性参数同步反演方法，其特征在于，岩相约束的贝叶斯线性反演获取的弹性参数反演结果为：
πμ＝πμT+(G·πΣm)T(G·πΣm·GT+σn 2I)-1(d-G·πμT)
其中， πI为纵波阻抗自然对数反演结果，πγ为纵横波速度比自然对数反演结果，πρ为密度自然对数反演结果；/> πIT、πγT和πρT分别表示利用前一次迭代更新的岩相π确定的纵波阻抗自然对数低频模型、纵横波速度比自然对数低频模型和密度自然对数低频模型；πΣm表示利用前一次迭代更新的岩相π确定的参数协方差矩阵；G表示由子波褶积矩阵、AVO近似公式系数以及差分矩阵确定的AVO正演矩阵；/>表示噪声的方差，I表示单位矩阵。
5.根据权利要求3所述的基于马尔科夫模型的岩相和弹性参数同步反演方法，其特征在于，弹性参数关于岩相的似然函数采用三变量高斯分布描述：
其中，表示第t个采样点的弹性参数反演结果关于第t个采样点的岩相的似然函数，/>表示第t个采样点的弹性参数反演结果向量，/>表示第t个采样点弹性参数低频趋势向量，其中tπ＝1,2,…K。
6.一种基于马尔科夫模型的岩相和弹性参数同步反演装置，其特征在于，包括：
初始化模块，初始化待反演工区的岩相；
计算模块，利用岩相约束的贝叶斯线性反演计算弹性参数，并计算弹性参数关于岩相的似然函数；
更新迭代模块，利用修正的Viterbi算法计算岩相关于弹性参数的后验概率最大值，并更新岩相，反复迭代直至收敛；
其中，所述修正的Viterbi算法包括：
(1)建立三维工区内所有道的随机访问顺序；
(2)按照访问顺序选择一道，根据周边相邻四道的岩相信息以及当前道的似然函数信息，利用常规Viterbi算法计算当前道岩相关于弹性参数反演结果后验概率的最大值，再进行回溯更新当前道的岩相；
(3)重复(2)更新工区内所有道的岩相；
(4)重复(1)至(3)，进行多次迭代，直至工区的岩相比例稳定；
其中，所述当前道岩相关于弹性参数反演结果的后验概率为：
其中表示当前道的后验概率，xπ表示当前道的岩相向量，/>表示当前道的弹性参数反演结果向量，δxπ表示与当前道相邻的四道岩相向量；p(x,tπ|x,t+1π)表示当前道的第t个采样点的岩相x,tπ关于当前道的第t+1个采样点的岩相x,t+1π的条件概率，为概率转移矩阵P的第t+1行、第t列元素；V(x,tπ,δx,tπ)表示势函数，表征同一采样点的相邻岩相之间的连续性，其中δx,tπ表示与当前道相邻四道的第t个采样点的岩相，β表示耦合系数，I(x,tπ＝i,tπ)表示指示函数；表示当前道的第t个采样点的弹性参数反演结果/>关于当前道的第t个采样点的岩相的似然函数。
7.一种电子设备，其特征在于，包括：
存储器，存储有可执行指令；
处理器，所述处理器运行所述存储器中的所述可执行指令，以实现权利要求1-5中任一项所述的基于马尔科夫模型的岩相和弹性参数同步反演方法。
8.一种计算机可读存储介质，其特征在于，该计算机可读存储介质存储有计算机程序，该计算机程序被处理器执行时实现权利要求1-5中任一项所述的基于马尔科夫模型的岩相和弹性参数同步反演方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010968412.XA CN114185090B (zh) | 2020-09-15 | 2020-09-15 | 岩相和弹性参数同步反演方法、装置、电子设备及介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010968412.XA CN114185090B (zh) | 2020-09-15 | 2020-09-15 | 岩相和弹性参数同步反演方法、装置、电子设备及介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114185090A CN114185090A (zh) | 2022-03-15 |
| CN114185090B true CN114185090B (zh) | 2024-04-05 |
Family
ID=80539753
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010968412.XA Active CN114185090B (zh) | 2020-09-15 | 2020-09-15 | 岩相和弹性参数同步反演方法、装置、电子设备及介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114185090B (zh) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5764515A (en) * | 1995-05-12 | 1998-06-09 | Institute Francais Du Petrole | Method for predicting, by means of an inversion technique, the evolution of the production of an underground reservoir |
| CN102508293A (zh) * | 2011-11-28 | 2012-06-20 | 中国石油大学(北京) | 一种叠前反演的薄层含油气性识别方法 |
| CN104297785A (zh) * | 2014-09-29 | 2015-01-21 | 中国石油天然气股份有限公司 | 岩相约束储层物性参数反演方法及装置 |
| WO2016041189A1 (zh) * | 2014-09-19 | 2016-03-24 | 杨顺伟 | 一种评价页岩气储层及寻找甜点区的方法 |
| CN107462927A (zh) * | 2016-06-03 | 2017-12-12 | 中国石油化工股份有限公司 | 基于朴素贝叶斯分类的地震岩相预测方法和装置 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CA3043310C (en) * | 2016-12-02 | 2022-07-26 | Ratnanabha SAIN | Method for estimating petrophysical properties for single or multiple scenarios from several spectrally variable seismic and full wavefield inversion products |
| US11163080B2 (en) * | 2018-05-18 | 2021-11-02 | Repsol Exploración, S.A | Computer implemented method for generating a subsurface rock and/or fluid model of a determined domain |
-
2020
- 2020-09-15 CN CN202010968412.XA patent/CN114185090B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5764515A (en) * | 1995-05-12 | 1998-06-09 | Institute Francais Du Petrole | Method for predicting, by means of an inversion technique, the evolution of the production of an underground reservoir |
| CN102508293A (zh) * | 2011-11-28 | 2012-06-20 | 中国石油大学(北京) | 一种叠前反演的薄层含油气性识别方法 |
| WO2016041189A1 (zh) * | 2014-09-19 | 2016-03-24 | 杨顺伟 | 一种评价页岩气储层及寻找甜点区的方法 |
| CN104297785A (zh) * | 2014-09-29 | 2015-01-21 | 中国石油天然气股份有限公司 | 岩相约束储层物性参数反演方法及装置 |
| CN107462927A (zh) * | 2016-06-03 | 2017-12-12 | 中国石油化工股份有限公司 | 基于朴素贝叶斯分类的地震岩相预测方法和装置 |
Non-Patent Citations (1)
| Title |
|---|
| 《基于弹性参数加权统计的地震岩相预测方法》;桂金咏等;地球物理学报;第63卷(第1期);第300-307页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114185090A (zh) | 2022-03-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| EP3171203B1 (en) | Adaptive ensemble-based method and device for highly-nonlinear problems | |
| CN110031896B (zh) | 基于多点地质统计学先验信息的地震随机反演方法及装置 | |
| de la Varga et al. | Structural geologic modeling as an inference problem: A Bayesian perspective | |
| CN112017289B (zh) | 一种基于深度学习的井震联合初始岩性模型构建方法 | |
| Azevedo et al. | Generative adversarial network as a stochastic subsurface model reconstruction | |
| CA2919543A1 (en) | System and method for estimating a reservoir parameter using joint stochastic inversion of multisource geophysical data | |
| US11055908B2 (en) | Inversion of large, nearly-homogeneous geobodies via ultra low-dimensional shape representation using orthogonal basis functions | |
| CN110927791B (zh) | 基于深度学习利用地震数据进行流体预测的方法及装置 | |
| CN112733449A (zh) | 一种cnn井震联合反演方法、系统、存储介质、设备及应用 | |
| CN112965103B (zh) | 一种多重孔隙储层叠前地震概率化多道反演方法 | |
| Chen et al. | Integration of principal-component-analysis and streamline information for the history matching of channelized reservoirs | |
| CN116047583A (zh) | 基于深度卷积神经网络的自适应波阻抗反演方法及系统 | |
| Yao et al. | Optimized algorithm for multipoint geostatistical facies modeling based on a deep feedforward neural network | |
| CN111273346B (zh) | 去除沉积背景的方法、装置、计算机设备及可读存储介质 | |
| CN114185090B (zh) | 岩相和弹性参数同步反演方法、装置、电子设备及介质 | |
| CN112444850A (zh) | 地震资料速度建模方法、存储介质及计算设备 | |
| US20230140656A1 (en) | Method and system for determining seismic processing parameters using machine learning | |
| Koochak et al. | A variability aware GAN for improving spatial representativeness of discrete geobodies | |
| CN113419278B (zh) | 一种基于状态空间模型与支持向量回归的井震联合多目标同时反演方法 | |
| CN116702577A (zh) | 基于高斯-马尔科夫先验约束的叠后随机地震反演方法 | |
| Zhang | Ensemble methods of data assimilation in porous media flow for non-Gaussian prior probability density | |
| CN115685314A (zh) | 一种地震储层物性参数预测方法及装置 | |
| CN113970789A (zh) | 全波形反演方法、装置、存储介质及电子设备 | |
| Xue | Novel stochastic inversion methods and workflow for reservoir characterization and monitoring | |
| Lee | Joint Integration of Time-Lapse Seismic, Electromagnetic and Production Data for Reservoir Monitoring and Management |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |