基于变分模态分解和同步提取变换识别时变结构瞬时频率

2022-03-27黄天立

振动与冲击 2022年6期

唐蕾，黄天立，万熹

(中南大学土木工程学院，湖南长沙 410075)

在实际运营过程中，由于环境和荷载的变化会导致结构出现时变特性，如斜拉桥拉索的索力变化、列车过桥引起的桥梁振动以及环境温度等都会导致结构动力特性随时间而改变。这使得结构具有时变特征，结构响应表现出非平稳特性。因此，开展时变结构的瞬时频率等特征参数识别方法研究，对于结构健康监测具有重要意义。目前国内外对时变结构瞬时频率识别的方法可以分为两类。第一，基于时频分析的方法，即从信号时频谱图中获取时变结构瞬时频率等参数，并进一步通过时频重排、能量重组等方法提高识别精度；第二，基于信号自适应分解的方法，即将多分量信号分解成多个单分量信号，再对各个单分量信号进行瞬时频率等参数识别。

时频分析方法能同时在时域和频域内分析信号的时频特征，是分析非平稳信号的有力工具[1]。传统的时频分析方法，如Wigner-Ville分布、Gabor变换、短时Fourier变换(short-time Fourier transform，STFT)[2]、小波变换[3-4]等，已在时变模态参数识别和非线性系统识别等领域得到了广泛的应用。为提高时频分析方法的精度，发展了针对时频谱进行能量重组的方法。Daubechies等[5]提出了同步压缩变换(synchrosqueezing transform，SST)，该方法将小波变换后的小波时频谱图进行重组，将小波系数压缩到中心频率附近获取更高精度的时频表示。刘景良等[6-7]采用同步压缩小波变换对Duffing系统自由振动响应和两层剪切型框架模型在地震激励下的响应进行了结构瞬时频率识别。王超等[8]采用同步压缩小波变换识别了移动车辆荷载作用下的桥梁时变参数。Oberlin等[9]在同步压缩小波变换的基础上，提出了同步压缩Fourier变换(fourier-based synchrosqueezing transform，FSST)。徐晓迪等[10]利用同步压缩短时Fourier变换，提取了高速列车轴箱振动加速度信号的瞬时频率。FSST提高了时频分辨率，但噪声同样也会被压缩，使得该方法的抗噪性能差。针对FSST存在的问题，于刚等[11]提出同步提取变换(synchroextrcting transform，SET)，SET提取STFT谱图在瞬时频率位置上的时频系数，具有较强的时频聚焦性。使用SET识别时变结构瞬时频率可以较好地剔除环境或温度等引起的噪声。应该指出，SET方法是针对STFT的改进，在处理多分量信号时，各分量成分信号的频谱之间须具有一定的距离，才能获得精度较高的时频谱图。因此，寻找适合的信号分解方法，有效分解近距离频率成分分量对于SET的应用至关重要。

信号的自适应分解方法中，Huang等[12]提出的经验模式分解(empirical mode decomposition，EMD)方法影响较大，通过EMD分解得到的固有模式函数(intrinsic mode function，IMF)，结合Hilbert变换可提取多分量信号的瞬时频率等参数。Shi等[13-14]基于EMD和Hilbert变换识别了刚度慢变、周期改变和突变3种类型单自由度时变系统在自由振动情况下和多自由度时变系统在受迫振动情况下的瞬时频率和阻尼系数。EMD通过包络线拟合的方法将多分量信号分解成多个固有模式函数之和，但是在拟合的过程中容易出现欠包络或过包络，从而导致模态混叠。许多学者利用EMD分解的思路发展了许多多分量信号分解方法，如局部均值分解[15]、局部特征尺度分解[16]等。相比EMD分解方法这些方法在一定程度上改善了分解效果，但是由于这些方法仍然都是基于包络线拟合的方法，不能完全解决模态混叠的问题。变分模态分解(variational mode decomposition，VMD)[17]是一种新近提出的信号自适应分解方法，它舍弃了包络线拟合，将基于极值点模态获取问题变为变分模型的构造，将分解方式变为变分模型的求解，解决了模态混叠的问题，同时VMD相当于自适应的维纳滤波器组，具有较强的抗噪性能。结合信号分解方法与单分量信号参数识别方法可以实现结构瞬时频率等参数的识别。赵亚军等[18]利用VMD和经验包络法对密集模态和瞬态系统的瞬时频率和瞬时阻尼比进行了识别。王超等[19]提出基于VMD和广义Morse小波相结合的时变结构瞬时频率识别方法，用具有时变特性的移动小车试验验证了所提方法的有效性。

