(武漢大學(xué)電子信息學(xué)院,湖北武漢430072)

外輻射源雷達(dá)滑窗擴(kuò)展相消算法并行實(shí)現(xiàn)

張 堅(jiān),萬(wàn)顯榮,劉玉琪

(武漢大學(xué)電子信息學(xué)院,湖北武漢430072)

雜波抑制是外輻射源雷達(dá)信號(hào)處理中的關(guān)鍵技術(shù)之一。ECA-S算法突破了ECA-B算法在實(shí)際數(shù)據(jù)處理中對(duì)分段數(shù)的限制,通過(guò)給每段分段信號(hào)增添滑窗,在增加分段數(shù)的同時(shí)保證了足夠的積累時(shí)間,具有更好的濾波效果,但這一改進(jìn)效果是以增大空間和時(shí)間復(fù)雜度為代價(jià)而得到。結(jié)合圖形處理器(GPU)數(shù)據(jù)吞吐量大、并行處理簡(jiǎn)單、適宜解決計(jì)算密集型問(wèn)題的特點(diǎn),提出一種適用于GPU處理的ECAS時(shí)域雜波抑制并行實(shí)現(xiàn)方法。實(shí)測(cè)數(shù)據(jù)驗(yàn)證了該算法的有效性,并滿足實(shí)時(shí)處理的需求。

外輻射源雷達(dá);時(shí)域雜波抑制;滑窗擴(kuò)展相消算法(ECA-S);并行實(shí)現(xiàn)

0 引 言

近年來(lái),外輻射源雷達(dá)[1-4]由于其巨大的潛在價(jià)值,在軍用和民用方面均引起了廣泛關(guān)注。它利用環(huán)境中已存在的照射源進(jìn)行目標(biāo)探測(cè),具有綠色環(huán)保、操作隱蔽的特點(diǎn)。

在外輻射源雷達(dá)中,監(jiān)測(cè)通道中不僅含有目標(biāo)的回波信息,同時(shí)還存在直達(dá)波和多徑雜波,后者的強(qiáng)度遠(yuǎn)遠(yuǎn)強(qiáng)于目標(biāo)回波,經(jīng)相干積累后雜波仍具有很高的距離副瓣和多普勒副瓣,使得距離多普勒譜的基底抬高,導(dǎo)致目標(biāo)被湮沒(méi)在回波譜中。因此,雜波抑制是外輻射源雷達(dá)信號(hào)處理中的關(guān)鍵技術(shù)之一。利用只含有直達(dá)波和多徑信息的參考信號(hào),可以從監(jiān)測(cè)信號(hào)中移除雜波信息。文獻(xiàn)[5-10]給出了雜波抑制的不同實(shí)現(xiàn)方法,它們各自具有不同的計(jì)算復(fù)雜度和處理效果。

在這些方法中,廣泛使用到的就是擴(kuò)展相消算法(Extensive Cancellation Algorithm,ECA)和它的分段版本——分段擴(kuò)展相消算法[6](Extensive Cancellation Algorithm Batches,ECA-B)。ECA算法利用參考信號(hào)構(gòu)建雜波子空間矩陣,通過(guò)消去監(jiān)測(cè)信號(hào)在雜波空間的投影分量來(lái)實(shí)現(xiàn)雜波抑制。ECA算法要求在整個(gè)積累時(shí)間內(nèi)計(jì)算自適應(yīng)濾波系數(shù),這使得算法構(gòu)建的雜波空間矩陣較大,最終導(dǎo)致算法的計(jì)算量和運(yùn)算存儲(chǔ)空間很大。相對(duì)地,ECA-B算法將監(jiān)測(cè)信號(hào)和參考信號(hào)分成多段,分別計(jì)算每一段的自適應(yīng)濾波系數(shù),最終濾波結(jié)果即為每一段計(jì)算結(jié)果的組合。ECA-B算法在計(jì)算量相當(dāng)?shù)那闆r下,降低雜波子空間的維度和最大運(yùn)算儲(chǔ)存空間消耗,增加了零陷寬度。并且通過(guò)分段處理,加快了自適應(yīng)濾波系數(shù)的更新速率,算法在非平穩(wěn)環(huán)境下具有更好的魯棒性。但是算法的零陷展寬由每段數(shù)據(jù)的ECA算法零陷決定,每段數(shù)據(jù)的積累時(shí)間不能過(guò)少,即ECA-B算法的最大分段數(shù)較少,制約了算法并行化實(shí)現(xiàn),不利于實(shí)時(shí)處理。為了保持ECA-B運(yùn)算存儲(chǔ)空間消耗較小以及在非平穩(wěn)環(huán)境下魯棒性強(qiáng)的優(yōu)點(diǎn),并消除分段帶來(lái)的影響,ECA-S[11]算法被提出。算法采用ECA-B算法的分段思想及相似的算法流程,但是在估計(jì)雜波子空間系數(shù)時(shí),在參考信號(hào)和監(jiān)測(cè)信號(hào)兩端分別各取一段滑窗信號(hào),用滑窗后參考信號(hào)構(gòu)建的雜波子空間和對(duì)應(yīng)的監(jiān)測(cè)信號(hào)來(lái)估計(jì)雜波子空間系數(shù)。算法在參數(shù)估計(jì)時(shí)需要作滑窗處理,使得整體運(yùn)算量大大增加。

