“低慢小”目标的探测是一个世界性难题，困难主要源于以下3个方面^[1-2]：①飞行高度“低”，一般在1000m以下，地面雷达观测时会有大量杂波进入接收机，尤其在建筑密集的城市环境中，探测难度更大。②飞行速度“慢”，一般小于200km/h，其回波信号处于杂波主瓣区，容易与鸟群等低速杂波混淆。为防止过多虚警，雷达往往设置速度阈值，导致低速的低空小目标被滤除。③雷达散射截面积(RCS)“小”，目标信噪比低，雷达反射信号忽隐忽现，目标不易被发现与识别。随着无人机技术的日渐成熟，其作为一种典型的“低慢小”目标，严重威胁航班安全运行^[3-4]。

雷达系统通常包括检测和跟踪2个部分。其中，检测部分将目标从噪声和杂波环境中提取出来，获得目标的速度、位置等状态信息；跟踪部分对目标下一刻的状态进行预测，进而形成稳定航迹。在常规雷达系统中，检测与跟踪部分通常分别考虑和设计^[5-6]。实际上，二者是相互影响的。检测是跟踪的前提，良好的检测是保证跟踪性能的基础；同时，跟踪能进一步完善检测结果，利用跟踪获得的目标动态特征能够提高检测能力^[7-8]。因此，本文提出一种基于Rao-Blackwellized蒙特卡罗数据关联的检测跟踪联合优化算法，对雷达目标检测与跟踪进行交互处理，将跟踪结果反馈给检测部分，提高雷达系统对低空小目标的检测能力，同时保证下一步的精确跟踪。

1 检测跟踪联合优化技术

雷达系统在完成目标检测之后，通常将波门内超过一定阈值的量测数据通过跟踪算法进行关联^[9]。如将跟踪结果反馈至检测中心，进而根据该反馈信息调整波门内的检测阈值，改善检测结果，此类算法即可称为检测跟踪联合优化算法。王云奇^[10]结合目标、环境等先验知识的辅助，量化分析检测性能与跟踪性能的关系，实现检测和跟踪的耦合。刘红亮、严俊坤等^[11-12]提出一种航迹恒虚警的目标检测跟踪一体化算法，根据帧虚警概率调整预测波门内的检测门限，完成目标检测与跟踪的联合优化过程。闫学昭^[13]采用DSP等硬件编程方法，搭建了雷达目标检测与跟踪交互处理模块，将跟踪结果反馈至检测系统，在提高检测准确率的前提下实现了精确跟踪。

综上所述，这些算法的共同点在于基于一定规则降低波门内的检测阈值，以牺牲波门内局部虚警率为代价提升局部目标的检测率^[11]。可见，波门的设置至关重要，如果波门面积(体积)过小，则可能将检测目标排除在外；如果波门面积(体积)过大，则可能引入过多虚警，导致计算机过载。针对此缺点，本文提出一种基于Rao-Blackwellized蒙特卡罗数据关联的检测跟踪联合优化算法。该算法采用序贯蒙特卡罗(Sequential Monte Carlo, SMC)(又称为粒子滤波)方法实现多目标跟踪，由于粒子能够用于表示多种不同的数据关联假设，基于SMC方法的多目标跟踪算法可视为多假设跟踪(Multiple Hypothesis Tracking, MHT)的推广^[14]。此外，Rao-Blackwellizion方法的应用能提高算法的准确率和效率，采用卡尔曼滤波或扩展卡尔曼滤波等方法估计目标状态，将SMC方法仅用于数据关联估计，使得联合后验分布由混合高斯表示，比纯粹的SMC方法具有更小的方差。同时，本文算法能够根据粒子的分布范围确定波门大小，在考虑粒子权重的前提下，利用检测单元与所有粒子的相对位置对检测门限进行修正。

2 基于SMC的检测跟踪联合优化算法

