欒鑄徵，俞成龍，顧 兵，趙先濤

(中國(guó)船舶重工集團(tuán)公司第723 研究所，江蘇 揚(yáng)州 225101)

0 引言

因?yàn)槟繕?biāo)受航路、動(dòng)力及環(huán)境等因素影響，目標(biāo)總是在做機(jī)動(dòng)運(yùn)動(dòng)，例如反艦導(dǎo)彈末端變軌，高空制導(dǎo)炸彈拋物線運(yùn)動(dòng)時(shí)受到空氣阻力和重力的作用，旋翼無(wú)人機(jī)受人為控制飛飛停停等。這種機(jī)動(dòng)性往往是不可預(yù)測(cè)的，使用單一固定的濾波模型很難準(zhǔn)確跟蹤機(jī)動(dòng)目標(biāo)狀態(tài)，濾波器結(jié)果會(huì)發(fā)散，導(dǎo)致跟蹤失敗。因此由Blom 和Bar-Shalom 提出交互多模型(IMM)算法，采用基于位置、勻速、勻加速、Singer 等多種濾波并存方式，目標(biāo)狀態(tài)估計(jì)是多個(gè)濾波器交互作用的結(jié)果，采用馬爾可夫(Markov)鏈控制模型間的交互，把各個(gè)模型上一時(shí)刻的濾波值進(jìn)行交互作用作為各模型的下一時(shí)刻的輸入，然后分別進(jìn)行濾波，得到的結(jié)果進(jìn)行模型概率加權(quán)輸出作為最終的結(jié)果，效果比單模型的好，從而IMM 算法廣泛應(yīng)用到各個(gè)領(lǐng)域[1-4]。但在常規(guī)IMM中馬爾可夫轉(zhuǎn)移概率矩陣是固定值，并且模型概率是通過(guò)卡爾曼濾波(Kalman)更新過(guò)程中產(chǎn)生的殘差來(lái)更新模型概率，模型概率更新及模型概率轉(zhuǎn)移沒(méi)有結(jié)合當(dāng)前的目標(biāo)狀態(tài)分布。所以本文提出了以模型間似然函數(shù)(Likelihood Function)及多模型貝葉斯后驗(yàn)估計(jì)(Bayesian Estimation)融合思想，采用當(dāng)前模型跟蹤結(jié)果更新模型交互概率和以貝葉斯估計(jì)融合多模型輸出作為目標(biāo)狀態(tài)更新值，與目標(biāo)實(shí)際機(jī)動(dòng)情況更加符合，本文對(duì)強(qiáng)機(jī)動(dòng)目標(biāo)和擾動(dòng)靜態(tài)目標(biāo)進(jìn)行了基于Kalman 濾波器的時(shí)變IMM 模型融合算法(TV-IMM)和常規(guī)IMM 方法(C-IMM)仿真，結(jié)果表明時(shí)變IMM 模型融合算法比常規(guī)IMM 方法更有效。

1 時(shí)變交互多模型融合理論分析

本方法仍然采用IMM 處理架構(gòu)，時(shí)變交互多模型融合濾波算法原理是同時(shí)使用多種濾波器對(duì)應(yīng)多個(gè)運(yùn)動(dòng)模型，基于貝葉斯后驗(yàn)估計(jì)的方法得到目標(biāo)當(dāng)前狀態(tài)的最小均方差估計(jì)。首先根據(jù)模型概率和模型交互馬爾可夫轉(zhuǎn)移概率完成各個(gè)模型之間的輸入交互作用，結(jié)果輸入給各個(gè)濾波器預(yù)測(cè)和更新?tīng)顟B(tài)，用各個(gè)濾波器目標(biāo)狀態(tài)分布求解模型間馬爾可夫轉(zhuǎn)移概率，同時(shí)輸入給多模型融合器，根據(jù)貝葉斯后驗(yàn)概率估計(jì)原理得到目標(biāo)狀態(tài)分布更新，再根據(jù)似然函數(shù)原理更新模型概率，從而完成多模型濾波的閉環(huán)跟蹤。

