王 飛

（中國民航大學空管學院，天津 300300）

空中交通流是大量航空器沿航路航線相繼飛行形成的航班流，是空中交通流量管理的主要對象.對空中交通流時空特性的剖析是分析空中交通流系統(tǒng)內(nèi)在特性、揭示空中交通流系統(tǒng)動態(tài)演化規(guī)律的重要基礎(chǔ)[1]，在空中交通流建模、預測和復雜性分析等領(lǐng)域具有廣泛應用前景.

當前，空中交通流的研究主要從3方面開展：基于數(shù)學模型的分析[2-3]，是對系統(tǒng)的高度抽象，難以真實反映空中交通流運行實際；基于系統(tǒng)仿真的研究[4-5]，考慮了不同研究對象交互作用對空中交通流的影響，但仍難以模擬不確定性因素；基于實測數(shù)據(jù)的數(shù)據(jù)挖掘分析[6-7]，可以直接從實際運行數(shù)據(jù)中挖掘空中交通流的特征，更為真實地反映空中交通流內(nèi)在規(guī)律性.

隨著研究的逐步深入，國內(nèi)外學者逐漸意識到空中交通流系統(tǒng)類似于地面交通流系統(tǒng)，也是一個復雜的非線性動力系統(tǒng)，并開始應用非線性動力學理論研究空中交通流.雖然非線性理論在道路交通方面得到廣泛研究和應用，但在空中交通領(lǐng)域的研究剛剛起步，研究成果較少.李善梅等提出應用小數(shù)據(jù)量法和小波去噪理論計算最大Lyapunov指數(shù)，并對飛行沖突時間序列的混沌特性進行了分析[8].Cong等以扇區(qū)交通流量為對象建立時間序列數(shù)據(jù)，并對其混沌特性進行了深入研究[9].鄭旭芳等從混沌與分形角度對交通流量時間序列的非線性特征進行了對比分析，并分析了時間尺度的影響[10].楊陽等利用混沌特性研究了扇區(qū)流量的短期預測方法[11].

傳統(tǒng)的空中交通流建模與預測方法均是基于統(tǒng)計意義上的分布函數(shù)，鑒于空中交通流系統(tǒng)已被證實具有混沌、分形等非線性特性，傳統(tǒng)方法過于理想化，難以準確把握空中交通流的內(nèi)在特征和演化規(guī)律.從現(xiàn)有研究可以看出，對空中交通流的混沌識別與預測研究較多，而對于分形特征的研究也僅停留在計算關(guān)聯(lián)維數(shù)這一表象，對其內(nèi)在的分形特征都尚未開展研究.基于實際運行數(shù)據(jù)構(gòu)建的時間序列是研究非線性系統(tǒng)的有效手段.研究空中交通流時間序列的分形特征，一方面為現(xiàn)有空中交通流建模方法是否合理提供判據(jù)，另一方面為應用分形理論實現(xiàn)空中交通流短期預測奠定基礎(chǔ).

本文以基于實際運行數(shù)據(jù)的管制扇區(qū)流量時間序列為研究對象，首先對時間序列進行非線性檢驗，然后系統(tǒng)性分析時間序列的分形特征，最后計算了時間序列的關(guān)聯(lián)維數(shù).

1 空中交通流時間序列

本文采用文獻[11]中的三亞飛行情報區(qū)2016年5月某天零時開始連續(xù)12 d的流量統(tǒng)計數(shù)據(jù).為全面分析空中交通流的非線性分形特征，以三亞管制區(qū)1號扇區(qū)為對象，將數(shù)據(jù)處理為時間尺度30、15、10 min和5 min的4個交通流時間序列，流量時間序列如圖1所示.

圖 1 不同時間尺度的空中交通流量時間序列Fig.1 Air traffic flow time series based on different time scales

2 空中交通流的非線性檢驗

檢驗空中交通流時間序列是否具備非線性是應用非線性理論的必要條件.替代數(shù)據(jù)法是檢驗時間序列非線性的常用檢驗方法，包括零假設(shè)、替代數(shù)據(jù)生成和計算檢驗統(tǒng)計量3個部分組成.該方法基本思路是：假設(shè)觀測的時間序列是線性的，在保持原有線性關(guān)系下根據(jù)原始分布的性質(zhì)（如均值、方差等）生成一組替代數(shù)據(jù)，然后計算原始數(shù)據(jù)和替代數(shù)據(jù)的檢驗統(tǒng)計量.若兩者檢驗統(tǒng)計量有顯著差異，則拒絕零假設(shè)，說明原始時間序列不太可能是從假設(shè)的線性系統(tǒng)中產(chǎn)生，即原始數(shù)據(jù)中存在確定的非線性成分.

