聯合光學和SAR遙感影像的山區公路滑坡易發性評價方法

余紹淮, 徐喬, 余飛

(中交第二公路勘察設計研究院有限公司,武漢 430056)

0 引言

隨著我國公路建設事業的發展,建設重點正逐步由東部向中西部地區、由平原微丘區向困難復雜的重丘山嶺區轉移,快速準確地對山區公路路線走廊區域進行滑坡易發性評價,可為后續道路的選線提供重要的參考資料; 但艱險山區地形、地質等基礎資料嚴重匱乏,且氣候條件復雜,地面數據獲取困難,野外作業風險高。遙感技術的快速發展,為復雜山區的滑坡易發性評價提供了新的技術手段,當前遙感技術已廣泛應用于區域性滑坡易發性評價中[1-4]。

滑坡災害易發性評價是對滑坡災害的時空分布和發生概率進行預測,可為滑坡災害風險管理提供重要的決策依據。當前應用較多的易發性評價方法有2大類: 第一類是基于知識驅動的評價方法,如層次分析法[5]、專家打分法[6]等經驗模型; 第二類是基于數據驅動的評價方法,包括邏輯回歸[7]、信息量法[8]、確定系數法[9]等統計分析模型,及人工神經網絡[10]、支持向量機[11]、隨機森林[12]等機器學習模型。與經驗模型和統計分析模型相比,基于數據驅動的機器學習方法能更好地分析評價因子和滑坡之間的非線性關系。隨著對地觀測技術的飛速發展,地表觀測數據越來越豐富,基于數據驅動的機器學習模型已廣泛應用到滑坡災害易發性評價模型中[13-14]。

基于長時間序列的SAR影像獲取的地表形變速率可直接反映滑坡運動狀態,有助于有效識別不穩定區域,對滑坡監測具有重要意義[15-17]; 但僅利用地表形變速率進行滑坡易發性評價,而忽略誘發滑坡災害的地形、地貌、地質等內在因素,會使評價結果易受地表形變誤差影響。目前,基于機器學習的滑坡易發性評價模型中較少使用SAR影像,更多使用高程、坡度、坡向等地面靜態數據,忽略了地表形變速率等地表動態數據,致使評價方法精度不高,實用性不強,難以推廣到其他區域[18-19]。因此,在評價方法中增加利用地表形變速率,可使評價結果與區域滑坡實際分布情況更相符,進一步提高結果的準確性與可靠性。

為此,針對上述不足,本文在綜合利用光學和SAR遙感影像基礎上,提出一種結合地表形變速率的滑坡易發性評價方法,該方法先利用光學遙感影像和地形數據提取多種滑坡災害靜態因子; 然后利用隨機森林(rondom forest, RF)模型計算滑坡易發性風險初始等級; 最后利用長時間序列SAR影像計算的地表形變速率對易發性風險初始等級進行修正,從而得到更加準確可靠的滑坡易發性評價分區圖,為后續山區公路地質選線提供數據支撐。

1 基于多源遙感數據的滑坡易發性評價模型構建

1.1 滑坡易發性評價因子提取

滑坡災害的發育是一個復雜的非線性過程,受多種因素綜合影響,包括地形地貌、地質構造、地層巖性、巖體結構及地下水等內在地質因素,也有降雨、人工開挖、地震等外在誘發因素; 因此,充分提取滑坡孕災環境信息對滑坡災害易發性評價至關重要。為快速準確建立路線區域的滑坡易發性評價模型,利用高分遙感影像及DEM數據提取多種滑坡災害因子; 其中,基于高分遙感影像的滑坡災害因子提取主要是通過對區域地質、地形地貌及生態景觀環境信息進行遙感解譯識別,提取因子包括地層巖性、地質構造、地形地貌等; 基于DEM數據的因子提取主要是通過GIS分析獲取,主要有坡度、坡向、坡長、水流流向、水流累積量等。

本文所用的滑坡易發性評價因子如表1所示。由于DEM數據和遙感影像的分辨率不同,致使基于二者所提取的滑坡災害因子尺度不同,考慮到可免費獲取的DEM數據的分辨率為30 m×30 m,本方法采用30 m×30 m的柵格單元作為滑坡易發性評價單元,并在此基礎上提取滑坡災害因子。

表1 滑坡易發性評價因子表

1.2 基于隨機森林的滑坡易發性風險等級劃分

隨機森林(random forests,RF)是由Breiman L于2001年提出的一種以決策樹為基礎分類器的集成機器學習算法。該算法具有分類表現優異、人工干預少、運算非?？斓葍烖c,在遙感圖像分類、特征選擇等方面應用廣泛[20-21]。

