隨機(jī)有限斷層法的俯沖帶板內(nèi)地震動(dòng)模擬

林德昕,馬 強(qiáng),陶冬旺,馬完君,解全才,劉名吉

(1.中國地震局工程力學(xué)研究所地震工程與工程振動(dòng)重點(diǎn)實(shí)驗(yàn)室,哈爾濱 150080;2.地震災(zāi)害防治應(yīng)急管理部重點(diǎn)實(shí)驗(yàn)室,哈爾濱 150080)

隨機(jī)有限斷層方法考慮了震源、傳播路徑和場地的全過程,使用簡單高效、輸入?yún)?shù)物理意義清晰,是地震工程中主流的強(qiáng)地震動(dòng)模擬方法。其模擬的地震動(dòng)記錄在頻譜、幅值以及持時(shí)上與觀測記錄對應(yīng)良好[1]。通過模擬已發(fā)生地震的地震動(dòng)可以驗(yàn)證相對穩(wěn)定的路徑和場地參數(shù),對于未來可能發(fā)生的地震利用現(xiàn)有合理輸入?yún)?shù)結(jié)合其震源信息,使用模擬方法快速得到地震動(dòng)影響場,進(jìn)而為后續(xù)的地震烈度速報(bào)、震后情景構(gòu)建和災(zāi)情評估等震后態(tài)勢感知系統(tǒng)提供地震動(dòng)輸入。

該方法被廣泛應(yīng)用于世界各地的地震動(dòng)模擬中[2-5],但這些地震多是發(fā)生在內(nèi)陸或震源深度不足25 km的俯沖帶淺層地殼區(qū)域,目前國內(nèi)針對淺層地殼以下俯沖帶地震引起的地面運(yùn)動(dòng)模擬研究較少。文獻(xiàn)[6]對俯沖帶板間地震的地震動(dòng)進(jìn)行了模擬,雖然淺層地殼地震、俯沖帶板間地震和板內(nèi)地震存在一定的關(guān)聯(lián),但俯沖帶板內(nèi)(subduction slab)事件的震源深度更深(通常超過50 km),其地震動(dòng)大于同條件俯沖帶板間事件的地震動(dòng),地震動(dòng)峰值衰減慢于淺層地殼地震,短周期譜相比于另外兩者會(huì)更高[7-10]?？紤]到俯沖帶板內(nèi)和板間地震的地震動(dòng)差異較大[11],隨機(jī)有限斷層法在俯沖帶板內(nèi)地震動(dòng)模擬中的適用性還有必要進(jìn)一步驗(yàn)證。作為地震多發(fā)國家,中國海域分布在大陸板塊與海洋板塊共同作用的俯沖帶地區(qū),南海北部及華南沿海歷史上曾發(fā)生7級以上地震,俯沖帶地震活動(dòng)性較強(qiáng)[12]。然而中國對俯沖帶地震的相關(guān)研究起步較晚,海域地震動(dòng)記錄稀缺,這對地震風(fēng)險(xiǎn)評估造成嚴(yán)重制約,由于中國臺(tái)灣地區(qū)俯沖帶與日本俯沖帶構(gòu)造相似[13],日本地區(qū)俯沖帶地震研究對中國俯沖帶地震的地震動(dòng)模擬以及地震工程研究有重要的參考價(jià)值。

北京時(shí)間2021年10月7日21時(shí)41分,日本千葉縣(35.6°N,140.1°E)發(fā)生Mj6.1地震,日本氣象廳及日本總務(wù)省消防廳發(fā)布的地震震源深度為60 km。根據(jù)美國地質(zhì)調(diào)查局(USGS)公開的震源機(jī)制解,此次地震斷層機(jī)制為逆斷層。本文利用Zhao等[14]的日本地區(qū)地震分類方法,進(jìn)一步將其確定為俯沖帶板內(nèi)地震,基于隨機(jī)有限斷層法對此次俯沖帶板內(nèi)地震引起的地面運(yùn)動(dòng)進(jìn)行了模擬。

1 隨機(jī)有限斷層法

