随着空中交通态势在许多因素的综合作用下变得异常复杂，许多问题在还原论思想的指导下得不出合理的解释。在这一时期，复杂性科学蓬勃发展，人们开始利用复杂性理论研究空中交通态势。Prandini等[2]考虑了飞机预测航迹的不确定性，提出了通过概率型的空域使用率评估3-D空域的空中交通复杂度算法，为战略冲突解脱提供机动策略。张进等[3]认为空中交通复杂性贯穿于空中交通管理的始终，已成为认识系统涌现的重要线索。张晨等[4]提出了基于连携效应的交通复杂性测度模型，改进了扇区范围内的改航航段和策略生成方法。
上述研究从宏观角度评估了空中交通的复杂性。宏观指标往往通过对航空器较多的繁忙时段采集获取，是对一段时间的统计信息，缺乏实时性。此外，尽管此类数据较易获取和应用，但是由于运行环境的差异导致数据缺乏普适性。在空域中，飞机之间的相对位置关系能够实时地反映当前空域中的冲突情况和调配难度即空中交通的复杂性。从这个角度出发，一些研究者建立了飞机状态网络，并用复杂网络理论评估空中交通的复杂性。陈新伟等[5]构建了基于复杂网络理论的航线飞行CD&R的基本框架。王红勇等[6-7]中以节点度、连通率和网络结构熵描述网络复杂性，提出以内聚力和稳定性指标反映飞机聚集成簇的程度。Zanin[8]以K-means聚类算法划分空域交通复杂情况。复杂网络分析方法为本文揭示飞行状态网络结构提供了重要的理论工具，为了描述网络结构特性，需要选择反映网络全局特性的复杂网络指标。
以上文献中构建的网络是无权网络，丢失了距离信息，且选择的指标对网络结构信息反映不全面。为了对空中交通复杂性进行综合分析，本文从复杂网络的角度描述了空中交通状况，并采用了一些新的指标进行识别。此外，评估方法也是复杂性识别中的关键问题。经典的评估方法往往需要主观设置阈值，聚类算法同样在很大程度上受类别数k的影响。为了解决上述问题，本文提出了一种空中交通复杂性识别的黑盒方法，通过对样本的学习直接找出复杂性识别的控制限，避免了阈值设置上的主观性。首先对空中交通状况进行飞行状态网络建模，选取了5个静态拓扑指标映射空中交通复杂性，然后引入独立主元分析(Independent Component Analysis, ICA)方法对空中交通复杂态势进行在线监控[9-11]。通过本文所提方法，可以实时获取空情信息，分析空中交通态势及其演化过程，为战略层飞行冲突预警提供决策支持。通过人工网络和实际飞行状态网络验证了所提方法的有效性。
1 飞行状态网络建模方法 1.1 网络的构建 以飞机为节点，当潜在冲突发生时，飞机通过机载防相撞系统(Airborne Collision Avoidance System，ACAS)与邻机通信，并立即采取行动解决冲突。如果ACAS通信距离的范围覆盖了邻近的飞机，则节点之间有一条边。这里，ACAS通信距离RACAS=40 km。飞行状态网络结构差异性如图 1所示。
 |
| 图 1 飞行状态网络结构差异性 Fig. 1 Difference of flight state network structures |
| 图选项 |
图 1中，相同的空域内分布着5架同型航空器，其ACAS通信距离均为40 km，3个空域按照复杂程度从小到大排序依次为：空域A、空域B、空域C。图 1(a)中，各飞机在空域中分布相对分散，反映在复杂网络中即为各节点相互孤立；图 1(b)中，部分飞机相互靠近聚集成簇，在彼此ACAS通信距离内，开始构成连边；图 1(c)中，所有飞机聚集在空域的中心区域，空中态势拥堵程度高，网络全连通甚至出现环边，极不利于管制员调配冲突。
此外，航空器间的距离直接影响飞行安全，本文设想将距离因素也考虑进飞行状态网络内。当两机间隔越小，冲突越突出，网络也越紧密。因此，本文将两机距离的倒数设置为网络边权，将上述模型拓展为加权网络模型，其权值设置如图 2所示。
 |
