一种脉冲雷达动目标检测方法

摘要

本发明涉及一种脉冲雷达动目标检测方法，包括：设置雷达系统参数以获取杂波图尺寸；根据所述雷达系统参数计算阈值；对所述杂波图进行二维DCT变换；根据所述阈值对经二维DCT变换的矩阵数据进行滤波；对经滤波后的矩阵数据进行IDCT变换，还原杂波图；及对所述还原后的杂波图检测。上述一种脉冲雷达动目标检测方法通过对杂波图进行二维DCT变换，DCT变换计算量较小，且具有明确的物理意义，精确度高；同时采用的阈值滤波，并不针对特定频率，具有更好的适应性。通过采用二维DCT滤波的脉冲雷达动目标检测方法可以简化计算同时具有更好的适应性。

权利要求

1.一种脉冲雷达动目标检测方法，其特征在于，包括下述步骤：设置雷达系统参数以获取杂波图尺寸；

根据所述雷达系统参数计算阈值；

对所述杂波图进行二维DCT变换；

根据所述阈值对经二维DCT变换的矩阵数据进行滤波；

对经滤波后的矩阵数据进行IDCT变换，还原杂波图；及对所述还原后的杂波图检测；

其中，对所述杂波图进行二维DCT变换是依据下述构造公式实现的：其中，式中，p表示变换域像素的纵坐标，q表示变换域像素的横坐标ap，aq是DCT变换的一组系数，m表示原图像像素的横坐标，n表示原图像像素的纵坐标，M表示原图像横向像素点的数量，N表示原图像纵向像素点的数量，Amn表示原图像中坐标为(m,n)的像素的值；

其中，对经滤波后的矩阵数据进行IDCT变换，还原杂波图是依据下述构造公式实现的：其中，式中，p表示变换域像素的纵坐标，q表示变换域像素的横坐标ap，aq是DCT变换的一组系数，m表示原图像像素的横坐标，n表示原图像像素的纵坐标，M表示原图像横向像素点的数量，N表示原图像纵向像素点的数量，Amn表示原图像中坐标为(m,n)的像素的值。

2.根据权利要求1所述的脉冲雷达动目标检测方法，其特征在于，其中，设置雷达系统参数以获取杂波图，包括下述步骤：根据雷达性能和检测范围选定帧大小；

根据待检测对象的运动特征设计采样系统；及根据所述采样系统需要的采样精度，获取杂波图尺寸。

3.根据权利要求1或2所述的脉冲雷达动目标检测方法，其特征在于，所述阈值包括：纵向低频阈值Tl、纵向高频阈值Th、横向低频阈值Tsl及横向高频阈值Tsh；

其中，所述纵向低频阈值Tl的计算方法为：Tl＝flow*M/(fV/2)

式中，flow为待检测对象运动频率，M为所述杂波图的行数，fV为慢时间采集频率；

其中，所述纵向高频阈值Th的计算方法为：Th＝fhi*M/(fV/2)

式中，fhi为待检测对象运动频率，M为所述杂波图的行数，fV为慢时间采集频率；

其中，所述横向低频阈值Tsl的计算方法为：Tsl＝fsl*N/(fH/2)

式中，fsl为实际脉冲波形，N为所述杂波图的列数，fH为快时间采集频率；

其中，所述横向高频阈值Tsh的计算方法为：Tsh＝fsh*N/(fH/2)

式中，fsh为实际脉冲波形，N为所述杂波图的列数，fH为快时间采集频率。

4.根据权利要求1所述的脉冲雷达动目标检测方法，其特征在于，其中，根据所述阈值对经二维DCT变换的矩阵数据进行滤波是依据下述构造公式实现的:式中，B(i，j)表示二维DCT变换后的变换域图像中坐标为(i,j)的像素的值。

5.根据权利要求3所述的脉冲雷达动目标检测方法，其特征在于，其中，对所述还原后的杂波图检测包括下述步骤：计算经还原后的杂波图上的每一列上的点的目标函数值；

记录超过幅度阈值的目标函数值对应的像素点位置；

统计上述像素点位置，并计算其概率密度；及选出所有密度大的点，即为目标位置。

