融合气象要素时空特征的深度学习水文模型

李 步，田富强，李钰坤，倪广恒

(清华大学水沙科学与水利水电工程国家重点实验室，北京 100084)

20世纪50年代以来，在气候变化与人类活动的叠加影响下，洪涝灾害日益加剧，严重威胁人民群众的生命财产安全[1-2]。高精度的径流预报与模拟是有效应对洪涝灾害和提高水资源管理水平的重要手段[3]。概念性或有物理基础的水文模型是常用的水文预报与模拟模型[4]，主要包括新安江和HBV等集总式模型，以及SWAT、MIKE和THREW等分布式模型。但是由于部分水文过程机制认识不足，目前物理水文模型精度仍难以满足实际洪水应对与高质量水资源管理的需求[5-6]。

深度学习近年来发展快速，从最初的多层感知器(即人工神经网络)发展到结构更加高级的卷积神经网络[7]和循环神经网络(RNN)，这些模型广泛应用于水文模拟中[8-9]。然而，由于模型结构限制，这些模型均难以准确刻画长时序过程规律，因此在水文模拟中难以准确反映诸如融雪、地下水等长时间尺度水文过程[10-11]。长短时记忆神经网络(LSTM)模型具有独特的“门”结构[12]，很好地解决了长时序过程训练梯度爆炸和梯度消失的问题，提高了长时序过程模拟精度。Kratzert等[11]将LSTM模型应用在美国241个小流域的径流模拟中，模拟精度优于物理水文模型和RNN。后续大量研究将LSTM模型成功应用在世界不同流域，例如Lees等[13]应用于英国669个小流域，Liu等[14]应用于中国汉江流域等。然而，目前在基于LSTM模型的径流模拟研究中，模型输入数据均为流域空间的平均值，损失了流域气象数据的空间信息。以往的研究表明，流域气象空间异质性是影响径流过程的重要因素[15]，如何将流域气象要素的空间特征降维输入至LSTM模型中，以提高LSTM模型径流模拟精度，值得进一步研究。此外，黄河源区是黄河流域重要的产流区和水源涵养区，径流变化将直接影响黄河中下游的用水安全。黄河源区属于青藏高原高寒流域，水分多相态变化频繁，高寒下垫面产汇流机制尚未完全厘清，给黄河源区高精度径流模拟预报与径流变化研究带来了挑战[16]。

本文将利用主成分分析(PCA)和LSTM搭建一种综合考虑气象要素时空特性的深度学习水文模型(PCA-LSTM)，以黄河源区为研究区域，以LSTM模型和THREW模型为比对模型，系统评估PCA-LSTM模型在水文模拟中的适用性和鲁棒性。

1 模型原理和搭建

1.1 PCA-LSTM模型

本文提出一种耦合PCA和LSTM的水文模型，如图1所示，基于PCA提取气象要素的空间信息，利用LSTM学习长时序过程规律。模型的输入数据为流域降水和气温等网格气象数据，输出为径流过程。耦合模型内部结构与比对模型简介如下。

1.1.1 主成分分析法

PCA是一种重要的降维方法，通过正交变换将原始变量变换为一组不相关的变量，实现利用较少变量来代替原始的较多变量，变换后的较少变量称之为主成分。主成分变量一般基于累计方差贡献率确定，当方差的累计贡献率大于或等于一定阈值时(以下简称为主成分阈值)，可利用这些成分来刻画原始变量。PCA具体计算步骤可参考文献[17]。本文利用PCA提取流域降水和气温等空间气象数据的主成分，作为LSTM模型的输入。主成分阈值设置为85%、90%和95%，基于径流模拟性能获取最优主成分阈值。

1.1.2 长短时记忆神经网络

LSTM模型作为一种RNN的变式[12]，引入了特殊的控制单元，其内部结构如图1(b)所示。其中，细胞状态(ct)和隐藏状态(ht)分别捕捉缓慢和快速的演变过程，输入门(i)、遗忘门(f)和输出门(o)分别用于控制信息的存储、丢失和传递。这些特殊单元解决了RNN模型在处理长序列数据时存在的梯度爆炸或者消失的问题，使得LSTM模型擅于处理长序列过程。LSTM模型的详细运算过程与在水文模拟中的应用可参考相关文献[11]。

1.1.3 PCA-LSTM模型

