CN107765302B - 不依赖震源子波的时间域单频波形走时反演方法 - Google Patents

Abstract

本发明涉及一种不依赖震源子波的时间域单频波形走时反演方法，针对地震数据缺失低频成分导致全波形反演周波跳跃现象：本发明提出时间域单频波形来进行波动方程走时反演，来为常规全波形反演构建一个高精度的初始速度模型，进而解决周波跳跃问题，最终得到高精度的全波形反演反演结果。本发明是为了降低全波形反演对低频成分的依赖性，同时避免震源子波的不准确导致的波形不匹配问题，并为全波形反演在实际应用提供更多的技术支持。该方法极大程度的缓解了全波形反演过程中的周波跳跃现象，改善了反演精度，并为常规全波形反演提供一个高精度初始速度模型。在解决全波形反演周波跳跃的同时，还能够避免震源子波不准确带来的影响。

Description

不依赖震源子波的时间域单频波形走时反演方法

技术领域

本发明涉及一种地震勘探的数据处理方法，尤其是利用单频波形的走时差信息来进行反演，并利用相位校正技术来消除震源子波不准确对反演过程中的影响。整个计算过程采用分炮并行策略，目的是加快整个波动方程走时反演的进程。

背景技术：

地震勘探全波形反演，是根据采集到的地震波波形的信息来得到地下的介质信息。传统的如偏移速度分析技术和走时层析成像技术都只能得到宏观的速度场，无法得到高精度的速度模型。随着石油工业的发展和勘探开发程度的不断深入，油气勘探开发从构造勘探阶段逐渐走向岩性勘探阶段，油气勘探急需一种反演方法能够得到复杂构造的高精度速度模型，为此，全波形反演迅速发展起来，并成为当今油气勘探开发的研究热点。全波形反演是一个基于地震全波场模拟的数据拟合过程，几乎使用了地震记录中所有的有效信息，该方法利用优化算法不断进行搜索，最终找到模拟数据与实际数据拟合差最小的速度模型。但是全波形反演是一个强非线性优化过程，同时野外采集的地震数据中常常缺失低频成分，导致了在反演过程中很容易出现周波跳跃问题。

为了克服全波形反演的周波跳跃现象，提出了许多解决方法。地震数据属于复频信号，因此可以通过简单的傅里叶变换以及反傅里叶变换提取到时间域地震数据的单频波形。将地震数据拆分开来进行反演能够有效避免因地震波形过于复杂而导致的周波跳跃现象。例如Bunks(1995)将时间域的波形分解到低中高三个尺度进行全波形反演，有效的解决了因复杂波形导致的周波跳跃现象，类似的还有Pratt(1999)提出的频域多尺度策略。这些分步分尺度的反演方法证明了将地震数据按照不同频率尺度拆解开来，分别进行反演是一种行之有效的反演方法。

震源子波是全波形反演过程中的一个重要参数，微小的误差都会导致观测记录与模拟记录不匹配问题，甚至导致整个反演向错误的方向更新。但实际地震数据中，震源函数往往是未知的，这成为全波形反演的又一巨大挑战。例如Ricket(2013)提出利用变量投影的方法可消除频率域全波形反演中震源子波对反演结果的影响，得到较好的反演结果；Lee(2003)利用反褶积的方法成功消除了震源函数对频域全波形反演的影响；Choi(2011)利用观测数据和模拟数据中的参考地震道信息进行交叉褶积，最终在时间域消除震源子波对全波形反演的影响，并推导了相应的目标函数和梯度，给出了详细的数学推导；Zhang等(2016)利用褶积与反褶积的方法在重构低频的同时消除了震源函数的影响，有效缓解了低频缺失和震源函数不准确带来的周波跳跃现象。

发明内容：

本发明的目的就在于针对上述现有技术的不足，提供一种不依赖震源子波的时间域单频波形走时反演方法。

