网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

扩展卡尔曼滤波对时变参数追踪性能的影响研究  PDF

  • 唐宏志
  • 姜金辉
南京航空航天大学机械结构力学及控制国家重点实验室,南京 210016

中图分类号: TB124

最近更新:2023-02-22

DOI:10.16356/j.1005-2615.2022.02.017

  • 全文
  • 图表
  • 参考文献
  • 作者
  • 出版信息
EN
目录contents

摘要

扩展卡尔曼滤波(Extended Kalman filter, EKF)理论由线性卡尔曼滤波理论发展而来。在结构的物理参数未知的情况下,扩展卡尔曼滤波理论将结构的物理参数与其状态向量组成增广状态向量,对增广状态向量进行识别,从而得到修正的结构参数。本文引入扩展卡尔曼滤波理论,对结构参数进行识别。为了验证该方法的有效性,引入1个时变参数三自由度系统作为仿真算例,分别探讨扩展卡尔曼滤波中的各项滤波参数(模型噪声的协方差矩阵Q,测量噪声的协方差矩阵R)对结构的时变参数追踪性能的影响。结果表明,合适的滤波参数能使得算法更迅速,准确地识别出结构中的时变参数变化趋势。还通过一个悬臂梁仿真算例,证明了即使在部分测点没有测量的情况下,EKF算法也能识别出结构参数。

参数识别是动力学的第一类逆问题。参数识别问题按识别参数的种类可以分为模态参数识别与物理参数识别;而如果按识别参数的方法,可以分为时域法参数识别和频域

