自適應(yīng)平方根無跡濾波的雷達(dá)組網(wǎng)狀態(tài)估計(jì)算法*

趙溫波，丁海龍

（解放軍陸軍軍官學(xué)院，合肥 230031）

自適應(yīng)平方根無跡濾波的雷達(dá)組網(wǎng)狀態(tài)估計(jì)算法*

趙溫波，丁海龍

（解放軍陸軍軍官學(xué)院，合肥 230031）

從最小方差的角度分析了雷達(dá)組網(wǎng)無跡濾波（Unscented Filter，UF）狀態(tài)估計(jì)的統(tǒng)計(jì)學(xué)本質(zhì)，并且針對UF的Cholesky分解遇到非半正定矩陣容易發(fā)散、不準(zhǔn)確濾波初值造成濾波發(fā)散以及異常擾動(dòng)影響濾波效果等問題，提出將自適應(yīng)平方根無跡濾波（Adaptive Square Root Unscented Filter，ASQUF）用于雷達(dá)組網(wǎng)狀態(tài)估計(jì)，結(jié)合合理的濾波初始化策略，提高了UF的工程可用性。仿真驗(yàn)證表明，提出的ASRUF算法用于雷達(dá)組網(wǎng)空域目標(biāo)狀態(tài)估計(jì)時(shí)，初始化平穩(wěn)無波動(dòng)，工程可用性好，狀態(tài)估計(jì)精度高，明顯優(yōu)于擴(kuò)展卡爾曼濾波（Extended Kalman Filter，EKF）算法。

雷達(dá)組網(wǎng)，擴(kuò)展卡爾曼濾波，無跡濾波，自適應(yīng)平方根無跡濾波，初始化

0 引言

雷達(dá)組網(wǎng)系統(tǒng)［1-3］為多源數(shù)據(jù)融合系統(tǒng)。雷達(dá)組網(wǎng)數(shù)據(jù)融合［3-4］，依據(jù)數(shù)據(jù)源類型，包括點(diǎn)跡融合［5］和航跡融合［6］兩種融合方式。其中，點(diǎn)跡融合方式目標(biāo)信息豐富，時(shí)延較小，但計(jì)算需求較大。本文主要針對雷達(dá)組網(wǎng)的點(diǎn)跡融合方式開展研究?？柭鼮V波（Kalman Filter，KF）為線性高斯系統(tǒng)的最優(yōu)估計(jì)［7-8］。而雷達(dá)組網(wǎng)系統(tǒng)中，每部組網(wǎng)雷達(dá)送達(dá)融合中心的測量數(shù)據(jù)，其格式為空域目標(biāo)在組網(wǎng)雷達(dá)極坐標(biāo)條件下的斜距、方位角和高低角（海拔高）測量參量。而目標(biāo)跟蹤過程通常在融合中心直角坐標(biāo)系中完成。顯然，組網(wǎng)雷達(dá)極坐標(biāo)測量數(shù)據(jù)與目標(biāo)跟蹤狀態(tài)呈現(xiàn)非線性關(guān)系，不滿足KF線性化使用要求。一般來說，針對非線性系統(tǒng)狀態(tài)估計(jì)問題，工程上通常使用擴(kuò)展卡爾曼濾波（Extended Kalmanfilter，EKF）［9-10］算法。EKF算法使用局部線性化方法，舍棄了非線性模型的二階以上高階項(xiàng)，必然導(dǎo)致目標(biāo)濾波估計(jì)精度下降。無跡（無味）濾波（Unscented Filter，UF）［11-14］算法基于Unscented變換（Unscented Transform，UT）原理，結(jié)合最小方差估計(jì)，達(dá)成逼近目標(biāo)狀態(tài)的后驗(yàn)概率分布目的。對于高斯噪聲系統(tǒng)，UT變換可以逼近四階矩精度；對于非高斯噪聲系統(tǒng)，UT變換也可以達(dá)到三階矩精度［11］。從狀態(tài)估計(jì)方式來看，UF算法為最小方差估計(jì)，其基本原理為基于估計(jì)誤差方差最小化為優(yōu)化目標(biāo)條件下的最優(yōu)狀態(tài)估計(jì)。

