基于包络目标函数的弹性波波形反演

王官超,杜启振

(中国石油大学(华东)地球科学与技术学院,山东青岛266580)



基于包络目标函数的弹性波波形反演

王官超,杜启振

(中国石油大学(华东)地球科学与技术学院,山东青岛266580)

摘要:弹性波场的复杂性致使反演问题非线性性增强,寻优过程极易陷入局部极值。多尺度反演能在一定程度上提高反演的稳定性及分辨率,但多尺度反演策略无法解决实际数据的低频缺失问题。针对该问题,提出了一种新的基于包络目标函数的弹性波多尺度反演策略。首先将包络反演方法应用到弹性波,恢复地下弹性介质参数中的低波数信息；然后将包络反演结果作为常规多尺度反演的输入,进一步提高反演的稳定性。Marmousi2模型试算结果表明该策略能够更加准确地反演纵、横波速度以及密度信息。

关键词：包络反演；低频缺失；多尺度；弹性波全波形反演

全波形反演技术(Full Waveform Inversion,FWI)是一种从全波场地震记录中定量提取地下介质参数的数据匹配技术,该技术利用波场的走时、振幅以及相位等多种信息,通过不同优化方法获得地下弹性介质参数,其分辨率明显高于传统基于射线理论的反演方法,反演结果不仅能为叠前深度偏移提供更为可靠的偏移速度,而且能直接提供精度较高的地下介质弹性参数信息。但是,该技术存在计算量大、非线性性强、初始速度场构建困难以及波场之间的不确定性等缺点,导致波形反演容易陷入循环跳动和局部极小值等问题,当涉及到弹性多参数反演时这些问题表现得尤为严重,其中初始模型的准确与否直接影响反演结果的正确性。

通常情况下,背景宏观速度模型可以利用初至走时层析[1]、立体层析[2-3]、反射/折射层析、叠加速度以及偏移速度分析等方法获得。当初始模型远离真实模型时,全波形反演目标函数存在多极值,往往导致最优化方法收敛到局部极小[4-5]。多尺度反演策略[6-8]能在一定程度上解决该问题。但是,由于采集技术的制约，实际地震勘探中很难记录到5Hz以下的频率成分,准确的含低波数信息的初始速度模型的构建仍然是FWI中亟需解决的关键问题。针对初始模型构建困难这一瓶颈,许多学者进行了大量研究,SHIN等[9-10]提出了利用Laplace域阻尼波场零频分量反演长波长宏观速度模型作为频率域全波形反演的初始速度模型。事实上,将Laplace域和Fourier域的波形反演集合成一个过程,既可以提高计算效率,同时也可以提高算法稳定性以及降低对初始速度模型的依赖性。XU[11]在反演过程中引入反偏移技术，进行了基于反射波的波形反演研究,并在实际资料应用中取得较好效果。胡光辉等[12]采用基于早至波的特征波波形反演建模方法解决收敛过程中由于初始模型不准确导致的局部极小问题,使基于特征波的全波形反演逐渐受到人们的关注。地震信号经过包络变换后的能量主要集中在极低频区域,这种含有极低频信息的包络信号可以有效反映地下介质的大尺度信息。因此,WU等[13-14]在声波假设下提出基于包络目标函数的全波形反演方法建立低频速度模型。通过对低频信息进行补偿、频率外推等方式,减少低频信息缺失带来的不利影响,同样是一种可行策略。

基于前人的研究成果,本文提出了一种基于弹性波包络反演方法的多尺度反演策略。首先将声波包络反演方法推广应用到弹性波,恢复地下弹性介质参数中的大尺度信息,降低弹性多参数反演的非线性性；其次将弹性波包络反演获得的含有低频信息的反演结果作为常规多尺度反演的初始模型,进一步解决反演过程中存在的不稳定性问题。Marmousi2模型试算结果表明,结合弹性波包络反演和多尺度反演的反演策略能够更加准确地反演纵、横波速度以及密度信息,实现弹性波多参数反演。

1基于包络目标函数的低频反演方法

常规全波形反演通过匹配观测数据与模拟数据,使其数据残差δu(t)=u(t)-d(t)达到最小,定量提取地下介质弹性参数m,目标函数定义为:

(1)

在弹性介质中,多分量地震数据可通过速度应力方程进行数值模拟合成,即:

式中:ρ(x)是密度；u(x,t)和T(x,t)分别为速度分量和应力分量；F(x,t)为震源函数；C(x)为弹性常数。利用高阶交错网格有限差分方法进行波动方程求解[15],采用PML吸收边界条件。

