鎮域尺度下秦巴山區堆積層滑坡易發性不同單元評價性能對比研究

李澤芝，王新剛

（大陸動力學國家重點實驗室，西北大學地質學系，陜西 西安 710069）

秦巴山區地質構造強烈，地形起伏大，地貌形態與巖土體結構類型復雜多樣，巖體節理裂隙發育，風化嚴重，降雨豐富，人類工程活動活躍，導致區內災害頻發，尤其是以第四紀堆積層滑坡為主的地質災害（隱患）具有數量多、分布廣、密度大、頻次高等特點，所造成的危害嚴重，是中國地質災害的重災區之一（強菲等，2015；趙力行等，2020；劉偉等，2021）。因此，秦巴山區堆積層滑坡具有孕災條件復雜多樣、部分災害評價數據獲取難度大等特征（張世林，2020），如斜坡堆積層厚度、斜坡結構類型。

滑坡易發性評價模型建立流程為：區域滑坡數據及孕災條件編錄→評價因子提取與選擇→評價模型選取→模型參數推定→評價結果精度驗證與分析等步驟，可見評價模型是十分重要的一環，關乎評價結果的準確性與可信度（黃發明等，2022）。目前滑坡易發性評價模型主要分為知識驅動模型、物理驅動模型、數據驅動模型3 類，其中，數據驅動模型以區域基礎環境因子為數據輸入，通過機器學習模型、數理統計模型等得到評價結果（Saro Lee et al.，2015；張林梵等，2022）。據研究表明，機器學習模型的預測性能明顯高于數理統計模型與傳統的知識驅動模型，其中隨機森林模型在預測精度、提高泛化能力等方面較其他機器學習模型有著更好的表現，如Merghadi 等（2018）以北非米拉盆地為例，比較了基于隨機森林、梯度提升機、邏輯回歸、神經網絡和支持向量機5種評價模型下滑坡易發性的預測結果，得出隨機森林模型預測性能更高的結論。

目前，滑坡易發性評價時所使用區劃單元類型主要為柵格單元與斜坡單元（田述軍等，2019；常志璐等，2023）。柵格單元區劃在空間上破壞了評價區域斜坡的自然屬性，但劃分時操作簡單，在給單元賦值時有更強的唯一性，信息損失較少；斜坡單元劃分時相較于柵格更復雜，給單位賦值時可能造成信息損失，但能夠很好的反應區域滑坡易發性的地形分布情況，對某些基于斜坡的特征因子有更好的表征，如斜坡結構、堆積層厚度等。

近年來，隨著1∶10 000 地質災害風險評價的開展，評價區域的尺度也轉向鎮域級，意味著對地質災害風險區的劃分有了相對縣域級評價更高的精度要求，主要表現在對全域內每個斜坡孕災地質條件的掌握（陜西省自然資源廳，2021）。不同地質環境背景下滑坡的主要孕災因子不盡相同（黃潤秋等，2007），如黃土地區由于土體的水敏性、灌溉等因素所致滑坡較為普遍（段釗等，2018），故土地利用類型為黃土地區滑坡易發性評價的重要因子；而秦巴山區多以堆積層碎石土-基巖的二元結構斜坡為主（范立民等，2004；唐睿旋，2017）層間的滲透性差異使該地區以堆積層滑坡為主，故堆積層厚度為堆積層滑坡易發性評價的重要因子。

基于上述分析，筆者選取秦巴山區柞水縣小嶺鎮作為典型研究區，結合堆積層滑坡特點，采取柵格單元、斜坡單元兩種單元類型，運用隨機森林模型方法對該區域進行了滑坡易發性評價與精度分析，得到鎮域滑坡空間概率，并探討易發性評價在兩種劃分方法下的差異，研究成果對當地政府部門的減災防災措施部署等具有重要意義。

1 研究區及模型特征因子分析

1.1 研究區概況

小嶺鎮地處秦巴山區腹地（圖1），區內地形陡峭、溝壑縱橫，其地質構造受東西展布的山陽—鳳鎮等深大斷裂影響較大，其滑坡孕災條件也較為復雜：區域內斜坡以二元結構的陡坡為主，豐富的礦產資源與稀少的平緩地帶，不可避免的形成了區內削坡建房頻繁、切坡修路密集、毀林造地普遍、礦山開發持續的現象，這些現象導致了斜坡開挖形成臨空面致使開挖處剪應力集中，加之堆積層大都為滲透性強的松散碎石土，礦區附近更有礦渣堆積使得堆積層自重增大，孔隙水壓力增強，導致堆積層滑坡在研究區內極易發育，經統計發現（表1），研究區地質災害以堆積層滑坡為主，占區內滑坡的92.86%，主要分布在道路兩側的陡坡處（圖2）。

表1 研究區滑坡規模分類Tab.1 Classification of landslide scale in study area

