一种基于矩阵低秩恢复的时间序列突变点检测方法与流程

本发明涉及特别适用于特定功能的数字计算设备或数据处理设备或数据处理方法的技术领域，特别涉及一种有效抑制稀疏噪声污染导致的检测虚警率上升、有效提高检测准确率的基于矩阵低秩恢复的时间序列突变点检测方法。

背景技术

时间序列的突变点检测在数据挖掘和信息安全领域应用广泛，其旨在从给定的序列当中发现性质突变的采样点。

现有技术中，常用的时间序列突变点检测方法是基于子空间的方法，该类方法使用时间片段中的采样点构造迹矩阵，并利用迹矩阵张成空间的差异性来描述相邻时间片段之间的性质变化。奇异谱变换(singularspectrumtransform,sst)是其中具有代表性的序列突变点检测方法，该方法运用主成份分析(principalcomponentanalysis,pca)将相邻时间序列采样片段之间的差异性转化为子空间的距离进行描述。

但是，pca是基于高斯分布假设的主成份分析方法，仅当数据存在小幅值高斯噪声时，pca可准确找出数据分布的主成份方向。换而言之，数据存在大幅值噪声(或者毛刺噪声)干扰时，pca的性能会受到很大影响，因此，基于pca的sst检测方法在面对此类噪声时，检测性能也会降低，表现为准确度的下降和虚警率的上升。

近来，针对观测数据中的稀疏噪声污染问题，还有一种鲁棒的pca方法(robustpca,rpca)被提出，该方法将原观测数据矩阵通过凸优化建模方法分解为低秩信号矩阵和稀疏噪声矩阵，rpca在矩阵低秩恢复上的优异性能使其广泛应用于图像和生物信息学处理当中。

技术实现要素：

本发明解决的技术问题是，现有技术中，观测序列当中存在稀疏大幅值噪声，而导致pca的性能受到很大影响，基于pca的sst检测方法在面对这类噪声时，检测性能降低，表现为准确度的下降和虚警率的上升，为此，本发明提供了一种优化的基于矩阵低秩恢复的时间序列突变点检测方法。

本发明所采用的技术方案是，一种基于矩阵低秩恢复的时间序列突变点检测方法，所述方法包括以下步骤：

步骤1：对待检测的时间序列进行预处理，得到m个迹矩阵，初始化计数变量m＝1；

步骤2：利用鲁棒主成分分析方法对第m个和第m+d个迹矩阵进行低秩恢复，其中，d为时延，d＞0；分别取低秩恢复后的低秩矩阵作为恢复后的矩阵；

步骤3：将恢复后的第m个低秩矩阵进行子空间构造，计算恢复后的第m+d个低秩矩阵的特征向量；

步骤4：求解步骤3中构造的子空间和特征向量之间的欧氏距离e；

步骤5：m＝m+1，若m+d>m，则执行下一步，否则返回步骤2；

步骤6：对生成的m-d个欧氏距离e进行归一化，并与阈值w比较，确定突变发生的位置，为时间序列突变点；0＜w＜1。

优选地，所述步骤1中，预处理包括以下步骤：

步骤1.1：令待截取的子序列长度为l，对时间序列的任一时间节点向前取l-1个时间节点，以所述l-1个时间节点和所述任一时间节点为子序列；

步骤1.2：针对该子序列构造迹矩阵；

步骤1.3：若迹矩阵个数小于m，则返回步骤1.1，否则得到m个迹矩阵。

优选地，所述子序列的长度在整个时间序列上固定。

优选地，m的值为时间序列的长度减去子序列的长度加一。

优选地，所述步骤3包括以下步骤：

步骤3.1：对第m个迹矩阵对应的恢复后的低秩矩阵做奇异值分解，将l个最大的奇异值对应的左奇异向量作为基向量构造该迹矩阵对应的子空间；

步骤3.2：对第m+d个迹矩阵对应的恢复后的低秩矩阵做奇异值分解，将最大的奇异值对应的左奇异值向量取出，作为描述第m+d个迹矩阵变化模式的特征向量。

优选地，所述步骤3.1中，低秩矩阵奇异值分解得到w个奇异值，将w个奇异值从大到小排列为{σ1,σ2,...,σw}，l取值为满足的最小值l0，1≤i≤w，1≤j≤w；90％≤η≤99％。

优选地，所述步骤6中，归一化为对生成的m-d个欧氏距离e归一化至区间[0,1]。

