基于長基線干涉儀相位差的多站無源定位方法

張 敏 張文俊 李 曦 郭福成

(國防科技大學電子科學學院 長沙 410073)

1 引言

通過多個觀測站接收非合作輻射源信號從而確定輻射源位置的多站無源定位技術，在偵察、監視、預警、人員搜救等軍事和民用領域都具有重要應用價值[1-3]，受到了廣泛關注和研究。通常利用輻射源信號的到達時差(Time Difference of Arrival,TDOA)、到達頻差(Frequency Difference of Arrival,FDOA)、到達方向(Direction of Arrival, DOA)等一種或多種觀測量實現多站無源定位[4-6]。

根據觀測量的特點，多站無源定位方法大致可分為兩大類：一類方法中的定位觀測量需要觀測站之間進行相關處理才能獲得，如TDOA, FDOA。基于這類觀測量的多站無源定位方法在每個觀測站上僅需要單個接收天線和接收通道，系統組成簡單。但對于應用越來越廣泛的超低旁瓣輻射源，如有源電掃陣列(Active Electronically Scanned Array, AESA)雷達[7]，可能會出現兩個或多個觀測站無法同時截獲輻射源信號的情況，這就會導致多站時差頻差定位方法失效。

另一類方法所需觀測量在各自觀測站上就能實現測量，如DOA、相位差變化率(Change Rate of Phase Difference, CRPD)、截獲時間(Time of Interception, TOI)等。采用這類觀測量的定位方法可適應超低旁瓣輻射源定位。其中基于DOA的多站定位應用最為廣泛[8]，但高精度定位依賴于高精度測向。常用的比幅測向精度差，從而導致定位精度差；干涉儀比相或陣列測向精度與天線之間的間距(即基線長度)成正比，基線長度越大測向精度越高，為了實現高精度測向通常要求基線長度大于信號半倍波長，但這就會導致出現干涉儀或陣列出現測向模糊，從而導致顯著的定位誤差。在實際應用中每個平臺上需要多個(典型的5～7個)接收天線和通道實現無模糊高精度測向，這就會帶來造價成本和系統復雜度高的問題。基于CRPD的多站定位[9]可以將每個觀測站所需的接收天線和通道降低到至少3個，但要求基線長度足夠長(數米甚至更長)，并且需要累積一段時間(如0.5 s)才能實現高精度CRPD測量，在觀測站與輻射源之間存在高階相對運動情況下，會引入顯著的CRPD測量誤差導致定位性能下降。基于TOI的多站無源定位方法[10]，各個平臺上只需單個接收天線和通道，但僅適用于輻射源信號具有某種周期規律的特殊輻射源(如具有勻速掃描周期的機械掃描雷達)，通常至少需要1個周期才能實現定位。

近年來無人機(Unmanned Aerial Vehicle,UAV)等平臺技術應用越來越廣泛，這類平臺具有低成本、小型化等特點，基于這種小型化平臺編隊的多站無源定位具有重要理論和應用價值[1,2]。而現有的多站無源定位技術在這類平臺上的應用可能面臨實現難度大、定位精度差、造價成本和系統復雜度高等缺點。針對這些問題，本文提出了一種基于相位差(Phase Difference of Arrival, PDOA)的多站無源定位新體制，在每個觀測站上僅需至少兩個接收天線和通道構成的長基線干涉儀(Long Baseline Interferometer, LBI)，通過測量輻射源信號到達LBI接收天線的PDOA實現定位，可在降低系統復雜度的同時實現高精度定位，并且能夠適應對超低旁瓣輻射源的定位應用。

為了實現高精度定位，觀測站上的PDOA基線長度通常需要遠大于信號半倍波長，這就會導致測量的PDOA出現 2π模糊。針對PDOA觀測量與輻射源位置之間的非線性和非連續性，提出了一種多假設迭代優化定位方法，首先利用兩個觀測站的PDOA獲得多個可能的輻射源位置初始值，然后采用高斯-牛頓(Gauss-Newton, GN)方法對每個位置初始值進行迭代優化并計算代價函數，最后選擇具有最小代價函數的優化結果作為最終的定位結果。計算機仿真結果表明該方法運算量適中，在測量噪聲為高斯分布且測量噪聲適中條件下，定位精度可達到克拉美羅下限(Cramer-Rao Lower Bound,CRLB)。

