小波多分辨率分析时变系统参数识别算法的鲁棒性研究

赵丽洁, 杜永峰,2, 李万润,2, 朱前坤,2

(1.兰州理工大学防震减灾研究所,甘肃 兰州 730050;2.兰州理工大学西部土木工程防灾减灾教育部工程研究中心,甘肃 兰州 730050)



小波多分辨率分析时变系统参数识别算法的鲁棒性研究

赵丽洁1, 杜永峰1,2, 李万润1,2, 朱前坤1,2

(1.兰州理工大学防震减灾研究所,甘肃 兰州 730050;2.兰州理工大学西部土木工程防灾减灾教育部工程研究中心,甘肃 兰州 730050)

针对多自由度时变系统参数识别问题,基于Daubechies小波多分辨率展开的时变参数辨识方法分析影响参数识别鲁棒性的各个因素。通过数值分析针对突变、线性慢变以及谐波快变的时变参数进行识别,研究结果表明:当基函数dbN一定时,在预先确立的分解尺度范围内,识别精度随分解尺度的增加而增加;待识别参数的频率特性对分解尺度的选择有很大影响,快时变参数比慢时变参数对分解尺度更为敏感;基函数dbN并不是影响识别精度的主要因素;在分解尺度相同的情况下,可以通过提高采样频率增加快时变参数识别精度。

时变系统; 参数识别; 小波多分辨率分析; 最优分解尺度; 算法鲁棒性

0 引言

时变参数的结构动力学问题,无论是正问题还是反问题一直都是研究的重要科学问题,也是实际工程问题中必须解决的关键问题[1]。随着工程科学朝高速、高精度及智能化方向发展,工程结构参数时域变化特点对系统动态特性的影响逐渐增大[2]。实际工程结构在服役期限内受到工作荷载或极端荷载作用时,其损伤不可避免且不断累积,本质上属于时变和非线性的结构系统[3-4]。时变物理参数的识别属于力学反问题研究范畴,时变参数识别问题主要是借鉴了控制理论、系统工程以及信号处理等领域的研究成果,逐渐形成各类方法的识别体系。研究人员一般从以下三种思路去考虑:一类是基于在线或递推技术的研究方法,常采用最小二乘类估计算法以及各种卡尔曼类估计算法,在此基础上引入常遗忘因子、变遗忘因子以及各种自适应因子矩阵来增强识别时变参数的时域跟踪能力[5-6];一类是子空间类方法,先后发展了基于集总数据的子空间方法和递推子空间方法,基本思想是通过跟踪系统矩阵的特征子空间来识别系统的时变参数[7];一类是基于信号的时频分析技术,从信号处理技术角度考虑,主要有HHT和小波变换等方法[8-11]。由于时变结构系统响应信号表现出非稳态特性,HHT和小波变换在处理非线性、非平稳信号方面具有独特优势。小波变换通过伸缩平移能够分析任意频率范围内的信号,具有良好的捕捉信号全局以及局部特性的能力。本文主要针对小波识别方法,即将每一个待识别参数在小波尺度函数上和小波函数上进行多分辨率展开,将时变系数转化成采用一系列小波系数表达的时不变系数的辨识问题[12]。

小波函数相比谐波、多项式等基函数更适于描述某些具有快变特征的参数。根据以上研究成果分析可知,基于小波理论对时变参数辨识逐渐形成基本理论体系。识别算法的鲁棒性是将其应用于实际结构进行系统识别的关键[13],但针对小波识别方法中影响参数识别鲁棒性的各个因素,如分解尺度的选择、小波基函数选择、抗噪声干扰能力,以及如何提高辨识精度等关键问题却鲜有报道。本文在已有研究的基础上,在Matlab环境下通过数值分析,针对突变、线性慢变以及谐波快变参数进行识别,为该识别方法在参数选取上给出一个定量、定性的评判规则,以期为应用于工程实际做出一些理论指导。

