基于IPv6 AS级复杂网络的特征分析与建模

刘 冰 杨 学 杨 琪 马永征

(中国互联网络信息中心 北京 100190)

0 引 言

随着互联网技术的迅猛发展和网络规模的日益扩大，IPv4协议在IP地址数量、安全性、移动性、服务质量等方面的不足日渐凸显，IPv6应运而生。伴随着IPv6网络规模的持续增长及各项相关技术的日趋完善，其重要性逐渐提高，网络结构也越发复杂。而全球互联网由自治系统(Autonomous System，AS)互相连接而构成，AS的工作方式、复杂的相互关系和自身属性都在很大程度上决定了整个网络的流量和行为[1]。互联网AS级拓扑关系能够反映出网络空间内各AS间的关联情况，据此可进一步分析得到其相互间的输入输出策略。因此，AS间的连接关系是网络拓扑结构中至关重要的一环，针对IPv6 AS级网络拓扑的分析和研究成为网络管理和网络安全控制的重要方向[2]。通过对IPv6网络拓扑的全面分析和建模，可以获得不同组织机构之间的关联关系和路由策略情况，为整体网络空间的优化提供了依据；能够加深对IPv6网络发展的全面理解，有利于对未来IPv6网络发展趋势作出更好的判断；能够进一步深入分析IPv6网络的内在机制，对更加有效地部署和规划下一代互联网有重要的指导作用，对国家的网络管理和维护具有深远的意义。

目前，IPv6 AS级网络结构已经非常复杂(如图1所示)，然而国内外对网络拓扑结构及模型的研究依然主要集中在IPv4层面，对IPv6网络拓扑的分析成果较少，仅有为数不多的基于IPv6地址构成的网络拓扑特征量的研究[3-6]。本文以时间为主线，通过对IPv6 AS级网络的静态特征和网络特性随时间的演化进行分析，并利用网络建模技术对未来IPv6 AS网络发展趋势作出分析和预测。

图1 2018年IPv6 AS复杂网络示例图

1 数据采集与预处理

1.1 数据采集

AS之间一般是借助专线或者公共网络的接入点实现连接的，而路由宣告和传送是利用AS间的路由协议——边界网关协议(Border Gateway Protocol，BGP)[7]确定的，BGP协议使AS可以根据其实际情况选择不同的路由策略来完成路由选择。本文基于RouteViews项目[8]获取BGP路由数据，提取出相关的AS信息做进一步分析。

RouteViews项目由俄勒冈大学(University of Oregon)先进网络技术中心创立，用户可以从互联网上的几个不同骨干和位置获取有关全球路由系统的实时BGP信息。RouteViews项目的初衷是能够向互联网服务提供商(ISP)提供查询其他网络前缀的服务帮助，使其能够很方便地调试他们接受到的网络访问，以便优化其自身网络情况。近几年内，该项目已经在学术和科研项目[9]等方面得到了广泛应用。1997年11月起，RouteViews项目就已经开始收集相关的BGP数据，截至2019年4月已有24个数据采集点，分别并行地采集全球互联网的BGP数据。我们通过具体数据比对发现，24个数据采集点中除位于肯尼亚、佩斯和贝尔格莱德的3个采集点外，其他采集点各自获取的AS总数基本一致，边总数(即AS节点连接数)相差不大，对网络结构状的分析影响有限。因此，为简单起见，本文在分析IPv6的AS级复杂网络特征与拓扑结构时，仅选取RouteViews项目位于俄勒冈的IPv6专用采集点[10](route-views6)2003年5月至2018年5月共15年的IPv6 BGP路由表原始数据。

1.2 数据预处理

