收稿日期：2020-08-12；网络首发时间：2021-01-25
网络首发地址：knski/kcms/detail/11.5020.TV.20210122.1708.003.html
基金项目：国家自然科学基金（51822907）；中国水利水电科学研究院“五大人才”计划（ID0145B742017）作者简介：韩信（1990-），博士生，主要从事灌区需耗水预测预报研究。E-mail ：通讯作者：张宝忠（1981-），教授级高级工程师，主要从事蒸散发尺度效应、现代灌区高效用水理论与技术研究。
E-mail ：文章编号：1672-3031（2021）01-0033-12
中国水利水电科学研究院学报
第19卷第1期
韩信1，2，张宝忠1，魏征1，李益农1，陈鹤1
（1.中国水利水电科学研究院流域水循环模拟与调控国家重点实验室，北京100038；
2.中国农业大学水利与土木工程学院，北京100083）
摘要：参考作物蒸散量（ET 0）的准确预测预报对于制定作物灌溉制度与实时灌溉调度具有重要意义，然而气象因子的不确定性极大的影响着ET 0的预测精度。因此本研究采用马尔科夫蒙特卡罗模拟与自适应采样算法相结合的方法（AM-MCMC ）对气象因子的不确定性进行修正，以气象站实测ET 0作为标准值，利用径向基神经网络（RBF ）模型建立气象因子与ET 0的映射关系，建立基于气象因子不确定性的ET 0不确定性预测模型（CU-RBF ），并以华北平原农田下垫面为例进行验证。结果表明，与传统的RBF 确定性预报结果相比，CU-RBF 预测结果的各精度评价效果均有所提高，纳什系数提高了10%，均方根误差和平均绝对误差分别降低了16.94%、17.05%，且单独修正平均风速的CU-RBF Ws 预测模型效果比分别单独修正最高温度、平均相对湿度的预测模型效果好。考虑气象因子不确定性开展ET 0的预报研究，减小了预测值与实际值的误差，可为农田下垫面的未来水分管理提供科学依据。关键词：参考作物蒸散量；气象因子；不确定性；贝叶斯预报系统；径向基神经网络中图分类号：P338 文献标识码：A
doi ：10.13244/jki.jiwhr.20200167
1研究背景
参考作物蒸腾量（ET 0）是水资源与灌溉管理的重要数据，是水分平衡和灌溉调度的关键因素
［1-3］
，
其计算和预测方法不仅成为农田生态系统水分循环和水平衡研究的重要领域，同时在农田灌溉制度优化和水土资源配置等方面具有重要作用［4］
。因此，为了更好地管理作物灌溉用水量以及提高作物水分利用效率，亟待对作物ET 0进行准确预测。
1998年FAO-56推荐的Penman-Monteith （PM ）模型是计算ET 0的标准计算方法，能够得到较精
确的结果［1］
。ET 0与气候变化等自然因素密切相关［5］
。近年来涌现出诸多的ET 0预测方法，如时间序列法
［6-7］
、灰模型法［8］和经验公式法
［9-10］
。由于ET 0与气象因子之间存在复杂的非线性关系，常规的
预测模型忽视其复杂的非线性关系造成预测精度不高，因此基于神经网络建立预测模型能够考虑众多因素对ET 0的影响
［11-12］
，具有良好的应用推广价值。且径向基人工神经网络（RBF ）较Back Propaga ⁃
tion （BP ）神经网络具有更好的适用性
［13-14］
。Ozgur
［15］
用Levenberge Marquardt （LM ）训练算法建立了基于
RBF 模型的ET 0预测，并得到较好的结果，Sudheer 等
［16］
利用有限的气象资料建立了估算水稻日蒸散
包钢大集体改制量的RBF 模型，模型预测性能可靠。但是影响ET 0的气象因子本身也具有不确定性
［17］
，因此，近年来
不确定性分析逐渐受到关注，但是考虑气象因子不确定性的影响对ET 0预测的研究较少。
本研究以华北平原农田下垫面为研究对象，以研究区1981—2018年气象因子和ET 0（PM 模型计
算）资料序列为基础，在贝叶斯理论框架下，根据马尔科夫蒙特卡罗模拟与自适应采样算法结合的方法修正ET 0影响因子，基于RBF 模型构建ET 0不确定性预测模型。
2预报方法简介
2.1研究区概况与数据来源北京市大兴区（39°26′—39°51′N，116°13′—116°43′E）位于华北平原北部永定河冲击平原，总面积1031km2，属于温带半湿润季风气候。季节变化明显，夏季短而炎热，冬季长而寒冷。多年平均气温12.1℃，平均最高气温14.8℃，最低气温6.7℃，平均相对湿度53%，年平均风速1.84m/s，年平均蒸散量1080mm，多年年平均降雨量540mm，7—9月降雨较多，降雨量占全年降雨总量的80%以上。北京市大兴区下垫面主要为农田，包括玉米/小麦、大豆，其中以玉米和小麦轮作为主，冬小麦整个生育期为260天左右（10月1日—次年6月30日），在正常年份冬小麦需补充灌溉，以保证作物对水分的需求。夏玉米生长期约90天左右（7月1日—9月30日），夏玉米全生育期不灌水［18］。
从中国气象科学共享服务网（http：//v）收集了北京市大兴区1981—2018年的逐日历史气象数据和2013—2018年的逐日天气预报数据，建模数据系列长度统一采用1981—2018年共38年的1—12月逐月数据。气象数据包括最高温度、最低温度、平均温度、平均相对湿度、最小相对湿度、日照时数、最大风速、极大风速、平均风速、参考作物蒸散量。天气预报数据及信息包括：最高温度、最低温度、相对湿度及天气情况等，为了保持前后因子的一致性，最终选取最高温度（T max）、最低温度（T min）、相对湿度（RH）和平均风速（Ws）建立ET0预测模型。天气预报信息中T max、T min、RH可以直接利用，Ws则需要根据研究区天气情况预报和气象标准进行解析［19］。根据风速对地面的影响而引起的各种现象，按照风速等级确定距离平面高处2m的平均风速（表1）。
风级0 1 2 3 4 5 6标准
无风
轻风
轻风
微风
和风
清风
强风
风速/(m/s)
0~0.1
0.2~1.1
1.2~
2.5
2.6~4.0
4.1~
5.9
6.0~8.0
8.1~10.3
均值/(m/s)
0.7
1.5
3.0
5.2
6.7
9.0
风级
7
8
9
10
11
脐点1988
12
标准
劲风
大风
烈风
狂风
暴风
飓风
风速/(m/s)
10.4~12.8
12.9~15.5
15.6~18.3
18.4~21.2
21.3~24.4
24.5~27.6
均值/(m/s)
12.0
14.2
17.2
19.4斯蒂文 斯皮尔伯格
23.2
26.2
表1风速等级
2.2模型方法
2.2.1参考作物蒸散量标准值以从中国气象科学共享服务网（http：//v）获取的研究区实测的逐日参考作物蒸散量（ET0）为标准值。
2.2.2基于气象因子不确定性的RBF模型构建利用径向基人工神经网络模型（RBF-ANN）建立气象因子与ET0的映射关系，并将其作为贝叶斯概率预报系统的似然函数，以ET0标准值为后验信息，对ET0影响因子先验信息进行贝叶斯修正，利用基于自适应算法的马尔可夫链蒙特卡罗随机模拟算法获取各气象因子的后验密度。最后将修正后的因子结合RBF-ANN模型构建考虑气象因子不确定性的ET0预报模型（CU-RBF）。
（1）因子不确定性修正。基于贝叶斯概率预报系统对气象因子的不确定性进行修正。贝叶斯概率预报系统由Krzysztofowicz［20］提出，是一种通用的概率水文预报理论，可以与确定性水文模型耦合，考虑因子不确定性的ET0预测是在贝叶斯推断的基础上提出的，计算式为：
h() y|x=
f()y p()
x|y
f()y p()
x|y d y
（1）
式中：h()
y|x为后验分布；p()
x|y为x对y的条件概率；f()y为y的边缘密度，即先验分布。
福柯知识考古学假设气象数据X （x 1为最高温度、x 2为最低温度、x 3为平均相对湿度、x 4为平均风速）的先验
密度为g ()x ；ET 0预测模型输出参考作物蒸散量S n 已知时，X 的似然函数为f ()S |X 。利用贝叶斯理论修正影响因子的不确定性，输出S n 的期望值为：
k ()S n |X =
-∞+∞
f ()S n |X
g ()X d X
（2）φ()
X |S n =
f ()S n |X
g ()X k ()
S n |X
=
f ()S n |X
g ()
X -∞+∞
f ()S n |X
g ()X d X
（3）
式中：φ为X ()
X ={x 1，x 2，}x 3的水文不确定性的定量表达；k 为归一化常数。
（2）修正后因子计算。AM 算法是由Haario 等［21］
提出的一种有效的MCMC 采样器，其主要特点是
随机抽样的参数空间初始协方差矩阵由先验密度确定的多维正态分布，AM 算法的主要步骤如下：
（a ）初始化，i =0；
（b ）计算转移密度的协方差矩阵：
C i =ìí
îïïC 0
，i ≤t 0s d
()
Cov ()X 0，X 1 X i -1+s d εI d （4）
式中：C 0为初始协方差；ε为一个较小的正数，避免C i 为奇异矩阵；s d 取值为2.42d ；I d 为d 维单位矩阵；t 0为初次抽样次数。
第一次迭代的协方差计算式为：
C i +1=i -1i C i +s d i
()
i ------X i -1X ˉT i -1-()i +1-X i ----X T i +X i X T拦阻索
i +εI d
（5）
式中：------X i -1、-X i 分别为前i -1次、i 次抽样的均值；X T i 为向量X i 的转置。
（c ）从转移密度N ()X iv ，C i 中产生新的样本X *i 。（d ）计算新样本X *i 的接受概率，计算式为：
α=min ìíî
ï
ï
1，
üý
þ禁娱令
ï
ï
f ()S n |X *i
g ()X *
i f (
)
S n |X i g ()X i （6）
（e ）产生一个均匀随机数μ~U []0，1（f ）若μ<α，接受X i +1=X *i ，反之，X i +1=X i
（g ）i =i +1，重复步骤（b ）—步骤（g ），直到达到预定精度。（3）CU-RBF 模型构建。径向基函数网络（RBF ）由Broomhead 提出［22］
。该模型具有较快的收敛速
度，在很多领域得到广泛应用［23-24］
。RBF 网络模型主要包括输入层、隐含层和输出层，RBF 模型的
表达式为：
F ()X =åi =1q
w i φi
X -c i （7）
φ为可逆高斯函数：
φ()t =exp æèçö
ø
÷
-t 2δ2（8）
式中：q 为隐含层节点个数即样本容量；w i 为隐藏层各节点到输出层F ()X 的权因子；c i 为第i 个径向基函数的中心值；δ为径向基函数的变量，决定了该径向基函数围绕中心点的宽度。
模型的建立通过Matlab2017b 的RBF 工具箱来实现，具体函数的调用格式为：
net =newrb ()
P ，T ，goal ，spread ，MN ，DF
（9）
式中：P 、T 分别为输入和输出变量；goal 为均方误差，取goal =0.001，避免过度拟合；spreads 为RBF 函数的宽度，通过隐含层神经元的中心距离确定；MN 为ANN 学习过程中神经元最大数目；
DF为2次显示之间添加神经元的个数。
仿真函数调用为：
Y=sim()
net，X（10）式中：Y为预测目标值；X为输入测试矩阵，net为式（9）中已经学习好的网络。
（4）训练集、验证集及测试集。选择1981—2018年1—12月444组修正后的历史气象月数据（T max、T min、RH和Ws）和以从中国气象科学共享服务网（http：//v）获取的研究区逐日ET0累加得到的月ET0作为训练集；选择2013—2015年36组修正后的历史天气预报月数据（根据天气预报日值累加所得）和月ET0数据作为验证集；选择2016—2018年36组修正后的历史天气预报月数据（根据天气预报日值累加所得）作为测试集。
2.3模型评价对于使用RBF神经网络模型预测值与PM计算标准值之间的相互比较，使用了5个评价指标，包括决定系数（R2）、纳什系数（NSE）、均方根误差（RMSE）、平均绝对误差（MAE）和相对误差（RE）来评估模型的预测性能［25］，计算如下：
R2=
å
i=1
n
()
x
i
-xˉ2()
y
i
-yˉ2
å
i=1
n
()
x
i
-xˉ2å
i=1
n
()
y
i
-yˉ2
（11）
NSE=1-å
i=1
n
()
y
i
-x
i
2
å
i=1
n
()
y
i
-yˉ2
（12）
RMSE=（13）
MAE=å
i=1
n|
|x
i
-y
i
n（14）
RE=x
i
-y
i
y
i
×100%（15）
式中：x i为ET0预报值（利用修正后的气象数据结合径向基神经网络预测所得）；y i为ET0标准值（从中国气象科学共享服务网（http：//v）获取的研究区逐日ET0）；y为标准值均值；i为预报样本序列，i=1，2，…，x；n为预报值的样本数。
3结果与分析
3.1数据集及先验密度图1和表2为华北平原农田下垫面1981—2018年气象因子和ET0的数据系列，气象因子包括最高温度（T max）、最低温度（T min）、相对湿度（RH）、平均风速（Ws）和ET0。T max、T min、RH、Ws和ET0多年平均值分别为15.15℃、6.74℃、4
4.20%、1.89m/s和1063mm，多年偏差系数分别为-0.36、-0.34、0.65、-0.36和-0.68，多年的气象因子和ET0都在均值周围有微小的波动，标准偏差分别为0.49℃、0.53℃、3%、0.16m/s和0.1mm。
图2为气象因子的正态概率分布。由图2可知，各气象因子均近似服从正态分布，通过SPSS软件分析所知各因子的正态性假设均通过5%的显著性水平。ET0影响因子的先验密度和标准方差见表3，由表3可以看出，各月Ws的分布差异性较大，变差系数（C V均值为0.13）比T max（C V均值为0.11）、T min （C V均值为0.06）和RH（C V均值为0.07）的大，表明Ws相对其他气象因子的不确定性较强。
3.2气象因子不确定性识别由于本研究所选气象因子为4个变量，即T max、T min、RH、Ws，因此C0矩阵对角线元素取为各变量的先验方差。设定AM-MCMC算法的初始抽样次数为2000，运行次数38次，每次采样10000个，这样共抽取有效模拟样本数为304000个，用于每月气象因子后验分布的统
计分析。气象因子在迭代2000次后已经收敛，所抽取的样本已具有其总体的统计特征。
分别以3月（作物生长早期）、7月（作物生长中期）和10月（作物生长后期）典型月为例，所抽取各影响因子的后验密度直方图见图3—图5。T max 、T min 、RH 、Ws 的后验密度分别为：最高温度x 1~N （10.46，0.982
），最低气温x 2~N （1.49，0.562
），相对湿度x 3~N （34.62，0.022），平均风速x 4~N （2.27，
0.132）；最高温度x 1~N （25.73，0.522），最低气温x 2~N （18.54，0.352
），相对湿度x 3~N （58.97，0.012
），平均风速x 4~N （1.65，0.062
）；最高温度x 1~N （16.03，0.682），最低气温x 2~N （7.03，0.412），相对湿度
x 3~N （49.69，0.012），平均风速x 4~N （1.60，0.072）。对于Ws 先验密度的修正较T max 、T min 和RH 大，在相同的显著性水平条件下，各影响因子的置信区间不同，其中Ws 比T max 、T min 和RH 的范围大，说明Ws 具有较强的不确定性，主要是因为Ws 受到外界的干扰较为敏感。因此，后验信息对Ws
的修正程
图1
研究区1981—2018年气象参数及ET 0资料系列
因子T max /℃T min /℃RH /%
Ws /(m/s)ET 0/(mm)
均值15.166.74
44.201.892.30中值15.246.7544
1.90
2.32RMSE 0.490.533
0.160.10变差系数-0.36-0.340.65
-0.36-0.68表2
气象参数和ET 0
图2各气象因子的正态概率图
T max /℃
T min /℃
RH /%
Ws /（m/s ）（a ）T max 正太概率（b ）T min 正太概率
（c ）RH 正太概率
（d ）Ws 正太概率
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议