一種基于GPU的SAR高效成像處理算法

孟大地* 胡玉新 丁赤飚

?

一種基于GPU的SAR高效成像處理算法

孟大地胡玉新 丁赤飚

(中國科學院電子學研究所 北京 100190)(中國科學院空間信息處理與應用系統技術重點實驗室 北京 100190)

合成孔徑雷達(SAR)成像處理是一項需要進行大量計算的處理任務。圖形處理器(GPU)具有數十倍于CPU的浮點計算能力以及傳輸帶寬，而CUDA技術的發展使得GPU能夠方便地進行通用計算。該文提出了一種在GPU上進行SAR成像的高效方法。與一般GPU處理方法相比，該方法使得處理過程中的CPU-GPU往返數據傳輸由4次減少到1次，而且同時利用了工作站上的CPU與GPU計算資源。實驗結果表明，該方法能夠帶來相對一般GPU處理方法2.3倍的處理效率提升，從而驗證了該方法的有效性。

SAR；CUDA；GPU；SAR成像處理

1 引言

合成孔徑雷達(Synthetic Aperture Radar, SAR)是一種全天時、全天候的微波遙感對地觀測手段。SAR圖像中包含豐富的地表特征信息，并具有穿透能力，因此在國民經濟以及軍事領域，具有光學遙感不可替代的重要作用。但與光學遙感相比，SAR圖像需要通過對SAR原始數據進行相干數據處理得到。常用的處理算法有以w-k算法為代表的頻域批處理算法(距離多普勒算法以及Chirp Scaling算法均屬此類)，以及反投影時域處理算法。

SAR圖像的方位向分辨能力有賴于對場景的長時間觀測期間的脈沖相干積累，SAR數據在兩個垂直維度上也存在耦合，因此在實際應用中，SAR數據量非常巨大；各SAR應用領域科學技術的飛速發展，對SAR分辨率、多波段、多極化等方面提出了更高的要求，這就使得SAR數據量成倍增加。因此，高效的SAR數據處理手段成為SAR技術發展的一項關鍵問題。

由于3D圖形領域的拉動作用，圖形處理器(Graphic Processing Unit, GPU)技術在近幾年來發展非常迅速。由圖1可見，在處理能力以及存儲器帶寬上，近七年間GPU所取得的進展與CPU相比超越約一個量級。GPU在運算以及傳輸速度上的飛速發展，在一定程度上應歸因于GPU主要專注于數據并行性任務的處理，而較少考慮了邏輯控制、條件判斷等CPU中不可或缺的重要功能。2006年11月，英偉達(NVIDIA)公司針對旗下主流型號GPU產品推出了一種通用并行計算架構CUDA(Compute Unified Device Architecture)，以及相關的并行編程模型和指令集，從而大大推動了GPU在高性能計算領域的廣泛應用。借助CUDA，在傳統CPU上運行的需要大量計算的代碼可以方便地移植到GPU上執行，而將代碼的邏輯控制以及條件判斷部分仍然由CPU執行。這種新的程序開發模式使得傳統上較為耗時的處理任務得到了數十倍甚至上百倍的速度提升。

利用CUDA這種高效的GPU編程模式，SAR數據處理軟件的速度有望得到大幅度提升。國內外已有少量該方面的研究成果。文獻[5-7]在GPU上實現了反投影處理算法，使得原本極其耗時的處理過程可以在可容忍的較短時間內完成；文獻[8-11]介紹了距離-多普勒算法在GPU上的實現方法。

當前支持CUDA編程的主流GPU產品的內存一般不大于4 GB，而在實際SAR數據處理過程中，所需處理的1景SAR數據往往超出單個GPU存儲能力。因此，在將現有SAR成像處理算法在GPU上實現時，除了要考慮原始算法中各大運算量代碼在GPU上的合理實現之外，還需要考慮如何高效地在CPU內存與GPU內存之間進行數據交互。

