晏 輝, 司偉建, 夏新凡

(1.哈爾濱工程大學(xué)信息與通信工程學(xué)院,黑龍江 哈爾濱 150001;2.哈爾濱工程大學(xué)先進(jìn)船舶通信與信息技術(shù)工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室,黑龍江 哈爾濱 150001;3.上海無(wú)線電設(shè)備研究所,上海 201109)

0 引言

波達(dá)方向(direction of arrival,DOA)估計(jì)作為陣列信號(hào)處理的關(guān)鍵技術(shù),已經(jīng)廣泛應(yīng)用于通信、聲納、雷達(dá)等領(lǐng)域[1]。均勻圓陣(uniform circular array,UCA)作為測(cè)向系統(tǒng)中常用的陣列之一,因其具備360°的方位角測(cè)向能力、容易共形和導(dǎo)向矢量共軛對(duì)稱等優(yōu)勢(shì)受到廣泛關(guān)注[2]。針對(duì)其復(fù)雜的導(dǎo)向矢量,文獻(xiàn)[3]提出模式空間變換的概念,將圓陣變?yōu)樘摂M均勻線陣(uniform linear array,ULA),使其導(dǎo)向矢量具備范德蒙結(jié)構(gòu)。文獻(xiàn)[4]通過(guò)構(gòu)造實(shí)值波束空間變換矩陣,結(jié)合多重信號(hào)分類(lèi)(MUSIC)算法和旋轉(zhuǎn)不變子空間(ESPRIT)算法提出基于均勻圓陣的實(shí)值波束空間多重信號(hào)分類(lèi)(UCA-RB-MUSIC)算法和基于均勻圓陣的旋轉(zhuǎn)不變子空間(UCA-ESPRIT)算法。但是波束空間變換會(huì)引入估計(jì)誤差。文獻(xiàn)[5-6]通過(guò)對(duì)誤差定量分析,采用迭代的方法消除誤差主要部分,算法性能得以提高,但是計(jì)算量也隨之增大。在陣列稀疏情況下,文獻(xiàn)[7]針對(duì)波束變換帶來(lái)的誤差問(wèn)題,對(duì)導(dǎo)向矢量進(jìn)行誤差補(bǔ)償,減少了陣元數(shù)需求。文獻(xiàn)[8]利用陣列流形分離技術(shù),結(jié)合傳播因子(PM)和求根MUSIC算法估計(jì)角度信息,但是會(huì)引入映射誤差。文獻(xiàn)[9]引入互質(zhì)稀疏陣列,不僅能夠減小互耦誤差,還能估計(jì)角度間隔很小的信源。但上述三種稀疏陣列算法都只能估計(jì)一維角度信息,無(wú)法估計(jì)二維角度信息。文獻(xiàn)[10-13]利用雙圓陣估計(jì)二維角度信息,但所需陣元數(shù)較多且沒(méi)有考慮陣列稀疏的情況。文獻(xiàn)[14]在稀疏陣列波束域誤差補(bǔ)償?shù)幕A(chǔ)上進(jìn)行相位校正,可以估計(jì)二維相干信號(hào),但是估計(jì)的信源數(shù)有限。

本文采用稀疏雙圓陣列,首先根據(jù)兩個(gè)子陣的旋轉(zhuǎn)不變性,用子陣接收數(shù)據(jù)協(xié)方差矩陣構(gòu)造波達(dá)矩陣,估計(jì)信號(hào)的俯仰角[15];其次利用波束變換理論,將UCA變?yōu)樘摂M的ULA,使導(dǎo)向矢量具備范德蒙結(jié)構(gòu);再次對(duì)稀疏陣列的導(dǎo)向矢量進(jìn)行誤差補(bǔ)償,減少稀疏陣元下波束空間變換帶來(lái)的誤差影響;最后利用求根MUSIC算法估計(jì)信號(hào)的方位角。所提方法可以應(yīng)用于陣列稀疏的實(shí)際工程場(chǎng)景中。

1 陣列結(jié)構(gòu)和信號(hào)模型

假設(shè)K個(gè)波長(zhǎng)為λ的遠(yuǎn)場(chǎng)窄帶非相干信號(hào)入射到各向同性的雙圓形稀疏陣列,如圖1所示。陣列子陣陣元數(shù)為M,相鄰陣元間距大于λ/2,陣列半徑為r,兩個(gè)子陣間距d=λ/2,入射信號(hào)俯仰角φ∈[0,π],方位角θ∈ [0,2π]。

