王 毅， 方志策， 牛瑞卿

(中國地質(zhì)大學(武漢),湖北 武漢 430074)

中國三峽庫區(qū)地質(zhì)結(jié)構(gòu)特殊,地形地貌復(fù)雜,人類工程活動較為頻繁,滑坡等地質(zhì)災(zāi)害廣泛發(fā)育,對人類生命安全構(gòu)成嚴重威脅,同時也給社會經(jīng)濟帶來巨大損失[1]?；聻?zāi)害易發(fā)性預(yù)測是通過分析歷史和現(xiàn)存滑坡災(zāi)害的空間分布,考慮研究區(qū)地形地貌、水文氣象、基礎(chǔ)地質(zhì)等各類斜坡不穩(wěn)定因素,對滑坡空間分布及其發(fā)生概率進行預(yù)測,能夠為滑坡防治提供關(guān)鍵的科學依據(jù)。

根據(jù)理論基礎(chǔ)的不同,滑坡災(zāi)害易發(fā)性預(yù)測方法分為定性和定量方法兩大類[2]。前者主要依賴專家經(jīng)驗和研究人員對滑坡災(zāi)害的評估與分析能力；后者主要通過分析滑坡空間分布與滑坡影響因子之間的關(guān)系來預(yù)測可能發(fā)生滑坡的位置,包括確定性方法和數(shù)據(jù)驅(qū)動法。確定性方法通?；跓o限斜坡模型開展滑坡易發(fā)性預(yù)測。Montgomery等[3]結(jié)合無限斜坡穩(wěn)定模型與坡面水文模型建立了SHLSTAB模型,并將其用于田納西河谷滑坡易發(fā)性預(yù)測。Pack等[4]將穩(wěn)定態(tài)水文學理論融入無限斜坡穩(wěn)定模型中,提出了穩(wěn)定指數(shù)制圖模型,并將其應(yīng)用于溫哥華流域易發(fā)性預(yù)測。在國內(nèi),蘭恒星等[5]在穩(wěn)定指數(shù)制圖模型計算過程中增加動水壓力項,將其應(yīng)用于云南小江流域滑坡易發(fā)性預(yù)測。

數(shù)據(jù)驅(qū)動法在近些年來備受研究人員關(guān)注,包括頻率比法[6]、證據(jù)權(quán)法[7]、邏輯回歸(Logistic Regression,LR)[8]、支持向量機[9]和人工神經(jīng)網(wǎng)絡(luò)[10]等。深度學習作為當前人工智能領(lǐng)域的研究熱點,在各個領(lǐng)域引起了廣泛關(guān)注。其中,卷積神經(jīng)網(wǎng)絡(luò)(Convolutional Neural Networks,CNN)與循環(huán)神經(jīng)網(wǎng)絡(luò)(Recurrent Neural Networks,RNN)作為深度學習中兩個最具代表性的方法,都具有強大的表征學習能力,已經(jīng)被成功應(yīng)用于圖像處理、場景理解和自動駕駛等領(lǐng)域中[11-13]。隨著三峽工程的建設(shè),長期的地質(zhì)災(zāi)害工作積累了大量地質(zhì)資料與災(zāi)害數(shù)據(jù),這類數(shù)據(jù)通常具有多源性、異構(gòu)性、時空性、模糊性和非線性等特征,如何有效、充分地利用這些數(shù)據(jù)開展災(zāi)害管理防治工作是亟待解決的難題。此外,研究人員已將單個深度學習方法成功應(yīng)用于滑坡災(zāi)害易發(fā)性預(yù)測,但其集成能力還未被挖掘。鑒于此,本文以長江三峽庫區(qū)——秭歸—巴東段(湖北省境內(nèi))為研究區(qū),利用Stacking集成學習技術(shù)融合CNN與RNN來構(gòu)建易發(fā)性評價模型,并將其應(yīng)用于研究區(qū)滑坡災(zāi)害預(yù)測。

1 數(shù)據(jù)準備與分析

1.1 研究區(qū)背景

