不同气候区参考作物需水量计算方法对比研究

崔伟敏，宋 妮，申孝军，粱悦萍，黄仲冬，宁慧峰，陈新国

(1. 中国农业科学院农田灌溉研究所/农业部作物需水与调控重点开放实验室，河南 新乡 453003；2. 石家庄市水利水电勘测设计研究院，石家庄 050001)

0 引 言

农田蒸散是表征农田土壤水分通过农田土壤表面蒸发以及作物植株蒸腾散发到大气系统过程的重要指标，也是气象、水文以及农田水利、灌溉管理的重要参数之一[1]。而农田蒸散量的定量化描述通常通过参考作物需水量(ET0)来间接估算，Allen等人[2]将参考作物需水量定义为一种理想作物的需水量，假设作物表面开阔、高度一致(0.12 m)、生长旺盛、完全遮盖地面(叶面阻力为70 sm-1，反射率为0.23)，无病虫害且没有水分胁迫。

准确估算ET0对于提高灌区灌水效率及灌溉水资源利用效率具有重要的意义。目前，国内外参考作物需水量的估算方法主要有蒸发皿法[3]和基于温度、基于辐射、基于能量平衡等数学方法，但这些方法均有一定的适用性[4]。Penman-Monteith公式是被推荐确定ET0的标准方法，经过不断的改进和发展，该方法已能准确计算不同地区和气候条件下的ET0，并在世界范围内被广泛采用，并作为其他估算模型的标准[5,6]；但该方法最大的不足表现在，所需的参数较多，而且部分参数在很多区域很难获得。特别是在一些发展中国家，诸如辐射、相对湿度、风速等气象数据资料还很不完善[4,7]。此外，安装和维护气象设备也需要投入大量的人力和物力。

当不能获得完整气象资料时，可以利用所需气象参数较少的估算模型确定ET0，但由于气候条件和区域等自然条件的差异，在利用这些估算模型确定ET0前，需要利用PM标准模型对简化公式进行率定和评估。国内外学者在模型适用性方面已经进行了大量研究，也取得了丰硕的成果，但这些研究大多集中于半干旱地区相关模型适用性评价，DehghaniSanij 等人研究表明[8]，半干旱条件下，与大型蒸渗仪观测结果相比，PM公式的计算结果最可靠；Turc模型更适宜估算冷湿环境条件下的ET0，在暖湿和半干旱气候区Hargreaves模型估算结果精度最高[9]。

但有关类似我国西北极端干旱荒漠绿洲、干旱区、半干旱区以及南方湿润特湿润气候区参考作物需水量计算方法综合比较评价的研究成果还鲜有报道。为此，笔者选取我国不同气候区9个典型站点建站以来的历史气象资料，分别利用Penman-Monteith公式、Hargreaves公式、McClound公式、Thornthwaite公式、Rohwer公式、Albrecht公式、Romanenko公式、Brockamp and Wenner公式和WMO公式计算了不同站点的ET0，以期筛选出我国不同气候区适宜的ET0估算模型，为气象资料缺失地区参考作物需水量及农田蒸散量的精确估算提供理论依据及数据支持。

1 材料与方法

1.1 数据来源

本研究中使用的数据集来自中国气象科学数据共享服务网(http:∥cdc.cma.gov.cn/home.do)，所选站点分别为阿克苏站(41.17° N，80.23° E，高程1 103.8 m)、乌鲁木齐站(43.78° N，87.65° E，海拔935.0 m)、兰州站(103.88° N，36.05° E，海拔1 517.2 m)、呼和浩特站(111.68° N，40.82° E，海拔1 063.0 m)、西安站(108.93° N，34.30° E，海拔397.5 m)、北京站(116.47° N，39.80° E，海拔31.3 m)、西宁站(101.75° N，36.73° E，海拔2 295.2 m)、广州站(113.33° N，23.16° E，海拔41.0 m)和昆明站(102.62° N，25.00° E，海拔1 888.1 m)，其中位于中国西北极端干旱荒漠绿洲区阿克苏和乌鲁木齐站以及位于干旱地区的兰州站和呼和浩特站可以代表温带大陆气候区(TCZ)，位于半干旱地区的西安站和北京站可以代表温带季风气候区(TMZ)，位于海拔2300 m的西宁站可以代表高原山地气候区(MPZ)，位于中国南部的昆明站可以代表副热带季风气候区(SMZ)，位于中国东南沿海的广州站可以代表热带季风气候区(TRMZ)[10]。

