瑪多MW7.3地震同震滑動分布的三維有限元模擬分析

張彩紅 譚 凱 魯小飛 李 琦 李承濤 李志才

1 中國地震局地震大地測量重點(diǎn)實(shí)驗室，武漢市洪山側(cè)路40號，430071 2 國家基礎(chǔ)地理信息中心測繪基準(zhǔn)部，北京市蓮花池西路28號，100830

北京時間2021-05-22 02:04，青海省果洛藏族自治州瑪多縣發(fā)生MW7.3地震[1-2]，震中(34.598°N、98.251°E)位于瑪多縣以南38 km處，處于青藏高原北部巴顏喀拉地塊內(nèi)部的北部邊緣，以及東昆侖南部約70 km的甘德-瑪多斷裂帶上。震源深度約10 km，斷層走向約92°，傾角約67°，朝向為S，破裂斷層走向近EW。美國地質(zhì)調(diào)查局USGS公布此次地震為具有顯著拉張分量的左旋走滑地震事件，而GCMT公布此次地震為接近純左旋走滑事件，二者存在較大差異[3-4]?，敹嗟卣鹗侵袊箨懤^汶川MW8.0地震后發(fā)生的最大地震，震級大、震源淺，造成震中附近房屋、道路、橋梁等基礎(chǔ)設(shè)施不同程度的隆起或坍塌。由于地震所處的巴顏喀拉塊體邊緣斷裂帶的構(gòu)造變形控制著地塊整體向東運(yùn)動，與邊緣地震的發(fā)生息息相關(guān)，因此對邊緣歷史地震的分布特征進(jìn)行研究，有助于深入了解瑪多地震的孕震機(jī)制以及巴顏喀拉塊體邊緣強(qiáng)震的活動特征。

隨著GNSS軟硬件的快速發(fā)展，科研工作者可在震后第一時間解算GNSS數(shù)據(jù)，對地震發(fā)生機(jī)制以及周圍地表影響作出快速響應(yīng)?，敹嗟卣鸢l(fā)生后，本課題組迅速響應(yīng)，組織野外考察隊趕赴現(xiàn)場進(jìn)行GNSS觀測。本文以此次觀測的GNSS站點(diǎn)、周邊CORS站和中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)(以下簡稱陸態(tài)網(wǎng)絡(luò))共40個連續(xù)運(yùn)行站的數(shù)據(jù)為約束，構(gòu)建三維有限元模型，分析和探討瑪多地震同震形變特征，評估周邊地區(qū)的地震危險性。

1 數(shù)據(jù)獲取與解算策略

本文使用的GNSS站均分布在距震中800 km范圍內(nèi)，均勻覆蓋震中區(qū)域(圖1)。距離震中最近的站點(diǎn)為青?，敹嗾?QHMD)，距離震中約38 km；最遠(yuǎn)的為青海茫崖站(QHMY)，距離震中約790 km。40個測站中12個測站的GNSS形變結(jié)果來自文獻(xiàn)[5]；10個測站的形變結(jié)果來自課題組對以往GNSS站進(jìn)行的重復(fù)性觀測，觀測時間為震前2020年doy180～186、doy239～250、doy254～256以及震后2021年doy162～172；18個測站的形變結(jié)果來自陸態(tài)網(wǎng)絡(luò)。

GNSS站和陸態(tài)網(wǎng)絡(luò)共28個測站的數(shù)據(jù)均由Bernese 5.2軟件[6]解算得到。由于Bernese 5.2軟件基線解算模式的精度高于單點(diǎn)定位精度，且定位結(jié)果與初始坐標(biāo)精度有關(guān)，因此為保證單日解精度的可靠性，首先用精密單點(diǎn)定位獲取cm或mm級(與數(shù)據(jù)質(zhì)量有關(guān))的單日解作為初始值，然后采用基線解算模型對中國及周邊IGS核心站進(jìn)行約束，解算得到ITRF2014框架下精度為mm級的單日解。在解算過程中，軌道和鐘差數(shù)據(jù)均采用IGS發(fā)布的最終產(chǎn)品。以CODE中心發(fā)布的全球電離層模型為基礎(chǔ)進(jìn)行高階電離層改正，相同頻率不同碼以及不同頻率之間存在的偏差采用CODE中心每月發(fā)布的碼差分偏差數(shù)據(jù)進(jìn)行校正，數(shù)據(jù)采樣率為30 s，衛(wèi)星截止高度角為5°。天線相位中心偏差采用絕對天線相位中心模型進(jìn)行改正。利用各站點(diǎn)的多天單日解組成三維時間序列，分別擬合地震前后的時間序列，得到高精度同震形變場(圖1)。

圖1 瑪多MW7.3地震同震形變場及青藏高原主要活動斷裂和歷史地震Fig.1 The coseismic deformation field of Madoi MW7.3 earthquake and major active faults and historical earthquakes in the Tibetan plateau

