基于突變理論的湟水河流域崩滑易發性評價

李 彬，周金喜，吳鈳橋，馬建全，祁漢文

（1.青海省青藏高原北部地質過程與礦產資源重點實驗室，青海 西寧 810000；2.青海省有色第一地質勘查院（青海省隱伏礦勘查重點實驗室），青海 西寧 810000；3.青海省有色第四地質勘查院，青海 西寧 810000；4.西安科技大學地質與環境學院，陜西 西安 710000）

0 引言

隨著社會的不斷發展，人類工程活動不斷加強，地質災害帶來的不良影響越來越多，如何科學準確的評判地質災害的危險性十分重要。崩滑災害的發生極其復雜，受地形地貌、巖性組合、地質構造、降雨、地震和人類活動等多種因素的影響[1]。2014年宋盛淵等[2]采用突變級數理論的方法對巖桑樹水電站庫區發育的潛在單體滑坡進行危險性評價，取得了較好的效果；尚志海等[3]通過對梅州市地質災害傷亡人數、直接經濟損失和發生頻率的統計，引入突變理論的突變級數法，建立了地質災害災度的評估指標體系和等級劃分標準，并對梅州市地質災害災度進行了評估；2020年，劉曉宇等[4]利用突變理論選取巖性、風化程度、坡度等7 個指標，對西部山區威脅公路的10 個典型滑坡進行危險性評價，評價結果與現場調查結果具有較高的吻合度。目前，突變理論在災害評價中的應用，僅限于單體災害危險程度的研究或是某一地區災情程度的研究上，對于區域易發性評價研究中尚未見相關報導[5]。本文旨在創建一種新的基于突變理論的崩滑災害易發性評價方法，開發一種基于突變級數理論的崩滑災害易發性評價平臺。區別于前人研究成果的是對整個研究區域柵格單元進行突變級數賦值，進而對整個研究區崩滑災害易發性進行評價，而不是單單對單體災害的評價。

1 突變級數理論基本原理及易發性評價平臺開發

1.1 突變級數法基本原理

將研究區內的評價因子按照倒立的樹枝狀結構進行多層次分解是突變級數法的基本思想[6]。按照“互補”或“非互補”原則，將最底層的控災因子進行無量綱處理（歸一化），將歸一化結果代入相應的突變模型中，并計算出該層的突變級數。然后逐層向上計算各層的突變級數，并根據最高層的突變級數將所研究的系統分級[7 -8]。

（1）建立逐層結構模型

建立逐層結構模型是突變級數法的關鍵，要充分考慮各層因子之間的關系，逐層分解，選擇相對應的突變函數，按先主后次排列成倒立的樹枝狀結構。常用的分解形式如圖1 所示。圖中，x為狀態變量；u、v、w、t為控制變量。

圖1 突變模型系統Fig.1 Mutation model system

（2）評價指標的無量綱化處理

各參評因子的取值范圍和單位量綱均不相同，無法對各參評因子進行比較。因此，就有必要將各參評因子轉變為0～1 之間的無量綱數據。值得注意的是，所有參評因子數據均應該轉化為“越大越好”型數據，對“越小越好”型數據可采用倒數法進行標準化處理。對于正指標，采用式（1）進行標準化；對于負指標，采用式（2）進行標準化[2]。

式中：xi——評價指標值；

xmax、xmin——各評價指標的最大值和最小值。

（3）突變級數的計算

各參評因子無量綱處理后就可以選擇對應的初等突變類型（表1）。

表1 初等突變類型Table 1 Primary mutation types

根據不同模型的歸一化公式和最底層因子的無量綱數據，可計算出對應因子的突變級數值。在計算時，無相關性因子取最小值，有相關性因子取平均值。逐層向上計算突變級數，最終求出總突變級數。突變級數法流程見圖2。

圖2 突變級數法流程圖Fig.2 Flow chart of catastrophe theory procedure

1.2 基于突變級數理論的易發性評價平臺開發

開發基于突變級數理論模型的易發性評價平臺采用的是MATLAB 語言，它是當今科研領域最常用的應用軟件之一，具有強大的符號運算、矩陣計算及數據處理等功能[9]。MATLAB 有包羅萬象的工具箱，“草稿紙式”的編程語言，易學易用，用戶不僅可以調用其內部函數進行計算，還可以根據自己的算法進行擴展編程[9]。對突變級數理論模型進行編程使其程序化，最終開發出基于突變級數理論模型的崩滑災害易發性評價平臺，其平臺開發流程圖（圖3）。

圖3 基于突變級數理論的易發性平臺開發流程圖Fig.3 Development flow chart of susceptibility platform based on catastrophe series theory

2 研究區概況及地質災害影響因素

2.1 研究區概況

