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

多目标条件下的CFAR计算

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

Neyman-Pearson检测

虚警概率是指雷达前方没有目标而雷达检测到了目标，显然雷达检测到了假的目标，这个概率称为$P_{fa}$ 。检测概率是指雷达前方有目标并且雷达也真实检测到了这个目标，这个概率称之为检测概率记为$P_D$ 。在实际检测中，我们总是希望$P_{fa}$越小越好，$P_D$越大越好，但事实上$P_{fa}$和$P_D$是相互矛盾的。虚警概率越小，则对应检测门限越高，那么对目标信号漏检的概率则越大。如果提高检测门限，虽然可以减低虚警概率，但是检测概率也被降低了，造成了大量虚警现象。在实际雷达检测系统中，往往将$P_{fa}$限制在某一特定水平，在给定虚警概率下，尽量提高检测概率，这就是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)=\frac{1}{u}{e}^{\frac{-x}{u}}$ x>0
$F(x)=1-e^{\frac{-x}{u}}$ x>0
首先我们计算参考单元的平均功率u，然后将平均功率uT作为门限，如果检测单元格功率高于此门限Tu，那么就可以认为检测单元格处有目标。当参考单元格数为N时，并且检测单元格处没有目标时，检测单元格虚警概率为： $P_{fa}^N=(1+T_N{)}^{-N}$， 因此基于上述公式，我们计算$T=(P_{fa}^N{)}^{\frac{-1}{N}}-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

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

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

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

$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

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