基于GPU加速的脈沖多普勒雷達信號處理

龔 昊,劉 瑩,馮建周,趙仁良,冷佳旭

(1.中國科學院大學計算機科學與技術學院，北京 100089；2.中國科學院大學數據挖掘與高性能計算實驗室，北京 100089；3.燕山大學信息科學與工程學院，河北 秦皇島 066004)

1 引言

雷達信號處理算法的高性能加速是現代雷達技術中至關重要的一環。傳統的雷達信號處理技術包括數字信號處理器DSP、現場可編程邏輯陣列FPGA和專用集成電路ASIC等[1 - 3]，雖然這些技術已經得到了大量的應用，但是它們都屬于專用設備，具有開發周期長、調試難度大等缺點，在開發和測試過程中需要耗費大量的資源。而GPU屬于通用設備，隨著NVIDIA、AMD等GPU廠商芯片技術的快速發展，CPU-GPU異構計算體系結構已經走進了千家萬戶，CUDA的出現更是大大簡化了程序員和科學家們對GPU進行編程的難度，基于CUDA的GPU開發環境逐步發展為成熟的貨架產品。

與CPU相比，GPU的眾核結構特別適合處理雷達信號這種大規模數據[4 - 6]，它的出現讓雷達信號處理領域煥發了生機。國內外關于使用GPU加速脈沖多普勒雷達這種目標探測雷達的研究比較少[7 - 19]，大部分研究成果都集中在SAR成像這種對實時性要求高的情景之中。文獻[8]針對脈沖多普勒雷達信號處理中的脈沖壓縮環節，采用FFTW(Fastest Fourier Transform in the West)和CUFFT(CUDA Fast Fourier Transform library)2種實現方式進行了性能基準測試，并分析了時域脈沖壓縮和頻域脈沖壓縮2種實現方式的性能差異。雖然實驗結果表明頻域脈沖壓縮可以取得更好的性能收益，但是并沒有針對脈沖壓縮模塊進行有針對性的性能提升與優化。文獻[11]研究了脈沖多普勒雷達中恒虛警率算法的GPU實現，雖然文中的方法降低了算法的復雜度，提高了并行性，但是總的計算量并沒有減少，提升空間有限。文獻[14]提出了一種基于GPU共享內存的恒虛警率算法，雖然使用共享內存可以優化程序的效率，但還是不可避免地引入了過多的重復運算。文獻[15]研究了脈沖壓縮模塊在不同數據處理規模和結構下，GPU硬件上的計算資源分布與基于CPU平臺的加速性能之間存在相關性，并基于實證數據提出了一種基于GPU的資源優化配置方法。但是，在大部分情況下，脈沖壓縮只是整個雷達信號處理系統中的一個流程，對其資源的過度傾向勢必會影響到其他模塊。

本文提出了一套完整的脈沖多普勒雷達信號處理GPU并行化方法。通過對脈沖多普勒雷達信號處理程序進行熱點測試來分析程序的性能瓶頸并進行針對性優化。本文的貢獻主要體現在以下幾點：(1)針對大規模雷達數據并行所造成的線程過度利用問題，提出了一種利用網格跨步并行技術[20]優化雷達信號處理的方法。該方法能有效地解決GPU創建和撤銷線程所造成的時間開銷，提高硬件資源利用率。(2)針對大規模雷達數據所造成的數據傳輸開銷，提出了一種利用多流異步處理技術加速脈沖多普勒雷達信號處理的方法。該方法能顯著降低大規模雷達數據的數據傳輸延遲，減小性能瓶頸。(3)針對恒虛警率CFAR(Constant False-Alarm Rate)算法中存在的冗余計算問題，提出了一種使用并行掃描來進行算法模塊數據預處理的方法。該方法能有效地避免重復計算，對于大規模的數據有明顯優勢。(4)設計了對比實驗，從性能測試和誤差分析等多角度來評估并行算法的準確性和實時性。

2 脈沖多普勒雷達信號處理算法分析

2.1 脈沖多普勒雷達信號處理流程

如圖1所示，脈沖多普勒雷達信號處理包含脈沖壓縮模塊、動目標顯示模塊、動目標檢測模塊和恒虛警率模塊4個基本功能模塊。

Figure 1 Flow chart of pulse Doppler radar signal processing圖1 脈沖多普勒雷達信號處理流程圖

2.2 熱點測試

