基于幅度信息的海上多目标稳健跟踪方法

柳超，张志国，孙进平

（1. 北京航空航天大学电子信息工程学院，北京 100191；2. 海军92853部队，辽宁 葫芦岛 125106）

1 引言

多目标跟踪是当前信息融合和计算机视觉领域的热点问题。检测的不确定性、量测的不确定性及关联的不确定性使多目标跟踪非常复杂[1]，尤其在雷达对海探测环境下，目标速度慢、杂波剩余多等因素使目标连续跟踪异常困难。目前，常见的多目标跟踪方法主要包括联合概率数据关联（JPDA,joint probabilistic data association）、多假设跟踪（MHT, multiple hypotheses tracking）、随机有限集（RFS, random finite set）等[2]。JPDA和MHT的实现过程是先进行数据关联，再进行单目标滤波。当量测数目增加时，数据关联运算量呈指数级增长。近年来，基于 RFS的多目标跟踪方法受到了很多关注，该方法将多目标状态和传感器量测分别建模为随机有限集，通过贝叶斯多目标滤波公式实现多目标后验概率密度的递归估计。与JPDA和MHT等传统跟踪方法不同的是，许多RFS方法不考虑数据关联也可以得到准确的多目标滤波结果，如概率假设密度（PHD, probability hypothesis density）滤波器[3]、势概率假设密度（CPHD, cardinalized PHD）滤波器[4]、多伯努利（MB, multi-Bernoulli）滤波器[5]。这3种方法由于实现过程简单、计算效率高等优点获得了大量应用，尤其是MB滤波器，其计算复杂度和 PHD相似，其序贯蒙特卡罗（SMC,sequential Monte Carlo）实现不需要聚类便能提取目标状态，因而在多目标滤波中更具优势。

在基于RFS的多目标滤波器中，杂波强度（即杂波率和杂波概率密度之积）和检测概率需要根据应用环境设定，但在雷达对海探测中，实时获取这2个参数非常困难。如果设定的参数值与真实值不一致，滤波器的性能便会下降。近年来，杂波强度和检测概率未知情况下的目标跟踪问题得到了深入研究。针对检测概率未知的情况，Guiseppe等[6]提出2种自适应跟踪方法；针对杂波强度未知的情况，学者们提出了多种杂波估计算法[7-10]；针对杂波强度和检测概率均未知的情况，Mahler等[11]提出一种在滤波的同时学习杂波率和检测概率的自适应CPHD滤波器，随后，Beard等[12]在此基础上提出一种自举滤波器[12]。文献[11-12]所提方法采用SMC实现时，需要进行聚类以提取目标状态，当目标的势估计误差较大时，状态估计性能也会变差。针对这一问题，Vo等[13]提出了稳健多伯努利滤波器（RMB, robust MB），因其SMC实现不需要聚类便能提取目标状态，该方法已应用于传感器选择[14]、生物学研究[15]、视频跟踪[16]等领域。文献[17]将稳健滤波思想[13]引入标签多伯努利滤波器，但没有给出杂波率的准确估计。

RMB滤波器在计算量测似然时，只利用目标和杂波的运动状态（如位置、速度等）信息，当目标和杂波相距较近时，多目标滤波效果欠佳。由于雷达回波不仅包含目标的运动状态信息，还包含幅度信息（AI, amplitude information），利用目标的幅度信息来改进量测似然，提升对目标和杂波的区分能力，实现更加准确的多目标滤波是一个可行的方法。

利用幅度信息的目标跟踪方法主要包括检测前跟踪（TBD, track-before-detect）[18-19]和幅度信息辅助目标跟踪（AIAOT, AI aided object tracking）[20-29]。前者对连续多帧回波进行非相参积累以提取可能的目标航迹，对低信杂比目标（SCR, signal to clutter ratio）跟踪效果好，但由于每一帧都要对扫描区域进行遍历搜索，因此运算量很大；后者是在传统跟踪方法的基础上增加幅度似然的计算，能够显著提升算法的跟踪效果，并且运算量小。

目前，AIAOT领域已取得许多成果。Lerro等[20-21]首先提出了 AIAOT的基本理论，随后Ehrman等[22]和 Brekke等[23-24]对 AIAOT 进行了深入研究。文献[20-24]只针对单目标跟踪场景。对于多目标跟踪的情况，Clark等[25-26]将Rayleigh杂波的幅度信息引入PHD滤波器，随后一些学者进行了更加深入的研究[27-28]。针对未知杂波强度的情况，Yuan[29]将 Rayleigh杂波和幅度恒定目标的AI引入RMB滤波器。