CUDA適合解決計(jì)算密集型問(wèn)題,為基于GPU的ECA-S算法并行實(shí)現(xiàn)提供了條件。

目前,已有一些研究人員在GPU平臺(tái)上并行實(shí)現(xiàn)外輻射源雷達(dá)信號(hào)處理算法,并取得了較好的加速效果。文獻(xiàn)[12-14]使用CUDA對(duì)外輻射源雷達(dá)雜波抑制算法進(jìn)行并行處理。同時(shí),基于GPU并行算法實(shí)現(xiàn)也被用于其他科學(xué)計(jì)算領(lǐng)域[15-18]。

1 算法描述

1.1 ECA和ECA-B算法

假設(shè)照射源發(fā)射信號(hào)為s(t),則外輻射源雷達(dá)參考通道和監(jiān)測(cè)通道接收到的回波信號(hào)可分別表示為sref(t)和ssurv(t):

式中,τ為參考通道回波相對(duì)發(fā)射信號(hào)的延時(shí),βi,τi分別為N條多徑中第i條回波衰減和時(shí)延,其中i=0時(shí)為該多徑為直達(dá)波的情況,αm,τm,f m分別為M個(gè)目標(biāo)中第m個(gè)目標(biāo)的衰減、時(shí)延和多普勒頻移,nref(t),nsurv(t)分別為參考通道和監(jiān)測(cè)通道中的噪聲,一般可認(rèn)為是高斯白噪聲。

利用最小二乘算法(Least Square,LS)進(jìn)行相消濾波,即求下式的最小殘留信號(hào)能量:

容易求得

可得ECA算法雜波抑制后的信號(hào)為

ECA-B算法將參考和監(jiān)測(cè)信號(hào)等分為B段,每一段長(zhǎng)度為n L的信號(hào)分別進(jìn)行ECA濾波,可得到第i段自適應(yīng)濾波權(quán)矢量為

濾波后的信號(hào)表示為

信號(hào)總積累時(shí)間為T(mén),則每段ECA濾波的積累時(shí)間T B=T/B。對(duì)于T B的選取,有兩方面的考慮:一方面希望T B盡量長(zhǎng)以獲得更小的多普勒分辨率和減少自適應(yīng)損失;另一方面又希望T B盡量短從而使算法更加具有魯棒性,更能應(yīng)對(duì)環(huán)境變量快速變化的情況。這些考慮使得分段數(shù)B的選取受到限制,兩方面的要求都不能得到較好的滿足。

1.2 ECA-S算法

為了克服上述的限制,ECA-S算法通過(guò)引入滑窗信號(hào)來(lái)使得兩方面的要求不再對(duì)立:每段信號(hào)加上滑窗信號(hào)使得足夠長(zhǎng)的積累時(shí)間用于自適應(yīng)估計(jì);增大分段數(shù)使得自適應(yīng)系數(shù)的更新速率足夠快以適應(yīng)復(fù)雜的環(huán)境情況。這樣可以滿足在分段數(shù)足夠大的同時(shí)保障每段信號(hào)的積累時(shí)間充分長(zhǎng)的需求。

