含峭壁V形峽谷對地震SH波散射的解析解

張寧， 潘家瑣， 代登輝， 高玉峰*

1 河海大學巖土力學與堤壩工程教育部重點實驗室， 南京 210098 2 河海大學土木與交通學院， 南京 210098

0 引言

當地震波遇到峽谷、高山等局部地形時會發生散射，散射波、入射波和反射波的相互作用會引起地震動的放大或衰減，稱為地形效應(周紅等，2010；田林等，2012).自從美國圣費爾南多地震記錄的異常放大現象(Trifunac and Hudson, 1971)得到關注之后，國內外大量的地震實測記錄和現場震害觀測進一步證實了地形效應的存在(王海云和謝禮立,2010；Hough et al., 2010；李平等, 2016).

為了揭示地形效應機理，地震學和土木工程領域學者開展了地表地形對地震地面運動影響規律的理論研究，包括數值模擬和解析分析兩類方法.數值研究眾多，例如有限差分法(Boore, 1972)、有限元法(廖振鵬和劉晶波,1992；章小龍等, 2017)、譜元法(李孝波等，2014；于彥彥等,2017；賀春暉等,2017)、邊界方法(林皋和關飛,1990；曹軍等,2004；巴振寧等,2017；Liu et al., 2018)及混合方法(Zhang and Zhao,1988；金峰等,1993；Zhou and Chen, 2006).數值方法可以模擬任意復雜地形情況，但其計算精度往往需要解析解的校準.解析方法主要是波函數展開法，Trifunac(1973)提出的半圓形峽谷對地震SH波散射的波函數級數解成為早期經典的解析解，隨后涌現出更多針對不同地形的解析解，如半橢圓形峽谷(Wong and Trifunac, 1974)、圓弧形峽谷(Yuan and Liao, 1994)、圓弧形飽和土峽谷(李偉華和趙成剛,2003)、圓弧形層狀峽谷(梁建文等,2003)、圓弧形層狀沉積谷(張郁山,2008，2010)、去底半圓形峽谷(Tsaur and Chang, 2009)、U形峽谷(Gao et al., 2012)、徑向非均質半圓形峽谷(Zhang et al., 2017, 2019)等.這些解析解為數值方法提供了驗證基準，并且引領了新型數值方法(Dai et al., 2019)的發展.

圖1 含峭壁V形峽谷示例：亞利桑那州北部Glen峽谷Fig.1 An example of V-shaped canyon with cliffs: Glen canyon in northern Arizona

由于V形峽谷分布廣泛，這類地形對地震波的散射問題得到了大量的數值(Sánchez-Sesma and Rosenblueth, 1979；趙崇斌等,1988)和解析(Tsaur and Chang, 2008；Tsaur et al., 2010；Zhang et al., 2012a，b；Gao and Zhang, 2013)研究.在長期的地質演變中，由于峽谷側壁的地應力集中，這種峽谷最終傾向于演變成含峭壁V形峽谷，典型的例子如亞利桑那州北部科羅拉多河的Glen峽谷(如圖1).然而，含峭壁V形峽谷地形效應研究尚未見報道.實際上，諸多峽谷具有峭壁，有必要開展含峭壁峽谷地形對地震波的散射特性研究.因此，本文采用波函數展開法提出了含峭壁V形峽谷對平面SH波散射的解析解，重點考察了上部峭壁對其地形效應的影響.

1 模型和理論推導

圖2 平面SH波作用下含峭壁V形峽谷的二維模型Fig.2 2D model of a V-shaped canyon with cliffs subjected to a plane SH waves

對于穩態反平面問題，可省略時間因子e-iω t，則模型兩個區域的穩態波場uj滿足亥姆霍茲方程(Pao and Mow, 1973)：

(1)

其中j=1,2,分別表示區域①和②，k=ω/cs表示波數.

除了亥姆霍茲方程，兩個區域的波場uj還需要滿足水平地表、峽谷表面和峭壁表面的應力自由條件：

(2)

(3)

(4)

另外，根據連續介質假設，由虛擬邊界S1分割而成的兩區域的位移和應力均須滿足連續性條件：

u1(r,θ)=u2(r,θ),r=b, 0≤|θ|≤β，

(5)

(6)

公式(5)和(6)可保證兩區域之間的位移和應力匹配.

區域①的波場包含兩部分(自由場和散射場)：

u1=uf+us,

(7)

自由場可以表示為

uf(r,θ)=exp[-ikrcos(θ+α)]