隨機森林的基本原理是先利用Bagging重抽樣方法從原始樣本中抽取多個樣本集; 然后對樣本集構建多個分類回歸樹(classification and regression tree,CART),其中每個CART樹的內部節點分裂時,對特征集進行一次隨機抽樣,并在抽取的特征集中進行最優選擇; 最后,利用所有CART樹的計算結果,通過均值或投票得出最終計算結果。隨機森林模型的構建過程見圖1。

圖1 隨機森林模型構建示意圖

由圖1可知,隨機森林模型中引入了樣本隨機抽樣及特征隨機抽樣,有效降低模型對樣本噪聲和異常值的敏感度,提高評價模型的準確率和穩定性。為此,本文基于滑坡災害因子集,利用滑坡樣本與非滑坡樣本,訓練基于隨機森林的滑坡易發性評價模型,完成滑坡易發性風險初始等級的劃分,主要包括以下步驟:

1)從N個原始滑坡訓練樣本中以有放回的方式重復取樣N次,得到一個訓練樣本集; 重復上述過程k次,得到k個訓練樣本集;

2)針對每一個訓練樣本集,通過隨機選取滑坡評價因子集作為分裂特征集,并選擇最優的分裂方式進行內部節點的迭代分裂,且分裂過程中不做減枝處理,迭代完成后得到k棵CART樹;

3)將生成的k棵CART樹組合成隨機森林,對輸入的滑坡評價因子集進行分類預測,統計每棵CART樹的分類結果,計算出滑坡易發性概率;

4)依據得到的滑坡易發性概率,利用自然斷點法進行分級處理; 將研究區域劃分低易發區、中易發區、高易發區和極高易發區4個等級,完成滑坡易發性風險初始等級的劃分。

1.3 結合形變特征的滑坡易發性風險等級修正

上述滑坡易發性風險初始等級的劃分是基于地面相關靜態參數如高程、坡度、坡向等計算得到,而地表形變速率反映了地表動態變形,可直接反映滑坡的運動狀態; 因此利用地表形變速率對滑坡易發性風險初始等級進行修正,可提高結果的準確性。為此,本文利用長時間序列SAR影像得到地表形變速率對易發性風險初始等級進行修正,從而得到更加準確有效的滑坡易發性評價分區圖。

PS-InSAR[22]與SABS-InSAR[23]是2種應用最為廣泛的時序InSAR技術,考慮到SABS-InSAR技術采用更加自由的時、空基線閾值來構成干涉對組合,削弱了空間失相關和大氣延遲的影響,可得到長時間緩慢地表的形變規律,更適用于自然地表形變速率的監測。因此,本文采用SABS-InSAR技術獲取該區域的形變速率,具體流程如圖2所示。

圖2 SABS-InSAR處理流程圖

考慮到SABS-InSAR技術只能探測沿雷達視線方向的形變信息,而滑坡多沿斜坡面進行滑動,本文利用SAR成像信息及DEM數據,將上一步獲取的形變速率轉化為沿坡度方向的形變速率,以反映斜坡面的地表形變信息。斜坡面的形變速率(Vslope)的計算表達式為:

Vslope=VLOS/C

,

(1)

式中:VLOS為雷達視線方向的形變速率;C為坡面位移和雷達視線位移的比例系數,公式為:

,

(2)

hLOS=cosθ

,

(3)

,

(4)

,

(5)

式中:φ為地形坡度;α為坡向;θ為SAR衛星成像的入射角;ω為衛星軌道方向和正北方向的夾角;δ為地理正西方向與SAR衛星成像視線向的夾角。

SBAS-InSAR技術獲取的坡向形變速率反映了區域坡面的動態變形狀態,將坡面形變速率加入到滑坡易發性評價中,可提高模型的敏感度與準確性。為此,按照評價單元大小,將所提取的坡面形變速率圖重采樣至30 m×30 m,并依據評價單元所對應的形變速率大小,將評價單元劃分為低形變區([0,15) mm/a)、中等形變區([15,30) mm/a)、高形變區([30,45) mm/a)及極高形變區(≥45 mm/a)4類; 然后,將滑坡易發性風險等級與形變等級關聯,建立風險等級更新矩陣,實現對滑坡易發性風險初始等級的修正,得到最終的滑坡易發性評價分區圖?；乱装l性風險等級更新矩陣如表2所示。

表2 滑坡易發性風險等級更新表

表2中的數值表示滑坡易發性風險初始等級修正的級別,其中0表示不修正,+1,+2,+3分別表示將風險等級提高1級、2級、3級。例如,某個評價單元易發性風險初始等級為低易發區,形變等級為高形變區,根據表2,該評價單元的風險等級需提高2級,則該評價單元的最終風險等級為高易發區。

2 工程實踐與分析

2.1 研究區概況

