郭杰 李云松 弋方 吳憲云

(西安電子科技大學 綜合業(yè)務網(wǎng)理論及關(guān)鍵技術(shù)國家重點實驗室，陜西 西安 710071)

高光譜遙感是一門集航空航天、傳感器、計算機等技術(shù)于一體的新興交叉學科［1］，涉及電磁波理論、光譜學與色度學、物理/幾何光學、電子工程、信息學、地球科學等多門學科，是20 世紀80年代以來人類在對地觀測方面所取得的重大技術(shù)突破之一.由于高光譜遙感成像技術(shù)所獲取的地球表面圖像包含了豐富的空間、輻射和光譜三重信息，因而受到國內(nèi)外學者的廣泛關(guān)注，并已廣泛應用于地質(zhì)制圖、植被調(diào)查、大氣研究、環(huán)境監(jiān)測等領域，在監(jiān)測地球環(huán)境和地物分類識別中發(fā)揮著越來越重要的作用.

像元是記錄高光譜遙感器所獲取的地面反射或輻射光譜信號的基本單位.由于遙感器的空間分辨率不夠精細，像元往往包含多種地表覆蓋類型，形成混合像元.根據(jù)光譜混合模型的相關(guān)背景，混合像元可以看作不同地表類型光譜響應特征的綜合.光譜混合模型從本質(zhì)上分為線性混合模型(LMM)和非線性混合模型(NLMM)［2］，將像元分解為不同的“基本組分單元”或者“端元”，并將這些端元成分定量分解，就是“光譜解混”過程［3］.由于LMM 實現(xiàn)簡單，且在大部分實驗中與真實場景較為吻合，故目前對光譜解混模型和算法的研究大多集中于線性光譜解混［4］.

在線性光譜解混領域，端元提取是關(guān)鍵步驟之一，目前國內(nèi)外已經(jīng)有很多成熟的端元提取方法［5］.比較著名的端元提取方法有:基于單形體向量投影性質(zhì)的純像元指數(shù)(PPI)［6］算法，利用光譜/空間信息進行端元搜索的AMEE［7］、SSEE［8］等算法，假設高光譜圖像中存在純光譜特征的OSP［9］、N-FINDR［10］、VCA［11］等算法，以及不采用純光譜特征假設的MVSA［12］、SISAL［13］等算法.同時也出現(xiàn)了與PPI［14-15］、N-FINDR［16］和VCA［17］等算法對應的硬件實現(xiàn)結(jié)構(gòu).PPI 算法因其開放性和實用性而廣泛應用于遙感分析和端元提取前處理中.目前商用的遙感圖像處理軟件ENVI 也已經(jīng)將PPI 算法作為標準化流程集成到端元提取模塊中.

現(xiàn)有的PPI 硬件實現(xiàn)結(jié)構(gòu)［14-15］因其投影矩陣的運算量非常大、沒有高效優(yōu)化的矩陣運算策略而導致硬件資源占用率較高，處理時間較長，故難以滿足高光譜解混實時處理的要求.為此，文中提出了一種基于PPI 算法的FPGA 實現(xiàn)方案，采用投影向量并行的矩陣運算策略和譜段并行的內(nèi)積運算電路，以減少接口的數(shù)據(jù)讀取和提高計算并行度，滿足星載/機載應用下的實時性要求.

1 PPI 算法

1.1 線性光譜混合模型的幾何學描述

假設遙感器可在B 個譜段上成像，所獲取的高光譜圖像有N 個像元xi(i=1，2，…，N)，則xi是一個包含B 個譜段的B 維列向量;每個像元由M 個端元ej(j=1，2，…，M)按照一定的比例混合而成，ei也是一個B 維列向量，且第j 個端元在第i 個像元中的比例為aij;第i 個像元的誤差項為wi.則存在如下關(guān)系:

設高光譜圖像矩陣X=［x1x2… xN］，端元矩陣E=［e1e2… eM］，比例系數(shù)所形成的豐度矩陣A=［a1a2… aN］，誤差矩陣W=［w1w2… wN］，則式(1)的矩陣形式可表示為

由于每個像元光譜均為B 維向量，故每個像元可以對應B 維空間中的一個點.在忽略誤差的前提下，依據(jù)式(2)，每個像元xi均是端元{ej}Mj=1 在下的凸組合，故和構(gòu)成了RB中的一個凸集.又因為端元向量一般是線性無關(guān)的，且譜段數(shù)B 一般遠大于端元數(shù)M，即B≥M－1 總是成立，所以RB中包含的最小凸集是一個M－1維空間的凸面單形體，這些端元就在凸面單形體的M 個端點上.具有3 譜段的線性混合模型的幾何學描述如圖1 所示.

圖1 具有3 譜段的線性混合模型的幾何學描述Fig.1 Geometrical description of linear mixing model with three spectral bands

1.2 PPI 算法描述

凸面單形體的形狀決定了其在空間中任意直線上的投影必為線段，且線段端點為單形體頂點的投影.PPI 算法利用這一性質(zhì)，在特征空間中隨機生成若干直線，并將所有像元投影到各條直線上，那么直線上所有投影點中最靠外的兩個投影點便是端元的投影.但由于模型誤差和噪聲的存在，也可能出現(xiàn)非端元被投影到線段端點的位置上，因此，需要生成大量隨機直線進行投影，并為每個像元定義一個純像元指數(shù)，用于記錄其被投影到線段端點的次數(shù).顯然，某個像元所對應的純像元指數(shù)越大，說明其是端元的可能性越大.經(jīng)過足夠多次投影后，可以根據(jù)每個像元的純像元指數(shù)判定端元.圖2 為PPI 算法的投影示意圖，每形成一個投影向量，會有兩個點的純像元指數(shù)增加，圖中點旁邊的數(shù)字即為該點的純像元指數(shù).

圖2 PPI 算法的投影示意圖Fig.2 Schematic diagram of projection of PPI algorithm

2 PPI 算法的并行硬件實現(xiàn)

2.1 內(nèi)積運算電路

空間中點到直線的投影可以通過向量內(nèi)積運算完成.假設高光譜圖像矩陣X 與隨機投影向量γ(γ是一個包含B 個分量的列向量)的內(nèi)積為X'=γTX，則X'中的最大和最小值分別對應投影線段的兩個端點(如圖2 所示).由1.2 節(jié)可知，PPI 算法中需要生成若干隨機投影向量，隨著投影向量數(shù)量K的增多，凸面單形體頂點的純像元指數(shù)不斷增加.經(jīng)過足夠多次投影后，純像元指數(shù)最多的那些像元便可以認為是端元所在的位置.

單個像元xi(i=1，2，…，N)與某一投影向量γk(k=1，2，…，K)的內(nèi)積可以表示為，因此一次內(nèi)積運算包含B 次乘法和B－1 次加法運算.由于高光譜圖像的譜段數(shù)B 一般大于100，若γ中的分量在二值集合{1，－1}中選取，則可以保證純像元指數(shù)計算的有效性，而且可以避免硬件實現(xiàn)中浮點數(shù)乘法對速度和資源帶來的負面影響.在文獻［14-15］的PPI 硬件實現(xiàn)中，計算純像元指數(shù)時所用的投影向量數(shù)K 一般取104，文中實現(xiàn)的硬件結(jié)構(gòu)中所用的投影向量數(shù)K 也為104.

筆者提出的一種譜段并行的內(nèi)積運算電路如圖3所示.將單個像元xi與某一投影向量γk的內(nèi)積按譜段展開進行并行計算.即對于某個像元，將其B個譜段的值分別與投影向量的B 個分量相乘并求和，由于投影向量各分量的取值僅為1 或－1，故在一個時鐘周期內(nèi)可以同時完成B 次×1 或×(－1)運算，因此并行度是順序完成B 次×1 或×(－1)運算電路的B 倍.

