地球物理学报  2015, Vol. 58 Issue (6): 1998-2010   PDF    
不依赖子波、基于包络的FWI初始模型建立方法研究
敖瑞德, 董良国, 迟本鑫    
同济大学海洋地质国家重点实验室, 上海 200092
摘要:地震全波形反演(FWI)从理论走向实际面临着诸多难题, 其中之一就是需要一个较高精度的初始模型, 另一个难题就是需要一个较为精确的震源子波, 初始模型和震源子波的准确程度严重影响着全波形反演的最终结果.为此, 本文提出了不依赖子波、基于包络的FWI初始模型建立的方法, 建立了相应的目标函数, 推导出了反演的梯度, 给出了伴随震源的表达式, 理论上分析了不依赖子波FWI的可行性.在数值试验中, 讨论了参考道的选取方式, 通过分析归一化目标函数收敛速率, 认为近偏移距参考道优于远偏移距参考道, 在地震数据含干扰噪音时, 平均道作为参考道要优于最小偏移距参考道.通过包络、包络对数、包络平方三种目标函数反演结果的比较, 发现包络对数目标函数对深层的反演效果最好.通过不同子波的试验进一步验证了本方法的正确性.
关键词全波形反演     子波     包络     初始模型     参考道     目标函数    
Source-independent envelope-based FWI to build an initial model
AO Rui-De, DONG Liang-Guo, CHI Ben-Xin    
State Key Laboratory of Marine Geology, Tongji University, Shanghai 200092, China
Abstract: Full waveform inversion (FWI) is facing many difficulties in practice. One of the issues is how to construct a relatively precise initial model. Another one is the treatment of source wavelet. Both seriously affect the final inversion result. To tackle this problem, we propose a source-independent envelope-based FWI approach to build an initial model for the traditional FWI.On one hand, the cross-convolved wavefields, including the convolution of observed wavefields with a reference trace from the modeled wavefields and the convolution of modeled wavefields with a reference trace from the observed wavefields, can eliminate source wavelet effect in the misfit function of waveform inversion.
On the other hand, the envelope of wavefields contains abundant low-frequency components, which are crucial for initial model inversion. Therefore, we use the envelope of cross-convolved wavefields to establish misfit functions, the feasibility of which is proved by theoretical analysis. The expressions of the corresponding gradient and adjoint source are derived. Furthermore, we use a Marmousi synthetic data example to test our method.
The proposed method is tested by using different reference traces, misfit functions and source wavelets. (1) The normalized misfit of long offset reference trace converges more slowly than that of short-offset reference trace for noise-free data and the average trace as reference trace is better than the minimum-offset reference trace for noisy data. (2) The logarithmic envelope is the best and the squared envelope is the worst misfit function to invert the deep part of the model. The initial model built by our method choosing the logarithmic misfit function and minimum-offset reference trace is much better than that of the traditional envelope-based FWI when the source wavelet is incorrect. (3) Using other kinds of source wavelets, such as 20 Hz Ricker wavelet without low frequency below 5 Hz and 12 Hz Gauss wavelet, both of which are quite different from the true 8 Hz Ricker wavelet, the inversion result of our method can provide a good initial model for the traditional FWI. In addition, the initial model can be used in seismic depth migration, and the imaging result is close to the RTM result with true velocity model and much better than that with gradient initial model.
In conclusion, first, we propose a source-independent envelope-based FWI approach to build an initial model for the traditional FWI without any prior information when the source wavelet is unknown and even the data lack low frequency. The initial model built by our method can also be used as the velocity model for pre-stack depth migration. Second, the short-offset reference trace is better than long-offset reference trace for noise-free data and average trace as reference trace is better than the minimum-offset reference trace for noisy data. Third, the logarithmic envelope misfit function is the optimal misfit function. Finally, by changing the source wavelet, it is further proved that our method is source-independent.
Key words: Full waveform inversion     Source wavelet     Envelope     Initial model     Reference trace     Misfit function    
1 引言

