基于最小二乘法的橢圓擬合實時航跡矢量化方法

基于最小二乘法的橢圓擬合實時航跡矢量化方法

民航業(yè)的快速發(fā)展使得各地空管部門壓力急劇增大，為了解決日益嚴峻的空中交通管制問題，本文提出實時航跡矢量化方法。該方法的原理為首先對實時雷達數據采用最小二乘法擬合航跡段，然后通過判定航跡段間的夾角，將航跡處理為直線段與曲線段的連接；最后通過橢圓擬合方法對實時航跡點預測，以及相應的偏差修正，以完成對于實時航跡的矢量化預測過程。對飛行航跡進行擬合的方法以最小二乘法最為普遍，最早的橢圓擬合方法源自最小二乘法，考慮隨機產生的誤差對數據的影響同時將這些影響以平方的形式最小化。通過實驗計算近千條真實飛行航跡數據，得出基于最小二乘法的橢圓擬合實時航跡矢量化方法的誤差在可控的范圍內，是一種有效，可行的實時航跡矢量化預測方法。

航空運輸行業(yè)的快速發(fā)展，為人們的出行提供了便利。但是隨著空中交通流量的持續(xù)增長，航空延誤、空域擁堵等問題時有發(fā)生，使得航空器沖突探測與解脫、進離場航班排序、 基于軌跡運行等空管自動化與智能化方法成為空中交通管理領域的研究重點，而如何快速準確的進行航空器飛行實時軌跡預測是實現上述方法的基礎與保障。

目前飛行軌跡預測的方法主要有兩種。第一種是基于混合估計理論實現的航跡預測，利用交互式多模型（IMM）算法通過對橫向、縱向以及垂直方向三種狀態(tài)估計加權求和實現隨機線性混雜系統的狀態(tài)估計，進而預測出航跡。第二種是基于航空器動力學及運動學模型， 利用各類機型的性能參數實現航跡的實時預測。如文獻提出根據飛行階段特點，用基本飛行模型構建水平航跡、高度剖面和速度剖面， 根據航跡特征點的飛行狀態(tài)信息擬合生成完整的 4D 航跡。

上述提出的兩類方法用于飛行軌跡的預測存在以下兩種問題：一是航空器的飛行模擬狀態(tài)需要隨實際情況的變化而變化，而IMM算法中假設的模擬轉換矩陣為固定值。二是在不考慮氣象信息或不能準確判斷航空器飛行意圖的情況下，很難保證航跡預測的準確度。

為了解決上述問題，本文中根據已知的飛機航行軌跡，判斷當前所處階段為直線航行或曲線航行，估測后續(xù)飛行軌跡點，進而擬合出行進軌跡的方法。最后，再對因軌跡預測可能產生的偏差分情況分別進行糾正，形成一條符合實時環(huán)境狀況、飛行器飛行特征和飛行意圖的平滑的行進軌跡。

方法原理

擬合原理

最小二乘法最早被稱為“回歸分析法”，由著名統計學家道爾頓于1806年所創(chuàng)，當時作為一種數學優(yōu)化技術，解決了日常生活中的計算問題。它通過最小化誤差的平方和尋找數據的最佳函數匹配。在現在，通過這種方法可以對未知數據進行預測，并且極大化地減小誤差，求得與真實數據誤差最小的預測數據。常用來解決曲線擬合、最佳逼近、物體運動軌跡等問題。其他一些優(yōu)化問題也可通過最小二乘法來解決。

當存在大量實驗數據時，不能絕對要求擬合函數?（x）在數據點（xi，yi） 處的誤差，即為0，但為了使曲線盡量逼近真實數據點的變化趨勢，將誤差的平方和最小，這就是最小二乘法的原理。

航跡的分段

就典型的航空器飛行航跡來說，航空器一般都是以直線的飛行軌跡從一個航跡點飛行至另一個航跡點，或者在空中進行了轉彎等飛行動作，以曲線的飛行軌跡從一個航跡點飛行至另一個航跡點，這些航路點信息可以從領航計劃報（filed flight plan message，FPL）中得知。即航跡可由一系列的直線飛行航跡段和曲線飛行航跡段組成。

直線飛行航跡段是飛行過程中的最簡單的部分，如圖1所示，P1（x1，y1，z1）點和P2（x2，y2，z2）點分別為一條直線飛行航跡段的起點和終點，在空間中兩點之間的直線距離為

曲線飛行航跡段是由調整航路或航班避讓以及進場盤旋等原因形成。通常，轉彎可分為三種模型：內切轉彎，末端轉彎和約束轉彎。

航路上最常使用的為內切轉彎，如圖2所示。其中點 p 1（x 1， y 1）為轉彎初始點 與p 2 （x 2 ， y 2） 轉彎結束點。

