CN107422379B - 基于局部自适应凸化方法的多尺度地震全波形反演方法 - Google Patents

Abstract

本发明涉及一种基于局部自适应凸化方法的多尺度地震全波形反演方法，其步骤：预处理，以零值序列作为初始值，初始速度模型作为起始值；截取直达波信息；得到正演模拟直达波；建立震源函数反演的目标函数；得到直达波残差、直达波残差反传波场；计算震源函数的更新梯度、更新方向并寻找步长；输出反演得到的高精度震源函数；进行衰减时窗处理；模拟时窗内的地震数据；进行局部凸化处理和分离处理；去除观测记录中的高频成分；建立最小二乘目标函数，得到波场残差；得到模型空间的残差反传波场；得到模型更新梯度；计算模型更新方向，并寻找步长；输出多尺度地震全波形反演的反演结果；输出最终反演结果。本发明在地震勘探技术领域中广泛应用。

Description

基于局部自适应凸化方法的多尺度地震全波形反演方法

技术领域

本发明涉及一种基于局部自适应凸化方法的多尺度地震全波形反演方法，属于地震勘探技术领域。

背景技术

20世纪90年代，Pratt将全波形反演推广到了频率域，提出只需要几个离散的频率就可以得到高精度的反演结果，且低频到高频的多尺度反演策略可以解决陷入局部极小值的问题。但全波形反演是一个强非线性问题，对噪声、初始模型、低频成分和震源函数非常敏感(Wang,Y.,Rao,Y.,2006；Virieux and Operto,2009)，其根本原因在于地震波传播的复杂性，即地震数据与地下介质物性参数之间的复杂变化关系 (Jannane et al.,1989；董良国，2013)。当地震数据中缺失低频成分，全波形反演结果都很容易出现周跳现象。如何缓解因地震数据缺失低频而带来的跳周问题，仍然是波形反演的研究热点。

面对波形反演的周跳问题，提出了许多解决方法。Bunks(1995)等提出基于低通滤波的多尺度全波形反演，该方法可在一定程度上缓解了全波形反演的周跳问题，但是当地震数据中缺失低频成分时，利用低通滤波滤除高频成分之后再进行全波形反演，不能很好的恢复模型的宏观构造。面向实际，也就是说，在波形反演中没有必要时刻追求全部地震信息的匹配，因为匹配部分信息就有可能解决目前地震成像中的宏观速度模型建立问题，还可以一定程度上规避匹配全部地震信息极易遇到的强烈非线性难题(董良国，2015)。

震源函数是全波形反演的一个重要参数,微小的误差都会导致观测记录与模拟记录存在不匹配问题，甚至导致整个反演向错误的方向更新,从而影响后续的岩性解释。当前求取震源函数最精确的方法是根据测井数据和井旁道数据来反演得到，但是实际工区不一定有测井数据，而且该方法工作量较大。工业上震源函数估计多是直接提取震源附近的直达波作为震源函数，该方法存在一定的激发延时误差和振幅不准确等问题。Rickett(2013)提出利用变量投影的方法可消除频率域全波形反演中震源函数的影响，得到较好的反演结果。Liu(2014)利用反褶积的方法成功消除了震源函数对频域全波形反演的影响。Ao(2015)利用交叉褶积方法消除震源函数对全波形反演的影响，并推导了相应的目标函数和梯度。Zhang(2016)利用褶积与反褶积的方法在重构低频的同时消除了震源函数的影响，有效缓解了低频缺失带来的周跳问题。在全波形反演的过程中消除震源函数的影响固然是一种好的方法，但是在消除震源函数的同时必然会去除一些细节信息，影响全波形反演的高精度成像，甚至可能会出现不真实的反演结果。

发明内容

针对上述问题，本发明的目的是提供一种基于局部自适应凸化方法的多尺度地震全波形反演方法，该方法能够解决震源函数反演不准确、地震数据中缺失低频以及低通滤波之后数据不匹配的问题。