圖3 譜段并行的內(nèi)積運算電路結(jié)構(gòu)Fig.3 Architecture of band-parallel inner-product computation circuit

2.2 矩陣運算策略

PPI 算法的實現(xiàn)實質(zhì)上轉(zhuǎn)化為隨機向量矩陣與高光譜圖像矩陣之間的乘法運算.矩陣的并行運算是降低處理時間的有效途徑.如圖4 所示的兩種并行運算策略:①像元向量并行策略，將各個像元向量并行與同一個投影向量進行內(nèi)積運算，可以計算出N 個像元在該投影向量下的結(jié)果;②投影向量并行策略，將各個投影向量同時與同一個像元向量進行內(nèi)積運算，可以計算出該像元在各個投影方向上的結(jié)果.比較這兩種運算策略，可以發(fā)現(xiàn):高光譜圖像數(shù)據(jù)非常龐大，按像元向量進行并行計算需要同時獲取所有像元數(shù)據(jù)，并且在與各投影向量計算時都需要保存圖像數(shù)據(jù)，這在硬件實現(xiàn)中對數(shù)據(jù)接口和片上存儲來說都是不可能實現(xiàn)的，因此按投影方向進行并行運算更適合硬件實現(xiàn).具體流程為:假設像元向量在各投影方向上的投影的初始最大和最小值均為零向量(即兩個K 維列向量的各分量均為0，分別對應最大和最小投影值向量);所有投影向量(對應于圖4(b)中投影向量矩陣的各行)與某個像元向量(對應于圖4(b)中高光譜圖像矩陣的某列)并行進行內(nèi)積運算，得到該像元向量在各投影方向上的投影值，更新最大和最小投影值向量，并記錄最大和最小值對應的像元向量位置;各投影向量與所有像元向量投影運算結(jié)束后，統(tǒng)計各像元對應最大和最小值出現(xiàn)的次數(shù)，即為該像元的純像元指數(shù).由于投影向量中的分量取值僅為1 或－1，因此可以用單個比特位的“0”和“1”表示，將隨機數(shù)算法生成的投影向量矩陣固化在硬件電路的存儲資源中，可以減少電路接口的數(shù)據(jù)讀取，有利于提高處理速度.

2.3 PPI 算法實現(xiàn)的硬件結(jié)構(gòu)

文中提出的PPI 算法實現(xiàn)的硬件結(jié)構(gòu)如圖5 所示.投影向量依次由投影向量矩陣中讀出，經(jīng)過數(shù)據(jù)同步及流水操作，格式轉(zhuǎn)化為按光譜方向排列;高光譜圖像數(shù)據(jù)經(jīng)過數(shù)據(jù)同步及流水操作后，單個像元也按光譜方向展開，像元數(shù)據(jù)和投影向量由內(nèi)積運算電路計算結(jié)果.所有投影向量順次與同一個像元向量進行運算完成投影向量并行的矩陣運算.數(shù)值比較及指數(shù)輸出采用投影向量并行的矩陣運算策略，即所有投影向量與單個像元向量依次進行并行的內(nèi)積運算，得到并更新像元向量在各投影方向上投影的最大和最小值，然后統(tǒng)計出各像元出現(xiàn)最大和最小值的次數(shù)，最終輸出各像元的純像元指數(shù).

圖5 PPI 算法實現(xiàn)的硬件結(jié)構(gòu)Fig.5 Hardware architecture of PPI algorithm implementation

3 實驗結(jié)果及實現(xiàn)性能

