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

1 應(yīng)急管理部國家自然災害防治研究院，北京 100085 2 中國科學院地質(zhì)與地球物理研究所，巖石圈演化國家重點實驗室，北京 100029 3 清華大學數(shù)學科學系，北京 100084 4 寧夏大學數(shù)學統(tǒng)計學院，銀川 750021 5 中國地質(zhì)大學地球物理與空間信息學院，武漢 430074

0 引言

數(shù)值求解地震波運動方程不僅有利于認識復雜介質(zhì)中地震波傳播規(guī)律，同時是全波形地震成像的基礎(chǔ).迄今為止，研究者已開發(fā)出多種算法用于地震波運動方程直接求解，如有限差分法（Finite Difference Method，F(xiàn)DM）（Madariaga，1976；董良國等，2000；Zhang et al.，2012；Zhang and Yao，2013；劉少林等，2013；Liang et al.，2018），偽譜法（Pseudo-Spectral Method，PSM）（Wang and Takenaka，2011；杜啟振等，2015；Li et al.，2021；Zou and Cheng，2018），有限元法（Finite Element Method，F(xiàn)EM）（Smith，1975；張懷等，2009）,譜元法（Spectral Element Method，SEM）（Komatitsch and Vilotte，1998；Komatitsch and Tromp，1999，2002；Yan et al.，2009；劉玉柱等，2020；劉少林等，2021；Shen et al.，2022）.FEM和SEM以其網(wǎng)格的靈活性以及能自動滿足自由地表邊界條件等優(yōu)點而引起眾多研究著的關(guān)注，并逐漸被應(yīng)用于實際地震波的正演模擬與成像之中（Wei et al.，2018；Lei et al.，2020）.

FEM采用變分原理將地震波運動方程轉(zhuǎn)化為積分形式（Marfurt，1984；Cohen，2002），并將求解域采用任意網(wǎng)格離散，最后在離散單元上對地震波方程數(shù)值離散并求解.不難發(fā)現(xiàn)FEM具有較好的網(wǎng)格靈活性，可適應(yīng)任意復雜地質(zhì)模型，同時可將FEM得到的地震波方程弱形式所包含的邊界積分項設(shè)置為零，即可滿足自由地表邊界條件.雖然FEM具有以上顯著優(yōu)點，但也存在明顯不足.其不足主要表現(xiàn)在三方面.首先，F(xiàn)EM求解地震波運動方程時計算量較大.FEM離散地震波運動方程以后得到的質(zhì)量矩陣為非對角矩陣（Marfurt,1984），因此FEM需要求解大型線性方程組，計算量較大.雖然可通過質(zhì)量集中算法將質(zhì)量矩陣對角化（Liu et al.，2014b；劉少林等，2014），集中質(zhì)量矩陣必然造成FEM精度損失（劉少林等，2014）.其次，F(xiàn)EM需要較大的計算機內(nèi)存空間.地震波運動方程離散以后，F(xiàn)EM在每個單元上形成的單元剛度矩陣為稠密矩陣，存儲稠密矩陣的非零元需要較大的計算機內(nèi)存.雖然可采用FEM的存儲改進算法，如核矩陣算法（Liu et al.，2014a）和單元矩陣重組算法（Meng and Fu，2017；蘇波等，2019），能有效改進FEM存儲量大的弊端，但這些算法往往會增加FEM的計算量.最后，F(xiàn)EM精度較低.在使用FEM模擬地震波場時，F(xiàn)EM在單元內(nèi)對求解函數(shù)插值逼近時常采用低階插值，雖然低階插值計算量較小、計算速度快，但低階插值面臨著精度低、數(shù)值頻散明顯等缺點（Fichtner，2011；劉少林等，2014）.然而采用高階插值不僅計算量增加而且Runge現(xiàn)象明顯.FEM所具有的缺點，阻礙了該方法在地震波模擬領(lǐng)域的廣泛應(yīng)用（Bielak et al.，2005）.