（1）零假設(shè)

觀測數(shù)據(jù)由具有原始數(shù)據(jù)的均值和方差的線性相關(guān)高斯過程產(chǎn)生[12].如果拒絕該零假設(shè)，則待檢驗時間序列包含非線性成分.

（2）生成替代數(shù)據(jù)

保留原時間序列功率譜增幅，通過改變功率譜中各頻率成份的相位來生成替代數(shù)據(jù).具體步驟如下：

步驟1將原始時間序列A進行傅里葉變換，生成新的序列X= {Xi}，i=1,2,···,n，n為原始序列長度.

步驟2隨機生成新的功率譜相位序列F={fi} ，fi∈ [0,2π]；

步驟3應用隨機相位按照式（1）改造X，構(gòu)建新的時間序列Y={Yi} ；

步驟4對時間序列Y進行逆傅里葉變換，得到新的時間序列，作為替代數(shù)據(jù).替代數(shù)據(jù)和原有數(shù)據(jù)具有相同的功率譜和自相關(guān)函數(shù).

（3）檢驗統(tǒng)計量

檢驗統(tǒng)計量的選擇將會直接影響判斷結(jié)果.本文選取KS檢驗作為檢驗統(tǒng)計量，該方法是一種有效、穩(wěn)定的非線性檢驗方法[13].

應用替代數(shù)據(jù)方法對4個交通流時間序列進行非線性檢驗.由于功率譜相位序列是隨機產(chǎn)生，每次產(chǎn)生的功率譜相位序列都不一樣，結(jié)合蒙特卡洛思想，每個時間序列均進行了1 000次檢驗（顯著性水平為0.05），具體檢驗結(jié)果見表1.

從表1可看出，不同時間尺度的時間序列非線性檢驗結(jié)果是不一樣的，即隨著時間尺度的增大，具備非線性的概率越小.需要說明的是，小概率非線性并不代表一定不是非線性，還應當結(jié)合其它指標（例如關(guān)聯(lián)維度）聯(lián)合檢驗.而時間尺度小于5 min的時間序列是否依然具備非線性，后續(xù)將另文研究.

表 1 時間序列非線性檢驗結(jié)果Tab.1 Nonlinear test results of time series

本文尺度為5 min的空中交通流時間序列可以確定具有非線性成分，但是否具備分形特征還需進一步分析.

3 空中交通流的分形特征分析

分形是指由各部分組成的形態(tài)，每個部分以某種方式與整體相似.交通流分形現(xiàn)象從整體來說包括自組織、自相似和無標度，從局部來說具有多重分形.根據(jù)這4類現(xiàn)象，歸納出空中交通流分形特征包括自相似、長相關(guān)、無標度和多重分形.

3.1 自相似性特征分析

空中交通流的自相關(guān)性并不是嚴格意義上的完全重合，而是在不同的空間或時間尺度的統(tǒng)計意義上的自相似.為了避免噪聲干擾，本文采用具有優(yōu)良的時域和頻域局域化能力的小波變換，對原始時間序列進行小波多尺度分解，分解尺度為3，提取其中低頻系數(shù)結(jié)果，如圖2所示.

圖 2 小波分解結(jié)果Fig.2 Results of wavelet decomposition

從圖2中可以看出，在不同日期里空中交通流呈現(xiàn)出類似“周期”的變化，顯示出較強的自相似性.當分解尺度為1、2時，可得到相同結(jié)論.

3.2 長相關(guān)性特征分析

長相關(guān)性也稱為長程相關(guān)性或記憶性，是產(chǎn)生分形的根源.長相關(guān)性反映的是過去的狀態(tài)對當前和未來的狀態(tài)的影響是長期的，這也是能夠根據(jù)過去狀態(tài)進行交通流預測的一個主要原因.R/S方法是分析時間序列長相關(guān)性的有效方法，通過計算Hurst指數(shù)H來判斷時間序列是否具有長程相關(guān)性.

