項(xiàng) 楊 楊晉生

(天津大學(xué)微電子學(xué)院，天津 300072)

1 引言

波達(dá)方向(direction of arrival， DOA)估計(jì)是陣列信號(hào)處理的重要內(nèi)容之一，受到了眾多學(xué)者廣泛的關(guān)注[1]。其中，二維DOA估計(jì)算法被廣泛應(yīng)用于各種形式的陣列，例如雙平行均勻線陣[2-3]，三平行均勻線陣[4-5]，L型均勻線陣[6-13]等。L型陣列具有形式簡(jiǎn)單且能夠提供較好的角度估計(jì)性能等優(yōu)點(diǎn)。因此，學(xué)者們提出了大量基于L型陣列的二維DOA估計(jì)算法。文獻(xiàn)[8]通過(guò)構(gòu)建匹配矩陣實(shí)現(xiàn)了方位角和俯仰角的自動(dòng)配對(duì)，但它需要類似譜峰搜索的運(yùn)算。文獻(xiàn)[11]解決了角度兼并問(wèn)題，在欠定條件下可實(shí)現(xiàn)角度估計(jì)且能夠自動(dòng)配對(duì)。

上述的算法有一個(gè)共同的特點(diǎn)，即陣元間距小于等于半波長(zhǎng)。擴(kuò)展孔徑可以有效地提高陣列的分辨率和角度估計(jì)精度，但會(huì)出現(xiàn)角度模糊的現(xiàn)象。文獻(xiàn)[15]提出了一種獲得高精度無(wú)模糊方向余弦的方法，因此有效地提高了DOA估計(jì)性能。當(dāng)包含方向余弦估計(jì)信息的DOA矩陣具有相同的對(duì)角線元素時(shí)，它不能解決奇異點(diǎn)問(wèn)題，即不能獲得正確的方向余弦估計(jì)值。在此基礎(chǔ)上，文獻(xiàn)[16]提出的算法將擴(kuò)展孔徑的陣列應(yīng)用到xoy平面的兩個(gè)雙平行陣列中，且使用傳播算子算法[17]降低了計(jì)算復(fù)雜度，但仍不能解決奇異點(diǎn)問(wèn)題。針對(duì)特征值相等時(shí)的角度配對(duì)失效現(xiàn)象，文獻(xiàn)[18]提出一種基于子空間正交的新配對(duì)方法，故能夠解決部分奇異點(diǎn)問(wèn)題。

為了解決擴(kuò)展孔徑的二維DOA估計(jì)存在的奇異點(diǎn)問(wèn)題，本文提出了非均勻L型陣列的聯(lián)合對(duì)角化二維DOA估計(jì)算法。首先，從延時(shí)互相關(guān)矩陣中提取信號(hào)子空間。然后，通過(guò)聯(lián)合對(duì)角化方法獲得自動(dòng)配對(duì)的兩種方向余弦，即高精度模糊的方向余弦和低精度無(wú)模糊的方向余弦。值得注意的是，即使DOA估計(jì)矩陣具有相同的對(duì)角元素，也可以獲得正確的方向余弦。最后，為了實(shí)現(xiàn)解模糊，以無(wú)模糊估計(jì)值為參考，從模糊估計(jì)值中得到高精度無(wú)模糊的估計(jì)值。因此，提出的算法解決了在擴(kuò)展孔徑的二維DOA估計(jì)中的奇異點(diǎn)問(wèn)題，并且在欠定條件下具有良好的表現(xiàn)性能。

符號(hào):(·)T,(·)*,(·)H和(·)+分別表示轉(zhuǎn)置,共軛,共軛轉(zhuǎn)置和偽逆運(yùn)算?！押?分別表示 Khatri-Rao (KR)積和Kronecker積。E[·]表示統(tǒng)計(jì)期望，arg(·)表示相位。IM是一個(gè)維數(shù)M×M單位矩陣。diag{·}是由列向量元素組成的對(duì)角矩陣。blkdiag{·}表示塊對(duì)角化。 circshift(·,m)是沿著行向右循環(huán)移動(dòng)m個(gè)單位。

2 陣列結(jié)構(gòu)和數(shù)據(jù)模型

圖1 陣列插圖Fig.1 Illustration of array

ρε(t)=Aεs(t)+nε(t)

(1)

(2)

此外，其他兩個(gè)子陣的陣列流型矩陣如下

(3)

3 提出的算法

所提出的算法主要包括三個(gè)部分。首先，根據(jù)KR運(yùn)算構(gòu)造延時(shí)互相關(guān)矩陣。其次，在延時(shí)互相關(guān)矩陣和選擇矩陣運(yùn)算得到對(duì)角矩陣的基礎(chǔ)上，通過(guò)聯(lián)合對(duì)角化方法得到自動(dòng)配對(duì)的兩種方向余弦估計(jì)值，即高精度模糊的方向余弦估計(jì)值和低精度無(wú)模糊的方向余弦估計(jì)值。最后，通過(guò)解模糊方法得到高精度無(wú)模糊方向余弦估計(jì)值，進(jìn)一步得到方位角和俯仰角。

