GP斷層破裂模型在青海門源MS6.9地震強地面運動模擬的適用性研究

羅超,曹曉雨,高陽,徐飛,馮懷平,王昊*

1 石家莊鐵道大學 省部共建交通工程結構力學行為與系統安全國家重點實驗室,石家莊 050043 2 石家莊鐵道大學 土木工程學院,石家莊 050043 3 石家莊鐵道大學 安全工程與應急管理學院,石家莊 050043

0 引言

北京時間2022年1月8日凌晨1時45分青海省海北州門源縣發生面波震級為6.9級的地震,震源深度約為10 km(尹曉菲等,2022).本次地震導致的人員傷亡較少,造成海北州門源等三個縣千余戶受災,9人受傷,超過四千間房屋損壞(王祖東等,2022).本次地震引發了工程結構的顯著破壞,靠近發震斷層的蘭新鐵路硫磺溝大橋和祁連山隧道破壞嚴重(韓帥等,2022).震害表明,近場地震動對工程結構的破壞力巨大,當類似的地震發生在人口密集和重大工程附近時后果不堪設想.

地震發生后,許多學者和研究機構都對門源地震區域斷層和震源參數開展了研究工作.李振洪等(2022)反演了本次地震的震源破裂過程,結果表明本次地震的類型為走滑型地震,子斷層最大滑動量為3.5 m,大約在地下4 km處.潘家偉等(2022)通過詳細野外考查分析得出,本次地震以左旋走滑運動為主,最大走滑位移約為3.7 m,形成了約27 km長的破裂帶,震源深度較淺.尹欣欣等(2022)通過雙差層析成像方法對2009—2022年之間的門源及周邊地區地震事件進行了精確定位處理,定位本次門源地震的震源深度為13.3 km,震中位置為(37.76°N,101.25°E).

在震后重建強地面運動場,是確定震后受災區域和估計未來強震影響場的重要手段.同時隨著一系列地震動輸入方法的提出和廣泛應用(Bielak et al.,2003; Luo et al.,2019,Wang et al.,2021a),對局部場地中密集地震動的獲取提出了要求.已有許多研究人員開展了強地面運動模擬方面的研究,其中隨機有限斷層法被廣泛使用.丁陽等(2018)使用改進的有限斷層法,對九寨溝MS7.0地震進行了模擬,結果顯示,部分模擬結果與強震觀測記錄存在一定差異.周紅等(2021)通過在隨機有限斷層法中考慮破裂過程,模擬了漾濞6.4級地震強地震動場,對地震動場的時空分布規律進行了分析.李春果等(2021)等考慮了瑪多MS7.4地震的多個震源破裂模型,基于隨機有限斷層法開展了強地震動模擬.通過與預測中位值對比,模擬記錄的PGA、PGV結果較為理想,但當與遠場觀測記錄相比時,基于模擬記錄的地震烈度值整體上偏低1～2度.強生銀等(2021)基于隨機有限斷層三維地震動模擬方法,給出了漾濞地震中28個強震動臺站的三分量加速度時程模擬結果.傅磊等(2023)通過隨機有限斷層法模擬了2022年蘆山MS6.1地震,分析了地震中9個距離斷層小于100 km的臺站的模擬記錄.結果顯示,模擬記錄與觀測記錄在加速度時程的S波部分符合良好,在8 s以下周期范圍內傅里葉振幅譜和反應譜形狀和幅值基本一致.綜上所述,隨機有限斷層法雖然計算效率較高,但是只能從宏觀上考慮地形和地質條件的影響,不能完全考慮速度結構、局部場地效應對地震動的影響.

