高光譜圖像貝葉斯分類算法GPU并行優化研究

趙海娜，吳遠峰，高建威，張兵

(1.中國科學院 遙感與數字地球研究所 數字地球重點實驗室，北京 100094；2.中國科學院大學，北京 100049)

1 引 言

當今遙感技術正朝著更高空間分辨率、更高時間分辨率和更高光譜分辨率方向發展，這給海量遙感圖像分類帶來了新的挑戰，尤其在災害監測、戰場目標探測以及智能衛星系統等時效性要求較高的應用領域的挑戰尤為突出[1-2]。

貝葉斯分類方法具有堅實的數學理論基礎以及綜合先驗信息和數據樣本信息的能力，是遙感圖像分類的實用方法[3-4]。但基于CPU的貝葉斯串行分類算法在處理較大圖像時非常耗時。

傳統遙感圖像分類加速處理多采用由CPU組成的分布式集群系統的處理模式[5-7]。但是集群系統價格昂貴、體積大、靈活性差，很難應用到戰場或星上等要求低重、高靈活性等場景。

近幾年來，國內外許多研究將GPU用在遙感影像融合、混合像元分解和實時目標探測等領域，對算法進行了加速處理，獲得高效的并行計算能力[8-10]。圖形處理器(Graphic Processing Unit，GPU)具有運算密集、高度并行、體積小和性價比高等特點，為解決數據密集型的計算提供了有利條件[11]。本文針對貝葉斯分類算法面臨大數據量圖像處理速度問題，提出了基于GPU的貝葉斯并行分類算法，并從數據傳輸和核函數設計兩方面對并行算法進行了優化設計。實驗結果表明，基于GPU的貝葉斯并行分類算法(以下簡稱并行算法)與基于CPU的串行算法(以下簡稱串行算法)相比，在保證分類精度的同時，分類效率明顯提高，能夠獲得25倍～54倍的加速比。

2 貝葉斯分類算法

2.1 算法流程

定義遙感圖像為X，band、sample和line分別代表波段數、列數和行數。x={a1,a2,…，am},(m=band)為其中的一個像元，ai為x的一個特征屬性；Y={y1,y2,…，yn}代表類別，n為類別數，D={d1,d2,…，dn}為訓練樣本，其中di代表每一類的訓練樣本的個數。貝葉斯分類方法[3-4]是一種基于統計模型的分類方法，并假設各對象之間是條件獨立的[12]，其原理是基于某地物的先驗概率，利用貝葉斯公式計算出其后驗概率，即該地物屬于某一類的概率，選擇具有最大后驗概率的類作為該地物所屬的類。計算步驟如下：

①計算訓練樣本中每個類別的先驗概率

(1)

②計算訓練樣本的均值μyi和標準差σyi；

③根據樣本均值、標準差計算每個類別條件下各個特征屬性條件概率

(2)

p(x|yi)p(yi)

=P(a1|yi)P(a2|yi)…P(am|yi)P(yi)

(3)

⑤用分類器進行分類，選擇概率最大的類作為X的所屬類別，對每一像元進行分類。

2.2 貝葉斯分類算法計算復雜度分析

2.1節中步驟①和步驟②計算了訓練樣本各類別的先驗概率以及均值和標準差，計算過程簡單并且計算量很小，因此在整個算法中這兩步運算的復雜度可以忽略不計。實驗中將先驗概率、訓練樣本的均值標準差作為已知。

表1是步驟③～步驟⑤對于處理分類類別為n，圖像尺寸為line×sample×band的圖像的時間復雜度。從以上的分析可以看出，數據量對算法復雜度的影響非常大，數據量的增加會嚴重影響算法的計算效率。

表1 各處理過程計算時間復雜度

3 基于GPU的分類算法并行優化

在CUDA架構下，程序計算分為在CPU上執行的host端和在GPU上執行的device端。Host端負責算法流程的控制以及簡單的計算；device端負責圖像的數據級并行計算[13]。根據上一節對貝葉斯算法復雜度的分析，基于CUDA對算法進行并行設計時，將計算流程拆分成兩部分：①數據讀寫、GPU初始化、函數調用等流程控制邏輯由CPU執行；②計算復雜度高的步驟③～步驟⑤由GPU執行，算法流程如圖1所示。

圖1 基于CUDA的貝葉斯分類并行算法流程

