一种基于局部流形结构与慢特征的分布式运行状态监测方法

：本发明涉及一种基于局部流形结构与慢特征的分布式运行状态监测方法。

背景技术

0、 背景技术：

1、在现代工业中，运行状态监测技术的研究与应用日益重要，尤其是在生产和设备维护中。数据驱动的运行状态监测方法凭借其高效性和准确性，逐渐成为主流。这类方法利用大量的历史数据和实时监测数据，通过机器学习和统计分析技术，自动识别和诊断设备故障。与传统基于模型的方法相比，数据驱动的方法不依赖于物理、数学模型，能更好地应对复杂系统中的不确定性。然而，数据驱动方法在处理高维度数据时，常常面临维度灾难和特征冗余等挑战，因此需要有效的特征提取与降维技术。

2、npe(neighborhood preserving embedding)和sfa(slow feature analysis)是近年来在数据降维和特征提取领域受到广泛关注的两种技术。npe主要通过保持邻域关系来实现特征的低维映射，能够有效保留数据的局部结构信息，适合处理非线性数据。sfa则关注于从数据中提取缓慢变化的特征，能够捕捉到数据的长期趋势与潜在模式。将这两种技术应用于运行状态监测领域，可以有效地提取和降低数据的维度，从而提高运行状态监测的准确性，为工业设备的智能化维护提供强有力的支持。

技术实现思路

0、 技术实现要素：

1、本发明要解决的技术问题是：传统运行状态监测方法在处理复杂工业过程变量时，尤其是当数据呈现出不同的统计特性时，现有方法往往无法有效适应多样化的变量特征，导致运行状态监测的准确性和可靠性不足。具体而言，过程变量可以根据其统计特性分为四类：高斯动态、高斯非动态、非高斯动态和非高斯非动态。然而，现有的运行状态监测技术通常对所有变量采用统一的方法，未能充分利用这些变量的独特特性，造成信息损失和性能下降。因此，本发明提出一种新颖的运行状态监测方法，通过将数据集细分为上述四类，一方面可以根据每类变量的特性灵活选择合适的故障检测技术(npe或sfa)，一方面可以求出各块的互信息，考虑块间变量的相关性，实现更精准的运行状态监测。npe能够有效保持高维数据的局部结构，而sfa则能捕捉缓慢变化的特征，通过结合这两种方法，并加入不同子块之间的互信息，本发明不仅提高了运行状态监测的准确性，还增强了对复杂动态过程的适应能力，为工业设备的稳定运行和智能化维护提供了有效支持，具体步骤如下：

2、步骤(1)：数据采集和数据预处理，具体实施过程如下：

3、步骤(1.1)：采集正常运行工况下n个采样时刻的m个传感器监测数据，得到数据集x＝[x1,x2,…,xn]t∈rn*m，其中n为数据集样本数，m为传感器数，r为实数集。

4、步骤(1.2)：将数据按照z-score标准化方法进行数据预处理，公式如下：

5、

6、其中xi是一个传感器采集到的样本，μi是样本的均值，σi为标准差。对每一列的数据进行标准化之后得到数据集xtrain∈rn*m；

7、步骤(2)：根据xtrain中不同变量是否具有动态性以及是否服从高斯分布进行子空间划分，将xtrain分为四个子块：高斯非动态块x1、非高斯非动态块x2、高斯动态块x3、非高斯动态块x4。具体步骤如下：

8、步骤(2.1)：动态、非动态子空间的划分依据变量自身的时序相关性进行，即变量有时序相关性，则划分为动态变量，反之为非动态变量。不同类型xi∈rn×1(i＝1,2,…,m)有着不同程度的时间序列相关性，采用自相关函数(acf)度量xi(t)与xi(t+kk)之间的相关性，其中kk为滞后。对于xi∈rn×1(i＝1,…,m)，滞后为kk的自相关公式为：

9、acf(kk)＝skk/s0 (2)

10、

