CN109389680A - 泥石流流域关键地形因子特征筛选方法 - Google Patents

Abstract

针对现有技术不能提供具体的判别指标和各因子的综合作用效果的缺陷，本发明公开一种泥石流流域关键地形因子特征筛选方法。本方法先从研究区域泥石流流域样本与非泥石流流域样本中采集可能关键地形因子特征集，然后建立Lasso惩罚函数，求解惩罚函数的最优解，最优解特征中β≠0回归系数所对应的特征即为所选关键地形因子特征。本发明还提供泥石流流域关键地形因子特征筛选方法的应用。本发明方法是基于Lasso方法对比分析筛选泥石流流域关键地形因子特征的方法，利用机器学习算法与程序完成定量分析，克服了现有技术因以统计方法为基础而只能定性地筛选关键地形因子特征的弊端。方法原理可靠，结果准确性高。

Description

泥石流流域关键地形因子特征筛选方法

技术领域

本发明涉及一种地形因子特征筛选方法，特别是涉及一种筛选泥石流流域区别于非泥石流流域的关键地形因子特征的方法，属于地质地貌测绘勘察技术、地质灾害防治技术领域。

背景技术

泥石流作为流域演化到一定阶段所出现的一种自然现象，物质和能量条件是决定其发育和分布的根本因素，而物质和能量条件最大程度上决定于流域地形条件。泥石流一般为小流域事件，因而流域的局部特征对水文过程能产生决定性的影响，许多小流域的滑坡、崩塌产生的物质，往往通过沟谷，以泥石流的形式输移到流域外部，形成坡-谷耦合的系统，而泥石流则是这种耦合关系的一种表现。凡是地形条件具备地区，都有发生泥石流的可能。通过比对分析典型泥石流流域和非泥石流流域的地形特征，去除不相关特征，能够确定决定泥石流发育的流域关键地形特征(以下称泥石流流域关键地形因子特征)。掌握了泥石流流域关键地形因子特征，一方面有利于构建更为简单和容易解释的泥石流流域地形研究模型，为基于形成机理的泥石流流域判识和监测研究奠定基础；另一方面可以对尚未发生过泥石流的流域实施全面调查，排查区分泥石流流域与非泥石流流域。两方面都对泥石流灾害防治具有重要意义。

现有技术对泥石流流域关键地形因子特征的筛选是采用数理统计的手段进行分析研究(Coe,J.A.and Godt,J.W.,2001,Debris flows triggered by the Elrainstorm of February 2-3,1998,Walpert Ridge and Vicinity；郑胥智.降雨及地文条件与土石流发生关系之研究[D].中国台湾：成功大学硕士论文,2003；曾奕超.土石流发生降雨及地文综合警戒指标之研究 [D].中国台湾：成功大学硕士论文,2004；Liu,C.N.,Huang H.F.,Dong J.J. Impacts of September 21,1999Chi-Chi earthquake on thecharacteristics of gully-type debris flows in central Taiwan.Natural Hazards,2008,47.)。这些研究受技术手段等的制约，当前对泥石流流域关键地形因子特征的筛选多凭借经验确定，属定性分析。同时采用传统数理统计分析的方法集中对各因素进行单独分析。现有技术没能提供具体的判别指标和各因子的综合作用效果，难以形成统一的具有代表性的结论。

发明内容

本发明的目的就是针对现有技术的不足，提供一种基于Lasso方法对比分析筛选泥石流流域关键地形因子特征的方法，该方法利用机器学习算法与程序完成定量分析，确定泥石流流域区别于非泥石流流域的关键地形因子特征。

为实现上述目的，本发明首先提供一种泥石流流域关键地形因子特征筛选方法，其技术方案如下：

一种泥石流流域关键地形因子特征筛选方法，用于确定泥石流流域区别于非泥石流流域的关键地形因子特征，其特征在于：

步骤S1、划定研究区域，调查获取研究区域内的O个泥石流流域样本与P个非泥石流流域样本，O等于或不等于P，分别为两类样本标记流域标签；

