周子勇

（油氣資源與探測國家重點實驗室 中國石油大學 地球科學學院，北京 102200）

0 前言

高光譜遙感技術近年發(fā)展迅速，并形成了一系列數(shù)據(jù)處理與信息提取方法，這些方法主要可以歸為三類：①繼承并發(fā)展多光譜遙感影像處理方法，如通過波段比值、PCA分析提取信息；②針對高光譜影像數(shù)據(jù)特點發(fā)展起來的一系列方法，如MNF（Minimum Noise Fraction）、SAM（Spectral Angle Mapper）、SAI（Spectral Absorption Index）、SFF（Spectral Feature Fitting）、MTMF（Mixture－Tuned Matched Filtering）、混合像元解混技術等；③借鑒信號處理領域的新方法對影像進行處理和信息提取，如小波分析方法、分形分析方法、神經(jīng)網(wǎng)絡方法等。

Hilbert－Huang（HHT）變換方法是1998年由Huang［3］首先提出來的一種信號處理方法，用于現(xiàn)實中非線性、非平穩(wěn)信號的處理與分析。HHT分析方法的前提是經(jīng)驗模態(tài)分解。所謂經(jīng)驗模態(tài)分解就是通過一種特定的濾波方法把原始數(shù)據(jù)分解成由高頻（細節(jié)成份）到低頻（基本成份）的一系列本征模態(tài)函數(shù)。由于不同的本征模態(tài)函數(shù)分別刻劃了原始信號中不同目標的信號成份，因此借助于EMD分析，有可能從遙感影像中分離出不同目標對應的信號成份，進而提取出有用信息。把EMD方法用于遙感影像的分析處理，目前的研究多是對單一圖像進行分析。如通過對影像的二維EMD分解（BEMD），提取影像的邊界信息［13］、紋理特征［10，12］或者通過 EMD 進行圖像融合［1］，或者通過EMD分析不同信號成份的空間分布規(guī)律［4］等。這種單波段影像的EMD分析，沒有發(fā)揮高光譜影像具有豐富的光譜信息的優(yōu)勢。因此作者提出基于像元高光譜曲線進行經(jīng)驗模態(tài)分解，構(gòu)建IMF影像，并對數(shù)據(jù)特征進行分析。

作者在本文以Hyperion高光譜影像為試驗數(shù)據(jù)，對每一像元的高光譜曲線進行經(jīng)驗模態(tài)分解，形成一系列IMF影像，通過分析IMF影像的特點，探討運用EMD方法進行高光譜影像處理與信息提取的可能性。作者首先簡述EMD算法的基本原理，然后對比分析不同IMF影像的特點，最后討論影響高光譜影像經(jīng)驗模態(tài)分解的主要因素以及可能的解決辦法。

1 遙感影像的經(jīng)驗模態(tài)分解

1.1 經(jīng)驗模態(tài)分解算法

經(jīng)驗模態(tài)分解算法的基本思想是，通過連續(xù)去除原始信號的本征信號要素IMF，從而達到信號分解的目的。其基本算法如下［3］：

給定一信號x（t），對信號x（t）按下列步驟進行處理：

（1）找出信號x（t）的所有極大值點Maxi（i＝1，2，…）和極小值點Mink（k＝1，2，…）。

（2）根據(jù)極大值點和極小值點進行三次樣條插值，分別得到信號的上包絡線Max（t）和下包絡線 Min（t）。

（3）計算每一點上包絡線和下包絡線的平均值e（t），即e（t）＝ （Max（t）＋Min（t））／2。

（4）從信號中減去e（t），即：x（t）＝ x（t）－e（t）。

（5）轉(zhuǎn)回到步驟（1），當x（t）不變化時，停止。這時所得到的x（t）稱為第一個IMF，記為IMF1。

（6）從原始信號移去IMF1，即x（t）＝x（t）－IMF1，然后把相減后的結(jié)果作為輸入數(shù)據(jù)重復步驟（1）至步驟（5）流程，可以得到IMF2、IMF3…。

（7）當最終分解的IMF沒有極值點時（單調(diào)上升或下降），停止EMD過程。

從上述算法原理可知，經(jīng)驗模態(tài)分解的實質(zhì)是把原始信號分解成一系列調(diào)幅／調(diào)頻（AM／FM）信號之和，即每一個IMF都是一個AM／FM信號。與傅立葉變換和小波變換相比，經(jīng)驗模態(tài)分解所得到的信號，更加直觀，易于分析。然而，該方法畢竟是經(jīng)驗分析方法，在實際操作過程中還存在很多的問題，信號的分解受多個因素的影響，后面將對這個問題進行詳細討論。

