高 偉， 康道馳

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

在中國民航近30年的持續(xù)高速發(fā)展過程中，經(jīng)空域規(guī)劃與管制技術(shù)的不斷升級完善，北京、上海、西安等中國繁忙機(jī)場的終端區(qū)空域通行能力大幅提升，機(jī)場容量的運(yùn)行瓶頸已逐漸轉(zhuǎn)移至終端區(qū)所在區(qū)域管制空域的接收能力。為了改善空域資源供給與利用，需要對區(qū)域管制空域的運(yùn)行服務(wù)能力即區(qū)域容量進(jìn)行評估。

涉及區(qū)域容量的研究主要分為空中交通流量管理、拓?fù)浣Y(jié)構(gòu)、交通流復(fù)雜度和危險(xiǎn)天氣影響4個方面。

空中交通流量管理主要以延誤、安全成本與容量為目標(biāo)對區(qū)域流量進(jìn)行優(yōu)化。Bayen[1]利用拉格朗日算法分析了航路容量和交通流延誤問題，并利用混合整數(shù)線性規(guī)劃(mixed integer linear programming，MILP)算法給出了優(yōu)化后的方案，提出繁忙空域交通流量歐拉網(wǎng)絡(luò)模型和基于空中交通流量管理的拉格朗日延遲模型來計(jì)算區(qū)域容量。Chen等[2]基于動態(tài)空中交通流構(gòu)建航路網(wǎng)絡(luò)模型，采用對偶分析法求解以運(yùn)行花費(fèi)最小作為目標(biāo)函數(shù)，以區(qū)域容量作為限制條件的混合整數(shù)線性規(guī)劃問題，從而實(shí)現(xiàn)對區(qū)域交通流的優(yōu)化。Hossain等[3]提出空域容量與安全風(fēng)險(xiǎn)折中計(jì)算模型，計(jì)算在可接受的安全風(fēng)險(xiǎn)下的區(qū)域扇區(qū)最大容量。

拓?fù)浣Y(jié)構(gòu)主要從航路網(wǎng)絡(luò)復(fù)雜度分析和航路網(wǎng)絡(luò)結(jié)構(gòu)兩方面進(jìn)行研究。網(wǎng)絡(luò)復(fù)雜度分析主要將航路網(wǎng)絡(luò)抽象為拓?fù)渚W(wǎng)絡(luò)，使用復(fù)雜網(wǎng)絡(luò)理論對航路網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)進(jìn)行分析。描述復(fù)雜網(wǎng)絡(luò)的特征量主要包括度、權(quán)、平均路徑長度、簇系數(shù)、度分布、介數(shù)等。陳才龍[4]提出一種空域容量的多目標(biāo)粒子群優(yōu)化方法，通過復(fù)雜網(wǎng)絡(luò)建模，給出了空域容量的一種定量描述方法，最后采用多目標(biāo)粒子群算法得到了航路匯聚點(diǎn)布局問題的非支配解。航路網(wǎng)絡(luò)基本結(jié)構(gòu)的研究主要包括航路交叉點(diǎn)個數(shù)及位置、航段長度、航路交叉角度等與航路網(wǎng)絡(luò)容量間的關(guān)系。Wang等[5]在航路網(wǎng)研究的基礎(chǔ)上考慮了限制區(qū)、危險(xiǎn)區(qū)、禁區(qū)的限制因素，通過對交叉點(diǎn)的合并與位置改變實(shí)現(xiàn)對航路網(wǎng)結(jié)構(gòu)優(yōu)化，并基于元胞自動機(jī)思想求解該模型。

交通流復(fù)雜度方面，Song等[6]提出了基于交通流模式識別的空域容量預(yù)測方法，考慮交通流屬性和管制員的管制策略因素，對未來30 min到幾小時內(nèi)的空域動態(tài)容量進(jìn)行預(yù)測，但并未探討交通模式與容量之間的關(guān)系，以及不同時間尺度內(nèi)的交通流模式的預(yù)測方法。

危險(xiǎn)天氣方面，Mitchell[7]研究了在確定的空域范圍和給定的惡劣天氣影響下的容量算法，并且分析了天氣對容量的影響；Krozel等[8]考慮惡劣天氣和各種運(yùn)行條件進(jìn)行了航路容量評估。

