基于3种可视图的进场航班流量波动特性适应性评估*

张 勰 肖恩媛 刘宏志 赵嶷飞 王梦琦

（1. 中国民航大学空中交通管理学院 天津 300300；2. 中国民航科学技术研究院民航发展规划研究院 北京 100028）

0 引 言

掌握空中交通流量涨落波动特性对于优化空域资源配置、提升运行效率、保障飞行安全方面发挥着先导性、基础性、关键性作用。而研究3 种可视图（可视图、水平可视图、有限穿越可视图）建网方法的适应能力和演化能力对于掌握航班流量波动规律和刻画流量波动动态演变模式是十分重要的一环，这有助于深刻理解整个空域交通系统的运作原理与演化机制，为航班流量演变建模仿真及预测提供理论基础与分析方法。

针对可视图（visibility graph，VG）对时间序列建网方面的研究，2017年，邢雪等[1]结合可视图理论来构建交通流时间序列的关联网络，从复杂网络角度实现交通流时间序列的特性分析。2018 年，刘宏志等[2]利用可视图和水平视图（horizontal visibility graph,HVG）刻画空中交通流波动。在此基础上，有限穿越可视图方法[3]（limited penetrable visibility graph，LPVG）也随之发展起来，并且显示出独特的优势，近些年被广泛应用于各个领域。2012年，周婷婷等[3]采用3 种可视图网络考察时间序列建网方法及其对周期、分形、混沌以及添加噪声信号后的适应能力。2016 年，Wang 等[4]采用了有限穿越可视图和相空间方法将EEG 序列映射成网络，研究AD 大脑潜在的混沌系统动态学。Ming等[5]对有限穿越可视图进行抗噪测试。2020年，Ren等[6]将有限穿越可视图与序模体方法结合，结果证明有限穿越可视图序模体方法有良好的鲁棒性。

对于空中交通流量时间序列的动态演变和波动特征方面，2018年，张勰等[7]将复杂网络理论引入到时间序列的应用当中，使用motif 理论提取波动模式，结果发现时间序列中存在显著的波动模式以及这些波动模式之间的转换循环。随后，2020年，围绕机场空中交通流量时间序列的复杂性问题，Liu 等[8]在引入了多个粗粒化的过程中，提出了1 种改进的MMPE 方法，有助于深入了解空中交通系统的演化机制，从而揭示时间序列多尺度复杂性的剧烈波动。2020 年，针对进离场空中交通流量时间序列，刘宏志[9]从交通流波动动态演化、波动多尺度分形特征和多尺度复杂性的角度探讨流量短期及长期的演变规律，实验结果证明该研究能够为流量波动的仿真建模、预测及管理提供科学的理论基础和分析工具。

基于分析航路网络特征方面，2020 年，王红勇等[10]从复杂网络的角度去构建网络模型。2021 年，隋东等[11]引入复杂网络修复理论和交通流分配理论，提出了1 种航路网络修复优化策略。针对机场终端区空域内进离场的空中交通流，2014—2016年，有关学者分别采用建模仿真工具[12-13]、跟驰模型[14]、元胞传输模型[15]和实证分析[16]等方法研究终端区空中交通流量、密度和速度参数之间的动态演化关系。随后，针对日前繁忙机场终端区常见的进离场航线结构，2020年，张洪海等[17]提出相应TBO模式下的进场交通流优化模型，并为此构建仿真运行环境。

围绕交通流时间序列内在特征的研究，2018年，王超等[18]研究了1 种基于改进加权一阶局域法的混沌交通流量时间序列预测模型。为了研究多尺度分析的时间序列与预测策略，2004 年，Weinreich 等[19]研究预测对于时间序列在不同时间尺度上所表现出来的依赖性。针对航空发动机时间序列预测中嵌入维数难于有效选取的问题，2010 年，张弦等[20]提出1 种基于嵌入维数自适应最小二乘法支持向量机的预测方法。针对交通系统内在特征，2020年，为了对空中交通系统自组织临界特性进行辨别和应用，王兴隆等[21]采用幂律特性直线的斜率判断航班延误、航班取消的管理水平；依据Hurst 指数大小，分析系统的航班延误、航班取消在自组织临界状态下的相关性与演化规律，结果表明空中交通系统呈现耗散结构。

