王增福 邵 毅 祁登亮 金術(shù)玲

①(西北工業(yè)大學(xué)自動化學(xué)院 西安 710072)

②(西安航天動力試驗技術(shù)研究所 西安 710100)

③(中國電子科技集團(tuán)公司第三十八研究所 合肥 230088)

1 引言

天基探測雷達(dá)具有全天候、全天時的戰(zhàn)略、戰(zhàn)術(shù)預(yù)警能力，且具有不受地球曲率限制、不易受攻擊等優(yōu)勢，在預(yù)警防御系統(tǒng)中具有巨大潛力。本文考慮面向大范圍自主搜索發(fā)現(xiàn)任務(wù)的天基雷達(dá)。設(shè)計過程中，為了實現(xiàn)較好的最小可檢測速度，通常方位向孔徑尺寸較大；而綜合考慮重量、收攏等因素，俯仰向口徑尺寸較小，俯仰向測角能力較弱。因此，單個天基探測雷達(dá)對空中目標(biāo)跟蹤定位時，相比于目標(biāo)徑向距量測及方位角量測誤差，目標(biāo)俯仰角量測誤差較大；僅利用單雷達(dá)徑向距量測、方位角量測、俯仰角量測難以有效估計空中目標(biāo)高度?？罩心繕?biāo)高度信息的缺失嚴(yán)重放大了目標(biāo)定位誤差，同時不利于目標(biāo)屬性判別與威脅估計。天基探測雷達(dá)空中目標(biāo)高度估計研究較為少見。Graves[1]提出了一種基于多路徑干擾模式的天基雷達(dá)空中目標(biāo)高度估計方法，理論上難以適用于低空目標(biāo)。而通過多個天基探測雷達(dá)組網(wǎng)，融合利用多個天基探測雷達(dá)徑向距、方位角量測，可提高空中目標(biāo)高度估計準(zhǔn)確度，對改善目標(biāo)定位精度及威脅估計能力等具有重要意義。

天基探測雷達(dá)目標(biāo)定位涉及多次坐標(biāo)轉(zhuǎn)換，具有較強(qiáng)的非線性，需要采用非線性濾波算法估計目標(biāo)運(yùn)動狀態(tài)。常見的非線性濾波方法有擴(kuò)展卡爾曼濾波(Extended Kalman Filter,EKF)、無跡卡爾曼濾波(Unscented Kalman Filter,UKF)、轉(zhuǎn)換量測的卡爾曼濾波(Converted Measurement Kalman Filter,CMKF)和粒子濾波(Particle Filter,PF)等。Kulikov等人[2]基于全局誤差控制常微分方程(Ordinary Differential Equations,ODE)求解器提出了一種精確的連續(xù)離散EKF方法，在高維雷達(dá)目標(biāo)跟蹤場景下對所提算法進(jìn)行了驗證，數(shù)值結(jié)果表明該方法可用于實際的目標(biāo)跟蹤，靈活度和魯棒性更高。Kulikov等人[3]和Zhang等人[4]采用無跡變換(Unscented Transform,UT)來近似非線性系統(tǒng)的狀態(tài)轉(zhuǎn)移過程和量測方程，結(jié)果表明其濾波精度要優(yōu)于EKF方法。然而，EKF和UKF濾波器忽略了狀態(tài)協(xié)方差矩陣的欠估計問題，導(dǎo)致濾波精度下降。張連仲等人[5]采用一種基于期望最大化的轉(zhuǎn)換量測卡爾曼濾波方法，將球坐標(biāo)系下的雷達(dá)量測值轉(zhuǎn)換到直角坐標(biāo)系下，利用期望最大化方法修正了量測噪聲協(xié)方差，最后利用標(biāo)準(zhǔn)卡爾曼濾波器進(jìn)行濾波，提高了野值干擾下的目標(biāo)跟蹤精度，但該方法需要利用真值計算量測轉(zhuǎn)換誤差矩陣，進(jìn)而導(dǎo)致兼容性問題。Wu等人[6]提出了一種重新加權(quán)的魯棒粒子濾波方法來降低量測野值對目標(biāo)跟蹤精度的影響，結(jié)果表明該方法具有較高的跟蹤精度，但是，由于使用蒙特卡羅采樣，粒子濾波所需計算時間成本過高。Ait-El-Fquih等人[7]提出了一種基于變分貝葉斯(Variational Bayes,VB)的粒子濾波方法，該方法將VB方法與隨機(jī)抽樣技術(shù)相結(jié)合以分割狀態(tài)空間，從而推導(dǎo)出一種新的計算效率高的粒子濾波算法。閆文旭等人[8]在VB框架下，將系統(tǒng)狀態(tài)的后驗分布近似于一個可解的變分分布，通過迭代優(yōu)化證據(jù)下界，并考慮懲罰因子調(diào)整迭代的步長，在非線性基礎(chǔ)上推導(dǎo)出閉合形式的迭代EKF。

