一种用于可靠且具有不确定性感知能力的干旱预测的智能加权自适应回归过程

《iScience》：An intelligent weighted adaptive regression process for reliable and uncertainty aware drought forecasting

　　

亮点
• 开发了一种混合机器学习模型WARP，用于精确的干旱预测
• 新的特征选择方法和多变体分解增强了预测精度
• WARP的预测表现优于现有方法，相关系数高达0.990
• 创新方法展示了可靠且低误差的干旱指数预测能力

摘要
准确的干旱预测有助于管理和减轻对水资源和农业系统的影响。因此，本研究提出了一种混合机器学习模型，称为加权自适应回归过程（WARP），基于高斯过程回归、广义岭回归、动态更新核矩阵的自适应机制，以及为输入变量分配权重的特征相关性机制。该模型使用增强型多算子Runge-Kutta算法进行优化。此外，开发了一种基于核岭回归和排列重要性的新型特征选择方法，并实施了一种增强型多变体变分模态分解来分解输入变量。WARP在伊朗三个气象站预测SPEI时，表现优于LSSVM、DRVFL和CFNN，相关系数高达0.99，最大绝对误差为0.49。不确定性分析表明其具有更低的不确定性和更窄的预测区间，证实WARP是一种可靠且有前景的干旱预测方法。

图文摘要

引言
干旱通常被定义为水文过程中长期显著缺水。这一现象是造成广泛破坏性影响的罕见自然灾害之一。由于难以准确预测其确切的开始和结束时间，干旱的影响和损害难以估算。按月度或季节预测未来干旱对于减轻其影响至关重要。干旱预测有助于供水系统更高效地运行。根据先前研究，预测干旱的最大挑战是开发能准确确定其起止时间的方法。减轻干旱影响的另一个挑战是无法预测未来几个月和几年的状况。

衡量气象干旱的两个重要参数是蒸散量和降水量。因此，干旱指数似乎通过包含参考蒸散率与降水量作为因素来考虑气候变化的后果。一些使用蒸散量和降水量变量监测干旱的指数包括标准化降水蒸散指数（SPEI）、侦察干旱指数（RDI）和帕默尔干旱指数（PDSI）。

尽管PDSI具有固有的优势，但它有几个缺点，其中之一是无法通过跟踪不同时间间隔同时产生干旱的短期和长期后果。关于RDI指数，在考虑降水量与蒸散量比率时，当蒸散量为零时RDI无法确定。其测量范围非常有限，无法充分体现温度对干旱指数值的影响。在众多干旱监测指标中，SPEI指数具有多项优势，最重要的是它能识别蒸散量变化和温度对全球变暖的影响。

干旱预测中使用的预测方法通常可分为传统机器学习和深度学习。传统机器学习模型，如高斯过程回归、自适应神经模糊推理系统、支持向量机、人工神经网络和极端梯度提升，具有较低的计算成本、在有限数据集上表现良好以及更容易模型解释等优势。然而，它们的性能高度依赖于特征质量、参数设置以及捕捉复杂非线性关系的能力。相比之下，深度学习模型能自动学习层次化特征并从数据中建模高度复杂的模式。然而，它们通常需要大量训练数据、广泛的计算资源和长时间训练，这限制了它们在许多水文研究中的适用性。

然而，近期水文时间序列预测通过使用机器学习方法引起了关注。例如，某研究者使用ANN、ANFIS、SVM和带外生变量的自回归积分滑动平均模型预测标准化降水指数参数。结果表明ARIMAX方法精度高于其他方法。Mokhtarzad等使用SVM、ANN和ANFIS方法预测博伊努尔德站的SPI。研究发现SVM方法比ANFIS和ANN方法结果更优。Kisi等使用粒子群优化算法、遗传算法、蚁群优化和蝴蝶优化算法优化ANFIS模型，以预测伊朗某些气象站的SPI参数。结果表明PSO-ANFIS预测SPI比其他方法更准确。Malik等使用PSO和黑鹰优化算法优化支持向量回归，以预测有效干旱指数。结果表明SVR-HHO算法在印度EDI参数预测中优于SVR-PSO。

在混合预测模型方面，有几项研究将机器学习模型与分解方法如变分模态分解、小波变换和经验模态分解相结合。例如，Belayneh等将小波变换方法与SVR、ANN和ARMIA结合，提出了W-ARMIA、W-ANN和W-SVR方法预测SPI参数。结果表明W-ANN模型优于其他模型。Soh等引入小波-ANFIS和小波-ARIMA-ANN方法预测马来西亚的SPEI。根据其结果，小波-ARIMA-ANN模型精确预测了SPEI3和SPEI6。Zhang等提供了一种基于小波和ANN模型（W-ANN）的混合模型，以预测中国12个月时间尺度上的SPEI。他们表明W-ANN模型表现优于ARIMA和SVM。Karbasi等将经验小波变换、离散小波变换和扩展卡尔曼滤波与多层感知机和群方法数据处理神经网络结合，预测伊朗大不里士和拉什的SPEI参数。结果表明EWT与MLP-EKT的结合能产生比其他模型更准确的预测。Rezaiy和Shabri开发了一种基于SVR和集成经验模态分解的混合方法，预测阿富汗中部的SPI参数。结果表明所提出的EEMD-SVR方法比传统SVR和ARIMA方法提供更准确的SPI预测。

Moghaddasi等提出了一种新的基于隐藏模式特征提取统计模态分解、递归特征消除和额外树回归器的混合模型，以改进1、6和12个月提前期的多步SPEI预测。该方法在伊朗两个具有湿润和超干旱气候的气象站进行了测试，并与典型基于分解的方法（包括TVF-EMD-ETR和VMD-ETR）进行了比较。实验结果表明HPFE-ETR模型在不同气候条件下提供了更好的预测精度。

Choursi等使用卫星数据、气象变量和气候指数，通过机器学习模型预测乌尔米耶湖盆地在多个时间尺度（3至12个月）的干旱。他们使用多个精度指标和可解释性工具比较了决策树、随机森林和极随机树模型。结果表明ERT模型表现最佳。分析还显示，短期干旱主要受局部因素影响，而长期干旱则更多由大尺度气候模式驱动。

根据先前的文献综述，先前研究的主要缺点可描述如下：
1. 先前研究中应用的传统方法如SVR、ANN和ARIMA不足以应对SPEI时间序列的复杂性和非线性。
2. 尽管优化算法如GA、PSO和HHO已在先前研究中应用，但它们的主要缺点是容易陷入局部最优解，尤其是在复杂的预测问题中。因此，使用稳健高效的优化方法至关重要。
3. 大多数先前研究未应用特征选择方法来识别数据集中最具影响力的特征。因此，它们可能包含无关的输入变量，这可能会负面影响模型性能。

根据先前研究的不足，本研究的主要贡献是开发一种新的有效混合机器学习模型来解决早期工作的关键弱点以预测SPEI干旱指数。在这方面，本研究开发了一种基于高斯过程回归、特征相关性机制、自适应机制和广义岭模型的新型混合模型。根据这些机制，所提出的方法称为加权自适应回归过程（WARP）模型。所提出的WARP方法使用特征相关性机制根据每个特征的重要性水平指定合适的权重。此外，实施自适应机制以根据噪声方差动态更新核矩阵。所提出的模型使用一种稳健的优化方法，即增强型多算子Runge-Kutta算法进行优化。此外，开发了一种基于核岭回归和排列重要性方法的新型特征选择方法来选择最具影响力的特征。鉴于所提出模型的上述特性，它能够有效地预测伊朗三个气象站（亚兹德、巴博尔斯尔和乌尔米耶）的SPEI，因为其先进的机制。

结果
GPR模型示意图、WARP模型示意图、三个气象站（亚兹德、巴博尔斯尔和乌尔米耶）的SPEI时间序列图以及三个气象站的位置分别如图1、2、3和4所示。

图1 GPR模型示意图
图2 WARP模型示意图
图3 SPEI-Y、SPEI-U和SPEI-B的时间序列图
图4 三个气象站（亚兹德、巴博尔斯尔和乌尔米耶）的位置

模型开发
本研究开发了一种新的混合机器学习模型，用于预测三个气象站（亚兹德、巴博尔斯尔和乌尔米耶）的SPEI。所提出的WARP模型基于高斯过程回归、广义岭回归、自适应机制和加权方法的组合。为评估所提出模型的能力，使用了两种深度学习模型（深度随机向量功能链接网络、深度极限学习机）和三种机器学习模型（LSSVM、级联前馈神经网络、最小绝对收缩和选择算子）来预测SPEI。此外，为提高预测精度，开发了一种基于KRidge模型和排列方法的新型特征选择方法。还实施了一种增强型MVMD模型来分解输入数据集。本研究的主要目的是利用SPEI的多个时间滞后来预测SPEI六个月提前量。数据集分为30%测试集和70%训练集。设计机器学习模型需要三个步骤，如下述子节所述。所提出的预测框架如图5所示。

图5 三个气象站SPEI预测框架

特征选择
特征选择是预处理阶段最重要的部分之一。它可以去除不必要的输入特征并降低预测问题的维度。本研究使用SPEI12的10个时间滞后作为输入特征，并使用所提出的特征选择方法来确定最有效的时间滞后。所提出的KRidge-PI-EMRUN特征选择方法被开发用于指定输入特征。KRidge-PI-EMRUN方法为SPEI12-B、SPEI12-U和SPEI12-Y计算的分数如图6所示。在本研究中，所提出的特征选择方法用于选择得分最高的前50%特征。从图中可知，选择了五个特征来预测三个气象站的SPEI12（t+6）。例如，为巴博尔斯尔气象站选择的特征是SPEI-B(t)、SPEI-B(t-6)、SPEI-B(t-1)、SPEI-B(t-7)和SPEI-B(t-4)。

图6 KRidge-PI特征选择结果（A）巴博尔斯尔、（B）乌尔米耶和（C）亚兹德

所选特征的分解
本研究开发了一种EMVMD方法来分解输入变量。如前一节所述，输入变量是SPEI12-B、SPEI12-U和SPEI12-Y的五个时间滞后，目标是SPEI12-B(t+6)、SPEI12-U(t+6)和SPEI12-Y(t+6)。在此方面，EMVMD被用于分解输入数据集，然后将其输入到机器学习方法中。

EMVMD通过捕获不同的频率模式来简化数据复杂性，提供高频和低频分量，并识别关键趋势。为了实施EMVMD方法，其主要参数（NIM和?）需要通过耦合MVMD和EMRUN算法来确定。表1列出了三个站点的最优NIM和?值。

参数 | NIM | ?
---|---|---
SPEI12-B | 12 | 81,310
SPEI12-U | 7 | 1,450
SPEI12-Y | 8 | 1,280

表1 所有站点的EMVMD参数

生成的本征模态函数（IMF）的数量由优化参数NIM直接决定。EMRUN算法分别确定了Babolsar、Urmia和Yazd的最优NIM值为8、7和8。因此，每个选定的输入特征被分解为相应数量的IMF。相应地，预测SPEI12-Babolsar的输入变量为8×5个IMF（等于40个输入），SPEI12-Urmia为7×5个IMF（等于35个输入），SPEI12-Yazd为8×5个IMF（等于40个输入）。

**模型调优**

机器学习模型需要精细调优以提高预测精度。本研究应用EMRUN算法为所有机器学习模型提供最优参数。种群规模（Np）和最大迭代次数（MIt）分别为50和100。需要注意的是，WARP的搜索范围为（0, 2.00E+10）；对于CFNN，层数范围为（1, 5），神经元数量范围为（1, 20）；对于DELM，神经元数量范围为（1, 6000），正则化因子范围为（0, 1）；对于LSSVM，搜索范围为（0, 2.00E+10）；对于LASSO，搜索范围为（0, 1）。参数范围的选取基于初步实验的经验。表2和表3分别列出了简单方法和基于EMVMD方法的六种机器学习模型的最优参数。

案例 | 模型 | 调优参数
---|---|---
SPEI12-B | WARP | a? = 1.98E+02, a? = 1.92E+06, a? = 5.56E+03, η? = 4.35E?03, η? = 1.21E+01
| CFNN | NuNo* =（2, 2, 1）
| DELM | NuNo =（400, 400）, AcF* = selu, ReF* = 0.11
| LSSVM | a? = 6.01E+02, a? = 9.47E+05, a? = 7.56E+05, sigma = 7.31E+05
| DRVFL | NuNo = 26, 层数 = 6, ReF = 81
| LASSO | Lambda = 0.011
SPEI12-U | WARP | a? = 1.50E+02, a? = 1.96E+04, a? = 3.37E+03, η? = 3.64E?08, η? = 1.48E+01
| CFNN | NuNo =（3, 2, 1）
| DELM | NuNo =（200, 200）, AcF = selu, ReF = 5.10E?03
| LSSVM | a? = 1.36E+08, a? = 1.96E+10, a? = 1.81E+10, sigma = 1.88E+09
| DRVFL | NuNo = 28, 层数 = 8, RFa* = 0.23
| LASSO | Lambda = 0.015
SPEI12-Y | WARP | a? = 1.66E+02, a? = 3.13E+01, a? = 3.45E+00, η? = 2.54E?01, η? = 1.47E+02
| CFNN | NuNo =（2, 2, 1）
| DELM | NuNo =（250, 250）, AcF = sigmoid, ReF = 1.10E?02
| LSSVM | a? = 2.83E+05, a? = 1.38E+12, a? = 2.00E+12, sigma = 6.90E+10
| DRVFL | NuNo = 18, 层数 = 10, RFa* = 0.31
| LASSO | Lambda = 0.041

