劉亞寧，張 秦，鄭桂妹

(空軍工程大學(xué)防空反導(dǎo)學(xué)院，陜西 西安 710051)

0 引言

近些年，二維DOA估計成為了研究陣列信號處理的熱點，在通信，雷達，探測等領(lǐng)域得到廣泛的應(yīng)用。二維DOA估計一般采用的陣列信號模型包括L型陣列，均勻平面陣列和平行陣列[1-3]，但也有一些特殊情況，如文獻[4]采用共形陣列，將二維DOA與極化狀態(tài)進行聯(lián)合估計。文獻[5]在稀疏陣列的基礎(chǔ)上，建立基于不動點迭代的空間譜函數(shù)，進行二維DOA估計。經(jīng)典的二維DOA估計算法包括多重信號分類(MUSIC)[6]，旋轉(zhuǎn)不變子空間(ESPRIT)[7]，傳播算子(PM)[8]，Capon算法[9]等。但這些經(jīng)典算法大多存在穩(wěn)定性差[10]，計算量復(fù)雜[11]的不足，因此人們在這些傳統(tǒng)方法的基礎(chǔ)上也提出了不少改進算法。文獻[12]提出了一種基于特征空間的MUSIC算法來最大限度的利用信號以及噪聲子空間。文獻[13]采用Unitary-ESPRIT簡化了運算復(fù)雜度。閆鋒剛等[14]提出了一種半實值Capon算法來降低運算量。

以上所提算法在入射信號相互獨立的條件下可以進行有效的估計，但在實際情況下，由于傳播環(huán)境的復(fù)雜性，直射信號和多徑傳播得到的反射信號會形成相干信號。信號的相干性會導(dǎo)致協(xié)方差矩陣不滿秩，上文所提算法的估計精度將明顯下降甚至失效。為了解決這一問題，學(xué)者們提出了不少對應(yīng)的解相干算法。常見方法有空間平滑算法以及toeplitz矩陣重構(gòu)法?？臻g平滑算法應(yīng)用最為廣泛，包括前向空間平滑算法(FSS)和前后向空間平滑算法(FBSS),如文獻[12]提出了一種基于特征空間MUSIC的空間平滑算法，文獻[15]實現(xiàn)了單快拍下利用空間平滑算法解相干。toeplitz矩陣重構(gòu)法利用方向矢量的特性，對其線性組合進行toeplitz重排，從而達到解相干的目的，計算精度較高。

本文為了解決相干信號二維DOA估計采用傳統(tǒng)方法計算量過大的問題，提出了基于降維Capon的相干信號二維DOA估計算法。

1 降維Capon算法

1.1 信號模型

根據(jù)文獻可知在常見陣列結(jié)構(gòu)中，L型陣列估計性能最優(yōu)。因此，本文在L型陣列基礎(chǔ)上進行相干信號二維DOA估計。L型陣列在x-y平面上，由分別位于x軸和y軸上的N個陣元構(gòu)成，陣元間隔為d。設(shè)遠場空間有k個窄帶信號源以不同二維波達方向入射，分別為(θ1,φ1)，(θ2,φ2)，…，(θk,φk)。其中，θk和φk分別表示第k個信號源的仰角和方位角。且θk∈[-90°,90°]，φk∈[-90°,90°]，αk和βk分別表示第k個信號源與x軸和y軸夾角。信號陣列模型如圖1所示。

圖1 信號陣列模型Fig.1 Signal array model

兩個信號的相關(guān)系數(shù)可表示為：

(1)

當ρij=0時，兩信號相互獨立；當0<|ρij|<1時，兩信號相關(guān)；當|ρij|=1時，兩信號為相干信號。若信號相干，則x軸和y軸上陣元接收的信號模型為：

X(t)=Ax(θ,φ)S(t)+Nx(t)

(2)

Y(t)=Ay(θ,φ)S(t)+Ny(t)

(3)

