基于空頻域稀疏表示的寬頻段DOA估計(jì)

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

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

(安徽省電子制約技術(shù)重點(diǎn)實(shí)驗(yàn)室 合肥 230037)

1 引言

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

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

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

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

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

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

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

c為電磁波傳播的速度。

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

2.2 基于稀疏表示的DOA估計(jì)模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟3 窄帶多信號(hào)頻率和DOA的估計(jì)

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

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

3.3 稀疏表示的求解

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

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

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

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

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

4 性能分析

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

4.1 有效性分析

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

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

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

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

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

表2 算法的性能比較

4.2 頻率估計(jì)誤差的影響

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

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

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

4.3 信噪比的影響

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

(1)算法在實(shí)現(xiàn)的過程中，利用N點(diǎn)FFT變換將信號(hào)從時(shí)域變換到頻域。因?yàn)閷?duì)于信號(hào)，該變換使信號(hào)相干積累，而噪聲在整個(gè)頻率均勻分布是非相干疊加，信噪比提高了N倍。

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

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

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

4.4 同頻信號(hào)功率不一致的影響

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

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

5 結(jié)束語

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

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

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

圖4 分辨概率同信號(hào)功率比之間關(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ī)的窄帶信號(hào)頻率和2維角度估計(jì)新方法[J]. 電子學(xué)報(bào), 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ì)算法[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.