步骤S2、分别调查获取两类样本的流域地形特征值，记录Q个两类流域样本地形因子特征为泥石流流域区别于非泥石流流域的可能关键地形因子特征，Q不大于O与P之和，对各可能关键地形因子特征值作中心标准化处理后构成可能关键地形因子特征集A；

步骤S3、依据Lasso准则，建立以集合A中可能关键地形因子特征为自变量X、流域标签为因变量Y的惩罚函数；

步骤S4、以惩罚函数最小为标准求解所述惩罚函数的最优解，得到惩罚函数的自变量活动集B；

步骤S5、将自变量活动集B中β≠0回归系数所对应的特征放入集合C，得到关键地形因子特征集C；集合C内的特征即为泥石流流域区别于非泥石流流域的关键地形因子特征。

在筛选泥石流流域关键地形因子特征中，面对的问题主要是泥石流流域地形的复杂性以及大量数据间存在错综复杂的关系，容易使构建关系矩阵成为对扰动敏感的病态阵，不能进行有效的特征筛选。本发明上述方法采用机器学习的Lasso方法实现。其基本思想是，依Lasso规则(在回归系数绝对值之和小于一个常数的约束条件下，使残差平方和最小化) 构建一个以可能关键地形因子特征为自变量，以流域属性为因变量的惩罚函数，通过训练，在约束条件下得到一个较为精炼的模型。通过该处理，能较好解决由于样本的共线性容易使矩阵成为对扰动敏感的病态阵，进而造成回归系数解析解的不稳定，解出的参数常常具有稀疏的特征，即很多特征对应的参数会为零的问题。这样使得特征选择成为可能，能够筛选出对结果最具影响的地形特征。

上述方法步骤S1中，可以利用研究区数字高程模型(DEM)划分研究区域内泥石流流域样本与非泥石流流域样本，也可围绕形成泥石流的物质与能量条件，通过现场调查、遥感影像解译手段划分泥石流流域样本与非泥石流流域样本。

上述方法步骤S3中，Lasso惩罚函数是：

t≥0 式3

其中，p—可能关键地形因子特征总数，t—给定的约束参数，N—泥石流流域样本与泥石流流域样本总数，y_i—因变量Y的第i个样本所对应的流域标签，x_ij—第i行第j列的自变量X值，β＝(β₁,…,β_j,…,β_p)^T—中心标准化处理后的各可能关键地形因子特征待解算的回归系数，β_j—第j个可能关键地形因子特征所对应的待解算的回归系数，β₀—常数项。

上述方法步骤S4中，依据Lasso规则，先将初始相关系数设定为0，先找出与残差值相关系数最大的变量，其标准是使惩罚函数最小，该变量设为x₁，将其加入到自变量活动集B，接着在x₁方向找到一个最长的步长，使得出现下一个变量x₂和残差的相关系数与x₁和残差的相关系数相等，将x₂加入活动集。用同样的方法，依次找到剩余变量中与前面加入活动集的变量和残差的相关系数相等的变量，加入活动集B。

利用上述方法筛选出的泥石流流域关键地形因子特征是从诸多可能特征错综交错的相互关系中经由定量分析确定，克服了传统从数量统计方法对各因素进行单独分析的缺点，因而更为准确。将本发明方法筛选出的关键地形因子特征付诸于对研究区域内流域的数学模型研究也能为得到研究区域地形更精确的数学结构提供基础。故，本发明还提供上述泥石流流域关键地形因子特征筛选方法在泥石流流域研究中的应用。

与现有技术相比，本发明的有益效果是：(1)提供了一种基于Lasso 方法对比分析筛选泥石流流域关键地形因子特征的方法，该方法利用机器学习算法与程序完成定量分析，确定出泥石流流域区别于非泥石流流域的关键地形因子特征，克服了现有技术因以统计方法为基础而只能定性地筛选关键地形因子特征的弊端。(2)本方法通过构建一个以可能关键地形因子特征为自变量，以流域属性为因变量的Lasso惩罚函数，得到一个较为精炼的模型，该模型压缩了部分系数并设定部分系数为零，从而较好地解决了由于样本的共线性容易使矩阵成为对扰动敏感的病态阵，进而造成回归系数解析解的不稳定，解出的参数常常具有稀疏特征，即很多特征对应的参数会为零的问题。保证了特征筛选的有效性。(3)提供了本发明泥石流流域关键地形因子特征筛选方法在泥石流流域研究中的应用。

