李利峰, 張曉虎, 鄧慧琳, 韓六平
(1.貴州工程應用技術學院土木建筑工程學院,畢節 551700;2.貴州工程應用技術學院理學院,畢節 551700)
滑坡常發生在山區與丘陵地帶,是最危險的地質災害之一[1]。滑坡每年造成大量的人員傷亡與財產損失。據統計,世界上每年因滑坡造成的死傷人數約為1 000人,經濟損失高達40億美元,并且呈現逐年遞增的趨勢[2]。中國是一個滑坡災害頻發的國家,特別是在西北地區,由于環境惡劣以及人口集中,滑坡對當地居民生命與財產形成巨大的威脅。對每個滑坡的治理需要投入大量的人力與財力,在當前經濟狀況下并不具備可執行性。而對區域滑坡發生進行預測可以起到防災減災、指導工程部署以及合理規劃土地等作用,具有重要的理論與實踐意義。
目前,常見的滑坡易發性分析方法主要分為兩大類:傳統的統計學分析方法與機器學習法。傳統的統計學分析法所建立的模型自變量的內部參數是由滑坡原始數據確定,可以避免人為因素的干擾,得到的參數結果更加客觀,但在選取評價因子的過程中,更容易受主觀因素的干擾。常見的傳統的統計學分析法有層次分析法[3-4]頻率比模型[5-6]、確定性系數模型[7]、證據權模型[8]、信息量模型[9]以及熵指數模型[10]等。機器學習模型可以避免各評價因子間互相依賴的問題,但由于模型的計算量大并且對模型參數的選取也較為困難,得到的評價結果很容易出現過擬合問題。最常見的機器學習算法有人工神經網絡法[11]、支持向量機[12]、隨機森林[13]、邏輯回歸算法[14]等。
針對以上兩類模型各自的優缺點,近年來,已有許多學者嘗試將這兩類模型進行耦合,形成混合模型應用于滑坡易發性評價中,并且取得了很好的預測結果。例如,謝洪斌等[15]利用模糊證據權法對汶川震區岷江流域雁門鄉至映秀段滑坡災害進行評價,得到較好的評價結果;安凱強等[16]將信息量模型與支持向量機模型相結合對三峽庫區的滑坡災害進行預測,將91.11%的滑坡預測為易發性高的區域。因此,將熵指數(index of entropy, IOE)模型與邏輯回歸(logistic regression, LR)模型相耦合,形成IOE-LR模型并用于藍田縣滑坡災害易發性評價中。該模型既可以避免評價因子在選取過程中受人為因素的干擾,又可以減少機器學習過程中出現過擬合的問題。但其不足之處是模型參數的選取是否恰當會直接影響滑坡預測結果可靠性。
藍田縣屬于西安市轄縣之一。地理位置為東經109°07′~109°49′,北緯33°50′~34°19′。東西長64 km,南北寬55 km,總面積2 007 km2。地形南北高、中間低,整體由東南向西北傾斜。區內最高海拔2 420 m,最低海拔397 m。按地貌單元將全縣劃分為河谷平原、山前洪積扇、黃土臺塬、黃土丘陵和基巖山地5種地貌。藍田縣屬暖溫帶半濕潤大陸性季風氣候,多年平均氣溫為13.1 ℃。境內河流眾多,屬黃河流域渭河水系,其中主要以灞河、浐河、零河及其支流水系為主。境內褶皺、斷裂較為發育,新構造運動十分活躍,垂直差異性升降運動明顯,進而增加了滑坡災害發生的幾率。人類工程活動主要有選址修房、開墾耕作、礦業開采、修建公路。全境在冊滑坡點共有227處,其中土質滑坡占80.6%,巖質滑坡僅占19.4%。區內較為典型滑坡點以及研究區地理位置以及滑坡點分布圖如圖1、圖2所示。

圖1 姜灣村滑坡Fig.1 The landslide of Jiangwan Village
熵指數模型是一種二元統計模型,它既可以單獨對數據進行分類預測,也可以作為構建混合模型的輸入數據。熵是指一個系統的不穩定或不確定性的程度。在滑坡易發性分析評價中,熵代表不同評價因子對滑坡災害發生的影響程度。熵指數值越大,代表評價因子的影響作用越大,反之亦然。各評價因子的權重值由以下公式計算得到:

圖2 地理位置以及滑坡災害點分布Fig.2 Geographical location and landslide points distribution of the study area

