2021年云南省漾濞Ms6.4地表二維形變分析

潘建平，趙瑞淇，袁雨馨，蔡卓言，付占寶

(重慶交通大學 智慧城市學院，重慶 400074)

據(jù)中國地震臺網(wǎng)消息，北京時間2021年5月21日T21:48云南省大理州漾濞縣(99.87°N，25.67°E)發(fā)生Ms6.4地震，震源深度8 km，地震致3人遇難，37人受傷，并造成不同程度損毀。漾濞地震發(fā)生在橫斷山滇西高山峽谷區(qū)，地質構造復雜，地震活動頻繁，震中附近斷裂構造復雜，斷裂眾多[1]。因此，確定本次地震的同震形變場，對認識該地區(qū)的地震孕震環(huán)境、同震滑動模型反演以及地震風險防范具有重要意義，過去對地表形變監(jiān)測大都基于水準儀、全站儀、GNSS等設備的外業(yè)手段，不能在較短時間內獲得大范圍的、監(jiān)測點分布密集的觀測數(shù)據(jù)[2],因而不適用于地震這種大范圍高烈度的自然災害形變監(jiān)測。

合成孔徑雷達差分干涉測量(differential synthetic aperture radar interferometry，D-InSAR)具有大范圍，高精度以及低成本的優(yōu)勢，可以監(jiān)測到cm級或mm級的形變[3]。1993年，測量美國加州Landers地震形變場首次應用D-InSAR技術[4]，從此D-InSAR技術被廣泛地應用到地震相關研究中，特別是在震間以及同震形變領域的應用[5]：張慶云[6]等采用D-InSAR技術，與ALOS-2和Sentinel-1A數(shù)據(jù)聯(lián)合使用得到了2017年伊朗Mw7.3地震的同震形變場，并在此基礎上反演出震發(fā)斷層的詳細幾何參數(shù)。Elias[7]等采用D-InSAR技術提取希臘Samos島近海M7.0級地震的同震形變，并聯(lián)合GNSS偏移量進行了斷層參數(shù)反演。Zhao L[8]等利用ALOS-2以及Sentinel-1A/B數(shù)據(jù)得到了菲律賓棉蘭老島4次Mw>6.0地震的InSAR同震形變場，并在InSAR監(jiān)測結果的基礎上對斷層幾何形狀和同震滑動分布進行了幾何學約束。Qingyun Zhang[9]等利用InSAR監(jiān)測結果與強震數(shù)據(jù)和GPS數(shù)據(jù)結合得到了2016日本熊本地震的三維形變場。本次地震發(fā)生后，我國學者立即采用InSAR技術開展研究，王紹俊[10]等基于D-InSAR技術獲得了此次地震的升降軌同震形變場，并以InSAR形變場和GNSS數(shù)據(jù)為約束，反演了斷層滑動分布。楊九元[11]等利用Sentinel-1A/B數(shù)據(jù)的InSAR同震形變反演了此次地震詳細的斷層幾何參數(shù)與滑動分布，并對該區(qū)域的地震危險性進行了評估。黎朕靈[12]等利用InSAR升降軌同震形變反演了該地震的震源滑動模型，并與GPS同震位移數(shù)據(jù)進行了聯(lián)合反演。此外，另有其他學者利用地震學相關理論對本次地震開展了震源機制反演以及深部孕震背景探討[13-14]。針對此次地震，文中利用歐空局Sentinel-1A升降軌數(shù)據(jù)，分別采用D-InSAR技術解算得到了視線向(line of sight,LOS)的同震形變場，并分別進行了形變方向和形變量級的分析。利用美國地質調查局USGS給出的震源機制解和魏云杰[15]等人基于梯度最速下降法得到的同震滑動分布結果模擬計算了本次地震并探討了本次地震引發(fā)的周邊斷裂庫侖應力變化并分析了發(fā)震危險性。針對本次地震的現(xiàn)有研究在InSAR數(shù)據(jù)源上僅利用LOS向的同震形變場，但是LOS向的形變場不能完全反應地表真實形變，因而不能為進一步開展震源參數(shù)和同震滑動分布反演提供較好的約束條件，因此文中聯(lián)合升降軌同震形變場，使用直接解算法得到了垂直向和東-西向的同震形變場,為進一步開展此次地震震源參數(shù)和同震滑動模型精細化反演提供更可靠的約束條件。