文中實驗所用的高光譜數(shù)據(jù)為美國國家航空航天局(NASA)發(fā)布的機載可見光/紅外成像光譜儀(AVIRIS)圖像(http:∥aviris.jpl.nasa.gov/html/aviris.freedata.html)，覆蓋光譜范圍為400～2500 nm，共224個譜段，光譜分辨率為10nm，地面分辨率為20 m.文中實驗選取Cuprite 場景中350×350 區(qū)域(見圖6(a))作為圖像數(shù)據(jù)，該區(qū)域下5 種典型礦物端元的光譜曲線(由美國地質(zhì)勘探局USGS 提供，http:∥speclab.cr.usgs.gov/spectral-lib.html)如圖6(b)所示.

采用光譜角距離(SAD)來評估實驗提取端元與USGS 光譜庫中礦物之間光譜曲線的近似程度，SAD值越小，表明兩條光譜曲線越接近.

圖6 實驗選取的高光譜圖像及其端元光譜Fig.6 Hyperspectral images selected in experiment and their endmember signatures

當K 為103、104、105和106時，Alunite 礦物的SAD 分別為0.297、0.084、0.084 和0.084.當K 從103變化到104時，SAD 從0.297 減少到0.084，說明提取的端元已明顯接近光譜庫中的參考端元;當K 繼續(xù)增加到105乃至106時，SAD 并無變化但計算量增大，表明在K=104時，端元提取精確程度和計算量之間已達到很好的平衡.K=104時采用ENVI軟件、文獻［14-15］實現(xiàn)結(jié)構(gòu)及文中提出的PPI 結(jié)構(gòu)所得到的SAD 結(jié)果見表1.由表可見，文中實現(xiàn)結(jié)構(gòu)所得到的SAD 值與文獻［15］實現(xiàn)結(jié)構(gòu)一致，且均優(yōu)于ENVI 軟件和文獻［14］實現(xiàn)結(jié)構(gòu)所得到的結(jié)果，說明文中實現(xiàn)結(jié)構(gòu)在端元提取精確性方面具有一定的優(yōu)勢，沒有因內(nèi)積運算電路中×1 或×(－1)的簡化結(jié)構(gòu)而降低精度.

表1 不同實現(xiàn)結(jié)構(gòu)獲取的端元與USGS 參考光譜庫之間的SAD 比較Table 1 Comparison of SAD between endmembers extracted by different implementations and the selected USGS reference signatures

使用Verilog HDL 語言，在型號為Xilinx Virtex-II PRO XC2VP30 的FPGA 開發(fā)板上設計并實現(xiàn)了文中提出的結(jié)構(gòu)，綜合工具為Xilinx ISE XST10.1.03.文中實現(xiàn)結(jié)構(gòu)所需要的邏輯單元僅占43%，時鐘頻率(89.174 MHz)接近90 MHz.在FPGA 開發(fā)板提供的50 MHz時鐘頻率下，由硬件結(jié)構(gòu)完成的純像元指數(shù)計算時間為9.68 s.文中實現(xiàn)結(jié)構(gòu)與文獻［14-15］實現(xiàn)結(jié)構(gòu)在對相同Cuprite 礦物區(qū)域高光譜圖像進行純像元指數(shù)計算時，所使用的Slice 資源和計算時間比較如表2所示，文中實現(xiàn)結(jié)構(gòu)所用的邏輯單元數(shù)最少、計算時間最短.

表2 不同實現(xiàn)結(jié)構(gòu)的硬件資源和計算性能對比Table 2 Comparison of hardware resource and computation performance among different implementation architectures

4 結(jié)語

文中提出了一種計算純像元指數(shù)的FPGA 實現(xiàn)結(jié)構(gòu)，旨在解決高光譜圖像中計算純像元指數(shù)運算量大的問題.該結(jié)構(gòu)采用投影向量并行的矩陣運算策略，以減少接口的數(shù)據(jù)讀取;通過譜段并行的內(nèi)積運算電路來提高投影計算的并行度，同時采用簡化的投影向量以避免硬件中難以實現(xiàn)的浮點乘法運算.實驗結(jié)果表明:該結(jié)構(gòu)的端元提取準確性優(yōu)于ENVI 軟件和其他已有結(jié)構(gòu)，對應的FPGA 實現(xiàn)方案可以滿足星載/機載現(xiàn)場實時處理的要求.隨著計算水平和硬件能力的不斷提升，在FPGA 平臺上實現(xiàn)其他復雜度更高的端元提取算法將是未來進一步研究的方向.

