沈笑云 王添瑩 張思遠(yuǎn) 焦衛(wèi)東 趙子璇

（中國(guó)民航大學(xué)天津市智能信號(hào)與圖像處理重點(diǎn)實(shí)驗(yàn)室，天津 300300）

1 引言

航路是由國(guó)家統(tǒng)一規(guī)定的通信導(dǎo)航設(shè)備比較齊全的空中通道，寬度通常為20 km［1］，航路網(wǎng)絡(luò)主要由多個(gè)航路點(diǎn)及連接航路點(diǎn)的航段構(gòu)成。近年來(lái)，隨著我國(guó)民航運(yùn)輸業(yè)快速發(fā)展，空中交通網(wǎng)絡(luò)的擁堵問(wèn)題日益嚴(yán)重，航路點(diǎn)的通行能力關(guān)系著整個(gè)航路網(wǎng)絡(luò)的運(yùn)輸能力，正確計(jì)算航路點(diǎn)通行能力對(duì)提高空中交通運(yùn)行效率，減少航班延誤具有重要意義。

通行能力在地面交通中定義為在一定道路幾何條件和交通管制條件下，單位小時(shí)內(nèi)車輛通過(guò)斷面的合理期望最大交通流［2-3］。擴(kuò)展到空域中，通行能力是指空中交通服務(wù)部門在正常運(yùn)行條件下能夠安全地提供空中交通管制服務(wù)的航空器數(shù)量，又稱聲明容量，通常指單位小時(shí)內(nèi)進(jìn)入空域的航空器數(shù)量［1］。

目前空域通行能力的評(píng)估主要有四種方法：基于計(jì)算機(jī)仿真模型的評(píng)估方法［4-5］，基于管制員工作負(fù)荷的評(píng)估方法［6-7］，基于數(shù)學(xué)計(jì)算模型的評(píng)估方法［8-9］，基于歷史統(tǒng)計(jì)數(shù)據(jù)的評(píng)估方法［10-12］。相比而言，基于歷史統(tǒng)計(jì)數(shù)據(jù)的方法更便于操作，結(jié)果更準(zhǔn)確。文獻(xiàn)［10］基于機(jī)場(chǎng)進(jìn)離港航班數(shù)，利用置信區(qū)間篩選處理數(shù)據(jù)并構(gòu)造容量包絡(luò)曲線評(píng)估跑道容量。文獻(xiàn)［11］基于天津機(jī)場(chǎng)小時(shí)進(jìn)離港架次統(tǒng)計(jì)數(shù)據(jù)，對(duì)樣本數(shù)據(jù)進(jìn)行正態(tài)檢驗(yàn)和轉(zhuǎn)換，并構(gòu)造區(qū)間估計(jì)評(píng)估跑道容量。文獻(xiàn)［12］將分位數(shù)回歸應(yīng)用到包絡(luò)曲線容量評(píng)估方法中，對(duì)機(jī)場(chǎng)容量進(jìn)行評(píng)估。文獻(xiàn)［10-12］都使用了基于歷史統(tǒng)計(jì)數(shù)據(jù)的方法，結(jié)果更貼近實(shí)際，但研究對(duì)象都集中于機(jī)場(chǎng)或跑道，沒(méi)有對(duì)航路點(diǎn)的通行能力進(jìn)行研究。在航路網(wǎng)絡(luò)規(guī)劃，優(yōu)化，資源分配等研究中，都以航路點(diǎn)的通行能力值為基礎(chǔ)，且通常利用基于歷史統(tǒng)計(jì)數(shù)據(jù)的飽和流率法來(lái)計(jì)算通行能力［13-14］。文獻(xiàn)［15］使用飽和流率法以一小時(shí)為時(shí)間間隔，統(tǒng)計(jì)航路交叉點(diǎn)小時(shí)流量，并以最大小時(shí)流量值作為航路點(diǎn)通行能力。飽和流率法是地面交通中常用的計(jì)算交叉路口通行能力的方法，但地面交通中又根據(jù)交叉路口的幾何形狀、各個(gè)方向流量的分配、信號(hào)配時(shí)等數(shù)據(jù)，通過(guò)對(duì)飽和流率進(jìn)行修正，計(jì)算得到通行能力［16］。然而空中交通也會(huì)由于飛行計(jì)劃，管制，天氣等因素出現(xiàn)超通行能力飛行的情況，實(shí)測(cè)的最大小時(shí)流率值并不能代表該點(diǎn)的通行能力。廣播式自動(dòng)相關(guān)監(jiān)視（ADS-B）是國(guó)際民航組織確定的未來(lái)主要監(jiān)視技術(shù)，ADS-B 報(bào)文數(shù)據(jù)解碼處理后包含航空器航班號(hào)，時(shí)間，速度，航向角，經(jīng)緯高位置等信息。通過(guò)對(duì)ADS-B 數(shù)據(jù)篩選統(tǒng)計(jì)可得到航路點(diǎn)小時(shí)流量數(shù)據(jù)。