2 定位模型

2.1 數學模型

假設地面有一固定輻射源向空中輻射信號。N個觀測站經過輻射源所在區域，接收輻射源輻射的信號，并利用其上的LBI測量得到PDOA，根據各個觀測站上的LBI基線矢量和觀測站位置，通過優化計算確定出輻射源位置。多站無源定位場景示意圖如圖1所示。

圖1 定位場景示意圖

在圖1所示定位場景中，假設信號輻射源位置為xt，輻射源信號載頻為fc，對應的信號波長為λ。第n個觀測站的位置為xn，其上的LBI基線矢量為bn。第n個觀測站的LBI測量得到的輻射源信號到達LBI兩個接收天線的無模糊PDOA為

其中，n=1,2,...,N,κ=2π/λ,rn=(xt-xn)/rn為輻射源到第n個觀測站的單位視線矢量，rn=||xtxn|| 為輻射源到第n個觀測站的距離，wn為PDOA測量噪聲，假設其滿足均值為零、方差為σn2的高斯分布。

將N個觀測站的無模糊PDOA表示為向量形式

其中，m=[?1,?2,...,?N]T,f(xt)=[f1(xt),f2(xt),...,fN(xt)]T,fn(xt)=κbTnrn,δ=[w1,w2,...,wN]T。

當LBI的基線長度大于輻射源信號半倍波長時，實際測量得到的PDOA存在 2π模糊，模糊PDOA為

其中， mod2π(·) 表示 2π取模運算，取模運算后的取值范圍為 -π～π。

將N個觀測站的模糊PDOA表示為向量形式

需要注意的是PDOA模型中基線矢量bn的計算通常會涉及不同的坐標系(如地固坐標系、北東地坐標系、本體坐標系等)之間的相互轉換，具體的坐標系定義和轉換可參見文獻[11]。

2.2 定位原理的幾何解釋

單個觀測站上的LBI測量的模糊PDOA，可以確定多條以該觀測站為頂點、以干涉儀基線矢量所在直線為旋轉軸的定位圓錐面(當LBI與輻射源在同一平面時，定位圓錐面退化為以觀測站為起點的定位射線)，定位圓錐面與地面相交得到一簇定位曲線，每條定位曲線對應模糊PDOA的一個模糊整數，其中一條定位曲線穿過輻射源所在位置，如圖2(a)所示。

圖2 定位原理示意圖

當存在多個觀測站時，由于觀測站之間的位置各不相同，觀測站上LBI的基線指向與信號入射方向之間的夾角也會隨之改變，這種變化就會引起不同位置觀測站上的LBI測量的模糊PDOA在地面上形成與其他觀測站不同的一簇定位曲線。不同觀測站之間的定位曲線相互交匯，不僅會在輻射源位置形成定位點，還會在其他位置上形成大量的模糊定位點(Ghost)，但只有輻射源位置上交匯的定位曲線的數量最大且等于觀測站數量，如圖2(b)所示。

對于圖2所示的定位場景，理論上至少需要3個觀測站的LBI測量的模糊PDOA才能實現無模糊定位。但受觀測噪聲的影響，往往需要大于3個觀測站測量的模糊PDOA，并通過一定的優化計算方法才能無模糊地(唯一地)確定輻射源位置。這就是本文提出的基于長基線干涉儀相位差的多站無源定位原理的幾何解釋。

2.3 PDOA的測量

各個觀測站的模糊PDOA通常可通過對LBI兩個天線接收的信號進行相關處理實現測量。圖3給出了一種典型的PDOA測量的基本原理[9]：

圖3 PDOA測量原理示意圖