本发明的发明思想是：针对地震数据缺失低频成分导致全波形反演周波跳跃现象：本发明提出时间域单频波形来进行波动方程走时反演，来为常规全波形反演构建一个高精度的初始速度模型，进而解决周波跳跃问题，最终得到高精度的全波形反演反演结果。利用时间域单频波形来进行波动方程走时反演方法不同于传统的波动方程走时反演方法，例如：单频波形中含有地震波的全波信息，不像常规波动方程走时反演只利用了初至波的信息，因而单频波形波动方程走时反演方法即使在地震数据缺失低频成分的情况下，依然能够得到高精度的反演结果。

针对震数据震源子波提取不准，导致反演错误问题：提取震源附近地震观测数据和模拟数据的直达波波形，并求解在单频情况两者的相位差，利用两者的相位差来对整个单炮记录进行校正，这样通过相对相位校正的方法能够避免因震源子波不准确产生的波形不匹配现象。观测数据和模拟数据的相位差异主要来自于震源子波的和速度模型的作用，当我们把震源子波的相位偏差因素消除了，剩下的波形差异就只是由于模型的因素导致的，这正是为反演提供了理论基础，与前期准备。

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

1).提取单频波形：利用傅里叶变换得到频域观测数据和模拟数据，提取目标频率的数据，并把其余数据赋值为0，再利用反傅里叶变换得到单频的时间域观测数据和模拟数据。

2).构建单频波形走时反演目标函数及梯度：地震数据反演实际上观测数据与模拟数据的匹配过程，本发明首先根据观测数据与模拟数据单频波形的走时差来构建一个L2范数目标函数，这样就可以用目标函数来衡量两个数据的走时接近程度，以及判断模型更新过程是否使得目标函数下降，确保速度模型向正确的方向更新。继而构建走时差与地震波场的连接方程，并利用用波动方程正演模拟以及伴随状态法来求得速度模型的更新梯度。

3).利用直达波形的相位校正技术，消除震源子波对地震反演的影响：观测数据和模拟数据的相位差异主要来自于震源子波的和速度模型的作用，因此必须把震源子波的相位偏差因素消除，才能确保目标函数的差异都是由于地下速度模型的差异所产生。通过提取震源附近的直达波地震数据求解在单频情况两者的相位差，并利用这相位差来对整个单炮记录进行校正。

4).利用循环频率多尺度策略进一步提高单频波形波动方程走时反演精度：单频波形波动方程走时反演具有频率域的优势，同时具有时间域的优势。例如：单频波形波动方程走时反演是利用离散的频率数据进行反演的方法，因此它具有频率域全波形反演从低频到高频的多尺度的策略，同时可以用几个离散的频率就可以得到高精度的反演结果。单频波形波动方程走时反演是利用时间域有限差分方法进行的数值模拟，有效避免了在三维情况下频域求解大型线性方程组Ax＝b，这样为波动方程类反演方法节约了大量的计算时间和机器内存。时间域单频波形波动方程走时反演方法结合了时间域全波形反演和频率域全波形反演的优点，同时很好地避开了这些方法中的缺点。

5).利用伪伴随震源进行反传得到反传波场：地震勘探的反演是根据目标函数实施的，不同的目标函数对应着不同的伴随震源求法。例如：波动方程走时反演方法需要对走时差进行求偏导数，需要连接方程把走时与波场连接起来，这样就可以通过反传波场与正传波场互相关来求得走时反演的更新梯度。因此首先需要求得该反演方法的伪伴随震源，利用新的伴随震源进行反传来求得更新梯度。通过伴随状态法求得的梯度很大程度上节约了反演的计算时间。

不依赖震源子波的时间域单频波形走时反演方法，包括以下步骤：

a、安装MATLAB2013a正版软件，并安装基于MATLAB语言编写的地震数据处理软件Crews工具包；

b、对实际采集的地震数据进行预处理，并将观测记录和观测系统输入计算机准备进行不依赖震源子波的时间域单频波形走时反演；

c、对观测数据进行傅里叶变换，在频率域挑选21个频率数据进行反演，从15Hz到25Hz，如图12所示。先利用低频的数据，再逐渐增加频率以实现多尺度的反演策略；