针对SET不能分离频率成分间隔相近的多分量信号的问题，本文结合VMD和SET，提出了高能量聚集性的时变结构瞬时频率识别方法。首先，通过傅里叶变换确定预设模态数量，利用VMD对多分量信号进行分解处理得到多个模态分量；然后，采用SET对每个模态分量进行时频分析获取瞬时频率；最后，将各个模态分量的时频谱图叠加得到完整的多分量信号时频谱图。该方法结合了VMD和SET的优势，解决了SET处理具有近距离频率成分多分量信号的不足，准确识别了时变结构的瞬时频率。多分量时变信号和两自由度时变结构自由振动响应信号数值算例和时变拉索试验验证了方法的正确性和适用性。

1 基本原理

1.1 变分模态分解

VMD是一种新近提出的信号自适应分解方法，它假设信号的每个模态分量具有不同中心频率且具有围绕各自中心频率最紧的带宽。VMD分解主要基于维纳滤波、希尔伯特变换、频率混合这3个概念，为了使得估计带宽之和最小，将待分解信号引入变分模型，求解变分模型的最优解获取信号的最佳分解。

1.1.1 变分模型构造

原始信号s(t)由K个分量信号构成，通过傅里叶变换等方法确定预设模态数量，VMD假定原始信号s(t)可被分解为K个IMF分量，定义各IMF分量为调频调幅信号uk(t)=Ak(t)cos(φk(t))，构造变分模型。

对uk(t)进行Hilbert变换，得到解析信号

(1)

式中：δ(t)为Dirac函数；j为虚数单位；*为卷积运算。

利用解析信号的频移特性，对其混合一个中心频率e-jωkt，将uk(t)的频谱移动到相应的基频带

(2)

计算频移后信号梯度范数的平方估计各模态分量的带宽，为使各IMF分量的带宽之和最小，建立约束变分模型如下

(3)

式中：{uk}和{ωk}为VMD分解得到的模态分量和对应的中心频率；s(t)为被分解信号。

1.1.2 变分模型求解

步骤1为了求解约束变分模型，引入拉格朗日乘子λ对约束严格控制，再引入二次惩罚因子α提高该变分问题的收敛性，将有约束变分问题转化为无约束变分问题

(4)

(5)

(6)

(7)

式中，τ为迭代步长。

步骤4设定判别精度e>0，若满足收敛条件式(8)，则分解结束，否则循环迭代步骤2和步骤3直到满足收敛条件

(8)

1.2 同步提取变换

SET是一种新近提出的高分辨率时频分析方法，其仅提取STFT谱图在瞬时频率位置上的时频系数，大大减少了噪声的影响，具有优越的抗噪性能；同时SET克服了同步压缩Fourier变换在能量压缩过程中将噪声分量压缩到时频谱图的缺陷。 SET基本步骤如下。

步骤1由STFT计算信号s(t)的时频谱图

(9)

令gω(u)=g(u-t)·ejωu，由于窗函数g通常采用实函数，故窗函数的复共轭等于它本身，根据Parseval定理，式(9)可以写成

(10)

(11)

(12)

(13)

步骤3提取STFT在瞬时频率位置的时频系数。为尽可能减轻噪声影响，采用δ函数仅提取ω=ω0处的时频系数，得到一个具有高分辨率的时频谱Te(t,ω)，

Te(t,ω)=Ge(t,ω)δ(ω-ω0(t,ω))

(14)

式中，δ(ω-ω0(t,ω))为同步提取算子(synchroextracting operator,SEO)。

考虑到计算误差，同时实际应用中要用到SEO的实部，建议用式(15)计算SEO

(15)

式中，Δω为离散频率之间的间隔。

对于多分量信号，STFT时频谱图为

(16)

根据相位信息估计的瞬时频率为

(17)

此时，SET的表达式为

Te(t,ω)=Ge(t,ω)δ(ω-φ′(t,ω))

(18)

1.3 时变结构瞬时频率识别流程

VMD可以有效地分解多分量信号，对于多自由时变结构也适用，解决了SET无法识别频率间隔较近多分量信号的不足，同时VMD+SET方法利用了SET高能量聚集性的优点，可以识别时变结构的瞬时频率。图1给出了基于VMD+SET方法识别时变结构瞬时频率的流程图，其基本步骤如下。

图1 VMD+SET方法识别时变结构瞬时频率流程图

步骤1基于待分析时变结构响应信号的傅里叶幅值谱图确定预设模态数量K。

步骤2采用VMD对多分量信号进行分解得到K个模态分量。

步骤3采用SET对每个模态分量进行时频分析获取瞬时频率。

步骤4将各个模态分量的时频谱图叠加得到完整的多分量信号时频谱图，从而揭示时变结构的时变特性。

2 数值算例

2.1 多分量时变信号

