李宇新,鄧念東*,馬建全,周 陽,崔陽陽
(1.西安科技大學(xué)地質(zhì)與環(huán)境學(xué)院,西安 710054;2.陜西省地質(zhì)調(diào)查院,西安 710054;3.陜西省水工環(huán)地質(zhì)調(diào)查中心,西安 710068)
隨著社會(huì)經(jīng)濟(jì)的快速發(fā)展,地質(zhì)環(huán)境逐漸惡化,衍生出一系列地質(zhì)災(zāi)害。漢中市漢臺(tái)區(qū)地處陜南漢中盆地中心,由于土地資源不合理利用、開墾山坡等人類工程活動(dòng),該區(qū)域地質(zhì)環(huán)境遭到嚴(yán)重破壞,特別是滑坡對(duì)人類生命財(cái)產(chǎn)造成巨大威脅。科學(xué)高效地進(jìn)行區(qū)域滑坡易發(fā)性評(píng)價(jià)是防災(zāi)減災(zāi)的有效措施,也是近些年來滑坡預(yù)防與土地規(guī)劃的熱點(diǎn)。
從20世紀(jì)80年代起,中外學(xué)者對(duì)滑坡易發(fā)性展開研究,其核心是評(píng)價(jià)指標(biāo)和評(píng)價(jià)模型的確定。根據(jù)前人研究,評(píng)價(jià)指標(biāo)主要根據(jù)區(qū)域地質(zhì)環(huán)境條件進(jìn)行選擇,尚未達(dá)成統(tǒng)一共識(shí)。評(píng)價(jià)方法主要可以劃分為確定性模型和非確定性模型[1]。確定性模型主要依靠專家經(jīng)驗(yàn)確定各類成災(zāi)因素的權(quán)重值后進(jìn)行疊加,取得了較好的預(yù)測(cè)效果,但過度依賴經(jīng)驗(yàn)判斷以及結(jié)論缺乏繼承性。隨著統(tǒng)計(jì)學(xué)習(xí)理論研究和計(jì)算機(jī)科學(xué)的發(fā)展,非確定性模型被廣泛應(yīng)用到滑坡易發(fā)性評(píng)價(jià)中,主要包括信息量[2-3]、確定性系數(shù)法[4]、證據(jù)權(quán)法[5]以及一系列數(shù)據(jù)挖掘與機(jī)器學(xué)習(xí)方法,比如邏輯回歸[6]、決策樹[7]、人工神經(jīng)網(wǎng)絡(luò)[8]、支持向量機(jī)[9]、隨機(jī)森林[10]、樸素貝葉斯[11]等。其中,以決策樹算法為基礎(chǔ)的一系列集成模型也被應(yīng)用至滑坡易發(fā)性評(píng)價(jià)中,如J48算法[12]、隨機(jī)森林模型、旋轉(zhuǎn)森林模型[13]、邏輯回歸樹模型[14]等,該類模型主要以決策樹算法或多種分類器為基礎(chǔ),采取不同的特征選擇和剪枝方式進(jìn)行模型構(gòu)建,提高模型分類正確率,同時(shí)一定程度上減少訓(xùn)練樣本過擬合。研究表明,集成模型往往比單一機(jī)器學(xué)習(xí)模型預(yù)測(cè)能力更高。為進(jìn)一步探索集成模型在我國滑坡易發(fā)性評(píng)價(jià)的適用性,分析漢中市漢臺(tái)區(qū)滑坡分布趨勢(shì)以及獲得較高精度的滑坡易發(fā)性圖,本文采用邏輯回歸樹和旋轉(zhuǎn)森林模型進(jìn)行該區(qū)域的滑坡易發(fā)性評(píng)價(jià)。
邏輯回歸樹模型(logistic model tree, LMT)是結(jié)合決策樹學(xué)習(xí)理論與邏輯回歸的集成模型[14],使用LogitBoost算法在樹的節(jié)點(diǎn)上生成邏輯回歸模型,同時(shí)使用分類與回歸樹算法對(duì)樹進(jìn)行剪枝。LMT模型通過交叉驗(yàn)證進(jìn)行大量迭代,采用遞增邏輯回歸計(jì)算分類LC,即是否為滑坡。最后使用線性邏輯回歸計(jì)算得到每個(gè)葉節(jié)點(diǎn)的后驗(yàn)概率值P,作為該模型滑坡易發(fā)性評(píng)價(jià)結(jié)果。LC、ρ的計(jì)算公式分別為
(1)
(2)
式中:βi為因子的系數(shù);xi為各評(píng)價(jià)因子;n為因子數(shù)量;D為分類數(shù)量,D=1,2。
旋轉(zhuǎn)森林(rotation forest, ROF)以決策樹為基分類器,構(gòu)建包含滑坡因子屬性以及類別標(biāo)簽的樣本集x=[x1,x2,…,xn]T,選擇L個(gè)基分類器{D1,D2,…,DL},隨機(jī)將特征集F劃分為K個(gè)不相交的子集,每個(gè)子集均有M個(gè)特征。對(duì)特征子集進(jìn)行主成分分析,旋轉(zhuǎn)換后得到大小均為M×1的向量子集,同時(shí)存儲(chǔ)主成分系數(shù)組成Ri矩陣[式(3)]。根據(jù)式(4)計(jì)算Di決策樹分類器的平均概率值[15-16],即代表評(píng)價(jià)單元滑坡發(fā)生的概率。
Ri=
(3)
(4)
研究區(qū)位于陜西省西南部漢中盆地中心,地理坐標(biāo)為東經(jīng)106°51′40″~107°10′25″,北緯33°01′40″~33°22′00″,總面積為556 km2。屬北亞熱帶濕潤季風(fēng)氣候區(qū),氣溫受地形影響差異明顯,降雨主要集中在6—9月。地形地貌從南到北依次為平原區(qū)、丘陵區(qū)、中低山區(qū)。其中,平原區(qū)地勢(shì)平坦,地質(zhì)環(huán)境問題較少;丘陵區(qū)膨脹土發(fā)育,發(fā)育膨脹土滑坡災(zāi)害;中低山區(qū)山坡陡峻,山背狹長,陡坡處基巖裸露,緩坡殘坡積土覆蓋,在降雨充沛時(shí)易發(fā)生殘坡積層滑坡。區(qū)內(nèi)出露地層為震旦系、寒武系、石炭系及第四系。受構(gòu)造作用影響,巖層褶皺變形強(qiáng)烈,巖石破碎,巖體內(nèi)部結(jié)合力較差。巖體以花崗巖、碳酸鹽巖、變質(zhì)巖為主。較堅(jiān)硬-較軟變質(zhì)巖易風(fēng)化產(chǎn)生破壞變形而產(chǎn)生滑坡。研究區(qū)人類工程活動(dòng)劇烈,例如開墾山坡、道路建設(shè)開挖邊坡、礦山棄渣不合理堆放等,為滑坡發(fā)育提供了條件。
在前人進(jìn)行區(qū)內(nèi)1∶50 000地質(zhì)災(zāi)害詳查的基礎(chǔ)上,筆者通過遙感影像解譯與實(shí)地調(diào)查進(jìn)行對(duì)比,共圈定40處滑坡。為更加便捷高效地獲取滑坡相關(guān)屬性信息,在劃分滑坡周界與對(duì)比核查屬性信息后,通過ArcGIS軟件提取其質(zhì)心作為滑坡屬性點(diǎn),得到區(qū)內(nèi)滑坡編目圖(圖1),為后續(xù)評(píng)價(jià)研究奠定基礎(chǔ)。
圖1 研究區(qū)位置與滑坡編錄圖
主要通過以下方式獲取研究區(qū)相關(guān)數(shù)據(jù)。
(1)從“地理空間數(shù)據(jù)云”平臺(tái)中ASTER GDEM獲取研究區(qū)數(shù)字高程模型(digital elevation model,DEM)數(shù)據(jù),獲得坡度、坡向、平面曲率、剖面曲率、歸一化植被利用指數(shù)(normalized difference vegetation index,NDVI)和地形濕度指數(shù)(topographic wetness index,TWI)因子。
(2)從1∶50 000地質(zhì)圖矢量化生成地形地貌和巖土體類型數(shù)據(jù)。
(3)根據(jù)研究區(qū)氣象觀測(cè)站數(shù)據(jù)生成研究區(qū)年均降雨量因子。
(4)通過Bigemap地圖軟件獲得研究區(qū)水系、道路矢量數(shù)據(jù)。
結(jié)合研究區(qū)地質(zhì)災(zāi)害詳查報(bào)告及相關(guān)資料,研究區(qū)內(nèi)滑坡受中低山地形貌控制,由于溝谷發(fā)育、斜坡高陡、降雨集中,以及強(qiáng)烈的人類工程活動(dòng)為滑坡的發(fā)育提供條件。因此選取高程、坡度、坡向、年均降雨量、地形地貌、巖土體類型、TWI、NDVI、距水系距離、距道路距離、平面曲率和剖面曲率共12類因子進(jìn)行評(píng)價(jià)。
根據(jù)湯國安經(jīng)驗(yàn)公式和研究區(qū)地質(zhì)圖比例尺大小,選擇30 m分辨率柵格單元進(jìn)行評(píng)價(jià),共計(jì)614 117個(gè)柵格[17]。依據(jù)上述數(shù)據(jù)源生成12類因子數(shù)據(jù),對(duì)連續(xù)型因子分別采用Jenks自然間斷法和等間距法進(jìn)行分級(jí);離散型因子根據(jù)二級(jí)因子類型進(jìn)行劃分,各類因子專題圖層見圖2。采用頻率比對(duì)各分級(jí)下滑坡進(jìn)行統(tǒng)計(jì),表1顯示區(qū)內(nèi)滑坡主要分布于高程在672~964 m、巖土體類型為石英片巖、坡度為7.37°~17.10°、平面曲率為-4.34~-0.50、年均降雨量在740~760 mm、坡向?yàn)闁|南以及NDVI為0.34~0.43的區(qū)域。
表1 各因子分級(jí)下滑坡頻率比
評(píng)價(jià)模型的數(shù)據(jù)之間往往存在一定的多重共線性關(guān)系,當(dāng)多重共線性較大時(shí),會(huì)嚴(yán)重影響到模型分類的辨識(shí)能力[18]。為有效避免這一影響,本文選擇方差膨脹因子(variance inflation factor,VIF)與容忍度(tolerance,TOL)進(jìn)行分析。VIF反映多元線性回歸模型中多重共線性嚴(yán)重程度,表示回歸系數(shù)估計(jì)量的方差與假設(shè)自變量間不線性相關(guān)時(shí)方差相比的比值,其中TOL為VIF的倒數(shù)。一般認(rèn)為VIF>10或TOL<0.1時(shí),數(shù)據(jù)之間存在著嚴(yán)重的共線性,需要進(jìn)行部分剔除與修正[19]。根據(jù)表2結(jié)果,所選12類因子之間共線性較弱。
因子貢獻(xiàn)率表征著與研究區(qū)滑坡發(fā)生的相關(guān)程度,一方面有利于滑坡發(fā)生機(jī)理分析;另一方面,貢獻(xiàn)率為0的因子對(duì)模型的分類預(yù)測(cè)不僅沒有幫助,反而容易會(huì)造成數(shù)據(jù)冗余,降低預(yù)測(cè)精度[20]。相關(guān)屬性評(píng)估(correlation attribute evaluation,CAE)通過計(jì)算影響因子與標(biāo)簽屬性之間的相關(guān)性來評(píng)估因子的重要程度,其結(jié)果包括平均貢獻(xiàn)率(average merit,AM)和標(biāo)準(zhǔn)差(standard deviation,SD),根據(jù)貢獻(xiàn)率降序排列見表2。結(jié)果顯示,12類因子對(duì)研究區(qū)滑坡發(fā)生均有促進(jìn)作用,其中貢獻(xiàn)率排列前三的影響因子分別為地形地貌(AM=0.258)、平面曲率(AM=0.223)以及巖土體類型(AM=0.225)。
表2 影響因子共線性分析與貢獻(xiàn)率
從滑坡范圍以外區(qū)域隨機(jī)提取等量的非滑坡點(diǎn)作為負(fù)樣本數(shù)據(jù)進(jìn)行模型數(shù)據(jù)構(gòu)建,按照7∶3比例隨機(jī)劃分為訓(xùn)練集(56處)和驗(yàn)證集(24處),并提取12類影響因子屬性值。本文通過Weka3.8軟件進(jìn)行LMT和ROF模型的構(gòu)建。通過十倍交叉驗(yàn)證進(jìn)行訓(xùn)練,然后代入驗(yàn)證集進(jìn)行測(cè)試,最后將整個(gè)研究區(qū)的屬性集代入上述兩種模型,得到分別基于LMT和ROF模型的滑坡易發(fā)性指數(shù)(landslide susceptibility index, LSI),通過ArcGIS的重分類工具將LSI值劃分為5類:極低易發(fā)區(qū)、低易發(fā)區(qū)、中等易發(fā)區(qū)、高易發(fā)區(qū)和極高易發(fā)區(qū),得到兩種模型的滑坡易發(fā)性分區(qū)圖(圖3)。
圖3 研究區(qū)滑坡易發(fā)性圖
圖3表明LMT與ROF模型分區(qū)結(jié)果基本一致,研究區(qū)滑坡高-極高易發(fā)區(qū)主要分布于中部丘陵區(qū)和北部中低山區(qū)。中部丘陵區(qū)巖性主要為第四系黏土,受降雨影響膨脹土變形加劇,黏性土滑坡最為發(fā)育。北部中低山區(qū)滑坡高易發(fā)區(qū)沿線狀分布,這是由于西北側(cè)人工擴(kuò)建公路,頻繁開挖坡腳,以及石英礦、磷礦等礦山開采程度高,一系列人類工程活動(dòng)破壞地質(zhì)環(huán)境、改變地形地貌所造成。中低山東側(cè)地勢(shì)高差較大,歷史滑坡主要沿河谷分布,地層巖性為震旦系千枚巖夾灰?guī)r,軟弱夾層多,較破碎、易風(fēng)化,在降雨充沛條件下易發(fā)生滑坡。研究區(qū)滑坡低-極低易發(fā)區(qū)主要分布于南側(cè)平原區(qū),該區(qū)域?yàn)閺V闊的一、二級(jí)階地區(qū)域,地勢(shì)平坦,地層巖性以第四系黏土為主,極少有滑坡發(fā)育。
進(jìn)行模型精度驗(yàn)證與比較是分類結(jié)果對(duì)滑坡易發(fā)性分區(qū)可靠程度驗(yàn)證的重要步驟。本文采用接受者工作特性曲線(receiver operating characteristic curve,ROC)與其線下面積(area under curve,AUC),以及滑坡頻率比對(duì)LMT模型和ROF模型進(jìn)行評(píng)估。ROC曲線于20世紀(jì)90年代起開始廣泛應(yīng)用于數(shù)據(jù)挖掘、機(jī)器學(xué)習(xí)等分類模型的評(píng)估,它以敏感度(即實(shí)際為滑坡,預(yù)測(cè)為滑坡)為縱坐標(biāo)、1-特異性(即實(shí)際為非滑坡,預(yù)測(cè)為滑坡)為橫坐標(biāo),通過動(dòng)態(tài)分類閾值避免界限值對(duì)結(jié)果的影響[21]。AUC取值范圍為0~1,值越大代表模型分類效果越好,通常認(rèn)為AUC>0.7時(shí),表明分類預(yù)測(cè)能力較強(qiáng)。圖4和圖5表明,ROF模型在訓(xùn)練集的正確率(77.4%)較LMT模型(75.5%)相比更高,同時(shí)驗(yàn)證集的預(yù)測(cè)率結(jié)果表明ROF模型(93.1%)優(yōu)于LMT模型(84.0%)。圖6表明兩者模型各易發(fā)性等級(jí)下,區(qū)內(nèi)滑坡集中分布在高-極高易發(fā)區(qū),而低-極低易發(fā)區(qū)很少或無滑坡分布,證明分區(qū)結(jié)果符合歷史滑坡分布規(guī)律。其中ROF模型高-極高易發(fā)區(qū)分布有37處歷史滑坡,多于LMT模型(31處);并且ROF模型低易發(fā)區(qū)的滑坡數(shù)量為0,少于LMT模型(4處)。將歷史滑坡密度與分區(qū)等級(jí)面積占比的比值作為頻率比,用來對(duì)比不同模型預(yù)測(cè)滑坡發(fā)生的敏感性。由圖7可看出,兩種模型頻率比總體呈上升趨勢(shì),ROF模型極高易發(fā)區(qū)的頻率比(6.52)高于LMT模型(2.07),說明ROF模型對(duì)滑坡分布更為敏感,預(yù)測(cè)結(jié)果更可靠。
圖4 訓(xùn)練集ROC曲線
圖5 驗(yàn)證集ROC曲線
圖6 分區(qū)結(jié)果與歷史滑坡數(shù)量對(duì)比
圖7 各易發(fā)性等級(jí)頻率比
機(jī)器學(xué)習(xí)模型的發(fā)展使得滑坡易發(fā)性評(píng)價(jià)更加快捷與高效,但仍存在如下不足。
(1)研究結(jié)果主要為對(duì)比LMT與ROF模型在中國滑坡易發(fā)性評(píng)價(jià)的適用性,后期宜展開模型參數(shù)調(diào)整對(duì)評(píng)價(jià)結(jié)果變化的動(dòng)態(tài)研究。
(2)本次研究?jī)烧吣P万?yàn)證集AUC值均高于訓(xùn)練集,這可能由于樣本數(shù)量因素以及樣本數(shù)據(jù)劃分的隨機(jī)性導(dǎo)致,后期宜探討滑坡樣本規(guī)模與評(píng)價(jià)模型適應(yīng)性的關(guān)系。
(3)由于滑坡發(fā)生機(jī)理復(fù)雜,滑坡易發(fā)性評(píng)價(jià)的因子選擇仍存在一定片面性或主觀性,后期宜對(duì)滑坡影響因子的選取與優(yōu)化進(jìn)行補(bǔ)充研究。
通過對(duì)漢臺(tái)區(qū)進(jìn)行基于LMT和ROF模型的滑坡易發(fā)性評(píng)價(jià),得到以下結(jié)論。
(1)通過相關(guān)資料與野外調(diào)查,選取研究區(qū)高程、坡度、坡向、年均降雨量、地形地貌、巖土體類型、TWI、NDVI、距水系距離、距道路距離、平面曲率和剖面曲率共12類因子進(jìn)行評(píng)價(jià),并通過CAE和VIF表明,因子屬性適合分類模型的構(gòu)建,并且各類因子對(duì)滑坡發(fā)育均有影響,其中關(guān)系最密切的是地形地貌、平面曲率、巖土體類型、NDVI和距道路距離。根據(jù)滑坡易發(fā)性分區(qū)圖,研究區(qū)滑坡主要為分布于中低山區(qū)和丘陵區(qū)黏性土滑坡。
(2)將數(shù)據(jù)代入Weka3.8軟件中構(gòu)建LMT和ROF模型,生成基于兩種模型的滑坡易發(fā)性分區(qū)圖。結(jié)果顯示兩種模型預(yù)測(cè)分布趨勢(shì)基本一致,且符合歷史滑坡分布規(guī)律。
(3)通過ROC曲線、AUC值以及頻率比進(jìn)行模型精度驗(yàn)證與比較,結(jié)果顯示兩種模型均適用于研究區(qū)滑坡易發(fā)性評(píng)價(jià)。ROF模型的訓(xùn)練集與驗(yàn)證集AUC分別為0.774和0.931,均高于LMT模型的0.755和0.840;滑坡頻率比顯示,ROF模型較LMT模型對(duì)研究區(qū)滑坡易發(fā)性更敏感,分區(qū)結(jié)果更為精確。