近年來,隨著計算能力的提高,研究人員也開展了有關確定性模擬方法的研究.?，摰?2012)在模擬汶川MS8.0地震的強地震動場時,通過對離散斷層面分級和對各級斷層面破裂上升時間分解的方法,增加了震源的高頻成分.結果表明,PGA模擬結果的誤差不超過50%,傅里葉譜在4 Hz以上頻段比觀測值小.趙宏陽和陳曉非(2017)基于遼寧省海城地區地質資料和已發表的海城MS7.3地震震源機制反演結果,模擬了1975年海城地震的強地震動場.計算結果表明,地震動場的分布特征受到包括震源、地形和地質條件等因素的影響.但是受到模型精度的影響,理論和實際等震線圖存在一定差異.賀春暉等(2017)提出了震源-介質-壩址峽谷場地的多尺度模擬方法,基于譜元法模擬了壩址區域的空間三維地震動場.Wang 等(2021b)在寬頻帶地震動模擬中,通過將震源分成多個層疊加模擬斷層破裂過程實現高頻信號的模擬,并以1994年北嶺地震為例進行了驗證.尹曉菲等(2022)使用曲線網格有限差分法,基于門源區域的地形、速度結構模型以及反演震源破裂過程的結果,模擬了門源地震的強地震動場.受模型精度和計算資源的制約,模擬的地震動場有效頻率不超過0.4 Hz,模擬結果的烈度較強震觀測記錄偏低.上述文獻中,均使用了震后反演得到的斷層破裂模型作為輸入,來模擬設定地震的強地面運動.但是不同地震的斷層破裂模型迥異,能否直接應用在其他有潛在發生地震區域的強地震動模擬中尚存疑問.

Graves和Pitarka提供了另一種思路,他們提出了預測可能發生的斷層破裂過程的GP法(Graves and Pitarka,2004,2010,2015,2016),可用于開展寬頻帶地面運動模擬.Pitarka和Graves(Pitarka et al.,2020)使用GP法和GP&IM混合法建立斷層破裂模型,對日本熊本MW7.0地震進行了寬頻帶地震動模擬;Pitarka等(2022)使用GP法對美國加利福尼亞州RidgecrestMW7.1地震進行了寬頻帶地震動模擬,兩篇文章通過將計算結果與強震觀測記錄對比,展現出了GP法在美國和日本地區優異的適用性.但是GP法是否適用于中國地區的強地震動模擬,目前僅巴振寧等(2023)使用GP模型對漾濞地震進行了模擬,研究工作還有待進一步開展.

研究表明,中國的城市50年內可能發生6級左右直下型地震數量占城市總數的80.4%,且有較多重要城市發生直下型地震的危險性較高(朱剛等,2007).此外,隨著我國重特大工程的建設向深地、深海和中西部地區的轉移,土木工程結構不可避免的需要跨越活動斷層(謝禮立等,2021).著眼城市建設和重大工程中對抗震韌性的需求,雖然已有多種提高結構韌性的減震構件和設計方法(Wang et al.,2023; Zhao et al.,2023; Hu et al.,2023),但是從地震動輸入的角度來看,如何合理地預測和模擬近/跨斷層地震動,對提高土木工程結構的抗震設防能力,實現城市和基礎設施的抗震韌性具有重要的意義和價值.

考慮到建設韌性城市和韌性交通系統中對潛在強震發生區域近/跨斷層地震動場模擬的迫切需求,本文基于Graves和Pitarka提出的GP法(Graves and Pitarka,2004,2010,2015,2016;Pitarka et al.,2020,2022),使用在震前可勘測的活動斷層的走向、傾角、滑動角以及可進行參數分析的震級、震源深度、斷層破裂長度與寬度等震源參數建立斷層破裂模型,然后使用“Seismic Waves,4th order”(簡稱SW4)有限差分計算程序(Petersson and Sj?green,2018)進行寬頻帶地震動模擬,重建門源區域的強地震動場.通過觀察分析模擬結果的衰減規律、模擬結果的峰值(PGA、PGV)和烈度分布與實測地震動記錄對比分析、不同方位臺站的觀測記錄和模擬記錄的波形對比,考察和驗證在中國區域內基于GP法的寬頻帶地震動模擬在震前預測和估計地震動場的能力,為近/跨斷層工程結構的規劃和設計提供參考和科學依據.

1 地震動模擬方法

地震的觀測記錄常被用來進行結構地震反應分析,但現有的地震觀測記錄只對地震動場的一小部分進行了采樣,地面運動記錄通常需要修正或按比例縮放,以此來適應目標場地的地震動特性,采用震后觀測記錄對地震進行預測和評估仍具有一定的局限性.

