擴(kuò)展卡爾曼濾波對時變參數(shù)追蹤性能的影響研究

唐宏志，姜金輝

（南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點實驗室，南京 210016）

參數(shù)識別是動力學(xué)的第一類逆問題。參數(shù)識別問題按識別參數(shù)的種類可以分為模態(tài)參數(shù)識別與物理參數(shù)識別；而如果按識別參數(shù)的方法，可以分為時域法參數(shù)識別和頻域法［1?2］參數(shù)識別。在時域內(nèi)識別結(jié)構(gòu)的物理參數(shù)，主要方法有擴(kuò)展卡爾曼濾波（Extended Kalman filter，EKF）、H∞濾波［3］、蒙特卡洛濾波［4］和最小二乘估計（Least square esti?mate，LSE）［5?6］。EKF 是在經(jīng)典卡爾曼濾波（Kal?man filter，KF）［7?8］的基礎(chǔ)上改進(jìn)得到的，文獻(xiàn)［9?11］使用擴(kuò)展卡爾曼濾波方法對多自由度線性結(jié)構(gòu)的參數(shù)進(jìn)行識別。文獻(xiàn)［12?13］為了實現(xiàn)對結(jié)構(gòu)時變參數(shù)的跟蹤，采用自適應(yīng)技術(shù)對EKF 方法進(jìn)行了改進(jìn)。文獻(xiàn)［14?15］在載荷未知的情況下，將最小二乘法與EKF 方法結(jié)合，識別時變的結(jié)構(gòu)參數(shù)，同時求出未知載荷。

這些文獻(xiàn)都對擴(kuò)展卡爾曼濾波方法進(jìn)行了改進(jìn)，但并沒有提及濾波參數(shù)的設(shè)置對識別效果的影響，尤其是結(jié)構(gòu)中參數(shù)會隨時間變化的情況。事實上，隨著算法的迭代計算，EKF 算法對測點響應(yīng)的利用程度會越來越小，此時識別結(jié)果無法追蹤到時變參數(shù)的變化情況。因此，必須選擇恰當(dāng)?shù)臑V波參數(shù)（模型噪聲的協(xié)方差矩陣Q，測量噪聲的協(xié)方差矩陣R）。本文首先介紹了一種結(jié)構(gòu)參數(shù)識別方法，適用于外界載荷已知的情況，并試圖通過仿真算例說明如何選擇最優(yōu)的濾波參數(shù)，來使得EKF算法對時變結(jié)構(gòu)參數(shù)的跟蹤性能達(dá)到最佳，即算法能快速、準(zhǔn)確地識別出結(jié)構(gòu)參數(shù)的變化情況。

1 結(jié)構(gòu)參數(shù)識別方法

本文使用擴(kuò)展卡爾曼濾波方法作為結(jié)構(gòu)參數(shù)的識別方法。對于1 個n自由度的動力學(xué)系統(tǒng)，當(dāng)結(jié)構(gòu)中的未知參數(shù)，例如質(zhì)量矩陣M，剛度矩陣K或阻尼矩陣C中的某些參數(shù)待確定時，參考卡爾曼濾波理論，將這些未知參數(shù)與狀態(tài)向量x一起組成擴(kuò)展?fàn)顟B(tài)向量z

式中：θ為待確定的未知參數(shù)向量，長度設(shè)為m；x(t)表示狀態(tài)向量，其長度為2n。令θ?=0，此時，卡爾曼濾波理論中的狀態(tài)更新方程可以重新寫為

而測量更新方程的形式與觀測量的類別有關(guān)，本文選擇的測量類別是加速度信號。此時，測量更新方程可以寫為

式中：M表示大小為n×n的質(zhì)量矩陣；C表示大小為n×n的阻尼矩陣；K表示大小為的n×n的剛度矩陣；p?(t)，p?(t)，p(t)分別代表結(jié)構(gòu)的加速度向量，速度向量，以及位移向量，向量的長度為n；Bu代表激勵的影響矩陣，與激勵在結(jié)構(gòu)上的作用位置有關(guān)，由0 和1 組成，它的矩陣大小為s×n；f(t)表示激勵向量，長度為s；H0為由0 和1 組成的位置矩陣。由于質(zhì)量矩陣M，剛度矩陣K或阻尼矩陣C包含了未知參數(shù)θ，這也就意味著關(guān)于擴(kuò)展?fàn)顟B(tài)向量z(t)的狀態(tài)更新方程fc(z(t)，f(t))是一個非線性方程，同樣測量更新方程h(z(t)，f(t))也是一個非線性方程。