d、建立线性递增初始模型，该速度模型作为反演的起始值，通过计算模型的更新方向，并加到初始速度模型上，来实现由线性初始模型逐渐向真实模型逼近的过程；

e、利用输入的观测系统、震源子波和线性初始速度模型进行波动方程正演模拟，得到模拟数据。并提取模拟数据震源附近直达波波形，用于震源波形相位校正，消除震源子波对全波形反演的影响；

f、将提取的观测数据直达波波形和模拟数据直达波波形进行傅里叶变换，并挑选对应的21个频率数据，求取这两个数据对应频率的相位差。并将求得的相位差用于单频观测数据的相位校正中；

g、对相位校正以后的单频观测数据进行傅里叶逆变换，得到时间域的单频观测数据；

h、对模拟数据和震源子波进行傅里叶变换，并提取对应21个频率数据，再对单频模拟数据和震源子波进行傅里叶逆变换，得到单频的时间域模拟数据和单频时间域震源子波；

i、利用单频时间域震源子波和速度模型进行波动方程有限差分正演模拟，得到单频的正传波场。并对正传波场进行存储，用于后续的梯度计算；

j、利用互相关方法计算单频观测数据和单频模拟数据的走时差；

k、根据最小二乘原理构造走时差目标函数：其中Δτ表示单频观测数据与单频模拟数据存在的走时差，v表示速度值，s表示震源数目，r表示检波器数目，x_s表示震源位置，x_r表示检波器位置。在求目标函数的梯度过程中，对目标函数两端关于速度v求导，得到梯度表达式：

则伪伴随震源的表达式可以表示为：

其中u(f₀,t,x_s)_cal表示单频正传波场，u(f₀,x_r,t,x_s)_cal单频模拟数据，f₀表示反演频率，表示单频模拟数据对时间的一阶导数，d(f₀,x_r,t+Δτ,x_s)_obs表示单频观测数据，表示单频观测数据对时间的一阶导数，(L^-1)^T表示对反传算子的转置；

l、计算伪伴随震源，并利用反传算子将伪伴随震源反传到模型空间，得到反传波场；

m、正传波场与反传残差波场做互相关得到模型更新梯度；

n、利用L-BFGS优化算法计算模型更新方向，并通过Wolfe收敛准则寻找最合适的步长，按照公式v_k+1＝v_k-α_kH_kg_k对速度模型进行更新，其中k表示迭代的次数，H_k表示当前步骤伪海森矩阵的逆，α_k表示当前步骤的更新步长，g_k表示当前步骤的更新梯度，v_k表示当前步骤的速度模型，v_k+1表示经过本次迭代以后的速度模型；

o、判断是够满足终止条件，若满足输出不依赖震源子波的时间域单频波形走时反演结果。若不满足终止条件，利用循环频率多尺度策略逐渐提高频率，将反演结果作为下一个循环的初始速度模型，返回e步骤；

p、将上一步的输出结果作为常规时间域全波形反演的初始速度模型，然后进行常规时间域全波形反演，并输出最终反演结果。

有益效果：提取观测数据与模拟数据在震源附近的直达波波形，求解在单频情况两者的相位差，并利用求得的相位差来对整个单炮单频观测记录进行相位校正，避免因震源子波不同而导致全波形反演过程中的相位不匹配问题。同时单频波形变化简单，差异性较小，并结合波动方程走时反演的思想能够有效缓解全波形反演因缺失低频带来的周波跳跃现象。本发明最终的目的是利用单频走时反演方法结合相位校正策略，对不含有低频成分和不知道震源子波的地震数据进行反演，得到一个高精度的初始速度模型；然后再利用不依赖震源子波的单频波形走时反演的结果作为常规全波形反演的初始速度模型，最终可以得到一个高精度的速度模型反演结果。本发明最基本的思想就是降低全波形反演对低频成分的依赖性，同时避免震源子波的不准确导致的波形不匹配问题，为全波形反演在实际应用提供更多的技术支持。