對于這種數據量超出GPU存儲能力的情況，文獻[8,9]通過減少處理脈沖數的方法強制減少1景SAR數據的數據量。由于處理所得每景SAR圖像的方位向兩端孔徑不全，因此不能達到全分辨率，在實際應用中需要截去，并且截去量與處理脈沖數無關，因此這種通過縮減方位景寬的方法將導致截去比例增加，降低處理效率。

當前GPU設備與主機之間主要通過PCI-E 2.0×16方式進行連接，根據在Tesla C1060上的實際測試結果，主機內存與GPU顯存之間進行1次2 GB數據的往返傳輸需要消耗時間約為1.1 s，而對2 GB float型存儲SAR數據利用單卡GPU1次成像(采用w-k成像處理算法，包括運動補償處理)時間約為2.5 s。由此可見，在成像處理過程中若有多次數據在主機內存與GPU顯存之間的往返傳輸，將會帶來效率的成倍下降。因此，需要盡可能減少主機內存與GPU顯存之間的數據傳輸次數，從而充分發揮GPU的高速浮點運算能力，提高運算速度。

本文對機載w-k成像處理算法在GPU上的實現進行了深入研究，提出了一種新的w-k算法在GPU上的實現方法方位時域分割法。該方法只需1次SAR原始數據導入GPU顯存以及1次float型SAR數據導出GPU顯存操作，從而大大降低了SAR數據在顯存傳輸上所消耗的時間。

另外，該方法將部分任務交由GPU處理，而數據導出后的剩余少量任務仍由CPU處理，這兩部分任務可并行執行。在配置4顆Intel Xeon X5550、16 G內存、1部Tesla C1060GPU的工作站上的實驗表明，在合理分配CPU核數與GPU設備數比例時，該方法在GPU端任務所花時間與CPU端相當。這時，1塊SAR數據的處理時間為GPU處理時間與CPU處理時間的最大值，從而在避免數據冗余傳輸的同時，充分利用了1臺主機上的所有計算資源。

本文第2節簡要介紹了w-k算法的基本步驟；第3節介紹了w-k算法的常規GPU實現方法；第4節提出了方位時域分割法，并給出了處理步驟和誤差分析；第5節利用點目標仿真SAR數據的處理結果驗證了算法的有效性和高效性；第6節對本文內容進行了總結。

2 w-k算法介紹

w-k算法是一種被廣泛采用的機載SAR成像處理算法。與其他常用算法(如距離-多普勒算法、Chirp Scaling算法)相比，由于其推導過程采用與SAR成像原理一致的雙曲距離歷程模型，并且除了載機速度不能沿距離向變化(在機載情況下該假設完全成立)以及插值誤差之外，未采取任何近似，因此w-k算法能夠適應各個波段的SAR數據處理，并能夠適應大方位波束角、大斜視角等情況。

另外，在實際SAR數據獲取過程中，由于載機不能嚴格保持勻速直線運動狀態，在成像處理過程中還需要加入運動補償處理環節(通過對距離壓縮后信號進行距離向重采樣以及相位補償實現)，以抵消載機運動的非理想性。

w-k算法的處理流程如圖2所示。該算法處理步驟如下：

(1) 首先對原始數據進行距離向壓縮(包括兩次距離向FFT以及距離向參考函數相乘)以及運動補償(包括插值以及相位生成、相位相乘)，并通過距離向FFT轉至距離向頻域；

(2) 通過方位向FFT(需要兩次額外轉置)轉至2維頻域；

(3) 在2維頻域進行距離遷移校正、2次距離壓縮、方位壓縮(通過距離向參考函數生成、相位相乘以及stolt插值實現)；

圖2 常規基于CPU的w-k算法處理流程圖

(4) 通過距離向IFFT轉至距離向時域；

(5) 最后通過方位向IFFT(需要1次額外轉置)得到時域SAR圖像。

3 w-k算法的常規GPU實現方法

