王 超， 毛 羚

(1.湖北工业大学土木建筑与环境学院 武汉， 430034) (2. 武汉理工大学土木工程与建筑学院 武汉，430070)

引 言

对于土木工程结构损伤诊断、健康监测与安全评估来说，模态参数识别是一个非常重要的研究课题。最近几十年，各国学者提出了许多模态参数识别技术，如频域分解法、时域分解法、随机子空间识别、基于响应传递比的参数识别法[1]和基于贝叶斯理论的模态参数识别法[2]等。这些参数识别方法都是针对时不变系统。然而，实际土木工程结构在运营过程中由于环境或者荷载变化可能导致结构发生时变，如不同时刻构件不均匀温度变换会引起结构模态频率随时间改变；桥上车辆荷载的变化引起结构拉索索力的变化，导致拉索刚度时变；车辆过桥时由于车辆位置不断改变引起车桥系统质量时变。因此许多学者提出了时变参数识别方法进行系统识别，如状态空间小波识别法[3]、S变换法[4]、HHT(Hilbert-Huang transform, 简称HHT)法[5]、自适应稀疏时频分析法[6]、同步挤压小波变换法[7]和解析模态分解法(analytical modal decomposition简称AMD)[8]等。

变分模态分解(variational mode decomposition,简称VMD)是最近提出的一种新的信号自适应分解方法[9],与HHT类似，它可以将非平稳的时变信号自适应的分解成多个具有有限带宽的模态分量信号之和，但是它没有HHT分解时遇到的模态混叠问题。由VMD分解得到模态分量需再采取合适方法提取信号的瞬时频率。HHT是采用Hilbert变换的方法进行分析得到Hilbert谱，但是它会产生负频率和端点效应问题。利用小波变换技术也可以提取信号瞬时频率，通常采用的复Morlet小波并不是完全解析的，当中心频率较小时可能产生负频率泄露，导致时频分析结果产生偏差。

因此，笔者采用VMD和广义Morse小波(generalized Morse wavelets，简称GMW)相结合的技术来识别时变结构瞬时频率。首先,利用信号在小波时频域具有稀疏性的特性，引入快速阈值迭代技术对信号进行去噪处理;然后,采用VMD将去噪信号分解为多个模态分量;最后,采用GMW对每个模态分量进行时频分析提取对应的瞬时频率。通过数值算例和一个移动小车试验验证了提出的方法。

1 快速阈值迭代降噪

现场实测信号通常会受到噪声污染,这里引入快速阈值迭代算法(fast iterative shrinkage-thresholding algorithm, 简称FISTA)进行降噪。含噪信号通过小波变换可以转化到时间-尺度域，噪声能量将分布到非常广泛的时间-尺度空间，而有效信号能量将集中在一个局部的时间-尺度域空间，具有稀疏特性，因此，通过最优稀疏性可以实现将噪声和信号尽量分离。假设实测的含噪信号y表示为

y=x+n=Iw(C)+n

(1)

其中:x为真实无噪信号;C为x小波变换系数矩阵;Iw(·)为小波逆变换;n为噪音。

稀疏性对应于l0范数最小化问题，对于非凸复杂性优化问题，可以采用l1范数来替代求解[10],因此，上述问题可以转换为求解优化问题

(2)

其中：s.t.为约束条件;ε为小量。

式(2)有约束优化问题可以通过施加拉格朗日乘子正则化方法转变为无约束优化问题来求解[11]

(3)

其中:λ为拉格朗日乘子。

这里采用Beck等[12]提出的快速阈值迭代算法(IST)通过下式快速迭代求解式(3)的解

(4)

其中：WT(·)为小波变换;Tλ(·)为软阈值函数;k为迭代数;t为迭代参数。

其中：c为小波系数矩阵C中的小波系数。

为了提高收敛效率，λ按指数衰减变化取值

(7)

其中：N为迭代求解次数;λmax=0.99max|C0|;λmin=max|C0|/q，q为系数C0的数量。

通过上述算法对响应数据进行分析，可以有效降低噪音干扰。

2 变分模态分解

VMD是一种新的信号自适应分解技术，将实信号分解为具有最优特定稀疏特性的多个分量信号(模态)之和。这里的稀疏特性采用分解信号的频谱带宽来衡量， VMD要求分解的每个模态分量具有围绕各自中心频率最紧的带宽。

为了估计每个模态分量的带宽，VMD利用希尔伯特变换计算各模态对应的解析信号, 得到其单边频谱, 然后通过频率平移将每个模态频谱平移到各自估计的中心频率位置，最后计算频移后信号梯度范数的平方来估计各模态分量的带宽。这样最优模态分解转化为如下有约束的变分优化问题来实现[13]

