混合二階全變分的抗核輻射圖像降噪方法

黃小莉，陳春梅,2，劉桂華,2

(1.西南科技大學(xué) 信息工程學(xué)院，四川 綿陽 621000；2.特殊環(huán)境機(jī)器人技術(shù)四川省重點實驗室，四川 綿陽 621000)

0 引 言

核環(huán)境中，圖像采集器采集到的圖像含有大量核噪聲亮斑，加大了后續(xù)對核輻射場景感知的難度。目前，對于核輻射干擾圖像的降噪研究相對較少，對核環(huán)境下的圖像進(jìn)行降噪處理具有極高的學(xué)術(shù)研究價值和較強(qiáng)的創(chuàng)新性。

傳統(tǒng)的圖像降噪方法，如中值濾波法[1]、均值濾波法[2]等在進(jìn)行降噪時會對圖像的細(xì)節(jié)信息造成損失。文獻(xiàn)[3]首次提出全變分(total variation，TV)算法，很好地解決了傳統(tǒng)降噪方法造成的圖像邊緣模糊的問題。文獻(xiàn)[4]將分裂布雷格曼方法應(yīng)用于全變分模型中，解決了全變分模型的正則化問題。文獻(xiàn)[5-6]提出了基于鄰近算子的全變分方法(proximity approximation operator total variation，PATV)，通過線性變換計算凸函數(shù)組成的近似算子提升降噪性能。

為了緩解全變分降噪引起的階梯偽影問題，很多文獻(xiàn)進(jìn)行了大量研究。文獻(xiàn)[7]提出了廣義全變分(generalized total variation，TGV)，在一定程度上能抑制階梯偽影的產(chǎn)生。文獻(xiàn)[8]提出了一種自適應(yīng)TGV模型的降噪算法，通過自適應(yīng)參數(shù)使得邊緣和區(qū)域平滑。文獻(xiàn)[9]提出了一種重疊組稀疏全變分(overlapping group sparse total variation，OGSTV)方法，通過突出圖像平滑區(qū)域與邊緣區(qū)域的差異性，有效抑制了階梯偽影。文獻(xiàn)[10]將高階非凸正則化加入到OGSTV上，使圖像紋理更為平滑。文獻(xiàn)[11-12]提出的高階方程法在抑制階梯偽影的同時很好地保留了圖像邊緣信息。

為了解決高階全變分造成的圖像過度平滑問題，很多學(xué)者提出了混合技術(shù)。文獻(xiàn)[13]將一階全變分同二階全變分相結(jié)合，既保持了圖像的邊緣又抑制了階梯效應(yīng)。文獻(xiàn)[14]提出利用交替方向乘子法(alternating direction method of multiplier, ADMM)作為快速求解器有效地解決L1-TV正則化問題。文獻(xiàn)[15]提出一種針對脈沖噪聲去噪的帶L1保真項的混合變分模型，增強(qiáng)模型的去噪性能。

受到混合算法的啟發(fā)，本文提出一種混合去噪算法，將重疊組稀疏正則化與非凸二階全變分方法相結(jié)合，采用ADMM和增廣拉格朗日乘子法進(jìn)行交替循環(huán)求解，將降噪后的圖像進(jìn)行相減，在初始圖像上對得到的差值進(jìn)行疊加操作與迭代處理，優(yōu)化最終降噪效果。與對比算法的結(jié)果相比，本算法保留了更多的邊緣細(xì)節(jié)信息，能更大程度地還原圖像。

1 混合二階全變分

圖像降噪可以看成一個線性逆問題的求解過程，即從待處理圖像中推導(dǎo)出原圖像，常見的去噪模型為

I=Kg+γ

(1)

(1)式中：g∈RMN×1為未受噪聲污染的原始圖像；γ表示噪聲；K為點擴(kuò)散函數(shù)；I∈RMN×1為含有噪聲的待處理圖像。

求解圖像降噪的關(guān)鍵在于對(1)式進(jìn)行正則化處理，即將原始圖像g中的先驗知識融入線性問題的求解過程中，由此獲得具有正則化解，同時抑制噪聲。

1.1 二維重疊組稀疏正則化

傳統(tǒng)的組稀疏正則化器[16]方程表示為

(2)

圖像g的二維M×M點群表示為

(3)

(4)

1.2 重疊組稀疏正則化與二階全變分混合

將非凸二階全變分與重疊組稀疏正則化混合構(gòu)成本文提出方法的前半部分，可以簡單表示為

(5)

通過假設(shè)圖像周期性邊界條件，已知?g與?2g為一階和二階差分算子，引入輔助變量將(5)式轉(zhuǎn)換為約束優(yōu)化問題，表示為

(6)

