徐燦， 田勇*， 牛科新， 鞏文健， 李桂芳

(1.南京航空航天大學(xué)民航學(xué)院， 南京 211106； 2.中國民用航空華東地區(qū)空中交通管理局， 上海 200335)

隨著民用航空的穩(wěn)步發(fā)展[1]，現(xiàn)有的空域容量無法滿足快速增長的交通需求[2-3]。終端區(qū)作為航空器飛行過程中的重要區(qū)域，運(yùn)行環(huán)境復(fù)雜，空域容量時(shí)常逼近閾值。目前大多終端區(qū)扇區(qū)是參照國際民航組織的相關(guān)文件，結(jié)合主觀經(jīng)驗(yàn)設(shè)計(jì)的。這顯然不能滿足日益增長的交通需求，易導(dǎo)致管制負(fù)荷不均勻，最終限制空域容量提升[4]。

對此，中外學(xué)者一直積極探索科學(xué)的扇區(qū)劃設(shè)方法。Kumar[5]以減少飛行碰撞為目標(biāo)劃設(shè)扇區(qū)。Sergeeva等[6]提出了一種啟發(fā)式算法優(yōu)化管制扇區(qū)。Zenlinski[7]總結(jié)了Voronoi圖、混合整數(shù)規(guī)劃和航跡聚類3種經(jīng)典扇區(qū)劃設(shè)方法的優(yōu)缺點(diǎn)。Yin等[8-9]對新加坡高空空域進(jìn)行了多目標(biāo)扇區(qū)劃設(shè)。Ghorpade[10]提出了一種基于加權(quán)圖的空域配置模型。Oktal等[11]對建立的多目標(biāo)扇區(qū)劃設(shè)模型進(jìn)行了標(biāo)量化。

近年來，學(xué)者們開始考慮三維扇區(qū)劃設(shè)。Tang等[12-13]提出了改進(jìn)的智能體模型，較好地滿足了相關(guān)約束。Min[14]提出了一種減少過載扇區(qū)管制負(fù)荷的三維扇區(qū)劃設(shè)方法。Chen等[15]通過改進(jìn)遺傳算法進(jìn)行了三維扇區(qū)劃設(shè)。王超等[16]使用二叉空間分割法劃設(shè)扇區(qū)。Wong等[17]提出了考慮下一小時(shí)流量的動態(tài)扇區(qū)配置方法。Wong等[18]提出了一種滾動優(yōu)化方法，根據(jù)流量的不同實(shí)現(xiàn)扇區(qū)的漸變。

相比于高空區(qū)域，終端區(qū)扇區(qū)劃設(shè)的研究較少。Wei等[19]對美國典型終端區(qū)進(jìn)行研究，先在垂直方向上進(jìn)行分割，而后對分割的高度層分別采用混合整數(shù)規(guī)劃方法劃設(shè)扇區(qū)。Granberg等[20]提出了新的終端區(qū)扇區(qū)劃分方法，能夠增強(qiáng)扇區(qū)形狀的凸性。Granberg等[21]以斯德哥爾摩阿蘭達(dá)機(jī)場為實(shí)例，為終端區(qū)動態(tài)扇區(qū)劃設(shè)提供了概念驗(yàn)證。

上述研究取得了較為豐富的成果，但主要針對航路扇區(qū)，終端區(qū)扇區(qū)劃設(shè)的研究較少，少數(shù)關(guān)于終端區(qū)扇區(qū)劃設(shè)的研究也未能有效考慮空域功能性。在管制負(fù)荷的衡量方面，研究大多基于靜態(tài)的交通密度，或只選取高峰小時(shí)內(nèi)的航跡數(shù)據(jù)進(jìn)行統(tǒng)計(jì)。《民航空管系統(tǒng)第十三個(gè)五年發(fā)展規(guī)劃》提出了“持續(xù)優(yōu)化管制扇區(qū)構(gòu)型，實(shí)施扇區(qū)功能化劃設(shè)”。為此，將終端區(qū)空域劃分為若干三維空域單元并賦予其確定的功能屬性，以管制負(fù)荷為優(yōu)化目標(biāo)，并考慮空域單元功能的相似性，對空域單元組合形成扇區(qū)。具體地，將管制負(fù)荷分為監(jiān)視負(fù)荷與協(xié)調(diào)負(fù)荷，監(jiān)視負(fù)荷基于動態(tài)密度計(jì)算，協(xié)調(diào)負(fù)荷基于移交航空器架次計(jì)算。

