一种基于动态时间规整与k-means聚类的边坡形变区域划分方法

本发明属于微波遥感和地质灾害预测预警的技术领域，尤其涉及一种基于动态时间规整与k-means聚类的边坡形变区域划分方法。

背景技术：

基于地基合成孔径雷达差分干涉测量技术，通过对同一位置、不同时刻获取的两幅雷达图像进行差分干涉处理，基于相位信息可以获得监测区域的位移信息。从大范围的目标区域中准确提取出各处形变区域，并结合各处形变区域的时间演变和空间分布信息等，有利于分析监测区域的整体稳定性，评估滑坡发生风险等级等。因此，边坡形变区域划分是滑坡稳定性监测中很重要的一步。

边坡形变区域划分方法，根据是否含有标签集，分为有监督学习算法和无监督学习算法两类。有监督学习算法主要通过搭建神经网络模型，进行监测区域内监测点的分类，但通常需要利用人为经验设置隐藏层个数、每层包含神经元个数等参数。无监督学习算法，如k-means聚类算法等，通常是以形变速率为单一指标进行各监测点的分类，由于形变速率无法评估监测点所处变形阶段，分类准确度较低。

理想位移-时间曲线包含初始变形阶段、等速变形阶段、加速变形阶段三个阶段，其滑坡演变阶段以累计加速度、累计加加速度为定量依据进行划分。位移时间序列能够描述一个监测点的变形演变情况及潜在变化，因此可引入各监测点所处形变状态作为边坡形变区域划分的评价标准。动态时间规整算法为一种模版匹配方法，通过遍历模版库内所有对象，选择出相似度足够高的结果，即认为其与该模版匹配。因此，可以通过动态时间规整算法计算不同监测点位移时间序列的相似度，来判断各监测点所处状态是否相近。

因此，针对常规边坡形变区域划分方法所存在的各项问题，有必要研究基于动态时间规整与k-means聚类的边坡形变区域划分方法。

技术实现要素：

为解决无监督边坡形变区域划分中计算量大、以位移速率为单一指标等问题，本发明提供一种基于动态时间规整与k-means聚类的边坡形变区域划分方法，能够实现边坡形变区域的有效划分。

一种基于动态时间规整与k-means聚类的边坡形变区域划分方法，包括以下步骤：

s1：在待监测边坡区域中，选取监测时间区间内累计位移量绝对值最大的监测点作为基准点；

s2：采用移动平均法分别从基准点与待监测边坡区域内其余监测点对应的位移时间序列中提取基准点与其余监测点的趋势项位移；

s3：分别对其余监测点的趋势项位移进行求导，得到其余监测点的变形速率序列，再将变形速率序列中变形速率零值的占比小于门限值的监测点作为备选监测点；

s4：采用动态时间规整算法分别获取各备选监测点的趋势项位移与基准点的趋势项位移之间的dtw相似度矩阵；

s5：将dtw相似度矩阵与监测时间区间内的累计位移量作为各备选监测点与基准点在k-means无监督聚类算法中对应的特征向量，并设定两个以上的类簇数目，再分别在不同类簇数目下，根据确定好特征向量的k-means无监督聚类算法对各备选监测点与基准点进行分类，得到各类簇数目对应的分类结果；

s6：按照设定的评价指标分别对各类簇数目对应的分类结果进行评估，将最大评估值对应的类簇数目作为边坡形变区域被划分的类别数，并将该类簇数目对应的分类结果作为边坡形变区域的划分结果。

进一步地，步骤s2中所述的基准点与其余监测点的趋势项位移的提取方法为：

s21：将基准点对应的位移时间序列表示为其余监测点对应的位移时间序列表示为其中，n为监测时间区间内的监测周期数量，m为监测点的序号；

s22：根据各监测点对应的位移时间序列来获取各监测点在第n个监测周期的移动平均结果得到各监测点的趋势项位移其中，移动平均结果的计算公式如下：

其中，t为移动平均法中设定的滑动周期，表示第m个监测点在第n个监测周期的累计位移量，且n＝t，t+1，…，n；

s23：根据基准点对应的位移时间序列来获取基准点在第n个监测周期的移动平均结果得到基准点的趋势项位移其中，移动平均结果的计算公式如下：

