考虑气象因子不确定性的参考作物蒸散量预报方法

韩 信，张宝忠，魏 征，李益农，陈 鹤

（1.中国水利水电科学研究院 流域水循环模拟与调控国家重点实验室，北京 100038；2.中国农业大学 水利与土木工程学院， 北京 100083）

1 研究背景

参考作物蒸腾量（ET0）是水资源与灌溉管理的重要数据，是水分平衡和灌溉调度的关键因素［1-3］，其计算和预测方法不仅成为农田生态系统水分循环和水平衡研究的重要领域，同时在农田灌溉制度优化和水土资源配置等方面具有重要作用［4］。因此，为了更好地管理作物灌溉用水量以及提高作物水分利用效率，亟待对作物ET0进行准确预测。

1998年FAO-56推荐的Penman-Monteith（PM）模型是计算ET0的标准计算方法，能够得到较精确的结果［1］。ET0与气候变化等自然因素密切相关［5］。近年来涌现出诸多的ET0预测方法，如时间序列法［6-7］、灰色模型法［8］和经验公式法［9-10］。由于ET0与气象因子之间存在复杂的非线性关系，常规的预测模型忽视其复杂的非线性关系造成预测精度不高，因此基于神经网络建立预测模型能够考虑众多因素对ET0的影响［11-12］，具有良好的应用推广价值。且径向基人工神经网络（RBF）较Back Propagation（BP）神经网络具有更好的适用性［13-14］。Ozgur［15］用Levenberge Marquardt（LM）训练算法建立了基于RBF模型的ET0预测，并得到较好的结果，Sudheer等［16］利用有限的气象资料建立了估算水稻日蒸散量的RBF模型，模型预测性能可靠。但是影响ET0的气象因子本身也具有不确定性［17］，因此，近年来不确定性分析逐渐受到关注，但是考虑气象因子不确定性的影响对ET0预测的研究较少。

本研究以华北平原农田下垫面为研究对象，以研究区1981—2018年气象因子和ET0（PM模型计算）资料序列为基础，在贝叶斯理论框架下，根据马尔科夫蒙特卡罗模拟与自适应采样算法结合的方法修正ET0影响因子，基于RBF模型构建ET0不确定性预测模型。

2 预报方法简介

2.1 研究区概况与数据来源 北京市大兴区（39°26′—39°51′N，116°13′—116°43′E）位于华北平原北部永定河冲击平原，总面积1031 km2，属于温带半湿润季风气候。季节变化明显，夏季短而炎热，冬季长而寒冷。多年平均气温12.1℃，平均最高气温14.8℃，最低气温6.7℃，平均相对湿度53%，年平均风速1.84 m/s，年平均蒸散量1080 mm，多年年平均降雨量540 mm，7—9月降雨较多，降雨量占全年降雨总量的80%以上。北京市大兴区下垫面主要为农田，包括玉米/小麦、大豆，其中以玉米和小麦轮作为主，冬小麦整个生育期为260天左右（10月1日—次年6月30日），在正常年份冬小麦需补充灌溉，以保证作物对水分的需求。夏玉米生长期约90天左右（7月1日—9月30日），夏玉米全生育期不灌水［18］。

从中国气象科学共享服务网（http：//cdc.cma.gov.cn）收集了北京市大兴区1981—2018年的逐日历史气象数据和2013—2018年的逐日天气预报数据，建模数据系列长度统一采用1981—2018年共38年的1—12月逐月数据。气象数据包括最高温度、最低温度、平均温度、平均相对湿度、最小相对湿度、日照时数、最大风速、极大风速、平均风速、参考作物蒸散量。天气预报数据及信息包括：最高温度、最低温度、相对湿度及天气情况等，为了保持前后因子的一致性，最终选取最高温度（Tmax）、最低温度（Tmin）、相对湿度（RH）和平均风速（Ws）建立ET0预测模型。天气预报信息中Tmax、Tmin、RH可以直接利用，Ws则需要根据研究区天气情况预报和气象标准进行解析［19］。根据风速对地面的影响而引起的各种现象，按照风速等级确定距离平面高处2 m的平均风速（表1）。

表1 风速等级

2.2 模型方法

2.2.1 参考作物蒸散量标准值 以从中国气象科学共享服务网（http：//cdc.cma.gov.cn）获取的研究区实测的逐日参考作物蒸散量（ET0）为标准值。