1 理论背景

1.1 Daubechies小波

Daubechies小波是法国学者Ingrid Daubechies提出构造的一系列二进制小波函数的总称,简写成dbN,N是小波的阶数。小波函数φ(t)和相应的尺度函数φ(t)的支撑区为2N-1,φ(t)的消失矩为N。dbN时域上是有限支撑的,即φ(t)长度有限,其高阶原点矩∫tpφ(t)dt=0,p=0～N;N值愈大,φ(t)的长度就越长;在频域上Ψ(ω),在ω=0处有N阶零点;φ(t)和它的整数位移正交归一,∫φ(t)φ(t-k)dt=δk;dbN小波具有较好的正则性,其中消失矩越高光滑性就越好,频域的局部化能力就越强,频带的划分效果越好。正是由于具有以上特点,Daubechies小波被Mallat应用到多分辨率分析框架体系中。dbN小波不具有对称性,且没有明确的解析表达式。图1为N=3、4时的db3小波函数与db4小波函数及对应的尺度函数曲线

1.2 小波多分辨率分析(WMRA)

基于Mallat多分辨率分析思想,把平方可积的函数f(t)∈L2(R)看成是某一逐级逼近的极限情况,也就是用不同分辨率基函数逐级逼近待分析函数f(t)[14]。

…,V0=V1⊕W1, V1=V2⊕W2, …, Vj=

Vj+1⊕Wj+1,…,

式中:Vj为尺度空间;Wj为小波空间;j∈Z,且j是从-∞到+∞的整数,j值愈大空间愈大。

对于一个能量有限且随时间变化的时间序列f(t),通过小波多分辨率可以近似表示为[14]:

(n=0,1,…,Nt-1)

(1)

图1 不同小波基函数Fig.1 Different wavelet basis function

其中:φj0,k(2j0t-k)和φj,k(2jt-k)分别是尺度函数φ(t)和母小波函数φ(t)的平移和伸缩的函数簇;cj0,k和dj,k分别是在尺度j0和j上的展开系数。式(1)中第一项给出了f(t)的低频分量或近似表达;第二项给出了对于不同分辨率尺度j的高频分量或细节表达。在式(1)中,确定对于不同分解尺度上展开系数cj0,k、dj,k取值个数,首先确定k0,Kj的取值范围,其与点数Nt、尺度函数支撑长度有关。消失矩为N尺度函数φ(t),φj,k(t)=2j/2φ(2jt-k)的支撑长度为[2-jk,2-j(k+2N-1)]。为保证平移伸缩之后φi,k(t)能够完全覆盖整个信号的长度,应保证最后一个基函数的初始时间大于信号的终止时间,且第一个基函数的终止时间小于信号的开始时刻。根据这个原则,k的取值范围应为k0=2-2N,Kj0=0,Kj=int(2jNt)-1。式(1)的离散形式矩阵表达式为:

F≈W(J)ξ(J)

(2)

2 识别方法

考虑一个多自由度体系的时变运动方程:

(3)

(4)式中:θ(t)为结构时变参数向量表达形式;恢复力向量R表达形式为R=[r1-r2,r2-r3,…,ri-ri+1,…,rn-1-rn]

(5)式中:ri(t)为i-1与i层之间的的恢复力,具体可以表达为:

(6)

R(t)=τ(t)θ(t)

(7)

由式(4)、(7)推知得:

(8)

y(t)=R(t)+e(t)=τ(t)θ(t)+e(t)

(9)

e(t)为测量噪声,写成矩阵的形式:

Y=ΓΘ+ε

(10)

式(10)可以看作一元线性回归模型。其中,测量输出向量Y=[Y1,Y2,…,Yn]T,而

Θ=[Θ1,Θ2,…,Θn]T;

Θi=[θi(0),θi(1),…,θi(Nt-1)]T

