基于空頻域稀疏表示的寬頻段DOA估計

李鵬飛 張 旻 鐘子發 羅 爭

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

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

1 引言

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

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

2 數學模型

2.1 陣列信號的數學模型

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

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

其中第i個導向矢量為

c為電磁波傳播的速度。

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

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

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

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

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

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

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

3.1 基于空間頻率的頻域稀疏模型的構建

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

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

為了對寬頻段范圍內的多個窄帶信號進行DOA估計，我們從頻域構建DOA估計的稀疏模型。假設信號與噪聲相互獨立，中心頻率為f的陣列協方差為

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

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

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

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

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

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

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

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

步驟1 寬頻段內窄帶信號頻率估計和分離

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

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

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

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

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

利用其最大的特征向量e，根據3.1節的變換構建DOA估計的稀疏表示模型：

其中D是利用空間頻率構建的過完備字典。

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

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

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

3.3 稀疏表示的求解

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

當上述模型含有噪聲時，式(21)的約束條件不再適用，問題轉化為最小化目標函數：

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

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

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

4 性能分析

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

4.1 有效性分析

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

圖1 寬頻段的頻率和DOA估計結果

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

表1 目標的真實值與估計值

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

表2 算法的性能比較

4.2 頻率估計誤差的影響

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

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

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

4.3 信噪比的影響

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

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

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

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

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

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

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

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

5 結束語

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

圖2 均方根誤差隨頻率估計誤差變化關系

圖3 均方根誤差隨信噪比變化關系

圖4 分辨概率同信號功率比之間關系

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

[2]劉聰鋒, 廖桂生. 寬帶接收機的窄帶信號頻率和2維角度估計新方法[J]. 電子學報, 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 ，France, 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]萬靈達, 楊曉光. 一種基于 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.