0 引言

泥石流灾害一直是威胁山区人民生产、生活安全的主要灾害之一。它是一种由降雨引发的、在沟谷限定的沟槽内及一定的坡度条件下, 由固液两相混合形成的快速由上游流向下游方向的重力流体, 呈现出黏性层流或稀性紊流的运动状态; 具有流速快、流量大、来势凶猛并且破坏力强、发生时间短及复发频率高的特点^[1-2]。

泥石流灾害广泛分布于世界各地, 据现有资料, 除南极洲外, 世界各大洲均有泥石流灾害爆发的痕迹, 全球超过60个国家饱受泥石流灾害的威胁^[1]。在亚洲, 中国、印度、日本、菲律宾、马来西亚等国都是泥石流灾害重灾区, 每年都因泥石流造成了巨大的财产损失和众多的人员伤亡。我国是一个多山的国家, 山区占全国陆地总面积的2/3^[3]; 亚热带气候的影响结合我国呈阶梯状展布的地貌轮廓, 造就了泥石流易发的客观条件, 使我国成为世界上泥石流灾害重灾区之一。据康志成等^[4]2004年统计:近年来全国29个省771个地区遭受过泥石流的危害, 平均每年每县的泥石流爆发次数为18次; 由泥石流灾害造成的人员死亡人数在近40年来高达14万人, 经济损失上万亿元^[5]。可见, 泥石流灾害对于我国易发灾害的山区人民来说是一个巨大的威胁, 泥石流灾害评价研究工作对当地的灾害预警及防治工作来说至关重要。

1 泥石流危险性评价方法

如今, 泥石流的风险评价分类方法逐渐由早期的定性分析转变为半定量评价。国外许多非线性研究方法的引进, 结合泥石流的风险评价手段, 亦使泥石流风险评价增加了许多新的方法。

1981年, 美国Hollingsworth等^[6]采用打分方法构建了泥石流危险性评价构架; 而因子叠加法的出现为以后的泥石流危险性评价指明了方向。现今的评价方法大多以此为基础, 从影响泥石流的启动、运移、堆积等因素入手, 采用不同数学方法确定因子权重及评价构架, 从而计算出泥石流危险性指标。而随着计算机技术的发展以及地理信息系统GIS平台的出现, 泥石流评价也会向着精度更高、速度更快、集成性能更高、更加便捷的方向发展^[7]。

根据研究范围, 泥石流的危险性评价可分为单沟危险性评价和区域危险性评价。其中:单沟危险性评价对待评价的泥石流沟认识明确、指标要求精度高、定量化程度高^[8]; 区域泥石流危险性评价更注重区域内各单元之间危险性大小的对比, 表示区域内相对各单元的危险性, 评价因子也有许多为相对评价指标。目前较为常见的泥石流危险性评价方法如图 1所示, 广泛应用的评价方法有可拓学、灰色关联度分析、信息熵、分形理论、模糊数学、神经网络等^[9-14]。这些方法中有些用于评价因子的选取及权重值计算, 有些在此基础上用于泥石流的危险性评价分类。

回归分析评价方法是建立在大量观测样本的前提下, 根据已知样本危险度进行回归分析, 推测未知样本危险度, 对泥石流沟谷危险性进行评价。模糊数学评价方法^[15]是先根据人为确定的各指标分类区间得出各指标模糊分级隶属度函数; 然后将需要评价的样本指标代入隶属度函数, 得出隶属度矩阵; 再由权重向量与隶属度矩阵运算后得出最终隶属度并得出评价结果。神经网络评价方法^[10]是根据观测样本选取神经网络结构, 训练得出评价模型。灰色关联度分析评价方法^[16]首先求得各因子与比较数列的关联系数; 然后通过关联系数与权重得出关联度, 关联度最大的类别即为评价结果。信息熵理论由香农提出, 我国学者艾南山将其应用于泥石流危险性评价。信息熵评价方法^{[13, 17]}是先通过地理系统信息熵与Strahler面积-高程积分值的关系得出评价流域地貌信息熵, 再根据艾南山提出的流域地貌信息熵与泥石流活动性的关系对泥石流危险性进行评价。地貌信息熵能够反映流域内地貌的变化能力, 在一定程度上评价泥石流危险性。分形理论评价方法^{[14, 18]}以盒维数为目标, 获取不断变换尺度下的研究目标数; 然后求得双对数坐标线性回归斜率, 即得到目标的盒维数, 进而研究泥石流危险性与盒维数相关性。

