段廣全 孫書利

近年來(lái),多傳感器信息融合技術(shù)得到廣泛關(guān)注.經(jīng)典Kalman濾波需要已知模型參數(shù)和噪聲統(tǒng)計(jì)特性[1].而在實(shí)際應(yīng)用中,系統(tǒng)可能會(huì)出現(xiàn)模型參數(shù)或噪聲統(tǒng)計(jì)特性未知情況.文獻(xiàn)[2?4]研究了帶未知噪聲方差系統(tǒng)的自校正估計(jì)問(wèn)題.文獻(xiàn)[5]和[6]研究了模型參數(shù)和噪聲方差未知系統(tǒng)的自校正融合估計(jì)問(wèn)題.文獻(xiàn)[7]針對(duì)帶未知模型參數(shù)和噪聲統(tǒng)計(jì)特性且噪聲相關(guān)的多傳感器系統(tǒng),分別應(yīng)用RELS算法、Gevers-Wouters算法對(duì)未知模型參數(shù)和噪聲統(tǒng)計(jì)特性進(jìn)行辨識(shí),并提出了自校正融合估值器.文獻(xiàn)[5?7]中在對(duì)未知模型參數(shù)進(jìn)行融合處理時(shí),采用加權(quán)平均方法獲得最終融合辨識(shí)器.該方法沒(méi)有考慮不同傳感器的局部參數(shù)辨識(shí)的差異,不能保證融合后的參數(shù)估計(jì)精度都優(yōu)于每個(gè)局部參數(shù)估計(jì).上述文獻(xiàn)所提出的辨識(shí)和估計(jì)算法都是基于完整的傳感器觀測(cè)數(shù)據(jù),而沒(méi)有考慮數(shù)據(jù)的不完整現(xiàn)象.

在實(shí)際的網(wǎng)絡(luò)化系統(tǒng)或傳感器網(wǎng)絡(luò)中,由于傳感器老化或故障、以及傳感器的觀測(cè)數(shù)據(jù)在通信傳輸過(guò)程中由于帶寬有限,可能會(huì)出現(xiàn)數(shù)據(jù)丟失、衰減、延遲等問(wèn)題,使得估值器所收到的傳感器數(shù)據(jù)具有不完整現(xiàn)象[8].文獻(xiàn)[9]研究了帶多丟包系統(tǒng)的最優(yōu)估計(jì)問(wèn)題.文獻(xiàn)[10]將文獻(xiàn)[9]的單傳感器系統(tǒng)推廣到了多傳感器系統(tǒng),對(duì)帶有不同丟包率的多傳感器系統(tǒng),提出了集中式和分布式融合估值器.文獻(xiàn)[11?14]考慮數(shù)據(jù)在傳輸過(guò)程中存在的丟失和延遲現(xiàn)象,設(shè)計(jì)了相應(yīng)的最優(yōu)估值器.其中文獻(xiàn)[13]在文獻(xiàn)[12]的基礎(chǔ)上研究了帶隨機(jī)乘性噪聲、多丟包和滯后系統(tǒng)的最優(yōu)估計(jì)問(wèn)題.文獻(xiàn)[14]考慮了過(guò)程噪聲和觀測(cè)噪聲具有一步自相關(guān)和互相關(guān)的情況.文獻(xiàn)[15]和[16]對(duì)帶有數(shù)據(jù)包丟失的多傳感器系統(tǒng),應(yīng)用協(xié)方差信息方法設(shè)計(jì)了分布式融合估值器.在上述文獻(xiàn)[9?17]中數(shù)據(jù)丟失現(xiàn)象均由一組滿足伯努利分布的隨機(jī)變量所描述,這種數(shù)據(jù)丟失可以看作是觀測(cè)數(shù)據(jù)發(fā)生衰減現(xiàn)象的一種特殊情況.文獻(xiàn)[18]和[19]考慮了帶隨機(jī)參數(shù)矩陣、相關(guān)噪聲和衰減觀測(cè)系統(tǒng)的最優(yōu)和次優(yōu)估計(jì)問(wèn)題.文獻(xiàn)[20]中研究了帶衰減觀測(cè)系統(tǒng)的Kalman 濾波估計(jì)問(wèn)題,同時(shí)分析了誤差協(xié)方差的有界性和穩(wěn)態(tài)特性.文獻(xiàn)[21?22]對(duì)帶衰減觀測(cè)的隨機(jī)不確定系統(tǒng),分別設(shè)計(jì)了多傳感器分布式和序貫融合估值器.上述文獻(xiàn)都是在假設(shè)觀測(cè)丟失率或衰減率已知的情況下得到的最優(yōu)估計(jì)結(jié)果.而沒(méi)有考慮觀測(cè)丟失率或衰減率未知的自校正估計(jì)問(wèn)題.目前,有關(guān)同時(shí)考慮帶有未知模型參數(shù)和未知衰減觀測(cè)率多傳感器系統(tǒng)的自校正融合估計(jì)問(wèn)題的報(bào)導(dǎo)甚少.