研究區(qū)位于湖北境內(nèi)(圖1),總面積達423.14 km2,地跨東經(jīng)110°18′～110°52′、北緯30°54′～31°03′,最高海拔達到1 968.1 m。屬亞熱帶季風氣候區(qū),年平均降雨量為1 100 mm,大部分降雨集中在5—9月,占年降雨量的70%。研究區(qū)地層發(fā)育較完整,除缺失泥盆系下統(tǒng)、志留系和石炭系上統(tǒng)和白堊系的大部分及第三系外,從震旦系至第四系皆有出露,總體上具有自東向西漸新展布的規(guī)律，巖性以灰?guī)r、白云巖等碳酸鹽巖和砂巖、泥巖、頁巖等碎屑巖為主。研究區(qū)大地構(gòu)造處于揚子準地臺的中西部,是川東褶皺帶、大巴山弧和淮陽山字形西翼反射弧的交匯區(qū)。褶皺和斷裂對研究區(qū)的構(gòu)造影響較大,其中褶皺發(fā)育有黃陵背斜、秭歸向斜等，斷裂主要有仙女山斷裂、九畹溪斷裂和香爐坪斷裂。

圖1 三峽庫區(qū)秭歸—巴東段滑坡災(zāi)害分布圖Fig.1 Distribution map of landslide locations in Zigui-Badong section of the Three Gorges Reservoir Region

近年來,人類工程活動對研究區(qū)的影響日益增大,包括移民城鎮(zhèn)建設(shè)、水庫建設(shè)、毀林開荒和礦山開采活動等。研究區(qū)主要的地質(zhì)災(zāi)害類型包括滑坡、崩塌、塌岸和泥石流等,滑坡是區(qū)內(nèi)發(fā)育數(shù)量最多的地質(zhì)災(zāi)害,占其總數(shù)的83%。

1.2 滑坡數(shù)據(jù)編錄

近年來,三峽庫區(qū)庫水位的周期性波動破壞了斜坡的穩(wěn)定結(jié)構(gòu),導致了滑坡災(zāi)害的頻繁發(fā)生。本文通過實地調(diào)查、Google Earth解譯以及歷史編錄數(shù)據(jù)統(tǒng)計得到研究區(qū)滑坡編錄數(shù)據(jù)(圖1)。研究區(qū)共有196處滑坡,其中滑坡最大的面積為1.51 km2,最小的為172 m2。典型滑坡主要有黃土滑坡、黃蠟石滑坡、范家坪滑坡、新灘滑坡和千將坪滑坡等。

1.3 滑坡評價因子