i从1到n,E为响应的测量噪声向量。根据多分辨率建模的思想,将时变参数矩阵Θ采用小波多分辨率展开写成矩阵形式:

Θ≈WΞ

(11)

W=diag[w1,w2,…,wn]是相应的小波展开的重构矩阵。ξi,wi分别为第i层时变参数的小波展开系数向量和小波基函数重构矩阵。将式(11)代入式(10)可以得到

Y=QΞ+ε

(12)

其中:Q=ΓW包含结构响应Γ和在不同分辨率下小波函数的采样值W。由式(12)可知,小波系数展开Ξ可以通过线性最小二乘伪逆解得到。

Ξ≈Q+Y

(13)

其中:Q+表示Q的伪逆,将Ξ回代到式(11)得到结构的时变物理参数。因此,采用多分辨率分析可以将结构的时变参数识别问题转化为一元回归模型中的时不变小波系数ξc、ξd的辨识问题。

3 数值分析

如图2所示,一个两自由度弹簧-质量-阻尼体系,M1=1 kg,M2=1 kg,k1=4 000 N/m,k2=4 000 N/m,c1=0.8 N·s/m,c2=0.8 N·s/m,假设振动过程中质量保持不变,分别考虑各参数线性慢变、突变工况来模拟参数变化,第一层结构的刚度参数k1在3 ～8 s退化25%;相对应的阻尼参数c1分别在3～8 s增加25%;k2、c2在10 s时相应退化了12.5%、25%。采样频率fs=50 Hz,点数Nt=1 001,以高斯白噪声作为输入激励,采用四阶Runge-Kutta法对运动方程进行迭代响应求解。为讨论不同因素对参数识别精度的影响,定义参数识别的平均绝对误差:

(14)

图2 两自由度弹簧-质量-阻尼体系模型Fig.2 A 2DOF spring-mass-damping model

3.1 分解尺度的选取

根据识别原理,对待时变参数进行小波展开之前,首先确定分解层数的选取范围[Pmin,Pmax]。长度为Nt的信号,由小波分解理论可知,P=log2Nt为信号的完全分解层数。由式(13)使方程有最小二乘伪逆解的最大的小波重构层数Pmin为

Pmin=ceil(1+lnq)

(15)

式中:q为每层识别参数的个数。对于Pmax的选择,参考文献[15]有:

(16)

图3为不同分解尺度下刚度、阻尼识别值与理论值的识别效果图,在预先确立的分解尺度范围内都能够准确地跟踪参数的变化趋势。由图4可知:(1)不管是阻尼时变参数还是刚度时变参数,相对误差均随分解尺度的增加而降低,识别精度越来越高。当分解尺度增加到一定程度时,反而MAPE增加,对于阻尼参数来说小波系数太多可能导致信号的过度拟合现象,阻尼识别时会出现震荡。(2)图4(a)中刚度的识别精度要比图4(b)阻尼的识别精度高,由于阻尼和刚度相差几个数量级,同时识别阻尼的精度稍差一些。(3)线性慢时变参数工况比突变工况时变精度高,突变参数相对于线性慢变参数对分解尺度的选择更加敏感。慢时变参数对分解尺度的选择范围更广,如图4(a)、(b)中所示J=-3～-7均能满足精度要求。通过大量的试算发现,待时变参数的频率特性对分解尺度的选择有很大的影响。

图3 不同分解尺度时变参数识别值与理论比较Fig.3 Comparison between identified values and theoretical values of time-varying parameters in different decomposition scales

图4 不同分解尺度下识别结果MAPE值Fig.4 MAPE of identification results in different decomposition scales

3.2 验证对噪声的敏感性

