王 飛

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

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

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

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

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

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

1 空中交通流時(shí)間序列

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

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

2 空中交通流的非線性檢驗(yàn)

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

（1）零假設(shè)

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

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

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

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

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

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

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

（3）檢驗(yàn)統(tǒng)計(jì)量

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

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

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

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

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

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

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

3.1 自相似性特征分析

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

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

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

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

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

3.2.1 Hurst指數(shù)

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

采用重標(biāo)級差法[14]先計(jì)算出重標(biāo)極差序列及其對應(yīng)的時(shí)間增量序列，繪制曲線圖，并運(yùn)用回歸分析計(jì)算出斜率，即為H.

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

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

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

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

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

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

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

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

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

3.3 無標(biāo)度性特征分析

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

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

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

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

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

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

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

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

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

以m= 14，τ1= 3為例，應(yīng)用GP算法繪制的雙對數(shù)關(guān)聯(lián)積分曲線和對應(yīng)的差分如圖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) 曲線計(jì)算得到區(qū)間為 l nr∈ [1.631 4,2.1910] ，從lnr-lnC′′(m,r)曲線可以看出，在0附近直線的lnr∈[1.949 8，2.260 0].綜合考慮兩個(gè)曲線圖，得到本時(shí)間序列的無標(biāo)度區(qū)間為 [1.949 8,2.191 0] .不同嵌入維度下的無標(biāo)度區(qū)間如表2所示.

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

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

3.4 多重分形特征分析

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

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

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

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

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

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

時(shí)間序列計(jì)算得到的多重分形譜如圖6所示.

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

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

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

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

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

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

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

本文采用GP算法來計(jì)算空中交通流時(shí)間序列的關(guān)聯(lián)維數(shù)D2.按照3.3節(jié)所述方法，首先設(shè)置m=2，繪制 lnr- l nC(r) 圖，在無標(biāo)度區(qū)間內(nèi)用最小二乘擬合得到斜率，即為此時(shí)的關(guān)聯(lián)維數(shù)估計(jì)D2(2) .然后，逐步增加嵌入維數(shù)，直到關(guān)聯(lián)維數(shù)不再隨著m的增加而增加.然后繪制m-D2(m) 圖來獲得關(guān)聯(lián)維數(shù).本文時(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可以看出，當(dāng)m= 10時(shí)，關(guān)聯(lián)維數(shù)D2達(dá)到最大值，此后D2隨m的增加而不再增加.因此，m= 10對應(yīng)的關(guān)聯(lián)維度是該時(shí)間序列的關(guān)聯(lián)維度，即D2=6.89 ，說明需要7個(gè)變量才能清晰描述該空中交通流.

5 結(jié) 論

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

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

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

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