地球物理学报  2015, Vol. 58 Issue (8): 2873-2885   PDF    
基于截断牛顿法的VTI介质声波多参数全波形反演
王义, 董良国    
同济大学海洋地质国家重点实验室, 上海 200092
摘要: 不同类别参数间的相互耦合使多参数地震全波形反演的非线性程度显著增加, 地震波速度与各向异性参数取值数量级的巨大差异也会使反演问题的性态变差.合理使用Hessian逆算子可以减弱这两类问题对反演的影响, 提高多参数反演的精度, 而截断牛顿法是一种可以比较准确地估计 Hessian 逆算子的优化方法.本文采用截断牛顿法在时间域进行了VTI介质的声波双参数同时反演的研究.不同模型的反演试验表明, 在VTI介质声波双参数同时反演中, 截断牛顿法比有限内存 BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno, L-BFGS)法能更准确地估计 Hessian 逆算子, 进而较好地平衡两类不同参数的同时更新, 得到了比较精确的反演结果.
关键词: 多参数     全波形反演     各向异性     耦合     Hessian逆算子     截断牛顿法    
Multi-Parameter full waveform inversion for acoustic VTI media using the truncated Newton method
WANG Yi, DONG Liang-Guo    
State Key Laboratory of Marine Geology, Tongji University, Shanghai 200092, China
Abstract: It is important to reconstruct both seismic velocities and anisotropy parameters by anisotropic full waveform inversion(FWI). But the trade-off between different parameter classes increases its nonlinearity severely. The significant difference of magnitude orders between seismic velocities and anisotropy parameters makes the inversion problem poorly conditioned. Judiciously incorporating the inverse Hessian operator within the FWI scheme can help to alleviate the negative impacts of these two problems on FWI. The truncated Newton(TN)method can better estimate the inverse Hessian operator. Therefore, we use the TN method to tackle these two problems and simultaneously reconstruct both seismic velocities and anisotropy parameters for acoustic VTI media.
The TN method has two advantages. On one hand, it does not need explicit Hessian matrix, which is prohibitive to build due to high computational cost and large memory requirements. To get the updated direction, the TN method solves Newton equation at each time of iteration with the conjugate gradient method(CG), which only uses Hessian-vector products. We compute these products by the finite difference method. On the other hand, the TN method can adjust flexibly the updated direction when the Hessian matrix is non-positive definite. The CG method at each step of the TN iteration is terminated as soon as a direction of negative curvature is generated. With this strategy the TN method may perform better than Limited-memory Broyden-Fletcher-Goldfarb-Shanno method(L-BFGS)for ill-conditioned problems.
We test the inversion capability of the TN method with two numerical experiments. In each experiment, the NMO velocity and anisotropy parameter eta are inverted simultaneously. The L-BFGS method is also tested for comparison. The first experiment is on a simple 2D VTI homogenous model with three different anomalies inside. The inversion is performed in the time domain with a full acquisition system. Inversion results show that the TN method can balance the updates of the two different parameters, and reconstruct more accurate models. In contrast, the L-BFGS method overestimates parameter eta but underestimates NMO velocity. The second experiment is on a more complex 2D model with a surface acquisition system. The two methods perform the same as they do in the first experiment.
The reason for their different performances may be found from the Hessian analysis part that we have finished before the two experiments. We design a small size(51 ×51)inclusion model and compute explicitly the Hessian matrix of the background initial model. The analyses show that the Hessian matrix has several negative eigenvalues and a large condition number(ill-conditioned matrix), which may result in the failure of the L-BFGS method. The TN method better overcomes these difficulties. Hence, in the context of simultaneous inversion of two parameters for acoustic VTI media, the TN method can account for more accurate inverse Hessian operator than the L-BFGS method and get more accurate inversion results. It is a promising method for multi-parameter FWI.
Key words: Multi-Parameter     Full waveform inversion     Anisotropy     Trade-off     Inverse Hessian operator     Truncated Newton method    
1 引言

