基于短基線InSAR技術的滑坡監測研究
——以四川茂縣為例

彭藝，吳云龍，楊玉，袁華清

(1.中國地震局地震研究所，中國地震局地震大地測量重點實驗室，湖北 武漢 430071；2.防災科技學院，河北省地震動力學重點實驗室，河北 三河 065201)

1 引 言

我國疆土遼闊、物質資源豐富，全國各地存在著各種各樣的地質條件，高原、山地、盆地、林地等各種地形遍布全國。這些各不相同的地質背景，加上各地不同的氣候影響，致使各種地質災害的發生，這大大影響著我國的民生建設和經濟發展。而其中，地面沉降作為一個典型的環境地質災害，不僅反映了沉降區域的地表變化趨勢，也表征了相應區域的地質變化規律，同時也可以在此基礎上探究其他災害，如地震、泥石流等的成因和發展趨勢。因此，合理利用現有的科學監測方法，選定較為典型的研究區域進行監察和預測，同時對監測結果進行分析并制訂相應的應對措施，是一項十分重要的研究。傳統滑坡沉降的監測主要依賴于人工勘測例如三角高程測量，即使是使用了GNSS等手段獲取數據，仍然存在著小范圍、低精度、高成本以及分辨率低和數據處理周期長等各種缺陷[1]。隨著衛星遙感領域的進步和發展，基于遙感手段的滑坡沉降監測越來越成為一種重要的方法。

合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar，InSAR)[2]技術是一種星載測量技術，相比于傳統的測量方式，將Sentinel-1衛星影像數據經相干處理得到沉降區域的地表形變信息，從而直觀地表征研究測區的地表沉降情況，具有全天候、全天時、覆蓋廣泛、分辨率高和精度高的優點，目前在山體滑坡、三角洲地區沉降[2]、地震火山、礦區沉降以及城市地下水開采引發的地面沉降監測等各方面已經得到廣泛應用。長安大學的趙延嶺利用干涉圖堆疊技術和SBAS技術對長江三峽區域的滑坡進行了探測與識別，獲取形變信息，分析滑坡成因，做出合理預測，整體效果比較理想[3]。InSAR技術主要分為兩種：傳統的差分干涉測量方法(D-InSAR)和時間序列差分干涉測量，其中D-InSAR技術的構建原理是對覆蓋同一地區的兩幅干涉圖像(一幅形變前的干涉圖像，另一幅形變后的干涉圖像)進行差分處理以獲取目標區域地表形變[4，5]。在早期的地表形變監測中，D-InSAR技術發揮著非常重要的作用，但是隨著科學技術的發展，其對時間和基線去相干過程中產生的誤差、干涉圖噪音、大氣效應和外部DEM誤差的消除越來越無法達到精度[6]，為了更好地提高影像的相干性，在此基礎上提出SBAS-InSAR小基線技術，選擇最短基線對影像對進行配準，提高了數據處理的精度，同時小基線技術對于DEM的數據精度要求相較于其他InSAR手段并不高，是目前應用較為廣泛的一種地表形變監測手段[7，8]。

為研究SBAS-InSAR技術在地面沉降監測中的應用，本文以四川省茂縣為例，利用SBAS-InSAR技術研究該地山體滑坡的發生以及對該地地形產生的影響，進而結合當地地質氣候條件推測我國西南地區的地形變化，再而分析我國全國范圍內山體滑坡的趨勢和影響，對經濟文化發展的影響也可以因地制宜地相應得到預測和解決。

2 研究區域概況

四川省茂縣是四川省阿壩藏族羌族自治州的一個下轄縣，該縣處于四川省的西北部，地質上處于青藏高原和川西平原的交界地，地勢變化大，總體呈現西北高東南低的趨勢。茂縣地貌以高山峽谷為主，海拔高低懸殊比較大，氣候受到西風環境和印度洋西南季風影響，局部氣候較為復雜，多種氣候狀況都有。常見各種氣候災害，再加上位于岷江上游，水資源豐富，春夏季節經常會有暴雨、冰雹、泥石流等自然災害發生，同時也位于地震多發區，極端天氣相對較多，出現大型滑坡的情況也相對較多。

