利用空間觀測技術研究青海瑪多7.4級地震孕育發(fā)生變形時空特征

姜衛(wèi)平，許才軍，李志偉，武艷強，譚凱，耿江輝，屈春燕，鄭剛，溫揚茂，賀克鋒，周曉慧，劉計洪，陳長云，馮光財，省天琛

1 武漢大學衛(wèi)星導航定位技術研究中心，武漢 430079 2 武漢大學測繪學院，武漢 430079 3 中南大學地球科學與信息物理學院，長沙 410012 4 中國地震局第一監(jiān)測中心，天津 300180 5 中國地震局地震研究所，武漢 430071 6 中國地震局地質研究所，北京 100029 7 青海省基礎測繪院，西寧 810101

0 引言

2021年5月22日，青?，敹喟l(fā)生了7.4級地震，震中位于北緯34.59°，東經(jīng)98.34°，此次地震延續(xù)了近20年中國大陸7級以上地震活動格局，即圍繞巴顏喀拉地塊的分布特征.作為青藏高原內部活動地塊之一的巴顏喀拉地塊(圖1)，是一個東西方向長近2000 km，而南北方向最窄處僅200 km的長條形地塊，東西橫跨青藏高原腹地(李平恩等，2019).在歐亞參考框架下，南部的羌塘地塊和川滇地塊的東-南東向運動速率明顯快于巴顏喀拉地塊，而巴顏喀拉地塊東向運動又明顯快于其北側的柴達木地塊，因此巴顏喀拉地塊南北邊界以及內部近東西向斷裂均表現(xiàn)出顯著的左旋走滑特征(張培震等，2003).巴顏喀拉地塊北邊界的東昆侖斷裂帶的滑動速率呈現(xiàn)明顯的自西向東衰減特征，其中西大灘—東大灘段全新世晚期滑動速率10～12 mm·a-1(任金衛(wèi)等，1993；青海省地震局和中國地震局地殼應力研究所，1999；Van Der Woerd et al., 2000)；托索湖段滑動速率6～7 mm·a-1(Harkins et al., 2010；Kirby et al., 2007)；下大武段滑動速率為6～7 mm·a-1(青海省地震局和中國地震局地殼應力研究所，1999)；瑪沁—瑪曲段全新世以來的滑動速率3～5 mm·a-1(青海省地震局和中國地震局地殼應力研究所，1999；何文貴等，2006；Kirby et al., 2007);若爾蓋段的滑動速率為1～3 mm·a-1(付俊東，2012)；塔藏段全新世右旋滑動速率0.8 mm·a-1，逆沖速率約為0.3 mm·a-1(Ren et al., 2013).

自20世紀90年代中葉起，青藏高原內部的一系列7級及以上強震都主要發(fā)生在巴顏喀拉地塊及周邊，包括2001年昆侖山口西8.1級、2008年汶川8.0級共2次8級以上地震，以及1997年瑪尼7.5級、2008年新疆于田7.3級、2010年玉樹7.1級、2013年蘆山7.0級、2014年新疆于田7.3級、2017年九寨溝7.0級、2021年瑪多7.4級等共7次7級以上地震(圖1).2021年青?，敹?.4級地震的發(fā)震斷裂為昆侖山口—江錯斷裂，是一條走向北西西的高傾角斷裂，全長約370 km；該斷裂發(fā)育于三疊紀砂巖、板巖中，形成于印支期，斷裂活動強烈，航片衛(wèi)片顯示斷裂在空間展布呈雁列狀分布，左旋位錯現(xiàn)象較明顯，震間GNSS(Global Navigation Satellite System)觀測結果給出江錯斷裂的滑動速率約為1.2±0.8 mm·a-1(朱亞戈等，2021)，古地震研究結果顯示該斷裂的地震復發(fā)周期約3500年(張裕明等，1996).

2021年瑪多7.4級地震是中國大陸繼1947年達日7.7級地震后(劉雷，2021)，又一次發(fā)生在活動地塊內部的地震.地震發(fā)生后，一般測震學最先可以得到其震源機制、余震分布(王未來等，2021)、以及簡單的地震破裂滑動模型.GNSS和InSAR(Interferometric Synthetic Aperture Radar)資料自20世紀90年代以來被廣泛應用于中國大陸地殼形變研究，可以有力推進強震變形與機理研究(單新建等，2004，2017；萬永革等，2008；張培震等，2009；江在森等，2009；Shen et al., 2009；Wang et al., 2011).因此有必要針對瑪多地震震前形變、同震位移分布和發(fā)震斷層位錯開展深入的研究工作.為此，本文系統(tǒng)收集了青海地區(qū)的GNSS觀測數(shù)據(jù)，以及Sentinel-1和ALOS-2(Advanced Land Observation Satellite-2)SAR數(shù)據(jù)，研究獲取了瑪多7.4級地震前GNSS形變場、同震InSAR和GNSS形變場，在此基礎上反演得到了發(fā)震斷層的位錯滑動分布.研究結果為認識板內7級以上強震(特別是活動地塊內部7級以上強震)的孕育發(fā)生機理提供了可靠的大地測量結果與模型基礎.

