一種子孔徑ωK聚束SAR成像算法的GPU實(shí)現(xiàn)

司 軍，陳家瑞

(中國船舶重工集團(tuán)公司第七二三研究所,江蘇 揚(yáng)州 225101)

0 引 言

合成孔徑雷達(dá)(SAR)是一種高分辨率成像雷達(dá)，相比光學(xué)和紅外遙感成像，其全天時、全天候和透射性強(qiáng)等特點(diǎn)在軍用和民用領(lǐng)域發(fā)揮著重要作用[1]。隨著SAR成像技術(shù)的不斷發(fā)展，對分辨率提出了更高的要求，因此SAR回波數(shù)據(jù)量更大，成像算法更復(fù)雜，傳統(tǒng)基于CPU的處理已經(jīng)無法滿足成像需求。

近幾年，圖形處理器(GPU)發(fā)展迅速，其擁有強(qiáng)大的浮點(diǎn)計(jì)算能力和很寬的存儲器帶寬，為大規(guī)模SAR數(shù)據(jù)處理與分析提供了新的運(yùn)算平臺。2007年，計(jì)算通用設(shè)備架構(gòu)(CUDA)發(fā)布。CUDA的出現(xiàn)為充分利用GPU的強(qiáng)大性能提供了條件。CUDA對圖形硬件單元和應(yīng)用程序編程接口(API)進(jìn)行了封裝[2]，開發(fā)人員不需要太過關(guān)注圖形學(xué)硬件和編程接口，只需在掌握了GPU架構(gòu)和并行算法的基礎(chǔ)上使用類C語言進(jìn)行開發(fā)，這大大簡化了編程過程。SAR成像處理的步驟主要包括快速傅里葉變換(FFT)/逆快速傅里葉變化(IFFT)、復(fù)數(shù)相乘和插值等，都具有高度并行性，因此可以利用GPU強(qiáng)大的并行處理能力來實(shí)現(xiàn)算法的加速。當(dāng)前支持 CUDA 技術(shù)的GPU 產(chǎn)品的顯存一般不超過6 GB，尚不能對大場景SAR數(shù)據(jù)進(jìn)行一次性處理，因此需要對其進(jìn)行分塊操作[3]。

本文以O(shè)mega-K算法為基礎(chǔ)，通過研究聚束模式SAR成像特點(diǎn)，提出了一種針對GPU加速處理的子孔徑成像處理方案。該方案將全孔徑劃分成若干個子孔徑，然后在GPU上分別成像，最后將子圖像相干融合以獲得高分辨率圖像。

1 子孔徑ωK算法

子孔徑的劃分主要是依據(jù)脈沖重復(fù)頻率(PRF)和多普勒帶寬之間的關(guān)系[4]。由于聚束式SAR是通過天線長時間照射成像區(qū)域來獲得更長的相干積累時間，因此其多普勒帶寬也相應(yīng)變大。為了避免方位頻譜出現(xiàn)混疊，必須滿足PRF大于多普勒帶寬。

如圖1所示為子孔徑ωK算法中每個子孔徑的處理流程。將回波數(shù)據(jù)根據(jù)PRF與多普勒帶寬之間的關(guān)系劃分為若干個子孔徑，并分別成像。子孔徑處理時，需要通過補(bǔ)零來擴(kuò)展方位處理長度。每個子孔徑的處理過程包括以下步驟：

(1) 通過二維FFT將SAR回波數(shù)據(jù)轉(zhuǎn)換到二維頻域。假設(shè)距離方程是雙曲形式，則距離R0處的點(diǎn)目標(biāo)的相位為：

θ2df(fτ,fη)=

(1)

式中：fτ和fη分別為距離和方位頻率變量；f0為雷達(dá)中心頻率；c為光速；Vr為雷達(dá)有效速度；Kr為距離chirp調(diào)頻率。

實(shí)際處理中，在方位向FFT后會乘以補(bǔ)償相位φ1(t)，使得后續(xù)處理中產(chǎn)生的移位引起的方位處理長度展寬最小化[3]。

φ1(t)=exp(-j2πfdcη)

(2)

式中：fdc為多普勒中心頻率；η為方位時間。

