一种改进的GM-C-CPHD 空间多目标跟踪算法

谢贝旭，张 艳，陈金涛，张任莉

（中山大学 航空航天学院，广东 深圳 518060）

0 引言

随着航天技术不断进步，太空探索活动逐渐增加，新发射的航天器及卫星增多，使得近地空间逐渐拥挤［1］；同时，空间的目标解体频发，产生的碎片数量逐渐增多，对卫星和航天器的运行造成巨大的威胁［2-3］。针对空间碎片，主要有以下4 个应对措施：监测和预警、碰撞规避和防护、离轨和弃置策略，以及主动清除［4-5］。从应对措施看，对轨道上目标进行精准识别［6］、定位与跟踪均为应对的必备前提，以此提前掌握可能发生碰撞的信息，并对应地操控卫星躲避或提前精确地清除碎片。

多目标跟踪（Multi-Target Tracking，MTT）［7-10］是指从含噪声的传感器量测中（存在虚警或杂波），估计监视区域内的目标数量和状态。MTT 算法主要分为数据关联和非数据关联2 类，如图卷积网络（Graph Convolutional Network，GCN）、联合概率数据关联（Joint Probabilistic Data Association，JPDA）和多假设跟踪（Multi-hypothesis Tracking，MHT）是基于数据关联的MTT 算法［11-12］，需先完成量测的分配步骤，再进行目标状态估计［13］。非数据关联的方法，是在随机有限集（Random Finite Set，RFS）的框架下进行，如概率假设密度（Probability Hypothesis Density，PHD）和势概率假设密度（Cardinalized Probability Hypothesis Density，CPHD）滤波算法，可以避开关联运算模块，并直接进行多目标状态估计，显著缓解了运算的复杂程度［14-19］。

R.MAHLER 从贝叶斯滤波和有限集理论出发，提出PHD 滤波算法［20］，该算法递归多目标后验Ⅰ阶近似矩，而不是整体地后验概率密度。随后又提出CPHD滤波算法，联合传递多目标概率密度的Ⅰ阶近似矩及目标的势分布，进一步提高跟踪精度。但实现PHD、CPHD 算法存在积分的困难；VO 等［21］在此基础上实现了高斯混合形式的PHD 及CPHD 滤波算法，降低了运算难度，并在二维仿真场景中应用。空间多目标跟踪问题的突出难点之一是由摄动力引起的运动方程不确定性，与其相关的不确定性参数会导致状态估计难度加大，进而导致环境对目标的影响不确定，参数识别效果变差，跟踪性能变差［22］。为了提高对空间目标的跟踪性能，MCCABE 和DEMARS 等［23-24］推导了一种考虑不确定性参数的高斯混合势概率假设密度（Gaussian-mixture Cardinalized Probability Hypothesis Density，GM-C-PHD）滤波器；YANG 等［25］在无迹卡尔曼滤波框架下推导了考虑面质比（Area-to-Mass Ratio，AMR）不确定性参数（Gaussian-mixture Considering Probability Hypothesis Density，GM-CPHD）的滤波器，在空间跟踪环境中使用。CPHD 滤波器同时传递多目标概率密度的Ⅰ阶近似矩以及目标的势分布，因此跟踪精度上，尤其是对目标数目的估计，其误差小于PHD。

针对空间目标跟踪的模型不确定性问题，本文提出一种考虑面质比（Area-to-Mass Ratio，AMR）不确定性参数（Gaussian-mixture Considering Cardinalized Probability Hypothesis Density，GM-CCPHD）的算法。该算法在轨道动力学模型中考虑了不确定参数AMR，在无迹卡尔曼滤波（Unscented Kalman Filter，UKF）滤波框架下通过协方差传递的方式传递参数AMR 对位置、速度状态估计的影响，以此减弱模型的不确定性影响程度，提高目标跟踪水平。仿真分析表明，相较于原来的GMCPHD 滤波器，目标的跟踪性能有所改善。

