天基红外预警卫星高精度协同目标跟踪方法

王 宇，马 成，彭 博，陶久亮，于一帆，赵 博

（北京宇航系统工程研究所，北京市 100076）

天基红外预警卫星常用于助推段的目标跟踪[1-5]，可被动探测目标热辐射，以多星协同的方式定位目标[6]。随着对定位精度的要求越来越高，传统假设天基红外预警卫星量测噪声的统计特性符合高斯分布，并基于高斯滤波方法[7,8]设计跟踪滤波器的方式已经无法满足要求，从而吸引了越来越多的学者开展高精度的协同跟踪方法研究。本文将针对天基红外预警卫星测量噪声的非高斯特性进行研究，建立相比于高斯分布更精确的量测噪声模型，并提出相应的状态估计方法。

天基红外预警相机的像平面是由离散像素组成的，投影在其上的连续光信号会转化为离散的电信号。在数字信号处理领域，用“量化”来描述连续信号被离散的过程，其一般由微处理器的有限计算精度或通信的带宽限制造成[9,10]。Widrow对量化的统计理论进行了深入的研究[11]，给出了量化器的定义，并对量化噪声进行了统计分析。Asmar等人[12]比较了几种用于垂直状态估计的滤波器，并说明了在量化噪声条件下传统Kalman滤波器性能较差。Rickard-Karlsson等人将量化噪声视为基于加性均匀噪声（AUN）假设的附加均匀噪声[13]，采用增大测量协方差的方法处理量化问题。Curry等人提出了一种通过计算以量化测量为条件的矩进行量化测量的方法[10]，这种方法包含复杂求积运算。为了提升算法效率，Duan[14]等人提出了一种量化测量下的近似MMSE（Minimum Mean Square Error）估计数值算法。可见，数字信号处理领域对连续信号被离散的过程进行了较深入的研究，本文将借鉴信号处理领域“量化”的概念，建立天基红外预警卫星的量测模型。

天基红外预警卫星目标跟踪系统是强非线性系统[15]。在各种非线性系统目标跟踪问题中，基于一阶线性化的扩展卡尔曼滤波（Extended Kalman Filter,EKF）算法应用最广。一些研究将容积卡尔曼滤波器[16,17]、粒子滤波[18-20]，应用于目标跟踪问题。但相对于传统EKF方法计算量较大。

本文主要工作包括：1）基于助推段运动微分方程，建立形式简单且符合飞行器运动特性的目标运动方程；2）从理论上分析天基红外预警卫星测量的特点,利用数字信号处理中的mid-riser量化器描述离散像平面，建立天基红外预警卫星的测量噪声模型；3）为解决基于天基红外预警卫星的高精度目标跟踪问题，提出了量化扩展卡尔曼滤波（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）为辅助参数。

求导可得：

定义状态向量x如下：

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

可见，该模型具备形式简单、不需要目标质量等先验信息的优点。

2 天基红外预警卫星量测模型建立

在实际应用中，至少有两颗卫星同时用于目标跟踪。为了表征天基红外预警卫星量测噪声的统计特性，高斯分布被广泛使用[6,7]。实际上，天基红外预警卫星相机的像平面是由离散的像素组成的，导致测量的量化。本节将分析天基红外预警卫星测量特性，建立量测模型。

2.1 理想量测方程

假设相机坐标系与卫星体坐标系重合，可表示为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颗卫星的理想量测为。由和组成，计算方式为：

其中，上标i表示第i个卫星；表示转换矩阵的第a行，第b列的元素。

2.2 量化量测模型

在天基红外预警卫星的离散像平面上，照射到像素内的光线均产生对应于该像素中心的坐标输出，见图2。目标在像素坐标系的坐标（,）受到量化噪声的干扰，干扰量最大为0.5像素。

图2 离散像平面Fig.2 Discrete image plane

从离散图像平面的输出，只能确定实际坐标所在的像素，而无法得到准确坐标。为了对天基红外预警卫星的量测进行精确建模，引入mid-riser量化器，将量测值坐标输出可以看作是mid-riser量化器的输出。mid-riser量化器的数学表达式如下：