利用泰勒展開，將上述兩個非線性方程（1，2）轉(zhuǎn)化為多項式之和的形式，并略去高階項，僅保留一階多項式，由此可將非線性方程近似看成線性方程

式中：zk-1|k-1為在tk-1時刻的擴(kuò)展?fàn)顟B(tài)向量；zk|k-1為在tk時刻的擴(kuò)展?fàn)顟B(tài)向量；?z fc k-1為狀態(tài)更新方程對擴(kuò)展?fàn)顟B(tài)向量z(t)的一階偏導(dǎo)；?zhk為測量更新方程對擴(kuò)展?fàn)顟B(tài)向量z(t)的一階偏導(dǎo)，又稱為雅各比矩陣。由泰勒展開定理可知，只要t時刻的擴(kuò)展?fàn)顟B(tài)向量z(t)與zk-1|k-1，zk|k-1之間大小差距很小，式（4，5）就能近似成立。而要使z(t)與zk-1|k-1，zk|k-1的差距很小，只要時刻t與時刻tk-1，tk差距很小即可。

得到離散化的近似線性等式（4，5）后，與卡爾曼濾波理論的方法類似。首先是狀態(tài)更新過程，其先驗估計的均值和估計誤差的方差可以寫為

由式（7）可知，擴(kuò)展?fàn)顟B(tài)向量中的參數(shù)θ在狀態(tài)更新過程中并不會變化，僅憑該過程并不能實現(xiàn)時變參數(shù)的追蹤。然后對于測量更新方程，后驗估計的均值和估計誤差的方差可以寫為

式中：Kk為卡爾曼濾波增益，它表示對測量響應(yīng)的利用程度。通過方程（8，9）可以使得擴(kuò)展?fàn)顟B(tài)向量zk|k中的參數(shù)θ發(fā)生變化。因此，可以認(rèn)為擴(kuò)展卡爾曼濾波算法是通過測量更新過程實現(xiàn)了對參數(shù)的追蹤，只要加大對測量響應(yīng)的利用程度，就能使zk|k較快地收斂到真實參數(shù)值。Kk其求法與卡爾曼理論中的相同，可以寫為

由此，可以得到擴(kuò)展卡爾曼濾波算法的流程如表1 所示。

表1 擴(kuò)展卡爾曼濾波法算法流程Table 1 Algorithm flow of extended Kalman filter

2 EKF 算法追蹤結(jié)構(gòu)時變參數(shù)的影響因素

在前文狀態(tài)更新方程的推導(dǎo)中，假定未知參數(shù)θ不會隨時間變化，即θ?=0。這也就意味著僅憑狀態(tài)更新方程，只能識別出時不變的參數(shù)，對時變參數(shù)或緩變參數(shù)是無法識別出來的。對于變化參數(shù)的追蹤，是通過測量更新方程來實現(xiàn)的。卡爾曼濾波增益Kk表示了對測量更新方程的利用程度。若希望能實現(xiàn)對時變參數(shù)的追蹤，則對測量更新方程的利用程度必須加大，總的宗旨是能盡量利用測量更新方程，即令Kk在迭代過程中仍保持較大值。

由卡爾曼濾波增益的計算公式可知，卡爾曼增益Kk由測量噪聲方差R和Pk|k-1共同決定，初始時刻Kk偏大，對狀態(tài)的估計主要依賴于測量更新方程。隨著迭代過程的進(jìn)行，測量對狀態(tài)的修正作用不斷減少，系統(tǒng)反而依賴于狀態(tài)更新方程，此時Kk會逐漸變小。