极大程度上缓解了全波形反演过程中的周波跳跃现象，改善了反演精度，并为常规全波形反演提供一个高精度初始速度模型。

在解决全波形反演周波跳跃的同时，还能够避免震源子波不准确带来的影响。单频波形的观测数据与模拟数据只存在走时信息的差别，走时信息相比于波形信息具有更强的线性关系，而且利用互相关函数求取走时信息差不存在周波跳跃现象，所以不依赖震源子波的时间域单频波形走时反演方法在未来的实际地震数据处理中有着很好的应用前景。

附图说明

图1不依赖震源子波的时间域单频波形走时反演方法流程图。

图2观测数据与模拟数据

(左)波形图，(右)频谱图。

图3观测数据与模拟数据单频波形图。

图4目标函数值

(左)复频波形目标函数图，(右)单频波形目标函数图。

图5单频单炮记录

(左)单频单炮观测数据，

(中)单频单炮模拟数据，

(右)单频单炮模拟数据与观测数据的差值，

图6不同震源子波单炮记录

(左)初始震源子波正演模拟得到的单炮记录，

(中)真实震源子波正演模拟得到的单炮记录，

(右)初始震源子波与真实震源子波正演模拟记录的差值。

图7地震子波与直达波波形

(左)真实震源子波与初始震源子波，

(右)真实震源子波与初始震源子波各自正演模拟得到的直达波波形图。

图8时间域单频波形图

(左)单频观测数据与模拟数据，

(右)相位校正以后单频观测数据与模拟数据。

图9单频单炮差值图

(左)单频单炮观测数据与模拟数据差剖面，

(右)相位校正以后单频单跑观测数据与模拟数据差剖面。

图10单频震源子波时间域正演模拟波场快照

(左)5Hz单频震源子波正演模拟波场快照，

(中)15Hz单频震源子波正演模拟波场快照，

(右)22Hz单频震源子波正演模拟波场快照。

图11循环频率多尺度频率分布图。

图12速度模型

(左)真实速度模型，

(右)初始速度模型。

图13震源子波

(左)波形图，

(右)频谱图。

图14常规全波形反演结果图(缺失低频测试)。

(左)常规全波形反演结果，

(右)常规多尺度全波形反演结果。

图15缺失低频对反演结果影响测试图。

(左)单频波形走时多尺度反演结果，

(右)单频波形走时多尺度反演+常规多尺度全波形反演结果。

图16震源子波对反演结果影响测试图。

(左)单频波形走时多尺度反演结果，

(右)不依赖震源子波的单频波形走时多尺度反演图。

具体实施方式

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

不依赖震源子波的时间域单频波形走时反演方法，包括以下步骤：

a、安装MATLAB2013a正版软件，并安装基于MATLAB语言编写的地震数据处理软件Crews工具包；

b、对实际采集的地震数据进行预处理，并将观测记录和观测系统输入计算机准备进行不依赖震源子波的时间域单频波形走时反演；

c、对观测数据进行傅里叶变换，在频率域挑选21个频率数据进行反演，从15Hz到25Hz，如图13所示。先利用低频的数据，再逐渐增加频率以实现多尺度的反演策略。

其中傅里叶变换可以表示为：

傅里叶逆变换可以表示为：

其中ω₀表示挑选的对应角频率，d(t)表示时间域观测数据，表示频率域观测数据，F[·]表示正傅里叶变换，F^-1[·]表示傅里叶逆变换；

d、建立线性递增初始模型，该速度模型作为反演的起始值，通过计算模型的更新方向，并加到初始速度模型上，来实现由线性初始模型逐渐向真实模型逼近的过程；

e、利用输入的观测系统、震源子波和线性初始速度模型进行波动方程正演模拟，得到模拟数据。并提取模拟数据震源附近直达波波形，用于震源波形相位校正，消除震源子波对全波形反演的影响。其中提取直达波波形的时窗可以表示为：

为了避免存在突变现象，给时窗加上圆滑处理，时窗可以表示为：

观测数据经过截取以后的数据可以表示为：

