基于流域单元的岷江上游泥石流危险性评价

李 炫，杨本勇，范建容，张建强，李磊磊

（1.中国科学院 水利部 成都山地灾害与环境研究所，成都610041；2.中国科学院大学，北京100049；3.国家测绘地理信息局 第三地理信息制图院，成都610100）

泥石流是在暴雨、融雪等外界条件的激发下产生的，来势凶猛并且过程短暂，具有运动速度快、暴发突然、能量巨大等特征，其携带的大石块具有冲击力强和破坏性大的特点［1］，泥石流通常给山区人民带来较大的经济财产损失和人员伤亡，如2010年8月7日甘肃省舟曲县城后山三眼沟及罗家沟突然暴发大规模泥石流，造成了1 744人死亡和失踪［2］。因此，做好泥石流的评估、预警工作，对防范泥石流灾害的发生，减少灾害带来的损失，以及区域的规划建设都具有重要的意义。

泥石流危险性评价是灾情评估、预测、防灾救灾的基础［3］。在国内好多区域泥石流危险性评价工作中大都采用栅格单元［4－6］，其结构简单计算高效，被国内学者广泛的应用于泥石流危险性评价中。但这种评价单元与地质地貌以及泥石流发生的其他环境条件缺乏联系，无法体现泥石流发生的实际情况。而泥石流沟谷的地形地貌特征、丰富的物质条件以及降雨是泥石流发生的必要条件［7］，基于流域单元的泥石流评价能充分考虑到与泥石流形成相关的内在因子的影响。目前流域划分方法大都是基于DEM数据［8－9］，提取流域时阈值的确定大都通过人工尝试，没有一个较为精确合理的方法，自动化程度不高，会产生的大量非闭合子流域，划分出的流域与实际的泥石流流域有所差别。本文选取岷江上游区域作为研究区，采用均值变点法和河网密度法探究研究区流域划分的合理阈值，并以划分的流域单元为基础，从泥石流的形成条件出发，选取地层岩性、断裂密度、平均坡度、形状系数、沟道比降、24h最大降雨量等六个指标，对岷江上游地区泥石流的危险性进行综合的评估。

1 研究区概况

岷江上游地区辖四川省阿坝藏族羌族自治州的汶川、茂县、理县、黑水和松潘五县的全部和都江堰小部分区域，流域面积为23 037km2，流域内地势起伏较大，最高海拔6 246m，最低海拔455m，该区地貌类型以山原和高山峡谷为主，山地灾害严重，泥石流灾害处于高风险区，是四川省乃至是全国泥石流的重灾区［10］，地质结构复杂，断裂发育，龙门山构造体系、岷江纵向构造体系共同作用，茂汶断裂、北川—映秀断裂、岷江断裂等深大断裂也分布于此，加上松潘地震带和龙门山地震带等活动地震带的影响，历史上发生过许多强大的地震如1976年的松潘地震和2008年的汶川地震等。受地震等灾害的影响，研究区内地质、地形条件都发生了不同程度的改变，地表植被损毁，松散固体物质增多，使得区域内地表抗蚀能力减弱，易损性变大，这些因素都为泥石流的形成创造了良好的条件。同时，该区在南北走向的特殊地貌和西南季风的共同作用下，焚风效应显著［11］，降雨季节性明显，5—8月集中了全年80%～90%的雨量，汛期时常会有暴雨出现，山洪和泥石流灾害频发。对岷江上游地区进行泥石流危险性评价，能为区域内泥石流灾害的预防提供参考，对减少人民经济财产损失和区域的长远发展都具有非常重要的意义。

2 流域单元划分

数字高程模型（Digital Elevation Model，DEM）是反映流域地形特征的数字模型，基于DEM的流域划分，大大节省了的人力、物力，从提取的效率和精度来看都是切实可行的［12］。基于DEM的流域划分需要对数据进行填洼处理，用网格坡度确定每个单元格的水流方向，得到对应的流向矩阵，计算汇流累积量，再给定一个汇水阈值，凡是汇水大于该阈值的网格，均为河网内的网格，将这些网格连接后，得出整个流域的排水网络。当流域的排水网络生成以后，就可以确定整个流域的界限。因此，阈值的确定是河网提取和流域划分的关键。

