趙磊， 許文斌， 方楠， 劉計洪， 馮光財

中南大學地球科學與信息物理學院， 長沙 410000

0 引言

2021年5月22日，在青海瑪多發(fā)生MW7.3地震，據(jù)中國地震臺網(wǎng)觀測，震源深度為17 km，震中(34.59°N，98.34°E)位于巴顏喀拉塊體內(nèi)部的江錯斷裂附近(圖1)，距離塊體北邊界的昆侖斷層南側(cè)約70 km(王未來等，2021；潘家偉等，2021).歐亞板塊和印度板塊的持續(xù)碰撞，在青藏高原形成了多條大型斷裂帶，這些斷裂帶把整個高原劃分為多個次級塊體.巴顏喀拉地塊為青藏高原中部的次級地塊，是一條狹長的三角構(gòu)造單元.巴顏喀拉地塊南北兩側(cè)分別受鮮水河斷裂和東昆侖斷裂控制.塊體向東擠出的過程中受到四川盆地的阻擋，形成了龍門山擠壓造山帶，即青藏高原東邊界.在過去20多年時間里，多次中強地震圍繞巴顏喀拉塊體發(fā)生(Klinger et al., 2005; 王衛(wèi)民等，2008；Xu et al., 2010; Huang et al., 2013; Liu et al., 2014; 季靈運等，2017；賀鵬超等，2018；Zhao et al., 2018)，但鮮有對于該塊體內(nèi)部活動斷層的研究.瑪多地震為我們研究塊體內(nèi)部斷裂的活動方式和構(gòu)造應力特征提供了很好的機會，而精確的發(fā)震斷層模型對于分析周緣斷層的活動方式和構(gòu)造應力起到關(guān)鍵作用.

圖1 區(qū)域構(gòu)造背景圖(a) 紫色圓點表示余震分布(Wang et al., 2021)，白色五角星代表震中位置，紅色沙灘球表示United States Geological Survey(USGS)瑪多地震震源機制解，綠色實線表示斷層破裂地表跡線，黑色實線表示活動斷層，藍色箭頭表示震間GPS速率(Wang and Shen, 2020)； (b) 紅色虛線表示巴顏喀拉塊體，黃色和綠色方框分別表示降軌和升軌數(shù)據(jù)范圍，紫色沙灘球表示歷史地震震源機制解(USGS).Fig.1 Regional tectonic setting(a) The purple dots denote the relocation aftershocks (from Wang et al., 2021), the white star represents the location of the epicenter, red beach balls show the focal mechanisms of the Madoi earthquake from United States Geological Survey (USGS), the green line indicates the rupture surface trace, the black lines denote the active faults, the blue arrows indicate the interseismic GPS velocity (Wang and Shen, 2020). (b) The red dotted line represents the Bayan Har block, the blue and green boxes represent the descending and ascending frame respectively, purple beach balls show the focal mechanism of historical earthquakes from the USGS.

已有研究成果表明，瑪多地震的余震分布和InSAR同震形變場呈現(xiàn)一個復雜的發(fā)震斷層幾何結(jié)構(gòu)，東西向延伸約150 km，發(fā)震斷層在東端分支，在西端改變斷層走向(Zhao et al., 2021; Wang et al., 2021; He et al., 2021).基于大地測量數(shù)據(jù)解算的同震三維形變場顯示垂直形變在斷層兩側(cè)正負交替，表明了斷層傾向在空間分布上的復雜性(Zhao et al., 2021; 姜衛(wèi)平等，2022).同震庫侖應力擾動表明此次地震對周緣斷層具有明顯的應力加載作用(祝愛玉等，2021).靜態(tài)庫侖應力結(jié)果顯示余震分布與應力加載區(qū)域有很好的空間相關(guān)性(王迪晉等，2022).雖然已有的研究提出了多種斷層模型(He et al., 2021; Jin and Fialko, 2021; Wang et al., 2021; Zhao et al., 2021; 王閱兵等，2022；余鵬飛等，2022He et al., 2022)，但大多基于余震和線性反演確定斷層幾何結(jié)構(gòu).對于余滑和余震的空間相關(guān)性及震后余滑的時空變化特征的研究較少.為了得到精確的斷層幾何結(jié)構(gòu)，本文聯(lián)合貝葉斯非線性反演方法和網(wǎng)格搜索反演方法估計斷層幾何和滑動分布模型.然后基于觀測數(shù)據(jù)和精確斷層模型反演同震和早期震后滑動分布并探討兩者之間的空間相關(guān)性，進而研究斷層的摩擦屬性.