1.3 ADMM和增廣拉格朗日乘子法優(yōu)化求解

ADMM是一種求解優(yōu)化問題的計算方法[19]，將受約束的全局問題分解為多個局部子問題。在本算法中，ADMM通過引入輔助變量將全局問題分裂為多個子問題，交替迭代優(yōu)化多個子問題的解從而得到全局問題的解。使用增廣拉格朗日函數(shù)和ADMM對(6)式進(jìn)行求解，增廣拉格朗日函數(shù)表示為

(7)

(8)

(8)式中，求解gm是一個最小二乘問題。

在圖像的周期邊界條件下，由于矩陣均為塊循環(huán)，因此，使用二維離散傅里葉變換對(8)式進(jìn)行求解，通過將時域乘法問題轉(zhuǎn)換為頻域問題，極大地減少計算復(fù)雜度，表示為

(9)

(9)式中，F(xiàn)-1為傅里葉逆變換。T的計算表達(dá)式為

T=

(10)

1.3.1 變量ζ1的優(yōu)化

變量ζ1的優(yōu)化是一個重疊組稀疏性問題?？梢员硎緸?/p>

(11)

使用優(yōu)化-最小化(MM)[20]方法解決重疊組稀疏的問題，在每次迭代中使用一個二次函數(shù)φ(ζ1)。在MM算法的每次迭代中，二階問題被最小化，(11)式的收斂式為

(12)

(12)式中，重疊組稀疏全變分正則化φ(ζ1)的定義為

(13)

(13)式中，M為重疊組大小。

1.3.2 變量ζ2的優(yōu)化

變量ζ2的優(yōu)化是一個非凸二階去噪問題，可以表示為

(14)

采用分離變量法將(14)式改寫為

(15)

采用迭代重加權(quán)算法(IRL1)[21]對(15)式進(jìn)行求解，每一次迭代，(15)式都可以表示為

(16)

(16)式中，每次迭代更新的權(quán)重t表達(dá)式為

(17)

(17)式中，δ是一個避免被除數(shù)為0的小數(shù)字。

利用收縮閾值函數(shù)對(16)式進(jìn)行計算，可得

(18)

(18)式中，shrink(·)表示一維收縮運算符。

1.3.3 變量ζ3的優(yōu)化

變量ζ3是一個確保像素值保持在[0,255]的投影，表示為

(19)

采用投影的方法來求解(19)式，使用更新規(guī)則，將迭代次數(shù)大大減小，最終的拉格朗日乘子迭代為

(20)

1.4 差值迭代

將上述混合非凸二階全變分算法降噪過程設(shè)為H，降噪后的圖像與原圖像進(jìn)行迭代降噪，多次操作后得到圖像的灰度值接近于原始圖像的灰度值，提高了降噪的效率。具體流程如下。

1)將圖像g通過上述提出的混合非凸二階全變分算法進(jìn)行降噪，得到處理后的圖像g1為

g1=H(g,1)

(21)

2)將處理后的圖像g1再次進(jìn)行混合非凸二階全變分降噪得到g2為

g2=H(g1,1)

(22)

將g1與g2相減得到gΔ，表示為

gΔ=g1-g2=H(g,1)-H(g1,1)

(23)

設(shè)定β為迭代次數(shù)的倒數(shù)，修改β的值來提高算法的降噪魯棒性。令β∈(0,1)，將差值gΔ進(jìn)行β倍的縮小并與原圖像相加，得到圖像g3，g3不僅包含有初始圖像的框架，還包含有前兩次降噪后的結(jié)果，得到的g3為

g3=g+β(g1-g2)=g+βgΔ

(24)

3)對g3進(jìn)行連續(xù)迭代降噪，設(shè)迭代次數(shù)為Ter，得到最終圖像gn為

gn=H(g3,Ter)

(25)

2 模擬核噪聲

2.1 核噪聲的特性

在核輻射場景下，將圖像采集鏡頭用鏡蓋蓋住，此時無可見光輸入，核輻射粒子直接作用于傳感器，生成噪聲亮斑圖像如圖1所示。從圖1a可以看出，核噪聲為多個塊狀的亮斑，每個亮斑的像素值相近。增強(qiáng)核輻射強(qiáng)度，得到圖1b所示圖像，對比圖1a、圖1b能看出，當(dāng)核輻射強(qiáng)度增強(qiáng)時，放射性粒子產(chǎn)生的脈沖電流的數(shù)量與強(qiáng)度也隨之增加，同時圖像上形成的噪聲亮斑的數(shù)量和面積也會隨之增加。由圖1c和圖1d得出，灰度值的幅值服從高斯分布且灰度級在50～150出現(xiàn)的概率最大。

圖1 無可見光下不同輻射強(qiáng)度采集的圖像及對應(yīng)灰度直方圖

