程 然，賀豐收，繆禮鋒

(中國航空工業(yè)集團(tuán)公司雷華電子技術(shù)研究所，江蘇無錫 214063)

1 引言

機(jī)載多平臺多傳感器信息融合技術(shù)能夠充分利用各類數(shù)據(jù)的特點，獲取不同視角、不同傳感器、不同特征的互補(bǔ)數(shù)據(jù)，大大擴(kuò)展了探測區(qū)域的覆蓋面積，提高了信息利用效率，加快了數(shù)據(jù)更新速率，從而能夠顯著提升飛機(jī)編隊對目標(biāo)的跟蹤能力，提供更快捷的航跡起始速度，更短的航跡收斂時間，更精確的航跡跟蹤精度和更穩(wěn)定的航跡連續(xù)性[1-2].

然而在艦載機(jī)編隊飛行等復(fù)雜的目標(biāo)跟蹤場景中，系統(tǒng)誤差的存在會導(dǎo)致航跡關(guān)聯(lián)混亂、融合精度降低等問題[2].導(dǎo)致產(chǎn)生這種現(xiàn)象的原因是不同機(jī)載平臺上的傳感器測量信息都是在各自的測量坐標(biāo)系下獲得的，當(dāng)飛機(jī)編隊對目標(biāo)信息進(jìn)行融合跟蹤時，需要將這些測量信息轉(zhuǎn)換到同一參考坐標(biāo)系中才能進(jìn)行關(guān)聯(lián)濾波.由于受到機(jī)載平臺自身定位誤差、傳感器測量誤差及安裝誤差等諸多系統(tǒng)誤差因素的影響，加之坐標(biāo)轉(zhuǎn)換又不可避免地會放大上述誤差，導(dǎo)致融合中心在公共參考坐標(biāo)系下難以將源于同一目標(biāo)的測量值關(guān)聯(lián)起來，嚴(yán)重影響目標(biāo)航跡的穩(wěn)定性和跟蹤精度[2].因此，系統(tǒng)誤差配準(zhǔn)是多源信息融合的前提，為保證后繼多平臺多傳感器信息融合的精度，有必要開展多平臺系統(tǒng)誤差配準(zhǔn)技術(shù)研究.

隨著對系統(tǒng)誤差配準(zhǔn)研究的不斷深入，一些消除系統(tǒng)誤差影響的方法逐漸發(fā)展起來，主要分為兩類:離線估計補(bǔ)償法和在線聯(lián)合估計法.離線估計補(bǔ)償法首先估計出系統(tǒng)誤差，然后對目標(biāo)量測進(jìn)行修正，并基于修正后的量測估計目標(biāo)的真實狀態(tài)，如最小二乘法[3]、廣義最小二乘法[4]、精確極大似然算法[5]等.此類算法都是假設(shè)系統(tǒng)誤差是緩慢變化的，從而對某一時間段內(nèi)的量測數(shù)據(jù)進(jìn)行批處理，但不能對系統(tǒng)誤差的變化做出實時地調(diào)整，缺乏實時性.在線聯(lián)合估計法考慮到系統(tǒng)誤差和目標(biāo)狀態(tài)的相互作用，將系統(tǒng)誤差作為目標(biāo)狀態(tài)的分量同時對二者進(jìn)行估計，如擴(kuò)維高斯濾波算法[6]等.但此類算法需要精確已知系統(tǒng)誤差的狀態(tài)空間模型，然而在實際工程應(yīng)用中，系統(tǒng)誤差的狀態(tài)空間模型可能是不準(zhǔn)確的或者是未知的，從而惡化擴(kuò)維高斯濾波器的性能，甚至導(dǎo)致濾波器發(fā)散.

