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

杜 鑫

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

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

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

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

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

全波形反演既可以在時間域實現也可以在頻率域或拉普拉斯域實現，具有在復雜地質條件下精確刻畫構造細節及巖性變化的潛力。

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

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

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

時間域各向同性介質聲波方程為：

(1)

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

(2)

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

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

(3)

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

1.3 時間域全波形反演

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

(4)

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

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

vk+1=vk+1+αkdk

(5)

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

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

2 譜共軛梯度法(SCG)

共軛梯度法的搜索方向表達式的構成為當前次目標函數的梯度與前一次搜索方向的線性組合，具體公式如下：

(6)

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

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

(7)

式中：θk為譜系數。

當時θk=1時，式(7)轉換為式(6)，實現譜共軛梯度法向對應的共軛梯度法的轉換，因此共軛梯度法的優點自然地出現在譜共軛梯度法上。譜共軛梯度法適合求復雜的無約束優化問題，能夠得到更好的數值計算結果。

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

(8)

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

(9)

(10)

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

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

3 模型測試

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

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

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

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

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

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

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

4 結論

(1)譜共軛梯度法具有共軛梯度法的全部優點，得益于譜系數的引入，利用了更多的目標函數信息，相對于共軛梯度法有更好的數值表現。

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