为实现上述目的，本发明采用以下技术方案：一种基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于包括以下步骤：1)将野外采集的地震数据进行预处理，去除表面多次波，交混回响，并进行保护低频去噪处理；2)以雷克子波作为真实震源函数，并以零值序列作为震源函数反演初始值；3)根据工区地质背景和速度分析结果建立初始速度模型，将此初始速度模型作为全波形反演的起始值；4)截取地震数据中的直达波信息；5)通过震源函数反演初始值、初始速度模型和常密度声波方程，通过有限差分进行数值模拟，得到正演模拟直达波数据；6)按照全波形反演方法，即观测数据与模拟数据残差最小的原则，建立震源函数反演的目标函数；7)将观测记录的直达波和模拟记录的直达波做差，得到直达波残差；8)将直达波残差数据作为伴随震源，用反传算子作用于伴随震源上，得到模型空间的直达波残差反传波场，并存储所有时刻的直达波反传残差波场；9)计算震源函数的更新梯度；10)采用L-BFGS 优化算法计算震源函数的更新方向，并通过Wolfe收敛准则寻找相应步长；11)判断是否满足预先设定的终止条件，若满足则输出反演得到的高精度震源函数；若不满足终止条件，将反演结果作为下一个循环的初始模型，返回步骤5)；12)对观测记录进行衰减时窗处理；13)利用反演得到的震源函数、初始速度模型和常密度声波方程，利用有限差分进行数值模拟，模拟时窗内的地震数据，并对模拟数据乘以相同的衰减函数；设定自适应低通滤波的截断频率，并存储所有时刻的波场快照；14)对观测记录和模拟数据进行局部凸化处理，得到地震数据的一个子集，对地震数据进行分离处理；15) 根据迭代次数自适应的调整低通滤波的截断频率，利用自适应低通滤波去除凸化之后观测记录中的高频成分，只保留低频成分；16)按照全波形反演方法，即观测数据与模拟数据残差最小的原则，建立最小二乘目标函数：17)将观测记录和模拟记录做差，得到波场残差；18)将波场残差作为伴随震源，用反传算子作用于伴随震源上，得到模型空间的残差反传波场，并存储所有时刻的反传残差波场；19)正传波场与反传残差波场做互相关，得到模型更新梯度；20)用L-BFGS优化算法计算模型更新方向，并通过Wolfe 收敛准则寻找步长；21)判断是否满足终止条件，若满足则输出局部自适应凸化方法的多尺度地震全波形反演的反演结果；若不满足终止条件，将反演结果作为下一个循环的初始模型，返回步骤13)，更换地震数据的局部凸子集；22)得到初始速度模型后，再进行常规时间域全波形反演，输出最终反演结果。

进一步，所述步骤6)中，震源函数反演的目标函数为：

其中，表示震源附近直达波模拟记录，表示震源附近直达波观测记录。

进一步，所述震源函数反演的目标函数对震源函数f的导数表示为：

时间域声波方程两边同时对震源函数f求导得到：

其中正演算子L与震源函数无关，有即波场u对震源函数f的导数：

得到震源函数反演的更新梯度：

其中震源函数伴随震源定义为：

进一步，所述步骤9)中，震源函数的更新梯度为：

进一步，所述步骤10)中，L-BFGS算法是一种迭代型的算法，迭代更新公式如下：

m_k+1＝m_k-α_kH_kg_k

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

进一步，所述步骤11)中，所述终止条件为预先设定的迭代次数及预设精度。

进一步，所述步骤12)中，衰减时窗公式为：

其中，T_win表示衰减时窗，t表示时间序列，t₀表示衰减的断点，N表示平滑的阶数。

进一步，所述步骤14)中，对地震数据进行分离处理，如下：

D(nt,nr)＝D⁺(nt,nr)+D^-(nt,nr)

其中，D表示地震数据，D⁺表示正值地震数据，D^-表示负值地震数据，nt表示时间采样点数，Nt表示总时间采样点数，nr表示检波器的个数，Nr表示总的检波器数；将地震数据的正值与负值分开，对不同的子集分别反演。