表2 Yazd、Babolsar和Urmia的简单机器学习方法控制参数值
NuNo*：神经元数量；AcF*：激活函数；ReF*：正则化因子。

案例 | 模型 | 调优参数
---|---|---
SPEI12-B | WARP | a? = 2.00E+02, a? = 8.82E+03, a? = 3.45E+02, η? = 7.26E?03, η? = 1.10E+00
| CFNN | NuNo* =（4, 4, 1, 1）
| DELM | NuNo =（4500, 4500）, AcF = selu, ReF = 8.00E?03
| LSSVM | a? = 1.34E+05, a? = 1.78E+10, a? = 1.01E+08, sigma = 2.00E+09
| DRVFL | NuNo = 200, 层数 = 10, ReF = 102
| LASSO | Lambda = 4.00E?04
SPEI12-U | WARP | a? = 1.52E+10, a? = 6.52E+09, a? = 1.93E+10, η? = 8.88, η? = 3.62E?3
| CFNN | NuNo =（3, 3, 1）
| DELM | NuNo =（5500, 5500）, AcF = selu, ReF = 5.20E?03
| LSSVM | a? = 1.12E+10, a? = 9.36E+11, a? = 4.09E+08, sigma = 2.00E+12
| DRVFL | NuNo = 100, 层数 = 20, ReF = 498
| LASSO | Lambda = 1.20E?05
SPEI12-Y | WARP | a? = 3.86E?01, a? = 1.07E+02, a? = 2.08E+01, η? = 7.96E?05, η? = 1.13E+01
| CFNN | NuNo =（3, 3, 1）
| DELM | NuNo = 5100, AcF = selu, ReF = 1.20E?02
| LSSVM | a? = 1.00E?10, a? = 8.52E+10, a? = 6.54E+09, sigma = 1.73E+12
| DRVFL | NuNo = 120, 层数 = 20, ReF = 210
| LASSO | Lambda = 2.50E?04

表3 Yazd、Babolsar和Urmia的EMVMD-ML方法控制参数值

需要注意的是，数据集在输入模型之前需要在（0, 1）范围内进行归一化处理。这一过程应用于使所有机器学习模型具有相同的尺度。为此，可用公式36对数据集进行归一化：

**公式1：** X_No = (X_M ? X_mi) / (X_ma ? X_mi)

其中，X_No表示归一化数据集，X_M为测量值，X_mi和X_ma分别为数据集的最小值和最大值。

**基于统计指标评估机器学习模型**

本节比较了简单方法和基于EMVMD的机器学习方法在预测SPEI12-Y、SPEI12-B和SPEI12-U方面的性能。在此背景下，为所有模型计算了六项指标（R、RMSE、MAPE、IA、NSE和MaxAE）。表4、表5和表6分别报告了所有模型预测SPEI12-Y、SPEI12-B和SPEI12-U的统计结果。

在SPEI12-Y预测方面，EMVMD-WARP模型在训练集和测试集上均表现出优异的性能（见表4）。与简单方法和基于EMVMD的机器学习方法相比，它在预测SPEI12-B(t+6)方面得出了更准确的结果（训练集：R = 0.988，RMSE = 0.159，MAPE = 60.601；测试集：R = 0.987，RMSE = 0.154，MAPE = 37.552）。此外，所提出的模型还产生了更高的IA（训练集0.994，测试集0.992）和NSE（训练集0.873，测试集0.854），以及更低的MaxAE（训练集1.039，测试集0.888）。从结果来看，EMVMD-WARP表现出比其他模型更好的泛化能力，因为它在训练和测试阶段都保持了高效率。因此，EMVMD-WARP模型是预测SPEI12-B(t+6)的最佳模型，归功于其新颖的机制。这些机制使模型能够在数据中识别复杂关系，同时呈现高精度和鲁棒性。其他模型，如EMVMD-LSSVM和EMVMD-LASSO，虽然表现出竞争性性能，但误差指标不及EMVMD-WARP。简单模型（WARP、DRVFL、DELM、LSSVM、LASSO和CFNN）的表现低于基于EMVMD的模型，这突显了基于EMVMD方法的有效性。

模型 | R | RMSE | MAPE | IA | NSE | MaxAE
---|---|---|---|---|---|---
EMVMD-WARP | 训练 0.988 | 0.159 | 60.601 | 0.994 | 0.873 | 1.039
| 测试 0.987 | 0.154 | 37.552 | 0.992 | 0.854 | 0.888
WARP | 训练 0.539 | 0.858 | 211.568 | 0.647 | 0.178 | 2.713
| 测试 0.671 | 0.702 | 197.599 | 0.690 | 0.247 | 2.143
EMVMD-DRVFL | 训练 0.993 | 0.116 | 80.753 | 0.997 | 0.896 | 0.578
| 测试 0.937 | 0.326 | 40.233 | 0.962 | 0.765 | 1.730
DRVFL | 训练 0.572 | 0.853 | 182.914 | 0.594 | 0.172 | 2.560
| 测试 0.618 | 0.762 | 186.951 | 0.589 | 0.188 | 2.401
EMVMD-DELM | 训练 0.998 | 0.075 | 37.747 | 0.999 | 0.928 | 0.214
| 测试 0.818 | 0.614 | 179.017 | 0.862 | 0.315 | 2.373
DELM | 训练 0.777 | 0.651 | 222.983 | 0.839 | 0.382 | 2.460
| 测试 0.572 | 0.816 | 190.093 | 0.692 | 0.099 | 2.076
EMVMD-LSSVM | 训练 0.991 | 0.139 | 71.057 | 0.995 | 0.879 | 0.680
| 测试 0.938 | 0.323 | 32.448 | 0.963 | 0.792 | 1.690
LSSVM | 训练 0.578 | 0.839 | 213.499 | 0.636 | 0.189 | 2.661
| 测试 0.599 | 0.754 | 216.249 | 0.624 | 0.190 | 2.301
EMVMD-LASSO | 训练 0.992 | 0.129 | 62.978 | 0.996 | 0.889 | 0.674
| 测试 0.935 | 0.330 | 36.865 | 0.962 | 0.797 | 1.598
LASSO | 训练 0.599 | 0.814 | 279.790 | 0.704 | 0.219 | 2.740
| 测试 0.573 | 0.753 | 235.542 | 0.661 | 0.183 | 2.147
EMVMD-CFNN | 训练 0.990 | 0.146 | 67.607 | 0.995 | 0.873 | 0.633
| 测试 0.928 | 0.346 | 62.190 | 0.957 | 0.722 | 1.802
CFNN | 训练 0.716 | 0.709 | 268.015 | 0.819 | 0.336 | 2.393
| 测试 0.479 | 0.805 | 220.920 | 0.611 | 0.122 | 2.058

表4 Yazd站点所有机器学习模型获得的统计指标

模型 | R | RMSE | MAPE | IA | NSE | MaxAE
---|---|---|---|---|---|---
EMVMD-WARP | 训练 0.994 | 0.100 | 45.677 | 0.997 | 0.900 | 0.373
| 测试 0.990 | 0.154 | 31.992 | 0.994 | 0.872 | 0.663
WARP | 训练 0.981 | 0.236 | 60.455 | 0.980 | 0.756 | 0.800
| 测试 0.677 | 0.793 | 108.691 | 0.767 | 0.285 | 2.359
EMVMD-DRVFL | 训练 0.993 | 0.112 | 48.646 | 0.996 | 0.884 | 0.390
| 测试 0.966 | 0.292 | 60.540 | 0.977 | 0.750 | 1.286
DRVFL | 训练 0.685 | 0.677 | 188.018 | 0.788 | 0.294 | 2.229
| 测试 0.524 | 0.880 | 127.131 | 0.681 | 0.207 | 2.375
EMVMD-DELM | 训练 0.998 | 0.069 | 39.141 | 0.999 | 0.928 | 0.188
| 测试 0.817 | 0.681 | 105.333 | 0.876 | 0.383 | 1.558
DELM | 训练 0.852 | 0.495 | 144.222 | 0.903 | 0.497 | 1.727
| 测试 0.575 | 0.922 | 111.718 | 0.705 | 0.189 | 2.243
EMVMD-LSSVM | 训练 0.989 | 0.138 | 46.311 | 0.994 | 0.862 | 0.597
| 测试 0.959 | 0.298 | 59.693 | 0.977 | 0.763 | 1.430
LSSVM | 训练 0.627 | 0.724 | 192.293 | 0.740 | 0.252 | 2.433
| 测试 0.671 | 0.803 | 106.977 | 0.757 | 0.280 | 2.329
EMVMD-LASSO | 训练 0.991 | 0.128 | 47.610 | 0.995 | 0.871 | 0.609
| 测试 0.965 | 0.274 | 51.922 | 0.981 | 0.782 | 1.380
LASSO | 训练 0.621 | 0.729 | 192.810 | 0.732 | 0.245 | 2.429
| 测试 0.643 | 0.821 | 120.119 | 0.740 | 0.254 | 2.291
EMVMD-CFNN | 训练 0.940 | 0.324 | 40.899 | 0.969 | 0.812 | 3.148
| 测试 0.935 | 0.367 | 79.867 | 0.967 | 0.714 | 1.964
CFNN | 训练 0.614 | 0.739 | 181.108 | 0.722 | 0.225 | 2.142
| 测试 0.651 | 0.856 | 95.254 | 0.725 | 0.235 | 2.085

表5 Babolsar站点所有机器学习模型获得的统计指标

模型 | R | RMSE | MAPE | IA | NSE | MaxAE
---|---|---|---|---|---|---
EMVMD-WARP | 训练 0.996 | 0.083 | 30.025 | 0.998 | 0.922 | 0.384
| 测试 0.989 | 0.127 | 30.150 | 0.994 | 0.852 | 0.486
WARP | 训练 0.944 | 0.641 | 158.245 | 0.931 | 0.354 | 1.849
| 测试 0.555 | 0.735 | 159.213 | 0.711 | 0.043 | 1.651
EMVMD-DRVFL | 训练 0.998 | 0.068 | 22.434 | 0.999 | 0.934 | 0.238
| 测试 0.943 | 0.290 | 31.213 | 0.965 | 0.771 | 1.897
DRVFL | 训练 0.786 | 0.608 | 142.813 | 0.869 | 0.385 | 1.602
| 测试 0.336 | 0.909 | 203.644 | 0.584 | ?0.227 | 2.022
EMVMD-DELM | 训练 0.998 | 0.092 | 27.994 | 0.998 | 0.902 | 0.239
| 测试 0.816 | 0.647 | 118.215 | 0.797 | 0.185 | 1.555
DELM | 训练 0.781 | 0.637 | 155.791 | 0.840 | 0.339 | 1.724
| 测试 0.436 | 0.860 | 159.307 | 0.621 | ?0.078 | 2.069
EMVMD-LSSVM | 训练 0.997 | 0.071 | 23.249 | 0.999 | 0.929 | 0.235
| 测试 0.958 | 0.249 | 31.828 | 0.975 | 0.807 | 1.655
LSSVM | 训练 0.770 | 0.628 | 148.093 | 0.854 | 0.349 | 1.727
| 测试 0.471 | 0.764 | 167.622 | 0.644 | ?0.043 | 1.596
EMVMD-LASSO | 训练 0.996 | 0.089 | 35.728 | 0.998 | 0.919 | 0.443
| 测试 0.942 | 0.288 | 31.063 | 0.966 | 0.766 | 1.763
LASSO | 训练 0.757 | 0.643 | 155.475 | 0.847 | 0.351 | 1.954
| 测试 0.419 | 0.837 | 185.477 | 0.641 | ?0.122 | 1.827
EMVMD-CFNN | 训练 0.989 | 0.148 | 29.699 | 0.994 | 0.907 | 1.745
| 测试 0.936 | 0.302 | 56.311 | 0.965 | 0.682 | 1.493
CFNN | 训练 0.787 | 0.609 | 158.960 | 0.870 | 0.390 | 1.677
| 测试 0.414 | 0.853 | 196.004 | 0.641 | ?0.126 | 1.942

表6 Urmia站点所有机器学习模型获得的统计指标

