全球气候变化与人类活动的加剧，导致水资源时空分布失衡问题日益严重. 径流预测作为水资源配置与管理以及生态流量保障的核心环节，其重要性愈发显著^[1-2]. 然而，受到气温、降水等多因素耦合作用，径流预测数据呈现非平稳性、非线性和季节性特征，须有效提取复杂径流特征并融合多特征数据，以提高预测精度.

水文预测模型包括物理模型^[3-4]和数据驱动模型^[5]. 物理模型基于水文过程机理，可解释性强但计算复杂；数据驱动模型则依据历史数据构建映射关系，无须了解物理过程^[6-7]，在多特征水文预测中表现优异. 数据驱动模型可分为统计学习模型、浅层机器学习模型和深度学习模型. 传统统计模型^[8-10]依赖平稳性假设，难以捕捉非线性关系^[11]. 机器学习模型^[12-17]如支持向量机（SVM）^[18-19]仍难应对水文数据的复杂性和不确定性.

近年来，深度学习模型在径流预测中取得了显著进展. 其中，长短期记忆网络（LSTM）^[20-23]通过引入门控机制有效解决了传统循环神经网络（recurrent neural network，RNN）梯度问题，能够出色地捕捉长短期依赖关系^[24-26]. 但在处理极长序列时仍存在稳定性问题. 为此，研究者们开发了多种混合架构以融合不同模型的优势. 王万良等^[27]提出结合时间卷积网络（temporal convolutional network，TCN）和LSTM的模型能够同时捕捉局部和长期依赖关系. Dai等^[28]提出基于双阶段注意力（dual-stage attention，DSA）的多模态深度学习（mutimodal deep learning，MDL）模型DSAMDL，该模型结合了一维卷积、双向长短期记忆网络（bidirectional long short-term memory，Bi-LSTM）和DSA的优势，提高了预测的可靠性和可解释性.

尽管现有模型能够处理非线性关系，但在提取非平稳径流序列特征时仍显不足. 且其训练过程需要耗费大量时间以寻找最优解，为此，信号分解方法与优化算法被广泛应用于深度学习模型. 变分模态分解（variational mode decomposition，VMD）^[29-30]通过频域分解有效地改善了水文序列非平稳问题^[31]；Yu等^[32]提出基于自适应傅里叶分解模型（adaptive Fourier decomposition model，AFDM）和多尺度时间卷积网络（multi-scale temporal convolutional network，MTCN）的AFDM-MTCN模型，并采用麻雀搜索算法（sparrow search algorithm，SSA)优化，但SSA存在收敛速度较慢、易陷入局部最优的问题，限制了模型预测精度的进一步提升.

针对上述问题，提出结合多特征融合（multi-feature fusion，MFF）和牛顿-拉夫逊优化（Newton-Raphson-based optimizer，NRBO）算法的LSTM模型MFF-NRBO-LSTM. 运用VMD分解历史径流序列，提取不同频率成分并与其他关键特征融合，并将NRBO应用于径流预测模型的超参数优化. 最终将所提模型应用于黄河流域径流预测，并通过对比和消融验证模型的有效性.

黄河作为中国第二长河，也是世界著名河流之一，起源于青海省巴颜喀拉山脉，干流全长约5 464 km，流经青海、宁夏、内蒙古、陕西、山西、河南和山东等省份，最终在山东省东营市垦利区流入渤海. 黄河流域面积约752 400 km²，支流繁多，含沙量很高. 黄河中下游地区的降雨量分布不均，容易引发洪涝和干旱灾害.

在黄河流域的众多研究站点中，花园口与利津具有特殊地位. 花园口站位于河南省郑州市，是黄河中下游的分界点，既承接了上游高原山地的来水来沙，又临近下游平原地区，是监测水文地形过渡变化的关键节点；而利津站位于山东省东营市，靠近黄河入海口，作为黄河入海前的最后一道“关卡”，其径流数据反映了整个黄河流域的水量及水沙输送的最终状态. 如图1所示，2个站点分别代表黄河中游与下游的典型水文地理特征，流域地形复杂多变，从高原、山地到平原，有着丰富的地貌特征，具有良好的代表性.