基于以上文獻(xiàn)分析,本文將對(duì)帶有未知模型參數(shù)和衰減觀測(cè)率的多傳感器隨機(jī)系統(tǒng),應(yīng)用相關(guān)函數(shù)和遞推增廣最小二乘算法,分別在線辨識(shí)衰減觀測(cè)的數(shù)學(xué)期望、方差和模型參數(shù).應(yīng)用線性無(wú)偏最小方差估計(jì)準(zhǔn)則,分別提出了分布式模型參數(shù)融合辨識(shí)器和自校正融合狀態(tài)濾波器,并分析算法的收斂性.

1 問(wèn)題闡述

考慮帶衰減觀測(cè)的多傳感器隨機(jī)系統(tǒng)

其中狀態(tài)x(t)∈Rn,觀測(cè)yi(t)∈R,過(guò)程噪聲www(t)∈Rr,觀測(cè)噪聲vi(t)∈R,下標(biāo)i表示第i個(gè)傳感器,L表示傳感器的個(gè)數(shù).{μi(t)}為一組在[0,1]區(qū)間取值的用來(lái)描述第i個(gè)傳感器衰減觀測(cè)現(xiàn)象的標(biāo)量隨機(jī)變量.其中E[μi(t)]=αi,Cov[μi(t)]=E為數(shù)學(xué)期望符號(hào),Cov 為協(xié)方差符號(hào).{μi(t)}與其他隨機(jī)變量不相關(guān).Φ,Γ,hi是適當(dāng)維數(shù)矩陣.

假設(shè)1.www(t)和vi(t)為零均值、方差陣分別為Qw和Qvi的不相關(guān)白噪聲.

假設(shè)2.初值x(0)不相關(guān)于w(t)和vi(t),且E{x(0)}=u0,E [x(0)?u0][x(0)?u0]T=P0.其中T為轉(zhuǎn)置號(hào).

假設(shè)3.Φ為穩(wěn)定矩陣,(Φ,hi)是完全可觀對(duì),(Φ,Γ)是完全可控對(duì).

假設(shè)4.Φ中部分參數(shù)未知,{μi(t)}的數(shù)學(xué)期望αi和方差未知.

問(wèn)題是基于觀測(cè)(yi(1),···,yi(t)),i=1,2,···,L,辨識(shí)Φ中未知參數(shù)、{μi(t)}的數(shù)學(xué)期望αi和方差并求Φ中未知參數(shù)的融合辨識(shí)器和狀態(tài)x(t)的自校正融合濾波器

注1.為了閱讀方便,這里對(duì)后文中經(jīng)常遇到一些術(shù)語(yǔ)給予解釋.局部濾波器,即基于單個(gè)傳感器的觀測(cè)數(shù)據(jù)所獲得的濾波器;最優(yōu)濾波器,即系統(tǒng)的模型參數(shù)以及衰減觀測(cè)期望和方差已知情況下,獲得的線性最小方差意義下的濾波器;自校正濾波器,即系統(tǒng)模型含有未知參數(shù),以及衰減觀測(cè)期望和方差未知情況下,通過(guò)辨識(shí)這些未知參數(shù),然后代入最優(yōu)濾波算法中獲得的濾波器.分布式融合濾波器,即基于各個(gè)傳感器的局部濾波器,應(yīng)用線性無(wú)偏最小方差矩陣加權(quán)融合估計(jì)算法[23]獲得的融合濾波器.

