聯(lián)合InSAR和高頻GNSS位移波形反演2022年青海門源M 6.9地震同震破裂過程

呂明哲，陳克杰*，柴海山，耿江輝，張生鵬，房立華

1 南方科技大學(xué)地球與空間科學(xué)系，深圳 518055 2 武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心，武漢 430079 3 青海師范大學(xué)地理科學(xué)學(xué)院，西寧 810016 4 青海省基礎(chǔ)測繪院，西寧 810101 5 中國地震局地球物理研究所，北京 100081

0 引言

據(jù)中國地震臺網(wǎng)中心（https:∥news.ceic.ac.cn）測定，北京時間2022年1月8日凌晨1時45分左右，青海省海北州門源縣發(fā)生M6.9地震，震源位置為101.26°E，37.77°N，深度10 km.此次地震震中海拔高且人煙稀少，有少量人員受傷，對當(dāng)?shù)夭糠指哞F線路和橋梁等造成了嚴(yán)重?fù)p壞（Yang et al.,2022）.野外現(xiàn)場考察結(jié)果顯示（潘家偉等,2022），門源地震形成了總長約27 km同震地表破裂帶，錯動方式為左旋走滑，與全球質(zhì)心矩張量（Global Centroid Moment Tensor,GCMT）給出的震源機(jī)制解一致（圖1c），發(fā)震斷層初步確定為冷龍嶺斷裂西段和托萊山斷裂東段（李智敏等,2022）.

冷龍嶺斷裂和托萊山斷裂位于海原斷裂帶西段.海原斷裂帶是一條長約1000 km的左旋走滑斷裂帶，具有部分由南向北的逆沖分量，是柴達(dá)木—祁連地塊北部邊界（圖1b），在調(diào)節(jié)印度板塊和歐亞板塊碰撞持續(xù)變形方面發(fā)揮著重要作用（袁道陽等,2004;Jolivet et al.,2012）.海原斷裂帶在20世紀(jì)發(fā)生了兩次大地震：1920年海原M8.7地震（Deng et al.,1986）和1927年古浪M8.0地震（Gaudemer et al.,1995）.在上述兩次強(qiáng)震之間存在一條260 km長的未破裂段，被稱為天祝地震空區(qū)（圖1c）.該段斷層在過去幾個世紀(jì)沒有記錄到7級以上地震，具有孕育大型地震的潛力（Gaudemer et al.,1995）.

2022年門源地震使天祝地震空區(qū)最西端冷龍嶺段和托萊山斷裂東端發(fā)生破裂（圖1c），其中冷龍嶺斷裂周邊區(qū)域在過去100年間發(fā)生過兩次較大地震：1986年門源MW6.0地震和2016年門源MW5.9地震（徐紀(jì)人等,1986;陳兵等,2003;Zhang et al.,2020），表明該區(qū)域有較強(qiáng)地震活動性.深入研究本次門源地震的破裂過程對于理解天祝地震空區(qū)西端冷龍嶺段發(fā)震特征和評估該區(qū)域未來地震危險性具有重要意義.

門源地震發(fā)生后，國內(nèi)學(xué)者利用InSAR和GNSS觀測資料研究了同震地表靜態(tài)形變和震源機(jī)制（李振洪等,2022;Yang et al.,2022;顏丙囤等,2022;馮萬鵬等,2022），但尚未對地震破裂時空演化過程進(jìn)行分析.本文首先采用差分干涉合成孔徑雷達(dá)（Differential InSAR,DInSAR）和像素偏移追蹤（Pixel Offset Tracking,POT）處理Sentinel-1數(shù)據(jù)獲取了同震地表形變場和地表破裂跡線，然后利用單歷元精密單點定位（Precise Point Positioning,PPP）獲取了震中附近16個GNSS測站（圖1）的1 Hz同震位移波形.最后，聯(lián)合上述觀測資料反演了本次地震同震滑動分布和時空過程，并對該區(qū)域靜態(tài)庫侖應(yīng)力擾動進(jìn)行分析.

1 同震形變獲取

1.1 InSAR同震形變

門源地震發(fā)生在海拔3500～4200 m的祁連山區(qū)，該區(qū)域幾乎沒有植被覆蓋（潘家偉等,2022），這使得歐空局（European Space Agency,ESA）C波段（波長5.54 cm）SAR衛(wèi)星Sentinel-1在該地區(qū)干涉效果較好.本研究選取覆蓋研究區(qū)域的4幅基于漸進(jìn)式掃描地形觀測（Terrain Observation with Progressive Scans,TOPS）升降軌Sentinel-1A干涉寬幅模式（Interferometric Wide,IW）SAR影像（T128A和T33D），利用開源軟件ISCE2進(jìn)行InSAR數(shù)據(jù)處理（Rosen et al.,2012）.SAR數(shù)據(jù)基本信息如表1所示，升降軌干涉影像標(biāo)準(zhǔn)差是從遠(yuǎn)離形變區(qū)域的像素估計得到（溫?fù)P茂等,2014）.