选择花园口与利津作为研究案例，既具备针对性又具备推广价值，所提模型可应用于其他具备类似水文条件或数据结构的站点，特别适用于非线性强的径流预测任务. 2个水文站对应的名称和经纬度坐标详见表1. 所选数据集涵盖2008.01.01~2022.12.30时间段的全年水文数据，还包括各个测站的地理位置坐标数据，共计5 478条数据. 该流域的水文数据从水利部黄河委员会^[33]获取，这些数据主要包括流量、含沙量、水位共3种水文特征，每天采集1次，水文数据的变量单位分别为m³/s、kg/m³和m，花园口的流量、含沙量和水文数据缺失个数均为16，利津的数据缺失个数均为12.

Spearman秩相关系数^[34]是一种非参数统计方法，用于衡量2个变量之间的单调关系. Spearman不要求变量之间存在线性关系，而是关注变量之间的单调性. 该方法通过将数据转换为秩（排名）来计算相关性，从而避免了异常值对结果的影响.

1) 排名：将2个变量的数据分别进行排序，并为每个数据点赋予秩. 如果数据中存在重复值，则为这些重复值分配相同的平均排名.

2) 计算秩差：对于每对数据点，计算它们在2个变量中的排名差值：

(1)$ d_i^2 = {\text{rank}}\,\left( {{X_i}} \right) - {\text{rank}}\,\left( {{Y_i}} \right){\text{ }}{\text{.}} $

3) 计算相关系数：

(2)$ \rho = 1 - \frac{{6\displaystyle\sum {d_i^2} }}{{n\left( {{n^2} - 1} \right)}}{\text{ }}{\text{.}} $

式中：$ d_i^{} $为每对数据的秩差，$ n $为数据点的数量. Spearman相关系数的取值范围为[−1,1]，其取值越接近1.0，表明变量之间的正相关性越强.

VMD是一种非递归信号分解方法，其核心是通过构建变分优化问题，将时序自适应分解为若干具有特定中心频率和带宽的固有模式函数（intrinsic mode function，IMF）：

(3)$ {u_k}(t) = {A_k}\left( t \right)\cos \left( {{\phi _k}(t)} \right). $

式中：$ {A_{{k}}}(t) $为信号的包络幅值，$ {\phi _{{k}}} $为瞬时相位.

该方法通过希尔伯特变换构造变分模型，并结合维纳滤波技术^[35]约束IMF带宽，其优化目标是最小化各IMF中心频率带宽之和，同时满足模态分量之和等于原始信号的约束. 在求解过程中，引入二次惩罚项及拉格朗日乘子，其数学模型可表示为

(4)$ \left. \begin{gathered} \min \left\{ {{u_k},{w_k}} \right\} = \\ {\text{ }} \sum\nolimits_k {||{\partial _t}\left[\left(\delta (t)+\frac{{\mathrm{j}}}{{{\text{π}}t}}\right) * {u_k}(t)\right]\exp\, \left( { - {\mathrm{j}}{w_k}t} \right)||_2^2} ,{\text{ }} \\ s.{\text{t}}.\;\sum\nolimits_k {{u_k} = Q.} {\text{ }} \\ \end{gathered} \right\} $

式中：$ {w_k} $为中心频率，$ \delta $为达拉克分布.

LSTM通过引入门控机制和记忆细胞改进RNN，有效捕捉序列中的长期依赖关系，其核心由输入门、遗忘门和输出门构成. 假设上一时刻的输出为${h_{t - 1}}$，当前$t$时刻的输入为${x_t}$，则遗忘门的输出$ {{\boldsymbol{f}}_t} $和输入门的输出$ {{\boldsymbol{i}}_t} $可表示为

(5)$ {{\boldsymbol{f}}_t} = \sigma \left( {{{\boldsymbol{W}}_f} \cdot [{h_{t - 1}},{x_t}]+{{\boldsymbol{b}}_f}} \right), $

(6)$ {{\boldsymbol{i}}_t} = \sigma \left( {{{\boldsymbol{W}}_i} \cdot [{h_{t - 1}},{x_t}]+{{\boldsymbol{b}}_i}} \right). $

式中：$ {{\boldsymbol{W}}_{\boldsymbol{f}}} $和$ {{\boldsymbol{b}}_{\boldsymbol{f}}} $分别为遗忘门的权重与偏置，$ {{\boldsymbol{W}}_{\boldsymbol{i}}} $和$ {{\boldsymbol{b}}_{\boldsymbol{i}}} $分别为输入门的权重与偏置，$ \sigma $为Sigmoid激活函数.