| 图 2 边权设置 Fig. 2 Edge weight setting |
| 图选项 |
当两机距离d=35 km时，两机ACi和ACj连边的边权为wij=1/35。由此看出，当空域中航空器数量相当时，空中交通复杂性受网络结构影响很大。
1.2 网络拓扑指标的选取 在一个无向飞行状态网络G={V, Ed, We}中，飞机节点表示为V={vi|i∈I, I=1, 2, …, N}，N为节点数，节点vi与vj之间存在连边Ed={edij=(vi, vj)|i, j∈I}V×V，Z={zij}n×n为网络的邻接矩阵，其中，元素zij可以表示为
 | (1) |
We={weij}n×n为网络的加权邻接矩阵，当节点之间不存在连边时，权值为无穷。选取以下指标，映射空中交通复杂性：
1) 环边数(LN)
环边数是一个能反映网络复杂程度的简单指标，环边数越多说明网络越紧密。本文采用深度优先搜索(Depth-first Search，DFS)算法计算网络环边数[12]。
2) 平均点强(ANS)
加权的节点度即是节点点强si，平均点强si的表达式为
 | (2) |
节点点强越大，说明飞机节点与周围飞机连接紧密。
3) 平均加权聚类系数(ACC)
平均加权聚类系数描绘了网络的聚集性特征[13]，单个节点的加权聚类系数表达式为
 | (3) |
式中：ki为节点度值。vj、vk为节点vi的2个相邻节点，zij反映节点对的连接状态，当vi和vj互连，zij=1，否则zij=0。为使c(i)∈[0, 1]，式(3)右侧除以si(ki-1)。平均加权聚类系数
4) 平均点介数(ABC)
介数中心性是级联故障现象与网络动力学行为的概念，反映了节点或连边在整体网络中的中心程度，点介数表达式为
 | (4) |
式中：σij(kv)为节点vi与节点vj间经由节点vkv的最短路径数目；σij为节点vi与节点vj间最短路径的数目。介数大的节点在网络中占据中心位置，对其他节点的行为产生更大影响[14-15]。
5) 网络效率(NE)
网络效率是全部节点间的距离倒数和的平均值，其表达式为
 | (5) |
式中：dij为节点vi与vj间的距离。网络效率可以反映网络信息传输的难易程度，Ew越大，信息传递越顺畅，网络鲁棒性就越强。
2 独立主元分析监测方法 2.1 独立主元分析方法 ICA方法是近年来发展起来的一种统计信号处理方法。对于给定的d个观测变量x(1), x(2), …, x(d)，它们可以表示为m(m≤d)维未知的独立成分s(1), s(2), …, s(m)的线性组合。其中，独立成分和观测变量都已进行零均值标准化处理，它们之间的关系可以表示为
 | (6) |
式中：X=[x(1), x(2), …, x(ns)]∈Rd×ns为观测变量；S=[s(1), s(2), …, s(ns)]∈Rm×ns为独立成分变量；A=[a(1), a(2), …, a(m)]∈Rd×m为未知的混合矩阵；E∈Rd×ns为残差矩阵，ns为样本数量。若d=m，残差矩阵E为一个零矩阵。
ICA方法要解决的主要问题，即是仅通过观测变量X，估计出混合矩阵A和独立成分变量S。本文定义一个解混矩阵W，实现对独立成分变量的估计：
 | (7) |
式中：?为独立成分变量S的估计量，当解混矩阵W是混合矩阵A的逆矩阵(W=A-1)时，?为S的最优估计。
ICA方法首先需进行白化处理，以消除随机变量之间的相关性。对于样本ky有d维随机变量x(ky)，其协方差为Rx=E(x(ky)xT(ky))。Rx可特征分解为
 | (8) |
白化处理可表示为
 | (9) |
式中：X(ky)为样本ky的观测变量; 白化矩阵Q=Λ-1/2UT。在这一白化处理下，Rz=E(z(ky)·zT(ky))为单位阵。因此在该处理下，式(9)可转化为
 | (10) |
式中：矩阵B为正交矩阵。因此，该问题从求解一个任意的满秩矩阵A简化为求解一个正交阵B。由式(10)可对s(ky)进行估计：
 | (11) |
由式(7)和式(11)，可表示出解混矩阵W和正交阵B的关系：
 | (12) |
本文采用Hyv?rinen提出的Fast-ICA算法[16-17]计算解混矩阵W，其主要步骤如下：
步骤1??随机选择初始化权值向量w。
步骤2??令wnew=E{xg(wTx)}-E{g′(wTx)}w。
步骤3??令w=wnew/||wnew||。
步骤4??wk=wk-
步骤5??如果结果不收敛，则转到步骤2。
步骤6??令ky=ky+1，如果ky小于独立成分(ICs)的数量，则转到步骤1。
在步骤2中，函数g通常取：
 | (13) |
 | (14) |
 | (15) |