本文從最小方差角度分析了UF的統(tǒng)計(jì)學(xué)本質(zhì)，UF在工程使用時(shí)容易受異常擾動(dòng)、不準(zhǔn)確初值和Cholesky對非半正定矩陣難以分解的影響而導(dǎo)致發(fā)散，本文擬使用自適應(yīng)平方根無跡濾波（Adaptive Square Root Unscented Filter，ASQUF）算法，解決雷達(dá)組網(wǎng)系統(tǒng)對空域目標(biāo)的狀態(tài)估計(jì)問題，ASQUF彌補(bǔ)了UF以上的幾個(gè)不足，能夠滿足工程應(yīng)用需求，推進(jìn)了雷達(dá)組網(wǎng)工程化水平。

2 雷達(dá)組網(wǎng)狀態(tài)線性最小方差估計(jì)

假定某時(shí)刻空域目標(biāo)T的真實(shí)狀態(tài)為X，組網(wǎng)雷達(dá)的相應(yīng)探測點(diǎn)跡為Z，探測點(diǎn)跡Z對目標(biāo)T的狀態(tài)估計(jì)（Z）為最優(yōu)化過程為［7］：

當(dāng)空域目標(biāo)T狀態(tài)噪聲、組網(wǎng)雷達(dá)探測噪聲均服從高斯白噪聲分布時(shí)，最小方差估計(jì)MV（Z）為以測量值為函數(shù)的線性最小方差估計(jì)［7］。狀態(tài)估計(jì)及其協(xié)方差表示如式（2）和式（3）所示。

其中，mX、mZ、PX、PXZ分別為狀態(tài)均值、測量均值、狀態(tài)自協(xié)方差矩陣、狀態(tài)-測量互協(xié)方差矩陣。若能獲取這些參量，即可解算空域目標(biāo)T的線性最小方差估計(jì)。

3 基于無跡變換的雷達(dá)組網(wǎng)狀態(tài)/測量的統(tǒng)計(jì)特性描述

一般來說，很難直接獲取式（2）、式（3）所需空域目標(biāo)狀態(tài)均值mX、測量均值mZ、狀態(tài)自協(xié)方差矩陣PX、狀態(tài)-測量互協(xié)方差矩陣PXZ解析表達(dá)。而通過無跡變換方法求解，為行之有效的解決途徑。

無跡變換為基于隨機(jī)分布逼近原理的非線性變換方法。雷達(dá)組網(wǎng)對空域目標(biāo)狀態(tài)估計(jì)，利用無跡變換方法，獲取線性最小方差估計(jì)所需的4個(gè)統(tǒng)計(jì)參量，進(jìn)而完成狀態(tài)濾波估計(jì)。

根據(jù)文獻(xiàn)［11］可知，k-1時(shí)刻目標(biāo)狀態(tài)的σ點(diǎn)集及其相應(yīng)權(quán)重表達(dá)式如式（4）和式（5）所示。

其中κ用于調(diào)節(jié)σ點(diǎn)集逼近狀態(tài)概率分布的高階矩精度，以減少逼近誤差，可取任意值，但要確保n+κ≠0，當(dāng)狀態(tài)分布是高斯分布時(shí)，取n+κ=3最有效［11］。

當(dāng)狀態(tài)維數(shù)增加時(shí)，按式（4）、式（5）方法采樣的σ點(diǎn)集會(huì)呈現(xiàn)離中心點(diǎn)越來越遠(yuǎn)現(xiàn)象［11］，盡管能確保與真實(shí)分布逼近（匹配）的均值和協(xié)方差，但可能會(huì)導(dǎo)致采樣σ點(diǎn)逼近概率分布時(shí)可能丟失局部區(qū)域信息，產(chǎn)生非半正定協(xié)方差矩陣，使矩陣Cholesky分解計(jì)算困難。為了解決這個(gè)問題，文獻(xiàn)［15］提出比例修正采樣，對式（4）的采樣點(diǎn)集進(jìn)行比例修正（修正過程如式（6）所示）。

比例修正系數(shù)α為小于1的正數(shù)，α取值要確保被逼近分布的方差矩陣為半正定矩陣，σ點(diǎn)集經(jīng)過比例修正后用于計(jì)算均值的權(quán)值如式（7）所示。