為了使得濾波過程中，一直保持對測量更新方程的利用，需要選擇合適的濾波參數(shù)。比如：模型噪聲的協(xié)方差矩陣Q，對應(yīng)狀態(tài)向量部分的取值較小。對應(yīng)參數(shù)部分的取值較大，這相當(dāng)于告知濾波系統(tǒng)對結(jié)構(gòu)未知參數(shù)的估計存在較大偏差，促使算法識別的結(jié)構(gòu)參數(shù)會向真實值靠攏；測量噪聲的協(xié)方差矩陣R的取值偏小，這相當(dāng)于承認(rèn)濾波系統(tǒng)所測量的響應(yīng)具有很高的置信度，加大系統(tǒng)對測量更新方程的利用程度，從而加強(qiáng)算法對結(jié)構(gòu)時變參數(shù)的追蹤能力；然后是外界的測量噪聲對參數(shù)識別效果的影響，算法能在響應(yīng)存在噪聲干擾的情況下識別出結(jié)構(gòu)未知參數(shù)，這是通過調(diào)整測量噪聲協(xié)方差矩陣R的方式來實現(xiàn)的；最后是觀測點數(shù)的影響，選擇懸臂梁作為仿真算例，由于懸臂梁既有轉(zhuǎn)角自由度，又有豎向位移自由度，選擇僅僅觀測豎向位移自由度可以看出算法仍能較好地識別出結(jié)構(gòu)參數(shù)變化的情況。

3 仿真算例

本文選擇一個三自由度系統(tǒng)，如圖1 所示，討論模型噪聲的協(xié)方差矩陣Q，測量噪聲的協(xié)方差矩陣R的影響，與觀測點數(shù)對算法跟蹤參數(shù)性能的影響。

圖1 三自由度系統(tǒng)Fig.1 A three degree of freedom system

在初始狀態(tài)下，k1=k4= 200 N / m，k2=k3= 100 N / m；m1=m2=m3= 1 kg。并假設(shè)這個三自由度系統(tǒng)的阻尼是比例阻尼，且有C=αM+βK，α=0.05，β=0.02。在1.5 s 時該結(jié)構(gòu)中的剛度k1從200 N / m 突然下降到150 N / m。假設(shè)激勵作用的位置和形式已知，質(zhì)量塊m1受到激勵f=sin(5π·t)+2·sin(2π·t)，由此可計算得到結(jié)構(gòu)各自由度的響應(yīng)。為了便于討論，選取加速度響應(yīng)作為測量向量，且不施加噪聲到測量向量中去。使用本文提出的擴(kuò)展卡爾曼濾波算法，此時激勵已知（作用位置，幅值大?。?，而未知參數(shù)僅選擇k1，則擴(kuò)展?fàn)顟B(tài)向量z=[p p?k1]T是一個7×1列向量。

由于參數(shù)是未知的，所以初始參數(shù)的選擇必然會有較大誤差，而結(jié)構(gòu)的初始位移和速度（即狀態(tài)向量）卻可以事先測得，所以可以設(shè)定成真實值。相應(yīng)的初始估計誤差協(xié)方差P0|0關(guān)于狀態(tài)向量的部分可以設(shè)定為一個較小的值，初始估計誤差協(xié)方差P0|0關(guān)于參數(shù)的部分必須設(shè)定為一個較大的值。令初始擴(kuò)展?fàn)顟B(tài)向量z0，初始誤差方差P0|0分別寫為

3.1 模型噪聲協(xié)方差矩陣Q 的影響

討論模型噪聲協(xié)方差矩陣Q的取值對時變參數(shù)識別的影響。Q為一個7 階對角矩陣，前6 階對應(yīng)的是狀態(tài)更新方程對結(jié)構(gòu)的位移和速度所產(chǎn)生的誤差，第7 階代表的是狀態(tài)更新方程對參數(shù)k1所產(chǎn)生的誤差。兩者是不一樣的，因此必須分開討論。令測量噪聲的協(xié)方差矩陣R為對角元素為10-4的三階對角矩陣。首先考慮模型噪聲的協(xié)方差矩陣Q參數(shù)部分取值的影響，假定Q關(guān)于狀態(tài)向量部分的取值為10-10，關(guān)于參數(shù)部分的取值為x，即

