若干新型群体智能算法优化高斯过程回归的年降水量预测

李 杰，崔东文

（1.云南省水利水电勘测设计研究院，昆明 650000；2.云南省文山州水务局，云南 文山 663000）

0 引 言

降水量预测是水资源合理开发利用的重要基础，同时是缓解区域水资源供需矛盾、提高防洪抗旱应急能力、保障区域水生态安全的重要支撑[1]。近年来，以神经网络为代表的机器学习技术因非线性问题处理能力强、无需问题解析函数、泛化能力好等优点，已在降水量预测研究中得到广泛应用，如BP 神经网络[2,3]、支持向量机（SVM）[4,5]、长短时记忆神经网络（LSTM）[6,7]等。然而，BP神经网络虽然能够以任意精度逼近函数，但存在模型结构复杂、易陷入局部最优、调节参数选取困难等不足；SVM 方法虽然需要的样本较少，但是其预测精度受限于核函数及其参数的准确选取，且易陷入局部极值；LSTM 网络虽然在时间序列预测方面具有较大优势，但LSTM 存在内部参数多、收敛速度慢、系统消耗资源大等不足。高斯过程回归（Gaussian process regression，GPR）是基于贝叶斯理论与统计学习理论的非参数监督学习方法，适于处理高维数、非线性等复杂的回归问题，已在时间序列预测分析、动态系统模型辨识、系统控制等多个领域得到广泛应用，并取得良好效果。然而，GPR 预测效果对超参数取值的依赖程度较高，目前主要采用实验试凑、经验选择等方法对GPR 超参数进行选取，存在随机性大、容易陷入局部最优等问题。针对这一问题，粒子群优化（PSO）算法[8,9]、北方苍鹰优化（NGO）算法[10]、鲸鱼优化算法（WOA）[11]、人工蚁群优化（ACO） 算法[12]等群体智能算法（swarm intelligence algorithms，SIA）偿试用于GPR 超参数优化，并取得较好的优化效果。

由于降水受地理位置、地形条件、大气环流、人类活动等诸多因素影响，其时间序列常表现出高噪声、非线性、非平衡性等特性[13]，若直接预测会使结果产生较大误差。为充分挖掘降水量时间序列所含信息，小波分解（WD）、变分模态分解（VMD）、完全集合经验模态分解（CEEMDAN）、总体经验模态分解（EEMD）、奇异谱分析（SSA）等分解技术已在降水量时间序列数据处理中得到应用，如王文川等[14]提出的基于小波分解（WD）和郊狼优化（COA）算法的长短期记忆神经网络（LSTM）降水量预测模型，徐冬梅等[15]提出VMD-时间卷积网络（TCN）月降水量组合预测模型，罗上学等[16]建立的CEEMDAN-LSTM 耦合模型，杨倩等[17]基于总体经验模态分解（EEMD） 和长短时记忆神经网络（LSTM） 构建的EEMDLSTM 耦合模型，张以晨[18]等建立的奇异谱分解（SSA）-支持向量回归机（SVR）耦合模型。然而，在实际应用中，上述分解技术存在一定不足，如WD只对低频部分进行再分解，缺乏对时序数据高频部分的分析，降低了分解效果；VMD 在非线性、非平稳数据分解中具有较高的精确性，但分解数值K 的合理选取对预测结果影响较大；CEEMDAN、EEMD 方法虽然在一定程度上解决了EMD 模型混叠和虚假分量问题，但人为经验添加噪声、分解分量过多等不足制约了其应用；SSA技术虽然具有较好的分解效果，但存在分解分量多、预测复杂度高、计算规模大等缺点。小波包变换（wavelet packet transform，WPT）衍生于WD，与之不同的是，WPT 同时将低频、高频信号再次分解，能将原始信号分解为更具规律的子序列分量，在各行业领域具有广泛的应用。

因此，为提高降水量预测精度，改进GPR 预测性能，同时拓展SIA 在GPR 超参数优化中的应用，本文基于WPT 分解技术、孔雀优化算法（peafowl optimization algorithm，POA）、沙猫优化（sand cat swarm optimization，SCSO）算法、猎豹优化（cheetah optimization，CO）算法和GPR 方法，建立WPTPOA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 模型，同时构建基于支持向量机（SVM）的WPT-POA-SVM、WPT-SCSOSVM、WPT-CO-SVM 模型、基于RBF 神经网络的WPT-POARBF、WPT-SCSO-RBF、WPT-CO-RBF 模型和未经优化的WPT-GPR 模型作对比分析模型，并利用云南省文山州1956-2021 年年降水量数据对10 种模型进行率定与验证，旨在检验WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 模型用于年降水量预测的可行性。