根據(jù)文獻(xiàn)［15］可知，為了使所逼近分布的協(xié)方差精確到前四階項(xiàng)，特別引入?yún)⒘喀隆?/p>

已經(jīng)證明，如果目標(biāo)狀態(tài)噪聲服從高斯白噪聲分布，β的最優(yōu)取值為2［15］。將式（5）代入式（7）、式（8），令，得到比例修正采樣k-1時(shí)刻的狀態(tài)σ點(diǎn)集及其相應(yīng)權(quán)值。

其中，W*m是均值加權(quán)所用的權(quán)值，W*c是協(xié)方差加權(quán)所用的權(quán)值，為矩陣（n+λ）Pk-1平方根矩陣的第j行或第j列，通常使用矩陣運(yùn)算的Cholesky分解來獲取。

4 雷達(dá)組網(wǎng)無跡濾波狀態(tài)估計(jì)建模

4.1 無跡濾波狀態(tài)/觀測建模

由于本文主要研究雷達(dá)組網(wǎng)系統(tǒng)對空域弱機(jī)動(dòng)目標(biāo)的跟蹤問題，擬采用無跡變換技術(shù)結(jié)合線性最小方差估計(jì)（稱為無跡濾波），針對組網(wǎng)雷達(dá)噪聲化測量值，實(shí)現(xiàn)對空域目標(biāo)狀態(tài)的最優(yōu)估計(jì)。

假定狀態(tài)噪聲和測量噪聲都為加性高斯白噪聲，空域目標(biāo)做勻速直線運(yùn)動(dòng)，則無跡濾波狀態(tài)方程和測量方程表達(dá)式如式（11）和式（12）所示。

其中，Φk|k-1為狀態(tài)轉(zhuǎn)移矩陣，Γk-1為策動(dòng)噪聲矩陣，Wk-1為系統(tǒng)策動(dòng)白噪聲，令策動(dòng)白噪聲方差為Qk-1，Xk為k時(shí)刻空域目標(biāo)的運(yùn)動(dòng)狀態(tài)（由融合中心慣性坐標(biāo)參量表示），本文只研究空域目標(biāo)在融合中心慣性坐標(biāo)系的位置和速度等狀態(tài)參量，則目標(biāo)狀態(tài)為。Zk為組網(wǎng)雷達(dá)測量點(diǎn)跡，由組網(wǎng)雷達(dá)對空域目標(biāo)的斜距、方位角和俯仰角（部分三坐標(biāo)雷達(dá)輸出海拔高測量參量，可以變換為雷達(dá)對目標(biāo)的俯仰角測量值，因此，本文對此不作針對性分析討論）測量參量構(gòu)成，即，Vk為組網(wǎng)雷達(dá)測量噪聲，其統(tǒng)計(jì)特性為Rk。h（g）為融合中心慣性坐標(biāo)系至組網(wǎng)雷達(dá)極坐標(biāo)系的坐標(biāo)變換矩陣［16］。式（11）和式（12）中的參量見文獻(xiàn)［17］。

4.2 無跡濾波狀態(tài)/觀測預(yù)測建模

4.2.1 狀態(tài)預(yù)測建模

針對k-1時(shí)刻空域目標(biāo)狀態(tài)σ點(diǎn)集，利用濾波器狀態(tài)轉(zhuǎn)移方程，生成目標(biāo)狀態(tài)預(yù)測σ點(diǎn)集，基于統(tǒng)計(jì)學(xué)理論方法，生成狀態(tài)預(yù)測統(tǒng)計(jì)特性。

式（13）為預(yù)測σ點(diǎn)集生成表達(dá)式，式（14）為狀態(tài)預(yù)測均值k|k-1表達(dá)式，式（15）為預(yù)測協(xié)方差陣Pk|k-1生成表達(dá)式。

4.2.2 觀測預(yù)測建模

利用式（13）生成的狀態(tài)預(yù)測σ點(diǎn)集，基于濾波器觀測方程，生成觀測σ點(diǎn)集，利用統(tǒng)計(jì)學(xué)方法，獲取觀測σ點(diǎn)集自身的統(tǒng)計(jì)特性及其與狀態(tài)預(yù)測σ點(diǎn)集的互協(xié)方差矩陣。

4.3 無跡濾波狀態(tài)估計(jì)建模

