基于联合降噪和经验模态分解的控制系统健康状态分析方法

摘要

本发明公开一种基于联合降噪和经验模态分解(EMD)的控制系统健康状态分析方法。通过引入结合中值滤波与采用尺度间相关性及改进阀值与阀值函数的小波阀值联合降噪的方法，实现了对控制系统状态信号包含的脉冲噪声和高斯随机噪声的有效抑制；针对降噪后的信号，提出采用端点延拓、集合经验模态分解(EEMD)和相关系数阀值比较相结合的方法，有效克服传统单纯依靠EMD分解在信号特征提取中存在的端点效应和模态混叠问题。本发明方法对控制系统易获得的常见状态信号进行分阶段处理，通过对计算得到的能量熵值进行分析，最终得到有效的状态特征信息，并与初始正常状态时的能量熵进行比较，可以快速实时判断出控制系统的健康状态，为控制系统的故障诊断、视情维修、容错控制等提供准确的判断依据。本发明用于高精度控制系统含噪状态信号的特征提取及健康状态实时检测。

主权项

对一类含有噪声干扰的控制系统状态信号的特征提取及健康状态实时分析方法，其特征在于：通过引入结合中值滤波与采用尺度间相关性的改进小波阀值联合降噪的方法，以及采用端点延拓、集合经验模态分解(EEMD)和相关系数阀值比较法相结合的方法，经过中值滤波、小波阀值降噪、经验模态分解、计算能量熵等步骤将控制系统状态特征提取出来。与初始检测的正常状态熵值比较，实时判断控制系统运行状态，包括如下具体步骤：步骤1)将离散含噪信号f(k)(信号长度为N)先进行中值滤波处理得到滤除可能的脉冲噪声；步骤2)优化分解层数，包括如下步骤：步骤2.1)对中值滤波处理后的进行第j层(从j＝1开始)小波分解；步骤2.2)对小波分解得到的低频近似系数a_j(k)予以保留，对高频细节系数d_j(k)进行式(1)自相关系数λ计算，若满足即各高频细节系数d_j的自相关系数满足自由度l的χ²分布，则对d_j继续进行j+1层分解；步骤2.3)直到第Δ+1层的细节系数d_Δ+1不能满足为止，确定优化分解层数为Δ层，同时得到各小波系数W_j(k)；d_j(k)高频细节系数，为d_j(k)的平均值；${\lambda }_i=\frac{\underset{k}{\overset{N-i}{\Sigma }}[d_j\left(k\right)-{\bar{d}}_j][d_j(k+i)-{\bar{d}}_j]}{\underset{k}{\overset{m-i}{\Sigma }}[d_j\left(k\right)-{\bar{d}}_j{]}^2},i=\mathrm{1,2},...,l$步骤3)通过使用式(2)的改进阀值以满足：①阀值随着分解尺度的递增而逐渐减小，经过小波分解后不同分解层的系数比例分布不同；②避免出现阀值较小而起不到去除竟可能多噪声的作用。σ为噪声标准差，N为信号采样长度，j为分解尺度；步骤4)通过使用式(3)的改进阀值函数以克服：①硬阀值函数的不连续性和引起信号的附加振荡；②软阀值函数存在恒定偏差，影响重构信号与真实信号的逼近程度。W为含噪信号小波变换后的小波系数，δ为阀值，W_δ为经过阀值降噪后的小波系数，μ(μ＞0)，v(v＞1)，p(p∈[0，1])，q(q≥0)均为可调参数；$\delta =\frac{\sigma \sqrt{2\ln \left(N\right)}}{\sqrt{j}}---\left(2\right)$$W_{\delta }=\left\{\begin{array}{cc}\mathrm{sgn}\left(W\right)\left(\right|W|-\frac{2p}{1+e^{q\left(\right|W|-\delta )}}\delta )& ,\left|W\right|\ge \delta \\ W·v^{-{\left(\mu \right|W|-\delta )}^2}& ,\left|W\right|\le \delta \end{array}\right.---\left(3\right)$步骤5)对处于阀值δ的[δ(1‑α)，δ(1+α)](α为调节因子，调节α可以改变估计算法的精度)邻域内的小波系数通过式(4)计算尺度间相关量θ(k)，将θ(k)∈[0，β]即具有较强相关性的小波系数记为(k)，其余记为(k)。W_j(k)为含噪信号f(k)在j层(j＝1，2，...，Δ)分解时的小波系数，maxW_j(k)为|W_j(k)|(j＝1，2，...，Δ)中的最大值，minW_j(k)为|W_j(k)|(j＝1，2，...，Δ)中的最小值；$\theta \left(k\right)=\ln \left(\frac{\max W_j\left(k\right)}{\min W_j\left(k\right)}\right)---\left(4\right)$步骤6)对不处于阀值δ的[δ(1‑α)，δ(1+α)]邻域内的其余所有小波系数W_j(k)进行判断，若|W_j(k)|≥δ，则将W_j(k)记为否则记为步骤7)对以|W|≥δ的情形带入式(3)加以处理，对以|W|＜δ的情形带入式(3)加以处理；步骤8)对处理后的各小波系数进行小波逆变换，得到降噪后的真实信号的估计步骤9)进行端点左右延拓(以右延拓为例)，包括如下步骤：步骤9.1)确定待延拓信号段：记的所有极大值点为{M₀，M₁，...，M_k}(从右向左标记)，所有极小值点为{m₀，m₁，...，m_k}(从右向左标记)，并记信号右端点为S；不妨假设从S向左先有极大值点后有极小值点，则记S到m₀的波形为待延拓信号段ω₀，待延拓长度为L，其中S到M₀的长度记为L′；步骤9.2)将ω₀以M₀为参考点在所有极大值点集合中依次移动，并使M₀始终与M_i重合，通过式(5)计算长度为L的波形ω₀和ω_i的匹配度系数ξ_i；步骤9.3)取ξ_i最小时的M_i为M_p，波形ω_i记为ω_p作为最佳匹配波形，则此时S所对应的坐标为X_p＝M_p+L′；步骤9.4)从X_p的后一点开始，将实际波形依次延拓到S后，一直延拓15个采样点。为S到m₀的波形为待延拓信号段ω₀与ω_i的协方差，D[·]为方差；${\xi }_i=\frac{d_i}{{\kappa }_i}=\frac{\underset{j=1}{\overset{L}{\Sigma }}|\hat{s}(M_i-L+L^{\prime }+j)-\hat{s}(M_0-L+L^{\prime }+j)|}{\frac{\mathrm{cov}(\stackrel{}{{{\hat{s}}_M}_i},{{\hat{s}}_M}_0)}{\sqrt{D[{{\hat{s}}_M}_i}]·\sqrt{D\left[{{\hat{s}}_M}_0\right]}}+2}---\left(5\right)$步骤10)对端点延拓后的信号进行EEMD分解，通过比较各本征模态函数(IMF)与信号的相关系数ρ与预设阀值ζ的大小，剔除“伪分量”。其中，ρ通过式(6)计算，σ[·]为标准差，ζ＝0.1max(ρ)；$\rho =\frac{\underset{k=1}{\overset{N}{\Sigma }}\hat{s}\left(k\right)·\mathrm{im}{f}_i\left(k\right)}{\sigma \left[\hat{s}\left(k\right)\right]·\sigma \left[\mathrm{im}{f}_i\left(k\right)\right]}---\left(6\right)$步骤11)对剔除“伪分量”后的η个IMF分量，计算各IMF分量c_i(k)的能量$E_i=\underset{k=1}{\overset{N}{\Sigma }}{c_i}^2\left(k\right),(i=\mathrm{1,2},...,\eta );$步骤12)通过式(7)计算能量熵；$H=-\underset{i=1}{\overset{\eta }{\Sigma }}\frac{E_i}{\underset{i=1}{\overset{\eta }{\Sigma }}{E}_i}·\ln \left(\frac{E_i}{\underset{i=1}{\overset{\eta }{\Sigma }}{E}_i}\right)---\left(7\right)$步骤13)通过得到的能量熵H判断该阶段的健康状态，具体为：在检测阶段中，若上位机则告警，提示控制系统出现故障；其中，H₀为初始控制系统健康状态检测得到的能量熵值，Thre为正常状态判断阈值，Thre取值为0.1～1。