(2)第1個關(guān)鍵聚焦步驟是參考函數(shù)相乘(RFM)。參考函數(shù)根據(jù)選定的距離Rref來計(jì)算。經(jīng)過參考函數(shù)相乘，參考距離處的目標(biāo)得到了完全聚焦。經(jīng)過RFM濾波后，二維頻域中的殘余相位近似為：

θRFM(fτ,fη)≈

(3)

(3) 第2個聚焦步驟是Stolt插值。它在距離頻域通過插值操作完成其它目標(biāo)的聚焦：

(4)

為了更好地理解Stolt插值，上式等號左邊的根號項(xiàng)可以根據(jù)泰勒公式展開得：

(5)

(6)

(4) 通過距離向IFFT轉(zhuǎn)換到距離多普勒域。在距離多普勒域，采用方位向Scaling與譜分析結(jié)合[5]的方式，乘以方位Scaling因子：

(7)

再通過方位向IFFT轉(zhuǎn)換到方位時域。

(5) 在方位時域乘以壓縮因子：

φ3(t)=exp(jπfηη2)

(8)

完成去調(diào)頻操作。最后經(jīng)過方位向FFT即可得到子孔徑圖像。

圖1 單個子孔徑處理流程

2 CUDA具體實(shí)現(xiàn)

子孔徑ωK算法實(shí)質(zhì)上是在CPU+GPU異構(gòu)平臺上完成的，CPU主要負(fù)責(zé)簡單的串行計(jì)算和邏輯事務(wù)處理，GPU則專注于大規(guī)模并行化數(shù)據(jù)運(yùn)算。如圖2所示為CUDA環(huán)境下子孔徑ωK算法的處理流程。CPU負(fù)責(zé)數(shù)據(jù)的讀寫和主機(jī)內(nèi)存與設(shè)備顯存之間的數(shù)據(jù)拷貝，GPU負(fù)責(zé)各個子孔徑的成像處理。最后將圖像數(shù)據(jù)拷貝到CPU進(jìn)行融合獲得高分辨率的圖像。

圖2 CUDA環(huán)境下子孔徑ωK算法處理流程

子孔徑ωK算法在GPU上的處理流程包括：5次FFT/IFFT，4次相位相乘和1次插值。其中5次FFT/IFFT可以通過CUDA自帶的CUFFT庫函數(shù)快速實(shí)現(xiàn)，該庫對FFT以及IFFT操作能夠達(dá)到很高的運(yùn)算性能，其余處理需要編寫對應(yīng)的Kernel函數(shù)[2]實(shí)現(xiàn)。

(1) 距離向FFT/IFFT

首先，使用cufftHandle制定1個FFT計(jì)劃plan；再利用cufftPlan1d函數(shù)創(chuàng)建1個CUFFT 句柄，對Na×Nr的數(shù)據(jù)分段批量實(shí)現(xiàn)FFT/IFFT，即cufftPlan1d(&plan,Nr,CUFFT_C2C,Na)，其中Nr為距離向采樣點(diǎn)數(shù)，Na為方位向采樣點(diǎn)數(shù)；最后調(diào)用cufftExecC2C函數(shù)完成距離向FFT/IFFT。

(2) 方位向FFT/IFFT

由于數(shù)據(jù)都是按先排距離向存儲的，因此進(jìn)行方位向處理前需要進(jìn)行轉(zhuǎn)置處理。轉(zhuǎn)置后的處理與步驟1基本相似。實(shí)現(xiàn)方位向FFT/IFFT后，再進(jìn)行轉(zhuǎn)置處理，將數(shù)據(jù)按先排距離向存儲。

(3) 復(fù)數(shù)相乘

子孔徑ωK算法中包含復(fù)數(shù)相乘的處理過程有：一致壓縮，方位向Scaling，方位向時移，方位向去調(diào)頻。

對于Na×Nr大小的數(shù)據(jù)，我們定義Na×Nr個并行線程來完成復(fù)乘操作。在CUDA中可以使用dim3類型的內(nèi)建變量threadIdx和blockIdx定義一維、二維和三維的索引來標(biāo)識線程。本文采用二維索引的方法，定義block的尺寸為(16,16)，即每個block包含256個線程。由于采樣點(diǎn)數(shù)通常為2的n次冪，因此不存在不能整除使得block數(shù)量為小數(shù)的情況。距離向處理時，定義grid的尺寸為(Nr/16，Na/16)，總的并行線程數(shù)為Na×Nr，每一個線程獨(dú)立完成1個采樣點(diǎn)的相位生成和復(fù)乘操作。方位向處理同上。

