李鵬飛 張 旻 鐘子發(fā) 羅 爭

(合肥電子工程學(xué)院309研究室 合肥 230037)

(安徽省電子制約技術(shù)重點實驗室 合肥 230037)

1 引言

在現(xiàn)代電子戰(zhàn)信號環(huán)境中，由于信號密度大，信號形勢復(fù)雜，使得電子偵察中的信號處理任務(wù)越來越困難。實際的測向裝備，特別是無線電通信測向面臨的是一個很寬的頻段(如頻率范圍在30～3000 MHz)。為了提高測向的搜索速度，測向接收機(jī)的瞬時帶寬朝著寬帶的趨勢發(fā)展，可以達(dá)到幾十MHz，甚至上百M(fèi)Hz。在如此寬的瞬時帶寬內(nèi)，必然包含了很多不同頻率的多種信號，如何對這些信號的方向進(jìn)行實時、高精度的估計是一個巨大的挑戰(zhàn)。直接利用傳統(tǒng)的高分辨率測向算法如多重信號分類法[1](MUltiple Signal Classification, MUSIC)將會產(chǎn)生嚴(yán)重的估計誤差，且無法獲得信號的準(zhǔn)確頻率[2]。并且 MUSIC算法不能直接對相干信號進(jìn)行處理，利用空間平滑的方法解相干將損失一定的陣列孔徑。利用稀疏分解的方法進(jìn)行 DOA估計具有很高的分辨率，不需要進(jìn)行任何預(yù)處理，可以直接應(yīng)用到相干信號上來，得到了國內(nèi)外學(xué)者的廣泛關(guān)注[3-6]。但是在進(jìn)行寬頻度DOA估計時字典的構(gòu)造是一個難題，Malioutov等[5]提出了根據(jù)頻率和方位角的組合構(gòu)建聯(lián)合字典的方法。隨著頻段的展寬，用于稀疏表示的過完備字典的長度隨之增加。一方面，稀疏分解的計算量將隨著字典的長度指數(shù)倍的增長[7]；另一方面，在寬頻段的范圍內(nèi)很難保證過完備字典沒有相同的元素，當(dāng)存在相同或者相近的元素時，字典的相關(guān)系數(shù)接近于 1，稀疏分解的精度急劇下降[8]；另外，該方法能分辨的信號數(shù)目受到陣元數(shù)目的限制。

針對以上問題，本文引入空間頻率的概念，將頻率和方位角的2維空間降為空間角的1維空間。利用空間角構(gòu)建稀疏表示的過完備字典，字典長度大大降低，減少了稀疏分解的計算量。為了提高算法在低信噪比的性能，摒棄了常用的在時域或頻域的稀疏模型，對頻域峰值協(xié)方差矩陣進(jìn)行特征分解，利用主特征向量建立稀疏模型進(jìn)行 DOA估計。算法從頻域?qū)π盘柕念l域分布進(jìn)行了預(yù)估計和分離，估計的信號數(shù)目能遠(yuǎn)大于陣元數(shù)目。

2 數(shù)學(xué)模型

2.1 陣列信號的數(shù)學(xué)模型

假設(shè)p個遠(yuǎn)場信號入射到陣元間距為r的M元均勻線陣上。第i個信號的方位角為θi，信號頻率為fi。M個陣元接收到的信號寫成一個矢量，可表示為

其中A(f,θ)是陣列流型矩陣，S(t)= [s(f1),s(f2),…,s(fp)]T為空間信號矢量，噪聲矩陣N(t)是獨立于信號的白高斯過程。陣列流型為

其中第i個導(dǎo)向矢量為

c為電磁波傳播的速度。

對于具有相同頻率的多個窄帶信號fi=f2=…fp, 信號的陣列模型與式(1)描述相同。常規(guī)的測向方法多建立在測向接收機(jī)是窄帶接收機(jī)的基礎(chǔ)上，接收的信號可以看作是在中心頻率附近很窄的頻率范圍內(nèi)，因此這些算法可以直接運(yùn)用。而寬帶測向接收機(jī)具有很寬的瞬時測向帶寬，在一幀的采樣數(shù)據(jù)內(nèi)，可能包含多個窄帶信號或者寬帶信號。MUSIC等算法處理起來難度就大大增加了，需要將整個寬頻段劃分為若干個窄帶子頻段，然后進(jìn)行處理。這樣在每個子頻段需要根據(jù)中心頻率的不同分別建立方向矢量，大大增加了存儲的空間和計算量。同時該方法不能直接對相干信號進(jìn)行處理，而利用空間平滑的方法解相干又將損失一定的陣列孔徑。