2 分布式最優(yōu)融合濾波

當(dāng)系統(tǒng)模型參數(shù)、{μi(t)}的數(shù)學(xué)期望αi和方差已知時(shí),應(yīng)用線性無(wú)偏最小方差意義下的矩陣加權(quán)融合估計(jì)算法[23]可獲得分布式最優(yōu)融合濾波器.下面給出實(shí)現(xiàn)過(guò)程.

由式(2)可得

其中

可計(jì)算其方差陣為

狀態(tài)二階矩X(t)=E[x(t)xT(t)]可遞推計(jì)算如下:

下面引理1 給出了最優(yōu)局部濾波器算法;引理2給出了互協(xié)方差計(jì)算公式;引理3給出了分布式最優(yōu)加權(quán)融合濾波算法.

引理1[24].在假設(shè)1～3 下,隨機(jī)系統(tǒng)(1)和(3)基于每個(gè)傳感器的觀測(cè)有最優(yōu)局部濾波器

初值為Pij(0|0)=P0.

引理3[23].基于引理1的各局部濾波器和引理2的任意兩個(gè)局部濾波誤差之間的互協(xié)方差陣,分布式最優(yōu)矩陣加權(quán)融合濾波器可計(jì)算如下:

加權(quán)矩陣計(jì)算為

其中e=[In,···,In]T,P(t|t)=[Pij(t|t)]nL×nL是以Pij(t|t)為第(i,j)元素的分塊矩陣.融合濾波器的估計(jì)誤差方差陣計(jì)算為

且有Po(t|t)≤Pi(t|t),i=1,···,L.

3 未知模型參數(shù)融合辨識(shí)器

上一節(jié)我們針對(duì)系統(tǒng)模型精確已知時(shí)給出了分布式最優(yōu)融合估計(jì)算法.而在實(shí)際應(yīng)用中,系統(tǒng)模型可能含有未知參數(shù).當(dāng)Φ中含有未知參數(shù)時(shí),本節(jié)采用RELS算法辨識(shí)未知模型參數(shù),并對(duì)辨識(shí)得到的L組參數(shù)估值進(jìn)行加權(quán)融合,獲得模型參數(shù)的分布式融合辨識(shí)器.下面給出具體實(shí)現(xiàn)過(guò)程.

由式(1)可得

式中q?1為單位滯后算子,即q?1x(t)=x(t?1).將式(16)代入式(3)得

將式(17)進(jìn)一步化簡(jiǎn)得

其中A(q?1)=det(In?q?1Φ),Bi(q?1)=hiadj(In?q?1Φ)q?1Γ,式中det和adj分別表示矩陣行列式和伴隨矩陣.A(q?1)和Bi(q?1)具有如下多項(xiàng)式形式:

其中ak,k=1,···,na和Bik,k=1,···,nBi是多項(xiàng)式系數(shù).nA,nBi分別為A(q?1)和Bi(q?1)的階次.式(18)等號(hào)右側(cè)兩個(gè)滑動(dòng)平均過(guò)程可以等價(jià)為一個(gè)穩(wěn)定的滑動(dòng)平均過(guò)程Di(q?1)εi(t)[24],即

其中εi(t)是零均值且?guī)в形粗肼暦讲瞀?εi的白噪聲,Di(q?1)具有如下多項(xiàng)式形式:

其中dik,k=1,···,nDi是多項(xiàng)式Di(q?1)的系數(shù),nDi是Di(q?1)的階次.

將式(18)重寫為

參數(shù)ak,k=1,···,na;dik,k=1,···,nDi未知.

基于每個(gè)單傳感器的觀測(cè)數(shù)據(jù),應(yīng)用RELS算法[24]可得到局部參數(shù)估計(jì)為

