基于優化數值積分的譜元法模擬地震波傳播

孟雪莉 劉少林 楊頂輝 汪文帥 徐錫偉 李小凡

(①應急管理部國家自然災害防治研究院，北京 100085； ②寧夏大學數學統計學院，寧夏銀川 750021；③清華大學數學科學系，北京 100084； ④中國地質大學(武漢)地球物理與空間信息學院，湖北武漢 430074)

0 引言

隨著計算機技術的發展，基于地震波運動方程的伴隨成像逐漸被應用于實際地震數據的成像[1-4]。在該類成像方法中，利用震源波場與接收器波場做互相關，從而構建目標函數關于模型參數的Fréchet導數。無論是震源波場的求解還是伴隨波場的求解，都需要數值求解地震波運動方程，因此地震波運動方程的數值求解也成為伴隨成像的關鍵。隨著研究的深入，伴隨成像的研究對象越來越復雜[5]，如成像模型具有強烈起伏界面[6]、模型含有固體與液體耦合[3]。針對復雜模型的伴隨成像，亟待研發適用性強、計算精度高的正演模擬方法。

關于地震波運動方程的直接求解，人們應用最早且最普遍的是有限差分方法(Finite Difference Method，FDM)[7]。FDM的主要優點包括數值格式簡單、程序易于實現、計算速度快等[8-9]，在實際伴隨成像中得到了較廣泛應用[10-11]。但FDM在模擬地震波時存在明顯不足，主要表現為網格不靈活、自由邊界條件難處理和數值頻散強。

針對FDM的不足，人們做了大量研究，也取得了具有重要意義的成果。為了增強FDM網格的靈活性，常采用非規則曲線網格將模型離散，再通過坐標變換方式將曲線網格變換成正交網格，最后在坐標變換后的規則網格上利用FDM求解地震波動方程[12-13]。雖然通過坐標變換方式能有效解決FDM網格不靈活的缺陷，但在坐標變換后的規則網格上使用的有限差分離散格式往往只能是二階中心差分，難以采用高階差分格式，以至于在粗網格條件下制約了FDM數值精度。雖然采用細網格能提升數值模擬精度，但此時需要較大的計算量和存儲量。此外，隨著模型復雜程度的增加，曲線網格在擬合復雜地形時易發生畸變，該網格畸變會嚴重降低曲線網格FDM的數值穩定性[5]。

為了近似自由邊界條件，可采用在自由地表之上設置虛擬層方式構造自由地表法向應力為零的自由邊界條件[13-14]。雖然通過設置虛擬層可方便實現FDM對自由邊界條件的近似，但該數值近似格式只具有二階精度[14]，導致地表處波場模擬精度不足，尤其是面波[15]。

為了減弱FDM的數值頻散，研究者常構造優化的FDM[8，16-21]。通過構造某種誤差泛函(如差分算法誤差的最小二乘泛函、相速度誤差的最小二乘泛函等)，利用優化方法選擇有限差分系數使誤差泛函極小。優化FDM相比于傳統FDM的優勢在于在不額外增加計算資源的條件下能有效壓制數值頻散。但在復雜模型中速度結構存在強烈間斷或模型存在固液耦合時，優化FDM表現出較大數值誤差[22]。壓制FDM數值頻散的另一種有效方法是利用網格點上地震波場值及其梯度共同近似空間微分算子，在這方面楊頂輝等[23-24]開展了卓有成效的研究，他們系統地發展了地震波模擬的近似解析離散方法。近似解析離散方法是FDM的推廣，該方法能有效壓制數值頻散緣于在近似空間微分算子時近似解析離散方法不僅利用了波場信息，也利用了波場梯度信息。相比于FDM，由于差分算子包含更多的波場信息，使該方法能在粗網格條件下有效壓制數值頻散。然后，與FDM類似，近似解析離散方法也較難使用非規則網格。

基于傅里葉變換的偽譜法(Pseudo-Spectral Method，PSM)可在較粗網格上模擬地震波傳播[25]， 其理論精度可達到FDM采用無窮階差分格式時的精度[26-28]。雖然PSM相比于FDM能在粗網格條件下有效壓制數值頻散[24]，但PSM仍存在網格不靈活、不易處理邊界條件以及難以適用于固液耦合模型等不足。

