基于GPU的MTD性能優化

楊千禾 袁子喬 扈月松

(西安電子工程研究所 西安 710100)

0 引言

雜波抑制是雷達信號處理中一個非常重要的環節。對于固定雜波，通常會使用結構簡單，運算量較低的動目標顯示(Moving Target Indication,MTI)實現對雜波的抑制。不過在實際情況下，硬件工藝水平的參差不齊和雷達工作環境的復雜多樣都會導致雷達目標回波出現起伏。這些現象將會影響MTI的雜波抑制效果，大部分的MTI對消器來說，其改善因子一般都被限制在20dB左右[1]，此外MTI對消器對于云雜波和雨雜波等氣象雜波的無法起到有效的抑制效果。現代雷達系統通常使用動目標檢測(Moving Target Detection,MTD)來克服MTI的不足，MTD通過多個帶通濾波器組對雷達回波進行脈沖多普勒處理[2]，再通過對濾波器組的輸出數據進行檢測，以致發現目標。相對于MTI，MTD能改善濾波器的頻率特性，使之更接近于最佳線性濾波，可以將改善因子提高至42dB，能夠檢測強雜波信號中的弱目標信號。

MTD濾波器組的實現方法一般分為兩類。第一類就是設計多階數的有限脈沖響應(Finite Impulse Response,FIR)濾波器組來對脈沖壓縮之后的數據進行濾波，這種方法能針對回波信號自適應的在雜波頻率處形成零陷，但是計算量過大，對于實時的信號處理系統的硬件計算能力的要求過高。另外一種方法就是用快速傅里葉變換(Fast Fourier Transform,FFT)來實現濾波器組，在計算資源有限的情況下能較好地完成MTD任務，但是需要更低副瓣的濾波器組來提升對于氣象雜波的抑制性能。更低副瓣會導致主瓣的展寬更加嚴重，可能導致主瓣發生明顯變形。

為了進一步提升雷達系統的計算能力，同時降低后續程序調試難度，實現雷達軟硬件的解耦，增強雷達的可重構性。使用圖形處理器(Graphics Processor Unit,GPU)代替DSP和FPGA作為雷達的計算核心，組成GPU+CPU的異構平臺作為新型的雷達信號處理系統架構。在使用GPU實現雷達信號處理過程中，MTD的優化效果相較脈沖壓縮(Pulse Compression,PC)和恒虛警檢測(Constant False Alarm Rate,CFAR)有明顯的不足，本文將著力解決這些問題并實現MTD在GPU計算平臺上的進一步優化。

1 MTD的實現方法及原理

1.1 基于MTI級聯FFT的MTD濾波器組

MTI的級聯FFT的MTD濾波器組是在FFT之前接一個二次對消器，它可以濾去最強的地物雜波，這樣就可以減少窄帶濾波器組所需要的動態范圍，并降低對濾波器副瓣的要求。由于DFT是一種特殊的橫向濾波器，所以濾波器組的系數可以按照DFT的定義來選擇，并采用快速算法FFT來實現MTD。

FFT的每點輸出，相當于N點數據在這個頻率上的積累，也可以說是以這個頻率為中心的一個帶通濾波器的輸出[3]。

根據DFT的定義，N組濾波器的權值為

wnk=e-j2πk/N,n=0,1,…,N-1;k=0,1,…,N-1

(1)

式(1)中，n表示第n個抽頭，k表示第k個濾波器，每一個k值決定一個獨立的濾波器相應，相應地對應于一個不同的多普勒濾波器響應。因此，第k個濾波器的頻率響應函數為

k=0,1,…,N-1

(2)

濾波器的幅頻特性為

(3)

每一個濾波器都有一定的副瓣，副瓣的大小決定著雜波抑制能力的大小。為了降低副瓣，通常在計算FFT之前，一般都需要使用一個窗函數w[k]對慢時間采樣數據進行加權來降低旁瓣的影響，加權之后的N點FFT可表示為

(4)

此時，第n個濾波器的幅頻響應可表示為

(5)

