基于NARX模型的参考作物蒸散发预测

武剑飞，康银红，宋鑫，梁友鹏

(四川农业大学水利水电学院，四川 雅安 625014)

作物蒸散量是参与大气水文循环的重要因子，同时也是农业灌溉需水量制定的重要参数.参考作物蒸散量(reference crop evapotranspiration，ET0)是国际上计算作物蒸散量的通用参数，准确计算及预报ET0具有重要意义.目前关于ET0的确定主要分为3个方面：试验测定、公式计算、数值模拟.运用蒸渗仪法测定的ET0精度好，但是仪器费用高且耗时长，只适用于一些特定研究，不适宜大范围推广[1-2].基于气象参数建立的ET0计算公式是许多学者常用的方法，目前已发展出50多种计算方法.在1998年，联合国粮农组织(food and agriculture organization，FAO)把Penman-Monteith(PM)公式作为计算ET0的唯一标准方法，并且在全世界范围内都得到了很好应用，但是该公式需要气象参数多，计算复杂，对于一些气象资料缺失的地区不具有适用性[3].因此，一些学者研究出了简化的ET0计算模型，例如温度法中的Hargreaves-Samani模型[4]，辐射法中的Makkink模型[5]、Priestley-Taylor模型[6]，经验法中的Irmak-Allen模型[7]等.这些简化模型虽然能够在气象资料缺失的情况下计算ET0，但是这些模型都是在特定的地区发展而来的，对于不同的环境需要筛选不同的计算模型.

随着计算机技术的发展以及仿真模拟的应用，基于机器学习下的ET0模拟计算以及预测成为国内外学者的研究热点.ANTONOPOULOS等[8]、JAHANBANI等[9]运用人工神经网络(artificial neural network，ANN)模拟了ET0，结果表明ANN模拟的精度比简化模型计算的精度高；ADAMALA[10]对印度不同地区建立广义小波神经网络(wavelet-neural network，WNN)以模拟ET0，结果发现WNN模型比传统的ANN模型、线性回归模型和Hargreaves模型的模拟误差小，具有良好的应用能力.崔宁博等[11]基于思维进化算法(mind evolutionary algorithm)优化的误差反向传波神经网络(back propagation neural network，BPNN)模型预测了西北干旱地区的ET0，发现该模型在相同气象参数输入的情况下，模拟精度远高于Hargreaves-Samani模型、Irmak模型和48-PM模型；冯禹等[12]基于只有温度资料建立广义回归神经网络(generalized regression neural network，GRNN)，通过与Hargreaves模型、改进的Hargreaves模型比较，发现该模型在时空尺度上能够很好地预测四川盆地的ET0；王升等[13]基于基因编程表达式算法(gene expression programming，GEP)以及径向基神经网络算法(radial basis function-artificial neural network，RBFNN)建立了ET0预报模型，发现在气象参数输入相同的情况下GEP模型比RBFNN模型的计算精度更高.

动态神经网络是一种具有动态反馈性的神经网络，能够有效实现复杂的动态系统建模.常见的动态神经网络有带外部输入非线性自回归(nonlinear autoregressive models with exogenous inputs，NARX)神经网络、非线性自回归(nonlinear autoregressive，NAR)神经网络等，目前广泛应用于电力负荷预测[14]、锅炉温度预测[15]等方面，对于ET0的预测还未有报道.

1 材料与方法

1.1 数据来源

选取四川西南山地的攀枝花站点作为研究站点，气象数据包括日最高温(Tmax)、日最低温(Tmin)、日照时数(n)、相对湿度(RH)和10 m高度处的风速(采用FAO风速轮廓线关系[3]换算为2 m高度处的风速u2)的2010—2018年的逐日气象资料.其中，2010—2017年的数据作为训练集，2018年的数据作为测试集；数据来自国家气象信息中心，且经过较好的质量控制.

1.2 参考作物计算模型

选取Penman-Monteith模型计算的ET0作为标准值，选用Hargreaves-Samani 模型、Irmak-Allen模型、Makkink模型以及Priestley-Taylor模型进行对比分析，具体计算公式见表1.