1 改进的GM-C-CPHD 滤波器

CPHD 滤波器在计算上存在积分困难，可通过高斯混合形式将递归过程中的目标强度函数及势分布函数进行高斯叠加和处理，进而有效避免积分困难的问题。本节首先简单分析不确定性参数AMR 对空间目标摄动力模型的影响，随后提出一种考虑不确定性参数AMR 的GM-C-CPHD 滤波器，并阐述该算法的具体计算流程。

1.1 AMR 与摄动力的影响关系

在地心惯性坐标系中，空间目标的摄动力模型如下：

式中：r为空间目标的位置矢量，m；为位置矢量加速度，m/s2；αarp为太阳光压加速度，m/s2；αdrag为大气阻力加速度，m/s2；αΘ为太阳引力加速度，m/s2；αΠ为月球引力加速度，m/s2；αJ2为J2摄动加速度，m/s2。

与AMR 参数相关的有大气阻力摄动和太阳光压摄动，本文主要考虑AMR 参数对太阳光压摄动力的影响。其影响关系表达式如下：

式中：S/m为面质比；rΘ为太阳的位置矢量，m；Cr为无量纲反射系数；Pr为地球表面处的光压强度；re，Θ为地日心平均距离。

1.2 状态增广的系统方程

本文采用高斯混合的实现形式，并在UKF 滤波框架下传递概率密度函数及势分布函数，以适应非线性的系统方程。目标的状态方程和量测方程如下：

式中：f(·)和g(·)为非线性函数；wk和vk均为高斯白噪声；c为面质比参数，即c=S/m；xk，c为考虑AMR 参数的增广状态矢量；zk为观测矢量。

xk，c和zk如下：

式中：A，E分别为方位角及俯仰角，rad；，分别为两者的角速度，rad/s；x、y、z分别为以观测站为原点的三维空间坐标，m；，，分别为对应的速度，m/s。

1.3 高斯混合实现形式

1）初始化

在非线性跟踪系统中，为了降低计算难度，利用高斯近似方法，将后验强度中的非高斯分量拆解成若干高斯分量形式。这里将满足均值为m、协方差为P的高斯概率密度记为N(·；m，P)。假设k时刻新生目标的强度函数是高斯混合形式：

2）预测步

记k-1 时刻真实目标的后验强度函数为

则k时刻的预测强度函数可计算如下：

式中：vS，k|k-1(xc)为存活目标的强度。

存活目标的强度计算如下：

式中：pS，k为新生目标存活概率。

记后验势分布为pk-1(·)，则k时刻的预测势分布函数如下：

式中：pΓ，k(·)为k时刻新生目标的势分布函数；pS，k为新生目标存活概率；为组合系数；Jk-1和Jk|k-1为高斯混合分量的个数；ω为分量权重。

3）更新步

k时刻目标的后验势分布及后验强度也是高斯混合形式，且两者的更新如下：

式中：vD(·)为检测概率下的更新强度函数。

1.4 滤波器实现

在使用高斯混合形式描述后验强度和后验势分布后，使用UKF 滤波公式传递各个高斯分量，更新目标强度函数和势分布，并实现GM-C-CPHD 滤波器。具体的滤波实现过程如下。

1）基于k-1 时刻的状态估计量开始Sigma 点采样，获取1 组Sigma 点及其相对应的权值wl：

2）将Sigma 点代进运动方程，获取Sigma 点的预测情况：

3）根据权值wl及各Sigma 点的预测值获取状态变量及协方差矩阵的预测情况如下：

4）根据上一步计算的状态变量预测值再次进行Sigma 点采样获取另一组Sigma 点及相对应的权值wl：

5）将上一步获取的Sigma 点代入观测方程，求得各Sigma 点的观测预测值为

6）根据权值wl及各Sigma 的观测预测值获取系统观测及协方差的预测值为

7）计算增益矩阵Kk为

8）计算状态变量及协方差矩阵的更新值：

