?

何 軍1 李滋剛2

（1.南京信息工程大學電子與信息工程學院 南京 210044）（2.東南大學儀器科學與工程學院 南京 210096）

1 引言

近年來，隨著傳感器技術和計算機技術的飛速發(fā)展，人們對于空間數(shù)據(jù)分析逐漸由全色波段、多光譜的柵格數(shù)據(jù)轉向更高維的高光譜影像。然而，光譜數(shù)據(jù)維度的提高卻給空間數(shù)據(jù)的分析帶來了極大的困難，這就是數(shù)據(jù)分析領域的維度災難問題（Curse of Dimensionality）［1］——即在保持一定取樣密度所需樣本的個數(shù)隨著數(shù)據(jù)空間維度的增加成指數(shù)級增加，而且散布在外圍的樣本數(shù)也隨著空間維數(shù)的增加成指數(shù)級增加，并且還表現(xiàn)為模型的復雜程度或模型的表示長度隨維數(shù)成指數(shù)級增加。為了克服維度災難給空間數(shù)據(jù)處理帶來的困難，我們希望能夠用低維的坐標以近似等價的方式去刻畫原來的高維空間中的數(shù)據(jù)。PCA（Principle components analysis）［2］和 MDS（Multi－dimensional scaling）［3］等線性算法能夠有效發(fā)現(xiàn)線性子空間，且在數(shù)據(jù)內插、外推、壓縮、去噪以及可視化等方面的應用已經被證明是非常有效的。由于PCA方法假設數(shù)據(jù)的固有結構是線性的，即PCA的本質是建立在全局線性假設之上，而高光譜數(shù)據(jù)中含有較多的非線性［4～6］，數(shù)據(jù)的分布呈現(xiàn)明顯的彎曲，不能再近似地認為高光譜數(shù)據(jù)取自線性子空間。因此，若要對高維空間數(shù)據(jù)進行更為準確的分析，需要研究發(fā)現(xiàn)非線性結構的方法。

近年來美國海軍研究實驗室的Bachmann博士等學者，首先基于非線性低維流形思想對于高光譜數(shù)據(jù)進 行 建 模［7～9］，并 取 得 了 非 常 顯 著 的 成果［10］。在Bachmann的工作啟發(fā)下，從擴散過程的角度出發(fā)，利用擴散映射思想［11～12］對高光譜數(shù)據(jù)的低維表示也展開了研究工作［13～14］。相對于傳統(tǒng)的線性表示技術，非線性的數(shù)據(jù)表示技術帶來的問題是較高的計算復雜性，難以處理的海量的數(shù)據(jù)。目前較為流行的一種方法是首先將大尺度的高光譜影像進行分片，對每個分片分別進行非線性降維，從而得到各自在流形上的內蘊坐標表示，然后選擇一個分片作為基準，并基于坐標系對準的思想，將其它分片的流形坐標映射到這個基準坐標系統(tǒng)中［9，13］。本文提出一種新的基于局部線性重構的擴散坐標系延拓算法，即對于海量的外采樣點（Out－of－sample），直接通過局部線性重構算法找到它們在擴散坐標系的坐標表示，這種新的方法相對于基于坐標系配準的延拓方法，高光譜影像數(shù)據(jù)在擴散坐標系下的表示更為一致。

2 擴散映射與局部線性重構算法

2.1 擴散映射

通過研究圖或流形上的擴散過程，進而發(fā)現(xiàn)蘊含在高維數(shù)據(jù)集中的信息，是Yale大學擴散幾何小組近年來發(fā)展起來的對于高維數(shù)據(jù)分析的理論［11］。擴散映射（Diffusion Maps）［12］算法簡單描述如下：

把數(shù)據(jù)集視為由數(shù)據(jù)作為頂點做構成的鄰接圖，數(shù)據(jù)與數(shù)據(jù)之間的邊的權重k（x，y）滿足如下的對稱半正定性質：

·k是對稱的：k（x，y）＝k（y，x）；

·k是正保持的：對于數(shù)據(jù)集X中所有的x和y，都有k（x，y）≥0；

·k是半正定的：即對于定義在X上所有的實有界函數(shù)f，均有：

這里μ是X上的概率測度。

對于對稱化后的擴散算子?A和擴散核?a（x，y），有如下的譜分解：

其中特征值滿足λ1＝1≥λ2≥λ3≥…。對于第m次擴散時刻的擴散算子?Am及對應的擴散核?a（m）（x，y），則有

