基于拟Hessian梯度预处理算子的勒夫波全波形反演研究

管建博，李 宇，殷裁云，杨 智，靳朝彬，赵 猛，杨 杭

基于拟Hessian梯度预处理算子的勒夫波全波形反演研究

管建博1,2，李 宇1,2，殷裁云3，杨 智1,2，靳朝彬1,2，赵 猛1,2，杨 杭1,2

(1. 长安大学 地质工程与测绘学院，陕西 西安 710054；2. 长安大学 西部矿产资源与地质工程教育部重点实验室，陕西 西安 710054；3. 华能煤炭技术研究有限公司，北京 101100)

构建近地表横波速度模型是煤田多分量地震资料处理的重要环节。相较于面波多道分析法，全波形反演在构建近地表横波模型中具有更高的分辨率。然而，在基于梯度的全波形反演中，由于地震记录频带有限、波场的非均匀覆盖以及双重散射等原因导致梯度算子不随深度的增加而缩放，模型深部参数得不到明显更新。目标函数的Hessian算子包含曲率信息，可清晰预测梯度算子中的焦散现象及双重散射产生的伪影，因此，逆Hessian算子则可作为反卷积算子实现对梯度的预处理，加强对模型深部的照明能力。然而Hessian算子具有巨大维度，对其显式计算十分困难。基于此，借鉴逆散射理论的思想，给出勒夫波全波形反演目标函数的拟Hessian算子的表达式，并提出一种梯度预处理的全波形反演方法。将该方法分别应用于断层模型、凹陷模型以及起伏界面模型的重构试验，反演结果表明：与传统的共轭梯度全波形反演方法相比，基于拟Hessian算子的预处理共轭梯度方法可加快收敛速度，提升成像质量。

全波形反演；勒夫波；Hessian算子；共轭梯度法；横波速度成像

静校正是复杂地区煤田地震勘探资料处理中的关键环节[1-2]，其中构建近地表横波速度模型较纵波速度模型更为复杂。利用地震数据中信噪比较高的“噪声”—面波建立横波速度模型是一种有效的方法。基于频散曲线反演的面波多道分析法(Multichannel Analysis of Surface Waves, MASW)已经成为当前浅地表横波速度成像最流行的方法之一[3-4]。孟小红等[5]指出该方法得到的是 P-SV 波横波静校正量的长波长趋势，无法预测横波静校正量的短波长分量，给出的成像结果也是通过插值形成的拟二维(2D)剖面，限制了该方法的横向分辨率；传统的偏移速度分析法[6]和走时层析反演[7-8]虽然能够直接实现横波速度2D成像，但只能得到宏观速度场，即速度的低频成分，难以得到高分辨率的速度模型。全波形反演(Full Waveform Inversion, FWI)[9-10]直接利用地震波形信息，可以得到高精度的地下速度模型，能够满足当前复杂构造成像方法对速度参数的要求。以往的FWI研究多集中于纵、横波反演[11-13]，而面波能量由于在近地表波场占据主导地位(垂直激发瑞雷波约70%，水平激发勒夫波约90%)使得发展面波全波形反演成为可能[14]；相比于瑞雷波，勒夫波由近自由地表多次反射SH波经干涉形成[15]，与纵波无关且对横波速度更为敏感，这决定了勒夫波FWI的框架更为简单且更为有效[16-17]。在经典的基于梯度的FWI方法中，由于双重散射、波场的非均匀覆盖以及子波的频带受限等原因导致由伴随状态法[18]求取的梯度算子存在奇异值且能量分布不均，模型深部的梯度振幅急剧减小，照明能力下降，模型参数得不到明显更新。O. Gauthier等[19]在计算各处梯度时，引入了与震源的距离参数，使远处的反演结果得到了改善，这种简单预处理对于复杂实际问题的效果一般；目标函数的二阶偏导–Hessian矩阵包含曲率信息，可清晰预测模型深部照明缺失现象，因此，Hessian矩阵的逆(逆Hessian算子)则可调节梯度算子比例并均衡其修正量以加强深部照明能力[13]。逆Hessian算子，或者称为梯度预处理算子，具有巨大维度，依靠当前的计算能力想要对其显式计算十分困难，通常采用估计或者构造的方式获得[13]。R. G. Pratt等[20]采用高斯–牛顿算法实现了频率–空间域的FWI，给出了Hessian矩阵的完整表达式并指出：在高频近似下逆Hessian矩阵是对角占优矩阵，对于包含个参数的反演问题，只需计算逆Hessian矩阵的个对角元素，这极大地简化了其的计算过程。C. Shin等[21]基于逆散射理论给出了叠前深度偏移中拟Hessian矩阵(用一对角占优矩阵近似代替Hessian矩阵)对角项的估计方法，但缺点是缺失几何扩散相关的因子。D. H. Sheen等[22]采用高斯–牛顿算法实现了时间域弹性波全波形反演，并给出了拟Hessian矩阵对角项的计算方法；还有一种思路是利用目标函数及其梯度来构建拟Hessian矩阵，其中以l-BFGS[23-24]最为常用。L. Métivier等[13]将FWI视为牛顿线性系统，采用二阶伴随状态法来估计更精确的对角项，但由于每次迭代都需要4个正演过程，大大增加了计算时间。为了进一步降低Hessian矩阵的计算成本并提高反演精度，修正的拟牛顿波形反演算法[25]结合梯度、模型信息及目标函数信息估计了一个投影Hessian矩阵并成功应用于浅地表的模型重构。Li Jing等[26]在勒夫波波动方程频散反演中通过对观测数据进行加权形成一个新的伴随源来对梯度进行修正。Liu Zhaolun等[27]利用该方法实现了瑞雷波的三维成像。Yan Yingwei等[28-30]基于叠前深度偏移理论中求取逆Hessian矩阵的方法求取预处理算子，并成功将其应用到勒夫波实际数据全波形反演中，取得较好效果。

