張 敏，劉金彥，郭福成

(1.國防科技大學(xué)電子科學(xué)學(xué)院，湖南 長沙 410073； 2.武警警官學(xué)院，四川 成都 610200)

0 引言

估計信號來波方向(AOA)，又稱之為測向，在電子偵察、監(jiān)視、預(yù)警、通信等領(lǐng)域受到廣泛關(guān)注和研究[1-4]。其中基于相位干涉儀的測向技術(shù)，由于其原理簡單、測向精度高，應(yīng)用越來越廣泛[3 ]。當基線長度大于信號半倍波長時，相位差測量可能會出現(xiàn)2π模糊。解相位差的2π模糊是干涉儀測向研究中的主要問題[4]之一。

一種常用的方法是采用長短基線干涉儀形成多基線實現(xiàn)解模糊。如文獻[5]基于長短組合基線、文獻[6]基于剩余定理互質(zhì)基線和文獻[7]虛擬基線的解模糊方法，但這都要求多組基線之間有一定的幾何約束關(guān)系，因此限制了干涉儀基線的布陣，且通常僅適用于線陣，較難推廣到任意平面陣。

文獻[8]的無模糊長基線解模糊方法，僅適用于特定的平面陣型，且該方法中對相位差的多次代數(shù)運算造成了較大的測向誤差。文獻[9]聚類法和文獻[10]立體基線解模糊方法通?？蛇m用于任意平面陣，但此類方法僅由特定的幾組基線得到多個待定的角度估計，利用剩余基線在此待定角度中選取一個值作為估計值，不能達到測向理論精度。

文獻[11]相關(guān)干涉儀測向方法可適用于任何陣型干涉儀，且使用了所有基線用于測向，但是基于網(wǎng)格搜索的相關(guān)干涉儀測向方法需要將網(wǎng)格劃分得足夠小，才能得到較好的結(jié)果，特別是二維測向條件下，計算量大，實時性差。粒子群優(yōu)化(PSO)等智能方法雖然運算量有所降低，但由于采用隨機搜索策略，在粒子數(shù)量較少的情況下，全局尋優(yōu)能力并不是十分可靠[12]。

文獻[13]中的最小二乘(LS)測向適用于任何陣型干涉儀，且使用了所有基線用于測向，但要求相位差無模糊。在主動雷達中，文獻[14]中直接利用多個觀測站上的干涉儀測量的模糊相位差，對運動目標進行定位跟蹤。文獻[15]中利用旋轉(zhuǎn)干涉儀模糊相位差直接對固定輻射源定位。這說明直接利用模糊相位差進行測向估計具有可行性。

在此思路下，本文提出了一種多假設(shè)非線性最小二乘(Multiple Hypothesis Taylor Series，MHTS)干涉儀測向算法。該算法使用所有基線用于測向，無需事先解相位差模糊，直接利用模糊相位差進行測向，對干涉儀幾何構(gòu)型無特殊要求。

首先給出后文一些數(shù)學(xué)符號定義：T表示矩陣轉(zhuǎn)置運算，mod表示取模運算，·表示向下取整運算，|·|表示取絕對值運算，‖·‖表示取2范數(shù)運算，IN表示N×N維單位矩陣，0N表示N×N維零矩陣。

1 數(shù)學(xué)模型

1.1 測向模型

本文中方位角α∈[-π,π]定義為入射信號在XOY平面上的投影射線OT′與Y軸正方向的夾角；俯仰角β∈[0,π/2]定義為入射信號與其在XOY平面上的投影射線OT′之間的夾角，如圖1所示。

圖1 入射角度定義示意圖

對應(yīng)的干涉儀到輻射源的視線矢量為：

(1)

接收天線An1和An2構(gòu)成的干涉儀基線n長度為dn，方位安裝角為θn，對于二維平面陣，俯仰安裝角εn=0，對應(yīng)的基線指向矢量為：

(2)

因此不妨將視線矢量和基線矢量寫為：

(3)

(4)

由此可得干涉儀測量得到的無模糊相位差為：

φn=gnθ+δn，n=1,2,…,N

(5)

當干涉儀基線長度dn>λ/2時，相位差測量可能會出現(xiàn)2π模糊，因此有：

φn=gnθ+δnmod2π

(6)

式中，φn∈-π,π。

將無模糊相位差寫成矩陣形式：

Θ=Gθ+E

(7)

將對應(yīng)的模糊相位差的矩陣形式為：

Φ=Gθ+Emod2π

(8)

模糊相位差中包含了輻射源角度信息，干涉儀測向問題即是如何利用N組干涉儀基線測量得到模糊相位差，估計出方位和俯仰角。

1.2 相位差噪聲模型

多通道干涉儀如圖2所示，接收通道對天線接收的信號進行低噪聲放大、混頻、濾波、A/D等處理，經(jīng)過FFT等處理得到接收信號到達不同天線的初相位。通過對不同基線的兩個天線的初相位相減得到相位差。

圖2 多通道干涉儀示意圖

由于各天線和接收通道相互獨立，因此接收信號處理后得到的初相位之間可認為近似獨立的：

E(ξiξk)=0，i,k=1,2,…,M,i≠k

(9)

在此條件下，相位差為：

φik=ψi-ψk+δikmod2π

(10)

式中，ψi和ψk為天線i和天線k對應(yīng)的初相，δik=ξi-ξk對應(yīng)基線的相位差測量誤差。

由此可得相位差噪聲特性為：

(11)

式中，l,m=1,2,…,M,l≠m。

本文將此條件下的相位差噪聲模型稱為相位差的相關(guān)噪聲模型。

通過機械旋轉(zhuǎn)[16]或電子切換[17]等方式形成時變基線維干涉儀進行測向，通常僅需要2個接收通道，對通道幅/相一致性要求和系統(tǒng)復(fù)雜度大為降低。如圖3所示的雙通道切換干涉儀，在觀測時刻tn，通過M選2切換開關(guān)，選通2根天線構(gòu)成干涉儀基線。由于測量的相位差是在不同觀測時刻tn，利用不同接收信號在不同基線上得到的，因此測量誤差之間可近似看作獨立的。在此條件下，相位差為：

φn=ψn1-ψn2+δnmod2π

(12)

式中，ψn1和ψn2為天線1和天線2在觀測時刻tn對應(yīng)的初相，ξn1和ξn2為對應(yīng)天線初相測量誤差，δn=ξn1-ξn2對應(yīng)基線的相位差測量誤差。

圖3 雙通道切換干涉儀示意圖

由此可得相位差噪聲特性為：

(13)

式中，c=1,2,…,N。

本文將此條件下的相位差噪聲模型稱為相位差的獨立噪聲模型。

2 測向方法

Gθ≈Gθm+Pmθ-θm

(14)

式中，Jacobi矩陣Pm=?Gθm/?θm：

Pm=HJm

(15)

式中，H=[κb1,κb2,…,κbN]T，Jm為視線矢量xm關(guān)于θm的Jacobi矩陣：

(16)

同時利用初始值得到相位差預(yù)測值：

Θm=Gθm

(17)

從而得到相位差測量值和預(yù)測值之間的殘差：

zm=Φ-Θmmod2π

(18)

利用最小二乘可得：

(19)

需要注意的是，在上述計算中，只有當初值在真值附近時，(18)式中的2π取模運算才能去掉。但由于無法事先獲知這一先驗信息，受文獻[15]中利用干涉儀模糊相位差定位方法的啟發(fā)，若能在觀測空間內(nèi)選取多個初值，對每個初值采用Taylor級數(shù)方法計算，只要有一個初值落在真值附近，通過選取最小代價函數(shù)即可確定方位和俯仰角估計值。此即本文提出的多假設(shè)Taylor級數(shù)(MHTS)測向方法的基本思路。

下面給出測向算法的主要流程：

Step1：確定輻射源來波方向初始值θm=[αm,βm]T，m=1,2,…,M。初始時刻令m=1。

Step2：利用(15)式計算相位差關(guān)于θm的Jacobi矩陣Pm。

Step3：利用(17)式計算相位差預(yù)測值Θm，利用(19)式計算測量值和預(yù)測值之間的殘差zm。

Step6：計算估計值對應(yīng)的檢驗量Cm：

(20)

Step7：判斷m若等于M則轉(zhuǎn)入Step8，否則令m=m+1，轉(zhuǎn)入Step2繼續(xù)處理剩余初始值。

Step8：將檢驗量C(m)最小值對應(yīng)的估計值作為視線矢量的估計值：

(21)

3 初值選取方法

上述算法中，關(guān)鍵是如何獲得角度初始值。針對二維干涉儀陣的特點，進行初始值選取。首先選擇利用具有一定夾角(如90°)和一定基線波長比(如10～20之間或最長的幾組基線中)的2個模糊相位差φi和φl，以及方位安裝角θi和θl。

根據(jù)φi和φl對應(yīng)的基線長度di和dl確定模糊相位差各自的取值范圍k∈[-di/λ,di/λ]，h∈[-dl/λ,dl/λ]。通過遍歷模糊數(shù)k和h，根據(jù)相位差與來波方向的關(guān)系，計算得到初始值。