式中：1≤a1≤2；a2≈1。
2.2 基于ICA的异常监测统计量 利用正常运行下的离线历史数据估计出统计量的“阈值”，称之为控制限，当在线监测样本的统计值超过该控制限时，则认为异常情况发生。
在正常运行下，有观测变量x(1), x(2), …, x(d)，在独立成分空间中标准化后的数据集为Xnormal。在观测变量维度与独立成分变量维度相等的情况下，通过Fast-ICA算法可解得解混矩阵W和独立成分的估计量?(?=WX)。同样地，在Fast-ICA运算过程中，矩阵A、B和Q也可求解到。将解混矩阵W各向量wi范数按大小排列，选取矩阵W前a列构成约化矩阵Wd(W的主要部分)，剩余列构成约化矩阵We(W的次要部分)。因此，新的主元空间和残差空间中的独立成分向量?new d(ky)和?new e(ky)，可由式(16)、式(17)求得
 | (16) |
 | (17) |
在主成分分析(Principal Component Analysis, PCA)方法中，有2类统计量用于过程的监控。D统计量针对过程变量的主元部分，而Q统计量针对过程变量的残差部分。类似地，本文将这2个统计量也用于ICA算法的过程监控中。D统计量，即I2统计量，对于样本ky而言，可由式(18)定义：
 | (18) |
Q统计量，也即SPE统计量，其定义如下：
 | (19) |
式中：
 | (20) |
此外，Simoglou等[18]基于典型关联分析(Canonical Variate Analysis, CVA)提出了监控系统状态的T2统计量。本文借鉴这一方法，采用基于残差空间的独立主元?new e(ky)的统计量Ie2监控异常状态。该统计量的优势在于能够补偿独立主元数量选择不恰当时的误差，Ie2统计量定义如下：
 | (21) |
Ie2与I2统计量互为补充，共同监控系统的运行状态。
2.3 控制限的确定 在离线获取正常数据后，样本m(m=1, 2, …, N)对应的统计值I2(m)、Ie2(m)和SPE(m)可由2.2节式(18)~式(21)计算求出。同样，在线监测过程中，采样时刻kt的样本统计量I2(kt)、Ie2(kt)和SPE(kt)也可实时计算得出。通过核密度估计(Kernel Density Estimation)方法计算正常样本统计特征的控制限(Control Limit)，绘制ICA监测图，比较在线监测样本的统计值与控制限的大小，若超过控制限，则认为该采样时刻的样本存在异常。
在PCA监测中，Hotelling的T2统计量分析和SPE统计图是提取数据关键特征的重要工具，而这一方法是在假设潜变量的概率密度函数服从多变量高斯分布的基础上进行的。然而在ICA监测方法中，本文无法预先知道独立成分的统计分布情况，因此，统计量的控制限不能直接通过某一具体分布求取。本文采用核密度估计方法计算I2、Ie2和SPE的控制限。对于核K的单变量核密度估计量
 | (22) |
式中：q为估计数据点; qi为观测值; h为平滑参数；no为观测值的数量；K为核函数。因此，核估计量是观测出“波动”的总和。核函数K定义了这种“波动”的模式，以及需要满足以下条件：
 | (23) |
核函数的选择对模型影响有限，而高斯核函数对于大多数模型都是适用的，选其作为核函数。求取ICA监测中的控制限主要有如下3个步骤：
步骤1??通过离线正常数据计算I2、Ie2和SPE统计量。
步骤2??使用单变量核密度估计量计算I2、Ie2和SPE统计量的密度函数。
步骤3??选择99%置信水平的点作为该统计量的控制限。
与通过Hotelling的T2统计量获取控制限相比，使用核密度估计方法求取控制限的主要优势在于与数据联系更紧密。
2.4 贡献图 通过I2、Ie2和SPE三个统计量的ICA监测图可以在线识别出数据的异常情况，然而，监测图只能反映异常是否存在，不能诊断出是何种故障或是哪些变量的问题导致了异常的发生。以空中交通复杂性识别为例，监测图仅能体现该状态下空域是否超出控制限，进入拥挤状态，而不能反映出拥挤产生的原因。因此，贡献图(Contribution Plots)方法是ICA异常诊断中分析故障原因的重要方法。当采样时刻kt发生异常，画出该时刻各统计量的监测变量贡献图，能够反映各监测变量对于模型或残差的影响大小。
在ICA监测中，变量x(kt)对于I2(kt)和Ie2(kt)的贡献可分别由式(24)和(25)表示：
 | (24) |
 | (25) |
