基于改进同步挤压小波变换识别信号瞬时频率*

2017-09-12刘景良郑锦仰郑文婷黄文金

振动、测试与诊断 2017年4期

关键词：理论值时频时变

刘景良，郑锦仰 ,郑文婷，黄文金

(1.福建农林大学交通与土木工程学院福州，350002)(2.福建工程学院土木工程学院福州，350118)

基于改进同步挤压小波变换识别信号瞬时频率*

刘景良1，郑锦仰1,郑文婷2，黄文金1

(1.福建农林大学交通与土木工程学院福州，350002)(2.福建工程学院土木工程学院福州，350118)

提出了改进同步挤压小波变换。首先，采用解析模态分解定理将非平稳响应信号分解为多个单分量信号；其次，通过频率转换将分解后的单分量信号从低频区域转换到高频区域以提高时间分辨率；然后，对频率重点关注区域的时频代表值进行局部挤压，提高时频曲线的频率精度；最后，通过两个数值算例和一个拉索试验验证了改进同步挤压小波变换方法的有效性和准确性。研究结果表明，改进同步挤压小波变换能够提高信号在特定区域内的瞬时频率识别精度，是一种可行的时变结构参数识别方法。

改进同步挤压小波变换；瞬时频率；时变结构；时频分析；非平稳信号

引言

实际土木工程结构在承受极限荷载或长期工作荷载时，本质上是时变和非线性结构系统，其响应信号呈现非平稳性。时频分析方法能够有效处理非平稳信号，因而在时变结构参数识别领域受到广泛关注，主要方法有短时傅里叶变换、二次型时频分布[1]、小波变换[2]、以经验模态分解为主的Hilbert-Huang变换(Hilbert-Huang transform，简称HHT)[3-4]、集合经验模态分解[5]和解析模态分解(analytical modal decomposition,简称AMD)[6]等，其中以HHT和小波变换方法研究最为深入。HHT本质上是一种经验的局域分析方法，它无法分离密集的模态响应，特别是模态频率叠混的信号[7]。Yan等[8]将HHT变换与连续小波变换进行对比，结果表明小波变换在识别密集模态时优于HHT。

小波变换能够自适应调整时窗跟频窗大小并实现多分辨率分析，十分适合非平稳信号的分析处理。孙鹏等[9]使用Morlet小波变换识别悬索桥模型的固有频率和阻尼比。许鑫等[10]基于状态空间和小波理论提出了时变系统的参数识别方法。王超等[11]提出一种基于动态规划提取信号小波脊线和瞬时频率(instantaneous frequency,简称IF)的方法，该方法具有较强的抗噪性和较高的计算效率。尽管小波变换在时变结构参数识别领域中取得成功，但如何提取清晰的时频曲线目前仍没有得到很好的解决。Daubechies等[12]提出同步挤压小波变换(synchrosqueezing wavelet transform,简称SST)，通过重组小波变换后的时频图获得较高频率精度的时频曲线。刘景良等[13]将同步挤压小波变换引入土木工程领域，识别了时变结构响应信号的瞬时频率。但是，同步挤压小波变换方法只能处理信号频率不变时尺度方向的扩散，对于时间维度上的扩散无能为力。Li等[14]针对标准同步挤压算法存在小波系数扩散的问题，提出了广义的同步挤压小波变换算法。该方法能够获得较高精度的时频曲线，但前提条件是响应信号的数学解析表达式必须已知，而实际工程中的响应信号并不满足这一先决条件。汪祥莉等[15]针对混沌干扰背景下多个谐波信号的提取问题，提出一种基于同步挤压小波变换的谐波信号抽取方法，但该方法无法对非平稳谐波信号进行有效处理。可见，同步挤压小波变换作为一种以小波变换为基础的全新的时频分析方法，目前虽有相关的理论研究，但其算法的改进和在实际工程中的应用还十分缺乏。因此，如何对同步挤压小波变换进行改进并使其能够同时处理时间和频率两个维度方向的扩散是一个急需解决的关键问题。

