地球物理学报  2019, Vol. 62 Issue (5): 1809-1823   PDF    
各向异性介质弹性波多参数全波形反演
刘玉柱, 黄鑫泉, 万先武, 孙敏傲, 董良国     
同济大学海洋地质国家重点实验室, 上海 200092
摘要:各向异性介质弹性波方程全波形反演过程中多参数之间的相互耦合,使得弱参数在反演过程中难得到理想的结果.本文以VTI介质为例,在各参数辐射模式分析的基础上,基于改进的散射积分算法实现目标函数梯度的直接求取,进一步构建高斯牛顿方向,实现Hessian矩阵的有效利用,以考虑Hessian矩阵非主对角线元素包含的各参数间的耦合效应,在不使用任何反演策略的情况下实现高精度的VTI介质弹性波方程多参数同步反演.同时,该方法在计算过程中无需存储庞大的核函数矩阵,且无需传统截断牛顿法中额外的正演计算,因此内存占用小,计算效率高.本文数值试验验证了该方法的有效性,为各向异性多参数全波形反演提供了一种新的解决方案.
关键词: 全波形反演      各向异性      弹性波方程      多参数      高斯-牛顿方向      改进的散射积分算法     
Elastic multi-parameter full-waveform inversion for anisotropic media
LIU YuZhu, HUANG XinQuan, WAN XianWu, SUN MinAo, DONG LiangGuo     
State Key Laboratory of Marine Geology, Tongji University, Shanghai 200092, China
Abstract: The trade-off effects between parameters in anisotropic media, which have been studied in radiation patterns analysis, makes it difficult to get acceptable weak parameter results during full-waveform inversion (FWI) based on an elastic wave equation. In this paper, to reduce the trade-off effects between parameters, we apply the truncated Gauss-Newton method to anisotropic FWI, taking VTI media for example. All parameters have been revealed simultaneously without using any strategy. Based on the improved scattering-integral (SI) algorithm, the steepest-descent direction can be directly calculated by vector operations, without storing the huge Fréchet kernels in memory beforehand. The Gauss-Newton descent direction iteration needs just two vector operations without any extra forward computations. Numerical experiments in the final part prove the effectiveness of the proposed method which opens a new approach on simultaneously constructing anisotropic parameters in FWI.
Keywords: Full-waveform inversion    Anisotropy    Elastic wave equation    Multi-parameter    Gauss-Newton method    Improved scattering-integral algorithm    
0 引言

地震全波形反演(FWI)是近三十年来发展起来的高分辨率地震数据处理方法,该方法利用了叠前地震数据的运动学和动力学信息重建地下属性,具有揭示复杂地质背景下的构造与岩性细节信息的能力.

自从Tarantola(1984)在最小二乘目标泛函框架下提出时间域FWI方法以来,FWI理论迅速发展.Pratt(1999)在Tarantola时间域FWI框架的基础上提出了频率域FWI方法.频率域FWI仅需要匹配观测数据中的若干有限离散频率成分就可以建立高精度的参数模型,推动了FWI向实际应用的发展.Sirgue等(2009)将频率域FWI应用于北海Valhall油田,获得了高精度地下速度结构,并在结果中清晰地展现了地下古河道分布和气云边界,揭示了FWI的巨大应用潜力.Luo等(2018)使用阻尼瞬时相位反演提供了良好的初始模型并帮助全波形反演取得更好的结果.基于声波近似的FWI已经进入实用化尝试阶段,但声波近似FWI的反演参数只包含纵波速度、密度及黏声参数.相对而言,弹性波全波形反演(EFWI)能更准确地描述实际地震波传播现象,理论上能够同时反演纵、横波速度及各向异性参数.最初的FWI并没有考虑地下介质的各向异性.但在实际应用中,由于忽略各向异性带来的运动学上的误差,偏移上表现为反射波归位不准,绕射波收敛不彻底,能量不聚焦,这种情况在大偏移距、宽方位数据的处理过程中表现得尤为明显(Barnes et al., 2008).通过考虑各向异性能够有效提高FWI的反演精度(Plessix,2009).因此,考虑各向异性的EFWI具有广阔的应用前景.

