基于同步提取和时频系数模极值的瞬时频率识别∗

2021-06-26刘景良王新宇郑锦仰骆勇鹏

振动、测试与诊断 2021年3期

刘景良，王新宇，郑锦仰，2，盛叶，骆勇鹏

（1.福建农林大学交通与土木工程学院福州，350002）（2.福建省建筑科学研究院福州，350025）

引言

环境激励下的在役土木工程结构本质上是时变和非线性结构，其响应信号呈现非平稳性。时频分析方法（time frequency analysis，简称TFA）能同时在时域和频域提取信号瞬时特征，是分析非平稳信号强有力的工具［1‐2］。其主要方法有Wigner‐Ville 分布［3］、短时傅里叶变换（short time Fourier trans‐form，简称STFT）［4］和小波变换［5］等。

采用时频分析方法识别时变结构响应信号瞬时频率需要解决的一个关键问题是脊线提取，而最常见的脊线提取方法是时频系数模极大值法。如Liu等［6］提出基于小波系数模局部极大值的小波脊线提取方法并将之应用于干涉图分析，该方法的主要优点是运算简单、计算量小，但易受噪声和端点效应的干扰［7］。为此，Carmona 等［8］提出基于随机走动的疯狂爬坡算法来提取信号各分量的小波脊线，但该方法的计算效率不高。Wang 等［9］提出基于动态规划的非平稳信号瞬时频率提取方法在一定程度上消除了端点效应，但其算法较为复杂。近年来，一些后处理方法被引入以解决噪声等外部因素对信号时频脊线提取的影响，如重新分配法［10‐11］、同步挤压小波变换［12‐14］、参数化TFA［15‐17］、解调TFA［18］和同步提取变换（synchroextracting transform，简称SET）［19］等。其中，SET 可以实现时频系数能量的高度聚集和响应信号的完整重构，但是目前尚未运用于土木工程领域。

通过对时变结构响应信号进行同步提取变换，瞬时频率被锁定在一个范围内，有效避免了随机噪声对瞬时频率的干扰。在此基础上，笔者将同步提取与时频系数模极大值脊线提取算法相结合，提出一种新的联合方法（combined method，简称CM）来精确识别响应信号的瞬时频率。首先，对响应信号进行STFT 得到时频图；其次，通过同步提取算法从时频图中提取一条或多条狭窄的频带，在降低噪声干扰的同时也减小了时频系数模极大值算法的计算量；最后，在限定的频带范围内逐点搜索时频系数模极大值点，并按时间顺序连接起来从而得到精确的时频脊线和瞬时频率曲线。

1 基本理论

1.1 同步提取变换

由于谐波信号经过时频变换得到的时频系数实部的频率与原信号的瞬时频率相等，因此同步提取变换通过计算各时频系数实部的频率来确定信号的瞬时频率（instantaneous frequency，简称IF）［19］。同步提取变换作为一种后处理技术，可以和STFT 等多种TFA 方法结合。以STFT 为例，原信号s(t)的STFT 可表示为

其中：τ为时移因子；ω为频率；t为时间；g(τ-t)为可移动的紧支高斯窗函数。

其中：σ为决定窗宽的尺度因子。

将式（2）代入式（1）可得

设定改进高斯窗φ(ω，t)=

则式（3）可改写为

根据傅里叶变换性质及Parseval 定理，式（4）可表示为

若以谐波信号s(t)=为例，其傅里叶变换可表示为

将式（6）代入式（5）得

对式（7）求偏导可得

至此，原信号s(t)的瞬时频率，即式（8）中的ω0可表示为

由于同步提取变换仅保留时频面中瞬时频率位置的时频系数而令其余位置的时频系数归零，因此可通过对时频系数矩阵乘以平移到ω0(τ，ω)的δ(ω)函数来实现同步提取，如式（10）所示

Te(τ，ω)=SSTFT(τ，ω)δ[ω-ω0(τ，ω)] （10）

其中：δ[ω-ω0(τ，ω)]为克罗内克函数。

考虑到计算误差的存在，将式（11）改写为

其中：Δω为采样频率间隔。

将式（12）代入式（10）得到同步提取的最后结果，如式（13）所示

由于同步提取变换算法涉及到求导、相除等运算，将造成一定的计算误差，而定义频带范围就是为了更好地考虑在瞬时频率提取过程中出现的识别误差。

1.2 基于同步提取和时频系数模极大值的联合方法提取时频脊线和瞬时频率

响应信号经时频变换后其时频图上的某些位置会出现亮点，这表明该处的时频系数模值较大。因此，可通过搜索时频系数模极大值位置来获取时频脊点，然后将脊点按照时间顺序连接起来即为时频脊线。然而，时频系数模极大值算法由于抗噪性较差，其提取的瞬时频率曲线含有大量毛刺，影响了识别精度。基于此，笔者引入同步提取算法并与时频系数模极大值法相结合来解决上述问题。该联合算法的具体步骤如下：

1）根据式（7）对信号进行STFT 得到时频系数矩阵；

