视场非完全重叠的分布式雷达多目标跟踪方法

达 凯 杨 烨 朱永锋 付 强

(国防科技大学电子科学学院ATR重点实验室 长沙 410073)

1 引言

分布式雷达多目标跟踪具有可扩展性强、跟踪精度高、抗干扰能力强等优点，且符合传感设备多站化、小型化、低功耗的发展趋势，目前已广泛应用于无人机群、监控预警、多弹协同等领域[1—5]。然而，受探测能力、波形设计以及天线指向等因素限制，不同雷达波束的覆盖范围可能不是完全重叠的(如图1)，这导致了网络中多目标检测信息的差异化，同时也给后续的分布式融合带来了巨大挑战。本文在利用高斯混合(Gaussian Mixture,GM)实现的集势概率假设密度(Cardinalized Probability Hypothesis Density,CPHD)滤波器[6]进行多目标跟踪的基础上，开展部分重叠视场下的分布式雷达多目标跟踪方法研究。

图1 两雷达的重叠视场Fig.1 The overlapping field of view of two radars

随机有限集多目标跟踪方法[7—12]在近些年获得了广泛的关注，其相比传统的多假设跟踪和概率数据互联方法，无须复杂的数据关联过程。CPHD滤波器作为随机集方法衍生的一种重要滤波器，相比于其他随机集滤波器有着独特的优势。例如，CPHD滤波器比PHD滤波器[13]的目标数目估计性能更好，比(标签)多伯努利滤波器[7,14]在分布式融合中的闭合性更好。目前，针对CPHD滤波器的分布式融合方法主要分为几何平均(Geometric Average,GA)和算术平均(Arithmetic Average,AA)融合，文献[15]对两种融合方法进行了详细的阐述。

由于分布式雷达网络中雷达视场仅是部分重叠，传统的分布式CPHD滤波器融合方法已经不再适用。目前，处理该问题的基本思想是：融合各雷达共同视场内的信息以提高跟踪性能，合并非共同视场内的信息以扩展感知范围。该思想的关键在于从表征全视场的多目标密度中找出代表部分共同视场的多目标密度。目前的方法大致分为3类：视场填充、视场切分以及聚类方法。

(1) 视场填充方法将所有雷达视场扩充至最大一致范围再进行融合。文献[16]和文献[17]在同时定位与构图的应用背景下，利用一个均匀密度将不同传感器的不同视场扩充为一个更大的相同视场，再利用标准的GA融合方法进行融合。这种方法的缺点是融合产生了大量的无效分量，这使得计算更为复杂。此外，由于非公共区域使用了概率质量较低的均匀分布进行填充，这导致融合结果会低估非公共区域目标的存在概率。

(2) 视场切分方法利用已知的雷达视场对多目标密度进行切分，进而融合切分结果以表示共同视场的多目标密度。文献[18]按照传感器视场的交叠对PHD进行切分，其将接收的粒子分类成共同视场内的公共粒子集和共同视场外的外部粒子集，然后通过GA方法融合公共粒子集的粒子，并直接合并外部粒子集。文献[19]同样利用视场交叠对粒子集进行切分，在有限视场下实现了PHD-AA 融合，但是它只融合了传感器共同视场中的信息，而将非共同视场的信息丢弃。文献[20]则通过多目标密度切分、融合以及合并的3步策略，开发了部分重叠视场下CPHD滤波器的分布式融合方法。

(3) 聚类方法根据不同度量寻找不同传感器对相同目标的估计匹配对，且仅在所选择的配对上进行融合。文献[21]计算不同传感器PHD高斯分量间的马氏距离，并通过GA方法融合马氏距离较近的高斯分量，实现了部分重叠视场下PHD滤波器的分布式GA融合。文献[22]采用和文献[21]类似的思想实现了部分重叠视场下CPHD滤波器的融合。文献[23]通过马氏距离来匹配多伯努利滤波器间的伯努利分量，实现了有限视场下多伯努利滤波器的分布式AA融合。文献[24]和文献[25]使用最优子模式分配(Optimal Subpattern Assignment,OSPA)距离作为度量来实现类似文献[21]中方法，不过文献[25]进一步考虑了标签匹配问题。文献[26]采用基于最高后验密度距离测度的方法实现了有限视场的 PHD滤波器融合。

总结以上，视场填充方法的缺点较为明显，会低估目标的存在概率。视场切分方法则要求传感器的视场精确已知，并按照视场的交叠进行多目标密度的切分。然而，一方面，在实际中尤其是在动态传感器网络中，视场的交叠非常复杂，基于视场来切分多目标密度十分困难；另一方面，要获得精确的传感器视场信息本身就难度较大。与之相比，基于聚类的方法原理上无须视场相关的先验信息，可以进行视场未知的分布式融合，但是聚类方法对高斯分量/粒子的聚集程度要求较高，容易造成目标信息错误的归类，此外聚类算法还会引入较高的计算复杂性。

