一种基于数据相似性的地震全波形反演速度建模方法

1.本发明涉及勘探地球物理技术领域，尤其涉及一种基于数据相似性的地震全波形反演速度建模方法。


背景技术：

2.全波形反演是一种非线性高分辨率速度建模方法，其通过求解最优化问题获得地下岩石速度参数信息。与传统的速度分析和射线类层析速度建模技术相比，全波形反演方法不仅能获得低波数背景速度信息，而且还能恢复中、高波数速度信息，最终得到全波数段高分辨率地下速度模型。其反演结果既可以用于后续的地震偏移成像，又可以直接用于地震解释来确定油气藏位置和内部精细构造。然而，在实际应用中，受数据采集频率、接收孔径和环境噪音等因素的影响，全波形反演的成功应用较为困难。特别是对于我国陆地低信噪比资料而言，实施全波形反演速度建模具有极大的挑战性。
3.在实际应用中，全波形反演面临两个主要问题：一个是受最大偏移距和数据有限频带以及初始模型精度低的限制，在数据拟合时容易出现周波跳跃问题，致使反演过程陷入局部极值，最终获得错误的反演结果；二是在观测数据信噪比较低的情况下，反演过程会拟合无关噪声，这些噪声最终会被投射到反演模型中，导致错误的结果和解释。对于周波跳跃问题，常用的解决方法是在野外采集时使用宽频检波器，记录大量低频有效信号，从较精确的初始模型(一般由层析成像获得)开始，然后进行分频段多尺度反演，从低频到高频逐步提高分辨率。然而，这种方法对野外施工要求较为严格，需要较大的偏移距和能够记录低频信号的高性能检波器。常规地震采集的数据往往缺失低频信息，而高性能低频检波器的价格一般是常规检波器的数十倍。因此，分频段多尺度全波形反演方法往往不能直接应用于低频缺失的常规采集数据，从而大大限制了其应用范围。对于数据信噪比低的问题，常用的解决办法是进行数据处理来减少非地震因素的干扰，主要包括对面波、线型噪声和背景噪声等进行滤波去除。但是，在进行降噪处理时需要在信号的保真度和抑制效果之间进行权衡，以便获得更好的反演结果。处理过度会损害有效信号，从而影响反演结果的准确性；而处理不足时残留的噪声仍会对反演结果造成干扰，可能造成错误的解释。
4.因此，亟需研发一种稳定实用的地震全波形反演速度建模方法，以解决其当前在实际资料应用中的难题。


技术实现要素：

