朱承元 晏楠欣

（中國民航大學(xué)空中交通管理學(xué)院 天津300300）

0 引 言

21世紀(jì)的空中交通管制需要應(yīng)用現(xiàn)代計算機科學(xué)與技術(shù)，以適應(yīng)未來更多流動人口的航空旅行需求。對航空旅行吞吐量的最大挑戰(zhàn)之一是惡劣天氣的存在，它導(dǎo)致飛機的延誤、取消和改道，而合理規(guī)劃和優(yōu)化惡劣天氣下的改航路徑是保障空域安全以及提升空域容量的有效途徑之一。

目前，國內(nèi)外已有不少學(xué)者研究了惡劣天氣下的改航問題[1-2]，例如：Joseph Prete等[3-4]研究了基于流量的航線規(guī)劃（flow-based route planner，F(xiàn)BRP）系統(tǒng)的算法，計算出避開時變危險天氣，供飛機安全通過的最短航線。Christine Taylor等[5]利用模擬退火方法確定受天氣影響的航班在運行上可接受的備選航線，實現(xiàn)更有效地使用空域。Jimmy Krozel等[6]將航線搜索算法應(yīng)用于避讓危險天氣的改航問題，得到最優(yōu)改航路徑。Jimmy Krozel等[7]還研究了根據(jù)最新天氣預(yù)報動態(tài)生成編碼起飛航線（CDRs）的方法，利用考慮天氣預(yù)報的算法生成這種天氣避讓路線。李雄等[8]在二維平面利用幾何算法進(jìn)行惡劣天氣下的改航路徑規(guī)劃。趙元棣等[9]利用A*算法構(gòu)建了靜態(tài)和動態(tài)條件下的改航路徑快速規(guī)劃模型，提出了動態(tài)危險天氣下的快速改航規(guī)劃的方法。張兆寧等[10]提出利用多普勒雷達(dá)識別散點狀分布危險天氣區(qū)域，采用平順曲線段和轉(zhuǎn)角曲線段組合的方式，生成終端區(qū)動態(tài)改航路徑。杜實等[11]通過構(gòu)建實時動態(tài)環(huán)境模型，并利用多目標(biāo)粒子群算法（multi-objective particle swarm optimization，MOPSO）對改航路徑進(jìn)行動態(tài)規(guī)劃。

上述研究主要通過建立數(shù)學(xué)模型，進(jìn)行約束限制，達(dá)到減少航班延誤，降低運行成本或生成最優(yōu)改航路徑等目的，其效果良好，但上述模型中存在著變量的關(guān)聯(lián)性復(fù)雜以及數(shù)學(xué)模型求解困難的問題，缺乏對管制員工作負(fù)荷的考慮，未將仿真模型和智能算法有效結(jié)合，不能有效的應(yīng)用于實際工程之中。本文在保證航空器安全運行的前提下，考慮惡劣天氣的動態(tài)變化和管制員工作負(fù)荷，借鑒FBRP方法的思想，考慮以隨機搜索為特征的粒子群算法[12]（particle swarm optimization，PSO）的靈活性以及可在較短時間內(nèi)找到最優(yōu)解的特點，利用全空域與機場模型（total airspace and airport modeller，TAAM）結(jié)合幾何算法和離散粒子群算法（discrete particle swarm optimization，DPSO），進(jìn)行改航路徑的規(guī)劃與優(yōu)化。

1 構(gòu)建改航環(huán)境模型

1.1 劃設(shè)飛行受限區(qū)

1）建立改航平面直角坐標(biāo)系。選取惡劣天氣左下方距離中心地帶較近的某點作為原點，磁北方向設(shè)為y軸正方向，磁北偏東90°設(shè)為x軸正方向。

2）劃設(shè)惡劣天氣區(qū)域。選取降雨量≥13.3 mm/h的區(qū)域，即將在雷達(dá)顯示圖中呈現(xiàn)為黃色、橙色、淡紅色及紅色的區(qū)域作為初始惡劣天氣區(qū)域，利用Matlab軟件和平面點集凸包（Graham）算法[13]，進(jìn)行灰度化處理和凸包劃設(shè)，得到散點狀分布的多個初始飛行受限區(qū)（flight forbidden area，F(xiàn)FA），見圖1。

圖1 劃設(shè)初始飛行受限區(qū)Fig.1 Designation of the initial flight-forbidden area

1.2 預(yù)測飛行受限區(qū)