為更準(zhǔn)確地計(jì)算航路點(diǎn)的通行能力，本文利用統(tǒng)計(jì)學(xué)的方法對(duì)歷史航路點(diǎn)實(shí)測(cè)小時(shí)流率進(jìn)行正態(tài)檢驗(yàn)和正態(tài)變換，利用置信區(qū)間修正飽和流率，最終得到航路點(diǎn)的通行能力值，并與基于數(shù)學(xué)模型計(jì)算出的理論通行能力值對(duì)比。

2 航路點(diǎn)通行能力計(jì)算

2.1 航路點(diǎn)理論通行能力數(shù)學(xué)模型

航路網(wǎng)絡(luò)中的節(jié)點(diǎn)通常連接兩條或多條航路，按照航路的交叉方式可將航路模型分為匯聚、發(fā)散、轉(zhuǎn)彎、交叉四種飛行模型［17］，如圖1所示。

航空器A1和A2分別以速度v1，v2（km/h）勻速飛行，α，β，θ，φ為兩相鄰航路的夾角，t1時(shí)刻，航空器A1位于點(diǎn)E，航空器A2位于點(diǎn)F，t2時(shí)刻，航空器A1位于點(diǎn)E'，航空器A2位于點(diǎn)F'。假設(shè)所有航空器均在同一高度層勻速飛行，T為航空器通過(guò)航路點(diǎn)的放行時(shí)間間隔，則通行能力表示為：

以匯聚飛行模型為例，航空器A1先通過(guò)O 點(diǎn)的情況下，設(shè)兩航空器距離為l（km），t=t2-t1，初始時(shí)兩航空器距離即放行間隔為D，在三角形OE'F'中：

令l最小值等于航空器最小飛行間隔L，整理可得最小放行間隔Dmin（km）：

最小放行時(shí)間間隔為Tmin（h）：

放行n架飛機(jī)之后，每?jī)杉茱w機(jī)之間的放行間隔的平均值為：

從而可以計(jì)算出該航路點(diǎn)的理論通行能力為：

同理，其他三種飛行模型的理論通行能力值也可計(jì)算出來(lái)。另外，當(dāng)有多條航路交叉時(shí)，將航路兩兩分組，分別計(jì)算該種模型下的通行能力，取最大值為該點(diǎn)的通行能力。

2.2 飽和流率法

用數(shù)學(xué)模型計(jì)算的結(jié)果偏理論，若基于實(shí)測(cè)數(shù)據(jù)則用飽和流率法計(jì)算航路點(diǎn)通行能力，即一小時(shí)為時(shí)間間隔，統(tǒng)計(jì)航路點(diǎn)每小時(shí)航空器經(jīng)過(guò)的數(shù)量，取最大的小時(shí)流率數(shù)據(jù)為該航路點(diǎn)的通行能力：

