基于CUBLAS和CUDA的MNF并行算法設計與優(yōu)化

周海芳+高暢+方民權

摘 要：為實現高光譜影像數據快速降維，基于nVidia 的圖像處理單元（graphic processing unit， GPU）研究最大噪聲分數變換（Maximum Noise Fraction Rotation，MNF Rotation）降維算法的并行設計與優(yōu)化，通過對加速熱點并行優(yōu)化，擇優(yōu)整合，設計并實現基于CUBLAS（CUDA Basic Linear Algebra Subprograms）庫的MNF-L（MNF-on-Library）算法和基于CPU/GPU異構系統的MNF-C（MNF-on-CUDA）算法.實驗結果顯示MNF-L算法加速11.5～60.6倍不等，MNF-C算法加速效果最好，加速46.5～92.9倍不等.研究結果表明了GPU在高光譜影像線性降維領域的巨大優(yōu)勢.

關鍵詞：圖像處理單元；GPU性能優(yōu)化；高光譜影像降維；最大噪聲分數變換；協方差矩陣計算

中圖分類號：TP391.4 文獻標志碼：A

Parallel Algorithm Design and Performance Optimization of Maximum Noise Fraction Rotation Based on CUBLAS and CUDA

ZHOU Haifang ，GAO Chang， FANG Minquan

（School of Computer， National University of Defense Technology， Changsha 410073， China ）

Abstract：To rapidly reduce the huge dimensions of hyperspectral image， this paper investigated the design and optimization of the parallel Maximum Noise Fraction （MNF） Rotation algorithm based on nVidia graphic processing units （GPUs）. In particular， a MNF-L （MNF-on-Library） algorithm based on the CUBLAS library functions and a MNF-C（MNF-on-CUDA） algorithm on the CPU/GPU heterogeneous system was presented by designing mapping schemes and parallel optimizing strategies. Experiment result shows that the MNF-L algorithm can obtain the speedups between 11.5 and 60.6， and the MNF-C algorithm can get the speedups between 46.5 and 92.9. Therefore， GPUs have a great advantage in reducing dimensions of hyperspectral images.

Key words： graphic processing unit； performance optimization for GPU； dimensionality reduction of hyperspectral image； maximum noise fraction rotation； covariance matrix calculation

高光譜遙感影像具有波段連續(xù)、光譜分辨率高的特點，能從其光譜空間中獲取豐富的地物特征信息，因其“圖譜合一”的優(yōu)勢，廣泛應用于農業(yè)、林業(yè)、軍事、環(huán)境科學、地質等各領域，具有很好的實用性和研究價值[1-2].在數據處理過程中，連續(xù)波段成像導致高光譜影像數據量龐大，且連續(xù)波段之間數據相關性強，信息冗余大，存儲處理困難，為應用和分析帶來不便，如產生維數災難、Hughes現象等[3]，因此，數據降維應運而生.

怎樣將高維空間數據映射到低維子空間，同時保持重要特征不被丟失是高光譜數據降維遵循的基本原則.高光譜影像降維主要涉及矩陣操作，如濾波、矩陣乘、協方差矩陣計算等，是典型的計算密集型和訪存密集型過程，傳統的降維過程一般采用串行方式進行，計算復雜度高，耗時長，無法滿足各領域對高光譜數據及時處理的需求[4-5].

2007年，nVidia公司發(fā)布統一計算設備架構（Compute Unified Device Architecture， CUDA），GPU并行計算開始在科學計算領域承擔重要角色[6].CPU+GPU異構系統在高性能計算領域表現突出，天河1A、泰坦等超級計算機相繼成為TOP500榜首[7].

