管建博，李 宇，殷裁云，楊 智，靳朝彬，趙 猛，楊 杭

基于擬Hessian梯度預(yù)處理算子的勒夫波全波形反演研究

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

(1. 長安大學(xué) 地質(zhì)工程與測繪學(xué)院，陜西 西安 710054；2. 長安大學(xué) 西部礦產(chǎn)資源與地質(zhì)工程教育部重點實驗室，陜西 西安 710054；3. 華能煤炭技術(shù)研究有限公司，北京 101100)

構(gòu)建近地表橫波速度模型是煤田多分量地震資料處理的重要環(huán)節(jié)。相較于面波多道分析法，全波形反演在構(gòu)建近地表橫波模型中具有更高的分辨率。然而，在基于梯度的全波形反演中，由于地震記錄頻帶有限、波場的非均勻覆蓋以及雙重散射等原因?qū)е绿荻人阕硬浑S深度的增加而縮放，模型深部參數(shù)得不到明顯更新。目標(biāo)函數(shù)的Hessian算子包含曲率信息，可清晰預(yù)測梯度算子中的焦散現(xiàn)象及雙重散射產(chǎn)生的偽影，因此，逆Hessian算子則可作為反卷積算子實現(xiàn)對梯度的預(yù)處理，加強對模型深部的照明能力。然而Hessian算子具有巨大維度，對其顯式計算十分困難。基于此，借鑒逆散射理論的思想，給出勒夫波全波形反演目標(biāo)函數(shù)的擬Hessian算子的表達式，并提出一種梯度預(yù)處理的全波形反演方法。將該方法分別應(yīng)用于斷層模型、凹陷模型以及起伏界面模型的重構(gòu)試驗，反演結(jié)果表明：與傳統(tǒng)的共軛梯度全波形反演方法相比，基于擬Hessian算子的預(yù)處理共軛梯度方法可加快收斂速度，提升成像質(zhì)量。

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