2 三維有限元模型的構(gòu)建

為模擬瑪多地震同震地殼形變，采用Abaqus軟件構(gòu)建彈性三維有限元模型[7]，該軟件能夠根據(jù)彈性斷層錯動、粘彈性松弛和孔隙回彈引起的地表形變進(jìn)行正演模擬，還可以對不同區(qū)域(細(xì)化到每個單元)賦予不同物質(zhì)屬性進(jìn)行模擬。首先在長2 000 km、寬2 000 km、深400 km的三維幾何模型的基礎(chǔ)上，建立長200 km、寬33.3 km、走向N103°E的斷層，斷層中心為34.575 7°N、98.251 7°E；然后對三維幾何模型進(jìn)行網(wǎng)格劃分，在保證解算精度的前提下提高計算效率。采用四面體單元劃分網(wǎng)格，根據(jù)距斷層距離的遠(yuǎn)近設(shè)置單元大小，單元長度由近及遠(yuǎn)逐漸增大，從斷層面及附近單元長度的4 km逐步增大至最遠(yuǎn)處的80 km，從而得到包括400個子斷層、95 571個單元、139 667個節(jié)點(diǎn)數(shù)的精細(xì)三維有限元模型(圖2)。為驗證模型的準(zhǔn)確性，將其與彈性半空間Okada模型[8]進(jìn)行比較。有限元模型的初始條件為零位移，在后續(xù)分析過程中，假設(shè)整個模型均處于彈性狀態(tài)(泊松比為0.25，楊氏模量為75 GPa)，即上表面不加任何應(yīng)力，處于自由狀態(tài)。將遠(yuǎn)場(即模型側(cè)面和底面)設(shè)置為零位移狀態(tài)，依次在每一個子斷層的走滑和傾滑方向上加載單位位移，解算由子斷層錯動引起的地表各節(jié)點(diǎn)位移。通過內(nèi)插方法得到GNSS站位移，從而獲得用于斷層滑動分布反演的格林函數(shù)：

(1)

式中，u為位移；x為走滑或傾滑方向位移；在三維模型中，指數(shù)i、j分別為x、y、z方向上的分量；k為3個分量方向的位移總和；G和λ分別為剪切模量和拉梅常數(shù)；δij為狄拉克函數(shù)。

圖2 瑪多地震三維有限元模型Fig.2 Three-dimensional finite element model for Madoi earthquake

為驗證本文模型和分析方法的準(zhǔn)確性，分別采用Okada模型和有限元模型，在斷層走滑和傾滑方向上加載單位位移，并將其與經(jīng)過斷層中心并垂直于斷層的剖面投影到地表的節(jié)點(diǎn)位移進(jìn)行對比(圖3)。由圖可見，無論是斷層單位走滑還是傾滑，采用有限元模型和Okada方法解算出的地表位移結(jié)果均相同，驗證了本文模型的準(zhǔn)確性。

圖3 斷層單位走滑/傾滑錯動引起的地表位移Fig.3 The surface displacements caused by unit strike/dip dislocation on the fault

在走滑和傾滑方向上依次對子斷層加載單位錯動，得到GPS站點(diǎn)的位移，即通過有限元模型獲取格林函數(shù)。斷層滑動分布公式為：

Gm=d

(2)

式中，G為綜合格林函數(shù)；m和d分別為斷層滑動分布和GNSS觀測到的同震位移。系數(shù)Gij為GNSS站點(diǎn)j上節(jié)點(diǎn)對i加載單位位錯而產(chǎn)生的位移，根據(jù)最小二乘原理可求得斷層滑動分布。為驗證利用最小二乘法反演有限元模型中斷層滑動分布的可靠性，本文對瑪多地震的斷層面賦予初始滑動分布值，如圖4(a)所示，其中白色代表滑動值為1 m，黑色代表滑動值為0 m。首先利用初始滑動分布模型(圖4(a))正演得到分布在震源近場和遠(yuǎn)場上的20個虛擬GNSS站點(diǎn)同震位移場，然后利用最小二乘原理和網(wǎng)格搜索法反復(fù)調(diào)試平滑因子，反演出斷層最優(yōu)滑動分布，如圖4(b)所示。結(jié)果表明，模擬得到的滑動分布與真實(shí)值(給予每個子斷層的單位錯動)在近地表的一致性較好，同震凹凸體主要分布在接近地面的斷層面上。該結(jié)果驗證了采用最小二乘原理和網(wǎng)格搜索法反演斷層滑動分布的準(zhǔn)確性。

圖4 斷層位錯模型和滑動分布Fig.4 The fault dislocation model and slip distribution

3 同震滑動分布分析