笔者基于逆散射理论中求取Hessian矩阵的方法求取勒夫波FWI目标函数的拟Hessian算子实现对梯度的预处理，以均衡整体梯度振幅加强算法对深部介质的照明能力，提升地层深部反演精度。为解决初始模型拙劣以及波场照明不足导致的反演精度丢失问题，采用伪井地联合反演，使地震记录中携带更多的深部介质信息，从本质上提升成像质量。

1 方法原理

1.1 正演问题

在二维各向同性介质中，假设质点在-平面上只有方向的非零位移(为水平距离，为深度)，结合虎克定律，得到伪保守形式的勒夫波弹性动力

学方程[31]:

式(1)的一般形式为：

式中：为模型参数；为式(1)的偏微分算子；为地震波场；为加载的震源项。

1.2 局部优化算法

FWI计算量巨大，使得全局优化和混合优化算法难以在FWI领域得到应用。因此，只有少数局部优化算法经过了测试[34]，即从一个初始模型出发逐步迭代寻优，迭代序列可以表示为：

对于基于梯度的算法，共轭梯度法(Conjugate Gradient method，CG)在FWI中得到了广泛的应用。其模型扰动为：

对于预处理算法，其核心是求解目标函数的逆Hessian矩阵，新的模型扰动为:

式中：上标T代表转置运算。对于预处理算法中的计算[35]，选用下式：

为了加快目标函数的收敛速度，令标量参数为非负值[36]：

1.3 目标函数梯度及预处理算子

全波形反演的本质思想就是将观测的地震数据与预测的地震数据进行匹配，使某种衡量二者误差关系的目标函数值达到最小，从而得到最佳的模型数据。目标函数可以定义为：

求取目标函数的梯度是局部优化算法的核心。伴随状态法[18]指出，梯度是正传波场的时间偏导与检波器位置处反传波场残差获得的反传波场的零延迟互相关，这在C. Castellanos等[37]的研究中有详细的推导说明。因此，结合链式法则[38]，可以得到勒夫波FWI中关于横波速度的梯度算子，其差分格式为：