高光譜影像并行處理已有廣泛研究，方民權[8-9]等人分別基于GPU和MIC研究了高光譜影像主成分分析降維算法，在2個GPU上獲得了128倍加速比，在3個MIC上加速133倍；Sergio[10]等人基于GPU集群實現了高光譜影像實時解混；Elmaghrbay[11]等人提出了高光譜影像端元提取的快速GPU算法；Chang[12]等人實現了一種并行模擬退火方法，并成功應用到高維遙感影像數據特征提取；Santos[13]等人基于GPU平臺實現了高光譜影像有損壓縮算法的并行加速，表明GPU在高效數據處理方面的巨大潛力.羅耀華[14]等人首次基于GPU實現了MNF并行方法研究，在數據規(guī)模為1 036 × 235 × 229時，最高取得了4倍的加速比，但該算法僅對加速熱點中的協方差矩陣運算進行了并行移植，且涉及的優(yōu)化分析較少，加速效果較差.MNF作為高光譜數據特征提取的主流算法，鮮有較全面的并行化研究，本文研究正以此為契機，對該算法的并行設計進行深入分析.

最大噪聲分數變換將原始數據中的噪聲分離，提取影像數據的主要特征，從而表征主要信息[14-15].該算法由兩層主成分變換構成，在主成分分析基礎上考慮了噪聲對圖像的影響，具有更好的效果，是目前主流的線性降維算法，實現其并行化在高光譜降維處理領域具有重要意義.

CUBLAS庫函數[16]和CUDA均可實現算法在GPU上并行加速，因其突出的加速效果而被廣泛使用.采用CUDA程序設計具有很大的靈活性，實現較為復雜，需要程序員熟練掌握GPU體系結構知識及其編程模型；CUBLAS函數庫內部整合了具體的并行和優(yōu)化細節(jié)，為使用者提供了方便的接口，但對加速算法存在限制.

本文以MNF降維算法為基礎，為實現較好的加速性能，分別基于CUBLAS庫和CUDA 進行GPU上的并行算法設計，并分析比較其性能.本文結構如下：第一部分對MNF算法進行熱點分析；第二部分提出并實現了基于CUBLAS庫的MNF GPU并行算法MNF-L（MNF-on-Library）；第三部分基于CUDA提出并實現了MNF的GPU映射和優(yōu)化算法MNF-C（MNF-on-CUDA）；第四部分通過實驗測試并分析比較MNF-L與MNF-C的并行性能；最后是總結和展望.

1 MNF算法及熱點分析

1.1 MNF算法概述

MNF算法實質是兩次層疊的主成分變換，第一層用于分離并重新調節(jié)數據中的噪聲，消除波段間的相關；第二層對噪聲白化數據進行標準主成分變換.

用X={X1，X2，X3，…，Xs}={Y1，Y2，Y3，…，YB}T表示高光譜影像數據，用S（W×H，寬W、高H）表示圖像像元數目，B表示波段數.高光譜影像實質是由B幅W×H圖像組成的三維立體模型.在實際處理中，X可用B行S列的二維矩陣表示.通過MNF變換，提取主要的m個特征波段（m

Step1. 對高光譜矩陣X濾波；

Step2. 計算濾波噪聲的協方差矩陣CN；

Step3. 對CN特征分解，

DN為降序排列的特征值，U為對應的特征向量，得到變換矩陣

Step4. 求原始高光譜X的協方差矩陣CD；

Step5. 對CD進行變換：

Step6. 對CD-adj特征值分解，

其中DD-adj為降序排列的特征值，V為所對應的特征向量；

1.2 加速熱點分析

實現高光譜影像串行MNF降維算法，對W=614，H=1 087，B=224數據降維，測試并統計各步驟占總計算時間的比例（圖1）.圖1數據顯示， Step2和Step4中協方差矩陣計算占比最大，超過80%，Step1的濾波和Step8的MNF變換耗時較為明顯，上述4個步驟是MNF算法并行設計和優(yōu)化的熱點.

2 基于CUBLAS的MNF并行算法

為簡化GPU移植編碼，nVidia引入經典BLAS（Basic Linear Algebra Subprograms）庫的CUDA實現版本——CUBLAS庫.CUBLAS庫充分利用了GPU體系結構，從底層匯編實現了BLAS庫函數的高效運算，是在GPU上進行矩陣運算的最佳選擇之一，憑借其在矩陣乘法中突出的性能表現而被廣泛應用[16-17].本節(jié)基于CUBLAS庫設計和實現協方差矩陣計算和MNF變換，并提出MNF-L算法.

2.1 基于 CUBLAS庫實現協方差矩陣計算

向量協方差計算及分解變形如公式（7）：

