基于粒子濾波算法修正節(jié)點(diǎn)隊(duì)列模型

邱智慧

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

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

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

針對上述問題，本文提出利用節(jié)點(diǎn)隊(duì)列模型得到多條飛機(jī)的下降路徑，并通過粒子濾波算法修正節(jié)點(diǎn)隊(duì)列模型的相關(guān)參數(shù)得到飛機(jī)的延誤時間，通過比對不同路徑的延誤時間得到時間較小的飛行路徑。實(shí)驗(yàn)結(jié)果表明，采用節(jié)點(diǎn)隊(duì)列模型可以準(zhǔn)確得到飛機(jī)的下降路徑，通過粒子濾波算法可以估計(jì)出進(jìn)近飛機(jī)的延誤時間，是一種有效的分析終端區(qū)流量控制的工具。

1 研究與方法

1.1 節(jié)點(diǎn)隊(duì)列模型定義和構(gòu)建

本文主要研究飛機(jī)多條路徑進(jìn)近過程中節(jié)點(diǎn)隊(duì)列模型的構(gòu)建。圖1所示為飛機(jī)到達(dá)終端區(qū)進(jìn)近的空域結(jié)構(gòu)圖[4]，飛機(jī)從起始進(jìn)近定位點(diǎn)(Initial approach fix，IAF)開始經(jīng)由中間進(jìn)近定位點(diǎn)(Intermediate approach fix，IF)，直到最后進(jìn)近定位點(diǎn)(Final approach fix，F(xiàn)AF)實(shí)施著陸，通常飛機(jī)在下降進(jìn)近過程中需要經(jīng)歷終端區(qū)內(nèi)規(guī)定的若干路徑點(diǎn)，當(dāng)經(jīng)過的路徑點(diǎn)間隔與標(biāo)準(zhǔn)進(jìn)近間隔相等時，路徑點(diǎn)與飛機(jī)重疊，飛機(jī)經(jīng)歷的路徑點(diǎn)集合構(gòu)成的下降節(jié)點(diǎn)隊(duì)列(即下降路徑)可等效為節(jié)點(diǎn)隊(duì)列模型。

圖1 終端區(qū)飛機(jī)進(jìn)近的空間結(jié)構(gòu)

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

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

(1)

圖2 終端區(qū)隊(duì)列模型中隊(duì)列狀態(tài)

飛機(jī)的運(yùn)動狀態(tài)由下一節(jié)點(diǎn)的忙閑狀態(tài)決定，Cij=1表示服務(wù)單元空閑，可以接受新服務(wù)，飛機(jī)可以進(jìn)入該服務(wù)單元接受服務(wù)，飛機(jī)狀態(tài)表示為Aij；Cij=0表示服務(wù)單元繁忙，不可以接受新服務(wù)，此時飛機(jī)不可以進(jìn)入服務(wù)單元，需實(shí)施路徑轉(zhuǎn)移，飛機(jī)狀態(tài)表示為Aij′。

根據(jù)上述方法，可以得到多條航跡上飛機(jī)的轉(zhuǎn)移狀態(tài)，當(dāng)一架飛機(jī)到達(dá)時(A表示一架飛機(jī)的狀態(tài))：

所有服務(wù)單元狀態(tài)Cij=1(i=1～M,j=1～N-1)，則該架飛機(jī)可以實(shí)施從服務(wù)單元S11到S1N進(jìn)近著陸。飛機(jī)狀態(tài)表達(dá)式

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

服務(wù)單元狀態(tài)CN=1, C1(N-1)=0,Cij=1(i=1～M,j=1～N-2)，則該架飛機(jī)從服務(wù)單元S11實(shí)施進(jìn)近著陸，在S1(N-2)處進(jìn)行轉(zhuǎn)移。飛機(jī)狀態(tài)表達(dá)式為

服務(wù)單元狀態(tài)CN=1,C1(N-2)=0,Cij=0(i=1～M,j=1～N-3)，則該架飛機(jī)從服務(wù)單元S11實(shí)施進(jìn)近著陸，在S1(N-3)處進(jìn)行轉(zhuǎn)移。飛機(jī)狀態(tài)表達(dá)式為

以此類推，可以得到當(dāng)C11=C12=1、Cij=0(i=1～M,j=3～N-1),則該架飛機(jī)從轉(zhuǎn)移S11到S12便不再繼續(xù)實(shí)施進(jìn)近著陸，需轉(zhuǎn)移到其他路徑點(diǎn)上。飛機(jī)狀態(tài)表達(dá)式為

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

服務(wù)單元狀態(tài)C1j=0,Cij=1(i=2～M,j=1～N-1),則該架飛機(jī)不能實(shí)施從服務(wù)單元S11進(jìn)近著陸，需從其他路徑點(diǎn)直接實(shí)施著陸。飛機(jī)狀態(tài)表達(dá)式為

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

當(dāng)有兩架飛機(jī)達(dá)到終端區(qū)時，仍按照上述方法實(shí)施進(jìn)近(令A(yù)表示第1架飛機(jī)的狀態(tài)、B表示第2架飛機(jī)的狀態(tài))。

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

服務(wù)單元狀態(tài)CN=1，C1N-1=0，其他Cij=1則兩架飛機(jī)從服務(wù)單元S11到S1N-2后便需要轉(zhuǎn)移，從而形成3條路徑。飛機(jī)狀態(tài)表達(dá)式為

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

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

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

上述為兩架飛機(jī)的轉(zhuǎn)移狀態(tài)，以此類推可以得到多架飛機(jī)達(dá)到終端區(qū)的轉(zhuǎn)移狀態(tài)和隊(duì)列結(jié)構(gòu)。

1.2 節(jié)點(diǎn)隊(duì)列模型參數(shù)估計(jì)

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