由于惡劣天氣的移動具有隨機性，因此，參考文獻(xiàn)[14]中灰色模型預(yù)測飛行受限區(qū)，簡化質(zhì)心位置和邊界坐標(biāo)的計算，優(yōu)化外推邊界的計算。目的是簡化計算，使惡劣天氣的呈現(xiàn)更加動態(tài)化，同時縮小最終飛行受限區(qū)的范圍。具體步驟如下。

1）明確觀測時刻Ti(i=0，1，…，s)

2）明確質(zhì)心位置

式中：v x=（xi-xi-1）/T，v y=(yi-yi-1)/T。

3）確定邊界坐標(biāo)

4）確定外推邊界

2 改航路徑的仿真優(yōu)化算法

2.1 仿真優(yōu)化算法的主體步驟

步驟1。依據(jù)初始FFA的范圍，明確優(yōu)化區(qū)域，確定初始方案，進(jìn)行仿真建模。

步驟2。設(shè)航空器到達(dá)改航起始點的時刻為當(dāng)前時刻，根據(jù)前一時刻和當(dāng)前時刻的惡劣天氣區(qū)域，以30 min為時間尺度，預(yù)測未來2 h的FFA"，在TAAM模型中建立相應(yīng)的FFA"。

步驟3。利用幾何算法規(guī)劃多條改航路徑，利用TAAM模型與Matlab軟件進(jìn)行建模和改航路徑選擇的優(yōu)化，從TAAM模型鏈接的MySQL數(shù)據(jù)庫中提取相關(guān)數(shù)據(jù)，計算多航空器對應(yīng)的管制員工作總負(fù)荷和改航總路徑距離。

步驟4。若滿足迭代次數(shù)，則得到優(yōu)化解，即惡劣天氣下的多航空器的優(yōu)化改航路徑；反之，運行改進(jìn)DPSO算法，進(jìn)行新一輪優(yōu)化計算，重復(fù)步驟3～4，直至得到多航空器的優(yōu)化改航路徑。

仿真優(yōu)化算法的主要流程見圖2。

2.2 改航路徑的規(guī)劃

規(guī)劃改航路徑的核心是設(shè)計多條可用的備選航線，并以改航路徑最短為優(yōu)先級，以增加既能避開FFA"又能減少交叉的路徑總數(shù)為原則。

圖2 仿真優(yōu)化算法的流程圖Fig.2 Flow of the simulated optimization algorithm

2.2.1 確定改航點

1）確定初始改航點。確定FFA"與原航線的交匯點（圖3中的ge和g l），選取距離交匯點最近的航路點作為改航的始末點。

2）確定中間轉(zhuǎn)彎點。

步驟1。將交匯點的連線作為分割線（圖3中的l g o g f），取l go g f的中點gm，分別連接對應(yīng)FFA"中距離g m最遠(yuǎn)的點，記為d r和d l。

步驟2。選取rmin=min{d r，d l}，過g m點作垂直于l g o g f的朝向rmin區(qū)域的射線l，以gm為圓心作圓弧，取圓弧和l的交點為中間轉(zhuǎn)彎點g r；然后，考慮偏航角度和航段距離的約束，過g r做平行于原航線的直線l r，再取2個新的改航點g"r，g"r作為中間轉(zhuǎn)彎點，而對應(yīng)的改航始末點依照規(guī)劃改航路徑的原則，內(nèi)移、外推或保持不變。其中，gr(xr，y r)確定見式（6）。

而g"r(x"r，y"r)、g"r(x"r，y"r)可參考文獻(xiàn)[8]中改航點修正的計算。

步驟3。考慮不觸發(fā)TCAS告警的最小距離間隔和航班潛在沖突[16]等因素，對航路進(jìn)行分流，選取rmax=max{d r，d l}，重復(fù)類似步驟2的操作，得到中間轉(zhuǎn)彎點g r1以及新改航點g"r1，g"r2，其坐標(biāo)計算同上。

改航點的確定見圖3。其中灰色區(qū)域表示初始FFA，黑色的虛線攘括的區(qū)域為外推邊界后的FFA"。

圖3 改航點的確定Fig.3 Determination of diverting points

2.2.2 改航約束

1）最大偏航角度限制。最大偏航角度不大于90°。

2）改航路徑有效性限制。改航路徑與FFA"相交，則改航路徑無效。反之，改航路徑有效。

3）中間轉(zhuǎn)彎點數(shù)量為定值限制。當(dāng)FFA"被航線一分為二時，有改航路徑的一側(cè)或兩側(cè)，其對應(yīng)的中間轉(zhuǎn)彎點數(shù)量為3。

4）最小航段距離限制 航空器為完成相應(yīng)的轉(zhuǎn)彎，飛行距離必須大于3.7 km；若要順利完成2次轉(zhuǎn)彎，則飛行距離不小于7.4 km，即圖3中

