閆少華 謝曉璇 張兆寧

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

0 引言

近年來，隨著我國空中交通流量不斷增加，空域容量趨近飽和，航班延誤甚至大面積航班延誤時常發(fā)生，空中交通管制壓力不斷增加[1]。為減緩空域壓力，降低航班延誤發(fā)生率，需進(jìn)一步提高空中交通管理能力，提高空中交通流量短時預(yù)測精度?？罩薪煌ǘ虝r流量預(yù)測一般是根據(jù)已有的交通數(shù)據(jù)對未來5，15，30 min或者1 h以上的流量變化情況預(yù)測，為管制員提供輔助決策信息，有利于實現(xiàn)動態(tài)交通控制與引導(dǎo)。

目前國內(nèi)外對空中交通短時流量預(yù)測的方法主要有：①基于飛行計劃的航跡預(yù)測；②考慮空中交通流量非線性特征；③機(jī)器學(xué)習(xí)；④基于數(shù)據(jù)挖掘與智能算法結(jié)合組合模型。

在基于飛行計劃的航跡預(yù)測研究中，陳愷等[2]將計劃航跡和實際航跡關(guān)聯(lián)對比分析，結(jié)合航向角的平面投影提高了航跡預(yù)測的準(zhǔn)確度；毛阿芳[3]通過歷史數(shù)據(jù)進(jìn)行初始航跡預(yù)測，建立航空器運動模型，之后隨著航空器運動狀態(tài)的變化不斷更新航空器運動模型，修正航跡預(yù)測結(jié)果，提高了預(yù)測精度；Pang等[4]考慮了天氣不確定性，提出了基于貝葉斯神經(jīng)網(wǎng)絡(luò)的飛機(jī)軌跡預(yù)測，這種方法能夠利用起飛前的最后1份存檔飛行計劃預(yù)測飛機(jī)軌跡，包括預(yù)測不確定性。在考慮空中交通流量非線性特征的研究中，國內(nèi)外許多學(xué)者通過對時間序列的分析研究中發(fā)現(xiàn)了空中交通流具有十分明顯的非線性特征，例如混沌性、分形性[5-9]；楊陽[10]利用空中交通流量時間序列具有混沌特性，建立了基于神經(jīng)網(wǎng)絡(luò)的流量預(yù)測算法。王飛等[11]利用Hurst指數(shù)驗證了空中交通流量時間序列的分形特性，并且利用相似日的概念對流量進(jìn)行短期預(yù)測，證明了利用分形特征對空中交通流量短期預(yù)測是可行的；在利用機(jī)器學(xué)習(xí)對空中交通流量預(yù)測中，常用的機(jī)器學(xué)習(xí)有：神經(jīng)網(wǎng)絡(luò)模型[12]、K近鄰算法[13]、支持向量機(jī)模型[14]等。其中，李桂毅等[15]在考慮航段相關(guān)性的前提下，利用神經(jīng)網(wǎng)絡(luò)建立了交通流參數(shù)融合預(yù)測模型，預(yù)測了航段流量以及密度參數(shù)；Pang等[16]將卷積層嵌入長短期記憶模型中，根據(jù)對流天氣條件以及飛機(jī)起飛前的飛行計劃預(yù)測飛行軌跡；在基于數(shù)據(jù)挖掘與智能算法結(jié)合組合模型研究方面，Guan等[17]通過ADS-B技術(shù)建立航空大數(shù)據(jù)平臺，將提取的信息映射到路線上，對不同城市之間的空中交通流量進(jìn)行計數(shù)和預(yù)測。尚然然[18]通過分析終端區(qū)流量特性，構(gòu)建基于SOM-k-means算法的相似日聚類模型，最后提出了1種結(jié)合聚類分析與LSTM-BP的終端區(qū)短期流量預(yù)測方法進(jìn)行未來時段的終端區(qū)流量值預(yù)測。

以上4種方法在空中交通流量短時預(yù)測中都有著廣泛的應(yīng)用，但是都存在一定的局限性，其中，基于飛行計劃的航跡預(yù)測的計算復(fù)雜度會隨著流量達(dá)到一定程度時迅速增大，而且航班在實際運行中，會受到各種隨機(jī)因素的干擾，影響航跡預(yù)測精度；考慮空中交通流量非線性特征的預(yù)測方法，參數(shù)設(shè)置對預(yù)測精度影響較大，主觀性較強(qiáng)；單純依賴機(jī)器學(xué)習(xí)的預(yù)測方法只能提取空中交通數(shù)據(jù)的局部特征；基于數(shù)據(jù)挖掘與智能算法結(jié)合組合模型，具有較強(qiáng)自我學(xué)習(xí)能力，但是需要大量數(shù)據(jù)才能挖掘出數(shù)據(jù)間的復(fù)雜關(guān)系，數(shù)據(jù)量少時預(yù)測效果不理想。