本文提出一种基于Rao-Blackwellized蒙特卡罗数据关联的检测跟踪联合优化算法，包括检测和跟踪2个部分，将粒子滤波跟踪获取的目标状态信息反馈至检测部分，以修正相关检测单元的检测门限，实现检测与跟踪的联合优化，算法流程如图 1所示，详细步骤如下：

步骤1  粒子滤波目标状态预估。

以一组N个粒子表示目标的当前状态{w_k－1⁽ⁱ⁾, c_k－1⁽ⁱ⁾, m_k－1⁽ⁱ⁾, P_k－1⁽ⁱ⁾:i=1, 2, …, N}，基于每个粒子i=1, 2, …, N在k-1时刻的预估值m_k－1⁽ⁱ⁾和协方差P_k－1⁽ⁱ⁾对k时刻的目标状态进行卡尔曼预估:

式中：c_k－1⁽ⁱ⁾为数据关联指数；A_k－1和Q_k－1分别为k-1时刻的转换矩阵和过程噪声矩阵；m_k|k－1⁽ⁱ⁾和P_k|k－1⁽ⁱ⁾分别为粒子i在k时刻的预估状态和协方差。

计算每个粒子新的非归一化权重为

式中：w_k－1⁽ⁱ⁾为粒子权重；y_k为k时刻的量测; y_1:k为1~k时刻的所有量测；π(m_k|m_k－1⁽ⁱ⁾, y_1:k)为前一时刻的重要性分布；p(m_k⁽ⁱ⁾|m_k－1⁽ⁱ⁾)为目标的动态模型，描述了目标的随机动态；p(y_k|m_k⁽ⁱ⁾)为测量模型，描述了基于目标当前状态的测量值分布。

进行权重归一化处理：

步骤2   检测单元相似度计算。

基于高斯概率分布函数计算检测单元u_k与每个粒子预估位置的空间相似度：

式中：H_k和R_k分别为k时刻的量测模型矩阵和量测噪声矩阵；N(·)表示检测单元与粒子预估位置的相似度。

结合每个粒子的权重，得到每个检测单元与所有粒子群预估位置的空间相似度为

步骤3  检测门限修正。

利用步骤2计算所得的每个检测单元与所有粒子群预估位置的空间相似度，对每个检测单元的固定检测门限θ进行修正^[15-16]，修正后的检测门限表示为

式中：γ为门限修正强度参数。

该模型的合理性在于：对于与粒子预估位置相距较远的检测单元，其为目标量测的概率较小，则p(u_k)值趋近于0，e^－γp(u_k)值趋近于1，以至于不对阈值进行修正；对于与粒子预估位置相距较近的检测单元，其为目标量测的概率增大，则p(u_k)值增大，且0 < e^－γp(u_k) < 1。参数γ能够调节该指数函数模型的非线性属性，较大的γ值能够强化模型对检测门限的影响; 相反，较小的γ值会弱化以上影响。在目标检测过程中，对于小弱目标，一般会设置较低的检测门限，通过跟踪算法剔除杂波和虚警，并通过本文算法进一步降低预估区域的检测门限；对于回波较强的目标，可适当提高检测门限，以减少虚警数量，提高跟踪效率。

步骤4  目标检测。

根据步骤3计算所得的检测门限，由式(7)判断检测单元中是否存在目标：

式中：rcs(u_k)为检测单元u_k处的量测回波强度值，若该值大于或等于检测门限，则数据确认为目标(E(u_k)=1)，由y_k表示，反之则无目标(E(u_k)=0)。

步骤5  数据关联。

假设马尔可夫模型为m次，则数据关联指数c_k取决于先前的m个关联结果c_k－m:k－1，假设虚警在量测空间V中均匀分布，卡尔曼滤波的量测相似度计算如下：

式中：j=1, 2，…, T，T为目标数目；KF_lh(·)代表卡尔曼滤波器测量相似度估计；H_{j, k}和R_{j, k}分别为目标j的量测模型矩阵和量测噪声矩阵。