除當地的地理氣候環境容易導致山體滑坡等的自然災害的發生，茂縣的地質條件特殊，這也決定著當地地質災害的發生。根據青藏高原的整體地勢走向，自西向東將青藏高原東緣分為青藏高原地貌區、龍門山高山地貌區和四川盆地地貌區[9]。而導致這些地貌產生的原因是地下的地質構造，青藏高原東緣地質構造自西向東分別是松潘-甘孜褶皺帶、龍門山造山帶和四川盆地[10]。本文選擇的四川省茂縣位于青藏高原地貌區和龍門山高山地貌區的交界處，即松潘-甘孜褶皺帶和龍門山造山帶的交界[11]，周圍發生過多次較大的歷史地震。地質活動活躍，地形起伏大，再加上岷江上游水源的沖刷，這里形成陡峭的峽谷地勢。此外，茂縣的地層巖石分布著砂巖、板巖、泥頁巖和碳酸鹽等，土質易變質變形，在地下構造變化的過程中極其容易受到受影響，巖石易碎，容易產生地面沉降，如果遇到雨水沖刷，山體滑坡和泥石流無法避免[11，12]。再加上地勢崎嶇不平、交通不便，在自然災害出現后經常無法第一時間進行救援，而這里土地、礦產和生物資源極其豐富，歷史上也是兵家必爭的要塞之地，所以四川茂縣的山體滑坡現象是我們必須時刻監測預警并制定合理措施的一項工作。

歷史上，四川茂縣發生過兩次十分嚴重的山體滑坡災害，對這里造成了非常大的影響。1933年8月25日，四川省茂縣疊溪發生的7.5級強烈大地震誘發山體滑坡，幾乎將這個位于四川省西部的羌族聚居區全部掩埋，只剩一座破敗不堪的城隍廟，據不完全統計，死亡人數約有 2 500人，造成嚴重損失。2017年6月24日，茂縣疊溪鎮新磨村的山體突發高位垮塌，被埋100多人，還造成了 2 km的堵塞河道，1 600余米的道路掩埋，以及光纜受損有 3 km長的損害，基本中斷基礎通信。經初步分析，這是一起由強降雨誘發的山體滑坡，再加上茂縣地處地震斷裂帶以及汶川地震對山體的影響，在連日降雨的情況下發生了山體滑坡。

3 研究方法

針對四川茂縣的山體滑坡監測，本文采用合成孔徑雷達(Interferometric Synthetic Aperture Radar，InSAR)的技術進行測量，相對于一些需要人工進行操作的三角高程測量、GNSS測量等，InSAR技術具有全天候、全天時、覆蓋廣泛、分辨率高和精度高的優點，是目前對山體滑坡進行監測精度最高也最合理的方法。在2002年，Berardino等人在InSAR技術的基礎上提出了小基線技術，獲取地表形變以研究大尺度上的形變[8]，后經Hopper等人逐步優化，逐步提高了傳統InSAR技術時間采樣率，減弱空間失相干造成的影響和地形對差分的影響，增強了形變監測的可靠性，在地表形變監測中發揮著非常重要的作用。

本文選取了2017年3月15日～6月19日期間四川茂縣以及周邊區域的8景Sentinel-1A數據，利用Hooper[13]提出的SBAS-InSAR數據處理模型進行處理，對影像數據進行匹配，然后進行干涉處理、濾波處理和解纏，得到最小基線構成的三角網后，再進行反演得到地表形變的速率，最后提取時間序列曲線，并進行時序分析，找到滑坡和疑似滑坡的區域。

3.1 算法模型

在進行反演之前，要先對影像的相位進行計算。根據衛星影像獲取的N+1幅SAR影像，形成M幅干涉圖，生成干涉相位，再利用外部數據去除地形相位影響，根據最小基線原則，利用三維相位解纏的方法進行解纏，之后進行濾波處理以消除各種影響。

在相應時間間隔內，使用同一傳感器獲取覆蓋同一地區的N+1幅影像，根據集合內基線距最小、集合間基線距最大的原則，通過給定閾值，可生成干涉圖M幅(M和N滿足(N+1)/2≤M≤[N(N+1)]/2)，第i幅干涉圖任意像元(x，y)處得到的干涉相位可表示為：

фobs，i(x，y)=фtopo，i(x，y)+фdefo，i(x，y)+фatm，i(x，y)+фnoise，i(x，y)

(1)

其中фobs表示纏繞相位，фtopo表示地形相位，фdefo表示DEM誤差引起的相位，фatm表示大氣噪聲引起的相位，фnoise表示一起熱噪聲等其他相位，這些相位都可以通過濾波來去除。

之后再經過圖像匹配、干涉處理、去除地形相位、解纏計算、濾波處理后進行最小二乘計算，得到地表形變速率。其中最小二乘計算的公式如下：

AVlos=φobs

(2)

Vlos=A+φobs

(3)

(4)

提取形變時序后，在各相位點上聯立方程組，利用最小二乘進行解算。利用式(2)求矩陣A+，使用的是矩陣的奇異值分解法(SVD)求其廣義逆，即A=USVT。再利用形變量的差和時間差得到各個時間點的形變速率矩陣即式(4)，帶入式(3)可得到形變速率的矩陣[1，2]，即通過A+=VS+UT求得A+，最終求得Vlos[6，8，13，14，15]。

3.2 計算流程

