基于投影矩陣搜索的相干源DOA 估計算法*

宗佳虹，戴旭初

（中國科學(xué)技術(shù)大學(xué) 電子工程與信息科學(xué)系 合肥 230027）

引 言

波達(dá)方向（DOA）估計是陣列信號處理的研究熱點，在雷達(dá)[1]、導(dǎo)航[2]等領(lǐng)域具有廣泛的應(yīng)用前景。傳統(tǒng)的DOA 估計算法中，多重信號分類算法MUSIC[3]、旋轉(zhuǎn)不變子空間算法ESPRIT[4]等子空間類算法具有良好的估計性能和超分辨能力，但是對于相干源，陣列輸出信號的協(xié)方差矩陣的秩虧缺導(dǎo)致子空間估計錯誤，從而使算法失效[5,6]。

針對相干源DOA 估計問題，許多學(xué)者提出了解相干算法，例如空間平滑法[7,8]、Toeplitz 矩陣重構(gòu)法[9]、基于SVD 的解相干算法[10]等。文獻(xiàn)[7]提出的空間平滑算法具有良好的解相干性能，且計算量小，被廣泛應(yīng)用于實際工程中。但是，空間平滑算法存在陣列孔徑的損失，從而降低估計的分辨率與精度，特別是，空間平滑算法只適合于特殊陣列，即需要子陣列結(jié)構(gòu)具有移不變性，例如均勻線陣。

對于一般的陣列，可以考慮采用虛擬內(nèi)插變換法[11,12]將其變換為均勻陣列，再進(jìn)行空間平滑等解相干處理。虛擬內(nèi)插變換是在某扇區(qū)范圍內(nèi)進(jìn)行內(nèi)插，將原陣列變換為一個虛擬陣列。不過，虛擬陣列變換只能針對某一入射角范圍，且會引入額外的誤差，內(nèi)插范圍越大，變換誤差越大，從而導(dǎo)致DOA 估計性能降低。

針對上述問題，本文提出了一種基于投影矩陣搜索的DOA 估計算法，其基本思想是由陣列流型導(dǎo)出噪聲子空間的投影矩陣，將接收信號樣本投影到噪聲子空間，獲得空間譜，進(jìn)而得到相干源的DOA估計。該算法適用于任意陣列結(jié)構(gòu)（如非均勻線陣），且無需進(jìn)行虛擬內(nèi)插、空間平滑等預(yù)處理，避免了變換誤差和孔徑損失，具有良好的估計性能。

1 信號模型

考慮K個遠(yuǎn)場窄帶平穩(wěn)信號入射到空間某窄帶陣列上，設(shè)陣列的陣元數(shù)為M(M＞K)，信號中心頻率為ω，第m個陣元在時刻t的輸出信號可表示為：

其中，si(t)為第i個窄帶源信號，θi為第i個信源到達(dá)陣列的入射角，τm(θi)為第i個源信號到達(dá)第m個陣元時相對于參考陣元的時延，v m(t)為第m個陣元的加性噪聲，其功率為σ2。將M個陣元t時刻的接收信號表示成矢量形式，則有：

設(shè)快拍數(shù)為N，則接收到的N個信號樣本用矩陣形式可表示為：

其中，X=[x(t1),…,x(tN)]，S=[s(t1),…,s(tN)]，V=[v(t1),…,v(tN)]。

定義陣列接收信號的協(xié)方差矩陣為：

其中，λi和ui分別為R的第i個特征向量和特征值，且M個特征值滿足如下關(guān)系：

基于特征空間理論，Us=[u1,u2,…,uK]構(gòu)成了接收信號的信號子空間，Un=[uK+1,uK+2,…,uM]構(gòu)成了接收信號的噪聲子空間，且Us與Un正交。子空間類算法就是利用該正交性進(jìn)行DOA 估計，因此傳統(tǒng)的子空間類算法需要對子空間進(jìn)行準(zhǔn)確的估計。

實際應(yīng)用環(huán)境中，源信號之間可能是相關(guān)或相干的。設(shè)兩平穩(wěn)信號si(t)、s k(t)，其相關(guān)系數(shù)定義為：

由Schwartz 不等式可知|ρik| ≤ 1，當(dāng)|ρik|=0時，稱si(t)與s k(t)為獨立信號源；當(dāng)0＜|ρik|＜1時，稱si(t)與s k(t)為相關(guān)信號源；當(dāng)|ρik|=1時，稱si(t)與s k(t)為相干信號源[13]。

對于獨立或弱相關(guān)信號源，信號子空間的維數(shù)等于信源數(shù)K；而對于強相關(guān)或相干信號源，信號子空間的維數(shù)小于K，也就是說，相干信號源會導(dǎo)致協(xié)方差矩陣R秩虧缺，從而使子空間Us與Un估計錯誤。因此，傳統(tǒng)的子空間類算法難以有效估計相干源DOA。

