載機(jī)平臺(tái)光電轉(zhuǎn)塔目標(biāo)定位的仿真算法

秦 川，吳玉敬，陶 忠，桑 蔚，安學(xué)智

（西安應(yīng)用光學(xué)研究所，陜西 西安 710065）

引言

機(jī)載光電目標(biāo)定位任務(wù)主要由載機(jī)、光電轉(zhuǎn)塔和GPS/慣導(dǎo)（INS）三大系統(tǒng)組成。其中GPS/INS提供載機(jī)的經(jīng)緯高信息和姿態(tài)信息（包括真航向、俯仰角、橫滾角），光電轉(zhuǎn)塔完成對(duì)目標(biāo)的搜索、捕獲、定位和跟蹤等功能，其瞄準(zhǔn)目標(biāo)后可提供瞄準(zhǔn)線（LOS）的方位、俯仰角等信息（涵蓋光電轉(zhuǎn)塔方位、俯仰角和目標(biāo)橫向、縱向像素偏差角），結(jié)合傳感器本身參數(shù)（包括視頻顯示的方位、俯仰視場(chǎng)角和視頻顯示的橫向、縱向分辨率）即可實(shí)現(xiàn)對(duì)海平面目標(biāo)的無源定位任務(wù)。若結(jié)合安裝于光電轉(zhuǎn)塔內(nèi)的激光測(cè)距機(jī)給出的測(cè)量信息，則可實(shí)現(xiàn)三維空間目標(biāo)的有源定位任務(wù)。

在海面目標(biāo)的無源定位過程中，目標(biāo)與載機(jī)的距離信息并非直接測(cè)量，因此需要：1）用載機(jī)高度的負(fù)值近似目標(biāo)在載機(jī)東北天坐標(biāo)系下的zG值，進(jìn)而求解距離R的近似算法；2）結(jié)合地球橢球模型，進(jìn)而求解距離R的嚴(yán)格解析算法。此外，無論GPS/INS 還是轉(zhuǎn)塔視軸穩(wěn)定度、視頻像素偏差讀取、激光測(cè)量距離等都存在一定誤差，這些誤差近似呈Gaussian 分布，所以每次定位結(jié)果都會(huì)或多或少偏離目標(biāo)真實(shí)位置，這在實(shí)際定位過程中是需要盡可能降低的，并且定位精度作為光電轉(zhuǎn)塔核心作戰(zhàn)技術(shù)指標(biāo)之一也需要預(yù)先估計(jì)出來。

文獻(xiàn)[1]對(duì)車載光電偵查系統(tǒng)的目標(biāo)定位模型及誤差源進(jìn)行了分析并通過外場(chǎng)試驗(yàn)加以驗(yàn)證；文獻(xiàn)[2]分析了多光譜圖像的定位與配準(zhǔn)，特別是對(duì)地理定位誤差源進(jìn)行了探討；文獻(xiàn)[3]對(duì)定位精度的評(píng)價(jià)指標(biāo)進(jìn)行研究，指出均方根誤差、平均誤差及置信概率區(qū)間方法的計(jì)算原則及使用條件；文獻(xiàn)[4]分析了測(cè)量平臺(tái)中的坐標(biāo)轉(zhuǎn)換誤差，給出載機(jī)姿態(tài)精度、位置精度提升對(duì)定位結(jié)果的影響大?。晃墨I(xiàn)[5]分析了光電偵查系統(tǒng)的定位精度影響因素，提出對(duì)INS 的誤差控制建議；文獻(xiàn)[6]分析了機(jī)載光電觀瞄系統(tǒng)的目標(biāo)定位原理，對(duì)算法進(jìn)行了試驗(yàn)研究；文獻(xiàn)[7]基于激光測(cè)距，對(duì)有源目標(biāo)定位技術(shù)進(jìn)行了研究；文獻(xiàn)[8]～[10]也對(duì)相關(guān)定位問題進(jìn)行了討論。

