2021年青?，敹郙W7.4地震GNSS同震形變場及其斷層滑動分布

王迪晉，王東振，趙斌*，李瑜，趙利江，王閱兵，聶兆生，喬學(xué)軍，王琪

1 中國地震局地震研究所，武漢 430071 2 中國地震局地震大地測量重點實驗室，武漢 430071 3 中國地震臺網(wǎng)中心，北京 100045 4 青海省地理信息與自然資源綜合調(diào)查中心，西寧 810001 5 中國地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，武漢 430074

0 引言

據(jù)中國地震臺網(wǎng)中心測定，北京時間2021年5月22日2時4分在我國青海省瑪多縣境內(nèi)發(fā)生MS7.4地震，震中位于34.59°N，98.34°E，震源深度17 km.美國地質(zhì)調(diào)查局(USGS)給出的震中位于瑪多縣黃河鄉(xiāng)東側(cè)(34.598°N，98.251°E)，震級為MW7.3，震源深度10 km，可能的發(fā)震斷層為一條近EW走向、帶有正斷分量的左旋走滑斷層(https:∥earthquake.usgs.gov/earthquakes/eventpage/us7000e54r/).地表破裂考察、余震分布、InSAR干涉影像等結(jié)果顯示，瑪多地震的發(fā)震斷層為昆侖山口—江錯斷裂，地表破裂長約160 km，可分成四個段落，自西向東分別為鄂陵湖南段、黃河鄉(xiāng)段、冬草阿龍湖段和昌麻河鄉(xiāng)段(李智敏等，2021；王未來等，2021；華俊等，2021).昆侖山口—江錯斷裂是巴顏喀拉塊體的內(nèi)部斷層，也是東昆侖斷裂帶的次級分支構(gòu)造之一(李陳俠等，2011).近30年來，中國大陸7級以上的地震事件都發(fā)生于巴顏喀拉塊體周緣斷裂帶上(鄧起東，2012；鄧起東等，2014)(圖1).此次地震是中國大陸地區(qū)繼2008年汶川地震以來震級最大的一次構(gòu)造地震，也是2017年九寨溝MS7.0地震后又一次震級大于7級的板內(nèi)地震.該地震的發(fā)生再次引起人們對巴顏喀拉塊體內(nèi)部及其邊界斷裂帶上的地震活動性的強烈關(guān)注.

巴顏喀拉塊體位于青藏高原北部，是青藏高原內(nèi)部物質(zhì)東向擠出的重要構(gòu)造單元之一.該塊體的北邊界為東昆侖斷裂；南邊界為甘孜—玉樹斷裂；東邊界為龍門山斷裂；西邊界與阿爾金斷裂帶交匯，這些邊界的構(gòu)造活動共同控制著塊體的現(xiàn)今變形.現(xiàn)今地殼速度場顯示青藏高原東北緣整體向北東方向運動(圖1)，巴顏喀拉塊體內(nèi)部區(qū)域形變相對較小.伴隨著巴顏喀拉塊體東向運移速率自西向東逐漸衰減，北邊界東昆侖斷裂帶的左旋走滑速率也隨之衰減(Kirby et al., 2007; Lin and Guo, 2008)，在該斷裂帶的東段發(fā)育了多條次級分支斷裂，如甘徳—瑪多斷裂、昆侖山口—江錯斷裂和達(dá)日斷裂等(圖1)，其中1937年在主斷裂托索湖段發(fā)生M7.5地震(Molnar and Deng, 1984；沈正康等，2003)、1947年在達(dá)日斷裂上發(fā)生MS7.7地震(Molnar and Deng, 1984；萬永革等，2007).地理位置上，瑪多地震同樣發(fā)生在靠近東昆侖斷裂帶的區(qū)域，發(fā)震斷層與相鄰的其他次級分支斷裂一起構(gòu)成了巴顏喀拉塊體北部寬闊的彌散型邊界帶(潘家偉等，2021).震間GNSS速度剖面的結(jié)果顯示巴顏喀拉塊體北部邊界斷裂帶兩側(cè)出現(xiàn)明顯的速度梯度變化，而本次地震發(fā)震斷層的滑動速率僅為1.2 ± 0.8 mm·a-1，比主邊界帶東昆侖斷裂托索湖段的滑動速率小(朱亞戈等，2021)，表明發(fā)震斷層昆侖山口—江錯斷裂應(yīng)力積累較慢，復(fù)發(fā)周期長，而這次地震很可能是長期震間構(gòu)造應(yīng)力積累下的自發(fā)破裂事件.