2 基于虛擬陣列變換和空間平滑的MUSIC 算法

本節(jié)簡要介紹已有的相干源DOA 估計方法。對于具有移不變性的陣列結(jié)構(gòu)，常用的相干源DOA估計方法是J.E.Evans 等人提出的空間平滑算法[7]；對于一般的非移不變的陣列結(jié)構(gòu)，考慮采用虛擬內(nèi)插變換法將其變換為具有移不變性的虛擬陣列[11]，再進(jìn)行空間平滑。

空間平滑算法是將陣列分成若干個相同結(jié)構(gòu)的子陣，對每個子陣的協(xié)方差矩陣進(jìn)行平均運算，這樣得到的平滑后的協(xié)方差矩陣，其秩得以恢復(fù)，從而實現(xiàn)解相干。但是，此方法是以犧牲陣列孔徑為代價的，且只適用于具有移不變性的陣列結(jié)構(gòu)。

虛擬內(nèi)插變換法最早由Friedlander.B 等人提出。假設(shè)源信號位于某個角度范圍Θ內(nèi)，將Θ區(qū)間離散化，即將Θ劃分為r個離散角度的集合，Θ=[θ1,θ2…,θr]，其中θ1、θr為角度范圍的左、右邊界，且r＞M，則與Θ=[θ1,θ2…,θr]相對應(yīng)的原陣列的陣列流型為：

通過對虛擬陣列進(jìn)行傳統(tǒng)的空間平滑處理，可實現(xiàn)解相干，然后可采用一般的空間譜估計方法進(jìn)行測角，如MUSIC 算法。

值得注意的是，虛擬內(nèi)插變換需要到達(dá)角的范圍，即需要確定信源到達(dá)角的大致位置。這里考慮采用Capon 算法[14]確定虛擬內(nèi)插變換的角度范圍。Capon 算法的基本思想是在期望方向信號功率不變的情況下，最小化總功率，也就是最小化噪聲及非信源方向干擾功率，該算法的空間譜表示為：

由于Capon 算法對相干性不敏感，但與子空間類算法相比，其性能較差，且分辨率受波束寬度限制，因此考慮采用Capon 算法對空間譜進(jìn)行一次搜索，取每個譜峰值處的一個波束寬度范圍，取它們的并集作為虛擬內(nèi)插的范圍。

綜上，基于虛擬陣列變換和空間平滑的MUSIC 算法流程總結(jié)如下。

輸入：信源數(shù)K，陣列流型A，接收信號樣本X；

輸出：各信源的DOA 估計。

①用Capon 算法進(jìn)行一次搜索，確定虛擬內(nèi)插范圍；

② 根據(jù)式（11）、（12）將原陣列變換為虛擬均勻陣列；

③針對虛擬均勻陣列，用空間平滑法進(jìn)行解相干處理；

④ 利用MUSIC 算法得到DOA 估計。

理論上，基于虛擬陣列變換的空間平滑MUSIC 算法可以進(jìn)行相干源DOA 估計，但是，虛擬陣列變換中引入的變換誤差會降低最后的DOA 估計性能，且內(nèi)插范圍越大，變換誤差越大。另外，空間平滑算法導(dǎo)致的孔徑損失也會降低DOA 估計的精度與分辨率。

3 本文算法

現(xiàn)有的基于空間譜DOA 估計方法，都是利用陣列接收到的數(shù)據(jù)來構(gòu)造信號子空間和噪聲子空間，然后將陣列的導(dǎo)向矢量投影到噪聲子空間來獲得空間譜，并通過對導(dǎo)向矢量的搜索得到DOA 估計。與現(xiàn)有的方法不同，本節(jié)將提出一種DOA 估計的新方法，即基于投影矩陣搜索的DOA 估計算法，其基本思想是利用陣列結(jié)構(gòu)來構(gòu)造信號子空間的投影矩陣，再根據(jù)信號子空間和噪聲子空間的正交性，得到噪聲子空間的投影矩陣，然后將陣列接收到的數(shù)據(jù)投影到噪聲子空間來構(gòu)造空間譜，從而通過投影矩陣的搜索得到空間譜和DOA 估計。

3.1 算法設(shè)計

3.1.1 噪聲子空間的構(gòu)造

設(shè)信源數(shù)為K，將K個信源的DOA 寫成矢量形式θ=[θ1,…,θK]，由于信號子空間與導(dǎo)向矢量張成的空間是同一個空間，因此，利用陣列流型A(θ)=[a(θ1),…,a(θK)]，可以得到信號子空間的投影矩陣[15]：

