面向高光譜影像分類的高性能計算及存儲優化

高 偉，李維良，林 妍

GAO Wei,LI Weiliang,LIN Yan

中國地質大學（武漢）信息工程學院，武漢430074

Department of Information Engineering,China University of Geosciences,Wuhan 430074,China

高光譜遙感影像分類是遙感定量信息提取的關鍵步驟之一，由于高光譜影像波段數目多，數據量往往較大，導致分類算法內存消耗大，計算時間較長。對此可以從計算機實現角度來進行算法改進，如進行并行化處理，以提高效率，滿足實際應用的需要。

針對遙感影像分類技術的并行化研究，國內外已經出現一些研究成果。蔣艷凰對遙感圖像高精度并行監督分類技術進行了研究，首先提出幾種學習算法用于提高監督分類的預測精度；然后采用并行處理技術，提高遙感圖像監督分類的處理速度[1]；楊靖宇等重點研究了基于GPU 的多光譜遙感影像分類的流程，通過實驗驗證了該技術方法的有效性[2]；Antonio Plaza 等利用異構網絡環境下的工作站實現了一種綜合影像光譜信息與空間信息的形態學分類算法，取得了較為理想的加速效果[3]；Dai Lijun 等提出了一種基于最大概率與加權平均的分類器，并采用計算機集群對影像分類過程進行了并行化研究[4]；Ujjwal Maulik 等提出了一種基于點對稱的遙感影像分類技術，并利用計算機集群對算法進行了并行化實現，取得了線性加速比[5]。

現有的研究多數是基于計算機集群、工作站來開展的，具有成本較高，部署困難等缺陷，不適合桌面遙感軟件的應用推廣；少數基于GPU 方式的研究主要是從流程的角度來論證該并行架構對提高算法效率的有效性，對基于GPU 的算法優化策略，尤其是存儲器優化策略，尚且涉及得不深入或不全面。

本文以高光譜遙感影像分類為例，研究GPU 架構下的桌面遙感影像處理系統高性能計算的設計思路、實現方法與關鍵技術，其中對算法存儲優化策略進行重點研究。在這個過程中，詳細分析了影像數據的組織形式及其對GPU 存儲器訪問效率的影響；改進通常所使用的2D blocks 映射方式為1D blocks 映射方式以實現全局顯存的最佳訪問形式；利用常量存儲器進行標準樣本光譜向量的存儲以減少GPU 對低速顯存的訪問頻率；利用共享存儲器進行影像數據的存儲以提升GPU 運算速度。

1 高光譜影像分類方法概述

經過幾十年的發展，高光譜遙感影像數據的處理與分析取得了長足的進步，形成了一系列面向高光譜影像特點的分類算法。這些算法可以總結為兩種思路：一種是基于地物物性的分類方法，主要是利用反映地物性質的光譜曲線來識別；另一種思路是基于圖像數據的分類方法，主要是利用數據的統計特征來建立分類模型。其中，基于物性，即基于地物的光譜反射或發射曲線的分類識別方法最具特色，這種方法的主要特點是利用光譜庫中已知地物的光譜數據，采用光譜匹配技術來識別圖像中的地物覆蓋類型。光譜匹配技術是通過對兩個光譜曲線進行比較與特定運算，來求解它們之間的相似度或者差異度。光譜匹配過程可以是全波長范圍內的比較，也可以是基于整波形特征的比較，或是利用感興趣波段的光譜，進行部分波長范圍的光譜組合參量的匹配。比較典型的匹配算法有：二值編碼匹配法、光譜波形匹配法以及光譜角填圖法等[6]。本文以光譜波形匹配法以及光譜角填圖法為例加以介紹。

1.1 高光譜影像光譜波形匹配法分類

光譜波形匹配包括兩種匹配方式，一種是將樣本光譜的全部或某一部分進行光譜曲線的特征函數擬合，通過計算像光譜與樣本光譜特征函數之間的擬合度來計算像元光譜隸屬于某一樣本的概率。常用的特征函數有植被倒高斯模型、光譜吸收谷函數模型等。另一種方式是直接計算影像中各像元光譜矢量與樣本光譜矢量之間的線性相似度來確定各像元對應地物所應歸屬的類別[7-8]。對于同一地物類型具有很高的線性相似度，而非同一地物類型則具有較低的線性相似度。Clark 等提出了一種波段擬合算法，對于n個波段的像元光譜與樣本光譜來說，兩者的光譜擬合程度可以用F表示：

式中，F越大，表示像元光譜與樣本光譜之間的擬合程度越高。本章以此式為代表進行采用CUDA 技術進行高光譜影像分類技術的實驗介紹。

