頻率域八階NAD有限差分模擬及全波形反演

韓如冰 郎 超

(①中國地質科學院地質研究所自然資源部深地動力學重點實驗室，北京 100037； ②中國地質大學(武漢)地球物理與空間信息學院，湖北武漢 430074)

0 引言

全波形反演(Full-waveform Inversion)以波動方程作為數學模型模擬地震波的傳播規律，能夠更充分地利用地震數據信息，具有成像更準確、分辨率更高的特點[1-2]。Tarantola[1]首先提出了基于廣義最小二乘的時間域全波形反演理論和方法，隨后Pratt[3-4]將其推廣到頻率域內。相比于時間域，頻率域全波形反演具有計算過程穩定、占用內存小、不存在累計誤差、緩和局部極小、降低反演的非線性、易于處理衰減頻散效應、易于并行計算等特點[5-10]，近些年得到廣泛發展與應用[11-12]。

在頻率域計算中，選擇合適的數值方法離散波動方程，處理區域邊界、求解大型線性方程組是必要的。常見的數值離散方法包括有限差分法[13-14]、有限元法[15-16]、譜元法[17-18]等。當網格不夠密集時，常規的有限差分(Ordinary Finite-difference，OFD)方法容易出現明顯的數值頻散現象，盡管通過加密離散網格可以在一定程度上壓制頻散，但是相應的計算量增大，所需時間和內存也會增加[19]。

針對以上問題，Yang等[20-21]引入了近似解析離散化(Nearly Analytic Discrete，NAD)方法，并應用于時間域地震波模擬，隨后Lang等[8]將其引入到頻率域正演模擬中，其基本思想是利用波場位移與其梯度值的聯合近似得到位移的高階偏導數，進而離散波動方程。相比于其他差分方法，NAD算子攜帶了更多的波場信息，特別是含有刻畫波場變化特征的梯度信息。前人研究表明，針對較粗的離散網格，NAD方法在壓制數值頻散方面具有很好的效果[8]，此外，與同階差分格式相比，NAD數值格式計算量小，易于并行計算。然而現有的頻率域四階NAD(NAD4)方法對于計算效率的提升還很有限，需要考慮在此基礎上構造更加有效的數值方法。

為提高全波形反演的計算效率，本文構造了頻率域八階NAD(NAD8)格式離散二維聲波方程，采用完全匹配層吸收邊界(Perfectly Matched Layer，PML)[22]條件。首先詳細介紹了NAD8格式的具體離散過程；然后使用一類不精確旋轉分塊三角預處理算子(IRBTP)加速廣義極小殘量方法(GMRES)[23]求解離散后的大型稀疏線性方程組；通過數值頻散分析以及波場模擬與六階NAD(NAD6)、八階OFD(OFD8)方法進行對比，結果表明NAD8方法在壓制頻散和提高計算效率等方面具有優勢；最后，基于共軛梯度法[24]，首次構造了基于NAD8方法的頻率域全波形反演算法，對兩種典型分層介質模型和經典Marmousi模型進行了反演試算，獲得了高保真、高分辨率的成像結果，進而驗證了所提方法的有效性和適用性。

1 方法原理

歸結起來，頻率域NAD方法的離散過程可概括為以下幾步：對頻率域波動方程關于各個方向求偏導，并在計算區域的邊界處施加吸收邊界條件，進而得到偏微分方程組，采用NAD網格差分模板離散其中的高階偏導數項，再按一定節點排序規則形成大型稀疏線性方程組。

假定地球為完全彈性介質時，可以用聲波方程研究地震波的傳播規律。常密度介質下的二維頻率域聲波方程為

(1)

式中：u(x，z，ω)表示位移或者壓力場；ω=2πf表示角頻率；c為地震波波速；s表示震源項。

對式(1)兩端分別關于x方向與z方向求偏導，有

(2)

