基于波场局部相关时移的时间域全波形反演方法

摘要

本发明涉及一种基于波场局部相关时移的时间域全波形反演方法，通过滑动时窗时移矫正将模拟及观测记录用时窗截断，根据时窗内模拟与观测记录互相关最大值的位置对模拟记录做时移，使互相关程度增大，在沿着采样时间轴移动窗口对下一个时窗内的波形做相关时移直到时窗移动覆盖所有采样点，即完成对所有模拟记录的波形时移校正，使模拟与观测记录匹配程度提高。由于对模拟记录的波形进行了位移，使各个采样点的相位信息更加接近观测记录，但对应的振幅信息并没有校正，为减小振幅错误对反演的影响，在进行全波形反演时采用对振幅信息依赖较小的全局互相关目标函数。本发明在不降低计算效率的情况下，减少全波形反演中跳周问题的发生，提高反演精度。

说明书

技术领域

本发明属于地震勘探技术领域，涉及一种地震勘探成像方法，具体涉及一种时间域全波形反演方法，特别涉及一种基于波场局部相关时移的时间域全波形反演方法。

背景技术

石油勘探开发领域目前已经由浅层油气藏转向深部复杂构造地质区域的油气藏。传统的地震勘探方法由于其勘探深度较浅，施工方法落后以及勘探精度较低，已经不能满足当今对石油能源勘探开发的需要。因此，许多新兴的工业技术和勘探方法逐渐发展开来。全波形反演作为一种高精度地下介质参数反演方法，近年来成为地球物理勘探领域的热门课题。它是通过对待重建的初始模型正演得到的模拟记录与实际勘探得到的观测记录进行匹配，得到模型参数更新量再通过最优化算法迭代更新最终得到地下介质参数的真实分布。全波形反演方法与传统地震反演成像方法最大的区别在于：其利用了叠前地震波场的全波信息包括直达波、折射波、反射波以及多次波甚至噪声等直接进行反演计算，因此包含了丰富的地下参数信息，是目前反演精度最高的方法。

20世纪80年代，Lailly首次提出全波形反演的概念，并在时间域尝试了全波形反演。但是由于反演梯度的计算中需要直接计算Jacobi矩阵，即得到地下介质每一点的梯度值都需要进行一次正演模拟，而一个模型中有成千上万个点，因此计算量是非常庞大的，这也极大地限制了全波形反演的发展。之后，Trantola提出伴随状态法求取梯度，使每次梯度的求取只需要两次正演计算，该方法的提出极大地推动了全波形反演的发展。20世纪90年代，Pratt在频率域实现了全波形反演，只需要从低到高几个离散的频率就可以实现全波形反演，使计算效率进一步提高。

全波形反演方法面临的主要问题就是跳周问题，即模拟记录和观测记录波形相位相差半个周期以上，导致波形匹配错误，从而导致反演结果出错。产生跳周的原因主要有两方面，一是反演的初始模型与真实地下模型相差太远，导致模拟记录波形与观测记录波形相差太大，从而导致跳周。1995年Bunks提出多尺度策略，由于低频波形波长较长，不容易跳周，所以先用模拟记录和观测记录中的低频信息反演出模型的大尺度构造，然后再以低频段反演结果为初始模型进行高频段反演，最终得到高精度地下介质参数分布。产生跳周的第二个原因是缺少低频信息，地震记录中的低频信号携带了大尺度构造的信息且不易跳周，然而实际观测记录中往往缺失低频信息，这时多尺度策略也无法将地下介质参数准确的反演出来。许多学者为了解决低频缺失导致的跳周问题做了很多研究，解决方法大致可以分为三类。

第一类是用走时层析成像方法为全波形反演建立更加精确的初始模型。由于传统基于最小二乘目标函数的全波形反演方法非常依赖于振幅信息，而缺失低频信息导致不能为反演提供一个好的初始模型，从而导致模拟记录与观测记录的振幅残差出现错误，导致反演失败。而基于走时的层析成像方法不需要振幅信息，而且走时信息对模型参数变化不敏感，可以为全波形反演提供一个好的初始模型。第二类是采用波场延拓类方法例如偏移速度分析为全波形反演提供一个较好的初始模型。第三类是对观测波场和模拟波场进行一些运算和变换产生人造低频信息，这些产生的低频信息与缺失的真实低频信息的振幅相位大致吻合，可以代替缺失的低频信息应用到全波形反演中。如Wu(2014)和Chi(2014)提出了基于希尔伯特变换的包络全波形反演方法，Hu(2014)提出了基于拍音理论的低频信息重构方法等。上述方法要么依赖于其他技术手段的辅助要么计算稳定性较差，在计算效率或者反演结果上都有一些不足之处。

