羅 健，郭 寧，楊 崢

(1.西昌學(xué)院資源與環(huán)境學(xué)院，四川 西昌 615013；2.西昌學(xué)院土木與水利工程學(xué)院，四川 西昌 615013)

0 引言

溝谷泥石流形成受地形、地質(zhì)、地貌、氣象、水文、植被、土壤等諸多不確定性因素的共同控制和影響，是典型的復(fù)雜系統(tǒng)。國(guó)外泥石流危險(xiǎn)度研究開展比較早[1-2]，而地理信息系統(tǒng)(GIS)技術(shù)也很早就應(yīng)用在了泥石流危險(xiǎn)度評(píng)價(jià)中[3]。20世紀(jì)80年代末90年代初，我國(guó)學(xué)者劉希林等先后提出了判定泥石流危險(xiǎn)度的指標(biāo)選取、權(quán)重確定[4]、以及區(qū)劃方法[5]，并指出了山區(qū)溝谷泥石流主要危險(xiǎn)因子[6]。唐川等[7]應(yīng)用泥石流二維非恒定流理論建立危險(xiǎn)度評(píng)價(jià)模型用于泥石流堆積泛濫區(qū)危險(xiǎn)區(qū)劃。

GIS憑借快速精準(zhǔn)獲取研究區(qū)范圍、數(shù)據(jù)矢量化和空間特征測(cè)度等方面空間分析技術(shù)優(yōu)勢(shì)廣泛應(yīng)用于泥石流危險(xiǎn)度評(píng)估研究之中[8]。早期的應(yīng)用中主要依據(jù)知識(shí)庫對(duì)因子進(jìn)行分級(jí)和賦權(quán)[9]，采用直接指標(biāo)法和間接指標(biāo)法[10]，易損性評(píng)價(jià)圖與危險(xiǎn)性評(píng)價(jià)圖疊加[11]等進(jìn)行分析，并集合各種評(píng)價(jià)模型建立評(píng)價(jià)決策系統(tǒng)[12]。隨著研究的進(jìn)一步深入，神經(jīng)網(wǎng)絡(luò)[13]、投影尋蹤聚類模型[14]、邏輯回歸模型[15]、分形維數(shù)[16]、信息熵[17]，主成分分析(PCA)[18]，層次分析法(AHP)[19]等各種評(píng)價(jià)模型和方法不斷應(yīng)用到泥石流危險(xiǎn)度評(píng)價(jià)中。為了克服研究者主觀賦權(quán)的影響，還出現(xiàn)了綜合兩種模型方法的組合賦權(quán)規(guī)則[20]。

西昌市位于安寧河流域中下游，地質(zhì)構(gòu)造復(fù)雜，境內(nèi)以牦牛山、螺髻山、小相嶺山地地貌，安寧河沖積平原，邛海斷陷盆地為主；屬亞熱帶季風(fēng)氣候，年平均降雨1 013.5 mm，降雨集中于6～9月，占全年降水量的76%，雨型多為暴雨，日最大降雨量達(dá)135.7 mm。據(jù)西昌市地質(zhì)災(zāi)害隱患調(diào)查評(píng)價(jià)項(xiàng)目，西昌境內(nèi)現(xiàn)發(fā)育地質(zhì)災(zāi)害點(diǎn)220處，其中泥石流100處、滑坡101處、不穩(wěn)定斜坡19處，每年雨季都有地質(zhì)災(zāi)害爆發(fā)并造成不同程度的災(zāi)情和損失，地質(zhì)災(zāi)害已經(jīng)成為西昌市自然災(zāi)害中破壞性最大的災(zāi)害之一。研究區(qū)地處安寧河斷裂帶地質(zhì)活動(dòng)較為活躍的東支，屬逆沖走滑斷裂帶[21]；同時(shí)該地區(qū)也是西昌市北部工業(yè)園區(qū)和未來西昌城市發(fā)展(北部新城)的重要拓展方向。對(duì)研究區(qū)溝谷泥石流危險(xiǎn)度進(jìn)行評(píng)價(jià)可以為該地區(qū)相關(guān)規(guī)劃建設(shè)和防災(zāi)減災(zāi)工作提供決策依據(jù)。