2.2 基于稀疏表示的DOA估計模型

Malioutov[5]提出了一種利用信號稀疏表示進(jìn)行DOA估計的方法。其主要思想是將信號的陣列流型矩陣擴(kuò)展成一個過完備的冗余字典D，它包含了所有可能的源位置信息。令θ= {θ1,θ2, …,θN}代表所有可能的源位置的一個采樣集合，θn代表信號的到達(dá)角度(DOA)，方向矩陣A可擴(kuò)展成如下的過完備字典：

定義N×1的信號向量H={h1,h2,…,hN}，當(dāng)且僅當(dāng)源信號位于角度θn處，hn有非零值，其他分量均為零，可以得到DOA估計的稀疏模型為

然后利用基追蹤(Basic Pursuit, BP)[9]等稀疏分解算法求得H，根據(jù)H中非零元素的位置可以得到信號的 DOA。為了進(jìn)行寬頻段 DOA估計，Malioutov[5]根據(jù)寬帶范圍內(nèi)若干個頻率點的導(dǎo)向矢量矩陣分別建立子字典，然后根據(jù)把這些子字典組合成一個新的字典。

從式(6)可以看出，隨著頻段的展寬，字典的長度極劇增長，稀疏分解的計算量也指數(shù)倍的增長?？紤]同子字典中的兩個元素a(fi,θi)、a(fk,θk)，當(dāng)fisin(θi)=fksin(θk)時，a(fi,θi)=a(fk,θk)。此時，字典的相干度為1，稀疏分解的精度將很難保證[8]。

3 基于空間頻率稀疏表示的寬頻段DOA估計

3.1 基于空間頻率的頻域稀疏模型的構(gòu)建

對于寬頻段的 DOA估計，簡單高效的冗余字典構(gòu)造是一個迫切需要解決的難題。定義空間頻率η=fircosθi/c，則式(3)可表示為

對于均勻線陣為了保證無模糊測向，要求fir/c≤0.5。因此空間頻率的取值范圍為 -0.5≤η≤0.5。根據(jù)空間頻率的取值范圍，設(shè)置N個采樣點，就可以構(gòu)建覆蓋整個無模糊頻段的過完備字典：

為了對寬頻段范圍內(nèi)的多個窄帶信號進(jìn)行DOA估計，我們從頻域構(gòu)建DOA估計的稀疏模型。假設(shè)信號與噪聲相互獨立，中心頻率為f的陣列協(xié)方差為

其中I為單位矩陣，RS(f)是頻率為f信號的協(xié)方差矩陣。常用的稀疏表示 DOA估計方法大多是在時域或頻域建立稀疏模型，容易受到噪聲的影響。本文從子空間域建立DOA估計的稀疏模型，對R做特征值分解得

如果信號協(xié)方差矩陣的秩為K(K≤p)，噪聲協(xié)方差矩陣RN=δ2I為滿秩矩陣，則有如下線性關(guān)系[10]：

其中1≤k≤K,ek為特征矢量，hk(i)為線性組合因子。由式(11)可知，當(dāng)噪聲協(xié)方差矩陣為理想白噪聲時，式(11)即簡化為

式(12)說明無論信號源是否相干，大特征值對應(yīng)的特征矢量是各信號源方向矢量的一個線性組合。令協(xié)方差矩陣的最大特征值為emax，則

根據(jù)式(8)的過完備字典，構(gòu)建式(13)的稀疏模型。

其中H={h1,h2,…,hN}為N×1的信號向量，當(dāng)且僅當(dāng)源信號的空間頻率位于ηn處，hn有非零值，其他分量均為零。對式(14)進(jìn)行稀疏分解得到H,H中非零元素的位置就代表了空間目標(biāo)的空間頻率信息。不同于基于快拍數(shù)據(jù)稀疏分解的 DOA估計方法，特征值分解減弱了噪聲的影響，而選取最大特征值對應(yīng)的特征向量進(jìn)一步減弱了噪聲的影響，因此本文方法在低信噪比情況下魯棒性較好。

