基于不同机器学习模型的川中丘陵区参考作物蒸散量模拟

黄滟淳，崔宁博,2,3，陈宣全，徐浩若，张艺璇

(1. 四川大学 水力学与山区河流开发保护国家重点实验室 水利水电学院，成都 610065；2. 南方丘区节水农业研究 四川省重点实验室，成都 610066；3. 西北农林科技大学旱区农业水土工程教育部重点实验室，陕西 杨凌 712100)

0 引 言

参考作物蒸散量(ET0)是指高0.12 m、冠层蒸散阻力为70 s/m、反照率为0.23的假想作物的蒸散速率[1]，ET0可以表征某一地区的大气蒸散能力，是计算作物需水量的基础，同时也是规划和设计农田水利工程的重要资料[2]，对作物需水量预测、水资源优化管理具有重大意义。

目前确定ET0的方法有实测法、公式法和数值模拟法，其中实测法操作繁杂，限制条件较多，不易推广[3]；公式法中的标准模型是由联合国粮食及农业组织(FAO)于1998年提出的Penman-Monteith (P-M) 模型[1]，该模型较全面地诠释了蒸散发过程的发生机制，但是P-M模型考虑的因素很多，对于不能获取完整气象数据的地区并不适用，因此，基于较少气象参数输入的简化预报模型尤为重要，目前已有60余种简化模型被相继提出，包括温度法中的Hargreaves-Samani模型[4,5]，辐射法中的Priestly-Taylor模型[6]、Irmak-Allen模型[7]，综合法中的Penman-Van Bavel模型[8]等；数值模拟法将既有气象数据作为样本，模型经过学习和训练，找到非线性关系的最优拟合，该方法高效、可移植性强，在区域范围广、气候条件复杂的地区[3]优势显著。

随着计算机技术及机器模型的发展，越来越多的数值模拟模型被用于ET0的预报。冯禹[9,10]等对中国西南部湿润地区的ET0预报研究表明，基于温度和地外辐射数据的随机森林(RF)和广义回归神经网络(GRNN)模型具有良好的预报效果[9]，极限学习机[10](ELM)的预报效果优于经验模型。Tabari[11]等将支持向量机(SVM)和自适应神经模糊推理系统(ANFIS)应用于伊朗半干旱区域的ET0预报，预报效果优于经验模型。Yassin[12]等利用人工神经网络(ANNs)和基因表达编程(GEP)对干旱条件下的ET0预报进行比较分析，结果表明ANNs模型预报精度高于GEP。

模型树常被用于数据挖掘和机器学习中，通常采用离散的标签进行决策分类，相关学者[13,14]对其进行改进，用分段线性函数替代离散分割，为模型树在连续型数据回归问题的应用提供了技术支撑。M5回归树不同于传统的神经网络，其精巧的分段式结构和灵活的模型尺度，使得它拓扑结构直观、收敛速度快，且能够得到模拟映射的线性表达，泛化能力更强[15]。目前M5回归树已成功应用于降雨—径流模拟中，在流量及区域降雨的预测问题上表现良好，但还没有将其利用于ET0预报的相关研究或报道。

本文以P-M模型作为标准，将M5回归树(M5-RT)与改进的传统神经网络BPNN、GRNN应用于川中丘陵区7个代表性站点的ET0预报，分析不同模型的ET0预报精度和泛化能力，并将其与精度较高的经验模型进行对比，为川中丘陵区的ET0简化模拟提供新思路。

1 材料与方法

1.1 基本资料

川中丘陵区位于四川省东部，属于典型的方山丘陵区，面积约8.4 万km2，占四川盆地约50%，海拔为250～600 m，地表起伏多但高差变化不大[16]。川中丘陵区是四川省农业生产的重要区域，主要经济作物有蚕桑、甘蔗、棉花等，区域内气候灾害频发且持续时间较长[17]，严重影响区域内农业发展，实现精准农业管理至关重要。考虑区域内地理、气候因素，选取7个代表性站点，具体分布见图1。

图1 川中丘陵区气象站点分布图Fig.1 Distribution of meteorological stations

1.2 数据获取及处理