現如今,譜元法、有限差分法等方法被廣泛應用于寬頻帶地震動模擬中,要實現寬頻帶地震動模擬,受到斷層破裂模型的頻帶范圍、速度結構模型的精細程度和網格尺寸的影響.本文使用SW4有限差分法軟件和具有寬頻特性的GP斷層破裂模型,通過控制網格尺寸不超過最小波長的1/8來實現0～10 Hz的寬頻帶地震動模擬.SW4是一種四階精度的有限差分計算程序(Sj?green and Petersson,2012; Petersson and Sj?green,2017a,b),它具有模擬3D材料屬性、各向異性、使用曲線網格近似地表地形(Petersson and Sj?green,2014)、沿地層深度細化網格(Wang and Petersson,2019)、設置吸收邊界條件最小化來自遠場邊界的人工反射波(Petersson and Sj?green,2015)的功能.

2 門源地震近場區域地震動模擬參數的確定

根據前兩節的介紹,本文研究工作的開展需要確定斷層破裂模型、三維地形地質結構等參數,下面逐一介紹.

2.1 斷層破裂模型

在門源地震發生后,多家相關機構迅速發布了這次地震的震源參數(GFZ,2022;CENC,2022;USGS,2022).本文選用美國地質調查局(USGS,2022)所提供的震源參數,斷層破裂面的走向、傾向和滑動角分別為104°、88°和15°,震中位于(37.828°N,101.290°E),震源深度為13 km.本次地震為走滑型地震,根據Wells統計的5.6

log(RLD)=-2.57+0.62×MW,

(1)

log(RW)=-0.76+0.27×MW.

(2)

經計算求得MW6.6地震斷層破裂的地下長度和寬度分別為33 km和11 km.考慮到斷層的長寬模擬值和最大滑移量應與USGS發布的值盡可能相近,本文模擬的斷層地下破裂長度為28 km,破裂寬度為8 km,最大滑移量為431 cm.綜上所述,根據USGS所提供的震源數據,作為GP法生成斷層破裂模型的參數,震源模型參數匯總見表1.

表1 門源地震震源模型參數

在GP法中,定義整個斷層面的平均上升時間與地震矩的關系如式(3)所示(Graves and Pitarka,2016).

(3)

其中τ為上升時間,αT是與震源機制相關的系數,rc為震源函數上升時間系數,M0為地震矩.在GP法中,默認的rc為1.6,對應的上升時間為0.72 s.

考慮到GP法是由美國西部地區的數據統計得出,而中國地區的地震動加速度反應譜相比于美國西部的地震動加速度反應譜在長周期(T>1.0 s)內幅值小(張齊,2016;張輝,2021),同時USGS提供的門源地震的震源函數上升時間也顯著大于0.72 s,所以在進行地震動模擬時需增加斷層破裂模型的上升時間.因此,為了研究rc對計算結果的影響,增加上升時間為4.2 s(USGS提供門源地震的震源時間函數的上升時間) 時的震源模型和上升時間為2.0 s時的震源模型.即本文將考慮三組上升時間,分別為0.72 s、2 s和4.2 s,所對應的rc為1.6、4.5和9.0.斷層破裂模型是使用南加州地震中心(The Southern California Earthquake Center,SCEC)開發的寬頻帶平臺(Broadband Platform,BBP)生成(Maechling et al.,2015),詳細的內容可以查看網站https:∥github.com/SCECcode/bbp/wiki.本文基于BBP的現有代碼,在僅改變參數rc的情況下生成斷層破裂模型,斷層的具體參數設置如表1所示,同時設置震源位于斷層破裂面的幾何中心,子斷層的大小為0.1 km×0.1 km,BBP生成的斷層破裂模型給出了每個子斷層的錯動速率時程,進而可計算斷層平面上滑移分布、起始破裂時間與滑移率等,斷層破裂模型文件詳細的內容和格式可以查看網站Standard Rupture Format-SCECpedia.在其他震源模型參數不變的情況下,改變rc時斷層滑移分布、起始破裂時間均不變,斷層滑移分布、起始破裂時間與滑移率如圖1所示.

