WPD-RSO-ESN和SSA-RSO-ESN模型在径流时间序列预测中应用比较

杨琼波，崔东文

（1.云南省水文水资源局红河分局，云南红河661100；2.云南省文山州水务局，云南文山663000）

月径流时间序列预测研究是水文预报研究领域的重要内容之一，提高月径流时间序列预测精度对于提高城市防洪抗旱能力、保障下游生态环境用水、优化配置水资源等具有重要意义。由于受地理位置、地形条件、大气环流、人类活动等诸多因素影响，其径流时间序列呈现出较强的非线性、非平稳性和多尺度等特征，当前尚无完善的模型可准确描述其变化过程，只有对径流时间序列进行细致分析，才能有效提高径流预测精度。当前“分解-重构”的时间序列预测方法已在降水［1］、径流［2-6］预测中得到应用，并获得较好的预测效果。时间序列分解方法有小波分解（Wavelet Decomposition，WD）、经验模态分解（Empirical Mode Decomposition，EMD）、总体经验模态分解（En⁃semble Empircal Mode Decomposition，EEMD）、变分模态分解（Variational Mode Decomposition,VMD）、奇异谱分析（Singular Spectrum Analysis，SSA）等；预测模型有ARMA 模型［3］、BP 神经网络［4-5］、广义回归神经网络（Generalized Regression Neural Net⁃work，GRNN）［6］、支持向量机（Support Vector Machine，SVM）［2］、相关向量机（Relevance Vector Machine，RVM）［7］、极限学习机（Extreme Learning Machine，ELM）［8］、LSTM 神经网络［9］、秩次集对分析法［10］等。

小波包分解（Wavelet Packet Decomposition，WPD）源于小波分解小波分解，它在分解信号低频子集同时，对高频子集继续分解，具有数学释意明确、使用简洁、能够自行设定分解层数和分解使用的小波函数等优点［11，12］，已在降水［1］时间序列分解中得到应用。奇异谱分析是近年来兴起的一种研究非线性时间序列数据的有效方法，它根据所观测到的时间序列构造出轨迹矩阵，并对轨迹矩阵进行分解、重构，从而提取出代表原时间序列不同成分的子序列。研究表明，SSA 对去除噪声信号效果优于经验模态分解和小波分解，目前SSA 已在径流时间序列分解［2］中得到应用。鼠群优化（Rat Swarm Optimization，RSO）算法是G.Dhiman等人于2020年基于自然界中老鼠追逐和攻击猎物而提出的一种新型仿生优化算法，具有实现简单、参数少、寻优精度效果好等特点［13］。回声状态网络（Echo State Network，ESN）是Jaeger 等提出的一种新型递归神经网络，特点是训练过程简单，拥有动态储备池，具备短时记忆功能，在处理时间序列等非线性系统辨识问题上表现突出，一般条件下可无限逼近动态系统。在实际应用中，ESN 网络的输出权值矩阵和储备池规模、稀疏度等超参数对ESN 网络性能有着重要影响，目前主要有优化输出权值矩阵［14-16］和优化储备池规模等超参数［17-21］两种途径。

为提高月径流时间序列预测精度，本文以云南省江边街水文站1957-2014年月径流时间序列预测为例，提出小波包分解（WPD）与奇异谱分解（SSA）-鼠群优化（RSO）算法-回声状态网络（ESN）相混合的径流时间序列预测方法。首先，分别利用WPD、SSA 将实例月径流时间序列数据分解为若干子序列，有效降低径流时间序列数据的复杂性；其次，介绍RSO 算法，在不同维度条件下选取6 个典型函数对RSO 算法进行仿真测试；最后，利用RSO 算法优化ESN 网络储备池规模等超参数，建立WPD-RSO-ESN、SSA-RSO-ESN 模型，并分别构建WPD-RSOSVM、WPD-ESN、WPD-SVM 和SSA-RSO-SVM、SSA-ESN、SSA-SVM 作对比分析模型，通过实例对8 种模型预测效果进行检验及对比分析。

1 研究方法

1.1 小波包分解（WPD）