基于以上算法，本文构建了耦合PCA和LSTM的水文模型，模型流程主要包括(图1(a))：

图1 PCA-LSTM模型结构与LSTM模型内部结构Fig.1 PCA-LSTM model architectures and the internals of LSTM cells

(1) 将降水、气温等气象空间数据分别重构为一维数据，并进行“最大值-最小值”归一化等预处理；

(2) 基于PCA和设定主成分阈值，分别筛选出降水和气温等气象数据的主成分变量，代表气象要素的空间特征；

(3) 多时刻(1～L)气象数据的主成分变量输入至LSTM模型，LSTM模型的输出变量引入全连接层(即线性回归模型)，获得待预测时刻(L)的径流变量。

1.2 比对模型

为了评估PCA-LSTM模型的适用性，本文搭建了2个比对模型：物理水文模型THREW和深度学习模型LSTM。其中LSTM模型的输入数据是流域空间平均的气象数据，其余设置与PCA-LSTM模型一致。

物理水文模型THREW采用代表性单元流域方法(REW)对流域进行空间离散，对每个REW划分为地上和地下2层，并基于下垫面和土壤特征划分为8种水文子流域类型，THREW模型具体描述可参考文献[18-20]。模型已成功应用于国内外多个流域，包括阿尔卑斯山以及青藏高原多个典型山区大流域[21-23]。THREW模型参数率定采用pySOT优化算法，该算法应用异步并行优化框架，充分利用计算资源，采用代理模型的方法寻找最优解，减少优化模型中的模型运行次数，提升计算效率。率定过程会在运行参数达到某一限定次数后停止，本文设置为3 000。

2 研究区域与数据

2.1 研究区域

如图2所示，研究区域黄河源区位于32°N—36°N、95°E—103°E之间，青藏高原东部，流域面积为12.3万km2，高程范围为2 656～6 253 m，下垫面以草地、林地为主。黄河源区属于典型的高原大陆性气候，冷热两季交替、干湿两季分明。黄河源区处于印度洋南亚季风和太平洋东亚季风控制区的边缘交错带上，降水同时受暖湿西南季风和东亚季风的支配，多年平均降水量约为510 mm，呈现由东南向西北递减的分布特征。流域出口为唐乃亥水文站。黄河源区是黄河流域的重要产流区和水源涵养区，高精度的径流模拟与预报可为黄河流域水旱灾害防御和水资源管理提供参考。

图2 黄河源区流域及其水文站位置Fig.2 Study area and the location of the hydrological station

2.2 数据

本文构建3个模型的输入数据主要包括：

(1) 降水。融合了GLDAS、TRMM 等多种数据集及国家气象站点数据的CMFD(China Meteorological Forcing Dataset)数据[24-25]，空间分辨率为0.1°。

(2) 气温。来源于全球再分析数据集ERA5的2m AGL温度数据[26]，空间分辨率为0.1°。

(3) 潜在蒸散发。来源于全球再分析数据集ERA5的潜在蒸散发数据，空间分辨率为0.1°。

(4) 高程。来源于SRTM的90 m分辨率的数字高程模型数据(http:∥srtm.csi.cgiar.org/)。

(5) 土壤。来源于联合国粮食和农业组织(FAO)和国际应用系统分析研究所构建的世界土壤数据库(HWSD)。

(6) 归一化植被指数(NDVI)和叶面积指数(LAI)分别为MODIS的MOD13A2和MOD15A2H数据。

其中，LSTM和PCA-LSTM模型输入均为降水和气温数据，LSTM模型的输入为流域空间平均数据；PCA-LSTM模型的输入为流域二维网格数据。THREW模型输入为以上所有数据。此外，用于深度学习模型训练与物理水文模型率定的目标为唐乃亥水文站径流实测数据(1982—2014年逐日数据，来源于西南源区重大研究计划集成项目)。气象与径流数据的时间长度为1982—2014年，时间分辨率为日。

2.3 模型参数设置

(1) PCA-LSTM模型。基于相关研究和前期实验，设置降水和气温为输入数据，实测径流为训练目标。其中，1983—2004年为训练期，2006—2009年为验证期，2011—2014年为测试期。为了评估不同主成分阈值对径流模拟的影响，主成分阈值分别设置为85%、90%和95%。PCA-LSTM模型中LSTM的主要参数设置为：LSTM层数为2，隐藏层数目为32，时间序列长度为120 d。模型损失函数为均方根误差，优化器采用Adaptive Moment Estimation (Adam)。模型批尺寸为32，迭代轮次为200。LSTM输入数据采用“最大值-最小值”方法进行归一化。

