結合Beltrami流和遞歸濾波的高光譜圖像分類方法

廖建尚 王立國 郝思媛

(1.廣東交通職業技術學院計算機工程學院, 廣州 510650； 2.哈爾濱工程大學信息與通信工程學院, 哈爾濱 150001;3.青島理工大學通信與電子工程學院, 青島 266520)

0 引言

通過成像光譜儀可以獲得達到幾百個波段光譜信息的高光譜遙感圖像，具有波段多、光譜分辨率高的優點，但是分類結果中椒鹽噪聲嚴重影響高光譜圖像分類的應用[1-2]。近年來，越來越多學者對在高光譜圖像分類中融入空間特征進行了深入研究，目前常見的空間特征提取方法有，形態濾波特征提取[3-4]、馬爾科夫隨機場特征提取[5-6]、圖像分割特征提取[7-8]，其中用濾波器提取高圖像紋理信息進行分類逐漸成為一個研究熱點。有部分學者用Gabor濾波器[9]提取紋理信息來輔助高光譜分類，HE等[10]通過識別低階Gabor濾波提取高光譜圖像的空間特征，提出了判別式低階Gabor濾波(Discriminative low-rank gabor filtering, DLRGF)的高光譜圖像分類方法；IMANI等[11]采用灰度共生矩陣、Gabor濾波器和形態學濾波器提取空間紋理信息、外形特征和像元鄰域信息，并通過不同的特征結合尋找最優的分類算法；JIA等[12]利用Gabor濾波器提取多維數據集特征，聯合稀疏表示提出一種有效地空譜結合分類方法GS-MTJSRC；馮逍等[13]按照紋理對不同地物響應的區分度逐一加入三維Gabor紋理特征并輔助光譜信息，用支持向量機(Support vector machine, SVM)進行分類。也有一些學者用雙邊濾波器[14]來提取高光譜空間紋理信息，如WANG等[15]先用概率SVM獲取初始分類結果，然后用雙邊濾波進行平滑處理，最后用圖像分割的方法獲取高光譜的分類結果，提高了分類精度。KANG等[16]先用SVM初始分類，再用雙邊濾波對分類結果進行優化，提出了邊緣保持濾波分類方法(Edge-preserving filtering, EPF)； SHEN等[17]用雙邊濾波器對空譜信息子集提取空間紋理信息，融合后再用極端學習機器對其進行有效分類。還有一些學者用遞歸濾波[18]提取空間紋理信息，如KANG等[19]用遞歸濾波對高光譜子集融合后提取空間信息，交由SVM實現分類，提出了圖像融合和遞歸濾波的高光譜圖像分類方法(Image fusion and recursive filtering, IFRF)，獲得較好的分類結果。另外有不少學者采用了Beltrami流濾波[20-24]進行了研究，Beltrami流應用廣泛，但目前沒有應用在高光譜圖像處理方面的研究報道。

過去在高光譜圖像空間紋理信息提取并分類的研究中取得了一定成效，但也存在一些不足：高光譜圖像空間分辨率低，空間紋理信息還有待進一步挖掘；地物分布往往具有統一性，在提取空間紋理信息時容易去除像元間的空間相關性。為提高高光譜圖像的分類性能，本文提出一種結合Beltrami流濾波和域轉換遞歸濾波的高光譜圖像分類算法(BFRF-SVM)，充分利用Beltrami流濾波提取的邊緣結構特征和域轉換遞歸濾波提取的空間相關信息，有效輔助SVM進行分類。

1 空間特征提取方法

1.1 Beltrami流濾波提取空間紋理信息

Beltrami流濾波是一種定義在圖像流形上的熱擴散方程，有清晰的幾何意義，可以獲取較好的圖像邊緣信息，對有n個波段的高光譜圖像，第i個波段圖像Ri的Beltrami流濾波在高光譜的定義為

(i=1,2,…,n)

(1)

其中

(2)

(3)

g——度量張量，表示圖像曲面的第一基本形式的行列式

