彈性波模擬中局部透射邊界的反射特征及參數優化

邢浩潔， 李小軍， 劉愛文， 李鴻晶， 周正華， 陳 蘇

(1. 中國地震局 地球物理研究所，北京 100081； 2. 北京工業大學 建筑工程學院，北京 100124；3. 南京工業大學 土木工程學院，南京 211816； 4. 南京工業大學 交通運輸工程學院，南京 211816)

地球物理、地震工程、計算聲學等領域的無限介質波動問題常常使用有限計算模型來模擬，此時需要設置適當的人工邊界條件來消除人為截取邊界所造成的虛假反射，以確保模擬結果的合理性[1-2]。提高復雜波動問題模擬中的邊界精度是地球物理正、反演模擬，強地震動場地效應分析，土與結構動力相互作用計算等重要技術的基礎共性需求。不同于力學領域傳統的Dirichlet邊界(固定邊界)和Neumann邊界(自由邊界)會將射向邊界的外行波完全反射回計算區域，人工邊界條件采用時空外推、應力平衡或區域衰減等方式，使計算得到的邊界節點運動與原來實際問題中該處運動基本一致，從而實現外行波在人工邊界上的無反射(或表述為透射、吸收、輻射、單向波動等)。

多次透射公式(multi-transmitting formula, MTF)[3-4]是一種重要的人工邊界條件，它基于直接模擬任意外行波的一般傳播過程而建立，具有形式簡單、精度可控、易于與各種內域離散格式相結合等優點。MTF已廣泛用于海洋地震工程所關心的兩相介質波動問題[5-7]、局部地形對彈性波的散射問題[8]、大壩-地基系統地震反應分析[9-10]、核電結構-場地土系統地震反應分析[11-12]等數值模擬當中。近幾年，MTF被用于高精度譜元法的地震波動模擬[13-16]，能夠與這種不等距節點的離散格式很好地結合，顯示出該邊界具有極強的適應能力。除了上述應用之外，透射邊界的穩定實現是一直以來該邊界研究的重點，研究結果主要包括：邊界失穩是隨著波在計算模型內部和邊界之間來回反射逐步積累的；一種小量修正的MTF公式能夠消除漂移失穩，在整個計算模型中注入小阻尼能夠消除振蕩失穩[17]；GKS(Gustafsson-Kreiss-Sundstr?m)定理的群速度解釋表明在集中質量有限元格式中MTF會出現振蕩失穩，而在二階中心差分格式中則不會[18]；在邊界附近添加彈簧和阻尼力學元件能夠抑制MTF的漂移失穩[19]；在一種修正的集中質量有限元格式中MTF能夠給出長持時穩定的模擬結果[20]；在波導問題模擬中可以通過控制單元長寬比來避免MTF出現振蕩失穩[21]。

為了提高透射邊界模擬復雜波動問題的精度，文獻[22]提出一種在多個不同方向上分別使用MTF再將結果加權平均的多向透射邊界；文獻[23]基于無窮大波速假定提出一種推廣的透射邊界，它在原MTF基礎上增加了對接近掠入射的外行波的模擬能力。最近，Xing等[24-25]基于已有研究揭示的MTF與Higdon邊界之間的密切聯系，將后者所具有的多種計算波速配置引入前者，提出一種基于多人工波速的新型局部透射邊界。然而，目前關于透射邊界精度嚴格的理論分析(這里指反射系數)主要基于單分量波動給出，而地震波作為一種彈性波，是典型的多分量波動。Shi等導出了現有MTF邊界對飽和多孔介質彈性波的反射系數，其中二階MTF對一些波動成分的反射比例能達到0.25～0.35左右，這大大高于單分量波動的反射系數所給出的結果。因此，需要專門研究多分量波動情形下透射邊界的反射特征并進行參數優化。本文將針對新提出的多人工波速局部透射邊界，導出其彈性波反射系數的嚴格表達式，據此詳細探討影響邊界反射特征的各方面因素，最后給出復雜波動情形的邊界參數取值建議。

1 一種新型局部透射邊界

1.1 多次透射公式

在圖1所示的局部復雜場地模型中，圍繞各個人工邊界節點統一定義一種局部坐標s0t，其空間軸s位于從邊界節點0指向內域的任意離散網格線上，正方向向內，時間軸t與整體模型一致。那么，用于推算任意時刻各個人工邊界節點運動的多次透射公式(MTF)[1]為

(1)