小波包分解（WPD）是在小波分解（WD）的基础上改进和发展而来的一种有效分解方法，其能同时解析信号的低频与高频部分，已在非平稳信号分解中得到广泛应用。根据分解的信号类别，小波包分解可分为离散小波转换和连续小波转换。理论上信号被小波包分解为3 层就能够提取信号中的有效信息，逼近任意非线性函数，从而解决实际实测问题［1，11，12］。连续小波转换数学描述见式（1）：

式中：x(t)表示原始信号；ψ*表示母小波函数；a表示比例因子；b表示转换参数；t表示原始信号对应时间点。

1.2 奇异谱分析（SSA）

奇异谱分解（SSA）是根据所观测到的时间序列构造出轨迹矩阵，并对轨迹矩阵进行分解、重构，从而提取出代表原时间序列不同成分子序列［22-24］。SSA包含分解与重构两个阶段：

（1）分解。选取窗口长度L构造轨迹矩阵X如下：

式中，M表示时间序列长度；L表示窗口长度，即嵌入维数，一般为不超过序列长度1/3的整数。

对矩阵X进行SVD（Singular Value Decomposition，SVD）分解：

式中：λi表示第i个特征值；Ui表示与第i个特征值对应的特征向量；Vi表示第i个主成分；d= rank(X)= max{i：λi＞ }0 。

（2）重构。根据给定阈值，取前p个分量作为原始序列的替代：记其中一个分量对应的矩阵为X*，采用对角平均的方式，将分解矩阵转换为序列，其表示为：

式中：向量数量K=M-L+ 1；L*= min(L，K)；K*= max(L，K)。

1.3 鼠群优化（RSO）算法

1.3.1 RSO算法简述

鼠群优化（RSO）算法是一种寻优精度高、全局搜索性能强的新型仿生群体智能算法，主要通过追逐猎物、攻击猎物两个过程实现待优化问题的求解［13］：

（1）追逐猎物。老鼠是群居动物，它们通过群居竞争行为来追逐猎物。为从数学上定义这种行为，假设最好老鼠搜索个体知道猎物的位置，其他老鼠体可以更新当前位置来获得最佳搜索位置。老鼠追逐猎物数学描述如下：

式中，表示当前猎物位置；表示第i只老鼠第x次迭代位置；A表示勘探参数，，R表示［1，2］范围内随机数，T表示最大迭代次数；C表示开发参数，C= 2 ·rand()，rand()表示［0，1］范围内随机数；表示当前迭代所处最佳老鼠个体位置。

（2）攻击猎物。为从数学上定义老鼠攻击猎物过程，提出以下数学表达式：

勘探和开发之间良好平衡是评估优化算法优化性能优劣的重要标志。勘探是在给定的搜索空间中探索有希望获得最优解的不同区域；而开发则是围绕有希望获得最优解的不同区域搜索最优解。RSO算法主要通过自动调整勘探参数A和开发参数C来获得勘探和开发之间的良好平衡，从而获得算法最优解。

1.3.2 RSO算法仿真验证

为测试RSO算法的寻优性能，选取Sphere等6个单、多峰目标函数在不同维度条件下对RSO 算法进行仿真测试。单峰目标函数用于评测RSO 算法的寻优精度，多峰目标函数用于评测RSO 算法的全局搜索能力。采用20 次寻优平均值对RSO 算法的寻优性能进行评估，见表1。设置RSO 算法老鼠种群规模n=50，最大迭代次数T=1000，其他参数采用RSO算法默认值。

表1 函数优化对比结果Tab.1 Function optimization comparison results

（1）对于单峰目标函数Sphere、Schwefel 2.22、Schwefel12，RSO 算法在不同维度条件下20 次寻优均获得理论最优值0，表现出较好的寻优精度。

（2）对于多峰目标函数Griewank、Rastrigin，RSO 算法在不同维度条件下20 次寻优均获得了理论最优值0；对于连续旋转不可分多峰目标函数Ackley，RSO 算法在不同维度条件下20 次寻优均获得相对理论最优值8.88e-16，表现出较好的全局搜索性能。

可见，对于上述6个测评函数，RSO 算法在不同维条件下均具有较好的寻优精度和全局搜索能力。

1.4 回声状态网络（SEN）

