2022年1月8日青海門源MS6.9地震強地面運動模擬及烈度分布估計＊

尹曉菲 王 芃， 張 偉 邵志剛曹文忠 熊仁偉

1） 中國北京 100036 中國地震局地震預(yù)測研究所

2） 中國廣東深圳 518055 南方科技大學地球與空間科學系

引言

2022年1月8日01時45分27秒青海省海北藏族自治州門源回族自治縣發(fā)生MS6.9地震，震中位置為（37.77°N，101.26°E），震源深度為10 km.本次地震造成了17 069人受災(zāi)，由于震中距離人口稠密地區(qū)較遠，尚無人員傷亡的相關(guān)報道（青海日報，2022），蘭新高鐵浩門至山丹軍馬場區(qū)間隧道群由于此次地震而發(fā)生局部塌方（央視網(wǎng)，2022）.門源MS6.9地震發(fā)生在青藏高原東北緣冷龍嶺斷裂、托萊山斷裂和肅南—祁連斷裂的階區(qū)部位.現(xiàn)場工作隊在冷龍嶺斷裂帶西段探測到長約22 km的地表破裂帶（青海省地震局，2022）.托萊山斷裂和冷龍嶺斷裂現(xiàn)為左旋兼擠壓斷裂（李強等，2013），歷史上無M7以上地震記載（姜文亮，2018），最近一次強震為2016年1月21日門源MS6.4地震，該地震的發(fā)生反映了青藏高原地塊向NE向不斷推擠生長的過程（胡朝忠等，2016）.門源MS6.9地震發(fā)生后，許英才等（2022）對門源地震早期序列（2022年1月8日至12日）進行了重定位和震源機制研究，其結(jié)果表明目前門源地區(qū)還存在一定的應(yīng)力積累且應(yīng)力尚未得到充分釋放，該地區(qū)仍有發(fā)生強震的危險.由于地震災(zāi)害主要是地震在地表產(chǎn)生的強地面運動造成的，因此為了獲得地震波場傳播過程及其引起的地表響應(yīng)，對門源MS6.9地震開展強地面運動初步模擬及烈度估計具有重要意義.

地震波場的正演模擬需要考慮復(fù)雜地表.曲線網(wǎng)格有限差分方法（Zhang，Chen，2006；Zhanget al，2012）適用于含起伏地形的地震波場模擬，該方法中地表的形狀用任意曲線網(wǎng)格近似，并采用牽引力鏡像法處理自由表面條件.曲線網(wǎng)格有限差分的強地面運動模擬方法在地震后的災(zāi)害評估中得到了非常廣泛的應(yīng)用.Zhang等（2008）結(jié)合三維介質(zhì)模型、震源破裂模型和地表地形數(shù)據(jù)模擬了2008年汶川MS8.0地震的強地面運動，研究表明斷層破裂方式和盆地構(gòu)造控制了地表峰值速度（peak ground velocity，縮寫為PGV）的分布，地表起伏劇烈地區(qū)往往對應(yīng)較大的PGV數(shù)值，應(yīng)關(guān)注其震害問題.張振國等（2014a，b）對2014年2月12日新疆于田MS7.3地震和2014年8月3日云南魯?shù)镸S6.5地震引起的強地面運動作了初步模擬和烈度預(yù)測，研究顯示地震動在山峰、山脊處具有較大幅值，該結(jié)果可指導(dǎo)震后的災(zāi)區(qū)重建工作.趙宏陽和陳曉非（2017）利用1975年2月4日遼寧海城MS7.3地震的地震地質(zhì)資料，模擬計算了海城地震的波場傳播過程，分析了強地面運動的方向性效應(yīng)、盆地效應(yīng)和近斷層效應(yīng)，得出理論烈度分布同震后調(diào)查烈度分布基本一致，驗證了研究所用的震源模型和速度結(jié)構(gòu)的合理性.

本文擬根據(jù)張勇①張勇2022年1月10日與作者的個人交流提供的門源地震震源破裂過程的初步結(jié)果，利用曲線網(wǎng)格有限差分方法模擬門源MS6.9地震的強地面運動，再結(jié)合地表峰值速度和烈度間的關(guān)系，計算地震烈度，并在此基礎(chǔ)上評估地震災(zāi)害分布特征，以期為門源地區(qū)的防震減災(zāi)提供科學依據(jù).

1 研究區(qū)概況及計算方法和模型

研究區(qū)范圍如圖1所示，可以看出，該區(qū)域的地貌特征變化大，地勢西南高東北低，高程介于1.3—5.0 km之間.門源MS6.9地震發(fā)生后至2022年1月20日，該地區(qū)已發(fā)生MS≥4.0余震23次，其中MS≥5.0余震2次.