用窗函數進行加權的主要目的是減少濾波器旁瓣的相對幅度，使信號盡量從主瓣進入，同時主瓣盡可能窄，使濾波器獲得較陡的過渡帶。一般來說，使用了非矩形窗對信號進行加窗之后，達成對旁瓣抑制的同時也會導致主瓣寬度地增加。主瓣電平會減少，信噪比也會相應地降低。因此要針對雷達信號體制和雜波的環境，選用合適的窗函數進行加權處理。圖1為加海明窗之后的濾波器響應，可以看出在降低副瓣的同時，主瓣有所展寬。

圖1 FFT濾波器組頻率特性

1.2 點最佳多普勒濾波器組

采用FIR濾波器對雜波進行抑制，同時對目標信號實現相參積累，這里的處理關鍵在于FIR雜波抑制濾波器的設計，在設計FIR濾波器時需要根據雜波感知得到的雜波幅度、雜波譜寬來確定濾波器的凹口深淺及凹口寬度。考慮到FIR濾波器設計實時性可能無法滿足的問題，通過設置默認的多種不同凹口深度和寬度的組合，提前優化設計多組FIR濾波器，使用時只需根據凹口深度和寬度選擇相應的濾波器系數即可。

點最佳多普勒濾波器組只在所需的多普勒處理頻段中的某一點上達到最佳，而在其他頻點都是不匹配的[3]。多普勒濾波器就是用許多濾波器填滿感興趣的多普勒區域。通常實際運用中的多普勒濾波器采用N個橫向濾波器填滿多普勒區域，N等于處理的相干脈沖輸。

點多普勒橫向濾波器復數信號表示為

(6)

式(6)中A是幅度，wd是多普勒的角頻率，N是相參脈沖數，Tr是雷達脈沖間的間距。信號適量則可表示為sT=(s1,s2,s3,…,sN),其中

sn=Aejwd(n-1)Trn=1,2,…,N

(7)

Rs為歸一化信號協方差矩陣,Rs=E(ss*T)

(8)

歸一化協方差的元素mskj為

mskl=ejwd(k-1)Te-jwi(j-1)T=ejwd(k-j)T;
k=1,2,..,n;j=1,2,…,n

(9)

式(9)中的k為行數，j為列數。

根據自適應濾波器原理，長度為N的濾波器中第k個濾波器的權適量為

wk=R-1a(fk)

(10)

式(10)中fk(k=1,2,...,N)為第k個濾波器的通帶中心頻率，a(fk)為導頻矢量,為

a(fk)=[1,ej2πfkTr,ej2πfk2Tr,…,ej2πfk(N-1)Tr]T

(11)

R=Rc+σ2I為雜波加噪聲協方差陣，RC為雜波協方差矩陣，I為單位矩陣，σ2為噪聲功率[3]。R-1的作用就是使濾波器自適應地在雜波頻率處形成零陷，從而抑制雜波[3]。

圖2 64階點最佳多普勒濾波器組頻率特性

2 GPU并行計算的特點

2.1 CPU和GPU的功能異同

CPU和GPU是兩個相互獨立的處理器，有著不同的內部架構。CPU和GPU的組成如圖3所示。

圖3 CPU與GPU的內部結構

圖3中Control是控制器、ALU為算術邏輯單元、Cache是CPU內部緩存、DRAM就是內存。可以看到GPU設計者將更多的晶體管用作執行單元，而不是像CPU結構包含大量的控制單元與緩存[5]。從實際來看，CPU芯片空間的5%是ALU，而GPU則占到了40%。物理結構上的差異導致GPU在面對大批次計算時的表現比CPU要更加優秀。

對于控制密集型的任務，特別是邏輯復雜、分支眾多或深度迭代的計算任務，CPU處理效果要優于GPU；而GPU更擅長處理易于并行化的計算密集型任務。在GPU+CPU的異構計算平臺中，需要根據數據特性和算法特性對一整個項目進行分解，再給CPU和GPU分配契合自身硬件特性與計算特性的任務，這樣能有效提高大規模計算項目的性能。