本文使用SBAS-InSAR模型，對衛星影像圖完成干涉處理、相位解纏和濾波處理，最后得到形變速率并分析其形變規律[16]。即先選擇的主影像，將副影像與其進行配準，并進行干涉處理，參照DEM高程數據去除地形相位影響，根據小基線原則，生成小基線集干涉圖，設定閾值后篩選SDFP候選點，通過估算其穩定性，選取穩定性高的點以去除估算誤差，在將此時得到的形變差利用三維相位解纏的方法進行解纏；最后進行濾波處理去除噪聲，提取形變時序，利用最小二乘法，同時利用矩陣的奇異值分解法建立方程組運算。具體流程如圖1所示：

圖1 SBAS處理數據流程圖

4 實驗結果

本文選取2017年3月15日～6月19日期間茂縣8幅Sentinel-1A影像數據，根據茂縣地理位置103°39′46″E，32°4′47″N，裁剪至合適大小，基于SBAS-InSAR模型，將各影像兩兩連接匹配，生成連接圖和時間基線，并自動定義超級主影像為20170315，再將各像對之間的匹配之后會對相距較遠的兩個數據進行解纏，得到最小基線的影像對配準三角網(圖2)。

圖2 影像對連接圖、時間基線分布圖、解纏后三角網

將解纏后的影像對進行干涉處理、濾波處理和3D解纏。參照拼接好的茂縣區域的DEM高程數據，選取閾值為0.2進行DEM法干涉處理，同時選擇高斯濾波去除噪聲提高影像對相干性，并使用最小費用法進行3D解纏。本文選擇相干性較好的影像對20170315-20170327的相干系數圖、干涉圖、解纏后影像圖進行展示(圖3)，并剔除相干性不好的影像對。

圖3 20170408-20170514相干系數圖、干涉圖、解纏后影像圖

選取相干性較好的影像對20170408-20170514，選取45個GCP點進行軌道精化和重去平以提高影像對的配準精度，選取GCP點時遵循盡量均勻分布，在未發生形變區域處人工選取的原則。經兩次反演，分別進行二次解纏，以達到優化干涉像對、去除DEM殘余誤差的目的，以此反演出地表形變速率和殘余地形，以及使用大氣濾波以去除大氣信息、估算形變組分。將兩次反演的SAR坐標系下地表信息轉換至WGS-84坐標系下，通過內插和平滑的方法進行地理編碼，得出地表形變平均速率，經密度分割后選取測區特征點繪制時間序列曲線分析其走勢，并結合地質氣候條件進行分析。

5 結果分析

將測區的地表形變速率進行密度分割(圖4)，該區域最大抬升速率大約為 50 mm/a，即圖中紅色的部分，而且分布在這里的點非常少，說明出現抬升的地方相對來說比較少。同時也可以看出，最大沉降速率可達到 -100 mm/a，如圖4所示藍色部分分布的點較多，說明這段時間內出現了比較嚴重的滑坡。

圖4 形變速率分布圖

同時選取4個地勢較高地特征點和1個地勢較低地特征點繪制時間序列曲線(圖5，圖6)，分析監測區域的地面沉降狀況。5個特征點的地表在3月～4月初整體形變十分微小，基本沒有變化，甚至會有微小的地表抬升，但是從4月～5月初，出現了明顯的高程值變化，特征點處呈現較陡的曲線走向，說此處發生了一定程度的山體滑坡，不過在5月中旬，5個特征點的高程變化值又集中在0的上下，說明5月中旬這里沒有發生大的地勢形變事件，即沒有影響較大的自然災害發生，而在5月末～6月中旬這段時間，時間序列曲線卻出現了驟變，如圖前4個特征點的時間序列曲線急速下降為負值，第5個點的時間序列曲線迅速上升，并且按照曲線的走向，在6月19日之后仍會按照這個趨勢下降和上升，說明此段時間內此地發生了較為嚴重的山體滑坡，并且隨著時間的增加，山體滑坡會繼續發生，而且甚至可能會有更加嚴重的自然災害發生，可能會對附近居民的生活以及道路產生十分嚴重的影響。

圖5 特征點1-4時間序列曲線

圖6 特征點5時間序列曲線

四川茂縣的山體形變發生上述地勢變化，與其所處的地勢條件和該處的氣候變化有著非常大的關系。這里地處四川盆地西北部，山地較多，四周分布著海拔為 1 000 m～3 000 m的青藏高原、云貴高原、巫山、大巴山和大婁山等幾座高山，周圍群山環繞，西北邊緣是很長的龍門山脈，中間地勢低矮，可明顯分為邊緣山地和中間盆地兩部分。同時，四川盆地也是我國著名的紅層盆地，其地表巖石主要是由紫紅色砂巖和頁巖組成，它們極易風化發育成紫色土，土壤中含有豐富的鈣、磷、鉀等營養元素，土壤極其肥沃，但是土層淺薄，母巖疏松，非常易于崩解，很容易在大雨天氣隨著雨水被沖走造成泥石流，嚴重的甚至可以造成山體滑坡。除此之外，這里屬于亞熱帶季風性濕潤氣候，夏季平均溫度24℃～28℃，極端高溫達到36℃～42℃，年平均降水量在 1 000 mm～1 300 mm之間，但是年內降水分布不均，整體氣候呈現冬干、春旱、夏澇、秋綿雨的分布趨勢，70%～75%的雨量都集中在6月～10月的夏季，高溫多雨的夏季很容易發生連續暴雨天氣，由此造成山體滑坡的大型自然災害現象。