SEM是一種極具潛力的地震波場數(shù)值模擬方法，該方法兼顧有限元法網(wǎng)格的靈活性與譜方法的高精度性.SEM最早在流體力學領(lǐng)域提出并得到發(fā)展和應(yīng)用（Patera，1984），流體動力學數(shù)值模擬結(jié)果表明SEM具有精度高、快速收斂等優(yōu)點.最早將SEM引入至地震學領(lǐng)域得益于Priolo和Seriani（1991）的研究.Priolo和Seriani（1991）利用SEM實現(xiàn)了二維聲波的數(shù)值模擬，研究表明采用高階SEM模擬聲波時，SEM一個波長內(nèi)的采樣點數(shù)為7即可獲得高精度模擬結(jié)果.隨后，Seriani（Seriani et al.,1992;Seriani,1998）將SEM推廣至彈性波的數(shù)值模擬.值得指出的是，早期的SEM在離散單元內(nèi)采用Chebyshev多項式插值（Priolo and Seriani，1991）（記該類SEM為Chebyshev-SEM），但由于Chebyshev多項式帶權(quán)正交（Wang et al.，2007），Chebyshev-SEM無法形成對角質(zhì)量矩陣，在地震波場模擬時，在每一個時間步更新過程中都需要通過迭代方式求解大型線性方程組（Seriani，1997），因此計算量較大.雖然可采用矩陣張量積算子將稠密的單元剛度矩陣分解成若干個子矩陣，并將單元剛度矩陣與解向量的乘積分解為子矩陣與解向量的乘積，該矩陣分解算法能有效提升Chebyshev-SEM的并行計算效率（Seriani，1997），但并不能降低其計算量，在模擬三維地震波場時計算效率依然較低（Seriani，1998）.

SEM除了可采用Chebyshev正交多項式以外，還可以采用Legendre正交多項式（記該類SEM為Legendre-SEM）.Legendre-SEM相比于Chebyshev-SEM具有三方面的優(yōu)點.首先，Legendre-SEM計算量較小.Legendre-SEM采用GLL點作為多項式的插值點，同時GLL點也為數(shù)值求積的積分點，由于積分點與插值點重合的特性，單元質(zhì)量矩陣為對角矩陣（Komatitsch and Vilotte，1998），對角質(zhì)量矩陣避免了迭代求解線性方程組，大幅減少了計算量.此外，插值點與數(shù)值積分點重合的特性可使得Legendre-SEM單元剛度矩陣具有稀疏的特點（劉少林等，2021），減少了矩陣與解向量乘積計算過程中的浮點運算次數(shù).分析表明Legendre-SEM的浮點運算量與FDM相當（劉少林等，2021）.其次，Legendre-SEM的存儲量較小.將Legendre-SEM的剛度矩陣采用矩陣張量積表示以后（劉少林等，2021），不難發(fā)現(xiàn)Legendre-SEM只用存儲GLL點處雅克比矩陣及其行列式的值，無需存儲稠密單元剛度矩陣，因此存儲量較小.事實上，在相同網(wǎng)格點規(guī)模條件下Legendre-SEM采用任意非規(guī)則網(wǎng)格所需要的存儲量與FDM相當（Komatitsch and Tromp，2002）.最后，Legendre-SEM單元內(nèi)插值所使用的GLL點實際上是Fekete點（Fichtner，2011），F(xiàn)ekete點能保證多項式插值具有最優(yōu)特性，插值多項式能對求解函數(shù)最優(yōu)逼近.由于Legendre-SEM具有以上顯著優(yōu)點，該方法已逐漸成為一種重要的地震波模擬工具（周紅，2018）.

少數(shù)研究者試圖發(fā)展SEM的優(yōu)化算法.例如Seriani和Oliveira（2007）將集中質(zhì)量矩陣與一致質(zhì)量矩陣優(yōu)化加權(quán)組合用于降低Chebyshev-SEM的數(shù)值頻散.然而該質(zhì)量矩陣加權(quán)優(yōu)化策略只在規(guī)則正方形單元網(wǎng)格中有效，而無法推廣至一般非規(guī)則網(wǎng)格.由于Legendre-SEM在正交多項式和插值點選取上保證了其具有較高的數(shù)值精度，以至于文獻中較少報道關(guān)于Legendre-SEM的優(yōu)化格式.