現(xiàn)基于區(qū)域航路網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)對研究區(qū)域容量采用細(xì)胞傳輸模型將區(qū)域航路網(wǎng)中每條航線分為若干細(xì)胞，將時間離散化為若干時段，以細(xì)胞容量作為限制條件，以容量最大作為目標(biāo)函數(shù)構(gòu)建MILP模型，通過遺傳算法對該模型進(jìn)行求解。通過與TAAM仿真軟件的對比來驗(yàn)證該模型的準(zhǔn)確性。

1 細(xì)胞傳輸模型

1.1 細(xì)胞傳輸模型概述

細(xì)胞傳輸模型(cell transmission model，CTM)是運(yùn)動學(xué)模型LWR(lighthill，whitham and richards)模型差分離散形式[9]。經(jīng)過十幾年來的發(fā)展，中外不少學(xué)者對其進(jìn)行了進(jìn)一步研究和完善，CTM 已經(jīng)成為模擬路網(wǎng)中的動態(tài)交通流傳播，研究交通流運(yùn)行規(guī)律的有效工具。將細(xì)胞傳輸模型引入空中交通，并對其進(jìn)行增維處理，加入高度維對區(qū)域空域網(wǎng)絡(luò)進(jìn)行建模。

1.2 定義

圖1所示為區(qū)域航路航線網(wǎng)絡(luò)圖，區(qū)域航線網(wǎng)絡(luò)分為若干高度層，每一高度層上的二維航線網(wǎng)絡(luò)由航線、區(qū)域邊界和航路點(diǎn)構(gòu)成。

1.2.1 點(diǎn)

點(diǎn)包括轉(zhuǎn)彎點(diǎn)、航線與航線的交叉點(diǎn)及區(qū)域邊界與航線的交點(diǎn)，用P表示。在圖1中共有10個點(diǎn)則P={P1,P2,…，P10},其中P2、P6表示轉(zhuǎn)彎點(diǎn)；P3、P9表示航線交叉點(diǎn)；P1、P4、P5、P7、P8、P10表示區(qū)域邊界與航線的交點(diǎn)。

1.2.2 邊

邊分為航段與區(qū)域邊界線。航段為航線上兩點(diǎn)之間的部分用L表示。在圖1中共有9個航段則L={L1,2,L2,3,L3,4,,L8,9,L9,10,L5,9L9,6L6,3,L3,7}。以L1,2為例表示點(diǎn)P1與點(diǎn)P2之間的航線部分，其方向?yàn)閺狞c(diǎn)P1到點(diǎn)P2，航段方向與所在航線方向一致。航線由航段組成，用R表示，航線分為經(jīng)由區(qū)域進(jìn)出終端區(qū)的航線和通過區(qū)域空域但不經(jīng)過終端區(qū)的飛越航線兩類。在圖1中共有3條航線，R={R1,R2,R3}。其中以航線R1為例，航線R1由航段L1,2、L2,3、L3,4組成，R1={L1,2,L2,3,L3,4}，其方向?yàn)镻1→P2→P3→P4。區(qū)域邊界為所選定區(qū)域空域的邊界線在圖1中以粗實(shí)線表示。

1.2.3 高度層

區(qū)域管制空域?yàn)槿S空間，分為若干高度層，每個高度層上的二維航線網(wǎng)絡(luò)均相同。高度用H表示，圖1中共有3個高度層，H={H1,H2,H3}。

圖1 區(qū)域空域航線網(wǎng)絡(luò)圖Fig.1 Regional airspace route network diagram

1.2.4 細(xì)胞

將每個航段按區(qū)域管制間隔分為若干細(xì)胞，用C表示。細(xì)胞作為組成航線網(wǎng)絡(luò)的最基本單元C={C1,C2，…，CN}，N表示航線網(wǎng)絡(luò)中細(xì)胞的總個數(shù)。圖2以航段L1,2為例將航段L1,2分為n個細(xì)胞每個細(xì)胞長度為區(qū)域管制間隔。由于航段長度與區(qū)域管制間隔并非整數(shù)對應(yīng)關(guān)系，所以要對每個航段的最后一個細(xì)胞長度進(jìn)行取整保留。

圖2 航線細(xì)胞結(jié)構(gòu)圖Fig.2 Diagram of route cell structure

(1)

式(1)中：n表示航段L1,2上細(xì)胞的個數(shù)；DL1,2表示航段L1,2的長度；Dcontrol表示區(qū)域管制間隔。

(2)

式(2)中：Clast表示位于航段L1,2上最后的細(xì)胞，其中Clast∈{C1,C2，…，CN}；DClast表示細(xì)胞Clast的長度。

(3)