本发明提供了一种优化的基于矩阵低秩恢复的时间序列突变点检测方法，通过对数据进行预处理，得到与时间序列采样片段对应的迹矩阵，基于鲁棒主成份分析从迹矩阵中恢复低秩矩阵，使用恢复的低秩矩阵构造子空间和特征向量，求解欧氏距离，得到变化得分，比较变化得分与设置阈值，最终确定突变发生的位置。

本发明区别于现有技术中rpca被广泛使用于解决多维数据的稀疏噪声问题，通过sst中的迹矩阵构造将rpca应用于一维时间序列当中，旨在通过对迹矩阵进行低秩恢复来滤除稀疏噪声。本发明主要针对观测序列当中存在的稀疏大幅值噪声问题，可以有效抑制稀疏噪声污染导致的检测虚警率上升的发生，同时可提升检测的准确率。

附图说明

图1为本发明的方法流程图；

图2为本发明的实施例的时间序列示意图，其中，t0为固定的时间节点，t0左侧有n个长度为w的子序列构造为迹矩阵h(t0)，以实线框表示，t0右侧为与迹矩阵h(t0)参与比较的、时延为d的同形状迹矩阵h(t0+d)，以虚线框表示。

具体实施方式

下面结合实施例对本发明做进一步的详细描述，但本发明的保护范围并不限于此。

本发明涉及一种基于矩阵低秩恢复的时间序列突变点检测方法，所述方法包括以下步骤。

步骤1：对待检测的时间序列进行预处理，得到m个迹矩阵，初始化计数变量m＝1。

所述步骤1中，预处理包括以下步骤：

步骤1.1：令待截取的子序列长度为l，对时间序列的任一时间节点向前取l-1个时间节点，以所述l-1个时间节点和所述任一时间节点为子序列；

步骤1.2：针对该子序列构造迹矩阵；

步骤1.3：若迹矩阵个数小于m，则返回步骤1.1，否则得到m个迹矩阵。

所述子序列的长度在整个时间序列上固定。

m的值为时间序列的长度减去子序列的长度加一。

本发明中，令给定时间序列为y(t)＝{y(1),y(2),…,y(t),…,y(n)}，子序列长度为l以及任意时间节点t0∈{l,l+1,...,n}，如图2所示。

本发明中，针对给定的y(t)，在t0处截取一个长度为w的时延向量v(t0-1)＝[y(t0-w),…,y(t0-1)]^t，其中，上标t表示矩阵的转置，基于n个这样的时延向量，构造迹矩阵h(t0)＝[v(t0-n),…,v(t0-2),v(t0-1)]，该迹矩阵h(t0)覆盖到时间序列的样本点为集合{y(t0-w-n+1),y(t0-w-n+2),...,y(t0-1)},该集合即作为时间节点t0对应的长度为l的子序列，故满足l＝w+n-1。

本发明中，步骤1.1中的子序列的长度l在整个时间序列上是保持固定的，m的值为时间序列的长度减去子序列的长度加一，在本实施例中子序列的长度固定为l，即m＝n-l+1。

本发明中，迹矩阵的详细描述已由参考文献knowledgediscoveryfromheterogeneousdynamicsystemsusingchange-pointcorrelations(tide，siaminternationalconferenceondatamining2005)公开，本领域技术人员可以获知相关技术。

步骤2：利用鲁棒主成分分析方法对第m个和第m+d个迹矩阵进行低秩恢复，其中，d为时延，d＞0；分别取低秩恢复后的低秩矩阵作为恢复后的矩阵。

本发明中，将第m个迹矩阵使用鲁棒主成分分析方法(rpca)分解为一个低秩矩阵和一个稀疏矩阵，将低秩矩阵作为恢复后的矩阵，同理对第m+d个迹矩阵做低秩恢复的操作。

本发明中，d为时延，人工设定，一般取值与迹矩阵的列数n相同。

本实施例，针对待处理的迹矩阵，以d∈r^w×n为例，rpca分解的目的是寻找一个低秩矩阵a和稀疏矩阵e，其中，低秩矩阵a即为待求的恢复后的迹矩阵，故针对迹矩阵的rpca求解问题可转化为如下优化问题：

其中，‖·‖*为核范数，‖·‖1为l1范数，λ为正则化参数。

本实施例中，上述优化问题通过不精确增广拉格朗日乘数(inexactaugmentedlagrangemultiplier，ialm)解决，步骤如下：