瑪多地震在巴顏喀拉塊體內(nèi)部形成了怎樣的形變特征？發(fā)震斷層的同震滑動分布如何？對周緣斷裂會造成哪些潛在的影響？這些都是亟待回答的科學(xué)問題.李智敏等(2021)通過詳細(xì)的野外地震考察，獲得了地震地表破裂組合特征和運動學(xué)性質(zhì)、破裂分段性、同震位移特征及分布規(guī)律，并根據(jù)沖溝、道路和拉張階區(qū)裂隙寬度確定地表同震位移量為1～2 m.潘家偉等(2021)通過經(jīng)驗公式估算該次地震的最大位移量約為4 m，并對發(fā)震構(gòu)造進(jìn)行了分析和探討，同時論述了發(fā)震斷裂與東昆侖斷裂帶的關(guān)系、巴顏喀拉地塊構(gòu)造變形的動力學(xué)機制以及區(qū)域未來強震活動性等問題.華俊等(2021)基于歐洲航空局Sentinel-1衛(wèi)星C波段的升、降軌數(shù)據(jù)獲取了瑪多地震InSAR同震形變場，確定了地震的地表破裂跡線，并結(jié)合庫侖應(yīng)力擾動分析了瑪多地震引起的周邊區(qū)域的地震危險性，反演得到的矩震級為MW7.45，斷層最大滑動量約為6 m，主體破裂集中在0～10 km深度.但由于近場失相干效應(yīng)明顯，利用InSAR資料得到的同震形變場缺少極震區(qū)的近場觀測，這也導(dǎo)致了由此反演得到的震級和最大滑動量均略大于地震波的結(jié)果.

圖1 (a) 青藏高原背景速度場； (b) 2021年瑪多地震的構(gòu)造背景圖(a)中的黑色粗實線給出了巴顏喀拉塊體的主邊界，黑色細(xì)實線給出了東昆侖斷裂附近發(fā)育的多條次級分支斷裂，紅色實線表示瑪多地震發(fā)震斷層的地表破裂軌跡，虛線矩形框指示了圖(b)的位置； 圖(b)中綠色圓點給出了余震序列的精定位結(jié)果(王未來等，2021).圖中的速度場結(jié)果引自文獻(xiàn)(Wang and Shen, 2020).Fig.1 (a) Background velocity field of the Tibetan plateau; (b) Tectonic setting of the 2021 Madoi earthquakeIn figure (a), the thick black solid line denotes the main boundary of the Bayan Har massif, the thin black solid lines denote secondary branch faults around the East Kunlun fault, and the red solid line denotes surface track of the rupture of the Madoi earthquake. The dashed rectangle box indicates the position of figure (b). The green circles in figure (b) represent relocated aftershocks (Wang et al., 2021). The velocity field in the figures is cited from (Wang and Shen, 2020).

隨著GNSS技術(shù)的迅速發(fā)展，該項技術(shù)已經(jīng)成為現(xiàn)代地震研究領(lǐng)域中的一個重要觀測手段(Bock and Melgar, 2016).尤其是近場GNSS觀測可以彌補InSAR影像在近場失相干區(qū)域的觀測不足，直觀詮釋發(fā)震機制，同時為地震動態(tài)破裂過程和靜態(tài)破裂分布模型提供重要的約束.李志才等(2021)展示了破裂帶附近9個GNSS基準(zhǔn)站的高頻(1 Hz采樣)動態(tài)形變波形，及12個GNSS連續(xù)站的同震形變場，觀測到的最大水平同震位移0.6 m，距離斷層20 km處.該研究雖然初步揭示了本次地震的形變特征，但由于點位稀疏，對斷層滑動分布的約束不足，破裂模型稍顯粗糙.

本文利用多家單位共同完成的GNSS觀測數(shù)據(jù)，在給出更為全面的同震形變場和高頻GNSS動態(tài)波形的基礎(chǔ)上，分析了本次地震的同震形變特征及其對周邊斷裂的影響，比較了GNSS觀測與InSAR衛(wèi)星影像的同震形變結(jié)果，反演得到精細(xì)的地震破裂模型，最后計算了地震在發(fā)震斷層面上產(chǎn)生的庫侖應(yīng)力變化，與余震精定位結(jié)果進(jìn)行了比較，并分析了巴顏喀拉塊體相關(guān)區(qū)域在未來的地震危險性.

1 變形觀測與數(shù)據(jù)處理

1.1 觀測資料

瑪多地震發(fā)生以后，中國地震局地震研究所第一時間派工作組趕赴震區(qū)開展GNSS應(yīng)急流動觀測.利用前期國內(nèi)不同機構(gòu)、不同研究團隊在震區(qū)布設(shè)的地殼形變觀測站點，在地震發(fā)生后十天內(nèi)完成了53個站點的流動觀測(每站觀測24～48 h).這些流動站主要包含陸態(tài)網(wǎng)絡(luò)工程的區(qū)域站，原國家測繪局B級網(wǎng)點，自然資源部于2019年新建的國家大地控制點和少量的三角點.此外，還有相關(guān)科研機構(gòu)在震中附近建立的流動GNSS站點，主要包括了東昆侖斷裂托索湖段近南北向GNSS剖線(正好位于發(fā)震斷層北側(cè)，Diao et al., 2019)和中國地震局地震研究所在該斷裂瑪沁—瑪曲段布設(shè)的一條GNSS剖線.