式(1)中的n為四舍五入后的取整值，當(dāng)航段長度與區(qū)域管制間隔相除的小數(shù)部分小于0.5時，航段上最后一個細(xì)胞Clast的長度如式(2)所示，反之取式(3)。每個高度層上的細(xì)胞分布均相同。

1.3 飛行速度分析

根據(jù)雷達(dá)數(shù)據(jù)分析，影響航空器飛行速度的主要因素為高度層與航段位置。如圖3所示，通過對同航段不同高度層上的航空器速度分布進(jìn)行對比，航空器速度每個高度層上呈正態(tài)分布，但高度層不同，速度的均值與方差也不同。如圖4所示，在相同高度層的不同航段上航空器的速度分布依然呈正態(tài)分布，航段不同，正態(tài)分布的參數(shù)也不同。

圖3 不同高度層速度分布Fig.3 Velocity distribution at different levels

圖4 不同航段速度分布Fig.4 Velocity distribution of different sections

通過以上分析可知，航空器的飛行速度在各航段各高度層上均呈正態(tài)分布，即各航段各高度層上航空器的飛行速度大部分落于某一速度區(qū)間之內(nèi)，以正態(tài)分布的均值作為航空器的飛行速度，不會造成過大的計(jì)算誤差，滿足精度的要求。由于正太分布均值與方差隨高度層和航段位置的變化而變化，速度變量將由航空器所在細(xì)胞的航段位置及高度所決定。

1.4 細(xì)胞傳輸模型構(gòu)建

1.4.1 模型假設(shè)

(1)航空器在細(xì)胞內(nèi)的飛行速度為恒定值，速度由細(xì)胞所在航段位置及高度決定。

(2)將時間離散化，以分鐘作為時間單元進(jìn)行計(jì)算。

(3)模型不考慮天氣等隨機(jī)因素的影響。

(4)區(qū)域飛越航線以平飛為主不發(fā)生高度層的改變，進(jìn)出終端區(qū)航線會發(fā)生高度層的改變。

(5)航空器需要改變高度層時，航空器從細(xì)胞的起始端開始高度變化在完成一個高度層跨越時其飛行的水平距離為一個細(xì)胞的長度，如圖5所示。航空器在改變高度時其水平速度取該細(xì)胞所在兩個高度層速度的平均值。航空器不在交叉點(diǎn)處發(fā)生高度層改變。

圖5 航空器高度層改變示意Fig.5 Elevation level change of aircraft

1.4.2 模型變量計(jì)算

1.4.2.1 飛越細(xì)胞時間

航空器飛越細(xì)胞所用時間由細(xì)胞所在航段位置及高度層決定。

(4)

1.4.2.2 進(jìn)入?yún)^(qū)域交通流

以1 min為時間單元對進(jìn)入?yún)^(qū)域的航空器進(jìn)行標(biāo)定。

(5)

1.4.2.3 高度層

每一架航空器都將生成一個高度層矩陣，高度層矩陣表示航空器在航線的每個細(xì)胞上飛行所使用的高度層。高度層矩陣的維度為高度層數(shù)×航線細(xì)胞數(shù)。

(6)

1.4.2.4 進(jìn)入細(xì)胞時間

當(dāng)式(5)的值為1時航空器進(jìn)入細(xì)胞時間取決于細(xì)胞所在位置及進(jìn)入?yún)^(qū)域空域時間。

(7)

式(6)中高度矩陣為0的元素進(jìn)入時間均為0，表示航空器不從該高度層進(jìn)入該細(xì)胞。

1.4.2.5 離開細(xì)胞時間

(8)

1.4.2.6 交叉點(diǎn)間隔

交叉航線可以看作由多個轉(zhuǎn)彎航線構(gòu)成，求解交叉點(diǎn)間隔則轉(zhuǎn)化為求解轉(zhuǎn)彎點(diǎn)間隔。轉(zhuǎn)彎點(diǎn)示意圖如圖6所示。航空器F1到達(dá)轉(zhuǎn)彎點(diǎn)P1時航空器F2距航空器F1的距離記為d。當(dāng)航空器F1通過點(diǎn)P1時兩航空器之間的距離可以表示為

圖6 轉(zhuǎn)彎航線示意圖Fig.6 Diagram of turning course

(9)

式(9)中：dF1-F2表示航空器飛過P1點(diǎn)后兩航空器之間的距離；時間t以航空器到達(dá)點(diǎn)P1時記作0；a為兩航段的夾角；vF1、vF2為兩航空器速度，其數(shù)值等于航空器所在細(xì)胞的速度。

