2021年雙湖MS5.8地震InSAR同震形變場及斷層滑動分布反演

李承濤 李 琦 譚 凱 魯小飛

1 中國地震局地震研究所，武漢市洪山側路40號，4300712 中國地震局地震大地測量重點實驗室，武漢市洪山側路40號，4300713 湖北省地震局，武漢市洪山側路48號，4300714 河北省地震動力學重點實驗室，河北省三河市學院街465號，065201

據中國地震臺網中心(CENC)測定，2021-03-30西藏那曲市雙湖縣(87.68°E, 34.38°N)發生MS5.8地震，震源深度10 km。美國地質調查局(USGS)提供的體波矩張量解顯示，該地震震級為MW5.56，震源深度4 km。該地震震中位于北羌塘盆地，鄰近琵琶湖-吐坡錯斷裂與琵琶湖-映天湖斷裂交匯處[1]。該地區1976年以來共發生38次MW≥5地震(https:∥www.globalcmt.org/CMTsearch.html)，其中大部分地震表現為正斷性質，近年發生的最大地震為2020-07-23西藏尼瑪MS6.6地震(距離本次地震約155 km)(圖1)。該地震震中附近5 km范圍內平均海拔約5 023 m，交通不便，缺少GNSS等觀測數據及地質資料。本文對Sentinel-1A數據進行DInSAR處理，獲取同震形變場，通過基于序列蒙特卡羅采樣的貝葉斯方法反演發震斷層的幾何參數[2]和同震滑動分布，可為快速、準確認識雙湖MS5.8地震的震中位置及同震特征提供參考。

圖1 2021年雙湖MS5.8地震構造背景Fig.1 Tectonic setting of the 2021 Shuanghu MS5.8 earthquake

1 InSAR數據處理

本文采用歐空局官方網站(https:∥sentinel.esa.int/)提供的Sentinel-1A干涉寬幅數據，具體信息見表1。

表1 Sentinel-1A數據相關參數

