空中交通流非线性分形特征

王 飞

（中国民航大学空管学院，天津 300300）

空中交通流是大量航空器沿航路航线相继飞行形成的航班流，是空中交通流量管理的主要对象.对空中交通流时空特性的剖析是分析空中交通流系统内在特性、揭示空中交通流系统动态演化规律的重要基础[1]，在空中交通流建模、预测和复杂性分析等领域具有广泛应用前景.

当前，空中交通流的研究主要从3方面开展：基于数学模型的分析[2-3]，是对系统的高度抽象，难以真实反映空中交通流运行实际；基于系统仿真的研究[4-5]，考虑了不同研究对象交互作用对空中交通流的影响，但仍难以模拟不确定性因素；基于实测数据的数据挖掘分析[6-7]，可以直接从实际运行数据中挖掘空中交通流的特征，更为真实地反映空中交通流内在规律性.

随着研究的逐步深入，国内外学者逐渐意识到空中交通流系统类似于地面交通流系统，也是一个复杂的非线性动力系统，并开始应用非线性动力学理论研究空中交通流.虽然非线性理论在道路交通方面得到广泛研究和应用，但在空中交通领域的研究刚刚起步，研究成果较少.李善梅等提出应用小数据量法和小波去噪理论计算最大Lyapunov指数，并对飞行冲突时间序列的混沌特性进行了分析[8].Cong等以扇区交通流量为对象建立时间序列数据，并对其混沌特性进行了深入研究[9].郑旭芳等从混沌与分形角度对交通流量时间序列的非线性特征进行了对比分析，并分析了时间尺度的影响[10].杨阳等利用混沌特性研究了扇区流量的短期预测方法[11].

传统的空中交通流建模与预测方法均是基于统计意义上的分布函数，鉴于空中交通流系统已被证实具有混沌、分形等非线性特性，传统方法过于理想化，难以准确把握空中交通流的内在特征和演化规律.从现有研究可以看出，对空中交通流的混沌识别与预测研究较多，而对于分形特征的研究也仅停留在计算关联维数这一表象，对其内在的分形特征都尚未开展研究.基于实际运行数据构建的时间序列是研究非线性系统的有效手段.研究空中交通流时间序列的分形特征，一方面为现有空中交通流建模方法是否合理提供判据，另一方面为应用分形理论实现空中交通流短期预测奠定基础.

本文以基于实际运行数据的管制扇区流量时间序列为研究对象，首先对时间序列进行非线性检验，然后系统性分析时间序列的分形特征，最后计算了时间序列的关联维数.

1 空中交通流时间序列

本文采用文献[11]中的三亚飞行情报区2016年5月某天零时开始连续12 d的流量统计数据.为全面分析空中交通流的非线性分形特征，以三亚管制区1号扇区为对象，将数据处理为时间尺度30、15、10 min和5 min的4个交通流时间序列，流量时间序列如图1所示.

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

2 空中交通流的非线性检验

检验空中交通流时间序列是否具备非线性是应用非线性理论的必要条件.替代数据法是检验时间序列非线性的常用检验方法，包括零假设、替代数据生成和计算检验统计量3个部分组成.该方法基本思路是：假设观测的时间序列是线性的，在保持原有线性关系下根据原始分布的性质（如均值、方差等）生成一组替代数据，然后计算原始数据和替代数据的检验统计量.若两者检验统计量有显著差异，则拒绝零假设，说明原始时间序列不太可能是从假设的线性系统中产生，即原始数据中存在确定的非线性成分.

（1）零假设

观测数据由具有原始数据的均值和方差的线性相关高斯过程产生[12].如果拒绝该零假设，则待检验时间序列包含非线性成分.

（2）生成替代数据

保留原时间序列功率谱增幅，通过改变功率谱中各频率成份的相位来生成替代数据.具体步骤如下：

步骤1将原始时间序列A进行傅里叶变换，生成新的序列X= {Xi}，i=1,2,···,n，n为原始序列长度.

步骤2随机生成新的功率谱相位序列F={fi} ，fi∈ [0,2π]；

步骤3应用随机相位按照式（1）改造X，构建新的时间序列Y={Yi} ；