1.2 高光譜遙感影像的經(jīng)驗模態(tài)分解與結(jié)果分析

1.2.1 數(shù)據(jù)來源與處理方法

本文實驗數(shù)據(jù)為2006年獲得的一景渤海灣地區(qū)的Hyperion影像數(shù)據(jù)。數(shù)據(jù)覆蓋區(qū)大部分為海域，作者截取了其中位于盤錦市附近的近海區(qū)域作為研究對象，地面覆蓋類型主要為水體、植被及裸露地表。Hyperion數(shù)據(jù)掃描帶寬度為7.5km，空間分辨率為30m，光譜分辨率約為10nm，覆蓋400nm～2 500nm的波長范圍，共有220個定標的波段。其中有波長重復波段以及噪聲特別大沒有有效信息的波段，去除這些波段后保留176個波段。表1為保留波段與原始波段對照表。本文中所指的波段號都是指176個波段的編號。

表1 保留波段與原始波段對照表Tab.1 Selected wavebands versus original wavebands

在進行EMD處理前，需要對影像進行預處理。由于作者主要目的是研究原始高光譜影像的經(jīng)驗模態(tài)分解效果，因此只對原始數(shù)據(jù)進行簡單的預處理，即把影像的DN值轉(zhuǎn)換為絕對輻射值。轉(zhuǎn)換方法很簡單，即把所有VNIR波段（1號～50號）除以“40”，所有SWIR波段（51號～176號）除以“80”，得到絕對輻射值。然后以絕對輻射值進行EMD處理。EMD處理以像元為基礎。由于每一像元都對應有一條高光譜曲線，因此可以按照類似于一維時間序列EMD分解方法進行高光譜曲線的分解，并分解得到一系列的IMF。對影像中的每一個像元都進行EMD處理，得到對應的IMF，然后把每一像元同一層次的IMF構(gòu)成IMF影像，得到IMF1、IMF2等一系列影像。由于不同地物像元高光譜曲線特征不同，分解得到的IMF數(shù)目也不一樣，實驗中首先對不同地物像元高光譜曲線進行經(jīng)驗模態(tài)分解試驗，分析其IMF特征。結(jié)果表明，每一像元可以分解的IMF個數(shù)在6個～7個之間。因此在對整個數(shù)據(jù)進行處理時，把每個像元的IMF個數(shù)設定為“8”，這樣既便于對不同像元的IMF進行對比，同時又盡可能減少計算量。整個計算過程通過ENVI／IDL編程實現(xiàn)。

1.2.2 結(jié)果分析

通過經(jīng)驗模態(tài)分解，可以得到一系列IMF影像，每一像元都有一系列的IMF曲線，同時每一個IMF也有一個對應的影像。這樣既可以分析每一個像元中每一個IMF的信號特征，又可以分析每一個IMF影像的空間特征。因此與原始高光譜影像數(shù)據(jù)相比，通過對不同IMF的對比分析，有可能突出更多的細節(jié)特征，得到一些其它信號分析方法無法得到的信息。下面通過對Hyperion經(jīng)驗模態(tài)分解結(jié)果的分析，來說明這一問題。

（1）光譜曲線的IMF特征。通過經(jīng)驗模態(tài)分解，把光譜曲線分解成一系列由高頻到低頻的IMF。不同層次的IMF反映了信號中不同的成份，或者說不同的信號來源，因此不同地物的高光譜曲線分解后的IMF也不同。圖1是影像中特定像元光譜曲線經(jīng)驗模態(tài)分解結(jié)果。在圖1中，IMF1反映了高光譜曲線的高頻特征，而IMF5則反映了曲線的總體趨勢。盡管從圖1上看，兩種地物其原始波譜曲線差異并不明顯，但分解的IMF差異卻相當明顯。根據(jù)IMF的特征差異，可以用于影像分類［11］。