(1)
(2)
(3)
Himax=log2S
(4)
(5)
(6)
Wi=Ii×Pi
(7)
式中:Pij為各評價因子分級的頻率比;a與b分別表示對應分級的面積百分比與滑坡百分比;(Pij)代表概率密度;Hi、Himax表示熵值;S為評價因子分級數;Ii為評價因子信息率;Wi為評價因子權重。
邏輯回歸是一種用于二分類或多分類的分類預測模型。在滑坡易發性評價分析中,各評價因子作為模型的輸入自變量,滑坡發生與否作為分類因變量(0代表滑坡不發生,1代表滑坡發生)。設P為滑坡發生的概率,1-P為滑坡不發生的概率,P的取值范圍為[0,1]。P的表達式為
(8)
式(8)中:x1,x2,…,xn代表各評價因子;β1,β2,…,βn為邏輯回歸系數;α為常數項,表示在不受任何有利或不利于滑坡發生因素影響的條件下,滑坡發生與不發生概率之比的對數值。
對P/(1-P)取自然對數,得到邏輯回歸方程為
(9)
式(9)中:Z為滑坡發生概率的目標函數。
上述IOE模型得到的Pij具有相同的維度,可以避免出現滑坡災害點與評價因子間出現線性相關的問題以及在建模過程減少噪聲的出現。將各評價因子按照Pij進行重分類,將Pij作為邏輯回歸模型中的輸入數據,構建LR-IOE模型,計算得到混合模型下的α與β。
用到的數據主要有滑坡點屬性數據、各評價因子屬性數據。其中滑坡點屬性數據主要是通過收集的歷史滑坡數據、野外調查以及衛星影像遙感解譯等手段獲取;各評價因子屬性數據是通過ArcGIS軟件對各評價因子圖層提取獲得。
在對藍田縣境內的地質環境條件、滑坡發育特征、形成條件等分析的基礎上,初步選取高程、坡度、坡向、曲率、降雨量、距水系距離、地層巖性、距斷層距離、距道路距離、歸一化植被指數(normalized difference vegetation index,NDVI)共10類影響因素作為評價因子。利用ArcGIS軟件提取各因子圖層,如圖3所示。其中,高程、坡度、坡向、曲率圖層主要提取自30 m×30 m分辨率的數字高程模型(digital elevation model,DEM)影像;降雨量圖層通過藍田縣氣象站點數據經克里金插值獲取;地層巖性、斷層、水系、道路圖層分別通過對1∶50 000地質圖、水系圖、路網圖進行矢量化,并轉為柵格得到。

圖3 評價因子圖層Fig.3 Evaluation factors layer
對研究區進行滑坡易發性評價時,評價因子選取的數量與質量都會影響評價結果。因此,并不是所有的評價因子對評價結果都能產生積極作用。當評價因子間存在多重共線問題時,不僅會增加模型的復雜度,而且會出現模型過擬合問題,大大降低模型的泛化能力,降低預測結果的準確度。因此,分別采用皮爾森相關系數(Pearson correlation coefficient,PCC)、方差膨脹因子(variance inflation factor,VIF)以及容忍度(tolerance,TOL)3種指標分析因子變量之間的多重共線性。
PCC常用來反映變量間的線性相關性,其絕對值越大,說明變量間的多重共線性越強。假設有樣本數據集(Xi,Yj)=(x1,y1), (x2,y2),…,(xn,yn),變量間的線性相關系數PCC的計算表達式為
(10)

一般認為,當0≤∣PCC∣≤0.3時,變量之間呈低度相關;當∣PCC∣>0.3時,變量之間相關度較高。
VIF與TOL也是判斷變量間是否多重共線的重要指標。VIF是將存在多重共線性時回歸系數估計量的方差與不存在多重共線性時回歸系數估計量的方差對比而得出的比值系數。TOL為VIF的倒數。通常認為,當0
對各評價因子的多重共線性進行分析,得到的結果如表1、表2所示。由表1可知,各評價因子之間相關系數的絕對值均小于或等于0.3,表明各因子間相對獨立;由表2可知,各評價因子的VIF值均小于2,TOL值均大于0.4,表明各因子之間多重共線性較弱。
對各因子分級是進行統計學分析的前提條件。根據滑坡災害點在區域分布的空間特征以及查閱相關區域文獻,對各評價因子進行分級,得到各評價因子的分級狀態,如表3所示。其中,坡度、坡向、距斷層距離、距水系距離以及距道路距離按照相等間隔法劃分;高程、曲率、NDVI、降雨量按照自然間斷點法劃分;地層巖性分類編號如下:1(太古界太華群)、2(下元古界陶灣組)、3(中元古界秦嶺群)、4(古近系與新近系砂泥巖)、5(第四系風成黃土)、6(第四系全新統沖積層)、7(第四系全新統洪積層)、8(加里東期花崗巖)、9(燕山期花崗巖)、10(喜山期花崗巖)。采用頻率比法對各評價因子與滑坡點的空間分布關系進行分析,得到各因子的頻率比如表3所示。
從表3中可以得出,高程在782~1 084 m,坡度在8°~24°,坡向為朝東或西南方向,曲率為0.44~0.12的凸型坡,其頻率比值相對較大,說明在此地形條件下滑坡災害最容易發生。斷層、水系與道路隨著距離的增加,中間稍有起伏波動,但總體呈現出頻率比值逐漸減小趨勢,說明距離水系、斷層、道路越近,越容易誘發滑坡。NDVI在0.02~0.07范圍內,頻率比最大,說明滑坡多發生于植被較少的區域。
將整個研究區按照30 m×30 m柵格大小進行劃分,共劃分為2 230 784個評價單元。利用ArcGIS軟件里的柵格轉點工具提取研究區各圖層的屬性值,建立研究區屬性數據庫。利用多值提取至點工具提取滑坡點各因子屬性值,建立滑坡屬性數據庫。
根據式(1)~式(7),利用Excel軟件求得的各評價因子權重值如表3所示。最終的滑坡易發性區劃圖由以下公式求得:
(11)
式(11)中:LSIIOE表示滑坡易發性指數;i表示評價因子數;e表示分級數最多的評價因子的分級數量;fi表示各評價因子的分級數;C為各評價因子分級后重分類值。