圖1 研究區位置圖Fig.1 Location of the study area

圖2 研究區已有堆積層滑坡分布圖Fig.2 The distribution of alluvial landslide in the study area

1.2 隨機森林模型簡介

隨機森林是一種集成Bagging 算法，可用于分類問題與回歸問題（方匡南等，2011）。在分類問題中，旨在生成眾多決策樹，并通過對建模數據集的樣本觀測和特征變量分別進行隨機抽樣，每次抽樣結果均為一棵樹，且每棵樹都會生成符合自身屬性的規則和分類結果，而森林最終集成所有決策樹的規則和分類結果，實現隨機森林算法的分類。其主要特點在于數據采樣隨機、特征選取隨機。

2 模型數據制備

2.1 數據來源及用途

將鎮域孕災因子數據賦予評價單元作為隨機森林分類模型的樣本數據，其中孕災因子為模型特征因子，區域滑坡環境特征因子數據源主要包括遙感影像數據、表面數字高程模型、地質圖、野外調查數據等，其分類及用途見表2。

表2 數據類型及用途Tab.2 Data types and uses

2.2 評價單元

①柵格單元：研究區柵格共計4 576 590 個柵格；提取出滑坡周界內所占柵格，共計6 192 個柵格單元（圖3）。

圖3 研究區柵格單元劃分Fig.3 Grid division of Study Area

②斜坡單元：對研究區DEM 用水文分析模型進行斜坡單元初步劃分，再結合研究區衛星影像圖與10 m 等高線進行人工校正（谷天峰等，2013；劉彬等，2021），最終將研究區劃分為729 個斜坡單元（圖4），單個面積概況見表3，對其中滑坡所在斜坡單元進行標記，共計28 個斜坡單元。

表3 斜坡單元面積概況Tab.3 Overview of slope unit area

圖4 研究區斜坡單元劃分Fig.4 Slope division of study area

2.3 特征因子

滑坡的發生是一個由多種因素作用下產生的非常復雜的非線性過程，要通過機器學習模型對全區發生滑坡的概率進行分級，特征因子的選取至關重要（Ahmed et al.，2016；吳潤澤等，2021），當根據研究區滑坡的孕災特點，有針對性的選取。在前人研究基礎上，結合區域孕災特征，初步選9 個孕災因子作為模型特征因子用于建模，并依據特征屬性，將特征因子進行離散化后做出屬性空間劃分（圖5），具體如下：

圖5 研究區滑坡地質災害易發性評價指標因子Fig.5 Index factors of landslide geological hazard Susceptibility assessment in the study area

①地形因子：坡度、坡高。

滑坡的主要動力來源于自重，坡度則控制著斜坡內部的應力分布，是滑坡的重要影響因素。將研究區斜坡坡度劃分為：0～15°、15°～25°、25°～35°、35°～45°和＞45°等5 個屬性空間。

在巖土類型大致相似的地質環境下，坡體高度影響著坡體自重大小以及斜坡范圍內地形變化，進而影響滑坡的規模與范圍（郭芳芳等，2008）。先將DEM重采樣為5 m×5 m 柵格，再利用分區統計工具得到坡體高度范圍，將其分為0～20 m、20～50 m、50～100 m、100～300 m、＞300 m 等5 個屬性空間。

②斜坡自身特征：斜坡坡面形態、斜坡結構、斜坡堆積層厚度。

斜坡形態以其剖面曲率劃分，剖面曲率定義為曲面在最陡斜坡方向上的曲率，當曲率值＜-0.5 時值判斷為凹型坡，曲率值＞0.5 時判斷為凸型坡，曲率值在-0.5～0.5 之間時判斷為直線型坡。

據前人研究，斜坡結構在堆積層滑坡的發生中也是一個重要的影響因素（邢林嘯等，2012；賈琳等，2021），從研究區已發生的堆積層滑坡來看，39.13%都發生在順向斜坡中，而逆向斜坡中無已發生滑坡，故根據野外調查數據，將區內斜坡結構劃分為：順向斜坡、斜向斜坡、橫向斜坡、近水平狀斜坡、切向斜坡、逆向斜坡等6 個屬性空間。

由于區內滑坡大都為堆積層滑坡，堆積層厚度在一定程度上決定了土地利用類型、邊坡開挖程度、地下水位埋深，對堆積層滑坡起著重要啟發作用，但對于這一指標，目前缺乏經濟又精準的測取方法，本研究方法如下：對每個坡體隨機抽點使用插桿進行測量得到平均數，結合依托項目的工程鉆探數據對調查結果進行修正并類比到全區，賦值給每個斜坡，將該因子劃分為：0～1 m、1～3 m、3～6 m、6～10 m、＞10 m等5 個屬性空間。

③人類工程相關特征：距道路距離、距河流距離、距礦區距離。