［1］童慶禧，張兵，鄭蘭芬.高光譜遙感:原理、技術(shù)與應用［M］.北京:高等教育出版社，2006:1-2.

［2］張兵，高連如.高光譜圖像分類與目標探測［M］.北京:科學出版社，2011:102-112.

［3］Keshava N，Mustard J F.Spectral unmixing ［J］.IEEE Signal Processing Magazine，2002，19(1):44-57.

［4］Bioucas-Dias J M，Plaza A，Dobigeon N，et al.Hyperspectral unmixing overview:geometrical，statistical，and sparse regression-based approaches［J］.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2012，5(2):354-379.

［5］Martín G，Plaza A.Spatial-spectral preprocessing prior to endmember identification and unmixing of remotely sensed hyperspectral data［J］.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2012，5(2):380-395.

［6］Boardman J W，Kruse F A，Green R O.Mapping target signatures via partial unmixing of AVIRIS data ［C］∥Proceedings of the 5th Annual JPL Airborne Earth Science Workshop.Pasadena:JPL Publication，1995:23-26.

［7］Plaza A，Martinez P，Perez R，et al.Spatial/spectral endmember extraction by multidimensional morphological operations ［J］.IEEE Transactions on Geoscience and Remote Sensing，2002，40(9):2025-2041.

［8］Rogge D M，Rivard B，Zhang J，et al.Integration of spatial-spectral information for the improved extraction of endmembers ［J］.Remote Sensing Environment，2007，110(3):287-303.

［9］Harsanyi J C，Chang C I.Hyperspectral image classification and dimensionality reduction:an orthogonal subspace projection［J］.IEEE Transactions on Geoscience and Remote Sensing，1994，32(4):779-785.

［10］Winter M E.N-FINDR:an algorithm for fast autonomous spectral endmember determination in hyperspectral data［C］∥Proceedings of SPIE Image Spectrometry V.Denver:SPIE，1999:266-277.

［11］Nascimento J M P，Bioucas-Dias J M.Vertex component analysis:a fast algorithm to unmix hyperspectral data［J］.IEEE Transactions on Geoscience and RemoteSensing，2005，43(4):898-910.

［12］Li J，Bioucas-Dias J M.Minimum volume simplex analysis:a fast algorithm to unmix hyperspectral data［C］∥Proceedings of IEEE International Geoscience and Remote Sensing Symposium.Boston:IEEE，2008:250-253.

［13］Bioucas-Dias J M.A variable splitting augmented lagragian approach to linear spectral unmixing ［C］∥Proceedings of IEEE Workshop on Hyperspectral Image Signal Process:Evolution in Remote Sensing.Grenoble:IEEE，2009:1-4.

［14］Plaza A，Plaza J，Vegas H.Improving the performance of hyperspectral image and signal processing algorithms using parallel，distributed and specialized hardwarebased systems ［J］.Journal of Signal Processing Systems，2010，61(3):293-315.

［15］Gonzalez C，Resano J，Mozos D，et al.FPGA implementation of the pixel purity index algorithm for remotely sensed hyperspectral image analysis ［J］.EURASIP Journal on Advances in Signal Processing，2010，2010:969806/1-13.

［16］Gonzalez C，Mozos D，Resano J，et al.FPGA implementation of the N-FINDR algorithm for remotely sensed hyperspectral image analysis ［J］.IEEE Transactions on Geoscience and Remote Sensing，2012，50(2):374-388.

［17］Lopez S，Horstrand P，Callico G M，et al.A novel architecture for hyperspectral endmember extraction by means of the modified vertex component analysis (MVCA)algorithm［J］.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2012，5(6):1837-1848.