步骤4对时间序列Y进行逆傅里叶变换，得到新的时间序列，作为替代数据.替代数据和原有数据具有相同的功率谱和自相关函数.

（3）检验统计量

检验统计量的选择将会直接影响判断结果.本文选取KS检验作为检验统计量，该方法是一种有效、稳定的非线性检验方法[13].

应用替代数据方法对4个交通流时间序列进行非线性检验.由于功率谱相位序列是随机产生，每次产生的功率谱相位序列都不一样，结合蒙特卡洛思想，每个时间序列均进行了1 000次检验（显著性水平为0.05），具体检验结果见表1.

从表1可看出，不同时间尺度的时间序列非线性检验结果是不一样的，即随着时间尺度的增大，具备非线性的概率越小.需要说明的是，小概率非线性并不代表一定不是非线性，还应当结合其它指标（例如关联维度）联合检验.而时间尺度小于5 min的时间序列是否依然具备非线性，后续将另文研究.

表 1 时间序列非线性检验结果Tab.1 Nonlinear test results of time series

本文尺度为5 min的空中交通流时间序列可以确定具有非线性成分，但是否具备分形特征还需进一步分析.

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

分形是指由各部分组成的形态，每个部分以某种方式与整体相似.交通流分形现象从整体来说包括自组织、自相似和无标度，从局部来说具有多重分形.根据这4类现象，归纳出空中交通流分形特征包括自相似、长相关、无标度和多重分形.

3.1 自相似性特征分析

空中交通流的自相关性并不是严格意义上的完全重合，而是在不同的空间或时间尺度的统计意义上的自相似.为了避免噪声干扰，本文采用具有优良的时域和频域局域化能力的小波变换，对原始时间序列进行小波多尺度分解，分解尺度为3，提取其中低频系数结果，如图2所示.

图 2 小波分解结果Fig.2 Results of wavelet decomposition

从图2中可以看出，在不同日期里空中交通流呈现出类似“周期”的变化，显示出较强的自相似性.当分解尺度为1、2时，可得到相同结论.

3.2 长相关性特征分析

长相关性也称为长程相关性或记忆性，是产生分形的根源.长相关性反映的是过去的状态对当前和未来的状态的影响是长期的，这也是能够根据过去状态进行交通流预测的一个主要原因.R/S方法是分析时间序列长相关性的有效方法，通过计算Hurst指数H来判断时间序列是否具有长程相关性.

3.2.1 Hurst指数

Hurst指数反映了时间序列A均值的累计离差随时间变化的范围.

采用重标级差法[14]先计算出重标极差序列及其对应的时间增量序列，绘制曲线图，并运用回归分析计算出斜率，即为H.

根据统计特性：H= 0.5表示观测序列是变化无规律，即随机的；H＜ 0.5表示观测序列是反持续的时间序列，即前一时点上升（下降），下一时点将有下降（上升）的趋势；H＞ 0.5表示所观测序列是持续性序列，即前一时点上升（下降），在下一阶段要继续这种正（负）的变化趋势.由于Hurst指数能反映序列的持续性或随机性，在实际中它得到广泛的应用.

根据时间序列的时间跨度不同，Hurst指数分为全局Hurst指数和局部Hurst指数.全局Hurst指数表示整个时间序列长相关性，在空中交通流中表示交通流保持之前趋势的强弱性.局部Hurst指数反映局部交通流的变化情况.两者的计算方法是相同的.

该时间序列的全局Hurst指数大于0.5，该时间序列并不是随机的，并且具有正相关性，即每一个观测值都保存着之前所有事情的记忆，未来的情况也与现在的情况有巨大关联性，说明该时间序列具有长相关性.

图 3 Hurst指数拟合曲线Fig.3 Fitting curve of Hurst exponent

为了检验R/S分析的可靠性，将原始时间序列的数值顺序进行随机打乱，形成新的时间序列.若原始时间序列是独立的，则新的时间序列的Hurst指数应基本保持不变，反映原始时间序列没有长相关性[14].按照重标级差法经过1 000次计算得到新时间序列的Hurst指数在（0.55，0.62）之间，相较于原始时间序列有大幅下降，说明原始时间序列不是独立随机序列，R/S分析结果是可靠的.

3.2.3 局部Hurst指数分析

为反映交通流随时间的变化规律，将原始时间序列按照288个数据（1 d）为一个子区间，以1个数据为间隔向后推，得到一系列时间序列子区间，计算各子时间序列的Hurst指数，如图4所示.

图 4 局部Hurst指数Fig.4 Local Hurst exponents

局部Hurst指数均大于0.5，表现出较强的长相关性，说明交通流能够保持之前的变化趋势.每一天周期内均会出现几个局部低值，表现出长相关性减弱、随机性增强，是交通流趋势最可能出现反转的点，反映的是交通流处于拥堵阶段.

3.3 无标度性特征分析

分形往往在一个最小标度和最大标度构成的区间内才存在分形规律，这个区间就是无标度区间.一旦逾越无标度区间，自相似性便不复存在，系统也就没有分形规律了.因此，无标度区间的确定一直是分形研究的热点问题.

无标度区间通常是与关联维数计算紧密联系的，GP (Grassbegr-Procaccia)算法是计算关联维数的常用方法，计算步骤如下：

步骤1对原始时间序列A采用时间差法重构相空间，形成新的时间序列X2为

式中：A= {A(i) }；m为嵌入维数，可人为设置数值，仅针对该嵌入维数进行相空间重构； τ1为时间序列的延迟时间，根据自相关系数法计算得到.

步骤2计算相空间中的每一个向量点A(i) 与其它向量点的距离，然后统计在以A(i) 为中心、r为半径的圆球中的向量点的个数，从而得到关联积分C(m,r)为

式中： ‖A(i)−A(j)‖ 表示相空间中两个向量点之间的距离；H(·) 表示Heaviside函数。

步骤3绘制lnr- l nC(m,r) 曲线图，寻找线性度最好的区间 [lnr1,lnr2] ，即无标度区间，lnr1、lnr2分别为该区间的端点.

目前，研究者大多采用视觉法识别无标度区间，简单易行，但主观性和随意性较强.文献[15]应用一阶差分数据k-means方法识别出无标度区间，但该方法对具有非线性动力方程的系统，如Lorenz等具有较好的识别效果，但对于诸如交通流时间序列的离散数据的识别效果不尽如人意.文献[16]提出应用双对数曲线点的二阶差分数据，识别出零值附近波动数据，再剔除粗大误差，用统计方法识别出无标度区间.本文采用文献[16]的方法，对空中交通流时间序列的 l nC(m,r) 差分数据进行分析.

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

以m= 14，τ1= 3为例，应用GP算法绘制的双对数关联积分曲线和对应的差分如图5所示.

图 5 双对数关联积分和二阶差分曲线Fig.5 Curves of double-logarithmic correlation integral and second-order difference

从图5可以看出，由lnr- l nC′(m,r) 曲线计算得到区间为 l nr∈ [1.631 4,2.1910] ，从lnr-lnC′′(m,r)曲线可以看出，在0附近直线的lnr∈[1.949 8，2.260 0].综合考虑两个曲线图，得到本时间序列的无标度区间为 [1.949 8,2.191 0] .不同嵌入维度下的无标度区间如表2所示.

表 2 不同嵌入维度下的无标度区间Tab.2 Non-scale range under different embedding dimensions

从表2可以看出，不同嵌入维数对应的无标度区间是有区别的，但并未发现无标度区间与嵌入维数具有明显规律性.

3.4 多重分形特征分析

若时间序列仅具有某种单一自相似特征，这样的分形称为简单分形，只需选择一个标度就可以刻画.多重分形则描述多标度复合分形的复杂体系.许多分形具有更复杂的尺度关系，需要一组参数去描述复杂系统特性，这就是多重分形.空中交通流时间序列的多重分形特征可以应用配分函数和多重分形谱来进行分析.具体步骤如下：

步骤1将原始时间序列A等分为N个长度为∆t的小区间，定义流量增量为

步骤2将配分函数定义为流量增量的q阶矩的加权和，表示为