圖1 2022年門源M 6.9地震測站分布及構(gòu)造背景（a） GNSS臺站分布（紅色三角形）和InSAR影像覆蓋區(qū)域（藍(lán)色矩形框）.藍(lán)色五角星為本次地震震中位置，黃色圓點為余震分布（Fan et al.,2022），黑色震源球表示GCMT提供的1976年到2022年MW>5.0地震震源機(jī)制解；（b） 研究區(qū)構(gòu)造背景，紅色矩形框為圖（a）所處區(qū)域；（c） 天祝地震空區(qū)在海原斷裂帶上所處位置（紅色斷層），綠色斷層表示1920年海原地震破裂帶，紅色、黑色震源球分別表示2022年門源M 6.9地震（GCMT提供）和1920年海原M 8.2地震（Ou et al.,2020）和1927年古浪M 8.0地震（Hou et al.,1999）震源機(jī)制解.Fig.1 Station distribution and tectonic background of the 2022 Menyuan M 6.9 earthquake（a） GNSS station distribution （red triangles） and InSAR coverage （blue rectangles）.The blue star denotes the epicenter,the yellow dots represent the distribution of aftershocks （Fan et al.,2022）,and the black beachball represents the focal mechanism solution of historical earthquakes with MW>5.0 provided by GCMT from 1976 to 2022;（b） Tectonic background of the study area,and red rectangle outlines area of figure （a）;（c） The location of the Tianzhu earthquake gap on the Haiyuan fault zone （red faults）,and the green faults are the rupture fault of the 1920 Haiyuan earthquake.Red and black beachballs show the focal mechanism solutions of the 2022 Menyuan M 6.9 earthquake provided by GCMT,the 1920 Haiyuan M 8.2 earthquake （Ou et al.,2020） and the 1927 Gulang M 8.0 earthquake （Hou et al.,1999）,respectively.

表1 InSAR觀測信息Table 1 InSAR data information

本文首先采用兩軌法對升降軌SAR數(shù)據(jù)進(jìn)行差分干涉處理，并使用ESA提供的精密軌道數(shù)據(jù)（https:∥s1qc.asf.alaska.edu）來調(diào)整SAR衛(wèi)星影像軌道參數(shù)，以實現(xiàn)SAR主輔影像精確配準(zhǔn)，然后采用美國宇航局發(fā)布的30 m分辨率SRTM DEM數(shù)據(jù)（Farr et al.,2007） 去除地形相位.為了抑制相位噪聲和提高信噪比，使用功率譜濾波器對干涉圖進(jìn)行濾波（Goldstein and Werner,1998），并進(jìn)行8×2（距離向×方位向）多視處理（Vaka et al.,2020）.使用SNAPHU軟件對濾波后的干涉圖進(jìn)行相位解纏（Chen and Zebker,2000）.最后，使用WGS84坐標(biāo)系對解纏后的干涉圖進(jìn)行地理編碼以獲取門源地震同震視線向（Line of Sight,LOS）地表形變場（圖2）.

SAR影像包含相干（干涉相位）和非相干信息（幅度特征）（Jiang et al.,2017）.在地震破裂近場區(qū)域，由于地表遭受劇烈破壞導(dǎo)致該區(qū)域位移梯度較高，可使SAR衛(wèi)星干涉相位失相干（Schmidt and Bürgmann,2006）.DInSAR在這些失相干區(qū)域無法進(jìn)行可靠相位解纏，因此LOS向形變場在該區(qū)域形變信息丟失.而近場形變信息在同震研究中能提取精確地表破裂跡線，從而有效約束斷層長度和走向等幾何信息.為了獲取本次地震地表破裂跡線，使用POT計算距離向像素偏移.

圖2 2022年門源M 6.9地震InSAR同震形變場Sentinel-1升軌和降軌經(jīng)DInSAR處理得到的LOS向形變場（a、c）和經(jīng)POT處理得到的距離向位移場（b、d）.藍(lán)色五角星表示震中，綠線為手動繪制的地表破裂跡線.Fig.2 Coseismic InSAR deformation field of the 2022 Menyuan M 6.9 earthquake（a），（c） are the LOS deformation field obtained by DInSAR and （b），（d） are the range displacement field obtained by POT of Sentinel-1 ascending and descending tracks,respectively.The blue star indicates the epicenter,and the green line is the manually drawn surface rupture trace.