將采樣率設置為1 GHz，產生40個雷達脈沖，總的采樣點數為9 600 000，對CPU程序利用Intel VTune進行熱點測試，結果如表1所示。

回波信號仿真模塊包括產生線性調頻信號和產生模擬回波信號2個部分，是程序中最耗時的部分。但是，它并不屬于真正的雷達信號處理系統的一部分，它是為了方便進行算法測試而實現的，在實際應用中這部分會被來自Rapid I/O等高速接口的真實數據所替代。

Table 1 CPU program hotspot test

初始化和內存釋放這2個步驟在整個程序的生命周期中只需要進行一次，而脈沖壓縮模塊、動目標顯示模塊、動目標檢測模塊和恒虛警率模塊這4部分是雷達信號處理系統的核心，需要進行針對性的加速。其中脈沖壓縮模塊和恒虛警率模塊是雷達信號處理系統中耗時最多的2個部分，本文將對這2個部分進行算法并行性分析。

2.3 算法并行性分析

2.3.1 脈沖壓縮模塊

脈沖壓縮模塊的作用是將信號的寬脈沖壓縮成窄脈沖，從而提高信號的探測距離和距離分辨率，其最終目的是提高信號的信噪比。脈沖壓縮的本質是計算回波信號對匹配濾波器的沖激響應。

脈沖壓縮模塊需要計算匹配濾波系數，它的計算方法是將原始線性調頻信號進行左右翻轉變換之后取共軛復數，其表達式如式(1)所示：

h(t)=s*(-t)

(1)

其中,h(t)為系統沖激響應函數,s*(t)為線性調頻信號復共軛。

計算沖激響應有2種方法，分別是時域處理和頻域處理。對于時域處理，直接計算回波信號對匹配濾波系數的離散卷積即可，如式(2)所示：

s0(t)=h(t)?sr(t)

(2)

其中,s0(t)為脈沖壓縮結果,sr(t)為回波信號。

對于頻域處理，先對回波信號和匹配濾波器的系數進行傅里葉變換，然后將這2種信號在頻域進行點乘，如式(3)所示：

S0(f)=FFT{h(t)?sr(t)}=H(f)·Sr(f)

(3)

其中,S0(f)為脈沖壓縮結果頻譜,H(f)為系統函數頻譜,Sr(f)為回波信號頻譜。

在這一過程中可以對H(f)進行頻域加窗，加窗之后的頻率響應輸出峰值會變小，信噪比也有一定的提升，將點乘的結果再進行傅里葉逆變換即可得到脈沖壓縮的結果，如式(4)所示：

s0(t)=IFFT(S0(f))

(4)

時域處理和頻域處理這2種方法都可以進行脈沖壓縮，但是頻域處理可以在計算過程中對頻域響應加窗。時域脈沖壓縮的計算效率不高，頻域脈沖壓縮相比較于時域脈沖壓縮具有更好的并行性。

在脈沖壓縮模塊中，多個脈沖可以同時進行脈沖壓縮，彼此之間不存在數據相關。對于單個脈沖來說，所涉及的操作也僅有FFT、IFFT和復數乘法這幾個基本算子。雖然脈沖壓縮算法是程序的瓶頸，但是在GPU并行化的基礎上可提升優化的程度有限。

2.3.2 恒虛警率模塊

恒虛警率模塊也叫做CFAR模塊，用于在信號處理的結果中判定目標。CFAR模塊會根據每一個采樣點周圍的幅度信息確定一個門限，超過這個門限的采樣點會被判決為目標。圖2展示了恒虛警率模塊的工作原理。

Figure 2 Schematic diagram of of CFAR module圖2 恒虛警率模塊原理圖

恒虛警率模塊主要由平方率檢波和雜波功率估計并判決這2個部分組成。其中平方律檢波部分是針對單個采樣點的變換過程，不存在與其他采樣點數據的相互依賴。而雜波功率估計并判決部分在計算每一個采樣點的門限值時需要用到相鄰采樣點的幅度信息，這就在并行之后的算法中引入了歸約與掃描計算。能否妥善處理算法中的歸約計算是制約算法性能提升的關鍵。

3 基于GPU的脈沖多普勒雷達信號處理算法并行化技術

3.1 基于網格跨步并行的細粒度并行化