本文基于載機(jī)平臺(tái)光電轉(zhuǎn)塔目標(biāo)定位需求，分析了有源、無源目標(biāo)定位的算法，其中對(duì)于海面目標(biāo)無源定位，目標(biāo)載機(jī)距離的計(jì)算采用嚴(yán)格的解析算法；對(duì)于定位精度的預(yù)估，本文提出轉(zhuǎn)塔視軸反演算法，并結(jié)合前者定位算法及Monte-Carlo 分析方法，從仿真角度得出一定測(cè)量精度條件下的海面目標(biāo)無源定位精度，即圓概率誤差（CEP）半徑約39.4 m@80%置信度及約18 km 距離；此外還做了載機(jī)位置、姿態(tài)、轉(zhuǎn)塔視軸不同時(shí)對(duì)目標(biāo)定位的80%置信度CEP 半徑、最大誤差距離等統(tǒng)計(jì)直方圖分析。本文的仿真方法及結(jié)果對(duì)于目標(biāo)定位研究具有一定的理論指導(dǎo)意義。

1 基本定義

1）地球橢球

在控制測(cè)量中，代表地球的橢球稱作地球橢球，它是地球的數(shù)學(xué)模型。地球橢球的幾何定義為：O是橢球中心，NS（南北極）為旋轉(zhuǎn)軸，a為長(zhǎng)半軸（取值6 378 137 m），b為短半軸（取值6 356 752 m）。赤道定義為：通過橢球中心的平行圈。緯圈定義為：垂直于旋轉(zhuǎn)軸的平面與橢球面相截所得的圓。子午圈定義為：包含旋轉(zhuǎn)軸的平面與橢球面相截所得的橢圓。其第一偏心率和第二偏心率分別為

2）載機(jī)機(jī)體直角坐標(biāo)系（OA-XAYAZA）

載機(jī)機(jī)體坐標(biāo)系為三維直角坐標(biāo)系，坐標(biāo)系的原點(diǎn)OA為載機(jī)質(zhì)心，XA軸為載機(jī)右翼正向，向外為正，YA軸為載機(jī)縱軸機(jī)頭正向，ZA軸通過右手螺旋定則確定，垂直XAYAZA平面朝向機(jī)身上方。

3）東北天直角坐標(biāo)系（OG-XGYGZG）

東北天坐標(biāo)系為三維直角坐標(biāo)系，坐標(biāo)系的原點(diǎn)OG為載機(jī)質(zhì)心，XG軸指向正東，YG軸指向正北，ZG軸沿地垂線指向天空。

4）地心直角坐標(biāo)系（Od-XdYdZd）

地心直角坐標(biāo)系為三維直角坐標(biāo)系，原點(diǎn)Od為參考橢球體中心，Zd軸與參考橢球體旋轉(zhuǎn)軸重合，向上為正，Xd軸為格林尼治子午面與地球赤道平面的交線，向外為正，Yd軸在赤道平面內(nèi)，并與Xd、Zd成右手直角坐標(biāo)系。

5）相對(duì)直角坐標(biāo)系（OR-XRYRZR）

相對(duì)坐標(biāo)系是為了實(shí)現(xiàn)東北天坐標(biāo)系到地心直角坐標(biāo)系的轉(zhuǎn)換，自主定義的三維直角坐標(biāo)系。其原點(diǎn)OR為載機(jī)質(zhì)心在地球表面的投影，XR軸與Od-Xd平行且方向一致，YR軸與OR-YR平行且方向一致，ZR軸與OR-ZR平行，且方向一致。

6）空間大地坐標(biāo)系（WGS-84，LBH）

空間大地坐標(biāo)系就是WGS-84 世界大地坐標(biāo)系，是美國國防部于1984年定義，通過大地經(jīng)度（L）、大地緯度（B）和大地高（H）來描述空間位置。經(jīng)度是空間中的點(diǎn)與參考橢球的自轉(zhuǎn)軸所在的面與參考橢球的起始子午面的夾角，緯度是空間的點(diǎn)與參考橢球面的法線與赤道面的夾角，大地高是空間點(diǎn)沿參考橢球的法線方向到參考橢球面的距離。

7）載機(jī)姿態(tài)（θεγ）

載機(jī)姿態(tài)包括真航向（θ）、俯仰角（ε）和橫滾角（γ），采用慣導(dǎo)進(jìn)行測(cè)量。真航向定義為飛機(jī)所在位置的真經(jīng)線北端順時(shí)針測(cè)量至航向線YA軸的夾角，北偏東為正。俯仰角定義為機(jī)體YA軸線與地平面間的夾角，以飛機(jī)抬頭為正。橫滾角定義為機(jī)體ZA軸與包含機(jī)體YA軸的鉛垂面間的夾角，以飛機(jī)向右傾斜為正（尾追觀察）。