(8)

其中：{si}和{ωi}分别为分解的模态分量和对应的中心频率;δ(t)为狄拉克函数;*为卷积;f为分解的信号。

通过引入二次罚函数项和拉格朗日乘子ζ来处理上式的重构约束条件，拉格朗日乘子ζ用来施加严格的约束条件，而二次罚函数项可以提高收敛性能，将式(8)的约束优化问题转变为式(9)无约束的函数优化问题

(9)

其中：α为带宽约束参数。

α增加可能导致分解的模态分量带宽减小，中心频率误差估计增大；α减小可能导致分解的模态带宽增加，包含更多的噪音成分。α建议取值区间为100～10 000, 文中取2 000。

通过交替方向乘子法(alternating direction method of multipliers,简称ADMM)可迭代求取原变分优化问题的解。步骤如下：

(10)

其中：F-1(·)表示傅里叶逆变换。

(11)

4) 按式(12)进行迭代更新

(12)

循环迭代步骤2～4直到满足收敛条件

(13)

其中:ζ为迭代收敛容差，文中取1e-7。

3 广义Morse小波瞬时频率识别

对于N自由度的时变结构系统，其振动响应z(t)可分解为N阶振型模态响应的叠加和

(14)

其中:zm(t)为第m阶振型模态响应。

每阶模态响应可表示为时变幅值Am(t)和相位φm(t)的乘积

xm(t)=Am(t)eiφm(t)

(15)

通过VMD对结构响应信号z(t)进行分析，可以将结构不同阶振型模态分离成各模态分量，然后分别对各模态分量采用连续小波变换进行时频分析，提取对应的小波脊线，由小波尺度与信号的频率关系可以得到结构的瞬时频率[7]。

连续小波时频分析通常采用复Morlet小波，然而，复Morlet小波只是当中心频率较大时近似解析，当中心频率较小时可能产生负频率泄露，导致分析结果产生一定偏差。文中引入GMW[13]进行时频分析, 该小波是完全解析的，不存在负频率泄露问题。GMW频域定义为

Ψβ, γ(ω)=H(ω)2(eγ/β)β/γ)ωβe-ω γ

(16)

其中：H(ω)为阶跃函数;β和γ为控制小波时域和频域衰减的参数。

当γ=3时，GMW接近高斯函数，具有很好的时频分辨率。文中采用的GMW参数为γ=3和β=27，对应的时域和频域波形如图1所示(示意图，幅值无量纲)。

图1 GMW时域和频域波形Fig.1 Wave form of GMW in time domainand frequency domain

4 数值算例

为验证提出的算法, 首先采用如下具有不同时变频率成分组合成的振动信号进行数值验证

(17)

信号采样频率为100 Hz，并添加10%的高斯白噪声。这里噪声强度10%定义为噪声的2-范数与信号的2-范数之比。

首先,利用FISTA对含噪信号进行降噪处理, 降噪前后信号的时频分布图(又称小波量图，尺度坐标为指数坐标，无量纲)如图2所示。

图2 信号小波量图Fig.2 Signal wavelet scalogram

图中颜色亮度表示小波系数幅值大小，越亮表示越大，对应能量就越大。可以看出，信号成分在小波域分布具有稀疏性，比较集中，能量较大，而噪声分布较广，通过FISTA处理后噪声明显减少。

随后,采用VMD对降噪后的信号进行分析，分解得到如图3所示的4个模态分量，原信号各成分也显示在图中用于比较。

图3 VMD分解的各模态分量(黑线)和对应的原信号xi(蓝线)Fig.3 The decomposed mode components by VMD (dark line) and corresponding signal xi (blue line)

从图可以看出,VMD将原信号不同成分完好的分解开来，其中第2个成分x2具有突变特性，VMD将其突变前后分解成两个部分，分解的各部分与原信号对应各成分吻合的非常好。

采用GMW对VMD分解的各模态分量进行连续小波分析提取对应的小波脊线，识别其时变频率。另外，也采用HHT的方法的识别结果用于对比，原信号的瞬时频率的理论值用来验证识别结果的准确性，分析结果如图4所示。

图4 信号各成分时变频率识别结果Fig.4 The identified time-varying frequency of each component

由图4可知，HHT方法识别的3个信号成分瞬时频率受噪音干扰波动较大，对于信号突变成分突变时刻存在一定的偏差，而且HHT方法的端点效应影响很大。提出的方法3个信号成分瞬时频率识别结果与理论值吻合的很好，只在端部由于小波端点效应稍有偏离，但是远小于HHT的影响。

