摘要
扩展卡尔曼滤波(Extended Kalman filter, EKF)理论由线性卡尔曼滤波理论发展而来。在结构的物理参数未知的情况下,扩展卡尔曼滤波理论将结构的物理参数与其状态向量组成增广状态向量,对增广状态向量进行识别,从而得到修正的结构参数。本文引入扩展卡尔曼滤波理论,对结构参数进行识别。为了验证该方法的有效性,引入1个时变参数三自由度系统作为仿真算例,分别探讨扩展卡尔曼滤波中的各项滤波参数(模型噪声的协方差矩阵Q,测量噪声的协方差矩阵R)对结构的时变参数追踪性能的影响。结果表明,合适的滤波参数能使得算法更迅速,准确地识别出结构中的时变参数变化趋势。还通过一个悬臂梁仿真算例,证明了即使在部分测点没有测量的情况下,EKF算法也能识别出结构参数。
参数识别是动力学的第一类逆问题。参数识别问题按识别参数的种类可以分为模态参数识别与物理参数识别;而如果按识别参数的方法,可以分为时域法参数识别和频域
这些文献都对扩展卡尔曼滤波方法进行了改进,但并没有提及滤波参数的设置对识别效果的影响,尤其是结构中参数会随时间变化的情况。事实上,随着算法的迭代计算,EKF算法对测点响应的利用程度会越来越小,此时识别结果无法追踪到时变参数的变化情况。因此,必须选择恰当的滤波参数(模型噪声的协方差矩阵Q,测量噪声的协方差矩阵R)。本文首先介绍了一种结构参数识别方法,适用于外界载荷已知的情况,并试图通过仿真算例说明如何选择最优的滤波参数,来使得EKF算法对时变结构参数的跟踪性能达到最佳,即算法能快速、准确地识别出结构参数的变化情况。
本文使用扩展卡尔曼滤波方法作为结构参数的识别方法。对于1个n自由度的动力学系统,当结构中的未知参数,例如质量矩阵M,刚度矩阵K或阻尼矩阵C中的某些参数待确定时,参考卡尔曼滤波理论,将这些未知参数与状态向量x一起组成扩展状态向量z
(1) |
式中: 为待确定的未知参数向量,长度设为m;表示状态向量,其长度为2n。令,此时,卡尔曼滤波理论中的状态更新方程可以重新写为
(2) |
而测量更新方程的形式与观测量的类别有关,本文选择的测量类别是加速度信号。此时,测量更新方程可以写为
(3) |
式中: M表示大小为的质量矩阵; C表示大小为的阻尼矩阵; K表示大小为的的刚度矩阵; ,,分别代表结构的加速度向量,速度向量,以及位移向量,向量的长度为n;代表激励的影响矩阵,与激励在结构上的作用位置有关,由0和1组成,它的矩阵大小为;表示激励向量,长度为s; 为由0和1组成的位置矩阵。由于质量矩阵M,刚度矩阵K或阻尼矩阵C包含了未知参数,这也就意味着关于扩展状态向量的状态更新方程是一个非线性方程,同样测量更新方程也是一个非线性方程。
利用泰勒展开,将上述两个非线性方程(1, 2)转化为多项式之和的形式,并略去高阶项,仅保留一阶多项式,由此可将非线性方程近似看成线性方程
(4) |
(5) |
式中: 为在时刻的扩展状态向量;为在时刻的扩展状态向量; 为状态更新方程对扩展状态向量的一阶偏导; 为测量更新方程对扩展状态向量的一阶偏导,又称为雅各比矩阵。由泰勒展开定理可知,只要时刻的扩展状态向量与,之间大小差距很小,式(
得到离散化的近似线性等式(
(6) |
(7) |
由
(8) |
(9) |
式中:为卡尔曼滤波增益,它表示对测量响应的利用程度。通过方程(8,9)可以使得扩展状态向量中的参数发生变化。因此,可以认为扩展卡尔曼滤波算法是通过测量更新过程实现了对参数的追踪,只要加大对测量响应的利用程度,就能使较快地收敛到真实参数值。其求法与卡尔曼理论中的相同,可以写为
(10) |
而、分别是非线性表达式,对扩展状态向量求偏导而得到的矩阵,即雅可比矩阵。若扩展状态向量所包含的未知参数仅为刚度或阻尼,则,可以分别写为
(11) |
(12) |
由此,可以得到扩展卡尔曼滤波算法的流程如
1.给定初始条件 和 2.时间更新
3.测量更新
|
在前文状态更新方程的推导中,假定未知参数不会随时间变化,即。这也就意味着仅凭状态更新方程,只能识别出时不变的参数,对时变参数或缓变参数是无法识别出来的。对于变化参数的追踪,是通过测量更新方程来实现的。卡尔曼滤波增益表示了对测量更新方程的利用程度。若希望能实现对时变参数的追踪,则对测量更新方程的利用程度必须加大,总的宗旨是能尽量利用测量更新方程,即令在迭代过程中仍保持较大值。
由卡尔曼滤波增益的计算公式可知,卡尔曼增益由测量噪声方差R和共同决定,初始时刻偏大,对状态的估计主要依赖于测量更新方程。随着迭代过程的进行,测量对状态的修正作用不断减少,系统反而依赖于状态更新方程,此时会逐渐变小。
为了使得滤波过程中,一直保持对测量更新方程的利用,需要选择合适的滤波参数。比如:模型噪声的协方差矩阵Q,对应状态向量部分的取值较小。对应参数部分的取值较大,这相当于告知滤波系统对结构未知参数的估计存在较大偏差,促使算法识别的结构参数会向真实值靠拢;测量噪声的协方差矩阵R的取值偏小,这相当于承认滤波系统所测量的响应具有很高的置信度,加大系统对测量更新方程的利用程度,从而加强算法对结构时变参数的追踪能力;然后是外界的测量噪声对参数识别效果的影响,算法能在响应存在噪声干扰的情况下识别出结构未知参数,这是通过调整测量噪声协方差矩阵R的方式来实现的;最后是观测点数的影响,选择悬臂梁作为仿真算例,由于悬臂梁既有转角自由度,又有竖向位移自由度,选择仅仅观测竖向位移自由度可以看出算法仍能较好地识别出结构参数变化的情况。
本文选择一个三自由度系统,如

图1 三自由度系统
Fig.1 A three degree of freedom system
在初始状态下,k1 = k4 = 200 N / m, k2 = k3 = 100 N / m;m1 = m2 = m3 = 1 kg。并假设这个三自由度系统的阻尼是比例阻尼,且有,=0.05, =0.02。在1.5 s时该结构中的刚度k1从200 N / m突然下降到150 N / m。假设激励作用的位置和形式已知,质量块m1受到激励,由此可计算得到结构各自由度的响应。为了便于讨论,选取加速度响应作为测量向量,且不施加噪声到测量向量中去。使用本文提出的扩展卡尔曼滤波算法,此时激励已知(作用位置,幅值大小),而未知参数仅选择k1,则扩展状态向量是一个7×1列向量。
由于参数是未知的,所以初始参数的选择必然会有较大误差,而结构的初始位移和速度(即状态向量)却可以事先测得,所以可以设定成真实值。相应的初始估计误差协方差关于状态向量的部分可以设定为一个较小的值,初始估计误差协方差关于参数的部分必须设定为一个较大的值。令初始扩展状态向量,初始误差方差分别写为
讨论模型噪声协方差矩阵Q的取值对时变参数识别的影响。Q为一个7阶对角矩阵,前6阶对应的是状态更新方程对结构的位移和速度所产生的误差,第7阶代表的是状态更新方程对参数k1所产生的误差。两者是不一样的,因此必须分开讨论。令测量噪声的协方差矩阵R为对角元素为的三阶对角矩阵。首先考虑模型噪声的协方差矩阵Q参数部分取值的影响,假定Q关于状态向量部分的取值为,关于参数部分的取值为x,即
x分别取值为1,,,时,讨论其参数部分的取值对参数识别效果的影响。同样,为了讨论状态向量部分取值的影响,令Q关于参数部分的取值为,关于状态向量部分的取值为x。EKF算法识别参数的效果如


图2 模型噪声的协方差矩阵Q对参数识别效果的影响
Fig.2 Influence of covariance matrix Q of model noise on parameter identification effect
从
对于测量噪声的协方差矩阵R,其最优取值与测量向量中的噪声有关。现讨论对参数识别效果的影响。令模型噪声的协方差矩阵Q状态向量部分取值,参数部分取值。令测量噪声的协方差矩阵R为对角元素为x的三阶对角矩阵,当x取不同值时,参数的识别结果如

图3 测量噪声的协方差矩阵R对参数识别效果的影响
Fig.3 Influence of covariance matrix R of measurement noise on parameter identification effect
测量噪声协方差矩阵R取一个较小值时,此时卡尔曼滤波增益保持一个较大的值,算法的追踪参数性能较好。但值得注意的是,R的取值不能太小,它需要与观测向量中的噪声相适应,如果选择过小的值,导致这些估计的误差超过了实际误差,识别的参数必然不会与实际情况相同,严重时会发生发散。因此,需要为测量噪声协方差矩阵R取一个合适的值,确保算法能较快地识别出结构参数值,识别结果不会发散。
当响应中存在噪声干扰时,本文提出的算法仍能正确识别出参数变化情况。在此分别讨论测量响应中不含噪声与含有5%噪声的情况。令模型噪声的协方差矩阵Q状态向量部分取值,参数部分取值。由前文分析可知,测量噪声的协方差矩阵R的最优取值与噪声水平的大小有关。对于响应中不含噪声时,设定测量噪声的协方差矩阵R是对角元素为的三阶对角矩阵,可以看出参数识别效果很好;而对于响应含有5%噪声时,此时R的取值必须大于,与观测噪声水平相适应。若选取R是对角元素为或的三阶对角矩阵,可以看出算法能较好地识别出参数变化情况。
由


图4 测量噪声对参数识别效果的影响
Fig.4 Effect of measurement noise on parameter identification
为说明观测点数的影响,取一个悬臂梁作为仿真对象,梁长度 L=1 m,矩形横截面尺寸为宽 b=0.1 m,高 h=0.01 m,材料密度ρ = 2 770 kg/

图5 悬臂梁结构
Fig.5 Cantilever beam structure
由
为了便于讨论,把梁单元的刚度特性作为单元发生刚度变化的指标。假定梁单元①的刚度特性在第1.5~3.5 s发生线性下降,下降至原来值的70%。其他单元的参数并不发生变化。在响应中添加1%的白噪声,利用本文所提出的算法,把该悬臂梁前3个单元的刚度特性作为未知参数向量,即扩展状态向量为,依次识别出各个刚度系数。
选择只观测第2,3,4,5节点上的竖直加速度,不测量转角自由度,因此这也算是部分观测。设模型噪声的协方差矩阵G,以及测量噪声的协方差矩阵R分别为
观测全部测点与部分测点的区别仅需改变测量方程中的位置矩阵。若选择只记录结构的竖向加速度响应,因此可写为
同理可求出既观测转角自由度,又观测竖向自由度情况下的参数识别效果,由此可得参数识别结果如图

图6 刚度系数E1I1识别效果图
Fig.6 Identification effect of stiffness coefficient E1I1

图7 刚度系数E2I2识别效果图
Fig.7 Identification effect of stiffness coefficient E2I2

图8 刚度系数E3I3识别效果图
Fig.8 Identification effect of stiffness coefficient E3I3
由图
在载荷已知的情况下,扩展卡尔曼滤波算法对结构参数的识别主要是通过测量更新方程实现的,而卡尔曼滤波增益表示了算法对测量更新方程的利用程度。为了实现算法对时变参数的追踪,本文提出了通过设置合适的滤波参数来保持较大的。这其中包括:模型误差Q关于参数部分取较大值,模型误差Q关于状态向量的部分取较小的值。还讨论了测量噪声对参数识别效果的影响,测量噪声的协方差矩阵R必须与测量噪声相适应,并尽量选取一个较小值。仿真结果表明,这样的滤波参数设置使得EKF算法在运行过程中能保持对测量更新方程的利用,从而能追踪到结构参数的变化情况,识别结构的时变参数。最后,通过悬臂梁仿真算例,证明了在部分测点观测的情况下,算法也能较好地识别出结构的参数变化情况,体现了EKF算法的优越性。
参考文献
BERNAL D, GUNES B. Flexibility based approach for damage characterization: Benchmark application[J]. Journal of Engineering Mechanics, 2004, 130(1):61-70. [百度学术]
FRISWELL M I , PENNY J E T. The choice of orthogonal polynomials in the rational fraction polynomial method[J]. Modal Analysis, 1993, 8(3):257-262. [百度学术]
SATO T,QI K.Adaptive H1 filter:Its application to structural identification[J]. Journal of Engineering Mechanics (ASCE), 1998, 124(11):1233-1240. [百度学术]
YOSHIDA I. Damage detection using Monte Carlo filter based on non-Gaussian noise[C]∥ Proceedings of Structural Safety and Reliability, ICOSSA 2001. Lisse: Swet & Zeitinger,2001: 8. [百度学术]
YANG J N, HUANG H. Sequential non-linear least-square estimation for damage identification of structures[J]. International Journal of Non-Linear Mechanics, 2006,41(1): 124-140. [百度学术]
LOH C H, LIN C Y, HUANG C C. Time domain identification of frames under earthquake loadings[J].Journal of Engineering Mechanics, 2000, 126(4):693-703. [百度学术]
黄小平,王岩.卡尔曼滤波原理及应用MATLAB仿真[M].北京: 电子工业出版社, 2015 [百度学术]
Huang Xiaoping,Wang Yan. Principle of Kalman filter and its application in MATLAB simulation[M]. Beijing: Electronic Industry Press, 2015. [百度学术]
KALMAN R E. A new approach to linear filtering and prediction problems[J]. Journal of Fluids Engineering, 1960, 82(1): 35-45. [百度学术]
HOSHIYA M, SAITO E. Structural identification by extended Kalman filter[J]. Journal of Engineering Mechanics, 1984, 110(12): 1757-1770. [百度学术]
尚久锉.卡尔曼滤波法在结构动态参数估计中的应用[J].地震工程与工程震动,1991, 11(2): 62-72. [百度学术]
Shang Jiucuo. Application of Kalman filter in structural dynamic parameter estimation[J]. Earthquake engineering and Engineering Vibration, 1991, 11(2): 62-72. [百度学术]
李杰.基于微分算子变换的广义卡尔曼估计方法[J].计算结构力学及其应用,1995, 11(4): 394-400. [百度学术]
Li Jie. Generalized Kalman estimation method based on differential operator transformation[J]. Computational Structural Mechanics and Its Application,1995, 11(4): 394-400. [百度学术]
张浩.基于强跟踪无迹卡尔曼滤波的结构时变参数识别[D].兰州: 兰州理工大学, 2016. [百度学术]
Zhang Hao. Time varying structural parameter identification based on strong tracking unscented Kalman filter[D]. Lanzhou: Lanzhou University of Technology, 2016. [百度学术]
穆腾飞.基于自适应追踪技术的迟滞非线性系统参数识别与损伤检测研究[D].南京: 南京航空航天大学, 2016. [百度学术]
Mu Tengfei. Parameter identification and damage detection of hysteretic nonlinear systems based on adaptive tracking technique[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2016. [百度学术]
YANG J N, PAN S, HUANG H. An adaptive extended Kalman filter for structural damage identifications II:Unknown inputs[J]. Struct Control Health Monit, 2007,14(3):497-521. [百度学术]
HUANG Q, XU Y L, LIU H J. An efficient algorithm for simultaneous identification of time-varying structural parameters and unknown excitations of a building structure[J].Engineering Structures, 2015, 98:29-37. [百度学术]