综上所述，针对可视图对时间序列建网方面的研究，仅基于单一可视图进行分析，缺乏对多种可视图多时间粒度角度的分析。对于空中交通流量时间序列的动态演变和波动特征方面，仅针对传统可视图，缺乏对改进可视图方法的引进和可视图理论的创新。基于分析航路网络特征方面，较少从波动演化角度结合复杂网络理论，缺乏对时间序列复杂网络波动演化的研究探索。围绕交通流时间序列内在特征的研究，多是针对时间序列分形、混沌等方面的改进与创新，对于时间序列波动特性的研究目前还较少。

为从时间维度考察3种可视图对于空中交通的动态演化适应能力，本文面向3个时间粒度，采用可视图建网方法将进场航班流量时间序列映射成复杂网络，分别从网络整体统计特性和网络局部结构模式2 个视角开展研究，考察多粒度多维度不同可视图对于流量时间序列的波动特性适应性以及波动模式的演化能力。

1 时间序列可视图复杂网络模型

1.1 时间序列

按照中国民航空管流量管理运行规则[22]扇区流量统计方法，本文分别基于2018年12月1日当日和及2018 年12 月1 日—31 日整月天津机场进场航班流量数据（06：30—23：50），天津机场进场航班流量数据，分别采用5，10，20 min时间间隔统计3个时间粒度的航班流量时间序列，3个粒度的时间序列长度分别为208，104，52，时间粒度5 min，见图1。航班当日流量时间序列见式（1）。

图1 天津滨海国际机场进场航班流量时间序列Fig.1 The arrival flight flow volume time series of ZBTJ

式中：xi为第i个时段内终端区进场流量（架次/min）；M为时间序列的长度（节点数量）。

1.2 3种可视图复杂网络

1.2.1 可视图复杂网络

Lacasa 等[23]提出了可视图建网方法（visibility graph，VG）。可视图将复杂网络与时间序列形成1种相互对应的关系，这使得从复杂网络分析得到的结论可以用来分析原始时间序列的性质。

可视图准则定义如下[8]：如果在ti和tj之间的任何观测值都满足如下标准，则在ti和tj得到的观测值xti和xtj彼此可见，并且可用可视线连接成为VG中的连接点对，定义见式（2）。

式中：xti为在时间周期t内第i个时段内终端区进场流量数值（架次/min）。

可视图利用柱状体来表示流量节点数值的大小，用柱状图表示时间周期内节点流量的起伏变化。若2个柱体顶端相互可见且无任何阻挡，则认为2个节点在网络中存在可视线。值得指出的是，节点与自己构成的自环不能当作可视线，可视线不得穿越任何柱体。例举5 min时间粒度航班流量时间序列见图2。

图2 空中交通流量时间序列映射得到的可视图Fig.2 Visibility graphs mapped from flight flow volume time series

1.2.2 水平可视图复杂网络

Luque 等[24]在可视图方法基础上提出了水平可视图（horizontal visibility graph，HVG），见图3，其节点和可视线的定义与可视图类似，并在此基础上要求2 个柱体之间的可视线必须为水平状态，期间不能有阻挡。值得注意的是，同一时间周期内时间序列的VG可视线能够包含HVG可视线，即HVG网络总是VG网络的1个子网络[3]。

其水平可视性准则[24]可以定义见式（3）。

式中：xti为在时间周期t内第i个时段内终端区进场流量数值。任何在ti和tj之间的xtk并没有穿越xti和xtj的直线，则该水平直线就是HVG中的1条有效可视线。例举5 min时间粒度航班流量时间序列，见图3。

图3 空中交通流量时间序列映射得到的水平可视图Fig.3 Horizontal visibility graphs mapped from flight flow volume time series

1.2.3 有限穿越可视图复杂网络

有限穿越可视图（limited penetrable visibility graph，LPVG）在可视图的基础上，定义了穿越距离N。若时间序列中2 节点连线截断的数目（即被其他节点截断的个数）满足d≤N，则说明2 个节点可连通[6]。这表明，当2个节点之间的连接视线被截断的次数不超过有限次数N时，2 个节点之间是存在有效可穿越可视线的。例举N=1时，5min时间粒度航班流量时间序列见图4。

图4 空中交通流量时间序列映射得到的有限穿越可视图Fig.4 Limited penetrable visibility graphs mapped from flight flow volume time series