步骤3针对每一个q值，绘制 l g(∆t) -lg(Sq(∆t))双对数曲线图，并利用最小二乘法进行线性拟合，其斜率为质量指数 τq，如式（8）所示.

步骤4将τq进行Legendre转换，即可得到多重分形谱为

式中： α为奇异性指数，用来描述各个子区间之间不同的奇异强度，f( α )为具有相同的子集的分数维度，通过分形谱函数可以刻画不同α值对应的分形特征.

时间序列计算得到的多重分形谱如图6所示.

多重分形谱构成的几何图形一般表现为钩状或者钟状，包含 6个参数.其中，参数 αmax和 αmin分别表示最大概率子集和最小概率子集的奇异指数；参数 ∆α=αmax−αmin表示多重分形谱的宽度，可以定量表征时间序列波动的强度；参数f(αmin) 和f(αmax) 分别表示最小概率子集和最大概率子集的分形维数；参数 ∆f=f(αmax)−f(αmin) 表示最大概率子集和最小概率子集的分形差异，反映时间序列的变化趋势.当∆f＞0时表明时间序列数据向高数值方向发展的可能性较大，当 ∆f＜0 时则表明时间序列数据向低数值方向发展的可能性较大.

从图6可以看出，原始时间序列的 α与f(α) 的关系呈现“钟形”，说明存在多重分形特征[16].∆α为0.042 6，说明交通流量并不是均匀分布，具有一定波动性.∆f为0.069 2，说明交通流量处于高流量的机会比处于低流量的机会大.

图 6 时间序列的多重分形谱Fig.6 Multifractal spectrum for time series

通常有两种类型的多重分形，一种是由时间序列的概率密度函数的“肥尾”分布引起的，另一种是由时间序列的长相关性引起的.判断一个时间序列产生多重分形的本质根源，可通过随机打乱原始时间序列元素顺序来分析.对于第一种情况，随机打乱顺序也无法消除序列的多重分形性质；对于第二种情况，随机打乱顺序将导致时间序列的长相关性被破坏，往往会呈现非多重分形性质.如果这两种情况同时存在，则打乱顺序的时间序列依旧呈现多重分形性质，但其多重分形性相较于原始序列较弱.

从图6可以看出，打乱顺序后的时间序列仍然具有“钟形”的多重分形特征，并且 ∆α和 ∆f都有所减弱，因此原始交通流的多重分形性质是由长相关性和概率密度的“肥尾”分布共同决定的.

4 空中交通流的分形维数计算

分形维数是定量刻画分形特征的参数，是分形的最基本特征.由于研究对象的不同，分形维数有多种定义，常用的包括Hausdorff维数、相似维数、信息维数、Lyapunov维数和关联维数等.虽然这些维数定义不同，但只要计算出的维数大于1，即可判定研究对象具有分形特征.

本文采用GP算法来计算空中交通流时间序列的关联维数D2.按照3.3节所述方法，首先设置m=2，绘制 lnr- l nC(r) 图，在无标度区间内用最小二乘拟合得到斜率，即为此时的关联维数估计D2(2) .然后，逐步增加嵌入维数，直到关联维数不再随着m的增加而增加.然后绘制m-D2(m) 图来获得关联维数.本文时间序列的lnr- l nC(r) 关系见图7，m-D2(m) 关系见图8.

图 7 双对数关联积分曲线Fig.7 Curves of double-logarithmic correlation integral

图 8 关联维数随嵌入维数的变化曲线Fig.8 Curve of the correlation dimensionD2with the increasing embedding dimensionm

从图8可以看出，当m= 10时，关联维数D2达到最大值，此后D2随m的增加而不再增加.因此，m= 10对应的关联维度是该时间序列的关联维度，即D2=6.89 ，说明需要7个变量才能清晰描述该空中交通流.

5 结 论

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

（2）本文空中交通流时间序列具有自相似、长相关、无标度和多重分形等典型的分形特征；其关联维度为6.89，需要7个参数才能清晰描述该空中交通流；

（3）后续将针对其他空中交通流时间序列开展非线性分形特征研究，以探索分形在空中交通流系统中的普适性.

致谢：中国民航大学省部级科研机构开放基金（KGJD201502）的资助.