基于變分原理的有限元法(Finite-Element Method，FEM)可采用任意非規則網格用于離散復雜地質模型[29-31]，因此FEM能有效模擬復雜模型中地震波傳播。除了模型適應性強以外，在FEM中直接令地震波動方程弱形式所包含的邊界積分項為零就可滿足自由地表邊界條件。數值算例表明FEM能較精確近似自由地表邊界條件，因此FEM能較精確模擬地震面波[15]。除以上優點外，FEM還可通過構造邊界條件的方式方便實現固液耦合模型中地震波數值模擬[30]。

FEM在模擬地震波傳播時具有優勢，也具有明顯缺點。傳統FEM的缺點主要表現在計算量大和存儲量大這兩個方面。離散地震波動方程所得到的質量矩陣非對角性以及單元剛度矩陣稠密性決定了FEM的計算量較大[32]。較大的計算量阻礙了FEM在大尺度乃至全球尺度地震波傳播模擬中的應用[33]。為了減少FEM的計算量，有效的方案是構造對角質量矩陣[34-37]，從而避免求解大型線性方程組。構造對角質量矩陣主要有兩種方案：第一種方案是將單元質量按一定比例分配至質量矩陣的對角線上，而將對角線之外的元素直接置為零[34-36]；第二種方案是按照一定規則在離散單元內選擇插值點使質量矩陣對角線上多項式的積分大于零，而對角線之外的積分為零[29，37]。第一種質量集中方案是在假設單元內質點間慣性力解耦的情況下直接對FEM的質量矩陣對角化，該處理方式必然在一定程度上造成FEM精度損失，但理論分析和實際算例都表明集中質量FEM與常規FEM具有相同的收斂階[36，38]，精度損失并不明顯。第二種質量集中方案可在不損失FEM精度的條件下獲得對角質量矩陣，但該質量集中方式需額外引入插值點，這必然增加FEM自由度。相比于傳統FEM，雖然第二種方案減少了計算量，但比第一種方案引入了額外計算量[37]。FEM將地震波動方程所包含的空間微分算子離散之后，在任意離散點處空間微分算子的近似需利用到周圍所有離散點信息，這就決定了FEM單元剛度矩陣的稠密性及在矩陣與向量乘積運算時計算量大的特點。單元剛度矩陣的稠密性也導致FEM需要較大的存儲空間。雖然采用FEM核矩陣的思路可有效地減少FEM的存儲量[16，39]，但存儲核矩陣的方案需要重建單元剛度矩陣而引入了額外計算量。

相對而言，譜元法(Spectral-Element Method，SEM)是一種較為高效的地震波模擬方法[40-42]，該方法結合了FEM的靈活性與譜方法的高精度性，且兼具計算量和存儲量小(與FDM相當)的優點[42]。SEM的靈活性得益于其采用與FEM相同的變分原理以及單元上利用分片插值方法； SEM的高精度性主要源于其借鑒了譜方法的思想而選擇正交插值多項式； SEM計算量和存儲量小的原因在于插值點與積分點重合，從而質量矩陣為對角矩陣、單元剛度矩陣為稀疏矩陣[42]。SEM以其顯著的優點已被運用至不同尺度地震波模擬與成像[2-3，43-44]。值得指出的是，SEM的質量矩陣與剛度矩陣包含的積分項中多項式的最高階數超過2N(其中N為空間插值階數)，在利用GLL(Gauss-Lobatto-Legendre)數值求積時數值積分的精度只有2N-1階，因此GLL數值求積無法精確求取多項式積分。從而積分精度不足而影響SEM的數值精度。為了提升SEM地震波數值模擬精度，本文構造關于GLL數值積分的目標函數，通過優化算法選取GLL數值積分權系數，以提高GLL數值積分精度，最終提升SEM地震波模擬精度，得到基于優化數值積分的譜元法(簡稱為優化譜元法(Optimal SEM，OSEM))。

1 譜元法基本原理

