李承濤 李 琦 譚 凱 魯小飛

1 中國(guó)地震局地震研究所，武漢市洪山側(cè)路40號(hào)，4300712 中國(guó)地震局地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室，武漢市洪山側(cè)路40號(hào)，4300713 湖北省地震局，武漢市洪山側(cè)路48號(hào)，4300714 河北省地震動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，河北省三河市學(xué)院街465號(hào)，065201

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

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

1 InSAR數(shù)據(jù)處理

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

表1 Sentinel-1A數(shù)據(jù)相關(guān)參數(shù)

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

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

2 斷層參數(shù)和滑動(dòng)分布反演

2.1 斷層參數(shù)反演

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

表2 反演的斷層幾何形狀參數(shù)

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

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

2.2 斷層滑動(dòng)分布反演

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

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

(1)

正則化后的后驗(yàn)概率密度函數(shù)PDF可表示為[11]：

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

(2)

式中，m為已知斷層參數(shù)，σK為第K個(gè)超參數(shù)，dobs,K為第K個(gè)觀測(cè)數(shù)據(jù)，平滑因子α可作為待求參數(shù)，s為子斷層待求的滑動(dòng)矢量。

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

在反演斷層滑動(dòng)分布時(shí)，參考表2結(jié)果對(duì)斷層參數(shù)進(jìn)行初始化設(shè)置，其中走向取47.810°，傾角取36.664°，延伸斷層長(zhǎng)、寬分別為15 km和19 km，各子斷層大小為1 km×1 km，平均滑動(dòng)角固定為-66.598°。

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

圖5 基于InSAR數(shù)據(jù)反演的同震滑動(dòng)分布Fig.5 Coseismic slip distribution constrained by InSAR

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

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

3 討 論

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

表3 不同機(jī)構(gòu)給出的震源參數(shù)

圖7 2.5D形變場(chǎng)Fig.7 2.5D deformation filed

4 結(jié) 語

本文基于Sentinel-1A干涉寬幅數(shù)據(jù)，利用DInSAR技術(shù)提取雙湖地震的同震形變場(chǎng)，并通過SMC方法反演發(fā)震斷層參數(shù)及同震滑動(dòng)分布，得到以下結(jié)論：

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

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

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

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