笔者提出的改进同步挤压小波变换方法首先通过解析模态分解定理将多分量响应信号分解为单分量信号；其次，通过频率转换将原始信号的瞬时频率值从低频区域转换到高频区域，达到了提高时间分辨率的目的；然后，对重点关注区域的时频代表值进行重组，提高了该区域瞬时频率识别值的频率精度。改进同步挤压小波变换摒弃了标准同步挤压小波变换沿整个频率轴对连续小波变换后的时频代表值进行挤压的做法，对重点区域的时频代表值进行时频重组[16]，使该区域的频率精度得到提高。因此，改进同步挤压小波变换同时提高了重点关注区域瞬时特征参数在频率和时间两个维度上的识别精度。为验证该方法的有效性，将其应用于两个多分量信号数值算例和一个拉索试验。结果表明，改进同步挤压小波变换能够有效提取非平稳响应信号的瞬时频率，是一种可行的时变结构参数识别方法。

1 标准同步挤压小波变换

一个典型的多分量信号可以表达为N个本征函数和一个余量之和。

(1)

其中：xi(t)=Ai(t)cos(φi(t))为第i个本征函数，其幅值的变化率A′(t)远小于相位对时间的变化率φ′(t)；r(t)为噪声或者观测误差。

对于给定的母小波函数ψ(t)，信号x(t)的连续小波变换定义为

(2)

(3)

通过对小波系数求导，初步求得式(4)的瞬时频率为

(4)

(5)

若尺度a及频率ω为连续变量，式(5)可改写为

(6)

2 改进同步挤压小波变换

2.1AMD定理提取多分量信号

(7)

(8)

si(t)=sin(ωbit)H[x(t)cos(ωbit)]-

cos(ωbit)H[x(t)sin(ωbit)]

(i=1,2,…,n-1)

(9)

其中：H[·]为希尔伯特算子。

在进行AMD分解之前，对信号进行连续小波变换，得到其小波量图并以此作为选取截止频率的依据，从而将两个截止频率之间的分量信号提取出来。

2.2 频率转换

标准同步挤压小波变换提高了信号瞬时频率的频率精度，却没有提高其时间精度。为了获取更高的时间分辨率，频率转换把信号的瞬时频率从低频区域转到高频区域。

单分量信号x(t)=cos[2π(f1+f0)t+sin2πf0t]可改写为

x(t)=cos[2πf1t+sin2πf2t]cos(2πf0t)-

sin[2πf1t+sin2πf2t]sin(2πf0t)

(10)

根据式(10)，给定信号x(t)的频率转换定义为

x*(t)=x(t)cos(2πf0t)-H[x(t)]sin(2πf0t)

(11)

其中：H[·]为希尔伯特算子。

信号x*(t)的瞬时频率ω*(a,b)为

ω*(a,b)=ω(a,b)+ω0

(12)

其中：ω(a,b)为信号x(t)通过时频分析方法得出的瞬时频率；ω0=2πf0。

在改进同步挤压小波变换算法中，为得到更好的时间分辨率，需预先设定f0>0，即正频率转换。若f0<0，则为负频率转换，目的是为了得到更好的频率分辨率，将信号的中心频率从高频转到低频。

2.3 时频重组

在时频重组阶段，改进同步挤压小波变换在重点关注频率区间[fm,fM]内对时频代表值进行挤压，提高了瞬时频率识别精度，同时也降低了计算量。定义两个中间变量lfm,lfM分别为

(13)

通过改进同步挤压小波变换得到的频率序列为

(14)

其中：n为重点关注区域离散频率点的个数。

得到信号的圆频率序列为

用ωis替换式(5)的ωl，得到同步挤压小波变换系数值Tx(ωis,b)为

(16)

2.4 改进同步挤压小波变换算法流程

改进同步挤压小波变换算法通过AMD定理分解多分量信号、频率转换和时频重组，同时提高了重点关注区域瞬时特征参数在频率和时间两个维度上的识别精度，其计算流程如图1所示。

图1 改进同步挤压小波变换算法流程图Fig.1 The flowchart of the improved synchrosqueezing wavelet transform algorithm