式中：变量贡献值的平方和分别等于该采样时刻I2、Ie2的值，即xcd(kt)xcdT(kt)=I2(kt)，xce(kt)xceT(kt)=Ie2(kt)。同样，SPE统计量的变量贡献值也可以计算出来。采样时刻kt，SPE统计量的定义由式(19)可知，有SPE(kt)=eT(kt)e(kt)。因此，向量e(kt)包含了每个监测变量特有的预测误差，绘制e(kt)各元素的柱状图作为SPE统计量的贡献图。
3 方法流程 基于此，本文提出了一种基于复杂网络理论的飞行状态网络建模方法，在此基础上，通过ICA异常诊断方法在线评估空中交通复杂性，将超出控制限的样本称为拥挤样本，低于控制限的，认为其较为流畅，以此监测空中交通拥挤。本文方法主要流程如图 3所示。
 |
| 图 3 本文方法主要流程 Fig. 3 Main flow of proposed method |
| 图选项 |
由图 3可知，提出的空中交通复杂性评估方法主要分为2个步骤：复杂网络建模与ICA在线监控。其中，ICA在线监控空中交通复杂性状况的步骤如下。
1) 离线训练阶段
步骤1??收集空中交通运行流畅状态(正常状态)下的数据作为训练集。
步骤2??对训练集进行预处理(中心化)，之后进行白化处理，以消除变量之间的相关性。
步骤3??利用Fast-ICA算法计算混合矩阵A和解混矩阵W。
步骤4??计算SPE、I2和Ie2统计量。
步骤5??利用核密度估计方法计算各统计量的控制限。
2) 在线监控阶段
步骤1??以待监测的数据为测试集，进行预处理。
步骤2??在线计算监测数据的SPE、I2和Ie2统计量，评估空中交通复杂性。
步骤3??若某统计量大于对应控制限，则认为该采样时刻空中交通拥挤；否则认为空中交通流畅。
4 仿真分析 为了验证本文方法的有效性，首先对空中交通态势进行仿真，考察在仿真环境中的监测效果；其次以ZPPP-昆明长水机场进近阶段数据为监测样本，对某日15:40—16:05期间的空中交通情况进行动态评估。
4.1 数据采集 在仿真场景中，空域大小设置为200 km×200 km。为了方便描述空域飞机分布情况，随机产生一批飞机但其位置服从特定分布，且从样本多样性角度考虑，本文假设飞机位置服从均匀分布、正态分布和指数分布，且x、y坐标值之间相互独立。首先，本文认为飞机(节点)数对空中交通复杂性影响最大，这在飞行状态网络结构中能够体现。以飞机位置服从均匀分布为例，节点数N=40, 80，其网络结构变化可如图 4所示。
 |
| 图 4 网络结构随节点数增多的变化情况 Fig. 4 Variation of network structure with increase of node number |
| 图选项 |
固定空域中，当节点数很少时，如图 4(a)所示，节点之间相互连接，但网络仍不是全连通的；当节点数量增加到80时，如图 4(b)所示，网络全连通，且随着飞机数量的增多，网络规模进一步扩大。经过大量测试发现，该空域中节点数N=40时，空中交通仍比较流畅，结论的依据是大部分样本中的飞行冲突数量都小于5。这里，二维平面中飞行冲突定义为两机间隔小于5 min。因此，本文采集节点数N=20, 30, 40且服从均匀分布的样本各50组，作为正常数据集。飞机位置服从均匀分布是一种空域利用率最高的理想状态，一般情况下，民航飞机按照机场-进近-区域飞行，军航飞机根据作训任务的不同，灵活选择空域或航线，鲜有文献对位于空域中的飞机位置分布规律进行统计分析。这里，为了方便描述，本文假设飞机位置服从正态分布和指数分布，且节点数均为N=80，网络结构如图 5所示。
 |
