何英杰 謝東海 鐘若飛
(首都師范大學(xué)資源環(huán)境與旅游學(xué)院, 北京 100048)
經(jīng)歷全色到多光譜掃描儀成像階段之后,在20世紀(jì)80年代初由于成像光譜概念的出現(xiàn),遙感技術(shù)進(jìn)入了高光譜遙感的嶄新的階段.高光譜遙感是當(dāng)前遙感的前沿技術(shù)之一,是將成像技術(shù)和光譜技術(shù)相結(jié)合的多維信息獲取的技術(shù).它是由多光譜遙感成像技術(shù)的基礎(chǔ)發(fā)展而來的,光譜分辨率得到了大幅的提高[1],為每個像元提供幾乎連續(xù)的地物光譜曲線,使我們利用高光譜反演陸地細(xì)節(jié)成為可能[2].
獲取較高的光譜分辨率的連續(xù)的圖像數(shù)據(jù)是高光譜遙感數(shù)據(jù)的閃耀之處,它能充分地反映地物反射率光譜的細(xì)微特征.然而,高光譜圖像成像儀器以及環(huán)境等因素對數(shù)據(jù)獲取造成的影響,使得在高光譜數(shù)據(jù)成像和預(yù)處理的各個環(huán)節(jié)中,都可能引入不同程度的噪聲,高光譜圖像具有空間圖像和地物光譜兩個方面的信息,因而噪聲對高光譜圖像的影響最終也會表現(xiàn)為空間域噪聲和光譜域噪聲兩方面[3],嚴(yán)重地限制了高光譜數(shù)據(jù)的高光譜分辨率的優(yōu)勢.
人們根據(jù)實際圖像的特點、噪聲的統(tǒng)計特征和頻譜分析的規(guī)律,發(fā)展了各式各樣的去噪方法.其中最直觀的方法是根據(jù)噪聲能量一般集中于高頻,而圖像頻譜分布于一個有限區(qū)間的這一特點,采用低通濾波方法.在基于三角網(wǎng)格模型的降噪算法研究方面,Yu等[4]提出一種兩階段特征保持網(wǎng)格去噪框架.此外,對圖像進(jìn)行平滑處理也是常用的方法.平滑的目的有兩個:改善圖像質(zhì)量和抽出對象特征.一個較好的平滑方法應(yīng)該既可以消掉噪聲影響又不使圖像的邊緣輪廓和線條變模糊.劉倩[5]基于非局部稀疏的圖像處理方法,充分利用自然圖像本身的有用信息,有效彌補(bǔ)了數(shù)學(xué)模型的缺陷,以數(shù)據(jù)驅(qū)動的方式大大提高了圖像平滑的效果.圖像平滑處理方法有空間域處理和頻率域處理兩大類.在空間域里一類方法是噪聲去除,即先判定某點是否為噪聲點,若是,重新賦值,若不是按原值輸出;另一類是平均,即不依賴于噪聲點的識別和去除,而對整個圖像進(jìn)行平均運(yùn)算.而在頻率域里是對圖像頻譜進(jìn)行修正,一般采用低通濾波方法,而不像在空間域里直接對圖像像素灰度級值進(jìn)行運(yùn)算.圖像的均值濾波模板和中值濾波模板是實現(xiàn)圖像的平滑去噪的主要方法.
由于高光譜圖像成像時自然光照明條件影響、地面地形影響、混合像元問題等原因引入各種噪聲,這些噪聲存在于光譜域中,同樣也是需要檢測并予以去除的[6].MNF(minimum noise fraction)變換通過降低數(shù)據(jù)的維數(shù)來達(dá)到隔離數(shù)據(jù)噪聲的目的,首先使用噪聲的協(xié)方差矩陣分解并將數(shù)據(jù)中的噪聲白化[7],此時噪聲具有單位方差,并且與原波段無相關(guān)性.然后針對上述噪聲白化的數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)的主成分變換;Boardman[8]針對這個問題提出了EFFORT方法,即通過統(tǒng)計高光譜圖像,尋找一種微調(diào)來減小誤差,從而提高反射率的精度;Gladden[9]同樣認(rèn)為反射率轉(zhuǎn)換后光譜中會產(chǎn)生很大的噪聲,并采用一種 SSE線性校正方法,即將轉(zhuǎn)換后反射率光譜擬合到實測地物光譜上,達(dá)到去除噪聲的目的.陳志剛和束炯[10]基于光譜二階導(dǎo)數(shù)檢測噪聲的基礎(chǔ)上,提出了基于經(jīng)驗?zāi)B(tài)分解EMDF(empirical mode decomposition based filter)的光譜域濾波方法及其修正,具有較好的光譜域濾波效果.由于多分辨率的分析特性和較好的時頻局部化,小波分析[11]依據(jù)小波基函數(shù),多層次地分解高光譜圖像光譜.然后針對分解后的小波系數(shù)選取合適的閾值進(jìn)行處理,最終獲得濾除噪聲后的光譜曲線,在圖像處理領(lǐng)域得到廣泛應(yīng)用.楊浩等[12]針對成像光譜儀采集的高光譜圖像數(shù)據(jù),提出一種小波閾值去噪方法去降低噪聲影響.遙感影像存在的光譜鋸齒噪聲[13],同樣也嚴(yán)重影響礦物光譜吸收參數(shù)的提取.
Savitzky和Golay[14]提出的最小二乘濾波算法,濾波器的選擇為Savitzky-Golay(SG) 濾波器,它是一種特殊的低通濾波器,最初由Savitzky和Golay于1964年提出.光譜曲線由于受各種噪聲的影響而呈現(xiàn)鋸齒狀,經(jīng)SG濾波處理后,平滑一些細(xì)微的鋸齒噪聲同時能保持整條曲線大的光譜特征不受影響[15-16].本文主要以航空高光譜遙感數(shù)據(jù)為研究對象,在SG濾波器的基礎(chǔ)之上進(jìn)行優(yōu)化,求取反射率光譜的二階導(dǎo)數(shù),以期使其能更加靈活且準(zhǔn)確地處理影像.
在平面坐標(biāo)系中,用一條曲線來擬合一組數(shù)據(jù),不妨假設(shè)這條曲線為y=a0+a1x+a2x2+a3x3+a4x4,當(dāng)每一個點的橫坐標(biāo)代入這個曲線方程后,所得值與該點的縱坐標(biāo)之差的平方之和最小時,這條曲線的擬合度最高,從而可以確定所有的系數(shù)ai(i=0,1,2,3,4).
從高光譜影像中每個像元提取出的反射率曲線本身有一定的物理意義,可以反映地物的細(xì)微特征,但是轉(zhuǎn)換后的像素光譜反射率曲線的噪聲一般呈鋸齒狀分布,這些鋸齒分布即為高光譜圖像光譜域噪聲,嚴(yán)重地限制了其高光譜分辨率的優(yōu)勢.為了去除光譜域噪聲,考慮采用多項式擬合的最小二乘平滑算法,剔除一些細(xì)微的鋸齒噪聲同時能保持整條曲線大的光譜特征不受影響.SG濾波法[14]是一種經(jīng)典的最小二乘平滑算法,它使用簡化的最小二乘擬合卷積方法對曲線進(jìn)行平滑處理并可計算平滑后曲線的各階導(dǎo)數(shù).
擴(kuò)展一下,通常假設(shè)曲線為p次多項式,即滿足公式(1)
Yi=a0+a1i+a2i2+…+apip
(1)
其中,Yi第i點平滑后的數(shù)值.
用上述多項式擬合曲線的誤差可以表示為:見公式(2)
(2)
其中,yj代表平滑前的數(shù)值,平滑窗口大小k=2m+1.為了讓誤差S最小,需要對S進(jìn)行偏微分,使得所求S偏微分各項為零.因此,可以得到以下各式:關(guān)系如公式(3)
…
(3)
Savitzky對整個求解過程進(jìn)行了推導(dǎo),給出了平滑窗中心點平滑后數(shù)值的最后公式及公式中系數(shù)的計算方法[14].Madden[17]基于Savitzky的計算公式,改正了原SG系數(shù)的一些錯誤,并且給出了各階導(dǎo)數(shù)平滑系數(shù)修正后的計算公式.根據(jù)其推導(dǎo)平滑后的曲線數(shù)值的最小二乘卷積方程可寫為:見公式(4)
(4)
(5)
式中,m為平滑窗寬度的一半[18].
SG濾波器有如下幾大優(yōu)點[19]:
1)利用最小二乘的多項式擬合方法非常清晰易懂,并且在計算上來說,多項式卷積的操作比最小二乘的計算可操作性更強(qiáng);
2)濾波系數(shù)只需要在對應(yīng)的卷積系數(shù)表中進(jìn)行查找,很容易獲得;
3)SG濾波器可以有任意的長度,因此有利于采樣頻率通常很低的生物學(xué)或者生物力學(xué)的數(shù)據(jù)處理.
SG算法的本質(zhì)是一種最小二乘法的卷積平滑,根據(jù)多項式擬合階數(shù)和平滑次數(shù)進(jìn)行綜合考量,但是多項式擬合曲線的誤差求取的是一階導(dǎo)數(shù).本文高光譜影像噪聲檢測SG算法是基于光譜二階導(dǎo)數(shù).光譜導(dǎo)數(shù)技術(shù)對光譜信噪比非常敏感,光譜的一階、二階導(dǎo)數(shù)計算公式如下[20-21]:
(6)
(7)
其中,近似選取相鄰兩點,其關(guān)系滿足公式(8)
Δλ=λk-λj=λj-λi,λk>λj>λi
(8)
光譜的低階導(dǎo)數(shù)處理對噪聲影響敏感性較低,而高階微分對噪聲影響敏感度高.相比較一階導(dǎo)數(shù),二階導(dǎo)數(shù)很好的反映了噪聲的實際分布情況(高階導(dǎo)數(shù)過于復(fù)雜).因此,可以考慮用光譜二階導(dǎo)數(shù)對光譜曲線進(jìn)行噪聲影響程度檢測.
對于光譜曲線進(jìn)行光譜二階導(dǎo)數(shù)計算,考慮將其與平均值進(jìn)行相減,然后一般想法是參考差值分布,與一個固定的值進(jìn)行比較,即設(shè)定一個閾值,判定噪聲與否.然后根據(jù)SG濾波法通過對反射率光譜二階導(dǎo)數(shù)的統(tǒng)計,得到光譜各波段中噪聲的判定函數(shù),在噪聲較嚴(yán)重的波段進(jìn)行較多點數(shù)的SG濾波,對噪聲不太嚴(yán)重的波段進(jìn)行較少點數(shù)的SG濾波,從而確定該波段濾波平滑窗的大小.根據(jù)計算表明,如果判別式滿足公式(9):
(9)
那么該波段反射率的值認(rèn)為是有較大噪聲的存在,其中σ表示光譜二階導(dǎo)數(shù)的標(biāo)準(zhǔn)差.
在噪聲檢測完成后,可以先采取較大窗口的SG濾波對較大噪聲波段進(jìn)行局部第一次平滑,其他波段保持不變,得到初步的濾波結(jié)果;然后采取較小窗口的SG濾波對初次濾波的結(jié)果進(jìn)行光譜曲線的整體平滑.這樣的好處是,既可以對光譜曲線進(jìn)行平滑濾波和噪聲去除,又可以最大程度的對光譜細(xì)節(jié)進(jìn)行保留.
本算法的具體步驟如下:
1)從高光譜圖像中逐像元的提取光譜特征,繪制反射率曲線;
2)計算曲線的二階導(dǎo)數(shù);
3)計算光譜二階導(dǎo)數(shù)的均值與方差,按照公式(9)判定噪聲;
4)進(jìn)行第一次平滑,對判定為噪聲的波段進(jìn)行例如11點SG濾波平滑,保持其他波段不變;
5)進(jìn)行第二次平滑,對第一次平滑后的光譜曲線進(jìn)行例如5點(一般選擇小窗口的值為大窗口值的一半取整)的SG平滑;
6)該點處理完畢,轉(zhuǎn)入下一像素點的處理,轉(zhuǎn)入步驟2,直到所有點處理完畢.
試驗數(shù)據(jù)為2012年武漢市某地區(qū)PHI高光譜遙感影像,包含80個波段,光譜范圍400~850 nm,大小為128行、128列.選取比較有代表性的波段進(jìn)行展示,例如49、51、53波段,如圖1所示.
圖1 49、51、53波段影像
圖像質(zhì)量評價的研究是圖像信息學(xué)科的基礎(chǔ)研究之一.對圖像處理或圖像通信系統(tǒng),其信息的主體是圖像,衡量這個系統(tǒng)的重要指標(biāo),就是圖像的質(zhì)量,而圖像去噪就是為了改變圖像的主觀視覺顯示圖像質(zhì)量[22].我們需要考慮的幾個因素總結(jié)如下:1)去噪后圖像應(yīng)盡量的平滑,不存在或有較少的噪聲痕跡.2)去噪結(jié)果不能使圖像過渡的失去結(jié)構(gòu)細(xì)節(jié)而變得模糊.3)沒有由于具體去噪方法產(chǎn)生的人工噪聲.4)均方差盡可能小.所有這些都要求有一個合理的圖像質(zhì)量評價方法.根據(jù)圖像質(zhì)量評價方法的思想,初步對所給數(shù)據(jù)進(jìn)行分析:波段存在劃痕狀條帶噪聲,以及分布不均勻的隨機(jī)噪聲,以49波段噪聲特點較為明顯.根據(jù)光譜域處理噪聲的結(jié)果與一般平滑去噪的不同,要盡量保證噪聲點附近平滑處理,而不影響非噪聲點光譜特性.基于高光譜影像應(yīng)用SG濾波算法,可以滿足預(yù)期設(shè)想.
根據(jù)已經(jīng)設(shè)計好的實驗步驟,編寫程序并在matlab環(huán)境下運(yùn)行,3個波段原始圖像、噪聲分布圖以及去噪后的對比圖像展示如圖2所示.顯然,如果噪聲分布不是明顯的劃痕式或者條帶式分布,那么很難直觀地對比去噪前后影像的差異.因此,提取算法輸出橫軸為波段數(shù),縱軸為亮度信息的曲線(如圖3(b)、圖4(b)、圖5(b)所示),曲線上每一點代表相同點在不同波段上的亮度表達(dá),從曲線走勢上可以清楚看出濾波前后差別.
選取第53波段中一點(39,7)如圖3,為第39行第7列的數(shù)據(jù)點,十字叉絲代表選點,下同.
圖2 濾波前后對比以及噪聲分布圖
圖3 第53波段處理示例
圖4 第49波段處理示例
圖5 第51波段處理示例
如圖3所示,前50波段曲線走勢變得平滑,但是由于波段窗口的選擇(此次實驗選取的大波段值為11),后面幾個數(shù)據(jù)沒有得到較好的處理,這是算法的不足之處;還有一種情況就是數(shù)據(jù)本身的原因,沒有濾掉噪聲.
圖4和圖5為第49波段和第51波段平滑處理后可以較好地保留圖像信息,曲線走勢保留較好,圖像信息較為完整.
從圖3(b)、圖4(b)、圖5(b)可以看出,算法對于高光譜圖像的最大作用體現(xiàn)在最大程度地保留曲線走勢;由于反射率變換后的圖像光譜曲線受各種噪聲影響而成鋸齒狀,經(jīng)過處理以后圖像光譜平滑了許多(從示例平滑數(shù)據(jù)圖像中可以明顯看出);同時保持了一些如吸收峰等細(xì)微的光譜特征.
為了與高光譜其他的去噪方法比較,應(yīng)用相似函數(shù)[23]:1)光譜角度距離:夾角越小,光譜間相似度越大;2)光譜歐幾里得距離:距離值越小,兩光譜越相似;3)光譜相關(guān)系數(shù):光譜間相關(guān)系與光譜相似度成正比,完全相似的光譜間的相關(guān)系數(shù)為1.通過計算得到如表1.
表1 三種方法處理后的相似函數(shù)
MNF變換通過降低數(shù)據(jù)的維數(shù)來達(dá)到隔離數(shù)據(jù)噪聲的目的,EFFORT方法通過統(tǒng)計高光譜圖像,尋找一種微調(diào)來減小誤差,從而提高反射率的精度.二者都是光譜域濾波的經(jīng)典算法,但是由于計算量和閾值的選取,本文的算法無需降維,閾值選取也較為簡單,在運(yùn)算效率方面有所改觀.由表1計算結(jié)果看,相似度越大,表明在最大程度上保留了原光譜的信息,因此在保證去除高光譜圖像中的噪聲的同時,相似度的判定可以比較算法的改進(jìn).基于光譜二階導(dǎo)數(shù)的實驗濾波方法都要優(yōu)于其他兩種方法,與實測光譜更加接近,濾波效果更佳.
根據(jù)光譜特征,基于高光譜影像提出一種利用光譜曲線二階導(dǎo)數(shù)的SG濾波去除噪聲的算法,得出以下結(jié)論:
1)基于光譜二階導(dǎo)數(shù)的濾波算法能最大程度保留曲線走勢,有效地去除光譜域噪聲,保留光譜原有的大部分光譜特征,是高光譜圖像預(yù)處理的一種簡單且有效的手段.
2)相比較其他算法,本文采用的算法花費的時間也較少,相似性函數(shù)反映濾波后光譜與實測光譜相似度最大.
當(dāng)然,算法也存在著不足.實驗算法存在著邊界效應(yīng),針對兩端波段不能進(jìn)行較好的處理;實驗僅對比了實測光譜(含噪聲)和濾波后光譜,并未在此基礎(chǔ)上添加不同類型噪聲,因此在適用性方面還有所欠缺.