在SPEI12-B的预测方面，简单方法和基于EMVMD的模型所获得的统计指标结果列于表5中。根据该表，EMVMD-WARP方法展现了最佳结果，具有最高的准确性。在测试阶段，它达到了最高的R值（0.990）、IA值（0.994）和NSE值（0.872），以及最低的RMSE（0.154）和MAPE（31.992）。这些发现表明，所提出的模型在应对数据集中复杂模式方面具有卓越能力。此外，它还产生了较低的误差率，并显示出与测量值的高度一致性。根据结果，EMVMD-LSSVM是第二好的模型，其表现良好但不如EMVMD-WARP模型出色。在测试阶段，它呈现了优于EMVMD-LASSO、EMVMD-DRVFL、EMVMD-CFNN和EMVMD-DELM的指标值（RMSE = 0.274，MAPE = 51.922，IA = 0.981，NSE = 0.782）。因此，这些结果强调了EMVMD-WARP模型作为预测SPEI12-B最佳模型的卓越能力。

对于最后一个站点的SPEI预测（SPEI-U），表6展示了机器学习模型所获得的误差指标。结果再次揭示了EMVMD-WARP在预测SPEI12方面达到了最高水平。所提出的模型能够实现更高的R（0.989）、IA（0.994）和NSE（0.852），以及更低的RMSE（0.127）、MaxAE（0.486）和MAPE（30.150），优于其他模型。EMVMD-DRVFL和EMVMD-LSSVM等其他模型也表现良好，但它们的测试指标不如EMVMD-WARP模型突出。因此，基于本节的发现，EMVMD-WARP模型脱颖而出，成为表现最优的方法，能够提供高度准确和可靠的SPEI12预测。

总之，根据测试阶段的指标，EMVMD-WARP作为表现最优的模型脱颖而出，提供了高度准确和可靠的预测。

**使用ARAS方法指定最佳模型**

本节使用加法比率评估法（ARAS）作为多准则决策方法，基于一组统计误差指标确定表现最佳的模型。通常，ARAS程序首先对指标值进行归一化处理，使所有准则具有可比的尺度，从而消除不同单位和量级带来的影响。归一化之后，为每个指标分配权重，以反映其在决策过程中的相对重要性。在本研究中，权重并非主观分配或基于专家判断，而是采用Gini系数方法。这是一种数据驱动的加权技术，其中准则权重根据每个指标在评估模型间的变异性或区分能力来确定。通过这种方式，能够更好地区分模型表现的指标被赋予更高的重要性，而信息量较少的指标则获得较低的权重。这种方法减少了主观性，提高了排序过程的客观性和鲁棒性。最后，基于加权归一化值为每个备选方案（机器学习模型）计算综合得分。

图7展示了三个站点和所有机器学习模型的ARAS结果。如图所示，EMVMD-WARP模型在三个站点中均获得了最高的综合排名，平均得分为0.987，其次是EMVMD-LSSVM（0.763）、EMVMD-LASSO（0.753）、EMVMD-DRVFL（0.732）、EMVMD-CFNN（0.626）和EMVMD-DELM（0.344）。这些排名为机器学习模型的比较评估提供了依据，并证明EMVMD-WARP方法在三个研究站点的SPEI12预测中始终优于其他模型。

图7 三个站点所有方法获得的ARAS得分

**使用散点图分析评估机器学习模型**

本节采用了一种独特的散点图分析，以更好地比较三个站点六个基于EMVMD的模型的SPEI12预测值与测量值。图8、图9和图10分别展示了SPEI12-Y、SPEI12-B和SPEI12-U的散点图。根据这些图，数据点围绕45度线（棕色线）分布。此外，上下界（蓝色和红色虚线）通过不确定性分析绘制，以评估机器学习方法的精度。为了确定边界值，需要预测的标准误差值（SEP）和t值方法。SEP衡量模型预测中的不确定性，并提供预测值偏离实际值的估计。关于计算SEP、上界和下界的主要细节见附录A。为了提高测量值-预测值散点图的可视化效果，使用高斯核密度估计（KDE）对观测值的局部密度进行了估计。标记的大小根据估计密度进行缩放，使得标记较大的区域是数据点集中的地方，反之亦然。这种技术减少了可视化中重叠点的影响，并提供了关于观测值分布的更多信息。需要注意的是，由于基于EMVMD的模型名称较长，从本节起至本文末尾将使用缩写模型名称。因此，省略"EMVMD"一词，在全文剩余部分仅使用缩写模型名称。

图8 亚兹德站点六种方法获得的散点图：(A) WARP, (B) DRVFL, (C) DELM, (D) LSSVM, (E) LASSO, (F) CFNN

图9 Babolsar站点六种方法获得的散点图：(A) WARP, (B) DRVFL, (C) DELM, (D) LSSVM, (E) LASSO, (F) CFNN

图10 Urmia站点六种方法获得的散点图：(A) WARP, (B) DRVFL, (C) DELM, (D) LSSVM, (E) LASSO, (F) CFNN

所提出的散点图通过定义上下界来展示每个机器学习模型的不确定性。根据该图，数据点越靠近45度线且落在上下范围内，模型的不确定性越小，置信度越高。本研究根据上下界之间的差异定义了预测区间（PI），以便更好地比较用于预测SPEI12的机器学习模型。PI展示了每个机器学习模型的平均不确定性范围。从图8、图9和图10可以看出，WARP模型在三个站点的PI值均最小（Yazd：PI = 0.54，Babolsar：PI = 0.54，Urmia：PI = 0.48），优于其他模型。根据结果，第二好的模型是DRVFL（Yazd：PI = 1.10，Babolsar：PI = 0.93，Urmia：PI = 0.96）。

本部分的结果表明，WARP模型在预测SPEI12方面具有最高的准确度和最低的预测不确定性。根据结果，LASSO和LSSVM方法具有中等程度的不确定性，而DELM和CFNN的不确定性较高。

利用相对误差分析评估ML模型
本部分利用相对误差（RE）分析，进一步评估ML模型在三个站点预测SPEI的性能。图11展示了WARP、LSSVM、LASSO、DRVFL和CFNN在三个站点的RE范围。该图使用带离群值的箱线图来展示RE的分布和变化。从图中可以看出，WARP模型在预测SPEI12-B（[-6.93, 6.86]）、SPEI12-Y（[-2.80, 2.88]）和SPEI12-U（[-1.26, 1.82]）时获得了较小的RE范围，这证明了所提模型相比LSSVM、LASSO、DRVFL和CFNN具有更优越的可靠性和准确度。根据结果，LSSVM、LASSO和DRVFL方法表现出中等程度的RE范围，而LASSO相对优于LSSVM和DRVFL模型。根据这些发现，WARP方法是预测三个站点SPEI12性能最佳的模型，具有最低的RE范围。

图11 各方法所得相对误差的箱线图：(A) Yazd，(B) Babolsar，(C) Urmia

利用相对误差密度分布评估ML模型
由于密度分布可以帮助识别预测值中的模式、偏差和不一致性，因此本部分利用相对误差（RE）的密度分布来评估ML模型的性能。图12展示了所有方法在三个站点的RE值的密度分布。根据图中结果，WARP模型相比其他ML方法具有更窄的置信区间（CI）（Yazd: CI = [-1.09, 1.18]，Babolsar: CI = [-1.76, 2.04]，Urmia: CI = [-0.88, 0.93]）。LASSO模型在Babolsar（CI = [-2.64, 2.96]）和Urmia（CI = [-0.91, 1.08]）获得了第二最佳的结果。此外，LSSVM方法在Yazd获得了第二最佳结果（CI = [-1.59, 1.85]）。因此，WARP模型在预测SPEI12时能够以最小的偏差获得高准确度。此外，所提方法的密度曲线峰值尖锐且集中在零附近，表明大部分误差接近于零。因此，WARP方法能够生成优越的结果来预测三个站点的SPEI12。

图12 六种方法所得相对误差的密度分布：(A) Yazd，(B) Babolsar，(C) Urmia

时间序列的不确定性分析
表7展示了四种模型（WARP、DRVFL、LSSVM和LASSO）在三个地点（Yazd、Babolsar和Urmia）的不确定性量化性能。模型的性能使用PICP、MPIW和PINAW进行衡量。一个好的模型应该具有较高的PICP值和较低的MPIW及PINAW值，这表明大多数观测值位于预测区间内，且边界更窄、更具信息量。对于Yazd，WARP获得了最高的PICP（0.982）和最低的MPIW（0.685）及PINAW（0.141）。其他模型的PICP值较低（0.955），且区间宽得多。因此，WARP在Yazd提供了更高的覆盖率和更窄的不确定性边界。对于Babolsar，所有模型达到了相同的PICP值（0.982）。然而，WARP产生了最小的MPIW（0.849）和PINAW（0.222）。DRVFL、LSSVM和LASSO生成的区间要宽得多。这表明在相同的覆盖水平下，WARP提供了最有效的不确定性估计。对于Urmia，所有模型再次达到了相同的PICP（0.982）。其中，WARP产生了最窄的区间，MPIW为0.857，PINAW为0.209。虽然LSSVM在区间宽度方面优于DRVFL和LASSO，但其MPIW和PINAW值仍大于WARP。总体而言，WARP在三个研究区域中都是最佳模型。它拥有最高或同等最高的覆盖概率，但预测区间最窄。MPIW和PINAW值更小表明WARP提供了更准确、更可靠的不确定性边界。因此，WARP可被认为是本研究中不确定性量化最有效的模型。

表7 基于PICP、MPIW和PINAW的模型比较

| 案例 | 模型 | PICP | MPIW | PINAW |
|------|------|------|------|-------|
| Yazd | WARP | 0.982 | 0.685 | 0.141 |
| Yazd | DRVFL | 0.955 | 1.666 | 0.344 |
| Yazd | LSSVM | 0.955 | 1.641 | 0.339 |
| Yazd | LASSO | 0.955 | 1.654 | 0.341 |
| Babolsar | WARP | 0.982 | 0.849 | 0.222 |
| Babolsar | DRVFL | 0.982 | 1.530 | 0.401 |
| Babolsar | LSSVM | 0.982 | 1.636 | 0.429 |
| Babolsar | LASSO | 0.982 | 1.521 | 0.399 |
| Urmia | WARP | 0.982 | 0.857 | 0.209 |
| Urmia | DRVFL | 0.982 | 1.642 | 0.399 |
| Urmia | LSSVM | 0.982 | 1.318 | 0.321 |
| Urmia | LASSO | 0.982 | 1.887 | 0.459 |

图13和图14分别展示了各ML模型在Yazd和Babolsar站预测SPEI12的不确定性边界。Urmia站的结果可在图15中观察到。图中，实际值（测量数据）和预测值分别用蓝色和橙色线表示。浅绿色阴影区域表示预测值的95%置信区间。图中显示，在WARP模型中，预测值与实际值非常接近，表明其在预测SPEI12方面具有很高的准确度。此外，所提模型中预测值周围最小的置信区间确认了相比其他方法其具有更小的不确定性和更高的可靠性。这些结果证明了利用具有低不确定性的ML模型的重要性，这对于有效的干旱条件管理和合适的决策至关重要。需要注意的是，在预测时程末端观察到的轻微偏差主要是由于目标序列的突然变化，而这些变化在训练数据中未能得到充分表示。这种效应在基准模型中更为明显，而所提模型由于增强了非线性学习能力，保持了更高的预测准确度，并表现出明显更小的边界偏差。

图13 ML模型时间序列与实际值的比较（含不确定性区间），Yazd站：(A) WARP，(B) DRVFL，(C) LSSVM，(D) LASSO

图14 ML模型时间序列与实际值的比较（含不确定性区间），Babolsar站：(A) WARP，(B) DRVFL，(C) LSSVM，(D) LASSO

图15 ML模型时间序列与实际值的比较（含不确定性区间），Urmia站：(A) WARP，(B) DRVFL，(C) LSSVM，(D) LASSO

Taylor图分析
本部分采用Taylor图作为有效的图形分析工具，评估ML方法预测SPEI12的能力。为比较ML方法，使用三个指标（SD、R和中心化均方根误差[cRMSE]）。根据图中结果，每种方法由一个特定点表示，越接近参考点的点表示性能越好。如图16所示，WARP模型在三个站点都是最接近参考点的模型。它能够以最高的R值和与实际数据集（参考点）相似的SD来预测SPEI。因此，WARP模型是预测SPEI12具有卓越精度和可靠性的最佳模型。此外，DRVFL在Yazd和Babolsar具有第二佳的效率，而LSSVM可作为Urmia站预测SPEI12的第二佳方法。因此，WARP模型被确定为预测干旱指数（SPEI）最有效的方法，使其成为环境监测和管理的重要工具。

图16 各方法所得Taylor图：(A) Yazd，(B) Babolsar，(C) Urmia

WARP方法与先进方法的比较
本部分将WARP模型的性能与两种先进的新ML模型进行比较：深度混合核极限学习机（DHKELM）和基于L1范数鲁棒正则化的非对称C损失极限学习机（L1-ACELM）。如前所述，基于EMVMD的模型优于简单模型。为此，本部分将EMVMD-WARP模型与两种先进方法的性能进行比较。

表8展示了EMVMD-WARP与EMVMD-DHKELM和EMVMD-L1-ACELM在Babolsar、Urmia和Yazd三个站点的比较结果。EMVMD-WARP在训练和测试阶段均表现最佳。它在几乎所有情况下都拥有最高的R和IA值，并且在所有情况下都具有最低的误差（RMSE、MAPE和MaxAE）。在Babolsar，EMVMD-WARP模型达到了非常高的准确度，测试R = 0.990，RMSE = 0.154，而其他模型表现较差。Urmia的情况也是如此，其测试R值等于0.989，RMSE为0.127。对于Yazd，EMVMD-WARP再次是最优模型，与其他模型相比具有最低的误差率和最高的准确度。总体而言，EMVMD-WARP是最佳模型，因为它为所有数据集提供了更稳定的结果、更高的准确度和更低的预测误差。

表8 EMVMD-WARP模型与先进ML模型的比较

| 模型 | R | RMSE | MAPE | IA | NSE | MaxAE |
|------|------|------|------|------|------|-------|
| Babolsar | EMVMD-WARP 训练 | 0.994 | 0.100 | 45.677 | 0.997 | 0.900 | 0.373 |
| Babolsar | EMVMD-WARP 测试 | 0.990 | 0.154 | 31.992 | 0.994 | 0.872 | 0.663 |
| Babolsar | EMVMD-DHKELM 训练 | 0.980 | 0.259 | 78.847 | 0.981 | 0.740 | 0.793 |
| Babolsar | EMVMD-DHKELM 测试 | 0.967 | 0.284 | 54.619 | 0.978 | 0.746 | 0.835 |
| Babolsar | EMVMD-L1-ACELM 训练 | 0.933 | 0.351 | 95.883 | 0.962 | 0.644 | 0.926 |
| Babolsar | EMVMD-L1-ACELM 测试 | 0.888 | 0.452 | 88.616 | 0.938 | 0.560 | 1.198 |
| Urmia | EMVMD-WARP 训练 | 0.996 | 0.083 | 30.025 | 0.998 | 0.922 | 0.384 |
| Urmia | EMVMD-WARP 测试 | 0.989 | 0.127 | 30.150 | 0.994 | 0.852 | 0.486 |
| Urmia | EMVMD-DHKELM 训练 | 0.986 | 0.207 | 76.311 | 0.988 | 0.790 | 1.012 |
| Urmia | EMVMD-DHKELM 测试 | 0.968 | 0.266 | 62.339 | 0.980 | 0.758 | 1.175 |
| Urmia | EMVMD-L1-ACELM 训练 | 0.963 | 0.262 | 73.265 | 0.980 | 0.731 | 0.988 |
| Urmia | EMVMD-L1-ACELM 测试 | 0.932 | 0.353 | 84.449 | 0.962 | 0.657 | 0.974 |
| Yazd | EMVMD-WARP 训练 | 0.988 | 0.159 | 60.601 | 0.994 | 0.873 | 1.039 |
| Yazd | EMVMD-WARP 测试 | 0.987 | 0.154 | 37.552 | 0.992 | 0.854 | 0.888 |
| Yazd | EMVMD-DHKELM 训练 | 0.966 | 0.285 | 144.568 | 0.976 | 0.732 | 1.884 |
| Yazd | EMVMD-DHKELM 测试 | 0.938 | 0.417 | 154.377 | 0.959 | 0.677 | 2.208 |
| Yazd | EMVMD-L1-ACELM 训练 | 0.933 | 0.348 | 126.350 | 0.962 | 0.658 | 2.118 |
| Yazd | EMVMD-L1-ACELM 测试 | 0.898 | 0.490 | 179.732 | 0.934 | 0.569 | 2.243 |

图17展示了三个ML模型在三个研究区域（Babolsar、Urmia和Yazd）的ARAS结果。最终得分绘制在y轴上，得分越高越好。从结果可以看出，EMVMD-WARP在所有三种情况下得分最高。其在Babolsar、Urmia和Yazd的得分均为1，表明它是所有地点的最佳模型。EMVMD-DHKELM是第二模型，其得分低于WARP：Babolsar为0.66，Urmia为0.50，Yazd为0.39。表现最差的是EMVMD-L1-ACELM，在Babolsar和Urmia的ARAS得分约为0.44，在Yazd的ARAS得分约为0.35。EMVMD-WARP在所有区域的综合ARAS得分最高，是最佳模型。

图17 三种情况下ML模型的ARAS得分

结论
干旱是最危险的自然灾害之一，因为它对环境造成严重破坏并带来巨大的经济损失。预测干旱是一项困难的工作，也是水资源管理者面临的最大挑战之一。因此，本研究开发了一种独特的混合ML模型，该模型基于GPR、GRM、一种动态更新GPR核矩阵的自适应机制以及一种为特征分配权重集的 特征相关性机制。所提模型称为加权自适应回归过程（WARP），通过一种高效的优化算法（EMRUN算法）进行优化。此外，本研究还开发了一种基于核岭回归和排列重要性（PI）技术的新特征选择方法，以识别有效的特征来预测SPEI参数。另外，利用EMRUN优化的MVMD方法增强版（EMVMD）被应用于分解输入数据集。最后，所提框架用于预测伊朗三个站点（Yazd、Babolsar和Urmia）的SPEI12。研究结果总结如下：

1. 所提WARP方法在三个站点（Yazd、Babolsar和Urmia）预测SPEI12时，相较于统计准则（R、MAPE、RMSE、NSE、IA和MaxAE）提供了更准确的结果。
2. WARP模型获得了最显著的R值，训练阶段变化范围为0.988至0.996，测试阶段为0.987至0.990。
3. 基于最低MaxAE，训练阶段变化范围为0.373至1.039，测试阶段为0.486至0.888，表明所提模型能够比LSSVM、DRVFL、DELM、CFNN和LASSO方法更好地降低极端误差。
4. 各种图形分析，如散点图、相对误差图、密度图和Taylor图，证明所提WARP模型比其他模型能更精确地预测SPEI12。
5. 关于不确定性分析，WARP模型在三个站点预测SPEI12时呈现了比其他方法更小的不确定性区间。所提模型的不确定性范围从-0.601到0.078。基于这一窄区间，它可以比其他方法更有效地降低不确定性。

因此，本研究证明所提框架可作为预测SPEI干旱指数的一种有用且高效的工具。该框架提供了准确可靠的预测，使其在干旱监测和管理方面具有很高的价值。该工具使决策者能够更好地规划和应对干旱情况。未来的研究可以将所提WARP框架应用于其他具有不同气候条件的地区，以评估其鲁棒性。此外，将模型与气候变化情景相结合可能有助于预测未来气候变化对干旱模式的影响。

研究局限性
本研究的一个局限性在于特征选择在信号分解之前执行，而未考察先分解后进行特征选择的替代工作流程。

此外，特征选择与EMVMD参数优化是依次进行的，而非同时进行。虽然这种方法降低了计算复杂度，但探索联合优化框架以识别更高效的特征子集和分解参数将很有意思。未来方向：未来工作应探索并比较其他特征提取技术，特别是HPFE技术[62]，以评估它们相对于所采用的MVMD框架的有效性和对预测性能的潜在影响。资源可用性：主要联系人：所有有关进一步信息和研究资源的请求均应联系主要联系人Iman Ahamdianfar（[e-mail]）。材料可用性：本研究未生成新的独特试剂。数据和代码可用性：• 本文中报告的数据将由主要联系人根据请求共享。• 本研究使用基于MATLAB和Python的实现进行机器学习建模。计算程序可根据合理请求向主要联系人获取。• 重现或重新分析本研究中所呈现数据所需的任何其他信息，均可根据合理请求向主要联系人获取。致谢：作者真诚感谢编辑和匿名审稿人提出的宝贵意见、建设性建议和富有洞察力的反馈，这些意见显著提升了本手稿的质量和清晰度。感谢加拿大自然科学和工程研究委员会-可持续农业研究计划资助金对本合作研究的支持。作者贡献：概念化，H.J.和I.A.；方法论，H.J.、I.A.和M.K.；软件，H.J.和I.A.；验证，H.J.、I.A.、M.K.和A.A.F.；形式分析，H.J.和I.A.；调查研究，H.J.；资源，H.J.和I.A.；撰写——初稿，H.J.、I.A.和M.K.；撰写——审阅与编辑，H.J.、I.A.、M.K.和A.A.F.；监督，H.J.、I.A.和A.A.F。所有作者均已阅读并同意所发表的手稿版本。利益声明：作者声明不存在竞争利益。生成式AI和AI辅助技术在手稿准备过程中的声明：作者在写作过程中使用生成式人工智能（AI）和AI辅助技术进行语法纠正。STAR★方法：关键资源表试剂或资源 | 来源 | 标识符：气象数据 | 本研究 | 报告于表9和表10；软件与算法：Python（v3.8或更高版本）| Python软件基金会 | https://www.python.org/；MATLAB R2018a | MathWorks | https://www.mathworks.com/；龙格-库塔算法 | 本文 | https://github.com/ImanAhmadianfar-AIOpen方法细节：加权自适应回归过程（WARP）：水文和水旱数据集通常包含非线性关系、变化的噪声水平和多重共线性。这些特征降低了传统预测模型的准确性。许多现有方法可以单独学习线性或非线性模式。高斯过程回归适用于非线性建模，但通常基于固定核函数和恒定噪声假设。它还假设所有输入变量同等重要，而不考虑它们作为预测变量的不同贡献。广义回归和岭回归模型能够处理线性结构，但可能无法描述非线性动态。此外，一些先前的混合模型使用简单平均方案组合预测[30]，这可能无法明确考虑单个模型之间的分歧或预测不确定性。为了解决这些问题，本文引入了加权自适应回归过程（WARP）。WARP采用特征相关性加权、自适应核更新和分歧感知融合来提升鲁棒性和预测准确性。高斯过程回归：高斯过程回归（GPR）是一种强大的非参数概率机器学习方法，广泛用于回归问题，特别是在水文和环境建模中，因为它能够量化预测中的不确定性。在GPR中，假设观测目标变量由一个受噪声干扰的潜在底层函数生成。令xi表示输入向量，yi表示大小为M的数据集中对应的观测输出。输入与输出之间的关系可表示为：??=???(??)+??（公式2）其中g表示指高斯过程的基本函数，即g～GaP(n,l)。l和n分别表示协方差和均值。高斯过程的假设是，有限数量的函数值呈多元高斯分布，从而为非线性关系提供了灵活的模型，而无需定义固定的函数形式。此外，假设噪声?呈正态分布，均值为零，方差为??2??，??～??????????????????(0,??2??)。这一假设便于从观测不确定性中区分真实的底层信号，特别是对于受测量误差和环境变化影响的水文数据。因此，y的分布可表示为公式3：??～??????????????????(0,??????(????,????)+??2???????)（公式3）其中Im表示单位矩阵，CM(xi,xj)表示协方差（核）函数，用于衡量输入向量xi和xj之间的相似性。这里，????=[?????1,?????2,…,???????]?和????=[?????1,?????2,…,???????]?分别是第i个和第j个观测对应的d维输入特征向量。对所有训练样本对评估CM(xi,xj)形成M×M协方差矩阵。这种协方差结构控制信息如何在观测之间共享，并决定回归模型的平滑度和灵活性。协方差函数定义为：??????(????,????)=cos?(2???·??1·????2)·exp?(???02·??3)（公式4）其中??0=∥?????????∥2，??=?????????（公式4.1）其中a1、a2和a3是协方差函数的超参数，a1控制周期分量的频率，a2控制余弦项的周期长度（缩放因子），a3通过指数分量控制平滑度和相关性衰减率。这些参数共同调节协方差结构的周期性、尺度和平滑度，使核函数能够有效捕获水文时间序列中的复杂时间模式。实际上，CM?(xi,xj)称为核函数，可表示为KrF?(xi,xj)。然后，测试数据集((X?,Y?))用于估计M?维向量Y。这里，X表示输入（训练）矩阵，定义为??=[??1,??2,…,????]?，其中??∈???×??。X的每一行对应一个d维输入向量????=[?????1,?????2,…,???????]?，其中i=1,2,…,M。M表示训练数据集中的样本数，而M?表示测试数据集中的样本数。相应地，KrF(X,X)是M×M协方差矩阵，KrF(X?,X?)是M?×M?协方差矩阵。y和Y的联合先验分布可通过公式5获得。然后公式变为[????]～???([???(??)???(???)],[?????????(??,??)+??2????????????????(??,???)?????????(???,??)?????????(???,???)])，其中?????????(??,??)=[?????????(????,????)]????,??=1，?????????(??,???)=[?????????(????,?????)]??=1,…,?????=1,…,??。（公式5）?????????(??,??)=????????????(??1,??1)??????????(??1,????)????????????(????,??1)??????????(????,????)???（公式5.1）?????????(??,???)=????????????(??1,???1)??????????(??1,??????)????????????(????,???1)??????????(????,??????)???（公式5.2）其中KrF(X?,X?)表示测试数据集的协方差矩阵。这种公式使模型能够捕获训练数据内部以及训练数据与未见测试样本之间的相关性。然后，使用贝叶斯回归计算Y的后验分布，定义为：???(??∣??,???,??)=???(???,???2)（公式6）???=???(???)+?????????(???,??)?[?????????(??,??)+??2?????]?1?(??????(??))假设零均值高斯过程，即???(??)=0，???(???)=0（公式6.1）公式6.1简化为???=?????????(???,??)?[?????????(??,??)+??2?????]?1???，预测协方差为???2=?????????(???,???)??????????(???,??)?[?????????(??,??)+??2?????]?1??????????(??,???)（公式6.2）其中ρ?表示预测值向量，也称为???。预测均值提供最终估计，而方差表示相关不确定性，使GPR特别适合不确定条件下的水文预测。因此，基于GPR的预测值可通过公式7获得。图1显示了GPR模型的示意图：???=?????????(???,??)?[?????????(??,??)+??2?????]?1???（公式7）其中[?????????(??,??)+??2?????]?1表示正则化训练协方差（核）矩阵的逆矩阵。正则化项??2???????通过考虑观测噪声方差来改善协方差矩阵的数值稳定性和可逆性。特征相关性机制：所提出的模型采用特征相关性机制（FRM）进行特征选择和训练过程中的权重确定。这是因为水文和环境数据的性质不同，同等使用所有输入变量会降低模型的有效性并导致变量冗余。上一节中的模型（GPR）可以通过核函数处理非线性关系。然而，它假设所有输入变量在构建协方差结构时同等重要。因此，提出FRM与GPR结合，引入自适应特征级加权策略，以在核学习之前改善输入表示。FRM基于正则化最小二乘法与核函数的结合，可以捕获输入变量与目标输出之间的线性和非线性关系。这种公式扩展了经典GPR，将重点从仅对样本间关系建模转移到额外评估单个特征的影响。该技术使用与GPR相同的相似性结构，通过注入核矩阵KrF，但在特征相关性层面，确保与第2.1节中建立的概率结构兼容。特征权重向量为：??=(????×????????×??+??1×??)?1?(????×(????????×y))（公式8）其中ω表示特征重要性向量。η1是FRM特有的正正则化参数。η1在[0, 100]范围内优化。用于确保(XTKrFX)求逆时的数值稳定性并防止特征权重向量ω估计中的过拟合。具体来说，η1控制应用于特征重要性解的收缩。这在处理多重共线和噪声较大的水文输入变量时很重要。调节η1使模型能够在对特征变化的敏感性和对噪声的鲁棒性之间找到平衡。所提出的模型使用权重向量ω对输入变量进行缩放。因此，具有缩放特征的输入矩阵(Xscaled)可表示为：???????????????????=??×??????????????????(????,??,1)（公式9）其中ωT表示向量ω的转置，repmat用于在整个数据集中复制ωT。实际上，这种缩放方法用于为机器学习模型创建更有效的输入特征。该变换在GPR核评估之前修改输入空间，确保信息量更大的变量对相似性测量具有更强的影响，而不相关的特征被抑制。广义回归模型：广义岭模型（GRM）作为所提出的混合框架的基本组件引入，以解决GPR的局限性。GPR能够高效学习非线性关系，但不通过基于核的学习显式建模数据中的线性结构。因此，建议GRM提供稳定的线性和半参数建模组件，通过捕获仅通过核函数无法有效建模的结构相互作用来提高所提出方法的总体灵活性。所提出的GRM利用广义线性模型（GLM）与岭回归正则化相结合，提供处理不同数据分布的能力，同时保证数值稳定性。GLM与岭回归正则化的结合在水文建模背景下特别有用，因为目标变量通常由多样化的统计行为特征化，例如正态、泊松或伽马分布。GRM确保分布假设不会限制所提出的框架，同时允许其通过灵活的建模结构适应不同的数据集特征。在本研究中，GRM在灵活的GLM框架内实现，可以通过适当的链接函数适应多个分布族（例如高斯、二项和泊松）。然而，对于本工作考虑的水文预测任务，模型使用具有恒等链接函数的高斯（正态）分布配置，定义为：??(??)=??，?′?(??)=1（公式10）线性预测子定义为：??=?????（公式11）其中X表示输入变量，β是系数向量，?表示线性预测子。系数向量使用加权最小二乘法（WLS）估计，公式为：??=(??????????+??2???)?1?(??????????′)（公式12）其中??=diag?(?′?(??)??????(??))（公式12.1）和??′=??+(?????)?′?(??)（公式12.2）这里，h′(?)表示链接函数关于?的导数，Vr(Y)是响应变量的方差函数。在本研究中使用的正态（高斯）假设下，方差函数简化为：??????(??)=??2（公式12.3）其中σ2表示误差项的方差，定义为σ2=Var(ε)，其中ε表示观测值与预测值之间的随机误差。在此背景下，σ2被假设为常数且不依赖于均值，反映了同方差误差结构。该参数量化了水文观测值围绕模型预测的噪声或不确定水平。这一假设对于连续水文预测问题是有效的，其中目标变量在预处理后近似正态分布。因此，公式10中的线性预测子和公式12中的WLS公式均在高斯GLM的设置下。这意味着模型简化为稳定的岭型回归框架，与更一般的GLM公式一致。需要注意的是，η1（公式8中）与GRM中使用的正则化参数η2不同。η1在基于核的学习之前进行特征加权阶段起作用，而η2用于估计回归系数β的阶段。因此，η1影响输入特征缩放，而η2控制最终回归模型的稳定性。最后，GRM预测值由以下公式给出：???????????=?????（公式13）其中???????????是GRM方法计算的预测值。WARP模型：在当前研究中，所提出的WARP模型通过将上述各小节中描述的GPR、FRM和GRM集成到统一的预测框架中而开发而成。