利用式（1）所示的線性最小方差估計(jì)，基于式（14）、式（15）、式（17）～式（19）獲取的狀態(tài)預(yù)測與觀測預(yù)測的統(tǒng)計(jì)特性，根據(jù)k時(shí)間的測量值Zk，獲取狀態(tài)值估計(jì)k及其協(xié)方差矩陣Pk。

5 雷達(dá)組網(wǎng)自適應(yīng)平方根無跡濾波狀態(tài)估計(jì)

將UF用于雷達(dá)組網(wǎng)狀態(tài)估計(jì)需要對估計(jì)協(xié)方差矩陣Pk進(jìn)行Cholesky分解，如果協(xié)方差矩陣Pk是負(fù)定矩陣，則Cholesky分解容易發(fā)散，有研究表明，通過QR分解和Cholesky因子更新方法，通過計(jì)算狀態(tài)預(yù)測方差、觀測預(yù)測方差的 Cholesky分解因子，用以替代Pk直接參與迭代。如式（22）所示，估計(jì)協(xié)方差矩陣Pk解算需要狀態(tài)預(yù)測協(xié)方差、狀態(tài)觀測預(yù)測互協(xié)方差和觀測預(yù)測方差，因此，Pk也可由和、綜合表示，避免了對負(fù)定Pk進(jìn)行Cholesky分解造成發(fā)散現(xiàn)象［16］。另外，不恰當(dāng)?shù)某踔颠x取甚至?xí)斐蔀V波發(fā)散，即使初值選取合理，但如果采樣的觀測預(yù)測點(diǎn)存在很大的誤差，也會(huì)影響濾波精度，甚至造成濾波發(fā)散。本文通過權(quán)衡觀測預(yù)測信息和觀測信息來設(shè)定自適應(yīng)因子，從而控制異常擾動(dòng)對觀測預(yù)測信息的影響，保證濾波能正常收斂。

如果ET=Q，則對ET進(jìn)行QR再分解可得與ET同維數(shù)的矩陣，矩陣的上三角部分就是Pk=EET的Cholesky分解因子，基于這個(gè)原理，可令：

再對ET進(jìn)行QR分解得到2n×n矩陣X

則n×n的上三角矩陣AX方陣就是EET的Cholesky分解因子，因此，狀態(tài)預(yù)測方差PXXk|k-1可表示為：

同理，令

對觀測預(yù)測方差進(jìn)行自適應(yīng)調(diào)節(jié)

對狀態(tài)預(yù)測與觀測預(yù)測協(xié)方差進(jìn)行自適應(yīng)調(diào)節(jié)，式（19）更新為：

直接用Cholesky分解因子進(jìn)行迭代，則式（20）、式（23）變換為：

其中Sk是估計(jì)誤差方差的Cholesky分解因子

式（39）表示進(jìn)行序列cholupdate，再回到第一步的點(diǎn)集采樣進(jìn)行迭代。

6 自適應(yīng)平方根無跡濾波初始化策略構(gòu)建

由于本文只估計(jì)空域目標(biāo)在融合中心坐標(biāo)系的運(yùn)動(dòng)位置和速度，因此，濾波初始化可采用兩點(diǎn)差分法，來確定ASQUF器的初始狀態(tài)估計(jì)和初始協(xié)方差矩陣。

則基于兩時(shí)刻的測量變換值，濾波器初始的估計(jì)誤差表達(dá)式經(jīng)整理如式（43）所示。

則濾波器初始協(xié)方差矩陣的二階中心矩表達(dá)式如式（44）所示。

其中，具體表達(dá)式參見文獻(xiàn)［17］。

7 仿真驗(yàn)證與結(jié)論分析

為了檢驗(yàn)本文ASQUF算法，用于雷達(dá)組網(wǎng)目標(biāo)狀態(tài)估計(jì)的性能，采用仿真方法進(jìn)行驗(yàn)證。仿真場景想定如下：組網(wǎng)雷達(dá)兩部（雷達(dá)1和雷達(dá)2），其掃描周期T都為10 s，雷達(dá)1的距離測量精度為100 m，方位角測量精度為0.2°，俯仰角測量精度為0.3°，配置位置為［119°30°150］T。雷達(dá)2的距離測量精度為80 m，方位測量精度為0.3°，俯仰角測量精度為0.2°，配置位置為［118°32°50］T。融合中心的配置位置為［118°31°170］T。仿真目標(biāo)航速為300 m/s，飛行航向角為60°，飛行高度7 000 m不變，飛行起點(diǎn)位置為［118°32.5°7 000］T，飛行地球大圓距離為300 km。Rk為k時(shí)刻探測點(diǎn)跡歸屬雷達(dá)的測量協(xié)方差矩陣，狀態(tài)噪聲設(shè)定為σx=0.05、σy=0.05、σz=0.05。

