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

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

0 引言

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

1 時變交互多模型融合理論分析

本方法仍然采用IMM 處理架構(gòu)，時變交互多模型融合濾波算法原理是同時使用多種濾波器對應(yīng)多個運(yùn)動模型，基于貝葉斯后驗(yàn)估計(jì)的方法得到目標(biāo)當(dāng)前狀態(tài)的最小均方差估計(jì)。首先根據(jù)模型概率和模型交互馬爾可夫轉(zhuǎn)移概率完成各個模型之間的輸入交互作用，結(jié)果輸入給各個濾波器預(yù)測和更新狀態(tài)，用各個濾波器目標(biāo)狀態(tài)分布求解模型間馬爾可夫轉(zhuǎn)移概率，同時輸入給多模型融合器，根據(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 個卡爾曼濾波器模型同時用于目標(biāo)跟蹤，模型間滿足相互獨(dú)立的多維高斯分布特性。TV-IMM 算法可總結(jié)為如下5步：模型輸入相互作用、模型濾波輸出、模型輸出融合、模型概率更新、模型間轉(zhuǎn)移概率更新。

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

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

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

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

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

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

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

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

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

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

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

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

狀態(tài)更新：

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

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

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

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

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

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

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

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

1.4 模型概率更新

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

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

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

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

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

傳統(tǒng)的IMM 算法中，馬爾可夫概率轉(zhuǎn)移矩陣是固定取值。本文提出了以各個模型濾波值為中心，各個模型濾波輸出結(jié)果代表了該處理方法對目標(biāo)狀態(tài)分布的估計(jì)，模型交互定義為各個模型處理的結(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時刻時變Markov 轉(zhuǎn)移概率矩陣。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 仿真結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

本文應(yīng)用似然函數(shù)理論實(shí)現(xiàn)了時變馬爾可夫概率轉(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ù)性問題，但是并沒有提高跟蹤精度，相反時變IMM 融合方法不但解決了目標(biāo)跟蹤的連續(xù)性問題，還提高了跟蹤精度。時變IMM 融合方法在航跡跟蹤連續(xù)性、航跡誤差和精度方面，統(tǒng)計(jì)結(jié)果都優(yōu)于傳統(tǒng)IMM 方法，原因是時變IMM 融合方法與目標(biāo)實(shí)際狀態(tài)更加吻合。所以理論和仿真結(jié)果都表明時變IMM 融合算法能更加準(zhǔn)確及時地跟蹤目標(biāo)，對跟蹤復(fù)雜機(jī)動目標(biāo)有現(xiàn)實(shí)意義。