表1中ET0为参考作物蒸发蒸腾量，mm/d；Δ为饱和水气压曲线斜率，kPa/℃；Rn为作物表面净辐射，MJ/(m2·d)；Ra为地外辐射，MJ/(m2·d)；Rs为太阳辐射，MJ/(m2·d)；G为土壤热通量，MJ/(m2·d)，以天计算的土壤热通量为0；γ为湿度常数，kPa/℃；Tmean为日平均温度，℃；es为饱和水气压，kPa；ea为实际水气压，kPa；λ为气化潜热，2.45 MJ/kg.

表1 参考作物蒸散量计算公式

1.3 NARX模型

NARX是一种非线性动态循环网络模型，由输入层、隐含层、输出层及输入和输出延迟构成[16]，模型结构见图1.

图1 NARX模型结构图

模型带有外部输入变量，能够接受来自输出神经元的反馈，反映过去输入与输出的历史信息，是一种具有记忆功能的动态神经网络[17-18].NARX同时引进外部输入和输出两种时间序列，运用输出y(t)和另一个时间序列x(t)的历史值以模拟预测未来时间内y(t)，模型数学表达形式为

y=f[x(t-1),…,x(t-m),y(t-1),y(t-2),…,(t-n)]，

(1)

式中：m，n分别为时间序列x(t)和y(t)的延迟阶数；t为时间.

从实现方式来说，耗材信息化建设的本质是信息技术与医疗行业的有机融合，其主导因素是信息技术的运用。但随着医院管理需求的变化，其对应的应用需求也要进行适时调整，这就决定了在信息化过程中，信息技术与信息管理是既独立又密切相连的特殊关系，充分保障信息化技术手段平稳发挥作用是医院管理部门的重要职责，通过技术手段实现医院运营管理一体化，能更好地发挥信息资源的高效利用［3］。

NARX模型的输入及输出延时阶数设为1，隐含层的个数为10，训练函数为“trainlm”，利用Matlab软件进行模拟预测.

1.4 NAR模型

NAR是一种基于时间序列且具有反馈和记忆功能的动态神经网络模型，模型的输入和输出并不是静态的映射方式，而是每个输出基于当前系统的动态结果综合得到[15].模型由时间序列y(t)的历史值预测未来时段内的y(t)值，数学表达式为

y(t)=f[y(t-1),…,y(t-n)].

(2)

模型由输入层、隐含层、输出层和滞后层组成.NAR模型的隐含层个数默认取10，取滞后阶数默认值1∶2，训练函数为“trainlm”.

1.5 模型评价方法

选取均方根误差(RMSE)、平均绝对偏差(MAE)、平均偏差(MBE)、纳什效率系数(NSE)、综合评价指标(GPI)评价各模型的精度，计算公式为

(3)

(4)

(5)

(6)

(7)

RMSE，MAE和MBE越接近0，表明模型的精度越好；NSE越接近1，表明模型效率高.GPI采用排名形式呈现，排名越靠前，模型的性能越好.

2 结果分析

2.1 不同气象因子组合下的NARX模型预测

基于Tmax，Tmin，RH，u2，n这5种与ET0相关的气象因子进行随机组合，建立11组NARX模型，分析不同气象因子组合下的模型模拟效果，具体结果见表2.

表2 不同气象因子组合下NARX模型的评价指标

当输入5个气象因子时，模型NARX-1的RMSE，MAE，MBE分别为0.425 mm/d，0.320 mm/d，0.069 mm/d，NES为0.920，GPI排名第11，模型模拟效果在所有组合中效果最差，这与冯禹等[19]、魏俊等[20]对ELM模型研究结果(输入5个气象因子的模型模拟效果最优)有所差异，这可能是由于某个因子对ET0的影响已包含在另一因子中(例如风速可通过温度反映对ET0的影响)，而该因子能够充分反映对ET0的影响，NARX模型包含了过多的冗余信息.过多的输入维度使得NARX模型趋于复杂，泛化能力下降，模拟精度较低.

当输入4个气象因子时，模型NARX精度相比模型NARX-1有所提高，但是NARX-2和NARX-3之间差异比较大.NARX-2(缺失RH)的RMSE，MBE分别为0.363 mm/d和0.035 mm/d，比NARX-3(缺失n)的RMSE，MBE分别减少了0.082 mm/d和0.032 mm/d，而NSE增幅为3.2%，两者GPI排名分别为第8、第9.说明日照时数对ET0的影响比相对湿度大，这与冯禹等[21]对四川地区蒸散发研究中得出日照时数是四川地区的主要影响因子的结论一致.