2.2 CPU和GPU協同工作方式

GPU自身不具備獨立處理進程的能力，要實現GPU的并行計算必須通過PCI-E總線與CPU相連，組成異構計算平臺才能完成對任務的處理。因為一般都是由CPU向GPU傳輸待計算的數據并發布計算任務，所以稱CPU所在位置為主機端(Host),稱GPU所在位置是設備端(Device)。

程序的啟動和初始化都在主機端完成，待計算數據是由CPU先處理加載到主機端內存中。在這一階段，會完成通信端口的建立，主機內存的開辟，數據的接收/加載、組包和重排，完成數據計算前的所有準備。在主機端的準備完成之后，設備端也要在顯存中(顯卡的內存)中開辟相應的數據空間，完成對數據的接收準備。圖4為GPU+CPU架構的基本工作原理，這種架構可以充分揮GPU的高速并行計算能力，同時兼顧CPU在邏輯運算、流程管理方面的能力，通過合理的任務分配，實現系統的高性能擇優。

圖4 CPU+GPU協同工作方式

3 MTD算法在GPU上的工程實現

在實際的計算中， R-D平面矩陣在設備的存儲器中被轉化成一維的向量，所有的計算都是直接針對這個一維向量進行操作。所以為了達成較好的計算性能并得到正確的計算結果，在計算前需要根據數據的存儲格式和數據的讀取方式來判斷對是否需要對輸入結果進行轉置。一般情況下，PC之后的距離點數都較大，在CPU上實現矩陣轉置耗費的時間較長并會占用大量的計算資源。

在GPU上實現MTD算法時，可以通過CUDA程序的函數和異構計算平臺主機端與設備端分別初始化的特點，將一些步驟同時進行甚至省略，如圖5所示。

圖5 基于GPU的MTD算法

圖中較深陰影部分可以在整個程序的初始化或者是在主存向顯存傳輸數據時獨立完成，可以理解為HOST與DEVICE并行處理，隱藏消耗時間；圖中無文字部分為省略部分，沒有計算操作，并不消耗時間；淺陰影部分則是消耗計算資源的部分，也是本文優化MTD過程的重點。

3.1 加窗與求模數據讀取優化

在GPU上對數據進行點乘和求模，需要使用自編的核函數來實現并行運算。GPU的計算核心在核函數內規劃時是抽象成了thread,block和grid三級，thread作為最基本的計算單元包含于block中，多個thread組成的block在組成grid。在核函數中，向量的單元數據被分配給啟動的thread。因為在求模和點乘的計算過程中，相鄰單元直接互不影響計算結果，所以可以實現多線程并行計算。常用的多線程并行計算的情況如圖6所示。

圖6 常用多線程并行情況示例

相對于計算數據的規模可能無限增大，GPU的計算核心卻是固定的，無法實現單個thread僅處理一個待計算單位理想情景。在數據規模較大的情況下，一個thread處理多個待計算單位時，圖6的線程分配方法將會造成線程阻塞。為了避免thread阻塞，提高數據讀取和計算的效率，應該使用如圖7所示線程讀取方式。

圖7 防止線程堵塞的數據讀取方式

3.2 使用CUFFT實現FFT

CUDA的CUFFT庫函數提供了一種高效優化的快速傅里葉變換，可以用來實現MTI級聯FFT的濾波器組。CUFFT的配置使用句柄(Plan)來實現，該句柄初始化定義了包括待變換信號的長度，數據類型，在內存中的存儲形式等信息。使用句柄來進行初始化的方式為

cufftPlanMany(&cufftHandleplanname,
rank,n,inembed,istride,idist,
onembed,ostride,odist,CUFFT_C2C,batch)

rank表示每個進行FFT運算的信號的維度；n是一個數組，表示進行FFT的每個信號的行數，列數和頁數值上一個參數rank相互呼應；inembed表示輸入數據的三維點數，[頁數,列數,行數]、[列數,行數]或[行數]，若數據是一維數據這個參數會被忽略；istride表示每個輸入信號的相鄰兩個原始直接的距離；idist為兩個連續輸入信號的起始元素之間的間隔；onembed為輸出數據的情況，同inembed；ostride和odist與此同理，最后batch為進行FFT運算的信號個數。