式中C為航路點(diǎn)通行能力，fn即第n個(gè)小時(shí)航路點(diǎn)航空器流量，n為小時(shí)流量數(shù)據(jù)即樣本數(shù)據(jù)的數(shù)量。

根據(jù)《中國(guó)航行資料匯編》搜集到航路及航路點(diǎn)信息，選取要研究的航路點(diǎn)，篩選以該航路點(diǎn)為圓心，10 km 為半徑的圓范圍內(nèi)，高度大于航路最低飛行高度的數(shù)據(jù)為該航路點(diǎn)的ADS-B 數(shù)據(jù)。統(tǒng)計(jì)每天每小時(shí)經(jīng)過(guò)該航路點(diǎn)的航空器架次，取最大值為該點(diǎn)的飽和流率參考值。

2.3 飽和流率的修正

由于空中交通飛行多種因素的影響，實(shí)測(cè)的飽和流率值并不能代表航路點(diǎn)的通行能力，所以參考民航局發(fā)布的容量評(píng)估方法，利用統(tǒng)計(jì)學(xué)的方法對(duì)飽和流率值進(jìn)行修正，將實(shí)測(cè)的歷史航路點(diǎn)小時(shí)流量數(shù)據(jù)作為樣本數(shù)據(jù)。在數(shù)理統(tǒng)計(jì)中由中心極限定理，隨機(jī)變量受多個(gè)互不相關(guān)的因素影響時(shí)即認(rèn)為其服從正態(tài)分布，從而構(gòu)建置信區(qū)間進(jìn)行假設(shè)檢驗(yàn)。由于航路點(diǎn)地理位置，重要程度及實(shí)際運(yùn)行情況的不同，會(huì)出現(xiàn)流量值過(guò)大或過(guò)小的情況，導(dǎo)致實(shí)際流量數(shù)據(jù)波動(dòng)過(guò)大，且不符合正態(tài)分布性質(zhì)，若直接進(jìn)行計(jì)算誤差較大。因此要首先對(duì)數(shù)據(jù)進(jìn)行正態(tài)檢驗(yàn)，看數(shù)據(jù)是否符合正態(tài)分布，如不符合則進(jìn)行轉(zhuǎn)換使數(shù)據(jù)具有正態(tài)性，再根據(jù)正態(tài)分布的性質(zhì)對(duì)數(shù)據(jù)進(jìn)行區(qū)間估計(jì)，最后通過(guò)正態(tài)逆變換得到最終的通行能力值。具體步驟如下：

（1）繪制Quantile-quantile 圖（Q-Q 圖）觀察數(shù)據(jù)分布。Q-Q圖展示的是樣本數(shù)據(jù)的分位數(shù)與所指定分布的分位數(shù)之間的關(guān)系，來(lái)檢驗(yàn)數(shù)據(jù)是否符合該分布。若圖中樣本數(shù)據(jù)散點(diǎn)分布在中心直線上，則該組樣本數(shù)據(jù)符合正態(tài)分布［18］，如圖2所示。

（2）Jarque-Bera 正態(tài)檢驗(yàn)。僅通過(guò)繪圖肉眼判斷是否符合正態(tài)分布過(guò)于主觀，為得到更精確的結(jié)果，需要進(jìn)行進(jìn)一步的假設(shè)檢驗(yàn)。Jarque-Bera 檢驗(yàn)適用于樣本數(shù)量大于30 的數(shù)據(jù)，它基于偏度和峰度，相對(duì)檢驗(yàn)更準(zhǔn)確［19］。設(shè)X1，…，Xn是從總體中抽取的樣本，令：

其中Br為r階樣本中心矩，Skew為樣本偏度，用于描述分布偏離對(duì)稱性程度。Kurt 是樣本峰度，是表征概率密度分布曲線在平均值處峰值高低的特征數(shù)。若總體是正態(tài)分布，則Skew應(yīng)近似為0，Kurt近似為3。那么可以構(gòu)造統(tǒng)計(jì)量J：