从评价方法发展情况来看:泥石流评价常用方法第一种是以可拓学为代表的, 即确定评价因子权重, 各因子打分并代入评价模型, 得出最终评价结果; 另一种是以神经网络为代表的, 即经过对已知评价结果的样本训练, 得出训练后模型, 再将各评价因子指标代入已确定模型得出评价结果; 还有一种是以聚类方法为代表的, 即根据评价因子自身特点评价区域内单沟相对危险度的评价方法。

泥石流评价分类方法与选取的评价因子密切相关, 在选取多个评价因子后, 综合考虑各项因子的权重, 计算各因子相对各分区的隶属度, 最后得到泥石流的评价分类; 然而, 一般的评价方法需要对选取的各评价因子的值域分区, 在未知单个影响因子的危险分区情况下, 评价方法的准确性难免不尽人意。

模糊聚类方法能够避开人为制定各评价因子的危险分区过程, 以评价因子的自身特性为基础, 考虑各沟的相似性, 得出各流域的泥石流危险性分类。聚类方法属于无监督机器学习理论, 也不需要训练样本; 通过得出模型进而评价泥石流分类的过程, 更适合样本缺乏情况下对泥石流沟谷的危险分类。支持向量机(support vector mechine, SVM)训练得出的模型只需要支持向量样本, 可以解决小样本、非线性及高维模式识别的问题, 适合小样本条件下的泥石流危险性分类, 同时能够明确划分聚类边界。因此, 本文结合模糊c均值聚类(FCM)方法与支持向量机(SVM)理论, 对研究区域的泥石流沟谷进行了危险性分类。

2 模糊c均值聚类方法

模糊c均值聚类(FCM)方法是由Bezdek等^[19]于1981年提出的, 作为早期硬c均值聚类(HCM)方法的一种改进, 是一种广泛应用的聚类方法。它能够根据样本数据自身特点进行分类, 无需人为给定各因子分类范围。它以目标函数的最小化为目标, 通过目标函数的变化小于阈值或者样本分类无变化为结束条件, 经过多次迭代, 得出样本数据的隶属度矩阵。根据隶属度矩阵中各样本(泥石流沟)对于所分类别的隶属情况, 得到泥石流危险性评价结果, 样本对应危险性类别隶属度最大的类别即为该样本所属类别。这种方法能够处理很难分类的样本并得出较为满意的结果^[20]。

FCM的计算原理^[21]如下:

设有样本集X=(x₁, x₂, …, x_k, …, x_n)', 其中x_k代表第k个泥石流沟的评价因子集, 隶属度矩阵

其中:u_ij为第j个样本属于第i类的隶属度; c为聚类的类别个数, 需要人为指定。

1)初始化隶属度矩阵U, u_ij为0到1的随机数并且满足条件:

2)计算聚类中心c_i:

式中:x_j为第j个泥石流沟的评价因子集; m为加权指数, m∈[1, ∞)。有研究表明:m越大, 聚类效果越差; m越小、越接近于1, 聚类越接近于硬聚类算法。因此, 参考前人文献[20-21], m一般取值为2。

3)计算目标函数J:

式中, d_ij为c_i与x_j之间的欧几里德距离, d_ij=‖c_i-x_j‖。

4)计算新的隶属度矩阵U。返回步骤2), 隶属度更新采用下式计算:

5)如果目标函数J两次变化值小于某个阈值ε, 或者分类结果不再变化, 结束循环, 得出最终隶属度矩阵。隶属度矩阵U内元素数值最大对应的类别即为所属类别。

3 支持向量机法

Vapnik提出的支持向量机是20世纪90年代发展起来的机器学习方法。支持向量机以训练误差作为优化问题的约束条件, 以置信范围值最小化作为优化目标, 即SVM是一种基于结构风险最小化准则的算法。它是基于统计学习理论产生的学习方法, 避免了传统方法容易发生的过学习问题, 适用于小样本及高维次分类情况, 其推广能力明显优于一些传统的学习方法^[22]。具体的计算方法如下:

设训练样本X={F_i, G_i}(i=1, 2, …, l), F_i∈R^D, G_i∈{-1, +1}。其中:F_i为各沟泥石流样本; G_i为样本标签; l为样本数; D为评价因子数。

线性判别函数的一般形式为

式中:ω、ω₀为待定变量, 可通过训练样本计算得出。sgn[g(F)]可判别样本类别, 通过归一化, 使样本满足|g(F)|=|ω·F+ω₀|≥1, 当取极值情况时, |g(F)|=1, 此时函数为类边界, 即有

g₁(F)与g₂(F)为两个类的边界, 此时, 两个类边界之间的分类间隔为。为使分类间隔最大, 将问题转变为求得最小‖ω‖, 进而为求解方便转为求解最小, 则分类函数的求解转为