(4) 插值

本文中的Stolt插值采用sinc插值核，卷積核為8點(diǎn)，實(shí)現(xiàn)插值的關(guān)鍵步驟是找出待插點(diǎn)位置。對于8點(diǎn)插值，每個脈沖左右兩端各4個點(diǎn)無法覆蓋8個點(diǎn)，因此在處理時令它們的值與最接近的可以實(shí)現(xiàn)插值的點(diǎn)相同，其余點(diǎn)根據(jù)二分法尋找待插點(diǎn)位置。在插值kernel函數(shù)內(nèi)為每個插值點(diǎn)分配1個獨(dú)立的線程，則8點(diǎn)sinc插值只需要在單個線程內(nèi)完成8次循環(huán)即可完成插值操作。

待子孔徑成像完成之后，在CPU端將子圖像相關(guān)融合以獲得高分辨率圖像。

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

本實(shí)驗(yàn)使用的GPU是NVIDIA Tesla C2075，硬件具體配置如表1所示。

表1 CPU與GPU型號及配置

本實(shí)驗(yàn)中采用的SAR仿真數(shù)據(jù)方位向大小Na為16 384，距離向大小Nr為32 768，表2是本實(shí)驗(yàn)用到的SAR實(shí)驗(yàn)數(shù)據(jù)的相關(guān)參數(shù)。

表2 聚束SAR實(shí)測數(shù)據(jù)

根據(jù)表2數(shù)據(jù)，經(jīng)計(jì)算可得仿真數(shù)據(jù)的大小為4 G。本實(shí)驗(yàn)劃分了4個子孔徑，則每個子孔徑的數(shù)據(jù)大小為1 G。子孔徑處理時，需要通過補(bǔ)零來擴(kuò)展方位處理長度。因此，在子孔徑的方位向左右兩邊各補(bǔ)2 048個零，數(shù)據(jù)大小增加為2 G。顯存上需要分配2塊同樣大小的存儲空間，一塊用來存儲從主機(jī)內(nèi)存拷貝的數(shù)據(jù)，一塊用來存儲運(yùn)算處理的中間結(jié)果。另外，利用CUFFT庫函數(shù)實(shí)現(xiàn)FFT/IFF時，會根據(jù)數(shù)據(jù)大小調(diào)用cuffftPlan1d函數(shù)來創(chuàng)建一個plan配置數(shù)據(jù)。以上三項(xiàng)的存儲空間之和未超過顯存大小6 G。

4個子孔徑的運(yùn)算時間取平均后可得，在GPU端的平均耗時為11.735 s，而同樣大小數(shù)據(jù)在CPU端的平均耗時為956.796 s，速度提高約81倍。

子孔徑相干融合后的點(diǎn)目標(biāo)仿真結(jié)果如圖3所示。

圖3 子孔徑融合后點(diǎn)目標(biāo)仿真結(jié)果

下面給出實(shí)測數(shù)據(jù)相干融合后的成像結(jié)果，成像大小為4 096×4 096，如圖4所示。

4 結(jié)束語

根據(jù)聚束SAR成像模式特點(diǎn)，本文提出了一種適合GPU加速的子孔徑成像方案。首先將回波數(shù)據(jù)劃分為若干個子孔徑，分別進(jìn)行子孔徑成像，最后在CPU上將子孔徑圖像相干融合，獲得高分辨率圖像。通過對實(shí)測數(shù)據(jù)的成像處理，可以看出本文提出的子孔徑成像方案可以有效地實(shí)現(xiàn)算法的加速。從Tesla C2075上的實(shí)驗(yàn)結(jié)果可見，本方案有效地解決了GPU顯存小，不足以容納大場景SAR數(shù)據(jù)的限制，且取得了良好的成像效果，與CPU上的處理速度相比，有81倍的效率提升。本文中的子圖像相干融合仍然是在CPU上實(shí)現(xiàn)的，未來將把這一處理放在GPU上實(shí)現(xiàn)，實(shí)現(xiàn)算法的加速。

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