所以正是由于四川茂縣處于這樣的地形和氣候條件中，才會出現如圖所示的時間序列的變化。圖中3月份5個特征點的高程變化值基本都集中在0的上下，會有微小的形變，因為四川盆地的地表巖石易于風化，而此時正值春天，相對來說，風比較大，再加上紫色土土質疏松，容易崩解，巖石的風化以及土壤的流失都是這里地勢發生變化的原因。根據2017年全年茂縣月降水量(圖7)分析，4月份開始進入夏天，雨水天氣逐漸增多，處在山脈地區再加上土壤疏松，泥土隨著雨水就從地勢高的地方流失到地勢較低的地方，也就產生了圖中所示的前4個特征點高程變化值是負值，而第5個特征點的高程變化值則為正的情況，且明顯從圖中分析得到，第5個特征點的地勢相較于前4個都比較低。5月中旬，根據圖中茂縣平均氣溫分布圖以及降水量分布圖，可以發現5月份茂縣的平均溫度并不算高，而且降水量也相對較少，因為5月17日～18日這段時間，茂縣遭受了大風、冰雹自然災害的襲擊，所以地勢沒有大規模的變化。

圖7 2017年茂縣月降水量和平均溫度圖

5月中旬的冰雹等極端天氣結束之后，正式進入夏季，在亞熱帶季風性濕潤氣候的影響下，氣溫持續上升，降水量也開始持續增加。此后降水集中，強降雨天氣持續，在雨水的沖刷下，疏松的土壤開始發生滑坡，所以5月中旬以后，地勢高的4個特征點的高程變化值發生驟變且為負值，而地勢較高的第5個特征點則發生了抬升。因為當地地勢變化大，歷史地震致使茂縣巖體結構發生破壞和裂縫，在持續的強降雨的沖刷下，具有單面山斜坡結構的后緣巖層被拉斷[12]，發生了順層滑動，誘發持續的泥土流失。截止到本次試驗研究的6月19日，特征點持續發生沉降和抬升，并且按照時間序列曲線，在未來的時間，該地區還會按照這個趨勢繼續沉降和抬升，而且曲線變化很陡，很有可能發生更為嚴重的自然災害，如大型山體滑坡。而事實也確實如此，2017年6月24日6時左右，茂縣疊溪鎮新磨村山體突發了高位垮塌，道路和村民被掩埋，是一起十分嚴重的自然災害事件，對當地的居民造成了十分嚴重的人身和財產傷害。

6 結 論

經過上述的實驗過程，對2017年3月～6月之間四川茂縣的衛星影像數據，基于小基線技術，對衛星影像進行了干涉，濾波和解纏等，再進行兩次反演得到其形變速率，最后對其進行地理編碼，找到特征點的時間序列，通過不同時間段的高程值變化來研究該處的形變情況，并同時結合當地全年月平均降水和氣溫數據分析其成因，最終得到的結論有：

(1)四川省茂縣在2017年3月15日～6月19日期間，地表發生了一定程度的沉降，最大沉降速率和抬升速率分別達到 -100 mm/a和 50 mm/a，且地面沉降隨著降水和溫度的變化而改變。

(2)3月茂縣附近整體地表形變不大，直至5月中旬，出現一定幅度的沉降，到5月中旬出現停止沉降的現象，之后到6月，沉降急速增加，并且有按照這個速度持續沉降的趨勢。

(3)結合2017年茂縣月平均降水和溫度數據，以及當地地質條件和構造，分析得出：處在地質構造帶交界處的茂縣地質條件特殊，巖體結構已遭受破壞，持續的降水容易誘發山體滑坡的發生。

以上就是基于小基線技術對3月～6月四川茂縣衛星影像數據的形變速率相干性分析，并結合當年月平均溫度和降水數據以及當地的地質構造背景分析，得出了此處可能會在將來發生大型山體滑坡的結論，整體吻合當地具體滑坡變化。依據SBAS-InSAR技術對地質活動活躍的地區進行定期監測，并結合當地地質條件和降水數據做出分析和預測，方法合理且結果可信，對城市的發展和規避一些自然災害帶來的損失具有重要意義。