隨機有限集推理動目標軌跡的相干累積檢測法

楊文彬, 李 旦, 張建秋

(復旦大學信息科學與工程學院電磁波信息科學教育部重點實驗室, 上海 200433)

0 引 言

在雷達信號處理過程中,基于多個脈沖多普勒濾波的動目標檢測(moving target detection, MTD)方法,在抑制強地面或氣象雜波、提高目標輸出信噪比等方面發(fā)揮著重要的作用[1-5]。然而,在應用MTD時,需保證目標在一個較長的相干處理間隔(coherent processing interval, CPI)內(nèi)距離和速度都不會發(fā)生變化,否則MTD的性能將因距離徙動(range migration, RM)[2,5]以及多普勒頻率徙動(Doppler frequency migration, DFM)[4-5]的存在而大幅降低。

Keystone變換[3]以及Radon-傅里葉變換(Radon Fourier transformation, RFT)[6]是最常見的RM處理方法。Keystone變換通過對接收信號的脈沖回波時間進行尺度變換來校正RM;而RFT通過對運動參數(shù)的聯(lián)合搜索來校正RM。但是,這兩種方法只適用于目標勻速運動的情況,不能有效地對同時存在RM和DFM的機動目標進行相干累積[4]。

為補償由目標高階運動導致的RM和DFM,文獻[7]在RFT的基礎上,提出了廣義RFT (generalized RFT, GRFT)。當目標的運動能以全局多項式模型準確描述時,該方法可通過對運動參數(shù)(距離、速度、加速度和加加速度等)的聯(lián)合搜索,實現(xiàn)理論上最優(yōu)的相干累積。但這樣復雜的搜索過程,在實際系統(tǒng)中幾乎不可實現(xiàn)[4]。為降低復雜度,文獻[8]提出了Radon改進高階LV分布(Radon modified high-order LV’s distribution, R-MHLV)法。該方法首先通過距離—速度的聯(lián)合搜索提取可能的目標軌跡,然后針對提取的軌跡,再利用LV分布來估計目標的運動參數(shù)。LV分布的引入減少了參數(shù)搜索的數(shù)目,從而在一定程度上減少了運算量的問題,但又會引入難以消除的交叉項[3]。

針對上述方法的運算復雜度較大的問題,相鄰交叉相關函數(shù)(adjacent cross correlated function,ACCF)法[4]、離散多相變換(discrete ploy phase transformation, DPPT)法[9]、對稱自相關函數(shù)(symmetric autocorrelation function,SAF)法[10]等被相繼提出。這些方法一般依賴于精心設計的某種自相關函數(shù),以期消除多項式運動模型的部分高階運動項(例如與加加速度有關的三階項)。雖然這些方法可在一定程度上減小相干累計的復雜度,但是自相關函數(shù)的非線性卻又使得它們在低信噪比條件下表現(xiàn)不佳[3-5]。

隨著技術的進步,非合作機動目標在一個較長的CPI內(nèi)已經(jīng)可以呈現(xiàn)出復雜的運動模式[11]。這樣,假定目標的運動可由全局多項式模型描述,在面對具有復雜運動的目標時,其性能將因模型的失配而惡化[11-12]。為解決上述問題,文獻[11]提出了短時GRFT(short time GRFT, STGRFT)法。該方法假設目標的運動可由有限個分段多項式描述,并對目標所有參數(shù)及它們的持續(xù)時間進行聯(lián)合搜索。然而,該方法需要進行比GRFT更高維的參數(shù)搜索,這就極大地限制了它的價值。除此之外,在電子波束掃描相控陣列(digital phased antenna array, DPAA)雷達中,高速機動目標可能會在一個較長CPI內(nèi)隨機地進入和/或離開波束的主瓣監(jiān)測區(qū)域[13],這種“跨波束”運動的現(xiàn)象將可能會進一步降低相干累積方法的性能。然而,現(xiàn)有文獻中未報道能有效處理這一現(xiàn)象的方案。

為解決上述問題,本文提出了一種隨機有限集推理動目標軌跡的相干累積檢測法。該方法首先通過距離頻率軸反轉變換,提取描述目標運動特性的信號。隨后,以該信號時頻變換結果中的局部極值構成測量的隨機有限集,并以狀態(tài)空間模型描述各目標的運動,便可將目標運動軌跡的推理轉化為貝葉斯多目標追蹤問題?；谕评淼玫降能壽E,可設計一組濾波器來補償RM和DFM。理論分析和仿真實驗結果表明,本文提出的方法適用于具有復雜甚至未知運動模式的機動目標。

