基于振幅衰减与线性插值的时间域全波形反演方法

摘要

本发明涉及一种基于振幅衰减与线性插值的时间域全波形反演方法，将产生跳周的原因分为两部分：第一部分是观测记录与模拟记录波形极性相反的波形；第二部分是观测记录与模拟记录极性相同但是相位差大于半个周期的波形。对第一部分波形采取振幅衰减的方法，即对这部分模拟记录的振幅乘以一个衰减函数使之减小，从而使跳周部分的波形对梯度的计算干扰减小；对第二部分波形采取线性插值的方法，使观测记录逐渐逼近模拟记录，提高博波形的相关性，从而减少跳周的发生。本发明采用全局互相关目标函数减小反演对振幅信息的依赖。基于振幅衰减与线性插值的时间域全波形反演方法在不降低计算效率的同时减少了跳周的发生，极大地提高了全波形反演的精度。

权利要求

1.一种基于振幅衰减与线性插值的时间域全波形反演方法，其特征在于，包括以下步骤：a、对实际地震观测记录进行子波估计、低频保护去噪、缺失地震道补偿多次波衰减、面波切除以及消除交混回响预处理；

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

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

d、实际地震勘探采集到的观测记录记为d，对初始模型进行正演模拟得到模拟记录记为u，首先对观测记录和模拟记录的振幅逐道进行编码：x为地表检波器个数，i为地震道数，ui为模拟记录中的第i道，t为采样点， 为模拟记录中第t个采样点的振幅极性， 为模拟记录中第t个采样点的振幅绝对值，同样，能够得到观测记录的振幅编码：di为观测记录中的第i道， 为观测记录中第t个采样点的振幅极性， 为观测记录中第t个采样点的振幅绝对值，取：Izer为观测记录与模拟记录振幅极性相反处的对应采样点t的集合；然后对模拟数据以观测数据为目标进行线性插值迭代：u*＝u+(d-u)/ε,(ε＞1)

u*为进行插值迭代后的模拟记录，ε为大于1的常数，设：th为迭代停止的阈值，即经过插值计算后的模拟数据与观测数据的归一化零延迟互相关，所以有th∈[-1,1]，对u*再进行振幅编码得到：为插值迭代后模拟记录的第i道， 为 中第t个采样点的振幅极性， 为 中第t个采样点的振幅绝对值，对 中Izer位置的采样点做振幅衰减：α为一个正整数，将做完插值迭代和振幅衰减后的模拟记录记为e、根据全局互相关原理建立目标函数：

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

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

Pf为时间域正传波场,Pb为反传波场；

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

说明书

基于振幅衰减与线性插值的时间域全波形反演方法

技术领域

[0001] 本发明属于地震勘探技术领域，涉及一种地震勘探领域的反演方法，具体涉及一种时间域全波形反演，特别涉及一种基于振幅衰减与线性插值的时间域全波形反演方法。

背景技术

[0002] 油气能源在当今世界中的作用是不可替代的，油气资源的勘探与开发甚至决定了许多国家的政治导向。因此，油气的勘探与开发迎来了前所未有的机遇和空前的发展。随着勘探与开发的不断深入，陆地和海洋浅层油气田越来越少，深地和深水油气变得愈发重要。地震勘探是深部油气资源勘探的一种方法，其以勘探深度大精度高成为实际生产工作中最主要的方法。它是通过对采集到的地震记录做一定的处理如速度分析、偏移以及反演等，可以得到丰富的地下介质信息如速度、密度、反射系数以及构造走向等。在诸多反演方法中，全波形反演近年来成为地球物理学家研究的热点课题，它是通过对待重建的初始模型正演所得到的模拟记录与实际勘探得到的观测记录进行匹配，得到模型参数更新量再通过优化算法迭代更新最终得到地下介质参数的真实分布。全波形反演方法与传统地震反演成像方法的最大区别在于：其利用了叠前地震波场的全波信息包括直达波、折射波、反射波、多次波甚至噪声等直接进行反演计算，因此包含了丰富的地下参数信息，是目前反演精度最高的方法。

[0003] 全波形反演方法最初在20世纪80年代由Tarantola在时间域实现的。其将地震勘探领域的地下介质参数反演问题，看成模拟地震记录与实际采集地震记录残差在最小二乘目标函数下的最优化问题。由于该问题是高度非线性以及病态的。因此，在Born近似的理论框架下，该优化问题可以采用局部优化的方法近似求解。反演的梯度用伴随状态法求取，即将目标函数对模型参数求导通过以记录残差为震源的反传波场与正传波场的零延迟互相关得到。伴随状态法的提出避免了直接计算Jacobi矩阵，即一次梯度的求取只需要计算两次波场模拟，从而大大减少计算量，奠定了全波形反演的基础。但是受计算机性能的限制以及反演过程需要很多次迭代，计算效率的问题很长时间内制约了该方法的发展。到了20世纪90年代，Pratt等将全波形反演方法在频率域实现，即采用从低到高的几个离散的频率就可以实现全波形反演，同时从低频到高频的多尺度策略避免了反演陷入局部极小值。这大大推动了之后学者对全波形反演的研究，使该方法发展的越来越快。

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