文獻(xiàn)[15]在隨機(jī)點(diǎn)源法[16-17]和靜拐角頻率的有限斷層方法[18-19]基礎(chǔ)上,提出了動(dòng)力學(xué)拐角頻率模型的隨機(jī)有限斷層方法,擺脫了子斷層尺寸劃分對地震動(dòng)模擬結(jié)果的影響,解決了同一個(gè)子斷層在地震中多次破裂的異常問題。此后為了保證隨機(jī)點(diǎn)源法和有限斷層法在遠(yuǎn)場模擬結(jié)果相匹配,文獻(xiàn)[20]進(jìn)一步做出重要改進(jìn),將子斷層拐角頻率的倒數(shù)作為震源上升持時(shí),在低頻段使用濾波器函數(shù)以保證低頻傅里葉譜的一致性,這使得隨機(jī)有限斷層法進(jìn)一步擴(kuò)大了適用震級范圍。

根據(jù)隨機(jī)有限斷層模型,將斷層面沿走向和傾向劃分為N個(gè)子斷層[21],頻域上可將第i行j列子斷層在距離為Rij的目標(biāo)場點(diǎn)加速度傅里葉譜描述為

(1)

式中:第1項(xiàng)為震源部分;M0ij、f0ij和R0ij分別為第i行j列子斷層的地震矩、拐角頻率和到目標(biāo)場點(diǎn)的震源距;C=RθφFV/(4πρβ3)為常數(shù),Rθφ為輻射圖型系數(shù),通常取剪切波的平均值0.55;F為自由表面放大,通常取2.0;V為2個(gè)水平分量的分配系數(shù),一般取0.71;ρ和β分別為震源附近的介質(zhì)密度和剪切波速。第i行j列子斷層的拐角頻率表示為

(2)

式中:Δσ為應(yīng)力降;M0為整個(gè)斷層的地震矩;N(t)為t時(shí)刻已經(jīng)破裂的子斷層的數(shù)量之和。

為了保證高頻輻射能守恒,文獻(xiàn)[15]引入高頻標(biāo)定因子Hij對低估的子斷層遠(yuǎn)場輻射能進(jìn)行補(bǔ)償,表達(dá)式如下:

(3)

將式(1)表示的所有子斷層在目標(biāo)場點(diǎn)產(chǎn)生的加速度譜經(jīng)過傅里葉逆變換,考慮斷層內(nèi)的破裂和地震波傳播造成的時(shí)間上的遲滯,在時(shí)域內(nèi)疊加即可得到該場點(diǎn)最終的地震動(dòng)加速度時(shí)程,可表示為

(4)

式中:nl和nw分別為沿?cái)鄬幼呦蚝蛢A向劃分的子斷層數(shù)目之和;aij為第i行j列子斷層的加速度時(shí)程,即式(1)中Aij的傅里葉逆變換值;Δtij為破裂傳播到第i行j列子斷層,并從該子斷層傳播到目標(biāo)場點(diǎn)的時(shí)間延遲。

2 數(shù)據(jù)處理和參數(shù)獲取

2.1 地震動(dòng)記錄篩選及處理

本次千葉縣Mj6.1俯沖帶板內(nèi)地震引起的地面運(yùn)動(dòng)模擬所用井下和井上加速度數(shù)據(jù)均來自日本防災(zāi)技術(shù)科學(xué)研究所的KiK-net強(qiáng)地震動(dòng)臺(tái)網(wǎng)。選取震中距≤100 km的臺(tái)站,利用文獻(xiàn)[22]基于信噪比的評估方法對記錄進(jìn)行篩選以確定記錄質(zhì)量,最終有25個(gè)臺(tái)站的記錄通過了閾值為5的篩選標(biāo)準(zhǔn),可被視為高質(zhì)量記錄,參考臺(tái)站及震中位置見圖1。