目前，RouteViews项目的BGP数据采集工作是以2小时为周期进行，每天数据共有12个压缩包，我们通过解析与合并，对该12个压缩包提取出的AS和对应边取合集，作为一天的BGP数据。由于从2003年至2018年长达15年的BGP数据反映了互联网在IPv6维度方面的网络结构动态演化过程，对每月或每年数据进行数据合并会屏蔽掉时间序列的特性(以月数据为例，将整个月的数据进行合并时，会屏蔽掉节点的撤销与新增等时间过程信息。例如：1月2日新增了节点AS1，25日撤销了该节点，当前网络结构中宣告节点应该不包含此AS，但将1月1日-1月31日的数据做合并后，AS列表中将会出现此节点，与实际网络状态不符)。因此，我们取每月最后一天的数据代表该月的情况，取每年最后一天的数据代表该年的情况，以此保证时序信息，使数据更加符合实际情况。

BGP原始数据包解压后为MRT格式，需要通过BGP数据包解析工具[11]将此二进制文件转换为可读格式。解析原始数据后提取AS_PATH字段的ASN相关数据，AS_PATH属性包含了一个有序的ASN列表，描述了到达目标网络所要经过的ASN序列。AS_PATH属性只有在BGP发言者向外部边界网关协议(EBGP)的邻居发布路由时，才会将自己本地系统的ASN作为最后一个元素添加到序列的最左边，而在向内部边界网关协议(IBGP)邻居发布路由时，并不会修改AS_PATH属性。因此，在提取AS网络结构的边数据时，需要根据AS_PATH属性中的ASN序列特性(右边为起始ASN，左边为目标ASN)成对提取。

AS_PATH具有4种类型：

1) AS_SEQUENCE(用于路由AS路径记录)。

2) AS_SET(用于聚合路由的明细路由AS集合)。

3) AS_CONFED_SEQUENCE(用于联盟路由AS路径记录)。

4) AS_CONFED_SET(用于联盟聚合路由)。

在BGP路由表中的显示格式[12]，如图2所示。

图2 AS_PATH显示格式

在对AS_PATH进行提取时，需进行如下处理：

1) 压缩处理。例如AS_PATH为“47065 1200 1200 5555”，可压缩为“47065 1200 5555”。由于BGP协议允许路由器将本地AS连续多次添加到AS_PATH属性中，使整体路由长度变大，干扰对等体路由的决策，因此在数据提取时需对此类路径做压缩处理。

2) 单节点路径。例如“47065”或“{47065}”，仅提取AS号加入节点数据列表中。

3) 多节点路径。例如“47065 1200 5555”，按照从右到左的顺序提取边的起始和终止节点。

4) 包含AS_SET多节点路径。例如“47065,{1200;3209;3320},5555”，拆分“{}”内的AS号与括号外的AS号分别组成新序列为“47065 1200 5555”、“47065 3209 5555”、“47065 3320 5555”，再重复步骤3)。

5) 过滤包含AS_CONFED_SEQUENCE或AS_CONFED_SET的路径。

6) 过滤存在回路的AS_PATH。一般情况下，回路是由管理员手动设置BGP规则而导致的人为失误，因此，本文在处理数据时对此类数据进行了清洗。

经过上述处理后提取出的AS网络结构的顶点和边数据，需要再分别进行去重和格式异常数据清洗操作，去除数据本身对算法的干扰，便于进一步分析该网络拓扑结构的演化情况。

2 特征与建模

针对大规模复杂网络的整体结构的研究，主要的途径是利用各种方法计算和分析其在拓扑结构上的特征，借助一些特征指标来实现对网络特征分析的量化研究。这些特征参数包括：节点总数、边总数、平均最短路径长度、网络直径、平均度、最大节点度及其对应节点、节点度大于平均度的节点所占比例、平均集聚系数、网络核数、社团数和社团模块度等。本文根据15年的IPv6 AS动态数据从以上特征维度进行了拓扑结构分析，并对富人俱乐部、自相似性、幂律分布等特性进行了验证。

2.1 静态特征分析

针对2003年5月-2018年5月的IPv6 BGP数据，统计分析AS级网络的静态特征，各维度特征值如表1、表2所示。

表1 IPv6 AS级网络静态特征汇总