1 問題描述

1.1 空域建模

將終端區(qū)空域分割為若干三維空域單元，每一空域單元內(nèi)包含部分雷達(dá)航跡，設(shè)在某一時(shí)間段內(nèi)的終端區(qū)運(yùn)行的航空器架次為m架，航跡數(shù)據(jù)的集合T可表示為

T={P1,P2,…,Pm}

(1)

每條航跡數(shù)據(jù)的集合可表示為

Pi={pi,1,pi,2,…,pi,n}

(2)

式(2)中：pi,j為第i條航跡的第j個(gè)航跡點(diǎn)，該航跡共有航跡點(diǎn)n個(gè)，按時(shí)間順序排列，每個(gè)航跡點(diǎn)可定義為一個(gè)六維的向量，即

pi,j={x,y,h,d,v,t}

(3)

式(3)中：x為經(jīng)度；y為緯度；h為高度；d為航向；v為速度；t為獲取雷達(dá)數(shù)據(jù)的時(shí)間。

終端區(qū)空域可表示為

(4)

P′={pi1,j1,pi2,j2,…,pin,jn}，

1≤it≤m；1≤jt≤nit

(5)

式中：L(P′)uvz為一個(gè)空域單元格；P′為該單元格中包含的航跡點(diǎn)的集合；m為終端區(qū)所有航跡數(shù)目；nit為航跡Pit中航跡點(diǎn)數(shù)目；u、v和z分別為單元Luvz中心點(diǎn)的經(jīng)度、緯度和高度。

1.2 優(yōu)化目標(biāo)

以最小化扇區(qū)的管制負(fù)荷之和，以及均衡扇區(qū)間的管制負(fù)荷為優(yōu)化目標(biāo)。管制負(fù)荷通常包括監(jiān)視負(fù)荷、協(xié)調(diào)負(fù)荷與沖突負(fù)荷。監(jiān)視負(fù)荷是指管制員檢查扇區(qū)內(nèi)每架航空器的飛行軌跡和飛行情況并對其做出具體管制提示的工作量；協(xié)調(diào)負(fù)荷是指當(dāng)航空器穿越扇區(qū)邊界時(shí)，管制員和飛行員、管制員與管制員之間的信息交換過程中管制員承受的工作量；沖突負(fù)荷為管制員解決扇區(qū)內(nèi)飛行沖突的各種行為的工作量。通過動態(tài)密度的概念將沖突負(fù)荷納入監(jiān)視負(fù)荷中，認(rèn)為管制員對所轄扇區(qū)內(nèi)航班的監(jiān)視包括了潛在沖突的解決等工作。

對于監(jiān)視負(fù)荷，早期研究中的監(jiān)視負(fù)荷通常只使用靜態(tài)密度，即一個(gè)扇區(qū)內(nèi)航空器的架次來表示。然而，同一時(shí)間段內(nèi)并非所有航空器所需要的監(jiān)視工作量都是相同的，飛行狀態(tài)改變的航空器通常需要耗費(fèi)管制員更多的精力，僅用航空器架次衡量監(jiān)視負(fù)荷并不夠準(zhǔn)確。對此，有學(xué)者提出了動態(tài)密度的概念[8-9]，采取各交通復(fù)雜度因子加權(quán)組合的動態(tài)密度反映交通流的復(fù)雜度。給定扇區(qū)i的動態(tài)密度ρi定義為

(6)

式(6)中：N為給定時(shí)間扇區(qū)內(nèi)的航空器架次；NH為航向變化大于15°的航空器架次；NS為速度變化大于18 km/h的航空器架次;NA為高度變化大于250 m的航空器架次;CP為橫向距離小于10 km、垂直距離小于300 m的航空器數(shù)量；W1、W2、W3、W4、W5為需要預(yù)先確定的主觀權(quán)重，分別取1、2.40、2.45、2.94、8。

