基于NGWarblet-WVD 的高质量时频分析方法

郝国成 , 冯思权 王巍 凌斯奇 谭淞元

时频分析方法采用时间和频率域联合函数来处理时序信号,获取信号频率随时间变化的细节信息,是现代信号处理领域的重要技术手段之一.1946 年,Gabor 首次对傅里叶变换加高斯窗函数,提出了著名的Gabor 变换,由此开启了时频联合分析的新思路[1].Ville 将量子力学的Wigner 分布用于信号分析与处理领域,提出了Wigner-Ville 分布(Wigner-Ville distribution,WVD)[2].WVD 具有良好的时频聚集特性,但存在严重的交叉项干扰问题.Cohen对WVD 进行时频二维卷积得到Cohen 类时频分布[3].Cohen 类时频分布通过构造核函数,达到消除或抑制交叉项的目的,缺点是降低了频率聚集度.20 世纪80 年代,Mallat[4]提出了多尺度分析思想和Mallat 算法,成功地统一了各种小波函数的构造模型,使用可调节的时间和频率窗口,有效提高了频率聚集度.1996 年,美国地球物理学家Stockwell等[5]对短时傅里叶变换和连续小波变换的思想进行延伸与推广,提出了S 变换.1998 年,Huang等[6]提出经验模态分解(Empirical mode decomposition,EMD),将信号分解为有限个固有模态函数(Intrinsic mode functions,IMF)的集合.EMD 方法可应用于多种类型信号的分解,在处理非平稳、非线性信号方面,效果良好[7−8].但相对来讲,S 变换和EMD的聚集度不够理想.2012 年,Yang等[9]提出了适用范围更广泛的参数化时频分析广义Warblet变换(Generalized Warblet transform,GWT),该方法具有真实反映信号频率分布的特点,但存在频率泄露现象,时频聚集度较差,需要结合其他算法加以改进.

WVD 作为一种优良的双线性时频分析方法具有其他方法不可替代的高锐化时频聚集特性,但在处理多分量信号时会出现交叉项干扰问题.文献[10−12]提出对信号进行WVD 处理前,分别采用变分模态分解 (Variational mode decomposition,VMD)、集成经验模态分解 (Ensemble empirical mode decomposition,EEMD)、自适应匹配追踪 (Adaptive matching pursuit,AMP)将多分量信号转换为单分量信号,有效避免产生交叉项;文献[13]通过对信号进行带通滤波和相位矫正的方法去除交叉项;文献[14]提出矩阵旋转变换的方法,将WVD 交叉项旋转至与频率轴平行,再通过滤波器滤除交叉项;文献[15]通过二重余弦信号的WVD,推导自项与交叉项位置关系以及振荡特性,滤除交叉项;文献[16−17]采用两种算法结合的思想(如Gabor-WVD、BGabor-NSPWVD、SPWVD-WVD 等),对比实验证明此方法有效抑制了WVD 交叉项.

本文在分析WVD 交叉项的产生原因基础上,将交叉项分为新产生的交叉项分量和混入自项的交叉项分量两种类型.利用GWT 较好还原信号真实频率分布的特性,将GWT 矩阵与WVD 矩阵联合处理以实现滤波效应,抑制WVD 的两种类型交叉项.使用两种定量评价方法将NGWT-WVD 算法与同类算法进行对比,检验算法有效性.最后将该算法用于处理金属破裂样本信号,获取破裂期间的时频分布图,找出滤波器的窗口门限近似频率,为声发射信号监测传感器采集卡提供门限设置依据.

1 NGWT-WVD 算法

1.1 WVD (Wigner-Ville 分布)

WVD 是一种双线性时频分布,设信号为z(t),其Wigner-Ville 分布表达式为

R(t,τ)=z(t+τ/2)z∗(t −τ/2)表示信号的自相关函数,WVD 可看作信号自相关函数的傅里叶变换.WVD 具有优良的时频聚集度,但处理多分量信号时会出现交叉项,如式(2)所示

Ra(t,τ)表 示自相关成分(自项),是有用信息;Rs(t,τ)表示互相关成分(交叉项),是干扰信息.对于同时间区间的多分量信号,交叉项具有两个特点: 任意两个信号分量均会产生一个交叉项;交叉项位于两频率分量中间.

选取三分量线性调频信号z1(t) 为例,说明信号各成分之间的关系.其表达式如式(3)所示

信号二维、三维WVD 时频分布如图1 所示.由图1(a)和图1(b)对比可见,三分量线性调频信号的Wigner-Ville分布引入了一些新分量,这些分量是新产生的交叉项分量;而由图1(c)的三维图可见,中间分量能量远高于两边分量,说明有部分交叉项混入进自项成分.本文需要解决的问题有两个: 一是抑制新产生的交叉项分量;二是抑制混入自项的交叉项分量.