1‑2参数识别。在时域内识别结构的物理参数,主要方法有扩展卡尔曼滤波(Extended Kalman filter, EKF)、H3、蒙特卡洛滤4和最小二乘估计(Least square estimate, LSE5‑6。EKF是在经典卡尔曼滤波(Kalman filter, KF7‑8的基础上改进得到的,文献[9‑11]使用扩展卡尔曼滤波方法对多自由度线性结构的参数进行识别。文献[12‑13]为了实现对结构时变参数的跟踪,采用自适应技术对EKF方法进行了改进。文献[14‑15]在载荷未知的情况下,将最小二乘法与EKF方法结合,识别时变的结构参数,同时求出未知载荷。

这些文献都对扩展卡尔曼滤波方法进行了改进,但并没有提及滤波参数的设置对识别效果的影响,尤其是结构中参数会随时间变化的情况。事实上,随着算法的迭代计算,EKF算法对测点响应的利用程度会越来越小,此时识别结果无法追踪到时变参数的变化情况。因此,必须选择恰当的滤波参数(模型噪声的协方差矩阵Q,测量噪声的协方差矩阵R)。本文首先介绍了一种结构参数识别方法,适用于外界载荷已知的情况,并试图通过仿真算例说明如何选择最优的滤波参数,来使得EKF算法对时变结构参数的跟踪性能达到最佳,即算法能快速、准确地识别出结构参数的变化情况。

1 结构参数识别方法

本文使用扩展卡尔曼滤波方法作为结构参数的识别方法。对于1个n自由度的动力学系统,当结构中的未知参数,例如质量矩阵M,刚度矩阵K或阻尼矩阵C中的某些参数待确定时,参考卡尔曼滤波理论,将这些未知参数与状态向量x一起组成扩展状态向量z

z(t)=x(t)θ=p(t)p˙(t)θ (1)

式中: θ为待确定的未知参数向量,长度设为mx(t)表示状态向量,其长度为2n。令θ˙=0,此时,卡尔曼滤波理论中的状态更新方程可以重新写为

z˙(t)=x˙(t)θ˙=p˙(t)p¨(t)[0]=p˙(t)-M-1Kp(t)-M-1Cp(t)+M-1Buf(t)[0]=              fc(z(t),f(t)) (2)

而测量更新方程的形式与观测量的类别有关,本文选择的测量类别是加速度信号。此时,测量更新方程可以写为

y(t)=H0M-1-Cp˙(t)-Kp(t)+Buf(t)=               hz(t),f(t) (3)

式中: M表示大小为n×n的质量矩阵; C表示大小为n×n的阻尼矩阵; K表示大小为的n×n的刚度矩阵; p¨(t)p˙(t)p(t)分别代表结构的加速度向量,速度向量,以及位移向量,向量的长度为nBu代表激励的影响矩阵,与激励在结构上的作用位置有关,由0和1组成,它的矩阵大小为s×nf(t)表示激励向量,长度为sH0为由0和1组成的位置矩阵。由于质量矩阵M,刚度矩阵K或阻尼矩阵C包含了未知参数θ,这也就意味着关于扩展状态向量z(t)的状态更新方程fc(z(t),f(t))是一个非线性方程,同样测量更新方程h(z(t),f(t))也是一个非线性方程。

利用泰勒展开,将上述两个非线性方程(1, 2)转化为多项式之和的形式,并略去高阶项,仅保留一阶多项式,由此可将非线性方程近似看成线性方程

fc(z(t),f(t))fc(zk-1|k-1,fk-1)+zfk-1c(z(t)-zk-1|k-1) (4)
h(z(t),f(t))h(zk|k-1,fk)+zhk(z(t)-zk|k-1) (5)

式中: zk-1|k-1为在tk-1时刻的扩展状态向量;zk|k-1为在tk时刻的扩展状态向量; zfk-1c为状态更新方程对扩展状态向量z(t)的一阶偏导; zhk为测量更新方程对扩展状态向量z(t)的一阶偏导,又称为雅各比矩阵。由泰勒展开定理可知,只要t时刻的扩展状态向量z(t)zk-1|k-1zk|k-1之间大小差距很小,式(45)就能近似成立。而要使z(t)zk-1|k-1zk|k-1的差距很小,只要时刻t与时刻tk-1tk差距很小即可。

得到离散化的近似线性等式(45)后,与卡尔曼滤波理论的方法类似。首先是状态更新过程,其先验估计的均值和估计误差的方差可以写为

zk|k-1=zk-1|k-1+tk-1tkfc(z(t),f(t))dt=zk-1|k-1+Δtfc(zk-1|k-1,f(tk-1),tk-1) (6)
Pk|k-1=Cov[zk-zk|k-1]=Cov[zk-zk|k-1]=       Cov[zk-1+tk-1tkfc(z(t),f(t))dt+w-zk-1|k-1-Δtfc(zk-1|k-1,fk-1)]=       (I+Δtzfk-1c)Pk-1|k-1(I+Δtzfk-1c)T+Q (7)

式(7)可知,扩展状态向量中的参数θ在状态更新过程中并不会变化,仅凭该过程并不能实现时变参数的追踪然后对于测量更新方程,后验估计的均值和估计误差的方差可以写为

zk|k=zk|k-1+Kk(yk-h(zk|k-1,fk)) (8)
Pk|k=Cov[zk-zk|k]=Cov[zk-zk|k]=     Cov[zk-zk|k-1-Kk(yk-h(zk|k-1,fk))]=      (I+Kkzhk)Pk|k-1(I+Kkzhk)T+KkRKkT (9)

式中:Kk为卡尔曼滤波增益,它表示对测量响应的利用程度。通过方程(8,9)可以使得扩展状态向量zk|k中的参数θ发生变化。因此,可以认为扩展卡尔曼滤波算法是通过测量更新过程实现了对参数的追踪,只要加大对测量响应的利用程度,就能使zk|k较快地收敛到真实参数值。Kk其求法与卡尔曼理论中的相同,可以写为

Kk=Pk|k-1(zhk)TzhkPk|k-1(zhk)T+Rk-1 (10)

zfk-1czhk分别是非线性表达式fc(z(t),f(t))hz(t),f(t)对扩展状态向量z(t)求偏导而得到的矩阵,即雅可比矩阵。若扩展状态向量z(t)所包含的未知参数θ仅为刚度或阻尼,则zfk-1czhk可以分别写为

                 zfk-1c=fcpfcp˙fcαz=zk-1k-1=0I00-M-1K-M-1C-M-1Kα1pk-1k-1-M-1Cα1p˙k-1k-1-M-1Kαmpk-1k-1-M-1Cαmp˙k-1k-10000 (11)
          zhk=hphp˙hαz=zkk-1=-H0M-1K-H0M-1C          -H0M-1Kα1pkk-1-H0M-1Cα1p˙kk-1-H0M-1Kαspkk-1     -H0M-1Cαsp˙kk-1 (12)

由此,可以得到扩展卡尔曼滤波算法的流程如表1所示。

表1  扩展卡尔曼滤波法算法流程
Table 1  Algorithm flow of extended Kalman filter

1.给定初始条件

x0|0P0|0

2.时间更新

zk|k-1=zk-1|k-1+tk-1tkfc(z(t),f(t))dt

Pk|k-1=(I+Δtzfk-1c)Pk-1|k-1(I+Δtzfk-1c)T+Q

3.测量更新

Kk=Pk|k-1(zhk)TzhkPk|k-1(zhk)T+Rk-1

zk|k=zk|k-1+Kk(yk-h(zk|k-1,fk))

Pk|k=(I+Kkzhk)Pk|k-1(I+Kkzhk)T+KkRKkT

2 EKF算法追踪结构时变参数的影响因素

在前文状态更新方程的推导中,假定未知参数θ不会随时间变化,即θ˙=0。这也就意味着仅凭状态更新方程,只能识别出时不变的参数,对时变参数或缓变参数是无法识别出来的。对于变化参数的追踪,是通过测量更新方程来实现的。卡尔曼滤波增益Kk表示了对测量更新方程的利用程度。若希望能实现对时变参数的追踪,则对测量更新方程的利用程度必须加大,总的宗旨是能尽量利用测量更新方程,即令Kk在迭代过程中仍保持较大值。

由卡尔曼滤波增益的计算公式可知,卡尔曼增益Kk由测量噪声方差RPk|k-1共同决定,初始时刻Kk偏大,对状态的估计主要依赖于测量更新方程。随着迭代过程的进行,测量对状态的修正作用不断减少,系统反而依赖于状态更新方程,此时Kk会逐渐变小。

为了使得滤波过程中,一直保持对测量更新方程的利用,需要选择合适的滤波参数。比如:模型噪声的协方差矩阵Q,对应状态向量部分的取值较小。对应参数部分的取值较大,这相当于告知滤波系统对结构未知参数的估计存在较大偏差,促使算法识别的结构参数会向真实值靠拢;测量噪声的协方差矩阵R的取值偏小,这相当于承认滤波系统所测量的响应具有很高的置信度,加大系统对测量更新方程的利用程度,从而加强算法对结构时变参数的追踪能力;然后是外界的测量噪声对参数识别效果的影响,算法能在响应存在噪声干扰的情况下识别出结构未知参数,这是通过调整测量噪声协方差矩阵R的方式来实现的;最后是观测点数的影响,选择悬臂梁作为仿真算例,由于悬臂梁既有转角自由度,又有竖向位移自由度,选择仅仅观测竖向位移自由度可以看出算法仍能较好地识别出结构参数变化的情况。

3 仿真算例

本文选择一个三自由度系统,如图1所示,讨论模型噪声的协方差矩阵Q,测量噪声的协方差矩阵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。并假设这个三自由度系统的阻尼是比例阻尼,且有C=αM+βKα=0.05, β=0.02。在1.5 s时该结构中的刚度k1从200 N / m突然下降到150 N / m。假设激励作用的位置和形式已知,质量块m1受到激励f=sin(5π·t)+2·sin(2π·t),由此可计算得到结构各自由度的响应。为了便于讨论,选取加速度响应作为测量向量,且不施加噪声到测量向量中去。使用本文提出的扩展卡尔曼滤波算法,此时激励已知(作用位置,幅值大小),而未知参数仅选择k1,则扩展状态向量z=[p  p˙  k1]T是一个7×1列向量。

由于参数是未知的,所以初始参数的选择必然会有较大误差,而结构的初始位移和速度(即状态向量)却可以事先测得,所以可以设定成真实值。相应的初始估计误差协方差P0|0关于状态向量的部分可以设定为一个较小的值,初始估计误差协方差P0|0关于参数的部分必须设定为一个较大的值。令初始扩展状态向量z0,初始误差方差P0|0分别写为

z0=000000150 P0|0=10-610-610-610-610-610-6103

3.1 模型噪声协方差矩阵Q的影响

讨论模型噪声协方差矩阵Q的取值对时变参数识别的影响。Q为一个7阶对角矩阵,前6阶对应的是状态更新方程对结构的位移和速度所产生的误差,第7阶代表的是状态更新方程对参数k1所产生的误差。两者是不一样的,因此必须分开讨论。令测量噪声的协方差矩阵R为对角元素为10-4的三阶对角矩阵。首先考虑模型噪声的协方差矩阵Q参数部分取值的影响,假定Q关于状态向量部分的取值为10-10,关于参数部分的取值为x,即

Q=10-1010-1010-1010-1010-1010-10x

x分别取值为1,10-110-210-3时,讨论其参数部分的取值对参数识别效果的影响。同样,为了讨论状态向量部分取值的影响,令Q关于参数部分的取值为10-1,关于状态向量部分的取值为x。EKF算法识别参数的效果如图2所示。

  

  

图2 模型噪声的协方差矩阵Q对参数识别效果的影响

Fig.2 Influence of covariance matrix Q of model noise on parameter identification effect

图2可以看出,模型误差Q关于参数部分倾向于取一个偏大的值,模型误差Q关于状态向量的部分倾向于取一个偏小的值。它本质上是使Kk保持一个较大的值,算法对参数的追踪性能也就会比较好;但Kk过大,会使得识别参数出现较大的波动。在本算例中,状态向量部分取值10-10是合适的;参数部分取值10-1是合适的。这样的参数设置能很好地追踪到参数的变化情况,且参数识别图也不会有太大的波动。

3.2 测量噪声协方差矩阵R的影响

对于测量噪声的协方差矩阵R,其最优取值与测量向量中的噪声有关。现讨论对参数识别效果的影响。令模型噪声的协方差矩阵Q状态向量部分取值10-10,参数部分取值10-1。令测量噪声的协方差矩阵R为对角元素为x的三阶对角矩阵,当x取不同值时,参数的识别结果如图3所示。

图3  测量噪声的协方差矩阵R对参数识别效果的影响

Fig.3  Influence of covariance matrix R of measurement noise on parameter identification effect

测量噪声协方差矩阵R取一个较小值时,此时卡尔曼滤波增益保持一个较大的值,算法的追踪参数性能较好。但值得注意的是,R的取值不能太小,它需要与观测向量中的噪声相适应,如果选择过小的值,导致这些估计的误差超过了实际误差,识别的参数必然不会与实际情况相同,严重时会发生发散。因此,需要为测量噪声协方差矩阵R取一个合适的值,确保算法能较快地识别出结构参数值,识别结果不会发散。

3.3 噪声影响

当响应中存在噪声干扰时,本文提出的算法仍能正确识别出参数变化情况。在此分别讨论测量响应中不含噪声与含有5%噪声的情况。令模型噪声的协方差矩阵Q状态向量部分取值10-10,参数部分取值10-1。由前文分析可知,测量噪声的协方差矩阵R的最优取值与噪声水平的大小有关。对于响应中不含噪声时,设定测量噪声的协方差矩阵R是对角元素为10-4的三阶对角矩阵,可以看出参数识别效果很好;而对于响应含有5%噪声时,此时R的取值必须大于10-4,与观测噪声水平相适应。若选取R是对角元素为10-210-3的三阶对角矩阵,可以看出算法能较好地识别出参数变化情况。

图4可以看出,无论测量向量中是否包含噪声,算法实践对时变参数的追踪是通过调整测量噪声协方差矩阵R的方式来实现的。

  

  

图4 测量噪声对参数识别效果的影响

Fig.4 Effect of measurement noise on parameter identification

3.4 观测点数影响

为说明观测点数的影响,取一个悬臂梁作为仿真对象,梁长度 L=1 m,矩形横截面尺寸为宽 b=0.1 m,高 h=0.01 m,材料密度ρ = 2 770 kg/m3,弹性模量 E=71e+9 Pa,设悬臂梁的阻尼特性为瑞利阻尼,且有α=1,β=8e-7。再由公式I=bh3/12A=bh可以推得该悬臂梁的截面惯性矩I与面积A。该悬臂梁的示意如图5所示。

图5  悬臂梁结构

Fig.5  Cantilever beam structure

图5可知,该悬臂梁被划分为4个单元,它的总体刚度矩阵和总体质量矩阵是由各个单元的刚度矩阵Ki与单元质量矩阵Mi组合而成的。因此该悬臂梁既包括转角自由度,又包括竖向位移自由度,即4个转角自由度和4个竖向位移自由度。其中KiMi可写为

Ki=EiIil3126l-126l6l4l2-6l2l2-12-6l12-6l6l2l2-6l4l2 
Mi=ρiAil420210000017l20000210000017l2

为了便于讨论,把梁单元的刚度特性EiIi作为单元发生刚度变化的指标。假定梁单元①的刚度特性E1I1在第1.5~3.5 s发生线性下降,下降至原来值的70%。其他单元的参数并不发生变化。在响应中添加1%的白噪声,利用本文所提出的算法,把该悬臂梁前3个单元的刚度特性作为未知参数向量,即扩展状态向量为z=[p,p˙,E1I1,E2I2,E3I3]T,依次识别出各个刚度系数。

选择只观测第2,3,4,5节点上的竖直加速度,不测量转角自由度,因此这也算是部分观测。设模型噪声的协方差矩阵G,以及测量噪声的协方差矩阵R分别为

G=10-510-510-5102102102
R=10-310-310-310-3

观测全部测点与部分测点的区别仅需改变测量方程中的位置矩阵H0。若选择只记录结构的竖向加速度响应,因此可写为

H0=1000001000000000   0000000010000010

同理可求出既观测转角自由度,又观测竖向自由度情况下的参数识别效果,由此可得参数识别结果如图6~8所示。

图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

由图6~8可以看出,在含有噪声的情况下,仅观测部分测点也能追踪到结构参数的变化情况,且效果比较好。这说明了扩展卡尔曼滤波算法能在测点部分观测的情况下识别参数变化。

4 结 论

在载荷已知的情况下,扩展卡尔曼滤波算法对结构参数的识别主要是通过测量更新方程实现的,而卡尔曼滤波增益Kk表示了算法对测量更新方程的利用程度。为了实现算法对时变参数的追踪,本文提出了通过设置合适的滤波参数来保持较大的Kk。这其中包括:模型误差Q关于参数部分取较大值,模型误差Q关于状态向量的部分取较小的值。还讨论了测量噪声对参数识别效果的影响,测量噪声的协方差矩阵R必须与测量噪声相适应,并尽量选取一个较小值。仿真结果表明,这样的滤波参数设置使得EKF算法在运行过程中能保持对测量更新方程的利用,从而能追踪到结构参数的变化情况,识别结构的时变参数。最后,通过悬臂梁仿真算例,证明了在部分测点观测的情况下,算法也能较好地识别出结构的参数变化情况,体现了EKF算法的优越性。

参考文献

1

BERNAL DGUNES B. Flexibility based approach for damage characterization: Benchmark application[J]. Journal of Engineering Mechanics20041301):61-70. [百度学术] 

2

FRISWELL M IPENNY J E T. The choice of orthogonal polynomials in the rational fraction polynomial method[J]. Modal Analysis199383):257-262. [百度学术] 

