陳飛,蔡超,李小雙,錢乾
(江西理工大學,a.資源與環(huán)境工程學院;b.江西省礦業(yè)工程重點實驗室,江西 贛州341000)
滑坡易發(fā)性評價是滑坡危險性評價的基礎[1],所以滑坡易發(fā)性的科學劃分與精度都至關重要。隨著GIS 與計算機算法的發(fā)展,常見滑坡易發(fā)性評價模型為統(tǒng)計分析模型,其中簡單的統(tǒng)計分析模型有信息量模型[2-4]、層次分析法[5-6]、證據(jù)權方法[7-8]等,但簡單的統(tǒng)計分析模型在數(shù)據(jù)分析中存在著人為的主觀性,并存在大量而繁重的工作,導致易發(fā)性評價精度較低。隨著大數(shù)據(jù)的發(fā)展,為了解決人為的過度干預,并減輕地質災害評價中的大量工作,大部分學者開始展開對復雜統(tǒng)計模型的研究。常見的復雜算法有邏輯回歸模型[9-10],神經(jīng)網(wǎng)絡模型[11]與支持向量機模型[12]等。
本文以上猶縣為研究區(qū),首先根據(jù)遙感圖像與實際調查獲取滑坡易發(fā)性評價的影響因子,通過相關性分析,確定弱相關性的評價因子[13],并將各個評價因子進行分級,得到分級后的評價因子與上猶縣歷史滑坡單元建立上猶縣災害評價數(shù)據(jù)庫。其次在滑坡單元數(shù)據(jù)500 m 以外選取500 個非滑坡單元,通過Arcgis 空間連接得到上猶縣滑坡單元與非滑坡單元的空間數(shù)據(jù)??紤]到存儲空間與制圖效果,確定柵格單元為30 m 分辨率,將上猶縣劃分為1717906 個柵格單元,得到每個柵格單元的空間數(shù)據(jù),并對上述空間數(shù)據(jù)進行歸一化處理[14]。最后將滑坡與非滑坡單元的空間數(shù)據(jù)隨機分為80%作為訓練集、20%作為測試集。上猶縣柵格單元的空間數(shù)據(jù)作為預測集,輸入邏輯回歸模型、神經(jīng)網(wǎng)絡模型、支持向量機模型進行數(shù)據(jù)的訓練、測試與預測,并利用自然斷點法進行滑坡易發(fā)性等級的劃分。為了提高易發(fā)性分區(qū)的安全性,將上猶縣柵格單元取最大的易發(fā)性指數(shù),即引入MAX(LSI (ANN),LSI(LG),LSI(SVM))將上述3 種模型的易發(fā)性指數(shù)取最大值計算,提高高危險區(qū)與較高危險所包含的滑坡數(shù)。其具體流程如圖1 所示。
邏輯回歸(LR)是一個二分類事件,將0 設為非滑坡單元,1 為滑坡單元[15]。邏輯回歸(LR)揭示一個因變量與多個自變量的關系[16],能很好地解決二分類問題。并且LR 模型的自變量分別為影響滑坡單元與非滑坡單元的評價因子。其公式如式(1)所示。式中:Y 代表非滑坡與滑坡單元,取值為0 或1;X1j,X2j,…,Xnj為各個評價因子的第j 級別;a0,a1,…,an代表邏輯回歸系數(shù)。若發(fā)生滑坡的概率用P 來表示,其計算如下:
在神經(jīng)網(wǎng)絡中,常見的算法有前饋神經(jīng)網(wǎng)絡[17](FNN)、反向傳播算法[18](B-P)、多層感知器[19](MLP)。由于前饋神經(jīng)網(wǎng)絡是單項的多層結構,沒有周期與循環(huán),而MLP 利用反向傳播作為訓練,又以前饋為連接,有多個隱藏層,算法復雜,對于數(shù)據(jù)的運算時間過長。所以本文選取反向傳播算法(B-P)作為上猶縣易發(fā)性評價的研究方法。假設各個評價因子用輸入神經(jīng)網(wǎng)絡中的xi(i=1,2,…,n)來表示,wi表示各個神經(jīng)單元之間的權重,取特征函數(shù)為雙向正切函數(shù),其公式為:
式中:s 代表神經(jīng)元的輸入總和;y 代表神經(jīng)元的輸出;θ 代表神經(jīng)網(wǎng)絡的閾值。
B-P 神經(jīng)網(wǎng)絡是采用反向傳播進行參數(shù)優(yōu)化,調整權重。目的是使損失函數(shù)最小。假設神經(jīng)網(wǎng)絡的誤差函數(shù)為Ep,用來衡量B-P 神經(jīng)網(wǎng)絡的能力。則其公式為:
SVM 是將輸入的數(shù)據(jù)映射在一個高維的空間,然后在此高維空間上進行線性分類,從而實現(xiàn)在原空間坐標系的非線性分類[20]。假設訓練數(shù)據(jù)有n 個,其表示為(xi,yi),其中xi表示評價因子,yi為表示是否為滑坡。只考慮一個變量,則支持向量機的超平面為:
式中:w 為斜率:b 為截距。當有n 個數(shù)據(jù)輸入,其超平面為:
式中:WTX 為wixi的累加。在滿足零均值與等方差的情況下,回歸方程的參數(shù)用最小二乘法來確定,在保證輸入變量的實際值與估測值的差值平方和達到最小的情況下,從而確定回歸方程的系數(shù),其公式如式(8)所示:
式中:yi為在第i 個因子觀測的輸出的預測值;(yi-y?i)2(i=1,2,…,n)是實際值與預測值相差的平方。設ε 為不敏感函數(shù),每個觀測的誤差函數(shù)都計入損失函數(shù),進行疊加,當誤差函數(shù)小于ε,即誤差可忽略不計[21]。本文將支持向量機的核函數(shù)選用線性核函數(shù)。
上猶縣位于江西省贛州市西部,坐標為東經(jīng)114°~114°40′,北緯25°42′~26°01′,屬于大陸性氣候,平均氣溫18°。區(qū)內地勢西高東低,常年多雨,災害發(fā)育頻繁。圖2 為上猶縣位置與災害分布圖。
通過調查在上猶縣有滑坡點500 個,其范圍覆蓋上猶縣全境,最高的滑坡位于五峰山鄉(xiāng),其高度為200 m,為花崗巖巖性。由于大部分災害圍繞縣城區(qū)域,在暴雨時期,這些滑坡可能釀成地質災害,嚴重損害上猶縣的經(jīng)濟發(fā)展。研究區(qū)內以構造侵蝕低山為主,主要有志留系、第四系、泥盆系、奧陶系、石炭系、三疊系與寒武系等地層組成,其巖性主要由粉砂巖、礫巖、細花崗巖、變輝長巖、砂巖、閃長巖、細粒輝長巖等組成。其次為構造中山,由花崗巖與組成,丘陵位于寒武系地層,由巖漿巖組成。侵蝕性崗地主要位于泥盆系,有巖漿巖等巖性組成。侵蝕性河谷主要有砂巖,有河流沖積形成。
文中滑坡易發(fā)性評價數(shù)據(jù)來源: ①1∶2000 上猶縣地質災害一覽表; ②地理數(shù)據(jù)云下載landsat5 數(shù)據(jù);③地理數(shù)據(jù)云下載的30 m 分辨率得DEM 數(shù)據(jù);④上猶縣1∶20 萬地質圖。
通過對上猶縣地質調查與水文分析,進行如下操作:①編寫上猶縣災害調查表,從上猶縣災害調查表得: 上猶縣災害點500 個,通過GIS 平臺投影在上猶縣圖層。并對災害點進行500 m 的緩沖區(qū)分析,在緩沖區(qū)之外選取500 個非滑坡單元作為機器學習的訓練集與測試集。②基于GIS 平臺,分別從DEM 數(shù)據(jù)獲得:坡度、高程、坡向、水系、平面曲率、剖面曲率、濕度指數(shù)(TWI)、地形起伏度(SPI)、地表粗糙度(TRI)等柵格圖。③將ladsat5TM 數(shù)據(jù)band5 與band4 輸入arcgis 中,先對兩波段柵格整型數(shù)據(jù)轉化為浮點型,然后運用arcgis 的地圖代數(shù)中的柵格計算器進行(band5-band4)/(band5+band4)運算,獲取上猶縣植被指數(shù)(NDVI)。④并將上述影響因子進行重分類,如圖3(a)~圖3(h)所示,將坡度分為0°~5°、5~9°、9°~12°、12°~16°、16°~20°、20°~25°、 25°~30°、30°~37°、>37°;坡向分為平面(-1)、北(0°~22.5°)、東北(22.5°~67.5°)、東(67.5°~112.5°)、東南(112.5°~157.5°),南(157.5°~202.5°)西 南(202.5°~247.5°)、西(247.5°~292.5°)、西北(292.5°~337.5°)、北(337.5°~360°);平面曲率分為0~9、9~17、17~25、25~34、34~43、43~53、53~63、63~73、>73; 剖面曲率分為0~2、2~4、4~6、6~8、8~10、10~12、12~15、15~20、>20;地形起伏度分為-71~-26、-26~-15、-15~-9、-9~-3、-3~1、1~6、6~13、13~22、>22;地面粗糙度分為1~1.01、1.01~1.04、1.04~1.07、1.07 ~1.11、1.11 ~1.15、1.15 ~1.21、1.21 ~1.31、1.31 ~1.46、>1.46; 濕度指數(shù)分為2~6、6~9、9~21、21~35、>35; 植被指數(shù)劃分為-0.34~0.17、0.17~0.40、0.40~0.53、0.53~0.62、>0.62;⑤根據(jù)上猶縣道路與水系進行歐氏距離分析,并利用重分類將道路緩沖區(qū)分為0~200 m、200~400 m、400~600 m、600~800 m、800~1000 m、>1000 m; 將水系緩沖區(qū)分為0~100 m、100~200 m、 200~300 m、300~400 m、400~500 m、>500 m。⑥根據(jù)1∶20 萬上猶縣地質圖,提取地層分布與斷層,分別在Arcgis 中做出如圖3(k)~圖3(l)所示,上猶縣的地層由志留系、第四系、泥盆系、奧陶系、石炭系、三疊系與寒武系組成。對地質構造進行緩沖區(qū)分析,并進行重分類分為三類,分別為0~1 km、1~2 km、>2 km。其最終評價因子與災害點分布如圖3 所示。
由于進行易發(fā)性評價過程中要保證數(shù)據(jù)之間不相關或弱相關,能有效地提高模型的效率。因此利用GIS 的波段集統(tǒng)計,將上述評價因子柵格圖層進行相關性分析。如表1 所列,由于高程與植被指數(shù)、水系、道路相關性大于0.3,且粗糙度與植被指數(shù)、濕度指數(shù)的相關性大于0.3,所以除去高程與粗糙度保證評價因子之間的不相關或弱相關。選取坡度、坡向、平面曲率、剖面曲率、地形起伏度、濕度指數(shù)、植被指數(shù)、距道路距離、距水系距離,地層分布與距斷層距離等11 個上猶縣滑坡易發(fā)性評價因子。
對上猶縣滑坡進行500 m 緩沖區(qū)分析,在緩沖區(qū)以外選取500 個非滑坡單元。將滑坡單元、非滑坡單元與11 個評價因子圖層進行空間連接,得到滑坡與非滑坡單元在上述評價因子空間數(shù)據(jù),并將數(shù)據(jù)進行歸一化。然后將滑坡單元與非滑坡單元的空間數(shù)據(jù)分為80%的訓練集和20%的測試集,輸入RapidMiner Studio 進行邏輯回歸、 神經(jīng)網(wǎng)絡與支持向量機模型的訓練與測試。為了保證模型的精度,采用交互驗證法,挑選神經(jīng)網(wǎng)絡的訓練次數(shù)120、學習率0.2、動量為0.1,支持向量機與邏輯回歸模型采用默認參數(shù)。并將上猶縣柵格單元所在空間歸一化數(shù)據(jù)輸入訓練好的模型中,得出上猶縣柵格單元的易發(fā)性指數(shù),最后利用自然斷點法將上猶縣柵格的滑坡易發(fā)性指數(shù)進行重分類。如圖4 所示,將上猶縣分為高易發(fā)區(qū)、較高易發(fā)區(qū)、中易發(fā)區(qū)、較低易發(fā)區(qū)、低易發(fā)區(qū)。
統(tǒng)計上述模型易發(fā)性分區(qū)的災害數(shù)與柵格數(shù),并計算出占總滑坡數(shù)比例、占總柵格數(shù)比例、滑坡比率,得表2、 表3、 表4。如表2、 表3 與表4 所示:LR、ANN、SVM 的易發(fā)性分區(qū)的高易發(fā)區(qū)與較高易發(fā)區(qū)所包含的滑坡單元占比分別為80.6%、74.6%、91%,SVM 的高易發(fā)區(qū)所包含的災害占比最高,針對該區(qū)域進行重點滑坡災害治理,可以更好地控制上猶縣滑坡災害。高易發(fā)區(qū)的面積占比能很好地衡量在易發(fā)性分區(qū)的科學性,根據(jù)統(tǒng)計表顯示:LR、ANN、SVM 的高易發(fā)區(qū)占比分別為23.1%、22.8%、32.0%,都占比較小,符合上猶縣的實際情況?;卤嚷适腔卤壤c易發(fā)區(qū)占比的比值,衡量分類最優(yōu)化,保證高易發(fā)區(qū)與較高易發(fā)區(qū)的滑坡密度最大。由統(tǒng)計表顯示:LR、ANN、SVM 模型的高易發(fā)區(qū)的滑坡比率分別為2.458、2.438、2.381,都有較大的滑坡比率,而且隨著易發(fā)性分區(qū)從高到低滑坡比率逐漸減小,符合滑坡易發(fā)性從高易發(fā)區(qū)到低易發(fā)區(qū)的分區(qū)結果。
表1 波段集統(tǒng)計各因子相關性Table 1 Correlation of various factors in band set statistics
表2 邏輯回歸易發(fā)性統(tǒng)計表Table 2 Statistics of logistic regression susceptibility
表3 神經(jīng)網(wǎng)絡易發(fā)性統(tǒng)計表Table 3 Statistics of neural network susceptibility
表4 支持向量機易發(fā)性統(tǒng)計表Table 4 Statistics of SVM susceptibility
文章通過受試者工作曲線(ROC)進行模型的精度評價,ROC 的評定模型精度的標準用曲線下的面積(AUC)來衡量。AUC>0.5 時,當AUC 越接近1,其精度越高。AUC 在0.5~0.7 之間說明其精度差,AUC在0.7~0.9 之間精度較好,AUC在0.9~1 之間精度很好,AUC<0.5 則不符合客觀事實[22]。將訓練所得到的易發(fā)性指數(shù)與訓練數(shù)據(jù)一起輸入SPSS 中,通過ROC分析得如圖5,橫軸為特異性,縱軸為敏感性。如圖5所示,LR、ANN、SVM 的AUC 值分別為0.897、0.939、0.884,都大于0.8,證明模型有較好的評價精度。其中神經(jīng)網(wǎng)絡模型有最高的評價精度。
為了上猶縣柵格單元的安全性,引入滑坡發(fā)生的預警值來衡量滑坡的易發(fā)性,即要在滑坡發(fā)生之前進行預警,所以選取3 種模型評價上猶縣柵格單元易發(fā)性的最大值來當柵格發(fā)生滑坡的可能性,其公式為:
式中,ILSI為易發(fā)性指數(shù),MAX 取LR、ANN、SVM 易發(fā)性指數(shù)的最大值。將式(9)所得的易發(fā)性指數(shù)導入上猶縣柵格,得到基于LR-ANN-SVM 模型的上猶縣柵格單元易發(fā)性指數(shù)。最后將上猶縣柵格單元滑坡易發(fā)性指數(shù)進行重分類,得如圖4(d)所示的上猶縣易發(fā)性評價分區(qū)圖。
高易發(fā)區(qū)所包含的災害比例能反映模型評價的科學性,將更多的災害單元包含在高與較高易發(fā)區(qū)更方便政府部門的治理。根據(jù)表2、表3、表4、表5 統(tǒng)計可知,LR、ANN、SVM、LG-ANN-SVM 的高易發(fā)區(qū)與較高易發(fā)區(qū)所包含的災害占總災害比例分別為:80.6%、74.6%、91%、93.2%。結果顯示LR-ANN-SVM易發(fā)性分區(qū)的高易發(fā)區(qū)有更高的滑坡災害比例,低危險區(qū)災害比例占比為1%,此分區(qū)更適用于上猶縣風險性的實際應用,如圖5 ROC 曲線下的面積(AUC)所示,LR-ANN-SVM 的AUC 值為0.815,說明LRANN-SVM 評價模型有很好的評價精度。
表5 LR-ANN-SVM 易發(fā)性統(tǒng)計表Table 5 Statistics of LR-ANN-SVM susceptibility
綜合上述4 種模型易發(fā)性分區(qū)結果分析,其高危險區(qū)多在水系、縣城附近,以及泥盆系與志留系地層上,河水的沖刷能帶走壓坡腳的土之外,還會增加土壤的含水率,使滑動面的摩擦力減小,導致滑坡更容易產(chǎn)生??h城附近是人類工程活動導致植被的破壞,植被指數(shù)較低,并且城鎮(zhèn)附近道路較多,在修建的過程中,難免會削坡,導致坡度的增加,破壞土體的穩(wěn)定性,從而產(chǎn)生了滑坡。最終分類結果與上猶縣歷史滑坡實際情況相符。
1) 以上猶縣為研究區(qū),對滑坡編錄數(shù)據(jù)進行分析,利用相關性分析確定評價因子。通過LR、ANN 與SVM 模型進行上猶縣滑坡易發(fā)性評價。并利用受試者工作曲線(ROC)對上述模型進行精度評價,結果顯示LR、ANN、SVM 的AUC 值 分 別 為0.897、0.939、0.884,3 個模型都有較好的易發(fā)性評價精度。
2) 為了上猶縣柵格單元的安全性,引入預警值來對上猶縣滑坡易發(fā)性進行分析,將預警值取為上述3 個模型得出的上猶縣柵格易發(fā)性指數(shù)的最大值,結果顯示LR-ANN-SVM 模型易發(fā)評價的AUC 值為0.815,有很好的評價精度。從高易發(fā)區(qū)與較高易發(fā)區(qū)滑坡占比來看,LR、ANN、SVM、LR-ANN-SVM 高易發(fā)區(qū)與較高易發(fā)區(qū)所包含的滑坡占總滑坡比例為80.6%、74.6%、91%,93.2%,LR-ANN-SVM 所包含的災害比例最高,且低易發(fā)區(qū)災害點占比僅含1%。通過預警能很好地進行滑坡災害治理,減少上猶縣滑坡災害發(fā)生的概率。
3) 高易發(fā)區(qū)與較高易發(fā)區(qū)位于道路與城鎮(zhèn)附近,植被發(fā)育相對薄弱的地帶,并隨著水系呈帶狀分布。與上猶縣歷史滑坡的實際情況相符,可作為上猶縣災害治理的科學依據(jù)。