根據(jù)KiK-net臺(tái)網(wǎng)發(fā)布的強(qiáng)震動(dòng)記錄信息和臺(tái)站鉆孔資料,將所用臺(tái)站的震源距、井上水平向峰值加速度(Apg)和井上臺(tái)場地分類整理見表1,其中場地分類按照文獻(xiàn)[23]基于場地周期的方法進(jìn)行劃分,對于給出完整鉆孔資料的臺(tái)站,使用土層厚度與對應(yīng)土層等效剪切波速之比計(jì)算出的特征周期(即T=4H/V)確定場地類別;未給出鉆孔資料的臺(tái)站,使用H/V平均反應(yīng)譜比的峰值特征大致判斷場地類別。

表1 臺(tái)站信息

由表1可以看出,KNGH10井上臺(tái)站同時(shí)記錄到了本次地震EW向和NS向地表Apg值,分別為140.6 cm·s-2和107.8 cm·s-2,地表水平向Apg>60 cm·s-2的臺(tái)站有4個(gè),分別為CHBH10、CHBH12、CHBH16和KNGH10臺(tái)站;井上臺(tái)場地情況復(fù)雜,涵蓋場地分類標(biāo)準(zhǔn)的全部類型。所有地震動(dòng)數(shù)據(jù)經(jīng)過基線校正后,使用頻帶范圍為0.1～25 Hz的四階巴通沃斯非因果濾波器進(jìn)行帶通濾波處理。

2.2 輸入?yún)?shù)的計(jì)算和選取

本文將USGS提供的震源機(jī)制信息作為本次模擬的震源參數(shù),矩震級Mw=5.9,震源深度取60.5 km,斷層走向、傾角和滑動(dòng)角取178°、33°和94°。在斷層幾何信息不明確的情況下,斷層尺寸使用文獻(xiàn)[24]的俯沖帶板內(nèi)地震震源尺寸與矩震級經(jīng)驗(yàn)統(tǒng)計(jì)關(guān)系估算,將其確定為10 km×8 km。結(jié)合上述斷層尺寸估算出的上緣埋深為58.5 km。應(yīng)力降利用文獻(xiàn)[25]給出的破裂面積、應(yīng)力降與地震矩的經(jīng)驗(yàn)關(guān)系估計(jì)出初始值為3 MPa,之后通過“試錯(cuò)法”最佳匹配擬反應(yīng)譜(Aps)殘差平均值以確定本次模擬使用的應(yīng)力降[26],增量設(shè)置為0.5 MPa,最終將應(yīng)力降取為6 MPa。本文使用隨機(jī)滑動(dòng)分布的方式對子斷層滑動(dòng)權(quán)重進(jìn)行分配。震源附近介質(zhì)密度和S波波速由Crust 1.0全球地殼模型確定,分別為3.17 g·cm-3和4.29 km·s-1。

對于地震動(dòng)持時(shí),文獻(xiàn)[4]使用90%能量持時(shí)作為地震動(dòng)持時(shí)模型,這使得模擬加速度記錄在持時(shí)上更加接近觀測記錄,本文認(rèn)為70%能量持時(shí)更能體現(xiàn)S波主導(dǎo)的強(qiáng)地震動(dòng)持續(xù)時(shí)間。因此,取地震動(dòng)能量處于5%和75%之間的持續(xù)時(shí)間(即70%能量持時(shí))作為本次模擬的地震動(dòng)持時(shí)模型,以距離為自變量,回歸得到二者的經(jīng)驗(yàn)關(guān)系式

T=a+b·R

(5)

式中:T可被視為包含了震源持時(shí)和路徑持時(shí)的地震動(dòng)持時(shí)模型(以下記為70%能量模型),a、b為模型系數(shù),R為震源距。由于場地對于地震動(dòng)持時(shí)的影響尚不明確,以上模型可能包含場地項(xiàng)對于持時(shí)的影響。本文計(jì)算了以上所選臺(tái)站記錄的能量持時(shí),去掉其中持時(shí)異常過大和過小的臺(tái)站(CHBH20和KNGN20)記錄,對篩選后的記錄使用最小二乘法進(jìn)行回歸,結(jié)果見圖2。圓點(diǎn)表示觀測記錄的能量持時(shí),點(diǎn)劃線表示地震動(dòng)持時(shí)模型,對應(yīng)井下和井上地震動(dòng)持時(shí)模型函數(shù)分別為Tdown=0.129 7R+0.797 95和Tup=0.090 8R+4.930 64。