1 數(shù)據(jù)來源及研究方法

1.1 研究區(qū)范圍及數(shù)據(jù)來源

泥石流危險(xiǎn)度評(píng)價(jià)單元通常有5種:柵格單元、地貌單元、均一條件單元、流域單元及地形單元[22]。鑒于研究區(qū)地形地貌特點(diǎn)，并考慮到泥石流發(fā)育形成過程與溝谷環(huán)境特征密切相關(guān)，本次評(píng)價(jià)采用流域作為評(píng)價(jià)單元。泥石流形成機(jī)理往往包括物源、流通、觸發(fā)3類因素，觸發(fā)因素除了降水外還包括地質(zhì)構(gòu)造運(yùn)動(dòng)和人類活動(dòng)，3類影響因素根據(jù)性質(zhì)的不同可細(xì)分為地質(zhì)、地形、氣象、水文、植被和人類活動(dòng)等。結(jié)合溝谷泥石流物質(zhì)形成、堆積、搬移和再次堆積的過程，以及我國(guó)西南山區(qū)地理環(huán)境特點(diǎn)本次評(píng)價(jià)選取了X1主溝長(zhǎng)度、X2平均坡度、X3流域高差、X4流域面積、X5近構(gòu)造線面積、X6巖性面積、X7森林覆蓋率、X8雨季月均降水量共8項(xiàng)評(píng)價(jià)指標(biāo)。其中，流域面積、近構(gòu)造線面積和巖性面積主要表征物質(zhì)形成來源；主溝長(zhǎng)度、平均坡度和流域高差表征松散物質(zhì)堆積和搬移條件；雨季月均降水量則是重要的觸發(fā)因素；而森林覆蓋率既體現(xiàn)了物質(zhì)搬移條件，又體現(xiàn)了物質(zhì)形成來源。研究區(qū)范圍基于DEM通過GIS水文分析生成分水嶺和集水區(qū)，然后根據(jù)水流匯集方向進(jìn)行合并(見圖1)，共計(jì)21個(gè)溝谷單元。由于安寧河流域西昌北段地處安寧河河谷開闊平坦處，寬度達(dá)到8～12 km，所以水文分析過程中研究區(qū)只選取盆地邊緣山區(qū)溝谷部分(見圖2)，主要是泥石流形成和通過區(qū)。本研究主要指標(biāo)數(shù)據(jù)由DEM數(shù)據(jù)生成，此外植被覆蓋數(shù)據(jù)由遙感解譯得到，地質(zhì)構(gòu)造線及巖性數(shù)據(jù)來自1∶50 000地質(zhì)圖，由涼山州國(guó)土局地質(zhì)環(huán)境監(jiān)測(cè)站(涼山州突發(fā)地質(zhì)災(zāi)害應(yīng)急會(huì)商指揮中心)提供。雨季月平均降水?dāng)?shù)據(jù)為國(guó)家生態(tài)系統(tǒng)觀測(cè)研究網(wǎng)絡(luò)1970年至2000年1 km×1 km降水柵格數(shù)據(jù)集[23-24]。

圖1 研究區(qū)域溝谷單元Fig.1 Catchment unit in the study area

圖2 研究區(qū)域DEM地形圖Fig.2 DEM in the study area

1.2 主成分分析方法

