基于波場分離的不規(guī)則地形下地震波輸入方法

張佳文，李明超，韓 帥，閆文鈺

(1.水利工程仿真與安全國家重點實驗室，天津大學，天津 300350；2.香港理工大學建筑與房地產學系智能建造實驗室，香港 999077)

不規(guī)則地形條件對地震波傳播路徑和地震烈度異常區(qū)的發(fā)生均有重大影響，主要體現在局部地形對地震波的散射和波型轉換[1-3]，國內外震害資料也反映了其對地表位移和加速度幅值放大效應等的巨大影響[4-6]。因此，場地反應分析一直以來都是地震學及工程抗震領域關注的重點。此外，近年來大量臺陣記錄均證實了地震波入射角度的存在[7]，斜入射地震波對結構物振動響應的危害性成為共識[8-9]，地震波的入射角度與行波效應等問題的引入使得不規(guī)則地形場地地震波場的求解難度進一步升級。

目前，場地反應的求解主要有解析法和數值模擬法[10]，解析方法包括波函數展開法、分離變量法、幾何射線理論等十幾種類型[11-14]，已成功應用在了規(guī)則局部地形中，例如：均質長方形場地[15]、半圓和半橢圓形峽谷[16]、任意圓弧凹陷場地[17]等。然而解析方法限制條件較多，針對場地介質本構關系復雜或形狀不規(guī)則地形條件，全域波場的解析求解難以實現。相比而言，數值模擬方法適用范圍更加廣泛，且能更好地應用到抗震分析中[18]，各類基于有限元法、有限差分法和邊界元法等的研究也正在蓬勃發(fā)展[19-21]。斜入射地震波作用下不規(guī)則地形地震波場的求解作為一個較難的課題，可通過解析推導與數值模擬相結合來實現。

ASHFORD 和SITAR[22]基于廣義一致透射邊界等分析了陡坎地形在傾斜SH 波作用下位移等的放大效應，將不規(guī)則場地分區(qū)域進行求解，并首次引入了此類問題的時域化解法，但波場計算的精度不高且并未研究平面內的波場；李山有和廖振鵬[23]應用顯式有限元法對波動方程進行離散化求解，研究了斷層臺階場地的反應特點，并指出了地表運動的差異性，有效考慮了不規(guī)則地形對地震波型的轉換問題，然而波場計算依賴于邊界條件且計算繁瑣；顧亮等[24]基于有限元數值模擬與透射邊界理論，研究了SV 波斜入射下陡坎地形地表峰值加速度的響應問題，文中對自由波場的計算較為簡便，但對局部地形的散射問題考慮不充分，因此精度不高；趙密等[25]基于黏性邊界條件推導了邊界上的等效節(jié)點力，研究了SV 波傾斜入射至階梯型場地時的地震動輸入方法，可將地表位移的峰值誤差控制在5%以內，但并未詳細介紹垂直入射與斜入射時地震波場計算的不同之處；BAO 等[26]提出了混合波場的思路，結合有限元法和波動輸入方法研究了SV 波垂直和傾斜入射下階梯型地形的地震動輸入方法，其精確度和效率有了較大的提升，然而地形條件和波型較為單一。以上幾類典型的方法不斷提升了計算精度與適用范圍，但是對局部地形的散射問題考慮還不夠充分，所采用的地震波波型、地震入射角度、場地地形條件等也有限，且在地震動輸入時普遍對人工邊界的限制條件較高，計算過程十分復雜；因此，亟需提出更為通用和高效的地震波場計算與輸入方法。

本文提出了一種不規(guī)則地形條件下的地震波輸入方法。依據地震波傳播規(guī)律進行波場分離，再將非外行波場轉換為人工邊界上的等效節(jié)點力形式以完成計算。該方法有效考慮了局部地形中散射波場的影響，計算過程明晰簡便，易于推廣至復雜場地條件。驗證部分選用了不同的地形條件，將計算結果分別與精確解、參考解、成熟方法解進行對比；并應用到多組斜入射P 波和SV 波作用下場地的振動的分析中。

1 方法原理

1.1 地震波場分離