當(dāng)i=j時(shí),P?ii(t)即為局部參數(shù)估計(jì)誤差方差陣P?i(t).εi(t)與εj(t)之間的互協(xié)方差可近似計(jì)算如下:

它可遞推地計(jì)算為

令?A=[a1,···,ana]T,則有?Ai=[Ina,0]?i.于是,我們有參數(shù)?A基于傳感器i的局部估計(jì)和估計(jì)誤差協(xié)方差陣如下:

當(dāng)i=j時(shí),P?Aii(t)即為局部參數(shù)?A的估計(jì)誤差方差陣P?Ai(t).

由式(18)可知,參數(shù)?A=[a1,···,ana]T是Φ中未知參數(shù)的函數(shù).假設(shè)Φ中未知模型參數(shù)組成的列向量為Λ[Φ]∈RnΦ,nΦ≤na,且可由?A唯一確定.設(shè)Λ[Φ]與?A之間滿足如下關(guān)系:

其中f(?A)為關(guān)于?A的線性或非線性函數(shù).

1)如果f(?A)是線性函數(shù),我們將式(35)重寫為

式中S,γ為適當(dāng)維數(shù)的系數(shù)陣.

那么,基于傳感器i的數(shù)據(jù)獲得的Φ中未知模型參數(shù)在t時(shí)刻的局部估計(jì)為

2)如果f(?A)是非線性函數(shù),我們將f(?A)在點(diǎn)處進(jìn)行線性化,有

那么,基于傳感器i的數(shù)據(jù)獲得的Φ中未知模型參數(shù)在t時(shí)刻的局部估計(jì)為

其中i,j=1,···,L.通過(guò)以上算法,基于L個(gè)傳感器的數(shù)據(jù)可獲得Φ中未知參數(shù)在時(shí)刻t處的局部估值i,j=1,···,L.由于對(duì)Φ中未知參數(shù)估計(jì)了L次,因此我們可應(yīng)用線性無(wú)偏最小方差加權(quán)融合估計(jì)算法[23]將它們進(jìn)行融合處理.下面定理1給出了未知參數(shù)分布式融合估計(jì)的結(jié)果.

定理1.基于局部參數(shù)估計(jì)局部參數(shù)估計(jì)誤差方差陣以及參數(shù)估計(jì)誤差互協(xié)方差可得在線性無(wú)偏最小方差意義下的矩陣加權(quán)參數(shù)融合辨識(shí)器如下:

參數(shù)融合加權(quán)矩陣計(jì)算為

參數(shù)融合辨識(shí)器的估計(jì)誤差方差陣計(jì)算為

注2.文獻(xiàn)[5?7,24]在進(jìn)行模型參數(shù)融合辨識(shí)時(shí),將各傳感器辨識(shí)得到的模型參數(shù)采用加權(quán)平均方法進(jìn)行融合處理.該方法不能保證所獲得的參數(shù)融合辨識(shí)器的估計(jì)精度不低于所有的局部參數(shù)估計(jì).而本文采用線性無(wú)偏最小方差分布式矩陣加權(quán)融合算法[23]對(duì)各傳感器辨識(shí)得到的參數(shù)進(jìn)行融合處理.所獲得的參數(shù)融合辨識(shí)器的估計(jì)精度不低于所有的局部參數(shù)估計(jì).因此,本文的矩陣加權(quán)融合的參數(shù)估計(jì)精度高于加權(quán)平均融合的參數(shù)估計(jì)精度.這在后面的仿真研究中也能看到.

4 μμi(t)的數(shù)學(xué)期望與方差辨識(shí)

當(dāng)各傳感器的衰減觀測(cè)率未知時(shí),為了能應(yīng)用第2節(jié)中的算法獲得狀態(tài)估計(jì),我們需要辨識(shí)描述衰減觀測(cè)的隨機(jī)變量{μi(t)}的均值和方差.下面我們采用相關(guān)函數(shù)來(lái)辨識(shí)它們.

由式(2)可計(jì)算零步相關(guān)函數(shù)Ri(t,0)=E[y2i(t)]為

