地球物理学报  2016, Vol. 59 Issue (10): 3765-3776   PDF    
引用本文
陈生昌, 陈国新. 2016. 时间二阶积分波场的全波形反演. 地球物理学报, 59(10): 3765-3776, doi: 10.6038/cjg20161021.  
CHEN S C, CHEN G X. 2016. Full waveform inversion of the second-order time integral wavefield. Chinese J. Geophys. (in Chinese), 59(10): 3765-3776, doi: 10.6038/cjg20161021.  
时间二阶积分波场的全波形反演
陈生昌 , 陈国新     
浙江大学地球科学学院, 杭州 310027
收稿日期 2015-11-10, 2016-03-18 收修定稿
基金项目: 国家自然科学基金项目（41074133，41374001）资助
作者简介: 陈生昌, 男, 1965年生, 教授, 博士生导师.主要从事勘探地球物理和计算地球物理研究.E-mail:chenshengc@zju.edu.cn
摘要: 通过对波场的时间二阶积分运算以增强地震数据中的低频成分，提出了一种可有效减小对初始速度模型依赖性的地震数据全波形反演方法-时间二阶积分波场的全波形反演方法.根据散射理论中的散射波场传播方程，推导出时间二阶积分散射波场的传播方程，再利用一阶Born近似对时间二阶积分散射波场传播方程进行线性化.在时间二阶积分散射波场传播方程的基础上，利用散射波场反演地下散射源分布，再利用波场模拟的方法构建地下入射波场，然后根据时间二阶积分散射波场线性传播方程中散射波场与入射波场、速度扰动间的线性关系，应用类似偏移成像的公式得到速度扰动的估计，以此建立时间二阶积分波场的全波形迭代反演方法.最后把时间二阶积分波场的全波形反演结果作为常规全波形反演的初始模型可有效地减小地震波场全波形反演对初始模型的依赖性.应用于Marmousi模型的全频带合成数据和缺失4 Hz以下频谱成分的缺低频合成数据验证所提出的全波形反演方法的正确性和有效性，数值试验显示缺失4 Hz以下频谱成分数据的反演结果与全频带数据的反演结果没有明显差异.
关键词: 散射波方程      时间二阶积分      增强低频      Born近似      全波形反演     
Full waveform inversion of the second-order time integral wavefield
CHEN Sheng-Chang, CHEN Guo-Xin     
School of Earth Sciences, Zhejiang University, Hangzhou 310027, China
Abstract: We proposed a new full waveform inversion (FWI) method, namely full waveform inversion of the second-order time integral wavefield, by enhancing low-frequency components of seismic data with a second-order time integration to seismic wavefield, which can efficiently reduce the initial model dependence of FWI. According to the propagation equation of scattering wavefield in scattering theory, we derived a propagation equation for the scattering wavefield with second-order time integral, and used the leading order Born approximation for the linearization to the propagation equation. Based on this equation, using the scattering wavefield to inverse the distribution of scattering sources in subsurface, and using wavefield modeling to construct the incident wavefield, and according to the linear relationship between the scattering wavefield and incident wavefield, velocity perturbation in the linear propagation equation for the scattering wavefield with second-order time integral, we applied a formula similar to the imaging formula of migration to estimate velocity perturbation, and established an iterative inversion method for the full waveform inversion of second-order integral wavefield. Applying the inversion result of full waveform inversion of second-order integral wavefield as the initial velocity model for the conventional FWI can efficiently reduce the initial model dependence of FWI. Numerical tests using synthetic data of the Marmousi model demonstrate the validity and feasibility of the proposed method. The final results of the new method can obtain much improved results than the conventional FWI. Furthermore, to test the independence to the seismic frequency-band, we used a low-cut source wavelet (cutting from 4Hz below) to generate the synthetic data. The inversion results by our new method show no appreciable difference from the full-band source results..
Key words: Propagation equation of scattering wave      Second-order time integral      Enhancement low-frequency components      Born approximation      Full waveform inversion     
1 引言

