基于高光譜影像的SG濾波算法的研究

何英杰 謝東海 鐘若飛

(首都師范大學(xué)資源環(huán)境與旅游學(xué)院， 北京 100048)

0 引 言

經(jīng)歷全色到多光譜掃描儀成像階段之后，在20世紀(jì)80年代初由于成像光譜概念的出現(xiàn)，遙感技術(shù)進(jìn)入了高光譜遙感的嶄新的階段.高光譜遙感是當(dāng)前遙感的前沿技術(shù)之一，是將成像技術(shù)和光譜技術(shù)相結(jié)合的多維信息獲取的技術(shù).它是由多光譜遙感成像技術(shù)的基礎(chǔ)發(fā)展而來的，光譜分辨率得到了大幅的提高[1]，為每個(gè)像元提供幾乎連續(xù)的地物光譜曲線,使我們利用高光譜反演陸地細(xì)節(jié)成為可能[2].

獲取較高的光譜分辨率的連續(xù)的圖像數(shù)據(jù)是高光譜遙感數(shù)據(jù)的閃耀之處,它能充分地反映地物反射率光譜的細(xì)微特征.然而，高光譜圖像成像儀器以及環(huán)境等因素對(duì)數(shù)據(jù)獲取造成的影響，使得在高光譜數(shù)據(jù)成像和預(yù)處理的各個(gè)環(huán)節(jié)中，都可能引入不同程度的噪聲，高光譜圖像具有空間圖像和地物光譜兩個(gè)方面的信息，因而噪聲對(duì)高光譜圖像的影響最終也會(huì)表現(xiàn)為空間域噪聲和光譜域噪聲兩方面[3]，嚴(yán)重地限制了高光譜數(shù)據(jù)的高光譜分辨率的優(yōu)勢.

人們根據(jù)實(shí)際圖像的特點(diǎn)、噪聲的統(tǒng)計(jì)特征和頻譜分析的規(guī)律，發(fā)展了各式各樣的去噪方法.其中最直觀的方法是根據(jù)噪聲能量一般集中于高頻，而圖像頻譜分布于一個(gè)有限區(qū)間的這一特點(diǎn)，采用低通濾波方法.在基于三角網(wǎng)格模型的降噪算法研究方面，Yu等[4]提出一種兩階段特征保持網(wǎng)格去噪框架.此外，對(duì)圖像進(jìn)行平滑處理也是常用的方法.平滑的目的有兩個(gè)：改善圖像質(zhì)量和抽出對(duì)象特征.一個(gè)較好的平滑方法應(yīng)該既可以消掉噪聲影響又不使圖像的邊緣輪廓和線條變模糊.劉倩[5]基于非局部稀疏的圖像處理方法，充分利用自然圖像本身的有用信息，有效彌補(bǔ)了數(shù)學(xué)模型的缺陷，以數(shù)據(jù)驅(qū)動(dòng)的方式大大提高了圖像平滑的效果.圖像平滑處理方法有空間域處理和頻率域處理兩大類.在空間域里一類方法是噪聲去除，即先判定某點(diǎn)是否為噪聲點(diǎn)，若是，重新賦值，若不是按原值輸出；另一類是平均，即不依賴于噪聲點(diǎn)的識(shí)別和去除，而對(duì)整個(gè)圖像進(jìn)行平均運(yùn)算.而在頻率域里是對(duì)圖像頻譜進(jìn)行修正，一般采用低通濾波方法，而不像在空間域里直接對(duì)圖像像素灰度級(jí)值進(jìn)行運(yùn)算.圖像的均值濾波模板和中值濾波模板是實(shí)現(xiàn)圖像的平滑去噪的主要方法.