由式（14）、（15）可以看出，本文算法的子空間是根據(jù)陣列流型構(gòu)造的，與信號無關(guān)。

3.1.2 數(shù)據(jù)樣本的降維處理

為了提高該算法的抗噪聲能力，同時降低算法的復(fù)雜度，這里先對信號樣本進(jìn)行時域降維預(yù)處理。將信號樣本X進(jìn)行奇異值分解，得到：

其中，U為左奇異矩陣，V為右奇異矩陣，D為一M×N維的對角陣；取V的前K列，記為V′，取D的前K列和前K行，記為D′，即保留信號空間的功率、去除噪聲空間的功率，則降維處理后的數(shù)據(jù)樣本為：

通常，N遠(yuǎn)大于K，通過上述處理將M×N維信號樣本X降維成M×K維的X′，減小了后續(xù)的投影運算的計算量。

3.1.3 數(shù)據(jù)樣本的投影和空間譜

將處理過的信號樣本X′投影到噪聲子空間，得到

由于X′始終位于信號空間，當(dāng)θ為真實的到達(dá)角時，PN(θ)對應(yīng)的噪聲子空間與真實的信號空間正交，此時PN(θ)與X′正交，Z=O；當(dāng)θ不是真實到達(dá)角時，該正交性不成立，即PN(θ)與X′不正交，Z≠O。

3.1.4 基于空間譜搜索的DOA 估計

因此，在可能的到達(dá)角范圍內(nèi)，通過對投影矩陣PN(θ)的搜索獲得空間譜，參考上節(jié)的算法一，可能的到達(dá)角范圍也可利用Capon 算法得到，本文算法的空間譜表示為：

最后通過峰值搜索，得到DOA 估計。

3.1.5 算法流程

基于上述分析和討論，本文提出算法的流程總結(jié)如下。

輸入：信源數(shù)K，陣列流型A，接收信號樣本X；

輸出：各信源的DOA 估計。

①用Capon 算法進(jìn)行一次搜索，確定可能的到達(dá)角范圍；

② 在可能的到達(dá)角范圍內(nèi)，根據(jù)式（14）和（15）計算噪聲子空間的投影矩陣；

③根據(jù)式（16）、（17）將接收信號進(jìn)行降維處理；

④ 根據(jù)式（18）、（19）計算空間譜，尋找峰值點，得到DOA 估計。

與算法一相比，本文算法是由陣列流型導(dǎo)出噪聲子空間，而不是由數(shù)據(jù)的協(xié)方差矩陣導(dǎo)出的噪聲子空間，因此理論上本文算法與信號的相干性無關(guān)，故對于相干源DOA 估計問題，本文算法無需進(jìn)行解相干處理，從而避免陣列孔徑的損失。此外，由于本文算法在滿足布陣要求的情況下對不同陣列均可適用，無需進(jìn)行虛擬變換處理，從而避免引入內(nèi)插誤差。

3.2 復(fù)雜度分析

設(shè)信源數(shù)為K，陣元數(shù)為M，快拍數(shù)為N，首先考慮基于虛擬陣列變換和空間平滑的MUSIC 算法（算法一）的復(fù)雜度，估計協(xié)方差矩陣的計算量為M2N，設(shè)虛擬內(nèi)插變換的范圍為Θ=[θ1,θ2…,θr]，虛擬陣列變換的計算量為3M3+2M2r；設(shè)空間平滑子陣列的陣元數(shù)為m，則空間平滑的計算量為2(m2M+M2m)(M-m+1)；最后采用MUSIC 算法在Θ上進(jìn)行一維搜索，計算量為2m2K+K3。因此算法一的總的計算量約為M2N+3M3+2M2r+2(m2M+M2m)(M-m+1) +(2m2K+K3)r，由于r?M,N?M，計算量近似為O(M2N)+O(2m2K+K3)r。

下面考慮本文算法（算法二）的復(fù)雜度，信號樣本矩陣奇異值分解的復(fù)雜度為O(MN2)，投影矩陣的計算量為2MK2+M2K+K3，投影計算的計算量為M2K，本文算法需進(jìn)行K維搜索，因此，算法二的計算量約為O(MN2)+O(2MK2+2M2K+K3)rK。

通過上述分析可知，在空間譜搜索部分，算法一是進(jìn)行一維搜索，而算法二需進(jìn)行K維搜索，因此算法二的計算復(fù)雜度較高。

4 仿真實驗及分析

本節(jié)通過仿真實驗對現(xiàn)有算法（算法一）和本文算法（算法二）的性能進(jìn)行分析和比較。