利用地震数据获取地下速度分布是地震勘探中的一个经典问题，同时又是地震数据反演方法理论研究中的难点.地震波场与地下速度模型之间是一种非线性的函数关系，利用观测的地震波场反演地下速度模型的分布属于典型的非线性反演问题.对于非线性反演问题的求解通常需要给定一个初始速度模型，由给定的速度模型和震源函数以及相应的初边值条件可计算出与速度模型对应的计算波场，进而获得观测波场与计算波场间的残差波场.如果残差波场小于设定的阈值，则认为给定的速度模型就是利用地震波场反演出的地下速度模型.如果残差波场大于设定的阈值，则需要对给定的速度模型进行修正.根据残差波场信息对初始速度模型进行修正主要两类方法，1)是非线性最优化方法；2)是线性化迭代反演方法.第1类的非线性最优化方法又包含全局最优化方法和局部最优化方法，由于得到计算波场的波场模拟计算量巨大，通常采用局部最优化方法对速度模型进行修正，这也就是目前广泛研究和应用的全波形反演方法.基于局部最优化方法的全波形反演自上世纪80年代初提出以来(Lailly，1983；Tarantola, 1984, 1986, 1987)，就一直是地震波场反演研究的热点，并逐渐成为了地震勘探中的地下速度建模方法(Virieux and Operto, 2009).本文把基于局部最优化方法的全波形反演称为常规全波形反演.第二类的全波形反演的线性化迭代反演方法主要基于地震波传播的散射理论(Chew，1990)，利用渐近展开的方法(如Taylor展开、一阶Born近似、Raytov近似等方法)把非线性的全波形反演问题线性化为迭代的线性反演问题，即在迭代中利用线性反演方法求解初始速度模型的修正.在数据拟合的最小平方准则下，基于局部最优化方法的全波形反演与基于线性化迭代反演方法的全波形反演是等价的(Vogel，2002).

常规全波形反演对初始速度模型具有很强的依赖性，这也是利用局部优化方法求解非线性反演问题的内在局限性，为了减小这种依赖性，人们提出了各种多尺度多层次的反演策略，如Bunks等提出的由低频到高频不同子集数据的多尺度时间域反演策略(Bunks et al., 1995)，Pratt等提出的在频率域对地震数据由低频到高频进行分频反演策略(Pratt et al., 1998; Pratt, 1999).通过对常规全波形反演目标函数的分析，认为目标函数在数据低频段具有比在高频段更大的凸性范围，全波形反演中起负作用的周期跳现象在低频段也可得到有效地压制，得到了低频数据的全波形反演对初始模型的要求相对较低的认识(Sirgue, 2006；Sirgue and Pratt, 2004).基于上述认识，应用各种不同数学变换挖掘地震数据中的低频信息，并将得到的低频信息用于全波形反演中的初始速度模型反演引起了不少研究者的重视，并提出了一些相应的全波形反演方法.Shin等提出的Laplace域、Lapace-Fourier域、对数域的全波形反演方法(Shin and Cha, 2008，2009；Shin and Min, 2006；Shin et al., 2007)，Donno等提出的归一化积分法全波形反演方法(Chauris et al., 2012；Donno et al., 2013)，Bozda等及其他研究者提出的包络全波形反演方法(Bozda et al., 2011；Wu et al., 2014；Chi et al., 2013, 2014)，上述这些反演方法都是通过对数据的数学变换来提取地震数据中的低频信息以反演速度模型中的长波长成分.这些方法虽然在全波形反演中取得了一定的效果，但是它们所用到的挖掘低频信息的方法多来自于信号变换，与描述地震波场传播的数学物理方程没有直接的关系，因此所得到的所谓低频信息的数学物理含义不明确.

常规全波形反演方法仅考虑地震波场与地下速度模型间的非线性关系，而没有考虑地震波场传播的物理过程，因此常规的全波形反演方法可视为一种纯粹基于计算数学而发展起来的方法.在地震波场的传播中，散射波场是最基本的波场，其它如反射、折射波场等是散射波场在高频近似下的退化.地震散射波场传播的物理过程包括，1)震源激发的入射波场在背景速度模型中传播；2)入射波场与相对于背景速度的速度扰动的相互作用产生散射波场，散射波场再与速度扰动的相互作用产生高次散射波场；3)散射波场在背景速度模型中传播.如果把常规全波形反演中的初始速度模型视为散射理论中散射波场传播方程的背景速度模型，并对散射波场的传播方程进行线性化，则上述的散射波场的产生与传播过程可为地震波场的全波形反演方法研究提供思路.

