興仁縣羅馬山滑坡體形變監測及成因分析

章 彭，劉文明

(1.貴州省水利水電勘測設計研究院，貴陽 550002；2.中國測繪科學研究院 地球觀測與時空信息科學國家測繪地理信息局重點實驗室，北京 100830)

1 引言

作為一種常見的地質災害，山體滑坡分布廣泛，往往極大地危害附近居民的安全和財產。2010—2016年的統計結果表明，山體滑坡占中國所有地質災害的70%以上(全國地質災害通報 2010—2016年)，而我國地質災害主要集中在南部及西南部地區。貴州省作為滑坡災害的主要集中區(呂剛等 2016)，發生于2017年8月的納雍縣普灑鎮塌方事件引起了國內外學者的廣泛關注。

目前已有多種方法被聯合用于滑坡體形變監測，如測地測量，基于全球定位系統的監測網絡，以及遙感圖像的解譯。盡管這些方法獲得了高度精確的結果，但由于人力和儀器成本，它們不容易獲得高密度的測量點。作為有源微波傳感器，合成孔徑雷達(SAR)特別適用于貴州省這種多陰天和雨季地區。利用地形的先驗知識，可以從在不同時間獲取的兩個圖像中檢測毫米變形。所謂的差分干涉SAR(D-InSAR)技術可以大規模地獲得地表形變量，精度可達到厘米級。因此在地震，火山和采礦變形的監測等多種應用中得到了廣泛的應用(Wang et al.2008；廖明生等 2013)。

本文選取貴州省黔西南布依族苗族自治州普安縣羅馬山作為研究區，共收集11景L波段的ALOS PALSAR2數據，生成10個干涉對，以形變時間序列的方式對研究區中的滑坡體進行形變監測分析。相對于X波段及C波段雷達信號，L波段具有較強的穿透能力，在黔西南地區這種多植被覆蓋區域能夠保證較小的時間去相干(田馨等 2013；Dong et al.2018；Zhao et al.2012)。文章最后結合該地區的降雨量信息和形變時間序列，發現降雨量于形變的高度相關性，這對滑坡災害的預防具有重要意義(程海琴等 2014)。

2 研究區基本特征

圖1 不穩定斜坡及礦區位置圖

3 理論及方法

3.1 數據源簡介

本文收集了2017年4月16日至2018年9月16日期間的11景ALOS PALSAR2雷達數據(見表1)，其工作模式為FBS(fine beam single polarization)，L波段(波長23.6 cm)，產品級別為Level1.1(單視復數影像數據)，HH極化，視角為32.8°，方位向像元大小為2.129 888 m，距離向像元大小為1.430 422 m，數據獲取日期信息如表1所示。

表1 實驗區ALOS PALSAR2數據

為了去除地形相位，在進行二軌差分法差分干涉測量時需要利用高精度DEM，本文采用美國地質調查局發布的30 m SRTM DEM數據。

3.2 二軌法D-InSAR技術

二軌法的基本思想是通過引入外部DEM去除干涉相位中的地形相位，實現過程中需要兩景覆蓋同一地區的SAR影像和輔助DEM(趙夢雪等 2017)。首先，需要將DEM與主影像進行配準處理；之后再利用DEM模擬成地形干涉條紋，得到地形相位；從主從影像生成的干涉相位中減去地形相位，即可得到形變相位；最后通過相高轉換，實現相位到視線向形變量的轉換。其公式可以表示為：

(1)

其數據處理的流程如圖2所示。

圖2 雙軌法D-InSAR數據處理流程

(1)影像配準

目前，常用的配準方法有2種：①基于強度和相關系數法的配準方法其包括兩景SAR影像之間距離向和方位向上初始偏移估計以及輔影像重采樣到主影像，分為粗配準和精配準。精配準需要對距離向和方位向的偏移量多項式進行亞像元級的精確估計，并最終將配準精度限制在0.1個像元內。②基于查找表的配準方法。本文采用第一種配準方法。

(2)干涉圖生成和相干性分析

主、輔影像精確配準后進行共軛相乘便得到了干涉圖。在干涉圖生成的同時，以計算相干系數作為評價干涉圖質量的依據，由干涉系數產生的相干圖直觀地顯示了干涉圖的質量。相干系數的計算公式如下：

(2)

其中M，N分別為相干性估計是的窗口大小。相干系數值越大，形變監測值越精確。

(3)干涉基線估計

目前常用的基線估計方法有3種：①采用軌道信息估計基線；②采用快速傅里葉變換(Fast Fourier Transform，FFT)根據干涉圖的條紋變化率估計基線；③基于地面控制點(Ground Control Point，GCP)估計基線。本文先根據方法1進行初始極限估計，然后使用方法2估計殘余基線并將其加入初始估計基線中，以達到基線精化的目的。

(4)相位解纏

相位解纏是恢復雷達影像的原始相位值。目前常用的方法有2種：①枝切法；②最小費用流法(Minimum Cost Flow，MCF)。因為MCF基本不需要人工參與，本文選用后者。

(5)形變圖生成及地理編碼

地理編碼是將相位解纏后的坐標轉化為橢球坐標系下的坐標。這1步的關鍵是確定解纏后影像上每個像元對應的三維坐標。根據相位和形變的對應關系，從解纏后的形變相位轉化為雷達視線向上的形變量。

4 實驗結果

使用二軌法對收集的10景ALOS PALSAR2數據進行差分干涉處理，共生成9組干涉對，干涉對垂直基線及時間基線情況如表2所示。

表2 干涉對基本參數

續表2

編號主影像輔影像時間基線/d垂直基線/m52017102920171126288.0616201711262018021884-149.8647201802182018051384132.0938201805132018062442-178.9589201806242018080542212.464

根據第3.2節的處理流程對9個干涉對進行處理，最終得到9幅視線向形變量圖，利用ENVI和ArcGIS等圖像處理軟件，將試驗區中一處典型滑坡——羅馬山滑坡制作成專題圖，結果如圖3所示。

圖3 羅馬山滑坡形變量圖

需要說明的是，圖中形變量的單位為m，而圖中最大值(0.25)及最小值(-0.25)為人為設置，目的是為了統一圖例。

同時，我們對羅馬山滑坡體上一點的形變歸算到型變速率，并進行了時間序列的統計。該點型變速率與降雨量疊加顯示結果如圖4所示。

圖4 形變速率與降雨量關系圖

5 結 論

從實驗結果中可以看出，滑坡體的形變速率往往伴隨著降雨量的增加而加快，形變速率與降雨量具有高度相關性。也就是說，強降雨是造成滑坡體形變的重要因素之一。同時，形變速率也是滑坡體穩定性評價的一個重要指標，滑坡發生前往往伴隨著滑坡體的形變加速度，使用D-InSAR技術能夠快速高效得獲得滑坡體的形變量，對滑坡災害的防治具有一定的應用價值。