5.为解决上述技术问题，本发明公开了一种基于数据相似性的地震全波形反演速度建模方法。本发明能够解决在地震资料信噪比低、初始模型不准确的情况下压制无关噪声干扰和缓解周波跳跃问题，促进全波形反演的成功实施，以获得准确可靠的反演结果；同时能够缩减野外采集的成本，简化反演步骤，节省大量时间和计算成本。
6.为实现上述目的，本发明采用下述技术方案：
7.一种基于数据相似性的地震全波形反演速度建模方法，包括如下步骤：
8.步骤s1、获取输入数据；
9.步骤s2、基于输入数据，使用地震子波和初始速度模型计算与观测数据相对应的合成数据u(xr,t)；
10.步骤s3、用目标泛函f(v)替换常规全波形反演的伴随源；
11.步骤s4、计算目标泛函f(v)中的相关系数；
12.步骤s5、计算目标泛函f(v)中对应的伴随源和梯度；
13.步骤s6、对梯度进行对角汉森预条件和平滑算子处理；
14.步骤s7、使用共轭梯度寻优算法进行寻优，得到最优解，并将更新的模型作为下一次反演的初始模型v1(x)；
15.步骤s8、判断是否达到精度要求或达到指定的迭代次数，若满足要求，则结束全波形反演过程；若不满足，依次重复步骤s2-s7，直到满足要求。
16.可选地，步骤s1中，所述输入数据包括：初始纵波速度模型v0(x)、地震子波、观测系统信息和野外低信噪比的观测数据d(xr,t)，其中，x代表反演模型空间的位置，xr为检波点的空间位置。
17.可选地，步骤s3中，目标函数f(v)公式如下：
[0018][0019]
式中，n为所有数据的总采样点数，s(xr,t)为由局部合成和观测数据计算的相关系数，取值范围为[-1,1]。
[0020]
可选地，步骤s4中，计算目标泛函f(v)中的相关系数的步骤，具体包括：
[0021]
步骤s41、利用二维/三维高斯窗g(kr,h)获取局部地震数据，公式如下：
[0022][0023]
式中，d(xr,t,kr,h)和u(xr,t,kr,h)分别代表局部窗口的观测与合成数据，kr和h分别代表在检波点空间位置和时间轴方向上的窗口大小；
[0024]
步骤s42、为将局部尺度数据恢复成真实尺度大小，定义相应逆变换为：
[0025][0026]
式中，∫∫∫d2krdh表示对局部窗口内的数据进行积分求和，此式表示将所有局部观测数据d(xr,t,kr,h)和局部合成数据u(xr,t,kr,h)分别与高斯窗g(kr,h)相乘之后进行积分，使局部数据恢复成全局数据大小；
[0027]
步骤s43、将局部数据去除均值，得到和公式如下：
[0028][0029]
[0030]
式中，n为局部窗口的总采样点数；
[0031]
步骤s44、求解相关系数s(xr,t)，由去除均值的局部数据求取归一化互相关得到，公式如下：
[0032][0033]
步骤s45、计算完相关系数，进行求和取均值，最后与负一相乘，得到目标泛函f(v)。
[0034]
可选地，步骤s5中，计算目标泛函f(v)中对应的伴随源和梯度的步骤，具体包括：
[0035]
步骤s51、通过对目标泛函f(v)求取相对于模型速度参数的导数，得到梯度g和相应的伴随源δf(xr,t)，公式如下：
[0036][0037]
式中，表示目标泛函相对于模型速度参数的导数，表示合成数据相对于模型速度参数的导数；
[0038]
步骤s52、伴随源δf(xr,t)，公式如下：
[0039][0040]
式中，为局部尺度伴随源a(xr,t,kr,h)去除均值后的结果，此式表示将去除均值的局部尺度伴随源与高斯窗g(kr,h)相乘之后进行积分，使其恢复成原始全局数据的大小；
[0041]
步骤s53、的计算方式如下：
[0042][0043]
步骤s54、局部尺度伴随源a(xr,t,kr,h)为：
[0044][0045]
步骤s55、通过反传伴随源δf(xr,t)得到反传波场，反传波场再与正传波场进行互相关得到梯度g。
[0046]
本发明的有益效果是，
[0047]
本发明提出的一种基于数据相似性的地震全波形反演速度建模方法，利用一个二维窗口来获得局部尺度的地震数据，然后应用零均值归一化互相关函数来测量局部合成数据和观测数据的相似性，并将其作为目标泛函来迭代更新地下速度模型；新的目标泛函增强了相位拟合的贡献，同时减少了振幅的约束，这有助于提高全波形反演对不准确初始模型的容忍度，缓解周波跳跃问题；通过观测数据和合成数据进行互相关，可以减轻噪声对反
演结果的影响；此外，通过计算归一化的局部波形相似性，可以自适应地平衡伴随源的强近偏移反射和弱远偏移折射的振幅，从而改善低波数更新，并在深层建模方面有着巨大的潜力；再结合分频多尺度反演策略进行迭代更新，最终可以获得高分辨率反演结果。本发明可以在地震资料信噪比低、初始模型不准确的情况下压制无关噪声干扰和缓解周波跳跃问题，促进全波形反演的成功实施，以获得准确可靠的反演结果；另一方面也缩减了野外采集的成本，简化了反演步骤，节省大量时间和计算成本。
附图说明
[0048]
图1为本发明一种基于数据相似性的地震全波形反演速度建模方法的流程示意图；
[0049]
图2为本发明一实施例示出的salt模型的真实速度场；
[0050]
图3为本发明一实施例示出的全波形反演输入的salt模型初始速度场；
[0051]
图4为本发明一实施例示出的一个观测数据炮集；
[0052]
图5为本发明一实施例示出的一个合成数据炮集；
[0053]
图6为本发明一实施例示出的由单个观测与合成数据炮集计算得到相关系数；
[0054]
图7为传统全波形反演第一次迭代伴随源(波形残差)；
[0055]
图8为本发明一实施例示出的第一次迭代伴随源；
[0056]
图9为常规计算所得的梯度；
[0057]
图10为本发明一实施例示出的进行对角汉森预条件和平滑算子处理后的梯度；
[0058]
图11为使用常规全波形反演方法获得的salt模型最终速度结果；
[0059]
图12为本发明一实施例示出的salt模型最终速度结果。
具体实施方式
[0060]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0061]
一种基于数据相似性的地震全波形反演速度建模方法，如图1所示，包括如下步骤：
[0062]
步骤s1、获取输入数据；
[0063]
所述输入数据包括：初始纵波速度模型v0(x)、地震子波、观测系统信息和野外低信噪比的观测数据d(xr,t)，其中，x代表反演模型空间的位置，xr为检波点的空间位置；
[0064]
步骤s2、基于输入数据，使用地震子波和初始速度模型计算与观测数据相对应的合成数据u(xr,t)；
[0065]
步骤s3、用目标泛函f(v)替换常规全波形反演的伴随源；
[0066]
目标函数f(v)公式如下：
[0067][0068]
式中，n为所有数据的总采样点数，s(xr,t)为由局部合成和观测数据计算的相关
系数，取值范围为[-1,1]；
[0069]
步骤s4、计算目标泛函f(v)中的相关系数；
[0070]
计算目标泛函f(v)中的相关系数的步骤，具体包括：
[0071]
步骤s41、利用二维/三维高斯窗g(kr,h)获取局部地震数据，公式如下：
[0072][0073]
式中，d(xr,t,kr,h)和u(xr,t,kr,h)分别代表局部窗口的观测与合成数据，kr和h分别代表在检波点空间位置和时间轴方向上的窗口大小；
[0074]
步骤s42、为将局部尺度数据恢复成真实尺度大小，定义相应逆变换为：
[0075][0076]
式中，∫∫∫d2krdh表示对局部窗口内的数据进行积分求和，此式表示将所有局部观测数据d(xr,t,kr,h)和局部合成数据u(xr,t,kr,h)分别与高斯窗g(kr,h)相乘之后进行积分，使局部数据恢复成全局数据大小；
[0077]
步骤s43、为了将数据的中心点移动到原点附近，以减少数据中的噪声和冗余信息，提高数据的可解释性，将局部数据去除均值，得到和公式如下：
[0078][0079][0080]
式中，n为局部窗口的总采样点数；
[0081]
步骤s44、求解相关系数s(xr,t)，由去除均值的局部数据求取归一化互相关得到，公式如下：
[0082][0083]
步骤s45、计算完相关系数，进行求和取均值，最后与负一相乘，得到目标泛函f(v)；
[0084]
步骤s5、计算目标泛函f(v)中对应的伴随源和梯度；
[0085]
计算目标泛函f(v)中对应的伴随源和梯度的步骤，具体包括：
[0086]
步骤s51、通过对目标泛函f(v)求取相对于模型速度参数的导数，得到梯度g和相应的伴随源δf(xr,t)，公式如下：
[0087][0088]
式中，表示目标泛函相对于模型速度参数的导数，表示合成数据相对
于模型速度参数的导数；
[0089]
步骤s52、伴随源δf(xr,t)，公式如下：
[0090][0091]
式中，为局部尺度伴随源a(xr,t,kr,h)去除均值后的结果，此式表示将去除均值的局部尺度伴随源与高斯窗g(kr,h)相乘之后进行积分，使其恢复成原始全局数据的大小；
[0092]
步骤s53、的计算方式如下：
[0093][0094]
步骤s54、局部尺度伴随源a(xr,t,kr,h)为：
[0095][0096]
步骤s55、通过反传伴随源δf(xr,t)得到反传波场，反传波场再与正传波场进行互相关得到梯度g；
[0097]
步骤s6、对梯度进行对角汉森预条件和平滑算子处理，用于增加梯度的平滑和深部照明，这样处理可以缓解周波跳跃问题，并且还可增加反演深度，具体而言，对角汉森预条件可以优化梯度的方向和增加深部照明，从而加快算法的收敛速度，而平滑算子处理可以减少梯度的高频分量，降低算法对噪声的敏感性，缓解周波跳跃问题以增加反演的稳定性和准确性；
[0098]
步骤s7、使用共轭梯度寻优算法进行寻优，得到最优解，并将更新的模型作为下一次反演的初始模型v1(x)；
[0099]
步骤s8、判断是否达到精度要求或达到指定的迭代次数，若满足要求，则结束全波形反演过程；若不满足，依次重复步骤s2-s7，直到满足要求。
[0100]
将本发明所提供的方法应用于seg/eage salt模型的速度反演建模中，取得了较好的反演结果。真实salt速度模型如图2所示，全波形反演输入的初始模型如图3所示，震源函数为雷克子波，观测数据集为34个炮集，均匀分布在地表，每个炮集包含645个检波点，间隔为12.5m。图4为其中一个观测数据炮集，图5为由初始模型获得的其中一个合成数据炮集。图6是本发明方法中由单个观测与合成数据炮集计算得到相关系数算子，迭代的目的是使所有相关系数算子的和最大化。图7为传统全波形反演第一次迭代中的伴随源，图8为本发明专利方法第一次迭代中的伴随源。新的伴随源不再是将合成数据与观测数据进行简单的波形相减，而是根据两者的相似性来自适应地平衡波形，优先拟合早至波，即先拟合模型的低波数成分，然后随着更新迭代不断自动添加高波数成分进行拟合，这样有效避免了周波跳跃问题，从而实现稳健的反演。图9为常规计算所得的梯度，图10为进行对角汉森预条件和平滑算子处理后的梯度，可以看出，经过处理后梯度的高频分量得到压制，整体变得平滑，而且增强了深部照明，这使得反演的收敛速度、稳定性和准确性都得到了显著提高。图
11为常规全波形反演方法获得的salt模型的最终速度结果，可以看到模型的浅部出现了很多点状、块状的更新，这是周波跳跃对反演结果造成的不良影响，这些虚假的结构会导致错误的解释。图12是使用本发明方法获得的salt模型最终速度结果，该模型的宏观速度结构和细节都得到了很好的恢复，得到了低、中、高全波数段高分辨率速度模型信息。值得注意的是，该模型的盐体轮廓特别是下边界和盐下结构，均得到了较好地恢复。
[0101]
本发明提出的一种基于数据相似性的地震全波形反演速度建模方法，主要解决全波形反演中由于初始模型精度低而引起的周期跳跃问题，避免反演算法陷入局部极值，同时可以减少无关噪声的干扰，自适应地平衡波形，实现自动从低波数向高波数拟合更新，使得反演算法更加稳健可靠，最终可以获得准确、高分辨率的反演结果。在实际应用中，全波形反演建模的成功往往需要大量高信噪比的低频(小于3hz)或超大偏移距(大于8km)的地震数据，而且需要高精度的初始模型。然而，低频检波器的费用通常是常规检波器价格的数十倍；大偏移距采集对地形有很大的要求，这在很多情况下都难以满足；对低信噪比的地震数据进行去噪处理时，为了不影响有效信号，处理的结果中会有残余的噪声，这仍会对反演造成干扰；高精度的初始模型一般需要用其他方法如层析成像进一步优化初始模型，这也增加了成本。为解决这些问题，本发明使用一个二维窗口来获得局部尺度的地震数据，然后应用零均值归一化互相关函数来测量局部合成数据和观测数据的相似性，并将其作为目标泛函来迭代更新地下速度模型；再将此方法与分频多尺度反演策略结合，可以进一步稳健反演算法并最终获得高分辨率反演结果。本发明利用地震数据的局部相似性实现了从低频到高频，从低波数到高波数对速度模型进行更新，从而解决了因初始模型精度低带来的周波跳跃问题，避免反演过程陷入局部极值，并可减少无关噪声的干扰，获得准确可靠的反演结果；另一方面降低了野外数据采集对高信噪比的低频数据和大偏移距数据的要求，缩减了野外采集的成本，而且此发明还简化了反演步骤，提高了反演的收敛速度和精度，节省了大量时间和计算成本。
[0102]
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

