一種子孔徑ωK聚束SAR成像算法的GPU實現

司 軍，陳家瑞

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

0 引 言

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

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

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

1 子孔徑ωK算法

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

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

(1) 通過二維FFT將SAR回波數據轉換到二維頻域。假設距離方程是雙曲形式，則距離R0處的點目標的相位為：

θ2df(fτ,fη)=

(1)

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

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

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

(2)

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

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

θRFM(fτ,fη)≈

(3)

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

(4)

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

(5)

(6)

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

(7)

再通過方位向IFFT轉換到方位時域。

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

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

(8)

完成去調頻操作。最后經過方位向FFT即可得到子孔徑圖像。

圖1 單個子孔徑處理流程

2 CUDA具體實現

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

圖2 CUDA環境下子孔徑ωK算法處理流程

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

(1) 距離向FFT/IFFT

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

(2) 方位向FFT/IFFT

由于數據都是按先排距離向存儲的，因此進行方位向處理前需要進行轉置處理。轉置后的處理與步驟1基本相似。實現方位向FFT/IFFT后，再進行轉置處理，將數據按先排距離向存儲。

(3) 復數相乘

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

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

(4) 插值

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

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

3 實驗結果

本實驗使用的GPU是NVIDIA Tesla C2075，硬件具體配置如表1所示。

表1 CPU與GPU型號及配置

本實驗中采用的SAR仿真數據方位向大小Na為16 384，距離向大小Nr為32 768，表2是本實驗用到的SAR實驗數據的相關參數。

表2 聚束SAR實測數據

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

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

子孔徑相干融合后的點目標仿真結果如圖3所示。

圖3 子孔徑融合后點目標仿真結果

下面給出實測數據相干融合后的成像結果，成像大小為4 096×4 096，如圖4所示。

4 結束語

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

圖4 實測數據成像