1 數(shù)據(jù)處理和建模方法

1.1 同震和震后數(shù)據(jù)處理

本文的InSAR同震形變場和偏移量來自Liu等(2022)，而震后4.5個月內(nèi)的形變場利用GAMMA軟件處理得到.采用震后降軌21景和升軌17景影像，通過雙軌差分干涉得到差分干涉圖.在數(shù)據(jù)處理過程中，距離向和方位向分別采用10和2進行多視處理提高信噪比，采用90 m分辨率的數(shù)字高程模型(SRTM DEM)用于校正地形相位.采用自適應濾波和最小費用流法(Goldstein and Werner, 1998)得到解纏干涉圖.之后掩膜相干性低于0.7的像素點和剔除垂直基線大于30 m的干涉對來提高觀測區(qū)域相干性.最終分別獲得77景降軌干涉圖和68景升軌干涉圖，使用開源軟件LICSBAS(Morishita, 2021)進行時序震后形變場解算.首先，將閉環(huán)相位值設(shè)為1.5 rad并且計算每三個干涉圖的閉環(huán)相位，剔除具有明顯解纏誤差的干涉圖，然后計算形變速率的標準差及掩膜低相干性和高標準差像素點，最后采用時空濾波得到時序震后形變場(Li et al., 2022).

1.2 同震和震后滑動分布反演

根據(jù)同震干涉圖和像素偏移量結(jié)果，將斷層面分成4個部分：發(fā)震斷層西端(斷層1)、中部的主斷裂段(斷層2)、東段分叉兩段(斷層3和斷層4).考慮觀測與模型誤差的不確定度，使用BEAT軟件(Vasyura-Bathke et al. 2020)采用貝葉斯算法來確定發(fā)震斷層主斷裂段的震源參數(shù).在反演過程中，升降軌數(shù)據(jù)采用相同權(quán)重.圖2顯示了貝葉斯搜索結(jié)果直方圖，在95%的置信區(qū)間內(nèi)，最優(yōu)傾角為北傾83°，走向為285°的左旋走滑斷層.在確定主斷裂段斷層參數(shù)后，我們采用SDM(Steepest Descent Method，最速下降法，Wang et al.,2013)基于最小二乘線性搜索其余三段斷層最優(yōu)傾向和傾角.基于模型的粗糙度和失配度搜索最優(yōu)平滑因子為0.13.然后將斷層面離散成2 km×2 km的斷層塊，將滑動角約束在[-45°，45°]范圍內(nèi)，假設(shè)泊松比為0.25來反演每個子斷層塊的滑動，并基于模擬觀測值的均方根誤差搜索最優(yōu)傾角.在整個反演過程中，將先搜索的最優(yōu)斷層段作為約束來搜索之后斷層段.首先，以斷層2作為約束來搜索斷層1的最優(yōu)傾角，其次以斷層1和斷層2作為約束來搜索斷層3的最優(yōu)傾角，最后以斷層1，2，3作為約束搜索斷層4的最優(yōu)傾角.因斷層3和斷層4位于發(fā)震斷層東側(cè)分叉處，反演時存在折衷關(guān)系，最后將斷層3和斷層4一起反演與單獨反演做對比保證反演的精確度.