+exp[ikrcos(θ-α)]，

(8)

利用Jacobi-Anger展開公式：

(9)

其中εn是紐曼因子(ε0=1，εn=2,n≥1)，Jn(·)是n階第一類貝塞函數，公式(8)中的自由場uf可以展開成：

(10)

模型物理意義明確，散射波由河谷處產生、向無窮遠處傳播，為了得到散射波場的唯一解，其應該滿足無窮遠處Sommerfeld輻射條件：

(11)

通過分離變量法求解公式(1)并排除不滿足輻射條件(11)的波場后可以得到：

(12)

在區域②中，滿足峽谷底部應力自由條件(公式(3))的駐波場u2可以表示為

(13)

其中Cn和Dn為待求系數，ν=π/(2β1).

為了便于問題的解決，需要用Graf加法公式將局部坐標系(r1,θ1)中的駐波場u2轉化到整體坐標系(r,θ)中.推導得到相應的變換公式為

(14)

其中h=h1+h2.

將公式(14)代入公式(13)，駐波場u2可以表示為

+1)ν+m]θ.

(15)

至此，利用前述區域分解策略完成了兩區域波場的構建，且兩區域波場分別自動滿足水平地表應力自由條件(公式(2))和峽谷表面應力自由條件(公式(3)).然而，波場的四組未知系數An、Bn、Cn和Dn仍然待定，需要利用余下的峭壁應力自由條件(公式(4))和虛擬邊界連續條件(公式(5)和(6))完成求解，即區域匹配技術：

首先，補充區域②在峭壁上的零應力條件，獲得區域②在與區域①的半圓形交界面上的分段應力函數：

(16)

據此可將峭壁應力自由條件(公式(4))和虛擬邊界應力連續條件(公式(6))合并為

(17)

對式(17)兩邊同時關于θ在區間[-π/2,π/2]積分，并利用三角函數在此區間的正交性，整理可得關系式：

(18)

(19)

(20)

(21)

然后，利用虛擬邊界S1上的位移連續條件(公式(5))，將兩區域位移場在區間[-β,β]上積分得：

(22)

(23)

最后，將公式(18)和(19)分別代入公式(22)和(23)，得到下面僅和Cn、Dn有關的兩個方程組：

(24)

(25)

將各無窮級數進行截斷、保留有限項之后即可進行未知系數計算(n、m和p分別截取N、M和P項).求解公式(24)和(25)可得到系數Cn和Dn，隨后通過公式(18)和(19)求得An和Bn.得到這些未知系數后，半空間中的波場即已知.

關于N、M和P的取值，其中N代表級數解的級數項數量，需要通過級數解收斂測試確定.本文選取了峽谷表面5個代表性位置進行了收斂測試.圖3給出了參數對應α=30°、h1/b=h2/b=1/3和η=4的收斂測試結果.從圖3中可以看出，隨著N的逐漸增大，地表位移逐漸收斂至一個固定值.通過多頻率下的收斂測試，本文采用N=25即可獲得精確的結果.M和P要保證駐波場u2和自由場uf的坐標變換精度，分別利用公式(14)和(9)的滿足情況確定其數值為M=150，P=150.

2 驗證

為了方便起見，首先定義本文位移幅值|u|和無量綱頻率η：

(26)

(27)

其中Re(·)和Im(·)分別表示復位移的實部和虛部，λ表示入射SH波波長.

當h1/b=0時，含峭壁的V形峽谷可以退化為V形峽谷；當h2/b=0時，可退化為去底半圓形峽谷.Tsaur和Chang(2008，2009)分別給出了V形峽谷和去底半圓形峽谷的解析解，因此可以用來驗證本文結果的準確性.

圖4將本文的結果與Tsaur和Chang(2008)進行了對比，對應參數為h1/b=0，h2/b=1/2，η=1.圖5將本文的結果與Tsaur和Chang(2009)進行了對比，參數為h1/b=1/2，h2/b=0，η=1.通過4個不同入射角度α=0°、30°、60°和90°下結果的良好吻合情況，驗證了本文方法的準確性.

圖3 峽谷表面五個代表性位置的地表位移收斂測試圖(α=30°, h1/b=h2/b=1/3, η=4)Fig.3 Convergence tests at five representative positions on the canyon surface (α=30°, h1/b=h2/b=1/3, η=4)