2 仿真分析

为了比较GM-C-CPHD 滤波器与GM-CPHD滤波器在不同空间目标跟踪环境下的跟踪效果，共设置3 组仿真实验，其变量控制见表1，3 个目标初始轨道根数见表2。其中目标1 存在时间为第1至第30 时刻；目标2 存在时间为第1 至第21 时刻；目标3 存在时间为第5 至第30 时刻。目标的滤波初始状态标准差设置为：x、y、z方向的位置标准差为σx=σy=σz=1 km；x、y、z方向的速度标准差为σvx=σvy=σvz=1 m/s。地心地固坐标系（Earthcentered，Earth-fixed，CEF）测站参数设置见表3。使用最优子模式分配（Optimal Subpattern Assignment，OSPA）一致性度量进行性能评估，其中惩罚参数c=50 km，阶次参数p=2。

表1 仿真实验参数设置Tab.1 Parameter settings of the simulation tests

表2 地球同步轨道目标的初始轨道根数Tab.2 Initial Keplerian elements of the geosynchronous orbit objects

表3 地心地固坐标系下测站参数说明Tab.3 Description of the ground station parameters in the coordinate system

2.1 仿真实验1

GM-C-CPHD 和GM-CPHD 滤波的位置、速度OSPA 误差及目标数目估计结果如图1 所示。其中，第30 时刻的位置、速度OSPA 误差统计值及目标数目错误估计的时刻数见表4。可以看到，随着时间的迭代，2 个滤波器最后均能收敛，但GM-CCPHD 滤波器的位置及速度OSPA 误差更小。GM-C-CPHD 滤波器的位置和速度OSPA 误差的最后收敛数值分别为130.311 和0.019，而GMCPHD 滤波器则为3951.62 和0.74。在目标数目估计方面，图1 显示2 个滤波器在大部分时间均可以准确估计，其中GM-C-CPHD 在第15 时刻目标数目估计错误，估计错误时刻数为1；而GM-CPHD滤波器在第15 和18 时刻均估计错误，错误估计时刻数为2。2 个滤波器在第30 时刻不同分位值的OSPA 误差见表5。由表5 可知，GM-C-CPHD 滤波器在第30 时刻的估计精度更高，显示出所提方法的有效性。

图1 OSPA 误差和目标数目估计Fig.1 Estimation of the OSPA errors and target numbers

表4 第30 时刻位置、速度OSPA 误差和目标数目错误估计时刻数Tab.4 OSPA errors of the position and velocity at the 30th moment and the total moment number of false estimation regarding the target number

表5 第30 时刻不同分位值的OSPA 误差Tab.5 OSPA errors regarding different quantile values at the 30th moment

2.2 仿真实验2

仿真实验2 进一步降低了检测概率，设置为0.95，其他实验参数设置与仿真实验1 一致。GM-C-CPHD 和GM-CPHD 滤波的仿真结果如图2、表6 和表7 所示。由图2 可知，2 个滤波器的跟踪误差值最后同样收敛，GM-C-CPHD 滤波器的位置及速度OSPA 误差收敛于更小的数值。由表6 可知，本文方法收敛数值分别为189.751 和0.025 3，而GM-CPHD 滤波器则为7 548.05 和1.463 8。在目标的数目估计方面，由于检测概率减小，2 个滤波器的估计效果相比于仿真1 有所降低。其中GM-CCPHD 滤波器在第11、15 和第22 时刻估计不准确；而GM-CPHD 滤波器在第11、12、15 和第22 时刻估计不准确。此外，表7 表明在5%分位值情况下，2 个滤波器的OSPA 误差较为接近，但其他分位值情况下，GM-C-CPHD 滤波器的位置及速度估计误差小于GM-CPHD 滤波器。

图2 GM-CPHD 和GM-C-CPHD 的OSPA 误差和目标数目估计Fig.2 OSPA errors and estimated target number

