江子特 趙遼英 鄒佳林

(杭州電子科技大學計算機學院 浙江 杭州 310018)

基于Fan模型非負矩陣分解的光譜解混并行計算

江子特 趙遼英 鄒佳林

(杭州電子科技大學計算機學院 浙江 杭州 310018)

基于FAN模型的廣義非負矩陣分解是一種非純像元假設下有效的高光譜圖像非線性光譜解混算法。針對基于FAN模型的廣義非負矩陣分解算法的快速實現問題，基于CUDA編程模型與存儲器模型設計并行優(yōu)化，對優(yōu)化后算法的串行與并行部分進行任務分配與線程映射，設計合理的核函數實現各關鍵步驟。通過真實高光譜數據的光譜解混實驗，結果表明CUDA并行優(yōu)化后的算法相比串行算法，能達到較高的加速比，驗證了其有效性。

廣義非負矩陣分解 非線性混合模型 并行計算 高光譜圖像

0 引 言

高光譜遙感技術以其圖譜合一的特點受到國內外的廣泛關注，被廣泛應用于環(huán)境監(jiān)測、礦產勘察、軍事偵察等領域?；旌舷裨钠毡榇嬖谑怯绊懜吖庾V圖像實際應用效果的一個重要因素，為此，混合像元分解一直是高光譜圖像處理研究的熱點。光譜混合從本質上分為線性混合和非線性混合兩種[1]，雙線性混合模型是非線性光譜解混中用得較多的一類非線性混合模型，包括Fan等提出的Fan模型[2]和Halimi等[3]提出的廣義雙線性混合模型等。自高光譜遙感出現后，各種線性和非線性的光譜解混技術不斷被提出[4-5]，其中非負矩陣分解NMF(Nonlinear matrix fraction)[6]是基于線性混合模型的非監(jiān)督光譜解混技術的典型代表。由于沒有考慮不同端元間的非線性相互作用，NMF不能有效地對非線性光譜混合數據進行光譜解混。針對該問題，文獻[7]基于Fan模型提出了稱為Fan_NMF的廣義NMF算法。

高光譜圖像的海量數據使得各種光譜解混算法不能很好地滿足高光譜圖像的實時應用。由于能夠為算法提供接近計算機集群的高計算能力， GPU/CUDA架構已經成為高性能并行高光譜圖像處理的一個重要選擇[8-9]。羅耀華等[10]提出基于GPU并行實現了高光譜遙感最小噪聲分離變換算法，使算法運行速度至少提高了4倍；文獻[11]提出了一種基于CUDA的高光譜端元投影向量法的實現方法，使其運算效率提高了十幾倍；葉舜[12]研究了基于GPU的稀疏性高光譜圖像混合像元分解的并行優(yōu)化算法，設計了CPU+GPU異構并行計算方法，在Telsa C2050平臺上進行實驗并取得了14倍的加速比?？梢?，基于GPU的并行計算可以顯著提高高光譜遙感圖像的處理效率。

分析Fan_NMF算法過程發(fā)現其具有很強的并行性，非常適合GPU加速。本文研究基于CUDA的Fan_NMF并行計算，設計了具體流程和關鍵步驟的并行計算方法，最后通過實驗進行了驗證。

1 Fan_NMF 算法

假設R=[r1,…,rN]∈L×N表示高光譜圖像數據矩陣，N為像素個數，L為圖像的光譜波段數，則rk的Fan光譜混合模型可以形式化表示為:

(1)

式中P為端元個數；M=[m1,m2,…,mP]∈L×P為端元光譜矩陣，sk=[sk1,sk2,…,skP]T∈P中的每個元素ski為該像元中端元mi對應的豐度；ε表示可能的誤差或噪聲項；⊙表示Hadamard積。

豐度向量sk受到如下兩個約束條件的限制：

1≥ski≥0i=1,2,…,P

(2)

(3)

記S=[s1,s2,…,sN]∈P×N，令Mb=[m1⊙m2,…,mP-1⊙mP]∈L×P(P-1)/2，Sb=[s1⊙s2,…,sP-1⊙sP]T∈P(P-1)/2×N，基于Fan_NMF的光譜解混問題可以描述為：給定非負矩陣R,尋找非負矩陣M和S，使得R≈MS+MbSb。

(4)

通過乘性迭代求解目標函數式(4)，得到式(5)和式(6)表示的更新規(guī)則：

(5)

(6)

3) 根據式(4)計算目標函數f(M(t),S(t))。

4) 根據式(5)更新M(t+1)，根據式(6)更新S(t+1)。

5)t=t+1，如果t

2 基于CUDA的Fan_NMF并行優(yōu)化

