天基紅外預(yù)警衛(wèi)星高精度協(xié)同目標(biāo)跟蹤方法

王 宇，馬 成，彭 博，陶久亮，于一帆，趙 博

（北京宇航系統(tǒng)工程研究所，北京市 100076）

天基紅外預(yù)警衛(wèi)星常用于助推段的目標(biāo)跟蹤[1-5]，可被動(dòng)探測(cè)目標(biāo)熱輻射，以多星協(xié)同的方式定位目標(biāo)[6]。隨著對(duì)定位精度的要求越來(lái)越高，傳統(tǒng)假設(shè)天基紅外預(yù)警衛(wèi)星量測(cè)噪聲的統(tǒng)計(jì)特性符合高斯分布，并基于高斯濾波方法[7,8]設(shè)計(jì)跟蹤濾波器的方式已經(jīng)無(wú)法滿足要求，從而吸引了越來(lái)越多的學(xué)者開展高精度的協(xié)同跟蹤方法研究。本文將針對(duì)天基紅外預(yù)警衛(wèi)星測(cè)量噪聲的非高斯特性進(jìn)行研究，建立相比于高斯分布更精確的量測(cè)噪聲模型，并提出相應(yīng)的狀態(tài)估計(jì)方法。

天基紅外預(yù)警相機(jī)的像平面是由離散像素組成的，投影在其上的連續(xù)光信號(hào)會(huì)轉(zhuǎn)化為離散的電信號(hào)。在數(shù)字信號(hào)處理領(lǐng)域，用“量化”來(lái)描述連續(xù)信號(hào)被離散的過程，其一般由微處理器的有限計(jì)算精度或通信的帶寬限制造成[9,10]。Widrow對(duì)量化的統(tǒng)計(jì)理論進(jìn)行了深入的研究[11]，給出了量化器的定義，并對(duì)量化噪聲進(jìn)行了統(tǒng)計(jì)分析。Asmar等人[12]比較了幾種用于垂直狀態(tài)估計(jì)的濾波器，并說(shuō)明了在量化噪聲條件下傳統(tǒng)Kalman濾波器性能較差。Rickard-Karlsson等人將量化噪聲視為基于加性均勻噪聲（AUN）假設(shè)的附加均勻噪聲[13]，采用增大測(cè)量協(xié)方差的方法處理量化問題。Curry等人提出了一種通過計(jì)算以量化測(cè)量為條件的矩進(jìn)行量化測(cè)量的方法[10]，這種方法包含復(fù)雜求積運(yùn)算。為了提升算法效率，Duan[14]等人提出了一種量化測(cè)量下的近似MMSE（Minimum Mean Square Error）估計(jì)數(shù)值算法?？梢?，數(shù)字信號(hào)處理領(lǐng)域?qū)B續(xù)信號(hào)被離散的過程進(jìn)行了較深入的研究，本文將借鑒信號(hào)處理領(lǐng)域“量化”的概念，建立天基紅外預(yù)警衛(wèi)星的量測(cè)模型。

天基紅外預(yù)警衛(wèi)星目標(biāo)跟蹤系統(tǒng)是強(qiáng)非線性系統(tǒng)[15]。在各種非線性系統(tǒng)目標(biāo)跟蹤問題中，基于一階線性化的擴(kuò)展卡爾曼濾波（Extended Kalman Filter,EKF）算法應(yīng)用最廣。一些研究將容積卡爾曼濾波器[16,17]、粒子濾波[18-20]，應(yīng)用于目標(biāo)跟蹤問題。但相對(duì)于傳統(tǒng)EKF方法計(jì)算量較大。