在雷达对海探测中，运动目标的姿态和视角相对雷达经常变化，导致目标的雷达散射截面积（RCS,radar cross section）起伏不定，从而造成目标回波幅度起伏。目前，表征目标 RCS起伏的模型包括Swerling模型和χ2分布模型等，工程上常用SwerlingI～SwerlingIV 模型，其中，Swerling I型表示RCS在扫描内的脉冲间相关，而两次扫描相互独立，为慢起伏，其概率密度服从瑞利分布。大量海杂波实测数据分析表明，通常海杂波相关时间为毫秒级[30]，服从慢起伏模型，且采用K分布模型能有效拟合较大入射余角范围内的海杂波幅度分布。目前，Swerling I型已广泛应用于描述海上目标RCS的起伏[23]。

本文提出了一种K分布海杂波环境下AI辅助的RMB滤波器（AI-RMB）。首先，将K分布海杂波及Swerling I型起伏目标的幅度似然函数引入稳健多伯努利滤波器，提升对海探测环境下雷达的多目标跟踪性能；其次，与经典 AIAOT方法[20-29]利用幅度信息计算理论检测概率和理论虚警概率不同，本文考虑到RMB滤波器的特点，采用贝塔分布迭代估计目标伯努利项的检测概率和杂波伯努利项的虚警概率；最后，为提升计算效率，采用积分表方法计算K分布杂波下的幅度似然函数。仿真实验表明，相较稳健多伯努利滤波器，本文所提方法在多目标状态估计、势估计以及杂波率估计方面性能更优，且运行时间没有显著增加。

2 基本模型

本节给出单个目标运动状态的时间演化模型和雷达的观测模型。

2.1 目标运动模型

假设k时刻第i个目标为近似匀速运动，其运动模型如式（1）所示。

其中，T和σv分别为采样间隔和加速度噪声的功率谱密度。

2.2 雷达观测模型

k时刻第i个量测矢量由运动状态相关的量测部分和附加的幅度部分构成。

其中，h（⋅）为非线性量测函数，为零均值白色高斯量测噪声，量测函数以及量测噪声协方差分别为

3 稳健多伯努利滤波器

为了在滤波的同时估计出环境中的杂波率，Vo等[13]提出的RMB滤波器将杂波建模为一类特殊类型的目标（又称为杂波发生器），并分别建立了目标和杂波的状态转移模型和量测模型，通过贝叶斯迭代滤波实现对目标状态和杂波率的估计。

定义 χ（Δ）= [ 0 ,1]为检测概率的取值空间，χ~=Rnx为目标的状态空间，nx为状态空间的维度，{0 ,1}为杂波发生器和目标的离散标签空间（本文用u =0表示杂波发生器，u=1表示真实目标），则目标和杂波的增广状态空间为目标和杂波的增广状态为其中，α为检测概率，为运动状态，u为标签。定义在增广状态空间上的任意函数，有

假设在k-1时刻，多目标概率密度可表示为如式（3）所示的多伯努利随机集。

其中，有

其中，fΔ,u和fχ~,u分别为检测概率和运动状态的转移密度。

将式（4）中预测的多目标密度表示为如式（8）的多伯努利随机集。

则k时刻后验多目标密度可表示为漏检多伯努利随机集和经量测更新的多伯努利随机集的并集。

且有

4 基于AI的稳健多伯努利滤波器

在上述RMB迭代中，多目标后验密度更新时仅利用了与运动状态相关的量测。本文引入目标和杂波的幅度信息，建立了目标和杂波的幅度量测似然，并将其加入更新过程，从而构成AI-RMB滤波器。由于AI-RMB和RMB的预测部分完全相同，为节省篇幅，这里只给出AI-RMB的更新部分。

4.1 幅度的似然函数

假设目标的幅度d和运动状态相互独立，则目标的量测似然函数 g（z|x）和杂波的量测似然函数c（z）可分别表示为[21]

其中，τ为检测门限。则经过门限检测后，目标和杂波的幅度似然函数分别为

4.2 AI-RMB滤波器

如果在k时刻，预测的多目标密度可表示为如式（22）所示的多伯努利随机集。

其中，有