而協(xié)調(diào)負(fù)荷的指標(biāo)則相對簡單明確，則以離開一個(gè)扇區(qū)的航班數(shù)量表示。綜上，扇區(qū)i的管制負(fù)荷可表示為

(7)

功能性扇區(qū)劃設(shè)的兩個(gè)目標(biāo)可分別表示為

(8)

Obj2=min[SD(Workload,i)]

(9)

式(9)中：SD(Workloadi)為各扇區(qū)間管制負(fù)荷的標(biāo)準(zhǔn)差。

1.3 約束條件

考慮到扇區(qū)劃設(shè)的實(shí)際可行性，在劃設(shè)扇區(qū)的過程中，考慮以下約束：①唯一性約束，即每一個(gè)空域單元只能被分配給一個(gè)扇區(qū)；②連通性約束，即同一扇區(qū)不能不相聯(lián)通，如圖1所示；③方便移交約束，即扇區(qū)邊界距離航路交叉點(diǎn)的距離不能小于10 km，如圖2所示。

圖1 扇區(qū)連通性約束Fig.1 Sector connectivity constraint

圖2 方便移交約束Fig.2 Facilitate handover constraints

2 劃設(shè)方法

2.1 高度分層

圖3 Delaunay三角剖分形成的Voronoi圖Fig.3 Voronoi diagram formed by Delaunay triangulation

目前終端區(qū)扇區(qū)高度上通常是劃為高低兩部分，高扇區(qū)用于處理初始進(jìn)場、最終離場和較高高度飛越航班，而低扇區(qū)用于處理五邊進(jìn)近、初始爬升離場航班。依照先垂直后水平的思路，首先對終端區(qū)內(nèi)的各高度層提取表征航空器運(yùn)動特性的參數(shù)，即平均高度變化率、平均航向變化率和平均速度變化率。接著，對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，采用Q型聚類的方法，將聚類數(shù)目分為2～4類，度量標(biāo)準(zhǔn)采用歐式距離。

2.2 三維空域單元的生成

在將終端區(qū)的高度層進(jìn)行合并的基礎(chǔ)上，分別對合并后的高度層劃設(shè)扇區(qū)。通過Voronoi圖將合并后的高度層水平劃分為若干子單元，各子單元覆蓋相應(yīng)高度，即為三維空域單元。設(shè)有空間目標(biāo)簇，可表示為

P={p1,p2,…，pn}?RN,V(pi)

={y|Dist(y,pi)≤Dist(y,pj)}

(10)

式(10)中：Dist(y,pi)為點(diǎn)y到pi的最小距離；V(pi)為以pi為生長源的Voronoi多邊形；pj為平面中其他生長源；V(pi)={V(p1),V(p2),…，V(pn)}，即為Voronoi圖。

生成Voronoi圖的關(guān)鍵在于生長源pi的選擇。在終端區(qū)空域中，可以將航路點(diǎn)作為自然的生長源，得到空域Voronoi劃分圖，如圖3所示。

在設(shè)計(jì)匹配網(wǎng)絡(luò)之前，根據(jù)頻率間隔和功率大小，做好設(shè)計(jì)規(guī)劃工作，共塔頻率之間的頻率比原理上不得小于1.1，實(shí)際上最好不要小于1.5，否則很難做到理想的匹配網(wǎng)絡(luò)。一個(gè)理想的天調(diào)網(wǎng)路應(yīng)具備穩(wěn)定性和可靠性，匹配狀態(tài)良好，駐波比小于1.1，帶寬足夠大。

結(jié)合高低扇的劃分，需確定各生長源的高度，將其分配至不同的高度層。通過提取各生長源100 m內(nèi)航跡點(diǎn)的高度，采用K-means算法對高度聚類，確定各生長源的高度。部分生長源被多條航線穿越，其附近航跡點(diǎn)具有多類高度。此外，考慮到方便移交約束，需合并生長源距離小于10 km的空域單元，以獲得蛻變的空域Voronoi圖作為扇區(qū)劃設(shè)的基本單元。

