秦世强,蒲黔辉,施 洲
(西南交通大学 土木工程学院,成都 610031)
模态参数识别是桥梁健康监测的重要组成部分,准确的识别模态参数,是进行有限元模型修正、结构损伤识别以及性能评定的前提。桥梁模态参数识别是指通过测得的动力响应,识别桥梁的固有振动频率、阻尼比和振型的过程。环境激励下的模态参数识别由于不需要外部激励设备且不影响正常交通,已经成为大型桥梁模态参数识别的主要方法。目前常用的方法有频域法、时域法和时频分析法[1]。
经验模态分解(Empirical mode decomposition,EMD)是由Huang等人提出的一种适用于非平稳、非线性信号分解的方法,将多成分的复杂信号分解成一系列的单成分的本征模态函数(IMF),突破了以傅里叶变换为基础的稳态分析;随机子空间识别(Stochastic subspace identification,SSI)是一种环境激励下结构模态参数识别的方法,由于只需要确定系统的阶次,且不涉及迭代和收敛性问题,因此受到广泛的关注。文献[2] 利用EMD结合随机减量技术识别了青马桥的模态参数,并验证了EMD在非平稳信号模态参数识别中的优势;文献[3] 利用HHT提取信号的瞬时特性,再结合自然激励技术识别结构模态参数;文献[4] 和文献[5] 研究了随机子空间方法在结构模态参数识别中的应用,并探讨了其系统阶次的确定、虚假模态等问题。大型桥梁自振频率低、模态密集,加之测试过程中的噪声影响,在EMD过程中往往会出现模态混叠;单纯的利用SSI方法识别模态参数又存在虚假模态的问题,因此尝试从解决这两个问题的角度出发,引入屏蔽信号以消除EMD中的模态混叠,将原始信号分解成一系列只包含结构某一阶固有振动信息的IMF,然后再利用SSI识别桥梁的模态参数。
EMD通过一种称为“筛”的算法[6]实现:① 找出原始信号x(t)所有极值点,用三次样条曲线分别拟合极大值点和极小值点,得到x(t)的上下包络线,取其平均值为m1(t),用原始信号x(t)减去m1(t),可得到一个新的数据序列h1(t):
判断h1(t)是否满足IMF的两个条件:① 在整个时间跨度内,极值点和过零点的个数相等或者至多相差一个;② 在任意一点处,由极大值构成的上包络线和由极小值构成的下包络线的均值为零。一般情况下,h1(t)不是一个IMF,此时将h1(t)视作原始信号,重复上述过程 n 次,直至 h1n(t)=h1,n-1(t)- m1n(t)满足IMF的条件,则h1n(t)为原始信号的第一个IMF,记作c1(t),表示原始信号中的最高成分;② 将原始信号x(t)减去c1(t),得到一个新的序列r1(t),对r1(t)重复步骤(1),可以得到一系列频率从高到低的cj(t),和一个不可分解的余量rn(t),它表示原始信号中的趋势项。至此,原始信号x(t)可以表示为:
线性振动系统的随机状态空间方程[4]为:
式中,xk表示系统的状态向量,yk为输出向量,A为n×n的系统的状态矩阵,n表示系统的阶次,C为输出矩阵;wk、vk分别表示过程噪声和测量噪声,均假定为均值为0的白噪声,相关函数为:
式中:E表示数学期望,δpq表示克罗内克函数。由式(3)、式(4)两式构成系统的随机状态空间模型,有多种方法可以实现振动系统识别,而SSI由于只需要确定系统阶次一个参数,且计算过程引入QR分解和SVD分解,一般不涉及解的收敛性问题,因而广受关注。SSI的主要目标是求解系统的状态矩阵A和输出矩阵C,为此,可以通过如下的算法[7]实现:
(1)确定系统的Hankel矩阵;桥梁环境振动测试时,由于测点较多而测试设备数量有限,一般会选取固定的几个参考点,进行连续的测试;而其它测点测试时间相对较短。
H∈R(m+n)i×j,m表示参考点的个数,n表示其他测点的个数。式中上标r表示参考点的输出,下标p,f分别表示“过去”和“将来”,是Hankel矩阵划分块行的一种方式。
(2)计算“将来”输入行空间在“过去”输入行空间上的投影,并通过QR分解,在保持系统原有信息的情况下缩减数据。
式中,(·)+表示·的广义逆。通过对H进行QR分解,Pi可以表示为:
投影的计算是随机子空间算法的核心,它表示可以利用“过去”行空间的信息预测“将来”;
(3)对投影进行奇异值分解,并结合卡尔曼滤波理论计算系统状态矩阵A和输出矩阵C;
(4)对A进行特征值分解,得到特征值和特征向量,求解系统的模态参数。
EMD能有效的处理非平稳数据,并且是自适应的。从提出至今,EMD方法已经在很多领域得到广泛的运用[8-10]。由于没有精确的数学模型,分解过程是一种以经验为基底的,EMD也存在一些问题限制了其进一步的推广[11]。模态混叠就是其中一种。模态混叠是指原始信号存在断续或畸变的情况下,不同频率分量被分解到同一阶IMF中或同一频率分量被分解到不同的IMF中。文献[12] 通过对信号进行傅里叶谱分析,在EMD过程中设置间断频率,使每一个IMF只包含频率大于间断频率的成分,这种方法简单易行,但FFT本身不适合分析非平稳信号;文献[13] 引入屏蔽信号消除了仿真信号的模态混叠;文献[14] 利用EMD的二进滤波特性,结合屏蔽信号的使用,提出一种限制带宽的EMD方法(Bandwidth Restricted EMD,BREMD)。大型桥梁结构尺寸大,模态通常较为密集,且测试过程受到噪声干扰严重,如果直接应用EMD,一般会产生模态混叠现象。为此,利用限制带宽的EMD消除模态混叠。为了使BREMD能更精确的提取实测信号的IMF,对其作部分改进,改进后算法如下:
(1)对原始信号x(t)做标准EMD,使其只分解出一阶IMF1以及一个余量r1(t):
(2)对IMF1进行Hilbert变换,计算待加入的屏蔽信号的频率,根据EMD的二进滤波特性,取带宽系数为 1.42:
式中,a1(i)、f1(i)分别表示IMF1的瞬时幅值和瞬时频率。
(3)构造屏蔽信号
(4)对信号y(t)=x(t)+s(t)进行标准的EMD分解,得到一阶IMFs1和一个余量;
(5)判断h1(t)=IMFs1-s(t)是否存在模态混叠,如不存在,则h1(t)即为原始信号的第一阶IMF;如果仍然存在模态混叠,则把h1(t)当作x(t)重复式(1)~式(4)步骤n次,直至 h1n(t)不存在模态混叠,记作c1(t);判断是否存在模态混叠的方法可以通过对h1(t)作Hilbert变换,查看其瞬时频率在时间轴上的分布;
(6)从原始信号中减去第一阶IMF,得到x1(t)=x(t)-c1(t),重复上述过程,直至所有的IMF都分解出来。
为了验证限制带宽EMD在消除模态混叠中的效果,以一仿真信号作为算例。仿真信号中含有三种频率成分,分别为 10 Hz、15 Hz和 18 Hz,y=sin(20πt)+sin(30πt)+sin(36πt)采样频率取为200 Hz,采样时间为0 s~2.5 s。利用标准EMD和限制带宽EMD分解的结果如图1所示,图中c1~c3表示标准EMD的前三阶IMF,BR-c1~BR-c3表示 BREMD的前三阶IMF,可以看出,c1和c3中存在较为严重的模态混叠,BR-c1~BR-c3基本能代表仿真信号中的三种频率成分,因此BREMD显著的改善了信号的模态混叠。
图1 仿真信号的EMD和BREMDFig.1 The results of EMD and BREMD of simulation signal
如何确定系统阶次是随机子空间识别的一个难点,通常的做法是对系统阶次分别取一个最小值nmin和最大值nmax,循环求解位于nmin和nmax之间所有阶次系统的模态参数,并将计算的结果绘在一个二维图中,以频率为横坐标,系统阶次为纵坐标。相邻频率点和阻尼点小于事先设定的容差的,合并为一点,形成系统的稳定图,在每一阶频率处会有较多的稳定点构成稳定轴。频率和阻尼稳定点的判据如下:
式中,fi、ξi分别表示识别的第 i阶频率和阻尼比;式(12)将阻尼比的判定相对容差放松至30%,这是由于目前对阻尼比的认知水平有限,计算中采取的假定与实际情况也不太一致;实桥识别的阻尼比通常偏差较大,因此适当放松阻尼比的判定标准[12]。受到测试过程中的噪声影响,稳定图中会出现虚假模态,如何剔除这些虚假模态是稳定图方法应用的关键。分析知虚假模态是由于环境激励不满足白噪声假设,使得原始信号中含有噪声模态,加之系统模态密集,在稳定图中难以区分系统真实模态的稳定轴和噪声模态的稳定轴。如果能将原始信号分解成只含有结构某一阶固有振动信息的IMF,再利用稳定图方法,此时稳定图中主要包含这一阶振动的稳定轴,噪声模态则离散成一些跳点,不会影响系统模态参数的识别。因此,相比利用原始信号作为SSI的输入,利用IMF作为SSI的系统输入会得到更为清晰的稳定图,后面将结合实例验证这一判断。
通过上面的分析知改进后的模态参数识别的关键问题是解决EMD的模态混叠。为此,首先对各个测点原始信号进行低通滤波,以消除信号中的高频噪声;对去噪后的信号进行BREMD,消除EMD过程中的模态混叠,使得每一阶IMF只包含结构一阶振动信息;然后组装各个测点中含有相同频率成分的IMF,利用随机子空间方法,作系统的稳定图,通过拾取稳定轴,识别结构的模态参数。具体的流程见图2。
图2 改进后的模态参数识别流程Fig.2 Process of improved method for modal parameters identification
为了验证本文方法在实桥模态参数识别中的有效性,结合赣龙铁路上某特大桥环境振动实验,识别了该桥的模态参数。
赣龙铁路某特大桥主桥为预应力混凝土连续梁桥,跨径布置为60 m+2×100+60 m,5号墩墩高100 m,属于典型的高墩大跨铁路桥梁。为了了解该类桥型的振动特性,对其进行了环境激励实验。实验共布置49个测点,60 m跨的测点布置在跨度的八分点处,100 m跨测点布置在十六分点处。测试了各个测点的横向和竖向位移,采样频率为20 Hz,采样时间为5 min。桥跨及测点的布置图见图3。
高墩大跨桥梁的横向振动效应明显,因此选取横向位移作为原始信号,识别桥梁前四阶横向振动模态。为了验证识别的效果,在模态参数识别时采用了三种方法:基于FFT谱分析的峰值拾取法(Peak Peaking,PP法)、随机子空间识别方法以及文中叙述的方法,并将三种分析方法的结果与有限元计算结果进行了对比。
图3 桥跨及测点布置图(cm)Fig.3 The configuration of bridge span and measuring points(cm)
图4 为测点5处的横向位移的傅里叶频谱图,图中幅值比为FFT变换后的实部与虚部表示的复数的模长与序列长度的比值,它代表各阶频率所占能量的相对比值;图5为单纯利用随机子空间识别得到的稳定图,二者均能识别桥梁横向前四阶振动模态,但FFT是一种总体平均的概念,在分析非平稳信号时缺乏物理意义,且模态较为密集;而SSI识别的稳定图中的虚假模态已经严重干扰到真实模态的拾取,因此利用FFT频谱图和SSI稳定图识别模态参数时均存在较多的人工干预。对于模态密集的大型桥梁,这种人工干预往往并不容易。另一方面,由于关心的频段在0~3Hz之间,对原始信号进行截止频率为3 Hz的低通滤波。并对每个测点的信号都作限制带宽的EMD,图6显示了测点5处的横向位移信号的BREMD结果,从图中可以看出,通过引入限制带宽的屏蔽信号,各阶IMF的模态混叠得到了很好的抑制。组装相同参考点的同频率成分的IMF,再利用随机子空间识别其模态参数,并作出如图7所示的稳定图。为了方便对比,将傅里叶谱以一定的比例放入稳定图中,可以看出稳定图中只有一个清晰的稳定轴,而噪声模态离散成一些跳点,为模态参数的提取提供了便利。
图4 测点5横向位移的傅里叶谱Fig.4 FFT spectrum of lateral displacement of measuring point 5
图5 随机子空间识别的稳定图Fig.5 Stabilization diagram of SSI
该桥横向前四阶模态识别结果如表1所示,表中列出了有限元计算结果、PP法识别结果、SSI识别结果以及本文方法识别的结果。理论振型特征是通过各个自由度上振型参与向量值所判定的。相比识别结果,有限元计算结果偏大,一方面是有限元计算中对约束和支撑刚度的模拟比实际中的大所导致,另一方面是由于未将实测结果对理论计算模型进行反馈修正,使得理论模型不能完全反映真实情况;而PP法、SSI与本文方法识别结果吻合良好。相比PP和SSI识别结果,利用本文方法作出的稳定图稳定轴清晰,在识别模态参数时减少了人工干预,更适合自振频率低、模态较为密集的大型桥梁的模态参数识别。
通过分析EMD和SSI稳定图在模态参数识别中存在的问题,建立一种了大型桥梁模态参数识别的方法,并成功应用于实桥环境振动实验模态参数识别中。该方法在EMD中引入了限制带宽的屏蔽信号,使得EMD的模态混叠得到显著抑制,每一阶IMF代表结构的某一阶固有振动,使得稳定图中只有一个较为清晰的稳定轴,能够更准确的识别结构模态参数。与PP法和SSI识别结果的对比知,该方法识别的结果具有一定精度,且不存在模态密集导致模态参数识别困难的问题,具有一定的应用前景。
图6 测点5横向位移的BREMD结果Fig.6 BREMD of lateral displacement at measuring point 5
表1 模态参数识别结果Tab.1 Results of modal parameters identification
图7 BREMD后随机子空间识别的稳定图Fig.7 Stabilization diagrams of SSI after BREMD
[1] 静 行,袁海庆,赵 毅.基于独立分量分析的结构模态参数识别[J] .振动与冲击,2010,29(3):137-141.
[2] 陈 隽,徐幼麟.HHT在结构模态参数识别中的应用[J] .振动工程学报,2003,16(3):383 -388.
[3] 韩建平,李达文.基于Hilbert-Huang变换和自然激励技术的模态参数识别[J] .工程力学,2010,27(8):54 -59.
[4] 常 军,张启伟,孙利民.稳定图方法在随机子空间模态参数识别中的应用[J] .工程力学,2007,24(2):39 -44.
[5] 肖 祥,任伟新.实时工作模态参数数据驱动随机子空间识别[J] .振动与冲击,2009,28(8):148 -153.
[6] Huang N E,Shen Z,Long S R,et al.The empirical mode decomposition and the Hilbert spectrum for non-linear and non-stationary time series analysis[J] .Proc.R.Soc.Lond.A,1998;454:903 -995.
[7] Van O P,Moor D B.Subspace identification for linear systems:theory implementation applications[M] .Dordrecht,The Netherlands:Kluwer Academic Publishers,1996.
[8] 傅光耀,段晨东,刘健新.经验模态分解剔除结构地震响应的位移飘移[J] .振动工程学报,2009,22(4):438-442.
[9] 罗洁思,于德介,彭富强.基于EMD的多尺度形态解调方法及其在机械故障诊断中的应用[J] .振动与冲击,2009,28(11):84-86.
[10] Yang Y B,Chang K C.Extraction of bridge frequencies from the dynamic responses of a passing vehicle enhanced by the EMD technique[J] .Journal of Sound and Vibration,2009,322(4-5):718-739.
[11] 刘小峰,秦树人,柏 林.EMD中存在的问题及解决方法探讨[C] .第九届全国振动理论及应用学术会议论文集,2007,194 -200.
[12] 禹丹江,任伟新.信号经验模式分解与间断频率[J] .福州大学学报(自然科学版),2005,33(5):638 -642.
[13] Deering R,Kaiser J F.The use of a masking signal to improve empirical mode decomposition[C] .USA:IEEE International Conference,2005,18-23.
[14] 何启源.基于现代时频分析的环境激励模态参数识别方法研究[D] .重庆:重庆大学,2009.