对于j=1, 2，…, T，有

式中：KF_p(·)代表卡尔曼滤波的预估部分；m_{j, k－1}⁽ⁱ⁾和P_{j, k－1}⁽ⁱ⁾分别为粒子i目标j的预估值和协方差；A_{j, k－1}和Q_{j, k－1}分别为动态模型转换矩阵和目标j的过程噪声矩阵。

数据关联结果由最优重要性分布决定，其概率分布结果计算如下：

1) 计算非归一化噪声关联概率：

式中：为粒子i的归一化噪声关联概率；p(··)表示条件概率。

2) 为每个目标j=1, 2，…, T计算非归一化目标关联概率：

式中：为粒子i归一化的目标j的关联概率。

3) 归一化重要性分布：

4) 概率代表每个粒子中的离散概率分布，其中数据关联指数c_k⁽ⁱ⁾(i=1, 2, …, N)采样如下：c_k⁽ⁱ⁾=0的概率为π₀⁽ⁱ⁾；c_k⁽ⁱ⁾=1的概率为π₁⁽ⁱ⁾；c_k⁽ⁱ⁾=2的概率为π₂⁽ⁱ⁾；…；c_k⁽ⁱ⁾=T的概率为π_T⁽ⁱ⁾。

步骤6  粒子滤波目标状态修正。

基于数据关联结果c_k⁽ⁱ⁾对每个粒子进行卡尔曼滤波更新：

式中：m_k⁽ⁱ⁾和P_k⁽ⁱ⁾分别为粒子i在k时刻更新后的估计值和协方差；v_k⁽ⁱ⁾、S_k⁽ⁱ⁾和K_k⁽ⁱ⁾为过程变量。

重新计算粒子权重并估计有效粒子数：

如有效粒子数过低(如n_eff < N/4)，则进行重采样。

3 实验结果分析

本节分别针对仿真数据和雷达实测数据，评价本文算法的有效性，评价指标包括目标数N_d、虚警数N_fa、检测率P_d、虚警率P_fa、工作特征(Receiver Operator Characteristic, ROC)曲线、均方根误差(Root-Mean-Square Error, RMSE)等。

3.1 仿真数据

根据雷达目标散射原理，RCS随目标与雷达的相对角度变化，具有起伏特性。雷达回波起伏与RCS值相关，雷达与目标间角度的变化导致雷达回波强度变化，最终决定雷达图像中像素灰度值的大小。同时，低空雷达具有较宽的垂直波束，接收的回波中容易混入地杂波而导致虚警。本节的仿真数据以瑞利分布模拟地杂波，以4类χ²分布模拟低空目标的RCS概率密度分布^[17]。

模拟低空目标的4类χ²分布概率密度函数为

式中：为RCS平均值；σ为RCS随机变量；k = 1, 2, 3, 4代表 4类分布的双自由度数值。

模拟地杂波的瑞利分布概率密度函数为

式中：b为瑞利系数。

以离散维纳过程速度模型建立杂波环境中的单目标运动模型，其中目标的状态可以写为

式中：(x_k, y_k)为二维直角坐标系中的目标位置；为速度。

离散动态可以表示为线性、时不变构造方程：

式中：q_k－1为离散高斯过程白噪声，矩特征为

其中：时间步长设定为Δt=0.1；过程噪声的功率谱密度设定为q=0.1。

为构造杂波测量环境，如果测量值为噪声，则数据关联指数c_k设定为0，如果为真实目标则设定为1。杂波测量值均匀地分布在空间[－5, 5]×[－4, 4]中。真实目标的测量模型与附加的高斯噪声成线性关系。因此，可以将联合测量相似度表示为

量测噪声矩阵为

数据关联的先验值完全独立，写为

说明量测为杂波的概率是p_c，量测为目标的概率是1-p_c。p_c代表了杂波在量测中所占的比例，该值越大，跟踪的难度也越大。