图1 三分量信号的WVD 时频图Fig.1 Time-frequency diagram of WVD of three-components signal

1.2 GWT (广义Warblet 变换)

GWT 是核函数以傅里叶级数为模型定义的参数化时频分析方法.其定义为

其中,t0∈R 表示时间窗滑动时的窗中心所在时间,wσ∈L2(R)定义了一个非负对称的标准化实窗,通常是高斯窗,(t) 计算式为

参数化时频分析方法选取合适的核函数对其分析效果有很大影响,采用傅里叶级数变换核的GWT能够分析具有周期性或非周期性时频特征的非平稳信号,以及具有强振荡时频特征的信号,使其适用范围更加广阔.

三分量线性调频信号的GWT 时频图如图2 所示.

图2 中GWT 算法采用傅里叶级数模型作为核函数,较好地还原了信号的真实频率分布,虽然时频聚集度较差,但能够保留真实频率分量,不会产生交叉项.由此,基于WVD 的高锐化特性和GWT真实还原信号时频分布的特性,将二者相结合,得到GWT-WVD 算法,既能抑制交叉项,又能保留高锐化时频聚集度的性能.

图2 三分量信号的GWT 时频图Fig.2 Time-frequency diagram of GWT of three-components signal

1.3 GWT-WVD 算法: 抑制新产生的交叉项分量

本文提出的GWT-WVD 算法,可有效抑制WVD 的交叉项分量.该算法通过对GWT 算法与WVD 算法得到的矩阵进行运算,在保留较好的时频聚集度的同时,能较好地消除或抑制交叉项.其表达式为

其中,GWarx(t,f)与W x(t,f) 分别代表GWT与WVD,p(x,y) 为联合处理函数.本文采用了3种不同的函数,得到GWT-WVD 算法的3种定义式

式(7)～(9)分别采用最小值法、二值化法、幂指数调节法对两矩阵进行处理,其各自实现思想如下.

1)最小值法.比较WVD 矩阵及GWT 矩阵对应位置元素,筛选其中较小的元素,按比例处理后组成GWT-WVD 矩阵,使WVD 矩阵新产生的交叉项所在位置的元素被GWT 矩阵元素取代.

2)幂指数调节法.调节幂指数,增强GWT与WVD矩阵中信号数据对应的元素,削弱交叉项数据对应的元素,再将两矩阵点乘,得到GWT-WVD矩阵.

3)二值化法.选取合适的阈值将GWT 矩阵二值化得到新矩阵,用新矩阵点乘WVD 矩阵得到GWTWVD 矩阵.WVD 矩阵有交叉项的元素位置对应的GWT 矩阵元素会小于阈值,二值化后为0,与WVD 矩阵相乘后可消除新产生的交叉项.

线性调频信号z1(t) 的3种GWT-WVD 时频分布如图3 所示.3种结合方法均能优化信号的时频分布,有效抑制交叉项,但幂指数调节法的时频聚集度比另外两种方法差(第3 节讨论时频聚集度问题).在处理复杂信号时,二值化法阈值难以确定.综合比较,最小值法的GWT-WVD 算法性能最佳.

图3 3种GWT-WVD 时频图Fig.3 Time-frequency diagram of three types of GWT-WVD

1.4 NGWT-WVD 算法: 抑制混入自项的交叉项分量

GWT-WVD 算法能够有效抑制新产生的交叉项分量,但无法抑制混入自项分量的交叉项.针对这个问题,本文对GWT-WVD 算法进行改进,又提出NGWT-WVD 算法,其主要实现思路如下.

1)采用广义Warblet 变换和WVD 分别对原信号进行处理,得到GWT 矩阵和WVD 矩阵;

2) 找出G W T 矩阵中元素数值的最大值GWTmax,并记录其所对应的位置 (i,j),将GWT矩阵中的各元素除以GWTmax,即对GWT 矩阵进行归一化,得到矩阵GWT-1;

3) 记录GWT-1 矩阵中元素数值的最小值GWTmin,最小值GWTmin 要求非零,并用GWTmin的数值替换掉矩阵GWT-1 中所有值为0 的元素;

4)找出WVD 矩阵中位置为 (i,j) 的元素,将其记为WVDmax,同时将WVD 矩阵中的各个元素除以WVDmax,得到矩阵WVD-1;

