后向投影成像算法的GPU優化方法研究?

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

2 后向投影算法的GPU設計

2.1 算法設計流程

圖2是后向投影算法在CPU+GPU上實現的流程圖。后向投影算法實質上是在CPU+GPU異構平臺上完成的,GPU負責完成數據級并行的操作,包括距離向脈沖壓縮以及反投影,CPU則主要負責簡單的串行操作。

2.2 脈沖壓縮

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

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

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

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

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

(3)編寫核函數完成FFT后的數據與參考函數的點乘,并將結果存入已經分配好的顯存內。

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

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

2.3 反投影

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

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

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

3 GPU算法優化

3.1 基于多流異步執行技術的脈沖壓縮優化

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

圖3 多流異步執行流程圖

3.2 反投影計算結構優化

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

圖4 優化之前執行流程圖

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

圖5 優化之后執行流程圖

3.3 計算精度與速度優化

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

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

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

4 CPU與GPU計算結果測試

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

表1 SAR成像參數

圖6為線性插值下算法優化前后中心點目標脈沖壓縮后的等高線圖。

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

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

反投影優化前后耗時對比圖如圖8所示。

圖6 線性插值點目標仿真圖

圖7 脈沖壓縮耗時圖

圖8 反投影耗時圖

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

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

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

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

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

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

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

圖9 實測數據成像圖

圖像對比度的公式如下：

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

由表3可知,GPU成像質量和CPU成像質量相當,具有很好的成像效果。

表3 CPU與GPU成像質量比較

5 結束語

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

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

[1]陳輝,雷霆,閆達亮.天波超視距雷達作戰效能綜合評估研究[J].雷達科學與技術,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].南京：南京航空航天大學,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設計與實現[J].雷達學報,2013,2(4)：481-491.

[6]盧媛.基于GPU的復雜場景實時SAR仿真[D].上海：上海交通大學,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實戰[M].北京：機械工業出版社,2011：40-45.