基于LSTM实时校正的WRF/WRF-Hydro耦合径流预报

刘昱辰，刘 佳，刘录三，李传哲，王 瑜

(1.中国环境科学研究院，北京 100012；2.中国水利水电科学研究院 流域水循环模拟与调控国家重点实验室，北京 100038)

1 研究背景

在气候变化的大背景下，极端降雨过程会更加频繁和强烈地发生，气象预报与水文预报的结合，有助于探索大气过程对陆面水文过程的影响，提高降雨和径流的模拟、预报能力，是水文预报发展的必然趋势[1]。相比其他陆气耦合系统，中尺度数值大气模式WRF(Weather Research and Forecasting Model)与其陆面水文模式WRF-Hydro的耦合可以实现大气与陆面水文过程的单向与双向交互，模拟水热通量在大气-陆面界面的运动与转化过程，其应用前景十分广泛[2-4]。但受到模型结构体系、初始状态和输入信息误差、参数辨识能力等限制，该耦合系统的不确定性较为显著且难以量化[5-6]，模拟结果难以直接应用于实际预报，对洪水预警的指导意义有限。由于降雨和径流形成的物理机制较为复杂，现有各类陆气耦合预报系统在模拟和预报精度方面往往都会存在一定的偏差。研究表明[7-8]，对陆气耦合系统的径流预报结果进行实时校正，能够有效提高耦合系统的降雨-径流预报能力。

随着计算机技术的发展，人工智能与大数据技术发展迅速，其中，作为人工智能核心技术之一的人工神经网络，广泛应用于水文预报的实时校正领域[9]。带有记忆单元的长短时记忆神经网络(Long Short Term Memory network，LSTM)[10]的提出，使时间序列上的记忆信息变得可控，也使LSTM在模拟和预报具有时间序列属性的过程方面发展的愈发成熟。Cho等[11]研究表明，将LSTM结合WRF-Hydro可以改进径流预报结果，有效提高WRF-Hydro在径流预报方面的潜力。Wang等[12]建立了基于LSTM的网格地表产流模型，结果证明LSTM在我国北方山区小流域地表径流预报中可表现出良好的性能。在LSTM应用于水文预报领域的研究中，大多是采用LSTM直接对径流等特征量进行预测[13-15]，缺乏与具有物理机制的水文模型相结合。在传统误差实时校正模型中，自回归滑动平均模型(Autoregressive Moving Average model，ARMA)是一种较为常用的函数模型，相对人工智能误差校正方法，ARMA模型参数简单，不易受环境变量的影响。这种方法在水文预报领域已被普遍认可，得到了广泛的应用研究[16-18]。

为充分提高耦合系统的径流预报能力，本研究在使用WRF模式数据同化改进预报降雨的同时，针对耦合系统仍然存在的径流预报偏差，采用LSTM和ARMA模型分别对预报径流进行实时校正。首先，通过同化雷达反射率与传统气象监测数据GTS(Global Telecommunication System)，改进WRF模式的降雨预报精度，在此基础上构建出WRF/WRF-Hydro耦合系统；其次，配合LSTM模型，实现预报径流的实时校正，提高陆气耦合系统的综合预报能力；最后，与传统实时校正模型ARMA结果进行对比验证。研究成果可为中小尺度流域的洪水风险管理与预报预警提供技术支撑。

2 研究区域概况

本研究选取位于大清河流域南支的阜平流域和北支的紫荆关流域作为研究区，流域面积分别为2210 km2和1760 km2，研究区地处太行山东麓，全境为山地，地势高低落差明显，陡峭的地形、裸露的地表和复杂的产汇流机制，极易造成严重的洪水灾害。流域面积研究区地理位置、地形地势以及雨量站的分布情况见图1。

图1 研究区位置图

