一種LDLT分解的外輻射源雷達(dá)雜波抑制并行處理技術(shù)

賈 東，溫 博，施 健，羅揚(yáng)靜，王海濤

(1.桂林電子科技大學(xué) 信息與通信學(xué)院，廣西 桂林 541002；2.中國電子科技集團(tuán)公司第五十四研究所，河北 石家莊 050081)

0 引言

外輻射源雷達(dá)是一種利用第三方輻射源來探測目標(biāo)的雙多基地雷達(dá),具有隱蔽性好、抗干擾能力強(qiáng)、反隱身能力和防低空突襲能力強(qiáng)等諸多優(yōu)勢,受到了國內(nèi)外的廣泛關(guān)注[1]。但外輻射源雷達(dá)接收信號中不僅包括目標(biāo)回波信號,還包括能量很強(qiáng)的直達(dá)波和多徑雜波,從而導(dǎo)致微弱目標(biāo)回波信號被雜波副瓣掩蓋。因此需要對其進(jìn)行雜波抑制,消除強(qiáng)雜波對目標(biāo)檢測的影響[2-4]。

目前外輻射源雷達(dá)[5]主流的雜波相消算法是擴(kuò)展相消批處理算法(Extensive Cancellation Algorithm Batches, ECA-B)[6-8],其分段利用參考信號的多普勒和時延矩陣來構(gòu)建多徑和雜波空間,然后將回波信號投影到此空間的正交子空間從而實現(xiàn)雜波消除,該方法不僅可實現(xiàn)靜止雜波和慢動目標(biāo)雜波的消除,同時具有很好的穩(wěn)健性。但近年來隨著可用的照射源信號帶寬越來越寬(如傳統(tǒng)調(diào)頻(Frequency Modulation, FM)信號的有效帶寬為80 kHz[9],目前第5代移動通信技術(shù)(5th Generation Mobile Communication Technology, 5G)信號帶寬信號已達(dá)到50 MHz[10]),一方面大幅度提升了外輻射源的探測能力,另一方面導(dǎo)致在對消相同距離雜波的情況下,雜波對消算法所需的階數(shù)越來越大和計算復(fù)雜度越來越高。因此近年來,在外輻射源雷達(dá)實時跟蹤系統(tǒng)中,如何提升ECA-B的實時處理性能成為研究熱點。

目前高性能并行計算中,中央處理器(Central Processing Unit, CPU)+圖形處理器(Graphic Processing Unit, GPU)并行異構(gòu)由于成本低、開發(fā)難度小、體積小、便于攜帶等特點被廣泛應(yīng)用于外輻射源雷達(dá)信號處理系統(tǒng)中[11-14]。

文獻(xiàn)[11]利用GPU求解失配濾波因子,其將所有數(shù)據(jù)分段進(jìn)行并行處理。文獻(xiàn)[12]對拓展的LS算法進(jìn)行GPU加速實現(xiàn),其中矩陣相乘采用CUBLAS庫函數(shù)完成。CUBLAS庫中的乘法函數(shù)計算行列差距較大的矩陣時,會產(chǎn)生極大的空間冗余。文獻(xiàn)[13]采用雙GPU對ECA-B進(jìn)行加速,每個GPU消除一半回?fù)苄盘栔械碾s波。文獻(xiàn)[14]在實現(xiàn)過程中不需要構(gòu)造參考矩陣,節(jié)省大量空間,且改進(jìn)了自相關(guān)矩陣的計算,減少了大量計算冗余。上述傳統(tǒng)文獻(xiàn)都是以高斯消元法為基礎(chǔ)原理對矩陣進(jìn)行求逆,且文獻(xiàn)[14]調(diào)用統(tǒng)一計算設(shè)備架構(gòu)(Compute Unified Device Architecture, CUDA)流實現(xiàn)并行計算。CUDA流并不能使所有分段中同一計算環(huán)節(jié)的核函數(shù)并行,使得算法計算耗時長。雖然高斯消元法具有計算時間短的特點,但是需要在GPU與CPU之間多次傳輸數(shù)據(jù),在數(shù)據(jù)量非常大且為雙精度的情況下,數(shù)據(jù)傳輸時間將會大于計算時間,不利于算法的實時處理[15-16]。因此隨著輻射源信號帶寬變大而帶來的巨大的數(shù)據(jù)量,且數(shù)據(jù)精度要求變高,傳統(tǒng)算法已經(jīng)難以滿足高精度實時處理需求。為此,本文基于GPU多線程并行處理技術(shù)并結(jié)合ECA-B各段相同子模塊特性,使得ECA-B的各子模塊之間并行;然后,針對傳統(tǒng)ECA-B求逆過程中數(shù)據(jù)傳輸耗時問題,利用自相關(guān)矩陣共軛對稱特性提出一種基于LDLT的并行迭代求逆方法[17],通過2個CUDA核函數(shù)實現(xiàn)求逆處理,解決了雙精度數(shù)據(jù)的矩陣求逆過程中數(shù)據(jù)傳輸耗時問題,進(jìn)一步提升ECA-B的實現(xiàn)效率。實驗結(jié)果表明,與傳統(tǒng)算法相比,本文提出的算法具有更高的時效性和有效性。