8）光電轉(zhuǎn)塔方位角（Tx）和俯仰角（Ty）

以載機(jī)坐標(biāo)系為基準(zhǔn)，光電轉(zhuǎn)塔光軸偏離XAYA平面的角度稱為俯仰角，光軸抬頭為正；光軸在XAYA平面的投影偏離YA軸的角度稱為方位角，投影向右偏離為正（俯視觀察）。

9）目標(biāo)橫向像素偏差角（σx）和縱向像素偏差角（σy）

目標(biāo)在視頻顯示屏上距離光軸的像素偏差，向右向上為正。

10）目標(biāo)的穩(wěn)瞄瞄準(zhǔn)線俯仰角（α）和方位角（β）

目標(biāo)的穩(wěn)瞄瞄準(zhǔn)線俯仰角和方位角分別定義為

關(guān)于2012年的水利規(guī)劃計(jì)劃工作，陳雷部長(zhǎng)提出了明確要求，學(xué)文同志進(jìn)行了具體的布置。這里，我再強(qiáng)調(diào)以下三項(xiàng)工作的落實(shí)。

式中：Vx、Vy分別為視頻顯示的方位和俯仰視場(chǎng)角；Rx、Ry分別為視頻顯示的橫向和縱向分辨率。

2 坐標(biāo)系轉(zhuǎn)換

2.1 目標(biāo)的瞄準(zhǔn)線俯仰角、方位角和載機(jī)機(jī)體直角坐標(biāo)系的轉(zhuǎn)換

顯然，目標(biāo)的穩(wěn)瞄瞄準(zhǔn)線俯仰角、方位角到載機(jī)機(jī)體直角坐標(biāo)系間轉(zhuǎn)換關(guān)系為[11]

式中：R表示目標(biāo)距載機(jī)距離，對(duì)于有源目標(biāo)定位，該值來自激光測(cè)距機(jī)等的信息。對(duì)于無源海平面目標(biāo)定位，該值可以采用兩種方法獲得：其一是近似方法，其二是嚴(yán)格解析方法。這兩種方法都先假設(shè)（5）式中R為未知變量，然后根據(jù)近似條件或者解析條件算出R，最后再帶入（5）式中重新計(jì)算目標(biāo)的載機(jī)機(jī)體直角坐標(biāo)值。

載機(jī)機(jī)體直角坐標(biāo)系到穩(wěn)瞄瞄準(zhǔn)線俯仰角、方位角間轉(zhuǎn)換關(guān)系為

2.2 載機(jī)機(jī)體直角坐標(biāo)系和東北天直角坐標(biāo)系的轉(zhuǎn)換

載機(jī)機(jī)體直角坐標(biāo)系到東北天直角坐標(biāo)系間轉(zhuǎn)換關(guān)系為[11]

其中轉(zhuǎn)換矩陣A為

對(duì)于海平面目標(biāo)，可近似假設(shè)zG=-H，從而求出R的近似表達(dá)式：

東北天直角坐標(biāo)系到載機(jī)機(jī)體直角坐標(biāo)系的轉(zhuǎn)換關(guān)系為

2.3 東北天直角坐標(biāo)系和相對(duì)直角坐標(biāo)系的轉(zhuǎn)換

東北天直角坐標(biāo)系到相對(duì)直角坐標(biāo)系的轉(zhuǎn)換關(guān)系為[11]

其中轉(zhuǎn)換矩陣C為

相對(duì)直角坐標(biāo)系到東北天直角坐標(biāo)系間轉(zhuǎn)換關(guān)系為

2.4 相對(duì)直角坐標(biāo)系和地心直角坐標(biāo)系的轉(zhuǎn)換

相對(duì)直角坐標(biāo)系到地心直角坐標(biāo)系間轉(zhuǎn)換關(guān)系為[11]

式中N為載機(jī)所在橢球卯酉圈曲率半徑，計(jì)算公式為

對(duì)于海平面目標(biāo)，其位于地球橢球面上，滿足：

