简介：面向雷达信号处理学习者与研究者的MATLAB恒虚警检测（CFAR）仿真程序，聚焦在噪声背景下如何设置恒定虚警门限以准确识别目标。代码精简，压缩包仅含1个m文件、大小约2KB，可直接运行或修改参数，适合初学雷达检测算法时对照原理逐行理解。资源已有464人学习下载。程序覆盖常见CFAR实现思路，如单元平均CA-CFAR、顺序统计OS-CFAR等，通过估计背景噪声功率并自适应调整门限，帮助读者清晰掌握保护单元与训练单元的选取方法，并能将频谱图中的目标检测流程快速落地为仿真实验。通过设置不同参考窗口长度可以观察门限变化对检测结果的影响，便于动手探索。作为练习素材，既能验证检测概率与虚警率的关系，也可在此基础上扩展至目标定位、多目标处理等更高级雷达信号处理场景。
1. 雷达恒虚警检测的核心矛盾：门限要跟着背景噪声滑动，而不是定死
雷达恒虚警检测（CFAR，Constant False Alarm Rate，对应标题里的 radar 恒虚警）不是把检测门限抬高，而是让门限跟着被检测单元周围的背景功率实时滑动。海杂波、地物杂波、降雨杂波的功率在距离维和帧间能差十几到几十 dB，固定门限要么埋掉目标，要么让虚警指数级上涨。CFAR 用 CUT 两侧的参考单元估计背景功率，再按锁定的虚警概率乘系数得门限，因此背景变化时单 bin 虚警率仍被压在同一个数量级。下文推导 CA-CFAR 阈值系数 α 的闭式解，给出 MATLAB 手写滑窗与 phased.CFARDetector 的对照，最后落到参数设置与蒙特卡洛验证。适合做雷达检测仿真、或拿到 RD 谱不知道怎么配 CFAR 窗口的工程师。
2. 恒虚警检测的判决理论：固定门限为何失效，CA-CFAR 的 α 怎么推
2.1 平方律检波下的固定门限：噪声功率一变，虚警率立刻失控
回波经过匹配滤波和正交检波后，I/Q 两路送给平方律检波器 z = I² + Q²。无目标时 I、Q 是零均值高斯，z 服从均值为 σ² 的指数分布，σ² 就是该距离单元内的噪声加杂波平均功率。如果用一个固定门限 T 判决，单 bin 虚警概率是 Pfa = exp(−T/σ²)。
这公式说明两件事：门限和 σ² 完全耦合；σ² 若翻一倍，Pfa 从 1e-6 变成 exp(−T/2σ²) = sqrt(1e-6) ≈ 1e-3，三个数量级没了。真实场景里海况变化、降雨区域、天线波束扫过地物，σ² 在同一个扫描帧内就可以差 20 dB 以上，固定门限根本没法选。所以恒虚警检测的本质是放弃"噪声功率已知"的假设，在判决的同时从数据里估计 σ²。估计窗口就是被检单元（CUT, cell under test）两侧的参考单元，滑窗每移动一个距离单元就重新估计一次，这就是雷达恒虚警在工程和书本上的标准模型。
2.2 CA-CFAR 的闭式解：Pfa=(1+α)^(−N)，α=Pfa^(−1/N)−1
单元平均 CFAR（CA-CFAR）把左右两侧共 N 个参考单元的平方律输出求和 Z = Σzᵢ。无目标时 zᵢ 独立同分布指数，Z 服从形状参数 N、尺度 σ² 的 Gamma 分布。给定 Z，目标单元判决 z₀ > αZ 的条件虚警概率是 exp(−αZ/σ²)，对 Z 求期望，用 Gamma 分布的矩母函数 (1−t)^(−N)，得到：
Pfa = E[exp(−αZ/σ²)] = (1+α)^(−N)
反解就是阈值系数 α = Pfa^(−1/N) − 1，注意这是对"求和量 Z"的乘子。如果习惯用平均形式 Z̄ = Z/N，等价的门限乘子是 N(Pfa^(−1/N) − 1)，工程书和 MATLAB 文档里两种写法都有，换算时别搞混，这是后面手写代码最容易出错的地方。
数值上，N=32、Pfa=1e-4 时 α = (1e4)^(1/32) − 1 ≈ 0.3335，门限约等于 0.3335×Z，等价于 10.67 倍平均噪声功率。理想情况下若 σ² 精确已知，同 Pfa 的门限是 −ln(1e-4)·σ² = 9.21σ²。多花的 10·log10(10.67/9.21) ≈ 0.64 dB 就是 CA-CFAR 在均匀背景里的 CFAR 检测损失，这个损失随 N 增大而减小，第 4 章给出不同 N 的对照表。
提示：CFAR 损失牺牲的是检测概率（Pd），不是虚警概率；Pfa 被设计值锁死，付出的代价是灵敏度随参考单元数变差。
推导过程还有一个容易被忽略的结论：Pfa 只依赖 α 和 N，不依赖 σ²，这正是"恒虚警"三个字的数学来源。
2.3 GO、SO、OS-CFAR：为多目标和杂波边缘准备的三种改型
CA-CFAR 在均匀背景里最优，但它假设所有参考单元都来自同一分布。目标进入参考窗（多目标场景）或者窗内一半是杂波一半是噪声（边缘场景）时，CA 的估计就错了。工程上按"哪种场景最致命"选改型：
GO-CFAR（选大）把左右两半窗分别求和，取大者当背景估计。杂波边缘处 CUT 在杂波区内、窗外是噪声时，CA 会把两侧平均导致门限偏低、虚警暴增，GO 取大窗能压住这类虚警；代价是两个邻近目标更容易被一起遮蔽，因为目标能量被算进大窗。
SO-CFAR（选小）恰好相反，两窗取小者，专门用于分辨相距很近的目标：强目标占住一侧参考窗时，CA 和 GO 的门限都被抬高，SO 可以选另一侧干净窗恢复门限。代价是杂波边缘处、CUT 位于杂波一侧时，SO 会选到噪声侧的小值，虚警明显上升。
OS-CFAR（有序统计）把全部参考单元排序，取第 k 个值作为背景估计。只要参考窗内的干扰目标数不超过 N−k，门限就基本不动，抗多目标能力最强；k 一般取 0.6N 到 0.75N，k 越大抗干扰目标越多、均匀背景损失也越大。常见取 k=3N/4，损失约 0.3~1 dB。实际工程里，海上目标稀疏时用 OS 的多，近程地物边缘多、目标密集时 SO/GO 配合使用。
| CFAR 类型 | 背景估计 | 均匀背景损失 | 邻近双目标 | 杂波边缘 |
|---|---|---|---|---|
| CA | 两窗求和 | 基准，N=32 约 0.64 dB | 弱目标被遮蔽 | 边缘内侧虚警偏高 |
| GO | 两窗选大 | 略增 0.1~0.3 dB | 遮蔽范围更大 | 抑制边缘虚警 |
| SO | 两窗选小 | 略增 0.2~0.4 dB | 两目标均可检出 | 边缘外侧虚警明显 |
| OS | 排序取第 k 个 | 与 k 相关，k=3N/4 约 0.3~1 dB | 可容忍 ≤N−k 个干扰 | 中等 |
MATLAB 的 phased.CFARDetector 只内置 CA/GO/SO，OS 需要自己写，第 4 章的多目标验证里会给出 SO 的手写版本。
3. 用 MATLAB 实现雷达恒虚警检测：手写滑窗与 phased.CFARDetector 对照
3.1 先造一份带目标的平方律检波数据
仿真从复高斯噪声开始。单单元噪声平均功率 σ² 时，实部虚部各占 σ²/2，平方律输出：
rng(202401); L = 4096; % 距离单元数，按一次脉冲的采样点算 sigma2 = 1; % 单单元噪声平均功率 x = sigma2 * (abs(randn(1, L) + 1i*randn(1, L)).^2) / 2; % randn 方差为 1，|I+jQ|^2 均值 2，除 2 后均值恰为 sigma2 A = 7; % 目标幅度 kTarget = 2500; x(kTarget) = x(kTarget) + A^2; % 功率域直接叠加目标 fprintf('目标 SNR = %.1f dB\n', 10*log10(A^2/sigma2));这段数据的物理含义是：I/Q 两路噪声方差各 0.5，合成后单 bin 指数分布均值为 1；目标功率 A² = 49 对应 16.9 dB SNR，在 CA-CFAR、Pfa=1e-4、N=32 的配置下检测概率很高，第 2500 号单元应当被稳定检出，方便后面肉眼确认结果。把目标直接加在功率域，等效于假设目标回波与噪声不相关，这是点目标仿真里的常用做法。
3.2 手写 CA-CFAR：循环版理解流程，向量化版应对 RD 谱
先写一个直白的循环实现，每个 CUT 单独取参考窗：
function det = myCACFAR(x, Nc, Ng, pfa) % x: 平方律检波输出（线性功率域），一行一个距离单元 % Nc: 单侧参考单元数，Ng: 单侧保护单元数 N = 2 * Nc; alpha = pfa^(-1/N) - 1; % 对求和量的阈值系数 L = numel(x); det = false(1, L); for k = Nc+Ng+1 : L-Nc-Ng refs = [x(k-Nc-Ng:k-Ng-1), x(k+Ng+1:k+Ng+Nc)]; Z = sum(refs); det(k) = x(k) > alpha * Z; end end循环版的好处是边界条件一目了然：k 从 Nc+Ng+1 起，保证左侧参考窗不越界；右侧同理，所以开头和结尾各有 Nc+Ng 个单元默认不判决。α 用的是求和量乘子 pfa^(−1/N)−1，不是平均形式的 N(...−1)，两种形式在 2.2 节已经说明，这里直接按代码能对上的来。
数据量大（比如 128×512 的 RD 谱逐距离门处理）时循环版太慢，换成 movsum 的向量化写法：
function det = cacfarVec(x, Nc, Ng, pfa) alpha = pfa^(-1/(2*Nc)) - 1; winSum = movsum(x, 2*(Nc+Ng)+1, 'Endpoints', 'discard'); guardSum = movsum(x, 2*Ng+1, 'Endpoints', 'discard'); refSum = winSum - guardSum - x; % 参考窗 = 全窗 - 保护带 - CUT det = x > alpha .* refSum; det(1:Nc+Ng) = false; % 边界单元不参与判决 det(end-Nc-Ng+1:end) = false; endmovsum 以每个单元为中心滑窗求和，全窗长度 2(Nc+Ng)+1 减去保护带长度 2Ng+1 再减 CUT 本身，正好剩下两侧共 2Nc 个参考单元，不用写索引循环。注意 'Endpoints','discard' 会让窗不满的位置补 0，如果不去掉边界，refSum 在边界处变成负数，x > α·refSum 恒真，就会在距离维两端刷出一排假目标——这个坑在实际代码里见过不止一次。
3.3 工具箱对照：phased.CFARDetector 的 NumTrainingCells 是总数
Phased Array System Toolbox 的 CFAR 检测器参数名和手写习惯不同，NumTrainingCells、NumGuardCells 都是"两侧总和"，而手写函数里的 Nc、Ng 是单侧值，映射时要乘 2：
Nc = 16; Ng = 2; hc = phased.CFARDetector( ... 'NumTrainingCells', 2*Nc, ... % 两侧参考单元总数 'NumGuardCells', 2*Ng, ... % 两侧保护单元总数 'ProbabilityFalseAlarm', 1e-4, ... 'ThresholdFactor', 'Auto', ... 'OutputFormat', 'Threshold'); % 先输出阈值，便于和手写版对账 idx = (1:L).'; thTool = hc(x(:), idx); % 手写版同样位置算一遍阈值 thManual = zeros(L, 1); for k = Nc+Ng+1 : L-Nc-Ng refs = [x(k-Nc-Ng:k-Ng-1), x(k+Ng+1:k+Ng+Nc)]; thManual(k) = (1e-4^(-1/(2*Nc)) - 1) * sum(refs); end disp(max(abs(thTool(Nc+Ng+1:L-Nc-Ng) - ... thManual(Nc+Ng+1:L-Nc-Ng)))); % 理想情况接近 1e-14ThresholdFactor 设为 'Auto' 时，对象内部按闭式公式由 Pfa 换算系数，等价于手写版里 alpha 乘在求和量上的做法，因此两者在远离边界的单元上数值一致。工具箱版本对边界单元有自己的处理策略，不一定和手写版相同，工程上直接约定"只统计有效区间"更干净，不要指望两边边界行为一致。要输出判决结果，把 OutputFormat 换成 'Detection result' 再调一次：
hc.OutputFormat = 'Detection result'; detTool = hc(x(:), idx).'; detManual = myCACFAR(x, Nc, Ng, 1e-4); fprintf('手写版检出目标: %d\n', detManual(kTarget)); fprintf('工具箱检出目标: %d\n', detTool(kTarget));| 手写版参数 | 工具箱属性 | 换算关系 |
|---|---|---|
| 单侧参考 Nc | NumTrainingCells | 2×Nc |
| 单侧保护 Ng | NumGuardCells | 2×Ng |
| pfa | ProbabilityFalseAlarm | 不变 |
| α（求和形式） | ThresholdFactor='Auto' | 由 Pfa 闭式计算 |
4. 恒虚警检测参数怎么设：保护单元、参考单元、Pfa 的联动与典型坑
4.1 保护单元数量由目标张角和脉压副瓣决定
保护单元的作用是把目标自身能量挡在参考窗之外。点目标经过匹配滤波/脉压后不是单根谱线，主瓣会占 2~3 个距离单元，脉压副瓣还能再外扩几个单元；如果目标能量漏进参考窗，Z 被抬高，门限跟着涨，第一个受害的就是目标自己。距离维上 Ng 取 2~4 比较常见，线性调频大时宽带宽积信号副瓣较高、或目标在距离维有明显展宽时取到 6~8。二维 RD 谱上，目标在速度维也有展宽，GuardBandSize 的多普勒分量必须一起给，只给距离维保护带、多普勒维不加，等于让旁瓣和多普勒泄漏直接污染参考单元。
4.2 参考单元数量：估计方差、CFAR 损失、均匀性三者的折中
参考单元 N 决定背景估计的精度。Z/σ² 服从 Gamma(N,1)，相对标准差 1/√N：N=16 时估计抖动 25%，N=64 时 12.5%。估计越抖，门限越不稳定，"恒虚警"就越差。另一方面 N 增大让阈值系数下降，CFAR 损失变小，但它同时拉长滑窗，窗外非均匀场景（杂波边缘、目标群）更容易毁掉估计。Pfa=1e-4 背景下不同 N 的对照：
| 参考单元总数 N | 均值形式阈值系数 | 相对理想门限的 CFAR 损失 |
|---|---|---|
| 16 | 12.4 | ≈1.3 dB |
| 32 | 10.7 | ≈0.64 dB |
| 64 | 9.9 | ≈0.34 dB |
| 128 | 9.5 | ≈0.16 dB |
损失用 10·log10(α_mean / (−ln Pfa)) 估算，α_mean = N(Pfa^(−1/N)−1)，理想门限 −ln(Pfa)·σ² = 9.21σ²。实际工程里取总数 32~64、即单侧 16~32 的最多：再多，均匀背景只省零点几 dB，非均匀代价却线性变大。Pfa 一般取 1e-4~1e-6，毫米波雷达量产点云里 Pfa 取得太松，一帧会出现几千个虚警点，后处理扛不住。
4.3 坑一：在 dB 域做单元平均，门限系统性偏低
把平方律输出转 dB 后再算参考单元均值，是入门最容易犯的错。指数分布 ln 的期望比 ln(均值) 小欧拉常数 γ≈0.577，换算成功率域相当于低估约 2.5 dB，且 N 越大这个偏差越稳定，不会随平均而消失。门限按设计 α 乘在这个偏低的估计上，实际虚警率会比名义 Pfa 高一到两个数量级。
CFAR 的闭式公式建立在平方律功率域，喂给检测器的数据也必须是功率（或幅度平方），不能是 dB、不能是幅度绝对值。如果你想基于幅度实现，闭式的 α 表达式要重新推导，不能直接套 pfa^(−1/N)−1。常见做法是全程线性功率域，只在显示和画图时转 dB。
提示：dB 域平均低估约 2.5 dB 的偏差不会随 N 增大而消失，这是系统性偏差，不是统计抖动。
4.4 坑二到坑四：多目标遮蔽、杂波边缘、边界裁剪
多目标遮蔽是 CA-CFAR 最经典的失效模式：两个目标相距小于参考窗半宽时，强目标进入弱目标的参考窗，Z 被抬高导致弱目标漏检。SO-CFAR 就是为这个设计的，手写只改一处：
function det = mySOCFAR(x, Nc, Ng, pfa) alpha = pfa^(-1/(2*Nc)) - 1; L = numel(x); det = false(1, L); for k = Nc+Ng+1 : L-Nc-Ng lead = x(k-Nc-Ng : k-Ng-1); trail = x(k+Ng+1 : k+Ng+Nc); Z = min(sum(lead), sum(trail)); % 取两窗小者 det(k) = x(k) > alpha * Z; end end验证：强目标放在 2009 号单元，弱目标放在 2025 号单元（间隔 16 < Nc+Ng=18），弱目标左侧参考窗会包含强目标：
L = 4096; x2 = (abs(randn(1,L) + 1i*randn(1,L))).^2 / 2; x2(2009) = x2(2009) + 12^2; % 强目标，功率 144 x2(2025) = x2(2025) + 4^2; % 弱目标，SNR 12 dB detCA = myCACFAR(x2, 16, 2, 1e-4); detSO = mySOCFAR(x2, 16, 2, 1e-4); fprintf('CA 检出弱目标: %d\n', detCA(2025)); % 0 fprintf('SO 检出弱目标: %d\n', detSO(2025)); % 1CA 会把 4²+12² 的功率算进背景，门限抬高后 12 dB 的弱目标直接漏掉；SO 选右侧干净窗，弱目标恢复检出。注意强目标要放到弱目标的左侧窗内而不是保护带内，间隔必须大于 Ng 但小于 Nc+Ng，参数不同结论会变，复现时先按这组数跑通再改。
杂波边缘的典型场景是云雨区或地物边界：窗一半在杂波内一半在外。处理顺序一般是先做一次 GO-CFAR 保住边缘内侧不虚警，再在均匀区段用 CA 或 SO 细分；或者用 OS-CFAR 一勺烩，但 k 值要对参考窗里可能出现的干扰目标数留余量。最后是边界裁剪：任何实现都要显式屏蔽前后 Nc+Ng 个单元，工具箱和手写版边界行为不一致时，以"统计时剔除边界"为准，否则蒙特卡洛统计虚警率时会被边界假目标污染。
5. 验证与进阶：蒙特卡洛核对虚警率，二维 CFAR 直接上距离-多普勒图
5.1 蒙特卡洛要跑多少样本，才敢说恒虚警达标
CFAR 写完后第一件事是验证实际虚警率贴近设计值。做法：生成纯噪声序列（无目标），跑检测，统计"虚警数 / 有效单元总数"：
rng(7); M = 300; L = 8192; Nc = 16; Ng = 2; pfa = 1e-4; valid = L - 2*(Nc+Ng); hit = 0; for m = 1:M xn = (abs(randn(1,L) + 1i*randn(1,L))).^2 / 2; hit = hit + nnz(myCACFAR(xn, Nc, Ng, pfa)); end pfaSim = hit / (M * valid); fprintf('模拟 Pfa = %.3g, 设计值 = %.3g\n', pfaSim, pfa);总共 300×8156≈2.45e6 个独立判决，按 pfa=1e-4 期望约 245 次虚警，标准差约 15.6，因此 pfaSim 相对偏差在 6% 上下是正常的，别看到 1.08e-4 就以为代码有问题。验证更低的 Pfa 时先算样本量：要让 95% 概率至少观察到一个虚警，需要的独立样本数约为 ln(0.05)/ln(1−Pfa)。
| 设计 Pfa | 至少出现 1 次虚警所需样本（95%） | 工程建议样本量 |
|---|---|---|
| 1e-4 | 约 3.0e4 | 1e6 以上 |
| 1e-5 | 约 3.0e5 | 1e7 以上 |
| 1e-6 | 约 3.0e6 | 1e8 以上 |
只跑 1e5 个样本就去判断 1e-5 的虚警率，即使 0 虚警也什么都说明不了，这是验证时最常见的统计功率不足问题。
5.2 二维 CFAR：把滑窗从一维距离换成距离-多普勒平面
毫米波雷达目标检测里普遍的做法是在距离-多普勒（RD）谱上做二维 CFAR，速度维和距离维同时估计背景，避免多普勒维上的强杂波干扰。工具箱对应对象是 phased.CFARDetector2D，TrainingBandSize 和 GuardBandSize 都是 [行, 列] 形式，数值含义是每个维度以 CUT 为中心的总窗宽与内层保护带：
cd2 = phased.CFARDetector2D( ... 'Method', 'CFAR', ... 'GuardBandSize', [2 2], ... % 距离维、多普勒维各 2 个保护单元 'TrainingBandSize', [8 8], ... % 外框总尺寸，已包含保护带 'ProbabilityFalseAlarm', 1e-5, ... 'ThresholdFactor', 'Auto', ... 'OutputFormat', 'Threshold'); [rr, dd] = ndgrid(1:size(rdMap,1), 1:size(rdMap,2)); cutIdx = [rr(:), dd(:)]; thMap = cd2(rdMap, cutIdx); % 列向量，reshape 成图便于画 thMap = reshape(thMap, size(rdMap));rdMap 必须是线性功率域数据，不能是 log 幅度图；FFT 窗函数造成的多普勒泄漏会推高邻近单元，所以多普勒维保护带不能省。整帧 RD 谱全图跑二维 CFAR 在实时循环里耗时可观，常见做法是先用较低 Pfa 或幅度门限做预筛选，只对候选单元计算二维 CFAR 阈值，把每帧的 CUT 数量从几万压到几千。边界行为同样要注意：TrainingBandSize 超出 RD 图边界的单元按规定应丢弃，别把 FFT 补零区域当成可判决单元。验证一维 CFAR 的蒙特卡洛流程加一层外循环即可套到二维，只是有效单元数要按"去除行、列边界后的面积"重新算。交付前把 α 的数值打印出来和手写版对一遍，确认 Pfa、参考单元数、阈值系数三者没有一处被"均值形式/求和形式"的表述带偏。
本文还有配套的精品资源，点击获取