選取各向同性彈性波動方程簡要介紹傳統譜元法(Conventional SEM，CSEM)。2D非均勻各向同性彈性波方程可表示為

(1)

式中：u為位移向量；T為應力張量；c為二階彈性張量；ρ為介質質量密度；f為震源項； ?為梯度算子。用任意測試向量v乘以式(1)的第一式，在2D空間積分，利用分部積分和格林公式可得

(2)

(3)

式中： (ξ,η)為等參單元坐標；Na(ξ,η)為關于第a個差值節點的形函數。由式(3)可將式(2)所包含的物理域上的積分變換至等參域上的積分，即有

(4)

為了計算式(4)，需離散各個等參單元，然后通過離散點上的u、v函數值構造多項式用于近似式(4)，在計算式(4)中積分時，采用GLL數值求積[36]。在等參單元中每個方向上的離散點的坐標是通過求解以下等式的零點而得到

(1-λ)P′N(λ)=0

(5)

式中PN(λ)為N階Legendre多項式，其上標“′”表示一階導數。通過單元上的離散點，可構造Lagrange多項式，利用Lagrange多項式可對u、v近似，即

(6)

式中：L為Lagrange多項式；uije、vije分別表示位移向量和測試向量在第e單元內第(i，j)標號離散點上的離散值。利用式(6-1)及式(3)可計算應力張量T。由已計算的T及式(3)和式(6)，式(4)的離散形式為

(7)

式中：Me為單元質量矩陣；Ke為單元剛度矩陣；Fe為單元上震源項的離散形式；Pe為離散應變張量與雅克比矩陣元素和坐標變換因子組合而成的向量。式(7)中矩陣和向量(包括震源項的離散形式)的具體形式已在文獻[42]中給出，本文不贅述。

值得指出的是，得益于插值點與積分點重合，Ke為稀疏矩陣，因此SEM的計算量和存儲量較小，幾乎與FDM處于同一量級[42]。

2 基于優化數值積分的譜元法(OSEM)

GLL數值求積的精度是2N-1階，而式(4)中被積多項式的最高階數超過2N階，由于GLL數值積分無法對式(4)中積分項精確估計，這必然損失SEM精度。在保證積分點與插值點重合的條件下，本文對GLL數值求積權優化，從而提升GLL數值求積精度，進而提高SEM的精度。

構造如下目標函數

(8)

式中：λj為式(5)的零點；wj為待優化的積分權。為了保證能量守恒和質量守恒，積分權需滿足

(9)

數值求積公式的積分權在對稱條件下能對任意奇數次冪函數λ2i+1精確求積，因此積分權還需滿足

wi=wN-i

(10)

利用共軛梯度算法(Conjugate Gradient)[45-46]求解帶約束條件(式(9)和式(10))的式(8)，可得到優化權系數(表1)?？紤]到對稱性(式(10))，表1只給出積分權的前半部分。由式(9)和式(10)可知，當N=1時積分權系數無法優化，故在表1中未列出。

表1 GLL數值求積優化權系數

3 數值頻散與穩定性分析

利用平面波分析方法[36，47]，可定量給出OSEM的數值頻散與穩定性。一方面，在相同網格采樣條件下S波的數值頻散往往大于P波[16]； 另一方面，S波波長小于P波。因此在空間網格一定的條件下，S波的采樣率(一個波長內采樣點數)大于P波。綜合這兩方面原因，在數值頻散分析時常只需分析S波的數值頻散。假設式(1)的平面波簡諧解為

(11)

式中： 下標j表示空間任意網格節點編號；Aj、Bj為任意常數；k為波數；ω為圓頻率。假設無限空間被離散成無窮多個正方形單元[47]，離散式(7)中的時間微分采用的是二階中心差分[48]，將式(11)代入式(7)，有關數值頻散的分析可退化成求解以下一般特征值問題

(12)

(13)

通過式(13)可計算出S波數值頻散結果(圖1)。由于SEM數值頻散與波傳播方向幾乎無關[47，50]，圖1展示q=0.25(保證2～9階時數值格式穩定)時地震波沿與水平方向成30°方向上傳播的數值頻散。值得指出的是，圖1中數值速度大于物理速度，這種現象由空間離散和時間離散共同決定。