圖1 區(qū)域構造背景與公元前780以來7級以上強震分布一、二級活動地塊邊界、綠色箭頭(數(shù)字代表地塊運動速率)依據(jù)張培震等(2003)，活動斷裂數(shù)據(jù)依據(jù)鄧起東等(2002).Fig.1 Regional tectonic background and distribution of strong earthquakes at/above magnitude 7 since 780 B.C.Boundaries of first- and second-order active blocks and green arrows (numbers represent movement rates of blocks) are from Zhang et al. (2003), and data of active faults are from Deng et al. (2002).

1 數(shù)據(jù)與方法

1.1 GNSS數(shù)據(jù)與處理方法

為了獲取高空間分辨率的GNSS形變結果，在前期獲取青海地區(qū)139個GNSS測站的基礎上(Zheng et al., 2017)，本文還收集了青海省衛(wèi)星導航定位基準站網(wǎng)(以下簡稱青海CORS)的GNSS觀測數(shù)據(jù)，共包括67個測站2017—2021年的觀測資料.本文采用高精度GNSS數(shù)據(jù)處理開源軟件PRIDE PPP-AR處理了這部分GNSS觀測數(shù)據(jù)，表1給出了具體的數(shù)據(jù)處理策略(Geng et al., 2019；Geng and Mao, 2021)，解算得到了67個測站在ITRF2014參考框架下的單天解坐標時間序列.

表1 青海CORS GNSS數(shù)據(jù)處理策略(Geng et al., 2019; Geng and Mao, 2021)Table 1 The processing strategies of GNSS data from Qinghai CORS (Geng et al., 2019; Geng and Mao, 2021)

為了提高青海CORS站坐標時間序列的可信度，本文對環(huán)境荷載非構造影響進行扣除.具體處理包括如下步驟：利用陸地水儲量模型(包括積雪深度和土壤濕度；https:∥hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH10_3H.2.1/；Rodell et al., 2004；Beaudoing and Rodell, 2020)計算了水文負載形變，時間分辨率為3 h，空間分辨率為1°×1°；采用美國國家環(huán)境預測中心(National Centers for Environmental Prediction, NCEP)發(fā)布的全球地表氣壓產(chǎn)品(https:∥psl.noaa.gov/data/gridded/data.ncep.reanalysis.surface.html；Kalnay et al., 1996)計算了非潮汐大氣壓負載形變，時間分辨率為6 h，空間分辨率為2.5°×2.5°；采用美國國家海洋合作計劃(National Oceanographic Partnership Program, NOPP)制定的海洋環(huán)流及氣候估值(Estimating the Circulation & Climate of the Ocean, ECCO)的全球海底壓力格網(wǎng)數(shù)據(jù)(https:∥ecco-group.org/products.htm；Gross et al., 2005；Gross, 2009)計算了非潮汐海洋負載形變，時間分辨率為12 h，空間分辨率為1°×(0.3°～1°).模型分析表明，負載形變周年信號振幅平均約為GNSS的E、N、U三分量時間序列周年振幅的36.5%、47.6%、63.7%.利用經(jīng)過模型處理后的測站坐標時間序列，本文計算了ITRF2014參考框架下的測站速度，然后采用Altamimi等(2017)給出的ITRF2014歐亞板塊歐拉矢量(Wx=-0.085,Wy=-0.531,Wz=0.770 mas·a-1)將測站速度轉換到穩(wěn)定歐亞參考框架下，圖2為基于青海CORS站的青海省現(xiàn)今地殼運動GNSS速度場.

同時，本文將青海CORS站在瑪多地震震前、震后各3天的數(shù)據(jù)均值做差，解算得到了地震在測站處產(chǎn)生的同震位移，結果如圖8藍色和紅色箭頭所示.結果表明：1)最大同震位移出現(xiàn)在距震中約47 km的JDUO站，約為28.5 cm；東方向最大同震位移也位于JDUO站，約為28.0 cm；北方向最大同震位移則位于距離震中約86 km的KANQ站，約為21.4 cm.2)地震同震影響范圍較大，在距離震中超過350 km的測站處，仍然可以探測到約3 mm的同震位移量.3)同震位移場顯示瑪多地震是一個典型的左旋走滑型地震，斷層上盤從東向西由西南向運動逐漸轉為西北向運動，下盤從西向東由東北向運動逐漸轉為東南向運動.

