孟 玉

(西北工業(yè)大學計算機學院，陜西 西安 710129)

0 引言

高光譜圖像處理技術(shù)被廣泛用于地質(zhì)勘探、氣象監(jiān)測、環(huán)境保護等領(lǐng)域，然而人們在獲取HS數(shù)據(jù)的同時，卻不可避免地引入了噪聲，如何有效地估計及去除噪聲［1］對HS圖像處理尤為重要。

空間光譜維去相關(guān)法［2］(Spectral and Spatial Decorrelation Method，SSDC)是一種經(jīng)典的評噪方法，雖然該方法較穩(wěn)定，但因強硬地將圖像分為等大的小塊，無法應(yīng)對具有復(fù)雜紋理的HS圖像，且評估結(jié)果誤差較大。局部均勻塊標準差法［3］(Local Homogeneous Block Standard Deviations Method，LHBSD)和同質(zhì)區(qū)域劃分光譜維去相關(guān)法［4-5］(Homogeneous Regions Division and Spectral Decorrelation Method，HRDSDC)則是基于同質(zhì)區(qū)域進行評噪，對復(fù)雜紋理的HS圖像依然適用，但HR的劃分結(jié)果并不穩(wěn)定，存在隨機性，影響了噪聲估計結(jié)果。Martin-Herrero［6］對HR劃分方法進行了評估和改進，但改進后的算法仍要高深度遞歸，計算負擔過重。

傳統(tǒng)方法［7］只適用于估計信號無關(guān)(Signal Independent，SI)噪聲，而電子器件敏感性的提升為HS圖像引入了光電噪聲，也即信號相關(guān)［8］(Signal Depend-ent，SD)噪聲，導(dǎo)致傳統(tǒng)方法難于進行估計。Qian［9-10］等人提出的三維非局部均值濾波方法評噪和去噪的效果顯著，而Acito［11］等人也提出了一種自動化的評噪方法，通過矩和基于梯度的多變量迭代方法求解噪聲模型的最大似然函數(shù)，繼而計算SD和SI的方差。這些方法具有較高的穩(wěn)定性和精確性，但復(fù)雜度較高，計算負擔較大。

為了解決傳統(tǒng)算法的混合噪聲估計問題，同時為了提高估計算法的性能，提出一種基于HR分割的混合噪聲估計算法。算法的主要思想是在HR分割區(qū)域內(nèi)進行MLR估計，從而得到含較強空間和光譜相關(guān)性的純信號估計，而該信號與原圖像間的殘差即為混合噪聲。得到混合噪聲后，可以利用信噪比和比例因子計算得到其內(nèi)部詳細的SD和SI的方差。該方法避免了復(fù)雜的預(yù)計算和后處理，不僅降低了算法復(fù)雜度，且能準確估計混合噪聲及其內(nèi)部參數(shù)。同時為了得到較好的HR分割并提升算法性能，提出一種最近鄰劃分的無監(jiān)督HR分割算法，該算法結(jié)合HS的空間和光譜特性對圖像進行掃描并分割，無需遞歸調(diào)用，算法高效穩(wěn)定且結(jié)果要優(yōu)于其它傳統(tǒng)方法。

1 同性質(zhì)區(qū)域分割

同性質(zhì)區(qū)域［5］定義為在空間上相鄰，且光譜特性相同或相似的大片連通區(qū)域，因此其內(nèi)部的像元具有較強的空間和光譜相關(guān)性，若除去HR內(nèi)像元間的相關(guān)性后仍存在誤差，則該誤差即由噪聲引起。通過在HR內(nèi)使用MLR去除強相關(guān)性，即可得到誤差最小的殘差，進而對殘差進行估計從而得到噪聲。理論上來看，HR的分割結(jié)果對噪聲的估計有較大影響。

此處使用光譜角距離(Spectral Angle Distance，SAD)［12-13］來度量像元間的光譜相似性。若將單個像元的光譜特性描述為一個n維向量，n是像元光譜的波段數(shù)，則每一個向量具有確定的長度和方向，分別代表了像元的亮度和光譜特征，由此得到2個像元間的SAD為:

其中t和r分別代表目標像元向量和參考像元向量。光譜角θ的單位為弧度，取值范圍為0～π/2，兩個像元間的光譜角越小，它們的光譜特征就越相似。由式(1)很容易看出，當t和r有尺度變化時，光譜角θ依然不變，這種尺度不變的特性可以很好地滿足刻畫高光譜圖像像元相似性的要求。

由于HR內(nèi)的像元具有相似的光譜特征，故可將HR內(nèi)所有像元向量的加權(quán)平均作為該區(qū)域的光譜特征。設(shè)參考區(qū)域R={r1，r2，…，rN}為包含N個像元的同質(zhì)區(qū)域，目標像元t與區(qū)域R的SAD為:

同性質(zhì)區(qū)域劃分算法步驟如下:

(1)創(chuàng)建一個與原始圖像同樣尺寸的劃分模板，并將左上角的像元標記為1，作為第一個區(qū)域。按照從左到右，從上到下的順序?qū)ふ椅礃擞浀南裨?，轉(zhuǎn)步驟(2)。

(2)尋找與目標像元上相鄰和左相鄰的參考像元(至少可以找到一個)。若只找到一個相鄰的參考像元，且找到的參考像元是單個像元，轉(zhuǎn)步驟(3);若參考像元屬于某個同質(zhì)區(qū)域，轉(zhuǎn)步驟(4)。若找到兩個相鄰的參考像元，轉(zhuǎn)步驟(5)。

(3)利用式(1)計算SAD，若SAD小于某一閾值μ則將目標像元與參考像元歸為同一同質(zhì)區(qū)域，并在劃分模板上目標像元的位置標記上參考像元的標號;否則將目標像元劃分為新的區(qū)域，劃分模板上標記為當前最大標號加1，轉(zhuǎn)步驟(6)。

(4)利用式(2)計算SAD，若SAD小于某一閾值μ則將目標像元歸入?yún)⒖嫉耐|(zhì)區(qū)域，并在劃分模板上目標像元位置處標記上參考區(qū)域的標號;否則將目標像元劃分為新的區(qū)域，劃分模板上標記為當前最大標號加1，轉(zhuǎn)步驟(6)。

(5)像元間或像元與區(qū)域間的SAD計算及處理參考步驟(3)和步驟(4)。當計算得到的2個SAD值都小于閾值μ時，將參與計算的兩個參考像元或區(qū)域和目標像元合并為一個區(qū)域，并更改標記。否則將目標像元劃分為新的區(qū)域，劃分模板上標記為當前最大標號加1，轉(zhuǎn)步驟(6)。

(6)繼續(xù)轉(zhuǎn)步驟(2)進行劃分，尋找未標記像元，直到所有像元劃分完畢，算法結(jié)束。

上述算法中的閾值μ是SAD的分割門限，需手動設(shè)置，關(guān)于它對算法的影響和討論將在后續(xù)的實驗部分給出。

2 混合噪聲估計算法

此處假設(shè)HS圖像僅受混合噪聲的影響，且混合噪聲是作為零均值的高斯隨機過程存在于圖像中的。一種普遍使用的混合噪聲模型［11］可表示為:

其中σ2(n)為混合噪聲的方差，μ(f)為圖像均值，σ2(u)和σ2(m)分別為SD和SI噪聲的方差。此處的目的就是估計得到除μ(f)之外的其它所有參數(shù)。

多元線性回歸(MLR)［11］是一種普遍使用的估計方法，它利用了HS圖像的高空間/光譜維相關(guān)性。假設(shè)HS圖像共L波段，每波段共有M個像元，若將第 l波段的圖像整理成行向量 Xl= ［，...，］，其中為第l波段第m個像元的DN值。再將某個像元m在其它L-1個波段上的DN值寫成列向量則某波段所有像元的向量可以組成矩陣Fl=由此得到第l波段的可用信號估計:

其中Wl為第l波段的線性混合系數(shù)。

對HS圖像的每個波段使用MLR進行估計，即可得到可用信號估計，其均值與圖像均值μ(f)相等。再對殘差rl=Xl－進行無偏估計，即可得到混合噪聲的標準差σ(n)。

為了對SD和SI的標準差進行估計，可在給定圖像混合噪聲信噪比SNRl=θ的情況下引入SD和SI的比例因子α，有:

其中第l波段SD和SI噪聲的信噪比為［11］:

結(jié)合式(3)和式(5)、式(6)即可解出:

綜上所述，給出混合噪聲的估計算法:

(1)首先對HS圖像進行分割得到HR，為了消除小區(qū)域和雜點的影響，實驗中可只保留像元個數(shù)大于某閾值的HR參與計算。

(2)對每個波段內(nèi)的HR分別使用MLR進行估計，并將所有計算結(jié)果的統(tǒng)計平均值作為本波段的最終估計，即可得到每個波段的混合噪聲參數(shù)σ2(n)和圖像均值μ(f)。

(3)給出混合噪聲的信噪比θ和比例因子α，通過式(7)計算得到每個波段的SD和SI噪聲的標準差 σSD和 σSI，算法結(jié)束。

3 實驗結(jié)果與分析

3.1 在仿真數(shù)據(jù)上實驗

從美國地質(zhì)勘探局(United States Geological Survey)的光譜庫中獲取5種物質(zhì)，其光譜特征如圖1所示。實驗選取其中350 nm到2500 nm(間隔10 nm)共216個波段，將其兩兩混合(從左到右，從上到下依次為:明礬石+方解石、黑松+方解石、柳樹+方解石、混凝土+方解石)，并依據(jù)HS線性混合模型產(chǎn)生200×200×216的HS仿真數(shù)據(jù)，如圖2所示。區(qū)域劃分模板如圖3所示。實驗在酷睿(R)E6500處理器(2.94 GHz)、2 GB 內(nèi)存、32位 Windows 7操作系統(tǒng)、Matlab環(huán)境下進行。

圖1 光譜特征圖

圖2 實驗數(shù)據(jù)第60波段(940 nm)的仿真圖像

圖3 區(qū)域劃分模板

圖4 α=3時混合噪聲(θ=30 dB)估計對比

圖5 α取不同值(θ=30 dB)SD和SI噪聲估計對比(實線為實際噪聲，點為估計噪聲)

表1 混合噪聲標準差估計RMSE(θ=30 dB)

表2 SD和SI噪聲標準差估計RMSE(θ=30 dB)

實驗加入θ為30 dB，α分別為1/3、1和3的3組混合噪聲到仿真數(shù)據(jù)中，然后將 SSDC［2］、LHBSD［3］、HRDSDC［5］和本文的方法進行對比，其中一組估計結(jié)果如圖4所示。此處使用相對均方誤差［11］(Relative Mean Square Error，RMSE)對估計方法進行定量評估:

圖5和表2給出本文方法在3組數(shù)據(jù)間對SD和SI噪聲的估計結(jié)果和RMSE對比，可以看出估計結(jié)果較優(yōu)，且α僅影響到SD和SI在混合噪聲中的比例關(guān)系，對估計結(jié)果的精度影響不大。

3.2 在AVIRIS數(shù)據(jù)上實驗

選取美國加利福尼亞州薩利納斯山谷的AVIRIS數(shù)據(jù)進行實驗，其圖像分辨率為512×217，空間分辨率為3.7 m，包含224個波段(已去除水扭曲波段)，其第100波段灰度圖如圖6(a)所示。

圖6 (a)AVIRIS數(shù)據(jù)第100波段灰度圖