由式(10)求取的梯度由于地震数据的频带受限、双重散射以及波场非均匀覆盖等原因会随着深度的加深振幅急剧减小，导致对深部地层的照明缺失，介质参数得不到明显更新，而目标函数的二阶导数(Hessian矩阵)由于包含曲率信息，可清晰预测这种现象。逆Hessian矩阵，也可称其为预处理算子，便可调节梯度算子比例并均衡其修正量，消除带限效应，以加强深部地层的照明能力，提高反演精度。由于其维度过于庞大，在当前的计算能力下显式计算这个矩阵几乎是不可能的，但是可以通过互易定理获得拟Hessian矩阵。本文参考叠前深度偏移中的拟Hessian矩阵的求解方法来求解时间域勒夫波全波形反演中的拟Hessian矩阵。这种可以压制由多次散射在梯度向量中产生的假象的求解 Hessian 矩阵的理论，被称为逆散射理论[21]。根据逆散射理论，拟Hessian矩阵的表达式为：

由此式即可得到勒夫波FWI中关于横波速度的拟Hessian矩阵，其差分格式为：

关于应力时间的导数在求取梯度时已经计算，求取拟Hessian矩阵时只需要相乘再叠加即可。由于逆Hessian矩阵可看作是对角占优矩阵，故只需计算出其对角项。同时，为了避免对角项中由于波场的几何扩散等原因出现近零值导致反演过程不稳定，故对其加上一个调节因子，第次迭代时的预处理算子表达式如下：

将式(14)代入式(5)，即形成了预处理的共轭梯度算法(Pre-conditioned Conjugate Gradient method，PCG)。

2 模型重构测试

为了验证上述PCG比CG更具优越性，我们将这2种FWI算法同时应用于断层模型、凹陷模型以及起伏界面模型进行反演试算。另外，在起伏界面模型试 算中加入井地联合反演，以证明其在提高反演精度、提升成像质量方面的有效性。在反演过程中，认为模型密度及地震子波已知。为了定量评价反演结果精度，对比2种算法的优劣，引入模型的对比误差(RMSE)，其表达式为：

2.1 断层模型

图1 断层模型重构测试中的真实模型、初始模型及反演结果

图1c、图1d分别为PCG与CG算法的反演结果，可以看出2种算法都较好地重建了地层异常，证明了程序的正确性。CG算法可以精确重建模型浅部形态且给出大致的异常构造及其边界，但不能重构模型深部形态。这是因为随着深度的加深，梯度算子振幅逐渐变小，模型深部参数没有得到明显更新，可从第22道和第35道处的1维(1-D)横波速度深度剖面(图3a、图3b)得到更为直观的反映，深层速度偏低。而PCG算法则可以更精确地重建模型，横波速度的误差终值仅为0.022 9。图3c为2种算法的归一化目标函数值变化曲线的对比图，PCG的收敛速度明显快于CG。在前5次迭代，PCG的函数值迅速收敛至极低水平，说明反演模型已更新到真实模型附近，随后逐渐逼近真实模型，迭代17次即满足精度要求，且目标函数终值也更小，仅为0.011 6。而CG前7次迭代的函数值收敛也较快，这是因为模型浅部的更新对函数值的贡献很大，之后由于模型深部的梯度振幅过小，参数更新不足，函数值收敛速度变慢，误差终值为0.047 2。由图3d两种算法反演结果的对比误差可以看出，在整个迭代过程中，PCG的对比误差一直低于CG，从0.080 1下降到0.022 9，而CG算法仅能下降到0.047 2(表1)，反演精度远远高于CG。图3e给出的波形拟合图说明2种方法的反演波形均能很好的拟合于真实波形。

图3 断层模型重构测试的评价曲线