圖2 斷層幾何參數(shù)的一維和二維后驗概率密度分布圖Fig.2 1-D and 2-D posterior probability density distribution plots of the fault geometry parameters

運動學余滑被廣泛用來觀測早期震后余滑特征，使用同震發(fā)震斷層模型來反演震后滑動分布，將震后斷層模型延伸至30 km寬度，其他斷層幾何參數(shù)與同震發(fā)震斷層相同.

2 結(jié)果和討論

2.1 同震和震后形變場

瑪多地震升降軌同震形變場整體呈東西向分布在斷層兩側(cè)，斷層兩側(cè)最大視線向(LOS)形變量約為1 m(圖3).升降軌相反的形變信號表明瑪多地震以走滑運動為主導.發(fā)震斷層東端分叉處的形變信號階躍表明發(fā)震斷層幾何結(jié)構(gòu)的復雜性，與昆侖斷裂在東端形成多條分支相似(Li et al., 2011; Ceylan et al., 2012).橫貫同震形變場的剖面BB′在斷層兩側(cè)信號突變，表明發(fā)震斷層破裂至地表.圖3g，h所示為Sentinel-1震后4.5個月升降軌震后形變場，升降軌形變信號相反表明震后形變依然以走滑運動為主.

圖4和圖5為震后降軌、升軌時序形變場，震后形變主要分布在斷層兩側(cè)近場區(qū)域，在遠場震后形變量級較小.早期震后形變主要分布在四個區(qū)域，區(qū)域2位于震中西側(cè)同震形變量級較小的區(qū)域，在震后階段呈現(xiàn)明顯的震后形變，震后4.5個月斷層兩側(cè)累積形變量達5 cm.在同震破裂西側(cè)延伸區(qū)域(區(qū)域1)震后最大形變量約為3 cm.而同震時破裂至地表且同震位移最大區(qū)域(區(qū)域3)存在明顯震后形變，最大LOS向形變量達5 cm.總體上，震后形變主要分布在同震形變較小的區(qū)域，而區(qū)域3震后形變與同震形變區(qū)域重疊，表明地震周期內(nèi)能量釋放的連續(xù)性和復雜性.

圖4 瑪多地震Sentinel-1降軌震后時序形變場Fig.4 Time-series of postseismic deformation after the 2021 Madoi earthquake in Sentinel-1 descending orbit

圖5 瑪多地震Sentinel-1升軌震后時序形變場Fig.5 Time-series of postseismic deformation after the 2021 Madoi earthquake in Sentinel-1 ascending orbit

2.2 同震滑動分布

基于非線性和線性反演，得到了瑪多地震發(fā)震斷層的精確幾何結(jié)構(gòu).在非線性反演得到的最優(yōu)主破裂段幾何約束下，斷層1最優(yōu)傾角為北傾75°(圖6a)，分開搜索結(jié)果顯示東端分叉處斷層3最優(yōu)傾角為81°(圖6b)，斷層4最優(yōu)傾角為84°(圖6c)，兩段都為南傾.斷層3和斷層4聯(lián)合反演結(jié)果(圖6d)與單獨搜索結(jié)果相同，確保了斷層傾角搜索結(jié)果的準確性.斷層1和斷層2北傾，斷層3和斷層4南傾，而震后短期余震分布也呈現(xiàn)在斷層中西端位于發(fā)震斷層北側(cè)，而在東端分段處位于斷層南側(cè).余震的空間分布特征驗證了發(fā)震斷層的幾何特征.

