宗佳虹，戴旭初

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

引 言

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

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

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

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

1 信號(hào)模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

輸出：各信源的DOA 估計(jì)。

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

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

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

④ 利用MUSIC 算法得到DOA 估計(jì)。

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

3 本文算法

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

3.1 算法設(shè)計(jì)

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

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

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

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

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

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

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

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

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

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

3.1.4 基于空間譜搜索的DOA 估計(jì)

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

最后通過峰值搜索，得到DOA 估計(jì)。

3.1.5 算法流程

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

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

輸出：各信源的DOA 估計(jì)。

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

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

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

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

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

3.2 復(fù)雜度分析

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

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

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

4 仿真實(shí)驗(yàn)及分析

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

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

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

4.1 實(shí)驗(yàn)一：相干源的估計(jì)精度比較

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

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

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

4.2 實(shí)驗(yàn)二：相干源的分辨率性能比較

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

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

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

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

從圖3 可以看出，對(duì)于相干、非相關(guān)混合源，算法二較算法一有更高的估計(jì)精度，在SNR 取值范圍為5～12dB 左右時(shí)，算法一基本失效，而算法二仍可以有效估計(jì)DOA；在SNR 較高時(shí)，兩種算法均可有效估計(jì)相干、非相關(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 估計(jì)，在已知信源數(shù)的情況下，無需進(jìn)行解相干預(yù)處理，且適用于任意陣列。而傳統(tǒng)的相干DOA 估計(jì)算法對(duì)于非均勻陣列，需要進(jìn)行虛擬變換后再進(jìn)行空間平滑，這樣會(huì)引入變換誤差，且空間平滑會(huì)損失陣列孔徑。仿真結(jié)果表明，相較于基于虛擬內(nèi)插變換的空間平滑MUSIC 算法，本文算法具有更優(yōu)的精度與分辨率，這與理論分析一致。但是，傳統(tǒng)DOA 估計(jì)算法的空間譜均為一維搜索，而本文算法的空間譜是K維搜索，因此本文算法的復(fù)雜度較高。后續(xù)工作將聚焦于如何降低本文算法的復(fù)雜度，可以考慮與其他算法如空域?yàn)V波[16]或STAP[17]相結(jié)合。