首先基于建立的三維有限元模型，利用400個子斷層的單位走滑和傾滑錯動引起的GPS點(diǎn)位移計算得到綜合格林函數(shù)；然后根據(jù)最小二乘原理加入平滑因子進(jìn)行網(wǎng)格搜索，反演斷層同震滑動分布模型；最后根據(jù)滑動分布估計GNSS站點(diǎn)的同震理論位移，將與觀測值最接近的理論值視為最優(yōu)滑動分布。因篇幅有限，本文僅給出最優(yōu)滑動分布以及根據(jù)最優(yōu)滑動分布解算得到的GNSS近場同震形變理論值(圖5)。由圖可見，遠(yuǎn)場GNSS同震形變觀測值不足1 cm，與理論值的差異為mm級。根據(jù)反演得到的滑動分布可以看出，本文構(gòu)建的瑪多地震斷層破裂區(qū)與余震精定位位置[9]一致，反演長度(200 km)大于實(shí)際破裂長度(約160 km)。地震引起的破裂主要分布在野馬灘和黃河鄉(xiāng)附近，破裂長度分別約為3.4 m和2.5 m，野馬灘大橋在地震中坍塌。本文模型與中國地震局地質(zhì)研究所采用InSAR構(gòu)建的滑動分布模型吻合較好，但與美國地質(zhì)調(diào)查局USGS以及北京大學(xué)張勇工作組的滑動分布模型相差較大[3]，可能是因為地質(zhì)所采用的InSAR數(shù)據(jù)以及本文采用的GNSS數(shù)據(jù)均覆蓋瑪多地震震源區(qū)域，可以更好地約束瑪多地震的斷層幾何形狀以及同震滑動分布；而其他2種模型則采用相對稀疏的全球地震波數(shù)據(jù)，且缺乏近場約束。

圖5 瑪多地震滑動分布及GNSS同震位移模擬形變場Fig.5 Coseismic slip distrubution and GNSS coseismic displacements field for Madoi earthquake

GNSS同震形變的水平和高程方向上的理論值與觀測值總體上保持一致，二者的差值稱為殘差。就GNSS水平位移場(圖5(a))而言，距離震源較近的站點(diǎn)QHAJ殘差最大(17 cm)，可能是存在觀測誤差或者幾何模型不夠精細(xì)所致，最小殘差不足1 cm。從同震水平形變的觀測值和理論值可以看出，瑪多地震是以左旋走滑為主的破裂事件，遠(yuǎn)場GNSS站同震水平位移和垂向位移相對較小，僅為0.1～3 cm，說明瑪多地震對周圍地表的影響隨震中距的增大迅速衰減。在GNSS同震形變高程方向上(圖5(b))，距離震源較近的GNSS站點(diǎn)理論值與觀測值相差較小，而距離震源較遠(yuǎn)的QSHE、2838和2819測站，模擬值接近于0，觀測值和模擬值相差較大，殘差最大為5.2 cm。根據(jù)瑪多地震震級和同震形變的水平位移場分布可知，地震對較遠(yuǎn)GNSS站點(diǎn)的位移影響非常有限，理論上垂向同震位移接近于0。2819和2838測站屬于流動觀測站，解算得到的垂向同震位移分別為1.9 cm和2.3 cm，誤差均為0.4 cm，該誤差可能與觀測天數(shù)較少有關(guān)。QSHE測站同震形變結(jié)果來自文獻(xiàn)[5]，垂向同震位移為6.9 cm，解算精度為0.8 cm，解算值和理論值不符，可能與觀測數(shù)據(jù)較少或者觀測環(huán)境不佳[5]有關(guān)。

4 結(jié) 語

根據(jù)GNSS數(shù)據(jù)獲取同震形變場，構(gòu)建瑪多地震三維有限元模型，并用彈性半空間Okada模型驗證其準(zhǔn)確性。通過調(diào)節(jié)滑動因子進(jìn)行網(wǎng)格搜索，反演瑪多地震同震滑動分布。結(jié)果表明，地震引起的破裂主要分布在野馬灘和黃河鄉(xiāng)附近，滑動值分別為3.4 m和2.5 m，與中國地震局地質(zhì)研究所得到的滑動分布結(jié)果相似。地質(zhì)所采用的InSAR數(shù)據(jù)和本文采用的GNSS數(shù)據(jù)均能較好地覆蓋震中附近區(qū)域，二者采用不同的大地測量觀測數(shù)據(jù)反演得到類似的斷層滑動模型，說明本文反演結(jié)果具有可靠性?；瑒臃植冀馑愕玫降腉NSS同震形變理論值與觀測值擬合較好，再次證明本文模型的合理性和準(zhǔn)確性。本文對于瑪多地震后續(xù)研究和該地區(qū)內(nèi)部結(jié)構(gòu)及其他地震的斷層滑動分布構(gòu)造等研究具有重要的理論價值和應(yīng)用意義，也可為相關(guān)研究提供重要的數(shù)據(jù)參考。