本文主要工作包括：1）基于助推段運(yùn)動(dòng)微分方程，建立形式簡(jiǎn)單且符合飛行器運(yùn)動(dòng)特性的目標(biāo)運(yùn)動(dòng)方程；2）從理論上分析天基紅外預(yù)警衛(wèi)星測(cè)量的特點(diǎn),利用數(shù)字信號(hào)處理中的mid-riser量化器描述離散像平面，建立天基紅外預(yù)警衛(wèi)星的測(cè)量噪聲模型；3）為解決基于天基紅外預(yù)警衛(wèi)星的高精度目標(biāo)跟蹤問題，提出了量化擴(kuò)展卡爾曼濾波（Quantized Extended Kalman Filter, QEKF）方法，該方法基于給定量化測(cè)量值的條件均值估計(jì)來(lái)進(jìn)行測(cè)量更新。

1 助推段目標(biāo)運(yùn)動(dòng)模型建立

在助推階段的目標(biāo)加速度可以表示為：

其中，aN表示凈加速度；aG表示重力加速度。

采用重力轉(zhuǎn)彎模型描述目標(biāo)的助推段運(yùn)動(dòng)[8]。由于助推段飛行攻角較小，在地球固連坐標(biāo)系中，假設(shè)推力和阻力平行于相對(duì)速度矢量，可得到：

定義[x,y,z]為目標(biāo)地球固連坐標(biāo)系的坐標(biāo)，可得：

其中，ωe為地球旋轉(zhuǎn)速度。且：

其中，α為待估計(jì)量，為推導(dǎo)其動(dòng)力學(xué)模型，假設(shè)非重力凈力的大小F（t）?T（t）-D（t）=F（t0）為常值；其中，t0為任意初始時(shí)間。假設(shè)目標(biāo)質(zhì)量m以常速率線性遞減，即m˙=-δm，其中，δm為常值。

對(duì)于t≥t0，有：

其中，β（t）為輔助參數(shù)。

求導(dǎo)可得：

定義狀態(tài)向量x如下：

根據(jù)式(7)-(10)，可整理得：

可見，該模型具備形式簡(jiǎn)單、不需要目標(biāo)質(zhì)量等先驗(yàn)信息的優(yōu)點(diǎn)。

2 天基紅外預(yù)警衛(wèi)星量測(cè)模型建立

在實(shí)際應(yīng)用中，至少有兩顆衛(wèi)星同時(shí)用于目標(biāo)跟蹤。為了表征天基紅外預(yù)警衛(wèi)星量測(cè)噪聲的統(tǒng)計(jì)特性，高斯分布被廣泛使用[6,7]。實(shí)際上，天基紅外預(yù)警衛(wèi)星相機(jī)的像平面是由離散的像素組成的，導(dǎo)致測(cè)量的量化。本節(jié)將分析天基紅外預(yù)警衛(wèi)星測(cè)量特性，建立量測(cè)模型。

2.1 理想量測(cè)方程

假設(shè)相機(jī)坐標(biāo)系與衛(wèi)星體坐標(biāo)系重合，可表示為Ob X bY b Zb；目標(biāo)空間點(diǎn)（x,y,z），相機(jī)透視中心Ob，像點(diǎn)（yi,zi）在同一條線上，見圖1；像平面坐標(biāo)系（單位為毫米）與像素坐標(biāo)系O pY pZp（單位為像素）重合。

圖1 像平面幾何關(guān)系Fig.1 Geometric relationship of image plane

定義k時(shí)刻第i顆衛(wèi)星的理想量測(cè)為。由和組成，計(jì)算方式為：

其中，上標(biāo)i表示第i個(gè)衛(wèi)星；表示轉(zhuǎn)換矩陣的第a行，第b列的元素。

2.2 量化量測(cè)模型

在天基紅外預(yù)警衛(wèi)星的離散像平面上，照射到像素內(nèi)的光線均產(chǎn)生對(duì)應(yīng)于該像素中心的坐標(biāo)輸出，見圖2。目標(biāo)在像素坐標(biāo)系的坐標(biāo)（,）受到量化噪聲的干擾，干擾量最大為0.5像素。

圖2 離散像平面Fig.2 Discrete image plane