表2 IPv6 AS级网络静态特征汇总

由表1可见，自2003年至2018年，IPv6 ASN总数和边总数增长迅猛，特别是2010年以后尤为突出。随着网络拓扑的复杂性越来越高，最大度值也随之逐年增大，通过分析发现自2008年之后，ASN为6 939的节点一直是持有最大度的节点，显示了其在整个网络拓扑中的重要性。进一步研究发现，AS6939是飓风电气有限责任公司(Hurricane Electric LLC)持有的主要AS编号，该公司是目前全球领先的原生IPv6互联网骨干网和主机托管服务商，根据欧洲互联网交换联盟(Euro-IX)的说法，该公司是世界上最大的交换中心参与者，运营了世界上以对等数目计算的最大IPv6网络，其中大多数是原生IPv6对等会话。该公司同时还提供了免费IPv6隧穿服务[13]，为IPv4用户或无法接入IPv6网络的用户通过隧道提供IPv6服务。

由表1可知，除网络直径、网络核数和社团数这些特征随网络规模的扩大而有较明显增长外，平均度、平均最短路径长度、平均集聚系数和社团模块度等静态特征基本稳定，只有轻微波动，说明网络容量虽然大幅度扩增，但其全局静态特征变化不大，网络环境相对稳定。从度的维度分析，由平均度、最大度的值和节点度大于平均度的占比可知，只有少数节点度很大(即“富节点”)，大部分节点的度很小，大致符合幂律分布[14]；同时，节点度大于平均度所占比例缓慢降低，说明随着网络规模的不断扩大，度值大的节点所占比例越来越小，新加入的节点更倾向于与度值大的节点连接，即这些富节点更倾向于选择彼此进行相连，则会形成“富人俱乐部”[15]。从时间的维度分析，网络特征在动态时间上表现出的稳定性很有可能是因为具有自相似性[16]。因此，针对IPv6 AS级网络拓扑可能存在的这些特性，我们将进行进一步的验证和分析。

2.2 网络特性验证

依据AS宣告的IPv6数量对AS进行排序，取前20%的AS节点及边数据重新计算全局网络特征，然后与全部AS节点的计算情况进行比对，结果如图3所示(由于平均最短路径长度、平均集聚系数、平均度等特征的变化曲线与网络直径特征的曲线走势完全一致，因此图3中只展示节点总数、边总数和网络直径三个特征曲线)，可以看出其发展趋势基本一致。由此可见，2003年-2018年IPv6 AS级网络拓扑结构随时间演化的过程符合二八定律，其中排名前20%的节点对网络结构的变化起着决定性影响，发挥着关键性作用。

图3 核心20%节点与全部节点特征对比

2.2.1幂律分布

幂律是指节点具有的连线数和节点数目的乘积为一个定值，即几何平均值是定值。也就是说节点数量和连线数量成反比，表现为在对数坐标上画出来会得到一条斜向下的直线。根据幂律分布公式：

Y=aX-b

(1)

通过对两边取以10为底的对数，得到：

logY=logaX-b=loga-blogX

(2)

令y=logY,x=logX,且c=loga为常数，公式变形为：

y=c-bx

(3)

即针对幂律公式的X和Y取双对数后，在坐标轴上应呈现线性方程图。

将上述概念应用到网络拓扑结构领域，节点间连线数量表现为各节点的度值，因此判断网络拓扑是否符合幂律分布，需对网络的度分布情况做分析验证。

我们根据2003年-2018年IPv6 AS数据，分别统计每年网络节点的度分布情况，如图4中各子图的左图所示(此处仅展示其中6年的图)，横坐标为度值，纵坐标为该度值对应的节点数在所有节点中的占比。由图可见，在整个网络结构中，度越大的节点数量越少，而度越小的节点数量越多，初步推断满足幂律分布的特性。为使结论更具准确性，进一步对该数据做线性拟合，得到图4各子图的右图，其中PE表示幂律分布的幂指数，SSR表示残差平方和的值。残差平方和是在线性模型中评估模型拟合程度的一个标准，其值越小说明线性模型拟合率越高，即越符合幂律分布特性。由图4可知，每年的IPv6 AS网络均满足幂律分布且幂指数几乎保持不变,线性模型拟合率约为0.01左右，拟合率极高。