发明内容

本发明的目的就在于针对上述现有技术的不足，提供一种解决跳周问题并且不增加额外计算量的基于波场局部相关时移的时间域全波形反演方法。

本发明的目是通过以下技术方案实现的：

本发明基于波场局部相关时移的时间域全波形反演方法，核心是提出一种滑动时窗时移矫正的方法，使模拟记录与观测记录之间匹配程度增大，达到减小跳周的效果。由于低频信息的缺失，在初始模型不准确的情况下，中高频信息在波形匹配的时候很容易产生跳周。采用时窗法，将模拟记录和观测记录用时窗截断，根据时窗内模拟记录与观测记录互相关最大值的位置对模拟记录做时移，使模拟记录和观测记录的互相关程度增大，然后沿着采样时间轴移动窗口对下一个时窗内的波形做相关时移直到时窗移动覆盖到所有采样点为止，即完成对所有模拟记录的波形时移校正，使模拟记录与观测记录的匹配程度提高。同时，由于对模拟记录的波形进行了位移，使各个采样点的相位信息更加接近观测记录，但是对应的振幅信息并没有校正，因此为了减小振幅错误对反演的影响，在进行全波形反演时采用对振幅信息依赖较小的全局互相关目标函数。

一种基于波场局部相关时移的时间域全波形反演方法，包括以下步骤：

a、对实际地震观测记录进行子波估计、低频保护去噪、缺失地震道补偿多次波衰减、面波切除、消除交混回响等预处理；

b、首先在预估速度范围建立线性递增初始模型，根据要求设定时间域全波形反演相关参数，包括地震子波主频f，低通滤波截断频率fluc，模型大小nz×nx，网格距dx,dz，采样总时间T，时间采样间隔dt，每个频段最大迭代次数iter max，最优化算法的迭代步长q，目标函数要求精度tol，模型速度估计的最大值vmax与最小值vmin；

c、用子波在初始模型上进行正演，得到模拟记录。对模拟记录和观测记录做低通滤波处理，得到低频段信号；

d、设置滑动时窗的长度l和窗口每次沿采样时间轴移动的距离s。对观测记录和模拟记录逐道做滑动时窗互相关计算，得到一个时窗内的互相关系数：

C_rw为互相关系数，t为时窗内的采样点，u为模拟记录，d为观测记录，x为检波器位置,取互相关系数最大时的τ值并对一个时窗内的模拟记录做时移可以得到：

表示一个时窗内做完时移后的模拟记录，ns表示时窗沿采样时间轴最大移动次数，将所有地震道做完滑动时窗时移校正的模拟记录记为u^*；

e、根据全局互相关原理建立目标函数：

J为目标函数，v为地下介质速度参数，令：

对目标函数两端对速度求导数可得梯度表达式为：

P_f为时间域正传波场,P_b为反传波场；

f、利用L-BFGS优化算法对速度模型进行迭代更新，先反演出模型的大尺度构造，再以低频段反演结果做为初始模型进行全频段反演，将模型的细节构造反演出来，最终得高精度的地下模型。

与现有技术相比，本发明的有益效果在于：本发明提出了一种基于波场局部相关时移的时间域全波形反演方法，采用时窗对模拟记录和观测记录进行截断,计算得到时窗内模拟记录与观测记录互相关的最大值，然后对模拟记录做时移校正使之与观测记录匹配程度提高。将时窗沿采样时间轴移动,逐个窗口对模拟记录做时移校正，最终使整个模拟记录与观测记录的匹配程度提高。由于设置了时窗沿采样时间轴移动的步长，只需要将时窗移动有限次就可将整个观测记录校正完毕，在计算效率上相比于传统时间域全波形反演方法并没有降低。本发明利用滑动时窗方法对模拟记录进行了时移校正，解决了以下问题：

1、利用滑动时窗对模拟记录和观测记录进行时移校正，使整个观测记录与模拟记录的匹配程度提高，消除了原始波形匹配时波形相位差大于半个周期的部分，有效的较少了跳周的发生。

2、采用了全局互相关目标函数减小反演对振幅信息的依赖。由于对模拟记录的波形进行了位移，使各个采样点的相位信息更加接近观测记录，但是对应的振幅信息并没有校正，因此采用全局互相关目标函数减小振幅错误对反演的影响，同时互相关计算使远偏移距记录的权重增加，而远偏移距记录往往携带了大尺度构造的信息，因此能进一步加强对模型大尺度构造的反演精度，为后续全频带反演提供更加精确的初始模型。

