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

Abstract

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

Description

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

技术领域

本发明涉及一种基于联合降噪和经验模态分解(EMD)的控制系统健康状态分析方法，属于控制系统信号处理与故障状态检测技术领域。

背景技术

由于现代控制系统组成复杂，常常需要在不同环境条件下长时间、高负荷工作，这就造成控制系统不可避免的会出现各类故障。特别是在航空航天、医疗、规模化机械生产等领域，细微的故障有时会造成极其严重的经济损失和人员伤害的后果，因此对设备运行的状态监测和故障诊断就成为了重要的研究课题，而保证状态监测和故障诊断准确性的前提，就是获取最能表征设备健康状况的信号特征信息。控制系统的各种状态信号，经常混杂大量的无用噪声信号，特别是在故障发生时的状态信号是典型的非平稳信号，同时由于传播介质的各向异性和多源性噪声的污染，加大了故障信号鉴别和特征提取的难度。

经验模态分解是一种非常适用于非线性、非平稳信号的分析方法。EMD把信号分解成一系列包含不同频率成分的本征模态函数(IMF)，从而使得信号的瞬时频率具有实际意义。现行通用方法对于信噪比高的信号是可行的，对干扰信号多、噪声强度大的信号，就无法分解出准确的IMF分量。而现场采集控制系统状态信号常常存在较强的随机噪声，因此在进行经验模态分解前，必须先对采样信号进行降噪处理。现阶段主要方法有：模极大值处理算法、空域相关去噪以及阀值去噪算法，其中使用最为广泛的就是小波阀值去噪算法。对于小波的选择和小波分解层数往往都是通过经验选取，同时在噪声抑制方面，只考虑了去除高斯随机噪声却没有能够很好地消除脉冲噪声的影响。

传统健康状态检测往往是设备停机后的人工检测故障方法，很难做到当控制系统刚刚出现一点故障的苗头但依然能够完成控制作用的情况下，在线实时准确地检测判断出控制系统已经处于不健康(故障)状态。这种大部分依赖操作者的工作经验来识别的方法，诊断效率低，且难以及时检测潜在故障。大量的停机检修不仅浪费大量人力与时间，还造成许多过度维修问题。

发明内容

发明目的：针对上述现有技术，提出一种基于联合降噪和经验模态分解(EMD)的控制系统健康状态分析方法，能够根据现场直接获取的控制系统状态信号，提取出能够表明信号信息的能量熵特征量且根据计算出的实时能量熵与正常状态时的熵值进行比较，在线实时判断控制系统健康状况。

技术方案：一种基于联合降噪和经验模态分解(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)的平均值；

${\mathrm{\λ}}_i=\frac{\underset{k}{\overset{N-i}{\mathrm{\Σ}}}[d_j\left(k\right)-{\stackrel{\‾}{d}}_j][d_j(k+i)-{\stackrel{\‾}{d}}_j]}{\underset{k}{\overset{m-i}{\mathrm{\Σ}}}{[d_j\left(k\right)-{\stackrel{\‾}{d}}_j]}^2},i=\mathrm{1,2},...,l---\left(1\right)$

步骤3)通过使用式(2)的改进阀值以满足：①阀值随着分解尺度的递增而逐渐减小，经过小波分解后不同分解层的系数比例分布不同；②避免出现阀值较小而起不到去除竟可能多噪声的作用。σ为噪声标准差，N为信号采样长度，j为分解尺度；

步骤4)通过使用式(3)的改进阀值函数以克服：①硬阀值函数的不连续性和引起信号的附加振荡；②软阀值函数存在恒定偏差，影响重构信号与真实信号的逼近程度。W为含噪信号小波变换后的小波系数，δ为阀值，W_δ为经过阀值降噪后的小波系数，μ(μ＞0)，v(v＞1)，p(p∈[0，1])，q(q≥0)均为可调参数；

$\mathrm{\δ}=\frac{\mathrm{\σ}\sqrt{2\ln \left(N\right)}}{\sqrt{j}}---\left(2\right)$