（2）通過IMF分析影像的噪聲特征。一般來說，IMF1對應是光譜曲線的高頻成份，反映信號中的噪聲或細節(jié)成份。因此通過分析IMF1，可以了解高光譜影像中每一波段的噪聲水平和性質(zhì)。為了評價各波段的噪聲水平，分別計算出了每一波段原始影像及對應IMF1影像的方差，并通過其比值來評價原始影像的噪聲特征。如圖2所示為每一波段原始影像、對應的IMF1影像的方差及其比值。由圖2可知，在不同波段其比值大小表現(xiàn)出不同的特征。在28號～57號（對應的波長范圍為：701.55nm～993.17nm，下 同 ），1 號 ～3 號（426.82nm ～ 447.17nm），70 號 ～ 71 號（1 124.28nm ～1 134.28nm），93 號 ～98 號（1 426.94nm～1 477.43nm），131 號 ～132 號（1 810.38nm～1 941.57nm），138 號 ～139 號（2 002.06nm～2 012.15nm）波段區(qū)間，其比值相對較小，其它波段區(qū)間比值相對較高。

由于IMF1反映信號中的噪聲或細節(jié)成份，因此該比值大小在某種程度上反映了各波段噪聲的性質(zhì)。比值較大的波段，IMF1反映的主要是影像中的系統(tǒng)噪聲；而比值較小的波段，既有可能是原始影像的信噪比較低，也有可能是IMF1反映了影像的細節(jié)。這點通過對比圖3和圖4可以看得很清楚。圖3為第16波段原始影像及其對應的IMF1影像，圖4為138波段原始影像及其對應的IMF1影像。顯然在圖3中，IMF1反映的是系統(tǒng)噪聲，而圖4中的IMF1則反映了原始信號中較強的背景噪聲，由于背景噪聲大，原始影像中目標信息被掩蓋，而在IMF1中，目標信息被突出了。因此通過對比原始影像及對應的IMF1影像，不但可以知道每一波段的噪聲大小，而且還可以了解噪聲的性質(zhì)。

（3）Smile效應。未經(jīng)處理的高光譜影像存在Smile效應。如果影像存在Smile效應，則通過MNF處理后的影像，其灰度在空間上有明顯的線性分帶特征。這種特征在IMF中也可以表現(xiàn)得非常明顯。如圖5所示分別為IMF1的28號和29號波段影像，由圖5可見明顯的灰度空間分帶現(xiàn)象。

（4）與小波分析或FFT分析不同，經(jīng)驗模態(tài)分解后的結(jié)果，與原始信號的波段存在一一對應關系，因此原始影像和IMF影像，對比更直觀，解釋更容易。

3 EMD中幾個問題的探討

影響高光譜影像經(jīng)驗模態(tài)分解結(jié)果的因素很多，其中影響最大的因素主要有以下四個方面：①極值點的位置和極值的確定；②插值算法；③端點效應；④IMF終止條件。下面對上述影響因素進行討論。

3.1 極值點的位置和極值的確定

確定極值點的位置和極值的大小是經(jīng)驗模態(tài)分解的基礎。尋找原始信號的極值點，可以采用下述算法：

極大值點 Maxi：If （xi≥xi－1and xi＞xi＋1）or （xi＞xi－1and xi≥xi＋1）

極小值點 Mink：If （xi＜xi－1and xi≤xi＋1）or （xi≤xi－1and xi＜xi＋1）

這種算法計算簡單，但是由于只是簡單的比較相鄰值，因此原始數(shù)據(jù)中的噪聲有可能產(chǎn)生一些不必要的極值點。特別是對于遙感影像，除系統(tǒng)本身的誤差外，還受大氣、下墊面等因素的影響，噪聲來源較多。而且在做EMD分解前，只是把DN值轉(zhuǎn)換成了輻射絕對值，并未做去噪處理，因此因噪聲產(chǎn)生的極值點不可避免，在進行EMD處理時，應該去除這類極值點。為了解決這一問題，在確定極值點時，先計算相鄰兩極值點的極差，然后設定一閾值，把極差小于閾值的極值點去除。也有文獻［6］提出通過相鄰三個極值點插值的方法來解決。這種方法的有效性有待于進一步驗證，而且計算量大，特別是對于本文的大數(shù)據(jù)量，計算量會增加很多，為此作者采用簡單的相鄰值比較方法。

3.2 插值算法