2.2 核噪聲的模擬

根據(jù)核噪聲的隨機(jī)性，采用隨機(jī)行走方法來模擬核噪聲斑塊的形狀。將圖像中的像素點比作核輻射粒子，隨機(jī)選擇一個像素點作為隨機(jī)行走的起點，以起點為中心向周圍8領(lǐng)域以概率P進(jìn)行擴(kuò)散。通過調(diào)整行走概率P和迭代次數(shù)來調(diào)整模擬的噪聲斑塊面積，噪聲斑塊的數(shù)量則由斑塊的分布率決定。將隨機(jī)行走得到的矩陣與高斯矩陣相乘得到模擬核噪聲斑塊圖。

綜上分析，模擬核噪聲斑塊的具體流程如下。

1)生成一個只含有0和1的二值隨機(jī)矩陣，矩陣中1的個數(shù)為泊松函數(shù)得到的泊松隨機(jī)數(shù)，將生成的泊松隨機(jī)數(shù)作為模擬生成的核噪聲個數(shù)；

2)以矩陣中的1為起點向周圍8領(lǐng)域以概率P進(jìn)行隨機(jī)行走，通過調(diào)整概率P和迭代次數(shù)來調(diào)整模擬的噪聲面積，接下來設(shè)置斑塊的分布率-斑塊噪聲的像素個數(shù)與圖像中所有像素個數(shù)的比值，從而模擬得到不同的斑塊數(shù)量與不同輻射強(qiáng)度中的噪聲斑塊；

3)將隨機(jī)行走得到的矩陣與高斯矩陣相乘得到最終模擬的核噪聲斑塊；

4)將模擬得到的核噪聲斑塊加入到不同場景的自然圖像中，得到不同場景下的核輻射噪聲圖像。

圖2所示為不同參數(shù)模擬的不同強(qiáng)度核噪聲斑塊，設(shè)定圖2a中的模擬弱核斑塊行走概率P為0.1，迭代次數(shù)為5，斑塊分布率為0.3；圖2b中的模擬中度核斑塊行走概率P為0.1，迭代次數(shù)為10，斑塊分布率為0.5；圖3c中的模擬輻射較強(qiáng)時的核斑塊行走概率P為0.1，迭代次數(shù)為25，斑塊分布率為0.8。由核噪聲特性可知，隨著核輻射強(qiáng)度的增加，產(chǎn)生的核輻射斑塊噪聲的面積變大，斑塊的密度也隨之增加。

圖2 不同參數(shù)模擬的不同強(qiáng)度核噪聲斑塊

將模擬的核斑塊加入到數(shù)字圖像處理常用的數(shù)據(jù)集set12中，圖3為其中一張圖像的效果。本文通過模擬得到不同場景下的核輻射噪聲圖片，為后續(xù)實驗提供條件。

圖3 加入不同強(qiáng)度模擬噪聲后的圖像

3 實驗參數(shù)的選取

3.1 差值迭代次數(shù)的確定

在實驗過程中需要設(shè)置差值迭代次數(shù)Ter，迭代次數(shù)過大會消耗過多的計算時間，迭代次數(shù)過小則無法達(dá)到最佳降噪效果。用不同迭代次數(shù)對獲得的模擬數(shù)據(jù)集進(jìn)行實驗，以峰值信噪比(peak signal-to-noise ratio, PSNR)和結(jié)構(gòu)相似性(structural similarity, SSIM)為指標(biāo)來獲得最佳迭代次數(shù)。圖4為不同迭代次數(shù)處理后的PSNR和SSIM折線圖。通過圖4看出，差值迭代次數(shù)Ter為15時，PSNR和SSIM數(shù)值達(dá)到最高，繼續(xù)增加迭代次數(shù)，PSNR和SSIM無明顯變化。

圖4 不同迭代次數(shù)處理后的PSNR和SSIM折線圖

3.2 重疊組M大小的確定

在實驗過程中需要設(shè)置重疊組M的大小。固定其他參數(shù)，選取模擬輻射強(qiáng)度為弱的數(shù)據(jù)集進(jìn)行實驗，以PSNR和SSIM為指標(biāo)來獲得最佳重疊組M的大小。圖5為不同重疊組M處理后的PSNR和SSIM折線圖。從圖5可知，當(dāng)M=4時，PSNR和SSIM數(shù)值最大。由此，設(shè)置重疊組M的大小為4。

圖5 不同重疊組M處理后的PSNR和SSIM折線圖

3.3 其他參數(shù)的調(diào)整