其中，表示基准点在第n个监测周期的累计位移量。

进一步地，采用max-min归一化方法对各备选监测点的趋势项位移与基准点的趋势项位移进行归一化后，再采用动态时间规整算法分别获取各备选监测点的趋势项位移与基准点的趋势项位移之间的dtw相似度矩阵。

进一步地，采用动态时间规整算法分别获取各备选监测点的趋势项位移与基准点的趋势项位移之间的dtw相似度矩阵具体为：

s41：将基准点的趋势项位移记为备选监测点的趋势项位移记为其中，为基准点在第t～n个监测周期的移动平均结果，为备选各监测点在第t～n个监测周期的移动平均结果，且n为监测时间区间内的监测周期数量，m为监测点的序号，t为移动平均法中设定的滑动周期；

s42：分别将各备选监测点作为当前监测点m执行以下步骤，得到各备选监测点的趋势项位移与基准点的趋势项位移之间的dtw相似度矩阵：

采用当前监测点m的趋势项位移与基准点的趋势项位移建立矩阵网格dm，其中，矩阵网格dm中的矩阵元素dm(i，j)表示趋势项位移第i个向量元素与趋势项位移第j个向量元素的欧氏距离其中，i＝1，2，…，n-t+1，j＝1，2，…，n-t+1；

根据当前监测点m对应的矩阵网格dm获取当前监测点m对应的dtw相似度矩阵rm，其中，dtw相似度矩阵rm第一列中各元素rm(i，1)的计算公式为：

(i＝1，…，n-t+1)

dtw相似度矩阵rm第一行中各元素rm(1，j)的计算公式为：

(j＝1，…，n-t+1)

dtw相似度矩阵rm其余行其余列各元素rm(i，j)的计算公式为：

(i＝2，…，n-t+1；j＝2，…，n-t+1)

其中，dtw相似度矩阵rm中的各元素rm(i，j)表示当前监测点m的趋势项位移的前i项构成的序列与基准点的趋势项位移前j项构成的序列的dtw相似度。

进一步地，所述根据确定好特征向量的k-means无监督聚类算法对各备选监测点与基准点进行分类具体为：

s51：构建目标函数如下：

其中，i⁽¹⁾～i^(r)为样本点对应的特征向量f(1)～f(r)所指定的簇的索引值，μv为对应簇的聚类中心，r＝1，2，3，…，r，r为基准点和备选监测点的总数，v＝1，2，3，…，v，v为聚类簇的个数，且所述样本点包括基准点和备选监测点；

s52：优化解算所述目标函数，以实现对各备选监测点与基准点的分类。

进一步地，步骤s6中所述的评价指标包括误差平方和点积pnsse、db指数以及ch分数，各类簇数目对应的分类结果的评估值的获取方法为：

s61：基于熵值法确定误差平方和点积pnsse、db指数以及ch分数的权值；

s62：分别获取各类簇数目对应的分类结果的误差平方和点积pnsse、db指数以及ch分数；

s63：分别将每一个类簇数目对应的分类结果的误差平方和点积pnsse、db指数以及ch分数按照对应的权值进行加权求和，所得和值为各类簇数目对应的分类结果的评估值。

进一步地，各类簇数目对应的分类结果的误差平方和点积pnsse的计算方法为：

其中，v为类簇数目，sse为类簇数目v对应的分类结果中，各样本点与其对应的聚类中心的欧式距离的平方和，其中，样本点包括备选监测点与基准点。

有益效果：

本发明提供一种基于动态时间规整与k-means聚类的边坡形变区域划分方法，首先提取基准点和其他监测点的形变数据的趋势项位移，然后采用动态时间调整算法计算监测点的趋势项位移与基准点的趋势项位移之间的相似度，最后基于k-means聚类实现形变区域划分，并通过计算多指标综合评价不同类簇数目下的分类结果，选取评价结果最优的类簇数目作为监测区域内各监测点的分类类别数量，从而获得监测区域内不同监测点位移时间序列分类信息，能够实现边坡形变区域划分结果的自动获取；由此可见，本发明有效弥补了以位移速率为单一指标的无监督学习算法在形变区域划分中的不足，为边坡形变区域划分提供了一种方便、准确的分类方法。

附图说明

图1为本发明实施例中边坡实例场景照片。

