李 宇，高教波，孟合民，張 磊，張茗璇

（西安應(yīng)用光學(xué)研究所，陜西 西安710065）

引言

成像光譜技術(shù)由于其不僅具有二維空間信息，還能獲取目標(biāo)的光譜維信息，對(duì)目標(biāo)識(shí)別有著極其重要的作用，因此近年來(lái)受到了各研究機(jī)構(gòu)的重視，文獻(xiàn)［1］研究得出了針對(duì)高光譜數(shù)據(jù)的快速重建技術(shù)，CUDA并行架構(gòu)非常適合于處理高光譜原始數(shù)據(jù)［1－2］，針對(duì)成像光譜技術(shù)而發(fā)展起來(lái)的系統(tǒng)分析、設(shè)備裝置研發(fā)、數(shù)據(jù)處理算法也不斷跟進(jìn)。成像光譜技術(shù)已經(jīng)在很多領(lǐng)域取得了較為成熟的應(yīng)用，如區(qū)域制圖、環(huán)境監(jiān)測(cè)、農(nóng)業(yè)遙感、公共安全等［3］。

時(shí)空調(diào)制型干涉光譜成像系統(tǒng)，相對(duì)于傳統(tǒng)的空間調(diào)制型和光柵分光型光譜成像系統(tǒng)具有高通量的優(yōu)勢(shì)，適合對(duì)微弱光譜輻射信號(hào)的探測(cè)，在氣體檢測(cè)方面具有一定的優(yōu)勢(shì)。針對(duì)大部分氣體在長(zhǎng)波紅外波段都有其獨(dú)特的吸收峰的特征，我們?cè)O(shè)計(jì)并研制了時(shí)空調(diào)制型長(zhǎng)波紅外干涉成像光譜系統(tǒng)，用于對(duì)氣體成分的檢測(cè)。由于氣體檢測(cè)屬于實(shí)時(shí)性要求較高的應(yīng)用場(chǎng)合，這就需要從推掃成像到最終的數(shù)據(jù)解算結(jié)果要有很高的處理效率，但是受干涉成像光譜的原理限制，采用傳統(tǒng)的CPU串行處理方法完全不能達(dá)到性能的要求。為此，利用CUDA平臺(tái)實(shí)現(xiàn)干涉光譜數(shù)據(jù)并行處理的研究，并取得了一定的成效，實(shí)驗(yàn)結(jié)果表明，采用GPU作為協(xié)處理器實(shí)現(xiàn)并行處理，比僅使用CPU作串行計(jì)算在計(jì)算效率上有了5到20倍的提升。

1 時(shí)空混合調(diào)制型干涉成像原理

自主研制的時(shí)空混合調(diào)制型長(zhǎng)波紅外干涉成像光譜儀基于分體式Sagnac干涉儀，樣機(jī)工作波段為8μm～12μm。

1.1 系統(tǒng)裝置

分體式Sagnac干涉儀采用三角共光路結(jié)構(gòu)，依靠反射鏡的非對(duì)稱(chēng)布置及分光鏡對(duì)光線的偏折作用產(chǎn)生橫向剪切量，分光干涉原理如圖1所示［4］。圖中，光學(xué)系統(tǒng)由分光干涉儀、長(zhǎng)波紅外成像物鏡及熱像儀組成，分光干涉原理如下：

1）遠(yuǎn)處物點(diǎn)的平行光入射分光干涉儀，分為2束，2束光以一定橫向剪切量平行出射；

2）2束光經(jīng)過(guò)成像物鏡在焦平面FPA上合束并干涉；

3）干涉光譜成像系統(tǒng)相對(duì)目標(biāo)運(yùn)動(dòng)，物點(diǎn)平行光進(jìn)入干涉系統(tǒng)角度不同，獲取不同光程差的干涉信息；

4）干涉條紋經(jīng)傅里葉逆變換，反演出目標(biāo)輻射的光譜信息［5］。

圖1 Sagnac干涉光譜成像儀Fig.1 Sagnac interferometer

1.2 干涉圖獲取模式

時(shí)空混合調(diào)制型干涉光譜儀獲取干涉圖的模式與其他干涉光譜儀有所區(qū)別，由于Sagnac采用前置光學(xué)系統(tǒng)，每次進(jìn)入視場(chǎng)的是當(dāng)前區(qū)域范圍內(nèi)的所有目標(biāo)像元，而且每個(gè)目標(biāo)像元在當(dāng)前幀上只有一個(gè)干涉數(shù)據(jù)點(diǎn)信息，要獲取單個(gè)目標(biāo)像元的完整干涉信息，就需要按照一定的規(guī)律，從不同的CCD幀中取得該像元的所有干涉信息后組合起來(lái)。被測(cè)目標(biāo)像元相對(duì)于儀器的入射角不斷變化，產(chǎn)生不斷變化的光程差，從而形成一系列干涉圖像。所以這就要求儀器和被測(cè)目標(biāo)之間必須是相對(duì)勻速運(yùn)動(dòng)，同時(shí)滿足采樣定理，方可獲得被測(cè)目標(biāo)均勻的干涉數(shù)據(jù)［6］。

2 CUDA架構(gòu)

圖形處理單元（graphic process unit，GPU）由于其眾核設(shè)計(jì)理念，使得它在浮點(diǎn)運(yùn)算能力上遠(yuǎn)超過(guò)CPU，因此非常適合數(shù)據(jù)密集型的復(fù)雜計(jì)算任務(wù)。到2009年，GPU浮點(diǎn)運(yùn)算的執(zhí)行速度可以達(dá)到1×1012次／s（1teraflops），而此時(shí)的多核CPU執(zhí)行速度僅是1×1011次／s，二者在峰值浮點(diǎn)運(yùn)算的吞吐量方面比率達(dá)到了10∶1，而這還僅僅是芯片中執(zhí)行資源可以支持的原始速度，不一定是應(yīng)用程序在運(yùn)行時(shí)能夠達(dá)到的最高速度［7］。

統(tǒng)一計(jì)算設(shè)備架構(gòu)（compute unified device architecture，CUDA）的出現(xiàn)，綜合了CPU的邏輯控制能力和GPU的數(shù)據(jù)密集型并行計(jì)算能力，并且使用對(duì)C語(yǔ)言的簡(jiǎn)單擴(kuò)展CUDA C作為編程語(yǔ)言，使得并行程序開(kāi)發(fā)變得相對(duì)容易［8］。在本研究中，我們采用CUDA架構(gòu)實(shí)現(xiàn)干涉數(shù)據(jù)的反演處理來(lái)解決處理效率的問(wèn)題。

2.1 CUDA并行原理

在CUDA中使用的是單指令多線程（SIMT）計(jì)算模型，即同一條指令可以處理多個(gè)數(shù)據(jù)對(duì)象。這相當(dāng)于使得用戶對(duì)程序有著線程級(jí)的控制，通過(guò)跳轉(zhuǎn)和分支，每個(gè)線程可以有獨(dú)立的、不同的行為，并且在GPU中線程是輕量級(jí)的，即其創(chuàng)建、切換和銷(xiāo)毀所耗費(fèi)的時(shí)間可以忽略不計(jì)。CUDA使用一種分層的編程模型來(lái)組織線程：用戶根據(jù)需要定制一定數(shù)量的線程準(zhǔn)備投入計(jì)算，并可以定制線程和數(shù)據(jù)的映射關(guān)系。

一個(gè)支持CUDA的GPU中至少包含一個(gè)流多處理器（SM），每個(gè)SM包含若干個(gè)標(biāo)量處理器（SP）。在實(shí)際運(yùn)行中，控制器將一個(gè)線程組分配給一個(gè)SM后，SM中的SP協(xié)同工作，并行處理所有的線程。

2.2 算法設(shè)計(jì)

時(shí)空混合調(diào)制型光譜成像系統(tǒng)，光譜特征信息以干涉條紋的形式疊加在景物所成的像上。從雙邊干涉圖像到光譜信息的光譜反演過(guò)程需要經(jīng)過(guò)以下步驟：

1）提取重排

同一物點(diǎn)以不同光程差干涉成像于各幀圖像的不同位置，因此需先將該物點(diǎn)的干涉強(qiáng)度信息從各幀圖像中提取出來(lái)，按照光程差順序重新排列。

2）基線校正

由于視場(chǎng)漸暈的影響，圖像邊緣的干涉條紋強(qiáng)度有所降低，為正確反演光譜，需對(duì)干涉條紋進(jìn)行基線校正處理以修正趨勢(shì)項(xiàng)及減去直流分量。

3）切趾

理想傅里葉變換的積分上下限是光程差從負(fù)無(wú)窮大到正無(wú)窮大，而實(shí)際儀器中的干涉圖函數(shù)在光程差達(dá)到最大值處截止，這使得儀器譜線函數(shù)變寬，反演光譜峰值周?chē)霈F(xiàn)較強(qiáng)的“旁瓣”，掩蓋了波峰旁邊的弱光譜信息，造成分辨光譜的困難。這個(gè)“旁瓣”出現(xiàn)的物理根源在于最大光程差附近干涉圖尖銳的不連續(xù)性，消除或抑制這些“旁瓣”的基本方法就是用一漸變權(quán)重函數(shù)來(lái)乘干涉圖，即通常所說(shuō)的“切趾”的方法進(jìn)行數(shù)學(xué)濾波，以緩和這種不連續(xù)性［9］。對(duì)條紋進(jìn)行切趾處理能有效降低旁瓣，但會(huì)使光譜展寬，不同的切趾函數(shù)有不同的特性［10］。

4）傅里葉逆變換

以提取重排、基線校正和切趾處理后的干涉條紋作傅里葉逆變換，設(shè)定采樣頻率，確定反演光譜的位置。利用傳統(tǒng)的CPU處理模式［11］，其計(jì)算流程如圖2所示。

圖2 傳統(tǒng)CPU處理流程Fig.2 Processing procedure using traditional CPU

而對(duì)于每個(gè)空間像元來(lái)說(shuō)，其所需要的處理都是上述的相同步驟，因此各像元之間獨(dú)立、可并行的特點(diǎn)非常適合使用CUDA實(shí)現(xiàn)，在本研究中基于CUDA架構(gòu)的干涉光譜反演并行算法對(duì)獲取的干涉圖進(jìn)行處理，處理流程如圖3所示。

圖3 采用CUDA架構(gòu)的處理流程Fig.3 Procedure using CUDA

3 實(shí)驗(yàn)結(jié)果及分析

文中采用的實(shí)驗(yàn)數(shù)據(jù)來(lái)自于自研干涉成像光譜儀采集的干涉數(shù)據(jù)，探測(cè)器視場(chǎng)大小為640×480像元，采用12bit量化，數(shù)據(jù)處理采用三角窗切趾和256點(diǎn)傅里葉變換。針對(duì)同樣的算法分別實(shí)現(xiàn)了CPU的串行計(jì)算版本和CUDA架構(gòu)的并行計(jì)算版本，并對(duì)比兩者的計(jì)算效率以及計(jì)算結(jié)果的正確性。

本實(shí)驗(yàn)采用的運(yùn)行環(huán)境是：

操作系統(tǒng)：Windows 7旗艦版64位SP1。

處理器：Intel Xeon X5690＠3.47GHz六核。

內(nèi)存：8GB

顯卡：Nvidia Quadro 4000，顯存大小為4GB，CUDA計(jì)算能力為2.0，CUDA驅(qū)動(dòng)版本為5.0。

開(kāi)發(fā)編譯環(huán)境：Microsoft Visual Studio 2008 SP1＋NVIDIA CUDA編譯器NVCC。