表1 評價因子之間的皮爾森相關系數Table 1 PCC between factors

表2 評價因子的方差膨脹因子與容忍度Table 2 VIF and TOL for factors

表3 評價因子與滑坡點空間分布關系Table 3 Spatial relationship between each factor and landslide
利用式(11),求得區域的易發性指數LSI范圍為2.41~9.48。LSI越大,表示滑坡越容易發生,LSI越小,表示滑坡發生的概率越小。采用自然間斷點法將LSI分為5個等級,分別為極低易發等級(2.41~4.63)、低易發等級(4.63~5.65)、中易發等級(5.61~6.71)、高易發等級(6.71~7.68)、極高易發等級(7.68~9.48),生成最終的滑坡易發性區劃圖,如圖4所示。其中,極高-高易發區占研究區面積的45.3%,滑坡占總滑坡數的86.5%,滑坡點密度為0.216個/km。

圖4 基于熵指數模型的研究區滑坡災害易發性區劃圖Fig.4 Landslide susceptibility mapping based on IOE model in study area
將滑坡點與隨機選取的相同數量的非滑坡點組成總的訓練樣本集,從總訓練樣本集中隨機選取70%的滑坡點與相等數量的非滑坡點,共計318個作為訓練樣本集,剩余的作為測試樣本集。采用R語言里的邏輯回歸對訓練樣本集進行學習,再將學習得到的模型用于測試樣本集,得到模型的預測正確率為82.5%;隨后將研究區屬性數據庫代入訓練好的模型中,得到研究區的滑坡災害易發性指數LSI,將LSI按自然間斷點法分為5個等級,分別為極低易發等級(0~0.16)、低易發等級(0.16~0.40)、中易發等級(0.40~0.52)、高易發等級(0.52~0.64)、極高易發等級(0.64~0.93),生成最終的滑坡易發性區劃圖,如圖5所示。其中,極高-高易發區占研究區面積的38.4%,滑坡占總滑坡數的84.3%,滑坡點密度為0.248個/km。

圖5 基于邏輯回歸模型的研究區滑坡災害易發性區劃圖Fig.5 Landslide susceptibility mapping based on LR model in study area
將式(2)求得的Pij作為輸入數據代入R語言的邏輯回歸模型進行學習,最終得到模型的預測正確率為85.6%。對整個研究區進行預測,得到研究區滑坡易發性指數LSI,同樣按照自然間斷點法將其劃分為5個等級,分別為極低易發等級(0~0.15)、低易發等級(0.15~0.31)、中易發等級(0.31~0.47)、高易發等級(0.47~0.65)、極高易發等級(0.65~0.97),生成最終的滑坡易發性區劃圖,如圖6所示。其中,極高-高易發區占研究區面積的30.4%,滑坡占總滑坡數的80.2%,滑坡點密度為0.298個/km。

圖6 基于熵指數與邏輯回歸耦合模型的研究區 滑坡災害易發性區劃圖Fig.6 Landslide susceptibility mapping based on IOE-LR model in study area
采用ROC曲線對3種不同模型的性能進行檢驗與比較,用ROC曲線下的面積(area under curve,AUC)分別表示訓練樣本與測試樣本的成功率與預測率,得到各模型成功率與預測率曲線如圖7、圖8所示。

圖7 成功率曲線Fig.7 The success rate curve

圖8 預測率曲線Fig.8 The prediction rate curve
IOE模型、LR模型、IOE-LR模型成功率曲線的AUC分別為0.735、0.742和0.805;預測率曲線的AUC分別為0.732、0.785和0.830,IOE-LR模型均有最高的成功率與預測率,表明IOE-LR模型整體要優于單一模型。
分別采用IOE模型、LR模型、IOE-LR模型對藍田縣滑坡災害易發性分析,取得以下結論。
(1)分別采用皮爾森相關系數、方差膨脹因子、容忍度3種指標對評價因子之間的多重共線性問題進行分析,結果表明,各選取因子之間多重共線性較低,可以認為各因子相對獨立。
(2)IOE、LR、IOE-LR3種模型在極高-高易發區內滑坡點密度分別為0.216、0.24、0.298個/km。說明IOE-LR模型在極高-高易發區內滑坡分布更加集中。
(3)采用ROC曲線對3種不同模型的性能進行檢驗與比較,IOE模型、LR模型、IOE-LR模型成功率曲線的AUC分別為0.735、0.742和0.805;預測率曲線的AUC分別為0.732、0.785和0.830,IOE-LR模型均有最高的成功率與預測率,說明IOE-模型整體上優于單一模型。