仿真實驗的條件設(shè)置為：五陣元的非均勻線陣(M=5)，陣元設(shè)置為0.25λ[0,1,3,6,8]，算法一的虛擬陣列設(shè)置為0.5λ[0,1,2,3,4]，λ為波長；源信號為窄帶高斯隨機信號，相干源的相關(guān)系數(shù)為1；實驗中，信噪比SNR 定義為陣列接收信號的總功率與噪聲總功率之比，即：

另外，DOA 估計精度由均方根誤差RMSE 來衡量，不失一般性，設(shè)獨立實驗次數(shù)為n，信源數(shù)為K，真實DOA 為 (θ1,…,θK)，第i次實驗的DOA 估計為，RMSE 定義為：

4.1 實驗一：相干源的估計精度比較

參數(shù)設(shè)定：信源數(shù)K=2，DOA 為(70°,80 °)，兩信源相干；快拍數(shù)N=200；獨立實驗次數(shù)n=500；算法一的空間平滑子陣列的陣元個數(shù)取3。圖1 為兩種算法的RMSE 與SNR 的關(guān)系：圖中藍(lán)色曲線表示算法一的RMSE 隨SNR 的變化關(guān)系；紅色曲線表示算法二的RMSE 隨SNR 的變化關(guān)系。

從圖1 可以看出，對于相干源，算法二較算法一有更高的估計精度，在SNR 取值范圍為2～12 dB左右時，算法一基本失效，而算法二仍可以有效估計DOA；在SNR 較高時，兩種算法均可有效估計相干源DOA，但算法二的RMSE 較低，相差約0.5°～0.9°。

圖1 兩種算法RMSE 隨SNR 變化曲線Fig.1 RMSE of two kinds of algorithm with various SNR

4.2 實驗二：相干源的分辨率性能比較

參數(shù)設(shè)定：信源數(shù)K=2，DOA 為(70°,70°)ψ+，兩信源相干；快拍數(shù)N=200；SNR=5 dB；獨立實驗次數(shù)n=500；算法一的空間平滑子陣列的陣元個數(shù)取3。圖2 為兩種算法的RMSE 與角度差ψ的關(guān)系：圖中藍(lán)色曲線表示算法一的RMSE 隨角度差ψ的變化關(guān)系；紅色曲線表示算法二的RMSE 與角度差ψ的變化關(guān)系。

從圖2 可以看出，對于相干源，在不同的信噪比下，算法二較算法一都有更高的分辨率。在實驗條件下，算法一在角度間隔大于14°左右時可有效進(jìn)行DOA 估計（RMSE ＜2 °），算法二在角度間隔大于9°左右時可有效進(jìn)行DOA 估計。

4.3 實驗三：相干、非相關(guān)混合信源時的性能比較

參數(shù)設(shè)定：信源數(shù)K=3，DOA 為(70°,80°,100°)，其中前兩個信源相干，第三個信源與前兩個信源不相關(guān)；快拍數(shù)N=200；獨立實驗次數(shù)n=500；算法一的空間平滑子陣列的陣元個數(shù)取4。圖3 為兩種算法的RMSE 與SNR 的關(guān)系：圖中藍(lán)色曲線表示算法一的RMSE 隨SNR 的變化關(guān)系；紅色曲線表示算法二的RMSE 隨SNR的變化關(guān)系。

從圖3 可以看出，對于相干、非相關(guān)混合源，算法二較算法一有更高的估計精度，在SNR 取值范圍為5～12dB 左右時，算法一基本失效，而算法二仍可以有效估計DOA；在SNR 較高時，兩種算法均可有效估計相干、非相關(guān)混合源DOA，但算法二的精度較高。

圖3 相干、非相關(guān)混合源情況下，兩種算法RMSE 隨SNR 變化曲線Fig.3 RMSE of two kinds of algorithm with various SNR for coherent and non-correlation sources

5 結(jié)束語

本文提出的算法可以有效進(jìn)行相干源DOA 估計，在已知信源數(shù)的情況下，無需進(jìn)行解相干預(yù)處理，且適用于任意陣列。而傳統(tǒng)的相干DOA 估計算法對于非均勻陣列，需要進(jìn)行虛擬變換后再進(jìn)行空間平滑，這樣會引入變換誤差，且空間平滑會損失陣列孔徑。仿真結(jié)果表明，相較于基于虛擬內(nèi)插變換的空間平滑MUSIC 算法，本文算法具有更優(yōu)的精度與分辨率，這與理論分析一致。但是，傳統(tǒng)DOA 估計算法的空間譜均為一維搜索，而本文算法的空間譜是K維搜索，因此本文算法的復(fù)雜度較高。后續(xù)工作將聚焦于如何降低本文算法的復(fù)雜度，可以考慮與其他算法如空域濾波[16]或STAP[17]相結(jié)合。