從圖1還可看出，OSEM與傳統SEM總體一致，都表現出弱頻散特點[44]。圖1a顯示N=2時OSEM的數值頻散明顯小于CSEM，當空間插值階數增加時，OSEM的數值頻散與CSEM差別變小。該現象的主要原因在于，隨著插值階數的增加，優化式(8)對數值積分的貢獻減弱。從圖1的局部放大圖可見OSEM的數值頻散小于CSEM，表明本文積分優化策略能提升SEM地震波模擬的精度。

由式(13)可知，為了保證OSEM迭代穩定，式(12)中的特征值應滿足

(14)

對于所有特征值，在式(14)的限定條件下，q能取到的最大值即為OSEM的穩定限。根據式(12)和式(14)，得到的OSEM穩定限如表2所示。便于對比，表2也列出了CSEM的穩定限?？梢姰敳逯惦A數為2或3時，OSEM的穩定限略小于CSEM，兩者的相對差別僅約為5/10000。當插值階數大于3時，OSEM與CSEM的穩定限幾乎無差別。

圖1 OSEM模擬地震波傳播S波的數值頻散(a)～(h)對應于空間插值階數分別為2～9時的數值頻散曲線； 局部放大圖展示OSEM和CSEM局部頻散特征； S波傳播方向與水平方向成30°，庫朗數為0.25

表2 空間插值階數為2～9時OSEM和CSEM的穩定限

4 數值算例

4.1 均勻模型

選擇如圖2所示的均勻介質模型用于檢驗OSEM地震波模擬的精度。模型尺寸為2000m×1000m，離散為邊長為10m的正方形單元。介質的P波、S波速度、密度參數如圖中所示。一個垂直集中力源(同時激發縱波和橫波)位于(1000m，50m)，其震源時間函數是主頻為20Hz的Ricker子波； 兩個接收器R1和R2分別位于(1500m，0)和(1800m，0)處。為了防止邊界截斷而引起的虛假反射，在模型的左、右和下邊界選用二階位移形式的PML吸收邊界條件(厚度為200m)接收反射波[16]，模型上界面為自由地表。模擬時長為1s。為了減少時間離散誤差，時間步長設定為0.05ms。

從表3可知OSEM和CSEM都能較精確模擬地震波傳播，當插值階數大于2時，OSEM和CSEM的最大相對誤差均小于1%。雖然階數為4～8時，OSEM計算得到的R1處相對誤差的最大值大于CSEM，但兩者的相對差別也只在1.0e-5左右，這種較小差別可能是由計算機浮點誤差或時間迭代誤差累積導致。對于其他情況，OSEM的誤差總體上明顯小于CSEM。

圖2 用于檢驗OSEM的數值精度的均勻介質模型震源位于(1000m，50m)； 兩接收器分列于(1500m，0)和(1800m，0)

表3 OSEM和CSEM模擬均勻介質中地震波傳播的誤差

均勻模型中的數值測試結果表明，相比于CSEM，OSEM總體上能有效提升地震波模擬的計算精度。值得指出的是，OSEM和CSEM的差別僅在積分權重，因此只用修改CSEM代碼中的積分權重值從而獲得OSEM； 在實際計算地震波場時，OSEM的計算時間和內存消耗量與CSEM一致。

4.2 傾斜界面模型

為了檢驗OSEM在具有非規則界面模型中的地震波模擬精度，設計了圖3所示的傾斜界面模型。模型尺寸為4000m×2000m。采用非結構化網格離散該模型，網格單元的平均尺寸為50m。傾斜界面上、下介質的P波和S波速度、密度如圖中所示。一個爆炸震源(只激發縱波)位于(2000m，300m)處，其震源時間函數是主頻為15Hz的Ricker子波。模型左、右和下邊界設置20個單元厚度的PML吸收邊界層，上界面為自由地表邊界。兩個接收器R1和R2分別位于(1000m，0)和(3000m，0)處。OSEM采用7階空間插值，時間離散步長為0.1ms； 參考解由10階CSEM采用0.05ms的時間步長計算得到。