全波形反演(Full waveform inversion,简称FWI)作为一种高分辨率地震反演成像方法,是近些年来勘探地球物理界研究的热点.20世纪80年代,Tarantola(1984)从波动方程出发,利用模拟波场与观测数据残差反传波场的互相关构造梯度方向,建立了一套时间域全波形反演的理论框架.20世纪90年代,Pratt等(1998)Pratt(1999)又将其发展到频率域,利用炮点和检波点的互易性原理推导出梯度的高效算法.这些研究是当前地震全波形反演的理论基础.但是,FWI目前在实际中还没得到很好的应用,其中一个问题就是如何得到一个较为精准的初始速度模型,另一个问题就是如何处理震源子波.

FWI普遍采用局部优化方法,但由于初始模型不够准确、低频成分的缺失、噪音的存在以及复杂波现象的近似描述等原因,局部优化方法极易导致目标函数收敛到局部极值(Virieux and Operto,2009).提供较准确的初始速度模型可以一定程度上缓解这一问题.要构造较准确的初始速度模型,观测地震数据中的低频成分非常关键.如果原始数据中含有丰富的低频信息,则可以在FWI过程中通过一些反演策略更好地反演速度的背景成分,如多尺度(Bunks et al.,1995;Sirgue and Pratt,2004)、复频率(Brossier et al.,2009)等反演策略.但由于接收数据的有限带宽和必要的去噪处理等原因,在实际地震数据中往往缺少可靠的低频成分.因此,如何获得地震数据中的低频信息,成为用FWI构建初始速度模型的关键.地震道的包络主要体现了地下模型的大尺度信息,用其反演地下速度变化的长波长分量很合适.基于包络目标函数的报道最早出现在2011年,但他们只讨论了相应的有限频核函数问题(Bozda et al.,2011),没有给出利用包络构造目标函数的FWI反演结果,也没有解释其可操作性.在2013年,基于声波包络的FWI的可行性得到了证实,用该方法获得的结果作为初始模型,再进行常规FWI反演的效果非常理想(Chi et al.,20132014;Wu et al.,2014).但在上述文献中,这些方法都建立在子波已知的基础上.

FWI要通过地震波正演模拟来实现,需要知道震源子波.但在实际地震数据中,子波往往是未知的,这成为全波形反演的又一大挑战.最直接的解决方法就是从地震数据中提取直达波作为震源子波,但直达波并不是真实的子波.Song等(1995)则将子 波的计算统一在速度模型的反演中,通过迭代的方式不断更新震源子波和速度模型,但如果速度模型与真实模型有较大的偏差,计算的子波也会不准确,最终造成反演的失败.另一种思路就是在反演过程中消除子 波的影响.Lee和Kim(2003)Zhou和Greenhalgh(2003)等提出了基于声波方程的不依赖子波的FWI,这种方法通过在频率域将地震数据与参考道波场相除,从而获得归一化的波场来消除震源子波的影响.该方法的缺点在于需要直接计算Jacobian矩阵,计算量大.Choi等(2005)在此基础上提出了利用地震数据的振幅谱消除子波的两种方法,即反褶积法与褶积法:将观测数据与模拟数据的振幅与各自参考道的振幅先相除后作差,以此构建目标函数的方法称为反褶积法;将观测数据与模拟数据的振幅与相互的参考道的振幅先相乘后作差,以此构建目标函数的方法称为褶积法.两种方法的梯度均可通过伴随状态法求取,减少了计算量,但这些方法 都是在频率域实现的.之后,Choi和Alkhalifah(2011)又提出了时间域不依赖子波的FWI,他们在时间域将模拟数据和观测数据互相与参考道波场褶积消除子波的影响,同时验证了该方法在加噪数据上的可行性,但是该方法仍需要提供较准确的初始模型.

本文在前人工作基础上,提出不依赖子波、基于包络的FWI初始速度模型的建立方法,该方法克服了FWI对子波的依赖性,降低了FWI的非线性程度,在观测地震数据缺少低频成分情况下也可以进行有效的FWI反演,从而为下一步FWI提供较为准确的初始模型.

2 方法原理 2.1 目标函数