2.1 阈值确定和流域划分

阈值的确定方法有多种，主要有水系分形分维法［13］、河网密度法［14］、均值变点法［15］等等。本文以空间分辨率为25m×25m的DEM数据作为数据源，综合运用均值变点法和河网密度法推求阈值，避免阈值选取的随意性。河网密度法认为阈值与河网密度（或河源密度）关系曲线趋于平缓时所对应的点为合理的阈值。如图1所示，为河网密度与阈值之间的变化关系。由图1可以看出：随着阈值的不断增大，总河网长度会逐渐减少，河网密度也逐渐减少；起初河网密度大小减少比较剧烈，随着阈值的慢慢增大，河网密度的减少变得越来越缓，当阈值为4 000的时候，二阶导数趋近于0，河网密度变化趋于稳定。

图1 河网密度随阈值变化趋势

均值变点法的离散数据模型为［15］：

a1＝…am1－1＝b1，am1＝…am2－1＝b2，amq＝…＝aN＝bq＋1，此处1＜m1＜m2＜…＜mq≤N，如果bj＋1≠bj，则mj就是一个变点。

计算公式为：

（1）令i＝2，…，N对每个i将样本分为2段：x1，x2，…，xi－1和xi，xi＋1，…，xN。

式中：Si——各段样本的方差和；Xi1，Xi2——每段样本的算术平均值。

（2）计算统计量。

式中：xm——样本均值；s——总体的方差；xt——单个样本值。

选取平均河网密度作为样本序列X，按照公式（3）计算样本总体平均值和方差由公式（2）计算出分段样本的算术平均值与方差，从而得到S与Si差值变化曲线如图2所示。变点的存在会使S和Si的差值增大，这个方法对恰有一个变点的检验最为有效［11］。由图2可以看出第5个点对应的差值最大，说明第5个点就是我们要找的变点，第5个点对应的阈值为4 000，这与之前用河网密度法所求的阈值相同。

图2 s与si差值变化曲线

根据计算的阈值提取出河网，利用GIS工具Hydrology Tools划分出子流域单元。将划分出的流域单元与实际的流域单元进行比较发现，大部分提取的流域单元与实际的流域单元都是比较接近的，但提取出的流域中包含大量的非闭合子流域，这些非闭合的子流域会对泥石流的危险性评价造成了干扰。同时，流域单元划分时只基于地形数据，提取出的流域单元在某些地势平坦的区域，流域边界并不准确，会出现整个流域被刨分成多个流域，某些细小的流域被划分成一个流域等情况，所以，应根据实际地形对流域单元进行修正。

2.2 流域修正

针对上述提取的流域边界与实际流域单元不吻合的情况，本文运用DEM生成的山体阴影图像结合高分辨率遥感影像，对部分集水单元进行手动修改，去除某些非闭合子流域，按照流域的实际形态合并或修正子流域边界，最终得到岷江上游地区的小流域分布（图3）。

图3 岷江上游地区流域单元划分

3 流域单元泥石流危险性评价

3.1 评价因子的选取与分级

影响泥石流形成的因素较多，根据形成机制来看，影响因素主要可以分为：地形地貌条件、物源条件、降雨等诱发条件。地形地貌条件制约着泥石流的形成、运动和规模等特征，地形因素包括流域面积、纵比降、坡度、形状因子、发育程度等［16］。韩用顺等人就选取了岩性、平均坡度、流域面积、相对高差、主沟长度等因子实现了汶川震后灾区泥石流危险度评价［17］。谭万沛等在研究区域暴雨泥石流预测时表明形状系数和沟道比降对泥石流的发生具有很大的制约作用，是影响泥石流流速和输送能力大小的重要条件，是评价泥石流流域危险度的重要因子［18］。孟国才等人在研究岷江上游泥石流特征时指出岷江上游地区地质构造复杂，断裂发育，新构造运动隆升强烈，茂汶断裂、北川—映秀断裂、岷江断裂等深大断裂分布于此，地层岩性和地质构造是岷江上游地区泥石流危险性评价的重要因素［11］。降雨是诱发泥石流灾害最直接的因素，对泥石流的形成起关键作用的是1h雨强和10min雨强，但是小时雨强和10min雨强不易观测和记录，而日降雨量的相对容易些，且三者之间存在正相关关系，一定程度上可以用日最大降水量来同时表示雨量和雨强。