2.3 功能性的確定

為了使最終劃設(shè)的功能性扇區(qū)內(nèi)的空域單元具有功能上的同一性或接近性，需確定各空域單元的功能性，各單元的功能取決于單元內(nèi)生長源的屬性。部分航路點(diǎn)本身具有多種功能，此外，Voronoi圖蛻變的過程中會導(dǎo)致有些空域單元包含多個(gè)航路點(diǎn)，因此會存在一個(gè)空域單元具有多個(gè)功能的情況。

在各有限單元進(jìn)行組合形成扇區(qū)的過程中，具有完全相同功能屬性的空域單元進(jìn)行組合時(shí)優(yōu)先級最高；兩空域單元屬性交集不為空集且不完全相同時(shí)優(yōu)先級次之；而功能屬性完全不同的空域單元在組合時(shí)優(yōu)先級最低。

圖4 上海終端區(qū)空域范圍示意圖Fig.4 Schematic diagram of the airspace range of Shanghai terminal area

2.4 算法設(shè)計(jì)

在空域單元組合時(shí)采用NSGA-II算法，具體流程如下。

(1)適應(yīng)度函數(shù)表示。適應(yīng)度函數(shù)即為目標(biāo)函數(shù)。

(2)初始解的產(chǎn)生。按照實(shí)數(shù)編碼規(guī)則產(chǎn)生100個(gè)初始可行解，對于種群中的每一個(gè)個(gè)體，編碼的索引代表空域單元的序號，編碼值表示Voronoi單元所處的扇區(qū)，根據(jù)Voronoi單元的性質(zhì)，在生成編碼后1.3節(jié)中的約束①和③自然滿足，為使得約束②得到滿足，通過以下步驟生成可行解。

步驟1規(guī)定已考慮的單元集A和未考慮的單元集U，初始的A為空集，U包含所有的單元。

步驟2在所有的Voronoi單元中隨機(jī)選擇n個(gè)(最終生成的扇區(qū)數(shù)量)，將其編碼為1～n，并將這n個(gè)單元所處集合修改為A。

步驟3在A中隨機(jī)選擇一個(gè)單元，并在U中隨機(jī)選擇一個(gè)與其相鄰的單元設(shè)置相同的編碼，將選中的單元所處集合修改為A。

步驟4循環(huán)步驟3直到U為空集。接著對初始解集中的每一個(gè)體分別計(jì)算適應(yīng)度值和擁擠度，并按適應(yīng)度的大小進(jìn)行非支配排序。

(3)選擇、產(chǎn)生新解。采用錦標(biāo)賽選擇方法挑選適合繁殖的父代，每次隨機(jī)選擇兩個(gè)個(gè)體，優(yōu)先選擇排序等級高的個(gè)體，如果排序等級一樣，則優(yōu)先選擇擁擠度大的個(gè)體。然后對挑選出的父代的基因進(jìn)行交叉和變異操作，產(chǎn)生新的子代。為保證產(chǎn)生的新解均為可行解，交叉操作選擇單點(diǎn)交叉方式，在隨機(jī)選擇的位置點(diǎn)上進(jìn)行分割并交換右側(cè)的部分，只需判斷交叉點(diǎn)位置的編碼是否滿足約束②。變異操作則將隨機(jī)選擇的位置點(diǎn)處的編碼隨機(jī)變?yōu)榕c該單元相鄰單元的編碼。在修改空域單元所處的扇區(qū)時(shí)，按照不同的優(yōu)先級次序進(jìn)行。

(4)合并子代和父代種群，對新的種群進(jìn)行適應(yīng)度計(jì)算，并選擇前50%的解作為新的解集。

(5)若未達(dá)到迭代次數(shù)，則返回(3)，直至達(dá)到迭代次數(shù)。

3 實(shí)例分析

3.1 上海終端區(qū)運(yùn)行環(huán)境

