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

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

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

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

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

1 一種新型局部透射邊界

1.1 多次透射公式

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

(1)

圖1 用于定義透射邊界的統(tǒng)一局部坐標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即為根據(jù)這些離散參考點的運動外推人工邊節(jié)點運動的表達式。

圖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)的含義及乘積規(guī)則，此時有

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)

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

(3)

為了直觀地體現(xiàn)新邊界與傳統(tǒng)單一人工波速MTF的聯(lián)系和區(qū)別，將邊界式(3)簡記為caj-MTF。

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

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

1.3 數(shù)值計算格式

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

(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的數(shù)值計算格式。

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

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

(9)

(10)

(11)

(T1)1×3=Qv

(12)

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

(13)

(15)

(16)

(17)

(18)

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

式(5)～式(18)組成了前3階caj-MTF實用的數(shù)值計算格式?？梢钥闯?，它與廖振鵬的研究中現(xiàn)有MTF的數(shù)值計算格式一樣簡潔。

2 彈性波反射系數(shù)

這一章將導(dǎo)出所提出的新型局部透射邊界對彈性波的反射系數(shù)，作為后文進行邊界反射特征與參數(shù)優(yōu)化分析的依據(jù)。

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

(19)

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

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

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

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

P波入射：P-P，P-S反射系數(shù)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反射系數(shù)Rsp，Rss

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

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

波勢函數(shù)ΦP，ΦS與數(shù)值模擬使用的波場分量u，w存在如下轉(zhuǎn)換關(guān)系

(22)

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

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

完成上述推導(dǎo)，最后得出P波入射產(chǎn)生P-P反射和P-S反射的反射系數(shù)Rpp和Rps為

(24)

(25)

S波入射產(chǎn)生S-P反射和S-S反射的反射系數(shù)Rsp和Rss為

(26)

(27)

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

(28)

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

3 邊界反射特征及參數(shù)優(yōu)化分析

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

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

圖6 彈性介質(zhì)參數(shù)(cp/cs)和邊界控制參數(shù)(N, caj) 不同取值情形下新型局部透射邊界的彈性波反射系數(shù)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兩種物理波速波動能量的同步吸收

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

Xing等給出了caj-MTF對單分量波動(如聲波或SH波動)的反射系數(shù)(基于等效表達式(19)導(dǎo)出)，通過將橫坐標設(shè)置為視波速，很好地解釋了新邊界對多種物理波速波動能量的同步吸收?；谝暡ㄋ俚膯畏至坎ǚ瓷湎禂?shù)表達式為

(29)

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

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

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

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

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

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

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

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

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

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

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

(30)

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

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

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

(31)

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

(32)

式(30)～式(32)中θp0,θs0和θsc的值可用數(shù)值方法求解。圖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都體現(xiàn)在圖6的反射系數(shù)曲線當(dāng)中，Rps和Rsp的兩個附加零反射角0°和90°亦在其中。當(dāng)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的增大而減小。了解這些角度值及其變化規(guī)律，有助于更好地掌握彈性波情形下局部透射邊界的反射特征。一句話總結(jié)，由彈性介質(zhì)自身特性帶來的附加零反射角是對邊界精度有利的因素。

3.3 S-P大反射問題

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

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

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

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

(33)

圖9繪出了式(33)在不同邊界參數(shù)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的基礎(chǔ)性放大貢獻，參數(shù)cp/cs是逐步增長的放大倍數(shù)，二者組合形成一個幅值約為2～10 (考慮cp/cs從1.2變化到5.0)的反射因子；左上子圖顯示，式(33)前半部分中各個反射因子的幅值隨著cp/cs的變化有可能達到遠超過1的水平，這是由于caj/cs位于分子caj/cp位于分母所造成的。S-P反射的峰值為這些反射因子的乘積，當(dāng)后者都遠超過1時，會導(dǎo)致峰值Rsp(θsc)出現(xiàn)極大的異常值，如幾十至數(shù)百。

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

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

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

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

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

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

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

4 算例驗證

4.1 均勻介質(zhì)波源問題

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

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

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

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

圖14繪出了cp/cs=4介質(zhì)模擬中(x,z)=(300 m, 500 m)處觀察點的u分量時程及其邊界反射誤差。此處提取了6種參數(shù)方案下的結(jié)果，已在圖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結(jié)果顯示，一些參數(shù)方案下u分量的邊界反射誤差會達到接近甚至超過入射波幅值的程度，這遠遠超出了人工邊界條件所能接受的合理誤差范圍。這表明在大縱橫波速比介質(zhì)(cp/cs=4)彈性波模擬中，S-P大反射對邊界精度具有嚴重破壞作用，必須采用本文所提出的方案來避開這個問題。

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

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

4.3 復(fù)雜不均勻介質(zhì)問題

圖17 復(fù)雜不均勻介質(zhì)模型Fig.17 A model of complicated heterogeneous media

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

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

5 結(jié) 論

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

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

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

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

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

Vol．41 No．12 2022