云南漾濞5.21地震誘發地表形變的D-InSAR快速探測分析

朱智富，甘淑，袁希平，張曉倫，黃志豪

(1.昆明理工大學國土資源工程學院，云南 昆明 650093；2.云南省高校高原山地空間信息測繪技術應用工程研究中心，云南 昆明 650093；3.滇西應用技術大學云南省高校山地實景點云數據處理及應用重點實驗室，云南 大理 671006)

1 引 言

根據中國地震臺網(http：//www.ceic.ac.cn/)報道，2021年5月21日，中國云南省大理州漾濞縣(北緯25.67°，東經99.87°)發生MS6.4級地震，震源深度 8 km。本次地震災害發生在夜間，全州十二縣市均有明顯震感，目前漾濞縣受災最嚴重，大量房屋遭到破壞，多處道路中斷，造成了嚴重的經濟損失。據歷史地震目錄記載，漾濞地區處于紅河斷裂帶[1]，屬于地震多發地區，在該斷裂帶的中段，即下關到漾濞這一段，歷史上發生的地震比較多，因此，針對漾濞震后的地表變化進行研究具有現實指導意義。

傳統的地表形變監測方式，由于地形條件和儀器自身局限性的影響，只能獲取到研究區域內有限的地表信息，不能在短時間內，得到大范圍、監測點連續的形變信息，并且這些監測點都是以離散點的形式存在，不能很好地反映出研究區域的地表變化。與之相比，差分合成孔徑雷達干涉測量技術(Differential Synthetic Aperture Radar Interferometry)具有全天時、全天候、大面積和高精度的優勢，能夠監測厘米級甚至毫米級的地表形變，可以為地震監測提供空間連續地表形變信息[2]。目前，D-InSAR技術已經成為監測地震地表形變的重要手段。文獻[3]利用合成孔徑干涉測量對九寨溝震后的滑坡隱患點進行識別與形變監測，共探測出13處滑坡隱患，驗證了基于InSAR技術對滑坡隱患點進行早期識別的可靠性。文獻[4]利用D-InSAR技術對2017年伊朗地震進行反演，獲取了發震斷層的精確幾何參數和最優斷層滑動分布。文獻[5]利用D-InSAR結合GIS技術對礦區開采的沉陷形變進行監測，很好地反映出煤炭在開采過程中的沉降區域位置和分布范圍。文獻[6]利用二軌法D-InSAR技術對北京城區地面沉降進行監測，發現北京市地面沉降較為明顯的兩個沉降中心，表明了D-InSAR技術在城市地區的地面不均勻沉降方面具有明顯優勢。文獻[7]利用升降軌D-InSAR技術對新疆伽師地震的同震形變場進行提取，成功獲取了此次地震的LOS向形變。基于D-InSAR的一系列優勢，本文將利用二軌法D-InSAR技術對漾濞震后的地表形變進行快速探測分析。

2 研究區及數據源

2.1 研究區概況

研究區位于漾濞彝族自治縣如圖1所示，隸屬于云南省大理州，位于大理州中部，東與大理市、巍山彝族回族自治縣毗鄰，西與永平、云龍二縣接壤，南交保山市昌寧縣，北連洱源縣。境跨北緯25°12′～25°54′、東經99°36′～100°07′，屬亞熱帶和溫帶高原季風氣候區，地勢北高南低，最高海拔 4 122 m，最低海拔 1 174 m。2021年5月21日，漾濞縣(北緯25.67°，東經99.87°)發生MS6.4級地震，震源深度 8 km。

圖1 研究區及其裁剪后的DEM數據

2.2 實驗數據

本次實驗數據來自開源的歐空局提供的Sentinel-1A C波段SAR數據(https：//scihub.copernicus.eu/)，模式為干涉測量寬幅模式(IW)，其極化方式為VV極化，距離向和方位向分辨率為 5 m×20 m。實驗數據的具體參數如表1所示。為了提高軌道精度，采用阿拉斯加(https：//s1qc.asf.alaska.edu/)提供的精密AUX_POEORB軌道數據文件對原始的SAR數據進行軌道精煉。采用美國太空總署(NASA)和國防部國家測繪局(NIMA)聯合測量的空間分辨率為 30 m的STRM1-DEM數據，作為參考高程數據，用于地形相位的去除。

Sentinel-1A干涉參數 表1

3 研究方法與技術

3.1 D-InSAR探測技術

InSAR技術是利用雷達系統獲取同一地區兩幅SAR影像所提供的相位信息進行干涉處理，來獲取地表的形變[8]。D-InSAR是在傳統InSAR基礎上發展起來，利用同一地區不同時刻獲取SAR數據進行差分干涉處理，通過差分處理去除兩次觀測相位中的平地效應、地形相位和大氣延遲等，得到形變相位，進而獲得地表形變信息[9]。目前的D-InSAR方法包括了二軌法，三軌法和四軌法，這主要是由所選影像的數量來決定的。綜合對比不同方法性能差異，二軌法干涉方法簡單而最為可靠，因此在本文實證研究中采用了二軌法。

研究采用的二軌法D-InSAR技術的基本原理及流程如圖2所示。在形變監測中，二軌法D-InSAR的主要技術核心是使用跨越形變期的兩幅SAR影像和該區域的外部DEM來完成差分干涉[10]。

圖2 二軌法D-InSAR原理示意圖

將兩幅SAR影像配準重采樣并共軛相乘形成一幅干涉圖，干涉圖中每個像元的涉相位包含該點相對于參考位置的高程和形變信息。考慮到大氣環境和地表的影響，干涉圖中的相位φint可以表示為：

φint=φref+φtopo+φdef+φatm+φnoise

(1)

在(1)式中，φref表示參考橢球相位，φtopo表示地形相位，φdef表示形變相位，φdef表示差分干涉相位，φatm表示為大氣相位，φnoise表示噪聲相位[11]。

特別需要說明的是，由于后面兩項的大氣和噪聲的影響是無法消除的，所以這里可以不再考慮大氣和噪聲的影響，從而相位φint可簡化為：

φint=φref+φtopo+φdef

(2)

通過軌道數據和外部DEM去除參考橢球相位和地形相位后，便可得最終的形變相位φdef，表示為：

(3)

其中λ為波長，△r為傳感器斜距方向的形變量。

3.2 地震前后D-InSAR數據處理

本次實驗采用了二軌法D-InSAR差分技術，利用SARscape軟件對震前和震后的SAR數據進行處理，其基本的處理技術流程如圖3所示，包括：①原始影像的預處理；②基線估算；③干涉相位圖生成；④濾波與相干性計算；⑤相位解纏；⑥軌道精煉與重去平；⑦相位轉形變及地理編碼。

圖3 D-InSAR處理流程

試驗研究中的具體操作步驟如下：

(1)原始影像的預處

每景Sentinel-1影像數據覆蓋范圍為 400 km×400 km，而研究區域遠遠小于一景影像的覆蓋范圍，為了節省時間和提高處理數據的工作效率，首先對原始影像進行裁剪，得到震前和震后的研究區影像。

(2)基線估算

這一步主要是估算主從數據的基線情況，包括時間基線、空間基線、多普勒偏移、一個相位周期變化代表的高程變化等信息[12]，其中空間基線的垂直分量不能超過臨界基線的閾值，否則就會丟失相位信息，從而無法進行干涉。因此，垂直基線是判斷所獲取的SAR數據能否做干涉的重要依據。

對裁剪后的研究區進行基線估算，得到空間基線為 49.061 m，遠遠小于臨界基線 6 262.828 m，滿足D-InSAR處理要求。

(3)干涉相位圖生成

干涉相位圖生成主要是基于軌道信息將兩景影像進行配準，然后用SLC數據所對應的像元進行復數共軛相乘，所得到的干涉圖像(圖4)包含了平地效應的影響，接下來引入外部DEM數據去除平地效應[13]，得到去平后的干涉圖像(圖5)。

圖4 干涉圖 圖5 去平后的干涉圖

(4)濾波與相干性計算

經過干涉處理后的圖像中還包含形變相位、大氣延遲相位和噪聲相位所帶來的影響，這些影響在差分干涉圖中表現為周期不明顯，條紋不清晰和斑點噪聲，為后續的相位解纏帶來困難[14]。通過相位濾波，削弱噪聲信號，提高相位解纏的準確性。常用的濾波方法為Goldstein自適應濾波，該方法可以在保證條紋清晰度的情況下有效的濾出噪聲，經過SARscape軟件處理后得到濾波后的干涉圖像(圖6)和相干性圖(圖7)。相干性系數的范圍為0-1，相干性系數越大，相干性越高，解纏效果也越好，反之亦然。

圖6 濾波后的干涉圖 圖7 相干性圖

(5)相位解纏

差分干涉圖中的形變相位被纏繞在[-π，π)區間內，不能直接反映出真實的地表形變信息。經過相位解纏處理可以恢復每一像素的相位整周數，得到連續的地形變化[15]。本實驗采用了SARscape軟件提供的Minimum Cost Flow(最小費用流法)進行相位解纏處理，得到解纏后的干涉圖像(圖8)。該方法可以考慮圖像上所有的像元，對相干性小于閾值的像元做掩膜處理。

圖8 解纏后的干涉圖

(6)軌道精煉與重去平

以濾波后和解纏后的干涉圖為參考，選取GCP地面控制點。利用選取的GCP點進行軌道精煉和相位偏移計算，消除可能的相位坡道和殘余的軌道相位，此時就只剩下形變的相位。在GCP控制點的選取過程中要遠離干涉條紋密集和解纏錯誤的區域，選擇相干性較高和平地區域。本實驗選取了11個GCP點，用Automatic Refinement(自動優化法)進行軌道精煉與重去平。

(7)相位轉形變及地理編碼

把經過絕對校準和解纏的相位，結合成合成相位，將僅包含形變信息的差分干涉相位轉為雷達視線方向的形變量[16]，然后通過地理編碼，把SAR坐標系下的形變圖轉化到WGS-84坐標下，最終得到LOS向的地理坐標形變圖(圖9)。

圖9 LOS方向上的地表形變圖

4 D-InSAR形變監測結果分析

4.1 形變監測結果的沉降凸顯斑塊識別分析

根據圖9的形變空間分異所示，首先可以看出此次地震造成4處形變比較明顯的區域，且4個區域呈現出以EE′為分界線的對稱分布，其中分界線EE′的西側區域的形變量為負值(遠離衛星方向)，表現為沉降，而東側區域的形變量為正值(靠近衛星方向)，表現為抬升。實驗中為了進行具體對比分析，對4個形變量較大集中區進行多邊形圈定如圖10所示，并分別命名為區域A、B、C、D，即得出區域A和區域B的最大沉降值為 8.0 cm，區域C和區域D的最大抬升值為 8.8 cm。區域A沉降區域面積較大，但沉降數值相對較小，而區域B雖然沉降區域面積較小，但它的沉降值相對較大，二者在沉降面積和沉降數值上形成鮮明對比。區域C和區域D的抬升面積和抬升數值基本一致，其中抬升面積較沉降區域面積而言，相對較小，說明此次地震造成的地表形變主要還是以沉降為主。最后結合4個形變區與最大震中位置的空間相關關系進行初步解讀分析，表明此次地震造成地表形變最大的區域并不是出現在最大震級發生的震源中心位置，而是發生在偏離最大震級震中的東南側，即形變區A、B、C、D4個區域。

圖10 典型形變區的多邊形圈定示意圖

4.2 形變監測結果與震源空間關聯分析

為了進一步對利用D-InSAR地震監測結果的形變空間分異特性開展深化探討，研究中收集整理了對中國地震臺網發布的關于漾濞縣在2021年5月21日～5月22日的4級以上地震記錄資料如表2。初步統計分析可知，在此期間，漾濞縣一共發生了14次4級以上的地震，而將這14次地震點的位置信息與D-InSAR監測的沉降形變圖進行空間疊加可得到圖11。分析圖示空間關聯可以發現：記錄的14次4級以上的余震中心點位置，有10次分布于沉降A區和B區，有2次分布于抬升C區，亦即有12/14次的頻度的4級以上的地震中心分布與利用D-InSAR技術探測出的地表形變較大區域吻合，該結果一定程度上較好地實證了利用D-InSAR技術進行定量、定位而連續探測地震誘發地表形變空間分異的可行性及可靠性。

4級以上地震記錄 表2

圖11 4級余震分布圖

4.3 典型形變監測剖面分析

最后，為了更加深入地揭示D-InSAR形變監測結果空間分異特性，進行不同方向的形變大小量化認識分析，針對D-InSAR形變監測結果空間圖層，研究中特別設置了與分界線EE′近似平行的，分別穿過沉降B區和抬升C區的兩條剖面線BB′和CC′，以及一條與分界線EE′近似垂直且穿過抬升D區中心的剖面線B′C′如圖12(a)所示，并繪制出對應剖面線的形變剖面圖(圖12(b)、(c)、(d))。根據三條典型形變剖面曲線圖的解析可以看出：

圖12 典型形變區剖面線圖

(1)剖面線BB′橫跨沉降區域A和B，主要是以沉降為主，剖面線整體呈現凹形，局部呈起伏變化狀態。主要集中在 15 km～20 km、22 km～25 km和26 km～28 km之間發生了3次劇烈沉降，形成3個較大的“漏斗型”沉降區，其中沉降的最大值達到 5.5 cm左右。

(2)剖面線CC′主要位于抬升區C，剖面線整體呈現凸形，局部呈起伏變化狀態，在剖面線 0 km～6 km之間變化迅速，抬升量達到了 7 cm左右，之后又相對變得平緩，在 11 km～15 km之間雖然曲線表現為下降趨勢，但形變值仍然沒有出現負值，說明該區域的主體形變仍然還是以抬升為主。

(3)剖面線B′C′橫跨抬升和沉降區，呈正負分布，在對稱線EE′的西側表現為沉降，東側表現為抬升，且抬升趨勢遠遠大于沉降趨勢，大約在 13 km位置剛好為沉降與抬升的分界點，可以大概推測它為分界線EE′的位置，剖面線B′C′經過分界線EE′從沉降區B進入抬升區C。剖面線在分界點處沒有出現斷裂現象，說明此次地震的破裂沒有到達地表。在 13 km～18 km之間上升趨勢劇烈，達到了抬升的峰值，之后的變化相對比較連續平穩，基本沒有“跳點”出現。

5 結 論

本文利用二軌法D-InSAR技術對5.21漾濞地震誘發的山地形變空間分異特征進行快速探測，識別出4個地表形變較大的區域，并且發現地表形變變化最大的位置并不是在最大震級的地震中心，而是分布在最大震級震中東南側的次級地震頻發的空間區域。結合中國地震臺網發布的大于4級余震發生的震中位置記錄資料整理，空間疊加D-InSAR形變探測結果進行初步關聯分析，發現這些點位與D-InSAR探測出的形變集中區域具有較大的空間關聯性，一定程度上實證了D-InSAR技術用于震后地表形變快速探測的實用性及可行性。最后，基于D-InSAR監測出的地震后形變空間分異圖及所識別出的4個地表形變集中地塊，通過三個典型剖面圖繪制及曲線特征發現，初步分析認為此次地震造成的主要地表形變變化以沉降為主。