本文所用的川中丘陵区1961-2016年逐日气象数据均来自国家气象信息中心(https:∥data.cma.cn/)，包括最高气温(Tmax)、最低气温(Tmin)、平均气温(Tmean)、日照时长(Tsun)、相对湿度RH(relative humidity)、10 m处风速(u10)等。针对达县、广元站缺失数据(约占0.15%)，采用线性插值[18]补全。并将1961-2016年的逐日气象数据按7：3的比例分为两部分，分别作为训练集(1961-2000年)和预测集(2001-2016年)。

选取FAO P-M模型[1]的计算值作为ET0标准值，参考P-M模型，风速项的条件为距地面2 m，故由风廓线[1]关系推出式(1)进行转化。

(1)

式中：z为地面至测点的垂直距离，m；uz为高度z处的风速，m/s。

1.3 ET0计算模型

P-M模型综合考虑了各种气象因素，基于水汽扩散理论推导，具有很强的普适性。另外，参考赵璐等[19]对川中丘陵区不同计算方法的比较和改进，以及冯禹等[20]对机器学习模型与Hargreaves模型在四川盆地的对比研究，选取3种在川中丘陵区精度较高的简化物理模型进行对比，分别是基于温度、日照时长和大气顶层辐射的Jensen-Haise模型[21]，李晨等[6]基于贝叶斯原理改进的Hargreaves-Li模型，基于温度和净辐射的Irmak-Allen模型[8]，计算公式见表1。

表1 ET0物理计算模型Tab.1 Calculation models of reference crop evapotranspiration

注：Δ为饱和水气压—温度曲线的斜率，kPa/℃；γ为湿度计常数，kPa/℃；Rn为净辐射，MJ/(m2·d)；Ra为大气顶层辐射，MJ/(m2·d)；ea为实际水气压，kPa；es为饱和水气压，kPa；G为土壤热通量，MJ/(m2·d)；N为白昼时长，h；λ为水的汽化潜热，取25 ℃下的标准值2.444 MJ/kg参数的具体计算见文献[1]，特别地，本文获得的气象数据以天为周期，可忽略土壤热通量的影响[1]，故取G=0。

1.4 机器学习算法

1.4.1 M5回归树模型M5-RT

模型树在数据挖掘和机器学习中应用广泛，其扩展结构如图2所示，像一棵倒置的树，叶片代表不同的分类或决策。传统的模型树采用离散的标签进行分类，Wang和Witten[13]对M5回归树的叶进行了改进，利用若干个线性函数替代离散标签，从而实现连续类数据的学习和预测回归，相关技术由Jekabsons等[14]于2010年继续完善，并利用Matlab工具箱研发相关程序(http:∥www.cs.rtu.lv/jekabsons/)。

M5回归树是一种多元的分段式线性回归模型，节点选择时采用贪心算法，理论上可使用足够多的节点完成任何复杂的非线性回归，实际建模中为提升收敛速度，会通过后续的剪枝操作[15]对树模型进行精简，分段函数的特性，使得M5回归树比其他数据模型更为灵活，且能得到映射关系的显式表达，目前已成功应用于降雨—径流模拟中。Goyal等[22]将M5回归树用于Pichola湖流域降雨量和水量预测，Bhattacharya等[23]将M5回归树与神经网络应用于河流流量预测，取得了优于传统模型的结果。

图2 M5回归树结构Fig.2 Structure diagram of M5 regression tree

1.4.2 交叉验证改良的广义回归神经模型CV-GRNN

广义回归神经网络(GRNN)是由Specht[24]提出的一种变形的径向基神经网络，以样本数据为后验条件，执行Parzen非参数估计[25]，具有很强的非线性映射能力。Ladlani等[26]将GRNN用于阿尔及利亚的ET0预报，并证明其精度优于传统的径向基神经网络。冯禹等[27]将基于温度资料的GRNN模型用于四川盆地的ET0预报，精度优于物理模型。GRNN模型由4层网络组成，如图3，依次为输入层、模式层、求和层和输出层。

图3 GRNN模型拓扑结构Fig.3 Topological structure of general regression neural network

(1)输入层为学习样本的直接输入，其神经元个数等于学习样本中自变量X=[x1,x2,…,xk]的维度k，深度为样本数量。