包络反演与常规时间域反演的区别在于其匹配的是地震数据的包络而非地震数据。根据声波包络反演的目标函数[12]可以直接给出弹性波包络反演目标函数:

(3)

式中:d(t)=[dx(t),dz(t)]和u(t)=[ux(t),uz(t)]分别为观测和合成的多分量地震记录；eobs(t)和esyn(t)分别为多分量地震记录d(t)和u(t)的包络；dh(t)和uh(t)为相应波形的希尔伯特变换；E为包络残差；参数p控制包络信号的能量,可以是任意正数,p值越大包络信号中深层反射的能量损失越严重,导致包络反演无法刻画深层结构,在本文的数值测试中p=2。

通过推导弹性波包络目标函数((3)式)相对于模型参数m的一阶偏导可得:

(4)

(4)式可表示为雅克比矩阵J以及有效残差向量η的乘积:

(5)

(6a)

(6b)

(6c)

式中:v=(vx,vz),w=(wx,wz,wxx,wzz,wzx)分别代表正传波场和反传残差波场。因为不同参数之间存在着相互关系,目标函数对于其它相关参数m=[vP,vS,ρ]的梯度可以通过链式法则进行求解。另外,迭代过程采用抛物线插值方法估计更新步长。

包络反演假设地震波形为调幅信号,包络算子是一种非线性的解调算子,可反映隐藏在调幅信号中的低频调制信号的波形变化。图1a和图1b分别展示了两个余弦衰减信号叠加的调幅信号和单道地震记录的波形及其包络,图1c和图1d为对应的频谱特征。从频谱上可看出,包络变换后信号的能量主要集中在低频区域,横坐标越接近0，振幅谱能量就越强(图1c和图1d),包络反演正是利用这些极低频信息来恢复地下介质参数的长波长分量。

图1　调幅信号及其包络(a),单道地震记录及其包络(b),对应图1a的频谱特征(c)以及对应图1b的频谱特征(d)

2基于包络目标函数的多尺度反演策略

BUNKS等[6]采用Hamming窗低通滤波器对地震子波和数据进行时间尺度分解,首次实现了时间域多尺度全波形反演策略。BOONYASIRIWAT等[19]指出利用Hamming窗进行滤波会有高频泄漏。为解决Hamming窗滤波的高频泄漏问题,可以选用其它无泄漏或少泄漏的滤波窗函数。本文中采用维纳滤波器进行低通滤波,实现多尺度反演策略,维纳滤波器表达式为:

(7)

式中:fWiener是维纳滤波器；WOriginal是原始信号；WTarget为期望得到的信号；ω是角频率；ε是为了避免不稳定现象发生而设定的一个很小的值。利用维纳滤波器将波场数据分解成不同频带,由低频至高频带进行逐频带反演，从而实现多尺度反演策略,降低反演过程陷入局部极小值的风险。图2a为利用20Hz Ricker子波生成的单炮Z分量记录,图2b 为利用5Hz Ricker子波对图2a的记录进行维纳滤波后的地震记录。由于实际地震数据往往缺失低频信息,基于低通滤波的常规多尺度反演策略无法反演模型参数中的大尺度信息,导致反演陷入局部极值或循环跳动。

图2　利用20Hz Ricker子波生成的单炮记录(a)以及利用5Hz Ricker子波对图2a进行维纳滤波后的单炮记录(b)

由于包络算子具有解调地震数据中隐含低频信息的作用，因此，利用包络算子解调出的低频信息可以有效恢复地下介质的长波长成分,进而解决常规多尺度反演无法克服的低频缺失或不足问题。同时,包络变换后的信号在富含低频信息的同时也存在相应的高频信息。图3是主频为7Hz(实线和虚线)、20Hz(点线和点划线)Ricker子波及其包络的频谱响应,可以看出包络信号携带的高频信息与原始信号主频带的高低是相对的。如果原信号主频较高，变换后的信号频率也相对较高，反之亦然。

图3　20Hz Ricker子波及其包络的频谱响应以及7Hz Ricker子波及其包络的频谱响应

我们通过数值算例来说明包络中高频信息对包络反演结果的影响。图4是以线性模型作为初始模型分别利用图3中两种震源子波进行弹性包络反演获得的纵波速度,可以看出,当高频信息较多时，包络反演结果不稳定(图4b黑色矩形框)。为提高包络反演的稳定性,在进行包络反演之前利用类似于时间域多尺度反演对原始地震数据进行低通滤波处理,这样就把弹性包络反演与常规多尺度反演有机结合在一起,自适应的实现了一种从极低频开始的多尺度反演策略。