近年來，概率圖模型及運(yùn)行其上的消息傳遞方法得到了極大發(fā)展，并被應(yīng)用于目標(biāo)跟蹤定位問題中[9–11]。概率圖模型結(jié)合概率論和圖論，利用圖模型表示變量的聯(lián)合概率分布，利用條件獨立性將高維推斷問題分解，從而構(gòu)造高效的學(xué)習(xí)和推斷算法。常見圖模型包括貝葉斯網(wǎng)絡(luò)、馬爾可夫隨機(jī)場以及因子圖，且它們之間可以相互轉(zhuǎn)化。同時，對于復(fù)雜推斷問題，具體的消息傳遞方法近年來也備受關(guān)注[12–15]，常見算法包括信度傳播(Belief Propagation,BP)、環(huán)信度傳播(Loopy BP,LBP)、變分消息傳遞(Variational Message Passing,VMP)、期望傳播(Expectation Propagation,EP)、非參數(shù)信度傳播(Nonparametric BP,NBP)、粒子信度傳播(Particle BP,PBP)、廣義信度傳播(Generalized BP,GBP)、平均場近似與信度傳播組合的消息傳遞等。

多天基雷達(dá)組網(wǎng)融合跟蹤本質(zhì)上是一個復(fù)雜的推斷問題。考慮到天基平臺計算資源的稀缺性等實際需求，需要準(zhǔn)確、低計算復(fù)雜度的融合跟蹤方法。鑒于因子圖模型對復(fù)雜推斷問題強(qiáng)大的表征能力，本文采用因子圖對多部雷達(dá)目標(biāo)跟蹤問題進(jìn)行建模，首先基于因子圖描述傳統(tǒng)的一致卡爾曼濾波，隨后引入耦合參數(shù)，并考慮到目標(biāo)跟蹤模型的非線性特性，推導(dǎo)了分布式一致擴(kuò)展卡爾曼濾波(Distributed Consistent Extended Kalman Filter,DCEKF)算法。隨后，針對EKF線性化非線性模型可能導(dǎo)致誤差變大的問題，在建立的因子圖的基礎(chǔ)上，利用分布式NBP (DCNBP)算法，有效解決了多雷達(dá)目標(biāo)跟蹤的非線性問題，相較于DCEKF算法跟蹤精度明顯提升。

2 系統(tǒng)建模與問題描述

2.1 天基雷達(dá)目標(biāo)跟蹤坐標(biāo)系轉(zhuǎn)換

天基雷達(dá)提供的目標(biāo)測量(徑向距、方位角和俯仰角)位于雷達(dá)測量坐標(biāo)系，而常見的空中目標(biāo)經(jīng)度、緯度、高度信息在大地坐標(biāo)系中描述。為了量測信息的統(tǒng)一，將大地坐標(biāo)系向雷達(dá)測量坐標(biāo)系轉(zhuǎn)換。首先是大地坐標(biāo)系到地心(第4)赤道坐標(biāo)系(簡稱ECEF坐標(biāo)系)的轉(zhuǎn)換，本文采用近似迭代方法[16]，該方法具有更好的轉(zhuǎn)換效率和精度，并且在地球半徑50 km以內(nèi)，不存在奇點和不收斂情況。下面介紹本文涉及的幾種坐標(biāo)系及坐標(biāo)轉(zhuǎn)換[8]。

