☰
FMCW雷达距离多普勒仿真与CA-CFAR目标检测完整实现
2026/10/22 5:14:37 网站建设 项目流程

这几年做雷达信号处理,被问得最多的不是波束赋形这类进阶方向,反而是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 MHz,chirp时长 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_chirp=128时,Δ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 m
chirp时长T_c20 μs决定最大速度约48.7 m/s
ADC采样率F_s20 MHz决定最大距离150 m
chirp数/帧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(:, n+1) = ifMat(:, n+1) + ... 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的全称是单元平均恒虚警率。每个待检测单元称作CUT(Cell Under Test)。在CUT周围取一个矩形邻域,紧挨着CUT的一圈是保护单元(guard cells),保护单元外面一圈是参考单元(reference cells)。保护单元的存在是为了防止目标主瓣或旁瓣能量污染噪声估计。一个目标在距离维和速度维都可能占据多个bin,如果这些bin落入参考单元,会把噪声均值拉高,阈值抬升,最终导致目标自己被"吃掉"。

二维CA-CFAR的窗口结构可以这样理解:设参考单元半径为 Rr,保护单元半径为 Gr,那么整个窗口大小是 (2Rr+1)×(2Rr+1),中心区域 (2Gr+1)×(2Gr+1) 不参与均值估计。用我前面那组参数计算,当 Rr=6、Gr=3时,窗口是13×13,保护区域7×7,参考单元总数:

Nref = (2*6+1)^2 - (2*3+1)^2 = 169 - 49 = 120

阈值系数α和虚警概率Pfa的关系为:

Pfa = (1 + alpha/Nref)^(-Nref)

反过来求α:

alpha = Nref * (Pfa^(-1/Nref) - 1)

取 Pfa = 1e-6,Nref = 120:

alpha = 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*Rr+1, 2*Rr+1); win(Rr+1-Gr:Rr+1+Gr, Rr+1-Gr:Rr+1+Gr) = 0; win(Rr+1, Rr+1) = 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(Rr+1:end-Rr, Rr+1: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阈值图、检测二值图。这三样放在一起才能判断问题是出在"目标本身太弱"还是"阈值选得不对"。我平时会针对同一组目标分别用 Pfa=1e-4、1e-5、1e-6 跑三次,看检测点的变化趋势。如果Pfa调低一个数量级虚警就消失,说明阈值方向选对了;如果Pfa调到1e-6仍然虚警满满,那基本是参考窗或者保护窗出了问题。

6.3 一点个人体会

做完整套仿真之后,最有价值的收获其实是"把每个环节的物理含义都过了一遍":波形参数决定分辨率,FFT决定频谱结构,CFAR决定检测门限,四者咬合紧密。跳过CFAR直接取峰值的做法,在单目标仿真里能出漂亮结果,但一换到多目标和噪声边缘就会露馅。我的习惯是仿真里也尽量模拟"噪声不均匀"的场景,故意让一侧底噪比另一侧高几个dB,这样CFAR的价值立刻体现出来。

最后分享一个小技巧:可以把CA-CFAR的阈值矩阵直接从函数里返回出来,画完RD图之后顺手画一张阈值剖面图,横轴是距离bin,纵轴是阈值幅度。你会发现阈值曲线并不是平的,而是跟随噪声底起伏,这一眼就能看出CFAR自适应的效果。看懂了这张曲线,也就看懂了恒虚警检测。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询