圖6 (a) 斷層1最優(yōu)傾角，藍綠色五角星表示最優(yōu)值； (b) 斷層3最優(yōu)傾角； (c) 斷層4最優(yōu)傾角； (d) 斷層3和斷層4聯(lián)合搜索最優(yōu)傾角； (e) 瑪多地震同震滑動分布； (f) 瑪多地震震后滑動分布，藍色圓點表示精定位余震，彩色直線表示斷層地表跡線，紅色五角星代表震中Fig.6 (a) The optimal dip angle of Seg.1, the cyan star represents the optimal value; (b) The optimal dip angle of Seg.3; (c) The optimal dip angle of Seg.4; (d) The optimal dip angle of jointly inversion of Seg.3 and Seg.4; (e) Coseismic slip distribution of the Madoi earthquake; (f) Postseismic slip distribution of the Madoi earthquake, the blue dots represent the aftershocks, the colored lines represent the rupture surface trace, the red star represents the epicenter

圖6e為瑪多地震同震滑動分布，整體呈現(xiàn)長條狀分布在上地殼淺部.同震滑動分布破裂至地表且主要分為四個區(qū)域：斷層2中東部和中西部的淺部滑動；斷層1和斷層3的淺部滑動.最大同震滑動量約為4 m位于斷層2中東部，平均滑動角為5°，同震滑動分布整體呈現(xiàn)左旋走滑運動.模擬殘差圖如圖7所示，較大的殘差值主要分布在斷層近場相干性較差的區(qū)域.假設(shè)剪切模量為30 GPa，瑪多地震同震釋放能量為～1.4446×1020Nm，等同于矩震級MW7.37，與Global Centroid-Moment-Tensor(GCMT)和USGS一致.

圖7 瑪多地震同震形變場(a)(d)(g)(j)分別為ALOS2降軌、ALOS2升軌、Sentinel-1降軌、Sentinel-1升軌觀測值，(b)(e)(h)(k)相應的模擬值，(c)(f)(i)(l)相應的殘差值，白色五角星代表震中.Fig.7 Coseismic deformation field of the Madoi earthquake(a), (d), (g), (j) ALOS2 descending, ALOS2 ascending, Sentinel-1 descending and Sentinel-1 ascending data, respectively; (b), (e), (h), (k) The corresponding modeled data; (c), (f), (i), (l) The corresponding residuals, white star represents the epicenter.

同震滑動分布在上地殼淺部不連續(xù)并存在多個滑動虧損區(qū)域，一處位于斷層2的兩個主要同震滑動區(qū)域之間，以及東端同震滑動未破裂至地表區(qū)域.同震滑動分布周圍未破裂以及滑動虧損區(qū)域通常為速度強化體或者同震時應力未達到破裂閾值的速度弱化體(Yamanakaand Kikuchi, 2004；Miyazaki et al., 2004).因此，滑動虧損區(qū)域以及未破裂區(qū)域在震后階段的活動特征對于區(qū)分它們的物理屬性具有重要意義.

采用Jackknife重采樣技術(shù)(Melgar et al., 2015)來評估斷層滑動模型的穩(wěn)定性.隨機從數(shù)據(jù)集中采取50%的數(shù)據(jù)進行反演，一次反演完成后，將數(shù)據(jù)放回數(shù)據(jù)集再重新隨機抽取50%數(shù)據(jù)進行反演，如此重復進行50次并計算每個斷層塊的標準差.圖8a所示為同震滑動分布的子斷層標準差，最大標準差大約為0.2 m，位于同震最大滑動處，約占同震滑動量級的4%，表明了數(shù)據(jù)與模型的強依賴關(guān)系.

圖8 (a) 同震滑動分布模型不確定度，彩色直線表示斷層地表跡線，紅色五角星代表震中； (b) 震后早期余滑分布模型不確定度； (c) 同震靜態(tài)庫侖應力變化； (d) 發(fā)震斷層速率狀態(tài)摩擦參數(shù)a-bFig.8 (a) Model coseismic slip parameter uncertainties, the colored lines represent the rupture surface trace, the red star represents the epicenter; (b) Model early afterslip parameter uncertainties; (c) Coseismic static Coulomb stress change; (d) The calculated rate and state friction parameter a-b of the seismogenic fault