全波形反演的目标函数一般定义为模拟数据与观测数据残差的 L2 泛函,形式如下:

其中 s和r 分别是炮号和道号,ui,j和di,j 分别是第 i 炮的第 j 个检波器处的模拟数据和观测数据.将 ui,j和di,j 写成格林函数和子波的褶积形式,目标函数变为

其中,g 是格林函数,w 是子波,*代表褶积,上标u和d分别代表模拟数据和观测数据.

为了避免在全波形反演中估算子波,可以通过交叉褶积波场的方法消除子波的影响(Choi and Alkhalifah,2011).对于每一个炮集,模拟数据与同一炮观测数据的第 k 个参考道进行褶积,获得模拟褶积波场;观测数据与同一炮模拟数据的第 k 个参考道进行褶积,获得观测褶积波场.此时,新的目标函数定义为

其中 ui,k和di,k 分别是第 i 炮的第 k 个参考道的模拟数据和观测数据. ui,j*di,k 为模拟褶积波场,di,j*ui,k 为观测褶积波场.

将(3)式中的模拟数据与观测数据都表示为格林函数与震源子波的褶积形式:

从(4)式可以发现,其中的减号左右两边都与共同的新震源子波(wu*wd)进行了褶积,这样在目标函数中就消除了实际子波与模拟子波差异的影响,而由(2)式可知,目标函数(1)却与子波差异有关.因此,基于目标函数(3)的全波形反演结果不再依赖于准确的子波(图 2b),而目标函数(1)依赖于子波,采用不准确的子波将严重影响反演结果的精度(图 2a).

进一步地,利用褶积波场的包络提取低频信息:

其中,Ecobs 就是观测褶积波场的包络,Ec 就是模拟褶积波场的包络,H(f)为实信号 f 的希尔伯特变换.将波场写成格林函数与震源子波的褶积形式后,EcobsEc 中都拥有共同的子波项(wu*wd). 所以此时建立的模拟褶积波场包络和观测褶积波场包络的泛函仍然不依赖子波差异:

除了(7)式中包络误差形式的目标函数,常用的目标函数还有下面的包络平方目标函数(Wu et al.,2014)和包络对数目标函数(Chi et al.,2014),它们都不受子波的影响:

2.2 梯度和伴随震源

在目前基于局部优化的全波形反演过程中,目标函数梯度的求取是关键.(10)式是包络对数目标函数(9)对模型参数 p 的梯度,本文中 p 代表纵波速度(详细的推导过程请参考附录A),

(10)式中,h 为希尔伯特算子.可以表示为上行波的格林函数 gx,ju和虚拟震源 vi,xu 的褶积,其中 x 表示地下模型空间的变化.对于常密度声波方程,根据Born近似,vi,xu 可以表示为或者 2pΔui,x. 虚拟震源的表达式与波动方程的类型有关.虚拟震源相似于下行波格林函数与震源子波的褶积(或者加权的正传震源波场)(Choi and Alkhalifah,2011).因此,上式仍为Born近似类型的反演.将改为 gx,juvi,xu 的褶积形式(同理),则(10)式变为

其中代表互相关运算.将伴随震源设为

格林函数 gx,ju 满足互易性准则,因此 gx,ju=gj,xu. 所以在(11)式中,互相关是任意 检波点 j 处伴随震源 r32 的反传波场,而互相关则是参考道 k 处伴随震源 r31 的反传波场.

当取参考道 di,k=ui,k=δ(t)时,褶积波场即为原始波场,参考道与空间位置无关,此时伴随震源只剩 r32. 应用 δ 函数的特殊性质,伴随震源进一步退化为

此时 EcobsEc等价于观测数据和模拟数据的包络.上式与依赖子波、基于包络的FWI对数目标函数的伴随震源表达式形式完全一致(Bozda et al.,2011),这说明不依赖子波、基于包络的FWI的伴随震源与依赖子波、基于包络的FWI的伴随震源拥有统一的形式.通过将参考道选为 δ(t),可以验证该方法的正确性.

当采用其他两种目标函数时,梯度的推导过程与上面类似.

对包络目标函数(7),在参考道 k 处的伴随震源为

在任意检波点 j 处的伴随震源为

对包络平方目标函数(8),在参考道 k 处的伴随震源为

在任意检波点 j 处的伴随震源为:

2.3 初始模型建立过程

利用本文提出的方法进行初始模型建立的步骤与时间域常规FWI反演的步骤基本一致,主要不同在于目标函数和伴随震源的计算,已在2.1和2.2节中分别给出.本文中的具体过程如下:(1)给定一个模型,正演计算每一炮的模拟波场;(2)选定合适的参考道,求取模拟褶积波场和观测褶积波场;(3)计算目标函数和伴随震源;(4)利用伴随震源计算伴随波场,并求取单炮梯度;(5)累加各单炮的梯度,从而获得全局梯度;(6)通过L-BFGS算法对梯度进行校正并求取步长,进行模型更新,直到满足终止条件,最终输出反演结果.

3 数值试验

采用的Marmousi速度模型(图 1a)参数如下:水平网格数nx=501,垂直网格数nz=140,水平网格大小dx=15 m,垂直网格大小dz=15 m,水层深30 m,水平总长度7500 m,垂直总长度2085 m.震源和检波器都位于地表处,共有80炮数据,第一炮位于195 m处,炮间距为90 m.第一个检波器位于0 m处,检波器间隔15 m.地震记录的采样间隔为0.0015 s,采样点数为2000,总记录时间为3 s.用高阶有限差分法(董良国等,2000)进行二维声波的正演.理论子波是主频为8 Hz的Ricker子波,延时0.15 s,子波最大振幅为1.0.

图 1 Marmousi准确模型(a)以及高斯平滑结果(b) Fig. 1 The Marmousi velocity model (a) and Gauss-window smoothed model (b)

即使采用相对精确的Gauss平滑初始模型(图 1b),利用错误的子波(主频3 Hz的Ricker子波,延时0.5 s,最大振幅2.0),常规FWI也会得到很差的反演结果(图 2a).而不依赖子波的FWI的反演结果对模型的高频成份进行了很好的恢复(图 2b).

图 2 子波错误时不同目标函数的FWI反演结果
(a) 目标函数(1)的常规FWI反演结果;(b) 目标函数(3)的不依赖子波的FWI反演结果.
Fig. 2 FWI results of different misfit functions with an incorrect wavelet
(a) The traditional FWI result corresponding to misfit function (1); (b) Source-independent FWI result corresponding to misfit function (3).

如果以梯度模型作为FWI的初始模型,常规FWI将不仅面对错误子波的问题,较差的初始模型还会使反演极易陷入局部极值,常规FWI的结果将变得更加糟糕.本文提出的不依赖子波、基于包络的FWI方法,就是要在不知道震源子波以及没有任何模型先验信息的条件下,为下一步的常规FWI反演提供一个含有可靠的低波数成分的初始模型.其中参考道的选择、方法的抗噪性、不同目标函数的比较以及不同子波的对比,是本文提出的方法中涉及的几个关键问题,下面通过数值试验和理论分析对这几个问题进行探讨.为了表达方便,将依赖于正确子波、基于包络的FWI称为传统的基于包络的FWI,在下面的试验中,它与本文提出的不依赖子波、基于包络的FWI均是在梯度模型(图 3)的基础上进行的.

图 3 梯度初始模型 Fig. 3 The gradient initial model
3.1 参考道的选择

本文提出的不依赖子波、基于包络的FWI方法中的一个关键问题,就是参考道应该如何选取.理论上,参考道的选择并不重要,然而实际中,参考道的选择对于反演结果有较大影响.关于参考道的选取,前人已经做了比较详细的分析.刘宝龙等(2014)比较了每个炮集中不同参考道的全波形反演结果,认为最小偏移距参考道为最优选择方案.但这种选取方案的FWI并没有涉及包络目标函数的问题.