5)用矩阵WVD-1 点除矩阵GWT-1,得到矩阵T,选取矩阵T中大于x的元素以及小于y的元素,x设置为5,y设置为2,将所对应的元素位置置1,并记录大于x的元素位置;

6)在WVD 矩阵中找出与上一步记录对应的元素位置,并将此位置上元素置为0,最后用WVD矩阵点除矩阵T,输出NGWT-WVD 矩阵.其流程如图4 所示.

图4 NGWT-WVD 算法流程图Fig.4 Algorithm flowchart of NGWT-WVD

线性调频信号的GWT-WVD 算法和NGWTWVD 算法三维时频分布如图5 所示.由图5(a)可见,中间分量的能量谱明显高于其他分量的能量谱,GWT-WVD 算法不能有效抑制混入自项的交叉项分量;而图5(b)中3 个信号分量能量一致,混入自项成分上的交叉项被有效抑制.图5 说明NGWTWVD 算法能有效抑制线性调频信号的两类交叉项.

图5 三分量信号三维时频图比较Fig.5 Three-dimensional time-frequency diagrams comparison of three-component signals

NGWT-WVD 算法中的阈值参数x是为了抑制存在于自项中的干扰项,阈值y是为了抑制新产生的干扰项并且去除发散的能量.阈值对NGWTWVD 算法性能起决定性作用,这里使用后文的两种定量评价方式对两个阈值的敏感性进行分析,其结果如图6 所示.

图6(a)是CM值(时频聚集度定量评价方法,数值越大代表聚集度越高,其计算式见第3.2 节)随阈值y变化折线图,此时阈值x设置为5,x1,x2,x3,x4 分别代表文中使用的4种仿真信号.图中,阈值y小于1 时,NGWT-WVD 算法的时频聚集度低于WVD 算法;随着阈值y的增加,CM值逐渐增加,WVD 中部分扩散的时频系数被滤除,时频聚集度提高.当阈值y超过2 时,时频聚集度CM值迅速提高,但这是以滤除部分WVD 有用信息分量为代价.因此阈值y取值为2 即可,此时新产生的交叉项基本消除干净(由图5(b)可见),且时频聚集度较理想.图6(b)是时变功率谱误差随阈值x变化折线图,此时阈值y设置为2 (基本消除了新产生的交叉项分量).阈值x是为滤除混入自项的交叉项(只有仿真信号1 含有混入自项的交叉项),当阈值x大于10 时,时变功率谱误差保持在25%左右,无法滤除混入自项的交叉项;阈值x在5～6 之间时,时变功率谱误差接近0,此时抑制效果较好;小于5 后,由于滤除了信息项,导致误差进一步增加,因此阈值x设置为5 较合适.实际应用信号由一系列单分量信号线性叠加而成,也可按照此参数设置相应阈值.

图6 阈值敏感性测试图Fig.6 The test chart of threshold sensitivity

NGWT-WVD 算法通过对GWT 矩阵进行归一化且做去零处理得到的GWT-1 矩阵,此矩阵为对照矩阵,通过设置两个阈值,实现了滤波效应,剔除了WVD 中发散能量,且抑制了混入自项的交叉项和新产生的交叉项,得到了更加理想的时频分布结果.

2 数值实验

为了评价NGWT-WVD 算法的处理效果,本文另外构造了3种具有代表性的信号,并通过与Garbor-WVD、VMD-WVD 两种较有代表性的方法进行对比,验证该算法的有效性.

2.1 分段信号

函数表达式为

分段信号z2(t) 的WVD、Gabor-WVD (文献[14]中采用的处理方法)、VMD-WVD (选取模数为4)、NGWT-WVD 算法时频分布如图7 所示.图7(a)中WVD 在端点处出现严重的交叉项干扰;图7(b)中Gabor-WVD 算法处理信号时仍然会出现少许交叉项;图7(c) 中,VMD 的参数设置准确的情况下,VMD-WVD 算法交叉项抑制效果较好;图7(d)中NGWT-WVD 算法得到的时频分布也基本没有交叉项.

图7 分段信号时频图Fig.7 Time-frequency diagram of segmented signal

2.2 交叉型信号

函数表达式为

交叉信号z3(t) 的WVD、Gabor-WVD、VMDWVD (选取模数为2)以及NGWT-WVD 算法的时频分布如图8 所示.此信号为两个线性调频信号交叉于一点,交叉点的频率分量混合在一起,分离困难.图8(b)中Gabor-WVD 算法在处理该信号时,不能十分有效地抑制交叉项,且时频聚集度较低.图8(c)中VMD-WVD 算法交叉项抑制效果较好,时频聚集度略差;图8(d)中NGWT-WVD 算法对此信号的频率分量刻画准确,时频聚集度高.

