一种滑坡稳定性时变规律分析的方法与流程
本发明涉及地质滑坡灾害监测技术领域，尤其涉及一种滑坡稳定性时变规律分析的方法。
背景技术：
滑坡在全球范围内分布广泛，危害极为严重，已成为仅次于地震的第二大地质灾害。一方面，滑坡的稳定性决定了其对外界触发因素的响应程度；另一方面，滑坡稳定性也在外界触发因素的作用下不断发生变化，因此，滑坡稳定性分析是危险度评价、提高位移预测、灾害预警的重要内容。
目前，滑坡状稳定性研究方法以极限平衡法和有限元法为主。尽管极限平衡法模型简单，计算简洁，但需要事先知道滑动面的大致位置和形状，所得安全系数只是假定滑裂面上的平均安全度，无法判断滑坡的变形模式。有限元方法不仅能够满足静力平衡等条件，同时能够将无限自由度结构转换为有限自由度的等价，具有非线性分析能力，因此得到了广泛的关注。但，有限元方法的适用性同样需要付出一定代价，比如计算量大，受物理参数影响大等。
滑坡的变形和破坏是一个渐进的过程，滑坡的稳定性也随着时间推移不断变化的。稳定性的变化不仅表示了滑坡对外界因素的响应程度，而且能够反映滑坡内在的状况及变形趋势。然而现有边坡稳定性分析方法多局限为静态分析方法。虽然已有学者探索滑坡实时稳定性的研究方法，但研究尚不充分，且涉及参数众多，计算复杂；分析结果局限为稳定性系数大于1表示稳定，反之则不稳定的两种状态，缺乏度的反映，进而导致不能准确给出预警。
技术实现要素：
有鉴于此，本发明的实施例提供了一种计算简单，能预测滑坡稳定性发展趋势，为滑坡高效预测预报提供依据的滑坡稳定性时变规律分析的方法。
本发明的实施例提供一种滑坡稳定性时变规律分析的方法，包括以下步骤：
(1)根据滑坡每月的变形速度和变形加速度构成每月的二维特征数据集；
(2)利用K均值算法对每月的二维特征数据集进行聚类分析得到类簇，并为每个类簇分配聚类标签；
(3)对各个类簇分别建立云模型，并计算二维特征数据集内二维特征数据相对于云模型的隶属度，根据隶属度对步骤(2)的聚类标签进行修正，并将修正后的聚类标签所代表的形变状态作为滑坡该月的形变状态，以此类推得到滑坡每月的形变状态；
(4)根据滑坡每月的形变状态对滑坡的形变状态进行统计，并据此计算滑坡每月形变状态的信息量模型，所述信息量模型包括状态出现信息量和状态转移信息量，将状态出现信息量、状态转移信息量和云模型中二维特征数据到对应云模型的隶属度乘积定义为马尔科夫信息量，据此计算每月的马尔科夫信息量；
(5)分析每月的马尔科夫信息量幅值及波动情况，初步判断滑坡下月的形变状态情况；
将各月的马尔科夫信息量累加得到累计马尔科夫信息量曲线，分析累计马尔科夫信息量曲线的整体变化趋势和阶段，据此判断滑坡稳定性时变规律；
对累计马尔科夫信息量曲线进行平滑处理，并进行历史最大值分析，据此判断滑坡稳定性时变规律的关键时间节点信息，预测滑坡稳定性发展趋势。
进一步，所述步骤(2)中，利用K均值算法对每月的二维特征数据集进行 聚类分析包括以下步骤：
预设聚类数为K；
在所述二维特征数据集中随机选取K个二维特征数据分别设为K个初始聚类中心；
计算除K个初始聚类中心外其他二维特征数据分别到K个初始聚类中心的距离，分别比较每个二维特征数据到K个初始聚类中心的距离，将每个二维特征数据和与其距离最小的初始聚类中心归为同一类簇，初次聚类后，计算每个类簇里所有二维特征数据的平均值，并将该平均值确定为新的聚类中心，重复计算每个二维特征数据分别到每个新的聚类中心的距离，并比较距离，重新得到类簇，依此重复直至损失函数的精度小于0.001或相邻损失函数的值不变，所述损失函数的值越小，各个类簇的区分度越大。
进一步，所述损失函数的计算方式为：
式中：L表示损失函数；Dintra表示一类簇里所有二维特征数据到对应聚类中心的距离总和，反映类内数据的一致性；Dinter表示K个聚类中心之间的距离总和，反映类间数据的差异性；xi表示二维特征数据i；cj为第j类簇的聚类中心。
进一步，利用K均值算法对每月的二维特征数据集重复进行N次聚类，N大于1，将N次聚类所得的类簇根据择多判决法确定最终的类簇。
进一步，所述步骤(3)中，对步骤(2)得到的每个类簇分别建立云模型，计算每个云模型的数字特征，所述数字特征包括均值、熵和超熵，然后计算并 比较所有二维特征数据相对于每个云模型数字特征的隶属度，将每个二维特征数据归入隶属度最高的云模型中。
进一步，所述数字特征通过逆向云发生器计算，所述云模型通过正向云发生器进行直观分析。
进一步，所述步骤(4)中，状态出现信息量计算公式如下：
式中：at为第t月滑坡的变形加速度，Si为本月滑坡的形变状态，Ni为滑坡形变状态Si出现的频数，pi为Si出现的概率，Ii为状态出现信息量。
进一步，所述步骤(4)中，状态转移信息量计算公式如下：
式中：Iij表示从Si转移到Sj的状态转移信息量，pij表示从Si转移到Sj出现的概率。
进一步，所述步骤(4)中，每月的马尔科夫信息量的计算公式如下：
式中：uij为马尔科夫信息量，且uij＝Ij*Iij*mbst，mbst表示第t月云模型中二维特征数据到对应云模型的隶属度。
与现有技术相比，本发明有效的运用马尔科夫信息量来表现滑坡演化变形过程中稳定性的时变规律，当滑坡进入加速变形阶段后，累计马尔科夫信息量将出现明显的波动上升趋势，从而能够实现对于滑坡的早期预警，而且计算简 便，同时，本发明引入的云模型不仅能分析数据的二维特征，以隶属度表达了各月数据与各形变状态典型数字特征的相似程度，采用云模型修正的标签也能够很好地反映监测数据的形变状态。
附图说明
图1为本发明一实施例的流程图。
图2为本发明一实施例的滑坡形变状态的类簇示意图。
图3为本发明一实施例的筛选比例0.6正向云模型生成结果图。
图4为本发明一实施例的K均值所得类簇与云模型修正结果的对比图。
图5为本发明一实施例的每月马尔科夫信息量示意图。
图6为本发明一实施例的累计马尔科夫信息量示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。
实施例
请参照图1，本实施例包括以下步骤：
(1)根据滑坡每月的变形速度和变形加速度构成每月的二维特征数据集；
(2)利用K均值算法对每月的二维特征数据集进行聚类分析得到类簇，并为每个类簇分配聚类标签；
预设聚类数为K；
在所述二维特征数据集中随机选取K个二维特征数据分别设为K个初始聚类中心；
计算除K个初始聚类中心外其他二维特征数据分别到K个初始聚类中心的距离，分别比较每个二维特征数据到K个初始聚类中心的距离，将每个二维特征数据和与其距离最小的初始聚类中心归为同一类簇，初次聚类后，计算每个类簇里所有二维特征数据的平均值，并将该平均值确定为新的聚类中心，重复计算每个二维特征数据分别到每个新的聚类中心的距离，并比较距离，重新得到类簇，依此重复直至损失函数的精度小于0.001或相邻损失函数的值不变，所述损失函数的值越小，各个类簇的区分布越大。
所述损失函数的计算方式为：
式中：L表示损失函数；Dintra表示一类簇里所有二维特征数据到对应聚类中心的距离总和，反映类内数据的一致性；Dinter表示K个聚类中心之间的距离总和，反映类间数据的差异性；xi表示二维特征数据i；cj为第j类簇的聚类中心。
利用K均值算法对每月的二维特征数据集重复进行N次聚类，N＝9，将N次聚类所得的类簇根据择多判决法确定最终的类簇。重复聚类有效减小类簇结果中存在的随机性。
(3)对每个类簇分别建立云模型，并计算二维特征数据集内二维特征数据相对于云模型的隶属度，根据隶属度对步骤(2)的聚类标签进行修正，并将修正后的聚类标签所代表的形变状态作为滑坡该月的形变状态，以此类推得到滑坡每月的形变状态；
对步骤(2)得到的每个类簇分别建立云模型，计算每个云模型的数字特征， 所述数字特征包括均值、熵和超熵，然后计算并比较所有二维特征数据相对于每个云模型数字特征的隶属度，将每个二维特征数据归入隶属度最高的云模型中。
所述数字特征通过逆向云发生器计算，所述云模型通过正向云发生器进行直观分析。
云模型是一种用于处理定性概念与定量描述的不确定转换模型。设U是一个用精确数值表示的定量论域，C是U上的定性概念，若定量值x∈U，且x是定性概念C的一次随机实现，x对C的确定度μ(x)∈[0,1]是具有稳定倾向的随机数
则x在论域U上的分布称为云，每一个x称为一个云滴，这里的随机实现是概率意义下的实现；确定度指模糊集意义下的隶属度，在云模型中，以均值(Ex)、熵(En)和超熵(He)描述云的整体特征，均值是云滴在论域空间分布的期望，是最能够代表定性概念的点；熵描述定性概念的不确定度，熵越大，概念越宏观；超熵是熵的熵，由熵的随机性和模糊性决定。
(4)根据滑坡每月的形变状态对滑坡的形变状态进行统计，并据此计算滑坡每月形变状态的信息量模型，所述信息量模型包括状态出现信息量和状态转移信息量，将状态出现信息量、状态转移信息量和云模型中二维特征数据到对应云模型的隶属度乘积定义为马尔科夫信息量，据此计算每月的马尔科夫信息量；
状态出现信息量计算公式如下：
式中：at为第t月滑坡的变形加速度，Si为本月滑坡的形变状态，Ni为滑坡形变状态Si出现的频数，pi为Si出现的概率，Ii为状态出现信息量。
状态转移信息量计算公式如下：
式中：Iij表示从Si转移到Sj的状态转移信息量，pij表示从Si转移到Sj出现的概率。
每月的马尔科夫信息量的计算公式如下：
式中：uij为马尔科夫信息量，且uij＝Ij*Iij*mbst，mbst表示第t月云模型中二维特征数据到对应云模型的隶属度。
滑坡危险度时变规律研究至少应包含两个方面：1)表征滑坡当前危险程度的指标；2)表征滑坡危险度变化趋势的指标。为了表征滑坡当前危险程度，从各形变状态出现概率的差异出发，构造状态出现信息量。为了表征滑坡稳定性变化趋势，充分考虑外界因素对滑坡变形的作用机制，引入马尔科夫链模型对形变状态跃迁规律进行分析，构造状态转移信息量。同时考虑到各月实际数据对形变状态不同的隶属度(云模型的输出结果)，最终将状态出现信息量、状态转移信息量和隶属度进行乘积运算，定义为马尔科夫信息量。
具体地，香农指出，信息量的大小与状态出现概率有关，且满足两个条件：①信息量是概率p的单调递降函数，即概率越大，信息量越小；②两个独立状态共同产生的信息量应等于其各自信息量的总和。同时满足这两个条件的函数即对数函数。设离散型随机变量X，p(x)表示X取值为x时的概率，则状态x的信 息量I(x)和信源X的信息熵H(X)如下所示：
I(x)＝-logp(x)
运用数理统计得到各形变状态的频数和概率，并计算各形变状态的信息量I(x)。考虑到任意形变状态的信息量I(x)是其出现多次的总贡献，因此将信息量I(x)除以对应的频数得到形变状态的平均信息量，以表达形变状态出现一次携带的信息量。与基础的信息量相比，平均信息量更加突出了小概率形变状态的重要性。然后对各形变状态的平均信息量进行归一化处理。为了使信息量能够反映滑坡的变形趋势，将变形加速度的符号作为信息量的符号，处理后的信息量定义为状态出现信息量，如下所示：
滑坡当月的形变状态将受到且仅受到上一个月形变状态影响，而与更早的形变状态无关。这种现象表明，滑坡形变状态是一个典型的马尔科夫链模型滑坡形变状态是一个典型的马尔科夫链模型。
马尔科夫链被提出用于描述这样一类状态序列，其每个状态只与前一个状态相关，而与其他状态无关，数学表达式如下所示：
P(xn＝in|xn-1＝in-1,…,x1＝i1)＝P(xn＝in|xn-1＝in-1)
式中：n为系统所处的时间，S为系统所有状态的集合，称为状态空间，in表示系统在n时刻的状态，且in∈S,n＝1,2,…N。
在马尔科夫链模型中，状态转移概率矩阵是描述状态跃迁规律的重要工具，其基本形式如下所示：
式中：第一列表示初始状态Si；第一行表示转移后的状态Sj；pij表示从状态Si经过一个时刻转移到状态Sj的概率，即pij＝P(xn＝Sj|xn-1＝Si)；N为总的状态数。
状态转移概率矩阵包含着滑坡形变状态跃迁的重要规律。一方面，严重变形出现的概率比较低，导致由其他形变状态转移到严重形变状态的概率较低。另一方面，严重形变状态表明滑坡具有较高的不稳定性，因此具有可持续性差的特点。显然，严重变形的持续时间越长，概率越低，滑坡的危险性越高，其对滑坡灾害防治的意义越大。因此，首先统计获得形变状态的状态转移矩阵，然后运用信息量描述各种转移情况的重要性。参照将归一化的加权信息量定义为状态转移信息量，如下所示：
式中：Iij表示从Si转移到Sj的状态转移信息量。
状态出现信息量和状态转移信息量从两个角度对滑坡变形规律进行了表达，虽然侧重点不同，但均从概率角度出发，以信息量作为度量。同时，信息量数值的增加均表明滑坡危险性的上升，具有一致的物理意义。同时考虑到相同标签的月份同样存在数据差异，因此引入云模型中的隶属度以示区分。最终，将状态出现信息量、状态转移信息量和隶属度的乘积定义为马尔科夫信息量，如下所示：
给定滑坡形变状态的序列，根据本月和上月的形变状态可直接读取本月的马尔科夫信息量。例如，若上月的形变状态为S1，本月的形变状态为S2，那么本月的马尔科夫信息量即为u12。
(5)分析每月的马尔科夫信息量幅值及波动情况，初步判断滑坡下月的形变状态情况；
将各月的马尔科夫信息量累加得到累计马尔科夫信息量曲线，分析累计马尔科夫信息量曲线的整体变化趋势和阶段，据此判断滑坡稳定性时变规律；
对累计马尔科夫信息量曲线进行平滑处理，并进行历史最大值分析，据此判断滑坡稳定性时变规律的关键时间节点信息，预测滑坡稳定性发展趋势。
根据每月马尔科夫信息量的变化情况初步判断滑坡稳定性变化情况，即如果在一段监测期内，每月马尔科夫信息量维持在较低水平，且存在正负交替情况，说明该滑坡在这一监测期内基本稳定，危险度较低，如果滑坡所处监测期内，与前期极值相当水平的数值频繁出现，则表明滑坡状态不稳定，危险度较高；对每月马尔科夫信息量进行累加得到累积马尔科夫信息量。
累计马尔科夫信息量曲线能够体现滑坡不稳定性的累积效应，以及滑坡稳定性在整个监测期内的时变规律。考虑到滑坡变形是一个渐变的过程，在发生破坏之前，滑坡危险度不可能发生突变，因此需要先对累计马尔科夫信息量曲线进行平滑处理。分析累计马尔科夫信息量时，历史最大值，即从开始监测到当前月份的最大值，是表明滑坡最危险状态的重要指标。每一次历史最大值的 更新，都表明滑坡变得更危险。当历史最大值频繁更新时，累计马尔科夫信息量曲线必将出现明显上升趋势，表明滑坡处于高度不稳定状态。此时应加强滑坡监测和灾害预警工作。
本方法验证如下：以三峡库区新滩滑坡为例，对三峡库区新滩滑坡稳定性时变规律进行分析。
首先获取新滩滑坡1978年1月到1985年5月的累计位移数据和变形速度，并计算变形加速度，变形速度体现了滑坡当前的变形程度，变形加速度则体现了变形趋势，因此，根据新滩滑坡每个月的变形速度和变形加速度构成每月的二维特征数据集(如图2所示)。
利用K均值算法对87个月的二维特征数据集进行聚类得到类簇，为所述类簇分配聚类标签；
预设聚类数为3，分别为聚类标签1、2、3；
利用K均值算法对每月的二维特征数据集重复进行9次聚类，将9次聚类所得的类簇根据择多判决法确定最终的类簇，所得结果如表1所示。重复聚类一定程度减小类簇中存在的随机性。
表1 二维特征数据集的类簇
由表1可以看出，三个聚类中心具有明显的数值差异，能够反映滑坡不同的变形程度。大部分数据分布在1类中，而1类的变形速度和变形加速度均处于较低的水平，能够证明滑坡变形过程中，轻微变形出现的概率高，严重变形出现概率低这一事实。
由于聚类分析属于无监督方法，所得类簇依然存在随机性，为了处理这种随机性和模糊性，本发明引入云模型进行分析。
利用云模型对聚类标签进行修正，根据变形速度和变形就速度的物理意义，在变形速度维度上将1类设置为右半云，2类为对称云，3类为左半云，在变形加速度维度上均设置为对称云，在筛选比例为0.6的条件下正向云生成结果可参见图3。图3中，三角形为筛选后的原始数据点，其余为3G生成的云滴。可以看出，云模型不仅考虑了数据的均值，同时还考虑了熵(方差)等属性，并以隶属度包络表达数据的二维特性。
K均值所得类簇与云模型修正结果的对比可参见图4。图4中，圆圈表示第1类，五角星表示第2类，正方形表示第3类。标记较大的表示K均值的聚类标签，标记较小的表示云模型的标签。可以看出，K均值的聚类标签与云模型的标签大部分一致，但存在1类被2类和3类侵占的情况，主要是因为1类数据密集，对数值偏差容忍度较小。
将所述云模型的标签作为修正后的聚类标签，进而代入信息量模型进行分析。统计各形变状态的出现频次，并计算得到各形变状态的状态出现信息量的绝对值(正负号由每月的加速度实际值决定)，分别为：1类0.1621，2类0.4980，3类0.3399：
最后计算得出的状态转移信息量矩阵如表2所示：
表2 状态转移信息量
由表2可以看出，S2跃迁到S1时的转移信息量为0。根据各类的数字特点易知，S2具有一定的正变形速度和变形加速度，S1的变形速度和变形加速度均 很小。当滑坡从S2跃迁到S1时，必定出现负的加速度，且变形速度介于S2变形速度和S1变形速度之间，这样的数字特征与S3的更为吻合，因此S2无法直接跃迁到S1。其次转移信息量最小的为S1到S1，即维持状态，这主要是因为S1代表轻微变形。转移信息量最大的出现在S2向S2的转移，表明了滑坡的变形速度在持续增加，滑坡的危险度也不断增加，符合实际意义。
最终计算得出各个月份的马尔科夫信息量，可参见图5。分析图5可以看出，前期月马尔科夫信息量维持在较低水平，且存在正负交替情况，说明新滩滑坡在监测期前段基本稳定，危险度较低。存在局部极值，即短时间的快速变形，但持续时间很短不长。在监测期后期，与前期极值相当水平的数值频繁出现，表明滑坡状态不稳定，危险度较高。
对月马尔科夫信息量进行累加得到累积马尔科夫信息量。同时对累计马尔科夫信息量进行均值平滑，得到累计马尔科夫信息量曲线，可参见图6。
图6中，柱状图为累计马尔科夫信息量，柱状图上的曲线为均值平滑结果，柱状图上的圆点为历史最大值。对图6进行历史最大值分析可以看出，新滩滑坡的累计马尔科夫信息量主要具有两种变化形式：零点波动型和波动上升型。1978年12月到1982年2月，累计马尔科夫信息量基本在零点上下波动，整体趋势略有降低。1979年8月(1号点)出现局部极值点，此时累计量约为0.32。1982年2月(2号点)，累计马尔科夫信息量出现监测期内最小值，约为-0.17。1982年2月之后，累计马尔科夫信息量开始具有明显波动上升趋势。由历史最大值点的分析可以看出，历史最大值同样主要存在两个阶段。1977年12月到1979年7月前后为第一阶段。该阶段处于监测初期，历史最大值处于较低水平，因此容易被更新。1982年6月前后到1985年4月为第二阶段，该阶段处于监测 期的后期，此时历史最大值被频繁更新表明滑坡的危险度在不断增加，需要对加强滑坡灾害预警。
根据已有新滩滑坡演化阶段的研究，新滩滑坡在1979年8月进入匀速变形阶段，于1982年7月进入加速变形阶段，并于1985年6月发生破坏。结合图6可以发现，1979年8月(1号点)对应于累计马尔科夫信息量前期的局部最大值，同时处于历史最大值第一阶段的末端，在其两侧历史最大值情况具有明显差异。1982年7月(2号点)位于累计马尔科夫信息量的波动上升型期间，同时处于历史最大值分析中第二阶段的起始点附近，是滑坡危险度不断上升过程的起点。1985年6月滑坡发生破坏前，累计马尔科夫信息量已经达到相当高的数值水平，且为新的历史最大值。
因此，本发明所述的马尔科夫信息量方法能够反映滑坡稳定性的变化情况。尤其在关键时间点上，与实际工程情况高度吻合，说明了该方法的有效性，能够反映滑坡稳定性时变规律，实现对于滑坡的早期预警。
- 还没有人留言评论。精彩留言会获得点赞！