如圖1所示，ECEF坐標(biāo)系以地球球心O為原點，OXE軸指向0經(jīng)度方向，OZE軸垂直赤道平面并指向地球自轉(zhuǎn)角速度方向，OXE,OYE與OZE軸構(gòu)成右手法則坐標(biāo)系。為了坐標(biāo)轉(zhuǎn)換需要，同樣在圖1中建立地心(第1)赤道坐標(biāo)系(簡稱ECI，又稱J2000.0坐標(biāo)系)，其以地球球心O為原點，OXJ軸指向春分點，OZJ軸垂直赤道平面并指向地球自轉(zhuǎn)角速度方向，OXJ,OYJ與OZJ軸構(gòu)成右手法則坐標(biāo)系。天基雷達(dá)通常安裝在衛(wèi)星上，故建立衛(wèi)星軌道坐標(biāo)系(簡稱LVLH坐標(biāo)系)，其以衛(wèi)星質(zhì)心o為原點，oxL軸為地心指向衛(wèi)星質(zhì)心的衛(wèi)星矢徑方向，ozL軸為衛(wèi)星運(yùn)動速度方向，oxL,oyL和ozL軸構(gòu)成右手法則坐標(biāo)系，下面介紹從ECEF坐標(biāo)系到雷達(dá)量測坐標(biāo)系的轉(zhuǎn)換步驟，坐標(biāo)轉(zhuǎn)換的流程如圖2所示。

圖1 坐標(biāo)系示意圖Fig.1 Coordinate systems

圖2 坐標(biāo)系轉(zhuǎn)換流程圖Fig.2 Procedure of coordinate systems transformation

為了得到較為精確的坐標(biāo)轉(zhuǎn)換，考慮地球的極移、自轉(zhuǎn)、章動、歲差的影響，從ECEF到ECI的坐標(biāo)變換為

其中，[A],[B],[C] 和 [D] 分別為t0時刻的極移、自轉(zhuǎn)、章動和歲差矩陣，且

其中，xp和yp為地極坐標(biāo)，θG為格林尼治恒星時角，εM,Δε和 Δ?分別為平黃赤交角、交角章動和黃經(jīng)章動，Zp,θp和ξp為歲差參數(shù)，t0為歷元起始時刻，Rx(·),Ry(·)和Rz(·) 分別為繞x軸、y軸和z軸旋轉(zhuǎn)的方向余弦矩陣。

從ECI坐標(biāo)系到衛(wèi)星LVLH坐標(biāo)系的旋轉(zhuǎn)矩陣為

其中，ν為緯度幅角，Ω和i為軌道傾角與升交點赤經(jīng)。

衛(wèi)星本體坐標(biāo)系oxbybzb以衛(wèi)星質(zhì)心為坐標(biāo)原點o，3軸為衛(wèi)星的3個慣性主軸，oyb軸垂直于星箭分離面，ozb軸垂直指向星體對地面，oxb,oyb與ozb構(gòu)成右手直角坐標(biāo)系。理想情況下，衛(wèi)星入軌姿態(tài)角全為 0°，衛(wèi)星LVLH坐標(biāo)系與衛(wèi)星本體坐標(biāo)系重合。衛(wèi)星LVLH坐標(biāo)系到衛(wèi)星本體坐標(biāo)系的轉(zhuǎn)換為

其中，γs,βs和αs分別為衛(wèi)星在本體坐標(biāo)系中姿態(tài)的滾動角、俯仰角和偏航角，M1,M2,M3分別為按照右手法則繞x,y,z軸逆時針旋轉(zhuǎn)一定角度的旋轉(zhuǎn)矩陣，xc,yc和zc為目標(biāo)在衛(wèi)星LVLH坐標(biāo)系下的位置，xb,yb和zb為目標(biāo)在衛(wèi)星本體坐標(biāo)系中的位置。

如圖3所示，天基雷達(dá)天線陣面坐標(biāo)系oxayaza坐標(biāo)原點o在天線中心，oxa軸指向天線橫軸，oya軸指向天線縱軸，oza軸指向天線陣面法向，oxa,oya與oza構(gòu)成右手直角坐標(biāo)系。一般天線安裝在衛(wèi)星的對地面(通常天線陣面坐標(biāo)系oxa的軸與衛(wèi)星本體坐標(biāo)系的oxb軸重合)，主要用于對地任務(wù)，入軌后轉(zhuǎn)動角度滿足對地工作要求，當(dāng)所有陣面姿態(tài)角均為 0°時，天線陣面坐標(biāo)系與衛(wèi)星本體坐標(biāo)系重合。衛(wèi)星本體坐標(biāo)系到天線陣面坐標(biāo)系的轉(zhuǎn)換為

圖3 天基雷達(dá)天線陣面坐標(biāo)系Fig.3 Space-based radar antenna array coordinate system