$W_{\mathrm{\δ}}=\left\{\begin{array}{cc}\mathrm{sgn}\left(W\right)\left(\right|W|-\frac{2p}{1+e^{q\left(\right|W|-\mathrm{\δ})}})& ,\left|W\right|\≥\mathrm{\δ}\\ W\·v^{-{\left(\mathrm{\μ}\right|W|-\mathrm{\δ})}^2}& ,\left|W\right|\≤\mathrm{\δ}\end{array}\right.---\left(3\right)$

步骤5)对处于阀值δ的[δ(1-α)，δ(1+α)](α为调节因子，调节α可以改变估计算法的精度)邻域内的小波系数通过式(4)计算尺度间相关量θ(k)，将θ(k)∈[0，β]即具有较强相关性的小波系数记为其余记为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，...，Δ)中的最小值；

$\mathrm{\θ}\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[·]为方差；

${\mathrm{\ξ}}_i=\frac{d_i}{{\mathrm{\κ}}_i}=\frac{\underset{j=1}{\overset{L}{\mathrm{\Σ}}}|\hat{s}(M_i-L+L^{\′}+j)-\hat{s}(M_0-L+L^{\′}+j)|}{\frac{\mathrm{cov}({\hat{s}}_{M_i},{\hat{s}}_{M_0})}{\sqrt{D\left[{\hat{s}}_{M_i}\right]}\·\sqrt{D\left[{\hat{s}}_{M_0}\right]}}+2}---\left(5\right)$

步骤10)对端点延拓后的信号进行EEMD分解，通过比较各本征模态函数(IMF)与信号的相关系数ρ与预设阀值ζ的大小，剔除“伪分量”。其中，ρ通过式(6)计算，σ[·]为标准差，ζ＝0.1max(ρ)；

$\mathrm{\ρ}=\frac{\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\hat{s}\left(k\right)\·\mathrm{im}{f}_i\left(k\right)}{\mathrm{\σ}\left[\hat{s}\left(k\right)\right]\·\mathrm{\σ}\left[\mathrm{im}{f}_i\left(k\right)\right]}---\left(6\right)$

步骤11)对剔除“伪分量”后的η个IMF分量，计算各IMF分量c_i(k)的能量 $\begin{array}{cc}{E}_i=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{c_i}^2\left(k\right)& (i=\mathrm{1,2},...,\mathrm{\η});\end{array}$

步骤12)通过式(7)计算能量熵；

$H=-\underset{i=1}{\overset{\mathrm{\η}}{\mathrm{\Σ}}}\frac{E_i}{\underset{i=1}{\overset{\mathrm{\η}}{\mathrm{\Σ}}}{E}_i}\·\ln \left(\frac{E_i}{\underset{i=1}{\overset{\mathrm{\η}}{\mathrm{\Σ}}}{E}_i}\right)---\left(7\right)$

步骤13)通过得到的能量熵H判断该阶段的健康状态，具体为：在检测阶段中，若上位机则告警，提示控制系统出现故障；其中，H₀为初始控制系统健康状态检测得到的能量熵值，Thre为正常状态判断阈值，Thre取值为0.1～1。

有益效果：一种基于联合降噪和经验模态分解的控制系统状态信号特征提取及健康状态实时分析方法，通过引入结合中值滤波与采用尺度间相关性的改进小波阀值联合降噪的方法，以及采用端点延拓、集合经验模态分解(EEMD)和相关系数阀值比较法相结合的方法，经过中值滤波、小波阀值降噪、经验模态分解、计算能量熵等步骤将控制系统状态特征提取出来。与初始检测的正常状态熵值比较，实时判断控制系统运行状态。具有如下具体优点：

①通过使用本发明能够对现场采集的控制系统状态信号进行较为有效地去除高斯随机噪声和脉冲噪声，且能够很好地抑制EMD分解中的端点效应与模态混叠的影响；并通过剔除“伪分量”，不仅增加了特征提取的准确性，而且减少了不必要的运算时间，有效地提取了控制系统状态信号特征；

②通过提取的能量熵的变化进行健康状况实时检测，而且具备预警潜在故障的能力。并且能够实现控制系统在工作过程中，通过实时监测到的状态数据进行健康状态评估；

