基于GPU的MTD性能優(yōu)化

楊千禾 袁子喬 扈月松

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

0 引言

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

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

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

1 MTD的實現(xiàn)方法及原理

1.1 基于MTI級聯(lián)FFT的MTD濾波器組

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

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

根據(jù)DFT的定義，N組濾波器的權(quán)值為

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

(1)

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

k=0,1,…,N-1

(2)

濾波器的幅頻特性為

(3)

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

(4)

此時，第n個濾波器的幅頻響應(yīng)可表示為

(5)

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

圖1 FFT濾波器組頻率特性

1.2 點最佳多普勒濾波器組

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

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

點多普勒橫向濾波器復(fù)數(shù)信號表示為

(6)

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

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

(7)

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

(8)

歸一化協(xié)方差的元素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為行數(shù)，j為列數(shù)。

根據(jù)自適應(yīng)濾波器原理，長度為N的濾波器中第k個濾波器的權(quán)適量為

wk=R-1a(fk)

(10)

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

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

(11)

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

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

2 GPU并行計算的特點

2.1 CPU和GPU的功能異同

CPU和GPU是兩個相互獨立的處理器，有著不同的內(nèi)部架構(gòu)。CPU和GPU的組成如圖3所示。

圖3 CPU與GPU的內(nèi)部結(jié)構(gòu)

圖3中Control是控制器、ALU為算術(shù)邏輯單元、Cache是CPU內(nèi)部緩存、DRAM就是內(nèi)存?？梢钥吹紾PU設(shè)計者將更多的晶體管用作執(zhí)行單元，而不是像CPU結(jié)構(gòu)包含大量的控制單元與緩存[5]。從實際來看，CPU芯片空間的5%是ALU，而GPU則占到了40%。物理結(jié)構(gòu)上的差異導(dǎo)致GPU在面對大批次計算時的表現(xiàn)比CPU要更加優(yōu)秀。

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

2.2 CPU和GPU協(xié)同工作方式

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

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

圖4 CPU+GPU協(xié)同工作方式

3 MTD算法在GPU上的工程實現(xiàn)

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

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

圖5 基于GPU的MTD算法

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

3.1 加窗與求模數(shù)據(jù)讀取優(yōu)化

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

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

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

圖7 防止線程堵塞的數(shù)據(jù)讀取方式

3.2 使用CUFFT實現(xiàn)FFT

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

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

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

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

圖8 FFT句柄初始化改變數(shù)據(jù)的排布格式

在完成對數(shù)據(jù)的初始化之后，調(diào)用CUFFT執(zhí)行函數(shù)即可。

cufftExecC2C(planname,datain,dataout,
CUFFT_FORWARD)

其中四個函數(shù)結(jié)構(gòu)分別對應(yīng)之前初始化的句柄，輸入輸出數(shù)據(jù)的指針和FFT運算的方向。CUFFT_FORWARD是FFT運算，CUFFT_INVERSE是IFFT運算。

3.3 使用CUBLAS實現(xiàn)FIR濾波器組運算

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

圖9 FIR濾波器組實現(xiàn)MTD

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

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

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

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

表1 不同數(shù)據(jù)排布結(jié)果的參數(shù)選擇

4 性能優(yōu)化分析

4.1 兩種MTD實現(xiàn)方案在GPU上的性能分析

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

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

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

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

圖10 MTD結(jié)果誤差

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

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

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

4.2 使用CUBLAS實現(xiàn)FIR濾波的進一步優(yōu)化猜想

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

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

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

5 結(jié)束語

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

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