與時間域不同，在頻率域施加PML邊界條件的過程中，不需要進行反褶積計算，實現更方便、準確。首先引入復坐標控制衰減，以x方向為例，有

(3)

將衰減函數取為

(4)

(5)

式中

(6)

與NAD4方法[19]不同，對于NAD8格式，任一點(i,j)的高階偏導數需要由(i-2，j)、(i-1，j)、(i，j)、(i+1，j)、(i+2，j)五個點處的波場及其空間梯度的加權組合近似。求解微分方程中高階偏導數的離散格式為

(7)

(8)

(9)

(10)

128ui-1,j+1-128ui+1,j-1+128ui-1,j-1-

7ui+2,j-2+7ui-2,j-2)+

(11)

1408ui+1,j+1-1408ui-1,j+1-62ui+2,j-

2816ui+1,j+2816ui-1,j+62ui-2,j+

31ui+2,j-2+1408ui+1,j-1-1408ui-1,j-1-

(12)

1408ui+1,j+1-2816ui,j+1+1408ui-1,j+1-

1408ui+1,j-1+2816ui,j-1-1408ui-1,j-1-

31ui+2,j-2+62ui,j-2-31ui-2,j-2)+

(13)

式中h表示空間離散步長。將離散結果依次代入式(5)，并寫成統一線性方程形式

(14)

(15a)

(15b)

(16a)

(16b)

(17a)

(17b)

由式(14)的離散結果發現，NAD8差分格式共涉及周圍17個節點，顯然算子長度小于同階精度的其他格式[25]，具體差分模板如圖1所示。

為快速有效求解離散后的大型稀疏線性方程組，本文采用一類不精確旋轉分塊三角(IRBT)預處理算子加速廣義極小殘量方法(GMRES)。該方法是一種專門用于非對稱線性方程組求解的Krylov子空間方法，充分利用了系數矩陣T的稀疏結構，結合不精確預處理的思想，在迭代次數并未顯著增加的前提下，可節省較多的運算時間與存儲量，從而能明顯提高計算效率[23]。

圖1 NAD8差分網格模板

2 正演計算

從數值頻散分析和波場模擬兩方面考查NAD8方法相對于NAD6和常規八階有限差分方法的模擬效率。震源選取雷克(Ricker)子波，根據Fourier變換的性質，其頻率域表達式[9]為

(18)

式中：A表示振幅，波場模擬時統一設置為1；f0是主頻。根據頻譜分析[8]，主要能量分布在0～3f0。定義網格頻率(Gf)為每個波長包含的離散網格點數，以度量網格剖分的大小[26]，其計算公式為

(19)

式中w表示波長。從式(19)可以看出，對于同一計算區域，Gf越大，每個波長內網格點數目會越多，從而計算精度越高，但是相應地增大了計算量。

2.1 理論頻散曲線分析

首先從理論上論證各種數值方法在壓制數值頻散方面的能力。基本思想為：通過計算數值速度與實際波速的偏差考察各種離散方法的準確性，偏差越小說明頻散越小。

為計算NAD方法的數值速度，首先引入均勻介質中的波動方程及其關于x、z方向的偏導

(20)

與式(5)的離散過程類似，通過NAD方法可將式(20)離散為矩陣形式

(21)

(22a)

(22b)

式中：a=exp(ikxh)；b=exp(ikzh)。

圖2 θ=0(a)和θ=π/6(b)時三種數值方法的頻散曲線

2.2 簡單模型的波場模擬

以均勻介質模型(M1)和雙層介質模型(M2)，進一步討論NAD8方法的計算精度與數值效率。

M1模型的速度v=4km/s； M2模型的速度v1=4km/s、v2=5km/s，界面位于z=4km處，其他參數設置如表1所示。

表1 M1和M2模型的基本參數