目标函数梯度的构建是FWI的关键.理论上,目标函数的梯度可以表达为数据残差沿着敏感核函数(Fréchet核函数)的投影.Woodward(1992)分别基于Born近似与Rytov近似给出了两种Fréchet核函数表达式,并将它们命名为Born波路径与Rytov波路径.但是,存储整个核函数以计算目标函数梯度甚至Hessian需要占用巨大的内存空间,现有的计算机软硬件条件无法满足要求.目前,伴随状态法被广泛用于构造FWI梯度,从而避免了核函数的计算与存储.所以,这两种核函数一直没能被直接应用于FWI研究中,而是被广泛应用于波动方程层析成像中,如Dahlen等(2000)提出的有限频层析成像,Spetzler和Snieder(2004)Liu等(2009)提出的菲涅尔体地震层析成像等.FWI考虑了波场的全部信息,但是由于弹性波传播的复杂性导致FWI具有强烈的非线性.在弹性波多参数反演过程各参数之间相互耦合,同时各参数对波场影响的敏感性也有强弱之分.通过使用不同的反演策略能一定程度上提高反演结果的精度(Prieux et al., 2013杨积忠等,2014刘玉柱等,2015),但是反演策略并没有从理论上考虑不同参数间的耦合问题,在反演过程中引入Hessian算子的逆,是解决这一问题的重要途径(Operto et al., 2013).

迭代方向计算中只利用目标函数一阶梯度的被称为一阶方向,同时利用了目标函数二阶梯度(Hessian)的称为二阶方向,即牛顿方向.精确Hessian矩阵由两项构成,利用精确Hessian的方向称为全牛顿方向.为简化计算同时考虑到多次散射能量相对较弱,精确Hessian的第一项(近似Hessian)被广泛应用,其对应的迭代方向称为高斯-牛顿方向.在多参数全波形反演问题中,Hessian矩阵的主对角线元素代表波传播过程的几何扩散效应,非主对角线元素包含了各参数间的耦合效应.因此,多参数FWI中Hessian的利用非常关键(Pan et al., 2016Sun et al., 2017).

Sheen等(2006)在时间域弹性波高斯-牛顿法FWI研究中将近似Hessian矩阵表达为列向量与行向量的乘积求和,这样就不需要存储核函数即可实现近似Hessian矩阵的计算,但存储该近似Hessian矩阵本身则需要更大的存储空间.因此,Choi等(2007)基于Hessian矩阵主对角线元素占优的特点,首先提出利用对角Hessian构造拟牛顿方向,但该拟牛顿方向误差较大且无法考虑参数耦合.Fichtner和Trampert(2011)采用二阶伴随状态法实现了精确Hessian的计算.Hessian矩阵作为庞大的病态矩阵,其逆矩阵无法直接求解.为了近似估计Hessian矩阵的逆矩阵,Sheen等(2006)利用l-BFGS方法实现了近似Hessian逆的迭代计算,其计算精度取决于迭代次数与初始Hessian的取值.Métivier等(2012)同样采用二阶伴随状态法实现了l-BFGS迭代算法中Hessian-向量乘计算,进而得到高斯-牛顿方向.Liu等(2015)提出了一种基于Born敏感核函数的改进的散射积分全波形反演方法,该方法通过对核函数-向量乘进行分解累加运算实现目标函数梯度的直接求取,且无需存储庞大的核函数矩阵;同时将Hessian-向量乘表达为两次核函数-向量乘,避免了Métivier等提出的截断牛顿法中额外的正演计算.

刘玉柱等(2015)基于改进的散射积分算法,并利用最速下降法实现了频率域拟声波方程VTI介质多参数全波形反演.相比于拟声波方程,各向异性弹性波方程能够更加准确地描述实际地震波传播的动力学特征,同时正演过程也更加稳定.Zhou和Greenhalgh(2009)推导并分析了一般情况下各向异性介质的敏感核函数.Lee等(2010)在伴随状态法的框架下实现了弹性波VTI介质多参数全波形反演,但未考虑多参数的耦合问题.基于上述研究背景,有必要尝试基于改进的散射积分法的VTI介质弹性波多参数全波形反演,并探讨Hessian矩阵在各向异性弹性波多参数反演中的作用.

1 方法原理

FWI利用了叠前地震数据的运动学和动力学信息重建地下属性,其反演过程为:求解能够使得理论合成数据和观测数据最佳匹配所对应的模型参数.FWI的反演过程用数学语言可以表达为求解特定目标泛函下的最优化问题.二维弹性波二范数意义下的目标函数可表达为:

(1)

(2)

其中,m为二维各向异性介质Thomsen参数表征下的vpvsεδ模型参数,Δux= uxobs- uxcal、Δuz= uzobs- uzcal(uxcaluzcal分别为理论合成波场的xz分量,uxobsuzobs分别为实测数据的xz分量)分别为xz分量的波场残差,†为共轭转置算子.