基于mid-riser量化器，目标在像素坐标系中的理想坐标与量化坐标之间的关系可表示为：

其中，wy和wz为弱照射条件、高温等引起的零均值高斯白噪声；ey和ez为与像素大小相关的量化噪声。

从式(17)可看出，天基红外预警卫星量测噪声由两部分组成：1) 高斯噪声wy和wz；2) 量化噪声ey和ez。在高斯噪声大于量化噪声的情况下，可延用传统EKF滤波方法；当量化噪声更大时，量测的非高斯特性会引起EKF方法精度差，甚至是滤波发散。

将式(15)代入式(17)，得到卫星i量测方程：

3 量化量测协同目标跟踪方法

考虑如下量化量测非线性离散系统

其中，xk∊Rn为状态向量；yk∊Rm为量化前的量测；vk和wk为相互独立的零均值高斯噪声，协方差分别为Qk和Rk；下标k代表离散时间tk；f和h分别为状态方程和离散方程；ykq为量化后的量测向量。需要指出的是，即使量测方程h为线性方程、wk为高斯噪声，由于量化噪声的存在，Q（yk）将是非线性的，ek是非高斯的。

定义

QEKF和EKF的主要区别在于量测更新，EKF量测更新公式为：

新息协方差Sk：

其中，Hk为雅克比矩阵。

滤波增益Kk：

状态及协方差更新方程为：

其中，

Kk可根据式(29)计算。

其中，

整理得：

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

协方差：

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

QEKF中的高维高斯函数积分没有解析解。本文结合Genz的转换方法以及伪蒙特卡洛积分方法提出一种数值求解方法，将积分转换为单位超立方体上的积分，从积分区间直接取点，可避免维度灾难。

涉及到的高斯积分为：

式(44)可以进一步写为：

之后，利用伪蒙特卡洛方法对式(43)进行积分，即：

至此，完成基于Genz转换和伪蒙特卡洛的多维高斯积分计算。对于式(31)和式(37)中的积分，仅需对以上积分函数进行简单改动即可。

4 数学仿真

考虑两颗天基红外预警卫星协同跟踪助推段目标的场景，利用本文提出的协同目标跟踪方法对目标进行跟踪。

4.1 仿真场景

目标助推段弹道如图3。

图3 目标助推段弹道Fig.3 Boost-phase Trajectory

仿真参数：两颗卫星分别位于105°E、145°E地球同步轨道；两卫星量测噪声为：高斯噪声σ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和传统EKF方法（EKF-Q），基于AUN假设的改进EKF方法（EKF-AUN）[13]，近似最小均方差滤波方法（NMMSE）[14]，以及纯高斯量测下的EKF方法进行对比。定义均方根误差（RMSE）衡量算法精度为：

其中，M表示蒙特卡洛打靶次数，表示k时刻第i次状态估计结果；xk为k时刻真实状态。

图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方法的精度最高，相对于传统基于AUN假设的方法位置估计精度提升了23%，表明了所提出算法在量化量测目标跟踪问题中的优越性。此外，可以看出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。可以看出，随着像素点尺寸减小，不同方法的估计精度都有所提升，且估计精度越来越接近。这是因为当像素点尺寸足够小，则量化量测可以等效为高斯量测，量化噪声的影响很小，不同方法之间区别较小。当像素点尺寸越大，QEKF相对于其他方法的估计精度优势越明显。

图8 像素点大小和位置、速度估计精度关系图Fig.8 The relationship between size of pixel and position,velocity estimation accuracy

5 结 论

本文建立了助推段目标运动模型、天基红外相机量化量测模型。从卡尔曼滤波框架出发，推导了量化量测条件下的量化扩展卡尔曼滤波方法。通过数学仿真，验证了所提方法在量化噪声条件下可提升轨迹跟踪精度。此外，本文还分析了像素点大小对位置、速度估计结果的影响。结果表明，随着像素点尺寸减小，估计精度逐渐提升。