当输入3个气象因子时，模型NARX精度差异明显.NARX-5(缺失RH，n)精度最高，它与NARX-4(缺失RH，u2)的RMSE均低于0.370 mm/d，NSE大于0.940，而NARX-6(缺失n，u2)的RMSE高于0.400 mm/d，NSE低于0.930，这说明n和u2对ET0有着重要的影响.除了NARX-6对ET0低估外(MBE<0),其他模型均高估了ET0.对比模型NARX-5(输入Tmax，Tmin，u2)和NARX-3(输入Tmax，Tmin，RH，u2)可以发现，在缺失RH时，模型NARX-5的RMSE下降幅度为25%，NSE增加了4.1%，GPI排名从第9升至第2，说明相对湿度对该地区ET0影响较小.

当输入2个气象因子时，NARX-7的RMSE比NARX-8低0.050 mm/d，MAE减少0.063 mm/d，MBE降低58%，NSE增加1.7%，GPI排名第4，精度明显优于NARX-8.对比NARX-8(缺失Tmax，Tmin，RH)和NARX-4(缺失RH，u2)和NARX-6(缺失n，u2)，减少输入参数后模型维度降低，精度并不一定提高.这主要与气象因子对ET0的贡献程度有关.

当输入1个气象因子时，模型NARX-11(输入RH)精度最低.NARX-9(输入u2)与NARX-10(输入n)在所有模型中精度均较高.其中，模型NARX-9的GPI排名第1，NSE为0.964，而RMSE和MAE均低于0.300 mm/d.进一步说明日照时数在该地区对ET0影响较大.而基于较少气象因子建立的BPNN模型[22]和WNN模型[23]的模拟精度较差，这说明NARX模型泛化能力更强，模拟效果更好.

2.2 模型NARX与NAR精度比较

为了比较模型NAR与NARX的精度，选取文中精度较差的模型NARX-1，分别绘制NARX-1和NAR的预测值与PM计算值的线性拟合图，如图2所示.由图可知，模型NARX-1的预测值与PM计算值拟合效果最好，趋势线斜率趋近于1，R2>0.900，与PM计算值有很好的相关性.模型NAR的预测值分布较为散乱，斜率为0.706 3，R2远小于0.700，模拟误差比较大.因此，模型NARX比NAR的模拟精度更高.

2.3 NARX模型与其他模型的精度比较

以PM模型计算的ET0作为标准值，比较基于Tmax和Tmin的NARX-7和Hargreaves-Samani模型以及基于Tmax，Tmin和n的NARX-4，Irmak-Allen，Makkink，Priestley-Taylor模型的精度.模拟精度如图3所示.各模型与PM模型的拟合效果之间差异具有统计学意义.模型Irmak-Allen和Makkink的ET0整体在1∶1线以下，其中Irmak-Allen拟合效果最差(拟合斜率为0.5168，R2为0.771)；模型Hargreaves-Samani和Priestley-Taylor的ET0在1∶1线附近分布较为散乱，R2在0.800左右；模型NARX-7和NARX-4的ET0均匀分布在1∶1线附近，且R2在0.950以上，拟合效果最优.

图2 NARX-1，NAR模型与PM模型的线性拟合

图3 PM模型与其他模型的线性拟合

进一步对各模型的精度进行统计分析.图4为NARX模型与其他计算模型精度比较，图中σ为精度比较误差.由图可以直观看出各模型误差情况.模型Irmak-Allen在所有模型中误差较高，RMSE大于1.000 mm/d，NSE小于0.200，精度最差；在输入气象因子相同的情况下，NARX模型的误差值均小于简化模型的误差值，模型Irmak-Allen和Makkink的MBE均小于-0.500 mm/d，表明计算的ET0值偏低较多.对比 NARX-7(输入Tmax，Tmin)和Hargreaves-Samani(输入Tmax，Tmin)模型，NARX-7的RMSE，MAE和MBE均低于0.500 mm/d，NSE接近于1，而Hargreaves-Samani的各项误差远大于0.600 mm/d，NSE低于0.700，比NARX-7的精度低很多.相比较输入相同气象因子的模型NARX-4，Irmak-Allen，Makkink和Priestley-Taylor，除了模型NARX-4和Priestley-Taylor的MBE比较接近，其他情况下模型NARX-4的MAE和RMSE比其他3个模型误差小得多，NARX-4的NSE在0.900以上，其他模型的NSE均低于0.700，说明NARX-4模拟精度远比其他3个模型精度高.