目标函数(1)是一个强烈的非线性反演问题,通常采用局部最优化方法求解.为了使目标函数取极小值,每一轮迭代的模型参数mn+1的更新为上一轮迭代的模型参数mn与模型参数扰动量Δm的和.Δm可以表达为更新方向p乘以步长:

(3)

其中α为步长.依据不同的优化算法选择不同的更新方向p.目标函数对模型参数求一阶偏导数并取负号,即可得到目标函数的最速下降方向p1,如公式(4)所示.目标函数对模型参数进行二阶泰勒展开,可推导得到牛顿方向p2的表达式,如公式(5)所示(Liu et al., 2015).

(4)

(5)

式中,∂uxcal/∂m、∂uzcal/∂m分别为模型参数扰动对波场的xz分量的影响,即波场的xz分量的敏感核函数KxKz.∂2∅(m)/ ∂m2为目标函数对模型参数m的二阶偏导数,即Hessian矩阵H.若不考虑高阶散射,舍弃Hessian矩阵的高阶项,可用Ha来进行近似,且满足

(6)

将式(2)代入式(4)、(5),并用核函数表达为:

(7)

(8)

其中,p1vpp1vsp1εp1δp2vpp2vsp2εp2δ分别为vpvsεδ的一阶方向和二阶方向. 分别为vpvsεδ在波场xz分量上的敏感核函数.

在牛顿优化算法中,Hessian矩阵作为庞大的病态矩阵,其逆矩阵很难直接求解.为避免Hessian矩阵的求逆运算,可以将方程Hap2= p1看作是一个新的最优化求解问题(Métivier et al., 2012),p2作为未知数,并建立新的目标函数:

(9)

Liu等(2015)利用共轭梯度法不完全求解(9)式,采用改进的散射积分法将内部迭代中涉及到的Hessian-向量乘转化为2次核函数-向量乘运算;数值试验表明,随着内部迭代次数的增加,二阶方向逐渐趋近于真解,且内循环迭代10次即可获得较好的二阶方向的反演结果.公式(10)展示了采用改进的散射积分算法计算一阶方向p1的计算过程.

(10)

式中,pj1代表一阶方向p1在第j个模型参数上的值,Δui代表第i个波场残差,kij表示第j个模型参数扰动对第i个观测数据的影响.公式(10)右端求和的每一项,代表着一个炮检对所对应波场残差在核函数上的投影,通过累加求和,即得到一阶方向(Liu et al., 2015).

因此,采用改进的散射积分算法进行VTI介质弹性波多参数全波形反演的核心是推导各参数的敏感核函数.附录A根据Born近似与表示定理最终导出了Thomson参数表征下的各参数敏感核函数,表达式如(A18)-(A21)所示.

2 辐射模式分析

地震数据对不同参数的敏感性在特定的传播角度上有相似性,这使得不同参数在多参数反演过程中相互干涉,并相互耦合,增加了反演的非线性程度.为探究Thomsen参数表征下弹性波各参数之间的耦合效应,设计如下实验来分析各参数的辐射模式.散射波场可以看成是背景模型下的入射波场在遇到散射点之后的散射波场的叠加,故在相同激发条件下,将存在一个异常点的模型下得到的正演波场与背景模型下得到的正演波场相减,即可构建该单一散射点的散射波场.通过求解散射波场在各个方位辐射能量的大小及方向性来分析各参数之间的耦合效应,以帮助解释各参数更新方向的特点.

vp的辐射模式为例.如图 1(1)所示,为求解坐标轴(0, 0)位置处vp散射点在正上方波场入射情况下散射波场的能量辐射情况,将P波占优的Ricker子波震源置于均匀的各向异性介质中(vp= 3000 m·s-1vs=2000 m·s-1ε=0.15、δ=0.05、ρ=2000 kg·m-3),正演得到背景波场.将背景模型vp在(1.0 km, 1.0 km)位置的数值设为3300 m·s-1用来构成散射点,在相同的激发条件下激发,正演得到vp散射点产生的散射波场.通过两个波场相减即可获得在正上方入射条件下的由单一的vp散射点所产生的散射波场,将模型空间任一位置的所有时刻的波场平方累加到该模型空间位置,得到散射波场能量在模型空间上的分布情况.该过程即可数值求解得到vp在该激发条件及背景模型下的辐射模式.由于VTI介质的空间对称性,可通过变换激发点的位置,获得vp在45°及90°入射角情况下的辐射模式,并由此推测在任意角度入射情况下的辐射模式.在相同的激发方式和背景模型情况下,vpvsεδ四个参数在不同的入射方向下的辐射模式图如图 1所示.