利用CUDA編程模型將w-k算法移至GPU實現時，首要考慮的是將SAR數據傳入GPU顯存。當GPU顯存容量足以容納所處理SAR數據以及處理所需額外內存時，可以將所有w-k算法處理流程交由GPU處理，最后將處理結果SAR圖像傳入主機內存并進行存儲。

當對SAR系統的分辨率、測繪帶提出更高要求時，GPU顯存存儲容量往往不能滿足1景SAR數據的存儲。這時需要在處理過程中多次將處理中間結果分塊導入GPU顯存進行處理，之后再將處理結果導出GPU顯存，在進行轉置之后，再進行后續處理。

常規GPU處理步驟如下：

(1) 對SAR原始數據在方位向進行分塊，各塊分別傳入GPU顯存進行距離壓縮、運動補償以及距離向FFT，并將結果傳入主機內存；

(2) 對SAR數據在距離向進行分塊，各塊分別傳入GPU顯存進行方位向FFT，并將結果傳入主機內存；

(3) 對SAR數據在方位向進行分塊，各塊分別傳入GPU顯存進行距離向參考相位相乘、stolt插值以及距離向IFFT，并將結果傳入主機內存；

(4) 對SAR數據在距離向進行分塊，各塊分別傳入GPU顯存進行方位向FFT，并將結果傳入主機內存，拼接之后得到時域SAR圖像。

按照這種思路，除了在GPU端進行應有的w-k算法處理操作以及SAR原始數據傳入GPU顯存、SAR圖像數據傳出GPU顯存之外，還需進行3次額外的GPU顯存與主機內存之間的往返數據傳輸。

另外，每次將數據傳入主機內存之前需要將數據進行轉置處理(共需3次)，以適應后續步驟的處理需求。

這種基于GPU的常規w-k算法處理流程如圖3所示。可見，這種GPU端的w-k成像處理方法在原理上與常規w-k算法完全一致，而將所有主要的處理步驟分別移至GPU端實現，但增加了3次額外的主機內存與GPU顯存之間的往返數據傳輸。

4 方位時域分割法

根據SAR成像處理原理可知，對w-k成像處理所得到的SAR圖像進行方位向逆壓縮可以得到距離遷移校正后、方位壓縮前的SAR數據，這時場景中各目標的回波信號都被校正到1個距離單元，沿方位向排列。如果在成像處理開始之前對SAR原始數據在方位向分塊，每塊單獨進行w-k成像處理以及其他輔助處理，得到方位壓縮前數據，則可以將各塊處理得到的數據在方位向進行拼接，拼接得到的SAR數據即近似相當于對所有SAR原始數據統一進行成像處理所得到的方位壓縮前數據(該等效的近似性將在5.1節進行分析)。再對該數據進行方位壓縮處理，即得到所有SAR原始數據所對應的SAR圖像。

圖3 常規基于GPU的w-k算法處理流程圖

基于這種成像處理方式，本文提出了一種w-k算法在GPU上新的實現方法方位時域分割法。該方法將SAR原始數據回波在方位向分為若干塊(每塊數據量小于GPU顯存)，每塊交由GPU一次性進行方位壓縮前的所有處理步驟，并將處理結果導入主機內存并在方位向拼接，再由主機端進行方位壓縮，得到最終的SAR圖像。

方位時域分割法的處理步驟如圖4所示，處理步驟如下：

(1) 將SAR原始數據在方位向進行分塊(每塊數據量小于顯存容量)之后，各塊分別傳入GPU并進行常規w-k算法處理步驟(無需分塊)；

(2) 在方位向IFFT生成SAR圖像之前通過方位向參考相位相乘得到方位壓縮前信號；

(3) 將各塊處理結果傳入主機內存，并拼合為全景SAR數據的方位壓縮前數據；

(4) 通過方位向FFT、方位向參考相位相乘、方位向IFFT得到最終SAR圖像。

由此可見，與常規方法相比，方位時域分割法在避免了3次往返數據傳輸的同時，在GPU端及CPU端各增加了一次參考相位生成及相乘，在CPU端增加了兩次FFT。