針對以上問題，本文提出基于小波優(yōu)化GRU-ARMA預(yù)測方法，首先通過對原始流量數(shù)據(jù)進(jìn)行多尺度不同頻率的小波分解，將原始流量時間序列數(shù)據(jù)看作不同頻率分量的組合，通過小波變換將原始流量時間序列的噪聲提取出來，作為噪聲項，剩余的數(shù)據(jù)作為趨勢項。其中，趨勢項真實反映了原始流量時間序列的隨時間演化的總體趨勢性，使用門控循環(huán)單元（gated recurrent unit，GRU）神經(jīng)網(wǎng)絡(luò)模型進(jìn)行預(yù)測，而噪聲項則是包含了過去一段時間由于危險天氣等突發(fā)狀況的隨機(jī)因素的干擾，使用自回歸移動平均（autoregressive moving average，ARMA）模型進(jìn)行噪聲項的預(yù)測，最后將2個部分結(jié)果疊加得到最終的預(yù)測結(jié)果。

1 預(yù)測模型

1.1 小波變換預(yù)處理原始流量數(shù)據(jù)

考慮到空中交通流量會受到危險天氣、軍方活動等各種時變性較強(qiáng)的隨機(jī)因素的影響，并且這些因素對于空中交通流量的影響程度不同，可將空中交通流量看成各種頻率的分量集合，認(rèn)為不同頻率的分量對空中交通流量的動態(tài)影響具有不同的重要性。小波變換（wavelet transform，WT）是1個專門用于時頻分析的工具，所以筆者先使用小波變換對原始流量數(shù)據(jù)進(jìn)行預(yù)處理，將其分解成低頻子分量的趨勢項和高頻子分量的噪聲項，再用不同模型分別對這些不同頻率的趨勢項和噪聲項進(jìn)行預(yù)測，降低隨機(jī)因素對預(yù)測精度的影響。

小波變換的本質(zhì)是抑制時間序列中的無用信息，增大有用信息占比的過程，本文中通過小波變換降低隨機(jī)因素的干擾，保留真實的流量變化趨勢。小波變換的過程如下。

1）選擇適當(dāng)?shù)男〔ɑ瘮?shù)，確定分解層數(shù)分解原始信號。小波基函數(shù)的種類主要有Haar，dbN，symN，biorN小波等，它們都是小波族，每個小波族里面還包含著具體的小波函數(shù)。在對原始流量數(shù)據(jù)進(jìn)行小波分解的過程中，不同的分解層數(shù)會導(dǎo)致原始流量數(shù)據(jù)表現(xiàn)出不同程度的趨勢性及隨機(jī)性，影響后續(xù)預(yù)測精度。

通過選擇合適的小波基函數(shù)及分解層數(shù)，將原始流量時間序列f(t)分解成多個部分的疊加。

式中：φj，k為近似分量；ψj，k為細(xì)節(jié)分量；cj，k與dj，k為近似分量和細(xì)節(jié)分量對應(yīng)的系數(shù)。近似分量代表著原始流量時間序列的低頻部分，細(xì)節(jié)分量代表著原始流量時間序列的高頻部分，但低頻和高頻是1種相對值，并不表示頻率的絕對大小。原始流量數(shù)據(jù)經(jīng)小波選定的分解層數(shù)進(jìn)行分解之后，每個分解層數(shù)都可以得到對應(yīng)的高頻和低頻信號。

2）設(shè)置合理的閾值，對小波系數(shù)進(jìn)行估算。本文選用啟發(fā)式閾值確定閾值，綜合固定閾值與無偏風(fēng)險估計閾值的優(yōu)勢。

如果crit＞eta，選擇固定閾值作為閾值；如果crit≤eta，則比較固定閾值和無偏風(fēng)險估計閾值，選擇較小的值。

其中，固定閾值為

式中：N為信號的長度；σ為噪聲方差，通過魯棒中值定理[19]計算得出。

無偏風(fēng)險估計閾值的計算分為3個步驟。

步驟1。把信號s進(jìn)行絕對值和升序處理，得到新的信號序列。

步驟2。如果設(shè)置閾值為f(k)的第k個值的平方根，即

則該閾值產(chǎn)生的風(fēng)險為