對(duì)正態(tài)分布而言，統(tǒng)計(jì)量J漸進(jìn)的服從自由度為2的卡方分布，即J=。在顯著性水平a=0.05條件下，可得決策值p：

（3）Johnson 正態(tài)轉(zhuǎn)換。若不符合正態(tài)分布，則需要將數(shù)據(jù)進(jìn)行可逆的正態(tài)變換使其具有正態(tài)性再進(jìn)行下步計(jì)算，根據(jù)統(tǒng)計(jì)學(xué)的相關(guān)研究表明，Johnson 轉(zhuǎn)換相比于Box-Cox 轉(zhuǎn)換，對(duì)數(shù)轉(zhuǎn)換等方法更有效且有更加廣泛的適用范圍，并在絕大多數(shù)參數(shù)取值范圍內(nèi)具有優(yōu)越性，因此本文使用Johnson轉(zhuǎn)換［20］。

Johnson 轉(zhuǎn)換法通過(guò)Johnson 曲線分布族進(jìn)行正態(tài)變化，三個(gè)分布類型分別表示為SB（bounded）、SL（lognormal）和SU（unbounded）［21］，選擇一個(gè)合適的k值。關(guān)于k值的選取，為了轉(zhuǎn)換后的數(shù)列具有較好的正態(tài)性，多次嘗試取值，使其分位值具有代表性，通常取值范圍在0.25 至1.25［22］。通過(guò)標(biāo)準(zhǔn)正態(tài)分布表找出對(duì)應(yīng)于{-3k，-k，k，3k}的分布概率{P-3k，P-k，Pk，P3k}，然后找到對(duì)應(yīng)分位數(shù){x-3k，x-k，xk，x3k}的樣本數(shù)據(jù)，同時(shí)，令M=x3k-xk，N=x-k-x-3k，Q=xk-x-k，并定義分位數(shù)比率QR=，當(dāng)QR<1 時(shí)，選取SB分布；當(dāng)QR=1 時(shí)，選取SL分布；當(dāng)QR>1時(shí)，選取SU分布，各類型的參數(shù)約束及正態(tài)轉(zhuǎn)換方程如表1所示。

表1 Johnson分布系統(tǒng)Tab.1 Johnson Distribution System

經(jīng)轉(zhuǎn)換后的樣本數(shù)據(jù)符合正態(tài)分布，對(duì)于正態(tài)分布Z～N(μ，σ2)，其樣本均值和樣本方差S2為：

a為顯著性水平，在置信水平為1-a的情況下，總體均值μ的置信區(qū)間為：

總體方差σ2的置信區(qū)間為：

設(shè)標(biāo)準(zhǔn)差系數(shù)為w，其值可以通過(guò)查標(biāo)準(zhǔn)正態(tài)分布表得到，覆蓋占比為g，即：

則通行能力的計(jì)算式為C=μ+wσ，可得通行能力的估計(jì)區(qū)間為(C1，C2)，其中：

為方便后續(xù)比較，取C1、C2的平均值Cz，再由表1中的正態(tài)逆變換公式求得Cx，即該航路點(diǎn)的通行能力為Cx。

3 仿真驗(yàn)證

根據(jù)國(guó)際民航數(shù)據(jù)公司OAG 公布的2019 年全球最繁忙航線前十名，我國(guó)的京滬航線位居世界最繁忙航線第七［23］，也是中國(guó)第一繁忙的航線，所以本文選取京滬航線構(gòu)建京滬航路網(wǎng)絡(luò)進(jìn)行通行能力的計(jì)算，根據(jù)中國(guó)航行資料匯編，搜集航路信息，將航路點(diǎn)及導(dǎo)航臺(tái)作為節(jié)點(diǎn)，連接兩航路點(diǎn)或?qū)Ш脚_(tái)的航路作為有向邊，航段長(zhǎng)度為邊權(quán)構(gòu)建有向加權(quán)航路網(wǎng)絡(luò)，如圖3所示。