圖1 用于定義透射邊界的統一局部坐標s0tFig.1 The unified local coordinates s0t that are used to define transmitting boundaries

如果用B(ca)表示時空移動算子，其含義為

B(ca)u(s,t) ≡B(caΔt, -Δt)u(s,t) =u(s+caΔt,t-Δt)

B2(ca)u(s,t) =B(ca)[B(ca)u(s,t)] =u(s+2caΔt,t-2Δt)

依次類推，則MTF可以簡潔地表示為如下離散算子形式

[I-B(ca)]Nu(s,t)=0

(2)

圖2繪出了二階和三階MTF的離散參考點在局部坐標s0t中所處的位置。可以看出，這些參考點依次由離散算子B(ca)及其乘方B2(ca),B3(ca), …所確定。由于這里只有一種人工波速ca, 故所有參考點等距離分布于一條傾斜直線(圖2中虛線)上。MTF即為根據這些離散參考點的運動外推人工邊節點運動的表達式。

圖2 二階和三階MTF邊界的離散參考點Fig.2 Discrete reference points of the 2nd-and 3rd-order MTF boundaries

1.2 一種新型局部透射邊界

將邊界表達式(2)中的N個相同的離散算子I-B(ca)替換為基于多個人工波速caj(j= 1, …,N)的N個不同的離散算子I-B(caj)。利用上述時空移動算子B(ca)的含義及乘積規則，此時有

B(ca1)u(s,t)≡B(ca1Δt,-Δt)u(s,t)=u(s+ca1Δt,t-Δt),

B(ca1)B(ca2)u(s,t) =u(s+ca1Δt+ca2Δt,t-2Δt),

B(ca1)B(ca2)B(ca3)u(s,t)=u(s+ca1Δt+ca2Δt+ca3Δt,t-3Δt)

依次類推。于是，導出一種基于多個人工波速的新型局部透射邊界，其表達式為

(3)

為了直觀地體現新邊界與傳統單一人工波速MTF的聯系和區別，將邊界式(3)簡記為caj-MTF。

圖3繪出了二階和三階caj-MTF的離散參考點在局部坐標s0t中所處的位置。由此可以看出，新邊界對MTF的優化是通過將原來單個參考點運動2倍、3倍分別替換為2個、3個不同參考點運動之和來實現的。如果從MTF“誤差波多次透射”的物理機制角度來分析caj-MTF，那么后者可以看作是保留了“多次透射”過程，但優化了對各階“誤差波”的描述。

圖3 二階和三階caj-MTF邊界的離散參考點Fig.3 Discrete reference points of the 2nd-and 3rd-order caj-MTF boundaries

1.3 數值計算格式

caj-MTF的數值計算格式可分為兩種情形：對于離散節點在人工邊界附近為等距離分布情形(如有限元或各種有限差分模型)，可以借鑒景立平研究的算子替換法。對于離散節點在人工邊界附近為非等距分布情形(如譜元離散模型)，則可以采用廖振鵬研究的6.1.1節的簡單內插方法。第一種情形最為常見，這里給出其數值計算格式的具體表達式，第二種情形將另文探討。

(4)

其中

sa1=ca1Δt/Δs,sa2=ca2Δt/Δs,sa3=ca3Δt/Δs

(5)

qx0=(sa1-1)(sa1-2)/2,qx1=-sa1(sa1-2),

qx2=sa1(sa1-1)/2

(6)

rx0=(sa2-1)(sa2-2)/2,rx1=-sa2(sa2-2),

rx2=sa2(sa2-1)/2

(7)

sx0=(sa3-1)(sa3-2)/2,sx1=-sa3(sa3-2),

sx2=sa3(sa3-1)/2

(8)

將式(4)中3個乘積項分別只保留一項、兩項，即為一階、二階caj-MTF的數值計算格式。

圖4 caj-MTF數值計算格式中離散節點的統一局部編號Fig.4 Unified local numbering of the discrete nodes involved in the numerical schemes of caj-MTF

逐步展開式(4)，并利用圖4中離散節點的統一局部編號，導出一階、二階、三階caj-MTF數值計算格式的最終形式如下

(9)

(10)

(11)

(T1)1×3=Qv

(12)

(T2)1×8=[T2_1-T2_2]

(13)

(15)

(16)

(17)

(18)

式中：帶下標的粗體變量表示由計算系數組成的行向量或由節點波場值組成的列向量；上標T表示轉置；括號下標( )1×3, ()1×8, ()15×1等僅是為了便于理解，用來標示出各個向量的維度，在實際編程時應忽略。式(15)中涉及的統一局部節點編號見圖4。