为了比较不同参考道对本方法反演结果的影响,我们进行了两组试验,分别利用不含噪数据和含有噪声的数据进行测试和对比.选取的子波都为主频3 Hz的Ricker子波,延时0.5 s,最大振幅2.0,显然与观测数据的子波(8 Hz的Ricker子波,延时0.15 s,子波最大振幅为1.0)差别很大.

首先比较数据不含噪声时不同偏移距的参考道.理论上,最小偏移距参考道的能量最大;考虑实际波场的传播,子波具有时变和空变特征,偏移距越小,畸变越小,因此当采用不依赖子波、基于包络的FWI建立初始模型时,最小偏移距参考道应优于远偏移距参考道.但直接比较不同偏移距的反演结果难以看出差别(图 4),而从图 5的归一化的目标函数随迭代次数的变化过程可以看出,随着参考道的偏移距增加,收敛速度减慢.这个试验说明:近偏移距参考道优于远偏移距参考道.因此在下面的试验中选取零偏移距参考道.

图 4 无噪时选取不同偏移距参考道迭代16次的反演结果
(a) 偏移距750 m;(b) 偏移距600 m;(b) 偏移距0 m.
Fig. 4 The inversion results of different offset reference traces after 16 iterations for noise-free data
(a) 750 m-offset;(b) 600 m-offset;(c) 0 m-offset.

图 5 无噪时不同偏移距参考道反演过程中目标函数收敛速度对比
红色:偏移距750 m参考道;绿色:偏移距600 m参考道;蓝色:0 m偏移距参考道.
Fig. 5 The history of normalized misfit functions of 750 m-offset reference trace(red), 600 m-offset reference trace (green) and 0 m-offset reference trace(blue) for noise-free data.

对于含有噪声的数据,如果只选取某一道作为参考道,噪声的影响会很大.而取平均道作参考道后,不同地震道的噪声可以相互抵消,从而可以更好地保留信号中的有效信息.所以在第二组试验中,对原始地震记录加入信噪比10 dB的噪声(图 6),其中的噪声是产生的随机噪音与正演模拟的子波褶积的结果,以保证噪音与有效信号频带一致.选取每个炮集炮点附近100道的平均道作为参考道.对比在无噪数据和加噪数据情况下,零偏移距参考道和平均道作参考道时不依赖子波、基于包络的FWI初始模型建立结果(图 7),可以发现,本方法在初始模型建立过程中对噪声并不敏感,均很好地反演出了模型的背景成份.从归一化的目标函数(图 8)中也可以看出,当数据中含有噪声后,平均道作参考道的收敛速度要优于零偏移距参考道的收敛速度.

图 6 加入噪声的一炮地震记录(信噪比10 dB) Fig. 6 One shot record contaminated by noise with 10 dB S/N ratio

图 7 不同参考道迭代16次反演结果
(a) 无噪时零偏移距参考道反演结果;(b) 加噪时零偏移距参考道反演结果;
(c) 无噪时平均道作参考道反演结果;(d) 加噪时平均道作参考道反演结果.
Fig. 7 The inversion results with different reference traces after 16 iterations
(a) 0 m-offset reference trace for noise-free data; (b) 0 m-offset reference trace for noisy data; (c) Average trace
as reference trace for noise-free data; (d) Average trace as reference trace for noisy data.

图 8 不同参考道反演时,归一化目标函数收敛曲线对比
红色:平均道作参考道;蓝色:零偏移距参考道. (a) 无噪时归一化目标函数收敛情况;(b)加噪时归一化目标函数收敛情况.
Fig. 8 The history of normalized misfit functions for different reference traces
(a) Noise-free data; (b) Noisy data. Red: average trace as reference trace; Blue: 0 m-offset reference trace.

由上可知,对于无噪数据,近偏移距参考道的目标函数收敛情况优于远偏移距参考道;对于加噪数据,平均道作参考道的目标函数收敛情况优于最小偏移距参考道.利用不依赖子波、基于包络的FWI方法,建立的初始模型效果对于不同参考道和数据噪声并不敏感,这是因为利用包络目标函数的反演方法仅对模型的低频成份敏感.最终的反演结果得到的是宏观的背景速度,可以为接下来的常规FWI提供一个较好的初始速度模型.