上海終端區(qū)目前為中國最繁忙的終端管制區(qū)，如圖4所示，其東西長約300 km，南北寬約200 km，管制高度為6 000 m(含)以下，主要負(fù)責(zé)向上海虹橋、浦東國際機(jī)場進(jìn)離場航班，轄區(qū)內(nèi)其他機(jī)場及周邊機(jī)場的飛越航班提供空中交通管制服務(wù)。

目前上海終端區(qū)劃分為11個(gè)扇區(qū)，其北向運(yùn)行時(shí)扇區(qū)劃設(shè)方案如圖5所示。

使用2021年6月14—20日上海終端區(qū)內(nèi)北向運(yùn)行時(shí)雷達(dá)記錄的航跡數(shù)據(jù)。對原始航跡數(shù)據(jù)進(jìn)行預(yù)處理，包括數(shù)據(jù)的清洗、分類及重采樣。

3.2 上海終端區(qū)三維空域單元的生成

根據(jù)2.1節(jié)所述高度分層方法，可得到上海終端區(qū)中各個(gè)高度層進(jìn)場、離場、飛越3種類型航班的運(yùn)行狀態(tài)。其中進(jìn)場航班各高度層內(nèi)的運(yùn)動參數(shù)絕對值如表1所示。

采用SPSS軟件，對各個(gè)高度層進(jìn)行聚類，結(jié)果如圖6所示?？梢钥闯?，對于進(jìn)場航班，高低扇區(qū)界限范圍為2 400～3 000 m。類似地，確定離場與飛越航班高低扇區(qū)界限范圍，最終以2 700 m為高低扇分界線。事實(shí)上，目前上海終端區(qū)高低扇的分界也正是2 700 m。

圖5 上海終端區(qū)北向運(yùn)行時(shí)現(xiàn)行扇區(qū)劃設(shè)方案Fig.5 Current sector plan for the northbound operation of the Shanghai terminal area

Y為標(biāo)準(zhǔn)高度層圖6 進(jìn)場航班垂直高度聚類樹狀圖Fig.6 Arrival flight vertical height clustering dendrogram

水平劃設(shè)時(shí)，以航路點(diǎn)為生長源將終端區(qū)分為若干個(gè)泰森多邊形，以此作為空域單元?？紤]到高低扇區(qū)，在K-means聚類時(shí)令聚類族的數(shù)量k=2。

表1 進(jìn)場航班各高度層參數(shù)統(tǒng)計(jì)Table 1 Parameter statistics of each altitude of the arriving flight

低扇區(qū)共被分配生長源51個(gè)，高扇區(qū)共被分配生長源55個(gè)。通過ArcGis軟件生成以不同生長源為中心的泰森多邊形，如圖7所示。

將生長源小于10 km的兩個(gè)相鄰空域單元合并，獲得蛻變的空域Voronoi圖，如圖8所示。部分生長源高度聚類結(jié)果如表2所示。

圖7 上海終端區(qū)初始空域Voronoi圖Fig.7 Voronoi map of initial airspace in Shanghai terminal area

按照機(jī)場及飛行性質(zhì)的不同將空域單元分為虹橋進(jìn)場、虹橋離場、浦東進(jìn)場、浦東離場以及飛越單元5種。表3為部分空域單元的功能屬性。

圖8 上海終端區(qū)蛻變空域Voronoi圖Fig.8 Voronoi map of transformed airspace in Shanghai terminal area

表2 部分生長源高度的K-means聚類結(jié)果Table 2 K-means clustering results of partial growth source heights

表3 部分空域單元功能屬性示例Table 3 Examples of functional properties of some airspace units

3.3 上海終端區(qū)功能性扇區(qū)的劃設(shè)

