搜索的圖像去噪MATLAB實現(xiàn))
圖像去噪做到后面瓶頸基本不在會不會調(diào)函數(shù)而在閾值怎么定、定在哪一層、定在哪個位置。我這套基于量子自適應變換 QAB 算法的圖像去噪系統(tǒng)本質(zhì)上就是把變換域去噪里那幾個最靠人肉試參的環(huán)節(jié)——分解層數(shù)、閾值縮放系數(shù)、收縮曲線的形狀、局部自適應的力度——打包成一個可搜索的參數(shù)空間再交給一個量子行為啟發(fā)的優(yōu)化器去自己找。整套東西用MATLAB實現(xiàn)純 CPU 跑不依賴任何特殊硬件也不需要顯卡一臺普通筆記本就能復現(xiàn)。它解決的是傳統(tǒng)小波軟閾值去噪里最讓人難受的三件事固定閾值對不同噪聲強度不友好、對紋理區(qū)域過度平滑、對平坦區(qū)域殘留噪聲。適合誰看手上有一批含噪圖像、想搭一條能自洽評估的去噪流水線的人正在寫相關方向論文、需要一個說得清為什么這么設計的完整方案的人以及已經(jīng)會wavedec2和idwt2、但發(fā)現(xiàn)調(diào)半天 PSNR 上不去、想找系統(tǒng)化調(diào)參思路的人。下面我會把選型邏輯、公式推導、可復制的代碼、踩過的坑全部攤開講。1. 標題里真正要解決的技術命題1.1 傳統(tǒng)變換域去噪的三處硬傷先別急著談量子。任何一套變換域去噪拆到最底層就是三步換一個域、對系數(shù)做收縮、再換回來。小波之所以長期占主流是因為它把圖像拆成低頻近似 高頻細節(jié)的多尺度金字塔噪聲在細節(jié)子帶里呈現(xiàn)為能量分散的小幅值系數(shù)而邊緣和紋理是大體量系數(shù)兩者在幅值上天然可分離——這就是收縮能工作的物理前提。但真拿它做工程問題立刻浮現(xiàn)。第一處硬傷是閾值是全局的。經(jīng)典的通用閾值 $T\sigma\sqrt{2\ln N}$ 只跟圖像總像素數(shù) $N$ 和噪聲標準差 $\sigma$ 有關它對整幅圖一視同仁。可現(xiàn)實中的圖像從來不是均勻的天空、墻面這種區(qū)域噪聲明顯草叢、毛發(fā)這種區(qū)域細節(jié)密集卻容易被誤殺。一刀切的閾值在平坦區(qū)會留下斑點在紋理區(qū)會把細節(jié)削平。第二處硬傷是收縮曲線太軟。標準軟閾值 $c\mathrm{sign}(c)\max(|c|-T,0)$ 對超過閾值的系數(shù)也統(tǒng)一減掉 $T$導致所有邊緣被系統(tǒng)性壓低圖像看上去發(fā)悶、發(fā)糊。硬閾值雖然保幅值但會在閾值附近產(chǎn)生跳變重構后出現(xiàn)振鈴和偽影。第三處硬傷是層間不一致。第 1 層細節(jié)子帶噪聲占比高第 3 層幾乎全是結構信息用同一個縮放系數(shù)的結果往往是低層沒削干凈、高層被削過頭。這三處硬傷指向同一個結論閾值不是一個標量而應該是一個隨尺度、隨空間位置變化的函數(shù)而描述這個函數(shù)的參數(shù)人的手調(diào)不過來。這就是我要引入自適應搜索器的動機。1.2 QAB 的骨架與選型理由QAB 是我給這套方案起的內(nèi)部代號全稱對應Quantum-behaved Adaptive Basis直白翻譯就是量子行為搜索 自適應變換基收縮。它的骨架拆成三塊第一塊是變換模塊用二維離散小波把圖像分解成多尺度子帶選小波是因為它在 MATLAB 里工具鏈最完整、可逆性最好、邊界處理成熟換成曲波或剪切波理論上前景更好但在沒有現(xiàn)成工具箱的情況下自己寫重構核的代價太高工程上不劃算。第二塊是自適應收縮核我不用固定閾值而是讓閾值隨局部鄰域方差聯(lián)動并且給收縮曲線加了一個冪次因子來控制保邊強度。第三塊是參數(shù)搜索器用帶量子行為的粒子群算法在連續(xù)參數(shù)空間里搜索最優(yōu)組合。為什么搜參器選 QPSO 而不是量子遺傳算法或者網(wǎng)格搜索網(wǎng)格搜索在三維連續(xù)參數(shù)上就是災難每維取 10 個點三維就是 1000 次完整去噪單次算力約幾百毫秒加起來半小時起步還只能拿到粗糙的格點。遺傳算法需要把連續(xù)量離散成比特串編碼精度受位長限制而且雜交、變異算子對這類平滑性較好的目標函數(shù)并不占優(yōu)。QPSO 用的是勢阱模型的采樣方式位置更新里帶一個 $\ln(1/u)$ 的對數(shù)項天然具備跳出局部極值的能力同時參數(shù)是實數(shù)編碼不需要離散化再反解精度上限高得多。我實測下來在 3 到 5 維、每維取值范圍可控的情況下QPSO 大約 40 到 60 代就能收斂到一個相當穩(wěn)的解比 GA 穩(wěn)定比網(wǎng)格快一到兩個數(shù)量級。這里必須提前說清楚一件事避免誤解QAB 里的量子是算法思想層面的指的是用概率幅、疊加態(tài)采樣這類機制來設計更新規(guī)則不是真的跑到量子計算機上執(zhí)行。整套代碼在經(jīng)典 CPU 上用雙精度浮點跑沒有任何硬件門檻。很多初學者一看到量子兩個字就以為要裝什么特殊環(huán)境實際上它跟量子力學實驗沒有任何關系只是把一個基于概率分布的啟發(fā)式搜索換了套數(shù)學外衣。2. 量子自適應變換的核心機制拆解2.1 概率幅表示與勢阱采樣到底在做什么傳統(tǒng)粒子群里粒子位置是確定的速度也是確定的更新公式 $v \leftarrow \omega v c_1 r_1 (p-x) c_2 r_2 (g-x)$ 本質(zhì)上是帶上限的梯度下降。它的問題在于一旦所有粒子都擠在局部最優(yōu)附近速度項趨零整個種群就僵住了。量子行為的粒子群不用速度和位置這套描述改用波函數(shù)粒子不再有確定軌跡而是以某個概率密度分布在搜索空間里這個密度由勢阱中心決定。具體到代碼上更新式寫成$$x_{i,d} p_d \pm \alpha,|m_d - x_{i,d}|,\ln!\left(\frac{1}{u}\right)$$其中 $p_d$ 是粒子自身最優(yōu)與全局最優(yōu)的隨機加權用一個隨機數(shù) $\phi$ 混合$m_d$ 是全體個體最優(yōu)的均值也就是種群的質(zhì)心$\alpha$ 是收縮擴張系數(shù)$u$ 是 $(0,1)$ 上的均勻隨機數(shù)。那個 $\ln(1/u)$ 是這套算法最精髓的地方$u$ 接近 1 時它趨近 0粒子幾乎不動做精細開發(fā)$u$ 接近 0 時它是個大數(shù)粒子一步跳很遠做全局探索。而且它是重尾分布跳遠的概率雖然小但不是零這保證了即使種群暫時陷入某個洼地也有機會被甩出去。$\alpha$ 的作用是控制搜索半徑的整體縮放。我一般讓它從 1.0 線性衰減到 0.2前期大步掃后期小步精修。這個衰減節(jié)奏不能太快衰減過快會導致早熟也不能太慢慢了最后幾代還在大范圍亂跳收斂曲線抖得看不清趨勢。2.2 自適應基閾值必須隨子帶和位置一起變量子那部分只是找參的工具真正決定去噪質(zhì)量的是收縮核怎么設計。我把閾值寫成兩層結構。第一層是尺度相關第 $j$ 層細節(jié)子帶的閾值基準取$$T_j \lambda,\hat{\sigma}\sqrt{2\ln N_j}$$其中 $N_j$ 是該子帶的系數(shù)個數(shù)$\hat{\sigma}$ 是噪聲標準差估計$\lambda$ 是一個全局縮放系數(shù)交給搜索器去調(diào)。為什么用 $N_j$ 而不是總像素數(shù) $N$因為子帶系數(shù)個數(shù)隨層數(shù)遞減低層子帶系數(shù)多、噪聲能量分散算出來的閾值相對高一些高層子帶系數(shù)少閾值自動降下來避免把結構信息當噪聲削掉。這個改動看著小實測對紋理豐富的圖像能穩(wěn)定貢獻接近 0.2 到 0.4 dB 的提升。第二層是空間相關在子帶內(nèi)部再按 3×3 鄰域的局部標準差對閾值做微調(diào)$$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 就完全跟隨局部方差。這個式子的直覺是局部方差明顯高于全局噪聲水平的區(qū)域大概率是有結構的邊緣紋理應當減弱收縮力度保護細節(jié)局部方差接近噪聲水平的區(qū)域就是純噪聲應當加強收縮。注意閾值調(diào)的是基準的高低不是直接決定保留多少所以 $\kappa$ 過大也會出問題——局部方差估計本身帶噪$\kappa$ 太大會把噪聲估計的抖動放大成閾值的抖動反而制造塊狀偽影。實操中我把 $\kappa$ 限制在 $[0, 0.8]$。2.3 冪次軟閾值解決保邊和去噪的矛盾標準軟閾值的問題是它把所有超閾值系數(shù)統(tǒng)一平移 $T$幅值越大的系數(shù)被壓得越多相對比例上邊緣就軟了。我用的改進形式是把收縮寫成$$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)$兩者量綱完全一致不會出現(xiàn)幅值小于 1 的系數(shù)被取冪后反而變大這種荒謬情況。$\gamma1$ 時超閾值系數(shù)被相對抬升邊緣更銳$\gamma1$ 時壓得更狠噪聲殘留更少但圖像更糊。這個參數(shù)交給搜索器去定比人工試快得多。順帶說一句邊界處理。做局部方差估計必須用conv2滑窗而conv2的same模式在圖像四邊會用零填充等于人為制造了一圈暗邊局部方差在邊緣處會虛高。解決辦法是先對系數(shù)矩陣做對稱延拓padarray(..., symmetric)卷積完再裁回來。這個細節(jié)不注意去噪結果的四周會有一圈明顯異常尤其在大圖上特別扎眼。3. 關鍵參數(shù)的計算過程與取值范圍3.1 噪聲標準差的三條估計路徑$\hat{\sigma}$ 估不準后面全是白搭。最常用的是 MAD 估計取最細尺度對角細節(jié)子帶 $cD_1$ 的中位數(shù)絕對值除以 0.6745。這個 0.6745 是標準正態(tài)分布的中位數(shù)絕對偏差與標準差的比值來源是 $\Phi^{-1}(0.75)\approx 0.6745$。為什么要用中位數(shù)而不是均值因為中位數(shù)對離群值穩(wěn)健即使子帶里混進了少量強邊緣系數(shù)估計值也不會被帶跑偏。function sigma estimate_sigma(y, wname) [~, ~, ~, cD1] dwt2(y, wname); sigma median(abs(cD1(:))) / 0.6745; end第二條路徑是分塊估計把圖像切成 32×32 的塊對每塊單獨做 MAD 估計然后取所有塊估計值的中位數(shù)或 10% 分位數(shù)。分位數(shù)比中位數(shù)更保守因為它偏向取噪聲小的塊能避免把紋理塊的高方差誤算成噪聲。第三條是局部方差最小法找方差最小的那一批塊認為它們近似平坦用這些塊的方差做噪聲估計。后兩條在噪聲非均勻的場景下更可靠代價是要多做一輪遍歷。注意如果圖像本身經(jīng)過了 JPEG 壓縮或者做過銳化最細尺度子帶里會殘留壓縮塊效應和銳化振鈴MAD 估計會明顯偏高。這種情況下建議先在稍粗的尺度上估或者直接對 $\lambda$ 的搜索下界放寬。3.2 搜索空間的設定邏輯搜索空間開得太大收斂慢且容易搜到退化解開得太小搜索器就形同虛設。我一般按下面的范圍設定這幾個上下界不是拍腦袋的是反復試出來的參數(shù)含義下界上界定界理由$\lambda$閾值全局縮放0.43.0低于 0.4 等于沒去噪高于 3.0 會削掉大量結構$\gamma$收縮曲線冪次0.71.3越界后要么振鈴明顯要么過度模糊$\kappa$局部自適應力度0.00.8超過 0.8 后局部估計噪聲被放大level小波分解層數(shù)24見下文說明分解層數(shù)這一維有個特殊處理它必須是整數(shù)。如果直接讓 QPSO 在 $[2,4]$ 的連續(xù)區(qū)間上搜再四舍五入會出現(xiàn)大片的平臺區(qū)——目標函數(shù)在 2.3 到 2.7 之間完全相同搜索器拿不到有效梯度信息白白浪費評估次數(shù)。我的做法是把層數(shù)提到外層做成網(wǎng)格枚舉內(nèi)層只搜三個連續(xù)參數(shù)最后比較三個層數(shù)各自的最優(yōu)結果。雖然評估次數(shù)變成三倍但總次數(shù)反而比在四維連續(xù)空間里瞎搜少而且結果可解釋性強得多。一般圖像邊長在 512 左右時3 層是最常見的贏家邊長 256 以內(nèi)時 2 層往往更好因為再往下分第 3 層的子帶只剩幾十個系數(shù)統(tǒng)計意義不足。3.3 適應度函數(shù)怎么設計才不過擬合最直覺的做法是拿 PSNR 當適應度。但這里有個陷阱如果整張圖都參與搜索搜出來的參數(shù)會把這幅圖特有的噪聲實現(xiàn)記住換一張圖立刻失效這叫對噪聲過擬合。我的處理是把訓練圖切成互不重疊的兩半一半用于搜索參數(shù)另一半用于驗證適應度也改成兩部分的加權$$f -\left(0.6\cdot\text{PSNR}{\text{train}} 0.4\cdot\text{PSNR}{\text{val}}\right)$$取負號是因為 QPSO 是最小化框架。更嚴格的做法是留出第三份做測試但圖像規(guī)模有限時兩份已經(jīng)夠用。另外我提醒一點如果你手上的數(shù)據(jù)根本沒有干凈參考圖那 PSNR 這條路徑直接斷了只能退而求其次用盲指標——比如估計殘差圖的噪聲水平是否接近目標值或者用局部方差的一致性作為代理。這時候搜索出來的參數(shù)可信度會明顯下降不要把結果當成金標準。4. MATLAB 實現(xiàn)逐模塊落地4.1 環(huán)境準備與目錄組織工具箱只需要兩個Image Processing Toolbox讀圖、padarray、imresize和小波工具箱dwt2、idwt2、wavedec2。MATLAB 版本建議 2019b 及以上主要是為了arguments塊的兼容性實際上老版本也能跑把參數(shù)校驗改成手動判斷即可。目錄我習慣這么分qab_denoise/ ├─ main_demo.m % 主腳本 ├─ qab_denoise.m % 去噪主函數(shù) ├─ qab_shrink.m % 收縮核 ├─ qab_fitness.m % 適應度封裝 ├─ qab_qpso.m % 搜索器 ├─ estimate_sigma.m % 噪聲估計 ├─ utils/my_psnr.m % 指標 └─ data/ % 測試圖這種拆法的好處是每個文件職責單一改收縮核的時候不用翻主函數(shù)調(diào)試搜索器的時候可以套一個假的目標函數(shù)比如二維球形函數(shù)單獨驗證收斂性。我強烈建議先用一個已知最優(yōu)點的測試函數(shù)把qab_qpso.m跑通再接到去噪上不然一旦結果不好你根本分不清是搜索器有問題還是收縮核有問題。4.2 去噪主函數(shù)核心思路是逐層dwt2分解記錄每層尺寸用于重構時對齊然后逆向逐層idwt2。function [xHat, info] qab_denoise(y, par) % y : 含噪圖像double 類型范圍 [0,1] % par : 參數(shù)結構體字段 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)); % 尺寸對齊避免奇數(shù)邊長誤差累積 end xHat min(max(A, 0), 1); info.sigma sigmaN; info.level level; end這里有個坑值得單說dwt2對奇數(shù)尺寸的輸入輸出會做上取整idwt2還原時可能多出一行一列。如果不在每層還原后立即按原始尺寸裁剪誤差會一層層累積三層之后可能差出好幾個像素后面算 PSNR 直接對不齊報錯。我第一次寫的時候就栽在這上面調(diào)試了半天以為是閾值算法寫錯了。小波基我默認用db4。sym4和db4在去噪任務上表現(xiàn)接近db4支撐長度稍短、計算快一點sym4對稱性更好、相位失真小一點。如果你的圖像里有大量接近對稱的結構建筑、機械零件可以試試sym4。別用haar雖然它快但分塊感太強重構后能看到明顯的方格。也別盲目上更高階的db8支撐變長意味著邊界效應擴散范圍更大對中等尺寸以下的圖像反而不利。4.3 收縮核的完整實現(xiàn)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兩個保護上下限是必須的。下限防止某些局部方差極低的區(qū)域閾值趨零那些區(qū)域的噪聲會被完整保留出現(xiàn)隨機亮斑上限防止強紋理邊緣把閾值抬得過高導致邊緣兩側被過度收縮形成光暈。這兩個界限看著粗暴但實測能顯著提升結果的視覺穩(wěn)定性。另外eps的位置要注意放在分母上是為了防止除零但如果Tloc本身很小eps不足以影響結果如果Tloc正常eps完全不起作用。不用擔心它引入偏差。4.4 QPSO 搜索器實現(xiàn)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參數(shù)取值上種群規(guī)模 20 到 30 就夠再大邊際收益很低迭代 60 到 100 代具體看收斂曲線的拐點。alpha0取 1.0alphaDecay取 0.97 到 0.98意味著 60 代后衰減到大約 0.2 到 0.3。我試過固定alpha不衰減結果最后十幾代還在大范圍跳動最優(yōu)值一直在變曲線沒法收斂。fitFun用匿名函數(shù)封裝把含噪圖和干凈參考圖綁進去fitFun (p) -qab_fitness(yTrain, p, par0, xTrainClean);qab_fitness里把連續(xù)參數(shù)映射過去注意層數(shù)是外層固定的所以這里只處理三個連續(xù)量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(迭代代數(shù)); ylabel(最優(yōu)適應度PSNR); grid on;跑一遍大概需要一到三分鐘取決于圖像尺寸和層數(shù)枚舉的個數(shù)。如果嫌慢可以把種群降到 16、迭代降到 40代價是最終 PSNR 可能掉 0.1 到 0.2 dB。5. 實驗結果與參數(shù)敏感性觀察5.1 測試構造與評價口徑評估這塊我不喜歡只報一個 PSNR太片面。我一般同時看三個量PSNR 反映整體誤差能量SSIM 反映結構相似性還有一個平坦區(qū)殘差標準差——在手工標注的若干塊純背景區(qū)域上計算去噪結果與原圖之差的標準差這個值能直接暴露噪聲有沒有削干凈。很多時候 PSNR 漲了但平坦區(qū)還是一片噪點說明算法把力氣都用去保細節(jié)了這種結果放到實際產(chǎn)品里是過不了驗收的。測試集構造上我一般選十張左右覆蓋不同內(nèi)容類型的圖一張大面積平滑的天空、墻面一張強紋理的草地、織物一張多邊緣的建筑、文字再加若干張自然場景圖作為綜合參考。噪聲類型默認加性高斯白噪聲$\sigma$ 取 10、20、30、50 四檔對應 $\sigma/255$ 的歸一化值。如果要測真實噪聲建議找?guī)О倒饣蚋?ISO 場景的圖但那種圖很難拿到干凈參考只能靠目視加盲指標。5.2 不同噪聲水平下的表現(xiàn)趨勢下面這組是量級參考具體數(shù)值隨圖像內(nèi)容、隨機種子、小波基都會浮動別拿它當硬指標對比。噪聲 $\sigma$最優(yōu) $\lambda$ 落點最優(yōu) $\gamma$ 落點最優(yōu) $\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這張表透露了一個很明確的規(guī)律噪聲越強最優(yōu)閾值越高收縮也應該越狠同時局部自適應的力度要更大。這完全符合直覺——噪聲大時子帶里噪聲和細節(jié)的重疊更嚴重需要更激進的閾值把噪聲推下去而噪聲強的時候局部方差估計相對更可靠信噪比低了但方差主要由噪聲貢獻所以可以放心加大 $\kappa$。反過來在低噪聲下$\gamma$ 和 $\kappa$ 都被壓到偏小的值說明這時候主要靠精細的幅度保留取勝過度收縮反而傷害細節(jié)。還有一個值得注意的現(xiàn)象$\lambda$ 和 $\gamma$ 之間存在明顯的負相關。$\lambda$ 增大時最優(yōu) $\gamma$ 傾向于也增大兩者在效果上具有補償關系。這意味著搜索空間里存在一條山谷沿谷底走目標函數(shù)變化很小這也是為什么網(wǎng)格搜索在這類問題上效率極低——大量評估花在了重復的谷底點上。5.3 消融實驗量子自適應到底貢獻了什么拆開看各部分的貢獻思路會更清楚。我設計了幾組對照第一組是純標準軟閾值$\lambda$ 用經(jīng)典公式算、$\gamma1$、$\kappa0$第二組開尺度自適應$\lambda$ 靠人工試三檔取最優(yōu)第三組加局部自適應$\kappa$ 人工試三檔第四組全開參數(shù)由 QPSO 搜索。配置相對基線提升主要改善區(qū)域標準軟閾值基線— 尺度自適應 $\lambda$0.4 ~ 0.6 dB整體 局部自適應 $\kappa$再 0.3 ~ 0.5 dB紋理區(qū)、邊緣 冪次收縮 $\gamma$再 0.1 ~ 0.3 dB邊緣銳度目視更明顯 QPSO 自動搜參再 0.2 ~ 0.4 dB各參數(shù)組合的協(xié)同從這個拆解能看出來收益的大頭來自尺度自適應 局部自適應這兩個結構性改動搜索器只是把剩下那點參數(shù)協(xié)同的空間榨出來。這點認知很重要——如果你時間有限先把局部自適應做好收益比死磕優(yōu)化器大得多。我見過一些方案把 90% 的篇幅放在優(yōu)化器上收縮核卻還是標準軟閾值那屬于用力的方向偏了。視覺上最直觀的差異在紋理區(qū)。標準軟閾值處理過的草地、織物紋理會被糊成一片均勻的灰細節(jié)全部丟失加上 $\kappa$ 之后紋理區(qū)按局部方差自動降低了閾值毛刺感被保留下來同時平坦區(qū)的噪點被更強的閾值壓下去同一幅圖里出現(xiàn)了該軟的地方軟、該硬的地方硬的效果。這種區(qū)域性差異是任何全局閾值都做不到的。6. 常見問題排查實錄6.1 PSNR 反降的幾種典型成因第一種噪聲估計偏高。如果圖像有壓縮痕跡或者做過銳化MAD 估計會把非噪聲的高頻成分算進去導致 $T_j$ 整體偏大把所有細節(jié)都削平。判斷方法是把估計出來的 $\hat{\sigma}$ 和真實加噪強度對比偏差超過 20% 就要警惕。解決辦法是改用分塊分位數(shù)估計或者對 $\hat{\sigma}$ 乘一個 0.8 左右的修正因子。第二種層數(shù)選太多。小波分解每深一層低頻近似就被不斷下采樣到了第 4、5 層一個系數(shù)對應原圖 16×16 甚至 32×32 的區(qū)域這個尺度上已經(jīng)談不上細節(jié)和噪聲的區(qū)分了強行收縮只會引入塊狀偽影。如果發(fā)現(xiàn)結果里出現(xiàn)了大塊的方格基本就是這個原因。判斷起來很簡單把層數(shù)逐次降低重跑一次看偽影是否消失。第三種局部自適應參數(shù)過大。$\kappa$ 超過 0.9 之后閾值會緊密跟隨局部方差而局部方差本身是從含噪系數(shù)算出來的本身就帶噪噪聲被放大回閾值里形成一種噪聲調(diào)制閾值、閾值放大噪聲的正反饋。表現(xiàn)是圖像上出現(xiàn)細密的隨機紋理看著更臟了。把 $\kappa$ 壓到 0.6 以下基本能解決。第四種尺寸對齊沒做。前面提過idwt2的尺寸和原圖差一個像素最后算指標時對不齊或者被imresize強行縮放PSNR 直接掉好幾個 dB。這種情況的特點是斷崖式下跌一眼就能看出來檢查一下重構循環(huán)里的裁剪語句就行。6.2 運行慢、內(nèi)存漲的排查路徑慢的原因九成在搜索循環(huán)里。每次評估都要跑一遍完整去噪包括兩層到四層小波分解和重構這是實打實的開銷。優(yōu)化手段有幾個按性價比排序第一降低搜索規(guī)模種群從 30 降到 20迭代從 80 降到 50評估次數(shù)直接砍一半多PSNR 損失通常在 0.1 dB 以內(nèi)。第二在縮略圖上搜參把圖降采樣到 256×256 上搜出參數(shù)再套回原尺寸圖因為最優(yōu)參數(shù)對尺度不敏感閾值公式里已經(jīng)含了 $N_j$ 的補償這招能省下 80% 以上的時間。第三并行化適應度評估parfor直接套在最內(nèi)層循環(huán)上四核能快三倍左右注意循環(huán)內(nèi)的變量不能有交叉依賴。內(nèi)存漲的典型場景是參數(shù)沒有及時清理。details這種 cell 數(shù)組在循環(huán)里被反復覆蓋MATLAB 的寫時復制機制在數(shù)據(jù)被別處引用時會保留舊副本。如果發(fā)現(xiàn)跑幾十輪之后內(nèi)存居高不下在每輪循環(huán)末尾顯式clear一下大變量或者把整個搜索過程包在一個函數(shù)里靠作用域自動回收。6.3 問題速查表現(xiàn)象最可能的成因快速驗證處理方式結果四周有一圈異常卷積邊界零填充放大看四邊padarray對稱延拓邊緣明顯發(fā)虛$\gamma$ 偏大把 $\gamma$ 手動設 0.8 重跑收緊 $\gamma$ 上界平坦區(qū)仍有噪點$\kappa$ 偏大或 $\lambda$ 偏小兩者分別單獨調(diào)整降 $\kappa$、升 $\lambda$出現(xiàn)塊狀偽影層數(shù)過多層數(shù)降到 2 重跑限制層數(shù)或換sym4搜索曲線一直抖alpha衰減太慢看最后十代是否有下降加快衰減或提前終止換圖后效果驟降對噪聲實現(xiàn)過擬合用另一張圖的驗證集測拆分訓練/驗證集運行報維度不匹配idwt2尺寸誤差打印每層size逐層裁剪對齊實操心得調(diào)試這套東西的順序很重要我建議按先固定參數(shù)跑通流程、再驗證搜索器、最后聯(lián)合調(diào)優(yōu)三步走。一上來就讓搜索器在沒有驗證過的主函數(shù)上跑出問題你根本定位不到是搜索器還是收縮核純屬浪費時間。7. 工程落地前的幾點提醒7.1 離線搜參、在線固定別在線上跑優(yōu)化器這一點我特別想強調(diào)因為它直接決定這套方案能不能上生產(chǎn)。QPSO 搜一次參數(shù)要幾十秒到幾分鐘而實際部署場景里一張圖進來你必須在一秒內(nèi)出結果甚至要去到幾十毫秒。所以正確的用法是在離線階段用一批有代表性的圖對每一檔噪聲水平、每一類內(nèi)容平滑型、紋理型、邊緣型各搜出一組最優(yōu)參數(shù)做成一張查找表存下來。在線階段只做一件事——估計輸入圖的噪聲水平然后從表里查最近的參數(shù)直接跑一次去噪。這條路走通的代價是泛化性的妥協(xié)查表參數(shù)一定不如現(xiàn)搜的最優(yōu)實測差距在 0.2 到 0.5 dB 之間。但從工程角度看這個代價換來的速度提升是幾百倍完全值得。而且只要你離線搜參時用的圖覆蓋足夠廣在線效果不會崩。還有個小技巧噪聲水平不需要估得太準按 10、20、30、50 分成幾檔實際值落到哪一檔就用哪一檔的參數(shù)中間地帶用線性插值。因為參數(shù)對噪聲水平的敏感度是連續(xù)的檔位之間差 5 個單位帶來的 PSNR 損失通常不到 0.1 dB。7.2 這套框架還能往哪兒長第一個方向是換變換。小波在處理各向異性邊緣斜線、曲線時效率不高換成剪切波或者非下采樣輪廓波邊緣表示能力會明顯更強尤其對細長的曲線結構。難點是重構核得自己實現(xiàn)需要仔細驗證完美重構條件我建議先在一個可逆性測試上跑通再接到去噪里。第二個方向是換成有監(jiān)督的參數(shù)預測?，F(xiàn)在參數(shù)是搜出來的能不能訓一個小網(wǎng)絡輸入含噪圖直接輸出四個參數(shù)這樣在線階段連估計噪聲那一步都省了。難度在于訓練數(shù)據(jù)的構造——你需要為每張訓練圖找到真正的最優(yōu)參數(shù)作為標簽這個標注成本很高但可以用現(xiàn)在的 QAB 流程離線生成正好形成一個閉環(huán)QAB 離線造標簽網(wǎng)絡在線快速預測。第三個方向是把收縮核換成非局部機制?，F(xiàn)在用的是 3×3 鄰域本質(zhì)上還是局部的。如果換成塊匹配、相似塊分組在組內(nèi)做協(xié)同收縮對重復紋理磚墻、織物、柵欄的去噪效果會有質(zhì)的提升。代價是計算量上去一到兩個數(shù)量級得配合積分圖或者快速傅里葉變換做加速。這條路的性價比很高尤其是如果你的應用場景里重復紋理多。第四個方向是處理彩色圖像?，F(xiàn)在整套流程是按單通道灰度寫的直接套到 RGB 三個通道上分別處理會有色度失真問題——三個通道的噪聲實現(xiàn)是獨立的分別去噪后顏色會偏移。正確做法是轉換到 YCbCr 或者小波域里做通道間聯(lián)合收縮對色度通道用更低的閾值人眼對色度噪聲不敏感但對色偏很敏感。我在這套東西上折騰了挺長時間最大的一個體會是優(yōu)化器再花哨也補不上收縮核設計的缺陷。早期我把大量精力花在調(diào) QPSO 的種群規(guī)模、衰減系數(shù)、隨機種子結果 PSNR 在一個范圍內(nèi)來回晃后來把 3×3 局部方差窗口改成 5×5把閾值下限保護加進去一次就往上跳了 0.3 dB。參數(shù)搜索能幫你挖到當前設計的天花板但它抬不高天花板本身。所以真要往深里做先把收縮核想清楚再考慮優(yōu)化器要不要換更復雜的。另外再分享一個調(diào)試時管用的小習慣把每次去噪的殘差圖含噪圖減去去噪圖存下來看一眼。理想情況下殘差應該接近純白噪聲看不見任何結構。如果殘差里能看出明顯的邊緣輪廓說明閾值削過頭了把結構當噪聲削掉了如果殘差里能看到大塊的斑駁區(qū)域說明局部自適應在那些地方過于激進了。這一眼比看 PSNR 數(shù)字直觀得多很多時候 PSNR 漲了但殘差里有結構那就是在過擬合。 基于QAB量子自適應變換與QPSO參數(shù)搜索的圖像去噪MATLAB實現(xiàn)
圖像去噪做到后面瓶頸基本不在會不會調(diào)函數(shù)而在閾值怎么定、定在哪一層、定在哪個位置。我這套基于量子自適應變換 QAB 算法的圖像去噪系統(tǒng)本質(zhì)上就是把變換域去噪里那幾個最靠人肉試參的環(huán)節(jié)——分解層數(shù)、閾值縮放系數(shù)、收縮曲線的形狀、局部自適應的力度——打包成一個可搜索的參數(shù)空間再交給一個量子行為啟發(fā)的優(yōu)化器去自己找。整套東西用MATLAB實現(xiàn)純 CPU 跑不依賴任何特殊硬件也不需要顯卡一臺普通筆記本就能復現(xiàn)。它解決的是傳統(tǒng)小波軟閾值去噪里最讓人難受的三件事固定閾值對不同噪聲強度不友好、對紋理區(qū)域過度平滑、對平坦區(qū)域殘留噪聲。適合誰看手上有一批含噪圖像、想搭一條能自洽評估的去噪流水線的人正在寫相關方向論文、需要一個說得清為什么這么設計的完整方案的人以及已經(jīng)會wavedec2和idwt2、但發(fā)現(xiàn)調(diào)半天 PSNR 上不去、想找系統(tǒng)化調(diào)參思路的人。下面我會把選型邏輯、公式推導、可復制的代碼、踩過的坑全部攤開講。1. 標題里真正要解決的技術命題1.1 傳統(tǒng)變換域去噪的三處硬傷先別急著談量子。任何一套變換域去噪拆到最底層就是三步換一個域、對系數(shù)做收縮、再換回來。小波之所以長期占主流是因為它把圖像拆成低頻近似 高頻細節(jié)的多尺度金字塔噪聲在細節(jié)子帶里呈現(xiàn)為能量分散的小幅值系數(shù)而邊緣和紋理是大體量系數(shù)兩者在幅值上天然可分離——這就是收縮能工作的物理前提。但真拿它做工程問題立刻浮現(xiàn)。第一處硬傷是閾值是全局的。經(jīng)典的通用閾值 $T\sigma\sqrt{2\ln N}$ 只跟圖像總像素數(shù) $N$ 和噪聲標準差 $\sigma$ 有關它對整幅圖一視同仁。可現(xiàn)實中的圖像從來不是均勻的天空、墻面這種區(qū)域噪聲明顯草叢、毛發(fā)這種區(qū)域細節(jié)密集卻容易被誤殺。一刀切的閾值在平坦區(qū)會留下斑點在紋理區(qū)會把細節(jié)削平。第二處硬傷是收縮曲線太軟。標準軟閾值 $c\mathrm{sign}(c)\max(|c|-T,0)$ 對超過閾值的系數(shù)也統(tǒng)一減掉 $T$導致所有邊緣被系統(tǒng)性壓低圖像看上去發(fā)悶、發(fā)糊。硬閾值雖然保幅值但會在閾值附近產(chǎn)生跳變重構后出現(xiàn)振鈴和偽影。第三處硬傷是層間不一致。第 1 層細節(jié)子帶噪聲占比高第 3 層幾乎全是結構信息用同一個縮放系數(shù)的結果往往是低層沒削干凈、高層被削過頭。這三處硬傷指向同一個結論閾值不是一個標量而應該是一個隨尺度、隨空間位置變化的函數(shù)而描述這個函數(shù)的參數(shù)人的手調(diào)不過來。這就是我要引入自適應搜索器的動機。1.2 QAB 的骨架與選型理由QAB 是我給這套方案起的內(nèi)部代號全稱對應Quantum-behaved Adaptive Basis直白翻譯就是量子行為搜索 自適應變換基收縮。它的骨架拆成三塊第一塊是變換模塊用二維離散小波把圖像分解成多尺度子帶選小波是因為它在 MATLAB 里工具鏈最完整、可逆性最好、邊界處理成熟換成曲波或剪切波理論上前景更好但在沒有現(xiàn)成工具箱的情況下自己寫重構核的代價太高工程上不劃算。第二塊是自適應收縮核我不用固定閾值而是讓閾值隨局部鄰域方差聯(lián)動并且給收縮曲線加了一個冪次因子來控制保邊強度。第三塊是參數(shù)搜索器用帶量子行為的粒子群算法在連續(xù)參數(shù)空間里搜索最優(yōu)組合。為什么搜參器選 QPSO 而不是量子遺傳算法或者網(wǎng)格搜索網(wǎng)格搜索在三維連續(xù)參數(shù)上就是災難每維取 10 個點三維就是 1000 次完整去噪單次算力約幾百毫秒加起來半小時起步還只能拿到粗糙的格點。遺傳算法需要把連續(xù)量離散成比特串編碼精度受位長限制而且雜交、變異算子對這類平滑性較好的目標函數(shù)并不占優(yōu)。QPSO 用的是勢阱模型的采樣方式位置更新里帶一個 $\ln(1/u)$ 的對數(shù)項天然具備跳出局部極值的能力同時參數(shù)是實數(shù)編碼不需要離散化再反解精度上限高得多。我實測下來在 3 到 5 維、每維取值范圍可控的情況下QPSO 大約 40 到 60 代就能收斂到一個相當穩(wěn)的解比 GA 穩(wěn)定比網(wǎng)格快一到兩個數(shù)量級。這里必須提前說清楚一件事避免誤解QAB 里的量子是算法思想層面的指的是用概率幅、疊加態(tài)采樣這類機制來設計更新規(guī)則不是真的跑到量子計算機上執(zhí)行。整套代碼在經(jīng)典 CPU 上用雙精度浮點跑沒有任何硬件門檻。很多初學者一看到量子兩個字就以為要裝什么特殊環(huán)境實際上它跟量子力學實驗沒有任何關系只是把一個基于概率分布的啟發(fā)式搜索換了套數(shù)學外衣。2. 量子自適應變換的核心機制拆解2.1 概率幅表示與勢阱采樣到底在做什么傳統(tǒng)粒子群里粒子位置是確定的速度也是確定的更新公式 $v \leftarrow \omega v c_1 r_1 (p-x) c_2 r_2 (g-x)$ 本質(zhì)上是帶上限的梯度下降。它的問題在于一旦所有粒子都擠在局部最優(yōu)附近速度項趨零整個種群就僵住了。量子行為的粒子群不用速度和位置這套描述改用波函數(shù)粒子不再有確定軌跡而是以某個概率密度分布在搜索空間里這個密度由勢阱中心決定。具體到代碼上更新式寫成$$x_{i,d} p_d \pm \alpha,|m_d - x_{i,d}|,\ln!\left(\frac{1}{u}\right)$$其中 $p_d$ 是粒子自身最優(yōu)與全局最優(yōu)的隨機加權用一個隨機數(shù) $\phi$ 混合$m_d$ 是全體個體最優(yōu)的均值也就是種群的質(zhì)心$\alpha$ 是收縮擴張系數(shù)$u$ 是 $(0,1)$ 上的均勻隨機數(shù)。那個 $\ln(1/u)$ 是這套算法最精髓的地方$u$ 接近 1 時它趨近 0粒子幾乎不動做精細開發(fā)$u$ 接近 0 時它是個大數(shù)粒子一步跳很遠做全局探索。而且它是重尾分布跳遠的概率雖然小但不是零這保證了即使種群暫時陷入某個洼地也有機會被甩出去。$\alpha$ 的作用是控制搜索半徑的整體縮放。我一般讓它從 1.0 線性衰減到 0.2前期大步掃后期小步精修。這個衰減節(jié)奏不能太快衰減過快會導致早熟也不能太慢慢了最后幾代還在大范圍亂跳收斂曲線抖得看不清趨勢。2.2 自適應基閾值必須隨子帶和位置一起變量子那部分只是找參的工具真正決定去噪質(zhì)量的是收縮核怎么設計。我把閾值寫成兩層結構。第一層是尺度相關第 $j$ 層細節(jié)子帶的閾值基準取$$T_j \lambda,\hat{\sigma}\sqrt{2\ln N_j}$$其中 $N_j$ 是該子帶的系數(shù)個數(shù)$\hat{\sigma}$ 是噪聲標準差估計$\lambda$ 是一個全局縮放系數(shù)交給搜索器去調(diào)。為什么用 $N_j$ 而不是總像素數(shù) $N$因為子帶系數(shù)個數(shù)隨層數(shù)遞減低層子帶系數(shù)多、噪聲能量分散算出來的閾值相對高一些高層子帶系數(shù)少閾值自動降下來避免把結構信息當噪聲削掉。這個改動看著小實測對紋理豐富的圖像能穩(wěn)定貢獻接近 0.2 到 0.4 dB 的提升。第二層是空間相關在子帶內(nèi)部再按 3×3 鄰域的局部標準差對閾值做微調(diào)$$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 就完全跟隨局部方差。這個式子的直覺是局部方差明顯高于全局噪聲水平的區(qū)域大概率是有結構的邊緣紋理應當減弱收縮力度保護細節(jié)局部方差接近噪聲水平的區(qū)域就是純噪聲應當加強收縮。注意閾值調(diào)的是基準的高低不是直接決定保留多少所以 $\kappa$ 過大也會出問題——局部方差估計本身帶噪$\kappa$ 太大會把噪聲估計的抖動放大成閾值的抖動反而制造塊狀偽影。實操中我把 $\kappa$ 限制在 $[0, 0.8]$。2.3 冪次軟閾值解決保邊和去噪的矛盾標準軟閾值的問題是它把所有超閾值系數(shù)統(tǒng)一平移 $T$幅值越大的系數(shù)被壓得越多相對比例上邊緣就軟了。我用的改進形式是把收縮寫成$$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)$兩者量綱完全一致不會出現(xiàn)幅值小于 1 的系數(shù)被取冪后反而變大這種荒謬情況。$\gamma1$ 時超閾值系數(shù)被相對抬升邊緣更銳$\gamma1$ 時壓得更狠噪聲殘留更少但圖像更糊。這個參數(shù)交給搜索器去定比人工試快得多。順帶說一句邊界處理。做局部方差估計必須用conv2滑窗而conv2的same模式在圖像四邊會用零填充等于人為制造了一圈暗邊局部方差在邊緣處會虛高。解決辦法是先對系數(shù)矩陣做對稱延拓padarray(..., symmetric)卷積完再裁回來。這個細節(jié)不注意去噪結果的四周會有一圈明顯異常尤其在大圖上特別扎眼。3. 關鍵參數(shù)的計算過程與取值范圍3.1 噪聲標準差的三條估計路徑$\hat{\sigma}$ 估不準后面全是白搭。最常用的是 MAD 估計取最細尺度對角細節(jié)子帶 $cD_1$ 的中位數(shù)絕對值除以 0.6745。這個 0.6745 是標準正態(tài)分布的中位數(shù)絕對偏差與標準差的比值來源是 $\Phi^{-1}(0.75)\approx 0.6745$。為什么要用中位數(shù)而不是均值因為中位數(shù)對離群值穩(wěn)健即使子帶里混進了少量強邊緣系數(shù)估計值也不會被帶跑偏。function sigma estimate_sigma(y, wname) [~, ~, ~, cD1] dwt2(y, wname); sigma median(abs(cD1(:))) / 0.6745; end第二條路徑是分塊估計把圖像切成 32×32 的塊對每塊單獨做 MAD 估計然后取所有塊估計值的中位數(shù)或 10% 分位數(shù)。分位數(shù)比中位數(shù)更保守因為它偏向取噪聲小的塊能避免把紋理塊的高方差誤算成噪聲。第三條是局部方差最小法找方差最小的那一批塊認為它們近似平坦用這些塊的方差做噪聲估計。后兩條在噪聲非均勻的場景下更可靠代價是要多做一輪遍歷。注意如果圖像本身經(jīng)過了 JPEG 壓縮或者做過銳化最細尺度子帶里會殘留壓縮塊效應和銳化振鈴MAD 估計會明顯偏高。這種情況下建議先在稍粗的尺度上估或者直接對 $\lambda$ 的搜索下界放寬。3.2 搜索空間的設定邏輯搜索空間開得太大收斂慢且容易搜到退化解開得太小搜索器就形同虛設。我一般按下面的范圍設定這幾個上下界不是拍腦袋的是反復試出來的參數(shù)含義下界上界定界理由$\lambda$閾值全局縮放0.43.0低于 0.4 等于沒去噪高于 3.0 會削掉大量結構$\gamma$收縮曲線冪次0.71.3越界后要么振鈴明顯要么過度模糊$\kappa$局部自適應力度0.00.8超過 0.8 后局部估計噪聲被放大level小波分解層數(shù)24見下文說明分解層數(shù)這一維有個特殊處理它必須是整數(shù)。如果直接讓 QPSO 在 $[2,4]$ 的連續(xù)區(qū)間上搜再四舍五入會出現(xiàn)大片的平臺區(qū)——目標函數(shù)在 2.3 到 2.7 之間完全相同搜索器拿不到有效梯度信息白白浪費評估次數(shù)。我的做法是把層數(shù)提到外層做成網(wǎng)格枚舉內(nèi)層只搜三個連續(xù)參數(shù)最后比較三個層數(shù)各自的最優(yōu)結果。雖然評估次數(shù)變成三倍但總次數(shù)反而比在四維連續(xù)空間里瞎搜少而且結果可解釋性強得多。一般圖像邊長在 512 左右時3 層是最常見的贏家邊長 256 以內(nèi)時 2 層往往更好因為再往下分第 3 層的子帶只剩幾十個系數(shù)統(tǒng)計意義不足。3.3 適應度函數(shù)怎么設計才不過擬合最直覺的做法是拿 PSNR 當適應度。但這里有個陷阱如果整張圖都參與搜索搜出來的參數(shù)會把這幅圖特有的噪聲實現(xiàn)記住換一張圖立刻失效這叫對噪聲過擬合。我的處理是把訓練圖切成互不重疊的兩半一半用于搜索參數(shù)另一半用于驗證適應度也改成兩部分的加權$$f -\left(0.6\cdot\text{PSNR}{\text{train}} 0.4\cdot\text{PSNR}{\text{val}}\right)$$取負號是因為 QPSO 是最小化框架。更嚴格的做法是留出第三份做測試但圖像規(guī)模有限時兩份已經(jīng)夠用。另外我提醒一點如果你手上的數(shù)據(jù)根本沒有干凈參考圖那 PSNR 這條路徑直接斷了只能退而求其次用盲指標——比如估計殘差圖的噪聲水平是否接近目標值或者用局部方差的一致性作為代理。這時候搜索出來的參數(shù)可信度會明顯下降不要把結果當成金標準。4. MATLAB 實現(xiàn)逐模塊落地4.1 環(huán)境準備與目錄組織工具箱只需要兩個Image Processing Toolbox讀圖、padarray、imresize和小波工具箱dwt2、idwt2、wavedec2。MATLAB 版本建議 2019b 及以上主要是為了arguments塊的兼容性實際上老版本也能跑把參數(shù)校驗改成手動判斷即可。目錄我習慣這么分qab_denoise/ ├─ main_demo.m % 主腳本 ├─ qab_denoise.m % 去噪主函數(shù) ├─ qab_shrink.m % 收縮核 ├─ qab_fitness.m % 適應度封裝 ├─ qab_qpso.m % 搜索器 ├─ estimate_sigma.m % 噪聲估計 ├─ utils/my_psnr.m % 指標 └─ data/ % 測試圖這種拆法的好處是每個文件職責單一改收縮核的時候不用翻主函數(shù)調(diào)試搜索器的時候可以套一個假的目標函數(shù)比如二維球形函數(shù)單獨驗證收斂性。我強烈建議先用一個已知最優(yōu)點的測試函數(shù)把qab_qpso.m跑通再接到去噪上不然一旦結果不好你根本分不清是搜索器有問題還是收縮核有問題。4.2 去噪主函數(shù)核心思路是逐層dwt2分解記錄每層尺寸用于重構時對齊然后逆向逐層idwt2。function [xHat, info] qab_denoise(y, par) % y : 含噪圖像double 類型范圍 [0,1] % par : 參數(shù)結構體字段 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)); % 尺寸對齊避免奇數(shù)邊長誤差累積 end xHat min(max(A, 0), 1); info.sigma sigmaN; info.level level; end這里有個坑值得單說dwt2對奇數(shù)尺寸的輸入輸出會做上取整idwt2還原時可能多出一行一列。如果不在每層還原后立即按原始尺寸裁剪誤差會一層層累積三層之后可能差出好幾個像素后面算 PSNR 直接對不齊報錯。我第一次寫的時候就栽在這上面調(diào)試了半天以為是閾值算法寫錯了。小波基我默認用db4。sym4和db4在去噪任務上表現(xiàn)接近db4支撐長度稍短、計算快一點sym4對稱性更好、相位失真小一點。如果你的圖像里有大量接近對稱的結構建筑、機械零件可以試試sym4。別用haar雖然它快但分塊感太強重構后能看到明顯的方格。也別盲目上更高階的db8支撐變長意味著邊界效應擴散范圍更大對中等尺寸以下的圖像反而不利。4.3 收縮核的完整實現(xiàn)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兩個保護上下限是必須的。下限防止某些局部方差極低的區(qū)域閾值趨零那些區(qū)域的噪聲會被完整保留出現(xiàn)隨機亮斑上限防止強紋理邊緣把閾值抬得過高導致邊緣兩側被過度收縮形成光暈。這兩個界限看著粗暴但實測能顯著提升結果的視覺穩(wěn)定性。另外eps的位置要注意放在分母上是為了防止除零但如果Tloc本身很小eps不足以影響結果如果Tloc正常eps完全不起作用。不用擔心它引入偏差。4.4 QPSO 搜索器實現(xiàn)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參數(shù)取值上種群規(guī)模 20 到 30 就夠再大邊際收益很低迭代 60 到 100 代具體看收斂曲線的拐點。alpha0取 1.0alphaDecay取 0.97 到 0.98意味著 60 代后衰減到大約 0.2 到 0.3。我試過固定alpha不衰減結果最后十幾代還在大范圍跳動最優(yōu)值一直在變曲線沒法收斂。fitFun用匿名函數(shù)封裝把含噪圖和干凈參考圖綁進去fitFun (p) -qab_fitness(yTrain, p, par0, xTrainClean);qab_fitness里把連續(xù)參數(shù)映射過去注意層數(shù)是外層固定的所以這里只處理三個連續(xù)量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(迭代代數(shù)); ylabel(最優(yōu)適應度PSNR); grid on;跑一遍大概需要一到三分鐘取決于圖像尺寸和層數(shù)枚舉的個數(shù)。如果嫌慢可以把種群降到 16、迭代降到 40代價是最終 PSNR 可能掉 0.1 到 0.2 dB。5. 實驗結果與參數(shù)敏感性觀察5.1 測試構造與評價口徑評估這塊我不喜歡只報一個 PSNR太片面。我一般同時看三個量PSNR 反映整體誤差能量SSIM 反映結構相似性還有一個平坦區(qū)殘差標準差——在手工標注的若干塊純背景區(qū)域上計算去噪結果與原圖之差的標準差這個值能直接暴露噪聲有沒有削干凈。很多時候 PSNR 漲了但平坦區(qū)還是一片噪點說明算法把力氣都用去保細節(jié)了這種結果放到實際產(chǎn)品里是過不了驗收的。測試集構造上我一般選十張左右覆蓋不同內(nèi)容類型的圖一張大面積平滑的天空、墻面一張強紋理的草地、織物一張多邊緣的建筑、文字再加若干張自然場景圖作為綜合參考。噪聲類型默認加性高斯白噪聲$\sigma$ 取 10、20、30、50 四檔對應 $\sigma/255$ 的歸一化值。如果要測真實噪聲建議找?guī)О倒饣蚋?ISO 場景的圖但那種圖很難拿到干凈參考只能靠目視加盲指標。5.2 不同噪聲水平下的表現(xiàn)趨勢下面這組是量級參考具體數(shù)值隨圖像內(nèi)容、隨機種子、小波基都會浮動別拿它當硬指標對比。噪聲 $\sigma$最優(yōu) $\lambda$ 落點最優(yōu) $\gamma$ 落點最優(yōu) $\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這張表透露了一個很明確的規(guī)律噪聲越強最優(yōu)閾值越高收縮也應該越狠同時局部自適應的力度要更大。這完全符合直覺——噪聲大時子帶里噪聲和細節(jié)的重疊更嚴重需要更激進的閾值把噪聲推下去而噪聲強的時候局部方差估計相對更可靠信噪比低了但方差主要由噪聲貢獻所以可以放心加大 $\kappa$。反過來在低噪聲下$\gamma$ 和 $\kappa$ 都被壓到偏小的值說明這時候主要靠精細的幅度保留取勝過度收縮反而傷害細節(jié)。還有一個值得注意的現(xiàn)象$\lambda$ 和 $\gamma$ 之間存在明顯的負相關。$\lambda$ 增大時最優(yōu) $\gamma$ 傾向于也增大兩者在效果上具有補償關系。這意味著搜索空間里存在一條山谷沿谷底走目標函數(shù)變化很小這也是為什么網(wǎng)格搜索在這類問題上效率極低——大量評估花在了重復的谷底點上。5.3 消融實驗量子自適應到底貢獻了什么拆開看各部分的貢獻思路會更清楚。我設計了幾組對照第一組是純標準軟閾值$\lambda$ 用經(jīng)典公式算、$\gamma1$、$\kappa0$第二組開尺度自適應$\lambda$ 靠人工試三檔取最優(yōu)第三組加局部自適應$\kappa$ 人工試三檔第四組全開參數(shù)由 QPSO 搜索。配置相對基線提升主要改善區(qū)域標準軟閾值基線— 尺度自適應 $\lambda$0.4 ~ 0.6 dB整體 局部自適應 $\kappa$再 0.3 ~ 0.5 dB紋理區(qū)、邊緣 冪次收縮 $\gamma$再 0.1 ~ 0.3 dB邊緣銳度目視更明顯 QPSO 自動搜參再 0.2 ~ 0.4 dB各參數(shù)組合的協(xié)同從這個拆解能看出來收益的大頭來自尺度自適應 局部自適應這兩個結構性改動搜索器只是把剩下那點參數(shù)協(xié)同的空間榨出來。這點認知很重要——如果你時間有限先把局部自適應做好收益比死磕優(yōu)化器大得多。我見過一些方案把 90% 的篇幅放在優(yōu)化器上收縮核卻還是標準軟閾值那屬于用力的方向偏了。視覺上最直觀的差異在紋理區(qū)。標準軟閾值處理過的草地、織物紋理會被糊成一片均勻的灰細節(jié)全部丟失加上 $\kappa$ 之后紋理區(qū)按局部方差自動降低了閾值毛刺感被保留下來同時平坦區(qū)的噪點被更強的閾值壓下去同一幅圖里出現(xiàn)了該軟的地方軟、該硬的地方硬的效果。這種區(qū)域性差異是任何全局閾值都做不到的。6. 常見問題排查實錄6.1 PSNR 反降的幾種典型成因第一種噪聲估計偏高。如果圖像有壓縮痕跡或者做過銳化MAD 估計會把非噪聲的高頻成分算進去導致 $T_j$ 整體偏大把所有細節(jié)都削平。判斷方法是把估計出來的 $\hat{\sigma}$ 和真實加噪強度對比偏差超過 20% 就要警惕。解決辦法是改用分塊分位數(shù)估計或者對 $\hat{\sigma}$ 乘一個 0.8 左右的修正因子。第二種層數(shù)選太多。小波分解每深一層低頻近似就被不斷下采樣到了第 4、5 層一個系數(shù)對應原圖 16×16 甚至 32×32 的區(qū)域這個尺度上已經(jīng)談不上細節(jié)和噪聲的區(qū)分了強行收縮只會引入塊狀偽影。如果發(fā)現(xiàn)結果里出現(xiàn)了大塊的方格基本就是這個原因。判斷起來很簡單把層數(shù)逐次降低重跑一次看偽影是否消失。第三種局部自適應參數(shù)過大。$\kappa$ 超過 0.9 之后閾值會緊密跟隨局部方差而局部方差本身是從含噪系數(shù)算出來的本身就帶噪噪聲被放大回閾值里形成一種噪聲調(diào)制閾值、閾值放大噪聲的正反饋。表現(xiàn)是圖像上出現(xiàn)細密的隨機紋理看著更臟了。把 $\kappa$ 壓到 0.6 以下基本能解決。第四種尺寸對齊沒做。前面提過idwt2的尺寸和原圖差一個像素最后算指標時對不齊或者被imresize強行縮放PSNR 直接掉好幾個 dB。這種情況的特點是斷崖式下跌一眼就能看出來檢查一下重構循環(huán)里的裁剪語句就行。6.2 運行慢、內(nèi)存漲的排查路徑慢的原因九成在搜索循環(huán)里。每次評估都要跑一遍完整去噪包括兩層到四層小波分解和重構這是實打實的開銷。優(yōu)化手段有幾個按性價比排序第一降低搜索規(guī)模種群從 30 降到 20迭代從 80 降到 50評估次數(shù)直接砍一半多PSNR 損失通常在 0.1 dB 以內(nèi)。第二在縮略圖上搜參把圖降采樣到 256×256 上搜出參數(shù)再套回原尺寸圖因為最優(yōu)參數(shù)對尺度不敏感閾值公式里已經(jīng)含了 $N_j$ 的補償這招能省下 80% 以上的時間。第三并行化適應度評估parfor直接套在最內(nèi)層循環(huán)上四核能快三倍左右注意循環(huán)內(nèi)的變量不能有交叉依賴。內(nèi)存漲的典型場景是參數(shù)沒有及時清理。details這種 cell 數(shù)組在循環(huán)里被反復覆蓋MATLAB 的寫時復制機制在數(shù)據(jù)被別處引用時會保留舊副本。如果發(fā)現(xiàn)跑幾十輪之后內(nèi)存居高不下在每輪循環(huán)末尾顯式clear一下大變量或者把整個搜索過程包在一個函數(shù)里靠作用域自動回收。6.3 問題速查表現(xiàn)象最可能的成因快速驗證處理方式結果四周有一圈異常卷積邊界零填充放大看四邊padarray對稱延拓邊緣明顯發(fā)虛$\gamma$ 偏大把 $\gamma$ 手動設 0.8 重跑收緊 $\gamma$ 上界平坦區(qū)仍有噪點$\kappa$ 偏大或 $\lambda$ 偏小兩者分別單獨調(diào)整降 $\kappa$、升 $\lambda$出現(xiàn)塊狀偽影層數(shù)過多層數(shù)降到 2 重跑限制層數(shù)或換sym4搜索曲線一直抖alpha衰減太慢看最后十代是否有下降加快衰減或提前終止換圖后效果驟降對噪聲實現(xiàn)過擬合用另一張圖的驗證集測拆分訓練/驗證集運行報維度不匹配idwt2尺寸誤差打印每層size逐層裁剪對齊實操心得調(diào)試這套東西的順序很重要我建議按先固定參數(shù)跑通流程、再驗證搜索器、最后聯(lián)合調(diào)優(yōu)三步走。一上來就讓搜索器在沒有驗證過的主函數(shù)上跑出問題你根本定位不到是搜索器還是收縮核純屬浪費時間。7. 工程落地前的幾點提醒7.1 離線搜參、在線固定別在線上跑優(yōu)化器這一點我特別想強調(diào)因為它直接決定這套方案能不能上生產(chǎn)。QPSO 搜一次參數(shù)要幾十秒到幾分鐘而實際部署場景里一張圖進來你必須在一秒內(nèi)出結果甚至要去到幾十毫秒。所以正確的用法是在離線階段用一批有代表性的圖對每一檔噪聲水平、每一類內(nèi)容平滑型、紋理型、邊緣型各搜出一組最優(yōu)參數(shù)做成一張查找表存下來。在線階段只做一件事——估計輸入圖的噪聲水平然后從表里查最近的參數(shù)直接跑一次去噪。這條路走通的代價是泛化性的妥協(xié)查表參數(shù)一定不如現(xiàn)搜的最優(yōu)實測差距在 0.2 到 0.5 dB 之間。但從工程角度看這個代價換來的速度提升是幾百倍完全值得。而且只要你離線搜參時用的圖覆蓋足夠廣在線效果不會崩。還有個小技巧噪聲水平不需要估得太準按 10、20、30、50 分成幾檔實際值落到哪一檔就用哪一檔的參數(shù)中間地帶用線性插值。因為參數(shù)對噪聲水平的敏感度是連續(xù)的檔位之間差 5 個單位帶來的 PSNR 損失通常不到 0.1 dB。7.2 這套框架還能往哪兒長第一個方向是換變換。小波在處理各向異性邊緣斜線、曲線時效率不高換成剪切波或者非下采樣輪廓波邊緣表示能力會明顯更強尤其對細長的曲線結構。難點是重構核得自己實現(xiàn)需要仔細驗證完美重構條件我建議先在一個可逆性測試上跑通再接到去噪里。第二個方向是換成有監(jiān)督的參數(shù)預測?，F(xiàn)在參數(shù)是搜出來的能不能訓一個小網(wǎng)絡輸入含噪圖直接輸出四個參數(shù)這樣在線階段連估計噪聲那一步都省了。難度在于訓練數(shù)據(jù)的構造——你需要為每張訓練圖找到真正的最優(yōu)參數(shù)作為標簽這個標注成本很高但可以用現(xiàn)在的 QAB 流程離線生成正好形成一個閉環(huán)QAB 離線造標簽網(wǎng)絡在線快速預測。第三個方向是把收縮核換成非局部機制?，F(xiàn)在用的是 3×3 鄰域本質(zhì)上還是局部的。如果換成塊匹配、相似塊分組在組內(nèi)做協(xié)同收縮對重復紋理磚墻、織物、柵欄的去噪效果會有質(zhì)的提升。代價是計算量上去一到兩個數(shù)量級得配合積分圖或者快速傅里葉變換做加速。這條路的性價比很高尤其是如果你的應用場景里重復紋理多。第四個方向是處理彩色圖像?，F(xiàn)在整套流程是按單通道灰度寫的直接套到 RGB 三個通道上分別處理會有色度失真問題——三個通道的噪聲實現(xiàn)是獨立的分別去噪后顏色會偏移。正確做法是轉換到 YCbCr 或者小波域里做通道間聯(lián)合收縮對色度通道用更低的閾值人眼對色度噪聲不敏感但對色偏很敏感。我在這套東西上折騰了挺長時間最大的一個體會是優(yōu)化器再花哨也補不上收縮核設計的缺陷。早期我把大量精力花在調(diào) QPSO 的種群規(guī)模、衰減系數(shù)、隨機種子結果 PSNR 在一個范圍內(nèi)來回晃后來把 3×3 局部方差窗口改成 5×5把閾值下限保護加進去一次就往上跳了 0.3 dB。參數(shù)搜索能幫你挖到當前設計的天花板但它抬不高天花板本身。所以真要往深里做先把收縮核想清楚再考慮優(yōu)化器要不要換更復雜的。另外再分享一個調(diào)試時管用的小習慣把每次去噪的殘差圖含噪圖減去去噪圖存下來看一眼。理想情況下殘差應該接近純白噪聲看不見任何結構。如果殘差里能看出明顯的邊緣輪廓說明閾值削過頭了把結構當噪聲削掉了如果殘差里能看到大塊的斑駁區(qū)域說明局部自適應在那些地方過于激進了。這一眼比看 PSNR 數(shù)字直觀得多很多時候 PSNR 漲了但殘差里有結構那就是在過擬合。