这是一个凸二次规划问题, 根据最优化理论, 这一问题存在唯一全局最小解。应用lagrange乘子法:

式中, α_i为拉格朗日乘子。

式中, x为向量。根据约束条件式(9)有:。因此, 原目标函数可转化为

根据目标函数式(11), 得出其对偶问题为, 当满足KKT条件(karush-kuhn-tucker最优化条件)时, 对偶问题与原问题有相同的极值d=p; 因此求解原问题转化为了求解其对偶问题。KKT条件为

原问题在满足KKT条件下变为

为了求得, 将式(9)分别对ω, ω₀偏微分, 可得

将式(14)代入式(9), 综合KKT条件整理后得到凸二次规划问题的对偶问题:

采用SMO优化算法求解上述对偶问题, SMO选用两个自由维度并控制其他变量, 通过不断反复迭代求得函数, 最终得出判别函数:

以上解决了线性可分条件下的判别函数解。对于线性不可分的情况, 如果提高样本的维次, 可将样本在高维空间分类; 通过引入核函数, 能将样本投射到高维空间, 能将原本线性可分条件下的解推广至线性不可分条件, 判别函数为

式中, K(F_i·x)即为核函数。常用的核函数有线性核函数、多项式核函数、高斯核函数及sigmoid函数, 其中, 高斯核函数K(x, z)=exp 最为常用, 本文使用核函数即为高斯核函数。

4 FCM-SVM联合方法及应用

本文使用FCM与SVM两种方法联合对泥石流危险性进行了评价。以各泥石流沟的9项评价因子为基础, 采用FCM对泥石流沟进行聚类分析; 找出聚类中心; 而后以聚类中心为训练样本, 进行了SVM模型训练, 得到了泥石流沟的危险性评价模型; 并对样本泥石流沟谷进行了评价。

4.1 工程概况

房山区隶属北京市, 位于北京西南, 地势西北高而东南低, 地貌复杂多样。地处暖温带半湿润季风大陆性气候区, 境内地貌复杂, 山区与平原间相对高差悬殊, 气候有明显差异。四季特征鲜明。岩性以凝灰质板岩、砂岩、砾岩、火山岩等为主。南窖沟位于北京市房山区西北部山区, 总面积24 km², 内有水峪村、北安村、南窖村、三合村等村落, 居住人口达七千人左右。若爆发泥石流将直接对沟内人员安全和财产造成威胁。本文选取南窖沟流域内7条支沟作为研究目标, 对各支沟进行泥石流危险性分类。

4.2 评价因子选取

影响泥石流危险性的主要有地质、地貌、气象、水文、植被及人类活动等情况, 这些因素均能影响流域内泥石流危险性的评判。泥石流研究区域内一次泥石流最大冲出量L₁、泥石流的发生频率L₂明显地影响流域内泥石流危险性, 可作为重要的影响因子; 此外综合各方面因素, 参考前人常用的泥石流影响因子选取方法^[23], 选取了包括流域面积S₁、主沟床长度S₂、流域最大相对高差S₃、流域切割密度S₄、主沟床弯曲系数S₅、泥沙补给段长度比S₆、流域内人口密度S₇为影响因子, 降水因素为重要的泥石流启动条件。因本文研究区域内降雨条件基本相同, 故没有将其列入影响因子。

4.3 FCM计算

以研究区域内7条泥石流沟谷为聚类样本, 各泥石流沟评价因子及参数见表 1, 根据上述FCM方法对其进行聚类。首先对研究流域内各项评价因子进行标准化, 随机选取初始聚类中心; 其次对每个循环后的目标函数值进行计算, 直到达到循环结束条件; 最后得出各聚类中心及样本隶属度矩阵, 将样本分为3类。计算得出各沟隶属度矩阵如表 2所示。

根据隶属度矩阵(表 2)分类结果得出:属于第Ⅰ类的有水峪沟; 属于第Ⅱ类的有南岗沟、猫儿山沟; 属于第Ⅲ类的有三合村沟、左支沟1、左支沟2、北安村沟。

4.4 FCM-SVM联合法计算及结果

使用FCM得出的聚类中心为训练样本, 经过SVM模型训练, 得出SVM分类模型, 最后代入各泥石流沟数据, 得出最终的泥石流危险性分类结果。具体计算过程如下。

1)首先采用FCM方法求出样本集的聚类中心。

2)将聚类中心集作为训练样本, 取其中一个聚类中心为第一类x_k={F_k, 1}, 其余中心为另一类x_i={F_i-1}(i≠k, i=1, 2, …, c)。

3)按照SVM训练方法得出模型, 代入泥石流沟谷样本, 得出与k聚类中心同类样本并记录。

