姜衛(wèi)平，許才軍，李志偉，武艷強(qiáng)，譚凱，耿江輝，屈春燕，鄭剛，溫?fù)P茂，賀克鋒，周曉慧，劉計(jì)洪，陳長(zhǎng)云，馮光財(cái)，省天琛

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

0 引言

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

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

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

圖1 區(qū)域構(gòu)造背景與公元前780以來(lái)7級(jí)以上強(qiáng)震分布一、二級(jí)活動(dòng)地塊邊界、綠色箭頭(數(shù)字代表地塊運(yùn)動(dòng)速率)依據(jù)張培震等(2003)，活動(dòng)斷裂數(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形變結(jié)果，在前期獲取青海地區(qū)139個(gè)GNSS測(cè)站的基礎(chǔ)上(Zheng et al., 2017)，本文還收集了青海省衛(wèi)星導(dǎo)航定位基準(zhǔn)站網(wǎng)(以下簡(jiǎn)稱青海CORS)的GNSS觀測(cè)數(shù)據(jù)，共包括67個(gè)測(cè)站2017—2021年的觀測(cè)資料.本文采用高精度GNSS數(shù)據(jù)處理開(kāi)源軟件PRIDE PPP-AR處理了這部分GNSS觀測(cè)數(shù)據(jù)，表1給出了具體的數(shù)據(jù)處理策略(Geng et al., 2019；Geng and Mao, 2021)，解算得到了67個(gè)測(cè)站在ITRF2014參考框架下的單天解坐標(biāo)時(shí)間序列.

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

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

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

為了獲取瑪多7.4級(jí)地震的高分辨率InSAR同震形變結(jié)果，本文使用了升降軌Sentinel-1衛(wèi)星和升降軌ALOS-2衛(wèi)星的4對(duì)SAR影像(見(jiàn)表2).通常情況下，Sentinel-1衛(wèi)星以TOPS模式獲取C波段的SAR數(shù)據(jù)，單景SAR數(shù)據(jù)幅寬可達(dá)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ù)的波長(zhǎng)更長(zhǎng)，因此在抵御失相關(guān)噪聲和測(cè)量大梯度形變方面具有較為明顯的優(yōu)勢(shì).

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

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

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

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

圖3 基于升降軌Sentinel-1衛(wèi)星/ALOS-2衛(wèi)星SAR數(shù)據(jù)以及DInSAR/POT/MAI/BOI獲取的青?，敹嗟卣鹧匦l(wèi)星視線向(line-of-sight，LOS)和方位向(azimuth，AZI)的地表形變觀測(cè)值紫色折線代表斷層線，黑色長(zhǎng)箭頭代表對(duì)應(yīng)衛(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)今地殼運(yùn)動(dòng)GNSS速度場(chǎng)紅色箭頭表示本文處理的青海CORS 67個(gè)GNSS測(cè)站的速度，藍(lán)色箭頭表示融合的Zheng 等(2017)中139個(gè)位于青海省的GNSS測(cè)站速度，誤差橢圓表示95%的置信區(qū)間，綠色粗線表示昆侖斷裂，紅色粗線表示昆侖山口—江錯(cuò)斷裂.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)今地殼應(yīng)變率場(chǎng)(a) 應(yīng)變率張量第二不變量分布; (b) 最大剪應(yīng)變率場(chǎng).紅白沙灘球表示近20年發(fā)生在青海省的3次大地震震源機(jī)制解，昆侖斷裂和昆侖山口—江錯(cuò)斷裂加粗顯示.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速度場(chǎng)，應(yīng)用球面近似方程(Savage et al., 2001；Wang and Wright, 2012)得到了應(yīng)變率張量第二不變量分布(圖5a)和最大剪應(yīng)變率場(chǎng)(圖5b).結(jié)果表明：1)應(yīng)變率張量第二不變量分布和最大剪應(yīng)變率場(chǎng)的高值區(qū)基本都位于2001年可可西里地震的發(fā)震斷層——昆侖斷裂西段及附近；2)近20年青海省發(fā)生的3次大地震(2001年可可西里7.8級(jí)地震、2010年玉樹(shù)6.9級(jí)地震、2021年瑪多7.4級(jí)地震)均位于應(yīng)變率高值區(qū)與低值區(qū)的轉(zhuǎn)換地帶.這里需要特別指出，由于本文采用的GNSS資料綜合了近20年觀測(cè)數(shù)據(jù)的平均結(jié)果，盡管在數(shù)據(jù)解算中已經(jīng)對(duì)2001年可可西里7.8級(jí)地震和2010年玉樹(shù)6.9級(jí)地震等強(qiáng)震影響做了扣除，但難免有殘留影響，因此關(guān)于這兩個(gè)地震震源區(qū)的震前變形特征僅供參考.另一方面，由于圖4和圖5涉及的全部GNSS觀測(cè)數(shù)據(jù)均來(lái)源于瑪多7.4級(jí)地震前，因此該地震的震前變形特征對(duì)強(qiáng)震孕育研究具有重要參考意義.

