非圓信號二維DOA和初始相位聯合估計方法

王 凌 李國林 謝 鑫 齊 率

(海軍航空工程學院 煙臺 264001)

1 引言

近年來利用目標空間中信源的信息冗余特性來提高空間譜估計的性能得到了廣泛關注。信號的非圓性作為時域上的信息擴充，Galy首先將其引入空間譜估計研究中，提出了非圓信號測向的MUSIC算法[1](MUSIC algorithm for NonCircular signals,NC-MUSIC),Abeida和Delmas進一步對非圓信號MUSIC類測向算法的Cramer-Rao Bound (CRB)進行了分析[2]。近年來國內外不少文獻對相關問題進行了研究[3－12]。目前，研究較多的一類非圓信號為直線信號，例如幅度鍵控信號和二進制相移鍵控信號等，其相位星座圖沿直線分布，非圓率ρ=1。文獻[3]中利用非圓信號為實值信號的特點，通過每次陣列快拍構造了一組數據矩陣，實現了非圓相干信號的完全解相干。文獻[4,5]在構造的 Toeplitz矩陣中均利用非圓信號的共軛信息，降低了陣列孔徑損失，提高了陣列自由度。文獻[6,7]利用非圓信號橢圓協方差E[ss]≠0特性，構造虛擬陣元，提高了角度估計的精度。但上述算法模型背景中均假設非圓信號的初始相位為零，這在實際環境中幾乎不可能滿足。容易證明，非圓信號s(t)和其任意旋轉s(t)ejφ具有不同的概率分布，因此對于非圓信號來說，初始相位不能忽略。如果考慮初始相位，上述基于非圓信號特性的算法均會存在估計性能下降的情況，特別是文獻[3－5]中直接利用非圓信號的實值特性構造數據矩陣，當考慮到初始相位時，信號模型將不再滿足實值特性，對算法的影響是嚴重的。近幾年，部分文獻對上述算法進行了修正[8,9]。

以上文獻都是針對1維DOA估計的非圓算法，針對非圓信號背景下的2維DOA估計問題，公開文獻中研究很少。僅有Haardt[13,14]等提出將2維酉ESPRIT用于非圓信號測向，以及劉劍[15,16]提出的基于雙平行線結構的擴展 MUSIC算法。但是上述算法運算量都較大，需要構造維數為33MM×的數據矩陣，并進行特征分解。2維DOA和初相的聯合估計可以看作是在3維參數空間進行點估計，基本思路是通過3維譜峰搜索實現，但運算量巨大。或者根據二次型的性質，將多維搜索轉化為1維搜索或求根。

針對上述方法計算量大，且算法復雜的問題，本文提出了一種基于垂直陣列結構的任意初始相位非圓信號2維DOA和初相聯合估計方法，利用陣列結構特點將3維參數估計轉化為3個子陣上的2維參數估計問題，然后在每一個子陣列上，同時利用信號子空間的旋轉不變性和噪聲子空間的正交特性將2維估計轉化為1維估計，且得到的2維DOA和初始相位能實現自動配對。該算法只需對 3個22MM×維的數據矩陣進行一次特征分解，且可在3個維度并行處理，相較于文獻[15]中算法，運算量得到降低。由于利用了非圓信號特性，本文算法相較于具有相同陣列結構的文獻[17]中算法，估計精度和可估計信源數都得到提升，且可估計多于陣元數的信源，最多可估計2(M?1)個信源。

2 陣列信號模型

考慮采用如圖1所示的傳感器陣列系統模型，該模型由3個相互垂直的均布線性子陣組成，記為，各子陣列陣元均為全向的，且陣列位于遠場。第1個子陣由陣元0,1,…,M?1組成，第2個子陣由陣元0,1,…,N?1組成，第3個子陣由陣元0,1,…,P?1組成，各個子陣的陣元間距分別為dx,dy,dz。空間有D個相同載頻的不相關窄帶信源入射該陣列系統，信源數D假設為已知或已估計得到，波達方向矢量角分別為 (θ1,θ2,…,θD)，其中θk=(αk,βk,γk),αk,βk,γk分別為第k個源信號入射方向與x軸，y軸和z軸的夾角。顯然有如下關系成立

可知3個角度變量中只有兩個獨立，因此，波達方向矢量角可以用1×2的矢量表示，即θk=(αk,βk)。稱αk,βk分別為方位角和俯仰角。以陣元0為參考陣元，則子陣列Xa中第i個陣元的輸出信號可表示為

