劉 佳，趙海軍，馬鳳山，郭 捷，孫琪皓，段學良

(1.中國科學院地質(zhì)與地球物理研究所/中國科學院頁巖氣與地質(zhì)工程重點實驗室，北京 100029；2.中國科學院地球科學研究院，北京 100029；3.中國科學院大學，北京 100049)

0 引言

青藏公路是世界上海拔最高、線路最長的公路，自建成后便加快了內(nèi)陸與高原地區(qū)的物資運輸速度，促進了西藏經(jīng)濟建設和人民生活改善。公路沿線泥石流災害是損毀路基、路面及橋梁、隧道等構造物的重要外在動力，也是造成公路長時間中斷的主要原因之一[1]。所以，利用數(shù)學方法和計算機技術，建立更高效的評價模型進行危險性分區(qū)，對泥石流的預測和防治至關重要。

國內(nèi)外對泥石流危險性評價的研究自20世紀伊始。HOLLINGSWORTH等[2]采用因子疊加法和打分方法構建滑坡泥石流危險度評價的基本框架，實現(xiàn)了滑坡危險度的評估；MANTOVANI等[3]基于RS技術對歐洲滑坡災害進行研究和分區(qū)；CHUNG等[4]針對調(diào)查區(qū)的滑坡災害情況，采用GIS的專業(yè)模型擴展分析功能進行了多指標綜合區(qū)劃；YALCIN[5]利用層次分析法和二元統(tǒng)計法研究了滑坡災害易發(fā)性；PRADHAN等[6]利用神經(jīng)網(wǎng)絡方法結合邏輯回歸模型對馬來西亞檳島地區(qū)滑坡災害進行危險性分區(qū)；YANG等[7]對比了兩種集對分析方法在泥石流危險性評價中的應用；殷坤龍等[8]應用GIS調(diào)查研究后發(fā)現(xiàn)滑坡災害發(fā)育具有區(qū)域性和重復性的規(guī)律；黃潤秋等[9]利用GIS技術和數(shù)學模型方法進行地災區(qū)劃；LIU[10]提出了泥石流災害等級與發(fā)生頻率之間指數(shù)關系的理論解；張春山等[11]研究了區(qū)域地質(zhì)災害風險評價要素權值的計算方法；薛東劍等[12]借助遙感數(shù)據(jù)解譯和地面調(diào)查以及GIS技術，運用信息量模型和邏輯回歸模型完成了對青川、平武縣的地質(zhì)災害危險性評價。REN等[13]建立了SPA-GRA模型并利用了變異系數(shù)法來確定礦山參數(shù)；由此看來，隨著科技進步、社會發(fā)展，GIS等計算機科學技術和多種數(shù)學方法逐漸被應用于泥石流危險性評價。

盡管評價方法不斷增加，但目前泥石流評價仍存在問題[14]：如評價指標賦權過程中，層次分析法主觀性過強，客觀變異系數(shù)法具有一定局限性、未考慮主觀因素的影響等等。本文以國道G109沿線拉薩至那曲地區(qū)路段為例，通過利用層次分析法來改進變異系數(shù)法，彌補了原方法適用性差的缺陷，并利用改進的變異系數(shù)法對研究區(qū)進行了泥石流危險性分區(qū)，從而確定危險區(qū)域。

1 評價模型構建

1.1 評價指標的選取

評價指標是整個評價體系最基本的組成元素。泥石流災害的形成與其所處的地貌形態(tài)、物質(zhì)、構造、水文、氣象和植被條件等多種因素密切相關[15]，在充分考慮調(diào)查區(qū)災害特點的基礎上，并著重考慮泥石流災害形成所需要的三個基本條件[16]，選取地層巖性、地形地貌、地質(zhì)構造(活動斷裂)、氣候氣象與地震五類因子作為公路沿線泥石流危險性評價的指標因子，以此作為評價體系中的準則層，選取依據(jù)如下：

(1)區(qū)內(nèi)主要分布有花崗巖、砂巖、頁巖及第四系堆積物等，易風化的巖漿巖和較軟弱巖層在溫度變化劇烈、日照強烈的青藏高原，更易于風化、崩解堆積于泥石流溝兩岸，往往能為泥石流提供更加豐富的物質(zhì)基礎。

(2)地形條件的影響主要體現(xiàn)在重力勢能以及能量轉化梯度方面，而重力勢能的高低和能量轉化率實質(zhì)上受坡體地形坡度和起伏度(泥石流坡體物源區(qū)至溝口的垂直距離)的影響：坡度不僅影響斜坡內(nèi)的應力分布，同時對斜坡上松散物質(zhì)的堆積厚度起著決定性作用，進而決定著斜坡的穩(wěn)定性[17]。相對高差大和地表切割強烈地區(qū)的巖土體剝蝕作用強烈，水流速度相對較快，具備了泥石流的賦存條件[18]。