典型的隊(duì)列分析需要計(jì)算出交通流量參數(shù)，例如延誤時間(等待的時間)，并且給出交通流量指標(biāo)和服務(wù)時間的分布狀態(tài)，其中交通流量指標(biāo)和服務(wù)時間可以通過測量實(shí)時數(shù)據(jù)來得到，而延誤時間是可以預(yù)測得到的，預(yù)測延誤時間是需要通過測量系統(tǒng)來劃分服務(wù)時間和延誤時間。

每次服務(wù)的服務(wù)時間表示為：

(2)

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

時間延遲估計(jì)模型為:

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)為沿著路徑的空速斜率，可以由飛機(jī)中的空速儀表盤來獲得。采樣時間：ΔT=t2-t1。其中tsvc(t2)可表示為：

(7)

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

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

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

(8)

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

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為服務(wù)的出場空速。由此可以得到代入空速加速度后的粒子濾波算法，從而實(shí)現(xiàn)延誤時間的預(yù)測。

2 結(jié)果與分析

2.1 路徑仿真分析

為了構(gòu)建節(jié)點(diǎn)隊(duì)列模型，需在終端區(qū)提取路徑點(diǎn)信息(經(jīng)度、緯度)并表示在坐標(biāo)系中，根據(jù)每個路徑點(diǎn)相關(guān)的緯度和經(jīng)度坐標(biāo),使用公式(13)和(14)，將路徑點(diǎn)的經(jīng)緯度信息轉(zhuǎn)換為笛卡爾坐標(biāo)表示出來。

x=Re·λ節(jié)點(diǎn)-Re·λ目的機(jī)場=Re·Δλ

(13)

y=Re·cosλ目的機(jī)場·(τ節(jié)點(diǎn)-τ目的機(jī)場)

(14)

式中：Re為地球半徑；τ和λ分別表示路徑點(diǎn)的經(jīng)緯度信息。由此，以美國舊金山終端區(qū)路徑信息為例，利用公式(13)和公式(14)就可以將路徑點(diǎn)表示在坐標(biāo)系中。可以得到圖3所示的舊金山終端區(qū)下劃線上的16個路徑點(diǎn)。 從不同方向進(jìn)入終端區(qū)的飛機(jī)，可選擇不同方向上的IAF點(diǎn)實(shí)施著陸，最后達(dá)到FAF點(diǎn)。

圖3 舊金山終端區(qū)著陸路徑點(diǎn)

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

圖4 以3海里劃分的路徑結(jié)構(gòu)

2.2 節(jié)點(diǎn)隊(duì)列模型參數(shù)仿真分析

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

圖5 舊金山機(jī)場飛機(jī)著陸路徑

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

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

由圖7可以看出路徑1延誤時間誤差較大，路徑4的時間均方誤差小，因此綜合考慮路徑4為最優(yōu)路徑。綜合所述，到達(dá)終端區(qū)的飛機(jī)選取路徑4進(jìn)行著陸更優(yōu)。

仍以舊金山機(jī)場節(jié)點(diǎn)隊(duì)列模型所得的5條路徑為基礎(chǔ)，假設(shè)飛機(jī)達(dá)到終端區(qū)后以250 km/h的速度實(shí)施勻減速著陸，得到圖8所示的5條路徑的延誤時間，從圖中可以看出路徑1的延誤時間最小，波動幅度也是最小的，最適合飛機(jī)實(shí)現(xiàn)進(jìn)近著陸。

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

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

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

將擴(kuò)展卡爾曼濾波算法和粒子濾波算法所測得的實(shí)驗(yàn)結(jié)果進(jìn)行對比，得到圖9的實(shí)驗(yàn)結(jié)果對比圖。

由圖9可以看出，飛機(jī)在相同條件的情況下，并以相同的速度和加速度通過同一條路徑時，運(yùn)用粒子濾波算法所測得的延誤時間更短，因此，粒子濾波算法對節(jié)點(diǎn)隊(duì)列模型的參數(shù)估計(jì)性能要比擴(kuò)展卡爾曼濾波算法的更優(yōu)。

圖9 不同算法實(shí)驗(yàn)結(jié)果對比

3 結(jié)論

(1)依據(jù)終端區(qū)的路徑點(diǎn)信息，若路徑點(diǎn)間距離為非標(biāo)準(zhǔn)進(jìn)近間隔，需將路徑點(diǎn)間距離劃分成每3海里一個節(jié)點(diǎn)的形式，再根據(jù)相鄰節(jié)點(diǎn)間只能容納一架飛機(jī)的限定條件，得到飛機(jī)的進(jìn)近路徑。

(2)通過粒子濾波算法預(yù)測模型參數(shù)得到飛機(jī)勻速狀態(tài)下隊(duì)列中的延誤時間分布狀態(tài)，通過對比不同路徑的延誤時間得到最優(yōu)的下降路徑。

(3)通過粒子濾波算法預(yù)測模型參數(shù)得到飛機(jī)變速狀態(tài)下隊(duì)列中的延誤時間分布狀態(tài)，通過對比勻速和變速所測得的延誤時間，可以看出飛機(jī)在變速情況下延誤時間更短，效果更好，在不發(fā)生沖突和飛機(jī)排序的情況下，使得節(jié)點(diǎn)隊(duì)列模型飛機(jī)數(shù)量達(dá)到最優(yōu)化，增大終端區(qū)的吞吐量。

(4)通過將擴(kuò)展卡爾曼濾波算法所得的延誤時間和粒子濾波算法所測的結(jié)果進(jìn)行對比，可以得到粒子濾波算法對于實(shí)現(xiàn)節(jié)點(diǎn)隊(duì)列模型延誤時間的估計(jì)效果更好。