表1 断层模型重构测试中PCG与CG算法的性能对比评价

在断层模型反演过程中我们发现，由于FWI中固有的参数过度估计与参数估计不足的缺陷，导致最终得到的反演结果中，地层深部速度上界高于500 m/s，地表速度下界低于300 m/s(图3a、图3b)。这是由于迭代过程中预处理算子对梯度算子过度调节，出现了照明过度的现象。因此，需要在迭代寻优过程中加入模型更新的约束控制条件来解决该问题。

2.2 凹陷模型

为进一步对比分析PCG比CG算法具有优越性，继续对凹陷模型进行反演试算。模型垂直与水平网格数目为81×101，垂直和水平网格步长均为0.5 m，即模型大小为40 m×50 m，如图4a所示。图4a模型分为4层：第一层厚度5 m，横波速度为200 m/s；第二层厚度10 m，横波速度为400 m/s，但在横向32～40 m间存在一个高度5 m的高速凸起；第三层厚度10～17.5 m，横波速度为400 m/s，在横向8～18 m间存在一个高度5 m的低速下陷。第四层厚度7.5～15 m，速度为400 m/s，密度为恒定值1 800 kg/m3。观测系统与断层模型基本一致，不同之处是激发点距改为5 m，保持激发点数目11不变(图4a中的白色箭头所示)，水平分量接收点为101个(图4a中的黄色圆点，已抽稀至51个显示)。波场正演采用断层模型参数，正演地震记录(第6炮)如图5所示。将初始模型划分为40层，同样参考地震记录中勒夫波的基阶频散曲线高频端、低频端相速度，设置初始模型速度范围为200～500 m/s，层速度增量为7.5 m/s，如图4b所示。

在该模型的试算中，仍存在参数过度估计与估计不足的现象(图6a、图6b)。由图6e的波形拟合图可以看出，初始波形与真实波形相差较大，说明初始模型比较拙劣，这是导致CG算法的反演波形没能与真实波形拟合的因素之一，同时也是图6c中CG算法的目标函数值终值较大，仅能收敛至0.052 1(表2)的原因之一。另一原因则是相较于断层模型，该模型深度变深，梯度振幅受限更为严重，导致反演精度进一步降低。但PCG算法的反演波形与真实波形完全拟合，同时目标函数值也近似收敛于0(表2)。图6d中CG算法的横波速度对比误差在前5次迭代中与PCG算法相当，这是因为模型浅部(10 m以上)的更新贡献较大。之后PCG算法由于对模型深部照明能力更强，使得深部介质参数得到更新，对比误差下降较快。图4c、图4d分别给出了PCG与CG算法的反演结果，可以清晰地看出，相较于断层模型，2种方法的成像质量差距明显加大。PCG算法较精确地重建了地层异常及特殊构造，清晰地刻画出了异常边界，横波速度对比误差终值由0.102 8下降至0.030 6(表2)；而CG算法只能较准确地给出模型浅层的形态，模糊地反演出模型中部存在高速凸起、低速凹陷，而模型底部的断层则完全没有显示出来，横波速度的对比误差下降量很小，终值为0.088 7(表2)。这一定程度上展现了预处理算子在提升深部照明强度、提高反演精度方面的作用，说明PCG算法在FWI中的有效性。

图4 凹陷模型重构测试中的真实模型、初始模型及反演结果

图5 凹陷模型重构测试的第6炮地震记录

2.3 起伏界面模型

最后对界面起伏模型进行反演试算，如图7a所示，真实模型由界面起伏的四层速度递增的地层组成，横波速度依次为200、300、400、500 m/s，差异已足够显著。模型尺寸、观测系统、波场正演参数与凹陷构造模型完全相同。图8为获得的第6炮

地震记录。将初始模型划分为20层，同样参考地震记录中勒夫波的基阶频散曲线高频端、低频端相速度，设置初始模型速度范围为200～500 m/s，层速度增量为15 m/s，如图7b所示。