步驟3。根據(jù)Risk(k)，取最小風(fēng)險點相對應(yīng)的值k(min)，無偏風(fēng)險估計閾值定義為

3）小波重構(gòu)。將原始流量時間序列f(t)分解成趨勢項trend(t)與噪聲項noise(t)之和。

1.2 GRU神經(jīng)網(wǎng)絡(luò)模型預(yù)測趨勢項

GRU神經(jīng)網(wǎng)絡(luò)模型是1種適用于時間序列預(yù)測的神經(jīng)網(wǎng)絡(luò)模型，包括控制上一時刻到當(dāng)前時刻信息的量的更新門（update gate），以及決定了上一時刻信息在當(dāng)前時刻的寫入程度的重置門（reset gate）。通過引入更新門和重置門決定提取歷史流量數(shù)據(jù)的量以及如何與新的流量數(shù)據(jù)結(jié)合，對趨勢項進(jìn)行預(yù)測，其算法原理圖見圖1。

圖1 GRU算法原理圖Fig.1 The structure of Gated Recurrent Unit（GRU）

GRU的計算見式（10）。

式中：Xt為t時刻該層的輸入變量；Rt為重置門，Zt為更新門；Ht為隱含狀態(tài)；Ht-1為候選隱含狀態(tài)；Wxr，Wxz，Wxh∈Rd×h和Whr，Whz，Whh∈Rh×h為權(quán)重參數(shù)；br，bz，bh∈R1×h為偏差參數(shù)；h為隱藏單元個數(shù)；d為輸入變量維度；σ(x)和tanhx都是激活函數(shù)；⊙表示矩陣乘法。

1.3 ARMA模型預(yù)測噪聲項

ARMA模型是時間序列預(yù)測的基本模型之一，利用ARMA模型預(yù)測噪聲項。令Xt為t時刻的真實值，假設(shè)Xt不僅與t時刻之前的真實值Xt-1，Xt-2，…，Xt-p有關(guān)，還與t時刻之前的干擾值εt-1，εt-2，…，εt-q相關(guān)，則Xt可以用式（11）表示。

式中：φ1，φ2，…，φp與θ1，θ2，…，θq，均為該線性組合的系數(shù)。

要求Xt是平穩(wěn)時間序列，εt是白噪聲序列，且滿足式（11）時，將{}Xt表示為p階自回歸與q階滑動的結(jié)合序列，記為ARMA(p，q)。

式中：E是Xt的數(shù)學(xué)期望，[Var]是Xt的方差。

1.4 基于小波優(yōu)化GRU-ARMA的預(yù)測模型

基于小波優(yōu)化GRU-ARMA模型的空中交通流量短時預(yù)測的流程圖見圖2。首先利用小波變換將原始流量數(shù)據(jù)分解為不同頻率的趨勢項和噪聲項，分別使用GRU神經(jīng)網(wǎng)絡(luò)模型和ARMA模型進(jìn)行預(yù)測。將2個模型的預(yù)測結(jié)果疊加得到流量預(yù)測結(jié)果。

圖2 基于小波優(yōu)化GRU-ARMA模型的空中交通流量短時預(yù)測流程Fig.2 Short term air traffic flow frediction process based on Wavelet-Optimized GRU-ARMAmodel

1.5 評價指標(biāo)

將小波優(yōu)化GRU-ARMA模型預(yù)測值與真實值對比，對誤差情況進(jìn)行分析。本文選用均方根誤差（RMSE）、平均絕對值誤差（MAE）、平均絕對百分比誤差（MAPE）作為誤差評價指標(biāo)，計算見式（13）～（15）。

2 實例分析

2.1 數(shù)據(jù)說明

本文的數(shù)據(jù)來源于中國某終端區(qū)2019年4月29日—6月3日的ADS-B數(shù)據(jù)，共計36 d的數(shù)據(jù)，時間間隔取5 min，共10 368組數(shù)據(jù)。

2.2 小波變換分析

小波基函數(shù)及分解層數(shù)的選擇沒有明確的規(guī)則，多靠經(jīng)驗及實驗結(jié)果確定。選擇bior2.2，db3和sym4作為小波基函數(shù)，對原始流量數(shù)據(jù)進(jìn)行3，4，5層的小波分解，信噪比計算結(jié)果見表1。由表1可見：在3種小波基函數(shù)不同分解層數(shù)的信噪比計算中，分解層數(shù)為3信噪比達(dá)到最大，信噪比越大，分解效果越好，說明分解層數(shù)設(shè)置為3時分解效果最好，所以在后續(xù)分析過程中，只考慮分解尺度為3的情況。