本文基于先进的CPHD滤波器框架，开展了部分重叠视场下的分布式雷达多目标跟踪研究。首先，文章利用多目标密度乘积来提取雷达共同观测目标信息的PHD函数，并采用多伯努利近似及解卷积方法计算切分的势分布。其次，对于共同观测目标信息采用标准的分布式CPHD滤波器方法进行融合，剩余部分则以补偿融合的方式进行合并。最终，利用势卷积及PHD合并得到雷达网络全视场下的多目标信息。本文方法无须视场先验信息，因此也可以适应雷达视场未知的分布式融合多目标跟踪场景。实验结果表明，本文所提出算法能够在时变未知雷达视场下较好的跟踪目标，且相比基于聚类的高斯匹配方法性能更好。

文章其余部分结构如下：第2节对问题背景及相关理论知识进行了介绍；第3节提出了基于多目标密度乘积切分共同视场信息的方法，并实现了雷达视场部分重叠时的分布式融合多目标跟踪；第4节利用仿真对所提出的方法进行了验证；第5节为结论。

2 问题背景

本节首先对随机集及CPHD滤波器进行简单阐述，然后介绍CPHD滤波器下的分布式融合方法。

2.1 随机集和CPHD滤波器

在空间 X中，随机集X表示F(X)中的一个随机变量，其中F(X)是 X所有有限子集的集合。随机集是一个元素互异、无序，且数目可变的有限集合，它可以很自然地表征多目标演化过程中的虚警、漏检、目标新生以及消亡。将多目标运动过程及量测过程建模为独立同分布群(Identical Independently Distributed Cluster,IIDC)随机集，则可以得到CPHD滤波器。IIDC随机集X的多目标密度以及PHD函数分别表示为

其中，s(x)表 示单目标空间分布，p(n)表示势分布。

相比于仅传递PHD函数的PHD滤波器，CPHD滤波器迭代传递PHD函数和势分布(指代目标数目的分布)，使其在目标数目估计的准确性和稳定性上拥有更加优越的性能。CPHD滤波器的迭代方程较为复杂，此处不再列出，其具体形式及相应的GM实现方法可查阅文献[27]。

2.2 基于KL平均的分布式融合

分布式随机集滤波器融合算法主要有GA/AA融合方法。由于它们均可通过最小化融合密度与融合结果之间的KL (Kullback-Leibler) 散度权重和得到，则都可视为KL平均融合。相应的CPHD滤波器的分布式融合形式以如下定理给出[15]。

定理1(AA-IIDC)给定来自传感器i∈Nj的IIDC分布表示传感器j的邻居传感器(包含j)，它们的AA融合的形式为且PH D和势分布分别为

其中，wi表示指向传感器i的权重且但是由于fAA(X) 不再是一个IIDC分布，无法在下一次CPHD滤波器迭代中使用。因此寻找一个与fAA(X)有 相同PHD及势分布的IIDC分布作为近似融合结果反馈至本地传感器，则

定理2(GA-IIDC)给定来自传感器i ∈Nj的IIDC分布表示传感器j的邻居传感器(包含j)，它们的GA融合仍然是IIDC分布，且形式为

根据以上可知，IIDC分布(对应CPHD滤波器)在GA融合下仍能得到IIDC分布，即该分布在GA融合下是闭合的。闭合性使得其无须寻找近似结果，能够得到精确的表达式。相比之下，(标签)多伯努利滤波器在分布式融合中的没有闭合性，需要经过一些启发式的方法来进行计算[15]。

以上分布式融合方法均假设不同传感器的视场相同。但在分布式雷达中各传感器视场不一致时，相应的PHD和势分布表征的信息也不一致，以上标准的AA/GA融合算法不适用于此场景。下面将针对分布式雷达视场部分重叠的情况，为GM-CPHD滤波器设计分布式融合算法。

3 雷达视场部分重叠的分布式融合

本节通过不同雷达多目标密度的乘积来确定视场间共同区域的非零分量，并通过切分、融合以及合并的步骤实现雷达视场部分重叠下的分布式融合。整体流程如表1所示。

表1 部分重叠视场下分布式雷达多目标跟踪算法Tab.1 Distributed multitarget tracking using radars with partially overlapping FoVs

3.1 PHD函数切分

假定雷达i和j在各自视场Fi和Fj(区域所对应状态空间分别为χi和χj)下使用CPHD滤波器，分别得到势 分布pi(n)和pj(n)，以及PHDDi(x)和Dj(x)。PHDD(x)是 单目标x的普通函数，因此其实际上也包含了雷达的空间覆盖信息。又PHD对空间外的目标密度无定义，这里首先对其进行扩展，即

