CN110244289B - 一种自适应粒子滤波地波雷达目标一体化探测方法 - Google Patents

Abstract

一种自适应粒子滤波地波雷达目标一体化探测方法，包括利用高频地波雷达获得的多个通道数据来构建距离‑多普勒‑时间三维数据，对三维数据进行杂波抑制、低门限CFAR预处理、疑似目标点初选，初始化目标粒子，计算状态转移后的粒子权重，根据终止采样条件进行粒子自适应选取，然后进行目标状态估计，最终得到目标的航迹集合。相对于固定粒子数粒子滤波方法，本发明针对地波雷达在复杂环境下不同信噪比的目标进行粒子个数的自适应选取。对于高信噪比目标，可以有效减少跟踪所需的粒子个数，提高算法效率；对于复杂情况下的目标，通过自适应增加粒子个数来提高目标跟踪性能。

Description

一种自适应粒子滤波地波雷达目标一体化探测方法

技术领域

本发明涉及一种地波雷达目标检测与跟踪一体化方法，具体涉及一种自适应粒子滤波地波雷达目标一体化探测方法。

背景技术

高频地波雷达(High Frequency Surface Wave Radar，HFSWR)是大范围海上船只目标监视监测的主要手段，利用高频电磁波(3～30MHz)沿海面爬行来实现超视距的目标(船只、低空飞机等)探测，可以实时提供目标的位置、航速航向等信息，探测距离最远可达300公里。

一般传统的高频地波雷达目标探测方法采用的是先检测后跟踪(Detect-before-Track，DBT)的思想，该方法对于回波能量较弱、信噪比较低的目标检测效果不佳，导致目标跟踪性能下降。检测与跟踪一体化方法可以解决弱目标检测困难的问题，而检测前跟踪方法是实现目标检测与跟踪一体化的重要方法，其思想是：对单帧雷达回波数据不进行目标有无判断，而是利用目标在时空上的关联特性和杂波噪声的随机性，先进行多帧数据积累，从而实现同一目标的回波能量的积累，由此提高目标信噪比，完成目标的检测和判决。检测前跟踪方法由于不设置检测门限，充分利用目标的原始信息，减少了先检测后跟踪过程中的航迹关联问题，降低了算法复杂度。

目前，国内外发展了各种TBD(Track-before-Detect，TBD)算法用来实现地波雷达目标检测与跟踪一体化探测。主要实现方法包含：三维匹配滤波(3-D matched filters)、投影变换、动态规划(Dynamic Programming，DP)、粒子滤波(particle filter，PF)等。粒子滤波相对于动态规划(Grossi E et al，2013)等方法，具有估计的目标状态在理论上最优，且适合于类似地波超视距雷达这种非线性、非高斯的系统，并且算法实现容易且因其递归结构，有存储量小等优点。

一般传统的粒子滤波算法确定粒子个数基本采用经验值，或者通过大量的蒙特卡洛实验，从而选取合适的粒子数目，粒子的数目基本上在粒子初始化阶段就固定下来。Rutten首次建立了背景噪声为瑞利噪声时的雷达TBD模型，并提出了基于优效重要性采样的PF-TBD算法，该算法提高了粒子使用效率，降低了粒子使用的个数。SU Hongtao等人采用粒子滤波检测前跟踪算法，针对复数量测、功率量测、目标幅度量测等不同量测构造检测似然比，实现了对天波超视距雷达的弱目标检测和跟踪。这些算法均采用固定数目的粒子去实现目标的检测与跟踪，而且大多基于仿真数据，缺少地波雷达实测环境应用。一般地波雷达背景噪声复杂、目标信噪比强弱差别较大。对于背景噪声简单，目标信噪比高的情况，只需要少量的粒子就可以满足系统的跟踪需求；对于复杂的目标环境，目标信噪比较低，为了保证跟踪精度，通常需要的粒子数目较大，而粒子滤波算法的复杂度与采样的粒子个数成正比。因此，如果统一采用固定的大量粒子这样就会导致计算量骤增，使得计算时间变长，影响地波雷达系统目标跟踪的实时性；如果采用较少的粒子，又无法满足跟踪的精度要求。所以，需要发展自适应的粒子选取方法解决这一问题。相关部分的参考文献如下：

[1]GrossiE，Lops M，Venturino L.A Novel Dynamic Programming Algorithmfor Track-Before-Detect in Radar Systems[J].IEEE Transactions on SignalProcessing，2013，61(10)：2608-2619.