圖4 本文結果(實線，對應參數：h2/b=0，h2/b=1/2，η=1)與Tsaur和Chang(2008)中無峭壁V形峽谷相應結果(點線)的比較(a) α=0°; (b) α=30°; (c) α=60°; (d) α=90°.Fig.4 Comparison between the solution of this study (solid lines) when h1/b=0, h2/b=1/2, η=1 and the corresponding results (dotted lines) of the V-shaped canyon without steep cliffs given by Tsaur and Chang (2008)

為了進一步檢驗本方法的正確性和實用性，繼續采用退化V形峽谷模型模擬中國臺灣翡翠峽谷臺陣的地震動并與實測結果進行對照.根據Huang和Chiu(1995)對1992年花蓮地震時翡翠河谷臺陣實測地震動的分析，編號為SC1的臺站地震記錄相比SC4加速度峰值放大了2.69倍.

翡翠臺陣所在的河谷剖面如圖6a所示，利用本文方法為之建立V形峽谷模型(b=563 m，d=300 m，h1=0 m).由于1992年花蓮地震入射波數據缺失，本文采用PGA=0.2g的集集地震波作為入射波，利用本文模型求得臺站SC1和SC4處的地震動頻域傳遞函數，然后據此得到臺站SC1和SC4的地震動加速度時程，如圖6b所示(頻域轉時域方法參見Dai et al.，2019).其中臺站SC1和SC4處的模擬地震動PGA分別為0.3324 g和0.1256 g，即SC1相對于SC4放大2.65倍，與前述實測2.69倍基本一致.可見，本方法可以較好地刻畫峽谷對地震動的差異放大規律.

圖5 本文結果(實線，對應參數h2/b=1/2，h2/b=0，η=1)與Tsaur和Chang(2009)中去底半圓形峽谷的相應結果(虛線)的比較(a) α=0°; (b) α=30°; (c) α=60°; (d) α=90°.Fig.5 Comparison between the solution of this study (solid lines) when h1/b=1/2, h2/b=0, η=1 and the corresponding results (dotted lines) of truncated semicircular canyon given by Tsaur and Chang (2009)

圖6 本文方法對中國臺灣翡翠河谷地震動的模擬(a) 翡翠河谷模型簡圖(據Huang and Chiu, 1995)； (b) 集集地震入射下SC1臺站和SC4臺站SH波地震動加速度時程模擬結果.Fig.6 Simulation of the ground motions for the Feitsui Canyon, Taiwan, China(a) Schematic diagram of Feitsui Canyon model (modified from Huang and Chiu, 1995)； (b) The calculated SH ground accelerations at positions SC1 and SC4 under the incidence of the Chi-Chi earthquake.

3 數值結果與討論

本文聚焦上部峭壁對峽谷地形效應的影響，故將參數h2/b設為固定值1/3，從而使得h1/b有較大的變化范圍.圖7—圖9分別為對應3個入射波無量綱頻率η=0.5、1和4的峽谷地表位移幅值結果，每個無量綱頻率結果都包含4個不同入射角度(α=0°,30°,60°,90°)，實線結果對應h1/b=0，虛線對應h1/b=1/3，虛實線對應h1/b=2/3.為了方便比較，x/b的計算范圍選取為-4到4.

圖7 對應η=0.5和不同入射角度的地表位移幅值結果(a) α=0°; (b) α=30°; (c) α=60°; (d) α=90°.黑色實線對應h1/b=0，h2/b=1/3；紅色虛線對應h1/b=1/3，h2/b=1/3；藍色點劃線對應h1/b=2/3，h2/b=1/3.Fig.7 Surface displacement amplitudes at four different incident angles and η=0.5(a) α=0°; (b) α=30°; (c) α=60°; (d) α=90°. h1/b=0, h2/b=1/3 (black solid lines), h1/b=1/3, h2/b=1/3 (red dashed lines), h1/b=2/3, and h2/b=1/3 (blue dotted lines).

圖8 對應η=1和不同入射角度的地表位移幅值結果(a) α=0°; (b) α=30°; (c) α=60°; (d) α=90°.黑色實線對應h1/b=0，h2/b=1/3；紅色虛線對應h1/b=1/3，h2/b=1/3；藍色點劃線對應h1/b=2/3，h2/b=1/3.Fig.8 Surface displacement amplitudes at four different incident angles and η=1(a) α=0°; (b) α=30°; (c) α=60°; (d) α=90°. h1/b=0, h2/b=1/3 (black solid lines), h1/b=1/3, h2/b=1/3 (red dashed lines), h1/b=2/3, and h2/b=1/3 (blue dotted lines).

