基于分解–集成方法的铁路客运量预测

Railway Passenger Transportation Volume Prediction Models Based on Decomposition-Aggregation Methods

摘要: 铁路客运量预测是铁路运输组织管理工作的重要基础和主要依据之一。本文建立多种分解–集成方法对全国铁路月度客运量进行预测分析。分别利用集合经验模态分解(EEMD)、奇异谱分解(SSA)和小波分解(WT)将原始序列进行分解，再分别使用季节差分移动自回归模型(SARIMA)模型和反向传播神经网络(BP)模型及其组合模型对分解后的子序列进行拟合、预测和集成。对比研究发现采用分解–集成方法有助于提高相关模型的预测准确性，且EEMD-SARIMA-BP组合模型在所有模型中预测效果最佳。

Abstract: The forecast of railroad passenger volume is one of the important foundations and main bases of railroad transportation organization and management. In this paper, we establish multiple decom-position-aggregation methods to forecast and analyze the monthly passenger volume of national railroads. We use the ensemble empirical modal decomposition (EEMD), the singular spectrum de-composition (SSA) and the wavelet analysis (WT) respectively to decompose the original data series into several sub-series, then we process the forecast by fitting, forecasting and aggregating by the seasonal difference moving autoregressive model (SARIMA) model and back propagation neural network (BP) and their combined model of the sub-series, respectively. We find that the use of de-composition-aggregation methods helps to improve the prediction accuracy, and the combined EEMD-SARIMA-BP model has the best prediction effect among all models.

1. 引言

铁路是国家重要交通基础设施及重大民生工程，关系着千千万万户的出行和国家经济发展 [1]。目前正处于我国铁路的高速发展期间，高速铁路服务网络覆盖范围不断扩大，越来越多人享受到安全快捷的高速铁路服务，高铁出行逐渐成为人们的长短途出行首选 [2]。铁路客运量作为反映铁路旅客运输的基本产量指标，可以直接反映铁路客运市场，体现铁路网络的客运服务规模。科学的铁路旅客运输量预测方法，可以建立高速铁路的需求分析，并结合建设标准和财政补贴，完善当前铁路发展水平，最终实现运营模式的盈利和更长足的发展 [3] [4]。

目前，我国许多学者对于铁路客运量预测的研究主要有传统计量和统计方法、人工智能方法和组合预测方法 [4] - [16]。传统的计量和统计方法是基于问题的历史数据并利用模式识别，参数估计，模型检测来建立问题的数学模型。王雷 [6] 等人根据铁路月客运量的趋势性和周期性，分别采用季节性指数平滑法和季节差分自回归移动平均法建立模型。缪巧芬 [7] 等人先用SARIMA模型对铁路客运量进行建模预测，再同时采用X-13A-S季节调整方法重新建模预测，得到了X-13A-S季节调整方法较优的结论。人工智能相关的预测方法也已取得了较好的预测效果，如将神经网络理论、灰色系统理论、支持向量机理论、遗传算法理论引入到客运量预测中 [8]。侯福均和吴祈宗 [4] 采用反向传播(BP)人工神经网络算法研究了铁路客运市场的时间序列。王卓 [9] 等人对客运量预测的BP神经网络模型进行优化，得到了改进后的BP神经网络模型的预测效果较标准BP模型较优的结论。汪健雄等人 [10] 为克服BP神经网络易计算效率和泛化能力低的问题，结合Gram-Schmidt正交化定理，提出了基于双层次BP神经网络模型的铁路客运量预测模型。针对客运量时间序列的特性和影响因素，许多学者也提出了新的组合预测方式来提高预测的精度。侯丽敏和马国峰 [13] 在灰色预测理论的基础上，建立了GM(1,1)与线性回归的组合模型。贺晓霞 [14] 等人为了充分考虑客运量的周期变动性，结合灰色预测理论GM(1,1)与周期扩展模型，得到了与实际值较吻合的预测值。刘琳玥 [15] 提出了利用主成分分析法去除和减少影响原始铁路客运量因素之间的相关性，再进行PCA-BP神经网络模型分析和预测，得到较好的预测效果。