靜校正是復(fù)雜地區(qū)煤田地震勘探資料處理中的關(guān)鍵環(huán)節(jié)[1-2]，其中構(gòu)建近地表橫波速度模型較縱波速度模型更為復(fù)雜。利用地震數(shù)據(jù)中信噪比較高的“噪聲”—面波建立橫波速度模型是一種有效的方法。基于頻散曲線反演的面波多道分析法(Multichannel Analysis of Surface Waves, MASW)已經(jīng)成為當(dāng)前淺地表橫波速度成像最流行的方法之一[3-4]。孟小紅等[5]指出該方法得到的是 P-SV 波橫波靜校正量的長波長趨勢，無法預(yù)測橫波靜校正量的短波長分量，給出的成像結(jié)果也是通過插值形成的擬二維(2D)剖面，限制了該方法的橫向分辨率；傳統(tǒng)的偏移速度分析法[6]和走時層析反演[7-8]雖然能夠直接實現(xiàn)橫波速度2D成像，但只能得到宏觀速度場，即速度的低頻成分，難以得到高分辨率的速度模型。全波形反演(Full Waveform Inversion, FWI)[9-10]直接利用地震波形信息，可以得到高精度的地下速度模型，能夠滿足當(dāng)前復(fù)雜構(gòu)造成像方法對速度參數(shù)的要求。以往的FWI研究多集中于縱、橫波反演[11-13]，而面波能量由于在近地表波場占據(jù)主導(dǎo)地位(垂直激發(fā)瑞雷波約70%，水平激發(fā)勒夫波約90%)使得發(fā)展面波全波形反演成為可能[14]；相比于瑞雷波，勒夫波由近自由地表多次反射SH波經(jīng)干涉形成[15]，與縱波無關(guān)且對橫波速度更為敏感，這決定了勒夫波FWI的框架更為簡單且更為有效[16-17]。在經(jīng)典的基于梯度的FWI方法中，由于雙重散射、波場的非均勻覆蓋以及子波的頻帶受限等原因?qū)е掠砂殡S狀態(tài)法[18]求取的梯度算子存在奇異值且能量分布不均，模型深部的梯度振幅急劇減小，照明能力下降，模型參數(shù)得不到明顯更新。O. Gauthier等[19]在計算各處梯度時，引入了與震源的距離參數(shù)，使遠處的反演結(jié)果得到了改善，這種簡單預(yù)處理對于復(fù)雜實際問題的效果一般；目標(biāo)函數(shù)的二階偏導(dǎo)–Hessian矩陣包含曲率信息，可清晰預(yù)測模型深部照明缺失現(xiàn)象，因此，Hessian矩陣的逆(逆Hessian算子)則可調(diào)節(jié)梯度算子比例并均衡其修正量以加強深部照明能力[13]。逆Hessian算子，或者稱為梯度預(yù)處理算子，具有巨大維度，依靠當(dāng)前的計算能力想要對其顯式計算十分困難，通常采用估計或者構(gòu)造的方式獲得[13]。R. G. Pratt等[20]采用高斯–牛頓算法實現(xiàn)了頻率–空間域的FWI，給出了Hessian矩陣的完整表達式并指出：在高頻近似下逆Hessian矩陣是對角占優(yōu)矩陣，對于包含個參數(shù)的反演問題，只需計算逆Hessian矩陣的個對角元素，這極大地簡化了其的計算過程。C. Shin等[21]基于逆散射理論給出了疊前深度偏移中擬Hessian矩陣(用一對角占優(yōu)矩陣近似代替Hessian矩陣)對角項的估計方法，但缺點是缺失幾何擴散相關(guān)的因子。D. H. Sheen等[22]采用高斯–牛頓算法實現(xiàn)了時間域彈性波全波形反演，并給出了擬Hessian矩陣對角項的計算方法；還有一種思路是利用目標(biāo)函數(shù)及其梯度來構(gòu)建擬Hessian矩陣，其中以l-BFGS[23-24]最為常用。L. Métivier等[13]將FWI視為牛頓線性系統(tǒng)，采用二階伴隨狀態(tài)法來估計更精確的對角項，但由于每次迭代都需要4個正演過程，大大增加了計算時間。為了進一步降低Hessian矩陣的計算成本并提高反演精度，修正的擬牛頓波形反演算法[25]結(jié)合梯度、模型信息及目標(biāo)函數(shù)信息估計了一個投影Hessian矩陣并成功應(yīng)用于淺地表的模型重構(gòu)。Li Jing等[26]在勒夫波波動方程頻散反演中通過對觀測數(shù)據(jù)進行加權(quán)形成一個新的伴隨源來對梯度進行修正。Liu Zhaolun等[27]利用該方法實現(xiàn)了瑞雷波的三維成像。Yan Yingwei等[28-30]基于疊前深度偏移理論中求取逆Hessian矩陣的方法求取預(yù)處理算子，并成功將其應(yīng)用到勒夫波實際數(shù)據(jù)全波形反演中，取得較好效果。

筆者基于逆散射理論中求取Hessian矩陣的方法求取勒夫波FWI目標(biāo)函數(shù)的擬Hessian算子實現(xiàn)對梯度的預(yù)處理，以均衡整體梯度振幅加強算法對深部介質(zhì)的照明能力，提升地層深部反演精度。為解決初始模型拙劣以及波場照明不足導(dǎo)致的反演精度丟失問題，采用偽井地聯(lián)合反演，使地震記錄中攜帶更多的深部介質(zhì)信息，從本質(zhì)上提升成像質(zhì)量。

1 方法原理

1.1 正演問題

在二維各向同性介質(zhì)中，假設(shè)質(zhì)點在-平面上只有方向的非零位移(為水平距離，為深度)，結(jié)合虎克定律，得到偽保守形式的勒夫波彈性動力

學(xué)方程[31]:

式(1)的一般形式為：

式中：為模型參數(shù)；為式(1)的偏微分算子；為地震波場；為加載的震源項。

1.2 局部優(yōu)化算法

FWI計算量巨大，使得全局優(yōu)化和混合優(yōu)化算法難以在FWI領(lǐng)域得到應(yīng)用。因此，只有少數(shù)局部優(yōu)化算法經(jīng)過了測試[34]，即從一個初始模型出發(fā)逐步迭代尋優(yōu)，迭代序列可以表示為：