附图说明

图1是凉山州境内典型泥石流流域和非泥石流流域分布图。

图2a是特征收敛速度图。

图2b是图2a的局部放大图。

图3是特征选择概率图。

具体实施方式

下面结合附图，对本发明的优选实施例作进一步的描述。

实施例一

如图1～图3所示，用本发明方法在四川省凉山彝族自治州境内筛选泥石流流域的关键地形因子特征。

1、筛选关键地形因子特征

步骤S1、前期流域调查

利用数字高程模型(DEM)，划分凉山州境内的流域，其中包含已知的典型泥石流流域与非泥石流流域(图1是凉山州境内典型泥石流流域与非泥石流流域分布图)。得到泥石流流域样本O＝679个，非泥石流流域样本P＝512个。为泥石流流域样本标记流域标签(+1)，为非泥石流流域样本标记流域标签(-1)。

步骤S2、分别提取两类样本的流域地形因子值

进一步利用DEM获取两类流域1191个样本的各种流域地形因子特征Q＝14个，记为可能关键地形因子特征。Q(＝14)不大于O(＝679) 与P(＝512)两者间的较小者，符合计算要求。14个地形因子特征分别是：平均高程(1)、坡度能量(2)、沟床比降(3)、高程方差(4)、相对高差(5)、最大坡度(6)、平均坡度(7)、坡度方差(8)、Melton ratio (9)、流域面积(10)、流域周长(11)、流域长度(12)、流域形状系数 (13)、流域形状因子(14)。提取可能地形因子特征时，对于只在一种流域样本中出现的地形因子特征，可在另一流域样本中设置为0。

对各可能关键地形因子特征值作中心标准化处理，将各特征值转化为处于0～1之间，使数据值原始特征转化为无量纲的量。标准化后的可能关键地形因子特征构成集合A。

步骤S3、依据Lasso准则建立惩罚函数

依据Lasso准则，建立以集合A中可能关键地形因子特征为自变量 X、流域标签为因变量Y的惩罚函数。

t≥0 式3

式中，p—可能关键地形因子特征总数，步骤S2确定，

t—给定的约束参数，在0～1间试算确定，

N—泥石流流域样本与非泥石流流域样本总数，步骤S2确定，

y_i—因变量Y的第i个样本所对应的流域标签，步骤S1确定，

x_ij—第i行第j列的自变量X值，步骤S2确定，

β＝(β₁,…,β_j,…,β_p)^T—中心标准化处理后的各可能关键地形因子特征待解算的回归系数，

β_j—第j个可能关键地形因子特征所对应的待解算的回归系数，

β₀—常数项，初值为0。

步骤S4、以Lasso惩罚函数最小为标准求惩罚函数的最优解，加入变量活动集B；

步骤S41、令β＝(β₁,…,β_j,…,β_p)^T各元素的初值为0，自变量活动集B的初始值为空；

