陳健,文成林

(1.吉林化工學(xué)院 信息與控制工程學(xué)院,吉林 132022;2.廣東石油化工學(xué)院 自動(dòng)化學(xué)院,廣東 茂名 510006;3.中國科學(xué)院深圳先進(jìn)技術(shù)研究院 廣東省機(jī)器人與智能系統(tǒng)重點(diǎn)實(shí)驗(yàn)室,廣東 深圳 518055)

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

1 線性系統(tǒng)加性濾波設(shè)計(jì)方法

1.1 問題陳述

考慮如下隨機(jī)系統(tǒng)

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為離散時(shí)間;A(k)為系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣;B(k+1)為觀測(cè)矩陣;D(k)和C(k+1)為系數(shù)矩陣;x(k)∈Rn為狀態(tài)向量;y(k+1)∈Rm為觀測(cè)向量;d(k)和r(k+1)為系統(tǒng)中的不確定因子。

假設(shè)系統(tǒng)建模誤差w(k)和v(k+1)為不相關(guān)的白噪聲序列,滿足

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積。

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

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

觀測(cè)模型y(k+1)沒有關(guān)于加性干擾d(k)的顯示信息,而d(k)的變化會(huì)反映在狀態(tài)x(k+1)中,因此,以隨機(jī)游走的形式建立一步滯后的動(dòng)態(tài)模型,可得

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

(5)

為了對(duì)加性干擾d(k)精確估計(jì),觀測(cè)模型簡(jiǎn)化為

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)

假設(shè)如下統(tǒng)計(jì)特性

(7)

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

(8)

根據(jù)標(biāo)準(zhǔn)卡爾曼濾波器設(shè)計(jì)原理,加性干擾d(k)的濾波器設(shè)計(jì)可以分為時(shí)間更新和測(cè)量更新兩個(gè)步驟,具體如下。

步驟一:時(shí)間更新

通過式(5),得到關(guān)于d(k)的預(yù)測(cè)值和預(yù)測(cè)誤差協(xié)方差矩陣

(9)

(10)

步驟二:測(cè)量更新

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

(11)

(12)

則d(k)的濾波器設(shè)計(jì)為

(13)

相應(yīng)的估計(jì)誤差為

(14)

結(jié)合正交原理得到濾波器中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)

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

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

(16)

至此,狀態(tài)模型中的加性干擾d(k)的濾波器設(shè)計(jì)就完成了,改寫測(cè)量模型,使模型中出現(xiàn)d(k)的顯示信息,本質(zhì)上是原始測(cè)量模型的等價(jià)改寫。

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

本節(jié)的目標(biāo)是設(shè)計(jì)觀測(cè)模型中加性干擾r(k+1)的Kalman濾波器,同樣以隨機(jī)游走的形式建立動(dòng)態(tài)模型即

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

(17)

假設(shè)如下統(tǒng)計(jì)特性

(18)

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

根據(jù)標(biāo)準(zhǔn)卡爾曼濾波器設(shè)計(jì)原理,加性干擾r(k+1)的濾波器設(shè)計(jì)過程分兩個(gè)步驟,具體如下。

步驟一:時(shí)間更新

基于r(k+1)的動(dòng)態(tài)模型,得到r(k+1)的預(yù)測(cè)值及預(yù)測(cè)誤差協(xié)方差矩陣

(19)

步驟二:測(cè)量更新

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

(20)

(21)

則r(k+1)的濾波器設(shè)計(jì)為

(22)

以及相應(yīng)的估計(jì)誤差為

(23)

最后,利用正交原理再結(jié)合估計(jì)誤差以及觀測(cè)模型可得增益矩陣,并根據(jù)式(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)

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

1.4 求解系統(tǒng)狀態(tài)變量x(k+1)估計(jì)的Kalman濾波器

本節(jié)的主要目的是在已獲得的加性干擾d(k)和r(k+1)估計(jì)值的基礎(chǔ)上,基于原始狀態(tài)模型和測(cè)量模型,求解系統(tǒng)狀態(tài)狀態(tài)變量x(k+1)的濾波器,即

(26)

根據(jù)標(biāo)準(zhǔn)卡爾曼濾波器設(shè)計(jì)原理,濾波器的設(shè)計(jì)過程如下。

步驟一:時(shí)間更新

根據(jù)式(1)可得x(k+1)的狀態(tài)預(yù)測(cè)值、狀態(tài)預(yù)測(cè)誤差及預(yù)測(cè)估計(jì)值誤差協(xié)方差矩陣,分別為

(27)

(28)

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

(29)

步驟二:測(cè)量更新

根據(jù)狀態(tài)預(yù)測(cè)值和觀測(cè)模型可得測(cè)量預(yù)測(cè)值和測(cè)量預(yù)測(cè)誤差,分別為

(30)

(31)

則狀態(tài)變量x(k+1)的濾波器設(shè)計(jì)為

(32)

以及相應(yīng)的狀態(tài)估計(jì)誤差

(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)

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

(35)

2 誤差分析

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

3 仿真驗(yàn)證

考慮如下隨機(jī)系統(tǒng)

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為各變量的估計(jì)曲線以及各變量的估計(jì)誤差曲線,由式(15)～式(18)、式(22)～式(29)、式(31)～式(35)得出。表1為各變量的估計(jì)誤差求和平均后的結(jié)果。

圖1 狀態(tài)x的估計(jì)值

圖3 參數(shù)a的估計(jì)值

圖5 參數(shù)r的估計(jì)值

表1 估計(jì)誤差平均值

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

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

4 結(jié)語

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