图2为本发明方法流程图。

图3为本发明实施例单一监测周期的累计位移量结果图。

图4为本发明实施例基准点位移时间序列变化结果图。

图5为本发明实施例基准点累积位移量变化及移动平均结果图。

图6为本发明实施例滤除点分布及监测区间内累计位移量结果图。

图7为本发明实施例监测区域dtw相似度结果图。

图8为本发明实施例多指标综合评价结果曲线变化结果图。

图9为本发明实施例边坡形变区域划分结果图。

具体实施方式

为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述。

本发明提出一种基于动态时间规整与k-means聚类的边坡形变区域划分方法，在滑坡灾害预测预警方面有重要用途。通过比较监测区间内累计位移量变化情况，选取基准点位移时间序列，经过归一化及移动平均预处理，提取各监测点累计位移量变化趋势；之后设置变形速率零值所占比门限值，选取需要进行后续分类的监测点集合后，采用动态时间规整算法计算各点位移时间序列与基准位移时间序列的相似度；最终选取累计位移量和对应相似度为特征，利用k-means聚类算法进行各监测点位移时间序列的分类，同时计算不同个数分类簇多指标综合评价结果得分，从而获得最优边坡形变区域划分结果。

以青海威斯特铜矿为例，说明本发明的具体实施过程。青海威斯特露天铜矿位于青海省果洛藏族自治州，监测场景照片如图1所示。实验采用地基合成孔径雷达(gb-sar)进行连续多天边坡形变监测，基于差分干涉技术实现了高精度形变监测。本发明给出青海威斯特铜矿的边坡形变区域划分方法，如图2所示，实施步骤如下：

s1：在待监测边坡区域中，选取监测时间区间内累计位移量绝对值最大的监测点作为基准点。

需要说明的是，步骤s1是为了获取基准点位移时间序列，具体的，本发明基于地基合成孔径雷达差分干涉测量技术，获取监测边坡区域的长时间形变信息。选取监测时间区间内累计位移量绝对值最大的监测点，以其位移时间序列作为基准点位移时间序列，设整个监测区间共包含等时间间隔的n个监测周期，基准时间序列表示为图3所示为累计位移量数据。基准点累积位移量变化如图4所示。

s2：采用移动平均法分别从基准点与待监测边坡区域内其余监测点对应的位移时间序列中提取基准点与其余监测点的趋势项位移。

需要说明的是，边坡的运动状态具有不确定性，因此监测区域内各点累计位移量通常是一个随时间变化的非稳定时间序列。在忽略地震、人类活动等随机扰动的情况下，通常将单点累计位移量分解为趋势项位移和周期项位移，其中趋势项位移可以表示该点累计位移量的变化趋势。

本发明采用移动平均法提取监测区域内各点及基准点位移时间序列的趋势项位移，包括以下步骤：

s21：将基准点对应的位移时间序列表示为其余监测点对应的位移时间序列表示为其中，n为监测时间区间内的监测周期数量，m为监测点的序号；

s22：根据各监测点对应的位移时间序列来获取各监测点在第n个监测周期的移动平均结果得到各监测点的趋势项位移其中，移动平均结果均计算公式如下：

其中，t为移动平均法中设定的滑动周期，表示第m个监测点在第n个监测周期的累计位移量，且n＝t，t+1，…，n；

s23：根据基准点对应的位移时间序列来获取基准点在第n个监测周期的移动平均结果得到基准点的趋势项位移其中，移动平均结果的计算公式如下：

其中，表示基准点在第n个监测周期的累计位移量，且基准点累积位移量变化及移动平均结果如图5所示。

s3：分别对其余监测点的趋势项位移进行求导，得到其余监测点的变形速率序列，再将变形速率序列中变形速率零值的占比小于门限值的监测点作为备选监测点。

需要说明的是，步骤s3实际为对待监测边坡区域中的监测点作剪枝预处理操作，其目的是在区域划分前，先初步检测出场景中的形变区域，滤除场景中在监测区间内基本保持稳定的监测点，从而降低模型运算量。设置变形速率零值所占比门限值α，筛选监测区间内各监测点变形速率零值点占比小于α的监测点集合用于形变区域划分，从而确定滑坡失稳关键区域，滤除点分布及监测区间内累计位移量如图6所示。