考虑一模拟多分量时变信号s(t)，由3个分量信号s1(t)、s2(t)和s3(t)组成，各分量瞬时频率分别随时间产生线性、突变和二次变化。取信号时长为60 s，采样频率为50 Hz，考虑到实际工作测试中噪声的影响，对信号添加10 dB的高斯白噪声，其时程曲线如图2所示。

图2 模拟的含噪多分量时变信号s(t)

(19)

采用高斯窗函数g(t)=e-πt2/0.322，时宽为2 s，对含噪信号s(t)分别进行STFT、FSST、SET处理得到时频谱图，如图3所示，其局部放大图如图4所示。从图3(a)和图4(a)可以看出，信号经STFT变换后，其时频谱图上的频率轨迹线比较模糊，虽然可以大致看出信号瞬时频率的变化趋势，但是不能准确地确定每个时刻各分量信号对应的瞬时频率。从图3(b)和图3(c)可以看出，信号经过FSST和SET变换后，其时频谱图上的频率轨迹线更为清晰，噪声影响明显减少。对比局部放大图4(b)和图4(c)可以发现SET的时频谱图能量较FSST更为集中，时频分辨率更高，FSST的谱图中的噪点较少，在能量压缩过程中噪声也不可避免地一同压入。由此可见，SET能较好地识别频率随时间线性、突变和二次变化的信号，且比STFT和FSST的能量聚集性强，SET适用于识别时变结构的瞬时频率。

图3 含噪声信号s(t)的时频谱图

图4 含噪声信号s(t)的时频谱图(局部放大图)

进一步采用VMD结合SET的方法识别该多分量信号的瞬时频率。根据傅里叶频谱图确定该多分量信号的模态数量K=4，用VMD对信号进行分解得到4个分量信号，如图5所示。从图5(b)和图5(c)可以看出，VMD分解出了突变信号，但是单独使用VMD无法判断各模态分量的具体成分，结合SET对VMD分解的各模态分量进行瞬时频率的识别，再将各个模态分量的时频谱图叠加得到完整的多分量时变信号的时频谱图，如图6所示。从图6可以看出，对比理论值(点线)，VMD结合SET的方法能准确识别频率随时间产生线性、突变和二次变化的3种时变信号的瞬时频率。

图5 VMD分解的各模态分量

图6 含噪声信号s(t)经VMD+SET处理的时频谱图

采用Rényi熵定量评价STFT、FSST、SET和VMD+SET 4种时频分析方法的时频能量聚集特性。Rényi熵是描述系统混乱度的一种衡量指标，熵值越低则代表时频谱图的能量聚集性更好，其计算公式为

(20)

式中：W(t,ω)为时频系数；α为Rényi熵的次序，本文采用3次Rényi熵进行评价。

对信号s(t)分别施加信噪比(signal-noise ratio ,SNR)从0～20 dB的噪声，计算分别采用STFT、FSST、SET和VMD+SET 4种时频分析方法获取的信号时频图的Rényi熵值，如图7所示。从图7可以看出，基于SET时频谱图计算得到的Rényi熵值较STFT和FSST时频谱图计算得到的Rényi熵值更低，验证了SET具有较好的时频聚焦能力。将信号先用VMD分解再采用SET处理，相当于对信号进行了降噪处理。因此，基于VMD+SET时频谱图计算得到的Rényi熵值最低。由此可见，在评价的4种方法中，VMD+SET方法具有最高的时频能量聚集特性和最优的抗噪性能，可以用于高精度的时变结构瞬时频率识别。

图7 噪声信号s(t)(信噪比0～20 dB)的STFT、FSST、SET和VMD+SET时频谱图的Rényi熵

2.2 两自由度时变结构

图8所示为一两自由度时变结构，调整其刚度使其两阶自振频率接近，不满足SET处理多分量信号的限制条件。在一定的初位移和初速度条件下，结构发生自由振动，振动过程中假设结构质量和阻尼不随时间变化，振动过程中由于损伤结构刚度随时间变化，结构的运动方程为

图8 两自由度时变结构

(21)

式中，m1=m2=1 kg，c1=c2=0.15 N·s/m，k2=300 N/m，刚度k1按式(21)随时间线性变化

(22)

设定初始条件x1(0)=0.3 m，x2(0)=-0.06 m，由四阶Runge-Kutta法求解结构的自由振动响应，其采样频率为100 Hz，信号时长为15 s。考虑测量噪声的影响，对质量块m1的位移响应信号添加10 dB的高斯白噪声，如图9所示。

图9 质量块m1的含噪声位移响应信号

