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

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

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

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

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

1 研究方法

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

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

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

(1)

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

1.2 二維形變提取方法

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

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

(2)

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

(3)

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

(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中黑色虛線框?yàn)槲闹型鹦巫兘馑愕膮⒖紖^(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í)驗(yàn)數(shù)據(jù)具體參數(shù)如表1所示。文中除了采用Sentinel-1A數(shù)據(jù)外，還采用歐空局提供的DORIS精密軌道數(shù)據(jù)AUX_POEORB進(jìn)行軌道誤差校正。升降軌影像覆蓋范圍在圖2中用藍(lán)色矩形框表示。

表1 影像信息

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

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

圖3 降軌同震形變場(chǎng)

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

圖4 余震時(shí)間分布及空間分布

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

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

圖5 剖面線形變曲線

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

圖6 升軌同震形變場(chǎng)

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

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

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

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

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

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

圖8 垂直向和東-西向同震形變場(chǎng)

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

4 結(jié)束語(yǔ)

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