发明内容

[0005] 本发明的目的就在于针对上述现有技术的不足，提供一种解决跳周问题的基于振幅衰减与线性插值的时间域全波形反演方法。

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

[0007] 本发明基于振幅衰减与线性插值的时间域全波形反演方法，核心是将产生跳周的原因分为两部分，第一部分是观测记录与模拟记录波形极性相反的波形，第二部分是观测记录与模拟记录极性相同但是相位差大于半个周期的波形。对第一部分波形采取振幅衰减的方法，即对这部分模拟记录的振幅乘以一个衰减函数使之减小，从而使跳周部分的波形对梯度的计算干扰减小；对第二部分波形采取线性插值的方法，使观测记录逐渐逼近模拟记录，提高了波形的相关性，从而减少跳周的发生。由于整个过程对波形的振幅做了较大的调整，而传统基于最小二乘目标函数的全波形反演方法对振幅信息非常敏感，很可能导致反演失败，因此本发明采用全局互相关目标函数减小反演对振幅信息的依赖。

[0008] 一种基于振幅衰减与线性插值的时间域全波形反演方法，包括以下步骤：

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

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

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

[0012] d、实际地震勘探采集到的观测记录记为d，对初始模型进行正演模拟得到模拟记录记为u。观测记录与模拟记录波形跳周的部分可分为两种情况，一种是二者波形振幅极性相反，另一种是波形极性相同。为了避免这两种跳周的情况，本发明分别采用两种方法对发生跳周的波形进行校正。首先对观测记录和模拟记录的振幅逐道进行编码：

[0013]

[0014]

[0015] x为地表检波器个数，i为地震道数，ui为模拟记录中的第i道，t为采样点， 为模拟记录中第t个采样点的振幅极性， 为模拟记录中第t个采样点的振幅绝对值，。同样，可以得到观测记录的振幅编码：

[0016]

[0017]

[0018] di为观测记录中的第i道， 为观测记录中第t个采样点的振幅极性， 为观测记录中第t个采样点的振幅绝对值，取：

[0019]

[0020] Izer为观测记录与模拟记录振幅极性相反处的对应采样点t的集合。然后对模拟数据以观测数据为目标进行线性插值迭代：

[0021] u*＝u+(d-u)/ε,(ε＞1)

[0022] u*为进行插值迭代后的模拟记录，ε为大于1的常数，设：

[0023]

[0024] th为迭代停止的阈值，即经过插值计算后的模拟数据与观测数据的归一化零延迟互相关，所以有th∈[-1,1]。通过设计阈值，使模拟记录与观测记录逐渐逼近，可以一定程度上消除波形相同部分的跳周问题。对u*再进行振幅编码得到：

[0025]

[0026]

[0027] 为插值迭代后模拟记录的第i道， 为 中第t个采样点的振幅极性， 为中第t个采样点的振幅绝对值。为了消除观测记录与模拟记录振幅极性相反造成的跳周问题，对 中Izer位置的采样点做振幅衰减：

[0028]

[0029] α为一个正整数，将做完插值迭代和振幅衰减后的模拟记录记为

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

[0031]

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

[0033]

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

[0035]

[0036] Pf为时间域正传波场,Pb为反传波场；

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

[0038] 与现有技术相比，本发明的有益效果在于：本发明提出了一种基于振幅衰减与线性插值的时间域全波形反演方法，将导致跳周问题产生的波形分为两部分分别进行处理：观测记录与模拟记录振幅极性相反的部分做振幅衰减；观测记录与模拟记录振幅极性相同的地方做线性插值，以此消除绝大部分波形跳周的情况，从而减少反演中速度更新量计算错误的情况，使反演结果精度提高。

[0039] 解决了以下问题：

[0040] 1、有效的减少了跳周的发生。本发明将产生跳周的原因分为两部分，第一部分是观测记录与模拟记录波形极性相反的波形，第二部分是观测记录与模拟记录极性相同但是相位差大于半个周期的波形。对第一部分波形采取振幅衰减的方法，即对这部分模拟记录的振幅乘以一个衰减函数使之减小，从而使跳周部分的波形对梯度的计算干扰减小；对第二部分波形采取线性插值的方法，使观测记录逐渐逼近模拟记录，提高了波形的相关性，从而减少跳周的发生。

[0041] 2、采用了全局互相关目标函数减小反演对振幅信息的依赖。由于对两部分模拟记录的波形分别进行了振幅衰减和线性插值，处理后的模拟记录振幅与真实幅值相差较大，而传统基于最小二乘目标函数的全波形反演方法对振幅信息非常敏感，很可能导致反演失败，因此本发明采用全局互相关目标函数减小反演对振幅信息的依赖。同时互相关计算使远偏移距记录的权重增加，而远偏移距记录往往携带了大尺度构造的信息，因此能进一步加强对模型大尺度构造的反演精度，为后续全频带反演提供更加精确的初始模型。