1 信號模型

1.1 接收信號的模型

假設脈沖多普勒雷達發(fā)射線性調(diào)頻(linear frequency modulation, LFM)信號[1,4]如下所示:

(1)

式中:τ表示快時間;Ts是脈沖持續(xù)時間;tm=(m-1)Tr表示慢時間;Tr是脈沖重復間隔(pulse repetition interval,PRI);m=1,2,…,M是脈沖編號且M∈Z+為一個CPI內(nèi)的脈沖數(shù);f0是雷達的工作頻率;B是帶寬,在窄帶條件下B?f0[1];rect(·)是矩形函數(shù),其取值在|x|≤1/2時為1,其余取值為0。

經(jīng)解調(diào)和脈沖壓縮(pulse compression, PC)處理后,K個點狀動目標的回波信號[1]可建模為

(2)

rk(tm)=rk,0+Δrk(tm)

(3)

1.2 速度模糊現(xiàn)象

式(2)中第k個目標的徑向瞬時速度可表示為vk(tm)=(Δrk(tm)-Δrk(tm-1))/Tr。若令fP=1/Tr表示雷達系統(tǒng)的脈沖重復頻率(pulse repetition frequency, PRF),那么當|vk(tm)|大于vmax=cfP/4f0時將發(fā)生速度模糊現(xiàn)象[6,10],此時有

Φ(vk(tm))

(4)

Φ(v)=v-round(v/vb)·vb

(5)

式中:round(·)表示取值到最近整數(shù)的運算。

2 動目標的運動估計

為補償RM和DFM,就需獲知目標的運動軌跡(增量)。為此,本節(jié)介紹一種隨機有限集推理機動目標軌跡的方法。

2.1 觀測信號

本節(jié)利用距離頻率軸反轉變換[14]提取出式(2)中與目標運動相關的部分。首先,對式(2)進行快時間軸傅里葉變換得

(6)

(7)

sself(τ,tm)+scross(τ,tm)

(8)

式中:IFT{·}表示逆傅里葉變換;sself(τ,tm)和scross(τ,tm)分別為

scross(τ,tm)=

(9)

式中:ak表示第k個目標項的幅度;bk,q表示第k個和第q個(q≠k)目標間交叉項幅度?？梢园l(fā)現(xiàn),與目標運動有關的信號就是sself(τ,tm)中τ=0附近的信號,即

ys(tm)=s2(0,tm)=

(10)

2.2 觀測的隨機有限集

對于式(10),可利用短時傅里葉變換(short time Fourier transform, STFT)[15]等時頻分析方法對其進行分析。利用STFT可得如下[15]分布:

TFD(tm,v)=

(11)

式中:速度變量v的取值范圍為[-vmax,vmax];y(tm)=ys(tm)+η(tm)表示信號ys(tm)疊加噪聲η(tm)后的觀測值;W(·)為窗函數(shù),其長度為0

(12)

|TFD(tm,v)|2刻畫了信號的能量分布,故在每一時刻,其局部極值點將與目標的瞬時模糊速度對應?？紤]到噪聲或雜波可能會造成誤檢測,因此如圖1所示,本文利用譜峰檢測技術[1]將|TFD(tm,v)|2每一時刻檢測到的局部極大值的速度值均視為測量。綜上,tm時構造的測量隨機有限集可表示為

(13)

圖1 局部峰值檢測示意圖Fig.1 Illustration of local peak detection

值得注意的是,雖然STFT是最常用的時頻分析手段,但其在面對模糊速度差異較小的回波信號時,可能并不理想。因此,一些時頻分析方法引入了描述頻譜稀疏性的約束,如l1范數(shù)、Cauchy范數(shù)等[16],以期提高局部頻譜的分辨率和能量集中度[17]。然而,這些方法需要進行多次復雜的迭代[16-17],計算復雜度高于STFT。

2.3 狀態(tài)空間模型

由STFT的時頻分析結果不能直接獲取目標瞬時速度的變化規(guī)律[18]。為解決上述問題,介紹一種軌跡推理的方法。首先,為了符號的簡潔,定義與運動有關的相位變量為

(14)

利用多項式預測模型(polynomial prediction model, PPM)[18-19]對式(14)的演化進行建模,可得如下狀態(tài)方程:

xk(tm)=Fxk(tm-1)+ξk(tm-1)

(15)

(16)