1.2 InSAR數(shù)據(jù)與處理方法

為了獲取瑪多7.4級地震的高分辨率InSAR同震形變結果，本文使用了升降軌Sentinel-1衛(wèi)星和升降軌ALOS-2衛(wèi)星的4對SAR影像(見表2).通常情況下，Sentinel-1衛(wèi)星以TOPS模式獲取C波段的SAR數(shù)據(jù)，單景SAR數(shù)據(jù)幅寬可達250 km，ALOS-2衛(wèi)星通?？色@取條帶模式(70 km幅寬)和寬幅掃描模式(350 km幅寬)的L波段SAR數(shù)據(jù).相較于C波段SAR數(shù)據(jù)，ALOS-2衛(wèi)星SAR數(shù)據(jù)的波長更長，因此在抵御失相關噪聲和測量大梯度形變方面具有較為明顯的優(yōu)勢.

在SAR數(shù)據(jù)處理方面，本文綜合利用差分干涉合成孔徑雷達測量方法(differential interferometric synthetic aperture radar，DInSAR；Gabriel et al., 1989)、像素偏移量追蹤方法(pixel offset-tracking，POT；Michel et al., 1999)、多孔徑干涉測量方法(multiple aperture interferometry，MAI；Bechor and Zebker, 2006)和子帶重疊區(qū)域干涉測量方法(burst overlap interferometry，BOI；Grandin et al., 2016)處理上述4對SAR影像.具體而言，針對Sentinel-1數(shù)據(jù)，升軌ALOS-2數(shù)據(jù)和降軌ALOS-2數(shù)據(jù)分別采用8×30、8×20和32×5(方位向×距離向)的多視因子，得到了形變觀測值的空間分辨率約為100 m×100 m的形變結果.在DInSAR數(shù)據(jù)處理過程中，采用了一種改進的Goldstein濾波方法(Li et al., 2008)來降低失相干噪聲對結果的影響，然后將相干性小于0.8的像素進行掩膜，利用最小費用流方法(Chen and Zebker, 2002)進行相位解纏.在POT數(shù)據(jù)處理過程中，采用了128×128的匹配窗口進行像素偏移量估計，其中過采樣因子為2；同時，在進行POT處理之前對Sentinel-1數(shù)據(jù)進行了去斜處理(Wegnüller et al., 2016).在ALOS-2數(shù)據(jù)的MAI處理過程中，對全孔徑的SAR數(shù)據(jù)進行前后視分頻后，對前視和后視孔徑的SAR數(shù)據(jù)分別進行標準的DInSAR數(shù)據(jù)處理，然后再對前后視的DInSAR干涉圖進行差分即可得到方位向形變(Bechor and Zebker, 2006).在Sentinel-1數(shù)據(jù)的BOI處理過程中，針對已配準的兩景影像，根據(jù)方位向每一行的時間即可提取出子帶重疊區(qū)域的相鄰子帶的兩個SAR數(shù)據(jù)，這兩個SAR數(shù)據(jù)可認為是MAI中的前后視SAR數(shù)據(jù)，進而類似MAI的數(shù)據(jù)處理過程即可得到BOI結果(Grandin et al., 2016).

圖2 基于青海CORS GNSS觀測數(shù)據(jù)的青海省現(xiàn)今地殼水平運動速度場(a) 負載改正前的速度場; (b) 負載改正后的速度場.Fig.2 Horizontal velocity field of present-day crustal movement of Qinghai Province derived from Qinghai CORS GNSS data(a) shows the velocity field with loading uncorrected and (b) shows the velocity field with loading corrected.

表2 本文所用的Sentinel-1衛(wèi)星和ALOS-2衛(wèi)星SAR數(shù)據(jù)的基本信息Table 2 The outline of the SAR data from Sentinel-1 and ALOS-2 satellites used in this study