根据阜平流域和紫荆关流域1990年以来发生的历史洪水情况，依据同流域相关研究[12，17，19-20]，按照不同的降雨类型及径流特征，选取了6场典型降雨过程开展陆气耦合径流预报研究。选取的降雨场次在气象成因、强度、历时、落区上均具有一定的代表性，其分别具有雨强小且持续时间较长、雨强较小但引发了明显的流量涨落、极短时间内雨强大引发洪水陡涨陡落等多种不同特点。为方便分析和比较，所选取6场降雨的降雨时长均为24 h，其降雨历时等详细信息见表1。其中，场次4在统计时段降雨发生前，已发生过较多降雨场次，自2016年5月1日至7月24日，阜平流域发生的累积面雨量达340.5 mm，并且在7月19日至7月22日，出现过一个降雨强度大且持续时间长的洪水过程，其中洪峰流量达到了581 m3/s。因此，在场次4发生时，流域内存在较大基流，且土壤湿度较为饱和，极易引发洪水，7月25日的最高洪峰流量达到1740.0 m3/s[21]。

表1 研究流域6场典型降雨场次信息

降雨场次的降雨强度、落区、历时的不同，会产生不同的径流形成结果。对流域而言，降雨时空分布的均匀程度是影响流域出口断面流量过程的最重要因素。为定量的描述降雨过程的时空分布情况，采用变差系数Cv[22]，分别从空间和时间尺度对典型降雨过程的均匀程度进行评价，并将其作为划分降雨类型的指标：

(1)

通过设置Cv临界值划分降雨类型，降雨空间分布均匀性的临界值Cv为0.40，降雨时间分布均匀性的临界值Cv为1.00，当降雨场次的Cv值小于Cv临界值时，判断降雨在时间或空间上是相对均匀的[17，21]。各降雨场次空间分布和时间分布的Cv值及降雨类型见表2。

表2 6场降雨的Cv值及降雨类型

3 实验方案

研究区域的地面观测数据包括来自阜平流域8个、紫荆关流域11个雨量站的实测小时降雨量和阜平水文站、紫荆关水文站的实测小时径流量，空间分布见图1，均由河北省水文水资源勘测局提供。本研究针对雷达反射率和GTS数据进行同化。雷达数据采用了位于石家庄市的SA波段多普勒雷达系统，该系统的初始数据来自中国气象局。雷达扫描范围见图2，覆盖半径为250 km，空间分辨率为1 km。GTS是传统的地面与高空气象数据的集合，由世界气象组织(WMO)提供，本研究通过程序仅截取覆盖研究区的GTS数据作为同化资料，并只针对该区域的数据转化为可被WRF-3DVAR同化的格式。GFS(Global Forecasting System)是由国家环境预报中心提供的实时数据，用于驱动WRF模式的初始场和边界场，空间分辨率为1.0°×1.0°，每隔6 h更新一次。

图2 阜平流域和紫荆关流域的网格嵌套方案

3.1 WRF/WRF-Hydro耦合预报系统构建考虑到两个研究流域位于山区，降雨的不确定性更为突出，WRF模式预报结果仍会受到空间误差的影响[23]。为了获取可靠的预报降雨，前期利用天气雷达对云雨系统的探测能力和获取瞬时降雨信息的优势，结合传统气象监测数据，基于WRF-3DVAR三维变分数据同化系统，开展了数据同化频率对WRF模式降雨预报的影响研究。前期成果[19]为本文研究提供了可靠的降雨输入。

为使设置的嵌套网格与雷达数据的空间分辨率尽可能相近，并提高优化效率，采用三层嵌套网格的方式逐层增加区域的水平分辨率，最外层网格范围的设置利于同化大尺度的GTS数据，使模式运行更加稳定。最内层网格完全覆盖研究区，相邻嵌套层的嵌套比例设为1∶3，网格尺寸分别为1、3和9 km。1.0°×1.0°全球预报系统为预报提供了初始和横向边界条件。由于GFS没有经过同化分析处理，因此在WRF模式中常作为预报降雨的初始驱动[24]，流域的嵌套网格和雷达观测范围见图2。

参数化方案的选择借鉴了Tian等[22]在相同研究流域的参数化方案敏感性试验结果。对同化试验中使用的降雨产生有重大影响的参数化细节见表3。其中，1和3 km网格的积云对流方案被关闭，WRF-Hydro模式仅在WRF最内层嵌套网格(即1 km网格)上运行，并采用100 m分辨率的次网格进行水文过程的模拟。陆面模式方案选择Noah，驱动数据输入步长和输出数据的时间间隔均为1 h。研究发现，在两个研究流域，数据同化对WRF模式预报降雨量及其时空分布的改进明显：数据同化前，WRF模式的预报降雨量远低于实际观测，数据同化后，各小时的预报降雨量均有所增加，同化雷达反射率结合GTS传统气象监测数据对WRF预报降雨有明显的改进效果，且1 h同化时间间隔和1 km分辨率的降雨输出可以使WRF模式更好地预报降雨的时间和空间分布情况。