同样也对Dj(x)进 行扩展得到(x)。由于增加的是零值，则扩展不影响势分布。两雷达的PHD的乘积为

假定雷达能够精确地获得其观测区域下的PHD，则必有

式(13)说明：在理想情况下，扩展PHD的乘积若不为零，则相应状态必位于两雷达视场的重叠部分。这一性质为利用乘积寻找共同观测目标信息提供了可能，即可以将乘积(x)大于零的部分视为两雷达的共同观测目标信息，并将其从原始密度中切分。

以上过程是在两个雷达完全无偏且对目标的检测估计完全一致的情况下进行的。但由于雷达的量测噪声，以及在检测过程中的漏检、虚警等因素，精确的点乘(如粒子滤波)是难以实现的，因此，本文将开展部分重叠视场下的高斯混合实现。为方便表述，以下均假设PHD已经过如式(11)的扩展，即不再区分(x)和D(x)。

对于雷达i，PHD的高斯混合形式可以表示为

其中，N(·;m,P)表 示均值为m协 方差为P的高斯分布。同理，雷达j的高斯混合参数为可将该PHD写成共同区域目标和单独区域目标之和，即

其中，两式中前一个高斯和表示两雷达共同区域下的目标估计，后一个高斯和表示仅在各自视场内的目标估计或者由杂波引起的分量。

根据高斯分布的乘积性质：两个高斯分布的乘积依然为一高斯分布，即

可以看出，当高斯分布的均值m1,m2距离较近时，乘积权重w12将 趋近于w1w2；当高斯分布的均值m1与m2距 离较远时，则乘积权重w12将趋近于0。因此，两雷达的高斯混合的乘积可近似为一高斯混合：

式(21)近似成立的的依据是：如式(17)，只有当两个高斯分量足够近，才能使得它们乘积的权重系数足够大。而对于两个高斯分量相距较远时，它们的乘积的权重系数将非常小。映射到视场空间下，若两雷达均能够检测到一目标，则它们对该目标的估计应当距离“较近”。反之，对于仅位于雷达i视场下的目标，由于雷达j未检测到它，也无法形成该目标的跟踪轨迹，则相应的高斯分量将会离雷达j中的所有高斯分量都比较远。对于杂波产生的高斯分量，其本身在空间和时间上没有相关性，因此杂波形成的分量和目标形成的分量一般也相距较远。

为了实现如式(21)的功能，需要滤除两高斯混合乘积中权重较小的分量，可以采用门限Γ作用于所有计算的权重，若Γ则保留相应的权重及其对应的高斯分量，否则丢弃该高斯分量，这类似于高斯混合的修剪步骤。

3.2 势分布切分

获得切分的共同观测目标信息PHDDi,I(x)之后，将权重进行以下处理，

其中，m od(·)是 求余操作表 示取小于x的最大整数。处理后的权重可以表示为集合：

以上处理保证了所有权重均不大于1，则此时可以将以上关于区域Fi,I的目标随机集近似为一个存在概率为vm ∈V的多伯努利随机集[20]，且对应的势分布可以计算为

同理，利用PHDDj,I(x)计 算势分布pj,I(n)。原理上也可以利用Di,O(x),Dj,O(x)以及多伯努利近似来计算剩余的势分布pi,O(n)和pj,O(n)。但是，这样就忽略了原始势分布pi(n)和pj(n)的信息，所有的过程仅对PHD操作，相当于PHD滤波器在有限视场下的融合，丧失了CPHD滤波器的优势。为此，这里利用总的势分布pi(n)和 切分的势分布pi,I(n)来计算雷达i剩余部分对应的势分布pi,O(n)。根据文献[20]的证明，总的势分布可以分解为若干个势分布的卷积，即

其中，*表示卷积。则在已知总势分布pi(n)和子势分布pi,I(n)的情况下，求取另一个势分布需要求解卷积操作。为此，基于数字信号处理中常用的卷积性质[29]，pi,O(n)可计算如下：

F(·)表 示离散傅里叶变换，F(-1)(·) 表示离散傅里叶逆变换。类似地，传感j非共同区域的势分布pj,O(n)也可以按照以上方法进行计算。由此，计算出了所切分的势分布pi,I(n),pi,O(n),pj,I(n)以 及pj,O(n)。

3.3 PHD及势分布的合并

获得以上切分的PHD及势分布后，再进行融合步骤。根据以上分析，Di,I(x)(pi,I(n)) 和Dj,I(x)·(pj,I(n))表示的是两个雷达共同观测的多目标密度信息，则按照标准的GA/AA融合方法(定理1、定理2)对它们融合可以得到PHDDI(x) 和势分布pI(n)。之后，计算全视场的势分布及PHD分别为

4 仿真及实验

