一种基于全局-局部散度的多元时间序列无监督降维方法

1.本发明属于计算机技术领域，具体涉及一种基于全局-局部散度的多元时间序列无监督降维方法。


背景技术：

2.时间序列通常是指通过各种传感器采集得到的一系列随时间变化的数据，其广泛存在于环境学、医学、金融学等领域，根据变量的数目时间序列可以分为一元时间序列(univariate time series，uts)和多元时间序列(multiple timeseries，mts)。mts可以看作是同一个系统中不同因素产生的相应uts的组合，其相比于uts不仅包含时间维度的高维性，还具有特征维度上的高维性，且特征之间存在相关性。因此，针对mts的数据挖掘相比于uts更为复杂。
3.目前mts的数据挖掘应用于工业故障监测中，由于故障监测中的数据大多是通过传感器获取，因此存在大量的噪声信息，且传感器采集到的监控数据存在特征维度和时间维度，在进行故障监测时由于高纬度容易导致数据处理困难，效率低下。因此，如何对故障监测中的多元时间序列数据进行有效降维是解决这一问题的关键所在。
4.时间序列数据挖掘通常包括聚类、分类、预测、异常检测、相关性分析等，这些挖掘任务通常与数据的规模和复杂性有关，mts由于同时具有两个维度的高维特性，因此在进行上述挖掘任务时通常要进行降维或特征表示来降低数据的复杂性以及减轻冗余信息造成的干扰，目前，现有技术主要分为针对特征维度的降维，针对时间维度的降维和同时针对两者的降维。而针对特征维度的降维无法解决不等长mts的时间长度不等的技术问题，在后续数据挖掘的相似性度量中存在困难；针对时间维度的降维则忽略了特征维度中变量相关性特征，降维后还可能存在信息冗余等问题，而针对两者的降维通常需要使用双向降维技术，计算成本高。


技术实现要素：