對于基于梯度的算法，共軛梯度法(Conjugate Gradient method，CG)在FWI中得到了廣泛的應(yīng)用。其模型擾動為：

對于預(yù)處理算法，其核心是求解目標(biāo)函數(shù)的逆Hessian矩陣，新的模型擾動為:

式中：上標(biāo)T代表轉(zhuǎn)置運算。對于預(yù)處理算法中的計算[35]，選用下式：

為了加快目標(biāo)函數(shù)的收斂速度，令標(biāo)量參數(shù)為非負值[36]：

1.3 目標(biāo)函數(shù)梯度及預(yù)處理算子

全波形反演的本質(zhì)思想就是將觀測的地震數(shù)據(jù)與預(yù)測的地震數(shù)據(jù)進行匹配，使某種衡量二者誤差關(guān)系的目標(biāo)函數(shù)值達到最小，從而得到最佳的模型數(shù)據(jù)。目標(biāo)函數(shù)可以定義為：

求取目標(biāo)函數(shù)的梯度是局部優(yōu)化算法的核心。伴隨狀態(tài)法[18]指出，梯度是正傳波場的時間偏導(dǎo)與檢波器位置處反傳波場殘差獲得的反傳波場的零延遲互相關(guān)，這在C. Castellanos等[37]的研究中有詳細的推導(dǎo)說明。因此，結(jié)合鏈?zhǔn)椒▌t[38]，可以得到勒夫波FWI中關(guān)于橫波速度的梯度算子，其差分格式為：

由式(10)求取的梯度由于地震數(shù)據(jù)的頻帶受限、雙重散射以及波場非均勻覆蓋等原因會隨著深度的加深振幅急劇減小，導(dǎo)致對深部地層的照明缺失，介質(zhì)參數(shù)得不到明顯更新，而目標(biāo)函數(shù)的二階導(dǎo)數(shù)(Hessian矩陣)由于包含曲率信息，可清晰預(yù)測這種現(xiàn)象。逆Hessian矩陣，也可稱其為預(yù)處理算子，便可調(diào)節(jié)梯度算子比例并均衡其修正量，消除帶限效應(yīng)，以加強深部地層的照明能力，提高反演精度。由于其維度過于龐大，在當(dāng)前的計算能力下顯式計算這個矩陣幾乎是不可能的，但是可以通過互易定理獲得擬Hessian矩陣。本文參考疊前深度偏移中的擬Hessian矩陣的求解方法來求解時間域勒夫波全波形反演中的擬Hessian矩陣。這種可以壓制由多次散射在梯度向量中產(chǎn)生的假象的求解 Hessian 矩陣的理論，被稱為逆散射理論[21]。根據(jù)逆散射理論，擬Hessian矩陣的表達式為：

由此式即可得到勒夫波FWI中關(guān)于橫波速度的擬Hessian矩陣，其差分格式為：

關(guān)于應(yīng)力時間的導(dǎo)數(shù)在求取梯度時已經(jīng)計算，求取擬Hessian矩陣時只需要相乘再疊加即可。由于逆Hessian矩陣可看作是對角占優(yōu)矩陣，故只需計算出其對角項。同時，為了避免對角項中由于波場的幾何擴散等原因出現(xiàn)近零值導(dǎo)致反演過程不穩(wěn)定，故對其加上一個調(diào)節(jié)因子，第次迭代時的預(yù)處理算子表達式如下：

將式(14)代入式(5)，即形成了預(yù)處理的共軛梯度算法(Pre-conditioned Conjugate Gradient method，PCG)。

2 模型重構(gòu)測試

為了驗證上述PCG比CG更具優(yōu)越性，我們將這2種FWI算法同時應(yīng)用于斷層模型、凹陷模型以及起伏界面模型進行反演試算。另外，在起伏界面模型試 算中加入井地聯(lián)合反演，以證明其在提高反演精度、提升成像質(zhì)量方面的有效性。在反演過程中，認(rèn)為模型密度及地震子波已知。為了定量評價反演結(jié)果精度，對比2種算法的優(yōu)劣，引入模型的對比誤差(RMSE)，其表達式為：

2.1 斷層模型

圖1 斷層模型重構(gòu)測試中的真實模型、初始模型及反演結(jié)果