故在實際航跡的預測中，需判斷得出一段航跡點中的直線段與曲線段，將相應段的航跡點進行劃分，從而得到相應的直線段和曲線段擬合方程。具體航跡分段方法如下：

在典型的水平飛行航跡中，若有A，B，C三個相鄰的航跡點，其中A（Xa，Ya，Za），B（Xb，Yb，Zb），C（Xc， Yc，Zc），則線段AB的直線方程為，斜率線段BC的直線方程為，斜率。故兩直線間的夾角θ有P=tanθ=，則θ=arctanP。若此時θ< 10°，則A，B，C三點間用直線連接，如圖3所示。若此時θ> 10°，則A，B，C三點間用曲線連接，如圖4所示。

實際航跡預測過程中的分段過程如圖5所示。若此時有a，b，c，d，e，f，g，h，8個航跡點。第一步首先判斷a，b，c三點，經過計算可得到直線ab與直線bc間的夾角大于170°，故判斷a到c點之間以直線連接；第二步，計算直線ac與直線cd之間的夾角，經計算得出，夾角小于170°，故此時的處理過程為，返回c點之前的點b，此時bd段可判斷為曲線；第三步，計算直線cd與直線de之間的夾角，經計算得出，夾角小于170°，此時ce段可判斷為曲線，又因為bd屬于曲線段，故此時be段屬于曲線部分；第四步，計算直線de與直線ef之間的夾角，經計算得出，夾角大于170°，此時e點為曲線和直線的連接點，故ef段屬于直線，采用相同的方法計算直線ef和直線fg之間的夾角，直線fg，和直線gh間的夾角，得出eh段屬于直線部分。此時如圖5所示，將整條航跡分為了三個部分，即直線部分ab段，曲線部分be段和直線部分eh段。

圖1　直線飛行航跡

圖2　內切轉彎示意圖

圖3　用直線連接

圖4　用曲線連接

圖5　航跡分段

圖6　直線航跡點的預測

實時航跡矢量化方法

在非實時航跡的矢量化中，可以通過對已知所有航跡點進行擬合，從而得到整條非實時的航跡；但是在實時航跡的矢量化中，僅通過預緩存的若干點進行最初階段航跡擬合，然后對下一時間點的實時航跡位置進行預測，在得到真實的該時間點航跡位置后，又采用相應偏差糾正方法對航跡進行平滑糾正。航跡點的預測主要分為直線航跡點的預測以及曲線航跡點的預測（此處我們用橢圓方程進行曲線擬合）兩類。

直線航跡點的預測

方法1：由于所預測的航跡點都是以4s為時間間隔，故在短時間內的運動皆近似為勻速運動處理。如圖6所示。

其中A、B、C、D分別為時間間隔4秒的真實航跡點，E為預測的下一個航跡點。具體預測方法為：采用最小二乘法擬合得到直線方程，其中以航跡點的緯度和經度分別定義坐標軸的橫縱坐標；因該段過程我們采用勻速運動處理，故（其中Xi表示i點的橫坐標，i為A、B、 C、D、E），將Xe的值代入擬合得到的直線方程中，得到E點的縱坐標Ye，綜上我們得到預測點E的坐標（Xe，Ye）。

曲線航跡點的預測

此處，我們采用橢圓方程對曲線進行擬合。

（1）焦點在坐標軸上

方法2：由于所預測的航跡點都是以4s為時間間隔，故在短時間內的運動皆近似為勻速運動處理。如圖7所示。

這是教師針對性批改作文的做法，是“授之以漁”的一個重要方法。因為每一個學生的習作，在選材、內容、結構都不一定相同，甚至相同的題材，側重點也不盡相同，好、中、差學生的習作也肯定不同。我們可以選擇優(yōu)良的習作面批面改，在全班范讀，指出好在哪里，起典范作用。若為了全面提高，也要選一些其他層次的學生習作，和學生面對面地指導，教給學生修改的方法。

其中，A、B、C、D分別為時間間隔4s的真實航跡點，E為預測的下一個航跡點。又因為該橢圓的焦點分別位于X軸Y軸上，故該橢圓的方程可采用參數方程：

其中a，b分別為橢圓的長短半軸，θi為i點的參數（i為A、B、C、D、E點）。原點O的坐標（0，0），設A點坐標（Xa，Ya），D點坐標（Xd，Yd），E點坐標（Xe，Ye）。若向量=（Xa，Ya），=（Xd，Yd）間的夾角為?，則，故將參數θe代入參數方程，即可得到E點的坐標。

（2）橢圓焦點不在坐標軸上