圖3為本文獲取的青?，敹嗟卣?6個SAR形變場結果，可以發(fā)現(xiàn)成像幾何相似的SAR形變場的形變空間特征較為相似.例如，降軌Sentinel-1數(shù)據(jù)的DInSAR形變場(圖3e)、降軌Sentinel-1數(shù)據(jù)的POT方法LOS向形變場(圖3f)和降軌ALOS-2數(shù)據(jù)的DInSAR形變場(圖3m)較為相似.由于ALOS-2數(shù)據(jù)比Sentinel-1數(shù)據(jù)的波長更長，因此，相對于Sentinel-1數(shù)據(jù)(圖3a, e)，基于ALOS-2數(shù)據(jù)的DInSAR方法(圖3i,m)可測量的形變梯度較高，即ALOS-2數(shù)據(jù)可獲得更為完整的近場同震形變結果.然而，從ALOS-2數(shù)據(jù)的方位向形變結果(圖3k,l,o,p)可以看出，ALOS-2數(shù)據(jù)受到了嚴重的電離層擾動影響，以至于相應的方位向觀測數(shù)據(jù)中的形變信號被電離層信號完全掩蓋，無法用于后續(xù)的形變分析和地震解譯.另外，電離層擾動也會影響LOS向的觀測結果.本文采用了距離向譜分解(range spectrum split，RSS)方法(Gomba et al., 2016)用來估計并去除ALOS-2數(shù)據(jù)中LOS向觀測值的電離層相位.對于POT形變觀測值而言，其測量精度明顯低于DInSAR和BOI方法.尤其對于Sentinel-1數(shù)據(jù)POT方法的方位向形變(圖3c,g)和降軌ALOS-2數(shù)據(jù)POT方法的LOS向形變(圖3n)，其相應的空間分辨率均較低(見表2)，因此相應的形變場包含了較為嚴重的噪聲信號.盡管BOI方法僅能獲取子帶重疊區(qū)域的方位向形變結果(圖3d,h)，但對于青海瑪多地震而言，BOI方法獲取的形變結果很好地揭示了斷層兩側的相對運動，為后續(xù)分析解譯提供了重要的數(shù)據(jù)基礎.

2 瑪多地震前地殼變形特征分析

在本文處理的67個青海CORS站速度場結果(圖2)的基礎上，進一步融合了前期得到的青海省139個GNSS測站速度(Zheng et al., 2017).青海CORS站東向和北向速度平均精度分別約為0.62 mm·a-1和0.75 mm·a-1，而前期GNSS測站東向和北向速度平均精度分別約為1.02 mm·a-1和0.86 mm·a-1，均略低于青海CORS站的精度.雖然前期青海地區(qū)GNSS數(shù)據(jù)的時間跨度為1998—2015年(Zheng et al., 2017)，跨度較長，而本文青海CORS數(shù)據(jù)的時間跨度為2017—2021年，跨度較短，但是前期GNSS數(shù)據(jù)大部分為流動觀測，平均間隔1～2年才有數(shù)天觀測，而本文青海CORS數(shù)據(jù)均為連續(xù)觀測，所以得到的測站速度精度會相對高一些.本文青海CORS站速度場和前期GNSS速度場均轉換到了穩(wěn)定歐亞參考框架下，但是不同數(shù)據(jù)源之間可能會存在參考框架差異，如果兩個數(shù)據(jù)源有適量的公共數(shù)據(jù)，則可以利用公共數(shù)據(jù)對差異進行擬合，但是本文兩個數(shù)據(jù)源之間沒有公共數(shù)據(jù)，而考慮到這種差異對速度場帶來的影響通常很小，一般小于速度場本身精度(Zheng et al., 2017)，所以本文忽略了參考框架差異對速度場融合帶來的影響.圖4給出了最終的穩(wěn)定歐亞參考框架下總共包含206個測站的瑪多地震震前GNSS速度場.

圖3 基于升降軌Sentinel-1衛(wèi)星/ALOS-2衛(wèi)星SAR數(shù)據(jù)以及DInSAR/POT/MAI/BOI獲取的青海瑪多地震沿衛(wèi)星視線向(line-of-sight，LOS)和方位向(azimuth，AZI)的地表形變觀測值紫色折線代表斷層線，黑色長箭頭代表對應衛(wèi)星飛行方向，紅色短箭頭代表LOS方向.Fig.3 The surface deformation fields of the Qinghai Madoi earthquake in the line-of-sight and azimuth directions derived from the ascending- and descending-track Sentinel-1 and ALOS-2 SAR data by DInSAR/POT/MAI/BOI methods, respectivelyThe purple line represents the fault; the black long arrow shows the azimuth direction; the red short arrow shows the line-of-sight direction.

