基于風集合預報的扇區概率擁堵預測

徐子玥，胡明華，張 穎，王 兵，謝 華，丁文浩

（南京航空航天大學民航學院，江蘇 南京 211106）

各類對實際運行產生干擾的因素中，天氣是最主要的不確定干擾因素，包括惡劣天氣和不確定性風[1]。 在航跡預測中風的不確定性直接影響航跡預測的準確性和可靠性，表現為航空器位置、過點時間、燃料消耗的離散[2]。 本文著重考慮不確定性風的影響，量化該不確定性，轉化為扇區需求預測的不確定性，以提高ATM 系統的穩定性和可預測性，減少運行成本，提高ATFM 策略的有效性。

如今，量化天氣不確定性的趨勢是使用集合預報系統（EPS）。 EPS 在ATM 問題中常見應用方法可分為以下兩類[2]：

1） 集合預報方法：將集合預報中的概率信息轉換為用戶相關信息，通過相關信息的統計特征值，量化天氣不確定性的大小。 Steiner 等[3]將對流天氣的不確定性轉化為可用空域容量的減少；Kicinger等[4]采用該方法量化天氣不確定性對飛行時間的影響，以衡量天氣不確定性在選擇最優路線方面的成本；Kim 等[5]將天氣不確定性轉化為軌跡的不確定性，軌跡集合提供了關于飛行參數的不確定性信息（飛行持續時間和航行燃料成本），將該信息用于軌跡選擇；Valenzuela 等[6]考慮風和溫度的不確定性，轉化為扇區需求的不確定性；Valenzuela 等[6-7]分析了降低軌跡不確定性對扇區需求的影響。

2） 概率轉換法：將天氣的概率分布函數（PDF）轉化為相關信息的概率分布函數。Rivas 等[8]，Franco等[9]考慮多航段，將風的PDF 轉化為飛行時間和燃料消耗的PDF；Rivas 等[8，10]考慮單航段，將風的PDF轉化為燃料消耗的PDF，預測了考慮風預報不確定性的飛行時間及飛行燃油消耗情況。

綜上，本文研究了考慮風不確定性的集合軌跡預測方法以及預測軌跡的不確定性分析方法，并在此基礎上建立了基于集合軌跡預測進行扇區交通需求集合預測， 以及扇區概率擁擠指標計算的方法， 進而獲得風不確定性影響下的軌跡預測集合。基于所預測的扇區概率擁堵，可以實施更為科學有效且具有魯棒性的流量管理策略。

1 考慮風不確定性的航跡預測

1.1 集合預報

為了表征和量化天氣預報固有的不確定性，常采用概率天氣預報，其制作方案[2]包括：基于歷史誤差的方案，預報員的主觀解釋，確定型預報的統計后處理，集合預報系統（EPS）。 其中，解決ATM 相關問題的趨勢是使用EPS。 EPS 通過多次運行初始條件稍微不同的確定性數值天氣預報（NWP）/稍微擾動的天氣模型[11]，時滯[12]，多模型[13]等方法構建集合，理想情況下實際天氣結果應在集合范圍內。

目前運行中常用的EPS 產品包括[5]：

1） 小規模大規模集合天氣預報（PEARP）。法國運行全球集合預報模型，由35 個成員組成，每天在0600UTC 和1800UTC 運行。 PEARP 的預報時間范圍為中短期（4～5 d），該模型的輸出時間間隔為6 h。利用初始條件的奇異向量和集合數據同化，構建34個擾動集合成員。

2）氣象局全球和區域集合預測系統（MOGREPS）。由12 個成員組成，每天在0000UTC、0600UTC、1200UTC 和1800UTC 運行。MOGREPS 的預報時間范圍為短期（1～2 d），輸出時間間隔為3 h。 該模型利用集合變換卡爾曼濾波器，構建11 個擾動集合成員的初始條件。

3） 歐洲中期天氣預報中心（ECMWF）。ECMWF的綜合預報系統（IFS）由51 個成員組成，每天在0000UTC 和1200UTC 運行。ECMWF 的預報時間范圍為中期（15 d），輸出時間間隔為6 h。 該模型利用初始條件的擾動（奇異向量技術）、集合數據量同化方法和隨機擾動模型的物理參數， 構建50 個擾動集合成員。

