这几年做雷达信号处理被问得最多的不是波束赋形这类进阶方向反而是FMCW雷达的基础链路波形生成、回波仿真、距离多普勒变换、目标检测。这几个模块拆开看都不算难可要串成一条能跑的MATLAB链路尤其是把CA-CFAR单元平均恒虚警率检测这一步做对、做得稳确实需要不少细节。这篇文章我就把一套从波形生成到CA-CFAR目标检测的完整仿真实现拿出来讲透包括我当时调参踩过的坑。适合刚接触雷达仿真、或者已经能画出距离多普勒图但不知道下一步怎么做的同学参考。1. 先搞明白FMCW雷达为什么能同时测距测速1.1 从线性调频chirp说起FMCW雷达发射的不是恒定频率的正弦波而是一个频率随时间线性变化的信号也就是常说的chirp。最常见的调制波形是锯齿波频率从 f_c 开始线性上升到 f_c B到达chirp末端后瞬间跳回起始频率进入下一个周期。发射信号可以写成s_tx(t) exp( j*2*pi*( f_c*t 0.5*K*t^2 ) )其中 K B / T_c 是调频斜率B是扫频带宽T_c是单个chirp的时长。为什么要让频率随时间变化因为连续波没有距离信息只有让发射信号带上时间标签才能从回波信号里反推出目标距离。逻辑上可以类比成两个人同时吹哨子一个人听到回声如果哨音是恒定频率我们只能知道回声存在却不知道它晚到了多久可如果哨音从低到高连续变化那晚到就直接体现在音调的高低差上这个高低差就是距离的度量。1.2 混频出来的拍频藏着距离发射信号遇到目标反射回来雷达接收到的回波是发射信号的延迟版本。假设目标静止距离为R那么回波延迟约等于tau 2R / c本振信号发射信号的一份副本与回波信号做混频再滤掉高频分量得到中频IF信号。混频的本质是两信号相乘取低频分量所以IF信号的频率恰好是发射频率和回波频率之差。由于频率是线性变化的这个频率差就是f_b K * tau K * 2R / c 2*B*R / (c*T_c)举个例子假设载频 f_c 77 GHz带宽 B 200 MHzchirp时长 T_c 20 μs采样率 F_s 20 MHz。目标距离R 50 m时K 200e6 / 20e-6 1e13 Hz/s tau 2*50 / 3e8 3.33e-7 s f_b 1e13 * 3.33e-7 3.33 MHz这个3.33 MHz的频率落在基带采样范围10 MHz以内可以在接收端直接用ADC采样。IF信号的频率在线性调制下和距离严格成正比所以只要对IF信号做FFT找到峰值对应的频率再乘一个比例系数就是目标的距离。1.3 运动目标慢时间维上出现了多普勒如果目标在径向运动回波延迟本身就在随时间变化。对于第n个chirp目标位置近似为R_n R0 v * n * T_c于是IF信号的相位里出现了关键的一项phi 2*pi*f_c * tau 4*pi*(R0 v*n*T_c) / lambda这里的 lambda c / f_c 是载波波长。相邻chirp之间的相位变化量是delta_phi 4*pi*v*T_c / lambda对应的多普勒频率为f_d 2*v / lambda所以问题就从一维变成了二维在一个chirp内部快时间维上IF频率只和距离有关在chirp与chirp之间的慢时间维上相位变化率只和多普勒也就是速度有关。一段完整的处理逻辑是先对每个chirp做距离FFT再把同一距离bin上所有chirp的复数输出排列起来做第二次FFT得到二维距离-多普勒矩阵。这就是所谓的2D-FFT。1.4 先算一算分辨率心里有个底仿真之前最好先把手算结果写在草稿上后面核对FFT输出用得上。距离分辨率ΔR c / (2B) 3e8 / (2*200e6) 0.75 m最大不模糊距离R_max (F_s/2) * c * T_c / (2B) 10e6 * 3e8 * 20e-6 / (2*200e6) 150 m速度分辨率Δv lambda / (2N_chirpT_c)。N_chirp128时Δv ≈ 0.0039 / (212820e-6) ≈ 0.76 m/s最大不模糊速度v_max lambda / (4*T_c) ≈ 48.7 m/s这些数字决定了目标在频谱上占几个bin也决定了后面CFAR窗口怎么选。2. 仿真系统搭建参数、回波模型与信号链路2.1 从指标倒推波形参数很多初学者拿到题目就写代码结果带宽、采样率都是随手填的最后检测效果对不上也不知道问题在哪。正确的做法是从指标倒推参数。设计目标无非几个距离分辨率要多少、最大探测距离多少、最大速度多少、速度分辨率多少。先定这四个指标再往回推波形参数。一套典型的车载雷达仿真参数如下参数符号取值备注载频f_c77 GHz常见车载雷达频段扫频带宽B200 MHz决定距离分辨率0.75 mchirp时长T_c20 μs决定最大速度约48.7 m/sADC采样率F_s20 MHz决定最大距离150 mchirp数/帧N_chirp128决定速度分辨率0.76 m/s每chirp采样点数N_s400T_c * F_s这套参数完全够做一个演示级的目标检测仿真距离分辨率0.75 m速度分辨率0.76 m/s最大距离150 m最大速度约48.7 m/s覆盖日常道路场景。2.2 中频信号建模的两种思路建回波的时候有两种做法。第一种是完整生成射频发射信号和回波信号再做混频第二种是直接跳过射频按照前面推导的中频相位表达式生成IF信号。仿真中我一般都选第二种因为射频混频带来的高载频分量在离散仿真里既费算力又没有额外信息。中频信号可以写成s_if(t, n) sum_k A_k * exp( j*2*pi*( f_c*tau_k K*tau_k*t ) )tau_k是第k个目标在第n个chirp的延迟。注意这个表达式已经包含了慢时间维的多普勒相位 4piR_n/lambda因为f_ctau_k展开后就有 vn*T_c项。下面这段MATLAB代码生成两个目标的回波数据并叠加高斯白噪声% FMCW雷达中频信号仿真 clear; close all; fc 77e9; % 载频 B 200e6; % 带宽 Tc 20e-6; % chirp时长 Fs 20e6; % ADC采样率 Nc 128; % chirp数 Ns round(Tc * Fs); % 单chirp采样点数 K B / Tc; % 调频斜率 c 3e8; lambda c / fc; % 目标: 距离(m), 径向速度(m/s), 幅度 targets [50, 10, 1; 80, -5, 0.35]; t_fast (0:Ns-1). / Fs; ifMat zeros(Ns, Nc); for n 0:Nc-1 for k 1:size(targets, 1) R targets(k,1) targets(k,2) * n * Tc; % 第n个chirp时目标距离 tau 2 * R / c; % 回波延迟 ifMat(:, n1) ifMat(:, n1) ... targets(k,3) * exp(1j*2*pi*(fc*tau K*tau*t_fast)); end end % 加噪: 噪声按纯噪声底归一化 amp sqrt(mean(abs(ifMat(:)).^2)); ifMat ifMat amp * 0.25 * (randn(Ns, Nc) 1j*randn(Ns, Nc));这里第二个目标幅度只有第一个的约三分之一模拟一个弱目标。跑通之后可以把这个幅度再压低观察CFAR漏检现象后面会专门讨论。2.3 混频后的处理起点有了ifMat后面所有处理都是标准的距离维FFT、多普勒维FFT、求功率谱、CFAR检测。需要提醒的是仿真信号是复数基带所以功率谱应该用abs(fft(...)).^2不要取实部。还有一个新手容易忽略的点混频后理论上应该滤除和频分量但在离散基带仿真里直接按中频表达式生成数据已经隐含了低通滤波不需要额外处理。如果希望噪声更贴近真实接收机也可以用awgn(ifMat, snr_ref, measured)来定量控制输入信噪比方便后面画检测性能随SNR变化的曲线。3. 距离-多普勒处理从一维FFT到二维矩阵3.1 距离维FFT加窗、补零、坐标映射距离维FFT是对每个chirp单独做的。原始数据有400个采样点直接做400点FFT的频率分辨率为50 kHz对应距离分辨率0.75 m这是理论极限。不过为了峰值更平滑、画图更漂亮、给后续插值留空间我习惯补零到512点再做FFT。注意补零不改变分辨率只改变插值密度别被补零能提高分辨率的说法误导。窗函数方面不加窗的矩形窗旁瓣很高强目标的旁瓣可能压在另一个目标的检测区域里。仿真里我通常用hann窗等效降低旁瓣约19 dB代价是主瓣略宽一点。对目标检测来说这点宽度完全可接受。NFFT_R 512; win_r hann(Ns); % 距离维FFT: 结果shape (NFFT_R, Nc) rd fft(ifMat .* win_r, NFFT_R, 1);这里不上fftshift因为距离轴从0开始正频率在前半段就够用了。注意加窗要在同一列上向量乘以整列MATLAB中ifMat .* win_r利用了隐式扩展每行乘同一个因子正确。3.2 速度维FFT慢时间维的第二次变换多普勒维FFT在慢时间维上进行。此时输入矩阵的每一行是一个固定距离bin上、跨128个chirp的复数序列对这个序列做FFT得到多普勒频率。因为多普勒有正有负靠近为正向远离为负向所以这次要对输出做fftshift把零频移到中间。NFFT_D 512; win_d hann(Nc).; rd rd .* win_d; % 多普勒维加窗 RD fftshift(fft(rd, NFFT_D, 2), 2);注意hann(Nc).转成行向量保证矩阵右乘时按列广播。3.3 坐标映射与RD图可视化FFT做完之后必须把bin索引换算成真实的距离和速度否则很难判断仿真是否正常。换算公式range_axis (0:NFFT_R-1) * Fs / NFFT_R; % 拍频频率轴 range_axis range_axis * c * Tc / (2*B); % 拍频频率 - 距离 doppler_axis (-NFFT_D/2:NFFT_D/2-1) * (1/(NFFT_D*Tc)); % 多普勒频率轴 velocity_axis doppler_axis * lambda / 2; % 多普勒频率 - 速度可视化时用dB功率谱powerMap abs(RD).^2; powerMap_dB 10*log10(powerMap / max(powerMap(:))); imagesc(velocity_axis, range_axis, powerMap_dB); xlabel(速度 (m/s)); ylabel(距离 (m)); colorbar; clim([-60 0]);仿真参数下50 m、10 m/s的目标应该出现在距离轴约50 m、速度轴约10 m/s的位置80 m、-5 m/s的目标出现在距离轴约80 m、速度轴约-5 m/s的位置。如果峰值位置和设计对不上先回头检查坐标公式这是排查问题最快的一步。另一个实用操作是静态杂波对消。如果在实际场景中存在固定强反射物如路面、护栏距离多普勒图上会有一条沿多普勒零频贯穿的亮线目标可能被它淹没。可以在距离FFT之后、多普勒FFT之前相邻chirp做一次差分即用后一个chirp减去前一个chirp把静止分量消掉。仿真演示里一般不加但画RD图时发现零频附近有异常高亮优先怀疑这个。4. CA-CFAR检测的工程实现阈值如何自适应噪声4.1 固定阈值为什么不够稳看过RD图的人都知道噪声底不是平的。靠近雷达的近距区可能有强反射远距区噪声功率下降加上窗函数旁瓣和杂波整体底噪起伏很大。如果用一个全局固定阈值阈值低了近距杂波区会爆出一堆虚警阈值高了远距弱目标又检不出来。CFAR的思路就是不再使用全局阈值而是就地取材用待检测单元周围一片区域的噪声平均值来估计当前噪声底再乘一个系数形成检测门限。这样阈值可以跟随局部噪声底波动远处的弱目标和近处的强目标得以用同一条逻辑处理。4.2 CUT、保护单元、参考单元CA-CFAR的全称是单元平均恒虚警率。每个待检测单元称作CUTCell Under Test。在CUT周围取一个矩形邻域紧挨着CUT的一圈是保护单元guard cells保护单元外面一圈是参考单元reference cells。保护单元的存在是为了防止目标主瓣或旁瓣能量污染噪声估计。一个目标在距离维和速度维都可能占据多个bin如果这些bin落入参考单元会把噪声均值拉高阈值抬升最终导致目标自己被吃掉。二维CA-CFAR的窗口结构可以这样理解设参考单元半径为 Rr保护单元半径为 Gr那么整个窗口大小是 (2Rr1)×(2Rr1)中心区域 (2Gr1)×(2Gr1) 不参与均值估计。用我前面那组参数计算当 Rr6、Gr3时窗口是13×13保护区域7×7参考单元总数Nref (2*61)^2 - (2*31)^2 169 - 49 120阈值系数α和虚警概率Pfa的关系为Pfa (1 alpha/Nref)^(-Nref)反过来求αalpha Nref * (Pfa^(-1/Nref) - 1)取 Pfa 1e-6Nref 120alpha 120 * ( (1e-6)^(-1/120) - 1 ) 120 * (10^(6/120) - 1) ≈ 120 * (1.122 - 1) ≈ 14.6意思是局部噪声均值大约要放大14.6倍作为检测门限。参考单元数越多α越接近一个小常数但估计方差越大参考单元数越少α越大但估计更离散。这个公式是CA-CFAR的把心建议记下来以后换参考窗大小可以直接算。4.3 用卷积快速实现二维CA-CFAR逐循环遍历每个CUT再求窗口内的均值在MATLAB里非常慢128×512的RD图就要几万次嵌套运算实测慢得没法用。工程上有个讨巧的办法参考单元求和本质上是一个box filter可以用卷积一次算完。思路是构造一个和参考窗同样大小的模板矩阵保护区域和CUT位置填0其余位置填1。把功率图与模板做卷积得到的就是每个CUT周围参考单元的和再除以模板中1的数量就是均值。function detMap caCFAR2D(powerMap, Rr, Gr, pfa) % 二维单元平均恒虚警率检测 % powerMap: 距离-多普勒功率谱 % Rr: 参考单元半径Gr: 保护单元半径pfa: 虚警概率 [nR, nD] size(powerMap); win ones(2*Rr1, 2*Rr1); win(Rr1-Gr:Rr1Gr, Rr1-Gr:Rr1Gr) 0; win(Rr1, Rr1) 0; % CUT本身不参与 Nref sum(win(:)); alpha Nref * (pfa^(-1/Nref) - 1); % 边界填充减少边缘效应 padMap padarray(powerMap, [Rr Rr], replicate); localSum conv2(padMap, win, valid); localMean localSum / Nref; threshold alpha * localMean; % 边界部分截掉 detMap powerMap(Rr1:end-Rr, Rr1:end-Rr) threshold; end边界处理是第一坑。直接对原图conv2用same边缘处只有部分参考单元参与平均算出来的localMean偏小阈值被压低边界会冒出一圈虚警。上面代码用padarray的replicate填充边缘再用valid卷积相当于把边缘当作邻居区域噪声继续延伸来处理能明显减少边缘虚警。4.4 参数选择与典型输出参考窗口大小的选择原则参考单元区域要足够包含噪声采样但又不能大到跨越RD图中的不均匀地形。保护区域则要覆盖目标可能占用的最大扩展范围。我常用的起点是场景参考半径Rr保护半径Gr参考单元数NrefPfaα稀疏目标631201e-614.6密集目标42641e-512.5强杂波842081e-410.7α随Nref和Pfa的变化规律Nref从64增加到208α反而从17.2降到14.2左右Pfa从1e-6放宽到1e-4α明显下降。仿真中Pfa选1e-6属于比较保守的设置能有效压制噪声毛刺。如果发现虚警太多先别急着调α检查参考单元里是不是混进了目标旁瓣如果漏检弱目标看看保护单元是不是太大了。在我的仿真参数下RD图功率谱上两个目标峰值都远高于局部噪声均值CA-CFAR检测结果会精确落在两个峰值位置附近。若把第二个目标幅度从0.35降到0.1它大概率会被漏检因为目标周围的参考单元可能已经包含了第一个目标的强能量这就是下一章要展开的遮蔽效应。5. 从CFAR检测点到目标量测峰值提取与性能局限5.1 检出来不等于量测好CA-CFAR输出的是一幅二值图像1代表检测到目标0代表没有。但检测到的点往往不止一个一个强目标的主瓣检测结果会在RD图上形成一小片连通区域如果直接取第一个1的位置作为目标距离速度会带来约半个bin的误差。更稳的做法是第一步用连通域标记区分不同目标簇cc bwconncomp(detMap);第二步对每个连通域用功率图加权求质心相当于在检测区域内按能量重心精确定位S regionprops(cc, powerMap_dB, WeightedCentroid);第三步把质心的像素坐标换算到距离和速度轴。这套检测聚类质心流程在密集目标场景下效果比取最大值稳定得多。曾经有个多目标场景两个目标相距不足两个距离bin直接取峰值的输出在目标1和目标2之间来回跳换成质心后就稳定了。如果想让精度再上一层可以对峰值附近做抛物线插值。距离维FFT在峰顶附近近似抛物线取峰值点及其左右各一个点的幅度代入三点抛物线公式求顶点位置可以把距离估计做到亚bin精度。实测中当目标静止时插值后距离抖动可以从0.6 m降到0.15 m左右代价只是几行代码。5.2 CA-CFAR的遮蔽效应是先天缺陷CA-CFAR有一个广为人知的缺陷当同一参考窗内存在两个功率相差很大的目标时强目标会把局部噪声均值抬高导致弱目标阈值被抬得过高而漏检这就是遮蔽效应。用前面仿真参数做实验第一个目标幅度1第二个目标幅度0.1两者在RD图上相距约20个距离bin恰好都在同一个参考窗范围内。检测结果显示强目标稳定被检出弱目标几乎无法通过检测。这被称为强目标掩盖弱目标。如果需要在这种情况下保证弱目标可检就要换用GO-CFAR或SO-CFAR。GO-CFAR把左右两侧参考窗分开估计取两者较大值作为阈值应对杂波边缘更稳但弱目标场景更吃亏SO-CFAR取较小值抗遮蔽能力强但均匀噪声下虚警率会比名义Pfa略高。选型口诀很简单杂波边缘多就用GO多目标密集就考虑SO均匀噪声且目标稀疏用CA就够了。5.3 边缘虚警的最后防线即使做了replicate填充RD图最外层仍然可能有零星虚警尤其是多普勒轴的两端。这是因为FFT边缘旁瓣叠加噪声后局部统计特性确实与非边缘不同。遇到这种情况我通常直接裁剪掉边缘区域只对距离轴4个bin和速度轴4个bin以内的有效区域做目标输出。仿真本来就是要看核心区域性能没必要让边缘的虚假点污染结果。代码里可以加一句话validRegion false(size(detMap)); validRegion(5:end-4, 5:end-4) true; detMap detMap validRegion;这个操作虽然粗暴但在工程输出目标列表之前过滤掉边缘虚警比事后人工剔点高效很多。6. 仿真之外参数调优顺序与我的几点体会6.1 调参顺序决定了排查效率我自己的经验是先固定波形参数把它当作不可动摇的地基再调CFAR参数。波形参数决定RD图上目标落在哪里、噪声底长什么样CFAR参数只是在这张图上画阈值。如果波形参数还没验过就先动CFAR参数常常是白忙活。具体顺序是先跑出RD图确认目标峰值位置和理论计算一致再固定Pfa调参考窗和保护窗观察检测结果的虚警点和漏检点最后如果检测结果依然有干扰回去检查是否需要在多普勒维做杂波对消或者调整窗函数。6.2 一套建议的验证流程每改一组参数建议存三样中间结果RD图功率谱、CFAR阈值图、检测二值图。这三样放在一起才能判断问题是出在目标本身太弱还是阈值选得不对。我平时会针对同一组目标分别用 Pfa1e-4、1e-5、1e-6 跑三次看检测点的变化趋势。如果Pfa调低一个数量级虚警就消失说明阈值方向选对了如果Pfa调到1e-6仍然虚警满满那基本是参考窗或者保护窗出了问题。6.3 一点个人体会做完整套仿真之后最有价值的收获其实是把每个环节的物理含义都过了一遍波形参数决定分辨率FFT决定频谱结构CFAR决定检测门限四者咬合紧密。跳过CFAR直接取峰值的做法在单目标仿真里能出漂亮结果但一换到多目标和噪声边缘就会露馅。我的习惯是仿真里也尽量模拟噪声不均匀的场景故意让一侧底噪比另一侧高几个dB这样CFAR的价值立刻体现出来。最后分享一个小技巧可以把CA-CFAR的阈值矩阵直接从函数里返回出来画完RD图之后顺手画一张阈值剖面图横轴是距离bin纵轴是阈值幅度。你会发现阈值曲线并不是平的而是跟随噪声底起伏这一眼就能看出CFAR自适应的效果。看懂了这张曲线也就看懂了恒虚警检测。