图 1 Thomsen参数在不同入射角情况下的辐射模式 Fig. 1 Radiation patterns of Thomsen parameters with different incident angles

图 1中,图(1)展示了炮点0°入射时与散射点的相对位置,图(2)和图(3)分别为45°入射和90°入射示意图.图(a1)、(a2)、(a3)分别为vp在0°、45°、90°入射时的辐射模式;图(b1)、(b2)、(b3)分别为vs在0°、45°、90°入射时的辐射模式;图(c1)、(c2)、(c3)分别为ε在0°、45°、90°入射时的辐射模式;图(d1)、(d2)、(d3)分别为δ在0°、45°、90°入射时的辐射模式.

通过分析图 1可以总结各参数辐射模式特征.vpvsδ的辐射模式的主能量辐射方位都与入射角度有关,但ε的主能量辐射方位基本保持不变.vp的辐射模式为椭圆形向四周辐射能量;vs主要沿垂直于入射的方向辐射能量;ε主要沿水平方向,即垂直于TI介质的对称轴方向辐射能量,且该特征不随入射角度变化;δ在正上方入射时主要沿水平方向辐射能量,90°角入射时主要沿垂直方向辐射能量,斜角度45°入射时分别主要向垂直和水平方向辐射能量,且在三个入射角度情况下,45°、135°、225°、315°方向处都有较明显的能量辐射.

由Thomsen参数辐射模式的特征得到各参数之间的耦合特性:vpvs作为强参数,在相同条件下产生的散射波场最强,且向四周都有能量辐射;εδ作为弱参数,对强参数反演的影响相对微弱,但当90°入射时,强参数vs对于δ的耦合则会对δ的反演结果造成影响,当εδ两参数主能量辐射方位重合时(如0°入射时),则反演过程中,两参数间表现为较强的参数耦合;当两参数的能量辐射方位不同时(如90°入射时),则有希望将两个参数分别反演出来.

由以上分析可见,四个参数的辐射模式差异较大,相互耦合关系非常复杂.必须采用一定的反演策略或合适的参数化手段才有可能将多个参数逐步反演出来.而本文探讨的则是如何有效利用Hessian这一包含多参数耦合信息的矩阵,而不采用任何策略实现多参数的同时反演.

3 数值试验

为了测试弹性波方程全波形反演方法对各参数的反演能力,依据观测方式的不同,设置了两组数值试验,同时反演vpvsεδ四个参数.实验一采用四周观测系统,对比分析了初始模型条件下的最速下降方向和高斯-牛顿方向,并在此基础上完成了最速下降法和高斯-牛顿法多参数反演.实验二为地表观测系统情况下进行的理论模型试验.

3.1 四周观测数值试验

实验一真实模型如图 2所示.模型被离散为201×201个网格点,网格间距10 m×10 m.76炮均匀分布在四边,炮点深度100 m.单边放炮,其他三边接收.炮点的起始点位置为(100 m,100 m)处,炮间距100 m.每一炮对应108个检波点并在其余三边均匀分布,检波点间距50 m.模型vp的背景速度为3000 m·s-1,球状异常体速度为3300 m·s-1.模型vs的背景速度为2000 m·s-1,球状异常体速度为2200 m·s-1.ε模型背景值为0.15,球状异常值为0.2.δ模型背景值为0.05,球状异常值为0.1.密度ρ=2000 kg·m-3且反演过程中保持不变.vpvsεδ球状异常体分别位于模型的左上方、右下方、左下方和右上方,直径为200 m.

图 2 (a) 为vp真实模型,(b)为vs真实模型,(c)为ε真实模型,(d)为δ真实模型 Fig. 2 The true models of (a) vp, (b) vs, (c) ε, (d) δ
3.1.1 更新方向对比

为了对比最速下降方向和高斯牛顿方向,在相同的初始模型和观测系统条件下,使用一个频率组(2.0、2.5、3.0、3.5、4.0、4.5、5.0、5.5、6.0、6.5、7.0、7.5、8.0、8.5、9.0、9.5、10.0、10.5、11.0、11.5、12.0、12.5、13.0、13.5、14.0、14.5、15.0),分别求解以真实背景为初始模型的最速下降方向和高斯牛顿方向,并进行归一化对比.求取高斯-牛顿方向的内部循环迭代次数为10次.vpvsεδ更新方向的对比如图 3所示.

图 3 四周观测系统下四参数联合反演各参数归一化的梯度图(a)、(c)、(e)、(g)分别为vpvsεδ的最速下降方向,(b)、(d)、(f)、(h)分别为vpvsεδ的高斯-牛顿方向. Fig. 3 Normalized gradients of different parameters for four-parameter simultaneous inversion with the acquisition surrounding the edges (a, c, e, g) are SD descent directions, (b, d, f, h) are GN descent directions.