1 材料与方法

1.1 数据来源

本文选用文山州1956-2021年年尺度降水量数据（数据来源于文山州历年水资源公报和文山州统计年鉴），过程线如图1 所示。从图1 可以看出，文山州年降水量序列波动性较大，复杂程度较高，呈现出较强的非线性和非平稳性，不利于直接预测。文山州位于云南省东南部，辖文山、砚山等7 县1市，国土面积31 456 km2，分属珠江、红河两大流域。珠江流域面积17 309 km2，占全州总面积的54.5%，主要有南盘江、清水江、驮娘江、西洋江等，多年平均水资源量75.36 亿m3；红河流域面积14 147 km2，主要有盘龙河、八布河、南利河、迷福河等，多年平均水资源量82.34 亿m3。文山州水资源总量丰富，多年平均降水量1 208 mm，水资源总量157.7 亿m3，属相对丰水地区。区域降水高度集中，汛期（5-10 月）降水量约占全年降水量的80%，其中约60%的降水集中在7-9月，降水量年内高度集中导致水旱灾害频发和水资源供需矛盾突出。

图1 文山州1956-2021年降水量Fig.1 Precipitation in Wenshan Prefecture from 1956 to 2021

1.2 研究方法

1.2.1 小波包变换（WPT）

WPT 能同时对信号低频部分和高频部分进行分解，更适用于降水量时间序列分解。WPT 对降水量原始信号进行分解的公式为[19-21]：

重构算法为：

式中各参数意见详见文献[19-21]。

1.2.2 孔雀优化算法（POA）

POA 是Jingbo Wang 等人于2022 年受绿孔雀群求偶、觅食和追逐行为启发而提出一种新型元启发式算法。该算法通过雄孔雀、雌孔雀和幼年孔雀来模拟种群在觅食过程中的群体行为，并基于孔雀种群角色划分建立雄孔雀、雌孔雀和幼年孔雀位置更新算子来求解待优化问题[22]。POA数学描述如下：

（1）角色划分。POA 将孔雀种群分为3 个角色：雄孔雀、雌孔雀和幼年孔雀。在优化过程中，所有个体都根据适应度值进行排名，并将适应度值最优的前五只孔雀视为雄孔雀Xpci（i=1,2,…,5），剩余的前30%的孔雀视为雌孔雀XPh，其余的孔雀视为幼年孔雀XPhC。

（2）雄孔雀位置更新。雄孔雀拥有华丽的羽毛，在其找到充足的食物后四处游走，并通过摇晃羽毛来求偶，雄孔雀越漂亮，吸引的雌孔雀就越多。POA 中，孔雀拥有的适应度值越高，其围绕食物源绕圈的概率就越大，圈半径越小；适应度值差的孔雀更容易原地旋转，圆半径较大。雄孔雀位置更新描述如下：

式中：Xpci为第i只雄孔雀位置，i=1,2,…,5；Rs为旋转半径；Xr为随机向量，描述为Xr= 2 · rand(1,Dim) - 1；‖Xr‖为Xr的模数；Dim 为问题维度；r1、r2、r3、r4为均匀分布在[0,1]范围内的随机数；t为当前迭代次数。

（3）雌孔雀接近行为。雌孔雀看到雄孔雀求偶舞蹈时，往往会先靠近雄孔雀再观察四周，雌孔雀被吸引的概率与雄孔雀的适应度值成正比。雌孔雀位置描述为：

式中：XPh为雌孔雀位置；r5为[0,1]范围内均匀分布的随机数；θ为雌孔雀接近系数，描述为θ=θ0+θ1-θ0(t/tmax)，t、tmax为当前迭代次数和最大迭代次数；θ0、θ1为接近系数初始值和最大值，分别取0.1和1。

（4）幼年孔雀搜索行为。除了接近具有良好食物资源（高适应度值）的雄性孔雀外，幼年孔雀还随机搜索，希望在搜索空间中找到更高质量的食物资源。幼年孔雀位置描述为：

式中：XPcC为幼年孔雀位置；α、δ为随迭代次数动态变化的因子；Levy 为莱维飞行，一种随机游走策略；XSPc、XPcC为随机选定的雄孔雀位置和幼年孔雀位置；其他参数意义同上。