5 誤差分析及實驗結果

5.1 方位時域分割法誤差的實驗分析

從原理上來講，在時帶積大于100時，駐定相位原理的近似性對SAR成像的影響可忽略不計。對于內存為4 G的GPU來說，所能容納的脈沖數一般能滿足該要求。這時，方位時域分割法對各塊處理所得的方位壓縮前時域數據拼接后近似等同于用w-k算法直接對SAR原始數據進行處理所得方位壓縮前時域數據。由于駐定相位原理近似性的理論分析較為困難，以下結合仿真對方位時域分割法的誤差特性進行研究。

圖4 方位時域分割法處理流程圖

對1塊距離、方位向分別為32768的仿真X波段SAR數據(有符號8位整型存儲，仿真中加入幅度約為3 m的運動誤差)進行成像及運動補償處理。仿真所用系統參數如表1所示。在該數據中，所設點目標回波的零多普勒位置位于中心距離單元的方位向中心，因此在利用方位時域分割法處理將仿真數據分為4塊進行處理時，該點目標回波信號的前后兩半將分別位于中間兩個分塊內。

表1仿真參數

Tab. 1 Simulation parameters

對表1仿真數據分別進行方位時域分割法處理以及基于CPU的w-k算法處理，點目標在相鄰兩個分塊內的方位壓縮前信號如圖5所示。這時信號總時帶積約為930，兩分塊內的局部信號時帶積約為465。由于駐定相位原理的幅度特性在信號邊緣有較大的震蕩以及低通效應，因此方位時域分割法的分塊操作導致各塊邊緣位置處有少量的信號泄漏，如圖5(a)所示；各段信號的邊緣泄漏在頻域也相應造成了較大的頻譜幅度震蕩，如圖5(b)所示。圖5中信號的相位特性以及方位壓縮后的方位波形如圖6所示，圖6(b)波形的各項壓縮指標如表2所示。由該仿真結果可見，方位時域分割法的分塊操作所引起的相位誤差以及幅度調制導致點目標的壓縮波形有所變化。這主要體現在方位向峰值旁瓣比(Peak Sidelobe Ratio, PSLR)以及積分旁瓣比(Integrated Sidelobe Ratio, ISLR)的惡化上(約0.1 dB)。在實際應用中，這樣的SAR圖像質量惡化程度通常是可以容忍的，而且可以通過加窗的手段進行抑制。

另外，由于運動誤差的方位空變性，兩個處理結果中方位向PSLR與理想值(-13.2 dB)相比惡化約0.1 dB。

5.2方位時域分割法性能分析

仿真所用工作站上配置了4顆4核Intel Xeon X5550 CPU以及1臺單卡NVIDIA Tesla C1060 GPU設備。基于CPU的w-k算法中使用了Intel MKL數學運算庫并通過Intel編譯器進行編譯，方位時域分割法使用CUDA 3.2編譯，其中主機代碼調用gcc編譯。在該工作站上，分別利用基于CPU的w-k算法、常規GPU方法、單CPU方位時域分割法、多CPU方位時域分割法4種方法對表1所示仿真數據進行成像及運動補償處理(方位時域分割法中將SAR原始數據分為4塊，每塊float型數據量為2 GB)。所花時間如表3所示(其中未包括SAR原始數據讀入主機內存時間以及SAR圖像文件輸出時間)。

圖5 方位壓縮前信號圖

圖6 方位時域分割法與w-k算法處理結果對比

表2 點目標仿真結果

Tab. 2 Simulation results of point target simulation

由仿真結果可見，在單核CPU配置情況下，與w-k算法的CPU實現相比，常規GPU實現方法效率提高將近30倍，但其中約60%時間花在數據在主機內存與GPU顯存之間的交互上，GPU的高速計算性能未能得到充分發揮。方位時域分割法避免了數據在主機內存與GPU顯存之間的冗余交互，但需要較長的CPU處理時間。