表1 不同分解層數(shù)信噪比對比Tab.1 Comparison of different decomposition layers of signal-to-noise ratio

分別選取小波基函數(shù)bior函數(shù)、db函數(shù)和sym函數(shù)對選取的終端區(qū)流量數(shù)據(jù)進(jìn)行小波分解，計算相應(yīng)的信噪比，根據(jù)信噪比計算結(jié)果選擇最優(yōu)的小波基函數(shù)，結(jié)果見圖3。

圖3 小波基函數(shù)信噪比Fig.3 Signal noise ratio of wavelet basis function

根據(jù)信噪比計算結(jié)果，選定小波基函數(shù)bior2.8，分解層數(shù)為3，采用啟發(fā)式閾值進(jìn)行閾值處理，對原始流量時間序列進(jìn)行小波分解處理，得到原始流量時間序列的近似序列和細(xì)節(jié)序列，結(jié)果見圖4。

圖4 原始流量時間序列小波分解圖Fig.4 Wavelet decomposition of original traffic time series

2.3 GRU神經(jīng)網(wǎng)絡(luò)模型對趨勢項進(jìn)行預(yù)測

基于Python的keras神經(jīng)網(wǎng)絡(luò)庫構(gòu)建GRU神經(jīng)網(wǎng)絡(luò)模型，通過loss曲線（見圖5）對模型進(jìn)行診斷，batch size設(shè)置為64，迭代次數(shù)為500。小波分解處理過后得到的趨勢項數(shù)據(jù)作為GRU神經(jīng)網(wǎng)絡(luò)模型的數(shù)據(jù)輸入，將前30 d的以5 min為間隔的數(shù)據(jù)作為訓(xùn)練集，31～35 d的以5 min為間隔的數(shù)據(jù)作為測試集，對第36 d的前2 h流量進(jìn)行預(yù)測，預(yù)測結(jié)果見圖6。

圖5 GRU神經(jīng)網(wǎng)絡(luò)loss曲線圖Fig.5 Loss graph of GRU

圖6 趨勢項預(yù)測結(jié)果Fig.6 Prediction results of trend items

2.4 ARMA（p，q）模型對噪聲項預(yù)測

利用ARMA模型對噪聲項進(jìn)行預(yù)測，分為4個步驟。

1）數(shù)據(jù)平穩(wěn)性檢驗。ARMA模型要求數(shù)據(jù)是平穩(wěn)的，首先要對噪聲項時間序列進(jìn)行平穩(wěn)性檢驗。本文使用增廣迪基-富勒檢驗（ADF）單位根進(jìn)行平穩(wěn)性檢驗，通過Python中的adftest函數(shù)計算得到該噪聲序列數(shù)據(jù)的ADF值為-5.360，見表2，該噪聲序列在90%，95%，99%的置信區(qū)間下通過平穩(wěn)性檢驗，說明可以使用ARMA模型。

表2 不同置信區(qū)間對應(yīng)的臨界ADF值Tab.2 Corresponding critical ADF value under different confidence intervals

2）模型定階。根據(jù)貝葉斯信息準(zhǔn)則（Bayesian information criterion，BIC）[20]繪制噪聲項時間序列的熱力圖，見圖7，當(dāng)AR系數(shù)p設(shè)為4、MA系數(shù)q設(shè)為3時，BIC的值最小，對應(yīng)方塊的顏色是最深的，所以確定模型為ARMA（4，3）。

圖7 噪聲項ARMA定階熱力圖Fig.7 ARMAthermal diagram of noise term

3）模型檢驗。為了確定ARMA模型的階數(shù)適用于噪聲項時間序列的預(yù)測，還需要對ARMA（4，3）模型進(jìn)行進(jìn)一步的殘差驗證。殘差就是原始信號去除模型根據(jù)原始信息擬合出的信號后剩下的信號。如果殘差是隨機(jī)而且是正態(tài)分布的、不存在自相關(guān)情況的，就說明殘差屬于白噪聲信號，也就代表了建立的ARMA（4，3）模型已經(jīng)將所有有用的信號都包含了。本文利用Durbin-Watson[21]對殘差進(jìn)行檢驗，計算出的結(jié)果為2.006 5，這個值與2接近程度越高越說明殘差通過檢驗，不存在相關(guān)性。

4）噪聲項預(yù)測。噪聲項的預(yù)測與趨勢項的預(yù)測相對應(yīng)，將10 368個噪聲項數(shù)據(jù)分為864組，每組12個，將前680組數(shù)據(jù)作為訓(xùn)練數(shù)據(jù)，將第681～862組數(shù)據(jù)作為測試數(shù)據(jù)，對后2組數(shù)據(jù)值進(jìn)行預(yù)測，得到噪聲項預(yù)測結(jié)果。

