CN110197211A - 一种面向大坝安全监测数据的相似性数据聚类方法 - Google Patents

Abstract

本发明公开了一种面向大坝安全监测数据的相似性数据聚类方法，步骤为：利用EMD算法分离单测点序列趋势项与高频噪声，采用累积和控制图的拐点检测方法对时间序列变化点进行检测，分裂得到所有子序列集合；在计算子序列距离问题上采用DTW距离度量方法，动态计算两条子序列数据之间的距离最小值；利用层次聚类将挖掘出的子时间序列聚类，通过得到的树状聚类图，动态分析不同聚类数下的时间序列聚类分布情况。本发明合理分析监测数据相似性，可以挖掘出同一监测点在时间序列上的相关性，同时能量化安全监测数据之间相关性大小。经相似性分析处理后的监测数据，能准确反映大坝在时间维度上的变化趋势，结合变化趋势规律能有效降低后续监测数据挖掘难度。

Description

一种面向大坝安全监测数据的相似性数据聚类方法

技术领域

本发明属于大坝安全监测时空序列预测技术领域，特别涉及一种面向大坝安全监测数据的相似性数据聚类方法。

背景技术

大坝安全监测数据是大坝运行状态的整体表现，在受相同或相似外界因素影响下，其各区域之间数据必定存在相关性，这种相关性主要表现在时间序列趋势变化相似或空间趋势变化相似上。如大坝顺河向位移与各个坝段压力应变存在相似变化趋势，大坝中各区域温度与应压力之间的存在相关性等。合理分析监测数据相似性，通过相似性分析不仅可以挖掘出哪些监测变量之间具有相关性，同时能量化安全监测数据之间相关性大小。经相似性分析处理后的监测数据，能准确反映大坝在时间与空间维度上的变化趋势，结合变化趋势规律能有效降低后续监测数据挖掘难度。

针对大坝安全监测中存在的大量时间序列数据，对序列的相似性分析能挖掘出其潜在的规律对后续研究提供可靠依据。在序列长度研究方面，根据序列长度是否一致，将序列相似性度量分为锁步距离与动态距离，分析锁步距离计算简单复杂度不高但可扩展性低，动态距离计算复杂但扩展性强。基于不同时间序列相似性度量方法的特点，序列相似性度量被广泛应用在经济、医学、交通等领域用来解决序列间的相似性计算、数据补全等问题。在医学领域，对医学上的心电图、脑电图、肌电图等时间序列进行相似性分析，提取不同类型序列变化模板，在未知疾病类型的情况下对疾病或症状进行分类。上述介绍通过序列整体分析相似性，即在判断两条序列相似考虑从一条序列变化为另一条序列整体最小代价。

发明内容

发明目的：为了克服现有技术中传统结构力学模型对数据处理单一，且受限于技术人员的工程认知和专业知识水平，监测异常甄别的质量和效率偏低的问题，本发明提供一种面向大坝安全监测数据的相似性数据聚类方法，合理分析监测数据相似性，不仅可以挖掘出哪些监测变量之间具有相关性，同时能量化安全监测数据之间相关性大小；经相似性分析处理后的监测数据，能准确反映大坝在时间维度上的变化趋势，结合变化趋势规律能有效降低后续监测数据挖掘难度。

技术方案：为实现上述目的，本发明提供一种面向大坝安全监测数据的相似性数据聚类方法，包括如下步骤：

(1)在对一整条长时间序列进行处理前，采用自相关性检验方法分析该序列是否适合子序列挖掘；为保证数据平滑，采用EMD算法处理数据依据数据自身的时间尺度特征来进行信号分解。

(2)在时间序列当中任意一点x_t，能够利用其相邻的两个序列时间点x_t-1,x_t+1表示其瞬间序列变化趋势。累积和(Cumulative Sum,CUSUM)控制图的拐点检测方法多用于数据中突变点即序列变化点检测，通过持续分裂，最终得到所有子序列集合。

(3)对时间序列形态挖掘过程当中，距离衡量方法决定了挖掘聚类效果。采用DTW方法计算子序列间的距离，解决了大坝子序列间长度不相等和数据伸缩变形无法比较的问题，对连续数据漂移、伸缩、平移问题都具有很好的适应性。