將解碼后的ADS-B 數(shù)據(jù)實(shí)時(shí)存入數(shù)據(jù)庫(kù)，篩選需要位置和高度范圍內(nèi)的數(shù)據(jù)，設(shè)計(jì)程序統(tǒng)計(jì)空域小時(shí)流量。為驗(yàn)證本文方法的可靠性，統(tǒng)計(jì)了2020年11 月01 日至12 月31 日9：00-21：00 京滬航路網(wǎng)絡(luò)所有航路點(diǎn)的ADS-B 數(shù)據(jù)，剔除航跡點(diǎn)缺失較多的時(shí)段數(shù)據(jù)，統(tǒng)計(jì)航路點(diǎn)小時(shí)航空器流量，處理后每個(gè)航路點(diǎn)有約350條小時(shí)流量數(shù)據(jù)。

以天津（CG）航路點(diǎn)為例，該航路點(diǎn)位于W34航路和V10，V30，W48，W55 航路的交叉點(diǎn)上，如圖4所示。

可以看出該航路點(diǎn)屬于交叉飛行模型，取實(shí)測(cè)航路點(diǎn)ADS-B數(shù)據(jù)每架航空器的平均速度為v，取航空器最小安全間隔為20 km，可得該航路點(diǎn)的理論通行能力值為28.12架次/小時(shí)。統(tǒng)計(jì)天津航路點(diǎn)的小時(shí)架次情況，經(jīng)篩選整理后該航路點(diǎn)共359 條小時(shí)流量樣本數(shù)據(jù)，其直方圖如圖5（a）所示，可得該航路點(diǎn)最大小時(shí)流量即飽和流率為32 架次/小時(shí)。繪制Q-Q 圖觀察該組數(shù)據(jù)，如圖5（b）所示，可以看出圖中頂部和底部的散點(diǎn)距直線有所偏離。

進(jìn)一步進(jìn)行Jarque-Bera 正態(tài)檢驗(yàn)，a為顯著性水平通常取0.05，h是假設(shè)檢驗(yàn)結(jié)果，返回值h=1 表示Jarque-Bera 檢驗(yàn)在默認(rèn)的5%顯著性水平上拒絕原假設(shè)。經(jīng)檢驗(yàn)所得h=1，p值為0.0132<0.05，說(shuō)明天津航路點(diǎn)小時(shí)流量樣本數(shù)據(jù)不符合正態(tài)分布。

利用Johnson 變換進(jìn)行正態(tài)變換，得到的最佳擬合k值為0.79，分位數(shù)比率QR 為0.44628<1，屬于SB分布類型，得到的轉(zhuǎn)換函數(shù)為：

轉(zhuǎn)換后的數(shù)據(jù)直方圖及Q-Q圖如圖6所示。

變換后的Z序列經(jīng)Jarque-Bera檢驗(yàn)所得P值為0.5>0.05符合正態(tài)分布，根據(jù)民航局空域容量評(píng)估指導(dǎo)辦法，本文取顯著性水平a=0.05情況下，95%～98%的覆蓋占比來(lái)構(gòu)造區(qū)間估計(jì)，查標(biāo)準(zhǔn)正態(tài)分布表可知標(biāo)準(zhǔn)差系數(shù)w分別為1.96 和2.33。由樣本數(shù)據(jù)符合SB分布時(shí)的轉(zhuǎn)化函數(shù)，利用表1 中的公式，可由Z值反推出原始樣本X的值，即可得通行能力的值。經(jīng)計(jì)算在95%置信水平下，變換后的Z序列總體均值置信區(qū)間為（-0.08，0.13），總體方差置信區(qū)間為（0.87，1.17），分別求得不同覆蓋占比時(shí)的Z值區(qū)間和通行能力區(qū)間如表2所示。

表2 不同覆蓋占比的通行能力區(qū)間Tab.2 Capacity range of different coverage ratios