通過式（7）變形，將協方差計算轉換為向量求和∑X與向量內積∑XY運算，其中向量內積可轉變?yōu)檩斎敫吖庾V矩陣與其轉置矩陣乘積過程.向量和與向量內積運算可用CUBLAS庫實現：

1）調用CUBLAS庫中的cublasSasum函數求各波段向量和，循環(huán)調用B次.

2）調用cublasSsyrk（單精度）求各波段向量內積.

cublasSsyrk 接口實現的功能為C = alpha*A*AT + beta*C，令alpha= 1.0f，beta = 0.0f，輸入矩陣A為高光譜影像矩陣X.

將1）和2）的計算結果傳回CPU端，并根據公式（7），在CPU端串行計算協方差矩陣中各元素值.

2.2 采用CUBLAS庫實現MNF變換

MNF變換Z=TTX，實質為矩陣乘法運算，因此可直接調用CUBLAS庫中矩陣乘函數cublasSgemm（單精度）加速MNF變換.

接口cublasSgemm 實現的功能為C = alpha*A*B + beta*C，為了完成Z=TTX 的功能，令alpha= 1.0f，beta = 0.0f，其中輸入參數A為變換矩陣TT，參數B為高光譜影像X.

2.3 基于CUBLAS庫的MNF-L （MNF-on-Library）算法

上述兩節(jié)分別介紹了使用CUBLAS庫函數加速協方差矩陣計算和MNF變換的過程，而串行算法的另一個加速熱點——濾波，在CUBLAS庫函數中，沒有對應的函數可直接加速.因此在MNF-L算法中，該步驟采用3.2節(jié)所述CUDA并行方案.根據上述方案，整合各熱點優(yōu)化過程，提出MNF-L并行算法，如圖2.

在該算法過程中，多次調用CUBLAS庫函數，由于函數接口只定義了浮點類型的數據運算，輸入的高光譜數據存儲類型為uchar，因此必須轉換為浮點數據類型后，才能進行計算，引入數據類型轉換開銷的同時，增加了數據傳輸時間.

3 基于CUDA的MNF并行算法

利用CUDA編程模型設計GPU并行算法時，需要設計者根據算法特點對GPU線程組織層次進行設計和存儲優(yōu)化等，具有很大的靈活性，其優(yōu)化效果很大程度上取決于設計者的映射方法，訪存設計等諸多細節(jié).本節(jié)基于CUDA C，對串行算法中的3個熱點分別進行了映射方案設計和優(yōu)化，實現MNF并行算法.

3.1 協方差矩陣計算并行方案與優(yōu)化

協方差矩陣中的元素表示各向量之間的協方差，兩個向量的協方差計算公式見2.1節(jié)公式（7）.

其中參與協方差計算的是各波段所有像元組成的向量，數據量較大，可通過拆分實現并行計算；另協方差矩陣中i行j列元素的計算只需波段i和波段j的數據，因此協方差矩陣中各元素的計算相互獨立，可并行計算；這就使得協方差矩陣計算存在兩層并行.

3.1.1 GPU上協方差矩陣計算的映射方案

GPU包含grid-block-thread 3個線程組織層次，本文針對協方差矩陣運算過程提出3種GPU映射方案.

方案1 GPU中每個thread負責協方差矩陣中一個元素的計算，如圖3所示，啟動B個block，每個block啟動B個thread，實質是將grid中所有線程組織為B×B大小的矩陣，以此對應協方差矩陣中每個元素的計算任務.由于協方差矩陣的對稱性，可先根據線程id和線程塊id的大小啟動下三角（上三角）位置的線程進行計算，最后將結果矩陣對稱位置補齊即可.

方案2 GPU上每個block完成協方差矩陣中一個元素的計算（圖4）.采用二維結構組織線程塊，即啟動B×B大小的block方陣，各block分別對應協方差矩陣中相同位置元素的計算任務.同方案1類似，采用對稱補償的方法減少實際參與計算的線程塊數目.

