近年来，随着信息技术的进步，网络的资源分配成为一个重要的研究课题.当网络过载或拥塞时，完善的资源分配机制能保证重要或高优先级的业务量不会发生延迟或丢弃，同时保证网络的高效运行^[1].网络流量预测技术的发展与成熟使得建立一种基于流量精确预测的动态资源分配成为可能^[2].网络流量预测根据应用场合不同，通常可分为短期预测、中期预测及长期预测^[3-4].一般而言，长期预测通常是以比较大的颗粒度即月、日的长时间段内的历史数据为参考来分析、建模、预测未来时刻流量变化，通常更注重趋势上的准确，而不是预测值的绝对准确.短期预测则更要求实时性，要求对未来秒级甚至更小尺度的网络流量进行预测.而中期预测一般是依据几分钟至几小时的历史数据，预测出未来几分钟或几小时内网络流量的变化，同时也要求具有较高的预测准确度.相对于长期预测，中短期预测更具有挑战性和难度.

目前，网络流量的中短期预测模型包括线性模型与非线性模型.线性预测模型有自回归滑动平均(ARMA, auto regressive moving average)模型^[5]、自回归积分滑动平均(ARIMA, auto regressive integrated moving average)模型^[6-7]及差分自回归求和滑动平均模型^[8-9]等.姜明等^[10]利用以上线性模型对不同时间尺度下的网络流量预测进行了研究，仿真结果表明线性模型的预测精度有限.同时，随着现今网络模型的动态化和复杂化，网络流量特性已经偏离了相关学者认为的泊松或马尔可夫分布^[11]，因此线性模型无法反映网络流量表现出的新特性.而非线性预测模型的典型代表包括支持向量机^[12-13]、最小二乘支持向量机^[14-15](LSSVM, least square support vector machine)、人工神经网络^[16-18]及灰色模型^[19]等.相对于线性模型，非线性模型的预测精度有了一定程度的提高，但也存在着各自的缺陷.同时，进一步的研究也指出非线性模型对于时间序列中线性成分的预测精度有限.对于1个时间序列，因为无法确定其是否包含线性或者非线性成分，所以采用单一的线性或非线性预测模型是不适合的.最合适的方法是针对时间序列中的线性成分采用线性预测模型，而非线性成分则采用非线性预测模型.

笔者通过BDS(Brock-Dechert-Scheinkman)统计量检验方法对实际采集的中期网络流量时间序列进行分析表明，网络中期流量既含有线性成分又含有非线性成分.因此，利用ARIMA模型对原始网络中期流量进行预测，得到网络中期流量的线性趋势；然后利用具有良好非线性逼近能力的高斯过程回归模型对预测误差进行预测，并通过人工蜂群(ABC, artificial bee colony)算法对高斯过程回归模型的参数进行优化，以提高预测精度；最后将ARIMA模型的线性预测值与高斯过程回归模型的非线性预测误差值进行相加而得到最终的预测值.通过仿真对比表明，所提出的网络流量中期预测方法具有较高的预测精度.

1 网络流量的线性与非线性判定

笔者利用BDS统计量检验法进行网络流量的线性/非线性判定. BDS统计量检验法是基于G-P相关积分创建的一种统计量，用于判定时间序列的非线性特征.对于长度为$n$的时间序列，构造嵌入维数为$m$、时延为$τ$的嵌入向量

$\begin{eqnarray} \boldsymbol{x}^{m}_{t}=(x_{1τ},x_{2τ},…,x_{(m-1)τ}) \end{eqnarray}$

(1)

其相关积分为

$\begin{eqnarray} C(m,n,r)=\frac{2}{M(M-1)}∑\limits_{t＜s}H(r-‖\boldsymbol{x}^{m}_{t}-\boldsymbol{x}^{m}_{x}‖) \end{eqnarray}$

(2)

其中：$M$=$n$-($m$-1)；$r$为包含时间序列点的$m$维球的半径；$s$为参考距离；$H$为Heaviside函数，即