泥石流形成和發(fā)育受到地質(zhì)、地形、植被、水文等多方面因素綜合影響作用，具備復(fù)雜系統(tǒng)特征，在分析評(píng)價(jià)時(shí)需要對(duì)相互關(guān)聯(lián)的原始變量進(jìn)行綜合化簡(jiǎn)，同時(shí)確保盡量多保留原有信息。所以，本次泥石流危險(xiǎn)度評(píng)價(jià)采用主成分分析法(PCA)[25]。主成分分析是一種適用的降維處理技術(shù)，通過對(duì)原始變量的線性組合，新生成主成分相互獨(dú)立且方差累積貢獻(xiàn)量足夠大。其基本分析步驟如下：

(1)計(jì)算相關(guān)系數(shù)矩陣

通過相關(guān)系數(shù)公式計(jì)算原始變量xi與xj之間的相關(guān)系數(shù)矩陣rij(i,j=1,2,...,p)得到：

(2)計(jì)算特征值和特征向量

通常用雅克比法求解特征值λi(i=1,2,…,p)，并按大小進(jìn)行排序，特征方程如下：

λI-R=0

然后分別求出對(duì)應(yīng)于特征值λi的特征向量ei(i=1,2,…,p)，計(jì)算公式為：

其中eij是向量ei的第j個(gè)分量。

(3)計(jì)算主成分累積貢獻(xiàn)率

一般取特征值大于1，且累積貢獻(xiàn)率大于80%的特征值λ1,λ2,…,λm所對(duì)應(yīng)的第一，第二，...，第m(m≤p)個(gè)主成分。

(4)計(jì)算主成分載荷及綜合得分

主成分載荷計(jì)算公式為

所以各主成分得分可以按下式計(jì)算

在計(jì)算過程中為了消除量綱的影響，原始變量xp往往采用Z-score標(biāo)準(zhǔn)化處理，使得標(biāo)準(zhǔn)化后的變量均值為0，標(biāo)準(zhǔn)差為1。但是，Z-score標(biāo)準(zhǔn)化方法參照系為樣本區(qū)的均值和標(biāo)準(zhǔn)差，這只能體現(xiàn)樣本區(qū)各指標(biāo)變量的相對(duì)離散程度，不同研究區(qū)之間不具有可比性。為了消除量綱并體現(xiàn)各指標(biāo)變量的絕對(duì)差異，本文采用10級(jí)分類賦值(見表1)。綜合得分Z由取得的前m個(gè)主成分加權(quán)得到，權(quán)重選取各主成分貢獻(xiàn)率占累積貢獻(xiàn)率的比重。

表1 參評(píng)變量分級(jí)賦值表

2 變量提取及危險(xiǎn)度評(píng)價(jià)

主溝矢量化后分別計(jì)算長(zhǎng)度、直線距離和比降；坡度由DEM數(shù)據(jù)直接生成(圖3)，再以溝谷單元統(tǒng)計(jì)各流域平均坡度；森林覆蓋率由遙感影像配準(zhǔn)后目視解譯獲取(圖4)。流域面積由水文分析生成的集水區(qū)結(jié)合河網(wǎng)和遙感影像進(jìn)行合并后測(cè)算得到；構(gòu)造線和巖性通過對(duì)1∶50 000地質(zhì)圖配準(zhǔn)拼接后矢量化生成。近構(gòu)造線面積由距離構(gòu)造線250 m、500 m，750 m、1 000 m、1 250 m多環(huán)緩沖區(qū)生成，并分別按系數(shù)1.0、0.8、0.6、0.4、0.2縮減后加和匯總(圖5)。巖性基于普氏系數(shù)分為5級(jí)(表2、圖6)，川西地區(qū)堅(jiān)固系數(shù)小于8的軟巖地區(qū)主要暴發(fā)中高頻率泥石流[26]，所以1至5級(jí)分別按系數(shù)0.2、0.4、0.6、0.8、1.0縮減后加和匯總。雨季月平均降水?dāng)?shù)據(jù)采用西昌雨季6月至9月1 km×1 km降水柵格數(shù)據(jù)集，按各溝谷分區(qū)統(tǒng)計(jì)均值(圖7)。

