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

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

1 中國地震局地震大地測量重點實驗室，武漢市洪山側路40號，430071 2 國家基礎地理信息中心測繪基準部，北京市蓮花池西路28號，100830

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

隨著GNSS軟硬件的快速發展，科研工作者可在震后第一時間解算GNSS數據，對地震發生機制以及周圍地表影響作出快速響應。瑪多地震發生后，本課題組迅速響應，組織野外考察隊趕赴現場進行GNSS觀測。本文以此次觀測的GNSS站點、周邊CORS站和中國大陸構造環境監測網絡(以下簡稱陸態網絡)共40個連續運行站的數據為約束，構建三維有限元模型，分析和探討瑪多地震同震形變特征，評估周邊地區的地震危險性。

1 數據獲取與解算策略

本文使用的GNSS站均分布在距震中800 km范圍內，均勻覆蓋震中區域(圖1)。距離震中最近的站點為青海瑪多站(QHMD)，距離震中約38 km；最遠的為青海茫崖站(QHMY)，距離震中約790 km。40個測站中12個測站的GNSS形變結果來自文獻[5]；10個測站的形變結果來自課題組對以往GNSS站進行的重復性觀測，觀測時間為震前2020年doy180～186、doy239～250、doy254～256以及震后2021年doy162～172；18個測站的形變結果來自陸態網絡。

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

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

(1)

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

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

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

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

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

Gm=d

(2)

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

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

3 同震滑動分布分析

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

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

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

4 結 語

根據GNSS數據獲取同震形變場，構建瑪多地震三維有限元模型，并用彈性半空間Okada模型驗證其準確性。通過調節滑動因子進行網格搜索，反演瑪多地震同震滑動分布。結果表明，地震引起的破裂主要分布在野馬灘和黃河鄉附近，滑動值分別為3.4 m和2.5 m，與中國地震局地質研究所得到的滑動分布結果相似。地質所采用的InSAR數據和本文采用的GNSS數據均能較好地覆蓋震中附近區域，二者采用不同的大地測量觀測數據反演得到類似的斷層滑動模型，說明本文反演結果具有可靠性。滑動分布解算得到的GNSS同震形變理論值與觀測值擬合較好，再次證明本文模型的合理性和準確性。本文對于瑪多地震后續研究和該地區內部結構及其他地震的斷層滑動分布構造等研究具有重要的理論價值和應用意義，也可為相關研究提供重要的數據參考。