與常規(guī)IMM 不同點(diǎn)在于馬爾可夫轉(zhuǎn)移概率矩陣計(jì)算方法、目標(biāo)狀態(tài)分布更新算法、模型概率更新算法上有不同。該濾波器實(shí)現(xiàn)流程如圖1 所示，假設(shè)有r 個(gè)卡爾曼濾波器模型同時(shí)用于目標(biāo)跟蹤，模型間滿足相互獨(dú)立的多維高斯分布特性。TV-IMM 算法可總結(jié)為如下5步：模型輸入相互作用、模型濾波輸出、模型輸出融合、模型概率更新、模型間轉(zhuǎn)移概率更新。

圖1 時(shí)變交互多模型融合濾波器結(jié)構(gòu)

1.1 模型輸入交互作用[1-5]

該過(guò)程和IMM 交互多模型一致，就是把各濾波器初始條件混合。利用上一步k-1 時(shí)刻得到的模型概率μj(k-1)和馬爾可夫(Markov)交互概率矩陣交互作用，產(chǎn)生新的模型交互概率，代表了模型間相互影響的程度。交互概率作用在每一個(gè)模型的濾波結(jié)果(k-1|k-1)，Pj(k-1|k-1)，其中j=1，…，r，得到輸入交互，馬爾可夫轉(zhuǎn)移矩陣πij表示：

Markov 矩陣第i 行表示第i 個(gè)模型轉(zhuǎn)換為其他模型的概率。第j 列表示由其他模型轉(zhuǎn)換到第j 個(gè)模型的概率。由于各個(gè)模型對(duì)目標(biāo)跟蹤影響程度不同，用模型概率ui表示模型i的影響程度，則第i 個(gè)模型轉(zhuǎn)化為第j 個(gè)模型的概率修正為：

可見(jiàn)模型概率ui影響模型交互概率uij(0)，模型間轉(zhuǎn)化概率滿足概率空間完備性，即第j 列概率總和為1：

概率空間由修正概率組成，所以模型交互概率由歸一化為：

第j 個(gè)模型協(xié)方差為：

模型輸入交互結(jié)果為r 個(gè)濾波器提供輸入。

1.2 模型濾波輸出[4-6]

濾波器是基于卡爾曼濾波原理，利用觀測(cè)空間得到的結(jié)果更新?tīng)顟B(tài)空間的目標(biāo)信息，是最小均方誤差估計(jì)，方程如下：

其中：j 表示第j 個(gè)模型；xj(k)是k 時(shí)刻系統(tǒng)狀態(tài)變量，是dj×1 維列向量，dj為目標(biāo)狀態(tài)維數(shù)；Zj(k)是k 時(shí)刻的系統(tǒng)量測(cè)變量，是nj×1 維列向量，nj為觀測(cè)向量維數(shù)；φj(k)是狀態(tài)轉(zhuǎn)移矩陣，是dj×dj矩陣；Hj(k)是測(cè)量矩陣，是nj×dj矩陣；Wj(k)是高斯型模型白噪聲，0 均值，協(xié)方差為Qj(k)，是dj×dj矩陣；Vj(k)是高斯型量測(cè)白噪聲，0 均值，協(xié)方差為Rj(k)，是nj×nj矩陣。濾波器輸入是基于混合初始狀態(tài)估計(jì)(k-1|k-1)和協(xié)方差(k-1|k-1)，應(yīng)用卡爾曼濾波計(jì)算k 時(shí)刻模型j的狀態(tài)估計(jì)(k|k)和協(xié)方差Pj(k|k)。

狀態(tài)預(yù)測(cè)：

狀態(tài)更新：

其中kj(k)是增益矩陣，vj(k)是殘差，Sj(k)是殘差協(xié)方差矩陣。從而，模型j的狀態(tài)概率分布可以表示為(dj為模型維數(shù))：