近十几年来，我国学者就铁路客运量预测展开了多维度、多角度的探讨与研究，取得了许多宝贵的成果。但是由于铁路客运量本身的季节性较强及外部影响因素较多，许多方法仍存在不足和缺陷。在此基础上，进一步探讨铁路旅客运输量预测模式，探索更多方法的结合是否能为预测带来更高的准确性，本文首先采用季节时间序列模型(SARIMA)和BP神经网络模型对铁路月客运量进行预测，再采用集合经验模态分解(EEMD)、奇异谱分解(SSA)和小波分解(WT)三种分解方法对月客运量数据进行分解，结合SARIMA模型和BP神经网络模型进行预测。最后，通过预测结果比较上述方法的有效性与精确性。

2. 季节性差分自回归移动平均模型(SARIMA)

2.1. 差分自回归移动平均模型基本原理(ARIMA)

具有下列结构的模型被称作差分自回归移动平均(Autoregressive Integrated Moving Average)模型，简称为ARIMA(p,d,q)模型 [17]：

$\{\begin{array}{l}\Phi \left(B\right){\nabla }^d{x}_t=\Theta \left(B\right){\epsilon }_t\\ E\left({\epsilon }_t\right)=0,Var\left({\epsilon }_t\right)={\sigma }_{\epsilon }^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s\ne t\\ E\left(x_s{\epsilon }_t\right)=0,\forall s<t\end{array}$ (2.1)

式中， ${\nabla }^d={\left(1-B\right)}^d$ ； $\Phi \left(B\right)=1-{\varphi }_{1}B-\cdots -{\varphi }_p{B}^p$ 为平稳可逆ARMA(p,q)模型的自回归系数多项式； $\Theta \left(B\right)=1-{\theta }_{1}B-\cdots -{\theta }_q{B}^q$ 为平稳可逆ARMA(p,q)模型的移动平均系数多项式。

2.2. 季节性差分自回归移动平均模型原理(SARIMA)

SARIMA模型来源于自差分回归移动平均模型(ARIMA)，又称为季节乘积模型。在短期相关性和季节性影响的乘积关系中，拟合模型实际上为ARMA(p,q)和ARMA(P,Q)_S相乘的结果。结合d阶趋势差分和以周期S为步长的D阶季节差分进行建模，乘法模型的构造如下列公式所示 [17]：

${\nabla }^d{\nabla }_S^D{x}_t=\frac{\Theta \left(B\right){\Theta }_S\left(B\right)}{\Phi \left(B\right){\Phi }_S\left(B\right)}{\epsilon }_t$ (2.2)

其中 $\Theta \left(B\right)=1-{\theta }_{1}B-\cdots -{\theta }_q{B}^q$， $\Phi \left(B\right)=1-{\varphi }_{1}B-\cdots -{\varphi }_p{B}^p$， ${\Theta }_S\left(B\right)=1-{\theta }_1{B}^S-\cdots -{\theta }_Q{B}^{QS}$， ${\Phi }_S\left(B\right)=1-{\varphi }_1{B}^S-\cdots -{\varphi }_P{B}^{PS}$，该乘法模型简记为 $\text{ARIMA}\left(p,d,q\right)\times {\left(P,D,Q\right)}_S$。

3. 集合经验模态分解理论(EEMD)

3.1. 经验模态分解基本原理

经验模态分解(EMD)，是美国华裔科学家黄诺登博士于1998年提出的一种新的自适应信号时频域处理技术 [18]。它根据数据本身的时间特性对信号进行分解为有限个本征模函数(IMF)，每个IMF分量分别为源信号的不同时间尺度的局部特征信号。

为了从原始序列中分解出IMF，EMD的分解过程如下：

先将数据序列分段筛选出极大值和极小值点，上下极值点的包络线 $e_{\max \left(t\right)}$ 和 $e_{\min \left(t\right)}$，用三次样条曲线拟合出来，并计算上下包络线的平均值 $m\left(t\right)$。设原始数据为 $x\left(t\right)$，在 $x\left(t\right)$ 中减去 $m\left(t\right)$ 得到：

$h\left(t\right)=x\left(t\right)-m\left(t\right)$ (3.1)

再根据预设判据判断 $h\left(t\right)$ 是否为IMF，重复以上过程直到 $h\left(t\right)$ 满足判据，则 $h\left(t\right)$ 就是需要提取的IMF： $C_k\left(t\right)$。每得到一次IMF分量后从原信号中扣除，直到剩余部分 $r_n\left(t\right)$ 为单调的序列或者常数序列。这样原始序列可分解为：

$x\left(t\right)={\displaystyle {\sum }_{i=1}^N{C}_i\left(t\right)}+r_n\left(t\right).$ (3.2)

3.2. 集合经验模态分解基本原理

集合经验模态分解(EEMD)，是为了解决EMD方法存在模态混叠等不足而提出的一种叠加高斯白噪声的多次经验模态分解。利用高斯白噪声具有频率均匀分布的统计特性，通过每次加入同等幅值的不同白噪声来改变信号的极值点特性，再对多次经验模态分解得到的相应IMF进行总体平均来抵消加入的白噪声，从而有效抑制模态混叠的产生。

EEMD分解的具体步骤 [19] [20] 如下：

第一步，设定总体平均次数(集合数) m；

第二步，将在原始序列 $x\left(t\right)$ 中加入高斯白噪声序列 ${\epsilon }_i\left(t\right),i=1,2,\cdots ,m$。

$x_i\left(t\right)=x\left(t\right)+{\epsilon }_i\left(t\right).$ (3.3)

第三步，对含噪序列 $x_i\left(t\right)$ 分别进行EMD分解，分解各自IMF分量 $C_{ij}\left(t\right),j=1,2,\cdots ,J$ 和残差分量 $r_i\left(t\right)$。

$x_i\left(t\right)={\displaystyle {\sum }_{j=1}^J{C}_{ij}\left(t\right)}+r_i\left(t\right).$ (3.4)

第四步，对于所得到对应的IMF分量求均值。从而得到原序列的第k个IMF分量 $c_k\left(t\right)$ 和剩余分量 $r\left(t\right)$ ：

$c_k\left(t\right)=\frac{1}{m}{\displaystyle {\sum }_{i=1}^m{C}_{ik}}\left(t\right),r\left(t\right)=\frac{1}{m}{\displaystyle {\sum }_{i=1}^m{r}_i}\left(t\right)$ (3.5)

第五步，原序列可分解为

$x\left(t\right)={\displaystyle {\sum }_{k=1}^m{c}_k}\left(t\right)+r\left(t\right)$ (3.6)

4. 奇异谱分析理论(SSA)

奇异谱分析(SSA) [21] [22] 是1978年Colebrook提出的一种用于非线性时序资料的新方法。建立了观察到的时间资料的轨迹矩阵，并将轨迹矩阵分解、重建，从表示原始时刻的序列中分离出各种分量信号，例如长期趋势、周期、噪声信号等，从所提取的信号中进一步分析和预测时间序列结构。

第一步，嵌入。将一维时间序列数据 $y=\left(y_1,y_2,\cdots ,y_N\right)$ 转化为其轨迹矩阵X：

$X={\left(x_{ij}\right)}_{i,j=1}^{L,K}=\left(\begin{array}{cccc}{y}_1& y_2& \cdots & y_K\\ y_2& y_3& \cdots & y_{K+1}\\ ⋮& ⋮& \ddots & ⋮\\ y_L& y_{L+1}& \cdots & y_N\end{array}\right)$ (4.1)

其中，L为选取的窗口长度， $1<L<N$, $K=N-L+1$。

第二步，奇异值分解。令 $S=X{X}^{\text{T}}$，对S进行奇异值分解后得到的其L个特征值 ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_L\ge 0$

和其所有的正交特征向量 $U_1,U_2,\cdots ,U_L$，令 $d=\max \left\{i,{\lambda }_i>0\right\}{U}_2$，记 $V_i=\frac{X^{\text{T}}{U}_i}{\sqrt{{\lambda }_i}},i=1,2,\cdots ,d$。则矩阵X的

奇异值分解可以写成：

$X=X_1+X_2+\cdots +X_d$ (4.2)

其中， $X_i=\sqrt{{\lambda }_i}{U}_i{V}_i^{\text{T}}$， $\sqrt{{\lambda }_i}$ 是矩阵X的奇异值， $U_i$ 是左特征向量， $V_i$ 是右特征向量， $\left(\sqrt{{\lambda }_i},U_i,V_i\right)$ 称为矩阵X的第i个三重特征向量。

第三步，分组。以不同的提取成分作为依据，将 $X_i$ 分为m个不同的组 $I_1,I_2,\cdots ,I_m$,并将每组的矩阵相加。如第 $I_i$ 组包含的子集为 $I_i=\left\{i_1,\cdots ,i_p\right\}$，则

$X_{I_i}=X_{i_1}+X_{i_2}+\cdots +X_{i_p}$ (4.3)

而X相应分解为

$X=X_{I_1}+X_{I_2}+\cdots +X_{I_m}$ (4.4)

其中， $X_{I_j}$ 的贡献率可表示为 ${\displaystyle {\sum }_{i\in I_j}{\lambda }_i}/{\displaystyle {\sum }_{1\le j\le m,i\in I_j}{\lambda }_i}$。

第四步，重构。设Y为 $L\times K$ 维矩阵，矩阵元为 $y_{ij}$， $1\le y_{ij}\le L$, $1\le j\le K$。定义 $L^{*}=\min \left(L,K\right)$， $K^{*}=\max \left(L,K\right)$， $N=L+K-1$， $y_{ij}^{*}=\{\begin{array}{l}{y}_{ij},若L<K\\ y_{ji},若L\ge K\end{array}$，则重构序列 $G=\left(g_0,g_1,\cdots ,g_{N-1}\right)$ 可通过下列式子计算获得：

$g_k=\{\begin{array}{l}\frac{1}{k+1}{\displaystyle {\sum }_{m=1}^{k+1}{y}_{m,k-m+2}^{*}}若0\le k\le L^{*}-1\\ \frac{1}{L^{*}}{\displaystyle {\sum }_{m=1}^{L^{*}}{y}_{m,k-m+2}^{*}}若L^{*}\le k\le K^{*}-1\\ \frac{1}{N-K^{*}}{\displaystyle {\sum }_{m=k-K^{*}+2}^{N-K^{*}+1}{y}_{m,k-m+2}^{*}}若K^{*}\le k\le N\end{array}$ (4.5)

式(4.5)本质上是对矩阵 $X_{ij}$ 在对角线方向上求出 $i+j=k+2$ 各单元的平均，并求得相应的 $g_k$ 值，从而获得 $X_{ij}$ 的重组序列G。

5. 小波变换理论(WT)

小波变换(Wavelet Transform，缩写为WT)，又称小波分析 [23] [24] [25]，是指信号通过有限长或者快速衰减的“母小波”的振荡波形来表出，它不仅传承和延续了短时傅里叶变换的局部化思想，也克服和改进短时傅里叶变换的窗口大小不随频率变化等缺点。小波变化旨在对信号进行多尺度细化分析，信号中的信息有效地通过伸缩和平移等运算的方式，以及时间和领域的局域变换被提取出来。小波分析主要包括分解、去噪和重建三个步骤。

小波分解是将某一小波基函数 $\psi$ 做位移b后，再在不同尺度a下，与尚未分析的信号y做内积。上述过程可逆，其逆过程称作小波重构。小波分解和小波重构的表达式分别如下：

$W_y\left(a,b\right)=\langle y\left(t\right),{\Psi }_{a,b}\rangle ={\left|a\right|}^{-1/2}{\displaystyle {\int }_{-\infty }^{+\infty }\overline{\Psi \left(\frac{t-b}{a}\right)}y\left(t\right)\text{d}t}$ (5.1)

$y\left(t\right)=K{\displaystyle \int {\displaystyle {\int }_{-\infty }^{+\infty }{W}_y\left(a,b\right){\left|a\right|}^{-\frac{1}{2}}\Psi \left(\frac{t-b}{a}\right)\frac{\text{d}a\text{d}b}{a^2}}}$ (5.2)

其中， $W_y\left(a,b\right)$ 为小波分解系数，a为伸缩因子，b为平移因子， ${\Psi }_{a,b}$ 为小波基函数， $\langle ,\rangle$ 为内积， $y\left(t\right)$ 为数据， $\bar{\Psi }$ 是 $\Psi$ 的共轭。

6. BP神经网络理论

反向传播神经网络(Back Propagation Neural Network，简称BP神经网络)运用输入数据的训练让网络存在联想记忆，从而使得联想记忆发挥其预测能力。模型的主要结构分为输入层，隐含层，输出层。BP神经网络的结构如图1所示。信号通过输入层、隐含层向输出层的正向传播，误差的逆向传播，以及对预测的错误部分的误差调节，使得预测结果与预期输出不断逼近 [26]。

各个层次的数学关系如下所示：

对于隐含层有：

$y_j=f\left(ne{t}_j\right),ne{t}_j={\displaystyle {\sum }_{i=1}^n{v}_{ij}{x}_i},j=1,2,\cdots ,m$ (6.1)

对于输出层有：

$o_k=f\left(ne{t}_k\right),ne{t}_k={\displaystyle {\sum }_{i=0}^m{w}_{ik}{y}_i},j=1,2,\cdots ,l$ (6.2)

一般情形下激活函数 $f\left(x\right)$ 取 $\frac{1}{1+{\text{e}}^{-x}}$，也可以根据需要取 $\frac{1-{\text{e}}^{-x}}{1+{\text{e}}^{-x}}$。

Figure 1. Structure of BP neural network

图1. BP神经网络的结构

7. 实证分析——全国铁路客运量预测分析

7.1. 数据描述

本文数据来自于国家统计局发布的月度数据报告中的铁路客运量当期值(万人)。本文选取1984年1月至2022年1月，共457条铁路月客运量数据进行建模和检验。图2中可以看出，铁路客运量当期量总体呈上升后近两年略有下降的趋势。由于新型冠状病毒疫情的爆发，2020年2月铁路客运量达到近十五年的最低值，铁路客运量仅有约3723万人。而后至2022年1月，铁路客运量的波动较大。

Figure 2. Time series of raw data

图2. 原始数据时序图

7.2. 原始序列SARIMA预测结果

本文使用R语言程序对原始序列构建SARIMA模型对铁路月客运量进行预测。一阶差分后通过平稳性和白噪声检验。首先确定短期相关模型，利用模型 $\text{ARMA}\left(\text{2},\text{4}\right)$ 作为原始序列差分后得到的短时自相关信息模型。考虑到季节自相关特性的情况下，使用 $\text{ARMA}{\left(0,2\right)}_{12}$ 模型拟合差分序列的季节自相关信息。因此，得到模型为 $\text{ARIMA}\left(\text{2},1,\text{4}\right)\times {\left(0,0,2\right)}_{12}$。对拟合模型进行检验，结果显示该模型通过残差白噪声检验，因此拟合函数可以写为：

$\nabla x_t{\nabla }^{12}=\frac{1+0.4000B+0.6108B^2-0.2056B^3}{1-0.4628B^2}\left(1-0.2849B^{12}-0.3373B^{24}\right){\epsilon }_t.$

将SARIMA模型拟合得到的结果与铁路月客运量的实际值进行对比，如图3所示。

Figure 3. Comparison of SARIMA predictions

图3. SARIMA预测结果对比图

7.3. 原始序列BP预测结果

初始化网络、训练和仿真是BP神经网络 [25] 在进行建模时的三个基本步骤。选取默认的函数线性函数(‘tansig’)，根据经验公式指定隐含层神经元个数为10。同时本文将最小均方误差(MSE)选作模型的评价指标，设置期望误差最小值为0.0001，最大训练步长为1000，学习率为0.05。其余设置则采用MATLAB神经网络工具箱中的默认设置。训练函数选择Levenberg_Marquardt算法。同时设定90%的训练样本数，5%的验证样本数和5%的测试样本数，即411个训练样本、23个验证样本和23个预测样本。

将原始序列未经任何分解后直接通过上述设定的BP神经网络模型进行训练和仿真，得到的预测结果与实际值的对比如图4所示。

Figure 4. Comparison of SARIMA predictions

图4. SARIMA预测结果对比图

7.4. EEMD分解预测模型

第一步，运用EEMD分解法对原始数据进行分解。Nstd是设置高斯白噪声的标准差，本文取0.2。NE是添加噪声的次数，本文设置为100。将加入白噪声序列的铁路月客运量的原始序列进行分解，得到7个IMF分量和1个剩余量，根据频率高到低的顺序将分解得到的分量排序，原始时间序列的波动特征成份由各自的波动代表，剩余量为无法分解的部分，如图5所示。

Figure 5. Time series of EEMD decomposed IMF components and residual components

图5. EEMD分解后的IMF分量及剩余分量时序图

第二步，对EEMD分解后的7个IMF分量和残差项分别进行SARIMA时间序列预测。SARIMA模型定阶全部由R语言程序中forecast包中的auto.arima()函数所选取的AIC最小的模型。得到7个IMF和残差项的模型参数如表1所示。

Table 1. Model order corresponding to each component

表1. 各分量对应的模型阶数

第三步，将EEMD分解后得到的7个IMF分量和残差项分别通过BP神经网络模型进行预测 [19]。隐含层神经元个数为10，网络参数设定与训练样本比例设定与前文相同，得到BP预测结果。

通过对EEMD分解后的7个IMF分量和残差项分别通过SARIMA模型和BP神经网络预测模型进行预测后，以MAPE为指标，得到预测效果如表2所示。从表2中可以看出，对于IMF1、IMF2、IMF4、IMF5、IM6、IM7，BP神经网络的预测效果较好。对于IMF3，SARIMA模型具有更高的精度。而残差项两者的预测误差都很小。

Table 2. Comparison of forecast results (MAPE)

表2. 预测结果对比(MAPE)

最后，进行EEMD-SARIMA-BP组合预测。根据表2，对IMF3和残差项选取SARIMA模型进行预测，对其他分量选取BP神经网络方法进行预测，再将得到的预测的分量值叠加后进行组合预测。得到组合预测模型的MAPE为3.96%，提高了预测精度。

取2020年3月至2022年1月共23个月的铁路月客运量的值作为预测组，将原始序列直接采用SARIMA模型预测、直接采用BP神经网络模型预测、EEMD分解后分别采用SARIMA模型预测、EEMD分解后分别采用BP模型预测和EEMD分解后SARIMA和BP模型组合预测得到结果的MAE、MAPE、RMSE进行对比，得到表3。

Table 3. Precision comparison of EEMD decomposition prediction model

表3. EEMD分解预测模型精度对比

7.5. SSA分解预测模型

第一步，采用Matlab实现奇异谱分解 [22]。首先，将数据标准化处理。再设定窗口长度L (嵌入维数)，一般设定的长度小于或者等于序列长度的一半，将原始序列延时排列成矩阵形式。本文的铁路月客运量数据为457条，将L确定为200，按非增的次序求出200个特征值。第一和第二个特征值被采用(贡献率为88.55%)，见图6。并用公式(4.5)计算其对应的重组序列，由此获得到低频段结果。剩余的特征值来计算延时矩阵，同用公式(4.5)求出对应的重构序列，从而得到高频部分。原始序列和分解成份的时序图如图7所示，可以看出铁路月客运量整体先缓慢上升后加速上升后开始逐步下降的趋势。而高频序列整体的波动情况较为稳定，随机波动性较强的部分可视为噪声成分，但由于内在蕴含着一定波动变化，不完全归于噪声，其中2020年左右的波动较其它时间更加剧烈，其代表了铁路月客运量序列中的随机性成分。

Figure 6. Low frequency sequence

图6. 低频序列

注：从上至下为原始序列、重构的低频序列、重构的高频序列，数据均标准化处理。

Figure 7. Original series and series diagram reconstructed after decomposition

图7. 原始序列与分解后重构的时序图

第二步，对分解得到的低频和高频分量分别进行SARIMA模型建模，得到高低频分量的模型参数如表4所示。取低频项和高频项分别预测的2020年3月至2022年1月共23个月的铁路客运量的值进行对比，得到平均绝对百分误差MAPE的值分别为11.11%和95.02%。

第三步，分别对低频部分和高频部分进行BP神经网络时间序列预测。模型参数设置与前文相同。图8为低频部分的预测结果，图9为高频部分的预测结果。同样选取2020年3月至2022年1月的客运量实际值和预测值进行平均相对误差的比较，得到低频分量的MAPE为0.04%，而高频分量的MAPE为86.05%。

Table 4. Model order determination parameters of high and low frequency components

表4. 高低频分量的模型定阶参数

Figure 8. BP prediction of low frequency series

图8. 低频部分的BP预测结果

Figure 9. BP prediction of high frequency series

图9. 高频部分的BP预测结果

相对于奇异谱分解的SARIMA模型预测，分解后采用BP神经网络模型预测得到的高低频分量的误差都较小。因此分解后采用BP模型进行预测可以有效地减小误差，不需要进行组合预测。

取2020年3月至2022年1月共23个月的铁路月客运量的值作为预测组，将原始序列直接采用SARIMA模型预测、直接采用BP神经网络模型预测、SSA分解后分别采用SARIMA模型预测、SSA分解后分别采用BP模型预测得到结果的MAE、MAPE、RMSE进行对比，结果如表5所示。

Table 5. Precision comparison of SSA decomposition prediction model

表5. SSA分解预测模型精度对比

7.6. 小波分解预测模型

第一步，进行序列的小波分解与重构。由于没有标准划一的方法来选定小波基函数，因此本文经过尝试后选择具有正交性和紧支撑性 [23] 的db5小波(即Daubechies小波)作为小波函数。若小波分解的层数太少，会影响到降噪后的平滑度和平稳性，而分解的层数太多，则导致了在分解时的运算错误较大 [24]，所以在综合考虑后，本文采用了三层小波分解法。图10显示了分解的结构。用Matlab软件进行分解和重组，其结果如图11所示。其中 $a_3$ 为低频序列 $c{a}_3$ 重构后的分量，基本保持了原始序列的形状，反映了铁路月客运量的长期变化趋势。 $d_1$ 、 $d_{\text{2}}$ 和 $d_{\text{3}}$ 为重构后的高频序列分量，其隐藏了原始序列的细节特征，反映了铁路客运量的波动特征和随机特征。

第二步，对小波分解及重构后得到的低频序列a3、高频序列d1、高频序列d2、高频序列d3分别进行SARIMA模型建模，得到4个分量的模型参数如表6所示。

第三步，将小波分解重构后得到的高低频数据BP神经网络训练和仿真。图12为经过小波分解重构后的序列进行分别进行BP神经网络预测叠加后的拟合图，可以看出拟合效果较好。

通过对小波分解重构后的4个序列分别通过SARIMA模型和BP模型进行预测后，同样选取2020年3月至2022年1月的客运量实际值和预测值进行平均相对误差的比较，得到预测效果如表7所示。从表7中可以看出，对于低频序列a1，BP模型和SARIMA模型预测精度差别不大，都有较小误差。而对于d1、d2，SARIMA模型具有更高的精度。对于d3，BP模型在预测准确性上有优势。

最后进行WT-SARIMA-BP组合预测。对高频序列d1和d2选取SARIMA模型进行预测，对低频序列a1和高频序列d3选取BP神经网络方法进行预测，再将得到的预测的分量值叠加后进行组合预测。

Figure 10. Structure of three-layer wavelet decomposition

图10. 三层小波分解结构图

Figure 11. Results of wavelet decomposition and reconstruction

图11. 小波分解与重构后的结果

Table 6. Model order determination parameters of high and low frequency components

表6. 高低频分量的模型定阶参数

Figure 12. WT-BP forecast result

图12. WT-BP预测结果图

Table 7. Comparison of wavelet decomposition prediction results

表7. 小波分解预测结果对比

取2020年3月至2022年1月共23个月的铁路月客运量的值作为预测组，将原始序列直接采用SARIMA模型预测、直接采用BP神经网络模型预测、WT分解后分别采用SARIMA模型预测、WT分解后分别采用BP模型和WT分解后采用组合的SARIMA和BP模型预测得到结果的MAE、MAPE、RMSE进行对比，得到表8。组合预测模型得到预测值的MAPE为9.72%，相较于WT-SARIMA模型的12.83%与WT-BP模型的9.80%有所下降。

Table 8. Precision comparison of wavelet decomposition prediction models

表8. 小波分解预测模型精度对比

7.7. 预测结果对比

7.7.1. 分解方法预测精度对比

对于BP神经网络模型，通过对比EEMD、SSA和小波分解的预测精度指标，得到表9。从表9中可以看出，三种分解方法都可以提高BP神经网络预测的精度。对于BP预测模型，综合三个指标的数据分析，精度从高到低的排序为：EEMD-BP模型、WT-BP模型、SSA-BP模型、BP模型。

对于SARIMA模型，通过对比三种分解方法的预测精度指标，得到表10。从表10可以看出，经过采用三种分解方法均可以提高SARIMA模型的预测精度。对于SARIMA预测模型，综合三个指标的数据分析，精度从高到低的排序为：WT-SARIMA模型、EEMD-SARIMA模型、SSA-SARIMA模型、SARIMA模型。

Table 9. Accuracy of BP prediction models with different decomposition methods

表9. 不同分解方法用法BP模型预测的精度指标

Table 10. Accuracy of SARIMA prediction models with different decomposition methods

表10. 不同分解方法用SARIMA模型预测的精度指标

7.7.2. 所有预测方法精度对比

对上述共10种不同的预测模型的精度指标从小到大进行大致排名，得到表11。从表11可以总结出：对于铁路月客运量数据，采用EEMD、SSA和WT三种分解方法都可以提高预测的精度。采用EEMD分解方法的预测效果普遍优于采用WT分解，而WT分解法的预测效果又优于SSA分解法。分解后采用BP模型进行预测的误差普遍小于分解后采用SARIMA模型预测得到的误差。其中，预测效果最好的是采用EEMD-SARIMA-BP组合模型进行预测，而预测效果最差的是直接使用BP模型进行预测。

Table 11. Accuracy of all prediction methods

表11. 所有预测方法的精度指标

8. 结论

铁路客运量具有较强的季节性，也具有趋势性，同时在内外部因素影响下还具有随机性和非线性等特点。因此，本文引入的具有季节效应的自回归移动平均模型和具有强非线性模拟、自学习能力的神经网络方法可对客运量进行有效的预测。基于改进预测精度的想法，本文引入集合经验模态分解、奇异谱分解和小波分解三种方法先对序列数据进行处理后再进行不同预测模型的建模。对1984年1月至2022年1月的全国铁路月客运量序列数据进行训练和预测，得到了以下结论：

1) 通过集合经验模态分解、奇异谱分解和小波分解三种分解方法处理过的数据，得到不同频率的分量再建立预测模型进行预测均可以有效地提高序列的预测精度。这说明分解数据可以为预测提取更有效的信息，使得预测模型的建模更具有针对性。

2) 在集合经验模态分解方法中，EEMD-SARIMA-BP组合模型的预测精度是本文中提到的十种预测模型中最好的。而EEMD-BP模型的精度略低于组合模型，这是因为组合模型中部分分量采用SARIMA可以得到更佳的预测效果，因此叠加后的组合模型可以提高整体预测的精确度。若不进行组合预测，将所有的IMF分量和残差项均进行BP模型预测，得到的预测结果的误差则是小于将所有分量全部进行SARIMA模型预测的误差。

3) WT-SARIMA-BP组合模型是所有小波分解预测方法中精度最高的，其次是WT-BP模型，预测效果最差的是WT-SARIMA模型。这说明，BP神经网络模型进行预测较SARIMA方法更适用于小波分解与重构后的铁路月客运量序列，可以提高预测的精确性。虽然WT-SARIMA模型的精度略低于上述两种方法，但还是高出原始序列直接采用SARIMA模型和BP模型预测的精度很多。

4) 虽然预测的精度低于上述两种分解方法的预测模型，但SSA-BP模型和SSA-SARIMA的预测精度也是高于直接用预测模型进行预测的。这表明通过奇异谱分解后的序列数据预测值较直接使用预测模型得到的预测值更加接近真实值，用该种分解方法进行分解后再预测同样可以提高预测的精度，但是对于铁路客运量预测值提升精度的效果有限。