期望最大化(expectation maximization，EM)算法是一種求解模型中不可觀測的未知參數(shù)最大似然估計的迭代算法，它能保證被估參數(shù)隨著迭代次數(shù)的増加而逐步收斂到一個局部極大值，在處理不完全數(shù)據(jù)中有著重要應(yīng)用[7-8].Henry等人為解決多傳感器信息融合中數(shù)據(jù)關(guān)聯(lián)與誤差配準(zhǔn)之間相互影響的問題，提出了一種卡爾曼濾波(Kalman filter，KF)和期望最大化算法相結(jié)合的算法EM-KF算法，使其能夠同時對系統(tǒng)狀態(tài)和空間誤差進(jìn)行估計，有效提高了信息融合的精度[9].然而該算法僅適用于線性系統(tǒng)的情況，無法推廣到非線性系統(tǒng).劉德浩等人采用EM與擴(kuò)展卡爾曼濾波(extended Kalman filter，EKF)相結(jié)合的方法，解決了量測維數(shù)不同的異類傳感器系統(tǒng)誤差配準(zhǔn)及航跡融合問題[10].然而該方法在計算非線性函數(shù)的雅可比矩陣時對非線性函數(shù)的一階線性化近似精度偏低，當(dāng)系統(tǒng)具有強(qiáng)非線性時，算法的數(shù)值穩(wěn)定性不佳，導(dǎo)致對系統(tǒng)狀態(tài)及誤差的估計精度嚴(yán)重下降，甚至引起濾波發(fā)散.Stefano等人將組網(wǎng)內(nèi)各傳感器的站址誤差、測量誤差及姿態(tài)誤差全部考慮在內(nèi)，提出了一種基于EM算法的多傳感器系統(tǒng)誤差配準(zhǔn)方法[11].然而該算法僅適用于解決固定平臺多傳感器的誤差配準(zhǔn)問題.

為了解決機(jī)載多平臺同質(zhì)多傳感器的系統(tǒng)誤差配準(zhǔn)問題，并主要考慮同質(zhì)傳感器量測系統(tǒng)誤差對目標(biāo)融合跟蹤精度的影響，本文借鑒最大似然估計算法的思想，提出了一種基于期望最大化與容積卡爾曼平滑器(cubature Kalman smoother，CKS)的機(jī)載多平臺同質(zhì)多傳感器系統(tǒng)誤差配準(zhǔn)算法.該算法將傳感器的量測系統(tǒng)誤差視為系統(tǒng)待估計的未知參數(shù)，構(gòu)建了新的傳感器量測方程.同時引入EM算法框架，在期望步驟(expectation step，E-step)中利用容積卡爾曼濾波器(cubature Kalman filter，CKF)和CKS近似計算完整對數(shù)似然函數(shù)的數(shù)學(xué)期望，在最大化步(maximization step，M-step)中對該數(shù)學(xué)期望進(jìn)行最大化處理，最后通過解析更新反復(fù)迭代的方式獲得各傳感器系統(tǒng)誤差的參數(shù)估計.最后，數(shù)值仿真驗證了本文提出算法的有效性.

2 問題陳述

2.1 EM算法

考慮如下狀態(tài)空間形式的離散非線性系統(tǒng)[12]:

其中:xk∈Rn和zk∈Rm分別表示k時刻系統(tǒng)狀態(tài)向量和量測向量;n和m分別表示系統(tǒng)狀態(tài)維數(shù)和量測維數(shù);fk-1(·)和hk(·)分別表示系統(tǒng)狀態(tài)轉(zhuǎn)移函數(shù)和量測函數(shù);wk∈Rn和vk∈Rm分別表示系統(tǒng)過程噪聲向量和量測噪聲向量，二者分別服從均值為零，方差為Qk和Rk的高斯分布，且互不相關(guān).初始系統(tǒng)狀態(tài)向量x0服從均值為方差為P0|0的高斯分布，且它與wk及vk互不相關(guān).η表示未知的系統(tǒng)誤差向量，本文的問題是如何估計出這個未知的系統(tǒng)誤差向量η，為此本文采用EM算法來解決這個問題.EM算法是一種基于貝葉斯推斷在概率模型中尋找未知參數(shù)向量的近似最大似然估計的迭代算法[13-14]，它通過最大化如下的聯(lián)合對數(shù)似然函數(shù)，從而得到未知參數(shù)向量η的一個近似最大似然解[15-16]，即

根據(jù)貝葉斯定理，聯(lián)合對數(shù)似然函數(shù)可以被分解如下:

對方程(6)等號兩邊同時取關(guān)于pηj(x0:N|z1:N)的數(shù)學(xué)期望，并代入到方程(5)中，可以得到

根據(jù)方程(7)，可以推出Lη(z1:N)和Q(η，ηj)的關(guān)系如下[15-16]:

從方程(9)可以看出，如果存在第j+1步η的一個合適近似ηj+1使得Q(ηj+1，ηj)≥Q(ηj，ηj)，那么就能保證成立.因此，EM算法能保證Lη(z1:N)在每一步的迭代過程中都是非遞減的.基于EM算法的未知參數(shù)向量η估計的具體實施步驟被總結(jié)在算法1中.

算法1基于EM算法的未知參數(shù)向量估計.

步驟1選擇一個未知參數(shù)的初始估計η0;

步驟2E-step:根據(jù)方程(5)計算Q(η，ηj);

步驟3M-step:最優(yōu)化計算:

步驟4重復(fù)步驟2-3直到滿足收斂要求.

2.2 容積卡爾曼濾波器和平滑器

在高斯近似濾波器和平滑器中涉及到求解帶高斯分布的非線性狀態(tài)函數(shù)和量測函數(shù)的多維積分，并且這些多維非線性積分全部集中在高斯近似濾波器中，在高斯近似平滑器中則不需計算多余的多維非線性積分.由于狀態(tài)函數(shù)和量測函數(shù)的非線性一般會導(dǎo)致計算多維積分難以實現(xiàn)，因此通常情況下只能采用數(shù)值近似的方法來獲得次優(yōu)的高斯近似濾波器.

基于不同的數(shù)值近似方法，目前已經(jīng)提出了許多次優(yōu)的高斯近似濾波器，其中主要包括基于無跡變換準(zhǔn)則的無跡卡爾曼濾波器(unscented Kalman filter，UKF)[17]，基于高斯埃爾米特求積準(zhǔn)則的高斯埃爾米特濾波器(Gauss Hermite filter，GHF)[18]和基于三階球徑容積準(zhǔn)則的CKF[19]等.UKF通過無跡變換產(chǎn)生確定性采樣點，并利用這些確定性采樣點經(jīng)非線性變換后的結(jié)果來計算系統(tǒng)狀態(tài)的后驗均值和協(xié)方差[17].然而在高維非線性系統(tǒng)條件下，基于無跡變換準(zhǔn)則所產(chǎn)生的確定性采樣點的權(quán)值極易出現(xiàn)負(fù)值情況，導(dǎo)致其矩積分中引入很大的截斷誤差，數(shù)值穩(wěn)定性不佳，濾波精度嚴(yán)重下降，甚至引起濾波發(fā)散.GHF利用高斯埃爾米特求積公式近似計算預(yù)測及預(yù)測協(xié)方差中的非線性函數(shù)多維積分，把多維積分看成是單重積分的嵌套序列，然后依次對每個變量應(yīng)用單重積分公式[18].然而，GHF受到非常嚴(yán)重的維數(shù)災(zāi)難的困擾，函數(shù)的計算次數(shù)隨著狀態(tài)維數(shù)的增長呈指數(shù)增長.盡管高斯埃爾米特公式的近似精度較高，但維數(shù)災(zāi)難所引起的高昂的計算代價限制了它的應(yīng)用推廣.而基于三階球徑容積準(zhǔn)則的CKF所產(chǎn)生的確定性采樣點具有相同的權(quán)值，且權(quán)值永遠(yuǎn)為正，在一定程度上緩解了濾波過程中由于舍入誤差的累計而導(dǎo)致的濾波發(fā)散和數(shù)值精度降低的問題.雖然近似精度比GHF略有降低，但因其計算簡單，易于實現(xiàn)，逐漸被人們接受，并被廣泛使用.因此，本文采用三階球徑容積準(zhǔn)則來計算這些高斯加權(quán)多維非線性函數(shù)積分，具體實施過程如下:

1) 選取容積采樣點ξ(i):

2) 計算含高斯權(quán)重的非線性函數(shù)積分:

其中:f(x)表示任意的非線性函數(shù)，m和P 分別表示高斯分布的均值和協(xié)方差矩陣，表示協(xié)方差矩陣P 的一個Cholesky分解因子.

3 基于EM-CKS的多傳感器系統(tǒng)誤差配準(zhǔn)方法

3.1 E-step

E-Step利用現(xiàn)有未知參數(shù)的估計值來計算完整數(shù)據(jù)對數(shù)似然函數(shù)的期望:

利用方程(1)-(2)給出的狀態(tài)空間模型的一階馬爾科夫性質(zhì)及貝葉斯定理，可以將聯(lián)合概率密度函數(shù)pη(x0:N，z1:N)分解如下:

其中:pη(xk+1|xk)表示系統(tǒng)狀態(tài)轉(zhuǎn)移概率密度函數(shù)，pη(zk+1|xk+1)表示似然概率密度函數(shù).將方程(14)代入到方程(13)中，可以得到

從方程(15)中可以看到，計算Q(η，ηj)需要用到平滑概率密度函數(shù)和聯(lián)合平滑概率密度函數(shù)為了解決這一問題，本文利用了CKS來近似計算這兩個平滑概率密度函數(shù).可被近似為[16]

根據(jù)方程(1)-(2)給出的狀態(tài)空間模型，狀態(tài)轉(zhuǎn)移概率密度分布函數(shù)pη(xk+1|xk)和似然概率密度分布函數(shù)pη(zk+1|xk+1)可以分別表示如下:

將方程(20)代入方程(18)中，Q1(η，ηj)可以被重新表示為

將方程(21)代入方程(19)中，Q2(η，ηj)可以被重新表示為

令關(guān)于任意概率密度函數(shù)p(x)的數(shù)學(xué)期望表示為

則方程(22)和方程(24)可以整理為

3.2 M-step

本節(jié)將解決方程(10)的最大化問題.定義

其中III表示與狀態(tài)同維數(shù)的單位矩陣，將式(31)帶入式(29)中，則I1(η，ηj)可以重新表示如下:

將hk+1(xk+1，η)在xk+1=和η=ηj處分別進(jìn)行一階泰勒級數(shù)展開，得

其中Hk+1和分別表示了系統(tǒng)非線性量測函數(shù)hk+1(xk+1，η)在xk+1=和η=ηj處的雅克比矩陣，即

則方程(30)可以重新表示如下:

將式(32)(38)分別代入式(27)-(28)中，整理得

利用向量求導(dǎo)公式:

將方程(39)等號兩邊對η=ηj+1求偏導(dǎo)，得

將方程(40)等號兩邊對η=ηj+1求偏導(dǎo)，得

而根據(jù)方程(10)中ηj+1的定義，需要滿足如下條件:

即需要滿足如下方程:

因此，可以推出

其中

本文提出的基于EM-CKS的多傳感器系統(tǒng)誤差配準(zhǔn)方法被詳細(xì)的總結(jié)在算法2中.

算法2基于EM-CKS的多傳感器系統(tǒng)誤差配準(zhǔn)方法.

步驟1選擇一個未知參數(shù)的初始估計η0，迭代次數(shù)M，并設(shè)置初始迭代步數(shù)i=0;

步驟2基于當(dāng)前時刻未知參數(shù)向量的估計值ηj運(yùn)行CKF和CKS，最后得到經(jīng)過平滑后的狀態(tài)估計值

步驟3計算方程(47)中的雅克比矩陣H，并將該雅克比矩陣H和步驟2得到的狀態(tài)平滑均值一同代入到方程(47)中更新ηj+1;

步驟4i+1 →i并檢查迭代終止條件i ＞M.如果這個條件滿足，那么迭代終止;否則迭代繼續(xù)，并返回步驟2.

3.3 收斂性討論

從方程(9)可以看出，若Q(ηj+1，ηj)≥Q(ηj，ηj)，那么

并且等號成立當(dāng)且僅當(dāng)滿足如下條件[15]:

根據(jù)方程(48)-(50)，當(dāng)Lη(z1:N)是一個關(guān)于參數(shù)向量η的凸函數(shù)時，并滿足[15]

