基于最大期望算法的多時變信號DOA估計方法

范保華, 左 樂,*, 唐 勇,2, 胡澤華

(1. 中國電子科技集團公司第二十九研究所, 四川 成都 610036; 2. 電子科技大學電子科學與工程學院, 四川 成都 611731)

0 引 言

輻射源參數估計問題在雷達、聲吶、麥克風陣列、射電天文以及移動通信中有廣泛應用,因此其在陣列信號處理領域具有重要的研究價值[1-5]。微波測向通過在不同空間位置布置天線接收來波的電磁信號進行。根據布置天線的方式可分為靜態陣列和時變合成陣列。靜態陣列每個單元均配有對應的接收通道,通過不同陣元同時采集空間中來波信號的場分布,并通過參數估計等方法完成測向[5]。靜態陣的優點是可以完成單脈沖測向,但缺點在于需通過較多的天線和相應的接收通道來解模糊及高精度測向。

與同時采集來波信號場分布的靜態陣不同,時變陣列在測量過程中移動采樣位置,通過時間積累和空間掃描,利用簡化的采樣硬件獲取來波的空間場分布信息。將同時采集用分時采集替代,在處理時間與硬件資源之間進行權衡。

無源合成陣列(passive synthetic array, PSA)既是典型的時變測向陣列。其通過較少的采樣單元數目,通過分時改變陣元位置來實現一維或二維場分布采集[6-7]。1996年,Friedlander首次將時變陣列用于測向,將靜態陣基于本征值的測向方法應用于時變陣[6]。

無源合成陣列由于其簡單的系統架構和較高的測向精度被廣泛應用,但在實際工程應用中易遇到以下難點:

(1) 信號分選問題。采用靜態陣測向時,由于采用多通道同時采樣,可同時測量每個信號的幅度、頻率、相位等信息,并以此計算得到該信號的入射角。然而采用時變陣測向時,每一次采樣僅能獲得入射場的幅度、相位、頻率等數據,無法建立多個數據間的顯性關聯。當多個輻射源在同一時間段內交替輻射信號時,采樣數據相互交錯,但又無法事先獲得單個采樣數據與其輻射源的對應關系。反觀入射角的估計,需首先完成信號分選,建立每個采樣信號與輻射源之間的關聯關系,才能采用來自同一個輻射源的數據進行該輻射源信號參數的精確估計。

(2) 采樣過程中信號幅度時變問題。如對于來波為波束掃描的雷達信號[7],在不同時刻收到信號幅度不同,導致其接收的信噪比不穩定,而不同信噪比(signal to noise ratio, SNR)的信號置信度不完全一樣,對于高精度測向需該考慮不同質量信號對測向結果的貢獻程度。

(3) 輻射信號頻率時變問題。對于頻率間隔大于信號帶寬的信號,可以從頻域上將其分為多個輻射源。但對于頻率間隔小于信號帶寬的信號,如跳頻信號或頻率調制信號,此時信號的頻率時變特性就必須謹慎處理,因為若將同一輻射源產生的變頻信號分為多個輻射源進行處理,則會造成信號分選錯誤。不僅如此,每組信號數據量的減少還會導致測向精度的下降。

(4) 信號的采樣位置取決于信號到達的時間,其分布并非均勻,空間采樣位置無周期性。

在眾多的波達方向(direction of arrival, DOA)估計方法中,最大似然(maximum likelihood, ML)方法被認為性能穩健且統計效率較高[5,8-9]。同時,ML方法容易結合最大期望(expectation maximization, EM)算法[10],進行信號的分選以及多信號的參數估計。EM算法可將多信號進行分離,其最大優勢在于能將多目標的優化問題分解為幾個獨立的單目標優化問題,通過反復迭代將問題降維,從而實現運算量的縮減。Weinstein和Feder首先將EM算法應用于測向領域,實現多徑時延和多輻射源定位[11]。隨后,Chung和Frenkel分別將EM算法應用于多目標點頻信號位置信息的獲取[12-13]。對于寬帶信號,Mada和Lu分別將EM算法應用于信號位置參數和波形的聯合估計[14-15]。

在時變測向陣列的移動形式方面,文獻[16-19]將無源合成線陣用于輻射源的測向和定位。雖然無源合成線陣通過陣列的直線移動形成一個大的測向物理尺寸對空間來波進行分時采集,但其存在移動路徑的左右無法分辨、移動姿態不易控制等問題。另一方面,由于陣列的旋轉容易控制,無源合成圓陣可通過兩傳感器的旋轉獲取各角度下的輻射源場分布信息,用于解算輻射源的二維入射角[20-32],并實現360°的全方位角覆蓋。文獻[21]提出采用兩天線旋轉形成一個無源合成圓陣,用于地球同步軌道衛星測向。但該方法未考慮多輻射源以及相位模糊問題,因此僅適用于單輻射源以及窄帶、小角度的測向系統。文獻[22]基于信號子空間測向方法,但要求頻率和采樣間隔固定。文獻[23]采用粒子群優化算法,進行目標參數計算,但是由于粒子群算法采用隨機搜索策略,其全局尋優能力并不十分可靠[24]。文獻[25]所采用的陣元數較大,且只適用于單信號。文獻[26-27]采用多假設偽線性迭代最小二乘方法,但未考慮信號分選問題,故僅能處理單輻射源情形。對于多信號問題,文獻[28]僅利用相位差信息,提出采用正交譜配對追蹤方法解決了多目標問題。但由于其譜值的計算過程要求頻率固定、采樣位置均勻,因此仍無法適應數據采集過程中信號頻率、幅度、位置間隔改變等問題。除此之外,文獻[29-31]雖然實現了多目標二維角度的無模糊估計,但均未考慮各采樣點信噪比各異的情況,會造成測向精度下降。

為同時解決上述4個無源合成陣列測向及分選所面臨的問題,本文提出采用EM算法,同時解決頻變多輻射源的測向以及信號分選問題。據筆者所知,尚未見將EM算法應用于無源合成陣列多目標測向的報道。對于每個輻射源的DOA估計,本文采用ML方法進行二維測向,并考慮每個接收信號采樣位置間隔、頻率、SNR均不同的情形。但直接采用復數形式的接收信號對來波方位角及仰角的ML估計是一個非線性、多維的極值問題,需要全局多維搜索,其計算量較大。由于不同DOA產生的信號相位不同,相位干涉儀同樣適用于時變測向陣列。注意到接收信號的相位差與二維入射角的投影呈線性關系,本文采用相位數據推導了不同SNR條件下二維入射角的閉合形式解析表達式。其次,無源合成陣列與靜態陣列相同,當涉及到采用干涉儀體制測量輻射源DOA時,會出現相位翻轉的模糊問題。為了解相位模糊的多值問題,通常需要增加傳感器數量來排除多值,這無疑會增加系統和算法的復雜度。通常的相位解模糊方法根據先計算出的一個粗角度來還原相位的真實值[32]。但已報道的相位解模糊方法只適合靜態線陣[32-34],不能直接用于無源合成圓陣。為了解決相位模糊問題,本文先采用無模糊的復數形式的ML方法,進行粗測向,以減小譜峰搜索的運算量,再采用解模糊后的相位數據進行入射角的解析運算。最后,數值實驗結果表明,本文提出的基于EM算法的能完成高精度測向,其測向誤差接近克拉米勞下限(Cramer Rao lower bound, CRLB),且測向與信號分選同時完成。

1 EM算法與模型建立

1.1 EM算法

EM算法最早于1977年由Dempster等人提出[10],其思路是將用于參數估計的數據進行分組,然后在每組內進行參數的分別估計,再根據新的估計參數重新進行數據的分組,以提高參數估計的質量,并以此交替迭代,進行數據的分解與參數的估計。EM算法的流程如圖1所示。隨著參數似然值的增加,分組結果最終也將收斂。EM算法的的一個顯著特點在于多個參數的估計不通過單次處理,而是每次參數估計后,用新參數重新劃分數據,再對劃分后的數據分別進行新一輪的參數估計。EM算法采用迭代的方式將一個高維問題分解成多個低維問題,通過降維實現運算量的精簡。

結合多目標測向與信號問題,由于輻射源的位置在旋轉周期內相對固定,因此可依據輻射源的入射角參數進行觀察數據的分解,即將來自于同一輻射源的接收數據分解入一組。而接收數據與輻射源位置的關系可通過ML方法求解的二維角度參數進行表征。基于此,首先建立接收信號的模型。

1.2 接收信號模型

設間隔為d的兩天線在與XOY面內旋轉,旋轉中心為坐標原點O,如圖2所示。在一個采樣周期內,該合成圓陣接收到來自P個位置固定的輻射源發出的N個信號。輻射源輻射的電磁信號頻率、幅度均不相同。

由于采樣時間極短,采樣周期內的旋轉角度可忽略不計。第p個輻射源的第n個脈沖信號的接收信號時域表達式為

Vn(t)=Anexp(jωnt+jφp(n))+en(t),n=1,2，…,N

(1)

φp(n)=kp(n)dsinθpcos(φ(tn)-φp)

(2)

式中：旋轉角度φ(tn)∈[0,2π)為旋轉臂與x軸的夾角;仰角φp∈[0,2π)為來波方向與x軸夾角;方位角θp∈[0,π)為來波方向與z軸夾角;kp(n)=ωn/c為波數,c為波的傳播速度。本文假設信號頻率測量無誤差。

此時,通過采樣點數為M的付氏變換,將時域信號轉化為頻域信號:

(3)

可得第無模糊相位測量值為

(4)

式中:ε(n)表示相位測量噪聲。由于不同采樣點間噪聲相互獨立,等效相位噪聲的方差[35]為

(5)

此即為頻域相位噪聲與時域高斯白噪聲的關系。其揭示了相位噪聲與信噪比的關系,即在采樣數一定的條件下,相位噪聲的方差與信噪比成反比。同時,等效相位噪聲服從高斯分布[36]。

當天線間距d大于半個工作波長時,會出現相位模糊問題。由于相位的測量值在[-π,π)之內,為相位的真實值對2π取模后的余數,故為

(6)

式中:h(n)為一整數,也被稱為模糊數。由此可見,相位的測量值為信號的空間相差疊加相位噪聲后,對2π取模后的余數。

當P個輻射源在一個采樣時間內交替輻射信號,該合成圓陣的采樣信號相位為P個輻射源所產生的信號相位疊加,即

(7)

2 二維入射角估計

2.1 基于相位的二維測向方法

(8)

對條件概率密度函數的對數值取偏微分,并令其為零,即

(9)

和

(10)

有

(11)

及

(12)

根據三角函數的和差化積關系,以及每個相位差采樣相位誤差方差與SNR的反比關系,可得二維入射角的簡潔表達式為

(13)

(14)

式中:b1和b2為矩陣b=[b1,b2]T的兩個元素,可通過下式計算得到:

b=(det(A))-1[A][c1,c2]T,

(15)

式中:det(A)為矩陣A的行列式值。矩陣A=[(a11,a12)T·(a21,a22)T]的元素為

(16)

以及

(17)

式(16)和式(17)即通過無模糊相位值計算來波二維入射角的解析表達式。為得到無模糊相位值,還需進行模糊解算。

2.2 基于復數響應的相位模糊解算

(18)

式中:

w(n)=exp(jε(n))-1=x(n)+jz(n)

(19)

當SNR較高時,w(n)的實部和虛部根據泰勒級數展開進行一階近似,有

(20)

z(n)=sin(ε(n))?ε(n)

(21)

因此,與式(8)相對應的yp=[yp(1),yp(2),…,yp(N)]的條件概率密度可寫為

(22)

(23)

(24)

式中:round為四舍五入運算。然后通過模糊數即可計算出相位真值,即

(25)

3 理論精度下限分析

采用信噪比分析時變陣對頻變目標的測向誤差下限,作為評價時變陣測向性能的指標。由式(1)可得關于方位角和仰角的Fisher矩陣元素為

(26)

式中：α1=θp;α2=φp;μn=cos(φp(n));νn=sin(φp(n))。再根據式(5)得Fisher矩陣的各元素為

(27)

求Fisher矩陣的逆,即可得到合成圓陣測量寬帶時變信號的方位角與仰角的克拉米勞精度下限,分別為

(28)

(29)

值得注意的是,若采樣點均勻,各采樣信號頻率、SNR相同,分別為k0c/2π和SNR0時,可得到仰角、方位角的克拉米勞下限表達式,分別為

var(Δφp)≥(NM×SNR0)-1(k0dsinθp)-2

(30)

var(Δθp)≥(NM×SNR0)-1(k0dcosθp)-2

(31)

可見,來波入射角遠離陣面法向,仰角測向誤差增加,而方位角測向誤差減小。同時,測量相位的采樣點和場分布的采樣點數越多、SNR越高,測向誤差越小。

4 EM算法求解多信號參數

在多輻射源入射條件下,直接求解法既求下式的ML解：

(32)

式中：dis[x]=|mod(x,2π)|定義相位的距離,mod(x,2π)表示x對2π取余數運算。注意到,在信號未分選的情況下,尚無法確定相位值與輻射源的對應關系,以及輻射源的入射角,因此對式(32)的直接求解是一個計算復雜度為非確定性多項式的問題[37]。

換個角度,若已知輻射源的入射角信息,則可依據角度信息計算得到該輻射源在合成圓陣處產生的理論相位差,并將與該相位差最接近的采樣數據來自該輻射源,以此進行信號分選,從而建立輻射源與采樣數據之間的對應關系。另一方面,若已知信號分選,即已建立各接收數據與其輻射源輻的對應關系,則該P維問題變為P個一維問題,通過接收數據進行入射角解算。基于此將高維復雜問題降維的思想,本文先將觀察數據進行隨機分組,根據每組頻率、幅度、相位差等數據分別進行單輻射源的入射角估計,然后根據估計的入射角對觀測數據進行重新歸類,并依次在估計和歸類間進行迭代,直至算法收斂。合成陣列求解時變多目標參數的EM算法流程如圖3所示,主要包含求期望值步驟和期望值最大化兩個步驟。

4.1 求期望值步驟

求期望值步驟的作用是進行信號分選,即根據輻射源入射角計算該位置在合成圓陣處產生的相位差,并以采樣相位差距此理論相位差最小為依據進行接收數據的重新分組,將的采樣數據歸為一組中。根據第p個輻射源第i次入射角估計值,其在合成圓陣處所產生的理論相位差為

(33)

根據每個觀測相位差與理論相位差的距離,將觀測相位差與理論相位差最接近的數據分為同一組,即

(34)

最終將觀測數據分為P組。

4.2 期望值最大化步驟

此步驟采用ML法,通過分組后的觀測數據分別計算每組數據對應的二維入射角,即將求期望值步驟中每組觀測數據分別求解下式的ML解:

(35)

5 數值仿真實驗

本節通過數值實驗仿真分析無源合成圓陣的信號分選和測向性能,并驗證求解時變多目標信號參數的EM算法的能力。

5.1 單輻射源情形

本節主要考察測向誤差的均方根值隨入射仰角的變化情況。兩天線間距為1.2 m,旋轉一周采樣點數為100。信號的中心頻率為1.5 GHz,頻率變化范圍為中心頻率的13%,即從1.4～1.6 GHz范圍內隨機變化。SNR在-10～10 dB內隨機分布。仿真采用的每個接收信號的SNR與來波頻率如圖4所示。

仰角從5°變化到75°,方位角固定為170°,分別進行500次獨立計算,并將測向誤差的均方根與CRLB進行對比,重點考察了考慮每個接收信號SNR不同的情形。計算結果對比如圖5所示。考慮不同采樣點的SNR能提高測向精度,且考慮SNR的測向精度接近理論誤差下限。同時,入射仰角越大,方位角測向誤差越小,仰角測向誤差越大,與分析結論相符。

5.2 多輻射源情形

6 結 論

本文重點討論了采用無源合成陣列在多輻射源環境下對時變信號的參數估計問題,包括對多個輻射源的二維測向以及對其產生交錯信號的分選。建立了無源合成圓陣測向的數學模型,考慮了目標頻變、測向位置非均勻、單次測量信噪比不同等因素。提出了采用無源合成陣列進行時變多信號參數估計與分選的EM算法。該算法通過將高維多參數估計問題分解成多個低維參數估計問題,從而降低運算復雜度。根據基于相位數據的ML估計方法,推導出閉合形式且性能接近最優估計的入射角計算解析解。此外,采用基于接收響應復數的ML方法,利用旋轉帶來的基線多樣性實現了相位模糊的解算。通過兩種不同似然方法的有機結合,兼顧了入射角初值計算運算量與最終入射角的求解精度。通過理論下限推導,分析了無源合成圓陣的測向精度,揭示了入射仰角、陣列口徑、采樣點數和信噪比對二維測向精度的影響。最后通過單輻射源和多輻射源的數值實驗驗證了分析結論和算法的有效性。