基于PHD滤波的雷达起伏目标检测前跟踪算法研究

吴孙勇，李东升，薛秋条，孙希延，蔡如华

1 引言

传统雷达目标跟踪算法［1］是对雷达天线扫描的单帧数据进行的，为了保证检测过程中较低的虚警率，通常采用较高的检测门限得到目标点迹信息后再进行目标航迹跟踪，而弱目标的回波通常低于门限，这可能会导致目标丢失，发生漏检；如果设置较低的门限，则会使得虚警率增高，目标的航迹无法维持. 为解决上述问题，方法之一是采用检测前跟踪（Track-Before-Detect，TBD）［2，4］算法. 该算法根据空间中目标运动的连续性和连续几帧目标回波数据时间上的关联性，对多帧数据进行联合处理，并通过多帧能量累积实现目标检测和跟踪. 传统的TBD 方法包括动态规划［5，6］、Hough 变换［7，8］等. 但这些方法一方面计算量大、算法复杂度高，另一方面只适用于线性高斯模型. 贝叶斯框架下的粒子滤波［9，12］TBD（Particle Filter TBD，PF-TBD）算法可以解决非线性和非高斯的问题，从而在微弱目标检测跟踪领域得到了快速发展.

PF-TBD方法的局限在于没有对目标的出现（新生）和消失（死亡）建模，因此在处理目标数目未知且变化的时候会导致算法复杂度急剧增大，滤波器性能存在一定局限性. 基于随机有限集（Random Finite Set，RFS）理论［13］的多目标跟踪算法可以解决目标时变，数据关联不确定性和检测不确定性的问题. 概率假设密度（Proba‑bility Hypothesis Density，PHD）滤波器［14］是一种基于随机有限集理论的多目标贝叶斯滤波算法，PHD滤波表示的是多目标后验概率密度的一阶统计矩，在递归过程中传递PHD，PHD 滤波具有计算复杂度低，实现简单等优点.PHD 滤波器一经提出就成为了研究热点，研究表明PHD 滤波器尤其适用于杂波密集且目标随时间变化的多目标场景；VO 等人给出了PHD 滤波器的序贯蒙特卡洛（Sequence Monte Carlo，SMC）实现［15］和高斯混合（Gaussian Mixture，GM）实现［16］. 文献［17］首次将SMCPHD滤波器用于红外图像多目标TBD中，利用序贯性蒙特卡洛技术实现对未知目标数的弱小多目标检测跟踪；文献［18］提出基于平滑的PHD-TBD 算法，减少了PHDTBD算法的目标数估计起伏；文献［19］解决了当目标数过多时PHD-TBD 杂波数不服从泊松分布的假设. 文献［20］对噪声进行泊松化，设计出一种能解决多目标TBD问题的PHD滤波器. 目前PHD-TBD算法常常是处理红外图像多目标检测跟踪，对于将PHD 理论运用在雷达TBD领域还存在诸多困难之处.

雷达量测方程中，来自目标的电磁能量散射的复杂性和未知时变的因素，目标的振幅可能在扫描之间突然改变，振幅起伏由均匀分布的相位和Swerling 模型模拟的幅度决定［12］. 对于幅度起伏目标数目不发生改变的检测跟踪常见的处理方式是在PF-TBD 下进行的，但由于振幅参数是未知的，似然函数并不能直接从量测方程中获得；求解未知振幅似然常见的解决方案［9，21，22］是对复量测进行取模处理；对于这种彻底抛弃相位的解决方法，有两种策略处理幅度起伏，一种是对幅度起伏密度的整体似然函数边缘化［21］；另一种是独立边缘化每个位置的似然［22］，后者可以求得Swerling 0，1，3 类型的封闭形式解［23］. 无论哪种方法，都彻底丢弃相位的空间相干性，文献［24］表明相位的丢失会导致滤波器灵敏度的损失. 为了避免这种损失，Davey 等人［24］提出在单目标Swerling 0 模型中，直接边缘化整个复量测的似然.Lepoutre等人［25］在其基础上推导了单目标Swerling 1，3 型复似然的封闭形式；对于目标数目不发生改变的多目标跟踪情况，给出了Swerling 1 下封闭形式，及Swerling 0，3 类型的近似解以降低计算成本.李渝等人利用复似然比改进粒子滤波下目标的实时检测前跟踪［26］. 对于幅度起伏目标数目发生改变的PHDTBD 算法通常是在平方量测下进行常幅值目标的检测与跟踪［27，29］，基于平方量测似然的Swerling 0 模型需多次进行贝塞尔函数运算，这会增加算法的复杂度. 平方量测忽略了相位信息会直接影响目标跟踪性能；在强度不同的目标同时存在且距离较近时，PHD-TBD 算法对较弱目标的检测会受到较强目标的影响，使得目标数目估计错误，目标航迹无法维持；传统的PHD-TBD滤波假定新生目标先验分布信息是已知的，而在真实情况下，目标可能出现在场景中的任何位置，此时传统的PHD-TBD滤波器将不再适用.

