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

鄒 強(qiáng)，王 青，劉延國(guó)

（1.西南科技大學(xué) 環(huán)境與資源學(xué)院，四川 綿陽(yáng)621010；2.中國(guó)科學(xué)院 水利部成都山地災(zāi)害與環(huán)境研究所 中國(guó)科學(xué)院 山地災(zāi)害與地表過(guò)程重點(diǎn)實(shí)驗(yàn)室，四川 成都610041）

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

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

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

1 研究區(qū)概況

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

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

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

2 研究方法

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

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

2.1 泥石流影響因子確定

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

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

2.2 因子量化方法

在泥石流易發(fā)性分析過(guò)程中，由于各個(gè)評(píng)價(jià)因子的量綱不同，各個(gè)評(píng)價(jià)因子作為回歸模型中的自變量首先需要進(jìn)行標(biāo)準(zhǔn)化處理，以獲取同一尺度的量化指標(biāo)值，便于進(jìn)行后續(xù)分析與計(jì)算。本文用每個(gè)影響因子的各個(gè)Ⅱ級(jí)分類(lèi)因子中泥石流災(zāi)害點(diǎn)的個(gè)數(shù)與各個(gè)Ⅱ級(jí)分類(lèi)因子的面積的比值作為各個(gè)Ⅱ級(jí)分類(lèi)因子的指標(biāo)值。計(jì)算表達(dá)式如下。

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

按照式（1）計(jì)算出一個(gè)初始指標(biāo)值xi，j，然后再把xi，j進(jìn)行標(biāo)準(zhǔn)化處理。標(biāo)準(zhǔn)化處理方法如下：

式中：Ii，j——各類(lèi)Ⅱ級(jí)因子的標(biāo)準(zhǔn)化值。

2.3 Logistic回歸分析方法

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

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

對(duì)（3）式進(jìn)行等價(jià)變換：

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

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

3 易發(fā)性評(píng)價(jià)模型

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

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

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

4 結(jié)果與分析

根據(jù)構(gòu)建的Logistic回歸模型，用式（6）計(jì)算每一網(wǎng)格單元各個(gè)影響因子的值，再用式（5）計(jì)算綜合評(píng)價(jià)圖層中每一單元的泥石流災(zāi)害易發(fā)性指數(shù)（發(fā)生概率）。依據(jù)泥石流易發(fā)性指數(shù)計(jì)算結(jié)果，按照自然斷點(diǎn)法將研究區(qū)分為微易發(fā)區(qū)、低易發(fā)區(qū)、中易發(fā)區(qū)、高易發(fā)區(qū)等4級(jí)，并制作泥石流易發(fā)性評(píng)價(jià)分區(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%。從公路全線(xiàn)范圍看，大渡河、瀾滄江、金沙江、怒江、帕隆藏布江等大江大河的峽谷區(qū)段，泥石流易發(fā)性指較大，而高原面與成都平原區(qū)的公路沿線(xiàn)泥石流易發(fā)性指數(shù)較小。

表1 各評(píng)價(jià)因子分級(jí)與量化指標(biāo)值

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

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

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

5 結(jié)論

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

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

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

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

［2］ 崔鵬，陳曉清，程尊蘭，等.西藏泥石流滑坡監(jiān)測(cè)與防治［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］ 鄒強(qiáng)，崔鵬，張建強(qiáng)，等.長(zhǎng)江上游地區(qū)泥石流災(zāi)害敏感性量化評(píng)價(jià)研究［J］.環(huán)境科學(xué)與技術(shù)，2012，35（3）：159－167.

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

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

［12］ 何易平，胡凱衡，韋方強(qiáng)，等.川藏公路迫隆藏布流域段泥石流活動(dòng)特征［J］.水土保持學(xué)報(bào)，2001，15（3）：76－80.

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

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

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

［16］ 朱阿興，裴韜，喬建平，等.基于專(zhuān)家知識(shí)的滑坡危險(xiǎn)性模糊評(píng)估方法［J］.地理科學(xué)進(jìn)展.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］ 蘭恒星，伍法權(quán)，周成虎，等.基于GIS的云南小江流域滑坡因子敏感性分析［J］.巖石力學(xué)與工程學(xué)報(bào)，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］ 王濟(jì)川，郭志剛.Logistic回歸模型：方法與應(yīng)用［M］.北京：高等教育出版社，2001.