技术领域
[0001] 本
发明涉及
自然灾害风险评价技术领域，具体的说是基于GIS与Logistic回归模型的泥石流易发性预测方法。
背景技术
[0002] 泥石流是山区常见的自然灾害现象，具有爆发突然、运动迅速、历时短暂、破坏性强等特点，会给当地居民的生命财产安全带来巨大的威胁。而泥石流的易发性评价研究一直是国内外研究的重要课题。泥石流易发性评价是指在现状自然地理、地质构造和
地层岩性等成灾背景条件下泥石流的易发程度，对可能发生的泥石流的概率进行定性和定量分析。易发性评价对泥石流灾害的管理和防治具有重要的指导意义。
[0003] 现有的泥石流易发性评价方法主要分为定性方法和定量方法，定性方法主要是指基于专家经验的现场分析和因子专题图
叠加分析，定量方法主要是指统计分析法。定性方法由于过多的涉及主观意愿及受专业人士
水平的影响，评价结果难以让人信服。而统计分析法中的Logistic回归由于方法简单，物理意义明确，很多学者利用Logistic回归模型来预测评价地质灾害的易发性，取得了不错成绩。但现有的Logistic回归模型存在评价
精度较低、无法处理变量间的共线性问题等不足。另外，过往大多关于泥石流易发性的研究最终停留在泥石流易发程度的判断，并没有识别出不同类型泥石流的主要影响因素，这给防治工作的开展造成不便。
发明内容
[0004] 针对
现有技术中存在的上述不足之处，本发明目的是提供一种评价精度更高、针对性更强的泥石流灾害易发性评价方法。
[0005] 本发明为实现上述目的所采用的技术方案是：基于GIS与Logistic回归模型的泥石流易发性预测方法，包括如下步骤：
[0006] 步骤1、确定评价单元：
[0007] 所述评价单元的确定方式为，基于数字高程模型数据，利用ArcGIS平台提取小流域，并结合遥感图像进行人为
修改，将修改后的小流域确定为所述评价单元，且相应的具有泥石流灾害发生的历史数据；
[0008] 步骤2、评价参数的选取与处理：
[0009] 选取评价参数，并对所选的评价参数进行无量纲化处理得到评价因子；
[0010] 步骤3、构建训练样本集：
[0011] 所述训练样本集依据所述评价单元内泥石流灾害发生的历史数据以及所述评价因子组成；
[0013] 采用模糊C均值聚类法对训练样本集进行聚类分析，使同类子样本间的相似度最大，而不同子样本间的差异性最明显；
[0014] 步骤5、因子分析：
[0015] 依据所述聚类分析结果，分别对每一组聚类分析后的子样本进行因子分析得出每一组的公因子及计算相应的因子得分；
[0016] 步骤6、建立Logistic回归模型：
[0017] 依据每组所述因子分析结果作为自变量以及所述评价单元内泥石流的发生情况作为因变量建立Logistic回归模型；
[0018] 步骤7、Logistic回归模型的分析：
[0019] 对所建立的Logistic回归模型进行拟合度分析、显著性分析以及回代法误判率分析；
[0020] 步骤8、评价结果的分类：
[0021] 采用自然断点法对评价结果进行分类，分为高易发性、中易发性、较低易发性及低易发性四类；
[0022] 步骤9、预测待评价区域内泥石流灾害的易发性：
[0023] 具体地，获取待评价区域内各评价因子的评价指标值输入已建立好的Logistic回归模型，计算出各评价单元发生泥石流灾害的概率。最后，根据已划分好的易发性等级进行分类。
[0024] 所述模糊C均值聚类法所涉及的方程式为：
[0025]
[0026] 所述模糊C均值聚类法分析步骤为，
[0027] 步骤1、设置计算参数：
[0028] 设置模糊指数m，最大
迭代次数maxTimes以及目标函数精度ε；
[0029] 步骤2、确定聚类数目：
[0030] 引入聚类有效性函数Vcs确定聚类中心数目；
[0031] 步骤3、初始化聚类中心；
[0032] 步骤4、更新隶属度矩阵和聚类中心：
[0033] 更新隶属度矩阵和聚类中心，若样本点和聚类中心的距离为0，则将该点与对应类的隶属度值设为1；
[0034] 步骤5、计算距离及目标函数值J：
[0035] 计算每个样本点到每个类的聚类中心的距离，并计算目标函数值J；
[0036] 步骤6、迭代计算：
[0037] 若计算达到最大迭代次数或前后两次J的绝对值差小于ε，则停止计算，否则返回步骤4；
[0038] 步骤7、按最大隶属度原理判断样本点归属。
[0039] 所述聚类有效性函数Vcs由紧致度和离散度组成，其中，紧致度的定义如下：
[0040]
[0041] 其中，
[0042] 所述因子分析的主要步骤为，
[0043] 步骤1、相关性分析：
[0044] 因子分析前，首先进行KMO检验和巴特利球体检验。在因子分析中，若拒绝原假设，则说明可以做因子分析，若不拒绝原假设，则说明这些变量可能独立提供一些信息，不适合做因子分析；
[0045] 步骤2、构造公因子F及因子荷载矩阵：
[0046] 计算标准化指标变量的协方差矩阵S和相关系数矩阵R，求矩阵R的特征值λk及
特征向量uk,并按公式 求出方差贡献率，将方差贡献率累加求出累积贡献率，选取m个公共因子。一般取累积贡献率达85％或者特征值大于1的λ1,λ2,…,λm所对应的第一，第二,…,第m(m≤p)个因子为公共因子。最后求出因子
载荷阵A，计算所需公式如下：
[0047] S＝Z`Z`/p-1；
[0048]
[0049] A＝(akj)p×m＝(ukj(λk)1/2)p×m；
[0050] X＝AF+ε；
[0051] 步骤3、公因子旋转：
[0052] 采用方差极大
正交旋转技术进行公因子旋转；
[0053] 步骤4、计算因子得分：
[0054] 计算因子得分采用是汤姆森回归法。
[0055] 所述评价参数包括流域面积、河网
密度、平均坡度、最大高程差、主沟比降、圆度、
断层密度、距断层的距离、距公路的距离、垦殖系数、归一化植被系数、7天最大雨量12类，且12个评价单元并不是相互独立的，而是具有一定的共线性，相互之间有影响。
[0056] 根据所述评价单元内泥石流灾害发生的历史数据，把曾经有泥石流灾害发生的评价单元设置为1；把没有泥石流灾害发生的评价单元设置为0。
[0057] 本发明的有益效果：利用评价单元内泥石流灾害发生的历史数据以及泥石流易发性评价参数，结合模糊C均值聚类和因子分析建立Logistic回归模型，并采用该Logistic回归模型对待评价区域内泥石流灾害的易发性进行评价。与现有的其他泥石流易发性评价方法相比，该方法通过将Logistic回归与模糊C均值聚类和因子分析法有机结合，不但很好地提高了模型的准确率，而且解决了评价因子间存在的共线性问题，进一步识别不同类型泥石流的主要影响因子，为防治工程的实施提供新的思路和依据。
附图说明
[0058] 图1为本发明结构示意图；
[0059] 图2为本发明结构示意图；
[0060] 图3为本发明结构示意图。
具体实施方式
[0061] 下面将结合本发明
实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0062] 实施例一
[0063] 请参阅图1，基于GIS与Logistic回归模型的泥石流易发性预测方法，包括如下步骤：
[0064] 步骤1、确定评价单元：
[0065] 评价单元的确定方式为，基于数字高程模型数据，利用ArcGIS平台提取小流域，并结合遥感图像进行人为修改，将修改后的小流域确定为评价单元，且相应的具有泥石流灾害发生的历史数据，根据评价单元内泥石流灾害发生的历史数据，把曾经有泥石流灾害发生的评价单元设置为1；把没有泥石流灾害发生的评价单元设置为0；
[0066] 步骤2、评价参数的选取与处理：
[0067] 选取评价参数，并对所选的评价参数进行无量纲化处理得到评价因子；
[0068] 步骤3、构建训练样本集：
[0069] 训练样本集依据评价单元内泥石流灾害发生的历史数据以及评价因子组成；
[0070] 步骤4、聚类分析：
[0071] 采用模糊C均值聚类法对训练样本集进行聚类分析，使同类子样本间的相似度最大，而不同子样本间的差异性最明显；
[0072] 步骤5、因子分析：
[0073] 依据聚类分析结果，分别对每一组聚类分析后的子样本进行因子分析得出每一组的公因子及计算相应的因子得分；
[0074] 步骤6、建立Logistic回归模型：
[0075] 依据每组因子分析结果作为自变量以及评价单元内泥石流的发生情况作为因变量建立Logistic回归模型；
[0076] 步骤7、Logistic回归模型的分析：
[0077] 对所建立的Logistic回归模型进行拟合度分析、显著性分析以及回代法误判率分析；
[0078] 步骤8、评价结果的分类：
[0079] 采用自然断点法对评价结果进行分类，分为高易发性、中易发性、较低易发性及低易发性四类；
[0080] 步骤9、预测待评价区域内泥石流灾害的易发性：
[0081] 具体地，获取待评价区域内各评价因子的评价指标值输入已建立好的Logistic回归模型，计算出各评价单元发生泥石流灾害的概率。最后，根据已划分好的易发性等级进行分类。
[0082] 综上，本发明得出的最终结果包括：1.根据因子分析结果中因子荷载矩阵来判断各个待评价样本的主要影响因素；2.易发性评价结果，为泥石流防治工作的开展提供新的思路。
[0083] 实施例二
[0084] 请参阅图1，基于GIS与Logistic回归模型的泥石流易发性预测方法，在实施例一中所描述的模糊C均值聚类法所涉及的方程式为：
[0085] 公式一：
[0086] 式中，m为大于1的任意实数(一般取m为2)，xi为测量数据的第i个分量，Pj是第j个聚类中心，||·||为距离测度。N表示样本数目，c表示聚类数目；uij表示xi在聚类j中的隶属度，满足uij∈[0,1],且 (j＝1，2，····，n)uij和Pj更新公式分别公式二和公式三，
[0087] 公式二：
[0088] 公式三：
[0089] 当满足 时，停止计算，
算法收敛于局部最小值。
[0090] 在实施例一中，模糊C均值聚类法分析步骤为，
[0091] 步骤1、设置计算参数：
[0092] 设置模糊指数m(本发明设置为2)，最大迭代次数maxTimes以及目标函数精度ε；
[0093] 步骤2、确定聚类数目：
[0094] 引入聚类有效性函数Vcs确定聚类中心数目；聚类有效性函数Vcs由紧致度和离散度组成，其中，紧致度的定义如下：
[0095] 公式四：
[0096] 其中，
[0097] Cij是第i类和第j类样本间的紧致度。当uij≥1/c时，说明第j类样本属于第i个类的隶属度较高，即第i类和第j类样本间的紧致度较高。当uij＜1/c时，对应的紧致度变为0，因为uij太小意味着第j类样本不太可能属于第i类，对紧致度测量Cij无意义。然后对所有的样本同所有的类别间的紧致度求和，这样当所有的样本明确属于某个类时，紧致度最大，即聚类结果紧致。
[0098] 而第i类和第j类之间的离散度定义为：
[0099] 公式五：Sij＝min(uik,ujk),，其中k＝1,2,···，n；
[0100] 即两个类之间的离散度定义为样本属于这两个类的隶属度中最小的值。
[0101] 总离散度的定义为：
[0102] 公式六：
[0103] 总离散度用于衡量界线最不明确的两个类之间的离散度。当离散度越小，说明两个类之间的差别越大，即所有的类划分时越明确，聚类结果越好。
[0104] 在此
基础上，定义基于隶属度的聚类有效性指数Vcs,即：
[0105] 公式七：
[0106] 综上所述，当紧致度越大、离散度越小时Vcs指数的值越大，说明聚类效果越好。
[0107] 步骤3、初始化聚类中心；
[0108] 步骤4、更新隶属度矩阵和聚类中心：
[0109] 根据公式二和公式三分别更新隶属度矩阵很聚类中心，若样本点和聚类中心的距离为0，则将该点与对应类的隶属度值设为1；
[0110] 步骤5、计算距离及目标函数值J：
[0111] 计算每个样本点到每个类的聚类中心的距离，并根据公式一计算目标函数值J；
[0112] 步骤6、迭代计算：
[0113] 若计算达到最大迭代次数或前后两次J的绝对值差小于ε，则停止计算，否则返回步骤4；
[0114] 步骤7、按最大隶属度原理判断样本点归属。
[0115] 实施例三
[0116] 请参阅图1，基于GIS与Logistic回归模型的泥石流易发性预测方法，在实施例一中，因子分析的主要步骤为，
[0117] 步骤1、相关性分析：
[0118] 因子分析前，首先进行KMO检验和巴特利球体检验。在因子分析中，若拒绝原假设，则说明可以做因子分析，若不拒绝原假设，则说明这些变量可能独立提供一些信息，不适合做因子分析；
[0119] 步骤2、构造公因子F及因子荷载矩阵：
[0120] 计算标准化指标变量的协方差矩阵S和相关系数矩阵R，求矩阵R的特征值λk及特征向量uk,并按公式 求出方差贡献率，将方差贡献率累加求出累积贡献率，选取m个公共因子。一般取累积贡献率达85％或者特征值大于1的λ1,λ2,…,λm所对应的第一，第二,…,第m(m≤p)个因子为公共因子。最后求出因子载荷阵A，计算所需公式如下：
[0121] 公式八：S＝Z`Z`/p-1；
[0122] 公式九：
[0123] 公式十：A＝(akj)p×m＝(ukj(λk)1/2)p×m；
[0124] 公式十一：X＝AF+ε；
[0125] 步骤3、公因子旋转：
[0126] 建立因子分析模型的目的不仅在于要找公因子，更重要的是知道每个公共因子的意义，以便对实际问题进行分析。本发明采用方差极大正交旋转技术进行公因子旋转；
[0127] 步骤4、计算因子得分：
[0128] 计算因子得分最常用的方法是汤姆森回归法，该方法假设公共因子可再对P个原始变量做回归分析，将公共因子表示为原始变量的线性组合，即：
[0129] 公式十二：
[0130] 公式中：W＝A`R-1称为因子得分系数矩阵。
[0131] 实施例四
[0132] 请参阅图1，基于GIS与Logistic回归模型的泥石流易发性预测方法，在实施例一种，评价参数包括流域面积、河网密度、平均坡度、最大高程差、主沟比降、圆度、断层密度、距断层的距离、距公路的距离、垦殖系数、归一化植被系数、7天最大雨量12类，且12个评价单元并不是相互独立的，而是具有一定的共线性，相互之间有影响，而评价参数具体如下：
[0133] 流域面积：利用ArcGIS中的几何运算得到流域面积，并结合遥感技术进行校正。
[0134] 河网密度：一定面积内河网长度，利用ArcGIS几何运算进行统计。
[0135] 平均坡度：利用ArcGIS坡度计算功能，并运用分区统计工具中的MEAN统计评价单元中的平均值作为该评价单元的值。
[0136] 最大高程差：每个流域的最高点和最低点的海拔高度由DEM确定，运用分区统计工具中的RANGE统计各单元的最大高程差。
[0137] 主沟比降：主沟的最大高差与其直线长度之比。
[0138] 圆度：它是指盆地面积与周长等于盆地周长的圆的面积之比。
[0139] 断层密度：一定面积内断层的长度。
[0140] 距断层的距离：本发明利用ArcGIS以200m为间距作多层缓冲分析，统计评价单元中的平均值作为该评价单元的值。
[0141] 距公路的距离：本发明利用ArcGIS以200m为间距作多层缓冲分析，统计评价单元中的平均值作为该评价单元的值。
[0142] 垦殖系数：一定区域内的耕地面积占土地总面积的比例。本发明结合Google影像及ArcGIS统计评价单元内耕地面积。
[0143] 归一化植被系数：植被归一化指数由MODIS产品提供栅格数据，本发明统计评价单元中的平均值作为评价单元值。
[0144] 7天最大雨量：以研究区附近12个降水站的资料为参考，在ArcGIS中采用普通克里格插值法进行插值。并统计评价单元中的平均值作为评价单元值。
[0145] 实施例五
[0146] 请参阅图1，基于GIS与Logistic回归模型的泥石流易发性预测方法，利用评价单元内泥石流灾害发生的历史数据以及泥石流易发性评价参数，结合模糊C均值聚类和因子分析建立Logistic回归模型，并采用该Logistic回归模型对待评价区域内泥石流灾害的易发性进行评价。与现有的其他泥石流易发性评价方法相比，该方法通过将Logistic回归与模糊C均值聚类和因子分析法有机结合，不但很好地提高了模型的准确率，而且解决了评价因子间存在的共线性问题，进一步识别不同类型泥石流的主要影响因子，为防治工程的实施提供新的思路和依据。
[0147] 对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附
权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
[0148] 此外，应当理解，虽然本
说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。