进一步，所述步骤16)中，根据常规时间域全波形反演目标函数，将局部自适应凸化方法的全波形反演的目标函数定义为：

其中表示经过局部自适应凸化处理的模拟记录，表示经过局部自适应凸化处理的观测记录；对速度模型v的导数表示为：

得到基于局部自适应凸化方法的全波形反演的近似梯度为：

其中L^-1表示伴随算子，表示残差反传波场；称为基于局部自适应凸化方法的伴随震源，则目标函数的梯度简化为：

其中表示由伴随震源向地下传播得到的残差反传波场，P_f表示正演模拟波场。

进一步，所述步骤19)中，模型更新梯度为：

其中表示由伴随震源向地下传播得到的残差反传波场，P_f表示正演模拟波场,v为速度模型，nr表示检波器的个数，ns为震源数目，T为总接收时间。

本发明由于采取以上技术方案，其具有以下优点：1、通过对实际数据中缺失低频问题和震源函数不准确等探究，本发明提出伴随状态法震源函数反演的策略得到精确的震源函数，同时提出局部凸化的处理方法，成功恢复出地震数据中的低频信号，利用该低频信号可为全波形反演提供一个好的初始速度模型。将衰减截断时窗多尺度和自适应低通滤波等技术应用到了全波形反演中，该反演策略一定程度上提高了全波形反演的精度，且更适合实际地震数据处理，同时本发明提出的伴随状态法反演震源函数，避免了震源函数不准确导致的跳周问题。2、伴随状态震源函数反演方法，本发明利用震源附近直达波信息就可以反演得到精确的震源函数，解决了震源函数估计不准确的问题。该震源函数反演方法具有不依赖初始速度模型、不依赖初始震源函数、反演精度高等优点。可为全波形反演以及逆时偏移提供一个精确的震源函数，有效缓解了因震源函数不准确带来的周跳问题。3、本发明利用局部凸化方法能够很好的恢复出地震数据中缺失的低频信息，将得到的低频信息应用到全波形反演中，可以得到一个更接近真实模型的宏观信息，利用该宏观信息作为下一步反演的初始速度模型，能够有效缓解全波形反演的跳周问题。4、本发明利用自适应低通滤波多尺度的反演策略，有效缓解了地震数据高频成分的影响，同时能够避免因滤波导致的数据不匹配现象。该多尺度反演策略使得整套数据处理流程。5、本发明利用衰减截断时窗多尺度反演策略，极大地缩短了正演所需要的计算时间。全波形反演是一个耗时的过程，每一次做匹配的过程都需要重新正演模拟，利用衰减截断时窗的策略，每一次只需要正演衰减截断时窗内的数据，大大提高了全波形反演的计算效率。同时在数据匹配的过程中只需要匹配衰减截断时窗内的数据，大大缩减了每一步反演的参数，一定程度上减弱了全波形反演的非线性。6、本发明将指数衰减函数引入到衰减截断时窗中，缓解了数据奇异点的影响，避免在模型空间中出现的画弧现象。直接利用时窗截取数据，可能会破坏一些波形的信息，例如直接利用0.8s的时窗进行截取，地震中的某个波形可能只被一半，这种不完整的波形在反演的过程中会反映到模型空间，出现画弧的现象。7、本发明利用L-BFGS优化算法，计算速度快，反演精度高，内存消耗小。L-BFGS优化算法有效的克服了显式求取Hessian矩阵及其逆的困难，其具有超线性收敛速度，计算效率高，占用内存小，精度高等优点，较适合求解大规模非线性优化问题。

附图说明

图1是本发明的整体流程示意图；

图2a是本发明雷克子波波形图；

图2b是本发明雷克子波频谱图；

图3a是本发明单炮观测记录；

图3b是本发明截取前0.15s的单炮观测记录；

图4是本发明真实震源函数与直达波示意图；

图5a是本发明真实震源函数示意图；

图5b是本发明线性初始震源函数示意图；

图6是本发明震源函数反演的目标函数值示意图；

