劉淳熙,王文亮,彭冬亮*,陳志坤,吳美嬋
(1.杭州電子科技大學 自動化學院,浙江 杭州 310018;2.中船(浙江)海洋科技有限公司,浙江 舟山 316021)
隨著電磁干擾和抗干擾技術的飛速發展,電磁環境日益復雜,電磁信號所蘊含的信息也越來越多[1-3]。普通陣列僅具有空域濾波能力,不能很好地區分期望信號和干擾信號,極化敏感陣列不僅能利用陣元間的相對空間位置進行空域濾波,還可以利用期望信號和干擾信號極化狀態的差異在極化域濾波[4-6],因此極化敏感陣列具有較強的抗干擾能力。
20世紀90年代以來,國內外學者對于極化敏感陣列的研究日益活躍,許多傳統標量陣列的經典參數估計算法被推廣到極化敏感陣列領域。
文獻[7-8]利用旋轉不變子空間(Estimation of Signal Parameters via Rotational Invariance Techniques,ESPRIT)算法實現了聯合參數估計,但該方法需要對所估計的空間角度和極化參數進行額外配對,增加了計算量。文獻[9]提出一種改進的ESPRIT算法,利用空域旋轉不變關系得到6個旋轉不變因子,通過數學運算得到一系列估計值,然后重構出多組對應的導向矢量,將其與噪聲子空間一一正交檢驗,提取出正確的一組無模糊估計值,但算法運算復雜度較高,且ESPRIT類算法均對陣列結構有嚴格要求,不利于實際應用推廣。文獻[10]將MUSIC算法直接拓展到極化敏感陣列,通過多維搜索即可獲得空域和極化域的參數信息。但是由于MUSIC算法需要進行譜峰搜索,同時對空域和極化域參數進行估計會使計算量急劇增大,因此有必要研究相關的降維算法。文獻[11]提出了一種模值約束降維求根MUSIC算法,通過數學運算對波達角和極化參數進行剝離,先估計DOA信息,然后根據模值約束條件構造代價函數,通過閉式解求出極化參數。文獻[12]提出一種恒模約束降維算法,對陣列接收矩陣用秩-1逼近恢復出信號的空域導向矢量和極化域導向矢量,然后估計出對應的參數。這些算法均能較為精確地估計出DOA和極化參數,但依然受到陣型的限制,僅適用于線陣或L陣。隨著四元數模型在空間定位等領域的研究愈發廣泛,一些學者將四元數模型與陣列參數估計算法相結合,提出了四元數MUSIC算法[13]、四元數求根MUSIC算法[14]以及降維四元數MUSIC算法[15],將其應用到極化敏感陣列中。該方法雖然可以適當地降低一部分運算量,但無法精確估計極化參數。
實際工程中普遍應用二維面陣,現有的針對線陣或L陣的參數估計算法實際應用價值有限。因此本文提出了一種適用于極化敏感面陣的不等式約束降維MUSIC算法。在強噪聲的多信源背景下,本算法能夠較好地識別出信源的空域-極化域參數。算法首先利用Kronecker積的數學性質將極化域矢量與空域矢量分離,然后將四維聯合譜估計問題轉化為二元函數極值的求解,之后利用極化矢量的模值有界性構建不等式約束方程,得到空域參數估計值,同時將空間角度與約束方程結合得到極化參數估計。本算法大大降低了運算量,同時能夠精確估計出DOA和極化參數。
本文研究極化敏感面陣,其陣元由雙正交偶極子對構成,圖1是由該陣元構成的M行N列的均勻面陣,圖中的“黑色圓點”代表雙正交偶極子。陣列擺放在XOY平面上,陣元間距d為半波長,即λ/2。定義俯仰角θ為與z軸正向的夾角,0≤θ≤π;定義方位角φ為順時針與x軸正向的夾角,0≤φ<2π。

圖1 極化敏感面陣結構Fig.1 Structure diagram of polarization sensitive array
假設空間中K個遠場窄帶信號入射到該面陣,其空間角度用(θk,φk),k=1,2,…,K表示,其中θk和φk分別表示第k個信號源的俯仰角和方位角。x軸和y軸上信號源的方向矢量分別為:
(1)
(2)
則陣列的空域導向矢量為[16]:
as(θk,φk)=ay(θk,φk)?ax(θk,φk),
(3)
式中,符號?表示Kronecker積。
電磁信號在空間中的傳播如圖2(a)所示,將信號的電場矢量E在兩正交單位矢量eθ,eφ方向上分解,如圖2(b)所示。