图 3中,(a)、(c)、(e)、(g)分别为vpvsεδ的最速下降方向,(b)、(d)、(f)、(h)分别为vpvsεδ的高斯-牛顿方向.(a)与(b)对比得,εδvp之间的参数耦合较弱,vp在(b)上与vs的参数耦合弱于(a);(c)与(d)对比得,εδvs之间的参数耦合较弱,vs在(d)上与vp的参数耦合弱于(c);(e)与(f)对比得,vpvsε之间的参数耦合在梯度上表现非常明显,且耦合方向水平,与ε的辐射模式一致,(f)上εvpvs之间的参数耦合显著减弱;(g)与(h)对比得,vpvsδ之间的参数耦合在δ更新方向上表现同样非常明显,由于δ能量辐射方向随入射角变化,vpvsδ之间的参数耦合在δ的更新方向上表现出各自的辐射方向特征.与(g)相比,(h)上δvpvs之间的参数耦合显著减弱.相比于最速下降方向,高斯牛顿方向上多参数之间耦合效应明显减弱.因此,在迭代更新的过程中,FWI能通过利用Hessian矩阵的信息有效提升多参数联合反演的精度.

3.1.2 多参数联合反演结果

分别使用最速下降法和高斯-牛顿法进行反演.最终的反演结果对比如图 4-图 7所示.

图 4 四周观测系统下的四参数联合反演vp结果对比(a)用最速下降方向的结果; (b)用高斯-牛顿方向的结果. Fig. 4 Inversion results of vp for four-parameter simultaneous inversion with the acquisition surrounding the edges (a) vp with SD; (b) vp with GN.
图 5 四周观测系统下的四参数联合反演vs结果对比(a)用最速下降方向的结果; (b)用高斯-牛顿方向的结果. Fig. 5 Inversion results of vs for four-parameter simultaneous inversion with the acquisition surrounding the edges (a) vs with SD; (b) vs with GN.
图 6 四周观测系统下的四参数联合反演ε结果对比(a)用最速下降方向的结果; (b)用高斯-牛顿方向的结果. Fig. 6 Inversion results of ε for four-parameter simultaneous inversion with the acquisition surrounding the edges (a) ε with SD; (b) ε with GN.
图 7 四周观测系统下的四参数联合反演δ结果对比(a)用最速下降方向的结果; (b)用高斯-牛顿方向的结果. Fig. 7 Inversion results of δ for four-parameter simultaneous inversion with the acquisition surrounding the edges (a) δ with SD; (b) δ with GN.

图 4(a, b)分别为使用最速下降方向(SD)和高斯-牛顿方向(GN)的vp反演结果.可以看出,两个方向的反演结果相似.高斯-牛顿方向的反演结果在vp异常的位置更接近于真实值,且vs的耦合效应削弱.

图 5(a, b)分别为使用最速下降方向(SD)和高斯-牛顿方向(GN)的vs反演结果.可以看出,最速下降方向和高斯-牛顿方向的反演结果相似,使用高斯-牛顿方向的反演结果在vs异常的位置更接近于真实值,且vp的耦合效应较小.

图 6(a, b)分别为使用最速下降方向(SD)和高斯-牛顿方向(GN)的ε反演结果.可以看出,高斯-牛顿方向的反演结果明显优于最速下降方向的反演结果,(b)在ε异常的位置更接近于真实值,且异常的形态刻画也更接近于真实.vpvs在(a)上表现出对ε明显的耦合效应,而在(b)上的参数耦合非常微弱.

图 7(a, b)分别为使用最速下降方向(SD)和高斯-牛顿方向(GN)的δ反演结果.可以看出,高斯-牛顿方向的反演结果同样明显优于最速下降方向的反演结果,(b)在δ异常的位置更接近于真实值.vpvs在(a)上表现出对δ明显的参数耦合效应,且在(a)上δ的异常产生的位置也错误,而在(b)上的参数耦合则非常微弱,能够更为真实的刻画出异常的位置和值.

通过对图 4-图 7的分析可得,最速下降法和高斯-牛顿法反演过程中,强参数vpvs都能获得很好的反演结果,其中使用高斯-牛顿方向的反演结果更佳;利用最速下降法无法获得良好的εδ反演结果,但利用Hessian矩阵的高斯-牛顿法反演得到的弱参数εδ的反演结果明显优于最速下降法.