回声状态网络（ESN）由输入层、动态储备池（Dynamic Res⁃ervoir，DR）、输出层组成。其中，动态储备池是网络的核心部分，它是由大量的、随机和以稀疏方式连接的隐含层神经元组成，故其具有短时记忆功能［16-18］。ESN网络结构如图1所示。

图1 ESN 神经网络结构图Fig.1 ESN neural network structure diagram

设网络输入为u(n)=[u1(n)，u2(n)，…，uK(n)]T，状态为x(n)=[x1(n)，x2(n)，…，xN(n)]T，输出为y（n）=[y1(n)，y2(n)，…，yL(n)]T，其中K表示输入维数，N表示内部神经元个数，L表示输出维数。则［11-13］：

式中：Win表示输入权值矩阵；W表示内部状态权值矩阵；Wb表示输出至内部状态权值矩阵；Wo表示内部状态至输出权值矩阵；f=(f1，f2，…，fN)表示内部神经元激活函数，一般采用sigmoid 函数或tanh 函数；fo表示输出函数，一般采用线性函数。ESN 神经网络参数中，Win、W和Wb均采用随机方式生成，在学习过程中保持不变。

1.5 WPD-RSO-SEN与SSA-RSO-SEN 建模流程

RSO 算法优化ESN 超参数的基本思想是：分别将ESN 储备池规模S、谱半径R、稀疏度SD、输入单元尺度IS映射为RSO 算法老鼠位置，设计RSO-ESN模型适应度函数，将RSO-ESN模型最优问题转化为求解适应度函数全局最小时对应的最佳老鼠个体位置，即全局最优解。根据最佳老鼠个体位置与ESN 超参数的映射关系，即可得到ESN 超参数值。预测实现流程图见图2，实现步骤如下：

图2 月径流时间序列预测流程图Fig.2 Monthly runoff time series forecast flowchart

步骤一：分别利用WPD、SSA 方法将原月径流序列分解为若干独立子序列；通过自相关函数法（AFM）确定各子序列输入向量，划分训练样本和预测样本。

步骤二：初始化ESN 神经网络结构、动态储备池DR 数量、谱半径、随机生成Win、W和Wb矩阵。构建训练样本相对误差绝对值之和作为RSO-ESN模型适应度函数：

步骤三：设置老鼠种群规模n，最大迭代次数T；初始化老鼠个体位置Pi，i= 1，2，…，n。

步骤四：计算每只老鼠适应度值，在给定搜索空间中确定最佳老鼠个体位置。令当前迭代次数x= 0。

步骤五：利用式（6）更新的老鼠个体位置位置。

步骤六：计算更新后的每只老鼠适应度值。若当前最佳位置优于前代最佳位置，则保存当前最佳位置；否则保存前代最佳位置。

步骤七：令x=x+ 1。判断x是否等于T，若是，输出最佳位置，算法结束；否则重复步骤五～步骤七。

步骤九：模型评估。利用平均绝对百分比误差（Mean Abso⁃lute Percentage Error，MAPE）、平均绝对误差（Mean Absolute Error，MAE）、均方根误差RMSE（Root Mean Square Error，RMSE）对各模型进行评估，见式（10）。

2 实例应用

江边街水文站位于云南省红河州弥勒县江边乡江边村，系南盘江干流控制站，控制径流面积25 116 km2，为国家重要水文站、中央报汛站、国家水质监测站。南盘江发源于曲靖市沾益马雄山，为珠江主源，河长2 214 km，径流面积42 800 km2，云南省境内长677 km，落差1 724 m，流经14 个县市，设有沾益、西桥、高古马、小龙潭、江边街、发蒙水文站，主要支流有曲江、泸江、甸溪河、清水江、黄泥河等。开展江边街水文站月径流时间序列预测研究，对于提高城市防洪抗旱能力、保障下游生态环境用水、优化配置水资源等具有重要意义。本文数据来源于江边街水文站1957-2014年55年共660 个实测月径流时间序列数据。

2.1 月径流时序数据

2.1.1 WPD分解

