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

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

1 引 言

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

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

2 研究區(qū)及數(shù)據(jù)源

2.1 研究區(qū)概況

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

圖1 研究區(qū)及其裁剪后的DEM數(shù)據(jù)

2.2 實驗數(shù)據(jù)

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

Sentinel-1A干涉參數(shù) 表1

3 研究方法與技術

3.1 D-InSAR探測技術

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

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

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

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

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

(1)

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

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

φint=φref+φtopo+φdef

(2)

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

(3)

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

3.2 地震前后D-InSAR數(shù)據(jù)處理

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

圖3 D-InSAR處理流程

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

(1)原始影像的預處

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

(2)基線估算

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

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

(3)干涉相位圖生成

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

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

(4)濾波與相干性計算

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

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

(5)相位解纏

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

圖8 解纏后的干涉圖

(6)軌道精煉與重去平

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

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

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

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

4 D-InSAR形變監(jiān)測結果分析

4.1 形變監(jiān)測結果的沉降凸顯斑塊識別分析

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

圖10 典型形變區(qū)的多邊形圈定示意圖

4.2 形變監(jiān)測結果與震源空間關聯(lián)分析

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

4級以上地震記錄 表2

圖11 4級余震分布圖

4.3 典型形變監(jiān)測剖面分析

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

圖12 典型形變區(qū)剖面線圖

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

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

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

5 結 論

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