圖3為頻率f=10和30Hz時， 本文構造的NAD8方法針對以上兩種模型得到的頻率域單頻波波場快照。圖4為NAD8、NAD6和OFD8三種方法模擬的t=0.5s時刻時間域波場快照。從圖中可以看出，在相同網格規模的情況下， NAD8方法有效地壓制了數值頻散，而后兩種方法則出現了明顯的頻散現象，OFD8方法頻散最為嚴重，表明在相同網格頻率下， NAD8方法較其他方法精度較高。根據Shannon采樣定理[27]，為了準確地刻畫波場，網格頻率至少大于2.0， 而NAD8方法突破了這個極限，Gf最小值可到達1.98，這是因為NAD方法除了利用位移信息，還利用了波場梯度信息，大大提高了波場模擬精度，降低了對網格頻率的要求。

圖3 由八階NAD方法計算的M1(左)和 M2(右)模型頻率域單頻波波場快照(a)10Hz； (b)30Hz

圖4 三種方法計算的M1(左)和M2(右) 模型t=0.5s時刻的時間域波場快照(a)NAD8； (b)NAD6； (c)OFD8

最后，針對以上兩個模型，在不產生數字頻散的情況下，統計三種方法運行時間(表2，10次計算結果的平均)，其中NAD8(Gf=1.98)用時最少，比NAD6方法(Gf=2.40)約省時12%，比OFD8方法(Gf=3.50)約省時25%。

表2 不同方法簡單模型波場模擬運行時間對比

2.3 復雜模型的波場模擬

選取經典的Marmousi模型[28](圖5)檢驗NAD8方法的適用性。模型尺寸為17.25km×5.63km；網格間距h=18.75m，網格頻率Gf=3.0；四周均為10層PML邊界，震源位于(8.63km，0.38km)，主頻為20Hz。

圖6為NAD8方法模擬的頻率f=10Hz和15Hz時Marmousi模型的頻率域單頻波波場快照，圖7為t=1.00、1.67s時刻的時間域波場快照。可以看出，對于復雜的Marmousi模型，沒有明顯的頻散，可見NAD8方法適用于復雜介質模型的正演模擬。

圖5 Marmousi速度模型

圖6 NAD8方法模擬的Marmousi模型單頻波波場快照(a)f=10Hz； (b)f=15Hz

圖7 NAD8方法模擬的Marmousi模型時間域波場快照(a)t=1.00s； (b)t=1.67s

3 頻率域全波形反演

3.1 基于NAD8方法的頻率域全波形反演

所謂反演即從初始模型出發，進行正演模擬，匹配計算波場與真實波場，然后不斷修正模型，最終達到特定的精度。本文基于NAD8正演算法對經典模型進行反演，反演采用非線性共軛梯度法，其優勢在于：與最速下降法相比，不用引入顯著的正演計算量便可以得到更加準確、穩健的計算結果，與牛頓類方法相比，不用計算反演目標函數的二階偏導數矩陣(Hessian矩陣)，收斂性較強。將離散后的頻率域聲波方程(式(14))改寫為

T(ω)u(ω)=s(ω)

(23)

式中：T表示系數矩陣；u表示波場項；s表示震源項，且三者都與角頻率ω有關。合成波場與觀測波場數據殘差的L2范數可表示為

(24)

式中：m為模型參數； “*”表示復共軛；ns為炮數；nr為每炮接收點數。δdij表示第i炮的第j個檢波點對應的模擬波場值u與觀測波場d的殘差。選擇合適的方法求解E(m)的最小值是問題的關鍵，基于梯度算法，通過推導計算可得相鄰步數之間的迭代基本格式為

(25)

通過對波動方程關于模型參數求偏導，進而可導出

(26)

式中v表示殘差反傳波場。可見計算反演目標函數梯度時候，不需要計算矩陣J的具體表達式，整個問題只需計算TTv=δd*。復雜的矩陣計算可以變成簡單的求解一次線性方程組的問題，大大降低了計算量。