根据上个时间点的输出$ {h_{t - 1}} $和当前时间的输入${x_t}$来计算当前记忆单元：

(7)$ \tilde {\boldsymbol{C}}_t= \tanh \left( {{{\boldsymbol{W}}_{\boldsymbol{C}}} \cdot [{h_{t - 1}},{x_t}]+{{\boldsymbol{b}}_{\boldsymbol{C}}}} \right). $

式中：$ {{\boldsymbol{W}}_{\boldsymbol{C}}} $和$ {{\boldsymbol{b}}_{\boldsymbol{C}}} $分别为记忆单元的权重与偏置，$ \tanh $为双曲正切函数. 之后更新记忆单元$ {{\boldsymbol{C}}_t} $：

(8)$ {{\boldsymbol{C}}_t} = {{\boldsymbol{f}}_t} \cdot {{\boldsymbol{C}}_{t - 1}}+{{\boldsymbol{i}}_t} \cdot \tilde {\boldsymbol{C}}_t . $

式中：$ {{\boldsymbol{C}}_{t - 1}} $为上一时刻的记忆细胞. 与遗忘门和输入门的计算方式相同，输出门$ {{\boldsymbol{o}}_t} $为

(9)$ {{\boldsymbol{o}}_t} = \sigma \left( {{{\boldsymbol{W}}_{\boldsymbol{o}}} \cdot [{h_{t - 1}},{x_t}]+{{\boldsymbol{b}}_{\boldsymbol{o}}}} \right). $

式中：$ {{\boldsymbol{W}}_{\boldsymbol{o}}} $和$ {{\boldsymbol{b}}_{\boldsymbol{o}}} $分别为输出门的权重与偏置. 当前输出${{\boldsymbol{h}}_t}$由记忆单元$ {{\boldsymbol{C}}_t} $和$ {{\boldsymbol{o}}_t} $共同决定：

(10)$ {{\boldsymbol{h}}_t} = {{\boldsymbol{o}}_t} \cdot \tanh \,\left( {{{\boldsymbol{C}}_t}} \right). $

Sowmya等^[36]提出基于牛顿-拉夫逊方法的新型优化算法NRBO. 该算法将种群分为探索者和开发者2种角色. 探索者负责在搜索空间中寻找新的潜在解，开发者利用这些发现来完善搜索过程. 探索者更新机制为

(11)$ {x_{n+1}} = {x_n} - \frac{{f'\left( {{x_n}} \right)}}{{f''\left( {{x_n}} \right)}}{\text{ }}. $

式中：$ f'\left( {{x_n}} \right) $为目标函数在当前位置的一阶导数，$ f''\left( {{x_n}} \right) $为二阶导数，$ {x_n} $为当前位置. 开发者更新规则为

(12)$ {x_{n+1}} = {x_n} - {\text{NRSR }}. $

式中：NRSR为牛顿-拉弗森搜索规则，用于指导开发者向更好的解迁移. 陷阱避免算子（trap avoidance operator，TAO）更新规则为

(13)$ \begin{split} {x_{n+1}} =& {x_n}+{\theta _1} \cdot \left( {{u_1} \cdot {x_{{\mathrm{best}}}} - {u_2} \cdot {x_n}} \right)+ \\ & {\theta _2} \cdot \partial \cdot \left( {{u_1} \cdot {\text{mean}}\left( {\boldsymbol{x}} \right) - {u_2} \cdot {x_n}} \right){\text{ }}. \end{split} $

式中：$ {x_{{\mathrm{best}}}} $为全局最佳位置；$ {x_n} $为当前位置；$ {\theta _2} $为服从正态分布的随机数；$ {u_1} $和$ {u_2} $分别为随机分配的1和−1的矩阵元素；$ \partial $为固定的参数，用于控制探索和开发的平衡.

采用相关性分析方法筛选模型的输入特征，并对历史径流序列进行VMD，将径流数据分解为多个IMF分量，捕捉径流在不同变化趋势下的信号特征，去除高频噪音，保留主要趋势. 将经过筛选的其他输入特征与IMF分量进行融合，构成多维特征矩阵. 为了进一步优化模型的性能，采用NRBO对LSTM的关键参数进行寻优，增强模型对径流变化规律的学习能力. 将各分量分别输入到优化的LSTM模型进行预测，并将各个预测结果进行叠加，得到最终的预测值.

1) 数据的缺失值和异常值处理：径流数据的完整性会影响预测的精度，为了提高数据质量，采用Python中的KNNImputer方法填充缺失值. 由于缺失比例极小，填充后的数据在统计特征与时序结构方面仍满足建模要求. 部分异常值可能反映极端气候或特殊水文事件，具有重要信息价值，予以保留，以提升模型对复杂特征的识别能力.

2) 数据的相关性分析：为了深入分析径流量的影响因素，采用Spearman方法对花园口和利津2个站点的水文数据进行相关性分析. 假定输入特征序列为$ {{\boldsymbol{M}}_t} $，模型预测得到的理想输出序列为$ {{\boldsymbol{y}}_{t+1}} $，为了预测$ t+1 $时刻的输出，模型应合理运用该时间段前的数据序列：