1.2 高光譜影像光譜角填圖（SAM）法分類

光譜角填圖將像元N個波段的光譜響應作為N維空間矢量，通過計算它與標準端元組分光譜矢量之間的廣義夾角來表征其匹配程度：夾角越小，說明越相似，從而確定每個像元所應歸屬的類別。設兩個n維空間矢量分別為T與R、為矢量R和矢量T的第i分量，則兩矢量廣義夾角可以通過公式（2）求解[6]：

式中，值越小，T與R的相似性越大。

對科萊恩新型催化劑的活性的分析中，測試了維持64%的轉化率所需的入口溫度（T64），且與對標催化劑（見圖3）進行比對。如果所需溫度降低，則表明催化劑活性更高。

2 高光譜影像分類算法的并行處理

2.1 CUDA 編程模型與存儲結構

CUDA是NVIDIA公司2007年提出的一種通用GPU計算模型。在CUDA 結構中，CPU 端被稱為主機端（host），GPU 端被稱為設備端（device），采用SIMT（single instruction multiple thread）模式執行程序。GPU 上并行執行的程序稱為內核函數（kernel）[9-10]。kernel 在執行時創建很多線程（thread），若干線程會被組織成線程塊（block），而若干block 再組成網格（grid）。每個thread 按照指定block ID 與thread ID 來唯一確定。

CUDA 的計算資源主要來源于GPU 內置的高度線程化的多核流處理器（Streaning Multiprocessor，SM）陣列，每個SM 又包含多個流處理器（Streaming Processor，SP）。通常，一個線程塊會被分配到一個SM，塊內的每個線程被映射到一個SP 上。程序通過大量線程來獲得并行性，這種并行主要有細粒度并行與粗粒度并行兩個級別：細粒度級并行是指block 內線程之間的并行，粗粒度級并行是指block 之間的并行。

圖1 計算任務分解過程

CUDA 的存儲器按線程訪問權限可以劃分為3 個層次。第1 層為單個線程的私有存儲器，包括寄存器文件（register files）、局部存儲器（local memory），第2 層為線程塊中所有線程可以共享的存儲器（shared memory），第3 層為網格內所有線程可以訪問的最外層存儲器，包括全局存儲器（global memory）、常量存儲器（const memory）和紋理存儲器（texture memory）。

2.2 高光譜影像分類算法的并行化

基于CUDA 的高光譜遙感影像分類，須首先確定網格中線程塊的數量以及每個線程塊中執行線程的數量。

線程塊數量的確定由待求解問題的規模來確定。在CUDA中，一個block必須被分配到一個SM中，但一個SM同一時刻可以有多個活動線程塊（active block）等待執行。一般來說，為便于理解，常常把影像網格劃分成維度大小為(IMG_WIDTH+X_BLOCK_DIM–1)/X_BLOCK_DIM，(IMG_HEIGHT-Y_BLOCK_DIM–1)/Y_BLOCK_DIM的2D block集合[11-12]，其中，IMG_WIDTH、IMG_HEIGHT分別為影像的寬度與高度；X_BLOCK_DIM、Y_BLOCK_DIM分別為每個block 內X方向上的線程數與Y方向上的線程數。

在確定了網格中的線程塊數量之后，每個block 待處理的子圖像塊就得以確定。此時，需要確定block 中每個執行線程的計算任務。在SAM 法高光譜影像分類過程中，子圖像塊中所有像素都要與各標準樣本進行光譜特征值的比較與計算，尋找SAM 值最小的標準樣本，這是分類過程中計算最密集的地方，也是線程計算任務劃分的關鍵。如圖1 所示，讓block 中的各線程T0、T1、T2、…、Tn去對應處理子圖塊內像素點p0、p1、p2、…、pn的光譜特征值，將影像像素映射到并行處理線程，線程的計算任務得到較好劃分。整個計算任務分解過程如圖1 所示。

完成計算任務分配后，就可以為內核函數提供影像數據并啟動內核函數。基于CUDA 的高光譜遙感影像分類計算過程如圖2 所示，可以總結為：將影像數據從內存拷貝到顯存；確定每個block 的維度與線程數量，建立block 與影像子圖塊之間的映射關系，并對block 中的線程進行計算任務分配，確定每個線程將要處理的像素；調用kernel 函數在GPU 上進行并行分類計算；kernel函數運行完畢后，將運算結果從顯存拷貝到內存。

圖2 基于CUDA 的并行分類流程圖

3 并行分類算法的存儲優化策略