$\begin{eqnarray} H(z)=\left\{\begin{array}{c} 0,&z≤0\\ 1,&z＞1\end{array} \right. \end{eqnarray}$

(3)

用BDS进行检查前，需要消除原时间序列线性相关的成分，通常对原时间序列拟合自回归模型AR($p$)，在寻找到合适的阶数$p$后，计算AR($p$)的残差序列并对该残差序列进行BDS检验. BDS检验的零假设：该残差序列是独立同分布的.如果结果拒绝零假设，则意味着原时间序列在某个显著水平下是非线性的.

采集了250组辽宁联通169核心网络接入路由器的网络中期流量数据，数据采集尺度为10 min，数据单位为MB. 250组网络流量数据如图 1所示.

令AR($p$)模型的$p$为2，得到残差序列，并对残差序列进行BDS检验，取嵌入维数$m$为2，取$r$为序列标准差的0.5、0.75、1、1.25及1.5倍数，检验结果如表 1所示.

标准正态分布95%的临界值为1.96，从表 1中可知，网络流量时间序列的BDS统计量检验结果远大于1.96，这就表明网络流量时间序列的残差序列不是独立同分布的时间序列，具有强烈的非线性特征.故而网络流量时间序列既含有线性又含有非线性成分.因此对残差序列需要用非线性预测模型进行预测.基于以上分析，笔者提出了一种线性与非线性模型相结合的预测方法.

2 ARIMA预测模型

ARIMA模型的基本思想是通过多次差分使非平稳序列平稳化，差分次数是参数$d$，然后用以$p$、$q$为参数的ARMA模型对该平稳序列建模，之后经过反变换得到原序列.以$p$、$d$、$q$为参数的ARIMA模型预测方程可表示为

$\begin{eqnarray} &y_{t}=θ_{0}+φ_{1}y_{t-1}+φ_{2}y_{t-2}+…+φ_{p}y_{t-p}+\\ &ε_{t}-θ_{1}ε_{t-1}-θ_{2}ε_{t-2}-…-θ_{q}ε_{t-q} \end{eqnarray}$

(4)

其中：$y_{t}$为时间序列的样本值，$φ_{i}$($i$=1, 2, …, $p$)和$θ_{i}$($i$=1, 2, …, $q$)为模型参数，$ε_{t}$为独立正态分布的白噪声.

平稳化时间序列后，首先计算原始序列的自相关函数和偏自相关函数.对于时间序列$y_{t}$，有自协方差为

$\begin{eqnarray} γ_{k}=\frac{1}{N}∑\limits^{N-k}_{j=1}y_{k}y_{j+k} \end{eqnarray}$

(5)

自相关函数为

$\begin{eqnarray} ρ=\frac{γ_{k}}{γ_{0}} \end{eqnarray}$

(6)

偏自相关函数为

$\begin{eqnarray} &α_{11}=ρ_{1}\\ &α_{k+1,k+1}=\left(ρ_{k+1}-∑ρ_{k+1-j}α_{kj}\right)\left(1-∑\limits^{k}_{j=1}ρ_{j}α_{kj}\right)^{-1}\\ &α_{k+1,j}=α_{kj}-α_{k+1,k+1}α_{k,k-j+1} \end{eqnarray}$

(7)

可通过$ρ_{k}$、$α_{k}$的截尾性初步确定模型的阶数.时间序列的参数辨识可通过最小二乘估计得到，即估计参数$φ_{1}$, $φ_{2}$, …, $φ_{p}$, $θ_{1}$, $θ_{2}$, …, $θ_{q}$使式(8)最小.

$\begin{eqnarray} ∑\limits^{N}_{t=1}α^{2}_{t}=∑\limits^{N}_{t=1}(θ^{-1}_{q}(Z)φ_{p}(Z)∇^{d}y_{t})^{2} \end{eqnarray}$

(8)

对不同的$p$、$d$、$q$参数进行组合，通过赤池信息量准则得到最优的参数模型.赤池信息量准则定义为

$\begin{eqnarray} A=-2\text{ln} L+2n \end{eqnarray}$

(9)

其中：$L$为模型的极大似然参数，$n$为模型的独立参数.

ARIMA模型的实质就是差分运算与ARMA模型的结合，即通过适当阶数的差分操作后，可实现任何非平稳时间序列的平稳，适合于线性时间序列的预测.

3 ABC算法优化的高斯过程回归

高斯过程回归是近年发展起来的一种机器学习回归方法，它有着严格的统计学习理论基础，方便预测未来事物的发展.其对处理高维数、小样本、非线性等复杂问题具有很好的适应性，且泛化能力强^[20].它的基本原理是：假设给定训练数据集$D$={($\boldsymbol{x}_{i}$, $y_{i}$)|$i$=1, 2, …, $n$}，其中$\boldsymbol{x}_{i}$表示训练数据集$D$中的第$i$个输入向量，$y_{i}$表示训练数据集$D$中的第$i$个目标输出，$n$表示训练数据集中样本的个数.

假定$f$是一个高斯过程，即$f$~GP($m$, $k$)，$f$是一个以$m$为均值函数、$k$为协方差函数的高斯过程.高斯过程是一个随机过程，与高斯分布类似，高斯过程完全由其均值函数与协方差函数确定.根据高斯过程的定义可知，$f$($\boldsymbol{x}_{1}$)、$f$($\boldsymbol{x}_{2}$)、…、$f$($\boldsymbol{x}_{n}$)服从多元高斯分布，且该多元高斯分布的均值向量为$m$($\boldsymbol{x}_{i}$)，协方差矩阵为$\boldsymbol{K}$，因此有

$\begin{eqnarray} f(\boldsymbol{x}_{i})～N[m(\boldsymbol{x}_{i}),\boldsymbol{K}],i=1,2,…,n \end{eqnarray}$

(10)

$\begin{eqnarray} D:y_{i}=f(\boldsymbol{x}_{i}),i=1,2,…,n \end{eqnarray}$

(11)

实际目标输出$\boldsymbol{y}$往往会包含一些噪声：

$\begin{eqnarray} \boldsymbol{y}=f(\boldsymbol{x}_{i})+ε_{i} \end{eqnarray}$

(12)

其中$ε$~$N$(0, $σ^{2}_{n}$).于是问题转化为已经观测到训练数据$D$:$y_{i}$=$f$($\boldsymbol{x}_{i}$)+$ε_{i}$($i$=1, 2, …, $n$)，需要在测试数据集$D_{*}$={($\boldsymbol{x}_{i}$, $y_{i}$)|$i$=$n$+1, $n$+2, …, $n$+$n_{*}$}预测对应的输出值$f_{*}$.训练数据集的输出向量$\boldsymbol{y}$和测试数据集的预测值$f_{*}$的多元高斯分布为

$\begin{eqnarray} \left[\begin{array}{c}\boldsymbol{y}\\f_{*}\end{array}\right]～N\left(0, \left[\begin{array}{c}\boldsymbol{K}+σ^{2}_{n}\boldsymbol{I} &\boldsymbol{K}^{\text{T}}_{*}\\ \boldsymbol{K}_{*}&\boldsymbol{K}_{**}\end{array}\right]\right) \end{eqnarray}$

(13)

其中

$\begin{eqnarray} &\boldsymbol{K}=\left[\begin{array}{c} k(\boldsymbol{x}_{1},\boldsymbol{x}_{1})&k(\boldsymbol{x}_{1},\boldsymbol{x}_{2})&…&k(\boldsymbol{x}_{1},\boldsymbol{x}_{n})\\ k(\boldsymbol{x}_{2},\boldsymbol{x}_{1})&k(\boldsymbol{x}_{2},\boldsymbol{x}_{2})&…&k(\boldsymbol{x}_{2},\boldsymbol{x}_{n})\\ ⋮&⋮&⋱&⋮\\ k(\boldsymbol{x}_{n},\boldsymbol{x}_{1})&k(\boldsymbol{x}_{n},\boldsymbol{x}_{2})&…&k(\boldsymbol{x}_{n},\boldsymbol{x}_{n})\end{array}\right]\\ &\boldsymbol{K}_{*}=[k(\boldsymbol{x}_{*}, \boldsymbol{x}_{1}), k(\boldsymbol{x}_{*}, \boldsymbol{x}_{2}), …, k(\boldsymbol{x}_{*}, \boldsymbol{x}_{n})]\\ &\boldsymbol{K}_{**}=k(\boldsymbol{x}_{*}, \boldsymbol{x}_{*}) \end{eqnarray}$

根据多元高斯分布的条件分布形式，可得出高斯过程预测方程的关键式为

$\begin{eqnarray} f_{*}|\boldsymbol{X},\boldsymbol{y},\boldsymbol{X}_{*}～\text{GP}[m(\boldsymbol{x}_{*}),\text{cov}(f_{*})] \end{eqnarray}$

(14)

其中：矩阵$\boldsymbol{X}$由训练数据输入$\boldsymbol{x}_{i}$的列向量组成，矩阵$\boldsymbol{X}_{*}$由测试集的输入$\boldsymbol{x}_{i*}$的列向量组成.

$\begin{eqnarray} m(\boldsymbol{x}_{*})=\boldsymbol{K}_{*}(\boldsymbol{K}+σ^{2}_{n}\boldsymbol{I})^{-1}\boldsymbol{y} \end{eqnarray}$

(15)

$\begin{eqnarray} \text{cov} (f_{*})=\boldsymbol{K}_{**}-\boldsymbol{K}_{*}(\boldsymbol{K}+σ^{2}_{n}\boldsymbol{I})^{-1}\boldsymbol{K}^{\text{T}}_{*} \end{eqnarray}$

(16)

高斯过程回归中协方差函数是一个满足Mercer条件的对称函数，其在有限输入集上正定，因此协方差函数等价于核函数，将式(15)改写为

$\begin{eqnarray} m(\boldsymbol{x}_{*})=∑\limits^{n}_{i=1}\mathbf{α}_{i}\boldsymbol{K}(\boldsymbol{x}_{i},\boldsymbol{x}_{*}) \end{eqnarray}$

(17)

其中

$\begin{eqnarray} \mathbf{α}_{i}=[\boldsymbol{K}(\boldsymbol{X}, \boldsymbol{X})+δ^{2}_{n}\boldsymbol{I}]^{-1}\boldsymbol{y} \end{eqnarray}$

这里选择平方指数函数作为核函数，即

$\begin{eqnarray} \boldsymbol{K}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})=δ^{2}_{f}\text{exp} \left(-\frac{1}{2l^{2}}(\boldsymbol{x}_{i}-\boldsymbol{x}_{j})^{2}\right)+δ^{2}_{n}δ_{ij} \end{eqnarray}$