[2]Mark G.Rutten，Neil J.Gordon，Maskell S.Particle-based track-before-detect in Rayleigh noise[C]//Signal&Data Processing of Small Targets.2004.

[3]Su H T，SHUIPL，Liu H W，et al.Particle Filter Based Track-before-detect Algorithm for Over-the-horizon Radar Target Detection and Tracking[J].Chinese Journal of Electronics.Vol.18，No.1，Jan.2009.

发明内容

本发明的目的是提供一种自适应粒子滤波地波雷达目标一体化探测方法，该方法能够充分利用粒子权重包含的信息，以自适应的方式确定粒子个数，实现高频地波雷达目标检测与跟踪一体化。

为了实现上述目的，本发明的技术方案如下：

一种自适应粒子滤波地波雷达目标一体化探测方法，包括以下步骤：

步骤0：利用高频地波雷达获得的多个通道数据来构建距离-多普勒-时间(Range-Doppler-Time，R-D-T)三维数据；

步骤1：对每一帧R-D谱数据依次进行杂波抑制、低门限CFAR预处理，得到预处理后的数据点；

其特征是还包括以下步骤：

步骤2：疑似目标点初选，具体如下：

(a)记录每一帧的预处理后的数据点在R-D谱上的位置并按幅值大小进行降序排序；

(b)对第一帧幅值最大的点，找出第一帧中与该幅值最大点的距离和多普勒均小于门限Δ的点，并删除；对于其他帧，以此类推；未被删除的点全部作为疑似目标点；所述门限Δ的距离维取值范围不超过3，多普勒维取值范围是不超过3；

步骤3：利用步骤2得到第一帧疑似目标的距离和多普勒位置，对第一个疑似目标进行随机粒子采样，作为目标的初始化粒子；

以疑似目标距离减1个单元格作为下限，加1个单元格作为上限，记为[R_min，R_max]；以疑似目标多普勒减2个单元格作为下限，加2个单元格作为上限，记为[D_min，D_max]；

在采样范围[R_min，R_max]、[D_min，D_max]内按照均匀分布函数

采样N个粒子，N为自然数，且100≤N≤5000，

表示第一帧中第1个疑似目标的状态，包括距离、速度、加速度信息；

表示第一帧中第1个疑似目标的均匀分布函数；

表示采样后的粒子，1≤j≤N；

步骤4：从第二帧起，挑选出第一帧第1个疑似目标的权重最大的粒子x_1，max，根据此粒子位置信息进行新粒子的随机采样；即按照均匀分布函数

进行新粒子随机采样，采样范围同上，采样数为1；

其中，x_1，max表示第一帧第一个目标的权重最大的粒子，U(x_1，max)表示粒子x_1，max的均匀分布函数，

表示采样后的粒子；

将采样得到的新粒子

按照如下公式进行状态转移，作为在第二帧的新粒子

计算公式如下；

其中，F为目标运动状态转移矩阵，为已知量；G为噪声协方差矩阵，为已知量；v_k-1为第k-1帧噪声向量，包含目标加速度变化引起的噪声和目标幅值波动引起的噪声，为可设置变量；该步骤中k＝2；

步骤5：完成状态转移后，结合现有的高频地波雷达目标观测模型计算第二帧的新粒子

的权重w_k，具体计算公式如下：

其中，z_k(r，d)为第k帧疑似目标在(r，d)位置的量测值，A_k(r，d)为疑似目标在(r，d)位置的回波幅度值，I₀(·)表示零阶修正的贝塞尔函数；δ是背景噪声的标准差，为已知量；C_r(x_k)＝{r-q，...r-1，r，r+1，...，r+q}，C_d(x_k)＝{d-p，...d-1，d，d+1，...，d+p}；r和d为已知量，分别表示x_k在测量数据中对应的距离单元和多普勒单元的坐标；p和q为小于等于3的自然数，分别表示受疑似目标点影响的距离单元和多普勒单元数目，该步骤中k＝2；

步骤6：对于第二帧，重复步骤4和5的通过均匀分布函数进行的粒子采样操作，将每一次操作所得的新粒子个数进行累加，将累加得粒子个数记为n；同时将每一次操作所得的新粒子的权重进行累加，将累加得到的粒子的权重之和记为S；

n＝n+1