地震全波形反演(Full waveform inversion,FWI)综合利用观测地震记录的波形信息来估计地下介质的物性参数(如密度、速度和各向异性参数等),在油气勘探、开发和地球动力学等领域都有广阔的应用前景(刘国峰等,2012石玉梅等,2014Fichtner et al.,2009).随着对复杂地下介质研究需求的增长和大规模科学计算能力的提升,FWI逐渐由单一纵波速度的反演向多个参数同时反演的方向发展(Virieux and Operto,2009).

地下介质的各向异性使地震波沿不同方位角以不同的速度传播,在地下存在各向异性情况时,使用各向同性FWI会得到错误的速度结构(Plessix and Cao,2011).各向异性参数在岩性识别和裂缝特征描述等方面均有重要作用.利用各向异性FWI能够更为精确地反演地下介质的速度和各向异性参数(张美根,2003Gholami,2012).但是,不同类别参数间存在的相互耦合效应增加了多参数FWI的非线性程度.同时,速度与各向异性参数取值数量级的巨大差异也会使反演问题的性态变差(Operto et al.,2013).

这两个问题都可以通过合理使用Hessian逆算子来部分解决.Hessian逆算子在FWI中起着重要作用.在单参数反演中,该算子可以对照明不足位置处的参数重新聚焦(Pratt et al.,1998).在多参数反演中,Hessian算子的非对角块矩阵中的元素描述了不同类别参数间的耦合效应(Operto et al.,2013),使用Hessian逆算子可以减弱这种耦合效应对FWI的影响(Métivier et al.,2014a).但由于巨大的计算和存储消耗,显式计算和存储Hessian及其逆算子是不现实的,这使得牛顿法和高斯牛顿法不能得到广泛应用(Virieux and Operto,2009).常规的梯度类方法如梯度法和非线性共轭梯度(Conjugate gradient,CG)法通过使用对角近似的Hessian逆算子对梯度进行预条件来加速反演的收敛过程(Shin et al.,2001Lee et al.,2010).L-BFGS法则更进一步,利用一定数目的模型和梯度信息构造正定矩阵来近似Hessian逆算子(Nocedal and Wright,2006).这种近似已不再是简单的对角近似,使得L-BFGS法的稳健性和收敛速度都优于常规的梯度类方法,在FWI中得到了广泛应用(Guitton and Symes,2003Brossier et al.,2009). 但是在很多病态反演问题中,Hessian矩阵的非正定 性使得L-BFGS法的性能变差,截断牛顿(Truncated Newton,TN)法则可能更加有效(Santosa and Symes,1988Nash and Nocedal,1991Métivier et al.,2013). 它在每步迭代中通过线性CG法求解牛顿方程来获取模型增量,以此来更加准确地利用Hessian逆算子的信息.

TN法在地震FWI领域逐渐得到研究和应用(Santosa and Symes,1988Akcelik,2002Epanomeritakis et al.,2008). 最近两年,Métivier等(2014b)在频率域进行了纵波速度的TN法全波形反演的研究,并用试验说明了存在强多次反射波情况下,TN法比截断高斯牛顿法、L-BFGS法以及非线性CG法能够更加准确地估计Hessian逆算子.在多参数方面,Bae等(2012)采用截断高斯牛顿法研究了各向同性声弹耦合介质的多参数反演,得到了纵横波速度和密度,但并未进一步研究TN法的反演性能.

在VTI介质FWI方面,Lee等(2010)用预条件梯度法反演了4个弹性系数,他们更注重采用分步反演的策略来克服多参数反演中的病态性.Plessix和Cao(2011)在声波VTI介质中分别对NMO速度和水平速度使用L-BFGS法进行单参数反演,前者的反演结果偏高,而用后者的结果转换出的等价各向异性参数的结果与真实模型相差较多.Gholami(2012)运用L-BFGS法进行了声波和弹性波的VTI 介质多参数同时反演.他使用参数标准化(Normalization)将不同类别的参数调整到相同的数量级上,使反演过程更加稳健,但得到的往往是过度估计的速度参数和低估的各向异性参数.Métivier等(2014b)运用TN法在频率域进行了声波VTI介质的纵波速度的单参数全波形反演.