在圖3中，天線A1和A2接收到信號s1(t)和s2(t)后送到對應的接收通道。每個接收通道利用低噪聲放大器(Low Noise Amplifier,LNA)提高處理增益，通過本振(Local Oscillator, LO)和下變頻器(Down Converter, D/C)進行混頻、濾波等處理后得到中頻信號x1(t)和x2(t)。將接收通道1和通道2中的中頻信號進行傅∫里葉變換

其中， angle(·)為對復數取相位函數，函數的取值范圍為 -π～π，ωmax為信號功率譜 |Y(ω)|的最大值對應的頻率。

3 定位方法

模糊PDOA不僅與輻射源位置之間存在非線性，還存在由于 2π取模運算導致的非連續性。這就使得基于模糊PDOA的多站定位解算面臨更為嚴峻的挑戰，難以獲得類似于多站測向定位中的偽線性最小二乘估計(Pseudolinear Least Square,PLS)或多站時差定位中的兩步法等閉合解析定位解[4-6,8]。網格搜索(Grid Search,GS)、粒子群優化(Particle Swarm Optimization,PSO)等智能優化方法對非線性非連續性優化問題求解具有通用性[12,13]，但GS方法要獲得接近CRLB的定位性能所需的運算量非常大，PSO方法難以保證收斂到最優結果。

通過對模糊PDOA進行正弦或余弦以及正余弦預處理[14]可消除LBI的PDOA中 2π取模運算的影響，但由于觀測噪聲的存在會在預處理觀測模型中引入模型誤差從而導致定位性能下降。本文直接利用模糊PDOA的概率密度函數建立最大似然估計，在此基礎之上提出了一種多假設迭代優化定位方法。

3.1 最大似然估計

其中，? =z-f(xt)為 相位差殘差向量，W=E(δTδ)=diag(σ12,σ22,...,σN2)為N個觀測站的PDOA測量誤差的協方差矩陣。

可通過最大化關于輻射源位置xt的對數似然函數lnp(z|xt) 得 到輻射源位置的ML估計x?t，即

其中，基于模糊PDOA的ML代價函數為

在第5節的仿真條件下，圖4給出了PDOA測量誤差為10°，在X軸范圍為–100～100 km, Y軸范圍為–100～100 km的區域內，利用式(13)得到的ML代價函數(圖中顯示的是代價函數的倒數)對應的網格圖和等值圖。

圖4 定位代價函數示意圖

從圖4(a)可以看到在輻射源位置真實值附近具有最高峰(對應于最小的代價函數)，這就意味著利用多站LBI的模糊PDOA可以唯一地確定出輻射源位置。但從圖4(b)可以看到在輻射源其他位置上也可能存在次高峰，這是由于PDOA的 2π模糊造成的模糊定位點。這就意味著在定位方法中需要考慮如何有效地消除 2π取模運算的影響。

3.2 必要條件

通常來講由ML代價函數得到式(13)的加權最小二乘(Weighted Least Square, WLS)定位優化目標函數，可通過采用高斯-牛頓(GN)方法等方法解算輻射源位置。而GN方法這類迭代優化方法，首先需要給出輻射源位置初始值xt0。下面分析輻射源位置初始值對ML代價函數中 2π取模運算的影響。將第n個觀測站上LBI測量的無模糊PDOA在輻射源位置初始值xt0附近進行Taylor級數展開并保留1階項有

當輻射源位置初始值xt0接近于輻射源位置真實值xt，且PDOA噪聲足夠小，即滿足式(15)必要條件

這就說明只要選擇合理的輻射源位置初始值，就可以通過模糊PDOA得到與無模糊PDOA等效的定位結果。

3.3 多假設迭代方法

由于輻射源是非合作的，因此在實際中很難直接獲得接近于輻射源位置真實值xt的初始值xt0。針對這個問題，本文提出了一種多假設迭代定位方法，首先在輻射源所在的區域內選擇多個可能的位置初始值，然后對每一個位置初始值假設其位于輻射源真實值附近，采用GN方法計算輻射源位置的ML估計，并將ML代價函數作為統計檢驗量，最后選擇具有最小代價函數的優化結果作為最終的定位結果。在滿足式(15)的必要條件下，只要存在一個輻射源初始值位于真實值附近，就可以得到與無模糊PDOA等價的且具有最小代價函數的定位結果。定位方法流程如圖5所示。