S＝S+w_k

当粒子个数n小于设置的最小粒子数N_min(为不小于100的自然数)或者当粒子个数n大于N_min且小于设置的最大粒子数N_max(大于100且不小于2000的自然数)，但粒子权重和S小于设定的阈值η(10⁶～10⁸)，均需要跳到步骤(4)继续采样粒子，其他情况下，停止采样；

步骤7：当步骤6停止采样后，将累加得到的粒子权重和S与设定的阈值η进行比较，判断目标是否存在，若粒子权重和S大于阈值η则说明目标存在，否则目标不存在；

步骤8：采样得到的所有粒子进行归一化处理，然后以现有的系统重采样方法进行粒子的重采样；

步骤9：根据步骤7判断目标存在后，根据停止采样后所得的粒子状态和权重对该存在的目标之状态进行估计，计算公式如下；

其中，x_k表示第k帧目标的状态，

表示第k帧第i个粒子的的权重，

表示第k帧第i个粒子的状态，该步骤中k＝2；

步骤10：对第三帧之后的所有帧，均开始执行步骤4-7的操作；

若经过步骤4-7的操作之后，连续三帧目标不存在，则终止步骤4-7的操作；否则，对

第三帧之后的所有帧，均执行步骤4-7的操作；

步骤11：跳到步骤(3)，对其他疑似目标完成步骤3-10的检测与跟踪操作；

步骤12：当所有疑似目标均完成检测与跟踪操作后，将得到的目标无向航迹集合作为最终输出。

与现有技术相比，本发明的创新之处体现在以下方面：

本发明根据粒子权重包含的信息，从重要性密度函数采样粒子来近似目标的真实状态。通过设置固定的粒子权重和，然后从随机产生的新粒子里自适应选取粒子，直到满足终止采样条件。终止采样条件包含固定粒子权重和、最小和最大粒子数。通过设置最小粒子数，防止因目标信噪比过高使得自适应选取的粒子过少，导致粒子多样性丧失，影响跟踪；通过设置最大粒子数，防止目标消失或信噪比过低使得自适应选取的粒子权重太小，浪费计算资源。

本发明提出的自适应粒子滤波方法，可以针对地波雷达这种复杂环境下不同信噪比目标进行自适应选取粒子。对于高信噪比目标，有效减少跟踪所需的粒子个数；对于复杂情况下的目标，通过自适应增加粒子个数来提高跟踪精度。

本发明针对固定粒子滤波计算量大、无法针对复杂情况自适应调整粒子个数的问题，通过提出自适应粒子滤波方法，实现了应用于地波雷达的自适应粒子滤波目标检测与跟踪一体化。

附图说明

图1为本发明的基本流程示意图。

图2为高频地波雷达R-D-T三维数据。

图3为粒子初始化示意图

图4为粒子状态转移与目标位置估计示意图

图5为使用本发明处理后的跟踪结果。

具体实施方式

下面结合附图，对本发明的方法做进一步说明：

一种自适应粒子滤波地波雷达目标一体化探测方法，如图1所示，包括以下步骤：

步骤0：利用高频地波雷达获得的多个通道数据来构建距离-多普勒-时间(Range-Doppler-Time，R-D-T)三维数据，作为自适应粒子滤波算法的输入数据；如图2所示；

步骤1：对每一帧R-D谱数据依次进行杂波抑制、低门限CFAR预处理，得到预处理后的数据点；

其特征是还包括以下步骤：

步骤2：疑似目标点初选，具体如下：

(a)记录每一帧的预处理后的数据点在R-D谱上的位置并按幅值大小进行降序排序；

(b)对第一帧幅值最大的点，找出第一帧中与该幅值最大点的距离和多普勒均小于门限Δ的点，并删除；对于其他帧，以此类推；未被删除的点全部作为疑似目标点；所述门限Δ的距离维取值范围不超过3，多普勒维取值范围是不超过3；

步骤3：利用步骤2得到第一帧疑似目标的距离和多普勒位置，对第一个疑似目标进行随机粒子采样，作为目标的初始化粒子；如图3所示；

以疑似目标距离减1个单元格作为下限，加1个单元格作为上限，记为[_Rmin，R_max]；以疑似目标多普勒减2个单元格作为下限，加2个单元格作为上限，记为[D_min，D_max]；

在采样范围[R_min，R_max]、[D_min，D_max]内按照均匀分布函数