CUDA 實質上將多個線程捆綁執行，每個block 被劃分成warp，warp 內線程執行同一條指令[13]，half-warp是存儲操作的基本單位[14]。本文重點從共享存儲器、全局存儲器、常量存儲器三個方面來探討高光譜影影像分類過程中GPU 的存儲優化策略，在存儲優化過程中，需重點關注half-warp 中的線程。

3.1 共享存儲器數據存儲與訪問優化

GPU 片內沒有對顯卡的全局存儲器進行緩存，導致全局存儲器的訪問具有較大延遲。因此，想要進一步提高分類算法的運算效率，必須最小化GPU 對全局存儲器的訪問。并行分類算法kernel 函數中，可以將block所對應的子圖塊定義為共享型。在kernel 函數執行之初，將指定的子圖塊從全局顯存載入到流多處理器的共享型存儲器中，以供后續分類過程的處理。由于每一個block 在圖像內獲取的子圖像塊不同，以塊索引號和線程索引號為判斷分支，可以將子圖塊載入共享存儲器這個過程作為kernel 函數的一部分在GPU 上并行進行。根據block 的索引號，每個block 找到與其對應的子圖像塊，再根據block 中thread 的索引號，每個thread 將其對應的子圖塊像素從全局存儲器中取出，傳入到共享存儲器中以供后續計算，如圖3 所示。

圖3 將子圖塊載入共享存儲器示意圖

共享存儲器被劃分為大小相等，能被同時訪問的存儲器模塊，稱為bank。Bank 的數據組織方式是：每個bank 的寬度為4 Byte，相鄰的4 Byte 數據被組織在相鄰的bank 上，每個bank 在每個時鐘周期可以提供4 Byte的帶寬，不同的bank 可以互不干擾地同時工作。如果half-warp 中的線程請求訪問的地址位于同一個bank中，存儲器模塊在同一時間無法響應多個請求，所有這些請求必須串行完成，降低了存儲器訪問效率。本文中，每個half-warp 中的線程對應處理一個像素的數據，為了使thread ID 相鄰的線程在同一時鐘周期內訪問相鄰的bank，共享存儲器中的子圖塊按照BSQ 方式進行存儲。

3.2 全局存儲器數據存儲與訪問優化

全局存儲器數據存儲與訪問優化的重點是獲取存儲器合并訪問條件，盡管不同計算能力的CUDA 設備對合并訪問條件的限制并不相同，其表現為隨著設備計算能力的提高，對合并訪問條件的限制相應地降低，但是，盡量遵循嚴格的合并訪問條件，預期可以較大地提升存儲器訪問效率。

當half-warp 中所有線程在執行同一條指令時，訪問經特定字長對齊的存儲段中的連續單元時，全局存儲器獲得最佳訪問形式，這種情況下，硬件把所有訪問請求結合成一個對全局存儲器連續單元的合并訪問[13]。一般而言，存儲段需要按照線程訪問字長的half-warp對齊[14]，本文采用float 類型對影像數據在顯存中進行存儲，因此，存儲段需要按照64 Byte進行對齊。

影像數據通常有BSQ、BIL、BIP三種存儲方式。halfwarp中每個線程對應處理一個像素，當影像數據按照BIL與BIP 兩種方式進行存儲時，half-warp 中相鄰線程所訪問的全局顯存地址不連續，無法滿足合并訪問條件，顯存訪問效率會受到影響，而當影像數據按照BSQ 方式進行存儲時，這種情況會得到一定程度的緩解，如圖4所示。

采用運行API cudaMalloc（）函數分配顯存，能夠保證其首地址按照256 Byte 對齊[15]。當影像各子圖像塊與2D blocks 建立映射后，首block 塊中首行線程滿足合并訪問條件，當選擇合適的block 塊大小（half-warp 的倍數）時，網格中首行block 塊中的首行線程滿足合并訪問條件，然而，網格中其他線程，如圖5 線程塊C 中的線程，可能由于待訪問的存儲段沒有對齊到指定地址而無法獲得最佳訪問形式。一種優化方法是將2D blocks 映射改進成1D blocks 映射，blocks 中的線程依照block ID 與thread ID 依次對應處理一個像素，當block 中線程數量為half-warp 的倍數時，可以獲得最佳訪問形式。

圖4 各數據存儲方式對應的顯存訪問示意圖

圖5 2D blocks映射所對應的全局顯存訪問示意圖

3.3 常量存儲器數據存儲優化