方案3 GPU中每個grid負責協方差矩陣中一個元素的計算，即每啟動一個kernel函數完成協方差矩陣中一個元素的計算，循環(huán) （1+B）*B/2次完成協方差矩陣中所有元素的計算.由于各波段數據量（S=W×H）龐大，協方差矩陣中單個元素計算涉及高維向量加與內積運算，將其映射到grid層可減小并行粒度，且成像波段數B一般為224，使循環(huán)次數控制在可接受范圍內.

上述3種方案的本質區(qū)別在于協方差矩陣中單個元素計算的映射層次不同，后續(xù)將針對方案1展開優(yōu)化研究，根據文獻[8]中相關內容對方案2，3進行優(yōu)化，后文3.1.3節(jié)將比較3種方案的執(zhí)行效果.

3.1.2 GPU上協方差矩陣計算的共享存儲優(yōu)化

CUDA線程可以訪問不同層次存儲空間的數據，如全局，本地，共享，常量或紋理等，不同層次的存儲器訪問速度不同.共享存儲器位于片上存儲，同一個線程塊內的線程可以共享訪問，其訪存速度比全局訪存快一個數量級.

1）方案1中，同一block中的線程存在數據復用，如線程塊i中的線程重復讀取第i波段的影像數據，圖5所示，將該波段數據存儲到共享存儲器中，block內所有線程直接訪問共享存儲，最多可減少（B-1）*S次全局訪存開銷.

基于2.1中式（7）對協方差計算的改寫，在方案1映射結構上予以實現，其中矩陣（與其轉置矩陣）相乘的過程可進行共享存儲優(yōu)化，如圖6顯示，將數據分塊讀取至共享存儲單元，進行分塊迭代.建立等同于線程塊大小（thx*thx）的共享存儲區(qū)，每次讀取對應塊到共享存儲，完成分塊相乘、相加；共享存儲向右滑動，重復上述過程，直到數據末尾.

2）從全局存儲讀取數據至共享存儲時，將右乘矩陣按照轉置后的位置進行保存，使矩陣行列相乘轉變?yōu)樾行邢喑耍苊庾x取共享存儲時出現bank conflict，提高矩陣相乘的效率.

根據上述思想在原始方案1基礎上進行優(yōu)化，并采用4組數據實驗對比不同優(yōu)化算法的執(zhí)行時間，圖7所示，優(yōu)化后的執(zhí)行時間比原始版本低1～2個數量級，同時使用兩類優(yōu)化方法的加速效果更為顯著；說明本文采取的共享存儲優(yōu)化策略是非常有效的.

3.1.3 協方差矩陣計算不同方案的性能比較

測定4組高光譜數據在3.1.1節(jié)中3種映射方案下（優(yōu)化后）的協方差矩陣計算時間（不包括數據通信時間），見圖8，其中方案1采用3.1.2節(jié)兩種優(yōu)化方法（記為G-cov），在3種方案中性能最佳，說明本文選擇的設計方案和進行的優(yōu)化改進具有顯著的成效.

3.2 噪聲估計（濾波）的并行映射與優(yōu)化

濾波是對像素點進行處理以得到相應點的噪聲估計值，該步驟需要目標點及周圍8個點參與運算.圖像各點的濾波計算相互獨立，不存在相互依賴，因此圖像中所有像元點都能并行計算，且高光譜圖像不同波段也能并行計算.在GPU中，每個線程可用來計算一個像元的噪聲估計值.

3.2.1 噪聲估計（濾波）的并行映射方案

GPU上每個線程完成高光譜矩陣中一個元素的濾波任務，如圖9（記為method1），邊界噪聲事先置零，不參與映射，啟動H-2個線程塊，每個線程塊內啟動W-2個thread，將噪聲矩陣非邊界元素的濾波任務映射到相應的線程中進行計算，每次完成一個波段的映射，每個線程循環(huán)計算B次，完成所有波段元素的濾波任務.

改變method1中線程塊及線程的組織方式，均采用二維分布，將噪聲矩陣分塊映射，得到映射圖10（記為method2）.

3.2.2 濾波的GPU細粒度優(yōu)化

針對濾波過程，主要采取以下兩種優(yōu)化方法：

1）利用寄存器替換本地存儲訪問.