图8 交叉型信号时频图Fig.8 Time-frequency diagram of cross-type signal

2.3 两分量调频信号

函数表达式为

信号z4(t) 是由一个抛物线调频信号和一个正弦调频信号组成,并存在交叉点.该信号的WVD、Gabor-WVD、VMD-WVD (选取模数为2)以及NGWT-WVD 算法的时频分布如图9 所示.由图可见,图9(b)中Garbor-WVD 不能完全抑制,图9(c)中VMD-WVD 在处理这一类复杂调频信号时,即使参数设置准确,也无法有效将单分量信号完全分开,交叉项干扰较为严重;而NGWT-WVD 算法对交叉项具有较好的抑制效果.由此,以上3 个仿真实验都能说明NGWT-WVD 算法具有较好的抑制信号交叉项效果,同时保留了高时频聚集度,提高了时频分析质量.

图9 两分量调频信号时频图Fig.9 Time-frequency diagram of two-component frequency modulated signal

3 算法性能评价

为定量说明NGWT−WVD 算法性能,本文选取交叉项抑制效果评价和时频聚集度评价两个指标评判该算法的时频分析效果.

3.1 交叉项抑制效果评价

交叉项抑制效果是WVD 改进算法性能评价的重要指标,以往均使用定性分析,即通过观察时频分布结果得出判断.本文在分析WVD 交叉项出现规律的基础上,提出一种定量的交叉项抑制效果评价方法.由式(2)知,信号经WVD 处理后,分量包括自项成分与交叉项成分,本文计算信号在时间−频率平面的时变功率谱Pi,并与信号标准时变功率谱IPi(不含交叉项的信号功率谱)作对比,计算出时变功率谱误差,以此评价交叉项抑制效果.选取信号z1(t),z2(t),z3(t),z4(t) 作为研究对象,计算各算法的功率谱与标准功率谱的平均相对误差,其计算式为

其中,N为信号采样点数.标准时变功率谱计算式为

其中,n表示信号含有单分量的个数,Pfj(i) 代表第j个单分量的时变功率谱.将各单分量分别求时频分布后相加,得到没有交叉项的标准时变功率谱.四种仿真信号平均时变功率谱误差对应的柱形图如图10 所示,数据如表1 所示,z1,z2,z3,z4 分别代表信号z1(t),z2(t),z3(t),z4(t).

图10 各算法的时变功率谱误差柱形图Fig.10 Time-varying power spectrum error column chart of each algorithm

表1 各算法的时变功率谱误差比较Table 1 Time-varying power spectrum error comparison of each algorithm

表1 中不同信号的NGWT-WVD 算法的时变功率谱误差均为最小,反映出算法抑制交叉项性能最佳.图11 为表1 的平均时变功率谱误差折线图形式.

图11 各算法的时变功率谱误差折线图Fig.11 Time-varying power spectrum error line chart of each algorithm

图11 中WVD 算法的误差值最大,说明时频分布中含有大量的交叉项,Gabor-WVD、GWTWVD 算法都能一定程度上抑制交叉项,但是无法有效抑制混入自项成分中的交叉项,其时变功率谱误差仍很大;VMD-WVD 算法会先分解信号,在处理恒频信号时,分解效果好,交叉项抑制效果也好,在处理复杂信号时,分解效果不佳,时变功率谱误差远高于NGWT-WVD 算法;而NGWT-WVD 算法可以较好地抑制两种交叉项,时变功率谱误差较小,且不需要设置初始参数,算法的适应性强.

3.2 时频聚集度评价

Shafi等[18]于2009 年提出了关于时频聚集度量化评定的方法,本文选取CM值作为时频聚集度量化标准,评价信号的时频图的聚集度,其计算式为

其中,n为时间窗长度,ω为信号在某点的频率,Q表示信号的时频分布.选取上述4种仿真信号作为实验对象,计算CM值.图12 为4种信号各时频分析方法的CM值柱形图,其数据如表2 所示.

图12 各算法的CM 值柱形图Fig.12 CM value column chart of each algorithm

表2 中NGWT-WVD 算法在评价不同信号的时频聚集度时,其CM值均为最大,反映出算法时频聚集度性能最佳,锐化程度最高.图13 时频聚集度CM值折线图.

表2 各算法的CM 值比较(×10−3)Table 2 CM value comparison of each algorithm(×10−3)