在该模型的试算中，仍存在参数过度估计与估计不足的现象(图9a、图9b)。图9c为2种算法的目标函数值变化曲线，CG算法由于对深部介质参数的更新不足，在第21次迭代便已满足寻优终止条件，目标函数终值过大，为0.117 3(表3)。PCG算法的收敛速度快于CG算法，且目标函数终值也更小，为0.011 5(表3)。图7c、图7d分别为PCG与CG算法的反演结果，可以看出CG算法仅仅能重建模型浅部形态，不能给出模型深部的地层边界，原因已在前2个模型试算中进行了详细分析。而PCG算法虽然由于预处理算子的补偿作用清楚地给出了各地层边界，反演精度有了一定程度的提高，但是由于“波峰”处的地层(图7c中箭头所指的部分)波场照明不足，地震波无法穿透，使地震记录中缺失这部分地层的波场信息，导致参数更新不足，由第22道处的1-D速度深度剖面(图9a)可以得到更为直观的反映。图9e给出了波形拟合结果，PCG算法的反演波形与真实波形已完全拟合，但横波速度对比误差变化曲线(图9d)表明2种算法的反演精度均有限，可能是由于短周期的周波跳跃导致的。另外，PCG算法在第6次迭代后，对比误差不降反升，主要原因是迭代后期对中深部地层照明过度，参数过度估计(图9b，第90道处的1-D速度深度剖面)。

表2 凹陷模型重构测试中PCG与CG算法性能对比评价

图7 起伏界面模型重构测试中的真实模型、初始模型及反演结果

Fig.7 The real model, initial model and inversion results in the undulating interface model reconstruction test

图8 起伏界面模型重构测试的第6炮地震记录

表3 界面起伏模型重构测试中PCG与CG算法性能对比评价

图9 界面起伏模型重构测试的评价曲线

2.4 伪井地联合反演

为了进一步提高FWI成像质量，本文开展了伪井地联合(Pseudo Borehole-Ground joint，PBGJ)的PCG算法测试。通过在地面和井中联合激发震源，加强模型深部的波场照明。由于地震数据包含了模型深部介质的丰富信息，反演精度将得到提高。将PBGJ-PCG算法应用于起伏界面模型测试中，除了地表激发11激发点外，在垂直方向上距离第1个激发点25 m处再增加1个激发点(图7a所示白色箭头)，其坐标为(0 m, 25 m)。

图10c给出了PBGJ-PCG算法反演结果的评价曲线，由于井地联合中地震记录包含的地层信息更为丰富，在前期的迭代寻优中便对模型整体进行了更新，而模型浅部的更新对目标函数值贡献很大，导致只能对模型浅层进行更新的PCG算法的目标函数收敛速度比PBGJ-PCG算法更快，但目标函数终止值更大(表4)，对比误差也一直更高(图10d)，但反演波形均完全拟合于真实波形(图10e)，推测是由于周波跳跃导致的。由图10a给出的反演结果及图10b所示的第22道处的1-D速度剖面可知，PBGJ-PCG算法提高了“波峰”处地层的反演精度，介质参数得到了明显更新，较精确地刻画了地层边界，重构了地层速度。对比误差由0.089 3下降至0.033 4，而PCG算法的对比误差只下降至0.078 5(表4)，展现了PBGJ-PCG算法在提高反演精度方面的有效性。

3 结论

a.拟Hessian梯度预处理算子可对梯度进行能量均衡，提高对深部介质的照明能力，对FWI具有重要意义。同时，由此衍生出的PCG算法相较于经典的CG算法，可加快收敛速度，即使初始模型比较拙劣，仍能比较准确地给出地层的异常结构及边界，提高深部介质的反演精度，提升成像质量。

b. PBGJ-PCG虽然无法解决时间域FWI固有的周波跳跃问题，但该方式获得的地震记录中携带了更多的深部介质信息，从本质上提高了对深部介质的照明能力，可进一步提升PCG算法的成像质量。

