本发明属于工业过程监测领域，具体涉及一种针对非平稳工业过程的故障监测方法和装置。
背景技术：
1、有效的过程监测是保障工业系统高效稳定运行的关键。一个完整的过程监测主要包括故障检测和故障诊断两个环节。通过及时发现故障并确定故障发生的位置，维护人员可以快速对其处理，降低发生重大事故的风险。然而，在实际的生产中，由于原料的变化、生产目标的调整以及未知的外部干扰等，工业过程往往呈现出非平稳特性，即采集到的观测变量的均值和协方差会随时间发生变化。
2、非平稳工业过程为传统的多元统计过程监测(mspm)带来两个方面的挑战，如图1所示。其一，传统的mpsm方法如主成分分析(pca),独立成分分析(ica)等仅对平稳过程有较好的故障检测效果。而对于非平稳过程，他们无法提取到正常训练数据的非平稳变化趋势，导致在线检测阶段，正常的数据常被误检为故障。而真正的故障效应很容易被非平稳趋势所掩盖，从而出现漏检。其二，在变量的非平稳变化以及故障的涂抹效应的双重影响下，很难准确评价各过程变量对故障的贡献，无法定位到故障变量。
3、针对非平稳过程，近些年已经有了一系列方法去解决故障检测时的误检、漏检以及变量冗余的干扰问题。在研究之初，一些方法通过计算非平稳时间序列的差值来消除非平稳特征的影响。然而，作差处理之后可能会对原始数据的动态信息和故障特征造成破坏，从而影响检测性能。后来，针对非平稳过程的故障检测方法主要有基于自适应策略，基于慢特征分析(sfa)，基于基于协整分析(ca)以及基于平稳子空间分析(ssa)等发展方向。基于自适应策略的方法通过模型的不断更新来获取数据的变化来解决非平稳的问题。文献[1]针对过程变化导致固定模型监控方法的误报问题，提出了递归pca的方法。基于sfa的方法通过学习数据不变或缓慢变化的特征用以故障检测。文献[2]提出基于sfa的过程监控策略能同时监测工作点偏差和工艺动力学异常。基于ca的方法试图找到非平稳变量之间的长期平衡关系，使用建立的协整模型对实时数据进行故障检测。文献[3]将协整分析应用于过程监测，通过协整模型描述装置的动态平衡关系。与ca类方法不同，基于ssa的方法不需要协整假设，对非平稳变量的积分顺序没有要求，因此可用于可靠的过程监测。ssa通过获取平稳和非平稳投影矩阵来分离平稳源和非平稳源。文献[4]用利用瓦瑟斯坦-2距离代替ssa的kullback-leibler(kl)散度提出了wasserstein ssa(wssa)。wssa获得了更高的计算效率，并对异常值更为鲁棒。
4、虽然上述方法一定程度解决了非平稳过程故障检测的系列问题，然而实现精准的故障诊断依旧是很大的挑战。故障检测和故障诊断是过程监测中密切相关且互为补充的两个重要阶段，但是关于非平稳过程的故障诊断的研究还并未广泛展开，目前只有一些初步的成果。文献[5]通过将最小绝对收缩和选择算子应用于协整分析，利用稀疏重构的策略实现在线识别故障。然而该方法依据每个变量被选择的频率确定的非平稳变量范围较为宽泛，受缺失重建比的阈值影响较大，可能导致诊断结果的不正确。文献[6]提出了基于主导趋势的逻辑回归的故障诊断方法来监测非平稳过程。该方法使用滑动窗口来收集过程数据和提取主要的趋势特征，用于区分故障并确定故障类型。虽然该方法可以实现较为准确的故障诊断效果，但是受到滑动窗口的大小和频率的影响，可能出现故障检测的延迟，而且高频的在线更新也会造成计算复杂度的提高。综上所述，目前非平稳过程依旧缺乏有效的故障诊断策略。
技术实现思路
1、针对非平稳工业过程的故障检测和诊断难题，本发明提供一种针对非平稳工业过程的故障监测方法和装置，实现高故障检测率以及精准确定故障位置。
2、为实现上述技术目的，本发明采用如下技术方案：
3、一种针对非平稳工业过程的故障监测方法，包括：
4、利用稀疏平稳投影矩阵对工业过程的多维观测变量进行稀疏线性组合构建多维度平稳源，并构建平稳源的监测统计指标；
5、其中，所述稀疏平稳投影矩阵，采用平稳子空间分析方法并在优化问题中引入l2,p范数作为稀疏约束，并利用由多维观测变量构成的训练数据求解得到；
6、根据每条训练数据，均计算平稳源的监测统计指标，并根据核密度估计计算平稳源的控制限；
7、在线监测阶段，根据实时观测数据计算平稳源的监测统计指标，并与控制限比较以完成故障检测，即判断工业过程是否发生异常；
8、再采用关键信息动态重构的方法，通过将每个观测变量对于监测统计指标的贡献值进行动态重构，通过比较重构贡献值的大小来确定故障变量，完成故障诊断。
9、进一步地，利用稀疏平稳投影矩阵对工业过程的多维观测变量进行稀疏线性组合构建平稳源的表达式为：
10、
11、其中，为t时刻的m个平稳源，为t时刻的d-m个非平稳源，x(t)为t时刻的d维观测变量，x(t)∈rd；bs和bn分别为平稳投影矩阵和非平稳投影矩阵，bs∈rm×d，bn∈r(d-m)×d。
12、进一步地，所述采用平稳子空间分析方法并在优化问题中引入l2,p范数作为稀疏约束，得到优化问题为：
13、
14、其中，bs为平稳投影矩阵；tr()表示括号内矩阵的迹；s为中间量，ui和∑i分别为第i个时期的观测训练数据的均值和协方差，和分别为所有n个时期的观测训练数据的平均均值和平均协方差，λ＞0是正则化系数，用以调节平稳性和稀疏性之间的平衡关系，p表示矩阵范数参数，右角标t表示转置运算，||·||2,p表示l2,p范数；im为m行m列的单位矩阵，m为平稳源的数目。
15、进一步地，求解优化问题以得到平稳投影矩阵具体为：
16、a1、将式(2)所示的优化问题转换为迹的形式：
17、
18、其中，d为中间量：
19、
20、其中，bj代表平稳投影矩阵bs的第j列向量；
21、a2、令各列向量：
22、
23、其中，是一个极小值；
24、a3、利用拉格朗日乘子方法来求解式(3)所示的优化问题，其中拉格朗日函数为：
25、l＝tr[bss(bs)t]+λtr[bsd(bs)t]-tr[λ(bs(bs)t-im)] (6)
26、其中，λ∈rm×m是拉格朗日乘子的矩阵；
27、a4、对bs求导得：
28、
29、a5、先保持d不变，更新平稳投影矩阵bs：
30、bs(s+λd)＝λbs (8)
31、式中，这是一个关于s+λd的特征值分解问题，bs由s+λd的m个最小特征值对应的m个特征向量给出平稳投影矩阵
32、然后根据更新的bs并按照式子(4)计算d；
33、a6、重复迭代优化步骤a5，直至收敛。
34、进一步地，采用马氏距离构建所有平稳源的整体监测统计指标：
35、
36、其中，μ和∑是所有训练数据的均值和协方差，xk为第k个采样瞬间观测得到的训练数据，xk∈rd；是训练数据xk对应的平稳源，tk表示平稳源的整体监测统计指标。
37、进一步地，故障诊断具体包括：
38、在训练阶段，针对每个平稳源单独计算监测统计指标：
39、
40、其中，j∈{1,2,...,m}，m是平稳源的数目，亦即平稳源的维度；代表第k个采样瞬间的第j个平稳源，bj代表平稳投影矩阵bs的第j行向量；μ是所有训练数据的均值；σj代表训练数据第j个平稳源的方差；表示第j个平稳源在第k个采样瞬间的监测统计指标；
41、针对每个平稳源的统计指标分别根据核密度估计计算对应的控制限；
42、在线监测阶段，根据实时观测数据xnew计算各平稳源单独的监测统计指标；若某个平稳源的监测统计指标大于对应的控制限，则认为故障变量潜在于该平稳源中；
43、根据潜有故障变量的平稳源构建实时观测数据xnew的动态选择矩阵wnew：
44、
45、其中，是潜有故障变量的h个平稳源的索引值，组成列向量为代表列向量τnew的第l个元素，是第个元素为1而其余元素为0的m维列向量；
46、基于动态选择矩阵重构监测统计量：
47、
48、其中，zi为沿着可能的故障方向ξi重构得到的变量，zi＝x-ξifi；ξi表示d维单位矩阵id的第i列，代表工业过程第i个非平稳变量故障的可能性；fi代表估计的ξi方向上故障的大小；为重构的监测统计量；中间量∑re＝(wnew)tbs∑(bs)twnew，xc＝xnew-μ；
49、求解最小时对应的fi的表达式，其最小二乘解为：
50、fi＝[mt(∑re)-1m]-1mt(∑re)-1wnewtbsxc (13)
51、其中，中间量m＝wnewtbsξi；则变量xi到监测统计量的重构贡献值为
52、
53、最终比较各个元素的大小，选择最大的所对应的非平稳变量xi为故障发生的源头。
54、本发明还提供一种基于上述故障监测方法的故障监测装置，包括：
55、平稳投影矩阵优化模块，用于：采用平稳子空间分析方法并在优化问题中引入l2,p范数作为稀疏约束，并利用由多维观测变量构成的训练数据求解得到平稳投影矩阵；
56、平稳源及指标构建模块，用于：利用平稳投影矩阵对工业过程的多维观测变量进行线性组合构建多维度平稳源，并构建平稳源的监测统计指标；
57、控制限计算模块，用于：根据每条训练数据，均计算平稳源的监测统计指标，并根据核密度估计计算平稳源的控制限；
58、在线监测模块，用于：根据实时观测数据计算平稳源的监测统计指标，并与控制限比较以完成故障检测，即判断工业过程是否发生异常；以及用于：采用关键信息动态重构的方法，通过将每个观测变量对于监测统计指标的贡献值进行动态重构，通过比较重构贡献值的大小来确定故障变量，完成故障诊断。
59、进一步地，所述采用平稳子空间分析方法并在优化问题中引入l2,p范数作为稀疏约束，得到优化问题为：
60、
61、其中，bs为平稳投影矩阵；tr()表示括号内矩阵的迹；s为中间量，ui和∑i分别为第i个时期的观测训练数据的均值和协方差，和分别为所有n个时期的观测训练数据的平均均值和平均协方差，λ＞0是正则化系数，用以调节平稳性和稀疏性之间的平衡关系，p表示矩阵范数参数，右角标t表示转置运算，||·||2,p表示l2,p范数；im为m行m列的单位矩阵，m为平稳源的数目。
62、进一步地，采用马氏距离构建所有平稳源的整体监测统计指标：
63、
64、其中，μ和∑是所有训练数据的均值和协方差，xk为第k个采样瞬间观测得到的训练数据，xk∈rd；是训练数据xk对应的平稳源，tk表示平稳源的整体监测统计指标。
65、进一步地，所述控制限计算模块还用于：在训练阶段，针对每个平稳源单独计算监测统计指标和根据核密度估计计算对应的控制限：
66、
67、其中，j∈{1,2,...,m}，m是平稳源的数目，亦即平稳源的维度；代表第k个采样瞬间的第j个平稳源，bj代表平稳投影矩阵bs的第j行向量；μ是所有训练数据的均值；σj代表训练数据第j个平稳源的方差；表示第j个平稳源在第k个采样瞬间的监测统计指标；
68、针对每个平稳源的统计指标分别根据核密度估计计算对应的控制限；
69、所述故障诊断具体包括：
70、(1)根据实时观测数据xnew计算各平稳源单独的监测统计指标；若某个平稳源的监测统计指标大于对应的控制限，则认为故障变量潜在于该平稳源中；
71、(2)根据潜有故障变量的平稳源构建实时观测数据xnew的动态选择矩阵wnew：
72、
73、其中，是潜有故障变量的h个平稳源的索引值，组成列向量为代表列向量τnew的第l个元素，是第个元素为1而其余元素为0的m维列向量；
74、(3)基于动态选择矩阵重构监测统计量：
75、
76、其中，zi为沿着可能的故障方向ξi重构得到的变量，zi＝x-ξifi；ξi表示d维单位矩阵id的第i列，代表工业过程第i个非平稳变量故障的可能性；fi代表估计的ξi方向上故障的大小；为重构的监测统计量；中间量∑re＝(wnew)tbs∑(bs)twnew，xc＝xnew-μ；
77、(4)求解最小时对应的fi的表达式，其最小二乘解为：
78、fi＝[mt(∑re)-1m]-1mt(∑re)-1wnewtbsxc (13)
79、其中，中间量m＝wnewtbsξi；则变量xi到监测统计量的重构贡献值为
80、
81、(5)最终比较各个元素的大小，选择最大的所对应的非平稳变量xi为故障发生的源头。
82、有益效果
83、针对非平稳工业过程带来的故障检测和诊断等方面的难题，本发明提出引入l2,p范数作为稀疏约束的swssa方法求解平稳投影矩阵，并基于马氏距离进行故障检测和基于cidr的故障诊断策略，有效实现对非平稳工业过程进行监测，而且故障检测率高且故障定位诊断精准。本发明发明可广泛用于现代工业过程的故障检测和诊断领域，所提的算法可以打包成可执行文件布置在工业现场，通过读取现场数据，进行实时过程监测。
1.一种针对非平稳工业过程的故障监测方法，其特征在于，包括：
2.根据权利要求1所述的针对非平稳工业过程的故障监测方法，其特征在于，利用稀疏平稳投影矩阵对工业过程的多维观测变量进行稀疏线性组合构建平稳源的表达式为：
3.根据权利要求1所述的针对非平稳工业过程的故障监测方法，其特征在于，所述采用平稳子空间分析方法并在优化问题中引入l2,p范数作为稀疏约束，得到优化问题为：
4.根据权利要求3所述的针对非平稳工业过程的故障监测方法，其特征在于，求解优化问题以得到平稳投影矩阵具体为：
5.根据权利要求3所述的针对非平稳工业过程的故障监测方法，其特征在于，采用马氏距离构建所有平稳源的整体监测统计指标：
6.根据权利要求3所述的针对非平稳工业过程的故障监测方法，其特征在于，故障诊断具体包括：
7.一种基于权利要求1所述故障监测方法的故障监测装置，其特征在于，包括：
8.根据权利要求7所述的故障监测装置，其特征在于，所述采用平稳子空间分析方法并在优化问题中引入l2,p范数作为稀疏约束，得到优化问题为：
9.根据权利要求7所述的故障监测装置，其特征在于，采用马氏距离构建所有平稳源的整体监测统计指标：
10.根据权利要求7所述的故障监测装置，其特征在于，所述控制限计算模块还用于：在训练阶段，针对每个平稳源单独计算监测统计指标和根据核密度估计计算对应的控制限：
