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

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

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

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

天基紅外預(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)跟蹤問題。但相對于傳統(tǒng)EKF方法計算量較大。

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

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

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

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

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

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

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

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

對于t≥t0，有：

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

求導(dǎo)可得：

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

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

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

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

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

2.1 理想量測方程

假設(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時刻第i顆衛(wèi)星的理想量測為。由和組成，計算方式為：

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

2.2 量化量測模型

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

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

從離散圖像平面的輸出，只能確定實(shí)際坐標(biāo)所在的像素，而無法得到準(zhǔn)確坐標(biāo)。為了對天基紅外預(yù)警衛(wèi)星的量測進(jìn)行精確建模，引入mid-riser量化器，將量測值坐標(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)星量測噪聲由兩部分組成：1) 高斯噪聲wy和wz；2) 量化噪聲ey和ez。在高斯噪聲大于量化噪聲的情況下，可延用傳統(tǒng)EKF濾波方法；當(dāng)量化噪聲更大時，量測的非高斯特性會引起EKF方法精度差，甚至是濾波發(fā)散。

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

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

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

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

定義

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

新息協(xié)方差Sk：

其中，Hk為雅克比矩陣。

濾波增益Kk：

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

其中，

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

其中，

整理得：

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

協(xié)方差：

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

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

涉及到的高斯積分為：

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

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

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

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

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

4.1 仿真場景

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

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

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

4.2 仿真結(jié)果

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

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

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

圖4 QEKF位置估計誤差Fig.4 Position estimation error of QEKF

圖5 QEKF速度估計誤差Fig.5 Velocity estimation error of QEKF

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

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

圖7 不同算法速度估計誤差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次蒙特卡洛仿真，得到平均位置、速度估計精度和像素點(diǎn)大小的關(guān)系如圖8?？梢钥闯?，隨著像素點(diǎn)尺寸減小，不同方法的估計精度都有所提升，且估計精度越來越接近。這是因為當(dāng)像素點(diǎn)尺寸足夠小，則量化量測可以等效為高斯量測，量化噪聲的影響很小，不同方法之間區(qū)別較小。當(dāng)像素點(diǎn)尺寸越大，QEKF相對于其他方法的估計精度優(yōu)勢越明顯。

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

5 結(jié) 論

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