图7a是本发明震源函数第一次反演更新梯度示意图；

图7b是震源函数最终反演结果示意图；

图8是本发明真实震源函数与缺失低频成分的直达波示意图；

图9a是本发明震源函数反演第一次更新梯度示意图；

图9b是本发明缺失低频成分的震源函数最终反演结果示意图；

图10是本发明Butterworth低通滤波器示意图；

图11a是本发明真实速度模型示意图；

图11b是本发明线性初始速度模型示意图；

图12a是本发明地震数据局部凸子集波形图示意图；

图12b是本发明地震数据局部凸子集频谱图；

图13a是本发明原始地震记录与高通滤波地震记录波形图示意图；

图13b是本发明原始地震记录与高通滤波地震记录频谱图；

图14a是本发明地震数据局部凸子集低通滤波记录波形图；

图14b是本发明地震数据局部凸子集低通滤波记录频谱图；

图15a1是原始混合震源观测记录波形图；

图15b1是凸化处理混合震源观测记录波形图；

图15c1是低通滤波处理之后混合震源观测记录波形图；

图15a2是原始混合震源观测记录频谱图；

图15b2是凸化处理混合震源观测记录频谱图；

图15c2是自适应低通滤波处理之后混合震源观测记录频谱图；

图16a是本发明Butterworth衰减单道函数；

图16b是本发明Butterworth二维衰减窗；

图17a是本发明常规截断时窗混合震源观测记录示意图；

图17b是本发明Butterworth衰减截断时窗混合震源观测记录示意；

图18a是本发明单道常规截断时窗观测记录示意图；

图18b是本发明单道Butterworth衰减截断时窗观测记录示意图；

图19a是本发明局部自适应凸化方法伴随震源；

图19b是本发明原始记录伴随震源；

图20a是本发明混合震源记录频谱图；

图20b是本发明混合震源原始缺失低频记录频谱图；

图20c是本发明混合震源局部自适应凸化频谱图；

图21a是本发明反演得到的精确震源函数全波形反演测试结果；

图21b是本发明直达波作为震源函数全波形反演测试结果；

图22a是本发明缺失低频常规时间域全波形反演结果；

图22b是本发明缺失低频自适应低通滤波多尺度全波形反演结果；

图23a是flcut＝15Hz的局部适应凸化方法全波形反演结果；

图23b是flcut＝25Hz的局部适应凸化方法全波形反演结果；

图23c是flcut＝20Hz的常规全波形反演结果；

图23d是常规全波形反演结果。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明提出了一种基于局部自适应凸化方法的多尺度地震全波形反演方法，该方法包括以下步骤：

1)将野外采集的地震数据进行预处理,预处理包括去除表面多次波，交混回响，并进行保护低频去噪处理；

2)以雷克子波作为真实震源函数(如图2a、图2b所示)，并以零值序列作为震源函数反演初始值(如图5a、图5b所示)；

3)采用已有方法根据工区地质背景和速度分析结果建立初始速度模型，将此初始速度模型作为全波形反演的起始值(如图11a、图11b所示)；4)截取地震数据中的直达波信息(如图3a、图3b，图4所示)，确保这一部分的直达波信息尽量不含有地下模型的信息；

5)通过震源函数反演初始值、初始速度模型和常密度声波方程，通过有限差分进行数值模拟，得到正演模拟直达波数据；

6)按照全波形反演方法，即观测数据与模拟数据残差最小的原则，建立震源函数反演的目标函数：(如图6所示)，其中表示震源附近直达波模拟记录，表示震源附近直达波观测记录；

震源函数反演的目标函数对震源函数f的导数可表示为：

时间域声波方程两边同时对震源函数f求导得到：

其中正演算子L与震源函数无关，有即波场u对震源函数f的导数：

得到震源函数反演的更新梯度：

其中震源函数伴随震源定义为：本发明利用伴随状态法计算Jacobian矩阵大大缩短了反演震源函数的时间。

7)将观测记录的直达波和模拟记录的直达波做差，得到直达波残差