4） 超級集合（SUPER）。 結合上述3 種EPS 構建的多模型集合系統。 該模型將PEARP 的35 名成員、MOGREPS 的12 名成員和IFS 的51 名成員混合在一起， 形成98 名成員。 該模型在1800UTC 運行。 SUPER 的預報時間范圍42 h，輸出時間間隔為6 h（所有組成模型的共同間隔）。

1.2 集合軌跡預測

EPS 應用于航跡預測相關問題的常見方法有兩種：

1） 集合預報法（EWF）：將集合中的天氣成員分別輸入確定性航跡預測器（dTP），輸出航跡集合。 航跡集合中的成員分別對應于相應的天氣成員場景，如圖1 所示。

圖1 集合預報法示意圖Fig.1 Schematic diagram of EWF

EWF 主張從集合每個成員中提取與航空相關的特征，而不是將集合預報提供的豐富信息總結成概率天氣描述。 強調每個集合成員的重要性，分析面向集合全體成員，不局限于集合成員的平均值。 不確定性通過集合成員的離散度（最大值-最小值）來表征。

2） 概率轉換法（PTM）：擬合集合中的天氣成員得到概率天氣分布，輸入概率性航跡預測器（pTP），輸出概率航跡，如圖2 所示。

圖2 概率轉換法示意圖Fig.2 Schematic diagram of PTM

PTM 允許通過給出的這些隨機變量的聯合概率密度函數，計算隨機變量轉換后的聯合概率密度函數。

本文擬采用的EWF 是目前應用最廣泛的概率預測方法，是一種能夠產生一系列未來天氣可能性的技術。 相較于PTM，EWF 計算成本更低，更適合于實際應用。

1.3 確定性航跡預測模型

本文采用基于動力學模型中應用最廣泛的點質量模型法（PMM）[14-17]，構建考慮風不確定性的航跡預測模型，用作EWF 法中的確定性航跡預測器。

1.3.1 構建模型的基本假設

已知航跡的描述有4 個維度： 經度、 緯度、高度、時間。 假設如下。

1） 航空器以恒定空速V=0.84 馬赫， 保持恒定高度H=9 800 m 水平飛行。巡航階段常采用固定馬赫數的方式勻速飛行， 該假設固定巡航高度層，避免高度變化導致速度不同，影響馬赫數的變化。

2） 航空器處于多航段的巡航飛行階段，航路點之間視作直線飛行。 該假設忽略航空器的轉彎過程，避免轉彎過程帶來的速度變化，保證航空器以恒定速度飛行。

3） 大氣環境：在國際標準大氣（ISA）模型基礎上疊加EPS 提供的集合預報風數據，將風的不確定性納入模型考慮范疇。

4） 運動學方程中不考慮橫向動力學， 將側風轉化為等效逆風。 本文中假設航空器沿固定航路水平飛行，不存在航跡偏離，忽略風的不確定性對水平航跡的影響，即經過的預定航路及航路點不發生變化。

1.3.2 模型的基本結構

本文著眼于預戰術階段，著重考慮風不確定性的影響，預測得到巡航飛行階段的航跡集合。 該模型的結構如圖3 所示。 輸入部分包括：

圖3 航跡預測模型結構示意圖Fig.3 Structure diagram of trajectory prediction model

1） 航空器意圖數據。為模型提供簡化駕駛員或飛行管理系統操縱航空器的抽象描述。 預戰術階段發生在運行前1～6 d，飛行計劃（FPLs）通常尚未可用，主要信息來源是航空器意圖（FIs）數據，包含：航班號、出發地和目的地機場、預計起飛時間、航空公司和飛機類型。 在本文中，該數據來源于航空公司運營數據以及導航數據。

2） 航空器性能參數。為模型航空器中動力學方程提供性能方面的參數值。 包含：飛行包線（最大速度、最小速度等），空氣動力學（機翼面積和阻力系數），發動機推力和油耗等參數。 本文中，該數據來源于航空器基礎資料（BADA）數據庫。

3） 氣象數據。 為模型提供相關環境信息。 本文考慮風不確定性的影響，風向風速的數據來源于集合預報。

本文僅研究航跡在時間維度上的不確定性——風的不確定性影響航空器的地速，表現為航空器過點時間的不確定性。

模型中航跡以航路點及其相應過點時間的形式來表示，模型輸出的航跡集合即為各航路點過點時間的集合列表。

1.3.3 數學模型

已知各網格點上的氣象數據，非網格點處航路點上的風值可通過三維線性插值得到（圖4）。 具體操作步驟如下：