1 研究方法

1.1 D-InSAR技術原理及數(shù)據(jù)處理

D-InSAR技術是在InSAR技術基礎上發(fā)展而來，利用衛(wèi)星重復軌道獲取的兩個時相的影像進行差分干涉，得到的干涉相位組成如式(1)所示。

φ=φref+φtop+φdef+φatm+φnoi.

(1)

式中:φref為參考橢球面相位;φtop為地形相位;φdef為形變相位;φatm為大氣延遲相位;φnoi為噪聲相位。文中利用瑞士SARMAP公司的Sarscape5.2.1軟件采用二軌法來提取形變相位，利用衛(wèi)星姿態(tài)參數(shù)以及美國宇航局發(fā)布的30 m分辨率的SRTM(shuttle radar topography mission)DEM分別去除參考橢球面相位和地形相位，采用10∶2的多視比抑制斑點噪聲并提高信噪比，采用GOLDSTEIN濾波方法消除部分噪聲相位的影響，這種濾波方法濾波器可變，提高了干涉條紋的清晰度。設置解纏閾值為0.2，使用最小費用流法(minimum cost flow, MCF)對干涉相位進行相位解纏[16]，最后通過相位轉形變和地理編碼得到了地理坐標系下的形變信息。

1.2 二維形變提取方法

利用一對升軌數(shù)據(jù)或者降軌數(shù)據(jù)提取到的形變信息只能反演單一的雷達視線向形變情況。單純的升軌或降軌數(shù)據(jù)獲取的沿LOS向的形變量相差可能較大，甚至完全相反，因此單一方向上的形變數(shù)據(jù)不能有效反映測區(qū)地表真實三維形變情況[17]。故本文利用升降軌數(shù)據(jù)提取到的形變信息進行聯(lián)合解算。升降軌側視成像幾何如圖1所示。利用常規(guī)D-InSAR技術提取到的形變是雷達視線方向上的一維形變，是地表東-西向、南-北向以及垂直向形變在視線向上的投影，由成像幾何與投影關系可以得到式(2)。

dLOS=dUcosθ-dNcoaαsinθ-dEsinαsinθ.

(2)

式中:dLOS代表視線向形變；dU、dN、dE分別代表地表形變在垂直、南-北和東-西向上的分量；θ代表雷達平均入射角;α代表衛(wèi)星平均飛行方位角。由式(2)可知要求解地表三維方向上的形變dU,dN和dE需要至少3個不同成像視角的LOS向形變信息，其聯(lián)合求解式為：

(3)

目前所有的SAR衛(wèi)星都是以近極地軌道的方式運行的，雷達視線方向對南-北向的形變不敏感，使獲取的南-北向形變質量很差[17],且受制于該區(qū)域Sentinel-1A數(shù)據(jù)的成像幾何，故文中研究內客忽略南-北向的形變，利用升軌和降軌兩個不同成像視角的LOS向形變信息對東-西向和垂直向的形變進行求解，其求解式為：

(4)

(5)

圖1 升降軌側視成像幾何圖

2 研究區(qū)概況及數(shù)據(jù)選取

2.1 研究區(qū)概況

漾濞地震震中距離漾濞縣城約7 km,距離大理市約40 km,其震中附近區(qū)域地質構造復雜，且發(fā)生過多次地質災害。漾濞地震震中位于川滇菱形塊體的西南邊界，是維西-喬后-巍山斷裂與紅河斷裂銜接處，震中區(qū)域由3條左階右旋斷裂組成，自西向東依次是維西-喬后-巍山斷裂、紅河斷裂北段和紅河斷裂中段[10]，震中區(qū)域高程信息以及主要斷裂帶分布如圖2所示，圖2中黑色虛線框為文中同震形變解算的參考區(qū)域。

圖2 震中主要斷裂帶分布及影像覆蓋范圍

2.2 數(shù)據(jù)選取

文中采用Sentinel-1A衛(wèi)星TOPSAR(terrain observation with progressive scan,SAR)條帶掃描模式獲取數(shù)據(jù)，數(shù)據(jù)覆蓋范圍為250 km，分辨率為5 m×20 m。文中實驗數(shù)據(jù)具體參數(shù)如表1所示。文中除了采用Sentinel-1A數(shù)據(jù)外，還采用歐空局提供的DORIS精密軌道數(shù)據(jù)AUX_POEORB進行軌道誤差校正。升降軌影像覆蓋范圍在圖2中用藍色矩形框表示。