| 图 5 相同节点数不同分布时网络结构的差异 Fig. 5 Difference of network structure with the same node number but different distribution of nodes |
| 图选项 |
与图 4(b)相比，同样是80个节点，节点之间紧密程度由强至弱依次为：指数分布，正态分布，均匀分布。特别是飞机坐标服从指数分布时，大量飞机聚集在中心位置，飞行冲突突出，空中交通复杂程度过高，拥挤严重。本文以这2种情况的样本为在线监测测试集，逐渐减少节点数，考察提出的方法对空中交通拥堵识别的效果。
在实际飞行场景中，本文以ZPPP-昆明长水机场进近阶段数据为样本，选择交通流畅的样本训练后，对某日15:40—16:05期间的空中交通情况进行了动态评估。每隔5 min，以快照(Snapshot)的形式记录下实时的空中交通情况，这样动态复杂网络被离散为6个静态子网络(Sub-networks)，通过提出的方法对空中交通状况进行在线监控。
4.2 仿真场景 首先，本文从静态子网络的角度出发，考察该方法的有效性。即在空域中随机产生一批飞机，记录下此时的飞机位置信息，收集若干样本作为测试集，在离线训练好的模型中在线监测。采集飞机数量N=20, 30, 40且服从均匀分布的样本各50组，作为正常数据集，且这些数据中飞行冲突数量均小于5。首先，本文探究所提方法在飞机数量增加时对空中交通拥堵的识别效果，测试集为服从均匀分布的50组节点数量均为N=50的样本，其中，飞行冲突数量均大于8。前150个样本为正常数据集，后50个样本为监测数据集，其结果如图 6所示。
 |
| 图 6 监测图(均匀分布且节点数为50) Fig. 6 Monitoring charts (uniformly distributed and node number equals to 50) |
| 图选项 |
从图 6中可以看出，提出的方法能够准确识别空中交通异常。I2和Ie2统计量从第151个样本开始监测到异常，不存在延迟和错误诊断。SPE统计量在训练过程中不稳定，导致控制限偏高；在监测阶段，从第151个样本起虽然SPE统计值明显高于训练阶段统计值，但大部分仍略低于控制限，未能判断为拥挤样本。
为了进一步说明是何种原因造成了交通拥挤，本文给出第175个样本的贡献图，变量1~5依次为环边数、平均点强、平均加权聚类系数、平均点介数和网络效率指标，结果如图 7所示。
 |
| 图 7 变量对偏差的贡献 Fig. 7 Contribution of variables to deviation |
| 图选项 |
从图 7中可以看出，网络效率对空中交通拥堵影响最大，其次为环边数、平均加权聚类系数和平均点强，平均点介数最弱。这是因为随着空域中飞机数量的不断增加，网络结构趋于复杂，这些从网络整体复杂度出发的指标最为敏感，而由于网络结构类似于随机网络，平均点介数和平均加权聚类系数指标相对变化较小。
另外，本文比较关心该方法在不同网络结构中的效果，这里以飞机坐标服从正态分布和指数分布的样本为测试集，训练集同上一试验。由于这2种情况下交通状态明显拥堵，本文设置空域中的飞机数量与训练集中相当，N=40, 30, 20，考察方法对复杂性的评估能力，监测结果如图 8所示。
 |
| 图 8 不同节点数和分布时的监测图 Fig. 8 Monitoring charts with different node number and distribution |
| 图选项 |
图 8(a)~图 8(c)服从正态分布，图 8(d)~图 8(f)服从指数分布，且各组样本中节点数量从40依次减少到20，监测样本中飞行冲突数量均大于8。当飞机位置服从正态分布且节点数量N=40和30时，如图 8(a)、图 8(b)所示，拥堵情况被部分统计量准确识别，验证了该方法对于不同结构的有效性。当节点数量减少到20时，如图 8(c)所示，I2统计量的控制限因训练阶段的一个错误样本而严重偏高，而SPE和Ie2统计量的监测样本与训练样本存在比较明显的差异，在控制限上下波动，空中交通已经比较流畅。当飞机位置服从指数分布时，如图 8(d)、图 8(e)和图 8(f)所示，拥挤情况均被所有或部分统计量识别出来。特别图 8(f)与图 8(c)相比，虽然节点数均为20，但指数分布的异常情况被SPE统计量识别出来，说明了相同飞机数量的情况下，指数分布比正态分布更加拥挤。上述试验验证了该方法对于静态网络的复杂性识别有效性，样本中各网络相互独立，没有时间上的先后顺序。
4.3 空管运行场景 在实际空管运行场景中，本文以ZPPP-昆明长水机场进近阶段数据为样本，对某日15:40—16:05期间的空中交通情况进行动态评估。每隔5 min，以快照的形式记录下实时的空中交通情况，这样动态复杂网络被离散为6个静态子网络，如图 9所示。
 |