1957-2014年江边街水文站逐月径流序列数据见图3。从图3 可以看出月径流数据为非平稳序列，波动较大，本文利用db4小波将1957-2014年逐月径流数据进行3层小波包分解，得到8 个分解空间的子序列数据［1］。其中，［3，0］、［3，1］、［3，2］、［3，3］低频空间数据和［3，4］、［3，5］、［3，6］、［3，7］高频空间数据波形见图4。低频空间数据振幅较大、波长较短，大致反映月径流时间序列的变化趋势；高频空间数据振幅逐渐减小、波长逐渐变长，反映月径流时间序列的波动情况。

图3 江边街站月径流变化曲线Fig.3 Monthly runoff variation curve of Jiangbianjie Station

图4 江边街站月径流WPD时间序列分解Fig.4 Decomposition of WPD time series of monthly runoff at Jiangbianjie Station

2.1.2 SSA分解

SSA方法是近年来兴起的一种研究非线性时间序列数据的有效方法，它根据所观测到的时间序列构造出轨迹矩阵，并对轨迹矩阵进行分解、重构，从而提取出代表原时间序列不同成分的子序列。由于月径流是以年为周期，因此本文设置SSA 窗口长度L= 12，即利用SSA 方法将原月径流时间序列数据分解为12 个独立的子序列IMF1～IMF12。SSA 分解与月径流随时间变化趋势如图5所示。从图5 可以看出，子序列IMF1 振幅最大、频率最高、波长最短，大致反映了原月径流时间序列的变化趋势；子序列IMF2—子序列IMF12 振幅逐渐减小、频率逐渐降低、波长逐渐变长，反映了原月径流时间序列的波动情况。

图5 江边街站月径流SSA时间序列分解Fig.5 Decomposition of SSA Time Series of Monthly Runoff at Jiangbianjie Station

2.2 时间序列建模

合理确定经分解后的各子序列的输入、输出向量是提高月径流时间序列预测精度的关键。本文采用自相关函数法（Auto⁃correlation Function，AFM）确定各子序列的输入、输出向量。确定原则：通过SPSS 软件计算各子序列的自相关系数，在滞后数H≥4 情况下，将自相关系数最大时所对应的滞后数H作为各子序列最优嵌入维数，即将预测月的前H个月径流数据作为输入向量，预测月作为输出向量，见表2、3。本文利用实例前533～536个月实测数据作为训练样本，后120个月作为预测样本。

表2 WPD分解各子序列自相关系数、嵌入维数及序列长度Tab.2 WPD decomposes the autocorrelation coefficient，embedding dimension and sequence length of each subsequence

2.3 参数设置及预测分析

2.3.1 参数设置

（1）WPD-RSO-ESN、SSA-RSO-ESN 模型：设置老鼠种群规模n=50，最大迭代次数T=100，其他参数采用算法默认值；ESN网络储备池规模S搜索范围设置为［2，20］，谱半径R、稀疏度SD、输入单元尺度IS搜索范围均设置为［0.01，0.99］，数据采用［0，1］进行归一化处理。

（2）WPD-RSO-SVM、SSA-RSO-SVM 模型：RSO 算法设置同上；SVM 惩罚因子、核函数参数搜索范围均设置为［2-10，210］，不敏感损失系数设置为0.01，交叉验证折数设置为3，数据采用［0，1］进行归一化处理。

（2）其他模型：通过反复试凑的方法确定WPD-ESN、SSAESN 模型的储备池规模、谱半径、稀疏度和输入单元尺度参数值；同法确定WPD-SVM、SSA-SVM 模型的惩罚因子和核函数参数值。

2.3.2 实例预测结果分析

利用所建立的WPD-RSO-ESN、WPD-RSO-SVM、WPDESN、WPD-SVM 模型对表2 中各子序列进行预测，将预测结果叠加即得到实例月径流最终预测结果；同样利用SSA-RSOESN、SSA-RSO-SVM、SSA-ESN、SSA-SVM模型对表3中各子序列进行预测，将预测结果叠加即得到实例月径流最终预测结果。各模型采用平均绝对百分比误差MAPE（%）、平均绝对误差MAE（m3/s）和均方根误差RMSE（m3/s）进行性能评估，结果见表4；预测效果见图6。

表3 SSA分解各子序列自相关系数、嵌入维数及序列长度Tab.3 SSA decomposes the autocorrelation coefficient，embedding dimension and sequence length of each subsequence

表4 实例120个月月径流时间序列预测精度对比Tab.4 Comparison of prediction accuracy of 120-month monthly runoff time series