式(2)和式(3)中，X(t)和Y(t)分別為x軸和y軸接收信號，Ax(θ,φ)=[ax(θ1,φ1),ax(θ2,φ2),…,ax(θk,φk)]為x軸陣元方向矩陣，ax(θk,φk)=(1,exp(j2πd/λ)cosθksinφk,…,exp(j2π(N-1)d/λ)·cosθksinφk)T，Nx(t)=[n1(t),n2(t),…,nN(t)]為x軸陣元接收到的噪聲矢量，Ay(θ,φ)=[ay(θ1,φ1),ay(θ2,φ2),…,ay(θk,φk)]為y軸陣元方向矩陣，ay(θk,φk)=(1,exp(j2πd/λ)sinθksinφk,…,exp(j2π (N-1)d/λ) sinθksinφk)T，Nx(t)=[n1(t),n2(t),…,nN(t)]為y軸陣元接收到的噪聲矢量，S(t)=[λ1,λ2,…,λk]s(t)表示相干信號(λk為復(fù)常數(shù)，s(t)為源信號)。陣元接收噪聲為均值為0，方差為σ2的高斯白噪聲，并與接收信號獨立。

1.2 解相干處理

相干信號的協(xié)方差矩陣不是滿秩矩陣，故不能直接采用常規(guī)二維DOA算法進行估計，本文采用前后向空間平滑算法進行解相干，恢復(fù)協(xié)方差矩陣的秩，再進行二維DOA估計。首先對x軸陣元進行解相干,將x軸看做由N個陣元組成的均勻線陣，將其劃分為Q個子陣，每個子陣有M個陣元，M=N-Q+1。以最左邊子陣為參考，則第q個前向子陣的輸出為：

AxmDq-1S(t)+Nq(t)

(4)

式(4)中，Axm為子陣對應(yīng)方向矩陣，D=diag[γ1,γ2,…γk]，γk=exp((j2πd/λ)sinαk)，Nq(t)為子陣對應(yīng)噪聲矢量。則此子陣的協(xié)方差矩陣為：

(5)

則前向空間平滑的協(xié)方差矩陣為：

(6)

同理，第q個后向子陣的輸出為：

AxmDq-1S(t)+Nq(t)

(7)

式(7)中，Axm為子陣對應(yīng)方向矩陣，D=diag[γ1,γ2,…，γk]，γk=exp((j2πd/λ)sinαk)，Nq(t)為子陣對應(yīng)噪聲矢量。則此子陣的協(xié)方差矩陣為：

(8)

后向空間平滑的協(xié)方差矩陣為：

(9)

(10)

同理，對y軸陣元進行處理后得到的協(xié)方差矩陣為：

(11)

因此，L型陣列解相干后的信號協(xié)方差矩陣為：

(12)

1.3 降維Capon算法

一般情況下二維Capon算法的空間譜為：

(13)

由圖1可推得：

(14)

由式(14)可得：

(15)

因此可得：

ax(θ,φ)=ax(α)=(1,exp(j2πd/λ)cosα,…, exp(j2π (M-1)d/λ)cosβ)T

(16)

ay(θ,φ)=ay(β)=(1,exp(j2πd/λ)cosβ,…, exp(j2π(M-1)d/λ)cosβ)T

(17)

式(13)可表示為：

(18)

根據(jù)Kronecker積性質(zhì)可將ax(α)?ay(β)轉(zhuǎn)換成[ax(α)?IM]ay(β)，則式(18)可表示為：

(19)

由于ax(α)∈CM，ay(β)∈CM，Rfb∈CM×M，Capon的運算均為復(fù)值運算，一次復(fù)值運算包括四次實值運算。如果只進行一次實值運算便可完成DOA估計，則計算量將減少75%，因此本文參考文獻方法引入一種半實值Capon算法(SRV-Capon)來完成二維DOA估計。

對協(xié)方差矩陣Rfb進行特征值分解可得：

Rfb=SUSSH+NUNNH

(19)

式(19)中，信號子空間S=[u1,u2,…,uk]，噪聲子空間N=[v1,v2,…,vM-k]，uk為Rfb中較大的k個特征值對應(yīng)的特征向量，vM-k為Rfb中較小的M-k個特征值對應(yīng)的特征向量。US，UN均為特征值矩陣，其主對角線分別為k個大特征值ξ1,ξ2,…,ξk，和M-k個小特征值。

因為信號子空間和噪聲子空間互補正交，因此容易證明：

(Rfb)-1=S(US)-1SH+σ-2NNH=σ-2(SCSNRSH+NNH)

(20)

式(20)中，CSNR=σ2(US)-1。

(21)

(Rfb)-1≈σ-2NNH

(22)

根據(jù)二維MUSIC空間譜函數(shù)和式(18)、式(22)可得：