技术特征：
1.一种基于数据相似性的地震全波形反演速度建模方法，其特征在于，包括如下步骤：步骤s1、获取输入数据；步骤s2、基于输入数据，使用地震子波和初始速度模型计算与观测数据相对应的合成数据u(x
r
,t)；步骤s3、用目标泛函f(v)替换常规全波形反演的伴随源；步骤s4、计算目标泛函f(v)中的相关系数；步骤s5、计算目标泛函f(v)中对应的伴随源和梯度；步骤s6、对梯度进行对角汉森预条件和平滑算子处理；步骤s7、使用共轭梯度寻优算法进行寻优，得到最优解，并将更新的模型作为下一次反演的初始模型v1(x)；步骤s8、判断是否达到精度要求或达到指定的迭代次数，若满足要求，则结束全波形反演过程；若不满足，依次重复步骤s2-s7，直到满足要求。2.如权利要求1所述的一种基于数据相似性的地震全波形反演速度建模方法，其特征在于，步骤s1中，所述输入数据包括：初始纵波速度模型v0(x)、地震子波、观测系统信息和野外低信噪比的观测数据d(x
r
,t)，其中，x代表反演模型空间的位置，x
r
为检波点的空间位置。3.如权利要求1所述的一种基于数据相似性的地震全波形反演速度建模方法，其特征在于，步骤s3中，目标函数f(v)公式如下：式中，n为所有数据的总采样点数，s(x
r
,t)为由局部合成和观测数据计算的相关系数，取值范围为[-1,1]。4.如权利要求1所述的一种基于数据相似性的地震全波形反演速度建模方法，其特征在于，步骤s4中，计算目标泛函f(v)中的相关系数的步骤，具体包括：步骤s41、利用二维/三维高斯窗g(k
r
,h)获取局部地震数据，公式如下：式中，d(x
r
,t,k
r
,h)和u(x
r
,t,k
r
,h)分别代表局部窗口的观测与合成数据，u(x
r
,t)代表步骤s2中计算所得的合成数据，k
r
和h分别代表在检波点空间位置和时间轴方向上的窗口大小；步骤s42、为将局部尺度数据恢复成真实尺度大小，定义相应逆变换为：式中，∫∫∫d2k
r
dh表示对局部窗口内的数据进行积分求和，此式表示将所有局部观测数据d(x
r
,t,k
r
,h)和局部合成数据u(x
r
,t,k
r
,h)分别与高斯窗g(k
r
,h)相乘之后进行积分，使局部数据恢复成全局数据大小；
步骤s43、将局部数据去除均值，得到和公式如下：式中，n为局部窗口的总采样点数；步骤s44、求解相关系数s(x
r
,t)，由去除均值的局部数据求取归一化互相关得到，公式如下：步骤s45、计算相关系数，进行求和取均值，最后与负一相乘，得到目标泛函f(v)。5.如权利要求1所述的一种基于数据相似性的地震全波形反演速度建模方法，其特征在于，步骤s5中，计算目标泛函f(v)中对应的伴随源和梯度的步骤，具体包括：步骤s51、通过对目标泛函f(v)求取相对于模型速度参数的导数，得到梯度g和相应的伴随源δf(x
r
,t)，公式如下：式中，表示目标泛函相对于模型速度参数的导数，表示合成数据相对于模型速度参数的导数；步骤s52、伴随源δf(x
r
,t)，公式如下：式中，为局部尺度伴随源a(x
r
,t,k
r
,h)去除均值后的结果，此式表示将去除均值的局部尺度伴随源与高斯窗g(k
r
,h)相乘之后进行积分，使其恢复成原始全局数据的大小；步骤s53、的计算方式如下：步骤s54、局部尺度伴随源a(x
r
,t,k
r
,h)为：步骤s55、通过反传伴随源δf(x
r
,t)得到反传波场，反传波场再与正传波场进行互相关得到梯度g。

技术总结
本发明公开了一种基于数据相似性的地震全波形反演速度建模方法，涉及勘探地球物理技术领域，包括：获取输入数据；计算与观测数据相对应的合成数据；用目标泛函f(v)替换常规全波形反演的伴随源；计算f(v)的相关系数；计算f(v)中对应的伴随源和梯度；对梯度进行对角汉森预条件和平滑算子处理；寻优，并将更新的模型作为下一次反演的初始模型；判断是否达到精度要求或达到指定的迭代次数，若不满足，依次重复上述步骤。本发明使获得局部尺度的地震数据，应用零均值归一化互相关函数来测量局部合成数据和观测数据的相似性，并将其作为目标泛函来迭代更新地下速度模型；再将此与分频多尺度反演策略结合，稳健反演算法并最终获得高分辨率反演结果。辨率反演结果。辨率反演结果。


技术研发人员：杨继东 于由财 黄建平 秦善源 张浩
受保护的技术使用者：中国石油大学（华东）
技术研发日：2023.05.06
技术公布日：2023/9/7