楊 華

(南充職業(yè)技術學院 電子信息工程系，四川 南充 637000)

1 引 言

數字圖像在成像、壓縮和傳輸過程中經常會引入噪聲，進而導致圖像可視程度降低和信息的丟失[1-2]。為了提高圖像成像質量，圖像處理領域的研究者提出了多種算法[3-6]。對很多圖像處理算法來說，在進行圖像處理時，需預先估計圖像噪聲的分布模型和噪聲的相關參數。

在大多數情況下，圖像噪聲服從高斯分布，為高斯白噪聲。高斯白噪聲一般為零均值。因此圖像噪聲估計參數一般為標準差或者方差。常用的噪聲方差估計方法包括奇異值分解[7]、主成分分析[8]和小波變換[9]等。圖像噪聲來源主要包括傳感器噪聲、電路失穩(wěn)、傳輸過程中的壓縮等[10]。Li等人[11]認為圖像中的信息可以被視為能量，被噪聲污染之后的圖像信息的部分損失可以歸結為能量損失。稀疏主成分分析對數字圖像這類矩陣數據的能量變化尤為敏感。在實驗中，我們觀察到噪聲圖像的高斯白噪聲標準差等級與稀疏主成分分析得到的部分主成分負載向量均值呈現較為良好的線性關系。基于該規(guī)律，本文提出了一種基于稀疏主成分分析的圖像高斯噪聲估計方法。該方法計算復雜度低且具有較強的魯棒性和精確性。

2 稀疏主成分分析

2.1 主成分分析

主成分分析(Principal Component Analysis, PCA)目的是在將每個樣本屬性值數量極大程度減小的前提下，保證原始數據矩陣的原始信息得以最大程度保留。

求解主成分(Principal Component)主要有兩種方法[12]。

第一種方法利用了協(xié)方差矩陣的特性：

設有包含n個樣本、p個屬性值的矩陣Xn×p。

將矩陣進行中心化和標準化，標準化即每一個屬性向量內的元素除以該向量內所有元素的標準差。其協(xié)方差矩陣為：

(1)

可以發(fā)現，協(xié)方差矩陣的對角線依次排列各屬性方差，非對角線上的元素是所有屬性兩兩之間的協(xié)方差。因此使用主成分分析得到的降維后的低維矩陣的協(xié)方差矩陣必須滿足這兩個條件：對角線元素盡可能大、非對角線元素為0。

設使用PCA對矩陣X降維得到矩陣Z，X、Z的協(xié)方差矩陣分別為C、D,P為一組正交基，有Z=XP。則有：

(2)

因此，優(yōu)化目標變?yōu)閷ふ沂沟迷紨祿仃嘪的協(xié)方差矩陣C對角化且對角線元素從上到下降序排列的矩陣P。由于協(xié)方差矩陣為實對稱矩陣，在線性代數上，實對稱矩陣有一系列極有實用價值的性質[13]。若某實對稱矩陣大小為p行p列，則一定可以找到p個互相正交的單位特征向量，設這p個特征向量組成矩陣E：

E=(e1,e2,…,ep)，

則對協(xié)方差矩陣C有：

取 1 mL 樣品溶液加入5 mL 試劑甲混勻，于20～25 ℃放置10 min，再加0.5 mL試劑乙(Folin-酚)，立即搖勻，記錄反應總體積為6.5 mL。在 20～25 ℃保溫30 min，然后于500 nm 處比色，以 0 mL標準蛋白代替作空白對照。

(3)

若將特征值向量按照降序排列，每個特征值對應的特征向量矩陣E也進行對應的重排。則式(3)得到的對角陣從上到下各元素的大小也為降序排列。

因此，主成分分析所求的P就是原始數據矩陣X的協(xié)方差矩陣C的特征向量E。按照特征值大小取前K個特征值對應的特征向量按列排序，構成PCA前k個主成分對應的負載矩陣βp×k，前K個主成分為：

Zn×k=Xn×pβp×k.

(4)

第二種方法是對Xn×p進行奇異值(Singular Value Decomposition, SVD)分解：

X=UDβT，

(5)

Zn×k=UD為X通過主成分分析降維后得到的前k個主成分矩陣,βp×k為前k個主成分對應的負載矩陣，Dp×p為對角陣，對角線元素為降序排列的原始矩陣的奇異值。

上述兩種方法適用不同類型數據矩陣。若以數據矩陣的每一行代表一個樣本，每一列代表一種屬性，第一種方法適用于樣本數量遠大于屬性數量的數據矩陣，第二種方法適用于樣本數量遠小于屬性數量的矩陣。

各個主成分對原數據矩陣的解釋程度不同，一般通過主成分的方差貢獻率(也稱方差解釋率)來衡量其在原始數據矩陣中所占的比重。對方法一、二，主成分的方差貢獻率分別使用式(6)、(7)計算，前k個主成分的累計方差貢獻率分別使用式(8)、(9)計算。

(6)