图13 中,GWT 算法CM值最小,时频聚集度最差.Gabor-WVD 算法CM值低于GWT-WVD算法,时频聚集度较GWT-WVD 算法低,其去除交叉项的同时,降低了时频聚集度.信号z1 因为是一种线性调频信号,其VMD-WVD、NGWT-WVD算法的CM值与WVD的CM值相差不大,但对于另外3种信号,VMD-WVD、NGWT-WVD 算法的时频聚集度CM值明显高于WVD 算法,且NGWT-WVD 算法略优于VMD-WVD 算法.通过比较4种信号各时频分析方法的时频聚集度CM值大小,验证了NGWT-WVD 算法具有高锐化的时频聚集度.

图13 各算法的CM 值折线图(×10−3)Fig.13 CM value line chart of each algorithm (×10−3)

4 金属破裂检测信号的时频分析应用

人造金刚石合成过程中,六面顶压机顶锤的破裂损坏是经常发生的生产事故.顶锤采用的钨钴类硬质合金,在高温高压环境中,长期处于超临界应力状态,会出现疲劳损伤,进而而导致顶锤破裂[19].若金刚石生产加工过程中未发现顶锤破裂前的异常情况,将会使六面顶压机顶锤出现不可逆性损坏,造成严重生产损失.六面顶压机和顶锤结构如图14所示.

图14 六面顶压机和硬质合金顶锤Fig.14 Cubic press and carbide anvil

声发射技术是一种有效的探伤检测手段.构件在外力或应变力的作用下会激发一定频谱的声发射信号,通过判断接收到的信号频谱强度来预判构件的缺陷严重程度[20].顶锤破裂大致分为: 裂纹成核、裂纹拓展、断裂3 个过程,在这3 个过程中应变能以弹性应力波的形式释放出来,会产生剧烈的声发射信号[20].

为了避免生产事故,需要在初步检测到钨钴合金出现破裂时立即终止加工过程.而采集到的声发射信号振幅微弱(如图15(a)所示),难于甄别判断.本文使用的疑似金属破裂信号数据采样频率为40 000 Hz,选取2 000 个采样点数据,由于篇幅限制,仅绘制时域图、WVD 算法时频分布图以及NGWT-WVD算法时频分布图,结果如图15 所示.两种时频分析方法的CM值对比如表3 所示,NGWT-WVD 的时频聚集度CM值最大,时频聚集度最好.

图15 疑似金属破裂样本时频分析Fig.15 Time-frequency analysis of suspected metal rupture samples

表3 六种算法的CM 值比较(×10−5)Table 3 CM value comparison of six algorithms(×10−5)

图15(a) 为某一疑似金属破裂样本信号时域图,该信号振幅微弱,不利于监测传感器报警阈值的设置;图15(b)和图15(c)分别为该样本信号WVD算法以及NGWT-WVD 算法处理得到的时频分布图.为便于更精确的分析,截取了频率发生剧烈波动的片段(采样点数为900～1 300 之间数据)进行分析,结果如图15(d)和图15(e)所示.图15(d)中WVD算法得到的时频分布图交叉项分量与信号分量混杂,整个时频分布图杂乱不清,难以确定出现金属破裂的精确时间节点,无法有效示警.图15(e)中可较为明显地看出抑制了图15(d) 中的交叉项(尤其是图15(d)方框中的主要交叉项),时频分布锐化聚集度有明显的提高,CM值验证了这一结果.图15(e)中椭圆框标记了金属破裂过程中频率较集中的频率分量区域,可作为监测滤波器组的通带上下限的选取范围.其中滤波器通带2 作为主要预警通带,滤波器通带1 由于频率较低,较容易受到外界噪声干扰,通带3 则由于频率阈值过高,容易遗漏低频预警信号,因此滤波器通带1、3 通常作为辅助预警通带,此外还可以在各通带间再设置部分滤波器通带,提高识别几率.实验结果表明,NGWTWVD 算法能够较精确地显示出各信号出现破裂的具体时间和频率窗口值,可为信号监测传感器和滤波器组提供可操作的判断阈值,提高设备的监测成功率.

5 结束语

本文分析了WVD 产生交叉项的原理,针对交叉项干扰和时频模糊问题,提出了NGWT-WVD算法.该算法不仅能够有效抑制新产生的交叉项分量,而且解决了Gabor-WVD 等算法无法消除混入自项成分的交叉项分量问题,在交叉项抑制效果评价和时频聚集度评价中表现良好.仿真结果表明,NGWT-WVD 算法能够实现保持高锐化聚集度的同时,有效抑制交叉项干扰,是一种高质量的时频分析方法.将该算法用于处理金属破裂样本信号,能够得到较为精确的信号时间和频率窗口值,为监测传感器报警阈值的设置和数据采集滤波器组的设计提供有效依据.