劉福臻，王 靈，肖東升，王軍朝

(西南石油大學 土木工程與測繪學院，四川 成都 610500)

滑坡易發(fā)性評價作為危險性評價和風險性評價的基礎，是防災減災過程中不可缺少的程序。近年來，隨著計算機技術和地理信息技術的快速發(fā)展，滑坡的易發(fā)性評價逐漸成為了研究熱點[1]。但由于滑坡發(fā)育的復雜性，目前還沒有一套完全成熟的理論能對其進行有效的防治。為此，國內(nèi)外學者對滑坡的易發(fā)性評價開展了一定程度的研究。李環(huán)禹等[2]和Kayastha等[3]從因子角度出發(fā)，分別探討了降雨因子和多因子之間的組合對滑坡易發(fā)性評價的影響。田述軍等[4]和Xiaohui Sun等[5]從評價單元角度出發(fā)，論述了斜坡單元和網(wǎng)格單元對評價精度的影響。李文彥等[6]和Bourenane H等[7]從多個模型對比的角度出發(fā)，探討了信息量、頻率比、邏輯回歸、層次分析法等模型在滑坡易發(fā)性評價中的精度表現(xiàn)。Nsengiyumva J B等[8]和孫德亮等[9]以機器學習為基礎，探討了黑箱模型在易發(fā)性評價中的應用。以上研究分別從不同角度對滑坡的易發(fā)性進行了探索，一定程度上促進了滑坡災害防治學科的發(fā)展，但鮮有學者對滑坡評價因子的量化和模糊綜合評價模型的隸屬函數(shù)作相應的研究。

鑒于此，本文以統(tǒng)計學理論為基礎，將歷史滑坡點與評價因子的二級屬性進行疊加分析，利用相對滑坡點密度作為本次實驗的量化值。該方法從數(shù)據(jù)本身的結構信息出發(fā)，實現(xiàn)了滑坡因子二級屬性的客觀量化。模型上，采用模糊綜合評判法作為評價模型，該模型以隸屬函數(shù)為基礎計算出每個評價單元在各分區(qū)的概率，通過選擇其中最大值所對應的易發(fā)性分區(qū)作為評價結果，從而很好的解決了評價等級的劃分問題[10]。模型的缺陷在于模糊隸屬函數(shù)的建立存在較強的主觀性，且現(xiàn)有的隸屬函數(shù)都不太適合滑坡的易發(fā)性評價。基于此，本文將反距離權重的思想用于確定因子屬性的隸屬度，探討該想法在評價過程中的可行性與有效性。

1 研究區(qū)概況及數(shù)據(jù)源

1.1 研究區(qū)概況

寧南縣位于我國四川省涼山州東南部，是連接川、滇兩地的重要通道。研究區(qū)內(nèi)崇山峻嶺，溝壑縱橫，黑水河貫穿其南北，于寧南縣葫蘆口處匯入金沙江。研究區(qū)氣候類型為南亞熱帶季風氣候，具有“一山分四季，十里不同天”的立體特征[11]。復雜多變的地貌和氣候特征為該地區(qū)提供了豐富的礦產(chǎn)資源，同時也導致了該地區(qū)地質(zhì)災害的頻發(fā)。

據(jù)寧南縣2019年1：5萬地質(zhì)災害詳查顯示，全縣共有482個地質(zhì)災害點。其中滑坡359個，泥石流100個，不穩(wěn)定斜坡12個，崩塌11個?；抡紦?jù)了所有災害點的74.5%，是寧南縣最主要的地質(zhì)災害。

圖1 研究區(qū)地理位置Fig.1 Geographical location of the study area

1.2 數(shù)據(jù)源

本文用到的數(shù)據(jù)主要有滑坡歷史災害點、寧南縣DEM、寧南縣土地利用覆蓋、1∶25萬地質(zhì)圖，詳細內(nèi)容見表1。其中滑坡災害點主要用于和各因子的疊加分析；DEM主要用于高程、坡度、坡向、坡型和水系的提??；1∶25萬地質(zhì)圖主要用于地層和斷層的提取。將所有圖層的投影坐標系設置為WGS_1984_UTM_Zone_47N，并重采樣為30*30的柵格單元。

表1 數(shù)據(jù)源Table 1 Data source

2 模型介紹

2.1 因子正向化(量化)