值得指出的是，Legendre-SEM質(zhì)量矩陣所包含的積分項中多項式的階數(shù)不低于2N階（其中N為單元插值多項式的階數(shù)），而GLL數(shù)值求積的精度為2N-1階，由于數(shù)值求積精度不足導致質(zhì)量矩陣元素存在誤差，質(zhì)量矩陣所包含的誤差在一定程度上會影響Legendre-SEM數(shù)值模擬精度.為了提升Legendre-SEM精度，需要對質(zhì)量矩陣改進.借鑒有限差分法優(yōu)化的思想（Zhang and Yao，2013；Wang et al.，2014；Liu,2020），構(gòu)建關(guān)于質(zhì)量矩陣的目標函數(shù)，將質(zhì)量矩陣優(yōu)化，減小質(zhì)量矩陣積分項數(shù)值近似而代入的誤差.通過數(shù)值頻散分析和數(shù)值算例證實了本文提出的質(zhì)量矩陣優(yōu)化方法能提升Legendre-SEM地震波數(shù)值模擬精度.

1 理論與方法

為了論述方便記優(yōu)化質(zhì)量矩陣Legendre-SEM為優(yōu)化譜元法（Optimal spectral-element method，OSEM），記常規(guī)的Legendre-SEM為傳統(tǒng)譜元法（Conventional spectral-element method，CSEM）.

1.1 CSEM基本原理

為了本文論述的完整性，以2D彈性波運動方程為例簡要介紹CSEM的基本原理，更詳細的原理介紹請見Fichtner（2011）和劉少林等（2021）.2D非均勻各向同性介質(zhì)中的彈性波運動方程可表示為

（1）

（2）

其中n為空間Ω邊界處的單位法向量.根據(jù)自由地表邊界條件和無窮遠處自由邊界條件，公式（2）中包含的邊界積分項為0，因此地震波方程弱形式可簡化為：

（3）

（4）

其中Ni（ξ1,ξ2）為關(guān)于單元控制點的形函數(shù).通過求解方程（1-λ2）P′N（λ）=0（其中PN（λ）為N階Legendre多項式）得到N+1個零點，通過這N+1個零點在每個四邊形單元內(nèi)可確定（N+1）2個插值點（等參坐標ξ1和ξ2每個坐標方向有N+1個插值點）.根據(jù)插值點以及插值點上u和v的離散值，在每個單元上的u和v可由以下公式近似：

（5）

（6）

（7）

（8）

（9）

（10）

其中K1和K2的具體表述形式如下：

（11）

其中?為矩陣運算的張量積算子（Seriani，1997）.假設(shè)震源為點源，公式（1）中震源項可由矩張量表示：

（12）

其中M表示矩張量，xs表示震源位置，且震源位于標號為es的單元中，S（t）為震源時間函數(shù).將（12）式以及公式（5）中ve的表示形式代入至（3）式等號右端第二項，關(guān)于震源項的離散形式有以下公式：

（13）

其中Fe具有以下形式：

（14）

其中Gij具有以下形式：

（15）

綜合公式（6）、（9）和（13），公式（3）的離散形式為

（16）

由于測試向量v的任意性，公式（16）可進一步簡化為

（17）

1.2 優(yōu)化質(zhì)量矩陣

從公式（7）和（A4）可知，質(zhì)量矩陣所包含的積分項中被積多項式的階數(shù)為2N+1階（在2D情況下，單元的控制點數(shù)為4個），雖然由于GLL積分點和積分權(quán)的對稱性，對于GLL數(shù)值求積公式能對任意奇數(shù)階多項式精確求積（積分值為0），但對N+1點GLL數(shù)值求積公式只能對階數(shù)不超過2N-2的偶數(shù)階多項式精確求積，而無法準確估計2N階多項式的積分值.由公式（7）可知，可通過調(diào)整數(shù)值積分權(quán)ωα從而使得質(zhì)量矩陣的近似誤差最小化.基于該考慮，構(gòu)造以下目標函數(shù):

（18）

其中變量c1和c2反映單元網(wǎng)格的變形程度，其具體形式見附錄A.總體而言c1和c2的變化范圍為[-0.2,0.2]（附錄A）.為了保證質(zhì)量守恒和能量守恒（Zhang and Yao,2013），需要對積分權(quán)系數(shù)增加以下限定條件：

（19）

圖1 利用CGM得到優(yōu)化系數(shù)的流程圖Fig.1 Flow chart for illustrating coefficient optimization using CGM