c. PCG算法与井地联合反演方式虽然在本文中仅应用于勒夫波的全波形反演中，但在瑞利波及体波的FWI中同样适用。

表4 PBGJ-PCG与PCG算法性能对比评价

d. PCG算法可能存在对梯度算子的过度调节，导致参数过度估计的现象出现。因此，迭代寻优过程中需要加入模型更新的约束控制条件，这需要进一步深入研究。

[1] 张胤彬，潘冬明，胡明顺，等. 复杂地表初至层析地震静校正技术及其应用[J]. 煤田地质与勘探，2012，40(4)：66–70. ZHANG Yinbin，PAN Dongming，HU Mingshun，et al. Research and application of first arrival tomographic static technique in complex surface area[J]. Coal Geology & Exploration，2012，40(4)：66–70.

[2] 周俊杰，王雨，侯玮. 黄土塬地区煤田三维地震综合处理技术[J]. 地球物理学进展，2016，31(5)：2299–2305. ZHOU Junjie，WANG Yu，HOU Wei. 3D seismic comprehensive processing technology of coalfield in loess tableland[J]. Progress in Geophysics，2016，31(5)：2299–2305.

[3] PARK C B，MILLER R D，XIA Jianghai. Multichannel analysis of surface waves[J]. Leading Edge，1999，18(3)：800–808.

[4] XIA Jianghai，XU Yixian，LUO Yinhe，et al. Advantages of using multichannel analysis of Love waves(MALW) to estimate near-surface shear-wave velocity[J]. Surveys in Geophysics，2012，33(5)：841–860.

[5] 孟小红，郭良辉. 利用地震瑞利波速度反演求取 P–SV波横波静校正量[J]. 石油地球物理勘探，2007，42(4)：448–453. MENG Xiaohong，GUO Lianghui. Using velocity inversion of seismic Rayleigh wave to compute S-wave statics of P-SV wave[J]. Oil Geophysical Prospecting，2007，42(4)：448–453.

[6] BIONDI B，SYMES W W. Angle-domain common-image gathers for migration velocity analysis by wavefield-continuation imaging[J]. Geophysics，2004，69(5)：1283–1298.

[7] ZHU Xianhuai，MCMECHAN G A. Estimation of a two-dimensional seismic compressional-wave velocity distribution by iterative tomographic imaging[J]. International Journal of Imaging Systems and Technology，1989，l(1)：13–17.

[8] ZHANG Jianzhong，SHI Taikun，ZHAO Yasheng，et al. Static corrections in mountainous areas using Fresnel-wavepath tomography[J]. Journal of Applied Geophysics，2014，111：242–249.

[9] TARANTOLA A. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics，1984，49(8)：1259–1266.

[10] TARANTOLA A. A strategy for nonlinear elastic inversion of seismic reflection data[J]. Geophysics，1986，51(10)：1893–1903.

[11] MÉTIVIER L，BROSSIER R，VIRIEUX J，et al. The truncated Newton method for full waveform inversion[C]//Las Vegas：SEG Annual Meeting，2012.

[12] MÉTIVIER L，BROSSIER R，VIRIEUX J，et al. Full waveform inversion and the truncated Newton method[J]. SIAM Journal on Scientific Computing，2013，35(2)：B401–B437.

[13] MÉTIVIER L，BRETAUDEAU F，BROSSIER R，et al. Full waveform inversion and the truncated Newton：Quantitative imaging of complex subsurface structures[J]. Geophysical Prospecting，2014，62(6)：1353–1375.

[14] 夏江海. 高频面波方法[M]. 武汉：中国地质大学出版社，2015. XIA Jianghai. High frequency surface wave method[M]. Wuhan：China University of Geosciences Press，2015.

[15] ESLICK R，TSOFLIAS G，STEEPLES D. Field investigation of Love waves in near-surface seismology[J]. Geophysics，2008，73(3)：1–6.