3 瑪多地震三維同震形變場(chǎng)與同震破裂模型

3.1 瑪多地震同震三維形變場(chǎng)分析

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

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

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

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

圖6 青海瑪多地震的同震三維形變場(chǎng)(a) 代表此次地震的同震三維形變場(chǎng)，其中底圖代表垂直向形變，箭頭代表降采樣后的水平形變； (b)和(c)分別代表此次地震的東西向形變和南北向形變，其中圓圈代表GNSS觀測(cè)站點(diǎn).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.

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

圖7 青?，敹嗟卣鹜鸹瑒?dòng)分布Fig.7 Co-seismic slip distribution of the Qinghai Madoi earthquake

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

圖9 InSAR導(dǎo)出和同震模擬的地表三維形變以及三方向殘差第一列到第三列分別為東西、南北和垂直方向形變分量；其中(a—c)為觀測(cè)值，(d—f)為模擬值； (g—i)為觀測(cè)值和模擬值之差.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級(jí)地震以來(lái)，巴顏喀拉地塊及周邊發(fā)生了一系列7級(jí)以上地震，其中2001年昆侖山口西8.1級(jí)和2008年汶川8級(jí)地震的相繼發(fā)生，巴顏喀拉地塊東向運(yùn)動(dòng)增強(qiáng)(江在森等，2009)，導(dǎo)致地塊北邊界以及與之鄰近的近東西向走滑斷裂(比如發(fā)生瑪多地震的昆侖山口—江錯(cuò)斷裂)表現(xiàn)出顯著的左旋運(yùn)動(dòng)加載效應(yīng)，該動(dòng)力效應(yīng)有利于瑪多地震的孕育發(fā)生.

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

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

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

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

5 結(jié)論

GNSS和多源InSAR資料顯示2021年瑪多7.4級(jí)地震導(dǎo)致了昆侖山口—江錯(cuò)斷裂顯著的左旋破裂，最大水平形變發(fā)生在斷層?xùn)|部、量值達(dá)到2.4 m，同時(shí)InSAR資料還表明垂直形變主要發(fā)生在斷層近場(chǎng)，在斷層兩側(cè)呈現(xiàn)正負(fù)交替現(xiàn)象，表明地震破裂的斷層傾角在空間上呈現(xiàn)非均一特征.GNSS和InSAR聯(lián)合反演結(jié)果顯示，此次瑪多地震的斷層滑動(dòng)主要分布在15 km以上，破裂延續(xù)到地表，最大滑動(dòng)量約4 m，矩震級(jí)為MW7.4，破裂沿?cái)鄬幼呦虺尸F(xiàn)明顯的分段特征，形成了多個(gè)顯著破裂區(qū).結(jié)果表明，發(fā)震的江錯(cuò)斷裂地震復(fù)發(fā)周期大約3000 a，青藏高原地塊內(nèi)部斷裂與邊界斷裂一起，分?jǐn)傁《取獨(dú)W亞板塊的匯聚變形.

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