3.2 基于頻域峰值協(xié)方差矩陣的頻率和角度估計方法

寬帶接收機(jī)采集的一幀數(shù)據(jù)中，可能包含多個不同頻率點的窄帶信號。首先需要對不同頻點上的窄帶信號進(jìn)行分離。如果利用等分的方法將寬頻段分成若干子頻段再針對每個子頻段進(jìn)行處理，將會對沒有信號的子頻段進(jìn)行運(yùn)算，白白增加了計算復(fù)雜度。因此本文根據(jù)信號在頻域的分布特點，自適應(yīng)地進(jìn)行信號的分離，然后針對有信號的頻點建立模型進(jìn)行DOA估計。具體算法步驟如下：

步驟1 寬頻段內(nèi)窄帶信號頻率估計和分離

首先對陣列接收數(shù)據(jù)進(jìn)行傅里葉變換，即對各通道的接收數(shù)據(jù)分別做N點FFT；其次，測頻并獲取頻域樣本，對任一通道FFT后的結(jié)果進(jìn)行譜峰搜索，并記錄譜峰值和相應(yīng)的頻點坐標(biāo)fk，其它通道根據(jù)頻點坐標(biāo)取出相應(yīng)位置的頻域復(fù)數(shù)值構(gòu)成相應(yīng)的數(shù)據(jù)矢量y(nk)，稱為頻域快拍矢量。

其中S(nk)= [s1(nk),s2(nk),…,sp(nk)]為p個信號在頻域nk的輸出矢量。W(nk)= [w1(nk),w2(nk),…,wM(nk)]為噪聲在頻域nk的輸出矢量。

步驟 2 基于頻域峰值協(xié)方差矩陣特征向量的稀疏DOA估計模型

對式(15)中的頻域峰值快拍矢量求協(xié)方差矩陣：

對式(16)進(jìn)行特征值分解：

利用其最大的特征向量e，根據(jù)3.1節(jié)的變換構(gòu)建DOA估計的稀疏表示模型：

其中D是利用空間頻率構(gòu)建的過完備字典。

步驟3 窄帶多信號頻率和DOA的估計

對式(18)進(jìn)行稀疏分解，得到頻域nk處信號對應(yīng)的空間頻率η= [η1,η2,…,ηp]。其中，p為頻域nk處信號的個數(shù)。則頻率nk處的信號的來波方向可以通過式(19)求得

重復(fù)步驟1到步驟3就可以得到所有寬頻段范圍內(nèi)所有信號的頻率和來波方向。

3.3 稀疏表示的求解

求解式(18)是一個稀疏分解的過程，找到最好的，也就是最稀疏的信號表示，等同于解決下述問題：

當(dāng)上述模型含有噪聲時，式(21)的約束條件不再適用，問題轉(zhuǎn)化為最小化目標(biāo)函數(shù)：

目標(biāo)函數(shù)中前一項反映失配程度，后一項反映稀疏性要求。由于H是復(fù)數(shù)，它的l1范數(shù)為

可以發(fā)現(xiàn)式(23)的兩邊平方仍不能消除平方根項，直接導(dǎo)致我們不能使用二次規(guī)劃的方式最小化目標(biāo)函數(shù)。為了解決這個問題，我們采用二階錐規(guī)劃的方法，其轉(zhuǎn)化形式為

利用內(nèi)點法可以有效解決上述優(yōu)化問題，得到式(14)中稀疏分解的系數(shù)，進(jìn)而根據(jù)非零元素的位置求得信源的DOA估計值。

4 性能分析

本節(jié)主要通過理論和仿真實驗對算法的性能進(jìn)行深入分析，以探討各種的因素對算法性能的影響。實驗中，仿真某型號測向設(shè)備的參數(shù)，天線陣列采用的是5元均勻線陣，陣元間距為0.6 m，測向接收機(jī)的瞬時帶寬為40 MHz，接收機(jī)的中頻為70 MHz，采樣頻率為92.16 MHz。

4.1 有效性分析

假設(shè) 9個窄帶信號入射到天線陣列，9個信號中有頻率重疊的。信噪比為10 dB, FFT變換的點數(shù)為8192點。表1為目標(biāo)的真實值和估計值；圖1為頻率和DOA估計結(jié)果。

