一種基于譜共軛梯度法的時間域全波形反演

杜 鑫

(中國石化石油工程地球物理有限公司華東分公司，江蘇南京 211112)

全波形反演一詞最早是由Pan等人在1986年提出的，全波形反演的核心思想是基于上世紀(jì)80年代初Tarantola[1]提出的利用正傳波場與反傳波場的廣義最小二乘擬合而形成的。1988年，Tarantola[2]在已有的全波形反演基礎(chǔ)上，給出了完整的時間域全波形反演方法的理論框架。隨后，Bunks[3]等引入多尺度全波形反演思想，使反演結(jié)果更加穩(wěn)定。目前，全波形反演已經(jīng)從時間域發(fā)展到頻率域并取得良好的效果[4-7]。

全波形反演目前應(yīng)用較多的求解方法是收斂速度較快的梯度類局部優(yōu)化迭代反演方法,主要包括最速下降法、共軛梯度法、擬牛頓法等[8-9]。Nocedal[10]提出了擬牛頓法(BFGS)，對收斂速度加速的同時能提高深層地層的分辨率，Pratt[11]等提出利用Hessian矩陣預(yù)處理梯度值以改善迭代的效果，張廣智[12]等借鑒Beck[13]等的思想，對共軛梯度法的迭代算法進行改造，提出快速共軛梯度法并將這種方法應(yīng)用到頻率域全波形反演中。

采用林穗華[14-17]等提出的方法，通過增加譜因子的方式對共軛梯度法進行改造，結(jié)合Hestenes和Stiefe提出的算法與Polak、Ribieren和Polyakn提出的算法優(yōu)點形成全新的修正譜共軛梯度法，并用這種方法來改善聲波全波形反演效果。Marmousi模型的測試結(jié)果表明，該方法能在改善收斂效果的基礎(chǔ)上加快迭代收斂速度。

1 時間域全波形反演基本原理

全波形反演既可以在時間域?qū)崿F(xiàn)也可以在頻率域或拉普拉斯域?qū)崿F(xiàn)，具有在復(fù)雜地質(zhì)條件下精確刻畫構(gòu)造細(xì)節(jié)及巖性變化的潛力。

1.1 時間域全波形反演基本框架

全波形反演是一種高精度速度建模方法，全波形反演包含初始模型的建立、地震波場正演數(shù)值模擬、目標(biāo)函數(shù)的建立、優(yōu)化方法的選擇等。全波形反演基本步驟為：①給定初始模型，并對其做正演模擬，設(shè)定迭代停止的精度要求及最大迭代次數(shù)；②利用正演模擬合成的地震記錄與實際觀測的地震記錄建立目標(biāo)函數(shù)；③判定目標(biāo)函數(shù)反演誤差值是否滿足終止條件或迭代次數(shù)達到最大設(shè)定值，滿足條件則輸出該模型并作為最終反演結(jié)果；若不滿足則轉(zhuǎn)至④；④利用優(yōu)化算法對目標(biāo)函數(shù)的最小二乘結(jié)果進行迭代求解，利用梯度類局部優(yōu)化算法進行求解，更新速度場；⑤利用模型更新方程和速度模型，轉(zhuǎn)至①。

1.2 時間域聲波方程正演計算

時間域各向同性介質(zhì)聲波方程為：

(1)

在使用基于交錯網(wǎng)格的有限差分算法進行數(shù)值模擬時，需要將聲波方程轉(zhuǎn)化為一階聲波方程：

(2)

式中：v為縱波速度，m/s；vx,vz分別為x和y方向速度分量，m/s；P為地層壓力分量，Pa；κ=ρv2為體積模量在空間的分布，Pa；t為傳播時間，s；x為平行于地表方向的距離，m；z為垂直與于地表的方向距離，m。

在有限差分地震波數(shù)值模擬過程中，需要引入具有吸收衰減屬性的邊界進行數(shù)值截斷，通常采用完全匹配層(PML)吸收邊界條件。其基本思想是將地震波場分裂為兩組:一組垂直于波場傳播方向，另一組平行于波場傳播方向。對到達PML邊界的波場，將平行界面法向傳播的平面波進行迅速衰減，而對于其他方向傳播的波并不進行衰減。在進行數(shù)值模擬時，衰減因子采用Collino推導(dǎo)的衰減公式如下：

(3)

式中：vmax為最大縱波速度，m/s；x為離PML內(nèi)層邊界計算點的距離，m；δ為PML匹配層厚度，m；R為介質(zhì)的反射系數(shù)。

1.3 時間域全波形反演

類似于其他反演問題，全波形反演的最小二乘意義下的目標(biāo)函數(shù)可寫為：

(4)

式中：dsyn為正演模擬得到的波場，dobs為實際地震數(shù)據(jù)的炮集波場。

通過迭代的方法更新初始速度模型，即：

vk+1=vk+1+αkdk

(5)

式中：α為迭代步長(可以通過線性搜索、拋物線插值等最優(yōu)化方法取得)，k為迭代次數(shù)，dk為搜索方向。

全波形反演多次迭代極小化實際觀測數(shù)據(jù)與合成地震數(shù)據(jù)的殘差即目標(biāo)函數(shù)，是實現(xiàn)反演速度場v的過程。目標(biāo)函數(shù)的極小化以初始速度v0為前進方向進行搜索，通過多次迭代使目標(biāo)函數(shù)在v0附近收斂到局部極小值。