(18)

其中：$l$为关联性测定超参数，$δ^{2}_{f}$为核函数的信号方差，$δ^{2}_{n}$为噪声方差，$δ_{ij}$为Kronecker符号.一般将$l$、$δ^{2}_{f}$、$δ^{2}_{n}$称为超参数.超参数直接决定了高斯过程回归的预测性能.目前一般通过共轭梯度法来确定高斯过程回归模型的参数，但是共轭梯度法通过梯度下降求取最优解，容易陷入局部最优，同时算法效果和收敛性依赖初始值.

ABC算法来源于蜜蜂采蜜的行为^[21]. ABC算法相对于遗传算法、差分进化算法以及粒子群优化算法具有极强的竞争力^[22].基于此，采用ABC算法对高斯过程回归的超参数进行优化.

ABC算法中蜜源的位置抽象为解空间内的点，蜜源($i$=1, 2, …, $S$)的质量为解的适应度.设求解问题的维数为$d$，对于第$t$次迭代蜜源$i$的位置可表示为$\boldsymbol{X}^{t}_{i}$=[$x^{t}_{i1}$ $x^{t}_{i2}$ $x^{t}_{id}$]，$x_{id}$∈($L_{d}$, $U_{d}$)，其中$L_{d}$为搜索空间下限，$U_{d}$为搜索空间上限，蜜源的位置按式(19)在搜索空间随机产生.