3.1 延時(shí)互相關(guān)矩陣的構(gòu)造

(4)

因此，根據(jù)KR運(yùn)算得到的延時(shí)互相關(guān)矩陣表示如下

(5)

(6)

(7)

(8)

根據(jù)式(8), 對(duì)應(yīng)的陣列形式如圖2所示。

圖2 陣列劃分Fig.2 Partition of array

3.2 方向余弦的估計(jì)

通過(guò)對(duì)R進(jìn)行奇異值分解(Singular Value Decomposition，SVD)，可以得到信號(hào)子空間Us以及具有K個(gè)較大奇異值的對(duì)角矩陣Λs

(9)

從式(8)和圖2(b)易知，平面陣列1與2、1與3、3與4以及2與4之間的間距均為d，且四個(gè)平面陣列內(nèi)部相鄰陣元間距為ds，故Us包含高精度模糊的方向余弦信息以及低精度無(wú)模糊的方向余弦信息。

平面陣列1與2，3與4之間均包含x軸低精度無(wú)模糊的方向余弦信息。因此，構(gòu)造選擇矩陣G1=[G01,G00,G02,G00]，用于選取與平面陣列1、3相對(duì)應(yīng)的Us中兩個(gè)子塊；G2=circshift(G1，M2)用于選取與平面陣列2、4相對(duì)應(yīng)的Us中的兩個(gè)子塊。其中，G00=02M2×M2，G01=[IM2;0M2]，G02=[0M2;IM2]。因此，包含x軸低精度無(wú)模糊的方向余弦對(duì)角矩陣可表示如下

(10)

式中，T=A+UsΛ1/2是一個(gè)酉矩陣。

(11)

四個(gè)平面陣列中均包含x軸高精度模糊的方向余弦信息，對(duì)應(yīng)對(duì)角矩陣表示如下

(12)

式中，G5=I4M?[IM-1,0(M-1)×1]，G6=circshift(G5,1)，Φ(θ)=diag{ej2πdscos θ1/λ,...,ej2πdscos θK/λ}。

(13)

式中，Φ(φ)=diag{ej2πdscos φ1/λ,...,ej2πdscos φK/λ}。

(14)

(15)

3.3 二維DOA估計(jì)的實(shí)現(xiàn)

(16)

(17)

根據(jù)以上的分析，第k個(gè)信號(hào)的方位角和俯仰角估計(jì)表達(dá)式如下

(18)

因此，我們可以得到自動(dòng)配對(duì)的俯仰角和方位角。所提出算法主要步驟概括如下：

(1)構(gòu)造式(7)所示的延時(shí)互相關(guān)矩陣。

(4)同理得到x軸上對(duì)應(yīng)的兩種方向余弦估計(jì)值。

3.4 算法分析

(1)角度估計(jì)

選擇矩陣的維數(shù)直觀地說(shuō)明了所提出的算法可實(shí)現(xiàn)多達(dá)2M2個(gè)信源估計(jì)。因此，所提出的DOA估計(jì)算法在欠定條件下仍能實(shí)現(xiàn)方位角和俯仰角的估計(jì)與自動(dòng)配對(duì)。同時(shí)，當(dāng)俯仰角和方位角滿足以下三種情況時(shí)，所提出的算法仍能獲得正確的方向余弦估計(jì):

情況 1θk=θp，cosφk≠cosφp±oλ/ds。

情況 2φk=φp，cosθk≠cosθp±oλ/ds。

情況 3θk≠θp，φk≠φp,cosθk=cosθp±oλ/ds，cosφk=cosφp±oλ/ds。

礦區(qū)生活區(qū)、選礦車間衛(wèi)生間等生活污水經(jīng)化糞池預(yù)處理后，與經(jīng)隔油池處理后的食堂含油廢水通過(guò)生活污水排水管道系統(tǒng)匯集，經(jīng)一體化生活污水處理設(shè)施處理，達(dá)到標(biāo)準(zhǔn)后回用至生產(chǎn)回水水池。

其中，o為正整數(shù)表示兩個(gè)來(lái)波信號(hào)在同一方向上的方向余弦估計(jì)值的差值為整數(shù)倍λ/ds。

證明如下：

由于R=ARss，故通過(guò)判斷A的階數(shù)即可得到R的階數(shù)。

1)首先，判斷A中互不相關(guān)行的個(gè)數(shù)。