其中，γc,βc和αc分別為雷達(dá)天線陣面的滾動角、俯仰角和偏航角，xa,ya和za為目標(biāo)在天線陣面坐標(biāo)系下的位置，T1,T2,T3分別為按照右手法則繞x,y和z軸逆時針旋轉(zhuǎn)一定角度的旋轉(zhuǎn)矩陣。

雷達(dá)陣面測量坐標(biāo)系為以天線中心為原點的球坐標(biāo)系，從天線陣面坐標(biāo)系到雷達(dá)量測坐標(biāo)系的轉(zhuǎn)換為

其中，rz,θz和?z分別為雷達(dá)陣面測量坐標(biāo)系下的空中目標(biāo)的徑向距、方位角、俯仰角。

2.2 目標(biāo)運(yùn)動模型

不失一般性，考慮空中弱機(jī)動目標(biāo)，此類目標(biāo)的運(yùn)動可以近似為勻速運(yùn)動，其運(yùn)動模型可以描述為

其中，xk為k時刻目標(biāo)在ECEF坐標(biāo)系下的位置速度矢量，，F(xiàn)為目標(biāo)狀態(tài)轉(zhuǎn)移矩陣

其中，T為采樣間隔，ωk-1為過程噪聲，且假設(shè)ωk-1～N(0,Q)。

2.3 雷達(dá)量測模型

假設(shè)在k時刻共有N部天基雷達(dá)對目標(biāo)進(jìn)行測量，每部雷達(dá)得到的目標(biāo)量測信息為徑向距、方位角。為了方便描述多雷達(dá)組網(wǎng)分布式融合中每個雷達(dá)的局部目標(biāo)跟蹤，重寫目標(biāo)運(yùn)動模型式(7)如下。在第n個雷達(dá)節(jié)點，目標(biāo)的運(yùn)動模型可以近似為

由于俯仰角誤差較大，僅利用雷達(dá)徑向距和方位角量測估計目標(biāo)高度及位置，第n個雷達(dá)在k時刻的量測方程表示為

2.4 基于因子圖的天基雷達(dá)組網(wǎng)目標(biāo)跟蹤建模

因子圖是概率圖模型的一種。一般情況下，高維復(fù)雜隨機(jī)系統(tǒng)的推斷問題可以通過帶有局部交互函數(shù)的概率圖模型來建模，即利用因子分解方式和條件獨立性假設(shè)緊湊表示多個變量的聯(lián)合分布。因子圖模型中包含變量節(jié)點和因子節(jié)點，其中每個變量節(jié)點表示一個隨機(jī)變量，每個因子節(jié)點表示一個局部函數(shù)，連接變量節(jié)點和因子節(jié)點的邊表示該因子是該變量的函數(shù)。因子圖模型結(jié)合消息傳遞方法為解決復(fù)雜環(huán)境下的目標(biāo)跟蹤與多源信息融合問題提供了統(tǒng)一框架，將復(fù)雜高維推斷問題轉(zhuǎn)換為優(yōu)化問題，在估計精度、計算復(fù)雜度和實現(xiàn)的靈活性方面具有顯著優(yōu)勢，能夠為目標(biāo)跟蹤問題提供效果好、效率高和可擴(kuò)展的統(tǒng)一解決方案。為此，在上述空中目標(biāo)與天基雷達(dá)量測模型的基礎(chǔ)上，本文采用因子圖方法對天基雷達(dá)組網(wǎng)目標(biāo)跟蹤進(jìn)行建模。

其中，p(xk |xk-1)為狀態(tài)轉(zhuǎn)移概率分布函數(shù)，p(yk |xk)為似然函數(shù)。式(12)可以用因子圖4表示，其中，fk和hk分別表示k時刻狀態(tài)轉(zhuǎn)移函數(shù)和似然函數(shù)對應(yīng)的因子節(jié)點。

圖4 目標(biāo)狀態(tài)與量測聯(lián)合概率分布函數(shù)的因子圖分解Fig.4 Factor graph decomposition of the joint probability distribution function of target states and measurements

圖5 k 時刻局部因子圖模型Fig.5 A local factor graph model at time k

3 基于因子圖的分布式濾波算法

3.1 分布式一致擴(kuò)展卡爾曼濾波

圖4所示因子圖模型中，k時刻，通過因子圖上的遞歸計算得到狀態(tài)xk的后驗邊緣分布和狀態(tài)xk+1的預(yù)測分布，即