(b) 2006年

(c) 2009年

(d) 2012年

(e) 2015年

(f) 2018年图4 2003年-2018年IPv6 AS网络度分布及幂律分布图

2.2.2富人俱乐部特性

富人俱乐部的连通性一般通过参数富人俱乐部系数φ(r/N)来度量，其计算方法是求取整个网络中前r个度值最大的节点所构成的网络中的边数L与这r个节点间两两互连能够存在的最大边数r(r-1)/2间的比值(节点总数为N)，即：

(4)

当φ(r/N)时，表示前r个度最大的节点组成的富人俱乐部为一个完全连通的子图。

利用式(4)，我们计算得到2003年-2018年IPv6 AS级网络富人俱乐部系数随时间的演化情况，结果如表3所示，其中，r在5到20每隔5取值一次，计算一次富人俱乐部系数。根据表可知：一方面，自2003年至2018年，同一r值下，AS网络每年的富人俱乐部系数值均有小幅的上下振动，但变化不大，相对稳定；另一方面，随着r取值越来越大，富人俱乐部系数会随之减小，说明富人俱乐部中的富节点数量越多其内部连通性就会相对越差，少量富节点构成的富人俱乐部其内部连通性越强，甚至接近于全连通状态。

表3 每年由度值排名前r的节点构成的富人俱乐部 的连通性统计

图5是选取2018年的IPv6 AS级真实网络数据，最大度为4 268，4个子图分别表示度大于500、300、100、50的网络结构图。可以明显看出：度大于50的节点已经屈指可数，度大的节点占整个网络规模的极少数，只有少数节点连接的边数很多，绝大多数的节点度很小；度大的节点间相互连接密切，体现了网络的富人俱乐部特性。

图5 2018年IPv6 AS网络度大的节点连接情况

2.2.3自相似性

自相似性作为分形理论的一个重要特性，越来越多地引起了广大研究者的兴趣和关注。在复杂网络中，对其自相似性的研究意义主要在于以自相似度研究为基础来准确进行链路预测、有效检测社团网络以及探索复杂网络的演化机制等[17]。目前，针对复杂网络自相似性的研究主要集中在对其节点或自身局部的相似性的研究，而在复杂网络全局拓扑特性层面上的相似度研究极少甚至没有。本文主要基于统计和分析全局的复杂网络拓扑结构特征来验证2003年-2018年动态网络的自相似程度。

1) Hurst指数。赫斯特指数(Hurst Exponent) 又称自相似参数，是用来衡量时间序列是否有长期记忆的一个指标，可以通过计算全局网络静态特征的Hurst指数来估计时间序列的自相关系数，以此评估网络拓扑结构的自相似性[18]。目前，针对时间序列的 Hurst指数的估计方法问题，国内外已经提出了 R/S分析法、绝对值法、周期图法等多种方法[19]。本文选用最常用的R/S分析方法(即重新标度的极差分析法，简称重标极差分析法)，针对AS网络主要特征维度的时间序列，计算各维度的Hurst指数。R/S分析法描述如下：

给定一时间序列xi，i=1,2,…,N，对于任意正整数n，定义累计和序列为：

(5)

计算相应的样本方差：

(6)

则R/S的统计值为：

n≥1

(7)

若存在如下关系:

• 若H=0.5，则表示该时间序列差分后的自相关系数等于0，即时间序列前后变化没有关联, 该时间序列是一个相互独立的随机序列。

• 若0.5

• 若0

2) 计算Hurst。针对2003年-2018年IPv6 AS网络的节点总数、边总数、平均度、平均最短路径长度和平均集聚系数等全局静态特征组成的时间序列，计算各单特征时间序列的Hurst指数，具体步骤如下：