式(5)～式(18)組成了前3階caj-MTF實用的數值計算格式。可以看出，它與廖振鵬的研究中現有MTF的數值計算格式一樣簡潔。

2 彈性波反射系數

這一章將導出所提出的新型局部透射邊界對彈性波的反射系數，作為后文進行邊界反射特征與參數優化分析的依據。

現有MTF的反射系數基于等距離散參考點上的振動解導出，但本文caj-MTF的離散參考點為不等距分布，難以沿用這個方法。不過，根據現有研究所揭示的離散形式MTF邊界與經典的連續形式Clayton-Engquist以及Higdon邊界之間的密切聯系，可以導出與caj-MTF等價的連續邊界形式，在后者基礎上導出所需的彈性波反射系數。刑浩潔等的研究證明了邊界式(3)中各個離散算子(I-B(caj))等價于連續形式的一維波動微分算子(?/?t-caj?/?s)，因此與邊界(3)等價的連續邊界形式為

(19)

式(19)仍然定義在圖1所示的統一局部坐標當中，空間軸s可以表示整體模型的x，y，z軸或者任意一條傾斜的離散網格線。

從文獻[26-28]得知，求解邊界反射系數的步驟為：先寫出入射波的諧波表達式；再將入射波與反射波(等于入射波乘以反射系數)組成的波場代入邊界條件；約去不必要的部分，即可得到反射系數。平直界面對彈性波的反射有P波入射和S波入射兩組模式，前者存在P-P，P-S反射系數Rpp，Rps，后者有S-P，S-S反射系數Rsp，Rss。不失一般性，考慮圖5所示的坐標系與外行波情形。

圖5 平面波的透射角度與視波速Fig.5 Incident angles and apparent propagation velocities of plane waves

ΦP和ΦS為P波和S波勢函數，則兩組反射模式下邊界附近的波場分別為

P波入射：P-P，P-S反射系數Rpp，Rps

ΦP=exp[i(ωt+kpxx+kpzz)]+Rppexp[i(ωt-kpxx+kpzz)]，

ΦS=Rpsexp[i(ωt-ksxx+kszz)]

(20)

S波入射：S-P，S-S反射系數Rsp，Rss

ΦP=Rspexp[i(ωt-kpxx+kpzz)],

ΦS=exp[i(ωt+ksxx+kszz)]+Rssexp[i(ωt-ksxx+kszz)](21)

波勢函數ΦP，ΦS與數值模擬使用的波場分量u，w存在如下轉換關系

(22)

此時caj-MTF的等價連續邊界式 (19)的具體形式為

將式(20)代入式(22)，再將結果代入(23)，求解聯立方程組，即可導出P波反射系數Rpp和Rps；同理，式(21)、式(22)和式(23)可以導出S波反射系數Rsp和Rss。

完成上述推導，最后得出P波入射產生P-P反射和P-S反射的反射系數Rpp和Rps為

(24)

(25)

S波入射產生S-P反射和S-S反射的反射系數Rsp和Rss為

(26)

(27)

式(24)～式(27)中同時含有θp和θs兩個角度。在分析中，需要將Rpp，Rps表示成只含有θp的形式，將Rsp，Rss表示成只含有θs的形式。根據Snell定律：入射波、透射波和反射波沿界面的視傳播速度相等，即cpz=csz(見圖5)。于是得到關系式cp/sinθp=cs/sinθs，即

(28)

將式(28)第二式代入式(24)和式(25)，將式(28)第一式代入式(26)和式(27)，即為分析所用的兩組(4種)彈性波反射系數。

3 邊界反射特征及參數優化分析

圖6繪出了16種情形下式(24)～式(28)各個彈性波反射系數的幅值曲線。這里考慮了兩種彈性介質性質cp/cs=2和4，后者為大縱橫波速比情形；兩種邊界階次N=2和3，二階邊界最為常用，三階邊界也有一定應用；3種單一人工波速方案和一種多人工波速方案，前3種方案與現有MTF邊界相同，最后一種方案是本文新型邊界的特色。

人工邊界的目標是實現對外行波的無反射，即所有反射系數在全部入射角范圍內為零，但觀察圖6可知這種理想狀況難以達到。不過，已有研究和數值試驗表明，接近掠入射附近角度(如70°～90°)的波動能量通常少到可以忽略不計，因此反射系數只需在小角度和中等角度范圍內(如0°～70°)足夠小，就能給出接近理想的模擬結果。以此為考察標準，接下來結合圖6，詳細探討彈性波情形下影響局部透射邊界反射特征的三方面因素并進行參數優化分析。

