登堂入室:毫米波雷达开发手册之信号模型

写在前面

本文档用于为实现基于AWR1243BOOST等单板毫米波雷达开发提供参考指南与解决方案，主要包括硬件配置、基础参数、信号模型、应用DEMO开发以及可深入研究方向思考等；为更好地匹配后续级联雷达应用的学习路线，在本手册中会尽可能同化单板雷达和级联雷达中的相关表述。

Xuliang，联系方式：xxxxxxxxx@xxx.edu.cn。未经本人允许，请勿用于商业和学术用途。

本章节为信号模型章节，主要解读时域信号模型、空域信号模型和信号处理栈。

信号模型分析

时域信号模型

通过比较雷达的TX和RX信号，可以估计出感知任务的基本参数：目标数目\(K_{tgt}\和每个目标\(k_{tgt}\的目标距离\(r_{k_{tgt}}\，目标径向速度\(v_{r_{k_{tgt}}}\，方位角\(\phi_{k}\和俯仰角\(\epsilon_k\。在每个测量周期中检测到的目标集合被称为目标列表，并且可以被解释为从雷达测量获得的点云。后续的毫米波雷达信号处理步骤主要有是聚类（即识别属于同一物理目标）、跟踪（即从目标列表的序列中导出物理上存在的对象的信息和统计数据）。本章节的核心是介绍TDMA-MIMO框架下的信号模型仿真及实现方法。

chirp组成的chirp序列构成的发射波形（绿色部分）和来自单个目标的接受波形信号（蓝色部分）。

chirp-chirp间隔的为\(T_{cc}\的\(Q\个线性调频脉冲波形\(s_c（t）\的序列给出：

\[s(t=\sum\limits_{q=1}^Q s_{\mathrm{c}}(t-q T_{\mathrm{cc}} \]

chirp的信号表示为：

\[s_\mathrm{c}(t=e^{j(2\pi f_1t+\pi\mu t^2+\varphi_0}\mathrm{rect}\left(\frac{t-T_\mathrm{c}/2}{T_\mathrm{c}}\right. \]

chirp斜率为\(\mu\，起始频率为\(f_1\，截止频率为\(f_h=f_l+W\，\(W=\mu T_c\表示chirp的带宽，chirp的中心频率\(f_c\定义为\(f_1+W/2\，一帧chirp序列的长度或称帧周期为\(T_{seq}=QT_{cc}-(T_{cc}-T_c\。

\[r(t=\sum_{k_{t g t}=1}^{K_{t g t}}a_{k_{t g t}}s(t-T_{\mathrm{d},k}(t \]

其中，目标衰减系数\(a_{k_{tgt}}\是由于自由空间传播和目标RCS引起的系数加权，往返时间\(T_{d,k}(t\是时间相关的，因为这考虑了目标的移动，通常帧周期为几毫秒或更小，与目标速度相比较短，无法形成积累效应。因此，\(T_{\mathrm{d},k}(t\可以近似表示为：

\[T_{\mathrm{d},k}(t=\frac{2r_{k_\mathrm{tgt}}+2v_{r,k_\mathrm{tgt}}t}{c_0}. \]

FMCW雷达的基本原理是，发射信号\(s（t）\和接收信号\(r（t）\相乘混频，混频后通过低通滤波器LPF以滤除在和频（约两倍载波频率）处振荡的和分量，保留低通滤波信号即中频信号IF Signal，然后与发射信号定时同步地ADC采样。对于每个chirp，ADC采样将生成P个样本的向量，这意味着\(Q\个线性调频脉冲的采样将生成\(Q\个向量，并且每个向量长度为\(P\，以二维矩阵\({\mathcal{B}}=(b_{p,q},{\mathcal{B}}\in{{\mathbb{C}}^{P\times Q}}\的形式存储，每个元素为\(b_{p,q}\。\(p\所在的维度称为快时间维，\(q\所在的维度称为慢时间维度，这在上一节内容中已经提及。

\[b[p,q]\approx\sum_{k_1\mathfrak{g}=1}^{k_1}a_{k_1\mathfrak{g}n}e^{j(\lambda_{k_1\mathfrak{g}n}p+\zeta_{k_1\mathfrak{g}n}q}e^{j\varphi_{k_1\mathfrak{g}n}}+\omega[p,q]. \]

其中，对应距离维度的归一化角频率可表示如下：

\[\lambda_{k_\mathrm{tgt}}=2\pi\left(\frac{2\mu r_{k_\mathrm{tgt}}}{c_0}+f_{\mathrm{d},k}\rightT_\mathrm{s}, \]

\[\xi_{k_{\mathrm{tgt}}}=2\pi\left(\frac{2f_{\mathrm{c}}v_{r,k_{\mathrm{tgt}}}}{c_{0}}\rightT_{\mathrm{cc}} \]

目标的初始相位则可表示如下：

\[\varphi_{k_\mathrm{tgt}}=\frac{2f_c r_{k_\mathrm{tgt}}}{c_0}. \]

空域信号模型

SIMO）配置；为了进一步提高角分辨率，在共同定位的相干MIMO配置中的若干发射和接收天线可用于现代毫米波雷达中，形成\(N=N_{TX}N_{RX}\的虚拟接收阵列。

MIMO后的虚拟阵列效果，左边天线标记为接收天线，右边天线标记为发射天线，可以很明显地看到TX和TX1/TX2在俯仰维度上存在偏移量，RX2和其他三根接收天线在俯仰维度上存在偏移量。为更好地理解MIMO机制，下面将手把手推导如何得到虚拟阵列图(假设每个方格代表阵元间距d。

TX1在四根接收天线上形成的水平相位为[0 w 2w 3w]，俯仰相位为[0 w 0 0]

TX2在四根接收天线上形成的水平相位为[w 2w 3w 4w]，俯仰相位为[0 w 0 0]
TX3在四根接收天线上形成的水平相位为[2w 3w 4w 5w]，俯仰相位为[2w 3w 2w 2w]

TI AWR2243Cascade具有12个发射天线和16个接收天线，实际在方位维度仅有86个虚拟阵元。为了排列准确，笔者写了以下程序来实现物理阵元到虚拟阵列的精确映射。

function cfgOut = ConfigureRTX(
	% 下面这个程序实现了不同级联雷达的物理发射阵元和接收阵元配置，并考虑了天线发射顺序（这对虚拟阵元的排列影响不大，但会影响adc数据中的数据排列）。
	% 实际的阵元间距可能与下面的存在误差，以下板子参数多为本人估算所得。
	% 具体代码结合注释和上下文语句进行分析，这里不再赘述。
	% By Xuliang
	
    cfgOut.Mode = 4; % 选择需要仿真的雷达模式 这里提供了6种毫米波雷达配置，分别为TI
    
    if cfgOut.Mode == 1 
        disp(strcat(['=====','已选择TI AWR1243-单板雷达模式','=====']; % 单板雷达模式3T4R
        cfgOut.numTx = 3;
        cfgOut.numRx = 4;  
        
        cfgOut.PosTX_X = [0 2 4]; % 发射天线水平方向阵元间距
        cfgOut.PosTX_Y = [0 1 0]; % 发射天线垂直方向阵元间距
        cfgOut.PosTX_BOARD_ID = [1 2 3]; % 发射天线板上顺序 后面没用到
        cfgOut.PosTX_Trans_ID = [1 3 2]; % 发射天线发射顺序 
        
        cfgOut.PosRX_X = [0 1 2 3]; % 接收天线水平方向阵元间距
        cfgOut.PosRX_Y = [0 0 0 0]; % 接收天线垂直方向阵元间距
        cfgOut.PosRX_BOARD_ID = [1 2 3 4]; % 接收天线板上顺序
        
    elseif cfgOut.Mode == 2
        disp(strcat(['=====','已选择TI AWR2243-级联雷达模式','=====']; % 级联雷达模式12T16R
        cfgOut.numTx = 12;
        cfgOut.numRx = 16;
        cfgOut.PosTX_X = [11 10 9 32 28 24 20 16 12 8 4 0];
        cfgOut.PosTX_Y = [6 4 1 0 0 0 0 0 0 0 0 0]; 
        cfgOut.PosTX_BOARD_ID = [12 11 10 3 2 1 9 8 7 6 5 4]; 
        cfgOut.PosTX_Trans_ID = [12:-1:1]; 
        
        cfgOut.PosRX_X = [11 12 13 14 50 51 52 53 46 47 48 49 0 1 2 3];
        cfgOut.PosRX_Y = zeros(1,cfgOut.numRx; 
        cfgOut.PosRX_BOARD_ID = [13 14 15 16 1 2 3 4 9 10 11 12 5 6 7 8]; 
        
    elseif cfgOut.Mode == 3  % 这个模式为特斯拉6T8R级联板
    	disp(strcat(['=====','已选择特斯拉6T8R-级联雷达模式','====='];
        cfgOut.numTx = 6;
        cfgOut.numRx = 8;
        cfgOut.PosTX_X = [0 17 34 41 48 55];
        cfgOut.PosTX_Y = [0 0 0 0 4 8 12]; 
        cfgOut.PosTX_BOARD_ID = [1 2 3 4 5 6]; 
        cfgOut.PosTX_Trans_ID = [1 2 3 4 5 6]; 
        
        cfgOut.PosRX_X = [0 2 3 5 8 11 14 17];
        cfgOut.PosRX_Y = zeros(1,cfgOut.numRx; 
        cfgOut.PosRX_BOARD_ID = [1 2 3 4 5 6 7 8]; 
        
    elseif cfgOut.Mode == 4 % 这个模式为加特兰 4T4R单板
    	disp(strcat(['=====','已选择加特兰4T4R-单板雷达模式','====='];
        cfgOut.numTx = 4;
        cfgOut.numRx = 4;
        cfgOut.PosTX_X = [0 5 12 18];
        cfgOut.PosTX_Y = [5 0 0 0]; 
        cfgOut.PosTX_BOARD_ID = [1 2 3 4]; 
        cfgOut.PosTX_Trans_ID = [4 3 2 1]; 
        
        cfgOut.PosRX_X = [0 2 5 7];
        cfgOut.PosRX_Y = [0 0 0 2]; 
        cfgOut.PosRX_BOARD_ID = [1 2 3 4]; 
        
    elseif cfgOut.Mode == 5 % 这个模式为福瑞泰克 6T8R级联板
    	disp(strcat(['=====','已选择福瑞泰克6T8R-级联雷达模式','====='];
        cfgOut.numTx = 6;
        cfgOut.numRx = 8;
        cfgOut.PosTX_X = [0 2 4 6 8 10];
        cfgOut.PosTX_Y = [0 0 0 0 0 0]; 
        cfgOut.PosTX_BOARD_ID = [1 2 3 4 5 6]; 
        cfgOut.PosTX_Trans_ID = [6 5 4 3 2 1]; 
        
        cfgOut.PosRX_X = [0 2 5 13 13+ones(1,4*1] + [0 0 0 0 0 2 5 13]; % 暂时不确定双级联间距
        cfgOut.PosRX_Y = [1 0 0 0 1 0 0 0 0]; 
        cfgOut.PosRX_BOARD_ID = [1 2 3 4 5 6 7 8]; 
        
    elseif cfgOut.Mode == 6 % 这个模式为TI双级联毫米波雷达
    	disp(strcat(['=====','已选择TI 6T8R-级联雷达模式','====='];
        cfgOut.numTx = 6;
        cfgOut.numRx = 8;
        cfgOut.PosTX_X = [0 2 4 7 11 15];
        cfgOut.PosTX_Y = [0 1.6 0 0 0 0]; 
        cfgOut.PosTX_BOARD_ID = [1 3 2 4 5 6]; 
        cfgOut.PosTX_Trans_ID = [6 5 4 3 2 1]; 
        
        cfgOut.PosRX_X = [0 1 2 3 19 20 21 22];
        cfgOut.PosRX_Y = [0 0 0 0 0 0 0 0 0]; 
        cfgOut.PosRX_BOARD_ID = [1 2 3 4 5 6 7 8]; 
        
    end
    
    virtual_azi_arr = repmat(cfgOut.PosTX_X(cfgOut.PosTX_Trans_ID, cfgOut.numRx, 1 + repmat(cfgOut.PosRX_X(cfgOut.PosRX_BOARD_ID, cfgOut.numTx, 1.';
    virtual_ele_arr = repmat(cfgOut.PosTX_Y(cfgOut.PosTX_Trans_ID, cfgOut.numRx, 1 + repmat(cfgOut.PosRX_Y(cfgOut.PosRX_BOARD_ID, cfgOut.numTx, 1.';
    
    virtual_array.azi_arr = reshape(virtual_azi_arr,[],1;
    virtual_array.ele_arr = reshape(virtual_ele_arr,[],1;
    virtual_array.tx_id = reshape(repmat(cfgOut.PosTX_Trans_ID, cfgOut.numRx, 1,[],1;
    virtual_array.rx_id = reshape(repmat(cfgOut.PosRX_BOARD_ID, cfgOut.numTx, 1.',[],1;
    virtual_array.virtual_arr = cat(2,virtual_array.azi_arr,virtual_array.ele_arr,virtual_array.rx_id,virtual_array.tx_id;
    
    cfgOut.virtual_array = virtual_array; % 将virtual_array参数存入更新
    
end

% 测试文件：可视化虚拟阵元排列
% By Xuliang

clc;clear;close all

cfgOut = ConfigureRTX(;
plot(cfgOut.virtual_array.virtual_arr(:,1, cfgOut.virtual_array.virtual_arr(:,2,'ro';
xlabel('方位维';ylabel('俯仰维'；grid minor；
% 这里不再进行可视化。

假设目标满足远场窄带信号假设，那么可进一步扩展时域信号至空域信号如下：$$b[p,q,n]\approx\sum_{k\text{g}s=1}^{{K_\text{g}t}a_{k\text{g}t}e}t}p+\zeta_{k\text{g}s}q}e^{{i\varphi_{k\text{g}t}}e}(\mathfrak{u}_{k\text{g}s}\mathfrak{p}_n}+\omega[p,q,n].$$

3D tensor或radar cube数据的元素，根据ISO8000:2标准定义的球面坐标系，目标\(k_{tgt}\对应的单位长度波达角矢量可以表示如下：

\[\mathbf{u}_k=\mathbf{r}_{k_\textrm{tgt}}/\|\mathbf{r}_{k_\textrm{tgt}}\|=\begin{bmatrix}\cos\phi_{k_\textrm{tgt}}\sin\theta_{k_\textrm{t}}\\ \sin\phi_{k_\textrm{gt}}\sin\theta_{k_\textrm{tgt}}\\ \cos\theta_{k_\textrm{tgt}}\end{bmatrix} \]

\[\mathbf{p}_n=\mathbf{p}_{\mathrm{RX},n_{\mathrm{RX}}}+\mathbf{p}_{\mathrm{TX},n_{\mathrm{TX}}},n_{\mathrm{RX}}=n\bmod N_{\mathrm{RX}},\quad n_{\mathrm{TX}}=n\bmod N_{\mathrm{TX}} \]

至此，我们可以很明确地建立目标角度和模型\(e^{(2\pi/\lambda\mathbf{u}_{k_\mathrm{tgt}}\mathbf{p}_n}\有关，上面提到的目标初始相位\(\varphi_{k_\mathrm{tgt}}=\frac{2f_c r_{k_\mathrm{tgt}}}{c_0}\可以衡量目标在参考坐标系原点处出现的相位，不大影响角度的估计。

TDMAI-MIMO下的radar cube数据文件，核心代码如下所示：

function [signal] = GenerateSigIQ(tarOut, cfgOut
    %% 本文件用于毫米波雷达的ADC-IQ信号
    %% By Xuliang
    tarNum = tarOut.nums; % 目标数目
    numTx = cfgOut.numTx; % 发射天线数目
    numRx = cfgOut.numRx; % 接收天线数目
    
    ADCNum = cfgOut.ADCNum; % ADC采样数目
    Frame = cfgOut.Frame; % 帧数
    ChirpNum = cfgOut.ChirpNum; % 每帧发射Chirp数目
    TotalChirpNum = Frame * ChirpNum; % 总Chirp数目
    
    % 信号参数
    fc = cfgOut.fc; % 载频 Hz
    fs = cfgOut.fs; % ADC采样频率 Hz 
    Ramptime = cfgOut.Ramptime; % 工作时间
    Idletime = cfgOut.Idletime; % 空闲时间
    Slope = cfgOut.Slope; % Chirp斜率
    validB = cfgOut.validB; % 实际有效带宽
    totalB = cfgOut.totalB; % 完整带宽
    Tc = cfgOut.Tc; % 单Chirp周期
    ValidTc = cfgOut.ValidTc; % 单Chirp内ADC有效采样时间
    
    % 硬件参数
    Pt = cfgOut.Pt; % dbm 发射功率
    Fn = cfgOut.Fn; % 噪声系数
    Ls = cfgOut.Ls;  % 系统损耗
    
    % 天线阵列
    antennaPhase = cfgOut.antennaPhase; % 天线相位
    virtual_array = cfgOut.virtual_array; % 虚拟阵列 这里的虚拟阵列是struct virtual_arr中包含了发射方位相位和俯仰相位 接收天线 发射天线顺序 
    arr = virtual_array.virtual_arr;
    
    % 初阶版本顺序发射可行
%     arr = cfgOut.array; % 虚拟阵元排列
    arrNum = size(arr, 1; % 虚拟阵列数目
    arrdelay = zeros(1, arrNum + Tc * reshape(repmat([0:numTx-1], numRx,1, 1, arrNum; % [0 0 0 0 1 1 1 1 2 2 2 2]
%     arrdelay = zeros(1, arrNum; % 这一项相当于常数项 本来考虑的是天线间时延
    
    delayTx = Tc * numTx; % 发射天线维度积累的时延
    arrdx = cfgOut.arrdx; % 归一化阵元间距
    arrdy = cfgOut.arrdy; 
    
    c = physconst('LightSpeed'; % 光速
    lambda = c / fc; % 波长
    Ts = 1 / fs; % ADC采样时间
    
    % 信号模型
    signal = zeros(ADCNum, TotalChirpNum, arrNum; % 信号
    
    noiseI = normrnd(0, 1, ADCNum, TotalChirpNum, arrNum; % 正态噪声
    noiseQ = normrnd(0, 1, ADCNum, TotalChirpNum, arrNum; % 正态噪声
%     noiseI = exprnd(1, ADCNum, TotalChirpNum, arrNum; % 指数分布噪声
%     noiseQ = exprnd(1, ADCNum, TotalChirpNum, arrNum; % 指数分布噪声
    for tarId = 1 : tarNum
        disp(strcat(['正在获取目标',num2str(tarId,'的数据'];
        % 目标参数
        targetR = tarOut.range(tarId;    % 距离 m
        targetV = tarOut.velocity(tarId; % 速度 m/s
        targetAzi = tarOut.Azi(tarId;    % 目标方位角 度
        targetEle = tarOut.Ele(tarId;    % 目标俯仰角 度
        targetRCS = tarOut.RCS(tarId;    % 目标RCS值
        targetGt  = tarOut.Gt(tarId;     % 发射增益 
        targetGr  = tarOut.Gr(tarId;     % 接收增益 
        [tarSNR] = CalculateSNR(targetR, targetRCS, targetGt, targetGr, lambda, Pt, Fn, Ls, validB; % 信噪比
        A = sqrt(2 * db2pow(tarSNR; % 信号幅度
        
        targPhi0 = rand * 2 * pi; % 产生随机相位[0 2*pi] 这里的目标初始相位也可以根据上面定义的形式来写 但影响不大

        tempSigI = zeros(ADCNum, TotalChirpNum, arrNum; % I路缓存信号
        tempSigQ = zeros(ADCNum, TotalChirpNum, arrNum; % Q路缓存信号
        tempSig = zeros(ADCNum, TotalChirpNum, arrNum; % I路缓存信号
        
        for channelId = 1 : arrNum
            for chirpId = 1 : TotalChirpNum
                for adcId = 1 : ADCNum
                    % 动目标时延
                    tarDelay = 2 * targetR / c + 2 * targetV * ((chirpId - 1 * delayTx + arrdelay(channelId + adcId * Ts / c; 
                    
                    % 控制来波方向
                    % 初阶版本的目标相位表达
%                     tarPhi = targPhi0 + 2 * pi * (arr(1, channelId * sind(targetAzi * arrdx + ...
%                             arr(2, channelId * sind(targetEle * arrdy + deg2rad(antennaPhase(channelId; % 考虑阵元初始、随机相位和目标波达角

%					tarPhi0 = 2 * targetR * fc  / c; % 目标初始相位 也可以用随机相位

                    tarPhi = targPhi0 + 2 * pi * (arr(channelId * sind(targetAzi * arrdx + ...
                            arr(arrNum+channelId * sind(targetEle * arrdy + deg2rad(antennaPhase(channelId; % 这里考虑了不同发射顺序天线相位

                    % 中频信号：exp[1j * 2 *pi * fc * tau + 2 * pi * S * t * tau - pi * tau * tau]  t的范围为0-Tc tau的范围为t+nTc
                    tempSigI(adcId, chirpId, channelId = A * cos(2 * pi * (fc * tarDelay + Slope * tarDelay * adcId * Ts - Slope * tarDelay * tarDelay / 2 + tarPhi;
                    tempSigQ(adcId, chirpId, channelId = A * sin(2 * pi * (fc * tarDelay + Slope * tarDelay * adcId * Ts - Slope * tarDelay * tarDelay / 2 + tarPhi;
                    tempSig(adcId, chirpId, channelId = tempSigI(adcId, chirpId, channelId + 1j * tempSigQ(adcId, chirpId, channelId;
                end
            end
        end
        signal = signal + tempSig; % 多目标信号相加
    end
    signal = signal - (noiseI + noiseQ; % 考虑噪声
%     signal = reshape(signal, size(signal,1, size(signal,2, numRx, numTx;
end

信号处理栈

在上一章节，我们已经得到了radar cube数据，这个数据和毫米波雷达ADC数据是异曲同工的。那么，针对现有数据可以展开信号处理栈的分析，简言之，传统的信号处理栈为RDA，如本文开篇的图所示。

radar cube数据的快时间和慢时间维度分别使用FFT（快速傅立叶变换）的到距离多普勒域数据RDD，这可以表示为如下形式：

\[x_n(\lambda,\zeta=\sum_{k_{tgt}=1}^{K_\text{tgt}}a_{k_{tgt}}{k_{tgt}}W_{\text{ft}}(\lambda-\lambda_{k_{tgt}}W_{\text{st}}(\zeta-\zeta_{k_{tgt}}e^{j\varphi_{k_{tgt}}}\cdot e^{j\frac{2\pi}{\lambda}(\mathbf{u}_{k\mathrm{tgt}}\cdot\mathbf{p}_n}+\Omega_n(\lambda,\zeta \]

为了进一步对距离-多普勒变换后的信号处理，将上述混合信号表示为如下形式：

radar cube的元素为\(x_{l,m,n}~~=~~x[l,m,n]\，快时间和慢时间维度的FFT点数对应\({P_{\mathrm{FFT}}}\和\({Q_{\mathrm{FFT}}}\

\[\lambda=\frac{2\pi}{P_{\mathrm{FFT}}}l,\quad\zeta=\frac{2\pi}{Q_{\mathrm{FFT}}}m,\quad x[l,m,n]=x_n\left(\frac{2\pi}{P_{\mathrm{FFT}}}l,\frac{2\pi}{Q_{\mathrm{FFT}}}m\right. \]

radar cube进行目标检测，这通常用恒虚警率检测算法(cfar, constant false alarm rate来实现，在正式进行检测前，我们需要对radar cube沿天线维度展开，通常可以选择一个天线通道对应的距离多普勒数据或使用所有天线通道的距离多普勒数据进行非相干积累处理（波束成形）得到cfar检测输入数据，这个过程可以表示如下：(下面第一式为方法1，第二式为方法2，\(N_{det}\为单天线通道，\(N^\prime\为全部天线通道

\[X^\prime=|x'[l,m]|^2=|x[l,m,n=N_{\det}]|^2 \]

\[X^\prime=|x'[l,m]|^2=\sum\limits_{n\in N'}|x[l,m,n]|^2 \]

cfar检测。通常，cfar检测针对数据维度可以分为1D-cfar和2D-cfar，针对处理思想可以分为均值类、有序统计类、自动筛选类和自适应类（这部分详细细节可以参考何友院士和关键老师的雷达目标检测与恒虚警处理书，后续会整理我的本科毕设中关于这方面的内容）。实际上，我们可以考虑两次1D-CFAR即对距离多普勒矩阵\(X^\prime\的距离和多普勒维度分别作cfar并关联匹配目标，或者考虑1次2d-cfar即直接在二维距离多普勒矩阵\(X^\prime\上作二维检测，这里我们主要用2d-cfar来检测目标。

2D-CFAR示意图来更好地理解原理。浅蓝色单元构成的大矩形区域为输入的距离-多普勒矩阵\(X^\prime\深红色单元为检测单元，浅红色单元为保护单元（防止旁瓣能量信号泄露主瓣形成掩蔽效应），浅绿色为参考单元，由深红色单元、浅红色单元和浅绿色单元构成的小矩形区域作为滑动窗窗；依次移动滑动窗在大矩形区域内遍历检测目标，移动过程中需要保证检测单元的中心从黄色的起始中心单元开始依次完成自左向右、自上而下的二维扫描，直至完成整个大矩形区域的遍历。【注意到，我们在这里强调了"从黄色的起始中心单元"，这部分需要在程序中重点处理，否则可能会导致部分区域是没有参与检测的。】

ca-cfar，cell-averaging cfar为例，首先需要根据预设虚警概率和滑动窗尺寸等先验参数计算门限因子\(\alpha\；在每次扫描中，需求求解参考窗内所有检测单元对应的信号电平的和并取平均值，这一估计值可以表示为\(\bar{X}\，将估计阈值\(\bar{X}\cdot \alpha\和检测单元下的信号电平进行比较，如果检测单元下信号电平大于估计阈值\(\bar{X}\cdot \alpha\，则认为检测单元出为目标信号，否则认为噪声信号。

2d-cfar过程：

function [cfarOut] = CFAR_2D(RDM, Pfa, TestCells, GuardCells
    % rd_map：输入数据
    % Pfa：虚警概率
    % TestCells：测试单元 [4 4]
    % GuardCells：保护单元 [4 4]
    % By Xuliang
    
    % CFAR检测器
    detector = phased.CFARDetector2D('TrainingBandSize',TestCells, ...
    'ThresholdFactor','Auto','GuardBandSize',GuardCells, ...
    'ProbabilityFalseAlarm',Pfa,'Method','SOCA','ThresholdOutputPort',true, 'NoisePowerOutputPort',true;

    N_x = size(RDM,1; % 距离单元
    N_y = size(RDM,2; % 多普勒单元
    % 获取二维滑动窗口的size
    Ngc = detector.GuardBandSize(2; % 保护窗列向
    Ngr = detector.GuardBandSize(1; % 保护窗行向
    Ntc = detector.TrainingBandSize(2; % 参考窗行向
    Ntr = detector.TrainingBandSize(1; % 参考窗行向
    cutidx = [];
    colstart = Ntc + Ngc + 1; % 列窗首
    colend = N_y + ( Ntc + Ngc; % 列窗尾
    rowstart = Ntr + Ngr + 1; % 行窗首
    rowend = N_x + ( Ntr + Ngr; % 行窗尾
    for m = colstart:colend
        for n = rowstart:rowend
            cutidx = [cutidx,[n;m]]; % 完整二维窗
        end
    end
    
    ncutcells = size(cutidx,2; % 获取窗口遍历区间
    
    rd_map_padding = repmat(RDM, 3, 3; % 对RDM补零
    chosen_rd_map = rd_map_padding(N_x+1-Ntr-Ngr:2*N_x+Ntr+Ngr,N_y+1-Ntc-Ngc:2*N_y+Ntc+Ngc;
    
    [dets, ~, noise] = detector(chosen_rd_map,cutidx; % 输出CFAR结果
    
    cfar_out = zeros(size(chosen_rd_map; % 检测点输出
    noise_out = zeros(size(chosen_rd_map; % 噪声
    snr_out = zeros(size(chosen_rd_map; % 信噪比
    for k = 1:size(dets,1
        if dets(k == 1
            cfar_out(cutidx(1,k,cutidx(2,k = dets(k; 
            noise_out(cutidx(1,k,cutidx(2,k = noise(k;
            snr_out(cutidx(1,k,cutidx(2,k = chosen_rd_map(cutidx(1,k,cutidx(2,k;
        end
    end
    
    cfarOut = {};
    cfarOut.cfarMap = cfar_out(Ntr+Ngr+1:Ntr+Ngr+N_x,Ntc+Ngc+1:Ntc+Ngc+N_y;
    cfarOut.snrOut = snr_out(Ntr+Ngr+1:Ntr+Ngr+N_x,Ntc+Ngc+1:Ntc+Ngc+N_y;
    cfarOut.noiseOut = noise_out(Ntr+Ngr+1:Ntr+Ngr+N_x,Ntc+Ngc+1:Ntc+Ngc+N_y;
    cfarOut.snrOut = cfarOut.snrOut ./ (eps + cfarOut.noiseOut;
end

经过cfar检测后，我们可以得到一组由\(K_{det}\个检测目标构成的集合：

\[\mathrm{cfar}\left(\left|x'[l,m]\right|^2\right=\{(r_1,v_{\mathrm{r},1},(r_2,v_{\mathrm{r},2},\ldots,(r_{k_{\det}},v_{\mathrm r,k_{\det}},\ldots,(r_{K_{\det}},v_{\mathrm r,K_{\det}}\} \]

下面，我们可以进入doa估计；不过我们需要准备好用于空间谱估计的输入信号，这需要用到由上面缓存的radar cube数据和经过cfar检测得到的目标距离-速度集合\(\mathrm{cfar}\left(\left|x'[l,m]\right|^2\right\来提取空域信号\(\mathbf{x}_{k\mathrm{det}}\，表示如下：（\(N\表示阵元数目）

\[\mathbf{x}_{k\mathrm{det}}=x[l_{k\mathrm{det}},m_{k\mathrm{det}},n]=[x_0,\ldots,x_{N-1}]^T \]

\[\mathbf{x}_{k_{\det}}=\sum_{k_{\text{doa}}=1}^{K_{\text{doa}}}s_{k_{\text{doa}}}\mathbf{a}(\phi_{k_{\text{doa}}},\epsilon_{k_{\text{doa}}}+\boldsymbol{\omega} \]

不难发现上述信号表示是标准的空间谱估计信号模型X=AS+N，其中A为流形，S为信号空间，N为噪声空间（上面用\(\omega\表示噪声项，通常假设满足复高斯白噪声分布），X为观测空间；在这里，我们也可以将A视为完备基或字典，字典原子或称导向矢量\(\mathbf{a}(\phi_{k_{\text{doa}}},\epsilon_{k_{\text{doa}}}\满足线性无关，那么目标的来波信号可以由字典中的原子通过复权重\(s_{k_{doa}}\加权得到，导向矢量和复权重的表示分布如下：

\[\mathbf{a}(\phi_{k_\text{doa}},\epsilon_{k_\text{doa}}=e^{j(2\pi/\lambda\mathbf{u}_{k_\text{doa}}\mathbf{p}_n} \]

\[s_{k\text{doa}}=a_{k\text{doa}}W_{\text{ft}}(\lambda-\lambda_{k\text{doa}}W_{\text{st}}(\zeta-\zeta_{k\text{doa}}e^{j\varphi_{k}\text{doa}} \]

doa估计的本意是在给定峰值目标空域信号\(\mathbf{x}_{k\mathrm{det}}\和阵列集合结构\(p_n\的基础上估计出目标数目\(K_{doa}\、目标来波方向\(\phi_{k_{doa}}\和\(\epsilon_{k_{doa}}\，即如下目标：\((r_{k_\mathrm{det}},v_{\mathrm{r},k_\mathrm{det}},\mathbf{x}_{k_\mathrm{det}}\to(r_{k_\mathrm{tgt}},v_{r,k_\mathrm{tgt}},\phi_{k_\mathrm{tgt}},\epsilon_{k_\mathrm{tgt}}\。在本处，我们可以对上述阵列信号完成对应维度的天线映射后，首先对方位天线维度做FFT处理，并认为方位维空间最大谱峰所在处为对应的目标方位角；其次，基于目标方位天线的峰值索引，我们沿其俯仰维度做FFT处理，并认为俯仰维空间最大谱峰所在处为对应的目标俯仰角，以此实现目标信号方位和俯仰维度的匹配关联。这种方法可以理解为“两步走”战略，和上面提及的两次1D-CFAR是类似的。

CBF、线性滤波Capon、子空间方法（MUSIC，ESPRIT，DML等）和压缩感知方法(OMP，IAA和LISTA等，这些内容我们将在后面的章节中讨论，在这里需要注意一点，我们要保证方位天线维度上的空间谱估计不能够破坏俯仰天线维度观测信号的相位信息，否则会导致俯仰角的估计错误。

RDA信号处理栈的基本物理量提取并最终得到点云表示\((r_{k_\mathrm{tgt}},v_{r,k_\mathrm{tgt}},\phi_{k_\mathrm{tgt}},\epsilon_{k_\mathrm{tgt}}\。所述过程仍存在一些细节性内容没有解释，这部分可以结合代码更好地去理解，下面给出核心信号处理实现部分。

for frame_id = 1 : Frame % Frame
    adcData = squeeze(RawData(:, :, frame_id, :;
    
    %% 距离维FFT和多普勒FFT
    disp(strcat(['=====','RD-MAP生成','=====']; 
    tic
    fftOut = rdFFT(adcData, IQFlag;
    toc
    rangeFFTOut = fftOut.rangeFFT;
    dopplerFFTOut = fftOut.dopplerFFT;
    if ShowRange
        figure(2;
        set(gcf,'unit','centimeters','position',[10,0,10,10]
        plot(rangeIndex, db(rangeFFTOut(:,1,1; 
        xlabel('Range(m';ylabel('Amplitude(dB'; 
        title(strcat(['第',num2str(frame_id,'帧-目标距离分布'];grid minor;  
        pause(0.1;
    end
    if ShowRD
        figure(3;
        set(gcf,'unit','centimeters','position',[20,12,10,10]
        imagesc(rangeIndex, velocityIndex, db(dopplerFFTOut(:,:,1.';
        xlabel('Range(m';ylabel('Velocity(m/s'; colormap('jet';
        title(strcat(['第',num2str(frame_id,'帧目标-距离多普勒分布'];
        grid minor; axis xy;
        pause(0.1;
    end

    %% 非相干积累
    disp(strcat(['=====','非相干积累','=====']; 
    RDM = dopplerFFTOut;
    tic
    [accumulateRD] = incoherent_accumulation(RDM;
    toc
    
    %% CFAR检测器
    disp(strcat(['=====','恒虚警检测','=====']; 
    Pfa = 1e-3; % 虚警概率
    TestCells = [8, 8]; % 参考窗
    GuardCells = [2, 2]; % 保护窗
    
    tic
    [cfarOut] = CFAR_2D(accumulateRD, Pfa, TestCells, GuardCells;
    toc
    cfarMap = cfarOut.cfarMap; % 检测点输出
    noiseOut = cfarOut.noiseOut; % 噪声输出
    snrOut = cfarOut.snrOut; % 信噪比输出
    
    if ShowCFAR
        figure(4;
        set(gcf,'unit','centimeters','position',[20,0,10,10]
        imagesc(rangeIndex, velocityIndex, cfarMap.';
        xlabel('Range(m';ylabel('Velocity(m/s'; colormap('jet';
        title(strcat(['第',num2str(frame_id,'帧目标-CFAR检测结果'];
        grid minor;axis xy;
        pause(0.1;
    end
    
    %% 峰值聚合-获取点目标
    disp(strcat(['=====','峰值聚合','=====']; 
    [range_idx, doppler_idx] = find(cfarMap;
    cfar_out_idx = [range_idx doppler_idx]; % 获取CFAR输出结果的行列索引
    tic
    [rd_peak_list, rd_peak] = peakFocus(db(accumulateRD, cfar_out_idx;
    toc
    peakMap = zeros(size(cfarMap; % 峰值聚合结果矩阵
    for peak_idx = 1 :size(rd_peak_list, 2
        peakMap(rd_peak_list(1,peak_idx, rd_peak_list(2,peak_idx = 1;
    end
    
    if ShowPeak
        figure(5;
        set(gcf,'unit','centimeters','position',[30,12,10,10]
        imagesc(rangeIndex, velocityIndex, peakMap.';
        xlabel('Range(m';ylabel('Velocity(m/s'; colormap('jet';
        title(strcat(['第',num2str(frame_id,'帧目标-峰值聚合结果'];
        grid minor;axis xy;
        pause(0.1;
    end

    %% DOA/AOA估计
    disp(strcat(['=====','DOA/AOA估计','=====']; 
    cfgDOA.FFTNum = 180; % FFT点数
    % 这里虽然封装了不同的DOA算法 但是需要注意的是 可用的算法有限 在本套代码里 建议使用的是FFT-FFT和FFT-MUSIC等
    % 因为MUSIC-MUSIC的使用会导致在方位维度空间谱估计时破坏相位信号 有损俯仰维度的相位估计
    
    cfgDOA.AziMethod = 'FFT'; % 方位维度DOA估计方法
    cfgDOA.EleMethod = 'MUSIC'; % 俯仰维度DOA估计方法
    
    cfgDOA.thetaGrids = linspace(-90, 90, cfgDOA.FFTNum; % 空间网格
    cfgDOA.AzisigNum = 1; % 约束每个RD-CELL上的信源数目
    cfgDOA.ElesigNum = 1; % 约束每个方位谱峰上的信源数目
    
    aziNum = length(virtual_array.noredundant_aziarr; % 方位天线数目
    
%     Aset = exp(-1j * 2 * pi * arrDx * [0:aziNum]' * sind(cfgDOA.thetaGrids; % 稀疏字典设计
    targetPerFrame = {}; 
    targetPerFrame.rangeSet = [];
    targetPerFrame.velocitySet = [];
    targetPerFrame.snrSet = [];
    targetPerFrame.azimuthSet = [];
    targetPerFrame.elevationSet = [];
    
    if ~isempty(rd_peak_list % 非空表示检测到目标
        rangeVal = (rd_peak_list(1, : - 1 * range_res; % 目标距离
        speedVal = (rd_peak_list(2, : - ChirpNum / 2 - 1 * doppler_res; % 目标速度
        
        doaInput = zeros(size(rd_peak_list, 2, arrNum;
        for tar_idx = 1 :size(rd_peak_list, 2
            doaInput(tar_idx, : = squeeze(dopplerFFTOut(rd_peak_list(1, tar_idx, rd_peak_list(2, tar_idx, :; % tarNum * arrNum
        end
        doaInput = reshape(doaInput, [], numRx, numTx;
        
        % 方位角估计前需要考虑多普勒补偿
        [com_dopplerFFTOut] = compensate_doppler(doaInput, cfgOut, rd_peak_list(2, :, speedVal, rangeVal; 
        
        tic
        for peak_idx = 1 : size(rd_peak_list, 2 % 遍历检测到的每个目标
            snrVal = mag2db(snrOut(rd_peak_list(1, peak_idx, rd_peak_list(2, peak_idx; % 信噪比的提升是由于chirpNum*ADCNum的积累
            tarData = squeeze(com_dopplerFFTOut(peak_idx, :,:;

    %         aziarr = unique(arr(1,arr(2,:==0; % 初阶版本获取方位维度天线排列
    %         aziArrData = arrData(aziarr+1; % 获取方位维度信号

            % 方位角解析
           sig = tarData;
           sig_space = zeros(max(virtual_array.azi_arr+1,max(virtual_array.ele_arr+1; % 初始化信号子空间
           for trx_id = 1 : size(cfgOut.sigIdx,2
               sig_space(cfgOut.sigSpaceIdx(1, trx_id, cfgOut.sigSpaceIdx(2,trx_id = sig(cfgOut.sigIdx(1,trx_id, cfgOut.sigIdx(2,trx_id; % 重排后的信号空间
           end
           % 至此我们生成的信号子空间维度为 方位虚拟天线数目 * 俯仰虚拟天线数目 

           eleArrData = zeros(cfgDOA.FFTNum, size(sig_space,2; % 俯仰维度数据
            for ele_idx = 1 : size(sig_space, 2 % 这里采取遍历是为了适配不同的空间谱估计方法
                tmpAziData = sig_space(:, ele_idx;
                [azidoaOut] = azimuthDOA(tmpAziData, cfgDOA; % 提取第一列方位维度天线信息进行方位角估计 
                eleArrData(:, ele_idx = azidoaOut.spectrum(:; % 复空间谱
            end

            for azi_peak_idx = 1 : length(azidoaOut.angleVal % 对方位维度检测的谱峰进行检索
                tmpEleData = eleArrData(azidoaOut.angleIdx(azi_peak_idx, :.'; % 获取与方位维目标关联的信号
                [eledoaOut] = elevationDOA(tmpEleData, cfgDOA; % 进行俯仰角估计

                % 关联目标的距离、多普勒、信噪比、方位和俯仰信息
                aziVal = azidoaOut.angleVal; 
                eleVal = eledoaOut.angleVal;
                targetPerFrame.rangeSet = [targetPerFrame.rangeSet, rangeVal(peak_idx];
                targetPerFrame.velocitySet = [targetPerFrame.velocitySet, speedVal(peak_idx];
                targetPerFrame.snrSet = [targetPerFrame.snrSet, snrVal];
                targetPerFrame.azimuthSet = [targetPerFrame.azimuthSet,aziVal];
                targetPerFrame.elevationSet = [targetPerFrame.elevationSet,eleVal];
            end
        end
        toc
    
        %% 点云生成 
        disp(strcat(['=====','点云生成','=====']; 
        tic
        pcd_x = targetPerFrame.rangeSet .* cosd(targetPerFrame.elevationSet .* sind(targetPerFrame.azimuthSet;
        pcd_y = targetPerFrame.rangeSet .* cosd(targetPerFrame.elevationSet .* cosd(targetPerFrame.azimuthSet;
        pcd_z = targetPerFrame.rangeSet .* sind(targetPerFrame.elevationSet;
        PointSet = [pcd_x.', pcd_y.', pcd_z.', targetPerFrame.velocitySet.', targetPerFrame.snrSet.'];
        toc

        %% 点云聚类
        eps = 1.1; % 邻域半径
        minPointsInCluster = 3; % 簇内最小点数阈值
        xFactor = 1;   % 变大控制距离变大，变小分类距离变小 椭圆
        yFactor = 1;   % 变大控制角度变大，变小分类距离变小 椭圆 
        figure(6;
        set(gcf,'unit','centimeters','position',[30,0,10,10]
        disp(strcat(['=====','点云聚类','=====']; 
        tic
        [sumClu] = dbscanClustering(eps,PointSet,xFactor,yFactor,minPointsInCluster,frame_id; % DBSCAN聚类
        toc
    end

end

在本文的尾巴，本文将展示本程序所实现的可视化功能。在以下仿真中，目标距离为[5 10 15 18]，目标速度为[0.4 -0.3 0 0.6]，目标方位角为[15 -2 30 -25]，目标俯仰角为[0 15 -5 8]，目标RCS为[20 20 10 20]，仿真实验平台为TI AWR1243BOOST，其他平台均可实现类似效果。

本完整程序可以通过邮箱形式向笔者咨询，希望读者以诚心求学心态咨询。

参考文献

[2] J. Fuchs, M. Gardill, M. Lübke, A. Dubey and F. Lurz, "A Machine Learning Perspective on Automotive Radar Direction of Arrival Estimation," in IEEE Access, vol. 10, pp. 6775-6797, 2022, doi: 10.1109/ACCESS.2022.3141587。

[3] 初出茅庐：毫米波雷达开发手册之基础参数 https://www.cnblogs.com/yuxuliang/p/MyRadar_4.html

编程笔记 » 登堂入室:毫米波雷达开发手册之信号模型