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

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

北京時(shí)間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°，朝向?yàn)镾，破裂斷層走向近EW。美國(guó)地質(zhì)調(diào)查局USGS公布此次地震為具有顯著拉張分量的左旋走滑地震事件，而GCMT公布此次地震為接近純左旋走滑事件，二者存在較大差異[3-4]?，敹嗟卣鹗侵袊?guó)大陸繼汶川MW8.0地震后發(fā)生的最大地震，震級(jí)大、震源淺，造成震中附近房屋、道路、橋梁等基礎(chǔ)設(shè)施不同程度的隆起或坍塌。由于地震所處的巴顏喀拉塊體邊緣斷裂帶的構(gòu)造變形控制著地塊整體向東運(yùn)動(dòng)，與邊緣地震的發(fā)生息息相關(guān)，因此對(duì)邊緣歷史地震的分布特征進(jìn)行研究，有助于深入了解瑪多地震的孕震機(jī)制以及巴顏喀拉塊體邊緣強(qiáng)震的活動(dòng)特征。

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

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

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

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

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

(1)

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

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

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

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

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

Gm=d

(2)

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

圖4 斷層位錯(cuò)模型和滑動(dòng)分布Fig.4 The fault dislocation model and slip distribution

3 同震滑動(dòng)分布分析

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

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

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

4 結(jié) 語(yǔ)

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