圖5 定位方法流程示意圖

假設可以選擇K個輻射源初始值，根據GN方法對第k(k=1,2,...,K)個輻射源位置初始值xt0,k進行優化計算得到迭代更新量為[16]

3.4 初始值選取方法

在本文提出的多假設迭代優化定位方法中，需要解決的一個關鍵問題是如何選擇多個可靠的輻射源位置初始值。如果選擇的位置初始值數量過多，則會造成計算量增大，算法效率降低；如果位置初始值數量過少，則可能造成位置初始值遠離位置真實值，從而導致無法收斂到全局最優。文獻[17]提出的一種基于觀測區域劃分的初始化方法，將每個子區域的中心點作為定位初始值，雖然與網格搜索方法相比這種方法獲得的初始點數量顯著降低，但沒有充分利用模糊PDOA包含的位置信息進一步降低定位初始點的數量。

從定位原理的幾何解釋中可以看到，對于單個觀測站上的LBI測量的模糊PDOA會產生多條可能的(模糊的)定位線，兩個觀測站上的LBI測量的模糊PDOA則會產生多個可能的(模糊的)定位點。雖然不能確定輻射源位于其中具體的哪一個定位點，但這已經將輻射源位置約束在數量有限的備選結果之中。因此利用兩個觀測站的模糊PDOA就可以高效地產生定位方法所需的輻射源位置初始值。對于3維場景下的2維定位(估計輻射源的xy坐標或經度緯度)問題，初始值選取方法如圖6所示。

圖6 初值選取方法示意圖

模糊PDOA對應的定位曲線退化為多條以LBI為起點的定位射線，其中一條定位射線穿過輻射源所在位置。假設選擇第p和q個觀測站，根據LBI基線長度可分別確定模糊PDOA的最大模糊數分別為

利用θap和θaq就可以根據僅測向定位中經典的偽線性最小二乘(PLS)等方法[8]得到輻射源位置初始值。對于3維場景下的3維定位(估計輻射源的xyz坐標或經度緯度高程)問題，可選擇3個觀測站的PDOA，采用類似的方法獲得輻射源位置初始值。

4 定位性能分析

4.1 定位誤差的CRLB

在滿足式(15)的必要條件下，模糊PDOA得到的定位結果等價于無模糊PDOA得到的定位結果，因此本文以無模糊相位差下的定位誤差CRLB作為評價基于模糊PDOA的多站無源定位精度指標。根據CRLB定義可得[19]

其中，p(xt|m) 為無模糊PDOA關于xt的條件概率密度函數，在高斯噪聲模型條件下有

其中，Jxt=J(xt)|xt為無模糊PDOA關于輻射源位置真實值xt的雅可比矩陣。

由此可得

4.2 定位方法運算量分析

首先從乘法運算次數的角度對所提方法單次迭代所需運算量進行分析；然后從相位差計算次數的角度分析本文方法與粒子群優化、網格搜索兩種典型定位方法的運算量。

本文方法中的GN方法單次迭代所需的乘法(含除法)運算主要涉及式(20)—式(22)。注意到加權矩陣W為對角矩陣，因此式(20)需要 12N+10次乘法運算；式(21)計算模糊PDOA需要 10N次乘法運算(不含N次開根號運算)；式(22)計算雅可比矩陣需要 21N次乘法運算(不含N次開根號運算)，因此單次GN迭代需要 43N+10次乘法運算。假設選取的輻射源位置初始值數量K=(2Aq+1)(2Ap+1)，假設每個輻射源位置初始值的最大迭代次數為I，則所提方法單次迭代所需的乘法運算次數為I(43N+10)(2Aq+1)(2Ap+1)。

