宋 聰，黎 楊*，2

1.武漢工程大學(xué)電氣信息學(xué)院，湖北 武漢430205；2.湖北省視頻圖像與高清投影工程技術(shù)研究中心，湖北 武漢430205

近年來，基于第三方非合作輻射源的外輻射源雷達(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á)中，天線接收微弱的目標(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ì)問題是一個(gè)極強(qiáng)干擾、極低信噪比情況下的DOA估計(jì)問題。

近幾十年來，出現(xiàn)了很多的算法用于估計(jì)信號(hào)DOA，例如：多重信號(hào)分類（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)景的正交線性預(yù)測(cè)傳播算子算法（orthogonal propagator method-linear prediction，OPM-LP）方法，該方法利用線性預(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天線陣列測(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ì)回波方向，但是該方法依賴于特定的Adcock天線陣列結(jié)構(gòu)。文獻(xiàn)［16］提出了一種適用于外輻射源雷達(dá)的比幅測(cè)角方法，該方法利用八單元均勻圓陣天線結(jié)合方向圖綜合技術(shù)形成覆蓋全空域的18個(gè)波束，用以對(duì)目標(biāo)進(jìn)行掃描。文獻(xiàn)［17］使用了2個(gè)監(jiān)視天線，通過計(jì)算距離-多普勒?qǐng)D上的相位差來估計(jì)目標(biāo)散射波。盡管上述方法在應(yīng)用中能獲得一些強(qiáng)目標(biāo)的方向，但是依賴于特定的天線配置，并且其精度和分辨率均很低。

針對(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)與天線各個(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)，最后通過頻率掃描的方式，用MUSIC算法進(jìn)行空間譜估計(jì)得到二維DOA-多普勒?qǐng)D，最終可獲得散射波的DOA，及其對(duì)應(yīng)的多普勒頻率。

1 問題描述

1.1 FM廣播信號(hào)

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］：

1.2 陣列信號(hào)模型

考慮一個(gè)均勻直線陣列（uniform linear array，ULA），陣元為全向性天線。假設(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)稱為“實(shí)際陣列信號(hào)”。

1.3 MUSIC空間譜

MUSIC空間譜可表示為：

其中譜的峰值對(duì)應(yīng)的角度θ為波達(dá)方向。

1.4 直達(dá)波參考信號(hào)獲取

波束形成技術(shù)用于空間濾波。參考通道中的直達(dá)波信號(hào)由yR(k)=wHx(k)獲得，其中wH是波束形成器指向直達(dá)波方向時(shí)的權(quán)矢量，可近似認(rèn)為經(jīng)過波束形成后yR(k)=s(k)。

1.5 模糊函數(shù)

雷達(dá)信號(hào)的模糊函數(shù)定義為：

其中τ表示時(shí)延，ν表示多普勒頻移。模糊函數(shù)的一個(gè)性質(zhì)是在A(0,0)處響應(yīng)最大。

2 新方法介紹

2.1 虛擬陣列信號(hào)構(gòu)建

參考信號(hào)yR(k)與第m個(gè)天線的接收信號(hào)xm(k)的互模糊函數(shù)表達(dá)式為：

其中m=1,2…M。式（8）可進(jìn)一步表示為矩陣形式的陣列互模糊函數(shù)F(τ,ν)：

矩陣F(τ,ν)具有與實(shí)際陣列信號(hào)相似的形式，在本文中稱為“虛擬陣列信號(hào)”。虛擬陣列信號(hào)F(τ,ν)具有3個(gè)虛擬信號(hào)分量：具有導(dǎo)向矢量aS的信號(hào)分量，稱為“虛擬直達(dá)波信號(hào)”；具有導(dǎo)向矢量ai的信號(hào)分量，稱為“虛擬散射波信號(hào)”；含n(k)的分量稱為“虛擬噪聲信號(hào)”。

假設(shè)一個(gè)期望散射波信號(hào)的導(dǎo)向矢量為a1，時(shí)差為τ1，多普勒頻移為ν1，且yR(k)=s(k)，則虛擬陣列信號(hào)F(τ1,ν1)可以表示為：

2.2 期望散射波信號(hào)的增強(qiáng)

因此，虛擬陣列信號(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)過式（11）處理后被增強(qiáng)了。從式（12）和式（13）可以看出δ和ε與模糊函數(shù)圖形上的峰值到底平面的差有關(guān)，γ與處理時(shí)長(zhǎng)有關(guān)，因此可以通過增加處理時(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ν矩陣被稱為“多個(gè)采樣點(diǎn)的虛擬陣列信號(hào)”。

2.3 散射波DOA估計(jì)

采用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-多普勒頻移圖。

2.4 實(shí)現(xiàn)步驟

算法流程如圖1所示，主要步驟總結(jié)如下：

步驟1：估計(jì)直達(dá)波信號(hào)的DOA：θS。

步驟2：對(duì)直達(dá)波方向進(jìn)行波束形成，獲得直達(dá)波信號(hào)yR。

步驟3：通過式（8）計(jì)算yR和每個(gè)天線接收信號(hào)的互模糊函數(shù)，然后將它們抽取組合構(gòu)造單個(gè) 采樣點(diǎn)的虛擬陣列信號(hào)F(τ,ν)，接著將多普勒頻移設(shè)定的區(qū)間范圍均勻劃分為L(zhǎng)ν個(gè)點(diǎn)v=linspace(-νmax,νmax,Lν)來構(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

3 仿真分析

3.1 FM輻射源信號(hào)

以采樣頻率為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)過下變頻后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

3.2 陣列信號(hào)參數(shù)

接收陣列為M=16的均勻直線陣，假設(shè)陣列流型已知，并忽略耦合效應(yīng)，2 s長(zhǎng)的信號(hào)經(jīng)天線陣接收后進(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

3.3 散射波DOA估計(jì)

在第一個(gè)仿真例子中，使用本文提出的方法，即步驟1～4，計(jì)算DOA-多普勒?qǐng)D。如圖3所示，圖3（a）結(jié)果表明：1）在DOA為30°時(shí)出現(xiàn)一條“脊線”，貫穿整個(gè)多普勒頻移區(qū)間。2）“脊線”的峰值表示多普勒頻移為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

3.4 性能對(duì)比

為了與其他方法進(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

4 結(jié) 論

本文提出了一種距離多普勒域陣列信號(hào)處理方法，用于估計(jì)FM外輻射源雷達(dá)散射波信號(hào)的DOA。通過計(jì)算直達(dá)波信號(hào)與天線各個(gè)陣元接收信號(hào)的互模糊函數(shù)，然后抽取數(shù)據(jù)組成多個(gè)采樣點(diǎn)的虛擬陣列信號(hào)。理論分析表明期望散射波信號(hào)在虛擬陣列信號(hào)中被顯著增強(qiáng)，接下來可使用MUSIC算法估計(jì)散射波信號(hào)的DOA。仿真結(jié)果表明，本文所提出的方法能夠在強(qiáng)直達(dá)波干擾情況下成功估計(jì)極低信噪比散射波的DOA，及其對(duì)應(yīng)的多普勒頻率。