后向投影成像算法的GPU優(yōu)化方法研究?

班陽陽,張勁東,陳家瑞,邱曉燕

(南京航空航天大學(xué)電子信息工程學(xué)院,江蘇南京210016)

0 引言

合成孔徑雷達(SAR)是一種全天時、全天候的微波成像雷達,能夠?qū)︼w機、艦船和導(dǎo)彈實現(xiàn)超視距檢測[1],高分辨率的特點使它在軍用和民用領(lǐng)域有著不可替代的作用。隨著合成孔徑雷達成像技術(shù)的發(fā)展,各種高分辨率成像算法應(yīng)運而生。然而高分辨率帶來巨大的計算量成為某些成像算法實際應(yīng)用的瓶頸,其中最為典型的就是后向投影(Back Projection,BP)算法。BP成像算法是一種時域成像算法,其成像過程就是計算各方位時刻雷達平臺位置與目標(biāo)點的雙程延時,再將不同方位時刻對應(yīng)的回波信號進行相干累加,最后得出目標(biāo)函數(shù)的過程。BP算法在原理上不存在任何理論近似,因此成像效果更好;但是逐點成像計算的過程,使得該算法運算量巨大[2]。

近年來,GPU技術(shù)的迅速發(fā)展為大規(guī)模數(shù)據(jù)的計算與分析提供了一種新的技術(shù)手段。與CPU相比,GPU明顯具有成本低、性能高的特點。目前Fermi與Kepler架構(gòu)的出現(xiàn),使得GPU通用計算的性能有了飛速的提升[3]。在此基礎(chǔ)上,各種復(fù)雜的SAR成像處理算法在GPU上得到加速實現(xiàn)。文獻[4]對高分辨率星載SAR成像中ECSA(Extended Chirp Scaling)算法提出了一種加速方案,并在NVIDIA Tesla C2075上面實現(xiàn)了6～8倍的有效加速;中科院電子所針對距離多普勒(Range Doppler,RD)算法、CS(Chirp Scaling)算法及ω-k算法提出了不同的加速方案,并在NVIDIA K20C上面實現(xiàn)了每秒約36 MB采樣點的實時處理[5]。文獻[6]對復(fù)雜場景的SAR成像進行實時仿真,把計算機的相關(guān)技術(shù)應(yīng)用到SAR圖像仿真中去,充分利用GPU的可編程圖形管線,實現(xiàn)了三維復(fù)雜場景的建立和實時陰影的生成。

本文以后向投影成像作為算法基礎(chǔ),研究了GPU設(shè)計中的多流異步執(zhí)行技術(shù)、數(shù)據(jù)傳輸模式和計算速度與精度,提出了一種在聚束模式下針對BP成像的GPU優(yōu)化方案。該方案充分利用后向投影算法內(nèi)在的并行性,最大程度地利用GPU設(shè)備的計算資源,使得成像速度與GPU非優(yōu)化方案的實現(xiàn)相比有了近一倍的提升。

1 后向投影算法的并行化分析

BP算法主要包括距離向壓縮和反投影兩個部分。它假設(shè)發(fā)射波是沖激球面,通過回波在時域的相干疊加實現(xiàn)高分辨率成像。慢時間域相干疊加的成像結(jié)果可表示為

因此,后向投影算法的過程是一個點對點的圖像重建過程,在反投影過程中需要數(shù)量巨大的插值操作,導(dǎo)致運算量巨大。以Na×Nr的回波數(shù)據(jù)為例(其中Na和Nr分別為方位向和距離向點數(shù)),成像網(wǎng)格方位向和距離向點數(shù)分別為Nia和Nir,假設(shè)Na、Nia和Nir都為N,那么后向投影中插值的計算量就為O(N3)。除了插值運算本身之外,斜距計算和相位補償?shù)炔僮髟贑PU計算中也是非常耗時的,這些因素制約了該算法在CPU平臺上的應(yīng)用和發(fā)展。

GPU最大的特點是采用了SIMT(Single Instruction Multiple Threads)組織模式管理硬件中同時存在的線程,實現(xiàn)對所有數(shù)據(jù)的并行處理[7]。同時,GPU在線程的組織管理中采用了單線程(Thread)、線程塊(Block)以及線程格(Grid)的方式,將線程與GPU硬件中的流多處理器(Streaming Multiprocessor,SM)和流處理器(Streaming Processor,SP)進行了映射。