實(shí)驗(yàn)針對(duì)黑體經(jīng)單色儀出射10μm單色光進(jìn)行推掃，數(shù)據(jù)處理過(guò)程中各步驟的結(jié)果見(jiàn)圖4～圖7。圖4所示是對(duì)采集到的數(shù)據(jù)進(jìn)行重排以后，在黑體表面選取的一個(gè)像素點(diǎn)的完整干涉序列，我們將所有像素點(diǎn)干涉序列組成的干涉立方體數(shù)據(jù)送至顯存，由GPU進(jìn)行并行計(jì)算的任務(wù)。圖5所示是為了校正視場(chǎng)漸暈而導(dǎo)致的邊緣能量畸變，將圖4中彎曲的譜線拉平。圖6是去直流及采用三角窗切趾得到的曲線。圖7是CUDA針對(duì)圖4中選取的干涉序列最終得到的反演結(jié)果，我們沒(méi)有采用CUDA自帶的cufft庫(kù)，而是在核函數(shù)中實(shí)現(xiàn)了頻域兩點(diǎn)IFFT算法，計(jì)算結(jié)果只進(jìn)行了波長(zhǎng)定標(biāo)，而未作輻射定標(biāo)。圖8是在CUDA中利用OpenGL渲染出的可視化結(jié)果，可以看到單色光的空間分布范圍很明顯。

圖4 10μm單色光原始采集干涉序列Fig.4 Original data sequence of 10μm homogeneous light

圖5 CUDA消除趨勢(shì)項(xiàng)Fig.5 Eliminating tendency using CUDA

圖6 CUDA三角切趾Fig.6 Triangle apodization using CUDA

圖7 10μm單色光CUDA反演結(jié)果Fig.7 Inversion result of 10μm homogeneous light using CUDA

圖8 10μm單色光反演結(jié)果顯示Fig.8 Inversion result visualization of 10μm homogeneous light

表1是采用CUDA和僅采用CPU對(duì)相同輸入數(shù)據(jù)和算法流程的計(jì)算結(jié)果對(duì)比（輸入數(shù)據(jù)量：515.625MB（640×480×880×2B）），從表中可以看出，對(duì)于同樣規(guī)模的數(shù)據(jù)量，采用CUDA計(jì)算只需要5.3s，而采用傳統(tǒng)的CPU計(jì)算則需要33s，說(shuō)明了CUDA完全能夠使用并行計(jì)算代替原本的串行計(jì)算，提高計(jì)算效率。

表1 CUDA和CPU計(jì)算時(shí)間對(duì)比Table 1 Calculation time performance contrast between CUDA and CPU

在近期的研究中，改為采用1024點(diǎn)傅里葉逆變換，使用CUDA計(jì)算320×480像元大小的光譜立方體僅需要13s，而采用CPU計(jì)算時(shí)間則達(dá)到了接近300s，加速比超過(guò)了20倍。表明對(duì)于復(fù)雜度高的算法來(lái)說(shuō)，CUDA的并行處理模式相對(duì)于CPU來(lái)說(shuō)具有很高的效率。

從圖9和圖10中還可以看出，采用CUDA計(jì)算和僅采用CPU計(jì)算，分別得到計(jì)算結(jié)果，其誤差基本在±0.005以?xún)?nèi)，所以采用CUDA進(jìn)行計(jì)算并不會(huì)對(duì)計(jì)算的精度有很大的影響。

通過(guò)分析工具Visual Profile的分析結(jié)果可以得知，應(yīng)用程序中仍然存在不足之處，影響了加速比的顯著提高，有幾方面的原因：一是采用的數(shù)據(jù)結(jié)構(gòu)不利于顯存的合并訪問(wèn)，全局內(nèi)存存／取和浮點(diǎn)運(yùn)算次數(shù)之比不理想；二是目前采用的并行方式所導(dǎo)致的，在每個(gè)線程中還是需要處理一維向量，而同一時(shí)刻并行執(zhí)行的線程數(shù)量是有限制的。

圖9 CUDA和CPU反演結(jié)果對(duì)比Fig.9 Inversion result contrast between CUDA and CPU

圖10 CUDA和CPU反演計(jì)算誤差Fig.10 Calculating deviation between CUDA and CPU

4 結(jié)束語(yǔ)

通過(guò)對(duì)時(shí)空混合調(diào)制型成像光譜儀干涉圖獲取模式的分析，結(jié)合CUDA平臺(tái)，實(shí)現(xiàn)了干涉光譜的并行反演算法，并通過(guò)實(shí)驗(yàn)驗(yàn)證了采用CUDA進(jìn)行并行計(jì)算，將干涉光譜數(shù)據(jù)的反演過(guò)程加速了5到20倍。