由于高光譜圖像成像時(shí)自然光照明條件影響、地面地形影響、混合像元問題等原因引入各種噪聲，這些噪聲存在于光譜域中，同樣也是需要檢測并予以去除的[6].MNF(minimum noise fraction)變換通過降低數(shù)據(jù)的維數(shù)來達(dá)到隔離數(shù)據(jù)噪聲的目的，首先使用噪聲的協(xié)方差矩陣分解并將數(shù)據(jù)中的噪聲白化[7],此時(shí)噪聲具有單位方差，并且與原波段無相關(guān)性.然后針對(duì)上述噪聲白化的數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)的主成分變換；Boardman[8]針對(duì)這個(gè)問題提出了EFFORT方法，即通過統(tǒng)計(jì)高光譜圖像，尋找一種微調(diào)來減小誤差，從而提高反射率的精度；Gladden[9]同樣認(rèn)為反射率轉(zhuǎn)換后光譜中會(huì)產(chǎn)生很大的噪聲，并采用一種 SSE線性校正方法，即將轉(zhuǎn)換后反射率光譜擬合到實(shí)測地物光譜上，達(dá)到去除噪聲的目的.陳志剛和束炯[10]基于光譜二階導(dǎo)數(shù)檢測噪聲的基礎(chǔ)上，提出了基于經(jīng)驗(yàn)?zāi)B(tài)分解EMDF(empirical mode decomposition based filter)的光譜域?yàn)V波方法及其修正，具有較好的光譜域?yàn)V波效果.由于多分辨率的分析特性和較好的時(shí)頻局部化，小波分析[11]依據(jù)小波基函數(shù)，多層次地分解高光譜圖像光譜.然后針對(duì)分解后的小波系數(shù)選取合適的閾值進(jìn)行處理，最終獲得濾除噪聲后的光譜曲線，在圖像處理領(lǐng)域得到廣泛應(yīng)用.楊浩等[12]針對(duì)成像光譜儀采集的高光譜圖像數(shù)據(jù)，提出一種小波閾值去噪方法去降低噪聲影響.遙感影像存在的光譜鋸齒噪聲[13]，同樣也嚴(yán)重影響礦物光譜吸收參數(shù)的提取.

Savitzky和Golay[14]提出的最小二乘濾波算法，濾波器的選擇為Savitzky-Golay(SG) 濾波器，它是一種特殊的低通濾波器，最初由Savitzky和Golay于1964年提出.光譜曲線由于受各種噪聲的影響而呈現(xiàn)鋸齒狀，經(jīng)SG濾波處理后，平滑一些細(xì)微的鋸齒噪聲同時(shí)能保持整條曲線大的光譜特征不受影響[15-16].本文主要以航空高光譜遙感數(shù)據(jù)為研究對(duì)象，在SG濾波器的基礎(chǔ)之上進(jìn)行優(yōu)化，求取反射率光譜的二階導(dǎo)數(shù)，以期使其能更加靈活且準(zhǔn)確地處理影像.

1 高光譜影像SG濾波法分析

在平面坐標(biāo)系中,用一條曲線來擬合一組數(shù)據(jù),不妨假設(shè)這條曲線為y=a0+a1x+a2x2+a3x3+a4x4，當(dāng)每一個(gè)點(diǎn)的橫坐標(biāo)代入這個(gè)曲線方程后，所得值與該點(diǎn)的縱坐標(biāo)之差的平方之和最小時(shí)，這條曲線的擬合度最高，從而可以確定所有的系數(shù)ai(i=0,1,2,3,4).

從高光譜影像中每個(gè)像元提取出的反射率曲線本身有一定的物理意義，可以反映地物的細(xì)微特征，但是轉(zhuǎn)換后的像素光譜反射率曲線的噪聲一般呈鋸齒狀分布，這些鋸齒分布即為高光譜圖像光譜域噪聲，嚴(yán)重地限制了其高光譜分辨率的優(yōu)勢.為了去除光譜域噪聲，考慮采用多項(xiàng)式擬合的最小二乘平滑算法，剔除一些細(xì)微的鋸齒噪聲同時(shí)能保持整條曲線大的光譜特征不受影響.SG濾波法[14]是一種經(jīng)典的最小二乘平滑算法，它使用簡化的最小二乘擬合卷積方法對(duì)曲線進(jìn)行平滑處理并可計(jì)算平滑后曲線的各階導(dǎo)數(shù).

