色噪聲下基于差分聚焦的寬帶DOA估計方法

曹司磊，曾維貴，王 磊

(海軍航空大學，山東 煙臺 264001)

波達方向(direction-of-arrival,DOA)估計是陣列信號處理中的重要內容，在雷達、通信及醫學領域有著廣泛的應用.空間譜方法是DOA估計的經典方法，如多信號分類(multiple signal classification,MUSIC)方法、旋轉子空間不變(estimating signal parameter via rotational invariance techniques,ESPRIT)算法等.但是，傳統的空間譜方法一般假設空間噪聲為高斯均勻白噪聲.在實際應用中，大多數陣列信號處理中面臨色噪聲環境.因此，有學者提出了色噪聲環境下的陣列DOA估計改進方法，主要可分為3類：參數化方法[1]、四階累積量算法[2-3]及協方差矩陣差分法[4-5].參數化方法一般需要預先知道噪聲協方差矩陣結構特征，以此來進行預白化處理；四階累積量算法通過構造觀測信號的四階累積量來消除非高斯色噪聲影響；協方差差分法一般假設色噪聲協方差矩陣滿足Toeplitz結構，并通過將觀測信號協方差矩陣與線性變換后的協方差矩陣進行差分處理來消除色噪聲影響.

上述針對色噪聲環境下DOA估計算法的大多數只適用于窄帶陣列，在寬帶陣列中的使用受限.目前，寬帶DOA估計主要解決途徑有：基于頻率分解的非相干信號子空間方法[6-7](incoherent signal subspace method, ISSM)及基于頻率聚焦的相干信號子空間方法[8-10](coherent signal subspace method,CSSM).ISSM方法對信噪比要求較高，低信噪比條件下性能急劇下降；CSSM方法在低信噪比條件下性能優于ISSM方法，但存在的問題是需要預估或者假設信號角度，以此來構建信號子空間，這種角度估計會導致聚焦矩陣出現偏差，最終導致DOA估計誤差.為避免CSSM方法中角度預估帶來的失配誤差，有學者提出了改進算法.文獻[11]提出基于頻域時延補償的DOA估計方法(DCF)，避免了角度預估的問題，但是算法復雜度較高，用于實時陣列信號處理難度較大.文獻[12]提出使用陣列自相關矩陣構造聚焦矩陣，避免了角度預估帶來的誤差，但是算法的估計性能有所下降.

結合來看，對于色噪聲下的寬帶陣列DOA估計問題，研究者提出解決一些方法.文獻[13]提出一種低復雜度寬帶色噪聲DOA估計方法，利用空間差分技術去除Toeplitz結構的色噪聲，然后通過RSS方法計算通用協方差矩陣，最后應用低復雜度的PM算法估計信號DOA；但是該算法對信源數較為敏感.文獻[14]提出了一種基于數據協方差矩陣修正及TOPS算法相結合的寬帶DOA估計方法(data covariance matrix correction with TOPS algorithm for DOA estimation, DCT-DE)，利用特殊構造矩陣將色噪聲協方差矩陣從數據協方差矩陣中剔除，然后利用平方TOPS算法實現寬帶DOA估計.文獻[15]提出基于FOCUSS的網格角度求解算法(MMV-FOCUSS)，通過差分算法去除色噪聲影響并構建新觀測模型，采用分布優化的方式進行DOA估計，但是算法復雜度較高.

針對色噪聲背景下寬帶DOA估計存在的上述問題，本文提出基于差分聚焦的寬帶陣列DOA估計方法.算法首先對信號協方差矩陣進行差分處理，消除色噪聲影響；然后對差分矩陣進行特征分解，以差分矩陣正特征值對應的信號子空間為基礎重構觀測信號模型；在此基礎上，分析了無需角度預估的特征向量空間聚焦DOA估計原理，并對特征向量矩陣進行重新分塊排列，簡化重構聚焦矩陣.理論分析及仿真結果表明，本文算法估計精度及穩健性較好，且不需角度預估、復雜度低.

1 陣列信號模型

陣列排布為一維線陣，陣元數為K，陣元間距為d.空間中有I個寬帶遠場信號s(t)=[s1(t),…,sI(t)]入射至陣列，且信號與噪聲相互獨立，入射角為θ=[θ1,…,θI]，則第k個陣元輸出信號xk(t)可以表征為

式中,τk,i為信號si(t)在第k個陣元相對于首個陣元的接收時延；nk(t)為第k個陣元接收非均勻噪聲矢量，并假設空間噪聲場為各向同性圓形或柱形場.

將足夠長時間內陣列接收到的信號分割成若干獨立快拍，對每個快拍內的數據進行離散傅里葉變換(discrete Fourier transform, DFT)，得到J個頻域窄帶信號.寬帶陣列接收信號的頻域表示形式如下:

X(fj)=A(fj,θ)S(fj)+N(fj),j=1,…,J;

X(fj)=[Xk(fj),…,XK(fj)];

A(fj,θ)=[a(fj,θ1),…,a(fj,θI)];

S(fj)=[Si(fj),…,SI(fj)];

N(fj)=[Nk(fj),…,NK(fj)].

式中,A(fj,θ)為頻點fj處的陣列流形；a(fj,θi)為入射角為θi的信號在頻點fj的方位導向矢量.

頻點fj處的陣列輸出信號自相關矩陣可表示為

R(fj,θ)=E[X(fj,θ)XH(fj,θ)]=
A(fj,θ)Rs(fj,θ)AH(fj,θ)+
E[N(fj)NH(fj)]=
A(fj,θ)Rs(fj,θ)AH(fj,θ)+QN.

式中QN表示噪聲信號協方差矩陣.高斯相關色噪聲協方差矩陣QN的第(n1,n2)個元素為[16]

QN(n1,n2)=σ2γ|n1-n2|e|jpi(n1-n2)/2|.

式中σ2為噪聲功率，用于調整信噪比參數；γ為回歸系數，用于調整噪聲相干性.

2 協方差矩陣差分理論

由于假設空間噪聲場為同性圓形或柱形場，因此色噪聲信號協方差矩陣滿足Toeplitz結構[4].可知，噪聲協方差矩陣分量旋轉不變，有

QN=ΦQNΦ.

式中Φ表示轉置矩陣.

對協方差矩陣變化前后求差得到

(1)

從上式可以看出，在差分協方差矩陣ΔR中，噪聲分量已被消除.對[A,ΦA]進行分解，有

(2)

其中(·)*表示矩陣共軛.將式(2)代入式(1)中，可得

從上式可以看出，A對應信源入射真實角度，而A*則對應信源入射負角度.因此，利用ΔR進行DOA估計的過程中，不可避免地會出現偽峰.

3 DOA估計方法

3.1 差分協方差矩陣特征分解

針對協方差矩陣差分過程導致的偽峰問題，本文對差分矩陣ΔR進行特征分解，取正特征值對應的信號子空間作為DOA估計基礎，有

ΔR=UΛUH=[U+,Un,U-]Λ[U+,Un,U-].

式中：U+是由I個正特征值對應的特征向量構成的特征矩陣；Un對應噪聲子空間；U-是由N個負特征值對應的特征向量構成的特征矩陣；Λ可分別表示為

從上述分析知，U+可張成導向矢量A，因此可建立新的觀測模型如下：

U+=Aq+ε.

式中：q為新觀測模型對應的信號矩陣；ε為新觀測模型對應的噪聲分量.新的觀測模型經過差分分解處理后，不再含有色噪聲分量，同時包含了全部目標方位信息，可較好用于寬帶DOA估計.

3.2 DOA估計方法

3.2.1 DOA估計基本原理

在構建新的觀測模型基礎上，本文提出特征向量信號子空間聚焦的方法進行寬帶DOA估計.

首先，對于新的觀測模型，信號自相關矩陣R+可表示為

(3)

其中P表示無噪聲信號自相關矩陣.P為Hermitian矩陣，因此可酉相似于對角矩陣，表示為

