基于位移观测序列聚类的边坡潜在滑移区早期辨识方法研究

白 宇，孙郡庆，刘明明，李卫鹏

(国能新疆托克逊能源有限责任公司，新疆维吾尔自治区吐鲁番市，830011)

0 引言

边坡安全稳定问题长期困扰着露天煤矿的高效开采，是露天煤矿安全生产不可回避的关键问题[1-2]。近年来，随着GPS、边坡雷达、柔性测斜仪等新型自动化监测设备的广泛应用，露天煤矿边坡变形监测技术得到长足发展，为主动防控边坡灾害提供越来越多的现场数据支持[3-6]。相比而言，露天煤矿边坡风险判识技术明显滞后，大量观测数据没有得到充分挖掘与利用[7-9]。目前，边坡位移观测数据主要用于临滑险情预警，难以在潜在风险区的早期识别方面发挥作用。滑坡预测缺乏多源监测数据融合分析手段，通常仅围绕某个监测点的观测数据进行分析，边坡状态判识存在主观性和盲目性。

从系统辨识的角度看，边坡立体监测网不同空间点位处观测到的监测数据集蕴含边坡系统在外部应力作用下，由无序稳态向有序非稳态转换过程中时空全域演化的信息[10-12]。为挖掘隐藏在多源监测数据背后的边坡系统演化规律，一些学者引入数据挖掘领域中的聚类分析概念，对边坡观测数据进行特征提取和知识挖掘。秦雨樵等[13]综合考虑边坡各点的位移以及点安全系数，提出基于K-means聚类算法的滑面搜索方法；雷从雨[14]将边坡节点位移聚类分析与变点分析相融合，提出改进的特征点位移突变判据；陈波等[15]采用划分聚类和投影聚类对边坡位移监测资料进行测点区域划分和特征提取，以达到边坡监测数据简约的目的。

作为按观测时间先后顺序排列的一组数据序列，边坡位移观测序列表现出十分明显的非线性、非平稳特征，不仅序列内部相邻观测值存在关联性，而且序列之间存在不等长、偏移性、时间轴扭曲等结构差异。然而，现有研究工作普遍忽略了位移观测序列结构的复杂性，或者仅对观测序列中特定时刻的观测值进行聚类[13-14]，或者虽然将整个观测序列看作聚类对象，但采用欧氏距离进行同步硬性度量[15]，聚类结果的准确性和可靠性无法得到保证。为此，笔者提出一种将动态时间弯曲(DTW)距离和K-means(K均值)聚类相结合的无监督数据挖掘方法，从多源监测数据集中聚类“观测域上相近、空间域上连续、变形趋势明显”的监测点子集，将其所在区域圈定为潜在滑移区。实验结果表明，上述聚类方法能够在滑动面尚未完全形成阶段，提前辨识出边坡潜在滑移区域。

1 边坡动力灾害演化规律

露天煤矿边坡属于复杂非线性系统范畴，其动力演化遵循一般非线性系统演化所具有的普遍规律，即从无序到有序的演化历程[10-11]。在边坡灾害的早期阶段，各子系统基本保持力学平衡，整个系统处于最无序的平衡态；随着边坡高度、暴露面积、暴露时间的增加，在大气降水、开采扰动等外部应力的作用下，局部区域岩土体单元力学平衡状态被打破，应力重新分布，表现出短程相关和局部有序；在边坡灾害的中后期阶段，随着微损伤从小尺度到大尺度的串级发展以及主控滑动面的逐渐形成，各子系统趋向协同作用，表现为长程相关和全局有序，直至整个边坡系统发生突变而形成滑坡。

当边坡子系统发展到有序阶段时，所处区域上的监测点位移观测序列必然表现出相似的演化行为，即有序区域具有“观测域上相近、空间域上连续、变形趋势明显”的特点。因此，通过边坡立体监测网各测点的时空域数据进行聚类分析，可以对露天煤矿边坡系统的有序程度进行定量化表征，为边坡灾害的早期辨识与晚期预警预报提供支撑。

