彈性波模擬中局部透射邊界的反射特征及參數(shù)優(yōu)化

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

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

1 一種新型局部透射邊界

1.1 多次透射公式

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

(1)

圖1 用于定義透射邊界的統(tǒng)一局部坐標(biāo)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的離散參考點在局部坐標(biāo)s0t中所處的位置。可以看出，這些參考點依次由離散算子B(ca)及其乘方B2(ca),B3(ca), …所確定。由于這里只有一種人工波速ca, 故所有參考點等距離分布于一條傾斜直線(圖2中虛線)上。MTF即為根據(jù)這些離散參考點的運動外推人工邊節(jié)點運動的表達(dá)式。

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

1.2 一種新型局部透射邊界

將邊界表達(dá)式(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)出一種基于多個人工波速的新型局部透射邊界，其表達(dá)式為

(3)

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

圖3繪出了二階和三階caj-MTF的離散參考點在局部坐標(biāo)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ù)值計算格式的具體表達(dá)式，第二種情形將另文探討。

(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)

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

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

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

這一章將導(dǎo)出所提出的新型局部透射邊界對彈性波的反射系數(shù)，作為后文進(jìn)行邊界反射特征與參數(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)一局部坐標(biāo)當(dāng)中，空間軸s可以表示整體模型的x，y，z軸或者任意一條傾斜的離散網(wǎng)格線。

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

圖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邊界相同，最后一種方案是本文新型邊界的特色。

人工邊界的目標(biāo)是實現(xiàn)對外行波的無反射，即所有反射系數(shù)在全部入射角范圍內(nèi)為零，但觀察圖6可知這種理想狀況難以達(dá)到。不過，已有研究和數(shù)值試驗表明，接近掠入射附近角度(如70°～90°)的波動能量通常少到可以忽略不計，因此反射系數(shù)只需在小角度和中等角度范圍內(nèi)(如0°～70°)足夠小，就能給出接近理想的模擬結(jié)果。以此為考察標(biāo)準(zhǔn)，接下來結(jié)合圖6，詳細(xì)探討彈性波情形下影響局部透射邊界反射特征的三方面因素并進(jìn)行參數(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ù)(基于等效表達(dá)式(19)導(dǎo)出)，通過將橫坐標(biāo)設(shè)置為視波速，很好地解釋了新邊界對多種物理波速波動能量的同步吸收?；谝暡ㄋ俚膯畏至坎ǚ瓷湎禂?shù)表達(dá)式為

(29)

式中，cn為外行波沿邊界計算方向的視波速，根據(jù)圖5，此時cn≡cx=ω/kx=v/cosθ，v為某種波動成分的物理波速，它可以是聲波波速c、P波波速cp、S波波速cs等，θ為每種波的透射角度。式(29)針對局部透射邊界模擬任意復(fù)雜波動問題的精度提供了一種無差別的度量標(biāo)準(zhǔn)，因此它比傳統(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波的吸收；三階邊界能夠進(jìn)一步降低二階邊界兩個零反射點附近的反射幅值，尤其是當(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ù)表達(dá)式。圖7給出了彈性波反射系數(shù)與單分量波反射系數(shù)繪制于同一坐標(biāo)系的結(jié)果，三幅子圖分別為一階、二階、三階caj-MTF邊界。圖7中，v為某個可以被約去的速度值，用于將各個波速無量綱化?？紤]cs=v，cp=3v，縱橫波速比cp/cs=3。橫坐標(biāo)中視波速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為灰色粗實線和虛線，位于坐標(biāo)范圍[v, 3v]，對應(yīng)于θs從0°～70.5°(即arcos(v/3v))；P波反射系數(shù)Rpp和Rps為黑色粗實線和虛線，位于坐標(biāo)范圍[3v, 6v]，對應(yīng)于θp為0°～60°(即arcos(3v/6v))；式(29)給出的單分量波反射系數(shù)R為細(xì)點線。總體上，彈性波曲線Rsp，Rss，Rpp，Rps所表示的整體反射特征與單分量波曲線R的變化規(guī)律基本一致，如：Rss與該段R完全重合；兩類反射曲線具有共同的零反射點cn=v(即ca1)和cn=3v(即ca2)；Rpp與該段R在零反射點附近比較接近。這證明了就局部透射邊界對于復(fù)雜波動問題中多種物理波速的處理而言，式(29)的確是一種簡單可靠的度量標(biāo)準(zhǔn)。

圖7中Rsp，Rss，Rpp，Rps的整體特征還存在幾處與R不一致的地方：Rsp在一定范圍內(nèi)出現(xiàn)異常大值；Rpp不像該段R那樣一直增大，而是先增大后下降；Rps的值要遠(yuǎn)小于該段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)異常大值，其幅值可能達(dá)到接近甚至遠(yuǎn)超過1的程度。一旦出現(xiàn)S-P大反射，在主要波動能量所處的小角度及中等入射角范圍內(nèi)，常用的二階或三階邊界的精度遠(yuǎn)不能令人滿意，因此必須深入分析其原因并加以解決。

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

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