原執(zhí)行函數采用中間值濾波，是對包括原濾波點在內的周圍9個點進行排序，選取中間的點與原濾波點做差，在這個過程中，將使用本地存儲器，而無法使用寄存器，本地存儲實際存儲在顯存中，訪存延遲大，因此性能較差.

改進濾波函數，采用均值濾波，即將周圍所有點相加求得的平均值與原濾波點比對，此時，無需本地存儲器，而僅需要寄存器參與運算，可大大加快濾波過程.經驗證，改進的濾波函數不影響處理精度，因此可采用改進方法替代原方法.

2）利用共享存儲減少復用數據的全局訪存.

在濾波過程中，非邊界的點會被相鄰的9個點濾波使用，因此存在9次復用.利用共享存儲訪存快、線程塊內共享的屬性，可以先將線程塊所需要的數據讀入共享存儲，各線程需要時從共享存儲讀取數據，此時只需訪問一次全局存儲，減少了8次全局存儲訪問.

圖9，圖10的方案中，由于block組織方式不同，采取的共享存儲數據也有所差異.

method1中每個線程塊負責一行數據的濾波，因此，將該行與前、后兩行讀取到共享存儲區(qū)（圖11）.共享存儲大小為3*S個字節(jié)，線程塊內所有線程并行讀取一次（個別線程讀取兩次）可完成本地到共享空間的轉儲.

圖12顯示，將method2中各矩陣塊及其周圍元素數據讀取到該數據塊濾波映射的線程塊中，使該block內所有線程可以直接從共享存儲區(qū)讀取計算數據.

3.2.3 不同優(yōu)化下的濾波計算性能對比

以W=614，H=1 087，B=224的高光譜數據為輸入，對比不同方法的執(zhí)行時間，圖13前兩組數據為method1分別采用中間值濾波和均值濾波的kernel執(zhí)行時間（不包括通信時間），結果表明改進后執(zhí)行時間減少了一半以上，均值濾波中寄存器訪存加速效果顯著.

method1，method2均同時采用均值濾波和共享存儲優(yōu)化進行實現，比較兩種方法優(yōu)化后的執(zhí)行時間，如圖13后兩組數據顯示，二者執(zhí)行時間十分接近，較均值濾波優(yōu)化方法提高3～4 ms，加速比例為11.7%～15.6%，說明3.2.2節(jié)討論的兩類優(yōu)化均能使算法取得加速效果.

3.3 MNF變換的GPU映射

2.2節(jié)采用CUBLAS庫中矩陣乘函數實現MNF變換，本節(jié)針對MNF變換中相乘矩陣規(guī)模的特殊性，在傳統矩陣乘并行方案基礎上進行并行設計、改進存儲策略，以充分發(fā)揮加速潛能.

3.3.1 MNF變換的映射方案

圖14為映射圖，根據分塊計算的思想，將結果矩陣縱向分割，每個線程塊完成一個分割塊的計算，即線程塊采用一維分布，各block負責blockDim.x列數據的計算.由于S一般遠遠超過block數量，因此，每個block執(zhí)行完一塊計算后，固定跳步到下一塊計算，直到所有分塊完成.

每個線程塊采用二維結構組織線程，每個線程負責相應位置元素的計算，如果m（波段數）大于blockDim.y，各線程最多需循環(huán)計算（m+blockDim.y-1）/blockDim.y次，通常數據降維后波段數小于32，設置blockDim.y為16，使每個block縱向最多映射2次.

3.3.2 MNF變換的優(yōu)化策略和性能比較

單個線程中實質進行的是向量乘法運算，需要不斷的讀取變換矩陣各行數據與高光譜影像各列數據，直接的global 訪存延遲大，效率低，因此考慮使用共享存儲降低讀取數據帶來的延遲，而考慮到變換矩陣規(guī)模較小（一般小于32×224），沒有超出常量存儲器的存儲容量（64 kB），因此可直接將變換矩陣整體保存在常量存儲中，進一步提升訪問速度的同時，避免各線程重復讀取數據到共享存儲.

因此在該步驟中進行了如下優(yōu)化：

1）采用常量存儲器存儲變換矩陣；

2）使用共享存儲器存儲原始圖像B×blockDim.x大小的塊數據.

根據不同的優(yōu)化方式實現了4類優(yōu)化版本：