圖1中，虛線框代表該部分計算在GPU上執行，分配到GPU上并行計算的部分由處理流程中的Kernel函數(核函數)完成，每個Kernel函數對應CUDA計算模型上的一個GRID，GRID包含多個Block計算單元，Block由多個Thread線程具體執行計算，將圖像像元映射到計算線程從而實現大規模的數據級并行計算。

在計算前，需要將圖像數據從主機內存拷貝到GPU的設備內存，計算完成后，再將分類后的數據從GPU的設備內存拷貝回主機內存。在這個過程中CPU與GPU的計算核心都處于空閑狀態。為了隱藏CPU和GPU間的數據傳輸時間，充分利用GPU的計算資源，本文通過CUDA創建流管理貝葉斯算法中的條件概率計算、后驗概率計算以及分類器分類過程。流內的操作是按順序執行的，而流與流之間則是亂序執行。實現一個流中數據傳輸與另一個流中核函數的執行同時進行，即GPU能夠持續地進行計算，而不必等待數據。以3個流為例，每個流處理一塊圖像數據，異步優化效果如圖2所示，不使用異步執行時貝葉斯并行模型的計算時間為(執行時間+傳輸時間)；而使用流和異步后，時間降為(執行時間+傳輸時間/3)。

圖2 異步優化示意圖

同時，由于Kernel的啟動會占一定的時間開銷，為了提高并行算法的計算效率，本文盡可能將多個計算過程合并到一個Kernel函數中執行。貝葉斯算法中的特征屬性條件概率計算、后驗概率計算以及分類器分類都是針對整個遙感圖像進行的，圖像像元與計算線程的映射相同，將以上三個過程合并到一個Kernel中進行計算，從而可以避免不必要的時間開銷。

4 實 驗

4.1 數據

所用的高光譜圖像是PHI航空成像光譜儀獲取的日本某地區的精細農業數據，共80個波段，原始像幅為570×350。圖3為原始圖像假彩色合成影像以及地面樣本圖。經主成分分析進行特征提取，選取變換后包含信息量最大的6個波段作為實驗數據。根據地面調查數據選取了中國白菜、日本白菜、蘿卜、生菜、牧場、薯、森林、無植被覆蓋區、地膜覆蓋區等9類地物作為分類的訓練樣本。另外，為了分析不同數據量對算法加速比的影響，將原始圖像空間維重采樣，生成5幅不同空間分辨率的模擬影像，像幅分別為256×256，512×512，1024×1024，2048×2048，4096×4096。

圖3 PHI日本某地區的精細農業數據

4.2 計算環境

本文的實驗測試平臺參數為，中央處理器CPU是2.8GHz、12個計算核心的Intel(R) Xeon(R) X5660；圖形處理器是NVIDIA Quadro 5000 GPU，2.5GB全局存儲器，352個計算單元。操作系統是Windows 7 64位。

貝葉斯分類算法分別基于串行和并行計算環境實現，前者采用Microsoft Visual Studio 2008 C/C++編譯器實現，后者采用GPU的統一計算架構CUDA4.0 C實現。

4.3 結果與分析

本文從分類精度和計算時間兩個方面對提出的并行算法的性能進行分析，并與對應的串行程序進行比較。首先，對串行、并行優化算法的分類精度進行比較。從多組實驗結果得出對于相同的實驗數據串并行算法得到的分類精度是一致的，圖4給出了分類結果。

圖4 PHI數據分類結果

為了定量評價分類精度，本文根據地面調查數據分別在影像的不同位置選取了4個～6個均勻區域作為驗證樣本，各類地物的驗證樣本分別包含約200個～1000個像元。通過混淆矩陣計算了貝葉斯算法的分類精度，如表2所示。

表2 貝葉斯算法分類精度

接下來對兩種算法的分類效率進行對比分析。實驗時間統計結果如圖5、圖6所示。其中并行計算時間是指包括GPU分配顯存時間、數據拷貝時間以及核函數執行時間在內的總時間；由于采用了流異步執行，核函數執行時間包括了從設備端向主機端的數據拷貝時間。加速比是指串行算法與并行算法計算時間的比值。

由圖5可以看出，當數據量比較小時，核函數執行時間與并行計算時間基本相同，但隨著數據量的增大，并行計算時間與核函數執行時間之間的差值不斷增大，原因是GPU顯存分配、主機端設備端數據拷貝時間隨著數據量的增大而不斷增加。

圖5 核函數執行與并行計算時間在不同數據量下的比較