$\begin{eqnarray} x_{id}=L_{d}+\text{rand}(0,1)(U_{d}-L_{d}) \end{eqnarray}$

(19)

搜索开始阶段，引领蜂在蜜源周围根据式(20)搜索而产生一个新的蜜源.

$\begin{eqnarray} v_{id}=x_{id}+φ(x_{id}-x_{jd}) \end{eqnarray}$

(20)

其中：$j$∈{1, 2, …, $S$}，$j$≠$i$，代表在蜜源中随机选择一个不同于$i$的蜜源；$φ$∈[-1, 1].当新蜜源$\boldsymbol{V}_{i}$=[$\begin{array}{c}v_{i1}&v_{i2}&…&v_{id}\end{array}$]的适应度优于$\boldsymbol{X}_{i}$时，利用贪婪算法将$\boldsymbol{V}_{i}$代替$\boldsymbol{X}_{i}$，否则保留$\boldsymbol{X}_{i}$.然后跟随蜂根据所有引领蜂的蜜源分享信息，按概率

$\begin{eqnarray} p_{i}=\frac{\text{fit}_{i}}{∑\limits^{S}_{i=1} }\text{fit}_{i} \end{eqnarray}$

(21)

进行跟随.其中fit$_{i}$为第$i$个蜜源的适应度值.也即跟随蜂产生一个[0, 1]间的随机数并与$p_{i}$进行比较，如果该随机数小于$p_{i}$则按式(20)产生一个新蜜源.搜索中，蜜源$\boldsymbol{X}_{i}$经过若干次迭代到达阈值而没有找到更好的蜜源，则$\boldsymbol{X}_{i}$将被抛弃，与之对应的引领蜂转变为侦查蜂，并按式(22)随机产生一个新的蜜源替代$\boldsymbol{X}_{i}$.