P=VMVH.

本文以子空間Γj為聚焦空間，設定參考聚焦矩陣Ψ0，則有

(4)

式中Ω為參考聚焦矩陣Ψ0在聚焦空間Γj中的系數矩陣.聚焦矩陣為T，則T滿足

按照文獻[10]推導，有

(5)

式中ξ為參考聚焦矩陣Ψ0的特征向量矩陣.若要實現完美聚焦，則應有

上式可化為

(6)

式中,V0及Vj分別為頻點f0及fj處的特征向量矩陣.觀察式(6)，為了方便計算合并，令

式中W為未知矩陣.則式(6)可進一步化為

式中B=Ω-WMjWH.由矩陣的范數特性可將上式進一步化為

(7)

可知，W為酉矩陣，Mj為對角矩陣，則WMjWH亦為對角矩陣.因此，若要式(7)為0，可行解之一為

式中Ω為對角矩陣.對照式(4)，可知V0即為參考聚焦矩陣Ψ0的特征向量矩陣.因此，根據式(5)，有

3.2.2 聚焦矩陣重構

為避免特征分解過程中零特征值所對應的特征向量對分辨門限的影響，對聚焦矩陣進行重構.目標信號的方位區間為β，對區間范圍內按角度間隔進行功率譜積分可得

(8)

對Π(fj)進行特征值分解，得出其顯著特征值個數ρ(即為信號數).

將V0及Vj按照對應特征值遞減的順序進行重新排列，表示為

式中V′0、V′j為K*ρ矩陣，即是由信號協方差矩陣特征分解后前ρ個顯著特征值對應的特征向量構成.Vj為酉矩陣，則可知

(V″j)HA(fj)=AH(fj)V″j=0.

因此，

V′0(V′j)HPjV′j(V′0)H.

可得重構后的聚焦矩陣

(9)

從上式看出，重構后的聚焦矩陣避免了協方差分解過程中非顯著特征值所對應的特征向量，提高了算法在整個信噪比范圍內的適應能力及分辨性能.

(10)

3.3 算法實現步驟

綜上，本文算法步驟總結為：

1)將陣列接收信號分割為獨立快拍片段，進行DFT處理，計算各頻點處的陣列輸出信號自相關矩陣R(fj,θ)；

2)根據式(1)計算差分協方差矩陣ΔR(fj,θ)，并對其進行特征分解得U+；

3)按照式(3) 求取R+(fj)，對R+(fj)進行特征分解，得V0與Vj；

4)按照式(8)求取Π(fj)，對Π(fj)進行特征分解，得ρ；

4 仿真分析

本小節通過仿真比較實驗，分析本文所提算法的測向準確性、角分辨能力、算法穩健性及復雜度特性.基本仿真參數為：陣元數為8，工作帶寬覆蓋X波段(8～12 GHz)，信號帶寬500 MHz，采樣頻率為2.4 GHz，陣元間距為2.5 cm.仿真實驗將本文算法與DCT-DE算法、MMV-FOCUSS算法及文獻[13]算法進行比較分析.其中，MMV-FOCUSS算法的網格間距取1°.

4.1 噪聲背景下算法的測向精確性比較

在高斯白噪聲與高斯色噪聲背景下，仿真比較本文算法與另外3種算法在不同噪聲背景下的測向準確性.仿真實驗中，設定2個獨立的寬帶輻射源，入射角度分別為-10°及20°，算法采用快拍數為200，白噪聲信噪比設定為5 dB，色噪聲信噪比設定為5 dB.4種算法的空間譜見圖1、2.

圖1 高斯白噪聲下空間譜

圖2 高斯色噪聲下空間譜

從圖1看出，白噪聲背景下本文算法、MMV-FOCUSS算法、DCT-DE算法及文獻[13]算法均具備較好的DOA估計性能.從圖2看出，高斯色噪聲背景下，文獻[13]算法及DCT-DE算法譜峰出現較大偏差，MMV-FOCUSS算法測向精度稍好，但估計方差仍在2°以上，而本文算法保持較好的估計性能，測向精度較高.