4)变换其他聚类中心为x_k={F_k, 1}并转入步骤2), 直至所有聚类中心都都经过轮换, 最后得出所有与各聚类中心同类的样本, 赋予原本聚类中心的类别标签, 即得出最后分类结果。

对比可拓学、灰色关联度、FCM及FCM-SVM分类方式得出的结果见表 3。由表 3可知:FCM算法结果中的Ⅰ、Ⅱ类相当于可拓学轻度危险指标, Ⅲ类为轻度至中度危险, 两种判别方式在趋势上相一致, 均能够较好地反映沟谷的泥石流危险性情况; 大部分沟谷的灰色关联度危险性判别情况与可拓学一致, 与FCM算法的结果也有一定的相关性, 但是由于灰色关联度分类类别只有轻度及中度两种, 难以明确分辨出变化情况; FCM-SVM危险性分类与FCM所得分类情况基本一致, 但是相比FCM结果, 前者经过SVM的运算, 得出了该地区的SVM泥石流危险性判别模型, 与FCM算法的概念不同, 确定了泥石流判别函数条件, 为该地区其他沟谷泥石流判别提供了一种方法。

FCM及FCM-SVM分类方法得出的结果充分反映了泥石流沟谷的危险性特征, 与实地调查结果相符。第Ⅰ类泥石流沟谷, 泥沙补给长度比例较低, 最大相对高差小, 综合活动强度弱; 泥石流发展阶段处于老年期, 水域沟地势相对平缓, 侵蚀作用不明显, 未见明显的大规模近期泥石流堆积物, 泥石流危险性相对较小。第Ⅲ类泥石流沟谷, 泥沙补给长度比例高, 相对高差较大, 平均坡度陡, 综合活动强度较强, 容易发生泥石流灾害; 现场可见泥石流沟谷内大量明显的崩坡积物及泥石流堆积体, 植被覆盖较差, 表明该沟谷的泥石流危险性高。第Ⅱ类泥石流沟谷爆发规模中等, 泥石流灾害危险性处于第Ⅰ类与第Ⅲ类之间; 沟谷内植被覆盖茂盛, 有少量的泥石流活动迹象。

为了进一步说明FCM-SVM方法的科学合理性, 本文另外选取金沙江研究区典型的20条泥石流沟谷作为研究对象^[24], 选取流域面积、主沟床长度、流域最大相对高差、植被覆盖率、24 h最大降雨强度、流域切割密度、主沟床弯曲系数、泥沙补给段长度比等12项影响因子为评价指标, 进行了FCM-SVM方法分类, 并与可拓学方法进行了对比(表 4)。计算结果表明:FCM-SVM算法所得的Ⅰ类泥石流沟一般为可拓学轻度危险, 少量出现中度危险; Ⅱ类泥石流沟大多处于中度至轻度危险; Ⅲ类泥石流沟仅有一条, 为重度危险性。评价结果基本相近, 反映了由评价因子客观因素确定的泥石流危险性, 证明了采用FCM-SVM方法对泥石流危险性进行评价具有可行性。

5 结语

选用模糊c均值聚类(FCM)理论及支持向量机(SVM)理论对泥石流危险性进行了分类, 得出以下结论:

1)现有泥石流评价方法按研究范围分为以单沟评价指标为基础的单沟泥石流评价方法和以一定区域内的区域性指标为评价标准的区域泥石流危险性评价方法, 各种不同的评价方法根据评价因子及自身算法有其自身的适应性。

2)应用FCM方法对流域泥石流危险性分类, 能够依据研究流域内的泥石流影响因子特征, 根据样本自身特点, 对研究区内泥石流危险性进行合理的分类, 克服了主观因素可能产生的不准确性。分类结果以隶属度矩阵的形式展示, 量化了研究目标对各分类级别的隶属情况, 提供了更为详细的信息。

3)SVM具有适合小样本、高维次、非线性分类模式的特点, 能够针对泥石流分类样本较少, 危险性影响因子多及分类非线性的特征, 适合泥石流危险性的分类, 基于统计学习理论的特点, 能够使得出的模型更具推广能力。

4)应用FCM理论和SVM理论时, 首先利用FCM算法得出各聚类中心, 然后使用SVM对聚类中心进行训练, 最后得出危险性分类模型; FCM-SVM模型能够对其他该地区泥石流沟谷进行评价, 具有单纯聚类算法所没有的优势。

5)相比传统方法的泥石流沟谷分类形式, FCM-SVM方法的分级边界随研究区域整体情况变化, 使危险性评价工作更加科学合理, 也有利于当地泥石流灾害预警及治理工作的开展。