如图 2所示，经过240步仿真，目标从(-4, -0.2)出发，在0~1.0 s以速度(1, 0)匀速运动，在1.1~3.0 s完成右转弯，在3.1~3.5 s以速度(0, -1)匀速运动，在3.6~5.5 s完成左转弯，在5.6~8.0 s以速度(1, 0)匀速运动，在8.1~10.5 s完成左转弯，在10.6~14.0 s以速度(0, 1)匀速运动，在14.1~16.5 s完成左转弯，在16.6~19.0 s以速度(-1, 0)匀速运动，在19.1~21.0 s完成左转弯，在21.1~21.5 s以速度(0, -1)匀速运动，在21.6~23.5 s完成右转弯，在23.6~24.0 s以速度(-1, 0)匀速运动，直至结束。杂波概率为0.3。针对4类目标RCS模型，按照图 2中的雷达目标运动轨迹，经过1000次蒙特卡罗仿真并取平均值，图 3给出不同杂波概率条件下(p_c=0.1, 0.3, 0.5)的ROC曲线和RMSE曲线，对比了采用本文算法对检测门限进行修正前后的跟踪结果，粒子数为10。

由图 3中4类模型的ROC曲线可见，本文算法均明显改善了4类目标模型的检测效果，其中，模型1的检测效果稍差，模型2~模型4的检测效果逐渐改善，这是由目标的散射特性决定的。具体而言，模型1的典型目标为前向观察的小型喷气飞机，模型2的典型目标为螺旋桨推进飞机、直升机等，模型3的典型目标为喷气飞机、大型民用客机，模型4的典型目标为侧向观测的导弹与高速飞行体等。显然，模型1的检测难度最大，模型2~模型4的检测难度逐渐降低。由于目标检测阈值在修正前后均采用了数据关联算法剔除杂波，p_c值对ROC曲线的影响不大。由图 3中4类模型的RMSE曲线可见，p_c值对RMSE曲线的作用明显，当p_c值提高时，RMSE明显提高，跟踪效果变差。同时发现，当杂波比例较低时(p_c=0.1, 0.3)，量测中的虚警较少，跟踪精度(RMSE)较高，此时，本文算法利用跟踪结果对检测门限进行反馈修正，能够进一步提高跟踪精度。但是，当杂波比例较高时(p_c=0.5)，杂波和虚警在量测中的比例达到50%，此时由于目标量测较少，导致跟踪精度(RMSE)降低。由实验结果可见，如果阈值设置过低(θ < 1.5)，检测门限经过修正后的跟踪结果反而可能进一步恶化，这可能是由于过低的阈值引入部分虚警导致的。因此，当虚警比例较高时，建议适当提高阈值(θ>1.5)，采用检测门限反馈修正算法能够在一定程度上改善跟踪结果，确保算法收敛。

针对4类目标RCS模型，表 1~表 4给出了粒子数N为10、50和100的情况下，本文算法给定不同分割阈值时的RMSE值，杂波比例均设定为p_c=0.1。通过对以上4组数据的分析，粒子数对检测结果影响不大，粒子数较多时的检测结果有时甚至稍逊于粒子数较少的时候。其原因在于本文算法采用Rao-Blackwellization方法将单目标跟踪与数据关联分开处理，将SMC(粒子滤波)方法用于数据关联，采用较少的粒子数，实现杂波与虚警量测中的多目标跟踪，可视为MHT的推广。Rao-Blackwellized粒子滤波的理论基础是某些滤波方程可以闭合的形式计算，其他采用蒙特卡罗采样，而不是对所有方程都采用采样方法。Rao-Blackwell的思想实现了较小的估计方差，其可视为用无穷集合去替代有限集合，往往能得到更为准确的结果。

3.2 实测数据