圖6為串并行計算時間以及加速比隨數據量變化圖。從圖6(a)中的串并行計算時間曲線可以看出，串并行計算時間隨數據量的增加均呈線性增長。對于傳統CPU環境下，縱坐標計算時間增量與橫坐標數據量增量之比(秒每兆，s/M)超過1∶5，在這種情況下，如果數據量比較小，其分類的處理時間還能夠容忍，如果數據量增大，則計算時間隨之成倍增加；對于GPU環境下，縱坐標計算時間的增量與橫坐標數據量的增量之比小于1∶200，即由數據量急劇增加而引起的計算時間增長十分緩慢。從圖6(b)中的加速比曲線來看，基于GPU并行計算的加速比隨數據量的增加呈增大趨勢，當數據量達到一定量時，計算效率不再有較大的提升，加速比基本趨于穩定，說明此時GPU線程達到滿載，即使增加數據量，加速比變化也不明顯；從多次實驗結果可知，本文提出的并行算法能夠獲得25倍～54倍的加速比。

更進一步，將本文提出的并行算法與ENVI軟件中的最大似然算法進行對比。在處理4096×4096大小的圖像時，本文的并行算法計算時間是1.78s，而ENVI軟件中最大似然算法的計算時間大約是34.2s，本文的并行計算優化方法比ENVI軟件同類功能的計算效率快19倍。

從以上實驗結果分析可以得出結論如下：(1)本文提出的貝葉斯并行算法與串行算法的分類結果一致；(2)在保證分類結果正確的同時，并行算法的計算效率要遠遠高于串行算法，并能獲得25倍～54倍的計算加速比；(3)串行計算時間增量與數據量增量之比大于1∶5，并行計算時間增量與數據量增量之比小于1∶200，數據量增加時，串行計算時間增長的更快；(4)加速比隨數據量的增大呈增大趨勢，當數據量達到一定時，加速比趨于穩定，此時GPU線程達到滿載。

圖6 串并行計算時間以及加速比隨數據量的變化

5 結束語

本文基于高光譜圖像貝葉斯分類算法的復雜度分析，給出了GPU并行優化計算方法，并且采用PHI航空成像光譜儀獲取了日本某地區的精細農業數據，在NVIDIA GPU上進行了實驗驗證。結果表明，該并行分類算法在保證分類精度的同時能大大提高算法的計算效率，獲得25倍～54倍的計算加速比。研究表明，基于GPU實現遙感圖像并行分類算法是提高遙感圖像分類處理計算效率的有效途徑，可以為遙感大數據處理提供技術保障。同時，如何利用多GPU系統進一步提高算法的計算效率也是下一步研究工作的重點。

參考文獻：

[1] ANTONIO P.Special issue on architectures and techniques for real-time processing of remotely sensed images[J].Journal of Real-Time Image Processing，2009，4(3)：191-193.

[2] 張兵.智能遙感衛星系統[J].遙感學報，2011，15(3)：415-431.

[3] DOMINGOS P，PAZZANI M.On the optimality of the simple bayesian classifier under zero-one loss[J].Machine Learning，1997，29(2-3)：103-130.

[4] 彭興媛.樸素貝葉斯分類改進算法的研究[D].重慶：重慶大學，2012.

[5] 劉曉云，康一梅，齊同軍，等.遙感圖像波譜角并行分類算法[J].計算機科學，2009，36(9)：267-270.

[6] LUO W，ZHANG B，JIA X.New improvements in parallel implementation of N-FINDR algorithm[J].Geoscience and Remote Sensing，IEEE Transactions on，2012，50(10)：3648-3659.

[7] PLAZA A J.Parallel techniques for information extraction from hyperspectral imagery using heterogeneous networks of workstations[J].Journal of Parallel and Distributed Computing，2008，68(1)：93-111.

[8] 盧俊，張保明，黃薇，等.基于 GPU 的遙感影像數據融合 IHS 變換算法[J].計算機工程，2009，35(7)：261-263.

[10] TARABALKA Y，HAAVARDSHOLM T V，KSEN I，et al.Real-time anomaly detection in hyperspectral images using multivariate normal mixture models and GPU processing[J].Journal of Real-Time Image Processing，2009，4(3)：287-300.

[11] 袁濤，馬艷，劉定生.GPU 在遙感圖像處理中的應用綜述[J].遙感信息，2012，27(6)：110-117.

[12] KIM S，PARK S，WON D.Proxy signatures，revisited[J].Information and Communications Security，1997：223-232.

[13] NVIDIA.NVIDIA CUDA compute unified device architecture—programming guide[EB/OL].http：//developer.nvidia.com/cuda.