基于GIS與Logistic模型的公路泥石流易發(fā)性分析

鄒 強，王 青，劉延國

（1.西南科技大學 環(huán)境與資源學院，四川 綿陽621010；2.中國科學院 水利部成都山地災害與環(huán)境研究所 中國科學院 山地災害與地表過程重點實驗室，四川 成都610041）

泥石流作為山區(qū)環(huán)境中一種常見的自然現(xiàn)象，由于其形成過程復雜，暴發(fā)突然，來勢兇猛，破壞力強，常給山區(qū)道路造成毀滅性的災害。隨著山區(qū)資源開發(fā)與路網(wǎng)建設力度的加大，泥石流災害對道路的危害日益突出。如何定量分析泥石流的發(fā)生可能性對道路的正常運營具有重大的現(xiàn)實意義。泥石流災害易發(fā)性分析是通過分析研究區(qū)泥石流致災因子的變化情況，估計其發(fā)生概率，并在此基礎上完成易發(fā)性評價與制圖，可為泥石流的預測與風險分析提供基礎信息。

自20世紀90年代，泥石流易發(fā)性評價成為泥石流災害研究的熱點，在全國范圍內(nèi)展開，涌現(xiàn)出許多突破性的研究成果，如全國泥石流災害分布與災害區(qū)劃圖［1］、西藏泥石流監(jiān)測方法［2］等。隨著泥石流易發(fā)性分析的理論與方法的不斷深入與完善，信息技術的發(fā)展及廣泛應用，特別是遙感（RS）與地理信息系統(tǒng)（GIS）技術為泥石流易發(fā)性研究提供了強有力的空間分析功能，使得空間數(shù)據(jù)集成分析更加快捷、精準，泥石流易發(fā)性評價也從過去單一、定性研究向綜合化、定量化發(fā)展［3－9］。目前，泥石流易發(fā)性評價工作的關鍵在于針對不同的評估目的、不同的區(qū)域尺度，對相關空間數(shù)據(jù)進行處理和分析，開發(fā)出可靠性強的分析評價模型。

本文擬探討基于GIS與Logistic模型的公路泥石流易發(fā)性評價的原理與方法，并選擇G318川藏公路段為研究區(qū)進行相關理論與技術的實證研究。

1 研究區(qū)概況

川藏公路東起四川省成都，西止西藏拉薩，是連接西藏與西南地區(qū)其他城市的主要交通命脈。受青藏高原強烈隆升的影響，川藏公路途徑地區(qū)，地勢西高東低，自東向西跨越不同的地貌單元，公路沿線地層較復雜，從上元古界震旦系到新生界第四系均有分布，全線地質(zhì)構造活動強烈，空間展布東西差異明顯。東部地區(qū)北東向構造帶主要以龍門山斷裂帶為主；中部地區(qū)近南北向構造帶以橫斷山褶皺帶為主，主要包括“三江斷裂”［10］——金沙江斷裂帶、瀾滄江大斷裂帶、怒江大斷裂帶，對公路沿線的地貌格局起著主要的控制作用；中西部地區(qū)為中段南北向構造與西段東西向構造的過渡區(qū)域，斷裂、褶皺及花崗巖體呈北北西向展布；西部地區(qū)的斷裂及褶皺呈東西向展布，為雅魯藏布江斷裂與三江斷裂之間的緩沖地帶。川藏公路沿線地區(qū)地處西藏東南部，地勢對氣候控制明顯，直接影響公路沿線的降水、氣溫、濕度等氣候因素。此外，受西南季風的控制，公路沿線地區(qū)降水較豐富。

川藏公路是中國泥石流災害最為嚴重的交通干線之一，僅帕隆藏布沿江線路長271km，災害點有399處，危害 長度達71.63km［10－11］。古 鄉(xiāng)溝、培 龍溝、米堆溝、瓦達溝等都是比較典型的泥石流溝，歷史上數(shù)次大規(guī)模泥石流［11－15］，給川藏公路造成嚴重的生命財產(chǎn)損失。G318川藏公路段泥石流分布狀況如圖1所示。隨著經(jīng)濟建設的發(fā)展，頻繁的泥石流災害，使公路線路路況日趨惡化，車毀人亡事件時有發(fā)生，阻車斷道十分嚴重。

圖1 G318川藏公路段泥石流分布

2 研究方法