(2)模式层与输入层间通过权重wp连接，模式层的神经元个数等于样本数量n，各神经元的传递函数为，

(2)

式中：Xi为第i个神经元的学习样本；σ称为光滑因子，是高斯分布的标准差，其取值影响GRNN网络的性能，需要优化。

(3)求和层使用两种计算法则对模式层的神经元进行处理。第一种对模式层所有神经元进行算术求和，连接权值为1，其求和函数如式(3)；第二种对模式层的神经元进行加权求和，本文ET0预报的输出变量为一维，设样本输出值为Y，第i个神经元与第i个输出变量的连接权值就是Yi，其求和函数如式(4)。

(3)

(4)

(4)输出层的神经元个数等于输出样本的维度，此处为1，模型输出的估计值是两类求和层的比值。

(5)

本文选取川中丘陵区7个代表性站点的逐日气象资料建立GRNN模型，代码见文献[28]，同时采用交叉验证(Cross-validation)算法对光滑因子σ及径向基扩展速度spread进行优化，将训练样本随机分为5个子样本，4对1重复验证5次，在参数的预设范围内循环求解，确定优化参数建立CV-GRNN模型。

1.4.3 双隐藏层优化的BPNN模型H-BPNN

1988年，Rumelhart[29]等提出了误差反向传播算法，即BP算法，应用于反向传播神经网络(BPNN)。BPNN属于前馈型神经网络的一种，会根据理想输出结果和实际预测结果的差值对模型中的传递权值进行修正，反复训练直至误差达到设定阈值以下或者迭代次数达到最大值[30]。

本文用动量批梯度下降函数(traingdm)替代默认的模型训练函数trainlm，动量项的引入提高了收敛速度，且能有效避免网络训练时[31]陷入局部极小值震荡。其权重更新规则改进为，

(6)

式中：w(t)为t时刻的权重值；η为学习速率；α为动量因子；E为网络训练的误差。

动量项的添加能加强权重变化方向的记忆效应，使得学习速率加快并有效拜托局部极小值区域。

隐含层的层数和节点个数对模型的预测性能和泛化能力有很大影响[32]，过拟合或欠拟合都会产生较大误差，因此，需要预试验确定模型的最佳结构[33]。本文采用双隐含层的BPNN网络，节点数取值范围初设为1～20，基于广元站的气象数据进行优化求解，代码基于Matlab 2018a，具体见文献[28]，以均方根误差为考察目标，设置最大迭代次数为4000，求解最优隐含层节点数，预试验结果见图4。

图4 BPNN网络性能与隐含层节点数的关系Fig.4 Performance of back propagation neural network in different hidden neurons

由图4，模型误差与隐含层节点数量n大致呈抛物线关系[34]，综合考察迭代次数及网络精度，本文建立隐含层节点数为7的双隐含层H-BPNN模型，模型结构见图5。

图5 H-BPNN模型结构图Fig.5 Structure diagram of H-BPNN

1.5 模型输入层

刘昌明等[35]对全国主要流域的地表潜在蒸散量敏感性研究表明，风速与日照时数的下降是造成ET0减小的主要原因。另一方面，太阳辐射是引起白昼温差的主要原因，Hargreaves[36]和Samani等[37]利用温差和大气顶层辐射(Extraterrestrial Radiation，Ra)对地表净辐射进行推演，弥补Tsun的缺失。因此，本文在输入项的气象因子组合中优先考虑温度和u2，并比较Tsun和Ra的作用，建立如表2的4种输入组合。

式(7)中Ra是与站点纬度(Latitude，单位rad)和日序数(J，每年1月1日起从1开始循环)相关的函数，具体计算方法见文献[1]。

Ra=F(Latitude,J)

(7)

1.6 模型评价因子

本文根据模型预测结果和P-M模型计算结果计算4个统计参数对模型预报精度进行评价，分别是纳什效率系数NSE，均方根误差RMSE，平均相对误差MRE和决定系数R2，并由对应的综合指标GPI[38](global performance indicator)对模型性能进行排序。MRE和RMSE越接近0，R2和NSE越接近1，表明模型预报精度越高，故GPI计算公式如式(8)，GPI越大，模型精度越高。