v0）全局訪存，無共享存儲優(yōu)化；

v1）僅使用共享存儲優(yōu)化；

v2）僅使用常量存儲優(yōu)化；

v3）同時使用共享存儲和常量存儲.

圖15展示了4類版本的性能對比，4組數據測試下的性能比較結果保持一致，即優(yōu)化效果v3>v1>v2>v0，v3性能最好，說明采用共享存儲和常量存儲均能取得優(yōu)化效果；隨著X數據量增加，v1較v2優(yōu)勢漸增，由于v2僅對變換矩陣使用常量存儲優(yōu)化，其全局訪存主要來自X，隨X數據量增大，v2全局訪存次數將顯著增加.

3.4 基于CUDA C的MNF-C（MNF-on-CUDA）

算法

參照前3節(jié)中關于協方差矩陣計算、濾波和MNF變換的GPU并行映射方案和優(yōu)化效果，以最優(yōu)方案替代串行算法中的相應步驟，整合提出MNF-C算法，圖16為其流程圖.

其中，第3個GPU kernel執(zhí)行期間，CPU端計算與GPU端協方差矩陣計算可同時進行，實現了Host（CPU）與device（GPU）計算重疊，充分發(fā)揮了異構系統的并行優(yōu)勢.

4 實驗結果

4.1 實驗平臺和測試數據集

本文使用CUDA C實現上述算法.實驗平臺為GPU微型超算平臺，包含2個8核的Intel（R） Xeon（R） CPU E5-2670和nVidia Kepler架構的Tesla K20c GPU，采用Intel（R） C Compiler XE 13.0.0.079和CUDA5.5工具包進行編譯.

實驗采用了4組AVIRIS高光譜數據，數據已經經過預處理，將光譜數據轉換為unsigned char類型存儲在二進制文件中.表1詳細列出了測試集的信息.

4.2 MNF-L與MNF-C性能對比

4.2.1 協方差矩陣計算性能對比

對比MNF-L中采用CUBLAS庫函數和MNF-C中程序語言實現協方差矩陣的執(zhí)行效率，采用4組數據計算總執(zhí)行時間，并分別標識運算、通信、創(chuàng)建銷毀等不同階段.如圖17所示，CUDA 實現版本（記為G-cov）同CUBLAS版本（記為L-cov）總執(zhí)行時間幾乎相同，各有兩次稍占優(yōu)勢，但整體相差細微.

高光譜數據集&實現方案

G-cov與L-cov執(zhí)行熱點主要集中在cov運算和數據傳輸，但最耗時階段形成鮮明對比，G-cov主要集中在運算階段，通信開銷較小，L-cov則與之相反.CUBLAS庫接口定義的數據類型均為浮點型，而高光譜數據初始為unsigned char，運算前需要進行類型轉換，數據量增加3倍，導致巨大通信開銷.

4.2.2 MNF變換性能對比

見表2，v3（3.3節(jié)）在MNF-C算法中采用的CUDA C實現的MNF變換，對比MNF-L中CUBLAS庫函數實現的執(zhí)行時間.

CUBLAS庫函數加速效果優(yōu)于v3版本，說明CUDA C的設計和優(yōu)化仍存在加速空間；但在實驗中發(fā)現CUBLAS首次啟動需要約200 ms的啟動開銷，大于CUDA C中的kernel的啟動開銷，如果只進行一次矩陣乘法，將難以發(fā)揮其運算優(yōu)勢.

同時隨數據規(guī)模不斷增加，CUBLAS的缺點逐漸顯現，如通信開銷、存儲限制等，4.4節(jié)將進行具體分析.

4.2.3 并行算法加速比

實驗分別測試了表1中4組數據的串行MNF降維、MNF-L和MNF-C并行降維時間（表3），計算MNF-L與MNF-C相對于串行MNF的加速比（圖18）.圖中數據顯示，基于CUBLAS庫實現的MNF-L算法可加速11.5～60.6倍不等；MNF-C算法加速效果最好，可加速46.5～92.9倍不等.