表3 WRF模式主要物理参数化方案

因此，在本研究中，WRF-Hydro的输入由逐小时循环同化后的WRF预报降雨提供，基于数据同化的陆气耦合流程见图3。

图3 基于数据同化的WRF/WRF-Hydro陆气耦合预报系统示意

3.2 径流预报偏差校正方法

3.2.1 LSTM实时校正模型 LSTM是一种特殊的循环神经网络，其允许先前输入的信息保存在网络的内部，从而影响网络的输出，相比于传统的循环神经网络，LSTM在隐藏层中加入了记忆细胞单元，可以较好地记住较长的历史信息。

记忆细胞单元可以对输入的数据进行选择性的记忆和遗忘，通过引入遗忘门、输入门、输出门，使神经网络长期传输记忆信息，LSTM结构见图4。在本研究中，LSTM模型建模在Matlab下运行，研究所使用的LSTM及其变体模块，均采用Matlab(版本R2021)框架中的LSTM模块设计。根据文献[26-28]，当流域较小或可获取的数据量不充足时，采用单层隐含层可获得更好的训练效率及运行结果，因而本研究选择了1层作为LSTM的隐含层数。通过试算得出，隐藏神经元个数n=32时可以使模型均方误差最小。模型训练前，还需要确定模型的学习率、批处理量参数和选代次数等参数。其中，学习率决定了模型权重更新的幅度，但目前还缺少确定最佳学习率的最优方法。在本研究依据经验将学习率设定为0.01，同时将学习率衰减因子设定为0.99，使学习率随着迭代次数增加缓慢变小，则模型收缩时震荡能得到缓解。批处理量参数即一次输入模型中的样本个数，参照相关研究[29-31]，批处理样本数与隐含层神经元数相同时处理效果较好，故将批处理量定为32。epoch是参数选代更新次数，根据本研究中的数据量大小，将epoch定为200。LSTM主要参数设置见表4。在本研究中，基于WRF-Hydro的逐小时径流预报构建了逐小时预报径流误差的历史序列，并将其作为LSTM输入变量，以1 h为时间步长预报不同预见期下的径流误差，最后将预测的误差序列叠加至径流预报结果中，完成耦合系统水文预报的实时校正。

表4 LSTM设置参数表

图4 LSTM模型结构示意图[25]

3.2.2 ARMA实时校正模型 ARMA模型由自回归模型(AR模型)与移动平均模型(MA模型)为基础“混合”构成，是可以反映水文变量在时序变化上统计特性的数学模型。本研究针对研究流域内实测径流量和WRF-Hydro预报径流量之间的逐小时误差序列{dQ(x)，x=1，2，…，t}，通过构建具有时间特征的误差序列，建立基于误差的自回归方程。设t时刻的预报值为Q′(t)，实测值为Q(t)，则dQ(x)=Q′(t)-Q(t)即为误差序列。假定误差序列具有一定的自相关性，则可以利用过去的误差预报未来的误差。ARMA(p，q)的结构取决于阶数p和q。在本研究流域，已有相关研究表明当p=3，q=1时可以取得较好的校正效果[17]，因此，本研究选用ARMA(3，1)模型对耦合系统的预报径流过程进行实时校正。

3.3 评价指标利用纳什效率系数(NSE)、均方根误差(RMSE)、洪峰流量的相对误差(Rf)、洪量的相对误差(Rv)以及峰现时间绝对误差(T)对WRF-Hydro模拟的径流量进行评估。对于评价指标，NSE用以验证水文模式模拟结果的可信程度。RMSE用来衡量预报值同观测值之间的偏差。Rf为洪峰的相对误差，反映模拟洪峰的可信程度，Rv为洪量的相对误差，反映模拟洪量的可信程度。T为实际峰现时间和预报峰现时间的时差，用于表示峰现时间预报的误差。

(2)

(3)

Rv=(r′v-rv)/rv

(4)