POT利用SAR影像中非相干信息（幅度特征）來提取位移信息，將震前震后兩幅SAR影像中M（通常為幾十）個像素進(jìn)行互相關(guān)，從而獲取亞像素精度的距離向（相當(dāng)于DInSAR中的LOS向）和方位向像素偏移（Simons and Rosen,2007）.這種方法不需要進(jìn)行相位解纏（Wang et al.,2015;He et al.,2019），因此不受SAR影像中相位相干性限制，可以獲取近場地表形變.但由于POT是基于像素集合實現(xiàn)，故空間分辨率相較于DInSAR顯著降低（Simons and Rosen,2007）.其中方位向偏移信噪比更低（Jin and Fialko,2021），因此本研究只采用距離向偏移.

圖2為分別利用DInSAR和POT獲取的門源地震同震形變場，升、降軌SAR影像均完整覆蓋了震中區(qū)域.具體而言，利用干涉相位獲取的升降軌同震形變（圖2a，2c）分布范圍為60～80 km，最大視線向形變?yōu)?.61 m，發(fā)震斷層附近部分失相干區(qū)域形變被掩膜.降軌干涉形變圖在斷層南北兩側(cè)形變大致呈對稱分布，由于升降軌SAR衛(wèi)星影像存在成像幾何差異，升軌LOS向形變場南側(cè)下沉區(qū)偏西，北側(cè)抬升區(qū)偏東，并大致以震中為中心呈對稱分布.而升軌和降軌LOS向形變圖上發(fā)震斷層兩側(cè)位移相反，進(jìn)一步證實了本事件因走滑斷層而造成的水平向運(yùn)動.

利用POT獲取的升降軌距離向偏移（圖2c，2d）顯示地表破裂最大距離向偏移為0.79 m，由于其在發(fā)震斷層附近也連續(xù)覆蓋，故大于在近斷層附近失相干的最大LOS向形變（0.61 m）.同時，距離向偏移表明發(fā)震斷層走向自西向東在～101.2°E處發(fā)生顯著變化，西段斷層走向角略小于90°，東段斷層走向角增大～30°，總地表破裂跡線長度為～30 km.由于POT利用SAR影像中的幅度特征信息而不是高精度干涉相位，所以距離向偏移圖較DInSAR的LOS向形變圖存在更多斑點噪聲（Hu et al.,2012;He et al.,2019）.

InSAR原始同震形變場數(shù)據(jù)量達(dá)百萬級，并且大部分低相位梯度區(qū)域的相鄰觀測數(shù)據(jù)之間空間相關(guān)性較高，考慮到反演效率，使用四叉樹降采樣法（Jónsson et al.,2002） 對原始觀測數(shù)據(jù)進(jìn)行降采樣，該降采樣方法可以在發(fā)震斷層附近高相位梯度區(qū)域更加頻繁地采樣，從而可以保留原始數(shù)據(jù)中主要形變特征，降采樣后的升、降軌LOS向形變場數(shù)據(jù)點為1305.

1.2 GNSS數(shù)據(jù)處理

本文獲取了震中附近來自陸態(tài)網(wǎng)絡(luò)和青海省連續(xù)運(yùn)行跟蹤站（Continuously Operating Reference Stations,CORS）共計16個GNSS臺站（分布見圖1a）的1 Hz GPS/GLONASS雙系統(tǒng)2022年1月7日（UTC）當(dāng)天觀測數(shù)據(jù)，其中來自陸態(tài)網(wǎng)絡(luò)測站接收機(jī)型號為Trimble Net R9，天線型號為TRM59800.00或TRM59900.00，青海省CORS站所用接收機(jī)型號為中海達(dá)VNet8U-I，天線類型為AT-2300，上述測站均裝備扼流圈以削弱多路徑效應(yīng).本文采用成熟的GPS+GLONASS雙系統(tǒng)組合方式（An et al.,2020），GPS和GLONASS兩系統(tǒng)定權(quán)設(shè)置為2∶1（Shi et al.,2013），系統(tǒng)間偏差在計算過程中估計，通過PRIDE PPP-AR軟件完成動態(tài)PPP解算（Zumberge et al.,1997;Geng et al.,2018,2019），使用卡爾曼濾波逐歷元求解并反向濾波，每個歷元的三維坐標(biāo)和接收機(jī)鐘差作為白噪聲隨機(jī)參數(shù)估計，以全天平均坐標(biāo)做站心轉(zhuǎn)換，得到東西、南北和垂向三個時間序列，具體處理策略如表2所示.

精密產(chǎn)品選取德國地學(xué)中心（GFZ）提供的30 s精密鐘差和15 min精密軌道.對于1 Hz數(shù)據(jù)解算來說，30 s時長內(nèi)鐘差高頻抖動并不明顯，使用30 s精密鐘差即可滿足計算需要（張小紅等,2010）.產(chǎn)品日期向前后各延長一天以確保當(dāng)天開始和結(jié)束時軌道內(nèi)插精度.對每一個連續(xù)弧段，設(shè)置一個待估模糊度參數(shù).為了減少多路徑效應(yīng)影響，選擇基于高度角的觀測值定權(quán)方式，且設(shè)置截止高度角為7°.電離層使用傳統(tǒng)L1+L2無電離層組合消除，以減少待估參數(shù)總數(shù).以Saastamoinen模型計算對流層延遲，并每兩小時估計一次對流層延遲殘余，投影函數(shù)選擇Global Mapping Function（GMF）.