表2 基于普氏系數(shù)的巖性堅(jiān)固(軟硬)分級(jí)

圖3 研究區(qū)域坡度分異Fig.3 Slope classification in the study area

圖4 研究區(qū)域森林覆蓋Fig.4 Forest cover classification in the study area

圖5 距離構(gòu)造線遠(yuǎn)近范圍Fig. 5 The distance classification from from the tectonic line

圖6 基于普氏系數(shù)分級(jí)的巖性類型Fig.6 The lithology classification in the study area

圖7 1971～2000年6月每平方公里平均降水量Fig.7 The average rainfall per square kilometer in June from 1971 to 2000

各溝谷單元評(píng)價(jià)指標(biāo)原始變量值見表3。主成分分析法通過對(duì)原始變量的線性組合進(jìn)行降維處理，但是泥石流災(zāi)害評(píng)價(jià)指標(biāo)中的坡度與危險(xiǎn)度卻不是簡(jiǎn)單的線性關(guān)系，因?yàn)槠露瘸^一定的限度之后松散物質(zhì)便很難堆積積累，所以對(duì)坡度指標(biāo)直接標(biāo)準(zhǔn)化處理會(huì)對(duì)載荷計(jì)算和多變量之間的共線性估計(jì)造成影響。為了避免上述問題，采用上文的各變量分級(jí)賦值進(jìn)行主成分分析；同時(shí)，為了解決同一賦值區(qū)間原始變量差異的消除，本文在分級(jí)賦值的基礎(chǔ)上再借助區(qū)間端點(diǎn)進(jìn)行線性插值。具體按下式進(jìn)行計(jì)算：

式中:Xij′——標(biāo)準(zhǔn)化之后的值;

Fl——所在分級(jí)區(qū)間賦值的下限值;

Fh——所在分級(jí)區(qū)間賦值的上限值;

Xij——指標(biāo)的原始值;

Xl——分級(jí)區(qū)間的下限值;

Xh——分級(jí)區(qū)間的上限值。

以賀家溝主溝長(zhǎng)度1.81 km為例，結(jié)合表1則有：

表3 各流域評(píng)價(jià)指標(biāo)原始值

原始變量分級(jí)線性插值后見表4，第10級(jí)區(qū)間上限值有正無窮大的在實(shí)際計(jì)算過程中替換為下限值的2倍值，如果實(shí)際值超過了下限值的兩倍，分級(jí)線性插值直接取值為10。

表4 各流域評(píng)價(jià)指標(biāo)原始變量分級(jí)線性插值

經(jīng)過線性插值處理的指標(biāo)兩兩計(jì)算相關(guān)系數(shù)，生成相關(guān)系數(shù)矩陣見表5。此外， KMO檢驗(yàn)值達(dá)0.838，大于0.6，并通過了Barrlett球形檢驗(yàn)，說明適合采用主成分分析方法進(jìn)行綜合評(píng)價(jià)。

表5 變量相關(guān)系數(shù)矩陣

累積方差貢獻(xiàn)率表中第一主成分和第二主成分特征值大于1(表6)，并且累積貢獻(xiàn)率達(dá)到87.671%，大于80%，符合提取條件，故提取前兩個(gè)主成分。

表6 累積方差貢獻(xiàn)率

在輸出的第一、二主成分向量矩陣基礎(chǔ)上通過主成分載荷公式計(jì)算各主成分載荷和特征向量。在此基礎(chǔ)上，通過各主成分貢獻(xiàn)率占累積貢獻(xiàn)率比重進(jìn)一步歸一化求取權(quán)重(表7)。

表7 載荷矩陣、特征向量及變量權(quán)重Table 7 The load matrix, characteristic vectors and variable weights

再按如下公式計(jì)算主成分得分：

Z=0.1573×ZX1+ 0.1289×ZX2+0.1413×ZX3+0.1530×ZX4+ 0.1550×ZX5+ 0.1513×ZX6-0.0164×ZX7+0.1295×ZX8

