简介资源包是一份面向雷达信号处理学习者与工程师的MATLAB仿真源码聚焦调频连续波FMCW雷达回波信号中的杂波抑制问题也可推广到其他频域杂波抑制场景。该仿真基于CLEAN算法在频域中对真实雷达回波数据进行处理通过相减函数抑制杂波分量从而还原目标信号相减函数当前采用余弦形式若信号为复数域可将其改成指数形式便于适配不同场景。整包仅含1个文件即MATLAB脚本源码体积约1KB代码短小精悍、无额外依赖复制到MATLAB环境中即可运行无需复杂配置很适合初学者对照学习CLEAN算法流程也可通过修改参数快速应用于自己的雷达回波数据。目前已有897人学习/下载说明其在雷达杂波抑制相关任务中具有一定参考价值。对正在做FMCW雷达仿真或频域杂波抑制项目的读者这份源码提供了一个可直接运行的算法模板既能帮助理解CLEAN算法的频域实现细节又能作为自定义开发的基础脚本尤其是采用真实回波数据时可作为快速验证的第一步。1. main()把雷达杂波仿真和CLEAN算法串起来的主心骨很多拿到雷达回波数据的人最先打开的是一个画图脚本而不是那个把数据读取、杂波仿真和CLEAN算法一步步串起来的main()。等你在距离-多普勒图上真正看到那一片杂波时才会意识到杂波不只是看起来“脏”它还会把弱目标死死按在强散射点的旁瓣里常规的恒虚警检测拿它没什么办法。这个工程方案的核心就是用main()搭一条完整链路——先在可控条件下仿真雷达杂波和回波数据再用CLEAN算法把杂波与旁瓣一层层剥掉最终输出目标列表。适合正在做雷达算法落地验证的工程师也适合想把仿真原型往Python/C工程里迁移的信号处理从业者。2. 先把杂波做出来从瑞利到K分布的雷达杂波仿真2.1 为什么仿真先于CLEAN你需要一块可控的“脏”回波直接拿实测数据调CLEAN最大的问题是“黑匣子”你不知道某个峰值到底是目标、杂波还是噪声更不知道算法是不是真的把杂波剥干净了。仿真杂波的价值在于目标位置、幅度、多普勒频移、杂波分布全部由自己指定CLEAN跑完一对比就知道对不对、差多少。后面要换场景只需要改几个参数不用扛着雷达外场跑一遍。杂波模型的选择一般看应用场景。均匀地面且雷达分辨率较低时大量独立散射点的叠加会让幅度趋近瑞利分布这也是教科书默认的模型。高分辨率雷达看海面地杂波的拖尾特别重强散射点密集这时候瑞利已经兜不住要用K分布或韦布尔分布。K分布的形状参数越小拖尾越重就越接近“少量强散射体加上大量微弱散射体”的海面实况。地杂波、海杂波在慢时间维上还有自己的多普勒谱通常建模成高斯谱中心频移对应杂波散射体的平均径向速度谱宽对应海浪起伏或天线扫描展宽。所以一套能用的杂波仿真至少要有三层可调幅度分布、多普勒谱形状、信杂比和信噪比。三层都控住生成的回波才能当作CLEAN算法的“测试靶场”。否则参数随机一拉出来的东西自己都不知道是什么那就是靠玄学调算法了。2.2 在main()里生成一帧含目标与杂波的IQ数据Python脚本与参数下面这段代码生成的是脉压后的距离-多普勒复图像横轴距离门纵轴慢时间脉冲每个元素是复数IQ。目标直接以“距离门多普勒频移信噪比”的形式插入杂波按CNR信杂比控制强度。import numpy as np def generate_frame(n_range1024, n_pulse128, prf1000.0, cnr_db30.0, clutter_typerayleigh, target_listNone): 生成一帧距离-多普勒域复数据维度 [n_pulse, n_range]。 target_list: [(range_bin, doppler_hz, snr_db), ...] 假定噪声功率归一为1。 rng np.random.default_rng(20240601) # 接收机噪声复高斯实部虚部方差各0.5总功率为1 noise (rng.standard_normal((n_pulse, n_range)) 1j * rng.standard_normal((n_pulse, n_range))) / np.sqrt(2) # 杂波基底复高斯得到瑞利幅度分布 clutter (rng.standard_normal((n_pulse, n_range)) 1j * rng.standard_normal((n_pulse, n_range))) / np.sqrt(2) # 慢时间维多普勒谱成形高斯谱中心-8Hz谱宽5Hz freq_axis np.fft.fftfreq(n_pulse, 1.0 / prf) sigma_f 5.0 center_f -8.0 filter_h np.exp(-0.5 * ((freq_axis - center_f) / sigma_f) ** 2) clutter np.fft.ifft(filter_h[:, None] * np.fft.fft(clutter, axis0), axis0) # 按信杂比归一化 noise_pow np.mean(np.abs(noise) ** 2) clutter_pow np.mean(np.abs(clutter) ** 2) clutter clutter * np.sqrt(10 ** (cnr_db / 10) * noise_pow / clutter_pow) frame noise clutter # 目标噪声功率为1幅度 10^(snr_db/20) if target_list: for rbin, df_hz, snr_db in target_list: dop_bin int(round(df_hz / prf * n_pulse)) amp 10 ** (snr_db / 20.0) frame[dop_bin % n_pulse, rbin] amp return frame这段代码的逻辑是先铺噪声底再叠加已经做过谱成形的杂波最后逐个放入目标。多普勒谱成形放在杂波之后、功率归一化之前是为了避免滤波改变杂波功率后还要二次校准。n_pulse和n_range决定帧尺寸prf决定多普勒轴刻度。target_list里的doppler_hz会被换成脉冲索引负数代表目标朝雷达运动正数远离。这个生成函数已经可以直接在main()里调用if __name__ __main__: # 两个距离接近的目标一个强一个弱用于验证CLEAN的旁瓣剥离能力 targets [(400, 0.0, 35.0), (403, 0.0, 20.0)] frame generate_frame(target_listtargets) # frame[0, :] 就是第一个脉冲对应的距离像两个目标距离只差3个门强目标脉压后的距离旁瓣正好压在弱目标的位置上这是后面CLEAN算法最典型的用武之地。3. CLEAN算法原理从射电天文的“脏图”到雷达的“脏回波”3.1 核心假设与迭代逻辑CLEAN算法最早来自射电天文Högbom在1974年提出用它重建综合孔径图像。观测到的“脏图”被认为是真实天空亮度分布与天线方向图点扩散函数PSF的卷积再加上噪声。CLEAN的做法是在脏图里找最亮的点认为它是真实源把这个亮点的PSF分量按一定比例减掉记录该源的幅度和位置然后再找下一个最亮点重复迭代直到剩余信号低到门限以下。这套逻辑放到雷达里完全成立。脉冲压缩后一个理想点目标在距离维上是一个sinc形状的响应主瓣旁边挂着一串距离旁瓣。如果同一个距离门附近有个强目标它的旁瓣会把弱目标盖住。雷达回波里的目标、强杂波散射体都可以看成稀疏分布的“点源”而CLEAN正是把这些点源的旁瓣逐个剥离出来的迭代算法。需要明确一点CLEAN不是超分辨算法。两个目标只要落进同一个主瓣宽度内峰值搜索就分不开它们CLEAN也一样无能为力。它的价值在于把“主瓣宽度之外被旁瓣掩盖”的目标找回来而不是突破瑞利分辨极限。理解了这一点你就不会拿它去硬解主瓣内多目标。3.2 增益、PSF与门限三个决定收敛行为的数字第一个数是循环增益loop gain记作g。标准CLEAN每次只减掉峰值幅度的g倍而不是全额减掉。g太小迭代次数多收敛慢g太大PSF与真实响应只要有一点不匹配误差就会被放大形成目标两边的虚假亮点。我一般取0.2~0.4。如果PSF是从标定数据测来的、质量很好可以放宽到0.5如果PSF只是理论近似就老实取0.25以下。第二个数是PSF本身。在雷达里它就是脉冲压缩后点目标响应的归一化复序列峰值恒为1。PSF必须覆盖到主瓣以外至少4~6个旁瓣否则每次减除都残留一段尾巴多次迭代后残留尾巴会累积成梳状结构。PSF的获取有两种从发射信号的自相关函数推导或者用强反射体的实测脉压数据截取。实测数据比理论模型可靠得多。第三个数是停止门限。CLEAN迭代不能无限跑否则会把噪声也当目标剥出来。门限通常设成“噪声底以上10~15dB”噪声底用最远距离门无目标区的幅度中位数来估计。这个门限要动态计算绝不能写死成某个绝对幅度值数据换一段就变。3.3 一维CLEAN的迭代式与伪代码一维CLEAN的迭代公式可以写成r^(k1)(t) r^(k)(t) - g · A_k · h(t - t_k)其中r为残差信号t_k为第k次迭代找到的峰值位置A_k是峰值位置的复幅度h是PSF。每次迭代把PSF平移到t_k处乘上g和A_k再从残差中减掉同时在目标列表里记录一个幅度为g·A_k的点目标。伪代码就是下面这个结构resid data.copy() peaks [] while True: idx np.argmax(np.abs(resid)) amp resid[idx] if np.abs(amp) threshold: break peaks.append((idx, amp)) # 以idx为中心将PSF缩放后减去 resid resid - gain * amp * shift_psf(psf, idx)注意amp是复数shift_psf把PSF数组平移到idx位置。复数运算必须全程保持不能只对幅度做操作否则相位信息丢得一干二净。4. 在main()里实现CLEAN从数据读取到目标列表输出4.1 main()的职责参数、文件读取与按帧循环main()在雷达处理链路里不是“随便写个函数入口”那么简单。它至少要负责三件事解析参数、按帧读取数据、把处理完的目标列表写出去。很多工程里CLEAN算法本身没问题翻车翻在数据路径写死、参数没传对、或者读到的IQ数据格式和算法预期不一致。main函数参数这块C/C工程里标准写法是int main(int argc, char *argv[])argc是参数个数argv是参数内容。我习惯把数据文件路径、PSF文件路径、迭代增益、停止门限全部通过参数传入而不是硬编码在源码里。下面这段是典型的框架import sys def main(): # argv[0]是程序名真正参数从argv[1]开始 if len(sys.argv) 3: print(usage: radar_clean iq_data psf.npy [gain] [thresh_db]) return 1 data read_iq(sys.argv[1]) psf np.load(sys.argv[2]) gain float(sys.argv[3]) if len(sys.argv) 3 else 0.3 thresh_db float(sys.argv[4]) if len(sys.argv) 4 else 12.0 hits, resid clean_1d(data, psf, gaingain, thresh_dbthresh_db) for idx, amp in hits: print(f{idx}\t{20 * np.log10(abs(amp)):.1f} dB) return 0 if __name__ __main__: sys.exit(main())read_iq这里没有展开常见做法是读int16交错IQ先拆实部虚部再转复数。关键点是argv处理一定要先判断长度否则用户忘了传参数时程序会直接崩这在自动化批处理里非常讨厌。数据路径、PSF路径、算法参数全部走参数传换数据就不需要重新编译。4.2 完整的CLEAN最小实现clean_1d与PSF构造现在给出可以直接抄走的clean_1d实现。PSF用sinc近似长度取主瓣宽度4个门、左右各5个旁瓣宽度这样能覆盖到远端旁瓣。def make_psf(width4, n_side5): 生成脉压点目标响应的sinc近似峰值归一化为1。 width主瓣过零点间隔对应的距离门数 n_side主瓣每侧覆盖的旁瓣周期数 x np.arange(-n_side * width, n_side * width 1) psf np.sinc(x / width) return psf / np.max(np.abs(psf)) def clean_1d(data, psf, gain0.3, thresh_db12.0, max_iter512): 一维CLEAN输入距离像输出目标列表和残差。 data: 一维复数距离像 psf: 一维复数PSF峰值已归一化 返回 (hits, resid)hits元素为 (距离门, 复幅度) pad psf.size # 两侧各补一个PSF长度的零省去边界索引判断 resid np.zeros(data.size 2 * pad, dtypecomplex) resid[pad:pad data.size] data # 噪声底取末尾5%距离门用幅度中位数估计 tail data[int(data.size * 0.95):] noise_level np.median(np.abs(tail)) / 0.6745 threshold noise_level * 10 ** (thresh_db / 20.0) half psf.size // 2 hits [] for _ in range(max_iter): peak_idx int(np.argmax(np.abs(resid))) amp resid[peak_idx] if abs(amp) threshold: break # 记录目标注意只记增益乘过的幅度 hits.append((peak_idx - pad, amp * gain)) # 以峰值位置为中心减掉缩放后的PSF start peak_idx - half resid[start:start psf.size] - gain * amp * psf clean_resid resid[pad:pad data.size] return hits, clean_residclean_1d的实现把PSF两侧补齐后再做减除这样peak靠近首尾时不需要写一堆边界分支。工程上对长数据会改成切片裁剪避免每次迭代处理整段补零信号但教学和验证场景这个补零版本更不容易写错。PSF构造里的width参数要和实际波形匹配线性调频信号时宽带宽积越大主瓣越窄width就越小。如果是从实测标定数据取的PSF直接np.load进来替换make_psf即可。提示clean_1d的噪声底估计依赖末尾5%距离门没有目标。雷达采集时专门留一段空回程就是为了给这类处理提供参考区。4.3 把仿真、脉压和CLEAN串进main()里的完整示例下面再把第2章的仿真函数和第4章的CLEAN函数接在一起跑通最小闭环import numpy as np def main_sim(): # 生成一帧含杂波和两个目标的距离-多普勒数据 targets [(400, 0.0, 35.0), (403, 0.0, 18.0)] frame generate_frame(n_range1024, n_pulse128, cnr_db30.0, target_listtargets) # 这里只取第一个脉冲的距离像做一维CLEAN演示 # 真实工程中会对每个多普勒门分别做一维CLEAN或直接做二维CLEAN data frame[0, :] psf make_psf(width4, n_side5) hits, resid clean_1d(data, psf, gain0.25, thresh_db12.0) for idx, amp in sorted(hits): print(frange_bin{idx:5d} amp{abs(amp):8.2f} f{20 * np.log10(abs(amp)):6.1f} dB) if __name__ __main__: main_sim()运行后会看到强目标在400被检出403附近那个弱目标如果在原距离像里被旁瓣淹没经过CLEAN后会被单独挑出来。排序输出是为了方便和target_list里的真值对比。main_sim这个入口保持了main()风格的参数化设计场景和目标列表都通过变量传入换实验条件不用改内部逻辑。真实工程里距离-多普勒图要完整处理不能只取第一个脉冲。常见做法是沿多普勒维逐行做一维CLEAN或者把PSF扩展成二维距离维用sinc多普勒维用主瓣宽度为1的理想窗然后直接在二维矩阵上搜峰、减二维PSF。二维CLEAN的计算量大约是每迭代一次做一次二维卷积级减除数据帧大时要考虑性能。5. CLEAN算法避坑与常见问题5个翻车现场的修复清单5.1 参数参考迭代增益、门限、PSF长度的初始推荐先给一张参数速查表新场景跑不通时按这组初始值试。参数推荐初值调节方向loop_gain0.25残差出现振荡就减小收敛慢可适当加大到0.4thresh_db12 dB虚警多就加大漏检多就减小PSF长度主瓣宽度5~8倍残差里出现梳状旁瓣就加长max_iter目标数估计值的3~5倍到上限还没收敛检查门限和增益噪声底估计区末尾5%距离门该区有目标时换到盲区/空回程段这几个参数互相牵制。最常见的问题是门限设太低导致迭代停不下来max_iter被耗尽后输出几百个假目标或者门限设太高弱目标还没被剥出来就终止了。调试时先把gain固定0.25只调thresh_db一次只动一个变量不然出了问题根本不知道是谁的锅。5.2 两个与main()直接相关的编译期问题先说C/C场景里很常见的“编译器未包含main类型”报错。现象是编译时报错说main没有返回类型或者链接时提示找不到main入口。原因通常有三类一是把main写成了void main()这在很多C标准里直接判错二是main被放进了某个没参与链接的源文件三是拼写错误main拼成了mian。解决方法是统一写int main(int argc, char *argv[])并在工程文件树里确认入口文件确实参与编译。第二个是main函数参数处理不到位。雷达程序里argv通常要接收数据文件路径、PSF路径和算法参数。如果代码里直接用argv[1]而不检查argc用户少传一个参数程序就段错误。工程上一定要先判断argc数量给足默认值。这算老生常谈但我在实际项目里见过太多因为这问题导致的启动崩溃尤其在同一套代码跑不同批次数据的时候。5.3 五个翻车场景现象、原因、解决第一个翻车场景迭代结束后残差里还有明显的梳状旁瓣结构。原因是PSF截断太短只覆盖了主瓣附近远端旁瓣每次都没减干净累积起来形成周期性的残留。解决是把PSF长度扩大到主瓣宽度的5倍以上必要时从标定数据里取实测PSF。第二个翻车场景门限写死成一个固定幅度换一段数据后输出目标数量剧烈变化。原因是数据换了之后噪声底变了但门限没跟着变。解决是用末尾无目标区的中位数动态估计噪声底门限始终取“噪声底thresh_db”。第三个翻车场景gain设成0.6结果强目标两边出现对称的假目标。原因是增益太大PSF和真实响应只要偏差一点点减除后的残差里就会留下反向的伪峰下一轮又被当成真实目标剥出来。解决是gain降到0.2~0.3如果已经产生了对称假目标把它们合并掉再用低增益重新迭代。第四个翻车场景复数数据用了实数sinc做PSF目标恢复出来的相位乱跳、幅度偏低。原因是复数信号只减实数PSF实部虚部被同步缩放破坏了原始相位关系。解决是全程使用复数PSF和复数减除即使理论sinc是实函数也要保证相位数据流不被拆散。第五个翻车场景两个目标靠得很近弱目标直接被当作强目标的旁瓣一起减没了。原因是CLEAN只做单峰搜索主瓣内多目标无法分辨。解决是先做一次CFAR预检测对主瓣展宽区域改用“多目标最小二乘”模型同时估计两个复幅度而不是让CLEAN盲剥。6. 用合成数据验证CLEAN性能目标恢复误差与蒙特卡洛检查验证CLEAN做得对不对最直接的办法是拿第2章的仿真器生成已知真值的数据跑完CLEAN后对比输出和真值。位置误差不超过2个距离门、幅度误差不超过3dB视为命中非真值位置的检出都算虚警。每轮换随机种子重复100次统计检测概率和虚警率。下面是验证脚本的核心片段def run_trial(n_trials100): total_true 0 total_hit 0 total_fa 0 psf make_psf(width4, n_side5) for seed in range(n_trials): rng np.random.default_rng(seed) targets [(400, 0.0, 30.0), (406, 0.0, 20.0)] frame generate_frame(n_range1024, n_pulse128, cnr_db25.0, target_listtargets) data frame[0, :] hits, _ clean_1d(data, psf, gain0.25, thresh_db10.0) truth [400, 406] hit_pos [idx for idx, _ in hits] for t in truth: total_true 1 if any(abs(t - idx) 2 for idx in hit_pos): total_hit 1 detected sum(1 for idx in hit_pos if all(abs(t - idx) 2 for t in truth)) total_fa detected print(fdetection rate{total_hit / total_true:.2f}, ffalse alarms per frame{total_fa / n_trials:.2f})我自己的血泪经验是第一次跑这组验证时把thresh_db设到18检测率只有0.5以为算法有问题。后来逐帧看残差才发现强目标早就被剥掉了弱目标幅度低于门限被提前终止。把门限降到10之后检测率直接上到0.9以上。所以验证时不要只看最终目标数还要把残差画出来看是不是真的压到底了——残差里如果还有一块明显高出噪声底的鼓包说明迭代没跑完或者门限偏高。这个蒙特卡洛框架换场景很容易改targets、改cnr_db、改杂波分布类型就能得到不同信杂比下的CLEAN性能曲线。以后接手这个方向的人第一件事应该是把这套验证跑一遍再碰实测数据否则永远不知道算法好坏是被什么参数卡住的。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/31 5:44:30 · 阅读 1.2k 基于CLEAN算法的雷达杂波仿真与目标检测工程实现
简介资源包是一份面向雷达信号处理学习者与工程师的MATLAB仿真源码聚焦调频连续波FMCW雷达回波信号中的杂波抑制问题也可推广到其他频域杂波抑制场景。该仿真基于CLEAN算法在频域中对真实雷达回波数据进行处理通过相减函数抑制杂波分量从而还原目标信号相减函数当前采用余弦形式若信号为复数域可将其改成指数形式便于适配不同场景。整包仅含1个文件即MATLAB脚本源码体积约1KB代码短小精悍、无额外依赖复制到MATLAB环境中即可运行无需复杂配置很适合初学者对照学习CLEAN算法流程也可通过修改参数快速应用于自己的雷达回波数据。目前已有897人学习/下载说明其在雷达杂波抑制相关任务中具有一定参考价值。对正在做FMCW雷达仿真或频域杂波抑制项目的读者这份源码提供了一个可直接运行的算法模板既能帮助理解CLEAN算法的频域实现细节又能作为自定义开发的基础脚本尤其是采用真实回波数据时可作为快速验证的第一步。1. main()把雷达杂波仿真和CLEAN算法串起来的主心骨很多拿到雷达回波数据的人最先打开的是一个画图脚本而不是那个把数据读取、杂波仿真和CLEAN算法一步步串起来的main()。等你在距离-多普勒图上真正看到那一片杂波时才会意识到杂波不只是看起来“脏”它还会把弱目标死死按在强散射点的旁瓣里常规的恒虚警检测拿它没什么办法。这个工程方案的核心就是用main()搭一条完整链路——先在可控条件下仿真雷达杂波和回波数据再用CLEAN算法把杂波与旁瓣一层层剥掉最终输出目标列表。适合正在做雷达算法落地验证的工程师也适合想把仿真原型往Python/C工程里迁移的信号处理从业者。2. 先把杂波做出来从瑞利到K分布的雷达杂波仿真2.1 为什么仿真先于CLEAN你需要一块可控的“脏”回波直接拿实测数据调CLEAN最大的问题是“黑匣子”你不知道某个峰值到底是目标、杂波还是噪声更不知道算法是不是真的把杂波剥干净了。仿真杂波的价值在于目标位置、幅度、多普勒频移、杂波分布全部由自己指定CLEAN跑完一对比就知道对不对、差多少。后面要换场景只需要改几个参数不用扛着雷达外场跑一遍。杂波模型的选择一般看应用场景。均匀地面且雷达分辨率较低时大量独立散射点的叠加会让幅度趋近瑞利分布这也是教科书默认的模型。高分辨率雷达看海面地杂波的拖尾特别重强散射点密集这时候瑞利已经兜不住要用K分布或韦布尔分布。K分布的形状参数越小拖尾越重就越接近“少量强散射体加上大量微弱散射体”的海面实况。地杂波、海杂波在慢时间维上还有自己的多普勒谱通常建模成高斯谱中心频移对应杂波散射体的平均径向速度谱宽对应海浪起伏或天线扫描展宽。所以一套能用的杂波仿真至少要有三层可调幅度分布、多普勒谱形状、信杂比和信噪比。三层都控住生成的回波才能当作CLEAN算法的“测试靶场”。否则参数随机一拉出来的东西自己都不知道是什么那就是靠玄学调算法了。2.2 在main()里生成一帧含目标与杂波的IQ数据Python脚本与参数下面这段代码生成的是脉压后的距离-多普勒复图像横轴距离门纵轴慢时间脉冲每个元素是复数IQ。目标直接以“距离门多普勒频移信噪比”的形式插入杂波按CNR信杂比控制强度。import numpy as np def generate_frame(n_range1024, n_pulse128, prf1000.0, cnr_db30.0, clutter_typerayleigh, target_listNone): 生成一帧距离-多普勒域复数据维度 [n_pulse, n_range]。 target_list: [(range_bin, doppler_hz, snr_db), ...] 假定噪声功率归一为1。 rng np.random.default_rng(20240601) # 接收机噪声复高斯实部虚部方差各0.5总功率为1 noise (rng.standard_normal((n_pulse, n_range)) 1j * rng.standard_normal((n_pulse, n_range))) / np.sqrt(2) # 杂波基底复高斯得到瑞利幅度分布 clutter (rng.standard_normal((n_pulse, n_range)) 1j * rng.standard_normal((n_pulse, n_range))) / np.sqrt(2) # 慢时间维多普勒谱成形高斯谱中心-8Hz谱宽5Hz freq_axis np.fft.fftfreq(n_pulse, 1.0 / prf) sigma_f 5.0 center_f -8.0 filter_h np.exp(-0.5 * ((freq_axis - center_f) / sigma_f) ** 2) clutter np.fft.ifft(filter_h[:, None] * np.fft.fft(clutter, axis0), axis0) # 按信杂比归一化 noise_pow np.mean(np.abs(noise) ** 2) clutter_pow np.mean(np.abs(clutter) ** 2) clutter clutter * np.sqrt(10 ** (cnr_db / 10) * noise_pow / clutter_pow) frame noise clutter # 目标噪声功率为1幅度 10^(snr_db/20) if target_list: for rbin, df_hz, snr_db in target_list: dop_bin int(round(df_hz / prf * n_pulse)) amp 10 ** (snr_db / 20.0) frame[dop_bin % n_pulse, rbin] amp return frame这段代码的逻辑是先铺噪声底再叠加已经做过谱成形的杂波最后逐个放入目标。多普勒谱成形放在杂波之后、功率归一化之前是为了避免滤波改变杂波功率后还要二次校准。n_pulse和n_range决定帧尺寸prf决定多普勒轴刻度。target_list里的doppler_hz会被换成脉冲索引负数代表目标朝雷达运动正数远离。这个生成函数已经可以直接在main()里调用if __name__ __main__: # 两个距离接近的目标一个强一个弱用于验证CLEAN的旁瓣剥离能力 targets [(400, 0.0, 35.0), (403, 0.0, 20.0)] frame generate_frame(target_listtargets) # frame[0, :] 就是第一个脉冲对应的距离像两个目标距离只差3个门强目标脉压后的距离旁瓣正好压在弱目标的位置上这是后面CLEAN算法最典型的用武之地。3. CLEAN算法原理从射电天文的“脏图”到雷达的“脏回波”3.1 核心假设与迭代逻辑CLEAN算法最早来自射电天文Högbom在1974年提出用它重建综合孔径图像。观测到的“脏图”被认为是真实天空亮度分布与天线方向图点扩散函数PSF的卷积再加上噪声。CLEAN的做法是在脏图里找最亮的点认为它是真实源把这个亮点的PSF分量按一定比例减掉记录该源的幅度和位置然后再找下一个最亮点重复迭代直到剩余信号低到门限以下。这套逻辑放到雷达里完全成立。脉冲压缩后一个理想点目标在距离维上是一个sinc形状的响应主瓣旁边挂着一串距离旁瓣。如果同一个距离门附近有个强目标它的旁瓣会把弱目标盖住。雷达回波里的目标、强杂波散射体都可以看成稀疏分布的“点源”而CLEAN正是把这些点源的旁瓣逐个剥离出来的迭代算法。需要明确一点CLEAN不是超分辨算法。两个目标只要落进同一个主瓣宽度内峰值搜索就分不开它们CLEAN也一样无能为力。它的价值在于把“主瓣宽度之外被旁瓣掩盖”的目标找回来而不是突破瑞利分辨极限。理解了这一点你就不会拿它去硬解主瓣内多目标。3.2 增益、PSF与门限三个决定收敛行为的数字第一个数是循环增益loop gain记作g。标准CLEAN每次只减掉峰值幅度的g倍而不是全额减掉。g太小迭代次数多收敛慢g太大PSF与真实响应只要有一点不匹配误差就会被放大形成目标两边的虚假亮点。我一般取0.2~0.4。如果PSF是从标定数据测来的、质量很好可以放宽到0.5如果PSF只是理论近似就老实取0.25以下。第二个数是PSF本身。在雷达里它就是脉冲压缩后点目标响应的归一化复序列峰值恒为1。PSF必须覆盖到主瓣以外至少4~6个旁瓣否则每次减除都残留一段尾巴多次迭代后残留尾巴会累积成梳状结构。PSF的获取有两种从发射信号的自相关函数推导或者用强反射体的实测脉压数据截取。实测数据比理论模型可靠得多。第三个数是停止门限。CLEAN迭代不能无限跑否则会把噪声也当目标剥出来。门限通常设成“噪声底以上10~15dB”噪声底用最远距离门无目标区的幅度中位数来估计。这个门限要动态计算绝不能写死成某个绝对幅度值数据换一段就变。3.3 一维CLEAN的迭代式与伪代码一维CLEAN的迭代公式可以写成r^(k1)(t) r^(k)(t) - g · A_k · h(t - t_k)其中r为残差信号t_k为第k次迭代找到的峰值位置A_k是峰值位置的复幅度h是PSF。每次迭代把PSF平移到t_k处乘上g和A_k再从残差中减掉同时在目标列表里记录一个幅度为g·A_k的点目标。伪代码就是下面这个结构resid data.copy() peaks [] while True: idx np.argmax(np.abs(resid)) amp resid[idx] if np.abs(amp) threshold: break peaks.append((idx, amp)) # 以idx为中心将PSF缩放后减去 resid resid - gain * amp * shift_psf(psf, idx)注意amp是复数shift_psf把PSF数组平移到idx位置。复数运算必须全程保持不能只对幅度做操作否则相位信息丢得一干二净。4. 在main()里实现CLEAN从数据读取到目标列表输出4.1 main()的职责参数、文件读取与按帧循环main()在雷达处理链路里不是“随便写个函数入口”那么简单。它至少要负责三件事解析参数、按帧读取数据、把处理完的目标列表写出去。很多工程里CLEAN算法本身没问题翻车翻在数据路径写死、参数没传对、或者读到的IQ数据格式和算法预期不一致。main函数参数这块C/C工程里标准写法是int main(int argc, char *argv[])argc是参数个数argv是参数内容。我习惯把数据文件路径、PSF文件路径、迭代增益、停止门限全部通过参数传入而不是硬编码在源码里。下面这段是典型的框架import sys def main(): # argv[0]是程序名真正参数从argv[1]开始 if len(sys.argv) 3: print(usage: radar_clean iq_data psf.npy [gain] [thresh_db]) return 1 data read_iq(sys.argv[1]) psf np.load(sys.argv[2]) gain float(sys.argv[3]) if len(sys.argv) 3 else 0.3 thresh_db float(sys.argv[4]) if len(sys.argv) 4 else 12.0 hits, resid clean_1d(data, psf, gaingain, thresh_dbthresh_db) for idx, amp in hits: print(f{idx}\t{20 * np.log10(abs(amp)):.1f} dB) return 0 if __name__ __main__: sys.exit(main())read_iq这里没有展开常见做法是读int16交错IQ先拆实部虚部再转复数。关键点是argv处理一定要先判断长度否则用户忘了传参数时程序会直接崩这在自动化批处理里非常讨厌。数据路径、PSF路径、算法参数全部走参数传换数据就不需要重新编译。4.2 完整的CLEAN最小实现clean_1d与PSF构造现在给出可以直接抄走的clean_1d实现。PSF用sinc近似长度取主瓣宽度4个门、左右各5个旁瓣宽度这样能覆盖到远端旁瓣。def make_psf(width4, n_side5): 生成脉压点目标响应的sinc近似峰值归一化为1。 width主瓣过零点间隔对应的距离门数 n_side主瓣每侧覆盖的旁瓣周期数 x np.arange(-n_side * width, n_side * width 1) psf np.sinc(x / width) return psf / np.max(np.abs(psf)) def clean_1d(data, psf, gain0.3, thresh_db12.0, max_iter512): 一维CLEAN输入距离像输出目标列表和残差。 data: 一维复数距离像 psf: 一维复数PSF峰值已归一化 返回 (hits, resid)hits元素为 (距离门, 复幅度) pad psf.size # 两侧各补一个PSF长度的零省去边界索引判断 resid np.zeros(data.size 2 * pad, dtypecomplex) resid[pad:pad data.size] data # 噪声底取末尾5%距离门用幅度中位数估计 tail data[int(data.size * 0.95):] noise_level np.median(np.abs(tail)) / 0.6745 threshold noise_level * 10 ** (thresh_db / 20.0) half psf.size // 2 hits [] for _ in range(max_iter): peak_idx int(np.argmax(np.abs(resid))) amp resid[peak_idx] if abs(amp) threshold: break # 记录目标注意只记增益乘过的幅度 hits.append((peak_idx - pad, amp * gain)) # 以峰值位置为中心减掉缩放后的PSF start peak_idx - half resid[start:start psf.size] - gain * amp * psf clean_resid resid[pad:pad data.size] return hits, clean_residclean_1d的实现把PSF两侧补齐后再做减除这样peak靠近首尾时不需要写一堆边界分支。工程上对长数据会改成切片裁剪避免每次迭代处理整段补零信号但教学和验证场景这个补零版本更不容易写错。PSF构造里的width参数要和实际波形匹配线性调频信号时宽带宽积越大主瓣越窄width就越小。如果是从实测标定数据取的PSF直接np.load进来替换make_psf即可。提示clean_1d的噪声底估计依赖末尾5%距离门没有目标。雷达采集时专门留一段空回程就是为了给这类处理提供参考区。4.3 把仿真、脉压和CLEAN串进main()里的完整示例下面再把第2章的仿真函数和第4章的CLEAN函数接在一起跑通最小闭环import numpy as np def main_sim(): # 生成一帧含杂波和两个目标的距离-多普勒数据 targets [(400, 0.0, 35.0), (403, 0.0, 18.0)] frame generate_frame(n_range1024, n_pulse128, cnr_db30.0, target_listtargets) # 这里只取第一个脉冲的距离像做一维CLEAN演示 # 真实工程中会对每个多普勒门分别做一维CLEAN或直接做二维CLEAN data frame[0, :] psf make_psf(width4, n_side5) hits, resid clean_1d(data, psf, gain0.25, thresh_db12.0) for idx, amp in sorted(hits): print(frange_bin{idx:5d} amp{abs(amp):8.2f} f{20 * np.log10(abs(amp)):6.1f} dB) if __name__ __main__: main_sim()运行后会看到强目标在400被检出403附近那个弱目标如果在原距离像里被旁瓣淹没经过CLEAN后会被单独挑出来。排序输出是为了方便和target_list里的真值对比。main_sim这个入口保持了main()风格的参数化设计场景和目标列表都通过变量传入换实验条件不用改内部逻辑。真实工程里距离-多普勒图要完整处理不能只取第一个脉冲。常见做法是沿多普勒维逐行做一维CLEAN或者把PSF扩展成二维距离维用sinc多普勒维用主瓣宽度为1的理想窗然后直接在二维矩阵上搜峰、减二维PSF。二维CLEAN的计算量大约是每迭代一次做一次二维卷积级减除数据帧大时要考虑性能。5. CLEAN算法避坑与常见问题5个翻车现场的修复清单5.1 参数参考迭代增益、门限、PSF长度的初始推荐先给一张参数速查表新场景跑不通时按这组初始值试。参数推荐初值调节方向loop_gain0.25残差出现振荡就减小收敛慢可适当加大到0.4thresh_db12 dB虚警多就加大漏检多就减小PSF长度主瓣宽度5~8倍残差里出现梳状旁瓣就加长max_iter目标数估计值的3~5倍到上限还没收敛检查门限和增益噪声底估计区末尾5%距离门该区有目标时换到盲区/空回程段这几个参数互相牵制。最常见的问题是门限设太低导致迭代停不下来max_iter被耗尽后输出几百个假目标或者门限设太高弱目标还没被剥出来就终止了。调试时先把gain固定0.25只调thresh_db一次只动一个变量不然出了问题根本不知道是谁的锅。5.2 两个与main()直接相关的编译期问题先说C/C场景里很常见的“编译器未包含main类型”报错。现象是编译时报错说main没有返回类型或者链接时提示找不到main入口。原因通常有三类一是把main写成了void main()这在很多C标准里直接判错二是main被放进了某个没参与链接的源文件三是拼写错误main拼成了mian。解决方法是统一写int main(int argc, char *argv[])并在工程文件树里确认入口文件确实参与编译。第二个是main函数参数处理不到位。雷达程序里argv通常要接收数据文件路径、PSF路径和算法参数。如果代码里直接用argv[1]而不检查argc用户少传一个参数程序就段错误。工程上一定要先判断argc数量给足默认值。这算老生常谈但我在实际项目里见过太多因为这问题导致的启动崩溃尤其在同一套代码跑不同批次数据的时候。5.3 五个翻车场景现象、原因、解决第一个翻车场景迭代结束后残差里还有明显的梳状旁瓣结构。原因是PSF截断太短只覆盖了主瓣附近远端旁瓣每次都没减干净累积起来形成周期性的残留。解决是把PSF长度扩大到主瓣宽度的5倍以上必要时从标定数据里取实测PSF。第二个翻车场景门限写死成一个固定幅度换一段数据后输出目标数量剧烈变化。原因是数据换了之后噪声底变了但门限没跟着变。解决是用末尾无目标区的中位数动态估计噪声底门限始终取“噪声底thresh_db”。第三个翻车场景gain设成0.6结果强目标两边出现对称的假目标。原因是增益太大PSF和真实响应只要偏差一点点减除后的残差里就会留下反向的伪峰下一轮又被当成真实目标剥出来。解决是gain降到0.2~0.3如果已经产生了对称假目标把它们合并掉再用低增益重新迭代。第四个翻车场景复数数据用了实数sinc做PSF目标恢复出来的相位乱跳、幅度偏低。原因是复数信号只减实数PSF实部虚部被同步缩放破坏了原始相位关系。解决是全程使用复数PSF和复数减除即使理论sinc是实函数也要保证相位数据流不被拆散。第五个翻车场景两个目标靠得很近弱目标直接被当作强目标的旁瓣一起减没了。原因是CLEAN只做单峰搜索主瓣内多目标无法分辨。解决是先做一次CFAR预检测对主瓣展宽区域改用“多目标最小二乘”模型同时估计两个复幅度而不是让CLEAN盲剥。6. 用合成数据验证CLEAN性能目标恢复误差与蒙特卡洛检查验证CLEAN做得对不对最直接的办法是拿第2章的仿真器生成已知真值的数据跑完CLEAN后对比输出和真值。位置误差不超过2个距离门、幅度误差不超过3dB视为命中非真值位置的检出都算虚警。每轮换随机种子重复100次统计检测概率和虚警率。下面是验证脚本的核心片段def run_trial(n_trials100): total_true 0 total_hit 0 total_fa 0 psf make_psf(width4, n_side5) for seed in range(n_trials): rng np.random.default_rng(seed) targets [(400, 0.0, 30.0), (406, 0.0, 20.0)] frame generate_frame(n_range1024, n_pulse128, cnr_db25.0, target_listtargets) data frame[0, :] hits, _ clean_1d(data, psf, gain0.25, thresh_db10.0) truth [400, 406] hit_pos [idx for idx, _ in hits] for t in truth: total_true 1 if any(abs(t - idx) 2 for idx in hit_pos): total_hit 1 detected sum(1 for idx in hit_pos if all(abs(t - idx) 2 for t in truth)) total_fa detected print(fdetection rate{total_hit / total_true:.2f}, ffalse alarms per frame{total_fa / n_trials:.2f})我自己的血泪经验是第一次跑这组验证时把thresh_db设到18检测率只有0.5以为算法有问题。后来逐帧看残差才发现强目标早就被剥掉了弱目标幅度低于门限被提前终止。把门限降到10之后检测率直接上到0.9以上。所以验证时不要只看最终目标数还要把残差画出来看是不是真的压到底了——残差里如果还有一块明显高出噪声底的鼓包说明迭代没跑完或者门限偏高。这个蒙特卡洛框架换场景很容易改targets、改cnr_db、改杂波分布类型就能得到不同信杂比下的CLEAN性能曲线。以后接手这个方向的人第一件事应该是把这套验证跑一遍再碰实测数据否则永远不知道算法好坏是被什么参数卡住的。希望帮到你。本文还有配套的精品资源点击获取