为了验证该方法在强各向异性和高异常介质下的应用效果,我们在相同的实验参数下进行实验,只将真实模型的背景值和异常值进行了修改,模型vp的背景速度为3000 m·s-1,球状异常体速度为3900 m·s-1.模型vs的背景速度为2000 m·s-1,球状异常体速度为2600 m·s-1.ε模型背景值为0.3,球状异常值为0.4.δ模型背景值为0.18,球状异常值为0.28.密度ρ=2000 kg·m-3且反演过程中保持不变.vpvsεδ球状异常体分别位于模型的左上方、右下方、左下方和右上方,直径为200 m.仍然以真实背景作为初始模型进行反演,结果如图 8所示.

图 8 强各向异性下GN方向的反演结果(a) vp的反演结果; (b) vs的反演结果; (c) ε的反演结果; (d) δ的反演结果. Fig. 8 Inversion results of (a) vp, (b) vs, (c) ε and (d) δ with GN in strong anisotropic media

图 8可以看出,在强各向异性下以及高异常的情况下,本文利用Hessian矩阵的高斯-牛顿法反演方法仍然有良好的反演效果.

3.2 地表观测数值试验

实验二真实模型如图 9a图 10a图 11a图 12a所示,模型被离散为501×160个网格点,网格间距10 m×10 m.84炮均匀分布在地表,炮点深度50 m.地表激发,地表接收.炮点的起始点位置为(10 m,50 m)处,炮间距60 m.每一炮对应84个检波点均匀分布在地表,检波点间距60 m.模型为Overthrust模型,密度ρ=2000 kg·m-3且反演过程中固定不变.反演使用的频率范围为2~15 Hz.1个频率组(2.0、2.5、3.0、3.5、4.0、4.5、5.0、5.5、6.0、6.5、7.0、7.5、8.0、8.5、9.0、9.5、10.0、10.5、11.0、11.5、12.0、12.5、13.0、13.5、14.0、14.5、15.0),频率间隔0.5 Hz.使用水平与垂直方向上相关长度为21个网格的高斯滤波器对真实模型进行平滑,作为vpvsεδ的初始模型.分别使用最速下降法和高斯-牛顿法进行四参数同时反演.最速下降方向的迭代次数为84次,高斯-牛顿方向外部循环迭代53次,内部循环迭代次数为10次.反演结果如图 9-图 12所示,水平位置x=2 km处的抽线对比如图 13所示.

图 9 地表观测系统下的四参数联合反演vp结果对比(a)真实模型; (b)初始模型; (c)最速下降方向反演结果; (d)高斯-牛顿方向反演结果. Fig. 9 Inversion results of vp for four-parameter simultaneous inversion with the acquisition at the surface (a) The true model; (b) The starting model; (c) The inversion result with SD; (d) The inversion result with GN.
图 10 地表观测系统下的四参数联合反演vs结果对比(a)真实模型; (b)初始模型; (c)最速下降方向反演结果; (d)高斯-牛顿方向反演结果. Fig. 10 Inversion results of vs for four-parameter simultaneous inversion with the acquisition at the surface (a) The true model; (b) The starting model; (c) The inversion result with SD; (d) The inversion result with GN.
图 11 地表观测系统下的四参数联合反演ε结果对比(a)真实模型; (b)初始模型; (c)最速下降方向反演结果; (d)高斯-牛顿方向反演结果. Fig. 11 Inversion results of ε for four-parameter simultaneous inversion with the acquisition at the surface (a) The true model; (b) The starting model; (c) The inversion result with SD; (d) The inversion result with GN.
图 12 地表观测系统下的四参数联合反演δ结果对比(a)真实模型; (b)初始模型; (c)最速下降方向反演结果; (d)高斯-牛顿方向反演结果. Fig. 12 Inversion results of δ for four-parameter simultaneous inversion with the acquisition at the surface (a) The true model; (b) The starting model; (c) The inversion result with SD; (d) The inversion result with GN.
图 13 地表观测系统下的四参数联合反演结果抽线对比 Fig. 13 Vertical profiles for four-parameter simultaneous inversion with the acquisition at the surface

图 9图 10中,(c)、(d)分别为vpvs最速下降方向(SD)和高斯-牛顿方向(GN)的反演结果.可以看出,最速下降方向和高斯-牛顿方向的反演结果相似,高斯-牛顿方向的反演结果在异常的位置更接近于真实值.