我們還收集了震區(qū)300 km范圍內(nèi)的53個GNSS連續(xù)觀測站，主要包括陸態(tài)網(wǎng)絡(luò)工程的連續(xù)觀測站和青海省的CORS站(衛(wèi)星導(dǎo)航與位置服務(wù)基準(zhǔn)站).其中部分連續(xù)站點同時還具有高頻采樣數(shù)據(jù)，因此我們也收集到了其中17個臺站1 Hz采樣的高頻GNSS觀測數(shù)據(jù).圖2展示了所有站點的分布圖.

1.2 靜態(tài)形變場處理

本文所有的GNSS數(shù)據(jù)處理均采用GAMIT/GLOBK10.6軟件 (Herring et al., 2015)，數(shù)據(jù)處理策略與流程與Zhao等(2015)一致，對固體潮、極潮、海潮進(jìn)行改正，并應(yīng)用最新的衛(wèi)星、天線絕對相位中心改正模型和全球氣壓溫度模型GPT2(Lagler et al., 2013)，估計得到包括測站坐標(biāo)、衛(wèi)星軌道、天頂對流層延遲的單日松弛解.接著用GLOBK將區(qū)域松弛解與SOPAC(Scripps Orbital and Permanent Array Center)解算的全球IGS(International GNSS Service)站的單日松弛解合并，得到一個包含全球IGS站和本文GNSS站的單日松弛解.最后在全球范圍內(nèi)選擇～60個用于實現(xiàn)參考框架轉(zhuǎn)換的參考站，以全球單日松弛解作為準(zhǔn)觀測值，利用GLOBK通過7參數(shù)的相似變換得到ITRF2014(Altamimi et al., 2016)下的單日坐標(biāo)解.

圖2 瑪多地震震區(qū)周圍GNSS點位分布圖中紅色實線表示發(fā)震斷層的地表破裂軌跡，黑色實線表示震中附近主要活動斷層.Fig.2 The distribution of GNSS sites around the Madoi earthquakeThe red solid line denotes surface track of the rupture of the Madoi earthquake and the black solid lines denote the main active faults around the epicenter in the figure.

貝葉斯后驗概率密度統(tǒng)計方法已被廣泛應(yīng)用于震間斷層運動參數(shù)、同震斷層滑動分布、震后形變估計等研究(Sun et al., 2013；周新, 2017；Ingleby et al., 2020).本文引入一個基于MCMC(Markov Chain Monte Carlo)采樣的貝葉斯方法來估計同震位移及其誤差.該算法可以充分顧及震間速率估計的不確定性，尤其適用于處理流動資料.

在地殼運動研究中，GNSS坐標(biāo)時間序列可表示為(Nikolaidis et al., 2001; Bock and Melgar, 2016)：

y(t)=a+bt+csin(2πt+φ1)+dsin(4πt+φ2)

(1)

式中a為初始形變量，b為長期線性速度，兩個正弦函數(shù)表示周年和半周年季節(jié)性變形，g為在T時刻的同震形變量，H為階躍函數(shù)(Heaviside function)，ε為觀測誤差.對于流動站點，由于采樣數(shù)據(jù)不足，不考慮周期項，同時忽略早期震后效應(yīng)；而對連續(xù)觀測站而言，由于計算同震形變只用到了地震前后幾天的數(shù)據(jù)，亦可忽略震后效應(yīng)及周期形變.因此，待估參數(shù)僅剩初始形變量a、長期速度b和同震形變量g.

依照貝葉斯理論，模型參數(shù)的后驗概率密度分布可表示為(Elliott et al., 2016):

(2)

式中d為GNSS觀測得到震前、震后兩期觀測結(jié)果，m為模型參數(shù)，P(m)為模型參數(shù)的先驗概率密度分布，Cd為數(shù)據(jù)的協(xié)方差矩陣，G(m)為將采樣的模型參數(shù)代入公式(1)得到的預(yù)測值.

首先利用沈正康等(Shen et al., 1996)提出的應(yīng)變率模擬方法(strain rate modeling method)進(jìn)行空間插值，得到長期速率初始值bmod.插值得到的長期速率必然存在誤差，我們利用均一分布(uniform distributions)給出其先驗概率密度分布，取值范圍為[bmod-2,bmod+2] mm·a-1.同震形變量g的取值范圍定為 [-2,2]m，同樣采用均一的先驗概率密度分布.采用Metropolis-Hastings采樣算法(Mosegaard and Tarantola, 1995)生成Markov鏈，包含6000個采樣，在計算時先燃燒(burn-in)前3000個采樣，以保證后一段進(jìn)入收斂區(qū)域, 最終求得估計模型參數(shù)并給出其誤差.由此得到瑪多地震106個GNSS測站的同震位移場，圖3給出水平方向的同震位移場(數(shù)值結(jié)果見電子附表1)，同震位移平均誤差為4 mm.

1.3 高頻GNSS動態(tài)形變數(shù)據(jù)處理

震時1 Hz高頻GNSS數(shù)據(jù)處理采用武漢大學(xué)PRIDE課題組發(fā)布的開源GNSS處理軟件PRIDE PPP-AR(Geng et al., 2019)進(jìn)行動態(tài)解算.首先，通過對數(shù)據(jù)的預(yù)處理剔除粗差和周跳，再利用最小二乘法進(jìn)行參數(shù)估計，以獲取浮點解，最后再利用相位鐘和相位偏差產(chǎn)品來完成整周模糊度的解算，獲得固定解.