區內平緩地帶稀少，道路開通處往往伴隨著邊坡開挖，形成陡立臨空面；河流對坡腳的侵蝕、依河而居帶來的切坡建房；礦山的開挖、震爆、不合理的礦渣堆積，這些因素都為滑坡的發生提供了有利條件。距河流、道路、礦區的距離均可反應人類工程活動的強弱程度，其各自屬性區間劃分如圖5 所示。

④地質條件：距斷裂距離。

區內受山陽—鳳鎮等深大斷裂影響較大，該斷裂呈近東西展布，經鳳鎮橫貫研究區，該斷裂對區內表層巖體破壞嚴重，加之研究區內節理裂隙極為發育，極大地破壞了巖土體的完整性，增加了滑坡的發生概率。對研究區地質圖進行矢量化，提取研究區內斷裂為矢量文件，對其做緩沖分析，劃分為0～100 m、100～400 m、400～700 m、700～1000 m、＞1000 m 等5 個屬性區。將上述特征因子賦值給兩種劃分單元：①柵格單元：在Arcgis10.7 中，將柵格轉化為點數據，用“多值提取至點”工具將各因子值提取至點，與初步選取的9 個因子構成矩陣A4576590×9，作為原始樣本1。②斜坡單元：除斜坡結構等因子在每個斜坡單元有唯一值外，對于坡面形態等離散型變量，取其眾數賦值；對于坡度等連續型變量，取其平均數賦值，729 個斜坡單元與初步選取的9 個因子構成矩陣B729×9，作為原始樣本2。

對兩種劃分單元下的原始樣本分別進行頻率比分析，得到各因子屬性區間的頻率比值。頻率比值的大小可反應該屬性空間滑坡發育的優勢程度（郭子正等，2019），若某屬性空間評價單元樣本數過少（如斜坡單元中近水平斜坡數量），則在選擇優勢空間時，剔除該屬性空間，最終統計得各因子的滑坡發育優勢空間（表4）。

表4 兩種評價單元下各因子的滑坡發育優勢空間Tab.4 Dominant space of landslide development of each factor under two evaluation units

隨機森林模型對特征因子的相關性有一定要求，因子之間相關性過高會降低模型的訓練精度。對研究區特征因子首先進行正態性檢驗，以確定相關性分析的計算方式。若數據滿足正態分布，則選用Pearson相關性分析；若數據不滿足正態分布，則用Speraman相關性分析；針對有序變量數據，則選用Kendall’staub 等級進行相關性分析。

由表5 可知，兩種評價單元下各因子顯著性P 值均為0.000***，水平呈現顯著性，拒絕原假設，因此數據不滿足正態分布。由于數據為有序數據，選用Kendall’s tau-b 等級進行相關性分析（表6、表7）。由結果可知，柵格單元下，研究區距河流距離與距道路距離特征因子的相關性系數為0.653＞0.3；斜坡單元下，研究區距河流距離與距道路距離特征因子的相關性系數為0.606＞0.3，說明兩種評價單元下道路與河流之間的相關性都較高。由于研究區堆積層滑坡主要分布在道路兩側，故剔除河流距離因子，最終采用其他剩余8 個因子進行建模評價。

表5 特征因子數據正態性檢驗結果Tab.5 Characteristic factor data Normality test results

表6 特征因子Kendall’s tau-b 等級相關系數矩陣（柵格單元）Tab.6 Characteristic factor Kendall’s tau-b rank correlation coefficient matrix (grid units)

表7 特征因子Kendall's tau-b 等級相關系數矩陣（斜坡單元）Tab.7 Characteristic factor Kendall's tau-brank correlation coefficient matrix (slope units)

3 評價及結果對比分析

3.1 評價結果

對于柵格單元，將4 576 590 個柵格單元轉化為連接有各因子數據的點后，將6 192 個滑坡柵格點作為滑坡樣本并將其易發性賦值為1，占總樣本的0.13%，在滑坡范圍外選取3 倍于滑坡樣本的點作為非滑坡樣本，以減少滑坡與非滑坡樣本之間的數量不平衡與空間相關性，并將其易發性賦值為0，兩類數據合并后以7∶3 的比例劃分為訓練集與測試集，在R 語言的RF 軟件包中構建隨機森林模型，將合并后數據輸入進行訓練。訓練結束后，將整個研究區點數據輸入，得到所有點的易發性值。將點數據轉為柵格，利用自然間斷法劃分出極低、低、中、高、極高等5 個易發區（圖6）。對于斜坡單元，選取28 個滑坡所在的斜坡單元作為滑坡樣本并將其易發性賦值為1，占總樣本的3.8%，在滑坡范圍外選取3 倍于滑坡樣本的點作為非滑坡樣本，并將其易發性賦值為0。剩余操作與柵格單元類似，結果見圖6。由圖6、表8 可知：①兩種評價單元下滑坡易發性分區大體相近，滑坡高、極高易發區主要位于鎮域北部礦區附近、西北部道路切坡密集區域、西南部盤山公路切坡密集道路附近，表明研究區內滑坡的發育與道路、礦區的影響有較強的相關性，這也與實際情況相符合。②柵格單元下預測的極低、低、中、高和極高易發區占研究區全區比例為：66.76%、19.50%、8.99%、3.04%和1.71%，已知滑坡樣本在高、極高易發樣本量中占95.49%；斜坡單元下預測的極低、低、中、高、極高易發區占研究區全區比例分別為61.92%、17.70%、10.09%、5.25% 和5.04%，已知滑坡樣本在高、極高易發樣本量中占96.42%。③柵格單元下高、極高易發區頻率比之和占總頻率比值為0.93，略低于斜坡單元下比值0.98，可知斜坡單元下的滑坡易發性可更好的判斷已知滑坡，且斜坡單元下的評價可很好的表現滑坡的地形分布特點。