對后向投影算法的分析可知,不同網(wǎng)格點的反向投影操作過程完全相同并相互獨立,這使得后向投影算法在GPU上并行實現(xiàn)成為可能。我們可以通過為每個成像網(wǎng)格點分配一個獨立工作的線程,每個線程都單獨執(zhí)行斜距計算、插值、相位補償?shù)炔僮?最終完成對每個像素點的反投影計算。后向投影算法的GPU并行計算分配示意圖如圖1所示,對于Nia×Nir的成像網(wǎng)格,一共分配Nia×Nir個并行線程,其相應(yīng)的并行度為Nia×Nir。

圖1 后向投影算法的GPU并行計算分配圖

值得注意的是,在GPU中所謂的“并行”執(zhí)行并不嚴謹,實際情況下是“并發(fā)”執(zhí)行的[8]。由于硬件資源的限制,在理想條件下,并行運算所能獲得的加速比等同于并行度,但由于并行處理所涉及的GPU資源分配及調(diào)度等原因,在實際情況下加速比要低于并行度。盡管如此,并行度為Nia×Nir的反投影操作經(jīng)過并行處理后依然能獲得相當(dāng)高的加速比。

2 后向投影算法的GPU設(shè)計

2.1 算法設(shè)計流程

圖2是后向投影算法在CPU+GPU上實現(xiàn)的流程圖。后向投影算法實質(zhì)上是在CPU+GPU異構(gòu)平臺上完成的,GPU負責(zé)完成數(shù)據(jù)級并行的操作,包括距離向脈沖壓縮以及反投影,CPU則主要負責(zé)簡單的串行操作。

2.2 脈沖壓縮

圖2 利用CUDA實現(xiàn)基于CPU+GPU的后向投影成像算法流程圖

后向投影算法的距離向壓縮和傳統(tǒng)SAR成像算法一樣,包含3個部分：回波數(shù)據(jù)距離向FFT、參考函數(shù)和FFT后的數(shù)據(jù)相乘、相乘后數(shù)據(jù)進行IFFT。CUDA有自帶的FFT庫,即CUFFT庫,它是一個基于CUDA編程環(huán)境下的FFT庫,支持多個分段數(shù)據(jù)批量進行一維FFT以及IFFT運算,該庫對FFT以及IFFT操作能夠達到很高的運算性能[9-10],為距離向脈壓提供了極大的方便。

如圖2所示,首先將Na×Nr的SAR回波數(shù)據(jù)從CPU拷貝到GPU,接著對Na×Nr的數(shù)據(jù)進行距離向脈沖壓縮處理,主要包括對回波數(shù)據(jù)分段批量進行一維FFT、在編寫的核函數(shù)內(nèi)進行參考函數(shù)相乘以及將相乘結(jié)果分段批量進行一維IFFT等操作,最后將脈沖壓縮后的數(shù)據(jù)從GPU拷貝到CPU。其具體實現(xiàn)步驟如下：

(1)在主機端申請Na×Nr大小的內(nèi)存來存儲距離壓縮后的數(shù)據(jù),同時在設(shè)備端也申請相同大小的顯存并完成數(shù)據(jù)從主機端到設(shè)備端的拷貝。

(2)使用CUFFT庫中的cufft Handle制定FFT計劃plan,并使用cufft Plan1d制定FFT方案,即cufftPlan1d(&plan,Nr,CUFFT_C2C,Na),對Na×Nr長度的數(shù)據(jù)分段批量作FFT;再調(diào)用cufftExecC2C函數(shù)分段批量完成回波數(shù)據(jù)的一維FFT。

(3)編寫核函數(shù)完成FFT后的數(shù)據(jù)與參考函數(shù)的點乘,并將結(jié)果存入已經(jīng)分配好的顯存內(nèi)。

(4)與步驟(2)類似,再次調(diào)用cufft ExecC2C函數(shù)對步驟(3)中顯存內(nèi)的數(shù)據(jù)并行批量完成一維IFFT,完成距離向壓縮。

(5)考慮到GPU顯存容量,將脈沖壓縮后的數(shù)據(jù)拷貝到主機端,并釋放掉在顯存上申請的空間。

2.3 反投影

如圖2所示,反投影操作沿著方位向依次進行,方位時刻na(0≤na≤Na-1)的反投影操作包括na時刻脈壓數(shù)據(jù)以及必要參數(shù)從CPU拷貝到GPU、反投影核函數(shù)執(zhí)行。在同一方位時刻,對于Nia×Nir大小的成像網(wǎng)格,為核函數(shù)分配Nia×Nir個并行線程,核函數(shù)中的每一個線程獨立完成網(wǎng)格中一個點的反投影操作,包括距離計算、插值計算、相位補償,不同方位時刻的反投影結(jié)果相疊加并存儲在設(shè)備端,方位向操作完成后將最終的圖像數(shù)據(jù)拷貝到主機。