圖4 三維線性插值示意圖Fig.4 Schematic diagram of 3D linear interpolation

1） 找到該航路點周圍的8 個網格點；

2） 進行高度上的線性插值；

3） 進行緯度上的線性插值；

4） 進行經度上的線性插值。

根據前文假設航空器保持恒定高度飛行，為簡化插值過程，航空器在航段內各位置處的風值由航段兩端航路點的緯/經向風值線性插值（取決于航空器在航段上的飛行距離）后矢量求和得到。航班i（i=1，2，…，I），集合成員m（m=1，2，…，M）在航段j（j=1，2，…，J）上任一點所受到的緯/經向風uij[m]（r）/vij[m]（r），r 為航空器在該航段上的飛行距離， 可由航段兩端航路點的u 分量和v 分量線性插值得到，數學表達式可寫作

式中：rij為航班i 航段j 的長度，m；us，ij[m]/vs，ij[m]（ue，ij[m]/ve，ij[m]）分別為航班i 集合成員m 在航段j 的航路起點（終點）的緯度和經度方向的風分量。

該點處的高空風wij[m]（r）可由緯度/經度方向風分量求和得到，數學表達式可寫作

圖5 風的矢量三角形示意圖Fig.5 Vector triangular diagram of wind

根據假設4，結合風的矢量三角形定理，航班i集合成員m 在航段j 上飛行距離為r 時的地速Vg，ij[m]（r）數學表達式為

航班i 集合成員m 經過各航路點p（p=0，1，…，J）的過點時間的數學表達式為

式中：p=0 時該點為航班i 的起時航路點；ti0為航班i 的預計起飛時間。

1.4 軌跡過點時間不確定性分析

1.4.1 定性分析

已知EWF 方法采用集合成員的離散度Spread（即最大值-最小值）來量化不確定性。

1） 風的不確定性：該不確定性具有時空相關性。

①空間相關性：風的不確定性常隨地理位置的變化而變化；

②時間相關性：風的不確定性常隨預報時間提前量（LAT）的增大而增大，呈正相關關系。

2） 航跡過點時間的不確定性： 表征為預測航跡集合的離散度， 即對于航班航跡上的任一航路點p，該點處航跡過點時間不確定性大小的數學表達式為

根據式（7）中航跡過點時間公式可知，過點時間是以航空器預計起飛時間為起始點，逐個航段推算獲得的， 風不確定性的影響沿著航跡不斷累積，即航跡的不確定性與航空器到該航路點的飛行距離呈正相關。

1.4.2 定量分析

采用逐步回歸法構建多元線性回歸模型，分析上述因素與扇區進入時間不確定性之間的定量關系。 已知定性分析中影響航跡不確定性的因素，結合扇區進入時間不確定性問題將其特殊化，則定量分析中， 影響航班的扇區進入時間不確定性的因素，即輸入多元線性回歸模型的解釋變量有

1） 扇區進入點pienter（經過獨熱編碼處理）；

2） 預報時間提前量Lienter；

3） 進入扇區時的飛行距離dienter。

逐步回歸的基本思想是將變量逐個引入模型，每次引入對因變量影響最顯著的解釋變量，當原來引入的解釋變量由于后面解釋變量的引入變得不再顯著時，則將其刪除。 這是一個反復的過程，直到既沒有顯著的解釋變量選入回歸方程，也沒有不顯著的解釋變量從回歸方程中剔除為止。 該方法保證將統計上不顯著的解釋變量剔除，最后保留在模型中的解釋變量之間多重共線性不明顯，而且對被解釋變量有較好的解釋貢獻。

具體步驟如下：

Step 1 將n 個解釋變量xn分別引入模型，同因變量y 構建n 個一元回歸模型。 計算各解釋變量相應的F 檢驗統計量的值。 當其最大值超過給定顯著水平對應的臨界值時，對應的解釋變量引入選入變量指標集合。

Step 2 將剩余n-1 個解釋變量分別與選入變量指標一起引入模型，同因變量y 構建n-1 個二元回歸模型，計算得到相應的n-1 個F 檢驗統計量的值， 當其最大值超過給定顯著水平對應的臨界值時，將新的解釋變量引入選入變量指標集合，否則終止引入過程。

Step 3 重復步驟2，構建n-h（h≤n）個（h+1）元回歸模型， 直到F 檢驗統計量的最大值小于給定顯著水平對應的臨界值或h=n 時， 終止引入過程。