Beltrami流濾波是一種最小化擴散的有效保持圖像邊緣、結構特點的噪聲濾波方法，在圖像的噪聲平滑與邊緣保持方面取得很好的平衡，式(1)的偏微分方程實現基于有限差分離散化，采用歐拉向前差分逼近It和中間偏差逼近空間導數，重新獲得Beltrami流濾波公式

(4)

圖1所示為印第安林高光譜圖像第10波段以及經PCA降維后的第1、2、3主成分濾波前后圖像，圖2所示為Beltrami流濾波對印第安林高光譜圖像第10波段以及經PCA降維后的第1、2、3主成分濾波后圖像，由圖1、2可以看出，Beltrami流濾波有較好的去噪效果，并能有效獲取高光譜邊緣結構特性。

圖1 印第安林高光譜圖像Fig.1 Hyperspectral images of Indian Pines

圖2 印第安林的Beltrami流濾波Fig.2 Beltrami flow filtering for Indian Pines

Beltrami流濾波具有較好的圖像邊緣保持特性，但是在保持邊緣的同時容易去除空間相關性信息，本文將用域轉換遞歸濾波[18]來獲取空間相關性信息，彌補Beltrami流濾波在圖像處理中的不足。

1.2 域轉換遞歸濾波提取空間相關信息

域轉換遞歸濾波是2011年由Gastal等提出的圖像特征提取算法，可以將二維的濾波轉換為一維的特征提取濾波，對第i個波段的高光譜圖像的域轉換遞歸濾波Di(n)在高光譜上的定義為

Di(n)=(1-ad)I[n]+adDi[n-1]

(5)

其中

d=f(yn)-f(yn-1)

(6)

(7)

(8)

(9)

(10)

式中ad——反饋系數I[n]——高光譜圖像Di[n-1]——第n-1次遞歸濾波d——兩個相鄰像元yn和yn-1在轉換域Ωw中的距離

f(yn)——某高光譜波段圖像Ri求偏微分后的累加求積分變換為遞增函數，將二維圖像濾波轉換為一維圖像濾波

r——濾波半徑

σs——空間標準偏差

σr——范圍標準偏差

σHt——第t次迭代的數值

N——總迭代次數

該濾波是一個指數衰減的無限沖擊響應濾波，當d增大即兩像元距離增大時，ad趨于0，阻止鏈接傳播，因此可以保持邊緣特性，式(5)并不對稱，依賴輸入和輸出，是一個因果濾波器，為了濾波對稱，式(5)需要執行2次，先從左邊到右邊，再從右邊到左邊，或者從頂部到底部，再從底部到頂部。

圖3所示為域轉換遞歸濾波對印第安林高光譜圖像第10波段以及經PCA降維后的第1、2、3主成分濾波后的圖像，由圖3可以看出，濾波后的圖像同類地物灰度相近，說明有一定的相關性。

圖3 印第安林的域轉換遞歸濾波Fig.3 Domain transform recursive filtering for Indian Pines

高光譜像元間有較強的空間相關性，為了分析Beltrami流濾波和域轉換遞歸濾波對空間相關性的保持特性，本文引入Moran于1950年提出的Moran’sI指數[25-26]來衡量其高光譜圖像濾波前后的空間自相關性，表達式為

(11)

I大代表空間相關性強，I小代表空間差異性大。為了對比兩類濾波器對空間自相關性的保持特性，分別對印第安林數據集濾波并計算濾波后的Moran’sI指數且求其均值。如圖4所示，域轉換遞歸濾波的大部分Moran’sI指數均值要高于Beltrami流濾波和光譜波段，因此域轉換遞歸濾波可以有效彌補Beltrami流濾波在圖像濾波中的不足。

圖4 印第安林數據集的Moran’s I指數的均值Fig.4 Average of Moran’s I for Indian Pines

1.3 兩種空間信息自適應融合分類

本文用Beltrami流濾波提取的空間紋理信息和域轉換遞歸濾波提取的空間相關性信息，兩種信息提取和融合后，由SVM完成分類，形成BFRF-SVM分類算法，實現過程如下:

(1)高光譜PCA降維：對于有l個波段的高光譜數據集R進行PCA降維，選擇前面維度的數據組成新的數據集B。