x分別取值為1，10-1，10-2，10-3時，討論其參數(shù)部分的取值對參數(shù)識別效果的影響。同樣，為了討論狀態(tài)向量部分取值的影響，令Q關(guān)于參數(shù)部分的取值為10-1，關(guān)于狀態(tài)向量部分的取值為x。EKF 算法識別參數(shù)的效果如圖2 所示。

從圖2 可以看出，模型誤差Q關(guān)于參數(shù)部分傾向于取一個偏大的值，模型誤差Q關(guān)于狀態(tài)向量的部分傾向于取一個偏小的值。它本質(zhì)上是使Kk保持一個較大的值，算法對參數(shù)的追蹤性能也就會比較好；但Kk過大，會使得識別參數(shù)出現(xiàn)較大的波動。在本算例中，狀態(tài)向量部分取值10-10是合適的；參數(shù)部分取值10-1是合適的。這樣的參數(shù)設(shè)置能很好地追蹤到參數(shù)的變化情況，且參數(shù)識別圖也不會有太大的波動。

圖2 模型噪聲的協(xié)方差矩陣Q 對參數(shù)識別效果的影響Fig.2 Influence of covariance matrix Q of model noise on parameter identification effect

3.2 測量噪聲協(xié)方差矩陣R 的影響

對于測量噪聲的協(xié)方差矩陣R，其最優(yōu)取值與測量向量中的噪聲有關(guān)。現(xiàn)討論對參數(shù)識別效果的影響。令模型噪聲的協(xié)方差矩陣Q狀態(tài)向量部分取值10-10，參數(shù)部分取值10-1。令測量噪聲的協(xié)方差矩陣R為對角元素為x的三階對角矩陣，當(dāng)x取不同值時，參數(shù)的識別結(jié)果如圖3 所示。

圖3 測量噪聲的協(xié)方差矩陣R 對參數(shù)識別效果的影響Fig.3 Influence of covariance matrix R of measurement noise on parameter identification effect

測量噪聲協(xié)方差矩陣R取一個較小值時，此時卡爾曼濾波增益保持一個較大的值，算法的追蹤參數(shù)性能較好。但值得注意的是，R的取值不能太小，它需要與觀測向量中的噪聲相適應(yīng)，如果選擇過小的值，導(dǎo)致這些估計的誤差超過了實際誤差，識別的參數(shù)必然不會與實際情況相同，嚴(yán)重時會發(fā)生發(fā)散。因此，需要為測量噪聲協(xié)方差矩陣R取一個合適的值，確保算法能較快地識別出結(jié)構(gòu)參數(shù)值，識別結(jié)果不會發(fā)散。

3.3 噪聲影響

當(dāng)響應(yīng)中存在噪聲干擾時，本文提出的算法仍能正確識別出參數(shù)變化情況。在此分別討論測量響應(yīng)中不含噪聲與含有5%噪聲的情況。令模型噪聲的協(xié)方差矩陣Q狀態(tài)向量部分取值10-10，參數(shù)部分取值10-1。由前文分析可知，測量噪聲的協(xié)方差矩陣R的最優(yōu)取值與噪聲水平的大小有關(guān)。對于響應(yīng)中不含噪聲時，設(shè)定測量噪聲的協(xié)方差矩陣R是對角元素為10-4的三階對角矩陣，可以看出參數(shù)識別效果很好；而對于響應(yīng)含有5%噪聲時，此時R的取值必須大于10-4，與觀測噪聲水平相適應(yīng)。若選取R是對角元素為10-2或10-3的三階對角矩陣，可以看出算法能較好地識別出參數(shù)變化情況。

由圖4 可以看出，無論測量向量中是否包含噪聲，算法實踐對時變參數(shù)的追蹤是通過調(diào)整測量噪聲協(xié)方差矩陣R的方式來實現(xiàn)的。

圖4 測量噪聲對參數(shù)識別效果的影響Fig.4 Effect of measurement noise on parameter identification

3.4 觀測點數(shù)影響

為說明觀測點數(shù)的影響，取一個懸臂梁作為仿真對象，梁長度L=1 m，矩形橫截面尺寸為寬b=0.1 m，高h(yuǎn)=0.01 m，材料密度ρ= 2 770 kg/m3，彈性模量E=71e+9 Pa，設(shè)懸臂梁的阻尼特性為瑞利阻尼，且有α=1，β=8e-7。再由公式I=b·h3/12，A=b·h可以推得該懸臂梁的截面慣性矩I與面積A。該懸臂梁的示意如圖5 所示。

圖5 懸臂梁結(jié)構(gòu)Fig.5 Cantilever beam structure

由圖5 可知，該懸臂梁被劃分為4 個單元，它的總體剛度矩陣和總體質(zhì)量矩陣是由各個單元的剛度矩陣Ki與單元質(zhì)量矩陣Mi組合而成的。因此該懸臂梁既包括轉(zhuǎn)角自由度，又包括豎向位移自由度，即4 個轉(zhuǎn)角自由度和4 個豎向位移自由度。其中Ki，Mi可寫為

為了便于討論，把梁單元的剛度特性Ei Ii作為單元發(fā)生剛度變化的指標(biāo)。假定梁單元①的剛度特性E1I1在第1.5～3.5 s 發(fā)生線性下降，下降至原來值的70%。其他單元的參數(shù)并不發(fā)生變化。在響應(yīng)中添加1% 的白噪聲，利用本文所提出的算法，把該懸臂梁前3 個單元的剛度特性作為未知參數(shù)向量，即擴(kuò)展?fàn)顟B(tài)向量為z=[p，p?，E1I1，E2I2，E3I3]T，依次識別出各個剛度系數(shù)。

選擇只觀測第2，3，4，5 節(jié)點上的豎直加速度，不測量轉(zhuǎn)角自由度，因此這也算是部分觀測。設(shè)模型噪聲的協(xié)方差矩陣G，以及測量噪聲的協(xié)方差矩陣R分別為

觀測全部測點與部分測點的區(qū)別僅需改變測量方程中的位置矩陣H0。若選擇只記錄結(jié)構(gòu)的豎向加速度響應(yīng)，因此可寫為

同理可求出既觀測轉(zhuǎn)角自由度，又觀測豎向自由度情況下的參數(shù)識別效果，由此可得參數(shù)識別結(jié)果如圖6～8 所示。

圖6 剛度系數(shù)E1I1識別效果圖Fig.6 Identification effect of stiffness coefficient E1I1

圖7 剛度系數(shù)E2I2識別效果圖Fig.7 Identification effect of stiffness coefficient E2I2

圖8 剛度系數(shù)E3I3識別效果圖Fig.8 Identification effect of stiffness coefficient E3I3

由圖6～8 可以看出，在含有噪聲的情況下，僅觀測部分測點也能追蹤到結(jié)構(gòu)參數(shù)的變化情況，且效果比較好。這說明了擴(kuò)展卡爾曼濾波算法能在測點部分觀測的情況下識別參數(shù)變化。

4 結(jié)論

在載荷已知的情況下，擴(kuò)展卡爾曼濾波算法對結(jié)構(gòu)參數(shù)的識別主要是通過測量更新方程實現(xiàn)的，而卡爾曼濾波增益Kk表示了算法對測量更新方程的利用程度。為了實現(xiàn)算法對時變參數(shù)的追蹤，本文提出了通過設(shè)置合適的濾波參數(shù)來保持較大的Kk。這其中包括：模型誤差Q關(guān)于參數(shù)部分取較大值，模型誤差Q關(guān)于狀態(tài)向量的部分取較小的值。還討論了測量噪聲對參數(shù)識別效果的影響，測量噪聲的協(xié)方差矩陣R必須與測量噪聲相適應(yīng)，并盡量選取一個較小值。仿真結(jié)果表明，這樣的濾波參數(shù)設(shè)置使得EKF 算法在運行過程中能保持對測量更新方程的利用，從而能追蹤到結(jié)構(gòu)參數(shù)的變化情況，識別結(jié)構(gòu)的時變參數(shù)。最后，通過懸臂梁仿真算例，證明了在部分測點觀測的情況下，算法也能較好地識別出結(jié)構(gòu)的參數(shù)變化情況，體現(xiàn)了EKF 算法的優(yōu)越性。