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

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

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

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

1 研究方法

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

D-InSAR技術(shù)是在InSAR技術(shù)基礎(chǔ)上發(fā)展而來，利用衛(wèi)星重復(fù)軌道獲取的兩個時相的影像進行差分干涉，得到的干涉相位組成如式(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濾波方法消除部分噪聲相位的影響，這種濾波方法濾波器可變，提高了干涉條紋的清晰度。設(shè)置解纏閾值為0.2，使用最小費用流法(minimum cost flow, MCF)對干涉相位進行相位解纏[16]，最后通過相位轉(zhuǎn)形變和地理編碼得到了地理坐標(biāo)系下的形變信息。

1.2 二維形變提取方法

利用一對升軌數(shù)據(jù)或者降軌數(shù)據(jù)提取到的形變信息只能反演單一的雷達視線向形變情況。單純的升軌或降軌數(shù)據(jù)獲取的沿LOS向的形變量相差可能較大，甚至完全相反，因此單一方向上的形變數(shù)據(jù)不能有效反映測區(qū)地表真實三維形變情況[17]。故本文利用升降軌數(shù)據(jù)提取到的形變信息進行聯(lián)合解算。升降軌側(cè)視成像幾何如圖1所示。利用常規(guī)D-InSAR技術(shù)提取到的形變是雷達視線方向上的一維形變，是地表東-西向、南-北向以及垂直向形變在視線向上的投影，由成像幾何與投影關(guān)系可以得到式(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)星都是以近極地軌道的方式運行的，雷達視線方向?qū)δ?北向的形變不敏感，使獲取的南-北向形變質(zhì)量很差[17],且受制于該區(qū)域Sentinel-1A數(shù)據(jù)的成像幾何，故文中研究內(nèi)客忽略南-北向的形變，利用升軌和降軌兩個不同成像視角的LOS向形變信息對東-西向和垂直向的形變進行求解，其求解式為：

(4)

(5)

圖1 升降軌側(cè)視成像幾何圖

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

2.1 研究區(qū)概況

漾濞地震震中距離漾濞縣城約7 km,距離大理市約40 km,其震中附近區(qū)域地質(zhì)構(gòu)造復(fù)雜，且發(fā)生過多次地質(zhì)災(zāi)害。漾濞地震震中位于川滇菱形塊體的西南邊界，是維西-喬后-巍山斷裂與紅河斷裂銜接處，震中區(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相位信息轉(zhuǎn)換為視線向的形變數(shù)據(jù)，得到了降軌和升軌兩個視線向形變場，其中降軌表現(xiàn)為兩個主要形變區(qū)域A和B，A和B大致呈橢圓形狀對稱分布，A表現(xiàn)為遠(yuǎn)離視線向的沉降，最大沉降量約為-7 cm，B區(qū)域表現(xiàn)為靠近雷達視線向的隆升，最大隆升量約為8 cm，這一結(jié)果與王紹俊[10]以及黎朕靈[12]等人的InSAR解算結(jié)果基本一致，降軌同震形變場如圖3所示，在A與B之間形成了一個近北-西走向的沒有明顯形變的長軸，如圖3中紅色虛線所示，有可能為本次地震發(fā)震斷層所在的空間位置，震中附近已知的維西-巍山斷裂位于川滇塊體西部邊緣，北起雪龍山東麓，止于點蒼山西南，走向北北西，但該北-西走向長軸位于維西-巍山斷裂以西，不在其主斷裂上，楊九元[11]等的研究認(rèn)為此次地震可能破裂在維西-巍山斷裂的隱伏分支斷層或一個獨立的未知的隱伏主斷層上，但是準(zhǔn)確的地震斷層還需要進一步研究確認(rèn)。

圖3 降軌同震形變場

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

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

為了更好地分析降軌同震形變場，繪制了aa*、bb*、cc*3條剖面線，其分別位于沉降區(qū)、隆升區(qū)、跨越沉降和隆升區(qū)。3條剖面線的位置已在圖3中標(biāo)明，其中剖面線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內(nèi)隆升量均超過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)星成像幾何的影響，噪聲污染嚴(yán)重，形變特征不明顯，但不難發(fā)現(xiàn)A和B兩個形變區(qū)域明顯擴大，且兩個區(qū)域的形變方向與降軌結(jié)果相反，其中A區(qū)域表現(xiàn)為隆升，最大形變量為7 cm，B區(qū)域表現(xiàn)為沉降，最大形變量為-7 cm，如圖6所示，相同區(qū)域表現(xiàn)出不同的形變方向，與升降軌不同的成像視角有關(guān)。自降軌輔影像成像時間2021年5月22日T23：14起至升軌輔影像成像時間6月1日止，震中附近區(qū)域發(fā)生多次余震，其中M3.0以上共7次(據(jù)中國地震臺網(wǎng))，故升軌解算結(jié)果應(yīng)包含了上述余震貢獻量，形變區(qū)域明顯擴大與余震形變累積有關(guān)。

圖5 剖面線形變曲線

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

圖6 升軌同震形變場

3.2 同震庫侖應(yīng)力變化計算

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

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

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

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

結(jié)合升降軌同震形變場解算結(jié)果，利用式(5)解算得到了東-西向和垂直向同震形變場，規(guī)定東-西向同震形變場形變量正值表示向東運動，負(fù)值表示向西運動，垂直向同震形變場形變量正值表示垂直方向的隆升，負(fù)值表示垂直方向的沉降。兩個方向的形變場呈現(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 結(jié)束語

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