考虑到实际测量过程中输出响应不可避免地存在噪声影响,为验证该方法对噪声的抗干扰能力,以第一层的时变参数k1(t),c1(t)为研究对象,分别向输出响应数据中添加不同信噪比的高斯白噪声模拟实测响应。高斯白噪声的特点是其频率成分均匀分布在每一个频段,因此每个展开小波系数中不仅携带了有用频率信息,而且携带了噪声频率信息。为了消除噪声的干扰,在保证时变参数的频段范围保留的情况下,尽可能用少的系数去表达重构时变参数,因此采用J=-5分解尺度来讨论其抗噪性能。

图5 不同信噪比的识别误差Fig.5 Error of identification results with different SNR

图5所示，随着SNR的减少识别误差越来越大。相对于阻尼来说,刚度的识别抗干扰能力更强一些,阻尼本身识别能力就稍差一些,再加之噪声的影响,掩盖了真实的阻尼值。如表1所列,当SNR=50dB时,阻尼的识别误差只能控制在10%以内;刚度的识别精度非常高,当SNR=30dB,误差仍然在1%以内。由此说明,在未经过任何去噪处理的情况下,采用该方法时时变刚度参数的抗噪性能很好,而时变阻尼参数对噪声很敏感,识别效果较差。这是由于在一般的实际工程中,刚度远大于阻尼,刚度变化时的高频能量泄露及噪声导致的刚度引起的恢复力的变化比阻尼突变时高频能量泄露及噪声引起的阻尼力变化要大,因此时变刚度识别的相对误差要比时变阻尼的相对误差要小的多。

表1 不同信噪比的识别MAPE值

3.3 dbN基函数的选取

1.1节中讨论的Daubechies系列小波基具有紧支撑、正交归一以及较好的正则性等优点,不同的N值对应不同dbN系列的小波基函数。以下讨论阶数基函数的选取、分解尺度与阶数N对识别精度的关系。

从图6中可以看出,当固定某一分解尺度时,不同阶数的小波识别精度有一定的影响,但影响不是很大。J=-6～-3的大部分的尺度范围内,阻尼识别相对误差值几乎都集中在0～0.005[图6(a)],识别精度相差不大,而刚度的识别误差更是如此,如图(b)所示;分解尺度与阶数N的存在一定联系,为了保证识别精度,N的阶数可以选择的较小一些或当分解尺度较小时,可以采用较高的阶数N。另外,小波基函数的选择并不是影响识别精度的主要因素,不同dbN仍然遵循随着分解尺度的增加识别精度越来越高的整体趋势。

3.4 采样频率的影响

图6 识别误差与不同分解尺度、不同dbN之间的关系Fig.6 Error of identification results with different J and dbN

由1.2节小波多分辨率分析理论和分解尺度的选取原则式(16)、(17)可知,分解层数p确定时(即分解尺度J确定),小于p层的小波系数将被忽略,剩下部分系数重构得到信号的最高频率为:

(17)其中:fmax为原时变参数的最高频率。在分解层数p确定时,如果时变参数的频率成分相对较高,可以通过提高采样频率来弥补小于p层被忽略的小波系数所占的高频成分,以提高时变参数的识别精度。如算例1所示,第二层的k2(t)为:

其余参数同算例1相同,假设都采用相同的分解尺度J=-6,采样频率fs分别为25 Hz、50 Hz时进行识别。

由图7可知,当采样频率fs=25 Hz,采用J=-6进行分解时,根据多分辨率分析理论,包含0.39～12.5的频段范围内的信息的小波系数被忽略,k2(t)中谐波震荡的0.5 Hz高频成分并没有识别出来;当fs=50 Hz时,包含0.78～12.5频率成分小波系数被忽略,去掉该部分系数并不影响参数的识别精度,0.5 Hz的频率成分很准确地识别出来,与理论值相接近。

图7 不同采样频率识别k2(t)的结果Fig.7 Identification results of k2(t) at different fs

4 结论

本文通过数值分析针对突变、线性慢变以及谐波快变的时变参数进行识别,分析了基于Daubechies小波多分辨率展开的时变参数辨识方法,以及影响参数识别鲁棒性的各个因素。