圖4 穩(wěn)定歐亞參考框架下的青海省現(xiàn)今地殼運動GNSS速度場紅色箭頭表示本文處理的青海CORS 67個GNSS測站的速度，藍色箭頭表示融合的Zheng 等(2017)中139個位于青海省的GNSS測站速度，誤差橢圓表示95%的置信區(qū)間，綠色粗線表示昆侖斷裂，紅色粗線表示昆侖山口—江錯斷裂.Fig.4 GNSS velocity field of present-day crustal movement of Qinghai Province with respect to stable EurasiaThe red arrows show the velocities at 67 Qinghai CORS GNSS stations processed in this study, and the blue ones show the velocities at 139 GNSS stations taken from Zheng et al. (2017). Error ellipses are 95% confidence interval. The green thick line shows the Kunlun fault and the red thick line shows the KunlunshanPass-Jiangcuo fault.

圖5 青海省現(xiàn)今地殼應變率場(a) 應變率張量第二不變量分布; (b) 最大剪應變率場.紅白沙灘球表示近20年發(fā)生在青海省的3次大地震震源機制解，昆侖斷裂和昆侖山口—江錯斷裂加粗顯示.Fig.5 Present-day strain rate fields in Qinghai Province(a) Distribution of second invariant of strain rate tensor; (b) Maximum shear strain rate field. The beach-balls show the focal mechanisms of three large earthquakes occurred in Qinghai Province in the recent 20 years. The Kunlun and KunlunshanPass-Jiangcuo faults are shown in bold.

為了直觀給出瑪多震源區(qū)附近的變形模式和量值，基于圖4所示GNSS速度場，應用球面近似方程(Savage et al., 2001；Wang and Wright, 2012)得到了應變率張量第二不變量分布(圖5a)和最大剪應變率場(圖5b).結果表明：1)應變率張量第二不變量分布和最大剪應變率場的高值區(qū)基本都位于2001年可可西里地震的發(fā)震斷層——昆侖斷裂西段及附近；2)近20年青海省發(fā)生的3次大地震(2001年可可西里7.8級地震、2010年玉樹6.9級地震、2021年瑪多7.4級地震)均位于應變率高值區(qū)與低值區(qū)的轉換地帶.這里需要特別指出，由于本文采用的GNSS資料綜合了近20年觀測數(shù)據(jù)的平均結果，盡管在數(shù)據(jù)解算中已經(jīng)對2001年可可西里7.8級地震和2010年玉樹6.9級地震等強震影響做了扣除，但難免有殘留影響，因此關于這兩個地震震源區(qū)的震前變形特征僅供參考.另一方面，由于圖4和圖5涉及的全部GNSS觀測數(shù)據(jù)均來源于瑪多7.4級地震前，因此該地震的震前變形特征對強震孕育研究具有重要參考意義.

3 瑪多地震三維同震形變場與同震破裂模型

3.1 瑪多地震同震三維形變場分析

本文用InSAR資料獲取了豐富的視線向(LOS)同震形變場，但只是一維形變.GNSS同震形變較少，且垂直形變精度較低.為了直觀，也為了獲得較豐富的三維同震形變，我們將瑪多地震InSAR和GNSS同震形變融合成高精度同震三維形變場.本文使用了一種基于地表應變模型(strain model，SM)和方差分量估計(variance component estimation，VCE)的InSAR三維地表形變分解方法(SM-VCE；Liu et al., 2018).鑒于ALOS-2數(shù)據(jù)的方位向形變觀測結果(圖3k,l,o,p)和降軌Sentinel-1數(shù)據(jù)POT方法的方位向形變觀測結果(圖3g)受電離層擾動和失相干噪聲的影響較為嚴重，本文僅利用剩余的11個SAR形變觀測數(shù)據(jù)進行分析，圖6給出了基于SM-VCE方法獲取的青?，敹嗟卣鸬耐鹑S形變場.結果顯示，此次地震的整體形變符合左旋走滑斷層的形變特征，最大水平形變發(fā)生在破裂斷層的東部，量值高達2.4 m.在震中附近的東西向形變量級較小，沿斷層軌跡從震中向東和向西區(qū)域的東西向形變量級較大，并且在斷層東端北側附近，水平形變方向從西南變?yōu)槲鞅?與水平形變相比，垂直向形變主要發(fā)生在近斷層區(qū)域，且形變量級較小(-0.52～0.64 m)，并且垂直向形變沿斷層走向在斷層兩側出現(xiàn)了正負交替的現(xiàn)象，顯示此次地震破裂的斷層傾角存在非均一性特征.總體上，瑪多7.4級地震的三維形變場結果顯示，大部分形變主要發(fā)生在巴顏喀拉地塊內部，地塊北邊界的東昆侖斷裂帶、南邊界的甘孜—玉樹以及鮮水河斷裂帶的形變不明顯.

