devbox
Monodyee
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

FMCW毫米波雷达中CFAR研究初探（附Python代码）_fmcw fft均值

作者：Monodyee | 2024-06-15 02:35:55

赞

踩

fmcw fft均值

多目标条件下的CFAR计算

汽车雷达的主要任务是探测前方区域内的所有目标，并计算目标的速度和位置信息。一般来讲，如果是在无噪声无杂波的背景下，目标检测会很容易，只需将雷达回波信号与一个信号固定门限比较，超过门限就会判定为目标。但在实际雷达探测应用中，由于地面，障碍物、雨云、箔条等干扰的存在，需要雷达在各种杂波中检测目标。恒虚警概率(CFAR)处理技术就是要在各种不同的杂波环境下，使雷达虚警概率保持在一个恒定范围内的信号处理方法。

Neyman-Pearson检测

虚警概率是指雷达前方没有目标而雷达检测到了目标，显然雷达检测到了假的目标，这个概率称为 P f a P_{fa} Pfa​ 。检测概率是指雷达前方有目标并且雷达也真实检测到了这个目标，这个概率称之为检测概率记为 P D P_{D} PD​ 。在实际检测中，我们总是希望 P f a P_{fa} Pfa​越小越好， P D P_{D} PD​越大越好，但事实上 P f a P_{fa} Pfa​和 P D P_{D} PD​是相互矛盾的。虚警概率越小，则对应检测门限越高，那么对目标信号漏检的概率则越大。如果提高检测门限，虽然可以减低虚警概率，但是检测概率也被降低了，造成了大量虚警现象。在实际雷达检测系统中，往往将 P f a P_{fa} Pfa​限制在某一特定水平，在给定虚警概率下，尽量提高检测概率，这就是Neyman-Pearson检测准则。

均值恒虚警（CA-CFAR）

一个FMCW雷达的典型信号处理流程如下图所示，ADC数据完成一次RangeFFT变换，然后在不同chirp的同一个RangeBin上在完成DopplerFFT，得到RDM（雷达数据矩阵）。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传在这里插入图片描述
(img-mfmV2ZU1-1655273377532)(Capture.JPG "雷达信号处理流程 ")]
恒虚警算法需要在RDM（雷达数据矩阵）上完成目标筛选，比如RDM中某一灰色单元我们称之为检测单元，周边单元称之为参考单元。一般可以认为参考单元和被检测单元的杂波模型都服从统一分布，如果我们RDM中的数据是代表的信号能量，也就是相当于雷达接收机使用平方律检波器。那么RDM上所有单元上信号符合指数分布。无目标时检测单元的概率密度函数（PDF）f(x)和累计分布函数（CDF）F(x)分别为：
f ( x ) = 1 u e − x u f(x)=\frac{1}{u}e^{\frac{-x}{u}} f(x)=u1​eu−x​ x>0
F ( x ) = 1 − e − x u F(x)=1-e^{\frac{-x}{u}} F(x)=1−eu−x​ x>0
首先我们计算参考单元的平均功率u，然后将平均功率uT作为门限，如果检测单元格功率高于此门限Tu，那么就可以认为检测单元格处有目标。当参考单元格数为N时，并且检测单元格处没有目标时，检测单元格虚警概率为： P f a N = ( 1 + T N ) − N P_{fa}^{N}=(1+T_{N})^{-N} PfaN​=(1+TN​)−N， 因此基于上述公式，我们计算 T = ( P f a N ) − 1 N − 1 T=(P_{fa}^{N})^{\frac{-1}{N}}-1 T=(PfaN​)N−1​−1
示例代码如下：

