基于優(yōu)化數(shù)值積分的譜元法模擬地震波傳播

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

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

0 引言

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

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

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

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

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

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

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

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

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

1 譜元法基本原理

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

(1)

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

(2)

(3)

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

(4)

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

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

(5)

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

(6)

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

(7)

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

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

2 基于優(yōu)化數(shù)值積分的譜元法(OSEM)

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

構(gòu)造如下目標(biāo)函數(shù)

(8)

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

(9)

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

wi=wN-i

(10)

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

表1 GLL數(shù)值求積優(yōu)化權(quán)系數(shù)

3 數(shù)值頻散與穩(wěn)定性分析

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

(11)

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

(12)

(13)

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

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

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

(14)

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

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

表2 空間插值階數(shù)為2～9時OSEM和CSEM的穩(wěn)定限

4 數(shù)值算例

4.1 均勻模型

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

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

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

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

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

4.2 傾斜界面模型

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

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

圖5顯示的是R1和R2兩處合成波形記錄的水平分量。觀察該圖，除了較強的直達波以外，還可看到較強的反射和多次反射震相。從其中的合成波形也可看出，OSEM計算的波形(藍線)與CSEM計算的波形(紅色虛線)高度一致。定量對比R1與R2兩處最大波形誤差，發(fā)現(xiàn)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在復(fù)雜起伏介質(zhì)中的模擬精度，設(shè)計了圖6所示的起伏界面模型。起伏地表的最大高程為1500m，模型水平尺度為10000m，最大深度為5000m。模型內(nèi)存在由正弦函數(shù)刻畫的分界面，且離散為平均單元尺寸為50m的網(wǎng)格，空間插值階數(shù)為5階。模型左、右和下邊界設(shè)置20個單元厚度的PML吸收邊界層，上界面為自由地表邊界。分界面上、下兩層介質(zhì)參數(shù)同第二個算例。起伏地表之上幾乎等距離分布有299個接收器，為了方便顯示，圖6中展示是稀疏提取后的接收器分布。一個爆炸震源(只激發(fā)縱波)位于(5000m，100m)處，其震源時間函數(shù)是主頻為15Hz的Ricker子波。

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

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

5 結(jié)論與討論

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

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

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