圖6 彈性介質參數(cp/cs)和邊界控制參數(N, caj) 不同取值情形下新型局部透射邊界的彈性波反射系數Fig.6 Elastic-wave reflection coefficients of the new local transmitting boundary with different values of the elastic medium parameter cp/cs and the boundary control parameters N and caj

3.1 P，S兩種物理波速波動能量的同步吸收

新型局部透射邊界所具有的多人工波速參數正是為了更好地解決存在多種物理波速的復雜波動問題，尤其是當不同物理波速之間差異較大時，新邊界比傳統單一人工波速邊界具有顯著優勢。

Xing等給出了caj-MTF對單分量波動(如聲波或SH波動)的反射系數(基于等效表達式(19)導出)，通過將橫坐標設置為視波速，很好地解釋了新邊界對多種物理波速波動能量的同步吸收。基于視波速的單分量波反射系數表達式為

(29)

式中，cn為外行波沿邊界計算方向的視波速，根據圖5，此時cn≡cx=ω/kx=v/cosθ，v為某種波動成分的物理波速，它可以是聲波波速c、P波波速cp、S波波速cs等，θ為每種波的透射角度。式(29)針對局部透射邊界模擬任意復雜波動問題的精度提供了一種無差別的度量標準，因此它比傳統基于透射角度的單分量波反射系數更具一般性。

新型局部透射邊界對P，S兩種物理波速波動能量的同步吸收，可以通過將其中兩個人工波速參數設定為ca1=cs和ca2=cp來實現(這對應于二階邊界)，對于更高階邊界的其余人工波速參數，可以在cs到略大于cp范圍內取值。這種做法的原理是：根據式(29)，此時在cn=cs和cn=cp處為兩個零反射點，零反射點附近的反射系數必然處于低值水平，因此在傳播主要波動能量的垂直入射及其附近一定入射角范圍內，P，S波都能夠被很好地吸收。一階邊界只有一個人工波速參數，無法較好地兼顧對P波和S波的吸收；三階邊界能夠進一步降低二階邊界兩個零反射點附近的反射幅值，尤其是當P，S波速差異較大時，有必要使用它；更高階邊界繼續提升精度的效果不明顯，卻更容易出現失穩問題，使用和討論得很少。

類似地，利用本文給出的彈性波反射系數可以得出同樣的結論。為證明這一點，可以將式(24)～式(27)改寫成基于視波速cn的形式，并將其與式(29)的曲線繪制于同一幅圖中，觀察二者相似度。根據上述視波速定義，對于Rpp和Rps，cn=cp/cosθp，有θp=arcos(cp/cn)；對于Rsp和Rss，cn=cs/cosθs，有θs=arcos(cs/cn)。于是，式(24)～式(28)均可改寫為關于變量caj,cn,cs,cp的形式。

為節省篇幅，這里不再寫出基于視波速的彈性波反射系數表達式。圖7給出了彈性波反射系數與單分量波反射系數繪制于同一坐標系的結果，三幅子圖分別為一階、二階、三階caj-MTF邊界。圖7中，v為某個可以被約去的速度值，用于將各個波速無量綱化。考慮cs=v，cp=3v，縱橫波速比cp/cs=3。橫坐標中視波速cn的值從v到6v，這涵蓋了所關心的S波和P波從垂直入射到以小到中等角度入射的視波速范圍。

圖7 基于視波速cn的彈性波與單分量波反射系數Fig.7 Elastic-wave and single-component-wave reflection coefficients in the context of the apparent velocity cn

在圖7中，S波反射系數Rsp和Rss為灰色粗實線和虛線，位于坐標范圍[v, 3v]，對應于θs從0°～70.5°(即arcos(v/3v))；P波反射系數Rpp和Rps為黑色粗實線和虛線，位于坐標范圍[3v, 6v]，對應于θp為0°～60°(即arcos(3v/6v))；式(29)給出的單分量波反射系數R為細點線。總體上，彈性波曲線Rsp，Rss，Rpp，Rps所表示的整體反射特征與單分量波曲線R的變化規律基本一致，如：Rss與該段R完全重合；兩類反射曲線具有共同的零反射點cn=v(即ca1)和cn=3v(即ca2)；Rpp與該段R在零反射點附近比較接近。這證明了就局部透射邊界對于復雜波動問題中多種物理波速的處理而言，式(29)的確是一種簡單可靠的度量標準。