由图1可知，改进同步挤压小波变换算法共分为5个步骤。

1) 判断初始响应信号是否为多分量信号。将响应信号进行连续小波变换并得到其小波量图。从小波量图中可以观察得出分量信号的数目及每个分量信号瞬时频率的大致范围，并以此作为选取截止频率的依据。

2) AMD定理提取分量信号。在选取合适的截止频率之后，可通过AMD及其拓展定理将多分量信号中的各个分量信号解析地提取出来。

3) 确定挤压范围。单分量信号的挤压范围可根据小波量图中的频率分布情况以及分量信号的截止频率进行选取。

4) 频率转换。将提取出的单分量信号进行周期延拓以减小端点效应对信号的影响，然后确定转换频率，并对周期延拓后的信号进行正频率转换以提高时域分辨率。

5) 时频重组。对完成频率转换步骤的单分量信号进行连续小波变换，将小波变换系数进行时频重组得到改进同步挤压小波变换系数并提取其瞬时频率，根据式(12)将瞬时频率转换到原来位置。

3 数值算例验证

为验证改进同步挤压小波变换的有效性和准确性，考虑多分量调频(frequency modulated,简称FM)信号，即

x(t)=x1(t)+x2(t)

(17)

其中：信号采样频率为1 kHz；采样时间为1 s；x1=sin[2π(100t+t2)]；x2=sin[2π(150t+20t3)]。

两个分量信号所对应的瞬时频率理论值分别为f1(t)=100+2tHz，f2(t)=150t+60t2Hz。为考虑信号噪声影响，对信号添加20%水平的高斯白噪声，噪声强度由信噪比定义(单位dB)。

(18)

其中：Asignal和Anoise分别为信号和噪声的均方根值；噪声水平是指Asignal与Anoise之间的比值。

添加20%水平噪声后的多分量调频信号如图2所示。选用复Morlet小波对含噪多分量调频信号进行连续小波变换，得到如图3所示的小波量图。由图3可知，分量信号x1(t)的瞬时频率集中在100 Hz附近，且呈线性变化，而分量信号x2(t)的瞬时频率分布范围为150～250 Hz，呈抛物线变化趋势。

图2 多分量调频信号Fig.2 Simulated multi-component FM signal with 20% Gauss white noise

图3 添加20%水平噪声的多分量调频信号小波量图Fig.3 Wavelet scalogram of the multi-component FM signal with 20% Gauss white noise

采用AMD定理提取的两个分量信号及其相应的理论值如图 4所示。由于20%水平噪声的影响，使提取的信号与理论值产生一定的偏差，但二者总体上仍十分吻合。设定f0=100Hz，对提取的分量信号进行改进同步挤压小波变换，识别的瞬时频率如图5所示。为验证改进同步挤压小波变换识别瞬时频率的准确性，图5同时给出了多分量调频信号瞬时频率理论值和基于标准同步挤压小波变换的识别值。可以看出，相比标准同步挤压小波变换，改进同步挤压小波变换识别的瞬时频率曲线较为光滑，且与理论值更加吻合，具有更好的分辨率。但是改进同步挤压小波变换以同步挤压小波变换为基础，因此端点效应依然存在，但所受影响要小得多。

图4 AMD提取的单分量信号Fig.4 Mono-components extracted by AMD theorem

图5 改进同步挤压小波变换识别的瞬时频率Fig.5 IF identified by improved synchrosqueezing wavelet transform

为量化瞬时频率的识别精度，采用瞬时频率在整个时间历程内的均方根值作为精度指标(index of accuracy,简称IA) 。IA值越小，说明识别值与理论值越接近。

(19)

其中：fd(t)为瞬时频率识别值；fe(t)为瞬时频率理论值。

表1为同步挤压小波变换法和改进挤压小波变换方法识别瞬时频率的精度指标IA。IAa和IAb分别为信号x1与x2的瞬时频率识别精度指标。由表1可知，针对多分量调频信号，改进同步挤压小波变换的识别效果优于基于标准同步挤压小波变换的瞬时频率识别方法。