6.根据权利要求5所述的脉冲雷达动目标检测方法，其特征在于，所述目标函数值为方差值，所述方差值的计算方法为：式中，xi(t)表示对应距离像上即还原后图像上的一纵列像素的值，mean(x)表示这一纵列像素所有样本的均值。

7.根据权利要求5所述的脉冲雷达动目标检测方法，其特征在于，所述目标函数值为每一列的峭度值，所述峭度值的计算方法为：式中，T为采样周期，μ取60*T，xi(t)表示对应距离像上即还原后图像上的一纵列像素的值， 表示峭度。

说明书

一种脉冲雷达动目标检测方法

技术领域

[0001] 本发明涉及脉冲雷达动目标检测技术，尤其涉及一种基于二维DCT(离散余弦变换)滤波的脉冲雷达动目标检测方法。

背景技术

[0002] 超宽带脉冲雷达已经广泛应用于社会生产的各个方面，其优越的定位性能和广适性能深受社会各界好评。然而雷达波在传输过程中，受到各种因素的影响；接收到的回波除携带有目标的运动信号外，还裹挟了非常多的干扰信号和噪声。这些干扰和噪声的存在严重影响了雷达的检测概率和恒虚警性能。

[0003] 通常算法仅对慢时间轴进行一维的滤波，没有用到快时间轴的信息量，效果不理想。采用传统的FIR(有限冲激响应)和FFT(快速傅里叶变换)的算法，运算量大，另外，通过加窗迭代计算二维滤波，本身就是一件非常消耗运算资源的事情，同时由于加窗截断数据，将导致时域信号被展宽，引入失真，尤其在需要保证一定探测范围的情况下，数据点数激增，运算量进一步加大。目前，虽然出现了基于DCT的雷达信号处理方法，但是均为一维DCT变换，这样仅利用了一维信号的信息量，精确度低且运算量大。

发明内容

[0004] 基于此，有必要针对上述脉冲雷达动目标检测技术存在的缺陷，提供一种精确度高同时运算量小的一种脉冲雷达动目标检测方法。

[0005] 一种脉冲雷达动目标检测方法，包括：设置雷达系统参数以获取杂波图尺寸；根据所述雷达系统参数计算阈值；对所述杂波图进行二维DCT变换；根据所述阈值对经二维DCT变换的矩阵数据进行滤波；对经滤波后的矩阵数据进行IDCT变换，还原杂波图；及对所述还原后的杂波图检测。

[0006] 在本实施例中，其中，设置雷达系统参数以获取杂波图，包括下述步骤：根据雷达性能和检测范围选定帧大小；根据待检测对象的运动特征设计采样系统；及根据所述采样系统需要的采样精度，获取杂波图尺寸。

[0007] 在本实施例中，所述阈值包括：纵向低频阈值Tl、纵向高频阈值Th、横向低频阈值Tsl及横向高频阈值Tsh；其中，所述纵向低频阈值Tl的计算方法为：Tl＝flow*M/(fV/2)[0008] 式中，flow为待检测对象运动频率，M为所述杂波图的行数，fV为慢时间采集频率；其中，所述纵向高频阈值Th的计算方法为：

[0009] Th＝fhi*M/(fV/2)

[0010] 式中，fhi为待检测对象运动频率，M为所述杂波图的行数，fV为慢时间采集频率；其中，所述横向低频阈值Tsl的计算方法为：

[0011] Tsl＝fsl*N/(fH/2)

[0012] 式中，fsl为实际脉冲波形，N为所述杂波图的列数，fH为快时间采集频率；

[0013] 其中，所述横向高频阈值Tsh的计算方法为：

[0014] Tsh＝fsh*N/(fH/2)

[0015] 式中，fsh为实际脉冲波形，N为所述杂波图的列数，fH为快时间采集频率。

[0016] 在本实施例中，其中，对所述杂波图进行二维DCT变换是依据下述构造公式实现的：

[0017]

[0018] 其中，式中，p表示像素的纵坐标，q表示像素的横坐标

[0019]