$\begin{eqnarray} &\boldsymbol{X}^{t+1}_{i}=\\ &\left\{ \begin{array}{l} L_{d}+\text{rand}(0,1)(U_{d}-L_{d}),~~尝试次数≥最大开采次数\\ \boldsymbol{X}^{t}_{i},~~尝试次数＜最大开采次数 \end{array} \right. \end{eqnarray}$

(22)

为了获取最小的均方根误差(RMSE, root mean square error)以及和ABC算法保持一致，采用式(23)作为ABC优化高斯过程回归参数的适应度函数.

$\begin{eqnarray} \text{fit}_{i}=\frac{1}{1+r_{i}} \end{eqnarray}$

(23)

其中$r_{i}$为实际值与预测值的RMSE：

$\begin{eqnarray} r_{i}=\sqrt{\frac{1}{S}∑\limits^{S}_{j=1}(y_{j}-\hat{y}_{j})^{2}} \end{eqnarray}$

(24)

ABC算法优化高斯模型回归预测模型的步骤如下：

步骤 1 生成训练数据样本集，同时参数初始化，包括蜜源的个数$S$，最大迭代次数$M$，蜜源最大开采次数为$L$，待优化的超参数$l$、$δ^{2}_{f}$、$δ^{2}_{n}$的取值范围，令$t$=1；

步骤 2 为蜜源分配一只引领蜂，按照式(20)进行搜索，产生新蜜源$\boldsymbol{V}_{i}$；

步骤 3 代入样本数据，根据高斯过程回归模型计算输出预测值，并根据式(23)计算适应度值，根据贪婪算法保留蜜源；

步骤 4 按式(21)计算蜜源被跟随的概率，跟随蜂搜索，根据贪婪算法保留蜜源；

步骤 5 判断蜜源$\boldsymbol{X}_{i}$是否应抛弃，若是，则引领蜂转变为侦查蜂，否则转步骤7；

步骤 6 侦查蜂按式(22)产生新的蜜源；

步骤 7 令$t$=$t$+1，如果满足终止条件，输出最优超参数，否则转到步骤2继续执行.

4 网络流量预测模型

设网络流量时间序列为$T$($k$)，$k$为数据采集时刻，高斯过程回归补偿ARIMA的网络流量预测模型如图 2所示.通过单步预测迭代循环的方式实现网络流量多步预测，即对$k$+1时刻的网络流量进行预测，然后将预测值看作实际值，通过滑窗机制代入网络流量序列，并剔除最旧的数据，然后按照预测算法进行$k$+2时刻网络流量的预测，直到满足最大预测步数.单步迭代预测方法较直接多步预测方法减少了模型训练时间，同时可利用拉伊达准则等算法对预测值数据进行合理性检测，结合插值操作，减少预测误差.

预测步骤描述如下：

步骤 1 预测模型训练，利用网络流量样本数据建立ARIMA预测模型，通过ARIMA的预测流量值与实际流量值误差建立高斯过程回归预测模型，模型中的超参数通过ABC算法得到；

步骤 2 利用$T$($k$), $T$($k$-1), …, $T$($k$-1-$p$)通过建立的ARIMA模型预测下一时刻的网络流量值$T$′($k$+1)；

步骤 3 利用网络流量实际值与ARIMA预测值生成预测误差序列$E$($k$), $E$($k$-1), …, $E$($k$-1-$p$)，通过高斯过程回归模型生成网络流量误差的预测值$E$′($k$+1)；

步骤 4 将$T$′($k$+1)与$E$′($k$+1)相加，得到最终的预测值$T_{p}$($k$+1)；

步骤 5 更新预测序列，将$T$′($k$+1)放入网络流量预测值队列中，将$T_{p}$($k$+1)放入网络流量实际值序列中，同时均剔除2个序列中最旧的数据，重复步骤2，直到预测结束.

5 仿真

仿真数据利用上文介绍的250组辽宁联通169核心网接入路由器的网络中期流量数据(数据采集尺度为10 min)，通过前200组数据进行ARIMA与高斯过程回归预测模型的建立，后50组数据用来进行预测模型精度的验证.

首先利用网络流量数据进行ARIMA预测模型的建立，训练完毕其模型为ARIMA(10, 1, 7)，利用该模型对50组网络流量测试样本进行了预测. 图 3为ARIMA模型的网络流量预测值与实际值的对比曲线.从图中可看出ARIMA模型的网络流量预测值在趋势上可以较好地拟合实际值，预测出了网络流量时间序列中的线性特征，但是存在预测误差过大的问题，同时该预测误差是非线性的，因此需要对其预测误差进行补偿修正.

利用得到的ARIMA预测模型对200组网络流量训练样本数据进行预测，得到200组网络流量样本预测误差，其曲线如图 4所示.从图 4可见，样本预测误差时间序列具有强烈的非线性特征，因此选择高斯过程回归模型是合适的.利用预测误差序列完成高斯过程回归预测模型的建立.

ABC算法参数设置：蜜源的个数$S$为20，最大迭代次数$M$为100，蜜源最大开采次数$L$为50.由于待优化的参数为3个，则最优解是3维向量.待优化的超参数取值范围为$l$∈[0, 10]、$δ^{2}_{f}$∈[0, 10 000]、$δ^{2}_{n}$∈[0, 100]. 图 5为ABC算法优化高斯过程回归适应度曲线，寻优之后的结果为$l$=5.61，$δ^{2}_{f}$=73.27，$δ^{2}_{n}$=13.26.利用ABC算法优化后的超参数，高斯过程回归模型对ARIMA模型网络流量预测误差的预测对比曲线如图 6所示，从图中可观察到高斯过程回归模型对于ARIMA模型的预测误差具有很好的预测能力.将ARIMA预测模型的网络流量预测值与高斯过程回归模型的网络流量误差预测值进行相加，即可得到最终的预测值，图 7为最后叠加后得到的预测值与实际值的对比曲线.

为了对比所提出预测方法的效果，图 8给出了多核LSSVM(优化后参数为$ρ$=0.89，$γ$=2.35，$σ^{2}$=1.05，$q$=4.26)^[15]、遗传算法优化回声状态网络(参数SR为0.617 3，$N$=58，IS为0.021 2，SD为0.419)^[16]、Elman神经网络(输入层数为30，最大迭代次数为3 000，迭代目标为0.000 1，隐层个数为15，输出层个数为1，利用迭代完成多步预测)^[17]、单一高斯过程回归(利用共轭梯度法确定模型的超参数为$l$=12.68，$δ^{2}_{f}$=27.23，$δ^{2}_{n}$=8.67)4种预测模型预测效果曲线对比.由图 7与图 8可观察到，在数据拟合程度上所提出的组合预测方法要优于其他几种方法.

图 9为几种预测方法的误差分布图，从图中可观察到所提出的预测方法在预测误差上要小于其他预测方法，同时误差分布更加均匀，即预测误差随着预测步长的增长变化相对平稳，受预测步长的影响较小.

为了进行预测效果的对比，引入如下4种性能指标.

1) RMSE $e_{\text{RMSE}}$

