本发明涉及石油地震勘探领域，尤其涉及绝对波阻抗反演成像方法、装置、设备及存储介质。
背景技术：
声阻抗(ai)是描述储层的关键弹性参数之一，是测井资料与地震观测资料之间的天然桥梁，宽频带ai反演成像是地震解释和储层描述的关键步骤。tarantola1984年《geophysics》发表的论文“inversionofseismicreflectiondataintheacousticapproximation”指出，在数学意义上，贝叶斯反演框架下的全波形反演(fwi)是估计地下宽带弹性参数反演阻抗最完美的方法。然而，对于实际地球物理数据反演，其客观条件是数据的物理采集不完整、正演方法不完备、地下初始参数模型不够精确等诸多因素，这使得fwi在实际高维弹性参数反演的应用仍然是有挑战的。
声阻抗估计的普遍方法是先获取地下反射系数，再通过一维递归反演方法获得声阻抗参数。berteussen和ursin在1983在《geophysics》发表的经典阻抗反演论文“approximatecomputationoftheacousticimpedancefromseismicdata”提出的一维递归法计算简便，目前在工业界仍有应用。但其一维实现对噪声敏感，在复杂介质中会带来严重的累积误差甚至反演不稳定。gholami于2015年在《geophysics》的论文“afastautomaticmultichannelblindseismicinversionforhigh-resolutionimpedancerecovery”提出了基于全变差的正则化连续模型反演方法，它可以克服一维递归法的不稳定的缺点，但在反演过程中没有融合背景阻抗，只能得到相对阻抗剖面。
因此，如何消除一维实现对噪声敏感，在复杂介质中会带来严重的累积误差甚至反演不稳定是亟待解决的技术问题。
上述内容仅用于辅助理解本发明的技术方案，并不代表承认上述内容是现有技术。
技术实现要素：
本发明的主要目的在于提供绝对波阻抗反演成像方法、装置、设备及存储介质，旨在解决现有技术中一维实现对噪声敏感，在复杂介质中会带来严重的累积误差甚至反演不稳定的技术问题，实现高维绝对波阻抗反演成像。
为实现上述目的，本发明提供绝对波阻抗反演成像方法，包括以下步骤：
采集勘探现场的井控数据，提取所述井控数据中阻抗的低波数部分alow；
获取地震数据，对所述地震数据进行高斯束偏移，获得反射系数r；
将所述反射系数r和所述阻抗的低波数部分alow输入绝对波阻抗反演模型进行约束计算，获得绝对波阻抗z；
通过所述绝对波阻抗z获得波阻抗剖面图像。
优选地，所述提取所述井控数据中阻抗的低波数部分alow，包括：
提取所述井控数据中的背景速度和背景密度；
将所述背景速度与所述背景密度相乘，获得阻抗的低波数部分alow。
优选地，所述对所述地震数据进行高斯束偏移，获得反射系数r，包括：
对所述地震数据进行常规高斯束偏移成像，获得偏移图像数据；
通过所述背景速度计算获得地下振幅补偿函数点扩散函数psf；
将所述psf和所述偏移图像数据通过反褶积计算，获得反射系数r；
优选地，所述将所述反射系数r和所述阻抗的低波数部分alow输入绝对波阻抗反演模型进行约束计算，获得绝对波阻抗z，包括：
根据所述反射系数r和所述阻抗的低波数部分alow通过声阻抗ai剖面计算公式进行循环计算，获得绝对波阻抗z。
优选地，所述声阻抗ai剖面计算公式具体为：
其中，alow为阻抗的低波数部分，l为低通滤波器，r为反射系数，||a||tv为一阶tv范数，a为阻抗的半自然对数，
优选地，所述根据所述反射系数r和所述阻抗的低波数部分alow通过声阻抗ai剖面计算公式进行循环计算，获得绝对波阻抗z，包括：
设置计算参数，包括：sk,
其中sk,
设置所述计算参数的初始值，并获取第一公式、第二公式和第三公式；
所述第一公式、所述第二公式和所述第三公式均从所述计算参数的初始值开始计算，令i和k的初始值均为0，
将所述计算参数带入所述第一公式进行计算，获得更新后的ak 1，所述第一公式具体为：
其中，
所述第一公式计算完成后将所述更新后的ak 1带入所述第二公式进行计算，获得更新后的sk,
所述第二公式计算完成后，k值加1，将数值加一后的k、所述更新后的sk,
循环更新a的值n次后，i的值加一，且令ai 1＝an，k的值归零；
将数值加一后的i带入所述第三公式进行计算，获得更新后的
将所述数值加一后的i、归零后的k、所述更新后的
重复更新i的值m次；
总共将所述a的值更新m*n次，获得绝对阻抗的半自然对数a′；
将所述a′带入第四公式进行计算获得绝对波阻抗z，所述第四公式具体为：
z＝exp(2*a′)
此外，为实现上述目的，本发明还提出绝对波阻抗反演成像设备，所述绝对波阻抗反演成像设备包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的绝对波阻抗反演成像程序，所述绝对波阻抗反演成像程序配置为实现如上文所述的绝对波阻抗反演成像方法的步骤。
此外，为实现上述目的，本发明还提出一种存储介质，所述存储介质上存储有绝对波阻抗反演成像程序，所述绝对波阻抗反演成像程序被处理器执行时实现如上文所述的绝对波阻抗反演成像方法的步骤。
此外，为实现上述目的，本发明还提出绝对波阻抗反演成像装置，包括：
第一数据提取模块，用于采集勘探现场的井控数据，提取所述井控数据中阻抗的低波数部分alow；
第二数据提取模块，用于获取地震数据，对所述地震数据进行高斯束偏移，获得反射系数r；
计算模块，用于将所述反射系数r和所述阻抗的低波数部分alow输入绝对波阻抗反演模型进行约束计算，获得绝对波阻抗z；
图像生成模块，用于通过所述绝对波阻抗z获得波阻抗剖面图像。
本发明中，通过井控数据提供背景阻抗和保幅成像获得反射系数，既分别得到阻抗的低波数部分和高波数部分，再融合背景阻抗和反射系数进行联合反演，开发出稳定高效的高维波绝对阻抗成像方法进行绝对波阻抗成像，指导地下岩性解释，进而推进油气勘探技术进步。
附图说明
图1为本发明实施例方案涉及的硬件运行环境的绝对波阻抗反演成像设备的结构示意图；
图2为本发明绝对波阻抗反演成像方法第一实施例的流程示意图；
图3为本发明绝对波阻抗反演成像方法第二实施例的流程示意图；
图4为本发明绝对波阻抗反演成像装置第一实施例的结构框图；
图5为本发明循环计算流程图；
图6为sigabee2a模型理论波阻抗；
图7为sigabee2a模型一维递归法反演得到的波阻抗；
图8为sigabee2a模型本发明提出的方法反演得到的波阻抗；
图9为真阻抗与道积分法反演得到的阻抗抽道对比；
图10为真阻抗与本发明提出的方法反演得到的阻抗抽道对比；
本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
参照图1，图1为本发明实施例方案涉及的硬件运行环境的绝对波阻抗反演成像设备结构示意图。
如图1所示，该绝对波阻抗反演成像设备可以包括：处理器1001，例如中央处理器(centralprocessingunit，cpu)，通信总线1002、用户接口1003，网络接口1004，存储器1005。其中，通信总线1002用于实现这些组件之间的连接通信。用户接口1003可以包括显示屏(display)，可选用户接口1003还可以包括标准的有线接口、无线接口，对于用户接口1003的有线接口在本发明中可为usb接口。网络接口1004可选的可以包括标准的有线接口、无线接口(如无线保真(wireless-fidelity，wi-fi)接口)。存储器1005可以是高速的随机存取存储器(randomaccessmemory，ram)，也可以是稳定的存储器(non-volatilememory，nvm)，例如磁盘存储器。存储器1005可选的还可以是独立于前述处理器1001的存储装置。
本领域技术人员可以理解，图1中示出的结构并不构成对绝对波阻抗反演成像设备的限定，可以包括比图示更多或更少的部件，或者组合某些部件，或者不同的部件布置。
如图1所示，作为一种计算机存储介质的存储器1005中可以包括操作系统、网络通信模块、用户接口模块以及绝对波阻抗反演成像程序。
在图1所示的绝对波阻抗反演成像设备中，网络接口1004主要用于连接后台服务器，与所述后台服务器进行数据通信；用户接口1003主要用于连接用户设备；所述绝对波阻抗反演成像设备通过处理器1001调用存储器1005中存储的绝对波阻抗反演成像程序，并执行本发明实施例提供的绝对波阻抗反演成像方法。
基于上述硬件结构，提出本发明绝对波阻抗反演成像方法的实施例。
参照图2，图2为本发明绝对波阻抗反演成像方法第一实施例的流程示意图，提出本发明绝对波阻抗反演成像方法第一实施例。
在第一实施例中，所述绝对波阻抗反演成像方法包括以下步骤：
s10：采集勘探现场的井控数据，提取所述井控数据中阻抗的低波数部分alow；
s20：获取地震数据，对所述地震数据进行高斯束偏移，获得反射系数r；
s30：将所述反射系数r和所述阻抗的低波数部分alow输入绝对波阻抗反演模型进行约束计算，获得绝对波阻抗z；
s40：通过所述绝对波阻抗z获得波阻抗剖面图像。
在具体实现中，本发明提出了一种融合背景阻抗的高维绝对声阻抗反演成像方法。首先，基于井控数据提供背景速度、密度，获得背景阻抗，即低波数部分。再通过真振幅偏移获得保真反射系数，即阻抗的高波数成分。其具体过程为：先采用gbm实现偏移成像和地下振幅补偿函数点扩散函数(psf)的计算，再使用psf对gbm得到的成像结果进行反褶积，以产生一个宽带、保真的地下反射系数。最后，基于估计出的反射系数和声阻抗之间的关系，提出基于一阶正则化的带约束的优化问题，通过求解约束最优化问题，从稀疏反射系数剖面和井控提供的背景阻抗中估计出绝对阻抗模型。
进一步地，步骤s10中所述提取所述井控数据中阻抗的低波数部分alow，包括：
s101：提取所述井控数据中的背景速度和背景密度；
s102：将所述背景速度与所述背景密度相乘，获得阻抗的低波数部分alow。
进一步地，步骤s20中所述对所述地震数据进行高斯束偏移，获得反射系数r，包括：
s201：对所述地震数据进行常规高斯束偏移成像，获得偏移图像数据；
s202：通过所述背景速度计算获得地下振幅补偿函数点扩散函数psf；
s203：将所述psf和所述偏移图像数据通过反褶积计算，获得反射系数r；
在具体实现中，所述地震数据可以由勘探地点采集获得，也可以有实验室模拟获得；
步骤s20具体过程为，先使用互相关成像条件，实现叠前高斯束偏移：
其中，g(x,xs,ω)和g(x,xr,ω)分别是炮点和检波点的格林函数，d(xr,xs,ω)为在频域中记录的地震数据。
利用式(1)计算的叠前gbm图像可以通过局部波数域中的psf与地下反射系数联系起来：
i(x,k)＝b(x,k)r(x,k)(2)
其中i(x,k)是地下图像，b(x,k)是psf，r(x,k)是反射系数。我们看到，psf作为一个滤波器，模糊反射系数，并将其转换成一个畸变的图像。如果可以获得局部psf，我们可以进行校正：
rb(x,k)＝i(x,k)/b(x,k)(3)
其中rb是校正后的图像。
本发明提出了一种利用gbm计算psf的新方法。方程(2)是局部波数域中的乘法运算，它等价于空间域的卷积：
i(x,x')＝b(x,x')*m(x,x')(4)
其中卷积是在局部坐标x'，如果式(4)的m(x,x')用一个单一的散射体代替，在数学上是脉冲函数，则有i(x,x')＝b(x,x')。也就是说，在背景速度下，通过对单个散射体成像可以直接获得psf。之后引入born模型(wu和aki,1985)来从一个点散射体生成合成数据：
ddelta(xr,xs,ω)＝2ω2∫dmdelta(x)g(x,xs,ω)g(x,xr,ω)dx(5)
其中mdelta(x)是由背景速度扰动和脉冲函数扰动组成的模型。由式(5)得到的合成数据ddelta代入式(1)计算散射点的图像，也就是psfr(x,x')。然后通过局部fft将空间域psf转换到波数域，用式(3)进行校正，然后通过反fft将校正后的成像返回到空间域。
参照图3，图3为本发明绝对波阻抗反演成像方法第二实施例的流程示意图，基于上述图2所示的第一实施例，提出本发明绝对波阻抗反演成像方法的第二实施例。
在第二实施例中，步骤s30中所述将所述反射系数r和所述阻抗的低波数部分alow输入绝对波阻抗反演模型进行约束计算，获得绝对波阻抗z，包括：
s301:根据所述反射系数r和所述阻抗的低波数部分alow通过声阻抗ai剖面计算公式进行循环计算，获得绝对波阻抗z。
在具体实现中，求得所述声阻抗ai剖面计算公式过程为：
一维反射系数可以用(声阻抗)ai表示为：
其中ρi 1,vi 1,zi 1＝ρi 1vi 1分别为第i 1层的密度，速度和阻抗，ri,ρi,vizi＝ρivi分别为第i层的反射系数，密度，速度和阻抗。当反射系数的绝对值较小时，反射系数与声阻抗之间的非线性关系可以重新表述为：
令a＝(lnz)/2，则一维反射系数和a之间的线性关系为：
其中，r是所有散射体的一维反射系数，a为对应的声阻抗的半自然对数。
为了考虑道与道之间的空间关系，对地震剖面数据进行处理，并采用同时求解多个地震道的正演模型。式(8)可推广到二维：
其中r＝(r1,r2,...,rm)代表表示正常入射反射系数成像剖面的矩阵，a＝(a1,a2,...,am)代表ai的半自然对数对应的剖面，
其中，alow为阻抗的低波数部分，l为低通滤波器，r为反射系数，||a||tv为一阶tv范数，a为阻抗的半自然对数，
对于2d模型，一阶tv范数为：
其中
进一步地，参照图5，所述根据所述反射系数r和所述阻抗的低波数部分alow通过声阻抗ai剖面计算公式进行循环计算，获得绝对波阻抗z，包括：
设置计算参数，包括：sk,
其中sk,
设置所述计算参数的初始值，并获取第一公式、第二公式和第三公式；
所述第一公式、所述第二公式和所述第三公式均从所述计算参数的初始值开始计算，令i和k的初始值均为0，
将所述计算参数带入所述第一公式进行计算，获得更新后的ak 1，所述第一公式具体为：
其中，
所述第一公式计算完成后将所述更新后的ak 1带入所述第二公式进行计算，获得更新后的sk,
所述第二公式计算完成后，k值加1，将数值加一后的k、所述更新后的sk,
循环更新a的值n次后，i的值加一，且令ai 1＝an，k的值归零；
将数值加一后的i带入所述第三公式进行计算，获得更新后的
将所述数值加一后的i、归零后的k、所述更新后的
重复更新i的值m次；
总共将所述a的值更新m*n次，获得绝对阻抗的半自然对数a′；
将所述a′带入第四公式进行计算获得绝对波阻抗z，所述第四公式具体为：
z＝exp(2*a′)(15)
进一步地，本发明可以估计出一个更为保真的反射系数，通过融合井控提供的背景阻抗，克服了一维递归反演方法的“挂面条”现象，可以更加高效稳定地进行高维绝对波阻抗反演成像。
为了测试本发明的有效性，选择了sigsbee2a模型进行测试。图6展示了sigabee2a模型理论波阻抗，图7展示了sigabee2a模型一维递归法反演得到的波阻抗，图8展示了sigabee2a模型用本发明提出的方法反演得到的波阻抗。对比图7，图8，可以看出用本发明提出的方法反演得到的波阻抗剖面更加干净，层界面比较清晰。
为了进一步证明该方法的有效性，将图7和图8的结果与图6的结果分别进行了抽道对比。图9展示了真阻抗与一维递归法反演得到的阻抗抽道对比，图10展示了真阻抗与本发明提出的方法反演得到的阻抗抽道对比。对比可以看出，利用本发明提出的方法反演得到的波阻抗与真实阻抗非常地接近。
此外，本发明实施例还提出一种存储介质，所述存储介质上存储有绝对波阻抗反演成像程序，所述绝对波阻抗反演成像程序被处理器执行时实现如上文所述的绝对波阻抗反演成像方法的步骤。
此外，参照图4，本发明实施例还提出绝对波阻抗反演成像装置，所述绝对波阻抗反演成像装置包括：
第一数据提取模块10，用于采集勘探现场的井控数据，提取所述井控数据中阻抗的低波数部分；
第二数据提取模块20，用于获取地震数据，对所述地震数据进行高斯束偏移，获得反射系数r；
计算模块30，用于将所述反射系数r和所述阻抗的低波数部分输入绝对波阻抗反演模型进行约束计算，获得绝对波阻抗z；
图像生成模块40，用于通过所述绝对波阻抗z获得波阻抗剖面图像。
本发明所述绝对波阻抗反演成像装置的其他实施例或具体实现方式可参照上述各方法实施例，此处不再赘述。
需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。
上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。词语第一、第二、以及第三等的使用不表示任何顺序，可将这些词语解释为标识。
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质(如只读存储器镜像(readonlymemoryimage，rom)/随机存取存储器(randomaccessmemory，ram)、磁碟、光盘)中，包括若干指令用以使得一台终端设备(可以是手机，计算机，服务器，空调器，或者网络设备等)执行本发明各个实施例所述的方法。
以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。
1.一种绝对波阻抗反演成像方法，其特征在于，包括以下步骤：
采集勘探现场的井控数据，提取所述井控数据中阻抗的低波数部分alow；
获取地震数据，对所述地震数据进行高斯束偏移，获得反射系数r；
将所述反射系数r和所述阻抗的低波数部分alow输入绝对波阻抗反演模型进行约束计算，获得绝对波阻抗z；
通过所述绝对波阻抗z获得波阻抗剖面图像。
2.根据权利要求1所述的一种绝对波阻抗反演成像方法，其特征在于，所述提取所述井控数据中阻抗的低波数部分alow，包括：
提取所述井控数据中的背景速度和背景密度；
将所述背景速度与所述背景密度相乘，获得阻抗的低波数部分alow。
3.根据权利要求2所述的一种绝对波阻抗反演成像方法，其特征在于，所述对所述地震数据进行高斯束偏移，获得反射系数r，包括：
对所述地震数据进行常规高斯束偏移成像，获得偏移图像数据；
通过所述背景速度计算获得地下振幅补偿函数点扩散函数psf；
将所述psf和所述偏移图像数据通过反褶积计算，获得反射系数r。
4.根据权利要求1所述的一种绝对波阻抗反演成像方法，其特征在于，所述将所述反射系数r和所述阻抗的低波数部分alow输入绝对波阻抗反演模型进行约束计算，获得绝对波阻抗z，包括：
根据所述反射系数r和所述阻抗的低波数部分alow通过声阻抗ai剖面计算公式进行循环计算，获得绝对波阻抗z。
5.根据权利要求4所述的一种绝对波阻抗反演成像方法，其特征在于，所述声阻抗ai剖面计算公式具体为：
其中，alow为阻抗的低波数部分，l为低通滤波器，r为反射系数，||a||tv为一阶tv范数，a为阻抗的半自然对数，
6.根据权利要求5所述的一种绝对波阻抗反演成像方法，其特征在于，所述根据所述反射系数r和所述阻抗的低波数部分alow通过声阻抗ai剖面计算公式进行循环计算，获得绝对波阻抗z，包括：
设置计算参数，包括：sk,
其中sk,
设置所述计算参数的初始值，并获取第一公式、第二公式和第三公式；
所述第一公式、所述第二公式和所述第三公式均从所述计算参数的初始值开始计算，令i和k的初始值均为0，
将所述计算参数带入所述第一公式进行计算，获得更新后的ak 1，所述第一公式具体为：
其中，
所述第一公式计算完成后将所述更新后的ak 1带入所述第二公式进行计算，获得更新后的sk,
所述第二公式计算完成后，k值加1，将数值加一后的k、所述更新后的sk,
循环更新a的值n次后，i的值加一，且令ai 1＝an，k的值归零；
将数值加一后的i带入所述第三公式进行计算，获得更新后的
将所述数值加一后的i、归零后的k、所述更新后的
重复更新i的值m次；
总共将所述a的值更新m*n次，获得绝对阻抗的半自然对数a′；
将所述a′带入第四公式进行计算获得绝对波阻抗z，所述第四公式具体为：
z＝exp(2*a′)。
7.一种绝对波阻抗反演成像装置，其特征在于，所述一种绝对波阻抗反演成像装置包括：
第一数据提取模块，用于采集勘探现场的井控数据，提取所述井控数据中阻抗的低波数部分alow；
第二数据提取模块，用于获取地震数据，对所述地震数据进行高斯束偏移，获得反射系数r；
计算模块，用于将所述反射系数r和所述阻抗的低波数部分alow输入绝对波阻抗反演模型进行约束计算，获得绝对波阻抗z；
图像生成模块，用于通过所述绝对波阻抗z获得波阻抗剖面图像。
8.一种绝对波阻抗反演成像设备，其特征在于，包括：存储器、处理器及存储在所述存储器上并可在所述处理器上运行的绝对波阻抗反演成像程序，所述绝对波阻抗反演成像程序被所述处理器执行时实现如权利要求1至6中任一项所述的一种绝对波阻抗反演成像方法的步骤。
9.一种存储介质，其特征在于，所述存储介质上存储有绝对波阻抗反演成像程序，所述绝对波阻抗反演成像程序被处理器执行时实现如权利要求1至6中任一项所述的一种绝对波阻抗反演成像方法的步骤。
技术总结