邢浩潔 劉愛文 李小軍 陳 蘇 傅 磊

1） 中國北京 100081 中國地震局地球物理研究所

2） 中國北京 100124 北京工業(yè)大學(xué)城市建設(shè)學(xué)部

引言

人工邊界條件是地震波動數(shù)值模擬技術(shù)的關(guān)鍵環(huán)節(jié)（廖振鵬，2002），是計算人工邊界節(jié)點運動的各種方法的統(tǒng)稱，其任務(wù)在于盡可能地“算準”外行波引起的邊界點運動，以達到消除不合理邊界虛假反射的目的.半個世紀以來，國內(nèi)外研究人員提出了基于各種直接或間接方法適用于不同波動問題的人工邊界條件，其中我國學(xué)者廖振鵬及其合作者所發(fā)展的多次透射公式（multi-transmitting formula，縮寫為MTF），本文稱之為廖氏透射人工邊界，具有形式簡單、易于實現(xiàn)、通用性好的優(yōu)點（廖振鵬等，1984；Liao，Wong，1984）.該邊界具有清晰明確的物理機制，即使用統(tǒng)一人工波速進行“誤差波多次透射”，因而受到研究人員和工程師的青睞.

三十多年來，MTF 人工邊界在場地地震反應(yīng)和土-結(jié)相互作用數(shù)值計算中發(fā)揮了重要作用（李小軍等，1995；Shiet al，2016；Huang，2018；陳少林等，2020），關(guān)于MTF 精度優(yōu)化（廖振鵬，李小軍，1995；Liao，1996）及數(shù)值穩(wěn)定性問題（Liao，Liu，1992；李小軍，廖振鵬，1996；廖振鵬等，2002；景立平，2004；章旭斌等，2015；Xie，Zhang，2017）的研究也隨之不斷深入.近幾年，MTF 開始被用于基于譜元法的地震波動模擬（戴志軍等，2015；邢浩潔，李鴻晶，2017a，b；于彥彥等，2017；Xinget al，2021a）.Xing 等（2021a）對二維SH 波動和P?SV 波動的研究結(jié)果表明MTF 邊界的精度略低于完美匹配層（perfectly matched layer，縮寫為PML）邊界（李寧等，2007；Liuet al，2014；謝志南，章旭斌，2017），但要優(yōu)于黏彈性邊界（劉晶波，李彬，2005；杜修力等，2006）.最近，作者（邢浩潔等，2021；Xinget al，2021b）成功地將多個人工波速引入MTF 表達式當中，提出了一種多人工波速優(yōu)化透射邊界，記為caj-MTF，其中caj（j＝1，2，···，N）表示與邊界階次N相對應(yīng)的多個人工波速.該邊界所具有的多個人工波速配置能夠與彈性波、兩相介質(zhì)波動、頻散波動等問題中的多種物理波速一一對應(yīng)（即各個人工波速參數(shù)可以分別取為各個物理波速），解決了傳統(tǒng)廖氏透射邊界中單一人工波速不便于選取的難題，使得復(fù)雜波動模擬中的邊界性能得到顯著提升.

與已有的廖氏透射邊界優(yōu)化形式（李小軍，廖振鵬，1996；Liao，1996）相比，多人工波速透射邊界與MTF 表達式更為接近，完全繼承了MTF 簡潔、精度可控和普遍適用性的優(yōu)點.廖振鵬（2002）曾經(jīng)提出一種含有多個人工波速的MTF 邊界的推廣形式，其解析表達式為

