第34卷第1期2023年1月㊀㊀
正则化是结构风险最小化策略的实现水科学进展ADVANCES IN WATER SCIENCE Vol.34,No.1Jan.2023
DOI:10.14042/jki.32.1309.2023.01.004基于变分贝叶斯深度学习的水文概率预报方法
李大洋1，姚㊀轶2，梁忠民3，周㊀艳4，李彬权3,5
(1.盐城工学院土木工程学院，江苏盐城㊀224051;2.南京水利科学研究院水文水资源与水利工程科学国家重点实验室,
江苏南京㊀210029;3.河海大学水文水资源学院，江苏南京㊀210098;4.盐城师范学院江苏省滩涂生物资源与环境保护重点实验室,
江苏盐城㊀224007;5.水安全与水科学协同创新中心，江苏南京㊀210024)摘要：目前水文领域关于深度学习的研究多集中于提高预测能力方面，与具有物理机制的水文模型相比，深度学习复杂的内部结构导致其不具备可解释性，预测结果难以被信任，因此发展可信赖的深度学习对于推进水科学发展具有重要意义㊂基于预报残差分析框架，构建具有物理机制的水文模型与深度学习融合的混合模型，以充分利用两者优势；引入变分贝叶斯理论，提出变分贝叶斯与深度学习耦合的概率预报模型VB-LSTM，以定量评估水文预报结果的不确定性㊁提高结果可靠度㊂以黄河源区1961 2015年的径流过程为研究对象，对VB-LSTM 模型
进行应用示例研究㊂结果表明：与长短时记忆网络(LSTM)相比,VB-LSTM 模型在验证期预报精度更高，结果更稳定；与传统基于 线性-正态 假设的水文概率预报方法相比,VB-LSTM 模型具有更高的预报精度，且不确定性更小㊁预报结果更可靠㊂
关键词：水文概率预报；深度学习；变分贝叶斯；长短时记忆网络；混合模型
中图分类号:TV124;P338㊀㊀㊀文献标志码:A㊀㊀㊀文章编号:1001-6791(2023)01-0033-09
收稿日期:2022-10-14；网络出版日期:2023-02-01
网络出版地址:https :ʊknski /kcms /detailʊ32.1309.P.20230131.1533.005.html
基金项目：国家自然科学基金资助项目(41730750;41877147)
作者简介：李大洋(1992 )，男，江苏徐州人，博士，主要从事水文水资源方面研究㊂
E-mail:dayangli_hhu@foxmail 通信作者：梁忠民,E-mail:zmliang@hhu.edu 1990年以来，机器学习在水文领域不断应用发展㊂人工神经网络(ANN)㊁循环神经网络(RNN)㊁支持向量机(SVM)㊁最近邻算法(KNN)等一系列模型被广泛应用于降水预报修正㊁径流预报㊁土壤水分预测㊁地下水位预测等水文领域的各个方面[1-2]㊂近些年，深度学习作为一种特定类型的机器学习日益受到水文学者的广泛关注㊂
与传统的 浅层 机器学习相比，深度学习通过复杂结构或多处理层进行多层次㊁多维度的数据抽象处理，一些研究显示深度学习模型如长短时记忆网络(LSTM)在水位㊁径流㊁土壤水等方面展现出了更强的预测能力[3-5]㊂然而，与具有物理机制的水文模型相比，深度学习无法真正理解水文过程及其物理规律，其结构复杂难以被解释，导致预测结果难以被充分信任，如何构建可信赖的深度学习是当前研究的热点与难点[6]㊂预报残差分析模型[7]作为一种在水文领域应用广泛的概率预报后处理模型，提供了将具有物理机制的模型与深度学习模型相结合从而充分发挥两者优势的途径，有望在提高深度学习预测结果可信赖程度方面发挥作用[8]㊂贝叶斯理论是进行水文概率预报的通用框架[9-10]，将贝叶斯理论应用于深度学习可以实现概率预报并对预测结果可靠度进行评估，然而其中常用于估计参数后验分布的抽样算法  马尔可夫链蒙特卡罗法(MCMC)的收敛速度偏慢，消耗计算资源较多，难以进行高维参数估计㊂Vrugt 等[11]指出采用MCMC 算法估计地下水模型MODFLOW 中的241个参数需花费数月才能完全收敛㊂深度神经网络的参数数量一般远高于这个量级，难以采用MCMC 算法在有效时间内(如1~3d 预见期)估计众多参数的不确定性并提供概率预报㊂
Fang 等[12]通过蒙特卡罗暂退算法(Monte Carlo Dropout,MC Dropout)估计深度神经网络LSTM 参数的不确定性，实现了土壤水的概率预报㊂MC Dropout 是一种将正则化算法Dropout 解释为高斯过程的贝叶斯近似，有
34㊀水科学进展第34卷㊀望在不增加计算成本的情况下实现深度学习的概率预报㊂但一些学者认为MC
Dropout算法并非基于贝叶斯理论，无法正确估计参数并提供合理的不确定性分析[13]；此外，该算法专为神经网络设计，无法用于水文领域其他高维参数估计问题㊂变分贝叶斯(Variational Bayes,VB)[14]是解决贝叶斯推断问题的新兴方法，与MC Dropout相比,VB方法理论基础与适用性更强；与MCMC算法相比,VB的计算效率更高㊂VB方法的本
质在于能快速寻参数化分布族的最佳近似分布，因此，对估计高维参数等问题具有显著优势㊂本文引入一种基于变分贝叶斯理论的最新算法  随机变分推断(Stochastic Variational Inference, SVI)[15]，提出一种耦合SVI与LSTM的概率预报方法(VB-LSTM)，作为水文模型的预报残差分析模型定量评估水文预报不确定性㊂为充分发挥物理模型与数据驱动模型的优势，采用分布式水文模型模拟径流过程，研究SVI算法的收敛性和时效性，分析VB-LSTM模型在减少水文预报不确定性和提高预报精度方面的优势㊂1㊀方法原理
1.1㊀预报残差分析模型
预报残差分析模型是一种后处理方法，用于构建预报残差与预报因子矩阵XɪDˑT之间的线性或者非线性回归关系，其中，为实数集,D为所选预报因子的个数,T为时间序列长度㊂此回归关系可写作
yε=f(X,θ)+ε(1)式中:yε=[yε1, ,yεT]ɪ1ˑT为实测流量(记为向量y o)与模拟流量(记为向量y m)之差，被称为残差;ε=[ε1, ,εT]ɪ1ˑT为随机误差项，通常被认为服从均值为0㊁方差为σ2的正态分布；参数向量θ= [θ1, ,
θN]ɪ1ˑN，其中,N为参数个数㊂水文预报不可避免地存在着误差[10,16]，在水文预报残差分析模型中，建模者认为水文预报的所有误差(如数据观测㊁模型初值㊁模型参数㊁模型结构)均集中在yε中，通过分析和模拟残差从而达到量化预报不确定性的目的㊂
贝叶斯学派认为θ是随机变量且可以通过概率分布描述㊂根据贝叶斯公式,θ的后验分布为
p(θ|y o)=p(y o|θ)p(θ)
p(y o)ɖp(y o|θ)p(θ)(2)式中：条件分布p(θ|y o)在贝叶斯中被称为后验分布;p(y o|θ)称为似然函数；边际分布p(θ)称为先验分布;p(y o)称为证据(Evidence)或归一化常数；符号ɖ表示 成正比例 ㊂贝叶斯推断的关键是估计参数后验分布，但这很难直接获得，主要原因在于证据p(y o)=ʏp(y o|θ)p(θ)dθ是一个积分，除个别理想情况外多数不可积㊂以MCMC为代表的贝叶斯抽样方法通过建立参数θ的马尔可夫链依次进行抽取样本㊂这意味着每次迭代过程中每个参数可获得1个抽样样本，但如果不符合一定条件，该样本将被舍弃㊂当参数过多时,MCMC抽样效率很低且难以收敛，需要更多的迭代次数[11]㊂变分贝叶斯方法采用已知的一组分布逼近真实后验的思路极大地提高了高维参数估计的计算效率㊂需要注意的是，变分贝叶斯方法的参数估计结果存在一定的偏差，因此，使用时需要在精度与速度之间做出权衡[17]㊂
1.2㊀变分贝叶斯方法
变分贝叶斯的主要思路是通过引入一组分布去逼近参数后验分布p(θ|y o)，这组分布被称为变分分布，记作q(θ|Λ)，其中Λ=[λ1, ,λN]ɪSˑN是与模型参数θ=[θ1, ,θN]对应的变分参数矩阵,S 意味着每个变分分布由S个变分参数控制㊂例如，假设变分分布为正态分布，则S=2，因为正态分布是由均值(μ)和方差(σ2)2个参数决定，该参数可以被称为变分参数，记作λ=[μ,σ2]㊂在变分贝叶斯计算中,Kullback-Leibler散度(D KL)，又称相对熵，被用来量化变分分布与真实后验之间的距离㊂D KL可推导后表示为:
D
KL[q(θ|Λ)p(θ|y o)]=-(Λ)+lg p(y o)(3)
㊀第1期李大洋，等：基于变分贝叶斯深度学习的水文概率预报方法35㊀
(Λ)=E
q[lg p(y o,θ)-lg q(θ|Λ)](4)式中:E
q(㊃)为关于q(θ|Λ)的期望;p(y o,θ)为y o与θ的联合分布;(Λ)为证据下界(Evidence Low-er BOund,ELBO)，相当于变分贝叶斯方法的目标函数㊂
调整式(3)各项位置后可推导出
log p(y
o)=D KL[q(θ|Λ)p(θ|y o)]+(Λ)ȡ(Λ)(5)㊀㊀从在变分贝叶斯中,D KL被用来计算变分分布与真实后验之间的距离，最小化D KL可以最终得到近似的
后验㊂最小化D KL的难点在于证据p(y o)难以被积分，导致D KL很难被计算㊂从式(3)可以看出，最小化D KL 等价于最大化ELBO，因此，变分贝叶斯可以通过最大化ELBO的方式来获得后验分布的近似㊂最大化EL-BO可以通过梯度优化算法获得，但ELBO中包含期望函数(见式(4))导致很难进行求导㊂Ranganath等[15]提出了一种目前最先进的变分贝叶斯方法  随机变分推断，该算法采用REINFORCE技巧解决此问题，详细介绍请见参考文献[15]㊂
1.3㊀变分贝叶斯与长短时记忆网络耦合模型  VB-LSTM
神经网络通常采用基于反向传播的梯度下降算法更新权重并最小化目标函数，得出最优的确定性预报结果，一般可称作确定性神经网络㊂然而，确定性神经网络中的众多参数使得预报结果具有明显的不确定性，预报结果难以被信任㊂贝叶斯理论提供了量化不确定性的基本框架，将贝叶斯理论用于神经网络参数不确定性量化并提供概率预报的方法可被称为贝叶斯神经网络[18]㊂图1比较了确定性神经网络与贝叶斯神经网络的结构差异㊂本文采用变分贝叶斯训练先进的深度学习网络LSTM，构建耦合模型  变分贝叶斯深度学习网络VB-LSTM，作为水文预报中的预报残差分析模型并提供概率预报结果㊂
VB-LSTM模型由深度神经网络和变分贝叶斯方法两部分构成，因此，在构建VB-LSTM模型时需考虑两者的特点㊂进行变分贝叶斯估计时，需要在构建LSTM模型的基础上，合理地设置参数的变分分布㊁先验分布以及似然函数㊂下面具体介绍VB-LSTM模型的各个部分㊂
图1㊀确定性神经网络与贝叶斯神经网络比较
Fig.1Comparison of deterministic neural networks and Bayesian neural networks
36㊀水科学进展第34卷㊀1.3.1㊀长短时记忆网络
长短时记忆网络是Hochreiter等[19]在1997年提出的一种特殊的循环神经网络㊂LSTM在RNN基础上引入了 门 ，在一定程度上解决了标准RNN存在的梯度消失或爆炸以及长时记忆问题㊂LSTM网络的架构包含若干LSTM层和1个全连接层㊂全连接层一般用线性函数表示，负责在最后输出最终结果㊂每一个LSTM 层由很多隐含单元(又叫神经元)组成，每个隐含单元又包括2个状态：隐含状态(又叫输出状态)和细胞状态㊂细胞状态是由3个门控制，分别是输入门㊁遗忘门㊁输出门㊂门负责决定信息的进入，其结构是一个
Sigmoid层和点乘操作的组合㊂Sigmoid层将输入值转化成0到1之间的数值，数值大小代表了输入值通过的程度，其中0代表不通过,1代表都通过㊂还有一个候选细胞为细胞状态添加额外信息㊂具体公式可见参考文献[19]㊂
1.3.2㊀先验分布、变分分布与似然函数
先验分布代表了进行推理前对该参数或数据的认知㊂由于LSTM参数没有明确的物理意义或可解释性，
因此，很难给出一个有意义的先验分布㊂在这种情况下，一般可以用无信息(Uninformative)的均匀分布或者数学性质较好的正态分布描述㊂这里假定LSTM参数的先验分布服从正态分布，误差项的方差为了保证非负则服从均匀分布，表达式如下:
θi~N(0,1)(6)
σ2~U(0,1)(7)式中:N(0,1)为均值为0㊁方差为1的正态分布;U(0,1)为下界为0㊁上界为1的均匀分布㊂在变分推断中，变分分布用来近似真实后验分布㊂本文假定LSTM参数的变分分布为正态分布，主要因为正态分布简单且具有良好的数学性质㊂为了保证误差项的方差估计为正，采用截断正态分布Nᶄ(㊃)作为变分分布㊂变分推断通过迭代不断更新变分参数(αi㊁γ2i㊁μ和2)从而控制变分分布逼近真实后验分布，迭代时需要设置初始变分分布㊂这里初始变分分布的均值设定为对应的先验分布均值，标准差为对应先验分布的1/10，写作:
θi|αi,γ2i~N(αi=0,γ2i=0.12)(8)
σ2|μ,2~Nᶄ(μ=0.5,2=0.0252)(9)似然函数用条件分布的形式衡量模拟值与实测值之间的吻合程度，似然函数的选择直接决定推断结果的合理性㊂为了使数据减少异方差性以及更服从正态分布，将Box-Cox转换应用于观测数据与模拟数据中㊂然后采用标准化减少变量的尺度差异，帮助SVI算法更快地收敛㊂在数据转换空间ϕ(㊃)中，假设随机误差独立同分布且服从均值为0㊁方差为常数σ2的正态分布㊂似然函数
可以表示为
l=-T2lg2πσ2+12ðT t=1-ε2tσ2+(λ1-1)ðT t=1lg(y t o+λ2)(10)其中:
εt=(y~~t m)-f(x㊆t,θ)(11)式中:y~t o和y~t m分别为t时刻经Box-Cox转换后的观测数据和水文模型模拟值;(y~~t m)表示(y~t o-y~t m)的标准化;x㊆t表示LSTM输入数据的标准化;λ1和λ2为Box-Cox转换的2个参数，根据文献建议固定λ1= 0.2和λ2=0[20]㊂
2㊀应用实例
2.1㊀研究流域与数据处理
本文选取黄河源区即唐乃亥以上流域进行应用研究㊂黄河源区的高程㊁水文站点㊁气象站点及水系分布
㊀第1期李大洋，等：基于变分贝叶斯深度学习的水文概率预报方法37㊀如图2所示㊂流域面积约12.2万km2，约占黄河流域面积的15%㊂本文所需的气象数据包括降水㊁平均气温㊁最高气温㊁最低气温㊁相对湿度㊁平均风速和日照时数等日系列资料，除降水作为水文模型直接输入外其余气象数据主要用于计算潜在蒸散发(PET)㊂考虑到资料的一致性和完整性，本文收集了图中所示的流域内7个气象站和1个水文站的实测日系列资料，时间范围为1961 2015年㊂此外，收集了空间分辨率为90m的数字高程模型(DEM)数据以及比例尺为1ʒ100万的土壤分类数据HWSD V1.1㊂采用反距离插值法将各类型数据插值到
对应计算网格㊂
图2㊀黄河源区流域
Fig.2Map of the source area of the Yellow River
2.2㊀模型设置
本文采用分布式水文模型MIKE SHE模拟黄河源区1961 2015年的日径流过程㊂将模拟时段分成3部分，即预热期(1961 1965年)㊁率定期(1966 2000年)和验证期(2001 2015年)㊂由于研究流域面积较大(约122000km2)，所以模型网格大小设置为5kmˑ5km，总共约5000个有效计算网格㊂采用数据中提取与人工率定相结合的方式率定MIKE SHE的参数，具体可见相关文献[21]㊂MIKE SHE模拟流量作为VB-LSTM的输入变量，将日模拟径流集成的月径流输入VB-LSTM构建月尺度的预报残差分析模型㊂选择月尺度主要有以下2点考虑:①月系列数据比日数据资料长度更短，可以减少计算时间，便于研究随机变分推断SVI的收敛性;②深度学习通常在大样本中表现优异，采用月尺度的小样本数据训练VB-LSTM，如果表现较好则可以在一定程度上说明VB-LSTM的通用性㊂
根据相关文献[5]以及人工试错,VB-LSTM设置为1个包含20个隐含单元的LSTM层，后面接1个全连接层㊂在SVI中，采用梯度下降优化算法Adam最大化证据下界ELBO，其相关超参数采用默认设置如下:ρ= 0.001,β1=0.9,β2=0.999以及 =10-8㊂VB-LSTM最大迭代次数设置为50000次㊂VB-LSTM在Python 环境中运行并基于开源机器学习库Pytorch与概率模型工具库Pyro构建，计算采用服务器集中英特尔CPU
中的4颗核心㊂
2.3㊀结果分析
图3给出了ELBO值(N ELBO)的迭代变化轨迹(NᶄELBO为N ELBO除以数据个数)㊂在SVI中，采用梯度下降优化算法Adam最大化证据下界来估计参数后验，设定最大迭代次数阈值为50000㊂SVI的收敛性可通过观察N ELBO的变化轨迹判断㊂如图3所示,N ELBO的轨迹图不够平滑，有不少噪声㊂为方便观察收敛性，引入窗口为20的滑动平均函数处理N ELBO(黑线条)㊂可以发现当迭代次数超过30000时平滑后的N ELBO稳定不变，这时可认为N ELBO已经收敛，并且收敛时的计算时间少于1h，能够满足一定的预报时效性需求㊂此时，本文

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。