def ca1d(x, guard_len=4, noise_len=8, mode='wrap', pfa=10**(-6)):
    """在一维雷达数据上，首先利用ca-cfar在保持恒虚警概率pfa不变，计算各个检测
    单元格的门限值，高于这个门限，则定义为目标，低于这个门限则定义为噪声

    Args:
        x (~numpy.ndarray): 输入信号，一维ndarray.
        guard_len (int): 保护单元格长度，由于FFT计算点数限制，会导致检测单元格功率泄漏到临近单元格，为了准确估算噪底，邻近单元格设置为保护单元格.
        noise_len (int): 噪声单元格，这些单元格中的功率将被用于噪声估算.
        mode (str): 怎样处理开始和结束位置的单元格，本例中建议采用“wrap”模式，即x[0]的左侧单元格为x[-1]也就是序列的尾部单元格.
        pfa (float or int): 恒虚警概率，即雷达前方没有目标而雷达报告检测到目标，即虚假目标的概率.

    Returns:
        Tuple [ndarray, ndarray]
            1. (ndarray): ret反应x各个单元格是否为信号，即高于门限值
            #. (ndarray): threshold 各个单元格的门限值

    Examples:
        >>> signal = np.random.randint(10, size=20)
        >>> signal[10]=1000  ##远远高于噪声
            array([   9,    0,    2,    7,    9,    5,    7,    4,    3,    4, 1000,
                      4,    4,    1,    1,    5,    7,    2,    8,    5])
        >>> ret,threshold = ca1d(signal, guard_len=1, noise_len=3, mode='wrap', pfa=10**(-6))
        >>>ret
            array([False, False, False, False, False, False, False, False, False,
                   False,  True, False, False, False, False, False, False, False,
                   False, False])
        >>> threshold
            array([  45.,   54.,   63.,   36.,   36.,   27., 1530., 1539., 1539.,
                     36.,   27.,   27., 1521., 1530., 1530.,   36.,   36.,   27.,
                     36.,   27.])


    """
    ####噪声单元格数量
    N=2*noise_len
    ####依据公式计算,平方律前提下，固定恒虚警下门限倍数值T
    T=pfa**(-1/N)-1
    if isinstance(x, list):
        x = np.array(x)
    assert type(x) == np.ndarray
    ##利用卷积计算参考单元的功率和
    kernel = np.ones(1 + (2 * guard_len) + (2 * noise_len), dtype=x.dtype) / (2 * noise_len)
    ##将保护单元格置零
    kernel[noise_len:noise_len + (2 * guard_len) + 1] = 0
    ##计算参考单元的噪声功率和
    noise_floor = convolve1d(x, kernel, mode=mode)
    ##计算门限值
    threshold = noise_floor *T
    ret = (x > threshold)

    return ret,threshold


  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52

有序统计恒虚警（OS-CFAR)

CA-CFAR在多目标环境下表现糟糕，一旦有目标落入参考单元格，那么将极大的提高检测单元格的门限值，如果检测单元格内有目标也将被遮蔽。因此在多目标环境下，建议使用OS-CFAR。有序统计恒虚警检测器(ordered statistics—CFAR)是一种对参考滑窗内的参考单元进行排序处理的恒虚警处理算法，OS类检测方法在多目标环境中具有良好的分辨能力，这一点优于前面的均值类恒虚警处理器。同时OS类检测方法在均匀背景中的性能下降也是适度的。在有序类恒虚警检测器中，首先对参考单元采样值按照幅度大小作排序处理，然后，选取第k个排序采样值作为总的背景杂波功率水平估计值。k一般取为所有参考单元长度N的( 1 2 \frac{1}{2} 21​~ 3 4 \frac{3}{4} 43​),推荐值为 3 4 \frac{3}{4} 43​。 恒虚警概率 P f a P_{fa} Pfa​= k ( k N ) ( k − 1 ) ! ( T + N − k ) ! ( T + N ) ! k(_{k}^{N})\frac{(k-1)!(T+N-k)!}{(T+N)!} k(kN​)(T+N)!(k−1)!(T+N−k)!​ python代码如下：

import math
##参考单元格数量
N=16
##选取第k值作为噪声参考
k=round(N*0.75)

##门限乘积因子
T=21
##计算恒虚警
Pfa=k*math.comb(N,k)*math.factorial(k-1)*math.factorial(T+N-k)/(math.factorial(T+N))

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10

一般我们选取恒虚警概率 P f a P_{fa} Pfa​=1* 1 0 − 6 10^{-6} 10−6对应的N，k和T值如下表格所示，

