CN109783889A - 基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法 - Google Patents

Abstract

本发明公开的基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法，首先通过安装在坡体上的激光位移传感器采集滑坡体的位移数据，并对采集的位移数据在时序方向上进行多状态划分；然后，将划分后的位移数据通过Baum‑welch算法进行训练，获得多状态MOG‑HMM模型，并对已划分的状态进行标记；最后，采用Viterbi算法对实时采集的位移数据进行当前状态估计，并将当前估计状态作为Dijkstra算法的输入，从而预测出滑坡的发生时间。本发明公开的方法根据滑坡演化过程中的位移数据，推演出当前滑坡演化状态，并计算出滑坡灾害发生的时间，提升了一般HMM方法在预测中的鲁棒性性能，提高了地质灾害预报的准确性。

Description

基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法

技术领域

本发明属于地质灾害预报方法技术领域，具体涉及一种基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法。

背景技术

滑坡不仅是一种自然灾害，同时也是一种严重的工程地质灾害。近几年来重大地质灾害频繁发生，常造成房屋破坏、通讯设施中断，道路崩塌、毁坏土地，乃至村毁人亡的事故。且事故常发生在地质结构复杂、地形陡峭的山区，对于预防的实施、灾后恢复的时间将被延长，造成居民生活不便，其规模和风险性大大超过了我们所能承受的范围。因此如何使用技术手段进行监测预报，进行实时提示受灾群众，减少灾害损失，成为我们关注的主要内容。

目前，关于滑坡预报主要集中在以下三个方面，第一、滑坡位移预报，中国矿业大学的高彩霞等提出了四种滑坡位移预测方法，即BP、RBF、新型智能算法和智能耦合模型，这四种方法均有较准确的位移预测精度，特别是对1步预测尤其有效；第二、在位移预测的基础上，并结合其他相关参数，进行滑坡稳定性预报，西南科技大学的蔡嘉伦等针对目前滑坡稳定性评价研究的不足，采用卫星实时监测的方法，利用表面位移规律对滑坡体进行稳定性分析，此种方法具有监测方法稳定直观、分析结果实时可靠的优势，三峡大学的胡安龙等基于相关系数理论描述滑坡稳定性影响因素对滑坡稳定状态影响的大小，再根据关联度大小筛选出影响滑坡稳定性的主要影响因素，其次，引入贝叶斯理论，滑坡稳定性的主要影响因素和滑坡的稳定状态建立了基于贝叶斯的滑坡稳定性预测模型；第三、滑坡发生时间预报，位移预测和滑坡稳定性预测能够对预估滑坡发生时间提供重要的参考作用，但如果需要准确计算出滑坡发生具体时间，还需要在位移预测和稳定性预测的基础上进一步分析，中国科学院力学研究所的王建锋利用Pearl曲线，结合滑坡位移变化数据，从运动学角度，判断出滑坡爆发的峰值点，为滑坡时间预报开辟出一条新的道路。现有方法存在的问题是虽然能够估计滑坡发生的时间点，但未能最终计算出多久后会发生滑坡灾害，不能真正称为滑坡发生时间预报。

从以上分析可以看出，滑坡灾害的预报呈现出越来越准确的趋势，如果滑坡灾害的时间预报能够像天气预报一样准确，那么势必会为防灾减灾工作提供更大的帮助。本发明建立了以多状态混合高斯隐马尔科夫模型(MixtureofGaussian-Hidden markov Model)MOG-HMM为模型的滑坡发生时间预测算法，能够根据滑坡演化过程中的位移数据，推演出当前滑坡演化状态，并进一步计算出生滑坡灾害发生的时间。

发明内容

本发明的目的是提供一种基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法，解决了现有滑坡发生预报算法准确率偏低的问题。

本发明所采用的技术方案是，基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法，具体过程包括如下步骤：

步骤1，通过安装在坡体上的激光位移传感器采集滑坡体的位移数据，并对采集的位移数据在时序方向上进行多状态划分；

步骤2，将划分后的位移数据通过Baum-welch算法进行训练，获得多状态MOG-HMM模型，并对已划分的状态进行标记；

步骤3，采用Viterbi算法对实时采集的位移数据进行当前状态估计，并将当前估计状态作为Dijkstra算法的输入，从而预测出滑坡的发生时间。

本发明的其他特点还在于，

步骤1中对采集的位移数据在时序方向上进行多状态划分原则为：根据滑坡的演化阶段细分。

