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