CN105094111B - 基于联合降噪和经验模态分解的控制***健康状态分析方法 - Google Patents
基于联合降噪和经验模态分解的控制***健康状态分析方法 Download PDFInfo
- Publication number
- CN105094111B CN105094111B CN201510172003.8A CN201510172003A CN105094111B CN 105094111 B CN105094111 B CN 105094111B CN 201510172003 A CN201510172003 A CN 201510172003A CN 105094111 B CN105094111 B CN 105094111B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- mover
- mfrac
- signal
- 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.)
- Expired - Fee Related
Links
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/0224—Process history based detection method, e.g. whereby history implies the availability of large amounts of data
- G05B23/0227—Qualitative history assessment, whereby the type of data acted upon, e.g. waveforms, images or patterns, is not relevant, e.g. rule based assessment; if-then decisions
- G05B23/0235—Qualitative history assessment, whereby the type of data acted upon, e.g. waveforms, images or patterns, is not relevant, e.g. rule based assessment; if-then decisions based on a comparison with predetermined threshold or range, e.g. "classical methods", carried out during normal operation; threshold adaptation or choice; when or how to compare with the threshold
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Automation & Control Theory (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了基于联合降噪和经验模态分解的控制***健康状态分析方法。通过引入改进阈值与阈值函数的小波阈值联合降噪的方法,实现对控制***状态信号包含的脉冲噪声和高斯随机噪声的有效抑制;针对降噪后的信号,提出端点延拓、集合经验模态分解(EEMD)和相关系数阈值比较相结合的方法,有效克服了端点效应和模态混叠问题。本发明方法对控制***常见状态信号进行分阶段处理,通过对计算得到的能量熵值进行分析,最终得到有效的状态特征信息,并与正常状态时的能量熵进行比较,快速实时判断出***的健康状态,为***的故障诊断、视情维修、容错控制等提供准确的判断依据。本发明用于高精度控制***含噪状态信号的特征提取及健康状态实时检测。
Description
技术领域
本发明涉及一种基于联合降噪和经验模态分解(EMD)的控制***健康状态分析方法,属于控制***信号处理与故障状态检测技术领域。
背景技术
由于现代控制***组成复杂,常常需要在不同环境条件下长时间、高负荷工作,这就造成控制***不可避免的会出现各类故障。特别是在航空航天、医疗、规模化机械生产等领域,细微的故障有时会造成极其严重的经济损失和人员伤害的后果,因此对设备运行的状态监测和故障诊断就成为了重要的研究课题,而保证状态监测和故障诊断准确性的前提,就是获取最能表征设备健康状况的信号特征信息。控制***的各种状态信号,经常混杂大量的无用噪声信号,特别是在故障发生时的状态信号是典型的非平稳信号,同时由于传播介质的各向异性和多源性噪声的污染,加大了故障信号鉴别和特征提取的难度。
经验模态分解是一种非常适用于非线性、非平稳信号的分析方法。EMD把信号分解成一系列包含不同频率成分的本征模态函数(IMF),从而使得信号的瞬时频率具有实际意义。现行通用方法对于信噪比高的信号是可行的,对干扰信号多、噪声强度大的信号,就无法分解出准确的IMF分量。而现场采集控制***状态信号常常存在较强的随机噪声,因此在进行经验模态分解前,必须先对采样信号进行降噪处理。现阶段主要方法有:模极大值处理算法、空域相关去噪以及阈值去噪算法,其中使用最为广泛的就是小波阈值去噪算法。对于小波的选择和小波分解层数往往都是通过经验选取,同时在噪声抑制方面,只考虑了去除高斯随机噪声却没有能够很好地消除脉冲噪声的影响。
传统健康状态检测往往是设备停机后的人工检测故障方法,很难做到当控制***刚刚出现一点故障的苗头但依然能够完成控制作用的情况下,在线实时准确地检测判断出控制***已经处于不健康(故障)状态。这种大部分依赖操作者的工作经验来识别的方法,诊断效率低,且难以及时检测潜在故障。大量的停机检修不仅浪费大量人力与时间,还造成许多过度维修问题。
发明内容
发明目的:针对上述现有技术,提出一种基于联合降噪和经验模态分解(EMD)的控制***健康状态分析方法,能够根据现场直接获取的控制***状态信号,提取出能够表明信号信息的能量熵特征量且根据计算出的实时能量熵与正常状态时的熵值进行比较,在线实时判断控制***健康状况。
技术方案:一种基于联合降噪和经验模态分解(EMD)的控制***健康状态分析方法,通过引入结合中值滤波与采用尺度间相关性的改进小波阈值联合降噪的方法,以及采用端点延拓、集合经验模态分解(EEMD)和相关系数阈值比较法相结合的方法,经过中值滤波、小波阈值降噪、经验模态分解、计算能量熵等步骤将控制***状态特征提取出来,与初始检测的正常状态熵值比较,实时判断控制***运行状态,包括如下具体步骤:
步骤1)将离散含噪信号f(k)先进行中值滤波处理得到滤除可能的脉冲噪声,其中,f(k)信号长度为N;
步骤2)优化分解层数,包括如下步骤:
步骤2.1)对中值滤波处理后的进行第j层小波分解,其中,从j=1开始分解;
步骤2.2)对小波分解得到的低频近似系数aj(k)予以保留,对高频细节系数dj(k)进行式(1)自相关系数λ计算,若满足(l)即各高频细节系数dj的自相关系数满足自由度l的χ2分布,则对dj继续进行j+1层分解;
步骤2.3)直到第Δ+1层的细节系数dΔ+1不能满足为止,确定优化分解层数为Δ层,同时得到各小波系数Wj(k);dj(k)为高频细节系数,为dj(k)的平均值;
步骤3)通过使用式(2)的改进阈值以满足:①阈值随着分解尺度的递增而逐渐减小,经过小波分解后不同分解层的系数比例分布不同;②避免出现阈值较小而起不到去除尽可能多噪声的作用;σ为噪声标准差,N为信号采样长度,j为分解尺度;
步骤4)通过使用式(3)的改进阈值函数以克服:①硬阈值函数的不连续性和引起信号的附加振荡;②软阈值函数存在恒定偏差,影响重构信号与真实信号的逼近程度;W为含噪信号小波变换后的小波系数,δ为阈值,Wδ为经过阈值降噪后的小波系数,μ,v,p,q均为可调参数,其中μ>0,v>1,p∈[0,1],q≥0;
步骤5)对处于阈值δ的[δ(1-α),δ(1+α)]邻域内的小波系数通过式(4)计算尺度间相关量θ(k),其中α为调节因子且调节α可以改变估计算法的精度;将θ(k)∈[0,β]即具有较强相关性的小波系数记为其余记为Wj(k)为含噪信号f(k)在j层分解时的小波系数,max Wj(k)为|Wj(k)|中的最大值,min Wj(k)为|Wj(k)|中的最小值,其中j=1,2,...,Δ;
步骤6)对不处于阈值δ的[δ(1-α),δ(1+α)]邻域内的其余所有小波系数Wj(k)进行判断,若|Wj(k)|≥δ,则将Wj(k)记为否则记为
步骤7)对以|W|≥δ的情形带入式(3)加以处理,对以|W|<δ的情形带入式(3)加以处理;
步骤8)对处理后的各小波系数进行小波逆变换,得到降噪后的真实信号的估计
步骤9)进行端点左右延拓,以右延拓为例进行介绍,包括如下步骤:
步骤9.1)确定待延拓信号段:记的所有极大值点为{M0,M1,...,Mk},所有极小值点为{m0,m1,...,mk},其中极大值点和极小值点分别从右向左标记,并记信号右端点为S;不妨假设从S向左先有极大值点后有极小值点,则记S到m0的波形为待延拓信号段ω0,待延拓长度为L,其中S到M0的长度记为L′;
步骤9.2)将ω0以M0为参考点在所有极大值点集合中依次移动,并使M0始终与Mi重合,通过式(5)计算长度为L的波形ω0和ωi的匹配度系数ξi;
步骤9.3)取ξi最小时的Mi为Mp,波形ωi记为ωp作为最佳匹配波形,则此时S所对应的坐标为Xp=Mp+L′;
步骤9.4)从Xp的后一点开始,将实际波形依次延拓到S后,一直延拓15个采样点,为S到m0的波形为待延拓信号段ω0与ωi的协方差,D[·]为方差;
步骤10)对端点延拓后的信号进行EEMD分解,通过比较各本征模态函数(IMF)与信号的相关系数ρ与预设阈值ζ的大小,剔除“伪分量”;其中,ρ通过式(6)计算,σ[·]为标准差,ζ=0.1max(ρ);
步骤11)对剔除“伪分量”后的η个IMF分量,计算各IMF分量ci(k)的能量其中i=1,2,…,η;
步骤12)通过式(7)计算能量熵;
步骤13)通过得到的能量熵H判断检测阶段的健康状态,具体为:在检测阶段中,若上位机则告警,提示控制***出现故障;其中,H0为初始控制***健康状态检测得到的能量熵值,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)先进行中值滤波处理得到滤除可能的脉冲噪声,其中,f(k)信号长度为N;
步骤2)优化分解层数,包括如下步骤:
步骤2.1)对中值滤波处理后的进行第j层小波分解,其中,从j=1开始分解;
步骤2.2)对小波分解得到的低频近似系数aj(k)予以保留,对高频细节系数dj(k)进行式(1)自相关系数λ计算,若满足即各高频细节系数dj的自相关系数满足自由度l的χ2分布,则对dj继续进行j+1层分解;
步骤2.3)直到第Δ+1层的细节系数dΔ+1不能满足为止,确定优化分解层数为Δ层,同时得到各小波系数Wj(k);dj(k)为高频细节系数,为dj(k)的平均值;
步骤3)通过使用式(2)的改进阈值以满足:①阈值随着分解尺度的递增而逐渐减小,经过小波分解后不同分解层的系数比例分布不同;②避免出现阈值较小而起不到去除尽可能多噪声的作用;σ为噪声标准差,N为信号采样长度,j为分解尺度;
步骤4)通过使用式(3)的改进阈值函数以克服:①硬阈值函数的不连续性和引起信号的附加振荡;②软阈值函数存在恒定偏差,影响重构信号与真实信号的逼近程度;W为含噪信号小波变换后的小波系数,δ为阈值,Wδ为经过阈值降噪后的小波系数,μ,v,p,q均为可调参数,其中μ>0,v>1,p∈[0,1],q≥0;
步骤5)对处于阈值δ的[δ(1-α),δ(1+α)]邻域内的小波系数通过式(4)计算尺度间相关量θ(k),其中α为调节因子且调节α可以改变估计算法的精度;将θ(k)∈[0,β]即具有较强相关性的小波系数记为其余记为Wj(k)为含噪信号f(k)在j层分解时的小波系数,max Wj(k)为|Wj(k)|中的最大值,min Wj(k)为|Wj(k)|中的最小值,其中j=1,2,...,Δ;
步骤6)对不处于阈值δ的[δ(1-α),δ(1+α)]邻域内的其余所有小波系数Wj(k)进行判断,若|Wj(k)|≥δ,则将Wj(k)记为否则记为
步骤7)对以|W|≥δ的情形带入式(3)加以处理,对以|W|<δ的情形带入式(3)加以处理;
步骤8)对处理后的各小波系数进行小波逆变换,得到降噪后的真实信号的估计
步骤9)进行端点左右延拓,以右延拓为例进行介绍,包括如下步骤:
步骤9.1)确定待延拓信号段:记的所有极大值点为{M0,Ml,…,Mk},所有极小值点为{m0,m1,...,mk},其中极大值点和极小值点分别从右向左标记,并记信号右端点为S;不妨假设从S向左先有极大值点后有极小值点,则记S到m0的波形为待延拓信号段ω0,待延拓长度为L,其中S到M0的长度记为L′;
步骤9.2)将ω0以M0为参考点在所有极大值点集合中依次移动,并使M0始终与Mi重合,通过式(5)计算长度为L的波形ω0和ωi的匹配度系数ξi;
步骤9.3)取ξi最小时的Mi为Mp,波形ωi记为ωp作为最佳匹配波形,则此时S所对应的坐标为Xp=Mp+L′;
步骤9.4)从Xp的后一点开始,将实际波形依次延拓到S后,一直延拓15个采样点,为S到m0的波形为待延拓信号段ω0与ωi的协方差,D[·]为方差;
步骤10)对端点延拓后的信号进行EEMD分解,通过比较各本征模态函数(IMF)与信号的相关系数ρ与预设阈值ζ的大小,剔除“伪分量”;其中,ρ通过式(6)计算,σ[·]为标准差,ζ=0.1max(ρ);
步骤11)对剔除“伪分量”后的η个IMF分量,计算各IMF分量ci(k)的能量其中i=1,2,…,η;
步骤12)通过式(7)计算能量熵;
步骤13)通过得到的能量熵H判断检测阶段的健康状态,具体为:在检测阶段中,若上位机则告警,提示控制***出现故障;其中,H0为初始控制***健康状态检测得到的能量熵值,Thre为正常状态判断阈值,Thre取值为0.1~1。
以上所述仅是本发明的优选实施方式,应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下,还可以做出若干改进和润饰,这些改进和润饰也应视为本发明的保护范围。
如图2、3所示,一种基于联合降噪和经验模态分解的控制***状态信号特征提取及健康状态实时分析方法的具体实现方式:就特征提取与健康分析模块而言,包括作为上位机的计算机1,计算机1中装有降噪模块4、特征量提取模块5和健康状态分析模块6,多通道数据采集卡2,告警装置/控制断路器3。需要对控制***的哪一个具体部分进行健康状态监测,就将特征提取与健康实时分析装置并入控制***的该部分旁。
下面以算例仿真和实际案例仿真说明实施方案的有效性。
1、算例仿真
构造一个调频调幅并加入噪声的非线性信号予以仿真。信号如下:
该信号由调幅频率为25Hz,基频为10Hz,调频频率为20Hz,和频率为100Hz的正弦信号以及加入的高斯随机噪声r(t)、脉冲噪声z(t)混合而成,采样频率为1kHz,采样点为1024。
对信号f(t)进行中值滤波,可以发现中值滤波对脉冲噪声有很明显的抑制作用,但对随机噪声的滤波效果不理想,因此需要进一步使用小波变换滤除信号中的随机噪声。
利用Matlab进行仿真,小波函数选择db5小波,分解层数确定为5层,对比硬阈值、软阈值和本文改进阈值处理方法,对信号进行小波阈值降噪。
为了更加直观地对比各降噪处理方法的效果,引入信噪比SNR和均方根误差MSE。
SNR定义为
MSE定义为
式中,为中值滤波后的含噪信号,为小波阈值降噪后的重构信号,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分量能量和能量熵
究其原因:当执行器处于正常状态时,能量分布比较均匀,所以熵值较大;而当执行器处于故障状态时,能量主要分布在故障频段,故而熵值较小。因此,通过能量熵的方法可以较为直观的判断出执行器是否处于故障状态,以及处于何种故障状态。而且,控制***中执行器、传感器以及被控对象发生故障时的能量熵都会处在不同的区段上,因此该方法还可以推广到对整个控制***不同组成部分的状态检测与辨识上。
当取正常时能量熵为初始值H0,Thre=0.1即在正常的±10%的范围,则不难发现,轻度故障和重度故障时均会告警。
Claims (1)
1.基于联合降噪和经验模态分解的控制***健康状态分析方法,其特征在于:通过引入结合中值滤波与采用尺度间相关性的改进小波阈值联合降噪的方法,以及采用端点延拓、集合经验模态分解(EEMD)和相关系数阈值比较法相结合的方法,经过中值滤波、小波阈值降噪、经验模态分解、计算能量熵四个步骤将控制***状态特征提取出来,与初始检测的正常状态熵值比较,实时判断控制***运行状态,包括如下具体步骤:
步骤1)将离散含噪信号f(k)先进行中值滤波处理得到滤除可能的脉冲噪声,其中,f(k)信号长度为N;
步骤2)优化分解层数,包括如下步骤:
步骤2.1)对中值滤波处理后的进行第j层小波分解,其中,从j=1开始分解;
步骤2.2)对小波分解得到的低频近似系数aj(k)予以保留,对高频细节系数dj(k)进行式(1)自相关系数λ计算,若满足即各高频细节系数dj的自相关系数满足自由度l的χ2分布,则对dj继续进行j+1层分解;
步骤2.3)直到第Δ+1层的细节系数dΔ+1不能满足为止,确定优化分解层数为Δ层,同时得到各小波系数Wj(k);dj(k)为高频细节系数,为dj(k)的平均值;
<mrow>
<msub>
<mi>&lambda;</mi>
<mi>i</mi>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mi>k</mi>
<mrow>
<mi>N</mi>
<mo>-</mo>
<mi>i</mi>
</mrow>
</munderover>
<mo>&lsqb;</mo>
<msub>
<mi>d</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mover>
<mi>d</mi>
<mo>&OverBar;</mo>
</mover>
<mi>j</mi>
</msub>
<mo>&rsqb;</mo>
<mo>&lsqb;</mo>
<msub>
<mi>d</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mi>i</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mover>
<mi>d</mi>
<mo>&OverBar;</mo>
</mover>
<mi>j</mi>
</msub>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mi>k</mi>
<mrow>
<mi>m</mi>
<mo>-</mo>
<mi>i</mi>
</mrow>
</munderover>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msub>
<mi>d</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mover>
<mi>d</mi>
<mo>&OverBar;</mo>
</mover>
<mi>j</mi>
</msub>
<mo>&rsqb;</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>,</mo>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
<mo>,</mo>
<mi>l</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
步骤3)通过使用式(2)的改进阈值以满足:①阈值随着分解尺度的递增而逐渐减小,经过小波分解后不同分解层的系数比例分布不同;②避免出现阈值较小而起不到去除尽可能多噪声的作用;σ为噪声标准差,N为信号采样长度,j为分解尺度;
步骤4)通过使用式(3)的改进阈值函数以克服:①硬阈值函数的不连续性和引起信号的附加振荡;②软阈值函数存在恒定偏差,影响重构信号与真实信号的逼近程度;W为含噪信号小波变换后的小波系数,δ为阈值,Wδ为经过阈值降噪后的小波系数,μ,v,p,q均为可调参数,其中μ>0,v>1,p∈[0,1],q≥0;
<mrow>
<mi>&delta;</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mi>&sigma;</mi>
<msqrt>
<mrow>
<mn>2</mn>
<mi>ln</mi>
<mrow>
<mo>(</mo>
<mi>N</mi>
<mo>)</mo>
</mrow>
</mrow>
</msqrt>
</mrow>
<msqrt>
<mi>j</mi>
</msqrt>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msub>
<mi>W</mi>
<mi>&delta;</mi>
</msub>
<mo>=</mo>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>sgn</mi>
<mrow>
<mo>(</mo>
<mi>W</mi>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<mo>|</mo>
<mi>W</mi>
<mo>|</mo>
<mo>-</mo>
<mfrac>
<mrow>
<mn>2</mn>
<mi>p</mi>
</mrow>
<mrow>
<mn>1</mn>
<mo>+</mo>
<msup>
<mi>e</mi>
<mrow>
<mi>q</mi>
<mrow>
<mo>(</mo>
<mo>|</mo>
<mi>W</mi>
<mo>|</mo>
<mo>-</mo>
<mi>&delta;</mi>
<mo>)</mo>
</mrow>
</mrow>
</msup>
</mrow>
</mfrac>
<mi>&delta;</mi>
<mo>)</mo>
</mrow>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mrow>
<mo>|</mo>
<mi>W</mi>
<mo>|</mo>
</mrow>
<mo>&GreaterEqual;</mo>
<mi>&delta;</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>W</mi>
<mo>&CenterDot;</mo>
<msup>
<mi>v</mi>
<mrow>
<mo>-</mo>
<msup>
<mrow>
<mo>(</mo>
<mi>&mu;</mi>
<mo>|</mo>
<mi>W</mi>
<mo>|</mo>
<mo>-</mo>
<mi>&delta;</mi>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</msup>
<mo>,</mo>
</mrow>
</mtd>
<mtd>
<mrow>
<mrow>
<mo>|</mo>
<mi>W</mi>
<mo>|</mo>
</mrow>
<mo>&le;</mo>
<mi>&delta;</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<mo>)</mo>
</mrow>
</mrow>
步骤5)对处于阈值δ的[δ(1-α),δ(1+α)]邻域内的小波系数通过式(4)计算尺度间相关量θ(k),其中α为调节因子且调节α可以改变估计算法的精度;将θ(k)∈[0,β]即具有较强相关性的小波系数记为其余记为Wj(k)为含噪信号f(k)在j层分解时的小波系数,max Wj(k)为|Wj(k)|中的最大值,min Wj(k)为|Wj(k)|中的最小值,其中j=1,2,...,Δ;
<mrow>
<mi>&theta;</mi>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mi>ln</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<mi>max</mi>
<mi> </mi>
<msub>
<mi>W</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>min</mi>
<mi> </mi>
<msub>
<mi>W</mi>
<mi>j</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>4</mn>
<mo>)</mo>
</mrow>
</mrow>
步骤6)对不处于阈值δ的[δ(1-α),δ(1+α)]邻域内的其余所有小波系数Wj(k)进行判断,若|Wj(k)|≥δ,则将Wj(k)记为否则记为
步骤7)对以|W|≥δ的情形带入式(3)加以处理,对以|W|<δ的情形带入式(3)加以处理;
步骤8)对处理后的各小波系数进行小波逆变换,得到降噪后的真实信号的估计
步骤9)进行端点左右延拓,以右延拓为例进行介绍,包括如下步骤:
步骤9.1)确定待延拓信号段:记的所有极大值点为{M0,M1,...,Mk},所有极小值点为{m0,m1,...,mk},其中极大值点和极小值点分别从右向左标记,并记信号右端点为S;不妨假设从S向左先有极大值点后有极小值点,则记S到m0的波形为待延拓信号段ω0,待延拓长度为L,其中S到M0的长度记为L′;
步骤9.2)将ω0以M0为参考点在所有极大值点集合中依次移动,并使M0始终与Mi重合,通过式(5)计算长度为L的波形ω0和ωi的匹配度系数ξi;
步骤9.3)取ξi最小时的Mi为Mp,波形ωi记为ωp作为最佳匹配波形,则此时S所对应的坐标为Xp=Mp+L′;
步骤9.4)从Xp的后一点开始,将实际波形依次延拓到S后,一直延拓15个采样点,为S到m0的波形为待延拓信号段ω0与ωi的协方差,D[·]为方差;
<mrow>
<msub>
<mi>&xi;</mi>
<mi>i</mi>
</msub>
<mo>=</mo>
<mfrac>
<msub>
<mi>d</mi>
<mi>i</mi>
</msub>
<msub>
<mi>&kappa;</mi>
<mi>i</mi>
</msub>
</mfrac>
<mo>=</mo>
<mfrac>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>L</mi>
</munderover>
<mo>|</mo>
<mover>
<mi>s</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<msub>
<mi>M</mi>
<mi>i</mi>
</msub>
<mo>-</mo>
<mi>L</mi>
<mo>+</mo>
<msup>
<mi>L</mi>
<mo>&prime;</mo>
</msup>
<mo>+</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mover>
<mi>s</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<msub>
<mi>M</mi>
<mn>0</mn>
</msub>
<mo>-</mo>
<mi>L</mi>
<mo>+</mo>
<msup>
<mi>L</mi>
<mo>&prime;</mo>
</msup>
<mo>+</mo>
<mi>j</mi>
<mo>)</mo>
</mrow>
<mo>|</mo>
</mrow>
<mrow>
<mfrac>
<mrow>
<mi>cov</mi>
<mrow>
<mo>(</mo>
<msub>
<mover>
<mi>s</mi>
<mo>^</mo>
</mover>
<msub>
<mi>M</mi>
<mi>i</mi>
</msub>
</msub>
<mo>,</mo>
<msub>
<mover>
<mi>s</mi>
<mo>^</mo>
</mover>
<msub>
<mi>M</mi>
<mn>0</mn>
</msub>
</msub>
<mo>)</mo>
</mrow>
</mrow>
<msqrt>
<mrow>
<mi>D</mi>
<mo>&lsqb;</mo>
<msub>
<mover>
<mi>s</mi>
<mo>^</mo>
</mover>
<msub>
<mi>M</mi>
<mi>i</mi>
</msub>
</msub>
<mo>&rsqb;</mo>
<mo>&CenterDot;</mo>
<msqrt>
<mrow>
<mi>D</mi>
<mo>&lsqb;</mo>
<msub>
<mover>
<mi>s</mi>
<mo>^</mo>
</mover>
<msub>
<mi>M</mi>
<mn>0</mn>
</msub>
</msub>
<mo>&rsqb;</mo>
</mrow>
</msqrt>
</mrow>
</msqrt>
</mfrac>
<mo>+</mo>
<mn>2</mn>
</mrow>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>5</mn>
<mo>)</mo>
</mrow>
</mrow>
步骤10)对端点延拓后的信号进行EEMD分解,通过比较各本征模态函数(IMF)与信号的相关系数ρ与预设阈值ζ的大小,剔除“伪分量”,其中,ρ通过式(6)计算,σ[·]为标准差,ζ=0.1max(ρ);
<mrow>
<mi>&rho;</mi>
<mo>=</mo>
<mfrac>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mover>
<mi>s</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>&CenterDot;</mo>
<msub>
<mi>imf</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>&sigma;</mi>
<mo>&lsqb;</mo>
<mover>
<mi>s</mi>
<mo>^</mo>
</mover>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
<mo>&CenterDot;</mo>
<mi>&sigma;</mi>
<mo>&lsqb;</mo>
<msub>
<mi>imf</mi>
<mi>i</mi>
</msub>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
</mfrac>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>6</mn>
<mo>)</mo>
</mrow>
</mrow>
步骤11)对剔除“伪分量”后的η个IMF分量,计算各IMF分量ci(k)的能量其中i=1,2,...,η;
步骤12)通过式(7)计算能量熵;
<mrow>
<mi>H</mi>
<mo>=</mo>
<mo>-</mo>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>&eta;</mi>
</munderover>
<mfrac>
<msub>
<mi>E</mi>
<mi>i</mi>
</msub>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>&eta;</mi>
</munderover>
<msub>
<mi>E</mi>
<mi>i</mi>
</msub>
</mrow>
</mfrac>
<mo>&CenterDot;</mo>
<mi>ln</mi>
<mrow>
<mo>(</mo>
<mfrac>
<msub>
<mi>E</mi>
<mi>i</mi>
</msub>
<mrow>
<munderover>
<mi>&Sigma;</mi>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>&eta;</mi>
</munderover>
<msub>
<mi>E</mi>
<mi>i</mi>
</msub>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>7</mn>
<mo>)</mo>
</mrow>
</mrow>
2
步骤13)通过得到的能量熵H判断检测阶段的健康状态,具体为:在检测阶段中,若上位机则告警,提示控制***出现故障;其中,H0为初始控制***健康状态检测得到的能量熵值,Thre为正常状态判断阈值,Thre取值为0.1~1。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510172003.8A CN105094111B (zh) | 2015-04-09 | 2015-04-09 | 基于联合降噪和经验模态分解的控制***健康状态分析方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510172003.8A CN105094111B (zh) | 2015-04-09 | 2015-04-09 | 基于联合降噪和经验模态分解的控制***健康状态分析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105094111A CN105094111A (zh) | 2015-11-25 |
CN105094111B true CN105094111B (zh) | 2017-12-05 |
Family
ID=54574789
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510172003.8A Expired - Fee Related CN105094111B (zh) | 2015-04-09 | 2015-04-09 | 基于联合降噪和经验模态分解的控制***健康状态分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105094111B (zh) |
Families Citing this family (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP6444851B2 (ja) * | 2015-12-24 | 2018-12-26 | ファナック株式会社 | ノイズの発生原因を検出する学習機能を有する制御装置 |
CN105699883B (zh) * | 2016-01-22 | 2018-08-14 | 合肥工业大学 | 一种模拟电路健康预测方法 |
CN105956514B (zh) * | 2016-04-15 | 2019-07-23 | 南京航空航天大学 | 一种振动数据驱动的直升机旋翼异常检测方法 |
CN106798554A (zh) * | 2017-01-12 | 2017-06-06 | 安徽大学 | 一种含噪imf分量及心电信号的去噪方法 |
CN106652456A (zh) * | 2017-01-20 | 2017-05-10 | 上海琛岫自控科技有限公司 | 用于路段拥堵报警的无人机检测*** |
CN110070874B (zh) * | 2018-01-23 | 2021-07-30 | 中国科学院声学研究所 | 一种针对声纹识别的语音降噪方法及装置 |
CN109580787B (zh) * | 2018-12-08 | 2021-02-02 | 国网四川省电力公司广安供电公司 | 用于变压器高压套管引线超声检测的超声回波去噪方法 |
CN112697887B (zh) * | 2020-12-08 | 2023-05-02 | 江苏科技大学 | 基于神经网络的超声检测缺陷定性识别方法 |
CN115811456B (zh) * | 2022-11-16 | 2024-02-23 | 南京源兴智达信息科技有限公司 | 一种基于ofdm的车载通信*** |
CN116703183B (zh) * | 2023-08-04 | 2023-09-29 | 国网吉林省电力有限公司经济技术研究院 | 一种新型配电网***用电负荷评估方法 |
CN116828354B (zh) * | 2023-08-30 | 2023-11-07 | 深圳市智纬科技有限公司 | 一种用于领夹无线麦克风的收音质量优化方法及*** |
CN117878971B (zh) * | 2024-03-12 | 2024-06-11 | 西安热工研究院有限公司 | 一种熔盐储能耦合火电机组的调频方法和*** |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7485095B2 (en) * | 2000-05-30 | 2009-02-03 | Vladimir Shusterman | Measurement and analysis of trends in physiological and/or health data |
KR100954817B1 (ko) * | 2009-03-13 | 2010-04-28 | (주)바이오센스크리에티브 | 맥파신호분석을 이용한 혈관건강 및 스트레스 검사시스템 및 방법 |
CN102436586A (zh) * | 2011-10-28 | 2012-05-02 | 哈尔滨工业大学 | 一种基于小波阈值降噪和经验模态分解的高光谱图像分类方法 |
CN102854874A (zh) * | 2012-06-18 | 2013-01-02 | 南京航空航天大学 | 一种基于联合多观测器的故障诊断与容错控制装置及方法 |
CN103674511A (zh) * | 2013-03-18 | 2014-03-26 | 北京航空航天大学 | 一种基于emd-svd与mts的机械磨损件性能评估与预测方法 |
-
2015
- 2015-04-09 CN CN201510172003.8A patent/CN105094111B/zh not_active Expired - Fee Related
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7485095B2 (en) * | 2000-05-30 | 2009-02-03 | Vladimir Shusterman | Measurement and analysis of trends in physiological and/or health data |
KR100954817B1 (ko) * | 2009-03-13 | 2010-04-28 | (주)바이오센스크리에티브 | 맥파신호분석을 이용한 혈관건강 및 스트레스 검사시스템 및 방법 |
CN102436586A (zh) * | 2011-10-28 | 2012-05-02 | 哈尔滨工业大学 | 一种基于小波阈值降噪和经验模态分解的高光谱图像分类方法 |
CN102854874A (zh) * | 2012-06-18 | 2013-01-02 | 南京航空航天大学 | 一种基于联合多观测器的故障诊断与容错控制装置及方法 |
CN103674511A (zh) * | 2013-03-18 | 2014-03-26 | 北京航空航天大学 | 一种基于emd-svd与mts的机械磨损件性能评估与预测方法 |
Non-Patent Citations (2)
Title |
---|
基于EEMD和自相关函数特性的自适应降噪方法;余发军 等;《计算机应用研究》;20150131;第32卷(第1期);第206-209页 * |
基于改进EEMD方法的数字滤波器;焦彦军 等;《电力自动化设备》;20111130;第31卷(第11期);第64-67页 * |
Also Published As
Publication number | Publication date |
---|---|
CN105094111A (zh) | 2015-11-25 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105094111B (zh) | 基于联合降噪和经验模态分解的控制***健康状态分析方法 | |
CN110135492B (zh) | 基于多高斯模型的设备故障诊断与异常检测方法及*** | |
CN105469138B (zh) | 基于粒子群和支持向量机的控制***执行器故障诊断方法 | |
CN104236911B (zh) | 一种列车转向架轴承服役过程监测与故障诊断***及方法 | |
CN103674511A (zh) | 一种基于emd-svd与mts的机械磨损件性能评估与预测方法 | |
Dambros et al. | Oscillation detection in process industries–Part I: Review of the detection methods | |
CN103033567B (zh) | 一种基于导波的管道缺陷信号识别方法 | |
Li et al. | Research on test bench bearing fault diagnosis of improved EEMD based on improved adaptive resonance technology | |
CN106796579A (zh) | 用于自动检测旋转轴上的缺陷的方法和*** | |
CN101936767B (zh) | 提取工程机械运行特征信号的方法 | |
CN109946075A (zh) | 一种轴承状态监测及故障诊断方法 | |
CN105092239A (zh) | 一种齿轮早期故障检测方法 | |
CN106771598B (zh) | 一种自适应谱峭度信号处理方法 | |
Yan et al. | Early chatter detection in thin-walled workpiece milling process based on multi-synchrosqueezing transform and feature selection | |
CN113177537B (zh) | 一种旋转机械设备的故障诊断方法及*** | |
CN105424366A (zh) | 基于eemd自适应消噪的轴承故障诊断方法 | |
CN102901630A (zh) | 基于自适应冗余提升小波降噪分析的轴承故障识别方法 | |
CN109190598A (zh) | 一种基于ses-lof的旋转机械监测数据噪点检测方法 | |
Han et al. | ESPRIT-and HMM-based real-time monitoring and suppression of machining chatter in smart CNC milling system | |
CN103175897A (zh) | 一种基于振动信号端点检测的高速道岔伤损识别方法 | |
Sun et al. | Application of wavelet soft threshold denoising algorithm based on EMD decomposition in vibration signals | |
CN107563403B (zh) | 一种高速列车运行的工况识别方法 | |
Patel et al. | Fault diagnostics of rolling bearing based on improve time and frequency domain features using artificial neural networks | |
CN106441843A (zh) | 一种旋转机械故障波形识别方法 | |
Wan et al. | Chatter detection for micro milling considering environment noises without the requirement of dominant frequency |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20171205 |
|
CF01 | Termination of patent right due to non-payment of annual fee |