③本发明提高了特征提取的针对性和准确性，实现了控制系统健康状态的自动化方式检测，对噪声干扰具有较强的鲁棒性。

本发明所提方法作为一种信号降噪、特征提取和控制系统健康状态实时检测的改进方法，具有一定的应用意义，易于实现，实时性好，能够有效提高控制系统安全性且可操作性强，节省时间，效率更高，可广泛应用于各种控制系统的特征提取与状态实时分析中。

附图说明

图1是本发明的方法的流程图；

图2是控制系统状态信号特征提取与健康状态分析框图；

图3是特征提取与健康状态分析模块框图；

图4是Quanser公司研制的用以研究四旋翼直升机控制的实验装置Qball-X4四旋翼直升机；

图5是Qball-X4控制系统状态信号特征提取与健康状态实时分析示意图；

图6是特征提取与健康状态分析模块示意图；

图7-图9是Qball-X4执行器正常、轻度故障和重度故障时的状态信号；

图10-图12是Qball-X4执行器正常、缓变故障和完全故障时的状态信号EMD分解图。

具体实施方式

下面结合附图对本发明做更进一步的解释。

如图1所示，一种基于联合降噪和经验模态分解(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)的平均值；

${\mathrm{\λ}}_i=\frac{\underset{k}{\overset{N-i}{\mathrm{\Σ}}}[d_j\left(k\right)-{\stackrel{\‾}{d}}_j][d_j(k+i)-{\stackrel{\‾}{d}}_j]}{\underset{k}{\overset{m-i}{\mathrm{\Σ}}}{[d_j\left(k\right)-{\stackrel{\‾}{d}}_j]}^2},i=\mathrm{1,2},...,l---\left(1\right)$

步骤3)通过使用式(2)的改进阀值以满足：①阀值随着分解尺度的递增而逐渐减小，经过小波分解后不同分解层的系数比例分布不同；②避免出现阀值较小而起不到去除竞可能多噪声的作用。σ为噪声标准差，N为信号采样长度，j为分解尺度；

步骤4)通过使用式(3)的改进阀值函数以克服：①硬阀值函数的不连续性和引起信号的附加振荡；②软阀值函数存在恒定偏差，影响重构信号与真实信号的逼近程度。W为含噪信号小波变换后的小波系数，δ为阀值，W_δ为经过阀值降噪后的小波系数，μ(μ＞0)，v(v＞1)，p(p∈[0，1])，q(q≥0)均为可调参数；

$\mathrm{\δ}=\frac{\mathrm{\σ}\sqrt{2\ln \left(N\right)}}{\sqrt{j}}---\left(2\right)$

$W_{\mathrm{\δ}}=\left\{\begin{array}{cc}\mathrm{sgn}\left(W\right)\left(\right|W|-\frac{2p}{1+e^{q\left(\right|W|-\mathrm{\δ})}})& ,\left|W\right|\≥\mathrm{\δ}\\ W\·v^{-{\left(\mathrm{\μ}\right|W|-\mathrm{\δ})}^2}& ,\left|W\right|\≤\mathrm{\δ}\end{array}\right.---\left(3\right)$

步骤5)对处于阀值δ的[δ(1-α)，δ(1+α)](α为调节因子，调节α可以改变估计算法的精度)邻域内的小波系数通过式(4)计算尺度间相关量θ(k)，将θ(k)∈[0，β]即具有较强相关性的小波系数记为其余记为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，...，Δ)中的最小值；

$\mathrm{\θ}\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)对处理后的各小波系数进行小波逆变换，得到降噪后的真实信号的估计 $\hat{s}\left(k\right);$

步骤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[·]为方差；

${\mathrm{\ξ}}_i=\frac{d_i}{{\mathrm{\κ}}_i}=\frac{\underset{j=1}{\overset{L}{\mathrm{\Σ}}}|\hat{s}(M_i-L+L^{\′}+j)-\hat{s}(M_0-L+L^{\′}+j)|}{\frac{\mathrm{cov}({\hat{s}}_{M_i},{\hat{s}}_{M_0})}{\sqrt{D\left[{\hat{s}}_{M_i}\right]}\·\sqrt{D\left[{\hat{s}}_{M_0}\right]}}+2}---\left(5\right)$