對每組模糊數(shù)k和h，由φi和φl以及對應(yīng)的方位安裝角θi和θl，計算無模糊相位差和系數(shù)矩陣：

(22)

(23)

根據(jù)最小二乘可得：

(24)

若‖xkh‖大于1，不能構(gòu)成單位視線矢量，因此舍去該組模糊數(shù)對應(yīng)的初值。否則計算得到對應(yīng)的方位和俯仰角的初始值：

(25)

(26)

由此可獲得M≤(2di/λ+1)(2dl/λ+1)組方位和俯仰角初值。該初始化方法可推廣到一維測向初值的選取。

4 仿真分析

通過計算機仿真，對不同干涉儀陣型以及一維和二維測向場景下，對本文測向方法性能進行驗證。為評估測向性能，采用Monte Carlo重復(fù)試驗統(tǒng)計定位誤差。仿真中，定義方位角和俯仰角的均方根(RMS)誤差為：

(27)

(28)

場景條件如下：一維測向中假設(shè)輻射源來波的方位角30°，信號頻率3GHz。二維測向場景中，假設(shè)輻射源來波的俯仰角32°，其余參數(shù)與一維測向場景中相同。本文測向算法最大迭代次數(shù)為3。

獨立噪聲模型下，兩個接收通道通過9選2電子切換開關(guān)在不同切換時刻接收輻射源信號，在每組基線停留時間為1ms。相關(guān)噪聲模型下，9個天線連接9個接收通道。各個通道測量得到初始相位，進而得到各組基線相位差。

仿真分析不同算法的二維測向性能。采用文獻[10]中聚類測向方法，仿真中采用了夾角為4π/9的4組最長基線用于聚類，其余基線用于解模糊，并將各組基線中與聚類結(jié)果最接近的角度取均值作為估計值。采用文獻[12]中的PSO方法用于相關(guān)干涉儀測向，利用本文方法中的代價函數(shù)進行，其中方位角搜索范圍為[-π,π]，俯仰角搜索范圍為[0,π/2]，粒子數(shù)量與本文方法初始點數(shù)量相同，最大迭代步數(shù)為400，慣性權(quán)重ω由0.9遞減到0.6，學(xué)習率c1和c2均為2.1。

仿真圖例中，Clstr表示聚類方法仿真結(jié)果； PSO表示PSO相關(guān)干涉儀方法仿真結(jié)果；MHTS表示本文算法仿真結(jié)果；CRLB表示測向誤差克拉美-羅下限，計算方法可參見文獻[18]；MSE表示測向理論均方誤差。仿真中，測向估計值與真值的誤差大于5倍理論誤差則判定為測向模糊。

1)仿真一：獨立噪聲模型下二維測向

圖4給出相位差獨立噪聲模型下，測量誤差3°～40°范圍內(nèi)，不同算法的二維測向解模糊概率。

圖4 獨立噪聲下不同算法二維測向解模糊概率

從圖4可以看出，PSO方法由于隨機搜索引起的個別局部極值，因此在相位差誤差3°～30°范圍內(nèi)，解模糊概率不能達到100%。聚類方法本質(zhì)上與相關(guān)方法相同，雖使用了多組基線，但仍未有效地使用所有基線。相位差誤差3°～40°條件下，本文的測向方法解模糊概率優(yōu)于聚類法和PSO方法。

圖5給出相位差獨立噪聲模型下，測量誤差3°～40°范圍內(nèi)，不同算法的二維測向精度。需要說明的是，圖5和后文中的圖7給出的統(tǒng)計誤差，是剔除野值后得到的RMS誤差。

圖5 獨立噪聲下不同算法二維測向性能

表1給出了二維測向方法仿真運算時間。

表1二維測向算法運行時間比較ms

算法相位差測量誤差/(°)51015Clstr9．189．189．17PSO212．8212．1214．7MHTS24．424．724．6

計算機操作系統(tǒng)為Microsoft Windows XP，硬件配置為Intel?CoreTMi5 CPU 750@2.67GHz , 2.66GHz，內(nèi)存3.24GB，使用Matlab?軟件進行計算。

從表1和圖5中可以看出， PSO方法運算時間高于其它方法一個數(shù)量級。聚類方法運算時間雖然優(yōu)于本文方法，但未有效地使用所有基線，因此測向性能不能達到CRLB。本文方法在相位差誤差3°～40°之間，方位角和俯仰角測量精度都可到達CRLB。

2)仿真二：相關(guān)噪聲模型下二維測向

圖6給出相位差相關(guān)噪聲模型下，測量誤差3°～30°范圍內(nèi)，不同方法的二維測向解模糊性能。

圖7給出相位差相關(guān)噪聲模型下，測量誤差3°～30°范圍內(nèi)，不同方法的二維測向精度。

從圖7可以看出，在相關(guān)噪聲模型下，在剔除野值后，聚類方法和本文方法的測向精度相當，且都可接近理論測向精度；PSO方法在初相測向誤差較小時，測向誤差略大于理論誤差。然而從圖6可以看出， PSO方法仍然由于隨機搜索引起的個別局部極值，導(dǎo)致解模糊概率不能達到100%；而本文方法的解模糊概率優(yōu)于其它兩種方法。

5 結(jié)束語

本文提出的一種適用于任意平面陣干涉儀的多假設(shè)Taylor級數(shù)(MHTS)測向方法，直接利用模糊相位差進行測向，無需事先解相位差模糊，對二維干涉儀幾何構(gòu)型無特殊要求。該方法利用相位差所包含的角度信息作為初值，采用Taylor級數(shù)方法進行計算，運算量適中。理論分析和數(shù)字仿真表明本文方法測向精度可接近理論測向精度?！?/p>

[1] Rong Yang, Pek Hui Foo, Boon Poh Ng, et al. RF emitter geolocation using amplitude comparison with auto-calibrated relative antenna gains [J]. IEEE Trans. on Aerospace and Electronic System, 2011, 47 (3): 2098-2110.

[2] Schmidt R. Multiple emitter location and signal parameter estimation[J]. IEEE Trans. on Antennas and Propagation, 1986, 34(3):276-280.

[3] Macphie RH, Tae Ho Yoon. On using the compound interferometer to obtain the power pattern of a conventional receiving array[J]. IEEE Trans. on Antennas and Propagation, 2009, 10(57):3356-3359.

[4] Ching-Sung Shieh, Chin-Teng Lin. Direction of arrival estimation based on phase differences using neural fuzzy network [J]. IEEE Trans. on Antennas and Propagation, 2000, 48(7): 1115-1124.

[5] Jacobs E，Ralston EW. Ambiguity resolution in interferometry [J]. IEEE Trans. on Aerospace and Electronic System, 1981, 17 (6): 766-780.

[6] 周亞強,皇甫堪. 噪擾條件下數(shù)字式多基線相位干涉儀解模糊問題[J].通信學(xué)報,2005,26 (8) :16-21.

[7] 曲志昱, 司錫才. 基于虛擬基線的寬帶被動導(dǎo)引頭測向方法[J]. 彈箭與制導(dǎo)學(xué)報,2007,27(4):92-95.

[8] 宋才水, 顧爾順. 無模糊長基線干涉儀測角的設(shè)計[J]. 現(xiàn)代防御技術(shù),2006,34(2):52-54.

[9] 張亮, 徐振海,熊子源,等. 基于圓陣干涉儀的被動導(dǎo)引頭寬帶測向方法[J]. 系統(tǒng)工程與電子技術(shù),2012,34(3):462-466.

[10] 張春杰, 李智東. 非均勻圓陣天線模型解模糊誤差研究[J]. 系統(tǒng)工程與電子技術(shù),2012,34(8): 1525-1529.

[11] Wei Hewen, Shi Yunguang. Performance analysis and comparison of correlative interferometers for direction finding [C]∥Proc. of the IEEE International Conference on Digital Object Identifier, 2010:393-396.

[12] Frans Van de Bergh. A new convergent particle swarm optimizer[C]∥IEEE International Conference on Systems, Man and Cybernetics, 2002:96-101.

[13] 翟慶偉. 最小二乘法測向技術(shù)研究[J].無線電工程,2008,38(3):55-57.

[14] Cheng Yongqiang, Wang Xuezhi, Caelli T. Tracking and localizing moving targets in the presence of phase measurement ambiguities[J]. IEEE Trans. on Signal Processing, 2011, 59(8): 3514-3525.

[15] 李騰,郭福成,姜文利. 基于旋轉(zhuǎn)干涉儀模糊相位差的多假設(shè)NLS定位算法[J].電子與信息學(xué)報,2012,34(4): 956-962.

[16] Kawase S. Radio interferometer for geosynchronous satellite direction finding[J]. IEEE Trans. on Aerospace and Electronic System, 2007, 43 (2):443-449.

[17] Wu Wei,Cooper CC, Goodman NA. Switched-element direction finding [J]. IEEE Trans. on Aerospace and Electronic System, 2009, 45(3): 1209-1217.

[18] Kay SM. Fundamentals of statistical signal processing, estimation theory [M]. Englewood Cliffs, New Jersey: Prentice Hall, 1993:30-36.