8)将直达波残差数据作为伴随震源，用反传算子作用于伴随震源上，得到模型空间的直达波残差反传波场P_b，并存储所有时刻的直达波反传残差波场；

9)计算震源函数的更新梯度：

10)采用L-BFGS优化算法计算震源函数的更新方向，并通过Wolfe收敛准则寻找合适的步长；

L-BFGS算法是一种迭代型的算法，迭代更新公式如下：

m_k+1＝m_k-α_kH_kg_k

其中，m_k为模型更新参数，α_k为Wolfe线性搜索得到的步长，H_k为近似Hessian 矩阵的逆矩阵，g_k为模型更新梯度，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_kg_k的乘积可以通过梯度g_k与向量对{s_k,y_k}之间一系列向量的内积与向量的和来获得的。在每一次更新完成后，上一步向量对将被当前新向量对 {s_k+1,y_k+1}取代。因此，向量对集合中包含最近n次迭代的曲率信号。在时实际中，当n≥5 时都能获得较满意的结果。L-BFGS优化算法的初始近似Hessian矩阵每一次迭代中都不同。近似Hessian矩阵的逆矩阵H_k需满足以下更新公式：

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

(1)令q＝g_k，则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矩阵及其逆的困难，其具有超线性收敛速度，计算效率高，占用内存小，精度高等优点，较适合求解大规模非线性优化问题。

11)判断是否满足预先设定的终止条件，若满足则输出反演得到的高精度震源函数 (如图7a～图9b所示)；若不满足终止条件，将反演结果作为下一个循环的初始模型，返回步骤5)；终止条件为预先设定的迭代次数及预设精度；

12)对观测记录进行衰减时窗处理，衰减时窗公式如下：

其中，T_win表示衰减时窗，t表示时间序列，t₀表示衰减的断点，N表示平滑的阶数。利用该衰减时窗避免波形被截断的现象(如图16a～图18b所示)。

13)利用反演得到的震源函数、初始速度模型和常密度声波方程，利用有限差分进行数值模拟，模拟时窗内的地震数据，并对模拟数据乘以相同的衰减函数。设定自适应低通滤波的截断频率，并存储所有时刻的波场快照；

14)对观测记录和模拟数据进行局部凸化处理，得到地震数据的一个子集，该子集数据能够减弱地震数据匹配过程中的强非线性。对地震数据进行分离处理(如图15a1～图15c2以及图20a、图20b所示)，如下：

D(nt,nr)＝D⁺(nt,nr)+D^-(nt,nr)

其中，D表示地震数据，D⁺表示正值地震数据，D^-表示负值地震数据，nt表示时间采样点数，Nt表示总时间采样点数，nr表示检波器的个数，Nr表示总的检波器数。该局部凸化处理的目的是将地震数据的正值与负值分开，对不同的子集分别反演，在此过程中地震数据中的低频成分得到了恢复，同时减弱了观测数据与模拟数据的匹配难度，降低了全波形反演的非线性程度(如图12a～图14b所示)。

15)根据迭代次数自适应的调整低通滤波的截断频率，利用自适应低通滤波去除凸化之后观测记录中的高频成分，只保留低频成分(如图14b所示)。

针对上述低通滤波之后观测记录与模拟记录不匹配问题，本发明利用自适应低通滤波多尺度反演策略，即用Butterworth低通滤波(如图10所示)对观测记录(D_obs) 进行低通滤波得到同时再对模拟记录(D_cal)进行低通滤波得到该方法利用相同的滤波器只对对观测记录和模拟记录进行低通滤波，减弱与因低通滤波产生的差异性。此时自适应低通滤波伴随震源定义为：

进而得到自适应低通滤波多尺度全波形反演的梯度：

式中，v为速度模型，ns为震源数目，T为总接收时间；

利用自适应低通滤波多尺度全波形反演方法与常规全波形反演方法的梯度计算公式基本相同，伴随状态法大大地缩短了Jacobian矩阵的计算时间。先用低主频的伴随震源恢复模型的宏观构造，逐渐增加滤波的截断频率来实现多尺度全波形反演，该多尺度策略能够有效缓解全波形反演周跳的问题。