表2 GNSS數(shù)據(jù)處理策略Table 2 Processing strategies of GNSS data

天線相位中心偏差（Phase Center Offsets）、天線相位中心變化（Phase Center Variations）以及天線相位纏繞均考慮并改正.固體潮、海潮和極潮通過國際地球自轉(zhuǎn)和參考框架服務(wù)（International Earth Rotation and Reference System Service）2010模型精確改正.接收機(jī)硬件碼偏差等會被吸收至接收機(jī)鐘差內(nèi)，精密產(chǎn)品估計所使用的是無電離層組合，故衛(wèi)星硬件碼偏差可忽略，相位偏差會被吸收至模糊度，導(dǎo)致其失去整數(shù)特性.由于動態(tài)PPP解算需要一定收斂過程，為保證結(jié)果不受影響，本文計算了地震發(fā)生當(dāng)日全天時間序列，最后截取各測站P波理論到時前10 s和到達(dá)后50 s共計60 s同震位移波.經(jīng)過上述處理，GNSS同震波拾取精度達(dá)到了1～2 cm.由于門源地震為左旋走滑運(yùn)動，垂向位移較小，并且GNSS水平向精度高于垂向精度，這導(dǎo)致垂向波形信噪比較低，本文僅采用水平向GNSS觀測波形進(jìn)行同震反演.

2 滑動分布反演

2.1 斷層幾何參數(shù)

根據(jù)POT得到的Sentinel-1升、降軌距離向偏移（圖2c，2d）揭示的斷層地表破裂跡線，本文使用兩個矩形斷層Fault 1（靠近托萊山斷裂東段）和Fault 2（靠近冷龍嶺斷裂西段）來擬合發(fā)震斷層（圖3a）.兩個矩形斷層長度分別設(shè)定為10 km和20 km，斷層寬度均設(shè)定為16 km，走向角分別為82°和112°.高分影像識別的門源地震地表破裂帶表明（劉璐等,2022），沿冷龍嶺斷裂西段破裂帶長度為21 km，沿托萊山斷裂東段破裂帶長度為5 km，說明本文選取的兩段斷層長度較為合理.根據(jù)GCMT震源機(jī)制解，門源地震矩形斷層傾角變化范圍初步設(shè)定為70°～90°，兩個斷層設(shè)置為等傾角，然后根據(jù)步長為1°搜索確定最佳傾角為80°.此外，InSAR形變揭示的門源地震地表破裂跡線與已知冷龍嶺斷裂和托萊山斷裂所處位置并不完全重合，其中Fault 2東側(cè)與冷龍嶺斷裂部分重合，F(xiàn)ault 2西側(cè)位于冷龍嶺斷裂和托萊山斷裂交界處，而Fault 1基本平行靠近于托萊山斷裂.

震中位置選取對于聯(lián)合反演中高頻GNSS位移波形擬合至關(guān)重要，本文根據(jù)中國地震臺網(wǎng)中心、USGS和Fan等（2022）提供的主震震中位置初步劃定格網(wǎng)搜索范圍，然后基于數(shù)據(jù)擬合程度在由地表破裂跡線構(gòu)建的斷層模型上反演搜索最優(yōu)震中位置.格網(wǎng)搜索結(jié)果如圖3所示，高頻GNSS位移波形擬合相較于InSAR對震中位置更為敏感，故本文基于高頻GNSS波形擬合選取最佳震中位置為101.26°E，37.79°N.該震中位于本文選取的矩形斷層Fault 2上，在已知斷層數(shù)據(jù)庫中位于冷龍嶺斷裂和托萊山斷裂交界處，并更接近于托萊山斷裂.

圖3 2022年門源M 6.9地震震中位置格網(wǎng)搜索（a） 高頻GNSS位移波形擬合方差降；（b） InSAR數(shù)據(jù)擬合方差降.Fig.3 Epicenter grid search of the 2022 Menyuan M 6.9 earthquake（a） Variance reduction of high-rate GNSS displacement waveforms;（b） Variance reduction of InSAR.

2.2 反演方法

在發(fā)震斷層幾何參數(shù)確定之后，為了細(xì)化斷層面滑動分布，將兩個斷層離散成2 km×2 km共計120個子斷層.然后采用非負(fù)最小二乘算法（Lawson and Hanson,1995）聯(lián)合InSAR和高頻GNSS位移波反演每個子斷層上的同震滑動分布，基本線性觀測方程可表示為：