在配置多顆CPU使得方位時域分割法的CPU處理時間與GPU處理時間相當時，可在CPU進行后期處理的同時讓GPU進行下1塊SAR數據的前期處理，這時總處理時間為GPU處理時間與CPU處理時間的最大值(表3實驗結果中為10.79 s)。按照表3實驗結果，在配置14核CPU的情況下(1個CPU核心用于GPU控制)，相比基于CPU的w-k算法，常規GPU實現方法的效率提高約2.5倍，而方位時域分割法的效率提高約5.8倍。

表3 w-k算法3種實現方法處理效率對比

由以上仿真分析可知，方位時域分割法由于避免了冗余的主機與顯存之間的數據傳輸，并結合CPU資源進行并行處理，可以在約10 s內完成8 GB SAR數據的成像及運動補償處理，相比常規GPU實現方法，性能提高約2.3倍。

6 結論

本文提出了一種w-k算法在GPU平臺上的新的實現方法方位時域分割法。與常規GPU實現方法相比，新方法避免了主機內存與GPU顯存之間的多次冗余數據傳輸。實驗結果表明，在常規GPU實現方法中，冗余數據傳輸占用了成像處理的大部分時間，而方位時域分割法只需要在成像初始階段將SAR原始數據傳入GPU顯存，在GPU處理結束之后將中間結果傳回主機內存，與常規GPU實現方法相比減少了3次主機內存與GPU顯存之間的往返數據傳輸，成倍縮減了GPU處理時間。

另外，方位時域分割法將SAR成像及運動補償處理任務劃分為前期GPU計算部分和后期CPU計算部分，避免了常規GPU實現方法對CPU計算資源的浪費。按照這種任務分配方式，在進行批量SAR數據處理時，可在后期CPU計算的同時，GPU對下一塊SAR數據進行前期計算，提高了SAR數據處理效率。在實驗所用工作站的性能條件下，利用Tesla C1060單卡GPU與14核CPU進行處理時，GPU端處理時間與CPU端處理時間相當，可以較為充分地利用工作站的計算資源。

[1] Cumming I G and Wong F H. Digital Processing of Synthetic Aperture Radar Data: Algorithms and Implementation[M]. Boston: Artech House, 2005.

[2] Bamler R. A comparison of range-doppler and wavenumber domain SAR focusing algorithms[J]., 1992, 30(4): 706-713.

[3] Lars M H U, Hans H, and Gunnar S. Synthetic-aperture radar processing using fast factorized back-projection[J]., 2003, 39(3): 760-776.

[4] Nvidia. NVIDIA CUDA C programming guide [OL]. http:// developer.download.nvidia.com/compute/cuda/3_2/toolkit/docs/CUDA_C_Programming_Guide.pdf.

[5] Fasih A R and Hartley T D R. GPU-accelerated synthetic aperture radar backprojection in CUDA[C]. IEEE International Radar Conference, Arlington, VA, USA, May 2010: 1408-1413.

[6] Blom M and Follo P. VHF SAR image formation implemented on a GPU[C]. International Geoscience and Remote Sensing Symposium, Seoul, Korea, July 2005: 3352-3356.

[7] Hartley T D R, Fasih A R, Berdanier C A,.. Investigating the use of GPU-accelerated nodes for SAR image formation[C]. Proceedings of the IEEE International Conference on Cluster Computing, New Orleans, LA, USA, Sept. 2009:1-8.

[8] Liu Bin, Wang Kai-zhi, Liu Xing-zhao,.. An efficient signal processor of synthetic aperture radar based on GPU[C]. European Conference on Synthetic Aperture Radar, Eurogress, Aachen, Germany, June 2010: 1054-1057.

[9] Ning Xia, Yeh Chun-mao, Zhou Bin,.. Multiple-GPU accelerated range-doppler algorithm for synthetic aperture radar imaging[C]. IEEE International Radar Conference, Kansas City, MO, USA, May 2011: 698-701.