基于上述的认识，本文根据散射理论中散射波场的传播方程，首先推导出时间二阶积分散射波场的传播方程，然后把一阶Born近似条件下的线性化迭代反演方法与散射波场的传播过程相结合，构建基于线性化迭代方式的时间二阶积分散射波场的全波形反演方法.对波场的时间二阶积分运算相当于一种低通滤波可以增强数据中已有的低频成分，不同于Laplace、Laplace-Fourier域全波形反演方法、归一化积分法和包络反演方法中利用各种不同变换来提取地震数据中的低频信息，利用时间二阶积分增强地震数据中低频成分的方法是来自于散射波场的传播方程，它具有明确的数学物理含义.在时间二阶积分散射波场全波形反演中，把常规全波形反演中的初始速度模型视为散射波场传播中的背景速度模型，把观测波场与由初始速度模型合成得到的计算波场之间的残差波场视为速度扰动产生的散射波场；在推导出的时间二阶积分散射波场传播方程的基础上，利用逆传播或伴随传播的方法由地表的时间二阶积分残差波场重建出地下散射源的分布，同时利用震源波场正向传播的方法获取地下入射波场分布，并根据时间二阶积分散射波场线性传播方程中散射波场与入射波场、速度扰动间的线性关系，再利用类似偏移中的成像公式获得初始速度模型的修正.最后再把时间二阶积分散射波场全波形反演得到的结果作为常规全波形反演的初始速度模型，对地震波场进行常规的全波形反演.通过Marmousi模型的全频带合成数据和缺低频合成数据的数值试验验证本文所提出的全波形反演方法的正确性和有效性.

2 方法原理

众所周知，正演是反演的基础，因此我们在这节首先推导出时间二阶积分散射波场的传播方程，并分析时间二阶积分波场的含义，然后再构建基于时间二阶积分散射波场传播方程的全波形反演方法.

2.1 时间二阶积分散射波场的传播方程

本文考虑单P波速度的全波形反演.给定初始(背景)速度模型v0(x)和震源函数s(t)，入射波场ui(x, t)的传播方程有

(1-1)

式中，xs为震源点位置；▽2为Laplace算子，有；ui(x, t)也称为与背景速度模型对应的背景波场；s(t)为震源函数.与方程(1-1)相应的Green函数g(x, x′, t, t′)，有

(1-2)

把速度模型表示为v(x)=v0(x)+δv(x)，δv(x)为速度扰动，也称为初始速度模型v0(x)的修正.根据散射理论，散射波场δu(x, t)的传播方程有

(2)

方程(2)右边的波场u(x, t)为与速度模型v(x)对应的总波场；δu(x, t)也称为与扰动速度模型δv(x)有关的扰动波场.下文把方程(2)中的近似等号写为等号，把方程(2)变换到频率域，有

(3)

(4)

则方程(3)可写为

(5-1)

把方程(5-1)变换到时间域，有

(5-2)

根据表达式(4)(对该式的详细解释在下一小节)可知，波场δui(x, t)是散射波场δu(x, t)的时间二阶积分，我们称之为时间二阶积分散射波场，方程(5-1)、(5-2)分别为时间二阶积分散射波场传播方程的频率域形式和时间域形式.

把一阶Born近似条件应用于方程(5-2)，得到

(6)

方程(6)就是一阶Born近似条件下的时间二阶积分散射波场δui(x, t)的线性传播方程，它也是本文全波形反演方法研究的基础方程.在方程(6)中，δui(x, t)与速度扰动δv(x)和入射波场ui(x, t)都成线性关系.为了后述公式推导的方便，利用方程(1-2)中的Green函数，把方程(6)写为积分形式，有

(7)

式中，T0为波场的最大记录时间；Ω为速度扰动δv(x)的分布空间.令，

(8)

m(x, t)可视为产生时间二阶积分散射波场的散射源，它是与δv(x)和ui(x, t)的分布、强弱有关的.利用式(8)，公式(7)可写为

(9)

式(9)的积分是时间-空间域的褶积.如果已知左边的δui(x, t)求右边积分式中的m(x′, t′)，则意味求解一个第一类Fredholm线性积分方程.对式(9)两边做时间Fourier变换，有

(10)

式中，和分别为g(x, x′, t, t′)和m(x′, t′)的时间Fourier变换.对于有

(11)

2.2 波场的时间二阶积分

我们知道，时间域的一阶微分运算，对应于频率域的乘积运算(-iω)，积分运算是微分运算的逆运算，(-iω)的逆运算对应于时间域的一阶积分运算∫dt，因此公式(4)所定义的运算实质上就是波场的时间域二阶积分在频率域的实现，即

(12)