f2D-cap(α,β)≈σ-2f2D-music(α,β)

(23)

利用Woodbury公式[12](U+V)-1=U-1-V-1·(I+VU-1)-1VU-1可得：

Re-1(Rfb) =2(Rfb+(Rfb)*)-1= 2(Rfb)-1-2(Rfb)-1(I+(Rfb)*·

(Rfb)-1)-1(Rfb)*(Rfb)-1

(24)

將式(22)代入式(24)得：

Re-1(Rfb) ≈2σ-2NNH-2σ-2HNNH=

2σ-2(I-H)NNH

(25)

式(25)中，H=(Rfb)-1(I+(Rfb)*(Rfb)-1)-1(Rfb)*。

根據(jù)矩陣性質(zhì)可得

Re-1(Rfb)∈R(N)

(26)

將式(24)中Rfb和(Rfb)*互換位置，再按照式(24)和式(25)推導(dǎo)可得：

Re-1(Rfb)∈R(N*)

(27)

結(jié)合式(26)和式(27)可得：

Re-1(Rfb)∈R(N)∩R(N*)

(28)

定義式(19)分母為Φ:

Φ=ay(β)H[ax(α)?IM]H·

(Rfb)-1[ax(α)?IM]ay(β)

(29)

此算法便轉(zhuǎn)化為求Φ最小值時對應(yīng)的α，β值。設(shè)G(α)=[ax(α)?IM]H(Rfb)-1[ax(α)?IM]。再將(Rfb)-1替換成Re-1(Rfb)，因為ax(α)∈CM，ay(β)∈CM，Re-1(Rfb)∈RM×M，所以此運算可看作半實值處理，可得P(α)=[ax(α)?IM]HRe-1(Rfb) [ax(α)?IM]，易得ay(β)Hay(β)=M，則將(29)轉(zhuǎn)化為：

(30)

利用Rayleigh-Ritz理論[16]可得，式(30)的最小值為M×λmin(α)，λmin(α)為P(α)的最小特征值，可表示為：

(31)

因此α的估計值可由下式推得：

(32)

由式(28)和式(32)可得：在α的對稱角度-α上，也可取到式(30)的最小值，因此，我們可以把搜索范圍縮小到原來的一半，但要通過比較G(α)和G(-α)的大小來判斷哪一個為正確的α估計值。

其判斷標準為：

1)若G(α)?G(-α)，則α為正確估計值；

2)若G(-α)?G(α)，則-α為正確估計值；

3)若G(-α)≈G(α)，則α和-α均為正確估計值。

ay(βk)=emin(P(αk))

(33)

定義φk=angle(ay(βk))，結(jié)合式(17)，可得：

φk=[0,(2πd/λ)cosβk,…,(2π(M-1)d/λ)cosβk]T

(34)

利用式(25)估計βk可采用最小二乘法[17]，令：

(35)

即求解min‖Cdk-φk‖2，可得其結(jié)果為：

(36)

式(36)中，dk(2)為結(jié)果值，則βk可由下式得：

βk=arcsin(dk(2)/(2πd/λ))

(37)

2 相干信號二維DOA估計算法

2.1 算法流程

本文所提基于降維Capon的相干信號二維DOA估計算法流程為：

步驟1)先采用FBSS解相干，得到陣元數(shù)為M的子陣的滿秩協(xié)方差矩陣。

步驟2)提取子陣協(xié)方差矩陣的實部，代入一種半實值Capon算法進行DOA估計。

步驟3)利用Kronecker積的性質(zhì)求出P(α)的最小特征值，代入式(32)得到α估計值，再根據(jù)G(α)的判斷標準篩選出正確的α估計值。

步驟4)根據(jù)式(33)—式(37)，得到β估計值。

步驟5)按照式(15)得到二維DOA估計值。

2.2 計算量分析

按照上述步驟分析本文算法的計算量，并與常規(guī)二維Capon算法(2D-Capon)進行對比，設(shè)快拍數(shù)為L，搜索點數(shù)為n,子陣個數(shù)為M。

步驟1)解相干求滿秩協(xié)方差矩陣的計算量為O{(N-M)L(2M+3)}。

步驟2)求最小特征值以及篩選正確估計值的計算量為O{n[(M2+M)-K(M+1)]+5M2K}。

步驟3)求β估計值的計算量為O{10M}。

步驟4)計算量可以忽略不計。