GPU的并行度是有限的，其并行度由硬件執行模型來決定。以Maxwell架構為例，每個SM(Streaming Multiprocessor)單元上最多可調度2個尺寸為1 024的線程塊，當線程塊數量遠大于所有SM單元上可調度的線程塊數量時，就需要頻繁地進行線程的激活和撤銷，造成了不必要的時間開銷。基于網格跨步并行的細粒度并行算法如算法1所示。

算法1基于網格跨步并行的細粒度并行算法

輸入:雷達回波信號。

輸出:雷達信號處理結果。

Step1初始化核函數網格尺寸和線程塊尺寸。

Step2主機端調用核函數。

Step3對每一個GPU線程：

Step3.1計算當前線程在網格中的索引；

Step3.2計算當前線程塊在網格中的跨度；

Step3.3對每一個分配給當前線程的任務:

Step3.3.1計算雷達回波信號處理結果。

當使用傳統方法時需要啟用更多的線程塊才能完成大規模的計算任務，但是當硬件資源不充裕時，這種方法會造成額外的時間開銷，線程塊會被頻繁地調用。因此，在CUDA編程中應盡可能增加每一個線程的計算量，減少核函數調用的次數。圖3展示了數據元素數量大于網格中線程數的情況。圖3中下半部分的深色方塊代表當前活躍的線程，在這種情況下要想完成計算任務需要喚醒足夠多的線程塊。

Figure 3 Schematic diagram of the number of data elements is greater than the number of threads in the grid圖3 數據元素數量大于網格中線程數示意圖

Figure 4 Schematic diagram of single-thread grid striding parallism圖4 單線程網格跨步并行示意圖

圖4所示的網格跨步并行可以很好地解決該問題。對于單個線程而言，在網格跨度循環中，線程計算的第1個元素使用threadIdx.x+blockIdx.x*blockDim.x得出。其中，threadIdx.x表示線程塊內線程索引，blockIdx.x表示網格內線程塊索引，blockDim.x表示線程塊內線程數。然后，線程會按網格中的線程數blockDim.x*gridDim.x向前推進，其中gridDim.x表示網格內線程塊數，直到其數據索引超出數據元素的數量為止。

圖5展示了當所有線程均按照網格跨步并行的方式運作時，所有計算元素均被覆蓋的情況。

Figure 5 Schematic diagram of all threads grid striding parallism圖5 所有線程網格跨步并行示意圖

網格跨步并行所使用的核函數尺寸是固定的，一般根據具體的硬件執行模型來設定而不會根據數據規模自動調整。對于每一個線程而言，它會持續進行運算，當所有的運算結束之后線程才會被撤銷。

網格跨步并行還可靈活分配核函數對GPU的利用率，當以多線程形式控制GPU完成工作時，GPU會被每一個線程以CUDA流的方式管制。GPU的計算資源是有限的，它并不能保證所有CUDA流的線程塊同時運行，必須減小同一時刻每一個CUDA流上核函數對GPU的占有率才能實現多流計算的真正并行。

3.2 基于多CUDA流的粗粒度并行化

GPU的計算資源是有限的，而核函數一般使用盡可能多的計算資源以達到最佳效率，這就意味著2個控制流的計算不會重疊，雖然網格跨步并行可以靈活分配每一個控制流的計算量，但是這樣做并不能使GPU持續保持較高的利用率。

CUDA并行計算由CPU和GPU配合完成，在GPU計算時CPU可以以同步的形式等待GPU的計算結果，也可以進行和GPU沒有數據相關的計算，這樣就可以重疊CPU和GPU的計算，提高程序運行效率。

CPU和GPU的數據流和控制流通信必須經過PCIE總線，這是CUDA程序最大的性能瓶頸。當內存以鎖頁內存的形式進行分配時，可以避免操作系統請求分頁管理將內存以頁面的形式換出到磁盤，從而為內存和顯存之間異步的數據傳輸創造條件，主機端與設備端之間異步的數據傳輸可以與計算重疊起來，以隱藏數據傳輸延遲。

如圖6所示的多CUDA流并發是粗粒度并行方案。圖6中H2D代表主機端到設備端的數據傳輸，D2H代表設備端到主機端的數據傳輸。多CUDA流并發的重要用途是重疊數據傳輸和計算的時間，無論是主機端向設備端的數據傳輸，還是設備端向主機端的數據傳輸，都可以進行重疊；無論是CPU計算的時間，還是GPU計算的時間，都可以進行重疊。重疊可以增加程序的吞吐量，更進一步地降低延遲。