3.2 不同目标函数的比较

由2.1节可知,不依赖子波、基于包络的FWI方法可以选择三种不同的目标函数.虽然三种目标函数均可以反演背景速度,但结果差别较大.图 9是选用同样的错误子波、在梯度模型的基础上、采用不同目标函数的反演结果,它们均是不依赖子波的,都对模型的背景成份进行了一定的恢复.但是,明显可以看出,包络对数目标函数对深层的恢复更好(图 9a),其次为包络目标函数(图 9b),而包络平方目标函数的效果最差(图 9c).在同样的迭代次数下,不依赖子波、基于包络平方目标函数的FWI方法对模型的深部速度基本不更新,收敛速度很慢.

图 9 不依赖子波、基于不同包络目标函数的FWI反演结果对比
(a) 包络对数目标函数; (b) 包络目标函数; (c) 包络平方目标函数.
Fig. 9 The inversion results of source-independent envelope-based FWI with different misfit functions
(a) Logarithmic envelope misfit function; (b) Envelope misfit function; (c) Squared envelope misfit function.

三种目标函数的不同,本质上是残差波场和伴随震源的不同,从而影响了求取的梯度,并最终影响了模型的更新.图 10图 11分别展示了不同目标函数在第一次迭代后第40炮的波场残差和相应的伴随震源,残差分别是 ln(Ecobs)-ln(Ec)、EcobsEc和(Ecobs)2-(Ec)2. 可以看出,采用包络对数目标函数时,来自深层的反射波对目标函数的贡献最大,并且这部分能量最终会体现在伴随震源中.因此,采用包络对数目标函数时,目标函数对深层更为敏感,并且这一部分能量将通过伴随震源反映到梯度中,从而提高模型深部的反演精度.相反,采用包络平方目标函数时,其残差波场中主要为初至波的能量,来自深层的反射波在其伴随震源中的贡献很弱,在求得的梯度中,模型深层的修改量很微弱,导致仅对模型的中浅部进行了较好的反演.包络目标函数则为介于两者之间.通过以上实验和分析可以得知,包络对数目标函数比较适用于地表观测的反射地震数据,可以对深层模型进行较好的反演;而包络平方目标函数更好地利用了初至波的信息,在井间观测数据的反演中可能会有更好的应用效果.

图 10 第一次迭代后,不同目标函数反演的波场残差
(a) 包络对数目标函数;(b) 包络目标函数;(c) 包络平方目标函数.
Fig. 10 The residual wavefields for different misfit functions after 1st iteration
(a) Logarithmic envelope misfit function; (b) Envelope misfit function; (c) Squared envelope misfit function.

图 11 第一次迭代后,不同目标函数反演的伴随震源
(a) 包络对数目标函数;(b) 包络目标函数;(c) 包络平方目标函数.
Fig. 11 The adjoint sources for different misfit functions after 1st iteration
(a) Logarithmic envelope misfit function; (b) Envelope misfit function; (c) Squared envelope misfit function.

最后,选取最小偏移距参考道和包络对数目标函数进行不依赖子波、基于包络的FWI建立初始速度模型,在此基础上进行常规的FWI,并与传统的基于包络的FWI反演结果对比(图 12).可以发现,错误的子波对于传统的基于包络的FWI反演影响是巨大的(图 12a),而以图 12a为初始模型的常规FWI结果当然非常糟糕(图 12b).而当采用了不依赖子波、基于包络的FWI方法后,反演的初始模型不再受错误子波的影响,较好地反映了模型速度变化的低波数成分(图 12c),而以图 12c为初始模型的常规FWI反演结果(图 12d)与准确速度模型(图 1a)非常接近,这一点也可以从沿深度抽线结果中清楚地看出(图 13).

图 12 子波错误时FWI反演结果对比
(a) 传统基于包络的FWI建立的初始模型;(b) 以(a)为初始模型的常规FWI反演结果;(c)不依赖子波、基于包络的FWI建立的初始模型;(d) 以(c)为初始模型的常规FWI反演结果.
Fig. 12 The comparison of different FWI results with an incorrect wavelet
(a) The initial model built by traditional envelope-based FWI; (b) The inversion result of traditional FWI using (a) as the initial model; (c) The initial model built by source-independent envelope-based FWI; (d) The inversion result of traditional FWI using (c) as the initial model.