(3)由于歐亞板塊與印度板塊的不斷擠壓，青藏高原構造活動劇烈。斷裂構造越發(fā)育則巖體中次生結構面越多、巖體越破碎。構造作用主要體現(xiàn)在斷層等宏觀表現(xiàn)上，主要根據(jù)采樣區(qū)內(nèi)各類構造線的密度分布狀況[19]，對所發(fā)育的泥石流溝的固體物儲備方量、補給形式以及泥石流的性質(zhì)有很大影響。地質(zhì)構造的復雜程度將直接影響到巖層的穩(wěn)定性[20]；進而為泥石流發(fā)育提供豐富的物源基礎。

(4)氣候也是泥石流產(chǎn)生的重要影響因素，它不僅影響風化作用，還會因為暴雨、連續(xù)降雨誘發(fā)泥石流滑坡等地質(zhì)災害。降雨過程中，孔隙水壓力不斷升高，直接影響滑裂面以下土體向泥石流轉化的速度[21]。

(5)地震烈度越高(地震加速度越大)則巖土體越易遭受破壞，為泥石流發(fā)育提供豐富物質(zhì)條件；地震通過為泥石流提供固體物質(zhì)來源、水源、動力條件和激發(fā)條件等影響泥石流的形成和發(fā)展[22]。

1.2 評價指標體系的建立

根據(jù)層次分析法的相關理論，評價指標體系一般包括目標層、準則層和指標層。目標層是指泥石流危險性的評價；準則層即評價危險性的指標因子，選取原則如上所述；指標層是對準則層的細化，綜合考慮調(diào)查區(qū)內(nèi)外部因素，將準則層細分為六個指標：巖石易風化程度、平均坡度、起伏度、斷層線密度(單位面積內(nèi)斷裂線的長度)、地震加速度、汛期(6～9月)月均降水量。對三個評價分析層次分別用字母A、B、C表示，最終建立危險性評價指標體系(圖1)。

1.3 改進變異系數(shù)法的提出

泥石流危險性評價因子權重的確定有主觀和客觀兩種方法。對于主觀方法來說：傳統(tǒng)層次分析法在計算權重系數(shù)時，相對權的最小平方法復雜一些，但其所需的一致性檢驗步驟保證了判斷矩陣的合理性和一致性。使用層次分析法主要是為了定性定量地研究所賦權重問題，具有試探性，所構造判斷矩陣具有很大主觀性，若是專家經(jīng)驗不足或者對所構造矩陣意見不一就會導致準確度較低甚至錯誤的結論。對于客觀變異系數(shù)法來說：此方法具有一定局限性，僅考慮數(shù)據(jù)的變異離散程度而確定權重系數(shù)，這點對于不可量化的評價指標因子不適用。所以本文對變異系數(shù)法靈活創(chuàng)新運用，利用層次分析法來修正可量化評價的指標因子所占權重，從而得到綜合考慮主觀、客觀兩方面因素的評價模型。

1.3.1層次分析法

本文將采用兩種層次分析法，即傳統(tǒng)層次分析法和權的最小平方法，要使上述兩種方法所得的權重值之間的差異程度與其相對應的分配系數(shù)間的差異程度相一致,引入距離函數(shù)的概念[23]，得到相應比例系數(shù)從而取得最終層次分析法所得權重值。

(1)傳統(tǒng)層次分析法

此方法大致可以分為四個步驟(詳細計算過程及一致性檢驗方法參見文獻[24])：在這一過程中根據(jù)專業(yè)人員建議，考慮西藏地區(qū)特殊的地質(zhì)環(huán)境條件，構造判斷矩陣見表1和表2。

表1 指標判斷矩陣(1)

表2 指標判斷矩陣(2)

(2)權的最小平方法

權的最小平方法是由A.T.W.Chu等于1979年提出，于1988年介紹到我國[25]。權的最小平方法的基本原理是在一定的約束條件下，建立目標函數(shù)，對其一致性判斷矩陣A求得最小值，即一致性判斷矩陣排序向量的最優(yōu)解[26]。由權的最小平方法構成的最優(yōu)化問題，用向量形式可表示為[27]：

(1)

②約束條件：

(2)

通過求解可得最優(yōu)解為：

(3)

式中：W*——一致性判斷矩陣排序權向量。

(4)

e=(1,1,1,…,1)T

(5)

