一种基于均值‑标准差的加权马尔可夫滑坡预警方法，属于滑坡灾害预报预警领域；该方法包括：收集滑坡监测点位移速度，构成位移速度集合；利用均值‑标准差法确定位移速度集合中各位移速度所属状态，并构建位移速度状态序列；使用加权马尔可夫链预测位移速度次日预测状态；进行滑坡预警判断；更新位移速度集合；对滑坡预警结果进行评价；本发明从滑坡演变过程中获得的位移参数的随机特性出发，基于马尔可夫链理论，综合运用滑坡演化三阶段理论和均值‑标准差分类法，超越了以往工程界对滑坡预警的认识；给出了具体预警过程，可用于从位移监测数据角度出发的滑坡预警。
1.一种基于均值-标准差的加权马尔可夫滑坡预警方法，其特征在于，包括如下步骤：步骤1：布置滑坡监测点，针对每一个滑坡监测点，执行步骤2到步骤7；
步骤2：每日收集一次滑坡监测点位移速度，构成位移速度集合{x1，x2，…，xi，…，xn}，其中，xi为第i日收集到的位移速度，i＝{1,2，…，n}；
步骤3：利用均值-标准差法确定位移速度集合中各位移速度所属状态，包括正常状态或异常状态，n个位移速度的状态组成位移速度状态序列；
步骤4：根据位移速度状态序列，使用加权马尔可夫链预测次日位移速度xn+1所属状态，即位移速度次日预测状态；
步骤5：进行滑坡预警判断：当满足位移速度前日状态为异常状态、位移速度当日状态为异常状态且位移速度次日预测状态为异常状态时，发出滑坡预警信号；
步骤6：判断是否有新收集到的位移速度数据，有，则将第一日位移速度从位移速度集合中删除，将第n+1日收集到的位移速度加入位移速度集合中，对新的位移速度集合重复步骤3到步骤5；否则，执行步骤7；
步骤3-3：计算位移速度的最大值max和最小值min；
步骤3-4：利用“D”“max”和“min”依次确定位移速度状态：若位移速度在区间[min，D)内，则位移速度状态为正常状态，用“1”表示；若位移速度在区间[D，max]内，则位移速度状态为异常状态，用“2”表示；n个位移速度状态组成位移速度状态序列。
2.根据权利要求1所述的基于均值-标准差的加权马尔可夫滑坡预警方法，其特征在于，步骤4具体步骤为：步骤4-1：使用马尔可夫链计算m阶位移速度次日预测状态的分布矩阵；
步骤4-1-1：根据位移速度状态序列，统计不同步长的状态转移频数矩阵q(k)：其中，k为步长，k＝{1,2，…，m}，m为阶数，m≤5；q11为位移速度状态序列中从正常状态间隔k-1个位移速度状态转移到下一个正常状态的转移次数；q12为位移速度状态序列中从正常状态间隔k-1个位移速度状态转移到异常状态的转移次数；q21为位移速度状态序列中从异常状态间隔k-1个位移速度状态转移到正常状态的转移次数；q22为位移速度状态序列中从异常状态间隔k-1个位移速度状态转移到下一个异常状态的转移次数；
步骤4-1-2：根据状态转移频数矩阵计算状态转移概率矩阵P(k)：其中，
步骤4-1-4：计算m阶位移速度次日预测状态的分布矩阵B：其中，
步骤4-2：使用加权马尔可夫链计算m阶位移速度次日预测状态：步骤4-2-1：计算位移速度的自相关系数rk：
其中，为位移速度均值，xt+k为第t+k日的位移速度，t＝[1，n-k]；
步骤4-2-2：计算不同步长权重wk，构成m阶权重向量wm＝(w1,w2,…,wk,…,wm)，其中，步骤4-2-3：计算m阶位移速度次日预测状态的概率分布向量C(n+1)：(n+1) m
步骤4-2-4：将向量C(n+1)中最大元素所在列对应的状态作为m阶加权马尔可夫链预测的次日位移速度xn+1所属状态，即位移速度次日预测状态；
步骤4-3：重复步骤4-1到步骤4-2，计算m取不同值时各阶位移速度次日预测状态，如果各阶位移速度次日预测状态均为“异常”，位移速度次日预测状态为“异常”；如果各阶位移速度次日预测状态不全为“异常”，位移速度次日预测状态为“正常”。
3.根据权利要求1所述的基于均值-标准差的加权马尔可夫滑坡预警方法，其特征在于，步骤7所述对滑坡预警结果进行评价的具体方法为：步骤7-1：评价滑坡预警的及时性：滑坡预警及时性以滑坡发生前滑坡前应急准备时间内滑坡预警密集度衡量；滑坡预警密集度越接近1，滑坡预警的及时性越强；滑坡预警密集度为：在滑坡发生之前R天进行的R次滑坡预警判断中，发出滑坡预警信号的次数为R*，则滑坡之前R天内的滑坡预警密集度步骤7-2：评价滑坡预警的抗干扰性：错误预警率越低，滑坡预警的抗干扰性越好；错误预警率为：如果在滑坡应急准备时间之前共做出M次滑坡预警判断，其中有M*次发出了滑坡预警信号，错误预警率为步骤7-3：评价滑坡预警的可信度：滑坡预警的可信度由预测一致率决定，预测一致率越大，滑坡预警的可信度越好；预测一致率为：在滑坡发生前一共对位移速度状态做Q次预*测，其中Q次与实际位移速度状态相符合，则预测一致率为
一种基于均值-标准差的加权马尔可夫滑坡预警方法
技术领域
[0001] 本发明属于滑坡灾害预报预警领域，具体涉及一种基于均值-标准差的加权马尔可夫滑坡预警方法。
背景技术
[0002] 滑坡作为一种频发的自然地质灾害，会造成巨大的生命财产损失，其危害已成为仅次于地震的第二大自然灾害。据初步统计，我国每年由崩塌、滑坡、泥石流等灾害所造成的直接经济损失约200亿元人民币，间接损失更是无法估量。
[0003] 监测预警作为规避滑坡等地质灾害的重要手段之一，越来越受到人们的重视。尽管我国近10年来在地质灾害监测预警方面取得了一定的研究成果，积累了一定的经验和教训。但是，在滑坡预警方面并不理想：一是在成功预警实例中专业预警比例过低；二是预警的滞后性，即预警不够及时；三是错误预警较多，即在并无滑坡危险时给出预警提示，对生产及生活造成一定干扰；此外，滑坡预警包括滑坡时间预警和滑坡区域预警，以往研究中比较侧重滑坡时间预警，而对滑坡区域预警的研究不足，而时间预警与区域预警应该同时实现。
发明内容
[0004] 针对上述现有技术存在的不足，本发明提供一种基于均值-标准差的加权马尔可夫滑坡预警方法。
[0005] 本发明的技术方案：
[0006] 一种基于均值-标准差的加权马尔可夫滑坡预警方法，包括如下步骤：
[0007] 步骤1：布置滑坡监测点，针对每一个滑坡监测点，执行步骤2到步骤7；
[0008] 步骤2：每日收集一次滑坡监测点位移速度，构成位移速度集合{x1，x2，…，xi，…，xn}，其中，xi为第i日收集到的位移速度，i＝{1,2，…，n}；
[0009] 步骤3：利用均值-标准差法确定位移速度集合中各位移速度所属状态，包括正常状态或异常状态，n个位移速度的状态组成位移速度状态序列；
[0010] 步骤3-1：计算位移速度的均值 和标准差s：
[0011]
[0012] 步骤3-2：确定状态划分的分界点D：
[0013]
[0014] 其中，α为可变参数；
[0015] 步骤3-3：计算位移速度的最大值max和最小值min；
[0016] 步骤3-4：利用“D”“max”和“min”依次确定位移速度状态：若位移速度在区间[min，D)内，则位移速度状态为正常状态，用“1”表示；若位移速度在区间[D，max]内，则位移速度状态为异常状态，用“2”表示；n个位移速度状态组成位移速度状态序列；
[0017] 步骤4：根据位移速度状态序列，使用加权马尔可夫链预测次日位移速度xn+1所属状态，即位移速度次日预测状态：
[0018] 步骤4-1：使用马尔可夫链计算位移速度次日预测状态的分布矩阵；
[0019] 步骤4-1-1：根据位移速度状态序列，统计不同步长的状态转移频数矩阵q(k)：
[0020]
[0021] 其中，k为步长，k＝{1,2，…，m}，m为阶数，m≤5；q11为位移速度状态序列中从正常状态间隔k-1个位移速度状态转移到下一个正常状态的转移次数；q12为位移速度状态序列中从正常状态间隔k-1个位移速度状态转移到异常状态的转移次数；q21为位移速度状态序列中从异常状态间隔k-1个位移速度状态转移到正常状态的转移次数；q22为位移速度状态序列中从异常状态间隔k-1个位移速度状态转移到下一个异常状态的转移次数；
[0022] 步骤4-1-2：根据状态转移频数矩阵计算状态转移概率矩阵P(k)：
[0023]
[0024] 其中，
[0025] 步骤4-1-3：计算初始状态概率向量
[0026]
[0027] 其中，
[0028] 步骤4-1-4：计算m阶位移速度次日预测状态的分布矩阵B：
[0029]
[0030] 其中，
[0031] 步骤4-2：使用加权马尔可夫链预测位移速度次日预测状态：
[0032] 步骤4-2-1：计算位移速度的自相关系数rk：
[0033]
[0034] 其中，为位移速度均值，xt+k为第t+k日的位移速度，t＝[1，n-k]；
[0035] 步骤4-2-2：计算不同步长权重wk，构成m阶权重向量wm＝(w1,w2,…,wk,…,wm)，其中，
[0036] 步骤4-2-3：计算m阶位移速度次日预测状态的概率分布向量C(n+1)：
[0037] C(n+1)＝(c1,c2,…,cr)＝wm*B
[0038] 步骤4-2-4：将向量C(n+1)中最大元素所在列对应的状态作为m阶加权马尔可夫链预测的次日位移速度xn+1所属状态，即位移速度次日预测状态；
[0039] 步骤4-3：重复步骤4-1到步骤4-2，计算m取不同值时各阶位移速度次日预测状态，如果各阶位移速度次日预测状态均为“异常”，位移速度次日预测状态为“异常”；如果各阶位移速度次日预测状态不全为“异常”，位移速度次日预测状态为“正常”；
[0040] 步骤5：进行滑坡预警判断：当满足位移速度前日状态为异常状态、位移速度当日状态为异常状态且位移速度次日预测状态为异常状态时，发出滑坡预警信号；
[0041] 步骤6：判断是否有新收集到的位移速度数据，有，则将第一日位移速度从位移速度集合中删除，将第n+1日收集到的位移速度加入位移速度集合中，对新的位移速度集合重复步骤3到步骤5；否则，执行步骤7；
[0042] 步骤7：对滑坡预警结果进行评价：
[0043] 步骤7-1：评价滑坡预警的及时性：滑坡预警及时性以滑坡发生前滑坡前应急准备时间内滑坡预警密集度衡量；滑坡预警密集度越接近1，滑坡预警的及时性越强；滑坡预警密集度为：在滑坡发生之前R天进行的R次滑坡预警判断中，发出滑坡预警信号的次数为R*，则滑坡之前R天内的滑坡预警密集度
[0044] 步骤7-2：评价滑坡预警的抗干扰性：错误预警率越低，滑坡预警的抗干扰性越好；错误预警率为：如果在滑坡应急准备时间之前共做出M次滑坡预警判断，其中有M*次发出了滑坡预警信号，错误预警率为
[0045] 步骤7-3：评价滑坡预警的可信度：滑坡预警的可信度由预测一致率决定，预测一致率越大，滑坡预警的可信度越好；预测一致率为：在滑坡发生前一共对位移速度状态做Q次预测，其中Q*次与实际位移速度状态相符合，则预测一致率为
[0046] 有益效果：一种基于均值-标准差的加权马尔可夫滑坡预警方法与现有技术相比，具有如下优势：
[0047] (1)从滑坡演变过程中获得的位移参数的随机特性出发，基于马尔可夫链理论，综合运用滑坡演化三阶段理论和均值-标准差分类法，超越了以往工程界对滑坡预警的认识；
[0048] (2)给出了具体预警过程，可用于从位移监测数据角度出发的滑坡预警。
附图说明
[0049] 图1为本发明一种实施方式的一种基于均值-标准差的加权马尔可夫滑坡预警方法流程图；
[0050] 图2为本发明一种实施方式的滑坡区及监测点布置平面图；
[0051] 图3为本发明一种实施方式的滑坡前各监测点预警情况示意图，(a)为滑坡前第8天及之前各监测点预警情况示意图，(b)滑坡前第7天各监测点预警情况示意图，(c)滑坡前第6天各监测点预警情况示意图，(d)滑坡前第5天各监测点预警情况示意图，(e)滑坡前第4天各监测点预警情况示意图，(f)滑坡前第3天各监测点预警情况示意图，(g)滑坡前第2天各监测点预警情况示意图，(h)滑坡前第1天各监测点预警情况示意图。
具体实施方式
[0052] 下面结合附图对本发明的一种实施方式作详细说明。
[0053] 本实施方式中，以露天矿为例，该露天矿于2013年4月17日发生了侏罗系地层滑坡。滑体高度近200m，走向长700多m，倾向宽度近400m，滑体平面投影面积约12万m2，滑面最3
大深度约50m，体积约300万m。
[0054] 如图1所示，基于均值-标准差的加权马尔可夫滑坡预警方法，包括如下步骤：
[0055] 步骤1：布置滑坡监测点，针对每一个滑坡监测点，执行步骤2到步骤7；如图2所示，滑坡区附近分别在▽438、▽462、▽512及▽548共布置了28个滑坡监测点，其中滑坡区外布置了7个滑坡监测点；
[0056] 步骤2：每日收集一次滑坡监测点位移速度，构成位移速度集合{x1，x2，…，xi，…，xn}，其中，xi为第i日收集到的位移速度，i＝{1,2，…，n}；本实施方式中，选择监测时间区段为2013年1月1日至2014年4月17日，选择样本容量为n＝20，将前20天获得的20个位移速度数据组成位移速度集合；鉴于该滑坡过程中的位移以水平方向为主，本实施方式只对水平位移数据进行研究；另一方面，由于数据量庞大，仅列出滑坡区内WY2300-512监测点数据。
[0057] 步骤3：利用均值-标准差法确定位移速度集合中各位移速度所属状态，包括正常状态或异常状态，n个位移速度的状态组成位移速度状态序列：
[0058] 步骤3-1：计算位移速度的均值 和标准差s：
[0059]
[0060] 本实施方式中，均值 标准差s＝0.6982；
[0061] 步骤3-2：确定状态划分的分界点D：
[0062]
[0063] 其中，α为可变参数；本实施方式中，α＝0.4，分界点D＝1.3008。
[0064] 步骤3-3：计算位移速度的最大值max和最小值min；
[0065] 步骤3-4：利用“D”“max”和“min”依次确定位移速度状态：若位移速度在区间[min，D)内，则位移速度状态为正常状态，用“1”表示；若位移速度在区间[D，max]内，则位移速度状态为异常状态，用“2”表示；n个位移速度状态组成位移速度状态序列；
[0066] 本实施方式中，位移速度的最小值min为-0.12，最大值max为2.39，若位移速度在区间[-0.12，1.30)内，则位移速度状态为异常状态，用“1”表示；若位移速度在区间[1.30，2.39]内，则位移速度状态为异常状态，用“2”表示。位移速度状态结果见表1：
[0067] 表1WY2300-512监测点位移速度状态表
[0068]
[0069] 步骤4：根据位移速度状态序列，使用加权马尔可夫链预测次日位移速度xn+1所属状态，即位移速度次日预测状态：
[0070] 步骤4-1：使用马尔可夫链计算m阶位移速度次日预测状态的分布矩阵；
[0071] 步骤4-1-1：根据位移速度状态序列，统计不同步长的状态转移频数矩阵q(k)：
[0072]
[0073] 其中，k为步长，k＝{1,2，…，m}，m为阶数，m≤5；q11为位移速度状态序列中从正常状态间隔k-1个位移速度状态转移到下一个正常状态的转移次数；q12为位移速度状态序列中从正常状态间隔k-1个位移速度状态转移到异常状态的转移次数；q21为位移速度状态序列中从异常状态间隔k-1个位移速度状态转移到正常状态的转移次数；q22为位移速度状态序列中从异常状态间隔k-1个位移速度状态转移到下一个异常状态的转移次数；
[0074] 步骤4-1-2：根据状态转移频数矩阵计算状态转移概率矩阵P(k)：
[0075]
[0076] 其中，
[0077] 本实施方式中，m＝5，步长k分别为1，2，3，4，5的状态转移频数矩阵和状态转移概率矩阵如下：
[0078]
[0079]
[0080]
[0081] 步骤4-1-3：计算初始状态概率向量
[0082]
[0083] 其中，
[0084] 本实施方式中，第20个位移速度状态为2，步长为1的初始状态概率向量同理，第19，18，17，16个位移速度状态为分别为1，1，1，2，步长2、3、4、5对应的初始状态概率向量分别为：
[0085] 步骤4-1-4：计算m阶位移速度次日预测状态的分布矩阵B：
[0086]
[0087] 其中，
[0088] 步骤4-2：使用加权马尔可夫链预测m阶位移速度次日预测状态：
[0089] 步骤4-2-1：计算位移速度的自相关系数rk：
[0090]
[0091] 其中，为位移速度均值，xt+k为第t+k日的位移速度，t＝[1，n-k]；本实施方式中，各步长自相关系数如表2：
[0092] 表2各步长自相关系数表
[0093]
[0094] 步骤4-2-2：计算不同步长权重wk，构成m阶权重向量wm＝(w1,w2,…,wk,…,wm)，其中， 本实施方式中，各步长的权重向量如表3：
[0095] 表3各步长权重向量表
[0096]
[0097] 步骤4-2-3：计算m阶位移速度次日预测状态的概率分布向量C(n+1)：
[0098] C(n+1)＝(c1,c2,…,cr)＝wm*B
[0099] 步骤4-2-4：将向量C(n+1)中最大元素所在列对应的状态作为m阶加权马尔可夫链预测的次日位移速度xn+1预测状态；
[0100] 步骤4-3：重复步骤4-1到步骤4-2，计算m取不同值时各阶位移速度次日预测状态，如果各阶位移速度次日预测状态均为“异常”，位移速度次日预测状态为“异常”；如果各阶位移速度次日预测状态不全为“异常”，位移速度次日预测状态为“正常”；
[0101] 本实施方式中，各阶位移速度次日预测状态的概率分布如表4：
[0102] 表4位移速度数据状态概率分布表
[0103]
[0104] 步骤5：进行滑坡预警判断：当满足位移速度前日状态为异常状态、位移速度当日状态为异常状态且位移速度次日预测状态为异常状态时，发出滑坡预警信号；
[0105] 未来位移速度的预测状态信息以及当前位移速度的真实状态信息，概括了马尔可夫链预测理论对边坡稳定性状况的评价，基于这些信息即可对是否给出预警做出判断。
[0106] 为确保预测的可信度，仅在各阶加权马尔可夫链的位移速度次日预测状态均为异常状态时，综合预测结果才确定为异常状态。由于监测数据不断更新，应及时纳入新数据，剔除历史数据，并不断以加权马尔可夫链预测未来位移速度状态。
[0107] 在给出每个次日位移速度预测状态的同时，也要做出预警与否的决定。当连续出现真实状态和预测状态均为异常状态时，应发出预警信号。这个连续出现异常状态的阶段对应于滑坡演化三阶段理论中位移加速度大于零的阶段，根据加速度a≥0这一判据的滑坡预警判断方法为：
[0108]
[0109] 步骤6：判断是否有新收集到的位移速度数据，有，则将第一日位移速度从位移速度集合中删除，将第n+1日收集到的实际位移速度加入位移速度集合中，对新的位移速度集合重复步骤3到步骤5；否则，执行步骤7；
[0110] 本实施方式中，监测点位移速度预测状态结果及预警情况如表5：
[0111] 表5WY2300-512监测点位移速度预测状态结果及预警情况表
[0112]
[0113]
[0114] 步骤7：对滑坡预警结果进行评价：
[0115] 步骤7-1：计算滑坡预警的及时性：滑坡预警及时性以滑坡发生前滑坡前应急准备时间内滑坡预警密集度衡量；本实施方式中，滑坡前应急准备时间为7天，故滑坡预警及时性以滑坡发生前7天内滑坡预警密集度，即7日滑坡预警密集度来衡量；对于发生滑坡的监测点，滑坡预警密集度应是一个动态变化逐渐接近1的变量；滑坡预警密集度越接近1，滑坡预警的及时性越强；而对于未发生滑坡的监测点，滑坡预警密集度则应是一个接近0的变量。滑坡预警密集度为：在滑坡发生之前R天进行的R次滑坡预警判断中，发出滑坡预警信号的次数为R*，则滑坡之前R天内的滑坡预警密集度
[0116] 步骤7-2：评价滑坡预警的抗干扰性：根据实际需要，错误预警率越低，滑坡预警的抗干扰性越好；错误预警率为：如果在滑坡应急准备时间之前共做出M次滑坡预警判断，其中有M*次发出了滑坡预警信号，错误预警率为
[0117] 步骤7-3：评价滑坡预警的可信度：滑坡预警的可信度由预测一致率决定，预测一致率越大，滑坡预警的可信度越好；预测一致率为：在滑坡发生前一共对位移速度状态做Q次预测，其中Q*次与实际位移速度状态相符合，则预测一致率为
[0118] 本实施方式中，由表5可知，对监测点的7日预警密集度为57％，错误预警率为0，预测一致率为58％。实际情况是2013/4/17发生滑坡，在83次预测中仅在滑坡前的第4、3、2、1日连续给出四次预警，满足预警及时性标准，而错误预警率为零则满足了预警抗干扰性标准，预测一致率为58％，具有一定的可信度。
[0119] 对滑坡区域内全部监测点的预警结果评价见表6：
[0120] 表6基于均值-标准差的加权马尔可夫滑坡预警方法预警评价表
[0121]
[0122] 滑坡区域内外各监测点临滑前15天的预警判断结果，见表7。表中“--”表示当日未给出预警信号，灰色框表示当日给出预警信号。
[0123] 表7监测点在滑坡前15天内的预警结果表
[0124]
[0125]
[0126] 监测区域内全部监测点在滑坡前的预警情况见图3(a)-图3(h)，从图中监测点的动态预警情况可以看出，本方法具有较好的时间预警和区域预警效果，能在一定程度上揭示滑坡的启动及演变过程。