表1 多分量调频信号瞬时频率识别精度指标IATab.1 IA of instantaneous frequency identification of multi-component FM signal %

为进一步验证改进同步挤压小波变换识别多分量信号的有效性和准确性，定义一个由两个瞬时频率呈余弦变化的调频分量组成的多分量非平稳信号y(t)为

y(t)=y1(t)+y2(t)

(20)

其中：y1(t)=5cos[2πt+0.8sin(0.4πt)]；y2(t)=10cos[3.6πt+0.5sin(0.5πt)]。

两个分量信号的瞬时频率分别为f1(t)=1.0+0.16cos(0.4πt)Hz，f2(t)=1.8+0.125×cos(0.5πt)Hz。采样频率为100 Hz，采样时间为50 s。添加20%水平噪声后的y(t)如图6所示。

图6 多分量非平稳信号Fig.6 Simulated multi-component non-stationary signal

选用复Morlet小波对含噪多分量非平稳信号进行连续小波变换，得到如图7所示的小波量图。可见，分量信号y1(t)和y2(t)的瞬时频率大致呈正弦变化。采用AMD定理提取的单分量信号如图8所示。图8中提取的单分量信号在整个时间历程内均与理论值基本一致，但是由于噪声的影响，使得提取值在某些峰值处出现不吻合的情况。

图7 添加20%水平噪声的多分量非平稳信号小波量图Fig.7 Wavelet scalogram of multi-component non-stationary signal with 20% Gauss white noise

图8 AMD提取的单分量信号

Fig.8 Mono-components extracted by AMD theorem

设定f0=30Hz，对AMD提取的单分量信号采用改进同步挤压小波变换算法识别，如图9所示。可见，改进同步挤压小波变换识别的瞬时频率与理论值基本保持一致，只是在峰值处产生些许误差。与标准同步挤压小波变换相比，改进同步挤压小波变换识别结果更为光滑和准确。改进同步挤压小波变换算法的端点效应依然存在，但影响较小。

图9 改进同步挤压小波变换识别的瞬时频率Fig.9 IF identified by improved synchrosqueezing wavelet Transform

基于两种方法的瞬时频率识别精度指标IA值如表2所示。由表2可知，改进同步挤压小波变换瞬时频率识别结果明显优于标准同步挤压小波变换识别结果。

表2 多分量非平稳信号瞬时频率识别精度指标IATab.2 IA of instantaneous frequency identification of multi-component non-stationary signal %

4 拉索试验验证

为验证改进挤压小波变换的有效性和准确性，设计一个时变拉索结构试验。试验拉索为一根7φs5的钢绞线。两锚固点间的索长为4.55 m，将加速度传感器竖向安装在拉索中部。索的一端用反力架锚固，另一端采用电液伺服加载系统的作动器施加拉力。首先，对索施加一定的预拉力，随后连续改变索的拉力，使索的刚度随时间发生变化，从而导致索的固有频率随之改变。在改变索力的同时，用力锤敲击拉索，采集索的竖向加速度冲击响应，采样频率为600 Hz。试验装置如图10所示。试验过程中采用了两个工况：即索拉力线性变化和正弦变化。为了与识别结果进行比较，笔者采用“冻结法[11]”，即假定在很小的时间间隔内结构参数保持不变，通过求解系统振动方程的特征值和特征向量近似得到拉索瞬时频率的理论值。本次试验分析中，频率主要采用一阶模态频率，即基频。

图10 试验装置图Fig.10 Cable test setup

4.1 拉力线性变化时索的瞬时频率识别

试验时设定索的拉力从20 kN开始，以1.67 kN/s的速率线性加载，同时采集索的冲击加速度响应，采样时间为7 s，采样频率为600 Hz。测得的加速度响应信号如图11所示。索的拉力变化曲线如图12所示。

图11 拉力线性变化时索的加速度响应Fig.11 Measured cable acceleration responses with linearly varying cable tension force

图12 实测的线性变化拉力Fig.12 Measured cable tension forces with linear variation