由式(8)可知,陣列流型矩陣A中第k列中的不同行之間可能存在的關(guān)系分別為ej2πdcos φk/λ，ej2πdcos φk/λ，ej2πdscos θk/λ，ej2πdscos φk/λ。

當(dāng)角度關(guān)系滿足情況 1時(shí)，ej2πdscos θk/λ=ej2πdscos θp/λ且ej2πdcos θk/λ=ej2πdcos θp/λ，但由于cosφk≠cosφp±oλ/ds，故ej2πdscos φk/λ≠ej2πdscos φp/λ且ej2πdcos φk/λ≠ej2πdcos φp/λ，結(jié)合3.2中關(guān)于高精度模糊的方向余弦信息的選擇矩陣的維數(shù)可知A中線性無(wú)關(guān)的行數(shù)至少為2M2。

當(dāng)角度關(guān)系滿足情況 2時(shí)，與情況1同理。

當(dāng)角度關(guān)系滿足情況3時(shí)，ej2πdcos θk/λ≠ej2πdcos θp/λ，ej2πdscos φk/λ≠ej2πdscos φp/λ且ej2πdscos θk/λ=ej2πdscos θp/λ，ej2πdcos φk/λ=ej2πdcos φp/λ，結(jié)合3.2中關(guān)于低精度無(wú)模糊的方向余弦信息的選擇矩陣的維數(shù)可知A中線性無(wú)關(guān)的行數(shù)為至少為2M2。

2)其次，判斷A中互不相關(guān)列的個(gè)數(shù)。

由式(8)可知, 陣列流型矩陣A中的各列之間并不存在線性關(guān)系，除非兩個(gè)信號(hào)的方位角與俯仰角完全相同。故對(duì)于以上三種情況，A中各列始終線性無(wú)關(guān)。

綜上所述，由于2M2>K，因此通過(guò)對(duì)R進(jìn)行SVD分解可以得到信號(hào)子空間Us以及K個(gè)較大奇異值。

3)最后，文獻(xiàn)[11]指出聯(lián)合對(duì)角化可以有效處理角度兼并問(wèn)題，因此所提出的算法在以上三種情況下仍能實(shí)現(xiàn)良好的角度估計(jì)。

(2)克拉美羅界

本文所提出算法的克拉美羅界(Cramer-Rao bound，CRB)[19]表達(dá)式如下

(19)

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

本文所提出算法的角度估計(jì)性能與文獻(xiàn)[11]、文獻(xiàn)[18]中的算法以及CRB[19]進(jìn)行比較。所有信源的聯(lián)合均方誤差定義為

(20)

在所有的仿真實(shí)驗(yàn)中，所提出算法中子陣的陣元個(gè)數(shù)M=3，文獻(xiàn)[11]中子陣的陣元個(gè)數(shù)為6，即L型陣列中共有11個(gè)陣元。由于文獻(xiàn)[18]中算法應(yīng)用的陣列結(jié)構(gòu)不同，令其共包含13個(gè)陣元。此外，所有的仿真實(shí)驗(yàn)均進(jìn)行1000次蒙特卡洛仿真。

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

假設(shè)在一種欠定條件下，其中快拍數(shù)N、數(shù)據(jù)幀數(shù)L、信源個(gè)數(shù)K、信噪比SNR和擴(kuò)展孔徑ds分別為1000，500，18，30 dB，5λ。圖3顯示了欠定條件下所提出算法仍具有良好的估計(jì)性能。

如圖3所示，當(dāng)方位角和俯仰角滿足情況1和情況2時(shí)，算法仍能夠較好地實(shí)現(xiàn)方位角和俯仰角估計(jì)。當(dāng)M=3時(shí)，本文提出的算法可實(shí)現(xiàn)多達(dá)18個(gè)信源估計(jì)。

圖3 欠定條件下的角度估計(jì)性能Fig.3 Angle estimation performance in the underdetermined case

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

當(dāng)方位角和俯仰角滿足情況3時(shí)，驗(yàn)證所提出算法的有效性。假設(shè)有K=2個(gè)信號(hào)入射到天線陣列，兩個(gè)信號(hào)分別為(90°,60°)，(120°,90°)或(65°,33°)，(85°,60°)。即o1=1，ds=2λ，cos 60°=cos 90°+1/2且cos 90°=cos 120°+1/2；o2=1，ds=3λ，cos 65°=cos 85°+1/3且cos 33°=cos 60°+1/3。其他實(shí)驗(yàn)條件均與仿真實(shí)驗(yàn)1相同。圖4為角度估計(jì)值分布散點(diǎn)圖。值得注意的是文獻(xiàn)[18]中的算法不能實(shí)現(xiàn)此種情況下的角度估計(jì)。

圖4 情況3下的角度估計(jì)性能Fig.4 Angle estimation performance for Case 3