Figure 6 Multi-CUDA stream concurrency圖6 多CUDA流并發

3.3 基于并行掃描的數據預處理

在恒虛警率模塊中，對于每一個采樣點都要完成一次求和運算。采用并行掃描的方法進行預處理可以避免重復運算。在每一次求和運算中，涉及到的數據元素數量非常少，只有幾十個數據元素參與歸約。而采樣點的數據量非常龐大，甚至能夠達到吉比特數量級。掃描是一種典型的以空間換時間的方法，當有內存空間可復用時，掃描能夠帶來很好的收益，如式(5)所示：

sum(i)=sum(i-1)+a(i)

(5)

其中,sum(i)為掃描處理結果,a(i)為采樣點數據。

掃描預處理的時間復雜度為O(n)，其中n表示參與并行掃描運算的數據元素個數。當需要獲取子區間之和時，可以使用式(6)通過O(1)的時間復雜度計算得出。

ans(l,r)=sum(r)-sum(l-1)

(6)

以上展示了串行掃描的方法。當使用CPU完成掃描預處理時，必然會造成主機端和設備端之間的內存傳輸。Blelloch算法[21]是一種在設備端計算掃描的方法，該算法分為Reduce和Down-Sweep 2個過程。

圖7展示了Reduce過程，該過程是使用相鄰配對實現的并行歸約，為Down-Sweep過程準備了數據。Reduce過程利用所有線程計算出了所有2的整數次冪索引值的和。

Figure 7 Reduce process圖7 Reduce過程

圖8展示了Down-Sweep過程，該過程是Reduce過程的逆過程。圖8中實線箭頭代表求和運算，虛線箭頭代表賦值運算。

Figure 8 Down-Sweep process圖8 Down-Sweep過程

該并行算法的階躍時間復雜度為O(2 logn)，工作時間復雜度為O(n)。

4 實驗與評估

本文所使用的測試數據為仿真線性調頻信號，具體參數如表2所示。

Table 2 Radar simulation signal parameters

4.1 性能測試

CPU實現中的傅里葉變換操作使用了FFTW3.5加速庫，GPU實現中的傅里葉變換操作使用了cuFFT10.2加速庫。測試所使用的CPU和GPU配置如下所示：

(1)CPU：Intel Core i7-6700HQ。Intel Core i7-6700HQ CPU基于Skylake架構，具有4個物理核心，超線程8線程。它含有2個FMA(Fused Multiply-Add)單元，可同時發射2條256 bit融合乘加指令。Skylake架構的單精度浮點執行單元數的計算公式為：

2(融合乘加指令)× 2(FMA單元數)× 8(1條指令可以處理256/32=8個單精度浮點數)= 32

Intel Core i7-6700HQ CPU的理論單精度峰值浮點運算能力計算公式為：

2.6 GHz (默認主頻，超頻3.5 GHz) × 32(單精度浮點執行單元數)× 4(物理核心數)= 332.8 GFlop/s = 0.3328 TFlop/s = 0.3328萬億次浮點計算

(2)GPU：NVIDIA GeForce GTX 950M。NVIDIA GeForce GTX 950M GPU基于Maxwell架構，它含有5個SM單元，每個SM單元包含了128個CUDA核心，其理論單精度峰值浮點運算能力計算公式為：

1124 MHz(GPU Boost主頻)× 640(CUDA核心數量)× 2(單個時鐘周期內能處理的浮點計算次數)= 1438.72 GFlop/s = 1.43872 TFlop/s = 1.43872萬億次浮點計算

不同型號的CPU的物理核心數差異很大。在實踐中往往通過比較單核CPU和GPU的性能來衡量加速效果的優劣。

4.1.1 網格跨步并行實驗結果

幾乎所有的核函數都可以使用網格跨步并行進行優化。為了體現網格跨步并行的性能優勢，本文以核函數運算量最大的CFAR模塊為例來進行網格跨步并行的對比測試。使用不同的策略在32個雷達脈沖，每個雷達脈沖8 192個采樣點的條件下分別進行3次測試，結果如表3所示。

Table 3 Experiment results of grid striding parallism

GPU運算時間分為數據傳輸時間和內核時間2部分，在本算法中數據傳輸時間對加速比有顯著的影響。本文所統計的時間為這2部分時間之和，網格跨步并行的優化效果只體現在內核時間上。本實驗采用了較小的數據量，以更好地體現網格跨步并行的性能提升結果。