步骤10)对端点延拓后的信号进行EEMD分解，通过比较各本征模态函数(IMF)与信号的相关系数ρ与预设阀值ζ的大小，剔除“伪分量”。其中，ρ通过式(6)计算，σ[·]为标准差，ζ＝0.1max(ρ)；

$\mathrm{\ρ}=\frac{\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\hat{s}\left(k\right)\·\mathrm{im}{f}_i\left(k\right)}{\mathrm{\σ}\left[\hat{s}\left(k\right)\right]\·\mathrm{\σ}\left[\mathrm{im}{f}_i\left(k\right)\right]}---\left(6\right)$

步骤11)对剔除“伪分量”后的η个IMF分量，计算各IMF分量c_i(k)的能量 $\begin{array}{cc}{E}_i=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{c_i}^2\left(k\right)& (i=\mathrm{1,2},...,\mathrm{\η});\end{array}$

步骤12)通过式(7)计算能量熵；

$H=-\underset{i=1}{\overset{\mathrm{\η}}{\mathrm{\Σ}}}\frac{E_i}{\underset{i=1}{\overset{\mathrm{\η}}{\mathrm{\Σ}}}{E}_i}\·\ln \left(\frac{E_i}{\underset{i=1}{\overset{\mathrm{\η}}{\mathrm{\Σ}}}{E}_i}\right)---\left(7\right)$

步骤13)通过得到的能量熵H判断该阶段的健康状态，具体为：在检测阶段中，若上位机则告警，提示控制系统出现故障；其中，H₀为初始控制系统健康状态检测得到的能量熵值，Thre为正常状态判断阈值，Thre取值为0.1～1。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

如图2、3所示，一种基于联合降噪和经验模态分解的控制系统状态信号特征提取及健康状态实时分析方法的具体实现方式：就特征提取与健康分析模块而言，包括作为上位机的计算机1，计算机1中装有降噪模块4、特征量提取模块5和健康状态分析模块6，多通道数据采集卡2，告警装置/控制断路器3。需要对控制系统的哪一个具体部分进行健康状态监测，就将特征提取与健康实时分析装置并入控制系统的该部分旁。

下面以算例仿真和实际案例仿真说明实施方案的有效性。

1、算例仿真

构造一个调频调幅并加入噪声的非线性信号予以仿真。信号如下：

f(t)＝s(t)+n(t)＝[1+sin(50πt)]·[cos(20πt)+0.2sin(40πt)](8)

+sin(200πt)+r(t)+z(t)

该信号由调幅频率为25Hz，基频为10Hz，调频频率为20Hz，和频率为100Hz的正弦信号以及加入的高斯随机噪声r(t)、脉冲噪声z(t)混合而成，采样频率为1kHz，采样点为1024。

对信号f(t)进行中值滤波，可以发现中值滤波对脉冲噪声有很明显的抑制作用，但对随机噪声的滤波效果不理想，因此需要进一步使用小波变换滤除信号中的随机噪声。

利用Matlab进行仿真，小波函数选择db5小波，分解层数确定为5层，对比硬阀值、软阀值和本文改进阀值处理方法，对信号进行小波阀值降噪。

为了更加直观地对比各降噪处理方法的效果，引入信噪比SNR和均方根误差MSE。

SNR定义为

$\mathrm{SNR}=10\×1g\left\{\frac{\underset{N}{\mathrm{\Σ}}{\tilde{f}}^2\left(k\right)}{\underset{N}{\mathrm{\Σ}}{[\tilde{f}\left(k\right)-\hat{s}\left(k\right)]}^2}\right\}---\left(9\right)$

MSE定义为

$\mathrm{MSE}=\frac{1}{N}\sqrt{\underset{N.}{\mathrm{\Σ}}{[\tilde{f}\left(k\right)-\hat{s}\left(k\right)]}^2}---\left(10\right)$

式中，为中值滤波后的含噪信号，为小波阀值降噪后的重构信号，N为采样点数。