不規(guī)則地形條件下場地的全域自由波場無法統(tǒng)一求解，需要在各個邊界上針對不同地形特征與輸入條件分開計算。李山有和廖振鵬[23]提出地震波的輸入應滿足非外行波的條件(非外行波即通過邊界向有限域內傳播的波和平行邊界的地震波)，并證明了地震波垂直入射情況下，底邊界的非外行波即入射波，而側邊界的非外行波可看作與遠場介質對應的自由波場。而本實驗進一步考慮到傾斜入射的情況，認為斜入射情況下(以左下角入射為例)，右邊界區(qū)域在初始階段便存在反射波場，因此非外行波場也只有入射波場。此外，杜修力和趙密[27]曾也指出，目前所有的人工邊界都只能模擬傳遞到無限域的外行波，而由結構物存在、不規(guī)則地形條件、局部地質等引發(fā)的散射問題，都需在邊界處進行波場分離。

因此，根據地震波入射、反射和散射等原理，首先在不規(guī)則地形場地的各個邊界面上進行波場分離。當地震波垂直入射時，在底面邊界上將波場分為入射波場和邊界外行場，左右邊界面上分離為自由波場和散射場；當地震波傾斜入射時(以從左下角入射為例)，在底面邊界和右面邊界將波場分為入射波場和邊界外行場，而左面界面上則分離為自由波場和散射場即可。進而輸入非外行波場，如圖1 所示：1) 垂直入射情況下在底部邊界輸入入射波場，兩側邊界輸入各自高度對應的自由波場；2) 傾斜入射情況下，在底部邊界和右側邊界輸入入射波場，在左邊界輸入自由波場。

圖1 不規(guī)則地形下各邊界地震波場輸入Fig.1 Input of seismic wave field at each boundary under irregular terrain

1.2 地震波場求解

入射波場和自由波場是決定場地響應的主要因素，需要分別確定；由局部地形引起的散射波場則通過有限元計算體現。

1) 入射波場：地震的入射波場求解可以直接通過幾何推導實現。如圖2 所示，利用人工邊界截斷形成一定的有限域空間，假設地震P 波從左下角(x0,y0)點進行入射且振動位移時程為uP(t)。由于P 波入射時質點振動方向與傳播方向一致，其入射波場的水平和豎向位移計算公式如下：

圖2 斜入射P 波作用下的地震波場Fig.2 Seismic wave field under obliquely incident P waves

式中：θ1為P 波入射的角度；Δt1即為P 波傳播到邊界上任一節(jié)點(x,y)的時間；cP為地震P 波的波速。

假設地震SV 波從左下角(x0,y0)點進行入射且振動位移時程為uSV(t)，由于SV 波入射時質點振動方向與傳播方向相互垂直，其入射波場的水平和豎向位移計算如下：

式中：θ2為SV 波入射的角度；Δt2即為SV 波傳播到邊界上任一節(jié)點(x,y)的時間；cSV為地震SV 波的波速。

2) 自由波場：如圖3 所示，自由波場包括入射場和反射場，例如P 波斜入射時又會在界面分層處同時產生反射P 波和反射SV 波，且幅值系數發(fā)生了變化，因此不能通過簡單的幾何推導求得。

圖3 斜入射地震波傳播規(guī)律Fig.3 Propagation law of obliquely incident seismic waves

目前在研究斜入射地震波下場地的自由波場時，主要有頻域和時域兩種方法，后者更適用于復雜的場地條件且與現代大型計算軟件匹配良好，因而廣泛使用。目前比較流行的有一維時域化方法[28]，其結合了集中質量有限元法和中心差分法以建立二維波動方程，計算出單一垂線上節(jié)點的位移，再通過幾何推導擴展到全域波場。

主要過程如下：圖4 為離散化后的半空間模型，橫向和豎向網格尺寸為Δx和Δy，地震P 波或SV 波從模型底部的人工邊界進行輸入。將y軸上的三類節(jié)點包括中間節(jié)點(0,n)、自由表面節(jié)點(0, 0)和人工邊界上的節(jié)點(0,N)在PΔt時刻的運動方程表示如下：

式中：M、C、K分別為節(jié)點的質量、阻尼和剛度矩陣，下標為節(jié)點的坐標信息；CB和FB分別為邊界的阻尼系數和需要輸入的集中力。在P 波和SV 波入射時分別有不同的FB計算方法，求解上述方程組，即可確定半空間內y軸上所有節(jié)點在(p+1)Δt時刻的位移，再結合幾何關系推廣到全域即可計算出整個半無限空間的自由波場。

