自适应雷达检测：第2部分

有没有想过雷达是如何生成点云的？现在不用再疑惑了，因为在本文关于雷达检测的内容中，我们将学习如何把雷达数据处理成点云！这将是一个实战教程，让我们开始吧！

这是关于自适应雷达系列文章中的第二篇。

• 第1部分：CFAR介绍
• 第2部分：密集雷达点云（本篇）
• 第3部分：稀疏雷达点云

照片来自 USGS，来自 Unsplash

我们将使用来自Radical数据集的数据，该数据集可以从这里下载。除非另有说明，所有照片均为作者拍摄。让我们检查数据并在单个帧上生成点云。主要代码片段如下，完整代码请参见笔记本。

    从 h5py 导入  

    # 设置数据路径  

    # 获取数据  
    data_dict = {}  
    with h5py.File(DATA_PATH, 'r') as h5_obj:  
        for key in h5_obj.keys():  
            data_dict.update({key : np.asarray(h5_obj[key])})  

    data_idx = 29  
    将 adc_data 设置为 data_dict['radar'][data_idx, ...]

现在我们将处理数据，具体的方法在这里不作详细介绍。如果你想知道更多细节，可以查看这个系列：系列。等我们处理完数据后，就可以开始检测了。

    # 执行快速时间FFT
    range_cube = np.fft.fft(adc_data, axis=2).transpose(2, 1, 0)

    # 执行慢时间多普勒变换
    range_doppler = np.fft.fftshift(np.fft.fft(range_cube, axis=2), axes=2)

    # 执行方位角处理
    steering_vector = compute_steering_vector(num_ant=8, angle_res=1, angle_rng=90)
    range_azimuth = capon_spectrum(range_doppler, steering_vector)

    # 计算功率谱密度
    range_azimuth_psd = np.abs(range_azimuth)**2

到目前为止，我们已经得到了范围和方位角的数据，让我们看看它画出来是什么样的。

图1. 范围-方位谱图。

高功率区域显示为更亮的颜色斑块，我们可以看到有几个显而易见的目标，看看设置一个静态阈值来检测这些目标会怎样。使用平均值似乎能取得不错的效果。

    rng_az_detections = range_azimuth_psd > range_azimuth_psd.mean() # `rng_az_detections` 表示 `range_azimuth_psd` 大于其平均值的那些值。

图2，左图是距离-方位谱；右图是粗略检测到的距离-方位目标。

这种方法确实能产生一些结果，但我们可以进一步清理数据，并且可能使用更稳健的阈值方法。不过我们会做得更好，我们将使用单元平均恒虚警率（CFAR，即恒虚警率），这种方法可以在未知特性的噪声环境中实现自适应检测，更多细节请参阅第一部分。现在我们将展示如何在距离-方位数据上实现这种方法。

我们将首先编写一个用于一维信号的CA-CFAR函数，然后将其应用于距离-方位矩阵的每一行和每一列中。

    从 scipy.ndimage 导入 convolve1d  

    def ca_cfar(x, prob_fa=0.05, num_train=8, num_guard=2):  
        """   
            用于根据指定虚警概率获得动态阈值的单元平均虚警函数。使用训练单元的双侧平均值来估计噪声水平。  
            输入：  
                x - 输入信号向量  
                prob_fa - 虚警概率（0-1）  
                num_train - 用于估计噪声的训练单元数  
                num_guard - 保护单元数  
            输出：  
                noise_level - 返回噪声水平和阈值  
                threshold - 动态阈值  
        """  
        # 计算动态阈值因子  
        a = num_train * (prob_fa ** (-1 / num_train) - 1)  

        # 获取核以高效计算信号向量上的 CFAR，其中核的中间部分是0，以排除保护单元的影响  
        cfar_kernel = np.ones((1 + 2 * num_guard + 2 * num_train), dtype=x.dtype) / (2 * num_train)  
        cfar_kernel[num_train : num_train + (2 * num_guard) + 1] = 0.  

        # 计算动态阈值  
        noise_level = convolve1d(x, cfar_kernel, mode='nearest')  
        # 其中 convolve1d 函数使用最近邻模式进行卷积  
        threshold = a * noise_level  

        return noise_level, threshold

我们可以这样使用CFAR：

    # 范围方向的CFAR
    noise_level_rng, threshold_rng = np.apply_along_axis(ca_cfar,   
                                                     axis=0,   
                                                     arr=range_azimuth_psd,  
                                                     prob_fa=0.05,   
                                                     num_train=10,   
                                                     num_guard=3)  

    # 方位角方向的CFAR
    noise_level_az, threshold_az = np.apply_along_axis(ca_cfar,   
                                                     axis=0,   
                                                     arr=range_azimuth_psd.T,  
                                                     prob_fa=0.05,   
                                                     num_train=10,   
                                                     num_guard=3)  

    # 获取CFAR峰值点
    range_peaks = range_azimuth_psd > threshold_rng  
    azimuth_peaks = range_azimuth_psd > threshold_az.T  

    峰值 = 范围峰值点 & 方位峰值点

看看这会怎样

图3展示了范围和方位数据的初步CFAR检测结果。

