圖切割快速生成扇區(qū)的蟻群算法

葉志堅，王建忠，張召悅，楊群亭，牟龍芳

中國民航大學 空中交通管理學院，天津 300300

管制扇區(qū)是指一塊劃定高度和水平范圍的管制空域。早期飛行活動主要在終端空域，劃分的扇區(qū)形狀接近扇形，而隨著交通越來越復雜，為降低管制員的工作負荷，航路空域也需要劃分成更多的扇區(qū)，但是航路扇區(qū)形狀大都不是扇形。扇區(qū)大小是根據(jù)交通流量、三維空間范圍和管制員的工作量來確定的。劃分扇區(qū)的目的是為了管制員在能力范圍內公平組織交通并保持航空器安全間隔，從而實現(xiàn)管制空域安全、高效和快速地運行。

需求與容量的管理一直是交通運輸領域的核心研究領域。一般來說，容量和基礎設施有密切關系，容量的提升需要較長的時間和巨大的投入。當天氣等因素造成容量下降或者預測的需求會超過容量時，通常的做法是控制需求；但隨著航空領域研究和實踐的發(fā)展，歐美等國的研究人員提出了通過動態(tài)改變管制員管轄范圍來提升空域通行能力的思路并付諸實施，這種做法稱為動態(tài)空域管理[1]。這與國內的開扇、合扇類似，但更加動態(tài)，更加強調空域劃分扇區(qū)適應交通需求，控制需求被視為需求與容量管理的次優(yōu)方案[2]。研究和推廣動態(tài)扇區(qū)劃分運行范式，對改善民航空管運行效率以滿足我國巨大的航空交通需求具有重要意義。

動態(tài)空域劃分研究的第一個熱點是管制員工作負荷或者說是空域的復雜度，管制員工作負荷是管制員主觀感受到的負荷，而空域復雜度是造成管制負荷的根本原因。通過回歸分析法或采用機器學習可以找到兩者之間的相關性聯(lián)系[3]。通過對操作者在工作中的表現(xiàn)[4~6]或借助眼動儀[7]等測試設備，建立任務負荷關聯(lián)模型來測試管制員工作負荷，也得到了國內學者的廣泛研究。但到目前為止，無論是復雜度模型還是管制員工作負荷模型都沒有統(tǒng)一標準。目前歐美采用的容量標準就是最簡單的航空器數(shù)量[8]。這正如奧卡姆剃刀理論在很多領域的運用一樣，如果你有兩個或多個原理，它們都能解釋觀測到的事實，那么你應該使用簡單或可證偽的那個，直到發(fā)現(xiàn)更多的證據(jù)[9]。另外，簡單有效的模型能減少扇區(qū)負荷的計算時間。

動態(tài)空域劃分研究另一個研究熱點是扇區(qū)劃分的模型和算法。由于動態(tài)扇區(qū)劃分核心就是要及時調整空域結構適應交通和環(huán)境的變化，因此算法的計算時間決定了模型能否用于實際運行。第一種是基于圖的模型，圖由邊和頂點構成，邊表示航段，其頂點表示現(xiàn)有航路的交點[10-12]。第二種是基于單元空域增長的模型。首先，將空域劃分為比所研究扇區(qū)更小的單元空域塊，然后對單元空域塊進行重組得到扇區(qū)。單元空域塊有采用正方形的，根據(jù)網(wǎng)絡流計算每個正方形空域塊的負荷，然后再重組[13]；也有用正六邊形作為單元空域塊的，每個空域塊的負荷包括協(xié)調負荷、監(jiān)控負荷、沖突負荷和高度改變工作負荷，單元空域塊的重組采用的是聚類算法，這種方法被用于TAAM（total airspace&airport modeller）系統(tǒng)[14]；有的單元空域塊是不規(guī)則圖形，采用混合整數(shù)規(guī)劃算法重組單元空域塊[15]，有的采用動態(tài)規(guī)劃算法重組單元空域塊[16-17]，有的采用遺傳算法的[18]、采用背包算法[19]以及人工神經(jīng)網(wǎng)絡算法的[20]。但這兩種方法都或多或少需要后處理以使得扇區(qū)形狀滿足特定要求[21]。第三種方法是利用Voronoi圖從上到下分割空域。該方法無需后處理，可一次成型，而且能充分保證最終扇區(qū)的連通性、凸性和壓縮性。這些特點使得這種方法特別有吸引力[22-25]。但這種方法需要的計算量大，計算效率不高是阻礙其有效使用的瓶頸。上述文獻中只有參考文獻[22]報道了計算時間，采用雙層規(guī)劃計算美國沃思堡中心生成16個扇區(qū)大約需要20 min時間。如果用在較為動態(tài)的環(huán)境中，該計算時間還是偏長。本文主要通過模型改進和算法改進來縮短這種自頂向下方法的計算時間。

扇區(qū)劃分涉及的問題非常復雜，有動態(tài)交通流多時段的，也有靜態(tài)交通流單時段的，還可以分為平面2維的，也有3維的。為了減少問題的復雜性，本文研究的是基于某個時段的歷史數(shù)據(jù)，推演今天同一時段的優(yōu)化扇區(qū)結構，是一個靜態(tài)2維的扇區(qū)劃分方法，本文的工作主要聚焦在2維模型的求解上。為未來的研究中結合交通流預測，采用計算機扇區(qū)自動動態(tài)輔助設計奠定基礎。

本文主要創(chuàng)新工作主要體現(xiàn)在以下幾點：采用航空器在所研究時空范圍內4D（4Dimensional）軌跡態(tài)勢判斷為基礎來計算工作負荷，大大縮短了計算工作負荷時間；采用自頂向下切割空域方法，一次成型，減少了預處理時間；通過動態(tài)步長結合蟻群搜索算法自身具有的個體和群體信息交互能力和并行性，生成滿足約束條件的可行解，提高了計算過程的收驗速度。這些工作為快速產(chǎn)生空域劃設預案，減少計算機扇區(qū)輔助設計的時間，服務動態(tài)空域運行模式，提升管制席位適應動態(tài)變化的交通流奠定了基礎。

1 問題描述

1.1 4D軌跡態(tài)勢與任務負荷的關系

為了求解飛機的四維軌跡在某個時間間隔[τi,τi+1)內穿越扇形邊界的準確次數(shù)，假設航空器f進入扇區(qū)時間表達成，航空器f飛出扇區(qū)的時間是。通過以下的例子來推導航空器的協(xié)調次數(shù)和在扇區(qū)中的飛行時間。

假定太原區(qū)調空域只有一個扇區(qū)，假定在τi=9：00時，有七架飛機（AC1、AC2、AC3、AC4、AC5、AC6、AC7），它們的位置在黃色圓圈指示的位置，τi+1=9：15時，它們沿4D軌跡移動到箭頭位置，如圖1所示。

圖1 7架飛機9：00—9：15的歷史軌跡Fig.1 Historical tracjectory of 7 aircrafts during 9：00—9：15

在09：00—09：15間隔時間內，由一架飛機引起的監(jiān)測負荷僅與該飛機在太原空域內（紅色部分）4D軌跡的飛行時間有關，由一架飛機引起的協(xié)調負荷與飛機穿越邊界的次數(shù)有關。有以下幾種情況。

情況1 AC2、AC3和AC4在間隔[τi,τi+1)開始時，τi=9：00，飛機已經(jīng)在扇區(qū)內，，在間隔結束時，飛機在扇區(qū)外，。而且所以當進出扇區(qū)時間和時間段有以下關系時：，則每個航空器與空中交通管制協(xié)調一次，在時間段[τi,τi+1)內航空器在扇區(qū)內的飛行持續(xù)時間為

情況5在屬于在時間段[τi,τi+1)內未曾經(jīng)過扇區(qū)，前面4種情況都是在時間段[τi,τi+1)內經(jīng)過扇區(qū)的。因此，只要判斷飛機軌跡在時間段[τi,τi+1)內的是否經(jīng)過扇區(qū)，可以判斷飛機是否進入扇區(qū)，進入扇區(qū)的航空器分以下4種情況計算協(xié)調次數(shù)和持續(xù)飛行時間，如表1。