表1 影像信息

3 同震形變分析3.1 視線向同震形變分析

將D-InSAR相位信息轉換為視線向的形變數(shù)據(jù)，得到了降軌和升軌兩個視線向形變場，其中降軌表現(xiàn)為兩個主要形變區(qū)域A和B，A和B大致呈橢圓形狀對稱分布，A表現(xiàn)為遠離視線向的沉降，最大沉降量約為-7 cm，B區(qū)域表現(xiàn)為靠近雷達視線向的隆升，最大隆升量約為8 cm，這一結果與王紹俊[10]以及黎朕靈[12]等人的InSAR解算結果基本一致，降軌同震形變場如圖3所示，在A與B之間形成了一個近北-西走向的沒有明顯形變的長軸，如圖3中紅色虛線所示，有可能為本次地震發(fā)震斷層所在的空間位置，震中附近已知的維西-巍山斷裂位于川滇塊體西部邊緣，北起雪龍山東麓，止于點蒼山西南，走向北北西，但該北-西走向長軸位于維西-巍山斷裂以西，不在其主斷裂上，楊九元[11]等的研究認為此次地震可能破裂在維西-巍山斷裂的隱伏分支斷層或一個獨立的未知的隱伏主斷層上，但是準確的地震斷層還需要進一步研究確認。

圖3 降軌同震形變場

根據(jù)中國地震臺網(wǎng)發(fā)布的相關數(shù)據(jù)顯示，自2021年5月21日T21：48主震發(fā)生后到降軌輔影像成像時間5月22日T23：14止，震區(qū)陸續(xù)共發(fā)生31次余震，其中M3.0以上的余震29次，該時間段內的余震時間序列如圖4(a)所示。將余震空間分布疊加到降軌形變場上，如圖4(b)所示，不難發(fā)現(xiàn)余震發(fā)震位置與降軌形變場高度吻合，故文中得到的InSAR降軌形變場包含上述余震的貢獻。

圖4 余震時間分布及空間分布

為了更好地分析降軌同震形變場，繪制了aa*、bb*、cc*3條剖面線，其分別位于沉降區(qū)、隆升區(qū)、跨越沉降和隆升區(qū)。3條剖面線的位置已在圖3中標明，其中剖面線aa*長度為12.7 km，剖面線bb*長度為10.3 km，剖面線cc*長度為12.5 km。沿3條剖面線的形變情況如圖5所示。由剖面線形變曲線可知:剖面線aa*在2.5～4.5 km中形成了1個約5 cm的大型沉降漏斗，而且在4.3 km處形變量達到了-5.6 cm。在7.6～8.7 km中形成了一個-4 cm左右的小型沉降漏斗，在8.1 km處沉降量達到了-4.3 cm；剖面線bb*在1～6.8 km內隆升量均超過7 cm，隆升量較大，其余部分隆升量均小于7 cm；剖面線cc*在6 km處為北-西走向長軸所在位置，在6 km之前表現(xiàn)為沉降，沉降量在-0.3～-3.8 cm，6 km之后表現(xiàn)為隆升，隆升量在0～8 cm。

升軌同震形變場因受山區(qū)獨特地理環(huán)境以及衛(wèi)星成像幾何的影響，噪聲污染嚴重，形變特征不明顯，但不難發(fā)現(xiàn)A和B兩個形變區(qū)域明顯擴大，且兩個區(qū)域的形變方向與降軌結果相反，其中A區(qū)域表現(xiàn)為隆升，最大形變量為7 cm，B區(qū)域表現(xiàn)為沉降，最大形變量為-7 cm，如圖6所示，相同區(qū)域表現(xiàn)出不同的形變方向，與升降軌不同的成像視角有關。自降軌輔影像成像時間2021年5月22日T23：14起至升軌輔影像成像時間6月1日止，震中附近區(qū)域發(fā)生多次余震，其中M3.0以上共7次(據(jù)中國地震臺網(wǎng))，故升軌解算結果應包含了上述余震貢獻量，形變區(qū)域明顯擴大與余震形變累積有關。

圖5 剖面線形變曲線

圖7 同震庫侖應力變化

圖6 升軌同震形變場