(14)$ {{\boldsymbol{M}}_t} = \left[ {\begin{array}{*{20}{c}} \begin{gathered} {x}_1^1 \\ {x}_1^5 \\ \end{gathered} &\begin{gathered} {x}_2^1 \\ {x}_2^5 \\ \end{gathered} &\begin{gathered} {x}_t^1 \\ {x}_t^5 \\ \end{gathered} \end{array}} \right]{\text{,}} $

(15)$ {\overset\frown{{{\boldsymbol{y}}}}_{{{t+1}}}} = f\left( {{{\boldsymbol{M}}_t}} \right){\text{ }}{\text{.}} $

式中：$ {{x}}_i^j $为第$ j $个特征在第$ i $时刻的值. 径流预测的目标就是建立数学函数使得函数拟合值$ \overset\frown{\boldsymbol{y}}_{{{t+1}}} $和真实值$ {{\boldsymbol{y}}_{t+1}} $最接近，一般使用损失函数来判定.

3) 数据集划分：数据经过预处理后，按时间序列顺序将数据集的70%划分为训练集，数据集的15%划分为验证集，数据集剩下的15%划分为测试集，避免随机划分导致时序信息泄露，如图2所示，图中Q表示流量，t表示时间.

4) 归一化处理：采用最小−最大归一化方法将数据缩放到[0,1]区间，消除特征间量纲差异，提升模型稳定性与收敛效率.

(16)$ {y^*} = \frac{{y - {y_{\min }}}}{{{y_{\max }} - {y_{\min }}}}. $

式中：$ y $为原始数据，$ {y^*} $为归一化之后的数据，$ {y_{\max }} $和$ {y_{\min }} $分别为训练集数据的最大值和最小值. 为了避免数据泄露，严格采用训练集计算得到的$ {y_{\max }} $和$ {y_{\min }} $对训练集、验证集及测试集进行归一化处理，确保模型评估的公正性与可靠性.

为了充分挖掘多特征数据对径流趋势的影响，提出基于VMD的MFF方法，具体流程如下.

1) VMD的镜像延拓：在VMD算法执行过程中，信号边界处易产生模态混叠与失真问题. 为了有效缓解该问题，采用镜像延拓策略对输入信号进行预处理. 具体步骤如下：针对长度为$ T $的原始信号，将其前半段与后半段分别对称复制至信号两端，构建总长度为2$ T $的周期性延拓信号；将该延拓信号输入VMD算法进行分解，通过频域变换与迭代优化获取多个模态分量. 由于延拓部分仅为消除边界效应的辅助数据，实际有效信息仍集中于原始信号区域，在分解完成后，选取延拓信号的中心区间$ \left[ {{T \mathord{\left/ {\vphantom {T 4}} \right. } 4}+1,{{3T} \mathord{\left/ {\vphantom {{3T} 4}} \right. } 4}} \right] $进行截取，将该部分的模态分量作为最终有效输出. 该方法通过对称扩展与精准截取，既保证了信号处理的完整性，又最大限度降低了边界效应对分解结果的干扰，确保了后续径流趋势分析的准确性.

2) VMD分解：原始径流序列$ Q\left( t \right) $通过VMD分解为$ K $个IMF分量：

(17)$ [{\text{IM}}{{\text{F}}_1},{\text{IM}}{{\text{F}}_2}, \cdots ,{\text{IM}}{{\text{F}}_K}] = {\text{VMD}}\left( {Q\left( t \right),\alpha ,\tau ,K} \right). $