5.本发明的目的在于克服现有技术的不足，提供一种基于全局-局部散度的多元时间序列无监督降维方法。
6.为实现上述目的，本发明采用的技术方案是：
7.一种基于全局-局部散度的多元时间序列无监督降维方法，其特征在于，包括以下步骤：
8.s1：通过传感器获取故障监测数据，将获取的故障监测数据多元时间序列组成多元时间序列原始数据集d，分别计算数据集中每个多元时间序列的协方差矩阵，提取协方差矩阵的上三角元素，将其组合为特征序列；并将所有特征序列构成特征序列集，且每一个特征序列的长度相同；
9.s2：基于得到的特征序列集，使用k近邻数和欧氏距离进行度量建立特征序列集中每个样本邻域集合；
10.s3：根据步骤s2得到的邻域集合，计算特征集中每个邻域集合的邻域中心序列；
11.s4：根据步骤s2得到的邻域集合，计算待投影后的样本点每个邻域集合的方差，再对这些方差累加求和，计算局部散度；
12.s5：根据步骤s3中得到的邻域中心点，计算邻域全局方差，得到全局散度；
13.s6：根据步骤s4和s5中得到的局部散度和全局散度，求解投影矩阵；
14.s7：根据步骤s6中得到的投影矩阵，对步骤s1中得到的特征序列集进行投影，从而得到降维特征序列：
15.yi＝w
tfi
ꢀꢀ
(6)；
16.s8：根据降维特征序列，得到故障监测数据的降维特征集d
′
＝{yi|i＝1,2,
…
,n}，其中，
17.s9：对降维后的故障监测数据进行处理，得到故障监测结果。
18.进一步地，步骤1的具体操作步骤包括：
19.s11：将获取的故障监测数据多元时间序列，组成多元时间序列原始数据集d＝{xi|i＝1,2,
…
,n}，其中n为样本个数，表示原始数据集中第i个mts样本，xi(i＝1,2,
…
,m)表示第i个变量的一系列观测值，m为变量数，ti是第i个多元时间序列的时间长度；并对每个多元时间序列进行零均值化处理：xi＝x
i-e(xi)；
20.s12：对每个零均值化处理后的多元时间序列进行协方差矩阵计算，且第i个mts的协方差矩阵为：
[0021][0022]
s13：将协方差矩阵的上三角阵元素进行提取，并按顺序组成行向量：
[0023]fi
＝(a
11
,a
12
,
…
,a
1m
,a
21
,
…
,a
2m
,
…
,a
mm
)
ꢀꢀ
(8)
[0024]
将fi作为第i个mts的特征序列，得到多元时间序列的特征集fea＝{fi|i＝1,2,
…
,n}。
[0025]
进一步地，步骤s3中计算邻域中心序列mi的计算公式为：
[0026][0027]
其中，k为近邻数，mi为fi与其k个近邻点的邻域中心序列。
[0028]
进一步地，s4的具体步骤包括：
[0029]
s41：基于所述邻域集合，计算投影后样本点每个邻域集合的方差，并对这些方差累加求和，得到局部散度计算公式：
[0030][0031]
其中，pi＝w
t
mi为邻域中心点的低维投影；yj是指样本点的低维投影，w为投影矩阵，下标l为英文local的缩写；
[0032]
s42：对式(2)进行变换可得：
[0033][0034]
其中，其中，为第i个特征序列fi的局部散度矩阵：
[0035][0036]
进一步地，s5的具体步骤包括：
[0037]
s51：根据邻域中心序列中的邻域中心点，计算邻域全局方差，得到全局散度计算公式：
[0038][0039]
其中，pi＝w
t
mi为邻域中心点的低维投影；下标g为英文global的缩写；
[0040]
s52：对式(3)进行变换可得：
[0041][0042]
其中，sg为全局散度矩阵：
[0043][0044]
s53：对式(12)进行化简，可得：
[0045]
[0046]
其中，为所有邻域中心的中心点，即全局邻域中心。
[0047]
进一步地，s6的具体步骤包括：
[0048]
s61：将公式(9)得到的局部散度和公式(11)得到的全局散度组合得到公式(4)：
[0049][0050]
s62：将式(4)转化为广义特征值求解问题，通过求解式(4)得到投影矩阵：
[0051]
sgω＝λs
l
ω
ꢀꢀ
(5)
[0052]
其中：ω为广义特征向量，λ为广义特征值；
[0053]
s63：求解式(5)得到前d个最大的特征值λ1,λ2,
…
,λd(d＜p)、对应的特征向量ω1,ω2,
…
,ωd以及投影矩阵w＝(ω1,ω2,
…
,ωd)。
[0054]
本发明与现有技术相比，其有益效果在于：
[0055]
本发明首先提出一种特征序列提取方法，提取多元时间序列协方差矩阵的上三角元素，将其组合为特征序列。然后以“局部散度最小、全局散度最大”为基本思想，提出了一种无监督降维模型，在保持局部近邻关系的同时，尽可能保留全局信息。将特征序列作为输入，最小化所有样本点邻域方差之和，最大化邻域中心点方差。求解模型得到的投影矩阵能够实现多元时间序列的降维；最终本发明通过20组公开数据集对本发明提出的降维方法及相关对比方法进行实验验证，实验结果表明，本发明所提出的降维方法能够有效对mts数据集进行降维，从而提高故障监测准确度。
具体实施方式
[0056]
下面将对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0057]
本发明公开了一种基于全局-局部散度的多元时间序列无监督降维方法，具体包括以下步骤：
[0058]
s1、通过传感器采集故障监测数，获取相应的多元时间序列，将其组成多元时间序列原始数据集d，将多元时间序列原始数据集d转化为等长特征序列集fea＝{fi|i＝1,2,
…
,n}；
[0059]
具体地，s1的具体操作步骤包括：
[0060]
s11、将获取的多元时间序列，组成多元时间序列原始数据集d＝{xi|i＝1,2,
…
,n}，其中n为样本个数，表示数据集中第i个mts样本，xi(i＝1,2,
…
,m)表示第i个变量的一系列观测值，m是变量数，ti是第i个mts的时间长度，对每个mts进行零均值化处理，即xi＝x
i-e(xi)；
[0061]
s12、对每个零均值化处理后的多元时间序列进行协方差矩阵计算，且第i个mts的协方差矩阵为：
[0062][0063]
其中，协方差矩阵为对称阵；
[0064]
将σi的上三角阵元素提取出来，按顺序组成行向量：
[0065]fi
＝(a
11
,a
12
,
…
,a
1m
,a
21
,
…
,a
2m
,
…
,a
mm
)
ꢀꢀ
(8)
[0066]
将fi作为第i个mts的特征序列，得到mts的特征集fea＝{fi|i＝1,2,
…
,n}。
[0067]
s2、基于得到的特征集fea，使用k近邻和欧氏距离(eucliddistance，ed)进行度量针对特征集得到特征集fea中每个样本的邻域集合nk(fi)＝{fj|j＝1,2,
…
,k}
[0068]
s3、根据步骤s2得到的邻域集合，计算特征集fea中每个邻域集合的邻域中心序列mi，mi的计算公式为：
[0069][0070]
其中，k为近邻数，mi为fi与其k个近邻点的邻域中心序列；
[0071]
s4、根据步骤s2得到的邻域集合，得到投影后样本点每个邻域集合的方差，再对这些方差累加求和，得到局部散度计算公式：
[0072][0073]
其中，pi＝w
t
mi为邻域中心点的低维投影；yj是指样本点的低维投影，w为投影矩阵，下标l为英文local(局部)的缩写；
[0074]
具体地，s4还包括：
[0075]
s41、对式(2)进行变换可得：
[0076][0077]
其中，s
li
为第i个特征序列fi的局部散度矩阵，即：
[0078][0079]
s5、根据步骤s3中得到的邻域中心点，计算邻域全局方差，得到全局散度：
[0080][0081]
其中，pi＝w
t
mi为邻域中心点的低维投影；下标g为英文global(全局)的缩写；
[0082]
具体地，s5还包括：
[0083]
s51、对式(3)进行变换可得
[0084][0085]
其中
[0086][0087]
s52、式(12)中sg为全局散度矩阵，对其进行化简可得：
[0088][0089]
其中，为所有邻域中心的中心点，称为全局邻域中心；
[0090]
s6、将步骤s4和s5中得到的局部散度(公式(9))和全局散度(公式(11))组合得到式(4)，通过求解式(4)得到投影矩阵w＝(ω1,ω2,
…
,ωd)：
[0091][0092]
其中，将式(4)转化为广义特征值求解问题：
[0093]
sgω＝λs
l
ω
ꢀꢀ
(5)
[0094]
其中，ω为广义特征向量，λ为广义特征值。
[0095]
通过求解式(5)，可以得到前d个最大的特征值λ1,λ2,
…
,λd(d＜p)以及对应的特征向量ω1,ω2,
…
,ωd和投影矩阵w＝(ω1,ω2,
…
,ωd)；
[0096]
s7、根据步骤s6中得到的投影矩阵，对步骤s1中得到的特征序列集中的样本点(即特征序列集fea中的元素)进行投影，得到降维特征序列：
[0097]
yi＝w
tfi
ꢀꢀ
(6)
[0098]
最终，得到降维特征集d
′
＝{yi|i＝1,2,
…
,n}，其中，
[0099]
实施例
[0100]
为了验证本发明所提降维方法(以下简称glsup)，进行相关实验。
[0101]
1、数据集选择
[0102]
设mts数据集d＝{xi|i＝1,2,
…
,n},其中n为样本个数，是数据集中第i个mts，xi(i＝1,2,
…
,m)是第i个变量的一系列观测值，m是变量数，特征序列长度为(m+1)m/2。假设数据集中mts的平均长度为t，选取的低维特征数为d(d＜＜m2)。
[0103]
2、算法复杂度
[0104]
本发明的训练过程主要包括三个部分：特征提取、模型求解以及投影。特征提取阶段，计算代价主要是mts协方差矩阵的计算，计算复杂度为o(nm2t)。模型求解阶段，计算代价主要是局部散度矩阵和广义特征值求解的计算，计算复杂度分别为o(nm4)和o(m6)。在求得投影矩阵w后，glsup方法对特征集fea中的特征序列进行投影，计算复杂度为o(dm2)。因此，glsup方法的复杂度为o(nm2t+nm4+m6)。
[0105]
3、mts数据降维实验过程
[0106]
本实施例选取了20个不同领域的多元时间序列数据集，见表1，分别为:lp1、lp2、lp3、lp4、lp5、daily and sportsactivities(dsa)、fingermovements(fm)、handmovementdirection(hmd)、natops、cricket、racketsports(rs)、epilepsy、basicmotions(bm)、lsst、articularywordrecognition(awr)、eegeye、wafer、walkvsrun(wr)、kickvspunch(kp)、australian sign language(asl)。前15个为等长数据集，后5个为不等长数据集。
[0107]
表1 mts数据集信息
[0108][0109]
降维有效性是指数据降维后，保留信息对mts特征的刻画程度。实验将降维结果输入knn(k＝1)分类器中，通过分类精度来评估降维有效性。因此实验过程可以描述为：使用
降维算法对原始mts数据集进行降维，得到降维数据集。从降维数据集中依次选取样本输入到分类器中，使用最近邻查询得到1个与被查询样本最相似的样本，将该样本标签值作为被查询样本所属类别，若与被查询样本标签值一致，则为正确分类，反之则为不正确分类。在对所有样本执行完操作后，得到分类精度：
[0110]
ε＝n
true
/n
[0111]
其中，ε为分类精度，n
true
为正确分类的样本数量，n为样本个数。
[0112]
选取pca、cpca、lpp、pblda、vpca等5种无监督降维方法作为对比方法。由于vpca方法只适用于等长数据集，只在15组等长数据集上对其进行实验。在不等长mts数据集上，pca、cpca方法的降维结果仍为不等长序列。欧氏距离只能度量等长序列。因此，在5组不等长数据集中，使用动态时间弯曲(dynamic time warping,dtw)距离对pca、cpca方法的降维结果进行度量，实现knn分类。
[0113]
pca、cpca、vpca方法涉及的参数有：方差贡献率σ；lpp方法涉及的参数有：近邻数k、热核参数t、低维特征数d；pblda方法涉及的参数有：特征维度pc、时间维度pr；glsup方法涉及的参数有：近邻数k，低维特征数d。实验中，pca、cpca和vpca方法的方差贡献率σ设为80％。近邻数k和热核参数t均设为1。时间维度pr与原时间序列长度保持不变，通过调整特征维度pc和低维特征数d来得到pblda、glsup、lpp三种降维算法的最优匹配精度。
[0114]
降维有效性实验结果如表2所示，表中每行的最高分类精度加粗显示。从结果来看，glsup方法在20个数据集中均具有较好的分类精度。glsup方法将mts转化为等长特征序列，保留了不同变量之间的相关性信息，之后同时考虑了数据集的全局和局部信息，对特征序列实现了降维。
[0115]
表2降维有效性实验结果
[0116][0117]
[0118]
4、实验结论
[0119]
对于其他几种方法，cpca方法相比于pca方法在降维有效性上有较大提升。原因在于前者将mts投影到公共低维子空间，后者则投影到不同低维子空间。但是，这两种方法仅对变量维度降维，序列长度没有缩减。vpca方法分类精度较高，但只能用于等长数据集。pblda方法从变量维度和时间维度进行降维，但在部分数据集上降维效果不佳，原因在于面对不等长数据集时，pblda方法采用截断方式将不等长序列转化为等长序列，造成信息损失。lpp方法从mts中提取特征序列并进行lpp降维，解决了不等长问题，但只考虑了数据集局部信息，忽略了全局信息。并且，该方法需要对每个mts进行奇异值分解，具有较高的计算复杂度。
[0120]
另外，在多类别数据集asl中，glsup方法相比其他方法具有明显优势。原因在于glsup方法考虑了数据集全局和局部信息，在无标签数据集中也能将样本投影形成可分性较好的簇。
[0121]
综上所述，本发明所提降维方法使用mts的协方差矩阵作为mts的特征序列，可以将时间长度不等的mts转化为等长特征序列；而且将等长特征序列投影到同一个共同的低维空间，获得的低维投影序列可代表原始mts，实现了较为明显的降维效果。
[0122]
以上实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在发明的保护范围之内。本发明的保护范围以权利要求书为准。

