时间二阶积分波场的全波形反演

陈生昌， 陈国新

浙江大学地球科学学院， 杭州　310027



时间二阶积分波场的全波形反演

陈生昌， 陈国新

浙江大学地球科学学院， 杭州310027

通过对波场的时间二阶积分运算以增强地震数据中的低频成分，提出了一种可有效减小对初始速度模型依赖性的地震数据全波形反演方法—时间二阶积分波场的全波形反演方法.根据散射理论中的散射波场传播方程，推导出时间二阶积分散射波场的传播方程，再利用一阶Born近似对时间二阶积分散射波场传播方程进行线性化.在时间二阶积分散射波场传播方程的基础上，利用散射波场反演地下散射源分布，再利用波场模拟的方法构建地下入射波场，然后根据时间二阶积分散射波场线性传播方程中散射波场与入射波场、速度扰动间的线性关系，应用类似偏移成像的公式得到速度扰动的估计，以此建立时间二阶积分波场的全波形迭代反演方法.最后把时间二阶积分波场的全波形反演结果作为常规全波形反演的初始模型可有效地减小地震波场全波形反演对初始模型的依赖性.应用于Marmousi模型的全频带合成数据和缺失4 Hz以下频谱成分的缺低频合成数据验证所提出的全波形反演方法的正确性和有效性，数值试验显示缺失4 Hz以下频谱成分数据的反演结果与全频带数据的反演结果没有明显差异.

散射波方程; 时间二阶积分; 增强低频; Born近似; 全波形反演

1　引言

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

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

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

2　方法原理

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

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

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

(1-1)

(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)

(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)

(11)

2.2波场的时间二阶积分

(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中的时间二阶积分波场信号是受波场传播方程控制，是一种具有明确数学物理意义的地震波场.

2.3时间二阶积分波场的全波形反演

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

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

图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

(14)

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

(15)

(16)

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

(17)

或

(18)

(19)

(20)

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

(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)

(30)

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

(31)

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

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

(32)

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

(33)

(33-1)

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

(34)

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

图6　Marmousi模型的速度分布Fig.6　The velocity distribution of Marmousi model

3　数值试验

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

3.1全频带地震数据试验

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

图7　用于反演的常梯度初始速度模型Fig.7　The initial velocity model with constant gradient for inversion

图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作初始速度模型的包络反演加常规全波形反演的反演结果.

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

图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

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

图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

图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

图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以下频谱成分数据和全频带数据得到的两个反演结果没有明显的差异.

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

Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion.Geophysics, 60(5): 1457-1473.

Chauris H, Donno D, Calandra H. 2012. Velocity estimation with the normalized integration method.∥ 74thEAGE Conference and Exhibition incorporating EUROPEC 2012. SPE, EAGE, W020.

Chew W C. 1990. Waves and Fields in Inhomogeneous Media. New York: Dover Publications.

Chi B, Dong L, Liu Y. 2013. Full waveform inversion based on envelope objective function.∥ 75thEAGE 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.JournalofAppliedGeophysics, 109: 36-46.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. Pratt R G, Shin C, Hicks G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 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.

Shin C, Min D J. 2006. Waveform inversion using a logarithmic wavefield.Geophysics, 71(3): R31-R42.

Shin C, Pyun S, Bednar J B. 2007. Comparison of waveform inversion, Part 1: Conventional wavefield vs. logarithmic wavefield.GeophysicalProspecting, 55(4): 449-464.

Shin C, Cha Y H. 2008. Waveform inversion in the Laplace domain.GeophysicalJournalInternational, 173(3): 922-931.

Shin C, Cha Y H. 2009. Waveform inversion in the Laplace-Fourier domain.GeophysicalJournalInternational, 177(3): 1067-1079.

Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies.Geophysics, 69(1): 231-248.

Sirgue L. 2006. The importance of low frequency and large offset in waveform inversion.∥ 68thEAGE Conference & Technical Exhibition. SPE, EAGE, A037.

Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation.Geophysics, 49(8): 1259-1266.

Tarantola A. 1986. A strategy for nonlinear elastic inversion of seismic reflection data.Geophysics, 51(10): 1893-1903.

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.

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.

(本文编辑刘少华)

Full waveform inversion of the second-order time integral wavefield

CHEN Sheng-Chang， CHEN Guo-Xin

SchoolofEarthSciences，ZhejiangUniversity，Hangzhou310027,China

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.

Propagation equation of scattering wave; Second-order time integral; Enhancement low-frequency components; Born approximation; Full waveform inversion

10.6038/cjg20161021.

国家自然科学基金项目(41074133，41374001)资助．

陈生昌，男,1965年生，教授，博士生导师．主要从事勘探地球物理和计算地球物理研究. E-mail: chenshengc@zju.edu.cn

10.6038/cjg20161021

P631

2015-11-10, 2016-03-18收修定稿

陈生昌， 陈国新. 2016. 时间二阶积分波场的全波形反演. 地球物理学报，59(10):3765-3776,

Chen S C， Chen G X. 2016. Full waveform inversion of the second-order time integral wavefield.ChineseJ.Geophys. (in Chinese),59(10):3765-3776,doi:10.6038/cjg20161021.