针对上述问题，本文基于文献［25］的思想提出了雷达起伏目标PHD-TBD 算法，检测和跟踪低信噪比下的多个不同慢起伏（Swerling 0，1，3）的微弱目标，对于PHD-TBD 算法实现过程中，最重要的是似然函数的计算，在以往情况下，假设似然函数服从高斯分布即可，但是目标起伏过程中需要对未知参数边缘化，所获得的似然函数不再是一般的高斯分布；通过边缘化处理未知时变的振幅信息，给出幅度起伏下PHD 滤波的SMC两种不同似然比的实现. 在起伏多目标中，强度较弱目标的检测会受到较强目标的影响，设置阈值使其在同一强度下进行估计；当起伏目标似然比出现无穷大时，为维持目标航迹跟踪，将预测先验信息代替更新后验. 对于场景内新生目标先验信息未知的问题，将量测区域进行划分，利用前一时刻的量测似然比来新生目标. 仿真实验表明，三类起伏目标在复似然比PHDTBD 算法下的检测和跟踪性能比幅度似然比有明显的提高，且运行速度更快.

2 模型与背景

2.1 随机有限集下的多目标跟踪问题

有N(k)个目标xk,1,…xk,N(k)，多目标状态Xk表示为RFS［13］：

Xk={xk,1,xk,2,…,xk,N(k)}∊F(X) （1）

其中，F(X)为X的空间有限子集. 多目标状态集Xk包括目标存活RFS，新生RFS以及衍生RFS.k时刻目标运动方程为：

xk,i=fk-1(xk-1,i) +vk-1,i（2）

假定在时刻k，在状态空间

Zk={zk,1,zk,2,…,zk,M(k)}∊F(Z) （3）

其中，F(Z)为Z的空间有限子集.

2.2 随机有限集下多目标贝叶斯滤波及PHD滤波

假设k-1 时刻多目标后验fk-1|k-1(Xk-1|Z1：k-1)已知，给定k时刻的新量测Zk及历史量测Z1：k-1={Z1,Z2,…,Zk-1}，则多目标贝叶斯滤波器预测和更新步骤如下［13］：

预测：

更新：

其中πk|k-1(Xk|Xk-1) 表示多目标转移概率密度，pk(Zk|Xk)表示多目标获取传感器量测的多目标似然函数，且：

是状态空间X的有限子集F(X)上函数f的集积分.

尽管多目标贝叶斯递归公式中涉及集积分等计算量大的问题，但基于RFS 的PHD 滤波算法通过传递多目标后验概率密度分布的一阶矩，将多目标状态空间简化为单目标状态空间，很大程度上减少计算复杂度，则PHD滤波的预测和更新方程如下［13］：

其中b(∙)和γ(∙)表示k时刻目标新生和衍生的PHD；κk(z)是虚警强度，pD,k(∙)是目标状态的检测概率.

在式（8）中，pD,k(∙)是与目标状态相关的检测概率；对于TBD 算法而言，在更新步骤之前没有检测过程，假设量测值包含所有的目标信息；因此pD,k(x) ≡1，所以更新式（8）的方程变为：

另外，TBD 算法中，目标不会影响场景内所有像素，直接使用p( |ZkXk)会导致计算效率很低；定义似然比为目标存在似然p1( |ZkXk)与目标不存在似然p0(Zk)的比值［24］：

故对于TBD 算法，式（8）中的似然利用式（10）的似然比来替代.

2.3 雷达量测模型

传统的雷达目标检测跟踪问题中，首先对每一帧回波数据进行阈值处理以形成点迹信息，然后对超过阈值的点迹信息进行关联、滤波等处理，最后获得目标的航迹，从而实现目标的跟踪. 对于这种先检测后跟踪的方法，适用于信噪比较高或目标回波幅度较大的情况. 如图1（a）所示，目标的强度远远高于杂波的强度，通过设置较大的阈值可以将目标与杂波分离. 对于低信噪比或目标回波信号较弱的情况下，目标回波湮没在噪声杂波中，如图1（b）所示，采用单帧过门限的检测方法，门限过高会造成目标漏检，而门限过低会造成虚警率增高，目标航迹无法维持.