本节仿真中利用3个运动的雷达平台观测4个运动的目标，仿真场景如图2所示。图中直线运动的目标轨迹由黑色粗实线画出，做小角度转弯运动的雷达轨迹由红色实线表示。图中显示了不同时刻k雷达视场覆盖目标的情况，雷达s1—s3的视场分别由红、绿和蓝色阴影部分表示。雷达传感器s1,s2和s3的波束中心在每一时刻分别指向目标t1,t2和t4，且波束宽度分别为30°、20°和30°。由于目标在雷达视场内的进出相当于目标的出现和消失，因此这里不再单独设置目标出现消失，即目标在整个监督区域中始终存在。

图2 不同时刻的雷达视场示意图(s1—s3为雷达平台，t1—t4为目标，k表示不同时刻，坐标轴单位为km)Fig.2 The illustration of the sensor field of views in different time(s1—s3 are radar platforms,t1—t4 are targets,k represents time index,and the unit of coordinate is km)

目标t1—t4的初始位置分别为[10,40] km,[13,35] km,[16,40] km和[20,35] km。雷达s1—s3的初始位置分别为[0,5] km,[12,1] km和[27,5] km。雷达使用GM-CPHD滤波器来跟踪多目标。本地雷达中以匀转弯运动模型来跟踪目标，运动状态表示为xk=[xp,k,xv,k,yp,k,yv,k,ωk]T，其中下标p和v分别表示位置和速度。状态转移方程表示为

雷达可以测得各目标在雷达坐标系下的方位角和斜距：

其中，vr,t和vφ,t是相互独立的零均值高斯白噪声且标准差分别为表示雷达i在k时刻的位置。

目前在随机集领域能够处理时变未知视场的分布式融合算法还非常少。基于聚类的方法原理上无须视场相关的先验信息，可以实现未知有限视场下的分布式融合。因此，本文仿真将所提出的算法与文献[21]中高斯匹配的方法进行对比，它对各传感器得到的GM-PHD利用马氏距离匹配，并对距离较近的高斯分量配对分别进行GA融合。为了使该方法能应用到CPHD滤波器中，在得到融合PHD后，利用多伯努利近似计算势分布。则仿真中对比了3种算法：“AA+本文算法”表示融合步骤使用AA融合的本文算法；“GA+本文算法”表示融合步骤使用GA融合的本文算法；“高斯匹配方法”表示文献[21]的算法。

仿真各项参数设置为：仿真结果由执行200次蒙特卡罗实验获得。目标的存活概率设为pS=0.90。各雷达的目标检测概率设置为pD=0.97。每个雷达检测到的杂波服从均值为λ=10的泊松分布，且在各自视场内均匀分布。目标新生过程服从均值为0.001的泊松分布，且新生目标的PHD由量测驱动的自适应方法生成[30,31]。高斯混合中修剪和合并的门限分别设置为 1×10-4和2 0 m。最终的性能由阶数p=2、截断因子c=400 m的OSPA误差[32]评估。分布式GA/AA融合采用等权重进行融合。为了防止融合后全视场的信息影响下一时刻基于乘积的信息提取，本仿真中融合的结果不反馈至本地滤波器，这一不反馈策略在文献[26]中也被使用。最终仿真结果如图3所示。

如图3(a)所示，本文提出的融合算法可以较为准确地估计目标数目，而高斯匹配方法会在一定程度上低估目标数目。引起这种现象的原因在于高斯匹配方法的最终融合仅在每两个权重高斯分量之间进行，这相当于进行了多组目标数为1的泊松密度的融合。根据GA融合性质可知，当目标数目为1时，泊松分布的GA融合存在估计下偏置[15]。因此，仿真中多个匹配对的GA融合会进一步加剧目标数目的低估现象。

如图3(b)所示，本文提出的算法在最终性能上明显优于高斯匹配的方法。而且根据如图3(c)和图3(d)的对OSPA度量的分解可知，这种优势主要来源于对目标数目的准确估计。

由图3(a)—图3(d)可以看出，在k=17,33,52,58s时刻，由于目标穿越雷达视场，目标在穿出视场下存在消亡现象，这使得所有融合方法的估计都会出现尖峰。相比之下，本文算法出现这种跳跃尖峰的值更大，这是由于未反馈融合结果所导致的，但所有尖峰会在很短的时间内消失。

图3 目标数目估计及OSPA误差Fig.3 Estimated number of targets and OSPA error

5 结论

本文基于GM实现的CPHD滤波器，在雷达视场部分重叠情况下实现了分布式雷达多目标跟踪。文章提出了使用多目标密度乘积来提取雷达共同观测目标的PHD函数，并通过多伯努利近似及卷积性质切分势分布的方法。共同观测目标信息仅通过简单的乘法和高斯修剪操作即可实现，且拥有较好的性能。仿真实验表明，所提出方法比高斯匹配方法拥有更好的性能。