式中，λi為原始矩陣Xn×p的協(xié)方差矩陣的第i個特征值(按降序排列)。

(7)

(8)

(9)

2.2 稀疏主成分分析

稀疏主成分分析是在主成分分析的基礎上，將主成分分析所得負載矩陣作為系數矩陣，與原始數據矩陣通過最小二乘法逼近主成分矩陣，同時引入L1正則化和L2正則化作為約束[14]。如式(10)所示:

(10)

3 圖像噪聲估計方法

以圖1所示的8幅尺寸均為512 pixel×512 pixel的灰度圖進行線性關系測試。對8幅圖分別添加標準差δ=5，10，20，30，…，70的高斯白噪聲，使得每一幅圖均對應生成8幅新的高斯白噪聲圖像。對噪聲標準差等級δ使用多個稀疏負載向量均值的加權平均進行線性擬合，如式(11)所示。

(11)

式中，δi代表第i幅噪聲圖像所對應的高斯白噪聲標準差等級。使用式(11)反復迭代計算之后確定r1=1,r2=5,ωj=1(j=1,2,…,5)時擬合效果最好。此時式(11)可重寫為式(12)。

(12)

此時各圖的擬合結果如圖2所示。

在實際應用中，對單個圖像而言，因未知量過多，因此式(12)無法求出唯一收斂解。為解決此問題，本文采用構造超定方程組的方式求解單個圖像的未知量。方法如下所述：

設原始圖像的像素矩陣為X0，添加了未知標準差等級δ0的高斯白噪聲n0之后得到的新圖像為X0+n0=Y0。對X0和Y0有：

D(Y0)=D(X0+n0)=D(X0)+δ02.

(13)

然后分別對Y0添加已知標準差等級δi(i=1,2,…,N)的高斯白噪聲ni(i=1,2,…,N)，分別得到Yi(i=1,2,…,N)。對Yi(i=1,2,…,N)有：

D(Yi)=D(X0+n0+ni)=D(X0)+δ02+δ2i，

(14)

則結合式(12)與(14)可得到：

(15)

(16)

圖1 用于線性關系測試的圖像Fig.1 Images for linear relationship testing

圖2 參數最優(yōu)時的擬合結果Fig.2 Fitting results of the optimal parameters

N)之后生成的Yi(i=1,2,…,N)進行稀疏主成分分析所得到的第j個稀疏負載向量均值。

4 實驗結果和分析

4.1 灰度圖像的噪聲估計

Lena標準灰度圖像被用于進行噪聲估計測試。文中所提出的噪聲估計方法與文獻[11,15]中提出的方法進行了比較。設δ0為圖像高斯白噪聲的真實標準差等級，δ1為使用估計方法得出的高斯白噪聲標準差等級的最優(yōu)估計。實驗結果如圖3和表1所示。

圖3 不同噪聲等級下3種噪聲估計方法的噪聲估計誤差Fig.3 Estimation error of three noise estimation methods under different noise levels

從圖4和表1可以看出，相比于其他兩種算法，本方法的精確度更高且更具有魯棒性。

4.2 彩色圖像的噪聲估計

彩色圖像包含R、G、B 3個像素矩陣。在實驗中，分別對R、G、B 3個像素矩陣采用所提出的噪聲估計方法進行噪聲估計。實驗所使用的Lena圖像的3個通道如圖4所示。與文獻[11,15]的比較結果如圖5和表2所示。

表1 不同噪聲等級下3種噪聲估計方法的估計結果

Tab.1 Estimation results of three noise estimation methods under different noise levels

δ0文獻[11]文獻[15]本文方法54.124.324.78109.609.2410.112021.0221.3120.573032.4231.5530.734041.5841.0840.365051.3652.7750.396063.8364.1961.027072.6671.8670.62

從圖5和表2可以看出，本文方法在3個通道都達到了最高的精確度，且估計效果不受噪聲等級影響。

(a)R通道(a) R channel

(b)G通道(b) G channel

(c)B通道(c) B channel圖4 Lena圖像的R、G、B 3個通道Fig.4 R, G, B channels of Lena image

(a)R通道(a) R channel

(b) G channel(b)G通道

(c)B通道 (c) B channel

5 結 論

高斯噪聲圖像的噪聲標準差等級與稀疏主成分分析得到的部分主成分負載向量均值呈現較為良好的線性關系，基于此特征，本文中提出了一種基于稀疏主成分分析的圖像高斯噪聲估計方法。在該方法中，經過實驗確定了前5個稀疏負載向量均值之和與噪聲等級的線性度最優(yōu);同時通過在圖像中依次添加已知噪聲等級的高斯白噪聲的方式構造超定方程組，在反復迭代之后求解出噪聲等級的最優(yōu)估計量。在標準圖像Lena的灰度圖和彩色圖像三個通道的測試中,該方法展現出了相比于其他兩種算法更高的精確度，在低噪聲(δ0=5)到高噪聲(δ0=70)條件下均具有較高的估計精度和較強的魯棒性，具有一定的工程實用價值。