| 图 9 进近阶段不同时刻雷达屏幕截图 Fig. 9 Radar screenshots at different moments in phase of approaching |
| 图选项 |
进近阶段空中交通态势如图 9所示，一个动态复杂网络按照相同的时间间隔5 min离散为6个静态子网络，客观反映了空中交通态势的动态变化过程。通过复杂网络建模，雷达管制实况可抽象为如图 10所示的网络结构图。
 |
| 图 10 进近阶段不同时刻离散子网络 Fig. 10 Discrete sub-network structures at different moments in phase of approaching |
| 图选项 |
从图 10可以比较直观地看出，随着将飞机从进近阶段向机场管制区域移交，空中交通态势的总体趋势从拥挤逐渐转变为流畅。然而，本文需要一种研究方法对空中交通复杂性进行定性的分析。因此，本文在该场景下验证提出方法的有效性。为了获得理想的监测效果，本文将管制区容量考虑在内，在大量的真实数据中选择了50个交通流畅的正常样本作为训练集。在表 1中列出了部分样本的5个拓扑指标。通过核密度估计方法计算出的SPE、I2和Ie2统计量的控制限分别为3.335 7×10-29、15.345 2和9.752 1。在异常监测阶段，6组样本的拓扑指标和对应的SPE、I2和Ie2统计值如表 2所示。
表 1 部分训练样本拓扑指标值 Table 1 Some topological indicator values of training samples
| 样本序号 | LN | NS | CC | BC | NE |
| 1 | 53 | 9.754 6 | 0.813 4 | 0.021 9 | 40.125 0 |
| 2 | 41 | 9.145 8 | 0.827 4 | 0.029 3 | 20.569 4 |
| 3 | 60 | 10.883 2 | 0.842 7 | 0.012 9 | 35.411 3 |
| 4 | 42 | 8.433 9 | 0.775 3 | 0.034 5 | 20.882 2 |
| 5 | 59 | 11.977 7 | 0.852 3 | 0.015 4 | 24.790 2 |
| 6 | 52 | 9.567 7 | 0.817 6 | 0.021 1 | 27.105 3 |
| 7 | 50 | 9.617 3 | 0.803 6 | 0.022 8 | 31.494 1 |
| 8 | 20 | 11.570 1 | 0.723 4 | 0.029 2 | 21.478 3 |
| 9 | 41 | 7.796 3 | 0.813 3 | 0.031 0 | 19.480 9 |
| ┇ | ┇ | ┇ | ┇ | ┇ | ┇ |
| 50 | 60 | 14.536 2 | 0.855 2 | 0.013 9 | 24.163 0 |
表选项
表 2 监测样本拓扑指标值和SPE、I2和Ie2统计值 Table 2 Topological indicator values of monitoring samples and statistic values of SPE, I2 and Ie2
| 时刻序号 | 时刻 | LN | NS | CC | BC | NE | SPE | I2 | Ie2 |
| 1 | 15:40:43 | 207 | 23.681 4 | 0.934 6 | 0.001 3 | 185 983 | 7.467 1×10-29 | 21.867 4 | 3.464 4 |
| 2 | 15:45:43 | 122 | 24.009 2 | 0.964 0 | 0.001 0 | 145.370 9 | 5.218 6×10-29 | 16.758 3 | 4.589 4 |
| 3 | 15:50:43 | 65 | 18.026 1 | 0.951 9 | 0.003 5 | 43.537 8 | 3.921 4×10-29 | 14.541 3 | 1.976 5 |
| 4 | 15:55:44 | 44 | 18 223 | 0.894 1 | 0.013 5 | 65.336 2 | 3.676 3×10-29 | 13.287 6 | 5.875 4 |
| 5 | 16:00:43 | 33 | 8.149 0 | 0.884 7 | 0.010 1 | 46.486 1 | 9.812 5×10-30 | 6.221 8 | 2.664 7 |
| 6 | 16:05:43 | 39 | 10.791 1 | 0.892 8 | 0.019 4 | 41.822 5 | 1.581 9×10-29 | 7.543 9 | 2.545 2 |
表选项
将监测样本统计值与控制限比较，本文发现SPE和I2统计量识别到了部分异常情况，而Ie2统计量总是低于控制限。SPE统计量判断前4个样本为拥挤态势，而I2统计量识别到前2个样本为拥挤态势，说明SPE统计量相较于I2统计量对指标的变化更加敏感。在这6个采样时刻中，根据SPE和I2统计量的评估结果，态势从拥挤到流畅的复杂性排序均为：1>2>3>4>6>5，差异仅在于2个统计量对于指标变化的敏感程度。因此，本文认为在进近区域的空中交通复杂性识别中，在判断空中态势是否拥挤时，SPE和I2统计量都应给予充分关注。为了验证提出方法的评估效果，利用文献[19]中的聚类算法测试同样的管制运行数据，与本文方法进行比较。文献[19]将不同维度和阶段的拓扑指标组合成一个综合复杂性指标向量，使用K-means算法将样本依据空中交通复杂性划分为3类。各指标和评估结果如表 3所示，表中：Ei, j为第i层第j级的边数，Ci, j为该层级的聚类系数。
表 3 文献[19]中K-mean算法对相同样本的复杂性识别结果 Table 3 Complexity recognition results of K-mean algorithm for the same sample in Ref.[19]
| 时刻序号 | 时刻 | N | E1, 1 | C1, 1 | E1, 2 | C1, 2 | E2, 1 | C2, 1 | E2, 2 | C2, 2 | 等级 |
| 1 | 15:40:43 | 24 | 87 | 0.934 6 | 71 | 0.923 8 | 66 | 0.954 7 | 72 | 0.955 8 | 高 |
| 2 | 15:45:43 | 18 | 72 | 0.964 0 | 81 | 0.962 7 | 45 | 0.943 1 | 64 | 0.942 5 | 高 |
| 3 | 15:50:43 | 14 | 66 | 0.931 9 | 51 | 0.925 4 | 41 | 0.922 1 | 35 | 0.927 6 | 低 |
| 4 | 15:55:44 | 13 | 58 | 0.914 1 | 52 | 0.912 4 | 45 | 0.935 0 | 33 | 0.898 8 | 中 |
| 5 | 16:00:43 | 11 | 35 | 0.884 7 | 34 | 0.882 1 | 18 | 0.884 6 | 16 | 0.885 6 | 低 |
| 6 | 16:05:43 | 12 | 41 | 0.892 8 | 40 | 0.881 7 | 23 | 0.883 2 | 25 | 0.853 0 | 低 |
表选项
从复杂性评估等级可以看出，前2个样本被识别为“高”等级，第4个样本为“中”，其余样本为“低”，复杂性的总体变化趋势与所提方法的结果几乎相同，但个别样本的结果存在差异。从方法的角度看，聚类算法缺乏训练过程，因此获得的复杂性等级是相对测试样本自身而言的。该模型受类别数k的影响；而在提出方法中，控制限是通过大量顺畅样本训练获得的，复杂性评级对指导飞行活动更具有参考意义。从仿真结果看，聚类算法认为第3个样本复杂性等级为“低”，第4个样本是“中”，但所提方法识别出前4个样本为复杂样本，且等级从高到低依次为：1>2>3>4>6>5。从图 10中可以看出第3个样本的复杂性明显比第4个更高。在类别数k=2的条件下再次测试，此时，前4个样本被评价为“高”，后2个样本被评价为“低”，结果与所提方法相同。此试验说明了聚类算法受类别数k的影响较大，进一步验证了此方法的有效性与鲁棒性。
5 结论 1) 针对战略层面冲突探测与解脱特点，本文提出了一种基于复杂网络分析和ICA异常监测算法的空中交通复杂性识别方法，弥补了中期冲突探测系统对空中交通态势整体评估的不足。
2) 通过构建飞行状态网络，利用拓扑指标与空中交通复杂性建立映射关系，选择的指标能够较全面地刻画网络结构和性能，为了客观地评估空中交通状况，引入了ICA监测的黑盒方法，通过训练交通顺畅样本计算控制限，将测试样本的SPE、Ie2和I2统计量与控制限比较来识别复杂情况，克服了其他评估方法阈值设置的主观性问题。
3) 在仿真场景和管制运行的真实场景中，3个统计指标从不同角度反映复杂性，互为补充，较好地监测空中交通态势，方法稳定性、普适性较强。
本文方法仅从空域的情况分析了空域中的交通复杂情况，空域情况是相对管制员而言的。因此，后期结合管制员的管制难度和工作负荷对空域复杂性进行分析，将能够更有针对性的体现空中交通情况。
参考文献
| [1] | ADAM R, JACEK S. The concept of initial air traffic situation assessment as a stage of medium-term conflict detection[J]. Procedia Engineering, 2017, 187: 420-424. DOI:10.1016/j.proeng.2017.04.395 |
| [2] | PRANDINI M, PUTTA V, HU J H. A probabilistic measure of air traffic complexity in 3-D airspace[J]. International Journal of Adaptive Control and Signal Processing, 2010, 24(10): 813-829. DOI:10.1002/acs.1192 |
| [3] | 张进, 胡明华, 张晨. 空中交通管理中的复杂性研究[J]. 航空学报, 2009, 30(11): 2132-2142. ZHANG J, HU M H, ZHANG C. Complexity research in air traffic management[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(11): 2132-2142. (in Chinese) |
| [4] | 张晨, 胡明华, 张进, 等. 基于交通复杂性的扇区资源管理[J]. 南京航空航天大学学报, 2010, 42(5): 607-613. ZHANG C, HU M H, ZHANG J, et al. Sector asset management based on air traffic complexity[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2010, 42(5): 607-613. (in Chinese) |
| [5] | CHEN X W, STEVEN J L, SHIMON Y N. A framework of enroute air traffic conflict detection and resolution through complex network analysis[J]. Computer in Industry, 2011, 62: 787-794. DOI:10.1016/j.compind.2011.05.006 |
| [6] | WANG H Y, WEN R Y, ZHAO Y F. Analysis of topological characteristic in air traffic situation networks[J]. Proceedings of the Institution of Mechanical Engineers, Part G:Journal of Aerospace Engineering, 2015, 229(13): 2497-2505. DOI:10.1177/0954410015578482 |
| [7] | WANG H Y, XU X H, ZHAO Y F. Empirical analysis of aircraft clusters in air traffic situation networks[J]. Proceedings of the Institution of Mechanical Engineers, Part G:Journal of Aerospace Engineering, 2017, 231(9): 1718-1731. DOI:10.1177/0954410016660870 |
| [8] | ZANIN M. Network analysis reveals patterns behind air safety events[J]. Physica A:Statistical Mechanics and Its Applications, 2014, 401: 201-206. DOI:10.1016/j.physa.2014.01.032 |
| [9] | HYVARNEN A, OJA E. A fast fixed-point algorithm for independent component analysis[J]. Neural Computation, 1997, 9: 1483-1492. DOI:10.1162/neco.1997.9.7.1483 |
| [10] | HYVARNEN A, OJA E. Independent component analysis:Algorithms and applications[J]. Neural Networks, 2000, 13(4-5): 411-430. DOI:10.1016/S0893-6080(00)00026-5 |
| [11] | LEE J M, YOO C K, LEE I B. Statistical process monitoring with independent component analysis[J]. Journal of Process Control, 2004, 14(5): 467-485. DOI:10.1016/j.jprocont.2003.09.004 |
| [12] | TANG J J, WANG Y H, LIU F. Characterizing traffic time series based on complex network theory[J]. Physica A:Statistical Mechanics and Its Applications, 2013, 392: 4192-4201. DOI:10.1016/j.physa.2013.05.012 |
| [13] | TARJAN R.Depth-first search and linear graph algorithms[C]//Symposium on Switching & Automata Theory.Piscataway: IEEE Press, 1971: 114-121. |
| [14] | FREEMAN L C. A set of measures of centrality based on betweenness[J]. Sociometry, 1997, 40(1): 35-41. |
| [15] | WANG X Y, LI J Q. Detecting communities by the core-vertex and intimate degree in complex networks[J]. Physica A:Statistical Mechanics and Its Applications, 2013, 392: 2555-2563. DOI:10.1016/j.physa.2013.01.039 |
| [16] | WANG X Y, CAO J Y, LI R, et al. A preferential attachment strategy for connectivity link addition strategy in improving the robustness of interdependent networks[J]. Physica A:Statistical Mechanics and Its Applications, 2017, 483: 412-422. DOI:10.1016/j.physa.2017.04.128 |
| [17] | PATERA R P. Space vehicle conflict probability for ellipsoidal conflict volumes[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(6): 1818-1821. |
| [18] | SIMOGLOU A, MARTIN E B, MORRIS A J. Statistical performance monitoring of dynamic multivariate processes using state space modelling[J]. Computers & Chemical Engineering, 2002, 26(6): 909-920. |
| [19] | WANG H, SONG Z, WEN R, et al. Study on evolution characteristics of air traffic situation complexity based on complex network theory[J]. Aerospace Science and Technology, 2016, 58: 518-528. DOI:10.1016/j.ast.2016.09.016 |
