川鄂褶皺山地溪洪-滑坡災害與主要自然因子的關系
——以香溪河流域為例

杜 俊, 丁文峰, 范仲杰, 李清溪

(1.長江科學院 水土保持研究所， 湖北 武漢 430010； 2.水利部 山洪地質災害防治工程技術研究中心， 湖北 武漢 430010)

三峽庫區是中國典型的生態環境脆弱區，2003年三峽水庫蓄水以來，庫區滑坡問題愈發突出，大量典型案例的調查、觀測和試驗表明，水庫水位的周期性升降會改變滑體靜/動水壓力、降低巖土體抗剪強度，進而促進滑體變形、滑帶貫通以及最終的累積釋放[1-4]。防汛部門的資料顯示，早在三峽水庫蓄水以前，庫區滑坡即已十分發育，是該區山洪地質災害(即溪河洪水以及由其引發的滑坡、泥石流等)[5]的主要表現形式。顯然水庫蓄水只是外因，庫區特有的自然地理和地質條件才是滑坡發育的實質，因此只有弄清當地滑坡發育格局與具體根植性要素的關系，才可能從根本上認知區域滑坡宏觀發育機制和開展更有針對性的防治工作。以往研究對于這一問題已有一定積累，鄧清祿[6]和陳劍[7]等指出庫區滑坡歷史上主要與新構造期以來的快速抬升和間冰期極大降水有關，因此許多大型滑坡在長江干流和主要支流深切方向上沿程分布；白世彪等[8]基于GIS技術和頻率統計，對三峽水庫175 m回水影響區的滑坡分布進行了統計，總結了滑坡分布與巖性巖組、高程、坡度、坡向、曲率等的敏感性特征；喬建平等[9]使用本底因子貢獻率法統計了地層、坡度、坡形、高差和坡向對三峽庫區云陽—巫山段滑坡分布的敏感性特征，并認為高差、坡度和地層的危險度權重更高；Wang和Niu[10]應用中巴地球資源衛星遙感影像，采集植被覆蓋、坡體結構、水庫水位、高程、工程巖組等20個指標信息，基于決策數方法對秭歸縣郭家壩鎮的滑坡分布進行了高精度預報；李雪等[11]分析了三峽庫首區滑坡分布與12項地形因子的相關性，認為高程、坡向、高程變異系數與滑坡的相關性較高；張俊等[12]總結的萬州區滑坡發育主要因子為地層巖性、地質構造、水系分布、坡度、坡向、坡體結構及土地利用。綜觀這些工作，大部分學者都十分重視統計滑坡發生頻率對地質、地形類因子的敏感性特征，而對各類因子之于滑坡分布的相對重要性缺乏興趣，更不會探究因子之間可能存在的交互影響；在因子選用方面也存在一定重復。此外，防汛部門所關注的溪洪—滑坡，主要指由于山區溪河洪水沖蝕坡岸所引發的各類急性或慢性滑坡，其收集的歷史災害信息可能與流域水文屬性具有更好的響應，這在以往研究中也較少關注。空間序列分析可以探究地形、土壤等短時間尺度內不隨時間變化因子對地理現象的影響，是因地制宜的開展目標對象關鍵控制因子研究的基本方法論。另一方面，地理探測器模型(Geodetector)是一種探測空間分異性并揭示其背后驅動力的一組統計學方法[13]，它不僅可以分析各自變量對因變量的影響程度，還可以探測自變量間交互作用于因變量的影響。據此，本文嘗試基于空間分析視角，以三峽庫區香溪河流域的溪洪—滑坡災害為研究對象，結合經典統計和Geodetector方法，系統分析研究對象的空間分異特征及其與主要自然因子的關系，探明以該流域為代表的庫區東部川鄂褶皺山地溪洪—滑坡災害空間分布的關鍵控制因子，對比分析地理探測器模型在本研究中所體現的優勢與局限。

1 研究區與數據方法

1.1 流域概況

香溪河位于長江上游三峽庫區左岸，發源于神農架南麓，由北向南流經湖北省興山縣和秭歸縣，于香溪鎮入注長江，干流全長約94 km，流域面積3 193 km2，是長江上游一級支流中距離三峽大壩壩首最近的中尺度河流。受亞熱帶大陸性季風氣候影響，本區四季分明，降水充沛，年均雨量1 000 mm以上，北部較南部稍高；流域地貌屬鄂西褶皺山地，為大巴山和巫山的余脈，地形以中低山地為主，山勢陡峻、河谷深切；流域植被曾于20世紀80年代遭到嚴重破壞，1989年起經過一系列水土保持項目治理后得到明顯改善，現森林覆蓋率達56%。流域涉及的主要構造形跡為近SN向展布的黃陵背斜和秭歸向斜，其中黃陵背斜一直處于上升趨勢，每年上升2～4 mm[14]；區內斷裂活動總體不顯著；出露巖層以碳酸鹽巖和砂、泥巖為主，主要包括三疊系中統嘉陵江組上部厚層白云制灰巖、嘉陵江組中部中—厚層灰巖、嘉陵江組下部薄層灰巖，侏羅系上統遂寧組磚紅色泥巖與砂巖互層，侏羅系上統沙溪廟組紅色泥巖和石英巖等。