Fan_NMF算法主要運算量在迭代過程，每次迭代要先對端元矩陣和豐度矩陣進行Hardmard積運算，求目標函數，進行端元和豐度矩陣的更新。因此，整個迭代過程由GPU處理，采用CPU+GPU的協同工作達到加速的目的。

并行設計流程圖如圖1所示。主機端主要負責原始高光譜圖像數據的讀取操作、初始端元矩陣和豐度矩陣的讀取以及端元矩陣及豐度矩陣的主機端存儲，設備端則完成算法迭代，最終將計算結果拷貝至主機端。為了簡化設計，只以迭代次數滿足設定值為迭代結束條件。

圖1 基于CUDA架構的FAN_NMF算法結構流程圖

下面介紹圖1所示的設備端7個步驟的具體設計過程。

Step1端元矩陣M的偽逆計算，利用公式(MTM)-1MT，調用CUBLAS庫函數cublasSgemm實現MTM，調用CULA函數culaDeviceSgetrf對MTM進行LU分解，再調用culaDeviceSgetri函數得到(MTM)-1，最后再調用cublasSgemm函數得M的偽逆。各函數啟動的線程數量以及每個塊的維度依據端元矩陣維度的大小而定。

Step2將高光譜數據及端元矩陣的偽逆作為輸入，通過cublasSgemm求解初始豐度，再設計Kernel函數創(chuàng)建N個線程，每個線程計算一個像素的歸一化豐度，將結果保存至全局存儲器。

Step3以初始端元矩陣為輸入，啟動點乘Kernel函數，求端元矩陣的Hadamard積。

Step4以豐度矩陣為輸入，啟動點乘Kernel函數，求豐度矩陣的Hadamard積。

Step5端元矩陣的更新，每個Kenrel啟動L×P個線程，具體過程如表1所示。

表1 端元矩陣更新步驟

Step6重復Step3的計算過程，重新計算端元矩陣的Hadamard積。

Step7豐度矩陣更新，每個Kenrel啟動L×N個線程，具體過程如表2所示。

表2 豐度矩陣更新步驟

按照Step4-Step7重復迭代，每次迭代結束，迭代次數加1，如果迭代次數大于設定的最大迭代次數，則終止迭代，得到最終的端元矩陣以及豐度矩陣。

3 實驗結果及分析

實驗數據1：HYDICE采集到的城市高光譜數據集URBAN，如圖2(a)所示。圖像大小為307×307像素，其包含210個波段。為了提升算法的性能，去除低信噪比和水蒸氣吸收波段(1～4，76，87，101～111，136～153和198～210)后剩余162個波段用于本實驗，該區(qū)域主要包含六種地物：道路、金屬、瀝青、草地、樹木和屋頂。

實驗數據2：農作物圖，如圖2(b)所示，大小為360×550，共124個波段。

圖2 真實高光譜圖像(第30波段)

實驗平臺是Windows 7 32位Intel core i5-3580p和NVIDIA GTX 660Ti及CUDA 5.5開發(fā)包。實驗中用FNSGA方法[14]提取初始端元，全約束最小二乘法[15]求初始豐度?；贑UDA的Fan_NMF并行計算迭代1 000次后，將CUDA處理的實驗結果寫入txt文件中，通過Matlab顯示結果，分別得到圖3和圖4所示的光譜曲線和圖5和圖6的豐度圖。從圖5和圖6，可以看出所求的各種地物的分布與圖2的目視效果基本一致。

圖3 URBAN圖端元光譜曲線

圖4 農作物圖端元光譜曲線

圖5 URBAN圖的豐度圖

圖6 農作物圖的豐度圖

最大迭代次數設置為1 000次時，兩幅真實高光譜圖像基于CPU和GPU實現Fan_NMF解混的時間比較如表3所示，其中，所有GPU耗時不考慮I/O數據傳輸以及圖像數據的初始化時間，CPU耗時不包括讀取圖像數據時間。設高光譜圖像數據大小為height×width，設置gridDim為[(width+15)/16]× [(height+15)/16]，blockDim(線程塊維度)為16，即線程塊的大小為256。從表3中可以看出，基于CUDA的并行Fan_NMF算法相比于基于串行CPU實現的Fan_NMF算法性能大幅的提升?；贑UDA的并行Fan_NMF算法處理農作物crops高光譜圖像數據時的加速比達到了54.72倍，對URBAN圖像進行并行處理加速比為80.96倍。

表3 GPU和CPU運行時間比較

表4給出了基于CUDA并行計算中各部分耗時。從表4中可以看出，在Fan_NMF算法并行實現中，端元及豐度更新操作因為其計算量大，所以占用了相當