圖1c、圖1d分別為PCG與CG算法的反演結(jié)果，可以看出2種算法都較好地重建了地層異常，證明了程序的正確性。CG算法可以精確重建模型淺部形態(tài)且給出大致的異常構(gòu)造及其邊界，但不能重構(gòu)模型深部形態(tài)。這是因為隨著深度的加深，梯度算子振幅逐漸變小，模型深部參數(shù)沒有得到明顯更新，可從第22道和第35道處的1維(1-D)橫波速度深度剖面(圖3a、圖3b)得到更為直觀的反映，深層速度偏低。而PCG算法則可以更精確地重建模型，橫波速度的誤差終值僅為0.022 9。圖3c為2種算法的歸一化目標(biāo)函數(shù)值變化曲線的對比圖，PCG的收斂速度明顯快于CG。在前5次迭代，PCG的函數(shù)值迅速收斂至極低水平，說明反演模型已更新到真實模型附近，隨后逐漸逼近真實模型，迭代17次即滿足精度要求，且目標(biāo)函數(shù)終值也更小，僅為0.011 6。而CG前7次迭代的函數(shù)值收斂也較快，這是因為模型淺部的更新對函數(shù)值的貢獻很大，之后由于模型深部的梯度振幅過小，參數(shù)更新不足，函數(shù)值收斂速度變慢，誤差終值為0.047 2。由圖3d兩種算法反演結(jié)果的對比誤差可以看出，在整個迭代過程中，PCG的對比誤差一直低于CG，從0.080 1下降到0.022 9，而CG算法僅能下降到0.047 2(表1)，反演精度遠遠高于CG。圖3e給出的波形擬合圖說明2種方法的反演波形均能很好的擬合于真實波形。

圖3 斷層模型重構(gòu)測試的評價曲線

表1 斷層模型重構(gòu)測試中PCG與CG算法的性能對比評價

在斷層模型反演過程中我們發(fā)現(xiàn)，由于FWI中固有的參數(shù)過度估計與參數(shù)估計不足的缺陷，導(dǎo)致最終得到的反演結(jié)果中，地層深部速度上界高于500 m/s，地表速度下界低于300 m/s(圖3a、圖3b)。這是由于迭代過程中預(yù)處理算子對梯度算子過度調(diào)節(jié)，出現(xiàn)了照明過度的現(xiàn)象。因此，需要在迭代尋優(yōu)過程中加入模型更新的約束控制條件來解決該問題。

2.2 凹陷模型

為進一步對比分析PCG比CG算法具有優(yōu)越性，繼續(xù)對凹陷模型進行反演試算。模型垂直與水平網(wǎng)格數(shù)目為81×101，垂直和水平網(wǎng)格步長均為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。觀測系統(tǒng)與斷層模型基本一致，不同之處是激發(fā)點距改為5 m，保持激發(fā)點數(shù)目11不變(圖4a中的白色箭頭所示)，水平分量接收點為101個(圖4a中的黃色圓點，已抽稀至51個顯示)。波場正演采用斷層模型參數(shù)，正演地震記錄(第6炮)如圖5所示。將初始模型劃分為40層，同樣參考地震記錄中勒夫波的基階頻散曲線高頻端、低頻端相速度，設(shè)置初始模型速度范圍為200～500 m/s，層速度增量為7.5 m/s，如圖4b所示。

