利用Sentinel-1A的時序干涉測量探測蒲縣滑坡形變

高二濤， 李 豪， 雍 琦， 付波霖， 范冬林

(1.桂林理工大學測繪地理信息學院，桂林 541006； 2.廣西空間信息與測繪重點實驗室，桂林 541006；3.山西省工業設備安裝集團有限公司，太原 030000)

中國地質環境復雜、氣候條件時空差異大，地震、滑坡、泥石流等地質災害類型多且分布廣，難以有效防范[1-2]。在上述自然災害中，滑坡是一種危害程度僅次于地震，排名第二的環境地質災害。滑坡廣泛發育于河流、湖泊、海洋和溝渠的岸坡，地形高差大的山區、峽谷區，道路的斜坡地段及強降雨地區等，具有隱蔽性、突發性、不確定性等特點，且破壞力巨大[3]。蒲縣屬于晉西黃土高原地質災害區，地貌復雜多樣，此地區應重點加強對居民居住區和采煤區的黃土滑坡、崩塌等地質災害的防范。滑坡在突發前一般伴隨長期緩慢的地表變形，在此背景下，基于變形監測為主的滑坡災害早期識別與持續監測可以為預報預警與防災減災工作提供科學依據，對于主動防范地質災害意義重大[4]。

水準測量、全球導航衛星系統(global navigation satellite system，GNSS)觀測、攝影測量法和位移計等傳統滑坡監測方法以其點觀測、高精度等特點，目前已被廣泛應用于滑坡監測中[5-7]。近年來，空間對地觀測技術迅速發展，以永久散射體雷達差分干涉測量技術(persistent scatterer interferometric synthetic aperture radar，PS-InSAR)與短基線集雷達干涉技術(small baseline subset，SBAS)為代表的時間序列合成孔徑雷達干涉(time series interferometric synthetic aperture radar，TS-InSAR)技術，可以遠程監測地球表面以獲取全球的高精度地形和地表微小形變信息，具有全天時、全天候、時空分辨率高、覆蓋范圍廣等特點，彌補了傳統滑坡形變監測方法計算量大且監測成本高的不足，在滑坡、地震、城市地面沉降等地質災害的普查與監測預警中展現出巨大的發展潛力[8-12]。現采用PS-InSAR與SBAS技術分別處理覆蓋山西蒲縣地區24景哨兵-1A(Sentinel-1A)數據，獲取研究區域2019年年平均形變速率和累計形變位移量，通過結合已有的資料和谷歌地球(Google Earth)等數據，進一步分析是否為滑坡隱患點及其可能性成因。

1 研究區概況與數據簡介

1.1 研究區概況

蒲縣(36°11′32″N～36°38′13″N，110°51′09″E～111°23′36″E)位于山西省西南部，臨汾盆地中心，呂梁山西南側，東部屬于呂梁山和太行山斷裂帶，西部屬于鄂爾多斯斷裂帶。蒲縣東接洪洞，西靠大寧，南依吉縣，北與隰縣、汾西毗鄰。境內由山地、丘陵、河谷等地形組成，地勢剖面為東高西低，北、東、南三面環山，植被茂盛。西部與中部為典型的黃土高原地貌，第四紀黃土發育良好，地表大多被不同厚度的黃土覆蓋，地形破碎、垂直節理發育、溝谷深切超過200 m，易引發崩塌、滑坡、泥石流等多種次生地質災害。蒲縣是山西煤炭開采的重要基地之一，含煤面積占全縣總面積的九成以上[13]。蒲縣位于北半球中緯度大陸內部，是典型的溫帶大陸性氣候。春季多風沙，降雨量少；夏季多暴雨，雨熱同期；秋季短暫、溫和涼爽；冬季漫長、寒冷干燥。蒲縣地理位置如圖1所示。

圖1 蒲縣地理位置

1.2 實驗數據

哨兵1號(Sentinel-1)衛星是歐洲航天局(ESA)發射的載有C波段合成孔徑雷達的地球監測衛星[14]，可全天時、全天候提供連續圖像。哨兵1號是哥白尼計劃的一個重要組成部分，具有1A和1B雙系統。2014年4月發射的Sentinel-1A衛星是太陽同步極地軌道衛星，衛星的重訪周期為12 d，軌道高度是693 km。它的主要用途是監測地球環境變化，存在四種成像模式、兩種極化方式，具有生產能力強、重訪周期短的優勢。現采用Sentinel-1A衛星獲取的24景干涉寬幅模式(interferometric wide swath，IW)的單視復數圖像數據(single look complex，SLC)。獲取影像的時間范圍是2019-01—2019-12，數據的詳細參數如表1所示。

表1 Sentinel-1A TOPS SAR影像數據參數

本文研究采用日本宇宙航空研究開發機構獲取的30 m分辨率的AW3D30(ALOS World 3D-30 m)作為外部DEM(digital elevation model)數據。DEM的范圍在35°N～38°N、108°E～112°E，運用其生成地質背景構造圖和去地形相位處理，以消除地形起伏的影響[15]。