采样N个粒子，N为自然数，且100≤N≤5000，

表示第一帧中第1个疑似目标的状态，包括距离、速度、加速度信息；

表示第一帧中第1个疑似目标的均匀分布函数；

表示采样后的粒子，1≤j≤N；

步骤4：从第二帧起，挑选出第一帧第1个疑似目标的权重最大的粒子x_1，max，根据此粒子位置信息进行新粒子的随机采样；即按照均匀分布函数

进行新粒子随机采样，采样范围同上，采样数为1；

其中，x_1，max表示第一帧第一个目标的权重最大的粒子，U(x_1，max)表示粒子x_1，max的均匀分布函数，

表示采样后的粒子；

将采样得到的新粒子

按照如下公式进行状态转移，作为在第二帧的新粒子

计算公式如下；

其中，F为目标运动状态转移矩阵，为已知量；G为噪声协方差矩阵，为已知量；v_k-1为第k-1帧噪声向量，包含目标加速度变化引起的噪声和目标幅值波动引起的噪声，为可设置变量；该步骤中k＝2；如图4所示；

步骤5：完成状态转移后，结合现有的高频地波雷达目标观测模型计算第二帧的新粒了

的权重w_k，具体计算公式如下：

其中，z_k(r，d)为第k帧疑似目标在(r，d)位置的量测值，A_k(r，d)为疑似目标在(r，d)位置的回波幅度值，I₀(·)表示零阶修正的贝塞尔函数；δ是背景噪声的标准差，为已知量；C_r(x_k)＝{r-q，...r-1，r，r+1，...，r+q}，C_d(x_k)＝{d-p，...d-1，d，d+1，...，d+p}；r和d为已知量，分别表示x_k在测量数据中对应的距离单元和多普勒单元的坐标；p和q为小于等于3的自然数，分别表示受疑似目标点影响的距离单元和多普勒单元数目，该步骤中k＝2；

步骤6：对于第二帧，重复步骤4和5的通过均匀分布函数进行的粒子采样操作，将每一次操作所得的新粒子个数进行累加，将累加得粒子个数记为n；同时将每一次操作所得的新粒子的权重进行累加，将累加得到的粒子的权重之和记为S；

n＝n+1

S＝S+w_k

当粒子个数n小于设置的最小粒子数N_min(为不小于100的自然数)或者当粒子个数n大于N_min且小于设置的最大粒子数N_max(大于100且不小于2000的自然数)，但粒子权重和S小于设定的阈值η(10⁶～10⁸)，均需要跳到步骤(4)继续采样粒子，其他情况下，停止采样；

步骤7：当步骤6停止采样后，将累加得到的粒子权重和S与设定的阈值η进行比较，判断目标是否存在，若粒子权重和S大于阈值η则说明目标存在，否则目标不存在；

步骤8：采样得到的所有粒子进行归一化处理，然后以现有的系统重采样方法进行粒子的重采样；

步骤9：根据步骤7判断目标存在后，根据停止采样后所得的粒子状态和权重对该存在的目标之状态进行估计，计算公式如下；

其中，x_k表示第k帧目标的状态，

表示第k帧第i个粒子的的权重，

表示第k帧第i个粒子的状态，该步骤中k＝2；

步骤10：对第三帧之后的所有帧，均开始执行步骤4-7的操作；

若经过步骤4-7的操作之后，连续三帧目标不存在，则终止步骤4-7的操作；否则，对

第三帧之后的所有帧，均执行步骤4-7的操作；

步骤11：跳到步骤(3)，对其他疑似目标完成步骤3-10的检测与跟踪操作；

步骤12：当所有疑似目标均完成检测与跟踪操作后，将得到的目标无向航迹集合作为最终输出，输出结果如图5所示。

Claims (1)

1.一种自适应粒子滤波地波雷达目标一体化探测方法，包括以下步骤：

步骤0：利用高频地波雷达获得的多个通道数据来构建距离-多普勒-时间三维数据，即R-D-T数据；

步骤1：对每一帧R-D谱数据依次进行杂波抑制、低门限CFAR预处理，得到预处理后的数据点；

其特征是还包括以下步骤：

步骤2：疑似目标点初选，具体如下：

(a)记录每一帧的预处理后的数据点在R-D谱上的位置并按幅值大小进行降序排序；