步骤2.1：初始化正则化参数λ、最大迭代次数kt、m维拉格朗日乘积向量y0＝d/j(d)、稀疏矩阵e0、标量μ0>0、迭代次数k＝0、阈值ε1和ε2，其中，j(d)＝max(‖y‖2,λ^-1‖d‖∞)，‖·‖∞是矩阵元素中的最大绝对值，‖·‖2为l2范数；

步骤2.2：求解当前循环中的低秩矩阵ak+1和稀疏矩阵ek+1。

具体地，首先进行奇异值分解如下：

其次进行求解如下：

其中，指的是收缩算子sε[x]，其中x具体为sε[x]定义如下：

参数ε>0。此收缩算子可扩展至矩阵，仅需使用该收缩算子处理矩阵中的每个元素。

步骤2.3：进行参数更新。

具体地，更新如下式：

yk+1＝yk+μk(d-ak+1-ek+1)(6)

其中，ρ>0为人为设定；k＝k+1。

步骤2.4：收敛判断。

具体为，若有不等式||d-ak+1-ek+1||f/||d||f<ε1满足，则判断为收敛，进行下一步；否则，跳转至步骤2.2。

步骤2.5：返回迭代结果(ak,ek)作为最终的矩阵分解结果。

步骤3：将恢复后的第m个低秩矩阵进行子空间构造，计算恢复后的第m+d个低秩矩阵的特征向量。

所述步骤3包括以下步骤：

步骤3.1：对第m个迹矩阵对应的恢复后的低秩矩阵做奇异值分解，将l个最大的奇异值对应的左奇异向量作为基向量构造该迹矩阵对应的子空间；

所述步骤3.1中，低秩矩阵奇异值分解得到w个奇异值，将w个奇异值从大到小排列为{σ1,σ2,...,σw}，l取值为满足的最小值l0，1≤i≤w，1≤j≤w；90％≤η≤99％。

步骤3.2：对第m+d个迹矩阵对应的恢复后的低秩矩阵做奇异值分解，将最大的奇异值对应的左奇异值向量取出，作为描述第m+d个迹矩阵变化模式的特征向量。

本发明中，对第m个迹矩阵对应的恢复矩阵做奇异值(svd)分解，将l个最大的奇异值对应的左奇异向量作为基向量构造该迹矩阵对应的子空间。在本实施例当中，将第m个迹矩阵对应的恢复矩阵表示为am，对矩阵am进行奇异值分解，将奇异值由大到小排列，取出前l个奇异值对应的左奇异向量u1,u2,…,ul，构造子空间hr≡span{u1,u2,…,ul}。

本发明中，对第m+d个迹矩阵对应的恢复矩阵做奇异值(svd)分解，将最大的奇异值对应的左奇异值向量取出，作为描述该迹矩阵变化模式的特征向量。具体地，在实施例当中，将第m+d个迹矩阵对应的恢复矩阵表示为am+g，对矩阵am+g进行奇异值分解，取出最大奇异值对应的左奇异向量β。

步骤4：求解步骤3中构造的子空间和特征向量之间的欧氏距离e。

本发明中，求解步骤3中的子空间和特征向量之间的欧氏距离e作为当前时间节点的“变化得分”(change-pointscore,cp)。

具体地，cp计算如下式：

其中，u＝[u1,u2,…,ul]。

本发明中，实际上，此欧氏距离e为向量至子空间的投影距离。

步骤5：m＝m+1，若m+d>m，则执行下一步，否则返回步骤2。

步骤6：对生成的m-d个欧氏距离e进行归一化，并与阈值w比较，确定突变发生的位置，为时间序列突变点；0＜w＜1。

所述步骤6中，归一化为对生成的m-d个欧氏距离e归一化至区间[0,1]。

本发明中，若归一化后的cp数值，即欧氏距离e，大于设置的阈值，则该cp对应的子序列将视为一个发生突变的片段。

本发明中，一个欧氏距离的计算涉及到两个矩阵，即第m个和第m+d个迹矩阵。因迹矩阵共有m个，故对应欧氏距离e的个数的上限ml满足ml+d＝m，则欧氏距离e的个数为m-d。

本发明通过对数据进行预处理，得到与时间序列采样片段对应的迹矩阵，基于鲁棒主成份分析从迹矩阵中恢复低秩矩阵，使用恢复的低秩矩阵构造子空间和特征向量，求解欧氏距离，得到变化得分，比较变化得分与设置阈值，最终确定突变发生的位置。本发明主要针对观测序列当中存在的稀疏大幅值噪声问题，可以有效抑制稀疏噪声污染导致的检测虚警率上升的发生，同时可提升检测的准确率。