图 11图 12中,(c)、(d)分别为εδ最速下降方向(SD)和高斯-牛顿方向(GN)的反演结果.可以看出,与最速下降方向(SD)相比,高斯-牛顿方向的ε反演结果在异常值的位置更接近于真实值.与以往的工作(刘玉柱等,2015Djebbi and Alkhalifah, 2017)相比,ε的反演精度得到明显的改善.但从图 13可以看出,δ作为最弱的参数,在复杂地下构造和仅有地表观测系统情况下难以反演得到真实值,借助于反演策略(如分步反演策略)或合适的参数化手段也许可以得到更好的结果.

4 结论

基于改进的散射积分算法的VTI介质弹性波全波形反演,在无需存储完整的Fréchet核函数和Hessian矩阵的情况下实现了最速下降方向和高斯-牛顿方向的灵活计算.通过分析Thomsen参数化情况下各参数的辐射模式及其随入射角的变化,说明了多参数耦合特征的复杂性,进而也说明了反演中有效利用Hessian矩阵的必要性.利用近似Hessian算子对梯度进行校正,有效降低更新方向上各参数之间的耦合效应,从而提升弹性波多参数联合反演的精度.数值试验表明,强参数vpvs在使用高斯-牛顿方向反演时,反演结果更接近真实值.在四周观测系统情况下,近似Hessian的利用能有效地减弱强、弱参数的耦合对反演结果的影响,进而εδ都能被很好地反演出来.在地表观测系统情况下,近似Hessian虽能降低反演中各参数之间的耦合作用,ε的反演精度得到了改善,但由于δ参数对数据的影响太小,仍不能被很好地反演出来,反演策略与参数化的应用仍是多参数反演中必要的辅助措施.

附录A 弹性波方程Thomsen参数核函数推导

二维VTI介质弹性波方程为:

(A1)

其中,uxuz为背景波场,c11c13c33c44ρ为背景模型参数.令 Δc11,且同样满足传播方程(A1),则:

(A2)

其中,Δc11为模型参数的微小扰动,Δux、Δuz是由模型扰动所引起的散射波场. 分别为x方向和z方向背景波场和散射波场之和,为背景模型与模型扰动之和.因为uxuzc11c13c33c44ρ同样满足传播方程(A1),化简方程(A2)得:

(A3)

根据散射积分原理,设GxxGzx为满足方程

(A4)

的解,其中(x′, z′)表示脉冲震源位置,GxxGzx分别表示在震源位置x方向上的集中力脉冲震源所对应的方程(A4)的x分量和z分量解.

同理,设GxzGzz为满足方程

(A5)

的解,GxzGzz分别表示在震源位置z方向上的集中力脉冲震源所对应的弹性波方程的x分量和z分量解.

根据表示定理,方程(A3)的解可用格林函数表示为:

(A6)

整理可获得由模型参数c11的微小扰动Δc11与其对应的扰动波场Δux、Δuz之间的关系式:

(A7)

将恒等式:

(A8)

代入方程(A7)中,其中的积分项满足:

(A9)

其中,xmaxxmin分别表示模型空间x方向的上边界和下边界.由于正演过程中模型边界位置存在吸收边界,反演过程中认为边界外的参数扰动Δc11为零,即方程(A9)中右端的积分项为零.由此可得:

(A10)

核函数反映了模型参数扰动与波场扰动之间的关系.在关系式(A10)中施加Born近似,用背景场uxuz来替代总场xz,得到模型参数c11的敏感核函数表达式:

(A11)

其中,kc11xkc11z分别表示c11x方向和z方向的敏感核函数.根据c11的核函数推导过程,同理可得c13c33c44ρ的敏感核函数表达式:

(A12)

(A13)

(A14)

(A15)

(A16)

利用Thomsen参数和刚度系数矩阵之间的关系式,通过链式法则可以得到Thomsen参数表征下任意模型参数m的核函数表达式:

(A17)

其中,kmi表示模型参数mi方向上的核函数(i=x, z).

将Thomsen参数与刚度系数之间的关系式代入公式(A17)中,并根据公式(A12)-(A16)可得vpvsεδ的核函数表达式:

(A18)

(A19)

(A20)

(A21)