这比使用静态阈值要好得多，我们可以实际地从这张图像中分离出各个峰值。我们可以通过一些非极大值抑制来清理这些内容，这是计算机视觉和信号处理领域中常用的一种技术。代码如下所示，所有 _cfar_nms() 函数所做的就是在 Range-Azimuth 矩阵中的每个像素上进行迭代，并将其与在半径为 2 的范围内所有其他像素进行比较。如果它是该半径内的最大像素，则设置检测为 1，否则设置为 0（无检测）。

    def get_nhood(arr, x, y, axis=None, nsize=2):  

        if not axis:  
            hood = arr[y-nsize:y+nsize+1, x-nsize:x+nsize+1]  
        elif axis == 0:  
            hood = arr[y-nsize:y+nsize+1, x]  
        elif axis == 1:  
            hood = arr[y, x-nsize:x+nsize+1]  

        return hood  

    def cfar_nms(cfar_peaks, signal_arr, nsize=2):  
        """ 对检测到的CFAR峰值执行非极大值抑制操作 """  
        nms_peaks = np.zeros(cfar_peaks.shape)  
        for y, x in zip(*np.where(cfar_peaks == 1)):  
            if np.all(signal_arr[y, x] >= get_nhood(signal_arr, x, y, None, nsize)):  
                nms_peaks[y, x] = 1  

        return nms_peaks  

    # 执行 NMS
    nms_peaks = cfar_nms(peaks, range_azimuth_psd, nsize=2)

图4. 非完全抑制CA-CFAR在距离-方位数据上的检测

这好多了，现在我们有了实际的点云数据，但我们还可以做一件事。

既然我们用的是雷达，为什么不利用我们测量到的多普勒速度信息呢？我们首先可以利用。

    ## 获取多普勒CFAR  
    range_doppler_mag = np.abs(range_doppler.sum(axis=1))**2  

    # 获取适当的缩放范围多普勒功率谱  
    range_doppler_psd = np.abs(range_doppler.sum(axis=1)/8)**2   
    range_doppler_psd = 10*np.log10(range_doppler_psd)  

    noise_level_dop, threshold_dop = np.apply_along_axis(ca_cfar,   
                                                 axis=0,   
                                                 arr=range_doppler_mag.T,  
                                                 prob_fa=0.0001,   
                                                 num_train=10,   
                                                 num_guard=2)  

    doppler_peaks = range_doppler_mag > threshold_dop.T

图5. 显示的多普勒速度测量（峰值部分）。

多普勒速度检测中的长线对应于零多普勒速度的分Bin。不在该线上的任何目标都表示检测到了运动。现在我们将根据我们已经检测到的距离目标，获取这些多普勒峰值的相关位置信息。

    doppler_locs = doppler_peaks[int(num_adc_samples) - range_locs, :]

接着我们将得到所有检测的距离、角度以及多普勒频移。

    # 获取范围、角度和多普勒数据  
    ranges = (num_adc_samples - range_locs) * range_resolution  
    angles = angle_vector[azimuth_locs]  
    dopplers = range_doppler_psd[num_adc_samples - range_locs, :]

最后我们可以把我们检测到的点云组合起来。

    # 获取检测到的范围和角度
    detected_ranges = ranges[np.where(doppler_locs)[0]]
    detected_angles = angles[np.where(doppler_locs)[0]]
    detected_dopplers = np.where(doppler_locs)[1]  # 多普勒值
    detected_doppler_mag = dopplers[np.where(doppler_locs)[0], np.where(doppler_locs)[1]]

    detections = np.stack((detected_ranges,
                           detected_angles,
                           detected_dopplers,
                           detected_doppler_mag)).T  # 将检测结果堆叠成一个数组

当我们显示数据时，我们可以根据检测到的多普勒速度来确定哪些检测对象是移动的。我们将点云叠加在范围-方位图上，并用洋红色标记移动点。

图6. 点云检测结果叠加在范围方位谱上。

这里的情况很难判断，所以我们来做个GIF，并把GIF和摄像头的画面放一起，这样我们就能清楚看到发生了什么。

图7. CFAR 目标检测 GIF。左：摄像头画面。右：检测结果叠加在距离-方位谱上。

在这一 GIF 中，我们可以看到哪些检测结果对应于移动的行人，并且可以看到速度计大多数时候是正确的。在这个 GIF 中，我们并没有显示他们的速度，只显示他们是否在移动。另一个需要注意的地方是，即使没有实际物体移动，检测结果也会波动（上下跳跃）。这是由于雷达信号中的噪声造成的，这是正常现象，因此我们的算法需要足够稳定来应对这种情况。最后要观察的是当行人通过后面的玻璃门进入时。随着行人的前进，他们的运动通过大量的检测结果来指示，这是因为玻璃门中的金属会为雷达产生强烈的反射，而打开门则会产生雷达信号的多普勒频移，从而让我们能够测量他们的相对速度。

这篇帖子到这里就结束了，然后在下次我们会看到如何高效生成稀疏的雷达点云数据，这更接近市面上汽车雷达实际生成的效果。

• 第1部分：CFAR介绍
• 第2部分：密集雷达点云（本文）
• 第3部分：稀疏雷达点云