為了顯示方便，表1只給出了優(yōu)化權(quán)系數(shù)的左半部分.應(yīng)該說明的是，當N=1時，在公式（19）的限定下，ωi無法進一步優(yōu)化，因此表1不再列出N=1時的優(yōu)化系數(shù).為了對比優(yōu)化積分權(quán)系數(shù)與常規(guī)GLL積分權(quán)系數(shù)的區(qū)別，表2給出了常規(guī)GLL積分權(quán)系數(shù).從表1和表2對比可知，優(yōu)化后的權(quán)系數(shù)與標準積分權(quán)相差較小，總體而言兩者的相對差別在1/1000以內(nèi).兩者的差別隨著N的增加而變小，產(chǎn)生這種現(xiàn)象的原因在于SEM采用高階多項式插值時標準的GLL數(shù)值求積得到的質(zhì)量矩陣離散誤差（公式（18））相比于低階時要小，因此高階插值時積分權(quán)系數(shù)無需過多調(diào)整即可使得公式（18）達到極小.

應(yīng)該指出的是，第一，本文的質(zhì)量矩陣優(yōu)化策略不會增加或減少SEM的計算量.本文構(gòu)造目標函數(shù)（18）式，通過優(yōu)化選取數(shù)值積分權(quán)系數(shù)，從而減小質(zhì)量矩陣離散誤差，該過程并未改變數(shù)值積分點和插值點的位置，數(shù)值積分點和插值點仍保持一致，因此質(zhì)量矩陣仍具有對稱性（公式（7））.第二，公式（11）中K矩陣包含的積分權(quán)仍采用標準的GLL數(shù)值求積積分權(quán).從K矩陣的推導過程可知（劉少林等，2021），由于空間梯度算子的作用（公式（1）和（2）），使得K矩陣所涉及的部分被積函數(shù)的階數(shù)低于質(zhì)量矩陣中被積函數(shù)的階數(shù)，以至于標準的GLL數(shù)值求積在用于K矩陣的形成過程中引入誤差相對較小，因此，本文不再討論K矩陣的優(yōu)化.第三，由于公式（18）不包含模型速度，因此優(yōu)化系數(shù)不依賴于模型速度結(jié)構(gòu)，當優(yōu)化系數(shù)選定以后，OSEM可以用于不同模型的地震波場模擬，而無需再次進行系數(shù)優(yōu)化.

表1 GLL數(shù)值求積的優(yōu)化權(quán)系數(shù)Table 1 Optimized weights for GLL quadrature rule

表2 GLL數(shù)值求積的標準權(quán)系數(shù)Table 2 Standard weights for the GLL quadrature rule

2 理論分析

本節(jié)分析OSEM的數(shù)值頻散特征和數(shù)值穩(wěn)定性.根據(jù)De Basabe和Sen（2007）的討論，可通過求解矩陣的一般特征值問題分析OSEM的數(shù)值頻散與穩(wěn)定性分析.本節(jié)只給出主要分析過程和結(jié)果，更加詳細的原理推導過程請見文獻De Basabe和Sen（2007）和劉少林等（2014）.

在理論分析中，空間網(wǎng)格由規(guī)則的正方形單元組成（Liu et al.,2017a），采用二階中心差分近似公式（17）中的時間微分（Liu et al.，2017b）.假設(shè)地震波運動方程（公式（1））的解可寫成如下簡諧波的形式：

（20）

（21）

（22）

SEM的數(shù)值頻散對地震波的傳播角度和波速比并不敏感（De Basabe and Sen，2007），為了討論方便，選取波數(shù)與x1坐標軸之間的夾角為30°和波速比r為1.6.為了減小時間離散誤差對結(jié)果的干擾，選擇較小的時間步長，使得q=0.005.OSEM和CSEM的數(shù)值頻散如圖2所示.圖2顯示，OSEM與CSEM的數(shù)值頻散總體相似，都表現(xiàn)出SEM弱頻散的特點.當N=2、3、5、7、8和9時，OSEM和CSEM的數(shù)值頻散曲線幾乎重合，為了對比需要，將頻散曲線局部放大.從圖2（a）、（b）、（d）、（f）、（g）和（h）中局部放大圖可看出，OSEM計算的S波數(shù)值速度與物理速度之比相比于CSEM更靠近1，表明OSEM比CSEM具有更弱的數(shù)值頻散.當N=4和6時（圖2（c）和（e）），OSEM的數(shù)值頻散明顯小于CSEM.當N=4和6時，數(shù)值頻散明顯得到改善的可能原因是這些階對應(yīng)的優(yōu)化系數(shù)相比與其他階的優(yōu)化系數(shù)更加接近全局最優(yōu)系數(shù).

