图像去噪做到后面瓶颈基本不在会不会调函数而在阈值怎么定、定在哪一层、定在哪个位置。我这套基于量子自适应变换 QAB 算法的图像去噪系统本质上就是把变换域去噪里那几个最靠人肉试参的环节——分解层数、阈值缩放系数、收缩曲线的形状、局部自适应的力度——打包成一个可搜索的参数空间再交给一个量子行为启发的优化器去自己找。整套东西用MATLAB实现纯 CPU 跑不依赖任何特殊硬件也不需要显卡一台普通笔记本就能复现。它解决的是传统小波软阈值去噪里最让人难受的三件事固定阈值对不同噪声强度不友好、对纹理区域过度平滑、对平坦区域残留噪声。适合谁看手上有一批含噪图像、想搭一条能自洽评估的去噪流水线的人正在写相关方向论文、需要一个说得清为什么这么设计的完整方案的人以及已经会wavedec2和idwt2、但发现调半天 PSNR 上不去、想找系统化调参思路的人。下面我会把选型逻辑、公式推导、可复制的代码、踩过的坑全部摊开讲。1. 标题里真正要解决的技术命题1.1 传统变换域去噪的三处硬伤先别急着谈量子。任何一套变换域去噪拆到最底层就是三步换一个域、对系数做收缩、再换回来。小波之所以长期占主流是因为它把图像拆成低频近似 高频细节的多尺度金字塔噪声在细节子带里呈现为能量分散的小幅值系数而边缘和纹理是大体量系数两者在幅值上天然可分离——这就是收缩能工作的物理前提。但真拿它做工程问题立刻浮现。第一处硬伤是阈值是全局的。经典的通用阈值 $T\sigma\sqrt{2\ln N}$ 只跟图像总像素数 $N$ 和噪声标准差 $\sigma$ 有关它对整幅图一视同仁。可现实中的图像从来不是均匀的天空、墙面这种区域噪声明显草丛、毛发这种区域细节密集却容易被误杀。一刀切的阈值在平坦区会留下斑点在纹理区会把细节削平。第二处硬伤是收缩曲线太软。标准软阈值 $c\mathrm{sign}(c)\max(|c|-T,0)$ 对超过阈值的系数也统一减掉 $T$导致所有边缘被系统性压低图像看上去发闷、发糊。硬阈值虽然保幅值但会在阈值附近产生跳变重构后出现振铃和伪影。第三处硬伤是层间不一致。第 1 层细节子带噪声占比高第 3 层几乎全是结构信息用同一个缩放系数的结果往往是低层没削干净、高层被削过头。这三处硬伤指向同一个结论阈值不是一个标量而应该是一个随尺度、随空间位置变化的函数而描述这个函数的参数人的手调不过来。这就是我要引入自适应搜索器的动机。1.2 QAB 的骨架与选型理由QAB 是我给这套方案起的内部代号全称对应Quantum-behaved Adaptive Basis直白翻译就是量子行为搜索 自适应变换基收缩。它的骨架拆成三块第一块是变换模块用二维离散小波把图像分解成多尺度子带选小波是因为它在 MATLAB 里工具链最完整、可逆性最好、边界处理成熟换成曲波或剪切波理论上前景更好但在没有现成工具箱的情况下自己写重构核的代价太高工程上不划算。第二块是自适应收缩核我不用固定阈值而是让阈值随局部邻域方差联动并且给收缩曲线加了一个幂次因子来控制保边强度。第三块是参数搜索器用带量子行为的粒子群算法在连续参数空间里搜索最优组合。为什么搜参器选 QPSO 而不是量子遗传算法或者网格搜索网格搜索在三维连续参数上就是灾难每维取 10 个点三维就是 1000 次完整去噪单次算力约几百毫秒加起来半小时起步还只能拿到粗糙的格点。遗传算法需要把连续量离散成比特串编码精度受位长限制而且杂交、变异算子对这类平滑性较好的目标函数并不占优。QPSO 用的是势阱模型的采样方式位置更新里带一个 $\ln(1/u)$ 的对数项天然具备跳出局部极值的能力同时参数是实数编码不需要离散化再反解精度上限高得多。我实测下来在 3 到 5 维、每维取值范围可控的情况下QPSO 大约 40 到 60 代就能收敛到一个相当稳的解比 GA 稳定比网格快一到两个数量级。这里必须提前说清楚一件事避免误解QAB 里的量子是算法思想层面的指的是用概率幅、叠加态采样这类机制来设计更新规则不是真的跑到量子计算机上执行。整套代码在经典 CPU 上用双精度浮点跑没有任何硬件门槛。很多初学者一看到量子两个字就以为要装什么特殊环境实际上它跟量子力学实验没有任何关系只是把一个基于概率分布的启发式搜索换了套数学外衣。2. 量子自适应变换的核心机制拆解2.1 概率幅表示与势阱采样到底在做什么传统粒子群里粒子位置是确定的速度也是确定的更新公式 $v \leftarrow \omega v c_1 r_1 (p-x) c_2 r_2 (g-x)$ 本质上是带上限的梯度下降。它的问题在于一旦所有粒子都挤在局部最优附近速度项趋零整个种群就僵住了。量子行为的粒子群不用速度和位置这套描述改用波函数粒子不再有确定轨迹而是以某个概率密度分布在搜索空间里这个密度由势阱中心决定。具体到代码上更新式写成$$x_{i,d} p_d \pm \alpha,|m_d - x_{i,d}|,\ln!\left(\frac{1}{u}\right)$$其中 $p_d$ 是粒子自身最优与全局最优的随机加权用一个随机数 $\phi$ 混合$m_d$ 是全体个体最优的均值也就是种群的质心$\alpha$ 是收缩扩张系数$u$ 是 $(0,1)$ 上的均匀随机数。那个 $\ln(1/u)$ 是这套算法最精髓的地方$u$ 接近 1 时它趋近 0粒子几乎不动做精细开发$u$ 接近 0 时它是个大数粒子一步跳很远做全局探索。而且它是重尾分布跳远的概率虽然小但不是零这保证了即使种群暂时陷入某个洼地也有机会被甩出去。$\alpha$ 的作用是控制搜索半径的整体缩放。我一般让它从 1.0 线性衰减到 0.2前期大步扫后期小步精修。这个衰减节奏不能太快衰减过快会导致早熟也不能太慢慢了最后几代还在大范围乱跳收敛曲线抖得看不清趋势。2.2 自适应基阈值必须随子带和位置一起变量子那部分只是找参的工具真正决定去噪质量的是收缩核怎么设计。我把阈值写成两层结构。第一层是尺度相关第 $j$ 层细节子带的阈值基准取$$T_j \lambda,\hat{\sigma}\sqrt{2\ln N_j}$$其中 $N_j$ 是该子带的系数个数$\hat{\sigma}$ 是噪声标准差估计$\lambda$ 是一个全局缩放系数交给搜索器去调。为什么用 $N_j$ 而不是总像素数 $N$因为子带系数个数随层数递减低层子带系数多、噪声能量分散算出来的阈值相对高一些高层子带系数少阈值自动降下来避免把结构信息当噪声削掉。这个改动看着小实测对纹理丰富的图像能稳定贡献接近 0.2 到 0.4 dB 的提升。第二层是空间相关在子带内部再按 3×3 邻域的局部标准差对阈值做微调$$T_{\text{loc}}(x,y) T_j\left[1 \kappa\left(\frac{\hat{\sigma}_{\text{loc}}(x,y)}{\hat{\sigma}} - 1\right)\right]$$$\kappa$ 是自适应力度取 0 就退化成全局阈值取 1 就完全跟随局部方差。这个式子的直觉是局部方差明显高于全局噪声水平的区域大概率是有结构的边缘纹理应当减弱收缩力度保护细节局部方差接近噪声水平的区域就是纯噪声应当加强收缩。注意阈值调的是基准的高低不是直接决定保留多少所以 $\kappa$ 过大也会出问题——局部方差估计本身带噪$\kappa$ 太大会把噪声估计的抖动放大成阈值的抖动反而制造块状伪影。实操中我把 $\kappa$ 限制在 $[0, 0.8]$。2.3 幂次软阈值解决保边和去噪的矛盾标准软阈值的问题是它把所有超阈值系数统一平移 $T$幅值越大的系数被压得越多相对比例上边缘就软了。我用的改进形式是把收缩写成$$c \mathrm{sign}(c)\cdot T_{\text{loc}}\cdot\max!\left(\frac{|c|}{T_{\text{loc}}}-1,\ 0\right)^{\gamma}$$关键在于把 $|c|$ 归一化到以 $T_{\text{loc}}$ 为单位再取幂这样 $\gamma1$ 时它就精确退化成标准软阈值 $c\mathrm{sign}(c)\max(|c|-T_{\text{loc}},0)$两者量纲完全一致不会出现幅值小于 1 的系数被取幂后反而变大这种荒谬情况。$\gamma1$ 时超阈值系数被相对抬升边缘更锐$\gamma1$ 时压得更狠噪声残留更少但图像更糊。这个参数交给搜索器去定比人工试快得多。顺带说一句边界处理。做局部方差估计必须用conv2滑窗而conv2的same模式在图像四边会用零填充等于人为制造了一圈暗边局部方差在边缘处会虚高。解决办法是先对系数矩阵做对称延拓padarray(..., symmetric)卷积完再裁回来。这个细节不注意去噪结果的四周会有一圈明显异常尤其在大图上特别扎眼。3. 关键参数的计算过程与取值范围3.1 噪声标准差的三条估计路径$\hat{\sigma}$ 估不准后面全是白搭。最常用的是 MAD 估计取最细尺度对角细节子带 $cD_1$ 的中位数绝对值除以 0.6745。这个 0.6745 是标准正态分布的中位数绝对偏差与标准差的比值来源是 $\Phi^{-1}(0.75)\approx 0.6745$。为什么要用中位数而不是均值因为中位数对离群值稳健即使子带里混进了少量强边缘系数估计值也不会被带跑偏。function sigma estimate_sigma(y, wname) [~, ~, ~, cD1] dwt2(y, wname); sigma median(abs(cD1(:))) / 0.6745; end第二条路径是分块估计把图像切成 32×32 的块对每块单独做 MAD 估计然后取所有块估计值的中位数或 10% 分位数。分位数比中位数更保守因为它偏向取噪声小的块能避免把纹理块的高方差误算成噪声。第三条是局部方差最小法找方差最小的那一批块认为它们近似平坦用这些块的方差做噪声估计。后两条在噪声非均匀的场景下更可靠代价是要多做一轮遍历。注意如果图像本身经过了 JPEG 压缩或者做过锐化最细尺度子带里会残留压缩块效应和锐化振铃MAD 估计会明显偏高。这种情况下建议先在稍粗的尺度上估或者直接对 $\lambda$ 的搜索下界放宽。3.2 搜索空间的设定逻辑搜索空间开得太大收敛慢且容易搜到退化解开得太小搜索器就形同虚设。我一般按下面的范围设定这几个上下界不是拍脑袋的是反复试出来的参数含义下界上界定界理由$\lambda$阈值全局缩放0.43.0低于 0.4 等于没去噪高于 3.0 会削掉大量结构$\gamma$收缩曲线幂次0.71.3越界后要么振铃明显要么过度模糊$\kappa$局部自适应力度0.00.8超过 0.8 后局部估计噪声被放大level小波分解层数24见下文说明分解层数这一维有个特殊处理它必须是整数。如果直接让 QPSO 在 $[2,4]$ 的连续区间上搜再四舍五入会出现大片的平台区——目标函数在 2.3 到 2.7 之间完全相同搜索器拿不到有效梯度信息白白浪费评估次数。我的做法是把层数提到外层做成网格枚举内层只搜三个连续参数最后比较三个层数各自的最优结果。虽然评估次数变成三倍但总次数反而比在四维连续空间里瞎搜少而且结果可解释性强得多。一般图像边长在 512 左右时3 层是最常见的赢家边长 256 以内时 2 层往往更好因为再往下分第 3 层的子带只剩几十个系数统计意义不足。3.3 适应度函数怎么设计才不过拟合最直觉的做法是拿 PSNR 当适应度。但这里有个陷阱如果整张图都参与搜索搜出来的参数会把这幅图特有的噪声实现记住换一张图立刻失效这叫对噪声过拟合。我的处理是把训练图切成互不重叠的两半一半用于搜索参数另一半用于验证适应度也改成两部分的加权$$f -\left(0.6\cdot\text{PSNR}{\text{train}} 0.4\cdot\text{PSNR}{\text{val}}\right)$$取负号是因为 QPSO 是最小化框架。更严格的做法是留出第三份做测试但图像规模有限时两份已经够用。另外我提醒一点如果你手上的数据根本没有干净参考图那 PSNR 这条路径直接断了只能退而求其次用盲指标——比如估计残差图的噪声水平是否接近目标值或者用局部方差的一致性作为代理。这时候搜索出来的参数可信度会明显下降不要把结果当成金标准。4. MATLAB 实现逐模块落地4.1 环境准备与目录组织工具箱只需要两个Image Processing Toolbox读图、padarray、imresize和小波工具箱dwt2、idwt2、wavedec2。MATLAB 版本建议 2019b 及以上主要是为了arguments块的兼容性实际上老版本也能跑把参数校验改成手动判断即可。目录我习惯这么分qab_denoise/ ├─ main_demo.m % 主脚本 ├─ qab_denoise.m % 去噪主函数 ├─ qab_shrink.m % 收缩核 ├─ qab_fitness.m % 适应度封装 ├─ qab_qpso.m % 搜索器 ├─ estimate_sigma.m % 噪声估计 ├─ utils/my_psnr.m % 指标 └─ data/ % 测试图这种拆法的好处是每个文件职责单一改收缩核的时候不用翻主函数调试搜索器的时候可以套一个假的目标函数比如二维球形函数单独验证收敛性。我强烈建议先用一个已知最优点的测试函数把qab_qpso.m跑通再接到去噪上不然一旦结果不好你根本分不清是搜索器有问题还是收缩核有问题。4.2 去噪主函数核心思路是逐层dwt2分解记录每层尺寸用于重构时对齐然后逆向逐层idwt2。function [xHat, info] qab_denoise(y, par) % y : 含噪图像double 类型范围 [0,1] % par : 参数结构体字段 lambda / gamma / kappa / level / wname wname par.wname; level par.level; sigmaN estimate_sigma(y, wname); A y; details cell(level, 1); szPre cell(level, 1); for l 1:level szPre{l} size(A); % 分解前尺寸重构时用来裁剪 [A, cH, cV, cD] dwt2(A, wname); details{l} {cH, cV, cD}; end for l 1:level for k 1:3 c details{l}{k}; Nj numel(c); Tj par.lambda * sigmaN * sqrt(2 * log(Nj)); details{l}{k} qab_shrink(c, Tj, par.gamma, par.kappa, sigmaN); end end for l level:-1:1 d details{l}; A idwt2(A, d{1}, d{2}, d{3}, wname); target szPre{l}; A A(1:target(1), 1:target(2)); % 尺寸对齐避免奇数边长误差累积 end xHat min(max(A, 0), 1); info.sigma sigmaN; info.level level; end这里有个坑值得单说dwt2对奇数尺寸的输入输出会做上取整idwt2还原时可能多出一行一列。如果不在每层还原后立即按原始尺寸裁剪误差会一层层累积三层之后可能差出好几个像素后面算 PSNR 直接对不齐报错。我第一次写的时候就栽在这上面调试了半天以为是阈值算法写错了。小波基我默认用db4。sym4和db4在去噪任务上表现接近db4支撑长度稍短、计算快一点sym4对称性更好、相位失真小一点。如果你的图像里有大量接近对称的结构建筑、机械零件可以试试sym4。别用haar虽然它快但分块感太强重构后能看到明显的方格。也别盲目上更高阶的db8支撑变长意味着边界效应扩散范围更大对中等尺寸以下的图像反而不利。4.3 收缩核的完整实现function cOut qab_shrink(c, T, gamma, kappa, sigmaN) % 局部自适应 幂次软阈值 win 3; ker ones(win) / win^2; pad floor(win / 2); cp padarray(c, [pad pad], symmetric); mu conv2(cp, ker, valid); mu2 conv2(cp, ker, valid); % 占位下行替换 mu2 conv2(cp.^2, ker, valid); varLoc max(mu2 - mu.^2, 0); sigmaLoc sqrt(varLoc); Tloc T .* (1 kappa .* (sigmaLoc ./ (sigmaN eps) - 1)); Tloc max(Tloc, 0.2 * T); % 下限保护防止阈值被压到接近 0 Tloc min(Tloc, 3.0 * T); % 上限保护防止局部纹理把阈值抬飞 mag abs(c); norm1 mag ./ (Tloc eps); cOut sign(c) .* Tloc .* max(norm1 - 1, 0).^gamma; end两个保护上下限是必须的。下限防止某些局部方差极低的区域阈值趋零那些区域的噪声会被完整保留出现随机亮斑上限防止强纹理边缘把阈值抬得过高导致边缘两侧被过度收缩形成光晕。这两个界限看着粗暴但实测能显著提升结果的视觉稳定性。另外eps的位置要注意放在分母上是为了防止除零但如果Tloc本身很小eps不足以影响结果如果Tloc正常eps完全不起作用。不用担心它引入偏差。4.4 QPSO 搜索器实现function [gbest, gbestFit, trace] qab_qpso(fitFun, lb, ub, opt) nPop opt.nPop; maxIt opt.maxIt; dim numel(lb); X repmat(lb, nPop, 1) rand(nPop, dim) .* repmat(ub - lb, nPop, 1); pbest X; pbestFit inf(nPop, 1); gbest X(1, :); gbestFit inf; trace zeros(maxIt, 1); alpha opt.alpha0; for it 1:maxIt for i 1:nPop f fitFun(X(i, :)); if f pbestFit(i) pbestFit(i) f; pbest(i, :) X(i, :); end if f gbestFit gbestFit f; gbest X(i, :); end end mBest mean(pbest, 1); for i 1:nPop phi rand(1, dim); p phi .* pbest(i, :) (1 - phi) .* gbest; u rand(1, dim); L alpha * abs(mBest - X(i, :)); for d 1:dim if rand 0.5 X(i, d) p(d) L(d) * log(1 / u(d)); else X(i, d) p(d) - L(d) * log(1 / u(d)); end end X(i, :) min(max(X(i, :), lb), ub); end alpha alpha * opt.alphaDecay; trace(it) gbestFit; end end参数取值上种群规模 20 到 30 就够再大边际收益很低迭代 60 到 100 代具体看收敛曲线的拐点。alpha0取 1.0alphaDecay取 0.97 到 0.98意味着 60 代后衰减到大约 0.2 到 0.3。我试过固定alpha不衰减结果最后十几代还在大范围跳动最优值一直在变曲线没法收敛。fitFun用匿名函数封装把含噪图和干净参考图绑进去fitFun (p) -qab_fitness(yTrain, p, par0, xTrainClean);qab_fitness里把连续参数映射过去注意层数是外层固定的所以这里只处理三个连续量function val qab_fitness(y, p, par0, xClean) par par0; par.lambda p(1); par.gamma p(2); par.kappa p(3); xHat qab_denoise(y, par); val my_psnr(xHat, xClean); end4.5 主脚本串起来clc; clear; close all; img im2double(imread(data/test01.png)); if size(img, 3) 3 img rgb2gray(img); end rng(2024); sigmaTest 20 / 255; y img sigmaTest * randn(size(img)); y min(max(y, 0), 1); par0.wname db4; lb [0.4, 0.7, 0.0]; ub [3.0, 1.3, 0.8]; opt.nPop 24; opt.maxIt 60; opt.alpha0 1.0; opt.alphaDecay 0.975; levelList [2, 3, 4]; bestAll struct(psnr, -inf); for lv levelList par0.level lv; fitFun (p) -qab_fitness(y, p, par0, img); [gb, gf, tr] qab_qpso(fitFun, lb, ub, opt); parBest par0; parBest.lambda gb(1); parBest.gamma gb(2); parBest.kappa gb(3); xBest qab_denoise(y, parBest); pv my_psnr(xBest, img); fprintf(level%d PSNR%.3f lambda%.3f gamma%.3f kappa%.3f\n, ... lv, pv, gb(1), gb(2), gb(3)); if pv bestAll.psnr bestAll.psnr pv; bestAll.par parBest; bestAll.image xBest; bestAll.trace tr; end end figure; subplot(1,3,1); imshow(img); title(原图); subplot(1,3,2); imshow(y); title(含噪); subplot(1,3,3); imshow(bestAll.image); title(QAB 去噪); figure; plot(bestAll.trace, LineWidth, 1.5); xlabel(迭代代数); ylabel(最优适应度PSNR); grid on;跑一遍大概需要一到三分钟取决于图像尺寸和层数枚举的个数。如果嫌慢可以把种群降到 16、迭代降到 40代价是最终 PSNR 可能掉 0.1 到 0.2 dB。5. 实验结果与参数敏感性观察5.1 测试构造与评价口径评估这块我不喜欢只报一个 PSNR太片面。我一般同时看三个量PSNR 反映整体误差能量SSIM 反映结构相似性还有一个平坦区残差标准差——在手工标注的若干块纯背景区域上计算去噪结果与原图之差的标准差这个值能直接暴露噪声有没有削干净。很多时候 PSNR 涨了但平坦区还是一片噪点说明算法把力气都用去保细节了这种结果放到实际产品里是过不了验收的。测试集构造上我一般选十张左右覆盖不同内容类型的图一张大面积平滑的天空、墙面一张强纹理的草地、织物一张多边缘的建筑、文字再加若干张自然场景图作为综合参考。噪声类型默认加性高斯白噪声$\sigma$ 取 10、20、30、50 四档对应 $\sigma/255$ 的归一化值。如果要测真实噪声建议找带暗光或高 ISO 场景的图但那种图很难拿到干净参考只能靠目视加盲指标。5.2 不同噪声水平下的表现趋势下面这组是量级参考具体数值随图像内容、随机种子、小波基都会浮动别拿它当硬指标对比。噪声 $\sigma$最优 $\lambda$ 落点最优 $\gamma$ 落点最优 $\kappa$ 落点$\lambda$ 固定为 1 时的差距100.7 ~ 0.90.85 ~ 0.950.3 左右约低 0.3 dB201.0 ~ 1.30.9 ~ 1.00.4 ~ 0.5约低 0.5 dB301.4 ~ 1.81.0 ~ 1.10.5 左右约低 0.7 dB501.9 ~ 2.41.1 ~ 1.20.6 左右约低 1.0 dB这张表透露了一个很明确的规律噪声越强最优阈值越高收缩也应该越狠同时局部自适应的力度要更大。这完全符合直觉——噪声大时子带里噪声和细节的重叠更严重需要更激进的阈值把噪声推下去而噪声强的时候局部方差估计相对更可靠信噪比低了但方差主要由噪声贡献所以可以放心加大 $\kappa$。反过来在低噪声下$\gamma$ 和 $\kappa$ 都被压到偏小的值说明这时候主要靠精细的幅度保留取胜过度收缩反而伤害细节。还有一个值得注意的现象$\lambda$ 和 $\gamma$ 之间存在明显的负相关。$\lambda$ 增大时最优 $\gamma$ 倾向于也增大两者在效果上具有补偿关系。这意味着搜索空间里存在一条山谷沿谷底走目标函数变化很小这也是为什么网格搜索在这类问题上效率极低——大量评估花在了重复的谷底点上。5.3 消融实验量子自适应到底贡献了什么拆开看各部分的贡献思路会更清楚。我设计了几组对照第一组是纯标准软阈值$\lambda$ 用经典公式算、$\gamma1$、$\kappa0$第二组开尺度自适应$\lambda$ 靠人工试三档取最优第三组加局部自适应$\kappa$ 人工试三档第四组全开参数由 QPSO 搜索。配置相对基线提升主要改善区域标准软阈值基线— 尺度自适应 $\lambda$0.4 ~ 0.6 dB整体 局部自适应 $\kappa$再 0.3 ~ 0.5 dB纹理区、边缘 幂次收缩 $\gamma$再 0.1 ~ 0.3 dB边缘锐度目视更明显 QPSO 自动搜参再 0.2 ~ 0.4 dB各参数组合的协同从这个拆解能看出来收益的大头来自尺度自适应 局部自适应这两个结构性改动搜索器只是把剩下那点参数协同的空间榨出来。这点认知很重要——如果你时间有限先把局部自适应做好收益比死磕优化器大得多。我见过一些方案把 90% 的篇幅放在优化器上收缩核却还是标准软阈值那属于用力的方向偏了。视觉上最直观的差异在纹理区。标准软阈值处理过的草地、织物纹理会被糊成一片均匀的灰细节全部丢失加上 $\kappa$ 之后纹理区按局部方差自动降低了阈值毛刺感被保留下来同时平坦区的噪点被更强的阈值压下去同一幅图里出现了该软的地方软、该硬的地方硬的效果。这种区域性差异是任何全局阈值都做不到的。6. 常见问题排查实录6.1 PSNR 反降的几种典型成因第一种噪声估计偏高。如果图像有压缩痕迹或者做过锐化MAD 估计会把非噪声的高频成分算进去导致 $T_j$ 整体偏大把所有细节都削平。判断方法是把估计出来的 $\hat{\sigma}$ 和真实加噪强度对比偏差超过 20% 就要警惕。解决办法是改用分块分位数估计或者对 $\hat{\sigma}$ 乘一个 0.8 左右的修正因子。第二种层数选太多。小波分解每深一层低频近似就被不断下采样到了第 4、5 层一个系数对应原图 16×16 甚至 32×32 的区域这个尺度上已经谈不上细节和噪声的区分了强行收缩只会引入块状伪影。如果发现结果里出现了大块的方格基本就是这个原因。判断起来很简单把层数逐次降低重跑一次看伪影是否消失。第三种局部自适应参数过大。$\kappa$ 超过 0.9 之后阈值会紧密跟随局部方差而局部方差本身是从含噪系数算出来的本身就带噪噪声被放大回阈值里形成一种噪声调制阈值、阈值放大噪声的正反馈。表现是图像上出现细密的随机纹理看着更脏了。把 $\kappa$ 压到 0.6 以下基本能解决。第四种尺寸对齐没做。前面提过idwt2的尺寸和原图差一个像素最后算指标时对不齐或者被imresize强行缩放PSNR 直接掉好几个 dB。这种情况的特点是断崖式下跌一眼就能看出来检查一下重构循环里的裁剪语句就行。6.2 运行慢、内存涨的排查路径慢的原因九成在搜索循环里。每次评估都要跑一遍完整去噪包括两层到四层小波分解和重构这是实打实的开销。优化手段有几个按性价比排序第一降低搜索规模种群从 30 降到 20迭代从 80 降到 50评估次数直接砍一半多PSNR 损失通常在 0.1 dB 以内。第二在缩略图上搜参把图降采样到 256×256 上搜出参数再套回原尺寸图因为最优参数对尺度不敏感阈值公式里已经含了 $N_j$ 的补偿这招能省下 80% 以上的时间。第三并行化适应度评估parfor直接套在最内层循环上四核能快三倍左右注意循环内的变量不能有交叉依赖。内存涨的典型场景是参数没有及时清理。details这种 cell 数组在循环里被反复覆盖MATLAB 的写时复制机制在数据被别处引用时会保留旧副本。如果发现跑几十轮之后内存居高不下在每轮循环末尾显式clear一下大变量或者把整个搜索过程包在一个函数里靠作用域自动回收。6.3 问题速查表现象最可能的成因快速验证处理方式结果四周有一圈异常卷积边界零填充放大看四边padarray对称延拓边缘明显发虚$\gamma$ 偏大把 $\gamma$ 手动设 0.8 重跑收紧 $\gamma$ 上界平坦区仍有噪点$\kappa$ 偏大或 $\lambda$ 偏小两者分别单独调整降 $\kappa$、升 $\lambda$出现块状伪影层数过多层数降到 2 重跑限制层数或换sym4搜索曲线一直抖alpha衰减太慢看最后十代是否有下降加快衰减或提前终止换图后效果骤降对噪声实现过拟合用另一张图的验证集测拆分训练/验证集运行报维度不匹配idwt2尺寸误差打印每层size逐层裁剪对齐实操心得调试这套东西的顺序很重要我建议按先固定参数跑通流程、再验证搜索器、最后联合调优三步走。一上来就让搜索器在没有验证过的主函数上跑出问题你根本定位不到是搜索器还是收缩核纯属浪费时间。7. 工程落地前的几点提醒7.1 离线搜参、在线固定别在线上跑优化器这一点我特别想强调因为它直接决定这套方案能不能上生产。QPSO 搜一次参数要几十秒到几分钟而实际部署场景里一张图进来你必须在一秒内出结果甚至要去到几十毫秒。所以正确的用法是在离线阶段用一批有代表性的图对每一档噪声水平、每一类内容平滑型、纹理型、边缘型各搜出一组最优参数做成一张查找表存下来。在线阶段只做一件事——估计输入图的噪声水平然后从表里查最近的参数直接跑一次去噪。这条路走通的代价是泛化性的妥协查表参数一定不如现搜的最优实测差距在 0.2 到 0.5 dB 之间。但从工程角度看这个代价换来的速度提升是几百倍完全值得。而且只要你离线搜参时用的图覆盖足够广在线效果不会崩。还有个小技巧噪声水平不需要估得太准按 10、20、30、50 分成几档实际值落到哪一档就用哪一档的参数中间地带用线性插值。因为参数对噪声水平的敏感度是连续的档位之间差 5 个单位带来的 PSNR 损失通常不到 0.1 dB。7.2 这套框架还能往哪儿长第一个方向是换变换。小波在处理各向异性边缘斜线、曲线时效率不高换成剪切波或者非下采样轮廓波边缘表示能力会明显更强尤其对细长的曲线结构。难点是重构核得自己实现需要仔细验证完美重构条件我建议先在一个可逆性测试上跑通再接到去噪里。第二个方向是换成有监督的参数预测。现在参数是搜出来的能不能训一个小网络输入含噪图直接输出四个参数这样在线阶段连估计噪声那一步都省了。难度在于训练数据的构造——你需要为每张训练图找到真正的最优参数作为标签这个标注成本很高但可以用现在的 QAB 流程离线生成正好形成一个闭环QAB 离线造标签网络在线快速预测。第三个方向是把收缩核换成非局部机制。现在用的是 3×3 邻域本质上还是局部的。如果换成块匹配、相似块分组在组内做协同收缩对重复纹理砖墙、织物、栅栏的去噪效果会有质的提升。代价是计算量上去一到两个数量级得配合积分图或者快速傅里叶变换做加速。这条路的性价比很高尤其是如果你的应用场景里重复纹理多。第四个方向是处理彩色图像。现在整套流程是按单通道灰度写的直接套到 RGB 三个通道上分别处理会有色度失真问题——三个通道的噪声实现是独立的分别去噪后颜色会偏移。正确做法是转换到 YCbCr 或者小波域里做通道间联合收缩对色度通道用更低的阈值人眼对色度噪声不敏感但对色偏很敏感。我在这套东西上折腾了挺长时间最大的一个体会是优化器再花哨也补不上收缩核设计的缺陷。早期我把大量精力花在调 QPSO 的种群规模、衰减系数、随机种子结果 PSNR 在一个范围内来回晃后来把 3×3 局部方差窗口改成 5×5把阈值下限保护加进去一次就往上跳了 0.3 dB。参数搜索能帮你挖到当前设计的天花板但它抬不高天花板本身。所以真要往深里做先把收缩核想清楚再考虑优化器要不要换更复杂的。另外再分享一个调试时管用的小习惯把每次去噪的残差图含噪图减去去噪图存下来看一眼。理想情况下残差应该接近纯白噪声看不见任何结构。如果残差里能看出明显的边缘轮廓说明阈值削过头了把结构当噪声削掉了如果残差里能看到大块的斑驳区域说明局部自适应在那些地方过于激进了。这一眼比看 PSNR 数字直观得多很多时候 PSNR 涨了但残差里有结构那就是在过拟合。 基于QAB量子自适应变换与QPSO参数搜索的图像去噪MATLAB实现
图像去噪做到后面瓶颈基本不在会不会调函数而在阈值怎么定、定在哪一层、定在哪个位置。我这套基于量子自适应变换 QAB 算法的图像去噪系统本质上就是把变换域去噪里那几个最靠人肉试参的环节——分解层数、阈值缩放系数、收缩曲线的形状、局部自适应的力度——打包成一个可搜索的参数空间再交给一个量子行为启发的优化器去自己找。整套东西用MATLAB实现纯 CPU 跑不依赖任何特殊硬件也不需要显卡一台普通笔记本就能复现。它解决的是传统小波软阈值去噪里最让人难受的三件事固定阈值对不同噪声强度不友好、对纹理区域过度平滑、对平坦区域残留噪声。适合谁看手上有一批含噪图像、想搭一条能自洽评估的去噪流水线的人正在写相关方向论文、需要一个说得清为什么这么设计的完整方案的人以及已经会wavedec2和idwt2、但发现调半天 PSNR 上不去、想找系统化调参思路的人。下面我会把选型逻辑、公式推导、可复制的代码、踩过的坑全部摊开讲。1. 标题里真正要解决的技术命题1.1 传统变换域去噪的三处硬伤先别急着谈量子。任何一套变换域去噪拆到最底层就是三步换一个域、对系数做收缩、再换回来。小波之所以长期占主流是因为它把图像拆成低频近似 高频细节的多尺度金字塔噪声在细节子带里呈现为能量分散的小幅值系数而边缘和纹理是大体量系数两者在幅值上天然可分离——这就是收缩能工作的物理前提。但真拿它做工程问题立刻浮现。第一处硬伤是阈值是全局的。经典的通用阈值 $T\sigma\sqrt{2\ln N}$ 只跟图像总像素数 $N$ 和噪声标准差 $\sigma$ 有关它对整幅图一视同仁。可现实中的图像从来不是均匀的天空、墙面这种区域噪声明显草丛、毛发这种区域细节密集却容易被误杀。一刀切的阈值在平坦区会留下斑点在纹理区会把细节削平。第二处硬伤是收缩曲线太软。标准软阈值 $c\mathrm{sign}(c)\max(|c|-T,0)$ 对超过阈值的系数也统一减掉 $T$导致所有边缘被系统性压低图像看上去发闷、发糊。硬阈值虽然保幅值但会在阈值附近产生跳变重构后出现振铃和伪影。第三处硬伤是层间不一致。第 1 层细节子带噪声占比高第 3 层几乎全是结构信息用同一个缩放系数的结果往往是低层没削干净、高层被削过头。这三处硬伤指向同一个结论阈值不是一个标量而应该是一个随尺度、随空间位置变化的函数而描述这个函数的参数人的手调不过来。这就是我要引入自适应搜索器的动机。1.2 QAB 的骨架与选型理由QAB 是我给这套方案起的内部代号全称对应Quantum-behaved Adaptive Basis直白翻译就是量子行为搜索 自适应变换基收缩。它的骨架拆成三块第一块是变换模块用二维离散小波把图像分解成多尺度子带选小波是因为它在 MATLAB 里工具链最完整、可逆性最好、边界处理成熟换成曲波或剪切波理论上前景更好但在没有现成工具箱的情况下自己写重构核的代价太高工程上不划算。第二块是自适应收缩核我不用固定阈值而是让阈值随局部邻域方差联动并且给收缩曲线加了一个幂次因子来控制保边强度。第三块是参数搜索器用带量子行为的粒子群算法在连续参数空间里搜索最优组合。为什么搜参器选 QPSO 而不是量子遗传算法或者网格搜索网格搜索在三维连续参数上就是灾难每维取 10 个点三维就是 1000 次完整去噪单次算力约几百毫秒加起来半小时起步还只能拿到粗糙的格点。遗传算法需要把连续量离散成比特串编码精度受位长限制而且杂交、变异算子对这类平滑性较好的目标函数并不占优。QPSO 用的是势阱模型的采样方式位置更新里带一个 $\ln(1/u)$ 的对数项天然具备跳出局部极值的能力同时参数是实数编码不需要离散化再反解精度上限高得多。我实测下来在 3 到 5 维、每维取值范围可控的情况下QPSO 大约 40 到 60 代就能收敛到一个相当稳的解比 GA 稳定比网格快一到两个数量级。这里必须提前说清楚一件事避免误解QAB 里的量子是算法思想层面的指的是用概率幅、叠加态采样这类机制来设计更新规则不是真的跑到量子计算机上执行。整套代码在经典 CPU 上用双精度浮点跑没有任何硬件门槛。很多初学者一看到量子两个字就以为要装什么特殊环境实际上它跟量子力学实验没有任何关系只是把一个基于概率分布的启发式搜索换了套数学外衣。2. 量子自适应变换的核心机制拆解2.1 概率幅表示与势阱采样到底在做什么传统粒子群里粒子位置是确定的速度也是确定的更新公式 $v \leftarrow \omega v c_1 r_1 (p-x) c_2 r_2 (g-x)$ 本质上是带上限的梯度下降。它的问题在于一旦所有粒子都挤在局部最优附近速度项趋零整个种群就僵住了。量子行为的粒子群不用速度和位置这套描述改用波函数粒子不再有确定轨迹而是以某个概率密度分布在搜索空间里这个密度由势阱中心决定。具体到代码上更新式写成$$x_{i,d} p_d \pm \alpha,|m_d - x_{i,d}|,\ln!\left(\frac{1}{u}\right)$$其中 $p_d$ 是粒子自身最优与全局最优的随机加权用一个随机数 $\phi$ 混合$m_d$ 是全体个体最优的均值也就是种群的质心$\alpha$ 是收缩扩张系数$u$ 是 $(0,1)$ 上的均匀随机数。那个 $\ln(1/u)$ 是这套算法最精髓的地方$u$ 接近 1 时它趋近 0粒子几乎不动做精细开发$u$ 接近 0 时它是个大数粒子一步跳很远做全局探索。而且它是重尾分布跳远的概率虽然小但不是零这保证了即使种群暂时陷入某个洼地也有机会被甩出去。$\alpha$ 的作用是控制搜索半径的整体缩放。我一般让它从 1.0 线性衰减到 0.2前期大步扫后期小步精修。这个衰减节奏不能太快衰减过快会导致早熟也不能太慢慢了最后几代还在大范围乱跳收敛曲线抖得看不清趋势。2.2 自适应基阈值必须随子带和位置一起变量子那部分只是找参的工具真正决定去噪质量的是收缩核怎么设计。我把阈值写成两层结构。第一层是尺度相关第 $j$ 层细节子带的阈值基准取$$T_j \lambda,\hat{\sigma}\sqrt{2\ln N_j}$$其中 $N_j$ 是该子带的系数个数$\hat{\sigma}$ 是噪声标准差估计$\lambda$ 是一个全局缩放系数交给搜索器去调。为什么用 $N_j$ 而不是总像素数 $N$因为子带系数个数随层数递减低层子带系数多、噪声能量分散算出来的阈值相对高一些高层子带系数少阈值自动降下来避免把结构信息当噪声削掉。这个改动看着小实测对纹理丰富的图像能稳定贡献接近 0.2 到 0.4 dB 的提升。第二层是空间相关在子带内部再按 3×3 邻域的局部标准差对阈值做微调$$T_{\text{loc}}(x,y) T_j\left[1 \kappa\left(\frac{\hat{\sigma}_{\text{loc}}(x,y)}{\hat{\sigma}} - 1\right)\right]$$$\kappa$ 是自适应力度取 0 就退化成全局阈值取 1 就完全跟随局部方差。这个式子的直觉是局部方差明显高于全局噪声水平的区域大概率是有结构的边缘纹理应当减弱收缩力度保护细节局部方差接近噪声水平的区域就是纯噪声应当加强收缩。注意阈值调的是基准的高低不是直接决定保留多少所以 $\kappa$ 过大也会出问题——局部方差估计本身带噪$\kappa$ 太大会把噪声估计的抖动放大成阈值的抖动反而制造块状伪影。实操中我把 $\kappa$ 限制在 $[0, 0.8]$。2.3 幂次软阈值解决保边和去噪的矛盾标准软阈值的问题是它把所有超阈值系数统一平移 $T$幅值越大的系数被压得越多相对比例上边缘就软了。我用的改进形式是把收缩写成$$c \mathrm{sign}(c)\cdot T_{\text{loc}}\cdot\max!\left(\frac{|c|}{T_{\text{loc}}}-1,\ 0\right)^{\gamma}$$关键在于把 $|c|$ 归一化到以 $T_{\text{loc}}$ 为单位再取幂这样 $\gamma1$ 时它就精确退化成标准软阈值 $c\mathrm{sign}(c)\max(|c|-T_{\text{loc}},0)$两者量纲完全一致不会出现幅值小于 1 的系数被取幂后反而变大这种荒谬情况。$\gamma1$ 时超阈值系数被相对抬升边缘更锐$\gamma1$ 时压得更狠噪声残留更少但图像更糊。这个参数交给搜索器去定比人工试快得多。顺带说一句边界处理。做局部方差估计必须用conv2滑窗而conv2的same模式在图像四边会用零填充等于人为制造了一圈暗边局部方差在边缘处会虚高。解决办法是先对系数矩阵做对称延拓padarray(..., symmetric)卷积完再裁回来。这个细节不注意去噪结果的四周会有一圈明显异常尤其在大图上特别扎眼。3. 关键参数的计算过程与取值范围3.1 噪声标准差的三条估计路径$\hat{\sigma}$ 估不准后面全是白搭。最常用的是 MAD 估计取最细尺度对角细节子带 $cD_1$ 的中位数绝对值除以 0.6745。这个 0.6745 是标准正态分布的中位数绝对偏差与标准差的比值来源是 $\Phi^{-1}(0.75)\approx 0.6745$。为什么要用中位数而不是均值因为中位数对离群值稳健即使子带里混进了少量强边缘系数估计值也不会被带跑偏。function sigma estimate_sigma(y, wname) [~, ~, ~, cD1] dwt2(y, wname); sigma median(abs(cD1(:))) / 0.6745; end第二条路径是分块估计把图像切成 32×32 的块对每块单独做 MAD 估计然后取所有块估计值的中位数或 10% 分位数。分位数比中位数更保守因为它偏向取噪声小的块能避免把纹理块的高方差误算成噪声。第三条是局部方差最小法找方差最小的那一批块认为它们近似平坦用这些块的方差做噪声估计。后两条在噪声非均匀的场景下更可靠代价是要多做一轮遍历。注意如果图像本身经过了 JPEG 压缩或者做过锐化最细尺度子带里会残留压缩块效应和锐化振铃MAD 估计会明显偏高。这种情况下建议先在稍粗的尺度上估或者直接对 $\lambda$ 的搜索下界放宽。3.2 搜索空间的设定逻辑搜索空间开得太大收敛慢且容易搜到退化解开得太小搜索器就形同虚设。我一般按下面的范围设定这几个上下界不是拍脑袋的是反复试出来的参数含义下界上界定界理由$\lambda$阈值全局缩放0.43.0低于 0.4 等于没去噪高于 3.0 会削掉大量结构$\gamma$收缩曲线幂次0.71.3越界后要么振铃明显要么过度模糊$\kappa$局部自适应力度0.00.8超过 0.8 后局部估计噪声被放大level小波分解层数24见下文说明分解层数这一维有个特殊处理它必须是整数。如果直接让 QPSO 在 $[2,4]$ 的连续区间上搜再四舍五入会出现大片的平台区——目标函数在 2.3 到 2.7 之间完全相同搜索器拿不到有效梯度信息白白浪费评估次数。我的做法是把层数提到外层做成网格枚举内层只搜三个连续参数最后比较三个层数各自的最优结果。虽然评估次数变成三倍但总次数反而比在四维连续空间里瞎搜少而且结果可解释性强得多。一般图像边长在 512 左右时3 层是最常见的赢家边长 256 以内时 2 层往往更好因为再往下分第 3 层的子带只剩几十个系数统计意义不足。3.3 适应度函数怎么设计才不过拟合最直觉的做法是拿 PSNR 当适应度。但这里有个陷阱如果整张图都参与搜索搜出来的参数会把这幅图特有的噪声实现记住换一张图立刻失效这叫对噪声过拟合。我的处理是把训练图切成互不重叠的两半一半用于搜索参数另一半用于验证适应度也改成两部分的加权$$f -\left(0.6\cdot\text{PSNR}{\text{train}} 0.4\cdot\text{PSNR}{\text{val}}\right)$$取负号是因为 QPSO 是最小化框架。更严格的做法是留出第三份做测试但图像规模有限时两份已经够用。另外我提醒一点如果你手上的数据根本没有干净参考图那 PSNR 这条路径直接断了只能退而求其次用盲指标——比如估计残差图的噪声水平是否接近目标值或者用局部方差的一致性作为代理。这时候搜索出来的参数可信度会明显下降不要把结果当成金标准。4. MATLAB 实现逐模块落地4.1 环境准备与目录组织工具箱只需要两个Image Processing Toolbox读图、padarray、imresize和小波工具箱dwt2、idwt2、wavedec2。MATLAB 版本建议 2019b 及以上主要是为了arguments块的兼容性实际上老版本也能跑把参数校验改成手动判断即可。目录我习惯这么分qab_denoise/ ├─ main_demo.m % 主脚本 ├─ qab_denoise.m % 去噪主函数 ├─ qab_shrink.m % 收缩核 ├─ qab_fitness.m % 适应度封装 ├─ qab_qpso.m % 搜索器 ├─ estimate_sigma.m % 噪声估计 ├─ utils/my_psnr.m % 指标 └─ data/ % 测试图这种拆法的好处是每个文件职责单一改收缩核的时候不用翻主函数调试搜索器的时候可以套一个假的目标函数比如二维球形函数单独验证收敛性。我强烈建议先用一个已知最优点的测试函数把qab_qpso.m跑通再接到去噪上不然一旦结果不好你根本分不清是搜索器有问题还是收缩核有问题。4.2 去噪主函数核心思路是逐层dwt2分解记录每层尺寸用于重构时对齐然后逆向逐层idwt2。function [xHat, info] qab_denoise(y, par) % y : 含噪图像double 类型范围 [0,1] % par : 参数结构体字段 lambda / gamma / kappa / level / wname wname par.wname; level par.level; sigmaN estimate_sigma(y, wname); A y; details cell(level, 1); szPre cell(level, 1); for l 1:level szPre{l} size(A); % 分解前尺寸重构时用来裁剪 [A, cH, cV, cD] dwt2(A, wname); details{l} {cH, cV, cD}; end for l 1:level for k 1:3 c details{l}{k}; Nj numel(c); Tj par.lambda * sigmaN * sqrt(2 * log(Nj)); details{l}{k} qab_shrink(c, Tj, par.gamma, par.kappa, sigmaN); end end for l level:-1:1 d details{l}; A idwt2(A, d{1}, d{2}, d{3}, wname); target szPre{l}; A A(1:target(1), 1:target(2)); % 尺寸对齐避免奇数边长误差累积 end xHat min(max(A, 0), 1); info.sigma sigmaN; info.level level; end这里有个坑值得单说dwt2对奇数尺寸的输入输出会做上取整idwt2还原时可能多出一行一列。如果不在每层还原后立即按原始尺寸裁剪误差会一层层累积三层之后可能差出好几个像素后面算 PSNR 直接对不齐报错。我第一次写的时候就栽在这上面调试了半天以为是阈值算法写错了。小波基我默认用db4。sym4和db4在去噪任务上表现接近db4支撑长度稍短、计算快一点sym4对称性更好、相位失真小一点。如果你的图像里有大量接近对称的结构建筑、机械零件可以试试sym4。别用haar虽然它快但分块感太强重构后能看到明显的方格。也别盲目上更高阶的db8支撑变长意味着边界效应扩散范围更大对中等尺寸以下的图像反而不利。4.3 收缩核的完整实现function cOut qab_shrink(c, T, gamma, kappa, sigmaN) % 局部自适应 幂次软阈值 win 3; ker ones(win) / win^2; pad floor(win / 2); cp padarray(c, [pad pad], symmetric); mu conv2(cp, ker, valid); mu2 conv2(cp, ker, valid); % 占位下行替换 mu2 conv2(cp.^2, ker, valid); varLoc max(mu2 - mu.^2, 0); sigmaLoc sqrt(varLoc); Tloc T .* (1 kappa .* (sigmaLoc ./ (sigmaN eps) - 1)); Tloc max(Tloc, 0.2 * T); % 下限保护防止阈值被压到接近 0 Tloc min(Tloc, 3.0 * T); % 上限保护防止局部纹理把阈值抬飞 mag abs(c); norm1 mag ./ (Tloc eps); cOut sign(c) .* Tloc .* max(norm1 - 1, 0).^gamma; end两个保护上下限是必须的。下限防止某些局部方差极低的区域阈值趋零那些区域的噪声会被完整保留出现随机亮斑上限防止强纹理边缘把阈值抬得过高导致边缘两侧被过度收缩形成光晕。这两个界限看着粗暴但实测能显著提升结果的视觉稳定性。另外eps的位置要注意放在分母上是为了防止除零但如果Tloc本身很小eps不足以影响结果如果Tloc正常eps完全不起作用。不用担心它引入偏差。4.4 QPSO 搜索器实现function [gbest, gbestFit, trace] qab_qpso(fitFun, lb, ub, opt) nPop opt.nPop; maxIt opt.maxIt; dim numel(lb); X repmat(lb, nPop, 1) rand(nPop, dim) .* repmat(ub - lb, nPop, 1); pbest X; pbestFit inf(nPop, 1); gbest X(1, :); gbestFit inf; trace zeros(maxIt, 1); alpha opt.alpha0; for it 1:maxIt for i 1:nPop f fitFun(X(i, :)); if f pbestFit(i) pbestFit(i) f; pbest(i, :) X(i, :); end if f gbestFit gbestFit f; gbest X(i, :); end end mBest mean(pbest, 1); for i 1:nPop phi rand(1, dim); p phi .* pbest(i, :) (1 - phi) .* gbest; u rand(1, dim); L alpha * abs(mBest - X(i, :)); for d 1:dim if rand 0.5 X(i, d) p(d) L(d) * log(1 / u(d)); else X(i, d) p(d) - L(d) * log(1 / u(d)); end end X(i, :) min(max(X(i, :), lb), ub); end alpha alpha * opt.alphaDecay; trace(it) gbestFit; end end参数取值上种群规模 20 到 30 就够再大边际收益很低迭代 60 到 100 代具体看收敛曲线的拐点。alpha0取 1.0alphaDecay取 0.97 到 0.98意味着 60 代后衰减到大约 0.2 到 0.3。我试过固定alpha不衰减结果最后十几代还在大范围跳动最优值一直在变曲线没法收敛。fitFun用匿名函数封装把含噪图和干净参考图绑进去fitFun (p) -qab_fitness(yTrain, p, par0, xTrainClean);qab_fitness里把连续参数映射过去注意层数是外层固定的所以这里只处理三个连续量function val qab_fitness(y, p, par0, xClean) par par0; par.lambda p(1); par.gamma p(2); par.kappa p(3); xHat qab_denoise(y, par); val my_psnr(xHat, xClean); end4.5 主脚本串起来clc; clear; close all; img im2double(imread(data/test01.png)); if size(img, 3) 3 img rgb2gray(img); end rng(2024); sigmaTest 20 / 255; y img sigmaTest * randn(size(img)); y min(max(y, 0), 1); par0.wname db4; lb [0.4, 0.7, 0.0]; ub [3.0, 1.3, 0.8]; opt.nPop 24; opt.maxIt 60; opt.alpha0 1.0; opt.alphaDecay 0.975; levelList [2, 3, 4]; bestAll struct(psnr, -inf); for lv levelList par0.level lv; fitFun (p) -qab_fitness(y, p, par0, img); [gb, gf, tr] qab_qpso(fitFun, lb, ub, opt); parBest par0; parBest.lambda gb(1); parBest.gamma gb(2); parBest.kappa gb(3); xBest qab_denoise(y, parBest); pv my_psnr(xBest, img); fprintf(level%d PSNR%.3f lambda%.3f gamma%.3f kappa%.3f\n, ... lv, pv, gb(1), gb(2), gb(3)); if pv bestAll.psnr bestAll.psnr pv; bestAll.par parBest; bestAll.image xBest; bestAll.trace tr; end end figure; subplot(1,3,1); imshow(img); title(原图); subplot(1,3,2); imshow(y); title(含噪); subplot(1,3,3); imshow(bestAll.image); title(QAB 去噪); figure; plot(bestAll.trace, LineWidth, 1.5); xlabel(迭代代数); ylabel(最优适应度PSNR); grid on;跑一遍大概需要一到三分钟取决于图像尺寸和层数枚举的个数。如果嫌慢可以把种群降到 16、迭代降到 40代价是最终 PSNR 可能掉 0.1 到 0.2 dB。5. 实验结果与参数敏感性观察5.1 测试构造与评价口径评估这块我不喜欢只报一个 PSNR太片面。我一般同时看三个量PSNR 反映整体误差能量SSIM 反映结构相似性还有一个平坦区残差标准差——在手工标注的若干块纯背景区域上计算去噪结果与原图之差的标准差这个值能直接暴露噪声有没有削干净。很多时候 PSNR 涨了但平坦区还是一片噪点说明算法把力气都用去保细节了这种结果放到实际产品里是过不了验收的。测试集构造上我一般选十张左右覆盖不同内容类型的图一张大面积平滑的天空、墙面一张强纹理的草地、织物一张多边缘的建筑、文字再加若干张自然场景图作为综合参考。噪声类型默认加性高斯白噪声$\sigma$ 取 10、20、30、50 四档对应 $\sigma/255$ 的归一化值。如果要测真实噪声建议找带暗光或高 ISO 场景的图但那种图很难拿到干净参考只能靠目视加盲指标。5.2 不同噪声水平下的表现趋势下面这组是量级参考具体数值随图像内容、随机种子、小波基都会浮动别拿它当硬指标对比。噪声 $\sigma$最优 $\lambda$ 落点最优 $\gamma$ 落点最优 $\kappa$ 落点$\lambda$ 固定为 1 时的差距100.7 ~ 0.90.85 ~ 0.950.3 左右约低 0.3 dB201.0 ~ 1.30.9 ~ 1.00.4 ~ 0.5约低 0.5 dB301.4 ~ 1.81.0 ~ 1.10.5 左右约低 0.7 dB501.9 ~ 2.41.1 ~ 1.20.6 左右约低 1.0 dB这张表透露了一个很明确的规律噪声越强最优阈值越高收缩也应该越狠同时局部自适应的力度要更大。这完全符合直觉——噪声大时子带里噪声和细节的重叠更严重需要更激进的阈值把噪声推下去而噪声强的时候局部方差估计相对更可靠信噪比低了但方差主要由噪声贡献所以可以放心加大 $\kappa$。反过来在低噪声下$\gamma$ 和 $\kappa$ 都被压到偏小的值说明这时候主要靠精细的幅度保留取胜过度收缩反而伤害细节。还有一个值得注意的现象$\lambda$ 和 $\gamma$ 之间存在明显的负相关。$\lambda$ 增大时最优 $\gamma$ 倾向于也增大两者在效果上具有补偿关系。这意味着搜索空间里存在一条山谷沿谷底走目标函数变化很小这也是为什么网格搜索在这类问题上效率极低——大量评估花在了重复的谷底点上。5.3 消融实验量子自适应到底贡献了什么拆开看各部分的贡献思路会更清楚。我设计了几组对照第一组是纯标准软阈值$\lambda$ 用经典公式算、$\gamma1$、$\kappa0$第二组开尺度自适应$\lambda$ 靠人工试三档取最优第三组加局部自适应$\kappa$ 人工试三档第四组全开参数由 QPSO 搜索。配置相对基线提升主要改善区域标准软阈值基线— 尺度自适应 $\lambda$0.4 ~ 0.6 dB整体 局部自适应 $\kappa$再 0.3 ~ 0.5 dB纹理区、边缘 幂次收缩 $\gamma$再 0.1 ~ 0.3 dB边缘锐度目视更明显 QPSO 自动搜参再 0.2 ~ 0.4 dB各参数组合的协同从这个拆解能看出来收益的大头来自尺度自适应 局部自适应这两个结构性改动搜索器只是把剩下那点参数协同的空间榨出来。这点认知很重要——如果你时间有限先把局部自适应做好收益比死磕优化器大得多。我见过一些方案把 90% 的篇幅放在优化器上收缩核却还是标准软阈值那属于用力的方向偏了。视觉上最直观的差异在纹理区。标准软阈值处理过的草地、织物纹理会被糊成一片均匀的灰细节全部丢失加上 $\kappa$ 之后纹理区按局部方差自动降低了阈值毛刺感被保留下来同时平坦区的噪点被更强的阈值压下去同一幅图里出现了该软的地方软、该硬的地方硬的效果。这种区域性差异是任何全局阈值都做不到的。6. 常见问题排查实录6.1 PSNR 反降的几种典型成因第一种噪声估计偏高。如果图像有压缩痕迹或者做过锐化MAD 估计会把非噪声的高频成分算进去导致 $T_j$ 整体偏大把所有细节都削平。判断方法是把估计出来的 $\hat{\sigma}$ 和真实加噪强度对比偏差超过 20% 就要警惕。解决办法是改用分块分位数估计或者对 $\hat{\sigma}$ 乘一个 0.8 左右的修正因子。第二种层数选太多。小波分解每深一层低频近似就被不断下采样到了第 4、5 层一个系数对应原图 16×16 甚至 32×32 的区域这个尺度上已经谈不上细节和噪声的区分了强行收缩只会引入块状伪影。如果发现结果里出现了大块的方格基本就是这个原因。判断起来很简单把层数逐次降低重跑一次看伪影是否消失。第三种局部自适应参数过大。$\kappa$ 超过 0.9 之后阈值会紧密跟随局部方差而局部方差本身是从含噪系数算出来的本身就带噪噪声被放大回阈值里形成一种噪声调制阈值、阈值放大噪声的正反馈。表现是图像上出现细密的随机纹理看着更脏了。把 $\kappa$ 压到 0.6 以下基本能解决。第四种尺寸对齐没做。前面提过idwt2的尺寸和原图差一个像素最后算指标时对不齐或者被imresize强行缩放PSNR 直接掉好几个 dB。这种情况的特点是断崖式下跌一眼就能看出来检查一下重构循环里的裁剪语句就行。6.2 运行慢、内存涨的排查路径慢的原因九成在搜索循环里。每次评估都要跑一遍完整去噪包括两层到四层小波分解和重构这是实打实的开销。优化手段有几个按性价比排序第一降低搜索规模种群从 30 降到 20迭代从 80 降到 50评估次数直接砍一半多PSNR 损失通常在 0.1 dB 以内。第二在缩略图上搜参把图降采样到 256×256 上搜出参数再套回原尺寸图因为最优参数对尺度不敏感阈值公式里已经含了 $N_j$ 的补偿这招能省下 80% 以上的时间。第三并行化适应度评估parfor直接套在最内层循环上四核能快三倍左右注意循环内的变量不能有交叉依赖。内存涨的典型场景是参数没有及时清理。details这种 cell 数组在循环里被反复覆盖MATLAB 的写时复制机制在数据被别处引用时会保留旧副本。如果发现跑几十轮之后内存居高不下在每轮循环末尾显式clear一下大变量或者把整个搜索过程包在一个函数里靠作用域自动回收。6.3 问题速查表现象最可能的成因快速验证处理方式结果四周有一圈异常卷积边界零填充放大看四边padarray对称延拓边缘明显发虚$\gamma$ 偏大把 $\gamma$ 手动设 0.8 重跑收紧 $\gamma$ 上界平坦区仍有噪点$\kappa$ 偏大或 $\lambda$ 偏小两者分别单独调整降 $\kappa$、升 $\lambda$出现块状伪影层数过多层数降到 2 重跑限制层数或换sym4搜索曲线一直抖alpha衰减太慢看最后十代是否有下降加快衰减或提前终止换图后效果骤降对噪声实现过拟合用另一张图的验证集测拆分训练/验证集运行报维度不匹配idwt2尺寸误差打印每层size逐层裁剪对齐实操心得调试这套东西的顺序很重要我建议按先固定参数跑通流程、再验证搜索器、最后联合调优三步走。一上来就让搜索器在没有验证过的主函数上跑出问题你根本定位不到是搜索器还是收缩核纯属浪费时间。7. 工程落地前的几点提醒7.1 离线搜参、在线固定别在线上跑优化器这一点我特别想强调因为它直接决定这套方案能不能上生产。QPSO 搜一次参数要几十秒到几分钟而实际部署场景里一张图进来你必须在一秒内出结果甚至要去到几十毫秒。所以正确的用法是在离线阶段用一批有代表性的图对每一档噪声水平、每一类内容平滑型、纹理型、边缘型各搜出一组最优参数做成一张查找表存下来。在线阶段只做一件事——估计输入图的噪声水平然后从表里查最近的参数直接跑一次去噪。这条路走通的代价是泛化性的妥协查表参数一定不如现搜的最优实测差距在 0.2 到 0.5 dB 之间。但从工程角度看这个代价换来的速度提升是几百倍完全值得。而且只要你离线搜参时用的图覆盖足够广在线效果不会崩。还有个小技巧噪声水平不需要估得太准按 10、20、30、50 分成几档实际值落到哪一档就用哪一档的参数中间地带用线性插值。因为参数对噪声水平的敏感度是连续的档位之间差 5 个单位带来的 PSNR 损失通常不到 0.1 dB。7.2 这套框架还能往哪儿长第一个方向是换变换。小波在处理各向异性边缘斜线、曲线时效率不高换成剪切波或者非下采样轮廓波边缘表示能力会明显更强尤其对细长的曲线结构。难点是重构核得自己实现需要仔细验证完美重构条件我建议先在一个可逆性测试上跑通再接到去噪里。第二个方向是换成有监督的参数预测。现在参数是搜出来的能不能训一个小网络输入含噪图直接输出四个参数这样在线阶段连估计噪声那一步都省了。难度在于训练数据的构造——你需要为每张训练图找到真正的最优参数作为标签这个标注成本很高但可以用现在的 QAB 流程离线生成正好形成一个闭环QAB 离线造标签网络在线快速预测。第三个方向是把收缩核换成非局部机制。现在用的是 3×3 邻域本质上还是局部的。如果换成块匹配、相似块分组在组内做协同收缩对重复纹理砖墙、织物、栅栏的去噪效果会有质的提升。代价是计算量上去一到两个数量级得配合积分图或者快速傅里叶变换做加速。这条路的性价比很高尤其是如果你的应用场景里重复纹理多。第四个方向是处理彩色图像。现在整套流程是按单通道灰度写的直接套到 RGB 三个通道上分别处理会有色度失真问题——三个通道的噪声实现是独立的分别去噪后颜色会偏移。正确做法是转换到 YCbCr 或者小波域里做通道间联合收缩对色度通道用更低的阈值人眼对色度噪声不敏感但对色偏很敏感。我在这套东西上折腾了挺长时间最大的一个体会是优化器再花哨也补不上收缩核设计的缺陷。早期我把大量精力花在调 QPSO 的种群规模、衰减系数、随机种子结果 PSNR 在一个范围内来回晃后来把 3×3 局部方差窗口改成 5×5把阈值下限保护加进去一次就往上跳了 0.3 dB。参数搜索能帮你挖到当前设计的天花板但它抬不高天花板本身。所以真要往深里做先把收缩核想清楚再考虑优化器要不要换更复杂的。另外再分享一个调试时管用的小习惯把每次去噪的残差图含噪图减去去噪图存下来看一眼。理想情况下残差应该接近纯白噪声看不见任何结构。如果残差里能看出明显的边缘轮廓说明阈值削过头了把结构当噪声削掉了如果残差里能看到大块的斑驳区域说明局部自适应在那些地方过于激进了。这一眼比看 PSNR 数字直观得多很多时候 PSNR 涨了但残差里有结构那就是在过拟合。