如图13所示，对实测加速度响应信号进行复Morlet小波变换，发现实测信号的小波脊线大致呈线性变化。如图14所示，采用AMD定理提取实测信号的基频分量。设定f0=30 Hz，对基频分量进行改进同步挤压小波变换，识别的拉索瞬时频率结果如图15所示。可见，相比标准同步挤压小波变换，改进同步挤压小波变换识别的瞬时频率曲线拥有更高的光滑程度，且与瞬时频率理论值更为接近，只是在信号的开始和结束阶段出现一些偏离。这主要是因为端点效应所致。此外，由于加速度响应信号末端幅值较小，因而信噪比低，在一定程度上影响了瞬时频率识别结果。表3中的IA结果也表明，改进同步挤压小波变换方法识别瞬时频率的准确性相对较高，优于标准同步挤压小波变换。

图13 实测响应信号小波量图Fig.13 Wavelet scalogram of measured cable acceleration responses

图14 AMD提取的基频分量信号Fig.14 The fundamental modal frequency component signal extracted by AMD theorem

图15 拉力线性变化时索的瞬时频率识别结果Fig.15 Identified instantaneous frequency with linearly varying cable tension force

表3 拉力线性变化时索的瞬时频率识别精度指标IATab.3 IA of instantaneous frequency identification with linearly varying cable tension force %

4.2 拉力正弦变化时索的瞬时频率识别

试验时索的拉力呈正弦曲线变化，变换幅度为±4 kN。实测加速度响应信号和索拉力变化曲线分别如图16，17所示。对实测加速度响应信号进行复Morlet小波变换，其小波量图如图18所示。由图18可知，实测信号的小波脊线大致呈正弦变化。通过AMD将实测信号基频分量提取出来，如图19所示。设定f0=20Hz，将改进同步挤压小波变换与标准同步挤压小波变换瞬时频率识别结果进行比较，如图20所示。从图20可知，由于端点效应以及低信噪比的影响，改进同步挤压小波变的瞬时频率识别结果出现了一定的偏差，但瞬时频率的正弦变化趋势与理论频率基本保持一致。两种方法的瞬时频率识别精度指标IA如表4所示。由表4可知，改进同步挤压小波变换精度指标值比同步挤压小波变换的相应值低1.2%，再次证明改进同步挤压小波变换提高了瞬时频率的识别精度。

图16 拉力正弦变化时索的加速度响应Fig.16 Measured cable acceleration responses with sinusoidal varying tension force

图17 实测的正弦变化拉力Fig.17 Measured cable tension force with sinusoidal variation

图18 实测响应信号小波量图Fig.18 Wavelet scalogram of measured cable acceleration responses

图19 AMD提取的基频分量信号Fig.19 The fundamental modal frequency component signal extracted by AMD theorem

图20 瞬时频率识别值比较Fig.20 Identified instantaneous frequency with sinusoidal varying cable tension force

表4 拉力正弦变化时索的瞬时频率识别精度指标IATab.4 IA of instantaneous frequency identification with sinusoidal varying cable tension force %

5 结论

1) 改进同步挤压小波变换能够有效地识别噪声干扰下多分量信号的瞬时频率，且识别效果优于标准同步挤压小波变换。

2) 改进同步挤压小波变换能够有效地识别实测响应信号的瞬时频率。当索的拉力线性变化时，识别出的瞬时频率与理论结果非常接近；当拉力正弦变化时，识别出的瞬时频率与理论识别结果也比较吻合。

3) 端点效应在一定程度上影响基于改进同步挤压小波变换的瞬时频率识别效果，但影响小于标准同步挤压小波变换。

[1] Wang Zhongren,Lin Jun.Time-frequency analysis for complex time-varying signals based on Wigner-Ville distribution[J].Acta Electronica Sinica,2005,33(12):2239-2241.

[2] Hou Z,Hera A,Shinde A.Wavelet-based structural health monitoring of earthquake excited structures[J].Computer-Aided Civil and Infrastructure Engineering,2006,21:268-279.