從公式（22）不難看出，數(shù)值算法穩(wěn)定性條件為

（23）

3 數(shù)值算例

本節(jié)通過三個數(shù)值算例說明OSEM在地震波場模擬時的高精度性.為了盡量減弱人工邊界處的虛假反射，在邊界截斷處加入約10個波長（地震波主頻對應(yīng)的P波波長）厚度的PML吸收層（Komatitsch and Tromp，2003；Liu et al.，2014a；Xie et al.，2016；謝志南和章旭斌，2017；謝志南等，2019）.為了減小時間誤差的干擾,數(shù)值模擬時采用遠小于穩(wěn)定限的時間步長.事實上，CSEM的高精度性已被大量文獻所證實（Komatitsch and Vilotte，1998；De Basabe and Sen，2007；Liu et al.，2014b），為了定量討論OSEM精度，本文選擇10階的CSEM作為參考.

3.1 均勻模型

選擇如圖3所示的各向同性均勻介質(zhì)模型檢驗OSEM的計算精度.模型大小為2000 m×1000 m.模型介質(zhì)的P波速度為2000 m·s-1，S波速度為1250 m·s-1，密度為2000 kg·m-3.一個垂直方向的集中力源位于（1000 m,50 m）處，其震源時間函數(shù)為主頻為10 Hz的Ricker子波.兩個接收器R1和R2分別位于（1500 m,0 m）和（1800 m,0 m）處.本次模擬空間網(wǎng)格單元尺寸為10 m，時間步長為0.05 ms.為了定量評估誤差，定義接收點的最大誤差max=‖Unum-Ur‖∞，其中U表示由歸一化位移水平或垂直分量的時間序列組成的向量，上標num和r分別表示接收器處數(shù)值解與參考解.表3給出了OSEM計算得到的位移垂直分量最大誤差（由于水平分量結(jié)果類似，因此不再討論）.為了對比需要，表3也給出了CSEM的計算結(jié)果.

表3給出了當空間插值為2-9階時OSEM和CSEM的最大誤差.從表3可知OSEM和CSEM模擬均勻介質(zhì)中地震波傳播時均具有高精度的特點.當階數(shù)為2時，OSEM和CSEM的相對誤差均在2%以內(nèi).隨著階數(shù)的增加，OSEM和CSEM的相對誤差均變小，當階數(shù)為9階時，相對誤差減少至0.2%以內(nèi).OSEM和CSEM雖然都能較精確模擬均勻介質(zhì)中彈性波，但OSEM的誤差要小于CSEM.

3.2 分層模型

為檢驗OSEM在分層模型中的計算精度，設(shè)計了如圖4所示的分層模型.模型的大小為4000 m×2000 m，模型內(nèi)部存在兩個分界面.上層界面由一傾斜面構(gòu)成；下層界面由一起伏面構(gòu)成.這兩個分界面將模型分層三層，從上至下分別編號為A、B和C，每一層的參數(shù)如表4所示.采用四邊形單元將模型離散，網(wǎng)格單元的平均尺寸為30 m.一個爆炸震源位于模型中心，其震源時間函數(shù)是主頻為10 Hz的Ricker子波.兩個接收器為R1和R2分別位于（3500 m,0 m）和（4500 m,0 m）處.空間采用8階多項式插值.時間離散步長為0.01 ms.

圖2 OSEM與CSEM數(shù)值頻散曲線紅色線是OSEM;黑色線是CSEM.（a）—（h）分別對應(yīng)的空間插值階數(shù)為2-9階.（a）、（b）、（d）、（f）、（g）和（h）中OSEM和CSEM數(shù)值頻散曲線幾乎重合,為了對比需要,將頻散曲線局部放大.Fig.2 The numerical dispersion curves of OSEM and CSEMThe red curves are calculated from OSEM;the black curves are computed from CSEM.Plots （a）—（h） are generated with spatial interpolations with second- to ninth-order,respectively.When N=2,3,5,7,8 and 9,the numerical dispersions of OSEM and CSEM are nearly consistent.For comparison purposes,portions of these dispersion curves in plots （a）,（b）,（d）,（f）,（g） and （h） are zoomed in.