Rf=(r′f-rf)/rf

(5)

T=|T1-T2|

(6)

4 结果分析

4.1 校正前的WRF/WRF-Hydro耦合预报结果同化后预报降雨驱动下的紫荆关流域和阜平流域的径流预报结果如图5所示。为更好的显示降雨引发的流量起涨过程，在24 h雨洪过程的基础上，对实测径流时间进行了延长。

图5 同化后预报降雨驱动下的WRF/WRF-Hydro耦合系统径流预报结果

总体上，通过WRF-3DVAR同化系统，可以使WRF模式预报降雨量增大，降雨预报整体与实际降雨量更为接近，其预报径流的起涨过程与实际流量过程线也较为吻合。同时，WRF-Hydro模式对量级小、陡涨陡落、退水较快的洪水过程有更好的模拟效果，但对量级较大且退水较慢的洪水过程流量过程线吻合程度较差。表5给出了WRF-Hydro径流预报结果。

表5 同化前后预报降雨动下的WRF-Hydro径流预报结果

从评价指标结果可以看出，WRF-Hydro径流预报精度主要取决于WRF模式的降雨预报误差，与驱动数据质量和降雨分布类型密切相关。如强对流降雨场次3，降雨在空间和时间尺度的变差系数Cv值分别达到了0.740和2.393，短时集中的强降雨过程导致其预报结果存在较大误差，不理想的降雨预报间接影响了洪水的预报精度。虽然数据同化技术为水文过程的模拟与预报提供了更精确的降雨输入，但径流与降雨并非简单的线性关系，模型结构的复杂程度和参数率定过程都会对径流预报产生影响，多因素导致WRF-Hydro在径流预报结果上仍存在一定的不确定性。

4.2 LSTM实时校正径流预报结果采用LSTM模型配合WRF/WRF-Hydro耦合预报系统，对6场径流的预报结果进行实时校正，预见期为1、2、3、6、9和12 h的径流预报结果见图6和表6。

表6 LSTM校正下WRF/WRF-Hydro耦合系统的径流预报指标

图6 LSTM校正后的WRF/WRF-Hydro径流预报结果

由图6可以看到，LSTM模型对WRF/WRF-Hydro耦合预报系统径流预报的实时校正具有较好的改进效果，尤其是对耦合系统在1 h预见期的预报结果，改进最为明显。预见期超过6 h后，LSTM模型对耦合系统径流预报的改进程度总体呈下降趋势，预报径流过程的波动增大，但洪峰流量和洪量没有明显的增加或降低。分析其原因，可能是由于随着预见期的延长，预报径流误差明显增加，径流序列数据之间的关联性降低。LSTM模型此时已较难准确学习到时间序列的数据特征，因而无法对预报径流序列进行有效的校正，具体体现在径流预报的波动也就越大。

将LSTM模型引入WRF/WRF-Hydro耦合预报系统之后(评价指标见表6)，场次洪水的预报精度有了明显的提高。整体来看，6场雨洪过程1 h预见期的NSE在0.75～0.98之间，平均NSE为0.90。随着预见期的延长，NSE逐渐降低。当预见期小于6 h时，预报精度下降较为缓慢，但仍有较好的预报校正效果；预见期大于6 h时，预报精度下降较快，其中表现最好的雨洪过程(场次1)在12 h预见期的NSE仍可达到0.88。从洪峰相对误差指标来看，1 h预见期的洪峰预报精度最高，径流过程的洪峰相对误差位于区间范围为[-3.65%，5.78%]，3 h预见期大多数场次洪峰相对误差小于20%。但由于预报流量过程线随着预见期的延长波动较大，洪峰流量和洪量在较长预见期内不呈现明显的规律。

4.3 ARMA实时校正径流预报结果以数据同化后的预报降雨驱动WRF-Hydro，采用ARMA模型对6场典型雨洪过程的径流预报结果进行实时校正，预见期为1、2、3、6、9和12 h的径流预报结果见图7和表7。

表7 ARMA校正下WRF/WRF-Hydro耦合系统的径流预报指标

图7 ARMA校正后的WRF/WRF-Hydro径流预报结果