参数N的设置能够有效地减少噪声对时间序列的影响[3]。参数N越大，允许穿越的节点数量就越多，时间序列节点之间的整体关联性就会越强，但是N选择过大会相对地弱化相邻节点的关联性，即N的选择应在合理提升抗噪能力的同时最大程度上保留可视图局部结构的原貌。

2 基于流量时间序列整体波动统计特性的3种可视图适应性评估

可视图作为可直接将时间序列关系映射成复杂网络的工具，能够不同程度地反映出时间序列的波动特性。为了考察3种可视图表征空中交通流时间序列复杂结构特征的能力及对其动态波动的适应性，采用矩阵从属阵和k-core 等方法对不同可视图方法流量时间序列整体波动整体统计特性适应症进行评估。

2.1 3种可视图映射的网络结构特征

2.1.1 网络结构特征实验分析

利用可视图可视线无方向性特征构造对等性可视化矩阵图，对3种可视图及其不同N取值LPVG进行块建模，利用块建模网络矩阵表征网络结构，有效地分析网络的结构特征。全部分为上下2个部分，上部分方块表征节点之间的连接关系，且沿对角线呈对称分布状态，下部分对应不同时间粒度时间序列。

连接关系矩阵-时间序列图例举5 min时间序列见图5～6。图中矩阵对角线上的矩阵点不代表自环，全图上下2 个部分利用不同灰度的从属阵作为子单位，3 种时间粒度从属阵个数分别为9，8，7，实线指向的节点表示从属阵的区分节点，区分节点的计算见式（4）。

图5 3种可视图的网络矩阵结构图Fig.5 Network matrix structure diagrams of three visibility graphs

图6 不同N取值的LPVG网络矩结构图Fig.6 Limited penetrable visibility graph network matrix structure diagrams with different N values

式中：i，j为节点数；xi为第i个节点（时段）内流量值（架次/min）；Ci为节点流量值大于左右节点流量值的最小次数；M为时间序列节点的总数（总长度），M=208，104，52；其中15表示经验值。

实线划分从属阵且指向从属阵中最高流量值节点，刻画了空中交通流量时间序列从属阵的流量峰值，全部从属阵节点连接关系动态变化表征时间序列整体波动动态特征，通过计算从属阵实际可视线比率来评估可视线密度的变化。见图5～6，时间序列X坐标排序方式表示从属阵可视线比率从大到小的排列，其计算见式（5）。

式中：l为第i个节点与n个（i＜n≤N）节点连接线段的数量；j为第j个从属阵；N为从属阵内部总节点数量。

对从属阵内所有节点的可视线比率按照式（5）进行计算，得到3种可视图可视线比率的统计结果，见表1，同时对不同N取值LPVG见表2。

表1 3 种可视图从属阵可视线比率Tab.1 Subordinate array visibility lines ratios of three visibility graphs

表2 不同N 取值LPVG 从属阵可视线比率Tab.2 Limited penetrable visibility graph subordinate array visibility lines ratios with different N values

3 种时间粒度可视线比率LPVG（N=1）比VG高出54.9%，37.0%，27.4%，VG与LPVG（N=1）可视线比率变化呈现正比关系。二者均能够表征节点之间的连接属性，具备刻画空中交通流时间序列复杂结构特征的能力。由图5 和表1 可见：HVG 从属阵可视线密度极低，因为HVG要求满足相同数值节点相连的连接条件，导致存在可视线间断的情形，不符合采用空中交通流量时间序列节点连接属性表征波动涨幅特性的基本要求。

在LPVG连接关系中，N值越大，允许穿越的节点数量越多，相应的可视线数量增多和可视线比率变大，但这种变化会削弱相邻节点原始的关联属性，因此定义细节损失率见式（6）。

式中：j为第j个从属阵；ρj为第j个从属阵的可视线密度；ρ0为当N=1 时相应的第j个从属阵的可视线比率。

利用式（5）～（6）绘制LPVG连接关系下不同N取值从属阵的细节损失率见图7，其中实线表示数值为0.5 的细节损失率，N=1 可视线比率为N=2～6的对比参数，所以N=1细节损失率为0。

图7 不同N 取值的LPVG从属阵细节损失率Fig.7 Limited penetrable visibility graph subordinate array detail loss rates with different N values