1 ECA-B原理

假定參考信道接收信號為Sref(t),目標(biāo)回波信號為Ssur(t),輻射源信號為S(t)。

Sref(t)=S(t-τ)+nref(t),

(1)

(2)

式中:τ為參考信道回波相對于發(fā)射信號的時延,βi、τi分別為M條多徑中第i條回波衰減和時延,其中i=0時這條多徑信號為直達(dá)波;αk、τk、fk分別為K個目標(biāo)中第k個目標(biāo)的衰減、時延和多普勒頻移,nref(t)、nsurv(t)分別為參考信道和目標(biāo)回波信道中的噪聲,一般情況可以看作高斯白噪聲。

第h段參考信號構(gòu)成的滑動矩陣為:

(3)

式中:L為數(shù)據(jù)長度,H為數(shù)據(jù)分段數(shù),h代表分段的第h段,a為雜波抑制的對消階數(shù)。

基于此,每段目標(biāo)回波信號自適應(yīng)濾波權(quán)值系數(shù)可以轉(zhuǎn)化為一個求解凸優(yōu)化問題:

(4)

求式(4)是個最小二乘優(yōu)化問題,其解析解為:

(5)

那么,第h段的剩余信號可以通過第h段的回波信號減去第h段的滑動矩陣與權(quán)值向量的積得到。

(6)

依照上述算法原理,可將ECA-B的數(shù)據(jù)處理流程總結(jié)如下:

① 對參考信號Sref(t)和回波信號Ssur(t)分段;

② 利用每段參考信號構(gòu)建滑動矩陣Vh;

③ 利用第h段回波信號Ssur,h和滑動矩陣Vh求解權(quán)值向量;

④ 進(jìn)行雜波抑制得到第h段的剩余信號。

算法具體流程如圖1所示。

圖1 ECA-B計算流程Fig.1 The calculation flowchart for ECA-B

各分段之間雜波相消是相互獨立且處理流程相同,因此段間并行性很高。但是在傳統(tǒng)計算自相關(guān)矩陣過程中需要構(gòu)造參考矩陣Vh,采用CUDA流實現(xiàn)各分段并行。而CUDA流只能使數(shù)據(jù)傳輸和核函數(shù)并行,并不能使各分段同一計算環(huán)節(jié)之間并行計算。基于ECA-B各分段計算流程相同,本文采用段間并行算法使各分段同一計算環(huán)節(jié)在同一個核函數(shù)中并行實現(xiàn)。此外所有分段一起計算導(dǎo)致算法中的自相關(guān)矩陣求逆環(huán)節(jié)數(shù)據(jù)量非常大,而傳統(tǒng)的自相關(guān)矩陣求逆方法是高斯消元法,該方法進(jìn)行N階矩陣求逆時,需要傳輸3×N次數(shù)據(jù)。在采用雙精度浮點數(shù)數(shù)據(jù)形式時,這種設(shè)備間傳輸大批量數(shù)據(jù)大大增加了雜波抑制的處理時間[15],基于此問題,本文采用LDLT分解法只需傳輸一次數(shù)據(jù),降低了數(shù)據(jù)傳輸時間。

2 算法的改進(jìn)和并行實現(xiàn)

2.1 EDA-B的分段并行實現(xiàn)

