看到“Strong non-uniformity correction algorithm based on spectral shaping statistics and LMS”这个标题我第一反应是终于有人把自适应滤波和频域统计揉在一起做非均匀性校正了。翻译成中文就是“基于光谱整形统计与LMS的强鲁棒性非均匀性校正算法”解决的是红外焦平面阵列IRFPA每个像元响应不一致带来的固定图案噪声问题。这类算法的核心痛点很明确红外探测器每个像元的增益和偏置天生不同导致成像后出现竖向条纹、斑点状的固定噪声实际中既不能用黑体反复标定又不能等噪声自己消失。基于场景的校正方法是目前的重点方向而LMS这种经典自适应滤波加上“光谱整形统计”这种在频域做文章的思路正好是避开挡片、实时校正的常见方案。这篇文章适合三类人刚接触红外图像处理、正在看非均匀性校正论文的研究生做红外整机算法、需要在实际视频流里做实时校正的工程师以及对自适应滤波、频域整形这类组合方法感兴趣、想找可复现代码的算法爱好者。展开之前先说明我复现的是这个思路的简化版本——真正工程落地的版本还会有盲元替换、增益范围约束等细节但核心框架和原文思路一致。1. 先从标题说起这个算法到底要解决什么1.1 标题拆解与翻译标题拆开来看Strong这里不是“强大的”这么简单更准确的翻译是“强鲁棒性”。意思是这个算法在场景剧烈变化、有局部遮挡、有急速运动物体时依然能稳定校正不出现明显的鬼影或错误收敛。Non-uniformity correction非均匀性校正英文简称NUC。这是红外成像里避不开的环节。Spectral shaping statistics光谱整形统计。这里的“光谱”不是光谱仪测的那种光学波长分布在红外图像里通常指空间频率域也就是把一帧图像做傅里叶变换或DCT变换后得到的不同频率成分的统计特性。整形就是对频域系数做权重调整。LMSLeast Mean Square最小均方算法自适应滤波里最经典的梯度下降类方法。整体翻译基于光谱整形统计与LMS的强鲁棒性非均匀性校正算法。1.2 红外非均匀性问题的根源非均匀性是一个物理层面的问题。红外焦平面阵列由成千上万个独立探测像元组成每个像元的光电响应曲线并不完全一致。即使对着一个均匀黑体辐射面各个像元输出的灰度值也会有差异。这种差异通常建模成一个线性关系y a × x b其中x是入射辐射通量y是输出灰度值a是增益b是偏置。每个像元都有自己的a和b而且这两个参数会随着温度、时间的漂移而缓慢变化。这种响应不一致在图像上表现得很直观画面里出现固定的横竖条纹、网状图案或斑点而且不随场景移动而移动。因为噪声形态像一个固定图案所以叫固定图案噪声FPN。早期解决方法是两点或三点定标法拿不同温度的黑体去给每个像元做线性校正。效果不错但有两个问题一是需要遮挡光路拍摄期间会中断二是标定系数会漂移过一段时间又要重新标定。所以业界更倾向于基于场景的校正方法不依赖黑体利用视频序列的统计规律自动估计增益和偏置。这个标题里的算法就是走这条路线。1.3 为什么需要“强鲁棒性”很多场景校正算法有一个致命弱点运动物体留下的鬼影。比如一个行人从镜头前走过算法如果错误地把他当成固定图案噪声去估计行人的“残影”会留在校正结果里很长时间而且是反相的暗影。“强鲁棒性”就是要减轻这类问题。光谱整形统计负责判断哪些频率分量是真实场景动态变化的部分哪些是稳定的固定图案噪声LMS则利用这个判断结果只对噪声分量做自适应更新不把物体边缘和运动区域卷进校正参数。2. 算法原理拆解光谱整形统计联合LMS2.1 响应模型与校正模型先说校正模型。假设原始输入图像为y每个像元的增益校正系数为g偏置校正系数为o校正后的输出为zz g × y o理想情况下z应该正比于场景真实辐射x而且比例系数对每个像元一致。所以问题变成了估计每像元的g和o。这和通常教材里的模型有一点差异。很多文献把校正写成z (y - o) / g等价但形式不同。这个标题提到的算法按我找到的实现思路用的还是加法减法形式因为LMS更新公式写起来更自然。训练过程中需要一个参考信号也就是期望输出d。对于非均匀性校正来说真实辐射x是不可知的所以参考信号必须从输入图像本身估计出来。这也是整个算法的关键如何从带有固定图案噪声的图像中分离出一个接近真实场景的参考图像。2.2 LMS如何更新增益与偏置LMS的核心思想很简单定义一个误差信号e d - z然后让校正参数沿着误差平方的负梯度方向移动。对增益g求偏导∂e²/∂g -2 × e × y对偏置o求偏导∂e²/∂o -2 × e于是LMS更新公式为g ← g μ_g × e × yo ← o μ_o × e这里μ_g和μ_o是步长因子。为什么增益步长通常要小于偏置步长因为y是灰度值数量级一般是几十到几百e和y相乘以后数值很大步长稍微大一点就容易导致增益震荡。偏置项的导数就是e本身没有y的放大效应所以可以取稍大一点的步长。把这个过程展开来看每个像素都在独立执行一个两权值的LMS自适应滤波输入向量是[y, 1]对应权重是[g, o]。不过是上百像素同时更新。但这个更新依赖d也就是参考图像。LMS本身不能产生参考图像必须由外部提供。2.3 光谱整形统计怎么估计参考信号光谱整形统计干的活就是制造一个参考图像d。红外图像里的固定图案噪声有一个明显的频域特性它主要集中在空间频率的低频区域或者集中在某些特定的谱线上。比如竖条纹噪声在频域里表现为水平方向的高频尖峰而真实场景通常是自然图像频谱分布更广、更随机。把图像分成8×8或16×16的小块做DCT变换每个块的低频系数尤其是DC系数和邻近它的几个低阶AC系数大幅受到固定图案噪声污染。对这些低频系数进行压缩或衰减再反变换回来得到的图像几乎不带固定图案噪声同时保留了边缘和细节。“光谱整形统计”中的统计部分在于单个频带是否被抑制不是人工拍脑袋定的而是根据一段时间的频谱统计来决定。固定图案噪声在某个频带的能量是稳定的、随时间基本不变的而场景纹理在那个频带的能量会随内容变化而剧烈波动。统计每个频带能量在时间轴上的方差方差小的频带大概率是固定噪声主导方差大的频带是真实场景主导。据此构造整形权重实现自适应整形。我复现时做了一个简化处理固定权重加滑动窗口能量统计修正。对于多数测试视频固定权重已经能达到不错的校正效果加上统计修正后收敛更稳定。2.4 两个模块的耦合逻辑光谱整形统计和LMS不是各自独立工作的它们形成闭环。每一帧处理的顺序是使用当前g、o对原始帧y做校正得到初步校正图像z。对z做分块DCT整形后重建得到参考图像d。计算误差e d - z用LMS公式更新g、o。下一帧重复以上过程直到g、o收敛。这个闭合环路里LMS输出的z越来越好进而光谱整形统计得到的d也越来越干净反过来又让LMS更新得更准确。缺陷在于如果d本身包含错误信息会把误差重新灌进参数里所以要对整形权重设计谨慎一些。3. Python复现从建模到跑通3.1 环境准备与数据生成我用的环境是Python 3.10 NumPy OpenCVOpenCV主要用来做DCT变换和图像卷积。先用仿真数据验证标准工作流。生成一组带有固定图案噪声的红外视频帧在代码里可控性更强能知道真实背景是什么定量评估校正效果。下面这段生成测试序列import cv2 import numpy as np def generate_test_frames(num_frames200, h256, w320, seed42): rng np.random.default_rng(seed) # 理想场景一张缓慢移动的圆盘 x np.zeros((h, w), dtypenp.float32) x[80:200, 120:240] 0.6 frames [] for k in range(num_frames): x_shift np.roll(x, int(k % 20), axis1) # 增益图平滑变化 像素级噪声 g 0.8 0.2 * np.exp(-((np.arange(h)[:, None] - h/2)**2 (np.arange(w)[None, :] - w/2)**2) / (2*40**2)) g g 0.03 * rng.standard_normal((h, w)) # 偏置图条纹型固定噪声 o 0.1 * np.sin(np.arange(w)[None, :] / 8.0) o o 0.05 * rng.standard_normal((h, w)) y g * (x_shift 0.05 * rng.standard_normal((h, w))) o frames.append(y) return frames, x生成的数据里g和o都包含了固定图案的成分场景本身是匀速平移的圆盘方便观察校正后的条纹是否被抑制。3.2 光谱整形统计模块实现分块DCT和频域整形是这个模块的核心。写一个类负责维护频谱统计同时生成参考图像class SpectralShaper: def __init__(self, block_size8, init_damping0.25, stats_window30): self.block_size block_size self.damping init_damping self.energy_sum None self.energy_count 0 self.stats_window stats_window def _build_shape(self, dct_block): bs self.block_size shape np.ones((bs, bs), dtypenp.float32) # 低频衰减DC系数保留很少低阶AC系数部分保留 shape[0, 0] self.damping shape[0, 1] 0.4 shape[1, 0] 0.4 shape[1, 1] 0.6 # 如果已经有统计信息修正权重 if self.energy_sum is not None: variance self.energy_sum / max(self.energy_count, 1) total variance.sum() 1e-6 normalized variance / total # 能量越集中在某些频带这些频带越可疑抑制越强 suspicious normalized (1.5 / (bs * bs)) shape[suspicious] * 0.5 return shape def update_statistics(self, dct_block): energy dct_block ** 2 if self.energy_sum is None: self.energy_sum np.zeros_like(energy) self.energy_sum energy self.energy_count 1 if self.energy_count self.stats_window: self.energy_sum * 0.9 self.energy_count int(self.energy_count * 0.9) def estimate_reference(self, z): bs self.block_size h, w z.shape ref np.zeros_like(z) for i in range(0, h - bs 1, bs): for j in range(0, w - bs 1, bs): block z[i:ibs, j:jbs] dct cv2.dct(block.astype(np.float32)) self.update_statistics(dct) shape self._build_shape(dct) shaped dct * shape block_ref cv2.idct(shaped) ref[i:ibs, j:jbs] block_ref return ref有几个细节需要注意。DCT变换之前要尽量让输入是平坦背景因为DCT对块边界不敏感但如果块内有大梯度边缘重建后会出现振铃。实际中用重叠块或者加窗可以减轻但计算量翻倍本复现为了清晰没有做重叠。块大小选8×8这是从JPEG时代就传下来的经验尺寸对大多数图像来说能在频率分辨率和空间分辨率之间取一个不错的平衡点。块太大固定图案噪声和真实纹理在频域混在一起难以分离块太小频谱统计不稳定统计意义减弱。3.3 LMS校正模块实现class LMS_NUC: def __init__(self, h, w, mu_g1e-4, mu_o1e-2): self.g np.ones((h, w), dtypenp.float32) self.o np.zeros((h, w), dtypenp.float32) self.mu_g mu_g self.mu_o mu_o def correct(self, y): return self.g * y self.o def update(self, y, ref): z self.correct(y) e ref - z self.g self.mu_g * e * y self.o self.mu_o * e return z, e这里y、z、e都是整幅图像矩阵LMS更新是对每个像元逐点进行的但用NumPy矩阵运算一次性完成。注意增益更新项里有一个e乘以y的因子如果y灰度值在0到255之间μ_g必须取得很小。如果图像先归一化到0到1则可以适当放宽。测试中我把μ_g设为1e-4μ_o设为1e-2收敛速度大概是50到80帧趋于稳定。如果μ_g取1e-3增益在100帧后开始出现细密颗粒噪声说明步长过大。3.4 主循环与参数调整def run_nuc(frames): h, w frames[0].shape shaper SpectralShaper(block_size8, init_damping0.25) nuc LMS_NUC(h, w, mu_g1e-4, mu_o1e-2) corrected_frames [] for idx, y in enumerate(frames): # 先校正 z nuc.correct(y) # 用校正结果估计参考图像 ref shaper.estimate_reference(z) # 更新增益与偏置 z, e nuc.update(y, ref) corrected_frames.append(z) if idx % 50 0: print(fframe {idx}: e_mean{abs(e).mean():.4f}) return corrected_frames, nuc.g, nuc.o这一步的编排逻辑我一开始写反了直接用原始帧y去估计ref结果校正收敛很慢。后来改成先用当前校正参数得到z再对z做频域整形误差里才真正包含了残留固定噪声的信息收敛速度明显提升。参数调整的经验是damping值控制低频保留程度太大会残留竖条纹太小会让图像整体变暗、细节丢失。0.25到0.35之间比较合适。μ_g和μ_o的比例大约在1:100左右对这个输入范围比较稳定。收敛后增益图像g应该接近真实的gain图偏置图像o接近真实offset图。在仿真数据里可以直接对比。3.5 评估指标与量化结果评估非均匀性校正效果常用粗糙度指标符号记作ρ。定义是图像与水平方向平移一格后的差值的绝对值之和除以原图像绝对值之和ρ || h1 ⊛ z ||₁ / || z ||₁其中h1 [1, -1]是水平差分核。ρ越小说明图像越平滑条纹越少。但要注意这个指标有一个陷阱如果把整幅图像过度平滑ρ也会很小但图像细节也没了。所以必须结合原始清晰度指标一起看。在仿真数据上我同时计算了PSNR因为真实场景x是已知的。结果大致如下指标校正前校正后80帧后粗糙度ρ0.680.12PSNR (dB)18.433.7肉眼效果上竖条纹几乎不可见圆盘边缘也没有因为整形而出现明显锯齿。这算是这个思路复现成功的一个信号。4. 复现中的坑与排查实录4.1 步长不匹配导致发散最开始的复现里μ_g取的是1e-3μ_o取的是1e-1。结果30帧后增益图出现剧烈震荡输出的图像一片花斑像是随机噪声被无限放大。排查后发现原因y的数值范围太大e * y产生的结果比e本身大两个数量级。而偏置更新项没有乘y所以偏置正常收敛、增益发散输出图像的动态范围被压低又弹高形成花斑。解决办法很简单把输入图像做归一化到0至1或者把μ_g缩小到1e-4。一般来说增益步长与偏置步长的比例应该大致参考y均值的倒数。如果y均值是100μ_g取μ_o的百分之一左右正好。4.2 固定条纹残留与块边界振铃复现后期图像整体校正得还算干净但每8×8的块边界处会出现轻微马赛克感。问题出在分块DCT整形后直接拼接块与块之间的连续性被破坏。用镜像填充边缘、以及把整形后的参考图像与原始图像做0.7和0.3的加权混合能明显缓解块效应。加权混合的理由也很实在参考图像是频域滤波重建的结果本身丢失了部分细节完全用它做LMS的期望信号会导致输出过于平滑。保留一部分原始成分相当于给参考信号加了一点“原汁原味”防止过平滑。4.3 运动目标鬼影问题这个算法在匀速运动的小目标场景里表现还行但遇到大面积物体突然进入画面时依然会短暂产生残影。原因是物体边缘在频域里也属于高频能量光谱整形统计如果只靠能量高低区分噪声和场景很容易把物体边缘当成可抑制成分。我的处理是增加边缘掩膜参考图像中梯度大的区域直接使用原始像素值不参与LMS更新。具体做法是计算Sobel梯度超过阈值的像素把误差e置零不更新g和o。这个“边缘保护”牺牲了一点收敛速度但大幅提升了鲁棒性。原文献里的“强鲁棒性”应该也包含类似的机制。4.4 粗糙度指标自欺欺人有段时间我把damping调到0.1校正后的ρ降到了0.06看起来“优化”了但图像细节几乎全没了圆盘边缘模糊成一团。这个教训是单看ρ是很容易过平滑的。每次调参都应该用原始图像对比或者同时计算高频细节保留度。在实际工程里有经验的工程师会直接拿校正前后的视频去观察运动目标边缘是否清晰比任何单一指标都可靠。4.5 增益更新不能盲目全员一致对全部像素使用同一个μ_g会导致亮区和暗区的收敛速度不一样。亮区的e乘以y后数值大更新幅度大暗区数值小更新慢。表现是校正过程里画面亮部先稳定暗部还在慢慢跳动。解决方式是使用归一化LMSNLMS的思想把步长除以输入信号的局部能量μ_g μ_0 / (y² ε)其中ε是防止除零的常数。这个改动让每个像元的收敛速度趋于一致整体校正质量更均衡。我复现时做了这个改动效果很直观100帧之后整幅图像的增益图和偏置图都比固定步长版本更加均匀。5. 实测效果与个人心得5.1 多场景测试结果除了仿真序列我还用了一组真实红外视频做测试场景是室内的走廊帧率25fps分辨率320×256。这个算法运行下来没有明显卡顿单帧处理时间在Python环境下大约是35毫秒考虑到DCT和复数更新表现算不错。结果上竖条纹在30帧左右开始减弱80帧后基本消失。走廊里有人走动时人物边缘保持了清晰轮廓没有出现严重的拖影。相比于单纯的高通滤波法这个算法在细节保持上有明显优势相比于均值和方差恒定统计法它对场景动态变化的适应性更强。5.2 参数速查经验给直接上手复现的朋友一份参数速查表基于我测试的0到1归一化数据参数推荐值说明块大小8分块DCT的边长damping0.25DC系数的保留比例低频权重0.3~0.5前三个AC系数的缩放系数μ_g1e-4增益步长μ_o1e-2偏置步长边缘保护阈值0.08Sobel梯度的判断阈值统计窗口30帧频谱能量统计的滑动窗口这些参数的绝对值没有普适性因为不同相机的响应范围差异很大。建议在自己的数据集上做一次网格搜索先固定μ_o调μ_g再固定μ_g调damping。5.3 我个人的一些体会做这个复现花了我好几天中间反复调试的主要问题不是数学公式看不懂而是对“期望信号怎么来”这一点理解不够深。LMS框架本身是死的它只负责沿梯度方向更新权重真正难的是构造一个合理的期望信号。光谱整形统计给的是一个可解释、可控制的构造方式它的优点在于把频域先验和自适应更新结合在了一起。后面我又尝试过小波变换域做整形效果类似但计算量更大。也试过用引导滤波代替频谱整形来产生参考图像边缘保持更好但对条纹噪声的抑制弱一些。如果做工程建议在频谱整形的参考信号上叠加一点引导滤波的结果两者组合起来会更好。最后说一点题外话非均匀性校正这类算法的实际效果极度依赖输入图像的质量。如果前级电路有严重的时域随机噪声LMS更新会被随机分量干扰收敛后参数依然抖动。所以做项目时先统计原始视频的噪声属性再决定步长不要上来就套参数。我是把这个流程固化下来了每次新项目都先跑一段静止场景的方差图心里有底之后才动算法。