[0042] 3、只需要对模拟记录做振幅衰减和线性插值两步简单操作就可以完成对模拟记录的处理，因此相比于传统时间域全波形反演几乎没有计算效率上的下降，在不损失计算效率的同时提高反演精度。

[0043] 基于振幅衰减与线性插值的时间域全波形反演方法原理简单易于编程实现，在保持计算效率不变的情况下有效的减少了跳周的发生，非常适用于实际地震勘探，该方法能使地震勘探反演效果显著提高。

附图说明

[0044] 图1基于振幅衰减与线性插值的时间域全波形反演方法流程图；

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

[0046] 图2b雷克子波频谱；

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

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

[0049] 图4真实模型图；

[0050] 图5线性递增初始模型图；

[0051] 图6a原始模拟记录与观测记录各地震道归一化互相关值(α＝64，th＝0.85)；

[0052] 图6b基于本发明方法处理后的模拟记录与观测记录各地震道归一化互相关值(α＝64，th＝0.85)；

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

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

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

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

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

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

具体实施方式

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

[0060] 如图1所示，本发明基于振幅衰减与线性插值的时间域全波形反演方法，将导致跳周问题产生的波形分为两部分分别进行处理：观测记录与模拟记录振幅极性相反的部分做振幅衰减；观测记录与模拟记录振幅极性相同的地方做线性插值，以此消除绝大部分波形跳周的情况，从而防止反演中速度更新量计算错误，使反演结果精度提高，包括以下步骤：

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

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

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

[0064] d、实际地震勘探采集到的观测记录设为d，对初始模型进行正演模拟得到模拟记录u。观测记录与模拟记录波形跳周的部分可分为两种情况，一种是二者波形振幅极性相反，另一种是波形极性相同。为了避免这两种跳周的情况，本发明分别采用两种方法对发生跳周的波形进行校正。首先对观测记录和模拟记录的振幅逐道进行编码：

[0065]

[0066]

[0067] x为地表检波器个数，i为地震道数，ui为模拟记录中的第i道，t为采样点， 为模拟记录中第t个采样点的振幅极性， 为模拟记录中第t个采样点的振幅绝对值，。同样，可以得到观测记录的振幅编码：

[0068]

[0069]

[0070] di为观测记录中的第i道， 为观测记录中第t个采样点的振幅极性， 为观测记录中第t个采样点的振幅绝对值。取：

[0071]

[0072] Izer为观测记录与模拟记录振幅极性相反处的对应采样点t的集合。然后对模拟数据以观测数据为目标进行线性插值迭代：

[0073] u*＝u+(d-u)/ε,(ε＞1)

[0074] u*为进行插值迭代后的模拟记录，ε为大于1的常数。设：

[0075]

[0076] th为迭代停止的阈值，即经过插值计算后的模拟数据与观测数据的归一化零延迟互相关，所以有th∈[-1,1]。通过设计阈值，使模拟记录与观测记录逐渐逼近，可以一定程度上消除波形相同部分的跳周问题。

[0077] 对u*再进行振幅编码得到：

[0078]

[0079]

[0080] 为插值迭代后模拟记录的第i道， 为 中第t个采样点的振幅极性， 为中第t个采样点的振幅绝对值。为了消除观测记录与模拟记录振幅极性相反造成的跳周问题，对 中Izer位置的采样点做振幅衰减：

[0081]

[0082] α为一个正整数。将做完插值迭代和振幅衰减后的模拟记录记为

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

[0084]

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

[0086]

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

[0088]

[0089] 简化后可表示为：

[0090]

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

[0092]

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

[0094]

[0095] 其中Pf为时间域正传波场，L-1λ表示反传波场，其中L-1为反传算子。

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

[0097] mk+1＝mk-αkHkgk

[0098] 其中，Hk为近似Hessian矩阵的逆矩阵，mk为模型更新参数，gk为模型更新梯度，αk为Wolfe线性搜索得到的步长，k表示迭代次数。

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

[0100] Hk+1＝VkTHkVk+ρkskskT

[0101] Vk＝I-ρkykskT

[0102] sk＝mk+1-mk，yk＝gk+1-gk

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

[0104]

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

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

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

[0108] (3)通过上述步骤得到更新量Hkgk＝r。

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

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

[0111] 本发明搭建了MATLAB并行工作库的安装环境，并安装MATLAB并行计算工具箱Parallel Computing Toolbox。

[0112] 实施例1

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

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

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

[0116]网格大小 网格距 测线长度 纵向深度 fluc 频带宽度
69×192 12.5m 2400m 862.5m 20Hz 11～20Hz

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

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

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

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

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