11、当滞后kk>0时，如果acf在给定的控制范围内，则认为该变量没有时间序列相关；否则，认为变量具有时间序列相关性。对于严格无时间序列相关的变量，acf＝0。考虑统计学意义，将acf的控制限设为0.05，从而将训练数据集xtrain∈rn×m划分为动态子空间和非动态子空间。

12、步骤(2.2)：通过jarque-bera检验对单个变量是否服从高斯分布进行判断，进而完成高斯和非高斯子块的划分：对于不同的变量xi∈rn×1(i＝1,2,…,m)，其偏度为sx，峰度为kx，当sx＝0以及kx＝3时，该变量严格服从高斯分布。jb统计量越小，表示变量越服从高斯分布。此外，jb统计量服从自由度为2的χ2分布，则可查得该统计量落在当前观测值右侧尾部的概率p，考虑统计学意义，将置信度设置为0.05，当p>0.05时认为变量服从高斯分布，反之不服从高斯分布，从而将训练数据集xtrain∈rn×m划分为高斯子空间和非高斯子空间。

13、

14、步骤(2.3)：根据以上划分，将数据集xtrain分为四块：高斯非动态块x1、非高斯非动态块x2、高斯动态块x3、非高斯动态块x4。

15、步骤(3)：对于高斯非动态块x1，设置用于局部保持的最近邻数量k1，采用npe方法对x1进行降维。假设x1∈rn*m1，最终得到降维后的特征降维映射矩阵p1。具体步骤如下：

16、步骤(3.1)：对于x1的数据x1i，找到其k1个最近邻，形成邻域ni，使用这些邻域来定义图的边，构建邻接矩阵w，表示样本点之间的连接程度，邻接矩阵w中的元素按下式构建。

17、

18、步骤(3.2)：散度矩阵s的定义为：

19、

20、步骤(3.3)：npe寻找映射矩阵p，使得样本点在低维空间中的表示能够保留原始样本点之间的距离关系，构造目标函数如下：

21、

22、该目标函数旨在最小化邻域内样本点在低维空间中的距离。

23、步骤(3.4)：通过引入拉格朗日乘子λ，构造拉格朗日函数以进行优化。其形式为：

24、l(p,λ)＝j(p)+λ(ptp-i) (12)

25、其中i是单位矩阵，确保降维后样本方差不变。

26、步骤(3.5)：通过对优化问题的求导并设置为零，得到特征值问题。将p作为特征向量进行求解：

27、sp＝λp (13)

28、其中s是散度矩阵。求解这个特征值问题可以得到映射矩阵p。

29、步骤(3.6)：将原始数据x1通过映射矩阵p映射到低维空间：

30、y1＝x1p (14)其中y1是降维后的特征数据。

31、步骤(4)：构建统计量gsd1，公式如下：

32、

33、其中，x1是x1中的数据样本，p是降维后的映射矩阵，cov(y1)是降维后数据y1的协方差矩阵。

34、步骤(5)：对于非高斯非动态子块x2，按照以上npe方法，得到统计量gsd2。

35、步骤(6)：对于高斯动态子块x3，采用sfa方法进行特征提取并降维，得到慢特征矩阵wd。具体步骤如下：

36、步骤(6.1)：对于数据集x3，求出其导数数据集

37、

38、步骤(6.2)：计算原始数据集和导数数据集的协方差矩阵：

39、导数数据集协方差矩阵如下：

40、

41、原始数据集协方差矩阵：

42、

43、步骤(6.3)：求解慢特征的映射矩阵wd，为了提取慢特征，目标是最小化慢特征的变化：

44、

45、构建广义特征值问题：

46、

47、解这个特征值问题，得到慢特征映射矩阵wd和对应的特征值λd。

48、步骤(7)：构建高斯动态子块的统计量gsd3，公式如下：

49、

50、其中，x3是x3中的数据样本，wd是慢特征矩阵。

51、步骤(8)：对于非高斯动态子块x4，仍按照以上sfa方法进行特征提取，构建统计量gsd4。

52、步骤(9)：对于x1、x2、x3、x4，通过滑动时间窗分别求出它们两两之间的互信息mm12、mm13、mm14、mm23、mm24、mm34。