式中，u（s，t）為人工邊界所在的指向內(nèi)域的離散網(wǎng)格線上的波場，N為邊界階次，I為單位算子，B（caj）為表示時間-空間移動的離散算子，caj（j＝1，2，···，N）為一組可調(diào)的人工波速.離散算子B（caj）對波場u（s，t）的作用是將其轉(zhuǎn)換為另一個時空節(jié)點（即參考點）處的波場，該點在時間上向前平移一個時刻Δt，空間上沿著人工邊界節(jié)點所在的離散網(wǎng)格線向計算區(qū)域內(nèi)部平移cajΔt距離；算子的乘積如B（ca1）B（ca2）表示時間-空間平移（2Δt，ca1Δt＋ca2Δt）距離，依此類推.當caj-MTF 邊界中各個人工波速取相同數(shù)值時，即退化為MTF 邊界的算子表達式［I?B（ca） ］Nu（s，t）＝0 （廖振鵬等，2002；景立平，2004）.邢浩潔等（2021）證明了這兩個邊界表達式具有等價性，雖然二者的表達形式和數(shù)值離散格式都不一樣，但在波動數(shù)值模擬中的精度和穩(wěn)定性幾乎完全相同.不過，離散邊界caj-MTF 比解析邊界更易于與不同的內(nèi)域離散方法相結(jié)合，前者只需將各個離散參考點的運動由臨近的內(nèi)域節(jié)點的運動進行插值表示即可.作者在前期工作中對caj-MTF 邊界模擬標量波和彈性波的性能進行了理論分析，并且在有限差分和有限元波動數(shù)值模擬中予以驗證，本文將進一步探討該邊界在譜元法地震波動模擬中的應(yīng)用.

本文將首先介紹多人工波速優(yōu)化透射邊界caj-MTF 的表達式，闡明該邊界提高邊界性能的原理；然后給出在該邊界譜元離散網(wǎng)格下的具體數(shù)值格式，并指出其與有限差分或有限元格式的差異；最后，采用簡單雷克（Ricker）子波激勵下的體波大角度透射和面波透射數(shù)值算例，研究該邊界對不同外行波的反射率，并與傳統(tǒng)MTF 邊界、完美匹配層邊界、黏彈性邊界、一階旁軸近似邊界（Clayton，Engquist，1977）及自由邊界模擬結(jié)果進行比較，以期獲得關(guān)于譜元法地震波動模擬中不同人工邊界基本性能的定量認識.

1 基于多人工波速的優(yōu)化透射邊界

多人工波速優(yōu)化透射邊界（caj-MTF）的具體原理詳見Xing 等（2021b）和邢浩潔等（2021），這里僅簡要介紹.常用的前三階優(yōu)化邊界表達式為

圖1 為caj-MTF 邊界示意圖.每個人工邊界節(jié)點0 的運動都由一系列參考點（圖中空心圓圈）的運動推算得到，這些參考點位于一條指向內(nèi)域的離散網(wǎng)格線上（圖中以局部坐標s表示）且分別處于不同時刻（t－Δt，t－2Δt，t－3Δt，···）.參考點的空間坐標ca1Δt，ca2Δt，ca1Δt＋ca2Δt，···由各個人工波速ca1，ca2，ca3和時間步長Δt共同決定.觀察圖1 并比較上述兩種邊界表達式可知，caj-MTF 邊界對MTF 邊界的優(yōu)化是通過將各個時刻的單一參考點替換為多個參考點來實現(xiàn)的.邢浩潔等（2021）將這種改進解釋為在保留“誤差波多次透射”物理機制的同時對各階誤差波的表述加以優(yōu)化.

圖1 多人工波速優(yōu)化透射邊界（ caj-MTF）示意圖（s0t 為定義人工邊界條件的局部坐標系）Fig. 1 Schematic diagram of optimized transmitting boundary with multiple artificial wave velocities （ caj-MTF）（s0t is a local coordinate system that is used to define the artificial boundary condition）

2 優(yōu)化透射邊界的譜元計算格式

由于caj-MTF 邊界具有與MTF 邊界非常相似的形式，其數(shù)值實現(xiàn)方式與后者也基本一致（邢浩潔，李鴻晶，2017a）.如圖2 所示，利用緊鄰人工邊界的譜單元節(jié)點來實現(xiàn)caj-MTF.對于任意人工邊界節(jié)點0，要求該節(jié)點位于一條指向內(nèi)域的離散網(wǎng)格線上，緊鄰人工邊界的譜單元在該離散網(wǎng)格線上的一組節(jié)點為0，1，···，M（M為單元插值函數(shù)的階次）.在采用caj-MTF 計算（p＋1）Δt時刻邊界節(jié)點0 的運動時，只需將式（3）中pΔt，（p－1）Δt，（p－2）Δt，···各個時刻的邊界參考點（圖中空心圓圈）的運動用譜單元節(jié)點0，1，···，M的運動插值表示即可.