傳輸到設備端的數據都是按照列優先的存儲規則來排布數據,可以通過改變參數istride與idist，或者ostride和odist來改變FFT的輸入輸出數據的格式，如圖8所示。通過句柄初始化直接根據需求對數據進行再次排序，在初始的過程中實現矩陣轉置，這一手段可以有效地優化MTD的運算時間。

圖8 FFT句柄初始化改變數據的排布格式

在完成對數據的初始化之后，調用CUFFT執行函數即可。

cufftExecC2C(planname,datain,dataout,
CUFFT_FORWARD)

其中四個函數結構分別對應之前初始化的句柄，輸入輸出數據的指針和FFT運算的方向。CUFFT_FORWARD是FFT運算，CUFFT_INVERSE是IFFT運算。

3.3 使用CUBLAS實現FIR濾波器組運算

在前文已經分析過，使用FIR濾波器組實現MTD能夠針對不同的回波信號，使濾波器自適應地在雜波頻率出形成零陷，從而靈活地抑制回波信號中的雜波。但是這種MTD算法的實現卻對硬件有著巨大的壓力，算法實現過程如圖9所示。

圖9 FIR濾波器組實現MTD

CUBALS庫對于這種規模較大的矩陣向量運算有著良好的加速效果，對于復數矩陣相乘C=alpha*A*B+beta*C有如函數

cublasCgemm(handle,transa,transb,
m,n,k,alpha,A,lda,ldb,beta,C,ldc)

transa和transb是表示矩陣A、B是否需要在計算前進行轉置;m,n,k就是計算矩陣的維度信息；lda,ldb,ldc是對應矩陣的主維度。

出于整個系統兼容性的考慮，CUBLAS也選擇使用列優先的存儲方式。如果在計算之前對矩陣做了轉置，那么矩陣的主維度就是矩陣的列數，此時的矩陣是按照C/C++以行優先的方式來存儲的；反之矩陣的主維度就是行數，保持了列有限的存儲方式。但是不論在計算之前是否對矩陣進行轉置操作，都只是改變了矩陣數據排列格式，結果矩陣C永遠都是按照列優先存儲方式，因此要獲得符合主機端矩陣排布的MTD結果，就需要對參與計算的矩陣進行轉置，得到行優先的CT。如果后續的CFAR和雜波圖計算仍然是在GPU是進行實現，則不需要此操作。表1總結了兩種情況下的參數設置規則，這是使用CUBLAS庫的關鍵。

表1 不同數據排布結果的參數選擇

4 性能優化分析

4.1 兩種MTD實現方案在GPU上的性能分析

在CPU+GPU異構平臺上對信號處理中的MTD算法進行設計，硬件使用兩塊6U-VPX板卡，其中一塊使用至強E5-2648L V4的14核CPU為主處理器，用于處理流程的控制和資源的調度；另一塊板卡為Tesla P6 GPU，其單精度浮點計算能力為6.1Tflops，傳輸帶寬為192 GB/s。

使用上述計算平臺實現某雷達的信號處理。該雷達測試數據包含5個波束(含一個匿影波束)，積累點數為32點，有效距離單元是800(對脈沖壓縮結果截取后產生)。將模擬產生的信號處理數據在GPU上運行100次之后取平均耗時，并在CPU上分別實現FFT的MTD和FIR的MTD，記錄下耗費時間，所得加速比如表2所示。

表2 MTD實現算法(開啟加速選項)

將GPU計算的MTD結果導出，與經過MatLab計算的結果進行比對，MTD結果誤差結果如圖10所示。需要注意的是在經過2次N點FFT之后的數據比原數據擴大了N倍，對結果進行誤差分析時需要將擴大的倍數考慮在內。最終求模輸出的MTD結果誤差在0.05%以內，考慮到在GPU上計算的數據類型為單精度(float)型，而MatLab默認的顯示數據類型為雙精度(double)，在進行數據轉換時會對結果產生微小的影響，因此這個級別的誤差基本不會影響后續的雜波圖建立和恒虛警檢測。