文獻[14]對MNF算法中協方差矩陣計算進行了熱點加速，使用3組高光譜數據進行測試，實驗顯示，隨高光譜圖像數據量增大，并行加速效果增加，當最大數據集為1 036 × 235 × 229時，加速比為4.58.而本文設計實現的兩種算法各有特色，較文獻[14]，采用了更全面的優(yōu)化策略和測試數據集，不僅協方差矩陣計算實現了更好的加速比，還發(fā)掘出濾波、MNF變換等步驟的加速潛能，加速效果提升了數十倍，可應用于高光譜等較大規(guī)模數據的特征提取，將大大提高執(zhí)行效率.

4.3 實驗討論

輸入不同測試數據，計算MNF-C和MNF-L算法執(zhí)行時間比值（圖19），實驗結果顯示，MNF-C算法加速效果優(yōu)于MNF-L算法，隨數據集不斷增大，兩種算法性能差距逐漸縮小.

高光譜數據集

根據算法特點，對MNF-L與MNF-C進行比較分析：

1）MNF-L中CUBLAS庫函數計算前需要先將unsigned char數據轉換為float型，引入了數據轉換和通信開銷.

2）MNF-C可顯式調整算法步驟，實現host與device的計算重疊，而MNF-L中的CUBLAS庫函數隱藏相關細節(jié)，無法實現不同設備間的協同計算.

3）在協方差矩陣計算部分，兩算法加速效果相近.雖然MNF變換中CUBLAS庫函數頗具優(yōu)勢，但該步驟所占時間比重較小，使得加速程度受限，加之CUBLAS庫函數的首次啟動開銷，一定程度上削弱了MNF-L算法的優(yōu)勢.

4）隨數據量增大，MNF-L加速效果逐漸提升，運算優(yōu)勢一定程度彌補了通信開銷，但對于更大規(guī)模的測試數據，所需存儲空間容量至少為高光譜數據的 4倍（使用雙精度接口函數達8倍），將導致顯存不足.

兩類算法各有特點，可應用于不同的加速需求領域.基于CUBLAS庫設計的并行算法，實現簡單，代碼簡潔，隱藏了相關算法的實現細節(jié)，對程序員是透明的，較容易掌握，可用于復雜和編碼量較大的算法.除了CUBLAS庫，CUDA還提供了一系列GPU加速庫函數，如針對計算機視覺和計算流體力學的cusolverDN庫、應用于化學反應動力學的cusolverSP、cusolverRF庫，應用于深度學習的cuDNN等，此類科學計算方法一般算法復雜，為簡化實現，采用GPU庫函數加速是最為直接和快速的方法.而對于GPU加速庫函數難以實現的功能或者對時效性要求較高的高端應用領域，需要獲得盡可能高的性能加速比， 那么CUDA并行設計則成為必要選擇.

總之，并行程序設計必須以運算結果正確為前提，根據算法特點、應用場景、性能要求等多方面合理選擇實現方案.

5 結 論

本文基于CPU/GPU異構系統研究了高光譜線性降維方法MNF的并行實現技術，提出了切實有效的MNF并行算法.分別針對協方差矩陣計算、濾波、MNF變換等熱點進行GPU并行加速和優(yōu)化研究，設計并實現了基于CUBLAS庫的MNF-L算法和基于CUDA的MNF-C并行降維算法.實驗結果顯示MNF-L算法加速11.5～60.6倍，MNF-C加速效果較好，獲得了46.5～92.9的加速比.最后對基于CUBLAS和CUDA實現的方法在性能、算法、應用領域等方面進行了分析討論.

針對本文應用領域——高光譜遙感數據降維，后續(xù)工作將圍繞可較好體現高維空間結構的非線性降維算法展開并行優(yōu)化的研究.

參考文獻

[1] 趙英時.遙感應用分析原理與方法[M].北京：科學出版社，2003：5-30.

ZHAO Yingshi.The principle and method of analysis of remote sensing application[M].Beijing：Science Press，2003：5-30.（In Chinese）

[2] 張良培，張立福.高光譜遙感[M].武漢：武漢大學出版社，2005：30-88.

ZHANG Liangpei，ZHANG Lifu. Hyperspectral remote sensing[M].Wuhan： Wuhan University Press，2005：30-88. （In Chinese）