圖9 對應η=4和不同入射角度的地表位移幅值結果(a) α=0°; (b) α=30°; (c) α=60°; (d) α=90°.黑色實線對應h1/b=0，h2/b=1/3；紅色虛線對應h1/b=1/3，h2/b=1/3；藍色點劃線對應h1/b=2/3，h2/b=1/3.Fig.9 Surface displacement amplitudes at four different incident angles and η=4(a) α=0°; (b) α=30°; (c) α=60°; (d) α=90°. h1/b=0, h2/b=1/3 (black solid lines), h1/b=1/3, h2/b=1/3 (red dashed lines), h1/b=2/3, and h2/b=1/3 (blue dotted lines).

從圖7—圖9的結果可以看出，含峭壁V形峽谷的地表位移幅值在自由場位移幅值2上下波動(|u|>2表示放大，|u|<2表示減小)，峽谷峭壁對地表位移幅值具有重要影響.一方面，峭壁使得地震放大效應更為顯著，通過比較h1/b=2/3和h1/b=0兩種情況，可以發現含峭壁峽谷的地表位移幅值可達不含峭壁峽谷的2倍.另一方面，隨著峭壁深度h1/b從1/3增大至2/3，位移振幅的最大值增大，而位移振幅的最小值減小.總之，峭壁對V形峽谷地形效應具有顯著的增強作用，這意味著在地震作用下峭壁峽谷周邊的建筑物將會遭受更為劇烈的振動，可能會造成更大的震害.

對于地震波豎向入射的情況(如圖7a、圖8a和圖9a所示)，由于波的聚焦，地表位移幅值的最大值往往發生在峽谷兩肩.不同于具有對稱性的豎向入射結果，地震波斜入射的情況(圖7(b—d)、圖8(b—d)和圖9(b—d))展示出左右不對稱的地表位移幅值.迎波側比背波側的位移幅值波動得更加頻繁，隨著入射角度的增加，這種不對稱表現得更加明顯.在小角度斜入射(α=30°)的情況下，迎波側的位移振幅僅略大于背波側，而對于大角度斜入射(α=60°, 90°)的情況，兩側差別更加明顯.隨著無量綱頻率η的增大，位移振幅的放大和減小越來越明顯，說明峽谷成為一個有效的障礙物，短波穿越峽谷到達對岸的能力越來越弱.在圖7(b—d)、圖8(c—d)和圖9d中，陰影區的位移幅值曲線趨于平穩，即體現了峭壁峽谷的濾波屏蔽效應.

圖10展示了在三種入射角度下V形峽谷地表位移幅值隨著無量綱位置x/b和無量綱頻率η的變化情況，其中圖10(a—c)是不含峭壁的情況(h1/b=0)、圖10(d—f)是含峭壁的情況(h1/b=1/3).圖中可見，隨著無量綱頻率的增大，含峭壁峽谷的地形放大形態趨于復雜.在同一入射角下，峭壁的存在使得峽谷兩側表現出更為明顯的放大效應，也就是說峭壁的出現對峽谷地形效應有加劇作用.

圖10 不同角度下含峭壁和不含峭壁峽谷地表位移幅值隨位置x/b和無量綱頻率η的變化圖(a) α=0°，h1=0，h2=1/3; (b) α=45°，h1=0，h2=1/3; (c) α=90°， h1=0，h2=1/3; (d) α=0°，h1=1/3，h2=1/3; (e) α=45°，h1=1/3，h2=1/3; (f) α=90°，h1=1/3，h2=1/3.Fig.10 Ground surface displacement amplitudes as a function of position x/b and dimensionless frequency η for canyons with and without steep cliffs at four different incident angles

4 結論

本文提出了含峭壁V形峽谷對地震SH波散射的波函數級數解，揭示了含峭壁V形峽谷對地震動的地形放大效應.本文解析模型可退化為不含峭壁的V形峽谷和去底半圓形峽谷兩個已有解析模型，通過退化計算驗證了本文方法.計算結果表明上部峭壁對峽谷地震動具有重要影響，峭壁峽谷(h1/b=2/3)較無峭壁峽谷(h1/b=0)對地表位移幅值的放大差異可達200%，且較陡的峭壁明顯加劇峽谷的地震放大效應，也就增加了周圍建筑物的震害風險，值得工程師的關注.本文結論基于二維峽谷場地SH波散射模型，關于峽谷場地的三維效應以及P波和SV波等其他類型地震波的散射特征，尚需進一步探索.