B=Pca(R)

(12)

(2)提取空間信息1：用Beltrami流濾波對PCA降維后的數據集H按照式(4)進行濾波，獲取空間信息Dmp。

(13)

(3)域轉換遞歸濾波提取空間信息2：用域轉換遞歸濾波對高光譜數據集R按照式(5)作域轉換標準卷積濾波，獲取空間信息Dts。

(14)

(4)空間信息線性疊加：Dmp和Dts按照

W=Dmp+Dts

(15)

結合成W。

(5)分類：①從空間信息數據集W中以一定比例隨機抽取訓練集Ws，其余部分作為測試集Wt。②用徑向基函數支持的SVM方法交叉驗證，尋找最優參數組合。③用徑向基函數支持的SVM對Ws進行訓練，獲取訓練模型。④獲取模型后，用徑向基函數支持的SVM對測試集Wt進行分類。

BFRF-SVM分類算法流程如圖5所示。

圖5 BFRF-SVM算法流程圖Fig.5 Flow chart of BFRF-SVM

2 實驗數據集

為驗證本文提出算法的有效性，采用印第安農林和帕維亞大學的2個高光譜數據集進行驗證，其特征如下：①印第安農林高光譜數據集來自機載可見紅外成像光譜儀(Airborne visible infrared imaging spectrometer)，是1992年在印第安納州西北部印第安農林收集到的高光譜遙感圖像，具有20 m的空間分辨率，包含144×144個像元，220個波段，由于噪聲和水吸收等因素除去其中的20個波段，剩余200個波段，共有16種地物。②帕維亞大學數據集來自反射式光學系統成像光譜儀(Reflective optics system imaging spectrometer)，拍攝于帕維亞大學的高光譜遙感圖像，包含610×340個像元，115個波段，由于噪聲等因素去除其中的12個波段，剩余103個波段，共有9個類別。

3 實驗設計與評價指標

為驗證本文所提算法在高光譜分類的優越性，使用7種方法進行比較：①利用文獻[27]提出的SVM算法。②利用PCA對光譜信息降維后，用SVM進行分類。③用Gabor濾波器、雙邊濾波器和導向濾波器分別對高光譜數據用PCA[27]降維后的前20個主成分提取空間信息，并將獲取的空間信息和光譜信息線性結合后，用SVM進行分類，形成SGB-SVM、SBL-SVM、SGD-SVM。④文獻[15]提出的EPF算法對高光譜圖像進行分類，有EPF-B-c和EPF-G-c。⑤文獻[19]提出的域轉換遞歸濾波方法IFRF。⑥文獻[29]提出的形態學濾波方法SMP-SVM。⑦BFRF-SVM。

本文用整體分類精度(Overall accuracy, OA)、平均分類精度(Average accuracy, AA)以及Kappa統計系數(Kappa statistic, Kappa)來衡量分類算法的精度，為了避免隨機偏差，每個實驗重復10次記錄平均結果，驗證平臺采用Matlab R2012b，i7-6700 CPU，8GB RAM實驗平臺。

3.1 參數尋優驗證

為了確定Beltrami流濾波參數迭代次數k和時間步長ht，本文采用窮舉法，通過前期實驗初步估算，限定取值范圍，分別為k∈[10,11,…,100]，ht∈[0.01,0.02,…,0.2]，采用印第安農林圖像作為實驗數據集，分別獲取最優分類結果的參數值。先對印第安林數據集進行PCA降維后取前30個主成分作為測試數據集，然后隨機抽取數據集的8%作為訓練樣本和數據集的92%作為測試樣本進行驗證，實驗結果如圖6所示。當k=36、ht=0.04，得到分類結果最優：整體分類精度為94.23%，因此本文Beltrami流濾波系數取k=36、ht=0.04。