圖7中Rsp，Rss，Rpp，Rps的整體特征還存在幾處與R不一致的地方：Rsp在一定范圍內出現異常大值；Rpp不像該段R那樣一直增大，而是先增大后下降；Rps的值要遠小于該段R的值。這些差異是由P波和S波耦合效應引起的，它們可歸結為附加零反射角和S-P大反射兩方面因素。

3.2 附加零反射角和S波臨界角

零反射角指邊界反射系數為零的透射角度，前文將其稱為零反射點(亦可稱為零反射中心)。從圖6和圖7不難看出，在所關心的透射角度或視波速范圍內，零反射點越多對邊界精度越有利。

彈性波反射系數Rpp，Rps，Rsp，Rss有兩種零反射角。一種是由邊界參數caj所確定的，根據式(24)～式(27)乘號之前的部分計算，它們用視波速表述為cn=caj，若用透射角表述，則P波為cp/cosθp=caj，S波為cs/cosθs=caj。另一種是由彈性介質自身特性即縱橫波速比cp/cs所確定的，根據式(24)～式(27)乘號之后的部分計算。本文將第二種角度稱為彈性介質自身特性帶來的附加零反射角，其值分別如下。

Rpp有一個附加零反射角θp0，其值根據式(30)確定

(30)

Rps有兩個附加零反射角，為θp0=0°和90°；

Rsp有兩個附加零反射角，為θs0=0°和90°；

Rss有一個附加零反射角θs0，其值根據式(31)確定

(31)

除了附加零反射角θp0和θs0之外，S-P反射中還有一種S波臨界角θsc，對應于反射P波角度為臨界狀態90°的情形。根據式(28)第二式可得θsc為

(32)

式(30)～式(32)中θp0,θs0和θsc的值可用數值方法求解。圖8繪出了它們的值隨cp/cs的變化曲線。

圖8 Rpp,Rss的附加零反射角θp0,θs0以及S波臨界角θscFig.8 The additional zero-reflection angles θp0, θs0 for Rpp, Rss and the S-wave critical angle θsc

圖8中cp/cs=2和4處的附加零反射角θp0,θs0和S波臨界角θsc都體現在圖6的反射系數曲線當中，Rps和Rsp的兩個附加零反射角0°和90°亦在其中。當cp/cs從1.42逐步變化到5時，Rpp的θp0從54.7°漸漸增大至78.7°，Rss的θs0從35.3°慢慢減小到11.3°，S波臨界角從45°逐步減小至11.5°。θsc的值略大于θs0，二者差異隨著cp/cs的增大而減小。了解這些角度值及其變化規律，有助于更好地掌握彈性波情形下局部透射邊界的反射特征。一句話總結，由彈性介質自身特性帶來的附加零反射角是對邊界精度有利的因素。

3.3 S-P大反射問題

圖6和圖7曲線顯示，S-P反射系數Rsp在S波臨界角θsc附近一定范圍內，常常會出現異常大值，其幅值可能達到接近甚至遠超過1的程度。一旦出現S-P大反射，在主要波動能量所處的小角度及中等入射角范圍內，常用的二階或三階邊界的精度遠不能令人滿意，因此必須深入分析其原因并加以解決。

首先觀察圖6各幅子圖，S-P大反射表現出以下特征：從左至右進行對比，人工波速caj全取為cs時未出現S-P大反射，其他參數情形均不同程度地出現了大反射，其中人工波速caj全取為cp時，異常大反射最為嚴重；S-P大反射在臨界角θsc處有一個幅值尖峰，此角度附近的幅值迅速降低；從上往下對比，彈性介質縱橫波速比cp/cs=4時的S-P大反射比cp/cs=2時嚴重得多，而三階邊界與二階邊界的區別在于新增一個反射因子(注：反射因子指式(24)～式(27)中相乘的每一項)所帶來的乘積效應。