标准GPR提供概率非线性回归模型，而FRM通过估计特征相关性并在核函数构建之前对输入变量进行缩放来增强GPR。与此同时，GRM作为互补的线性预测模型，用于捕捉基于核函数的GPR可能无法充分表示的线性关系。基于这些组件，所提出的WARP通过基于残差噪声估计的自适应核函数细化和结合改进后GPR与GRM预测结果的自适应融合策略，进一步改进了预测过程。因此，WARP提供了一个统一的框架，能够同时在水文干旱预测中建模非线性关系、线性结构、变化的噪声水平以及预测不确定性。所提出的WARP框架从三个主要方面解决了传统水文预测模型的局限性：1.基于残差的自适应核函数更新；2.集成在GPR内部的特征相关性加权；3.GPR与GRM预测的不一致性感知融合。

WARP的实现遵循一个顺序程序。首先，通过FRM根据公式8估计特征重要性向量，所获得的权重用于根据公式9构建缩放后的输入矩阵。然后，利用缩放后的输入变量构建GPR模型的协方差矩阵，使得核函数赋予信息量更大的预测因子更大的权重。GPR模型使用加权输入进行训练，其预测结果根据公式7获得。随后，观测值与预测值之间的残差用于估计噪声方差，该方差通过自适应核函数更新策略纳入协方差矩阵，该策略的公式如下：

K_RF^WARP = K_RF + NoVar × I （公式14）

在公式14中，NoVar表示从GPR模型的训练残差中估计的噪声方差。其计算方式为观测值与预测值之间残差的均方误差（MSE）：

NoVar = (1/N) × Σ_{i=1}^{N} (y_i - ?_i)2 （公式14.1）

其中y_i表示观测值，?_i是从GPR模型获得的对应预测值，N为训练样本数量。该公式基于这样的假设：噪声方差可以通过残差误差的方差来近似，假设噪声服从高斯分布。估计得到的NoVar随后代入公式14的自适应核函数更新中，以增强对水文数据噪声和不确定性的鲁棒性。

同时，使用相同的加权输入变量训练GRM，以提供捕捉输入输出关系中线性成分的额外预测结果。最后，所提出的WARP模型通过两个缩放参数（a和b）将改进后的GPR预测值与GRM预测值进行组合，其公式如下：

y^WARP = a × y^GPR + (1-a) × y^GRM + b × |y^GPR - y^GRM| （公式15）

在公式15中，a和b是用于组合GPR和GRM模型预测结果的自适应缩放系数。这些参数通过本研究采用的EMRUN优化算法进行超参数优化。系数a控制GPR和GRM输出的相对贡献，而b调节不一致项|y^GPR - y^GRM|的影响，该不一致项考虑了两种模型之间的预测不确定性。两个参数均约束在范围a, b ∈ [0, 1]内。这种自适应公式使所提出的WARP框架能够动态地平衡预测精度与不确定性表示。

可以得出结论，基于残差的噪声方差估计、加权预测和更新后的核矩阵三者的结合，使所提出的模型对不同数据集具有高度适应性。它可以处理具有以下特征的数据集：
• 具有不同噪声水平的数据集
• 同时具有线性和非线性关系的数据集
• 具有复杂模式、需要在GPR与GRM之间取得平衡的数据集

WARP模型的示意图如图2所示。

**EMRUN优化方法**

在本研究中，提出了一种改进的多算子龙格-库塔（EMRUN）优化方法，用于有效识别机器学习模型的最优参数。所建议的EMRUN使用多种高级算子和自适应机制来提高解的质量、保持种群多样性并加快收敛过程。这些机制包括一种基于排名的选择方法来挑选感兴趣的候选解、用于动态参数修改的自适应更新规则，以及一种改进的增强解质量（IESQ）算子以提高优化器的探索与开发能力。同时，序列二次规划（SQP）和协方差矩阵适应进化策略（CMA-ES）方法被嵌入框架中，以提高局部搜索性能和优化精度。EMRUN算法的主要组成部分和操作程序在附录B中展示。

**特征选择**

本研究开发了一种高效的特征选择（FS）方法，该方法基于核岭回归（KRidge）与置换重要性（PI）方法的集成。这种组合同时利用了KRidge模型在建模复杂关系方面的优势以及PI在评估特征重要性方面的鲁棒性。基于上述说明，所提出的特征选择方法被称为KRidge-PI模型，并通过EMRUN算法进行优化。以下子节将描述所提出的FS方法及其背后的逻辑。

**使用KRidge模型进行训练阶段**

KRidge模型使用所有特征进行训练。此步骤涉及将模型拟合到训练数据，并基于核函数（例如小波核）确定系数a。系数a根据公式16计算：

a = (K_RF + ? × I)?1 X? y （公式16）

其中

K_RF(x_i, x_j) = cos(2π · a? · (x_i - x_j) / a?) · exp(-‖x_i - x_j‖2 / (2a?)) （公式16.1）

其中?表示正则化因子。K_RF是具有三个控制参数（a?、a?和a?）的小波核函数。

**使用EMRUN进行超参数优化**

KRidge模型的超参数使用EMRUN算法进行优化。考虑了一个四维决策向量x = [?, a?, a?, a?]，其中每个变量均约束在范围[0, 200]内。对于EMRUN生成的每个候选解，训练KRidge模型并在测试数据集上进行评估。优化过程的目标是最小化均方根误差（RMSE），计算如下：

RMSE = √((1/N) × Σ_{i=1}^{N} (y_i - ?_i)2) （公式17）

其中N表示测试样本数量，y_i是观测值，?_i是对应的模型预测值。EMRUN算法以种群规模20和最大迭代次数100实施。EMRUN获得的最优参数向量随后用于训练最终的KRidge模型。优化后，应用置换重要性分析对输入特征进行排序，并保留前50%的特征。

**基线性能评估**

在训练KRidge模型后，使用测试数据集评估基线性能，以建立特征重要性的参考基准。为建立基线性能，采用RMSE指标。基线性能定义为：

Per_baseline = RMSE(y_m, y_p) （公式18）

其中Per_baseline表示基线性能，y_m和y_p分别为SPEI12的测量值和预测值。

**置换数据上的性能评估**

对每个特征应用PI方法。对于数据集中的每个特征，该方法对该特征在所有数据中的值进行置换。由于时间序列数据具有顺序性，因此需要注意如何处理时间序列值。因此，应考虑以下步骤。在置换特征后，使用RMSE指标评估模型性能（公式19）：

Per_permuted = RMSE(y_m, y_{p,permuted}) （公式19）

其中y_{p,permuted}表示置换后输入数据集的预测值。

**计算重要性分数**

每个特征的重要性分数（IS）可以通过基线性能与置换该特征后的性能之差来计算，其公式为：

IS = Per_baseline - Per_permuted （公式20）

基于上述关系，较大的IS值表明置换该特征会导致模型性能产生更大的变化，意味着该特征对模型预测具有更强的影响，更为重要。

为了获得鲁棒的特征重要性估计，对每个特征的置换过程重复多次，最终的重要性分数计算为所得IS值的平均值。然后根据重要性分数对特征进行排序，得分最高的特征被认为是最具影响力的预测因子。

将PI应用于时间序列数据时需要仔细保留时间依赖性。与独立观测不同，连续的时间序列样本之间存在相关性，随机置换整个数据集中的某个特征可能会破坏时间结构并产生不真实的输入模式。这种破坏可能导致特征重要性估计出现偏差。