步骤2中将划分后的位移数据通过Baum-welch算法进行训练并获得MOG-HMM模型的具体过程如下：

步骤2.1通过Baum-welch算法计算MOG-HMM模型λ的6个参数，如公式(1)所示：

λ＝(P,A,B,μ(D(S_t)),σ(D(S_t)),S_final) (1)

其中，λ为基于MOG-HMM的滑坡预报模型，等式右边为该模型的6个相关参数，P：初始状态分布；A：状态转移概率矩阵；B：表示观察值概率分布矩阵；μ(D(S_t))：MOG函数中的均值，表征每个状态持续的平均时间，其中，D(.)表示访问时间，S_t表示状态序列；σ(D(S_t))：MOG函数中的方差；S_final：最终状态；

步骤2.2训练MOG-HMM模型。

步骤2.2中训练MOG-HMM模型的具体过程如下：

(1)记录坡体位移原始数据X_i，i＝1，2，...，n；

(2)初始化X_i的多状态个数N，其中，T为整个滑坡演化过程持续时间，t为滑坡灾害发生状态持续时间；设置每个状态的高斯混合数M；

(3)初始化P，A和B，将观测序列X_i平均分为N份，每份称为一个状态；采用聚类方法将X_i的每个状态聚类为M类，并将每类按照高斯混合序号(1,2,……,M)进行编号，得到时间序列O_t(t＝1，2，...，N)；

(4)训练参数P和A，通过Baum-welch算法求解使期望P(O_t|λ)最大时对应的P和A；参数B由滑坡位移时序通过混合高斯分布参数估计如公式(2)和公式(3)求得：

其中，b_j(O_t)是B中的具体元素值，j为B中元素的索引值，N为B中元素的个数，M为加权系数个数，C_j，m为加权系数，ξ(O，μ_j，m，σ_j，m)为混合高斯函数；

(5)然后利用Viterbi方法解码，得到整个时序O_t对应的状态序列S_t；

(6)通过公式(4)～(6)计算中间参数均值μ(D(S_t))、方差σ(D(S_t))和最终状态S_final；

其均值和方差计算过程如公式(4)和(5)所示。

S_final＝S_N (6)

其中，μ(D(S_t))：MOG函数中的均值，表征每个状态持续的平均时间，S_t表示状态序列；D(S_t)表示第t个状态S_t的访问时间，D(S_t,l)表示第l次访问状态S_t的访问时间，S_final：最终状态，S_N表示滑坡序列划分出的第N个状态。

步骤3的具体过程如下：

步骤3.1采用Viterbi算法对实时采集的数据进行状态解码，从而判断当前滑坡所处的演化状态，即估计出测试数据属于滑坡演化过程中的哪个阶段，具体过程如下：

(1)输入：观测时间序列O_t和滑坡发生时间预报模型λ；

输出：最短路径I^*＝(i₁ ^*，i₂ ^*，...，i_T ^*)；

(2)t＝1时，δ₁(i)＝P_ib_i(o₁)，i＝1，2，…，N；ψ₁(i)＝0；

其中，δ₁(i)和ψ₁(i)为计算的中间变量，δ₁(i)由P和B的第i个元素相乘得到，其中O₁为O_t中的第1个元素；

(3)t＝2，3，…，T时，

ψ_t(i)＝argmax[δ_t-1(j)a_j，i]；

其中，δ_t(i)为上一时刻所有δ_t-1(j)与A中元素a_j,i和B中元素b_i(O_t)乘积的最大值，其中，1≤j≤N，ψ_t(i)为使δ_t-1(j)a_j,i最大对应的索引；

(4)终止

其中C*为所有δ_t(i)中元素的最大值，i_T ^*为使δ_t(i)最大对应的索引；

(5)最优路径回溯：对t＝T-1，T-2，…，1，利用i_t ^*＝ψ_t+1(i_i+1 ^*)，i_t ^*为I^*中的具体元素，求得S_t＝(S₁，S₂，…，S_T)；判断出当前状态后，截取新采集数据对应状态序列S_t的最后一段序列，如公式(7)所示：

S_end＝(S_t-l，…，S_t-2，S_t-1，S_T) (7)

其中，l为最后一段序列的截取位置，T为当前时刻；

步骤3.2将S_end中数量最多的值作为当前状态S_current，对当前状态到达滑坡发生状态的最短时间RUL进行估计：