那么

非線性系統(tǒng)辨識的收斂性是一個非常重要但非常困難的問題，因為其包含的多維非線性積分一般無法解析求解，從而無法獲得Q(η，ηj)和pη(x0:N|z1:N)的解析解.此外在實際應(yīng)用中，似然函數(shù)Lη(z1:N)有可能是一個多峰函數(shù)，存在多個局部解，而EM算法無法避開這個問題.所以通常情況下是無法分析EM算法的收斂性或給出EM算法的收斂條件的.為此，我們將定性地討論參數(shù)選擇對所提出算法收斂性的影響，并說明所提出算法的適用領(lǐng)域和適用條件.

當(dāng)似然函數(shù)Lη(z1:N)是一個多峰函數(shù)時，初始參數(shù)估計η0對所提出算法的收斂性有顯著的影響.如果η0與最優(yōu)解很接近，那么隨著迭代的進(jìn)行，最終有方程(53)成立.如果η0與最優(yōu)解相差很大，那么隨著迭代的進(jìn)行，ηj可能會偏離最優(yōu)解或者陷入局部解.因此，通過選取與最優(yōu)解相近的初值可以改善提出算法的收斂性.在很多實際應(yīng)用中，用戶對模型中的未知參數(shù)的大致取值范圍都有一定的了解，這種信息可以幫助用戶選擇合適的初值，以更好地應(yīng)用本文所提出的算法.

在本文所提出的算法中，狀態(tài)的濾波概率密度和平滑概率密度都被近似為高斯分布，并且它們的高斯近似精度會隨著噪聲方差的增大而降低，這是因為系統(tǒng)的非線性度會隨著噪聲方差Qk和Rk的增大而增大.狀態(tài)的濾波概率密度和平滑概率密度的高斯近似誤差會增大參數(shù)η的估計誤差.而參數(shù)η的估計誤差反過來又會增大狀態(tài)的濾波概率密度和平滑概率密度的高斯近似誤差，從而使得參數(shù)η的估計誤差隨著迭代的進(jìn)行而積累，并且噪聲方差越大這種積累誤差也越大.因此，本文所提出算法對系統(tǒng)誤差的辨識精度會隨著噪聲方差增大而降低.此外，本文所提出算法只有在具有較低或者適中噪聲方差的非線性系統(tǒng)中才會收斂，而對于具有大噪聲方差的非線性系統(tǒng)，提出的算法可能會出現(xiàn)發(fā)散的情況，這是因為大的噪聲方差可能會導(dǎo)致巨大的累積誤差.

3.4 計算復(fù)雜度分析

算法的計算復(fù)雜度是工程應(yīng)用中必須要考慮的一個問題，因此，本節(jié)將對所提出算法的計算復(fù)雜度進(jìn)行分析.分析算法計算復(fù)雜度最為有效的方法是對其浮點操作數(shù)進(jìn)行準(zhǔn)確統(tǒng)計.一次浮點操作數(shù)定義為兩個浮點數(shù)進(jìn)行一次加、減、乘或除法運(yùn)算，下面給出幾種常用的基本代數(shù)運(yùn)算所需的浮點運(yùn)算數(shù)[20].

1) 矩陣加減.

若A∈Rn×m，B∈Rn×m，計算A±B所需的浮點運(yùn)算數(shù)為nm次;

2) 矩陣相乘.

若A∈Rn×m，B∈Rm×l，計算AB所需的浮點運(yùn)算數(shù)為2mnl-nl次;

3) 矩陣求逆.

若A∈Rn×n，計算A-1所需的浮點運(yùn)算數(shù)為n3次;

4) Cholesky分解.

若A∈Rn×n，計算A的Cholesky因子所需的浮點運(yùn)算數(shù)為次;

由第3.2節(jié)的算法2中可以看出，本文所提出的算法實際上是一個不斷迭代更新的過程，且每次迭代都需要依次經(jīng)過CKF、CKS、雅克比矩陣求解及參數(shù)更新這幾個步驟.下面，本文結(jié)合前面給出的基本代數(shù)運(yùn)算，并參照文獻(xiàn)[20]首先給出CKF和CKS所需的浮點運(yùn)算數(shù)，整理如表1所示.