经ARMA模型校正后，在1 h预见期，耦合系统能准确预报出所有径流过程的发展趋势，NSE在0.85～0.98之间。3 h预见期的预报径流比较稳定，平均NSE为0.64，但3 h预见期后的NSE开始出现较明显的下降。除预见期外，ARMA的实时校正效果与校正前的径流预报结果也有一定关系。例如，校正后场次1在3 h预见期的径流预报NSE为0.89，12 h预见期的径流预报NSE仍能达到0.56，实时校正结果十分理想。对比之下场次3的预报结果和校正效果则相对较差，预报精度衰减速度也相对较快。校正前，场次1的径流预报NSE为0.85，洪峰误差为5.15%，场次3的NSE为0.39，洪峰误差为-25.43%。由于校正前场次3的径流预报精度较低，ARMA作为误差序列的实时校正模型，是根据以往的历史径流误差和扰动来推测未来时刻的误差值，所以当误差序列相关性较差时，误差校正效果也不明显。

4.4 实时校正结果对比将LSTM与ARMA模型校正下的耦合系统径流预报结果进行深度对比，进一步分析LSTM对WRF/WRF-Hydro陆气耦合系统的径流预报改进效果。图8给出了校正前后6场典型雨洪过程不同预见期的NSE指标变化情况。

图8 LSTM和ARMA对6场典型雨洪过程实时校正后的NSE指标对比图

由图8可以看出，与校正前的径流预报结果相比，实时校正模型使预报径流的精度明显提高。根据不同预见期与1 h预见期的NSE指标变化百分比计算预报径流的NSE指标衰减速度，可以发现，LSTM校正下耦合系统预报径流的衰减速度明显慢于ARMA模型。例如，场次1在LSTM模型校正下，6 h预见期径流预报的衰减速度为5.10%，12 h预见期的衰减速度为10.20%，而ARMA模型校正下6 h预见期的衰减速度为18.37%，12 h预见期的衰减速度为42.86%。除场次3以外，其他场次的径流预报都显示了相似的衰减规律。这与校正前耦合系统的径流预报结果有关，误差序列之间的相关性越强，径流误差的校正效果则更为显著。

5 讨论与结论

5.1 讨论为了获得较高的水文预报精度，通常采用实时校正方法，针对历史预报误差的相关性对预报径流进行实时校正，但该方法受预见期的影响较大。就本研究流域而言，由于流域汇流时间较短，实时校正后的预报径流过程线出现了不同程度的波动(如图6、图7)，且波动随着预见期的延长而增大。可见，LSTM与ARMA模型对中小流域在6 h以上的预见期下预报径流的实时校正能力有限。考虑到大尺度流域汇水时间较长，流量起涨平稳缓慢，预报径流误差即使在6 h或12 h预见期以上也具有较高的相关性，LSTM与ARMA模型的校正潜力则可以体现在更长预见期的径流预报。

5.2 结论为进一步提高WRF/WRF-Hydro耦合系统的径流预报能力，在数据同化预报降雨的同时，针对耦合系统仍然存在的径流预报偏差，研究采用LSTM模型对预报径流进行实时校正，并与ARMA模型的径流校正结果进行了对比，结论如下：

(1)ARMA(3，1)与LSTM都能进一步改进WRF/WRF-Hydro耦合系统的径流预报能力，尤其在1 h预见期对预报径流的改进效果最好。校正后，研究区6场雨洪场次的1 h预见期NSE区间分别为0.75～0.98(LSTM)、0.85～0.98(ARMA)，但随着预见期的延长，两个模型对径流预报的改进程度逐渐降低。

(2)随着预见期的延长，LSTM模型对耦合系统径流预报的整体改进效果优于ARMA模型，耦合系统预报精度的衰减速度慢于ARMA模型。当预见期达到9 h时，LSTM在研究区大部分雨洪场次的NSE指标在0.50以上。

本研究采用LSTM构建了预报径流的误差校正模型，实现了对WRF/WRF-Hydro耦合系统预报径流的实时校正。研究还存在一些不足，如LSTM模型在预见期超出6 h时的预报精度仍有待提高等。在进一步的研究中，将采用LSTM与小波变换或BP神经网络等方法结合，利用径流序列的频率域误差特征，结合多种校正方法，开展多模型联合优化方法，以更好地提高陆气耦合预报系统的准确性。