其中t表示时间，t_c表示截断的时间，n表示圆滑的阶数，W表示时窗函数，G_d表示观测数据的格林函数，x_s表示震源位置，r_d表示真实震源子波；

f、将提取的观测数据直达波波形和模拟数据直达波波形进行傅里叶变换，并挑选对应的21个频率数据，求取这两个数据对应频率的相位差。并将求得的相位差用于单频观测数据的相位校正中，测试效果如图7,8,9。首先需要对观测数据和模拟数据进行归一化，在频率域归一化可以表示为：

其中φ(ω₀,x_s)_u表示模拟数据的相位，φ(ω₀,x_s)_d表示观测数据的相位，i表示虚部单位。经过归一化的观测数据与模拟数据可以写成如下形式：

假设海洋勘探中水层的速度已知，观测数据的格林函数和模拟数据的格林函数相等：

G_u(ω₀,x_s,W)＝G_d(ω₀,x_s,W)

则观测数据与模拟数据的相位差可以表示为：

其中Im[·]表示取数据的虚部。经过相位校正以后的观测数据可以表示为：

g、对相位校正以后的单频观测数据进行傅里叶逆变换，得到时间域的单频观测数据。

h、对模拟数据和震源子波进行傅里叶变换，并提取对应21个频率数据，再对单频模拟数据和震源子波进行傅里叶逆变换，得到单频的时间域模拟数据和单频时间域震源子波。

i、利用单频时间域震源子波和速度模型进行波动方程有限差分正演模拟，得到单频的正传波场。并对正传波场进行存储，用于后续的梯度计算；

j、利用互相关方法计算单频观测数据和单频模拟数据的走时差；

k、根据最小二乘原理构造走时差目标函数：其中Δτ表示单频观测数据与单频模拟数据存在的走时差，v表示速度值，s表示震源数目，r表示检波器数目，x_s表示震源位置，x_r表示检波器位置。建立单频波形走时信息与波场之间的连接方程：

通过互相关连解方程来求取单频波形的走时差：

C(f₀,x_r,Δτ,x_s)＝max{C(f₀,x_r,τ,x_s)|τ∈[-T,T]}

对互相关连解方程对时间求导数，建立隐函数方程：

其中A(f₀,x_r,x_s)_obs表示观测数据的最大值，A(f₀,x_r,x_s)_cal表示模拟数据的最大值。则目标函数对速度的导数可以表示为：

其中本发明只讨论常密度声波方程，可以表示为：

其中u代表正传波场，s代表震源函数。将声波方程写成离散的矩阵形式：

L[u(f₀,t,x_s)]＝s

表示正演算子，正演算子对速度参数的导数可以表示为:

则地震波场对速度模型的导数可以表示为:

其中u表示地震波场的模拟数据，将公式代入得到模拟数据对速度模型的偏导数：

其中(L^-1)^T表示伴随算子，用于将伴随震源反传到模型空间，对于二阶声波方程如公式(11)存在(L^-1)^T＝(L^-1)。T表示矩阵的转置。根据波场对速度的偏导数，可以得出走时信息对速度的偏导数：

对于单频波形只存在走时上的差别，因此存在：

不依赖震源子波的时间域单频波形走时反演方法的梯度可以表示为：

则伪伴随震源的表达式可以表示为：

多炮地震数据的梯度可以表示：

其中u(f₀,t,x_s)_cal表示单频正传波场，u(f₀,x_r,t,x_s)_cal单频模拟数据，f₀表示反演频率，表示单频模拟数据对时间的一阶导数，d(f₀,x_r,t+Δτ,x_s)_obs表示单频观测数据，表示单频观测数据对时间的一阶导数，(L^-1)^T表示对反传算子的转置；

l、计算伪伴随震源，并利用反传算子将伪伴随震源反传到模型空间，得到反传波场；

m、正传波场与反传残差波场做互相关得到模型更新梯度；