本文獲取了陸態(tài)網(wǎng)絡(luò)連續(xù)站和青海CORS站部分站點的高頻波形數(shù)據(jù)，圖4給出了其中17個臺站(9個陸態(tài)網(wǎng)絡(luò)連續(xù)站、8個青海CORS站)在瑪多地震期間的三分量高頻波形記錄.陸態(tài)網(wǎng)絡(luò)連續(xù)站中距離斷層最近的青?，敹嗾?QHMD，距離斷層34 km)形成了30 cm的永久形變(最大振幅南北向達(dá)到約16 cm，東西向達(dá)到約28 cm)，而青海CORS站中有多個臺站均形成了永久形變，包括靠近斷層的JDUO站(距離斷層41 km)，該站也形成了近30 cm的永久形變.另外距離斷層不到80 km的KANQ站和HSHX站均形成了明顯的永久形變，其他臺站也均記錄到了明顯的動態(tài)波形.

2 同震變形特征分析

在圖3中我們利用三種比例尺的箭頭，給出了不同震中距范圍內(nèi)的觀測站點的同震形變結(jié)果.從圖3中可以看出，遠(yuǎn)場連續(xù)GNSS觀測站表現(xiàn)出明顯的四象限分布特征(圖3紅色箭頭所示)，斷層西北、東南兩側(cè)測站同震位移均向遠(yuǎn)離斷層的方向運動，而東北、西南兩側(cè)測站的同震位移則指向斷層方向.而在近場范圍內(nèi)(圖3藍(lán)色箭頭所示)，GNSS觀測到的最大水平同震形變?yōu)?.2 m，位于野馬灘附近(大野馬嶺站，編號：4454)，距離斷層約100 m.另一個同震水平形變大于1.1 m的站點位于昌麻河鄉(xiāng)附近(昌麻河鄉(xiāng)站，編號：4499)，距離斷層近5 km.

觀測結(jié)果顯示，沿地表破裂兩側(cè)區(qū)域50 km范圍內(nèi)有20 cm以上的地表永久變形，震中范圍300 km仍可觀測到同震形變，影響范圍最遠(yuǎn)波及至東南方向的四川甘孜爐霍一帶.在地表破裂跡線西北角距離斷層20 km的鄂陵湖沿岸，有多個GNSS測站顯示出有20～30 cm的NW向運動，而更靠近震中的瑪多縣城附近的站點則表現(xiàn)為NWW向運動.地表破裂跡線東南角距斷層25 km處的陸態(tài)區(qū)域站J406有幅度約30 cm的SE向運動.

圖3 瑪多地震水平方向同震位移場不同顏色、不同比例尺箭頭表示同震水平位移，紅色實線表示發(fā)震斷層的地表破裂軌跡，淺黑色實線表示震中附近主要活動斷層.兩個灰色虛線框標(biāo)示了兩條跨斷層的GNSS剖面：P1，跨江錯斷裂剖面；P2，跨東昆侖斷裂瑪沁—瑪曲段的剖面.Fig.3 Horizontal coseismic displacements of the Madoi earthquakeThe arrows with different color and different scale denote the horizontal coseismic displacements. Red solid line denotes surface track of the rupture and shallow black solid lines denote main active faults around the epicenter. Two gray dashed boxes denote two GNSS profiles, which cross the Jiangcuo fault (P1) and the Maqên-Maqu segment of East Kunlun fault (P2), respectively.

圖4 連續(xù)GNSS站記錄到的1 Hz高頻波形數(shù)據(jù)上行給出了8個青海CORS站的記錄，下行給出了9個陸態(tài)網(wǎng)絡(luò)工程的連續(xù)站的記錄.左圖右側(cè)給出的公里數(shù)表示臺站的震中距.灰色虛線標(biāo)示了發(fā)震時刻.Fig.4 1-Hz high-frequency waveforms recorded by the continuous GNSS sitesThe top row plots the waveforms recorded by 8 Qinghai CORS sites and the bottom row plots the waveforms recorded by 9 CMONOC continuous sites.The kilometers at the right sides of two left subfigures denotes epicentral distance of every sites.The gray dashed line in every subfigure denotes the original time of the event.

圖5 兩條GNSS剖面水平同震形變的投影結(jié)果(a)和(b)為跨發(fā)震斷層江錯斷裂的GNSS剖面結(jié)果，(c)和(d)為跨東昆侖斷裂瑪沁—瑪曲段的GNSS剖面結(jié)果.上行為平行于斷層位移，下行為垂直于斷層位移.F1: 東昆侖斷裂，F(xiàn)2：久治斷裂，F(xiàn)3: 瑪多—甘德斷裂，F(xiàn)4: 達(dá)日斷裂.Fig.5 The projection of horizontal coseismic deformation of two GNSS profiles(a) and (b) are the projection of the GNSS profile which cross the seismogenic fault—the Jiangcuo fault; (c) and (d) are the projection of the GNSS profile which cross the Maqên-Maqu segment of East Kunlun fault.The top row is the projection parallel to the strike of faults and the bottom row is the projection perpendicular to the strike of faults.F1: East Kunlun fault, F2: Jiuzhi fault, F3: Madoi-Gande fault, F4: Dari fault.