假設(shè)ECA-S算法將參考信號(hào)和監(jiān)測(cè)信號(hào)分為B段,每一段參考信號(hào)和監(jiān)測(cè)信號(hào)都向前和向后多取Ns/2個(gè)點(diǎn)作為滑窗。則第i段信號(hào)濾波后為

圖1為ECA-S算法的流程圖。對(duì)第i段信號(hào),首先利用加入滑窗后的求得自適應(yīng)權(quán)矢量α(i),然后利用構(gòu)建雜波空間矩陣X(i),最后依據(jù)式(8)計(jì)算得到濾波處理結(jié)果

圖1 ECA-S算法處理流程

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

2.1 R x的改進(jìn)計(jì)算

考慮到矩陣X(i)和的特殊結(jié)構(gòu),矩陣中元素與參考信號(hào)的對(duì)應(yīng)關(guān)系為

式中,X(i,j)表示矩陣X中第i行第j列的元素。

容易得到矩陣X的共軛轉(zhuǎn)置矩陣XH為

式中,?表示取共軛。

則有

由式(14)可知,關(guān)于R x(i,j)的計(jì)算有極大的計(jì)算冗余,直接計(jì)算需要L次復(fù)數(shù)浮點(diǎn)乘法和L-1次復(fù)數(shù)加法,而若利用對(duì)角元素的前一計(jì)算值R x(i-1,j-1),則只需要兩次復(fù)數(shù)乘法和兩次復(fù)數(shù)加法即可得到。從而可以通過(guò)計(jì)算矩陣第1行的值,然后迭代計(jì)算剩余K-1行的值,并最終得到R x,顯然此算法可以極大降低計(jì)算量。對(duì)比直接矩陣相乘,兩種方法計(jì)算量之比近似為K。改進(jìn)方法從原理上大大減少了計(jì)算量[12]。

在ECA及其改進(jìn)算法的濾波過(guò)程中,都可通過(guò)式(10)和式(11)從參考信號(hào)得到對(duì)應(yīng)雜波空間矩陣X中各元素的值,從而不需要直接構(gòu)建X,節(jié)省了存儲(chǔ)空間。特別是在ECA-S算法中,免去了構(gòu)建的空間,兩者都可從參考信號(hào)中提取出各元素的值。

2.2 分段ECA-S濾波并行計(jì)算

在使用CPU進(jìn)行ECA-S濾波時(shí),每一段數(shù)據(jù)都是循環(huán)讀取且依次進(jìn)行處理的。這是一個(gè)串行處理的過(guò)程,總處理時(shí)間是各段數(shù)據(jù)處理時(shí)間之和。而在CUDA對(duì)數(shù)據(jù)進(jìn)行并行實(shí)現(xiàn)中,可以將各段數(shù)據(jù)一起讀取同時(shí)進(jìn)行處理,這樣總處理時(shí)間只是各段處理時(shí)間中最長(zhǎng)的部分,理論上加速比近似為B。實(shí)際處理過(guò)程中由于GPU對(duì)于雙精度浮點(diǎn)數(shù)據(jù)的處理能力較弱,數(shù)據(jù)的傳輸速率也慢于CPU,同時(shí)由于GPU的內(nèi)核數(shù)有限,最大可并行度受不同顯卡內(nèi)核的限制,無(wú)法做到B段數(shù)據(jù)全部并行一起處理,使得加速比有所降低。這些硬件上的阻礙可通過(guò)更換更先進(jìn)的設(shè)備來(lái)優(yōu)化。

在CUDA編程模型中,CPU作為主機(jī)端,只負(fù)責(zé)數(shù)據(jù)的傳遞分配、運(yùn)行參數(shù)的配置等,而將大規(guī)模的數(shù)據(jù)計(jì)算交給GPU設(shè)備端進(jìn)行處理。GPU按照粒度粗細(xì)分為Grid,Block和Thread,Grid內(nèi)部Block間粗粒度并行,數(shù)據(jù)間無(wú)法直接通信,Block內(nèi)部Thread間細(xì)粒度并行,通過(guò)__syncthread()函數(shù)保證數(shù)據(jù)同步。在實(shí)際編程實(shí)現(xiàn)過(guò)程中,將每一步計(jì)算中可以并行計(jì)算的部分寫(xiě)成一個(gè)或多個(gè)核函數(shù),CPU控制整體的計(jì)算流程,GPU運(yùn)行Kernel函數(shù)并行計(jì)算實(shí)現(xiàn)。