积分运算相当于一种低通滤波，因此对波场进行式(12)所表示的频率域运算有助于增强波场的低频成分.把式(12)进行时间Fourier反变换，有

(13)

F-1表示时间Fourier反变换.因此时间二阶积分波场ui(t)相对于原始波场u(t)其低频成分得到了有效增强.

假定波场u(t)为主频等于15 Hz的Ricker子波，如图 1中的蓝线所示，对应的时间二阶积分波场ui(t)如图 1中的红线所示，它们的振幅谱如图 2所示，图 2中的蓝线为Ricker子波的振幅谱，红线为Ricker子波时间二阶积分的振幅谱.由图 1的时间信号和图 2的振幅谱可看出，Ricker子波的二阶积分信号相对于Ricker子波表现出明显的低频化，通过二阶积分运算使Ricker子波中的低频成分得到了增强，二阶积分信号的能量主要集中在低频段，强化了低频信息在数据中的作用.

图 1 主频为15 Hz的Ricker子波及其时间二阶积分：蓝线为Ricker子波；红线为Ricker子波的时间二阶积分 Fig. 1 The ricker wavelet with peak-frequency 15Hz and its second-order time integral. The blue-line denotes ricker wavelet. The red-line denotes the second-order time integral of ricker wavelet
图 2 主频为15 Hz的Ricker子波及其时间二阶积分的振幅谱；蓝线为Ricker子波的振幅谱，红线为Ricker子波时间二阶积分的振幅谱 Fig. 2 The amplitude spectra of ricker wavelet with peak-frequency 15 Hz and its second-order time integral. The blue-line denotes the amplitude spectra of ricker wavelet. The red-line denotes the amplitude spectra of second-order time integral of ricker wavelet

图 3为Marmousi模型的单炮合成记录，图 4为对图 3中的地震波场进行时间二阶积分运算得到的时间二阶积分波场，对比图 3、4可看出，图 4中波场的波形变化相对于图 3中的波形表现缓慢，图 4中的波场和图 3中的波场相比表现出明显的低频化.图 5为图 3中波场的包络信号图，由于信号的包络处理和所得到的包络信号的数学物理意义不明确，图 5所示的包络信号仅是一种能量集中在浅部、变化非常缓慢的非负值信号.图 4中的时间二阶积分波场信号是受波场传播方程控制，是一种具有明确数学物理意义的地震波场.

图 3 Marmousi模型的单炮合成记录 Fig. 3 The synthetic single-shot record of Marmousi model
图 4 图 3所示的单炮记录的时间二阶积分 Fig. 4 The second-order time integral of the synthetic single-shot record in Fig.3
图 5 图 3所示的单炮记录的包络 Fig. 5 The envelope of the synthetic single-shot record in Fig.3
2.3 时间二阶积分波场的全波形反演

由于地震波场是速度模型的非线性函数，利用观测的地震波场数据uobs(x, t)反演地下速度模型v(x)属于非线性反演问题，本文采用线性化迭代反演方法进行求解.

给定一个初始速度模型v0(x)，并计算出与之对应的计算波场ucal(x, t)，把观测波场uobs(x, t)与

ucal(x, t)间的残差波场，即

(14)

视为由地下真实速度模型与给定速度模型之间速度差异，δv(x)=v(x)-v0(x)，引起的散射波场.计算散射波场δu(x, t)的频率域、时间域二阶积分散射波场，即

(15)

(16)

式(15)中的F表示时间Fourier变换.

根据入射波场的传播方程和前面推导出的时间二阶积分散射波场线性传播方程以及公式(8)或(11)，我们就可以在时间域利用δui(x, t)或在频率域利用反演初始速度模型v0(x)的修正量δv(x).其过程为，先利用已知的δui(x, t)或，通过求解线性积分方程(9)或(10)，得到时间域的m(x, t)或频率域的，结合通过波场模拟得到的地下入射波场，再利用公式(8)或(11)得到δv(x)的估计，即

(17)

或

(18)

式中，Re表示取实部运算；符号*表示复共轭运算；ui(x, t)和分别表示时间域和频率域的地下入射波场.为了计算的稳定性，利用类似于偏移成像公式的处理方式，公式(17)和(18)可分别写为，

(19)

(20)

得到了速度扰动δv(x)的解估计后，我们就可以对给定的初始速度模型v0(x)进行修正，得到新的速度模型v1(x)，即v1(x)=v0(x)+αδv(x)，其中α为修正步长，然后再把v1(x)作为下一步反演的初始速度模型，进行迭代反演.