所述利用自相关性检验方法和EMD方法对数据处理的具体步骤如下：

(1.1)时间序列自相关性检验：长度为n的时间序列数据S＝{x₁,x₂,...,x_t}，sub_i＝{x_i,...,x_i+k}为i到i+k连续时间范围内的观测序列，sub_j＝{x_j,.....,x_j+l}为j到j+l连续时间范围的观测序列。其中1≤i＜i+k≤n，1≤j＜j+l≤n，i≠j，k＞0，l＞0。若存在函数A，当A(sub_i)≈A(sub_j)则有sub_i与sub_i为一组自相似对，记作I＝<sub_i,sub_j>。；

(1.2)自相关性方法分析序列自相关性：自相关性方法(Auto CorrelationFunction ACF)在时间序列分析时，被用来分析同一条序列中的滞后自相关性规律：

E代表期望，k表示滞后的时间戳，σ代表数据标准差，S(t)表示序列中t时刻监测值，S(t+k)表示在t+k时刻的时间序列监测值，μ_t,μ_t+k表示同一条时间序列中在t，t+k时刻的期望值；

(1.3)EMD模态分解：EMD处理数据依据数据自身的时间尺度特征来进行信号分解，即局部平稳化，而无须预先设定基函数，将含有噪声的信号分解为一组变化单一的本征模函数(IMF Intrinsic Mode Function)，且IMF须满足以下两个性质：函数在整个时间范围内，局部极值点和过零点的数目必须相等，或最多相差一个；在任意时刻，局部最大值的包络(上包络线)和局部最小值的包络(下包络线)平均必须为零；

(1.4)计算均值包络线：

(1.5)计算剩余项：h₁(t)＝r₁(t)-m₁(t)，r₁(t)为原始序列S(t)，e₊(t)为r₁(t)的上包络线，e_-(t)为r₁(t)下包络线，m₁(t)表示r₁(t)序列的均值包络线；

(1.6)若h₁(t)满足上述IMF两个条件，则h₁(t)为原数据的第一个本征函数记作c₁(t)＝h₁(t)并令r₁(t)剩余项r₂(t)＝r₁(t)-c₁(t)继续寻找下一个IMF，若h₁(t)不满足条件，则h₂(t)＝h₁(t)-m₁(t)不断计算剩余项的上下包络线，减去上下包络线平均值，直到所有数据满足上述条件为止。最终得到原始时间序列n个本征函数c_i(t)与趋势项r_n(t)；

(1.7)根据IMF序列的特征频率变化区间，序列分割时能够依据频率范围选择指定误差范围内的IMF组合最终达到CUSUM分裂所要求的平滑数据。

所述内容利用累积和控制图的拐点检测方法分裂整条序列数据的具体步骤如下：

(2.1)求序列均值；

(2.2)设定累积和初值为0，计算各点累积和；

(2.3)求得该序列数据最大累积和：

split_point＝max(abs(B(i)))

B(i)函数表示每个时间点累积和结果，x_i为时间点i值，∑x_i/n为整个序列的平均值。式中split_point表示累积和序列中绝对值最大值点，通过得到绝对最大值下标，找出原始序列变化的累积最大点，对于大坝监测点顺河向位移数据中突变点也正表示序列模式转换点；

(2.4)迭代分割：基于已求得的拐点，将序列分割为两段，各段序列继续采用如上步骤找出累积和最大的拐点，进行迭代分割，直至分割序列长度均小于设定的最短序列，最终得到所有子序列集合。

所述内容采用标准化DTW算法计算子序列间距离；层次聚类将子序列进行聚类的具体步骤如下：

(3.1)计算两个子序列间的距离：对于两个子序列sub₁，sub₂，长度分别为l₁，l₂,两者DTW距离计算：

D(i,j)＝d(x_i,x_j)+min{D(i-1,j),D(i,j-1),D(i-1,j-1)}

D(i,j)表示在以l₁，l₂构造的l₁xl₂大小矩阵中下标为i，j位置的最大相似度。d(x_i,x_j)为x_i，x_j两点之间欧式距离。DTW计算过程中，在所有路径中找出最优路径使得两条时间序列数据之间距离最小；