2.3 程序流程

程序整體實(shí)現(xiàn)流程如圖2所示。首先讀取數(shù)據(jù),將參考和監(jiān)測(cè)通道的數(shù)據(jù)拷貝到設(shè)備端?？梢詫lock Dim.y設(shè)置為分段數(shù)B,第i段信號(hào)可以通過(guò)blockIdx.y來(lái)分別讀取。依據(jù)與第i段參考信號(hào)中數(shù)據(jù)的一一對(duì)應(yīng)關(guān)系,可以免去構(gòu)建雜波矩陣,節(jié)省了存儲(chǔ)空間。由于構(gòu)建X(i)和的參考信號(hào)都是從原始參考信號(hào)相同的位置截取的一段,不同的是前者比后者多截取了長(zhǎng)度為Ns的滑窗信號(hào),為計(jì)算方便,可以將滑窗信號(hào)的選取方式改為只在分段信號(hào)后截取,這樣可以用相同的指針表示兩者首元素的位置。實(shí)際處理結(jié)果表明這種滑窗選取方式對(duì)濾波效果無(wú)影響。的求逆采用了文獻(xiàn)[18]介紹的Gauss-Jordan原地求逆算法,相比基礎(chǔ)的Gauss-Jordan順序消去法,節(jié)省了一半的顯存空間和計(jì)算量。由于GPU的硬件結(jié)構(gòu)限制,Thread的最大數(shù)目不能超過(guò)規(guī)定數(shù)值,所以在計(jì)算矩陣相乘時(shí),一個(gè)Block計(jì)算一個(gè)元素,若需要累加的項(xiàng)數(shù)超過(guò)限制,則每個(gè)Kernel計(jì)算多個(gè)相乘的結(jié)果并累加,最后對(duì)Block內(nèi)所有線程進(jìn)行規(guī)約求和。同時(shí)由于的計(jì)算使用的是改進(jìn)方法,整個(gè)計(jì)算過(guò)程中的矩陣相乘部分都是矩陣乘上列向量,結(jié)果仍為列向量。可將block Dim.x設(shè)為列向量的行數(shù),Grid內(nèi)的x維用來(lái)計(jì)算相乘結(jié)果。

每一段監(jiān)測(cè)信號(hào)雜波抑制后的計(jì)算結(jié)果保存在sECA-S的對(duì)應(yīng)位置,所有線程計(jì)算完后即可得到ECA-S算法的結(jié)果。

圖2 程序流程圖

3 實(shí)測(cè)數(shù)據(jù)處理結(jié)果

2013年11月,武漢大學(xué)外輻射源雷達(dá)實(shí)驗(yàn)基地利用FM廣播信號(hào)進(jìn)行了目標(biāo)探測(cè)的外場(chǎng)實(shí)驗(yàn)。信號(hào)處理的軟硬件配置如表1所示。圖3為選取其中一段FM信號(hào)繪制的雜波抑制前的距離多普勒譜,其中信號(hào)中心頻率為103.8 MHz,帶寬為500 k Hz。采用ECA-S算法進(jìn)行并行計(jì)算,抑制距離元數(shù)為500,滑窗信號(hào)長(zhǎng)度取50 000,監(jiān)測(cè)通道數(shù)據(jù)總長(zhǎng)度為500 000,采用雙精度浮點(diǎn)計(jì)算,分段數(shù)為100時(shí),總耗時(shí)為0.861 s,相同設(shè)備下串行處理平均時(shí)間為18.531 s,總體加速比為21.5。圖4為ECA-S雜波抑制前后距離多普勒譜,雜波抑制前目標(biāo)被直達(dá)波和多徑雜波旁瓣掩蓋,雜波抑制后目標(biāo)凸顯。同時(shí)對(duì)相同的數(shù)據(jù)使用ECA-B算法濾波處理,分段數(shù)取10,圖5為ECA-B濾波后的距離多普勒?qǐng)D。對(duì)比可以看到,圖5中的3個(gè)目標(biāo)都能在圖4中找到對(duì)應(yīng),然而后者比前者多觀測(cè)到一個(gè)低速目標(biāo)。該實(shí)驗(yàn)現(xiàn)象驗(yàn)證了ECA-S算法相比ECA-B算法對(duì)低速目標(biāo)具有更好的觀測(cè)效果。圖6為ECA-S算法使用Matlab和GPU計(jì)算所得結(jié)果的絕對(duì)誤差,10-5量級(jí)的誤差完全可以忽略不計(jì)。