$\begin{eqnarray} e_{\text{RMSE}}=\sqrt{\frac{1}{N}∑\limits^{N}_{i=1}(T_{i}-\hat{T}_{i})^{2}} \end{eqnarray}$

(25)

其中：$N$为样本数量，$T_{i}$为网络流量实际值，$\hat{T}$$_{i}$为预测模型的网络流量预测值.

2) 平均绝对误差(MAE, mean absolute error)$e_{\text{MAE}}$

$\begin{eqnarray} e_{\text{MAE}}=\frac{1}{N}∑\limits^{N}_{i=1}|T_{i}-\hat{T}_{i}| \end{eqnarray}$

(26)

3) MAPE平均绝对百分误差(MAPE, mean absolute percentage error) $e_{\text{MAPE}}$

$\begin{eqnarray} e_{\text{MAPE}}=\frac{1}{N}∑\limits^{N}_{i=1}\frac{|T_{i}-\hat{T}_{i}|}{T_{i}}×100 \end{eqnarray}$

(27)

4) 可靠性

$\begin{eqnarray} R^{(1-a)}=\left[\frac{ξ^{(1-a)}}{N}-(1-a)\right]×100\% \end{eqnarray}$

(28)

其中$ξ^{(1-a)}$为在置信度(1-$a$)下实际值落入预测置信区间的个数.

表 2给出了几种预测方法的RMSE、MAE、MAPE的性能指标对比. 图 10是这些预测方法的RMSE、MAE以及MAPE随着预测步数增加的变化曲线. 表 2与图 10中的结果同样显示，所提出的预测方法在误差的性能指标上也优于其他预测方法.

图 11为几种预测方法的可靠性值与置信度的分布图，从图中可看出所提出的方法在可靠性上要优于其他几种方法.

6 结束语

提出了一种高斯过程回归补偿ARIMA的网络流量中期预测方法.利用ARIMA模型进行网络流量线性成分的预测，通过高斯过程回归模型进行具有非线性特征的网络流量预测误差的预测，同时采用ABC算法对高斯过程回归模型的参数进行优化.两种方法预测值累加得到最终的预测值，通过仿真实验对比表明所提出的方法具有很好的预测效果，对解决未来设计网络资源分配原则、拥塞控制策略、网络管理优化等问题提供了良好的基础.

致谢: 中国联合网络通信集团有限公司辽宁省分公司集团客户响应中心的顾斌工程师为笔者提供了网络流量数据，在此表示感谢！