2.3 震后滑動分布

震后滑動分布采用的升降軌影像對為2021-05-26—2021-10-05，升降軌震后模擬數(shù)據(jù)都能較好的擬合觀測數(shù)據(jù)(圖9).震后余滑分布在多個區(qū)域(圖6f)，在上地殼淺部，余滑整體分布在同震滑動虧損區(qū)域內(nèi)，最大滑動量約為20 cm.在中下地殼存在30 km×10 km的深部余滑，震后4.5個月累積余滑量約為14 cm.最大余滑位于斷層2兩個主要同震滑動區(qū)域之間的滑動虧損區(qū)域內(nèi).而多個余滑區(qū)域顯示較大的同震滑動區(qū)域在震后階段持續(xù)滑動.在發(fā)震斷層東端和西端延伸區(qū)域，余滑分布在斷層14～15 km深度范圍內(nèi)，在斷層3表現(xiàn)為深部余滑.早期震后余滑釋放能量約為5.0417 × 1018Nm，等同于矩震級MW6.4，震后短期內(nèi)的快速余滑是同震應力擾動釋放的主要方式.Heki和Tamura(1997)發(fā)現(xiàn)震后24 h內(nèi)的余滑約占同震滑動的30%，多項研究也表明短期余滑在整個震后周期內(nèi)占據(jù)較大比重(Kim and Dreger, 2008; Yano et al., 2014; Liu and Xu, 2019；Ragon et al., 2019 ).震后第一景影像距離同震發(fā)生時間為震后4天，缺失了短期內(nèi)的震后形變，因此我們得到的震后矩震級低于實際量級.

圖9 瑪多地震震后形變場(a)(d) 分別為Sentinel-1降軌、升軌觀測值; (b)(e) 相應的模擬值；(c)(f)相應的殘差值，白色五角星代表震中.Fig.9 Postseismic deformation field of the Madoi earthquake(a) (d) represent Sentinel-1 descending and ascending data, respectively; (b) (e) The corresponding modeled data; (c) (f) The corresponding residuals, white star represents the epicenter.

采用相同的方法來評估震后數(shù)據(jù)與模型的匹配度，圖8b所示為震后滑動分布的子斷層的均方根誤差，較大的均方根誤差主要分布在斷層淺部震后余滑較小的區(qū)域，整體均方根誤差～1 cm，表明了震后數(shù)據(jù)與模型的強依賴關(guān)系.

2.4 余滑空間分布揭示的摩擦性質(zhì)

震后余滑一般分布在同震滑動區(qū)域周圍，在同震階段未發(fā)生破裂的障礙體上，以及沿發(fā)震斷層走向和傾向延伸的區(qū)域(孫權(quán)等，2022).余滑一般分布在同震庫侖應力加載區(qū)域內(nèi)(Dogan et al., 2014; Wallace et al., 2018).運動學余滑通常被用來觀測短期震后余滑的時空分布特征.本文使用Coulomb 3.3軟件計算了瑪多地震斷層面上的同震庫侖應力變化(圖8c)，余滑整體分布在同震庫侖應力加載區(qū)域內(nèi).