圖2 優(yōu)化透射邊界的譜元計算格式Fig. 2 Spectral-element computational scheme of the optimized transmitting boundary

根據(jù)上述插值方法得到優(yōu)化透射邊界caj-MTF 的譜元計算格式為

式中：u1，u2和u3表示單元節(jié)點運動組成的列向量；T1，T21，···，T33表示插值系數(shù)組成的行向量；sk（k＝0，1，···，M）表示從邊界節(jié)點0 到編號為k的內(nèi)域節(jié)點的距離；T31，T32和T33中的插值系數(shù)與T21，T22類似，為節(jié)省篇幅，不再具體列出.式（5）—（6）所表示的caj-MTF 邊界的譜元計算格式與已有的MTF 邊界的譜元計算格式（戴志軍等，2015；邢浩潔，李鴻晶，2017a，b；于彥彥等，2017；Xinget al，2021a）非常相似，二者都簡單地將各個邊界參考點的運動用該方向上譜單元節(jié)點的運動插值表示.這種簡單實用的施加方式為廖氏透射邊界的重要特色，本文邊界繼承了這個優(yōu)點.不過，本文邊界也具有Xing 等（2021a）所指出的譜元法中MTF 數(shù)值格式的缺點，即它不能像傳統(tǒng)MTF 有限差分或有限元格式那樣由一階邊界格式的乘方來得到高階邊界格式，這使得譜元法中高階透射邊界的精度和穩(wěn)定性受到數(shù)值離散誤差的影響很大.因此，我們建議在譜元法地震波動模擬中使用MTF 或caj-MTF 邊界時，邊界階次取N＝2.

3 波動模擬結(jié)果分析

3.1 計算模型和用于對比的幾種經(jīng)典人工邊界

本文以譜元法模擬脈沖型地震動輸入的波動模擬結(jié)果來檢驗多人工波速優(yōu)化透射邊界（caj-MTF）的實用性，并通過與幾種經(jīng)典人工邊界以及自由邊界（物理邊界）進行比較分析，得到關(guān)于caj-MTF 邊界誤差量級及其相對優(yōu)勢的定量認識.采用如圖3 所示的計算模型，分別考慮體波大角度透射問題和面波透射問題兩種情形.

圖3 人工邊界性能分析模型（a） 體波大角度透射問題；（b） 面波透射問題Fig. 3 Analytical models of artificial boundary performance（a） The problem of body wave transmission over large incident angles；（b） Surface wave transmission problem

對于體波問題，波速設(shè)定為cP＝1 000 m/s，cS＝440 m/s，對應(yīng)于泊松比ν＝0.38.模型尺寸為400 m×1 200 m，波源位于坐標（120 m，600 m），三個觀察點位于左邊界上，A（0，500 m），B（0，400 m）和C（0，200 m）.這樣大的縱橫波速比（或泊松比）和透射角度對邊界性能提出了很高的要求.對于面波問題，波速cP＝1 000 m/s，cS＝577.4 m/s，對應(yīng)的泊松比為ν＝0.25.模型尺寸為800 m×600 m，波源位于坐標（500 m，20 m），三個觀察點位于右邊界上，從地表向下依次為D（800 m，0），E（800 m，20 m）和F（800 m，80 m）.兩個問題的波源均采用雷克子波，施加在節(jié)點的z方向位移uz之上，波源函數(shù)表達式為

式中，t0為脈沖起始時刻，f0為中心頻率，與之對應(yīng)的脈沖時間寬度為2/f0，對應(yīng)的峰值頻率fmax約為2.5f0.本文取t0＝0，f0＝10 Hz，則脈沖時間寬度為0.2 s，峰值頻率約為25 Hz.輸入波的時程及其傅里葉譜如圖4 所示.

圖4 輸入波時程及其傅里葉幅值譜Fig. 4 Time history and Fourier spectrum of the incident wave