Gm=d，

（1）

式中m表示每個子斷層上的滑動參數(shù)，G表示斷層位錯和地殼形變之間的線性響應(yīng)（格林函數(shù)），d表示InSAR和高頻GNSS位移波的觀測矩陣.每個子斷層用5個50%時間重疊窗口進(jìn)行參數(shù)化，每個時間窗口由2 s持續(xù)時間的三角形滑移率函數(shù)描述.基于門源地區(qū)一維地殼速度分層模型（表3）（尹欣欣等,2022），使用頻率-波數(shù)積分法（Zhu and Rivera,2002）計算InSAR和高頻GNSS位移波的格林函數(shù).

表3 震中附近的地殼速度模型Table 3 The crustal velocity model near the epicenter

此外，通過離散化子斷層可以將需要反演的無限數(shù)量參數(shù)減少到有限數(shù)量，雖然此時需要反演的參數(shù)數(shù)量小于觀測值數(shù)量，但由于在反演過程中位于深部的子斷層約束不佳以及部分區(qū)域觀測值分布稀疏等原因，滑動分布反演問題總是一個不適定問題，因此需要施加正則化約束（Wang et al.,2016）.聯(lián)合InSAR靜態(tài)位移和高頻GNSS動態(tài)位移波形進(jìn)行運(yùn)動學(xué)反演需要同時對空間和時間施加正則化約束，本研究采用一階拉普拉斯空間正則化矩陣Ls和時間一階導(dǎo)數(shù)正則化矩陣Lt來對該反演問題進(jìn)行約束（Melgar Moctezuma,2014），式（1）可以擴(kuò)展為：

（2）

式中λs和λt分別表示空間和時間正則化的平滑參數(shù)，最佳平滑參數(shù)由ABIC（Akaike′s Bayesian Information Criterion,Akaike貝葉斯信息準(zhǔn)則）方法（Akaike,1980;Fukahata et al.,2004）確定.

本研究嘗試?yán)酶哳lGNSS波形擬合殘差來確定地震最佳破裂速度，但由于GNSS臺站距離震中較遠(yuǎn)且該地震破裂尺度較小，高頻GNSS位移波形對破裂速度并不敏感，故本次地震破裂速度取經(jīng)驗值（0.8倍S波速度）2.8 km·s-1（Heaton,1990）.確定不同種類數(shù)據(jù)相對權(quán)重比以便每個數(shù)據(jù)集都能較好擬合對于多源數(shù)據(jù)聯(lián)合反演來說至關(guān)重要（Chen et al.,2020）.我們根據(jù)經(jīng)驗定權(quán)法（Feng et al.,2010）將InSAR和高頻GNSS權(quán)重比設(shè)置為1∶1，以確保能對每個數(shù)據(jù)集都能較好地解釋和擬合（Zhang et al.,2009）.

2.3 結(jié)果分析

最佳斷層滑動分布模型如圖3所示.結(jié)果表明，門源地震引起的斷層破裂長度約為26 km，與野外現(xiàn)場考察結(jié)果基本一致（潘家偉等,2022）.兩條發(fā)震斷層Fault 1和Fault 2均以左旋走滑運(yùn)動為主，其中Fault 2中部兼具少量逆沖成分.門源地震主要滑動沿深度主要分布在震源（5 km）上部，沿走向則主要集中在Fault 2的冷龍嶺斷裂部分，其中最大滑動也位于Fault 2，最大滑動量為～4 m，所處深度為～4 km.Fault 1存在一個較小的高滑區(qū)，最大滑動量為～2 m，對應(yīng)深度為～4 km.

圖4 2022年門源M 6.9地震同震滑動分布模型及其不確定性（a） 滑動分布在地表的投影，子圖為震源時間函數(shù)；（b） 同震滑動模型，藍(lán)色五角星為震中，藍(lán)色箭頭表示滑動方向；（c） Jackknife測試的方差系數(shù).Fig.4 Coseismic slip distribution model of the 2022 Menyuan M 6.9 earthquake and its uncertainty（a） Surface projection of the slip distribution,inset plot shows the source time function;（b） Coseismic slip model,the blue star is the epicenter,and the blue arrows indicate the slip direction; （c） The coefficient of variance from Jackknife test.

破裂演化圖和震源時間函數(shù)表明，門源地震為雙側(cè)不對稱破裂事件（圖4a，圖5）.地震成核后，破裂在Fault 2由震中向兩側(cè)傳播，并在3.8 s時達(dá)到地震矩釋放速率峰值.在5～8 s，地震破裂在Fault 2傳播終止，但是在8 s后又繼續(xù)釋放，推測是由于地震破裂在震源東南側(cè)遇到了某些阻礙導(dǎo)致破裂在該處幾乎中斷.西側(cè)Fault 1在5 s后地震矩繼續(xù)釋放，但速率快速下降.門源地震總共在～10 s內(nèi)釋放了8.65×1018Nm的地震矩，對應(yīng)矩震級為MW6.56.