4.2 噪聲背景下算法的分辨能力比較

本次仿真實驗主要比較本文算法與其他3種算法在高斯色噪聲背景下的角度分辨能力.仿真條件為：空間環境中設定兩個獨立的等功率輻射源，快拍數為500，噪聲參數同上.入射角度間隔Δθ分別設定為Δθ=[2°,4°,8°,15°]，即入射角度θ為

θ=[(θ1,θ2)]=

[(-1°,1°),(-2°,2°),(-4°,4°),(-7°,8°)].

仿真實驗中，通過200次蒙特卡洛仿真實驗，統計各類算法的平均分辨能力，比較分析4種算法在色噪聲背景下的角度分辨性能.本文對角度分辨率η的定義為

Qmin=min[Q(θ)],θ∈[θ1,θ2].

可知，η代表兩譜峰與峰間谷值之間的平均差，η越高，表示分辨能力越強.統計4種算法的平均分辨率見表1.

表1 算法分辨能力比較

從表1看出，橫向比較表格數據，兩輻射源信號入射角度的間隔越大，則各算法的角度分辨率越高；縱向比較發現，相同的空間色噪聲環境下，本文算法的角度分辨能力相較于MMV-FOCUSS及DCT-DE等算法更好.

4.3 算法穩健性分析

本節對4種算法的穩健性進行仿真分析，主要包括估計方差及快拍數影響等.仿真實驗中，單次估計誤差<3°認為是有效實驗，而估計方差及快拍數影響是對所有有效實驗結果的統計.快拍數影響分析實驗中，信噪比取5 dB.背景噪聲為高斯色噪聲，參數不變.仿真結果為200次蒙特卡洛仿真實驗的統計平均.

從圖3看出，色噪聲背景下，隨著輸入信噪比的不斷提高，本文算法與MMV-FOCUSS始終保持較好的估計性能；而文獻[13]算法在低信噪比條件下的估計性能較差.從圖4看出，隨著快拍數的提升，文獻[13]算法的性能提升較大，但總體性能差于本文算法及MMV-FOCUSS算法；在仿真所取的快拍數范圍內，本文算法性能始終優于其他算法.

圖3 均方根誤差隨信噪比變化曲線

圖4 均方根誤差隨快拍數變化曲線

4.4 算法復雜度分析

仿真條件設置為：設定兩輻射源入射場景及三輻射源入射場景，入射角度分別為θ=[0°,10°]及θ=[-10°,0°,10°]，快拍數取200，信噪比取5 dB.仿真平臺為MATLAB 2014a，運行于聯想P1移動工作站，進行500次蒙特卡洛實驗，求取算法平均運行時間，結果見表2.

表2 算法運行時間比較

分析表2，本文算法的運行時間最短，MMV-FOCUSS算法的運行時間最長，而文獻[13]算法與DCT_DE算法運行時間相當.從算法實現過程看，MMV-FOCUSS算法需要在每個角度間隔內進行迭代運算，造成的計算量較大；文獻[13]算法在估計單頻點參考聚焦矩陣的過程中，需要進行2次矩陣逆運算及4次矩陣乘法；DCT-DE算法需要進行3次特征值分解及2次矩陣乘法；本文算法需要進行3次特征值分解及1次矩陣乘法.因此，本文算法的運算復雜度相對較小.

5 結 論

針對色噪聲背景下寬帶陣列DOA估計精度差及穩健性不足的問題，本文提出基于差分聚焦的寬帶陣列DOA估計方法.利用信號自相關矩陣差分去除色噪聲影響并以其正特征值對應的信號子空間為基礎構建新的觀測模型；通過對新模型信號協方差矩陣特征向量子空間的分析，提出無需角度預估的聚焦矩陣計算方法.理論分析及仿真實驗證明，色噪聲背景下，本文方法相較于對比算法具有更好的估計精度及分辨率，且算法穩健性好、復雜度低.