地坐標意義下三維單站無源定位迭代算法

王妙妙，江怡帆，趙 悅，姬利海

(1.北京工業大學,北京 100190；2.北京應用物理與計算數學研究所,北京 100094；3.多倫多大學，多倫多 M5S 2E8；4.東北師范大學,吉林 長春 130024)

0 引 言

無源定位技術能在自身不輻射的條件下,隱蔽確定輻射源位置,具有作用距離遠、抗干擾能力強的特點, 對于提高系統在電子戰環境下的生存能力和作戰效能具有十分重要的作用，因此無源定位已經成為電子對抗的一個相當重要、不可或缺的技術。雷達輻射源無源定位方面現有成果主要針對地面直角坐標系意義下開展的工作[1]。眾所周知, 在過去的幾十年中，無源定位技術的研究不斷深入,主要包括定位精度和效率的分析。其中, 最流行的方法有到達角(AOA),到達時間差(TDOA)和到達頻率差(FDOA),可參考文獻[2]～ [6]。文獻[7]引入了相位差變化率定位算法(定位速度較快, 但相位差變化率測量精度需提高), 但該方法對硬件要求比較高,需要很高的時差測量精度，同時該方法僅適用于大樣本定位, 因此對于當前現代相控陣雷達的無源偵察很難有效，對于功率比定位，現有技術功率測量精度很有限。由于空對地偵察過程中受到諸多硬件條件限制, 而基于AOA的單站定位是利用一個平臺上單個或多個接收機(絕大部分為單個)在不同時刻對信號方向角測量和處理獲得目標方位經緯度信息, 具有經濟性、機動性、靈活性, 同時不需要對外數據通信等優點, 所以基于單偵察站(AOA)無源定位是空對地偵察中應用最廣的定位方法[8-10]。遺憾的是, 現有算法基本都是基于直角坐標系推導而來, 本文核心算法解決了現有體制采樣設備與無源定位算法直接對接問題, 在實際地坐標條件下討論問題, 同時充分考慮了地球橢球曲率對定位模型的影響。本文給出了實際偵察環境地坐標條件下小樣本和大樣本無源定位算法, 具體包括三維情形下, 兩點定位及理論誤差估計算法、大樣本廣義最小二乘迭代定位算法, 從理論上小樣本兩點定位確定的解可以作為大樣本迭代定位的迭代初值。

下一節給出三維地坐標與地面直角坐標間的轉換關系，第2節給出了地坐標條件下小樣本定位算法, 該算法可以為后面大樣本條件下迭代算法提供基本的迭代公式, 同時可以為迭代算法提供一個比較好的初值, 兩點定位誤差估計采用等概率曲線來確定誤差協方差矩陣所對應的誤差橢圓所覆蓋的定位概率。直角坐標和地坐標條件下廣義最小二乘迭代算法在第3節中給出，數值模擬在第5節中給出, 本節從數值角度分析了第2節中兩點交叉定位和第3節中2種迭代定位的定位精度, 同時分析了針對快速移動目標的定位收斂速度、定位精度與目標移動速度之間的依賴關系。第5節給出了總結。

1 三維地坐標與地面直角坐標間的變換

地球的實際形狀為橢球, 以地球的球心o1為坐標原點,B點為零度經緯線的交點之一。以o1B所在直線為x1軸的方向, 赤道所在平面上與x1軸垂直的方向為y1軸所在方向,o1z是地心軸所在的直線。赤道接近圓,且半徑為r=a=6 378 137 m。y1o1z1所在平面對應的橢圓短半軸為：b=6 356 752 m。

圖1 地坐標與地面直角坐標關系示意圖

(1)

如圖1所示,設以A點為坐標圓點的三維直角坐標架為x2y2z2o2(o2與A重合),其中o2y2為過A點的正北方向,o2x2為過A點的正東方向,o2z2過A點鉛直于地面。結合公式(1),若僅考慮直角坐標系x1y1z1o1與x2y2z2o2之間的平移關系, 則：

(2)

本文中其它地方也用類似的表示方法。另外,顯而易見,兩坐標架之間還存在著依賴于α,γ的旋轉變換,結合公式(2)與旋轉變換,一般地,直角坐標系x1y1z1o1與x2y2z2o2之間的坐標變換關系為：

(3)

(4)

2 三維地坐標小樣本定位算法

設兩偵察點位P1(α1,γ1,h1),P2(α2,γ2,h2)及其相應的2個點位的偵察目標到達方位角和俯仰角為(θ1,β1),(θ2,β2),基于這些已知條件求偵察目標的經緯度(αR,γR,hR)。三維地坐標小樣本定位算法基于直角坐標系下的交叉定位算法,因此需借助上一節中的坐標變換，先把經緯度坐標變換為直角坐標下的形式。由式(1)知,P1(α1,γ1,h1)在x1y1z1o1坐標架中的坐標為:

(5)

(6)

類似地,P2(α2,γ2,h2)在x1y1z1o1坐標架中的坐標為：

(7)

(8)

(9)

(10)

由公式(9)～(10)可知：

(11)

(12)

在直角坐標系下三維交叉定位模型為:

(13)

由公式(6)～(8),(11)～(13)可聯立進行兩點交叉定位計算(只需用公式(13)中的3個表達式,這時不失一般性, 假設用前3個表達式來交叉定位)。

3 三維地坐標大樣本定位算法

接下來討論三維地坐標條件下廣義最小二乘迭代定位算法：

(14)

(15)

根據公式(14),有：

(16)

從式(12)可得到：

(17)

(18)

因此,公式(16)和(17)變成：

(19)

hR)cosαRcosγR-

(20)

(21)

經一階泰勒展開，上式可簡化為：

(22)

(23)

經一階泰勒展開，上式可簡化為：

-NPicosθicos2γi-(1-

2e2)NPitanβicosγisinγi+

(NR+hR)(1-e2)sinγRtanβicosγi-

(24)

(25)

若樣本量m較大時, 可通過下面推廣形式的最小二乘估計, 滿足公式(21)的待估參數使得采樣點擬合殘差最小:

(26)

(27)

結合傳統的最小二乘估計表達式可給出下面的迭代估計算法：

(4) 循環迭代至收斂。

4 數值模擬

例1(靜態偵察目標情形)，設偵察機飛行速度為(300,200)km/h, 初始坐標為(28.5°E,37.3°N),初始到達角為θ=π/3,偵察機初始高度為h=8 km。設每隔dt=0.01 s采集1組數據, 數據量為200, 共采集時間長度為T=1 s。并設采集到樣本的誤差服從N(0,0.12)的正態分布。假設最小二乘方法初始迭代步數為200, 迭代初值選取為(28.5°E,37.3°N,0.8 km), 迭代收斂誤差限制為10-6，則不同目標位置所需要的迭代步數、迭代結果以及相對誤差如表1～表3所示。

在以下的數值實驗中,定義相對誤差如下:

(28)

(29)

表1 目標位置為(126°E,37°N,1 km)的估計結果

表2 目標位置為(130.49°E,30.38°N,0.5 km)的估計結果

另外, 我們選取2種不同的迭代初值, 定位同一個目標(30°E, 40°N, 1 km)，定位收斂速度、定位精度如表4和表5所示。

表3 目標位置為(127°E,40°N,0.3 km)的估計結果

表4 初值為(20°E,30.30°N,0.8 km)的估計結果

表5 初值為(10°E,60°N,4 km)的估計結果

同時, 我們仍然選擇同一個定位目標(30°E, 40°N, 1 km), 選擇的初值分別為小樣本交叉定位坐標與任意坐標, 從定位結果表5和表6可以看到, 利用小樣本交叉定位的解作為迭代算法的初值收斂速度更快, 迭代所需要的樣本更少, 因此在針對現代相控陣雷達的無源定位偵察過程中很有必要引入該迭代定位思路。

從表4～表6結果可以發現, 本文提出的廣義最小二乘算法對迭代初值的依賴不強, 在另一個層面上表明本文所提算法具有很好的穩定性。

表6 初值選取為到達角(AOA)計算得到的值

例 2(動態情形)，假設偵察機的初始參數設置和例1中一樣, 不失一般性，假設目標初始位置坐標為(30°E,40°N,1 km), 并設其勻速直線運動速度為VR，則在不同移動速度、不同時刻、相同的迭代初值(20°E,30°N,0.8 km)情況下, 利用最小二乘方法計算結果如表7～表10所示。

表7 速度VR=30 km/h的估計結果

表8 速度VR=45 km/h的估計結果

表9 速度VR=72 km/h的估計結果

表10 速度VR=100 km/h的估計結果

從表7～表10的數值模擬結果可以看出。本文算法對地面移動目標的定位性能。算法的定位收斂速度很快, 收斂誤差及相對誤差都比較小，定位精度很高, 由此可見本文算法可以成功應用于實際應用中。

5 結束語

本文算法主要應用于對靜態目標或慢速移動目標的定位,并且能夠達到較高的定位精度。對于移動目標定位可以從數值分析得到滿足不同精度需求時本文算法能定位的目標速度的大小, 基本可以滿足實際應用中對靜態目標和地面快速移動目標的定位需求。缺點是計算太復雜, 影響計算速度，硬件化不太方便等,以后可進一步通過大量例子進行模擬和改進。