通過以上過程計算出梯度之后，進而可以確定搜索方向。本文選用計算效率較高的拋物線擬合方法選擇步長，即近似認為目標函數是關于步長的二次函數，確定一個試探步長r1，使其滿足在目標函數的下降方向上，不斷重復搜索步驟，直到目標函數不降反增，確定另一步長r2； 結合初始步長r0，三點確定拋物線，求取極值點便為搜索步長。最后，反演頻率的選擇參照非等距頻率選擇策略[29]，這樣頻率間距隨著頻率的增加變得越來越大，從而整體上越來越松散，但不損失反演結果的分辨率，使用更少的頻率參與反演，可節省計算時間。

此外，定義標準化殘差的概念來衡量反演誤差的大小，即

(27)

式中：mt表示真實模型； ‖·‖2表示L2范數。顯然給定某一頻率，κ是一個隨迭代步數變化的曲線，曲線越趨近于0表示反演結果越好。

3.2 簡單模型反演

基于NAD8方法的正演程序，依據上文所述的頻率域反演原理，首先對兩個簡單的速度模型進行反演與分析。

模型一(圖8)參數設置如下：模型尺寸為2.5km×2.5km，背景速度為4.0km/s，中央有一速度為4.5km/s、邊長為0.5km的異常體，網格間距h=25m； 將四周設置為10層PML邊界，主頻設置為10Hz，震源振幅為1.0×105(為了避免波場值太小而引起的機器出錯)，初始速度設置為4.0km/s。炮點與接收點布設兩組，位于計算區域內。在頻率域接收點的多少并不會增大反演的計算量，而震源數目的增多會使計算量成倍增大。為保證計算效率，需要選擇適量的震源數目，而盡量多布設接收點，橫向間距為一個網格距離。為了得到深部更高分辨率的成像結果，分別設置兩排震源和接收器，設置震源數目為上、下各11個，共計22個，縱向上分別位于0.275km和2.225km處，橫向上為0.25～2.25km，間距為0.2km。上、下接收點同時接收傳播信號。

圖8 速度模型一及震源與接收點布設示意圖

首先選取5個頻率(1、7、13、19、25Hz)，低頻勾勒大致輪廓，高頻細化內部結構，每個頻率的最大迭代步數為30步，將反演過程記為S1，最終f=25Hz的反演結果如圖9a 所示。可以看出反演的速度模型與理論模型整體輪廓基本一致，但是分辨率相對較低，特別在靠近速度異常體周邊的數值與理論值差異較大。進一步加密震源，每行放置21個，即2×21個，增加反演頻數到30個(1～30Hz，間隔1Hz)，每個頻率最大的迭代步數為50，將反演過程記為S2，最終f=30Hz的反演結果如圖9b所示。與圖9a相比，圖9b的分辨率明顯提高，細節刻畫更清晰，對異常體的刻畫更準確。

圖9 模型一反演結果(a)S1； (b)S2

圖10中藍色曲線是頻率取25Hz時S1反演的誤差曲線，可以看出，當步數到達最終30時，誤差曲線還有下降趨勢；紅色曲線為頻率取25Hz時S2反演的誤差曲線，表現為前幾步急劇下降，后期逐步平緩，停滯在小于0.1范圍內且基本不再下降，1～30Hz中其他頻率的誤差曲線類似。

圖10 模型一反演誤差曲線

模型二(圖11)為雙層模型，尺寸為2.5km×2.5km，上層速度為4.00km/s，下層速度為4.50km/s在分界面處有一速度為4.25km/s、邊長為0.5km的異常體，網格間距為h=25m。將四周設置10層的PML邊界，震源主頻設為10Hz，震源振幅設置為1.0×105。反演的初始速度設置為4.00km/s。接收器盡可能多的放置，間距為一個網格距離，震源數目為上、下各布設41個，共計82個，縱向上分別位于0.275km和2.225km處，橫向上從0.25km到2.25km，間距為0.05km。上、下接收點同時接收傳播信號。