擴(kuò)展一下，通常假設(shè)曲線為p次多項(xiàng)式，即滿足公式(1)

Yi=a0+a1i+a2i2+…+apip

(1)

其中，Yi第i點(diǎn)平滑后的數(shù)值.

用上述多項(xiàng)式擬合曲線的誤差可以表示為：見公式(2)

(2)

其中，yj代表平滑前的數(shù)值，平滑窗口大小k=2m+1.為了讓誤差S最小，需要對(duì)S進(jìn)行偏微分，使得所求S偏微分各項(xiàng)為零.因此，可以得到以下各式：關(guān)系如公式(3)

…

(3)

Savitzky對(duì)整個(gè)求解過程進(jìn)行了推導(dǎo)，給出了平滑窗中心點(diǎn)平滑后數(shù)值的最后公式及公式中系數(shù)的計(jì)算方法[14].Madden[17]基于Savitzky的計(jì)算公式，改正了原SG系數(shù)的一些錯(cuò)誤，并且給出了各階導(dǎo)數(shù)平滑系數(shù)修正后的計(jì)算公式.根據(jù)其推導(dǎo)平滑后的曲線數(shù)值的最小二乘卷積方程可寫為：見公式(4)

(4)

(5)

式中，m為平滑窗寬度的一半[18].

SG濾波器有如下幾大優(yōu)點(diǎn)[19]：

1)利用最小二乘的多項(xiàng)式擬合方法非常清晰易懂，并且在計(jì)算上來說，多項(xiàng)式卷積的操作比最小二乘的計(jì)算可操作性更強(qiáng)；

2)濾波系數(shù)只需要在對(duì)應(yīng)的卷積系數(shù)表中進(jìn)行查找，很容易獲得；

3)SG濾波器可以有任意的長度，因此有利于采樣頻率通常很低的生物學(xué)或者生物力學(xué)的數(shù)據(jù)處理.

SG算法的本質(zhì)是一種最小二乘法的卷積平滑，根據(jù)多項(xiàng)式擬合階數(shù)和平滑次數(shù)進(jìn)行綜合考量，但是多項(xiàng)式擬合曲線的誤差求取的是一階導(dǎo)數(shù).本文高光譜影像噪聲檢測SG算法是基于光譜二階導(dǎo)數(shù).光譜導(dǎo)數(shù)技術(shù)對(duì)光譜信噪比非常敏感，光譜的一階、二階導(dǎo)數(shù)計(jì)算公式如下[20-21]：

(6)

(7)

其中，近似選取相鄰兩點(diǎn)，其關(guān)系滿足公式(8)

Δλ=λk-λj=λj-λi,λk>λj>λi

(8)

光譜的低階導(dǎo)數(shù)處理對(duì)噪聲影響敏感性較低，而高階微分對(duì)噪聲影響敏感度高.相比較一階導(dǎo)數(shù)，二階導(dǎo)數(shù)很好的反映了噪聲的實(shí)際分布情況(高階導(dǎo)數(shù)過于復(fù)雜).因此，可以考慮用光譜二階導(dǎo)數(shù)對(duì)光譜曲線進(jìn)行噪聲影響程度檢測.

對(duì)于光譜曲線進(jìn)行光譜二階導(dǎo)數(shù)計(jì)算，考慮將其與平均值進(jìn)行相減，然后一般想法是參考差值分布，與一個(gè)固定的值進(jìn)行比較，即設(shè)定一個(gè)閾值，判定噪聲與否.然后根據(jù)SG濾波法通過對(duì)反射率光譜二階導(dǎo)數(shù)的統(tǒng)計(jì)，得到光譜各波段中噪聲的判定函數(shù)，在噪聲較嚴(yán)重的波段進(jìn)行較多點(diǎn)數(shù)的SG濾波，對(duì)噪聲不太嚴(yán)重的波段進(jìn)行較少點(diǎn)數(shù)的SG濾波，從而確定該波段濾波平滑窗的大小.根據(jù)計(jì)算表明，如果判別式滿足公式(9):

(9)