因子正向化是指將量化的因子轉換成與滑坡發(fā)生成正相關形式的過程。由于數(shù)據(jù)的多源性，原數(shù)據(jù)通常不能直接為我們所用。以連續(xù)因子中的坡度為例，寧南縣的坡度范圍大約為0～75°，然而有利于滑坡發(fā)育的范圍并不是隨著坡度的增加而變大。同理，對于類似坡向、土地利用、巖組等分類數(shù)據(jù)，都需要將之量化為能表現(xiàn)滑坡易發(fā)性指數(shù)的形式。因此，本文從統(tǒng)計學原理出發(fā)，利用歷史滑坡點數(shù)據(jù)，選取相對滑坡點密度作為各因子的正向化值。相對滑坡點密度計算方法見式(1)。

(1)

式中Xi為因子X在二級屬性i下的相對滑坡點密度；ni為因子X在二級屬性i下的滑坡個數(shù)；si為研究區(qū)內(nèi)因子X在二級屬性i下的柵格個數(shù)；N為研究區(qū)總滑坡個數(shù)；S為研究區(qū)總柵格數(shù)。

2.2 熵權法

Shannon在1948年提出了熵值理論[12]，熵值是對不確定信息的一種度量，對系統(tǒng)或事件離散性和無序性的概述。其主要思想為：指標的變異程度越大，對應的權值也越大。在本文中體現(xiàn)為滑坡因子二級屬性值的變異程度。也就是說，對于所有的樣本而言，倘若某個因子正向化后的二級屬性值都相同，則認為這個指標的權值為0，即該因子對評價系統(tǒng)起不到任何幫助。同理，某個正向化后的因子的屬性值變異越大，則該因子對滑坡的影響也就越大[13]。具體計算過程如下：

(1)假設x表示事件X可能發(fā)生的某種情況，P(x)表示這種情況發(fā)生的概率。

我們可以定義：

I(x)=-ln(P(x)).

(2)

(2)如果事件X可能發(fā)生的情況分別為：x1,x2,x3…xm，事件X的信息熵值H(X)可用式(3)計算。

(3)

(3)最后，指標的熵權計算方法見式(4)。

(4)

式中wi為第i個因子的權重；Hi為第i個因子的信息熵；n為評價因子的類別個數(shù)。

2.3 模糊綜合評判法模型

2.3.1 模糊綜合評判法模型介紹

模糊綜合評判模型通過對因子的屬性做隸屬度賦值，以此來達到評判的目的。其核心在于用隸屬程度代替屬于和不屬于，即你有黑的可能性，也有白的可能性[14]。評價過程以隸屬函數(shù)和指標權重為基準，通過唯一確定評價單元的評判等級，從而解決了指標等級劃分問題。具體步驟如下：

(1)設因素集U={u1,u2,u3,…,um}，評價集V={v1,v2,v3, …,vn}。因素論域和評價論域之間的模糊關系用評價矩陣R來表示：

(5)

其中rij=μ(ui,vj)表示從因素ui的角度來判定事件評為vj的隸屬度。

(2)給定出影響因子的權重。

W=(w1,w2,w3,…,wi)，其中wi為第i個因子的權重，該權重通過熵權法計算獲得。

(3)利用模糊矩陣R和權重矩陣W計算評價集Z。

Z=(z1,z2,z3,…,zn)=W1×m⊙Rm×n.

(6)

(4)最后根據(jù)每個評價單元的最大隸屬度b來確定評價單元的等級vi。

b=max(z1,z2,z3,…,zn).

(7)

2.3.2 隸屬函數(shù)的確定

隸屬函數(shù)的確立，是模糊數(shù)學及其應用在理論上和實踐中的關鍵問題[15]。常用的隸屬函數(shù)有模糊統(tǒng)計法和指派法[16]。然而這些常規(guī)的方法都不太適合滑坡的易發(fā)性評價?；诨掠绊懸蜃拥膹碗s性，本文將反距離權重法的思想用于模糊隸屬度的計算。過程如下：

(1)利用相對滑坡點密度給每個影響因子的二級屬性賦值，并選出與評價等級相匹配的屬性值作為相應等級的基準。

(2)對于評價集V={v1,v2,…,v5}，因子ui對應的隸屬度rij=μ(ui,vj)的計算方法見式(8)。

ri1=1-|xi-xi1|
ri2=1-|xi-xi2|
ri3=1-|xi-xi3|
ri4=1-|xi-xi4|
ri5=1-|xi-xi5|.

(8)

式中xi是評價因子ui的量化值，xi1,xi2…xi5是評價集v1,v2,…,v5的基準。xi1,xi2…xi5的選擇一般由ui量化后，根據(jù)從小到大的順序而定。

3 滑坡分布特征及評價因子量化結果

3.1 坡度