表3 CSEM和OSEM模擬均勻介質(zhì)中地震波傳播的誤差對比Table 3 A comparison of the maximum errors from CSEM and OSEM when modeling seismic waves in homogeneous media

圖3 均勻介質(zhì)模型用于檢驗OSEM的地震波數(shù)值模擬精度Fig.3 Homogeneous model is employed to investigate the accuracy of OSEM

OSEM模擬結(jié)果如圖5和表5所示.從圖5可看出，總體上看OSEM和CSEM得到的合成波形與參考波形幾乎重合，兩者在模擬層狀介質(zhì)中地震波傳播時都具有高精度性.從圖5可看到，在振幅較強的直達P波之后可看到明顯的透射、反射震相.為了更清楚對比OSEM和CSEM的波形模擬精度，圖6展示了波形誤差.從圖6可看出，總體上OSEM的波形誤差（青色線）小于CSEM的波形誤差（藍色線）.為了定量對比OSEM和CSEM的計算精度，表5給出了OSEM和CSEM在R1和R2處合成波形水平分量和垂直分量的最大誤差.從最大誤差對比可看出，OSEM的精度要明顯高于CSEM.

圖4 分層介質(zhì)模型兩個分層界面將模型分成三層，從上往下分別記為A、B和C，每一層的介質(zhì)參數(shù)如表4所示.圖中五角星表示震源；三角形表示接收器.Fig.4 A model with three layersTwo interfaces divide the model into three layers,which is marked by A,B,and C,respectively.The parameters of each layer are shown in Table 4.The upper interface is an inclined line;the lower interface is an undulating line.

圖5 合成波形與參考波形對比（a）與（c）分別為R1處波形的水平分量與垂直分量；（b）與（d）分別為R2處波形的水平分量與垂直分量.圖中青線和藍線分別由OSEM和CSEM計算得到；紅線為參考波形.Fig.5 A comparison between the synthetic and reference waveformsPlots （a） and （c） respectively show the horizontal and vertical components of waveforms at R1.Plots （b） and （d） respectively show the horizontal and vertical components of waveforms at R2.Cyan lines are generated by OSEM;blue lines are computed by CSEM.Red lines represent reference waveforms.

表4 分層模型介質(zhì)參數(shù)Table 4 Media parameters of the layered model in Fig.3

圖6 OSEM和CSEM得到的合成波形誤差對比（a）和（b）為R1處的合成波形誤差；（c）和（d）為R2處的合成波形誤差.（a）和（c）為水平分量誤差；（c）和（d）為垂直方向誤差.Fig.6 A comparison between errors of synthetic waveformsPlots （a） and （b） are at the receiver R1;plots （c） and （d） are at the receiver R2;Plots （a） and （c） are for the horizontal components;plots （b） and （d） are for the vertical components.Cyan lines are generated by OSEM;blue lines are computed by CSEM.

表5 層狀介質(zhì)中OSEM與CSEM模擬精度對比Table 5 A comparison of the accuracy of OSEM and CSEM in layered model

3.3 起伏模型

設(shè)計如圖7的所示模型用于檢驗OSEM在起伏地表模型中地震波模擬精度.模型水平方向延伸10000 m，垂直向下延伸至地下5000 m處.如圖7所示，模型除了地表具有強烈起伏外，模型內(nèi)部還存在分層界面.模型上下兩層的P波、S波以及密度參數(shù)顯示在圖7中.一個爆炸源位于（5000 m,1000 m）處，其震源時間函數(shù)是主頻為15 Hz的Ricker子波.5個接收器（圖7中黑色三角形）從上往下分別編號為1-5.采用四邊形單元離散模型，四邊形單元的平均尺寸為30 m，該網(wǎng)格尺寸能保證最短波長內(nèi)平均單元采樣數(shù)約為2.OSEM空間插值多項式的階數(shù)為4.時間步長為0.01 ms.模擬結(jié)果如圖8所示.

圖7 地表起伏模型模型地表處起伏強烈，此外模型內(nèi)部具有一傾斜界面.模型介質(zhì)參數(shù)在圖中給出.五角星代表震源；地表處的三角形表示接收器.Fig.7 A model with rugged topographyBeside the topography,an incline interface is inside the model.The parameters of the model are presented in the plot.The star at （5000 m,1000 m） denotes the source;the triangles represent the receivers.