n、利用L-BFGS优化算法计算模型更新方向，并通过Wolfe收敛准则寻找最合适的步长，按照公式v_k+1＝v_k-α_kH_kg_k对速度模型进行更新，其中k表示迭代的次数，H_k表示当前步骤伪海森矩阵的逆，α_k表示当前步骤的更新步长，g_k表示当前步骤的更新梯度，v_k表示当前步骤的速度模型，v_k+1表示经过本次迭代以后的速度模型。

o、判断是够满足终止条件，若满足输出不依赖震源子波的时间域单频波形走时反演结果。若不满足终止条件，利用循环频率多尺度策略逐渐提高频率，将反演结果作为下一个循环的初始速度模型，返回e步骤。

p、将上一步的输出结果作为常规时间域全波形反演的初始速度模型，然后进行常规时间域全波形反演，并输出最终反演结果。

实施例1

根据勘探要求，将MATLAB2013a在Windows 7旗舰版系统下进行安装，并安装基于MATLAB语言编写的地震数据处理软件Crews工具包。

本发明利用Marmousi进行测试，Marmousi具有复杂的构造、细微的层理和深部储油藏构造，很适合进行理论方法测试。原始Marmousi模型较大，考虑到硬件内存以及计算时间因素，本发明对原始模型进行抽稀处理，利用抽稀之后的Marmousi模型进行不依赖震源子波的时间域单频波形走时反演方法测试。真实模型(图12左)和线性递增初始模型(图12右)。

模型参数如下：

表1：不依赖震源子波的时间域单频波形走时反演方法测试参数(时间域)

模型大小	网格间距	横向距离	纵向深度	速度范围	地震子波主频
						69*384	12.5m	2.4km	0.8625km	1.5-4km/s	22Hz

时间域频率多尺度策略参数如下：

对观测数据和模拟数据进行抽频处理，并得到对应的单频波形。低频段对应着模型的宏观信息，高频段对应着模型的细节信息，只有先得到模型的宏观结构以后，才能在宏观结构的框架上得到模型准确的细节信息。本发明采用的频率范围为freq＝15:0.5:25Hz进行反演。从低频到高频不断增加反演频率来达到高精度的反演结果。每一个频率迭代30次，上一个频率得到的反演结果用于下一个频率的初始速度模型，同时按照图11所示的示意图来进行选择频率，“Z”字型交错向前反演。时间域单频波形波动方程走时反演方法，由宏观到细微的多尺度反演策略，极大地减小了每一次反演过程中的非线性程度，有效避免了周跳现象，同时增加了反演程序的稳定性。

基于以上参数，在表2所示的测试环境下进行不依赖震源子波的时间域单频波形走时反演方法。

表2性能测试环境

从图14中可以看出，Marmousi速度模型的反演结果左侧出现了明显的错误，这是由于地震数据的波形匹配过程中出现了错位的现象，导致这一主要原因是地震数据中波形变化复杂，且不存在低频成分。从图14(右)可以看出，当地震数据中缺失低频成分时，即使利用低通滤波多尺度策略仍然不能有效的解决全波形反演过程中的周波跳跃现象。图15(左)是利用单频波形波动方程走时反演反演所得到的，在反演的过程中利用波形的走时信息，相比于利用波形信息匹配具有更弱的非线性关系，所以只得到了Marmousi速度模型的宏观结构，利用这个速度模型作为初始速度模型，进行常规全波形反演得到图15(右)所示结果，反演精度较常规多尺度全波形反演高了很多，同时基本不存在反演错误的现象。证明了时间域单频波形走时反演方法不依赖于初始速度模型，即使在缺失低频成分的数据中，也能得到高精度的反演结果。

从图16(左)可以看出，当震源子波不准确的时候对反演结果有着很严重的影响，甚至得不到地下的任何信息。相比于图16(右)利用了相位校正方法，很好消除了震源子波不准确对单频波形波动方程走时反演结果的影响，同时证明了基于相位校正方法在消除震源子波方面的可靠性。