因此，可通過求解含有R的一元二次方程得到目標(biāo)與載機(jī)的距離信息解析值。

地心直角坐標(biāo)系到相對(duì)直角坐標(biāo)系間轉(zhuǎn)換關(guān)系為

2.5 空間大地坐標(biāo)系和地心直角坐標(biāo)系的轉(zhuǎn)換

空間大地坐標(biāo)系向地心直角坐標(biāo)系的轉(zhuǎn)換關(guān)系為[12]

式中：LT、BT、HT分別為目標(biāo)的經(jīng)度、緯度和高度；NT為目標(biāo)所在橢球卯酉圈曲率半徑。

地心直角坐標(biāo)系向空間大地坐標(biāo)系的轉(zhuǎn)換關(guān)系為[13]

式中：BT和HT的求解含有交叉項(xiàng)，可以通過zd/sqrt（xd2+yd2）先求出BT的初值，然后利用BT的初值結(jié)合（15）式和（19）式，即可以求出NT、HT的初值，接著再次求解BT的值，循環(huán)往復(fù)，直至滿足精度要求[13]。BT、HT也可以采用解一元四次方程的方法[12]，或者Bowring 解法[14]：

3 仿真分析

3.1 算法流程

對(duì)于目標(biāo)定位，即是已知載機(jī)經(jīng)、緯、高（LBH）、載機(jī)姿態(tài)（θεγ）、光電轉(zhuǎn)塔瞄準(zhǔn)線方位角（Tx）和俯仰角（Ty）、目標(biāo)橫向像素偏差角（σx）和縱向像素偏差角（σy）、視頻顯示的方位和俯仰視場(chǎng)角（Vx和Vy）、視頻顯示的橫向和縱向分辨率（Rx和Ry），求解目標(biāo)的空間大地坐標(biāo)（LTBTHT）。若目標(biāo)與載機(jī)的距離R已知，則是有源定位；若未知，則是無源定位。對(duì)于有源定位，目標(biāo)不必限定位于地球表面；而對(duì)于無源定位，若只通過載機(jī)和光電轉(zhuǎn)塔一組測(cè)量數(shù)據(jù)進(jìn)行解算，必須限定目標(biāo)位于地球表面。

與定位過程相反的過程，即是光電轉(zhuǎn)塔視軸反演過程，即已知載機(jī)LBH、θεγ，光電轉(zhuǎn)塔σxσy、VxVy、RxRy，目標(biāo)LTBTHT等信息，求解光電轉(zhuǎn)塔的TxTy。

本文采用R的解析求解方法和BT、HT的直接求解方法，其算法流程如圖1所示。其中圖1（a）過程對(duì)應(yīng)光電轉(zhuǎn)塔視軸反演，圖1（b）過程對(duì)應(yīng)有源目標(biāo)定位，圖1（c）過程對(duì)應(yīng)無源目標(biāo)定位。

圖1 目標(biāo)定位與反演算法流程Fig.1 Flow diagram of target location and inversion algorithms

3.2 光電轉(zhuǎn)塔視軸反演

假設(shè)目標(biāo)在（東經(jīng)108°53′53″，北緯34°12′45″，高度426.6 m），傳感器視場(chǎng)0.45°×0.25°，分辨率1 920×1 080 像素，載機(jī)飛行速度為400 km/h，運(yùn)動(dòng)軌跡如圖2所示，即載機(jī)在6 000 m 高度先向北平飛，然后逆時(shí)針做3/4 圓周盤旋，接著開始向東爬升高度至7 000 m 后平飛，姿態(tài)依照飛行航跡相應(yīng)調(diào)整。

圖2 載機(jī)航跡與姿態(tài)Fig.2 Aircraft tracks and attitudes

采用圖1中反演算法，在10 min 時(shí)間內(nèi)，目標(biāo)距離載機(jī)的距離R如圖3（a）所示，光電轉(zhuǎn)塔的方位角和俯仰角分別如圖3（b）和3（c）所示。

圖3 反演結(jié)果Fig.3 Inversion results

3.3 誤差源分析