為了確定域轉換遞歸濾波空間標準偏差σs、范圍標準偏差σr、總迭代次數N，使得分類效果達到最優，采用窮舉法，限定域轉換遞歸濾波3個參數的取值范圍，分別為σs∈[10,14,…,250]，σr∈[0.10,0.11,…,0.20]，N∈[3,4,…,10]，同樣采用印第安農林圖像作為實驗數據集，獲取最優分類結果的參數值。隨機抽取實驗數據集的8%作為訓練樣本和實驗數據集的92%作為測試樣本進行驗證，實驗結果如圖7所示，當σs=170，σr=0.18，N=6時，得到的分類結果最優：整體分類精度為94.94%，因此域轉換遞歸濾波系數取σs=170，σr=0.18，N=6。

圖6 Beltrami流濾波分類參數尋優Fig.6 Parameter optimization of Beltrami flow filtering for classification

圖7 域轉換遞歸濾波分類參數尋優Fig.7 Parameter optimization of domain transform recursive filtering for classification

圖8 印第安農林數據集分類Fig.8 Classification for data set of Indian Pines

3.2 實驗驗證

用7種方法分別對印第安林和帕維亞大學數據集進行分類對比驗證，其中驗證方法如下：①印第安農林數據集圖像，地物分布如圖8a所示，選取全部16個類別，每類隨機選取7%的樣本組成有標簽訓練集，其余93%作為測試集，地物數量較少的3類隨機抽取21%作為訓練集。表1為各種分類方法對印第安農林數據集的分類精度統計，分類效果如圖8所示。②帕維亞大學數據集圖像，地物如圖9a所示，選取全部9個類別，每類隨機選取3%樣本組成訓練集，其余97%作為測試集,表2為各種分類方法對帕維亞大學數據集的分類精度統計，分類效果如圖9所示。

表1 印第安農林數據集圖像分類數據統計Tab.1 Classification statistics for data set of Indian Pines

圖9 帕維亞大學數據集分類Fig.9 Classification for data set of Pavia University

地物樣本總數訓練/%測試/%SVMSVM-PCASGB-SVMSBL-SVMSGD-SVMEPF-B-cEPF-G-cIFRFSMP-SVMBFRF-SVM柏油馬路664139788.6888.7088.4590.1196.2497.8298.1790.3992.2597.36草地1864939795.9695.8495.7997.8698.73100.0099.9299.5298.8299.75砂礫209939759.2058.9565.0164.5580.9969.7174.2886.2881.8487.64樹木306439788.1887.3790.4894.8896.9994.8193.7488.9697.2697.77金屬板134539799.0998.8198.7298.4598.2999.8799.9299.2197.7498.89裸土502939758.5260.9762.3083.4292.3465.9467.5499.8491.2798.19瀝青屋頂133039768.9466.8076.1172.9089.3187.5487.9194.1381.5591.54地磚368239781.1277.2179.4882.3091.1098.8096.6580.4189.6191.83陰影94739788.7985.7392.6085.9592.5194.5094.2465.4486.5098.39OA/%85.8985.5886.8390.7595.5093.1893.3694.1694.3397.46AA/%80.9480.0483.2285.6092.9489.8990.2689.3590.7695.71Kappa80.9780.5882.3087.6694.0290.7891.0392.2592.4696.64

3.3 實驗分析

(1)用BFRF-SVM對2種數據集進行分類，其中印第安農林數據集的OA為96.01%，帕維亞大學數據集的OA為97.46%，比SVM整體分類精度高出12～15個百分點，比PCA-SVM整體分類精度高出12～16個百分點，比SGB-SVM、SBL-SVM、SGD-SVM高出2～12個百分點，比EPF算法高出4～5個百分點，比IFRF算法高出1～3個百分點，比SMP-SVM算法高出2～6個百分點，充分驗證了BFRF-SVM算法的有效性。

(2)印第安農林數據集的訓練樣本僅為7%，分類精度OA超過96%，可知BFRF-SVM對分布復雜的高光譜數據有較好的效果；帕維亞大學數據集分類訓練樣本僅為3%，分類精度OA超過97%，由圖8、9可以看出，算法對數據集的分類效果提升明顯，有效去除了椒鹽現象。