然后分析反射系數Rsp表達式(26)，并與Rpp，Rps，Rss的式(24)、式(25)、式(27)進行比較。最容易看出的是式(26)中存在一個放大因子cp/cs，縱橫波速比越大，其對Rsp幅值的放大作用越顯著；與之相反，式(25)中這個因子為cs/cp，對Rps幅值有縮小作用；式(24)和式(27)中沒有類似因子。另一個比較直觀的觀察是caj/cs的值要大于caj/cp，而在Rsp中出現了前者位于分子后者位于分母的組合，導致其反射因子的幅值較大，可能達到接近甚至遠超過1的水平；與之相比，Rpp，Rps，Rss的各個反射因子中均沒有類似組合。這兩個直觀觀察粗略地解釋了Rsp有時會出現異常大反射，而Rpp，Rps，Rss均沒有類似問題的原因。

為嚴謹起見，進一步對Rsp表達式(26)進行定量分析。S-P大反射的幅值在S波臨界角θsc(見式(32)和圖8)處有一個尖峰，圖6未充分顯示其幅值變化規律，這里對其進行探究。將θsc代入式(26)，得到此處Rsp的值為

(33)

圖9繪出了式(33)在不同邊界參數caj(j=1,…,N)取值下，各個反射因子的幅值及其組合得到的S-P反射的峰值Rsp(θsc)隨著縱橫波速比cp/cs的變化。

圖9 S-P大反射峰值分析Fig.9 Analysis of the peak value of S-P large reflection

圖9從定量角度嚴格證明了上述基于直觀觀察對S-P大反射原因做出的解釋。具體為：圖9右上角子圖顯示，sin(2θsc)/cos(π/2-θsc)提供了幅值接近于2的基礎性放大貢獻，參數cp/cs是逐步增長的放大倍數，二者組合形成一個幅值約為2～10 (考慮cp/cs從1.2變化到5.0)的反射因子；左上子圖顯示，式(33)前半部分中各個反射因子的幅值隨著cp/cs的變化有可能達到遠超過1的水平，這是由于caj/cs位于分子caj/cp位于分母所造成的。S-P反射的峰值為這些反射因子的乘積，當后者都遠超過1時，會導致峰值Rsp(θsc)出現極大的異常值，如幾十至數百。

解決S-P大反射問題的方法已經由基于圖6的直觀觀察和圖9的定量分析所暗示：只要將所有人工波速參數嚴格設定為S波波速，即caj=cs(j=1, …,N)，就不會出現S-P大反射。其原理為：圖9左上子圖顯示當caj=cs時，式(33)前半部分中反射因子的幅值接近于零，若干個這樣因子的乘積能夠有效地抵消后半部分天然具有的約為2～10的大幅值，從而使得S-P反射峰值Rsp(θsc)降至人工邊界精度所要求的低值水平(如小于0.1)。

但是，上述解決S-P大反射問題的邊界參數取值方案與3.1節所給出的解決P和S兩種物理波速問題的方案存在沖突，后者要求至少有兩個caj參數取值為caj=(cs,cp)。從圖9左上子圖可知，caj=cp的反射因子幅值能夠達到1～3左右(考慮cp/cs位于2～4)。此時，為了抵消這個因子和后半部分因子的大幅值，需要至少兩個caj=cs的反射因子，這對應于控制參數為caj=(cs,cs,cp)的三階局部透射邊界。這種參數配置在圖6中沒有出現，表明Xing的研究基于單分量波反射系數式(29)所建議的三階邊界參數caj= (cs,cs,cp)在彈性波情形下并不是好的選擇。圖10給出了本文建議的caj=(cs,cs,cp)三階邊界的彈性波反射系數。

圖10 參數為(cs, cs, cp)的三階邊界的彈性波反射系數Fig.10 Elastic-wave reflection coefficients of the 3rd-order boundary with parameters of (cs, cs, cp)

與圖6各幅子圖比較可知，圖10邊界方案效果更好，它既實現了對P和S波的同步吸收，又不存在S-P大反射問題。進一步研究發現，避免S-P大反射的方法可以退化為：只需將大部分(如1/1, 2/2, 2/3, 3/4, 3/5, …)人工波速參數嚴格設定為S波波速，就可以避開這個問題。

最后需要指出，現有局部透射邊界最常使用的是二階形式，但各種參數方案下的二階邊界都無法同時兼顧對P,S波的同步吸收與避開S-P大反射兩個方面。不過，根據上述理論分析并觀察圖6和圖9曲線，可以給出二階邊界的使用建議如下：可用于縱橫波速比不大(如cp/cs不超過2.5)的彈性介質；優先考慮采用caj=(cs,cp) 作為邊界參數；嚴格避免使用參數方案caj=(cp,cp)；理論上無法實現對各種波動成分的完美吸收，對于縱橫波速比不大的問題，能給出較高精度模擬結果。除此之外，對于縱橫波速較大(如cp/cs>2.5)的問題，則建議使用圖10給出的(cs,cs,cp)三階邊界。由于巖土介質的力學性質常常以泊松比υ來表述，為便于了解大縱橫波速比cp/cs所對應的巖土介質類型，圖11繪出了υ與cp/cs的關系曲線。