本文将截断牛顿法应用到VTI介质声波的双参数同时反演中,并分析了它在反演中的性能.在介绍了TN法的基本原理以及分析了Hessian矩阵的性态之后,通过时间域不同模型的合成数据数值反演试验,对比研究了TN法与L-BFGS法的性能.

2 截断牛顿法原理

设最小二乘目标函数为

其中,m 为模型参数,d obsf(m)分别为观测和模拟记录,f是联系模型参数和模拟记录的算子,在时间域*代表转置.求解(1)式常采用局部迭代优化方法,其模型更新公式为
其中,p kαk分别为第k次迭代的更新方向和步长.对目标函数的二阶Taylor展开进行最小化,可以得到牛顿法的更新方向所满足的方程:
其中,▽ m E(m k)和Hessian矩阵 H(m k)分别为目标函数对参数 m k的梯度和二阶导数.

由于显式Hessian及其逆矩阵的计算和存储消耗巨大,在大规模FWI中求取精确的牛顿更新方向是不现实的,使得具有局部快速收敛优势的牛顿法无法获得广泛应用.TN法则可以克服这个困难,它在每步迭代中使用线性CG法求解(3)式来获取比较准确的牛顿更新方向,而线性CG法只需整体使用Hessian矩阵和已知向量的乘积.每个Hessian向量乘积都可以利用有限差分法(Nocedal and Wright,2006)或者二阶伴随状态法(Fichtner and Trampert,2011Métivier et al.,2013)通过四次正演计算得到.本文使用有限差分法计算Hessian向量乘积,计算公式为

其中,v 为已知向量,h为一个很小的正数.相比于牛顿法,TN法的计算量已经可以承受,在大规模反演问题中获得了应用并得到了较好的反演结果(Nash and Nocedal,1991Métivier et al.,2014b).

下面的算法1给出了TN法的流程,内层while循环代表了线性CG迭代过程.为避免计算量过大,CG迭代需要设定合理的停止条件(Wang等,2013). 当CG迭代中遇到负曲率方向(即算法1中β1 <0),表示当前的Hessian矩阵为非正定矩阵,继续迭代会得到不下降的更新方向,因而将CG迭代停止.通过这种CG过程,TN法主要利用了Hessian逆算子 的大特征值的作用,因而对模型更新的计算有内在的正则化作用(Kaltenbacher et al.,2008).

本文采用的VTI介质声波方程的模型参数为(Vnη,δ),具体方程和梯度公式见附录A.其中,η=(ε-δ)/(1+2δ),NMO速度Vn=Vp$\sqrt{1+2\delta }$,Vp为纵波速度,ε,δ为各向异性参数(Thomsen,1986).Alkhalifah和Plessix(2014)通过辐射模式分析指出该组参数间存在着有限的耦合效应,更有利于反演.由于δ对数据的敏感性很低,反演中将其固定为真实值(Gholami et al.,2011Plessix and Cao,2011).因而,下文中反演的模型参数为Vnη.

算法1 截断牛顿法

给定初始模型 m 0,令k=0

while 反演停止条件未满足 do

计算梯度 ▽ m E(m k);

l=0,d =0,r = ▽ m E(m k),x =- r ;

while CG法停止条件未满足 do

计算 H(m k)x ;

β1=(H(m k)xx); β2=‖ r2;

d = d +(β2/β1)x ;

r = r +(β2/β1)H(m k)x ;

x =- r +(‖ r2/β2)x ;

l=l+1;

end while

更新方向 p k= x ;

计算步长αk;

更新模型 m k+1= m k+αk p k;

k=k+1;

end while

3 Hessian矩阵性态分析

Hessian矩阵的性态会影响到局部迭代优化方法的性能.当Hessian矩阵的条件数很大时(病态矩阵),会导致局部迭代优化方法收敛缓慢甚至不收敛(Nocedal and Wright,2006).下面我们采用一个简单的模型来分析参数Vnη的Hessian矩阵的结构和性态.

