基于降維Capon的相干信號二維DOA估計(jì)算法

劉亞寧，張 秦，鄭桂妹

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

0 引言

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

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

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

1 降維Capon算法

1.1 信號模型

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

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

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

(1)

當(dāng)ρij=0時(shí)，兩信號相互獨(dú)立；當(dāng)0<|ρij|<1時(shí)，兩信號相關(guān)；當(dāng)|ρij|=1時(shí)，兩信號為相干信號。若信號相干，則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的高斯白噪聲，并與接收信號獨(dú)立。

1.2 解相干處理

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

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個(gè)后向子陣的輸出為：

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軸陣元進(jìn)行處理后得到的協(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的運(yùn)算均為復(fù)值運(yùn)算，一次復(fù)值運(yùn)算包括四次實(shí)值運(yùn)算。如果只進(jìn)行一次實(shí)值運(yùn)算便可完成DOA估計(jì)，則計(jì)算量將減少75%，因此本文參考文獻(xiàn)方法引入一種半實(shí)值Capon算法(SRV-Capon)來完成二維DOA估計(jì)。

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

Rfb=SUSSH+NUNNH

(19)

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

因?yàn)樾盘栕涌臻g和噪聲子空間互補(bǔ)正交，因此容易證明：

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

(30)

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

(31)

因此α的估計(jì)值可由下式推得：

(32)

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

其判斷標(biāo)準(zhǔn)為：

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

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

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

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)估計(jì)βk可采用最小二乘法[17]，令：

(35)

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

(36)

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

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

(37)

2 相干信號二維DOA估計(jì)算法

2.1 算法流程

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

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

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

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

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

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

2.2 計(jì)算量分析

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

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

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

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

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

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

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

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

3 仿真驗(yàn)證

3.1 算法的正確性驗(yàn)證

(38)

式(38)中，Eθ為方位角的標(biāo)準(zhǔn)誤差，Eφ為仰角的標(biāo)準(zhǔn)誤差。

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

圖3 目標(biāo)二維DOA估計(jì)值Fig.3 Estimation of target 2D DOA

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

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

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

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

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

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

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

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

4 結(jié)論

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