16)按照全波形反演方法，即观测数据与模拟数据残差最小的原则，建立最小二乘目标函数：

根据常规时间域全波形反演目标函数，可将局部自适应凸化方法的全波形反演的目标函数定义为：

其中表示经过局部自适应凸化处理的模拟记录，表示经过局部自适应凸化处理的观测记录。对速度模型v的导数可表示为：

得到基于局部自适应凸化方法的全波形反演的近似梯度为：

其中L^-1表示伴随算子，表示残差反传波场。称为基于局部自适应凸化方法的伴随震源。目标函数的梯度可以简化为：

其中表示由伴随震源向地下传播得到的残差反传波场，P_f表示正演模拟波场。

17)将观测记录和模拟记录做差，得到波场残差(如图19a、图19b所示)；

18)将波场残差作为伴随震源，用反传算子作用于伴随震源上，得到模型空间的残差反传波场，并存储所有时刻的反传残差波场；

19)正传波场与反传残差波场做互相关，得到模型更新梯度：

20)用L-BFGS优化算法计算模型更新方向，并通过Wolfe收敛准则寻找合适的步长；

21)判断是否满足终止条件，若满足则输出局部自适应凸化方法的多尺度地震全波形反演的反演结果；若不满足终止条件，将反演结果作为下一个循环的初始模型，返回步骤13)，并更换地震数据的局部凸子集；

22)得到高精度的初始速度模型后，再进行常规时间域全波形反演，输出最终反演结果(如图21a～23d所示)。

综上所述，本发明利用国际通认的速度模型进行局部自适应凸化方法的多尺度地震全波形反演测试，根据说明书附图可以看出，当地震数据中缺失10Hz以下的低频信息时(如图8所示)，仍然能得到较好的反演结果，证明了本发明提出的地震数据局部凸化方法恢复出的低频信号真实有效。当地震数据缺失低频成分和初始模型与真实模型相差较远的时候，局部自适应凸化方法的多尺度地震全波形反演能够弥补常规全波形反演方法不足。再加上衰减截断时窗多尺度反演策略和自适应低通滤波多尺度反演策略一定程度上增加了反演流程的稳定性，同时很好的压制了全波形反演的跳周问题。

上述各实施例仅用于说明本发明，其中方法的实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (10)

1.一种基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于包括以下步骤：

1)将野外采集的地震数据进行预处理，去除表面多次波，交混回响，并进行保护低频去噪处理；

2)以雷克子波作为真实震源函数，并以零值序列作为震源函数反演初始值；

3)根据工区地质背景和速度分析结果建立初始速度模型，将此初始速度模型作为全波形反演的起始值；

4)截取地震数据中的直达波信息；

5)通过震源函数反演初始值、初始速度模型和常密度声波方程，通过有限差分进行数值模拟，得到正演模拟直达波数据；

6)按照全波形反演方法，即观测数据与模拟数据残差最小的原则，建立震源函数反演的目标函数；

7)将观测数据的直达波和模拟数据的直达波做差，得到直达波残差；

8)将直达波残差数据作为伴随震源，用反传算子作用于伴随震源上，得到初始速度模型空间的直达波残差反传波场，并存储所有时刻的直达波残差反传波场；

9)计算震源函数的更新梯度；

10)采用L-BFGS优化算法计算震源函数的更新方向，并通过Wolfe收敛准则寻找相应步长；

11)判断是否满足预先设定的终止条件，若满足则输出反演得到的高精度震源函数；若不满足终止条件，将反演结果作为下一个循环的初始速度模型，返回步骤5)；

12)对观测数据进行衰减时窗处理；

13)利用反演得到的震源函数、初始速度模型和常密度声波方程，利用有限差分进行数值模拟，模拟时窗内的地震数据，并对模拟数据乘以相同的衰减函数；设定自适应低通滤波的截断频率，并存储所有时刻的波场快照；