式中：$ \alpha $为带宽约束，$ \tau $为噪声容忍度，$ K $为分解个数. 高频分量反映短期波动，低频分量反映长期趋势.

3) 滑动窗口构建时序特征矩阵与直接多步预测：为了实现1、2、3 d的预见期的直接多步预测，采用滑动窗口技术对历史步长$ k $内的多特征横向拼接，对于每个IMF分量d，输入特征矩阵$ {\boldsymbol{X}}_t^{(d)} $为

(18)$ {\boldsymbol{X}}_t^{(d)} = \left[ {{H_{t - 1}}, \cdots ,{H_{t - k}}} \right] \oplus \left[ {{\text{IM}}{{\text{F}}_d}(t - 1), \cdots ,{\text{IM}}{{\text{F}}_d}(t - k)} \right]. $

式中：$ H $为水位，$ {\text{IM}}{{\text{F}}_d} $为第$ d $个VMD分量，$ \oplus $为向量拼接操作.

4) 独立建模与结果融合：每个IMF分量的特征矩阵$ {\boldsymbol{X}}_t^{(d)} $输入独立的LSTM预测子结果，最终预测值为各分量输出之和：

(19)$ Q\left( {t+1} \right) = \sum\limits_{d = 1}^K {{Q^{\left( d \right)}}\left( {t+1} \right)} . $

上述过程将水文数据和IMF分量融合在同一特征空间中，构建多维信息特征矩阵. 该矩阵作为模型的输入，不仅增强了模型理解时序数据复杂结构的能力，还能够为后续预测提供更完整的信息. 通过融合多特征数据，模型能够更准确地捕捉在不同的预见期内径流量与水文因子之间的复杂关系，从而显著提升预测精度和泛化能力. 此外，利用多维信息可以增强模型的可解释性，使预测过程更加透明，有助于更好地理解模型决策所依据的物理规律与数据特征.

构建深度学习模型，该模型由LSTM、ReLU激活函数、全连接层和回归层组成. LSTM层被用于捕捉径流时序数据的长期依赖关系；ReLU激活函数提供非线性变换以增强模型的表达能力；全连接层通过线性组合提取高阶特征；回归层则将模型输出映射到预测结果. 如图3所示.

在MFF的基础上，将NRBO应用于LSTM的关键超参数（隐藏层单元数、学习率和训练epoch数）优化. NRBO具备强大的全局搜索能力，能够有效地克服传统优化算法易陷入局部最优的问题，从而提升模型的超参数寻优效率. 与此同时，LSTM以优异的时序建模能力，能够充分捕捉径流时间序列的长期依赖性和复杂非线性特征. 两者结合可以提高径流预测模型的泛化能力和预测精度，实现对流域径流变化的更准确模拟和预估.

具体优化过程如图4所示. 将每个粒子编码为1组LSTM超参数组合，通过计算该组合下LSTM在验证集上的均方根误差（root mean square error，RMSE）作为适应度函数，算法以预设的最大迭代次数MaxIt作为停止条件，在每次迭代中，NRBO根据适应度值更新粒子位置，动态调整超参数组合，并记录当前全局最优适应度及其对应的超参数. 当迭代次数达到MaxIt时，算法终止，输出最优超参数组合用于LSTM训练. 通过该策略，既能保证算法在可控计算资源下完成寻优，又可避免因过早收敛导致的次优解问题. 此阶段中，NRBO在LSTM优化上采用多种方法.

传统的优化算法通常采用固定参数或线性衰减的方法来平衡搜索和开发，易导致模型陷入局部最优解或收敛速度过慢，影响泛化能力. 引入基于NRBO的非线性动态衰减权重因子$ \delta $：

(20)$ \delta = {\left(1 - \frac{{2 \cdot {\text{it}}}}{{{\mathrm{MaxIt}}}}\right)^5}. $

式中：it为当前迭代次数，MaxIt为最大迭代次数. $ \delta \in (0,1.0) $.

根据当前迭代次数与最大迭代次数的关系动态调整. 权重因子大时算法侧重全局搜索，权重因子小时转向局部微调. 通过自适应调节搜索范围，平衡模型搜索与开发权重，满足训练不同阶段的需求，实现LSTM模型最佳参数寻优，提升收敛速度与泛化能力.