为解决此问题，可以在预定义的时间窗口内而非在整个数据集范围内对特征值进行置换。这种方法在扰动选定特征与目标变量之间关系的同时保留了局部时间特征。每次置换后，重新评估KRidge模型性能，并使用RMSE的变化量来量化特征重要性。通过在置换过程中保持时间结构，所提出的PI程序为时间序列预测问题提供了更可靠的重要性估计，并有助于识别对模型性能贡献最大的预测因子。

**分解方法**

在本研究中，采用多变量变分模态分解（MVMD）方法将输入变量分解为多个具有不同频率特性的本征分量，从而提高了预测模型的特征提取能力。采用MVMD方法是因为干旱相关时间序列的特性——这些序列通常是非线性、非平稳的，并受到不同时间尺度上多种相互关联的气象和水文现象的影响。与EMD和EEMD相比，MVMD具有更优越的模态分离能力和更低的模态混叠，同时能够处理多个相关变量。与需要预先指定母小波的WT不同，MVMD可以直接从数据进行自适应分解。此外，MVMD在分解过程中保持了多变量信号之间的相互依赖性，有助于更好地刻画干旱驱动因素之间的物理联系。

这些优势使MVMD特别适合非线性、非平稳时间序列的干旱预测及其他水文应用。MVMD在水文、环境和预测研究中也被报道具有类似的优势。然而，MVMD的性能高度依赖于其控制参数的适当选择。因此，将EMRUN算法与MVMD框架集成，以自适应优化分解参数并增强搜索过程中探索与开发的平衡。EMRUN优化器用于自动确定MVMD主要参数的最优值：本征模态函数数量（NIM）和二次惩罚因子（?）。因此，所提出的分解方法被称为增强型多变量变分模态分解（EMVMD）。

为评估分解质量，采用基于包络熵（EnE）的适应度函数。包络熵量化了分解后本征模态函数（IMF）的复杂性和分布特征。首先，使用希尔伯特变换获取每个IMF的包络信号。将包络幅度归一化为概率分布，然后从该分布计算香农熵。较低的包络熵表明信号中包含的信息集中在较少数量的主导结构内，从而产生更有组织和信息更丰富的分解，减少了噪声和模态混叠。相反，较高的熵值表明信号特征分布更为分散和不规则，表明分解效果较差。

在优化过程中，EMRUN通过最小化所有分解IMF的平均包络熵来搜索最优的本征模态函数数量（NIM）和二次惩罚因子（?）。因此，优化的目的是获得复杂度更低、模态分离更好且信号表示更优的分解。这样的分解预计能够保留原始信号的关键信息内容并减少噪声和模态混叠的影响，从而提高后续机器学习模型的预测性能。相应地，适应度函数定义为：

F = (1/M) × Σ_{j=1}^{M} EnE(j) （公式21）

其中

EnE_prds(i) = -Σ_{j=1}^{N} prds_j × log?(prds_j) （公式22）

prds_j = SeS(j) / Σ_{j=1}^{N} SeS(j) （公式23）

SeS(j) = √[x(j)]2 + {HiT[x(j)]}2 （公式24）

其中F表示适应度函数，M为分解得到的IMF总数，IMF的包络熵由EnE_prds表示，N表示信号长度，prds表示信号的概率分布，SeS(i)是通过希尔伯特变换（HiT）获得的包络信号序列。HiT[•]表示希尔伯特变换。希尔伯特变换用于构建每个IMF的解析表示并提取其瞬时包络幅度。对于信号x(j)，希尔伯特变换定义为：

HiT[x(t)] = x?(t) = (1/π) × P.V. ∫_{-∞}^{∞} x(τ)/(t-τ) dτ （公式24.1）

其中P.V.表示柯西主值，由于被积函数在τ = t处奇异，因此采用主值。对于希尔伯特变换，主值定义为：

P.V. ∫_{-∞}^{∞} x(τ)/(t-τ) dτ = lim_{ε→0} [∫_{-∞}^{t-ε} x(τ)/(t-τ) dτ + ∫_{t+ε}^{∞} x(τ)/(t-τ) dτ] （公式24.2）

其中ε是一个小的正数，用于排除奇点τ = t，极限ε→0确保积分的适当计算。

因此，完整的希尔伯特变换（HiT）可以表示为：

$$\text{HiT}[x(t)] = \frac{1}{\pi} \lim_{\varepsilon \to 0} \left[ \int_{-\infty}^{t-\varepsilon} \frac{x(\tau)}{t-\tau} d\tau + \int_{t+\varepsilon}^{\infty} \frac{x(\tau)}{t-\tau} d\tau \right] \quad \text{(公式24.3)}$$

因此，公式22用于计算单个本征模态函数（IMF）的包络熵，而公式21则通过对所有IMFs的包络熵值取平均来评估整体分解质量。（F）的值越低，表明分解越有效，具有更好的模式分离度和更低的复杂性。以下是用于分解输入变量的EMVMD过程的描述：

1. NIM的搜索范围在区间[2, 12]内初始化，二次惩罚因子（?）定义在[100, 3000]范围内。相应地，优化问题包含两个与EMVMD分解过程相关的决策变量。EMRUN算法的种群规模（Np）设为20，最大迭代次数（MIt）设为50，以在计算效率和优化精度之间取得平衡。

2. EMRUN算法随机生成初始候选解，对应于EMVMD控制参数（NIM和?）。然后，MVMD方法应用这些参数值将输入变量分解为具有不同频率特征的若干本征模态函数。

3. 使用公式21计算适应度函数，用于评估分解信号的质量。包络熵（EnE）用于衡量所生成的本征模态函数的复杂性、稳定性和有效性。

4. 然后检查停止准则，即优化过程是否已达到最大迭代次数。如果当前迭代次数（It）等于最大迭代次数（MIt），则停止优化过程并转到步骤（7）。否则，算法转到步骤（5）以更新当前候选解并创建新解。

5. 通过EMRUN优化算子更新候选解，包括增强的搜索机制、自适应参数更新、交叉策略和IESQ算子，以提高优化器的探索和开发性能。

6. 然后，算法返回步骤（3），重新计算新解的适应度函数并评估其分解性能。

7. 最后，获得NIM和?的优化值。EMVMD框架使用这些优化参数有效地分解输入变量，并为后续预测模型提供高质量的分解信号。

## 计算标准化降水蒸发指数

标准化降水蒸发指数（SPEI）是最广泛使用的干旱指数之一，用于确定干旱。SPEI使用降雨量和潜在蒸散发来计算干旱期。SPEI可获得1至24个月的多时间步长。

在这方面，为了确定干旱严重程度，根据公式25计算潜在蒸散发（PET）与降水量（P）之间的差值：

$$DP_i = PET_i - P_i \quad \text{(公式25)}$$

其中DPi表示PET与P参数之间的差值。本研究使用Hargreaves和Samani（1985）方法来确定蒸散发量。SPEI通过将DP序列拟合到三参数对数逻辑概率分布函数（LPDF）获得，该函数由于能够表示偏态水文气候数据而在水文和气候学研究中得到广泛应用。公式26用于定义对数逻辑分布的概率密度函数（PDF），从中可推导出累积分布函数（CDF）。

$$z(x) = \frac{\vartheta}{\sigma} \left(\frac{x-\varsigma}{\sigma}\right)^{\vartheta-1} \left[1 + \left(\frac{x-\varsigma}{\sigma}\right)^{\vartheta}\right]^{-2} \quad \text{(公式26)}$$

其中σ、?和?表示LPDF的主要参数，x表示累积序列。LPDF的参数使用L矩方法估计，该方法为水文时间序列提供稳健且稳定的参数估计，并在原始SPEI公式中推荐采用。为了确定CDF的主要参数，使用L矩方法：

$$\vartheta = \frac{2w_1 - w_0}{6w_1 - w_0 - 6w_2} \quad \text{(公式27)}$$

$$\sigma = \frac{(w_0 - 2w_1) \cdot \vartheta}{\Gamma(1+1/\vartheta) \cdot \Gamma(1-1/\vartheta)} \quad \text{(公式28)}$$

$$\varsigma = w_0 - \frac{\sigma \cdot \Gamma(1+1/\vartheta) \cdot \Gamma(1-1/\vartheta)}{\vartheta} \quad \text{(公式29)}$$

其中w0、w1和w2表示三个概率矩权重，Γ(?)表示伽马函数。基于LPDF，DP的概率密度函数表示如下：

$$Z(x) = \left[1 + \left(\frac{x-\varsigma}{\sigma}\right)^{\vartheta}\right]^{-1} \quad \text{(公式30)}$$

根据PDF（Z(x)），通过Z(x)的标准化量可获得不同时间步长的SPEI，表示如下：

$$SPEI = \frac{Q - \varepsilon_0 + \varepsilon_1 Q + \varepsilon_2 Q^2}{1 + r_1 Q + r_2 Q^2 + r_3 Q^3} \quad \text{(公式31)}$$

其中

$$Q = \sqrt{-2 \ln(EP)}, \quad \text{当} \; EP \leq 0.5 \quad \text{(公式32)}$$

$$EP = 1 - Z(x), \quad \text{当} \; EP > 0.5 \quad \text{(公式33)}$$

其中EP表示DP序列的超越概率。ε0、ε1、ε2、r1、r2和r3的参数值分别为2.515517、0.802853、0.010328、1.432788、0.189269和0.001308。本研究使用SPEI12来确定伊朗三个监测站的月度干旱变化。表11展示了SPEI的常用分类标准。

| SPEI范围 | 干旱状态 |
|---------|---------|
| SPEI > 2.0 | 极度湿润 |
| 1.5 < SPEI ≤ 2.0 | 中度湿润 |
| 1.0 < SPEI ≤ 1.5 | 接近正常 |
| ?1.5 < SPEI ≤ ?1.0 | 中度干旱 |
| ?2.0 < SPEI ≤ ?1.5 | 严重干旱 |
| SPEI ≤ ?2.0 | 极度干旱 |

表11 SPEI分类

## 基于平稳自助法的不确定性量化

预测模型中的不确定性通过平稳自助法对模型残差进行估计。该方法适用于时间序列数据，因为它考虑了模型误差通常在时间上具有相关性，而非从一个观测到下一个观测相互独立的事实。许多环境时间序列的残差表现出持续性，即某一时间步长偏离模型后，后续步骤往往也会出现类似的偏离。标准自助法在重采样残差时不考虑这种依赖性结构，从而打乱了时间顺序，可能导致对不确定性的不真实表现。

平稳自助法通过以任意大小的连续观测块重采样残差来解决这一局限性，同时保持误差结构中的局部自相关性。块长度按概率方式生成，提供了灵活性以记录不同时期的依赖性，同时保持整体平稳性。该重采样程序保留了原始误差序列的主要时间特性，有助于生成更真实的预测误差合成实现，从而更准确地估计预测不确定性。

设yt表示时间t的观测SPEI值，$\hat{y}_t$表示相应的模型预测值。残差定义为：

$$e_t = y_t - \hat{y}_t, \quad t = 1, 2, \ldots, n \quad \text{(公式34)}$$

其中et表示时间t的模型误差，n是观测总数。残差序列表示为E = {e1, e2, …, en}。

### 平稳自助法程序

平稳自助法用于考虑残差中的自相关性。该方法通过以随机大小的连续残差块重采样生成自助样本，从而保留局部时间依赖性。开始新块的概率定义为：

$$p = \frac{1}{L} \quad \text{(公式35)}$$

其中L为平均块长度。在本研究中，选取L=12以反映月度数据并保持残差结构中的年际依赖性。使用该方案，总共生成B=10,000个自助重复样本。对于每个自助重复样本b∈{1, …, B}，使用平稳自助法获得重采样的残差序列：

$$E^*(b) = \{e^*_1(b), \; e^*_2(b), \; \ldots, \; e^*_n(b)\} \quad \text{(公式36)}$$

每个序列通过随机块重采样构建，同时保持局部时间依赖性。

### 自助预测集合

对于每个自助重复样本b，通过将重采样的残差加到原始模型预测上来生成合成预测序列：

$$\hat{y}^{*(b)}_t = \hat{y}_t + e^{*(b)}_t \quad \text{(公式37)}$$

其中$\hat{y}_t$为模型预测值，$e^{*(b)}_t$为从第b个平稳自助样本中获得的残差。对B=10,000个自助重复样本重复此过程，在每个时间步产生一组可能的预测集合。该集合代表与模型误差相关的不确定性，同时保留其时间依赖性结构。自助预测的分布为估计预测区间和量化预测不确定性提供了基础。

### 预测区间

预测区间通过对自助预测集合取百分位数来构建。在每个时间t，95%预测区间的下界和上界分别为：

$$\text{LowerBound}_t = Q_{0.025}(\hat{y}^{*(1)}_t, \; \hat{y}^{*(2)}_t, \; \ldots, \; \hat{y}^{*(B)}_t) \quad \text{(公式38)}$$