考虑圆形异常模型,如图 1所示.模型网格剖分为51×51,间距均为10 m.η真实模型的四周是宽度为100 m的椭圆各向异性层,震源和检波器置于此层内可以压制S波的影响(Plessix and Cao,2011).在模型的每条边上分别均匀布置12个炮点,炮间距为40 m,每炮都采用四周观测方式接收数据.震源使用主频为8 Hz的Ricker子波,时间采样间隔为1 ms.初始模型为没有异常体的背景模型,我们使用有限差分法近似计算出初始模型中心部分(31×31)的Hessian矩阵(Nocedal and Wright,2006),如图 2a所示.

图 1 真实模型(a)Vn(b)η Fig. 1 True model(a)Vn(b)η

图 2 (a)初始模型的Hessian矩阵及其(c)特征值分布,(b)为单独显示的Hessian矩阵的右上角块矩阵 Fig. 2 (a)Hessian matrix of initial model and its(c)Eigenvalue distribution, (b)Display for the right-upper block matrix of the Hessian

由于考虑了两类参数,Hessian矩阵由4个961×961 大小的分块矩阵组成,在图 2a中用黄色虚线进行了划分.左上角块矩阵中的非对角元素代表两个不同位置处的Vn对目标函数的影响,对角线上的元素则代表该位置的Vn对目标函数的二阶影响,右下角块矩阵类似地对应着参数η的作用.而非对角块矩阵(右上和左下)中的元素则描述了不同类别参数间的相互耦合效应.图 2b单独给出了Hessian矩阵的右上角块矩阵,它的主对角元素比较突出,说明耦合效应主要集中在相同位置的Vnη上.从图 2a可以看出,右下角块非常突出.经计算Hessian矩阵的条件数为3.58×1010,病态比较严重.图 2c给出了该矩阵的特征值分布,该矩阵没有0特征值(绝对值均大于1×10-5),是非奇异矩阵.可以看出有负特征值存在且部分幅度较大,因而该Hessian矩阵不是正定矩阵.若此时采用L-BFGS法进行参数Vnη的同时反演,由于它构造的近似Hessian逆算 子为正定矩阵(Nocedal and Wright,2006),可能会收敛缓慢甚至失败.如果使用TN法则可能较好地处理这类病态问题.下面我们通过数值试验进行检验.

4 数值试验

本节采用两个不同的二维模型测试TN法的反演性能.每个试验中都将TN法与L-BFGS法(使用最近8代的模型和梯度向量)的收敛速度和反演结果进行了对比分析.

4.1 简单模型试验

真实模型包含了圆形(半径400 m)、正方形(边长800 m)和三角形(底边长为800 m的等腰三角形)三个异常体,如图 3所示.η真实模型的四周各有厚度为100 m的椭圆各向异性层,反演时固定不变.网格剖分为401×401,间距均为10 m.为便于叙 述,记该模型为CST模型.在模型的每条边上分别 均匀布置16个炮点,炮间距为240 m,每炮都采用四周观测方式接收记录.震源采用主频为8 Hz的Ricker子波,时间采样间隔为1 ms.两类参数的初始模型均为没有异常体的背景模型.

图 3 CST真实模型(a)Vn(b)η Fig. 3 True CST model(a)Vn(b)η

两种方法的归一化的目标函数的下降曲线如图 4所示,其横坐标为迭代次数.可以看出,TN法目标函数(红色实线)随迭代次数的下降速度明显快于L-BFGS法(黑色实线),说明了TN法在单次迭代中能利用更多的Hessian矩阵信息.但从第2部分的原理介绍可知,下降速度的加快需要以增加模拟计算次数为代价.表 1给出了两种方法在反演过程中的计算数据对比.表中第2行数据表明在第25次迭代时,TN法的目标函数下降到0.0039,远低于L-BFGS法的0.1762,但使用的正演次数也增加了很多.这反应出TN法的总体计算消耗较大.同时表中第1行给出了TN法第6代和L-BFGS法第12代的对比,它们都使用了48次正演.TN法的目标函数已经下降到0.052,但L-BFGS法只下降到0.2332,因而在同等计算消耗下,TN法下降更快.