2.2.2 基于气象因子不确定性的RBF模型构建 利用径向基人工神经网络模型（RBF-ANN）建立气象因子与ET0的映射关系，并将其作为贝叶斯概率预报系统的似然函数，以ET0标准值为后验信息，对ET0影响因子先验信息进行贝叶斯修正，利用基于自适应算法的马尔可夫链蒙特卡罗随机模拟算法获取各气象因子的后验密度。最后将修正后的因子结合RBF-ANN模型构建考虑气象因子不确定性的ET0预报模型（CU-RBF）。

（1）因子不确定性修正。基于贝叶斯概率预报系统对气象因子的不确定性进行修正。贝叶斯概率预报系统由Krzysztofowicz［20］提出，是一种通用的概率水文预报理论，可以与确定性水文模型耦合，考虑因子不确定性的ET0预测是在贝叶斯推断的基础上提出的，计算式为：

式中：h(y|x)为后验分布；p(x|y)为x对y的条件概率；f(y)为y的边缘密度，即先验分布。

假设气象数据X（x1为最高温度、x2为最低温度、x3为平均相对湿度、x4为平均风速）的先验密度为g(x)；ET0预测模型输出参考作物蒸散量Sn已知时，X的似然函数为 f(S |X )。利用贝叶斯理论修正影响因子的不确定性，输出Sn的期望值为：

式中：φ为X(X ={x1，x2， x3})的水文不确定性的定量表达；k为归一化常数。

（2）修正后因子计算。AM算法是由Haario等［21］提出的一种有效的MCMC采样器，其主要特点是随机抽样的参数空间初始协方差矩阵由先验密度确定的多维正态分布，AM算法的主要步骤如下：

（a）初始化，i=0；

（b）计算转移密度的协方差矩阵：

（3）CU-RBF模型构建。径向基函数网络（RBF）由Broomhead提出［22］。该模型具有较快的收敛速度，在很多领域得到广泛应用［23-24］。RBF网络模型主要包括输入层、隐含层和输出层，RBF模型的表达式为：

φ为可逆高斯函数：

式中：q为隐含层节点个数即样本容量；wi为隐藏层各节点到输出层F(X ) 的权因子；ci为第i个径向基函数的中心值；δ为径向基函数的变量，决定了该径向基函数围绕中心点的宽度。

模型的建立通过Matlab2017b的RBF工具箱来实现，具体函数的调用格式为：

式中：P、T分别为输入和输出变量；goal为均方误差，取goal=0.001，避免过度拟合；spreads为RBF函数的宽度，通过隐含层神经元的中心距离确定；MN为ANN学习过程中神经元最大数目；DF为2次显示之间添加神经元的个数。

仿真函数调用为：

式中：Y为预测目标值；X为输入测试矩阵，net为式（9）中已经学习好的网络。

（4）训练集、验证集及测试集。选择1981—2018年1—12月444组修正后的历史气象月数据（Tmax、Tmin、RH和Ws）和以从中国气象科学共享服务网（http：//cdc.cma.gov.cn）获取的研究区逐日ET0累加得到的月ET0作为训练集；选择2013—2015年36组修正后的历史天气预报月数据（根据天气预报日值累加所得）和月ET0数据作为验证集；选择2016—2018年36组修正后的历史天气预报月数据（根据天气预报日值累加所得）作为测试集。

2.3 模型评价 对于使用RBF神经网络模型预测值与PM计算标准值之间的相互比较，使用了5个评价指标，包括决定系数（R2）、纳什系数（NSE）、均方根误差（RMSE）、平均绝对误差（MAE）和相对误差（RE）来评估模型的预测性能［25］，计算如下：

式中：xi为ET0预报值（利用修正后的气象数据结合径向基神经网络预测所得）；yi为ET0标准值（从中国气象科学共享服务网（http：//cdc.cma.gov.cn）获取的研究区逐日ET0）；y为标准值均值；i为预报样本序列，i=1，2，…，x；n为预报值的样本数。

3 结果与分析

3.1 数据集及先验密度 图1和表2为华北平原农田下垫面1981—2018年气象因子和ET0的数据系列，气象因子包括最高温度（Tmax）、 最低温度（Tmin）、相对湿度（RH）、平均风速（Ws）和ET0。Tmax、Tmin、RH、Ws和ET0多年平均值分别为15.15℃、6.74℃、44.20%、1.89 m/s和1063 mm，多年偏差系数分别为-0.36、-0.34、0.65、-0.36和-0.68，多年的气象因子和ET0都在均值周围有微小的波动，标准偏差分别为0.49℃、0.53℃、3%、0.16 m/s和0.1 mm。

表2 气象参数和ET0

图1 研究区1981—2018年气象参数及ET0资料系列