因此，可以引入如下的一族擴散映射｛Φm｝：

對于高光譜影像在擴散幾何坐標下的可視化表示而言，該族擴散映射能夠將高光譜數(shù)據(jù)集Γ中任一x嵌入至低維的擴散空間：

其中j0是在不引入明顯誤差情況下，所嵌入擴散空間的維度。由擴散映射理論可知，擴散空間把握了特征空間內蘊幾何結構，因此我們將此低維坐標表示稱之為高光譜數(shù)據(jù)的擴散幾何坐標［13］。

2.2 局部線性重構算法

對于規(guī)模諸如O（106）的全景高光譜影像構造擴散映射，直接計算無論是運行內存還是運算時間均不可行。因此一種有效的解決方法是首先在少量采樣（In－Sample）的高光譜數(shù)據(jù)集上構造擴散映射，然后建立起由外采樣點（Out－of－Sample）到擴散幾何坐標的映射關系，利用該映射能夠近似得到外采樣點的擴散坐標。這里我們引入局部線性重構算法（Local Linear Reconstruction）能夠建立起外采樣點和已經構造好的擴散坐標系之間的映射。局部線性重構算法的示意圖如圖1所示。

對于高維數(shù)據(jù)集合Γ?RD，假定Γ分布在RD中的低維流形M 之上或附近，將Γ分為兩個子集合Γ1和Γ2，其中Γ1＝｛x1，x2，…，xn｝，Γ2＝｛xn＋1，xn＋2，…，xN｝∈RD，由Γ1構造出擴散映射Φ：RD→Rd，其中d≤D，使得yi＝Φ（xi），i＝1，2，…，n，這里yi對應于xi在擴散空間中的擴散坐標。相對于x1，x2，…，xn，xn＋1，xn＋2，…，xN稱之為外采樣點。對于Γ2的數(shù)據(jù)點xj，j＝n＋1，…，N，無法直接套用Φ得到對應的yj，因為Φ并沒有顯示的表達式，它蘊含在特征值求解過程。由于Γ1和Γ2分布在相同的低維流形M之上或附近，因此借助于局部線性假設，我們可以把Γ2中任一數(shù)據(jù)點xj放入Γ1中，并找到xj在Γ1中的局部鄰域｛xj1，xj2，…，xjk｝，在該局部鄰域內xj可由｛xj1，xj2，…，xjk｝線性表示，而線性表示的系數(shù)｛Wj1，Wj2，…，Wjk｝在嵌入空間保持不變，于是yj能夠利用此重構系數(shù)計算得出。具體算法描述見算法1。

圖1 局部線性重構算法

算法1 局部線性重構算法

Step 1.對集合Γ1構造擴散映射Φ：RD→Rd，得到Γ1在Φ 作用下的像｛y1，y2，…，yn｝；

Step 2.將Γ2中每一數(shù)據(jù)點xj在Γ1中進行最近鄰居搜索，如kNN搜索，得到Γ1中與xj在歐式距離意義下的若干最近鄰居｛xj1，xj2，…，xjk｝；

Step 3.根據(jù)xj在Γ1中的近鄰｛xj1，xj2，…，xjk｝，計算線性重構的權值Wjk，使得

Step 4.由計算所得的重構權值Wjk和鄰域對應的擴散坐標｛yj1，yj2，…，yjk｝，重構出yj，即：

3 基于局部線性重構的擴散坐標系延拓算法

在以往的工作中［13］，利用局部線性重構方法可以得到分片中的部分點在骨干（Backbone）的擴散坐標系下的近似坐標表示，進而構造了分片擴散坐標系與Backbone擴散坐標系之間的坐標變換矩陣，但是由于局部線性重構方法構造的兩個坐標系之間的近似對應關系，因此坐標變換的方法并不能保證把各個分片精確地統(tǒng)一到Backbone的擴散坐標系，所以在最終擬合的全景高光譜擴散坐標中分片與分片之間存在不連續(xù)性。

