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

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

0 引言

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

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

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

圖1 (a) 青藏高原背景速度場(chǎng)； (b) 2021年瑪多地震的構(gòu)造背景圖(a)中的黑色粗實(shí)線給出了巴顏喀拉塊體的主邊界，黑色細(xì)實(shí)線給出了東昆侖斷裂附近發(fā)育的多條次級(jí)分支斷裂，紅色實(shí)線表示瑪多地震發(fā)震斷層的地表破裂軌跡，虛線矩形框指示了圖(b)的位置； 圖(b)中綠色圓點(diǎn)給出了余震序列的精定位結(jié)果(王未來(lái)等，2021).圖中的速度場(chǎng)結(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ā)展，該項(xiàng)技術(shù)已經(jīng)成為現(xiàn)代地震研究領(lǐng)域中的一個(gè)重要觀測(cè)手段(Bock and Melgar, 2016).尤其是近場(chǎng)GNSS觀測(cè)可以彌補(bǔ)InSAR影像在近場(chǎng)失相干區(qū)域的觀測(cè)不足，直觀詮釋發(fā)震機(jī)制，同時(shí)為地震動(dòng)態(tài)破裂過(guò)程和靜態(tài)破裂分布模型提供重要的約束.李志才等(2021)展示了破裂帶附近9個(gè)GNSS基準(zhǔn)站的高頻(1 Hz采樣)動(dòng)態(tài)形變波形，及12個(gè)GNSS連續(xù)站的同震形變場(chǎng)，觀測(cè)到的最大水平同震位移0.6 m，距離斷層20 km處.該研究雖然初步揭示了本次地震的形變特征，但由于點(diǎn)位稀疏，對(duì)斷層滑動(dòng)分布的約束不足，破裂模型稍顯粗糙.

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

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

1.1 觀測(cè)資料

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

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

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

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

圖2 瑪多地震震區(qū)周?chē)鶪NSS點(diǎn)位分布圖中紅色實(shí)線表示發(fā)震斷層的地表破裂軌跡，黑色實(shí)線表示震中附近主要活動(dòng)斷層.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.

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

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

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

(1)

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

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

(2)

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

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

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

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

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

2 同震變形特征分析

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

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

圖3 瑪多地震水平方向同震位移場(chǎng)不同顏色、不同比例尺箭頭表示同震水平位移，紅色實(shí)線表示發(fā)震斷層的地表破裂軌跡，淺黑色實(shí)線表示震中附近主要活動(dòng)斷層.兩個(gè)灰色虛線框標(biāo)示了兩條跨斷層的GNSS剖面：P1，跨江錯(cuò)斷裂剖面；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個(gè)青海CORS站的記錄，下行給出了9個(gè)陸態(tài)網(wǎng)絡(luò)工程的連續(xù)站的記錄.左圖右側(cè)給出的公里數(shù)表示臺(tái)站的震中距.灰色虛線標(biāo)示了發(fā)震時(shí)刻.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ā)震斷層江錯(cuò)斷裂的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.

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

利用近場(chǎng)GNSS測(cè)站，本文構(gòu)建了兩條GNSS剖面(圖3中的兩個(gè)灰色虛線框)，分別為跨發(fā)震斷層的剖面P1(長(zhǎng)寬分別為400 km和180 km)和一條跨東昆侖斷裂瑪沁—瑪曲段的剖面P2(長(zhǎng)寬分別為300 km和80 km).圖5給出了這兩條GNSS剖面水平同震形變?cè)谄叫杏跀鄬雍痛怪庇跀鄬臃较蛏系耐队敖Y(jié)果.

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

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

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

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

3 與InSAR圖像的對(duì)比

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

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

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

(3)

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

圖6 歐空局哨兵1號(hào)衛(wèi)星給出的瑪多地震InSAR同震形變場(chǎng)(a)為升軌觀測(cè)結(jié)果，(b)為降軌觀測(cè)結(jié)果.兩者形變量符號(hào)正好相反.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向位移觀測(cè)上圖為升軌觀測(cè)，下圖為降軌觀測(cè).圖中同時(shí)給出了前文跨江錯(cuò)斷裂的GNSS剖面的位移結(jié)果，為方便與InSAR觀測(cè)比較，GNSS位移已經(jīng)通過(guò)方向矢量轉(zhuǎn)化為L(zhǎng)OS方向的投影.灰色粗實(shí)線標(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剖線與對(duì)應(yīng)的GNSS剖面結(jié)果的比較圓點(diǎn)為InSAR剖線結(jié)果，三角形為GNSS剖面的結(jié)果(已利用升、降軌方向矢量轉(zhuǎn)化為對(duì)應(yīng)的LOS向位移).上行為升軌結(jié)果，下行為降軌結(jié)果.紅色實(shí)線為擬合結(jié)果，灰色豎直粗實(shí)線標(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向位移，以及對(duì)應(yīng)的GNSS剖面位移的LOS向投影作為觀測(cè)數(shù)據(jù)，估計(jì)模型參數(shù).考慮到模型參數(shù)與觀測(cè)數(shù)據(jù)之間的非線性關(guān)系，采用格網(wǎng)搜索方法來(lái)確定s和d1的最佳取值.通過(guò)調(diào)節(jié)兩個(gè)參數(shù)的取值使得觀測(cè)數(shù)據(jù)與模擬結(jié)果之間的殘差最小，以達(dá)到最佳的數(shù)據(jù)擬合效果.圖8每列圖形的行首給出了各條剖線的模型參數(shù)最佳取值，擬合結(jié)果由每個(gè)子圖中的紅色實(shí)線表示.參數(shù)估計(jì)結(jié)果顯示，最大平均滑動(dòng)量出現(xiàn)在冬草阿龍湖段，達(dá)到了3.765 m，其斷層破裂寬度為6.8 km；而最大的斷層破裂寬度則出現(xiàn)在昌麻河段，達(dá)到13.5 km，平均滑動(dòng)量為1.565 m.

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

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

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

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

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

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

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

5 討論

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

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

圖10 (a) GNSS觀測(cè)的瑪多地震水平方向同震位移與破裂模型模擬結(jié)果的比較； (b) 模擬值和觀測(cè)值的擬合殘差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)生的庫(kù)侖應(yīng)力變化圖中灰色圓圈標(biāo)示了精定位得到的余震位置(王未來(lái)等, 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觀測(cè)約束的滑動(dòng)速率，朱亞戈等(2021)估計(jì)得到相同破裂規(guī)模的地震復(fù)發(fā)周期(T)為1100～5500年.震后大地測(cè)量約束的青藏高原不同區(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)，由此計(jì)算的震后松弛時(shí)間(Tr=2μ/η,μ為剪切模量)遠(yuǎn)小于地震復(fù)發(fā)周期.在此條件下，震間的滑動(dòng)速率強(qiáng)烈依賴(lài)于上次地震的離逝時(shí)間(te q), 震后早期滑動(dòng)速率大于整個(gè)地震周期平均速率，而在晚期則低于平均速率.由此推測(cè)，朱亞戈等(2021)估計(jì)的結(jié)果一定程度上低估了滑動(dòng)速率，高估了地震復(fù)發(fā)周期.為此，可以結(jié)合瑪多地震的震后資料開(kāi)展更深入的研究.

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

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

6 結(jié)論

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

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

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

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

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