2）按照式（9）初步估算出信号的瞬时频率；

3）根据式（13）将时频面上非瞬时频率范围位置上的时频系数归零；

4）令i，j分别表示时频系数矩阵的列和行，初始时令i=1，搜索第i列在限定频带范围内的时频系数最大值，最后令j(i)=j，i=i+1，并开始搜索下一列在限定频带范围内的时频系数最大值；

5）重复步骤4，直至求出时频系数矩阵所有列中限定频带范围内的时频系数模极大值，并将极大值点连接成线，即为待求时频脊线和瞬时频率曲线。

2 数值算例验证

2.1 单分量信号数值算例

考虑如式（14）所示的单分量调幅调频信号

根据瞬时频率定义，x(t)的理论频率为f=25+5cos(πt)。设定信号采样频率为600 Hz，采样时长为10 s。为模拟真实情况中噪声对信号的影响，对上述信号添加50%水平的高斯白噪声，其噪声强度由信噪比（signal noise ratio，简称SNR）定义

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

添加50%水平高斯白噪声后的信号如图1 所示。首先，对信号进行STFT 得到如图2 所示的时频图；其次，对信号x(t)进行同步提取变换，得到的时频曲带如图3 所示；最后，采用时频系数模极大值法在锁定的瞬时频率曲带范围内搜索时频系数模极大值。通过比较图2 和图3 可以看出，同步提取变换将瞬时频率锁定在一个曲带范围内，同时也去除了图2 中所含的毛刺。为方便比较，图4 同时给出了基于CM 法、小波系数模极大值法［6］和同步挤压小波变换方法［5］的瞬时频率识别值以及理论值。图5 为上述3 种方法瞬时频率识别结果的局部放大效果图。图6 给出了3 种方法的相对误差。由图4和图5 可知：笔者提出的CM 方法识别结果与理论值吻合较好，且具有较高的抗噪性；同步挤压小波变换的识别效果次之，但是同步挤压小波变换是在每一中心频率附近区间固定一个频率值来细化时频曲线，同时最终提取的瞬时频率曲线也是时频图中能量最高或曲率最小的部分；该方法虽然提高了频率分辨率，但也降低了时间分辨率［20］；小波系数模极大值法识别的瞬时频率曲线因含有大量毛刺而与理论频率值存在较大的偏差，整体识别效果较差。

图1 添加50%水平高斯白噪声的x(t)Fig.1 x(t)with 50% Gaussian white noise

图2 x(t)的时频图Fig.2 The time-frequency spectrogram of x(t)

图3 同步提取变换识别x(t)的瞬时频率曲带Fig.3 The IF curve band of x(t)extracted by SET

图4 x(t)的瞬时频率识别结果对比Fig.4 Comparison of IF identification results of x(t)

图5 x(t)的瞬时频率识别结果局部放大效果Fig.5 Local magnification of IF identification results of x(t)

图6 x(t)的瞬时频率识别相对误差Fig.6 Relative errors of IF identification of x(t)

2.2 多分量信号数值算例

考虑如式（16）所示多分量调幅调频信号

设定采样频率为600 Hz，采样时间为10 s。瞬时频率理论值分别为f1=10+2.5 cos(0.5πt)和f2=30+2.8 cos(1.4πt)。对信号添加50%水平的高斯白噪声，其含噪信号如图7 所示。对信号进行STFT 得到如图8 所示的时频图。利用同步提取变换去除噪声干扰并锁定瞬时频率的限定范围，结果如图9 所示。图10 同时给出了基于CM、小波系数模极大值法和同步挤压小波变换提取的瞬时频率识别结果以及理论瞬时频率曲线。为更清楚地比较3种方法的瞬时频率识别效果，图11 给出了局部放大结果，图12 则给出了3 种方法的相对误差。通过对比图11（a）和图11（b）可以看出：在识别f1时，由于f1变化缓慢，该信号为慢变信号，CM 和同步挤压小波变换识别的瞬时频率结果均与理论值十分吻合；在识别f2时即瞬时频率相对快变的情况下，CM 法仍具有很好的稳定性且识别结果与理论值吻合较好，但是同步挤压小波变换识别结果由于时间分辨率下降的原因其识别效果相对较差；小波系数模极大值法识别的瞬时频率曲线则含有大量毛刺，识别效果不佳。总之，CM 识别的瞬时频率曲线不仅最平滑而且基本没有毛刺，其识别值也最接近理论值，在3 种瞬时频率识别方法中表现最优。

图7 添加50%水平高斯白噪声的y(t)Fig.7 y(t)with 50% Gaussian white noise

图8 y(t)的时频图Fig.8 The time-frequency spectrogram of y(t)

图9 同步提取变换识别y(t)的瞬时频率曲带Fig.9 The IF curve band of y(t)extracted by SET

图10 y(t)的瞬时频率识别结果对比Fig.10 Comparison of IF identification results of y(t)

图11 y(t)的瞬时频率识别结果局部放大效果Fig.11 Local magnification of IF identification results of y(t)