（5）孔雀间互动行为。由于雄孔雀Xpc1拥有最好的食物资源，其余四只雄孔雀将被吸引逐渐向Xpc1移动。

1.2.3 沙猫优化（SCSO）算法

SCSO 是Amir Seyyedabbasi 等人于2022 年提出一种新型优化算法。该算法灵感来自于沙猫特殊的低频噪声检测能力，该能力有助于沙猫在地面或地下快速找到并捕捉猎物。SCSO算法中，沙猫觅食行为分搜索猎物（勘探）和攻击猎物（开发）两个阶段，并通过自适应机制在勘探和开发间保持平衡[23]。SCSO数学描述简述如下：

（1）初始化。设置沙猫种群规模N，利用式（1）初始化沙猫个体位置Xi。

式中：Xi表示第i只沙猫个体初始位置；ub、lb 分别表示搜索空间上、下限值；rand(0,1)表示介于0 和1 之间均匀分布的随机数。

（2）搜索猎物（探索）。在该阶段，每只沙猫根据最佳候选位置、当前位置及其敏感度范围r更新自己位置，使其能在搜索区域中获得局部最优值。沙猫位置更新描述如下：

式中：Xi(t+ 1)为第i只沙猫第t+ 1 次迭代位置；Xbc(t)为第t次迭代沙猫最佳候选位置；Xc(t)为第t次迭代沙猫位置，即沙猫当前位置；r为沙猫敏感度范围，描述为r=rG· rand(0,1)，其中rG为灵敏度控制参数，描述为rG=SM-(2SM·t)/(T+t)，SM为沙猫听觉特性参数，取值为2；t、T分别为当前迭代次数和最大迭代次数；其他参数意义同上。

SCSO 通过定义探索与开发之间转换的平衡参数R来保持探索与开发之间的平衡，数学描述如下：

式中：R为探索与开发之间转换的平衡参数；其他参数意义同上。

（3）攻击猎物（开发）。在该阶段，SCSO 通过轮盘选择算法为每只沙猫选择一个随机角度进行位置更新以接近猎物。沙猫位置更新描述如下：

式中：Xb(t)为迄今为止沙猫最佳位置；θ为随机角度，θ介于0和360之间；其他参数意义同上。

1.2.4 猎豹优化（CO）算法

CO 算法是Mohammad AminAkbari 等人于2022 年受自然界猎豹狩猎启发而提出一种新型群体智能优化算法。该算法通过模拟猎豹在狩猎过程中搜索、坐等和攻击3种策略来实现位置更新，即待优化问题的求解[24]。CO算法数学描述如下：

（1）初始化。与其他群体智能优化算法类似，CO 算法也是从种群初始化开始。设在d维搜索空间中，猎豹初始化位置描述为：

式中：Xi,j为第i头猎豹第j维位置；UBj、LBj为第j维搜索空间上、下限值；rand 为介于0 和1 之间的随机数；n为猎豹种群规模；d为问题维度。

（2）搜索策略。猎豹在其领地（搜索空间）或周围区域进行全范围扫描或主动搜索，以找到猎物。该策略数学描述为：

（3）坐等策略。在搜索模式下，猎物可能会暴露在猎豹视野中，在这种情况下，猎豹的每一个动作都可能会导致猎物逃跑。为避免该情况发生，猎豹采取坐等伏击策略（躺在地上或躲进灌木丛）以接近猎物。该策略数学描述为：

式（12）各参数意义同上。该策略不但提高狩猎成功率（获得取优解），而且避免CO过早收敛。

（4）攻击策略。在CO 算法中，每头猎豹都可以根据逃跑猎物、领头猎豹或附近猎豹的位置来调整自己的位置，以获得最佳攻击位置。该策略数学描述为：

1.2.5 高斯过程回归（GPR）

GPR 是通过有限的高维数据来拟合出相应的高斯过程，从而来预测任意随机变量下的函数值。设训练集(X,y)={(Xi,yi)|i=1,2,…,n}，其中X=[x1,x2,…,xn]代表d维的输入向量矩阵，y=[y1,y2,…,yn]表示输出值。高斯过程（Gaussian process，GP）f(x)的有限维随机变量都服从一个多元高斯分布，其性质由均值函数m(x)和协方差函数（核函数）k(x,x＇)决定，因此高斯过程可定义为[8-12]：

考虑噪声后，GPR回归模型可表示为：

式中：ε为高斯噪声，且满足ε～N(0,σ2)。