2 扇區交通需求及扇區擁堵預測

2.1 扇區交通需求集合預測及分析

將第1 節得到的航跡集合，轉化為扇區需求集合，并根據扇區需求集合的特征統計量，分析考慮風不確定性的扇區需求。

本文將扇區需求定義為：選定時間段Pk內進入扇區的航空器數量。 Pk的數學表達式如下

對所預測的扇區需求集合中的各預測值進行統計分析，計算如下統計指標值：最大值（Max）、最小值（Min）、平均值（Avg）、離散度（Spread），計算式如下

2.2 扇區擁堵指標

扇區擁堵定義為：某時間段內，扇區需求超過扇區聲明容量（Capacity）。 本文中扇區聲明容量假設為定值，扇區需求存在不確定性，通過匹配確定性容量和不確定性需求判斷扇區擁堵。

由于ATFM 措施的嚴重性取決于扇區擁堵的可能性及嚴重性，本文定義扇區擁堵指標，扇區擁堵概率，預期扇區容量缺失值如下

式中：G[m]（Pk）為選定時間段Pk內，集合成員m 的扇區容量缺失值，定義為扇區需求與扇區聲明容量的差值，為正值時表示扇區擁堵；C 為扇區容量。

引入單位階躍函數的概念用于指示扇區是否擁堵，數學表達式為

3 實例計算

3.1 數據

1） 扇區：本文選取ZUUUAR04 扇區作為對象，分析扇區概率需求， 已知該扇區的聲明容量為28架/h，7 架/15 min。

2） 航班：北京時間2019 年6 月8 日00:00—24:00 途經扇區ZUUUAR04 的航班，共計429 架。扇區 進 入 點 有：CDX （229 架/53.4%）、GYN （1 架/0.2%）、OMBON（192 架/44.8%）、P247（7 架/1.6%）。

3） 集合天氣預報：ECMWF 在UTC 時間2019年6 月7 日00：00 的預報。

3.2 確定性航跡與不確定性航跡集合的對比

本文采用確定性模型進行考慮風預報不確定性的航跡預測，風預報的不確定性通過EWF 方法來量化。 根據前文所述ECMWF 的集合構建方法，可知風集合預報中成員0 為初始天氣預報， 成員1～50 為擾動成員。 現將成員0 得到的確定性航跡， 與考慮風預報不確定性得到的不確定性航跡（成員0～50 得到的確定性航跡集合） 對比如圖6所示。

圖6 確定性航跡與不確定性航跡集合Fig.6 Determined trajectory vs undetermined trajectory set

由圖6 可知， 橫坐標表示該航跡的航路點編號，縱坐標表示該航跡的預計過點時間。 紅色實線表示集合成員0 預測得到的確定性過點時間，紫色/藍色直線分別代表集合成員預測得到的最小/最大過點時間，紫色藍色直線范圍內即為航跡的可能過點時間區間，相較于確定性過點時間，增加了預測航跡的魯棒性。

3.3 不確定性分析

3.3.1 軌跡過點時間不確定性分析

1） 定性分析。繪制各航班進入扇區時的飛行距離、預報時間提前量與扇區進入時間離散度的散點圖，分別如圖7，圖8 所示。 藍色實心點分別為各航班進入扇區時的飛行距離、LAT 及其相應的進入扇區時間的離散度， 紅色實線是散點擬合出的趨勢線。 可以看出扇區進入時間離散度與航班進入扇區時的飛行距離、預報時間提前量成正相關關系。

圖7 扇區進入時間不確定性與飛行距離關系圖Fig.7 Dispersion of the entry time vs flight distance to the entry point

圖8 扇區進入時間不確定性與LAT 關系圖Fig.8 Dispersion of the entry time vs LAT

2） 定量分析。 于SPSS 平臺上，基于逐步回歸分析法構建多元線性回歸模型，分析各因素對扇區進入時間離散度的影響情況。

本文用于判定是否引入該解釋變量的臨界值設為0.1，表示當候選變量中最大F 值的P 值小于或等于0.1 時，引入相關變量。在引入方程的變量中，最小F 值的P 值大于或等于0.1 時，則剔除該變量。

逐步回歸過程中每一步驟輸入/除去解釋變量，構建多元線性回歸模型的結果摘要如表1 所示。

表1 模型摘要Tab.1 Model summary