圖10 MTD結果誤差

可以看到使用FIR濾波器組的MTD性能要比使用FFT的MTD提升了42.4%，原本最費時的FIR濾波器組在GPU計算平臺上實現了性能的超越。

GPU的結構特性利于完成大規模的矩陣運算，針對數據讀取的優化與高速庫函數的調用進一步提升了濾波器組的計算速度。FIR濾波器組能夠讓雷達系統自適應的對雜波產生抑制，增強了系統的靈活性；但是這樣設計的得到的自適應的MTD濾波器會有較高的旁瓣電平，雖然不會造成主瓣的過度展寬以致主瓣變形，卻會造成不同濾波器之間的目標相互影響，帶來虛警。FFT的MTD對于氣象雜波的抑制效果不佳，需要使用更低副瓣的濾波器組；但是加窗獲得了更低的副瓣電平卻造成了主瓣的進一步展寬。

不過經過優化過后的兩種MTD方法耗時差異較小，經過論證分析，已在某雷達的信號處理程序中都進行了實現，根據具體情況選取更加合適的MTD方法。

4.2 使用CUBLAS實現FIR濾波的進一步優化猜想

在GPU高性能計算中，核函數與庫函數的啟用、初始化和計算數據的導入在整個計算過程中占據的時間花銷不能忽視。所以在使用GPU進行高性能計算時，為了達到最優的計算速度，整個應用的結構和數據的傳輸要根據實際情況作出相對應的調整。計算數據規模的線性增加不一定導致計算時間的線性增加，所以我們將使用固定的總計算量，不同的循環次數和啟動函數次數來測試MTD的計算時間。為了避免開啟加速優化項對多次測試的影響，在未加速情況下計算128點積累，R-D平面總數據量為128×100000個復數的FIR濾波過程，分別才去多段計算和一批次計算的方法，測試時間如表3所示。

表3 各種情況下的FIR濾波器組計算時間(未開啟加速選項)

從得出結果的后三項可以看出，并不是每次計算的矩陣維度小就一定能節省每次的運算時間，單次函數啟動和數據導入的開銷據決定著要選則合適的數據維度才能達到最好的計算優化；其次每次計算的數據維度越大，在計算的總耗時上就會有提升。所以在信號處理過程中，在不造成系統其它分系統阻塞或者使計算資源發生空閑的情況下，盡量在GPU計算前堆積足夠大量的數據，單次直接計算大量數據，減少函數啟動的次數，優化整體計算的時間。

5 結束語

本文從實現MTD的兩種方法入手，分析了FIR濾波器組和FFT計算實現MTD的優缺點。針對在較大數據情況下，使用FIR濾波器組實現MTD會導致計算資源的緊張和影響系統實時性的問題，提出了使用CPU+GPU異構架構來完成計算過程優化的方式。結果表明在GPU上實現MTD算法的處理結果與在CPU上所得結果誤差在0.05%內。經過算法優化和GPU加速選項優化之后，MTI級聯FFT的濾波器組在GPU上的計算效率相比在CPU上的串行計算提高了近32.64倍，而本身以運算量較大的FIR濾波器組算法相比CPU提高了近254.58倍。

FIR濾波器組能夠實現自適應的雜波抑制，增加了系統的靈活性，但副瓣較高；FFT濾波器組具有更低的副瓣電平，卻因為加窗抑制電平造成主瓣的進一步展寬。在經過優化之后，兩種MTD在GPU上的計算耗時來到了同一個水平。所以在以GPU為核心的雷達信號處理平臺中，可以根據具體需求擇優選擇兩種MTD實現方法。最后對于CUBALS庫在雷達信號處理過程中的進一步優化，提出了單次計算大量數據，減少函數啟動次數的建議。