目前CUDA并行主要用CUDA流實現(xiàn)并行計算,CUDA流的優(yōu)點在于數(shù)據(jù)傳輸?shù)耐瑫r調(diào)用核函數(shù),不能使核函數(shù)之間直接并行計算。此外,隨著輻射源數(shù)據(jù)量不斷增大,且高精度要求,使得ECA-B各子模塊處理時長變長,從而整體時長也會大大增大,因此對于ECA-B來說,并行度不高[18]。針對此問題,本文將算法中各分段的同一環(huán)節(jié)一起計算,提高了各分段間的并行度,減小數(shù)據(jù)量變大對算法實時處理產(chǎn)生影響。具體來說,此方法利用CUDA編程特點,通過線程索引控制每個線程訪問全局內(nèi)存中具體位置的數(shù)據(jù),利用這一特點實現(xiàn)線程塊分開計算各分段數(shù)據(jù)[19]。根據(jù)上述特點,把ECA-B所有分段的同一變量存儲在同一段連續(xù)存儲器內(nèi)。例如將數(shù)據(jù)分為10段計算,在存儲自相關(guān)矩陣Rh時,R1、R2、…、Rh、…、R10可按行優(yōu)先順序存儲,而矩陣Ch、Wh、Sh、Vh在計算機(jī)中按列優(yōu)先順序存儲,Rh的公式為:

(7)

在實現(xiàn)ECA-B的計算過程中,用線程塊獨立計算不同段的計算環(huán)節(jié)。例如以Wh計算為例,線程網(wǎng)分配10個線程塊,每個線程塊開辟10個線程。10個線程塊對應(yīng)計算ECA-B 10段中Wh,利用x維度索引矩陣的行數(shù)和列數(shù),y維度索引具體某一段中的矩陣,即線程塊Block0(線程號0～255)計算矩陣W0,線程塊Block0內(nèi)線程中迭代計算W0(i,i+k),線程塊Blockh(線程號256～511)計算矩陣Wh,線程塊Blockh內(nèi)線程中迭代計算Wh(i,i+k),以此類推。以為Wh例,一個核函數(shù)同時計算10個塊的示意如圖2所示。

圖2 GPU計算矩陣C示意Fig.2 Schematic diagram of matrix Cfor GPU calculation

同理,圖1中的矩陣Rh、Wh、Vh、Sh與圖2類似,也是將數(shù)據(jù)按一定規(guī)則存儲,用不同的線程塊去讀取和計算對應(yīng)的地址數(shù)據(jù),這樣即可實現(xiàn)塊間并行。此外,此并行框架由于訪問內(nèi)存并不完全為順序訪問,所以計算時間并不是每段單獨計算的k倍,而是各分段計算時長中最長的部分。算法實現(xiàn)流程如圖3所示。圖3中Vh為參考矩陣,Rh為自相關(guān)矩陣,Ssur,h為第h段回波信號,Wh為式(4)中的權(quán)值向量。ECA-B的數(shù)據(jù)處理流程總結(jié)如下:

圖3 算法實現(xiàn)流程Fig.3 Algorithm implementation flowchart

① 在各分段中計算自相關(guān)矩陣Rh;

② 對自相關(guān)矩陣Rh進(jìn)行求逆;

③ 利用第h段回波信號Ssur,h,Vh和求逆后的自相關(guān)矩陣Rh求解權(quán)值向量Wh;

④ 進(jìn)行雜波抑制得到第h段的剩余信號。

2.2 LDLT分解求逆

傳統(tǒng)的ECA-B實現(xiàn)過程中對自相關(guān)矩陣求逆通常采用高斯消元法。此方法需要多次進(jìn)行GPU與CPU之間的數(shù)據(jù)傳輸,導(dǎo)致數(shù)據(jù)傳輸時間大于數(shù)據(jù)計算時間。而外輻射源雷達(dá)的數(shù)據(jù)量非常大,使得高斯消元法的數(shù)據(jù)傳輸時間更大,因此在CPU與GPU傳輸帶寬小的設(shè)備中,高斯消元法不適合處理ECA-B。針對這一問題,本文根據(jù)厄米特矩陣Rh的對稱特性,采用LDLT分解法,通過2個GPU核函數(shù)實現(xiàn)自相關(guān)矩陣求逆。其中,一個核函數(shù)實現(xiàn)LDLT分解,另一個核函數(shù)實現(xiàn)矩陣求逆,并將其應(yīng)用到段間并行實現(xiàn)中,從而減少數(shù)據(jù)傳輸耗時。