(1) 将时间序列[x1,x2,…,xN]按不同片段长度分别进行均分。比如，按照下列4种片段长度划分，其中N表示整个时间序列的总长度：

① 片段长度n=N，均分后共1个片段；

(2) 计算每个片段的均值E、累积离差序列及其极差R和标准差S。

(3) 计算每个片段的R/S值，并求其在步骤(1)的不同分割方法下对应的平均值。

(4) 计算差分序列的自相关系数H：

通过上述步骤，得到R/S双对数，分别计算其斜率得到Hurst指数结果，如表4所示。结果显示，各特征维度时间序列Hurst指数均在0.8以上，表明IPv6 AS级网络在时间维度上是自相似的，过去和未来具有正相关性，且自相似程度很高，具有持续性。因此，可以通过过去和当前的网络特征情况推测未来的发展趋势和状态。

表4 各全局静态特征时间序列Hurst指数

2.3 网络建模

为了更准确地预测和评估未来IPv6 AS级互联网的发展趋势，基于对IPv6 AS网络历年来拓扑结构及特征的上述分析结果，考虑到网络演化过程中出现的变化，本文在分析网络拓扑特性的基础上，利用差分自回归移动平均模型(ARIMA(p,d,q))对主要特征的时间序列进行了建模，为未来的拓扑变化趋势作出预测。

2.3.1基础模型

ARMA(p,q)模型是目前已经相对成熟和完善的用于分析时间序列的重要模型之一，它分为两个主要步骤，其中AR代表p阶自回归过程，MA代表q阶移动平均过程，其公式如下：

Zt=φ1Zt-1+φ2Zt-2+…+φpZt-p+at-

θ1at-1-…-θqat-q

(8)

简化后得到：

φp(B)Zt=θq(B)at

(9)

式中：

φp(B)=1-φ1B-φ2B2-…-φpBpφq(B)=1-θ1B-θ2B2-…-θqBq

而ARIMA(p,d,q)模型是在ARMA模型的基础上增加了差分的操作，用来得到平稳序列，保证数据的稳定性，其中d是差分的阶数。

2.3.2平稳性检验及处理

序列平稳性是进行时间序列分析的前提条件，对时间序列进行建模的过程是基于大数定理和中心极限定理的，而其要求的样本同分布原则与序列平稳性异曲同工。为了保证结论的可靠性和准确性，必须满足此原则。例如，在输入和输出变量均平稳的情况下，可以通过t统计量来对标准化系数的显著性进行验证；而在输入和输出变量均不平稳的情况下，该标准化系数不再满足t分布，如果继续用t对其显著性进行验证和分析，则会增加拒绝原假设的概率，最终导致错误结论的产生[19]。针对非平稳时间序列直接建立回归，很容易产生伪回归。所以，在进行时间序列分析前，必须进行平稳性检验和处理，确保序列是平稳的，以避免回归分析中存在伪回归。

以节点总数的时间序列为例，自2003年至2018年IPv6 AS网络的节点总数是逐年增加的，如图6(a)所示，结合2.2.3自相似性的计算结果可知，该时间序列的自相关系数是基本稳定的，并没有跟随时间的增加而变化，由于平稳序列的自相关系数会快速衰减，因此该时间序列并非平稳序列。为得到平稳序列，对其进行1阶差分，得到图6(b)仍然不平稳，继续进行2阶差分得到图6(c)，经二次差分后的时间序列在均值和方差上趋于平稳。

(a) 时间序列原图 (b) 1阶差分处理

(c) 2阶差分处理图6 节点总数时间序列差分过程示意图