式中:gl(l=1,2,…,L)為一有限沖激響應(finite impact respose, FIR)濾波器的系數(shù);I(L-1)×(L-1)和0(L-1)×1分別表示大小為(L-1)×(L-1)的單位陣和(L-1)×1的全0向量。分析表明,FIR的系數(shù)可由L和局部多項式的階數(shù)唯一確定[19-20]。

利用xk(tm)改寫式(12),就能為第k個目標建立如下的觀測模型:

(17)

不失一般性地,可假設各目標的演化過程相互獨立,且產(chǎn)生觀測的過程互不影響[21-22],這樣由式(15)和式(17),可得到如下目標狀態(tài)空間模型:

(18)

式中:N(·;μ,Σ)表示高斯概率密度函數(shù),其均值和方差分別為μ和Σ?；谑?18),并以Zm為測量,便可將軌跡的推理轉化為貝葉斯框架下的多目標追蹤問題。這樣,多假設追蹤(multi hypotheses tracker, MHT)[21]、高斯混合概率假設密度(Gaussian mixture probability hypothesis density, GM-PHD)濾波[22]等方法便可用于軌跡的推理。

2.4 濾波方法

設tm-1時的狀態(tài)概率密度,可借由高斯混合概率密度[22]的形式描述為

P(x|Z1:m-1)=

(19)

Pk(xk(tm)|Z1:m-1)=

N(xk(tm);Fμk(tm-1),FPk(tm-1)FT+Qk)

(20)

(21)

(22)

(23)

在濾波過程中,對于tm-1時就已經(jīng)存在的目標,其存在概率可由如下方式進行更新:

(24)

式中:

0<β<1則表示目標存在概率的衰減因子。在式(20)高斯線性模型的假設下,對于tm-1時存在的目標,其均值及方差的更新過程詳見文獻[21]和文獻[22]。

隨著觀測的到來,假設目標數(shù)將逐漸增加,為提高運行效率,就需要進行“剪枝”[22]?？紤]到噪聲/雜波是隨機的,而目標回波則是連續(xù)且有規(guī)律的,這樣,隨著觀測的到來,真實目標的存在概率將接近1,而由噪聲/雜波產(chǎn)生的虛假目標的存在概率將趨近于0。設閾值δ2>δ1>0,則:① 當目標的存在概率大于δ2時,該假設目標被確認為真實目標;② 當目標的存在概率大于δ1而小于δ2時,該目標處于待定狀態(tài);③ 否則,該目標被認定是來自雜波/噪聲的虛警目標或是正在動態(tài)消失的目標,應將其存在概率置0,并不再對其進行更新。對應地,新生目標的初始存在概率可被假設為δ1和δ2之間的一個常數(shù)。

θi=[θi(t1),θi(t2),…,θi(tM)]T=

(25)

(26)

3 距離聚焦

根據(jù)式(26),第i個假設目標在給定的模糊常數(shù)p下的匹配濾波器[6]可設計為

Ξi(f,tm,p)=

(27)

為方便敘述,不妨設第i個假設目標恰與式(6)中的第k個目標對應,那么將Ξi(f,tm,p)與信號S(f,tm)相乘后,得到解調(diào)信號

(28)

式中:

(29)

(30)

式中:若δrk(tm,p)≈0(m=1,2,…,M)則有且僅有第k個目標回波信號的能量集中在fd=0的附近。

(31)

特別地,若模糊常數(shù)p=pk,則有

Ci(τ,pk)≈

(32)

綜上,本文提出的方法可總結為圖2所示的流程圖。

圖2 所提方法流程圖Fig.2 Flowchart of the proposed method

4 分 析

本文方法在進行反轉變換時,需沿著快時間軸進行M次傅里葉變換和逆傅里葉變換,以及MN次乘運算,計算復雜度為O(2MN(log2N+1))[4,10]。軌跡推理階段,STFT的計算復雜度不大于O(MIlog2I);濾波方法復雜度為O(MKL3)[3,23-24]。每個目標在每個模糊常數(shù)下進行MN次乘運算,N次慢時間維傅里葉變換和1次快時間維的逆傅里葉變換,計算復雜度為O((2P+1)K(MN+MN·log2M+Nlog2N))。在一般應用場景中,有l(wèi)og2M>1和log2N>1,若再取I=M,那么本文方法計算復雜度便為O((2P+1)KMNlog2M+2MNlog2N+M2log2M+MKL3)。