s4：采用动态时间规整算法分别获取各备选监测点的趋势项位移与基准点的趋势项位移之间的dtw相似度矩阵。

需要说明的是，由于监测区间内累计位移量大小与位移时间序列的变化趋势不直接相关，为了使得各点位移时间序列变化趋势具有可比性，因此可以采用max-min归一化方法对各备选监测点的趋势项位移与基准点的趋势项位移进行归一化，使得各点累计位移量处于同一数量级后，再采用动态时间规整算法(dtw)计算监测区域内位移时间序列的相似度，具体为：

s41：将基准点的趋势项位移记为备选监测点的趋势项位移记为其中，为基准点在第t～n个监测周期的移动平均结果，为备选各监测点在第t～n个监测周期的移动平均结果，且n为监测时间区间内的监测周期数量，m为监测点的序号，t为移动平均法中设定的滑动周期；

s42：分别将各备选监测点作为当前监测点m执行以下步骤，得到各备选监测点的趋势项位移与基准点的趋势项位移之间的dtw相似度矩阵：

采用当前监测点m的趋势项位移与基准点的趋势项位移建立矩阵网格dm，其中，矩阵网格dm中的矩阵元素dm(i，j)表示趋势项位移第i个向量元素与趋势项位移第j个向量元素的欧氏距离其中，i＝1，2，…，n-t+1，j＝1，2，…，n-t+1；

由于通过动态规划寻找从矩阵元素dm(1，1)至dm(n-t+1，n-t+1)距离累计值最小的路径即为最优路径，且时间序列相似度即为最优路径所得距离累计值，则可以根据当前监测点m对应的矩阵网格dm获取当前监测点m对应的dtw相似度矩阵rm，其中，dtw相似度矩阵rm第一列中各元素rm(i，1)的计算公式为：

(i＝1，…，n-t+1)

dtw相似度矩阵rm第一行中各元素rm(1，j)的计算公式为：

(j＝1，…，n-t+1)

dtw相似度矩阵rm其余行其余列各元素rm(i，j)的计算公式为：

(i＝2，…，n-t+1；j＝2，…，n-t+1)

其中，dtw相似度矩阵rm中的各元素rm(i，j)表示当前监测点m的趋势项位移的前i项构成的序列与基准点的趋势项位移前j项构成的序列的dtw相似度，也即从起始点dm(1，1)至dm(i，j)最优路径的距离累计值，因此第m个监测点的趋势项位移与基准点趋势项位移的dtw相似度即为矩阵元素rm(n-t+1，n-t+1)。

也就是说，动态时间规整算法即通过计算最优路径的距离累计值，比较监测点位移时间序列与基准点位移时间序列的相关性，监测区域内所有监测点的dtw相似度结果图如图7所示。

s5：将dtw相似度矩阵与监测时间区间内的累计位移量作为各备选监测点与基准点在k-means无监督聚类算法中对应的特征向量，并设定两个以上的类簇数目，再分别在不同类簇数目下，根据确定好特征向量的k-means无监督聚类算法对各备选监测点与基准点进行分类，得到各类簇数目对应的分类结果。

下面详细说明如何根据确定好特征向量的k-means无监督聚类算法对各备选监测点与基准点进行分类，实现边坡形变区域的划分。

选取监测区间内r个监测点(包括备选监测点和基准点)的第n个监测周期的累计位移量(r＝1，2，…，r)和dtw相似度矩阵rr(r＝1，2，…，r)为特征，构成单一监测点的特征向量f(r)(r＝1，2，…，r)，从而构建用于边坡形变区域划分的样本集

k-means聚类为本发明用于边坡形变区域划分的无监督算法，单一样本点通过比较对应特征向量与各聚类中心的欧氏距离进行簇分配，从而基于贪心算法实现监测区域内所有监测点的分类，则优化目标函数j可以表示为：

其中，i⁽¹⁾～i^(r)为样本点对应的特征向量f(1)～f(r)所指定的簇的索引值，μv为对应簇的聚类中心，r＝1，2，3，…，r，r为基准点和备选监测点的总数，v＝1，2，3，…，v，v为聚类簇的个数，且所述样本点包括基准点和备选监测点；