[0020] 0≤p≤M-1，0≤q≤N-1，ap，aq是DCT变换的一组系数，m表示原图像像素的横坐标，n表示原图像像素的纵坐标，M表示原图像横向像素点的数量，N表示原图像纵向像素点的数量，Amn表示原图像中坐标为(m，n)的像素的值。

[0021] 在本实施例中，其中，根据所述阈值对经二维DCT变换的矩阵数据进行滤波是依据下述构造公式实现的：

[0022]

[0023] 式中，B(i，j)表示二维DCT变换后的变换域图像中坐标为(i，j)的像素的值。

[0024] 在本实施例中，其中，对经滤波后的矩阵数据进行IDCT变换，还原杂波图是依据下述构造公式实现的：

[0025]

[0026] 其中，式中，p表示变换域像素的纵坐标，q表示变换域像素的横坐标[0027]

[0028] 0≤p≤M-1，0≤q≤N-1，ap·aq是DCT变换的一组系数。m表示原图像像素的横坐标，n表示原图像像素的纵坐标，M表示原图像横向像素点的数量，N表示原图像纵向像素点的数量，Amn表示原图像中坐标为(m，n)的像素的值。

[0029] 在本实施例中，其中，对所述还原后的杂波图检测包括下述步骤：计算经还原后的杂波图上的每一列上的点的目标函数值；记录超过幅度阈值的目标函数值对应的像素点位置；统计上述像素点位置，并计算其概率密度；及选出所有密度大的点，即为目标位置。

[0030] 在本实施例中，所述目标函数值为方差值，所述方差值的计算方法为：

[0031]

[0032] 式中，xi(t)表示对应距离像上即还原后图像上的一纵列像素的值，mean(x)表示这一纵列像素所有样本的均值。

[0033] 所述目标函数值为每一列的峭度值，所述峭度值的计算方法为：

[0034]

[0035] 其中：T为采样周期，μ通常取60*T，xi(t)表示对应距离像上即还原后图像上的一纵列像素的值， 表示峭度。

[0036] 上述一种脉冲雷达动目标检测方法通过对杂波图进行二维DCT变换，DCT变换计算量较小，且具有明确的物理意义，精确度高；同时采用的阈值滤波，并不针对特定频率，具有更好的适应性。通过采用二维DCT滤波的脉冲雷达动目标检测方法可以简化计算同时具有更好的适应性。

附图说明

[0037] 图1为本发明实施例提供的一种脉冲雷达动目标检测方法的流程示意图。

[0038] 图2为本发明实施例提供的设置雷达系统参数以获取杂波图的流程示意图。

[0039] 图3为本发明实施例提供的对还原后的杂波图检测的流程示意图。

[0040] 图4为本实施例提供的对每一帧进行去均值处理得到新的杂波图。

[0041] 图5为本发明提供的对杂波图进行二维DCT变换得到的杂波图B。

[0042] 图6为本发明实施例提供的对经滤波后的矩阵数据进行IDCT变换得到的杂波图C。

[0043] 图7为本发明实施例提供的对4个窗口进行联合分析的图样。

[0044] 图8为本发明实施例提供的最终的判决图。

具体实施方式

[0045] 请参阅图1，为本发明实施例提供的一种脉冲雷达动目标检测方法的流程示意图，包括下述步骤：

[0046] 步骤S10：设置雷达系统参数以获取杂波图尺寸。

[0047] 请参阅图2，为本发明实施例提供的设置雷达系统参数以获取杂波图的流程示意图，步骤S10具体为：

[0048] 步骤S11：根据雷达性能和检测范围选定帧大小。

[0049] 雷达性能指雷达本身对信号的采样速率和对信息的存储能力。假设雷达的最小距离分辨率为R，存储深度为M，则一帧数据代表的距离为R*M；若R和M可调，则根据需要调整R和M的值，使其满足具体需求。

[0050] 步骤S12：根据待检测对象的运动特征设计采样系统。

[0051] 采样系统应根据实际应用进行设计。对于指定的应用，若已知待检测对象的运动特征引起的慢时间最大频率偏移为f，则选取采样频率为f的3～5倍即可，而无需追求精确采样雷达中频信号；而对于快时间采样频率由雷达本身的性能决定，与雷达硬件能够做出反应的最快时间相对应。