3.2.1 Hurst指數(shù)

Hurst指數(shù)反映了時間序列A均值的累計離差隨時間變化的范圍.

采用重標級差法[14]先計算出重標極差序列及其對應的時間增量序列，繪制曲線圖，并運用回歸分析計算出斜率，即為H.

根據(jù)統(tǒng)計特性：H= 0.5表示觀測序列是變化無規(guī)律，即隨機的；H＜ 0.5表示觀測序列是反持續(xù)的時間序列，即前一時點上升（下降），下一時點將有下降（上升）的趨勢；H＞ 0.5表示所觀測序列是持續(xù)性序列，即前一時點上升（下降），在下一階段要繼續(xù)這種正（負）的變化趨勢.由于Hurst指數(shù)能反映序列的持續(xù)性或隨機性，在實際中它得到廣泛的應用.

根據(jù)時間序列的時間跨度不同，Hurst指數(shù)分為全局Hurst指數(shù)和局部Hurst指數(shù).全局Hurst指數(shù)表示整個時間序列長相關(guān)性，在空中交通流中表示交通流保持之前趨勢的強弱性.局部Hurst指數(shù)反映局部交通流的變化情況.兩者的計算方法是相同的.

該時間序列的全局Hurst指數(shù)大于0.5，該時間序列并不是隨機的，并且具有正相關(guān)性，即每一個觀測值都保存著之前所有事情的記憶，未來的情況也與現(xiàn)在的情況有巨大關(guān)聯(lián)性，說明該時間序列具有長相關(guān)性.

圖 3 Hurst指數(shù)擬合曲線Fig.3 Fitting curve of Hurst exponent

為了檢驗R/S分析的可靠性，將原始時間序列的數(shù)值順序進行隨機打亂，形成新的時間序列.若原始時間序列是獨立的，則新的時間序列的Hurst指數(shù)應基本保持不變，反映原始時間序列沒有長相關(guān)性[14].按照重標級差法經(jīng)過1 000次計算得到新時間序列的Hurst指數(shù)在（0.55，0.62）之間，相較于原始時間序列有大幅下降，說明原始時間序列不是獨立隨機序列，R/S分析結(jié)果是可靠的.

3.2.3 局部Hurst指數(shù)分析

為反映交通流隨時間的變化規(guī)律，將原始時間序列按照288個數(shù)據(jù)（1 d）為一個子區(qū)間，以1個數(shù)據(jù)為間隔向后推，得到一系列時間序列子區(qū)間，計算各子時間序列的Hurst指數(shù)，如圖4所示.

圖 4 局部Hurst指數(shù)Fig.4 Local Hurst exponents

局部Hurst指數(shù)均大于0.5，表現(xiàn)出較強的長相關(guān)性，說明交通流能夠保持之前的變化趨勢.每一天周期內(nèi)均會出現(xiàn)幾個局部低值，表現(xiàn)出長相關(guān)性減弱、隨機性增強，是交通流趨勢最可能出現(xiàn)反轉(zhuǎn)的點，反映的是交通流處于擁堵階段.

3.3 無標度性特征分析

分形往往在一個最小標度和最大標度構(gòu)成的區(qū)間內(nèi)才存在分形規(guī)律，這個區(qū)間就是無標度區(qū)間.一旦逾越無標度區(qū)間，自相似性便不復存在，系統(tǒng)也就沒有分形規(guī)律了.因此，無標度區(qū)間的確定一直是分形研究的熱點問題.

無標度區(qū)間通常是與關(guān)聯(lián)維數(shù)計算緊密聯(lián)系的，GP (Grassbegr-Procaccia)算法是計算關(guān)聯(lián)維數(shù)的常用方法，計算步驟如下：

步驟1對原始時間序列A采用時間差法重構(gòu)相空間，形成新的時間序列X2為

式中：A= {A(i) }；m為嵌入維數(shù)，可人為設(shè)置數(shù)值，僅針對該嵌入維數(shù)進行相空間重構(gòu)； τ1為時間序列的延遲時間，根據(jù)自相關(guān)系數(shù)法計算得到.

步驟2計算相空間中的每一個向量點A(i) 與其它向量點的距離，然后統(tǒng)計在以A(i) 為中心、r為半徑的圓球中的向量點的個數(shù)，從而得到關(guān)聯(lián)積分C(m,r)為