由1∶5萬地質(zhì)災害詳查結果得知，全縣76.3%的滑坡在5°到25°的平緩地帶。其中5°—15°為量化的峰值。量化結果整體呈現(xiàn)的趨勢為：隨坡度的增加先變大后減小。研究表明，雨水在緩坡地段容易匯集，使得雨水的滲透作用加強，導致土石層飽和，從而增加了滑體的重量，降低了土石層的抗剪能力[17]。所以寧南縣的滑坡主要聚集在平緩地帶。坡度較陡的地帶大部分為巖質(zhì)斜坡，通常情況下巖質(zhì)斜坡在抗剪力上遠遠超過土質(zhì)斜坡，這些斜坡往往發(fā)育為崩塌，所以在坡度大于45°的研究區(qū)內(nèi)沒有滑坡災害點的分布。本文將研究區(qū)坡度分為了[0，5°]、[ 5°，15°]、[ 15°，25°]、[ 25°，35°]、[ 35°，45°]、[ ≥45°]六大類，圖2為滑坡坡度分布統(tǒng)計及其量化圖。

圖2 坡度統(tǒng)計量化結果Fig.2 Statistics and quantitative results of slope

3.2 坡向

利用ARCGIS中的坡向工具將DEM數(shù)據(jù)轉化為坡向數(shù)據(jù)。寧南縣的坡向分為平地、北、東北、東、東南、南、西南、西和西北9大類型。經(jīng)實地調(diào)查發(fā)現(xiàn)，在寧南縣范圍內(nèi)，除了平地沒有滑坡外，各坡向量化的相對滑坡點密度值起伏變化不大。整體的數(shù)據(jù)走勢表明坡向因子在研究區(qū)范圍內(nèi)對滑坡的影響較小。圖3為坡向分布統(tǒng)計及其量化圖。

圖3 坡向統(tǒng)計量化結果Fig.3 Statistics and quantitative results of aspect

3.3 坡型曲率

通常情況下，山體的坡型特征可以這樣描述：曲率大于0的地方為凸型地貌，曲率小于0的地方為凹形地貌。然而這樣簡單的分類方法并不能很好的展現(xiàn)斜坡復雜的形態(tài)。為了更為細致的突出滑坡分布特征，本文沒有采用傳統(tǒng)的凹凸型分類，而是根據(jù)坡型曲率將研究區(qū)分為了六類：[≤-1.3]、[-1.3，-0.56]、[-0.56，0]、[0，0.45]、[0.45，1.2]、[≥1.2]。圖4為坡型曲率分布統(tǒng)計及其量化圖。

圖4 坡型統(tǒng)計量化結果Fig.4 Statistics and quantitative results of slope type

3.4 水系

寧南縣溝谷地貌發(fā)育，受河流沖刷影響，河流兩岸斜坡風化侵蝕嚴重。由統(tǒng)計數(shù)據(jù)可知，在研究區(qū)域內(nèi)，河流緩沖距離在300～600 m之間的量化值最大。量化的整體趨勢可描述為：離河道越近，發(fā)生滑坡的數(shù)量越多，相對滑坡點密度也越大。圖5為水系分布統(tǒng)計及其量化圖。

圖5 水系統(tǒng)計量化結果Fig.5 Statistics and quantitative results of water system

3.5 巖組

根據(jù)寧南縣地層分布情況，將研究區(qū)內(nèi)分布的地層巖組劃分為9大類，分別為：Z、Pz1、P-T、Qb-Z、T3C、J、O-S、Qal、-O。45%的滑坡點分布于P-T和J這兩組地層，其中P-T類型的量化值最大，高達3.82。圖6為巖組分布統(tǒng)計及其量化圖。

圖6 巖組統(tǒng)計量化結果Fig.6 Statistics and quantitative results of rock group

3.6 土地利用覆蓋

根據(jù)國際分類標準，寧南縣范圍內(nèi)的土地利用覆蓋可分為耕地、森林、草原、灌木叢、水、防滲透表面8個類型。其中有182個滑坡發(fā)生在耕地范圍內(nèi)，占據(jù)了總滑坡的51%左右，而在灌木叢和防滲水面一共只有5個滑坡災害點。由此可見，寧南縣范圍內(nèi)土地利用覆蓋因子在滑坡的發(fā)育過程中起著重大作用。圖7為土地利用覆蓋分布統(tǒng)計及其量化圖。

圖7 土地利用覆蓋統(tǒng)計量化結果Fig.7 Statistics and quantitative results of land use cover

3.7 斷層