湟水流域位于青藏高原與黃土高原的過渡地帶，總體為高原地貌。地勢起伏較大，西北高、東南低，亦有南北高、中間低的特點（圖4）。最高海拔位于西北部的湟水干流源區，約4 350 m，最低海拔位于東南部的湟水干流河谷，約1 650 m。根據流域內實際，又可細分為中高山區、丘陵區和河谷區。由于湟水流域地形地貌以溝壑縱橫的丘陵為主，氣候氣象具有干旱少雨但多集中于汛期的特點，出露巖性以黃土、泥質巖等軟弱巖土體多見，植被覆蓋較為稀疏且森林覆蓋率低，斷裂、褶皺等構造發育，因而湟水流域地質環境脆弱，水土流失嚴重，崩塌、滑坡等地質災害多發易發；加之，流域內人口密集、人類工程活動強烈，對地質環境的改造明顯和突出，人為因素導致的崩塌、滑坡等也越來越多。因此，湟水流域崩滑災害多發，特別是中下游地區的丘陵區是青海省乃至全國地質災害高易發區。

圖4 湟水河流域崩滑災害分布圖Fig.4 Distribution map of landslide disasters in the Huangshui River Basin

2.2 因子選取及建立逐層機構模型

本文參考國內外研究學者經常選取的崩滑災害影響因子，同時考慮所選評價因子容易量化程度[10]，最終選擇具有代表意義的12 個因子作為崩滑災害易發性評價指標。選取因子分別為：工程地質巖組、距離城鎮、平面曲率、高程、剖面曲率、距離水系、地形地貌、距離公路、坡度、距離斷層、地形濕度指數（TWI）、多年平均降水量。

將12 個控災因子作為最底層因子指標進行多層次分解，先主后次排列成倒立的樹枝狀結構。逐層向下分解，直至分解到將災害點危險度作為目標層，以地質因素和環境因素兩個變量為準則層，地貌、地形、工程地質巖組、地質構造、水文氣象、人類工程活動6 個變量作為中間層，高程、坡度、平面曲率、剖面曲率、工程地質巖組、地形濕度指數（TWI）、距離斷層、年均降水量、地貌類型、距離水系、距離公路、距離城鎮12 個控災因子作為最底層變量，建立逐層結構模型。其中，目標層災害點危險度即災害點突變級數，各最底層變量即為控災因子，各層之間根據變量的不同選用不同的突變模型。表2 為湟水河流域崩滑災害危險度評價逐層結構模型。

表2 湟水河流域崩滑災害危險度評價逐層結構模型Table 2 Layer-by-layer structure model for landslide hazard assessment in the Huangshui River Basin

3 研究區崩滑災害評價

3.1 評價指標無量綱化

由于最底層評價指標之間的取值范圍和單位量綱均不相同，各指標之間無法進行比較，需將其轉變為0～1 之間的無量綱數據。利用歸一化原則，將湟水河流域內1 171 處崩滑災害進行無量綱處理。表3 為評價因子無量綱處理結果表。

表3 評價因子無量綱處理結果表Table 3 Dimensionless processing results of evaluation factors

3.2 突變級數的計算

在對流域內1 171 處崩滑災害點的各類因子進行無量綱化后，利用突變模型的突變勢函數計算各層評價指標的突變級數。按照表2 的突變模型，其中間層、準則層及目標層的突變級數計算公式為：

中間層：

根據式（3）—（11）計算得到湟水河流域崩滑災害點的突變級數，即各災害點危險度（表4）。

表4 湟水河流域崩滑災害點危險度計算結果表Table 4 Calculation results of landslide hazard at various points in the Huangshui River Basin

3.3 崩滑災害危險性評價

對于崩滑災害危險性等級的劃分，尚未形成統一的標準。因此，選用ArcGIS 重分類自然間斷點法將研究區內1 171 處崩滑災害點危險度進行分級，結合研究區崩滑災害易發性實際特征，將崩滑災害危險性劃分為輕度危險、中度危險、重度危險和極度危險4 個等級，崩滑災害危險性評價指標劃分標準見表5。

表5 崩滑災害危險性評價指標劃分標準Table 5 Division standard of landslide hazard evaluation index

對表4 中研究區1 171 處崩滑災害點突變級數按照劃分標準表5 進行分類，其評價結果如表6 所示。

表6 崩滑災害危險性評價結果Table 6 Risk assessment results of rockfall disasters

根據表6 危險性評價結果可知，研究區內1 171 處崩滑災害中輕度危險點75 處、中度危險點201 處、重度危險點476 處、極度危險點419 處。根據評價結果，研究區內災害點以重度和極度危險為主，與實際災情狀況一致。