图 13 本文方法的反演抽线结果
(a) 1000 m处速度抽线; (b) 5800 m处速度抽线. 绿色:准确模型;蓝色:本文方法建立的初始模型;红色:最终FWI反演结果.
Fig. 13 Depth velocity slices of inversion results
(a) X=1000 m; (b) X=5800 m. Green: true model; Blue: the inverted initial model by our method; Red: the final inversion result.
3.3 子波试验

为了进一步说明本文提出的方法在为常规FWI建立初始模型时并不依赖于子波,我们选取了不同参数的子波(均不同于原始观测数据中的子波)来进行测试,选用的子波参数如下:(1)主频20 Hz、 滤去5 Hz以下低频成分的Ricker子波,延时0.5 s,最大振幅2.0;(2)主频为12 Hz的高斯子波,延时0.15 s,最大振幅1.0.这里需要再次说明的是,合成地震数据时所用的是主频为8 Hz、延时0.15 s、最大振幅为1.0的Ricker子波.

图 14a图 14c是本文方法采用上述两种错误子波时的反演结果,可以看出,即使采用错误的子波,甚至在地震数据缺少低频信息的情况下,反演结果也都较好地恢复了模型中的背景成分.用图 14a图 14c作为初始速度模型,常规FWI的反演结果都比较好地刻画了真实模型的高低波数成分(图 14b图 14d).

图 14 本文方法采用两个错误子波时建立的初始模型以及相应的常规FWI反演结果
(a) 采用主频为20 Hz滤去5 Hz以下的Ricker子波建立的初始模型;(b)以(a)为初始模型的常规FWI反演结果;(c) 采用主频为12 Hz的高斯子波建立的初始模型;(d) 以(c)为初始模型的常规FWI反演结果.
Fig. 14 The inverted initial models using two incorrect wavelets by our method and the corresponding traditional FWI results
(a) The inverted initial model using 20 Hz Ricker wavelet without low frequency below 5 Hz; (b) The corresponding traditional FWI result using (a) as the initial model; (c) The inverted initial model using 12 Hz Gauss wavelet; (d) The corresponding traditional FWI result using (c) as the initial model.

为了验证本文方法反演的背景速度的可靠性,我们分别用四种速度模型对地震数据进行了逆时偏移,结果如图 15所示.偏移所用的四种速度模型分别是:准确模型的Gauss平滑结果(图 1b)、上述两种错误子波情况下本文方法的反演结果(图 14a图 14c)以及梯度模型(图 3).相比于梯度模型的偏移结果(图 15d),本文提出的方法反演的模型的成像结果(图 15bc)大大改善,与Gauss平滑模型的偏移结果(图 15a)基本一致,这说明了本方法在子波未知情况下,仍对模型中的低波数成分进行了很好的恢复,这样的反演结果当然就为接下来的常规FWI提供了很好的初始模型.

图 15 不同速度模型的逆时偏移结果对比
(a) 高斯平滑初始模型;(b) 主频为20 Hz滤去5 Hz以下低频Ricker子波建立的初始模型;(c) 12 Hz高斯子波初始模型;(d) 梯度初始模型.
Fig. 15 The comparison of RTM results with different velocities
(a) Gauss-window smoothed initial model; (b) The inverted initial model using 20 Hz Ricker wavelet without low frequency below 5 Hz; (c) The inverted initial model using 12 Hz Gauss wavelet; (d) Gradient initial model.
4 结论

(1)为进一步推动FWI的实用性,本文在前人 研究的基础上,提出了不依赖子波、基于包络的FWI反演方法.在子波未知、没有任何模型先验信息、以及地震数据缺少低频信息的情况下,这种方法可以为常规FWI提供一个可靠的初始速度模型.利用本文方法反演的初始模型,也可以作为叠前深度偏移的速度模型.