表2 模型输入项参数组合Tab.2 Input combinations and parameters

注：InputI(I=1, 2, 3, 4 )表示4种不同的输入组合，下同；按输入项维度递减排列，其中大气顶层辐射Ra可由日序数和纬度直接计算[1]，非观测量。

(8)

式中：Ckm为模型m的第k个参数；Midk为所有模型中参数k的中位数；GPIm为模型m的综合得分。

对于R2和NSE，ηk=1，对于MRE和RMSE，ηk=-1。

2 结果分析

2.1 各模型ET0日值预报结果及精度分析

根据1.6的模型评价方法对各模型的ET0预报结果进行分析，计算川中丘陵区不同模型ET0预报的统计参数，7个站点的参数范围见表3。由表3可知，输入组合为Input 1(Ra、Tmax、Tmin、RH、Tsun和u2)时，M5回归树(M5-RT)表现出极高的精度，NSE、RMSE、MRE和R2的变化范围分别为0.998 0～0.999 2、0.037 4～0.058 0 mm/d、0.88%～1.62% 和0.998 2～0.999 2，GPI排名第1，CV-GRNN和H-BPNN模型也有较好的预报精度，其NSE和R2均大于0.940 2。在3种数值模拟模型中，M5-RT最优，H-BPNN次之，CV-GRNN相对较差，但3种模型的NSE均大于0.843 8、R2均大于0.893 3、RMSE均小于0.5 mm/d、MRE均小于13.59%，都能在一定程度上呈现气象参数与ET0间复杂非线性映射关系，且模拟精度均高于所选经验模型，尤其是M5回归树，在4种输入项组合下都有最优的表现，从图6可以看出M5-RT模拟预报误差的异常波动最小，在川中丘陵区的ET0预报研究中有很高的实用价值。

表3 川中丘陵区不同模型的逐日ET0预报精度Tab.3 Daily statistic performance evaluation results for different models in hilly area of central Sichuan

注：XXXI(I=1, 2, 3, 4 )分别表示4种不同输入下的XXX模型，下同；**表示在1%的水平上显著相关；“x±y”表示参数范围是[x-y,x+y],下同。

图6 2001-2016年川中丘陵区不同模型ET0日值预报误差箱线图Fig.6 Box-plot of ET0 forecasting accuracy for different models in hilly area of central Sichuan from 2001 to 2016

Input 1包含了最多的气象因子，从图6的整体变化趋势不难发现，Input 1下模型具有精度最高，其他3种输入组合下模型的精度随输入参数减少而降低，但M5-RT、CV-GRNN和H-BPNN的ET0预报精度仍处于较高水平。在缺省气象资料的情况下，模型的可靠性以及实用性评价需要考量减少输入项对模型精度产生影响的程度，只要总体精度达标，相比全输入模式下轻微的精度降低是可以被接受的，这一点在Input 2和Input 3上表现极为突出。

Input 2相比于Input 1去掉了Tsun和RH，模型预报精度降低幅度低于0.4%，CV-GRNN模型反而略有升高，其精度波动可以归于模型训练的随机性误差，即Tsun和RH的缺失对数值模拟模型几乎没有影响。正如Hargreaves[36]和Samani[37]等研究发现，利用温差和Ra对地表净辐射进行推演，能弥补Tsun在模型中的作用，刘昌明[36]等的研究也表明风速对ET0的影响力远大于RH。另一方面，Input 3在Input 2的基础上(去掉了Ra)约有5%的精度降低，进一步体现了大气顶层辐射Ra的重要性。

Input 4(Ra和u2)输入下，M5和H-BPNN均表现出良好的预报精度，NSE和R2均大于0.9、RMSE均小于0.4 mm/d、MRE均小于10%，模型精度高于其他在川中丘陵区适应能力较好的物理模型。Input 4保留了最为重要的气象因子，其中Ra是温差与辐射的体现，风速与其他气象参数(例如温度和辐射)的关系比较微弱，这种不可替代性使得对模型性能有着重要的影响，这一点与张皓杰[39]和冯禹[20]的研究结论一致，在地理环境、气候条件各异的不同流域内，风速对地表蒸散发的影响均不可忽视[35]。