首先，對變量λ的范圍進(jìn)行調(diào)整。變量λ控制數(shù)據(jù)保真項，在圖像去噪上具有重要作用。由于在實驗中使用不同模擬噪聲強(qiáng)度的數(shù)據(jù)集，因此，根據(jù)經(jīng)驗設(shè)置λ∈[1×10-4，4]。

然后，對混合算法的內(nèi)部迭代次數(shù)進(jìn)行調(diào)整。設(shè)置差值迭代次數(shù)分別為1、5、10、15、20，對不同情況下的PSNR進(jìn)行對比，以此得到最佳的內(nèi)部迭代次數(shù)。選取模擬輻射強(qiáng)度為弱的數(shù)據(jù)集“house”和“woman”進(jìn)行實驗，結(jié)果如表1所示。當(dāng)?shù)螖?shù)為10的時候，PSNR數(shù)值達(dá)到最大值。由于迭代次數(shù)過高會導(dǎo)致計算時間增加，綜合考慮，選取內(nèi)部迭代最佳次數(shù)為10。

表1 不同迭代次數(shù)的降噪結(jié)果

4 實驗結(jié)果

為了驗證算法的有效性和可靠性，在兩類不同數(shù)據(jù)集上進(jìn)行實驗，數(shù)據(jù)集分為模擬得到的核噪聲數(shù)據(jù)集及真實核環(huán)境下的數(shù)據(jù)集。本算法與傳統(tǒng)中值濾波、鄰近算子全變分方法(PATV)[5]、廣義全變分方法(TGV)[8]、重疊組稀疏的全變分方法(OGSTV)[9]、L1全變分方法[14]等多個算法進(jìn)行比較，采用峰值信噪比、結(jié)構(gòu)相似性以及圖像粗糙度作為客觀評判算法性能的指標(biāo)，同時從主觀角度出發(fā)，通過圖像細(xì)節(jié)部分對圖像最終處理結(jié)果進(jìn)行評定。本實驗在處理器為AMD Ryzen 72700X Eight-Core Processor at 3.70 GHz，內(nèi)存為16.0 GByte，系統(tǒng)為Windows 10專業(yè)版，MATLAB為R2015b的環(huán)境下進(jìn)行實驗。

4.1 對模擬得到的圖像進(jìn)行實驗

對獲得的多張不同輻射強(qiáng)度的圖片進(jìn)行處理，表2為不同去噪算法的結(jié)果對比。由表2可見，在不同的模擬核輻射強(qiáng)度中，本文算法得到的PSNR和SSIM值最高。將不同算法得到的最終結(jié)果的紅框處進(jìn)行放大，如圖6—圖7所示。由圖6—圖7能看出，本算法在降噪的同時能更好地保留圖像的邊緣信息。

表2 不同去噪算法的結(jié)果對比

圖6 “peppers”處理后的局部放大圖

圖7 “woman”處理后的局部放大圖

4.2 對真實核環(huán)境下的核圖像進(jìn)行實驗

為了進(jìn)一步驗證算法的有效性，對真實核環(huán)境下采集到的噪聲圖像進(jìn)行降噪處理，同時將處理結(jié)果與不同算法進(jìn)行對比。

圖8為不同算法對核環(huán)境下噪聲圖片進(jìn)行處理后的結(jié)果。圖8a為廣州番禺事故及河南杞縣事故處理時在真實核環(huán)境下采集到的圖像，整個圖像被噪聲包圍且能見度極低，嚴(yán)重影響了核環(huán)境的識別及后續(xù)操作。將圖8b—圖8g紅框處部分進(jìn)行放大，可得到“鐵柱”處理后的局部放大圖如圖9所示，得到“小車”處理后的局部放大圖如圖10所示。從圖9可知，本算法結(jié)果的細(xì)節(jié)處更加清晰，引起的圖像缺陷也更少，裂縫銜接處紋理更加流暢。從圖10可知，本算法比其他算法的處理結(jié)果損失更少，結(jié)果中含有的噪聲含量較少，同時物體邊緣更加清晰。

圖8 不同算法對核環(huán)境下噪聲圖片進(jìn)行處理的結(jié)果

圖9 “鐵柱”處理后的局部放大圖

圖10 “小車”處理后的局部放大圖

5 結(jié)束語

本文通過將非凸二階全變分與重疊組稀疏結(jié)合起來，在降噪的同時抑制階梯效應(yīng)，從而極大地保留了圖像的細(xì)節(jié)信息。在此混合算法的基礎(chǔ)上，將處理后得到的圖像進(jìn)行差值迭代，使得圖像的灰度更加偏向原圖灰度，進(jìn)一步提升了圖像的質(zhì)量。實驗結(jié)果表明，本混合算法在視覺主觀效果上和客觀數(shù)據(jù)評價指標(biāo)上優(yōu)于對比算法，證實了算法的優(yōu)越性和可行性。