5 移动荷载试验

5.1 试验模型

为进一步验证提出的方法，设计了一个移动小车试验，采用马达牵引小车通过一个简支钢板梁来模拟移动荷载作用下的结构时变特性。板梁尺寸为1 000 mm×150 mm×9 mm，质量为10.555 kg，采用圆钢杆支撑模拟的简支支撑。为避免小车进入和离开主梁时产生大的振动与冲击，在板梁首尾分别连接了一段过渡的引梁和尾梁，但它们与主梁断开。小车设有两轴，轴距为0.117 m，质量为2.308 kg，小车内通过砝码调整重量, 本试验时小车总重为5.296 kg。主梁上设置了一根与主梁中线平行的穿过小车的的钢丝，使小车能沿板梁中线移动。试验模型如图5所示。

图5 试验模型图Fig. 5 The experimental mode

在梁底板起点和终点，1/4,1/2,3/4跨处布设了应变片，在节点3,5,7,9处(见图6)布设了加速度计，测点布置如图6所示。

图6 测点布置图Fig. 6 The layout of measure points

5.2 模态试验

首先进行模态试验。在节点4处用力锤锤击主梁，采集结构的加速度响应，采样频率为1 000 Hz。进行模态分析并对结构的有限元模型进行校正。随后将小车当做移动质量块，利用校正后的有限元模型计算小车重心通过不同节点位置时结构的第1阶模态频率，结果如表1所示。

表1 小车重心通过不同节点时结构第1阶频率(FEA)Tab.1 The first frequency of structure when the vehicle cross different node (FEA)

5.3 移动小车试验

通过马达牵引小车匀速通过试验梁，小车移动时导致主梁上质量分布发生改变, 引起结构的时变。采集响应的应变和加速度响应，采样频率为200 Hz。节点1和11处测得的应变响应如图7所示。

图7 节点1和节点11处测点应变响应Fig.7 Response of strain by sensor at node 1 and node11

通过首尾节点应变测点响应、主梁长度、小车轴距可以计算小车通过主梁时速度。应变峰值对应小车前后轮正好经过测点位置，由图7测得前轮通过节点1的时刻为4.925 s，后轮通过节点11的时刻为7.895 s。由此估算得到小车速度为0.376 m/s。

采用提出的方法对节点9处加速度测点测得的响应进行分析识别结构的第一阶瞬时频率。为进行比较，也分别采用HHT和复Morlet连续小波变换方法(continuous wavelet transform,简称CWT)进行时变频率识别，结果如图8所示。

图8 节点9处测点识别的瞬时频率Fig.8 Identified instantaneous frequency by sensor at node 9

从识别的结果可以看出，采用提出了VMD和GMW方法识别的结果比较光滑，性能较好，通常的CWT识别结构波动稍大，主要存在不同模态间干扰和噪音干扰，HHT识别结果最差，端点效应影响很大。

采用GMW方法提取小波脊线时，对识别精度的影响主要是信号小波脊线对应的各阶成分是否被充分激起，另外信号各阶成分频率间距太小，间距过密可能会造成能量的相互干扰。

随后对所有测点(节点3,5,7,9)测得的加速度响应进行分析，首先用FISTA算法进行降噪，然后采用VMD和GMW相结合的方法识别结构的第1阶瞬时频率。由应变估计的小车移动速度和通过主梁的起止时间，可以推导出小车在不同时刻在主梁上的位置。结合表1中有限元计算的结果，可以得到小车通过主梁不同时刻的第1阶瞬时频率有限元计算结果，与提出的方法进行对比分析，如图9所示。

图9 所有测点识别的结构瞬时频率Fig.9 Identified instantaneous frequency by all sensors

由图可见，所有测点数据的识别结果与有限元计算结果吻合的较好，主要在两端位置稍有偏差，分析偏差的原因:a.尾梁跟主梁连接位置还是存在一定冲击振动干扰;b.小波端点效应有一定影响。另外，所有测点识别结果相差很小，表明提出的方法可以有效识别结构时变参数。

6 结束语

土木工程结构在运营过程中通常表现出时变特性，笔者提出了基于VMD和GMW的方法识别时变结构瞬时频率。对于实测信号的噪声污染问题，提出采用快速阈值迭代算法(FISTA)进行降噪处理。数值模拟结果表明，提出的FISTA算法通过将信号转换到小波域利用信号稀疏特性可以较好地对信号进行降噪，VMD和GMW法可以有效准确地识别时变结构瞬时频率。移动小车试验进一步验证了方法的有效性和精度，提出的方法比HHT和CWT具有更好的效果。