結(jié)合表1所給出的相關(guān)參數(shù)的運(yùn)算復(fù)雜度，單次遞推CKF算法的計算復(fù)雜度可以表示如下[20]:

單次遞推CKS算法的計算復(fù)雜度可以表示如下:

由式(54)-(55)可以看出，CKF和CKS的計算復(fù)雜度均為O(n3)和O(m3)的量級，其中O(·)表示浮點運(yùn)算的階數(shù)，即與狀態(tài)維數(shù)和量測維數(shù)均成三次冪關(guān)系.最后根據(jù)方程(36)(46)-(47)，本文所提出算法的浮點運(yùn)算總數(shù)可以表示如下:

其中:M表示迭代次數(shù)，Mh表示計算量測函數(shù)需要的浮點運(yùn)算數(shù).由式(56)可以很直觀的得出，本文所提出算法的計算復(fù)雜度主要取決于迭代次數(shù)、狀態(tài)維數(shù)及量測維數(shù).因此，在實際工程應(yīng)用中，迭代次數(shù)的選取需要權(quán)衡工程計算量和所需精度.

表1 單次遞推過程計算量Table 1 Calculation amount of single recursive process

4 仿真

本節(jié)將通過模擬雙機(jī)共同跟蹤同一非合作目標(biāo)的仿真試驗來驗證本文所提出的機(jī)載多平臺多傳感器系統(tǒng)誤差配準(zhǔn)算法的有效性和性能.具體仿真場景設(shè)置為由兩架艦載直升機(jī)平臺共同探測跟蹤一架戰(zhàn)斗機(jī)目標(biāo)，其中每架艦載直升機(jī)各攜載一部三維傳感器.戰(zhàn)斗機(jī)飛行速度較快，飛行海拔高度為2 km，艦載直升機(jī)飛行速度較慢，飛行海拔高度為1 km，所有飛行器均沿經(jīng)線飛行，三維傳感器掃描時間間隔為1 s，大地坐標(biāo)系下艦載直升機(jī)平臺與戰(zhàn)斗機(jī)目標(biāo)真實運(yùn)動信息的具體仿真場景如圖1所示.

各傳感器分別測量得到目標(biāo)的距離、方位和俯仰角三維信息，且它們的測量噪聲均為零均值高斯白噪聲，對應(yīng)的各傳感器量測噪聲協(xié)方差矩陣可以表示為

各傳感器系統(tǒng)誤差均設(shè)置為η=[1000 m，0.5°，0.5°]T.

圖1 仿真環(huán)境Fig.1 Diagrams of simulation environment

為了實施本文所提出的算法，在接下來的雙機(jī)目標(biāo)跟蹤仿真試驗中，本文將使用CKF和CKS來近似計算濾波概率密度和平滑概率密度.此外，為了對比說明本文所提出的算法相比于現(xiàn)有的解決多傳感器系統(tǒng)誤差配準(zhǔn)問題的算法所具有的優(yōu)越性，本文在仿真中重點比較文獻(xiàn)[10]所提出的基于EM-EKF的系統(tǒng)誤差配準(zhǔn)方法和本文所提出的基于EM-CKS的系統(tǒng)誤差配準(zhǔn)方法.為了公平比較，所有濾波器都設(shè)置為相同的初始條件，即都從以真實參數(shù)為中心的對稱區(qū)間中隨機(jī)選取未知參數(shù)的初值η0.且都執(zhí)行50次獨立的蒙特卡洛仿真試驗，每次蒙特卡洛仿真試驗的迭代次數(shù)設(shè)置為100次，使用系統(tǒng)誤差在每次迭代的均方根誤差(root mean square error，RMSE)作為性能評價指標(biāo)