1.2 参考作物需水量计算方法

1.2.1 FAO推荐的基于Penman-Monteith公式

Rn=Rns-Rnl

(2)

Rns=(1-0.23)Rs

(3)

Rs=(as+bsn/M)Ra

(4)

N=24×ωs/π

(5)

Ra=24×60/πGscD[ωssinφsinδ+cosφcosδsinωs]

(6)

D=1+0.033cos(DOY2π/365)

(7)

ωs=arccos[-tanφtanδ]

(8)

δ=0.409sin(2π/365DOY-1.39)

(9)

Rso=(as+bs)Ra

(11)

ea=RHmean/100 [e0(Tmax)+e0(Tmin)]/2

(12)

e0(T)=0.610 8 exp[17.27T/(T+237.3)]

(13)

Δ=4 098e0(Tmean)/(Tmean+237.3)2

(14)

es=e0(Tmean)

(15)

γ=CpP/(ελ)

(16)

P=101.3 [(293-0.006 5z)/293]5.26

(17)

式中：ET0为参考作物需水量，mm/d；△为饱和水汽压与温度关系曲线的斜率；Rn为净辐射量，MJ/(m2·d)；G为土壤热通量，MJ/(m2·d)；γ为湿度计常数，kPa/℃；T为空气平均温度，℃；u2为地面2m高处的风速，m/s；es为空气饱和水汽压，kPa；ea为空气实际水汽压，kPa;Rns是净短波辐射通量；Rnl是净长波辐射通量；Rs是入射到地表的太阳短波辐射通量；Ra是入射到大气外层太阳短波辐射；n是实际日照时数；N是最大可能日照时数；as，bs是随季节而变化的常数；δ是太阳赤纬；Gsc是太阳常数，其值为0.082 0MJ/(m2·min)；D是以日地平均距离为单位的日地之间的距离，无量纲；ωs是日落时角；ψ是纬度；DOY是一年内的日序数，即一年内的第几天；σ是斯特潘—玻耳兹曼常数，5.67×10-8W·m2·K4；Tmax,k是日内最大绝对温度；Tmin,k是日内最小绝对温度；Rso是晴空入射到地表的太阳短波辐射；RHmean是日平均相对湿度；e0(T)是温度为T时的饱和蒸汽压；Tmax，Tmin和Tmean分别是日最大、最小和平均气温；P是气压；Cp是空气定压比热，其值为1.013×10-3MJ/(kg·℃)；ε是水蒸汽与干燥空气的分子量比率，其值为0.622；λ是水汽化潜热，2.45MJ/kg；z是高程。

1.2.2PMT方法

适应于气象资料不全，仅有日平均最高气温和日平均最低气温的情况下，风速取多年平均风速；在计算净辐射项中，太阳短波辐射Rs计算：

式中：KRs是调节系数，内陆取值常为0.17；Tmax、Tmin分别是日平均最高气温和日平均最低气温，℃；平均相对湿度RHmean缺失时，实际水汽压ea可用下式进行计算：

1.2.3Hargreaves公式

该公式由Hargreaves和Allen等人于1985年构建，利用Hargreaves公式估算参考作物需水量只需要当地日最高和最低温度，外空辐射数据可根据当地纬度和日序数进行估算。其计算公式：

式中：ET0,HS是用Hargreaves公式计算得到的ET0值，mm/d；Ra是大气顶太阳辐射，可根据日序数和站点地理纬度计算，MJ/(m2·d)；C、E、Toff为公式中的3个参数，推荐取值分别为0.002 3，0.5和17.8[11]。

1.2.4McClound公式

ET0,M=KW1.8T

(21)