2 位移观测序列相似性度量方法

边坡位移观测序列是按观测时间先后顺序排列的数据序列，可以采用时间序列相似性度量方法进行相似性评价。时间序列相似性度量方法主要分为欧氏距离和DTW距离2种。欧式距离和DTW距离数据匹配如图1所示。

图1 欧式距离和DTW距离数据匹配

图1中蓝色曲线与红色曲线代表2个时间序列，灰色曲线代表2个时间序列元素之间的相似匹配关系，曲线两端代表2个元素。欧式距离是一对一元素相似匹配(蓝色曲线上的一个元素匹配红色曲线上的一个元素)，DTW距离是一对多元素相似匹配(蓝色曲线上的一个元素有可能匹配红色曲线上的多个元素)。

2.1 欧氏距离

欧氏距离(ED)是一种对时间序列进行同步硬性度量(“一对一”匹配)的方法。设A=(a1，a2，…，an)和B=(b1，b2，…，bn)分别代表2个具有相同长度的位移观测序列，那么它们之间的欧氏距离定义为：

(1)

式中：d(A，B)——欧氏距离；

n——位移观测序列的长度，即观测时间的总次数；

ai——位移观测序列A中第i个观测值；

bi——位移观测序列B中第i个观测值。

d值越小说明序列A和B相似度越高，当d=0时，可认为A=B。

2.2 DTW距离

DTW距离是一种对时间序列进行异步相似形态度量(“一对多”匹配)的方法[16]。假定2个位移观测序列A=(a1，a2，…，an)和B=(b1，b2，…，bn)，长度分别为n和m，则序列A和B的DTW距离计算过程如下。

2.2.1 距离矩阵

构造大小为n×m的距离矩阵：

(2)

式中：d(ai，bj)——序列A中第i个观测值ai到序列B中第j个观测值bj的距离。

2.2.2 弯曲路径

定义弯曲路径：

(3)

式中：wr——序列A中第i个点与序列B中第j个点之间的匹配关系。

并且，弯曲路径需要满足3个约束条件：

(1)边界约束。弯曲路径起于距离矩阵D的左下角，止于距离矩阵D的右上角，即w1=(1，1)，wk=(n，m)。

(2)单调性约束。对于相邻元素wr=(i，j)和wr+1=(i′，j′)，需要满足i≤i′和j≤j′，即弯矩路径只能沿着时间轴单向移动。

(3)连续性约束。对于相邻元素wr=(i，j)和wr+1=(i′，j′)，需要满足i′≤i+1和j′≤j+1，即弯曲路径上任意2个相邻的元素在距离矩阵D中也相邻。

2.2.3 DTW距离

定义序列A和B之间的DTW距离：

DTW(A，B)=min(∑(i，j)∈Wd(ai，bj))

(4)

由式(4)可知，DTW(A，B)对应的弯曲路径为最优路径，在此路径上的累积距离最小。基于动态规划，求解式(2)等价于构造累积距离矩阵Γ，其中的元素服从如下递归关系：

(5)

式中：γ(i，j)——累积距离。

当完成累积距离矩阵Γ计算后，元素γ(n，m)便为序列A和B的DTW距离，即：

DTW(A，B)=γ(n，m)

(6)

进一步，以Γ的右上角为起点，即wk=(n，m)，以累积距离最小为选取原则，逆序寻找W中各元素，直至Γ的左下角，即w1=(1，1)，可以得到完整的最优弯曲路径。

DTW距离不仅能实现不等长时间序列的度量，而且对时间序列的偏移、形状缩放、振幅变化、时间点异常数据具有较强的鲁棒性。因此，笔者采用DTW距离对位移观测序列的相似性进行分析。

3 位移观测序列相似性聚类方法

聚类分析一直是数据挖掘和机器学习研究的热点之一，其目的是将数据集中的数据点按各自的相似度划分为不同的簇类，并确保簇内的数据点相似程度最大、簇间数据点相似程度最小，从而判识数据集合的结构，提取需要的重要信息[17]。为深度挖掘位移观测序列背后蕴含的边坡系统演化规律，采用基于DTW距离的K-means聚类算法，对位移观测序列进行相似性聚类分析。

