登堂入室：毫米波雷達開發手冊之訊號模型
寫在前面
深知新手在接觸毫米波雷達板硬體時需要花費的沉沒成本，因此在行將告別毫米波雷達之際，總結這兩年以來在毫米波雷達上的一些經驗和教訓。
本檔案用於為實現基於AWR1243BOOST等單板毫米波雷達開發提供參考指南與解決方案，主要包括硬體設定、基礎引數、訊號模型、應用DEMO開發以及可深入研究方向思考等；為更好地匹配後續級聯雷達應用的學習路線，在本手冊中會盡可能同化單板雷達和級聯雷達中的相關表述。
本指南作者資訊：Xuliang，聯絡方式：[email protected]。未經本人允許，請勿用於商業和學術用途。
希望後者在使用本指南時可以考慮參照作者在毫米波雷達旅途中的相關工作，如本文參考文獻[1].
本章節為訊號模型章節，主要解讀時域訊號模型、空域訊號模型和訊號處理棧。
往期內容：
眼觀四海：自動駕駛&4D成像毫米波雷達 如今幾何？
揚帆起航：毫米波雷達開發手冊之硬體設定
初出茅廬：毫米波雷達開發手冊之基礎引數
訊號模型分析
時域訊號模型
如圖所示，雷達收發器（TRX）目前通常被實現為單個單片微波積體電路（MMIC）或甚至包括數位功能的片上系統（SoC），其生成雷達發射（TX）波形並經由連線到其TX通道的天線輻射該波形。電磁波通過空間傳播，直到它命中一個或多個目標。假設目標的各向同性輻射，與目標雷達截面（RCS）成比例的入射波通量密度的一部分被捕獲並朝向雷達再輻射。根據經典理論，目標通常被假設為孤立的點狀散射體即點目標，如汽車、行人或環境之類的物理物件由若干點狀目標組成，並且各個目標響應的總和則是總體接收（RX）訊號。
通過比較雷達的TX和RX訊號，可以估計出感知任務的基本引數：目標數目\(K_{tgt}\)和每個目標\(k_{tgt}\)的目標距離\(r_{k_{tgt}}\)，目標徑向速度\(v_{r_{k_{tgt}}}\)，方位角\(\phi_{k}\)和俯仰角\(\epsilon_k\)。在每個測量週期中檢測到的目標集合被稱為目標列表，並且可以被解釋為從雷達測量獲得的點雲。後續的毫米波雷達訊號處理步驟主要有是聚類（即識別屬於同一物理目標）、跟蹤（即從目標列表的序列中匯出物理上存在的物件的資訊和統計資料）。本章節的核心是介紹TDMA-MIMO框架下的訊號模型模擬及實現方法。
本文所研究的模型與快速線性調頻調頻連續波（FMCW）雷達相關，這是目前智慧駕駛、智慧家居領域中所使用的主要雷達波形。如下圖所示為典型的由Q個chirp組成的chirp序列構成的發射波形（綠色部分）和來自單個目標的接受波形訊號（藍色部分）。
發射訊號\(s（t）\)由以chirp-chirp間隔的為\(T_{cc}\)的\(Q\)個線性調頻脈衝波形\(s_c（t）\)的序列給出：
其中脈衝持續時間為\(T_c\)的每個chirp的訊號表示為：
其中，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)\)。
每個點目標作為發射訊號\(s(t)\)的時延副本以加權形式參與到接收訊號\(r(t)\)中，其表示形式如下：
其中，目標衰減係數\(a_{k_{tgt}}\)是由於自由空間傳播和目標RCS引起的係數加權，往返時間\(T_{d,k}(t)\)是時間相關的，因為這考慮了目標的移動，通常幀週期為幾毫秒或更小，與目標速度相比較短，無法形成積累效應。因此，\(T_{\mathrm{d},k}(t))\)可以近似表示為：
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\)所在的維度稱為慢時間維度，這在上一節內容中已經提及。
基於所給的調變引數，可以將二維矩陣\({\mathcal{B}}\)的每個元素\(b_{p,q}\)展開如下：
其中，對應距離維度的歸一化角頻率可表示如下：
對應多普勒維度的歸一化角頻率可表示如下：
目標的初始相位則可表示如下：
空域訊號模型
為了建立與目標方位角和仰角的關係，可以使用\(N_{RX}\)個接收天線的陣列，這將被稱為單輸入多輸出（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]\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].\)$
假設目標滿足遠場窄頻訊號假設，那麼可進一步擴充套件時域訊號至空域訊號如下：
上述訊號對應於雷達3D tensor或radar cube資料的元素，根據ISO8000:2標準定義的球面座標系，目標\(k_{tgt}\)對應的單位長度波達角向量可以表示如下：
其中，\(\phi_{k_\textrm{tgt}}\)和\(\theta_{k_\textrm{tgt}}\)分別表示目標的方位角和傾斜角，俯仰角則表示為\(\epsilon_{k_{tgt}}=\pi-\theta_{k_{tgt}}\)；虛擬天線的位置\(p_n\)可以表示如下（\(n=0....N-1\)）
至此，我們可以很明確地建立目標角度和模型\(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，這可以表示為如下形式：
其中\(W_{ft}(\lambda)\)和\(W_{st}(\zeta)\)分別表示在快時間和慢時間維度上應用的窗函數傅立葉變換，這實現了訊號移位至目標對應的距離維和多普勒峰值頻率處：\(\lambda_{k_{tgt}}\)和\(\zeta_{k_{tgt}}\)。目標初始相位項\(e^{j\varphi_{k_{tgt}}}\)和相位分佈項\(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}}}\)
其次，我們需要對距離-多普勒域的radar cube進行目標檢測，這通常用恆虛警率檢測演演算法(cfar, constant false alarm rate)來實現，在正式進行檢測前，我們需要對radar cube沿天線維度展開，通常可以選擇一個天線通道對應的距離多普勒資料或使用所有天線通道的距離多普勒資料進行非相干積累處理（波束成形）得到cfar檢測輸入資料，這個過程可以表示如下：(下面第一式為方法1，第二式為方法2，\(N_{det}\)為單天線通道，\(N^\prime\)為全部天線通道)
那麼，現在我們可以對處理後的距離多普勒矩陣\(X^\prime\)進行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}\)個檢測目標構成的集合：
其中，每個檢測目標由與其關聯的引數元組進行表示，每個引數元組\((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\)表示陣元數目）
基於經過距離多普勒FFT的混合訊號模型\(x_n(\lambda,\zeta)\)，我們可以建立空域訊號模型為：
不難發現上述訊號表示是標準的空間譜估計訊號模型X=AS+N，其中A為流形，S為訊號空間，N為噪聲空間（上面用\(\omega\)表示噪聲項，通常假設滿足復高斯白噪聲分佈），X為觀測空間；在這裡，我們也可以將A視為完備基或字典，字典原子或稱導向向量\(\mathbf{a}(\phi_{k_{\text{doa}}},\epsilon_{k_{\text{doa}}})\)滿足線性無關，那麼目標的來波訊號可以由字典中的原子通過復權重\(s_{k_{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，其他平臺均可實現類似效果。
本完整程式可以通過郵箱形式向筆者諮詢，希望讀者以誠心求學心態諮詢。
參考文獻
[1] X. Yu, Z. Cao, Z. Wu, C. Song, J. Zhu and Z. Xu, "A Novel Potential Drowning Detection System Based on Millimeter-Wave Radar," 2022 17th International Conference on Control, Automation, Robotics and Vision (ICARCV), Singapore, Singapore, 2022, pp. 659-664, doi: 10.1109/ICARCV57592.2022.10004245.
[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