模型的空間離散采用25 節(jié)點正方形譜單元，單元尺寸為20 m×20 m.時間步長在滿足內(nèi)域時間積分穩(wěn)定性的前提下，隨著不同邊界條件的影響進行一定調(diào)整.用于與本文caj-MTF 邊界相比較的其它幾種邊界分別為：MTF 邊界（Xinget al，2021a）、PML 邊界（Zenget al，2011；Liuet al，2014）、黏彈性邊界（杜修力等，2006）、一階旁軸近似邊界（Clayton，Engquist，1977）以及自由邊界（物理邊界，譜元法自動滿足該條件）.采用擴展計算區(qū)域的數(shù)值解作為參考解，來衡量各個人工邊界的反射誤差.

3.2 體波大角度透射問題

不同邊界條件模擬體波大角度透射問題所得到的z方向位移uz的波場快照如圖5 所示，各子圖的快照時刻分別為0.3 s，0.4 s，0.5 s，0.7 s.從圖5 可以看出：自由邊界這種物理邊界會造成大量虛假反射，且邊界反射存在波型轉(zhuǎn)換（四種反射模式：P?P，P?S 和S?P，S?S）；所有人工邊界都能夠極大地降低不必要的虛假反射，其模擬結(jié)果遠遠優(yōu)于自由邊界，這體現(xiàn)出人工邊界條件的重要性；PML 邊界精度最高，已看不出邊界反射；本文caj-MTF 邊界的精度明顯好于MTF、黏彈性邊界和一階旁軸近似邊界，雖然不及PML 那樣完美，但也達到了很好的模擬效果.由此可見，波場快照定性地反映了不同人工邊界的性能差異.

圖6 給出了觀察點A，B和C處z方向位移uz的時程及其相對于參考解的誤差曲線.為便于比較，圖中未顯示誤差很大的自由邊界結(jié)果以及與MTF 和黏彈性邊界精度相當?shù)囊浑A旁軸近似邊界結(jié)果，后文表格將列出全部結(jié)果.通過觀察點時程曲線可以看到，黏彈性邊界結(jié)果總是小于參考解，表明該邊界給計算模型提供的約束“偏剛”，另外幾種邊界結(jié)果則看不出明顯規(guī)律.從誤差曲線可知，除PML 之外的幾種人工邊界均產(chǎn)生了顯著的誤差.不過，caj-MTF 邊界相對于MTF 邊界和黏彈性邊界具有明顯優(yōu)勢，不僅其誤差幅值大大小于后兩者，還對大角度透射波（見觀察點B和C的誤差曲線）實現(xiàn)了較高精度的模擬.前文所指出的caj-MTF 邊界能夠較傳統(tǒng)邊界更好地適應(yīng)含有不同物理波速或存在大角度透射（大角度對應(yīng)于大的視傳播速度）的復(fù)雜波動問題，在這里得到驗證.

為定量地比較不同邊界條件對體波模擬效果的差異，表1 詳細列出了觀察點A，B和C處波場分量ux和uz的數(shù)值解相對于參考解（擴展模型解）的誤差曲線（uz的誤差曲線見圖6 右側(cè)，ux的未畫出）的峰值（絕對值），以及這些誤差峰值相對于參考解時程峰值的百分比（對于每種人工邊界條件，將各個觀察點的結(jié)果進行平均）.

圖6 觀察點A （a），B （b）和C （c）處z 方向位移 uz的時程及誤差曲線Fig. 6 Time histories and error curves of the displacement uz at the observation points A （a），B （b） and C （c）

由表1 可知，各個觀察點的最大誤差絕對值存在明顯差異，從A點到B點再到C點，誤差逐步降低，這體現(xiàn)了波傳播過程中幾何衰減所導(dǎo)致的幅值下降.不過，各個觀察點的最大誤差總體上隨著不同邊界條件呈現(xiàn)出一致的變化趨勢，這個趨勢最終體現(xiàn)在邊界的峰值百分誤差上.不同邊界的峰值百分誤差顯示出邊界精度由高到低依次為：PML 邊界最高，caj-MTF 邊界較高，MTF 邊界、黏彈性邊界和一階旁軸近似邊界中等，自由邊界最低.各個峰值百分誤差數(shù)值存在著數(shù)量級的差異，PML 邊界低至3.2%，caj-MTF 邊界為15.3%，MTF 邊界、黏彈性邊界和一階旁軸近似邊界分別為37.0%，62.1%，33.7%，自由邊界則高達119.2%.由于峰值百分誤差是對各個邊界誤差最大值的度量，而從圖5 所示的波場快照來看，總體上人工邊界對波動能量的反射比例要遠小于這些峰值.因此可以得出結(jié)論，在對大縱橫波速比和大角度透射的體波進行模擬時，除PML 邊界具有近乎完美的性能之外，本文的caj-MTF 邊界亦具有很高的精度.

