賈 東，溫 博，施 健，羅揚靜，王海濤

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

0 引言

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

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

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

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

1 ECA-B原理

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

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

(1)

(2)

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

第h段參考信號構成的滑動矩陣為:

(3)

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

基于此,每段目標回波信號自適應濾波權值系數可以轉化為一個求解凸優(yōu)化問題:

(4)

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

(5)

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

(6)

依照上述算法原理,可將ECA-B的數據處理流程總結如下:

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

② 利用每段參考信號構建滑動矩陣Vh;

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

④ 進行雜波抑制得到第h段的剩余信號。

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

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

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

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

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

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

(7)

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

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

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

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

① 在各分段中計算自相關矩陣Rh;

② 對自相關矩陣Rh進行求逆;

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

④ 進行雜波抑制得到第h段的剩余信號。

2.2 LDLT分解求逆

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

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

(8)

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

(9)

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

根據逆矩陣定義LB=E,推導可得:

(10)

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

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

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

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

④ 重復計算步驟②和③。

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

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

3 算法性能分析

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

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

為了驗證本文算法的雜波抑制效果,分別用Matlab和GPU進行雜波消除。圖5是Matlab對信號進行雜波抑制的距離-多普勒頻譜。圖6是GPU對信號做雜波抑制的距離-多普勒頻譜圖,雜波被消除,目標清晰可見?？梢钥闯龆邘缀跻粯?都得到了目標的距離-多普勒頻譜。圖7給出了Matlab和GPU計算結果的絕對誤差,誤差低于7×10-15。實驗結果說明了本文算法的有效性。

圖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

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

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

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

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

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

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

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

4 結束語

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