為了對比，圖8也展示了CSEM波形模擬結(jié)果.從模擬結(jié)果可看出，OSEM和CSEM均能準確模擬地表強烈起伏模型中地震波傳播.從合成地震記錄中可看到除了振幅較強的直達波以外，還在直達波之后發(fā)現(xiàn)明顯的反射波，以及可看到來自地表起伏界面的多次反射波.從圖8難以直接看出OSEM和CSEM兩者精度差別.為了定量對比，計算OSEM和CSEM的水平和垂直分量最大誤差，OSEM的水平和垂直分量最大誤差分別為3.881776×10-2和4.219775×10-2，CSEM的水平和垂直分量最大誤差分別為4.048275×10-2和4.297425×10-2.由此對比可以知，OSEM的精度要高于CSEM.

4 討論與結(jié)論

本文提出了一種質(zhì)量矩陣優(yōu)化策略用于減小SEM質(zhì)量矩陣中包含的誤差.本文核心思想在于構(gòu)造數(shù)值積分與真實積分的最小二乘目標函數(shù)，通過優(yōu)化方法選取數(shù)值積分權(quán)用于減少質(zhì)量矩陣誤差.通過理論分析和實際算例證實了該優(yōu)化策略能有效提升SEM地震波模擬的計算精度.本文對SEM的優(yōu)化并不增加計算量，而只用修改SEM質(zhì)量矩陣中所包含的積分權(quán)便可獲得精度提升.應(yīng)該指出的是，隨著插值階數(shù)的提高，OSEM的精度提升越來越小.在實際地震波場模擬時，為了在計算效率與精度之間達到平衡，往往空間插值階數(shù)不超過5階.因此，當采用OSEM用于實際地震波場模擬時，建議空間采用2-5階插值.

應(yīng)該指出的是，本文提出的優(yōu)化策略并未破壞SEM正交多項式的正交性，因此SEM的高精度性以及指數(shù)收斂性仍得到了保持.實事上，傳統(tǒng)的SEM通過選取正交多項式和Fekete點用于多項式插值和逼近待求解函數(shù)，使其具有較高的數(shù)值精度.本文的貢獻在于，可通過優(yōu)化策略將SEM的精度進一步提升.本文在優(yōu)化目標函數(shù)的過程中，只采用CGM用于系數(shù)的優(yōu)化，以至于只能得到局部最優(yōu)的優(yōu)化系數(shù).如果采用全局優(yōu)化方法（如模擬退化方法、蒙特卡洛方法），有可能獲得全局最優(yōu)的優(yōu)化系數(shù)，從而有望進一步提升SEM精度.本文提出的質(zhì)量矩陣優(yōu)化策略可方便推廣至三維SEM，通過構(gòu)造三維情況下關(guān)于質(zhì)量矩陣的目標函數(shù)（與公式（18）類似），可得到適用于三維SEM的優(yōu)化權(quán)系數(shù)和優(yōu)化質(zhì)量矩陣.

圖8 接收器記錄的波形對比（a）和（b）分別為水平方向位移和垂直方向位移；青色線、藍色線和紅色線分別為OSEM、CSEM和參考解得到的合成波形.Fig.8 Comparison of synthetic waveformsPlots （a） and （b） show horizontal and vertical components of displacements,respectively;Cyan,blue and red lines are generated by OSEM,CSEM and reference solution,respectively.

致謝感謝兩名匿名審稿專家對本文提出的建設(shè)性意見.感謝與中國科學院地質(zhì)與地球物理研究所劉有山副研究員的討論.本文的計算模擬在中國地質(zhì)大學（武漢）地球物理與空間信息學院的TS10000服務(wù)器上完成.

附錄A

（A1）

由公式（A1），雅克比矩陣可表示為

（A2）

其中a1-a4和b1-b2的表達形式如下：

（A3）

由公式（A2）可知雅克比矩陣行列式的值可表示為

附圖A1 四種典型四邊形單元用于討論c1和c2（a）—（d）對應(yīng)的（c1,c2）值分別為（0,0）、（0.2,0.2）、（-0.2,-0.2）和（-0.2,0.2）.Fig.A1 Four typical quadrilateral elements to discuss the values of c1 and c2The elements in plots （a）—（d） are with （c1,c2） of （0,0）,（0.2,0.2）,（-0.2,-0.2） and （-0.2,0.2）,respectively.