53、具体步骤如下：

54、步骤(9.1)：以x1块和x2块为例，x1和x2的变量数量分别为m1和m2，设定滑动窗口长度为kk，对于x1中任意一个窗口位置w，以及任意一个变量xj，提取时间窗口内的变量数据为：x1,j(w)＝{xj(w),xj(w+1),…,xj(w+kk-1)},j＝1,2,…,m1，类似地，提取x2中窗口位置w，以及变量xk中的时间窗口变量数据为x2,k(w)。

55、步骤(9.2)：求取x1,j(w)与x2,k(w)的互信息，公式如下：

56、

57、其中，p(x1,j,x2,k)是x1,j和x2,k的联合概率密度，p(x1,j)和p(x2,k)分别是x1,j和x2,k的边缘概率密度。

58、步骤(9.3)：对每一窗口位置w，计算所有可能变量对的互信息值，形成互信息矩阵：

59、

60、步骤(9.4)：对矩阵内所有元素求和，得到第w个窗口内的互信息总量mm12(w)，随后根据滑动窗口在时间序列上的移动，依次计算得到互信息总量的时间序列mm12，公式如下：

61、

62、mm12＝{mm12(1),mm12(2),…,mm12(t)} (25)

63、其中t为滑动窗口的有效位置数量。

64、步骤(9.5)：根据以上步骤，求出其它互信息mm13、mm14、mm23、mm24、mm34。

65、步骤(10)：通过lof(local outlier factor)方法，考虑各块统计量gsd1～4以及互相关信息mm12、mm13、mm14、mm23、mm24、mm34，建立综合性的统计量lof，并求出其统计阈值lofm。具体步骤如下：

66、步骤(10.1)：构造检测指标矩阵：

67、ss＝[gsd1,gsd2,gsd3,gsd4,mm12,mm13,mm14,mm23,mm24,mm34] (16)

68、步骤(10.2)：ss是ss中的一个样本，根据其k个近邻值构造其邻域空间n(ss)。定义ssf为ss的第f个近邻，d(ss,ssf)为它们的欧式距离，令fdis(s)＝d(ss,ssk)，其中ssk为ss距离最远的近邻，那么fdis(ssf)即为ssf与其最远近邻之间的欧式距离。取d(ss,ssf)与fdis(ssf)的较大值作为reach(ss,ssf)的值，则局部可达密度以及lof值为：

69、

70、步骤(10.3)通过核密度估计kde得到lof统计量的阈值lofm。

71、上述步骤(1)至步骤(10)为本发明的离线建模阶段，如下的步骤(11)至步骤(18)为本发明方法的在线运行状态监测过程：

72、步骤(11)：采集m个传感器的监测数据xt。

73、步骤(12)：将数据xt基于数据集x＝[x1,x2,…,xn]t∈rn*m的均值以及标准差进行数据预处理，得到标准化之后的数据xt。

74、步骤(13)：将数据xt对应分为四块：高斯非动态xt1、非高斯非动态xt2、高斯动态xt3、非高斯动态xt4。

75、步骤(16)：根据离线建模步骤，分别构建数据xt各块的统计量gsdt1～4。

76、步骤(17)：根据离线建模步骤，构建数据xt各块两两之间的互信息mmt12、mmt13、mmt14、mmt23、mmt24、mmt34。

77、步骤(18)：通过lof方法，建立数据xt的综合性的统计量lof，与统计阈值lofm进行比较，以判断运行状态。

78、与传统方法相比，本发明方法的优势在于：

79、将数据集按照高斯动态、高斯非动态、非高斯动态和非高斯非动态进行分块，针对不同类型变量的分别使用适当的特征提取方法，并加入了不同子块之间的互信息，从而提高了模型的针对性以及精确性。分块处理能够更好地保留每类变量的重要信息，与直接构建整个模型相比，减少了不同变量特性之间的干扰，提高了运行状态监测模型的精度和敏感性，使得运行状态监测结果更加准确和可靠。