式(13)與式(14)為目標(biāo)全局狀態(tài)變量的更新與預(yù)測。根據(jù)2.3節(jié)的雷達(dá)觀測模型，在局部狀態(tài)變量中，引入一個耦合因子節(jié)點gji來表示全局變量xk的兩個復(fù)制狀態(tài)變量和在相鄰節(jié)點i和j上的關(guān)系，其中，i ∈Dj，Dj表示j的鄰居集合。將gji定義為一個指數(shù)函數(shù)[17]

因此，k時刻局部變量的后驗邊緣分布為

因子節(jié)點gji到變量節(jié)點的消息為

將式(15)和式(21)代入到式(20)中，可以得到迭代方程。初始化后，第l次迭代時的方程為

至此，完成了DCEKF算法推導(dǎo)，具體流程總結(jié)如算法1所示。

3.2 分布式非參數(shù)信度傳播算法

在DCEKF算法中，將量測方程式(10)中的非線性進(jìn)行準(zhǔn)線性化近似，以達(dá)到采用均值和方差表示消息的目的。但線性化近似可能產(chǎn)生較大損失，因此針對因子圖模型中的非線性、非高斯概率分布，結(jié)合信度傳播與粒子濾波，采用分布式非參數(shù)信度傳播算法(DCNBP)計算因子圖中的消息傳遞。在2.4節(jié)基礎(chǔ)上加入耦合函數(shù)后，k時刻的因子圖如圖6所示。

算法 1 天基雷達(dá)組網(wǎng)目標(biāo)跟蹤DCEKF算法Alg.1 DCEKF algorithm for target tracking in space-based radar networks

圖6 k時刻因子圖模型Fig.6 The factor graph at time k

其中，g-1(·)表示從雷達(dá)陣面坐標(biāo)系到ECEF坐標(biāo)系的轉(zhuǎn)換函數(shù)，且?z～U(0,2π)。因此，從因子節(jié)點傳遞至變量節(jié)點的消息為，其粒子化形式為

其中，N(x;xi,s,Pi)表示歸一化的高斯密度函數(shù)，均值和協(xié)方差分別為xi,s和Pi。式(41)中，S個高斯密度函數(shù)選取相同的協(xié)方差

因此，D個輸入消息的乘積可以表示為

易知，D個高斯密度函數(shù)的乘積仍然是一個高斯密度函數(shù)

因此，式(45)所示的高斯密度函數(shù)N(x;)的權(quán)重可以表示為

由式(40)可知，p(x)表示D個高斯混合式的乘積，每個高斯混合又含有S個高斯分量，則p(x)一共含有SD個高斯分量，直接采樣的計算復(fù)雜度為O(SD) 。因此，對D個高斯混合式的乘積采用重要性采樣[18]，引入輔助變量用來表示抽取樣本的高斯分量

且ζ1:D={ζ1,ζ2,...,ζD}。離散隨機(jī)變量ζi=li表示從第i個高斯混合式中的第li個高斯分量中抽取一個樣本，則式(40)可寫成

令重要性采樣的建議性分布為p(ζ1:D)，且

綜上，以線圖拓?fù)湎碌奶旎走_(dá)網(wǎng)絡(luò)為例，基于DCNBP的天基雷達(dá)組網(wǎng)分布式目標(biāo)融合跟蹤結(jié)構(gòu)圖如圖7所示。每個雷達(dá)節(jié)點包含利用本地量測的目標(biāo)跟蹤模塊以及利用其他節(jié)點估計值的融合模塊。初始化后，在k時刻，每部雷達(dá)在目標(biāo)跟蹤模塊基于本地量測進(jìn)行本地跟蹤。本地跟蹤后，每部雷達(dá)再根據(jù)相鄰雷達(dá)傳遞的本地跟蹤值與本地量測進(jìn)行消息融合更新，經(jīng)過多次迭代更新，實現(xiàn)多部雷達(dá)目標(biāo)跟蹤的一致性融合。

算法 2 基于重要性采樣的消息乘積粒子化算法Alg.2 Message product particleization algorithm based on importance sampling

3.3 算法復(fù)雜度分析

本節(jié)通過計算執(zhí)行的基本運(yùn)算，比較分析DCEKF算法與DCNBP算法的計算復(fù)雜度。記目標(biāo)的局部狀態(tài)向量大小為s，局部量測向量大小為M，DCNBP算法中粒子數(shù)為S，消息迭代次數(shù)為L。