3.1 算法原理

假定数据集X由n个位移观测序列组成，Χ=(Χ1，Χ2，…，Χn)，每个序列含m个观测值，Χi=(χi1，χi2，…，χim)。设定聚类数k，则聚类算法目标是找到一组簇类集合{S1，S2，…，Sk}，使得簇间观测序列尽可能不同，而簇内观测序列尽可能相似。其中，每个簇的聚类中心由cj=(cj1，cj2，…，cjm)表示，1≤j≤k。具体K-means算法框架如下。

(1)初始化阶段。

步骤1：在数据集X中，随机选取k个观测序列作为初始聚类中心。

(2)聚簇阶段。

步骤2：遍历数据集X中所有的观测序列，根据各观测序列到每个聚类中心的DTW距离，将观测序列分配到距离最近的聚类中心所代表的簇中，构建k个子集{S1，S2，…，Sk}。

步骤3：计算每个子集的位移观测均值，作为新的聚类中心。

重复步骤2～3，直到k个聚类中心不再变化，或达到指定的迭代次数。在此过程中，每个观测序列有且仅属于一个簇，并且与该簇聚类中心之间的距离最小。

3.2 聚类中心计算

每个簇的聚类中心由该簇内观测序列各个时间点观测值的均值构成，相应计算公式如下：

(7)

式中：cjt——簇Sj聚类中心cj第t个元素；

χlt——簇Sj内观测序列Χl第t个时间点的观测值；

|Sj|——簇Sj中观测序列的个数。

4 边坡潜在滑移区智能辨识

依据露天煤矿边坡系统所表现出的从无序到有序的动力演化特征，提出一种基于多源数据融合的边坡潜在滑移区辨识方法，通过对监测数据集的聚类分析，将“观测域上相近、空间域上连续、变形趋势明显”的监测点选定为有序测点，将有序测点所在区域圈定为潜在滑移区。具体辨识流程如下。

(8)

(2)对监测数据集F进行数据清洗，消除数据集中的缺失值、重复值、异常值。

(3)设定聚类数k，采用基于DTW距离的K-means聚类算法，对F观测域上的位移序列进行相似性度量及聚类分析，将F分解为k个簇，每个簇中测点位移观测序列行为相近。计算同簇测点位移观测序列的均值，获得每个簇的聚类中心。

(4)对各簇中心进行排序，选出[t0，tn]区间位移变化最大的簇中心，并对该簇测点的空间连续性进行评估。如果存在连续性，则将该簇测点选定为有序测点，将该簇测点所在区域圈定为潜在滑移区。滑坡滑动面智能辨识的流程如图2所示。

5 实验研究与结果分析

采用有限元数值建模方式仿真模拟典型边坡渐进失稳全过程，输出边坡不同演化阶段的位移场作为实验测试样本，对提出的边坡潜在滑移区辨识方法展开适用性研究。

5.1 数值建模与计算

采用文献[18-19]中的边坡算例：边坡高度H=20 m，坡角β=26.57°，杨氏模量E=100 MPa，泊松比ν=0.3，容重ϒ=20 kN/m3，黏聚力c=10 kPa，内摩擦角φ=20°。土体材料为理想弹塑性材料，强度准则采用平面应变条件下莫尔-库仑结合DP准则(DP4准则)。边坡尺寸及网格划分如图3所示，边坡左、右两侧边界为法向约束，底边为双向固定约束，共划分四边形单元5 686个、节点5 878个。采用有限元强度折减方法，按照c/ω，tanφ/ω方式不断折减土体材料强度参数，模拟边坡渐进破坏过程及其位移场演化信息，直至边坡发生整体失稳。其中，ω为折减系数，可看作广义时间，其初始值ω0=1.0，步长△ω=0.002，第m步的折减系数ωm=ωm-1+△ω。

图3 边坡尺寸及网格划分