2.2 M5-RT、CV-GRNN、H-BPNN模型可移植性分析

Input 2(Ra、Tmin、Tsun和u2)在较少输入的情况下取得了较高的ET0预报精度，输入项中的Ra使得模型学习到更多关于站点地理位置的信息。本文为检验模型在川中丘陵区不同地区间的适应能力，在Input 2的基础上，将7个代表性站点按照纬度分为3类(偏北部的广元、绵阳、巴中，偏中部的乐山、达县、遂宁，偏南部的叙永)，以叙永站2001-2016年逐日气象资料及P-M模型计算的ET0作为输出，分别从北、中各3个站点中随机抽取20年日值数据进行组合，即建立“北-南”、“中-南”两类训练-预测组合，训练集数据量为60年，预测集数据量为16年，按照1.6节的模型评价方法计算相关参数，并将预报结果与P-M模型标准值进行比较，图7为P-M模型计算结果与模型预测输出的散点图。

可移植性分析结果表明，训练、预测站点交叉组合下，3种模型的预报精度都较高，NSE和R2均大于0.9。同时，其精度与相同输入下的独立站点模拟精度相比都有一定程度的降低，特别是均方根误差RMSE，增大了近50%～100%。

由图7可知，M5-RT模型两种站点交叉组合下模型的ET0预报能力接近，R2为0.955 7～0.956 0，NSE为0.953 3～0.955 2，模拟结果平均相对误差约为7.5%～7.6%，精度降低约5%，RMSE为0.29～0. 31 mm/d，ET0模拟输出较为稳定，精度波动较小，整体效果好；H-BPNN模型和CV-GRNN模型出现明显的截断误差，当ET0大于6 mm/d时，预报值普遍过小，具有一定局限性，模型可靠性降低；CV-GRNN模型在整体上相关性较好，但预测值与标准值的趋势线斜率较小(0.819和0.913)，ET0预报值较低，效果相对较差。对比发现M5-RT模型泛化能力最强，在气候条件接近的地区间预报能力优且稳定，可作为川中丘陵区预报ET0的理想模型。

3 结 论

(1)日尺度误差分析表明，4类输入项组合下，M5回归树模型最优，H-BPNN模型次之，CV-GRNN模型相对较差，且对比不同输入组合下模型精度发现，Tsun和RH对ET0预报的影响很小，而Ra和u2对模型的精度影响较大。另外，3类机器学习模型的ET0预报精度都较高， 其中R2为0.89～0.99，NSE为0.84～0.99，MRE为1.25%～12.45%，RMSE为0.05～0.49 mm/d，优于经验模型Jensen-Haise、Hargreaves-Li和Irmak-Allen。

(2)当输入项气象因子减少时，3类机器学习模型的预报精度均有一定程度的降低，但M5-RT模型表现最优，特别是气象因子数量最少的Input 4输入下，M5-RT仍表现出良好的预报效果，其R2为0.917 9～0.932 3，NSE为0.917 8～0.931 0，RMSE小于0.386 7 mm/d，MRE小于9.48%，在川中丘陵区的适应性最广。M5-RT3和M5-RT4均可作为气象数据缺失下川中丘陵区ET0预报的推荐模型。

(3)通过Input 2下各模型的可移植性分析发现，M5回归树泛化能力最强，预报能力优且稳定。整体上，ET0大于6 mm/d时，H-BPNN模型和CV-GRNN模型出现明显的截断误差，且CV-GRNN模型的ET0预报值普遍偏小， M5回归树更适合作为川中丘陵区预报ET0的推荐模型。M5回归树良好的泛化能力使得当某地区气象数据在时间尺度上残缺时，可以适当采用附近地区的气象数据用于M5回归树的训练，仍可获得精度较高的ET0预报结果。

M5-RT、CV-GRNN和H-BPNN等数值模拟模型能辨析出气象因子间固有的隐含关系，完成部分气象因子的替代和剔除，从而减少观测层面的工作量，这正是一般经验模型所不具备的能力。为进一步分析各个气象因子对ET0预报精度的影响，应设置更多的输入组合进行对比分析。

□