N=1 时数值所占节点数量（208，104，52）百分比分别为0.48%，0.96%，1.92%。利用LPVG 连接关系N=1 可视线比率数值作为对比参量，N值的选择需要确保细节损失率至少在0.5以内，保证相邻节点的信息失真率在此范围内变化达到较为真实反映节点的连接属性的条件。

2.1.2 网络结构特征实验分析结果

综上所述，其计算结果如下。

1）时间粒度为5 min时（208节点）N=2 和N=3的从属阵细节损失率均值分别为0.270 和0.461，N值数值所占节点数量百分比为0.962%和1.442%；

2）当时间粒度为10 min时（104节点），N=2～6细节损失率均值范围在0.145～0.337，N=2～6 数值所占节点数量百分比为1.923%～5.770%；

3）当时间粒度为20 min 时（52 节点），N=2～6细节损失率均值范围在0.072～0.130，N=2～6 数值所占节点数量百分比为3.846%～11.538%。

时间粒度设置为10，20 min 能够保证N=2～6情况下细节损失率都小于0.5，这是因为时间粒度选择过大导致统计后的流量时间序列缩短，原始时间序列整体节点连接特征被削弱，所以N=2～6 不同N取值情况下，被削弱整体统计特征的流量时间序列不能够表达出更多的细节特性，以至于从属阵细节损失率能够保持在0.5以内。

基于上述对3 种可视图网络及其不同N 取值LPVG网络整体结构的分析，结论如下。

1）VG与LPVG能够表征节点之间的连接属性，能够具备刻画空中交通流时间序列复杂结构特征的能力，HVG不符合采用空中交通流量时间序列节点连接属性表征波动涨幅特性的基本要求。

2）对于LPVG而言，N选取数值与节点数量占比在0.48%～1.442%区间既能够保证细节损失率在0.5范围内，又能确保在不损失时间序列相邻节点连接信息属性的前提下，最大程度表征该时间序列的动态波动特征，针对该交通流时间序列研究中推荐使用208（时间粒度为5 min）节点N=1～3 的LPVG网络。

2.2 k-core 评估3 种可视图对于流量波动特征的适应性

2.2.1 流量波动特征实验分析

在可视图中，节点的流量值差异越大，涨落起伏越剧烈，可视线越密集。因此，可视线的密度反映了流量时间序列的波动强度。为了刻画流量时间序列的整体波动强度，采用k-core算法[25]对所得到的网络进行节点聚类分析。该算法利用节点的度刻画单个节点的连接关系，度越大，节点的可视线越多，节点阶数越大，因此可用k阶核量化流量的波动强度。

k阶核定义：对于图G=（V，L）（V是节点，L是连边），定义n=|V|，m=|L|，W 是节点集合，则定义子图H=（W，L|W）当且仅当∀Degree_V（V∈H）≥k,且H是1 个具有这样属性的最大子网络时，称H 为节点集合W诱导的1个k阶核。各节点所属核的阶数为其所能参与的最高的k阶。该方法能够得到节点度不小于k的最大子网络。

基于上节结论，后续内容不涉及水平可视图，同时LPVG连接关系时使用208（时间粒度为5 min）节点N=1～3。采用k-core 波动强度算法考察VG 与LPVG对波动动态的适应性。2种可视图采用k-core算法分类的可视化网络见图8，表3 为k-core 统计特性。

表3 k-core 统计特性Tab.3 The statistic characteristics of k-core

图8 3种可视图k-core网络图Fig.8 Three visibility graphs k-core network graphs

VG，LPVG（N=1），LPVG（N=2），LPVG（N=3）网络k-core 算法中core 分类数量分别为8，13，17，20，同时最高频数k-core 分别为5 阶核（44次）、10 阶核（53 次）、16 阶核（44 次）和20 阶核（58次）。允许穿越节点（N）数量增加1，最高频数k阶核分别能够增长100%（5阶核）、60%（6阶核）和25%（4 阶核）。可视图可穿越距离（N）越大，节点k阶核种类越多，最高频次k阶核阶数越高，刻画流量时间序列的波动的动态特征就越细致。

依据以下计算步骤，进一步考察2 种可视图时对时间序列整体波动态势的适应性，见图9～10。