圖11 彈性介質泊松比與縱橫波速比的關系Fig.11 The relation of Poisson’s ratio and the ratio of longitudinal and transverse velocities in elastic media

圖11表明，局部透射邊界S-P大反射問題比較突出的大縱橫波速比彈性介質，比如cp/cs>2.5，其泊松比范圍為υ∈(0.405, 0.500)。河谷、湖泊、沉積盆地、海洋沉積物等軟土介質經常能夠達到這么大的泊松比。因此，當基于彈性波方程來模擬這類介質中的地震波傳播時，需要采用本文提出的理論與方法來提高局部透射邊界的精度，確保模擬結果合理性。在當前大型沉積盆地上的城市抗震安全、海洋地震工程等重要問題受到廣泛關注的背景下，本文工作有著巨大應用和推廣價值。

4 算例驗證

4.1 均勻介質波源問題

圖12、圖13分別為cp/cs=2,cp/cs=4兩種彈性介質在不同邊界參數下的模擬結果。圖12、圖13給出了水平分量u在兩個典型時刻的波場快照，其中1.5 s，0.85 s為P波穿過邊界，2.6 s，2.4 s為S波穿過邊界。

圖12 不同邊界參數對cp/cs=2彈性介質的模擬結果Fig.12 Modeling results of elastic waves (cp/cs=2) using different boundary parameters

圖13 不同邊界參數對cp/cs=4彈性介質的模擬結果Fig.13 Modeling results of elastic waves (cp/cs=4) using different boundary parameters

圖12、圖13給出的邊界反射規律與圖6、圖10的反射系數曲線及第3章理論分析結果完全相符。主要包括：所有單一人工波速方案都無法很好地實現對P波和S波的同步吸收，而含有caj=(cs,cp)的多人工波速方案能夠做到。小縱橫波速比(cp/cs=2)彈性介質僅在邊界參數為(cp,cp)和(cp,cp,cp)時出現了輕度的S-P反射問題，但大縱橫波速比(cp/cs=4)彈性介質出現了極其嚴重的S-P大反射。(cp,cp)和(cp,cp,cp)方案最容易出現S-P大反射且幅值最大；(cs,cs)、(cs,cs,cs)，(cs,cs,cp)這幾種caj=cs數目占優勢的方案都沒有出現S-P大反射。值得注意的是，在圖12中表現良好的(cs,cp,cp)方案(Xing等推薦的三階邊界方案)，在圖13中也出現了顯著的S-P反射誤差。在圖12和圖13中，本文所建議的(cs,cs,cp)方案三階邊界均給出了非常滿意的結果。

圖14繪出了cp/cs=4介質模擬中(x,z)=(300 m, 500 m)處觀察點的u分量時程及其邊界反射誤差。此處提取了6種參數方案下的結果，已在圖13中用倒三角形示出。在計算邊界反射誤差時，以(cs,cs,cp)方案作為參考解。

圖14 觀察點的u分量時程及邊界反射誤差(cp/cs=4)Fig.14 Time histories of the component u at the observation point and their boundary-reflection errors (cp/cs=4)

圖14結果顯示，一些參數方案下u分量的邊界反射誤差會達到接近甚至超過入射波幅值的程度，這遠遠超出了人工邊界條件所能接受的合理誤差范圍。這表明在大縱橫波速比介質(cp/cs=4)彈性波模擬中，S-P大反射對邊界精度具有嚴重破壞作用，必須采用本文所提出的方案來避開這個問題。

另外需要指出，圖14給出的u分量反射誤差并未達到如圖9所示的數十倍以上的程度。這種差異具有合理性，因為：反射系數Rsp是關于波勢函數ΦP和ΦS的，u分量與它們存在轉換關系式(22)；反射系數基于平面諧波給出，而本試驗為暫態的短持時脈沖波；最重要的原因是圖9給出的是S-P反射在臨界角θsc處的峰值Rsp(θsc)，臨界角附近的S-P大反射幅值會迅速降低(參考圖6)，本試驗中的S波入射角與臨界角不同。總之，第3章的理論內容很好地預測了此處數值試驗結果。