最優(yōu)同震滑動模型對應(yīng)的InSAR形變場、高頻GNSS位移波如圖6、7所示.從圖6可以看出，反演模型預(yù)測的InSAR升降軌LOS向形變與觀測值吻合較好，總體InSAR擬合方差降（Variance Reductions,VR）為90.2%.從InSAR降軌LOS向擬合殘差分布可以看出，在近斷層區(qū)域仍有一些殘差分布，推測是由于簡化矩形斷層模型并不能反映真實斷層幾何的小規(guī)模細(xì)節(jié)，以及在最上層地殼中存在沉積層導(dǎo)致出現(xiàn)非彈性響應(yīng)等多種因素導(dǎo)致（Fialko et al.,2005;Xu et al.,2010）.

從圖7可以看出，GNSS水平向預(yù)測波形與觀測波形整體趨勢擬合較好，但部分高頻信息擬合較差，這些高頻信號可能屬于GNSS觀測噪聲，簡化的矩形斷層模型也可能導(dǎo)致了部分臺站擬合較差.

圖5 2022年門源M 6.9地震破裂1s間隔快照藍(lán)色五角星為震源，黑色虛線為破裂速度為2.8 km·s-1的破裂前緣等時線.Fig.5 Snapshot of the 2022 Menyuan M 6.9 earthquake rupture history with 1 s intervalThe blue star is the epicenter,the black dotted line is the isochron of the rupture front with a rupture velocity of 2.8 km·s-1.

圖6 InSAR同震形變場的觀測、擬合和殘差第一行和第二行分別表示升軌和降軌數(shù)據(jù)，藍(lán)色五角星表示震中.Fig.6 Coseismic InSAR deformation field,synthetics and residualsThe top and bottom rows show data from ascending and descending tracks,respectively.The blue star represents the epicenter.

圖7 水平向高頻GNSS觀測（黑色）和模擬（紅色）波形Fig.7 Horizontal high-rate GNSS recorded （black） and predicted （red） waveforms

為了評估同震滑動模型的不確定性，本研究對滑動模型中每個子斷層進(jìn)行了Jackknife測試（Kim and Dreger,2008）.首先隨機(jī)刪除20%反演數(shù)據(jù)集（InSAR升、降軌聯(lián)合數(shù)據(jù)隨機(jī)刪除260個點，高頻GNSS位移波數(shù)據(jù)隨機(jī)剔除3個臺站），然后利用每次刪除之后的數(shù)據(jù)集進(jìn)行100次重復(fù)反演，模型參數(shù)保持不變.根據(jù)這100次的反演結(jié)果，可以計算得到平均滑動模型以及每一個子斷層上的滑動標(biāo)準(zhǔn)差，利用平均滑動模型和滑動標(biāo)準(zhǔn)差可以計算方差系數(shù)（Coefficient of Variation,CV）.方差系數(shù)為滑動標(biāo)準(zhǔn)差和平均滑動的比值，可用來評估每個子斷層上滑動量的變化是否穩(wěn)健，方差系數(shù)越小，滑動量反演越可靠.結(jié)果如圖4c所示，主滑移區(qū)方差系數(shù)均小于0.2，表明斷層面上這部分區(qū)域滑動穩(wěn)定.

3 同震庫侖應(yīng)力變化

地震會改變周圍斷層剪應(yīng)力和正應(yīng)力，使得應(yīng)力增加區(qū)域地震活動率上升，應(yīng)力下降區(qū)域地震活動率下降（Stein,1999）.基于庫侖破裂準(zhǔn)則（King et al.,1994），庫侖應(yīng)力變化ΔCFS可表示為：

ΔCFS=Δτ+μ′Δσn，

（3）

式中Δτ和Δσn分別是接收斷層面上剪應(yīng)力和正應(yīng)力變化，μ′是有效摩擦系數(shù).李振洪等（2022）利用該原理研究了2016年門源地震和2022年門源地震之間的觸發(fā)關(guān)系，其結(jié)果表明2016年門源MW5.9地震增加了2022年地震發(fā)震斷層上的庫侖應(yīng)力變化，對2022年門源地震有一定觸發(fā)作用.類似地，本文利用庫侖破裂準(zhǔn)則研究2022年門源地震對余震分布和周圍區(qū)域庫侖應(yīng)力狀態(tài)的影響.根據(jù)King等（1994）的結(jié)果，本文將有效摩擦系數(shù)μ′設(shè)置為0.4，選擇10 km深度為研究對象.

圖8 2022年門源M 6.9地震庫侖應(yīng)力變化藍(lán)色五角星表示震中，白色圓點表示余震分布.Fig.8 Coulomb stress change caused by the 2022 Menyuan M 6.9 earthquakeThe blue star is the epicenter,and the white dots indicate aftershocks.

