CN115047853B - 基于递推规范变量残差和核主元分析的微小故障检测方法 - Google Patents

基于递推规范变量残差和核主元分析的微小故障检测方法 Download PDF

Info

Publication number
CN115047853B
CN115047853B CN202210737859.5A CN202210737859A CN115047853B CN 115047853 B CN115047853 B CN 115047853B CN 202210737859 A CN202210737859 A CN 202210737859A CN 115047853 B CN115047853 B CN 115047853B
Authority
CN
China
Prior art keywords
matrix
kernel
fault
canonical variable
vector
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210737859.5A
Other languages
English (en)
Other versions
CN115047853A (zh
Inventor
史贤俊
秦玉峰
秦亮
聂新华
王朕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Naval Aeronautical University
Original Assignee
Naval Aeronautical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Naval Aeronautical University filed Critical Naval Aeronautical University
Priority to CN202210737859.5A priority Critical patent/CN115047853B/zh
Publication of CN115047853A publication Critical patent/CN115047853A/zh
Application granted granted Critical
Publication of CN115047853B publication Critical patent/CN115047853B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00—Testing or monitoring of control systems or parts thereof
    • G05B23/02—Electric testing or monitoring
    • G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
    • G05B23/0218—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults
    • G05B23/0243—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults model based detection method, e.g. first-principles knowledge model
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B2219/00—Program-control systems
    • G05B2219/20—Pc systems
    • G05B2219/24—Pc safety
    • G05B2219/24065—Real time diagnostics
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/02—Total factory control, e.g. smart factories, flexible manufacturing systems [FMS] or integrated manufacturing systems [IMS]

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Automation & Control Theory (AREA)
  • Tests Of Electronic Circuits (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及故障诊断技术领域，尤其涉及一种基于递推规范变量残差和核主元分析的微小故障检测方法，首先构造规范变量残差,从中提取数据的线性特征，利用指数加权滑动平均法对规范变量残差进行递推滤波处理,提高规范变量残差对微小故障的敏感程度，然后使用KPCA提取规范变量残差中的非线性主成分作为非线性特征,根据提取的特征提出了两个新的故障检测统计量；此外,利用核密度估计确定故障检测统计量的控制限。由于同时提取了过程数据的线性和非线性特征,有效地提高了非线性动态过程中微小故障的可检测性。

Description

基于递推规范变量残差和核主元分析的微小故障检测方法
技术领域
本发明涉及故障诊断技术领域，尤其涉及一种基于递推规范变量残差和核主元分析的微小故障检测方法。
背景技术
随着***日益大型化和复杂化,故障诊断技术已成为保证***安全运行的一种重要手段。按照故障引起的征兆大小,可分为显著故障和微小故障。微小故障在早期的主要特征是发展变化缓慢,故障征兆不明显,容易被噪声所淹没；随着时间累积,故障幅值缓慢增加进而发展成为显著故障,如果不及早发现,可能会导致***失效从而引发严重后果。
微小故障的特点导致了在早期检测到微小故障的发生是非常困难的。现有的微小故障诊断方法主要包括基于知识的故障诊断方法、基于解析模型的故障诊断方法以及基于数据驱动的故障诊断方法。考虑到***结构组成及功能的复杂性,基于知识和基于解析模型的故障诊断方法往往难以实施。因此,基于数据驱动的微小故障诊断方法成为了研究热点。其优点在于不需要完备的结构功能等先验知识,也不需要构建精确的物理模型。大量多元统计分析技术被用于微小故障检测领域。如:Harmouche等针对传统PCA方法的T2统计量存在对微小故障不敏感的问题,提出了一种基于KL散度(Kullback-Leiblerdivergence,KLD)的微小故障检测方法；Zhang等将PCA和KLD相结合,并对PCA所得到的投影向量进行优化,使得投影向量对于KLD故障检测方法是局部最优的；Chen等基于KLD对非高斯电驱动***的早期微小故障进行了检测,并分析了该方法在较宽信噪比范围内的鲁棒性能；Cai等考虑到核主元分析模型不能敏感地检测微小故障初期的变化,引入KLD来度量核主成分的变化程度,提出了一种基于KLD-KPCA的微小故障诊断方法；陶松兵等建立了基于协方差矩阵特征值变化与KLD变化的微小故障幅值估计模型；Gautam等利用扩展卡尔曼滤波器建立故障检测指标和故障特征,然后基于KLD设计了故障决策统计量。与KLD类似,Zhang等基于JS散度(Jensen-Shannondivergence,JSD)对早期微小故障进行检测和估计。上述文献利用概率密度函数对微小故障较为敏感的特点,通过衡量故障发生前后检测数据概率密度函数之间的差异,实现对微小故障的检测。这类方法通常需要假设检测数据服从正态分布,但是实际***可能不满足该要求,因此应用范围受到一定限制。
另一类方法旨在提高残差对微小故障的灵敏度。如:Ruiz-Crcel等提出了基于规范变量分析(canonicalvariateanalysis,CVA)的微小故障检测方法,并验证了该方法的故障检测效果优于传统的PCA方法；Wu等将深度神经网络引入CVA中,利用贝叶斯推理分类器对故障进行分类；商亮亮等在CVA中引入了一阶干扰理论,显著降低了计算负荷；此外,Pilario等根据过去和未来规范变量之间的差异,构造规范变量残差,通过规范变量残差分析(canonicalvariatedissimilarityanalysis,CVDA)来处理早期微小故障检测问题；Shang等基于CVDA提出了一种加权平均统计量,提高了故障检测率；Li等基于CVDA提出了一个新的故障检测统计量,对微小故障的发展变化更敏感,同时仍然保持低的虚警率,并改进了传统的贡献图方法,提高了故障可识别性；Chen和Luo提出了一种新的多变量q-sigma规则来监测规范变量残差,并为每个变量都设置控制限,降低了检测延时和虚警率；肖姝君将CVDA推广到非线性过程,提出了一种基于核规范变量残差分析(kernelcanonicalvariatedissimilarityanalysis,KCVDA)的故障检测方法；Pilario等将不同的核函数进行组合,提出了一种基于混合KCVDA的非线性动态过程早期微小故障方法。
虽然CVDA以及KCVDA在微小故障检测方面具有一定的有效性,但单一采用一种模型并不是最佳选择:CVDA仅提取数据中的线性特征,无法提取数据的非线性特征,这些非线性特征通常出现在线性模型的残差空间中；KCVDA将原始数据映射到高维空间,从而忽略了原有空间的信息。因此,若仅用线性或非线性模型进行故障检测,微小故障的可检测性是相对较低的。
发明内容
为此，本发明提供一种基于递推规范变量残差和核主元分析的微小故障检测方法，用以克服现有技术中仅用线性或非线性模型进行故障检测,微小故障的可检测性是相对较低的的问题。
为实现上述目的，本发明提供一种基于递推规范变量残差和核主元分析的微小故障检测方法,包括：
步骤S1，离线训练，
步骤S11：获取正常运行状态下的检测数据Y0并对其进行标准化,得到标准化后的检测数据Y；
步骤S12：构造过去观测矩阵Yp和未来观测矩阵Yf，并分别计算Yp和Yf的协方差和互协方差矩阵；
步骤S13：对矩阵H执行奇异值分解并确定主导奇异值个数q；
步骤S14：计算规范变量残差d并对d进行EWMA滤波以得到滤波后的规范变量残差
Figure BDA0003714775950000031
步骤S15：构造核矩阵K并均值中心化，求解其特征值和特征向量；
步骤S16：计算主元得分向量
Figure BDA0003714775950000032
步骤S17：构造故障检测统计量
Figure BDA0003714775950000033
和Qck；
步骤S18：计算对应的控制限
Figure BDA0003714775950000034
QUCL，ck；
步骤S2，在线检测，
步骤S21：获取实际检测数据,并使用Y0的均值和协方差对其进行标准化,得到标准化后的检测数据Ytest；
步骤S22：构造过去观测矩阵Yp,test和未来观测矩阵Yf,test；
步骤S23：采用所述步骤S14的方法计算规范变量残差dtest并对dtest进行EWMA滤波以得到滤波后的规范变量残差
Figure BDA0003714775950000035
步骤S24：利用
Figure BDA0003714775950000036
和/>
Figure BDA0003714775950000037
构造核矩阵Ktest并均值中心化,求解其特征值和特征向量；/>
步骤S25：采用所述步骤S16的方法计算主元得分向量
Figure BDA0003714775950000038
步骤S26：采用所述步骤S17的方法计算故障检测统计量
Figure BDA0003714775950000039
和Qck，test；
步骤S27：若
Figure BDA0003714775950000041
或Qck，test＞QUCL，ck，则检测到故障发生。
进一步地，在所述步骤S11中，对正常运行状态下的检测数据进行标准化处理的方法为，设
Figure BDA0003714775950000042
为原始检测数据，其中，/>
Figure BDA0003714775950000043
n为样本个数，m为变量个数，/>
Figure BDA0003714775950000044
为第k个样本，对原始检测数据进行标准化处理，设定
Figure BDA0003714775950000045
式中，
Figure BDA0003714775950000046
为第j个变量的均值，sj为第j个变量的标准差，j＝1，2，。。。，m，则原始检测数据矩阵Y0转化为：Y＝[y1…yn]T。
进一步地，在所述步骤S12中，构造过去观测矩阵Yp和未来观测矩阵Yf，并分别计算Yp和Yf的协方差和互协方差矩阵的方法为，当原始检测数据矩阵Y0转化为：Y＝[y1…yn]T时，对于第k个检测样本，设定其过去观测向量yp(k)和未来观测向量yf(k)的表达式分别为：
Figure BDA0003714775950000047
Figure BDA0003714775950000048
式中，p和f分别表示过去观测向量和未来观测向量的窗口长度；
设定过去观测矩阵
Figure BDA0003714775950000049
和未来观测矩阵/>
Figure BDA00037147759500000410
的表达式分别为：
Yp＝[yp(p+1)yp(p+2)…yp(p+N)] (4)
Yf＝[yf(p+1)yf(p+2)…yf(p+N)] (5)
式中，N＝n-f-p+1，
设定Yp的协方差的表达式为：
Figure BDA0003714775950000051
设定Yf的协方差的表达式为：
Figure BDA0003714775950000052
设定Yp和Yf的互协方差的表达式为：
Figure BDA0003714775950000053
式中，参数p和f满足{mp，mf}＜N。
进一步地，在所述步骤S13中，对矩阵H执行奇异值分解并确定主导奇异值个数q的方法为，在完成所述步骤S12后，对矩阵H执行奇异值分解计算，设定：
Figure BDA0003714775950000054
式中，U和V分别是由左右奇异向量组成的矩阵，对角矩阵S由有序奇异值组成,设定S＝diag(∑1，…，∑γ，0，…，0),γ为矩阵H的秩，取U和V具有最大相关性的前q列，其中q＜mp,得到降维后的矩阵Uq和Vq，则投影矩阵J和L的表达式分别为：
Figure BDA0003714775950000055
Figure BDA0003714775950000056
进一步地，在所述步骤S14中计算规范变量残差d并对d进行EWMA滤波以得到滤波后的规范变量残差
Figure BDA0003714775950000057
的方法为，完成投影矩阵J和L的计算以使得JYp(k)和LYf(k)之间的相关性最大化,其中JYp(k)和LYf(k)为规范变量，对于第k个检测样本,设定其状态向量x(k)和残差向量e(k)的表达式分别为：
x(k)＝JYp(k) (12)
Figure BDA0003714775950000058
式中，I为适维单位阵，利用x(k)和e(k)构造如下故障检测统计量：
T2(k)＝x(k)Tx(k) (14)
Q(k)＝e(k)Te(k) (15)
式中，T2(k)统计量来度量状态向量x(k)的变化，Q(k)统计量度量残差向量e(k)的变化，则第k个检测样本的规范变量残差d(k)的表达式为：
d(k)＝Lyf(k)-SqJyp(k) (16)
式中，Sq＝diag(Σ1，…，Σq)记所有样本的规范变量残差组成的矩阵为Yd，其协方差矩阵的表达式为：
Figure BDA0003714775950000061
基于马氏距离的相关定义构造故障检测统计量D，设定：
Figure BDA0003714775950000062
采用指数加权滑动平均法EWMA对规范变量残差d进行滤波处理，滤波后的数据
Figure BDA0003714775950000063
的表达式为：
Figure BDA0003714775950000064
式中，
Figure BDA0003714775950000065
为权重因子，/>
Figure BDA0003714775950000066
取值越大，/>
Figure BDA0003714775950000067
越能反应实际数据信息，/>
Figure BDA0003714775950000068
取值越小，/>
Figure BDA0003714775950000069
对数据的微小变化越敏感。
进一步地，在所述步骤S15中构造核矩阵K并均值中心化，求解其特征值和特征向量的方法为，在得到滤波后的规范变量残差
Figure BDA00037147759500000610
后，采用KPCA算法进一步提取/>
Figure BDA00037147759500000611
中的非线性特征，假设/>
Figure BDA00037147759500000612
通过非线性函数/>
Figure BDA00037147759500000613
隐式映射到高维特征空间/>
Figure BDA00037147759500000614
且/>
Figure BDA00037147759500000615
然后求解样本协方差矩阵的特征值，设定：
Figure BDA00037147759500000616
Cξ＝λξ (21)
式中，C为空间
Figure BDA00037147759500000617
中的样本协方差矩阵，λ为特征值，ξ为特征向量且被包含于
Figure BDA0003714775950000071
所张成的子空间中,因此存在向量η＝[η1…ηN]T,使得ξ被表示为
Figure BDA0003714775950000072
的线性组合，设定：
Figure BDA0003714775950000073
则式(21)改写为：
Figure BDA0003714775950000074
式中，<·，·>表示内积，将式(20)(22)代入到式(23)中可得：
Figure BDA0003714775950000075
定义核矩阵
Figure BDA0003714775950000076
K中元素满足：
Figure BDA0003714775950000077
式中，k(·，·)为核函数,i，j＝1，…，N，选取高斯径向基函数作为核函数，则：
Figure BDA0003714775950000078
式中，h为核宽度，则公式(24)改写为：
Figure BDA0003714775950000079
即：
λNη＝Kη (28)
根据公式(28)即可确定特征向量η1，…，ηN及其对应的特征值λ1，…，λN，另外,在计算前需要对矩阵K进行均值中心化：
Figure BDA00037147759500000710
式中，
Figure BDA0003714775950000081
且其中每个元素都为1/N，取累积方差贡献率前95％的r个特征值λ1，…，λr和其对应的特征向量η1，…，ηr。
进一步地，在所述步骤S16中，计算主元得分向量
Figure BDA0003714775950000082
的方法为，在完成所述步骤S15后，对于第k个规范变量残差/>
Figure BDA0003714775950000083
通过以下投影获得其主元得分向量/>
Figure BDA0003714775950000084
设定：
Figure BDA0003714775950000085
Figure BDA0003714775950000086
式中，
Figure BDA0003714775950000087
为特征向量ηi中的第j个元素,i＝1，…，r。
进一步地，在所述步骤S17中构造故障检测统计量
Figure BDA0003714775950000088
和Qck的方法为，在完成所述步骤S16后，构造故障检测统计量/>
Figure BDA0003714775950000089
和Qck，设定：
Figure BDA00037147759500000810
Figure BDA00037147759500000811
式中，Λ＝diag(λ1，…，λr),
Figure BDA00037147759500000812
进一步地，在所述步骤S18中计算控制限的方法为，采用核密度估计KDE确定故障检测统计量的控制限，假设x为一个随机变量,p(x)为x的概率密度函数，则：
Figure BDA00037147759500000813
通过高斯核函数可得：
Figure BDA00037147759500000814
通过高斯核函数估计x的概率密度函数：
Figure BDA00037147759500000815
式中:ψ为带宽,x(i),i＝1，2，…，N为x中第i个样本，设定故障检测统计量为J,其控制限为JUCL，给定一个显著性水平α,则JUCL的计算表达式为：
Figure BDA0003714775950000091
进一步地，若J≤JUCL，则判定没有发生故障；若J＞JUCL，则判定检测到故障。
与现有技术相比，本发明的有益效果在于，本发明提出了一种基于RCVD-KPCA的微小故障检测方法，其主要贡献在于，提出了一种混合统计建模方法,同时提取过程数据的线性和非线性特征,改进了非线性动态过程的早期微小故障检测性能,提高了微小故障的可检测性。
进一步地，仿真结果表明,本发明所述方法与传统的CVDA和KCVDA方法相比,本发明所述方法所得到的故障检测统计量不仅能够更快地检测到微小故障,而且虚警率和漏检率较低,验证了本发明所述方法具有较好的故障检测性能。
附图说明
图1为本发明实施例公开的基于递推规范变量残差和核主元分析的微小故障检测方法的流程图。
具体实施方式
为了使本发明的目的和优点更加清楚明白，下面结合实施例对本发明作进一步描述；应当理解，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。
下面参照附图来描述本发明的优选实施方式。本领域技术人员应当理解的是，这些实施方式仅仅用于解释本发明的技术原理，并非在限制本发明的保护范围。
需要说明的是，在本发明的描述中，术语“上”、“下”、“左”、“右”、“内”、“外”等指示的方向或位置关系的术语是基于附图所示的方向或位置关系，这仅仅是为了便于描述，而不是指示或暗示所述装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
此外，还需要说明的是，在本发明的描述中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域技术人员而言，可根据具体情况理解上述术语在本发明中的具体含义。
请参阅图1所示，其为本发明实施例公开的基于递推规范变量残差和核主元分析的微小故障检测方法的流程图，包括：
步骤S1，离线训练，
步骤S11：获取正常运行状态下的检测数据Y0并对其进行标准化,得到标准化后的检测数据Y；
步骤S12：构造过去观测矩阵Yp和未来观测矩阵Yf，并分别计算Yp和Yf的协方差和互协方差矩阵；
步骤S13：对矩阵H执行奇异值分解并确定主导奇异值个数q；
步骤S14：计算规范变量残差d并对d进行EWMA滤波以得到滤波后的规范变量残差
Figure BDA0003714775950000101
步骤S15：构造核矩阵K并均值中心化，求解其特征值和特征向量；
步骤S16：计算主元得分向量
Figure BDA0003714775950000102
步骤S17：构造故障检测统计量
Figure BDA0003714775950000103
和Qck；/>
步骤S18：计算对应的控制限
Figure BDA0003714775950000104
QUCL，ck；
步骤S2，在线检测，
步骤S21：获取实际检测数据,并使用Y0的均值和协方差对其进行标准化,得到标准化后的检测数据Ytest；
步骤S22：构造过去观测矩阵Yp,test和未来观测矩阵Yf,test；
步骤S23：采用所述步骤S14的方法计算规范变量残差dtest并对dtest进行EWMA滤波以得到滤波后的规范变量残差
Figure BDA0003714775950000105
步骤S24：利用
Figure BDA0003714775950000111
和/>
Figure BDA0003714775950000112
构造核矩阵Ktest并均值中心化,求解其特征值和特征向量；
步骤S25：采用所述步骤S16的方法计算主元得分向量
Figure BDA0003714775950000113
步骤S26：采用所述步骤S17的方法计算故障检测统计量
Figure BDA0003714775950000114
和Qck，test；
步骤S27：若
Figure BDA0003714775950000115
或Qck，test＞QUCL，ck,则检测到故障发生。
在本实施例中，基于CVDA的规范变量残差构建：
作为一种线性降维技术,CVA能够最大程度地关联过去和未来数据集。因此,可以根据过去数据集对未来数据集的预测程度检测出数据变化。CVDA在CVA的基础上,利用过去和未来规范变量之间的差异构造出规范变量残差,通过规范变量残差检测数据变化。下面介绍CVDA的具体实现方法。
具体而言，在所述步骤S11中，对正常运行状态下的检测数据进行标准化处理的方法为，设
Figure BDA0003714775950000116
为原始检测数据,其中，/>
Figure BDA0003714775950000117
n为样本个数，m为变量个数,/>
Figure BDA0003714775950000118
为第k个样本，对原始检测数据进行标准化处理，设定
Figure BDA0003714775950000119
式中，
Figure BDA00037147759500001110
为第j个变量的均值,sj为第j个变量的标准差,j＝1,2，。。。,m，则原始检测数据矩阵Y0转化为：Y＝[y1…yn]T。
具体而言，在所述步骤S12中，构造过去观测矩阵Yp和未来观测矩阵Yf，并分别计算Yp和Yf的协方差和互协方差矩阵的方法为，当原始检测数据矩阵Y0转化为：Y＝[y1…yn]T时，对于第k个检测样本,设定其过去观测向量yp(k)和未来观测向量yf(k)的表达式分别为：
Figure BDA00037147759500001111
/>
Figure BDA0003714775950000121
式中，p和f分别表示过去观测向量和未来观测向量的窗口长度；
设定过去观测矩阵
Figure BDA0003714775950000122
和未来观测矩阵/>
Figure BDA0003714775950000123
的表达式分别为：
Yp＝[yp(p+1)yp(p+2)…yp(p+N)] (4)
Yf＝[yf(p+1)yf(p+2)…yf(p+N)] (5)
式中，N＝n-f-p+1，
设定Yp的协方差的表达式为：
Figure BDA0003714775950000124
设定Yf的协方差的表达式为：
Figure BDA0003714775950000125
设定Yp和Yf的互协方差的表达式为：
Figure BDA0003714775950000126
式中，为了避免Σpp和Σff是奇异的，参数p和f需要满足：{mp，mf}＜N。
具体而言，在所述步骤S13中，对矩阵H执行奇异值分解并确定主导奇异值个数q的方法为，在完成所述步骤S12后，对矩阵H执行奇异值分解计算，设定：
Figure BDA0003714775950000127
式中，U和V分别是由左右奇异向量组成的矩阵，对角矩阵S由有序奇异值组成，设定S＝diag(∑1，…，∑γ，0，…，0)，γ为矩阵H的秩，由于只有q(q＜mp)个主导奇异值描述***的动态特性，因此取U和V具有最大相关性的前q列，得到降维后的矩阵Uq和Vq，则投影矩阵J和L的表达式分别为：
Figure BDA0003714775950000131
Figure BDA0003714775950000132
具体而言，在所述步骤S14中计算规范变量残差d并对d进行EWMA滤波以得到滤波后的规范变量残差
Figure BDA0003714775950000133
的方法为，CVA的目标是得到投影矩阵J和L，使得JYp(k)和LYf(k)之间的相关性最大化，其中JYp(k)和LYf(k)称为规范变量，对于第k个检测样本，设定其状态向量x(k)和残差向量e(k)的表达式分别为：
x(k)＝JYp(k) (12)
Figure BDA0003714775950000134
式中，I为适维单位阵，利用x(k)和e(k)构造如下故障检测统计量：
T2(k)＝x(k)Tx(k) (14)
Q(k)＝e(k)Te(k) (15)
式中，T2(k)统计量来度量状态向量x(k)的变化，Q(k)统计量度量残差向量e(k)的变化，过去观测向量对未来观测向量的可预测性可以有效地检测数据的微小变化，因此第k个检测样本的规范变量残差d(k)的表达式为：
d(k)＝Lyf(k)-SqJyp(k) (16)
式中，Sq＝diag(∑1，…，∑q)记所有样本的规范变量残差组成的矩阵为Yd，其协方差矩阵的表达式为：
Figure BDA0003714775950000135
基于马氏距离的相关定义构造故障检测统计量D，设定：
Figure BDA0003714775950000136
现有技术中证明了统计量D(k)对微小故障检测的有效性，但它只能评估过程数据中线性特征的变化。由于线性模型的残差通常具有非线性特征，其影响不能与其他不确定性相分离，使得模型具有更高的控制限，从而降低了微小故障的可检测性，为了提高故障可检测性，有必要进一步提取规范变量残差d中的非线性特征。考虑到目前核方法技术成熟，因此本实施例应用KPCA方法实现非线性特征提取。
在本实施例中，基于RCVD-KPCA的故障检测：
为了提高规范变量残差d对数据变化的敏感程度,首先采用指数加权滑动平均法(exponentiallyweightedmovingaverage,EWMA)对规范变量残差d进行滤波处理，EWMA是工程***过程测量中一种常用的数据处理方法,其求解过程实际上是一个递推过程，滤波后的数据
Figure BDA0003714775950000141
的表达式为：
Figure BDA0003714775950000142
式中，
Figure BDA0003714775950000143
为权重因子,/>
Figure BDA0003714775950000144
取值越大,/>
Figure BDA0003714775950000145
越能反应实际数据信息，/>
Figure BDA0003714775950000146
取值越小,/>
Figure BDA0003714775950000147
对数据的微小变化越敏感。
具体而言，在所述步骤S15中构造核矩阵K并均值中心化，求解其特征值和特征向量的方法为，在得到滤波后的规范变量残差
Figure BDA0003714775950000148
后，采用KPCA算法进一步提取/>
Figure BDA0003714775950000149
中的非线性特征，假设/>
Figure BDA00037147759500001410
通过非线性函数/>
Figure BDA00037147759500001411
隐式映射到高维特征空间/>
Figure BDA00037147759500001412
且/>
Figure BDA00037147759500001413
然后求解样本协方差矩阵的特征值，设定：
Figure BDA00037147759500001414
Cξ＝λξ (21)
式中，C为空间
Figure BDA00037147759500001415
中的样本协方差矩阵,λ为特征值,ξ为特征向量，由于/>
Figure BDA00037147759500001416
无法显式表示,因此公式(21)不能够直接进行求解，考虑到ξ被包含于/>
Figure BDA00037147759500001417
所张成的子空间中,因此存在向量：η＝[η1…ηN]T,使得ξ被表示为/>
Figure BDA00037147759500001418
的线性组合，设定：
Figure BDA00037147759500001419
则式(21)改写为：
Figure BDA00037147759500001420
式中，<·，·>表示内积，将式(20)(22)代入到式(23)中可得：
Figure BDA0003714775950000151
定义核矩阵
Figure BDA0003714775950000152
K中元素满足：
Figure BDA0003714775950000153
式中，k(·，·)为核函数,i，j＝1，…，N，本实施例选取高斯径向基函数作为核函数，则：
Figure BDA0003714775950000154
式中，h为核宽度，则公式(24)改写为：
Figure BDA0003714775950000155
即：
λNη＝Kη (28)
根据公式(28)即可确定特征向量η1，…，ηN及其对应的特征值λ1，…，λN，另外,在计算前需要对矩阵K进行均值中心化：
Figure BDA0003714775950000156
式中，
Figure BDA0003714775950000157
且其中每个元素都为1/N，取累积方差贡献率前95％的r个特征值λ1，…，λr和其对应的特征向量η1，…，ηr。
具体而言，在所述步骤S16中，计算主元得分向量
Figure BDA0003714775950000158
的方法为，在完成所述步骤S15后，对于第k个规范变量残差/>
Figure BDA0003714775950000159
通过以下投影获得其主元得分向量/>
Figure BDA00037147759500001510
设定：
Figure BDA00037147759500001511
/>
Figure BDA0003714775950000161
式中，
Figure BDA0003714775950000162
为特征向量ηi中的第j个元素,i＝1，…，r。
具体而言，在所述步骤S17中构造故障检测统计量
Figure BDA0003714775950000163
和Qck的方法为，在完成所述步骤S16后，构造故障检测统计量/>
Figure BDA0003714775950000164
和Qck，设定：
Figure BDA0003714775950000165
Figure BDA0003714775950000166
式中，Λ＝diag(λ1，…，λr),
Figure BDA0003714775950000167
在本实施例中，基于核密度估计的控制限设计：
核密度估计(kerneldensityestimation,KDE)是一种确定控制上限的常用方法,尤其适用于非线性或非高斯分布过程数据。考虑到实际检测数据不一定服从正态分布,本实施例利用KDE确定故障检测统计量的控制限。
具体而言，在所述步骤S18中计算控制限的方法为，假设x为一个随机变量,p(x)为x的概率密度函数，则：
Figure BDA0003714775950000168
通过高斯核函数可得：
Figure BDA0003714775950000169
通过高斯核函数估计x的概率密度函数：
Figure BDA00037147759500001610
式中:ψ为带宽,x(i),i＝1，2，…，N为x中第i个样本，设定故障检测统计量为J,其控制限为JUCL，给定一个显著性水平α,则JUCL的计算表达式为：
Figure BDA00037147759500001611
若J≤JUCL，则判定没有发生故障；若J＞JUCL，则判定检测到故障。
根据公式(37)计算对应的控制限
Figure BDA0003714775950000171
和QUCL，ck，在所述步骤S27中，采用所述步骤S17的方法计算故障检测统计量/>
Figure BDA0003714775950000172
和Qck，test，
若
Figure BDA0003714775950000173
或Qck，test＞QUCL，ck，则检测到故障发生。
至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征做出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。
以上所述仅为本发明的优选实施例，并不用于限制本发明；对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于递推规范变量残差和核主元分析的微小故障检测方法，其特征在于，包括：
步骤S1，离线训练，
步骤S11：获取正常运行状态下的检测数据Y0并对其进行标准化,得到标准化后的检测数据Y；
步骤S12：构造过去观测矩阵Yp和未来观测矩阵Yf，并分别计算Yp和Yf的协方差和互协方差矩阵；
步骤S13：对矩阵H执行奇异值分解并确定主导奇异值个数q；
步骤S14：计算规范变量残差d并对d进行EWMA滤波以得到滤波后的规范变量残差
Figure FDA0004143176050000011
步骤S15：构造核矩阵K并均值中心化，求解其特征值和特征向量；
步骤S16：计算主元得分向量
Figure FDA0004143176050000012
步骤S17：构造故障检测统计量
Figure FDA0004143176050000013
和Qck；
步骤S18：计算对应的控制限
Figure FDA0004143176050000014
QUCL，ck；
步骤S2，在线检测，
步骤S21：获取实际检测数据,并使用Y0的均值和协方差对其进行标准化,得到标准化后的检测数据Ytest；
步骤S22：构造过去观测矩阵Yp,test和未来观测矩阵Yf,test；
步骤S23：采用所述步骤S14的方法计算规范变量残差dtest并对dtest进行EWMA滤波以得到滤波后的规范变量残差
Figure FDA0004143176050000015
步骤S24：利用
Figure FDA0004143176050000016
和/>
Figure FDA0004143176050000017
构造核矩阵Ktest并均值中心化,求解其特征值和特征向量；
步骤S25：采用所述步骤S16的方法计算主元得分向量
Figure FDA0004143176050000018
步骤S26：采用所述步骤S17的方法计算故障检测统计量
Figure FDA0004143176050000019
和Qck，test；
步骤S27：若
Figure FDA0004143176050000021
或Qck，test＞QUCL，ck,则检测到故障发生；
在所述步骤S17中构造故障检测统计量
Figure FDA0004143176050000022
和Qck的方法为，在完成所述步骤S16后，构造故障检测统计量/>
Figure FDA0004143176050000023
和Qck，设定：
Figure FDA0004143176050000024
Figure FDA0004143176050000025
式中，Λ＝diag(λ1，…，λr),
Figure FDA0004143176050000026
在所述步骤S18中计算控制限的方法为，采用核密度估计KDE确定故障检测统计量的控制限，假设x为一个随机变量,p(x)为x的概率密度函数，则：
Figure FDA0004143176050000027
通过高斯核函数可得：
Figure FDA0004143176050000028
通过高斯核函数估计x的概率密度函数：
Figure FDA0004143176050000029
式中:ψ为带宽,x(i),i＝1，2，…，N为x中第i个样本，设定故障检测统计量为J,其控制限为JUCL，给定一个显著性水平α,则JUCL的计算表达式为：
Figure FDA00041431760500000210
2.根据权利要求1所述的基于递推规范变量残差和核主元分析的微小故障检测方法，其特征在于，在所述步骤S11中，对正常运行状态下的检测数据进行标准化处理的方法为，设
Figure FDA00041431760500000211
为原始检测数据,其中，/>
Figure FDA00041431760500000212
n为样本个数，m为变量个数,/>
Figure FDA00041431760500000213
为第k个样本，对原始检测数据进行标准化处理，设定
Figure FDA0004143176050000031
式中，
Figure FDA0004143176050000032
为第j个变量的均值,sj为第j个变量的标准差,j＝1,2，。。。,m，则原始检测数据矩阵Y0转化为：Y＝[y1…yn]T。
3.根据权利要求2所述的基于递推规范变量残差和核主元分析的微小故障检测方法，其特征在于，在所述步骤S12中，构造过去观测矩阵Yp和未来观测矩阵Yf，并分别计算Yp和Yf的协方差和互协方差矩阵的方法为，当原始检测数据矩阵Y0转化为：Y＝[y1…yn]T时，对于第k个检测样本,设定其过去观测向量yp(k)和未来观测向量yf(k)的表达式分别为：
Figure FDA0004143176050000033
Figure FDA0004143176050000034
式中，p和f分别表示过去观测向量和未来观测向量的窗口长度；
设定过去观测矩阵
Figure FDA0004143176050000035
和未来观测矩阵/>
Figure FDA0004143176050000036
的表达式分别为：
Yp＝[yp(p+1) yp(p+2) … yp(p+N)] (4)
Yf＝[yf(p+1) yf(p+2) … yf(p+N)] (5)
式中，N＝n-f-p+1，
设定Yp的协方差的表达式为：
Figure FDA0004143176050000037
设定Yf的协方差的表达式为：
Figure FDA0004143176050000041
设定Yp和Yf的互协方差的表达式为：
Figure FDA0004143176050000042
式中，参数p和f满足{mp，mf}＜N。
4.根据权利要求3所述的基于递推规范变量残差和核主元分析的微小故障检测方法，其特征在于，在所述步骤S13中，对矩阵H执行奇异值分解并确定主导奇异值个数q的方法为，在完成所述步骤S12后，对矩阵H执行奇异值分解计算，设定：
Figure FDA0004143176050000043
式中，U和V分别是由左右奇异向量组成的矩阵，对角矩阵S由有序奇异值组成,设定S＝diag(∑1，…，∑γ，0，…，0),γ为矩阵H的秩，取U和V具有最大相关性的前q列，其中q＜mp,得到降维后的矩阵Uq和Vq，则投影矩阵J和L的表达式分别为：
Figure FDA0004143176050000044
Figure FDA0004143176050000045
5.根据权利要求4所述的基于递推规范变量残差和核主元分析的微小故障检测方法，其特征在于，在所述步骤S14中计算规范变量残差d并对d进行EWMA滤波以得到滤波后的规范变量残差
Figure FDA0004143176050000046
的方法为，完成投影矩阵J和L的计算以使得JYp(k)和LYf(k)之间的相关性最大化,其中JYp(k)和LYf(k)为规范变量，对于第k个检测样本,设定其状态向量x(k)和残差向量e(k)的表达式分别为：
x(k)＝JYp(k) (12)
Figure FDA0004143176050000047
式中，I为适维单位阵，利用x(k)和e(k)构造如下故障检测统计量：
T2(k)＝x(k)Tx(k) (14)
Q(k)＝e(k)Te(k) (15)
式中，T2(k)统计量来度量状态向量x(k)的变化，Q(k)统计量度量残差向量e(k)的变化，则第k个检测样本的规范变量残差d(k)的表达式为：
d(k)＝Lyf(k)-SqJyp(k) (16)
式中，Sq＝diag(∑1，…，∑q)记所有样本的规范变量残差组成的矩阵为Yd，其协方差矩阵的表达式为：
Figure FDA0004143176050000051
基于马氏距离的相关定义构造故障检测统计量D，设定：
Figure FDA0004143176050000052
采用指数加权滑动平均法EWMA对规范变量残差d进行滤波处理，滤波后的数据
Figure FDA0004143176050000053
的表达式为：
Figure FDA0004143176050000054
式中，
Figure FDA0004143176050000055
为权重因子，/>
Figure FDA0004143176050000056
取值越大，/>
Figure FDA0004143176050000057
越能反应实际数据信息，/>
Figure FDA0004143176050000058
取值越小，/>
Figure FDA0004143176050000059
对数据的微小变化越敏感。
6.根据权利要求5所述的基于递推规范变量残差和核主元分析的微小故障检测方法，其特征在于，在所述步骤S15中构造核矩阵K并均值中心化，求解其特征值和特征向量的方法为，在得到滤波后的规范变量残差
Figure FDA00041431760500000510
后，采用KPCA算法进一步提取/>
Figure FDA00041431760500000511
中的非线性特征，假设/>
Figure FDA00041431760500000512
通过非线性函数/>
Figure FDA00041431760500000513
隐式映射到高维特征空间/>
Figure FDA00041431760500000514
且/>
Figure FDA00041431760500000515
然后求解样本协方差矩阵的特征值，设定：
Figure FDA00041431760500000516
Cξ＝λξ (21)
式中，C为空间
Figure FDA00041431760500000517
中的样本协方差矩阵，λ为特征值，ξ为特征向量且被包含于
Figure FDA0004143176050000061
所张成的子空间中,因此存在向量η＝[η1 … ηN]T,使得ξ被表示为
Figure FDA0004143176050000062
的线性组合，设定：
Figure FDA0004143176050000063
则式(21)改写为：
Figure FDA0004143176050000064
式中，<·，·>表示内积，将式(20)(22)代入到式(23)中可得：
Figure FDA0004143176050000065
定义核矩阵
Figure FDA0004143176050000066
K中元素满足：
Figure FDA0004143176050000067
式中，k(·，·)为核函数,i，j＝1，…，N，选取高斯径向基函数作为核函数，则：
Figure FDA0004143176050000068
/>
式中，h为核宽度，则公式(24)改写为：
Figure FDA0004143176050000069
即：
λNη＝Kη (28)
根据公式(28)即可确定特征向量η1，…，ηN及其对应的特征值λ1，…，λN，另外,在计算前需要对矩阵K进行均值中心化：
Figure FDA00041431760500000610
式中，
Figure FDA0004143176050000071
且其中每个元素都为1/N，取累积方差贡献率前95％的r个特征值λ1，…,λr和其对应的特征向量η1，…，ηr。
7.根据权利要求6所述的基于递推规范变量残差和核主元分析的微小故障检测方法，其特征在于，在所述步骤S16中，计算主元得分向量
Figure FDA0004143176050000072
的方法为，在完成所述步骤S15后，对于第k个规范变量残差/>
Figure FDA0004143176050000073
通过以下投影获得其主元得分向量/>
Figure FDA0004143176050000074
设定：
Figure FDA0004143176050000075
Figure FDA0004143176050000076
式中，
Figure FDA0004143176050000077
为特征向量ηi中的第j个元素,i＝1，…，r。
8.根据权利要求7所述的基于递推规范变量残差和核主元分析的微小故障检测方法，其特征在于，
若J≤JUCL，则判定没有发生故障；
若J＞JUCL，则判定检测到故障。
CN202210737859.5A 2022-06-27 2022-06-27 基于递推规范变量残差和核主元分析的微小故障检测方法 Active CN115047853B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210737859.5A CN115047853B (zh) 2022-06-27 2022-06-27 基于递推规范变量残差和核主元分析的微小故障检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210737859.5A CN115047853B (zh) 2022-06-27 2022-06-27 基于递推规范变量残差和核主元分析的微小故障检测方法

Publications (2)

Publication Number Publication Date
CN115047853A CN115047853A (zh) 2022-09-13
CN115047853B true CN115047853B (zh) 2023-06-06

Family

ID=83163931

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210737859.5A Active CN115047853B (zh) 2022-06-27 2022-06-27 基于递推规范变量残差和核主元分析的微小故障检测方法

Country Status (1)

Country Link
CN (1) CN115047853B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115500829A (zh) * 2022-11-24 2022-12-23 广东美赛尔细胞生物科技有限公司 一种应用于神经内科的抑郁症检测分析***

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106959397B (zh) * 2017-04-24 2019-07-19 南京航空航天大学 一种用于高铁逆变器的微小故障诊断***的设计方法
CN107544477B (zh) * 2017-10-23 2019-05-31 中国石油大学(华东) 基于核主元分析的非线性工业过程故障检测方法
CN109145256B (zh) * 2018-11-14 2022-09-16 保控(南通)物联科技有限公司 一种基于规范变量非线性主成分分析的过程监测方法
JP7115644B2 (ja) * 2020-01-24 2022-08-09 日本精工株式会社 ノイズキャンセラ、異常診断装置、およびノイズキャンセル方法
CN111610021B (zh) * 2020-06-04 2022-05-10 沈阳科网通信息技术有限公司 一种齿轮箱早期故障检测方法
CN112906158A (zh) * 2021-02-25 2021-06-04 武汉科技大学 一种基于多传感器多元数据融合的机械故障诊断方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
一种新的早期微小故障检测方法;王佳瑜;周福娜;宋洋;;科技创新与应用(第02期) *
基于DPCA和KL散度的微小故障检测方法;周伟;潘海鹏;吴平;陈亮;;传感器与微***(第03期) *

Also Published As

Publication number Publication date
CN115047853A (zh) 2022-09-13

Similar Documents

Publication Publication Date Title
Ge et al. Multivariate statistical process control: Process monitoring methods and applications
CN107632592B (zh) 基于高效递推核主元分析的非线性时变过程故障监测方法
Harmouche et al. Incipient fault detection and diagnosis based on Kullback–Leibler divergence using principal component analysis: Part I
CN108062565B (zh) 基于化工te过程的双主元-动态核主元分析故障诊断方法
CN109739214B (zh) 工业过程间歇故障的检测方法
Ammiche et al. A combined monitoring scheme with fuzzy logic filter for plant-wide Tennessee Eastman Process fault detection
CN113642754B (zh) 一种基于rf降噪自编码信息重构和时间卷积网络的复杂工业过程故障预测方法
CN112904810B (zh) 基于有效特征选择的流程工业非线性过程监测方法
Monroy et al. Fault diagnosis of a benchmark fermentation process: a comparative study of feature extraction and classification techniques
Fan et al. AutoEncoder based high-dimensional data fault detection system
CN115047853B (zh) 基于递推规范变量残差和核主元分析的微小故障检测方法
CN112000081B (zh) 基于多块信息提取和马氏距离的故障监测方法及***
Lahdhiri et al. Interval valued data driven approach for sensor fault detection of nonlinear uncertain process
CN111368428A (zh) 一种基于监控二阶统计量的传感器精度下降故障检测方法
Sun et al. Fault detection for closed-loop control systems based on parity space transformation
Wang Enhanced fault detection for nonlinear processes using modified kernel partial least squares and the statistical local approach
CN105718733B (zh) 基于模糊贴近度和粒子滤波的故障预报方法
CN116627116B (zh) 一种流程工业故障定位方法、***及电子设备
CN110244690B (zh) 一种多变量工业过程故障辨识方法及***
CN116304823A (zh) 一种垃圾焚烧过程在线诊断方法
CN114611606A (zh) 基于核混合空间投影的故障检测方法
CN114200914A (zh) 一种基于mw-occa的质量相关早期故障检测方法
Harkat et al. Uncertain Dynamic Process Monitoring Using Moving Window PCA for Interval-Valued Data.
Ma et al. Process monitoring of the pneumatic control valve using canonical variate analysis
CN117499199B (zh) 一种基于vae的信息增强解耦网络故障诊断方法及***

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant