擴展卡爾曼濾波對時變參數追蹤性能的影響研究

2022-04-27 01:45:38唐宏志姜金輝

南京航空航天大學學報 2022年2期

唐宏志，姜金輝

（南京航空航天大學機械結構力學及控制國家重點實驗室，南京 210016）

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

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

1 結構參數識別方法

本文使用擴展卡爾曼濾波方法作為結構參數的識別方法。對于1 個n自由度的動力學系統，當結構中的未知參數，例如質量矩陣M，剛度矩陣K或阻尼矩陣C中的某些參數待確定時，參考卡爾曼濾波理論，將這些未知參數與狀態向量x一起組成擴展狀態向量z

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

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

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

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

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

由式（7）可知，擴展狀態向量中的參數θ在狀態更新過程中并不會變化，僅憑該過程并不能實現時變參數的追蹤。然后對于測量更新方程，后驗估計的均值和估計誤差的方差可以寫為

式中：Kk為卡爾曼濾波增益，它表示對測量響應的利用程度。通過方程（8，9）可以使得擴展狀態向量zk|k中的參數θ發生變化。因此，可以認為擴展卡爾曼濾波算法是通過測量更新過程實現了對參數的追蹤，只要加大對測量響應的利用程度，就能使zk|k較快地收斂到真實參數值。Kk其求法與卡爾曼理論中的相同，可以寫為

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

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

2 EKF 算法追蹤結構時變參數的影響因素

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

由卡爾曼濾波增益的計算公式可知，卡爾曼增益Kk由測量噪聲方差R和Pk|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分別寫為

3.1 模型噪聲協方差矩陣Q 的影響

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

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

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

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

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-2或10-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=b·h3/12，A=b·h可以推得該懸臂梁的截面慣性矩I與面積A。該懸臂梁的示意如圖5 所示。

圖5 懸臂梁結構Fig.5 Cantilever beam structure

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

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

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

觀測全部測點與部分測點的區別僅需改變測量方程中的位置矩陣H0。若選擇只記錄結構的豎向加速度響應，因此可寫為

同理可求出既觀測轉角自由度，又觀測豎向自由度情況下的參數識別效果，由此可得參數識別結果如圖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 算法的優越性。