通過對式(9)求導(dǎo)求解兩航空器之間的最小值，令最小值等于區(qū)域管制間隔。通過對導(dǎo)數(shù)為零與最小值等于區(qū)域管制間隔兩式聯(lián)立解出d。

(10)

過點(diǎn)P1的時間間隔為

(11)

交叉航線可以看作由多個轉(zhuǎn)彎航線構(gòu)成(圖7)，其中交叉航線由n×n種轉(zhuǎn)彎航線組合而成則交叉點(diǎn)間隔取轉(zhuǎn)彎點(diǎn)間隔的最大值，如式(12)所示。

圖7 交叉航線示意圖Fig.7 Diagram of crossing route

(12)

1.4.2.7 細(xì)胞通過量

1.4.2.8 扇區(qū)通過量

扇區(qū)作為區(qū)域空域的組成部分，將扇區(qū)通過量定義為給定時間內(nèi)通過扇區(qū)的航空器數(shù)量，用Sx表示。

(13)

式(13)中：Sx表示給定時間內(nèi)通過扇區(qū)x的航空器數(shù)；C表示扇區(qū)入口處細(xì)胞的集合；K為給定時間集合；H為高度層集合。

1.4.2.9 區(qū)域通過量

區(qū)域通過量與扇區(qū)通過量類似，表示給定時間內(nèi)航空器通過區(qū)域空域的數(shù)量，用A來表示。

(14)

1.4.3 區(qū)域容量

由于以求解區(qū)域最大容量為目的，所以目標(biāo)函數(shù)為一定時間內(nèi)進(jìn)入?yún)^(qū)域航空器數(shù)量最大。限制條件為細(xì)胞容量與交叉點(diǎn)間隔。由于細(xì)胞長度為區(qū)域管制間隔，所以所有細(xì)胞容量均為1。

(15)

限制條件：

(16)

(17)

式(17)中：Pi為航路交叉點(diǎn)的總集合；Cx1、Cx2表示以交叉點(diǎn)Pi為細(xì)胞終點(diǎn)的細(xì)胞集合；Kx1、Kx2分別表示細(xì)胞Cx1、Cx2進(jìn)入?yún)^(qū)域空域的時間。

式(16)表示相鄰兩架航空器后機(jī)進(jìn)入細(xì)胞C前機(jī)進(jìn)入細(xì)胞C的時間差大于細(xì)胞C的通過時間確保了細(xì)胞C的容量為1；式(17)則表示通過交叉點(diǎn)的時間差大于允許值。

2 遺傳算法求解

2.1 遺傳算法概述

遺傳算法是模擬達(dá)爾文遺傳選擇和自然淘汰生物進(jìn)化過程的計(jì)算模型[10]。該算法是一種基于生物自然選擇與遺傳機(jī)理的隨機(jī)搜索算法，起源于對生物系統(tǒng)進(jìn)行的計(jì)算機(jī)模擬研究，是自然遺傳學(xué)和計(jì)算機(jī)科學(xué)相互結(jié)合滲透而成的算法。

與傳統(tǒng)搜索算法不同，遺傳算法[10-11]從一組隨機(jī)產(chǎn)生的初始解開始搜索過程，該初始解稱為種群(population)。種群中的每個個體是問題的一個解，每個個體由攜帶基因編碼的“染色體”組成。遺傳算法通過對初始種群的選擇、交叉、變異運(yùn)算篩選出最優(yōu)個體，即為所求的最優(yōu)解。

2.2 遺傳算法構(gòu)建

2.2.1 初始種群生成

(18)

式(18)中：G(R1)表示個體中航線R1基因序列，矩陣G(R1)中的每一行為個體的一個染色體，染色體的長度為給定時間段的長度Kn。

(19)

式(19)中：G表示個體的基因序列由區(qū)域空域中的每一條航線的基因序列組成。

(20)

式(20)中：Y表示染色體。

個體中的染色體Y的數(shù)目為矩陣G的行數(shù)，即Hn×Rn。對生成個體進(jìn)行檢驗(yàn)使生成的個體需要滿足式(16)、式(17)的限制條件，生成多個滿足條件的個體構(gòu)成初始種群。

2.2.2 適應(yīng)度

適應(yīng)度是在研究自然界中生物的遺傳和進(jìn)化現(xiàn)象時，用來度量某個物種對于生存環(huán)境的適應(yīng)程度。對生存環(huán)境適應(yīng)度較高的物種將獲得更多的繁殖機(jī)會,而對生存環(huán)境適應(yīng)程度較低的物種，其繁殖機(jī)會就會相對較少，甚至逐漸滅絕。本文中個體以給定時間內(nèi)進(jìn)入?yún)^(qū)域空域航班數(shù)的總和為適應(yīng)度，如式(21)所示。