LDLT分解法由Cholesky分解(Cholesky Factorization)法改進(jìn)而來,Cholesky分解法通過對三角矩陣L的列向量歸一化可以得到改進(jìn)后的Cholesky分解法為A=LDLH,其中D為對角矩陣,L為對角線上都是1的下三角矩陣,LH為L的復(fù)共軛轉(zhuǎn)置矩陣。當(dāng)對矩陣A求逆,即對角矩陣D-1的對角元素矩陣D對角元素倒數(shù),計算量很小,所以對矩陣A求逆主要就是對下三角矩陣L求逆。

(8)

根據(jù)矩陣乘法計算,再化簡可以得到:

(9)

式中:gij為過渡矩陣G的元素,lij為下三角矩陣L的元素,dj為對角矩陣D的元素。過渡矩陣G也是下三角矩陣,輔助計算下三角矩陣L和對角矩陣D。

根據(jù)逆矩陣定義LB=E,推導(dǎo)可得:

(10)

根據(jù)式(9),當(dāng)j=1時,矩陣G的第一列等于A的第一列,即gi1=ai1。當(dāng)計算矩陣G的第j列元素時,需要用到L的第j行的所有非零元素,而每行中的單個元素計算與其他行的元素互不影響,因此改進(jìn)型Cholesky分解法可以并行迭代實現(xiàn)。一個線程迭代分解得到矩陣G的一行數(shù)據(jù),線程中每次迭代計算需要用到上一次的結(jié)果,所以需要所有線程同步。而線程同步函數(shù)__syncthreads()僅保持同一個線程塊內(nèi)所有線程同步,當(dāng)矩陣維度大于線程塊內(nèi)最大線程數(shù)時,迭代放在CPU端控制,即CPU控制核函數(shù)調(diào)用,核函數(shù)迭代最大線程數(shù)次。具體GPU實現(xiàn)流程如下:

① 計算矩陣L的第一列,即gi1=ai1,1≤i≤C;

② 利用G和L前j-1行的值計算矩陣G第j列的值,線程同步;

③ 計算矩陣D和L的第j列的值;

④ 重復(fù)計算步驟②和③。

LDLT分解求逆實現(xiàn)示意如圖4所示。

圖4 LDLT分解求逆示意Fig.4 Inversion process based on LDLT decomposition

3 算法性能分析

為了評估本文算法的性能,在以下環(huán)境中進(jìn)行性能測試:CPU為Intel(R)Core(TM)i5-6200F CPU@2.30 GHz,內(nèi)存4 GB,圖形處理器為NVIDIA GeForce GT 940M顯卡,2 048 MB顯存;操作系統(tǒng)為Windows 10,配置了Matlab R2018a、Microsoft Visual Studio 2019和CUDA 11.0。

實驗中數(shù)據(jù)為實測接收的8路GSM信號,每路信號數(shù)據(jù)點數(shù)為80 000,取第一組作為參考信號,給每路信號加上6組不同的多徑干擾,再加上同一目標(biāo)回波信號,對8路信號進(jìn)行波束形成處理之后再消除雜波。采樣率fs為250 kHz,數(shù)據(jù)分段數(shù)為10,采用雙精度浮點數(shù)運算,重復(fù)實驗50次。為了驗證算法的性能,矩陣求逆模塊的驗證,將本文的LDLT算法與文獻(xiàn)[14]中的求逆算法進(jìn)行對比,ECA-B的整體性能與文獻(xiàn)[14]進(jìn)行對比。性能評價所采用的指標(biāo)為算法加速比。算法加速比為50次重復(fù)實驗下對比文獻(xiàn)[14]算法與本文算法的平均處理時間之比。

為了驗證本文算法的雜波抑制效果,分別用Matlab和GPU進(jìn)行雜波消除。圖5是Matlab對信號進(jìn)行雜波抑制的距離-多普勒頻譜。圖6是GPU對信號做雜波抑制的距離-多普勒頻譜圖,雜波被消除,目標(biāo)清晰可見。可以看出二者幾乎一樣,都得到了目標(biāo)的距離-多普勒頻譜。圖7給出了Matlab和GPU計算結(jié)果的絕對誤差,誤差低于7×10-15。實驗結(jié)果說明了本文算法的有效性。

圖5 Matlab雜波抑制距離-多普勒頻譜Fig.5 Clutter suppression range-Doppler spectrum in Matlab

圖6 GPU雜波抑制距離-多普勒頻譜Fig.6 Clutter suppression range-Doppler spectrum by GPU

圖7 GPU和Matlab計算的絕對誤差Fig.7 Absolute error calculated by GPU and Matlab