在我們最近的研究工作中，發(fā)現(xiàn)只要構造了Backbone并計算出Backbone的擴散坐標系，那么對于全景高光譜影像中的任何數(shù)據(jù)點，均能用局部線性重構的方法在Backbone的擴散坐標系下重構出該點對應的擴散幾何坐標，所以如果對Backbone之外所有的數(shù)據(jù)點做同樣的計算，不僅能夠成功地得到全景高光譜影像的擴散幾何坐標表示，而且解決了不連續(xù)性的問題。因此，本節(jié)提出了基于局部線性重構的擴散坐標系延拓算法，該算法的關鍵之處在于兩點，第一是Backbone的構造，第二則是把除了Backbone之外的數(shù)據(jù)點，即外采樣點，“插入”至Backbone的擴散坐標系下，“插入”的方法就是“局部線性重構算法”。

3.1 骨干坐標系（Backbone）的構造

在本擴散坐標延拓算法中，僅需直接計算出Backbone的擴散幾何坐標，外采樣點的擴散幾何坐標不需直接算出，這些外采樣點將由局部線性重構算法被“插入”至Backbone的擴散坐標系。因此Backbone的選擇非常重要。

從局部線性重構算法知道Backbone的選擇取決于兩個因素，即Backbone對高光譜特征空間的覆蓋性和計算可行性。這里引入一個假設，即全景高光譜影像的特征空間相對樣本空間的稀疏性假設。我們知道，全景高光譜影像是遙感設備對一個地區(qū)地表物質光譜吸收和反射情況的測量。樣本空間的大小size（S）對應于遙感設備的空間分辨率，而特征空間的大小size（F）取決于觀測地區(qū)地表成分的復雜程度。很顯然，對于O（106）分辨率而言，樣本空間的大小size（S）＝O（106），而絕不會出現(xiàn)該地區(qū)的地表任何一處都是不同的成分，即size（F）?size（S）。

為了使得采樣S′能夠較好的覆蓋S，我們引入ε－net的概念，即定義在S中的子集S′，所有點均滿足ε－net要求，也就是說S′中的任何兩個點xi、xj之間的距離均大于ε，而對于S中的其它的點y，在S′中存在至少一個點x到y(tǒng)的距離小于ε。ε－net的概念在圖2做了直觀的解釋?；讦牛璶et可以設計出更為有效的隨機采樣算法，即算法2隨機εnet采樣方法。

圖2 ε－net示意圖

算法2 隨機ε－net采樣方法

采樣算法執(zhí)行的過程就是構造ε－net的過程。隨機選擇x0∈S；由歸納法，對于k≥1，假設已經選擇了x0，x1，…，xk∈S，使得所選擇的任何兩個點之間的距離均不小于ε。令

表示采樣k個點后，在S中挖去k個ε－球后的剩余空間。如果Rk為空集，則算法停止，否則繼續(xù)從Rk中選擇點xk＋1。

由Rk的定義可知，xk＋1與x0，x1，…，xk之間的距離必然不小于ε。因此當構造ε－net的過程結束時，也就是說當?shù)趉＊個點被選中后，那么對于任意的y∈S我們都能在S′中找到一個點x，使得這兩個點之間的距離不超過ε，否則的話則有y∈Rk，上述構造過程不會停止。因此，算法2保證了隨機采樣很好地覆蓋了原樣本空間。

3.2 外采樣點在骨干坐標系上的局部線性重構

基于局部線性重構算法，依據(jù)高維數(shù)據(jù)集中嵌入低維流形結構的假設和局部線性思想，能夠得到外采樣點在低維流形上的參數(shù)化表示。另外，由于我們采用ε－net隨機采樣所構成的Backbone能夠較好地逼近全景高光譜影像的特征空間，因此Backbone同樣能夠較好逼近全景高光譜影像的低維流形，計算得出的Backbone擴散幾何坐標則是對該低維流形結構的參數(shù)化表示。于是對于Backbone之外的外采樣點，每一個點均可以利用局部線性重構算法得到該點在Backbone擴散坐標系下的坐標表示，算法的示意圖見圖3，具體算法步驟見算法3。

圖3 外采樣點在Backbone上的局部線性重構

算法3 基于局部線性重構的擴散坐標系延拓算法

Step 1.使用算法2介紹的ε－net隨機采樣方法對規(guī)模為｜S｜的全景高光譜影像S進行采樣，構造出S的Backbone S′，并使用二維數(shù)組B記下S′在S 中的標號（i，j）；

Step 2.對Backbone S′構造擴散映射ΦS′，并由此得到Backbone的擴散坐標系｛φ1，φ2，…，φd｝S′；

Step 3.對于外采樣點SO∶＝S＼S′中每一個點x，執(zhí)行如下子過程，并同時將SO在S中的標號記入二維數(shù)據(jù)B；