[10] Clemente C, Bisceglie M D, Santo M D,.. Processing of synthetic aperture radar data with GPGPU[C]. IEEE Workshop on Signal Processing Systems, Tampere, Finland, Oct. 2009: 309-314.

[11] 俞驚雷, 柳彬, 王開志, 等. 一種基于GPU 的高效合成孔徑雷達信號處理器[J]. 信息與電子工程, 2010, 8(4): 415-418. Yu J L, Liu B, Wang K Z,.. A highly efficient GPU- based signal processor of Synthetic Aperture Radar[J]., 2010, 8(4): 415-418.

[12] 張舒, 褚艷利. GPU高性能運算之CUDA[M]. 北京: 中國水利水電出版社, 2009: 120-124.

Zhang Shu and Zhe Yan-li. CUDA-High Perfermance Computing on GPU[M]. Beijing: China WaterPower Press, 2009: 120-124.

[13] Nvidia. Tesla C1060 specifications[OL]. http://www.nvidia. com/docs/IO/43395/BD-04111-001_v06.pdf.

[14] Intel. Intel xeon processor X5550 specifications[OL]. http://ark.intel.com/Product.aspx?id=37106.

[15] Fornaro G, Franceschetti G, and Perna S. On center-beam approximation in SAR motion compensation[J]., 2006, 3(2): 276-280.

[16] Prats P, Macedo K A C, Reigber A,..Comparison of topography-and aperture-dependent motion compensation algorithms for Airborne SAR[J]., 2007, 4(3): 349-353.

Efficient Algorithm for Processing SAR Data Based on GPU

Meng Da-di Hu Yu-xin Ding Chi-biao

(Institute of Electronics, Chinese Academy of Sciences, Beijing 100190, China)(Key Laboratory of Technology in Geospatial Information Processing and Application System,Chinese Academy of Sciences, Beijing 100190, China)

Data processing is time-consuming in the field of Synthetic Aperture Radar (SAR). Graphics Processing Units (GPUs) have tremendous float-point computational ability and a very high memory bandwidth, and the developing Compute Unified Device Architecture (CUDA) technology has enabled the application of GPUs to general-purpose parallel computing. A new method for processing SAR data using GPUs is presented in this paper. Compared with the nominal GPU-based SAR processing method, the number of data transfers between the CPUs and a GPU is reduced from 4 to 1, and the CPUs are exploited to cooperate with the GPU synchronously. By using the proposed method, we can speed up the data processing by 2.3 times, which is verified by the testing with simulated SAR data.

Synthetic Aperture Radar (SAR); Compute Unified Device Architecture (CUDA); Graphics Processing Unit (GPU); SAR data processing

TN957.52

A

2095-283X(2013)02-0210-08

10.3724/SP.J.1300.2013.20098

孟大地(1979-)，男，陜西渭南人，2006年于中國科學院電子學研究所獲得工學博士學位，2001年于西安交通大學獲得工學學士學位，現任職于中國科學院電子學研究所，副研究員，研究方向為機載合成孔徑雷達信號處理。E-mail: mengdadi@hotmail.com

胡玉新(1981-)，男，內蒙古赤峰人，2007年于中國科學院電子學研究所獲得工學博士學位，2002年于內蒙古大學獲得工學學士學位，現任職于中國科學院電子學研究所，副研究員，研究方向為合成孔徑雷達信號處理系統設計。E-mail: yxhu@mail.ie.ac.cn

丁赤飚(1969-)，男，研究員，博士生導師，現任中國科學院電子學研究所副所長，主要從事合成孔徑雷達、遙感信息處理和應用系統等領域的研究工作，先后主持多項國家863重點項目和國家級遙感衛星地面系統工程建設項目，曾獲國家科技進步一等獎、二等獎各一項。E-mail: cbding@mail.ie.ac.cn

2012-12-17收到，2013-04-07改回；2013-04-19網絡優先出版

國家大科學工程航空遙感系統地面數據處理與管理分系統項目資助課題

孟大地 mengdadi@hotmail.com