插值是EMD的一個重要組成部分。Huang［3］提出采用三次樣條方法進行插值，根據(jù)極小值和極大值求得最大最小包絡線。由于不同插值方法得到的最大最小包絡線可能不完全一樣，因此分解的IMF也可能存在一定的差異。許多文獻討論了插值方法對EMD結(jié)果的影響［5，9］，如何根據(jù)實際數(shù)據(jù)選擇最優(yōu)插值方法，目前并沒有更多的研究和判別標準，而且也沒有一種普適的算法。但目前大多數(shù)的文獻都采用三次樣條方法來處理不同的數(shù)據(jù)，因此作者仍然采用三次樣條方法。

3.3 分解結(jié)束判別標準

EMD整個過程涉及兩個停止條件，①IMF的判別條件；②整個分解過程的結(jié)束條件。對于前者，Huang給出了兩個條件：一是IMF 中極值點的數(shù)目與曲線穿越零值點的個數(shù)相等或相差1個；二是IMF的上包絡線和下包絡線各點的平均值為零。顯然，第一個條件包含在第二個條件中。所以確定IMF的關鍵是如何給出第二個條件，由于實際計算中第二個條件很難精確確定，Huang建議通過前后兩次計算的均方差來判斷IMF的合理性，但這種方法有時難以達到好的效果［2］。

對于EMD的結(jié)束條件，可以根據(jù)最終分解的結(jié)果是否為單調(diào)曲線來確定。但是這一判別條件，有可能出現(xiàn)IMF冗余現(xiàn)象，即分解的IMF比實際的要多，從而不能反映信號的真實成份。為此，有文獻提出采用能量變化作為判別準則。無論哪種方法，都只是針對單一信號曲線而言。對于高光譜影像，由于需要對每一個像元的高光譜曲線進行一維經(jīng)驗模態(tài)分解，并最終形成IMF影像，因此必須保證每一個像元對應高光譜曲線分解后的IMF具有可比性。為了解決多信號IMF可比性問題，Rato提出通過一個所謂的分辨率因子的辦法，文中所說的分辨率因子，指的是上下包絡線平均值的能量與原始信號能量的比值。顯然分辨率因子值越小，分解的IMF個數(shù)越少；反之越大。因此通過調(diào)整每一個信號分解過程中的分辨率因子，就可以保證各信號分解的IMF個數(shù)相同，使相同級的IMF具有可比性。但文中并沒有具體給出分辨率因子的確定方法。實際上，當對未知信號進行分解時，分辨率因子的選擇本身就存著在判別標準的問題。

還有一種簡單的方法，就是把IMF的個數(shù)設定為一常數(shù)。當大致知道分解的IMF個數(shù)時，可通過這種方法簡單設定，而且也便于對比。通過對試驗數(shù)據(jù)進行經(jīng)驗模態(tài)分解，發(fā)現(xiàn)每一像元的高光譜曲線分解的IMF個數(shù)都是“6”或“7”，基于上述考慮，作者把IMF分解個數(shù)設定為“8”，經(jīng)計算結(jié)果表明，此結(jié)束判別條件不會對結(jié)果產(chǎn)生太大的影響。

3.4 端點效應

高光譜曲線的兩個端點一般都不是極值點，對曲線兩個端點的處理，會對分解的IMF產(chǎn)生影響，即端點效應問題。很多文獻討論了減少端點效應的解決辦法。常用方法主要有以下幾種：①依據(jù)端點與鄰近值的關系直接指定極值點，即如果與端點相鄰點為極小值點，則把該端點作為極大值；否則作為極小值。該方法簡單，但誤差較大；②根據(jù)與端點相鄰的兩個極值點向相反方向延展極值點［6］?；镜乃悸啡缦拢涸O信號第一個極大值和極小值分別為Max1、Min1，其對應的位置分別為tmax1和tmin1，則在－tmin1處插入極大值 Max1，在－tmax1處插入極小值Min1，按同樣的方法在信號末端拓展兩個極值，然后進行插值處理；③與第二種類似，只是采用更多的極值（總極值數(shù)的一半）進行拓展［8］；④分別以信號首端和尾端為對稱點把原始數(shù)據(jù)進行鏡像取值，形成一個三倍于原始數(shù)據(jù)的數(shù)據(jù)量，插值完成后，截取中間部分進行計算。除此之外，還有一些其它方法，但是這些方法都是基于實際數(shù)據(jù)的比較，并沒有一個判別標準。對于高光譜曲線，也可以采用上述方法消除端點效應。本文作者分別試驗了上述四種方法，結(jié)果認為第二種方法效果較好。

4 結(jié)論