图 4 CST模型的归一化的目标函数的下降曲线 Fig. 4 The decrease curve of normalized misfit function for CST model

表 1 CST模型试验中TN法和L-BFGS法的对比 Table 1 Comparison of TN method and L-BFGS method for the CST model testing

反演结果的好坏是评判优化方法的关键因素.图 5给出了两种方法各自对应于表 1第1行的反演结果.可以看出,两组反演结果差异较大:在3个异常体位置,L-BFGS法对模型η的更新过度,对Vn的更新则很小,有可能收敛不到全局最优模型;而TN法对两类参数都有较好更新,反演结果在异常体结构和幅度上都已经比较接近真实模型.接下来进一步对比两种方法的第25代反演结果,如图 6所示.在TN法的结果中,Vn模型(图 6c)的3个异常体的结构比第6代时更加清晰,幅度更逼近真实模型,η模型(图 6d)的幅度比第6代也有细微改善.对于L-BFGS法,模型Vn的更新仍旧非常缓慢,而模型η偏离真实模型的程度加重.

图 5 L-BFGS法第12代的(a)Vn(b)η 与TN法第6代的(c)Vn(d)η Fig. 5 Inversion results at the 12th iteration of L-BFGS method(a)Vn(b)η and those at the 6th iteration of TN method(c)Vn(d)η

图 6 第25代反演结果. L-BFGS法的(a)Vn(b)η,TN法的(c)Vn(d)η Fig. 6 Inversion results at the 25th iteration . L-BFGS method(a)Vn(b)η and TN method(c)Vn(d)η

为了更清晰地观察两种方法的差别,我们对比了反演结果在不同位置(图 6c中白色点线所示)的模型参数曲线.图 7给出了对应于图 5的反演结果在深度2 km的横向对比曲线,图 8则是相应反演结果在地表距离1 km,2 km和3 km处的垂向对比曲线.这些曲线也表明TN法第6代的反演结果,尤其是η模型,已经比较接近真实模型.而L-BFGS法第12代所反演的η模型已较大幅度偏离真实模型. 两种方法第25代反演结果(即图 6中模型)的对比曲线分别由图 9图 10提供.可以看出,TN法的Vn模型曲线已经非常接近真实模型,η模型曲线虽然在异常体的区分上仍有不足,但在异常体处的幅度与真实模型基本一致.对于L-BFGS法,迭代次数的增多并未使反演结果好转:η模型与真实模型的偏差增大,Vn模型几乎没有更新.这说明L-BFGS法很可能收敛到远离真实模型的极小值点,甚至迭代发散.

图 7 L-BFGS法第12代和TN法第6代反演结果在深度2 km的横向曲线对比黑色和红色实线分别为真实和初始模型,蓝色和绿色虚线分别为L-BFGS法和TN法的结果. Fig. 7 Horizontal profile comparisons at depth 2 km between inversion results at the 12th iteration of L-BFGS method(blue dash line) and those at the 6th iteration of TN method(green dash line)The black and red solid lines represent true and initial models,respectively.

图 8 L-BFGS法第12代和TN法第6代反演结果在地表不同距离的垂向曲线对比黑色和红色实线分别为真实和初始模型,蓝色和绿色虚线分别为L-BFGS和TN法的结果. Fig. 8 Vertical profile comparisons at different surface distances between inversion results at the 12th iteration of L-BFGS method(blue dash line) and those at the 6th iteration of TN method(green dash line)The black and red solid lines represent true and initial models,respectively.

图 9 两种方法第25代反演结果在深度2 km处的横向曲线对比黑色和红色实线分别为真实和初始模型,蓝色和绿色虚线分别为L-BFGS法和TN法的结果. Fig. 9 Horizontal profile comparisons at depth 2 km between inversion results at the 25th iteration of both L-BFGS(blue dash line) and TN method(green dash line)The black and red solid lines represent true and initial models,respectively.