評價因子的選取對于滑坡易發(fā)性評價至關(guān)重要。目前,研究人員在三峽地區(qū)開展了一系列滑坡災(zāi)害易發(fā)性預(yù)測研究[14-18]。因此,本文通過野外調(diào)查研究及已有相關(guān)研究,綜合考慮研究區(qū)的地形地貌、基礎(chǔ)地質(zhì)、水文條件等因素,最終選取20個滑坡影響因子(圖2)。為了保證這些因子的空間一致性,將這些因子圖層重采樣至相同空間分辨率。在此基礎(chǔ)上,根據(jù)歷史文獻及專家經(jīng)驗對連續(xù)型圖層進行重分類,再對所有圖層子類別進行賦值。相關(guān)因子數(shù)據(jù)獲取源如下：數(shù)字高程模型(digital elevation model,DEM)數(shù)據(jù)來自地理空間數(shù)據(jù)云ASTER GDEM v2版本數(shù)據(jù)(http://www.gscloud.cn/)；坡向、坡度、流域面積、流域坡度、曲率和斜坡形態(tài)因子通過ArcGIS軟件從DEM數(shù)據(jù)中提?。坏匦挝恢弥笖?shù)(Terrain Position Index,TPI)、地形粗糙指數(shù)(Terrain Ruggedness Index,TRI)、地形表面曲率(Terrain Surface Convexity,TSC)、地形表面紋理(Terrain Surface Texture,TST)和地形濕度指數(shù)(Topographic Wetness Index,TWI)因子通過SAGA軟件從DEM數(shù)據(jù)中提?。粠r性和斷層距離因子提取于湖北省地質(zhì)局提供的地質(zhì)圖。降雨量因子來自2003—2010年雨量站點數(shù)據(jù),震級因子來自1970年以來研究區(qū)及附近區(qū)域地震資料。通過地形圖和遙感影像數(shù)據(jù)提取研究區(qū)水系網(wǎng),再利用歐氏距離工具獲得水系距離因子；利用支持向量機對2010年獲取的Landsat-7衛(wèi)星影像進行分類獲得土地利用因子；利用ENVI軟件波段計算工具獲得植被歸一化指數(shù)(Normalized Difference Vegetation Index,NDVI)和歸一化水指數(shù)(Normalized Difference Water Index,NDWI)。

圖2 滑坡影響因子專題圖Fig.2 Thematic maps of landslide influencing factors

2 滑坡易發(fā)性評價模型

2.1 卷積神經(jīng)網(wǎng)絡(luò)

CNN最初是由LeCun等人針對手寫數(shù)字識別而提出的一種特殊的神經(jīng)網(wǎng)絡(luò)模型[19],可用于處理具有已知網(wǎng)格拓撲結(jié)構(gòu)的數(shù)據(jù)。CNN除了擁有與普通神經(jīng)網(wǎng)絡(luò)共有的輸入層和輸出層之外,其結(jié)構(gòu)還包含卷積層、池化層和全連接層(圖3)。

圖3 典型CNN網(wǎng)絡(luò)結(jié)構(gòu)圖[20]Fig.3 Typical network structure of CNN

如圖3所示,卷積層中的卷積核可看作規(guī)則形狀的濾波器,卷積窗口從輸入數(shù)據(jù)的最左上方開始,按照從左往右、從上往下的順序依次滑動,將其固定的權(quán)重與窗口內(nèi)的數(shù)據(jù)做內(nèi)積運算。卷積核在整個輸入數(shù)據(jù)上完成一次卷積運算后會輸出卷積特征圖,即卷積操作所提取出的初級特征。以二維卷積為例,其計算公式如下：

(1)

式中：k為卷積核的個數(shù)；Cj表示第j個卷積核的輸出；f表示非線性激活函數(shù)；i為卷積操作的位置；xi表示卷積窗口對應(yīng)的輸入數(shù)據(jù)；wj和bj分別表示權(quán)重和偏置。在卷積操作中,輸入數(shù)據(jù)通常為多維數(shù)組,而卷積核則是通過學習算法不斷調(diào)整固定大小的多維數(shù)組,這些多維數(shù)組被稱為張量。

池化層中池化操作不僅能夠通過降低卷積層輸出特征圖的維度來提高計算效率,還能使提取的特征保持平移不變性[21]。最常見的池化操作為最大池化,其計算公式如下：

(2)

全連接層也是CNN隱藏層中的一種,一般位于多個卷積層和池化層之后,其神經(jīng)單元都與上一層的神經(jīng)元相連。全連接層可看作為特殊的分類器,其目的是將卷積、池化操作提取出來的高維特征映射到低維特征空間。

2.2 循環(huán)神經(jīng)網(wǎng)絡(luò)

RNN考慮了輸入數(shù)據(jù)之間的序列關(guān)系,并將“時間序列”引入了神經(jīng)網(wǎng)絡(luò)模型[22]。相比傳統(tǒng)人工神經(jīng)網(wǎng)絡(luò)模型,RNN隱藏層中點的信息都能傳遞到下一個時間節(jié)點,因此數(shù)據(jù)中所蘊含的信息能夠一直傳遞下去。本文構(gòu)建了基于因子重要性排序的序列數(shù)據(jù),利用因子選擇方法對滑坡評價因子進行重要性賦值與排序,再將其輸入到RNN模型中。如圖4所示。假設(shè)x為輸入的序列數(shù)據(jù),時間步長t時的輸出層為yt,它由當前節(jié)點的輸入和上一節(jié)點共同決定,其計算公式如下：

圖4 典型RNN結(jié)構(gòu)圖[23]Fig.4 Typical structure of RNN

ht=σ(Whxt+Uhht-1+bh)

(3)

yt=σ(Wyht+by)

(4)

式中：σ為激活函數(shù)；Wh和Uh為權(quán)重矩陣；b為偏置向量。

2.3 Stacking集成方法

Stacking集成是一種流行的異構(gòu)集成學習方法,它通過使用元學習器將不同的基分類器組合在一起以獲得更準確的預(yù)測結(jié)果[24]。與常用的Bagging和Boosting兩種集成方法不同,Stacking用于集成不同類型的模型,而不是多個相同類型的模型。圖5展示了Stacking集成的一般結(jié)構(gòu),該集成算法首先通過交叉驗證的方式,使用訓練基分類器未使用的樣本來產(chǎn)生元學習器的訓練樣本。以5折交叉驗證為例,初始訓練集被劃分為5個大小相同的集合D1,D2,…,D5。令其中一個集合為初級測試集,剩余的4個集合為初級訓練集訓練初級分類器,再利用該初級分類器對初級測試集進行預(yù)測,預(yù)測值作為元學習器訓練集的一部分,重復(fù)5次得到元學習器的訓練集。最后,利用該訓練集訓練元學習器并進行預(yù)測。Stacking集成方法可以通過初級學習步驟共同估計所有初級分類器的誤差,并使用元學習步驟減少預(yù)測殘差。

圖5 Stacking集成結(jié)構(gòu)圖Fig.5 Structure of the Stacking ensemble

2.4 模型評價指標

本文選用總體精度(Overall Accuracy,OA)、精確率(Precision)、召回率(Recall)、馬修斯相關(guān)系數(shù)(Matthews Correlation Coefficient,MCC)和一致性系數(shù)(Index of Agreement,IOA)這5個統(tǒng)計學指標評價模型。其中,精確率表示預(yù)測為滑坡的樣本中有多少是真實的滑坡樣本,召回率表示樣本集中的滑坡樣本有多少被準確預(yù)測[25]。MCC是一個較為均衡的評價指標,能夠描述實際類別與預(yù)測類別之間的相關(guān)系數(shù)[26]。IOA是一種模型預(yù)測誤差程度的標準化度量,它是由均方誤差和潛在誤差之比計算得到[27]。上述評價指標的公式如下所示。這5個評價指標的取值范圍在0～1之間,其中取值越接近于1表示對應(yīng)的模型性能越好。本文使用的另外一個評價指標為受試者工作特征(Receiver Operating Characteristic,ROC)曲線,它通過設(shè)定多個不同的臨界值計算一系列敏感性和特異性,并基于這兩個值繪制成曲線[8]。其曲線下面積(Area Under the Curve,AUC)可用來衡量模型的優(yōu)劣。AUC值越高,表明該模型性能更優(yōu)異。

(5)

(6)

(7)

(8)

(9)

3 實驗結(jié)果與分析

3.1 滑坡評價因子分析

根據(jù)歷史研究及專家經(jīng)驗,初步選取的滑坡評價因子之間難免存在統(tǒng)計學上的共線性關(guān)系,這會導致易發(fā)性模型無法準確地分析評價因子與滑坡之間的真實關(guān)系。此外,影響滑坡因素眾多,評價因子中存在的冗余信息可能會影響易發(fā)性結(jié)果的精度。因此,本文從共線性分析及重要性分析兩個角度對滑坡評價因子進行定量評價與篩選。本研究區(qū)評價因子多重共線性分析結(jié)果如表1所示,其中VIF和TOL分別表示方差膨脹系數(shù)和容忍度[28],結(jié)果顯示所有的滑坡影響因子之間不存在多重共線性。

表1 滑坡影響因子多重共線性分析結(jié)果Table 1 Multicollinearity analysis results of landslide influencing factors

此外,本文選用隨機邏輯回歸方法來定量分析滑坡評價因子的重要性。該算法對訓練數(shù)據(jù)重采樣并對每個樣本應(yīng)用L1懲罰邏輯回歸模型以選擇重要特征[29]。在多次迭代重復(fù)該過程后,如果始終選擇某個特征,則該特征對于建模至關(guān)重要,而如果從未選擇該特征則無關(guān)緊要。滑坡評價因子的重要性排序如圖6所示,其中TWI相對其他因子而言對滑坡貢獻最大,而坡向、斜坡形態(tài)和TPI因子重要性統(tǒng)計值為0,可視為無關(guān)因子而被剔除。因此,將剩余的17個評價因子用于后續(xù)的模型構(gòu)建。

圖6 滑坡影響因子重要性Fig.6 Importance of landslide influencing factors

3.2 模型構(gòu)建

在進行滑坡易發(fā)性預(yù)測之前,需要選取訓練樣本和測試樣本用于構(gòu)建滑坡易發(fā)性模型,且訓練集與測試集均由正樣本(滑坡樣本)和負樣本(非滑坡樣本)共同組成。由文中1.2節(jié)可知,研究區(qū)共有196處滑坡、25 380個滑坡柵格單元,從非滑坡區(qū)域隨機選取25 380 個柵格單元作為負樣本。為了探索滑坡易發(fā)性模型在滑坡樣本非常有限的條件下的預(yù)測能力,從正樣本和負樣本中只隨機選取10%的數(shù)據(jù)進行模型訓練,而剩余的90%樣本用于模型測試。

Stacking集成方法通過使用訓練集集成CNN和RNN。首先,選擇訓練集的10%作為驗證集來訓練模型，CNN和RNN模型基于梯度下降算法迭代訓練直到損失函數(shù)收斂(圖7)。隨著訓練次數(shù)的增加,損失值減小直到達到一個低水平,這表明訓練過程是令人滿意的。然后,LR結(jié)合CNN和RNN的預(yù)測輸出最終的分類結(jié)果。所使用的深度學習模型結(jié)構(gòu)如下：CNN包含一個有著20個卷積核的卷積層、一個最大池化層和一個包含50個神經(jīng)元的全連接層；RNN僅包含一個具有17個隱藏節(jié)點的循環(huán)層。此外,在實驗中,相關(guān)的超參數(shù)是通過先前的研究和驗證集的微調(diào)確定的[25,30],相關(guān)超參數(shù)如表2所示。

圖7 訓練過程模型損失值變化Fig.7 Change of loss value in training process model

表2 模型超參數(shù)設(shè)置Table 2 Hyperparameter settings of the models

3.3 滑坡易發(fā)性預(yù)測圖

為了評估Stacking方法的有效性,選擇了CNN、RNN和LR這三種模型進行比較。本文使用自然斷點法[6，9，17]將研究區(qū)劃分為5個易發(fā)性等級,分別為極低、低、中、高和極高易發(fā)區(qū)。四種模型的滑坡易發(fā)性預(yù)測分區(qū)圖如圖8所示，可以看出四個分區(qū)圖在空間分布上具有相似性，例如極高易發(fā)區(qū)主要分布在長江沿岸,大多數(shù)滑坡也都位于滑坡易發(fā)性極高區(qū)域。此外,四個分區(qū)圖也存在少許差異，例如LR的中易發(fā)區(qū)分布(黃色區(qū)域)最多。

圖8 滑坡災(zāi)害易發(fā)性預(yù)測分區(qū)圖Fig.8 Prediction zoning map of landslide susceptibility

為了定量分析易發(fā)性分區(qū)圖,本報告使用滑坡密度來衡量結(jié)果的可靠性和準確性?；旅芏葹橐装l(fā)性分區(qū)中滑坡比例與柵格個數(shù)比例的比值[31]。如表3所示,四種模型對應(yīng)的易發(fā)性分區(qū)圖中滑坡密度最高的為極高易發(fā)區(qū),其次為高易發(fā)區(qū)、中易發(fā)區(qū)、低易發(fā)區(qū)和極低易發(fā)區(qū)。這與實際情況較為吻合,即滑坡往往更容易出現(xiàn)在極高易發(fā)區(qū),而極低易發(fā)區(qū)滑坡出現(xiàn)的概率非常低。同時,從表3可以發(fā)現(xiàn),Stacking的極高易發(fā)區(qū)滑坡密度在4種模型中最高,其次為CNN、RNN和LR，這說明Stacking模型預(yù)測的滑坡災(zāi)害高危區(qū)能夠?qū)v史滑坡包含在內(nèi),也證明了Stacking對研究區(qū)的易發(fā)性預(yù)測更為合理與準確。