依據(jù)第2 節(jié)和第3.1 節(jié)的分析，誤差主要來自GPS 測(cè)量載機(jī)經(jīng)緯高的誤差、INS 測(cè)量載機(jī)姿態(tài)的誤差、光電轉(zhuǎn)塔視軸穩(wěn)定精度帶來的誤差、目標(biāo)橫向與縱向像素偏差角的視頻判別誤差、測(cè)距機(jī)測(cè)量誤差等。上述誤差，可以看成許多量的疊加，根據(jù)中心極值定理，從統(tǒng)計(jì)意義上講是滿足Gaussian分布的隨機(jī)誤差，其大小如表1所示。此外，還存在GPS、INS 和光電轉(zhuǎn)塔的安裝誤差，在載機(jī)結(jié)構(gòu)剛性假設(shè)下，這三者是固有誤差，可通過地面校軸加以減小，本文暫不列入分析。對(duì)于傳感器視場(chǎng)角和分辨率，也可通過地面檢驗(yàn)進(jìn)行測(cè)量，本文暫將其列為固定值。

3.4 有源目標(biāo)定位

當(dāng)目標(biāo)與載機(jī)距離用測(cè)距機(jī)測(cè)量得到時(shí)，采用圖1中有源目標(biāo)定位算法，計(jì)算得到的目標(biāo)經(jīng)度、緯度和高度分別如圖4（a）、4（b）和4（c）所示。其中藍(lán)色曲線是10 min 內(nèi)單點(diǎn)有源定位直接結(jié)果，綠色曲線是上述直接結(jié)果在1 s 內(nèi)做均值濾波得到，紅色直線代表目標(biāo)真實(shí)位置。從結(jié)果中可以看到：1）經(jīng)度、緯度和高度直接定位結(jié)果的最大-最小誤差分別為0.001 304°、0.001 630°、64.119 4 m；2）做1 s 濾波后，經(jīng)度、緯度和高度定位結(jié)果的最大-最小誤差分別為0.000 093°、0.000 118°、5.083 5 m；3）目標(biāo)定位的經(jīng)緯高是概率收斂到其真實(shí)位置。

表1 誤差源及其大小Table1 Error source and its intensity

圖4 有源目標(biāo)定位結(jié)果Fig.4 Active target location results

3.5 無源目標(biāo)定位

由于單點(diǎn)無源目標(biāo)定位理論所限，其要求目標(biāo)位于海平面上，即假設(shè)第3.2 節(jié)中目標(biāo)的高度為0，接著，類似第3.2 節(jié)反演過程計(jì)算光電轉(zhuǎn)塔的視軸，然后對(duì)加入噪聲的觀測(cè)量采用圖1中無源目標(biāo)定位算法，計(jì)算得到目標(biāo)經(jīng)度、緯度，分別如圖5（a）、5（b）所示，目標(biāo)定位位置與真實(shí)位置之間的距離如圖5（c）所示。其中藍(lán)色曲線是10 min 內(nèi)單點(diǎn)無源定位直接結(jié)果，綠色曲線是上述直接結(jié)果在1 s 內(nèi)做均值濾波得到的，圖5（a）、5（b）中紅色直線代表目標(biāo)真實(shí)位置。從結(jié)果中可以看到：1）經(jīng)度、緯度直接定位結(jié)果的最大-最小誤差分別為0.002 466°、0.001 966°；2）做1 s 濾波后，經(jīng)度、緯度定位結(jié)果的最大-最小誤差分別為0.000 170°、0.000 139°；3）目標(biāo)定位位置與真實(shí)位置之間的距離1 s 均值濾波前后分別為129.969 7 m 和44.495 5 m；4）目標(biāo)無源定位結(jié)果比有源定位結(jié)果差。

圖5 無源目標(biāo)定位結(jié)果Fig.5 Passive target location results

4 目標(biāo)定位精度分析

4.1 Monte-Carlo 方法

蒙特卡洛（Monte-Carlo）方法亦稱為隨機(jī)模擬（random simulation）方法，有時(shí)也稱作隨機(jī)抽樣（random sampling）技術(shù)或統(tǒng)計(jì)試驗(yàn)（statistical testing）方法。該方法是先建立一個(gè)概率模型或隨機(jī)過程，然后通過對(duì)模型或過程的抽樣試驗(yàn)來計(jì)算統(tǒng)計(jì)特征，最后給出所求解的近似值[15]。概率論中的大數(shù)定律和中心極限定理是蒙特卡洛方法的理論基礎(chǔ)。