基于一组S波段非相参雷达数据和一组S波段相参雷达数据，在前期研究成果的基础上，采用本文算法进行处理，并与现有经典算法进行对比分析，2组数据的采样周期均为2.5s。第1组测试数据为S波段非相参雷达采集的图像序列，图中目标为一架波音737客机，分辨率1024×1024，量程22km，共120帧。图 4给出了某帧图像的检测跟踪结果。其中，图 4(a)在原始图像中标定了目标位置，分别采用固定阈值(Fixed Threshold, FT)、单元平均恒虚警(Cell Average-ConstantFalseAlarmRate, CA-CFAR)和本文算法进行目标检测；图 4(b)为目标跟踪结果。实验结果表明，目标检测阈值设置为θ=70时，FT算法在引入228个虚警的前提下能够检测到目标，CA-CFAR算法不仅能检测到目标，且将虚警数降为25个，基于前期研究成果^[18]，结合本文算法，粒子数为10，能够检测到小弱目标，且将虚警数减少到2个。图 4(b)为检测阈值设定为θ=120时，本文算法的目标跟踪结果，此时120帧图像的目标检测累积之和为N_d=119，N_fa=3，所有信息叠加在雷达背景图像上，3个虚警位置疑为未形成轨迹的低空未知目标。

为比较各类算法的鲁棒性，表 5给出了设置不同分割阈值的情况下，各类算法处理120帧雷达图像获取的目标数和虚警数总和，包括FT、CA-CFAR、最小选择恒虚警(Small of-Constant False Alarm Rate, SO-CFAR)和本文算法。形态学处理(Morphological Processing, MP)是一种经典的图像处理算法，对于剔除雷达图像中由单个像素组成的杂波信息，效果良好^[19]。此处，将MP算法作为所有算法的后处理，将“腐蚀”和“膨胀”处理相结合，既能剔除虚警，又能将属于同一目标的多个量测区域重新联通，避免将目标信息误认为虚警。由表 5可见，当检测到的目标数为114时，FT的虚警数为3539，CA-CFAR的虚警数为2469，而SO-CFAR的虚警数为3541。本文算法在设定不同阈值时，几乎能检测到所有目标，虚警数最小值仅为1。总体而言，本文算法最优，CA-CFAR算法优于FT算法，SO-CFAR算法甚至略逊于FT算法，其原因在于SO-CFAR算法通常将局部检测阈值设置过低，容易导致较高虚警率。

第2组测试数据在某雷达测试外场，由S波段相参雷达采集的轻小型无人机目标数据，该无人机目标沿正北偏西方向逐渐远离雷达飞行，量程6km。门限θ为目标RCS值，典型无人机目标的RCS值很小，属于小弱目标。

图 5中，采样周期为360个。表 6对比了设定不同θ值时的检测结果。本文算法将粒子数设定为10，能够在检测到339个无人机目标的情况下将虚警数降低到17。相参雷达已经剔除了静止的背景地物回波，且对虚警具有一定的抑制作用，本文算法的作用主要在于跟踪目标的同时提高对小弱目标的检测能力。

4 结论

本文结合SMC方法计算的粒子范围确定波门大小，在考虑粒子权重的前提下，利用检测单元与所有粒子的相对位置对检测门限进行修正，提高对小弱目标的检测能力。通过将本文算法与其他经典算法进行对比，得出以下结论：

1) 由仿真结果可见，本文算法对于各种类型雷达目标的检测结果均有明显改善，在杂波比例较低时，对跟踪精度也有一定提高。

2) 本文算法分别对每个目标进行卡尔曼滤波预估、更新和测量相似度估计，比纯粹的粒子滤波方法效率更高，以较少的粒子数就能实现对目标的精确跟踪。

3) 将本文算法与前期研究成果相结合，并应用于S波段非相参雷达采集的图像序列，检测结果表明，本文算法能够改善局部区域内的目标检测结果，优于其他经典算法。

4) 将本文算法应用于S波段相参雷达数据，能够提高对低空无人机目标的检测跟踪能力。