因此,基于包络目标函数的多尺度反演策略的具体实现过程为:首先,对地震数据进行时间尺度分解；其次,采用低通滤波后的地震数据实施包络反演,恢复地下弹性介质参数中的大尺度信息；最后,将弹性波包络反演结果作为初始模型进行常规多尺度反演以获取精确的最终反演结果。

图4　弹性波包络反演得到的纵波速度a 7Hz震源子波; b 20Hz震源子波

3Marmousi2模型测试

在反演测试过程中采用两种不同类型的震源子波合成多分量地震数据(主频为20Hz的全频带Ricker子波和低频缺失(<7Hz)Ricker子波)来验证弹性波包络反演恢复地下模型参数大尺度信息的能力及其在地震数据缺失低频情况下的适用性。为保证包络反演的稳定性,首先利用7Hz Ricker子波对两种类型震源子波生成的地震记录进行维纳低通滤波,然后进行包络反演。图6和图7 为以线性模型作为初始模型进行包络反演迭代30次后得到的平滑背景场。可以看到,针对两种不同震源子波包络反演获得的背景模型构造非常相近,这说明了包络反演对低频信息的依赖性很弱,即使在低频缺失的情况下同样可以反演得到模型中的低频信息,为常规多尺度波形反演提供含有低频信息的初始速度模型,降低全波形反演的非线性性。

图5　真实Marmousi2纵波速度模型(a)及其线性模型(b)

图6　线性模型作为初始模型的弹性波包络反演30次迭代后的结果(全频带震源子波)a 纵波速度; b 横波速度; c 密度

图7　线性模型作为初始模型的弹性波包络反演30次迭代后的结果(低频缺失震源子波)a 纵波速度; b 横波速度; c 密度

图8和图9分别为以图6和图7中包络反演结果作为初始模型,选取5,10,20Hz的Ricker子波对地震记录进行维纳滤波进而实施多尺度反演得到的最终结果。相比于直接利用线性模型作为初始模型进行多尺度反演的结果(图10和图11),利用包络反演加多尺度反演得到的最终结果在构造(黑色椭圆处)上得到了明显提高,构造假象明显减少。对比图10和图11可以看出,利用低频缺失震源子波的反演结果相比于全频带震源子波反演结果的构造假象更为严重,这说明传统全波形反演对震源的频谱特性更为敏感。

图8　以图6所示包络反演结果作为初始模型得到的多尺度反演结果(全频带震源子波)a 纵波速度; b 横波速度; c 密度

图9　以图7所示包络反演结果作为初始模型得到的多尺度反演结果(低频缺失震源子波)a 纵波速度; b 横波速度; c 密度

图10　以线性模型作为初始模型得到的多尺度反演结果(全频带震源子波)a 纵波速度; b 横波速度; c 密度

图11　以线性模型作为初始模型得到的多尺度反演结果(低频缺失震源子波)a 纵波速度; b 横波速度; c 密度

4结论

通过理论研究与模型试算得到以下认识:①包络算子可以解调出地震数据中隐含的低频信息。包络反演利用解调出的低频信息能有效恢复模型中的大尺度信息,在地震数据缺失低频的情况下同样适用。②包络信号在富含低频信息的同时也存在相应的高频信息。高频信息较多时会导致包络反演不稳定,通过将包络反演与基于低通滤波的多尺度反演相结合,自适应的实现了一种从极低频开始的多尺度反演策略,最大限度地降低了反演对初始模型的依赖性。

实际应用是全波形反演的挑战,野外数据低频缺失或不足是影响波形反演获得全波数带参数估计结果的主要因素。包络反演技术可恢复隐藏在地震资料中的低频信息,建立低频初始模型,具有很好的发展前景。

参考文献

[1]NOLET G.Seismic tomography:with applications in global seismology and exploration geophysics [M].Bonston:Kluwer Academic Publishers,1987,1-24

[2]WOODWARD M J.Wave-equation tomography[J].Geophysics,1992,57(1):15-26

[3]PRIEUX V,OPERTO S,BROSSIER R,et al.Application of acoustic full waveform inversion to the synthetic Valhall model[J].Expanded Abstracts of 77thAnnual Internat SEG,Mtg,2007:2268-2272

[4]VIRIEUX J,OPERTO S.An overview of full-waveform inversion in exploration geophysics[J].Geophysics,2009,74(6):WCC1-WCC26

[5]杨勤勇,胡光辉,王立歆.全波形反演研究现状及发展趋势[J].石油物探,2014,53(1):77-83

YANG Q Y,HU G H,WANG L X.Research status and development trend of full waveform inversion [J].Geophysical Prospecting for Petroleum,2014,53 (1):77-83

[6]BUNKS C,SALECK F M,ZALESKI S,et al.Multiscale seismic waveform inversion[J].Geophysics,1995,60(5):1457-1473