對比青海瑪多地震近場GNSS的水平同震位移結果(圖6b—c的圓圈處)，InSAR與GNSS東西向和南北向形變差值的均方根誤差(root mean square error，RMSE)分別為5.2 cm和12.2 cm，其中InSAR和GNSS在震中最近處測點南北向形變之差高達35 cm，直接導致了南北向形變的RMSE值較大.原因可能是在斷層附近破裂嚴重，導致GNSS站點發(fā)生了局部劇烈形變，也可能是因為InSAR結果的南北向形變幾乎全部由精度較低的方位向形變觀測數(shù)據(jù)得到所致，SAR南北向形變本身的精度較低.在剔除此GNSS站點之后，重新計算的東西向和南北向RMSE數(shù)值分別為5.2 cm和6.4 cm，說明本文獲取的青?，敹嗟卣鹑S形變場是可靠的.

3.2 三維同震形變約束下的同震破裂模型

為獲得瑪多地震同震破裂的精細幾何結構和同震滑動分布，本文首先基于InSAR測量、現(xiàn)場勘察(潘家偉等，2021)和重定位余震(王未來等，2021)構建了一個傾向沿走向變化的三維同震破裂幾何模型，斷層整體為高傾角特征，從西北端的西南傾沿走向向東逐漸轉為高傾角的東北傾.隨后，將斷層劃分為3 km×3 km的均勻三角單元，以InSAR和GNSS觀測合成的三維同震形變數(shù)據(jù)為約束，進行了瑪多地震同震破裂的滑動分布反演.三維同震形變具有直觀、水平形變和垂直形變精度都很高的優(yōu)點，為減少冗余和提高效率，采用基于斷層模型的自適應降采樣方法(Lohman and Simons, 2005)降采樣InSAR三維形變.同時在反演過程中為避免得到一個力學不可信且過于粗糙的滑動分布結果，對滑動分布反演施加了一個改進的平滑約束(Maerten et al., 2005)，既考慮了主-分支斷層的幾何關系，又有效避免了斷層淺部由于傳統(tǒng)平滑約束所導致的滑動虧損.反演在剪切模量為30 GPa、泊松比為0.25的均勻彈性半空間內進行(Meade, 2007).此外，考慮到瑪多地震走滑為主的破裂特征(潘家偉等，2021)以及余震分布特點(王未來等，2021)，反演中僅允許滑動向量的滑動角在-50°～50°范圍內變化，整個反演過程依據(jù)模擬殘余和模型粗糙度的L曲線來確定一個合理的平滑因子，以至于使滑動分布模型在不至于過分粗糙的情況下，能夠很好地擬合觀測數(shù)據(jù).

圖6 青海瑪多地震的同震三維形變場(a) 代表此次地震的同震三維形變場，其中底圖代表垂直向形變，箭頭代表降采樣后的水平形變； (b)和(c)分別代表此次地震的東西向形變和南北向形變，其中圓圈代表GNSS觀測站點.Fig.6 The 3-D co-seismic deformation field of the Qinghai Madoi earthquake(a)shows the 3-D co-seismic deformation field, in which the base map shows the vertical deformation, and the arrows show the down-sampling horizontal deformation. (b) and (c) show the deformation fields in the E-W and N-S directions, respectively. The circles in (b) and (c) represent GNSS stations.

模擬結果顯示瑪多地震破裂至地表，滑動主要分布在地表至15 km深度，最大滑動位于近地表且在地表沒有明顯的弱化(圖7)，最大滑動量值約4 m，釋放累積矩1.78×1020Nm，對應矩震級MW7.4，滑動分布沿斷層呈現(xiàn)峰-谷交替的特征，顯示斷層破裂具有明顯的橫向不均勻性.傾滑分量主要發(fā)生在斷層走向轉換區(qū)，且這些轉換區(qū)主要對應于拉分階區(qū)，表明斷層幾何和滑動分布有很好的空間對應關系.余震主要分布于同震滑動分布高值區(qū)周緣，和高滑動分布位置呈現(xiàn)空間互補的特征.正演GNSS和InSAR三維形變(圖8和圖9)，可以發(fā)現(xiàn)模擬值和觀測值有很好的對應性，對應的加權殘余誤差分別為4.9 mm和15.9 cm.表明滑動分布結果可以有效再現(xiàn)GNSS和InSAR兩個數(shù)據(jù)集的主形變特征.兩個臺站GPS形變信號擬合較差的主要原因可能是由于這兩個站點建立在非基巖上，因此信號受到場地效應的影響，致使其偏離了同震形變主趨勢.InSAR三維形變擬合相對更差的原因可能是三維形變導出誤差和InSAR觀測原始數(shù)據(jù)對應解纏和大氣擾動等因素共同作用所致.InSAR的殘差較大測點主要集中在近斷層位置，推測可能是由于彈性模型無法解釋近斷層的非彈性形變或主斷層形變模式無法模擬近斷層形變分布所致(潘家偉等，2021).此外，對比原始InSAR LOS觀測和GPS聯(lián)合反演結果，發(fā)現(xiàn)兩者具有相似的滑動分布模式，意味著這兩個模型具有很好的一致性，均能較好的解釋相應的觀測資料.相較于原始InSAR LOS觀測直接用于滑動分布反演，InSAR三維形變數(shù)據(jù)由于引入了物理約束(應變)，可以有效避免原始觀測數(shù)據(jù)集中局部噪聲對滑動分布結果的影響，雖然這可能會導致部分區(qū)域信號被平滑.此外，由于InSAR三維形變數(shù)據(jù)有效避免了多數(shù)據(jù)集之間的定權問題，因此也就降低了不合理定權對反演結果的影響.