式中： ‖A(i)?A(j)‖ 表示相空間中兩個向量點之間的距離；H(·) 表示Heaviside函數(shù)。

步驟3繪制lnr- l nC(m,r) 曲線圖，尋找線性度最好的區(qū)間 [lnr1,lnr2] ，即無標度區(qū)間，lnr1、lnr2分別為該區(qū)間的端點.

目前，研究者大多采用視覺法識別無標度區(qū)間，簡單易行，但主觀性和隨意性較強.文獻[15]應用一階差分數(shù)據(jù)k-means方法識別出無標度區(qū)間，但該方法對具有非線性動力方程的系統(tǒng)，如Lorenz等具有較好的識別效果，但對于諸如交通流時間序列的離散數(shù)據(jù)的識別效果不盡如人意.文獻[16]提出應用雙對數(shù)曲線點的二階差分數(shù)據(jù)，識別出零值附近波動數(shù)據(jù)，再剔除粗大誤差，用統(tǒng)計方法識別出無標度區(qū)間.本文采用文獻[16]的方法，對空中交通流時間序列的 l nC(m,r) 差分數(shù)據(jù)進行分析.

式中：lnC′(m,ri)和lnC′′(m,ri)分別表示lnC(m,r)的一階差分和二階差分.

以m= 14，τ1= 3為例，應用GP算法繪制的雙對數(shù)關(guān)聯(lián)積分曲線和對應的差分如圖5所示.

圖 5 雙對數(shù)關(guān)聯(lián)積分和二階差分曲線Fig.5 Curves of double-logarithmic correlation integral and second-order difference

從圖5可以看出，由lnr- l nC′(m,r) 曲線計算得到區(qū)間為 l nr∈ [1.631 4,2.1910] ，從lnr-lnC′′(m,r)曲線可以看出，在0附近直線的lnr∈[1.949 8，2.260 0].綜合考慮兩個曲線圖，得到本時間序列的無標度區(qū)間為 [1.949 8,2.191 0] .不同嵌入維度下的無標度區(qū)間如表2所示.

表 2 不同嵌入維度下的無標度區(qū)間Tab.2 Non-scale range under different embedding dimensions

從表2可以看出，不同嵌入維數(shù)對應的無標度區(qū)間是有區(qū)別的，但并未發(fā)現(xiàn)無標度區(qū)間與嵌入維數(shù)具有明顯規(guī)律性.

3.4 多重分形特征分析

若時間序列僅具有某種單一自相似特征，這樣的分形稱為簡單分形，只需選擇一個標度就可以刻畫.多重分形則描述多標度復合分形的復雜體系.許多分形具有更復雜的尺度關(guān)系，需要一組參數(shù)去描述復雜系統(tǒng)特性，這就是多重分形.空中交通流時間序列的多重分形特征可以應用配分函數(shù)和多重分形譜來進行分析.具體步驟如下：

步驟1將原始時間序列A等分為N個長度為?t的小區(qū)間，定義流量增量為

步驟2將配分函數(shù)定義為流量增量的q階矩的加權(quán)和，表示為

步驟3針對每一個q值，繪制 l g(?t) -lg(Sq(?t))雙對數(shù)曲線圖，并利用最小二乘法進行線性擬合，其斜率為質(zhì)量指數(shù) τq，如式（8）所示.

步驟4將τq進行Legendre轉(zhuǎn)換，即可得到多重分形譜為

式中： α為奇異性指數(shù)，用來描述各個子區(qū)間之間不同的奇異強度，f( α )為具有相同的子集的分數(shù)維度，通過分形譜函數(shù)可以刻畫不同α值對應的分形特征.

時間序列計算得到的多重分形譜如圖6所示.

多重分形譜構(gòu)成的幾何圖形一般表現(xiàn)為鉤狀或者鐘狀，包含 6個參數(shù).其中，參數(shù) αmax和 αmin分別表示最大概率子集和最小概率子集的奇異指數(shù)；參數(shù) ?α=αmax?αmin表示多重分形譜的寬度，可以定量表征時間序列波動的強度；參數(shù)f(αmin) 和f(αmax) 分別表示最小概率子集和最大概率子集的分形維數(shù)；參數(shù) ?f=f(αmax)?f(αmin) 表示最大概率子集和最小概率子集的分形差異，反映時間序列的變化趨勢.當?f＞0時表明時間序列數(shù)據(jù)向高數(shù)值方向發(fā)展的可能性較大，當 ?f＜0 時則表明時間序列數(shù)據(jù)向低數(shù)值方向發(fā)展的可能性較大.