折减系数分别为ω90=1.18、ω150=1.30、ω185=1.36和ω192=1.384时，边坡等效塑性应变云图和局部化带扩展路径如图4所示。由图4可以看出，ω90、ω150、ω185和ω192时刻分别对应于边坡渐进演化的早期、中期、晚期和失稳阶段；当边坡处于前3个阶段时，虽然滑带自坡脚向坡体内部逐渐扩展，但整个滑带尚未贯通，滑移区也没有完全形成，无法通过塑性云图和滑带扩展路径提前辨识潜在滑移区；当边坡进入失稳阶段、滑带完全贯通后，通过塑性云图和滑带扩展路径可以准确获取滑移区范围。

图4 不同折减系数工况下边坡内部等效塑性应变云图及局部化带扩展路径

5.2 监测样本构建

将边坡模型中的所有节点选做监测点，以ωi0为观测起始时间，以ωn为观测截止时间，从数值模型计算结果中提取观测数据，构建边坡监测数据集。

(9)

5.3 辨识算法测试

以观测起始时间为分类标准，将监测数据集细分为自边坡渐进演化的早期、中期、晚期启动观测的监测数据集，按照前文所述方法构建监测数据集。设定聚类数k=3，开展监测数据集观测域聚类。聚类结果通过测点空间分布图和簇中心曲线图可视化展示。其中，同簇测点及其簇中心曲线使用相同颜色。

5.3.1 边坡渐进演化早期阶段启动的观测数据集(简称“早期监测数据集”)

图5 同簇测点空间分布与簇中心曲线(早期监测数据集聚类结果)

5.3.2 边坡渐进演化中期阶段启动的观测数据集(简称“中期监测数据集”)

图6 同簇测点空间分布与簇中心曲线(中期监测数据集聚类结果)

5.3.3 边坡渐进演化晚期阶段启动的观测数据集(简称“晚期监测数据集”)

5.3.4 辨识结果分析

通过对早期、中期和晚期监测数据集的聚类挖掘，可以得到如下结论：

(2)在监测点数量足够多的前提下，边坡潜在滑移区辨识效果将主要取决于观测的起始时间和观测的累积时长。观测起始时间越早，滑移区可辨识的时间也越早，但所需的观测时长也越多。以ω90为观测起始时间，于ω155时刻可辨识出潜在滑移区，需要65个时间步；以ω150为观测起始时间，于ω175时刻可辨识出潜在滑移区，需要25个时间步；以ω185为观测起始时间，滑移区辨识仅需要5个时间步，但辨识时间也最晚。

(3)现场位移监测数据中的本底噪声对辨识效果有很大影响，只有当观测累积位移值明显大于本底噪声时，才有可能对潜在滑移区进行聚类辨识。因此，选用观测精度高的监测设备对于潜在滑移区的早期辨识十分重要。

6 结论

(1)边坡立体监测网不同空间点位处观测到的监测数据集蕴含边坡系统在外部应力作用下，由无序稳态向有序非稳态转换过程中时空全域演化信息。为深度挖掘隐藏在多源监测数据背后的边坡系统演化规律，提出一种将DTW距离和K-means聚类相结合的无监督数据挖掘方法，从多源监测数据集中聚类“观测域上相近、空间域上连续、变形趋势明显”的监测点子集，将其所在区域圈定为潜在滑移区。

(2)采用强度折减有限元数值模型仿真模拟典型边坡渐进失稳全过程，输出边坡不同演化阶段的位移场作为实验测试样本，对潜在滑移区辨识方法展开适用性研究。实验结果表明，该聚类方法能够捕捉边坡系统从无序到有序的渐进演化特征，可以在滑动面尚未完全形成阶段，提前辨识出边坡潜在滑移区域。

(3)需要指出的是，采用数值模型构建监测样本与边坡现场观测数据之间存在差异，并且将所有节点都作为监测点是理想情况。因此，加快构建露天煤矿边坡地表-深部立体监测网，开展基于现场有限测点数据的潜在滑移区早期识别，将成为未来研究的重点。