圖7 青海瑪多地震同震滑動分布Fig.7 Co-seismic slip distribution of the Qinghai Madoi earthquake

圖8 瑪多地震GNSS同震形變觀測值和模擬值Fig.8 Observed and modeled values of GNSS co-seismic deformation of the Madoi earthquake

圖9 InSAR導出和同震模擬的地表三維形變以及三方向殘差第一列到第三列分別為東西、南北和垂直方向形變分量；其中(a—c)為觀測值，(d—f)為模擬值； (g—i)為觀測值和模擬值之差.Fig.9 3-D InSAR-derived and modeled surface deformation fields and residuals in three directionsThe subfigures in the first, second and third columns show the deformation in EW, NS and UP directions, respectively; (a—c) show the observed values; (d—f) show the modeled values; (g—i) show the differences between the observed and modeled values.

4 討論

自1997年瑪尼7.5級地震以來，巴顏喀拉地塊及周邊發(fā)生了一系列7級以上地震，其中2001年昆侖山口西8.1級和2008年汶川8級地震的相繼發(fā)生，巴顏喀拉地塊東向運動增強(江在森等，2009)，導致地塊北邊界以及與之鄰近的近東西向走滑斷裂(比如發(fā)生瑪多地震的昆侖山口—江錯斷裂)表現(xiàn)出顯著的左旋運動加載效應，該動力效應有利于瑪多地震的孕育發(fā)生.

從長期地震危險性研究角度，瑪多7.4級地震發(fā)生在中國地震局M7專項工作組給出的昆侖山—瑪多危險區(qū)內(M7專項工作組，2012)，潛在的發(fā)震斷裂主要為東昆侖斷裂帶西大灘段和昆侖山口—達日斷裂等，判定依據(jù)主要為大型活動斷裂帶上的地震空區(qū)、2001年昆侖山口西地震破裂演化以及歷史大地震的時間-序次關系等.從形變場特征與強震發(fā)生地點角度來看，此次瑪多7.4級地震與2001年以來其他4次7級以上強震(2001年昆侖山口西8.1級、2010年玉樹7.1級、2014年新疆于田7.3級、2017年九寨溝7.0級)具有一致性，均發(fā)生在最大剪應變率高值區(qū)邊緣(武艷強等，2020).具體到發(fā)震斷層尺度，GNSS速度剖面結果顯示瑪多地震發(fā)震斷層的滑動速率明顯小于其北側東昆侖斷裂帶托索湖段(朱亞戈等，2021)，同時瑪多地震發(fā)震斷層兩側連續(xù)變形更緩慢，反映了一定程度的孕震晚期特征(Savage and Burford,1973；Meade and Hager, 2005).

地震發(fā)生后，測震學最先用地震波數(shù)據(jù)約束反演獲得震源機制以及破裂滑動分布模型，為地震危險性評估提供快速參考.USGS基于震源機制的走向106°、傾角76°構建了一個平板斷層面，用地震波約束反演了有限斷層滑動分布，矩震級MW7.4、最大滑動量4 m的數(shù)值與本文結果一致，但其最大滑動量位于本文最大滑動量西側約50 km處，從最大滑動處向兩側分散.Wang等(2021)雖然在地震波的基礎加了Sentinel InSAR數(shù)據(jù)約束，但還是沒有考慮斷裂東側的破裂情況，其最大滑動達7 m，位于本文最大滑動東側幾十公里處，埋深較大，可能更多的受到遠場地震波的影響.Chen等(2021)用升降軌Sentinel數(shù)據(jù)約束的模型最大滑動5 m，最大滑動位置和深度與本文比較接近.本文用近場更多的升降軌Sentinel和ALOS InSAR數(shù)據(jù)和GNSS數(shù)據(jù)分辨出5個以上滑動峰值區(qū)，揭示了更多的破裂細節(jié)和斷裂東側的分叉破裂，與余震活動分布一致，可能更多的受到近場大地測量形變的影響.