由于仅通过图像观察存在一定程度的误差性，因此本文通过单位根检验法(Augmented Dickey-Fuller,ADF)进一步验证二次差分后时间序列的平稳性。由于该序列在年周期性和长期趋势上表现明显，我们通过设置窗口为12的移动平均来处理其年周期成分，并通过差分的方法来处理其长期趋势：假设该时间序列具有单位根，为非平稳序列，对于一个平稳的时序数据，就需要在给定的置信水平上显著，以拒绝原假设。根据节点总数时间序列原数据，进行ADF检验，得到p值为0.998 6，大于99%，说明不能拒绝原假设，即时间序列为非平稳序列；经过一次差分后，进行ADF检验，得到p值为0.601 1，即60%左右，说明并不能完全拒绝原假设，即该时间序列具有一定的非平稳性；经过二次差分后，进行ADF检验，得到p值为0.001 8，小于1%，说明可以拒绝原假设，即经过二次差分的时间序列为平稳序列，满足继续建模分析的条件。

2.3.3模型识别

数据平稳后，需要对模型定阶，即确定p、q的阶数。首先检查平稳时间序列的自相关图(见图7)和偏自相关图[21](见图8)。观察图7和图8，发现置信区间被设置为默认值95%的情况下，滞后阶数在1到30之间时自相关和偏相系数都存在明显的拖尾(自相关12阶拖尾，偏自相关13阶拖尾)，符合ARMA模型的特点。利用贝叶斯信息准则(BIC)来衡量模型拟合度，获得最优ARIMA模型。通过计算得到BIC为-1 544.227 2，p为2，q为4时模型拟合度最高，结合差分情况，即完整模型为ARIMA(2,2,4)。

图8 节点总数时间序列二次差分后的偏自相关图

2.3.4样本拟合

确定模型后，对时间序列进行预测。由于模型的拟合过程中进行了差分等预处理，所以通过模型预测得出的值也需要通过反向处理进行数据还原，对还原后的预测序列做可视化处理[22]，得到图9所示结果(其中，黑线表示模型预测结果，灰线表示原时间序列)。观察图9可见，模型拟合率很高，预测效果较为可观。依据同样的方法，对其他主要特征进行了同样的建模与预测，得到图10。

图9 节点总数时间序列与模型预测结果对比图

(a) 边总数 (b) 最大度

(c) 节点度大于平均度占比(d) 平均最短路径长度

(e) 最大核数 (f) 社团数图10 主要特征的时间序列与模型预测结果对比图

2.4 分析与评价

根据IPv6 AS级网络动态数据的静态特征、网络特性及建模结果可知，2003年-2018年网络节点总数、边总数持续增长，最大度呈明显上升趋势，网络直径、网络核数和社团数呈现缓慢增长趋势，平均节点度、平均集聚系数、平均最短路径长度和社团模块度均相对稳定，振幅很小，随着时间的推移几乎不变；节点度大于平均度所占比例随时间推移呈现缓慢降低趋势，说明新增节点更倾向于与度值大的节点建立连接关系；动态网络的节点度分布情况满足幂律分布，说明大多数节点的度值很小，只有少数的节点度值较大；动态网络存在富人俱乐部现象，说明度较大的节点更倾向于与其他度较大的节点建立连接，形成富节点组成的富人俱乐部；动态网络具有自相似性，说明过去的网络特征与未来的网络特征成正相关关系，可以通过过去的网络状态预测未来的网络发展趋势；利用ARIMA模型可以很好地拟合IPv6 AS级网络状态，通过该模型可以较准确地预测未来网络拓扑的变化情况，对未来网络的发展趋势分析有指导作用和长远意义。

3 结 语

本文通过对2003年-2018年15年间的IPv6 AS级网络静态特征及网络特性的分析，观察网络规模扩增过程中各特征量的变化，总结IPv6 AS级网络拓扑的演化规律，并通过建模分析网络拓扑状态，预测发展趋势。根据分析结果可以推测，未来IPv6 AS级网络规模总体会继续呈上升趋势增长，但平均度、平均最短路径长度、平均集聚系数和社团模块度等全局静态特征变化不大，未来将依然处于相对稳定的状态，这符合IPv6 AS级网络特性和发展规律。