基于同質(zhì)區(qū)域分割的高光譜圖像混合噪聲估計(jì)

孟 玉

(西北工業(yè)大學(xué)計(jì)算機(jī)學(xué)院，陜西 西安 710129)

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

(3)利用式(1)計(jì)算SAD，若SAD小于某一閾值μ則將目標(biāo)像元與參考像元?dú)w為同一同質(zhì)區(qū)域，并在劃分模板上目標(biāo)像元的位置標(biāo)記上參考像元的標(biāo)號(hào);否則將目標(biāo)像元?jiǎng)澐譃樾碌膮^(qū)域，劃分模板上標(biāo)記為當(dāng)前最大標(biāo)號(hào)加1，轉(zhuǎn)步驟(6)。

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

(5)像元間或像元與區(qū)域間的SAD計(jì)算及處理參考步驟(3)和步驟(4)。當(dāng)計(jì)算得到的2個(gè)SAD值都小于閾值μ時(shí)，將參與計(jì)算的兩個(gè)參考像元或區(qū)域和目標(biāo)像元合并為一個(gè)區(qū)域，并更改標(biāo)記。否則將目標(biāo)像元?jiǎng)澐譃樾碌膮^(qū)域，劃分模板上標(biāo)記為當(dāng)前最大標(biāo)號(hào)加1，轉(zhuǎn)步驟(6)。

(6)繼續(xù)轉(zhuǎn)步驟(2)進(jìn)行劃分，尋找未標(biāo)記像元，直到所有像元?jiǎng)澐滞戤?，算法結(jié)束。

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

2 混合噪聲估計(jì)算法

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

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

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

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

對(duì)HS圖像的每個(gè)波段使用MLR進(jìn)行估計(jì)，即可得到可用信號(hào)估計(jì)，其均值與圖像均值μ(f)相等。再對(duì)殘差rl=Xl－進(jìn)行無(wú)偏估計(jì)，即可得到混合噪聲的標(biāo)準(zhǔn)差σ(n)。

為了對(duì)SD和SI的標(biāo)準(zhǔn)差進(jìn)行估計(jì)，可在給定圖像混合噪聲信噪比SNRl=θ的情況下引入SD和SI的比例因子α，有:

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

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

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

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

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

(3)給出混合噪聲的信噪比θ和比例因子α，通過(guò)式(7)計(jì)算得到每個(gè)波段的SD和SI噪聲的標(biāo)準(zhǔn)差 σSD和 σSI，算法結(jié)束。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 在仿真數(shù)據(jù)上實(shí)驗(yàn)

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

圖1 光譜特征圖

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

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

圖4 α=3時(shí)混合噪聲(θ=30 dB)估計(jì)對(duì)比

圖5 α取不同值(θ=30 dB)SD和SI噪聲估計(jì)對(duì)比(實(shí)線為實(shí)際噪聲，點(diǎn)為估計(jì)噪聲)

表1 混合噪聲標(biāo)準(zhǔn)差估計(jì)RMSE(θ=30 dB)

表2 SD和SI噪聲標(biāo)準(zhǔn)差估計(jì)RMSE(θ=30 dB)

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

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

3.2 在AVIRIS數(shù)據(jù)上實(shí)驗(yàn)

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

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

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

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

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

4 結(jié)束語(yǔ)

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

［1］張兵，高連如.高光譜圖像分類(lèi)與目標(biāo)探測(cè)［M］.北京:科學(xué)出版社，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］高連如，張兵，張霞，等.基于局部標(biāo)準(zhǔn)差的遙感圖像噪聲評(píng)估方法研究［J］.遙感學(xué)報(bào)，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］.光譜學(xué)與光譜分析，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.