式中：ET0,M是利用McClound计算的参考作物需水量，mm/d；T是日平均气温，℃；K=1.243，W=1.03[10]。

1.2.5Rohwer公式

ET0=0.44 (1+0.27μ) (es-ea)

(22)

1.2.6WMO公式

ET0=(0.129 8+0.093 4μ) (es-ea)

(23)

1.2.7Romanenko公式

Et0=0.001 8(Ta+25)2(100-RH)

(24)

1.2.8BrockampandWenner公式

ET0=0.543μ0.456(es-ea)

(25)

1.3 评价标准

本文以Penman-Monteith公式计算结果为依据，利用均方根误差(RMSE)、误差百分比(PE)、平均偏差(MBE)和确定系数(R2)来评价不同计算方法在不同气候区的适用性。

2 结果与分析

2.1 不同区域ET0计算结果

根据中国不同典型区域代表站点多年气象资料，利用8种常用ET0估算模型计算了各站点逐日参考作物需水量，图1～图5分别给出了温带大陆气候区(TCZ)、温带季风气候区(TMZ)、高原山地气候区(MPZ)、副热带季风气候区(SMZ)和热带季风气候区(TRMZ)典型站点计算结果的多年平均值。从图1可以看到，阿克苏站、乌鲁木齐站和呼和浩特站计算结果相似，Rohwer和WMO公式计算结果明显高于其他6种公式计算结果，而且其他6种公式计算结果比较接近，其计算结果在年内呈现先增加后减小的趋势。温带季风气候区计算结果(图2)表明，北京站Rohwer和WMO公式计算结果普遍高于其他6种公式计算结果，特别是5月份Rohwer和WMO公式计算结果平均值分别比PM公式计算值高出13.45和9.22mm/d，西安站的差值相对较小。高原山地气候区、副热带季风气候区和热带季风气候区不同公式计算结果差异相对较小，而且年内变化相对较为平缓。

2.2 不同计算方法对ET0计算结果的影响

以PM公式逐日ET0计算结果为标准对不同气候区利用不同方法计算结果进行评价，统计分析结果(表1)表明，与PM公式计算结果相比，不同计算方法的计算结果随着气候区域的不同而存在不同程度的差异。就温带大陆气候区(TCZ)而言，Hargreaves公式的表现最好(RMSE、PE、MBE和R2值分别为0.47mm/d、12.93%、0.16mm/d和0.99)，其次是McClound公式(RMSE、PE、MBE和R2值分别为0.79mm/d、10.28%、-0.28mm/d和0.91)和BrockampandWenner公式(RMSE、PE、MBE和R2值分别为0.84mm/d、17.93%、0.26mm/d和0.97)，其他4种公式的计算结果均不理想，PE值均大于30%。就温带季风气候区(TMZ)而言，BrockampandWenner公式(RMSE、PE、MBE和R2值分别为0.52mm/d、7.62%、0.06mm/d和0.91)的表现最好，其次是Hargreaves公式(RMSE、PE、MBE和R2值分别为0.54mm/d、12.18%、0.30mm/d和0.98)和McClound公式(RMSE、PE、MBE和R2值分别为0.71mm/d、10.88%、0.23mm/d和0.83)，其他4种公式的计算结果均不理想，PE值均超过了22.28%。就高原山地气候区(MPZ)而言，Hargreaves公式(RMSE、PE、MBE和R2值分别为0.32mm/d、7.41%、0.19mm/d和0.99)的表现最好，其次是Romanenko公式(RMSE、PE、MBE和R2值分别为0.40mm/d、7.47%、0.19mm/d和0.96)，其他5种公式的计算结果均不理想，PE值均超过了21.08%。副热带季风气候区(SMZ)计算结果分析表明，Hargreaves公式(RMSE、PE、MBE和R2值分别为0.46mm/d、13.66%、0.40mm/d和0.94)的表现较为理想，McClound公式和Romanenko公式虽然PE值均较小，但R2值过低。热带季风气候区(TRMZ)考虑所有的分析结果表明，Hargreaves公式表现相对较好(RMSE、PE、MBE和R2值分别为0.77mm/d、24.03%、0.67mm/d和0.78)。