CUDA 提供了常量存儲器，常量存儲器位于全局存儲器中，但采用了緩存提高了訪問效率[15]。在影像分類過程中，各個標準樣本的光譜特征向量也會十分頻繁地與影像各個像素的光譜特征向量進行比較與計算。標準樣本的光譜特征向量數據在并行分類開始之前從主機內存拷貝到GPU 全局存儲器，且在整個分類過程中維持不變。在GPU 分類運算開始之前，將標準樣本的光譜特征向量存放在常量存儲器中，可以在一定程度上減少對全局存儲器的訪問次數。

4 實驗結果與分析

根據以上研究思路，對基于CUDA 的高光譜影像分類存儲優化策略進行了實驗。本文所使用的硬件實驗環境為：NVIDIA GEFORCE GT 425M 顯卡，該顯卡所使用的GPU 內置兩個SM，每個SM 擁有48 個SP，共96 個SP，顯卡的顯存容量為1 GB；軟件實驗環境為：Windows 7 操作系統及VC++2005 編譯器。

4.1 高光譜遙感影像分類結果

實驗首先對比利用CPU 與利用GPU 的高光譜影像分類結果。實驗數據采用AVIRIS 高光譜影像，原始影像波段數量為60，在本次實驗中，為加大相鄰波段之間的差異性，減少相關性，便于計算結果的觀察，采用等間隔取樣的方法選擇其中25 個波段的數據形成新的原始影像作為實驗數據，影像尺寸為350×350，影像數據類型的單精度浮點型，如圖6 所示。

圖6 分類結果實驗原始影像

高光譜遙感影像分類結果如圖7 所示，其中圖7（a）為CPU 串行計算情況之下的光譜角匹配法分類結果，圖7（b）為基于CUDA 的光譜角匹配法分類結果，圖7（c）為CPU 串行計算情況之下的SAM 法分類結果，圖7（d）為基于CUDA 的SAM 法分類結果。

4.2 全局存儲器、共享存儲器、常量存儲器數據存儲與訪問優化

接下來，進行全局存儲器、共享存儲優化、常量存儲的優化實驗。實驗所用的高光譜影像與4.1 節中所使用的數據相同。為滿足實驗所用CUDA 設備共享存儲器容量大小的限制，block 尺寸不宜過大，在本次實驗中設置為16×8。實驗分4 組進行，第1 組實驗不進行存儲優化，第2 組為全局存儲器優化實驗，采用3.2 節中所示的略仍能在一定程度上減少GPU 對低速全局存儲器的過度訪問，算法計算效率仍有一定程度的提高。1D blocks 與子圖塊進行映射，每個block 中包含16×8個線程，第3 組與第4 組分別為共享存儲器優化實驗與常量存儲器優化實驗，均采用2D blocks 與子圖塊進行映射。

圖7 分類結果

表1 全局存儲器、共享存儲器、常量存儲器數據存儲與訪問優化實驗結果 ms

表2 實驗影像信息

從表1 第2 列與第3 列的對照實驗結果可見，由于影像尺寸（350×350）的制約，在進行2D blocks 映射時，線程在對影像第2 列以及以后各列數據進行訪問時，無法獲得存儲器的最佳訪問形式，對程序運行效率造成影響，這種情形在采用1D blocks 映射優化后得到一定程度的改善。從表1 第4 列與第2 列的對照實驗結果可見，共享存儲器優化有效地減少了GPU 對全局存儲器的訪問頻率，較大程度地提升了程序計算的效率。從表1 第4 列與第5 列以及第5 列與第2 列的對照實驗結果可見，在僅使用常量存儲器優化策略的情況之下，算法優化率不及僅使用共享存儲器優化策略的情況，但是，相比不采取任何存儲優化措施的情況，常量存儲器優化策

4.3 高性能計算運算效率實驗

接下來，采用多幅不同尺寸與數據量的高光譜遙感影像對比GPU 高性能計算情況之下與CPU 串行計算情況下的分類算法用時，以證明GPU 高性能計算對高光譜遙感影像分類算法效率提升的有效性。表2 顯示了這些影像的尺寸與數據量信息，各影像的波段信息與4.1 節中所述的相同。鑒于顯存容量的限制，在將影像從內存拷貝到顯存的過程中，采用分塊拷貝的策略，每次拷貝整副影像中適宜尺寸的影像塊進行并行處理，這里影像塊大小采用256×256 像素；同樣為滿足實驗所用CUDA 設備共享存儲器容量大小的限制，block尺寸不宜過大，這里設置為16×8。實驗結果如表3 與表4 所示。