圖1 雙圓形稀疏陣列結(jié)構(gòu)圖

在不考慮陣列誤差、通道不一致和互耦影響的前提下,子陣1和子陣2的接收信號(hào)矢量表達(dá)式分別為

式中:A為M×K維陣列導(dǎo)向矢量矩陣;S(t)=[s1(t),s2(t),...,sK(t)]T為K×1維入射信號(hào)矢量,T 表示轉(zhuǎn)置運(yùn)算;n(t)=[n0(t),n1(t),…,nM-1(t)]T為M×1維噪聲矩陣;Φ為相位差矩陣。本文假設(shè)噪聲為加性高斯白噪聲,彼此獨(dú)立且與信號(hào)不相關(guān),噪聲功率為δ2。導(dǎo)向矢量矩陣A可以表示為

其中

式中:a(θi,φi)為信號(hào)i的導(dǎo)向矢量,其中,φi和θi分別為信號(hào)i的俯仰角和方位角;k=2π/λ為信號(hào)波數(shù);γn=2πn/M為陣元n的位置。

相位差矩陣Φ可表示為

式中:diag(·)為矩陣對(duì)角化函數(shù)。

2 基于稀疏陣列的二維DOA估計(jì)

2.1 俯仰角估計(jì)

文獻(xiàn)[15]首次提出波達(dá)方向矩陣法,用于雙線性平行陣列二維參量估計(jì),文獻(xiàn)[10]將其擴(kuò)展到雙圓陣中。該方法首先根據(jù)接收信號(hào)矢量得到兩個(gè)子陣的自協(xié)方差矩陣RXX和互協(xié)方差矩陣RYX,表達(dá)式為

其中

定義波達(dá)方向矩陣

文獻(xiàn)[15]已證明,當(dāng)導(dǎo)向矢量矩陣A和信號(hào)協(xié)方差矩陣RSS滿秩時(shí),波達(dá)方向矩陣R滿足

即通過(guò)對(duì)R的特征分解得到特征值和特征向量,其中非零特征值和對(duì)應(yīng)的特征向量分別與對(duì)角陣Φ和導(dǎo)向矢量A對(duì)應(yīng)相等,通過(guò)此對(duì)等關(guān)系即可求出俯仰角和方位角。本文采用文獻(xiàn)[10]中的方法對(duì)俯仰角求解,計(jì)算公式為

式中:Arg(·)為求取復(fù)數(shù)輻角函數(shù);ηi為R特征分解得到的非零特征值。但是在求解方位角時(shí),首先特征分解操作并不能保證特征向量唯一(存在比例關(guān)系),導(dǎo)向矢量矩陣A的對(duì)等關(guān)系不一定成立;其次由于圓陣方位角θ∈[0,2π],容易出現(xiàn)測(cè)角模糊,想要得到準(zhǔn)確的角度信息必須進(jìn)行額外處理,增大了計(jì)算難度。

2.2 方位角估計(jì)

針對(duì)2.1節(jié)方位角估計(jì)所遇到的問(wèn)題,引入稀疏陣列下的波束空間變換法。首先利用波束空間變換,將圓陣轉(zhuǎn)化為導(dǎo)向矢量具備范德蒙結(jié)構(gòu)的虛擬線陣;然后進(jìn)行導(dǎo)向矢量補(bǔ)償,消除稀疏陣列情況下的誤差影響;最后通過(guò)求根MUSIC算法估計(jì)出信號(hào)的方位角。

在得到俯仰角后,利用波束空間變換理論[3],直接用波束變換矩陣W乘以去噪后的接收信號(hào)協(xié)方差矩陣RXX,0,使UCA的導(dǎo)向矢量具備范德蒙結(jié)構(gòu)。波束變換矩陣的表達(dá)式為

其中

式中:ε=ceil(kr)為最大相位模式數(shù),其中ceil(·)表示向上取整。當(dāng)M＞2ε+1時(shí),波束域協(xié)方差矩陣

其中

式 中:Γ(φi)=diag(V-ε(φi),…,Vm(φi),…,Vε(φi))為波束域?qū)蚴噶康母┭鼋遣糠?其中Vm(φi)=jmJm(kr sinφi),Jm·（）表示階數(shù)為m的第一類(lèi)貝塞爾函數(shù);av(θi)為虛擬線陣導(dǎo)向矢量。由此看出經(jīng)過(guò)波束空間變換后,導(dǎo)向矢量中的方位角和俯仰角已經(jīng)分離。