青海省沿黃公路共和至大河家段公路地處青藏高原東北緣,青海省東南部,起點位于青海省共和縣,終點位于甘肅省積石山縣,是青海省第一條沿黃河公路; 建設高質量的沿黃河公路,對于開發黃河上游水能資源,促進黃河谷地開發利用,完善公路網建設,適應地區經濟快速發展的交通需求,具有十分重要的現實意義和社會意義。

工程地處青海省拉脊山斷裂南側,受拉脊山斷裂帶的控制作用,并隨著青藏高原的急劇抬升和黃河不斷下切,導致路線經過的龍羊峽至拉西瓦、李家峽庫區、隆務峽至公伯峽等黃河峽谷段兩岸地勢陡峭,深切河谷發育,地形地貌復雜多樣,泥石流、滑坡、崩塌等不良地質現象發育。為此,本文選取選擇隆務峽至公伯峽段進行工程實踐,以驗證本文提出的聯合光學和SAR遙感影像的山區公路滑坡易發性評價方法的準確性和有效性。研究區的范圍與地形如圖3所示。

圖3 研究區范圍與地形示意圖

2.2 實驗數據及處理

為準確評價研究區的滑坡易發性風險等級,選擇QuickBird衛星影像、Sentinel-1A衛星影像及ASTER GDEM數據作為評價模型的基礎數據。其中,QuickBird衛星影像用于提取滑坡易發性評價因子與歷史滑坡區的解譯,其分辨率為全色0.61 m,多光譜2.44 m?？紤]到研究區內的地形起伏較大,先對QuickBird影像進行幾何糾正與空間配準,以消除影像上存在的幾何畸變; 然后對QuickBird衛星的多光譜與全色波段進行融合處理,以提高QuickBird影像的空間分辨率的同時保留其多光譜特性。在此基礎上,利用處理后QuickBird影像提取出植被指數、地表覆蓋等多種評價因子; 然后對QuickBird影像進行滑坡遙感解譯,以得到研究區的歷史滑坡區,并在此基礎上選擇訓練樣本。圖4為研究區的歷史滑坡遙感解譯圖。

圖4 研究區歷史滑坡遙感解譯圖

Sentinel-1A衛星影像用于提取研究區內的坡面形變速率,作為后續滑坡易發性風險初始等級調整的依據。實驗選擇覆蓋研究區范圍的37景2018年1月—2020年6月的Sentinel-1A升軌影像,影像的成像方式為干涉寬幅模式,選用VV極化方式的影像進行干涉處理,并利用精密軌道文件進行軌道誤差校正,所用Sentinel-1A的影像信息如表3所示。

表3 Sentinel-1A 影像信息表

本文利用ENVI 5.3中的SARScape模塊,對上述Sentinel-1A影像進行裁剪、配準、干涉、去平、濾波、解纏及地理編碼等列干涉處理; 其中,所選擇的超級主影像日期為2018-05-14,所設置的時間基線閾值與最大空間基線閾值分別為120 d與5%,所利用的干涉相對數目為165,平均絕對時間基線與絕對空間基線分別為70 d與53.88 m; 在此基礎上,利用SBAS-InSAR方法得到研究區沿雷達視線方向的形變速率,結果如圖5所示。

圖5 研究區2018—2020年沿雷達視線方向形變速率分布圖

利用上述結果,結合Sentinel-1A衛星成像信息及ASTER GDEM數據,利用前文提及的轉換方法,將沿雷達視線方向的形變速率轉換為沿坡面方向的速率。由于斜坡面和雷達視線方向的夾角接近90°時,會導致斜坡面形變速率值趨近無窮大,結合區域的實際情況與其他學者的研究經驗[24-25],本文將比例系數C的閾值固定為±0.3,即:

。

(6)

另外,考慮到斜坡體的形變應沿著斜坡面向下運動,只有形變速率為負值的形變點才能反應斜坡體的實際形變狀態,故剔除形變速率為正值的形變點,從而得到研究區沿斜坡面的形變速率分布圖,結果如圖6所示。

圖6 研究區2018—2020年沿斜坡面方向形變速率分布圖

2.3 結果分析與討論

利用基于高分辨率QuickBird影像解譯的歷史滑坡區域,從中選取257個樣本點作為訓練樣本,并從影像其他區域隨機選擇非滑坡點,二者比例保持1∶1,提取出各樣本點的評價因子數據作為隨機森林模型的訓練樣本。本文利用隨機森林算法得到的研究區滑坡易發性初始評價分區圖如圖7所示。

圖7 研究區滑坡易發性初始評價分區圖