圖1 斷層破裂隨機滑移分布圖(a); 起始破裂時間(b); 圖(c)—(e) rc值依次為1.6、4.5、9.0時斷層的峰值滑移率

2.2 地震動模擬研究區域

本文分為區域1、區域2、區域3三個區域進行寬頻帶地震動模擬,區域范圍和尺寸如表2和圖2所示.受計算能力所限,區域1中難以開展(0～10 Hz)地震動模擬,主要用于模擬速度場和估計區域尺度的烈度分布,地震動模擬的頻帶范圍為(0～2 Hz).區域2用于模擬PGA和PGV隨震中距的衰減分布規律,地震動模擬的頻帶范圍為(0～10 Hz).區域3用于模擬近斷層區域寬頻帶地震動場,對比分析近場區域臺站觀測記錄和模擬記錄的波形和峰值,地震動模擬的頻帶范圍為(0～10 Hz).本文使用SW4軟件進行地震動模擬時,需要用到輸入腳本文件、地形模型、速度結構模型和斷層破裂模型.輸入腳本文件內的代碼在計算時會讀取地形模型、速度結構模型和斷層破裂模型,同時會定義計算區域的原點位置、計算區域的總尺寸、網格尺寸、模擬的總時間、模擬網格點的輸出等關鍵運行參數,如表2所示.

圖2 門源地震震中位置、研究區域地形分布圖

表2 SW4關鍵運行參數

2.3 門源地區的地形模型與速度結構模型

地殼介質的波速、密度和品質因子等特性在空間中存在顯著的差異,建立和使用合理的地殼速度結構模型是開展地震動模擬的基礎.本次門源地震發生在青海省門源縣,在本文研究區域內(36.6°N—39.0°N,100.0°E—103.0°E),使用美國地質調查局發布的具有30 m精度的全球地形數據STRM1(Farr et al.,2007)建立研究區域內的地形模型.使用模型USTClitho2.0(Han et al.,2022)來建立速度結構模型,由于模型USTClitho2.0中沒有地下巖層的密度,使用全球地殼模型CRUST1.0(Laske et al.,2013)中的地下巖層密度進行補充.P波與S波的深度切片分布圖如圖3所示,P波與S波豎向剖面分布圖如圖4、圖5所示.地震波的非彈性衰減采用粗粒度方法,品質因子根據Pitarka和Graves給出的經驗關系公式計算(Pitarka and Graves,2016):

圖3 P波與S波隨深度變化切片分布

圖4 P波在垂直方向上變化切片分布

圖5 S波在垂直方向上變化切片分布

QS=vS×50,

(4)

QP=2×QS,

(5)

其中,QS為P波品質因子,QP為S波品質因子,vS為S波波速.

3 計算結果與討論

根據前文介紹的模擬方法、震源參數、地形結構模型與地殼結構模型,按照上升時間0.72 s、2.0 s、4.2 s,將計算工況分為三組,對應的rc分別為1.6、4.5、9.0,對目標區域的地震動場開展強地震動模擬.為了驗證本文計算模型的合理性,選取區域2內觀測點的模擬記錄的東西向、南北向的PGA和PGV,與俞言祥等(2013)編制建立的我國青藏地震區基巖場地水平地震動峰值預測方程進行對比,地震動峰值預測方程的計算公式為

lgY=A+BM+Clg(R+DeEM)±σ,

(6)

其中Y代表地震動加速度和速度峰值,A、B、C、D、E代表回歸分析得到的系數,M、R和σ分別代表面波震級、震中距和標準差.預測方程與PGA模擬記錄對比時,預測方程的回歸系數A、B、C、D、E分別取值為2.457、0.388、-1.854、0.612、0.457,M和σ分別取6.9和0.236;預測方程與PGV模擬記錄對比時,預測方程的回歸系數A、B、C、D、E分別取值為0.443、0.474、-1.696、0.612、0.457,M和σ分別取6.9和0.271.圖6繪制了基于不同rc的模擬記錄峰值地面加速度(PGA)和峰值地面速度(PGV)隨震中距的衰減規律,其中黑色實線代表預測中位值,黑色虛線表示±1倍標準差.