[16] PAN Yudi，XIA Jianghai，XU Yixian，et al. Love-wave waveform inversion in time domain for shallow shear-wave velocity[J]. Geophysics，2015，81(1)：1–14.

[17] DOKTER E，KOHN D，WILKEN D，et al. Full waveform inversion of SH- and Love-wave data in near-surface prospecting[J]. Geophysical Prospecting，2017，65(Sup. 1)：216–236.

[18] 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.

[19] GAUTHIER O，VIRIEUX J，TARANTOLA Albert. Two-dimensional nonlinear inversion of seismic waveforms：Numerical results[J]. Geophysics， 1986，51(7)：1387–1403.

[20] PRATT R G，SHIN C，HICKS G J. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J]. Geophysical Journal International，1998，133：341–362.

[21] SHIN C，JANG S，MIN D J. Improved amplitude preservation for prestack depth migration by inverse scattering theory[J]. Geophysical Prospecting，2001，49：592–606.

[22] SHEEN D H，TUNCAY K，BAAG C E，et al. Time domain Gauss：Newton seismic waveform inversion in elastic media[J]. Geophysical Journal International，2006，167(3)：1373–1384.

[23] BROSSIER R，OPERTO S，VIRIEUX J. Seismic imaging of complex on shore structures by 2D elastic frequency-domain full-waveform inversion[J]. Geophysics，2009，74：63–76.

[24] RAO Ying，WANG Yanghua，HAN Dehao. Seismic waveform tomography with simplified restarting scheme[J]. IEEE Geoscience and Remote Sensing Letters，2018，16(1)：135–139.

[25] 刘璐，刘洪，张衡，等. 基于修正拟牛顿公式的全波形反演[J]. 地球物理学报，2013，56(7)：2447–2451. LIU Lu，LIU Hong，ZHANG Heng，et al. Full waveform inversion based on modified quasi-Newton equation[J]. Chinese Journal of Geophysics，2013，56(7)：2447–2451.

[26] LI Jing，HANAFY S，LIU Zhaolun，et al. Wave equation dispersion inversion of love waves[J]. Geophysics，2019，84(5)：693–705.

[27] LIU Zhaolun，LI Jing，HANAFY S M，et al. 3D wave-equation dispersion inversion of Rayleigh waves[J]. Geophysics，2019，84(5)：673–691.

[28] YAN Yingwei，WANG Zhejiang，LI Jing，et al. A preconditioned technique for SH- and Love-wave full-waveform inversion in time domain and crosstalk analysis[J]. Journal of Geophysics and Engineering，2019，17(1)：160–174.

[29] YAN Yingwei，WANG Zhejiang，LI Jing，et al. Elastic SH- and Love-wave Full-Waveform Inversion for shallow shear wave velocity with a pre-conditioned technique[J]. Journal of Applied Geophysics，2020，173(2)：103947.

[30] 闫英伟. 浅地表高频面波成像技术研究[D]. 长春：吉林大学，2019. YAN Yingwei. Resear on high-frequency surface wave imaging technique for the shallow subsurface[D]. Changchun：Jilin University，2019.

[31] VIRIEUX J. SH wave propagation in heterogeneous media：Velocity stress finit-difference method[J]. Geophysics，1984，49：1933–1942.

[32] GRAVES R W. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences[J]. Bulletin of the Seismological Society of America，1996，86(4)：1091–1106.

[33] MEZA-FAJARDO K C，PAPAGEPRGIOS A S. A nonconvalutional，spllit-field，perfectly matched layer for wave propagation in isotropic and anisotropic elastic media：Stability analysis[J]. Bulletin of Seismological Society of America，2008，98(4)：1811–1836.

[34] RODI W，MACKIE R L. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion[J]. Geophysics，2001，66(1)：174–187.