在式（25）、式（26）、式（28）和式（29）中，目标和杂波的检测概率pD,u,k虽然可以通过式（18）和式（19）求得理论值，但该值不能直接应用于AI-RMB的更新过程。其原因如下：一方面，与仅用于估计多目标状态的MB滤波器不同，RMB滤波器不仅需要估计目标状态，还要在滤波的过程中估计出杂波率，由于通常情况下虚警概率比检测概率低得多（相差几个数量级），如果采用式（18）和式（19）计算出的理论检测概率，将会导致杂波伯努利项的权重过低，经过航迹修剪和合并，这些项很可能丢失，从而无法实现对杂波率的估计；另一方面，由式（18）计算出的理论检测概率并不能准确反映RMB中每个目标伯努利项的检测概率，这是因为只有在RMB滤波器中的目标伯努利项数与实际目标数完全一致的情况下才能准确表示目标伯努利项的检测概率，同样，由式（19）计算出的理论虚警概率也不能准确反映RMB中每个杂波伯努利项的检测概率。基于上述两方面原因，本文在计算检测概率和虚警概率时，摒弃了常规AIAOT[20-29]中所采用的方法，而是采用贝塔分布来迭代估计每个伯努利项的检测概率和虚警概率。

4.3 K分布海杂波及目标的幅度似然函数

幅度服从K分布的海杂波具有如式（31）所示的概率密度函数（PDF, probability density function）。

其中，G（⋅）表示Gamma函数，Kv表示第二类修正贝塞尔函数，v是形状参数，b是尺度参数。则根据式（19）可得海杂波的虚警概率为

其中，τ为检测门限。

通过K分布杂波的复合解释[23]，可得Swerling I型起伏目标加杂波的PDF为

则根据式（18）可得目标的检测概率为

将式（31）和式（32）代入式（21）可得

将式（33）和式（34）代入式（20）可得

其中，有

式（36）中分子和分母的积分都没有闭合解。针对这一问题，Brekke等[23]提出了一种采用网格的数值积分方法，以对g（η）的积分为例说明如下。

首先进行变量代换，如式（37）所示。

则式（36）中分子部分的积分可以表示为

其中，有

q（u）的极值up可通过对 q（u）求一阶导数得到，如式（39）所示。

当u→∞时，q（u）与u-1-2v成正比。为了估计q（u）的积分，使用一个变分辨率的积分网格，该网格包含低区和高区两部分，分别对应 q（u）极大值附近的积分和u→∞时的积分。对低区网格，采样点为

对高区网格，采样点为

这一数值积分方法准确性很高，但是计算量很大，不利于实时应用。为实现式（36）中积分的快速计算，本文制作了积分表（InT, integral table）。当需要计算式（36）的积分时，只需要根据K分布杂波的参数和量测幅度到积分表中查找相应的积分值即可。由于积分表是在算法执行前离线制作的，不消耗运行时间，因而可使海杂波中目标幅度似然的计算效率大大提升。具体操作步骤如下。

步骤1根据杂波参数b、v以及可能的信杂比，设定量测幅度a的可能取值。为了保证积分表在较大的信杂比范围内可用，a的取值范围要尽可能大。

步骤2利用Brekke提出的网格法计算杂波参数b、v以及量测a取不同值时的积分值，制作成积分表。

步骤3根据参数b、v以及a的取值，从积分表中提取相应的积分值，完成式（36）中积分的计算。

5 仿真

为检验所提AI-RMB的性能，通过仿真实验将其与MB滤波器和RMB滤波器进行对比。对海雷达观测区域如下：方位为0°～180°，距离为0～2 000 m。在总计100 s的观测时间内先后有10个目标进入观测区域，每个目标都为匀速运动。目标的真实轨迹如图1所示，其中，○表示航迹起点，△表示航迹终点。

图1 目标真实轨迹

每一帧中杂波数服从均值为1 000的泊松分布，海杂波幅度服从参数为v=4，b=0.25的K分布，空间分布为方位-距离上的均匀分布，但在直角坐标系中为非均匀分布。所有目标的平均信杂比为

目标的幅度测量值服从Swerling I型起伏。虚警概率设为 0.02，则根据式（32）可求得该信杂比下的检测门限。当杂波和目标的量测幅度超过检测门限时生成量测点迹。理论上检测后杂波的数目服从均值为20的泊松分布，且在方位-距离上服从均匀分布，但在平面直角坐标系上为非均匀分布。检测后包含目标和杂波的量测点迹如图2所示。

图2 量测点迹

目标的存在概率为 pS,1,k=0.99，虚警的存在概率为 pS,0,k=0.9。目标和杂波的检测概率的转移密度为如下的贝塔分布。

其中，有

并给定均值 μa,1,k|k-1=ak-1，对于目标，标准差σa,1,k|k-1=0.01，而对于杂波，标准差为σa,0,k|k-1=0.07。

目标的新生项为多伯努利随机集，参数为

新生目标贝塔分布的参数为 sΓ,1,k=85，tΓ,1,k|=15，运动状态的概率密度为

杂波的新生项也是多伯努利随机集，其参数为

本文通过计算每种算法的 OSPA位置误差和OSPA势误差[31]来评价算法的跟踪性能。