危險(xiǎn)度綜合得分見表8。泥石流爆發(fā)危險(xiǎn)度最大值為10，研究區(qū)中約一半的溝谷單元危險(xiǎn)度得分在5分以上，形勢(shì)較為嚴(yán)峻，其中熱水河分值高達(dá)9.24。

表8 各流域單元主成分得分及加權(quán)綜合得分

3 泥石流危險(xiǎn)度評(píng)價(jià)結(jié)果與實(shí)際災(zāi)害點(diǎn)對(duì)比分析

危險(xiǎn)度評(píng)價(jià)得分結(jié)果常見分級(jí)方法中，自然斷裂點(diǎn)和分位數(shù)都表征的是數(shù)值的相對(duì)離散程度，不同時(shí)間節(jié)點(diǎn)和不同研究區(qū)域之間也都沒有可比性。所以本次評(píng)價(jià)危險(xiǎn)度綜合得分結(jié)果采用絕對(duì)閾值劃分危險(xiǎn)程度，并使得每個(gè)危險(xiǎn)度等級(jí)數(shù)值等距為2(分級(jí)結(jié)果見表9)。

表9 泥石流危險(xiǎn)度等級(jí)表

危險(xiǎn)度分級(jí)結(jié)果見圖8，結(jié)果圖中疊置了西昌市地質(zhì)災(zāi)害隱患調(diào)查評(píng)價(jià)項(xiàng)目成果表中泥石流空間分布，并按照DD 2008-02(1∶50 000滑坡崩塌泥石流災(zāi)害調(diào)查規(guī)范)分為小型、中型、大型3級(jí)。較危險(xiǎn)和極危險(xiǎn)區(qū)溝谷內(nèi)部和溝口堆積區(qū)都發(fā)生了泥石流災(zāi)害，其中熱水河、大塘河與下黑砂河流域內(nèi)部分布有一處大型泥石流。7個(gè)臨界危險(xiǎn)區(qū)中有4個(gè)流域內(nèi)部和匯流通道前端堆積區(qū)形成了泥石流災(zāi)害。其中，比較特殊是上黑砂河，其內(nèi)部有一處大型泥石流。通過實(shí)地踏勘調(diào)查，該流域曾是我國(guó)典型的泥石流災(zāi)害點(diǎn)，中國(guó)科學(xué)院第二地質(zhì)調(diào)查大隊(duì)做過長(zhǎng)期觀測(cè)，并由成都地理所設(shè)計(jì)支持，西昌市和涼山州政府相關(guān)部門經(jīng)過進(jìn)行治理，通過水庫、攔沙壩、谷坊、順?biāo)畨巍⒆o(hù)岸、導(dǎo)流堤、分洪捻、排洪溝等工程措施和植樹造林使得災(zāi)情得到了很好的控制[27]。此外，僅馬家溝一處評(píng)價(jià)等級(jí)為較安全，仍在其堆積區(qū)形成了一處泥石流災(zāi)害點(diǎn)。綜上，研究區(qū)泥石流危險(xiǎn)度評(píng)價(jià)結(jié)果和實(shí)際災(zāi)害點(diǎn)分布較為符合。

圖8 研究區(qū)域泥石流災(zāi)害點(diǎn)分布及危險(xiǎn)度比對(duì)Fig.8 The hierarchical results of debris flow hazardous degree and distribution of debris flow hazard points in the study area備注：泥石流災(zāi)害點(diǎn)圖示符號(hào)中含*號(hào)的是指已經(jīng)采用工程治理措施的災(zāi)害點(diǎn)。

4 結(jié)論與討論

