宋 聰,黎 楊*,2
1.武漢工程大學(xué)電氣信息學(xué)院,湖北 武漢430205;2.湖北省視頻圖像與高清投影工程技術(shù)研究中心,湖北 武漢430205
近年來(lái),基于第三方非合作輻射源的外輻射源雷達(dá)系統(tǒng),由于具有潛在的隱蔽、反隱身目標(biāo)和抗干擾等特性,受到了特別關(guān)注[1-2]。外輻射源雷達(dá)利用第三方非協(xié)作輻射源進(jìn)行目標(biāo)探測(cè),由于其本身不發(fā)射信號(hào),因此相對(duì)于傳統(tǒng)的有源雷達(dá),具有很多的優(yōu)勢(shì),如成本低、體積小、對(duì)電磁環(huán)境無(wú)污染和具有探測(cè)隱身目標(biāo)的潛力等[3]。
目前,多種外輻射源信號(hào)被用做無(wú)源探測(cè),如調(diào)頻(frequency modulation,F(xiàn)M)廣播信號(hào)[4]、移動(dòng)通信基站[6]、數(shù)字/模擬電視信號(hào)[5]等。FM廣播信號(hào)具有站點(diǎn)多、覆蓋范圍廣、組網(wǎng)潛力大、模糊函數(shù)距離分辨率高等優(yōu)點(diǎn),被認(rèn)為是最適合用于無(wú)源探測(cè)的輻射源之一[7]。
在外輻射源雷達(dá)中,天線(xiàn)接收微弱的目標(biāo)散射波信號(hào),同時(shí)還會(huì)接收到很強(qiáng)的直達(dá)波、多徑和雜波干擾。由于散射波的信噪比(signal-to-noise ratio,SNR)通常極低(接收到散射波的SNR典型值約為-40 dB左右),同時(shí)還存在很強(qiáng)的直達(dá)波干擾(直達(dá)波與散射波的功率之比的典型值大于90 dB)[8],散射波的波達(dá)方向(direction-of-arrival,DOA)估計(jì)問(wèn)題是一個(gè)極強(qiáng)干擾、極低信噪比情況下的DOA估計(jì)問(wèn)題。
近幾十年來(lái),出現(xiàn)了很多的算法用于估計(jì)信號(hào)DOA,例如:多重信號(hào)分類(lèi)(multiple signal classification,MUSIC)[9]、旋轉(zhuǎn)不變子空 間(estimation of signal parameters via rotational invariance technique,ESPRIT)[10],以及各種特殊應(yīng)用場(chǎng)景的方法[11]等。然而,現(xiàn)有方法很少能適用于微弱信號(hào)DOA估計(jì)。在微弱信號(hào)DOA估計(jì)方面,文獻(xiàn)[12]提出了一種基于軟稀疏表示的DOA估計(jì)方法,該方法基于稀疏性的前提下,應(yīng)用迭代加權(quán)最小方差法求解軟稀疏解,最后估計(jì)目標(biāo)方位。文獻(xiàn)[13]提出了一種協(xié)方差矩陣的重構(gòu)方法,該方法能夠明顯地提高協(xié)方差矩陣的信噪比,將新的協(xié)方差矩陣應(yīng)用到最小方差無(wú)畸變響應(yīng)算法進(jìn)行DOA估計(jì),改善了低信噪比條件下的DOA估計(jì)性能。文獻(xiàn)[14]提出了一種適用于低信噪比場(chǎng)景的正交線(xiàn)性預(yù)測(cè)傳播算子算法(orthogonal propagator method-linear prediction,OPM-LP)方法,該方法利用線(xiàn)性預(yù)測(cè)原理,重構(gòu)降噪的數(shù)據(jù)協(xié)方差矩陣,然后使用正交傳播算法進(jìn)行DOA估計(jì),但是該方法不適用于-40 dB極低信噪比情況。盡管上述方法一定程度上能夠估計(jì)微弱信號(hào)的DOA,但并不適用于本文信噪比極低,同時(shí)還存在很強(qiáng)的直達(dá)波干擾的應(yīng)用場(chǎng)景。
在外輻射源雷達(dá)目標(biāo)散射波估計(jì)方面,文獻(xiàn)[15]提出了一種利用四元Adcock天線(xiàn)陣列測(cè)量目標(biāo)散射波DOA估計(jì)方法,該方法預(yù)先計(jì)算每一個(gè)通道的接收信號(hào)與參考信號(hào)的互模糊函數(shù),在距離多普勒?qǐng)D上消除了直達(dá)波、多徑和雜波干擾之后,再利用成對(duì)監(jiān)視信號(hào)之間的角度關(guān)系估計(jì)回波方向,但是該方法依賴(lài)于特定的Adcock天線(xiàn)陣列結(jié)構(gòu)。文獻(xiàn)[16]提出了一種適用于外輻射源雷達(dá)的比幅測(cè)角方法,該方法利用八單元均勻圓陣天線(xiàn)結(jié)合方向圖綜合技術(shù)形成覆蓋全空域的18個(gè)波束,用以對(duì)目標(biāo)進(jìn)行掃描。文獻(xiàn)[17]使用了2個(gè)監(jiān)視天線(xiàn),通過(guò)計(jì)算距離-多普勒?qǐng)D上的相位差來(lái)估計(jì)目標(biāo)散射波。盡管上述方法在應(yīng)用中能獲得一些強(qiáng)目標(biāo)的方向,但是依賴(lài)于特定的天線(xiàn)配置,并且其精度和分辨率均很低。
針對(duì)目標(biāo)散射波的DOA估計(jì)難題,本文提出一種距離多普勒域陣列信號(hào)處理(range-doppler domain array signal processing,RD-ASP)方法,在估計(jì)其散射波的DOA之前增強(qiáng)期望散射波信號(hào)?!霸鰪?qiáng)”一詞意味著分別極大地提高了期望散射波信號(hào)的功率與直達(dá)波信號(hào)的功率、與干擾散射波信號(hào)的功率、以及與噪聲的功率之比。該方法首先計(jì)算直達(dá)波信號(hào)與天線(xiàn)各個(gè)陣元接收信號(hào)之間的互模糊函數(shù),然后在互模糊函數(shù)上抽取特定時(shí)差和多普勒的數(shù)據(jù)構(gòu)建單個(gè)采樣點(diǎn)虛擬陣列信號(hào),再將單個(gè)采樣點(diǎn)虛擬陣列信號(hào)擴(kuò)展為多個(gè)采樣點(diǎn)虛擬陣列信號(hào),最后通過(guò)頻率掃描的方式,用MUSIC算法進(jìn)行空間譜估計(jì)得到二維DOA-多普勒?qǐng)D,最終可獲得散射波的DOA,及其對(duì)應(yīng)的多普勒頻率。
FM廣播信號(hào)一般由主信道信號(hào)、副信道信號(hào)和導(dǎo)頻信號(hào)組成的立體聲復(fù)合信號(hào)。主信道信號(hào)是左聲道信號(hào)和右聲道信號(hào)的和信號(hào),副信道信號(hào)是左聲道和右聲道信號(hào)的差信號(hào),導(dǎo)頻信號(hào)是為了接收而傳送的輔助信號(hào)。t時(shí)刻的立體聲復(fù)合信號(hào)m(t)可表示為[18]:
考慮一個(gè)均勻直線(xiàn)陣列(uniform linear array,ULA),陣元為全向性天線(xiàn)。假設(shè)無(wú)同頻干擾信號(hào),接收信號(hào)由一個(gè)直達(dá)波和多個(gè)散射波組成。接收到的陣列信號(hào)x (k)可表示為:
其中k是采樣點(diǎn)數(shù);i=S表示直達(dá)波,i≥1表示散射波;ai,i=S,1,…,L是第i個(gè)信號(hào)的導(dǎo)向矢量;s(k)是時(shí)域中的直達(dá)波;ηi,τi和νi是第i個(gè)信號(hào)(相對(duì)于直達(dá)波)的幅度衰減,時(shí)差和多普勒頻移;FS是采樣率;n(k)表示噪聲。x(k),ai和n(k)的維數(shù)為M×1維,其中M表示陣元數(shù)。假定噪聲是高斯白噪聲。在本文中,陣列信號(hào)x(k)稱(chēng)為“實(shí)際陣列信號(hào)”。
MUSIC空間譜可表示為:
其中譜的峰值對(duì)應(yīng)的角度θ為波達(dá)方向。
波束形成技術(shù)用于空間濾波。參考通道中的直達(dá)波信號(hào)由yR(k)=wHx(k)獲得,其中wH是波束形成器指向直達(dá)波方向時(shí)的權(quán)矢量,可近似認(rèn)為經(jīng)過(guò)波束形成后yR(k)=s(k)。
雷達(dá)信號(hào)的模糊函數(shù)定義為:
其中τ表示時(shí)延,ν表示多普勒頻移。模糊函數(shù)的一個(gè)性質(zhì)是在A(0,0)處響應(yīng)最大。
參考信號(hào)yR(k)與第m個(gè)天線(xiàn)的接收信號(hào)xm(k)的互模糊函數(shù)表達(dá)式為:
其中m=1,2…M。式(8)可進(jìn)一步表示為矩陣形式的陣列互模糊函數(shù)F(τ,ν):
矩陣F(τ,ν)具有與實(shí)際陣列信號(hào)相似的形式,在本文中稱(chēng)為“虛擬陣列信號(hào)”。虛擬陣列信號(hào)F(τ,ν)具有3個(gè)虛擬信號(hào)分量:具有導(dǎo)向矢量aS的信號(hào)分量,稱(chēng)為“虛擬直達(dá)波信號(hào)”;具有導(dǎo)向矢量ai的信號(hào)分量,稱(chēng)為“虛擬散射波信號(hào)”;含n(k)的分量稱(chēng)為“虛擬噪聲信號(hào)”。
假設(shè)一個(gè)期望散射波信號(hào)的導(dǎo)向矢量為a1,時(shí)差為τ1,多普勒頻移為ν1,且yR(k)=s(k),則虛擬陣列信號(hào)F(τ1,ν1)可以表示為:
因此,虛擬陣列信號(hào)δ的值比實(shí)際陣列信號(hào)大。
將虛擬期望散射波信號(hào)與虛擬干擾散射波信號(hào)的功率的比值定義為ε,則虛擬陣列信號(hào)F(τ1,ν1)的ε為:
因此,虛擬陣列信號(hào)F(τ1,ν1)的ε比實(shí)際陣列信號(hào)x(k)大。
對(duì)于噪聲,將虛擬期望散射波信號(hào)與虛擬噪聲信號(hào)的功率比定義為γ,使用柯西-施瓦茨不等式可證明虛擬陣列信號(hào)F(τ1,ν1)的γ比實(shí)際陣列信號(hào)x(k)大。
其中PS是直達(dá)波信號(hào)的功率,PN是噪聲的功率。
由式(12)~式(14)可知,導(dǎo)向矢量為a1的虛擬期望散射波信號(hào)經(jīng)過(guò)式(11)處理后被增強(qiáng)了。從式(12)和式(13)可以看出δ和ε與模糊函數(shù)圖形上的峰值到底平面的差有關(guān),γ與處理時(shí)長(zhǎng)有關(guān),因此可以通過(guò)增加處理時(shí)間,進(jìn)一步增強(qiáng)期望散射波信號(hào)。
然而式(11)中的虛擬陣列信號(hào)F(τ1,ν1)只有單個(gè)采樣點(diǎn),很難直接使用傳統(tǒng)的DOA估計(jì)方法。下面將單個(gè)采樣點(diǎn)的虛擬陣列信號(hào)擴(kuò)展為多個(gè)采樣點(diǎn)的虛擬陣列信號(hào)。將多普勒頻移設(shè)定的區(qū)間范圍(-νmax,νmax)均勻劃分為L(zhǎng)ν個(gè)點(diǎn)ν=linspace(-νmax,νmax,Lν),對(duì)于頻率νFix∈(-νmax,νmax),可 以 獲 得Lt個(gè) 矩 陣:F[t(1),νFix],F[t(2),νFix],…,F[t(Lt),νFix],每個(gè)矩陣的維數(shù)為M×1,構(gòu)建一個(gè)M×Lt維矩陣Gν:
Gν矩陣被稱(chēng)為“多個(gè)采樣點(diǎn)的虛擬陣列信號(hào)”。
采用MUSIC算法進(jìn)行空間譜估計(jì)。對(duì)于多個(gè)采樣點(diǎn)的虛擬陣列信號(hào)Gν,其采樣數(shù)據(jù)協(xié)方差矩陣可表示為:
其中:ri和λi是Rν的特征值和對(duì)應(yīng)的特征向量;特征值按降序排列,即r1>…>rd>>rd+1=…=rM;Rν的前d個(gè)大特征值對(duì)應(yīng)的特征向量組成虛擬的信號(hào)子空間ES=[ λ1,…,λd];剩下的M-d個(gè)小特征值對(duì)應(yīng)的特征向量組成虛擬的噪聲子空間EN=[λd+1,…,λM]。
最后用MUSIC算法進(jìn)行空間譜估計(jì):
其中P (θ,νFix)的峰值表示具有νFix頻移信號(hào)的DOA。如果掃描所有多普勒頻移區(qū)間,則可以獲得所有散射波和直達(dá)波信號(hào)的DOA,以及它們對(duì)應(yīng)的多普勒頻移,最后可以組合得到DOA-多普勒頻移圖。
算法流程如圖1所示,主要步驟總結(jié)如下:
步驟1:估計(jì)直達(dá)波信號(hào)的DOA:θS。
步驟2:對(duì)直達(dá)波方向進(jìn)行波束形成,獲得直達(dá)波信號(hào)yR。
步驟3:通過(guò)式(8)計(jì)算yR和每個(gè)天線(xiàn)接收信號(hào)的互模糊函數(shù),然后將它們抽取組合構(gòu)造單個(gè) 采樣點(diǎn)的虛擬陣列信號(hào)F(τ,ν),接著將多普勒頻移設(shè)定的區(qū)間范圍均勻劃分為L(zhǎng)ν個(gè)點(diǎn)v=linspace(-νmax,νmax,Lν)來(lái)構(gòu)建多個(gè)采樣點(diǎn)的虛擬陣列信號(hào)Gν。
步驟4:對(duì)多個(gè)采樣點(diǎn)的虛擬陣列信號(hào)Gν應(yīng)用MUSIC算法進(jìn)行空間譜估計(jì)。
圖1 算法流程圖Fig.1 Flowchart of algorithm
以采樣頻率為44.1 kHz從一段語(yǔ)音信號(hào)中截取時(shí)間為2 s的語(yǔ)音信號(hào),將此語(yǔ)音信號(hào)作為基帶信號(hào)進(jìn)行FM調(diào)制,將此調(diào)制信號(hào)進(jìn)行下變頻,下變頻后的采樣率為300 kHz。圖2(a)和圖2(b)分別為經(jīng)過(guò)下變頻后2 s時(shí)長(zhǎng)的FM信號(hào)的功率譜圖和模糊函數(shù)圖。其中P表示FM信號(hào)功率,f表示FM信號(hào)頻率,SNR的值為κ,R表示距離差,ν表示多普勒頻移,χ=20log10||A()
τ,ν。
圖2 下變頻后的FM信號(hào):(a)功率譜圖,(b)模糊函數(shù)圖Fig.2 FM signal after down-conversion:(a)power spectrum chart,(b)ambiguity function diagram
接收陣列為M=16的均勻直線(xiàn)陣,假設(shè)陣列流型已知,并忽略耦合效應(yīng),2 s長(zhǎng)的信號(hào)經(jīng)天線(xiàn)陣接收后進(jìn)行下變頻,得到接收陣列信號(hào)。表1給出了直達(dá)波信號(hào)和5個(gè)散射波信號(hào)的仿真參數(shù)。用于掃描的多普勒頻移區(qū)間從-40 Hz變化到40 Hz,步長(zhǎng)為0.25 Hz,距離差區(qū)間從-100 km變化到100 km,步長(zhǎng)為1 km。
表1 仿真參數(shù)Tab.1 Simulation parameters
在第一個(gè)仿真例子中,使用本文提出的方法,即步驟1~4,計(jì)算DOA-多普勒?qǐng)D。如圖3所示,圖3(a)結(jié)果表明:1)在DOA為30°時(shí)出現(xiàn)一條“脊線(xiàn)”,貫穿整個(gè)多普勒頻移區(qū)間。2)“脊線(xiàn)”的峰值表示多普勒頻移為0 Hz和DOA為30°的直達(dá)波信號(hào)。3)在圖3(a)中可觀測(cè)到5個(gè)以上具有不同的DOA和對(duì)應(yīng)的多普勒頻移的峰。
圖3(b)給出的是移除直達(dá)波信號(hào)所在區(qū)域后的DOA-多普勒頻移圖,可清晰地觀察到5個(gè)峰值,分別對(duì)應(yīng)于5個(gè)散射波信號(hào)。其中SNR的值為-40 dB的散射波信號(hào)3的DOA也被成功檢測(cè)到。
圖3 DOA-多普勒頻移圖:(a)所有角度,(b)部分角度Fig.3 DOA-Doppler-shift spectra:(a)all directions,(b)some directions
為了與其他方法進(jìn)行性能比較,將二維的DOA-多普勒?qǐng)D沿著多普勒維投影成一維的空間譜。其中ξ=10log10[P(θ,νFix)]。
圖4(a)給出的是低信噪比情況下的3種算法性能比較圖,從圖4(a)中可以看到MUSIC算法、OPM-LP算法、RD-ASP算法在低信噪比情況下都能夠分辨直達(dá)波和5個(gè)散射波信號(hào)的DOA。
圖4(b)給出的是極低信噪比情況下的3種算法性能比較圖,從圖4(b)中可以觀察到OPM-LP和MUSIC均不能分辨極低SNR的散射波,而本文提出的RD-ASP方法能夠清晰地分辨5個(gè)散射波的DOA。
圖4 3種算法的性能比較:(a)低信噪比,(b)極低信噪比Fig.4 Performance comparison of threealgorithms:(a)low signal-to-noiseratio,(b)extremely low signal-to-noiseratio
本文提出了一種距離多普勒域陣列信號(hào)處理方法,用于估計(jì)FM外輻射源雷達(dá)散射波信號(hào)的DOA。通過(guò)計(jì)算直達(dá)波信號(hào)與天線(xiàn)各個(gè)陣元接收信號(hào)的互模糊函數(shù),然后抽取數(shù)據(jù)組成多個(gè)采樣點(diǎn)的虛擬陣列信號(hào)。理論分析表明期望散射波信號(hào)在虛擬陣列信號(hào)中被顯著增強(qiáng),接下來(lái)可使用MUSIC算法估計(jì)散射波信號(hào)的DOA。仿真結(jié)果表明,本文所提出的方法能夠在強(qiáng)直達(dá)波干擾情況下成功估計(jì)極低信噪比散射波的DOA,及其對(duì)應(yīng)的多普勒頻率。