一种网络流量估计方法技术领域
本发明属于计算机网络通信技术领域，特别涉及其中的网络流量估计方法。
背景技术
随着互联网技术高速发展，网络类型多样化、规模大型化、结构复杂化已成为一种趋
势，网络中各种性能参数的变化也越来越复杂。通信网络中各种性能参数的测量是网络行
为学的基本的研究内容，也是网络的优化设计和管理、网络监测和评估、流量异常监测和
分析等工作的基础。为了更好地对网络进行设计、监控、管理，需要了解和感知这些网络
状态参数的动态变化，包括链路延时、网络拓扑结构、网络OD对流量等参数。网络操作
人员或管理人员通过了解网络流量在全网络中的分布情况，可充分利用网络宽带资源提高
网络资源的利用率。流量矩阵是网络流量在全网络中的分布的具体表示，通过实时获取准
确的流量矩阵，网络操作人员或网络管理人员能更好的制定出路由策略，使网络总流量负
载均衡，而且能实时检测网络中流量的异常，防止外来攻击进而保证网络的安全。
实际网络情况很复杂，直接测试流量矩阵很困难，最终可能覆盖不到需要测量的网络
节点。但是网络中链路带宽的利用率、网络拓扑和路由配置状况等信息是可以获取的，通
过建立估计模型或算法来获得流量矩阵是一种通用的做法。
流量矩阵(Traffic Matrix,TM)是网络流量的具体描述，由源节点到目的节点的流量
(Origin-Destination,OD)组成，是网络中指定节点对间流量大小的具体值，给出了网络流量
在全网中各个OD对间流量的分布情况。网络流量工程需要基于流量矩阵进行分析，流量
矩阵也是容量规划、可靠性分析，异常检测等网络管理方面的重要输入项。因此，在实际
网络中具有重要的意义。求解流量矩阵的方法有很多，主要如下两种：
主成分分析(Principal Component Analysis,PCA)算法需要一段时间的所有OD对流量
值作为先验信息，基于主成分分析方法和矩阵理论分析，通过使用OD对流量的主要成分(特
征流向量)来表示OD对流量集合，将流量矩阵估计问题转化为气逆特征流向量的问题。由
于线性的PCA算法计算量大，并且需要的存储空间要求比较大，它忽略了数据在高维上的
关系。
空时压缩感知算法是利用直接测量的流量矩阵中由于链路或者是路由器损坏的原因丢
失的一些OD对流量值进行插值运算，充分利用稀疏正则化矩阵分解及空间与时间域的约
束进行流量矩阵估算，而后利用KNN算法进行插值对估算结果进行修正的一种算法,得到
最终的估算结果。在实际的流量矩阵的估计中，由于要估计的网络节点很多，时间点比较
长，所以TM非常大，不进行降维计算会出现计算量很大，计算很复杂，需要用于计算的
环境配置要求很高。
发明内容
本发明的目的是为了解决现有的网络流量估计方法存在的上述问题，本发明提出了一
种网络流量估计方法。
本发明的技术方案为：一种网络流量估计方法，包括如下步骤：
步骤1：获取路由矩阵及链路负载；
步骤2：剔除步骤1获取的链路负载中部分已知OD对流量所占负载，并通过网络总节
点及已知流量的OD对标记未知流量的OD对；
步骤3：利用LP算法估算步骤2获得的未知的OD对流量；
步骤4：用已知流量的OD对与步骤3估算出来的OD对流量对总的OD对进行重构；
步骤5：基于KPCA分簇对步骤2获得的未知流量的OD对进行插值；
步骤6：用步骤5插值出来的OD对流量对步骤3通过LP估算出来的OD对流量进行
校正，得到最优解。
进一步的步骤5具体包括如下分步骤：
步骤51：由流量矩阵TM确定丢失模型LM；
步骤52：根据丢失模型LM，从流量矩阵TM中取出固定个数OD对的小矩阵，其中，
无缺失OD对小矩阵与有缺失OD对小矩阵的大小一致;
步骤53：对取出来的小矩阵用KPCA算法映射到特征空间;
步骤54：对当前的无缺失OD对的小矩阵进行分簇;
步骤55：对当前的缺失OD对小矩阵选择最优的簇作为参照进行插值，得到最优的插
值结果。
本发明的有益效果：本发明的网络流量估计方法，实现了在已知部分OD对流量的情
况下，对IP网络流量矩阵进行精确的估计，采用了KPCA分簇对缺失OD对进行插值，与
传统的PCA算法相比可以降低计算复杂度，用于处理大规模的TM；而LP算法能满足链路
流量的约束，特别是在链路条数比较多，网络节点数比较少情况下，能够精确估计，降低
均方根相对误差。本发明的方法通过LP算法估算出来的值经过基于KPCA插值估算出来的
结果校正以后，使估算的OD对流量更加接近真实的OD对流量。
附图说明
图1为本发明的网络流量估计方法的流程示意图。
图2为本发明中用基于KPCA分簇方法对OD对流量进行插值方法的流程示意图。
具体实施方式
下面结合附图和具体实施例对本发明做进一步的说明。
本发明的网络流量估计方法的流程示意图如图1所示，具体步骤如下:
步骤1：获取路由矩阵A及链路负载Ytotal；
路由矩阵的获取主要涉及基于metric域内路由查找算法(用最短路径算法—Dijkstra算
法)，获得流量矩阵估计方程：
Ytotal=AXtotal (1)
其中，Xtotal为总的OD对流量。
步骤2：剔除步骤1获取的链路负载Ytotal中部分已知OD对流量Xpart所占的负载，获得
未知的OD对在链路中的负载Yleft，并通过网络总节点及已知流量的OD对标记未知流量的
OD对；
由于已经知道部分OD对的流量，在用线性规划算法(Linear Programming,LP)计算的时
候为了减少计算量及提高计算精确度，对已经知道的部分OD对流量不再进行估计，如下
所示：
X=XtotalΘXpart (2)
其中，Xpart表示部分已知OD对流量，X为未知OD对流量，Θ表示从总Xtotal中排除
Xpart项。
于是可以根据总链路负载中已知OD对流量所占负载，得到未知OD对流量在链路中
的负载，表达式如下：
Yleft=Ytotal-A*Xpart (3)
其中，Ytotal为总的链路负载，Yleft为剔除部分已知OD对流量后的链路负载。
步骤3：利用LP估算步骤2剔除部分已知OD对流量后的未知的OD对流量。
首先由路由矩阵A及剩余链路负载Yleft获得OD对流量的上确界和下确界。由于链路流
量为Yleft，未知OD对流量为X，于是OD对流量估计的系统方程可以定义为：
Yleft＝AX (4)
估计上确界Xup，目标函数为：
max X up subject to Y left = A X up X up < ub - - - ( 5 ) ]]>
其中，ub为X的上限约束，为向量，在没有其它约束信息的时候，采用ub=max(Yleft)，
也就是当前的最大链路值。在实际应用中可以加更精细的约束，比如部分OD对流量可以
获得更小范围的上限值等其它约束条件。
估计下确界为Xlow的目标函数：
max X low subject to Y left = AX low lb < X low - - - ( 6 ) ]]>
其中，lb为X的下限约束，在没有其它约束信息条件下，常采用lb=0。当然也可
以像上限约束ub一样，采用其它的约束信息，以提高求解精度，那么线性规划的解：
X = X up + X low 2 - - - ( 7 ) ]]>
该算法的核心为线性规划算法的选择，线性规划算法选择主要能够满足大尺度大规模
的方程求解，此处采用文献：“Yin Zhang.Solving Large-Scale Linear Programs by Interior-Point
Methods Under the MATLAB Environment.Department of Mathematics and Statistics University
of Maryland Baltimore County，Technical Report TR96-01,February，1996”中提到用于处理大
规模病态线性方程的primal-dual Infeasible-Interior-Point algorithm(原-对偶不可行内点算法)
进行估算。
步骤4：用已知OD对流量与步骤3估算出来的未知OD对流量对总的TM进行重构；
在LP估算出未知OD对流量之后，就与原来的已知OD对合并重构为最初总的流量矩
阵TM，如表1所示：
表1
已知OD
OD已知
OD缺失
OD已知
...
OD已知
OD缺失
OD已知
估算OD
OD缺失
OD估算
OD缺失
...
OD缺失
OD估算
OD缺失
总OD
OD已知
OD估算
OD已知
...
OD已知
OD估算
OD已知
步骤5：基于KPCA分簇对缺失OD对进行插值，具体过程如图2所示，展开如下：
步骤51：由流量矩阵TM确定丢失模型LM。
由于原来的TM部分OD对是知道流量的，可以把它看作一个丢失模型LM(Lost Model)，
LM是一个N2×m大小的矩阵，其中N表示网络节点的个数，N2表示OD对的个数，m表
示路由中总的链路条数。表示形式如下：
其中，X(i,j)表示经过第j条链路的第i个OD对，这可以用来确定小矩阵内需要对哪
些缺失OD对进行插值。
步骤52：根据LM，从TM中取固定个数OD对的小矩阵，其中，无缺失OD对小矩
阵与有缺失OD对小矩阵的大小一致；
根据LM从TM中随意取出l个小矩阵f，大小为(w×h，其中w与h可以随意定义大
小，但不能太大，分别表示取出矩阵的行数与列数，即OD对的个数与链路的条数)，这些
小矩阵不包含有缺失的OD对（由LM中得知），再取一个有缺失OD对的小矩阵(w×h)，
在OD对流量不知道的地方记作在OD对流量知道的地方记作Ω。
步骤53：对取出来的小矩阵用KPCA算法映射到特征空间。
对于上面的取出来的小矩阵fi(其中，i=1,2,3...,l)，可以获取两个向量
其中，xi与yi表示OD对的流量，不知道流量的OD对初始化
为0，表示不知道OD对流量的集合，表示知道OD对流量的集合，NΩ表示知
道OD对流量的个数，为知道的OD对流量，为不知道的OD对
流量，初始化为0。通过KPCA把它映射到特征空间φx(xi)与φy(yi)，具体如下：
TM数据X(N2×m)通过点积的形式把它从数据空间(实际OD对流量，低维空间)映射
到特征空间(数学上，高维空间)φ(x1),φ(x2),......,φ(xl)，映射定义为φ，如下面式子所示：
其中，表示映射到特征空间的过程，ξ表示映射到的特征空间。由于非线性映射函
数φ非常复杂，实现起来很困难，在运算中用核函数代替会简单很多。核函数如下，可以表
示为特征空间的点积的形式：
K(xi,xj)=φ(xi)·φ(xj) (10)
核函数使用高斯函数，则(10)可以写成下面的式子:
K ( x i , x j ) = exp ( - | | x i - x j | | 2 σ 2 ) - - - ( 11 ) ]]>
其中，σ2为高斯核函数的参数，通过这个公式计算出的核函数将在后面用到。
步骤54：对当前的无缺失OD对的小矩阵进行分簇。
由于非线性特征空间有许多结构，要选择最优结构的特征空间作为参照目标用于后面
的插值，所以必须进行簇分类。
先定义一个向量为 φ i = φ x ( x i ) φ y ( y i ) , ]]>向量为步骤53中OD对流量映射到特征空间的向量，
后面将用核函数代替。使用下面的分簇准则公式来进行簇分类，分簇准则公式为：
E = Σ k = 1 M Σ j = 1 M k { | | x j k - x ‾ j k | | 2 + | | y j k - y ‾ j k | | 2 } - - - ( 12 ) ]]>
其中，||·||表示求二范数，与分别是属于簇k的xi和yi，即为当前簇的OD对流量,
Mk是簇k的元素个数,K为总的簇个数，与分别表示当前分类小矩阵的OD对流量。
给定 φ j k = φ x ( x j k ) φ y ( y j k ) ]]>和 φ ~ j k = ζ j k ξ j k , ]]>满足：
φ ~ j k = U k U k ′ ( φ j k - φ ‾ k ) + φ ‾ k - - - ( 13 ) ]]>
其中，Uk′为Uk的转置，和满足和与分
别是当前分类小矩阵的已知OD对流量与未知OD对流量映射到特征空间后的值。在式(13)
中
是ΞkHkHkΞk′的特征向量矩阵，其中，
Ik是一个Mk×Mk的单位矩阵，1k=[1,1,...,1]'是一个
Mk×1的向量的特征向量矩阵，ΞkHk满足下面的奇异值分解：
Ξ k H k ≅ U k Λ k V k ′ - - - ( 14 ) ]]>
其中，Λk与Vk分别表示特征值及向量矩阵，Vk′为Vk的转置，UkUk还可以写成：
U k ≅ Ξ k H k V k Λ - 1 - - - ( 15 ) ]]>
是簇k的一个中心向量，可以把式(13)写成：
φ ~ j k ≅ Ξ k T k Ξ k ′ φ j k - 1 M k Ξ k ( T k Ξ k ′ Ξ k - I k ) 1 k - - - ( 16 ) ]]>
其中，最终式(12)中的可以由如下推断得到：
| | x j k - x ‾ j k | | = - σ x 2 log { φ x ( x j k ) ′ φ x ( x ‾ j k ) } = - σ x 2 log { ( 1 - ψ j k ′ Ξ x k ′ Ξ x k ψ j k + 2 φ x ( x j k ) ′ Ξ x k ψ j k ) 2 } - - - ( 17 ) ]]>
其中，是高斯核函数的参数。
可以下式得到：
| | y j k - y ‾ j k | | ≅ - σ y 2 log { ( 1 - ψ j k ′ Ξ y k ′ Ξ y k ψ j k + 2 φ y ( y j k ) ′ Ξ y k ψ j k ) 2 } - - - ( 18 ) ]]>
其中，是高斯核函数的参数。
使用上面的方法，可以对提取的OD对小矩阵进行分簇，如果当前小矩阵使式(12)的值
最小，它就属于当前簇。
分簇的步骤可以归纳为：
步骤541：定义分簇准则公式(12)，通过它划分OD对小矩阵属于哪个簇。
步骤542：由KPCA映射的特征空间向量满足式(13)。
步骤543：由簇K的中心向量，及式(14)、(15)，可以把式(13)重新写为式(16)。
步骤544：由式(17)和式(18)，可以求出式(12)中的与由此就可以计
算出式(12)中的E。
步骤545：当使式(12)最小时，就属于当前簇。
步骤55：对当前的缺失OD对小矩阵选择最优的簇作为参照进行插值。
在步骤54中已经把提取的OD对小矩阵进行了分簇，下一步要选取最优的特征空间作
为缺失OD块的参考簇，下面就是具体选择最优参考簇及插值的操作：
根据步骤53把目标矩阵f映射到特征空间φf， φ f = φ x ( x ) φ y ( y ) , ]]>其中，φx(x)与φy(y)分
别为目标矩阵f中知道的OD对流量与不知道的OD对流量在特征空间的映射，由前面的分
簇公式计算出一个新的向量 φ f k ^ = ζ k ξ k ]]>满足下面的两个约束来获取估计未知的OD对流量。
①φy(y)是知道的OD对流量的特征空间，它固定为空间向量
②在特征空间，目标向量在非线性特征空间通过簇k的特征向量组成。
因此，满足
φ ~ f k = U k U k ′ ( φ f k - φ k ‾ ) + φ k ‾ - - - ( 19 ) ]]>
由式(16)可以把式(19)写成
φ ~ f k ≅ Ξ k T k Ξ k ′ φ f k - 1 M k Ξ k ( T k Ξ k ′ Ξ k - I k ) 1 k - - - ( 20 ) ]]>
如果定义ξk=φy(y)，那么由式(20)可以得到：
ζ t k ≅ Ξ x k T k ( Ξ x k ′ ζ t - 1 k + Ξ y k ′ φ y ( y ) ) - 1 M k Ξ x k { T k ( Ξ x k ′ Ξ x k + Ξ y k ′ Ξ y k ) - I k } 1 k - - - ( 21 ) ]]>
其中，是一个任意的初值向量，可以定义一个初始的OD对流量。
但是现在不知道插值所需要的最优特征空间，为了求得它，用另外一个准则公式代替
式(12)，如下所示：
E k ~ = | | y - y k ^ | | 2 wh - N Ω - - - ( 22 ) ]]>
其中，满足：
φ y ( y k ^ ) ≅ Ξ y k T k ( Ξ x k ′ ζ t - 1 k + Ξ y k ′ φ y ( y ) ) - 1 M k Ξ y k { T k ( Ξ x k ′ Ξ x k + Ξ y k ′ Ξ y k ) - I k } 1 k - - - ( 23 ) ]]>
与式(18)类似，能够通过下面的式子得到：
E k ~ ≅ - σ y 2 log { ( 1 - ψ k ′ Ξ y k ′ Ξ y k ψ k + 2 φ y ( y ) ′ Ξ y k ψ k ) 2 } wh - N Ω - - - ( 24 ) ]]>
其中，wh-NΩ为小矩阵中不知道流量的OD对的个数，ψk为：
ψ k = T k ( Ξ x k ′ ζ k + Ξ y k ′ φ y ( y ) ) - 1 M k { T k ( Ξ x k ′ Ξ x k + Ξ y k ′ Ξ y k ) - I k } 1 k - - - ( 25 ) ]]>
通过式(24)可以选择得到最优的簇kopt，使在最小的情况下得到φf。
通过上面的方法在簇kopt下得到在满足的条件下，把映射到空间的向量
对应的数据当作缺失OD对流量的估计值，为最优插值的估计值。
可以将选择最优簇归纳为如下分步骤：
步骤551：在特征空间，建立目标向量，如式(19)。
步骤552：由式(20)，可以得到插值OD对流量目标向量式(21)。
步骤553：建立最优准则公式(22)。
步骤554：求解使最优准则最小时的的簇k即为最优簇。
步骤555：计算出插值OD对流量目标向量对应的就是最优插值的估计值。
步骤6：用插值出来的OD对流量对LP估算出来的OD对流量进行校正，得到最优解。
通过校正LP估算出来的缺失的部分OD对流量X，得到最优的解，由于计算出
来的值不止一个，又由于LP估算出来的值满足链路流量的约束，所以只需使
其中φ(X)为LP算法估算的未知OD对流量X
映射到特征空间的值，从而得到对应的最优解yopt，这样对LP算法估算出来的值进行了校
正优化（它满足LP约束，同时又满足插值要求），即完成了校正，计算出了最终的结果。
综上可以看出，本发明的方法实现了在已知部分OD对流量的情况下，对IP网络流量
矩阵进行精确的估计，具有如下优点：
（1）KPCA算法是一种基于非线性特征空间的算法，结合它对OD对进行分簇，然
后对缺失OD对流量进行插值，适用于处理大规模的TM有缺失OD对流量的情况。
（2）LP算法能满足链路流量的约束，特别是在链路条数比较多，网络节点数比较少
情况下，本身估计就比较精确了，RMSR(Mean Square Of Relative Error，均方根相对误差)
能达到15%以下；通过LP估算出来的OD流量值经过基于KPCA插值估算出来的结果校正
以后，使更加接近真实的OD对流量值。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原
理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术
人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和
组合，这些变形和组合仍然在本发明的保护范围内。