式中：aij——判斷矩陣中的數(shù)值；

wi或wj——權重系數(shù)。

1.3.2對變異系數(shù)法進行改進

變異系數(shù)又稱為標準差率，它是利用求取一組數(shù)據(jù)的標準差與平均值之比，來作為評價體系指標權重的客觀賦權方法。該方法可以直接利用ArcGIS軟件中關于指標因子的柵格數(shù)據(jù)文件的平均值和標準差，通過簡單計算就可以得到各權重系數(shù)，具體計算公式如下：

(6)

(7)

式中：S——某組指標因子數(shù)據(jù)的標準差；

R——平均值；

CVi——某組指標因子的變異系數(shù)；

Mi——指標因子所占的客觀權重。

此處，首先選取傳統(tǒng)層次分析法和權的最小平方法進行組合賦權，各自的權重系數(shù)依據(jù)距離函數(shù)所得，目的是減少誤差、增大可信度。

W12=αW1+βW2

(8)

α=0.511，β=0.489

(9)

式中：W12——聚合賦權法所確定的權重；

W1——傳統(tǒng)層次分析法所確定的權重；

W2——權的最小平方法所得結果。

其次，利用組合賦權后的權重系數(shù)，對變異系數(shù)法進行改進。對于不可量化的評價指標，權重系數(shù)取層次分析法所得結果，這一過程可以避免變異系數(shù)法不能衡量不可量化指標權重的缺點；對于可量化評價指標，本著總權重和為1的原則，權重系數(shù)為由層次分析法所得的總權重系數(shù)與變異系數(shù)法所得的各個評價指標權重系數(shù)的乘積：

(10)

式中：Vi——評價指標改良后最終權重系數(shù)；

W12i——利用上述組合賦權法確定的權重系數(shù)；

Mi——變異系數(shù)法確定的指標權重。

1.4 綜合評價模型

確定各個評價指標因子的權重，僅僅是完成了該指標因子在致災評價中作用高低的評價，而要想獲得最終的綜合評價結果以及建立綜合風險性分區(qū)，需要利用各個指標因子的權重與分級分數(shù)的乘積確定。本文采用加權綜合評分模型進行綜合評價，公式如下：

(11)

式中：Q——泥石流危險性系數(shù)；

xi——第i個指標因子在某一評價單元的分級賦值分數(shù)；

vi——該指標因子的權重系數(shù)。

2 研究區(qū)域危險性分析

2.1 研究區(qū)概況

研究區(qū)域位于西藏自治區(qū)拉薩市至那曲市的國道G109(毗鄰青藏鐵路)附近10 km范圍內(nèi)，地處青藏高原東南部，北鄰唐古拉山脈，南靠喜馬拉雅山脈，地勢總體北高南低，西高東低；構造活動劇烈，多高山河谷地貌；區(qū)域內(nèi)包括西藏拉薩市和那曲市的眾多重要地區(qū)：拉薩、當雄、那曲、安多和唐古拉山口(其為西藏自治區(qū)和青海省的分界)。氣候?qū)儆诟咴瓪夂?，平均氣溫低，晝夜溫差大。?jù)已有調(diào)查資料顯示：那曲地區(qū)已發(fā)生地質(zhì)災害110多處，造成幾千億元經(jīng)濟損失。

2.2 計算過程

2.2.1數(shù)據(jù)來源