表8 柵格單元與斜坡單元下評價結果頻率比Tab.8 Frequency ratio of evaluation results under grid unit and slope unit

圖6 研究區滑坡易發性區劃Fig.6 The division of landslide susceptibility in the study area

由模型得到兩種評價單元下各特征因子的重要性比例（圖7），可知道路開挖、斜坡結構對滑坡的發生起主要作用，事實上，現有滑坡多發生在道路兩旁及順向結構斜坡中。

圖7 各特征因子貢獻值Fig.7 Contribution value of each characteristic factor

3.2 模型精度及評價結果對比

結合前人研究（劉堅等，2018；鄭迎凱等，2020），筆者選取ROC 曲線來進行模型精度驗證。將全部樣本滑坡易發性概率遞減排序，閾值從1 至0 變更，計算各閾值下對應的（FPR、TPR）數值對，將數值對繪制于直角坐標系中。橫坐標FPR 代表非滑坡預測為滑坡的數量/真實非滑坡數量；縱坐標TPR 代表滑坡預測為滑坡的數量/真實滑坡數量。曲線下面積（AUC）越接近1，則說明模型預測性能越高。根據兩種評價單元下預測結果，繪制各自ROC 曲線（圖8），可以看出，基于柵格單元、斜坡單元下隨機森林模型的AUC值分別為0.864、0.921。表明5 m 高精度柵格單元下隨機森林模型的預測精度小于基于斜坡單元下模型的預測精度。

圖8 ROC 曲線Fig.8 ROC curve

滑坡的發生相對于整個區域是一個小概率事件，整個區域內滑坡易發性概率在預測模型精度較高的情況下，均值越低，表明預測的不確定性越低；標準差越低，表明預測區分度更優，由表9 可知，柵格單元下預測滑坡易發性概率均值為0.10 小于斜坡單元下滑坡易發性概率均值0.13；柵格單元下預測滑坡易發性概率標準差為0.13，小于坡單元下滑坡易發性概率標準差0.18。說明柵格單元作為評價單元，預測結果具有較差的區分度，斜坡單元反之。

表9 不同評價單元下易發性概率均值與標準差Tab.9 Mean and standard deviation of probability of Susceptibility under different evaluation units

4 結論

（1）結合研究區孕災特征，選取了坡度、坡高、坡面形態、斜坡結構類型、堆積層厚度、道路、礦區、區域地質斷裂等8 個孕災因子作為模型特征因子，柵格單元、斜坡單元兩種評價單元的預測結果，均發現模型主要根據道路、斜坡結構等因素區分滑坡、非滑坡區域。據各因子頻率比值可知，研究區內，在距道路100 m 內的順向淺層堆積層斜坡中，最易發生滑坡。

（2）在柵格單元、斜坡單元兩種評價單元下，研究區95.49%、96.42% 的滑坡落在隨機森林模型預測分區出的高、極高易發區內；通過AUC 值、易發性概率均值、易發性概率標準差對比表明：隨機森林是一種可靠、全面的滑坡易發性預測模型，評價單元為斜坡單元時預測精度更高、準確性更強。

（3）滑坡易發性評價中，柵格單元與斜坡單元在評價結果上，兩種評價單元的預測結果都有良好的表現。柵格單元的樣本數遠高于斜坡單元，意味著樣本數據中滑坡樣本比例減小，樣本的不均衡性增大，使得隨機森林模型的預測結果性能變弱，但以柵格單元為評價單元可減少特征因子的數據損失，在災害防治具體空間部署上有著更精細的參考；斜坡單元在特征因子賦值時，會造成部分數據損失，但其能更好體現滑坡災害的地形空間屬性，相較于柵格單元其樣本的不均衡性更小，使得其基于隨機森林模型的預測性能高于柵格單元。