項 楊 楊晉生

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

1 引言

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

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

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

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

2 陣列結構和數據模型

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

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

(1)

(2)

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

(3)

3 提出的算法

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

3.1 延時互相關矩陣的構造

(4)

因此，根據KR運算得到的延時互相關矩陣表示如下

(5)

(6)

(7)

(8)

根據式(8), 對應的陣列形式如圖2所示。

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

3.2 方向余弦的估計

通過對R進行奇異值分解(Singular Value Decomposition，SVD)，可以得到信號子空間Us以及具有K個較大奇異值的對角矩陣Λs

(9)

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

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

(10)

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

(11)

四個平面陣列中均包含x軸高精度模糊的方向余弦信息，對應對角矩陣表示如下

(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估計的實現(xiàn)

(16)

(17)

根據以上的分析，第k個信號的方位角和俯仰角估計表達式如下

(18)

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

(1)構造式(7)所示的延時互相關矩陣。

(4)同理得到x軸上對應的兩種方向余弦估計值。

3.4 算法分析

(1)角度估計

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

情況 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)生間等生活污水經化糞池預處理后，與經隔油池處理后的食堂含油廢水通過生活污水排水管道系統(tǒng)匯集，經一體化生活污水處理設施處理，達到標準后回用至生產回水水池。

其中，o為正整數表示兩個來波信號在同一方向上的方向余弦估計值的差值為整數倍λ/ds。

證明如下：

由于R=ARss，故通過判斷A的階數即可得到R的階數。

1)首先，判斷A中互不相關行的個數。

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

當角度關系滿足情況 1時，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/λ，結合3.2中關于高精度模糊的方向余弦信息的選擇矩陣的維數可知A中線性無關的行數至少為2M2。

當角度關系滿足情況 2時，與情況1同理。

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

2)其次，判斷A中互不相關列的個數。

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

綜上所述，由于2M2>K，因此通過對R進行SVD分解可以得到信號子空間Us以及K個較大奇異值。

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

(2)克拉美羅界

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

(19)

4 仿真實驗

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

(20)

在所有的仿真實驗中，所提出算法中子陣的陣元個數M=3，文獻[11]中子陣的陣元個數為6，即L型陣列中共有11個陣元。由于文獻[18]中算法應用的陣列結構不同，令其共包含13個陣元。此外，所有的仿真實驗均進行1000次蒙特卡洛仿真。

仿真實驗1

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

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

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

仿真實驗2

當方位角和俯仰角滿足情況3時，驗證所提出算法的有效性。假設有K=2個信號入射到天線陣列，兩個信號分別為(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。其他實驗條件均與仿真實驗1相同。圖4為角度估計值分布散點圖。值得注意的是文獻[18]中的算法不能實現(xiàn)此種情況下的角度估計。

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

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

仿真實驗3

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

仿真實驗4

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

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

仿真實驗5

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

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

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

仿真實驗6

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

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

仿真實驗7

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

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

5 結論

本文提出了一種擴展孔徑的二維DOA估計算法。基于非均勻L型陣列，所提出的算法首先構造了四個延時互相關矩陣并對延時互相關矩陣進行奇異值分解得到了信號子空間。其次，利用選擇矩陣從信號子空間中獲得對角矩陣。再次通過聯(lián)合對角化方法得到自動配對的精確和模糊的方向余弦估計值。最后，以低精度無模糊的方向余弦估計值為參考，從高精度模糊的方向余弦估計值中獲取高精度無模糊的方向余弦估計值。所提出的算法能夠實現(xiàn)方位角和俯仰角的自動配對。仿真結果表明所提出的算法有效地解決了擴展孔徑的二維DOA估計算法中的奇異點問題。