1.2 數據來源

本研究著眼于目標流域溪洪—滑坡災害與主要自然因子的關系，因此在資料選取上力圖排除2002年之后水庫蓄水的影響，涉及的主要數據和來源詳見表1。

表1 研究區基礎數據與來源

1.3 指標體系建立與出圖計算

一般認為，溪河洪水的形成與強降水、大起伏地形、植被生長、土壤屬性、流域形狀等自然要素密切相關，但對于由溪洪引發的(及一般的)滑坡、泥石流，巖性、斷層(裂)、松散物質分布等地質要素的作用也非常重要。

基于學界對山洪、滑坡成因的基本認知[8-12,14-15]，結合資料收集情況，選取暴雨、土壤、巖性、斷層、歷史災害等指標作為分析基礎構建指標體系，各指標圖層如圖1所示。

圖1 香溪河流域主要自然因子的過程柵格圖

1.3.1 暴雨〔多年平均年暴雨極值(A1)與變異系數Cv(A2)〕 多年平均暴雨極值反映區域暴雨的絕對能量，暴雨極值的變異系數反映其不確定性，與溪洪、滑坡災害均存在理論上的聯系。統計收集的10個站點[16]1971—2002多年平均年暴雨極值量和Cv值，在GIS中進行反距離權重(IDW)插值，以結果的均方根誤差(RMSE)最小為宜，得到流域面上的A1和A2。

1.3.2 一般下墊面〔土地利用產匯流能力(A3)、土壤入滲性(A4)與NDVI(A5)〕 地表綜合產匯流能力不僅與溪洪、滑坡發育有直接關聯，也不同程度影響著地下水分布及孔隙水壓力情況。研究表明黏粒占比較少、有機質含量較高的土壤類型，或土地利用中林草覆蓋較多的類型，下墊面入滲性能較好[17-19]，不易產流形成洪水，但理論上有增加滑坡體重力勢能的可能，參考文獻[15]的工作，對相關指標進行單向打分；NDVI直接采用2000—2002年4—9月研究區NDVI值進行算術平均得到。

1.3.3 地形〔坡度(A6)與坡向(A7)〕 依據文獻[8-9]對庫區敏感性坡度的總結，應用GIS中的“slope”和“raster calculator”工具提取流域10°～45°之間的坡度，其余坡度賦予低值；坡向由GIS中“aspect”工具對DEM處理后，再對方向分級得到。

1.3.4 地質巖性〔脆弱性(A8)、巖性軟硬程度(A9)與斷層緩沖區(A10)〕 地層巖性是滑坡發育的重要考量因素，本研究結合巖性類別、表面破碎程度、產狀水平和完整性，使用專家打分法得到巖性脆弱性指標[20](表2)。此外，依據巖性軟硬分級，得到軟硬程度指標；汶川震區的研究表明[21]，大部分大型滑坡點分布在距離斷層5 km范圍內，庫區斷層活躍程度相對較弱，故緩沖區依距離斷層0.5,1,2,3,10,30 km設置。

表2 香溪河巖性脆弱性賦值

1.3.5 小流域水文屬性〔小流域面積(A11)、河網密度(A12)、主溝坡降(A13)和形狀系數(A14)〕 發育溪洪的小流域一般面積較小，基于DEM，利用GIS的水文分析工具，得到香溪河流域43個面積在6.5～183.6 km2的小流域，并利用GIS的地統計分析功能分別提取這些小流域大于2 km2的河網密度、主溝坡降和形狀系數。

1.3.6 歷史災害〔歷史溪洪—滑坡災害的綜合災度(B)〕 由于很多歷史災害點分布在小流域的邊緣地帶，不能充分反映流域特性的影響，且歷史災害點本身只能表現已經發生過的事件，對潛在的災害易發范圍沒有明確指示。基于“距離歷史災害點越近的地區孕災環境越相似”這一假設，使用空間插值的方法對收集的109個歷史災害點進行補點插值，具體方法是結合實地環境調查，對流域內及周邊無明顯災害記錄的地區補0點，對發災記錄較多地區的災害點合并為高值點，使網點在空間分布上相對均勻，最后使用IDW插值得到表征香溪河流域溪洪—滑坡災害潛在易發程度的綜合災度圖(圖2)。