當(dāng)M＜2ε+1(陣元間距大于λ/2)時(shí),陣列變?yōu)橄∈桕嚵?波束空間變換帶來(lái)的殘差不可忽略。為了消除誤差影響,本文采用文獻(xiàn)[7]中的方法,重新定義稀疏下的陣元數(shù)為2N+1(N＜ε),則新的波束變換矩陣可表示為

波束域虛擬線陣長(zhǎng)度變?yōu)?N+1,稀疏陣列第i個(gè)信號(hào)源的導(dǎo)向矢量可表示為

其中

式中:a0(θi,φi)為稀疏波束域中不含誤差的導(dǎo)向矢量部分;Δa(θi,φi)為相位模式數(shù)在 [N,ε]之間的波束變換估計(jì)誤差;Γs(φi)為稀疏陣列導(dǎo)向矢量的俯仰角部分;avs(θi)為稀疏情況下虛擬線陣的導(dǎo)向矢量。新的波束域協(xié)方差矩陣表達(dá)式為

式中:I為(2N+1)×(2N+1)維單位矩陣;Jl是單位矩陣I的后ε-N列(2N+1)×(ε-N)維矩陣(ε≤3N+1);Jr是單位矩陣I的前ε-N列(2N+1)×(ε-N)維矩陣。則補(bǔ)償后的波束域?qū)蚴噶?/p>

式中:Δa′(θi,φi)為補(bǔ)償后的殘差項(xiàng),可以通過(guò)選擇合適的模式數(shù)達(dá)到任意小值。誤差消除后利用求根MUSIC算法可以估計(jì)出方位角,第i個(gè)信號(hào)的求根多項(xiàng)式為

則本文所提方法的二維DOA估計(jì)具體步驟為:

a)獲取兩個(gè)陣列的接收數(shù)據(jù)矢量X(t),Y(t);

b)分別求取子陣1的自協(xié)方差矩陣RXX和兩個(gè)子陣的互協(xié)方差矩陣RYX,然后根據(jù)式(11)構(gòu)造波達(dá)方向矩陣R;

c)對(duì)R進(jìn)行特征分解,然后利用式(13)估計(jì)俯仰角φi;

f)根據(jù)式(24),利用求根MUSIC算法估計(jì)方位角θi。

3 仿真實(shí)驗(yàn)

為了驗(yàn)證所提方法的正確性,做如下仿真實(shí)驗(yàn):實(shí)驗(yàn)1驗(yàn)證本文方法的有效性;實(shí)驗(yàn)2驗(yàn)證本文方法隨信噪比變化的情況,對(duì)比算法為文獻(xiàn)[4]中的UCA-ESPRIT算法和UCA-RB-MUSIC算法;實(shí)驗(yàn)3驗(yàn)證高度相關(guān)情況下信噪比對(duì)本文方法的影響;實(shí)驗(yàn)4驗(yàn)證快拍數(shù)對(duì)本文方法的影響。在構(gòu)造補(bǔ)償矩陣P時(shí),本文通過(guò)大量仿真實(shí)驗(yàn)證明,當(dāng)ε=3N+1時(shí),估計(jì)性能較好。所以本文所有仿真實(shí)驗(yàn)均取N=(ε-1)/3。

(1)實(shí)驗(yàn)1

假設(shè)4個(gè)獨(dú)立信號(hào)分別以入射角(100°,35°),(110°,38°),(190°,40°),(220°,45°)入射到圖1所示陣列。取r=1.5λ,則N=3,M=7,信噪比為25 dB。仿真的快拍數(shù)為128,獨(dú)立進(jìn)行100次蒙特卡羅實(shí)驗(yàn)。仿真結(jié)果如圖2所示。對(duì)于間隔較近的信號(hào),本文所提方法能夠準(zhǔn)確估計(jì)出結(jié)果。

圖2 4個(gè)信號(hào)二維DOA估計(jì)結(jié)果

針對(duì)3個(gè)信源情況,可以進(jìn)一步減少陣元數(shù)。取r=λ,則N=2,M=5,其他條件不變,得到入射角分別為(50°,20°),(100°,80°),(150°,50°)的3個(gè)獨(dú)立信號(hào)的估計(jì)結(jié)果,如圖3所示。

圖3 3個(gè)信號(hào)二維DOA估計(jì)結(jié)果

(2)實(shí)驗(yàn)2