考虑噪声观测值的先验分布为：

式中：K(X,X)为n×n阶协方差矩阵；In为单位矩阵。

观测值y和预测值f*的联合先验分布为：

测试过程中根据Bayes原理求得后验分布为：

协方差函数k(x,x＇)（又称“核函数”）体现了输入样本之间的相似关系，决定了GPR 的预测精度。GPR 核函数中超参数选取对模型精度有重要影响，本文基于二次有理式协方差函数（RQ）构建GPR，并采用上述POA、SCSO、CO 算法对GPR核函数参数、噪声方差、协方差3个超参数进行调优，以改善GPR模型预测性能。

1.3 建模流程

WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 模型的预测步骤归纳如下（其他模型参考实现）：

步骤一：利用WPT 将实例年降水量时序数据进行2 层分解，得到1 个周期项分量[2,4]和3 个波动项分量[2,1]～[2,3]，见图2。周期项分量可有效反映原始年降水量数据的周期性，波动项分量可反映原始年降水量数据的震荡变化。本文选取1956-2015 年降水量时序数据作为训练样本，2016-2021 年降水量时序数据作为预测样本。

图2 文山州年降水量时序数据分解Fig.2 Decomposition of Annual Precipitation Time Series Data in Wenshan Prefecture

步骤二：采用Cao 方法[19-21]确定图1 中周期项分量[2,4]和波动项分量[2,1]～[2,3]的输入步长a，并利用前a年降水量预测当年降水量，即输入层节点数为a，输出层节点数为1。经Cao方法计算，图1中周期项分量[2,4]和波动项分量[2,1]～[2,3]的输入步长a分别为14、13、10、11。

因此，WPT-POA-GPR 等模型的输入、输出可表述为：

式中：M为样本数量；a为输入步长。

步骤三：利用周期项、波动项分量训练样本的预测值与实际值构建均方误差（MSE）作为优化GPR 超参数的适应度函数：

步骤四：设置POA、SCSO、CO种群规模为30，最大迭代次数为100，其他采用算法默认值。初始化蜣螂、珍鲹、猎豹个体位置。

本文基于二次有理式协方差函数（RQ）构建GPR 模型，GPR 核函数、噪声方差、协方差搜索空间分别设置为[0.01,10]、[0.01,10]、[0.1,100]；SVM 选择RBF 径向基核函数，惩罚因子、核函数参数、不敏感损失系数搜索空间分别设置为[0.01,100]、[0.01,100]、[0.000 1,0.1]；RBF 扩展速度搜索空间设置为[0.001,1 000]。为验证优化效果，WPT-GPR模型超参数采用试算法确定；所有模型的输入数据均采用[0,1]进行归一化处理。

步骤五：基于式（20）计算孔雀、沙猫、猎豹个体适应度值，找到并保存最佳个体位置。令t= 1。

步骤六：分别利用上述POA、SCSO、CO 算法位置更新算子更新个体新位置。

步骤七：利用位置更新后的孔雀、沙猫、猎豹个体计算适应度值，比较并保存当前最佳个体位置。

步骤八：重复步骤六～步骤八直至满足算法最大迭代次数。

步骤九：输出最佳个体位置，该位置即为GPR 最佳超参数向量。利用该向量建立WPT-POA-GPR 等模型对周期项分量[2,4]和波动项分量[2,1]～[2,3]进行预测和加和重构。

步骤十：利用平均绝对百分比误差MAPE、平均绝对误差MAE、均方根误差RMSE和决定系数DC对模型进行评价。其中，MAPE、MAE用于反映模型预测误差大小，RMSE用于衡量观测值同真值之间的偏差，MAPE、MAE、RMSE越小，说明模型性能越优，预测精度越高；DC反映变量之间相关关系的密切程度，其值越大，说明模型性能越好。

2 结果与分析

构建WPT-POA-GPR 等10 种模型对实例年降水量进行训练及预测，结果见图3 和图4；预测相对误差、绝对误差见图5。

图3 年降水量训练误差Fig.3 Annual precipitation training error

图4 年降水量预测误差Fig.4 Annual precipitation prediction error

图5 年降水量训练-预测误差效果对比图（后6年为预测样本）Fig.5 Comparison Chart of Annual Precipitation Training and Prediction Error Effects (the next 6 years are prediction samples)

依据图3～图5可以得出以下结论：