泥石流易發(fā)性評價的目的在于分析泥石流在一定的致災因子作用下發(fā)生的可能性或概率。目前，泥石流易發(fā)性分析方法主要包括定性與定量兩類方法。定性方法是借助于泥石流的理論知識、經(jīng)驗認識對泥石流發(fā)展變化規(guī)律進行科學分析與判斷，一般用于區(qū)域性的泥石流研究，如專家知識法［16］、層次分析法（AHP）［17］、加權線性組合法（WLC）［18］。定量方法是根據(jù)數(shù)據(jù)資料，應用科學的方法建立反映地質(zhì)災害與其控制因素之間關系的數(shù)學模型，進而實現(xiàn)定量化評價的一類方法，如確定性方法［19］、人工智能法［20］和多元統(tǒng)計法［21］。泥石流易發(fā)程度的特點是因變量的取值范圍介于“1”（泥石流發(fā)生）和“0”（泥石流不發(fā)生）這兩個值之間，這要求建立的模型必須保證因變量的取值范圍為［0，1］。二項Logistic模型所需的假設簡單，不要求滿足誤差分布趨于正態(tài)分布的假設，也不要求自變量符合正態(tài)分布的條件，模型對識別變量的分布未作任何要求［22］，可以用來預測具有二項特點的因變量概率問題，因此，該模型符合泥石流易發(fā)性評價的建模要求。

基于GIS與Logistic回歸模型的泥石流易發(fā)性分析方法包括三方面的內(nèi)容：（1）關鍵影響因子確定；（2）因子指標量化；（3）Logistic回歸模型構建。

2.1 泥石流影響因子確定

泥石流的發(fā)生是一個復雜的系統(tǒng)過程，影響因素很多，包括地質(zhì)條件、地形地貌條件、氣候條件、水文條件、植被條件和人類活動等，各因素的影響作用也不盡相同。道路泥石流影響因子的選擇的主要依據(jù)研究區(qū)內(nèi)泥石流的形成條件、發(fā)育現(xiàn)狀和災害分布，從泥石流形成與演化的物質(zhì)條件、能量條件和環(huán)境條件3個方面出發(fā)，選取坡度、高程、坡面粗糙度，地層巖性、距斷裂距離，地震烈度，土地利用類型等7個相對穩(wěn)定的評價因子構建指標體系，進而研究泥石流易發(fā)性。

針對川藏公路沿線泥石流發(fā)生的實際情況，充分考慮資料的可獲取性與研究尺度，高程、坡度和坡面粗糙度可通過地形數(shù)據(jù)獲取，地層、巖性、斷裂等地質(zhì)信息可通過研究區(qū)地質(zhì)圖提取，土地利用類型可依據(jù)研究區(qū)土地利用數(shù)據(jù)，參考第二次全國土地調(diào)查技術規(guī)程（TDT 1014—2007）進行劃分。

2.2 因子量化方法

在泥石流易發(fā)性分析過程中，由于各個評價因子的量綱不同，各個評價因子作為回歸模型中的自變量首先需要進行標準化處理，以獲取同一尺度的量化指標值，便于進行后續(xù)分析與計算。本文用每個影響因子的各個Ⅱ級分類因子中泥石流災害點的個數(shù)與各個Ⅱ級分類因子的面積的比值作為各個Ⅱ級分類因子的指標值。計算表達式如下。

式中：xi，j——各 類 Ⅱ 級 因 子 的 指 標 值 （條／km2）；Si，j——各類Ⅱ級因子的面積（km2）；Ni，j——各類Ⅱ級因子中泥石流數(shù)（條），i＝1，2，…，n（影響因子序號），j＝1，2，…，m （各影響因子Ⅱ級分類因子序號）。

按照式（1）計算出一個初始指標值xi，j，然后再把xi，j進行標準化處理。標準化處理方法如下：

式中：Ii，j——各類Ⅱ級因子的標準化值。

2.3 Logistic回歸分析方法

Logistic回歸模型是二值響應變量（因變量）對回歸變量（自變量）的分析，其結(jié)果是預測某種事件的發(fā)生概率。在泥石流易發(fā)性評價過程中，各因子數(shù)據(jù)作為自變量，泥石流的發(fā)生與否作為分類因變量（0代表泥石流不發(fā)生，1代表泥石流發(fā)生），實現(xiàn)泥石流易發(fā)性評價模型的構建。

Logistic回歸模型中的響應變量為二值變量，記為Y，只取2個值：0和1。令P為泥石流發(fā)生的概率，取值范圍（0，1），則1－P為泥石流不發(fā)生的概率。對于n個泥石流影響因子X1，X2，…，Xn，當X1＝x1，X2＝x2，…，Xn＝xn時，Y＝1的泥石流發(fā)生概率為P（YH＝1│X1＝x1，X2＝x2，…，Xn＝xn）?？紤]Logistic變換［19］，以災害發(fā)生的概率P為因變量，影響因子集X1，X2，…，Xn為自變量，建立模型：

對（3）式進行等價變換：

模型中C0是常數(shù)項，C1，C2，…，Cn是回歸系數(shù)。通過對已知樣本的求解，即可確定回歸系數(shù)。

