基于粒子濾波算法修正節點隊列模型

邱智慧

(中國民航大學電子信息與自動化學院，天津 300300)

民航空中交通運行區域可分為航路、終端區及機場等三個部分，其中終端區是航路起始、終止的必經區域，也是整個空域系統運行過程中的瓶頸區域[1]。由于終端區的空域容量有限，進近飛機數量多，因此如何有效地實現飛機進近排序是終端區交通流量控制的重要任務。

針對終端區飛機進近排序問題，國內外學者提出了很多創新性方案。來自NASA的艾姆斯研究中心的Gregory C. Carr等人為了不改變飛機順序實現終端區飛機進近排序，提出了時間提前(TimeAdvance，TA)算法[2],通過提前首架飛機的到達時間來解決后續飛機的延誤時間。該算法要求在飛機隊列中必須有兩架飛機滿足最小時間間隔約束條件，若無滿足條件的飛機，該算法將不能適用。我國的劉洪等人提出了一種針對動態飛機流的優化調度算法[3]。該算法結合了各種空管限制條件,引入飛機的優先級,并通過實驗驗證了該算法的可行性，但是實現管制限定條件的獲取上具有一定難度。

針對上述問題，本文提出利用節點隊列模型得到多條飛機的下降路徑，并通過粒子濾波算法修正節點隊列模型的相關參數得到飛機的延誤時間，通過比對不同路徑的延誤時間得到時間較小的飛行路徑。實驗結果表明，采用節點隊列模型可以準確得到飛機的下降路徑，通過粒子濾波算法可以估計出進近飛機的延誤時間，是一種有效的分析終端區流量控制的工具。

1 研究與方法

1.1 節點隊列模型定義和構建

本文主要研究飛機多條路徑進近過程中節點隊列模型的構建。圖1所示為飛機到達終端區進近的空域結構圖[4]，飛機從起始進近定位點(Initial approach fix，IAF)開始經由中間進近定位點(Intermediate approach fix，IF)，直到最后進近定位點(Final approach fix，FAF)實施著陸，通常飛機在下降進近過程中需要經歷終端區內規定的若干路徑點，當經過的路徑點間隔與標準進近間隔相等時，路徑點與飛機重疊，飛機經歷的路徑點集合構成的下降節點隊列(即下降路徑)可等效為節點隊列模型。

圖1 終端區飛機進近的空間結構

飛機在終端區運行時，路徑點間距離通常是非等間隔且具有多條路徑，因此利用公式(1)計算得出((xij,yij)表示路徑點Sij的坐標)路徑點間間隔，再將路徑點距離離散分割成標準進近間隔進行分析。若飛機在下一路徑點產生沖突，則應在當前路徑點進行轉移，且應保障在距離FAF前5海里處完成飛機航路的轉移。如圖2所示，為保障飛機安全著陸，應在服務單元SN之前調整飛機狀態使之對正下降著陸航段。巡航飛機通過不同航路及高度層路徑到達目的機場終端區，且均通過服務單元S11開始實施進近著陸，經過S12至SN點到著陸形成一條下降路徑，當一條路徑無法滿足飛機進行著陸時，就會轉移到Sij,i=2～M路徑點上，形成M條路徑。其中Sij表示服務單元(i=1～M,j=1～N)，M表示飛機下降所形成的路徑數目，N為終端區飛機制定飛行計劃時所經歷的路徑點中除FAF外的其他路徑點數目，相鄰服務單元間的距離為3海里。

(i=1～M,j=1～N-1)

(1)

圖2 終端區隊列模型中隊列狀態

飛機的運動狀態由下一節點的忙閑狀態決定，Cij=1表示服務單元空閑，可以接受新服務，飛機可以進入該服務單元接受服務，飛機狀態表示為Aij；Cij=0表示服務單元繁忙，不可以接受新服務，此時飛機不可以進入服務單元，需實施路徑轉移，飛機狀態表示為Aij′。

根據上述方法，可以得到多條航跡上飛機的轉移狀態，當一架飛機到達時(A表示一架飛機的狀態)：

所有服務單元狀態Cij=1(i=1～M,j=1～N-1)，則該架飛機可以實施從服務單元S11到S1N進近著陸。飛機狀態表達式

[AS11→S12，AS12→S13，…，AS1(N-1)→SN,ASN→SFAF]

服務單元狀態CN=1, C1(N-1)=0,Cij=1(i=1～M,j=1～N-2)，則該架飛機從服務單元S11實施進近著陸，在S1(N-2)處進行轉移。飛機狀態表達式為

服務單元狀態CN=1,C1(N-2)=0,Cij=0(i=1～M,j=1～N-3)，則該架飛機從服務單元S11實施進近著陸，在S1(N-3)處進行轉移。飛機狀態表達式為

以此類推，可以得到當C11=C12=1、Cij=0(i=1～M,j=3～N-1),則該架飛機從轉移S11到S12便不再繼續實施進近著陸，需轉移到其他路徑點上。飛機狀態表達式為