[7]SIRGUE L,PRATT R G.Efficient waveform inversion and imaging:a strategy for selecting temporal frequencies[J].Geophysics,2004,69(1):231-248

[8]李媛媛,李振春,张凯.频率域多尺度弹性波全波形反演[J].石油物探,2015,54(3):317-323

LI Y Y,LI Z C,ZHANG K.Multi-scale elastic full waveform inversion in frequency domain [J].Geophysical Prospecting for Petroleum,2015,54 (3):317-323

[9]SHIN C,HA W.A comparison between the behavior of objective functions for waveform inversion in the frequency and Laplace domains[J].Geophysics,2008,73(5):VE119-VE133

[10]SHIN C,CHA Y H.Waveform inversion in the Laplace domain [J].Geophysical Journal International,2008,173(3):922-931

[11]XU S,WANG D,CHEN F,et al.Full waveform inversion for reflected seismic data[J].Expanded Abstracts of 74thEAGE Annual Conference,2012:114-117

[12]胡光辉,王立歆,王杰,等.基于早至波的特征波波形反演建模方法[J].石油物探,2015,54 (1):71-76

HU G H,WANG L X,WANG J,et al.Characteristics waveform inversion based on early arrival waves [J].Geophysical Prospecting for Petroleum,2015,54 (1):71-76

[13]WU R S,LUO J R,WU B.Seismic envelope inversion and modulation signal model[J].Geophysics,2014,79(3):WA13-WA24

[14]LUO J R,WU R S.Seismic envelope inversion:reductionoflocal minima and noise resistance [J].Geophysical Prospecting,2015,63(3):597-614

[15]董良国,马在田,曹景忠.一阶弹性波交错网格高阶有限差分解法[J].地球物理学报,2000,43(3):411-419

DONG L G,MA Z T,CAO J Z.Higher order difference method for first order elastic wave equation with staggered grid [J].Chinese Journal of Geophysics,2000,43(3):411-419

[16]TARANTOLA A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[17]PLESSIX R E.A review of the adjoint-state method for computing the gradient of a functional with geophysical applications [J].Geophysical Journal International,2006,167(2):495-503

[18]WANG J,ZHOU H,TIAN Y,et al.A new scheme for elastic full waveform inversion based on velocity-stress wave equations in time domain[J].Expanded Abstracts of 82ndAnnual Internat SEG,Mtg,2012:1-5

[19]BOONYASIRIWAT C,VALASEK P,ROUTH P,et al.An efficient multiscale method for time-domain waveform tomography[J].Geophysics,2009,74(6):WCC59-WCC68

(编辑：朱文杰)

Elastic full waveform inversion based on envelope objective function

WANG Guanchao,DU Qizhen

(SchoolofGeosciences,ChinaUniversityofPetroleum,Qingdao266580,China)

Abstract:The complexity of elastic wavefield will increase the nonlinearity of inversion and make the inverse problem fall into local minima.Multi-scale inversion algorithms can improve the stability and resolution of inversion to a certain extent,but it cannot solve the problem of low frequency missing in real data.To alleviate this problem,we propose a new multi-scale scheme based on envelope inversion to obtain the elastic parameters accurately.The main procedures are as follows:Firstly,we extend the envelope inversion method to elastic medium to estimate the low wavenumber components of elastic parameters; then,the envelope inversion results are used as the initial model of conventional multi-scale inversion to further improve the stability of the inversion.The Marmousi2 model test results show that the new multi-scale scheme can estimate P- and S- wave velocities more accurately,as well as the density information.

Keywords:envelope inversion,low frequency missing,multi-scale,elastic full waveform inversion

文章编号：1000-1441(2016)01-0133-09

DOI:10.3969/j.issn.1000-1441.2016.01.017

中图分类号：P631

文献标识码：A

基金项目：国家自然科学基金(41174100和41074087)和中国石化地球物理重点实验室开放基金联合资助。

作者简介：王官超(1990—),男,硕士在读,主要从事全波形反演研究工作。通讯作者：杜启振(1969—),男,教授,博士生导师,主要从事地震弹性波传播理论与成像等研究工作。

收稿日期：2015-07-07;改回日期：2015-08-24。

王官超,杜启振.基于包络目标函数的弹性波波形反演[J].石油物探,2016,55(1):-141

WANG Guanchao,DU Qizhen.Elastic full waveform inversion based on envelope objective function[J].Geophysical Prospecting for Petroleum,2016,55(1):-141

This research is financially supported by the National Science Foundation of China (Grant Nos.41174100,41074087) and the SINOPEC Key Laboratory of Geophysics Project.