传统的LSTM模型训练时多使用固定学习率，无法适应不同训练阶段的需求，导致训练效率低下. 通过引入基于NRSR的自适应学习率动态规划方法，计算当前位置与最佳位置的差异，并结合随机数生成动态步长，使LSTM在不同训练阶段自动调整学习率，有效提升训练效率.

传统优化算法在迭代过程中容易陷入局部最优解，缺乏有效的机制来帮助粒子逃离局部最优解. 为此，引入基于NRBO的双随机扰动机制，该机制通过在方向和强度2个方面引入随机性，并结合$ \delta $，在迭代早期提供较大扰动，帮助粒子逃离局部最优解；在迭代后期逐渐减少扰动，以增强算法的开发能力. 通过设置触发机制，让算法仅在部分迭代中激活TAO，避免了每次迭代都进行复杂的随机扰动计算，从而平衡了计算开销和算法性能. 这种概率触发机制在保证算法效率的同时，可以有效地避免过早收敛.

在高维搜索空间中，传统的粒子位置更新方法依赖逐元素操作，计算效率低，在NRBO算法中引入向量化操作和矩阵运算，将多个粒子更新合并为单个矩阵运算，减少了循环和条件判断开销. 在LSTM应用中，这种高效的计算方式可以缩短训练时间.

选取4种常用的评价指标：RMSE、平均绝对百分比误差（mean absolute percentage error，MAPE）、平均绝对误差（mean absolute error，MAE）和纳什效率系数（Nash-Sutcliffe efficiency coefficient，NSE）.

(21)$ {\text{RMSE}} = \frac{1}{n}\sqrt {\sum\limits_{i = 1}^n {{{({y_i} - {{\hat y}_i})}^2}} {\text{ }}} , $

(22)$ {\text{MAPE}} = \frac{1}{n}\sum\limits_{i = 1}^n {\frac{{|{y_i} - {{\hat y}_i}|}}{{|{y_i}|}}} {\text{ ,}} $

(23)$ {\text{MAE}} = \frac{1}{n}\sum\limits_{i = 1}^n {\left| {{y_i} - {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{y} }_i}} \right|} {\text{ ,}} $

(24)$ {\text{NSE}} = 1 - {{{{\displaystyle\sum\limits_{i = 1}^n {{{({y_i} - {{\hat y}_i})}^2}} }}} \mathord{\left/ {\vphantom {{{{\displaystyle\sum\limits_{i = 1}^n {{{({y_i} - {{\hat y}_i})}^2}} }}} {{{\displaystyle\sum\limits_{i = 1}^n {{{({y_i} - {{\bar y}_i})}^2}} }}}}} \right. \kern-\nulldelimiterspace} {{{\displaystyle\sum\limits_{i = 1}^n {{{({y_i} - {{\bar y}_i})}^2}} }}}} . $

式中：$ {y_i} $为第$ i $天径流实测值，$ {\hat y_i} $为第$ i $天径流预测值，$ {\bar y_i} $为实测径流平均值，$ n $为径流时间序列的长度.

基于NRBO优化LSTM的损失函数采用均方误差（mean square error，MSE）：

(25)$ {\text{MSE}} = \frac{1}{n}\sum\limits_{i = 1}^n {{{({y_i} - {{\hat y}_i})}^2}} . $

在LSTM的优化过程中，所使用的目标函数$ \mathcal{J}({\boldsymbol{\theta }}) $不仅包含基础损失函数项$ L({\boldsymbol{\theta }}) $，还考虑了模型复杂度的约束（例如引入正则项防止过拟合）：

(26)$ \mathcal{J}({\boldsymbol{\theta }}) = L({\boldsymbol{\theta }})+\lambda \cdot {\left\| {\boldsymbol{\theta }}\right\|^2}. $

式中：$ L({\boldsymbol{\theta }}) $为MSE；$ {\left\|{\boldsymbol{ \theta }}\right\|^2} $为L2正则项；$ \lambda $为正则化系数，设置为$ 1 \times {10^{ - 7}} $.

模型训练环境如下：处理器为Inter Core i9-10900K CPU和NVIDIA GeForce RTX 1660 Laptop GPU，内存为64 GB. 模型的初始参数如下：VMD的分解层数为5层；CNN模型的卷积核大小为2，输出通道数为16；TCN模型包含2个一维卷积层，每个卷积层使用大小为2的卷积核，输出通道数为16. CNN、TCN和LSTM的隐藏层层数为70，初始学习率为0.01，最大训练轮数为70，L2正则化参数为1×10⁻⁷，梯度阈值设置为1，学习率调整策略为分段常数，学习率降低周期为60. 优化算法的种群数量为100，最大迭代次数为6，变量为最佳隐藏单元数目、最佳训练周期和最佳初始学习率，参数的上限分别为500、500、0.01，下限分别为15、50、0.000 1.