(3.2)对距离进行标准化处理：DTW序列相似性衡量结果仅能够表示两条序列间转化所需要的最小距离。当有多条序列距离比较时会出现无法比较的情况，因此对聚类中多个子序列之间距离D_dtw大小进行标准化处理：

D_dtw(sub_i,sub_j)为动态时间弯曲距离，max(len(sub_i),len(sub_j))为两条子序列中较长序列长度。由DTW计算过程知D_dtw(sub_i,sub_j)为对齐两条序列所花费的最小整体代价。D_dtw(sub_i,sub_j)包含两条序列的长度因素，直接使用D_dtw作为聚类依据缺乏合理性。利用D_dtw(sub_i,sub_j)除以max(len(sub_i),len(sub_j))保证所有D_dtw(sub_i,sub_j)距离为均值，因此使用作为聚类判别标准；

(3.3)计算相似度距离矩阵：由于DTW距离满足对称性，因此关于SIM计算矩阵同样属于对称矩阵，计算过程中仅考虑计算矩阵上三角矩阵或下三角矩阵；

(3.4)利用距离矩阵进行聚类：聚类过程中利用SIM矩阵不断计算两个聚类之间最小平均距离，当两个聚类属于最小距离时则合并两聚类形成较大的聚类，最终将所有序列归并到同一聚类中。

有益效果：本发明与现有技术相比具有以下优点：

1.EMD算法将原本含有噪声及趋势的序列被分解为各个本征函数以及趋势序列，根据IMF序列的特征频率变化区间，序列分割时能够依据频率范围选择指定误差范围内的IMF组合，最终达到CUSUM分裂所要求的平滑数据，去除了原始数据中趋势项与高频噪声。

2.累积和(Cumulative Sum,CUSUM)控制图的拐点检测方法通过持续分裂，最终得到所有子序列集合，解决了长序列的分割问题，以便进行后续子序列的层次聚类。

3.采用DTW方法计算子序列间的距离，解决了大坝子序列间长度不相等和数据伸缩变形无法比较的问题，对连续数据漂移、伸缩、平移问题都具有很好的适应性。

附图说明

图1为具体实施例中子序列挖掘实例数据图；

图2为具体实施例中子序列挖掘总体思路图；

图3为具体实施例中累积和控制图拐点检测法分割序列示意图；

图4为具体实施例中监测序列自相关性检验图；

图5为具体实施例中EMD算法降噪分解示意图；

图6为具体实施例中监测序列层次聚类分析示意图。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

已知有1200高程29坝段中某测点顺河向位移监测数据。数据时间跨度为2012-01-01至2018-08-01共10158条监测值。此时间段内数据依据工程经验被先验分割标注为泄洪期、枯水期、蓄水期等10类序列段，如图1所示序列。

图2给出了本发明提供的面向大坝安全监测数据的相似性数据聚类方法的总体思路，它介绍本发明的主要工作流程，具体实施如下：

(1)在对一整条长时间序列进行处理前，采用自相关性检验方法分析该序列是否适合子序列挖掘；为保证数据平滑，采用EMD算法处理数据依据数据自身的时间尺度特征来进行信号分解，即局部平稳化，而无须预先设定基函数。

(1.1)时间序列自相关性检验：长度为n的时间序列数据S＝{x₁,x₂,…,x_t}，sub_i＝{x_i,…,x_i+k}为i到i+k连续时间范围内的观测序列，sub_j＝{x_j,…..,x_j+l}为j到j+l连续时间范围的观测序列。其中1≤i＜i+k≤n，1≤j＜j+l≤n，i≠j，k＞0，l＞0。若存在函数A，当A(sub_i)≈A(sub_j)则有sub_i与sub_i为一组自相似对，记作I＝<sub_i,sub_j>。；

(1.2)自相关性方法分析序列自相关性：自相关性方法(Auto CorrelationFunction ACF)在时间序列分析时，被用来分析同一条序列中的滞后自相关性规律：

