動態系統加性干擾與狀態估計的實時濾波方法

陳健,文成林

(1.吉林化工學院 信息與控制工程學院,吉林 132022;2.廣東石油化工學院 自動化學院,廣東 茂名 510006;3.中國科學院深圳先進技術研究院 廣東省機器人與智能系統重點實驗室,廣東 深圳 518055)

隨著濾波器的不斷發展與改善,濾波器設計和狀態估計問題逐漸成為研究的熱點之一[1]。實際應用中,最具有代表性的是線性白噪聲系統下的卡爾曼濾波,它是基于均方誤差最小準則下的遞歸最優狀態估計方法[2,3]。當系統呈現為弱非線性時,一般將非線性系統線性化,轉化為符合線性卡爾曼濾波器的標準形式,再進行相應的濾波器設計[4-8]。模型中誤差變動的影響,大大增加了模型的不確定性,導致濾波器的性能下降,容易引起濾波的發散[9]。因此,針對模型的不確定而造成的濾波器性能下降,進而引起濾波器性能失效的問題,本文建立了相應的多不確定性下的新型濾波器。本文主要圍繞一類線性系統中因建模誤差引起模型不確定情況進行研究,首先將不確定因子視為加性干擾[10-12],通過設計出求取狀態模型中不確定因子當前值以及未來值的輔助方程,再基于未來觀測信息的順序求出狀態模型中不確定因子的當前值和測量模型中不確定因子的未來值,進一步優化求解未來狀態中基于測量模型和狀態模型在線調整的自適應濾波器[13]。

1 線性系統加性濾波設計方法

1.1 問題陳述

考慮如下隨機系統

x(k+1)=A(k)x(k)+D(k)d(k)+w(k)

(1)

y(k+1)=B(k+1)x(k+1)+C(k+1)r(k+1)+v(k+1)

(2)

式中:k為離散時間;A(k)為系統狀態轉移矩陣;B(k+1)為觀測矩陣;D(k)和C(k+1)為系數矩陣;x(k)∈Rn為狀態向量;y(k+1)∈Rm為觀測向量;d(k)和r(k+1)為系統中的不確定因子。

假設系統建模誤差w(k)和v(k+1)為不相關的白噪聲序列,滿足

E{w(k)}=0,E{v(k)}=0

(3)

E{w(k)wT(j)}=Q(k)δkj,E{v(k)vT(j)}=R(k+1)δkj

(4)

式中:Q(k)為噪聲w(k)的方差;R(k+1)為噪聲v(k+1)的方差;δkj為Kronecker積。

根據式(1)和式(2)可以看出,狀態變量x(k+1)的變化受加性干擾d(k)和r(k+1)的影響。因此,需要利用d(k)和r(k+1)的估計值對狀態x(k+1)進行估計。

1.2 加性干擾d(k)的Kalman濾波器設計

觀測模型y(k+1)沒有關于加性干擾d(k)的顯示信息,而d(k)的變化會反映在狀態x(k+1)中,因此,以隨機游走的形式建立一步滯后的動態模型,可得

d(k)=Ad(k-1)d(k-1)+wd(k-1)

(5)

為了對加性干擾d(k)精確估計,觀測模型簡化為

y(k+1)=B(k+1)(A(k)x(k)+D(k)d(k)+w(k))+C(k+1)r(k+1)+v(k+1)

(6)

假設如下統計特性

(7)

式中:Ad(k-1)為狀態轉移矩陣;wd(k-1)為零均值的高斯白噪聲;Qd(k-1)為噪聲wd(k-1)的方差。

(8)

根據標準卡爾曼濾波器設計原理,加性干擾d(k)的濾波器設計可以分為時間更新和測量更新兩個步驟,具體如下。

步驟一:時間更新

通過式(5),得到關于d(k)的預測值和預測誤差協方差矩陣

(9)

(10)

步驟二:測量更新

根據式 (6),可得加性干擾d(k)的測量預測值及測量預測誤差分別為

(11)

(12)

則d(k)的濾波器設計為

(13)

相應的估計誤差為

(14)

結合正交原理得到濾波器中Kd(k)

Kd=Pd(k|k-1)DT(k+1)BT(k+1)[B(k+1)A(k)Px(k|k)AT(k)BT(k+1)

+B(k+1)Q(k)BT(k+1)+C(k+1)Pr(k+1|k)CT(k+1)

+R(k)+B(k+1)D(k)Pd(k|k-1)DT(k)BT(k+1)]-1

(15)

最后,根據式(14)求得估計誤差協方差矩陣為

Pd(k|k)=[I-Kd(k)B(k+1)D(k)]Pd(k|k-1)

(16)

至此,狀態模型中的加性干擾d(k)的濾波器設計就完成了,改寫測量模型,使模型中出現d(k)的顯示信息,本質上是原始測量模型的等價改寫。

1.3 加性干擾r(k+1)的Kalman濾波器設計