隨著輻射源信號的不斷增高,外輻射源處理的數(shù)據(jù)量也越來越大,且數(shù)據(jù)精度要求越來越高。針對此,文獻(xiàn)[14]中的矩陣求逆算法傳輸耗時問題使其不再適用于ECA-B中的自相關(guān)矩陣求逆。本文將LDLT算法與文獻(xiàn)[14]中的矩陣求逆算法進(jìn)行了對比,實驗數(shù)據(jù)由Matlab生成的雙精度Hermitian矩陣。實驗結(jié)果如圖8所示,由圖8可以看出,矩陣階數(shù)越大,LDLT算法耗時比文獻(xiàn)[14]中的矩陣求逆算法越少。矩陣階數(shù)從64階增長至512階,LDLT算法時間增大了162倍,而文獻(xiàn)[14]中的矩陣求逆算法增大了424倍,相比之下,LDLT算法更適用于處理大批高精度數(shù)據(jù)。

圖8 求逆實驗結(jié)果對比Fig.8 Comparison of inversion experiment results

本文段間并行實現(xiàn)將所有分段一起計算,而文獻(xiàn)[14]采用CUDA流的方法實現(xiàn)了ECA-B的段間串行算法。所有分段一起計算也導(dǎo)致了矩陣求逆環(huán)節(jié)的數(shù)據(jù)量增大了h倍,同時文獻(xiàn)[14]中的矩陣求逆算法的數(shù)據(jù)傳輸時間增大。基于此,本文在段間并行實現(xiàn)中將2種求逆算法再次進(jìn)行對比,當(dāng)對消階數(shù)為128時,段間并行算法中的矩陣求逆采用文獻(xiàn)[14]中的矩陣求逆算法,ECA-B處理時間為133 ms,而當(dāng)求逆方法采用LDLT分解法時,ECA-B處理時間則為45 ms。很明顯,在段間并行算法中,LDLT分解法優(yōu)于獻(xiàn)[14]中的矩陣求逆算法。本文算法與文獻(xiàn)[14]仿真結(jié)果對比,結(jié)果如圖9所示,最大加速比為3.5倍,大大縮短了算法運算時間。從圖9中可以看出,當(dāng)對消階數(shù)在128階以內(nèi)時,段間并行算法的處理時間維持在45 ms以內(nèi),而當(dāng)對消階數(shù)增大到256時,段間并行算法也僅僅只有296 ms,可以看出對消階數(shù)的增大對段間并行算法的影響較小,而隨著對消階數(shù)的不斷變大,對文獻(xiàn)[14]的段間串行算法影響較大。由此可以看出,本文的段間并行算法和LDLT求逆算法有著非常好的時效性能。

圖9 雜波抑制實驗結(jié)果對比Fig.9 Comparison of clutter suppression experiment results

為了驗證ECA-B消除雜波性能,本文將ECA-B與傳統(tǒng)算法中的SMI算法性能進(jìn)行對比,將目標(biāo)回波信號經(jīng)過雜波抑制處理之后進(jìn)行距離多普勒處理。圖10為ECA-B的距離-多普勒處理中多普勒維,圖11為SMI算法的距離-多普勒處理中多普勒維。從圖中可以看出,二者都可以清晰看出目標(biāo)信號,并且ECA-B在0多普勒通道產(chǎn)生明顯零陷,直達(dá)波和多徑雜波被完全抑制,雜波抑制效果較好。ECA-B和SMI算法的區(qū)別在于ECA-B將信號進(jìn)行分段處理,因此適合GPU將所有分段并行處理。

圖10 ECA-B的多普勒維Fig.10 Doppler dimension of ECA-B

圖11 SMI算法的多普勒維Fig.11 Doppler dimension of SMI algorithm

4 結(jié)束語

隨著輻射源帶寬的不斷增加,ECA-B對消階數(shù)也在不斷增加,使得算法所需處理的數(shù)據(jù)量不斷變大,對數(shù)據(jù)精度的要求也在不斷變高。ECA-B的分塊處理過程需要極大的時間和空間復(fù)雜度,本文將ECA-B的所有分段在GPU上統(tǒng)一進(jìn)行段間并行計算,并在GPU上實現(xiàn)矩陣的LDLT分解求逆,矩陣的LDLT分解求逆相對于Gauss-Jordan順序消去法有著很好的加速效果。ECA-B很好地滿足了外輻射源雷達(dá)對實時性的要求。