本步驟完成了各濾波器對(duì)目標(biāo)狀態(tài)分布估計(jì)。

1.3 模型結(jié)果融合輸出[6-8]

傳統(tǒng)IMM 算法是用各個(gè)模型濾波結(jié)果與模型概率加權(quán)乘積求和后得到目標(biāo)狀態(tài)分布更新，所以IMM 是一種趨勢(shì)控制正確的近似估計(jì)，而不是統(tǒng)計(jì)意義下的最優(yōu)估計(jì)。本文采用貝葉斯后驗(yàn)估計(jì)融合方法得到最優(yōu)后驗(yàn)概率分布估計(jì)，在第2 步中得到各個(gè)模型的卡爾曼濾波狀態(tài)矢量，每一個(gè)結(jié)果代表了模型后驗(yàn)概率目標(biāo)分布估計(jì)，根據(jù)式(9)得到第j 個(gè)模型的目標(biāo)狀態(tài)后驗(yàn)概率分布為，Pj(k|k))，那么根據(jù)貝葉斯估計(jì)原理，由r個(gè)濾波模型同時(shí)產(chǎn)生的目標(biāo)狀態(tài)后驗(yàn)分布可表示為：P(x(k|k)|x1(k|k)，…，xr(k|k))，由于各個(gè)模型間是相互獨(dú)立多維高斯分布，因此有：

根據(jù)貝葉斯估計(jì)原理，最佳估計(jì)就是后驗(yàn)概率分布的期望，即：

而對(duì)于模型j，由式(10)得到概率表達(dá)式為：

從而結(jié)合式(10)，得到融合后狀態(tài)分布顯式表達(dá)，得到融合期望均值和協(xié)方差如下：

本步驟完成了目標(biāo)狀態(tài)融合，獲得了貝葉斯后驗(yàn)概率分布的最優(yōu)估計(jì)。

1.4 模型概率更新

傳統(tǒng)的IMM 算法中，模型概率更新是通過(guò)殘差及協(xié)方差，計(jì)算殘差的似然度，再結(jié)合模型輸入交互概率來(lái)計(jì)算模型概率，沒(méi)有考慮多模型交互更新后的目標(biāo)分布。本文用多模型融合后的目標(biāo)狀態(tài)分布計(jì)算各個(gè)模型似然函數(shù)，可以消除部分狀態(tài)預(yù)測(cè)誤差和觀測(cè)值誤差影響，提高模型準(zhǔn)確度。在第3 步，多模型融合結(jié)果輸出代表了最優(yōu)貝葉斯后驗(yàn)估計(jì)的目標(biāo)狀態(tài)分布，以各個(gè)模型濾波結(jié)果在融合輸出分布的似然函數(shù)，作為模型概率更新的似然度，來(lái)更新模型概率，根據(jù)似然函數(shù)公式，模型j的似然函數(shù)計(jì)算為：

以u(píng)j表示模型概率，根據(jù)概率空間完備性總和為1，即：

概率空間由似然函數(shù)構(gòu)成，所以概率空間根據(jù)模型似然度歸一化為：

本步驟完成了模型概率更新計(jì)算，為1.1 節(jié)k+1 時(shí)刻計(jì)算模型轉(zhuǎn)移概率提供了輸入。

1.5 模型間轉(zhuǎn)移概率計(jì)算

傳統(tǒng)的IMM 算法中，馬爾可夫概率轉(zhuǎn)移矩陣是固定取值。本文提出了以各個(gè)模型濾波值為中心，各個(gè)模型濾波輸出結(jié)果代表了該處理方法對(duì)目標(biāo)狀態(tài)分布的估計(jì)，模型交互定義為各個(gè)模型處理的結(jié)果在其他模型輸出分布的似然函數(shù)，作為模型轉(zhuǎn)化概率更新值。即模型i的濾波值在模型j的目標(biāo)狀態(tài)分布似然函數(shù)，作為模型i 在模型j的概率分布；模型j 自身的似然度由自身模型計(jì)算，從而模型i 到模型j的似然函數(shù)為：