在該模型的試算中，仍存在參數(shù)過度估計與估計不足的現(xiàn)象(圖6a、圖6b)。由圖6e的波形擬合圖可以看出，初始波形與真實波形相差較大，說明初始模型比較拙劣，這是導(dǎo)致CG算法的反演波形沒能與真實波形擬合的因素之一，同時也是圖6c中CG算法的目標(biāo)函數(shù)值終值較大，僅能收斂至0.052 1(表2)的原因之一。另一原因則是相較于斷層模型，該模型深度變深，梯度振幅受限更為嚴(yán)重，導(dǎo)致反演精度進一步降低。但PCG算法的反演波形與真實波形完全擬合，同時目標(biāo)函數(shù)值也近似收斂于0(表2)。圖6d中CG算法的橫波速度對比誤差在前5次迭代中與PCG算法相當(dāng)，這是因為模型淺部(10 m以上)的更新貢獻較大。之后PCG算法由于對模型深部照明能力更強，使得深部介質(zhì)參數(shù)得到更新，對比誤差下降較快。圖4c、圖4d分別給出了PCG與CG算法的反演結(jié)果，可以清晰地看出，相較于斷層模型，2種方法的成像質(zhì)量差距明顯加大。PCG算法較精確地重建了地層異常及特殊構(gòu)造，清晰地刻畫出了異常邊界，橫波速度對比誤差終值由0.102 8下降至0.030 6(表2)；而CG算法只能較準(zhǔn)確地給出模型淺層的形態(tài)，模糊地反演出模型中部存在高速凸起、低速凹陷，而模型底部的斷層則完全沒有顯示出來，橫波速度的對比誤差下降量很小，終值為0.088 7(表2)。這一定程度上展現(xiàn)了預(yù)處理算子在提升深部照明強度、提高反演精度方面的作用，說明PCG算法在FWI中的有效性。

圖4 凹陷模型重構(gòu)測試中的真實模型、初始模型及反演結(jié)果

圖5 凹陷模型重構(gòu)測試的第6炮地震記錄

2.3 起伏界面模型

最后對界面起伏模型進行反演試算，如圖7a所示，真實模型由界面起伏的四層速度遞增的地層組成，橫波速度依次為200、300、400、500 m/s，差異已足夠顯著。模型尺寸、觀測系統(tǒng)、波場正演參數(shù)與凹陷構(gòu)造模型完全相同。圖8為獲得的第6炮

地震記錄。將初始模型劃分為20層，同樣參考地震記錄中勒夫波的基階頻散曲線高頻端、低頻端相速度，設(shè)置初始模型速度范圍為200～500 m/s，層速度增量為15 m/s，如圖7b所示。

在該模型的試算中，仍存在參數(shù)過度估計與估計不足的現(xiàn)象(圖9a、圖9b)。圖9c為2種算法的目標(biāo)函數(shù)值變化曲線，CG算法由于對深部介質(zhì)參數(shù)的更新不足，在第21次迭代便已滿足尋優(yōu)終止條件，目標(biāo)函數(shù)終值過大，為0.117 3(表3)。PCG算法的收斂速度快于CG算法，且目標(biāo)函數(shù)終值也更小，為0.011 5(表3)。圖7c、圖7d分別為PCG與CG算法的反演結(jié)果，可以看出CG算法僅僅能重建模型淺部形態(tài)，不能給出模型深部的地層邊界，原因已在前2個模型試算中進行了詳細分析。而PCG算法雖然由于預(yù)處理算子的補償作用清楚地給出了各地層邊界，反演精度有了一定程度的提高，但是由于“波峰”處的地層(圖7c中箭頭所指的部分)波場照明不足，地震波無法穿透，使地震記錄中缺失這部分地層的波場信息，導(dǎo)致參數(shù)更新不足，由第22道處的1-D速度深度剖面(圖9a)可以得到更為直觀的反映。圖9e給出了波形擬合結(jié)果，PCG算法的反演波形與真實波形已完全擬合，但橫波速度對比誤差變化曲線(圖9d)表明2種算法的反演精度均有限，可能是由于短周期的周波跳躍導(dǎo)致的。另外，PCG算法在第6次迭代后，對比誤差不降反升，主要原因是迭代后期對中深部地層照明過度，參數(shù)過度估計(圖9b，第90道處的1-D速度深度剖面)。

表2 凹陷模型重構(gòu)測試中PCG與CG算法性能對比評價

圖7 起伏界面模型重構(gòu)測試中的真實模型、初始模型及反演結(jié)果

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

圖8 起伏界面模型重構(gòu)測試的第6炮地震記錄

表3 界面起伏模型重構(gòu)測試中PCG與CG算法性能對比評價

圖9 界面起伏模型重構(gòu)測試的評價曲線

2.4 偽井地聯(lián)合反演