假設(shè)4個(gè)獨(dú)立信號(hào)分別以入射角(60°,20°),(120°,30°),(32°,35°),(200°,40°)入射到圖1所示陣列,其他仿真條件與實(shí)驗(yàn)1相同。為了定量分析,本文算法陣元數(shù)M取7,其他兩種算法陣元數(shù)均取14。均方根誤差定義為

式中:n為蒙特卡羅實(shí)驗(yàn)次數(shù);(θi,φi)為信號(hào)i入射角的真實(shí)值;(θij,φij)為信號(hào)i入射角的第j次估計(jì)值。

信噪比從5 dB增加到29 dB,步長(zhǎng)3 d B,進(jìn)行100次蒙特卡羅實(shí)驗(yàn),仿真得到本文算法、文獻(xiàn)[4]中UCA-ESPRIT算法和UCA-RB-MUSIC算法估計(jì)結(jié)果的均方根誤差隨信噪比變化情況,如圖4所示。

圖4 實(shí)驗(yàn)2均方根誤差隨信噪比變化情況

從圖4中可以看出,本文算法性能介于UCA-RB-MUSIC算法和UCA-ESPRIT算法之間。相比二維譜峰搜索的UCA-RB-MUSIC算法,本文算法計(jì)算量大大減少,運(yùn)行時(shí)間顯著提高;相比UCA-ESPRIT算法,本文算法由于多次特征分解,計(jì)算量稍大,但是分辨率更高,估計(jì)性能更好。

(3)實(shí)驗(yàn)3

由于多徑傳播和障礙物遮擋等原因,實(shí)際環(huán)境中充斥著大量高度相關(guān)信號(hào),需驗(yàn)證高度相關(guān)情況下信噪比對(duì)本文方法的影響。假設(shè)4個(gè)相關(guān)信號(hào)分別以入射角(120°,20°),(60°,40°),(320°,60°),(200°,80°)入射到圖1所示陣列,各信號(hào)間的相關(guān)系數(shù)為0.9,其他條件與實(shí)驗(yàn)1相同。信噪比從5 d B增加到29 dB,步長(zhǎng)3 d B,進(jìn)行100次蒙特卡羅實(shí)驗(yàn),仿真得到本文算法估計(jì)結(jié)果的均方根誤差隨信噪比變化情況,如圖5所示?？芍岱椒▽?duì)于高度相關(guān)信號(hào)具備較好的分辨性能,當(dāng)信噪比高于15 dB時(shí),均方根誤差顯著減小。

圖5 實(shí)驗(yàn)3均方根誤差隨信噪比變化情況

(4)實(shí)驗(yàn)4

假設(shè)4個(gè)獨(dú)立信號(hào)分別以入射角(120°,20°),(50°,40°),(200°,60°),(270°,50°)入射到圖1所示陣列,其他仿真條件與實(shí)驗(yàn)1相同。進(jìn)行100次蒙特卡羅實(shí)驗(yàn),改變快拍數(shù),仿真得到本文方法估計(jì)結(jié)果的均方根誤差隨快拍數(shù)變化情況,如圖6所示。當(dāng)快拍數(shù)大于50時(shí),所提方法估計(jì)結(jié)果的均方根誤差已然很小,算法實(shí)時(shí)性較好,基本符合實(shí)際工程應(yīng)用中小快拍數(shù)要求。

圖6 實(shí)驗(yàn)4均方根誤差隨快拍數(shù)變化情況

4 結(jié)束語(yǔ)

本文利用稀疏雙圓陣接收數(shù)據(jù)矢量的協(xié)方差矩陣構(gòu)造波達(dá)方向矩陣,通過(guò)波達(dá)方向矩陣特征分解估計(jì)入射信號(hào)俯仰角;然后對(duì)子陣1進(jìn)行波束空間變換,使圓陣變?yōu)閷?dǎo)向矢量具備范德蒙結(jié)構(gòu)的虛擬線陣;再對(duì)稀疏情況下的導(dǎo)向矢量進(jìn)行補(bǔ)償,消除了波束空間變換引起的誤差;最后利用求根MUSIC算法估計(jì)入射信號(hào)的方位角。所提方法無(wú)需進(jìn)行復(fù)雜的二維譜峰搜索,方位角和俯仰角自動(dòng)配對(duì),可以應(yīng)用在陣列稀疏情況的二維測(cè)向系統(tǒng)中。雖然雙圓陣列損失了一個(gè)子陣的陣元,但實(shí)際上稀疏情況下總的陣元數(shù)并不多,滿足工程實(shí)際的要求。未來(lái)主要的工作是將本文方法應(yīng)用在信號(hào)相干情況下。