本文研究數(shù)據(jù)主要涉及泥石流評價指標因子的數(shù)據(jù)。其中，地層巖性和斷層線密度數(shù)據(jù)主要來自于1∶50 萬地質(zhì)圖資料；地形坡度和起伏度的DEM數(shù)據(jù)來源于中國科學院計算機網(wǎng)絡信息中心地理空間數(shù)據(jù)云平臺(http://www.gscloud.cn)，借鑒前人泥石流危險性評價經(jīng)驗[28-29]采用30 m精度(考慮到數(shù)據(jù)獲取的難易程度和最大利用效果)；汛期月均降雨量數(shù)據(jù)來源于全國降雨統(tǒng)計資料；地震加速度數(shù)據(jù)來自于《中國地震動峰值加速度區(qū)劃圖》(GB 18306—2015)。將這些數(shù)據(jù)導入ArcGIS軟件并劃分網(wǎng)格，共有11 917 個面積為1 km×1 km的網(wǎng)格，形成柵格數(shù)據(jù)文件，并將其分為五級。

2.2.2權重系數(shù)的計算

根據(jù)前面章節(jié)內(nèi)容提到的評價指標權重的確定方法，代入數(shù)據(jù)整理計算后可得結果。傳統(tǒng)層次分析法(W1)和權的最小平方法(W2)、組合后的層次分析法(W12)、變異系數(shù)法(M)以及改進變異系數(shù)法(V)的計算結果見表3。

表3 各種方法所得權重

圖2 評價指標分級圖Fig.2 Classifying picture of evaluation index

2.2.3確定評價指標分級分數(shù)

泥石流的形成離不開特定的自然地理、地質(zhì)和氣候條件：利用ArcGIS軟件，根據(jù)上述所選評價指標因子，對西藏地區(qū)拉薩至那曲國道G109附近10 km范圍(前人研究主張距離鐵路5 km以上的泥石流均不會對青藏鐵路交通干線產(chǎn)生破壞性影響[30]，為了綜合評價更廣泛的區(qū)域災害特點，本文決定擴大為10 km范圍)內(nèi)的區(qū)域進行指標因子等級劃分(表4)，結合“《中國公路綜合自然區(qū)劃》體系框架研究”和文獻[31]的研究成果，將六類評價指標分別劃分為五個等級(圖2)；為便于更直觀地了解等級劃分結果，對其中部分指標數(shù)據(jù)采取適當取整處理，分別賦值1、3、5、7、9；與大部分研究中都選擇15作為分級數(shù)值相比，此方法使得分級分數(shù)間拉開差距，更具明顯性(g為地震加速度，取值10 m/s2)。

2.3 結果驗證分析

利用ArcGIS數(shù)據(jù)處理功能，得到研究區(qū)域泥石流災害各個指標因子的分級圖，在此基礎上，利用軟件的加權疊加功能，形成泥石流綜合評價危險性分區(qū)圖(圖3)。借助其數(shù)據(jù)統(tǒng)計功能，對改進變異系數(shù)法所建模型進行分析：泥石流危險度中高區(qū)域(中等、較高、高)面積約2.3×103km2，占調(diào)查區(qū)總面積的20%，主要分布在：國道G109拉薩至羊八井路段約1 000 km2范圍內(nèi)；羊八井—當雄路段道路右側部分區(qū)域；當雄縣東北方向約30 km處的道路左側拐角處約170 km2范圍；以及安多縣周圍部分區(qū)域。泥石流危險度低區(qū)域(較低、低)約占總調(diào)查面積的80%：大部分分布在那曲縣城北部。將實際災害數(shù)量分布特點(圖3)和危險性分區(qū)圖對照，結合災害點數(shù)量在各個等級危險區(qū)的分布情況(圖4)，綜合來看評價可靠度較高。根據(jù)現(xiàn)場調(diào)查結果及四川省廣漢地質(zhì)工程勘察院所撰寫的地災調(diào)查與區(qū)劃報告資料得出，23處(約占泥石流總量的79%)泥石流位于評價模型中高危險區(qū)，其余6處(約占泥石流總量的21%)處于低、較低危險區(qū)。根據(jù)現(xiàn)場調(diào)查情況分析，圖3中高危險區(qū)往往處于高山峽谷段，地勢較陡，而低危險區(qū)往往地勢平坦，這就驗證了改進變異系數(shù)方法具有較好的準確度和合理性。

表4 泥石流評價指標等級劃分

圖3 危險性分布圖Fig.3 The distribution map of hazards

圖4 災害數(shù)量分布圖Fig.4 The quantity distribution map of disaster

3 結論

本文針對西藏拉薩至那曲地區(qū)的國道G109附近10 km范圍內(nèi)的泥石流危險性評價問題，利用改進變異系數(shù)法借助ArcGIS的數(shù)據(jù)統(tǒng)計和加權疊加等功能，構建了高效合理的泥石流危險性評價模型，并進行泥石流危險性分區(qū)及評價，所得主要結論包括：

(1)提出了一種改進變異系數(shù)法，該方法克服了層次分析法的主觀性及客觀變異系數(shù)法的局限性，并將其應用于泥石流危險性綜合評價模型，以此確定的中高危險區(qū)與泥石流高發(fā)區(qū)重合度較大，此方法可靠性較高。

(2)研究區(qū)泥石流中高危險區(qū)域面積約2.3×103km2(僅占調(diào)查區(qū)總面積的20%)，泥石流發(fā)生數(shù)量占總量的79%，國道G109東側比西側危險性更高，拉薩至當雄路段尤其明顯；泥石流低危險區(qū)約9.6×103km2，占總調(diào)查面積的80%，但災害數(shù)量僅占21%。

(3)變異系數(shù)法在地質(zhì)災害危險度研究中還不是很廣泛，利用改進變異系數(shù)法，可克服原方法只能計算可量化指標的缺點，具有更廣泛的應用空間。