4.1.2 多CUDA流并發實驗結果

以整個雷達信號處理流程來進行多CUDA流并行的對比測試，在測試中使用1 GHz的采樣率，產生40個雷達脈沖，總采樣點數為9 600 000，調用脈沖壓縮、MTI(Moving Target Indication)、MTD(Moving Target Detection)和CFAR模塊100次來模擬算法模塊實時工作的情況。由于顯存容量的限制，本實驗中CUDA流的最大數量為8。在不同的CUDA流數量下進行實驗得到的實驗結果如表4所示。

在理想情況下多CUDA流并行最多可取得近3倍的速度提升，那是從主機端到設備端的數據傳輸、從設備端到主機端的數據傳輸和GPU運算這三者的時間完全重疊的情況。進行不同的任務數據傳輸和計算所占的比重不同，在本實驗中數據傳輸所占的比重不到1/10，并不能達到理想的時間重疊效果，多CUDA流并行的性能提升只有12%。

Table 4 Experiment results of multi-CUDA stream parallism

4.1.3 并行掃描預處理實驗結果

通過調整采樣率來增大數據規模，針對CFAR模塊采用非并行掃描和并行掃描2種策略，分別設計程序并進行實驗，結果如表5所示。

Table 5 Experiment results of parallel scan preprocessing

并行掃描通過對數據的預處理減少了CFAR模塊中門限計算時的重復計算，從根本上提升了算法的工作效率。

4.1.4 不同硬件平臺性能測試加速比

表6整理了不同硬件平臺上和不同采樣率下算法CPU實現和GPU實現所消耗的時間之比。實驗結果顯示，程序在PC機上達到了10倍左右的加速比，在嵌入式設備中達到了42倍左右的加速比。對于PC機來說，GPU運算能力是CPU運算能力的4.3倍；對于嵌入式設備來說，GPU運算能力是CPU運算能力的2.5倍，無論是哪種情況，加速比都超過了計算能力差距，這說明針對GPU的性能優化取得了很好的效果。

Table 6 Acceleration ratio of performance test

本次測試使用的PC機是HP Pavilion Gaming Notebook，其CPU和GPU的性能差距較小，其中NVIDIA GeForce GTX 950M是移動定制版，顯存帶寬為28.8 GB/s，較低的峰值浮點運算能力和顯存帶寬使得其在性能方面無法與桌面級顯卡相媲美。

在生產環境部署時使用的CPU為雙路Intel Xeon E5-2683 v3，每一路CPU都具有14個物理核心，其理論峰值浮點運算能力為896 GFlop/s。生產環境使用的GPU為NVIDIA GeForce GTX 1080 Ti，其理論峰值浮點運算能力為11 339.78 GFlop/s，顯存帶寬為484.44 GB/s，較大的性能差距使得算法在生產環境上能夠獲得更加出色的加速比。

Figure 9 Stage result error analysis圖9 階段結果誤差分析

4.2 誤差分析

圖9展示了算法CPU實現和GPU實現的相對誤差，其中脈沖壓縮模塊的絕對誤差小于3.5×10-5，動目標顯示模塊絕對誤差小于4.5×10-5，動目標檢測模塊絕對誤差小于6×10-4。從雷達信號處理中間結果的誤差分析圖中可以看出，相比于CPU實現，GPU實現具有極高的精度。

圖10展示了CFAR模塊的結果對比情況，圖10中橫縱坐標分別表示目標的距離信息和速度信息。從圖10中可以看到，使用CPU和GPU處理的結果其距離信息和速度信息完全相同，表明GPU實現產生了正確的結果。

Figure 10 Comparison of CFAR results圖10 CFAR結果對比

5 結束語

本文提出了一種基于GPU加速的脈沖多普勒雷達信號處理方法，通過網格跨步并行、多CUDA流并發和并行掃描等多種優化策略來實現加速。本文提出的方法既滿足了雷達信號處理大吞吐量和高實時性的要求，又充分發揮了GPU設備眾核并行的優勢。本文所闡述的方法最終在生產環境上達到了300倍的加速比。

致謝：

本文研究得到國家自然科學基金(71671178)、中國科學院大學優秀青年教師科研能力提升重點項目以及總裝備部裝備預先研究基金項目的支持。同時，感謝北京雷鷹科技有限公司對本文研究提供的幫助和支持。