图 10 两种方法第25代反演结果在地表不同距离的垂向曲线对比黑色和红色实线分别为真实和初始模型,蓝色和绿色虚线分别为L-BFGS和TN法的结果. Fig. 10 The black and red solid lines represent true and initial models,respectively.The black and red solid lines represent true and initial models,respectively.

综合以上对比可知,TN法能够较好地平衡Vn和η的同时更新,得到了比较准确的反演结果.而L-BFGS法则过度集中在一类参数的更新上,反演结果较差.由第3部分的分析可知导致两种方法性能差异的原因可能在于此时Hessian矩阵的性态比较差.L-BFGS法未能对Hessian逆矩阵进行比较准确的近似,从而没能准确解释两类参数间的耦合效应,将Vn引起的数据误差转移到了η上,造成其过度更新.这也反应了多参数同时反演的复杂性.TN法则比L-BFGS法准确地估计了Hessian逆矩阵,展现出良好性能.

4.2 复杂模型试验

真实模型如图 11(a,b)所示,网格剖分为301×71,间距均为20 m.在地表激发56炮,炮间距为100 m,在地表所有网格点处接收地震记录.采用主频为9 Hz的Ricker子波,时间采样间隔为2.5 ms.对真实模型进行高斯平滑得到两类参数的初始模型,如图 11(c,d)所示.试验中也进行了L-BFGS法的反演.

图 11 真实模型(a)Vn(b)η和初始模型(c)Vn(d)η Fig. 11 True models(a)Vn(b)η and initial modes(c)Vn(d)η

两种方法的归一化的目标函数随迭代次数的下 降曲线由图 12给出.类似于CST模型试验,TN法 的目标函数(红色)随迭代次数的下降同样明显快于 L-BFGS法(黑色).图 13分别给出了TN法的最终反演结果(第46代,正演次数超过1000,目标函数值为0.0112)与L-BFGS法在相同迭代次数时(第46代,正演次数为184,目标函数值为0.35)的反演结果.可以看出,TN法反演的Vn模型(图 13c)已经与真实模型比较接近,η模型中(图 13d)也得到了比较准确的界面位置.L-BFSG法则主要对η模型进行了过度更新,对Vn模型改善很小.图 14分别给出了两种方法的反演结果在地表 1.5 km,3 km和4.5 km处(具体位置如图 13c中的白色点线所示)的垂向曲线对比.从中可以更清晰地看出,TN法真实地反演出了Vn模型的结构和幅度.虽然η模型的幅度与真实模型仍有差距,但能够得到比较准确的界面结构.L-BFSG法的模型曲线表明,它对Vn模型的更新非常有限,而主要对η模型的浅部进行了过度更新,因此该方法未能收敛到全局最优模型.

图 12 复杂模型的归一化的目标函数下降曲线 Fig. 12 The decrease curve of normalized misfit function for complex model

图 13 第46代反演结果 L-BFGS法的(a)Vn(b)η,TN法的(c)Vn(d)η. Fig. 13 Inversion results at the 46th iteration L-BFGS method(a)Vn(b)η and TN method(c)Vn(d)η.

图 14 两种方法第46代反演结果在地表不同距离的垂向曲线对比黑色和红色实线分别为真实和初始模型,蓝色点线和绿色虚线分别为L-BFGS和TN法的结果. Fig. 14 Vertical profile comparisons at different surface distances between inversion results at the 46th iteration of both L-BFGS(blue dotted line) and TN method(green dash line)The black and red solid lines represent true and initial models,respectively.

以上复杂模型试验同样说明了TN法能够比L-BFGS法更为准确地估计Hessian逆矩阵,较好地平衡了两类不同参数的同时更新.

5 结论与讨论