(a) 空間電磁信號傳播

(b) 電場矢量分解圖2 電磁信號傳播與分解Fig.2 Electromagnetic signal propagation and decomposetion diagram
圖2(b)中Vθ和Vφ為電場矢量E的2個分量,其中:
Vθ=sinγejη,Vφ=cosγ,
(4)
式中,γ∈[0,π/2]和η∈[-π,π]分別表示極化輔角和極化相位差,γ的正切值表示了eθ方向電場幅度與eφ方向電場幅度的比,η表示eθ方向電場與eφ方向電場的相位差。
在空間直角坐標系中,電場矢量E可以分解為:
E=Vθeθ+Vφeφ=
(Vθcosθcosφ-Vφsinφ)ex+
(Vθcosθsinφ+Vφcosφ)ey-(Vθsinθ)ez。
(5)
將矢量寫成矩陣形式:
(6)
式中,γk,ηk分別表示第k個信源的極化輔角和極化相位差。本文所提陣列陣元為雙正交偶極子對,僅能接收到與x軸和y軸平行的電磁波,則陣列的極化域導向矢量為:
(7)
因此陣列的極化域-空域聯合導向矢量為:
a(θk,φk,γk,ηk)=ap(θ,φ,γ,η)?as(θk,φk)。
(8)
當K個遠場窄帶完全極化波入射到該極化敏感面陣時,陣列的接收數據x(t)可以寫成下述形式:
AS(t)+N(t),
(9)
式中,S(t)∈CK×1為入射信號矢量矩陣;N(t)∈C2MN×1為空-時-極化域白噪聲矩陣;A∈C2MN×K為陣列流行矩陣:
A=[a1a2…ak…aK],
(10)
式中,
ak=a(θk,φk,γk,ηk)。
(11)
根據式(9)所構建的信號模型,對K個目標的二維空間角度和二維極化角度進行估計的過程稱為極化域-空域聯合譜估計。

(12)
式中,X=AS+N為P次快拍的接收數據矩陣。
對其進行特征值分解得到:
(13)
式中,US表示K個較大特征值所對應的特征向量構成的信號子空間;UN表示2MN-K個較小特征值所對應的特征向量構成的噪聲子空間。
通過分析可知,信號子空間與陣列的導向矢量張成的空間是同一空間,信號子空間與噪聲子空間相互正交[17],則陣列的導向矢量與噪聲子空間也正交,即:
aH(θk,φk,γk,ηk)UN=0,k=1,2,…,K。
(14)
考慮到實際接收數據矩陣是有限長的,且實際接收數據矩陣中混有噪聲,因此a(θ,φ,γ,η)與UN不完全正交,即式(14)不完全成立,實際上求DOA是以最小優化搜索實現的,即四維聯合譜函數定義為:
(15)
遍歷θ,φ,γ,η四個參數,找出譜函數最大的K個峰值對應的方位角、俯仰角、極化輔角和極化相位差即為K個來波信號的參數估計值。
該傳統MUSIC算法需要進行四維全局譜峰搜索,運算量極大,難以實現。本文利用極化參量模值的有界性,引入不等式約束將四維譜峰搜索降低至二維,先行估計二維空域參數,進而得到極化參數。
結合式(7)和式(8),把陣列導向矢量寫成僅包含空域參數的矢量和僅包含極化域參數的矢量的乘積形式:
Γ(θk,φk)σ(γk,ηk),
(16)
式中,
(17)
(18)
則式(15)譜函數可以改寫成:
PMUSIC(θ,φ,γ,η)=
(19)
分析式(18)可知極化參數矢量σ(γ,η)的模值有如下關系:
(20)
令:
(21)
則四維譜函數的譜峰搜索過程可以轉化為在不等式約束條件下的極值求解問題:
(22)
從四維譜峰搜索轉換到不等式約束條件下的極值求解過程的證明過程可參見文獻[18],其中定義檢測量為:
(23)
令:
(24)
以am(β)為待求解變量,對式(23)進行角度估計可轉化為求解以下優化方程:
(25)
式中,e1=[1,0,…,0]T。本文中用σ(γk,ηk)和Γ(θk,φk)分別替換式(23)中的am(β)和[an(α)?IM],以σ(γk,ηk)為待求解變量,相應改變對其模值的約束條件即可將式(19)的四維譜峰搜索問題轉化為式(22)所示的優化方程求解問題。
構造代價函數:
L(θ,φ,γ,η,λ,ω)=σ(γk,ηk)HV(θ,φ)σ(γk,ηk)-
λ(‖σ(γk,ηk)‖2-3+ω2),
(26)
式中,λ是拉格朗日乘子;ω是新引入的松弛變量,目的是為了將不等式約束經過松弛后變為等式約束。對式(26)求導,并令結果為零:
(27)
進一步化簡得:
V(θ,φ)σ(γk,ηk)=λσ(γk,ηk)。
(28)