［1］ 楊智雄，廖寧放，呂航，等.基于CUDA的傅里葉變換成像光譜儀數(shù)據(jù)重建［J］.光學(xué)技術(shù)，2013，39（2）：128-132.YANG Zhi-xiong，LIAO Ning-fang，Lü Hang，et.al.Fast data reconstruction method of Fourier transform imaging spectrometer based on CUDA［J］.Optical Technique，2013，39（2）：128-132.（in Chinese with an English abstract）

［2］ 王建宇，舒嶸，劉銀年，等.成像光譜技術(shù)導(dǎo)論［M］.北京：科學(xué)出版社，2011.WANG Jian-yu，SHU Rong，LIU Yin-nian，et al.Introduction of imaging spectrum tecniques［M］.Beijing：Science Press，2011.（in Chinese）

［3］ SCHOWENGERDT R A.遙感圖像處理模型與方法［M］.北京：電子工業(yè)出版社，2010.SCHOWENGGERDT R A.Remote sensing：modeling methods for image processing［M］.Beijing：Publishing House of Electronics Industry，2010.（in Chinese）

［4］ 張茗璇，高教波，孟合民，等.基于傅里葉變 換光譜技術(shù)的Zoom-FFT 算法研究 ［J］.應(yīng) 用 光 學(xué)，2013，34（3）：452-456.ZHANG Ming-xuan，GAO Jiao-bo，MENG He-ming，et al.Zoom-FFT based on Fourier transform spectroscopy［J］.Journal of Applied Optics，2013，34（3）：452-456.（in Chinese with an English abstract）

［5］ 孟合民，高教波，肖相國(guó)，等.紅外高通量干涉成像光譜儀的設(shè)計(jì)與驗(yàn)證［J］.紅外與激光工程，2012，41（11）：2903-2908.MENG He-min，GAO Jiao-bo，XIAO Xiang－guo，et al.Design and validation of infrared interferential imaging spectrometer with high flux［J］.Infrared and Laser Engineering，2012，41（11）：2903-2908.（in Chinese with an English abstract）

［6］ 嚴(yán)新革，張淳民，趙葆常.時(shí)空混合調(diào)制型偏振干涉成像光譜儀干涉圖獲取模式研究［J］.物理學(xué)報(bào)，2010，59（5）：3123-3129.YAN Xin-Ge，ZHANG Chun-Min，ZHAO Baochang.Research on the mode of obtaining interferograms based on the temporally and spatially mixd modulated polarization interference imaging spectrometer［J］.Acta Physica Sinica，2010，59（5）：3123-3129.（in Chinese with an English abstract）

［7］ KIRK D B，HU Wen-mei.大規(guī)模并行處理器編程實(shí)戰(zhàn)［M］.陳曙輝，熊淑華，譯.北京：清華大學(xué)出版社，2010.KIRK D B，HU Wen-mei.Programming Massively Parallel Processors［M］.translated by CHEN Shuhui，XIONG Shu-hua.Beijing：Tsinghua University Press，2010.（in Chinese）

［8］ FARBER R.高性能CUDA應(yīng)用設(shè)計(jì)與開(kāi)發(fā)：方法與最佳實(shí)踐［M］.北京：機(jī)械工業(yè)出版社，2013.FARBER R.CUDA application design and development［M］.Beijing：Machinery Industry Press：2013.（in Chinese）

［9］ 張淳民.干涉成像光譜技術(shù) ［M］.北京：科學(xué)出版社，2010.ZHANG Chun-min.Interferential imaging spectroscopy ［M］.Beijing：Science Press，2010.（in Chinese）

［10］ 楊露.干涉成像光譜技術(shù)中干涉圖處理的研究［D］.南京：南京理工大學(xué)，2007.YANG Lu.Study on interferometry processing of imaging spectroscopy interference［D］.Nanjing：Nanjing University of Science and Technology，2007.（in Chinese with an English abstract）

［11］ 韓剛.干涉光譜復(fù)原算法研究［D］.西安：中國(guó)科學(xué)院西安光學(xué)精密機(jī)械研究所，2012.HAN Gang.Research on the algorithm of inteferogram spectrum reconstruction［D］.Xi'an：Xi'an Institute of Optics and Precision Mechanics，2012.（in Chinese with an English abstract）