經過逐步回歸， 模型3 中修正后的R2=0.874，大于0.8，表示模型擬合效果較好；F 值為28.022，伴隨回歸的顯著性檢驗0.000＜0.05， 說明逐步回歸法調整后的方程聯合顯著性檢驗通過。 模型3 擬合得到的多元線性回歸方程為

式中：y 為航班進入扇區時間的離散度；x1為航班進入扇區時的飛行距離dienter；x2為航班進入扇區時的預報提前時間Lienter；x3為航班進入扇區時的航路點OMBON。

3.3.2 扇區需求不確定性分析

所預測扇區需求的不確定性特征影響扇區概率擁堵的計算，本節重點分析不同的統計時間段周期長度T 下的扇區需求的不確定性特征。

當選定時間段周期T=60 min，扇區進入計數的統計特征量（Min、Max、Avg）如圖9 所示。

圖9 T=60 min 時各時間段內扇區需求Fig.9 Sector demand for T=60 min

扇區進入計數的平均值在22:00—23:00 處取最大值，考慮到扇區聲明容量為28 架/h（圖中用紅線表示）， 可知在21:00—22:00 和22:00—23:00 這兩個時間段內，累積持續時間120 min，扇區需求大于扇區聲明容量，存在擁堵的可能性。

各選定時間段內， 扇區進入計數的不確定性（離散度）如圖10 所示。

圖10 T=60 min 時各時間段內扇區需求離散度Fig.10 Dispersion of the sector demand for T=60 min

離散度在22:00—23:00 時間段處取最大值，為3 架。 相較于扇區的聲明容量28 架/h，大約占比11%。

當選定時間段周期時，扇區進入計數的統計特征量（Min、Max、Avg）如圖11 所示。

圖11 T=15 min 時各時間段內扇區需求Fig.11 Sector demand for T=15 min

扇區聲明容量隨選定時間段周期成比例減小為7 架/15 min，如圖12 所示，共計11 個時間段內，累積持續時間165 min， 扇區進入計數大于扇區聲明容量，存在擁堵的可能性。

各選定時間段內， 扇區進入計數的不確定性（離散度） 如圖12 所示。 離散度在02:45—3:00、22:15—22:30、22:30—22:45、23:00—23:15 這4 個時間段處取最大值，為3 架。

圖12 T=15 min 時各時間段內扇區需求離散度Fig.12 Dispersion of the sector demand for T=15 min

比較上述兩種情況，選定時間段周期T 減小時，扇區進入計數成比例縮小，但離散度變化不大，即扇區需求成比例縮小， 需求不確定性的重要性顯著提高。 與此同時，根據扇區平均需求計算的超出容量事件的發生頻率提高，累積超出容量持續時間變長。

3.4 扇區概率擁堵指標分析

分別計算得出選定時間段周期T=60，15 min時，各時間段內扇區擁堵概率及其相應的預期容量缺失值，并將其可視化，如圖13 所示。

圖13 各時間段內扇區擁堵概率及預期容量缺失值Fig.13 Sector congestion probability and expected capacity deficit values for each time period

圖13 中， 藍色柱狀圖表示該時間段內扇區擁堵概率，紅色折線圖表示該時間段內預期的扇區容量缺失值。

從圖13 可以看出，當T=60 min 時，在22：00—23：00 時間段內； 當T=15 min 時， 在23：00—23：15時間段內，扇區擁堵概率最大，同時扇區容量缺失值最大，扇區概率擁堵風險最大，ATFCM 策略最嚴格。

4 結論

1） 研究了考慮風預報不確定性的集合軌跡預測方法，并分析了軌跡預測的不確定性隨預報時間提前量以及飛行距離的變化趨勢。 實例驗證可知，軌跡預測的不確定性與預報時間提前量、飛行距離均成正相關關系。

2） 研究了航班過點時間預測不確定性量化方法， 基于逐步回歸法建立了多元線性回歸方程，經檢驗該模型具有顯著性，能量化各個因素對航跡過點時間預測不確定性的貢獻度。

3） 研究了基于集合軌跡預測進行扇區交通需求集合預測以及扇區概率擁堵指標計算的方法，并分析了扇區概率擁堵相關的兩個指標隨不同統計時間段粒度的變化特征。 實例計算結果表明：當選定時間段周期減小時，扇區擁堵的概率增大，扇區擁堵概率達到閾值的累積持續時間增大。