References
Barnes C, Charara M, Tsuchiya T. 2008. Feasibility study for an anisotropic full waveform inversion of cross-well seismic data. Geophysical Prospecting, 56(6): 897-906. DOI:10.1111/j.1365-2478.2008.00702.x
Choi Y, Shin C, Min D J. 2007. Frequency domain elastic full waveform inversion using the new pseudo-Hessian matrix: elastic Marmousi-2 synthetic test.//2007 SEG Annual Meeting. San Antonio, Texas: SEG.
Dahlen F A, Hung S H, Nolet G. 2000. Fréchet kernels for finite-frequency traveltimes-I. Theory. Geophysical Journal International, 141(1): 157-174. DOI:10.1046/j.1365-246X.2000.00070.x
Djebbi R, Alkhalifah T. 2017. Frequency domain multi-parameter full waveform inversion for acoustic VTI media.//79th EAGE Conference and Exhibition. EAGE.
Fichtner A, Trampert J. 2011. Hessian kernels of seismic data functionals based upon adjoint techniques. Geophysical Journal International, 185(2): 775-798. DOI:10.1111/j.1365-246X.2011.04966.x
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. DOI:10.1111/j.1365-246X.2010.04767.x
Liu Y Z, Dong L G, Wang Y W, et al. 2009. Sensitivity kernels for seismic Fresnel volume tomography. Geophysics, 74(5): U35-U46. DOI:10.1190/1.3169600
Liu Y Z, Wang G Y, Yang J Z, et al. 2015. Multi-parameter full-waveform inversion for VTI media based on Born sensitivity kernels. Chinese Journal of Geophysics (in Chinese), 58(4): 1305-1316. DOI:10.6038/cjg20150418
Liu Y Z, Yang J Z, Chi B X, et al. 2015. An improved scattering-integral approach for frequency-domain full waveform inversion. Geophysical Journal International, 202(3): 1827-1842. DOI:10.1093/gji/ggv254
Luo J R, Wu R S, Gao F C. 2018. Time-domain full waveform inversion using instantaneous phase information with damping. Journal of Geophysics and Engineering, 15(3): 1032-1041. DOI:10.1088/1742-2140/aaa984
Métivier L, Brossier R, Virieux J, et al. 2012. Toward gauss-newton and exact newton optimization for full waveform inversion.//74th EAGE Conference and Exhibition. SPE, EAGE.
Operto S, Gholami Y, Prieux V, 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. DOI:10.1190/tle32091040.1
Pan W Y, Innanen K A, Margrave G F, et al. 2016. Estimation of elastic constants for HTI media using Gauss-Newton and full-Newton multiparameter full-waveform inversion. Geophysics, 81(5): R275-R291. DOI:10.1190/geo2015-0594.1
Plessix R É. 2009. Three-dimensional frequency-domain full-waveform inversion with an iterative solver. Geophysics, 74(6): WCC149-WCC157. DOI:10.1190/1.3211198
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. DOI:10.1190/1.1444597
Prieux V, Brossier R, Operto S, et al. 2013. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field. Part 1:imaging compressional wave speed, density and attenuation. Geophysical Journal International, 194(3): 1640-1664. DOI:10.1093/gji/ggt177
Sheen D H, Tuncay K, Baag C E, et al. 2006. Time domain Gauss-Newton seismic waveform inversion in elastic media. Geophysical Journal International, 167(3): 1373-1384. DOI:10.1111/j.1365-246X.2006.03162.x
Sirgue L, Barkved O I, van Gestel J P, et al. 2009.3D waveform inversion on Valhall wide-azimuth OBC.//71st EAGE Conference and Exhibition. SPE, EAGE.
Spetzler J, Snieder R. 2004. The Fresnel volume and transmitted waves. Geophysics, 69(3): 653-663. DOI:10.1190/1.1759451
Sun M A, Yang J Z, Dong L G, et al. 2017. Density reconstruction in multiparameter elastic full-waveform inversion. Journal of Geophysics and Engineering, 14(6): 1445-1462. DOI:10.1088/1742-2140/aa93b0
Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8): 1259-1266. DOI:10.1190/1.1441754
Woodward M J. 1992. Wave-equation tomography. Geophysics, 57(1): 15-26. DOI:10.1190/1.1443179
Yang J Z, Liu Y Z, Dong L G. 2014. A multi-parameter full waveform inversion strategy for acoustic media with variable density. Chinese Journal of Geophysics (in Chinese), 57(2): 628-643. DOI:10.6038/cjg20140226
Zhou B, Greenhalgh S. 2009. On the computation of the Fréchet derivatives for seismic waveform inversion in 3D general anisotropic, heterogeneous media. Geophysics, 74(5): WB153-WB163. DOI:10.1190/1.3123766
刘玉柱, 王光银, 杨积忠, 等. 2015. 基于Born敏感核函数的VTI介质多参数全波形反演. 地球物理学报, 58(4): 1305-1316. DOI:10.6038/cjg20150418
杨积忠, 刘玉柱, 董良国. 2014. 变密度声波方程多参数全波形反演策略. 地球物理学报, 57(2): 628-643. DOI:10.6038/cjg20140226