為了進一步提高FWI成像質(zhì)量，本文開展了偽井地聯(lián)合(Pseudo Borehole-Ground joint，PBGJ)的PCG算法測試。通過在地面和井中聯(lián)合激發(fā)震源，加強模型深部的波場照明。由于地震數(shù)據(jù)包含了模型深部介質(zhì)的豐富信息，反演精度將得到提高。將PBGJ-PCG算法應(yīng)用于起伏界面模型測試中，除了地表激發(fā)11激發(fā)點外，在垂直方向上距離第1個激發(fā)點25 m處再增加1個激發(fā)點(圖7a所示白色箭頭)，其坐標(biāo)為(0 m, 25 m)。

圖10c給出了PBGJ-PCG算法反演結(jié)果的評價曲線，由于井地聯(lián)合中地震記錄包含的地層信息更為豐富，在前期的迭代尋優(yōu)中便對模型整體進行了更新，而模型淺部的更新對目標(biāo)函數(shù)值貢獻很大，導(dǎo)致只能對模型淺層進行更新的PCG算法的目標(biāo)函數(shù)收斂速度比PBGJ-PCG算法更快，但目標(biāo)函數(shù)終止值更大(表4)，對比誤差也一直更高(圖10d)，但反演波形均完全擬合于真實波形(圖10e)，推測是由于周波跳躍導(dǎo)致的。由圖10a給出的反演結(jié)果及圖10b所示的第22道處的1-D速度剖面可知，PBGJ-PCG算法提高了“波峰”處地層的反演精度，介質(zhì)參數(shù)得到了明顯更新，較精確地刻畫了地層邊界，重構(gòu)了地層速度。對比誤差由0.089 3下降至0.033 4，而PCG算法的對比誤差只下降至0.078 5(表4)，展現(xiàn)了PBGJ-PCG算法在提高反演精度方面的有效性。

3 結(jié)論

a.擬Hessian梯度預(yù)處理算子可對梯度進行能量均衡，提高對深部介質(zhì)的照明能力，對FWI具有重要意義。同時，由此衍生出的PCG算法相較于經(jīng)典的CG算法，可加快收斂速度，即使初始模型比較拙劣，仍能比較準(zhǔn)確地給出地層的異常結(jié)構(gòu)及邊界，提高深部介質(zhì)的反演精度，提升成像質(zhì)量。

b. PBGJ-PCG雖然無法解決時間域FWI固有的周波跳躍問題，但該方式獲得的地震記錄中攜帶了更多的深部介質(zhì)信息，從本質(zhì)上提高了對深部介質(zhì)的照明能力，可進一步提升PCG算法的成像質(zhì)量。

c. PCG算法與井地聯(lián)合反演方式雖然在本文中僅應(yīng)用于勒夫波的全波形反演中，但在瑞利波及體波的FWI中同樣適用。

表4 PBGJ-PCG與PCG算法性能對比評價

d. PCG算法可能存在對梯度算子的過度調(diào)節(jié)，導(dǎo)致參數(shù)過度估計的現(xiàn)象出現(xiàn)。因此，迭代尋優(yōu)過程中需要加入模型更新的約束控制條件，這需要進一步深入研究。

[1] 張胤彬，潘冬明，胡明順，等. 復(fù)雜地表初至層析地震靜校正技術(shù)及其應(yīng)用[J]. 煤田地質(zhì)與勘探，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] 周俊杰，王雨，侯瑋. 黃土塬地區(qū)煤田三維地震綜合處理技術(shù)[J]. 地球物理學(xué)進展，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]. 武漢：中國地質(zhì)大學(xué)出版社，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]. 地球物理學(xué)報，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] 閆英偉. 淺地表高頻面波成像技術(shù)研究[D]. 長春：吉林大學(xué)，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] 劉聰，王者江，閆英偉. 基于伴隨狀態(tài)法二維時間域勒夫波全波形反演研究[J]. 地球物理學(xué)進展，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

國家重點研發(fā)計劃課題(2018YFC0807803)；陜西省自然科學(xué)基礎(chǔ)研究計劃項目(2019JLM8)；長安大學(xué)中央高校基本科研業(yè)務(wù)費專項基金項目(300102260203)

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

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

管建博，李宇，殷裁云，等. 基于擬Hessian梯度預(yù)處理算子的勒夫波全波形反演研究[J]. 煤田地質(zhì)與勘探，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

(責(zé)任編輯 聶愛蘭)