基于安寧河流域西昌北段東岸溝谷環(huán)境特征以及相關(guān)研究成果，采用流域作為評(píng)價(jià)單元，并選取了X1主溝長(zhǎng)度、X2平均坡度、X3流域高差、X4流域面積、X5近構(gòu)造線面積、X6巖性面積、X7森林覆蓋率、X8雨季月均降水量共8項(xiàng)評(píng)價(jià)指標(biāo)。評(píng)價(jià)結(jié)果顯示，安寧河流域西昌北段東岸中一半的溝谷單元危險(xiǎn)度得分在5分以上，形勢(shì)較為嚴(yán)峻，其中熱水河分值高達(dá)9.24。疊置研究區(qū)地質(zhì)災(zāi)害隱患調(diào)查中的泥石流災(zāi)害點(diǎn)空間分布后發(fā)現(xiàn)，較危險(xiǎn)和極危險(xiǎn)區(qū)溝谷內(nèi)部和溝口堆積區(qū)都發(fā)生了泥石流災(zāi)害，其中熱水河、大塘河與下黑砂河流域內(nèi)部分布有一處大型泥石流。7個(gè)臨界危險(xiǎn)區(qū)中有4個(gè)流域內(nèi)部和匯流通道前端堆積區(qū)形成了泥石流災(zāi)害。通過實(shí)地踏勘驗(yàn)證，上黑砂河內(nèi)部有一處大型泥石流，但是評(píng)價(jià)等級(jí)為較安全區(qū)，主要得益于長(zhǎng)期實(shí)施的工程和生物治理措施。此外，僅馬家溝一處評(píng)價(jià)等級(jí)為較安全，仍在其堆積區(qū)形成了一處泥石流災(zāi)害點(diǎn)。綜上，研究區(qū)泥石流危險(xiǎn)度評(píng)價(jià)結(jié)果和實(shí)際災(zāi)害點(diǎn)分布較為符合。

在溝谷泥石流危險(xiǎn)度評(píng)價(jià)中，主成分分析法通過對(duì)原始變量的線性組合進(jìn)行降維處理，一定程度上解決了變量間的多重共線性；同時(shí)，絕對(duì)閾值分級(jí)賦值刻畫了各變量絕對(duì)差異，可以跨越區(qū)域和時(shí)間節(jié)點(diǎn)進(jìn)行比較。當(dāng)涉及與危險(xiǎn)度不是簡(jiǎn)單的線性關(guān)系的指標(biāo)，比如坡度，直接標(biāo)準(zhǔn)化處理會(huì)對(duì)載荷計(jì)算和多變量之間的共線性估計(jì)造成影響。為了解決這個(gè)問題，本文在分級(jí)賦值的基礎(chǔ)上再借助區(qū)間端點(diǎn)進(jìn)行線性插值。需要注意的是，第10級(jí)區(qū)間上限值有正無窮大的在實(shí)際計(jì)算過程中替換為下限值的2倍值，如果實(shí)際值超過了下限值的兩倍，分級(jí)線性插值直接取值為10。此外，負(fù)相關(guān)指標(biāo)在進(jìn)行分級(jí)線性插值的時(shí)候賦值區(qū)間的上下限值應(yīng)當(dāng)互換，比如植被覆蓋率。最后，危險(xiǎn)度評(píng)價(jià)結(jié)果與泥石流災(zāi)害點(diǎn)調(diào)查分級(jí)結(jié)果的疊置分析有助于兩者相互驗(yàn)證。對(duì)于危險(xiǎn)度等級(jí)較低但是發(fā)生了較大規(guī)模泥石流災(zāi)害的要通過實(shí)地踏勘是否是治理效果的體現(xiàn)，或是流域集水區(qū)劃分有誤；對(duì)于危險(xiǎn)度等級(jí)較高但是沒有發(fā)現(xiàn)災(zāi)害點(diǎn)應(yīng)當(dāng)引起重視，實(shí)地踏勘驗(yàn)證是否有泥石流災(zāi)害將要發(fā)生的跡象，盡早采用監(jiān)測(cè)和防治措施。