表1给出了采用中值滤波前后不同阀值函数对含噪信号进行降噪处理后的信噪比和均方根误差。对比数据不难发现，使用改进阀值函数的方法较其他两种传统降噪法，信噪比增加均方根误差减小；采用联合降噪后的信噪比明显提高，均方根误差明显降低。

通过本算例仿真证明，选择本发明的联合降噪方法既可以有效地消除脉冲噪声的干扰，也对高斯随机噪声有很好地抑制作用。

表1不同降噪方法性能指标比较

性能指标	SNR	MSE
			小波降噪(硬阀值函数)	3.4678	1.0578
小波降噪(软阀值函数)	2.4469	1.3381
			小波降噪(改进阀值函数)	4.1411	0.9059
联合降噪(硬阀值函数)	5.4368	0.5114
			联合降噪(软阀值函数)	3.6170	0.7776
联合降噪(改进阀值函数)	6.9831	0.3582

对联合降噪后的信号进行EMD分解，不难发现从IMF1开始就存在明显的端点效应和模态混叠问题。因此，使用本发明提出的改进端点延拓方法抑制端点效应，EEMD和相关系数预设阀值比较方法消除模态混叠与剔除“伪分量”。

对经过延拓后的信号进行EMD分解，发现端点效应被较好地抑制。再利用EEMD方法对模态混叠问题进行处理。为进一步判断“伪分量”的存在问题，计算相关系数ρ，如表2所示。比较相关系数ρ和预设阀值ζ＝0.1max(ρ)＝0.07086，因IMF7、IMF8和IMF9的相关系数均小于预设阀值，所以可认为IMF7、IMF8和IMF9均为“伪分量”并予以剔除。表3为剔除“伪分量”后的各IMF能量及特征提取出的能量熵。

表2相关系数ρ

IMF	1	2	3	4	5
						ρ	0.6328	0.5114	0.7086	0.6790	0.2191
IMF	6	7	8	9
						ρ	0.1046	0.0239	0.0571	0.0552

表3能量与能量熵

2、实际案例仿真

采用由Quanser公司研制的用以研究四旋翼直升机控制的实验装置Qball-X4四旋翼直升机飞行控制系统作为应用研究对象，Qball-X4如图4。Qball-X4四旋翼直升机，存在六维度变量即(X，Y，Z，ψ，θ，φ)，其中X，Y，Z为位置变量，ψ为偏航角，θ为俯仰角，φ为滚转角。此处以X轴方向通道信号作为研究对象，分别在正常状态、轻度故障状态和重度故障状态下，对执行器的状态信号进行特征提取。如图5、6所示，为本发明在Qball-X4四旋翼直升机X方向通道的具体实现方式。图7、8、9分别为执行器正常工作、轻度故障和重度故障时动态变量的状态信号。

使用本发明提出的联合降噪方法，对执行器状态信号进行降噪处理，然后采用改进EMD分解方法，对降噪后信号进行经验模态分解，分别得到执行器正常工作、轻度故障和重度故障时动态变量的状态信号的EMD分解图，如图10、11、12所示。

分别计算执行器正常工作、轻度故障和重度故障时各IMF分量与其各自降噪后信号的相关系数ρ，如表4、5、6所示。并比较预设阀值剔除改进EMD分解过程的“伪分量”。可以认为：在执行器正常运行过程中，不存在“伪分量”；在发生轻度故障以及处于重度故障状态时，IMF6-IMF8为“伪分量”。

表4执行器正常时的状态信号的相关系数ρ

IMF	1	2	3	4
					ρ	0.1880	0.3227	0.5291	0.0746
IMF	5	6	7	8
					ρ	0.1690	0.4168	0.9428	0.7131

表5执行器轻度故障时的状态信号的相关系数ρ

IMF	1	2	3	4
					ρ	-0.0036	0.0382	0.8573	0.8742
IMF	5	6	7	8
					ρ	0.2249	0.0511	0.0032	-0.0215

表6执行器重度故障时的状态信号的相关系数ρ

IMF	1	2	3	4
					ρ	0.0182	0.0011	0.8383	0.9765