C11=1，C1j=0，Cij=1(i=2～M,j=1～N-1)，則該架飛機進入S11后便不能再實施進近著陸，需進行轉移到其他路徑點上。飛機狀態表達式為

服務單元狀態C1j=0,Cij=1(i=2～M,j=1～N-1),則該架飛機不能實施從服務單元S11進近著陸，需從其他路徑點直接實施著陸。飛機狀態表達式為

[ASi1→Si2,ASi2→Si3,…,ASi(N-1)→SN,ASN→SFAF]

當有兩架飛機達到終端區時，仍按照上述方法實施進近(令A表示第1架飛機的狀態、B表示第2架飛機的狀態)。

所有服務單元狀態Cij=1(i=1～M,j=1～N-1)，且CN=1則兩飛機可以實施從服務單元S11到S1N進近著陸，但是為了不產生沖突，第2架飛機需在SN-1處轉移，形成兩條航跡。飛機狀態表達式為

服務單元狀態CN=1，C1N-1=0，其他Cij=1則兩架飛機從服務單元S11到S1N-2后便需要轉移，從而形成3條路徑。飛機狀態表達式為

以此類推可得，當CN=1，C11=C12=1、Cij=0(i=1～M,j=3-N-1)，則兩架飛機從轉移S11到S12便不再繼續實施進近著陸，需轉移到其他路徑點上，從而形成3條路徑。飛機狀態表達式為

服務單元狀態CN=1,C11=1,C1j=0,Cij=1(i=2～M,j=1～N-1)，則兩架飛機進入S11后便不能再實施進近著陸，需進行轉移到其他路徑點上，從而形成3條路徑。飛機狀態表達式為

所有服務單元狀態CN=1,C1j=0,Cij=1(i=2～M,j=1～N-1)，則兩架飛機不能實施從服務單元S11進近著陸，需從其他路徑點直接實施著陸。飛機狀態表達式為：

上述為兩架飛機的轉移狀態，以此類推可以得到多架飛機達到終端區的轉移狀態和隊列結構。

1.2 節點隊列模型參數估計

由于節點隊列模型中參數是非線性、非高斯分布的[5-6]，飛機的到達時間分布和服務時間分布都是大于零，其次，這些分布不是高斯線性分布，通常服從Poisson、Erlang分布描述。因此需要考慮更先進的技術，選擇可適用于節點隊列模型的算法，所以本文采用粒子濾波算法實現對模型參數的估計[7-8]，粒子濾波算法和擴展卡爾曼相比可以不依賴于任何局部的線性化技術，也不使用任何非線性函數近似逼近方法，它是通過有限數量的狀態樣本(粒子云)來近似后驗概率，每個樣本可大致對應于狀態空間中的區域,其精度可以逼近最優估計。

典型的隊列分析需要計算出交通流量參數，例如延誤時間(等待的時間)，并且給出交通流量指標和服務時間的分布狀態，其中交通流量指標和服務時間可以通過測量實時數據來得到，而延誤時間是可以預測得到的，預測延誤時間是需要通過測量系統來劃分服務時間和延誤時間。

每次服務的服務時間表示為：

(2)

式中：s代表服務的長度，在本文中為由若干3海里的區間組成的路徑點之間的距離；V(s)代表沿著路徑長度的空速參數；L代表從起始進近定位點到最終進近定位的直徑距離。

時間延遲估計模型為:

a(t2)=a(t1)+wap(t1)*ΔT

(3)

x1(t2)=x1(t1)+x2(t2)*ΔT+1/2*ΔT*wtp(t1)

(4)

x2(t2)=x2(t1)+ΔTwtp(t1)

(5)

y1(t2)=x1(t2)+tSVC(t2)+ws(t2)

(6)

式中：x1、x2為延誤時間，wap(t1)是變化的空速斜率的噪聲；wtp(t1)是延誤時間的噪聲；ws(t2)是傳感器噪聲；ɑ(t)為沿著路徑的空速斜率，可以由飛機中的空速儀表盤來獲得。采樣時間：ΔT=t2-t1。其中tsvc(t2)可表示為：

(7)

式中：Ventry為飛機在終端區的出場空速，在上述粒子濾波算法中，飛機需保持勻速實施進近，但是飛機在終端區運行時，速度是勻減速變化的，因此需對上式中tsvc進行修正，代入飛機空速變化參數，進行分析。

飛機在終端區實施進近著陸時，飛機的空速不是一個定值，而是不斷變化的，因此可以引入空速變化斜率來修正粒子濾波方程的測量方程，從而使得觀測結果更加準確。假設空速的變化是沿著路徑線性變化的：V(s)=as+b,上式的積分可以表示為

tsvc=1/a*ln(a/b*L+1)

(8)

式中：ɑ為沿著路徑的空速斜率；b為進入終端區時的入場速度。由此可得修正后的粒子濾波算法方程為：