表1 體波情形下A，B 和C 點誤差曲線峰值（絕對值）以及每種邊界條件下誤差峰值與參考解峰值的百分比Table 1 The peaks （absolute value） of error curves of the points A，B，C，and the percentages between those error peaks and reference solution peaks on different boundary conditions：Body wave case

圖5 不同邊界條件模擬體波大角度透射問題的波場快照（時刻依次為0.3 s，0.4 s，0.5 s，0.7 s）（a） caj-MTF 邊界；（b） MTF 邊界；（c） PML 邊界；（d） 黏彈性邊界；（e） 一階旁軸近似邊界；（f） 自由邊界（物理邊界）Fig. 5 Wave field snapshots of body wave propagating over large incident angles obtained from different boundary conditions （time instants are 0.3 s，0.4 s，0.5 s，0.7 s）（a） caj-MTF boundary；（b） MTF boundary；（c） PML boundary；（d） Viscous-spring boundary；（e） The first-order paraxial-approximation boundary；（f） Free boundary （physical boundary）

3.3 面波透射問題

不同邊界條件模擬面波透射問題所得到的z方向位移uz的波場快照如圖7 所示，各子圖的快照時刻分別為0.5 s，0.9 s，1.2 s.面波快照結(jié)果與上述體波情形相似，主要為：與自由邊界存在的大量虛假反射相比，所有人工邊界都能夠消除絕大部分邊界反射；PML 邊界模擬效果近乎完美；本文caj-MTF 僅存在微量面波反射，精度很高；MTF、黏彈性邊界和一階旁軸近似邊界則有著較為顯著的面波反射，精度一般.

圖7 不同邊界條件模擬面波問題的波場快照（時刻依次為0.5 s，0.9 s，1.2 s）（a） caj -MTF 邊界；（b） MTF 邊界；（c） PML 邊界；（d） 黏彈性邊界；（e） 一階旁軸近似邊界；（f） 自由邊界（物理邊界）Fig. 7 Wave field snapshots of surface waves obtained by different boundary conditions（time instants are 0.5 s，0.9 s，1.2 s）（a） caj -MTF boundary；（b） MTF boundary；（c） PML boundary；（d） Viscous-spring boundary；（e） The first-order paraxial-approximation boundary；（f） Free boundary （physical boundary）

圖8 為面波觀察點D，E和F處uz的時程及其相對于參考解的誤差曲線.同樣地，為了便于比較，圖中未給出誤差很大的自由邊界結(jié)果以及與MTF 和黏彈性邊界精度相當?shù)囊浑A旁軸近似邊界結(jié)果，后文表格將全部列出.D點位于地表，E點和F點分別距離地表20 m和80 m.三個觀察點處時程曲線的幅值逐步降低，表現(xiàn)出面波幅值從地表向下逐步衰減的規(guī)律；從誤差曲線來看，三個觀察點的結(jié)果均表明PML 邊界誤差最小，caj-MTF 有少量誤差，MTF 邊界和黏彈性邊界的誤差進一步增大.綜上可見，面波問題的模擬結(jié)果再次證明caj-MTF 邊界相對于MTF 和黏彈性邊界具有非常明顯的優(yōu)勢.

為定量比較不同邊界條件對面波模擬效果的差異，表2 詳細列出了觀察點D，E和F處波場分量ux和uz的數(shù)值解相對于參考解（擴展模型解）的誤差曲線（uz的誤差曲線見圖8 右側(cè)，ux的未畫出）的峰值（絕對值），以及這些誤差峰值相對于參考解時程峰值的百分比（對于每種人工邊界條件，將各個觀察點的結(jié)果進行平均）.