綜上可得總計算量為O{(N-M)L(2M+3)+n[(M2+M)-K(M+1)]+5M2K+10M}。而常規(guī)二維Capon的計算量為O{(N-M)L(2M+3)+4n(M2+M)+4M2K+10M}。

為直觀表示，給出具體數(shù)值作圖2來對比兩種算法的計算量大小。設(shè)N-M=2，快拍L=100，信號個數(shù)K為3，搜索點數(shù)n=3×103，由圖2可得，本文算法相比常規(guī)二維Capon算法計算量顯著減小，可比常規(guī)算法減少接近75%。

圖2 兩種算法計算量對比圖Fig.2 Comparison chart of two algorithms’ computation complexity

3 仿真驗證

3.1 算法的正確性驗證

(38)

式(38)中，Eθ為方位角的標準誤差，Eφ為仰角的標準誤差。

為了驗證本算法，進行如下仿真。假設(shè)陣列為平面陣列，有3個窄帶信號入射到陣列上，其入射角度分別為(15°，10°)，(25°，35°)，(40°，50°)。x軸，y軸方向陣元個數(shù)均為8，子陣個數(shù)為6，陣元間隔d=λ/2，角度搜索精度為0.01°?？炫臄?shù)L為50，信噪比SNR為20 dB，進行500次蒙特卡羅實驗。圖3目標二維DOA估計值顯示了三個信源條件下本算法的估計結(jié)果。從圖中可看出，對三個目標的入射角度估計基本保持在真實值周圍，波動較小，表明本算法能正確對二維DOA進行估計，并且精度較高。

圖3 目標二維DOA估計值Fig.3 Estimation of target 2D DOA

3.2 信噪比與估計性能的關(guān)系

采用同4.1節(jié)相同的實驗驗參數(shù)，蒙特卡羅實驗次數(shù)為500，快拍數(shù)為50，x軸，y軸方向陣元個數(shù)均為8，子陣個數(shù)為6，信噪比分別取-4 dB，0 dB，4 dB，8 dB，12 dB，16 dB，20 dB。將本文算法與2D-Capon和常規(guī)二維MUSIC算法(2D-MUSIC)進行比較，結(jié)果如圖4所示。

圖4 RMSE與SNR關(guān)系圖Fig.4 The relationship between RMSE and SNR

由圖3可得，本文算法和2D-Capon和2D-MUSIC算法在信噪比較小的時候有一定差距，但隨著噪比增大，差距差距逐漸縮小，估計精度基本相同，且兩個角度的估計值RMSE隨信噪比的增大逐漸減小，說明估計性能逐漸變好，但減小的趨勢越來越緩。

3.3 快拍數(shù)與估計性能的關(guān)系

采用同4.1節(jié)相同的實驗參數(shù)，蒙特卡羅實驗次數(shù)為500，x軸，y軸方向陣元個數(shù)均為8，子陣個數(shù)為6，信噪比取20 dB，快拍數(shù)分別為10，20，40，60，80，100，120，140，160。將本文算法與2D-Capon和2D-MUSIC進行比較，結(jié)果如圖5所示。

圖5 RMSE與快拍數(shù)關(guān)系圖Fig.5 The relationship between RMSE and snapshot

由圖5可得，當快拍數(shù)小于40的時候，本文算法RMSE小于2D-Capon，本文給出的分析是在小快拍下，信號相關(guān)性較大，復(fù)數(shù)域范圍的協(xié)方差矩陣相比于實數(shù)域范圍的協(xié)方差矩陣秩虧更大，因此有可能導(dǎo)致2D-Capon估計性能不如SRV-Capon。當快拍數(shù)較大時，本文算法相比2D-Capon和2D-MUSIC算法基本相同。隨著快拍數(shù)的增大，RMSE值減小，但RMSE值減小的趨勢越來越緩。說明雖然估計性能越來越好，但變好的趨勢是越來越緩。由于快拍數(shù)的選擇影響算法的計算復(fù)雜度，因此并非快拍數(shù)越大越好。

4 結(jié)論

本文提出了基于降維Capon的相干信號二維DOA估計算法，該算法采用空間平滑和半實值降維Capon算法，相比常規(guī)二維Capon算法，計算量顯著減小，仿真結(jié)果表明本文算法的估計性能與常規(guī)二維Capon算法和二維MUSIC算法基本相同。