方法3：由于所預測的航跡點都是以4s為時間間隔，故在短時間內的運動皆近似為勻速運動處理。如圖8所示。

其中，A、B、C、D分別為時間間隔4s的真實航跡點，E為預測的下一個航跡點，O'點坐標為O'（m，n）。其中，橢圓的焦點不在坐標系XOY上。虛線所表示的為原始坐標系XOY，實線所表示的是經旋轉和平移后的坐標系X'O'Y'，兩坐標系橫軸的夾角（即旋轉角）為ν。故，兩坐標系之間的轉化關系為

又因為該橢圓的焦點分別位于旋轉平移后的X' 軸Y'軸上，故該橢圓的方程可采用參數方程

分別擬合。其中a，b分別為橢圓的長短半軸，θi為i點的參數（i為A、B、C、D、E點）。原點為O'，若向量間的夾角為?，則故將參數θ代入參數方程，

e即可得到E點的坐標。

圖7　焦點在坐標軸上曲線航跡點的預測

圖8　焦點不在坐標軸上曲線航跡點的預測

圖9　直線段到直線段的過渡偏差糾正

圖10　曲線段到曲線段的過渡偏差糾正

實時航跡糾正偏差方法

通過實驗發(fā)現當飛機行駛方向發(fā)生偏轉時，由于預測具有延遲性，預測點與真實點存在偏差，如果立即按照新預測的航跡線進行飛行，就會造成整體的航跡線在轉彎處轉彎角過大。為了解決上述問題，從直線段到直線段的過渡、曲線段到曲線段的過渡、直線段到曲線段的過渡及直線段到曲線段的過渡以上四種情況對航跡進行糾偏，使得經過糾偏處理后的數據擬合出的航跡更加平滑。

（1）實時航跡直線段到直線段過渡偏差修正

如圖9所示，按照擬合完成的直線A'B'方程，飛機由A'點飛至B'點，將A'B'段的航跡預測看做理想狀態(tài)，故此處的實際航跡AB段與預測航跡A'B'重合。在B'點之后，按照本文第三部分中所介紹的實時航跡點預測方法1，預測得到4s后的航跡點C'，但飛機經過4s后由B'點飛行至C'點，此時雷達傳來真實的航跡點C點。為了將模擬航跡平滑過渡至真實航跡，采用以下方法：用C點以及該時刻向后倒推4s，8s，12s的航跡點，采用最小二乘法擬合出直線方程“f1=k1x +b1”；采用相同的方法，擬合出關于C'的直線方程“f2=k2x +b2”，當時，令當時，令所預測的下一時刻實時航跡點D'坐標，由此時所得到的直線D'B'的方程“f=kx +b”按照實時航跡點預測方法1得出。不斷地迭代此過程，直至預測航跡與實際航跡接近乃至重合。故整個預測模擬航跡為點A'B'C'D'E'F'的連線。

（2）曲線段到曲線段的過渡偏差修正

如圖10所示，按照擬合完成的橢圓1方“A1x2+B1xy +C1y2+D1x+E1y+F1=0”，飛機由A點飛至D點，將A'D'段的航跡預測看做理想狀態(tài)，故此處的實際航跡AD段與預測航跡A'D'重合。在D'點之后，按照本文第三部分中所介紹的實時航跡點預測方法2或3，預測得到4s后的航跡點E'，但飛機經過4s后由D'點飛行至E'點，此時雷達傳來真實的航跡點E點。為了將模擬航跡平滑過渡至真實航跡，采用以下方法：用E點以及該時刻向后倒推4s，8s，12s的航跡點，采用最小二乘法擬合出橢圓2方程時，令當時，令其中a=A，B，…，F 。預測的下一時刻實時航跡水平面坐標點F'由橢圓3的方程得出。不斷地迭代此過程，直至預測航跡與實際航跡接近乃至重合。故整個預測模擬航跡為點A'B'C'D'E'F'G'的連線。

（3）曲線段到直線段的過渡偏差修正

如圖11所示，按照擬合完成的曲線AD方程，飛機由A點飛至D點，將AD段的航跡預測看做理想狀態(tài)，故此處的實際航跡AD段與預測航跡A'D'重合。在D'點之后，按照本文第三部分中所介紹的實時航跡點預測方法2或3，預測得到4s后的航跡點E'，但飛機經過4s后由D'點飛行至E'點，此時雷達傳來真實的航跡點E點。為了將模擬航跡平滑過渡至真實航跡，采用以下方法：用E點以及該時刻向后倒推4s，8s，12s的航跡點，采用最小二乘法擬合出直線方程“f1=k1x +b1”；采用相同的方法，擬合出關于E'的直線方程“f2=k2x +b2”，當k1-k2＜0.1時，令k=k1，當時，令所預測的下一時刻實時航跡點F'坐標，由此時所得到的直線C'F'的方程f=kx +b按照實時航跡點預測方法1得出。不斷地迭代此過程，直至預測航跡與實際航跡接近乃至重合。故整個預測模擬航跡為點A'B'C'D'E'F'G'的連線。