[0052] 步骤S13：根据采样系统需要的采样精度，获取杂波图尺寸。

[0053] 根据步骤S12得到慢时间采样频率，而快时间采样精度由实际应用决定，综合考虑距离分辨率R，快时间轴上表示的帧数num-M，慢时间轴上的频率分辨率Rm等参数。假如慢时间上的帧率为nr，则为了达到慢时间轴上的频率分辨率Rm，需要的成帧数N＝nr/Rm；则杂波图尺寸为：M×num-M×N，其中M为快时间一帧的点数，N为图像纵向尺寸。

[0054] 步骤S20：根据雷达系统参数计算阈值。

[0055] 雷达系统参数可包括：快时间采样频率为fH，慢时间采样频率为fV及截取的杂波图尺寸，记为M行N列；阈值包括：纵向低频阈值Tl、纵向高频阈值Th、横向低频阈值Tsl及横向高频阈值Tsh；这里取待检测对象运动频率为0.2Hz～5Hz并分别记为flow～fhi，实际脉冲波形fsl～fsh，则：

[0056] Tl＝flow*M/(fV/2)

[0057] Th＝fhi*M/(fV/2)

[0058] Tsl＝fsl*N/(fH/2)

[0059] Tsh＝fsh*N/(fH/2)。

[0060] 步骤S30：对杂波图进行二维DCT变换。

[0061] DCT是一种实数域变换，其变换核为实数余弦函数；对一幅图像进行离散余弦变换后，许多有关图像的重要信息都集中在DCT变换的一小部分系数中。在本发明提供的实施例中，对杂波图进行二维DCT变换是依据下述构造公式实现的：

[0062]

[0063] 其中，式中，p表示变换域像素的纵坐标，q表示变换域像素的横坐标[0064]

[0065] 0≤p≤M-1，0≤q≤N-1，

[0066] ap，aq是DCT变换的一组系数，m表示原图像像素的横坐标，n表示原图像像素的纵坐标，M表示原图像横向像素点的数量，N表示原图像纵向像素点的数量，Amn表示原图像中坐标为(m，n)的像素的值。

[0067] 步骤S40：根据阈值对经二维DCT变换的矩阵数据进行滤波。

[0068] DCT是一种实数域变换，其变换核为实数余弦函数。对一幅图像进行离散余弦变换后，许多有关图像的重要可视信息都集中在DCT变换的一小部分系数中。根据DCT变换的意义可知，整帧图像的能量向左上角聚集，越靠近左上角表示频率越低。因此根据雷达系统的参数可以对将要设计的滤波器进行计算。记快时间采样频率为fH，慢时间采样频率为fV，截取的杂波图尺寸为M行N列；这里取待检测对象的运动频率0.2Hz～5Hz记为flow～fhi，实际脉冲波形fsl～fsh，则：

[0069] Tl＝flow*M/(fV/2)

[0070] Th＝fhi*M/(fV/2)

[0071] Tsl＝fsl*N/(fH/2)

[0072] Tsh-fsh*N/(fH/2)

[0073]

[0074] 式中，Tl为纵向低频阈值，Th为纵向高频阈值，Tsl为横向低频阈值，Tsh为横向高频阈值，B(i，j)表示二维DCT变换后的变换域图像中坐标为(i，j)的像素的值。根据上式对杂波图的DCT系数构成的新图进行滤波，指定面积内的DCT系数将被清零，而保留了我们需要的频段内的数据。当然如果要进行更加复杂的滤波只需依据上式计算出要滤除的频段对应的面积大小和位置，然后对相应的DCT系数进行清零或者特定的修改即可形成新的算法，该方法适应性很广。

[0075] 步骤S50：对经滤波后的矩阵数据进行IDCT变换，还原杂波图。

[0076] 在本发明提供的实施例中，对经滤波后的矩阵数据进行IDCT变换，还原杂波图是依据下述构造公式实现的：

[0077]

[0078] 其中，式中，p表示变换域像素的纵坐标，q表示变换域像素的横坐标[0079]