斷層是構造運動中最為廣泛的形態(tài)之一。地區(qū)的地形地貌以及巖層的連續(xù)性通常和斷層的發(fā)育有著極強的聯(lián)系[18]。根據(jù)寧南縣1∶25萬地質(zhì)圖和野外調(diào)查結果顯示，滑坡發(fā)育的地帶通常距斷層較近，而距離斷層較遠的區(qū)域很少發(fā)生滑坡。本文通過統(tǒng)計分析，將斷層因子分為了七大類，分別為:[≤500 m]、[ 500 m，1 000 m]、[ 1 000 m，1 500 m]、[ 1 500 m，2 000 m]、[ 2 000 m，2 500 m]、[ 2 500 m-3 000 m]、[ ≥3 000 m]。圖8為斷層分布統(tǒng)計及其量化圖。

圖8 斷層統(tǒng)計量化結果Fig.8 Statistics and quantitative results of fault

3.8 高程

寧南縣地貌復雜，地形起伏大，最低海拔526 m，最高海拔3 869 m。統(tǒng)計顯示，有119個滑坡位于海拔1 094 m到1 466 m范圍之間，在海拔大于2 570 m后基本沒有滑坡歷史點?；诖耍米匀粩帱c法將高程分為了7大類，分別為：[526 m，1 094 m]、[1 094 m，1 466 m]、[1 466 m，1 824 m]、[1 824 m，2 190 m]、[2 190 m，2 572 m]、[2 572 m，2 987 m]、[2 987 m，3 869 m]。圖9為高程分布統(tǒng)計及其量化圖。

圖9 高程統(tǒng)計量化結果Fig.9 Statistics and quantitative results of elevation

4 研究區(qū)易發(fā)性評價

4.1 因子相關性和權值計算

4.1.1 因子相關性計算

考慮到因子之間可能存在的相關性會使結果產(chǎn)生偏差，實驗對量化的因子進行了相關性分析，結果表明8個因子之間的相關性都小于0.33，因子之間沒有較強的相關性。計算結果中斷層因子和坡型因子的相關性為0，但卻沒有通過顯著性檢驗，其余因子均通過了5%的顯著性檢驗。實驗目的是想證明因子之間沒有相關性，因此，可以利用所選的8個因子進行滑坡易發(fā)性評價。詳細內(nèi)容見表2。

表2 因子相關性分析Table 2 Factor correlation analysis

4.1.2 熵權法計算權重

在計算因子權重之前，先將各因子的二級屬性標準化，標準化計算方法見式(9)。

(9)

將標準化后的相對滑坡點密度看作概率P(x)，利用式(3)計算每個因子的信息熵。例如，坡度的信息熵值計算方法見式(10)。

(10)

最后計算得出坡度、坡向、坡型、水系、巖組、土地利用、斷層、高程的信息熵分別為0.814 3、0.929 8、0.869 2、0.949 3、0.861 0、0.788 3、0.956 0、0.770 8。

再利用計算的信息熵代入式(4)，計算得出每個因子的權值，結果見表3。

表3 因子統(tǒng)計計算結果Table 3 Calculation results of factor statistics

從計算結果中我們可以分析得出，在南江縣范圍內(nèi)，高程和土地利用類型對滑坡的影響程度最大，而斷層的影響相對較小。

4.2 模糊綜合評價

4.2.1 反距離權重法計算隸屬度

由于滑坡發(fā)育的復雜性，傳統(tǒng)的隸屬函數(shù)不能凸顯滑坡因子屬性在各易發(fā)性級別的隸屬度，基于此，本文將反距離權重法的思想用于模糊隸屬度的計算。

以坡度因子為例，具體步驟如下：

(1)利用量化模型給坡度二級屬性賦值，其中[0，5°]、[ 5°，15°]、[ 15°，25°]、[ 25°，35°]、[ 35°，45°]、[ ≥45°]的屬性值分別為0.981 930 058、1.615 744 912、1.396 320 587、0.638 616 124、0.232 246 663、0。

(2)從量化模型的函數(shù)表達式得知，量化的值越大，發(fā)生滑坡的概率也越大。量化結果顯示坡度在[ 5°，15°]這一范圍內(nèi)的量化值最大，從而推斷出[ 5°，15°]這一范圍屬于高易發(fā)區(qū)的隸屬度也最大。在坡度[ ≥45°]這一范圍內(nèi)量化的屬性值為0，相應的這一范圍內(nèi)發(fā)生滑坡的概率小，屬于不易發(fā)區(qū)的概率也就大。本文將易發(fā)區(qū)分為了五大類，分別為不易發(fā)、低易發(fā)、中易發(fā)、高易發(fā)、和極高易發(fā)。由于屬性值1.616是坡度中最大的量化值，則選取該值作為極高易發(fā)區(qū)的一個基準，量化值離1.616越近，隸屬高易發(fā)區(qū)的值也越大，離1.615 7越遠，隸屬高易發(fā)區(qū)的值也就越小。計算方法見式(11)。

