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

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

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

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

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

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

1 助推段目標運動模型建立

在助推階段的目標加速度可以表示為：

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

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

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

其中，ωe為地球旋轉速度。且：

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

對于t≥t0，有：

其中，β（t）為輔助參數。

求導可得：

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

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

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

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

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

2.1 理想量測方程

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

圖1 像平面幾何關系Fig.1 Geometric relationship of image plane

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

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

2.2 量化量測模型

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

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

從離散圖像平面的輸出，只能確定實際坐標所在的像素，而無法得到準確坐標。為了對天基紅外預警衛(wèi)星的量測進行精確建模，引入mid-riser量化器，將量測值坐標輸出可以看作是mid-riser量化器的輸出。mid-riser量化器的數學表達式如下：

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

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

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

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

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

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

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

定義

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

新息協(xié)方差Sk：

其中，Hk為雅克比矩陣。

濾波增益Kk：

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

其中，

Kk可根據式(29)計算。

其中，

整理得：

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

協(xié)方差：

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

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

涉及到的高斯積分為：

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

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

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

4 數學仿真

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

4.1 仿真場景

目標助推段彈道如圖3。

圖3 目標助推段彈道Fig.3 Boost-phase Trajectory

仿真參數：兩顆衛(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 仿真結果

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

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

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

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

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

開展500次蒙特卡洛仿真，得到如圖6和圖7的不同算法位置、速度估計誤差對比結果。將不同算法10 s～60 s位置、速度的RMSE取平均，得到結果如表1。可以看出，EKF-Q的位置、速度估計精度最差；EKF-AUN精度比EKF-Q高，但比NMMSE、QEKF精度差。這是因為EKF-AUN算法基于AUN假設，并不完全符合量化噪聲特性，從而導致精度較差。相比之下，本文提出的QEKF方法的精度最高，相對于傳統(tǒng)基于AUN假設的方法位置估計精度提升了23%，表明了所提出算法在量化量測目標跟蹤問題中的優(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

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

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

5 結 論

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