图2为气象因子的正态概率分布。由图2可知，各气象因子均近似服从正态分布，通过SPSS软件分析所知各因子的正态性假设均通过5%的显著性水平。ET0影响因子的先验密度和标准方差见表3，由表3可以看出，各月Ws的分布差异性较大，变差系数（CV均值为0.13）比Tmax（CV均值为0.11）、Tmin（CV均值为0.06）和RH（CV均值为0.07）的大，表明Ws相对其他气象因子的不确定性较强。

表3 各月ET0影响因子的先验密度

图2 各气象因子的正态概率图

3.2 气象因子不确定性识别 由于本研究所选气象因子为4个变量，即Tmax、Tmin、RH、Ws，因此C0矩阵对角线元素取为各变量的先验方差。设定AM-MCMC算法的初始抽样次数为2000，运行次数38次，每次采样10 000个，这样共抽取有效模拟样本数为304 000个，用于每月气象因子后验分布的统计分析。气象因子在迭代2000次后已经收敛，所抽取的样本已具有其总体的统计特征。

分别以3月（作物生长早期）、7月（作物生长中期）和10月（作物生长后期）典型月为例，所抽取各影响因子的后验密度直方图见图3—图5。Tmax、Tmin、RH、Ws的后验密度分别为：最高温度x1～N（10.46，0.982），最低气温x2～N（1.49，0.562），相对湿度x3～N（34.62，0.022），平均风速x4～N（2.27，0.132）；最高温度x1～N（25.73，0.522），最低气温x2～N（18.54，0.352），相对湿度x3～N（58.97，0.012），平均风速x4～N（1.65，0.062）；最高温度x1～N（16.03，0.682），最低气温x2～N（7.03，0.412），相对湿度x3～N（49.69，0.012），平均风速x4～N（1.60，0.072）。对于Ws先验密度的修正较Tmax、Tmin和RH大，在相同的显著性水平条件下，各影响因子的置信区间不同，其中Ws比Tmax、Tmin和RH的范围大，说明Ws具有较强的不确定性，主要是因为Ws受到外界的干扰较为敏感。因此，后验信息对Ws的修正程度较其他3个因子大。Pereira等［4］的研究结果也表明为了避免因Ws的不确定性带来的作物耗水预测不确定性，在进行耗水预测时将不考虑Ws，而仅考虑T和RH等因子。由此可见，考虑Ws进行ET0预测时，对其不确定进行修正是必要的。

图3 3月份频率直方图和密度曲线

图4 7月份频率直方图和密度曲线

图5 10月份频率直方图和密度曲线

3.3 参考作物蒸散量预测

3.3.1 单气象因子不确定性修正对ET0预测影响 由3.2节可知，在相同的显著性水平条件下，各气象因子的均值置信区间范围不同，其中Ws比Tmax、Tmin和RH的范围大，说明Ws表现出较强的不确定性。图6和表4为预测期单独考虑修正某一因子（Tmax、Ws和RH）后RBF预测ET0的效果。

由图6和表4可知，当单独对Tmax修正后，基于RBF模型的ET0的预测值与标准值的变化趋势基本一致，且效果较未修正前好，模型评价参数效果均有所提高，R2、NSE分别增加0.04、0.01，RMSE、MAE和RE分别减少0.72 mm、0.02 mm和0.22%，当分别单独考虑Ws和RH修正后，基于RBF模型的ET0的预测值与标准值的变化趋势基本一致，效果均较未修正前好，模型评价参数效果均有所提高。R2、NSE分别增加0.06、0.01和0.02、0.01，RMSE、MAE和RE分别减少6.08 mm、0.17 mm、1.71%和0.21 mm、0.01 mm、0.14%。由图6对比可知，单独对Ws修正，相较于分别单独修正Tmax和RH时的预测效果较好，可能是因为风速与气象因子以外的农田下垫面、观测高度等有关，修正后的Ws能够削弱外界条件对其的不确定性影响。而单独对RH进行修正时预测效果较差，R2、NSE分别相差0.05、0.02，RMSE、MAE和RE分别相差5.87mm、0.16mm、1.57%，可能是因为RH与降雨有直接关系，该研究结果与Mohmmad、刑贞相等研究结果一致［26-27］，在降雨条件下，各个气象因子的不确定性增加，未对其他气象因子进行修正必然会造成预测的误差。

图6 最高温度、相对湿度和平均风速分别修正前后ET0预测效果

表4 分别单独修正气象因子不确定性修正前后ET0预测效果