[0080] 0≤p≤M-1，0≤q≤N-1；ap，aq是DCT变换的一组系数。m表示原图像像素的横坐标，n表示原图像像素的纵坐标，M表示原图像横向像素点的数量，N表示原图像纵向像素点的数量，Amn表示原图像中坐标为(m，n)的像素的值。

[0081] 步骤S50：对还原后的杂波图检测。

[0082] 请参阅图3，为本发明实施例提供的对还原后的杂波图检测的流程示意图，包括下述步骤：

[0083] 步骤S51：计算经还原后的杂波图上的每一列上的点的目标函数值。

[0084] 对于经过DCT滤波并反变换的杂波图，其尺寸仍为N×M，N行M列。对每一列上的N个点进行统计计算目标函数的值，这里以协方差为例进行说明：设第i列上的像素点为xi(t)t∈[1，N]，则该距离像上的方差

[0085] 式中，xi(t)表示对应距离像上即还原后图像上的一纵列像素的值，mean(x)表示这一纵列像素所有样本的均值。

[0086] 可以理解，除协方差外，还可以各列的峭度作为目标函数进行计算。峭度代表的是一组数据分布的集中程度，当数据样本很大，或者未知时，应该满足正态分布，其峭度为3；若因某些原因导致数据的分布发生改变，则峭度必然偏离3，因此我们也可以通过计算其峭度判断是否有目标。峭度的计算式可采用简化的迭代算法：

[0087]

[0088] 其中：T为采样周期；μ通常取60*T，xi(t)表示对应距离像上即还原后图像上的一纵列像素的值， 表示峭度。这是一个迭代式，通过该式遍历样本即可求得数值。。

[0089] 步骤S52：记录超过幅度阈值的目标函数值对应的像素点位置。

[0090] 对于采用方差作为目标函数的情况选取最大值的80％左右阈值，即计算所有M个距离像上的方差值，取其最大值的80％为阈值，方差大于阈值的距离像点i，i∈[1，M]被记录到一个新的数组Index里。

[0091] 对于采用峭度作为目标函数的情况，将峭度值与3做减法，然后取其绝对值，选出所有距离像上大于最大值80％的距离像，表示有目标出现，而小于80％的数据进行丢弃。

[0092] 步骤S53：统计上述像素点位置，并计算其概率密度。

[0093] 在本发明提供的实施例中，为保证鲁棒性，多次重复上述像素点位置的方法；重复次数以4～10次之间为宜。统计Index里每个i出现的次数，计算其概率(出现次数除以重复次数)。

[0094] 步骤S54：选出所有密度大的点，即为目标位置。

[0095] 选出密度大的N个点，即为目标位置：对Index中的概率密度取峰值，有多少个峰就有多少个运动目标，且其对应的距离像i，即为目标的距离位置。

[0096] 可以理解，密度的大小取决于人为的选择，在本实施例中，优选为大于最大密度值的80％为大的密度。

[0097] 以下通过具体实施例对本发明的进行进一步说明。

[0098] 为描述方便，假设快时间上每帧存储深度为512个点，点与点之间的对应的距离约为8毫米；慢时间轴采用20Hz的采样频率。我们选取200帧为一个判决窗口，则每隔10秒钟可以得到一个200行，512列的杂波图，记为Frame_org，其中每个像素点记为f(i，n)，其中i表示快时间轴上的位置，i的取值范围是[1，512]；n表示慢时间轴上的位置，n的取值范围是[1，200]。

[0099] 请参阅图4，为本实施例提供的对每一帧进行去均值处理得到新的杂波图，记为Frame。该图是雷达取得的原始图像，横轴表示快时间轴，这里横轴上共有512个点，纵坐标上我们采用256点为一个窗口，共用了4个窗口进行拼接，一共1024点。

[0100] 请参阅图5，为本发明提供的对杂波图进行二维DCT变换得到的杂波图B。离散余弦变换(DCT)是一种实数域变换，其变换核为实数余弦函数；对一幅图像进行离散余弦变换后，许多有关图像的重要信息都集中在DCT变换的一小部分系数中。一个矩阵的二维DCT定义如下：