近場GNSS觀測結(jié)果較為清晰地勾勒出了發(fā)震斷層兩側(cè)的運動模式.斷層北側(cè)臺站位移方向由東到西，從西南向逐漸轉(zhuǎn)為北西西方向、北西方向；南側(cè)臺站位移則呈現(xiàn)出一種順時針變化的變形形態(tài).近場形變在西北—東南方向上表現(xiàn)出明顯拉張變形特征，在東北—西南方向表現(xiàn)出擠壓變形特征，近場形變模式與遠(yuǎn)場的四象限變形特征相符.發(fā)震斷層兩側(cè)的近場同震形變表現(xiàn)為明顯的左旋走滑特征，與震源機制一致.

利用近場GNSS測站，本文構(gòu)建了兩條GNSS剖面(圖3中的兩個灰色虛線框)，分別為跨發(fā)震斷層的剖面P1(長寬分別為400 km和180 km)和一條跨東昆侖斷裂瑪沁—瑪曲段的剖面P2(長寬分別為300 km和80 km).圖5給出了這兩條GNSS剖面水平同震形變在平行于斷層和垂直于斷層方向上的投影結(jié)果.

跨江錯斷裂剖面P1平行于斷層方向上的位移呈現(xiàn)出明顯的雙曲線特征(如圖5a所示)，近斷層同震形變較大，平行于斷層的同震形變隨斷層距離的增加而逐漸衰減.GNSS觀測到的同震形變場直觀地展示了發(fā)震斷層的左旋走滑性質(zhì).垂直于斷層方向的同震位移相對較小(如圖5b所示).

在復(fù)測的GNSS站點中，距斷層最近的兩個站點為位于野馬灘附近的4454(距離斷層～100 m)和位于黃河鄉(xiāng)附近的2836(距離斷層～1 km)，分別位于同震破裂的北側(cè)和南側(cè).上述兩個站點觀測到的平行于斷層的同震位移分別為～1.2 m和～0.6 m，因此點位間位移差為～1.8 m.這一量值與地表調(diào)查(潘家偉等，2021)獲得的～2 m的地表破裂量值相當(dāng).

東昆侖斷裂托索湖段距離瑪多地震震中僅～70 km.跨東昆侖斷裂的同震GNSS位移剖面顯示(如圖5a、b藍(lán)色粗實線指示了東昆侖斷裂F1的位置)，東昆侖斷裂兩側(cè)站點同震響應(yīng)顯著.這些站均向西運動且東昆侖斷裂南側(cè)的站點形變量較大，這與東昆侖斷裂左旋走滑的震間形變特征相反，因此瑪多地震的發(fā)生對該斷裂段的震間應(yīng)力積累有緩解作用.

跨東昆侖斷裂瑪沁—瑪曲段的剖面P2顯示，斷層兩側(cè)水平運動差異同樣較為明顯(見圖5c、d)，南側(cè)站點向遠(yuǎn)離地震震中的東南方向變形，最大值可達(dá)6 cm，斷層北側(cè)點位朝西南方向運動，最大值不到2 cm.東昆侖斷裂瑪沁—瑪曲段兩側(cè)的地殼差異運動將造成該斷層的庫侖應(yīng)力加載，未來地震風(fēng)險性將進(jìn)一步增強.

3 與InSAR圖像的對比

基于合成孔徑雷達(dá)差分干涉測量(D-InSAR)技術(shù)，利用歐空局(European Space Agency, ESA)C波段升、降軌哨兵1號衛(wèi)星數(shù)據(jù)可獲取青海瑪多地震的InSAR同震形變場.InSAR資料顯示瑪多地震產(chǎn)生的同震形變影響空間范圍廣，形變場的長軸呈NWW向，升、降軌觀測到的形變量符號相反(見圖6).總體而言，瑪多地震所在區(qū)域干涉條件很好，InSAR觀測形成的干涉條紋連續(xù)，所獲的同震形變清晰，為與GNSS觀測相互驗證提供了良好的資料.然而，由于受到黃河水系、湖泊、沼澤等地貌影響，在斷層近處InSAR影像存在失相干區(qū)域.因此，近場的GNSS觀測為失相關(guān)區(qū)域的形變特征提供了關(guān)鍵約束.

本文計算了4條垂直于斷層InSAR剖線(圖6中綠色實線)的同震觀測結(jié)果，連同在圖3中剖面P1的GNSS位移觀測一并展示在圖7中.這4條剖線分別處于四段地表破裂的中間位置.以圖6中每條InSAR剖線為中軸線，東西各展開25 km的距離，形成4條GNSS剖面.結(jié)合InSAR升、降軌方向矢量，將各GNSS剖面中的觀測站點的三維同震位移轉(zhuǎn)化為LOS(Line-of-Sight)向形變，將其與對應(yīng)的InSAR剖線升、降軌LOS向位移進(jìn)行比較(見圖8).利用升、降軌InSAR數(shù)據(jù)得到的最大LOS向形變量約為0.8 m，與由GNSS三維位移導(dǎo)出的LOS向形變量大體一致(見圖7、圖8).