表3 模型各易發(fā)性分區(qū)滑坡密度Table 3 Landslide density of different susceptibility classes

3.4 模型比較與驗證

四種模型的精度如表4所示，可以看出Stacking的精度最高,OA高達83%,比深度學習基分類器高出0.87%～1.61%；LR的精度最低,OA為80.11%。模型ROC曲線如圖9所示，可以看出無論是基于訓練集還是測試集,Stacking都能取得最高的AUC值,其次為CNN、RNN和LR。

圖9 模型ROC曲線Fig.9 ROC curves of different models

表4 模型精度評價Table 4 Performance of different models

4 結(jié)論

本文以長江三峽庫首區(qū)——秭歸—巴東段為例,選取高程、坡度、坡向、曲率、巖性等20個滑坡影響因子,利用Stacking集成策略融合CNN和RNN模型開展

滑坡災(zāi)害易發(fā)性預(yù)測。研究表明,Stacking集成學習可以充分利用基分類器的優(yōu)點,擁有較好泛化能力和較高的預(yù)測精度。實驗結(jié)果表明,Stacking集成方法在使用極其有限的樣本進行建模時,能獲得更為可靠的滑坡災(zāi)害易發(fā)性預(yù)測圖,其預(yù)測總體精度為83%,皆高于CNN、RNN和LR。綜上所述,本文不僅證明了Stacking集成在滑坡易發(fā)性預(yù)測有著巨大的應(yīng)用潛力,也為深度學習集成模型的應(yīng)用打下了基礎(chǔ)。此外,深度學習在滑坡災(zāi)害易發(fā)性預(yù)測中的應(yīng)用也存在超參數(shù)尋優(yōu)較復(fù)雜的問題,解決這些問題并探索更有效的深度學習模型將是未來的研究目標。

致謝：感謝三峽庫區(qū)地質(zhì)災(zāi)害防治工作指揮部提供的相關(guān)實驗數(shù)據(jù)和資料。