[35] 刘聪，王者江，闫英伟. 基于伴随状态法二维时间域勒夫波全波形反演研究[J]. 地球物理学进展，2019，34(1)：136–143. LIU Cong，WANG Zhejiang，YAN Yingwei. Research on two-dimensional Love-wave full waveform inversion in time domain based on adjoint state method[J]. Progress in Geophysics，2019，34(1)：136–143.

[36] GILBERT J C，NOCEDAL J. Global convergence properties of conjugate gradient methods for optimization[J]. Siam Journal on Optimization，1992，2(1)：21–42.

[37] CASTELLANOS C，ETIENNE V，HU Guanghui，et al. Algorithmic and methodological developments towards full waveform inversion in 3D elastic media[C]//San Antonio：SEG Annual Meeting，2011.

[38] MORA P. Nonlinear two-dimensional elastic inversion of multioffset seismic data[J]. Geophysics，1987，52(9)：1211–1228.

Love wave full waveform inversion via Pseudo-Hessian gradient pre-conditioning operator

GUAN Jianbo1,2, LI Yu1,2, YIN Caiyun3, YANG Zhi1,2, JIN Chaobin1,2, ZHAO Meng1,2, YANG Hang1,2

(1. School of Geological Engineering and Geomatics, Chang’an University, Xi’an 710054, China; 2. Key Laboratory of Western China Mineral Resources and Geological Engineering, Chang’an University, Xi’an 710054, China; 3. Huaneng Coal Technology Research Co., Ltd., Beijing 101100, China)

The construction of near surface shear wave velocity is an important step in multi-component seismic data processing in coalfield. Compared with the multichannel analysis of surface wave, the full waveform inversion(FWI) has higher resolution in the construction of near surface shear wave velocity model. However, in the gradient-based FWI, the gradient operator is not scaled with increasing depth due to the narrow frequency band of seismic records, the non-uniform coverage of the wavefield, and the double scattering. The parameters of the deep model cannot be updated significantly. The Hessen operator of the objective function contains curvature information, which can clearly predict the defocusing phenomenon and the artifacts generated by double scattering in the gradient operator. The inverse Hessen operator can be used as a deconvolution operator to realize gradient pre-conditioning and enhance the illumination ability of the deep model. However, the explicit calculation of Hessian operator is very difficult because it has huge dimensions. Based on this, inverse scattering theory is referred to, the expression of the pseudo-Hessian operator of the objective function of full-waveform inversion is given, and a pre-conditioned gradient-based FWI method is developed. The proposed method was applied to the reconstruction tests of the fault model, subsidence model, and undulating interface model, respectively. The inversion results show that, compared with the classic conjugate gradient-based FWI, the pre-conditioned conjugate gradient method based on the pseudo-Hessian operator can accelerate the convergence rate and improve the inversion accuracy.

full waveform inversion; Love wave; Hessian operator; conjugate gradient method; shear-wave velocity imaging

P315.9

A

1001-1986(2021)04-0049-11

2021-04-02；

2021-04-21

国家重点研发计划课题(2018YFC0807803)；陕西省自然科学基础研究计划项目(2019JLM8)；长安大学中央高校基本科研业务费专项基金项目(300102260203)

管建博，1997年生，男，河北承德人，硕士研究生，研究方向为浅地表面波成像. E-mail：Jianbo_Guan @126.com

李宇，1983年生，男，湖北孝感人，博士，讲师，研究方向为浅地表地震勘探.E-mail：liyupa@chd.edu.cn

管建博，李宇，殷裁云，等. 基于拟Hessian梯度预处理算子的勒夫波全波形反演研究[J]. 煤田地质与勘探，2021，49(4)：49–59. doi: 10.3969/j.issn.1001-1986.2021.04.007

GUAN Jianbo，LI Yu，YIN Caiyun，et al. Love wave full waveform inversion via Pseudo-Hessian gradient pre- conditioning operator[J]. Coal Geology & Exploration，2021，49(4)：49–59. doi: 10.3969/j.issn.1001-1986.2021. 04.007

(责任编辑 聂爱兰)