$$\text{UpperBound}_t = Q_{0.975}(\hat{y}^{*(1)}_t, \; \hat{y}^{*(2)}_t, \; \ldots, \; \hat{y}^{*(B)}_t) \quad \text{(公式39)}$$

其中Q0.025和Q0.975分别为自助预测集合的2.5%和97.5%百分位数。这些百分位数边界代表自助生成的预测的中间95%部分。由此产生的区间反映了残差变异性所导致的不确定性，提供了未来观测值可能出现的区间估计，依据为拟合模型和观测到的误差结构。

### 区间评估指标

预测区间的质量通过三个公认指标进行评估，即预测区间覆盖概率（PICP）、预测区间平均宽度（MPIW）和预测区间归一化平均宽度（PINAW）。这些指标用于评估所创建预测区间的可靠性和尖锐度。PICP为观测值落在预测区间内的比例，定义为：

$$PICP = \frac{1}{n} \sum_{t=1}^{n} c_t \quad \text{(公式40)}$$

其中

$$c_t = \begin{cases} 1, & \text{若} \; y_t \in [L_t, \; U_t] \\ 0, & \text{否则} \end{cases} \quad \text{(公式41)}$$

其中Lt和Ut分别为时间t预测区间的下界和上界。PICP值越高，表明区间的可靠性越好。MPIW反映预测区间的宽度，定义为：

$$MPIW = \frac{1}{n} \sum_{t=1}^{n} (U_t - L_t) \quad \text{(公式42)}$$

MPIW值越小，表明预测区间越尖锐、信息越丰富。为了在不同量级的数据集之间进行更易比较，将MPIW除以观测值的范围以得到PINAW：

$$PINAW = \frac{MPIW}{y_{\max} - y_{\min}} \quad \text{(公式43)}$$

其中ymax和ymin分别为观测值的最大值和最小值。PINAW值越小，意味着相对于数据范围，预测区间越小，因此预测精度越高。PICP用于评估区间的可靠性，MPIW和PINAW用于评估区间的尖锐度和效率。理想的预测区间应具有高PICP值和低MPIW及PINAW值。

## 计算散点图的上界和下界

预测标准误差（SEP）衡量模型的预测不确定性，通过将每个数据点的影响度（杠杆）和模型的残差误差结合计算得出。杠杆表示每个观测值在模型中的相对影响程度，基于其与输入数据集平均值的距离；离平均值越远的数据点影响越大，因此对预测不确定性的贡献也越大。每个数据点的杠杆为：

$$Leverage_i = \frac{(X_i - \bar{X})^2}{\sum_{i=1}^{N} (X_i - \bar{X})^2} \quad \text{(公式44)}$$

其中Xi表示第i个数据点，$\bar{X}$表示数据集的均值。Leveragei表示第i个数据点的杠杆值。

模型残差误差计算为观测值与预测值之间的均方误差（MSE），定义为：

$$MSE = \frac{1}{N} \sum_{i=1}^{N} (y_i - y_{p_i})^2 \quad \text{(公式45)}$$

其中yi表示实测值，ypi表示预测值。

基于杠杆和残差误差，每个数据点的预测标准误差计算如下：

$$SEP_i = SR \times \sqrt{1 + LV_i} \quad \text{(公式46)}$$

为了构建预测区间（上界和下界），首先定义观测值-预测值关系的45度参考线（y45）。然后使用t分布计算上下界如下：

$$\text{Upper} = y_{45} + t_{\alpha/2, \; df} \times SEP_i \quad \text{(公式47)}$$

$$\text{Lower} = y_{45} - t_{\alpha/2, \; df} \times SEP_i \quad \text{(公式48)}$$

其中$t_{\alpha/2, \; df}$表示t值，α表示显著性水平，df表示自由度。

## EMRUN算法

EMRUN算法的主要算子在以下各节中描述。

### 搜索机制

EMRUN算法使用一种新的增强搜索机制（ESM），通过有效的基于排名的选择策略来有效平衡全局探索与局部开发阶段。该机制增强了优化器寻找有前景搜索区域的能力，同时保持种群多样性并加速收敛至最优解。所提出的ESM数学表达如下：

$$Y_1 = \begin{cases} (x_{c1} + r \times A \times x_c) + A \times (SM + 0.5 \times \rho \times (x_{c2} - x_{c1})), & \text{若} \; rand > \theta(ID_k) \\ (x_{c2} + r \times A \times x_{c2}) + A \times (SM + 0.5 \times \rho \times (x_{a1} - x_{a2})), & \text{否则} \end{cases} \quad \text{(公式49)}$$

其中r表示取值为1或?1的随机整数，xa1和xa2为从[1, Np]范围内随机选择的两个解，Np表示种群规模，A表示自适应参数，ρ表示服从正态分布的随机值。SM在第a.1.2节中有详细描述。xc1和xc2由公式49.1和49.2给出：

$$x_{c1} = \nu \cdot x_M + (1 - \nu) \cdot x_{a1} \quad \text{(公式49.1)}$$

$$x_{c2} = \nu \cdot x_{bs} + (1 - \nu) \cdot x_{pbs} \quad \text{(公式49.2)}$$

其中

$$\rho = 0.5 + 0.1 \cdot randn \quad \text{(公式49.3)}$$

$$A = 4 \cdot (0.5 - rand) \times e \quad \text{(公式49.4)}$$

$$e = 2 \times \exp\left(-8 \cdot rand \cdot (I_{term}^{max} - I_{term})\right) \quad \text{(公式49.5)}$$

以及

$$\theta(ID_k) = 1 - \cos\left(Rank(ID_k) \cdot \frac{\pi}{2}\right) \quad \text{(公式49.6)}$$

其中

$$Rank(ID_k) = \left[\frac{N_p - k}{N_p}\right]^2 \quad \text{(公式49.7)}$$

其中θ表示每个解的排名，IDk表示解根据适应度函数排序后的索引，Rank表示种群中每个成员的总体排名。

此外，为促进解的多样性，还实施了交叉算子。

公式50定义了通过交叉算子生成解 (xck) 的过程：