注：图中数值为ET0月平均值(即首先计算出逐日ET0，然后根据计算结果求得逐日ET0多年平均值，最后计算出ET0月平均值)，下同；阿克苏(1953年6月1日-2014年12月31日)，乌鲁木齐(1951年1月1日-2014年12月31日)，兰州(1951年1月1日-2014年12月31日)，呼和浩特(1951年2月1日-2014年12月31日)

注：北京(1951年1月1日-2014年12月31日)，西安(1951年1月1日-2014年12月31日)

图3 高原山地气候区(西宁)不同计算方法ET0计算结果(1954年1月1日-2014年12月31日)Fig.3 The result of ET0 calculated by different models in MPZ(1954.01.01-2014.12.31)

2.3 不同气候区ET0适宜计算方法及相关参数率定

综合考虑不同公式计算结果精度、所需参数等因素，可以初步确定Hargreaves、McClound和BrockampandWenner公式比较适合温带大陆气候区(TCZ)ET0计算；Brockamp and Wenner、Hargreaves和McClound比较合适温带季风气候区(TMZ)ET0计算；Hargreaves和Romanenko公式比较适合高原山地气候区(MPZ)ET0计算；Hargreaves公式比较副热带季风气候区(SMZ)和热带季风气候区(TRMZ)ET0计算。

图4 副热带季风气候区(昆明)不同计算方法ET0计算结果(1951年1月1日-2014年12月31日)Fig.4 The result of ET0 calculated by different models in SMZ(1951.01.01-2014.12.31)

图5 热带季风气候区(广州)不同计算方法ET0计算结果(1951年7月1日-2014年12月31日)Fig.5 The result of ET0 calculated by different models in TRMZ(1951.07.01-2014.12.31)

表1 不同方法计算结果与PM计算结果统计分析Tab.1 Statistical performance of the different methods versus the PM 56 model for estimating ET0

为了进一步提高不同气候区ET0计算精度，利用不同气候区典型气象站2013年12月31号以前的历史气象资料以及得用PM公式计算的ET0结果对Hargreaves和McClound公式中的参数进行率定。提出了不同气候区适宜的参数(表2)。从表2可以看到，不同地区Hargreaves和McClound公式的参数存在较大的差异，即使是同一地区不同站点也存在不同程度的差异，C值介于0.015 1～0.004 21，E值介于0.253 8～0.668 3；Toff值介于9.83～36.26。McClound公式的参数在不同气候区间的变异相对较小，K值介于0.825 0～1.483 1；W值介于1.025 3～1.036 9之间。

表2 不同气候区ET0计算方法及相关参数Tab.2 The parameters of models for the ET0 in different climatic region

根据2014年实测气象资料，利用PM公式计算结果对不同气候区适宜ET0计算公式计算结果进行验证，试验结果(表3)表明，与PM公式计算结果相比，Hargreaves公式计算的RMSE、PE和R2值分别介于0.42～0.95 mm/d、0.81%～10.69%和0.66～0.94，McClound公式计算的RMSE、PE和R2值分别介于0.80～1.23 mm/d、1.36%～14.52%和0.60～0.87，在温带大陆气候区(TCZ)Brockamp and Wenner公式计算结果的PE值普遍偏高，而在TMZ和MPZ地区的计算结果较为理想。

3 结 语

(1)与PM公式计算结果相比，8种常用ET0计算结果在不同区域的精度存在不同程度差异，即使同一公式计算结果在不同气候区的计算精度也存在较大差异。

表3 不同气候区适宜方法计算结果与PM计算结果统计分析Tab.3 Statistical performance of the different methods versus the PM 56 model for estimating ET0

(2)综合考虑计算精度、参数获得难易程度以及可操作性，初步筛选出了适宜不同气候区的简易ET0计算模型。

(3)利用历史气象资料率定了不同气候区典型站点的Hargreaves和McClound公式的参数，并利用2014年实测气象资料对其进行验证，结果表明，通过对参数优化可以进一步优化筛选的ET0计算公式计算精度。