圖6 基于不同rc的模擬記錄 PGA 和 PGV 隨距離的衰減分布

由圖6可知:三組不同rc的模擬記錄和預測中位值的PGA與PGV隨距離的變化總體衰減速率基本一致.當rc為1.6時,模擬結果PGA和PGV較預測中位值偏大.當rc為4.5和9.0時,模擬結果PGA和PGV總體上與預測中位值較符合.隨著rc的增加,在相同震中距的情況下,三組模擬結果的PGA和PGV均逐漸降低.這是因為周期為1～2 s的地面運動受rc的影響顯著,即rc的值越大,上升時間越大,峰值滑移率越小,地震動速度和加速度的幅值越小.綜上所述,當rc為4.5與9.0時,模擬結果與預測中位值符合較好.

表3給出了距離震中最近的兩個臺站C0028、C0029的觀測記錄、模擬記錄在10 Hz低通濾波后的PGA和PGV.從表3可以看出,當rc為1.6時,臺站C0028模擬記錄的PGA和PGV在南北向較觀測記錄偏大,東西向和豎向與觀測記錄較符合;臺站C0029模擬記錄的PGA和PGV在三個方向與觀測記錄較符合.當rc為4.5和9.0時,臺站C0028與臺站C0029模擬記錄PGA與PGV的三個方向與觀測記錄均較符合.綜上所述,當rc為4.5與9.0時,模擬結果與觀測記錄符合較好.

表3 臺站C0028、C0029觀測與模擬記錄的PGA、PGV

根據《中國地震烈度表》(國家市場監督管理總局和國家標準化管理委員會,2020),計算研究區域地震烈度.結合模擬所得加速度時程與速度時程,計算研究區域地震烈度.本次門源震后國家預警工程自動產出結果(王士成等,2017)見圖7a,rc為1.6、4.5和9.0對應的研究區域的地震烈度空間分布分別由圖7b、7c和7d表示.圖8為不同rc參數下PGA和PGV地表分布圖.

圖7 國家預警工程自動產出烈度圖與地震動模擬產出烈度圖

圖8 不同rc參數下PGA和PGV相應模擬量的地表分布圖

從圖7a中可知,斷層附近區域極震區的最大烈度可達9度,本文模擬的三幅圖與國家預警工程的最大烈度Ⅸ相一致.對比圖7a和圖7d可知,圖7d高烈度區(Ⅵ、Ⅶ、Ⅷ、Ⅸ度區)分布的經緯度范圍與圖7a較符合,具體的經緯度范圍如表4所示;圖7a中高烈度區在平行和垂直于斷層方向的烈度較其他方向衰減慢,與圖7d中高烈度區呈垂直斷層和平行斷層分布特征相似;從表5所示的各烈度區面積可以看出圖7d中高烈度區與圖7a符合度較好.

表4 自動產出結果與地震動模擬結果Ⅵ、Ⅶ、Ⅷ與Ⅸ度區范圍

表5 自動產出結果與地震動模擬結果Ⅶ、Ⅷ與Ⅸ度區面積

從表5和圖7中可以看出,隨著rc的增大,沿著斷層走向的高烈度區面積逐漸減小,垂直于斷層的高烈度區面積基本不變.當rc為1.6時,模擬結果的Ⅶ度區面積大于自動產出結果,Ⅷ、Ⅸ度區面積遠大于自動產出結果;當rc為4.5時,模擬結果Ⅶ度區面積與自動產出結果吻合較好,Ⅷ度區面積遠大于自動產出結果;當rc為9.0時,與自動產出結果相比,模擬結果Ⅶ、Ⅷ度區面積略偏大.綜上所述,隨著rc的增加,模擬的各烈度區域面積整體呈現減小趨勢,這是由于rc對1～2 s周期地面運動的影響顯著,即rc的值越小,速度幅值越大,烈度區的面積越大.