對于DCNBP算法，每個節(jié)點的消息均為單獨計算，因此其狀態(tài)向量和量測向量維度為s和L。每個時刻，首先采用式(28)和式(29)生成帶權(quán)值粒子進(jìn)行粒子初始化，生成粒子和權(quán)值的計算復(fù)雜度分別為s2S和s2+sS+s，從而該步每個節(jié)點的計算復(fù)雜度為O(s2S)。然后，采用算法2計算算法3中的所有消息乘積的粒子化形式。對于算法2，記節(jié)點的輸入消息個數(shù)為D，計算所有高斯混合消息的均值與協(xié)方差的復(fù)雜度為sS和s2DS2。根據(jù)設(shè)置輔助變量采樣粒子，式(46)和式(47)的計算復(fù)雜度分別為s3D+s2D+s2和sD+s2+2s。因此，D個輸入消息乘積的粒子化形式的計算復(fù)雜度為DS(s3D+s2D+sD+2s2+2s)。從而，算法2的計算復(fù)雜度為O(s2DS2+s3D2S)。因此，經(jīng)過多次迭代計算每個節(jié)點消息的計算復(fù)雜度為O(s2LDS2+s3LD2S)。最后通過算法2計算每個節(jié)點的后驗分布。因此DCNBP算法的計算復(fù)雜度為O(s2LDS2+s3LD2S)。

算法 3 天基雷達(dá)組網(wǎng)目標(biāo)跟蹤DCNBP算法Alg.3 DCNBP algorithm for target tracking in space-based radar networks

本文選擇3顆天基雷達(dá)測量空中目標(biāo)的徑向距和方位角，單顆雷達(dá)的局部量測向量維度M與消息輸入個數(shù)大小D均較小。因此，DCEKF和DCNBP算法的計算復(fù)雜度分別為O(s3L)和O(s2S2L+s3LS)?？梢钥闯觯珼CNBP算法的粒子數(shù)越多，計算復(fù)雜度相較DCEKF算法越大。

3.4 濾波初始狀態(tài)的選取與俯仰角先驗信息

在實際跟蹤中，濾波器初值的選取會對跟蹤性能造成極大的影響，空中目標(biāo)俯仰角誤差過大是天基雷達(dá)目標(biāo)定位的一個關(guān)鍵問題。

雷達(dá)得到的俯仰角誤差越大，得到的濾波器初值與真實值偏差就越大。本文利用幾何法[19]先求解目標(biāo)高度h，再結(jié)合地球半徑re、徑向距rz和方位角θz估計目標(biāo)俯仰角?z，由雷達(dá)與目標(biāo)空間位置關(guān)系可知

其中，rs為雷達(dá)到地心的距離，(xc,yc,zc)為目標(biāo)在LVLH坐標(biāo)系下的坐標(biāo)。實際天基雷達(dá)天線較多安裝在對地面，即陣面姿態(tài)角中僅滾動角變化較多，為簡化推導(dǎo)，僅假設(shè)滾動角γc非零，則

其中，(xa,ya,za)為目標(biāo)在雷達(dá)陣面坐標(biāo)系下的坐標(biāo)，將式(52)代入到式(51)中，得

整理得

再將式(6)代入到式(54)，整理得

通過上述方法可估計出雷達(dá)量測的俯仰角信息，并估計出目標(biāo)高度作為濾波算法初值。一方面提高了濾波算法初值精度，另一方面對于單部雷達(dá)將俯仰角估計值初值加入到量測方程中，可以部分解決單部雷達(dá)因為俯仰角缺失或誤差太大無法跟蹤目標(biāo)的問題，但目標(biāo)高度估計誤差仍然較大。

4 仿真與分析

本節(jié)通過仿真驗證DCNBP算法的性能，與幾何法、傳統(tǒng)集中式EKF (CEKF)、DCEKF算法、集中式粒子濾波(CPF)算法進(jìn)行對比，基于高度求解結(jié)果對俯仰角進(jìn)行估計并將其加入到量測，仿真驗證高度估計誤差對目標(biāo)定位精度的影響。

4.1 仿真場景設(shè)置