3.2 同震庫侖應力變化計算

漾濞地震發(fā)生在維西-巍山斷裂和紅河斷裂交會處，該區(qū)域的地震背景活動性較強，處于背景地震預測模型的低值區(qū)，地震風險較高[10]。本次地震的發(fā)震斷層位于維西-巍山斷裂以西，且周圍區(qū)域斷層分布密集，主要有紅河斷裂北段和紅河斷裂中段，且這些斷層均為右旋走滑斷層，因此有必要計算本次地震引起的同震庫侖應力變化，對周圍斷層的地震危險性做出評估。

地震發(fā)生后，USGS給出了此次地震的震源機制解：走向135°、傾角82°、滑動角-165°；魏云杰等[15]給出了此次地震的平均同震滑動量0.19 m，發(fā)震斷層以右旋走滑為主。選擇上述參數(shù)作為同震庫侖應力計算的輸入?yún)?shù)，設置摩擦系數(shù)為0.4(走滑型斷層的典型值)?；趶椥园肟臻g位錯模型，利用Coulomb3.3軟件分別計算了深度為7.5 km和10 km的同震庫侖應力變化，如圖7所示。

同震庫侖應力變化計算結果顯示，此次地震對周邊斷層主要表現(xiàn)為庫侖應力卸載作用，僅在斷層東北部和東南部表現(xiàn)為加載作用，東北部對應為紅河斷裂北段，東南部沒有已知斷裂分布，因此此次地震引起的同震庫侖應力變化對紅河斷裂北段具有加載作用，但是沒有達到觸發(fā)值(0.1 bar)。盡管該結果可以為地震危險性初步評估提供一定的參考，但是精確的震源參數(shù)以及同震滑動分布可以為同震庫侖應力計算提供更可靠的數(shù)據(jù)源，進而可以實現(xiàn)對周圍斷裂發(fā)震危險性的精確評估，因此精確的GNSS數(shù)據(jù)和多維度的同震形變場可以為此提供更好的約束條件。

3.3 東-西向和垂直向同震形變分析

結合升降軌同震形變場解算結果，利用式(5)解算得到了東-西向和垂直向同震形變場，規(guī)定東-西向同震形變場形變量正值表示向東運動，負值表示向西運動，垂直向同震形變場形變量正值表示垂直方向的隆升，負值表示垂直方向的沉降。兩個方向的形變場呈現(xiàn)出C和D兩個主要形變區(qū)域，其中C區(qū)域面積大約為60 km2，D區(qū)域面積大約為40 km2，如圖8所示。

圖8 垂直向和東-西向同震形變場

圖8(a)顯示，C區(qū)域表現(xiàn)為向東運動，最大運動量達到了5.7 cm，D區(qū)域表現(xiàn)為向西運動，最大運動量達到了-7 cm。圖8(b)顯示，C區(qū)域表現(xiàn)為隆升，最大隆升量為3.9 cm，D區(qū)域表現(xiàn)為沉降，最大沉降量為-8 cm。

4 結束語

文中利用Sentinel-1A數(shù)據(jù)，采用二軌差分D-InSAR技術解算得到了2021年5月21日云南漾濞地震升軌和降軌的視線向同震形變場，并計算了此次地震引起的同震庫侖應力變化。最后結合升降軌同震形變場直接解算得到了本次地震的東-西向和垂直向同震形變場，可以對進一步開展此次地震震源參數(shù)和同震滑動模型精細化反演提供更可靠的約束條件，主要形成了以下幾點結論：①視線向同震形變場中，降軌的最大隆升量為8 cm，最大沉降量為-7 cm。升軌的最大隆升量為7 cm，最大沉降量為7 cm。②同震庫侖應力變化計算結果顯示此次地震僅對地震斷層東北部的紅河斷裂北段有加載作用，但加載量較小，發(fā)震危險性較低。③東-西向和垂直向同震形變場呈現(xiàn)兩個主要形變區(qū)域C和D，形變區(qū)域C整體表現(xiàn)為垂直向向上隆升，水平方向向東運動，水平方向最大形變量為5.7 cm，垂直向最大隆升量為3.9 cm。形變區(qū)域D整體表現(xiàn)為垂直向向下沉降，水平方向向西運動，水平方向最大形變量為-7 cm，垂直向最大沉降量為-8 cm。