為嚴(yán)謹(jǐn)起見，進(jìn)一步對Rsp表達(dá)式(26)進(jìn)行定量分析。S-P大反射的幅值在S波臨界角θsc(見式(32)和圖8)處有一個尖峰，圖6未充分顯示其幅值變化規(guī)律，這里對其進(jìn)行探究。將θ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從定量角度嚴(yán)格證明了上述基于直觀觀察對S-P大反射原因做出的解釋。具體為：圖9右上角子圖顯示，sin(2θsc)/cos(π/2-θsc)提供了幅值接近于2的基礎(chǔ)性放大貢獻(xiàn)，參數(shù)cp/cs是逐步增長的放大倍數(shù)，二者組合形成一個幅值約為2～10 (考慮cp/cs從1.2變化到5.0)的反射因子；左上子圖顯示，式(33)前半部分中各個反射因子的幅值隨著cp/cs的變化有可能達(dá)到遠(yuǎn)超過1的水平，這是由于caj/cs位于分子caj/cp位于分母所造成的。S-P反射的峰值為這些反射因子的乘積，當(dāng)后者都遠(yuǎn)超過1時，會導(dǎo)致峰值Rsp(θsc)出現(xiàn)極大的異常值，如幾十至數(shù)百。

解決S-P大反射問題的方法已經(jīng)由基于圖6的直觀觀察和圖9的定量分析所暗示：只要將所有人工波速參數(shù)嚴(yán)格設(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的反射因子幅值能夠達(dá)到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大反射問題。進(jìn)一步研究發(fā)現(xiàn)，避免S-P大反射的方法可以退化為：只需將大部分(如1/1, 2/2, 2/3, 3/4, 3/5, …)人工波速參數(shù)嚴(yán)格設(shè)定為S波波速，就可以避開這個問題。

最后需要指出，現(xiàn)有局部透射邊界最常使用的是二階形式，但各種參數(shù)方案下的二階邊界都無法同時兼顧對P,S波的同步吸收與避開S-P大反射兩個方面。不過，根據(jù)上述理論分析并觀察圖6和圖9曲線，可以給出二階邊界的使用建議如下：可用于縱橫波速比不大(如cp/cs不超過2.5)的彈性介質(zhì)；優(yōu)先考慮采用caj=(cs,cp) 作為邊界參數(shù)；嚴(yán)格避免使用參數(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á)到這么大的泊松比。因此，當(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)了極其嚴(yá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分量的邊界反射誤差會達(dá)到接近甚至超過入射波幅值的程度，這遠(yuǎn)遠(yuǎn)超出了人工邊界條件所能接受的合理誤差范圍。這表明在大縱橫波速比介質(zhì)(cp/cs=4)彈性波模擬中，S-P大反射對邊界精度具有嚴(yán)重破壞作用，必須采用本文所提出的方案來避開這個問題。

另外需要指出，圖14給出的u分量反射誤差并未達(dá)到如圖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波入射角與臨界角不同?？傊?章的理論內(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反射的幅值要遠(yuǎn)遠(yuǎn)大于后者。從圖16可以進(jìn)一步看出，(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ì)中的波傳播過程相去甚遠(yuǎn)。(cp,cp,cp)邊界在左側(cè)和右側(cè)邊界上都出現(xiàn)顯著的大反射(根據(jù)前文理論可知此為S-P大反射)，嚴(yán)重影響了模擬結(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ù)引入廖氏透射邊界的算子表達(dá)式，提出一種優(yōu)化的新型局部透射邊界，并給出簡潔實用的數(shù)值計算格式。新邊界比廖氏透射邊界更具一般性，后者為前者多個人工波速取相同值的特例。局部透射邊界對彈性波的反射特征受到P波與S波耦合效應(yīng)的重要影響，但現(xiàn)有透射邊界理論所給出的基于單分量波(如聲波或SH波)的反射系數(shù)或“誤差波多次透射”的物理機制并未考慮這個效應(yīng)。本文給出了彈性波情形下局部透射邊界的兩組(4種)反射系數(shù)Rpp，Rps和Rsp，Rss的明確定義，推導(dǎo)出反射系數(shù)的嚴(yán)格表達(dá)式，據(jù)此詳細(xì)探討了邊界反射特征并進(jìn)行參數(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)幅值接近甚至遠(yuǎn)超過1的異常放大，這會嚴(yán)重破壞邊界精度。不過，當(dāng)大部分(如1/1, 2/2, 2/3, 3/4, 3/5, …)人工波速參數(shù)被嚴(yán)格設(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