(2) LSTM模型。模型输入为流域空间平均的降水和气温数据，其余设置与PCA-LSTM模型的LSTM部分保持一致。

(3) THREW模型。黄河源区划分83个REWs，率定参数包括上、下层土壤饱和导水率平均值等17个参数[19]，率定期设置为1983—2004年，测试期设置为2011—2014年，预热期为2 a。

(4) LSTM和PCA-LSTM模型对输入数据不确定性的鲁棒性评估。为模型中LSTM的输入数据添加不同高斯噪音，评估模型性能随着噪音增加的变化趋势，若模型性能随着噪音平滑增大平滑下降，表明模型未出现数据过拟合的现象，模型具有很好的鲁棒性。参考Kratzert等[27]的相关研究，设置10组符合正态分布的噪音，均值为0,标准差分别为0.02、0.04、…、0.18、0.20。对于每组噪音，生成50条随机噪音向量，以此降低噪音自身随机性对模型性能的影响。

以上所有深度学习程序均由Python的Pytorch库编译实现。

2.4 评价指标

本文选取纳什效率系数[28](ENS)及其3个分解变量[29](相关系数(r)，标准差偏差(α)和均值偏差(β))用于系统评估模型性能。r用于评估模拟值与实测值的趋势一致性，α与β分别评估模拟值标准差与平均值的偏差。

(1)

(2)

(3)

(4)

式中：Qsim，l(Qobs，l)分别为时刻l的模拟(观测)径流;L为模拟期长度;σsim(σobs)和Qsim(Qobs)分别为模拟(观测)的径流的标准值和均值。此外，基于观测值的超越概率设置洪峰偏差(Bp)、中值偏差(Bm)和低值偏差(Bl)，评估模型在不同流量条件下的性能表现[30]。

(5)

式中：i分别代表p(洪峰，超越概率为0～0.02)、m(中值，超越概率为0.3～0.7)和l(低值，超越概率为0.7～1);Li为模拟期中以上不同超越概率对应的数据长度。

3 结果与讨论

3.1 PCA-LSTM模型适用性评估

由于气温数据的空间相关性较高，第一主成分变量的方差累计贡献率达96%，本文设置的3个主成分阈值对应的气温主成分数量均为1。在主成分阈值设置为85%、90%和95%时，降水主成分数量分别为9、14和29。图3展示了降水和气温的第一主成分变量与空间平均值的相关关系，二者在降水与气温的相关系数r2>0.99，表明降水与气温的第一主成分变量包含信息与其流域空间平均值基本一致。

图3 黄河源区降水与气温的第一主成分变量与流域空间平均值比较Fig.3 Comparisons between precipitation and temperature data of first principal component and basin spatially average

图4展示了PCA-LSTM模型在主成分阈值分别设置为85%、90%和95%时唐乃亥水文站径流的模拟与实测径流结果，评价指标结果如表1所示。3种不同主成分阈值的PCA-LSTM模型模拟ENS分别为0.92、0.92和0.93，r分别为0.96、0.96与0.97，模拟径流与实测径流具有很好的一致性；此外，均值偏差结果说明3个主成分阈值对应的模拟径流均存在略微低估，低估小于5%；其中，洪峰、中值与低值分别出现低估、低估与高估，偏差绝对值均低于20%。说明模型偏差较小，能够准确模拟径流量，且不同主成分阈值的模拟结果之间差异较小。综合以上指标表明本文提出的PCA-LSTM模型能够准确捕捉水文气象要素的时空特征，在黄河源区径流模拟表现优异。此外，主成分阈值分别为85%、90%和95%时，ENS均超过0.90，模拟表现无明显差异，说明主成分阈值为85%时PCA方法能够提取满足黄河源区降水—径流过程模拟的气象空间特征，进一步增加主成分阈值对径流模拟影响较小。

图4 不同主成分阈值下PCA-LSTM模型模拟和实测径流过程线Fig.4 Performance of streamflow simulations of PCA-LSTM models with various thresholds in the test period

表1 不同主成分阈值下PCA-LSTM模型评价指标结果Table 1Evaluation metrics results of PCA-LSTM models with different PCA thresholds