從圖4可以看出，當(dāng)方位角和俯仰角滿足情況3時(shí)，所提出的算法能夠清晰地分辨這兩個(gè)來(lái)波信號(hào)。根據(jù)文獻(xiàn)[20]中的定理1可知，式(10)～(13)所示的四個(gè)對(duì)角矩陣中沒(méi)有相等的對(duì)角元素時(shí)，算法才可獲得正確的方向余弦估計(jì)值。在擴(kuò)展孔徑的二維 DOA估計(jì)算法中隨著ds的增大，存在o，ds使得cosθ1(cosφ1)=cosθ2(cosφ2)+oλ/ds成立。因此，實(shí)驗(yàn)1和實(shí)驗(yàn)2證明了所提出的算法已解決了在擴(kuò)展孔徑的二維DOA估計(jì)算法中的奇異點(diǎn)問(wèn)題，且角度估計(jì)性能較好。

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

圖5 RMSE隨陣元間距的變化Fig.5 RMSE versus element spacing

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

本實(shí)驗(yàn)研究所提出算法的估計(jì)性能隨信噪比的變化情況。其中，N=200，L=50，ds=8λ。假設(shè)有K=4個(gè)等功率非相關(guān)信號(hào)入射到天線陣列，信號(hào)的方位角和俯仰角分別為(60°,90°)，(60°,110°)，(80°,90°)和(80°,110°)。算法的角度估計(jì)性能隨信噪比的變化情況如圖6所示。

圖6 RMSE隨SNR的變化Fig.6 RMSE versus SNR

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

本實(shí)驗(yàn)研究所提出算法的估計(jì)性能隨快拍數(shù)的變化情況。其中SNR=0 dB，其他實(shí)驗(yàn)條件均與仿真實(shí)驗(yàn)4相同。算法的角度估計(jì)性能隨快拍數(shù)的變化情況如圖7所示。

圖7 RMSE隨快拍數(shù)的變化Fig.7 RMSE versus the number of snapshots

從圖6和圖7可以看出，隨著SNR和快拍數(shù)的增加，兩種算法的RMSE均減小,證明了所提出算法的有效性。值得注意的是，所提出的算法比文獻(xiàn)[11]中的算法具有更好的角度估計(jì)性能。這是由于所提出的算法使用了擴(kuò)展孔徑的非均勻L型陣列，因此陣列孔徑更大。

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

本實(shí)驗(yàn)研究所提出算法對(duì)應(yīng)不同角度的估計(jì)情況。其中快拍數(shù)N、數(shù)據(jù)幀數(shù)L、信噪比SNR和擴(kuò)展孔徑ds分別為200，10，15 dB，8λ。信號(hào)的方位角和俯仰角均在10°～80°之間以2°的步長(zhǎng)變化。圖8為對(duì)應(yīng)不同角度的角度估計(jì)聯(lián)合均方誤差。

圖8 不同方位角俯仰角的RMSEFig.8 RMSE for different azimuth-elevation angles

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

由于文獻(xiàn)[18]中的算法不能實(shí)現(xiàn)欠定條件下的角度估計(jì)，故研究其算法實(shí)現(xiàn)對(duì)兩個(gè)信號(hào)源(60°,70°)，(60°,80°)的角度估計(jì)。其中，信噪比SNR=15 dB。其他實(shí)驗(yàn)條件均與仿真實(shí)驗(yàn)4相同。圖9為兩種算法角度估計(jì)值分布對(duì)比散點(diǎn)圖。由圖9可知，與文獻(xiàn)[18]中的算法相比，本文所提出的算法在解決奇異點(diǎn)問(wèn)題時(shí)的估計(jì)性能更好。

圖9 情況1下的角度估計(jì)性能Fig.9 Angle estimation performance for Case 1

5 結(jié)論

本文提出了一種擴(kuò)展孔徑的二維DOA估計(jì)算法?；诜蔷鶆騆型陣列，所提出的算法首先構(gòu)造了四個(gè)延時(shí)互相關(guān)矩陣并對(duì)延時(shí)互相關(guān)矩陣進(jìn)行奇異值分解得到了信號(hào)子空間。其次，利用選擇矩陣從信號(hào)子空間中獲得對(duì)角矩陣。再次通過(guò)聯(lián)合對(duì)角化方法得到自動(dòng)配對(duì)的精確和模糊的方向余弦估計(jì)值。最后，以低精度無(wú)模糊的方向余弦估計(jì)值為參考，從高精度模糊的方向余弦估計(jì)值中獲取高精度無(wú)模糊的方向余弦估計(jì)值。所提出的算法能夠?qū)崿F(xiàn)方位角和俯仰角的自動(dòng)配對(duì)。仿真結(jié)果表明所提出的算法有效地解決了擴(kuò)展孔徑的二維DOA估計(jì)算法中的奇異點(diǎn)問(wèn)題。