πij表示模型i 到模型j的轉(zhuǎn)移概率，j 取值1 到r，根據(jù)概率空間完備性概率總和為1。

概率空間由模型似然函數(shù)構(gòu)成，所以模型轉(zhuǎn)化概率由似然度歸一化，從而有：

本步驟得到模型間轉(zhuǎn)移概率，為1.1 節(jié)提供了k+1時(shí)刻時(shí)變Markov 轉(zhuǎn)移概率矩陣。

2 時(shí)變IMM與常規(guī)IMM 方法仿真比較

2.1 仿真場(chǎng)景設(shè)計(jì)

時(shí)變IMM 融合方法與常規(guī)IMM 方法進(jìn)行了蒙特卡羅仿真對(duì)比，模擬兩類不易跟蹤的目標(biāo)場(chǎng)景，強(qiáng)機(jī)動(dòng)目標(biāo)場(chǎng)景和強(qiáng)擾動(dòng)靜態(tài)目標(biāo)場(chǎng)景，進(jìn)行了50 次隨機(jī)航路的仿真。首先根據(jù)Kalman 濾波原理，對(duì)目標(biāo)狀態(tài)空間和目標(biāo)觀測(cè)空間建模。

(1)目標(biāo)狀態(tài)空間建模

機(jī)動(dòng)目標(biāo)的運(yùn)動(dòng)模型可以通過(guò)具有加性加速度高斯噪聲的統(tǒng)計(jì)來(lái)描述，所以本文采用勻速目標(biāo)狀態(tài)疊加加速度噪聲模型來(lái)建模，通過(guò)控制參數(shù)和狀態(tài)初始化，模擬產(chǎn)生強(qiáng)機(jī)動(dòng)目標(biāo)和靜態(tài)擾動(dòng)目標(biāo)。

目標(biāo)狀態(tài)方程為：

w(k)的機(jī)動(dòng)目標(biāo)協(xié)方差矩陣為：

其中a 為加速度標(biāo)準(zhǔn)差，取值為0.2。

擾動(dòng)靜態(tài)目標(biāo)協(xié)方差矩陣為：

其中a 為速度標(biāo)準(zhǔn)差,取值為0.2?！鱰 取值為0.5。

(2)目標(biāo)觀測(cè)空間建模

采用兩種模型建立觀測(cè)空間，模型1 是位置模型，模型2 是常速度模型對(duì)目標(biāo)跟蹤。觀測(cè)目標(biāo)狀態(tài)為：

其中：模型1 觀測(cè)矩陣H1=，模型1 觀測(cè)噪聲V1(k)的協(xié)方差矩陣為R1=b2，其中b 為模型1 觀測(cè)誤差，取值為0.5。

模型2 觀測(cè)矩陣H2=，模型2 觀測(cè)噪聲V2(k)的協(xié)方差矩陣為R2=b2，其中b 為模型2觀測(cè)誤差，取值為0.5。

(3)常規(guī)IMM 算法馬爾可夫轉(zhuǎn)移概率矩陣：

2.2 仿真結(jié)果

2.2.1 強(qiáng)機(jī)動(dòng)目標(biāo)仿真結(jié)果

根據(jù)2.1 節(jié)的建模進(jìn)行仿真，由于加速度較大，因此目標(biāo)始終處于強(qiáng)機(jī)動(dòng)狀態(tài)，單次航路仿真跟蹤曲線見(jiàn)圖2，仿真了目標(biāo)真實(shí)位置、模型1 跟蹤、模型2 跟蹤、時(shí)變IMM 融合跟蹤、常規(guī)IMM 跟蹤，共計(jì)5 條曲線，從圖可見(jiàn)，4種方法都可以跟蹤目標(biāo)。圖3 是航路誤差統(tǒng)計(jì)。