E代表期望，k表示滞后的时间戳，σ代表数据标准差，S(t)表示序列中t时刻监测值，S(t+k)表示在t+k时刻的时间序列监测值，μ_t,μ_t+k表示同一条时间序列中在t，t+k时刻的期望值，相关性检验结果如图4所示，对于t时刻大坝顺河向位移值与滞后500、900、2400的顺河向位移值相关性最大，但随着时间的推移，顺河向位移自相关性在时间序列上逐渐降低；

(1.3)EMD模态分解：EMD处理数据依据数据自身的时间尺度特征来进行信号分解，即局部平稳化，而无须预先设定基函数，将性能不好的信号分解为一组性能较好的本征模函数(IMFIntrinsic Mode Function)，且IMF须满足以下两个性质：函数在整个时间范围内，局部极值点和过零点的数目必须相等，或最多相差一个；在任意时刻，局部最大值的包络(上包络线)和局部最小值的包络(下包络线)平均必须为零；

(1.4)计算均值包络线：

(1.5)计算剩余项：h₁(t)＝r₁(t)-m₁(t)，r₁(t)为原始序列S(t)，e₊(t)为r₁(t)的上包络线，e_-(t)为r₁(t)下包络线，m₁(t)表示r₁(t)序列的均值包络线；

(1.6)若h₁(t)满足上述IMF两个条件，则h₁(t)为原数据的第一个本征函数记作c₁(t)＝h₁(t)并令r₁(t)剩余项r₂(t)＝r₁(t)-c₁(t)继续寻找下一个IMF，若h₁(t)不满足条件，则h₂(t)＝h₁(t)-m₁(t)不断计算剩余项的上下包络线，减去上下包络线平均值，直到所有数据满足上述条件为止。最终得到原始时间序列n个本征函数c_i(t)与趋势项r_n(t)；

(1.7)EMD分解C4-A29-IP-01顺河向位移分解结果如图5所示，原本含有噪声及趋势的序列被分解为各个本征函数以及趋势序列。根据IMF序列的特征频率变化区间，序列分割时能够依据频率范围选择指定误差范围内的IMF组合最终达到CUSUM分裂所要求的平滑数据。

(2)在时间序列当中任意一点x_t，能够利用其相邻的两个序列时间点x_t-1,x_t+1表示其瞬间序列变化趋势。累积和(Cumulative Sum,CUSUM)控制图的拐点检测方法多用于数据中突变点即序列变化点检测，通过持续分裂，最终得到所有子序列集合。

(2.1)求序列均值；

(2.2)设定累积和初值为0，计算各点累积和；

(2.3)求得该序列数据最大累积和：

split_point＝max(abs(B(i)))

B(i)函数表示每个时间点累积和结果，x_i为时间点i值，∑x_i/n为整个序列的平均值。式(3-2)中split_point表示累积和序列中绝对值最大值点，通过得到绝对最大值下标，找出原始序列变化的累积最大点，对于大坝监测点顺河向位移数据中突变点也正表示序列模式转换点，如图3所示，C4-A29-IP-01监测点顺河向位移首次计算累积和得出的分裂点split_point过程；

(2.4)迭代分割：通过图3累积求和得出原始数据在时间点6480位置，累积和整体趋势发生变化，因此第一次时间序列分裂点在t＝6480。由t＝6480位置将原始顺河向位移数据分割成两条子序列sub₁，sub₂，设置子序列最小长度阈值为l_min，若sub₁或sub₂子序列长度大于阈值l_min且累积和存在极值点则继续分裂，最终得到所有子序列集合。

(3)对时间序列形态挖掘过程当中，距离衡量方法决定了挖掘聚类效果。采用DTW方法计算子序列间的距离，解决了大坝子序列间长度不相等和数据伸缩变形无法比较的问题，对连续数据漂移、伸缩、平移问题都具有很好的适应性。

(3.1)计算两个子序列间的距离：对于两个子序列sub₁，sub₂，长度分别为l₁，l₂,两者DTW距离计算：

D(i,j)＝d(x_i,x_j)+min{D(i-1,j),D(i,j-1),D(i-1,j-1)}

D(i,j)表示在以l₁，l₂构造的l₁xl₂大小矩阵中下标为i，j位置的最大相似度。d(x_i,x_j)为x_i，x_j两点之间欧式距离。DTW计算过程中，在所有路径中找出最优路径使得两条时间序列数据之间距离最小；