1）将2018 年12 月1 日—31 日相同时段的时间序列采用VG与LPVG连接关系，利用k-core算法得出相应阶数。

2）1个月节点阶数见式（7）。

式中：Ki为1 个月的求和k阶核；i为时间节点，1 ≤i≤M；M为时间序列长度208；j为日期（1—31日）。

3）1个月的流量值见式（8）。

式中：Xi为1 个月流量平均值；架次/min；i为时间节点，1 ≤i≤M；M为3 个时间序列长度208；j为日期（1—31日）。

4）绘制计算结果，并且用不同灰度标注。再利用方差评估1个月k-core浮动变化波动适应性评估值。

图9 为时间序列1 个月k-core 整体波动趋势，4种情形都能够通过k-core算法表现时间序列整体波动趋势，k-core 在N=3 情形下呈现最大值。在此基础上，考察4种情形1个月k-core浮动变化波动适应性评估值见图10。

图9 时间序列整体波动态势Fig.9 Overall fluctuation trend of the time series

图10 为4 种情形下1 个月的k-core 浮动变化波动适应性评估值。VG 均值为2.665，LPVG 关系N=1 为4.810、N=2 为6.973、N=3 为9.883。随着N值的增大，波动适应性评估值提高80.5%（2.145）、45.0%（2.163）和41.7%（2.910）。

图10 方差评估波动适应性Fig.10 Variance assessment fluctuation adaptability

2.2.2 流量波动特征实验分析结果

综上所述，基于上述对3种可视图网络及其不同N取值LPVG流量波动特性实验的分析，结论如下。

1）可视图可穿越距离（N=1～3）越大，k-core浮动变化波动适应性评估值越高，对于时间序列动态波动的适应性能力越强，越能够有效精准刻画时间序列动态波动特性。

2）基于上述对VG及其N=1～3取值LPVG网络波动适应能力的分析可知，VG 与LPVG（N=1～3）都能够表征空中交通流量时间序列的波动强度，刻画时间序列的动态波动特征。

3）VG 与LPVG（N=1～3）网络对时间序列波动适应性评估值为2.665，4.810，6.973，9.883。

3 2 种可视图网络对流量局部结构波动动态演化能力评估

为了揭示VG 与LPVG（N=1～3）网络对于流量时间序列局部动态演化能力，采用序模体方法，从网络局部结构角度分析研究二者对于交通流的动态演化能力。

Shen-Orr 等[26]提出了模体（motif）概念，Iacovacci等[27]受到启发提出了序模体概念，在此基础上Ren等[6]提出了有限穿越可视图序模体。在可视图节点按时间顺序排列的情况下，构造1个滑动窗口，依序滑动，在每个窗口中得到的网络子结构（即子网络）就是1个有可视图序模体。

针对VG与LPVG（N=1～3）映射得到的4个可视图网络，此时4 种网络能够存在的最小序模体节点数量为3，4，5，6。因为当N=1 时，时间序列2 个节点之间只相隔1 个节点的可视线必定存在，这表明LPVG（N=1）网络要刻画波动模式的差异需要4个节点以上长度的序模体，N=2～3 以此类推。

3.1 序模体类型

在得到的VG 网络中3，4，5 节点序模体出现的类型总数2种、5种和22种，LPVG（N=1）网络4，5，6 节 点 序 模 体 分 别 有2 种、6 种 和32 种，LPVG（N=2）网络5，6，7 节点序模体分别有2 种、6 种和34种，LPVG（N=3）网络6，7，8节点序模体分别有2种、6 种和34 种。例举LPVG（N=1）的所有序模体类型见图11。

图11 LPVG（N=1）流量时间序列中出现的序模体类型Fig.11 Sequential motif types of limited penetrable visibility graph（N=1）presented in the flow time series

3.2 4 种情形下表征流量波动模式动态演化的能力评估

将流量时间序列中每1个时间点上的空中交通流量视为1个节点，而节点之间的连接关系由2种可视图中的可视线准则确定，分别按照不同节点长度宽度设置滑动窗口，移动步长为1，当滑动窗口移动经过所有的节点，则得到时间序列中出现的所有序模体。将不同序模体类型作为波动模式动态转化图中的不同节点，按照滑动窗口的移动所体现的类型转移关系，绘制不同序模体类型之间的转移方向，即波动模式动态转化图中相应节点之间的有向边，当滑动窗口完成在整个时间序列上的移动后，即可得到空中交通流量波动模式动态转化图。表4描述了在不同单位尺度下每种模式的统计特征，包括每种类型出现次数和转移次数的均值（μ）和标准差（σ）。