根据当前状态S_current和滑坡发生时间预报模型λ得出的状态S_final，将当前估计状态作为Dijkstra算法的输入，从而得出从S_current到达S_final最短时间，如公式(8)所示：

本发明的有益效果是，建立了以多状态MOG-HMM为模型的滑坡发生时间预测算法，能够根据滑坡演化过程中的位移数据，推演出当前滑坡演化状态，并进一步计算出从当前时刻开始多久后会发生滑坡灾害。本发明最大程度上提升了一般HMM方法在预测中的鲁棒性性能，提高了地质灾害预报的准确性，为进一步更及时有效的预报滑坡发生奠定理论基础。

附图说明

图1是本发明的基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法流程图；

图2是实施例中测试数据进行多状态划分的曲线图；

图3是实施例中测试数据的测试结果；

图4是实施例中测试数据的测试结果的绝对误差。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明的基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法，具体过程包括如下步骤：

步骤1，通过安装在坡体上的激光位移传感器采集滑坡体的位移数据，并对采集的位移数据在时序方向上进行多状态划分；

步骤1中对采集的位移数据在时序方向上进行多状态划分原则为：根据滑坡的演化阶段细分；

步骤2，将划分后的位移数据通过Baum-welch算法进行训练，获得多状态MOG-HMM模型，并对已划分的状态进行标记；

步骤2中将划分后的位移数据通过Baum-welch算法进行训练并获得MOG-HMM模型的具体过程如下：

步骤2.1通过Baum-welch算法计算MOG-HMM模型λ的6个参数，如公式(1)所示：

λ＝(P，A，B，μ(D(S_t))，σ(D(S_t))，S_final) (1)

其中，λ为基于MOG-HMM的滑坡预报模型，等式右边为该模型的6个相关参数，P：初始状态分布；A：状态转移概率矩阵；B：表示观察值概率分布矩阵；μ(D(S_t))：MOG函数中的均值，表征每个状态持续的平均时间，其中，D(.)表示访问时间，S_t表示状态序列；σ(D(S_t))：MOG函数中的方差；S_final：最终状态；

步骤2.2训练MOG-HMM模型。

步骤2.2中训练MOG-HMM模型的具体过程如下：

(1)记录坡体位移原始数据X_i，i＝1，2，...，n；

(2)初始化X_i的多状态个数N，其中，T为整个滑坡演化过程持续时间，t为滑坡灾害发生状态持续时间；设置每个状态的高斯混合数M；

(3)初始化P，A和B，将观测序列X_i平均分为N份，每份称为一个状态；采用聚类方法将X_i的每个状态聚类为M类，并将每类按照高斯混合序号(1,2,……,M)进行编号，得到时间序列O_t(t＝1，2，...，N)；

(4)训练参数P和A，通过Baum-welch算法求解使期望P(O_t|λ)最大时对应的P和A；参数B由滑坡位移时序通过混合高斯分布参数估计如公式(2)和公式(3)求得：

其中，b_j(O_t)是B中的具体元素值，j为B中元素的索引值，N为B中元素的个数，M为加权系数个数，C_j，m为加权系数，ξ(O，μ_j，m，σ_j，m)为混合高斯函数；

(5)然后利用Viterbi方法解码，得到整个时序O_t对应的状态序列S_t；

(6)通过公式(4)～(6)计算中间参数均值μ(D(S_t))、方差σ(D(S_t))和最终状态S_final；

其均值和方差计算过程如公式(4)和(5)所示。

S_final＝S_N (6)

其中，μ(D(S_t))：MOG函数中的均值，表征每个状态持续的平均时间，S_t表示状态序列，D(S_t)表示第t个状态S_t的访问时间，D(S_t,l)表示第l次访问状态S_t的访问时间，S_final：最终状态，S_N表示滑坡序列划分出的第N个状态。

步骤3，采用Viterbi算法对实时采集的位移数据进行当前状态估计，并将当前估计状态作为Dijkstra算法的输入，从而预测出滑坡的发生时间。

步骤3的具体过程如下：

步骤3.1采用Viterbi算法对实时采集的数据进行状态解码，从而判断当前滑坡所处的演化状态，即估计出测试数据属于滑坡演化过程中的哪个阶段，具体过程如下：

(1)输入：观测时间序列O_t和滑坡发生时间预报模型λ；

输出：最短路径I^*＝(i₁ ^*,i₂ ^*,...,i_T ^*)，；