那么該波段反射率的值認(rèn)為是有較大噪聲的存在，其中σ表示光譜二階導(dǎo)數(shù)的標(biāo)準(zhǔn)差.

在噪聲檢測完成后，可以先采取較大窗口的SG濾波對(duì)較大噪聲波段進(jìn)行局部第一次平滑，其他波段保持不變，得到初步的濾波結(jié)果；然后采取較小窗口的SG濾波對(duì)初次濾波的結(jié)果進(jìn)行光譜曲線的整體平滑.這樣的好處是，既可以對(duì)光譜曲線進(jìn)行平滑濾波和噪聲去除，又可以最大程度的對(duì)光譜細(xì)節(jié)進(jìn)行保留.

本算法的具體步驟如下：

1)從高光譜圖像中逐像元的提取光譜特征，繪制反射率曲線；

2)計(jì)算曲線的二階導(dǎo)數(shù)；

3)計(jì)算光譜二階導(dǎo)數(shù)的均值與方差，按照公式(9)判定噪聲；

4)進(jìn)行第一次平滑，對(duì)判定為噪聲的波段進(jìn)行例如11點(diǎn)SG濾波平滑，保持其他波段不變；

5)進(jìn)行第二次平滑，對(duì)第一次平滑后的光譜曲線進(jìn)行例如5點(diǎn)(一般選擇小窗口的值為大窗口值的一半取整)的SG平滑;

6)該點(diǎn)處理完畢，轉(zhuǎn)入下一像素點(diǎn)的處理，轉(zhuǎn)入步驟2，直到所有點(diǎn)處理完畢.

2 實(shí)驗(yàn)驗(yàn)證與分析

2.1 影像數(shù)據(jù)選擇

試驗(yàn)數(shù)據(jù)為2012年武漢市某地區(qū)PHI高光譜遙感影像，包含80個(gè)波段，光譜范圍400～850 nm，大小為128行、128列.選取比較有代表性的波段進(jìn)行展示，例如49、51、53波段，如圖1所示.

圖1 49、51、53波段影像

圖像質(zhì)量評(píng)價(jià)的研究是圖像信息學(xué)科的基礎(chǔ)研究之一.對(duì)圖像處理或圖像通信系統(tǒng)，其信息的主體是圖像，衡量這個(gè)系統(tǒng)的重要指標(biāo)，就是圖像的質(zhì)量，而圖像去噪就是為了改變圖像的主觀視覺顯示圖像質(zhì)量[22].我們需要考慮的幾個(gè)因素總結(jié)如下：1)去噪后圖像應(yīng)盡量的平滑，不存在或有較少的噪聲痕跡.2)去噪結(jié)果不能使圖像過渡的失去結(jié)構(gòu)細(xì)節(jié)而變得模糊.3)沒有由于具體去噪方法產(chǎn)生的人工噪聲.4)均方差盡可能小.所有這些都要求有一個(gè)合理的圖像質(zhì)量評(píng)價(jià)方法.根據(jù)圖像質(zhì)量評(píng)價(jià)方法的思想，初步對(duì)所給數(shù)據(jù)進(jìn)行分析：波段存在劃痕狀條帶噪聲，以及分布不均勻的隨機(jī)噪聲，以49波段噪聲特點(diǎn)較為明顯.根據(jù)光譜域處理噪聲的結(jié)果與一般平滑去噪的不同，要盡量保證噪聲點(diǎn)附近平滑處理，而不影響非噪聲點(diǎn)光譜特性.基于高光譜影像應(yīng)用SG濾波算法，可以滿足預(yù)期設(shè)想.

2.2 結(jié)果分析

根據(jù)已經(jīng)設(shè)計(jì)好的實(shí)驗(yàn)步驟，編寫程序并在matlab環(huán)境下運(yùn)行，3個(gè)波段原始圖像、噪聲分布圖以及去噪后的對(duì)比圖像展示如圖2所示.顯然，如果噪聲分布不是明顯的劃痕式或者條帶式分布，那么很難直觀地對(duì)比去噪前后影像的差異.因此，提取算法輸出橫軸為波段數(shù)，縱軸為亮度信息的曲線(如圖3(b)、圖4(b)、圖5(b)所示)，曲線上每一點(diǎn)代表相同點(diǎn)在不同波段上的亮度表達(dá)，從曲線走勢上可以清楚看出濾波前后差別.