其中sk(t)為第k個直線信號，,λ為其中心波長，nxi(t)為子陣列Xa中第i個陣元中的復圓高斯白噪聲，噪聲功率為σ2，且噪聲與信源不相關。直線信號，其中φk為第k個直線信號入射到參考陣元的初始相位。文獻[3－7]中均假定φk為零，然后利用實值性進行處理。本文考慮較為合理的情形，即各直線信號的非圓相位不全為零。

則子陣列Xa的輸出信號矢量可表示為

式中

為M×D維陣列流型矩陣，a(αk)為對應的方向向量，且有

同理可以得到子陣列Ya和Za的輸出信號矢量

式中

3 本文算法描述

本文算法主線是利用空間信號入射到垂直陣列后的陣列流型特點，將2維DOA和初相聯合估計問題轉化為3個可并行處理的2維參數估計，然后在每一子陣列上，同時利用信號子空間的旋轉不變性和噪聲子空間的正交性，通過一次特征分解得到的信號和噪聲子空間將1維DOA和初始相位的2維搜索轉換為1維估計問題，這樣就把3維參數估計問題全部轉化為1維估計，使算法運算量得到大幅的降低。

3.1 算法描述

受文獻[13,18]啟發，并利用非圓信號橢圓協方差非零特性，首先在x軸向構造

式中，ΠΜ為M×M維置換矩陣，其反對角線上元素為1，其余元素為0。,從bxk(αk,φk)的表達式中可以看出，此時已將(α,β,φ)的3維估計轉化為(α,φ)的2維估計問題。改寫bxk如下式

式中⊙為哈達瑪積。從bxk的表達式不難發現，通過利用信號非圓特性構造式(14)等效于在x軸向的負向端擴展了一倍的陣列孔徑，且擴展虛擬陣元的幅相增益為e?2jφk。按照上文思想，在y軸和z軸向同樣構造擴展的觀測矢量Gy(t)和Gz(t)，則經過擴展后的虛擬陣列結構如圖2所示。

圖2 擴展后的虛擬陣列系統模型

先計算x軸向的擴展協方差矩陣

式中UxS為D個大特征值對應的特征矢量張成的信號子空間，UxN為2M?D個小特征值對應特征矢量張成的噪聲子空間。如果直接利用噪聲子空間的正交性進行估計，即x軸向的方向向量bxk與噪聲子空間UxN正交，得到

由于噪聲影響，式(18)并不完全接近于0，此時初始相位和1維DOA的估計就轉化為最小化下式來實現

此時需要對α,φ進行2維搜索，為了避免復雜搜索，本文考慮在利用噪聲子空間特性的同時使用信號子空間的旋轉不變性來降維。以x軸向為例，擴展后的虛擬陣列如圖3所示。經典的子陣劃分方法為圖中虛線所示。但由于非圓信號初始相位的影響，導致按圖中虛線劃分的子陣間將不再滿足旋轉不變性。由于虛擬陣元幅相增益的存在，可以考慮按實線方式進行子陣劃分，此時

圖3x 軸向子陣劃分

其中V1和V2為選擇矩陣，同樣令

因為Bx2=Bx1Ψx，其中Ψx=diag[u(α1),u(α2),…,u(αD)]，根據信號子空間旋轉不變原理，易得：，其中Tx為一個非奇異轉換矩陣，可得

利用最小二乘或總體最小二乘方法求解式(24)得到兩子陣信號子空間之間的旋轉不變關系矩陣Θx，通過對Θx特征分解得到的Ψx即可求得空間信源的方位角αk，然后將得到的αk代入式(19)，通過1維搜索即可求得對應于第k個信源的初始相位φk，這樣就把(α,φ)的2維估計問題轉化為兩個1維估計。

同理，對另外兩個軸向的旋轉不變關系Θy和Θz進行特征分解可得到Ψy和Ψz，進而可計算出各個信號與y軸和z軸之間的夾角βi和γj。2維角度參數的配對可以通過計算d(k,i,j)的最小值來實現，d(k,i,j)的定義式為

3.2 算法步驟

由上文分析，可得到本文算法步驟如下：

(1)通過陣列天線獲取垂直陣列的觀測矢量X(t),Y(t)和Z(t)；

(2)根據式(14)形式構造擴展觀測矢量Gx(t),Gy(t)和Gz(t)；

(3)計算擴展協方差矩陣Rxx,Ryy和Rzz；

(4)對擴展協方差矩陣進行特征值分解，得到各維度上相應的信號和噪聲子空間；

(5)利用式(24)中信號子空間的旋轉不變性計算信源與各軸之間夾角α,β和γ，并用式(25)進行配對；

(6)利用上步得到結果代入式(19)中，根據噪聲子空間的正交性得到對應的初始相位φ。

4 數值仿真