一步相關(guān)函數(shù)Ri(t,1)=E[yi(t)yi(t?1)]為

零階和一階相關(guān)函數(shù)陣Ri(t,r),r=0,1可通過(guò)如下采樣相關(guān)函數(shù)逼近:

最后,由式(47)可求出隨機(jī)變量{μi(t)}的數(shù)學(xué)期望為

將式(49)代入式(46)可得隨機(jī)變量{μi(t)}的方差為

5 自校正濾波器的收斂性分析

下面的引理4給出了DESA方法,應(yīng)用DESA方法可證明自校正融合狀態(tài)濾波器的收斂性.

引理4[24].考慮動(dòng)態(tài)誤差系統(tǒng)

其中t≥0,輸出δ(t)∈Rn,輸入u(t)∈Rn,并且矩陣T(t)∈Rn×n是一致漸近穩(wěn)定的.若u(t)是有界的,則δ(t)是有界的.當(dāng)t→∞時(shí),若u(t)→0,則δ(t)→0.

定理2.在假設(shè)1～4下,自校正局部濾波器收斂于最優(yōu)局部濾波器即

證明.由式(7)可得自校正濾波器為

6 仿真例子

在仿真中假設(shè)未知模型參數(shù)a11=0.6,a12=?0.2,μi(t),i=1,2,3,的概率分布分別為P{μ1(t)=0.3}=0.3,P{μ1(t)=0.5}=0.2,P{μ1(t)=1}=0.5,P{μ2(t)=0.4}=0.4,P{μ2(t)=0.7}=0.3,P{μ2(t)=0.9}=0.3,P{μ3(t)=0.1}=0.2,P{μ3(t)=0.6}=0.6,P{μ3(t)=0.9}=0.2.我們可以計(jì)算μi(t),i=1,2,3的數(shù)學(xué)期望和方差分別為

為了與文獻(xiàn)[5?7,24]中的參數(shù)加權(quán)平均融合算法相比較.如下給出了模型參數(shù)a11和a12的局部估計(jì)誤差方差、加權(quán)平均估計(jì)誤差方差以及分布式加權(quán)融合估計(jì)誤差方差算法.

a)局部參數(shù)估計(jì)和分布式加權(quán)融合估計(jì):

根據(jù)前面的第3節(jié),可知未知模型參數(shù)與?A=[a1,a2]T有如下線性關(guān)系

b)參數(shù)的加權(quán)平均融合估計(jì):

Λ[Φ]的加權(quán)平均估值為

圖1給出了應(yīng)用部分3參數(shù)融合辨識(shí)算法獲得的未知模型參數(shù)融合辨識(shí)結(jié)果.由圖可知隨著時(shí)間的增長(zhǎng),辨識(shí)結(jié)果收斂于真值.

圖2和圖3分別給出了未知模型參數(shù)a11和a12的局部估計(jì)誤差方差、加權(quán)平均估計(jì)誤差方差和分布式融合估計(jì)誤差方差的比較結(jié)果.由圖可知分布式加權(quán)融合辨識(shí)誤差方差小于各局部辨識(shí)誤差方差和加權(quán)平均辨識(shí)誤差方差.圖中Si,i=1,2,3表示第i個(gè)傳感器的局部辨識(shí)的誤差方差,DWF表示分布式加權(quán)融合辨識(shí)的誤差方差,WAEV表示加權(quán)平均融合辨識(shí)的誤差方差.

圖1 Φ中未知參數(shù)估計(jì)Fig.1 Identif ication of parameters of Φ

圖2 a11估計(jì)誤差方差Fig.2 Estimation error variance of a11

圖3 a12估計(jì)誤差方差Fig.3 Estimation error variance of a12

圖4和圖5給出了應(yīng)用部分4中辨識(shí)算法分別對(duì)不同傳感器的隨機(jī)變量{μi(t)},i=1,2,3的數(shù)學(xué)期望和方差進(jìn)行辨識(shí)的結(jié)果.曲線表示辨識(shí)結(jié)果,直線表示相應(yīng)的真值.由圖可知隨著時(shí)間的增長(zhǎng),辨識(shí)結(jié)果收斂于真值.圖6和圖7 給出了自校正融合狀態(tài)濾波器,可見(jiàn)自校正融合估計(jì)具有有效性.