表1 軟硬件配置情況

圖3 雜波抑制前的距離多普勒譜

圖4 ECA-S抑制后的距離多普勒?qǐng)D

圖5 ECA-B抑制后的距離多普勒?qǐng)D

圖6 GPU與Matlab計(jì)算的絕對(duì)誤差

為了探究ECA-S算法處理時(shí)間與分段數(shù)的關(guān)系,選取不同的分段數(shù),分別使用C語(yǔ)言和CUDA進(jìn)行數(shù)據(jù)處理,處理時(shí)間如圖7所示?？梢钥闯?隨著分段數(shù)的增加,串行處理時(shí)間線性增長(zhǎng),而CUDA處理時(shí)間基本不變,對(duì)應(yīng)的加速比也隨之線性增加。這是由于滑窗信號(hào)的長(zhǎng)度相對(duì)分段信號(hào)較長(zhǎng),分段信號(hào)長(zhǎng)度變化對(duì)整體影響不大,從而每段信號(hào)的濾波處理時(shí)間基本不變,串行處理時(shí)間隨著分段數(shù)的變化而線性增加,而在GPU的最大并行處理限度內(nèi),并行處理時(shí)間不變。

圖7 不同分段數(shù)下ECA-S算法處理時(shí)間

4 結(jié)束語(yǔ)

本文針對(duì)ECA-S算法處理時(shí)間過(guò)長(zhǎng)而不利于實(shí)時(shí)化的問(wèn)題,使用CUDA在GPU上對(duì)算法進(jìn)行并行加速實(shí)現(xiàn)。實(shí)測(cè)數(shù)據(jù)驗(yàn)證了該算法的有效性。相對(duì)于常規(guī)方法,該方法能極大減少顯存需求,縮短計(jì)算時(shí)間,使得ECA-S算法能夠滿足外輻射源雷達(dá)雜波抑制的實(shí)時(shí)處理要求。

[1]萬(wàn)顯榮.基于低頻段數(shù)字廣播電視信號(hào)的外輻射源雷達(dá)發(fā)展現(xiàn)狀與趨勢(shì)[J].雷達(dá)學(xué)報(bào),2012,1(2):109-123.

[2]KUSCHEL H,O’HAGAN D.Passive Radar from History to Future[C]∥11th International Radar Symposium,Vilnius,Lithuania:IEEE,2010:1-4.

[3]HOWLAND P E,GRIFFITHS H D,BAKER C J.Passive Bistatic Radar Systems[M]∥Cherniakov M.Bistatic Radar:Emerging Technology.Weinheim:Wiley,2008:247-311.

[4]HOWLAND P E,MAKSIMIUK D,REITSMA G.FM Radio Based Bistatic Radar[J].IEE Proceedings:Radar,Sonar and Navigation,2005,152(3):107-115.

[5]CARDINALI R,COLONE F,FERRETTI C,et al.Comparison of Clutter and Multipath Cancellation Techniques for Passive Radar[C]∥IEEE Radar Conference,Boston,MA:IEEE,2007:469-474.

[6]COLONE F,O’HAGAN D W,LOMBARDO P,et al.A Multistage Processing Algorithm for Disturbance Removal and Target Detection in Passive Bistatic Radar[J].IEEE Trans on Aerospace and Electronic Systems,2009,45(2):698-722.

[7]PALMER J E,SEARLE S J.Evaluation of Adaptive Filter Algorithms for Clutter Cancellation in Passive Bistatic Radar[C]∥IEEE Radar Conference,Atlanta,GA:IEEE,2012:493-498.

[8]MELLER M,TUJAKA S.Processing of Noise Radar Waveforms Using Block Least Mean Squares Algorithm[J].IEEE Trans on Aerospace and Electronic Systems,2012,48(1):749-761.