圖2 地震動(dòng)持時(shí)模型

井下臺(tái)的場地效應(yīng)考慮了包括地殼放大和表示局部場地高頻快速衰減的κ0值,在此基礎(chǔ)上,井上臺(tái)還考慮了局部場地效應(yīng),其中包含局部土層放大效應(yīng)以及對井上臺(tái)κ0值的修正。

近地表地殼放大通過改進(jìn)后的四分之一波長法[28]計(jì)算各個(gè)臺(tái)站所在地區(qū)的波速模型獲得,KiK-net提供了部分臺(tái)站對應(yīng)井下土層結(jié)構(gòu)的數(shù)據(jù)資料,但不足以使用四分之一波長法計(jì)算全部臺(tái)站的近地表地殼放大。對于沒有提供井下波速結(jié)構(gòu)的臺(tái)站,本文結(jié)合文獻(xiàn)[29]提出的日本綜合速度模型作為補(bǔ)充,選取臺(tái)站所在劃分區(qū)塊的地下結(jié)構(gòu)速度模型,計(jì)算出每個(gè)場點(diǎn)對應(yīng)的近地表地殼放大因子,全部結(jié)果見圖3。

圖3 地殼介質(zhì)放大因子

與場地效應(yīng)相關(guān)的高頻衰減部分,本文使用文獻(xiàn)[30]提出的κ0值作為加速度幅值譜在高頻部分的衰減,計(jì)算κ0所用地震動(dòng)數(shù)據(jù)為2018—2021年內(nèi)目標(biāo)區(qū)域及鄰近臺(tái)站記錄到的90組井下水平向地震動(dòng)加速度記錄,所用記錄均來自KiK-net臺(tái)網(wǎng),且經(jīng)過基線校正和濾波處理。按照文獻(xiàn)[31-32]的處理流程,取每條記錄加速度時(shí)程S波到時(shí)到能量持時(shí)80%的部分進(jìn)行傅里葉變換。根據(jù)文獻(xiàn)[33-34]的觀點(diǎn),將傅里葉幅值譜的衰減起始頻率fE控制在5 Hz≤fE<10 Hz,衰減終止頻率fmax控制在10 Hz≤fmax≤20 Hz,且起始頻率與終止頻率的間隔>10 Hz。本文所用起始頻率為5 Hz,終止頻率為20 Hz,對每組傅里葉譜中間的高頻衰減部分使用最小二乘擬合,得到EW分量和NS分量的κ關(guān)于震源距R的關(guān)系式κ=0.047 5+0.000 033 7R和κ=0.048 3+0.000 016 9R(見圖4,圓點(diǎn)表示根據(jù)觀測記錄計(jì)算出的κ值)。當(dāng)R=0 km時(shí)κ記為κ0,表示目標(biāo)研究區(qū)域的局部場地高頻衰減效應(yīng),本文使用2個(gè)水平分量κ0的算數(shù)平均值0.048作為此次模擬的局部場地高頻衰減輸入?yún)?shù)。

圖4 EW和NS向高頻衰減κ值隨震源距R的分布

由于KiK-net具有豐富的井上井下觀測記錄,本文采用物理意義明確的井上井下譜比法對各個(gè)臺(tái)站的局部場地效應(yīng)進(jìn)行估計(jì)。首先,對所有基線校正和濾波后的井上井下水平向地震動(dòng)記錄進(jìn)行傅里葉變換,并使用Konno-Ohmachi 窗對其幅值譜進(jìn)行平滑處理。對平滑后的兩分量傅里葉幅值譜求幾何平均值,再用同一個(gè)臺(tái)站的井上和井下幅值譜均值做譜比,最終得到25個(gè)臺(tái)站的井上井下譜比,將其作為本次模擬的局部場地效應(yīng)參數(shù),全部結(jié)果見圖5。

圖5 井上/井下譜比

本次隨機(jī)有限斷層法模擬所用震源、路徑和場地部分的輸入?yún)?shù)見表2。