表1 情況判斷標準、協(xié)調次數(shù)和持續(xù)飛行時間Table 1 Situation judgment criteria，coordination times and duration of flight

1.2 任務負荷測量

將任務負荷分為監(jiān)視和協(xié)調兩種負荷，監(jiān)視負荷包括了沖突負荷等其他類型負荷，而且和航空器在扇區(qū)內的飛行時間成正比，而協(xié)調負荷取決于協(xié)調一次的時間和協(xié)調次數(shù)。每架次航空器在扇區(qū)內飛行10 min（600 s）需要的監(jiān)控時間和協(xié)調一次所用的時間，都是根據(jù)在太原區(qū)調空域現(xiàn)場測試統(tǒng)計的平均數(shù)。以下描述如何根據(jù)這兩個參數(shù)，以及1.1節(jié)描述的在所研究的時空范圍內航空器4D航跡的態(tài)勢，來計算管制員的任務負荷。

不失一般性，假設空域A包含k個扇區(qū)，s j,j∈{1 ,2,…,k}，這些扇區(qū)預計在時間間隔[τi,τi+1）期間開放。

扇區(qū)s j的任務負荷wl s j包括監(jiān)控和協(xié)調任務負荷，如下式所示：

1.2.1 監(jiān)控任務負荷測量

假設有m架次飛機通過扇區(qū)s j，在間隔[τi,τi+1)內扇區(qū)s j的監(jiān)視任務負荷計算如下：

ls f j是表1計算的飛機f通過扇區(qū)s j的持續(xù)時間（單位：s）。αmot是通過現(xiàn)場測試統(tǒng)計的，每架次航空器在扇區(qū)內飛行10 min（600 s）需要的監(jiān)控時間。

1.2.2 協(xié)調任務負荷測量

協(xié)調任務負荷計算公式如下：

γsfj是根據(jù)表1計算的飛機f在間隔[τi,τi+1)內與s j協(xié)扇區(qū)管制員協(xié)調的次數(shù)。αcod是管制員與飛機協(xié)調一次所用的時間（單位：s）。

1.3 目標函數(shù)

有兩個目標函數(shù)被用于評估解決方案。第一個目標函數(shù)是首要目標，第二個是次要目標。

1.3.1 第一個目標函數(shù)

最小化任務負荷不平衡：

使用以下公式計算任務負荷不平衡程度：

其中，k為需要開的扇區(qū)數(shù)；j為扇區(qū)標號；μ為平均任務負荷為扇區(qū)s j中的總任務負荷。扇區(qū)s j的總任務負荷包括監(jiān)視負荷wlmons j和協(xié)調負荷wlcods j。

1.3.2 第二個目標函數(shù)

最小化總協(xié)調任務負荷：

協(xié)調任務負荷是影響任務總負荷的關鍵因素。由于假設監(jiān)測工作量與飛行器在空域的飛行時間成正比，因此幾乎不可能降低總的監(jiān)測任務負荷。因此，一個好的空域配置方案必須是最小協(xié)調任務負荷的方案。一種扇區(qū)劃分方案的總的協(xié)調負荷是公式（3）表示的各個扇區(qū)協(xié)調負荷的累計。但不能忽略監(jiān)視負荷，因為監(jiān)視負荷涉及到工作負荷平衡，這體現(xiàn)在第一個目標函數(shù)中?？偟谋O(jiān)視任務負荷的計算公式如下：

1.4 約束條件

考慮了三個約束條件。首先，每個扇區(qū)的平均停留時間必須大于一定閾值。其次，必須限制每個扇區(qū)的總負荷量，不能太大，也不能太小。最后，從扇區(qū)邊界到任何關鍵交叉點的距離必須滿足最小距離約束。

1.4.1 平均停留時間約束

美國聯(lián)邦航空局的扇區(qū)容量計算是基于監(jiān)控警報參數(shù)（monitor alert parameter，MAP）計算的[26]。MAP大約是扇區(qū)平均駐留時間（average dwell time，ADT，單位：min）的5/3。這意味著更大的平均駐留時間將有助于提高扇區(qū)容量。平均駐留時間測量如下：

各扇區(qū)的平均駐留時間的約束可以表示為：

adt0是用戶定義的平均駐留時間的最小允許值。為了保證能計算每架次航空器的平均駐留時間，將數(shù)據(jù)提取時間由[τi,τi+1)擴展到[τi-1,τi+2)，保證在時間段[τi,τi+1)空域中出現(xiàn)的飛機，在時間段[τi-1,τi+2)都能完成穿越扇區(qū)。

1.4.2 任務負荷約束

理論上，管制員可以有效使用的最大可接受任務負荷(wlmax)必須小于可用時間。假設在計劃期間[τi,τi+1)，每個管制員的工作負荷的上限是：

為了公平起見，管制員在各扇區(qū)的負荷不應低于平均負荷μ太多。因此，每個扇區(qū)的工作負荷約束如下所示：

α1和α2是用戶定義的介于0和1之間的值，大多介于0.6和0.9之間。μ是平均工作負荷，每個管制員的負荷可以在平均工作負荷附近波動，但是必須小于α1(τi+1-τi)，(τi+1-τi)是工作周期，是一個時間段，總的用時間表示的工作負荷不能大于α1(τi+1-τi)。實際應用場景中會有高工作負荷（大交通流）情況或者低工作負荷情況，低交通流強度情況下，管制員一般不會超負荷，α1和α2可以設置比較小的值，而高強度交通流情況下，α2盡量接近1比較好，α1的設置要根據(jù)實際管制員的水平盡量大一些，強度越大，需要大家共同分擔的任務越多。

1.4.3 最小距離約束

扇區(qū)邊界和任何關鍵交叉點之間的距離不得小于給定距離。此約束確保管制員有足夠的時間解決此節(jié)點上可能發(fā)生的沖突。將由至少兩條路線的交叉口形成的交叉點，且交通流量達到2架次/15 min的定義為臨界交叉口。如果至少存在Q個臨界交點，則從每個臨界交點到任何一個扇區(qū)邊界（多邊形）的水平距離不能小于用戶定義的距離D0。

從臨界交叉點p i:(φi,λi)到任何扇區(qū)s j的邊界的水平距離可以表示為如下：

使用文獻[27]中的方法來求解這個距離。那么，最小距離約束可以表示為：

2 切割空域的蟻群智能算法

如前所述，為了充分利用Voronoi圖能滿足凸性、連通性和壓縮性的特點，來切割空域成扇區(qū)，前面描述了具體的問題和模型，這一章主要研究如何采用蟻群智能算法實現(xiàn)問題的求解。

國外采用Voronoi圖切割空域的文獻中的模型求解基本都用遺傳算法[22-25]，而且計算時間偏長；國內有文獻采用了蟻群算法求解，但從最終解的圖形結果來看，是多個Voronoi圖組合成扇區(qū)，并沒有實現(xiàn)利用Voronoi圖直接切割成具有Voronoi圖顯著特征的完整扇區(qū)，理論上還是屬于元素扇區(qū)重組的范疇[28-30]。本章主要描述解的表達形式和采用蟻群算法的求解過程。

2.1 解的表達形式

用Voronoi圖切割空域，有一個顯著的特點：在空域邊界確定的情況下，最終扇區(qū)的形成主要取決于Voronoi圖的產(chǎn)生點的位置，即Voronoi圖的產(chǎn)生點和最終解有一一對應關系。因此，假設要劃分成k個扇區(qū)解的表達式可以寫成如表2形式。

表2 解的表達形式Table 2 Expression of solution

2.2 蟻群算法的解與蟻群個體的信息映射

為了理解中間解和最終解的關系，以及解與蟻群個體間的關系，下面以父體產(chǎn)生初始種群為例來做說明。不失一般性，假設空域要劃分成k個扇區(qū)，蟻群算法的種群數(shù)量為g。種群的產(chǎn)生過程，解及其個體信息的映射關系如下：

（1）將空域邊界頂點的坐標轉換為平面坐標。

（2）如圖2所示，在所研究空域生成均衡分布的點。

圖2 均衡分布的點及聚類中心Fig.2 Evenly distributed points and clustering centers

（3）在第一代扇區(qū)生成之前，將上一步生成的這些點的K均值聚類中心（圖2中“+”號顯示位置）作為Voronoi圖的生成點。算法開始只有一個父代，尚未產(chǎn)生種群，需要在其周圍鄰域范圍內構造種群，種群中的個體至少有一部分足所在位置一定和聚類中心位置不同。如果種群數(shù)量為g，生成k個扇區(qū)，則一個螞蟻個體有k個足，而且為了保證多樣性g只螞蟻的足應該不完全相同。為了讓個體有差異，需要用圖4所示的聚類中心（圓圈中心）周圍圓圈上的點（聚類中心的鄰域）來替換聚類中心，生成其他個體。聚類中心作為父體，生成種群個體時，為了保證個體的多樣性，需要位置替換；在每次個體更新的時候也要替換，每次替換，相當于個體的位移（傳統(tǒng)蟻群算法，是螞蟻只有一個位置，本文的個體有k個足，k個位置（在2維平面上），相當于k×2維下的一個位置）。個體只有發(fā)生位移才能去探索尋優(yōu)。位置替換以后，個體就可以在尋優(yōu)過程中交換信息，提高并行性能。如果僅僅一個聚類中心，相當于種群只有一個個體。

（4）使用上一步獲得的某個個體的生成點（k只足），即表2所述的中間解，產(chǎn)生扇區(qū)，即表2所述最終解（由每個扇區(qū)邊界頂點構成），例如，如圖3所示，扇區(qū)是利用上一步獲得的Voronoi圖生成點生成的Voronoi圖切割所研究空域產(chǎn)生的。圖3只是一個個體切割空域的例子，種群數(shù)量有g個個體，則有這樣的空域切割方案g個。

圖3 按生成點生成扇區(qū)示意圖Fig.3 Generating sectors schematic by generating points

（5）在圖3的基礎上就可以計算某個個體對應的信息：每個扇區(qū)的工作負荷、個體中最小工作負荷、個體中最大工作負荷、個體中平均工作負荷、不平衡性；航空器平均駐留時間，扇區(qū)邊界到重要點的距離。

解對應的映射信息，是蟻群個體間傳遞信息的基礎，也是解是否可行判斷以及優(yōu)選的基礎。

2.3 蟻群算法的中間解產(chǎn)生機制

不失一般性，假定蟻群個體數(shù)量是g個，ps代表個體上標，則，代表一個中間解。如圖4所示，空域要分成k=6個扇區(qū)，6個扇區(qū)的初始中間解位置在6個Voronoi圖中的圓圈中心位置。

圖4 當前生成點的位置和下一代備用生成點的位置Fig.4 Location of current generation point and next generation standby generation point

定義個體最優(yōu)解為個體在歷代尋優(yōu)過程中的最優(yōu)值，表達式如下：

等式（23）和（24）中的參數(shù)δnc-1被用于限制位置調整量。如果一個扇區(qū)的工作負荷接近個體的平均工作量，則調整幅度會更小。否則，將進行更多的調整。算法初期，扇區(qū)可能不平衡性比較大，在當前位置附近可能沒有接近平衡的點集存在，位置調整量就大些，以便搜索到更平衡的解；在已經(jīng)接近平衡的情況下，最優(yōu)解可能就在附近，需要減少移動步長，增加搜索精度。該參數(shù)計算公式如下：

等式（23）和（24）中的還有三個步長調整參數(shù)，在參考文獻[31]中步長都是固定的（γ1=0.7，γ2和γ3都是1.4）。調整這三個參數(shù)，使它們隨著迭代次數(shù)的增加而減少，如圖5所示。

圖5 γ1、γ2和γ3隨著迭代次數(shù)的變化Fig.5 γ1，γ2 and γ3 changing with number of iterations

三個參數(shù)（γ1，γ2，γ3）隨著迭代次數(shù)的增加而減少可以使后期位置變化更小，扇區(qū)的負荷值變化也更小，精度更高，更容易收斂到負荷平衡點，這三個參數(shù)隨著迭代次數(shù)變化的公式如等式（26）、（27）所示：

其中，nc是當前迭代次數(shù)。

2.4 滿足目標的帕累托前沿選擇機制和最優(yōu)解產(chǎn)生流程

如前所述，主要目標是扇區(qū)負荷平衡，次要目標是總工作負荷。為了得到盡可能平衡的扇區(qū)結果，在優(yōu)選可行解的時候，只選擇扇區(qū)負荷不平衡越來越小的方案，而允許選擇總工作負荷高于前一代的可行解。隨著代數(shù)的增加，允許高于前一代總工作負荷的值越來越小。優(yōu)選條件表達如下：

式中，F(xiàn)2是前一代最佳可行解的總工作負荷，F(xiàn)′2是當前一代某個可行解的的總工作負荷。γ1的表達式見等式（26）。在種群中個體數(shù)量為g，生成k個扇區(qū)的情況下，最優(yōu)解產(chǎn)生流程如圖6所示。

圖6 蟻群算法流程圖Fig.6 Flow chart of ant colony algorithm

2.5 算法終止條件根據(jù)場景設置分析

算法的解首先要滿足約束，但探索過程允許不滿足約束的中間解存在，滿足約束的解才能去參與評優(yōu)，評優(yōu)考慮兩個目標，一個是總工作負荷，一個是扇區(qū)間的不平衡性，兩個目標都是越小越好，而縮小扇區(qū)間的不平衡性是主要目標，因此將“扇區(qū)間的不平衡性”的目標值作為終止條件。如果約束條件越多、越苛刻，比如在重要點分布比較密集的情況下，由于要求扇區(qū)邊界距離重要點必須一定距離，則產(chǎn)生可行解就會變得困難，因此，在這種情況下，要通過Voronoi圖切割空域得到負荷完全平衡的解就會變得非常困難，這時候就要放低終止條件要求，增加“扇區(qū)間的不平衡性”的目標值。

當然也可以用重復迭代多次，目標值沒有改進作為條件，對于要求快速得到結果且對平衡性能滿足允許誤差要求的情況下，可能會消耗更多時間。因此，本文采用“扇區(qū)間的不平衡性”的預設值作為終止條件。

在其他場景下，采用扇區(qū)總工作負荷重復迭代多次無提升，或者扇區(qū)總工作負荷小于某個預設值也是可以的。但扇區(qū)劃分中，總工作負荷最小化，往往扇區(qū)間負荷極度不平衡，除非交通流在時空上均勻分布。因此，很少以此為主要目標，或者在與平衡性構建線性目標函數(shù)時，其權重也會設置相對較小。

因此，算法終止條件的設置要綜合實際應用場景，權衡優(yōu)化目標、得到可行解的難易程度、算法效率、計算時間要求，來根據(jù)具體場景設置終止條件。

3 實例驗證

3.1 實驗設計

實驗目的是為了驗證了模型的正確性和測試算法的計算效率。實驗采用的數(shù)據(jù)是文獻[21]中的山西區(qū)調空域的交通流數(shù)據(jù)。第一個測試是計算質量測試：通過計算模型中的兩個目標（不平衡程度和總工作負荷）來驗證解決方案的質量，并與相關的工作進行了比較。第二個測試是蟻群算法與文獻[21]中MC-CLFV算法重復100次情況下，用箱線圖統(tǒng)計計算質量和計算效率。測試中為了比較計算質量和計算效率，將程序終止條件設置為第一個目標函數(shù)(costimb)達到120，同時第二個目標函數(shù)達到當前最小值。為了測試算法在扇區(qū)數(shù)量增加時的計算效率，松弛了關鍵點到扇區(qū)邊界的平均停留時間約束和最小距離約束。這是由于用于研究的案例空域太小，目前只有四個扇區(qū)，如果保留這兩個約束，就沒有可行的解決方案，或者需要降低性能評估標準。表3列出了算法中使用的參數(shù)。

表3 算法中使用的參數(shù)Table 3 Parameters used in ant colony algorithm

在以下測試中，使用的計算機設備的處理器是英特爾?Xeon?銀牌4215 CPU@2.50 GHz，使用的計算機內存是16 GB RAM。軟件平臺是Matlab 2018a。

3.2 測試1及結果分析

測試1在所有約束條件下，測試扇區(qū)間任務負荷的不平衡程度和所有扇區(qū)的總任務負荷，以評估解決方案的質量。兩個值越低，解的質量越高。圖7是由本文算法生成4個扇區(qū)的扇區(qū)結構圖。

圖7 蟻群算法生成的扇區(qū)Fig.7 Sectors generated by ant colony algorithm

表4 蟻群算法與以往算法結果的比較Table 4 Comparison of results between ant colony algorithm and algorithm in previous work

顯然，扇區(qū)之間的平衡程度和總任務負荷，蟻群算法都優(yōu)于文獻[21]兩種算法的結果和山西區(qū)調管制扇區(qū)當前構型下的值。

圖8是3種算法隨機一次生成4個扇區(qū)的不平衡性隨迭代次數(shù)的變化情況對比。由于有最小距離約束，以及解的鄰域不平衡目標下降，使得優(yōu)化目標值不是隨迭代次數(shù)連續(xù)下降，而是平臺式下降，可以推斷，扇區(qū)劃分的解存在許多鞍點。蟻群算法具有更強的跳出鞍點的能力，MC-CLFV次之，MC-RC跳出鞍點能力最差。從圖8中也可以明顯觀察到蟻群算法的迭代次數(shù)明顯少于另外兩種算法。啟發(fā)式算法的結果有一定隨機性，通過單次的結果比較很難有說服力，因此在測試2中做了重復100次的實驗，以測試改進后的算法的性能統(tǒng)計優(yōu)勢。

圖8 不平衡性隨迭代次數(shù)的變化對比Fig.8 Comparison of variation of imbalance with number of iterations

3.3 測試2及結果分析

在測試2中，測試太原空域劃分為3~8個扇區(qū)時，當終止條件為δ=120時，且松弛了最小平均停留時間約束和關鍵點到扇區(qū)邊界的最小距離約束情況下，蟻群算法和MC-CLFV算法在100次重復計算中的解質量和求解效率統(tǒng)計。

統(tǒng)計結果的箱線圖如圖9所示。粉色箱線圖為MC-CLFV算法運行100次的任務負荷不平衡統(tǒng)計；藍色箱線圖為蟻群算法運行100次的任務負荷不平衡統(tǒng)計。在圖9中的每個框上，凹陷處中心標記表示中位數(shù)值，框的底部和頂部邊緣分別表示下四分位數(shù)Q1（25%）和上四分位數(shù)Q2（75%）。上下邊緣用“－”表示，用“o”符號繪制異常值。

在圖9中，可以看到藍色框（蟻群算法）的中位數(shù)值、Q1和Q2小于粉色框（MC-CLFV算法）。這表明蟻群算法比MC-CLFV算法更有可能得到更好的解。同時，還觀察到兩種算法的三個參數(shù)隨扇區(qū)數(shù)的增加而增加。這意味著，隨著扇區(qū)數(shù)量的增加，解的質量都會下降。

圖9 兩種算法不平衡性的箱線圖對比Fig.9 Comparison of box line diagram with imbalance between two algorithms

兩種算法的總工作負荷統(tǒng)計圖10所示是蟻群算法（藍色）與MC-CLFV算法（粉色）的總工作負荷比較。兩種算法的總工作負荷都隨扇區(qū)數(shù)量增加而增加，這是由于扇區(qū)數(shù)量增加必然增加航空器穿越扇區(qū)邊界的次數(shù)，從而總工作負荷增加。兩種算法都以平衡性為主要目標，因此差異不大。

圖10 兩種算法的總工作負荷比較Fig.10 Comparison of total workload of two algorithms

圖11是蟻群算法（藍色）與MC-CLFV算法（粉色）的計算效率比較圖，很明顯，藍框（蟻群算法）的中位數(shù)、Q1和Q2都小于粉色框（MC-CLFV算法）。這說明，蟻群算法的計算效率高于MC-CLFV算法。也可以看出，隨著扇區(qū)數(shù)量的增加，兩種算法的三個參數(shù)之間的差異越來越大。在扇區(qū)數(shù)量為8的時候，MC-CLFV算法的中位數(shù)值（9 856）約為蟻群算法的中位數(shù)值（1 006）的10倍。另外MC-CLFV算法需要前期元素空域重組，也需要大量的時間。因此，從計算效率上自頂向下地用Voronoi圖切割空域，極大地提高了分扇區(qū)的效率。

圖11 兩種算法的計算時間比較Fig.11 Comparison of computing time between two algorithms

3.4 多層規(guī)劃求解大規(guī)模扇區(qū)劃分的構想

直接用Voronoi圖切割空域，至少生成3個扇區(qū)，因此2分扇區(qū)需要借助其他方法。當需要劃分的扇區(qū)數(shù)量達到9個以上時，采用單層規(guī)劃，一次成型耗費的時間會比較長，采用蟻群算法大約需要1 h，顯然這個計算時間太長。而如果采用雙層規(guī)劃，先把空域分成3個，每一個再分成3個，總共消耗時間相當于4個1分3的時間，大約是4×20=80 s。如果要分成16個扇區(qū)，同樣第一層采用1分4，第二層再采用1分4，采用蟻群算法總消耗時間是5×30=150 s。雖然研究的空域、算法和文獻[22]有很大的不同，但是這篇文章中報道的采用雙層規(guī)劃計算16個扇區(qū)大約需要20 min，遠遠高于推測的計算時間150 s。

如果目標扇區(qū)數(shù)可分解，多級規(guī)劃引理如下：

引理1如果空域被劃分為n個扇區(qū)，并且最多有n個備選方案，且扇區(qū)數(shù)目可以表示為n=ai×b i×ci，i=1,2,…,m，則三層規(guī)劃的計算時間計算公式必是：

m是可能的數(shù)字組合方案數(shù)量。在備選方案中，時間最短的方案必定滿足以下條件：

t a i是一個空域被分割成a i份空域的計算時間。t b i、t c i和t a i的定義一樣。公式中的“|”符號表示條件符號。

對于10個扇區(qū)可以雙層規(guī)劃，先借助其他方法2分，再5分。而11、13、17、19這樣的數(shù)字，可以通過先切割部分扇區(qū)或者增加虛擬扇區(qū)的形式，變成可分和數(shù)來解決。由于篇幅的關系，這個問題留在未來的研究中解決。

4 結論

本文開發(fā)了一種用Voronoi圖直接切割空域形成管制扇區(qū)的方法，模型效率高，能自動滿足凸約束、連通性約束和壓縮性約束，不需要預處理與后處理。基于航空器在所研究時空范圍的4維航跡態(tài)勢來計算工作負荷計算模型能準確估計管制員的工作任務時間，節(jié)約計算時間。本文開發(fā)的蟻群算法比基于蒙特卡洛的MC-CLFV算法探索能力強、求解質量高，而且計算時間大大縮短，計算效率得到顯著提升。

對于劃分9個以上數(shù)量的扇區(qū)，本文提出了采用雙層或者多層規(guī)劃的解決方案，并推斷采用多層規(guī)劃計算時間會比采用單層規(guī)劃顯著減少，對于多層規(guī)劃較大奇數(shù)個扇區(qū)的劃分和設計，需要在未來繼續(xù)研究。