圖4 μi(t)的數(shù)學(xué)期望辨識(shí)Fig.4 Identif ication of Mathematical expectation of μi(t)

圖5 μi(t)的方差辨識(shí)Fig.5 Identif ication of variance of μi(t)

圖6 自校正狀態(tài)分量1融合濾波器Fig.6 The first state component of self-tuning fusion filter

圖8和圖9給出了局部和融合的最優(yōu)與自校正狀態(tài)估計(jì)誤差方差圖.由圖可見(jiàn),各局部自校正誤差方差收斂于局部最優(yōu)誤差方差,自校正融合誤差方差收斂于最優(yōu)融合誤差方差,即自校正濾波器具有漸近最優(yōu)性.而且自校正融合濾波器比各局部自校正濾波器具有更高精度.圖中Si,i=1,2,3表示第i個(gè)傳感器的局部自校正估計(jì)誤差方差,SF表示自校正融合估計(jì)誤差方差,直線表示相應(yīng)的最優(yōu)方差.

圖7 自校正狀態(tài)分量2融合濾波器Fig.7 The second state component of self-tuning fusion filter

圖8 局部、融合最優(yōu)與自校正狀態(tài)分量1的濾波誤差方差Fig.8 Variance of the first state component of local,fusion optimal and self-tuning flters

目前參考文獻(xiàn)[5?7,24]中的自校正濾波算法大都沒(méi)有考慮傳感器的衰減觀測(cè)現(xiàn)象.圖10給出了傳感器存在衰減觀測(cè)而沒(méi)有給予考慮的自校正融合濾波器與本文考慮衰減觀測(cè)的自校正融合濾波器在30次蒙特卡洛實(shí)驗(yàn)下均方誤差跡的比較.可見(jiàn),在傳感器存在衰減觀測(cè)時(shí),本文考慮衰減觀測(cè)的自校正融合濾波器具有更高的精度.

7 結(jié)論

對(duì)帶未知模型參數(shù)和衰減觀測(cè)率的多傳感器隨機(jī)系統(tǒng),應(yīng)用RELS算法和相關(guān)函數(shù)分別對(duì)未知模型參數(shù)、描述衰減觀測(cè)現(xiàn)象的隨機(jī)變量的數(shù)學(xué)期望和方差進(jìn)行在線實(shí)時(shí)辨識(shí),提出了線性無(wú)偏最小方差矩陣加權(quán)融合模型參數(shù)辨識(shí)器.與已有文獻(xiàn)的加權(quán)平均融合模型參數(shù)辨識(shí)算法相比,本文所提出的線性無(wú)偏最小方差矩陣加權(quán)融合參數(shù)辨識(shí)算法具有更高的估計(jì)精度.將實(shí)時(shí)辨識(shí)的模型參數(shù)、數(shù)學(xué)期望和方差代入到最優(yōu)局部和融合狀態(tài)估計(jì)算法中獲得了相應(yīng)的自校正狀態(tài)濾波算法.利用DESA方法證明了自校正狀態(tài)濾波器收斂于最優(yōu)狀態(tài)濾波器.與現(xiàn)有文獻(xiàn)的帶未知模型參數(shù)的自校正估計(jì)算法相比,本文還考慮了傳感器的衰減觀測(cè)現(xiàn)象,并給出了采用相關(guān)函數(shù)辨識(shí)衰減觀測(cè)的數(shù)學(xué)期望和方差的算法.

圖9 局部、融合最優(yōu)與自校正狀態(tài)分量2的濾波誤差方差Fig.9 Variance of the second state component of local,fusion optimal and self-tuning filters

圖10 考慮衰減觀測(cè)與沒(méi)有考慮衰減觀測(cè)自校正融合濾波器的均方誤差的跡Fig.10 Trace of mean square error of the self-tuning fusion flters with/without considering fading measurements