通过图5的相关性分析，可以发现，2个站点的水位和流量的相关性最高，而其他特征值和流量的相关系数不足0.5. 因此，将水位和流量作为后续研究的输入特征.

为了全面评估模型的性能，将MFF-NRBO-LSTM与随机森林（random forest，RF）^[37]、卷积神经网络（CNN）^[38-39]、TCN、TCN-LSTM、LSTM、MFF-LSTM、MFF-SSA-LSTM和MFF-黑翅鸢算法（black kite algorithm，BKA）^[40]-LSTM这8种模型进行对比. 当预见期为1 d时，结果如表2所示，MFF-NRBO-LSTM模型在所有评价指标上均表现最优，说明该模型的预测精度最高. 与单一的LSTM模型相比，MFF-NRBO-LSTM模型的RMSE、MAE和MAPE分别下降了74.68%、74.59%、70.98%，NSE提升了5.56%.

MFF-LSTM模型在引入MFF后，增强了模型的时域分解能力，能够有效地提取径流序列中复杂的时序特征，精准捕捉径流变化趋势，提升模型对径流序列动态变化的学习与预测能力. 与原始的LSTM相比，MFF-LSTM的RMSE、MAE和MAPE分别下降了45.70%、44.68%、37.51%，而NSE提升了4.09%. 验证了MFF方法在特征提取和噪声抑制方面的有效性. 进一步地，MFF-NRBO-LSTM的RMSE、MAE和MAPE相较于MFF-SSA-LSTM和MFF-BKA-LSTM分别降低了29.12%、17.37%、4.63%，说明NRBO算法在径流预测方面优于SSA和BKA算法. 如表2和3所示，当预见期为1、2、3 d时，MFF-NRBO-LSTM模型的精度均为最高.

在5种单一模型中，LSTM的拟合效果最好，在后面几种混合模型中，MFF-NRBO-LSTM对于变化趋势和峰值的拟合效果均最好. 图6和7的散点图的v表示实际值，p表示预测值. 2个站点的实际值和预测值的散点图也进一步表明，在单一模型中的LSTM和在混合模型中的MFF-NRBO-LSTM的预测点最为密集地分布在拟合线附近，说明其具有更强的预测一致性.

消融实验从单独引入NRBO和MFF这2个关键组件出发，进一步验证各部分对模型性能的贡献. 如表4所示，NRBO-LSTM相比LSTM的RMSE、MAE和MAPE分别下降了9.81%、6.37%和3.21%，NSE则有所提高，表明NRBO-LSTM具备较强的参数优化能力. 引入MFF后，MFF-NRBO-LSTM比NRBO-LSTM的RMSE、MAE和MAPE分别下降了71.91%、72.81%和70.12%，NSE提升了4.04%，说明MFF能够聚焦多特征数据关键信息，并能够显著提高预测精度.

为了克服现有预测模型在提取多特征数据特征方面的局限性，并提升日径流预测的精度，创新性地提出MFF-NRBO-LSTM多步预测模型，该模型集成了MFF和NRBO优化2大阶段. 在MFF阶段，通过MFF方法丰富模型的信息来源，使MFF-NRBO-LSTM能够更全面地理解数据的结构，从而显著增强了对序列特征的捕捉能力. 同时，NRBO对LSTM的超参数进行了精细调整，不仅增强了模型的全局搜索能力，也显著提升了预测性能. 为了验证该模型在日径流多步预测中的适用性，对黄河流域的花园口和利津水文站这2个具有代表性的站点进行预见期为1~3 d的多步预测，实验结果显示，MFF-NRBO-LSTM与其他8种对比模型相比，预测精度最高，误差最小. 有效验证了该模型的实际应用价值.

尽管本研究重点突出了NRBO在提升LSTM预测精度方面的显著优势，但未来的研究可以进一步探索该模型在收敛速度和计算复杂度方面的表现. 通过全面评估其在水文预测任务中的适用性和效率，可以为模型的选择和优化提供更为科学的依据，从而推动水文预测技术的进步和发展.