從圖4所示的位移水平分量的波場快照可見：當t=0.42s時，由于模型存在分層界面，在上層介質中除了直達波、地表反射波和轉換波以外，還可看到界面的反射波，在下層介質中可看到透射P波、轉換的S波(圖4a和圖4b)； 當t=0.60s時，由于地表反射波和轉換波到達傾斜界面而進一步透射和反射，以至于圖4c和圖4d中，波場成分更復雜。從該圖還可看出OSEM計算結果(圖4a和圖4c)與CSEM計算結果(圖4b和圖4d)無明顯差別。

圖5顯示的是R1和R2兩處合成波形記錄的水平分量。觀察該圖，除了較強的直達波以外，還可看到較強的反射和多次反射震相。從其中的合成波形也可看出，OSEM計算的波形(藍線)與CSEM計算的波形(紅色虛線)高度一致。定量對比R1與R2兩處最大波形誤差，發現OSEM計算波形的水平分量和垂直分量誤差均小于CSEM的，最大相對誤差降低了約0.3%。說明OSEM算法的精度略高于CSEM。

圖3 傾斜界面模型震源位于(2000m，300m)，接收器R1和R2分列于(1000m，0)和(3000m，0)

圖4 傾斜界面模型中地震波場快照(a)、(b)t=0.42s時的波場快照； (c)、(d)t=0.60s時的波場快照； (a)、(c)OSEM計算得到； (b)、(d)CSEM計算得到

圖5 傾斜界面模型中接收器記錄的波形(a)R1處波形； (b)R2處波形

4.3 起伏界面模型

為了檢驗OSEM在復雜起伏介質中的模擬精度，設計了圖6所示的起伏界面模型。起伏地表的最大高程為1500m，模型水平尺度為10000m，最大深度為5000m。模型內存在由正弦函數刻畫的分界面，且離散為平均單元尺寸為50m的網格，空間插值階數為5階。模型左、右和下邊界設置20個單元厚度的PML吸收邊界層，上界面為自由地表邊界。分界面上、下兩層介質參數同第二個算例。起伏地表之上幾乎等距離分布有299個接收器，為了方便顯示，圖6中展示是稀疏提取后的接收器分布。一個爆炸震源(只激發縱波)位于(5000m，100m)處，其震源時間函數是主頻為15Hz的Ricker子波。

圖6 起伏界面模型起伏地表最大高程為1500m，模型內分界面形狀是正弦函數震源位于(2000m，100m)，接收器幾乎等距分布于起伏地表

模擬結果圖7顯示的是OSEM和CSEM兩種方法合成位移水平分量的共炮記錄。從該共炮記錄中可看出明顯的直達波同相軸。在2s以下，來自模型內起伏界面的反射波顯得較明顯。對比圖7a與圖7b可知，OSEM與CSEM模擬結果幾乎一致。但對比所有接收器的誤差總和，發現OSEM比CSEM的約小3%。顯然，OSEM在復雜模型中的計算精度也高于CSEM。

5 結論與討論

針對傳統SEM數值積分精度上存在的不足，通過構造關于數值積分的誤差函數，采用共軛梯度方法優化數值積分的積分權，最終提升SEM中GLL數值積分的精度，從而提升SEM的地震波模擬精度。數值頻散分析證實了OSEM比CSEM更具有弱頻散的特點。數值算例表明，OSEM精度總體上高于CSEM。值得指出的是，本文發展的OSEM只優化了CSEM的數值求積公式，而未改變CSEM對邊界條件的適應性，因此OSEM與CSEM一樣能方便地對各類邊界條件(如PML吸收邊界條件[30]、透射邊界條件[52]和固液耦合邊界條件[30]等)進行數值離散。

圖7 OSEM(a)和CSEM(b)合成共炮記錄

構造FDM優化格式提升其地震波數值模擬精度較為常見，但優化的FEM和SEM較少。雖然本文通過優化數值積分，總體上將SEM的精度提升約3%，提升幅度較小，但本文貢獻在于提出了可構造優化數值積分格式提升FEM或SEM的數值精度，這為FEM或SEM的優化提供了新思路。