x1(t2)=x1(t1)+x2(t2)*ΔT+1/2*ΔT*wtp(t1)

(9)

x2(t2)=x2(t1)+ΔTwtp(t1)

(10)

y1(t2)=x1(t2)+tsvc(t2)+ws(t2)

(11)

(12)

式中：Ventry為服務的出場空速。由此可以得到代入空速加速度后的粒子濾波算法，從而實現延誤時間的預測。

2 結果與分析

2.1 路徑仿真分析

為了構建節點隊列模型，需在終端區提取路徑點信息(經度、緯度)并表示在坐標系中，根據每個路徑點相關的緯度和經度坐標,使用公式(13)和(14)，將路徑點的經緯度信息轉換為笛卡爾坐標表示出來。

x=Re·λ節點-Re·λ目的機場=Re·Δλ

(13)

y=Re·cosλ目的機場·(τ節點-τ目的機場)

(14)

式中：Re為地球半徑；τ和λ分別表示路徑點的經緯度信息。由此，以美國舊金山終端區路徑信息為例，利用公式(13)和公式(14)就可以將路徑點表示在坐標系中?？梢缘玫綀D3所示的舊金山終端區下劃線上的16個路徑點。 從不同方向進入終端區的飛機，可選擇不同方向上的IAF點實施著陸，最后達到FAF點。

圖3 舊金山終端區著陸路徑點

對于終端區路徑點間距離為非標準進近的間距時，需要將相鄰路徑點間距離劃分成標準進近距離，得到飛機的排隊隊列，由此得到圖4中的飛機節點隊列，且相鄰節點間僅能容納一架飛機。得到圖5所示的舊金山飛機著陸路徑信息，由該圖可以看出共形成5條下降路徑。接下來在此基礎上利用粒子濾波算法對延誤時間進行預測分析。

圖4 以3海里劃分的路徑結構

2.2 節點隊列模型參數仿真分析

以圖5中舊金山機場節點隊列模型所得的5條路徑為基礎，假設飛機達到終端區后保持230 km/h的速度勻速著陸，通過粒子濾波算法在采樣100次的條件下，得到5條路徑的延遲時間，如圖6所示，從圖中可以看出路徑3和路徑4的延遲時間相對其他路徑較小，波動幅度也是最小的，最適合飛機實現進近著陸。

圖5 舊金山機場飛機著陸路徑

圖6 粒子濾波算法得到5條路徑延誤時間

將5條路徑所得的仿真時間與真實時間進行比對，得到仿真時間的誤差，圖7所示為典型路徑的誤差時間。

由圖7可以看出路徑1延誤時間誤差較大，路徑4的時間均方誤差小，因此綜合考慮路徑4為最優路徑。綜合所述，到達終端區的飛機選取路徑4進行著陸更優。

仍以舊金山機場節點隊列模型所得的5條路徑為基礎，假設飛機達到終端區后以250 km/h的速度實施勻減速著陸，得到圖8所示的5條路徑的延誤時間，從圖中可以看出路徑1的延誤時間最小，波動幅度也是最小的，最適合飛機實現進近著陸。

通過以圖6和圖8做對比可以看出，引入空速斜率后的粒子濾波算法準確值更高，極大地縮短了延誤時間，若到達機場飛機采用該算法實現進近著陸，將會極大的減少終端區的延誤現象，使得終端區空域充分利用。

圖7 2條典型路徑的延誤時間均方誤差

圖8 引入空速斜率粒子濾波的延誤時間

將擴展卡爾曼濾波算法和粒子濾波算法所測得的實驗結果進行對比，得到圖9的實驗結果對比圖。

由圖9可以看出，飛機在相同條件的情況下，并以相同的速度和加速度通過同一條路徑時，運用粒子濾波算法所測得的延誤時間更短，因此，粒子濾波算法對節點隊列模型的參數估計性能要比擴展卡爾曼濾波算法的更優。

圖9 不同算法實驗結果對比

3 結論

(1)依據終端區的路徑點信息，若路徑點間距離為非標準進近間隔，需將路徑點間距離劃分成每3海里一個節點的形式，再根據相鄰節點間只能容納一架飛機的限定條件，得到飛機的進近路徑。

(2)通過粒子濾波算法預測模型參數得到飛機勻速狀態下隊列中的延誤時間分布狀態，通過對比不同路徑的延誤時間得到最優的下降路徑。

(3)通過粒子濾波算法預測模型參數得到飛機變速狀態下隊列中的延誤時間分布狀態，通過對比勻速和變速所測得的延誤時間，可以看出飛機在變速情況下延誤時間更短，效果更好，在不發生沖突和飛機排序的情況下，使得節點隊列模型飛機數量達到最優化，增大終端區的吞吐量。

(4)通過將擴展卡爾曼濾波算法所得的延誤時間和粒子濾波算法所測的結果進行對比，可以得到粒子濾波算法對于實現節點隊列模型延誤時間的估計效果更好。