图1是整个反演过程的流程图，从流程图可以看出首先通过相位校正技术来消除震源子波对单频波形波动方程走时反演的影响，然后再结合单频多尺度的优势，从低频逐渐向高频进行反演，最后将不依赖震源子波的时间域单频波形走时反演结果作为常规时间域全波形反演的初始模型，得到高精度的反演结果。证明了本发明提出的不依赖震源子波的时间域单频波形走时反演在缓解周波跳跃问题上有着一定的优势。

Claims (1)

1.一种不依赖震源子波的时间域单频波形走时反演方法，其特征在于，包括以下步骤：

a、安装MATLAB2013a正版软件，并安装基于MATLAB语言编写的地震数据处理软件Crews工具包；

b、对实际采集的地震数据进行预处理，并将观测数据和观测系统输入计算机准备进行不依赖震源子波的时间域单频波形走时反演；

c、对观测数据进行傅里叶变换，在频率域挑选21个频率数据进行反演，从15Hz到25Hz，先利用低频的数据，再逐渐增加频率以实现多尺度的反演策略；

d、建立线性递增初始速度模型，该速度模型作为反演的起始值，通过计算模型的更新方向，并加到初始速度模型上，来实现由线性初始速度模型逐渐向真实模型逼近的过程；

e、利用输入的观测系统、震源子波和线性初始速度模型进行波动方程正演模拟，得到模拟数据，并提取模拟数据震源附近直达波波形，用于震源波形相位校正，消除震源子波对全波形反演的影响；

f、将提取的观测数据直达波波形和模拟数据直达波波形进行傅里叶变换，并挑选对应的21个频率数据，求取这两个数据对应频率的相位差，并将求得的相位差用于单频观测数据的相位校正中；

g、对相位校正以后的单频观测数据进行傅里叶逆变换，得到时间域的单频观测数据；

h、对模拟数据和单频震源子波进行傅里叶变换，并提取对应21个频率数据，再对单频模拟数据和震源子波进行傅里叶逆变换，得到单频的时间域模拟数据和单频时间域震源子波；

i、利用单频时间域震源子波和速度模型进行波动方程有限差分正演模拟，得到单频的正传波场，并对正传波场进行存储，用于后续的梯度计算；

j、利用互相关方法计算单频观测数据和单频模拟数据的走时差；

k、根据最小二乘原理构造走时差目标函数：其中Δτ表示单频观测数据与单频模拟数据存在的走时差，v表示速度值，s表示震源数目，r表示检波器数目，x_s表示震源位置，x_r表示检波器位置，在求目标函数的梯度过程中，对目标函数两端关于速度v求导，得到梯度表达式：

则伪伴随震源的表达式可以表示为：

其中u(f₀,t,x_s)_cal表示单频正传波场，u(f₀,x_r,t,x_s)_cal表示单频模拟数据，f₀表示反演频率，表示单频模拟数据对时间的一阶导数，d(f₀,x_r,t+Δτ,x_s)_obs表示单频观测数据，表示单频观测数据对时间的一阶导数，(L^-1)^T表示对反传算子的转置；

l、计算伪伴随震源，并利用反传算子将伪伴随震源反传到模型空间，得到反传波场；

m、正传波场与反传残差波场做互相关得到模型更新梯度；

n、利用L-BFGS优化算法计算模型更新方向，并通过Wolfe收敛准则寻找最合适的步长，按照公式v_k+1＝v_k-α_kH_kg_k对速度模型进行更新，其中k表示迭代的次数，H_k表示当前步骤伪海森矩阵的逆，α_k表示当前步骤的更新步长，g_k表示当前步骤的更新梯度，v_k表示当前步骤的速度模型，v_k+1表示经过本次迭代以后的速度模型；

o、判断是否满足终止条件，若满足则不依赖震源子波的时间域单频波形走时反演结果，若不满足终止条件，利用循环频率多尺度策略逐渐提高频率，将反演结果作为下一个循环的初始速度模型，返回e步骤；

p、将上一步的输出结果作为常规时间域全波形反演的初始速度模型，然后进行常规时间域全波形反演，并输出最终反演结果。