圖1 寬頻段的頻率和DOA估計結(jié)果

從表1和圖1可以看出，本文方法對寬頻段內(nèi)的多個信號的頻率和 DOA可以達(dá)到很高的估計精度。頻率的估計誤差小于50 kHz，角度的估計誤差在1°以內(nèi)。具有優(yōu)越的工程應(yīng)用價值。

表1 目標(biāo)的真實值與估計值

表2是本文方法和文獻(xiàn)[5]方法的性能比較。仿真環(huán)境：Matlab7.5 平臺，Pentium 4處理器，1 G內(nèi)存，接收機(jī)的瞬時帶寬為 40 MHz。文獻(xiàn)[5]方法中，每隔10 MHz抽取一個頻點建立過完備字典，字典長度為1795；本文方法利用空間頻率建立字典，字典長度僅為 201，卻能覆蓋整個頻段。本文方法的角度和頻率估計精度遠(yuǎn)優(yōu)于文獻(xiàn)[5]的方法，而運(yùn)算時間卻是后者的三十分之一。

表2 算法的性能比較

4.2 頻率估計誤差的影響

為了降低寬帶稀疏分解時字典的長度，引入了空間頻率的概念。定義為η=frs inθ/c，可以理解為信號頻率和陣元位置引起時延的乘積。因此頻率估計的誤差會對來波方向的估計精度產(chǎn)生影響。信號的波達(dá)方向和頻率存在如下的關(guān)系：

從式(26)可以看出，隨著頻率的升高，測頻誤差對 DOA估計的影響越小。假設(shè)有同頻相干兩個信號入射到天線陣列，信號的角度在0°～90°范圍內(nèi)隨機(jī)產(chǎn)生，分析頻率估計誤差對 DOA估計精度的影響，圖2為均方根誤差隨頻率估計誤差的變化曲線。每個頻率估計誤差下進(jìn)行1000次獨立實驗。

從圖2可以看出，隨頻率估計誤差的增加，角度估計誤差也隨之增加。但是在頻率估計誤差小于1 MHz時，角度估計的均方根誤差小于0.5°。目前基于 FFT變換的測頻技術(shù)在采樣點數(shù)足夠的情況下，可以保證測頻誤差只有幾Hz[11]。因此本文提出的方法可以忽略測頻誤差的影響。

4.3 信噪比的影響

本文提出的方法具有良好的抗噪聲干擾能力，這與算法中采取的兩個運(yùn)算操作有關(guān)：

(1)算法在實現(xiàn)的過程中，利用N點FFT變換將信號從時域變換到頻域。因為對于信號，該變換使信號相干積累，而噪聲在整個頻率均勻分布是非相干疊加，信噪比提高了N倍。

(2)在構(gòu)建信號的稀疏模型中，采用的頻域峰值協(xié)方差矩陣的最大特征向量。特征值分解減弱了噪聲的影響，而選取最大特征值對應(yīng)的特征向量進(jìn)一步減弱了噪聲的影響。

假設(shè)有同頻相干兩個信號入射到天線陣列，信號的角度在0°～90°范圍內(nèi)隨機(jī)產(chǎn)生，分析信噪比對DOA估計精度的影響，圖3為均方根誤差隨信噪比變化關(guān)系。每個信噪比下進(jìn)行1000次獨立實驗。

從圖3可以看出，本文方法對噪聲具有良好的抑制能力。在信噪比-10 dB時，均方根誤差達(dá)到0.6°；信噪比大于0 dB時，均方根誤差小于0.2°。算法的良好抗噪聲性能與理論分析一致。

4.4 同頻信號功率不一致的影響

本文方法采用的是頻域峰值協(xié)方差矩陣的最大特征向量作為稀疏分解的輸入矢量，削弱了噪聲的影響。但是對于同頻率的兩個信號，信號功率的不同會對其造成一定的影響。假設(shè)有兩個同頻但功率不同的信號入射到天線陣列，信號的角度在0°～90°范圍內(nèi)隨機(jī)產(chǎn)生，分析分辨概率同信號功率比之間關(guān)系(如圖4所示)。每個功率比下進(jìn)行1000次獨立實驗。