同震滑動區(qū)域是速度弱化體成核區(qū)域，余滑一般圍繞同震滑動區(qū)域分布，主要位于周邊的速度強化體區(qū)域.前面所提及的，上地殼淺部震后余滑部分位于同震滑動量值最大區(qū)域.Yabe和Ida(2018)通過摩擦屬性非均勻分布的斷層面模擬實驗表明，如果同震階段斷層未能充分釋放積累的能量，殘余應力會在震后階段以余滑和余震的方式進行釋放.地下介質(zhì)在空間分布上是非均勻的，同震破裂時快速的滑動可以有效的穿過局部的速度強化體.此外，斷裂帶的摩擦屬性在地震周期中也會發(fā)生變化，當摩擦屬性沿斷層發(fā)生變化并且同時存在速度弱化體和速度強化體，這些斷層區(qū)域會同時參與同震和震后滑動.綜上，我們認為上地殼淺部同震和震后滑動重疊區(qū)域主要是由于殘余應力的釋放驅(qū)動了震后余滑在這些區(qū)域的滑動.相似的震后余滑位于同震滑動區(qū)域在很多地震事件中可以觀測到(Diao et al., 2014; Iinuma et al., 2016; Liu and Xu, 2019).震后余滑的傳播也會對斷層面的應力加載產(chǎn)生影響，進而引發(fā)余震的發(fā)生.余震與余滑在空間上形成很好的互補關(guān)系，基本填充了同震庫侖應力加載區(qū)域.位于深部中下地殼的滑動分布，震后呈現(xiàn)線性滑動狀態(tài)(圖10b)，其中包含了下地殼黏塑性形變機制黏彈性松弛的影響，可能意味著中下地殼較為軟弱的物質(zhì)層的存在，因此之后的研究中應重點關(guān)注震后遠場區(qū)域應力驅(qū)動余滑和黏彈性松弛的折衷關(guān)系，分析青藏高原東北部下地殼和上地幔的流變性.

圖10 (a) 同震、震后早期余滑和余震歸一化累積地震矩隨深度變化； (b) 震后早期余滑地震矩隨深度的時序變化Fig.10 (a) The normalized accumulated seismic moment of coseismic, early afterslip and aftershocks with depth; (b) Temporal evolution of early afterslip seismic moment with depth

2.5 震后余滑的時空分布

時序震后余滑模型顯示余滑量級隨時間逐漸增加.基于震后滑動分布計算的震后時序地震矩(圖10b)顯示，上地殼淺部余滑在2～8 km深度范圍內(nèi)快速衰減，震后3個月后趨于穩(wěn)定.歸一化的同震和震后滑動分布顯示(圖10a)，震后滑動分布在同震滑動量較小的區(qū)域，兩者在空間上形成良好的互補關(guān)系.Gualandi等(2014)表明長期的余滑位于距離同震破裂較遠的區(qū)域.實驗室模擬實驗也表明淺部和深部余滑速率不同，有效正應力是與深度相關(guān)的，深部余滑速率小于淺部余滑(Ariyoshi et al., 2007).這表明同震庫侖應力擾動觸發(fā)了上地殼淺部的速度強化體產(chǎn)生快速余滑，深部余滑隨著距離同震破裂區(qū)域的距離增加趨于穩(wěn)定滑動.對于大地震，余滑可以持續(xù)幾年到幾十年(Barbot et al., 2009)，淺部余滑在震后4.5個月基本釋放，深部將是未來主要余滑區(qū)域.

3 結(jié)論

本文基于ALOS-2和Sentinel-1同震形變場反演獲得了瑪多地震主破裂段斷層幾何參數(shù)，并確定了發(fā)震斷層在西端和東端分叉段的傾角.同震滑動模型顯示斷層位錯主要發(fā)生在上地殼淺部并且破裂至地表，最大滑動量為4 m，位于2 km深處.利用N-SBAS技術(shù)獲得震后4.5個月升降軌時序震后形變場，最大震后形變量近場約為5 cm，遠場約為2 cm，上地殼淺部最大余滑量約為20 cm.震后余滑在上地殼淺部分布在多個區(qū)域，且在深部存在長條狀余滑.時空分布的震后余滑表明上傾角和下傾角區(qū)域不同的摩擦屬性，上地殼淺部在震后階段快速余滑，而深部余滑可能持續(xù)更長時間.

致謝感謝European Space Agency (ESA) 提供的免費Sentinel-1雷達影像數(shù)據(jù)，日本Japan Aerospace Exploration Agency(JAXA)提供的ALOS-2雷達影像數(shù)據(jù)(ER3A2N521).本文所有圖片采用GMT和MATLAB繪制.