為了擬合瑪多地震同震InSAR觀測的升、降軌LOS向形變和GNSS觀測的位移，本文采用一種半無限空間淺表位錯模型(shallow dislocation model)來表征地表位移與同震滑動的關(guān)系.在同震階段，地表位移與發(fā)震斷層的寬度(傾角方向)和平均滑動量可以簡單的表示為(Segall, 2010)：

(3)

式中u3為平行斷層走向的位移，s為斷層平均滑動量，d1為斷層寬度，x1為測站離開斷層的距離，其中s和d1為待求模型參數(shù).假設(shè)走滑斷層同震破裂中平行于斷層走向的位移遠(yuǎn)大于垂向位移，即忽略垂向位移，則可將平行斷層走向的同震位移u3轉(zhuǎn)換成對應(yīng)的InSAR觀測LOS向位移.

圖6 歐空局哨兵1號衛(wèi)星給出的瑪多地震InSAR同震形變場(a)為升軌觀測結(jié)果，(b)為降軌觀測結(jié)果.兩者形變量符號正好相反.Fig.6 InSAR coseismic deformation field caused by the Madoi earthquake from Sentinal-1A satellite of the ESA(a) is LOS displacements along the ascending orbit and (b) is LOS displacements along the descending orbit.The sign of the two displacement values are opposite.

圖7 4條InSAR跨斷層剖線的LOS向位移觀測上圖為升軌觀測，下圖為降軌觀測.圖中同時給出了前文跨江錯斷裂的GNSS剖面的位移結(jié)果，為方便與InSAR觀測比較，GNSS位移已經(jīng)通過方向矢量轉(zhuǎn)化為LOS方向的投影.灰色粗實線標(biāo)示了發(fā)震斷層的位置.Fig.7 The LOS displacements of 4 InSAR profiles crossing the seismogenic faultThe top figure is LOS displacements obtained from the ascending orbit and the bottom figure is LOS displacements obtained from the descending orbit.The coseismic deformation of the GNSS profile crossing the Jiangcuo fault is also shown in two subfigures and the displacements had been converted to the projection along LOS with the direction vector due to the comparison with InSAR observation.The thick gray solid line in every subfigure denotes the seismogenic fault.

圖8 4條InSAR剖線與對應(yīng)的GNSS剖面結(jié)果的比較圓點為InSAR剖線結(jié)果，三角形為GNSS剖面的結(jié)果(已利用升、降軌方向矢量轉(zhuǎn)化為對應(yīng)的LOS向位移).上行為升軌結(jié)果，下行為降軌結(jié)果.紅色實線為擬合結(jié)果，灰色豎直粗實線標(biāo)示了發(fā)震斷層的位置.Fig.8 The comparison between 4 InSAR profiles and the corresponding GNSS displacementsThe solid dots denote the LOS displacements of the InSAR profiles and the triangles denote the corresponding GNSS displacements in the LOS direction converted by the direction vectors of the ascending and descending orbit.The top row is the results along the ascending orbit and the bottom row is the results along the descending orbit.The red solid line is the fitting result and the vertical gray solid line denotes the seismogenic fault in every subfigure.

將圖8中各條InSAR剖線的升、降軌LOS向位移，以及對應(yīng)的GNSS剖面位移的LOS向投影作為觀測數(shù)據(jù)，估計模型參數(shù).考慮到模型參數(shù)與觀測數(shù)據(jù)之間的非線性關(guān)系，采用格網(wǎng)搜索方法來確定s和d1的最佳取值.通過調(diào)節(jié)兩個參數(shù)的取值使得觀測數(shù)據(jù)與模擬結(jié)果之間的殘差最小，以達(dá)到最佳的數(shù)據(jù)擬合效果.圖8每列圖形的行首給出了各條剖線的模型參數(shù)最佳取值，擬合結(jié)果由每個子圖中的紅色實線表示.參數(shù)估計結(jié)果顯示，最大平均滑動量出現(xiàn)在冬草阿龍湖段，達(dá)到了3.765 m，其斷層破裂寬度為6.8 km；而最大的斷層破裂寬度則出現(xiàn)在昌麻河段，達(dá)到13.5 km，平均滑動量為1.565 m.

擬合曲線總體上較好地擬合了各條剖線的LOS向位移，尤其在斷層近場附近，多數(shù)剖線上的擬合結(jié)果與實際觀測的LOS向位移均有較好的符合.而昌麻河段InSAR剖線P4在斷層北側(cè)近場附近的升、降軌LOS向位移均與擬合曲線有相對較大的偏差.昌麻河段真實斷層分布較為復(fù)雜，地表破裂野外考察(潘家偉等，2021)和InSAR影像的結(jié)果，均證實了在地表破裂昌麻河段存在分支斷層.在該段的擬合過程中，未考慮分支斷層，因此可能造成了擬合曲線在該段與實際觀測出現(xiàn)相對較大的偏差.此外，走滑斷裂端部的垂向位移相對較大，與假設(shè)不符，這是出現(xiàn)較大偏差的另一因素.