(3)為驗證監督數據對算法的影響，選擇不同的訓練樣本測試算法的分類精度。由圖10a可以看出，印第安農林總體分類精度在4%的訓練樣本比例時，OA超過了91%，10%的訓練樣本時，OA超過了97%；由圖10b可以看出，帕維亞大學總體分類精度OA在訓練樣本為0.5%時超過了90%，在4.5%的訓練樣本時OA超過了98%，驗證了BFRF-SVM算法在低訓練樣本的情況下也能得到較優的總體分類精度，且有一定的穩定性。

圖10 不同訓練樣本比例分類后OA、AA和Kappa折線圖Fig.10 Curves of OA、AA and Kappa with different training samples for classification

4 結束語

提出了用Beltrami流濾波和域轉換遞歸濾波結合的高光譜圖像分類算法(BFRF-SVM)，通過一系列實驗，驗證了BFRF-SVM算法優越性。實驗結果表明，本文提出的BFRF-SVM算法與單純光譜信息SVM算法、空譜結合分類算法、邊緣保持分類算法、遞歸濾波分類算法、形態學濾波算法相比，OA、Kappa和AA都有較大幅度提高，說明了用Beltrami流濾波和域轉換遞歸濾波提取空間信息結合后，能有效提高高光譜圖像的分類精度，該算法有以下特點：①用兩種濾波器將高光譜圖像光譜信息提取為空間信息并進行分類，有效利用了邊緣結構信息和空間相關性信息。②Beltrami流濾波提取高光譜圖像的邊緣結構信息，有效輔助SVM進行像元分類。③域轉換遞歸濾波獲得高光譜圖像的空間相關信息，彌補Beltrami流濾波中的不足，結合后有效提高了高光譜圖像的分類精度。

BFRF-SVM算法為高光譜遙感提供了一種有效的分類方法，可以應用于農作物的精準分類與識別領域。

1 童慶禧, 張兵, 張立福. 中國高光譜遙感的前沿進展[J]. 遙感學報, 2016, 20(5): 689-707. TONG Qingxi, ZHANG Bing, ZHANG Lifu. Current progress of hyperspectral remote sensing in China[J]. Journal of Remote Sensing, 2016,20(5): 689-707. (in Chinese)