1.4 數學方法

1.4.1 逐步回歸 逐步回歸方法各類文獻中均有詳細介紹，這里不再累述。本研究使用GIS提取的香溪河43個小流域的綜合災度(B)和除坡向(A7)外14項影響因子的平均值，以前者為因變量、后者為自變量構筑逐步回歸模型，分析主要自然因子對溪洪—滑坡災害的影響。

1.4.2 地理探測器模型(Geodetector model) 地理探測器模型是圍繞地理現象的分層(類)異質性或空間分異性研發的一種統計學方法[13]，它基于“兩變量空間分布越相似則關聯性越大”的假設，通過對因變量的再分類計算因變量類內方差與類間方差的關系，以此分析自變量對因變量的影響程度，其基本公式為：

(1)

除了上述基本的因子探測功能外，Geodetector還具有生態探測和交互作用探測功能，前者主要使用F檢驗判斷兩自變量對因變量的影響是否存在顯著差異；后者的原理為比較兩自變量在GIS中疊加后生成的q值與相應獨立自變量q值的關系，主要的交互作用方式有：

①q(A∩B)

② min(q(A),q(B))

③q(A∩B)>max〔q(A),q(B)〕 (雙因子增強)

④q(A∩B)=q(A)+q(B) (獨立)

⑤q(A∩B)>q(A)+q(B) (非線性增強)

由于本研究自變量擬分類數較多，僅提取43個小流域的數值建立關聯可能導致失真，因此這里除去反映小流域統計屬性的A11—A14，對余下的11項變量的柵格圖統一采樣精度(100 m)并轉化為點矢量圖，每張柵格圖得到3 123個矢量點，將自變量矢量點值分類后與因變量矢量點值聯立，采用因子和交互作用探測功能分析各自然因子對綜合災度的影響。

2 結果與討論

2.1 香溪河流域溪洪-滑坡綜合災度分布格局

依據上述方法，經自然斷點法分類后，得到香溪河流域溪洪—滑坡綜合災度顯示圖(圖2)，綜合災度是從宏觀格局上對區域潛在災害易發程度的綜合反映，并不是實際發災情況，因此不會具體到坡岸尺度。依據綜合災度分布，香溪河流域的溪洪—滑坡災害高易發區主要分布在興山以下的流域西南部，東部和西北也有少量分布，北部上游的大起伏山地反而發災較少，可能與該區人口密度不高、植被覆蓋較好有關。

圖2 香溪河流域溪洪-滑坡綜合災度分布

2.2 香溪河流域溪洪-滑坡綜合災度與主要自然因子的定量關系

通過逐步回歸和地理探測器模型分析，認為下墊面條件，特別是地層巖性與植被發育，是控制香溪河流域溪洪—滑坡綜合災度空間分異的主要方面，降水和土壤因素的影響有限，小流域水文屬性要素的作用不明顯。逐步回歸結果顯示進入模型的5項因子可以解釋因變量80.6%的變化，依據標準化回歸系數(表3)，各因子的貢獻率依次為巖性脆弱性(A8)28.6%，斷層緩沖區(A10)19.8%，NDVI(A5)18.9%，暴雨極值(A1)16.4%，以及巖性軟硬程度(A9)16.3%，巖層多傾斜、破碎、巖性較軟弱、距離斷層較近、植被覆蓋少的地區，發災潛力更大。

表3 香溪河流域溪洪-滑坡綜合災度與主要自然因子的逐步回歸結果

地理探測器模型分析的結果與回歸模型稍有不同(表4)，依據各因子對因變量的解釋度(q值)排序，僅土壤入滲性(A4)取代了暴雨極值(A1)。巖性脆弱性(A8)不僅在因子探測分析中具有最高的解釋度，與其他4項因子還分別存在非線性增強的交互影響。實際上與地質環境關系密切的巖性脆弱性(A8)、斷層緩沖區(A10)和巖性軟硬程度(A9)在兩類模型中均有穩定的表現，體現了地質要素對區域滑坡分布絕對的控制力。有學者[22]認為正是由于侏羅系黃陵背斜的褶皺隆升，引起秭歸向斜發生陸相繞曲沉積，導致沉積巖層產狀嚴重傾斜，加上沉積物質來源不同，沉積相軟硬互層，才有了當前香溪河中下游的滑坡集中區，本文的結果在一定程度上支持這一觀點。此外，NDVI(A5)在兩類模型的影響力僅次于地層巖性與斷層類因子，且與巖性脆弱性(A8)在交互作用下可以解釋因變量68.8%的變化，顯示了其對控制滑坡的積極影響。

表4 香溪河流域溪洪-滑坡綜合災度與主要