考慮一個空中勻速直線飛行的目標(biāo)，高度為9 km，采用3顆天基探測雷達(dá)對此空中目標(biāo)進(jìn)行跟蹤。3顆衛(wèi)星軌道六根數(shù)，即軌道半長軸a，偏心率e，軌道傾角η，近地點幅角ω，升交點赤經(jīng)Ω,平近點角M參數(shù)設(shè)置如下。衛(wèi)星1：6978.14 km,1.369e–15,20°,0°,0°,50°；衛(wèi)星2：7078.14 km,1.369e–15,20°,0°,0°,60°；衛(wèi)星3：7178.14 km,1.369e–15,20°,0°,0°,70°；3顆衛(wèi)星姿態(tài)角均為0°。3部天基雷達(dá)陣面姿態(tài)角均為滾動角 30°、俯仰角 0°、偏航角 0°。雷達(dá)量測噪聲標(biāo)準(zhǔn)差相同，根據(jù)典型的雷達(dá)與目標(biāo)參數(shù)(脈沖寬度、方位向波束寬度等)取值，設(shè)置為徑向距0.09 km，方位角 0.03°。

設(shè)置迭代次數(shù)L=4 ；系統(tǒng)噪聲協(xié)方差均為Q=diag{10-3,10-3,10-3,10-3/T,10-3/T,10-3/T}，T=15 s；仿真得到目標(biāo)在ECEF坐標(biāo)系下的位置，并將其轉(zhuǎn)換到大地坐標(biāo)系，得到目標(biāo)高度估計。算法性能評估指標(biāo)為空中目標(biāo)高度的均方根誤差，即

其中，N表示蒙特卡羅仿真次數(shù)，取N=100，和hk分別表示k時刻目標(biāo)估計高度和真實高度。在CPF算法中，令粒子數(shù)S=800。

4.2 結(jié)果分析

首先令耦合參數(shù)κ=500，粒子數(shù)S=800。CEKF、DCEKF、DCNBP、幾何法估計以及CPF的目標(biāo)高度RMSE結(jié)果在圖8中給出，從圖中可以看出，幾何法的目標(biāo)高度估計誤差最大，在1.5～3.0 km，這是因為僅采用了量測信息。集中式算法由于能夠利用全局信息，一般情況下，相比于分布式估計方法，具有更高的估計精度。因此，可以看到，傳統(tǒng)CEKF算法和CPF算法的精度最高，目標(biāo)高度估計的RMSE為0.50 km和0.45 km。但實際工作時，集中式算法需要將所有量測信息與狀態(tài)信息發(fā)送到一顆衛(wèi)星上進(jìn)行數(shù)據(jù)處理，對該單星的通信負(fù)載與計算速率要求高，并且如果該星發(fā)生故障后，可能會導(dǎo)致整個系統(tǒng)失效，從而，對于天基雷達(dá)應(yīng)用，包括CEKF,CPF算法在內(nèi)的傳統(tǒng)集中式濾波器魯棒性與靈活性較差。而DCEKF,DCNBP等分布式濾波算法能更好地克服此類問題。圖8中，DCEKF算法、DCNBP算法高度估計RMSE分別為0.85 km和0.55 km。相比于DCEKF算法，DCNBP算法估計精度提升了35.3%，同時曲線波動更小，可見DCNBP算法不僅具有較高的估計精度，也具有較好的靈活性和魯棒性。本文未實現(xiàn)分布式PF算法，理論上，與DCNBP算法相比，分布式PF算法具有更高的估計精度，但計算量也更大；天基平臺計算資源較為稀缺，有必要在目標(biāo)高度估計精度與計算量之間做出權(quán)衡。

圖8 空中目標(biāo)高度估計RMSEFig.8 RMSE of aerial target altitude estimation

DCNBP算法中，耦合參數(shù)κ和采樣粒子數(shù)S為關(guān)鍵參數(shù)，能夠影響算法估計性能。圖9和圖10分別為不同κ和S下DCNBP算法的目標(biāo)高度估計結(jié)果，表1中列出了不同參數(shù)下算法估計的RMSE平均值。可以看出，隨著κ和S增大，DCNBP算法高度估計的RMSE變小，精度也越高。這是由于隨著κ的增加，相鄰變量的相關(guān)性越強(qiáng)，鄰居節(jié)點量測信息利用越充分，估計精度越好。而根據(jù)3.1節(jié)，當(dāng)κ →∞時，局部變量與全局變量分布相同，高度估計精度并不會隨著κ的增加一直提高，其最終會趨于某一定值。同時當(dāng)粒子數(shù)S越多，對非線性分布的近似誤差越小，最終的估計精度也更高，但S的增加也會提高算法的計算復(fù)雜度與通信負(fù)載，在應(yīng)用當(dāng)中需要根據(jù)實際情況進(jìn)行權(quán)衡。