通過計(jì)算得到各空域單元內(nèi)的動態(tài)密度分布，可以看出高扇相比低扇的監(jiān)視負(fù)荷分布更加均勻，監(jiān)視密度分布如圖9所示。高扇內(nèi)幾乎每一個(gè)空域單元都分布有一定數(shù)量的監(jiān)視負(fù)荷。而低扇的監(jiān)視負(fù)荷密度顯然分布不均勻，靠近中心的單元密度較大，而外圍單元?jiǎng)t十分稀疏?？紤]到空域容量限制，根據(jù)上海終端區(qū)現(xiàn)行扇區(qū)數(shù)量設(shè)定高扇劃設(shè)7個(gè)扇區(qū)，低扇劃設(shè)4個(gè)扇區(qū)，最后一代種群生成的Pareto解集即為各空域單元組合的優(yōu)解集，如圖10所示。

3.4 結(jié)果分析

由圖9可知，降低各扇區(qū)管制負(fù)荷之和與管制負(fù)荷的標(biāo)準(zhǔn)差出現(xiàn)了此消彼長的現(xiàn)象。經(jīng)計(jì)算，原始扇區(qū)劃設(shè)方案各扇區(qū)管制負(fù)荷之和為183 230.93架次，扇區(qū)間管制負(fù)荷的標(biāo)準(zhǔn)差為6 429.24架次。在Pareto優(yōu)解集中，使得各扇區(qū)管制負(fù)荷之和最低的解可將管制負(fù)荷之和降低7.54%，管制負(fù)荷標(biāo)準(zhǔn)差降低11.19%。而使得各扇區(qū)管制負(fù)荷標(biāo)準(zhǔn)差最小的解可將管制負(fù)荷之和降低5.23%，管制負(fù)荷標(biāo)準(zhǔn)差降低20.61%，如表4所示。

從優(yōu)化效果來看，相比原扇區(qū)劃設(shè)方案，優(yōu)化后各扇區(qū)管制負(fù)荷之和可減少5%～8%；標(biāo)準(zhǔn)差可減少11%～20%，均衡了各扇區(qū)的管制負(fù)荷。這說明原扇區(qū)劃設(shè)方案下各扇區(qū)間管制負(fù)荷之和不大，但各扇區(qū)間管制負(fù)荷卻相差較大，即管制負(fù)荷標(biāo)準(zhǔn)差可優(yōu)化空間更大。因此，選擇各扇區(qū)管制負(fù)荷之和小于原始扇區(qū)的情況下，降低所有扇區(qū)的管制負(fù)荷之和效果更好的解作為上海終端區(qū)優(yōu)化后的扇區(qū)劃設(shè)方案，如圖11所示。

圖9 上海終端區(qū)空域單元監(jiān)視密度分布Fig.9 Distribution of airspace unit surveillance density in Shanghai terminal area

圖10 NSGA-II算法得出的Pareto優(yōu)解集Fig.10 Pareto optimal solution set obtained by NSGA-II algorithm

表4 各扇區(qū)劃設(shè)方案對比Table 4 Comparison of various sector planning schemes

圖11 優(yōu)化后的上海終端區(qū)扇區(qū)劃設(shè)方案Fig.11 The optimized sector plan for Shanghai terminal area

4 結(jié)論

結(jié)合空域的功能性研究了終端區(qū)三維扇區(qū)劃設(shè)問題，以管制扇區(qū)負(fù)荷為優(yōu)化目標(biāo)，考慮了包括連通性約束等在內(nèi)的約束條件。在劃設(shè)方法上，針對不同航班類型通過層次聚類的方法在高度上將終端區(qū)分為高低扇，然后通過泰森多邊形生成若干空域單元，并賦予各空域單元以功能性。采用NSGA-II算法對各空域單元組合生成扇區(qū)。選取了上海終端區(qū)進(jìn)行實(shí)例驗(yàn)證，結(jié)果表明，相比于原扇區(qū)劃設(shè)方案，優(yōu)化后的扇區(qū)劃設(shè)方案最高可將所有扇區(qū)管制負(fù)荷之和降低7.54%，各扇區(qū)管制負(fù)荷的標(biāo)準(zhǔn)差降低20.61%，因此，提出的方法顯示了其在實(shí)際扇區(qū)劃設(shè)中的巨大潛力，未來可以進(jìn)一步研究如何將空域功能性細(xì)分并動態(tài)組合扇區(qū)。