在提取慢時間域描述目標運動軌跡的時頻信號階段,可先使用Radon變換作為非相干累積的手段對軌跡進行粗略搜索,以較高的虛警閾值將回波信號sPC(τ,tm)中不存在目標信號的單元置0。隨后,對處理后的信號再進行距離頻率軸反轉變換,提取描述目標運動的時頻信號。這種兩階段的策略將以增加計算復雜度為代價,提高本文方法在低信噪比情況下的性能[25]。

表1給出了MTD、RFT[6]、R-MHLV[8]、GRFT[4]以及本文方法的計算復雜度。由表1可知,基于全局運動參數(shù)多維度搜索的GRFT的計算復雜度遠高于本文方法。

表1 不同方法計算復雜度的對比

5 數(shù)值仿真

雷達系統(tǒng)仿真參數(shù)如表2所示。設雷達所需探測的速度范圍為[-1 000 m/s,+1 000 m/s],則雷達參數(shù)的模糊整數(shù)搜索范圍為[-8,8]。數(shù)值仿真中,距離門的編號nr∈N和徑向距離r∈R+之間的關系為

nr=round((r-rmin)/ρr)

(33)

式中:rmin表示檢測場景中所考慮的最小徑向距離,本文中取5.5×103m;ρr=c/2fs是距離門的寬度。

表2 雷達系統(tǒng)參數(shù)

5.1 有效性驗證

本仿真實驗中,兩個目標的運動分別為

(34)

(35)

式中:β=3.0。本組實驗中,脈沖壓縮后目標1和目標2的信噪比均為5 dB[4]。由式(5)可知,將目標1和目標2的速度映射到區(qū)間[-vmax,vmax]內(nèi),分別得到-55 m/s和31.38 m/s。

綜合考慮復雜度與濾波的性能,時頻分析方法選用STFT,取長度為51的高斯窗,頻點數(shù)I為512。在式(15)的狀態(tài)模型中,FIR濾波器的階數(shù)為2,抽頭數(shù)L為3。假設所有目標具有相同的過程噪聲協(xié)方差,且協(xié)方差對角線上的元素為[1e-5, 1e-5, 1e-5];觀測噪聲的方差由文獻[5]中的方法估計。多目標追蹤方法為MHT;PS=0.95;δ1=0.1,δ2=0.5,δd=0.05,β=0.85;新生初始存在概率為0.3。

圖3給出了本文方法對回波信號的處理結果。圖3(a)為脈沖壓縮后的結果,由圖3(a)可以看到目標有明顯的距離徙動。圖3(b)為提取的時頻信號的STFT譜。在圖3(c)中,本文方法很好地估計了目標速度的變化。圖3(d)為方法執(zhí)行過程中假設的目標生成數(shù)和確認的目標數(shù)?？梢钥吹?每次迭代過程中新生假設目標數(shù)量是隨機的。圖3(e)～圖3(f)分別為兩個目標具有最大響應值的距離聚焦輸出信號,圖中將輸出信號模的最大值進行了歸一化。上述結果表明,本文方法能應對回波信號中多個機動目標距離和速度變化軌跡相互交叉的情況。值得注意的是,上述仿真過程沒有復雜的二維及二維以上的運動參數(shù)搜索,即能以較低的計算復雜度,實現(xiàn)對具有復雜運動形式的機動目標的相干累積。

圖3 本文方法處理結果Fig.3 Processing results of the proposed method

接下來將繼續(xù)研究速度估計的相對均方根誤差(relative root mean square error, R-RMSE)和檢測率之間的關系。速度的R-RMSE定義為

R-RMSEk=

(36)

實驗中,使用單元平均—恒虛警率(cell averaging constant false alarm rate, CA-CFAR)[1]檢測目標。在檢測時,將進行如下二元檢測:

(37)

式中:Γi表示|Ci(r)|的最大值;ηT,i表示第i個假設目標的自適應檢測閾值,其值可在給定虛警率下由式(31)中的|Ci(r)|的統(tǒng)計量進行估計。具體而言,將Γi附近若干個單元作為保護單元(本文中取20個左右距離單元),將其他所有距離單元作為參考單元[1,6],并計算如下統(tǒng)計量:

(38)

(39)

本實驗中,設置虛警率為Pfa=10-4,PC后信噪比變化范圍為-20 dB到20 dB。圖4為目標2的運動以式(35)中的方程進行描述且取不同值時的探測結果?？梢钥闯?隨著PC后信噪比的逐漸升高,速度相對誤差逐漸減小且探測率逐漸增加;β的絕對值越小,即DFM的變化范圍越小,方法的表現(xiàn)越好。