(3.2)对距离进行标准化处理：DTW序列相似性衡量结果仅能够表示两条序列间转化所需要的最小距离。当有多条序列距离比较时会出现无法比较的情况，因此对聚类中多个子序列之间距离D_dtw大小进行标准化处理：

D_dtw(sub_i,sub_j)为动态时间弯曲距离，max(len(sub_i),len(sub_j))为两条子序列中较长序列长度。由DTW计算过程知D_dtw(sub_i,sub_j)为对齐两条序列所花费的最小整体代价。D_dtw(sub_i,sub_j)包含两条序列的长度因素，直接使用D_dtw作为聚类依据缺乏合理性。利用D_dtw(sub_i,sub_j)除以max(len(sub_i),len(sub_j))保证所有D_dtw(sub_i,sub_j)距离为均值，因此使用作为聚类判别标准；

(3.3)计算相似度距离矩阵：由于DTW距离满足对称性，因此关于SIM计算矩阵同样属于对称矩阵，计算过程中仅考虑计算矩阵上三角矩阵或下三角矩阵；

(3.4)利用距离矩阵进行聚类：聚类过程中利用SIM矩阵不断计算两个聚类之间最小平均距离，当两个聚类属于最小距离时则合并两聚类形成较大的聚类，最终将所有序列归并到同一聚类中。C4-A29-IP-01序列进行子序列聚类部分结果如图6所示，根据子时间序列特点，选择自底向上方式层次聚类方式，聚类方法最初将所有子时间序列本身作为一个聚类簇，然后找出距离最近的两个簇，将相似的子序列合并为较大的簇，不断重复以上步骤直到所有的时间序列归纳到一个簇中。如图6中所示，对大坝时间序列相似分析中则是采用相似性度量方法不断的将相似时间子序列聚类。层次聚类应用在子时间序列聚类中不仅能直观得出聚类情况，同时能够通过得到的树状聚类图，动态分析不同聚类数下的时间序列聚类分布情况。

Claims (5)

1.一种面向大坝安全监测数据的相似性数据聚类方法，其特征在于，包括如下步骤：

(1)对待处理的完整的长时间大坝监测序列数据，采用自相关性检验方法分析该序列是否适合子序列挖掘；采用EMD算法处理数据依据数据自身的时间尺度特征来进行信号分解，保证数据平滑；

(2)采用累积和控制图的拐点检测方法，检测序列数据中的突变点，通过持续分裂，最终得到所有子序列集合；

(3)对时间序列形态挖掘过程当中，采用DTW方法计算子序列间的距离，比较大坝监测数不同子序列间距离，评价挖掘聚类效果。

2.根据权利要求1所述的一种面向大坝安全监测数据的相似性数据聚类方法，其特征在于，所述步骤(1)中采用自相关性检验方法和EMD方法对数据处理的具体步骤如下：

(1.1)时间序列自相关性检验：长度为n的时间序列数据S＝{x₁，x₂，…，x_t}，sub_i＝{x_i，…，x_i+k}为i到i+k连续时间范围内的观测序列，sub_j＝{x_j，…..，x_j+l}为j到j+l连续时间范围的观测序列；其中1≤i＜i+k≤n，1≤j＜j+l≤n，i≠j，k＞0，l＞0；若存在函数A，当A(sub_i)≈A(sub_j)则有sub_i与sub_i为一组自相似对，记作I＝<sub_i，sub_j>；

(1.2)自相关性方法分析序列自相关性：自相关性方法在时间序列分析时，被用来分析同一条序列中的滞后自相关性规律：

E代表期望，k表示滞后的时间戳，σ代表数据标准差，S(t)表示序列中t时刻监测值，S(t+k)表示在t+k时刻的时间序列监测值，μ_t，μ_t+k表示同一条时间序列中在t，t+k时刻的期望值；