以上分析表明，在日照时数、湿度、风速气象资料缺失的情况下，NARX模型的模拟精度比简化模型更具有优势，可作为该地区预测ET0的优良模型.

图4 NARX模型与其他计算模型精度比较

3 讨 论

基于气象因子建立的Hargreaves-Samani和Makkink等各种ET0简化模型虽然能够很好地解决某些气象资料缺失下的ET0获取问题，但是这些简化模型普遍存在精度较差且模型具有区域性问题，需要进行筛选且修正后才能应用到该地区.随着机器模拟技术的发展，基于各种机器学习条件下的ET0模型不需要具体的公式，仅建立气象因子与ET0之间的关系进行模拟，极大地提高了ET0的计算精度.NARX模型与静态神经网络结构相比，输入层增加了历史值y(t)，每次y(t)反馈到输入层，从而达到长期记忆的能力，而且对y(t)影响因素较大的输入变量能够很好反映变量间的非线性关系.

研究表明，在气象因子缺失条件相同的情况下，NARX模型比简化的ET0计算模型的精度更高.通过建立不同气象因子下的NARX模型发现，基于单个气象因子的NARX模型的精度(RMSE低于0.400 mm/d，NSE在0.940以上)高于其他一些机器学习模型，例如张皓杰等[24]对西北地区建立极限学习机(extreme learning machine，ELM)模型时，基于单个气象因子的ELM模型精度最低(RMSE高于1.300 mm/d，NSE在0.500以下)；KISI等[25]对Ankara和Kirikkale地区比较了模型ELM，WANN和ANN在不同气象因子组合下的精度，结果发现基于少量气象因子建立的这些模型精度均较低(RMSE均高于1.200 mm/d，NSE低于0.600).这主要是由于模型结构不同导致，模型NARX比模型ELM和WANN在输入层增加了一个外部输入反馈，使得模型可以根据历史值调整输出，从而逼近实际值.

文中应用NARX模型仅进行了单站点短期ET0的预测，但对于多站点长序列的ET0预测以及和其他神经网络模型对比，还有待后续进一步地研究.

4 结 论

以攀枝花为研究站点，基于5个气象因子作为输入参数，PM模型计算的ET0作为输出，构建了动态神经网络中NARX和NAR模型，并且与其他模型进行对比分析，得到以下主要结论.

1) 不同气象因子组合下的NARX模型对比发现，输入单个气象因子的模型精度普遍高于其他组合，基于u2的NARX-9模拟精度最高，其RMSE为0.285 mm/d，MAE和MBE分别为0.237 mm/d和0.019 mm/d，NSE为0.964，GPI排名第1，能够很好地实现该地区缺失气象资料下的ET0预测.

2) 将模型NARX-1与NAR对比发现，带有外部输入的NARX-1模型精度高于无外部输入的NAR模型.

3) 通过模型NARX-4，NARX-7和模型Hargreaves-Samani，Irmak-Allen，Makkink以及Priestley-Taylor对比分析发现，基于Tmax，Tmin，n时，模型NARX-4的RMSE，MAE，MBE均在0.500 mm/d以下，NSE高于0.950，而模型Irmak-Allen，Makkink，Priestley-Taylor的RMSE，MAE均在0.600 mm/d以上，NSE均在0.700以下，精度相对较差；基于Tmax，Tmin时，模型NARX-7(RMSE<0.500 mm/d，NSE>0.940)精度高于模型Hargreaves-Samani(RMSE>0.600 mm/d，NSE<0.700).因此，在气象资料缺失条件下，NARX模型相较于其他简化计算模型的精度高，更适合用于模拟计算ET0.