郭杰 李云松 弋方 吳憲云

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

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

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

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

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

1 PPI 算法

1.1 線性光譜混合模型的幾何學(xué)描述

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

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

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

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

1.2 PPI 算法描述

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

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

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

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

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

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

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

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

2.2 矩陣運(yùn)算策略

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

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

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

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

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

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

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

圖6 實(shí)驗(yàn)選取的高光譜圖像及其端元光譜Fig.6 Hyperspectral images selected in experiment and their endmember signatures

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

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

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

4 結(jié)語

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

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

［2］張兵，高連如.高光譜圖像分類與目標(biāo)探測［M］.北京:科學(xué)出版社，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.