圖4 目標2速度估計的R-RMSE及探測率間的關系Fig.4 Relationship between R-RMSE velocity estimation and detection probability of target 2

5.2 方法對比

5.2.1 輸出結果對比

本節(jié)給出MTD、RFT、R-MHLV、GRFT以及STGRFT的輸出結果。為直觀地比較各方法的性能,GRFT會先分別搜索目標的加速度、加加速度等參數(shù)。隨后,再利用估計的參數(shù)補償?shù)裟繕?的DFM,再由RFT計算被補償后的信號的距離—速度圖。R-MHLV和STGRFT使用和GRFT相同的策略。除此之外,STGRFT在進行參數(shù)搜索時,假定每個目標都包含兩個運動階段,且第二運動階段起始脈沖的搜索范圍均為[230,280]。為直觀地進行比較,本節(jié)方法也先利用推理的軌跡補償RM和DFM,再對補償后的信號使用RFT計算其距離—速度圖。實驗中,目標2的信噪比為-2 dB,其余設置同第5.1節(jié)。

圖5給出了各方法的輸出結果對比。在MTD(圖5(a))和RFT(圖5(b))的結果中,目標的能量完全淹沒在噪聲中。圖5(c)和圖5(d)中R-MHLV和GRFT只能對目標2的進行能量集中度較低的相干累計,這是因為R-MHLV和GRFT中的全局三階多項式運動模型無法很好地描述目標2的運動。圖5(e)為STGRFT的結果,可以發(fā)現(xiàn)其優(yōu)于R-MHLV和GRFT,但仍不太理想,這是因為STGRFT所假設的兩分段二階運動模型同樣與目標2的運動模式存在一定程度上的失配。相比其他方法(圖5(a)至圖5(d)),本文方法的距離—速度域輸出結果(如圖5(f)所示)更優(yōu)。

5.2.2 探測性能對比

本節(jié)比較各方法的檢測性能。在仿真中,使用CA-CFAR檢測器檢測目標,虛警率為Pfa=10-4。信噪比的變化范圍為-30 dB到20 dB。在每種實驗條件下,都進行500次蒙特卡羅仿真。其余設置同第6.1節(jié)。圖6(a)展示了各方法在不同信噪比下對目標1的探測率。對于目標1,本文方法略遜于STGRFT方法,但優(yōu)于其余方法。圖6(b)給出了各方法對目標2的探測率。目標2的運動模型更為復雜,因此本文方法具有更加明顯的優(yōu)勢。

圖6 不同信噪比下的目標探測率Fig.6 Detection probability via different input signal to noise ratio

5.3 復雜情況的處理結果

考慮動態(tài)出現(xiàn)和/或消失情況,設計出的二維匹配濾波器為

(40)

本節(jié)分析本文方法對動態(tài)目標的有效性。設某CPI內(nèi),目標1自120個脈沖開始出現(xiàn),并持續(xù)至CPI結束,其速度如圖7(a)所示,但具體表達式未知,目標1的初始徑向距離為6.8×103m;目標2從第1個脈沖時起就存在,但在第400個脈沖時消失。其運動表達式為

Δr2(tm)=6.5×103+3.0cos 2.5πtm+340tm,

0

(41)

其余設置同第5.1節(jié)。從圖7(b)PC后的結果可以看到,目標存在動態(tài)的出現(xiàn)/消失情況。圖7(c)為本文方法對目標瞬時速度的估計結果;圖7(d)為本文方法對目標數(shù)的估計結果。圖7(e)為不同信噪比下目標2的探測性能,蒙特卡羅仿真設定同第5.2節(jié)。以上仿真實驗結果表明,對于存在動態(tài)出現(xiàn)和/或消失現(xiàn)象的目標,本文方法能取得較優(yōu)的結果。

圖7 動態(tài)目標的處理結果Fig.7 Processing results of the dynamic targets

6 結 論

針對相干處理間隔內(nèi),由RM和DFM導致的相干累積性能下降問題,本文提出了一種隨機有限集推理動目標軌跡的相干累積檢測法。該方法首先通過距離頻率軸反轉變換,提取出可描述目標運動特性的時頻信號。然后,通過時頻分析方法將該信號轉換為觀測的隨機有限集合,結合描述目標運動軌跡演化的狀態(tài)空間模型,利用隨機集方法實現(xiàn)運動軌跡的估計?；诘玫降倪\動軌跡,設計出合適的匹配濾波器來補償目標的RM和DFM。理論分析和仿真結果表明,本文提出的方法,能有效地處理具有復雜運動形式的機動目標。