基于ISCE(InSAR scientific computing environment)軟件對Sentinel-1A數據進行DInSAR處理[3]。在具體處理流程中，距離向和方位向按5∶1進行多視處理，生成干涉圖；利用精密軌道文件(https:∥scihub.copernicus.eu/gnss/#/home)和SRTM(shuttle radar topography mission)數字高程模型[4]，分別消除軌道和地形影響；為提高信噪比，采用Goldstein方法進行濾波[5]；相位解纏采用軟件SNAPHU(statistical-cost, network-flow phase-unwrapping algorithm)[6]；通過地理編碼，獲取WGS84坐標系下雙湖地震的LOS同震形變場，其中對相關性低值區進行掩膜處理。從雙湖地震的同震干涉圖(圖2)可以看出，形變長軸呈近NS向展布，影響范圍約10 km×10 km。

圖2 雙湖地震同震干涉圖Fig.2 Coseismic interferograms of Shuanghuearthquake

2 斷層參數和滑動分布反演

2.1 斷層參數反演

為獲取發震斷層相關的幾何形狀參數，利用基于序列蒙特卡羅采樣的貝葉斯方法進行反演[2,7-8]。本文采用四叉樹算法對LOS同震形變場進行降采樣處理[9]，最終得到426個數據點，其中T114軌道228個，T121軌道198個，計算過程共迭代1.9×106次。圖3為斷層幾何參數1D和2D的后驗概率密度函數(posterior probability density function, PDF)，可表征幾何參數的離散程度。選取最大后驗解(maximum a posterior solutions, MAP)作為最佳斷層參數：東偏移量-0.452 km、北偏移量7.025 km，表示的是參考GCMT震中(87.7° E, 34.32° N)的偏移量，代表斷層上邊界中心點所處的位置，經坐標轉換后為87.695° E、34.383° N，相應深度為3.662 km，斷層長度、寬度分別為4.271 km和6.461 km，斷層走向、傾角、滑動角分別為47.810°±2.218°、36.664°±2.499、-66.598°±3.258°。具體參數見表2。

表2 反演的斷層幾何形狀參數

基于上述斷層模型得到的LOS殘差分布如圖4所示，圖中淺灰色散點表示采樣點，黑色直線與矩形分別表示斷層上邊界和斷層面的位置。從圖中可以看出，升降軌均顯示以沉降為主，最大沉降量分別約5.8 cm和4.8 cm。對于T114升軌，殘差范圍為-0.9～0.9 cm，殘差均方根(RMS)為0.4 cm；對于T121降軌，殘差范圍為-0.6～0.7 cm，殘差RMS為0.3 cm，表明反演獲取的斷層幾何參數較為合理。

圖4 斷層幾何參數反演獲得的LOS形變殘差Fig.4 LOS deformation residuals of fault geometry inversion

2.2 斷層滑動分布反演

在已知斷層幾何參數的基礎上，將斷層離散為若干子斷層，斷層滑動分布反演可簡化為求解每個子斷層對應的滑動。利用基于Python語言的Pyrocko-GF軟件生成格林函數庫，然后通過基于序列蒙特卡羅采樣的貝葉斯方法反演滑動分布[2,10]，子斷層滑動分布與相應格林函數庫的乘積即為地表位移。每個子斷層的滑動量搜索范圍設置為0～1.0 m。

在貝葉斯方法中，正則化可通過高斯先驗p(s,α)實現，其包含一個非對角項的協方差矩陣，Cholesky分解等價于一個大小為M×M的拉普拉斯有限差分算子L(M為子斷層個數)。p(s,α)可表示為[2]：

(1)

正則化后的后驗概率密度函數PDF可表示為[11]：

p(m,σ1,σ2,…,σK|dobs,1,dobs,2,…,dobs,K)

(2)

式中，m為已知斷層參數，σK為第K個超參數，dobs,K為第K個觀測數據，平滑因子α可作為待求參數，s為子斷層待求的滑動矢量。

直方圖中紅色直線和2D相關圖中的紅點表示最大后驗解，藍色表示高概率區域圖3 雙湖地震矩形震源的1D和2D后驗概率密度分布Fig.3 1D and 2D PDF of a rectangular source for Shuanghu earthquake

在反演斷層滑動分布時，參考表2結果對斷層參數進行初始化設置，其中走向取47.810°，傾角取36.664°，延伸斷層長、寬分別為15 km和19 km，各子斷層大小為1 km×1 km，平均滑動角固定為-66.598°。

反演采用的數據共426個，與前文反演斷層幾何參數所使用的數據個數一致。從同震滑動分布(圖5)可以看出，主要滑動區集中在深度3.9～7.5 km，最大滑動量約0.42 m，平均滑動角結果表明，雙湖地震同震破裂主要表現為正斷拉張兼少量左旋走滑性質。反演得到的震中位置為87.715° E、34.365° N，震源深度為5.763 km。假設地殼的剪切模量為30 GPa，計算得到地震釋放的地震矩約為2.189×1017Nm，矩震級為MW5.5，略小于USGS提供的震級。

圖5 基于InSAR數據反演的同震滑動分布Fig.5 Coseismic slip distribution constrained by InSAR

圖6為同震滑動分布模型的殘差分布，圖中小矩形表示斷層位置，大矩形表示拓展長寬后得到的斷層面。從圖中可以看出，T114升軌的殘差范圍為-1.2～1.5 cm， T121降軌的殘差范圍為-0.5～0.9 cm，升降軌殘差RMS分別為0.7 cm和0.3 cm，表明滑動分布結果具有合理性與可靠性。

圖6 斷層滑動分布模型的LOS形變殘差Fig.6 The LOS deformation residuals of thefault slip distribution model

3 討 論

表3為不同機構提供的震源參數，由于采用的數據和方法不同，獲得的參數存在差異。反演結果表明，發震斷層屬于具有正斷性質的隱伏斷層，震中位于琵琶湖-吐坡錯斷裂與琵琶湖-映天湖斷裂的交匯處，地震可能受這2個斷裂共同作用。Sentinel-1衛星對南北方向不敏感，當南北向分量近似忽略時，通過升降軌數據可以獲取近東西向和垂直向的形變分量，稱為2.5D形變場[12]，結果如圖7所示，圖中灰色矩形框表示滑動分布的斷層面，黑色和紅色五角星分別表示CENC和InSAR計算的震中位置，后者更靠近形變中心。從圖7可以看出，東向運動最大值約2.8 cm，西向運動最大值約2.2 cm；最大沉降量約6.4 cm，位于斷層上盤附近。2.5D形變場結果表明，雙湖地震以正斷特性為主，發震斷層屬于東傾斷層。此次地震反演得到的斷層傾角較小，表明地震的成因不僅受自身重力影響，還可能受到區域拉張作用的影響，具體情況需要作進一步研究。雙湖地區是羌塘盆地最具油氣資源潛力的地區之一[13]，本次地震確定的地下隱伏斷層相關參數有利于區域內部結構的探測，可為相關資源的開采提供參考。

表3 不同機構給出的震源參數

圖7 2.5D形變場Fig.7 2.5D deformation filed

4 結 語

本文基于Sentinel-1A干涉寬幅數據，利用DInSAR技術提取雙湖地震的同震形變場，并通過SMC方法反演發震斷層參數及同震滑動分布，得到以下結論：

1)T114升軌、T121降軌的LOS形變場均顯示雙湖地震以沉降為主，升降軌LOS最大沉降量分別約5.8 cm和4.8 cm。

2)反演得到發震斷層走向約47.810°±2.218°，傾角約36.664°±2.499。

3)斷層滑動分布結果表明，雙湖地震主要破裂區集中在3.9～7.5 km深度，最大滑動量約0.42 m，反演確定的震中位置為87.715° E、34.365° N，震源深度為5.673 km，平均滑動角約-66.598°±3.258°，表明本次地震以正斷拉張為主，兼少量左旋走滑性質。假設地殼的剪切模量為30 GPa，計算得出的矩震級為MW5.5，略小于USGS的結果。

致謝：感謝歐洲空間局(ESA)提供Sentinel-1衛星升降軌SAR數據。