本文在上述分析和前人研究的基础上，从影响泥石流的三大因素出发，筛选出对岷江上游地区泥石流发生起一定主导作用的地层岩性、断裂密度、平均坡度、形状系数、沟道比降、24h最大降雨量等六个指标建立泥石流危险性的评价模型。通过对研究区978个流域的6个指标进行量化，并进行相关分析，以避免影响因子的重复引入，相关矩阵如表1所示。从相关矩阵可以看出各个因子之间没有明显的相关关系，均从不同角度反映了流域的特征。

表1 评价因子之间的相关矩阵

由于各个评价因子的计量单位不同，取值范围变幅很大，会影响评价的结果，所以需要对上述指标进行归一化处理。本文对各因子采用5级量化指标，通过不同级别反映各因子对泥石流发育影响程度的差异，等级越高，说明其对泥石流的影响程度越大，泥石流危险性越高。对于地层岩性因子，由获得的研究区地质岩性图，综合考虑本区岩土体类型、物理力学性质，并结合岩土体结构特征，参考工程岩体分级标准，将研究区内岩组分为坚硬岩组（花岗岩、正长岩、石英等）、较坚硬岩组（白云岩、大理岩等）、软硬相间岩组（千枚岩、砂质泥岩等）、软弱岩组（含煤层、半固结岩）、极软岩组（泥岩和煤系地层）五类，危险度分别赋予1，2，3，4，5。对于断裂密度、平均坡度、形状系数、沟道比降、24h最大降雨量等评价因子，采用自然间断点分级法，划分为五个等级，并按照各指标值的大小分别赋值为5，4，3，2，1。具体归一化方案见表2。

表2 评价因子等级划分

3.2 评价因子权重的确定

在评价指标体系中，不同影响因子对泥石流灾害的贡献程度是不同的［6］，这种贡献程度的不同通过各评价因子的权重系数加以体现。目前泥石流评价中权重确定的方法一般有两种，主观赋权法，和客观赋权法。主观赋权法如层次分析法、专家打分法等，权重的结果受到人为因素影响较大。而客观赋权法如主成分分析法和因子分析法等是根据各指标间的相关关系或各项指标值的变异程度来确定权重系数，可以避免人为因素的影响。本文采用主成分分析法来确定个因子的权重，计算方法如下［19］：

（1）对各因子进行相关分析，确定相关系数矩阵R（又称协方差矩阵）；

（2）计算矩阵R的特征值和特征向量，即通过正交变换将矩阵化为对角矩阵，对矩阵中的对角元素即为所求特征值，按其大小排列，分别称为第一、二特征值（λ1＞λ2＞…λn）；

（3）计算主成分的累积方差贡献率，并确定出累积贡献率大于85%的主成分数；

（4）计算主成分中载荷系数，即特征值的方根与相应的特征向量的乘积，并将其归一化，得到各参评因子的权重。

将量化的978个流域数据作主成分分析，以累积方差大于93%确定主成分数（4个），确定出各个评价因子的权重值，见表3。

表3 评价因子权重值

3.3 流域单元的危险度评价

4 结果与验证

通过加入研究区内已查明的205条泥石流沟，与评价结果相叠加，结果如图4所示。研究区内有9条泥石流沟位于极低危险区内，占总数的4.39%，21条泥石流沟位于低度危险区域内，占总数的10.24%，52条泥石流沟位于中度危险区域内，占总数的25.36%，84条泥石流沟位于高度危险区内，占总数的40.98%，39条沟位于极度危险区域内，占总数的19.03%。已查明的泥石流沟一共有175条位于中度危险、高度危险、极度危险区域内，占总数的85.37%。结果分析表明，评价的结果能较好的反映出岷江上游的泥石流分布情况。

图4 岷江上游泥石流危险分布

5 结 论