分別選取α=0.1，β=2，κ=-3，為了計(jì)算的簡便性，特別假定兩組網(wǎng)雷達(dá)對目標(biāo)探測時(shí)間間隔均為5 s，兩雷達(dá)間隔探測，則300 km仿真航程每1.5 km進(jìn)行一次采樣，全程共采樣201個(gè)觀測點(diǎn)集。初始狀態(tài)估計(jì)X^0和初始協(xié)方差矩陣P0按式（22）、式（26）解算。為了與本文算法形成對比，對于同一仿真測試數(shù)據(jù)，使用EKF進(jìn)行測試對比分析。文中所有仿真計(jì)算均重復(fù)100次。

圖1 雷達(dá)組網(wǎng)仿真場景

圖2 ASQUF位置估計(jì)、噪聲化探測點(diǎn)跡、位置真值對比圖（近區(qū)）

圖1 為兩部雷達(dá)組網(wǎng)仿真場景設(shè)置。圖2是本文提出的ASQUF位置估計(jì)與噪聲化測量點(diǎn)跡、位置真值的對比圖（近區(qū)放大圖），“+”形虛線表示真值點(diǎn)跡，“.”形點(diǎn)劃線表示噪聲化測量點(diǎn)跡，“*”形實(shí)線表示ASQUF估計(jì)。從圖2可以看出，本文構(gòu)造的初始化策略用于ASQUF算法收斂性很好，沒有造成初始波動(dòng)［18］。由于組網(wǎng)雷達(dá)測量噪聲原因，測量點(diǎn)跡對比真實(shí)點(diǎn)跡，位置波動(dòng)較大，但使用本文的ASQUF算法用于雷達(dá)組網(wǎng)目標(biāo)狀態(tài)估計(jì)時(shí)，估計(jì)值能夠逼近目標(biāo)真值。圖3和圖4為本文ASQUF與EKF的位置估計(jì)針對同一測試數(shù)據(jù)集的比較結(jié)果。可以看出，當(dāng)空域目標(biāo)距離兩組網(wǎng)雷達(dá)較近時(shí)，ASQUF與EKF估計(jì)誤差相差無幾。隨著空域目標(biāo)與組網(wǎng)雷達(dá)的空間距離越來越遠(yuǎn)，EKF估計(jì)誤差越來越大，漸趨發(fā)散。由此可見，EKF用于雷達(dá)組網(wǎng)目標(biāo)狀態(tài)估計(jì)時(shí)，在估計(jì)精度方面，要明顯劣于本文構(gòu)造的ASQUF效果。

圖3 X方向ASQUF位置估計(jì)誤差VsEKF位置估計(jì)誤差（100次仿真結(jié)果統(tǒng)計(jì)）

圖4 Y方向ASQUF位置估計(jì)誤差VsEKF位置估計(jì)誤差（100次仿真結(jié)果統(tǒng)計(jì)）

［1］李雪松，江旭敏.雷達(dá)組網(wǎng)技術(shù)研究［J］.探測與定位，2011（6）：9-14.

［2］韓曉威.雷達(dá)組網(wǎng)顯控技術(shù)研究［J］.現(xiàn)代導(dǎo)航，2011，30（10）：367-370.

［3］金虎，王玉松，劉江波.陸上陣地雷達(dá)組網(wǎng)技術(shù)［J］.電子對抗，2011，28（3）：42-45.

［4］王晨熙，王曉博，朱靖，等.雷達(dá)與AIS信息融合綜述［J］.指揮控制與仿真，2009，31（2）:1-4.

［5］方青.雷達(dá)組網(wǎng)數(shù)據(jù)融合處理中的點(diǎn)跡融合技術(shù)［J］.現(xiàn)代電子，2002，19（4）：5-12.