在第r(r＝1，2，3，…，r)次迭代中，前r-1次簇分配所得聚类中心为{μ1，μ2，…，μv}，则第r个样本点f(r)的分类过程可分为以下两个步骤：

i^(m)＝argmini∈{1，2，...，v}||fm-μi||2

其中，r＝1，2，3，…，r，f(r)为第r个监测点的特征向量，μi(i＝1，2，…，v)为前r-1次迭代所得聚类中心，同时将样本点f(r)划入相应的簇即为第i^(r)个簇更新后所包含的样本点特征向量集合；

其中，i＝1，2，3，…，v，μ′i为更新后所得聚类中心，ci为第i个聚类簇所包含样本点的特征向量集合，即通过移动聚类中心使得当前目标函数取得最小值。

s6：按照设定的评价指标分别对各类簇数目对应的分类结果进行评估，将最大评估值对应的类簇数目作为边坡形变区域被划分的类别数，并将该类簇数目对应的分类结果作为边坡形变区域的划分结果。

需要说明的是，对于形变区域内的单一监测点，k-means无监督聚类算法分别通过改变索引值和移动聚类中心进行分类，从而优化目标函数；进一步地，为了获得簇聚合程度较高的分类结果，本发明利用误差平方和点积pnsse、db指数以及ch分数评价无监督边坡形变区域划分结果，其中，本发明获得的边坡形变区域划分结果如图9所示。

下面详细介绍误差平方和点积pnsse、db指数以及ch分数的计算方法。

假设共选择k个不同的分类簇个数vi(i＝1，2，…，k)，则一共需要进行k次k-means预分类。设当前分类簇个数为vi，为了考虑分类簇个数vi及误差平方与sse结果值曲线变化的相关性对于分类结果的综合影响，构建评价指标分类误差平方和点积pnsse，可以表示为：

其中，sse为分类后各样本点与其对应聚类中心的欧式距离的平方和。根据k次计算结果构建pnsse结果值向量r1＝[pnsse1，pnsse2，…，pnssek]，pnssei为第i个分类簇个数对应的分类误差平方和点积。

为了衡量样本点分类结果内部的距离，引入db指数为评价指标。设簇cq内样本点间的平均距离为avg(cq)，可以表示为：

其中，|cq|为簇cq包含的样本点个数，基于上式可计算db指数所得结果值，为：

则db结果值向量可以表示为r2＝[db1，db2，…，dbk]。

为了评估各类之间的离散程度，引入ch分数为评价指标，可以表示为：

其中，m为监测点样本总数，trb(vi)表示类间离差矩阵的迹，trw(vi)表示类内离差矩阵的迹，则ch结果值向量可以表示为r3＝[ch1，ch2，…，chk]，因此可建立评价指标结果值矩阵

为了获取合理的多指标综合评价结果得分，本发明利用熵值法计算各评价指标权重分配，首先，计算各项评价指标在各次分类结果中所占比重，则当前更新结果值矩阵r′的矩阵元素可以表示为：

其中，i＝1，2，…，k，j＝1，2，3；然后，计算各项指标的熵值，可以表示为：

即可求得各项指标的信息熵冗余度dj，可以表示为：

dj＝1-ej

最终，可求得各项指标权值，即为：

则各次分类结果所得多指标综合评价结果得分si(i＝1，2，…，k)可以表示为：

因此，当前数据集设置下多指标综合评价结果得分变化曲线如图8所示，获得最优边坡形变区域划分结果如图9所示。

由此可见，本发明通过比较监测区间各监测点累计位移量变化情况，选取基准点位移时间序列；之后采用归一化及移动平均平滑预处理，提取监测区域内各监测点变化趋势，通过设置门限值进行剪枝预处理，筛选进行后续分类的监测点坐标位置集合，从而利用动态时间规整算法计算集合内位移时间序列相似度；最后以监测区间内累计位移量和dtw相似度矩阵为输入特征，建立以k-means无监督聚类算法的监测点位移时间序列分类模型，从而获得边坡形变区域划分结果；本发明能仅利用监测区域内各监测点位移时间序列得到该点趋势项位移变化情况，并通过结合剪枝预处理、动态时间规整算法、k-means无监督聚类算法得到更有效的监测点分类模型，从而为边坡形变区域划分提供了一种简洁、高效的分类方法。

当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当然可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。