1.3 地震波動輸入

地震屬于振動中的外源問題，在進行波場推求后還需進行地震動輸入才能準確計算場地的動力響應，人工邊界和地震動輸入方式的選擇是計算中的兩個關鍵環(huán)節(jié)[29]。數值模擬過程中通常在無限域地基中截取感興趣的部分作為有限域，并在截斷處賦予各類人工邊界，其主要作用是吸收外行波以達到對無限地基輻射阻尼的模擬。地震動輸入則主要分為振動法和波動法，波動法最早由劉晶波和呂彥東[30]于1998 年提出，可有效克服振動法中對于場地材料、激勵方式和邊界處理上的簡化造成的誤差，可以處理地震波斜入射問題和復雜場地的輸入問題，也因其時空解耦性的特點而與各類大型通用軟件適配性極高，是目前模擬精度最高且最為流行的地震動輸入方法。波動法的核心原理是將地震動輸入轉換為人工邊界上的等效荷載力，需要引入黏性邊界和黏彈性邊界等人工邊界；因此需要同時計算由自由場產生的應力和引入邊界后的附加應力，針對不同網格尺寸和方向不同的應力條件也需逐一判斷，計算過程十分復雜；尤其是考慮復雜場地條件時，等效荷載力計算的難度顯著增大，尤其是所要探討的不規(guī)則地形條件，傳統(tǒng)的波動方法已不具備高效處理的能力，亟需進行優(yōu)化。

劉晶波等[31]和寶鑫等[19]對波動法進行了改進，從一維土柱模型入手證明了在邊界節(jié)點附近的所有相鄰節(jié)點同時施加入射波場位移或自由波場位移，再將輸入位移后計算得到的荷載力反向輸入回模型中可完成地震波輸入，進而提出了基于人工邊界子結構的地震波動輸入方法，如圖5所示。圖5(a)為不規(guī)則地形的有限元模型；圖5(b)中將含有邊界節(jié)點的所有單元稱為邊界子結構，將地震波場位移時程輸入到子結構中可得到邊界節(jié)點的反力，經過證明此反力在數值上恰好等于需要的等效荷載力；圖5(c)表示重新在邊界節(jié)點上輸回此反力即可實現波動輸入。該方法保留了波動法的諸多優(yōu)勢并簡化了等效荷載力的計算流程，避免了復雜等效節(jié)點力的過程，可成功地將不規(guī)則地形場地的地震波輸入問題聚焦為地震波場求解問題。

圖5 基于人工邊界子結構的波動輸入方法Fig.5 Wave input method based on substructure of artificial boundaries

2 算例驗證

2.1 有限元模型

為了驗證所提出的方法在各類不規(guī)則地形條件下的適配性，分別以階梯型、V 字河谷型和梯形河谷型這三類工程中常見的不規(guī)則地形條件為例進行驗證。圖6 展示了三類地形有限元模型尺寸信息，在模型的左側、右側和底面邊界分別采用了黏彈性邊界條件進行約束。由于斜入射地震波場在不規(guī)則地形條件下的頻域解析解不能求得，可擴大截取的有限域場地范圍，采用遠置邊界的方法保證地震波不會通過人工邊界反射后再回到需要觀測的部位，并將遠置邊界法得到的結果作為有限元近似精確解(簡稱參考解)。圖7 展示了上述三類場地相對應的遠置邊界有限元模型信息，邊界均無任何約束。圖6 和圖7 中的規(guī)則部分為正方體網格，尺寸為1 m×1 m，A點、B點、C點為自由表面觀測點，利用觀測點的位移信息來進行對比驗證。三類場地的材料參數及地震波速信息見表1。

表1 場地模型參數Table 1 Model parameters of the site

圖6 不規(guī)則地形場地的有限元模型 /mFig.6 Finite element model of irregular terrain site

以平面P 波和SV 波為例，傾斜入射時假設地震波統(tǒng)一從場地的左下角進入，若考慮從右下角只需將左右邊界的波場分離方式互換即可，此處不再詳細展開。選用脈沖地震波的時程曲線及傅里葉譜曲線如圖8 所示。兩種場地模型參數條件下滿足一維時域化方法自由場計算精度的時間間隔分別為0.0029 s 和0.0020 s，統(tǒng)一選取0.002 s 作為地震波的輸入以及后續(xù)有限元計算的時間間隔。不規(guī)則地形條件下各個邊界輸入波場的計算方法參考第1 節(jié)；而在遠置邊界法中，各邊界均采用對應高度的自由波場即可。