???????={??1,?? ???? ????????
其中：

????????=τ??????,??+0.1×?????????? （公式50.1）

其中，crsk 表示交叉概率，krd 表示区间 [1, D] 内的随机整数，τcrs,k 表示初始迭代中等于 0.5 的数值。τcrs,k 可通过公式50.2获得：

τ??????,??=σ·τ??????,??+(1?σ)·∑????????(????????)??=1 ?????????? / ????????(????????) （公式50.2）

其中 σ 表示 0.9 的值；Scsr 表示每次迭代中 crsk 的所有成功值。

**搜索机制（SM）**

SM 可表示为：

????=16·(????????)·Δ?? （公式51）

????????=??1+2×??2+2×??3+??4 （公式51.1）

??1=1/(2Δ??)·(????????×???????λ·??????) （公式51.2）

λ=??????????(1+????????)·(1?????????) （公式51.3）

??2=1/(2Δ??)·(????????×(??????+??1·??1·Δ??)?(λ·??????+??2·??1·Δ??)) （公式51.4）

??3=1/(2Δ??)·(????????×(??????+??1·(1/(2??2))·Δ??)?(λ·??????+??2·(1/(2??2))·Δ??)) （公式51.5）

??4=1/(2Δ??)·(????????×(??????+??1·??3·Δ??)?(λ·??????+??2·??3·Δ??)) （公式51.6）

其中，?????? 和 ?????? 分别表示最劣解和最优解，??1 和 ??2 表示 [0, 1] 范围内的两个随机数。Δ?? 定义为：

Δ??=2·????????·|??| （公式51.7）

??=????????×[(???????????????·????)+χ] （公式51.8）

χ=????????×(?????????????×(?????????))·exp(?4·????/??????) （公式51.9）

其中，?? 表示步长，???? 和 ???? 分别表示可行域的上界和下界，???? 表示解的均值，???? 表示迭代次数，?????? 表示总迭代次数。?????? 和 ?????? 基于公式51.10确定：

???? ??(????) ????????=????
????????=????????
else
????????=????????
????????=????
end （公式51.10）

其中，???????? 表示在三个随机解（????1, ????2, 和 ????3）中获取的最优解。???????? 和 ???????? 分别表示局部最优解和局部最劣解。

**IESQ算子**

IESQ算子被嵌入到EMRUN算法中，以提高候选解的质量并降低过早收敛到局部最优解的可能性。所提出的EMRUN方法采用从种群平均解和随机选取的解中生成的两个差分向量，以增强探索和开发能力。因此，利用IESQ机制根据公式52生成新候选解??2：

???? ????????<(1?????/??????)2 或 ????<0.1×??????
??2={????1+??·(?????????)·?? 如果 ????????<ε； ????+????·(?????????)+0.5·????·(????????1?????????2) 其他情况} end （公式52） 其中： ????=??×????????1+(1???)×???? （公式52.1） ??={1, 如果><ε；0 其他情况}；ε=0.5+0.1×?????????? （公式52.2）

其中，?? 是一个逻辑参数，用于指定向量 (????-????) 的哪些部分可用于生成 ??2。???? 表示由 ????????1 和 ???? 组合创建的解。????????1 和 ????????2 分别表示较好的随机解和无能随机解。为确定这些解，首先将种群分为两个相等的部分，第一部分包含排名较高的解，第二部分包含排名较低的解。因此，????????1 和 ????????2 分别从第一部分和第二部分中随机选取。???? 表示自适应参数。为指定 ????，需要计算解 ????2,??。相应地，解 ????2,?? 可计算为：

????2,??={??2,?? 如果 ????????
其中：

??????2,??=τ??????2,??+0.1·?????????? （公式53.1）

其中，??????2,?? 表示交叉概率，τ??????2,?? 表示初始迭代中等于 0.5 的数值。τ??????2,?? 由公式53.2确定：

τ??????2,??=σ·τ??????2,??+(1?σ)·∑????????(????????2)??=1 ????????2,?? / ????????(????????2) （公式53.2）

其中，????????2 表示每次迭代中所有成功的 ??????2,?? 值。

自适应参数 ???? 定义如下：

????=?????,??+0.1×tan(0.5·π·??????????) （公式53.3）

其中：

?????,??=σ·?????,??+(1?σ)·∑????????(??????)??=1 ω′·??????,??2 / ∑????????(??????)??=1 ω′·??????,?? （公式53.4）

ω′=|??(????2,??)???(????)| / ∑????????(??????)??=1 |??(????2,??)???(????,??)| （公式53.5）

其中，?????? 表示每次迭代中所有成功的 ???? 值，ω′ 表示加权因子。

由于解 ????2,?? 可能没有比 ???? 更好的适应度函数值 (??(????))，因此生成解 ?? 以给优化器第二次生成合适解的机会，其公式为：

???? ????????<0.3
??=????+????·(???????????)+0.5·????·(????????1?????????2)
end （公式54）

emrun算法采用动态种群规模策略，以提高早期迭代的探索能力并增强优化过程后期的开发性能。这种自适应机制使优化器能够在搜索开始时保持种群多样性，并逐渐聚焦于有前景的区域，以加快向最优解的收敛速度。所提出的动态种群调整机制可表示为：

??????=??????????[????????+(?????????????????)·(????/??????)^0.5] （公式55）

其中，???????? 和 ???????? 分别表示 ???? 的最小值和最大值。

**模型性能评估准则**

本研究使用六项统计准则评估所有机器学习模型，包括均方根误差（rmse）、相关系数（r）、平均绝对百分比误差（mape）、nash-sutcliffe效率系数（nse）、最大绝对误差（maxae）以及willmott一致性指数（ia）。这些指标公式如下：

??=∑????=1(??????????,??????????????)×(??????????,??????????????) / √[∑????=1(??????????,??????????????)2×∑????=1(??????????,??????????????)2] （公式56）

????????=(1/??)·∑????=1 |(??????????,?????????????,??)/??????????,??|×100 （公式57）

??????????=????????=1,…,?? |??????????,?????????????,??| （公式58）

????????=√[(1/??)·∑????=1(??????????,?????????????,??)2] （公式59）

??????=1?∑????=1(??????????,?????????????,??)2 / ∑????=1(??????????,??????????????)2 （公式60）

????=1?∑????=1(??????????,?????????????,??)2 / ∑????=1(|??????????,??????????????|+|??????????,??????????????|)2, 0
其中，??????????,?? 和 ??????????,?? 分别为实测和预测的spei值。??????????? 和 ??????????? 分别为实测和预测spei值的平均值。

**案例研究**

本研究使用伊朗不同水文气候区的三个气象站——亚兹德、巴博尔斯尔和乌鲁米耶——的气象数据，开发并评估了六种用于干旱预测的机器学习模型。亚兹德站（东经54°20′，北纬31°52′，海拔1230米）位于伊朗中部，具有炎热干旱的气候；巴博尔斯尔站（东经52°39′，北纬36°42′，海拔-21米）位于伊朗北部，受里海影响的湿润气候；乌鲁米耶站（东经45°03′，北纬37°32′，海拔1328米）位于伊朗西北部，具有半干旱气候。这些气象站的选择旨在确保所开发的模型能够在不同的水文气候条件下进行评估。

数据集包含从1987年到2019年共33年的连续月度观测数据，用于计算12个月时间尺度的标准化降水蒸散发指数（spei12）。亚兹德、巴博尔斯尔和乌鲁米耶的spei12指数序列分别表示为spei-y、spei-b和spei-u。数据集包含不同形式的干旱和湿润条件，为训练和测试机器学习模型提供了足够的变异性。表9显示了所选气象站的气候特征，表10总结了spei12序列的统计属性。此外，spei序列的时间变化如图3所示，三个气象站的地理位置如图4所示。

**表9 三个气象站的气候信息**

| 气象站 | 纬度 | 经度 | 海拔（米） | 平均降水量（p）（毫米） | 平均气温（t）（°c） | 平均相对湿度 | de martonne干旱指数（p/t+10） |
|--------|------|------|-----------|----------------------|-------------------|------------|------------------------------|
| 乌鲁米耶 | 37°32′n | 45°03′e | 1328 | 310 | 11.7 | 58 | 14.5 |
| 巴博尔斯尔 | 36°42′n | 52°39′e | -21 | 911 | 17.8 | 80 | 33 |
| 亚兹德 | 31°52′n | 54°20′e | 1230 | 50 | 20.2 | 28 | 1.7 |

**表10 三个气象站spei12的统计信息**

| 序列 | 最大值 | 最小值 | 平均值 | 标准差 | 偏度 | 峰度 | 25%分位数 | 50%分位数 | 75%分位数 |
|------|--------|--------|--------|--------|------|------|----------|----------|----------|
| spei12-b | 2.466 | -2.184 | -3.06e-04 | 0.979 | 0.034 | -0.627 | -0.711 | -0.017 | 0.728 |
| spei12-u | 2.172 | -2.466 | 1.47e-02 | 0.981 | 0.097 | -0.565 | -0.739 | 0.047 | 0.622 |
| spei12-y | 2.097 | -3.090 | 5.69e-03 | 0.986 | -0.061 | -0.244 | -0.625 | -0.047 | 0.710 |

spei-12指标因其与水文干旱指标显著相关而被采用。每个数据集使用九项统计指标进行描述，包括最大值（max）、最小值（min）、平均值（mean）、标准差（sd）、偏度（skew）、峰度（kurt）以及分位数（25th q、50th q和75th q）。从表中可以看出，spei12-b具有最高的最大值（2.466），其次是spei12-u（2.172）和spei12-y（2.097）。这表明与spei12-u和spei12-y相比，spei12-b具有最大的正极端值（exv）。此外，spei12-y具有最低的最小值（-3.090），其次是spei12-u（-2.466）和spei12-b（-2.184）。因此，spei12-y具有最负的极端值。

在偏度方面，spei12-b（0.034）和spei12-u（0.097）具有轻微的正偏度，表明数据略微倾向于在右侧（正侧）具有较长的尾部。同样，spei12-y（-0.061）具有轻微的负偏度，表明略微倾向于在左侧（负侧）具有较长的尾部。在分位数分析方面，spei12-u具有最低的25%分位数（-0.739），表明在较低范围内的负值最多，而spei12-b的50%分位数最高（-0.017），表明中心分布更加平衡。对于75%分位数，spei12-b（0.728）也显示出最大的正值，其次是spei12-y（0.710）和spei12-u（0.622）。

**量化与统计分析**

在预测亚兹德、巴博尔斯尔和乌鲁米耶站的spei12时，所提出的warp框架始终优于lssvm、drvfl、delm、cfnn和lasso，基于各种统计指标显示出最高的精度。其具有非常高的相关系数（训练集r=0.988–0.996，测试集r=0.987–0.990）和最低的最大误差，表明在减少极端预测误差方面具有更好的能力。此外，图形评估和不确定性分析也证实了其稳健性，与所有其他竞争性模型相比，其呈现出不确定性区间最窄（-0.601至0.078）。总之，warp是一种精确可靠的干旱预测工具，在干旱管理和未来气候变化影响估计方面具有很高的潜力。

**附加资源**

本研究没有相关的附加资源。 其他情况}；ε=0.5+0.1×?????????? （公式52.2） 其中，?? 是一个逻辑参数，用于指定向量 (????-????) 的哪些部分可用于生成 ??2。???? 表示由 ????????1 和 ???? 组合创建的解。????????1 和 ????????2 分别表示较好的随机解和无能随机解。为确定这些解，首先将种群分为两个相等的部分，第一部分包含排名较高的解，第二部分包含排名较低的解。因此，????????1 和 ????????2 分别从第一部分和第二部分中随机选取。???? 表示自适应参数。为指定 ????，需要计算解 ????2,??。相应地，解 ????2,?? 可计算为： ????2,??={??2,?? 如果 ????????
其中，?? 是一个逻辑参数，用于指定向量 (????-????) 的哪些部分可用于生成 ??2。???? 表示由 ????????1 和 ???? 组合创建的解。????????1 和 ????????2 分别表示较好的随机解和无能随机解。为确定这些解，首先将种群分为两个相等的部分，第一部分包含排名较高的解，第二部分包含排名较低的解。因此，????????1 和 ????????2 分别从第一部分和第二部分中随机选取。???? 表示自适应参数。为指定 ????，需要计算解 ????2,??。相应地，解 ????2,?? 可计算为：

????2,??={??2,?? 如果 ????????
其中：

??????2,??=τ??????2,??+0.1·?????????? （公式53.1）

其中，??????2,?? 表示交叉概率，τ??????2,?? 表示初始迭代中等于 0.5 的数值。τ??????2,?? 由公式53.2确定：

τ??????2,??=σ·τ??????2,??+(1?σ)·∑????????(????????2)??=1 ????????2,?? / ????????(????????2) （公式53.2）

其中，????????2 表示每次迭代中所有成功的 ??????2,?? 值。

自适应参数 ???? 定义如下：

????=?????,??+0.1×tan(0.5·π·??????????) （公式53.3）

其中：

?????,??=σ·?????,??+(1?σ)·∑????????(??????)??=1 ω′·??????,??2 / ∑????????(??????)??=1 ω′·??????,?? （公式53.4）

ω′=|??(????2,??)???(????)| / ∑????????(??????)??=1 |??(????2,??)???(????,??)| （公式53.5）

其中，?????? 表示每次迭代中所有成功的 ???? 值，ω′ 表示加权因子。

由于解 ????2,?? 可能没有比 ???? 更好的适应度函数值 (??(????))，因此生成解 ?? 以给优化器第二次生成合适解的机会，其公式为：

???? ????????<0.3
??=????+????·(???????????)+0.5·????·(????????1?????????2)
end （公式54）

emrun算法采用动态种群规模策略，以提高早期迭代的探索能力并增强优化过程后期的开发性能。这种自适应机制使优化器能够在搜索开始时保持种群多样性，并逐渐聚焦于有前景的区域，以加快向最优解的收敛速度。所提出的动态种群调整机制可表示为：

??????=??????????[????????+(?????????????????)·(????/??????)^0.5] （公式55）

其中，???????? 和 ???????? 分别表示 ???? 的最小值和最大值。

**模型性能评估准则**

本研究使用六项统计准则评估所有机器学习模型，包括均方根误差（rmse）、相关系数（r）、平均绝对百分比误差（mape）、nash-sutcliffe效率系数（nse）、最大绝对误差（maxae）以及willmott一致性指数（ia）。这些指标公式如下：

??=∑????=1(??????????,??????????????)×(??????????,??????????????) / √[∑????=1(??????????,??????????????)2×∑????=1(??????????,??????????????)2] （公式56）

????????=(1/??)·∑????=1 |(??????????,?????????????,??)/??????????,??|×100 （公式57）

??????????=????????=1,…,?? |??????????,?????????????,??| （公式58）

????????=√[(1/??)·∑????=1(??????????,?????????????,??)2] （公式59）

??????=1?∑????=1(??????????,?????????????,??)2 / ∑????=1(??????????,??????????????)2 （公式60）

????=1?∑????=1(??????????,?????????????,??)2 / ∑????=1(|??????????,??????????????|+|??????????,??????????????|)2, 0
其中，??????????,?? 和 ??????????,?? 分别为实测和预测的spei值。??????????? 和 ??????????? 分别为实测和预测spei值的平均值。

**案例研究**

本研究使用伊朗不同水文气候区的三个气象站——亚兹德、巴博尔斯尔和乌鲁米耶——的气象数据，开发并评估了六种用于干旱预测的机器学习模型。亚兹德站（东经54°20′，北纬31°52′，海拔1230米）位于伊朗中部，具有炎热干旱的气候；巴博尔斯尔站（东经52°39′，北纬36°42′，海拔-21米）位于伊朗北部，受里海影响的湿润气候；乌鲁米耶站（东经45°03′，北纬37°32′，海拔1328米）位于伊朗西北部，具有半干旱气候。这些气象站的选择旨在确保所开发的模型能够在不同的水文气候条件下进行评估。

数据集包含从1987年到2019年共33年的连续月度观测数据，用于计算12个月时间尺度的标准化降水蒸散发指数（spei12）。亚兹德、巴博尔斯尔和乌鲁米耶的spei12指数序列分别表示为spei-y、spei-b和spei-u。数据集包含不同形式的干旱和湿润条件，为训练和测试机器学习模型提供了足够的变异性。表9显示了所选气象站的气候特征，表10总结了spei12序列的统计属性。此外，spei序列的时间变化如图3所示，三个气象站的地理位置如图4所示。

**表9 三个气象站的气候信息**

| 气象站 | 纬度 | 经度 | 海拔（米） | 平均降水量（p）（毫米） | 平均气温（t）（°c） | 平均相对湿度 | de martonne干旱指数（p/t+10） |
|--------|------|------|-----------|----------------------|-------------------|------------|------------------------------|
| 乌鲁米耶 | 37°32′n | 45°03′e | 1328 | 310 | 11.7 | 58 | 14.5 |
| 巴博尔斯尔 | 36°42′n | 52°39′e | -21 | 911 | 17.8 | 80 | 33 |
| 亚兹德 | 31°52′n | 54°20′e | 1230 | 50 | 20.2 | 28 | 1.7 |

**表10 三个气象站spei12的统计信息**

| 序列 | 最大值 | 最小值 | 平均值 | 标准差 | 偏度 | 峰度 | 25%分位数 | 50%分位数 | 75%分位数 |
|------|--------|--------|--------|--------|------|------|----------|----------|----------|
| spei12-b | 2.466 | -2.184 | -3.06e-04 | 0.979 | 0.034 | -0.627 | -0.711 | -0.017 | 0.728 |
| spei12-u | 2.172 | -2.466 | 1.47e-02 | 0.981 | 0.097 | -0.565 | -0.739 | 0.047 | 0.622 |
| spei12-y | 2.097 | -3.090 | 5.69e-03 | 0.986 | -0.061 | -0.244 | -0.625 | -0.047 | 0.710 |

spei-12指标因其与水文干旱指标显著相关而被采用。每个数据集使用九项统计指标进行描述，包括最大值（max）、最小值（min）、平均值（mean）、标准差（sd）、偏度（skew）、峰度（kurt）以及分位数（25th q、50th q和75th q）。从表中可以看出，spei12-b具有最高的最大值（2.466），其次是spei12-u（2.172）和spei12-y（2.097）。这表明与spei12-u和spei12-y相比，spei12-b具有最大的正极端值（exv）。此外，spei12-y具有最低的最小值（-3.090），其次是spei12-u（-2.466）和spei12-b（-2.184）。因此，spei12-y具有最负的极端值。

在偏度方面，spei12-b（0.034）和spei12-u（0.097）具有轻微的正偏度，表明数据略微倾向于在右侧（正侧）具有较长的尾部。同样，spei12-y（-0.061）具有轻微的负偏度，表明略微倾向于在左侧（负侧）具有较长的尾部。在分位数分析方面，spei12-u具有最低的25%分位数（-0.739），表明在较低范围内的负值最多，而spei12-b的50%分位数最高（-0.017），表明中心分布更加平衡。对于75%分位数，spei12-b（0.728）也显示出最大的正值，其次是spei12-y（0.710）和spei12-u（0.622）。

**量化与统计分析**

在预测亚兹德、巴博尔斯尔和乌鲁米耶站的spei12时，所提出的warp框架始终优于lssvm、drvfl、delm、cfnn和lasso，基于各种统计指标显示出最高的精度。其具有非常高的相关系数（训练集r=0.988–0.996，测试集r=0.987–0.990）和最低的最大误差，表明在减少极端预测误差方面具有更好的能力。此外，图形评估和不确定性分析也证实了其稳健性，与所有其他竞争性模型相比，其呈现出不确定性区间最窄（-0.601至0.078）。总之，warp是一种精确可靠的干旱预测工具，在干旱管理和未来气候变化影响估计方面具有很高的潜力。

**附加资源**

本研究没有相关的附加资源。>