從圖6可以看出，原始時間序列的 α與f(α) 的關(guān)系呈現(xiàn)“鐘形”，說明存在多重分形特征[16].?α為0.042 6，說明交通流量并不是均勻分布，具有一定波動性.?f為0.069 2，說明交通流量處于高流量的機會比處于低流量的機會大.

圖 6 時間序列的多重分形譜Fig.6 Multifractal spectrum for time series

通常有兩種類型的多重分形，一種是由時間序列的概率密度函數(shù)的“肥尾”分布引起的，另一種是由時間序列的長相關(guān)性引起的.判斷一個時間序列產(chǎn)生多重分形的本質(zhì)根源，可通過隨機打亂原始時間序列元素順序來分析.對于第一種情況，隨機打亂順序也無法消除序列的多重分形性質(zhì)；對于第二種情況，隨機打亂順序?qū)е聲r間序列的長相關(guān)性被破壞，往往會呈現(xiàn)非多重分形性質(zhì).如果這兩種情況同時存在，則打亂順序的時間序列依舊呈現(xiàn)多重分形性質(zhì)，但其多重分形性相較于原始序列較弱.

從圖6可以看出，打亂順序后的時間序列仍然具有“鐘形”的多重分形特征，并且 ?α和 ?f都有所減弱，因此原始交通流的多重分形性質(zhì)是由長相關(guān)性和概率密度的“肥尾”分布共同決定的.

4 空中交通流的分形維數(shù)計算

分形維數(shù)是定量刻畫分形特征的參數(shù)，是分形的最基本特征.由于研究對象的不同，分形維數(shù)有多種定義，常用的包括Hausdorff維數(shù)、相似維數(shù)、信息維數(shù)、Lyapunov維數(shù)和關(guān)聯(lián)維數(shù)等.雖然這些維數(shù)定義不同，但只要計算出的維數(shù)大于1，即可判定研究對象具有分形特征.

本文采用GP算法來計算空中交通流時間序列的關(guān)聯(lián)維數(shù)D2.按照3.3節(jié)所述方法，首先設(shè)置m=2，繪制 lnr- l nC(r) 圖，在無標度區(qū)間內(nèi)用最小二乘擬合得到斜率，即為此時的關(guān)聯(lián)維數(shù)估計D2(2) .然后，逐步增加嵌入維數(shù)，直到關(guān)聯(lián)維數(shù)不再隨著m的增加而增加.然后繪制m-D2(m) 圖來獲得關(guān)聯(lián)維數(shù).本文時間序列的lnr- l nC(r) 關(guān)系見圖7，m-D2(m) 關(guān)系見圖8.

圖 7 雙對數(shù)關(guān)聯(lián)積分曲線Fig.7 Curves of double-logarithmic correlation integral

圖 8 關(guān)聯(lián)維數(shù)隨嵌入維數(shù)的變化曲線Fig.8 Curve of the correlation dimensionD2with the increasing embedding dimensionm

從圖8可以看出，當m= 10時，關(guān)聯(lián)維數(shù)D2達到最大值，此后D2隨m的增加而不再增加.因此，m= 10對應的關(guān)聯(lián)維度是該時間序列的關(guān)聯(lián)維度，即D2=6.89 ，說明需要7個變量才能清晰描述該空中交通流.

5 結(jié) 論

（1）并非所有時間尺度的空中交通流時間序列都具有非線性，本文5 min時間尺度的空中交通流時間序列是非線性的，可以進一步分析其分形特征；

（2）本文空中交通流時間序列具有自相似、長相關(guān)、無標度和多重分形等典型的分形特征；其關(guān)聯(lián)維度為6.89，需要7個參數(shù)才能清晰描述該空中交通流；

（3）后續(xù)將針對其他空中交通流時間序列開展非線性分形特征研究，以探索分形在空中交通流系統(tǒng)中的普適性.

致謝：中國民航大學省部級科研機構(gòu)開放基金（KGJD201502）的資助.