表2 隨機(jī)有限斷層法模擬輸入?yún)?shù)

3 模擬結(jié)果及分析

考慮到模擬結(jié)果的穩(wěn)定性和隨機(jī)方法的不確定性,使用有限斷層法對全部井上井下臺(tái)站記錄各進(jìn)行了30次模擬,得到相應(yīng)阻尼比5%的Aps和加速度時(shí)程,并將每組30次模擬Aps的幾何平均值作為目標(biāo)臺(tái)站的Aps模擬結(jié)果。分別選取地表水平向Apg> 60 cm·s-2的井上臺(tái)觀測記錄和其對應(yīng)井下臺(tái)記錄,將其Aps值與模擬值進(jìn)行比較,并在井上部分的對比中加入地震動(dòng)預(yù)測方程(GMPEs)Aps預(yù)測平均值(以下簡稱預(yù)測值),其中GMPEs為文獻(xiàn)[7]提出的日本地區(qū)俯沖帶板內(nèi)水平向地震動(dòng)衰減關(guān)系(以下簡稱Zhao16),Zhao16所用場地見表1,比較結(jié)果見圖6。其中,圖6(a)、(b)分別為井下和井上的比較結(jié)果。整體來看,井上和井下臺(tái)模擬記錄與觀測記錄得到的Aps譜形高度相似,低頻部分模擬稍微高估。由圖6(b)中可以看出,Zhao16也可以較好地描述代表地面運(yùn)動(dòng)的井上臺(tái)記錄頻譜特征,但Zhao16的Aps預(yù)測值在低頻上略低于觀測值和模擬值,模擬得到的Aps值更加接近觀測值,由于以上臺(tái)站處于關(guān)東盆地,盆地效應(yīng)及厚覆蓋層可能對實(shí)際地震動(dòng)的低頻造成一定程度的放大。為驗(yàn)證這種放大與盆地相關(guān),本文在山地區(qū)臺(tái)站中選取海拔>300 m的臺(tái)站(KNGH21、KNGH18、TKYH13、KNGH20)加入比較,見圖6(c)。可以看出,對于山地區(qū)臺(tái)站,Zhao16的Aps值并沒有出現(xiàn)低頻段低估的現(xiàn)象。

對模擬Aps的評價(jià)可用觀測記錄與模擬記錄的Aps值對數(shù)殘差表示,在0.1～10 Hz的頻帶范圍內(nèi),全部井上井下臺(tái)站的模擬記錄與水平向觀測記錄Aps對數(shù)殘差結(jié)果見圖7,圖7(a)為井下結(jié)果,圖7(b)為井上結(jié)果,實(shí)線為對數(shù)殘差平均值,越接近0表示模擬結(jié)果越接近觀測值,陰影區(qū)域?yàn)閷?shù)殘差的一倍標(biāo)準(zhǔn)差范圍?？梢钥闯鲈谒x頻段內(nèi),所有臺(tái)站模擬記錄的Aps對數(shù)殘差平均值均在±0.4以內(nèi),其中0.6～10 Hz頻段基本在±0.2以內(nèi);Aps對數(shù)殘差的一倍標(biāo)準(zhǔn)差主要在±0.5范圍內(nèi)。由此可以看出,模擬記錄能很好地反映本次地震動(dòng)的頻域特征。

圖7 地震動(dòng)模擬記錄與觀測記錄Aps殘差

檢驗(yàn)?zāi)M與觀測記錄加速度時(shí)程的匹配程度,同樣選取地表水平向Apg>60 cm·s-2的臺(tái)站,從模擬30次得到的水平向加速度時(shí)程中挑出Apg與觀測記錄Apg幾何平均值相近的井上記錄進(jìn)行對比,同時(shí)為了檢驗(yàn)70%能量模型的表現(xiàn),本文將廣泛使用的由文獻(xiàn)[35]開發(fā)的地震動(dòng)持時(shí)模型(以下簡稱AB98)加入對比,見圖8?？梢钥闯鱿啾扔贏B98,70%能量模型得到的模擬加速度時(shí)程與觀測記錄在強(qiáng)震動(dòng)段更相似,其持時(shí)也更接近地震動(dòng)觀測記錄,模擬地震動(dòng)記錄在時(shí)域上的表現(xiàn)良好。