注：僅列出通過顯著性檢驗且q值大于0.1的變量。

一些被認為是滑坡發育基本因子的指標，如降水、地形類指標在兩類模型中的期望低于預期，這一方面是因為空間分析本身的“瓶頸效應”，即區內某類指標的數值趨于均化，或早已滿足因變量變化的要求，不構成發育瓶頸，則此類指標的作用就很難體現，如降水在香溪河乃至整個庫區均有非常充分的供給[16,23]，不是滑坡發育的瓶頸，因此本研究中相關指標的貢獻較低；地形與滑坡關系的特殊性也是原因之一，地形對滑坡發育的影響需與其他因子配合，并不是充分條件，因此Geodetector因子探測的q值較低(A6, 0.012;A7, 0.008)，但交互探測顯示，坡度與巖性脆弱性的交互作用強度達62.0%，表明坡度對滑坡的影響需聯合地質因素。考慮到此前也有案例反映坡度、坡向等指標對滑坡發育的影響并不大[24-25]，本研究認為單純的坡度、坡向對流域滑坡格局的影響有限，它們的作用更多體現在協同方面。

2.3 Geodetector的應用比較

Geodetector通過計算比較類別化了的自變量對因變量的類內與類間方差，得出因變量與自變量間的關聯大小或相似程度，以此判斷自變量對因變量的影響程度，為定量分析非線性變量間的關聯性提供了新的視角。本例中該模型探測了各單一因子對溪洪—滑坡綜合災度分布的影響，得到的各因子解釋度大小和排序與回歸結果相似，區別體現在土壤和暴雨因子的地位。直觀比較，土壤入滲性較之于暴雨極值，擁有與流域綜合災度更好的空間對應關系(圖1—2)。此外，即便使用Geodetector的矢量點數據進行回歸，土壤入滲性的貢獻仍然較低(表5)，Geodetector在處理全局性問題時對從屬因子擁有更高的敏感性和顯示度。然而，正是由于Geodetector對自變量的類別化處理，消除了自變量的數值屬性，使其無法探知自變量對因變量的影響方向，進而判斷這種影響是否存在假相關。此外，雖然“生態探測”功能可以對各自變量于因變量的影響格局進行判斷，但這種分析只能通過兩兩比較進行，且結果僅能說明比較因子在影響格局上的相似性，而影響格局相似的因子并不一定是重復的。回歸模型可以近似表達溪洪—滑坡綜合災度與自然因子的定量關系和作用方向，偏相關或逐步回歸等方法也能從統計上剔除可能存在的假相關或重復性變量，與Geodetector具有一定的互補性。

表5 香溪河流域溪洪-滑坡綜合災度與主要自然因子的矢量點數據逐步回歸結果

3 結 論

結合逐步回歸和Geodetector方法，從空間分異視角對三峽庫區香溪河流域的溪洪—滑坡災害與其自然影響因子的關系進行了定量研究。結果表明，進入逐步回歸模型的5項因子總解釋度80.6%，貢獻率依次為：巖性脆弱性28.6%，斷層緩沖區19.8%，NDVI 18.9%，暴雨極值16.4%，巖性軟硬程度16.3%；地理探測中對因變量解釋度較大的5項因子分別為：性脆弱性24.5%，斷層緩沖區20.9%，NDVI 18.3%，土壤入滲性16.3%，巖性軟硬程度14.8%；巖性脆弱性與NDVI，坡度的交互作用強度分別達68.8%和62.0%；基于定量分析、定性理解的原則，認為地質要素該本地區溪洪—滑坡綜合災度的空間分布處于支配地位，地層巖性越傾斜、破碎、軟弱，距離斷層越近，發災可能性越大；植被生長對綜合災度分布也有重要影響，降水、土壤屬性的控制作用較小，小流域水文屬性與綜合災度的關系不明顯，地形的影響體現在協同方面。

需要強調的是，本文的結果僅明確了研究區內影響溪洪—滑坡災害分布的瓶頸性因子，未進入模型或貢獻低于預期的因子，未否認它們對目標現象的機理性影響，僅表明這些因子并非這一地區溪洪—滑坡災害分布的瓶頸。Geodetector為定量分析非線性變量間的關聯性提供了新的視角，但它消除了自變量的數值屬性，無法探知自變量對因變量的影響方向，多變量分析時也很難判斷各自變量對因變量可能存在的重復性問題，因此建議結合經典統計方法，對結果進行綜合分析與判斷。

致謝：感謝中科院水利部成都山地災害與環境研究所的張信寶研究員在巖性脆弱性賦值方面給予的支持與幫助；感謝長江科學院重慶分院的盧陽高級工程師和秦維工程師提供的歷史災害數據支持。