CN112099079B - 自适应分频串联反射率反演方法及系统 - Google Patents

自适应分频串联反射率反演方法及系统 Download PDF

Info

Publication number
CN112099079B
CN112099079B CN201910527924.XA CN201910527924A CN112099079B CN 112099079 B CN112099079 B CN 112099079B CN 201910527924 A CN201910527924 A CN 201910527924A CN 112099079 B CN112099079 B CN 112099079B
Authority
CN
China
Prior art keywords
frequency
seismic
gather
low
inversion
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
Application number
CN201910527924.XA
Other languages
English (en)
Other versions
CN112099079A (zh
Inventor
刘俊州
时磊
刘兰锋
温立峰
王震宇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Sinopec Exploration and Production Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Exploration and Production Research Institute
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, Sinopec Exploration and Production Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201910527924.XA priority Critical patent/CN112099079B/zh
Publication of CN112099079A publication Critical patent/CN112099079A/zh
Application granted granted Critical
Publication of CN112099079B publication Critical patent/CN112099079B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. analysis, for interpretation, for correction
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/282Application of seismic models, synthetic seismograms
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/30Analysis
    • G01V1/306Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles

Abstract

公开了一种自适应分频串联反射率反演方法及系统。该方法可以包括:针对地震叠前道集,通过自适应广义S变换,分别获得低频地震叠前道集、中频地震叠前道集、高频地震叠前道集;根据已知井数据及层位数据,建立低频地震初始模型;根据低频地震初始模型、低频地震叠前道集、中频地震叠前道集、高频地震叠前道集,计算高频弹性参数,即为最终反演结果。本发明通过自适应广义S变换与串联反演,实现反演结果与实际井曲线基本吻合,提高反演结果的分辨率、精度与效率。

Description

自适应分频串联反射率反演方法及系统
技术领域
本发明涉及地球物理勘探技术领域,更具体地,涉及一种自适应分频串联反射率反演方法及系统。
背景技术
AVO技术自提出以来逐步发展成为储层识别和烃类检测最主要、最有效的方法。目前,地震反演研究主要涉及基于精确Zoeppritz方程及其近似式的反演方法和全波形反演方法。
受储层非均质强以及调谐效应等因素的影响,薄互储层的地震响应非常复杂,具有多次波、转换波,且透射损失严重等特征。基于Zoeppritz及其近似方程的常规反演方法均是假设接收的地震数据只包含一次反射波,忽略了波传播效应,例如多次波、多种类型的转换波及透射损失等,没有考虑储层厚度对AVO响应的影响。同时,Zoeppritz方程近似式假设界面两侧弹性参数变化小并且小角度入射,当界面两侧阻抗差异或者地震数据偏移距比较大时,基于近似式的AVO技术误差较大。因此,基于精确Zoeppritz方程或者其近似式来进行薄互储层的反演需要考虑大偏移距数据对反演精度的影响,且在反演前需要进行精细的振幅校正,例如透射损失补偿,多次波压制等。然而,大偏移距数据有利于提高密度反演的精度,而透射损失补偿在薄互层地震反射中难以实现,且薄互层多次波(尤其是层间多次波)的压制,在不损伤初始反射的情况下,难以实现。因此,面向薄互储层,基于精确Zoeppritz及其近似方程的传统AVO反演方法存在计算精度低的问题。
理论上,弹性波动方程反演是最为适合的薄互储层弹性参数的方法,但受计算效率、效益、反演稳定性和地震资料质量等因素制约,目前三维弹性波动方程反演还难以在实际生产中推广。
总之,基于地震数据的储层弹性参数反演方法存在以下几方面的问题:
1、基于精确Zoeppritz方程及其近似式的传统AVO反演方法,假设地震数据只包含一次反射波,忽略了波传播效应的影响,面向薄互层储层,该方法存在计算精度低的问题;
2、油气藏储层特别是针对页岩气储层最常见的结构是泥砂岩薄互层结构,地震响应多次波发育,各层反射波之间相互干涉叠加,改变了常规AVO响应特征,传统的AVO反演方法不再适用;
3、全波形反演虽然能够利用全波场的信息预测储层参数,但其计算量大,在反演尺度和计算效率上不能满足实际油藏储层精细表征要求;
4、针对低频能量强、高频能量弱的地震资料,常规反演策略反演结果无法体现高频细节信息,其反演结果分辨率有待进一步提高。
因此,如何开发出一种新的AVO储层反演处理策略,有效获得地下介质弹性参数模型信息作为先验信息,以降低反演的不确定性,提高反演的精度以及反演结果分辨率体现高频细节信息是本领域亟待解决的技术难题。
利用反射率法来求解弹性波动方程是1D假设下的全波场模拟方法。它能够模拟多次波、转换波,考虑了薄层厚度以及透射损失对地震数据的影响,克服传统AVO技术存在的问题,适合于薄互层的地震波场模拟分析与反演。反射率法首先由Thomson(1950)提出,后来由Fuchs and Müller(1971),Kennett(1979,1983),Kennett and Kerry(1979),Fryer(1980),Muller(1985)这些学者进行研究及改进。其中应用最多的反射率法是Kennett(1983)提出的递归矩阵算法,并且后来大多数应用反射率法的反演都是基于此理论或者对其改进得来,所以该类反射率法被称为KRM(Kennett Reflectivity Method)。然而,KRM的递归算法相对复杂,计算量大且耗时。所以Phinney等(1987)提出了一种向量化计算递归矩阵的反射率法,该方法很大程度上简化了计算过程,节约计算时间。但是该方法求导困难且正演求取的道集是不经过动校正的。因此,本技术方案对该向量化的反射率法进行改进并引入后面的反演理论中。
虽然利用反射率法求解弹性波动方程的研究很多,但其反演问题的研究却很少。1994年,Huasheng和Ursin等基于Kennett的反射率法,采用Levenberg-Marquardt算法进行反演。1998年,Gouveia和Scales基于Fuchs和Muller(1971)的反射率法求解波动方程,引入贝叶斯理论并用共轭梯度法寻找最优解。2003年,Sen和Roy基于Kennett的反射率法进行正演,用梯度下降法求解目标函数,并用共轭梯度迭代寻优。2015年,Mallick和Adhikari采用反射率法求解波动方程,然后基于贝叶斯框架用遗传算法进行反演。虽然这些反演方法都很成功,但是都是基于KRM的反演,其计算过程相对复杂、耗时。为此,刘洪星(2016)基于高斯先验分布利用向量化的反射率法来反演三参数。虽然该方法计算效率大大高于KRM系列的计算方法,但该反演利用高斯先验分布约束反演,其反演的分辨率有待提高。因此,有必要开发一种自适应分频串联反射率反演方法及系统。
公开于本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解,而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。
发明内容
本发明提出了一种自适应分频串联反射率反演方法及系统,其能够通过自适应广义S变换与串联反演,实现反演结果与实际井曲线基本吻合,提高反演结果的分辨率、精度与效率。
根据本发明的一方面,提出了一种自适应分频串联反射率反演方法。所述方法可以包括:针对地震叠前道集,通过自适应广义S变换,分别获得低频地震叠前道集、中频地震叠前道集、高频地震叠前道集;根据所述低频、中频、高频地震叠前道集,建立低频地震初始模型;根据所述低频地震初始模型、所述低频地震叠前道集、所述中频地震叠前道集、所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果。
优选地,所述根据所述低频地震初始模型、所述低频地震叠前道集、所述中频地震叠前道集、所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果包括:根据所述低频地震初始模型与所述低频地震叠前道集,计算中频地震初始模型;根据所述中频地震初始模型与所述中频地震叠前道集,计算高频地震初始模型;根据所述高频地震初始模型与所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果。
优选地,进一步包括:根据对应频带的地震初始模型与测井数据,获得模型参数;根据所述模型参数,构建协方差矩阵;构建所述模型参数的先验分布函数;根据反演似然函数与所述先验分布函数,获得后验概率分布函数,进而确定反演的目标函数;根据所述协方差矩阵与所述目标函数,计算对应频带的弹性参数,作为高一级别频带的地震初始模型。
优选地,所述根据所述模型参数,构建协方差矩阵为:所述模型参数包括横波速度、纵波速度与密度,求取各个参数的自相关系数与互相关系数,构建所述协方差矩阵为:
Figure BDA0002098815190000041
其中,Cm为协方差矩阵,δvp为纵波速度方差,δvs为横波速度方差,δrho为密度方差。
优选地,所述模型参数的先验分布函数为:
Figure BDA0002098815190000051
其中,n为噪音,Cd是噪声的协方差矩阵,
Figure BDA0002098815190000052
Figure BDA0002098815190000054
为噪声方差,I为单位矩阵。
优选地,所述后验概率分布函数为:
P(m|d)∝P(d|m)P(m) (3)
其中,P(m|d)为后验概率分布,P(d|m)是似然函数,P(m)是先验分布。
优选地,所述目标函数为:
Figure BDA0002098815190000053
其中,d为观测的地震叠前道集,d=G(m)+n,G为正演算子,G(m)为合成角度道集。
根据本发明的另一方面,提出了一种自适应分频串联反射率反演系统,其特征在于,该系统包括:存储器,存储有计算机可执行指令;处理器,所述处理器运行所述存储器中的计算机可执行指令,执行以下步骤:针对地震叠前道集,通过自适应广义S变换,分别获得低频地震叠前道集、中频地震叠前道集、高频地震叠前道集;根据已知井数据及层位数据,建立低频地震初始模型;根据所述低频地震初始模型、所述低频地震叠前道集、所述中频地震叠前道集、所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果。
优选地,所述根据所述低频地震初始模型、所述低频地震叠前道集、所述中频地震叠前道集、所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果包括:根据所述低频地震初始模型与所述低频地震叠前道集,计算中频地震初始模型;根据所述中频地震初始模型与所述中频地震叠前道集,计算高频地震初始模型;根据所述高频地震初始模型与所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果。
优选地,进一步包括:根据对应频带的地震初始模型与测井数据,获得模型参数;根据所述模型参数,构建协方差矩阵;构建所述模型参数的先验分布函数;根据反演似然函数与所述先验分布函数,获得后验概率分布函数,进而确定反演的目标函数;根据所述协方差矩阵与所述目标函数,计算对应频带的弹性参数,作为高一级别频带的地震初始模型。
优选地,所述根据所述模型参数,构建协方差矩阵为:所述模型参数包括横波速度、纵波速度与密度,求取各个参数的自相关系数与互相关系数,构建所述协方差矩阵为:
Figure BDA0002098815190000061
其中,Cm为协方差矩阵,δvp为纵波速度方差,δvs为横波速度方差,δrho为密度方差。
优选地,所述模型参数的先验分布函数为:
Figure BDA0002098815190000062
其中,n为噪音,Cd是噪声的协方差矩阵,
Figure BDA0002098815190000063
Figure BDA0002098815190000065
为噪声方差,I为单位矩阵。
优选地,所述后验概率分布函数为:
P(m|d)∝P(d|m)P(m) (3)
其中,P(m|d)为后验概率分布,P(d|m)是似然函数,P(m)是先验分布。
优选地,所述目标函数为:
Figure BDA0002098815190000064
其中,d为观测的地震叠前道集,d=G(m)+n,G为正演算子,G(m)为合成角度道集。
其有益效果在于:
(1)相比于基于Zoeppritz及其近似式的传统反演方法,反射率法除了可以模拟一次反射,还能够模拟透射波、多种转换波和多次波,考虑了地震波相位变化、地层厚度以及透射损失的影响,反演精度得到很好的提高;
(2)相比于受计算效率、经济效益、反演稳定性和地震资料质量等因素制约的弹性波动方程反演,反射率法计算具有成本低,更高的计算效率等特点;
(3)通过自适应广义S变换将地震叠前道集分成不同频带范围,然后进行分频串联反演,相较传统的反演策略进一步提高了反演结果的分辨率和精度;
(4)通过反射率反演过程中矩阵存储策略,进一步提高了反演效率;
(5)通过服从微分拉普拉斯的垂向块约束,可以很好的突出薄互储层的边界特征,提高反演精度。
本发明的方法和装置具有其它的特性和优点,这些特性和优点从并入本文中的附图和随后的具体实施方式中将是显而易见的,或者将在并入本文中的附图和随后的具体实施方式中进行详细陈述,这些附图和具体实施方式共同用于解释本发明的特定原理。
附图说明
通过结合附图对本发明示例性实施例进行更详细的描述,本发明的上述以及其它目的、特征和优势将变得更加明显,其中,在本发明示例性实施例中,相同的参考标号通常代表相同部件。
图1示出了根据本发明的自适应分频串联反射率反演方法的步骤的流程图。
图2示出了根据本发明的一个实施例的低频、中频、高频地震子波的示意图。
图3示出了根据本发明的一个实施例的地震叠前AVA道集的示意图。
图4a、图4b和图4c分别示出了根据本发明的一个实施例的低频、中频、高频地震叠前道集示意图。
图5a、图5b和图5c分别示出了根据本发明的一个实施例的低频纵波速度、横波速度、密度的示意图。
图6a、图6b和图6c分别示出了根据本发明的一个实施例的中频纵波速度、横波速度、密度的示意图。
图7a、图7b和图7c分别示出了根据本发明的一个实施例的高频纵波速度、横波速度、密度的示意图。
具体实施方式
下面将参照附图更详细地描述本发明。虽然附图中显示了本发明的优选实施例,然而应该理解,可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反,提供这些实施例是为了使本发明更加透彻和完整,并且能够将本发明的范围完整地传达给本领域的技术人员。
图1示出了根据本发明的自适应分频串联反射率反演方法的步骤的流程图。
在该实施例中,根据本发明的自适应分频串联反射率反演方法可以包括:步骤101,针对地震叠前道集,通过自适应广义S变换,分别获得低频地震叠前道集、中频地震叠前道集、高频地震叠前道集;步骤102,根据低频、中频、高频地震叠前道集,建立低频地震初始模型;步骤103,根据低频地震初始模型、低频地震叠前道集、中频地震叠前道集、高频地震叠前道集,计算高频弹性参数,即为最终反演结果。
在一个示例中,根据低频地震初始模型、低频地震叠前道集、中频地震叠前道集、高频地震叠前道集,计算高频弹性参数,即为最终反演结果包括:根据低频地震初始模型与低频地震叠前道集,计算中频地震初始模型;根据中频地震初始模型与中频地震叠前道集,计算高频地震初始模型;根据高频地震初始模型与高频地震叠前道集,计算高频弹性参数,即为最终反演结果。
在一个示例中,进一步包括:根据对应频带的地震初始模型与测井数据,获得模型参数;根据模型参数,构建协方差矩阵;构建模型参数的先验分布函数;根据反演似然函数与先验分布函数,获得后验概率分布函数,进而确定反演的目标函数;根据协方差矩阵与目标函数,计算对应频带的弹性参数,作为高一级别频带的地震初始模型。
在一个示例中,根据模型参数,构建协方差矩阵为:模型参数包括横波速度、纵波速度与密度,求取各个参数的自相关系数与互相关系数,构建协方差矩阵为:
Figure BDA0002098815190000091
其中,Cm为协方差矩阵,δvp为纵波速度方差,δvs为横波速度方差,δrho为密度方差。
在一个示例中,模型参数的先验分布函数为:
Figure BDA0002098815190000092
其中,n为噪音,Cd是噪声的协方差矩阵,
Figure BDA0002098815190000093
Figure BDA0002098815190000094
为噪声方差,I为单位矩阵。
在一个示例中,后验概率分布函数为:
P(m|d)∝P(d|m)P(m) (3)
其中,P(m|d)为后验概率分布,P(d|m)是似然函数,P(m)是先验分布。
在一个示例中,目标函数为:
Figure BDA0002098815190000101
其中,d为观测的地震叠前道集,d=G(m)+n,G为正演算子,G(m)为合成角度道集。
具体地,根据本发明的自适应分频串联反射率反演方法可以包括:
针对地震叠前道集,通过公式(5)进行自适应广义S正变换:
Figure BDA0002098815190000102
进而通过公式(6)进行反变换,获得低频、中频、高频地震叠前道集:
Figure BDA0002098815190000103
图2示出了根据本发明的一个实施例的低频、中频、高频地震子波的示意图。
反演过程中,不同地震记录对应不同的地震子波,只有子波对应的情况下,反演合成的地震记录才可以很好的与实际地震记录吻合。因此,根据低频、中频、高频地震叠前道集,获得低频、中频、高频地震子波,如图2所示,将提取的子波和已知测井数据作为反射率正演的输入,正演得到地震记录与井旁地震道对比,调节振幅比例因子,调整输入子波振幅,再正演对比井旁地震道,确保得到合适的振幅比例因子。通常情况下,子波包括波形和振幅两个主要信息,提取的子波一般振幅与实际地震记录振幅存在差异,所以需要振幅比例因子对子波振幅进行调节,使得子波与反演结果作用得到的正演结果与实际地震记录达到最佳匹配程度。只有子波信息与实际地震记录匹配的情况下,才可以反演得到一个较好的结果。
根据低频、中频、高频地震叠前道集,建立低频地震初始模型。根据低频地震初始模型与低频地震叠前道集,计算中频地震初始模型;根据中频地震初始模型与中频地震叠前道集,计算高频地震初始模型;根据高频地震初始模型与高频地震叠前道集,计算高频弹性参数,即为最终反演结果。
在每一步计算中均为:
根据对应频带的地震初始模型与测井数据,获得模型参数,基于贝叶斯理论引入模型参数的先验信息并建立了通用的反演框架,后验概率分布函数为公式(3),有P(m)=P0mexp{-F(m)},其中P0m视为常数项,F(m)是反演正则化项:
Figure BDA0002098815190000111
其中m=[α1…αN1…βN1…ρN]T为模型参数,μ为模型参数背景值,
Figure BDA0002098815190000115
为正则化因子,N为数据长度。
如果地震叠前道集中所含的噪声相互独立且满足高斯分布,那么模型参数的先验分布函数为公式(2),联合d=G(m)+n和式(3)可求取到似然函数,其表达式如下:
Figure BDA0002098815190000112
其中d=G(m)+n是观测的地震叠前道集,G为正演算子,G(m)为合成角度道集,P0为:
Figure BDA0002098815190000113
将先验分布P(m)乘上似然函数P(d|m)就能够求取后验概率分布函数,其表示式为:
Figure BDA0002098815190000114
在求得后验概率分布之后,求取后验概率分布函数关于模型参数的导数,并令其等于零,计算最大后验概率解,即求取目标函数公式(4)的极小值。
利用高斯-牛顿迭代来求解目标函数是最常见且最有效的方法,高斯-牛顿迭代公式为:
mk+1=mk-ηH(mk)-1γ(mk) (11)
其中,η为反演迭代步长,γ是S0(m)对模型参数的一阶偏导,H是S0(m)对模型参数的二阶偏导,称为海森矩阵。矩阵γ和矩阵H的表达形式如下:
Figure BDA0002098815190000121
由于海森矩阵的求取比较复杂,计算耗时较长,本论文为了降低计算成本,在反演过程中省略了海森矩阵的二阶导数项。
利用微分拉普拉斯先验分布来约束反演过程,并兼顾了反演的分辨率和数值稳定性,其指数项表式为:
Figure BDA0002098815190000122
其中D为一阶导数算子矩阵,表示为:
Figure BDA0002098815190000123
式中最后一项中的下角标j表示时间点,κl代表缩放因子,参数l=1,2,3分别代表了纵波速度Vp、横波速度Vs和密度Rho。而垂直约束算子C(x)是一个正则化项,其增强了先验约束对块数据的反演能力,其表达式为:
Figure BDA0002098815190000124
模型参数包括横波速度、纵波速度与密度,求取各个参数的自相关系数与互相关系数,构建协方差矩阵为公式(1),根据协方差矩阵与目标函数,计算对应频带的弹性参数,作为高一级别频带的地震初始模型。
由此得微分拉普拉斯先验分布对模型参数求一阶导数,可表示为:
Figure BDA0002098815190000131
其中,
Figure BDA0002098815190000132
微分拉普拉斯分布是凸函数,不仅反演分辨率高,且数值稳定性好。尤其是在反演界面两侧物性差异比较大的情况下,基于微分拉普拉斯分布的反演能够很好贴合实际数据。
建立一个3M×1的模型参数m:
m=[α1…αM1…βM1…ρM]T (18)
根据反射率法的正演理论计算τ-p域的反射系数:
Figure BDA0002098815190000133
根据褶积理论,将计算的反射系数R褶积上子波W,就能获得合成地震记录:
G(m)=WR (20)
根据前面贝叶斯反演理论可知,求解反演问题的关键就是要求取合成地震记录G(m)关于模型参数m的导数:
Figure BDA0002098815190000134
由于子波与模型参数无关,所以合成地震记录关于模型参数的导数只和反射系数关于模型参数的导数有关。而反射系数关于模型参数的导数可表达为:
Figure BDA0002098815190000141
其中
Figure BDA0002098815190000142
分别为第n层的纵横波速度以及密度。认真辨析向量ν0可知,除了Qn
Figure BDA0002098815190000143
相关,其余项都与
Figure BDA0002098815190000144
无关。所以ν0和Qn关于模型参数的一阶偏导可以表示为:
Figure BDA0002098815190000145
其中En
Figure BDA0002098815190000146
Figure BDA0002098815190000147
关于模型参数的一阶导数表达式。
在反演编程实现过程中,存储每层的中间算子,虽然增加了内存开销,但是可以提高运行效率。
在获得频率-慢度域的导数之后,时间-慢度域的偏导可以通过频率积分获得:
Figure BDA0002098815190000148
将计算获得的矩阵G(m)、矩阵
Figure BDA0002098815190000149
带入反演迭代式(11)中,就可建立基于反射率法的贝叶斯反演方法。
本方法通过自适应广义S变换与串联反演,实现反演结果与实际井曲线基本吻合,提高反演结果的分辨率、精度与效率。
应用示例
为便于理解本发明实施例的方案及其效果,以下给出一个具体应用示例。本领域技术人员应理解,该示例仅为了便于理解本发明,其任何具体细节并非意在以任何方式限制本发明。
图3示出了根据本发明的一个实施例的地震叠前AVA道集的示意图。
图4a、图4b和图4c分别示出了根据本发明的一个实施例的低频、中频、高频地震叠前道集示意图。
根据本发明的自适应分频串联反射率反演方法可以包括:
针对如图3所示的地震叠前道集,通过自适应广义S变换,分别获得低频地震叠前道集、中频地震叠前道集、高频地震叠前道集,如图4a-4c所示。
根据已知井数据及层位数据,建立低频地震初始模型。
图5a、图5b和图5c分别示出了根据本发明的一个实施例的低频纵波速度、横波速度、密度的示意图。
图6a、图6b和图6c分别示出了根据本发明的一个实施例的中频纵波速度、横波速度、密度的示意图。
图7a、图7b和图7c分别示出了根据本发明的一个实施例的高频纵波速度、横波速度、密度的示意图。
根据低频地震初始模型与低频地震叠前道集,计算获得低频弹性参数,即为中频地震初始模型,如图5a-5c所示;根据中频地震初始模型与中频地震叠前道集,计算获得中频弹性参数,即为高频地震初始模型,如图6a-6c所示;根据高频地震初始模型与高频地震叠前道集,计算高频弹性参数,即为最终反演结果,如图7a-7c所示,在每一步计算中均为:根据对应频带的地震初始模型与测井数据,获得模型参数;模型参数包括横波速度、纵波速度与密度,求取各个参数的自相关系数与互相关系数,构建协方差矩阵为公式(1);根据反演似然函数与先验分布函数公式(2),获得后验概率分布函数为公式(3),进而确定反演的目标函数为公式(4);根据协方差矩阵与目标函数,计算对应频带的弹性参数,作为高一级别频带的地震初始模型。
综上所述,本发明通过自适应广义S变换与串联反演,实现反演结果与实际井曲线基本吻合,提高反演结果的分辨率、精度与效率。
本领域技术人员应理解,上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果,并不意在将本发明的实施例限制于所给出的任何示例。
根据本发明的实施例,提供了一种自适应分频串联反射率反演系统,其特征在于,该系统包括:存储器,存储有计算机可执行指令;处理器,所述处理器运行所述存储器中的计算机可执行指令,执行以下步骤:针对地震叠前道集,通过自适应广义S变换,分别获得低频地震叠前道集、中频地震叠前道集、高频地震叠前道集;根据低频、中频、高频地震叠前道集,建立低频地震初始模型;根据低频地震初始模型、低频地震叠前道集、中频地震叠前道集、高频地震叠前道集,计算高频弹性参数,即为最终反演结果。
在一个示例中,根据低频地震初始模型、低频地震叠前道集、中频地震叠前道集、高频地震叠前道集,计算高频弹性参数,即为最终反演结果包括:根据低频地震初始模型与低频地震叠前道集,计算中频地震初始模型;根据中频地震初始模型与中频地震叠前道集,计算高频地震初始模型;根据高频地震初始模型与高频地震叠前道集,计算高频弹性参数,即为最终反演结果。
在一个示例中,进一步包括:根据对应频带的地震初始模型与测井数据,获得模型参数;根据模型参数,构建协方差矩阵;构建模型参数的先验分布函数;根据反演似然函数与先验分布函数,获得后验概率分布函数,进而确定反演的目标函数;根据协方差矩阵与目标函数,计算对应频带的弹性参数,作为高一级别频带的地震初始模型。
在一个示例中,根据模型参数,构建协方差矩阵为:模型参数包括横波速度、纵波速度与密度,求取各个参数的自相关系数与互相关系数,构建协方差矩阵为:
Figure BDA0002098815190000171
其中,Cm为协方差矩阵,δvp为纵波速度方差,δvs为横波速度方差,δrho为密度方差。
在一个示例中,模型参数的先验分布函数为:
Figure BDA0002098815190000172
其中,n为噪音,Cd是噪声的协方差矩阵,
Figure BDA0002098815190000173
Figure BDA0002098815190000175
为噪声方差,I为单位矩阵。
在一个示例中,后验概率分布函数为:
P(m|d)∝P(d|m)P(m) (3)
其中,P(m|d)为后验概率分布,P(d|m)是似然函数,P(m)是先验分布。
在一个示例中,目标函数为:
Figure BDA0002098815190000174
其中,d为观测的地震叠前道集,d=G(m)+n,G为正演算子,G(m)为合成角度道集。
本系统通过自适应广义S变换与串联反演,实现反演结果与实际井曲线基本吻合,提高反演结果的分辨率、精度与效率。
本领域技术人员应理解,上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果,并不意在将本发明的实施例限制于所给出的任何示例。
以上已经描述了本发明的各实施例,上述说明是示例性的,并非穷尽性的,并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下,对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。

Claims (6)

1.一种自适应分频串联反射率反演方法,其特征在于,包括:
针对地震叠前道集,通过自适应广义S变换,分别获得低频地震叠前道集、中频地震叠前道集、高频地震叠前道集;
根据已知井数据及层位数据,建立低频地震初始模型;
根据所述低频地震初始模型、所述低频地震叠前道集、所述中频地震叠前道集、所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果;
其中,所述根据所述低频地震初始模型、所述低频地震叠前道集、所述中频地震叠前道集、所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果包括:
根据所述低频地震初始模型与所述低频地震叠前道集,计算中频地震初始模型;
根据所述中频地震初始模型与所述中频地震叠前道集,计算高频地震初始模型;
根据所述高频地震初始模型与所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果;
其中,进一步包括:
根据对应频带的地震初始模型与测井数据,获得模型参数;
根据所述模型参数,构建协方差矩阵;
构建所述模型参数的先验分布函数;
根据反演似然函数与所述先验分布函数,获得后验概率分布函数,进而确定反演的目标函数;
根据所述协方差矩阵与所述目标函数,计算对应频带的弹性参数,作为高一级别频带的地震初始模型。
2.根据权利要求1所述的自适应分频串联反射率反演方法,其中,所述根据所述模型参数,构建协方差矩阵为:
所述模型参数包括横波速度、纵波速度与密度,求取各个参数的自相关系数与互相关系数,构建所述协方差矩阵为:
Figure FDA0003596389700000021
其中,Cm为协方差矩阵,δvp为纵波速度方差,δvs为横波速度方差,δrho为密度方差。
3.根据权利要求1所述的自适应分频串联反射率反演方法,其中,所述模型参数的先验分布函数为:
Figure FDA0003596389700000022
其中,n为噪音,Cd是噪声的协方差矩阵,
Figure FDA0003596389700000023
为噪声方差,I为单位矩阵。
4.根据权利要求1所述的自适应分频串联反射率反演方法,其中,所述后验概率分布函数为:
P(m|d)∝P(d|m)P(m) (3)
其中,P(m|d)为后验概率分布,P(d|m)是似然函数,P(m)是先验分布。
5.根据权利要求1所述的自适应分频串联反射率反演方法,其中,所述目标函数为:
Figure FDA0003596389700000031
其中,d为观测的地震叠前道集,d=G(m)+n,G为正演算子,G(m)为合成角度道集。
6.一种自适应分频串联反射率反演系统,其特征在于,该系统包括:
存储器,存储有计算机可执行指令;
处理器,所述处理器运行所述存储器中的计算机可执行指令,执行以下步骤:
针对地震叠前道集,通过自适应广义S变换,分别获得低频地震叠前道集、中频地震叠前道集、高频地震叠前道集;
根据已知井数据及层位数据,建立低频地震初始模型;
根据所述低频地震初始模型、所述低频地震叠前道集、所述中频地震叠前道集、所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果;
其中,所述根据所述低频地震初始模型、所述低频地震叠前道集、所述中频地震叠前道集、所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果包括:
根据所述低频地震初始模型与所述低频地震叠前道集,计算中频地震初始模型;
根据所述中频地震初始模型与所述中频地震叠前道集,计算高频地震初始模型;
根据所述高频地震初始模型与所述高频地震叠前道集,计算高频弹性参数,即为最终反演结果;
其中,进一步包括:
根据对应频带的地震初始模型与测井数据,获得模型参数;
根据所述模型参数,构建协方差矩阵;
构建所述模型参数的先验分布函数;
根据反演似然函数与所述先验分布函数,获得后验概率分布函数,进而确定反演的目标函数;
根据所述协方差矩阵与所述目标函数,计算对应频带的弹性参数,作为高一级别频带的地震初始模型。
CN201910527924.XA 2019-06-18 2019-06-18 自适应分频串联反射率反演方法及系统 Active CN112099079B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910527924.XA CN112099079B (zh) 2019-06-18 2019-06-18 自适应分频串联反射率反演方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910527924.XA CN112099079B (zh) 2019-06-18 2019-06-18 自适应分频串联反射率反演方法及系统

Publications (2)

Publication Number Publication Date
CN112099079A CN112099079A (zh) 2020-12-18
CN112099079B true CN112099079B (zh) 2022-08-05

Family

ID=73748448

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910527924.XA Active CN112099079B (zh) 2019-06-18 2019-06-18 自适应分频串联反射率反演方法及系统

Country Status (1)

Country Link
CN (1) CN112099079B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103245970A (zh) * 2012-02-08 2013-08-14 中国石油化工股份有限公司 叠前地震宽角度反演方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9201155B2 (en) * 2013-06-12 2015-12-01 Halliburton Energy Services, Inc. Systems and methods for downhole electromagnetic field measurement
CN104614763B (zh) * 2015-01-19 2017-06-06 中国石油大学(北京) 基于反射率法的多波avo储层弹性参数反演方法及系统
CN105182416A (zh) * 2015-09-06 2015-12-23 中国石油天然气股份有限公司 基于分频数据的地震反演方法及其装置
CN107179547A (zh) * 2017-06-06 2017-09-19 中海石油(中国)有限公司 一种地震波阻抗反演低频模型建立方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103245970A (zh) * 2012-02-08 2013-08-14 中国石油化工股份有限公司 叠前地震宽角度反演方法

Also Published As

Publication number Publication date
CN112099079A (zh) 2020-12-18

Similar Documents

Publication Publication Date Title
Kaur et al. Seismic ground‐roll noise attenuation using deep learning
Yu et al. Attenuation of noise and simultaneous source interference using wavelet denoising
EP2335093B1 (en) Estimation of soil properties using waveforms of seismic surface waves
US9625593B2 (en) Seismic data processing
US7937224B2 (en) Diplet-based seismic processing
US20120275267A1 (en) Seismic Data Processing
US20150293245A1 (en) Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media
US20120253758A1 (en) Method of Wavelet Estimation and Multiple Prediction In Full Wavefield Inversion
CN110456417B (zh) 一种地震数据多次波压制方法
Wang et al. Seismic, waveform modeling and tomography
CN110261897B (zh) 基于组稀疏的叠前四参数反演方法
CN110780351A (zh) 纵波和转换波叠前联合反演方法及系统
CN111077567B (zh) 一种基于矩阵乘法的双程波叠前深度偏移的方法
US11340366B2 (en) Accurate velocity model estimation and imaging in the presence of localized attenuation (Q) anomalies
Aghamiry et al. Complex-valued imaging with total variation regularization: an application to full-waveform inversion in visco-acoustic media
CN111665556B (zh) 地层声波传播速度模型构建方法
CN110687597A (zh) 一种基于联合字典的波阻抗反演方法
CN107621654A (zh) 一种基于最大相关熵的地震叠后波阻抗反演方法
Métivier et al. A review of the use of optimal transport distances for high resolution seismic imaging based on the full waveform
Sun et al. Multiple attenuation using λ-f domain high-order and high-resolution Radon transform based on SL0 norm
CN112099079B (zh) 自适应分频串联反射率反演方法及系统
Zhong et al. Elastic reverse time migration method in vertical transversely isotropic media including surface topography
Kazei et al. Mapping full seismic waveforms to vertical velocity profiles by deep learning
Kontakis et al. Efficient 1.5 D full waveform inversion in the Laplace-Fourier domain
Mozayan et al. Blocky inversion of multichannel elastic impedance for elastic parameters

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