(b)对第一帧幅值最大的点，找出第一帧中与该幅值最大点的距离和多普勒均小于门限Δ的点，并删除；对于其他帧，以此类推；未被删除的点全部作为疑似目标点；所述门限Δ的距离维取值范围不超过3，多普勒维取值范围是不超过3；

步骤3：利用步骤2得到第一帧疑似目标的距离和多普勒位置，对第一个疑似目标进行粒子随机采样：

以疑似目标距离减1个单元格作为下限，加1个单元格作为上限，记为[R_min,R_max]；以疑似目标多普勒减2个单元格作为下限，加2个单元格作为上限，记为[D_min,D_max]；

在采样范围[R_min,R_max]、[D_min,D_max]内按照均匀分布函数

采样N个粒子，N为自然数，且100≤N≤5000，

表示第一帧中第1个疑似目标的状态，包括距离、速度、加速度信息；

表示第一帧中第1个疑似目标的均匀分布函数；

表示采样后的粒子，1≤j≤N；

步骤4：挑选出上一帧第一个疑似目标的权重最大的粒子x_1,max，根据此粒子位置信息进行新的粒子随机采样：

对于第二帧，挑选出上一帧第一个疑似目标的权重最大的粒子x_1,max，根据此粒子位置信息进行粒子随机采样，即按照均匀分布函数

进行粒子随机采样，采样范围同上，采样数为1；

其中，U(x_1,max)表示粒子x_1,max的均匀分布函数，

表示采样后的粒子；

将采样得到的粒子

按照如下公式进行状态转移，作为下一帧的粒子

计算公式如下；

其中，F为目标运动状态转移矩阵，为已知量；G为噪声协方差矩阵，为已知量；v_k-1为第k-1帧噪声向量，包含目标加速度变化引起的噪声和目标幅值波动引起的噪声，为可设置变量；k为帧数，k≥2；

步骤5：完成状态转移后，结合现有的高频地波雷达目标观测模型计算下一帧的粒子的权重w_k，具体计算公式如下：

其中，z_k(r,d)为第k帧疑似目标在(r,d)位置的量测值，A_k(r,d)为疑似目标在(r,d)位置的回波幅度值，I₀(·)表示零阶修正的贝塞尔函数；δ是背景噪声的标准差，为已知量；

C_r(x_k)＝{r-q,…r-1,r,r+1,…,r+d},C_d(x_k)＝{d-p,…d-1,d,d+1,…,d+p}；r和d为已知量，分别表示x_k在测量数据中对应的距离单元和多普勒单元的坐标；p和q为小于等于3的自然数，分别表示受疑似目标点影响的距离单元和多普勒单元数目，k为帧数，k≥2；

步骤6：对于后续帧，按照步骤4和5的通过均匀分布函数进行的粒子随机采样操作，将每一次操作所得的粒子个数进行累加，将累加得粒子个数记为n；同时将每一次操作所得的粒子的权重进行累加，将累加得到的粒子的权重之和记为S；

n＝n+1

S＝S+w_k

k为帧数，k≥2，

当粒子个数n小于设置的最小粒子数N_min，N_min为不小于100的自然数，或者当粒子个数n大于N_min且小于设置的最大粒子数N_max，N_max为大于100且不小于2000的自然数，但粒子权重和S小于设定的阈值η，10⁶≤η≤10⁸，均需要跳到步骤4继续采样粒子，其他情况下，停止采样；

步骤7：当步骤6停止采样后，将累加得到的粒子权重和S与设定的阈值η进行比较，判断目标是否存在，若粒子权重和S大于阈值η则说明目标存在，否则目标不存在；

步骤8：采样得到的所有粒子进行归一化处理，然后以现有的系统重采样方法进行粒子的重采样；

步骤9：根据步骤7判断目标存在后，根据停止采样后所得的粒子状态和权重对该存在的目标之状态进行估计，计算公式如下；

其中，x_k表示第k帧目标的状态，

表示第k帧第i个粒子的的权重，

表示第k帧第i个粒子，k为帧数，k≥2；

步骤10：对第二帧及之后的所有帧，均以步骤4-7的方式进行操作；

若以步骤4-7的方式进行操作之后，连续三帧目标不存在，则终止操作；

步骤11：对所有疑似目标按照步骤3-10的方式完成检测与跟踪操作；

步骤12：当所有疑似目标均完成检测与跟踪操作后，将得到的目标无向航迹集合作为最终输出。

Images