14)对观测数据和模拟数据进行局部凸化处理，得到地震数据的一个子集，对地震数据进行分离处理；

15)根据迭代次数自适应的调整低通滤波的截断频率，利用自适应低通滤波去除凸化之后观测数据中的高频成分，只保留低频成分；

16)按照全波形反演方法，即观测数据与模拟数据残差最小的原则，建立最小二乘目标函数：

17)将观测数据和模拟数据做差，得到波场残差；

18)将波场残差作为伴随震源，用反传算子作用于伴随震源上，得到初始速度模型空间的残差反传波场，并存储所有时刻的残差反传波场；

19)正传波场与残差反传波场做互相关，得到初始速度模型更新梯度；

20)用L-BFGS优化算法计算初始速度模型更新方向，并通过Wolfe收敛准则寻找步长；

21)判断是否满足终止条件，若满足则输出局部自适应凸化方法的多尺度地震全波形反演的反演结果；若不满足终止条件，将反演结果作为下一个循环的初始速度模型，返回步骤13)，更换地震数据的局部凸子集；

22)得到初始速度模型后，再进行常规时间域全波形反演，输出最终反演结果。

2.如权利要求1所述的基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于：所述步骤6)中，震源函数反演的目标函数为：

其中，表示震源附近直达波模拟数据，表示震源附近直达波观测数据。

3.如权利要求2所述的基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于：所述震源函数反演的目标函数对震源函数f的导数表示为：

时间域声波方程两边同时对震源函数f求导得到：

其中正演算子L与震源函数无关，有即波场u对震源函数f的导数：

得到震源函数反演的更新梯度：

其中震源函数伴随震源定义为：

4.如权利要求1所述的基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于：所述步骤9)中，震源函数的更新梯度为：

式中，f为震源函数，L为正演算子，表示震源附近直达波模拟数据，表示震源附近直达波观测数据。

5.如权利要求1所述的基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于：所述步骤10)中，L-BFGS优化算法是一种迭代型的算法，迭代更新公式如下：

m_k+1＝m_k-α_kH_kg_k

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

6.如权利要求1所述的基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于：所述步骤11)中，所述终止条件为预先设定的迭代次数及预设精度。

7.如权利要求1所述的基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于：所述步骤12)中，衰减时窗公式为：

其中，T_win表示衰减时窗，t表示时间序列，t₀表示衰减的断点，N表示平滑的阶数。

8.如权利要求1所述的基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于：所述步骤14)中，对地震数据进行分离处理，如下：

D(nt,nr)＝D⁺(nt,nr)+D^-(nt,nr)

其中，D表示地震数据，D⁺表示正值地震数据，D^-表示负值地震数据，nt表示时间采样点数，Nt表示总时间采样点数，nr表示检波器的个数，Nr表示总的检波器数；将地震数据的正值与负值分开，对不同的子集分别反演。

9.如权利要求1所述的基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于：所述步骤16)中，根据常规时间域全波形反演目标函数，将局部自适应凸化方法的全波形反演的目标函数定义为：

其中表示经过局部自适应凸化处理的模拟数据，表示经过局部自适应凸化处理的观测数据，nr表示检波器的个数，ns为震源数目，T为总接收时间；对初始速度模型v的导数表示为：

得到基于局部自适应凸化方法的全波形反演的近似梯度为：

其中L^-1表示伴随算子，表示残差反传波场；称为基于局部自适应凸化方法的伴随震源，则目标函数的梯度简化为：

其中表示由伴随震源f_s ^conv向地下传播得到的残差反传波场，P_f表示正演模拟波场，u为波场。

10.如权利要求1所述的基于局部自适应凸化方法的多尺度地震全波形反演方法，其特征在于：所述步骤19)中，初始速度模型更新梯度为：

其中表示由伴随震源f_s ^conv向地下传播得到的残差反传波场，P_f表示正演模拟波场,v为初始速度模型，nr表示检波器的个数，ns为震源数目，T为总接收时间。