5）流量控制固定點的數(shù)量限制 流量控制固定點定義為飛行受限區(qū)對應(yīng)的改航始（末）點。數(shù)量限制是指FFA"對應(yīng)流量控制固定點的數(shù)值只能為2，3或4，即2≤改航始（末）點的數(shù)量≤4。

6）改航路徑與FFA"最小側(cè)向間隔限制 改航路徑與FFA"保持至少14 km（航路兼最小空中走廊寬度的1/2）的側(cè)向間隔。

其中，最大偏航角度限制和改航路徑有效性限制的約束，可通過改航始（末）點的內(nèi)移、外推修正。

2.3 改航路徑優(yōu)化的數(shù)學(xué)模型

以多架航空器的改航總路徑（L）最短和改航路徑對應(yīng)扇區(qū)s管制員工作總負(fù)荷最小為目標(biāo)，建立多目標(biāo)數(shù)學(xué)模型，見式（7）～（9）。

式中：αi(i=1，2，…，5）為各類管制員工作負(fù)荷的權(quán)重系數(shù)。

針對上述數(shù)學(xué)模型，有如下約束條件。

1）改航時間限制。根據(jù)惡劣天氣的影響范圍，預(yù)估改航所需要總時間T。設(shè)D sf為航空器對應(yīng)改航始末點之間的距離；vc為航空器的巡航速度；α為繞飛FFA"所需增加航程的最大百分比。

2）航空器與惡劣天氣安全間隔限制。即航空器與惡劣天氣的間隔至少大于8 km的最小空中走廊寬度。

3）流量控制固定點的流量限制。流量限制為在滿足我國民航區(qū)域雷達(dá)管制間隔（10 km）的前提下，3 min內(nèi)過固定點的飛機流≤2架。

4）低復(fù)雜度路徑優(yōu)先選擇限制。此處的復(fù)雜度定義為規(guī)劃改航路徑的航路點數(shù)量。數(shù)量越少對應(yīng)復(fù)雜度越低，反之，越高。

2.4 改進(jìn)DPSO算法

本文針對基于幾何算法規(guī)劃的改航路徑和原航線中部分固定航段的組合進(jìn)行優(yōu)化，該問題屬于離散問題，也屬于NP問題。因此，借鑒相關(guān)文獻(xiàn)[18]求解旅行商的DPSO算法的設(shè)計思想，重新定義了DPSO算法中狀態(tài)表示和運算規(guī)則。以貴陽區(qū)域管制區(qū)為例，多航空器改航路徑的見圖4。

圖4 多航空器改航路徑示意圖Fig.4 Diverting routes of multi aircrafts

2.4.1 運算規(guī)則的定義

1）粒子位置的說明。令第i(i=1，2，…，N)個粒子的位置向量為

式中：Xi1，Xi2，…，Xiq，…，XiW為W個航班的航班走向。第i個粒子的第q個航班的航班走向的可表示為

2）粒子速度和位置的更新公式。第i個粒子在DPSO算法中速度和位置的向量公式為

式中：V i為粒子i的速度向量，i=1，2，…，m；k為當(dāng)前迭代步數(shù)；r1，r2為學(xué)習(xí)因子，取[ ]0,1之間的隨機數(shù)；pb為粒子i的個體最優(yōu)位置；g b為粒子i的群體最優(yōu)位置。

第i個粒子的第q個航班，在第j個飛行受限區(qū)中可用改航路徑的速度和位置的向量公式，見式（16）～（17）。

2.4.2 算法的設(shè)計

1）確定粒子編碼。為簡化算法的運算，對改航路徑中的關(guān)鍵改航點進(jìn)行編碼，對式（13）中的中間轉(zhuǎn)彎點進(jìn)行編碼，則第i個粒子的第q個航班改航路徑編碼為以圖4中，1號飛行受限區(qū)的第i個粒子的第q個航班的編碼為例，k的取值為1，2，3和4，其具體編碼如下。

2）確定初始粒子種群。將FFA"下的幾何算法規(guī)劃的改航路徑的集合內(nèi)可隨機選擇的航班走向作為改進(jìn)DPSO算法的初始粒子種群。

3）確定適應(yīng)度函數(shù)。令fit(x)=minf。其中，f為改航路徑優(yōu)化的多目標(biāo)函數(shù)。

3 實例仿真分析

我國西南地區(qū)為雷暴多發(fā)地帶，由中國氣象數(shù)據(jù)網(wǎng)（http：//data.cma.cn）可知，2020年7月18日18：00—19日10：00，貴州境內(nèi)持續(xù)有惡劣天氣，選取貴州區(qū)域管制區(qū)作為改航對象。首先，根據(jù)中國民航國內(nèi)航空資料匯編，得到關(guān)于貴陽管制區(qū)（ZUGY）內(nèi)航路航線、航路點和扇區(qū)等相關(guān)信息；其次，結(jié)合從貴陽管制單位獲得的相關(guān)數(shù)據(jù)，選取受惡劣天氣影響的航路（W180，W181，W182，A581和H24）和經(jīng)過相關(guān)航路航段的航空器，進(jìn)行改航的規(guī)劃和優(yōu)化。

其中，經(jīng)過ZUGY區(qū)域并以W180，W181，W182，A581和H24航路為主的航班走向有百余次，參考貴陽管制區(qū)2019年高峰日高峰時段的實際運行數(shù)據(jù)，選取18：30—20：30的飛越、起飛和落地貴陽區(qū)域管制區(qū)的11個航班走向，45個航班，進(jìn)行初始建模。

3.1 相關(guān)系數(shù)的確定

1）參考文獻(xiàn)[15]，確定環(huán)境模型中FFA"的安全裕度：σmin=5 km，σmax=25 km。

2）參考文獻(xiàn)[11]，繞飛FFA"所需增加航程的最大百分比α取20%；

3）確定DPSO算法中的各類參數(shù)。參考臨時航線的運行特點，固定值B取0.6；迭代次數(shù)取100，粒子群種群數(shù)取20；使用基于群組AHP與交叉熵的組合賦權(quán)法[17]，并進(jìn)行規(guī)范化處理，最終確定的多目標(biāo)函數(shù)權(quán)重為：w1=0.584，w2=0.023。

4）TAAM模型中管制員工作負(fù)荷的各類負(fù)荷權(quán)重αi(i=1，2，…，5)取默認(rèn)設(shè)置，依次為0.2，0.1，0.1，0.3，0.3。

3.2 仿真數(shù)據(jù)分析

參考文獻(xiàn)[14]，以像素為單位建立惡劣天氣的坐標(biāo)系，明確惡劣天氣的坐標(biāo)和范圍，計算精確度Aa和冗余度Ab。其中，以18：24—18：36（T0，T1）內(nèi)H24航路惡劣天氣的雷達(dá)實況為例，預(yù)測18：36—18：54即Ti(i=2，3，4)的飛行受限區(qū)，則18：24—18：54的飛行受限區(qū)見圖5。此外，計算得到A a=76.43%，Ab=2.68%。

圖5 惡劣天氣對應(yīng)的飛行受限區(qū)圖示Fig.5 Flight-forbidden area corresponding to severe weather

1）構(gòu)建改航環(huán)境。采用Graham算法對零星并且相近的惡劣天氣進(jìn)行凸殼劃設(shè)，并基于滾動時間窗（天氣預(yù)報的雷達(dá)更新周期為6 min）進(jìn)行惡劣天氣的預(yù)測以及合并，得到貴陽區(qū)域管制區(qū)內(nèi)的散點狀分布的FFA"，見圖6（a）。

2）生成改航路徑。根據(jù)環(huán)境模型，考慮實際運行和相關(guān)約束，利用幾何算法規(guī)劃多條改航路徑，并設(shè)置空域規(guī)則以貼近實際的航班運行，改航路徑在TAAM模型中的情況，見圖6（b）。

圖6 改航環(huán)境和改航路徑的圖示Fig.6 Diagrammatic representation of the redirected environment and diverting routes

3）管制工作負(fù)荷和改航路徑的匯總。惡劣天氣主要影響的是AR01，AR03和AR04扇，改航數(shù)據(jù)匯總見表1，管制員小時負(fù)荷數(shù)據(jù)見圖7。其中Origin代表原始（改航前）方案，PSO代表PSO算法對應(yīng)的最優(yōu)方案，DPSO代表改進(jìn)DPSO算法對應(yīng)的最優(yōu)方案。而改航環(huán)境下算法適應(yīng)度的變化見圖8，其中改航距離通過規(guī)范化處理，量級上與管制員工作負(fù)荷基本保持一致，故適應(yīng)度值的單位取為當(dāng)量架次。

圖7 管制員工作小時負(fù)荷圖示Fig.7 Graphical representation of workhour load of the controller

表1 改航數(shù)據(jù)匯總Tab.1 Summary of diverting data

圖8 改航環(huán)境下算法適應(yīng)度的變化Fig.8 Variation of algorithm"s adaptation in diversion

由圖8可知，迭代50次左右，改進(jìn)DPSO和PSO算法適應(yīng)度值均趨于平穩(wěn)，改進(jìn)DPSO算法的適應(yīng)度值更低，并且其最優(yōu)方案的管制員工作負(fù)荷和改航距離較原方案分別增加了9.86%和18.66%，而對比PSO算法最優(yōu)方案的管制員工作負(fù)荷和改航距離則分別降低了7.52%和4.48%，改航效果較貼合實際。以16：48：00從杭州蕭山機場起飛，18：54：34落地畢節(jié)飛雄機場（貴陽管制區(qū)內(nèi)的小機場）的CDC8823航班數(shù)據(jù)為例，其機型為A320，巡航速度（地速）為814.88 km/h，具體改航數(shù)據(jù)見表2。其中A1，A2，A3，B1，B2和B3均為改航點，表中距離為與跑道末端的距離，n mile。

由表2可知，CDC8823航班的改航滿足式（11）改航時間限制的約束條件。

3.3 改航算法的對比和驗證

為驗證惡劣天氣下多航空器改航路徑的仿真優(yōu)化算法的有效性、科學(xué)性，對比基于惡劣天氣的多目標(biāo)粒子群算法（MOPSO）和非支配排序遺傳算法（NSGA-II）的改航路徑規(guī)劃。在2種對比算法中對應(yīng)構(gòu)建30×30的柵格環(huán)境，每個柵格設(shè)置為13 km×13 km[11]，并以ZUGY區(qū)域中受惡劣天氣影響的H24航路中的航段（P159-ZHJ-P293-XONID-UBDID-MASRO）為例，此航路段長183 n mile（338.916 km），其與FFA"在柵格環(huán)境中的情況，見圖9（a）。

其一，采用MOPSO算法，設(shè)置相關(guān)參數(shù)[11]；其二，采用NSGA-II算法，設(shè)置適當(dāng)參數(shù)[19]，其中種群數(shù)取40。2種算法的改航路徑結(jié)果見圖9。

其中，MOPSO算法的改航路徑長為350.032 km；NSGA-II算法的改航路徑長為347.849 km；以CDC8823航班數(shù)據(jù)為例，考慮W182航路上的FFA"，在XONID-UBDID有2次改航，采用改進(jìn)DPSO算法的仿真優(yōu)化算法的改航路徑長為386.698 km（含2次改航增加的34.077 km），而采用PSO算法的仿真優(yōu)化算法的改航路徑長為394.291 km（含2次改航增加的40.188 km）。以A320巡航速度（地速）814.88 km/h為例，上述4種方法均滿足偏航角度限制和改航時間限制的約束。

表2 CDC8823航班的改航數(shù)據(jù)Tab.2 Date of the diversion of Flight CDC8823

圖9 MOPSO算法和NSGA-II算法的改航路徑結(jié)果Fig.9 Results of the MOPSOalgorithm and the NSGA-II algorithm for diversion

此外，雖然MOPSO算法與NSGA-II算法和本文方法相比，其改航路徑距離均略小，但是其算法的運算時間較長，未考慮管制員工作負(fù)荷，當(dāng)出現(xiàn)空中交通擁堵時，可能存在改航路徑與其他航班飛行路徑的潛在沖突，且在Matlab程序中難以實現(xiàn)管制員的工作模擬和實時仿真。

而本文方法規(guī)劃了多條改航路徑，考慮了管制員工作負(fù)荷，利用了改進(jìn)DPSO算法和TAAM軟件的仿真優(yōu)化算法，可在縮短改航路徑距離的同時減少管制員工作負(fù)荷，較之其他算法，更加貼近實際管制工作，更具有工程上的可行性。

4 結(jié)束語

本篇針對典型區(qū)域管制區(qū)進(jìn)行改航規(guī)劃和優(yōu)化，從航班安全運行和管制員有效調(diào)配的角度出發(fā)，結(jié)合約束條件，改進(jìn)DPSO算法中粒子速度和位置的運算規(guī)則，使其適合于改航路徑的優(yōu)化與求解，通過Matlab軟件和TAAM模型實現(xiàn)算法的多次迭代，快速得到多個航空器的優(yōu)化改航路徑。結(jié)果證明，本文研究的仿真優(yōu)化算法，可達(dá)到航班分流，減少管制員工作負(fù)荷和縮短改航路徑的目的。

未來工作包括優(yōu)化環(huán)境模型，建立更精細(xì)的動態(tài)環(huán)境模型，改進(jìn)惡劣天氣的度量，以便更好地說明哪些天氣類型（形狀、大小、嚴(yán)重程度）會對改航產(chǎn)生不利影響，以及考慮改航路徑的三維建模。