(2)t＝1时，δ₁(i)＝P_ib_i(o₁)，i＝1，2，…，N；ψ₁(i)＝0；

其中，δ₁(i)和ψ₁(i)为计算的中间变量，δ₁(i)由P和B的第i个元素相乘得到，其中O₁为O_t中的第1个元素。

(3)t＝2，3，…，T时，

ψ_t(i)＝argmax[δ_t-1(j)a_j，i]；

其中，δ_t(i)为上一时刻所有δ_t-1(j)与A中元素a_j,i和B中元素b_i(O_t)乘积的最大值，其中，1≤j≤N，ψ_t(i)为使δ_t-1(j)a_j,i最大对应的索引；

(4)终止：

其中C*为所有δ_t(i)中元素的最大值，i_T ^*为使δ_t(i)最大对应的索引；

(5)最优路径回溯：对t＝T-1，T-2，…，1，利用i_t ^*＝ψ_t+1(i_i+1 ^*)，i_t ^*为I^*中的具体元素，求得S_t＝(S₁，S₂，…，S_T)；判断出当前状态后，截取新采集数据对应状态序列S_t的最后一段序列，如公式(7)所示：

S_end＝(S_t-l，…，S_t-2，S_t-1，S_T) (7)

其中，l为最后一段序列的截取位置，T为当前时刻；

步骤3.2将S_end中数量最多的值作为当前状态S_current，对当前状态到达滑坡发生状态的最短时间(Remaining useful life,RUL)进行估计：

根据当前状态S_current和滑坡发生时间预报模型λ得出的状态S_final，将当前估计状态作为Dijkstra算法的输入，从而得出从S_current到达S_final最短时间，如公式(8)所示：

实施例

测试数据：卧龙寺滑坡5#监测点

卧龙寺滑坡位于我国西部黄土地区，属于典型的塬边黄土滑坡，由于频繁的地震活动和地壳的振荡性上升，为滑坡的形成和分布提供了有利的地形、岩性条件。

预测过程：

(1)多状态划分

将数据测试数据1的所有时序划分为6个状态，如图2所示；

(2)预报模型训练

对以上划分过状态的滑坡位移数据用Baum-welch算法进行训练。设置状态数N为8，每个状态中的高斯混合数M为2。首先随机初始化P和A以及B，由于本文中采用的HMM结构为左右型，所以P不需要更新，经过Baum-welch算法训练后得到高斯混合M、P以及A的值。

(3)RUL计算

对测试数据进行测试，测试结果如图3所示；从该图可以看出，本文方法能够比较准确的预报滑坡发生时间，如第30天时，计算得到的RUL为29天，误差较小。

(4)预测结果评价

直观评价预测的一个指标是绝对误差，如图4所示；

从图4的结果可以分析得出该预测是有效的，为了进一步评价该预测的性能，根据较权威的预测评价指标对预测结果进行了评价。表1为利用预测标准中的评价指标对上述滑坡发生时间预测结果进行评价的结果，如果计算出来的各个预测评价指标都在范围内，说明该预测算法针对滑坡发生时间是可行且有效的。

表1基于MOG-HMM的滑坡发生时间预测算法的预测效果评价结果

通过表1的结果可以看出，基于测试数据的预测结果，三项指标的结果均在范围内，说明本发明方法对于解决HMM用于滑坡灾害发生时间鲁棒性差的问题是有效的。

Claims (5)

1.基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法，其特征在于，具体过程包括如下步骤：

步骤1，通过安装在坡体上的激光位移传感器采集滑坡体的位移数据，并对采集的位移数据在时序方向上进行多状态划分；

步骤2，将划分后的位移数据通过Baum-welch算法进行训练，获得多状态MOG-HMM模型，并对已划分的状态进行标记；

步骤3，采用Viterbi算法对实时采集的位移数据进行当前状态估计，并将当前估计状态作为Dijkstra算法的输入，从而预测出滑坡的发生时间。

2.如权利要求1所述的基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法，其特征在于，所述步骤1中对采集的位移数据在时序方向上进行多状态划分原则为：根据滑坡的演化阶段细分。

3.如权利要求1所述的基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法，其特征在于，所述步骤2中将划分后的位移数据通过Baum-welch算法进行训练并获得MOG-HMM模型的具体过程如下：

步骤2.1通过Baum-welch算法计算MOG-HMM模型λ的6个参数，如公式(1)所示：

λ＝(P，A，B，μ(D(S_t))，σ(D(S_t))，S_final) (1)