在表4 中，VG 与LPVG（N=1～3）允许穿越的节点数量（N取值）越大序模体长度越长，序模体越长导致序模体类型出现次数的均值（μ）以及序模体类型转移频次（μ）变小，说明N取值较大会影响LPVG网络对于时间序列局部结构波动动态的演化结果。例举LPVG（N=1）多元序模体转移矩阵图以及4种情形下序模体转移概率图见图12～13。

图12 序模体矩阵图Fig.12 Sequential motif matrix diagram

表4 序模体动态演化的统计特征Tab.4 Statistical characteristics of sequential motifs dynamic evolutions

通过对多元序模体转移态势进行描绘，VG 以及LPVG（N=1～3）的3，4，5，6节点序模体只有2种序模体类型，并且转移概率几乎相同，所以不具备研究价值。而VG以及LPVG（N=1～3）的5，6，7，8节点序模体所包含的序模体种类为22，32，34，34，由于序模体转移类型过多，而主要的转移模式次数过少，导致序模体类型转移的识别不够明显，转移概率过低。故对4种情形可视图序模体长度优先选择4，5，6，7 节点。4 种情形可视图序模体长度为4，5，6，7节点的序模体转移概率分布图见图14。

通过图14对4种情形序模体转移概率分布情况的描绘，VG网络序模体转移概率方差为0.079 67，且大于LPVG（N=1～3）的3组方差（0.049 8，0.051 89，0.071 47），4 种情形可视图序模体方差最小值为LPVG 的N=1 情况，且能够低于其他3 组网络平均0.0179。

图14 序模体转移概率分布图Fig.14 Transition probability distribution curve of sequential motifs

图13 序模体转移概率图Fig.13 Transition probability graphs of sequential motifs

3.3 4 种情形下表征流量波动模式动态演化的能力评估结果

综上所述，对VG及其N=1～3 取值LPVG网络局部波动动态演变能力分析，结论如下。

1）由于序模体过长在空中交通混沌特性的影响下对于预测没有意义，VG 及其N=1～3 取值LPVG 网络序模体长度推荐使用选择4，5，6，7 节点长度；

2）LPVG（N=1）网络5-node 序模体能够确保序模体类型数量以及序模体类型出现次数和序模体类型转移频次以及处于1 个合理数值范围内，同时序模体转移概率方差低于其他3 组网络平均0.017 9，能维持1个稳定均匀变化的水平。

4 结束语

1）针对天津机场进场航班流量数据，LPVG 网络N 选取数值与节点数量占比在0.48%～1.442%区间能够保证细节损失率在0.5范围内，即使用LPVG网络208（时间粒度为5 min）节点N=1～3。

2）VG 与LPVG（N=1～3）均能有效刻画航班流量时间序列的波动强度，其对时间序列波动动态特性的适应性评估值分别为2.665，4.810，6.973，9.883。

3）由于序模体过小导致类型转移概率趋于相同，过长在交通流混沌特性的影响下对于预测没有意义，故VG及其N=1～3 取值LPVG网络序模体长度推荐使用选择4，5，6，7节点长度。LPVG（N=1）网络5-node序模体能够确保序模体类型数量以及序模体类型出现次数和序模体类型转移频次以及处于1个合理数值范围内，同时序模体转移概率方差平均低于其他3组网络0.017 9。

综上所述，针对以日或月为单位的航班流量时间序列整体波动动态特性的分析，以及网络局部波动动态演化的研究，使用较小时间粒度映射成时间序列，并且选择使用LPVG 网络，N取值为1，深入分析波动模式转移概率矩阵和准确揭示空中交通时间维度的演变规律可以对航班的延误预测提供依据，并对航班实际运行管理提供先导性作用。

本文研究只涉及1 d 及1 个月的进场航班流量数据，所以后续以此进行航班流量时间序列波动特性预测所得到的精度将会受到影响。后续考虑使用1年或者更长时段的航班交通流量数据，并且增加其他相关信息，如天气条件、进离场排序等。