（4）直線段到曲線段的過渡偏差修正

如圖12所示按照擬合完成的直線AC方程，飛機由A點飛至C點，將A'C'段的航跡預測看做理想狀態(tài)，故此處的實際航跡AC段與預測航跡A'C'重合。在C'點之后，按照本文第三部分中所介紹的實時航跡點預測方法1，預測得到4s后的航跡點D'，但飛機經過4s后由C'點飛行至D'點，此時雷達傳來真實的航跡點D點。為了將模擬航跡平滑過渡至真實航跡，采用以下方法：用E點以及該時刻向后倒推4s，8s，12s的航跡點，采用最小二乘法擬合出曲線段C'E'的橢圓方程，采用相同的方法擬合出曲線段CE的橢圓方程時，令當時，令其中a=A，B，…，F。預測的下一時刻實時航跡水平面坐標點E'由曲線段CE'（兩條曲線中所夾的曲線）橢圓的方程得出。不斷地迭代此過程，直至預測航跡與實際航跡接近乃至重合。故整個預測模擬航跡為點A'B'C'D'E'F'的連線。

圖11　曲線段到直線段的過渡偏差糾正

圖12　直線段到曲線段的過渡偏差糾正

圖13　實時航跡直線段到直線段的偏差糾正

圖14　實時航跡曲線段到曲線段的偏差糾正

實驗驗證

本實驗以首都機場為例，采用首都機場提供的2013年6、7月的飛機飛行數據及機場地理位置信息，隨機抽取1000條離場航跡進行實驗。首先，依據實時航跡預測方法對數據進行處理預測出實時航跡，在平面坐標系內進行效果對比。利用預測出的航跡與實際航跡之間數據的剩余方差，判斷預測航跡與實際航跡之間的誤差，從而驗證方法的有效性。

通過實驗前后對比圖可看出，該實時航跡預測方法預測效果較好，可滿足實時航跡預測的要求。

航跡的預測精度即航跡的預測位置與實際位置之間的差值 。通過計算該預測航跡與實際航跡的穩(wěn)定程度來近似估算預測方法的精度 具體來說主要步驟如下：

1）對預測航跡或實際航跡S，假設按時間順序該航跡由點P1，...，Pm組成，對這些點的緯度x，經度y值分別做差分得到兩組序列X=｛X1，...，Xm-1｝，Y=｛Y1，...，Ym-1｝，其中X1=Pi+1（x）-Pi（x），Yi=Pi+1（y）-Pi（y），通過序列X，Y把目標在x，y軸上的運動自動分解成一系列的勻速運動和勻加速運動區(qū)間U1，...，Uk；

2） 對每個區(qū)間Ui根據目標在該區(qū)間x，y 軸上的運動模型用對應的線性或二次函數對該區(qū)間的運動軌跡進行回歸分析 ，得到目標在該區(qū)間的運行軌跡函數 x（t）和y（t）；

4）對航跡S上各區(qū)間的剩余方差取均值得到各航跡的剩余方差S（x2）和S（y2）；

5）對計算得到的近1000條真實航跡的剩余方差取平均值。

圖15　實時航跡曲線段到直線段的偏差糾正

圖16　實時航跡直線段到曲線段的偏差糾正

表1　435條直線段直線段航跡剩余方差的平均值

表2　386條曲線段曲線段航跡剩余方差的平均值

表3　236條曲線段直線段航跡剩余方差的平均值

表4　215條直線段曲線段航跡剩余方差的平均值

通過實驗結果得到，四種情況下實際航跡與預測航跡的最大剩余方差平均值與最小剩余方差平均值的差值為0.03211，該數量級的誤差，在航跡矢量化處理后顯示在雷達，將會是2至6個像素的誤差，故在允許的誤差范圍之內。故說明該實時航跡矢量化方法是一種正確可靠，有效的實時航跡矢量化方法。

結語

本文提出了一種基于最小二乘法橢圓擬合實時航跡矢量化方法。新方法的可行性已在多次實驗中得到驗證。計算得到實際航跡與預測航跡的最大剩余方差平均值與最小剩余方差平均值的差值誤差，實驗結果理想，驗證了實時航跡預測方法的有效性。新方法不僅可以解決實時航跡的矢量化問題，提供一個可行的實時航跡預測方案，還可以為場間雷達與空中雷達連接方案的制定提供有利參考。

10.3969/j.issn.1001- 8972.2016.18.021