依据表4及图6可以得出以下结论：

图6 8种模型月径流时间序列预测效果Fig.6 Prediction of monthly runoff time series of 8 models

（1）从小波包分解的4 种模型预测效果来看，WPD-RSOESN 模型对实例月径流时间序列预测的MAPE、MAE、RMSE分别为2.73%、2.71 m3/s、3.99 m3/s，MAPE分别较WPD-RSO-SVM、WPD-ESN、WPD-SVM 模型降低27.0%、29.6%、45.6%，MAE分别降低14.8%、21.9%、30.9%，RMSE分别降低4.10%、18.9%、20.4%；从奇异谱分解的4种模型预测效果来看，SSA-RSO-ESN模型对实例月径流时间序列预测的MAPE、MAE、RMSE分别为3.90%、4.46 m3/s、6.83 m3/s，MAPE分别较SSA-RSO-SVM、SSAESN、SSA-SVM 模型降低16.7%、20.7%、30.5%，MAE分别降低5.30%、8.40%、9.30%，RMSE分别降低2.70%、9.50%、3.90%。WPD-RSO-ESN、SSA-RSO-ESN 模型均具有较高预测精度和较好的泛化性能，将其用于月径流时间序列预测是可行的。

（2）从不同分解方法的同一模型预测效果对比来看，WPDRSO-ESN 预测精度优于SSA-RSO-ESN，WPD-RSO-SVM 优于SSA-RSO-SVM，WPD-ESN 优于SSA-ESN，WPD-SVM 优于SSA-SVM，表明WPD 分解可有效降低时间序列数据的复杂性，大大提高模型的预测精度，其对月径流时间序列的分解效果要优于SSA分解方法。

（3）从同一模型预测效果对比来看，WPD-RSO-ESN 预测精度优于RSO-ESN，WPD-RSO-SVM 优于RSO-SVM，SSARSO-ESN 优于SSA-ESN，SSA-RSO-SVM 优于SSA-SVM，表明RSO 算法能有效优化ESN 网络储备池规模等超参数及SVM 惩罚因子和核函数参数。

（4）从表4可以看出，WPD-RSO-ESN 模型的预测精度优于WPD-RSO-SVM、WPD-ESN、SSA-RSO-ESN 模型，明显优于WPD-SVM、SSA-RSO-SVM、SSA-ESN、SSA-SVM 模型，具有更佳的预测效果。

（5）从图6来看，WPD-RSO-ESN 模型预测的最大绝对百分比误差仅为17.9%，最大绝对误差仅为17.3m3/s，除少数样本外，绝大多数样本的绝对百分比误差均在0 值附近波动，具有更高预测精度和更好的预测效果。

3 结论

本文基于小波包分解（WPD）、奇异谱分解（SSA）与鼠群优化（RSO）算法、回声状态网络（ESN）建立WPD-RSO-ESN、SSARSO-ESN月径流时间序列预测模型，并分别构建WPD-RSOSVM、WPD-ESN、WPD-SVM 和SSA-RSO-SVM、SSA-ESN、SSA-SVM 作对比分析模型，利用云南省江边街水文站1957-2014年月径流时间序列数据对8 种模型进行验证，得出以下结论。

（1）鼠群优化（RSO）算法在5维、10维、30维、50维和100维条件下均具有较好的寻优精度和全局搜索能力，将RSO 算法用于ESN网络储备池规模、谱半径等超参数优化是可靠的。

（2）WPD-RSO-ESN、SSA-RSO-ESN 模型对实例月径流时间序列均具有较高预测精度和较好的泛化性能，将其用于月径流时间序列预测是可行的。模型有效地提高了径流时间序列的预测精度，为水文时间序列预测研究提供新的途径。其中，WPD-RSO-ESN 模型的预测效果要优于SSA-RSO-ESN模型。

（3）实例验证表明，小波包分解（WPD）可将复杂的径流时间序列分解为更加稳定、更具规律的空间时间序列，有效降低了时间序列数据的复杂性，大大提高了模型的预测精度，分解效果优于奇异谱分解（SSA）方法；RSO 算法能有效优化ESN 网络超参数，大大提高ESN网络的预测性能。□