從表3 與表4 所示的CPU 計算時間與GPU 運行時間中的有效計算時間對比可見，針對各種尺寸與數據量的遙感影像，無論是光譜波形匹配分類算法還是SAM分類算法，在CUDA 架構下的計算時間相對于CPU 而言都有大幅度的縮減。盡管由于內存與顯存之間數據傳輸的影響，使得在GPU 運行總時間中有效計算時間所占的比例降低，在一定程度上影響了CUDA 架構下分類算法的運算效率，但是，由于GPU 使得有效計算時間大幅減少，且內存與顯存之間數據傳輸時間并不長，從表3 以及表4 最后一列所示加速比可知，針對不同尺寸與數據量的影像，CUDA 架構下的高光譜遙感影像分類算法的運行效率都得到了顯著的提升。

表3 光譜波形匹配分類算法效率對照

表4 SAM 分類算法效率對照

表5 分類精度對照

4.4 高性能計算運算精度實驗

接下來評價高性能計算的精度。在下面的實驗中，計算GPU 高性能計算下與CPU 串行計算下結果影像對應位置像素值之間的差異，并將差異率從大到小劃分若干個區間，計算落在各個區間的像素數量及所占的比例。實驗采用4.3 節中所示的影像作為數據源，實驗結果如表5 所示。

表5 第2 大列顯示了GPU 與CPU 計算結果在各差異范圍內的像素數量及所占比例，從中可以看出，無論SAM 算法還是光譜波形匹配算法，GPU 高性能計算與CPU 串行計算的結果中除極少數像素的像素值差異在10%之上，超過99.99%的像素的差異率在10-8之下。基于CUDA 的高光譜影像高性能計算在大幅度提升計算效率的同時維護了計算結果的準確性。

5 結束語

本文以高光譜遙感影像分類算法為例，介紹利用多核CUDA 技術實現該算法并行化的設計思路與實現方法，重點對包括共享存儲器、常量存儲器、全局存儲器在內一系列GPU 存儲優化策略進行研究。實驗取得良好效果。為探索利用GPU 技術，充分挖掘現有單機計算機系統的計算資源，提升桌面遙感影像處理系統的運行效率，提供了一條具有借鑒意義的途徑。

當然，本文在以高光譜影像監督分類為例探討利用CUDA 技術進行桌面遙感影像處理系統高性能計算時，重點從CUDA 設備各類存儲器的利用研究算法優化策略，對于其他策略如：指令選擇、指令展開、最小化主機與設備間數據傳輸次數等則沒有涉及。為進一步提升算法運行效率，這些都有待后續研究與探討。

[1] 蔣艷凰.遙感影像高精度并行監督分類技術研究[D].長沙：國防科學技術大學，2004.

[2] 楊靖宇，張永生，董廣軍.基于GPU 的遙感影SAM 分類算法并行化研究[J].測繪科學，2010，35（3）：9-11.

[3] Plaza A，Plaza J，Valencia D.AMEEPAR：parallel morphological algorithm for hyperspectral image classification on heterogeneous.networks of workstations[C]//Lecture Notes in Computer Science，2006，3993：24-31.

[4] Dai Lijun，Liu Chuang.Research on remote sensing image of land cover classification based on multiple classifier combination[J].Wuhan University Journal of Nature Science，2011，16（4）：363-368.

[5] Maulik U，Sarkar A.Efficient parallel algorithm for pixel classification in remote sensing imagery[J].GeoInformatica，2012，16（2）：391-407.

[6] 童慶禧，張兵，鄭蘭芬.高光譜遙感——原理、技術與應用[M].北京：高等教育出版社，2006.

[7] 楊哲海.高光譜影像分類若干關鍵技術的研究[D].鄭州：解放軍信息工程大學，2006.

[8] 楊國鵬.基于核方法的高光譜影像分類與特征提取[D].鄭州：解放軍信息工程大學，2007.

[9] 沈玉琳.通用GPU 計算技術在高性能計算平臺上的應用研究[D].蘭州：蘭州大學，2012.

[10] 趙進.基于GPU 的遙感圖像并行處理算法及其優化技術研究[D].長沙：國防科學技術大學，2011.

[11] 肖漢，張祖勛.基于GPGPU 的并行影像匹配算法[J].測繪學報，2010，39（1）：46-51.

[12] 羅耀華.高性能計算在高光譜遙感數據處理中的應用研究[D].成都：成都理工大學，2013.

[13] KIRK D B，W HWU Wen-mei.大規模并行處理器編程實戰[M].陳曙暉，熊淑華，譯.北京：清華大學出版社，2010.

[14] 張舒，褚艷利.高性能運算之CUDA[M].北京：中國水利水電出版社，2009.

[15] NVIDIA.NVIDIA CUDA C programming guide[EB/OL].（2012）.http：//developer.download.nvidia.com/compute/Dev-Zone/docs/html/C/doc/CUDA_C_Programming_Guide.pdf.