本文采用TN法开展了时间域声波VTI介质的双参数同时反演的研究.地震波速度与各向异性参数间存在着耦合效应,增加了多参数FWI的复杂程度.同时,地震波速度和各向异性参数间取值数量级的巨大差异也使VTI介质多参数反演问题的性态变差.合理使用Hessian逆算子能够减弱这两类问题对FWI的影响.不同模型的反演试验表明,在地震波速度与各向异性参数的同时反演中,L-BFGS法对各向异性参数进行了过度更新,对地震波速度的更新量很小,反演效果较差,这也验证了多参数同时反演问题的复杂性.而TN法则较好地平衡了这两类不同参数的同时更新,得到了比较准确的反演结果,这说明了TN法比L-BFGS法能更准确地估计Hessian逆算子.因而,TN法在多参数反演中的应用潜力较大.

从文中的反演试验可以看出,虽然TN法得到了比较好的反演结果,对多参数FWI有重要的理论价值,但由于每步迭代中的CG循环都需要增加额外的模拟过程,也使方法本身付出了较大的计算代价,这也是限制TN法大规模实际应用的主要原因. 所以,采用适当的手段(如对CG循环进行预条件,或与其他方法混合迭代等)来提高TN法的计算效率需要进一步研究.

多参数FWI是一项复杂的研究工作.开发高效稳健的优化方法是其中非常重要的一个方面.同时,也需要注重优化方法与反演策略的有效结合.例如,将不同类型参数的数量级调整到相同取值范围的参数标准化策略可以改善Hessian矩阵的性态,因而可能会加速TN法的收敛.此外,VTI介质声波方程有多种不同的参数化方式,每种方式的Hessian算子的性态可能不同.因此,TN法在不同参数化方式下的性能表现也需要具体研究.

附录A

文中采用二维声波VTI方程(Plessix and Cao,2011):

其中pn=pv$\sqrt{1+2\delta }$,phpv分别为负的水平和垂直应力,s为震源.为压制横波假象,设定震源和检波器附近为椭圆各向异性介质,即η=0,从而pn=ph.记p*为观测记录,合成记录取为p=S(α1ph+α2pn),满足α1+α2=1,其中S为取样算子.采用有限差分法进行模拟,四周均采用NPML边界条件(McGarry and Moghaddam,2009).

目标函数对Vnη的梯度分别为(Plessix and Cao,2011):

其中λhλn满足以下的伴随状态方程:

参考文献
[1] Akcelik V. 2002. Multiscale Newton-Krylov methods for inverse acoustic wave propagation. Pittsburgh:Civil and Environmental Engineering, Carnegie Mellon University.
[2] Alkhalifah T, Plessix R.2014. A recipe for practical full-waveform inversion in anisotropic media:An analytical parameter resolution study. Geophysics, 79(3):R91-R101.
[3] Bae H S, Pyun S, Chung W, et al. 2012. Frequency-domain acoustic-elastic coupled waveform inversion using the Gauss-Newton conjugate gradient method. Geophysical Prospecting, 60(3):413-432.
[4] 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.
[5] Epanomeritakis I, Akelik V, Ghattas O, et al. 2008. A Newton-CG method for large-scale three-dimensional elastic full-waveform seismic inversion. Inverse Problems, 24(3):034015.
[6] Fichtner A, Kennett B L N, Igel H, et al. 2009. Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods. Geophysical Journal International, 179(3):1703-1725.
[7] Fichtner A, Trampert J. 2011. Hessian kernels of seismic data functionals based upon adjoint techniques. Geophysical Journal International, 185(2):775-798.
[8] Gholami Y, Brossier R, Operto S, et al. 2011. Acoustic VTI full waveform inversion:sensitivity analysis and realistic synthetic examples.//81th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2465-2470.
[9] Gholami Y. 2012. Two-dimensional seismic imaging of anisotropic media by full waveform inversion. Nice:Université de Nice-Sophia Antipolis.
[10] Guitton A, Symes W W. 2003. Robust inversion of seismic data using the Huber norm. Geophysics, 68(4):1310-1319.
[11] Kaltenbacher B, Neubauer A, Scherzer O. 2008. Iterative Regularization Methods for Nonlinear Ill-posed Problems. New York:Walter de Gruyter.
[12] Lee H -Y, Koo J M, Min D -J, et al. 2010. Frequency-domain elastic full waveform inversion for VTI media. Geophysical Journal International, 183(2):884-904.
[13] Liu G F, Liu H, Meng X H, et al. 2012. Frequency-related factors analysis in frequency domain waveform inversion. Chinese J. Geophys.(in Chinese), 55(4):1345-1353, doi:10.6038/j.issn.0001-5733.2012.04.030.
[14] McGarry R, Moghaddam P. 2009. NPML boundary conditions for second-order wave equations.//79th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3590-3594.
[15] Métivier L, Brossier R, Virieux J, et al. 2013. Full waveform inversion and the truncated Newton method. J. Sci. Comput., 35(2):B401-B437.
[16] Métivier L, Brossier R, Operto S, et al. 2014a. Multi-parameter FWI-an illustration of the Hessian operator role for mitigating trade-off between parameter classes. 6th EAGE Saint Petersburg International Conference & Exhibition. Saint Petersburg, Russia.
[17] Métivier L, Bretaudeau F, Brossier R, et al. 2014b. Full waveform inversion and the truncated newton method:quantitative imaging of complex subsurface structures. Geophysical Prospecting, 62(6):1353-1375.
[18] Nash S G, Nocedal J. 1991. A numerical study of the limited memory BFGS method and the truncated Newton method for large-scale optimization. SIAM J. Optim., 1(3):358-372.
[19] Nocedal J, Wright S J. 2006. Numerical optimization. New York:Springer.
[20] Operto S, Brossier R, Gholami Y, et al. 2013. A guided tour of multiparameter full waveform inversion for multicomponent data:from theory to practice. The Leading Edge,32(9):1040-1054.
[21] Plessix R.2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophysical Journal International, 167(2):495-503.
[22] Plessix R, Cao Q. 2011. A parametrization study for surface seismic full waveform inversion in an acoustic vertical transversely isotropic medium. Geophysical Journal International, 185(1):539-556.
[23] Pratt R G, Shin C, Hicks. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion. Geophysical Journal International, 133(2):341-362.
[24] Santosa F, Symes W W. 1988. Computation of the Hessian for least-squares solutions of inverse problems of reflection seismology. Inverse Problems, 4(1):211-233.
[25] Shi Y M, Zhang Y, Yao F C, et al. 2014. Methodology of seismic imaging for hydrocarbon reservoirs based on acoustic full waveform inversion. Chinese J. Geophys.(in Chinese), 57(2):607-617, doi:10.6038/cjg20140224.
[26] Shin C, Jang S, Min D -J. 2001. Improved amplitude preservation for prestack depth migration by inverse scattering theory. Geophysical Prospecting, 49(5):592-606.
[27] Thomsen L. 1986. Weak elastic anisotropy. Geophysics, 51(10): 1954-1966.
[28] Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6):WCC1-WCC26.
[29] Wang Y, Dong L G, Liu Y Z. 2013. Improved hybrid iterative optimization method for seismic full waveform inversion. Applied Geophysics, 10(3):265-277.
[30] Zhang M G, Wang M Y, Li X F, et al. 2003. Full wavefield inversion of anisotropic elastic parameters in the time domain. Chinese J. Geophys.(in Chinese), 46(1):94-100.
[31] 刘国峰,刘洪,孟小红等.2012.频率域波形反演中与频率相关的影响因素分析.地球物理学报,55(4):1345-1353,doi:10.6038/j.issn.0001-5733.2012.04.030.
[32] 石玉梅,张研,姚逢昌等.2014.基于声学全波形反演的油气藏地震成像方法.地球物理学报,57(2):607-617,doi:10.6038/cjg20140224.
[33] 张美根,王妙月,李小凡等.2003.时间域全波场各向异性弹性参数反演.地球物理学报,46(1):94-100.