式（4）中，如果邏輯回歸系C1，C2，…，Cn確定，則根據(jù)不同的指標值xi，j，計算研究區(qū)泥石流發(fā)生概率P值。根據(jù)P值大小，劃分泥石流發(fā)生可能性等級，完成泥石流易發(fā)性分區(qū)。

3 易發(fā)性評價模型

依據(jù)研究區(qū)1∶5萬DEM 數(shù)據(jù)，1∶25萬地質(zhì)圖，全國土地利用類型數(shù)據(jù)以及沿線泥石流數(shù)據(jù)，分析各個Ⅱ級分類因子的面積和泥石流分布數(shù)量，采用公式（1）計算Ⅱ級分類因子的指標值，再依據(jù)公式（2）對各類評價因子進行標準化處理，獲取統(tǒng)一尺度的量化指標值。各個影響因子的分類分級情況和計算結(jié)果如表1所示。

采用各關鍵因子數(shù)據(jù)類的標準化值作為邏輯回歸模型的自變量，泥石流易發(fā)概率作為因變量，依據(jù)Binary Logistic Regression分析，確定各個關鍵因子與泥石流發(fā)生的相關系數(shù)，即為二項邏輯回歸模型中各個關鍵因子參數(shù)的系數(shù)。Logistic回歸分析結(jié)果中，除坡度的系數(shù)所對應的sig值（0.076）略大于0.05外，其余均小于0.05，可以認為分析結(jié)果通過5%顯著性水平檢驗。因此，依據(jù)各個因子系數(shù)值，構建Logistic回歸模型：

式中：IELE，ISLO，IROU，ILIT，IFAU，IEAR，ILAN——高程、坡度、坡面粗糙度，地層巖性、距斷層距離，地震烈度，土地利用的因子參數(shù)。

4 結(jié)果與分析

根據(jù)構建的Logistic回歸模型，用式（6）計算每一網(wǎng)格單元各個影響因子的值，再用式（5）計算綜合評價圖層中每一單元的泥石流災害易發(fā)性指數(shù)（發(fā)生概率）。依據(jù)泥石流易發(fā)性指數(shù)計算結(jié)果，按照自然斷點法將研究區(qū)分為微易發(fā)區(qū)、低易發(fā)區(qū)、中易發(fā)區(qū)、高易發(fā)區(qū)等4級，并制作泥石流易發(fā)性評價分區(qū)圖（圖2）。

分析結(jié)果表明，泥石流高易發(fā)區(qū)的面積為12 655.51km2，占總面積的36.22%；泥石流中易發(fā)區(qū)的面積為14 812.66km2，占總面積的42.39%；泥石流低易發(fā)區(qū)的面積為6 071.65km2，占總面積的17.38%；泥石流微易發(fā)區(qū)的面積為1 402.92km2，占總面積的4.01%。從公路全線范圍看，大渡河、瀾滄江、金沙江、怒江、帕隆藏布江等大江大河的峽谷區(qū)段，泥石流易發(fā)性指較大，而高原面與成都平原區(qū)的公路沿線泥石流易發(fā)性指數(shù)較小。

表1 各評價因子分級與量化指標值

圖2 G318川藏公路段泥石流易發(fā)性評價分區(qū)

將研究區(qū)內(nèi)934條歷史泥石流災害位置與泥石流易發(fā)性區(qū)劃結(jié)果進行對比，發(fā)現(xiàn)共有493條泥石流位于高易發(fā)區(qū)，占災害總數(shù)的52.8%；中易發(fā)區(qū)有385條泥石流，占災害總數(shù)的41.2%；微易發(fā)區(qū)沒有泥石流發(fā)生（表2）。上述分析說明，本文的易發(fā)性評價結(jié)果具有可靠性。

表2 泥石流易發(fā)性評價結(jié)果

5 結(jié)論

分析了泥石流形成與演化過程中孕災因子對泥石流易發(fā)性的空間關系，建立了一套結(jié)合GIS和Logistic回歸模型的交通干線泥石流易發(fā)性定量分析與制圖流程包括：（1）建立基礎數(shù)據(jù)庫，選取關鍵影響因子；（2）劃分評價單元，計算各評價因子的易發(fā)性指標值；（3）建立評價模型，計算綜合易發(fā)指數(shù)；（4）進行泥石流易發(fā)性分區(qū)。

本文選取坡度、地面高程、坡面粗糙度、地層巖性、距斷裂距離、地震烈度、土地利用類型為泥石流易發(fā)性評價因子，引入Logistic回歸分析方法，構建泥石流易發(fā)性評價模型。定量分析G318川藏公路段泥石流易發(fā)性，結(jié)果表明，泥石流高易發(fā)區(qū)與中易發(fā)區(qū)占總面積的78.61%，主要分布在大渡河、瀾滄江、金沙江、怒江、帕隆藏布江等大江大河的峽谷區(qū)段；而微易發(fā)區(qū)僅占總面積的4.01%，主要分布在高原面與成都平原地區(qū)。