考慮到本文采用了兩個(gè)月的樣本數(shù)據(jù)，且空中航空器有穿越航路飛行，或有由于因天氣流量等管制原因臨時(shí)改變航線飛行的情況，為保守估計(jì)航路點(diǎn)通行能力，取置信水平為95%，標(biāo)準(zhǔn)差系數(shù)為1.96 時(shí)的數(shù)值區(qū)間為航路點(diǎn)通行能力。為更準(zhǔn)確的進(jìn)行下一步比較和計(jì)算，取通行能力區(qū)間上下界兩值的平均值作為最終評(píng)估的航路點(diǎn)通行能力，將利用數(shù)學(xué)模型計(jì)算出的航路點(diǎn)理論通行能力作為參考值。

經(jīng)正態(tài)檢驗(yàn)得出鹽城，濟(jì)南，PIMUS 和半塔集航路點(diǎn)的數(shù)據(jù)符合正態(tài)分布，因此這些數(shù)據(jù)無(wú)需進(jìn)行Johnson 變換。其余24個(gè)航路點(diǎn)均不屬于正態(tài)分布，需將其轉(zhuǎn)換為正態(tài)分布。京滬航路網(wǎng)絡(luò)所有航路點(diǎn)修正前后的飽和流率值及數(shù)學(xué)模型法的結(jié)果如表3所示。

表3 京滬航路網(wǎng)絡(luò)航路點(diǎn)各算法通行能力值對(duì)比（單位：架次/小時(shí)）Tab.3 Comparison of the capacity values of each algorithm of the Beijing-Shanghai airway network waypoints（unit：sorties/hour）

分別計(jì)算修正飽和流率值與參考值的誤差以及實(shí)測(cè)飽和流率值與參考值的誤差，并將兩組誤差值進(jìn)行對(duì)比，如圖7所示。

飽和流率值與參考值的平均誤差為6.08 架次/小時(shí)，修正后的平均誤差為3.24 架次/小時(shí)，誤差區(qū)間由［0.45，14.72］縮小到［0.42，7.69］，可以看出修正后的飽和流率相比修正前的更加接近參考值，平均誤差率降低了9%。有時(shí)空域繁忙，會(huì)出現(xiàn)流量大于通行能力值的情況，所以若僅用流量最大值作為該點(diǎn)的通行能力值過(guò)于草率，結(jié)果可能會(huì)偏大。因此，從整體上看，用本文方法在95%置信水平下修正后得到的值更加準(zhǔn)確可靠。

4 結(jié)論

（1）本文對(duì)由ADS-B數(shù)據(jù)統(tǒng)計(jì)得到京滬航路網(wǎng)絡(luò)航路點(diǎn)小時(shí)航空器架次進(jìn)行分析，針對(duì)傳統(tǒng)飽和流率法的不足，參考機(jī)場(chǎng)容量評(píng)估技術(shù)，綜合統(tǒng)計(jì)學(xué)中的方法，提出基于置信區(qū)間的航路點(diǎn)通行能力計(jì)算方法，經(jīng)過(guò)實(shí)例驗(yàn)證，將飽和流率修正前的通行能力值和飽和流率修正后的通行能力值分別與基于數(shù)學(xué)模型計(jì)算出的理論通行能力值對(duì)比，修正前后的結(jié)果與理論值的平均誤差分別為6.08 架次/小時(shí)和3.24 架次/小時(shí)。本文方法相比傳統(tǒng)的飽和流率法結(jié)果更加可靠，魯棒性更好。

（2）另外，通過(guò)實(shí)測(cè)數(shù)據(jù)可以看出部分航路點(diǎn)實(shí)際流量超出理論通行能力很多，這反映了京滬航線空域繁忙，可能因此出現(xiàn)空域擁擠的情況。今后可以在此基礎(chǔ)上進(jìn)行空域通行能力優(yōu)化，空域資源分配等方面的研究。