圖1 研究區(qū)及周邊的地質(zhì)概況和2022年1月8—20日震源區(qū)MS≥3.0 地震分布Fig. 1 Geological overview of the study area and its surrounding regions，spatial distribution of MS≥3.0 earthquakes in source zone during the period of 8 to 20 January，2022

采用曲線網(wǎng)格有限差分方法（Zhang，Chen，2006；Zhanget al，2012）對門源MS6.9地震的強地面運動進行模擬.計算中需要設(shè)置描述地形起伏的網(wǎng)格模型、反映地下物質(zhì)屬性的介質(zhì)模型，以及表示地震破裂過程的震源模型.

地形選取GTOPO30地形數(shù)據(jù)，其水平分辨率大約為1 km.整個計算區(qū)域尺度為350 km×220 km，深度為60 km.將研究區(qū)域離散成700×440×120個網(wǎng)格，垂直方向采用等間距排列，單位網(wǎng)格為邊長500 m的立方體.

考慮到面波對縱波速度的靈敏度和起伏地形的影響，Han等（2022）提出一種改進的體波和面波數(shù)據(jù)聯(lián)合反演方法獲得中國大陸地殼和上地幔水平分辨率為0.5°的速度結(jié)構(gòu)（USTClitho2.0模型）.該模型對于青藏高原地塊內(nèi)體波射線覆蓋相對稀疏的區(qū)域，加入面波頻散數(shù)據(jù)可更好地提高縱波速度和橫波速度的準確度，因此本文以USTClitho2.0模型為基礎(chǔ)建立研究區(qū)域的網(wǎng)格化橫波速度和縱波速度模型.每個網(wǎng)格的東向、北向和垂直方向的長度分別為43.75 ，55 和5 km，同時依據(jù)Brocher （2005）提出的密度和縱波速度的經(jīng)驗關(guān)系確定各網(wǎng)格的密度.模擬采用的橫波速度結(jié)構(gòu)如圖2所示，可以看出，2—12 km深度存在高速層，沿北東方向地殼厚度逐漸減薄.雖然青海湖處于研究區(qū)域內(nèi)（圖1），但是由于其位于模擬計算區(qū)域的邊緣，且水深較淺，因此在建立介質(zhì)模型時忽略其影響.

圖2 橫波速度結(jié)構(gòu)Fig. 2 S-wave velocity structure

張勇①張勇2022年1月10日與作者的個人交流反演了該地震的震源破裂過程，結(jié)果顯示斷層滑動出露地表，走向為103°，傾角為88°，滑動角為-6°，共有31×11＝341個子斷層，子斷層的空間分辨率為2 km×2 km，采樣點為80個，每間隔0.25 s給出一個滑動速率，破裂過程持續(xù)時間約20 s，最大滑動量為1.8 m.基于這一結(jié)果，本文設(shè)置了相應(yīng)的震源模型，如圖3所示.

圖3 本文所用的門源MS6.9地震震源模型Fig. 3 Source model of Menyuan MS6.9 earthquake used in this study

2 強地面運動模擬結(jié)果

強地面運動模擬計算的時間步長為0.01 s，總步數(shù)為1萬步，模擬時長共100 s，使用640個計算核心進行計算.門源地震x分量的模擬速度場快照如圖4所示，可見：地震的主要能量由位于震中附近下方的位錯產(chǎn)生，且以水平走向錯動為主，因此在斷層垂向上產(chǎn)生了能量較強的S波；大約第2.36 s時，地震波到達地表，在初 始破裂時刻，速度的最大值集中在斷層破裂的前鋒上；隨著遠離發(fā)震斷層，地震波場能量逐漸減小.

圖4 門源MS6.9地震x分量粒子速度波場快照Fig. 4 Wavefield snapshots of the x component particle velocity of Menyuan MS6.9 earthquake