雷达传感器接收的量测值由距离匹配滤波和自适应波束形成后的距离和方位组成. 给定目标状态xk,i，量测值z,k(l,m)由以下非线性方程给出［25］：

其中，h()l,m(xk,i)表示以目标位置(xk,i,yk,i)为中心的第i个目标的模糊函数，nk是均值为0，协方差为复高斯Γ的量测噪声；φk,i是[0,2π)上均匀分布的相位，ρk,i是由Swerling模型表示的幅度：

（1）在Swerling 0中，ρk,i无起伏，等于未知参数ρi：

pSw0(ρk,i)=ρi（12）

（2）对于Swerling 1，幅度ρk,i服从瑞利分布，且

（3）起伏类型为Swerling 3，ρk,i服从四个自由度的

图1 不同信噪比下目标回波

从而k时刻的量测可以表示为RFS：Zk=1,…,Nr,m=1,…,Nθ}. 其中，Nr为距离分辨单元的个数，Nθ为方位分辨单元的个数.

雷达在极坐标中定义区域[rmin,rmax]×[θmin,θmax].考虑对于距离，假定发射的脉冲是带宽B和持续时间Te的线性调频信号，距离模糊函数由文献［20］给出：

其中，|τl|=2(rk-r)l/c，c是电磁波速度；rk=

在接收端，雷达由线性相控阵组成，其中Na根天线的间距是是载波频率的波长. 则方位模糊函数由文献［25］给出：

其中：表示角度分辨率

距离方位单元(l,m)中的总体模糊函数由乘积

3 不同幅度起伏下似然比的计算

3.1 幅度起伏情况下幅度似然比的计算

量测方程式（11）不仅与目标扩散有关，还取决于未知参数相位φk,i和幅度ρk,i，故不能直接计算目标存在的似然p(Zk|xk). 常用的策略是考虑复量测的模平方量|Zk|2，噪声协方差矩阵为Γ=2σ2INc，INc是Nc维的单位方阵；

其中，I0(∙)是第一类修正的贝塞尔函数.γl为非中心参数：

对于目标不存在时的似然p(|Zk|2)有：

对式（17）积分边缘化幅度和相位后，除以式（19）即可得幅度似然比. 下面给出不同幅度起伏类型下目标存在的幅度似然比.

对于Swerling 0：

在Swerling 1的情况下，修正文献［25］可得

当目标幅度起伏类型是Swerling 3 时，结合文献［23］可推导出Swerling 3幅度似然比为：

3.2 幅度起伏情况下复似然比的计算

利用量测平方计算似然比并不是最佳的，因为它没有考虑复振幅的空间相干性以及该方法需假设噪声协方差为Γ=2σ2INc. 为了避免这些缺点，文献［24，25］考虑了相位的空间相干性，并且可以处理空间相关噪声. 复似然比的计算不仅提高了目标检测跟踪的性能，而且有效降低了计算复杂度.

通过量测方程式（11）可计算p(Zk|Xk,ρk,1：Nk,φk,1：Nk)［25］：

其中

对于目标不存在时：

似然比L(Zk|xk,ρk,1：Nk,φk,1：Nk)由式（23）除以式（24）可得；再通过对参数ρk,1：Nk,φk,1：Nk上边缘化，有：

在Swerling 0中：

其中，ρk是常数.

当目标幅度起伏类型是Swerling 1时：

对于Swerling 3：

4 算法实现

4.1 自适应新生幅度起伏PHD-TBD粒子滤波实现

对于k-1 时刻，用一组带权重的粒子代表PHD的后验密度，即：

4.1.1 预测

预测的粒子：

其中，

ps,k(是状态为的目标从时间k-1 到时间k的存活概率；sk|k-1()是状态为xk-1的衍生概率密度；是新生概率密度.

4.1.2 自适应新生

一般来说，PHD 滤波器的标准实现假设新生分布γk(x(i)k)是已知的，但当先验知识未知或低信噪比环境下，先验信息已知的PHD-TBD算法并不适用，PHD-TBD可能由于连续的误检而不能初始化目标. 通常的策略是将新生分布覆盖整个状态空间，然而这种方式需要大量的粒子来表示具有SMC实现的新生模型，这是极其低效的，尤其对于场景很大的时候［30］. 为克服新生先验的不足，提出一种基于划分量测区域似然比的自适应新生分布.

首先，对于目标跟踪区域进行划分，将目标场景等分为几个小场景；其次在不同的小场景上进行量测挑选；然后对于每个挑选的量测在其附近生成目标状态的位置信息，由于本文讨论的传感器接收的是目标的距离和方位信息，故目标状态的速度信息在一定的范围内均匀生成；在低信噪比的情况下，目标信号较强的位置不一定是真实值产生的，为降低杂波对初始粒子选取的影响，在每次新生粒子后都进行似然比的计算，即赋予每个粒子权重；最后重采样挑选出真实目标附近的粒子. 记自适应新生分布为这里给出其算法步骤：

4.1.3 更新

PHD的更新用SMC可以表示为：

其中，Zkxk)不一样. 注意，由量测方程可知，目标的回波信号为sinc 函数，此时的目标存在位置的似然比会很高；又当目标幅度起伏时，回波信号的强度各不相同，为解决这种情况下的多目标的跟踪问题，文献［31］提出对数似然比的方法将目标存在位置的似然比降低. 在本文中，为降低算法的复杂度，同时解决不同强度的目标匹配跟踪，将似然比大于阈值ThL的都赋值为最高似然比.