本文以DEM数据为基础进行流域划分，综合运用均值变点法和河网密度法确定出岷江上游地区流域划分的合理阈值为4 000，将自动划分的流域单元，参考实际的遥感影像，对子流域单元进行合并或修正，最终得到岷江上游地区大小流域共978个。通过遥感影像的修正，划分的流域与实际的流域比较吻合。

从影响泥石流形成的三大条件出发，选取地层岩性、断裂密度、平均坡度、形状系数、沟道比降、24h最大降雨量等六个评价指标，运用主成分分析法确定各影响因子的权重值，对流域的量化数据进行归一化、分级处理，在ARCGIS支持下完成了岷江上游地区泥石流危险性评价，得到危险性评价图，其中极度危险度流域97个，占总流域的9.92%，高度危险流域205个，占总流域的20.96%，中度危险流域286个，占总流域的29.24%，低度危险流域240个，占总流域的24.54%，极低危险流域150个，占总流域的15.34%。用已有灾害数据验证表明，评价结果能较好的反映出岷江上游地区泥石流的分布规律，为该区泥石流灾害的防治提供了参考和依据。

［1］ 张怀珍，范建容，郭芬芬，等.中国单沟泥石流危险度评价模型比较研究［J］.水土保持研究，2011，18（4）：20－26.

［2］ 余斌，杨永红，苏永超，等.甘肃省舟曲8.7特大泥石流调查研究［J］.工程地质学报，2010，18（4）：437－444.

［3］ 孙绍骋.灾害评估研究内容与方法探讨［J］.地理科学进展，2001，20（2）：122－130.

［4］ 唐川，朱大奎.基于GIS技术的泥石流风险评价研究［J］.地理科学，2002，22（3）：300－304.

［5］ 宁娜，马金珠，张鹏，等.基于GIS和信息量法的甘肃南部白龙江流域泥石流灾害危险性评价［J］.资源科学，2013，35（4）：892－899.

［6］ 王欢，丁明涛，陈廷方.基于GIS的三江并流区泥石流危险性评价［J］.水土保持通报，2011，31（5）：167－170.

［7］ 王晓朋，潘懋，徐岳仁.基于流域单元的泥石流区域危险性评价［J］.山地学报，2006，24（2）：177－180.

［8］ 李丽.舟曲县泥石流危险性评价［D］.北京：中国地质大学，2012.

［9］ 王琼.基于流域尺度的震后汶川县潜在泥石流危险性评价［D］.成都：成都理工大学，2012.

［10］ 刘希林，苏鹏程.四川省泥石流风险评价［J］.灾害学，2004，19（2）：23－28.

［11］ 孟国才，王士革，谢洪，等.岷江上游泥石流灾害特征分析［J］.灾害学，2005，20（3）：94－98.

［12］ 陈加兵，励惠国，郑达贤，等.基于DEM的福建省小流域划分研究［J］.地球信息科学，2007，9（2）：74－77.

［13］ 杨邦，任立良.集水面积阈值确定方法的比较研究［J］.水电能源科学，2009，27（5）：11－14.

［14］ 孔凡哲，李莉莉.利用DEM提取河网时集水面积阈值的确定［J］.水电能源科学，2006，23（4）：65－67.

［15］ 赖晗，芮小平，梁汉东，等.基于均值变点分析的三峡库区河网提取研究［J］.测绘科学，2012，5（37）：056.

［16］ Liu C N，Dong J J，Peng Y F，et al.Effects of strong ground motion on the susceptibility of gully type debris flows［J］.Engineering Geology，2009，104（3）：241－253.

［17］ 韩用顺，李龙伟，朱颖彦，等.汶川地震灾区泥石流危险性评估：以都江堰—汶川公路为例［J］.中国水土保持科学，2012，10（1）：6－11.

［18］ 谭万沛，王成华，姚令侃，等.暴雨泥石流滑波的区域预测与预报：以攀西地区为例［M］.成都：四川科学技术出版社，1994.

［19］ 韩小孩，张耀辉，孙福军，等.基于主成分分析的指标权重确定方法［J］.四川兵工学报，2012，33（10）：124－126.