實驗對比了 HRDSDC［5］、Martin-Herrero［6］和本文方法的HR分割結(jié)果，并給出本文方法在不同閾值時的HR分割結(jié)果。對比圖6(c)、(d)、(e)在(b)中的黑框部分可以看出，閾值相同時本文方法的分割結(jié)果較優(yōu)。對比圖6(e)、(f)、(g)可看出，在閾值過小時易出現(xiàn)過分割現(xiàn)象，閾值過大時則易造成錯分少分的現(xiàn)象，可見閾值的選擇對HR分割有較大影響。經(jīng)過多次實驗，最終選擇閾值為0.05時可得到較滿意的HR分割。

圖7 閾值不同時多種方法間混合噪聲估計結(jié)果對比

圖7顯示了取不同閾值時各種方法間的混合噪聲估計結(jié)果對比。可以看出，HRDSDC受閾值影響較大，其改進方法受影響較小，而本文方法在不同閾值下都能獲得較一致的估計結(jié)果，可見算法相當穩(wěn)定，對HR分割結(jié)果的依賴性較低。

4 結(jié)束語

本文提出了一種基于同性質(zhì)區(qū)域分割的高光譜圖像混合噪聲估計方法。在HR分割方面，解決了HRDSDC具有的隨機性問題，同時在算法上避免了大量的遞歸調(diào)用。在評噪方面，解決了傳統(tǒng)評噪方法不能估計SD噪聲的問題，同時降低了算法復(fù)雜度。實驗表明，本文提出的方法在穩(wěn)定性和精確性上都要優(yōu)于一般傳統(tǒng)方法，具有一定的工程應(yīng)用價值。

［1］張兵，高連如.高光譜圖像分類與目標探測［M］.北京:科學出版社，2011:25-46.

［2］Roger R E，Arnold J F.Reliably estimating the noise in AVIRIS hyperspectral images［J］.International Journal of Remote Sensing，1996，17(10):1951-1962.

［3］高連如，張兵，張霞，等.基于局部標準差的遙感圖像噪聲評估方法研究［J］.遙感學報，2007，11(2):201-208.

［4］Martin-Herrero J.Hybrid object labelling in digital images［J］.Machine Vision and Applications，2007，18(1):1-15.

［5］Gao L R，Zhang B，Zhang X，et al.A new operational method for estimating noise in hyperspectral images［J］.IEEE Geoscience and Remote Sensing Letters，2008，5(1):83-87.

［6］Martin-Herrero J.Comments on“a new operational method for estimating noise in hyperspectral images”［J］.IEEE Geoscience and Remote Sensing Letters，2008，5(4):705-709.

［7］Rakwatin P，Takeuchi W，Yasuoka Y.Stripe noise reduction in MODIS data by combining histogram matching with facet filter［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45(6):1844-1856.

［8］Uss M L，Vozel B，Lukin V V，et al.Local signal-dependent noise variance estimation from hyperspectral textural images［J］.IEEE Journal of Selected Topics in Signal Processing，2011，5(3):469-486.

［9］Qian Yuntao，Shen Yanhao，Ye Minchao，et al.3-D nonlocal means filter with noise estimation for hyperspectral imagery denoising［C］//Proceedings of the 2012 IEEE International Geoscience and Remote Sensing Symposium.2012:1345-1348.

［10］Qian Yuntao，Ye Minchao.Hyperspectral imagery restoration using nonlocal spectral-spatial structured sparse representation with noise estimation［J］.IEEE Journal of Select-ed Topics in Applied Earth Observations and Remote Sensing，2013，6(2):499-515.

［11］Acito N，Diani M，Corsini G.Signal-dependent noise modeling and model parameter estimation in hyperspectral images［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49(8):2957-2971.

［12］何中海，何彬彬.基于權(quán)重光譜角制圖的高光譜礦物填圖方法［J］.光譜學與光譜分析，2011，31(8):2200-2204.

［13］Andreou C，Karathanassi V.A novel multiple endmember spectral mixture analysis using spectral angle distance［C］//Proceedings of the 2012 IEEE International Geoscience and Remote Sensing Symposium.2012:4110-4113.