由上述的反演步骤、公式可知，利用二阶积分散射波场和线性积分方程(9)或(10)求解出时间域的m(x, t)或频率域的是关键.把方程(9)写为算子形式，有

(21)

利用最小二乘方法得到算子方程(21)的最小二乘解，即

(22)

式中G-1为散射波场传播算子G的最小二乘逆算子，有

(23)

式(22)所表示的运算为在时间-空间域对波场δui进行逆传播.由于这种逆传播不仅计算量巨大，而且还很不稳定，常用G的伴随算子Gt代替逆算子，即G-1≈Gt，因此式(22)改写为

(24)

式(24)所表示的运算为在时间-空间域对波场δui进行伴随(逆时)传播.这种伴随传播不仅计算量较小，而且还稳定.

把式(23)代入式(22)，再把式(22)代入式(17)就可得到速度扰动的解估计，并把各个共炮道集数据的反演结果叠加，有

(25)

也可以把式(23)代入式(22)，再把式(22)代入式(19)得到速度扰动的解估计，并把各个共炮道集数据的反演结果叠加，有

(25-1)

把式(24)代入式(19)，同样也把各个共炮道集数据的反演结果叠加，有

(26)

在上述式中，xs表示共炮道集的炮点位置；ui(xs, x, t)表示初始速度模型下正向传播的震源波场.式(25)、(25-1)、(26)为在时间域求解δv(x)的公式.

上述的时间域反演工作也同样可在频率域实现，把频率域的方程(10)写为算子形式，有

(27)

同样利用最小二乘方法得到算子方程(27)的最小二乘解，即

(28)

式中为散射波场传播算子的最小二乘逆算子，有

(29)

式(28)所表示的运算为在频率-空间域对波场进行逆传播.同样由于这种逆传播不仅计算量巨大，而且还很不稳定，常用的伴随算子代替逆算子，即，因此式(28)改写为

(30)

式(30)所表示的运算为在频率-空间域对波场进行伴随(逆时)传播.这种伴随传播不仅计算量较小，而且还稳定.

把式(29)代入式(28)，再把式(28)代入式(18)就可得到速度扰动的解估计，并把各个共炮道集数据的反演结果叠加，有

(31)

把式(29)代入式(28)，再把式(28)代入式(20)得到速度扰动的解估计，并把各个共炮道集数据的反演结果叠加，有

(31-1)

把式(30)代入式(20)，同样也把各个共炮道集数据的反演结果叠加，有

(32)

在上述式中表示背景速度模型下正向传播的震源波场.式(31)、(31-1)、(32)为在频率域求解δv(x)的公式.

把正向传播的震源波场和二阶积分波场逆传播的具体形式代入式(31)、(31-1)，分别有，

(33)

(33-1)

把正向传播的震源波场和二阶积分波场逆时传播的具体形式代入式(32)，有，

(34)

在上述式中为震源函数s(t)的频谱.

利用上述的时间二阶积分波场的全波形反演虽然可以有效地减小对初始速度模型的依赖性，但时间二阶积分波场相对于原始波场存在明显的低频化，因此时间二阶积分波场全波形反演得到的反演结果分辨率受到一定的限制，为了充分地利用地震数据中的全部信息，以得到高分辨率的全波形反演结果，我们把时间二阶积分波场全波形反演得到的反演结果作为常规全波形反演的初始速度模型，再进行常规的全波形反演.

3 数值试验

为了验证本文所提出的时间二阶积分波场全波形反演(英文简写为IntI)的正确性和有效性，我们用Marmousi模型的合成地震数据进行试验.试验中Marmousi模型的网格化参数：nx=369，nz=125，dx=25 m，dz=24 m.图 6为Marmousi模型的速度分布图.试验用的合成地震数据共有92炮，炮间距100 m，每炮369道，道间距25 m.时间采样长度5 s，采样间隔2 ms.震源为主频8 Hz的Ricker子波.我们合成了2套地震数据用于试验，一套是全频带地震数据；另一套是完全缺失4 Hz以下频率成分，并以4~5 Hz为过渡带的缺低频地震数据.图 7为试验用的常梯度初始速度模型.

图 6 Marmousi模型的速度分布 Fig. 6 The velocity distribution of Marmousi model
图 7 用于反演的常梯度初始速度模型 Fig. 7 The initial velocity model with constant gradient for inversion
3.1 全频带地震数据试验