4.2 圓概率誤差仿真

目標(biāo)定位的圓概率誤差定義為：以目標(biāo)真實(shí)位置為圓心劃一圓圈，若定位點(diǎn)位于圓圈內(nèi)的概率為P，圓圈的半徑為r，則稱定位結(jié)果在P置信度下的CEP 半徑為r。

假設(shè)載機(jī)在（東經(jīng)108°58′，北緯34°5′，高度6 000 m），姿態(tài)水平，光電轉(zhuǎn)塔方位北偏西30°，俯仰水平向下20°，像素偏差為0，傳感器視場(chǎng)及分辨率同前文，各誤差源的大小同表1，采用無源定位算法，設(shè)定1 000 次樣本實(shí)驗(yàn)，結(jié)果如圖6所示。圖中，紅色圓圈代表目標(biāo)真實(shí)位置，此時(shí)目標(biāo)距離載機(jī)17.606 km，綠色和藍(lán)色十字代表1000 次無源定位計(jì)算得到的目標(biāo)位置，其中綠色十字代表接近目標(biāo)真實(shí)位置占比80%的定位結(jié)果，其半徑為39.404 m，亦即CEP 半徑=39.404 m@80%置信度，最差定位結(jié)果與目標(biāo)真實(shí)位置相距80.463 m。

圖6 無源目標(biāo)定位蒙特卡洛分析結(jié)果Fig.6 Passive target location results by Monte-Carlo analysis

4.3 統(tǒng)計(jì)直方圖仿真

假設(shè)載機(jī)在（東經(jīng)108°58′，北緯34°5′），載機(jī)高度范圍[1 000，8 000]，航向范圍（-180，180]，俯仰范圍[-45，45]，橫滾范圍[-45，45]，轉(zhuǎn)塔方位范圍（-180，180]，俯仰范圍[-75，-15]，其他參數(shù)同前。限定目標(biāo)至載機(jī)的距離在20 km 內(nèi)為合理取值，選取1 000 個(gè)樣本點(diǎn)，則有810 個(gè)樣本的載機(jī)-目標(biāo)距離在合理范圍，其分布如圖7（a）所示。目標(biāo)推算位置與真實(shí)位置的距離落于80%范圍內(nèi)的半徑最小值為8.59 m，最大值為117.82 m，其大多數(shù)情況位于40 m 內(nèi)，位于30 m 內(nèi)有726 個(gè)樣本，占有效樣本數(shù)的89.63%，其分布如圖7（b）所示。目標(biāo)推算位置與真實(shí)位置的距離最大半徑最小值為18.25 m，最大值為372.28 m，其大多數(shù)情況位于100 m 內(nèi)，其分布如圖7（c）所示。

圖7 統(tǒng)計(jì)直方圖分析結(jié)果Fig.7 Analysis results of statistical histogram

5 結(jié)論

本文對(duì)于光電目標(biāo)定位從理論到仿真進(jìn)行了研究：首先分析了目標(biāo)定位中用到的坐標(biāo)系及其相互轉(zhuǎn)換關(guān)系，然后給出光電轉(zhuǎn)塔視軸反演、有源目標(biāo)定位、無源目標(biāo)定位的算法流程，接著通過仿真實(shí)驗(yàn)加以驗(yàn)證，其間考慮了定位過程中可能的隨機(jī)誤差來源，并分析了是否采用均值濾波對(duì)定位結(jié)果的影響，最后通過Monte-Carlo 分析計(jì)算了定位精度。通過給定條件的仿真實(shí)驗(yàn)，得到以下結(jié)論：1）有源定位比無源定位的精度高（精度約提高1 倍）；2）均值濾波后，定位精度有較大提升（約提高15 倍）；3）約18 km 距離典型無源定位精度在80%置信度條件下約為39.4 m；4）載機(jī)高度、姿態(tài)、轉(zhuǎn)塔視軸等在一定范圍內(nèi)采樣統(tǒng)計(jì)，限定目標(biāo)載機(jī)距離在20 km 內(nèi)為合理取值，其80%置信度CEP 半徑從8.59 m～117.82 m，最大誤差為18.25 m～372.28m。