4.2 自由地表散射問題

圖15 散射問題波場快照(u分量)Fig.15 Snapshots of a scattering problem (u component)

圖16 (x, z)=(300 m, 240 m)處觀察點的u分量時程Fig.16 Time histories of the component u at the observation point (x, z)=(300 m, 240 m)

圖15中右側邊界對于兩道P波和一道S波實現了完美的同步吸收，表明對于這種視波速可以事先確定的情形，局部透射邊界效果極佳。在底部邊界上，(cs,cs)邊界很好地吸收了S波而對P波有明顯反射，(cp,cp)邊界很好吸收了P波而對S波有很大反射，基于視波速的邊界精確解則實現了對兩種波的同步吸收。觀察(cp,cp)邊界的反射，S波反射分為S-P和S-S兩道，S-P反射的幅值要遠遠大于后者。從圖16可以進一步看出，(cp,cp)方案S-P反射造成的誤差很大，是影響模擬結果的最重要因素。在短持時脈沖波模擬結果中，不合理的邊界方案尚且出現如此顯著的誤差，對于長持時的地震波模擬，更應注意采用合理的邊界方案來確保邊界精度。

4.3 復雜不均勻介質問題

圖17 復雜不均勻介質模型Fig.17 A model of complicated heterogeneous media

圖18為3組模擬得到的波場快照。結果表明：自由邊界會將外行波完全反射回計算區域，其模擬結果與實際無限介質中的波傳播過程相去甚遠。(cp,cp,cp)邊界在左側和右側邊界上都出現顯著的大反射(根據前文理論可知此為S-P大反射)，嚴重影響了模擬結果的可用性。據此以推斷對于長持時的實際地震波問題，邊界反射誤差不斷疊加累積，會導致更加不合理的結果。(cs,cs,cp)邊界很好地消除了人工邊界上的虛假反射波，說明它所推算的各個人工邊界節點的運動比較符合實際無限介質問題中該處的運動，所以有限模型的模擬結果再現了實際波動過程。

圖18 不均勻介質彈性波模擬結果Fig.18 Modeling results of elastic waves in heterogeneous media

5 結 論

通過將多個人工波速參數引入廖氏透射邊界的算子表達式，提出一種優化的新型局部透射邊界，并給出簡潔實用的數值計算格式。新邊界比廖氏透射邊界更具一般性，后者為前者多個人工波速取相同值的特例。局部透射邊界對彈性波的反射特征受到P波與S波耦合效應的重要影響，但現有透射邊界理論所給出的基于單分量波(如聲波或SH波)的反射系數或“誤差波多次透射”的物理機制并未考慮這個效應。本文給出了彈性波情形下局部透射邊界的兩組(4種)反射系數Rpp，Rps和Rsp，Rss的明確定義，推導出反射系數的嚴格表達式，據此詳細探討了邊界反射特征并進行參數優化分析。理論與數值試驗相結合，得到以下研究結論：

(1) 新型局部透射邊界具有的多個人工波速配置能夠比現有基于單一人工波速的透射邊界更好地處理含有多種物理波速的復雜波動問題。當一組人工波速取值中同時含有P波和S波波速時，新邊界能夠很好地實現對這兩種波的同步吸收。

(2) 彈性波反射系數中除了由邊界參數決定的零反射角之外，還存在由P波與S波耦合效應(彈性介質自身特性)帶來的附加零反射角。零反射角越多，在0°～90°內分布越均勻，邊界精度越好。

(3) 在縱橫波速比較大(如cp/cs>2.5)的彈性介質中，S-P反射系數Rsp有時會出現幅值接近甚至遠超過1的異常放大，這會嚴重破壞邊界精度。不過，當大部分(如1/1, 2/2, 2/3, 3/4, 3/5, …)人工波速參數被嚴格設定為S波波速時，就不會出現S-P大反射問題。

(4) 局部透射邊界S-P大反射問題比較突出的大縱橫波速比介質(如cp/cs>2.5)，其泊松比范圍為v∈(0.405, 0.500)，河谷場地、沉積盆地、海洋沉積物等軟土介質通常具有較大泊松比。因此，當基于彈性波方程來模擬這類介質中的地震波傳播時，需要采用本文方法來提高邊界精度，確保模擬結果的合理性。在目前大型沉積盆地上的城市抗震安全、海洋地震工程等受到廣泛關注的背景下，本文工作有著巨大應用和推廣價值。

Vol．41 No．12 2022