圖8 入射脈沖波的位移時程和傅里葉譜Fig.8 Displacement and Fourier spectrum of the incident seismic wave

共進行了4 組驗證實驗，選取了入射角度為0°和30°的平面脈沖P 波和入射角度為0°和15°的平面脈沖SV 波。其中，P 波和SV 波垂直入射計算時選用階梯型場地，P 波30°入射選用梯形河谷型場地，SV 波15°入射時選用V 字河谷型場地，可充分驗證方法對不同條件下的適配性。波場計算過程均利用python 編程，有限元計算過程在大型通用軟件ABAQUS 中進行。

2.2 P 波入射情況驗證

如圖9 所示，將節(jié)點等效荷載力輸入到有限元模型的人工邊界后，得到了P 波垂直和傾斜入射下不規(guī)則場地的自由表面點A、點B、點C的橫向與縱向位移情況，圖例中的Ux和Uy分別代表觀測點橫向和縱向位移，為將結果顯示清楚模擬解部分的位移時間間隔擴大了15 倍?？梢钥闯霰疚姆椒ǖ慕Y果(模擬解)與參考解吻合很好。其中，需要注意的是由于遠置邊界模型尺寸和原有的不規(guī)則地形尺寸不同，在比較位移結果時需要去除前面一段時間(地震波從遠置邊界到原有邊界處的時間)。此外，觀察結果發(fā)現：1) 觀測點的豎向高度越高，其位移幅值越大；2) 雖然兩個案例的場地參數不同，但由于輸入的脈沖波位移一樣，可以在位移幅值上做出比較，P 波傾斜入射相較于垂直入射時，橫向位移幅值增大，豎向位移幅值略有減小。這兩點與實際情況相符。

圖9 P 波入射下自由表面點位移Fig.9 Displacement of free surface under incident P waves

2.3 SV 波入射情況驗證

圖10 為SV 波垂直和傾斜入射下自由表面點A、點B、點C的水平與豎向位移情況，結果同樣與參考解幾乎吻合，說明方法對P 波和SV 波都是通用的。結果發(fā)現：1) 觀測點的豎向高度越高，其位移幅值越大，規(guī)律同P 波相同；2) SV 波傾斜入射相較于垂直入射時，橫向位移幅值減小，豎向位移幅值增大。這兩點也符合實際規(guī)律。

圖10 SV 波入射下自由表面點位移Fig.10 Displacement of free surface under incident SV waves

2.4 結果與應用討論

2.4.1 驗證結果對比

表2 展示了本文方法與遠置邊界法的對比，表2 中的計算時長和文件均指第2 次節(jié)點力輸入計算；實驗使用的電腦型號為Intel(R) Core(TM) i7-8700 CPU @ 3.20 GHz，機帶RAM 為8 G，ABAQUS為2017 版本?？梢钥闯鍪褂帽疚姆椒稍诒ＷC高精度的情況下有效縮短計算時長并減小結果文件存儲空間，當研究范圍擴展到三維空間后優(yōu)勢將更加明顯。

表2 兩種方法的結果對比Table 2 Comparison of the two methods

2.4.2 應用分析

由于本文方法涉及的理論推導和步驟較多，可以通過每次計算后的關鍵點位移和波場快照等進行分步檢驗，確保其合理性后再繼續(xù)進行。主要流程分兩次輸入與計算，統(tǒng)一以P 波入射為例進行說明。

1) 非外行波場計算及第1 次位移輸入：圖11(a)為P 波垂直入射時采用波場分離技術計算后，三個邊界上典型位置點的輸入位移時程(具體位置如圖5，G1為底部邊界中點)。兩側自由表面的E1、F1點的橫向位移幅值分別為底部G1點的兩倍，這是因為自由波場包括了入射波場和反射波場，縱向位移均為0，結果是符合規(guī)律的。圖11(b)顯示P 波入射角度為30°時，左邊界E3點的橫、豎向位移均為F3點和G3點的兩倍，也是由于自由波場同時包括了入射場和反射場造成的位移幅值加倍；F3點的位移曲線同底部G3點變化趨勢一樣，但稍有延遲，這是由行波效應引起的，整體符合波場傳播的規(guī)律。