其中，X={ x1,… ,xm}和Y={ y1,… ,yn}分别表示目标的真实状态和估计状态，m、n∈ℕ0={ 0 ,1,2…}。本文中设置p=1和c=300。

图3～图5为MB、RMB和AI-RMB在某次实验中的滤波效果。从图3～图5可以看到，参数设置完全正确的MB性能最优，只产生了很少的漏检和虚假航迹；AI-RMB次之，产生了更多的漏检和虚假航迹；未采用幅度信息的RMB性能最差，产生了大量漏检和虚假航迹。

图3 单次实验中MB滤波结果

图4 单次实验中RMB滤波结果

图5 单次实验中AI-RMB滤波结果

图6给出了MB、RMB和AI-RMB进行30次仿真后的平均 OSPA误差，其中，图 6（a）是 OSPA位置误差，图6（b）是OSPA势误差。从图6可以看到，参数配置正确的MB位置误差最小；AI-RMB的位置估计和势估计误差都明显小于未采用幅度信息的RMB。这是因为，RMB在计算量测似然时只利用了目标和杂波的运动信息，当目标和杂波距离较近时难以区分，从而影响目标粒子和杂波粒子权重的计算，造成目标状态估计性能下降。而采用幅度信息可以更好地区分目标和杂波，从而提升目标的状态估计性能。

图7是MB、RMB和AI-RMB在30次仿真后的平均势估计效果。从图7可以看到，这3种算法都能比较准确地估计出目标数目。但是，当目标数增多时（即第60 s～80 s），未采用幅度信息的RMB产生了比较明显的欠估计，而采用幅度信息辅助的AI-RMB以及参数配置完全正确的MB对目标数目的估计更加准确。

图6 30次实验中平均OSPA位置与势误差

图7 30次实验中平均势估计结果

图8是RMB和AI-RMB对杂波数目的估计结果，其中，图8（a）为某次仿真中每一帧的真实杂波数目和2种算法估计的杂波数目，图8（b）为30次仿真后2种算法对杂波数目估计的均方根误差对比。从图8可以看出，在开始的大约20帧中，2种算法对杂波数目的估计都存在比较大的误差，且误差大小一致，但误差随滤波时间的增加而迅速降低。在20帧之后，2种算法的误差基本上都达到了一个稳定的状态，容易看到采用幅度信息辅助的AI-RMB的误差小于未采用幅度信息的RMB。

表1为通过30次仿真实验统计出的MB、RMB、未采用积分表的AI-RMB和采用积分表的AI-RMB（InT）的平均运行时间。MB、RMB和AI-RMB的计算复杂度均与当前时刻的量测数和目标数成正比，即O（|Zk|×|Xk|）。在本文仿真中，这 3种滤波器处理的量测数是相同的，区别仅在于目标数。从表1可以看到，MB的运行时间最短，这是因为该滤波器仅用于估计真实目标，因而每个时刻的|Xk|都比较小。RMB的运行时间约为MB的6倍，这是因为该滤波器需要同时估计目标和杂波，因此每个时刻的|Xk|都比较大。从仿真设置可以看到，在每个时刻，MB的新生目标项数为4，而RMB为24，为MB的6倍，因此其计算复杂度也应为MB的6倍，这与仿真结果是一致的。理论上，AI-RMB的计算复杂度与RMB相同，但由于幅度似然函数的计算中需要对K分布求积分，而该积分没有解析表达式，只能通过数值积分求解，因此大大消耗了运行时间。从表1还可以看到，未采用积分表的AI-RMB的运行时间远远大于 RMB，而采用积分表的AI-RMB（InT）的运行时间则降低到接近RMB的水平。

图8 30次实验中平均杂波数估计和均方误差

表1 30次实验中平均运行时间对比

6 结束语

稳健多伯努利滤波器适用于杂波强度和检测概率未知的场景，因而在雷达对海探测环境下具有重要的应用价值。但该滤波器仅利用与目标运动状态相关的量测计算似然函数，当目标和杂波距离较近时性能不佳。针对这一问题，本文基于K分布海杂波和 Swerling I型起伏目标建立了幅度似然函数，并将其引入该滤波器的更新过程，以提升算法对目标和杂波的区分能力；在配置检测概率时，舍弃常规方法中利用目标和杂波幅度计算理论检测概率和理论虚警概率的思路，采用贝塔分布迭代估计每个目标伯努利项的检测概率和每个杂波伯努利项的虚警概率。仿真实验表明，本文所提方法在目标状态估计、势估计及杂波数估计方面均优于稳健多伯努利滤波器。且通过使用积分表，所提方法的计算效率与稳健多伯努利滤波器相当。