為驗證本文算法的有效性，從3個方面進行數值仿真驗證，仿真中對比了具有相同陣列結構的文獻[17]中算法。第1步驗證本文算法對空間信源2維DOA和初始相位的聯合估計性能；第2步驗證本文算法對快拍數和信噪比的敏感度。最后驗證本文算法的抗過載能力，即多于陣元數的信源估計能力。

4.1 2維DOA和初始相位聯合估計驗證

仿真1利用圖1中的陣列系統模型，取陣元數M,N和P均為8，為避免估值模糊，陣元間距dx,dy和dz均取為λ/2。3個空間信號的2維到達方向分別為(30°,90°,60°),(60°,60°,45°),(85°,75°,1 5 .9°)。其附加初始相位分別為ejπ/6,ej2π/3和ejπ/3。信噪比為15 dB，快拍數為500，進行100次Monte-Carlo試驗。圖4(a)和4(b)就是采用本文算法得到2維DOA估計的星座圖和初始相位的 MUSIC譜圖。從仿真結果可知，本文算法能夠準確地估計得到2維DOA和初始相位，并能實現自動配對。

圖4 2維DOA和初相聯合估計

4.2 信噪比和快拍數對本文算法的影響

仿真1只是定性說明了本文算法的估計性能，下面通過兩個仿真試驗驗證本文算法的估計均方根誤差大小隨信噪比和快拍數的變化趨勢，并加入文獻[17]中算法作為比較。

仿真2陣列和信號模型與仿真1相同，文獻算法采用同樣模型參數。快拍數固定為1000，改變信噪比，在0 dB到20 dB分別進行50次的Monte-Carlo統計試驗，得到如圖5所示的2維DOA估計性能與信噪比變化曲線。從圖中可以看出，本文算法對每個信號的估計均方根誤差都優于文獻[17]中算法，并且隨著信噪比的提升，算法估計誤差逐漸減小。定義第k個信號的均方根誤差為

仿真3仿真參數同上。固定信噪比為5 dB，改變快拍數，在400次到1400次分別進行50次的Monte-Carlo統計試驗，得到如圖6所示的估計誤差與快拍數變化曲線。從仿真結果可得，隨著快拍數的增加，均方根誤差逐漸減小，在快拍數高于900次后，誤差趨于收斂。總體來看，本文算法在低信噪比和短快拍環境下，對2維DOA估計的均方根誤差能控制在1°以內，且優于文獻中算法。

4.3 陣列抗過載性能檢驗

由于使用了利用非圓信號特性構造的擴展觀測矢量，因此本文能夠估計的空間信源數為2(M? 1 )，能夠對多于陣元數的空間信號進行估計，此時文獻[17]中算法已經失效。

圖5 估計誤差隨信噪比變化曲線

仿真4取陣元數M,N和P均為4，4個空間信號的2維到達方向分別為(30°,90°,60°)，(60°,60°,4 5 °),(85°,7 5 °,1 5 .9°),(70°,2 5°,7 5 .6°)。其附加初始相位分別為ejπ/6,ej2π/3和ejπ/3,和ej5π/6。信噪比為20 dB，快拍數為500，進行100次Monte-Carlo試驗。圖7(a)和7(b)就是采用本文算法得到2維DOA估計的星座圖和初始相位的MUSIC譜圖。此時陣列處于過載狀態，傳統的子空間估計方法失效，本文算法仍能準確得到結果。圖7(a)和7(b)中不同來波方向信源的2維DOA和初始相位估計誤差有區別的原因是仿真中采用了不同的信號形式。

5 結論

以往的2維DOA和初始相位聯合估計問題，一方面需要復雜的3維搜索，另一方面需要對高維度數據矩陣進行特征分解運算，針對傳統算法的復雜性，本文提出了一種基于垂直陣列結構的任意初始相位非圓信號2維DOA和初相聯合估計方法，利用垂直陣列特點，將3維參數估計問題轉化為可并行處理的3個2維參數估計，在每一個子陣上，同時使用噪聲子空間正交性和信號子空間旋轉不變性，將2維參數估計進一步轉化為1維估計問題，最終只需要對協方差矩陣進行一次特征分解即可實現 2維DOA和初相的聯合估計及自動配對。由于使用擴展觀測矢量構造的協方差矩陣，本文算法適用于空間信源處于過載的情形和低信噪比、短快拍環境，因此適合多目標和干擾下對2維DOA和初相估計要求較高的機載、彈載陣列天線系統和小型地面陣列天線系統。數值仿真驗證了本文算法的上述優良性能。

圖6 估計誤差隨快拍數變化曲線

圖7 過載環境下聯合估計