对于合成得到的全频带地震数据全波形反演试验，我们首先用图 7所示的常梯度速度模型作为时间二阶积分波场全波形反演的初始速度模型，在反演中应用时间域公式(26)进行速度修正量的计算.图 8为应用时间二阶积分波场全波形反演得到的反演结果，对比图 6、8可看出，由时间二阶积分波场反演得到的速度模型大体结构正确，但分辨率较低.我们再把图 8所示的速度模型作为常规全波形反演(英文简写为FWI)的初始速度模型，图 9为用图 8作初始速度模型的常规全波形反演结果.作为比较，我们还用图 7所示的常梯度速度模型作为常规全波形反演的初始速度模型，图 10为用图 7作初始速度模型的常规全波形反演结果.对比图 6、9、10可看出，对于同样的常梯度初始速度模型，时间二阶积分波场全波形反演加常规全波形反演(IntI+FWI)的反演结果要明显优于单纯的常规全波形反演的反演结果.为了更仔细对比上述两种反演结果的细节，图 11为Marmousi模型中第85、235网格点两种反演结果的速度-深度曲线对比图，图中绿线为真实速度，黑线为初始速度，红线为常规全波形反演的结果，蓝线为时间二阶积分波场全波形反演加常规全波形反演的结果.由图 11中的速度-深度曲线可看出，时间二阶积分波场全波形反演加常规全波形反演能很好地重建出真实速度模型.

图 8 时间二阶积分波场全波形反演(IntI)的反演结果 Fig. 8 The inversion result by full waveform inversion of the second-order time integral wavefield
图 9 时间二阶积分波场全波形反演加常规全波形反演(IntI+FWI)的反演结果 Fig. 9 The inversion result by full waveform inversion of the second-order time integral wavefield plus conventional full waveform inversion
图 10 常规全波形反演(FWI)的反演结果 Fig. 10 The inversion result by conventional full waveform inversion
图 11 第85、235网格点两种反演结果的速度-深度曲线对比：上图为第85网格点；下图为235网格点.图中绿线为真实速度，黑线为初始速度，红线为常规全波形反演的结果，蓝线为时间二阶积分波场全波形反演加常规全波形反演的结果 Fig. 11 Comparisons of velocity-depth curves of the two inversion results at 85th grid-point and 235th grid-point: the up-panel is for 85th grid-point; the bottom-panel is for 235th grid-point. The green line denotes the true velocity, the black line denotes the initial velocity. The red line denotes the velocity from conventional FWI. The blue line denotes the velocity from full waveform inversion of the second-order time integral wavefield plus conventional FWI
3.2 缺低频地震数据试验

为了进一步验证本文所提出的全波形反演方法对于缺低频地震数据的反演效果，我们进行了对完全缺失4 Hz以下频率成分，并以4~5 Hz为过渡带的缺低频地震数据全波形反演试验.在试验中，我们首先用图 7所示的常梯度速度模型作为时间二阶积分波场全波形反演的初始速度模型，在反演中仍然应用时间域公式(26)进行速度修正量的计算.图 12为应用时间二阶积分波场全波形反演得到的反演结果.我们再把图 12所示的速度模型作为常规全波形反演的初始速度模型，图 13为用图 12作初始速度模型的常规全波形反演结果.作为比较，我们还用图 7所示的常梯度速度模型作为初始速度模型对缺低频地震数据进行了单纯的常规全波形反演和包络反演加常规全波形反演(EI+FWI).图 14为用图 7作初始速度模型的常规全波形反演结果，图 15为用图 7作初始速度模型的包络反演加常规全波形反演的反演结果.

图 12 时间二阶积分波场全波形反演(IntI)的反演结果 Fig. 12 The inversion result by full waveform inversion of the second-order time integral wavefield
图 13 时间二阶积分波场全波形反演加常规全波形反演(IntI+FWI)的反演结果 Fig. 13 The inversion result by full waveform inversion of the second-order time integral wavefield plus conventional full waveform inversion
图 14 常规全波形反演(FWI)的反演结果 Fig. 14 The inversion result by conventional full waveform inversion
图 15 包络反演加常规全波形反演(EI+FWI)的反演结果 Fig. 15 The inversion result by envelope inversion plus conventional full waveform inversion