當rc為9.0時,模擬的各烈度區面積與國家預警工程較符合,但是Ⅶ、Ⅷ度區面積較自動產出記錄仍然略有偏大.一般來說,臺站多布置在地形較為平坦的區域,而地震動幅值在山地區域會有顯著放大,從而導致了模擬結果的Ⅶ、Ⅷ度區面積偏大.圖9繪制了臺站C0029附近區域的地形圖和南北方向速度分布云圖.可以看出,圖中北部和東北部山地區域出現了顯著的地震動放大效應.考慮到模擬區域山區分布廣泛,本文模擬的Ⅶ、Ⅷ度區面積較自動產出列度結果偏大.

圖9 臺站C0029周圍地形圖(a)與臺站C0029附近區域南北方向PGV(b)

當rc=9.0時,臺站C0028和C0029在10 Hz低通濾波下的觀測記錄波形和模擬記錄波形如圖10所示.其中圖10a為臺站C0028觀測記錄和模擬記錄的加速度波形對比圖;圖10b為臺站C0028觀測記錄和模擬記錄的速度波形對比圖.圖10c為臺站C0029觀測記錄和模擬記錄的加速度波形對比圖;圖10d為臺站C0029觀測記錄和模擬記錄的速度波形對比圖.從圖10中可以看出,除臺站C0029在南北向模擬記錄的速度峰值較小,兩臺站觀測記錄和模擬記錄的加速度和速度在峰值和波形上均較為符合.

圖10 臺站C0028和C0029在10 Hz低通濾波下的觀測記錄波形和模擬記錄波形

綜合模擬記錄峰值與青藏地區的地震動峰值預測方程、強震觀測記錄和震后國家預警工程自動產出烈度區分布的對比結果,同時對比分析臺站C0028和C0029在10 Hz低通濾波下的觀測記錄波形和模擬記錄波形.可以看出,當rc為9.0時,地震動模擬結果與觀測記錄的符合度較好.

4 結論

本文采用Graves和Pitarka斷層破裂模型(GP法)和“Seismic Waves,4th order” (SW4)有限差分計算程序,對2022年門源MS6.9地震開展了強地面地運動模擬,開展了以下工作:(1)對比了模擬記錄峰值與青藏地區的地震動峰值預測中位值隨震中距的衰減規律;(2)將C0028和C0029兩個臺站的模擬記錄與觀測記錄的PGA和PGV進行對比分析;(3)根據模擬記錄給出烈度分布圖,并與國家預警工程自動產出烈度圖進行對比分析;(4)對比分析C0028和C0029兩個臺站的模擬記錄與觀測記錄的波形.可以得到如下結論:

(1) 對于三種不同上升時間系數(rc)的震源模型,rc為4.5和9.0時模擬記錄的峰值(PGA、PGV)與青藏地區的地震動預測方程的預測中位值符合程度較好.

(2) 在近斷層區域寬頻帶地震動場的計算結果中,對于臺站C0028和C0029,rc為4.5和9.0時模擬記錄與強震觀測記錄的PGA和PGV基本吻合.

(3) 本文模擬的各烈度區面積隨著rc的增大而減小.當rc為9.0時烈度區走向、范圍和烈度區面積與國家預警工程自動產出結果基本吻合.

(4) 當rc為9.0時,臺站C0028和C0029的模擬記錄和觀測記錄的加速度和速度在峰值和波形上均較為符合.

(5) 綜上所述,本次強地震動模擬所得結果較好的重現了門源地震,在合理選取rc參數的情況下,GP法較好的重現了門源地震動場,但是影響強地震動模擬的因素還很多,未來還需開展破裂速度、震源時間函數等參數的影響和基于中國大陸強震觀測記錄的震源破裂模型研究.

致謝感謝中國地震局工程力學研究所“國家強震動臺網中心”為本研究提供數據支持.感謝審稿專家和編輯對文稿提出了寶貴的修改意見和建議.同時感謝開源繪圖軟件GMT(Wessel et al.,2019)為本文提供了方便高效的繪圖平臺.