Step 3.1 利用近似最近鄰居搜索算法在S′中搜索x的k個最近鄰居x1，x2，…，xk，求解如下方程

從而得到在k－鄰域內x的局部線性表示；

Step 4.將重構出的SO的擴散幾何坐標與Backbone S′的擴散幾何坐標按照其在全景高光譜影像S中的標號順序進行排列，即完成該算法。

4 實驗效果

針對NASA所公開的AVIRIS系統(tǒng)所采集的高光譜影像進行實驗。這里選擇AVIRIS系統(tǒng)在1997年在Moffet地區(qū)采集的高光譜影像的第二景（Scene）數(shù)據(jù)。該景影像的分辨率是614×512，含有224個近似連續(xù)的波段，波段范圍從0.369μm～2.506μm，根據(jù)AVIRIS的說明，在這224個波段中存在大約30個完全吸收和噪聲波段，在此實驗中將此30個波段排除。實驗平臺是Intel Core i5CPU 2.3G，8G內存，Mac OSX 10.6操作系統(tǒng)，基于 Matlab2010b所實現(xiàn)的算法。

4.1 骨干（Backbone）坐標系的實驗確定

圖4 AVIRIS高光譜影像ε－net隨機采樣ε選擇與采樣比率的曲線

圖5 ε設置為0.48時Backbone在波段1和波段60的投影

首先需要確定的是Backbone。這里采用ε－net隨機采樣方法，因為ε－net在較小的采樣率情況下，能夠很好地覆蓋原樣本空間，所以對AVIRIS高光譜影像進行ε－net隨機采樣實驗，確定合適的ε。對于該高光譜影像，ε取值與所采樣數(shù)據(jù)點的規(guī)模曲線如圖4所示，可以看到當我們將ε減小到0.6時，采樣規(guī)模迅速上升，對應于采樣的密度越來越大。這里我們將ε設置為0.48，圖5展示了采樣所得Backbone在波段1和波段60的分布情況，能夠很好地逼近原高光譜特征空間在這兩個波段的投影。

4.2 AVIRIS高光譜影像的全景可視化

當ε設置為0.48時我們采樣所得Backbone由6827個點組成，由擴散映射算法［12］，計算Backbone的擴散坐標系。由于一次性讀入全景高光譜影像會造成龐大的內存開銷，因此將高光譜影像劃分為7×7個分片，分別逐次地應用局部線性重構算法，求出所有外采樣點在Backbone擴散坐標系下的坐標表示，這里同樣設置最近鄰居搜索的參數(shù)k＝10，擴散空間的維度為15。基于算法XX將Backbone擴散坐標系延拓至全部的外采樣點，得到全景高光譜影像的擴散幾何坐標表示，其中進行ε－net隨機采樣的耗時為10.5s，計算Backbone的擴散幾何坐標的耗時為3.1s，通過局部線性重構算法將外采樣點延拓至Backbone擴散坐標系的耗時為57.6s。由以上參數(shù)設置計算得出在擴散坐標系下高光譜影像的表示如圖6和圖7所示，其中圖6是由擴散幾何坐標（φ2，φ3，φ4）所構成的RGB彩色影像，圖7由擴散幾何坐標（φ5，φ6，φ7）所構成的RGB彩色影像。

圖6 由擴散幾何坐標（φ2，φ3，φ4）所構成的RGB彩色影像

實驗表明，基于局部線性重構的擴散坐標延拓算法能快速地計算出全景高光譜影像的擴散幾何坐標，解決了擴散幾何坐標難以處理海量規(guī)模的全景高光譜影像問題。然而，與我們之前基于坐標系配準的方法對比［13］，發(fā)現(xiàn)同樣是由擴散幾何坐標（φ5，φ6，φ7）構成的 RGB彩色影像，表示效果卻大相徑庭，造成這樣結果的原因是擴散幾何坐標的計算本質上是數(shù)據(jù)驅動（Data－driven）的，在構造擴散映射時選擇不同的數(shù)據(jù)點就會產生不同的擴散幾何坐標。因此，從實際應用角度而言，更傾向于本文所提出的基于局部線性重構的擴散坐標系延拓算法。因為算法3構造的基礎是經過精心選擇的Backbone，一旦Backbone確定，那么所延拓的擴散坐標系也能夠保持一致，在實際應用中就可以依據(jù)Backbone的擴散坐標系對原有的基于線性的高光譜分析方法做相應的改進。同時，算法3本質上是流形半監(jiān)督學習的一個特例，隨著機器學習領域流形半監(jiān)督學習方法的不斷完善，對于海量高光譜影像能夠得出更為精確的擴散幾何坐標。