4 同震滑動分布及其應(yīng)力擾動

利用本文獲取的遠(yuǎn)近場GNSS同震形變場結(jié)果，我們構(gòu)建了簡單的斷層幾何模型，采用約束最小二乘算法反演了斷層滑動分布.InSAR成像及野外地質(zhì)考察(李智敏等，2021；潘佳偉等，2021)顯示瑪多地震地表破裂呈現(xiàn)明顯的分段特征，斷層走向總體為285°，但每段存在差異.如前文所述，在東側(cè)昌麻河段還存在分支斷層破裂.余震精定位結(jié)果(王未來等，2021)也顯示出較為明顯的分段性，但兩者在地表并不重合，可能說明斷層傾角的變化，也可能與余震定位的偏差有關(guān).為此我們依據(jù)InSAR影像識別的破裂跡線構(gòu)建斷層面，并假設(shè)斷層傾角為90°，寬度為30 km.斷層傾角的微弱變化并不會改變采用GNSS資料反演的滑動分布基本特征.我們將帶有分支斷層的曲面離散化為1574個三角形子塊，三角形的平均邊長由近地表2 km隨深度逐漸增加至5.5 km.

基于GNSS同震形變，本文反演得到的瑪多地震破裂模型見圖9.結(jié)果顯示，同震破裂具有明顯的分段性，并且同震滑動出露地表，與野外地表破裂考察和余震分布吻合.最大的破裂區(qū)域處于冬草阿龍湖段，黃河鄉(xiāng)段和鄂陵湖南段均有凹凸體存在，滑動量略小.反演結(jié)果顯示地震主體破裂位于斷層面0～10 km的淺部區(qū)域，主破裂面最大滑動量達(dá)到了4.6 m，地震矩1.63×1020N·m，矩震級為MW7.4.昌麻河段分支斷層最大滑動量可達(dá)2.3 m.

圖9所示的瑪多地震破裂模型可以很好地擬合GNSS觀測的同震形變場(見圖10)，擬合殘差為3.9 cm.通過單一手段給出的斷層滑動分布對原始觀測量的恢復(fù)程度較好，說明斷層幾何參數(shù)的選取相對可靠.近場的GNSS同震形變結(jié)果可以有效填補InSAR觀測在極震區(qū)由于失相干造成的觀測空白，對斷層滑動破裂提供有效約束.

華俊等(2021)基于InSAR觀測反演得到的瑪多地震斷層滑動分布結(jié)果顯示，斷層最大滑動量約為6 m，矩震級為MW7.45.與之相比，本文基于GNSS觀測反演得到的最大滑動量和矩震級都略小，但同震破裂的分段性特征一致，且均破裂到地表，主體破裂的深度范圍也較一致，均集中于斷層面0～10 km的淺部區(qū)域.

根據(jù)圖9所示的滑動分布模型，假設(shè)摩擦系數(shù)為0.4，計算了斷層面上的靜態(tài)庫侖應(yīng)力變化(見圖11).結(jié)果顯示庫侖應(yīng)力增加區(qū)域主要集中在同震滑移下傾方向，在黃河鄉(xiāng)段中部和冬草阿龍湖的西段存在淺部應(yīng)力增強的區(qū)域.靜態(tài)庫侖應(yīng)力增加區(qū)域與余震分布具有一致性，說明余震主要是由靜態(tài)庫侖應(yīng)力加載而觸發(fā)的.

5 討論

觀測到的同震形變顯示黃河鄉(xiāng)段和昌麻河段同震滑動顯著，這一結(jié)論與InSAR結(jié)果相符(華俊等，2021).但地表調(diào)查結(jié)果顯示鄂陵湖南段的地表破裂最為顯著，黃河鄉(xiāng)段和昌麻河段地表破裂相對較小(潘家偉等，2021).我們推測造成上述差異的原因可能有兩方面.其一是鄂陵湖南段水系發(fā)達(dá)，地表土層含水量較高，在同震破裂中可能發(fā)生較強的非彈性形變，導(dǎo)致相對顯著的地表破裂.其二是鄂陵湖南段的GNSS站點較少，且缺乏近斷層站點約束，導(dǎo)致我們的同震模型對該區(qū)域的同震滑動分辨能力不足.

圖9 基于GNSS同震形變反演得到的瑪多地震破裂模型Fig.9 The rupture model of Madoi earthquake inverted from the GNSS observed coseismic displacements

圖10 (a) GNSS觀測的瑪多地震水平方向同震位移與破裂模型模擬結(jié)果的比較； (b) 模擬值和觀測值的擬合殘差Fig.10 (a) The comparison of GNSS observed horizontal coseismic displacements of Madoi earthquake with simulated results by the rupture model; (b) The residuals between the observation and the simulations