[9]ZHAO Y D,ZHAO Y K,LU X D,et al.Block NLMS Cancellation Algorithm and Its Real-Time Implementation for Passive Radar[C]∥IET International Radar Conference,Xi’an:IET,2013:1-5.

[10]GUAN X,HU D H,ZHONG L H,et al.Strong Echo Cancellation Based on Adaptive Block Notch Filter in Passive Radar[J].IEEE Geoscience and Remote Sensing Letters,2015,12(2):339-343.

[11]COLONE F,PALMARINI C,MARTELLI T.Sliding Extensive Cancellation Algorithm for Disturbance Removal in Passive Radar[J].IEEE Trans on Aerospace and Electronic Systems,2016,52(3):1309-1326.

[12]陳偉,萬(wàn)顯榮,張勛,等.外輻射源雷達(dá)多通道時(shí)域雜波抑制算法并行實(shí)現(xiàn)[J].雷達(dá)學(xué)報(bào),2014,3(6):686-693.

[13]武勇,王俊,張培川,等.CUDA架構(gòu)下外輻射源雷達(dá)雜波抑制并行算法[J].西安電子科技大學(xué)學(xué)報(bào),2015,42(1):104-111.

[14]李曉波,關(guān)欣,仲利華,等.基于GPU的外輻射源雷達(dá)信號(hào)處理實(shí)時(shí)實(shí)現(xiàn)方法[J].系統(tǒng)工程與電子技術(shù),2014,36(11):2192-2198.

[15]鄧婕,張興浦,陳世友.基于GPU的信息融合并行方法研究[J].艦船電子工程,2016,36(6):35-37.

[16]沈聰,高火濤.使用GPU加速計(jì)算矩陣的Cholesky分解[J].計(jì)算機(jī)應(yīng)用與軟件,2016,33(9):284-287,305.

[17]賈春剛,郭立新,劉偉.基于GPU的并行FDTD方法在二維粗糙面散射中的應(yīng)用[J].電波科學(xué)學(xué)報(bào),2016,31(4):683-687.

[18]劉麗,沈杰,李洪林.基于GPU的矩陣求逆性能測(cè)試和分析[J].華東理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,36(6):812-817.

Parallel Implementation of Sliding Extensive Cancellation Algorithm for Passive Radar System

ZHANG Jian,WAN Xianrong,LIU Yuqi
(School of Electronic Information,Wuhan University,Wuhan430072,China)

Cancellation of clutter is one of the key signal processing techniques in passive radar.Sliding extensive cancellation algorithm(ECA-S)has broken through the limitations of extensive cancellation algorithm batches(ECA-B)in the number of batches during real data processing.By adding a sliding window to each segmented signal,a sufficient integrated time is ensured while the number of batches is increased,and a batter filtering effect is obtained.However,this improvement is achieved at the cost of increasing the space and time complexity.Considering the advantages of graphic processing unit(GPU)in high memory throughput,parallel processing,and computationally intensive problem,this paper proposes a parallel realization of ECA-S algorithm based on GPUs.The experimental results verify the effectiveness of the proposed algorithm.It also meets the demands of real-time processing.

passive radar;time-domain clutter suppression;sliding extensive cancellation algorithm(ECA-S);parallel implementation

TN958.97

A

1672-2337(2017)02-0115-05

10.3969/j.issn.1672-2337.2017.02.001

2016-11-02;

2016-12-02

國(guó)家重點(diǎn)研發(fā)計(jì)劃(No.2016YFB0502403);國(guó)家自然科學(xué)基金(No.61331012,61371197,U1333106,61271400);湖北省科技支撐項(xiàng)目(No.2015BCE075)

張 堅(jiān)男,1992年生,湖北孝感人,武漢大學(xué)電子信息學(xué)院電波傳播實(shí)驗(yàn)室碩士研究生,主要研究方向?yàn)槔走_(dá)信號(hào)處理。

E-mail:zhangjian6215@126.com

萬(wàn)顯榮男,1975年生,博士,教授、博士生導(dǎo)師,主要研究方向?yàn)橥廨椛湓蠢走_(dá)系統(tǒng)、高頻雷達(dá)系統(tǒng)及雷達(dá)信號(hào)處理。

劉玉琪男,1990年生,博士研究生,主要研究方向?yàn)槔走_(dá)系統(tǒng)、雷達(dá)信號(hào)處理。