圖7 由擴散幾何坐標（φ5，φ6，φ7）所構成的RGB彩色影像

5 結語

直接用擴散幾何坐標表示全景高光譜影像存在計算不可行的困難，由于全景高光譜影像通常達到計算機系統(tǒng)難以直接處理的龐大規(guī)模，使得構造擴散映射的幾個關鍵步驟難以進行。本文通過借鑒局部線性嵌入LLE算法［15］的思想，設計了局部線性重構算法，并將該算法應用于全景高光譜影像擴散坐標表示算法。

從應用價值而言，全景高光譜影像利用擴散幾何坐標進行表示，能夠有效地發(fā)現(xiàn)某一地區(qū)環(huán)境的物種分布情況，對于實施大范圍的環(huán)境監(jiān)測具有重要的意義。如果能夠得到某一地區(qū)充分多的高光譜影像，利用ε－net隨機采樣就能夠得到該地區(qū)準確的Backbone，進而就可以將此Backbone的擴散幾何坐標作為該地區(qū)各個物種的“指紋”特征，將對日后的應用帶來極大的便利。

［1］Bellman R E.Dynamic Programming［M］.New Jersey：Princeton University Press，1957.

［2］Jolliffe T.Principal Components Analysis［M］.Berlin：Springer，1986.

［3］Cox T F，Cox M A A.Multidimensional Scaling［M］.London：Chapman and Hall，2001.

［4］Goodin D G，Gao J C，Henebry G M.The effect of solar illumination angle and sensor view angle on observed patterns of spatial structure in tallgrass prairie［J］.IEEE Transactions on Geoscience and Remote Sensing，2004，42（1）：154－65.

［5］Sandmeier S R，Middleton E M，Deering D W，et al.The potential of hyperspectral bidirectional reflectance distribution function data for grass canopy characterization［J］.Journal of Geophysical Research－Atmospheres，1999，104（D8）：9547－9560.

［6］Keshava N，Mustard J F.Spectral unmixing［J］.IEEE Signal Processing Magazine，2002，19（1）：44－57.

［7］Bachmann C M，Ainsworth T L，F(xiàn)usina R A.Modeling data manifold geometry in hyperspectral imagery［C］／／IGARSS 2004：IEEE International Geoscience and Remote Sensing Symposium Proceedings，Vols 1－7－Science for Society：Exploring and Managing a Changing Planet，2004：3203－3206.

［8］Bachmann C M，Ainsworth T L，F(xiàn)usina R A.Exploiting manifold geometry in hyperspectral imagery［J］.IEEE Transactions on Geoscience and Remote Sensing，2005，43（3）：441－454.

［9］Bachmann C M，Ainsworth T L，F(xiàn)usina R A.Im－proved manifold coordinate representations of hyperspectral imagery［C］／／Proceedings of the IGARSS 2005：IEEE International Geoscience and Remote Sensing Symposium，Vols 1－8，Proceedings，2005：4307－4310.

［10］Bachmann C M，Ainsworth T L，F(xiàn)usina R A，et al.Bathymetric Retrieval From Hyperspectral Imagery Using Manifold Coordinate Representations［J］.IEEE Transactions on Geoscience and Remote Sensing，2009，47（3）：884－897.

［11］Coifmfman R R，Lafon S，Lee A B，et al.Geometric diffusions as a tool for harmonic analysis and structure definition of data：Diffusion maps［J］.Proceedings of the National Academy of Sciences of the United States of America，2005，102（21）：7426－7431.

［12］Coifman R R，Lafon S.Diffusion maps［J］.Applied and Computational Harmonic Analysis，2006，21（1）：5－30.

［13］HE J，Zhang L，WANG Q，et al.Using Diffusion Geometric Coordinates for Hyperspectral Imagery Representation［J］.IEEE Geoscience and Remote Sensing Letters，2009，6（4）：767－771.

［14］何軍.空間數(shù)據(jù)低維幾何表示及其應用研究［D］.南京：東南大學儀器科學與工程學院，2009.

［15］Roweis S T，Saul L K.Nonlinear dimensionality reduction by locally linear embedding［J］.Science，2000，290（5500）：2323－2325.