部分的并行處理時間，端元擴展以及豐度擴展在計算中耗費的處理時間相對較少，各個Kernel之間進行數據傳輸也耗費了相當部分的時間，對不同步驟的耗時分析可進一步找出性能優(yōu)化點。

表4 CUDA實現Fan_NMF算法中各部分耗時 s

4 結 語

本文實現了基于CUDA的Fan_NMF算法并行計算，除了添加虛擬波段操作在主機端Host完成外，其余計算都在設備端完成，減少了主機端與設備端之間的數據交互。在Windows 7 32位的Intel core i5-3580p，NVIDIA GeForce GTX 660Ti實驗平臺對兩幅真實高光譜圖像進行Fan_NMF光譜解混，分別取得了54.72倍和80.96倍的加速比。這充分證明了GPU并行計算針對數據量大的高光譜數據解混的優(yōu)勢。本文方法在實現上仍然存在一定的加速空間，例如可以使用共享存儲器來提高數據存取效率從而提高加速比。因此，后續(xù)工作中還需進一步對程序進行存儲器方面的優(yōu)化。

[1] 童慶禧,張兵,鄭蘭芬.高光譜遙感——原理、技術與應用[M].北京:高等教育出版社,2006:38-44.

[2] Fan W,Hu B,Miller J,et al.Comparative study between a new nonlinear model and common linear model for analysing laboratory simulated forest hyperspectral data[J].International Journal of Remote Sensing,2009,30(11):2951-2962.

[3] Halimi A,Altmann Y,Dobigeon N,et al.Nonlinear unmixing of hyperspectral images using a generalized bilinear model[J].IEEE Transactions on Geoscience and Remote Sensing,2011,49(11):4153-4162.

[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 & Remote Sensing,2012,5(2):354-379.

[5] Dobigeon N,Tourneret J Y,Richard C,et al.Nonlinear Unmixing of Hyperspectral Images:Models and Algorithms[J].IEEE Signal Processing Magazine,2013,31(1):82-94.

[6] 賈森.非監(jiān)督的高光譜圖像解混技術研究[D].杭州:浙江大學,2007.

[7] Eches O,Guillaume M.A Bilinear-Bilinear Nonnegative Matrix Factorization Method for Hyperspectral Unmixing[J].IEEE Geoscience & Remote Sensing Letters,2013,11(4):778-782.

[8] 張舒.GPU高性能運算之CUDA[M].中國水利水電出版社,2009:13-21.

[9] 趙穎輝,蔣從鋒.遙感影像的高性能并行處理技術研究[J].計算機技術與發(fā)展,2014(7):201-205.

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

[11] 羅耀華,郭科,趙仕波.基于線性光譜模型的混合像元分解高性能計算研究[J].物探化探計算技術,2013(3):344-348.

[12] 葉舜.基于GPU的高光譜圖像混合像元分解并行優(yōu)化研究[D].南京理工大學,2014.

[13] Bertsekas D P.Constrained optimization and Lagrange multiplier methods[M].New York:Academic,1982.

[14] 王麗姣.基于單形體體積增長的高光譜圖像端元提取及快速實現[D].浙江大學,2015.

[15] Heinz D C,Chein-I-Chang.Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[J].IEEE Transactions on Geoscience & Remote Sensing,2002,39(3):529-545.

SPECTRALUNMIXINGPARALLELCOMPUTINGOFNONNEGATIVEMATRIXDECOMPOSITIONBASEDONFANMODEL

Jiang Zite Zhao Liaoying Zou Jialin

(SchoolofComputerScienceandTechnology,HangzhouDianziUniversity,Hangzhou310018,Zhejiang,China)

The generalized non-negative matrix factorization algorithm based on Fan model is an effective nonlinear hyperspectral unmixing algorithm without the pure pixel assumption. First, the parallel optimization algorithm based on CUDA programming model and memory model was designed for the fast implementation of the generalized nonnegative matrix factorization algorithm based on FAN model. Then, task assignment and thread mapping were performed for the serial and parallel parts of the optimized algorithm and reasonable kernel functions were designed for the implementation of the key steps. The spectral unmixing experiment on real hyperspectral data shows that the CUDA parallel optimization algorithm can achieve higher speedup than the serial algorithm, which verifies the validity of the proposed algorithm.

Generalized non-negative matrix factorization Nonlinear mixture model Parallel computing Hyperspectral image

2017-01-17。國家自然科學基金項目(61571170)；浙江省自然科學基金項目(LZ14F030004)。江子特，本科生，主研領域：計算機應用與軟件。趙遼英，教授。鄒佳林，碩士生。

TP751.1

A

10.3969/j.issn.1000-386x.2017.12.018