对比图 6、13、14、15可看出，对于同样的常梯度初始速度模型，时间二阶积分波场全波形反演加常规全波形反演的反演结果要明显优于单纯的常规全波形反演的反演结果和包络反演加常规全波形反演的反演结果，包络反演加常规全波形反演的反演结果次之，单纯的常规全波形反演的反演结果最差.为了更仔细地对比上述三种反演结果的细节，图 16为Marmousi模型中第85、235网格点三种反演结果的速度-深度曲线对比图，图中绿线为真实速度，黑线为初始速度，深红线为常规全波形反演的结果，浅红线为包络反演加常规全波形反演的结果，蓝线为时间二阶积分波场全波形反演加常规全波形反演的结果.由图 16中的速度-深度曲线可看出，即使对于缺低频的地震数据，时间二阶积分波场全波形反演加常规全波形反演也能很好地重建出真实速度模型.

图 16 第85、235网格点三种反演结果的速度-深度曲线对比：上图为第85网格点；下图为235网格点.图中绿线为真实速度，黑线为初始速度，深红线为常规全波形反演的结果，浅红线为包络反演加常规全波形反演的结果，蓝线为时间二阶积分波场全波形反演加常规全波形反演的结果 Fig. 16 Comparisons of velocity-depth curves of the three inversion results at 85th grid-point and 235th grid-point: the up-panel is for 85th grid-point; the bottom-panel is for 235th grid-point. The green line denotes the true velocity. The black line denotes the initial velocity, the dark-red line denotes the velocity from conventional FWI. The light-red line denotes the velocity from envelope inversion plus conventional FWI. The blue line denotes the velocity from full waveform inversion of the second-order time integral wavefield plus conventional FWI

由全频带地震数据的时间二阶积分波场全波形反演得到的结果图 8与缺低频地震数据的时间二阶积分波场全波形反演得到的结果图 12的对比，以及全频带地震数据的时间二阶积分波场全波形反演加常规全波形反演得到的结果图 9与缺低频地震数据的时间二阶积分波场全波形反演加常规全波形反演得到的结果图 13的对比，我们可看出对于缺低频地震数据应用时间二阶积分波场全波形反演方法进行反演可得到与全频带地震数据应用时间二阶积分波场全波形反演方法进行反演基本一致的反演结果，这表明时间二阶积分波场全波形反演方法对缺失低频地震数据的有效性.图 17为反演结果图 9、13中第85、235网格点的速度-深度曲线对比图，图中绿线为真实速度，黑线为初始速度，红线为缺低频地震数据的反演结果，蓝线为全频带地震数据的反演结果，由图 17的对比可看出缺失4 Hz以下频谱成分数据的反演结果与全频带数据的反演结果没有明显差异.

图 17 第85、235网格点全频带地震数据与缺低频地震数据反演结果的速度-深度曲线对比：上图为第85网格点；下图为235网格点.图中绿线为真实速度，黑线为初始速度，红线为缺低频地震数据的反演结果，蓝线为全频带地震数据的反演结果 Fig. 17 Comparisons of velocity-depth curves of inversion result by full band seismic data and inversion result by lack low frequencies seismic data at 85th grid-point and 235th grid-point. The up-panel is for 85th grid-point; the bottom-panel is for 235th grid-point. The green line denotes the true velocity, the black line denotes the initial velocity, the red line denotes the velocity frominversion result by lack low frequencies seismic. The blue line denotes the velocity from inversion result by full band seismic data

同时给出时间域和频率域的公式推导，主要是便于对本文方法的理解和与现有一些FWI公式的对比.本文反演计算是在时间-空间域进行的，由给出的公式可知同样可以在频率-空间域进行反演计算.

4 结论

基于局部极值问题和周期跳现象在低频地震数据全波形反演中相对不敏感的认识，由散射波场传播方程推导出了时间二阶积分散射波场的传播方程，以及一阶Born近似下的时间二阶积分散射波场的线性传播方程.把推导出的时间二阶积分散射波场线性传播方程应用于非线性反演问题的线性化迭代反演，并结合时间二阶积分散射波场线性传播方程中散射波场与入射波场和速度扰动间的线性关系，建立了一种首先反演地下散射源分布，然后再求解速度扰动的时间二阶积分波场全波形反演方法.对地震波场的时间二阶积分可有效地增强波场中已有的低频成分，强化了地震数据中的低频信息在反演中的作用，这是本文方法与由低频到高频的分频反演方法的不同之处，也是其长处，使本文提出的时间二阶积分波场全波形反演方法对初始模型的依赖性减弱.不同于包络全波形反演和归一化积分全波形反演等方法中的数据变换方法，本文的波场时间二阶积分运算是来自于描述散射波场传播的数学物理方程，因此得到的时间二阶积分波场具有清楚的数学物理含义.把时间二阶积分波场全波形反演方法与常规全波形反演方法相结合，不仅可以得到高分率的反演结果，还能在缺低频地震数据的全波形反演中发挥作用.本文的数值试验显示时间二阶积分波场全波形反演方法应用于缺失4 Hz以下频谱成分数据和全频带数据得到的两个反演结果没有明显的差异.