(1) 在预先确立的分解尺度范围内都能够准确地跟踪参数的变化趋势,识别精度随着分解尺度的增加而增加;刚度参数的识别精度要比阻尼参数的识别精度高。

(2) 待时变参数的频率特性对分解尺度的选择有很大影响,线性慢时变参数工况比突变工况时变精度高,突变参数相对于线性慢变参数对分解尺度的选择更加敏感。

(3) 在未经过任何去噪处理的情况下,采用该方法时时变刚度参数的抗噪性能很好,而时变阻尼参数对噪声很敏感,识别效果较差些。

(4) 对某一待识别时变参数,在相同的分解尺度条件下,情况允许的条件可以提高采样频率以提高参数识别精度。

References)

[1] 于开平,庞世伟,赵婕.时变线性/非线性结构参数识别及系统辨识方法研究进展[J].科学通报,2009,54:3147-3156.

YU Kai-ping,PANG Shi-wei,ZHAO Jie.Advances in Method of Time-varying Linear/Nonlinear Structural System Identification and Parameter Estimate[J].Chinese Science Bulletin,2009,54:3147-3156.(in Chinese)

[2] 陈恩伟,陆益民,刘正士,等.Taylor展开的线性时变系统参数辨识及误差分析[J].机械工程学报,2011,47(7):90-96.

CHEN En-wei,LU Yi-min,LIU Zheng-shi.Parameter Identification and Error Analysis of Linear Time Varying System Based on Taylor Expansion[J].Journalof Mechanical Engineering,2011,47(7):90-96.(in Chinese)

[3] 郑山锁,代旷宇,孙龙飞.钢框架结构的地震损伤研究[J].地震工程学报,2015,37(2):290-296.

ZHENG Shan-suo,DAI Kuang-yu,SUN Long-fei.Research on the Seismic Damage of Steel Freame Structure[J].China Earthquake Engineering Journal,2015,37 (2):290-296.(in Chinese)

[4] 裴强,王丽,全厚辉.钢筋混凝土框架结构参数时变特性的研究[J].地震工程与工程振动,2013,33(1):41-46.

PEI Qiang,WANG Li,QUAN Hou-hui.Study on Characteristics of Time-varying Parameters of Reinforced Concrete Frame Structure[J].Earthquake Engineering and Engineering Vibration,2013,33 (1):41-46.(in Chinese)

[5] Cooper J E,Worden K.On-line Physical Parameter Estimation with Adaptive Forgetting Factors[J].Mechanical Systems and Signal Processing,2000,14:705-730.

[6] Yang J N,Lin S.Identification of Parametric Variations of Structures Based on Least Square Estimation and Adaptive Tracking Technique[J].Journal of Engineering Mechanics ASCE,2005,131(3):290-298.(in Chinese)

[7] 庞世伟,于开平,邹经湘.用于时变系统辨识的递推子空间方法[J].振动工程学报,2005,18:233-237.

PANG Shi-wei,YU Kai-ping,ZOU Jing-xiang.Time-varying System Identification Using Recursive Subspace Method Based on Free Response Data[J].Journal of Vibration Engineering,2005,18:233-237.(in Chinese)

[8] 黄天立,邱发强,楼梦麟.基于改进HHT方法的密集模态结构参数识别[J].中南大学学报:自然科学版,2011,42(7):2055-2062.HUANG Tian-li,QIU Fa-qiang,LOU Meng-lin.Application of an Improved HHT Method for Modal Parameters Identification of Structures with Closely Spaced Modes[J].Journal of Central South University:Science and Technology,2011,42(7):2055-2062.(in Chinese)

[9] 许鑫,史治宇,WieslawJstaszewski.利用加速度响应连续小波变换的时变系统物理参数识别[J].振动工程学报,2013,26(1):8-13.