首先來看本文提出算法的有效性.圖2-4給出了各傳感器系統(tǒng)誤差的估計結(jié)果.此外，為了進(jìn)一步展示本文所提出算法的精度，表2中給出了各傳感器系統(tǒng)誤差的最終估計.從圖2-4中，可以看到本文所提出的算法在每一次蒙特卡洛運(yùn)行中都是收斂的，并且這種收斂性獨立于初始的傳感器系統(tǒng)誤差參數(shù)值和系統(tǒng)量測輸出.從圖2-4中，也可以看到所有傳感器系統(tǒng)誤差的參數(shù)估計經(jīng)20步迭代后都收斂了，因而本文所提出的算法具有較快的收斂速度.從表2中還可以看到本文所提出的算法具有良好的估計精度.以上這些仿真結(jié)果說明了本文所提出算法可以用于解決機(jī)載多平臺多傳感器系統(tǒng)誤差的非線性系統(tǒng)辨識問題.

圖2 各傳感器徑向距離系統(tǒng)誤差估計結(jié)果Fig.2 Estimated results of range system errors

圖3 各傳感器方位角系統(tǒng)誤差估計結(jié)果Fig.3 Estimated results of azimuth system errors

圖4 各傳感器俯仰角系統(tǒng)誤差估計結(jié)果Fig.4 Estimated results of elevation system errors

表2 真實的參數(shù)值和估計的參數(shù)值(均值±標(biāo)準(zhǔn)差)Table 2 True parameter values and estimated parameter values(mean values±standard deviations)

然后再來看本文提出算法與現(xiàn)有方法的仿真對比.圖5-6分別展示了提出方法和現(xiàn)有方法對傳感器1的方位角系統(tǒng)誤差RMSE和傳感器2的距離系統(tǒng)誤差估計結(jié)果的RMSE，表3給出了兩種方法對傳感器1方位角系統(tǒng)誤差和傳感器2距離系統(tǒng)誤差估計RMSE均值的統(tǒng)計結(jié)果以及單次蒙特卡洛仿真的運(yùn)行時間.

圖5 方位角系統(tǒng)誤差RMSE估計結(jié)果Fig.5 Estimated results of azimuth system errors RMSE

圖6 距離系統(tǒng)誤差RMSE估計結(jié)果Fig.6 Estimated results of distance system errors RMSE

從圖5-6及表3中可以看出，本文提出的基于EMCKS的系統(tǒng)誤差配準(zhǔn)方法比現(xiàn)有的基于EM-EKF的系統(tǒng)誤差配準(zhǔn)方法具有更小的RMSE，但EM-CKS的運(yùn)行時間要慢于EM-EKF.仿真結(jié)果分析如下:本文提出的基于EM-CKS的系統(tǒng)誤差配準(zhǔn)方法采用的是3階球徑容積準(zhǔn)則的確定性采樣方法來近似計算含高斯分布的多維非線性函數(shù)積分，理論上，它能至少以二階泰勒精度逼近任何非線性系統(tǒng)狀態(tài)的后驗均值和協(xié)方差.而EM-EKF采用的是一階泰勒展開線性化的方法，在計算雅克比矩陣時存在截斷誤差，無法保證較高的估計精度.在計算量方面，由于EM-CKS在每一次迭代過程都需要產(chǎn)生一定數(shù)量的采樣點，并經(jīng)過非線性函數(shù)傳播這些采樣點，這一計算過程會隨著系統(tǒng)狀態(tài)及量測維數(shù)的增加而逐漸變慢，因此，EMCKS的單步運(yùn)行時間會多于EM-EKF.雖然本文提出方法在精度上得到了保證，但其卻以犧牲計算量為代價.在實際工程應(yīng)用中，往往需要權(quán)衡工程計算量和所需精度，因此，本文所提出方法適用于對計算量要求不高但對估計精度要求較高的場合.

表3 RMSE均值比較Table 3 Averaged RMSEs

5 結(jié)論

本文提出了一種新的基于最大似然判據(jù)的機(jī)載多平臺多傳感器系統(tǒng)誤差的非線性系統(tǒng)辨識方法.仿真結(jié)果表明本文所提出的算法適用于解決機(jī)載多平臺多傳感器系統(tǒng)誤差的配準(zhǔn)問題，并且所提出的算法具有較快的收斂速度和良好的估計精度.