3

SATO TQI K.Adaptive H1 filter:Its application to structural identification[J]. Journal of Engineering Mechanics (ASCE)199812411):1233-1240. [百度学术] 

4

YOSHIDA I. Damage detection using Monte Carlo filter based on non-Gaussian noise[C]∥ Proceedings of Structural Safety and ReliabilityICOSSA 2001. LisseSwet & Zeitinger,2001: 8. [百度学术] 

5

YANG J NHUANG H. Sequential non-linear least-square estimation for damage identification of structures[J]. International Journal of Non-Linear Mechanics2006411): 124-140. [百度学术] 

6

LOH C HLIN C YHUANG C C. Time domain identification of frames under earthquake loadings[J].Journal of Engineering Mechanics20001264):693-703. [百度学术] 

7

黄小平王岩.卡尔曼滤波原理及应用MATLAB仿真[M].北京电子工业出版社2015 [百度学术] 

Huang XiaopingWang Yan. Principle of Kalman filter and its application in MATLAB simulation[M]. BeijingElectronic Industry Press2015. [百度学术] 

8

KALMAN R E. A new approach to linear filtering and prediction problems[J]. Journal of Fluids Engineering1960821): 35-45. [百度学术] 

9

HOSHIYA MSAITO E. Structural identification by extended Kalman filter[J]. Journal of Engineering Mechanics198411012): 1757-1770. [百度学术] 