(1.3)EMD模态分解：EMD处理数据依据数据自身的时间尺度特征来进行信号分解，即局部平稳化，而无须预先设定基函数，将含有噪声的信号分解为一组变化单一的本征模函数IMF，且IMF须满足以下两个性质：函数在整个时间范围内，局部极值点和过零点的数目必须相等，或最多相差一个；在任意时刻，局部最大值的包络即上包络线和局部最小值的包络即下包络线平均必须为零；

(1.4)计算均值包络线：

(1.5)计算剩余项：h₁(t)＝r₁(t)-m₁(t)，r₁(t)为原始序列S(t)，e₊(t)为r₁(t)的上包络线，e_{_}(t)为r₁(t)下包络线，m₁(t)表示r₁(t)序列的均值包络线；

(1.6)不断计算剩余项的上下包络线，减去上下包络线平均值，直到所有数据满足上述条件为止；最终得到原始时间序列n个本征函数c_i(t)与趋势项r_n(t)；

(1.7)根据IMF序列的特征频率变化区间，序列分割时能够依据频率范围选择指定误差范围内的IMF组合最终达到CUSUM分裂所要求的平滑数据。

3.根据权利要求1所述的一种面向大坝安全监测数据的相似性数据聚类方法，其特征在于，所述步骤(2)中采用累积和控制图的拐点检测方法分裂整条序列数据的具体步骤如下：

(2.1)求序列均值；

(2.2)设定累积和初值为0，计算各点累积和；

(2.3)求得该序列数据最大累积和，并记录该处的拐点x_i；

(2.4)迭代分割：基于已得的x_i，将序列分割为两段，各段序列继续采用如上步骤找出累积和最大的拐点，进行迭代分割，直至分割序列长度均小于设定的最短序列，最终得到所有子序列集合。

4.根据权利要求1所述的一种面向大坝安全监测数据的相似性数据聚类方法，其特征在于，所述步骤(3)中采用标准化DTW算法计算子序列间距离；层次聚类将子序列进行聚类的具体步骤如下：

(3.1)计算两个子序列间的距离：对于两个子序列sub₁，sub₂，长度分别为l₁，l₂，两者DTW距离计算：

D(i，j)＝d(x_i，x_j)+min{D(i-1，j)，D(i，j-1)，D(i-1，j-1)}

D(i，j)表示在以l₁，l₂构造的l₁xl₂大小矩阵中下标为i，j位置的最大相似度；d(x_i，x_j)为x_i，x_j两点之间欧式距离；DTW计算过程中，在所有路径中找出最优路径使得两条时间序列数据之间距离最小；

(3.2)对距离进行标准化处理：DTW序列相似性衡量结果仅能够表示两条序列间转化所需要的最小距离；当有多条序列距离比较时会出现无法比较的情况，因此对聚类中多个子序列之间距离D_dtw大小进行标准化处理：

D_dtw(sub_i，sub_j)为动态时间弯曲距离，max(len(sub_i)，len(sub_j))为两条子序列中较长序列长度；由DTW计算过程知D_dtw(sub_i，sub_j)为对齐两条序列所花费的最小整体代价；D_dtw(sub_i，sub_j)包含两条序列的长度因素，直接使用D_dtw作为聚类依据缺乏合理性；利用D_dtw(sub_i，sub_j)除以max(len(sub_i)，len(sub_j))保证所有D_dtw(sub_i，sub_j)距离为均值，因此使用作为聚类判别标准；

(3.3)计算相似度距离矩阵；

(3.4)利用距离矩阵进行聚类：聚类过程中利用SIM矩阵不断计算两个聚类之间最小平均距离，当两个聚类属于最小距离时则合并两聚类形成较大的聚类，最终将所有序列归并到同一聚类中。

5.根据权利要求3所述的一种面向大坝安全监测数据的相似性数据聚类方法，其特征在于，所述步骤(2.3)中求得该序列数据最大累积和，并记录该处的拐点x_i的具体步骤如下：

求得该序列数据最大累积和：

split_point＝max(abs(B(i)))

B(i)函数表示每个时间点累积和结果，x_i为时间点i值，∑x_i/n为整个序列的平均值；其中split_point表示累积和序列中绝对值最大值点，通过得到绝对最大值下标，找出原始序列变化的累积最大点，对于大坝监测点顺河向位移数据中突变点也正表示序列模式转换点。