川鄂褶皱山地溪洪-滑坡灾害与主要自然因子的关系——以香溪河流域为例

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

(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为定量分析非线性变量间的关联性提供了新的视角，但它消除了自变量的数值属性，无法探知自变量对因变量的影响方向，多变量分析时也很难判断各自变量对因变量可能存在的重复性问题，因此建议结合经典统计方法，对结果进行综合分析与判断。

致谢：感谢中科院水利部成都山地灾害与环境研究所的张信宝研究员在岩性脆弱性赋值方面给予的支持与帮助；感谢长江科学院重庆分院的卢阳高级工程师和秦维工程师提供的历史灾害数据支持。