［6］陳科文，張祖平，龍軍.多源信息融合關(guān)鍵問題、研究進(jìn)展與新動(dòng)向［J］.計(jì)算機(jī)科學(xué)，2013，40（8）:8-13.

［7］秦永元，張洪鉞，汪叔華.卡爾曼濾波與組合導(dǎo)航原理［M］.西安:西北工業(yè)大學(xué)出版社，1998：1-43.

［8］Renato Z，Manoranjan M，Robert H B，et al.Norm-Constrained Kalman Filtering［J］.Acta Astronautica.Journal of Guidance ControlandDynamics，2009，32（5）:1458-1465.

［9］Wang Z D，Liu X H，Liu Y R，et al.An Extended Kalman Filtering Approach to Modeling Nonlinear Dynamic Gene Regulatory Networks via Short Gene Expression Time Series［J］.IEEE/ACM Transactions on Computational Biology and Bioinformatics，2009，6（3）:410-419.

［10］李明鎖，井亮，鄒杰，等.結(jié)合擴(kuò)展卡爾曼濾波的CamShift移動(dòng)目標(biāo)跟蹤算法［J］.電光與控制，2011，18（4）:1-5.

［11］Simon Julier，Jeffrey U，Hugh F D W.A New Method for the Nonlinear Transformation of Means and Covariances in Filters and Estimators［J］.IEEE Transactions on Automatic Control，2000，45（3）:477-482.

［12］Orlando X.Diaz.Analysis and Comparison of Extended and Unscented Kalman Filtering Methods for Spacecraft Attitude Determination［D］.California：Naval Postgraduate School，2010.

［13］Li P H，Zhang T W，Ma B.Unscented Kalman Filter for Visual Curve Tracking［J］.Image and Vision Computing，2004（22）:157-164.

［14］王璐，李春光.基于極大似然準(zhǔn)則和最大期望算法的自適應(yīng)UKF算法［J］.自動(dòng)化學(xué)報(bào)，2012，38（7）：1200-1210.

［15］黃鵬，季新源，邱一文.某型火控雷達(dá)組網(wǎng)誘偏AEM研究［J］.四川兵工學(xué)報(bào)，2013（1）：29-32.

［16］Simon J J.The Scaled Unscented Transformation［C］// Proceedings of the American Control Conference.Jefferson City，2002：4555-4559.

［17］Rudolph V M，Eric A W.The Square-root unscented Kalman Filter for State and Parameter Estimation［J］.Proc. IEEE Ine'l Conf.Acoustics，Speech，and Signal Processing，2001（6）：3461-3464.

［18］趙溫波，都基焱，盧達(dá)，等.組網(wǎng)雷達(dá)測量噪聲在慣性坐標(biāo)系中的統(tǒng)計(jì)特性研究［J］.炮兵學(xué)院學(xué)報(bào)，2010（5）：91-95.

Study on State Estimation Algorithm in Radar Networking Based on Adaptive Square Root Unscented Filter

ZHAO Wen-bo，DING Hai-long
（Army Officer Academy of PLA，Hefei 230031，China）

This paper analyzes the statistical nature of Unscented Filter（UF）used in Radar Networking System （RNS）for state estimation at from the perspective of minimum variance.And for problems that Cholesky decomposition in UF can't deal with non-positive semidefinite matrix，inaccurate initial value may cause filter divergence，and disturdance may interfere with filter result in RNS state estimation，this paper proposes Adaptive Square Root Unscented Filter（ASQUF）.Combining rational covariance matrix modeling in initialization，ASQUF improves engineering usability of UF.The simulating verification demonstrates that ASQUF algorithm proposed in this paper has smooth initialization and high precision of estimation and good engineering usability and better performance than Extended KFA（EKFA）.

radarnetworking，extended Kalman filter，unscented filter，adaptive squareroot unscented filter，initialization

TN95

A

1002-0640（2015）05-0043-06

2014-03-25

2014-05-07

國家自然科學(xué)基金（61273001）；安徽省自然科學(xué)基金資助項(xiàng)目（11040606M130）

趙溫波（1972- ），男，吉林省吉林人，副教授。研究方向：數(shù)據(jù)融合，目標(biāo)跟蹤，統(tǒng)計(jì)信號處理。