2 時序InSAR數據處理

2.1 基于PS-InSAR技術的蒲縣滑坡形變監測

首先，利用PS-InSAR技術監測蒲縣地區滑坡形變過程，主要處理步驟如圖2所示。實驗程序自動選取2019-06-30的影像數據為公共主影像，其余23景從影像逐個與之配準并建立連接關系，共得到23個干涉像對。SAR數據時空基線分布如圖3所示，具體時空基線如表2所示。可知，時間基線范圍為-180～168 d，空間基線范圍為-135～110 m。

表2 Sentinel-1A影像數據獲取時間及基線參數信息

圖2 PS-InSAR處理流程

圖3 PS-InSAR時空基線分布

使用AW3D30 DEM對干涉相位進行去平處理，實驗采用的DEM精度和分辨率越高，去除地形相位的效果越理想。打開Goldstein濾波器，將距離向多視數設置為5，方位向多視數設置為1，輸出差分干涉圖。選擇永久散射體(PS點)時，應選取有穩定散射特性的點，如橋墩、大壩、房頂等人造地物或者裸巖等自然地物。通過密集分布的PS點消除由于信號在大氣中傳播延遲造成的波動，從而有效改善實驗精度。為提高穩定散射體點的準確性，采用振幅離差指數法確定PS點，將相干系數閾值設置為0.75，對研究區域內的穩定散射目標進行有效識別[16]，如圖4所示，在圖4中用綠點標注，PS點目標主要分布在裸露的巖石、公路上，而植被覆蓋率高的地區幾乎沒有PS點。

圖4 PS-InSAR技術永久散射體(PS點)分布

利用重去平的干涉圖建立線性模型，估算變形率、殘余高度等。大氣效應的估算和清除通過兩種濾波程序進行，分別是大氣高通濾波和大氣低通濾波。對大氣校正后的結果進行地理編碼，即可得到蒲縣地區2019-01—2019-12年平均沉降速率(圖5)。對PS-InSAR技術所得結果圖進行觀察分析可知，研究區域內形變速率范圍為-53～48 mm/a。研究區域絕大部分地區處于穩定狀態，形變速率范圍為-15～10 mm/a。

沿雷達視線方向LOS (line of sight)地表抬升的地區用紅色表示，發生地面沉降的區域顯示為藍紫色

2.2 基于SBAS技術的蒲縣滑坡形變監測

為驗證PS-InSAR方法監測的準確度與可靠性，采用SBAS方法在同一數據源的基礎上對山西蒲縣地區進行滑坡監測，SBAS處理流程如圖6所示。

圖6 SBAS處理流程

采用蒲縣地區2019-01—2019-12的24景Sentinel-1A數據，設置臨界基線最小、最大百分比分別為0和2，時間基線閾值設定為60 d，系統自主選取2019-06-18的影像為超級主影像。根據設定得到閾值，共生成170個干涉像對，圖7為SBAS技術時空基線分布圖。

圖7 SBAS時空基線分布

為使實驗結果在距離向和方位向的分辨率盡可能保持相同的效果，實驗將方位向和距離向多視數分別設置為1和4，與超級主影像保持一致。選用常被應用于植被密集或環境潮濕地區的最小費用流法(minimum cost flow， MCF)進行相位解纏，將解纏相關系數閾值設置為0.23，利用Goldstein濾波方法以提高研究區中低相干位置的濾波強度。

打開存放差分干涉成果的文件夾，仔細查看每一個像對的相干性圖片以及解纏后的圖片，用連接圖編輯工具剔除相干性低和含有明顯大氣相位的像對。軌道精煉與重去平是為了估算差分后依舊殘存的恒定相位和解纏后還留有的相位坡道，利用GCP點對所有的數據對執行重去平以消除地形相位。在實驗過程中選取了遠離形變區域、孤立相位、殘余地形條紋等位置，位于相對穩定地點的70余個GCP點。在第一次反演中，對相干點建模并構建方程組矩陣，由于線性模型的穩定性最佳，因此，實驗選用線性模型求解所有影像對的變形率和新的DEM等。基于第一次反演所得形變速率，第二次反演是通過設置大氣高通和大氣低通來估計和清除實驗中的大氣影響，獲得更為準確純凈的時間序列中的位移結果并進行地理編碼，將其轉換到地理坐標系，獲得形變信息和年平均位移速率，如圖8所示。研究區域形變速率范圍為-70～52 mm/a，有4處地表沉降較為明顯的疑似滑坡隱患地點，大部分地表形變在-10～10 mm/a，趨于穩定狀態。通過SBAS方法采集的蒲縣地區2019年清晰明了的時序地表變化過程如圖9所示。

紅色表示雷達視線方向存在抬升，藍色表示下沉

圖9 基于SBAS技術的蒲縣地區2019年地表累計形變

3 時序形變結果分析

3.1 結果對比