圖11 場地關鍵點的位移輸入情況Fig.11 Input displacement of key points in the site

2) 第2 次等效荷載力輸入：即第1 次輸入波場位移后進行有限元計算，提出各邊界上的等效節(jié)點力再反向輸回的過程。圖12 和圖13 為第2 次地震波輸入后的波場快照，可以看出：開始階段波場基本沿平面形狀進行入射，證明各個邊界等效節(jié)點力計算及輸入的相對時間關系準確；在傳播到不規(guī)則地表后波陣面形狀會發(fā)生變化，可看出局部地形的影響；最后隨著時間增加波場位移逐漸消失，說明黏彈性邊界有效吸收了散射波。因此，不需要將邊界進行遠置即可獲得較準確的結果。

圖12 P 波垂直入射時的波場快照Fig.12 Wave field snapshot of site under vertical incident P waves

圖13 P 波入射角度30°時的波場快照Fig.13 Wave field snapshot of site under P waves with the incident angle of 30°

2.5 方法對比驗證與拓展

2.5.1 與頻域精確解的對比

由于不規(guī)則地形的地震波場無頻域精確解，故前文結果均與參考解進行對比。實際上，本方法針對大部分半空間場地是通用的，為進一步證明其正確性與精度，建立規(guī)則地形模型進行驗證，規(guī)則地形地震波場的建立基于較多的假設與簡化，因而有精確的解析解。

參考了文獻[12]選取相應的解析方法，圖14為規(guī)則地形場地的有限元模型，網格劃分方式、場地模型參數、邊界約束、荷載施加方式和地震波形式均與2.1 節(jié)的模型一致。選取自由表面中點C點為觀測點，圖15 為P 波30°入射和SV 波垂直入射下C點的位移?？梢钥闯霰疚姆椒ㄋ贸龅慕Y果與精確解幾乎吻合，驗證了所提出方法的可靠性。

圖14 規(guī)則地形場地的有限元模型 /mFig.14 Finite element model of regular terrain site

圖15 地震波入射下規(guī)則場地自由表面點位移Fig.15 Displacement of free surface of regular site under seismic waves

2.5.2 與已有方法的對比

不規(guī)則地形場地的波場求解已有一些成熟的方法，選取文獻[25]中的方法進行進一步對比。文獻中將波場在人工邊界處進行分解，用黏性邊界模擬輸入波場外的其余波場，并根據地基參數確定邊界阻尼參數等。將此方法應用到圖6 中的階梯型場地與梯形河谷型場地模型中，分別選取了自由表面C1和C3點為觀測點，圖16 展示了其在P 波30°入射和SV 波垂直入射下的位移。

圖16 地震波入射下階梯型場地自由表面點位移Fig.16 Displacement of free surface of stepped site under seismic waves

從圖16 中可以看出，本方法和文獻[25]中方法得出的位移結果基本一致，證實了其在垂直入射和斜入射時波場計算均較為高效，同時本文方法相比于成熟方法而言計算節(jié)點力過程更加簡單，對人工邊界的限制也較低。

2.5.3 拓展應用

波場分離方法同樣適用復雜的場地條件，只是當介質參數較為復雜時，入射波場和自由波場的求解難度會顯著增大。為使方法推廣性更強，這里以最為常見的成層場地為例，補充一下其計算思路。

成層介質自由波場的推求已有較多頻域和時域的成熟方法[32-33]。而大部分文獻中沒有詳細介紹成層地基內入射波場求解，這里推薦2 種方案：1) 頻域法求解可結合傳遞函數法和剛度矩陣而推求；2)時域方法求解可根據幾何關系[34-36]，截取地震波還未入射到自由表面前的位移時程，即可從自由波場中去除反射波場的影響，得到入射波場。

此外，更為復雜的三維真實地形條件下地震波場的構建也是后續(xù)研究的重點。

3 響應結果分析

3.1 入射角度對結果的影響