(29)
對應求解即可得到第k個信源的極化參數估計值。
本節對4D-MUSIC算法和本文所提的RD-MUSIC算法的計算量進行比較分析,2種算法均選用圖1所示的陣列結構,豎直和水平陣列參數分別為M和N,入射信號數為K,快拍數為L,搜索的角度點數為n。
對于4D-MUSIC算法,求陣列接收信號自相關矩陣計算量為O{L(2MN)2},矩陣特征值分解計算量為O{(4MN)3},四維角度搜索計算量為O{n4(4MN+1)(4MN-K)},總計算量為O{L(2MN)2+(4MN)3+n4(4MN+1)(4MN-K)}。
本文所提降維MUSIC算法計算自相關矩陣和獲得特征值分解的計算量與4D-MUSIC算法相同,角度搜索時先對二維空間角度進行搜索,后對二維極化角度搜索,計算量為o{n2(16MN+4)(4MN-K)} ,總計算量為o{L(2MN)2+(4MN)3+n2(16MN+4)×(4MN-K)}。
圖3給出了信源數L=2,分別采用1°粗搜索和0.1°精搜索時,4D-MUSIC算法和本文算法計算量隨著陣元數變化的對比關系。

圖3 2種算法復雜度對比Fig.3 Complexity comparison of two algorithms
為了更直觀地看出本文算法在運行時間上的優勢,表1給出了4D-MUSIC算法和本文算法的運行時間隨著陣元數變化的對比關系。

表1 4D-MUSIC算法和本文算法的運行時間隨著 陣元數變化的對比關系
由圖3和表1可見,本文提出算法的復雜度實際上遠低于MUSIC算法的復雜度。這是因為在算法運行過程中M,N,K和L遠遠小于n,傳統MUSIC算法的運算量為n4量級,而本算法運算量為n2量級,運算量大大降低。
本節用計算機仿真來驗證本文所提算法的準確性及有效性。仿真所用陣列結構如圖1所示,陣元數M=4,N=2,陣元間距為λ/2。
仿真1目標DOA與極化參數估計
考慮2個非相干信號入射到陣列,信號俯仰角為θ=[5°,20°];方位角φ=[10°,50°];極化輔角γ=[30°,45°];極化相位差η=[60°,80°]。快拍數snap=1 000,信噪比SNR=20 dB。圖4給出了空域參數的MUSIC譜圖像,圖上顯示出2個清晰可辨的峰值,說明本算法能夠正確估計出目標的二維空間參數。圖5和圖6分別給出了本算法所估計出的二維空間角度和二維極化角度的坐標圖。從圖6看出,本算法也能夠對極化參數進行較為精準的估計。

圖4 空域MUSIC譜Fig.4 Spatial MUSIC spectrum

圖5 方位角與俯仰角坐標Fig.5 Coordinates of azimuth and elevation angle

圖6 極化輔角和極化相位差坐標Fig.6 Coordinates of polarization auxiliary angle and polarization phase difference
仿真2算法性能與信噪比關系
本仿真對比了傳統的4D-MUSIC算法和本文算法的空域、極化域參數估計性能。蒙特卡羅次數L為500,信噪比從-10 dB變化至20 dB,用均方根誤差(Root Mean Square Error,RMSE)評估算法的空域和極化域參數估計性能,實驗中RMSE計算如下:
(30)


圖7 DOA參數的RMSE隨SNR變化曲線Fig.7 RMSE curve of DOA parameters varied with SNR

圖8 極化參數的RMSE隨SNR變化曲線Fig.8 RMSE curve of polarization parameters variedwith SNR
本文首先建立了一種適用于均勻平面陣列的長矢量接收信號模型,然后提出了一種適用于該陣列的空域-極化域聯合估計算法。該算法利用極化矢量范數的有界性,將四維聯合譜搜索過程轉化為二維的約束優化問題,在估計精度基本不變的情況下顯著降低了計算量。仿真實驗證明算法在較低信噪比條件下仍能得到較高精確度的估計值。