同震庫侖應(yīng)力變化如圖8所示.2022年門源地震對絕大多數(shù)余震分布在10 km深度處造成的庫侖應(yīng)力增加均超過了0.01 MPa，這說明主震對余震有一定的觸發(fā)作用.具體而言，本文根據(jù)InSAR地表破裂跡線選取的兩個矩形斷層與余震序列顯示的破裂面總長度（32 km）基本一致（Fan et al.,2022），但東西兩個斷層分段長度以及空間位置略有差異，推測是由于兩條發(fā)震斷層均為南傾，且主震使斷層～4 km處發(fā)生破裂，而余震則主要分布在～10 km，故淺部主震的發(fā)生可能觸發(fā)了深部斷層幾何結(jié)構(gòu)更為復(fù)雜的余震.此外，2022年門源地震也使門源斷裂西段和民樂—大馬營斷裂東段部分區(qū)域庫侖應(yīng)力增加超過0.01 MPa，使得這些區(qū)域未來發(fā)生地震概率增大.

4 討論

4.1 分段破裂特征

位于大陸內(nèi)部的活動斷層通常具有分段特征，斷層長度受限于當(dāng)?shù)氐貧ず穸龋虼嗣恳欢位顒訑鄬娱L度通常在10～30 km，而這種固有最大斷層尺寸限制了大陸地震震級，每一段活動斷層上潛在最大發(fā)震震級被限制在MW6～7之間（Pacheco et al.,1992;Scholz,1997;Stock and Smith,2000;Klinger,2010）.超過該震級閾值的大陸地震通常會有多個斷層發(fā)生破裂，比如2016年新西蘭凱庫拉M7.8地震至少有12個斷層發(fā)生破裂（Hamling et al.,2017），2021年青?，敹郙7.4地震發(fā)震斷層被分為了4段（Chen et al.,2022）.2022年門源地震產(chǎn)生了兩段斷層破裂，其中Fault 1長10 km，F(xiàn)ault 2長20 km，符合一般活動斷層分段長度.

地震破裂過程可以用多個具有阻礙或觸發(fā)地震破裂的障礙體來描述（Aki,1979）.這種地震破裂障礙體一般可分為兩類：幾何障礙體（geometric barriers）和松弛障礙體（relaxation barriers），幾何障礙體與斷層方向發(fā)生變化的位置相關(guān)，而松弛障礙體則是低應(yīng)力積累區(qū)域，震間淺層蠕變或者已有地震發(fā)生均會導(dǎo)致活動斷層應(yīng)力降低（King,1986;Lauer et al.,2020）.因此，幾何障礙體的出現(xiàn)可以阻礙或觸發(fā)在另一條斷層上的地震破裂，而松弛障礙體則一般會阻礙地震破裂.同震滑動分布模型（圖4a,4b）顯示在震源西北側(cè)發(fā)生明顯的斷層走向變化，滑動分布在該處也發(fā)生明顯變化：地震破裂在Fault 2向西北側(cè)傳播時破裂深度逐漸變小，但是當(dāng)破裂傳遞到Fault 1時，破裂深度有所增大.同時由動態(tài)破裂演化（圖5）可以看出，地震破裂在5 s之后主要發(fā)生在Fault 1上，這表明斷層走向改變和分段特征導(dǎo)致了地震在Fault 2上釋放完大部分能量之后觸發(fā)了相鄰斷層Fault 1的后續(xù)破裂，同時也說明Fault 1沒有阻礙地震的進(jìn)一步破裂，因此不是低應(yīng)力積累松弛障礙體區(qū).因此綜合斷層幾何學(xué)和地震破裂運(yùn)動學(xué)，推測該斷層走向變化導(dǎo)致的斷層分段為幾何障礙，它將震中西北側(cè)兩個走向不同的斷層分開，這使得地震在震源兩側(cè)呈非對稱破裂特征.

4.2 無震蠕滑與淺層滑動虧損

活動斷層中淺層無震蠕滑段的存在限制了該部分?jǐn)鄬釉杏蟮卣鸬臐摿Γˋvouac,2015）.海原斷裂帶上260 km長的天祝地震空區(qū)因其具有孕育大地震風(fēng)險而一直備受關(guān)注（Gaudemer et al.,1995）.該段斷層自西向東可以分為4個長度相當(dāng)?shù)牟糠郑↙iu-Zeng et al.,2007）：冷龍嶺段、金強(qiáng)河段、毛毛山段和老虎山段.其中利用InSAR觀測發(fā)現(xiàn)老虎山段存在一條長約35 km的無震淺層蠕滑段（Cavalié et al.,2008;Jolivet et al.,2012;Li et al.,2021），其主要特征是在該段區(qū)域斷層兩側(cè)速度梯度急劇變化，而西部冷龍嶺斷裂和托萊山斷裂的相關(guān)研究較少.Huang等（2022）利用Sentinel-1數(shù)據(jù)繪制的2014—2021年InSAR形變速度圖首次覆蓋了海原斷裂帶西部的冷龍嶺斷裂和托萊山斷裂，升降軌LOS向速度場顯示這兩條斷裂兩側(cè)速度梯度變化較為平緩，表明在近年來該段沒有明顯淺層蠕變發(fā)生.