從離散圖像平面的輸出，只能確定實(shí)際坐標(biāo)所在的像素，而無(wú)法得到準(zhǔn)確坐標(biāo)。為了對(duì)天基紅外預(yù)警衛(wèi)星的量測(cè)進(jìn)行精確建模，引入mid-riser量化器，將量測(cè)值坐標(biāo)輸出可以看作是mid-riser量化器的輸出。mid-riser量化器的數(shù)學(xué)表達(dá)式如下：

基于mid-riser量化器，目標(biāo)在像素坐標(biāo)系中的理想坐標(biāo)與量化坐標(biāo)之間的關(guān)系可表示為：

其中，wy和wz為弱照射條件、高溫等引起的零均值高斯白噪聲；ey和ez為與像素大小相關(guān)的量化噪聲。

從式(17)可看出，天基紅外預(yù)警衛(wèi)星量測(cè)噪聲由兩部分組成：1) 高斯噪聲wy和wz；2) 量化噪聲ey和ez。在高斯噪聲大于量化噪聲的情況下，可延用傳統(tǒng)EKF濾波方法；當(dāng)量化噪聲更大時(shí)，量測(cè)的非高斯特性會(huì)引起EKF方法精度差，甚至是濾波發(fā)散。

將式(15)代入式(17)，得到衛(wèi)星i量測(cè)方程：

3 量化量測(cè)協(xié)同目標(biāo)跟蹤方法

考慮如下量化量測(cè)非線性離散系統(tǒng)

其中，xk?Rn為狀態(tài)向量；yk?Rm為量化前的量測(cè)；vk和wk為相互獨(dú)立的零均值高斯噪聲，協(xié)方差分別為Qk和Rk；下標(biāo)k代表離散時(shí)間tk；f和h分別為狀態(tài)方程和離散方程；ykq為量化后的量測(cè)向量。需要指出的是，即使量測(cè)方程h為線性方程、wk為高斯噪聲，由于量化噪聲的存在，Q（yk）將是非線性的，ek是非高斯的。

定義

QEKF和EKF的主要區(qū)別在于量測(cè)更新，EKF量測(cè)更新公式為：

新息協(xié)方差Sk：

其中，Hk為雅克比矩陣。

濾波增益Kk：

狀態(tài)及協(xié)方差更新方程為：

其中，

Kk可根據(jù)式(29)計(jì)算。

其中，

整理得：

其中，MSE（x k|,yk）可根據(jù)公式(29)計(jì)算，即：

協(xié)方差：

因此，可得量測(cè)更新如式(31)(32)(35)(37)。

QEKF中的高維高斯函數(shù)積分沒有解析解。本文結(jié)合Genz的轉(zhuǎn)換方法以及偽蒙特卡洛積分方法提出一種數(shù)值求解方法，將積分轉(zhuǎn)換為單位超立方體上的積分，從積分區(qū)間直接取點(diǎn)，可避免維度災(zāi)難。

涉及到的高斯積分為：

式(44)可以進(jìn)一步寫為：

之后，利用偽蒙特卡洛方法對(duì)式(43)進(jìn)行積分，即：

至此，完成基于Genz轉(zhuǎn)換和偽蒙特卡洛的多維高斯積分計(jì)算。對(duì)于式(31)和式(37)中的積分，僅需對(duì)以上積分函數(shù)進(jìn)行簡(jiǎn)單改動(dòng)即可。

4 數(shù)學(xué)仿真

考慮兩顆天基紅外預(yù)警衛(wèi)星協(xié)同跟蹤助推段目標(biāo)的場(chǎng)景，利用本文提出的協(xié)同目標(biāo)跟蹤方法對(duì)目標(biāo)進(jìn)行跟蹤。

4.1 仿真場(chǎng)景

目標(biāo)助推段彈道如圖3。

圖3 目標(biāo)助推段彈道Fig.3 Boost-phase Trajectory

仿真參數(shù)：兩顆衛(wèi)星分別位于105°E、145°E地球同步軌道；兩衛(wèi)星量測(cè)噪聲為：高斯噪聲σ1=0.1，量化噪聲Δ=1；初始位置誤差[5000 5000 5000] m；初始速度誤差[500 500 500] m/s；像素大小εy×εz=50 μm×50 μm；焦距長(zhǎng)度f(wàn)=1 m。