P f a P_{fa} Pfa​NKT
1* 1 0 − 6 10^{-6} 10−68646
1* 1 0 − 6 10^{-6} 10−612927
1* 1 0 − 6 10^{-6} 10−6161221
1* 1 0 − 6 10^{-6} 10−6201518
1* 1 0 − 6 10^{-6} 10−6241816
1* 1 0 − 6 10^{-6} 10−6322415

OS-CFAR实现代码如下:

def os1d(x, guard_len=1, noise_len=4):
    """在一维雷达数据上，首先利用os-cfar在保持恒虚警概率pfa不变，计算各个检测
    单元格的门限值，高于这个门限，则定义为目标，低于这个门限则定义为噪声,在本函数中pfa固定为1e-6
            | $P_{fa}$ | N | K | T |
            |:---:|:---:|:---:|:---:|
            | 1*$10^{-6}$| 8 | 6 | 46 |
            | 1*$10^{-6}$| 12 | 9 | 27 |
            | 1*$10^{-6}$| 16 | 12 | 21 |
            | 1*$10^{-6}$| 20 | 15 | 18 |
            | 1*$10^{-6}$| 24 | 18 | 16 |
            | 1*$10^{-6}$| 32 | 24 | 15|

    Args:
        x (~numpy.ndarray): 输入信号，一维ndarray.
        guard_len (int): 保护单元格长度，由于FFT计算点数限制，会导致检测单元格功率泄漏到临近单元格，为了准确估算噪底，邻近单元格设置为保护单元格.
        noise_len (int): 噪声单元格，这些单元格中的功率将被用于噪声估算,如果noise_len小于4将被重置为4.
       

    Returns:
        Tuple [ndarray, ndarray]
            1. (ndarray): ret反应x各个单元格是否为信号，即高于门限值
            #. (ndarray): threshold 各个单元格的门限值

    Examples: >>>x=np.random.randint(4,size=32)
              >>>x[10]=100
              >>>x[14]=100
              >>>ret,thres=os1d(x,guard_len=1,noise_len=8)
              >>>ret
              >>>array([False, False, False, False, False, False, False, False, False,
                     False,  True, False, False, False,  True, False, False, False,
                     False, False, False, False, False, False, False, False, False,
                     False, False, False, False, False])



    """
    Tdict=dict([(8,46),(12,27),(16,21),(20,18),(24,16),(32,15)])
    if noise_len<4:
        noise_len=4
    elif noise_len>16:
        noise_len=16
    ####噪声单元格数量
    N=2*noise_len//4*4
    ####依据公式计算,平方律前提下，固定恒虚警下门限倍数值T
    k=round(N*3/4)
    T=Tdict[N]
    xlen=len(x)
    noise_floor=np.zeros_like(x)
    threshold=np.zeros_like(x)
    for i in range(xlen):
        LwinLeftIndex=np.mod(i-noise_len-guard_len,xlen)
        LwinRightIndex=np.mod(i-guard_len,xlen)
        if LwinLeftIndex>LwinRightIndex: #说明有部分噪声数据从队尾取，有部分从开头取
            Lwin=np.concatenate((x[LwinLeftIndex:],x[:LwinRightIndex]))
        else:
            Lwin=x[LwinLeftIndex:LwinRightIndex]
              
        RwinLeftIndex=np.mod(i+guard_len+1,xlen)
        RwinRightIndex=np.mod(i+guard_len+noise_len+1,xlen)
        if RwinLeftIndex>RwinRightIndex:  ##说明有部分数据需要从对头取
            Rwin=np.concatenate((x[RwinLeftIndex:], x[:RwinRightIndex]))
        else:
            Rwin=x[RwinLeftIndex:RwinRightIndex]
        

        window = np.concatenate((Lwin,Rwin))
        window.partition(k)
        noise_floor[i] = window[k]
        threshold[i] = noise_floor[i] * T
   
    ##计算门限值
    threshold = noise_floor *T
    ret = (x > threshold)

    return ret,threshold

  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70
  • 71
  • 72
  • 73
  • 74
  • 75
声明：本文内容由网友自发贡献，转载请注明出处：【wpsshop博客】
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号