基于PS-InSAR技術和SBAS方法，利用哨兵升軌SAR數據分別獲取蒲縣地區年平均形變速率，通過目視解譯選取出兩種技術所得結果中共同存在的4處較為明顯和典型的形變位置進行對比分析，所選位置分布如圖10所示，包括太林鄉半溝村羅克線沿線、水泉窊015鄉道兩側、郭家山周邊、井子洼村趙克路附近，具體地理位置如表3所示。圖11為4處形變速率放大圖，從放大圖中可以清楚地看到4個疑似隱患點的沉降信息。除圖11所示4處較明顯的疑似滑坡隱患點，研究區域仍存在部分升降不一的形變點，為使結果更加精準可靠，只對上述4處共同存在的疑似滑坡點做進一步對比分析。

圖11 4處疑似滑坡點SBAS技術形變速率放大圖

表3 4處疑似滑坡點詳細地理位置

圖10 4處疑似滑坡點位置

為了更加精準地分析PS-InSAR技術與SBAS技術對于同一數據源的結果精度對比情況，基于4處疑似滑坡點累計形變量的詳細信息，計算每一期數據的差值，并用兩種數據所得結果的平均值代替真值計算其均方根誤差(RMSE)，RMSE代表每個數據與真值間的偏離程度，均方根誤差越大，表示兩種技術獲取的結果差異越大，可靠性越低；反之，則表示測量的精度越高，進一步說明時序InSAR技術識別并監測滑坡的有效性，如表4所示。

由表4可得，一號疑似滑坡點PS-InSAR技術與SBAS技術所得最終累計沉降量分別為-51.23 mm 和-52.32 mm，差值為1.09 mm，RMSE范圍在0～2.25 mm，RMSE均值為0.75 mm；二號疑似滑坡點PS-InSAR技術與SBAS技術所得最終累計沉降量分別為-48.03 mm和-49.51 mm，差值為1.48 mm，RMSE范圍在0～2.06 mm，RMSE均值為0.77 mm；三號疑似滑坡點PS-InSAR技術與SBAS技術所得最終累計沉降量分別為-40.87 mm和-42.42 mm，差值為1.55 mm，RMSE范圍在0～3.25 mm，RMSE均值為1.22 mm；四號疑似滑坡點PS-InSAR技術與SBAS技術所得最終累計沉降量分別為-46.91 mm和-39.80 mm，差值為7.11 mm，RMSE范圍在0～3.55 mm，RMSE均值為1.05 mm。由上述數據可得PS-InSAR技術與SBAS技術監測結果的均方根誤差在0～4 mm，均值在0.5～1.5 mm，充分說明了PS-InSAR技術與SBAS技術的一致性較高。

表4 4處疑似滑坡點累計形變量的詳細信息

為了使兩種時序InSAR技術所得監測結果對比的更為直觀簡明，基于上述4處疑似滑坡點的累計形變量數據，對選取出的4處點繪制累計形變折線圖進行對比分析。圖12為識別出的4處疑似滑坡點2019-01-01—2019-12-15的LOS向累積形變位移量對比。總體看來，兩種技術所得結果的沉降趨勢較為一致，說明2019年蒲縣地區部分地區發生地表沉降，有疑似滑坡的可能，更進一步說明了時序InSAR識別并監測滑坡的可行性與準確性。1號疑似滑坡點在1—5月地表相對穩定，6月出現小幅沉降至20 mm，7月地面出現略微抬升，8—10月緩慢沉降，并在11、12月急劇下沉至-50 mm左右。2號疑似滑坡點1—6月緩慢沉降至-21 mm左右，并在7、8月出現小幅抬升，隨后逐漸沉降至-49 mm左右。3號疑似滑坡點沉降成線性變化，4號疑似滑坡點1—4月地表處于穩定狀態，隨后呈線性下沉模式。

圖12 4處疑似滑坡點累計形變位移量PS-InSAR與SBAS監測結果對比

3.2 可能性誘因分析

蒲縣是典型的溫帶大陸性氣候，雨熱同期。由上述分析可知，疑似滑坡點7、8月地面出現略微抬升，可能是由于雨季降雨量驟增的緣故，地表水滲入地下，導致地下含水層增厚。11、12月地表下沉明顯，可能是由于冬季降水量過少，當地居民大量使用地下水的緣故導致地下水位下降。3號疑似滑坡隱患點位于郭家山周邊，郭家山地區存在采煤區，該點的沉降疑似受到煤炭開采等人類活動的影響。

4 結論

利用時序InSAR技術探測研究蒲縣地區的滑坡形變，獲取研究區域內2019年年平均形變速率和累計形變位移量，確定太林鄉半溝村羅克線沿線、水泉窊015鄉道兩側、郭家山周邊、井子洼村趙克路附近等4處位置為疑似滑坡點，并結合煤炭開采、人類活動、降雨量信息等，對該地區分布的疑似滑坡隱患進行早期識別與監測預警，為該地區滑坡誘因提供科學的依據，主動規避自然地質災害風險，對未來滑坡災害的監測與預警提供思路與方法。