2 王俊淑, 江南, 張國明,等. 高光譜遙感圖像DE-self-training半監督分類算法[J/OL]. 農業機械學報, 2015, 46(5):239-244. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20150534&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2015.05.034. WANG Junshu, JIANG Nan, ZHANG Guoming, et al. Semi-supervised classification algorithm for hyperspectral remote sensing image based on DE-self-training [J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2015,46(5):239-244. (in Chinese)

3 IMANI M, GHASSEMIAN H. Edge patch image-based morphological profiles for classification of multispectral and hyperspectral data[J]. IET Image Processing, 2017, 11(3):164-172.

4 BAO R, XIA J, MURA M D, et al. Combining morphological attribute profiles via an ensemble method for hyperspectral image classification[J]. IEEE Geoscience & Remote Sensing Letters, 2017, 13(3):359-363.

5 CHEN P, NELSON J D B, TOURNERET J Y. Toward a sparse bayesian markov random field approach to hyperspectral unmixing and classification[J]. IEEE Transactions on Image Processing, 2017, 26(1): 426-438.

6 CAO X, XU L, MENG D, et al. Integration of 3-dimensional discrete wavelet transform and Markov random field for hyperspectral image classification[J]. Neurocomputing, 2017, 226:90-100.

7 ZHANG Z, PASOLLI E, CRAWFORD M M, et al. An active learning framework for hyperspectral image classification using hierarchical segmentation[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2016, 9(2): 640-654.

8 CUI B, MA X, XIE X, et al. Classification of visible and infrared hyperspectral images based on image segmentation and edge-preserving filtering[J]. Infrared Physics & Technology, 2017, 81:79-88.

9 SHEN L, BAI L. Mutual Boost learning for selecting Gabor features for face recognition[J]. Pattern Recognition Letters, 2006, 27(15):1758-1767.

10 HE L, LI J, PLAZA A, et al. Discriminative low-rank gabor filtering for spectral-spatial hyperspectral image classification[J]. IEEE Transactions on Geoscience & Remote Sensing, 2017, 55(3):1381-1395.

11 IMANI M, GHASSEMIAN H. GLCM, gabor, and morphology profiles fusion for hyperspectral image classification[C]∥ Iranian Conference on Electrical Engineering, 2016.

12 JIA S, HU J, XIE Y, et al. Gabor cube selection based multitask joint sparse representation for hyperspectral image classification[J]. IEEE Transactions on Geoscience & Remote Sensing, 2016, 54(6):3174-3187.

13 馮逍, 肖鵬峰, 李琦,等. 三維Gabor濾波器與支持向量機的高光譜遙感圖像分類[J]. 光譜學與光譜分析, 2014, 34(8):2218-2224. FENG Xiao, XIAO Pengfeng, LI Qi, et al. Hyperspectral image classification based on 3-D Gabor filter and support vector machines[J].Spectroscopy and Spectral Analysis,2014,34(8):2218-2224.(in Chinese)

14 TOMASI C, MANDUCHI R. Bilateral filtering for gray and color images[C]∥ICCV, 1998:839-846.

15 WANG Y, SONG H, ZHANG Y. Spectral-spatial classification of hyperspectral images using joint bilateral filter and graph cut based model[J]. Remote Sensing, 2016, 8(9):748.

16 KANG X, LI S, BENEDIKTSSON J A. Spectral-spatial hyperspectral image classification with edge-preserving filtering[J]. IEEE Transactions on Geoscience & Remote Sensing, 2014, 52(5):2666-2677.

17 SHEN Y, XU J, LI H, et al. ELM-based spectral-spatial classification of hyperspectral images using bilateral filtering information on spectral band-subsets[C]∥IGARSS 2016—2016 IEEE International Geoscience and Remote Sensing Symposium, 2016:497-500.

18 GASTAL E S L, OLIVEIRA M M. Domain transform for edge-aware image and video processing[C]∥ACM Transactions on Graphics (TOG), 2011, 30(4): 69：1-69：12.

19 KANG X, LI S, BENEDIKTSSON J A. Feature extraction of hyperspectral images with image fusion and recursive filtering[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(6): 3742-3752.

20 SOCHEN N, DERICHE R, PEREZ L L. The Beltrami flow over implicit manifolds[C]∥IEEE International Conference on Computer Vision, 2003:832.

21 SHI Y, LAI R, WANG D J, et al. Metric optimization for surface analysis in the Laplace-Beltrami embedding space[J]. IEEE Transactions on Medical Imaging, 2014, 33(7):1447-1463.

22 FERNANDEZ J J, MARTINEZ J A. Three-dimensional feature-preserving noise reduction for real-time electron tomography[J]. Digital Signal Processing, 2010, 20(4):1162-1172.

23 MISHRA K N. An efficient Laplace-Beltrami Spectra based technique for online iris image compression and identification[C]∥2016 International Conference on Control, Computing, Communication and Materials (ICCCCM), 2016: 1-6.

24 TAN M, QIU A. Spectral laplace-beltrami wavelets with applications in medical images[J]. IEEE Transactions on Medical Imaging, 2014, 34(5):1005-1017.

25 MORAN P A P. The interpretation of statistical maps[J]. Journal of the Royal Statistical Society. Series B (Methodological), 1948, 10(2): 243-251.

26 MORAN P A P. Notes on continuous stochastic phenomena[J]. Biometrika, 1950, 37(1/2): 17-23.

27 MELGANI F, BRUZZONE L. Classification of hyperspectral remote sensing images with support vector machines[J]. IEEE Transactions on Geoscience and Remote Sensing, 2004, 42(8): 1778-1790.

28 CAI S, DU Q, MOORHEAD R, et al. Noise-adjusted principal component analysis for hyperspectral remotely sensed imagery visualization[C]∥16th IEEE Visualization Conference(Vis2005),2005:119-120.

29 FAUVEL M, BENEDIKTSSON J A, CHANUSSOT J, et al. Spectral and spatial classification of hyperspectral data using SVMs and morphological profiles[J]. IEEE Transactions on Geoscience & Remote Sensing, 2008, 46(11):3804-3814.