按照突變級數法可以快速評價單體災害危險性，其優點是操作簡單，評價結果可靠，但其缺點是只能進行單體災害危險性分析，不能進行全流域易發性評價。為了對湟水河流域全域易發性進行評價，必須開發一套基于突變級數理論的易發性平臺，既滿足單體災害快速且準確率高的優點，也要能滿足全流域易發性評價的要求。因此，利用MATLAB 編程，對突變級數理論模型進行升級，在單體災害危險性評價基礎上開發一套基于突變級數理論模型的易發性評價平臺。

4 突變級數理論易發性平臺開發

4.1 平臺設計流程

（1）柵格單元賦值

通過突變級數計算獲得了湟水河流域內1 171 處崩滑災害的危險度，將崩滑災害各指標中間層、準則層及目標層突變級數值導入ArcGIS，利用ArcGIS 柵格轉出工具將數據轉出為ASCII 格式。湟水河流域總面積為16 100.65 km2，按照30 m×30 m 柵格單元劃分，可以將研究區劃分為17 889 612 個柵格單元，利用MATLAB編程，對每個柵格單元進行突變級數賦值。該方法的優點是將整個流域離散化，具有比常規方法精度高、操作簡單、評價結果與實際情況符合程度高等特點。具體步驟如下：

①構建評價因子柵格圖。（加載研究區評價底圖）

②對所有評價因子進行歸一化處理，得到最底層D（運用式（1）（2）對選取的高程、坡度等12 個控災因子進行歸一化處理，使其可以在同一維度內進行對比，得到12 個控災因子的歸一化結果D）。

③根據論文中表2 中逐層結構模型及突變模型，最底層D的各個層進行柵格計算并疊加得到中間層C（選取不同的突變模型，運用式（3）-（8）對D進行柵格計算，得到地貌、地形、工程地質巖組、地質構造、水文氣象、人類工程活動6 個中間層的突變級數C）。

④根據論文中表2 中逐層結構模型及突變模型，中間層C的各個層進行柵格計算并疊加得到中間層B（選取蝴蝶型和尖點型突變模型，運用式（9）（10），對C進行柵格疊加計算，得到地質和環境2 個準則層的突變級數B）。

⑤根據論文中表2 中逐層結構模型及突變模型，準則層B的各個層進行柵格計算并疊加得到目標層A，以此得到整個流域的災害危險度突變級數（選取尖點型突變模型，運用式（11）對B進行柵格疊加計算，得到目標層災害點危險度的突變級數A，即得到整個流域的災害柵格單元突變級數）。

⑥根據論文中表5，將計算得到的突變級數柵格圖進行重分類，得到易發性分區圖5（基于表5 給出的災害危險性評價指標劃分標準，對流域內災害危險度進行重分類并繪圖，得到流域災害易發性分區圖，見圖5）。

圖5 湟水河流域崩滑災害易發分區圖Fig.5 Landslide susceptibility zoning map in the Huangshui River Basin

（2）預定義平臺界面及屬性

新建一個平臺、定義平臺的結構、對加載數據定義到函數句柄、定義輸出及平臺的調用句柄、定義及更新平臺的調用句柄，輸出界面見圖6。

圖6 預定義平臺界面及屬性Fig.6 Predefined platform interface and properties

（3）讀取數據并繪圖

（以“D1地貌”按鈕為例）讀取數據、剔除異常值、保存數據的文本格式、構建句柄定義的坐標區窗口、展示圖像、設置顏色、保存數據、將更新后的數據和文本數據加載到句柄中、保存句柄到指定對象中、彈出確認提示，輸出界面見圖7。

圖7 讀取數據并繪圖Fig.7 Data retrieval and plotting

（4）數據歸一化并繪圖

（以D1地貌為例“歸一化”按鈕）恢復對象中的句柄數據、定義“地貌歸一化”值、構建句柄定義的坐標區、展示圖像、設置顏色、設置透明度、設置圖名、保存數據指定對象句柄中。設置背景色，選擇對應突變函數，對所有因子進行歸一化處理，對D1—D12重復進行步驟,輸出界面見圖8。

圖8 湟水河流域柵格單元突變級數計算圖Fig.8 Calculation of catastrophe series for grid units in the Huangshui River Basin

（5）突變級數計算并繪圖

讀取各因子歸一化數據、分層計算輸出目標值A、構建坐標區、目標值A歸一化、展示A歸一化值的圖像、設置顏色、設置透明度、保存數據到指定對象句柄。

（6）計算結果輸出

調用對象中保存得句柄信息、創建文件“YF.txt”、將dm 中的文本數據轉換為字符串數組、將HC 中所有行數據輸入YF.txt、將Ag 所有行數據輸入YF.txt、保存結束，輸出結果。

4.2 湟水河流域危險性評價結果