中國的長大交通干線，如川藏公路，青藏鐵路，沿途穿過多個差異巨大的地貌單元，地勢高差和氣候差異顯著，地質(zhì)條件復雜。道路安全不僅面臨高山峽谷區(qū)、山地災害頻發(fā)區(qū)等災害環(huán)境，還面臨著全球變化誘發(fā)超常規(guī)模地表災害，具有巨大的工程風險。為此，確定交通干線泥石流災害主控因素，精確判定泥石流發(fā)生與危害程度，建立交通干線泥石流災害風險分析與預警平臺意義重大，此方面的研究工作還需進一步完善。

［1］ 唐邦興，柳素清，劉世建.中國泥石流分布及其災害危險區(qū)劃圖［M］.成都：成都地圖出版社，1991.

［2］ 崔鵬，陳曉清，程尊蘭，等.西藏泥石流滑坡監(jiān)測與防治［J］.自然雜志，2010，32（1）：19－25.

［3］ Carrara A.Multivariate models for landslide hazard evaluation［J］.Mathematical Geology，1983，15（3）：403－426.

［4］ Montgomery D R，Dietrich W E.A physically based model for the topographic control on shallow landsliding［J］.Water Resource Research，1994，30（4）：1157－1171.

［5］ Chung C F，F(xiàn)abbri A G.Probabilistic Prediction models for landslide hazard mapping［J］.Photogrammetric Engineering＆ Remote Sensing，1999，65（12）：1388－1399.

［6］ Cascini L.Applicability of landslide susceptibility and hazard zoning at different scales［J］.Engineering Geology，2008，102（3）：164－177.

［7］ Rapolla A，Nocera S D，Matano F，et al.Susceptibility regional zonation of earthquake－induced landslides in Campania，Southern Italy［J］.Natural Hazards，2012，61（1）：115－126.

［8］ 唐川.云南怒江流域泥石流敏感性空間分析［J］.地理研究，2005，24（2）：178－185.

［9］ 鄒強，崔鵬，張建強，等.長江上游地區(qū)泥石流災害敏感性量化評價研究［J］.環(huán)境科學與技術，2012，35（3）：159－167.

［10］ 羅德富，朱平一，陳瑞，等.川藏公路南線（西藏境內(nèi)）山地災害及防治對策［M］.北京：科學出版社，1995.

［11］ 蔣忠信.西藏帕隆藏布河谷崩塌滑坡、泥石流的分布規(guī)律［J］.地理研究，2002，21（4）：495－503.

［12］ 何易平，胡凱衡，韋方強，等.川藏公路迫隆藏布流域段泥石流活動特征［J］.水土保持學報，2001，15（3）：76－80.

［13］ 游勇，程尊蘭，胡平華，等.西藏古鄉(xiāng)溝泥石流模型試驗研究［J］.自然災害學報，1997，6（1）：52－58.

［14］ 朱平一，程尊蘭.川藏公路培龍溝泥石流輸砂堵江成因探討［J］.自然災害學報，2000，9（1）：80－83.

［15］ 李德基，游勇.西藏波密米堆冰湖潰決淺議［J］.山地學報，1992，10（4）：219－224.

［16］ 朱阿興，裴韜，喬建平，等.基于專家知識的滑坡危險性模糊評估方法［J］.地理科學進展.2006，25（4）：1－12.

［17］ Ouri A E，Amirian S.Landslide hazard zonation using MR and AHP methods and GIS techniques in Langan watershed，Ardabil，Iran［C］∥International conference on ACRS.Beijing，China，2009：19－22.

［18］ Ayalew L，Yamagishi H，Ugawa N.Landslide susceptibility mapping using GIS－based weighted linear combination：The case in Tsugawa area of Agano River，Niigata Prefecture，Japan［J］.Landslides，2004（1）：73－81.

［19］ 蘭恒星，伍法權，周成虎，等.基于GIS的云南小江流域滑坡因子敏感性分析［J］.巖石力學與工程學報，2002，21（10）：1500－1506.

［20］ Lee S，Ryu J H，Lee M J，et al.Use of an artificial neural network for analysis of the susceptibility to landslides at Boun，Korea［J］.Environmental Geology，2003，44（7）：820－833.

［21］ Ayalew L，Yamagishi H.The application of GIS－based logistic regression for landslide susceptibility mapping in the Kakuda—Yahiko Mountains，Central Japan［J］.Geomorphology，2005，65（1）：15－31.

［22］ 王濟川，郭志剛.Logistic回歸模型：方法與應用［M］.北京：高等教育出版社，2001.