(2)比较了不同参考道的反演结果和归一化目标函数的收敛速度,认为近偏移距参考道优于远偏移距参考道.当数据含有噪音时,平均道作为参考道优于最小偏移距参考道.

(3)在不依赖子波、基于包络的FWI的不同目标函数中,包络对数目标函数为最优选择,对深层的反演效果最好.

(4)不依赖子波、基于包络的FWI反演方法能够克服传统FWI方法对震源子波的依赖.

附录A 不依赖子波、基于包络对数目标函数的梯度公式推导

不依赖子波、基于包络对数目标函数为公式(9),它对模型参数 p 求导,得

其中

h 为希尔伯特变换的滤波因子.

将(A2)和(A3)带入方程(A1)中并展开,得到如下四项:

其中第一项(A4)和第三项(A6),第二项(A5)和第四项(A7)推导过程类似.

下面推导(A4).(A4)可变为

(A8)进一步变换为

其中代表相关运算.将(A10)中的 ui,k p 写成褶积形式,得

(A11)进一步变换为

作类似的变换,第二项(A5)化为

对(A14)括号内的部分应用相关褶积定理,同时注意滤波因子 h 为奇函数,得

三、四项变换过程与上面类似.将变换后的结果带入原式,合并后得到最终结果:

对于包络目标函数、包络平方目标函数,它们对参数的梯度的推导过程与此类似,这里不再赘述.

参考文献
[1] 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.
[2] Brossier R, Operto S, Virieux J. 2009. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion. Geophysics, 74(6): WCC105-WCC118.
[3] Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion. Geophysics, 60(5): 1457-1473.
[4] Chi B X, Dong L G, Liu Y Z. 2013. Full waveform inversion based on envelope objective function. // 75th Annual International Conference and Exhibition, EAGE, Extended Abstracts.
[5] 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.
[6] Choi Y, Shin C, Min D J, et al. 2005. Efficient calculation of the steepest descent direction for source-independent seismic waveform inversion: An amplitude approach. Journal of Computational Physics, 208(2): 455-468.
[7] Choi Y, Alkhalifah T. 2011. Source-independent time-domain waveform inversion using convolved wavefields: Application to the encoded multisource waveform inversion. Geophysics, 76(5): R125-R134.
[8] Dong L G, Ma Z T, Cao J Z, et al. 2000. A staggered-grid high-order difference method of one-order elasticwave equation. Chinese Journal of Geophysics (in Chinese), 43(3): 411-419, doi: 10.3321/j.issn:0001-5733.2000.03.015.
[9] Lee K H, Kim H J. 2003. Source-independent full-waveform inversion of seismic data. Geophysics, 68(6): 2010-2015.
[10] Liu B L, Dong L G, Chi B X. 2014. Method of Source-independent Full Waveform Inversion[Ph. D. thesis] (in Chinese). Shanghai: Tongji University.
[11] Pratt R G, Shin C, Hick G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2): 341-362.
[12] Pratt R G. 1999. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model. Geophysics, 64(3): 888-901.
[13] Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies. Geophysics, 69(1): 231-248.
[14] Song Z M, Williamson P R, Pratt R G. 1995. Frequency-domain acoustic-wave modeling and inversion of crosshole data: Part II-Inversion method, synthetic experiments and real-data results. Geophysics, 60(3): 796-809.
[15] Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266.
[16] Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26.
[17] Wu R S, Luo J R, Wu B Y. 2014. Seismic envelope inversion and modulation signal model. Geophysics, 79(3): WA13-WA24.
[18] Zhou B, Greenhalgh S A. 2003. Crosshole seismic inversion with normalized full-waveform amplitude data. Geophysics, 68(4): 1320-1330.
[19] 董良国, 马在田, 曹景忠等. 2000. 一阶弹性波方程交错网格高阶差分解法. 地球物理学报, 43(3): 411-419, doi: 10.3321/j.issn:0001-5733.2000.03.015.
[20] 刘宝龙, 董良国, 迟本鑫. 2014. 不依赖于子波的全波形反演方法研究[硕士论文]. 上海: 同济大学.