作者在本文采用EMD方法對Hyperion影像進行處理，并對結(jié)果進行了初步分析。實驗結(jié)果表明，通過分析IMF，可以更好地描述遙感影像，了解影像的固有特征，因而可以基于IMF進行特征提取。但是，由于EMD方法本身存在不足，在基于EMD的高光譜影像處理時，還存在許多問題有待于進一步研究。另外對于大量的高光譜數(shù)據(jù)，其計算量也是一個必須考慮的問題。盡管如此，EMD方法還是有它的獨特之處，特別是經(jīng)驗模態(tài)分解結(jié)果，與原始數(shù)據(jù)可比性好，與其它分析（如小波分析、分形分析）相比，結(jié)果更直觀，因此有值得進一步研究的必要。作者只是初步討論了運用EMD分解來了解Hyperion影像的特征，下一步將綜合高光譜影像的一維和二維經(jīng)驗模態(tài)分解結(jié)果，研究不同地物與IMF的對應關系，以及基于IMF的高光譜影像特征提取方法。

［1］CHEN SHAOHUI，SU HONGBO，ZHANG RENHUA，et al.Fusing remote sensing images using atrous wavelet transform and empirical mode decomposition［J］.Pattern Recognition Letters.2008，29：330－342.

［2］CHENG J S，YU D J，YANG Y.Research on intrinsic mode function（IMF ）criterion in EMD method［J］.Mechanical Systems and Signal Processing，2006，20：817－824.

［3］HUANG N E，SHEN Z，LONG S R，et al.The em－pirical mode decomposition and Hilbert spectrum for nonlinear and non－stationary time series analysis［J］.Proceedings of the Royal Society of London A.1998，454：903－995.

［4］NUNES J C，BOUAOUNE Y，DELECHELLE E.Image analysis by bidimensional empirical mode decomposition［J］.Image and Vision Computing，2003，21，12：1019－1026.

［5］QIN S R，ZHONG Y M.A new envelope algorithm of Hilbert Huang transform［J］.Mechanical Systems and Signal Processing，2006，20：1941－1952.

［6］RATO R T，ORTOGIEORA M D，BATISTA A G.On the HHT，its problems，and some solutions［J］.Mechanical Systems and Signal Processing.2008，22：1374－1394.

［7］RILLING G，F(xiàn)LANDRIN P.On the influence of sampling on the empirical mode decomposition，in：Proceedings of the International Conference on A－coustics［J］.Speech and Signal Processing，ICASSP 2006：III－444－III－447.

［8］RILLING G，F(xiàn)LANDRIN P，GONcALVES P.On empirical mode decomposition and its algorithms，in：IEEE－EURASIP Workshop on Nonlinear Signal and Image Processing NSIP－03，Grado［C］.2003.

［9］XU Z G，HUANG B X，LI K W.An alternative envelope approach for empirical mode decomposition［J］.Digital Signal Processing，2010，20：77－84.

［10］ZHANG Y，SUN Z X，LI W H.Texture synthesis based on Direction Empirical Mode Decomposition［J］.Computers ＆ Graphics，2008，32：75－86.

［11］董士偉，周子勇，基于EMD與神經(jīng)網(wǎng)絡的油膜高光譜數(shù)據(jù)特征提?。跩］.遙感技術與應用，2010，25（2）：221－226.Dong Shiwei，Zhou Ziyong，Wen Baihong.Feature Extraction of Offshore Oil Slick from Hyperspectral Data Based on EMD and Neural Network［J］.Remote sensing technology and application，2010，25（2）：221－226.

［12］沈濱，崔峰，彭思龍.二維EMD的紋理分析及圖像瞬時頻率估計［J］.計算機輔助設計與圖形學學報，2005，17（10）：2345－2352.Shen Bin，Cui Feng，Peng Silong.Bidimensional EMD for Texture Analysis and Estimation of the Instantaneous Frequencies of an Image［J］.Journal of computer2aided design ＆computer graphics，2005，17（10）：2345－2352.

［13］萬建，任龍濤，趙春暉，二維EMD應用在圖像邊緣特征提取中的仿真研究［J］.系統(tǒng)仿真學報，2009，21（3）：799－822.Wan Jian，Ren Longtao，Zhao Chunhui.Simulation Research of Applying Two－dimensional Empirical Mode Decomposition on Image Feature Extraction［J］.Journal of System Simulation，2009，21（3）：799－822.