技术特征：
1.一种基于全局-局部散度的多元时间序列无监督降维方法，其特征在于，包括以下步骤：s1：通过传感器获取故障监测数据，将获取的故障监测数据多元时间序列组成多元时间序列原始数据集d，分别计算数据集中每个多元时间序列的协方差矩阵，提取协方差矩阵的上三角元素，将其组合为特征序列；并将所有特征序列构成特征序列集，且每一个特征序列的长度相同；s2：基于得到的特征序列集，使用k近邻数和欧氏距离进行度量建立特征序列集中每个样本邻域集合；s3：根据步骤s2得到的邻域集合，计算特征集中每个邻域集合的邻域中心序列；s4：根据步骤s2得到的邻域集合，计算待投影后的样本点每个邻域集合的方差，再对这些方差累加求和，计算局部散度；s5：根据步骤s3中得到的邻域中心点，计算邻域全局方差，得到全局散度；s6：根据步骤s4和s5中得到的局部散度和全局散度，求解投影矩阵；s7：根据步骤s6中得到的投影矩阵，对步骤s1中得到的特征序列集进行投影，从而得到降维特征序列：y
i
＝w
t
f
i
ꢀꢀꢀꢀꢀꢀꢀ
(6)；s8：根据降维特征序列，得到故障监测数据的降维特征集d
′
＝{y
i
|i＝1,2,
…
,n}，其中，s9：对降维后的故障监测数据进行处理，得到故障监测结果。2.根据权利要求1所述的一种基于全局-局部散度的多元时间序列无监督降维方法，其特征在于，步骤1的具体操作步骤包括：s11：将获取的故障监测数据多元时间序列，组成多元时间序列原始数据集d＝{x
i
|i＝1,2,
…
,n}，其中n为样本个数，表示原始数据集中第i个mts样本，x
i
(i＝1,2,
…
,m)表示第i个变量的一系列观测值，m为变量数，t
i
是第i个多元时间序列的时间长度；并对每个多元时间序列进行零均值化处理：x
i
＝x
i-e(x
i
)；s12：对每个零均值化处理后的多元时间序列进行协方差矩阵计算，且第i个mts的协方差矩阵为：s13：将协方差矩阵的上三角阵元素进行提取，并按顺序组成行向量：f
i
＝(a
11
,a
12
,
…
,a
1m
,a
21
,
…
,a
2m
,
…
,a
mm
)
ꢀꢀꢀꢀꢀꢀꢀ
(8)将f
i
作为第i个mts的特征序列，得到多元时间序列的特征集fea＝{f
i
|i＝1,2,
…
,n}。3.根据权利要求2所述一种基于全局-局部散度的多元时间序列无监督降维方法，其特征在于：步骤s3中计算邻域中心序列m
i
的计算公式为：
其中，k为近邻数，m
i
为f
i
与其k个近邻点的邻域中心序列。4.根据权利要求3所述一种基于全局-局部散度的多元时间序列无监督降维方法，其特征在于，s4的具体步骤包括：s41：基于所述邻域集合，计算投影后样本点每个邻域集合的方差，并对这些方差累加求和，得到局部散度计算公式：其中，p
i
＝w
t
m
i
为邻域中心点的低维投影；y
j
是指样本点的低维投影，w为投影矩阵，下标l为英文local的缩写；s42：对式(2)进行变换可得：其中，其中，为第i个特征序列f
i
的局部散度矩阵：5.根据权利要求4所述一种基于全局-局部散度的多元时间序列无监督降维方法，其特征在于，s5的具体步骤包括：s51：根据邻域中心序列中的邻域中心点，计算邻域全局方差，得到全局散度计算公式：其中，p
i
＝w
t
m
i
为邻域中心点的低维投影；下标g为英文global的缩写；s52：对式(3)进行变换可得：其中，s
g
为全局散度矩阵：s53：对式(12)进行化简，可得：
其中，为所有邻域中心的中心点，即全局邻域中心。6.根据权利要求5所述一种基于全局-局部散度的多元时间序列无监督降维方法，其特征在于，s6的具体步骤包括：s61：将公式(9)得到的局部散度和公式(11)得到的全局散度组合得到公式(4)：s62：将式(4)转化为广义特征值求解问题，通过求解式(4)得到投影矩阵：s
g
ω＝λs
l
ω
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)其中：ω为广义特征向量，λ为广义特征值；s63：求解式(5)得到前d个最大的特征值λ1,λ2,
…
,λ
d
(d＜p)、对应的特征向量ω1,ω2,
…
,ω
d
以及投影矩阵w＝(ω1,ω2,
…
,ω
d
)。

技术总结
本发明公开了一种基于全局-局部散度的多元时间序列无监督降维方法，其包括以下步骤：S1、计算多元时间序列的协方差矩阵,提取协方差矩阵的上三角元素,将其组合为特征序列，得到多元时间序列的特征集Fea＝{f


技术研发人员：李正欣 胡钢 刘嘉 吴虎胜 吴丹阳 刘斌 周漩 杨波 吴诗辉
受保护的技术使用者：中国人民解放军空军工程大学
技术研发日：2023.03.21
技术公布日：2023/7/12