圖8 觀察點D （a），E （b）和F （c）處z 方向位移 uz的時程及誤差曲線Fig. 8 Time history and error curves of the displacement uz at the observation points D （a），E （b） and F （c）

由表2 可知，各個觀察點的最大誤差總體上隨著不同邊界條件呈現(xiàn)出一致的變化趨勢，這種趨勢最終體現(xiàn)在峰值百分誤差上：PML 邊界誤差最小，caj-MTF 邊界次之，MTF 邊界、黏彈性邊界、一階旁軸近似邊界誤差中等，自由邊界則誤差極大.不同邊界的面波反射誤差存在數(shù)量級差異，PML 邊界為4.8%，caj-MTF 邊界為11.7%，MTF 邊界、黏彈性邊界、一階旁軸近似邊界分別為59.7%，30.4%，41.0%，自由邊界高達425.4%.因此，面波模擬結(jié)果同樣表明caj-MTF 邊界是一種精度很高的人工邊界，雖然未達到PML 邊界那樣近乎完美的精度，但是顯著優(yōu)于MTF 邊界、黏彈性邊界和一階旁軸近似邊界.

表2 面波情形下D，E 和F 點誤差曲線峰值（絕對值）以及每種邊界條件下誤差峰值與參考解峰值的百分比Table 2 The peaks （absolute value） of error curves of the points D， E，F(xiàn)， and the percentages between those error peaks and reference solution peaks under different boundary conditions：Surface wave case

4 討論與結(jié)論

在使用有限的數(shù)值離散模型研究無限域波動問題時，如何施加人工邊界條件（即吸收邊界條件）是尚未得到圓滿解決的關(guān)鍵環(huán)節(jié).本文將作者（邢浩潔等，2021；Xinget al，2021b）最近發(fā)展的多人工波速優(yōu)化透射邊界（caj-MTF）應(yīng)用于高精度譜元法的地震波動模擬，并通過與現(xiàn)有的廖氏透射邊界（MTF）、完美匹配層邊界（PML）、黏彈性邊界以及一階旁軸近似邊界進行比較分析，證明了這是一種便捷、高效的人工邊界條件實現(xiàn)方法.

caj-MTF 邊界對傳統(tǒng)MTF 邊界的計算波速參數(shù)進行優(yōu)化，但保持了后者“誤差波多次透射”的基本原理以及簡單通用的表達形式，因此caj-MTF 邊界像MTF 邊界一樣對于不同類型的波動問題和數(shù)值方法具有普適性.研究結(jié)果顯示，在譜元法中使用caj-MTF 邊界或傳統(tǒng)MTF 邊界時，二階邊界效果最好，而三階及以上階次的邊界則會受到數(shù)值離散誤差的影響反而降低精度，這是目前尚難以解決的問題.MTF 邊界或傳統(tǒng)MTF 邊界時，二階邊界效果最好，而三階及以上階次的邊界則會受到數(shù)值離散誤差的影響反而降低精度，這是目前尚難以解決的問題.

與現(xiàn)有的PML 邊界相比，caj-MTF 邊界的精度相差不大，而在復(fù)雜程度、計算量以及通用性方面則有著巨大優(yōu)勢.與MTF 邊界、黏彈性邊界或一階旁軸近似邊界相比，caj-MTF 邊界在精度方面具有顯著優(yōu)勢，尤其是對于外行波以大角度射向人工邊界和面波透射情形.

本文研究結(jié)果對于使用高精度譜元法進行地震波動模擬時人工邊界條件的優(yōu)化選取具有重要參考價值.雖然文中只給出了二維彈性波的數(shù)值算例，但caj-MTF 邊界或MTF 邊界的施加方式?jīng)Q定了它們可以無差別地應(yīng)用于一維、二維及三維波動離散模型，而其它人工邊界條件則不具備這樣的簡單通用性.下一步擬將該方法拓展應(yīng)用到一些更具實際意義的復(fù)雜波動問題當中，如局部場地的地震波散射問題、飽和介質(zhì)彈性波傳播問題、低頻震源問題等，以探討提高邊界精度的方法以及確保長持時計算穩(wěn)定性所要采取的措施.