(21)

式(21)中：AG表示個體G的適應(yīng)度。

2.2.3 選擇

選擇指決定以一定的概率從種群中選擇若干個體的操作。一般而言，選擇的過程是一種基于適應(yīng)度的優(yōu)勝劣汰的過程，本文中選擇過程將采用基于適應(yīng)度選擇的輪盤賭算法，即適應(yīng)度越大選擇概率越高，原理類似于博彩游戲中的輪盤賭。個體選擇概率為

(22)

2.2.4 交叉

有性生殖生物在繁殖下一代時2個同源染色體之間通過交叉而重組，即在2個染色體的某一相同位置處DNA被切斷，兩染色體互相交換被分離的染色體片段后組成新的染色體。本文中采用單點(diǎn)交叉，由于單個染色體的長度為給定時間段的長度，所以交叉點(diǎn)的選取范圍為(2,Kn-1)。圖8所示為個體G1與個體G2的Y1進(jìn)行交叉重組。

圖8 染色體交叉示意圖Fig.8 Diagram of chromosome crossing

2.2.5 變異

變異指染色體復(fù)制時可能(很小的概率)產(chǎn)生某些復(fù)制差錯，變異產(chǎn)生新的染色體，表現(xiàn)出新的性狀。本文中則表示個體中染色體上為0的基因以一定概率突變成1。

2.2.6 遺傳算法總流程

3 算例

選取覆蓋北京終端區(qū)的12個區(qū)域扇區(qū)作為研究對象，如圖10所示，共有80個航路點(diǎn)、184個航段、636個細(xì)胞和153條航線其中134條通過區(qū)域航線，10條由區(qū)域進(jìn)入終端區(qū)航線，9條從終端區(qū)進(jìn)入?yún)^(qū)域航線。北京區(qū)域扇區(qū)的管制高度為7 800～12 500 m，按照東單西雙原則向東的航線分為8 100、8 900、9 500、10 100、10 700、11 300、11 900、12 500 m共8個高度層，向西的航線分為7 800、8 400、9 200、9 800、10 400、11 000、11 600、12 200 m共8個高度層。

圖10 選取區(qū)域空域示意圖Fig.10 Diagram of selected regional airspace

3.1 計(jì)算結(jié)果

選取1 h為總時長，計(jì)算區(qū)域空域1 h的最大容量。所需參數(shù)：區(qū)域管制間隔為20 km，初始種群中的個體數(shù)為100交叉概率0.6，變異概率0.01。遺傳算法經(jīng)50次迭代后結(jié)果如圖11所示。最終收斂的均值為884.52，即選定區(qū)域空域的小時容量為884架次。進(jìn)入終端區(qū)和從終端區(qū)進(jìn)入?yún)^(qū)域空域的小時架次116.54，通過區(qū)域空域的小時架次為767.98。

圖11 適應(yīng)度迭代曲線Fig.11 Fitness iteration curve

3.2 TAAM仿真驗(yàn)證

選擇收斂后群體的一個個體作為樣本生成時刻表，通過TAAM對該時刻表的仿真與CTM模型計(jì)算結(jié)果進(jìn)行對比。表1為各扇區(qū)小時通過架次的比較，可以看出，交通流主要集中在12扇區(qū)、16扇區(qū)、17扇區(qū)和19扇區(qū)，總體上計(jì)算值要略高于仿真值。

表1 扇區(qū)小時通過架次統(tǒng)計(jì)Table 1 Sector hour statistics

4 結(jié)論

采用增維細(xì)胞傳輸模型構(gòu)建區(qū)域交通流MILP模型，利用遺傳算法進(jìn)行求解得到區(qū)域空域容量，取得了與TAAM基本一致的結(jié)果。這說明CTM適合區(qū)域交通流這種中觀尺度的空域容量分析，且不需要購置昂貴的仿真平臺，建模簡單，只要給定區(qū)域空域結(jié)構(gòu)，即可快速分析其容量，有助于分析區(qū)域空域?qū)K端區(qū)進(jìn)入?yún)^(qū)域航班的接收能力，從而對合理分析與規(guī)劃機(jī)場容量提供重要決策依據(jù)。該模型只要將有關(guān)場景涉及的細(xì)胞部分稍作處理，還可以得到災(zāi)害性氣象條件或軍事活動影響下的區(qū)域管制空域容量。