反投影操作主要包括斜距計算、插值和相位補償。對于Nia×Nir大小的成像網(wǎng)格,首先在設(shè)備上為成像網(wǎng)格分配一個Nia×Nir的存儲空間H0,并初始化為0,以存儲GPU上反投影的疊加結(jié)果,同時在主機上也分配一個同樣大小的內(nèi)存空間Hl以存儲從設(shè)備上拷貝回來的反投影數(shù)據(jù)。計算任意方位時刻雷達平臺與網(wǎng)格點之間的斜距R(t)求得雙程延時td,由BP算法原理可知,利用此延時構(gòu)造多普勒相位補償因子ej2πfctd來進行相位補償。在SAR成像處理過程中,距離壓縮后的數(shù)據(jù)都為離散信號,上面求得的雙程延時很難剛好匹配于采樣點,所以需要對距離壓縮后的回波數(shù)據(jù)進行插值重采樣,以實現(xiàn)精確的相位補償。每個方位時刻的反投影結(jié)果均與H0疊加并更新H0中的數(shù)據(jù),一直到方位向操作全部結(jié)束,H0圖像數(shù)據(jù)更新完成,并將H0數(shù)據(jù)拷貝回主機上事先分配好的內(nèi)存Hl里。

插值重采樣一般采用線性插值、最鄰近點插值以及8點sinc插值等。無論哪一種插值方式,都要通過斜距求出待插點位置,然后根據(jù)位置進行插值。由2.1節(jié)可知,在插值核內(nèi)為每個網(wǎng)格點的插值分配一個獨立的線程,線性插值和最鄰近插值只需要根據(jù)待插點位置便可一次性求出待插點的值,而8點sinc插值則需要在單個線程內(nèi)完成8次循環(huán)才能完成插值操作。目前主流的CPU主頻超過1 GHz、2 GHz,甚至3 GHz,而主流GPU的主頻在500～600 MHz,性能好的在1 GHz左右。從主頻來看,GPU執(zhí)行每個數(shù)值計算的速度并沒有CPU快,所以在GPU上編程要盡量避免單線程循環(huán)操作。基于以上分析,在GPU上實現(xiàn)8點sinc插值要比線性插值和最鄰近插值耗時長久。

3 GPU算法優(yōu)化

3.1 基于多流異步執(zhí)行技術(shù)的脈沖壓縮優(yōu)化

如圖2所示,距離向脈沖壓縮過程是先將SAR原始數(shù)據(jù)從CPU拷貝到GPU,完成距離向脈沖壓縮之后再將數(shù)據(jù)由GPU拷貝到CPU。為了隱藏數(shù)據(jù)在主機內(nèi)存與設(shè)備顯存之間的交互傳輸,充分利用GPU設(shè)備的計算資源,本文將利用CUDA多流異步執(zhí)行技術(shù)[11]來優(yōu)化此過程。首先創(chuàng)建4個流,并將SAR數(shù)據(jù)按方位向分塊并分別交由流0、流1、流2和流3處理,每個流單獨完成本塊數(shù)據(jù)的脈沖壓縮操作,包括“CPU→GPU數(shù)據(jù)拷貝、脈壓處理、GPU→CPU數(shù)據(jù)拷貝”三項任務(wù)。流之間同類處理任務(wù)不能并行執(zhí)行,不同種類的處理任務(wù)可以并行執(zhí)行。在整個任務(wù)處理過程中,除了首位部分時間4個流未完全并行外,中間處理過程都處于并行執(zhí)行狀態(tài)。由于GPU上的計算單元一直處于工作狀態(tài),從而有效隱藏了數(shù)據(jù)在主機和設(shè)備之間的傳輸時間。具體執(zhí)行過程如圖3所示。

圖3 多流異步執(zhí)行流程圖

3.2 反投影計算結(jié)構(gòu)優(yōu)化

如圖2所示,在反投影過程中習(xí)慣將脈壓數(shù)據(jù)沿方位向依次拷貝到GPU,即每個方位時刻的反投影過程都包括CPU→GPU數(shù)據(jù)拷貝以及核函數(shù)執(zhí)行兩個部分。反投影計算過程優(yōu)化之前的流程圖如圖4所示。

圖4 優(yōu)化之前執(zhí)行流程圖

在CUDA架構(gòu)中,只有使計算單元一直處于工作狀態(tài)才能充分利用GPU資源,而由傳統(tǒng)設(shè)計方案來看,在數(shù)據(jù)拷貝到GPU之前,GPU計算單元會一直處于等待狀態(tài),使得GPU資源沒有充分得到利用。為了解決以上問題,在反投影操作前將待用的數(shù)據(jù)拷貝到GPU內(nèi)存H2,然后在核函數(shù)中根據(jù)na計算出的偏移量在H2中選擇當(dāng)前方位時刻需要的數(shù)據(jù)進行反投影操作。這樣在沿方位向進行反投影操作的時候,不同方位時刻的核函數(shù)就會連續(xù)執(zhí)行,充分利用了GPU的計算資源。反投影計算過程優(yōu)化之前的流程圖如圖5所示。

圖5 優(yōu)化之后執(zhí)行流程圖

3.3 計算精度與速度優(yōu)化

本文中的GPU實驗平臺為NVIDIA Tesla C2075,共有448個CUDA cores,6 GB顯存以及1.15 GHz的時鐘頻率。Tesla C2075是NVIDIA公司為高性能計算量身打造,通過加入一顆Tesla協(xié)處理器,使其擁有高速的浮點計算能力,更適合用于以浮點運算為主的計算任務(wù),并且可以執(zhí)行符合IEEE-754標(biāo)準(zhǔn)的硬件雙精度計算。

成像質(zhì)量和成像速度是本次實驗的兩個主要指標(biāo),由于后向投影成像算法對精度要求相對較高,所以在滿足精度要求的情況下,最大限度的提高計算速度是本實驗追求的主要目標(biāo)。雙精度雖然能取得很好的成像效果,但其實現(xiàn)效率卻不太理想。反投影操作中要進行相位補償,需要計算sin和cos值,sinf(x)、cosf(x)及相應(yīng)的雙精度指令開銷非常昂貴,尤其是x絕對值較大時更是如此,而用CUDA自帶的__sinf(x)、__cosf(x)函數(shù)吞吐量卻能達到每個時鐘周期一個操作;考慮到sin和cos函數(shù)對小數(shù)位的精度要求不高,并且當(dāng)使用單精度時可以用__sinf和__cosf函數(shù)進行加速,所以為了提高計算效率,本實驗將使用混合精度來計算,即相位補償計算用單精度進行替換。本實驗中的波長λ為0.03 m,雷達與網(wǎng)格點之間的斜距R為10 000 m左右,相位補償因子4πR/λ的量級為107,此情況下數(shù)值較大不能直接使用__sinf和__cosf函數(shù)進行加速。基于余弦的周期性,先對相位補償因子進行如下處理：

處理完之后數(shù)量級減小,則可以直接用-sinf和-cosf函數(shù)進行加速。

4 CPU與GPU計算結(jié)果測試

本文中所用實測數(shù)據(jù)來自某型機載雷達所錄取,該雷達工作在X波段,以聚束模式采集數(shù)據(jù);點目標(biāo)仿真數(shù)據(jù)參數(shù)和實測數(shù)據(jù)參數(shù)相同,均為聚束模式數(shù)據(jù)。SAR仿真數(shù)據(jù)方位向大小Na為4 096,距離向大小Nr為32768,成像網(wǎng)格大小Nia×Nir為512×512,網(wǎng)格間隔大小為0.1 m×0.1 m。表1是本文用到的SAR實驗數(shù)據(jù)的相關(guān)參數(shù)。

表1 SAR成像參數(shù)

圖6為線性插值下算法優(yōu)化前后中心點目標(biāo)脈沖壓縮后的等高線圖。

線性插值下脈沖壓縮優(yōu)化前后耗時對比圖如圖7所示。

在多次計算求出平均耗時時間(不包括內(nèi)存/顯存的分配和釋放時間),(a)情況下脈沖壓縮耗時507.880 1 ms,(b)情況下耗時267.892 1 ms,優(yōu)化后的計算效率比優(yōu)化前提升了將近一倍。

反投影優(yōu)化前后耗時對比圖如圖8所示。

圖6 線性插值點目標(biāo)仿真圖

圖7 脈沖壓縮耗時圖

圖8 反投影耗時圖

在多次計算求出平均耗時時間,(a)情況下脈沖壓縮耗時1 936.326 8 ms,(b)情況下耗時1 100.356 1 ms,優(yōu)化后的計算效率比優(yōu)化前提升了將近一倍。

表2是基于CPU/GPU不同插值方法下成像速度的比較。

表2 基于CPU/GPU下不同插值方法成像時間對比(混合精度) ms

綜上結(jié)果分析可知,BP成像算法在GPU和CPU成像質(zhì)量相當(dāng),并且該算法在GPU上得到了150倍左右的加速,GPU優(yōu)化后的算法比優(yōu)化前速度提升了約一倍。另外由成像時間可知,線性插值速度最快,其次是最近鄰點插值,8點sinc插值速度最慢,這也驗證了之前的理論分析。

使用以上方案對實測數(shù)據(jù)進行成像處理,下面給出線性插值下實測數(shù)據(jù)分別在CPU和GPU上的成像結(jié)果,成像網(wǎng)格大小為512×512,網(wǎng)格間隔大小為0.8 m×0.8 m,如圖9所示。

本文通過圖像熵和圖像對比度準(zhǔn)則來判定圖像質(zhì)量。圖像熵的公式如下：

式中,0≤q≤Na-1,q為方位向索引,Na為方位向脈沖個數(shù);0≤k≤Nr-1,k為距離單元索引,Nr為距離單元個數(shù);D(q,k)為圖像的散射強度密度,且

圖9 實測數(shù)據(jù)成像圖

圖像對比度的公式如下：

式中,σ(·)表示方差,E(·)表示均值,|I(q,k)|2表示圖像各點的像素強度。對于給定的SAR圖像,E[|I(q,k)|2]表示信號能量,是一常數(shù)。

由表3可知,GPU成像質(zhì)量和CPU成像質(zhì)量相當(dāng),具有很好的成像效果。

表3 CPU與GPU成像質(zhì)量比較

5 結(jié)束語

本文主要研究了GPU設(shè)計中的多流異步執(zhí)行技術(shù)、數(shù)據(jù)傳輸模式和計算速度與精度,提出了一種針對BP成像的GPU優(yōu)化方案。該方案利用全新的GPU通用計算開發(fā)理念和編程模型,在低成本硬件的條件下,最大程度地加速了該算法的實現(xiàn),與GPU非優(yōu)化方案相比有了近一倍的速度提升。

CUDA環(huán)境下基于GPU的成像算法實現(xiàn),解決了單核CPU甚至是多核CPU計算能力不足的難題,最大限度地加速了復(fù)雜算法的實現(xiàn)。將算法的并行性最大限度地與GPU的編程框架相結(jié)合,充分利用GPU強大的計算能力和低廉的計算成本,開發(fā)出新型的SAR成像處理系統(tǒng),以應(yīng)對SAR成像領(lǐng)域目前面臨的新挑戰(zhàn)。

[1]陳輝,雷霆,閆達亮.天波超視距雷達作戰(zhàn)效能綜合評估研究[J].雷達科學(xué)與技術(shù),2014,12(2)：127-132.CHEN Hui,LEI Ting,YAN Da-liang.Comprehensive Evaluation of Combat Efficiency of Sky Wave OTHR[J].Radar Science and Technology,2014,12(2)：127-132.(in Chinese)

[2]聶鑫.SAR超高分辨率成像算法研究[D].南京：南京航空航天大學(xué),2010：14-17.

[3]LAI J,SEZNEC A.Performance Upper Bound Analysis and Optimization of SGEMM on Fermi and Kepler GPUs[C]∥2013 IEEE/ACM International Symposium on Code Generation and Optimization(CGO),[S.l.]：IEEE,2013：1-10.

[4]侯明輝.基于GPU的高分辨率星載SAR成像處理研究[J].電子科技,2013,26(10)：29-32,35.

[5]孟大地,胡玉新,石濤,等.基于NVIDIA GPU的機載SAR實時成像處理算法CUDA設(shè)計與實現(xiàn)[J].雷達學(xué)報,2013,2(4)：481-491.

[6]盧媛.基于GPU的復(fù)雜場景實時SAR仿真[D].上海：上海交通大學(xué),2010：27-54.

[7]張舒,褚艷利.GPU高性能運算-CUDA[M].北京：中國水利水電出版社,2009：20-24.

[8]VOLKOV V,KAZIAN B.Fitting FFT onto the G80 Architecture[D].Berkeley：University of California,2008：40.

[9]GOVINDARAJU N K,LLOYD B,DOTSENKO Y,et al.High Performance Discrete Fourier Transforms on Graphics Processors[C]∥Proceedings of the 2008 ACM/IEEE Conference on Supercomputing,[S.l.]：IEEE Press,2008：1-12.

[10]SANDERS J,KANDROT E,聶雪軍.GPU高性能編程CUDA實戰(zhàn)[M].北京：機械工業(yè)出版社,2011：40-45.