對比圖6與圖7可知,部分形變速率較大的區域被劃分為低易發區,而長時間存在形變的區域往往發生滑坡的風險較高,因此僅考慮高程、坡度、坡向等地面靜態因子對研究區的滑坡易發性進行評價,會導致評價結果的可靠性不高,與區域實際情況存在較大差異。為此,本文利用圖6所示的斜坡面方向的形變速率對圖7的滑坡易發性風險初始等級進行修正,得到最終的滑坡易發性評價分區圖,結果如圖8所示。經統計,圖8所示的滑坡易發性評價分區圖中,極高易發區、高易發區、中易發區及低易發區的面積占比分別為7.73%,13.67%,20.34%和58.26%,其中高易發區及極高易發區占比超過20%,表明研究區內大片區域受到滑坡的潛在威脅。通過前文解譯的歷史滑坡區可知,大部分滑坡區位于黃河峽谷及黃河兩岸,與圖7所示的高易發區、極高易發區的分布基本一致。此外,選取部分重點區域,將形變速率圖、初始評價分區圖與最終的評價分區圖進行對比分析,以驗證本文方法的有效性,對比圖見圖9。

圖8 研究區滑坡易發性評價分區圖

(a) 形變速率圖 (b) 滑坡易發性初始評價分區圖 (c) 結合形變特征的易發性評價分區圖

由圖9(a)可知,該區域存在大量形變速率較高的區域,其中A和B區域形變速率最高,具有較大的滑坡潛在風險,應具有極高的易發性風險等級; C,D,E區域也存在明顯形變,應具有較高的易發性風險等級。但圖9(b)所示的初始評價分區圖將上述區域主要劃分為低易發區,少量區域劃分為中易發區與高易發區,與研究區形變速率較高的現狀不符,而圖9(c)所示的結合形變特征的評價分區圖將上述區域主要劃分為極高易發區與高易發區,與研究區的實際狀況一致。對比分析結果表明,僅利用地表靜態因子得到的滑坡易發性評價分區圖精度有限,增加利用地表形變特征可顯著提高評價結果的準確性與實用性。

頻率比是指滑坡易發性等級內滑坡樣本分級占比與易發性分級面積占比的比值,可用于定量評價滑坡易發性分級的準確性[26-27]。為進一步驗證本文方法的準確性,利用研究區的歷史滑坡遙感解譯圖,采用頻率比檢驗的方式,對所得到的滑坡易發性評價分區圖進行定量評價,評價結果如表4所示。

表4 滑坡易發性評價精度統計表

由表4可知,上述2種方式所得到的滑坡易發性評價分區圖中,高易發區與極高易發區的分級面積占比均不超過研究區總面積的22%,但卻分布有 37%以上的滑坡樣本,遠超過其他易發區等級,且頻率比值從低易發區到極高易發區均顯著增大,極高易發區頻率比值也遠大于其他易發性等級,表明上述2種方式能有效地評價研究區的滑坡易發性。

利用形變特征對滑坡易發性初始評級結果進行修正后,有26.73%的滑坡樣本分布于極高易發區,22.28%的滑坡樣本分布于高易發區,均高于初始分級圖中的21.88%與15.18%,改正后分區圖中的極高易發區的面積占比從6.99%提高到7.73%,頻率比值從3.130 2提高到3.458 0,精度提升明顯,表明增加利用地表形變特征可顯著提高評價結果的準確性,尤其是具有較高風險等級的區域。修正后的分區圖中低易發區的滑坡樣本比例也顯著降低,頻率比值從0.749 1降低到0.555 6,但仍有32.37%的滑坡樣本區域位于低風險區,這主要是由于使用的DEM數據精度與分辨率不高及滑坡樣本點包含部分處于穩定狀態的古滑坡所造成?？偟膩碚f,本方法得到的滑坡易發性評價分區圖與研究區實際情況基本一致,驗證了本文聯合光學和SAR遙感影像對山區公路走廊進行滑坡易發性評價的有效性。

3 結論與展望

本文采用聯合光學和SAR遙感影像進行山區公路滑坡易發性評價的方法,準確提取了多種靜態與動態滑坡災害因子; 基于隨機森林算法綜合利用地質、水文、地形地貌及地表覆蓋等多種因子,實現了滑坡易發性風險初始等級的快速評估; 利用SABS-InSAR技術獲取的地表形變因子,對滑坡易發性風險初始等級進行精準修正,并有效開展艱險山區公路走廊范圍內的滑坡易發性評價工作。通過在青海省沿黃公路隆務峽至公伯峽段中的應用實踐,表明本文方法可快速準確地得到山區公路走廊范圍內的滑坡危險性評價分區圖,適用于山區公路遙感地質勘察,可為后續公路地質選線提供數據參考。

但是,本文提出的滑坡易發性評價方法仍存在以下不足: ①滑坡樣本需人工選擇,影響工作效率與評價精度; ②采用的DEM數據精度與分辨率有限,降低了滑坡易發性評價分區圖的準確性。因此,下一步工作應改進滑坡易發性評價算法,減少人工干預,采用高精度的DEM數據,提高滑坡易發性評價工作的效率與準確性。