3.2 不同模型结果的对比

图5和表2展示了PCA-LSTM模型(主成分阈值为85%)与比对模型THREW和LSTM的模拟径流过程和评价指标结果。THREW、LSTM与PCA-LSTM模型的ENS分别为0.75、0.88和0.92，r分别为0.91、0.95和0.96，表明3个模型均能较好模拟径流趋势，PCA-LSTM模型的表现最为优异，LSTM模型次之，THREW模型最差。此外，PCA-LSTM模型径流量偏差(-2%)远低于LSTM模型(-11%)与THREW模型(-23%)；比较不同阶段的径流偏差可以发现，PCA-LSTM与LSTM模型的低值和中值偏差远低于THREW模型，PCA-LSTM模型峰值与低值偏差在3个模型中最低，但3个模型在径流峰值模拟中均存在高于10%的低估。

表2 不同模型评价指标结果Table 2Evaluation metrics results of various models

图5 不同模型模拟和实测径流过程线Fig.5 Performance of streamflow simulations of various models in the test period

物理水文模型由于对下垫面特征、含水层结构与部分水文过程刻画不足，径流过程模拟总体存在一定误差。径流低水过程阶段较小的绝对误差会产生较大的相对误差，使得THREW模型对于低水过程模拟的相对偏差显著高于其他阶段。而深度学习模型是直接基于历史气象与径流观测数据训练拟合降水—径流过程[6]，由于强大的数据挖掘和拟合能力，LSTM和PCA-LSTM模型径流模拟精度高于THREW模型。深度学习模型的优异性能证明了基于现有水文气象数据可以挖掘更为准确的降水—径流关系，解析深度学习水文模型内部结构，揭示隐藏的降水—径流过程机理规律，将有助于推动对水文过程的认识与提高物理水文模型的模拟精度。另一方面，对于深度学习模型，研究区枯水期的气象与径流年际差异较小，降水—径流过程关系较为单一，而在丰水期降水与径流的年际差异显著，降水—径流过程关系更为复杂，使得丰水期降水—径流过程训练难度更大。此外，研究区下垫面空间异质性明显，产汇流机制差异大，因此降水的空间分布对径流高水过程影响显著。LSTM模型只利用黄河源区降水与气温的流域平均值作为输入，忽略了降水的空间特征，导致径流峰值模拟存在较大偏差。本文提出的PCA-LSTM模型考虑了降水与气温的空间特征，与LSTM模型相比，ENS、径流偏差与峰值偏差分别提高了5%、9%和7%，表明气象要素空间特征可以提升深度学习水文模拟表现。

3.3 模型鲁棒性评估

图6展示了LSTM和PCA-LSTM模型在添加不同标准差的高斯噪音后的模拟结果。随着高斯噪音标准差的增加，ENS降低，变化范围增大，表明模拟表现逐渐降低，噪音自身的随机性对模拟结果的影响增加。但是模拟表现随着噪音光滑增加呈现光滑下降而不是骤降趋势，表明PCA-LSTM模型未对作为LSTM模型输入的主成分变量过拟合，具有良好的鲁棒性。

图6 基于高斯噪音的LSTM与PCA-LSTM模型鲁棒性评估结果Fig.6 Performance of LSTM and PCA-LSTM models with increasing Gaussian noises added to LSTM inputs

4 结 论

(1) 本文提出了融合气象要素时空特性的深度学习水文模型PCA-LSTM，利用主成分分析提取了气象要素的空间特征，基于长短时记忆神经网络学习了长时序过程规律，提高了深度学习水文模型对气象要素空间异质性的刻画能力。

(2) 利用比对模型THREW和LSTM在黄河源区对PCA-LSTM模型进行系统评估，结果表明PCA-LSTM模型径流模拟纳什效率系数达到0.92，较比对模型THREW和LSTM分别提高了23%和5%，证明PCA-LSTM模型具有良好的适用性与鲁棒性，在深度学习模型中引入气象要素空间特征能够提升径流模拟精度。

本文主旨是在LSTM水文模型中考虑气象要素空间信息，但气象要素空间特征提取方法较多，后续将利用图像识别领域神经网络开展进一步研究。此外，物理水文模型由于具有明确的物理机制，仍然是实际水文预报中最常用的工具，如何深入耦合物理水文模型与深度学习模型，同时发挥二者优势，值得未来深入研究。