表6 第30 时刻位置、速度OSPA 误差和目标数目错误估计时刻数Tab.6 OSPA errors of the position and velocity at the 30th moment and the total moment number of false estimation regarding the target number

表7 第30 时刻不同分位值的OSPA 误差Tab.7 OSPA errors regarding different quantile values at the 30th moment

2.3 仿真实验3

在仿真2 的基础上，进一步考虑杂波的情况，杂波的产生服从泊松分布，泊松分布参数λ=1，其他参数设置与仿真实验2 保持一致。GM-C-CPHD 和GMCPHD 滤波的仿真结果如图3、表8 和表9 所示。

图3 GM-CPHD 和GM-C-CPHD 的OSPA 误差和目标数目估计Fig.3 OSPA errors and the estimated target number by GM-CPHD and GM-C-CPHD

表8 第30 时刻位置、速度OSPA 误差和目标数目错误估计时刻数Tab.8 OSPA errors of the position and velocity at the 30th moment and the total moment number of false estimation regarding the target number

表9 第30 时刻不同分位值的OSPA 误差Tab.9 OSPA errors regarding different quantile values at the 30th moment

如图3 所示，GM-C-CPHD 滤波器的位置及速度OSPA 误差更小，跟踪精度更优。由表8 可知，本文方法位置及速度OSPA 误差收敛数值分别为154.473 0 和0.023 2，而GM-CPHD 滤波器则为13 843.3 和2.737 4。在目标数目估计方面，由于跟踪场景存在杂波，估计效果相比于仿真2 均有所减小。其中GM-CPHD 滤波器在第11、14、15、19、20和21 时刻有估计错误情况，而GM-C-CPHD 滤波器只在第14 和15 时刻估计存在偏差。此外，表9显示，2 个滤波器在第30 时刻不同分位值的OSPA误差，可知GM-C-CPHD 滤波器在各分位值情况下的位置及速度估计误差均小于GM-CPHD 滤波器。

2.4 其他指标比较及分析

除了使用OSPA 指标评价滤波器的性能外，还统计了2 个滤波器在所有时刻的OSPA 误差均值、标准差及均方根误差，统计结果见表10。由表10 可知，GM-C-CPHD 滤波器的OSPA 误差标准差比GM-CPHD 滤波器的标准差大，说明后期误差减小幅度更大，而平均值及均方根误差均比GM-CPHD滤波器小，表明了GM-C-CPHD 滤波器对目标状态估计的整体误差水平更低。

表10 OSPA 误差统计值Tab.10 Statistical OSPA errors for all the tests

综合仿真结果及其他指标进行分析，由于不确定性参数AMR 对轨道动力学模型存在影响，使得滤波过程的协方差矩阵不能清晰地反映多目标跟踪的精度。其中，GM-CPHD 滤波器在递归过程中将AMR 参数假设为常数，参数的影响效果就会确定化，那么滤波器估计参数的不确定性也会受到影响，进而在协方差中表现出来。而GM-C-CPHD 滤波器通过协方差增广的方式将AMR 参数考虑到协方差传递过程中，提高该不确定性参数的可观性，进而降低其对滤波估计的影响，使得目标跟踪精度得到改善。

3 结束语

本文针对空间目标运动模型的不确定性问题，提出一种考虑面质比不确定性参数的GM-CCPHD 多目标跟踪算法，该算法在轨道动力学模型中考虑了面质比不确定参数，并在UKF 滤波框架下，通过协方差传递的方式传递该参数对位置、速度状态估计的影响，以此缓解空间目标运动模型的不确定性影响。通过考虑降低检测概率、增加杂波率等条件设置3 组仿真实验，仿真结果表明，与GMCPHD 滤波器相比，GM-C-CPHD 滤波器在OSPA误差及其均值、均方根误差等指标方面均有所改善。未来，主要从以下2 个方面开展工作：1）跟踪更多数量且密集的空间多目标；2）进一步优化算法，考虑检测概率的不确定性影响。