基于时频分析的非线性结构参数识别方法对比分析

袁平平,张 健,吴慧山

(1.江苏科技大学 土木工程与建筑学院, 镇江 212005) (2.中设设计集团股份有限公司, 南京 210014)

土木工程结构所处的复杂自然环境和外部激励对结构的特征参数和力学行为有着显著的影响．结构在地震、台风等较强荷载作用下会表现出一定的非线性行为,因此准确地识别出非线性结构动力系统的瞬时特性,对结构的状态评定和安全评估具有重要的工程价值意义．作为非线性响应信号分析方法之一,时频分析法将一维的时域信号扩展到二维的时频平面上,它可以清晰地描述信号的能量随时间和频率的变化规律．目前,此类方法主要包括短时傅里叶变换法、Wigner-Ville分布法、希尔伯特-黄变换法(Hilbert-Huang transform,HHT)和小波变换法等[1]．

文献[2]利用短时傅里叶变换方法对高层框架结构的响应数据进行了分析,得到了结构的模态参数．结合最优窗自适应傅里叶变换方法,文献[3]对具有非平稳特征的悬索桥模型以及钢框架模型进行了参数识别．通过对HHT的研究,文献[4]开发了一套环境激励下结构模态参数识别的程序．结合带通滤波技术,文献[5]将HHT应用到桥梁的模态参数识别中．由于HHT存在分解能力不足以及各本征函数(intrinsic mode function,IMF)分量非正交等问题,文献[6]提出正交化经验模式分解(empirtia mode deconpostion，EMD)方法,并成功地将其应用于密集模态结构的参数识别．文献[7]通过离散小波变换将结构的时变参数在多尺度上展开,并利用最小二乘法来识别低频尺度的展开系数,从而得到了原始的时变参数．文献[8]采用罚函数来降低噪声影响,并引入动态规划来提高小波脊线的精度．基于小波变换,文献[9]对电力系统谐波监测进行研究．文献[10]联合小波变换和重组法提出同步挤压小波变换时频分析方法,通过重组小波变换的时频图来获得较高精度的时频曲线．为了解决同步挤压算法中小波系数的发散问题,文献[11-12]提出广义的同步挤压小波变换算法,并利用齿轮箱算例验证其可行性．文献[13]利用同步挤压小波变换成功地识别了时变结构和非线性结构的瞬时频率．

文中基于HHT、小波变换和同步挤压小波变换等时频分析方法，通过一个具有非线性刚度的框架结构数值算例对各方法的特点、适用范围、识别精度进行了对比分析,进一步为非线性结构的参数识别提供了参考依据．

1 基于时频分析的非线性结构参数识别

1.1 基于HHT的瞬时参数识别

对于任意一个单分量信号x(t),可以根据HHT得到其解析形式为:

z(t)=x(t)+jH[x(t)]=a(t)ejθ(t)

(1)

对公式(1)进行微分,得到x(t)的瞬时频率:

(2)

对于多分量信号,为了使通过HHT获得的瞬时频率能够真实地反映结构响应信号的物理机制,文献[14]提出IMF的概念和EMD方法．

假设X(t)为原始信号,则对其进行EMD的具体步骤如下:首先,获取X(t)中的局部极大值序列和局部极小值序列,然后采用三次样条插值曲线分别拟合局部极大值序列和局部极小值序列,得到原始信号的上下包络线,并计算其平均值m1(t)，将X(t)减去m1(t)得到分量h1(t)，即:

h1(t)=X(t)-m1(t)

(3)

此过程称为一次“筛选”．重复上述筛选过程,以获得最终的IMF分量．在第二次的筛选过程中,可以将h1(t)作为原始信号,同样地计算其三次样条上下包络线的平均值m2(t)，然后得到:

h2(t)=h1(t)-m2(t)

(4)

c1(t)=hk(t)

(5)

文献[14]给出了一种类似柯西收敛准则的“筛选停止准则”,即:

(6)

一般情况下,当SDk的取值在0.2～0.3,可结束筛选过程．

通过式(5)得到的第一阶IMF分量c1(t)是原始信号X(t)中的最高频的分量．将第一阶IMF分量c1(t)从原始信号X(t)中分离出来,即可得到相应的残余信号r1(t),即:

r1(t)=X(t)-c1(t)

(7)

将残余信号r1(t)作为新的信号,对r1(t)重复上述的“筛选”处理便可得到第二阶IMF分量c2(t),将r1(t)减去c2(t)后可以得到r2(t),如此重复循环下去,得到一系列的IMF分量cj(t)和残余信号rj(t):

rj(t)=r(j-1)(t)-cj(t)

(8)

当满足以下两个条件中的一个时,整个筛选过程结束,EMD过程完成:(1) 得到的IMF分量cn(t)或者残余信号rn(t)小于预先假定值;(2) 残余信号rn(t)变成时间的单调函数．经过EMD,原始信号X(t)可以表示为n阶IMF分量和余量rn(t)的和,即:

基体钢片：直径25 mm、厚度3 mm的普碳钢Q235圆形钢片.先用粒径为10μm，88μm，65μm碳化硅砂纸打磨Q235普碳钢圆片至其表面无明显划痕，然后用乙醇溶液清洗掉试样表面的铁屑及灰尘，干燥后将其置于密封袋内封存备用.Q235钢片主要化学元素C，Si，Mn，S和P的质量分数分别为0.135%，0.178%，0.409%，0.037%，0.017%.

(9)

结合上述EMD方法将非线性结构的响应信号分解成一系列IMF分量,然后利用分量信号的HHT便可得到其相应的瞬时幅值和瞬时频率变化曲线,从而识别出非线性结构系统的瞬时参数．

1.2 基于连续小波变换的瞬时参数识别

渐进单分量实信号x(t)，其解析信号表示为Zx(t)=Ax(t)exp[jφx(t)]．采用连续小波函数ψ(t)=Aψ(t)exp[jφψ(t)]对其进行小波变换,得到:

(10)

相应的小波系数的相位为:

(11)

由式(11)可知,其时间-尺度平面上小波脊点上的小波尺度为:

(12)

将平稳相位点ts(a,b)=b带入式(12),可得:

(13)

由式(13),其分母为渐进单分量信号x(t)的瞬时角频率ω(b)．将整个时间-尺度平面上的小波脊点(r(b),b)连成曲线,可得到小波脊线．由于小波函数的尺度与频率一一对应且成反比关系,因此小波脊线可以表征信号的瞬时频率,即:

(14)

提前给定小波函数,便可得到小波中心频率φψ′(0)．因此,根据公式(14)可知,只需求得小波脊线上的小波尺度r(b),便可以求得相应的瞬时频率．非线性结构的响应信号可以分解成一系列IMF．根据连续小波变换的线性叠加性质,多分量信号的小波变换系数可以表示为各分量的小波系数的叠加,即:

(15)

选择合理的小波母函数和小波参数对响应信号进行连续小波变换,可以得到多条互不干扰的小波脊线,从而将响应信号的各分量有效地分离出来,并提取出各分量信号的瞬时频率．

1.3 基于同步挤压小波变换的瞬时参数识别

根据Plancherel定理,连续小波变换的频域表达式为:

(16)

对于单一简谐波X(t)=Acos(ωt),其连续小波变换为:

(17)

(18)

(19)

若a和w为连续变量,则式(19)转换为:

(20)

(21)

如果Δi=1,则公式(19)中的同步挤压小波系数可表示成:

(22)

2 两层非线性框架模型数值算例

为了验证上述3种方法对非线性结构进行参数识别的可行性和准确性,采用如图1的两层非线性剪切框架结构模型进行数值模拟计算．框架结构系统的质量、层间线性刚度系数、非线性刚度系数以及阻尼系数如表1．

表1 两层剪切框架模型结构参数Table 1 Structural parameters of the two-layer shear frame model

图1 非线性框架结构系统Fig.1 Nonlinear frame structure system

其相应的运动微分方程为:

(23)

给定结构的初始条件x1=0.1 m，x2=-0.05 m,采用四阶Runge-Kutta法求解结构的自由振动响应．图2为非线性框架结构的第一层位移响应信号,其采样频率为100 Hz,采样时间为10 s．

对第一层位移响应进行EMD分解,得到的前两阶IMF分量如图3．

同时,分别对前两阶IMF分量进行HHT便可得到相应的瞬时频率．但由于EMD分解能力的不足,残余信号中仍包含低频分量,从而导致瞬时频率的识别精度受到不利影响．

采用复Morlet小波对位移响应进行连续小波变换,小波中心频率fc取1Hz,带宽fb为2,得到小波量如图4，图可以看出，存在两条小波脊线,其分别对应位移响应的两个振动成分．通过小波系数模的极大值提取小波脊线，可以识别出结构响应的瞬时频率．同时,在小波量图的基础上进行同步挤压，得到瞬时频率的变化曲线如图5．

图5 非线性框架结构的瞬时频率Fig.5 Instantaneous frequency of thenonlinear frame structure

3 结论

3种时频分析方法都能识别出非线性框架结构的瞬时频率,且计算复杂难易程度、计算效率相当．对于HHT方法,由于EMD分解的不足和HHT的端点效应,其识别的瞬时频率在两端偏差较大,且不适用于密集模态的识别;对于小波变换方法,提取的小波脊线存在许多毛刺,从而影响了瞬时频率的识别精度;而同步挤压小波变换对小波系数进行挤压重组,其计算较小波变换方法耗时略长,但识别的瞬时频率则大体上消除了毛刺的影响,能够较为准确地反映频率的瞬时变化情况．数值结果表明:同步挤压小波变换识别效果好于HHT和小波变换的识别结果．