（1） WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR模型对实例年降水量预测的MAPE在0.46%～0.52%、MAE在5.25～5.80 mm、RMSE在7.72～8.20 mm 之间，确定性系数DC≧0.996 7，预测精度优于WPT-POA-RBF、WPT-SCSO-RBF、WPT-CO-RBF、WPT-GPR 模型，远优化WPT-POA-SVM、WPT-SCSO-SVM、WPT-CO-SVM 模型；模型对实例年降水量训练精度同样优化其他7 种模型。可见，WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 模型具有更高的预测精度和更好的泛化能力，将其用于年降水量时间序列预测是可行的。

（2） 与WPT-POA-RBF、WPT-SCSO-RBF、WPT-CORBF 模型和WPT-POA-SVM、WPT-SCSO-SVM、WPT-COSVM 模型相比，WPT-POA-GPR、WPT-SCSO-GPR、WPTCO-GPR 模型预测的MAPE分别提高64.9%、76.7%以上，表明GPR 在处理高维度、小样本、非线性等复杂回归问题中表现出色，具有较好的拟合精度及泛化性能。

（3） 与WPT-GPR 模型相比，WPT-POA-GPR、WPTSCSO-GPR、WPT-CO-GPR 模型预测的MAPE 提高72.8%以上，表明POA、SCSO、CO能有效优化GPR超参数，提高GPR预测性能。

3 讨 论

提高年降水量预测精度是水文预报研究的热点和难点。然而，由于年降水量影响因素众多，往往呈现出非线性、多尺度等特征，传统单一模型难以获得较好的预测精度。当前，基于“分解算法+预测模型”的预测方法已在年降水量预测研究中得到应用，并取得较好的预测效果，但也存在一定的不足，如文献[14]建立的模型中，WD分解方法存在分解不彻底，缺乏对高频数据的分析，LSTM 存在计算量大、耗时长等不足；文献[15]中，VMD 分解方法存在分解数值K选取困难，TCN 存在运行速度慢、超参数选取依靠人工调试等不足；文献[16]中，CEEMDAN 分解方法存在计算量大、复杂度高等问题；文献[17]中，EEMD 分解方法存在人为经验添加噪声、分解分量过多等不足；文献[18]中，SSA 分解方法存在分解过程复杂，分解分量多、计算规模大等不足，SVR 存在超参数选取困难、处理大规模数据时效果不佳等不足。

本文融合了WPT 分解方法、POA/SCSO/CO 3 种新型群体智能算法和GPR 预测器，建立了WPT-POA-GPR、WPTSCSO-GPR、WPT-CO-GP 年降水量时间序列预测模型，并构建 WPT-POA-SVM、WPT-SCSO-SVM、WPT-CO-SVM、WPT-POA-RBF、WPT-SCSO-RBF、WPT-CO-RBF、WPTGPR 对比分析模型，通过文山州年降水量预测实例验证了所构建的3种模型具有较好的预测效果，模型及方法可在类似降水量预测研究中进一步研究与推广。今后可在以下方面作进一步研究：

（1） 通过实例对比验证WPT 与上述WD、VMD、CEEMDAN、EEMD、SSA分解方法的分解效果。

（2）通过与传统粒子群优化算法、遗传算法优化GPR 超参数的对比，进一步验证POA/SCSO/CO 算法在调优GPR 超参数中的优势。

（3）本文仅以年降水量数据作为模型输入，未考虑气象、地形等因素，因此，该模型及预测精度仍有进一步提升的空间。

4 结 论

为提高年降水量预测精度，本文融合了WPT 方法、POA/SCSO/CO 三种新型群体智能算法和GPR 模型，提出WPTPOA-GPR、WPT-SCSO-GPR、WPT-CO-GPR 预测模型，并构建若干对比模型，结合具体算例，主要结论有：

（1） WPT-POA-GPR、WPT-SCSO-GPR、WPT-CO-GPR模型预测误差小于其他对比模型，具有较高的预测精度和较好的泛化能力，将其用于年降水量时间序列预测是可行的。模型及方法可为相关降水量时间序列预测研究提供参考。

（2）针对GRP 在超参数寻优时依赖参数初始值、容易陷入局部最优等问题，利用POA、SCSO、CO 优化GRP 超参数，不但改善了GPR 预测性能，而且有效提升了模型的智能化水平，具有较好的适用性。

（3）针对非线性、非平稳较强的年降水量时间序列，利用WPT 对其进行平稳化处理，可得到更具规律的周期项分量和波动项分量，显著提高了年降水量的预测精度。