10

尚久锉.卡尔曼滤波法在结构动态参数估计中的应用[J].地震工程与工程震动1991112): 62-72. [百度学术] 

Shang Jiucuo. Application of Kalman filter in structural dynamic parameter estimation[J]. Earthquake engineering and Engineering Vibration1991112): 62-72. [百度学术] 

11

李杰.基于微分算子变换的广义卡尔曼估计方法[J].计算结构力学及其应用1995114): 394-400. [百度学术] 

Li Jie. Generalized Kalman estimation method based on differential operator transformation[J]. Computational Structural Mechanics and Its Application1995114): 394-400. [百度学术] 

12

张浩.基于强跟踪无迹卡尔曼滤波的结构时变参数识别[D].兰州兰州理工大学2016. [百度学术] 

Zhang Hao. Time varying structural parameter identification based on strong tracking unscented Kalman filter[D]. LanzhouLanzhou University of Technology2016. [百度学术] 

13

穆腾飞.基于自适应追踪技术的迟滞非线性系统参数识别与损伤检测研究[D].南京南京航空航天大学2016. [百度学术] 

Mu Tengfei. Parameter identification and damage detection of hysteretic nonlinear systems based on adaptive tracking technique[D]. NanjingNanjing University of Aeronautics and Astronautics2016. [百度学术] 

14

YANG J NPAN SHUANG H. An adaptive extended Kalman filter for structural damage identifications II:Unknown inputs[J]. Struct Control Health Monit2007143):497-521. [百度学术] 

15

HUANG QXU Y LLIU H J. An efficient algorithm for simultaneous identification of time-varying structural parameters and unknown excitations of a building structure[J].Engineering Structures20159829-37. [百度学术] 

您是第位访问者
网站版权 © 南京航空航天大学学报
技术支持:北京勤云科技发展有限公司
请使用 Firefox、Chrome、IE10、IE11、360极速模式、搜狗极速模式、QQ极速模式等浏览器,其他浏览器不建议使用!