利用K-means 方法对重采样后的粒子进行聚类划分，考虑到本文低信噪比和幅度起伏的情况，每次聚类后删除权重和低于门限Thk‑means的粒子群；将权重和高于门限Thk‑means的粒子群视为估计的目标状态.

对于不同幅度起伏类型和似然比，LSw( |

4.2 幅度起伏PHD-TBD SMC实现

算法2 给出幅度起伏PHD-TBD SMC 实现的伪代码.

5 仿真实验及结果分析

为验证本文算法的有效性，我们将评估不同PHDTBD 策略对Swerling 0，1，3 类型目标的检测性能. 考虑一个二维运动场景，每个目标的状态定义为xk=[xk,ẋk,yk,ẏk]T，其中(xk,yk)和(ẋk,ẏk)分别是笛卡尔坐标系下目标的位置和速度；位置(xk,yk)在极坐标p=[rmin,rmax]×[θmin,θmax]场景内，rmin,rmax和θmin,θmax分别为最小和最大目标范围和方位；目标状态按匀速直线运动演化：

xk=Fxk-1+vk-1（36）

传感器扫描时间间隔T=1 s接收100帧图像，其中：

过程噪声vk-1服从高斯分布，其协方差为

噪声标准差σv=5 m，目标存活的概率为ps,k(x)=0.98.

对于目标速度(ẋk,ẏk)在区域C上：

其中，vmin=200 m/s，vmax=800 m/s.

表3给出所需雷达和模拟场景具体的参数.

表3 模拟数值

评估算法性能采用最优子模式分配距离［33］（Opti‑mal Sub Pattern Assignment，OSPA），OSPA 度量可以评估多目标滤波器的目标个数估计误差和目标位置估计误差，给定两个有限集X={x1,x2,…,xm} 和Y={y1,y2,…,yn}，OSPA定义如下：

其中 ，dC(X,Y)=min{C,db(X,Y)}，db(X,Y)=p，C>0为截断参数，用于惩罚目标个数的估计偏差，p为阶数，用于惩罚多目标状态估计偏差.

在本文仿真实验中，设置p=3，C=1000. OSPA 值越小，表明目标数目和状态估计越准确. 假设总共有五个目标，初始状态见表4，其中图2 为目标真实运行轨迹图.

表4 目标初始状态

图2 目标运动的真实轨迹

5.1 同一信噪比下不同幅度起伏的自适应新生算法性能仿真

在SNR=7 dB，每帧新生粒子数2000，利用50 次蒙特卡洛仿真，在不同幅度起伏下，对比算法仿真. 结果如图3~5所示.

图3 Swerling 0下复似然比和幅度似然比算法对比图

图3（a）对比了Swerling 0 下幅度似然比和复似然比对应的OSPA，由图可知，从目标出现到目标消失，复似然比的OSPA 比幅度似然比低. 在新生目标出现的时候，OSPA 会突然增大，这是因为自适应新生算法是利用前一时刻量测导致的，具有一定的延后性. 由图3（b）可知，两种方法下的平均势估计比较接近. 复似然比与幅度似然比在常幅值下的效果差异不是很明显，但是复似然比平均运行一次的时间为139 s，幅度似然比平均运行一次的时间为1281 s，幅度似然比的运行时间差不多是复似然比的9.22倍.

图4（a）表示Swerling 1 下PHD 滤波下复似然比与幅度似然比的OSPA 对比图，由图可知复似然比对目标位置信息的估计比幅度似然比效果好，由图4（b）可知，造成幅度似然比OSPA 过高的原因主要是对于目标势的估计，幅度似然比在起始过程中会估计较多的错误位置信息，且这些误检大多在雷达边缘，幅度似然比损失的目标相位信息导致在估计时刻无法达到目标最大势；在Swerling 1 下的幅度似然比估计效果不佳. 复似然比平均运行一次的时间为132 s，而幅度似然比平均运行一次的时间为347 s.

图5表示，在Swerling 3条件下，PHD-TBD复似然比对比幅度似然比损失最小.Swerling 3 描述的是由很多较弱的散射体和一个特别强的散射体构成的目标特性. 对于Swerling 3来说，如果粒子权重和无穷大，则保持粒子状态不变，用一个尽可能小的数代替粒子权重，这样做虽然解决了粒子和无穷大而导致目标航迹无法延续的问题，但舍弃了目标状态的更新，对于目标数目的估计以及状态的估计都会造成影响. PHD 下Swer‑ling 3复似然比平均运行一次的时间为177 s，而幅度似然比平均运行一次的时间为452 s.

图4 Swerling 1下复似然比和幅度似然比算法对比图

图5 Swerling 3下复似然比和幅度似然比算法对比图

综合图2~图5，PHD-TBD 幅度起伏算法需要经过多帧跟踪估计，才能检测到目标；同时前一时刻的自适应新生，目标有延后效应，因此目标在首次出现的时候，OSPA 会增大，但随着多帧积累，PHD-TBD 幅度起伏能实现对目标数目的估计. 在信噪比是7dB 的情况下，基于PHD 幅度起伏的检测前跟踪算法复似然比方法比幅度似然比方法更有优势，且复似然比比幅度似然比方法能更快地估计到目标.

5.2 不同信噪比下幅度起伏的自适应新生算法性能仿真

在SNR=7 dB和SNR=5 dB，每帧新生粒子数2000，在不同幅度起伏下，对比算法仿真. 结果如图6~8所示.

图6 Swerling 0类型不同信噪比下复似然比和幅度似然比算法对比图

由图6 可知，Swerling 0 下信噪比为5 dB 的幅度似然比是几乎完全无法正确的估计，而5 dB 的复似然比跟踪性能良好；复似然比的跟踪效果在信噪比降低的情况所受影响不大.

由图7 可知，起伏类型是Swerling 1 的情况下，降低信噪比会导致幅度似然比的目标个数估计在初始时刻增大，而相比于复似然比，目标的个数估计与目标位置的估计对于整个时刻而言所受影响不大.

由图8 可知，Swerling 3 下复似然比目标个数估计会随着信噪比的降低而导致目标个数估计不准确，会产生一定的漏检，但相比于幅度似然比所引起的误差，复似然比对于信噪比降低产生的影响是较小的.

图7 Swerling 1类型不同信噪比下复似然比和幅度似然比算法对比图

图8 Swerling 3类型不同信噪比下复似然比和幅度似然比算法对比图

目标姿态相对于雷达发生变化时，雷达到目标的距离和时间也将发生变化，运动目标雷达横截面（Ra‑dar Cross Section，RCS）随着时间推移而起伏，RCS 的起伏会降低检测概率，RCS 表示目标返回到雷达的回波信号的幅度. 而目标幅度起伏流行的表示方法是Swer‑ling 描述的统计模型，对于Swerling 0，RCS 是一个常数，典型类型如金属圆球；对于Swerling 1，目标RCS 在一次天线波束扫描期间时完全相干，扫描间的脉冲不相干，典型的类型如前向观察的小型喷气飞机等；对于Swerling 3，情况与Swerling 1 一样，只是概率密度不同，典型的类型如螺旋桨推进飞机，直升机等. Swerling 2和Swerling 4 类型分别考虑了与Swerling 1 和Swerling 3相同的概率密度，但起伏速度较快，在本文不予考虑.在实际情况下，运动目标姿态随时会发生改变，本文考虑不同幅度起伏类型，模拟不同情况下的多目标运动，并给出算法实现，体现了PHD-TBD 滤波中考虑相位信息的复似然比优势，而对于先验信息未知时，所提自适应新生算法也能很好地跟踪迭代.

6 结论

本文对雷达幅度起伏弱小目标问题进行了深入分析和研究，将PHD 滤波的思想应用在检测前跟踪上，给出PHD-TBD幅度起伏两种似然比的序贯性蒙特卡洛算法实现，完善了基于PHD 滤波的检测前跟踪算法；同时利用贝叶斯思想，考虑先验未知的情形，提出一种量测划分区域自适应新生的算法. 针对幅度起伏在PHD 滤波下可能出现的问题也做了进一步详细的说明. 仿真实验表明，在PHD-TBD目标幅度起伏的情形下，复似然比比幅度似然比能更准确更快地发现和估计目标.