简介:全局法图像垂直条纹去除是一份面向遥感图像处理、计算机视觉方向开发者的Matlab实用程序,用于对高光谱图像或普通RGB图像中出现的规则垂直条纹噪声进行全局建模与消除。该思路适用于成像传感器像元响应不一致导致的条带伪影,还针对倾斜条纹场景提示先旋转再处理,并提供了无循环高速版与循环低速版两套实现,经高光谱数据对比,高速版运行效率可提升约20倍。资源压缩包共2个文件,均为Matlab脚本(.m),体积仅2KB,结构精简,注释详尽,便于直接运行、核对算法流程,也可作为理解全局条纹去除原理的入门范例。当前已有886人学习下载,适合希望快速验证算法效果、对比不同实现版本或在此基础上二次开发的工程人员与研究人员。
1. 垂直条纹为什么难缠:空域滤波失效背后的全局结构问题
做图像算法的人,几乎都在某个时刻被「图像垂直条纹」折磨过:红外热像仪画面上一道道竖向亮暗纹,工业线扫相机采集的图每隔几列就有一条,遥感影像里经常整景都是条带噪声。这类噪声最烦人的地方在于,你拿中值滤波、高斯滤波、双边滤波去试,要么滤不干净,要么把边缘和纹理一起抹花了,属于典型的「越修越糟」。
原因在于垂直条纹是覆盖整幅图像的周期结构,而空域滤波只能看到局部窗口——窗口小了滤不掉,窗口大了伤细节,这是结构性矛盾。处理这类问题必须换到全局视角:把图像变换到频率域,在频谱上把条纹对应的能量精准剔除,或者在全图范围内拟合列方向的全局趋势做修正。这份资源就是围绕「全局法」整理的一套可执行方案,包含频域陷波滤波与列回归校正两条完整路线,以及我踩过的各种坑。适合做红外图像处理、工业视觉检测、遥感影像预处理和低剂量 CT 重建的工程师,也适合刚接触傅里叶变换但想直接把它落地到真实图像的初学者。
2. 频域视角看条纹:傅里叶变换如何定位周期噪声,陷波器又该怎么设计
2.1 垂直条纹在频谱上留下的记号:一眼认出噪声峰
要理解全局法,先得知道垂直条纹在频域里长什么样。这里直接说结论:一幅 M×N 的图像 I(x,y) 上叠加了垂直条纹,条纹沿 y 方向延伸、灰度沿 x 方向呈周期变化,可以写成 n(x) = A·sin(2πu₀x + φ) 的形式。对它做二维傅里叶变换后,这个正弦分量的能量会集中到频域平面上的两个点 (±u₀, 0) 附近——注意第二个坐标是 0,所以噪声峰全部落在频谱图的水平中轴线上,也就是 v=0 的那一行。
在机器视觉和工业相机采集场景里,图像傅里叶变换是分析周期噪声的第一选择,原因就在于此:空间域里横跨整幅图的条纹,在频域里只是两个孤立的峰,处理起来成本极低。用 fftshift 把零频移到中心后,这两个峰就对称出现在原点左右两侧的水平轴线上。峰位到原点的距离和条纹周期直接相关:图像宽度为 W 像素,条纹周期为 T 像素,那么峰位距原点的距离大约是 W/T 个频域像素。举个例子,512 像素宽的图像上有周期为 16 像素的条纹,峰位就在距原点 32 个像素的位置。
实际图像往往不只一组条纹。传感器读出噪声、电源纹波、扫描机构抖动可能叠加出多个周期分量,频谱水平轴线上会依次排开多个亮峰,二次谐波、三次谐波也可能出现。这些峰的高度通常远高于周围频谱成分的平均水平,这就是后面自动诊断算法的依据。需要特别说明的是,空域滤波处理不了这类噪声的根本原因也在这里:条纹是全域周期信号,能量高度集中在特定频率,你在空间域无论选多大窗口都只是在跟一个「无处不在」的信号较劲,而频域只需动几个点。
2.2 陷波器的设计:为什么不能只置零一个像素
明确了噪声峰的位置后,常规做法是在频谱上做一个带阻滤波器,把峰及其邻域的能量压掉,叫陷波滤波(Notch Filter)。但这里有一个新手最容易犯的错:以为把峰值那个像素置零就行了。实际做出来你会发现条纹确实淡了一点,但整幅图像出现一圈圈水波纹一样的振铃,而且条纹并没有彻底干净。
原因在于离散傅里叶变换里,周期噪声的能量并不是理想地集中在一个点上,而是分布在以峰位为中心的一个小邻域内,形状近似一个窄峰。只置零单个像素相当于用一根极细的针去戳一个面,能量从旁边漏过去,同时这个锐利的置零操作在反变换时又引入了新的高频振荡,也就是振铃。正确做法是用一个平滑过渡的掩模,把峰位周围的一小片区域压下去。
我一般用高斯型陷波,它的传递函数形式是:H(u,v) = 1 - exp(-(D²(u,v) - D₀²)² / (2·W²·D₀²)),其中 D(u,v) 是当前像素到原点的距离,D₀ 是噪声峰到原点的距离,W 控制陷波的带宽。这个公式的好处是掩模从 1 到 0 是渐变过渡的,不会在反变换时引入新的振铃。下面这段代码演示如何生成一个这样的陷波掩模:
import numpy as np import cv2 def gaussian_notch_mask(shape, center, radius=2, bandwidth=5): """ 生成高斯型陷波掩模 shape: (H, W) 频谱尺寸 center: (u0, v0) 噪声峰位置 radius: 需要压制的邻域半径 bandwidth: 陷波过渡带宽 """ H, W = shape u = np.arange(H) - center[1] # 行方向相对中心 v = np.arange(W) - center[0] # 列方向相对中心 U, V = np.meshgrid(v, u) D = np.sqrt(U**2 + V**2) D0 = np.sqrt((center[0] - W//2)**2 + (center[1] - H//2)**2) # 高斯带阻:1 表示保留, 0 表示滤除 mask = 1 - np.exp(-((D**2 - D0**2)**2) / (2 * bandwidth**2 * D0**2 + 1e-8)) # 只压制以 center 为中心的小范围, 远处恢复为 1 mask[D > radius] = 1 mask[center[1] - radius:center[1] + radius, center[0] - radius:center[0] + radius] = mask[center[1] - radius:center[1] + radius, center[0] - radius:center[0] + radius] return mask这里几个参数需要讲清楚。radius 是压制范围,经验上取 2~4 个频域像素就够,因为能量扩散范围有限;bandwidth 控制过渡的陡峭程度,bandwidth 越小过渡越陡,越接近二值掩模,bandwidth 越大过渡越平缓,但可能波及旁边真实的频谱分量。D0 的计算是为了让带阻的中心正好落在噪声峰上,公式里减去的 W//2 和 H//2 是因为 fftshift 后原点在频谱中心。实际使用中,我会对频谱峰值点周围的几个像素逐一测试,看反变换后条纹能量下降多少、细节损失多少,以此决定 radius 的取值,这比纠结理论公式更快。
3. 直接能跑的处理流水线:频谱自动诊断、陷波生成与图像重建
3.1 完整处理流程:从读图到输出干净的图像
理论铺垫完之后,这里给出一套可以直接跑的完整方案。它的流程是:读取灰度图 → 计算中心化傅里叶变换 → 在频谱水平轴线上自动寻找噪声峰 → 生成陷波掩模 → 频谱滤波 → 反变换回空间域。核心脚本如下:
import numpy as np import cv2 def remove_vertical_stripes(image, peak_thresh=8, notch_radius=3, border_width=3): """ 频域陷波法去除垂直条纹 image: 输入灰度图, uint8 或 float, 建议先转 float peak_thresh: 噪声峰判定阈值, 相对频谱均值的倍数 notch_radius: 陷波压制半径, 单位是频域像素 border_width: 沿水平轴线检测峰的半宽 """ # 1. 转 float 并计算中心化频谱 img_float = image.astype(np.float32) F = np.fft.fft2(img_float) F_shift = np.fft.fftshift(F) mag = np.abs(F_shift) # 2. 在水平中轴线上找噪声峰 H, W = mag.shape cy, cx = H // 2, W // 2 # 取中轴线附近一小条带, 降低噪声干扰 strip = mag[cy - border_width: cy + border_width + 1, :] profile = strip.mean(axis=0) # 沿竖直方向平均, 得到一维频谱能量 baseline = profile.mean() # 排除原点附近的低频区, 避免把图像本身亮度误判为条纹 profile[:cx - 40] = 0 profile[cx + 41:] = 0 profile[cx - 10: cx + 11] = 0 # 挖掉 DC 区域不要检测 # 找所有超过阈值的局部极大值点 threshold = baseline * peak_thresh peaks = [] for i in range(1, W - 1): if profile[i] > threshold and profile[i] >= profile[i-1] and profile[i] >= profile[i+1]: peaks.append(i) # 合并距离过近的峰, 只保留能量更高者 peaks = _merge_close_peaks(peaks, profile, min_dist=6) # 3. 构造陷波掩模并滤掉所有峰 mask = np.ones_like(mag) for peak_x in peaks: mask = _apply_notch(mask, (peak_x, cy), notch_radius) # 4. 滤波并反变换 F_shift_filtered = F_shift * mask img_filtered = np.fft.ifft2(np.fft.ifftshift(F_shift_filtered)).real return np.clip(img_filtered, 0, 255).astype(np.uint8), peaks这里_merge_close_peaks和_apply_notch是两个辅助函数:前者把相互距离小于 6 个像素的峰合并,只保留谱线能量更高的那个,因为条纹周期非常接近时会出多个相邻峰,全滤会导致带宽过大;后者在指定位置生成高斯陷波并乘到掩模上。peak_thresh 是最难调的参数,取 8 意味着「峰的能量至少是全频谱均值的 8 倍」才认为是条带噪声。这个值对大部分红外图像和遥感影像都适用,但如果你处理的图像对比度特别高、细节纹理很丰富,建议先看频谱剖面再决定。
3.2 自动找峰阈值怎么定,掩模参数怎么改
上面代码里 peak_thresh 和 notch_radius 是两个最关键的手工参数。不少第一次用的人会直接把 peak_thresh 设得很低,想「宁可错杀不可放过」,结果把图像本身的横向纹理也滤掉了。我一般会让用户先跑一段频谱诊断脚本,把水平轴线上的能量分布打印成数值列表,肉眼看清楚噪声峰和正常频谱分量的高度差,再定阈值。
阈值合理性有个经验区间:对典型的红外焦平面条纹噪声,峰值通常是频谱均值的 15~30 倍,peak_thresh 取 8~15 都行;对比较微弱的光学条纹,峰值可能只高 5~8 倍,阈值就得降到 3~4。但注意,阈值降到 3 以下时,图像自身横向边缘(比如地平线、建筑轮廓)也会被当成峰,这种情况我会改用第 4 章的列回归法,而不是硬调阈值。notch_radius 建议从 3 开始,检查反变换结果后按需增减 1。
另外要注意处理顺序:如果图像里叠加了多组不同周期的条纹,一次把检测到的所有峰全部滤波往往比逐次处理效果更好,因为各峰之间的陷波区域不会重叠。真正需要迭代处理的是那种峰位恰好落在另一个峰的过渡带里的情况,这时先滤能量大的峰,再重新检测剩余的小峰,通常两轮就干净了。
3.3 滤波完怎么验收:三分钟看三个信号
代码跑完别急着存图,先做三件事确认效果。第一步看列均值曲线:对去条纹后的图像逐列求平均,一条平滑的曲线说明列间亮度连续,如果再出现锯齿状抖动说明滤得不彻底。第二步看残差图:把去条纹前后的图像相减,残差图上应该只看到竖条纹的痕迹,如果残差里出现明显的物体轮廓,说明陷波带宽把真实细节也压掉了,需要减小 notch_radius。第三步看频谱:对去条纹后的图再做一次 FFT,水平轴线上应该看不到突出的峰值。
这三步可以用一个简单的评估函数串起来,我通常把列均值差分绝对值之和作为量化指标,这个指标在本文第 6 章还会展开讲。如果你在这个阶段发现条纹反而变得更明显,多半是掩模中心位置算错了——确认一下你的代码里有没有用 fftshift,以及峰值坐标是否和掩模生成时使用的坐标系一致,这是最常翻车的地方。
4. 频域解决不了的部分场景:列回归与分块校正的思路
4.1 非周期条纹:为什么频谱里找不到明显的峰
频域陷波不是万能的。红外焦平面阵列的非均匀性噪声、线扫相机的列暗电平不一致、某些 CMOS 传感器的列读出噪声,这类条纹的特征是周期不固定、列与列之间的偏差是随机或缓变的,频谱上不会出现干净利落的尖峰,而是弥散在低频区域,陷波滤波器无从下手。
这类条纹在空间域的特征反而更明显:逐列统计灰度均值后,正常图像的列均值曲线是平滑变化的(因为场景亮度在横向是渐变的),而带条纹图像的列均值曲线上叠加了明显的逐列跳变,像锯齿一样。于是有了另一条全局法思路——列回归校正:拟合出列均值曲线中的低频趋势部分,把它当作「场景真实亮度」,用原始列均值减去这条趋势线,剩下的残差就是列方向的噪声偏差,最后把残差从每一列中减掉。
这个思路在工业界叫列均一化或 flat-field correction 的简化版,它对非周期性列噪声的效果比频域法可靠得多,而且计算量小,适合批量处理。
4.2 列均值拟合与逐列修正:代码与参数细节
核心实现不复杂,就是统计、拟合、修正三步。下面这段代码用多项式拟合列均值趋势,并完成逐列修正:
import numpy as np def column_regression_correction(image, degree=4, block=1): """ 列回归去垂直条纹 image: 输入灰度图 float degree: 多项式拟合阶数, 推荐 3~5 block: 分块列数, 1 表示全局拟合; 图像亮度分布复杂时设为 128 """ H, W = image.shape img_out = image.copy() # 分块处理, 每 block 列作为一个拟合区间, 避免场景亮度突变干扰 for start in range(0, W, block): end = min(start + block, W) col_means = img_out[:, start:end].mean(axis=0) # 每列均值 col_positions = np.arange(start, end) # 多项式拟合低频趋势 coeffs = np.polyfit(col_positions, col_means, degree) trend = np.polyval(coeffs, col_positions) # 残差 = 实际列均值 - 低频趋势, 即列噪声 residual = col_means - trend # 逐列减去残差 img_out[:, start:end] -= residual.reshape(1, -1) return np.clip(img_out, 0, 255) corr_img = column_regression_correction(img.astype(np.float32), degree=4, block=256)多项式阶数 degree 是这里最需要小心的参数。阶数取 1~2 时拟合的是全局线性渐变,适合亮度均匀的图像;取 3~5 时可以表达场景中横向的明暗过渡,比如红外图像中间亮四周暗;但如果超过 6,拟合曲线会开始跟着条纹本身的抖动走,趋势线里混入噪声,残差减小,条纹就去除不干净了。分块 block 参数处理的是场景里有明显亮度分区的情况,比如一行图里有天空有地面,亮度相差很大,全局一条多项式曲线无法同时拟合两段不同亮度,分块后每段单独拟合就自然多了。我一般先用 block=256 试,如果分块边界处出现亮度跳变,就减小到 128 或 64。
这里还有一个容易忽略的点:列均值用的是算数平均,如果图像里有高亮饱和点或坏点,会污染该列的均值估计,导致这一列修正过头。做之前先用中位滤波或百分位截断把极端值压一下,代价小、收益明显。这也解释了为什么有些场景下用列中位数代替均值会更稳——中位数对离群点不敏感,但计算量大得多,按需选择。
4.3 频域法与列回归法怎么选:一个对比表
实际项目中两条路线可以互补,选型主要看条纹形态。这里给一个对比表,方便你根据现象快速决策:
| 对比维度 | 频域陷波法 | 列回归法 |
|---|---|---|
| 适用噪声 | 周期性条纹、固定周期条带 | 非周期列噪声、增益不均 |
| 频谱表现 | 水平轴线有明显尖峰 | 无明显尖峰, 低频弥散 |
| 细节保留 | 只压制特定频率, 细节保留好 | 依赖拟合阶数, 阶数高会伤横向渐变 |
| 参数数量 | 阈值+半径, 需要看频谱诊断 | 阶数+分块, 需要看列均值曲线 |
| 典型场景 | 遥感影像条带、扫描仪周期噪声 | 红外非均匀性、线扫暗电平不一致 |
| 风险点 | 误伤横向纹理、振铃 | 误把真实场景渐变当条纹抹掉 |
一句话总结我的选型习惯:先做频谱诊断,水平轴线上有清晰尖峰的走频域陷波;频谱上看不出峰的走列回归。遥感影像里的传感器条带噪声大多是周期性的,低剂量 CT 重建出现的条状伪影也偏周期结构,这两类用频域法更干净;红外热像仪和工业线扫相机的列噪声则更接近随机偏差,直接上列回归法更合适。
5. 避坑与常见问题:五个翻车现场与排查顺序
5.1 忘了 fftshift,掩模位置全错,图像越滤越花
现象:滤波后的图像整个变暗或出现大范围条纹状残留,和原图相比变得模糊不清,检测到的峰值坐标看起来「完全对不上」。
原因:fft2 的输出中零频在左上角 (0,0),如果不做 fftshift 就去找峰、生成掩模,所有坐标都偏移了半个图像尺寸。掩模没有对准任何真实噪声峰,滤波相当于把低频能量和高频噪声同时削弱了。
解决:生成频谱后立刻 fftshift,得到以中心为原点的坐标;构造掩模时同样以中心为原点;反变换前再用 ifftshift 把频谱还原。记住这个铁律:fftshift 和 ifftshift 必须成对出现,中间所有坐标运算都基于中心化之后的坐标,否则等于白做。
5.2 陷波只置零一个像素:条纹没下去,还出了振铃
现象:条纹颜色淡了一点但轮廓还在,图像边缘处出现沿着轮廓扩散的水波纹,高频细节看起来脏兮兮的。
原因:周期噪声的能量在离散频谱里分布在峰位周围至少 2~3 个像素的邻域内,单点置零只去掉了很少一部分能量;同时锐利的单点截断在反变换时产生吉布斯现象,也就是振铃。
解决:改用高斯型或巴特沃斯型陷波掩模,把峰位周围覆盖起来。notch_radius 从 3 起步,逐步增加 1,观察残差图。如果振铃严重而条纹没净,还要检查是不是掩模只有 0 和 1 两个值——二值掩模很容易引发振铃,换成渐变过渡的掩模立刻好转。
5.3 陷波带宽太宽:条纹清了,图像也糊了
现象:条纹确实消失了,但整幅图像像蒙了一层纱,边缘不再锐利,横向纹理(比如云层、水面波纹)明显变淡。
原因:带宽设太大时,陷波器把条纹峰附近的正常频谱分量也压制了。尤其当图像自身有横向结构时,它的频谱能量也分布在水平轴线附近,和噪声峰重叠,一起被滤掉了。
解决:把 notch_radius 降回 2~3,或者把带宽参数调窄,让陷波只覆盖峰的主瓣。用第 3.3 节的残差检查法确认:残差图里如果出现物体轮廓,说明带宽过宽。遇到图像自身横向纹理很强的情况,建议换列回归法,不要在频域里硬抠。
5.4 拿一张图的频谱掩模去滤一组图:每张都出问题
现象:同一批采集的图像,A 图用这套参数去条纹效果很好,B 图用了同一套掩模后出现新的条纹或者图像变糊。
原因:条纹周期不是恒定不变的。传感器温度变化、积分时间调整、电源波动都会让条纹频率发生轻微漂移,峰位在频谱上会移动几个像素。固定掩模对不上新图的峰位,自然滤不干净,还可能压偏。
解决:对每帧图像重新做峰值检测,至少做一次频谱诊断。如果帧率很高、逐帧诊断成本大,就用多帧合并的方式生成一个公共掩模——把一组帧的频谱按像素取中位数再找峰,这样得到的掩模对整组图都有适用性,具体做法见下一章。
5.5 对 uint8 图像原地做减法:输出出现新的方块噪声
现象:列回归法跑完后,图像里出现一行行规律排列的亮暗块,或者原先的条纹没去掉又多了雪花点。
原因:uint8 是无符号 8 位整数,范围 0~255,做减法时低于 0 的值会溢出变成 255 附近的值,形成亮斑;最后如果用 np.clip 处理不当,这些溢出就是新的噪声源。列回归修正时残差里有正有负,直接在 uint8 上运算必然出问题。
解决:全流程用 float32 计算,输入先转浮点,所有减法在浮点域进行,最后统一 clip 到 [0, 255] 再转回 uint8。这里没有捷径,我第一次用 uint8 直接做列回归时就吃了这个亏,输出图里多出一堆方块噪声,排查了半天才发现是数据类型的问题。
6. 进阶:条纹能量指标与多帧自动掩模生成
6.1 一个能写进验收报告的条纹能量指标
去条纹做完了,怎么跟上下游交代效果?PSNR 和 SSIM 对整幅图敏感,但条纹是细结构,它们未必能灵敏地反映条纹能量的变化。我习惯用一个专门指标:列均值曲线的相邻差分绝对值之和。它衡量的是「列与列之间亮度跳变的剧烈程度」,条纹越重,这个值越大;去干净了,这个值会显著下降。下面这段代码可以算出来:
def stripe_energy(img): """ 条纹能量: 列均值曲线相邻差分绝对值之和 值越大说明列间亮度跳变越剧烈, 条纹越重 """ col_means = img.astype(np.float32).mean(axis=0) diff = np.abs(np.diff(col_means)) return float(diff.sum()) before = stripe_energy(img) after = stripe_energy(result) print(f"去条纹前: {before:.2f}, 去条纹后: {after:.2f}, 下降: {(before-after)/before*100:.1f}%")这个指标对周期性和非周期性条纹都有效,而且计算一次只要几毫秒。实践中,红外条纹图像去完后这个值通常会下降 60%~85%,如果下降不到 30%,说明滤波没滤干净或者参数没找到点上。把它写进验收报告,比截图更有说服力。需要注意的是,这个指标要和残差检查配合使用——防止为了降低指标而把图像整体模糊化,那种做法指标很漂亮但细节全没了。
6.2 多帧自动掩模:让一组图像共享同一个陷波器
批量处理一组图像时,逐帧诊断很烦。我的做法是先取 10~20 帧有代表性的图,分别做 FFT,取幅度谱的逐像素中位数合成一张「平均频谱」。中位数操作可以压掉单帧图像自身的纹理结构,只留下所有帧共有的周期噪声峰。然后在这个合成频谱上做峰值检测,生成统一的陷波掩模,应用整组图像。这样参数只需要调一次,后续全部自动。如果传感器工作状态改变(比如换了积分时间),重新采集一组帧再算一次掩模就行。
从这个角度说,频域陷波法的核心其实不在「滤波」本身,而在「定位噪声峰的准确性」。定位准了,剩下就是把峰值区域轻轻盖住。我处理过的项目里,凡是效果翻车的基本都是定位环节偷懒了——没有做频谱诊断就套参数。从那以后,我每次拿到一张条纹图,都会先花三分钟跑一遍频谱诊断,看一眼水平轴线上的峰分布再决定走频域还是列回归,这个习惯帮我省下了大量反复调参的时间。希望帮到你。
本文还有配套的精品资源,点击获取