設置反演頻數為30個(1～30Hz，間隔為1Hz)，每個頻點的最大迭代次數為100。圖12a為f=30Hz的反演結果，可以看出，整體分辨率較高，細節刻畫清晰，特別是在速度分界面與異常體的形態反映較為真實，整體效果較為理想；圖12b為f=30Hz的反演誤差曲線，開始時下降明顯，40次后穩定在0.1以內，符合高精度反演的誤差曲線形態。可見，本文的基于NAD8的頻率域反演方法對于簡單的模型一和模型二取得了較為理想的結果。

圖11 速度模型二及震源與接收器布設示意圖

3.3 Marmousi模型反演

針對Marmousi模型(圖5)的反演，參數設置如下：網格間距h=15m，反演區域為3.47km×1.14km，PML層數為10，主頻設置為10Hz，震源振幅設為1.0×105。初始模型(圖13)是對真實模型進行Gauss平滑后的結果，只具有大致的分層結構。震源布設44個，縱向上位于0.17km處，橫向上從0.15km到3.35km， 間距為0.075km。僅在地表放置一排。接收器排列與震源的間隔相同，縱向位置上位于0.15km處，橫向上往外多擴了兩個網格距離。

選擇30個頻點(1～30Hz，間隔為1Hz)、每個頻點迭代200次進行反演，圖14為10和30Hz的反演結果。由圖可以看出，10Hz低頻反演結果基本給出了模型的大體輪廓，但是整體分辨率較低，出現了多處異常亮點，特別是深部高速異常體的精度較低。30Hz反演結果分辨率明顯提高、細節更加清晰，特別是深部高速異常體的刻畫更加精確。

圖12 模型二反演結果(a)及反演誤差曲線(b)

圖13 初始Marmousi速度模型及震源與接收點布設示意圖

圖14 Marmousi模型10(a)和30Hz(b)的速度反演結果

圖15是x=2.0km處、頻率為30Hz的速度反演曲線與真實、初始速度曲線的對比，可以看出，整體的反演結果較為理想，特別是淺部低速體保真度較高，但是深層精度略有下降，導致反演速度與真實模型速度之差的L2范數略高(約為0.8)，這是因為震源位于地表處，波傳播到深層時振幅降低、誤差累積，從而分辨率下降。

圖16為頻率為30Hz的反演誤差曲線，可以看出，整體上一開始下降明顯，并逐漸穩定在0.1以內，反演效果較好，其他頻率反演誤差曲線基本相似。可見，基于NAD8的頻率域全波形反演方法針對復雜模型依然適用，且分辨率、保真度都較高，效果良好。

圖15 x=2.0km處反演曲線與真實模型、初始模型的對比

圖16 Marmousi模型反演誤差曲線

4 結論

為進一步提高頻率域全波形反演的計算效率，本文首次構造了頻率域NAD8差分格式離散波動方程，給出了PML邊界條件與頻率域NAD方法相結合的正演算法，并詳細推導總結了整個數值離散過程，得到了大型線性代數方程組，根據方程組系數矩陣的數學結構，采用一類不精確旋轉分塊三角預處理算子加速Krylov迭代方法，提高了正演計算效率。數值頻散分析以及均勻介質和雙層介質模型的波場模擬結果表明：①與OFD8方法、NAD6方法相比，NAD8方法計算精度最高，在壓制頻散方面效果最明顯； ②在有效壓制頻散的情況下，NAD8方法所需網格頻率最小(Gf=1.98)，突破了Shannon采樣定理可達到的理論最小值(Gf=2.0)，所用的計算時間最短，大約比NAD6方法(Gf=2.4)省時12%，比OFD8方法(Gf=3.5)省時25%。

對于兩種典型的分層介質模型及經典Mar-mousi模型，本文基于NAD8格式的頻率域全波形反演算法獲得了較高分辨率和保真度的結果，結合反演誤差曲線驗證了方法的正確性和適用性。