[3] Huang N E,Shen Z,Long S R.A new view of nonlinear water waves:the Hilbert spectrum[J].Annual Review of Fluid Mechanics,2003,31(1):417-457.

[4] 熊飞,程远胜,刘均.基于HHT方法的时变多自由度系统的损伤识别[J].振动、测试与诊断，2008，28(2)：122-125.

Xiong Fei,Cheng Yuansheng,Liu Jun.Damage identification of time varying multi-degrees of freedom systems based on Hilbert-Huang transform[J].Journal of Vibration,Measurement &Diagnosis,2008,28(2):122-125.(in Chinese)

[5] Wu Zhaohua,Huang N E.Ensemble empirical mode decomposition:a noise-assisted data analysis method[J].Advances in Adaptive Data Analysis,2011,1(1):1-41.

[6] Wang Zuocai,Ren Weixin,Liu Jingliang.A synchrosqueezed wavelet transform enhanced by extended analytical mode decomposition method for dynamic signal reconstruction[J].Journal of Sound and Vibration,2013,332:6016-6028.

[7] Chen Genda,Wang Zuocai.A signal decomposition theorem with Hilbert transform and its application to narrowband time series with closely spaced frequency components[J].Mechanical Systems and Signal Processing,2012,28(2):258-279.

[8] Yan B,Ayaho M.A comparative study of modal parameter identification based on wavelet and Hilbert-Huang transforms [J].Computer-Aided Civil and Infrastructure Engineering,2006,21:9-23.

[9] 孙鹏,丁幼亮,李爱群，等.利用Morlet小波变换识别悬索桥模型模态参数[J].振动、测试与诊断，2012， 32(2)：238-243.

Sun Peng,Ding Youliang,Li Aiqun,et al.Modal identification for suspension-bridge model using morlet wavelet transform[J].Journal of Vibration,Measurement &Diagnosis,2012,32(2):238-243.(in Chinese)

[10] 许鑫，史治宇.状态空间下基于小波变换的时变系统参数识别[J].振动工程学报，2010，23(4)：415-419.

Xu Xin,Shi Zhiyu.Parameter identification of time-varying system based on state space and wavelet transform method [J].Journal of Vibration Engineering,2010,23(4):415-419.(in Chinese)

[11] 王超，任伟新，黄天立.基于复小波变换的结构瞬时频率识别[J].振动工程学报,2009,22(5):492-496.

Wang Chao,Ren Weixin,Huang Tianli.Instantaneous frequency identification of a structure based on complex wavelet transform[J].Journal of Vibration Engineering,2009,22(5):492-496.(in Chinese)

[12] Daubechies I,Lu Jianfeng,Wu H T.Synchrosqueezed wavelet transforms:an empirical mode decomposition-like tool[J].Applied and Computational Harmonic Analysis,2011,30:243-261.

[13] 刘景良，任伟新，王佐才，等.基于同步挤压小波变换的结构瞬时频率识别[J].振动与冲击，2013，32(18)：37-42.

Liu Jingliang,Ren Weixin,Wang Zuocai,et al.Instantaneous frequency identification based on synchrosqueezing wavelet transformation[J].Journal of Vibration and Shock,2013,32(18):37-42.(in Chinese)

[14] Li Chuan,Liang Ming.Time-frequency signal analysis for gearbox fault diagnosis using a generalized synchrosqueezing transform[J].Mechanical Systems and Signal Processing,2012,26 (1):205-217.

[15] 汪祥莉，王斌，王文波，等.混沌干扰中基于同步挤压小波变换的谐波信号提取方法[J].物理学报，2015，64(10)：11-20.

Wang Xiangli,Wang Bin,Wang Wenbo,et al.Harmonic signal extraction from chaotic interference based on synchrosqueezed wavelet transform[J].Acta Physica Sinica,2015,64(10):11-20.( in Chinese)

[16] Cao Hongrui,Xi Songtao,Chen Xuefeng,et al.Zoom synchrosqueezing transform and iterative demodulation:Methods with application [J].Mechanical Systems and Signal Processing,2016,72-73:695-711.