在地質災害評價中地質災害易發分區是一項重要內容，地質災害易發分區是根據災害的易發程度對研究區進行劃分等級[11-13]。根據平臺計算結果，將“YF.txt”文件導入ArcGIS，讀取柵格單元突變級數，按照自然間斷分級法將研究區柵格單元突變級數分為五個等級，即湟水河流域危險指數。分別是極低、低、中等、高和極高，得到基于突變級數理論的湟水河流域崩滑災害易發分區圖和易發分區統計表。圖5 為湟水河流域崩滑災害易發分區圖。

根據表7 統計結果可知，湟水河流域崩滑災害的極高易發區、高易發區主要分布于湟水河南北兩岸以及一級支流兩側的低山丘陵區，面積分別為4 661.79 km2和5 511.44 km2，占整個研究區總面積的28.95%和34.23%；中易發區分布于湟水河南北兩岸以及各大支溝兩側的構造侵蝕中山區，面積占研究區總面積的23.99 %；低易發區、極低易發區主要分布在湟水河河谷平原地帶及全流域中高山區，面積占比分別占總面積的11.80 %和0.90%。根據評價結果，湟水河流域以中高易發區為主，占研究區總面積的63.18%（圖9）。

表7 突變理論模型危險性分區統計表Table 7 Statistical table of hazard zoning using catastrophe theory model

圖9 易發性評價因子歸一化圖Fig.9 Normalized evaluation factors for susceptibility assessment

4.3 評價結果精度驗證

ROC 曲線為二分類模型的預測精度提供了定量化的研究方法[14-15]。評價過程簡單，評價結果非常直觀，以至于能僅憑肉眼進行預測效果初步判斷（圖10）。ROC 以曲線下面積判斷預測精度好壞（AUC）。對于任何預測試驗中，AUC值在0.5～1。特征曲線的基本原理是在崩滑危險性區劃結果中隨機選取了等量的崩滑點和隨機點，來繪制ROC 曲線，并確定曲線下面積AUC值[16-17]。其評判標準如表8 所示。

表8 AUC 值評判依據Table 8 Criteria for AUC value evaluation

圖10 基于突變理論的湟水河流域易發性模型ROC 精度驗證結果Fig.10 ROC accuracy validation results for the susceptibility model in the Huangshui River Basin based on catastrophe theory

湟水河流域崩滑災害分布特征與基于突變級數理論的易發性評價結果基本一致。為了驗證這一結果，在ArcGIS 隨機抽取30 %的災害點并創建同等數量的隨機點，通過SPSS 軟件中的分析工具ROC 曲線，獲得基于突變級數理論的湟水河流域易發性評價ROC 曲線。根據結果，基于突變級數理論的湟水河流域易發性模型AUC值為0.843，根據評判標準，表明該預測模型精度很好，這與流域內崩滑災害災情狀況基本相一致，說明該模型滿足評價要求。

5 結論

（1）本文以地質因素和環境因素兩個變量為準則層，地貌、地形、工程地質巖組、地質構造、水文氣象、人類工程活動 6 個變量作為中間層，高程、坡度、平面曲率、剖面曲率、工程地質巖組、地形濕度指數（TWI）、距離斷層、年均降水量、地貌類型、距離水系、距離公路、距離城鎮 12 個控災因子作為最底層變量，建立了基于突變理論湟水河流域崩滑災害危險度評價逐層結構模型。

（2）采用ArcGIS 重分類自然間斷點法將研究區內1 171 處崩滑災害點危險度計算結果劃分為輕 度危險、中度危險、重度危險和極度危險4 個等級，結果顯示，研究區內1 171 處崩滑災害中輕度危險點75 處、中度危險點201 處、重度危險點476 處、極度危險點419處。研究區內災害點以重度和極度危險為主，與實際災情狀況一致。

（3）開發了一種新的崩滑災害易發性評價平臺——突變級數理論平臺對湟水河流域崩滑災害易發性進行評價。結果表明，湟水河流域崩滑災害的極高易發區、高易發區主要分布于湟水河南北兩岸以及一級支流兩側的低山丘陵區，面積分別為4 661.79 km2和5 511.44 km2，占整個研究區總面積的 28.95% 和 34.23%；中易發區分布于湟水河南北兩岸以及各大支溝兩側的構造侵蝕中山區，面積占研究區總面積的 23.99 %；低易發區、極低易發區主要分布在湟水河河谷平原地帶及全流域中高山區，面積占比分別占總面積的 11.80 % 和 0.90%。

（4）通過ROC 曲線驗證，基于突變級數理論的湟水河流域易發性模型ACU值為 0.843，表明該預測模型精度很好，其優點是操作簡單，評價結果可靠，且彌補了突變級數只能進行單體災害危險性分析的短板，為崩滑災害防治提供理論依據。