3、通过设置时窗沿采样时间轴移动的步长，使时窗通过有限次的移动就能将整个模拟记录校正完毕，每个时窗内只进行互相关运算和时移两步操作，因此相比于传统时间域全波形反演几乎没有计算效率上的下降，在不损失计算效率的同时提高反演精度。

基于波场局部相关时移的时间域全波形反演方法原理简单易于编程实现，在保持计算效率不变的情况下有效的减少了跳周的发生，非常适用于实际地震勘探，该方法能使地震勘探反演效果显著提高。

附图说明

图1基于波场局部相关时移的时间域全波形反演方法流程图；

图2a 20Hz主频的雷克子波；

图2b雷克子波频谱；

图3a去除低频成分11Hz以下信息的雷克子波；

图3b去除除低频成分11Hz以下信息的雷克子波频谱；

图4真实模型图；

图5线性递增初始模型图；

图6a原始模拟记录与观测记录各地震道归一化互相关值(l＝50,s＝20)；

图6b基于本发明方法处理后的模拟记录与观测记录各地震道归一化互相关值(l＝50,s＝20)；

图7a基于传统方法的20Hz低通滤波时间域全波形反演结果图；

图7b基于本发明方法的20Hz低通滤波时间域全波形反演结果图；

图8a以图7a为初始模型进行基于传统方法的全频带全波形反演结果图；

图8b以图7b为初始模型进行基于本发明方法的全频带全波形反演结果图；

图9a图8b反演结果60道速度对比图；

图9b图8b反演结果150道速度对比图。

具体实施方式

下面结合附图和实例对本发明进一步的详细说明

如图1所示，基于波场局部相关时移的时间域全波形反演方法，用滑动时窗的方法提高地震观测记录与模拟记录之间的匹配程度，较少跳周的发生，使反演结果更加精确，包括以下步骤：

a、对实际地震观测记录进行子波估计、低频保护去噪、缺失地震道补偿多次波衰减、面波切除、消除交混回响等预处理；

b、首先在预估速度范围建立线性递增初始模型，根据要求设定时间域全波形反演相关参数，包括地震子波主频f，低通滤波截断频率fluc，模型大小nz×nx，网格距dx,dz，采样总时间T，时间采样间隔dt，每个频段最大迭代次数iter max，最优化算法的迭代步长q，目标函数要求精度tol，模型速度估计的最大值vmax与最小值vmin；

c、用子波在初始模型上进行正演，得到模拟记录。对模拟记录和观测记录做低通滤波处理，得到低频段信号；

d、设置滑动时窗的长度l和窗口每次沿采样时间轴移动的距离s。对观测记录和模拟记录逐道做滑动时窗互相关计算，得到一个时窗内的互相关系数：

C_rw为互相关系数，t为时窗内的采样点，u为模拟记录，d为观测记录，x为检波器位置,取互相关系数最大时的τ值并对一个时窗内的模拟记录做时移可以得到：

表示一个时窗内做完时移后的模拟记录，ns表示时窗沿采样时间轴最大移动次数。将所有地震道做完滑动时窗时移校正的模拟记录记为u*；

e、根据全局互相关理论建立目标函数：

J为目标函数，v为地下介质速度参数。令：

对目标函数两端对速度求导数可得梯度表达式为：

简化后可表示为：

其中，λ为伴随源，其表达式为：

通过链式法则可以得到基于伴随状态法的时间域全波形反演梯度公式：

其中，P_f为时间域正传波场，L^-1λ表示反传波场，其中L^-1为反传算子；

f、反演基于L-BFGS优化算法进行迭代更新，并通过Wolfe收敛准则寻找合适的步长；其迭代公式表示为：

m_k+1＝m_k-α_kH_kg_k

其中，H_k为近似Hessian矩阵的逆矩阵，m_k为模型更新参数，g_k为模型更新梯度，α_k为Wolfe线性搜索得到的步长，k表示迭代次数。

L-BFGS优化算法在迭代计算过程中需要保存的矩阵个数很少，对计算机内存要求小，同时，对用于更新Hessian矩阵，其迭代公式如下：

H_k+1＝V_k^TH_kV_k+ρ_ks_ks_k^T

s_k＝m_k+1-m_k，y_k＝g_k+1-g_k

其中，V_k为计算过程的中间变量，s_k为计算过程的中间变量，ρ_k为中间变量，I表示单位矩阵，H_k+1是根据向量对{s_k,y_k}和H_k计算得到，只储存n个向量对来隐式表达Hessian矩阵的逆矩阵。H_kgk的乘积可以通过梯度g_k与向量对{s_k,y_k}之间一系列向量的内积与向量的和来获得的。在每一次更新完成后，上一步向量对将被当前新向量对{s_k+1,y_k+1}取代。因此，向量对集合中包含最近n次迭代的曲率信号。在时实际中，当n≥5时都能获得较满意的结果。L-BFGS优化算法的初始近似Hessian矩阵每一次迭代中都不同。近似Hessian矩阵的逆矩阵H_k需满足以下更新公式：