本節的目標是設計觀測模型中加性干擾r(k+1)的Kalman濾波器,同樣以隨機游走的形式建立動態模型即

r(k+1)=Br(k)r(k)+wr(k)

(17)

假設如下統計特性

(18)

式中:Br(k)為狀態轉移矩陣;wr(k)為零均值的高斯白噪聲;Qr(k)為噪聲wr(k)的方差。

根據標準卡爾曼濾波器設計原理,加性干擾r(k+1)的濾波器設計過程分兩個步驟,具體如下。

步驟一:時間更新

基于r(k+1)的動態模型,得到r(k+1)的預測值及預測誤差協方差矩陣

(19)

步驟二:測量更新

根據式(2)求得加性干擾r(k+1)的測量預測值和測量預測誤差分別為

(20)

(21)

則r(k+1)的濾波器設計為

(22)

以及相應的估計誤差為

(23)

最后,利用正交原理再結合估計誤差以及觀測模型可得增益矩陣,并根據式(27)得到Pr(k+1|k+1)

Kr(k+1)=P(k+1|k)CT(k+1)[B(k+1)Px(k+1|k)BT(k+1)

+C(k+1)Pr(k+1|k)CT(k+1)+R(k+1)]-1

(24)

(25)

至此,我們完成了狀態模型和觀測模型中的加性干擾r(k+1)和d(k)的Kalman設計過程。

1.4 求解系統狀態變量x(k+1)估計的Kalman濾波器

本節的主要目的是在已獲得的加性干擾d(k)和r(k+1)估計值的基礎上,基于原始狀態模型和測量模型,求解系統狀態狀態變量x(k+1)的濾波器,即

(26)

根據標準卡爾曼濾波器設計原理,濾波器的設計過程如下。

步驟一:時間更新

根據式(1)可得x(k+1)的狀態預測值、狀態預測誤差及預測估計值誤差協方差矩陣,分別為

(27)

(28)

Px(k+1|k)=A(k)Px(k|k)AT(k)+D(k)Pd(k|k)DT(k)+Q(k)

(29)

步驟二:測量更新

根據狀態預測值和觀測模型可得測量預測值和測量預測誤差,分別為

(30)

(31)

則狀態變量x(k+1)的濾波器設計為

(32)

以及相應的狀態估計誤差

(33)

式中:Kx(k+1)為增益矩陣。利用正交原理求解得增益矩陣可得

Kx(k+1)=Px(k+1|k)BT(k+1)[B(k+1)Px(k+1|k)BT(k+1)

+C(k+1)Pr(k+1|k+1)CT(k+1)+R(k+1)]-1

(34)

最后,根據式(33)求期望得到Px(k+1|k+1)

(35)

2 誤差分析

從濾波更新過程中可以看出,加性干擾d(k)和r(k+1)以及狀態x(k+1)三者互為影響,對于已經估計出來的變量,在測量方程中顯示為k+1的估計值,對于未估計出的變量,顯示為預測值,而對比現有的方法統一視為預測值。因此,與STF相比,本文的方法的實時性更好。

3 仿真驗證

考慮如下隨機系統

x1(k+1)=x1(k)+a1(k)+w1(k),y1(k+1)=x1(k+1)+0.01r1(k+1)+v1(k+1)

x2(k+1)=x2(k)+a2(k)+w2(k),y2(k+1)=x2(k+1)+0.02r2(k+1)+v2(k+1)

(36)

圖1～圖6為各變量的估計曲線以及各變量的估計誤差曲線,由式(15)～式(18)、式(22)～式(29)、式(31)～式(35)得出。表1為各變量的估計誤差求和平均后的結果。

圖1 狀態x的估計值

圖3 參數a的估計值

圖5 參數r的估計值

表1 估計誤差平均值

由圖1～圖6和表1可看出,對比STF,本文所提出的實時濾波方法具有更好的濾波性能,加性干擾a1(k)、a2(k)、r1(k+1)和r2(k+1)的估計準確率分別提升24.85%、29.01%、37.52%、31.63%,狀態變量x1(k+1) 和x2(k+1)的估計準確率分別提升16.49%和18.43%。

三者之間存在著相互影響的關系,若加性干擾的濾波準確率提高,則加性干擾對模型的影響就會減弱,從而提高原始狀態變量的估計精度,證明了本文所提方法的有效性。

4 結語

本文針對具有加性干擾的不確定系統設計了一種三階段卡爾曼濾波器。本文提出建立加性干擾d(k)當前值的動態輔助方程和建立測量模型中r(k+1)的輔助方程,結合原始測量模型進行濾波;視加性干擾d(k)和r(k+1)的估計值為測量,設計原始狀態變量x(k+1)的卡爾曼濾波器。仿真結果表明,該方法可以有效應對模型中的加性干擾,減少加性干擾對原始狀態變量的干擾,能夠有效實現濾波估計。