[3] HUGHES G. On the mean accuracy of statistical pattern recognizers[J]. IEEE Transactions on Information Theory， 1968， 14（1）：55-63.

[4] UTO K， KOSUGI Y， SAITO G.Semi-supervised hyperspectral subspace learning based on a generalized eigenvalue problem for regression and dimensionality reduction[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2014，7（6）： 2583-2599.

[5] JIA Xiuping， JOHN A Richards. Efficient transmission and classification of hyperspectral image data[J]. IEEE Transactions on Geoscience and Remote Sensing，2003，41（5）：1129-1131.

[6] 張見明，余列祥，劉路平.基于GPU 加速的邊界面法正則積分的研究[J].湖南大學學報：自然科學版，2013，40（3）：41-45.

ZHANG Jianming，YU Liexiang，LIU Luping. Research on regular integration in boundary face method based on GPU-acceleration[J]. Journal of Hunan University：Natural Sciences，2013，40（3）：41-45. （In Chinese）

[7] TOP500. TOP500 Supercomputer Sites[EB/OL].http：//www.top500.org.

[8] 方民權， 周海芳， 申小龍. 基于GPU的高光譜PCA降維多級并行算法[J]. 東北大學學報：自然科學版， 2014，35（S1）： 238-243.

FANG Minquan， ZHOU Haifang， SHEN Xiaolong. Multilevel parallel algorithm of PCA dimensionality reduction for hyperspectral image on GPU[J]. Journal of Northeastern University：Natural Science， 2014， 35（S1）： 238-243. （In Chinese）

[9] 方民權， 張衛(wèi)民， 張理論， 等. 面向MIC的高光譜影像降維多級并行算法及性能優(yōu)化[J]. 東北大學學報：自然科學版， 2014， 35（S1）： 25-30，36.

FANG Minquan， ZHANG Weimin， ZHANG Lilun， et al. Multilevel parallel algorithms and performance optimization of hyperspectral image dimensionality reduction on MIC[J]. Journal of Northeastern University：Natural Science， 2014， 35（S1）： 25-30，36. （In Chinese）

[10]SERGIO Sánchez，RUI Ramalho，LEONEL Sousa. Antonio Plaza，real-time implementation of remotely sensed hyperspectral image unmixing on GPUs[J].Journal of Real-Time Image Processing，2015，10（3）：469-483.

[11]ELMAGHRBAY M， AMMAR R， RAJASEKARAN S. Fast GPU algorithms for endmember extraction from hyperspectral images[C]// Proceedings of the 2012 IEEE Symposium on Computers and Communications （ISCC 12）.Cappadocia，2012： 000631-000636.

[12]CHANG Y L ， CHEN K S ， HUANG B， et al. A parallel simulated annealing approach to band selection for high-dimensional remote sensing images[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2011，4（3）：579-590.

[13]SANTOS L， MAGLIE， VITULLI R，et al. Highly-parallel GPU architecture for lossy hyperspectral image compression[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing， 2013，6（2）：670-681.

[14]羅耀華，郭科，趙仕波.基于GPU 的高光譜遙感MNF并行方法研究[J]. 四川師范大學學報：自然科學版， 2013，36（3）：476-479.

LUO Yaohua，GUO Ke，ZHAO Shibo.Minimum noise fraction of hyperspectral remote sensing in parallel computing based on GPU[J]. Journal of Sichuan Normal University：Natural Science，2013， 36（3）：476-479. （In Chinese）

[15]LIU Xiang，ZHANG Bing，GAO Lianru，et al. A maximum noise fraction transform with improved noise estimation for hyperspectral images[J].Science in China， Series F： Information Sciences，2009（9）：1578-1587.

[16]MANUEL Carcenac.From tile algorithm to stripe algorithm： a CUBLAS-based parallel implementation on GPUs of Gauss method for the resolution of extremely large dense linear systems stored on an array of solid state devices[J].The Journal of Supercomputing，2014，68（1）：365-413.

[17]ZHANG Bo，YANG Xiang，YANG Fei，et al. The CUBLAS and CULA based GPU acceleration of adaptive finite element framework for bioluminescence tomography[J].Optics Express，2010， 18（19）：20201-20214.