從圖4可以看出，強(qiáng)信號會對弱信號的估計產(chǎn)生一定的影響，但是這種影響對一般的應(yīng)用場合影響不大。在強(qiáng)弱信號功率比小于200時，分辨概率大于80%，適用于大多數(shù)應(yīng)用場合。對特殊應(yīng)用場合，例如強(qiáng)干擾下的微弱信號 DOA估計，我們采取了歸一化信號子空間的處理方法，將另文討論。

5 結(jié)束語

對寬頻段范圍內(nèi)的多個窄帶信號進(jìn)行頻率和方位角的估計對于無線電測向特別是通信信號測向具有十分重要的意義。本文提出一種基于空間頻率稀疏分解的寬頻段估計方法。理論分析和仿真實驗驗證了本文方法對寬頻段內(nèi)的多個窄帶信號可以達(dá)到很高的估計精度。該方法利用空間頻率代替頻率和方位角的2維組合構(gòu)建過完備字典，字典的長度僅相當(dāng)于窄帶信號 DOA估計的字典長度，卻可以覆蓋整個測向頻段；利用FFT變換和頻域峰值協(xié)方差矩陣的特征向量大大降低了噪聲的影響，估計的信號數(shù)目能遠(yuǎn)大于陣元數(shù)目。論文主要討論對寬頻度內(nèi)的多個窄帶信號進(jìn)行來波方位估計，下一步將對多個相干或不相干寬帶信號進(jìn)行估計。

圖2 均方根誤差隨頻率估計誤差變化關(guān)系

圖3 均方根誤差隨信噪比變化關(guān)系

圖4 分辨概率同信號功率比之間關(guān)系

[1]Schmidt R O. Multiple emitter location and signal parameter estimation [J].IEEE Transactions on Antennas Propagation,1986, 34(3): 276-280.

[2]劉聰鋒, 廖桂生. 寬帶接收機(jī)的窄帶信號頻率和2維角度估計新方法[J]. 電子學(xué)報, 2009, 37(3): 523-528.Liu Cong-feng and Liao Gui-sheng. Novel method of narrow band signal frequency and 2D angle estimation for wide band receiver [J].Acta Electronica Sinica, 2009, 37(3): 523-528.

[3]Guo Xian-sheng, Wan Qun, Chang Chun-qi,et al.. Source localization using a sparse representation framework to achieve superresolution [J].Mulitidimensional Systems and Signal Processing, 2010, 21(4): 391-402.

[4]Md Mashud Hyder and Kaushik Mahata. Direction-of-Arrival estimation using a mixedL2,0norm approximation [J].IEEE Transactions on Signal Processing,2010, 58(9):4646-4655.

[5]Malioutov D, Cetin M, and Willsky A S. A sparse signal reconstruction perspective for source localization with sensor arrays [J].IEEE Transactions on Signal Processing, 2005,53(8): 3010-3022.

[6]Zhang Shu, Li Yu-mei, and Song Jun-cai. A novel method for DOA estimation based on generalized-prior distribution [C].2010 International Conference on Measuring Technology and Mechatronics Automation, Changsha, China, 2010: 244-247.

[7]Duarte M F, Wakin M B, and Baraniuk R G. Fast reconstruction of piecewise smooth signals from incoherent projections. Workshop on Signal Processing with Adaptative Sparse. Structured Representations-SPARS’05, Rennes ，F(xiàn)rance, 2005. http://spars05.isisa.fr/ ACTES/TS5-3.pdf.

[8]Donoho D and Huo X. Uncertainty principles and ideal atomic decompositions [J].IEEE Transactions on Information Theory, 2001, 47(7): 2845-2862.

[9]Chen Hao-bing, Donoho D L, and Saunders M A. Atomic decomposition by basis pursuit[J].SIAM Joumal on Scientific Compting, 2001, 20(1): 33-66.

[10]Cadzow J A, Kim Y S, and Shiue D C. General direction-ofarrival estimation: a signal subspace approach.IEEE Transactions on Aerospace and Electronic Systems, 1989,25(1): 31-47.

[11]萬靈達(dá), 楊曉光. 一種基于 FFT的高精度頻率估計算法[J].電子科技, 2010, 23(10): 79-81.Wan Ling-da, Yang Xiao-guang. A high-accuracy frequency estimation algorithm based on FFT [J].Electronic Science&Technique, 2010, 23(10): 79-81.