本文的滑動分布反演結果顯示2021年青?，敹郙W7.4地震發(fā)生在昆侖山口—江錯斷裂帶，破裂長度～160 km，最大破裂位于近地表，且滑動分布在地表未見明顯的滑動衰減.這一滑動特征在合成3D形變中也有所反映，形變最大值位于毗連斷層的兩側，顯示斷層滑動分布最大值至少位于近地表，否則由于淺地表的空間濾波效應，地表形變將更加連續(xù)平滑且形變峰值將偏離斷層一定距離.這一特征區(qū)別于常見的伴隨非成熟斷層破裂產(chǎn)生的淺部滑動虧損(shallow slip deficiency，SSD；例如：2019年Ridgecrest地震)，表明瑪多地震發(fā)震斷層(昆侖山口—江錯斷裂)是一條發(fā)育成熟的斷裂帶.由于合成3D形變可能導致形變在空間上相較于多數(shù)據(jù)集原始形變更加連續(xù)平滑，但基本的破裂滑動分布特征仍然可以假定其穩(wěn)定性，如最大凹凸體位置和大小，滑動分布沿走向變化的整體空間特征.結合地質現(xiàn)場勘察(潘家偉等，2021)和動態(tài)破裂反演(USGS, 2021)，我們可以發(fā)現(xiàn)瑪多地震成核于近階區(qū)位置，而進一步沿直且連續(xù)的簡單斷層幾何傳播，且在該簡單幾何區(qū)間形成最大破裂，最后破裂終止于兩個釋放階區(qū)，符合典型的走滑地震成核-破裂-傳播-終止特征(如：1999年Hector Mine地震).此外，瑪多地震在斷層兩側釋放階區(qū)產(chǎn)生大的滑動集中區(qū)可能是由于釋放階區(qū)斷層一般相對更弱，更容易產(chǎn)生分支斷裂和發(fā)生中小能量的集中釋放(Scholz, 2019).

瑪多地震同震變形結果顯示，此次地震的應變釋放主要集中在巴顏喀拉地塊內部的狹長區(qū)域.地震斷層破裂滑動最大達4 m，若以斷層滑動1.2 mm·a-1計算，則地震復發(fā)周期是3333 a，與古地震研究結果接近.而據(jù)GNSS速度場剖面，巴顏喀拉地塊的邊界斷層滑動速率一般在3～12 mm·a-1，地震復發(fā)周期更短.發(fā)震斷裂北側的東昆侖斷裂帶托索湖段和瑪沁—瑪曲段對地震的形變響應不明顯，但瑪多地震對該斷裂、龍日壩斷裂帶以及龍門山斷裂帶南段呈現(xiàn)明顯的應力加載效應.另一方面，庫侖應力計算結果同樣表明，東昆侖斷裂帶的托索湖段與瑪沁—瑪曲段、龍日壩斷裂帶為庫侖應力增加段(Li et al., 2021).因此，綜合斷裂帶對瑪多地震的同震響應特征以及庫侖應力影響，需重點關注上述斷裂帶的強震緊迫性.

5 結論

GNSS和多源InSAR資料顯示2021年瑪多7.4級地震導致了昆侖山口—江錯斷裂顯著的左旋破裂，最大水平形變發(fā)生在斷層東部、量值達到2.4 m，同時InSAR資料還表明垂直形變主要發(fā)生在斷層近場，在斷層兩側呈現(xiàn)正負交替現(xiàn)象，表明地震破裂的斷層傾角在空間上呈現(xiàn)非均一特征.GNSS和InSAR聯(lián)合反演結果顯示，此次瑪多地震的斷層滑動主要分布在15 km以上，破裂延續(xù)到地表，最大滑動量約4 m，矩震級為MW7.4，破裂沿斷層走向呈現(xiàn)明顯的分段特征，形成了多個顯著破裂區(qū).結果表明，發(fā)震的江錯斷裂地震復發(fā)周期大約3000 a，青藏高原地塊內部斷裂與邊界斷裂一起，分攤消化印度—歐亞板塊的匯聚變形.

致謝感謝“中國地殼運動觀測網(wǎng)絡”和青海省基礎測繪院提供GNSS數(shù)據(jù).感謝日本宇航局提供的ALOS-2數(shù)據(jù)(PER2A2N038)，感謝歐空局提供的Sentinel-1數(shù)據(jù).