其中，λ为基于MOG-HMM的滑坡预报模型，等式右边为该模型的6个相关参数，P：初始状态分布；A：状态转移概率矩阵；B：表示观察值概率分布矩阵；μ(D(S_t))：MOG函数中的均值，表征每个状态持续的平均时间，其中，D(.)表示访问时间，S_t表示状态序列；σ(D(S_t))：MOG函数中的方差；S_final：最终状态；

步骤2.2训练MOG-HMM模型。

4.如权利要求3所述的基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法，其特征在于，所述步骤2.2中训练MOG-HMM模型的具体过程如下：

(1)记录坡体位移原始数据X_i，i＝1，2，...，n；

(2)初始化X_i的多状态个数N，其中，T为整个滑坡演化过程持续时间，t为滑坡灾害发生状态持续时间；设置每个状态的高斯混合数M；

(3)初始化P，A和B，将观测序列X_i平均分为N份，每份称为一个状态；采用聚类方法将X_i的每个状态聚类为M类，并将每类按照高斯混合序号(1,2,……,M)进行编号，得到时间序列O_t(t＝1，2，...，N)；

(4)训练参数P和A，通过Baum-welch算法求解使期望P(O_t|λ)最大时对应的P和A；参数B由滑坡位移时序通过混合高斯分布参数估计如公式(2)和公式(3)求得：

其中，b_j(O_t)是B中的具体元素值，j为B中元素的索引值，N为B中元素的个数，M为加权系数个数，C_j，m为加权系数，ξ(O，μ_j，m，σ_j，m)为混合高斯函数；

(5)然后利用Viterbi方法解码，得到整个时序O_t对应的状态序列S_t；

(6)通过公式(4)～(6)计算中间参数均值μ(D(S_t))、方差σ(D(S_t))和最终状态S_final；

其均值和方差计算过程如公式(4)和(5)所示。

S_final＝S_N (6)

其中，μ(D(S_t))：MOG函数中的均值，表征每个状态持续的平均时间，S_t表示状态序列；D(S_t)表示第t个状态S_t的访问时间，D(S_t,l)表示第l次访问状态S_t的访问时间，S_final：最终状态，S_N表示滑坡序列划分出的第N个状态。

5.如权利要求4所述的基于混合高斯隐马尔科夫模型的滑坡发生时间预测方法，其特征在于，所述步骤3的具体过程如下：

步骤3.1采用Viterbi算法对实时采集的数据进行状态解码，从而判断当前滑坡所处的演化状态，即估计出测试数据属于滑坡演化过程中的哪个阶段，具体过程如下：

(1)输入：观测时间序列O_t和滑坡发生时间预报模型λ；

输出：最短路径I^*＝(i₁ ^*，i₂ ^*，...，i_T ^*)；

(2)t＝1时，δ₁(i)＝P_ib_i(o₁)，i＝1，2，…，N；ψ₁(i)＝0；

其中，δ₁(i)和ψ₁(i)为计算的中间变量，δ₁(i)由P和B的第i个元素相乘得到，其中O₁为O_t中的第1个元素；

(3)t＝2，3，…，T时，

ψ_t(i)＝argmax[δ_t-1(j)a_j，i]；

其中，δ_t(i)为上一时刻所有δ_t-1(j)与A中元素a_j,i和B中元素b_i(O_t)乘积的最大值，其中，1≤j≤N，ψ_t(i)为使δ_t-1(j)a_j,i最大对应的索引；

(4)终止

其中C^*为所有δ_t(i)中元素的最大值，i_T ^*为使δ_t(i)最大对应的索引；

(5)最优路径回溯：对t＝T-1，T-2，…，1，利用i_t ^*＝ψ_t+1(i_i+1 ^*)，i_t ^*为I^*中的具体元素，求得S_t＝(S₁，S₂，…，S_T)；判断出当前状态后，截取新采集数据对应状态序列S_t的最后一段序列，如公式(7)所示：

S_end＝(S_t-l，…，S_t-2，S_t-1，S_T) (7)

其中，l为最后一段序列的截取位置，T为当前时刻；

步骤3.2将S_end中数量最多的值作为当前状态S_current，对当前状态到达滑坡发生状态的最短时间RUL进行估计：

根据当前状态S_current和滑坡发生时间预报模型λ得出的状态S_final，将当前估计状态作为Dijkstra算法的输入，从而得出从S_current到达S_final最短时间，如公式(8)所示：