3.3.2 考虑气象因子不确定性修正的ET0预测 以AM-MCMC算法随机模拟每月8000组ET0影响因子相应的RBF预报值为对应月份ET0的预测集合，将均值作为该月的ET0预测值，并得出指定概率（85%）的置信区间，实现了考虑了RBF模型输入不确定性对ET0预报不确定性的影响（CU-RBF）。其中2013—2015年为验证期，2016—2018年为测试期，验证期和测试期的各月ET0计算与实测过程以及RBF的确定性和不确定性预报过程分别见图7和图8。

图7 验证期和测试期RBF、CU-RBF与标准值相关性分析

图8 验证期和测试期RBF、CU-RBF预测值与标准值对比

由图7（a）、图8（a）和表5可以看出，验证期（2013—2015）和测试期（2016—2018）在未考虑气象因子不确定性下建立的ET0预测模型（RBF）预测值与基于PM模型计算的标准值的相关性R2分别为0.87和0.89，拟合回归系数分别为0.83和0.81，NSE分别为0.88、0.93。RBF确定性模型预测得到的ET0与PM模型计算的标准值趋势基本一致，月均值分别为91.28 mm和88.98 mm，均方根RMSE和MAE分别为50.91 mm、46.27 mm和1.41 mm、1.29 mm；均值RE分别为6.99%、3.34%。由图7（b）、图8（b）和表5可以看出，验证期和测试期在考虑气象因子不确定性下建立的ET0概率预测模型（CURBF）预测值与基于PM模型计算的标准值的相关性的R2分别为0.94和0.97，拟合回归系数分别为1.18和0.99，NSE分别为0.90、0.99。CU-RBF确定性模型预测得到的ET0与PM模型计算的标准值趋势基本一致，相较于单独考虑Tmax、RH、Ws的预测模型效果较好，NSE增加0.08；月均值分别为92.18 mm和91.05 mm，RMSE和MAE分别为47.28 mm、38.43 mm和1.31 mm、1.07 mm；均值RE分别为1.33%、-0.02%。综上分析可知，无论在验证期和测试期，CU-RBF模型预测的ET0与标准值的R2、NSE均有所提高，RMSE、MAE和RE均有所下降。R2、NSE分别提高0.07、0.08和0.09、0.06，RMSE、MAE和RE分别下降3.36 mm、7.84 mm，0.1 mm、0.22 mm和5.36%、3.36%。且CU-RBF能给出预报值的85%置信区间，其中ET0标准值包含在指定的概率区间内的比例较高。Sudheer、邢贞相和德佳硕等［16，27-28］研究同样表明气象因子的不确定性必然影响作物需水预测的精度，在进行基于神经网络模型对参考作物蒸散量预测前对气象因子的不确定性修正能够较大程度的提升其预测精度。本研究验证了气象因子不确定性修正在参考作物蒸散量预测方面实际应用的可行性，可为农田下垫面的未来水分管理提供科学依据。

表5 CU-RBF和RBF逐月ET0预测精度评价

4 结论

本研究通过分析华北平原区农田下垫面气象因子和ET0的变化情况，结合贝叶斯理论，利用基于自适应采样算法的AM-MCMC求解ET0影响因子（最高温度、最低温度、相对湿度、平均风速）的后验概率密度，应用RBF-ANN建立ET0与气象因子的函数关系，进行ET0预报，得到以下结论：

（1）各气象因子均近似服从正态分布，且Ws的变差系数比Tmax、Tmin和RH的大，相对其他主要气象因子的不确定性较强。

（2）典型月份（3月、7月和 10月）的后验密度分别为：Tmax～N（10.46，0.982），Tmin～N（1.49，0.562），RH～N（34.62，0.022），Ws～N（2.27，0.132）；Tmax～N（25.73，0.522），Tmin～N（18.54，0.352），RH～N（58.97， 0.012）， Ws～N（1.65， 0.062）； Tmax～N（16.03， 0.682）， Tmin～N（7.03， 0.412）， RH～N（49.69，0.012），Ws～N（1.60，0.072）。

（3）单独修正Ws不确定性的CU-RBFWs模型比分别单独修正Tmax和RH的CU-RBF预测模型好，R2、NSE分别为0.95、0.92，RMSE、MAE和RE分别为40.19 mm，1.12 mm和1.63%。

（4）在验证期和测试期，考虑气象因子不确定性构建的CU-RBF模型的预报精度均有很大程度的提高。R2、NSE分别提高0.07、0.08和0.09、0.06，RMSE、MAE和RE分别下降3.36 mm、7.84 mm，0.1 mm、0.22 mm和5.36%、3.36%，且CU-RBF能得到85%置信区间的预报值。