2.5 總的流量預(yù)測

將GRU神經(jīng)網(wǎng)絡(luò)模型預(yù)測的趨勢項數(shù)據(jù)和ARMA模型預(yù)測的噪聲項數(shù)據(jù)疊加得到總的預(yù)測值，計算出每個預(yù)測點的相對誤差，見圖8。

圖8 小波優(yōu)化GRU-ARMA預(yù)測結(jié)果Fig.8 Prediction results of the wavelet-optimized GRU-ARMAmodel

1）不同模型對比。為了驗證小波優(yōu)化GRU-ARMA模型的預(yù)測效果及優(yōu)越性，本文將小波優(yōu)化GRU-ARMA模型、單一的GRU神經(jīng)網(wǎng)絡(luò)模型、BiLSTM模型、CNN-LSTM模型和ARMA模型的預(yù)測值與真實值進(jìn)行對比，結(jié)果見圖9?？梢钥闯龌谛〔▋?yōu)化GRU-ARMA模型預(yù)測結(jié)果更接近真實值。

圖9 不同模型預(yù)測結(jié)果Fig.9 Prediction results of different models

2）誤差分析。為了進(jìn)一步表明小波優(yōu)化GRU-ARMA模型在真實值與預(yù)測值的接近程度相對于傳統(tǒng)的神經(jīng)網(wǎng)絡(luò)模型和ARMA模型的優(yōu)越性，本文對比了5種模型在預(yù)測每個點的相對誤差，見圖10。由圖10可見：基于小波優(yōu)化GRU-ARMA模型在每個點的相對誤差基本都保持在2%左右，其他4種模型由于直接使用原始流量數(shù)據(jù)進(jìn)行預(yù)測，每個點的誤差都比較大，基本保持在5%～15%，GRU神經(jīng)網(wǎng)絡(luò)模型最大誤差達(dá)到了28.57%，BiLSTM最大誤差達(dá)到了37.14%，CNN-LSTM模型和ARMA模型最大誤差達(dá)到34.29%，說明基于小波優(yōu)化GRU-ARMA模型的組合模型預(yù)測效果相對于直接使用原始數(shù)據(jù)進(jìn)行預(yù)測的模型更穩(wěn)定。

圖10 5種模型誤差對比Fig.10 Error comparison of five models

3）模型評價。計算5種模型的RMSE，MAE和MAPE，見表3。由表3可見：基于小波優(yōu)化GRU-ARMA模型的預(yù)測精度最高，MAPE值為1.74%，與直接使用原始數(shù)據(jù)進(jìn)行預(yù)測的GRU，BiLSTM，CNN-LSTM，ARMA模型相比，預(yù)測精度分別提高了3.02%，5.39%，5.05%，4.30%，更適合空中交通流量短時預(yù)測。

表3 5種模型的評價指標(biāo)Tab.3 Evaluation indexes of four models

3 結(jié)束語

實驗結(jié)果表明，基于小波優(yōu)化GRU-ARMA模型的空中交通流量短時預(yù)測方法通過小波變換將原始流量數(shù)據(jù)分成不同頻率的趨勢項和噪聲項2個部分，去除噪聲的趨勢項能夠更好的反映空中交通流量的變化規(guī)律，使用了GRU神經(jīng)網(wǎng)絡(luò)模型對趨勢項時間序列及ARMA模型對噪聲項時間序列進(jìn)行預(yù)測，結(jié)合2個模型的預(yù)測結(jié)果得到最終的預(yù)測值。實例分析表明：相比于直接采用原始數(shù)據(jù)進(jìn)行預(yù)測的傳統(tǒng)神經(jīng)網(wǎng)絡(luò)模型和單一的ARMA模型，基于小波優(yōu)化的GRU-ARMA模型預(yù)測誤差更小，解決了空中交通流量預(yù)測受隨機(jī)因素干擾的影響，與直接使用原始數(shù)據(jù)進(jìn)行預(yù)測的GRU，BiLSTM，CNN-LSTM，ARMA模型相比，預(yù)測效果更加理想，對于空中交通短時流量預(yù)測效果有所提升。

本文研究是對空中交通流量進(jìn)行短時預(yù)測，未來將進(jìn)一步，研究在多方面因素干擾下交通量的動態(tài)性預(yù)測，為空中交通動態(tài)管理提供依據(jù)。