ri5=1-|xi-1.616|.

(11)

其中xi為滑坡二級因子的量化值，ri5為xi所對應的二級因子隸屬極高易發(fā)區(qū)的隸屬度。

同理，針對另外4個易發(fā)分區(qū)，根據(jù)量化值大小選取相應的基準。坡度上，本文選擇了0.232、0.639、0.982、1.396、1.616分別作為不易發(fā)、低易發(fā)、中易發(fā)、高易發(fā)、和極高易發(fā)的基準。相應的隸屬結果如表4。

表4 坡度隸屬度計算結果Table 4 Calculation results of slope membership

從表中我們可以得出，量化值越大的屬性，對應的柵格隸屬于高易發(fā)的可能性越大；量化值偏小的屬性，對應的柵格隸屬于低易發(fā)的可能性越大。以上信息說明該方法符合滑坡易發(fā)的基本情況，且利用反距離權重法能夠很好的過度各屬性在各分區(qū)上的隸屬值，避免了人為干涉對數(shù)據(jù)造成的影響，同時也保留了數(shù)據(jù)本身的離散度等基本信息。實驗證明，該方法能夠更科學、更客觀的對研究區(qū)進行易發(fā)區(qū)分區(qū)制圖。值得注意的是，由于某些因子在量化時最大值和最小值差距太大(巖層因子)，導致在使用反距離權重法賦隸屬度時出現(xiàn)不合理的情況，基于此，實驗選擇將這些異變的隸屬度根據(jù)實際情況進行了微調(diào)，從而使計算的結果更加符合實際情況。表5為具體的計算結果。

表5 各因子隸屬度計算結果Table 5 Calculation results of membership degree of each factor

4.2.2 模糊評價分區(qū)結果

通過上述步驟獲取了權重矩陣W和模糊矩陣R，再由式(12)計算評價集結果。

Z=(z1,z2,z3,z4,z5)=W1×8⊙R8×5.

(12)

對于每一個柵格，通過式(12)計算后會獲取5個值，取其最大值所對應的評價分區(qū)作為最終結果，計算方法見式(13)。

b=max(z1,z2,z3,z4,z5).

(13)

利用arcgis平臺，本文將研究區(qū)的柵格分別賦值為1、2、3、4、5表示對應的不易發(fā)、低易發(fā)、中易發(fā)、高易發(fā)、和極高易發(fā)區(qū)域，結果如圖10。

圖10 寧南縣滑坡易發(fā)分區(qū)圖Fig.10 Zoning map of landslide susceptibility in Ningnan County

將所得分區(qū)圖與滑坡歷史點位圖層進行疊加分析，結果顯示落在不易發(fā)、低易發(fā)、中易發(fā)、高易發(fā)和極高易發(fā)的滑坡點個數(shù)分別為2、18、22、144、173。采用相對滑坡比來驗證模型可行性。計算方法見式(14)。

(14)

其中ni為易發(fā)性分區(qū)i下的滑坡個數(shù)；si為易發(fā)性分區(qū)i的柵格個數(shù)；N為研究區(qū)總滑坡個數(shù)；S為研究區(qū)總柵格數(shù)。結果如表6。

表6 分區(qū)統(tǒng)計結果Table 6 District statistics results

5 結論

(1) 將各因子的統(tǒng)計特征量化為相對滑坡點密度，以此為基礎，使用熵權法計算獲得了各評價因子的權值，結果顯示高程和土地利用是研究區(qū)主要的滑坡易發(fā)性影響因子。

(2) 將反距離權重法的思想用于模糊隸屬度的計算。實驗證明該方法能夠很好的過度因子二級屬性在各分區(qū)上的模糊隸屬值，同時能夠更科學，精確的對研究區(qū)進行易發(fā)性分區(qū)制圖。

(3) 寧南縣滑坡易發(fā)性分區(qū)圖中，從不易發(fā)區(qū)到極高易發(fā)區(qū)所占的比例分別為18.7%、25.8%、11.3%、31.4%、和12.7%。其中高易發(fā)區(qū)和極高易發(fā)區(qū)多位于寧南縣東南側和黑河水系兩側。將易發(fā)性分區(qū)圖與歷史災害點圖層進行疊加分析，結果顯示有88.3%的滑坡分布于極高易發(fā)區(qū)和高易發(fā)區(qū)，從低易發(fā)到極高易發(fā)的相對滑坡比分別為0.03、0.19、0.54、1.28和3.78。說明該模型計算結果與實際情況高度相似，反距離權重法的思想在模糊隸屬度的計算上具有較高的適用性。