根据“冻结法”的原理，假定短时间间隔内时变结构为时不变结构，通过结构动力学方法求解得到各时刻时不变结构的固有频率，作为此时变结构的瞬时频率理论值。图10(a)～图10(c)给出了对图9所示含噪声位移响应分别进行STFT、FSST、SET处理后的时频谱图，采用高斯窗函数g(t)=e-πt2/0.322，时宽为1 s。从图10可以看出，STFT时频谱图中能量主要集中在4～6 Hz内，两阶频率离得比较近，发散部分产生了交叉重叠，不能满足间隔距离条件，导致两阶模态频率没有清晰地分离开；而FSST和SET作为对STFT 时频谱图的后处理手段，也不能将两阶模态分量清晰地分离开，由此使得低能量的一阶模态频率无法识别。采用VMD+SET方法对该两自由度时变结构的自由振动位移响应进行分析，首先由VMD分解信号得到振幅逐渐衰减的两阶模态分量信号，如图11所示，再利用SET获得时频谱图，如图12所示(图12中实线为理论频率值)。从图12可以看出，VMD+SET方法很好地分离了时变两自由度结构的两阶模态分量，且识别的瞬时频率值与理论频率值基本一致。

图10 位移响应的时频谱图

图11 VMD分解的两阶模态分量

图12 两自由度时变结构自由振动位移响应的VMD+SET时频谱图

3 时变拉索试验验证

为进一步验证VMD+SET方法识别时变结构瞬时频率的有效性，设计一个时变拉索装置，通过改变拉索在不同时刻的张拉力，模拟拉索的刚度时变特性。

3.1 拉索试验装置

图13给出了拉索试验装置的示意图和现场布置图[21]。试验拉索为7φS5一根的钢绞线，水平放置，弹性模量为E=1.95×105MPa，截面积为1.374×10-4m2，线质量密度为1.1 kg/m，两锚固点间拉索总长度为4.55 m。拉索一端锚固于反力架上，一端用电液伺服加载系统(electro hydraulic servo,MTS)作动器对其施加拉力，如图13所示。试验开始时，首先对拉索施加20 kN的预拉力，然后通过MTS作动器调整索力的大小，使得拉索刚度随时间发生变化，从而导致拉索固有频率随之改变。在索力变化过程中，通过力锤敲击拉索，并利用在拉索中部安装的一个加速度传感器，采集拉索的竖向加速度响应。

图13 拉索试验装置

试验研究索力线性和正弦变化两种工况，对比基于MTS实测索力换算得到的拉索理论频率，验证基于VMD+SET的时变结构瞬时频率识别方法识别拉索瞬时频率的正确性。

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

试验开始时，拉索预拉力为20 kN，拉索拉力通过MTS作动器以1.67 kN/s的速率线性增长，MTS作动器同步记录拉力变化曲线，如图14所示。索力变化过程中，采用力锤敲击拉索后，同时利用拉索中部安装的加速度传感器采集拉索竖向冲击加速度响应，如图15所示。响应时长6 s，采样频率为600 Hz。

图14 MTS实测拉力(线性)

图15 索力线性变化时拉索的加速度响应

对实测加速度响应信号采用VMD+SET方法识别拉索的瞬时频率，采用时宽为0.5 s的高斯窗函数g(t)=e-πt2/0.322，图16给出了识别得到的拉索第1阶模态瞬时频率以及线性拉力作用下，基于MTS实测拉力数据换算得到的拉索第1阶模态瞬时频率理论值(虚线)。从图16中可以看出，拉索的瞬时频率线性增加，识别的瞬时频率与理论瞬时频率非常接近，识别准确度高。应该指出，由于振动衰减，振幅减小，噪声影响较大，加速度响应信号尾部的VMD+SET时频谱图不太明显，因此识别的瞬时频率值稍有偏差。此外，信号两端由于“端点效应”影响，识别结果与理论结果之间也存在一定的误差。

图16 理论频率与识别频率结果(线性)

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

试验开始时，拉索预拉力为20 kN，拉索拉力通过MTS作动器按正弦规律变化，变化幅度为±4 kN，MTS作动器同步记录拉力变化曲线，如图17所示。索力变化过程中，采用力锤敲击拉索后，同时利用拉索中部安装的加速度传感器采集拉索竖向冲击加速度响应，如图18所示。响应时长6 s，采样频率为600 Hz。

图17 MTS实测拉力(正弦)

图18 索力正弦变化时拉索的加速度响应

对实测加速度响应信号采用VMD+SET方法识别拉索的瞬时频率，采用时宽为0.5 s的高斯窗函数g(t)=e-πt2/0.322，图19给出了识别得到的拉索第1阶模态瞬时频率以及正弦变化拉力作用下，基于MTS实测拉力数据换算得到的拉索瞬时频率理论值(虚线)。从图19中可以看出，识别的拉索瞬时频率略低于理论频率值，与理论值稍有偏差，能明显看出具有正弦变化规律。同样由于振动衰减，以及“端点效应”影响，识别结果与理论结果之间也存在一定的误差。