圖2 強(qiáng)機(jī)動(dòng)目標(biāo)航路跟蹤仿真結(jié)果圖

圖3 強(qiáng)機(jī)動(dòng)目標(biāo)誤差跟蹤仿真結(jié)果圖

強(qiáng)機(jī)動(dòng)目標(biāo)50 次隨機(jī)航路蒙特卡羅仿真結(jié)果如圖4 所示，可見(jiàn)時(shí)變IMM 融合方法效果最好，航跡誤差曲線在最下方，而常規(guī)IMM 方法誤差曲線介于兩個(gè)模型之間。航跡誤差統(tǒng)計(jì)分布結(jié)果見(jiàn)表1。

表1 強(qiáng)機(jī)動(dòng)目標(biāo)蒙特卡羅仿真結(jié)果分析表

圖4 強(qiáng)機(jī)動(dòng)目標(biāo)50 次蒙特卡羅誤差仿真結(jié)果

2.2.2 強(qiáng)擾動(dòng)靜態(tài)目標(biāo)仿真結(jié)果

強(qiáng)擾動(dòng)靜態(tài)目標(biāo)一次航路仿真跟蹤曲線見(jiàn)圖5，從圖5 可見(jiàn)，模型2 跟蹤發(fā)散，模型1、時(shí)變IMM 融合和常規(guī)IMM 方法可以跟蹤目標(biāo)。圖6 是航路誤差統(tǒng)計(jì)，可見(jiàn)時(shí)變IMM 融合精度最好，常規(guī)IMM 精度在兩個(gè)單模型之間。

圖5 強(qiáng)擾動(dòng)靜態(tài)目標(biāo)航路跟蹤仿真圖

圖6 強(qiáng)擾動(dòng)靜態(tài)目標(biāo)誤差仿真圖

強(qiáng)擾動(dòng)靜態(tài)目標(biāo)50 次蒙特卡羅仿真精度結(jié)果如圖7所示，時(shí)變IMM 融合方法航跡誤差曲線精度最好，分布在下方。統(tǒng)計(jì)的航跡誤差分布結(jié)果見(jiàn)表2，時(shí)變IMM 融合效果顯著。

表2 強(qiáng)擾動(dòng)靜態(tài)目標(biāo)蒙特卡羅仿真結(jié)果分析表

圖7 強(qiáng)擾動(dòng)靜態(tài)目標(biāo)50 次蒙特卡羅誤差仿真結(jié)果

3 結(jié)論

本文應(yīng)用似然函數(shù)理論實(shí)現(xiàn)了時(shí)變馬爾可夫概率轉(zhuǎn)移矩陣和目標(biāo)模型概率更新方法，應(yīng)用貝葉斯估計(jì)理論實(shí)現(xiàn)了多模型融合，得到目標(biāo)最優(yōu)貝葉斯后驗(yàn)準(zhǔn)確估計(jì)。蒙特卡羅仿真結(jié)果驗(yàn)證表明，傳統(tǒng)IMM 方法可以解決目標(biāo)跟蹤的連續(xù)性問(wèn)題，但是并沒(méi)有提高跟蹤精度，相反時(shí)變IMM 融合方法不但解決了目標(biāo)跟蹤的連續(xù)性問(wèn)題，還提高了跟蹤精度。時(shí)變IMM 融合方法在航跡跟蹤連續(xù)性、航跡誤差和精度方面，統(tǒng)計(jì)結(jié)果都優(yōu)于傳統(tǒng)IMM 方法，原因是時(shí)變IMM 融合方法與目標(biāo)實(shí)際狀態(tài)更加吻合。所以理論和仿真結(jié)果都表明時(shí)變IMM 融合算法能更加準(zhǔn)確及時(shí)地跟蹤目標(biāo)，對(duì)跟蹤復(fù)雜機(jī)動(dòng)目標(biāo)有現(xiàn)實(shí)意義。