XU Xin,SHI Zhi-yu,WieslawJstaszewski.Time-varying System Physical Parameters Identification Using the Continuous Wavelet Transform of Acceleration Response[J].Journal of Vibration Engineering,2013,26(1):8-13.(in Chinese)

[10] 黄东梅,周实,任伟新.基于小波变换的时变及典型非线性振动系统识别[J].振动与冲击,2014,33(13):124-129.

HUANG Dong-mei,ZHOU Shi,REN Wei-xin.Parameter Identification of Time-varying and Typical Nonlinear Vibration System Based on Wavelet Transform[J].Journal of Vibration and Shok,2014,33(13):124-129.(in Chinese)

[11] Dziedziech K,Staszewski W J,Uhl T.Wavelet-based Modal Analysis for Time-variant Systems[J].Mechanical Systems and Signal Processing,2015,50-51:323-337.

[12] Chang C C,Shi Y.Sub-structural Time-Varying Parameter Identification Using Wavelet Multiresolution Approximation[J].Journal of Engineering Mechanics,2012,138:50-59.

[13] 朱旭东,吕西林.多自由度非线性结构参数识别的鲁棒性研究[J].中南大学学报:自然科学版,2013,44(1):303-308.

ZHU Xu-dong,LV Xi-lin.Robust Study on Parametric Identification of Multi-degree-of-freedom Nonlinear Structure[J].Journal of Central South University:Science and Technology,2013,44(1):303-308.(in Chinese)

[14] Mallat S.A Wavelet Tour of Signal Processing[M].London:Academic,1999.

[15] Wei H L,Billings S A.Identification of Time-varying Systems Using Multiresolution Wavelet Models[J].International Journal of Systems Science,2002,33(15):1217-1228.

[16] 杜永峰,赵丽洁,党星海.基于最优复Morlet小波的结构密集模态参数识别[J].振动与冲击,2015,34(5):136-139.

DU Yong-feng,ZHAO Li-jie,DANG Xing-hai.Identiication of Structural Closely Modes Parameters Based on Optimized Complex Morletwavelet[J].Journal of Vibration and Shok,2015,34(5):136-139.(in Chinese)

Robustness of a Parametric Identification Algorithm for Time-varying Systems Based on Wavelet Multi-resolution Analysis

ZHAO Li-jie1, DU Yong-feng1,2, LI Wan-run1,2, ZHU Qian-kun1,2

(1.InstituteofEarthquakeProtectionandDisasterMitigation,LanzhouUniversityofTechnology,Lanzhou730050,Gansu,China; 2.WesternCenterofDisasterMitigationinCivilEngineeringofMinistryofEducation,LanzhouUniversityofTechnology,Lanzhou73005,Gansu,China)

Research,based on Daubechies wavelet multi-resolution analysis,was carried out to solve parameter identification problems in multiple degrees of freedom time-varying systems.In order to improve identification efficiency and accuracy,numerical experiments,based on the above method,were conducted to study the various factors that affect performance.The results show that when the basic function dbNwas fixed in the preset decomposition scale,identification accuracy increased with an increase in the decomposition scale.The frequency component of the time-varying parameters had great influence on the choice of decomposition scale,and the fast time-varying parameters were more sensitive than the slow.The choice of the basic function dbNaffects the identification accuracy,but is not a key factor; an increase in sampling rate can improve the identification accuracy of fast time-varying parameters under the same decomposition scale.

time-varying system; parameter identification; wavelet multi-resolution analysis; optimal decomposition scale; robustness

2015-06-22

国家自然科学基金项目(51178211,51578274);甘肃省青年科技基金计划(148RJYA004)

赵丽洁(1988-),女,博士生,主要从事结构健康监测研究。E-mail:ljzhaocz@126.com。

杜永峰(1962-),男,博士,教授,博导,主要从事结构减震控制、结构健康监测研究。

TB12; TU352.1+2

A

1000-0844(2016)05-0720-08

10.3969/j.issn.1000-0844.2016.05.0720