從相位差計算次數的角度來看，本文方法所需的相位差計算次數為NI(2Aq+1)(2Ap+1)。對于粒子數量為NPSO，最大迭代步長為IPSO的PSO方法，最大的相位差計算次數為NPSOIPSON。在2D問題的GS方法中，兩個維度上的搜索范圍分別為W1和W2，需要達到的精度為C1和C2，網格數量分別為NG1=W1/C1和NG2=W2/C2，對應的相位差計算次數為NG1NG2N。根據輻射源載頻、LBI基線長度、觀測站數量，以及算法收斂準則等條件的不同，3種方法對應的運算量也可能隨之變化。

5 仿真分析

采用3維場景下的2維定位應用對所提方法性能進行仿真分析。假設16個觀測站均勻分布在邊長12 km的正方形區域內，區域中心位于坐標原點正上方100 m，每個觀測站上的LBI基線長度為1 m；輻射源位于 [56, 60]Tkm，輻射源信號載頻為5 GHz。觀測站和輻射源的在地面上投影的幾何構型如圖7所示。圖中觀測站上的箭頭表示平臺運動方向(也是基線指向)。

圖7 仿真場景示意圖

首先分析了PDOA噪聲標準差(Standard Deviation, STD)為5°～35°條件下[17]基于模糊PDOA、余弦PDOA、正弦PDOA、正余弦PDOA模型[14]下的定位均方根誤差(Root Mean Square Error,RMSE)，蒙特卡羅仿真次數為1 000次，GN最大迭代次數為3次。為了進行性能對比，還給出了無模糊PDOA模型的RMSE以及CRLB。圖8(a)給出了各個觀測站PDOA噪聲具有相同STD(σn2=σ2，σ為仿真中設定的PDOA噪聲的STD)條件下的仿真結果。圖8(b)給出了各個觀測站PDOA具有不同STD條件下的仿真結果，其中假設第n個觀測站PDOA測 量 誤 差 的S T D 為σn=σ(1+0.5 cos(2πn/N))。

圖8 不同PDOA模型下定位RMSE

從圖8(a)可以看到在各站PDOA噪聲具有相同STD條件下，正余弦PDOA、模糊PDOA的定位誤差RMSE與無模糊PDOA的RMSE相同且接近CRLB；從圖8(b)可以看到在各站PDOA噪聲具有不同STD條件下，僅有模糊PDOA的RMSE與無模糊PDOA的相同并接近于CRLB，這驗證了本文提出的定位數學模型的有效性。同時注意到在圖8(a)和圖8(b)中正弦PDOA和余弦PDOA的RMSE都偏離CRLB，且各站PDOA的STD差異越大偏離程度越高。這是由于正弦PDOA、余弦PDOA以及正余弦PDOA預處理方式雖然消除了模糊PDOA中的不連續性，同時也不同程度地引入了模型誤差，從而導致定位性能下降。

從圖9(a)和圖9(b)可以看到本文提出的定位方法的RMSE可接近CRLB，且算法耗時接近粒子數量為100的PSO定位方法；PSO定位方法的RMSE在PDOA噪聲的STD較小時難以達到CRLB，粒子數量越多RMSE越接近CLRB，但算法耗時會隨粒子數量的增大而增長；網格搜索方法的精度受網格寬度的影響難以達到CRLB，但算法耗時隨著PDOA噪聲的增大而降低；隨著PDOA噪聲的增大，模糊定位點代價函數可能低于真實位置對應的代價函數，從而導致定位RMSE偏離CRLB。

圖9 不同定位方法的RMSE

6 結語

本文提出了一種基于PDOA的多站無源定位新體制，利用每個觀測站上長基線干涉儀測量輻射源信號到達觀測站的PDOA實現定位。針對模糊PDOA的非線性和非連續性，提出了一種基于多假設迭代優化定位方法，通過對選取的多個可能的輻射源位置，直接利用模糊PDOA進行GN優化計算確定出輻射源位置。仿真結果表明所提算法運算量適中，在高斯觀測噪聲條件下定位精度可達到CRLB。