表1 不同參數(shù)下DCNBP算法估計的RMSE平均值(km)Tab.1 Average RMSE obtained by DCNBP with different parameters (km)

圖9 空中目標(biāo)高度估計RMSE(粒子數(shù) S=400)Fig.9 RMSE of aerial target altitude estimation(particle number S=400)

圖10 空中目標(biāo)高度估計RMSE(粒子數(shù) S=800)Fig.10 RMSE of aerial target altitude estimation(particle number S=800)

在天基雷達(dá)跟蹤目標(biāo)過程中，雷達(dá)測量誤差與雷達(dá)到目標(biāo)距離相關(guān)，因此分別取雷達(dá)對空中目標(biāo)量測標(biāo)準(zhǔn)差為：徑向距0.07 km，方位角 0.02°；徑向距0.09 km，方位角 0.03°；徑向距0.11 km，方位角 0.04°。圖11所示為不同雷達(dá)測量誤差下DCNBP算法空中目標(biāo)高度估計結(jié)果。由圖可知，當(dāng)目標(biāo)與雷達(dá)距離越近時，雷達(dá)量測誤差越小，目標(biāo)最終的高度估計結(jié)果越精確。

前文中推導(dǎo)了基于高度估計的雷達(dá)俯仰角計算方法，該計算結(jié)果可以被認(rèn)為是一個雷達(dá)偽測量值，高度估計越準(zhǔn)確，該俯仰角偽量測值誤差就越小。為了進(jìn)一步提高目標(biāo)定位精度，將該偽測量值當(dāng)作量測信息加入到濾波中，得到雷達(dá)徑向距、方位角與俯仰角信息，對空中目標(biāo)進(jìn)行定位估計。圖12為不同高度估計誤差下空中目標(biāo)定位結(jié)果的RMSE?？梢钥闯觯叨裙烙嬚`差越大，俯仰角偽量測誤差越大，目標(biāo)定位結(jié)果的RMSE也就越大；當(dāng)高度估計誤差很大時，可等效為俯仰角信息缺失，此時，目標(biāo)定位RMSE最大；隨著高度估計誤差的減小，俯仰角偽量測誤差逐漸減小，目標(biāo)定位RMSE也逐漸減小，可知對目標(biāo)高度估計越準(zhǔn)確，加入俯仰角估計后的目標(biāo)位置估計精度越高。由此可得，優(yōu)先對目標(biāo)高度進(jìn)行估計進(jìn)而得到雷達(dá)俯仰角偽測量值能夠彌補(bǔ)雷達(dá)在實際情況中俯仰角量測誤差大的缺點，進(jìn)一步提高目標(biāo)定位精度。

圖12 不同高度誤差空中目標(biāo)定位估計RMSEFig.12 RMSE of aerial target position estimation with different altitude errors

5 結(jié)語

本文針對天基探測雷達(dá)空中目標(biāo)定位時面臨的俯仰角信息缺失、量測非線性等問題，提出了一種基于一致性的分布式天基雷達(dá)組網(wǎng)空中目標(biāo)高度估計與定位方法。首先建立了空中目標(biāo)運(yùn)動模型與天基雷達(dá)量測模型，并基于因子圖對傳感器目標(biāo)跟蹤進(jìn)行建模?？紤]到雷達(dá)量測非線性，推導(dǎo)了基于因子圖消息傳遞的DCEKF算法；針對DCEKF對高非線性問題近似誤差大的問題，提出了基于采樣的DCNBP算法。仿真驗證表明，相較于DCEKFF、幾何法，DCNBP算法對目標(biāo)高度估計精度更高。同時，分析了DCNBP算法參數(shù)對估計性能的影響。最后，根據(jù)得到的目標(biāo)高度值估計出目標(biāo)的俯仰角，并將該估計值當(dāng)作偽測量估計目標(biāo)位置。仿真結(jié)果顯示，目標(biāo)高度估計越準(zhǔn)確，加入俯仰角估計后的目標(biāo)位置估計精度越高。

利益沖突所有作者均聲明不存在利益沖突

Conflict of Interests The authors declare that there is no conflict of interests