4.2 仿真結(jié)果

將本文提出的QEKF和傳統(tǒng)EKF方法（EKF-Q），基于AUN假設(shè)的改進(jìn)EKF方法（EKF-AUN）[13]，近似最小均方差濾波方法（NMMSE）[14]，以及純高斯量測(cè)下的EKF方法進(jìn)行對(duì)比。定義均方根誤差（RMSE）衡量算法精度為：

其中，M表示蒙特卡洛打靶次數(shù)，表示k時(shí)刻第i次狀態(tài)估計(jì)結(jié)果；xk為k時(shí)刻真實(shí)狀態(tài)。

圖4-5給出了QEKF單次仿真跟蹤結(jié)果?？梢?，QEKF可以很快收斂，且位置、速度估計(jì)誤差在3σ邊界以內(nèi)。

圖4 QEKF位置估計(jì)誤差Fig.4 Position estimation error of QEKF

圖5 QEKF速度估計(jì)誤差Fig.5 Velocity estimation error of QEKF

開展500次蒙特卡洛仿真，得到如圖6和圖7的不同算法位置、速度估計(jì)誤差對(duì)比結(jié)果。將不同算法10 s～60 s位置、速度的RMSE取平均，得到結(jié)果如表1。可以看出，EKF-Q的位置、速度估計(jì)精度最差；EKF-AUN精度比EKF-Q高，但比NMMSE、QEKF精度差。這是因?yàn)镋KF-AUN算法基于AUN假設(shè)，并不完全符合量化噪聲特性，從而導(dǎo)致精度較差。相比之下，本文提出的QEKF方法的精度最高，相對(duì)于傳統(tǒng)基于AUN假設(shè)的方法位置估計(jì)精度提升了23%，表明了所提出算法在量化量測(cè)目標(biāo)跟蹤問題中的優(yōu)越性。此外，可以看出QEKF和NMMSE的誤差曲線比EKF-Q和EKF-AUN更加平滑，說(shuō)明QEKF和NMMSE方法可以減少更多的不確定性。

圖6 不同算法位置估計(jì)誤差Fig.6 Position estimation error of different algorithms

圖7 不同算法速度估計(jì)誤差Fig.7 Velocity estimation error of different algorithms

表1 不同方法10-60秒平均RMSETab.1 Average RMSE of different algorithms between 10～60 s

設(shè)置像素大小由10 μm到70 μm，進(jìn)行500次蒙特卡洛仿真，得到平均位置、速度估計(jì)精度和像素點(diǎn)大小的關(guān)系如圖8?？梢钥闯觯S著像素點(diǎn)尺寸減小，不同方法的估計(jì)精度都有所提升，且估計(jì)精度越來(lái)越接近。這是因?yàn)楫?dāng)像素點(diǎn)尺寸足夠小，則量化量測(cè)可以等效為高斯量測(cè)，量化噪聲的影響很小，不同方法之間區(qū)別較小。當(dāng)像素點(diǎn)尺寸越大，QEKF相對(duì)于其他方法的估計(jì)精度優(yōu)勢(shì)越明顯。

圖8 像素點(diǎn)大小和位置、速度估計(jì)精度關(guān)系圖Fig.8 The relationship between size of pixel and position,velocity estimation accuracy

5 結(jié) 論

本文建立了助推段目標(biāo)運(yùn)動(dòng)模型、天基紅外相機(jī)量化量測(cè)模型。從卡爾曼濾波框架出發(fā)，推導(dǎo)了量化量測(cè)條件下的量化擴(kuò)展卡爾曼濾波方法。通過數(shù)學(xué)仿真，驗(yàn)證了所提方法在量化噪聲條件下可提升軌跡跟蹤精度。此外，本文還分析了像素點(diǎn)大小對(duì)位置、速度估計(jì)結(jié)果的影響。結(jié)果表明，隨著像素點(diǎn)尺寸減小，估計(jì)精度逐漸提升。