以圖5(a)中的階梯型模型為例，分別計算P 波和SV 波不同入射角度下自由表面A1點的位移響應情況。因為平面SV 波入射時存在臨界入射角，一般在35°左右，因此選取SV 波入射角度范圍為0°～30°，P 波入射角度為0°～90°。如圖17 和圖18所示，隨著P 波入射角度的增大，橫向位移幅值先增大后減小，大約在60°左右時達到最大；豎向位移呈現逐漸減小的趨勢。但當P 波入射角度為90°時，位移呈現波動變化趨勢并未完全收斂，說明本文方法對于超大角度的波動輸入問題仍有一定問題需要處理。隨著SV 波入射角度的增大，橫向位移呈現逐漸減小的趨勢但變化幅度較小，縱向位移幅值呈現逐漸增大的趨勢。對于不規(guī)則地形來說，入射角度對振動響應的影響也是符合基本變化規(guī)律，且不同角度情況下位移變化幅度較大，說明對于地震波入射角度的考慮是必要的。

圖17 P 波不同入射角度下階梯型場地A1 點位移Fig.17 Displacement of A1 on the stepped site under P waves with different incident angles

圖18 SV 波不同入射角度下階梯型場地A1 點位移Fig.18 Displacement of A1 on the stepped site under SV waves with different incident angles

3.2 局部地形對結果的影響

選用圖14 中的規(guī)則地形場地作為參考，即100 m×50 m 的長方形場地，其動力參數均取為相同。表3 和表4 分別展示了P 波和SV 波入射條件下兩種地形條件下自由表面A1點的位移幅值對比；其中，相差百分比的計算方法統(tǒng)一為(規(guī)則地形下位移-階梯型地形下位移)/階梯型地形下位移。由表3 和表4 可見，在兩種場地條件下，位移隨著地震波入射角度的變化均呈現相同的變化規(guī)律，與3.1 節(jié)的規(guī)律一致。同時可以看出局部地形條件對場地位移的影響較大，在P 波入射條件下大部分情況階梯型場地下位移較大，SV 波入射下則相反。說明地形對地震波散射等傳播規(guī)律的影響復雜多變，會隨著波型和入射角度的變化而變化，在地震計算中對于地形因素的考慮也是必要的。

表3 P 波入射時不同地形條件下自由表面A1 點位移Table 3 Displacement of A1 under different terrain conditions when P waves incident

表4 SV 波入射時不同地形條件下自由表面A1 點位移Table 4 Displacement of A1 under different terrain conditions when SV waves incident

4 結論

合理的地震動輸入是動力分析的必要基礎，采取了解析推導結合數值模擬的手段，提出了一種基于波場分離技術和應力型人工邊界的地震波場計算與輸入方法，分別詳細推導了地震P 波和SV 波在垂直入射和傾斜入射下的地震波場，在充分考慮到局部地形產生的散射波場的同時，采用了改進的波動輸入方法實現地震波的高效輸入。此外，對比分析了多組不同入射角度下的規(guī)則和不規(guī)則地形的場地反應情況，主要得到了以下結論：

(1) 通過與遠置邊界法得到的自由表面位移的參考解進行對比，證明了方法的精確度與高效性，且適用于多類不規(guī)則地形條件的場地以及不同的地震波型和入射角度。

(2) 所采用的改進的波動輸入法，極大簡化了等效節(jié)點力計算，將地震波輸入的重難點聚焦于場地自由波場和入射波場的計算，使得所提出的方法易于推廣至三維場地以及成層地基情況。

(3) 地震波的入射角度變化對不規(guī)則場地地表位移響應影響較大，自由表面點隨著P 波入射角度的增大，橫向位移幅值先增大后減小，豎向位移呈現逐漸減小的趨勢；隨著SV 波入射角度增大，橫向位移呈現逐漸減小的趨勢，縱向位移規(guī)律相反。整體變換規(guī)律與規(guī)則場地條件下的相一致且符合實際。

(4) 不規(guī)則地形條件對場地的地表位移響應也有較大影響，地震波通過局部地形時傳播路徑會發(fā)生變化并產生額外的散射波場，因此，其振動反應相較于規(guī)則場地有明顯不同，并同時受到地震波型和入射角度等的影響。

所提出的方法和得出的規(guī)律性結論可為不規(guī)則地形上的振動響應分析提供手段，后續(xù)研究將進一步聚焦于三維真實地形條件下地震波場的構建與輸入方法改進提升。