IMF	5	6	7	8
					ρ	0.2198	0.0464	0.0497	0.0290

计算各真实IMF分量能量和能量熵，如表7、8、9。比较分析表7、8、9，不难得出：基于EMD分解的能量熵在正常状态是最大，随着故障的程度加深能量熵逐渐减小。

表7执行器正常时的状态信号的IMF分量能量和能量熵

表8执行器轻度故障时的状态信号的IMF分量能量和能量熵

表9执行器重度故障时的状态信号的IMF分量能量和能量熵

究其原因：当执行器处于正常状态时，能量分布比较均匀，所以熵值较大；而当执行器处于故障状态时，能量主要分布在故障频段，故而熵值较小。因此，通过能量熵的方法可以较为直观的判断出执行器是否处于故障状态，以及处于何种故障状态。而且，控制系统中执行器、传感器以及被控对象发生故障时的能量熵都会处在不同的区段上，因此该方法还可以推广到对整个控制系统不同组成部分的状态检测与辨识上。

当取正常时能量熵为初始值H₀，Thre＝0.1即在正常的±10％的范围，则不难发现，轻度故障和重度故障时均会告警。

Claims (1)

1.对一类含有噪声干扰的控制系统状态信号的特征提取及健康状态实时分析方法，其特征在于：通过引入结合中值滤波与采用尺度间相关性的改进小波阀值联合降噪的方法，以及采用端点延拓、集合经验模态分解(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)的平均值；

${\mathrm{\λ}}_i=\frac{\underset{k}{\overset{N-i}{\mathrm{\Σ}}}[d_j\left(k\right)-{\stackrel{\‾}{d}}_j][d_j(k+i)-{\stackrel{\‾}{d}}_j]}{\underset{k}{\overset{m-i}{\mathrm{\Σ}}}[d_j\left(k\right)-{\stackrel{\‾}{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)均为可调参数；

$\mathrm{\δ}=\frac{\mathrm{\σ}\sqrt{2\ln \left(N\right)}}{\sqrt{j}}---\left(2\right)$

$W_{\mathrm{\δ}}=\left\{\begin{array}{cc}\mathrm{sgn}\left(W\right)\left(\right|W|-\frac{2p}{1+e^{q\left(\right|W|-\mathrm{\δ})}}\mathrm{\δ})& ,\left|W\right|\≥\mathrm{\δ}\\ W\·v^{-{\left(\mathrm{\μ}\right|W|-\mathrm{\δ})}^2}& ,\left|W\right|\≤\mathrm{\δ}\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，...，Δ)中的最小值；

$\mathrm{\θ}\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[·]为方差；

${\mathrm{\ξ}}_i=\frac{d_i}{{\mathrm{\κ}}_i}=\frac{\underset{j=1}{\overset{L}{\mathrm{\Σ}}}|\hat{s}(M_i-L+L^{\′}+j)-\hat{s}(M_0-L+L^{\′}+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(ρ)；

$\mathrm{\ρ}=\frac{\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\hat{s}\left(k\right)\·\mathrm{im}{f}_i\left(k\right)}{\mathrm{\σ}\left[\hat{s}\left(k\right)\right]\·\mathrm{\σ}\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}{\mathrm{\Σ}}}{c_i}^2\left(k\right),(i=\mathrm{1,2},...,\mathrm{\η});$

步骤12)通过式(7)计算能量熵；

$H=-\underset{i=1}{\overset{\mathrm{\η}}{\mathrm{\Σ}}}\frac{E_i}{\underset{i=1}{\overset{\mathrm{\η}}{\mathrm{\Σ}}}{E}_i}\·\ln \left(\frac{E_i}{\underset{i=1}{\overset{\mathrm{\η}}{\mathrm{\Σ}}}{E}_i}\right)---\left(7\right)$

步骤13)通过得到的能量熵H判断该阶段的健康状态，具体为：在检测阶段中，若上位机则告警，提示控制系统出现故障；其中，H₀为初始控制系统健康状态检测得到的能量熵值，Thre为正常状态判断阈值，Thre取值为0.1～1。