致谢

感谢审稿专家和编辑部的支持.

参考文献
Bozdağ E, Trampert J, Tromp J. 2011. Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements. Geophysical Journal International , 185 (2) : 845-870. DOI:10.1111/gji.2011.185.issue-2
Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion. Geophysics , 60 (5) : 1457-1473. DOI:10.1190/1.1443880
Chauris H, Donno D, Calandra H. 2012. Velocity estimation with the normalized integration method.//74th EAGE Conference and Exhibition incorporating EUROPEC 2012. SPE, EAGE, W020.
Chew W C. Waves and Fields in Inhomogeneous Media. New York: Dover Publications, 1990 .
Chi B, Dong L, Liu Y. 2013. Full waveform inversion based on envelope objective function.//75th EAGE Conference and Exhibition incorporating SPE EUROPEC 2013. EAGE, TU-P04-09.
Chi B X, Dong L G, Liu Y Z. 2014. Full waveform inversion method using envelope objective function without low frequency data. Journal of Applied Geophysics , 109 : 36-46. DOI:10.1016/j.jappgeo.2014.07.010
Donno D, Chauris H, Calandra H. 2013. Estimating the background velocity model with the normalized integration method.//75thEAGE Conference and Exhibition incorporating SPE EUROPEC 2013. EAGE, Tu-07-04.
Lailly P. 1983. The seismic inverse problem as a sequence of before stack migrations.//Conference on Inverse Scattering, Theory and Applications. Philadelphia, PA:Society of Industrial and Applied Mathematics, 206-220.
PrattR G, ShinC, HicksG J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International , 133 (2) : 341–362.
Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part I:Theory and verification in a physical scale model. Geophysics , 64 (3) : 888-901. DOI:10.1190/1.1444597
Shin C, Min D J. 2006. Waveform inversion using a logarithmic wavefield. Geophysics , 71 (3) : R31-R42. DOI:10.1190/1.2194523
Shin C, Pyun S, Bednar J B. 2007. Comparison of waveform inversion, Part 1:Conventional wavefield vs. logarithmic wavefield. Geophysical Prospecting , 55 (4) : 449-464. DOI:10.1111/gpr.2007.55.issue-4
Shin C, Cha Y H. 2008. Waveform inversion in the Laplace domain. Geophysical Journal International , 173 (3) : 922-931. DOI:10.1111/gji.2008.173.issue-3
Shin C, Cha Y H. 2009. Waveform inversion in the Laplace-Fourier domain. Geophysical Journal International , 177 (3) : 1067-1079. DOI:10.1111/gji.2009.177.issue-3
Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging:A strategy for selecting temporal frequencies. Geophysics , 69 (1) : 231-248. DOI:10.1190/1.1649391
Sirgue L. 2006. The importance of low frequency and large offset in waveform inversion.//68th EAGE Conference & Technical Exhibition. SPE, EAGE, A037.
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics , 49 (8) : 1259-1266. DOI:10.1190/1.1441754
Tarantola A. 1986. A strategy for nonlinear elastic inversion of seismic reflection data. Geophysics , 51 (10) : 1893-1903. DOI:10.1190/1.1442046
Tarantola A. 1987. Inverse problem theory:Methods for data fitting and model parameter estimation. Amsterdam:Elsevier Science Publ. Co., Inc.
Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics , 74 (6) : WCC1-WCC26. DOI:10.1190/1.3238367
Vogel C R. 2002. Computational methods for inverse problems. Philadelphia, PA, US:Society for Industrial and Applied Mathematics.
Wu R S, Luo J R, Wu B Y. 2014. Seismic envelope inversion and modulation signal model. Geophysics , 79 (3) : WA13-WA24. DOI:10.1190/geo2013-0294.1