2 譜共軛梯度法(SCG)

共軛梯度法的搜索方向表達式的構(gòu)成為當(dāng)前次目標(biāo)函數(shù)的梯度與前一次搜索方向的線性組合，具體公式如下：

(6)

式中：βk為共軛梯度修正因子；gk為梯度方向。

為了提高共軛梯度法的計算效率，2001年Birgin和Martinez[18]結(jié)合譜梯度方法和CS，提出了SCG，其搜索方向如下：

(7)

式中：θk為譜系數(shù)。

當(dāng)時θk=1時，式(7)轉(zhuǎn)換為式(6)，實現(xiàn)譜共軛梯度法向?qū)?yīng)的共軛梯度法的轉(zhuǎn)換，因此共軛梯度法的優(yōu)點自然地出現(xiàn)在譜共軛梯度法上。譜共軛梯度法適合求復(fù)雜的無約束優(yōu)化問題，能夠得到更好的數(shù)值計算結(jié)果。

林穗華提出的譜共軛梯度法如下：

(8)

利用譜梯度法進行全波形反演，需要計算梯度方向譜因子，還需要確定迭代步長ak。常用的步長求取方法有線搜索方法和拋物線搜索方法。本文采用的是改進弱Wolfe非精確線搜索方法。其具體條件為：

(9)

(10)

利用線搜索計算步長時，在給定搜索區(qū)間[0,αmax]求取能夠使式(9)和(10)成立的最大步長αk。其中式(9)的目的是要求函數(shù)下降值大于沿切線下降。式(10)可以防止步長選取過小，確保函數(shù)值有足夠的下降量。

共軛梯度算子取值為譜PRP算法以及譜HS算法的最小值，該方法結(jié)合了PRP共軛梯度法，能較好地克服小步長時收斂緩慢以及HS共軛梯度法不依賴于步長的優(yōu)點，都滿足共扼關(guān)系，具有充分的下降性。具體流程如下：①給定初始速度模型v1,終止條件ε>0，k=1，令初始方向d1=-g1；②判別gk與ε的關(guān)系，決定是否終止計算；③利用線搜索方法計算初始步長αk；④由式(8)計算θk，βk,再由式(7)計算搜索方向dk+1(④對算法的效果影響最大)；⑤令k=k+1,計算vk+1。轉(zhuǎn)回②。

3 模型測試

為了測試譜共軛梯度法全波形反演效果，利用Marmousi模型進行測試。不同算法選用相同的反演參數(shù)，雷克子波主頻為5 Hz，空間采樣間隔為dx=dz=25 m，正演時間采樣為2 ms，單道記錄長度8 s，迭代40次。本次模型實驗采用PML吸收邊界條件，其橫向與縱向的網(wǎng)格數(shù)均設(shè)置為50。采用地表放炮，共41炮，炮間距為250 m，共481個檢波器，檢波點道距為25 m。計算平臺主要配置為：處理器i5 3570，8 GB內(nèi)存，500 GB硬盤存儲，Ubuntu18操作系統(tǒng)。圖1是Marmousi模型和初始速度模型。算法目前只支持二維地震數(shù)據(jù)以及速度模型輸入。

圖2是Marmousi模型兩種方法迭代反演的速度模型對比。由圖1和圖2可以看出，兩種方法都能得到基本符合目標(biāo)模型的反演結(jié)果。譜共軛梯度法(用時56 h)用時略少于共軛梯度法(用時60 h)，深層反演結(jié)果的準(zhǔn)確性及分辨率更高，更接近實際速度模型(圖2圈住部分)，驗證了本文所提方法的優(yōu)勢。

圖1 Marmousi模型(上)與初始速度模型(下)對比

圖2 譜共軛梯度法反演結(jié)果(上)與共軛梯度法反演結(jié)果(下)對比

在進行全波形反演迭代時，初始搜索方向均為負(fù)梯度方向。將每次迭代誤差值除以初始誤差值得到歸一化誤差，通過對比CG與SCG方法的歸一化誤差曲線(圖3)，可見SCG方法在22次迭代時的誤差與傳統(tǒng)CG方法在35次迭代時的誤差相近，SCG方法收斂速度更快且收斂過程相對穩(wěn)定。

圖3 兩種方法反演誤差曲線對比

相比于CG方法，SCG方法的額外計算量為譜因子θk的計算，其計算公式中的共軛因子βk已知，而gkTdk-1及范數(shù)求解在共軛因子計算的過程中可以得到。因此，譜因子計算過程中的額外計算量主要為數(shù)乘運算，與整個全波形反演計算過程中的矩陣乘法相比完全可以忽略。

4 結(jié)論

(1)譜共軛梯度法具有共軛梯度法的全部優(yōu)點，得益于譜系數(shù)的引入，利用了更多的目標(biāo)函數(shù)信息，相對于共軛梯度法有更好的數(shù)值表現(xiàn)。

(2)Marmousi模型的反演計算結(jié)果在提高速度模型反演結(jié)果準(zhǔn)確性的基礎(chǔ)上能夠加快迭代收斂速度。