[0101]

[0102] 其中，式中，p表示像素的纵坐标，q表示像素的横坐标

[0103]

[0104] 0≤p≤M-1，0≤q≤N-1，ap，aq是DCT变换的一组系数，m表示原图像像素的横坐标，n表示原图像像素的纵坐标，M表示原图像横向像素点的数量，N表示原图像纵向像素点的数量，Amn表示原图像中坐标为(m，n)的像素的值。。

[0105] 杂波图B中，横坐标和纵坐标是原始图像经DCT变换后的系数坐标，表示为余弦函数的频率。从杂波图B看出，采用上述构造公式(1)对Frame进行二维DCT，得到与原杂波图具有同样大小的DCT系数，且图像重要的信息都集中在图像的左上部。

[0106] 在得到上述图样后，根据阈值对经二维DCT变换的矩阵数据进行滤波。具体如下：记快时间采样频率为fH，慢时间采样频率为fV，截取的杂波图尺寸为M行N列；这里取人体运动频率0.2Hz～5Hz记为flow～fhi，实际脉冲波形fsl～fsh，则：

[0107] Tl＝flow*M/(fV/2)

[0108] Th＝fhi*M/(fV/2)

[0109] Tsl＝fsi*N/(fH/2)

[0110] Tsh＝fsh*N/(fH/2)

[0111]

[0112] 式中，T1为纵向低频阈值，Th为纵向高频阈值，Tsl为横向低频阈值，Tsh为横向高频阈值。根据上式对杂波图的DCT系数构成的新图B进行滤波，指定面积内的DCT系数将被清零，而保留了我们需要的频段内的数据。

[0113] 请参阅图6，为本发明实施例提供的对经滤波后的矩阵数据进行IDCT(DCT反变换)变换得到的杂波图C。

[0114] 在本实施例中，对经滤波后的矩阵数据进行IDCT变换，还原杂波图是依据下述构造公式实现的：

[0115]

[0116] 其中，式中，p表示变换域像素的纵坐标，q表示变换域像素的横坐标[0117]

[0118] 0≤p≤M-1，0≤q≤N-1；ap·aq是DCT变换的一组系数。m表示原图像像素的横坐标，n表示原图像像素的纵坐标，M表示原图像横向像素点的数量，N表示原图像纵向像素点的数量，Amn表示原图像中坐标为(m，n)的像素的值。

[0119] 杂波图C中，表示滤波还原后的时域图，从杂波图C可以看出在几个距离像上有运动的痕迹。依据得到的杂波图C，对每一列数据进行统计，这里选择以能量为目标函数进行判断，计算其协方差(或者交流能量)，大于最大值80％的峰值可称为明显较大，将满足这个条件的点定位目标位置。除协方差外，还可以各列的峭度作为目标函数进行计算。峭度代表的是一组数据分布的集中程度，当数据样本很大，或者未知时，应该满足正态分布，其峭度为3；若因某些原因导致数据的分布发生改变，则峭度必然偏离3，因此我们也可以通过计算其峭度判断是否有目标。峭度的计算式可采用简化的迭代算法：

[0120]

[0121] 其中：T为采样周期，μ通常取60*T，xi(t)表示对应距离像上即还原后图像上的一纵列像素的值， 表示峭度。

[0122] 请参阅图7，为本发明实施例提供的对4个窗口进行联合分析的图样。在图7表示的是对各个距离像求方差后的结果，具体为先对每个窗口进行计算，然后将4个窗口的值拼接起来得到最后的判据。为保险起见，可对多个这样的窗口进行多次计算，然后对结果进行蒙特卡洛分析，选取出现概率最大的几个，作为最终结果。从图7中表明，本发明实施例提供的方法可以有效提取运动及其距离信息。

[0123] 请参阅图8，为本发明实施例提供的最终的判决图。图8中，横坐标表示距离像，纵坐标表示对应的距离像上，大于最大方差值80％的点的个数。从图5中清晰地看出有3个运动目标，其位置分别为0.3米，1米，3.4米，与实验情况吻合。

[0124] 以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。