圖11 瑪多地震在發(fā)震斷層上產(chǎn)生的庫侖應(yīng)力變化圖中灰色圓圈標(biāo)示了精定位得到的余震位置(王未來等, 2021).Fig.11 The Coulomb stress changes on the modeled fault imposed by the Madoi earthquake Open gray circles represent relocated aftershocks (Wang et al., 2021).

根據(jù)震間GNSS觀測約束的滑動速率，朱亞戈等(2021)估計得到相同破裂規(guī)模的地震復(fù)發(fā)周期(T)為1100～5500年.震后大地測量約束的青藏高原不同區(qū)域的下地殼穩(wěn)態(tài)黏滯系數(shù)η≈1×1019Pa·s(Huang et al., 2014; Zhao et al., 2017; Diao et al., 2018; Liu et al., 2019)，由此計算的震后松弛時間(Tr=2μ/η,μ為剪切模量)遠(yuǎn)小于地震復(fù)發(fā)周期.在此條件下，震間的滑動速率強烈依賴于上次地震的離逝時間(te q), 震后早期滑動速率大于整個地震周期平均速率，而在晚期則低于平均速率.由此推測，朱亞戈等(2021)估計的結(jié)果一定程度上低估了滑動速率，高估了地震復(fù)發(fā)周期.為此，可以結(jié)合瑪多地震的震后資料開展更深入的研究.

巴顏喀拉塊體內(nèi)部的GNSS觀測站點分布較為稀疏，難以約束塊體內(nèi)部斷裂的震間閉鎖深度(朱亞戈等，2021).2021年瑪多地震的同震破裂集中于斷層面0～10 km淺部，本文第3節(jié)給出的4條剖線的平均斷層寬度為9.6 km，亦小于10 km，表明本次地震的主要破裂并未向深部擴展，因此孕育本次地震的昆侖山口—江錯斷裂的震間閉鎖深度很可能在10 km左右.前人利用震間大地測量觀測反演的青藏高原活動斷層的平均閉鎖深度多為15～20 km(如Wang et al., 2011; 李煜航等，2014).昆侖山口—江錯斷裂“低于預(yù)期”的閉鎖深度表明，震間大地測量觀測獲取的地表形變與常用的彈性位錯模型(Savage and Burford，1973)可能容易高估活動斷層(特別是滑動速率較低的塊體內(nèi)部斷層)的閉鎖深度.

長期以來，我們對類似的板塊內(nèi)次級分支斷裂的關(guān)注度一向偏低，如龍日壩斷裂、發(fā)生2017年九寨溝地震的虎牙斷裂都屬于這種類型.瑪多地震的發(fā)生引起了我們對次級斷裂活動性及其孕震能力的重新審視和思考.

6 結(jié)論

本文使用GNSS觀測技術(shù)獲取了瑪多地震較大范圍的同震形變場，揭示了瑪多地震在巴顏喀拉地塊內(nèi)部的形變特征及對周邊斷裂造成的潛在影響.

(1)同震形變結(jié)果顯示，瑪多地震為左旋走滑型地震，同震變形場呈明顯的四象限分布，觀測到的最大同震位移達(dá)到1.2 m.GNSS和InSAR數(shù)據(jù)較為一致地展示了瑪多地震同震形變的表現(xiàn)方式和變化幅度，近場GNSS為InSAR近斷層失相干區(qū)域提供了必要的補充和檢核，直觀詮釋了昆侖山口—江錯斷裂左旋走滑的發(fā)震機制.

(2)采用GNSS數(shù)據(jù)反演同震滑動分布的結(jié)果顯示，發(fā)震斷層的滑動破裂存在多個凹凸體，破裂分段特征明顯且出露地表，與野外地表破裂考察和余震精定位結(jié)果一致，主體破裂位于斷層面0～10 km的淺部區(qū)域，最大滑動量達(dá)到4.6 m，地震矩1.63×1020N·m，矩震級為MW7.4.同震破裂引起發(fā)震斷層面上的應(yīng)力擾動，余震主要集中在同震滑移下傾方向庫侖應(yīng)力增加的區(qū)域.

(3)從瑪多地震的同震形變場空間特征來看，形變主要集中在巴顏喀拉塊體內(nèi)部，瑪多地震的發(fā)生對東昆侖斷裂托索湖段的震間應(yīng)力積累有緩解作用；而東端瑪沁—瑪曲段的同震響應(yīng)形成了明顯的斷層南北兩側(cè)的差異運動，因此瑪多地震對東昆侖斷裂瑪沁—瑪曲段有應(yīng)力加載作用，未來地震風(fēng)險性將進(jìn)一步增強，需要加強監(jiān)測.巴顏喀拉地塊周緣強震具有明顯的跳躍性特征，長期的大區(qū)域GNSS觀測將為今后監(jiān)測地塊邊界及其內(nèi)部斷層活動性、劃定地震危險性區(qū)域提供可行的技術(shù)途徑.

致謝感謝中國地質(zhì)大學(xué)(武漢)熊熊教授提供相關(guān)點位的站點信息.