模型的更新方向，通过以下方法实现：

(1)令则q＝q-α_iy_i，其中，i＝k-1,k-2,…k-n；α_i为计算过程中的中间变量；

(2)令则r＝r-s_i(α_i-β)，其中，i＝k-n,k-n+1,…,k-1；为初始近似Hessian矩阵；

(3)通过上述步骤得到更新量H_kg_k＝r。

通过上述方法求得模型的更新量H_kg_k，然后再通过Wolfe线性搜索获得合适的步长α_k进行更新迭代。L-BFGS优化算法有效的克服了显式求取Hessian矩阵及其逆的困难，其具有超线性收敛速度，计算效率高，占用内存小，精度高等优点，较适合求解大规模非线性优化问题。

最后，判断反演结果是否满足设置的终止条件，即反演结果与真实模型相差很小，若果是则停止迭代，如果不是，则继续计算直到满足终止条件为止，结束计算。

本发明的程序是在MATLAB2016b软件框架下编写完成，根据相应的并行计算要求，搭建MATLAB并行工作库的安装环境，并安装MATLAB并行计算工具箱(Parallel ComputingToolbox)。

实施例1

根据勘探要求，将Parallel Computing Toolbox和MATLAB DistributedComputing Server(R2016b)在Windows 10专业版系统下进行安装，进行MATLAB并行平台的搭建。

利用Marmousi进行缺失低频信息全波形反演测试。真实模型(图4)和线性递增的初始模型(图5)。

表1缺失低频信息全波形反演测试参数

网格大小网格距测线长度纵向深度fluc频带宽度69×19212.5m2400m862.5m20Hz11～20Hz

模型网格大小为69×192，网格距12.5，横向距离为2400m，纵向深度为862.5m，模型中地震波速度范围从1500m/s到4000m/s，地震检波器安置在模型表面，每一个网格点都是一个检波器，且检波器之间的间隔为12.5m，一共激发12个震源。震源选用缺失11Hz以下信息20Hz主频的雷克子波(图3a)，采样间隔为0.001s，实际采样总长度为2s，频率范围从11Hz到20Hz。滑动时窗的窗长l取50个采样点，步长s取20个采样点。

缺失低频信息的低频带反演效果对比见图7a和图7b。图7a为传统方法的全波形反演结果受低频信息缺失的影响较大，多处发生了明显的跳周。对观测记录做低通滤波后有效信息的频带在11～20Hz，对于线性递增初始模型，该初始模型与真实模型相差较大，在波形匹配过程中，低频信息的缺失导致多处观测记录与模拟记录相位相差半个周期以上，发生跳周，导致相应位置的梯度更新量计算错误。随着反演过程中迭代次数的增加，错误的速度更新量不断累积最终使反演结果与真实模型相差较大。图7b为基于本发明方法的低频段时间域全波形反演结果。从图中可以看出反演结果基本能将真实模型的大尺度构造反演出来且没有速度异常区域。经过滑动时窗波形时移校正后，模拟记录与观测记录的匹配程度提高，原始波形中跳周的部分被校正，使模型速度值向正确的方向更新，最终得到比较精确的反演结果，且没有速度异常区域。

以图7a所示反演结果作为初始模型，得到基于传统全波形反演方法的全频带反演结果(图8a)。由于图7a中存在跳周导致的速度异常区域，以之为初始模型继续进行全频带反演过程中，速度异常区域并不能被校正，而且该区域的错误速度更新量会随着迭代次数的增加继续累积，最终使反演结果严重偏离真实模型。

以图7b所示反演结果作为初始模型，得到基于本发明方法的全频带反演结果(图8b)。由于低频段反演结果(图7b)提供了一个比较准确的初始模型，再以此结果为基础进行全频带全波形反演，最终得到一个非常精确的反演结果，基本将真实模型准确的反演出来，且没有速度异常区域。

图9a和图9b为基于本发明方法的全频带反演结果与初始模型和真实模型的单道对比图，抽取了第60和第150两道。从图中可以看出，在初始模型非常不精确的情况下，基于本发明的全波形反演结果速度变化曲线与真实模型速度曲线大致相同，且没有速度异常值，可见本发明方法在初始速度模型不精确的情况下依然能将地下速度参数准确的反演出来。