通過強地面運動模擬得到各網(wǎng)格點x，y，z三個方向不同時刻的速度，然后對三個方向的速度分量求矢量和獲得各網(wǎng)格點的運動速度隨時間的變化，根據(jù)各網(wǎng)格點的速度最大值確定研究區(qū)域的地表峰值速度（PGV）分布，并結(jié)合國家市場監(jiān)督管理總局和國家標準化管理委員會（2021）給出的PGV與烈度之間的關(guān)系得到對應(yīng)的地震烈度，如圖5所示，結(jié)果表明：沿平行斷層走向方向的地震動衰減明顯小于垂直斷層走向方向；地震的最大烈度為Ⅷ度（PGV＝37.63 cm/s），位于震源破裂起始點附近區(qū)域；Ⅷ度區(qū)主要涉及門源縣、祁連縣和肅南縣的部分區(qū)域；Ⅶ度區(qū)主要涉及大通縣、永昌縣、民樂縣等部分區(qū)域.中國地震局（2022）野外調(diào)查的烈度分布顯示等震線長軸呈WNW走向，同理論模擬顯示的高烈度區(qū)主要沿WNW方向延伸的結(jié)果一致.理論烈度分布與中國地震局工程力學研究所強震動觀測組（2022）給出的儀器烈度分布（圖6）較為接近，即高烈度區(qū)主要沿斷層走向展布，且斷層南側(cè)的地震烈度大于北側(cè).受模型分辨率和計算成本的限制，本文模擬的地震波場的最高有效頻率為0.4 Hz（2.5 s），缺少高頻成分，因此獲得的烈度較強震儀和烈度儀觀測的結(jié)果偏低，但相較于使用低分辨率介質(zhì)模型的模擬結(jié)果（徐劍俠等，2015），地震動模擬的最高頻率得到了一定的提升.

圖5 門源MS6.9地震烈度分布圖Fig. 5 The simulated intensity of Menyuan MS6.9 earthquake

圖6 門源MS6.9地震儀器觀測的烈度分布（引自中國地震局工程力學研究所強震動觀測組，2022）Fig. 6 Instrumental seismic intensity distribution of Menyuan MS6.9 earthquake （after Strong Motion Observation Group，Institute of Engineering Mechanics，China Earthquake Administration，2022）

地震災(zāi)害區(qū)域集中在發(fā)震斷層附近，并向WNW方向和ESE方向延伸，這主要由兩個原因造成：一方面，該地震為高傾角的走滑型地震，地震波能量主要沿走向方向傳播，表現(xiàn)出強地面運動的方向性效應(yīng)；另一方面，起伏地表對地震波傳播具有重要影響（Zhanget al，2008），發(fā)震斷層的WNW方向和ESE方向分布有托萊山、大通山、達坂山和冷龍嶺等山脈（圖1），屬于山地地貌，因此該區(qū)域的地震動高值可能與山脊地區(qū)地震波多次反射有關(guān).

3 討論與結(jié)論

強地面運動模擬結(jié)果由地形數(shù)據(jù)、介質(zhì)模型和震源模型共同決定，因此數(shù)據(jù)選擇對結(jié)果的可靠性有較大影響，為使結(jié)果更為可信，本文參考美國加州綜合地震破裂預(yù)測模型（Fieldet al，2014）中提出的“使用可獲得最優(yōu)解”原則對介質(zhì)模型進行選擇.常用的介質(zhì)模型包括CRUST1.0模型（Laskeet al，2012）和CRUST2.0模型（Bassinet al，2000），二者均基于水平層狀介質(zhì)的假設(shè)，水平分辨率分別為1.0°和2.0°，而Han等（2022）提出的USTClitho2.0模型的水平分辨率更高，因此研究中選擇USTClitho2.0模型作為介質(zhì)模型.本文基于門源地震震源區(qū)附近的地形數(shù)據(jù)、介質(zhì)模型和震源模型，使用曲線網(wǎng)格有限差分方法計算了該地震的近場地震波傳播過程，得到了理論的速度場快照和地震烈度分布.研究結(jié)果表明：總體上垂直斷層走向方向的地震動衰減大于平行斷層走向方向的地震動衰減，震中最大烈度為Ⅷ度；模擬得到的理論烈度同野外調(diào)查的地震烈度分布基本一致.受強地面運動方向性效應(yīng)和起伏地表的影響，地震災(zāi)害主要沿斷層的WNW方向和ESE方向分布.由于地震發(fā)生區(qū)域主要以山地地貌為主，該地震的發(fā)生造成了局部邊坡崩塌、滾石和凍土開裂以及蘭新高鐵大橋橋面受損等次生災(zāi)害（頡滿斌，2022）.考慮到研究區(qū)域仍有發(fā)生強震的危險（許英才等，2022），今后的防震減災(zāi)工作中有必要加強地表起伏劇烈區(qū)域特別是山脊地貌的震害防御工作.

北京大學張勇教授為本文提供了震源破裂過程的初步結(jié)果，中國科學技術(shù)大學張海江教授為本文提供了三維速度結(jié)構(gòu)，中國地震局工程力學研究所馬強研究員為本文提供了儀器觀測的地震烈度分布圖，中國地震局地震預(yù)測研究所徐岳仁研究員與作者就青藏高原東北緣的地質(zhì)概況進行了討論，審稿專家為本文提出了寶貴意見，作者在此一并表示感謝.