图12 y(t)的瞬时频率识别相对误差Fig.12 Relative errors of IF identification of y(t)

3 试验验证

笔者从刚度时变和质量时变两个角度来设计时变结构动力试验以进一步验证所提CM 方法的有效性和准确性，并将瞬时频率识别结果与基于小波系数模极大值和同步挤压小波变换的识别值进行了对比。

3.1 刚度线性变化拉索试验

刚度线性变化拉索试验装置如图13 所示。试验所用拉索为1 根4.55 m，7Ø5 的预应力钢绞线拉索，一端用反力架固定，另一端固定于电液伺服加载系统（mechanical testing &simulation，简称MTS）上。在拉索中部竖向安装加速度传感器并采集加速度冲击响应。在初始阶段对拉索施加20 kN 的预拉力后，通过调整MTS 的作动器使索的拉力以1.67 kN/s 的速率线性增加，拉力的变化将引起拉索刚度的改变。在改变拉力的同时，用力锤敲击拉索以采集拉索的加速度冲击响应，采样频率设为600 Hz，采样时间为6 s，所得响应信号如图14 所示。

图13 试验装置图Fig.13 The cable test setup

图14 拉索线性变化时的加速度响应Fig.14 The measured cable acceleration responses with lin‐early varying cable tension force

首先，采用解析模态分解方法［21］对实测响应信号进行分解，得到1 阶分量信号并对其进行STFT得，到如图15 所示的时频图；其次，对信号进行同步提取变换以去除噪声影响并将瞬时频率锁定在一定范围内，结果如图16 所示；最后，在限定范围内提取信号的时频系数模极大值并识别其瞬时频率。为方便比较，图17 给出了CM 法、小波系数模极大值法、同步挤压小波变换的瞬时频率识别结果以及基于“冻结法”［22］求解的瞬时频率理论值。图18 和图19则分别给出了相对误差曲线图和相对误差曲线局部放大图。由图17 可知，CM 法识别的拉索瞬时频率变化趋势为线性，而且最接近理论值。在信号的末端由于随机噪声掩盖了真实响应信号，此时同步挤压小波变换容易出现端点效应，其瞬时频率识别结果偏离理论值较大［5］；而CM 法在信号末端并没有出现毛刺现象和明显的端点效应，其针对实测信号的识别效果较佳。

图15 拉索响应信号的时频图Fig.15 The time-frequency spectrogram of the response from the cable

图16 同步提取变换识别拉索响应信号的瞬时频率曲带Fig.16 The IF curve band of the response from the cable ex‐tracted by SET

图17 拉索响应信号的瞬时频率识别结果对比Fig.17 Comparison of IF identification results of the re‐sponse from the cable

图18 拉索信号的瞬时频率识别相对误差Fig.18 Relative errors of IF identification of the response from the cable

图19 拉索响应信号的瞬时频率相对误差局部放大效果Fig.19 Local magnification of IF relative errors of the re‐sponse from the cable

3.2 质量突变悬臂梁试验

质量突变悬臂梁试验装置如图20 所示。铝合金悬臂梁的质量为0.81 kg，截面尺寸为40 mm×15 mm，长度为500 mm。悬臂端的质量块为1 kg，在质量块上方放置一块悬挂的磁铁。在梁上每隔100 mm 设置一锤击点，然后通过力锤敲击悬臂梁的自由端，并在某个时刻放下细绳使永磁铁垂直靠近质量块并将其吸起以实现结构质量的改变。设定采样频率为2 kHz，在梁的跨中位置安装加速度传感器以采集加速度冲击响应。试验前，基于“冻结法”测得带质量块与不带质量块时悬臂梁的基频分别为21.24 和47.06 Hz。试验过程中，先在悬臂梁末端用力锤施加激振力，并在2 s 附近利用磁铁吸引质量块以改变悬臂梁质量，最后在2.3 s 再次锤击悬臂梁以避免响应信号衰减过快。

图20 铝合金悬臂梁试验装置图(单位：mm)Fig.20 The aluminum cantilever beam test rig (unit:mm)

加速度传感器记录的响应信号如图21 所示。首先，采用解析模态分解定理对实测响应信号进行分解得到1 阶分量信号，再对1 阶分量信号进行ST‐FT 得到如图22 所示的时频图，由图可知频率在2s附近发生突变，与实际情况相符；其次，对响应信号进行同步提取变换，得到如图23 所示的瞬时频率曲带；最后，采用时频系数模极大值法提取实测响应信号的瞬时频率。为方便比较，图24 同时给出了CM法、小波系数模极大值法和同步挤压小波变换法提取的瞬时频率曲线以及理论值。图25 给出了瞬时频率识别结果的局部放大图。图26 和27 分别给出了相对误差曲线和局部相对误差曲线。由图24 和图25 可知：小波系数模极大值法识别的瞬时频率曲线存在大量毛刺；虽然同步挤压小波变换识别的曲线比小波系数模极大值法平滑，但其识别效果仍然不如CM 法，且在2.0 s 附近出现了明显波动。