圖8 模擬記錄與觀測記錄加速度時(shí)程比較

峰值地面加速度(Apg)是計(jì)算儀器烈度的常用參數(shù),與烈度分布有很強(qiáng)的相關(guān)性。為研究地表Apg特征及其隨距離的衰減特征,本文將25組井上臺(tái)站水平向Apg觀測值、Apg模擬值與Zhao16的Apg預(yù)測值進(jìn)行比較。地表Apg觀測值使用井上水平向Apg的幾何平均值,將模擬30次井上記錄獲得的Apg值進(jìn)行幾何平均作為地表Apg模擬值,Zhao16預(yù)測值所用場地見表1。比較結(jié)果見圖9?？梢钥闯?觀測值和模擬值基本落在Zhao16一倍標(biāo)準(zhǔn)差范圍內(nèi),表明模擬記錄與Zhao16均能很好地描述本次俯沖帶板內(nèi)地震的Apg衰減規(guī)律。此外,Zhao16的Apg預(yù)測值普遍高于地表Apg觀測值,地表模擬值更接近觀測值,大多數(shù)模擬值與觀測值匹配良好。由于Apg主要受地震動(dòng)高頻成分的影響,關(guān)東盆地的厚軟土覆蓋層可能會(huì)抑制地震動(dòng)的高頻成分,從而導(dǎo)致Zhao16在Apg預(yù)測上的高估。相比于Zhao16模型,隨機(jī)有限斷層法可以更靈活地考慮研究區(qū)域的局部場地效應(yīng)。

進(jìn)一步研究本次地震動(dòng)的Apg空間分布,基于井上水平向觀測記錄和模擬記錄,利用連續(xù)曲率樣條插值法方法[36]繪制地表Apg等值線,其中地表Apg模擬值為30次井上模擬值的幾何平均值,結(jié)果見圖10,可以看出觀測與模擬的Apg空間分布相似,模擬記錄得到的Apg在東京地區(qū)略有低估,整體來看具有與觀測記錄較為一致的空間分布特征。

4 結(jié) 論

1)井上井下臺(tái)模擬記錄與觀測記錄的Aps幅值整體對應(yīng)較好,雖然在0.1～0.6 Hz的低頻段內(nèi),模擬Aps普遍存在高估現(xiàn)象,但在中高頻段上(0.6～10 Hz),Aps模擬值與觀測值高度吻合。此外,Zhao16也能很好地描述本次地震動(dòng)的頻譜特征,但盆地區(qū)臺(tái)站的預(yù)測值在低頻出現(xiàn)了一定程度的低估,而山地區(qū)臺(tái)站并未出現(xiàn)上述情況,這可能是由于盆地效應(yīng)及厚覆蓋層對實(shí)際地震動(dòng)造成了低頻的放大。

2)模擬加速度時(shí)程與觀測記錄時(shí)程波形相似,基于70%能量持時(shí)回歸得到的地震動(dòng)持時(shí)模型可以有效改善模擬記錄與觀測記錄在強(qiáng)震動(dòng)段不匹配的情況,模擬加速度記錄的持時(shí)也更接近真實(shí)記錄。

3)地表Apg模擬值與觀測值匹配良好,隨機(jī)有限斷層法和Zhao16均可以很好地描述本次俯沖帶板內(nèi)地震動(dòng)的衰減特性。其中,Zhao16的Apg預(yù)測值相比于模擬值和觀測值普遍出現(xiàn)高估,這可能是目標(biāo)研究區(qū)域軟土層對地震動(dòng)高頻抑制引起的。

4)井上模擬記錄與觀測記錄獲得的Apg等值線相似度較高,通過隨機(jī)有限斷層法獲得的模擬地震動(dòng)在Apg空間分布上與此次觀測到的俯沖帶板內(nèi)地震引起的地震動(dòng)有很好的適配性。