在多個地震周期中，斷層不同深度處的滑動量（同震、震后和震間滑動的總和）應(yīng)保持一致（Tse and Rice,1986;Xu et al.,2016）.然而利用同震地表形變數(shù)據(jù)反演大型（MW～7）走滑地震結(jié)果表明，地殼最上部的同震滑動通常系統(tǒng)性地小于發(fā)震深度（4～10 km）的同震滑動（如，Simons et al.,2002;Fialko et al.,2005;Dolan and Haravitch,2014），這種現(xiàn)象被稱為“淺層滑動虧損”（Shallow Slip Deficit,SSD）.聯(lián)合滑動分布模型顯示門源地震的滑動集中分布在震源上方（5 km以上），計算可得SSD為6.3%，小于之前發(fā)生的大多數(shù)大型走滑地震（圖9）.造成SSD的原因目前仍有較大爭議，一種可能的解釋是由于淺部斷層在震間的無震蠕滑導(dǎo)致（Hussain et al.,2016;Xu et al.,2016）.而冷龍嶺斷裂和托萊山斷裂在近8年來并沒有觀測到明顯的淺層蠕變（Huang et al.,2022）.結(jié)合2022年門源地震相對較小的SSD，可以推測這兩條斷裂淺部斷層（0～2 km）應(yīng)力積累大部分仍以地震的形式釋放，在更長震間期內(nèi)發(fā)生淺層蠕變的概率較小.同時由4.1節(jié)討論的地震分段破裂特征可知，西側(cè)與Fault 1相近的托萊山段不屬于低應(yīng)力松弛障礙體區(qū)，而是一直處于高應(yīng)力積累狀態(tài)，這也表明托萊山斷裂在震間期幾乎不存在淺層蠕滑.同時在該區(qū)域也沒有現(xiàn)代地震活動記錄，因此推測托萊山斷裂未來具有孕育大型地震的潛力，需要對其進(jìn)行持續(xù)的地震危險性跟蹤監(jiān)測.

圖9 走滑地震累計同震滑動分布，數(shù)據(jù)來源：Fialko et al.（2005）,Jin and Fialko （2021）Fig.9 Cumulative coseismic slip distribution of strike-slip earthquakes,data source from Fialko et al.（2005） and Jin and Fialko （2021）

5 結(jié)論

本文利用InSAR和高頻GNSS位移波形數(shù)據(jù)研究了2022年1月8日門源M6.9地震的同震形變和運(yùn)動學(xué)破裂過程，主要結(jié)論如下：

（1）使用DInSAR和POT從Sentinel-1 SAR衛(wèi)星影像中提取的形變場顯示，最大LOS向（距離向）同震形變達(dá)0.79 m，地表破裂跡線揭示了走向角分別為82°和112°的兩個發(fā)震斷層.

（2）門源地震為雙側(cè)不對稱破裂，主要滑動發(fā)生在Fault 2的冷龍嶺段，最大滑動量為～4 m，所處深度為～4 km，西北側(cè)Fault 1也存在一個最大滑動量為～2 m的小范圍高滑區(qū).震源時間函數(shù)表明地震在～10 s內(nèi)釋放了8.65×1018Nm地震矩，對應(yīng)矩震級為MW6.56.

（3）本次地震也使門源斷裂西段和民樂—大馬營斷裂東段部分區(qū)域的庫侖應(yīng)力增加超過了0.01 MPa，應(yīng)密切關(guān)注這些區(qū)域未來發(fā)生強(qiáng)震的可能性.

（4）門源地震SSD僅6.3%，且冷龍嶺斷裂和托萊山斷裂在近8年沒有明顯淺層蠕變，推測這兩條斷裂在震間更長時間尺度范圍內(nèi)發(fā)生淺層蠕變概率較小.另外，震源西北側(cè)幾何障礙體導(dǎo)致門源地震分段破裂，表明托萊山斷裂屬于高應(yīng)力積累區(qū)域，未來有孕育中強(qiáng)地震的風(fēng)險.

致謝感謝歐空局（ESA）提供Sentinel-1 A衛(wèi)星影像和中國地震局地球物理研究所范莉蘋博士提供余震重定位目錄.開源軟件ISCE2的使用得到了北京大學(xué)梁存任教授指導(dǎo)，論文部分圖件由GMT繪制.