步骤S42、以使惩罚函数最小为标准，确定与残差相关系数最大的自变量X的第i行元素(x_i＝(x_i1,…,x_ij,…,x_ip)^T，将该自变量 x_i的地形因子值加入自变量活动集B；残差中，为回归拟合值，初值为0；

沿着自变量活动集B中变量前进，逐点增大自变量x_i对应的系数β_i，以降低自变量x_i与残差的相关系数，此时回归拟合值残差为

继续计算新的残差与各自变量X之间的相关系数，随着系数β_i的逐点增大，自变量x_i与新的残差的相关系数逐渐减小，直到出现一个新的自变量x_k满足内积然后沿矢量x_i与x_k的角平分线方向前进，对应地逐点增大β_i与β_k，相应残差为

如此循环找到第三个自变量，确定三个自变量中第一自变量与第二自变量夹角的第一角平分线，第二自变量与第三自变量夹角的第二角平分线，沿两条角平分线组成夹角的角平分线方向逐点移动，调整三个自变量的系数，使得在移动过程中三个自变量与当前残差的相关系数均减小且一直相等；

步骤S43、经上述依次逐点调整自变量活动集B中各个自变量X对应的回归系数β，使自变量活动集B中各自变量与当前残差的相关系数均减小且一直相等，残差为其中n为自变量活动集B 中的自变量个数；

步骤S44、确定自变量活动集B中各自变量对应的非零系数是否变为0，若是，将该系数为0的自变量从自变量活动集B中剔除，并继续执行S45；若否，继续执行S45；

步骤S45、确定是否存在自变量活动集B外自变量X各行元素与当前残差的相关系数等于自变量活动集B中各自变量与当前残差的相关系数的情况若是，则继续执行S46；否则，返回执行S43；

步骤S46、将自变量活动集B外的自变量X中的该行元素加入自变量活动集B；

步骤S47、确定自变量X中是否有未被处理的行元素，若是，则返回执行步骤S43；若否，结束操作。

步骤S5、将自变量活动集B中符合预设筛选规则的因子放入关键地形因子特征集C，预设筛选规则是：对于步骤S4中解出各可能的关键地形因子特征待解算的回归系数β＝(β₁,…,β_j,…,β_p)^T，β中非零元素对应的特征放入集合C；集合C内的特征即为泥石流流域区别于非泥石流流域的关键地形因子特征。

经本发明方法处理后，根据各流域地形特征收敛速度(图2a是特征收敛速度图，图2b是图2a的局部放大图)以及关键特征的选择概率(图 3是特征选择概率图)，确定泥石流流域的关键地形特征4个，分别是：平均高程(1)、高程方差(4)、流域形状因子(8)、Meltonratio(14)。

2、SVM及k-fold交叉检验

将选取的泥石流和非泥石流流域地形特征数据分作两类：一类是经过本发明方法选择的特征，另一类未经选择的原始特征。

分别将两类特征均匀分成K组，对每类特征数据，利用SVM分类方法将每个地形特征子集数据分别做一次检验集，其余的K-1组子集数据作为训练集，得到K个模型，用这K个模型最终的检验集的分类准确率平均数作为此K-fold下分类器的性能指标。对比两种分类下的性能指标，从而确定特征选择的可靠性。

分别对各流域的原始特征和选择后的特征利用SVM分类器进行k-folk交叉检验，经过测试，k值取7，结果如下：

以经过选择后的4个特征为输入，分类结果准确率为57.9818％；以未经选择的14个原始特征为输入，分类结果准确率为57.9405％。根据以上结果检验，可以看出，两种结果差别不大，因此经过选择的4个特征应该就是泥石流流域所具有的关键地形特征，分别为：平均高程(1)、高程方差(4)、流域形状因子(8)、Melton ratio(14)。

Claims (9)

1.泥石流流域关键地形因子特征筛选方法，用于确定泥石流流域区别于非泥石流流域的关键地形因子特征，其特征在于：

步骤S1、划定研究区域，调查获取研究区域内的O个泥石流流域样本与P个非泥石流流域样本，O等于或不等于P，分别为两类流域样本标记流域标签；

步骤S2、分别提取两类流域样本的流域地形特征值，记录两类流域样本的Q个地形因子特征为泥石流流域区别于非泥石流流域的可能关键地形因子特征，Q不大于O与P之和，对各可能关键地形因子特征值作中心标准化处理后构成可能关键地形因子特征集A；

步骤S3、依据Lasso准则，建立以集合A中可能关键地形因子特征为自变量X、流域标签为因变量Y的惩罚函数；

步骤S4、以惩罚函数最小为标准求解所述惩罚函数的最优解，得到惩罚函数的自变量活动集B；

步骤S5、将自变量活动集B中特征的回归系数不为0的特征放入集合C，得到关键地形因子特征集C；集合C内的特征即为泥石流流域区别于非泥石流流域的关键地形因子特征。

2.根据权利要求1所述的方法，其特征在于：

所述步骤S1中，泥石流流域样本标记的流域标签为+1，非泥石流流域样本的流域标签为-1；

所述步骤S3中，构建的惩罚函数是：

t≥0 式3

式中，p—可能关键地形因子特征总数，步骤S2确定，

t—给定的约束参数，在0～1间试算确定，

N—泥石流流域样本与非泥石流流域样本总数，步骤S2确定，

y_i—因变量Y的第i个样本所对应的流域标签，步骤S1确定，

x_ij—第i行第j列的自变量X值，步骤S2确定，

β＝(β₁,…,β_j,…,β_p)^T—中心标准化处理后的各可能关键地形因子特征待解算的回归系数，

β_j—第j个可能关键地形因子特征所对应的待解算的回归系数，

β₀—常数项，初值为0。

3.根据权利要求2所述的方法，其特征在于：所述步骤S4中，依如下步骤实施：

步骤S41、令β＝(β₁,…,β_j,…,β_p)^T各元素的初值为0，自变量活动集B的初始值为空；

步骤S42、以使惩罚函数最小为标准，确定与残差相关系数最大的自变量X的第i行元素(x_i＝(x_i1,…,x_ij,…,x_ip)^T，将该自变量x_i的地形因子值加入自变量活动集B，所述残差中，为回归拟合值，初值为零；

沿着自变量活动集B中变量前进，逐点增大自变量x_i对应的系数β_i，以降低自变量x_i与残差的相关系数，此时回归拟合值残差为

继续计算新的残差与各自变量X之间的相关系数，随着系数β_i的逐点增大，自变量x_i与新的残差的相关系数逐渐减小，直到出现一个新的自变量x_k满足内积然后沿矢量x_i与x_k的角平分线方向前进，对应地逐点增大β_i与β_k，相应残差为

如此循环找到第三个自变量，确定三个自变量中第一自变量与第二自变量夹角的第一角平分线，第二自变量与第三自变量夹角的第二角平分线，沿两条角平分线组成夹角的角平分线方向逐点移动，调整三个自变量的系数，使得在移动过程中三个自变量与当前残差的相关系数均减小且一直相等；

步骤S43、经上述依次逐点调整自变量活动集B中各个自变量X对应的回归系数β，使自变量活动集B中各自变量与当前残差的相关系数均减小且一直相等，残差为其中n为自变量活动集B中的自变量个数；步骤S44、确定自变量活动集B中各自变量对应的非零系数是否变为0，若是，将该系数为0的自变量从自变量活动集B中剔除，并继续执行S45；若否，继续执行S45；

步骤S45、确定是否存在自变量活动集B外自变量X各行元素与当前残差的相关系数等于自变量活动集B中各自变量与当前残差的相关系数的情况若是，则继续执行S46；否则，返回执行S43；

步骤S46、将自变量活动集B外的自变量X中的该行元素加入自变量活动集B；

步骤S47、确定自变量X中是否有未被处理的行元素，若是，则返回执行步骤S43；若否，结束操作。

4.根据权利要求3所述的方法，其特征在于：所述步骤S5中，对于步骤S4中解出各可能的关键地形因子特征待解算的回归系数β＝(β₁,…,β_j,…,β_p)^T，β中非零元素对应的特征即为关键地形因子特征。

5.根据权利要求1所述的方法，其特征在于：所述步骤S4中，求解所述惩罚函数的最优解采用LASSO算法对目标函数进行有约束的最小化计算，获得回归向量。

6.根据权利要求5所述的方法，其特征在于：通过修正的LARS求解所述惩罚函数，得到惩罚函数的自变量活动集B。

7.根据权利要求1所述的方法，其特征在于：所述步骤S1中，可能的关键地形因子特征包括平均高程、坡度能量、沟床比降、高程方差、最大坡度、相对高差、流域形状系数、流域形状因子、坡度方差、平均坡度、流域周长、流域面积、流域长度、Melton ratio。

8.根据权利要求1所述的方法，其特征在于：所述步骤S1中，利用数字高程模型划分研究区域内泥石流流域样本与非泥石流流域样本，或者围绕形成泥石流的物质与能量条件，通过现场调查、遥感影像解译手段划分泥石流流域与非泥石流流域，获取流域样本。

9.权利要求1～8所述泥石流流域关键地形因子特征筛选方法在泥石流流域研究中的应用。