選取第53波段中一點(diǎn)(39,7)如圖3，為第39行第7列的數(shù)據(jù)點(diǎn)，十字叉絲代表選點(diǎn)，下同.

圖2 濾波前后對(duì)比以及噪聲分布圖

圖3 第53波段處理示例

圖4 第49波段處理示例

圖5 第51波段處理示例

如圖3所示，前50波段曲線走勢變得平滑，但是由于波段窗口的選擇(此次實(shí)驗(yàn)選取的大波段值為11)，后面幾個(gè)數(shù)據(jù)沒有得到較好的處理，這是算法的不足之處；還有一種情況就是數(shù)據(jù)本身的原因，沒有濾掉噪聲.

圖4和圖5為第49波段和第51波段平滑處理后可以較好地保留圖像信息，曲線走勢保留較好，圖像信息較為完整.

從圖3(b)、圖4(b)、圖5(b)可以看出，算法對(duì)于高光譜圖像的最大作用體現(xiàn)在最大程度地保留曲線走勢；由于反射率變換后的圖像光譜曲線受各種噪聲影響而成鋸齒狀，經(jīng)過處理以后圖像光譜平滑了許多(從示例平滑數(shù)據(jù)圖像中可以明顯看出)；同時(shí)保持了一些如吸收峰等細(xì)微的光譜特征.

2.3 精度評(píng)價(jià)

為了與高光譜其他的去噪方法比較，應(yīng)用相似函數(shù)[23]：1)光譜角度距離：夾角越小，光譜間相似度越大；2)光譜歐幾里得距離：距離值越小，兩光譜越相似；3)光譜相關(guān)系數(shù)：光譜間相關(guān)系與光譜相似度成正比，完全相似的光譜間的相關(guān)系數(shù)為1.通過計(jì)算得到如表1.

表1 三種方法處理后的相似函數(shù)

MNF變換通過降低數(shù)據(jù)的維數(shù)來達(dá)到隔離數(shù)據(jù)噪聲的目的，EFFORT方法通過統(tǒng)計(jì)高光譜圖像，尋找一種微調(diào)來減小誤差，從而提高反射率的精度.二者都是光譜域?yàn)V波的經(jīng)典算法，但是由于計(jì)算量和閾值的選取，本文的算法無需降維，閾值選取也較為簡單，在運(yùn)算效率方面有所改觀.由表1計(jì)算結(jié)果看，相似度越大，表明在最大程度上保留了原光譜的信息，因此在保證去除高光譜圖像中的噪聲的同時(shí)，相似度的判定可以比較算法的改進(jìn).基于光譜二階導(dǎo)數(shù)的實(shí)驗(yàn)濾波方法都要優(yōu)于其他兩種方法，與實(shí)測光譜更加接近，濾波效果更佳.

3 結(jié) 論

根據(jù)光譜特征，基于高光譜影像提出一種利用光譜曲線二階導(dǎo)數(shù)的SG濾波去除噪聲的算法，得出以下結(jié)論：

1)基于光譜二階導(dǎo)數(shù)的濾波算法能最大程度保留曲線走勢，有效地去除光譜域噪聲，保留光譜原有的大部分光譜特征，是高光譜圖像預(yù)處理的一種簡單且有效的手段.

2)相比較其他算法，本文采用的算法花費(fèi)的時(shí)間也較少，相似性函數(shù)反映濾波后光譜與實(shí)測光譜相似度最大.

當(dāng)然，算法也存在著不足.實(shí)驗(yàn)算法存在著邊界效應(yīng)，針對(duì)兩端波段不能進(jìn)行較好的處理；實(shí)驗(yàn)僅對(duì)比了實(shí)測光譜(含噪聲)和濾波后光譜，并未在此基礎(chǔ)上添加不同類型噪聲，因此在適用性方面還有所欠缺.