基于气象-生理的夏玉米作物系数及蒸散估算

王振龙，范 月，吕海深，许莹莹，董国强

（1. 安徽省（水利部淮委）水利科学研究院 水利水资源安徽省重点实验室，蚌埠 233000；2. 河海大学理学院，南京 210098；3. 河海大学水文水资源学院，南京 210098）

0 引 言

蒸散是植物蒸腾量与植株间土壤蒸发量之和，是土壤-植物-大气系统（Soil-Plant-Atmosphere Cotinuum，SPAC）中重要的水交换过程[1]。其对作物的生长发育情况有重要影响，在大气循环中扮演重要角色，还能进一步调节气候[2-3]。构建精准的作物系数模型对准确预测蒸散量及进一步制定精准的灌溉计划有重要意义。

目前针对作物蒸散量计算方法，主要方法是作物系数法。世界粮农组织（Food and Agriculture Organization of the United Nations，FAO）提出了一种估算作物蒸散量的标准方法（参考作物蒸散量乘以作物系数），其中参考蒸散通过修正的彭曼（Penman-Menteith，P-M）模型计算得到，目前该方法已得到广泛运用[4]，而其中的关键问题就是确定作物系数。然而，根据FAO 推荐的4 个生长阶段的作物系数计算得到的蒸散量与实际蒸散可能存在较大偏差[4-5]，因此国内外学者建立了多种作物系数模型。如Steele 等[6]利用回归分析建立了作物系数与播后天数间的数学关系式，生成了玉米生长期内作物系数的动态变化曲线；王娟等[7]利用气象数据和P-M 模型构建了夏玉米蒸散模型，并得出了不同生育期的作物系数；雷志栋等[8]采用气象系列数据，通过P-M 模型和非水分胁迫条件下作物腾发（Non-Water Stress Evapotranspiration，NWSE）模型模拟计算了冬小麦返青-收获期逐日作物系数；刘海军和康跃虎[9]在冬小麦拔节抽穗期探讨了农田实际蒸散量和参考作物蒸散量的关系，利用回归分析建立了叶面积指数和株高与作物系数的数学关系；王振龙等[10]基于温度效应对作物生长发育的影响从而计算作物系数。尽管这些模型取得了一定成果，但综合考虑作物系数影响要素的具体计算模型较少，且大多都是通过传统方法计算拟合。本研究尝试采用地温及叶面积指数构建了气象-生理双函数乘法模型，并融合了机器学习算法（涉及多领域交叉学科），进一步提高了作物系数计算精度。

本研究采用五道沟水文实验站大型称重式蒸渗仪实测蒸散数据，结合系列气象资料计算出实际作物系数，将作物系数与叶面积指数利用米氏方程进行拟合，得出受生物因子影响的作物系数，进而利用多元回归筛选出影响作物系数最密切的气象因子，即地温（0 cm），进行了指数拟合，最后通过梯度下降法优化求取模型参数，采用实测数据对计算结果进行评估，对掌握夏玉米生育期内作物系数的动态变化特征，精确估算作物实际蒸散量具有重要意义。

1 材料与方法

1.1 试验区概况

试验研究在安徽淮北平原五道沟水文实验站进行。该地区（33°09’N，117°21’E）属暖温带半湿润季风气候区，四季分明，雨热同期。地下水位变化范围为1～3 m，种植作物主要有玉米、小麦、花生和大豆。试验区土壤主要包括2 种：砂姜黑土和黄潮土，分别约占淮北平原土壤的54%和33%。试验区土壤田间持水率为28%～30%，适宜作物生长的土壤含水率18%～25%，凋萎含水率为10%～13%。试验区设有气象观测场，可获取地温、空气温度及湿度、降雨量、风速及风向、水面蒸发量和日照时数等60 余年不间断长系列气象要素。其中地温有0、5、10、15、20、40、80、160 和320 cm 长系列资料。另外还设有自动高精度气象站，可获取空气温度、空气湿度、风速、净辐射、太阳辐射、土壤热通量等，可用于计算参考作物蒸散量。

1.2 试验设计与方法

本研究以淮北平原分布较为广泛的砂姜黑土区为研究对象，作物为夏玉米（登海618），于2018 年6 月22日播种，10 月8 日收获。玉米实际蒸散量由大型称重式蒸渗仪测得，蒸渗仪土柱口径为4.0 m2，高度为4.0 m，蒸渗仪地下水埋深设为1 m，土壤为分层回填土，10、30、50 和100 cm 埋深处分别设有土壤水分、温度、电导3 个参数传感器，数据每10 min 获取1 次。蒸渗仪自动采集数据资料时段选取2018 年6 月22 日至10 月8 日。选取同期气象观测场水文气象要素数据。利用叶面积仪测定蒸渗仪内玉米不同生育阶段的叶面积，结合实际种植密度计算叶面积指数。

1.3 蒸渗仪数据处理

由于蒸渗仪会受天气、仪器自身因素等的影响，需对数据进行预处理：1）对蒸渗仪蒸散量异常值进行剔除和插值处理；2）对1 h 内实际蒸发筒重和10 min 差值分别平均再求和计算整点蒸发筒重，计算出时变化量，再通过移动加权平均法修正1 h 变化量，得出时蒸散量；3）对小时蒸散量累计计算日蒸散量。

1.4 参考作物蒸散量

本研究采用高精度气象站数据，利用FAO 修正的P-M 公式计算逐日参考作物蒸散量（Reference Evapotranspiration，ET0，mm/d），如式（1）[11-14]所示

式中Rn为地表净辐射，MJ/(m2·d)；G 为土壤热通量，MJ/(m2·d)；U2为地面上方2 m 处风速，m/s；es为饱和水汽压，kPa；ea为实际水汽压，kPa；Δ 为饱和水汽压随温度变化曲线的斜率，kPa/℃；γ 为干湿球常数，kPa/℃；T为2 m 高处日平均气温，℃，由日最高气温和最低气温平均值计算得出。

1.5 作物系数计算

本研究所用测筒控制地下水埋深为1 m，属于不充分供水，这种情况下作物实际蒸散量（Evapotranspiration，ET，mm/d）除了受环境因素和作物自身特性的影响外，还受到土壤含水率的制约，因此，本研究引入了土壤水分胁迫系数Ks，并进行实际作物系数的计算[15]，如式（2）所示

式中Kc为作物系数。

1.6 土壤水分胁迫系数计算

土壤水分胁迫系数反映了土壤水分对蒸散量的影响，通过式（3）[16]确定为

式中θ(0-40)为0～40 cm 的平均土壤质量含水率，由蒸渗仪内10 cm 和30 cm 体积含水率平均计算所得，土壤容重为1.4 g/cm3；θw为凋萎含水率，θw= 10%；θf为田间持水率，θf= 28%。

2 作物系数模型建立

2.1 模型构建

2.1.1 基本模型

作物系数可以反映作物生长状况。当前作物系数的计算一般基于作物生理特性，但以往研究表明：作物系数受到作物生理特性及气象因子的共同影响。因此，本研究综合考虑两种类型因素对作物系数（Crop Coefficient，Kc）的影响[17]，通过建立气象-生理双函数乘法模型计算作物系数，模型式如式（4）所示

式中Ke为由气象因素构成的作物系数函数，KLAI为由作物生理因素构成的作物系数函数，其中LAI（Leaf Area Index)为叶面积指数，m2/m2。

当前，不同研究者针对作物系数与叶面积指数间关系采取了不同拟合方式，如高斯拟合、指数拟合等。在生长初期，作物快速生长，蒸散随作物蒸腾量的增加而增加，作物系数也增大，当增大到一定程度时，LAI 与作物系数的关系不再明显，作物系数主要受环境因子影响，本研究根据作物系数对叶面积指数的响应特征，遂采用米氏（Michaelis-Mentenn）方程拟合作物系数与LAI[18-19]，即KLAI为式（5）所示

式中m 和n 为待拟合参数；KLAI为仅受生物因子影响时的作物系数。

将实际Kc去除生物因素影响的处理，即用实际Kc除以拟合得到的KLAI，得Ke，分析其与同期各气象因子包括日平均气温、日最高气温、日最低气温、日照时数、风速、绝对湿度、相对湿度、饱和差、水汽压力差和0、5、10、20、40、80、160、320 cm 地温的相关性，采用多元逐步回归法筛选出解释程度较高的因子，并对各要素进行共线性诊断，根据SPSS 分析结果可知各个系数的方差膨胀因子（Variance Inflation Factor，VIF）较小（0

式中a、b、c 为未知参数；D0为地表温度，℃。

表1 各函数拟合优度和均方根误差比较

2.1.2 梯度下降法

针对上述模型中的未知参数，本研究利用梯度下降法确定参数，梯度下降算法是机器学习中较为广泛的优化算法[20-22]，本研究主要以线性回归算法损失函数求极小值确定参数a、b、c，因此需将式（6）转化为线性形式，其步骤如下：

1）首先对式（6）进行转化，对其两边取对数，得式（7）为

式中θ0、θ1和θ2为未知参数；x0为特征常量设为1；x1为地温，℃；x2为地温的平方，℃2。

为了算法能更快收敛，要确保不同特征的取值在相近的范围内，需要对自变量进行特征归一化。式（8）中D0和D02取值范围相差较大，需要对其进行特征缩放，即为式（10）和式（11）所示

式中x1'为均值归一化后的地温，℃；x1av为地温平均值，℃；x1max为地温最大值，℃；x1min为地温最小值，℃；x2'为均值归一化后的地温平方，℃2；x2av为地温平方平均值，℃2；x2max为地温平方最大值，℃2；x2min为地温平方最小值，℃2。

因此将式（9）转化为式（12）所示

2）对式（12）利用梯度下降法求得参数；在线性回归算法中，目标损失函数定义为式（13）所示

式中J(θ)为损失函数；θ 为由θ0、θ1和θ2构成的参数向量；m 为样本总数；xi为由x0、x1'及x2'构成的特征向量；yi为实际作物系数构成的向量；hθ(xi)为对应预测值构成的函数。具体过程如下：对θ 赋初值，这个值可以随机选取，也可以让其为一个全0 的向量，本研究令其初始值为0 向量；更新θ 的值，使得目标损失函数按梯度下降的方向减少，通过式（14）迭代

式中α 为学习率，需人为指定，本研究取值0.001，α 若过大会导致震荡无法收敛，若过小收敛速度会很慢；对多元函数的参数求偏导，将求得的各个参数的偏导数以向量的形式表示就是对应的梯度，本式中为对应的梯度，梯度计算结果如式（15）所示

3）当每次更新后的结构都能让损失函数变小，最终达到最小则停止下降。

2.2 模型评价指标

本研究评价指标包括平均绝对误差（Mean Absolute Deviation，MAE）、均方误差（Mean-Square Error，MSE）、均方根误差（Root Mean Square Error，RMSE）和相关系数（Correlation Coefficient，r），用来评价计算结果与实测数据间的误差情况；采用准确率（Precision，P）对模型精度进行评价，对于作物系数，定义绝对误差在0.2、0.3 和0.4 以内的数据个数占总数据集的比例为相应准确率，分别记为P0.2、P0.3和P0.4；对于蒸散量，定义绝对误差在2、3 和4 mm 以内的数据个数占总数据集的比例为相应准确率，分别记为P2、P3和P4。各指标通过式（16）～式（19）计算

式中yi为实际值；i为模型计算值；i 为样本序数，i=1,2,…m；为实际值均值；为计算值均值；m 为样本个数。

3 结果与讨论

3.1 作物系数模型

夏玉米全生育期地温、叶面积指数及实际作物系数变化（图1）。地温、LAI 与Kc的r 分别为0.75 和0.52，相关性较高。作物根系生长对地温变化非常敏感，地温越高，作物的生长发育越快，作物蒸腾作用增强，土壤水的移动越频繁；地温低，根系的生长和代谢都会受到抑制，作物根系吸水缓慢，蒸腾作用较弱。且作物蒸腾主要通过叶片蒸腾作用，与叶面积关系密切[23-26]。

图1 全生育期地温和作物系数Kc 变化过程线图 Fig.1 Process of changes in ground temperature and crop coefficient Kc during the whole growth period

利用蒸渗仪实测蒸散数据及气象系列资料对作物系数模型中各参数进行估计。首先对于KLAI，选取2018 年玉米生长季（6 月22 日至10 月8 日）的实测叶面积指数数据进行研究，每隔1 d 抽取1 d 数据作为训练样本（55 d），其余作为检验样本（54 d）。将根据式（2）计算得到的作物系数Kc与实测叶面积指数LAI 进行拟合，利用Matlab拟合为式（5）的形式，所得方程如式（20）所示

式中0.158 4 为米氏常数，数值越小，作物系数对LAI 变化敏感性越强；由于时间尺度为日尺度，故LAI 均对应逐日作物系数Kc，故R2较小，但两者之间具有显著相关性（P<0.01）。

其次针对Ke，利用梯度下降法（如式13～式15）解得式（12）中的参数θ0、θ1和θ2分别为-0.063 4、0.542 8和0.478 8，再根据所得θ0、θ1和θ2的值求得式（6）中a、b 及c 的值分别为0.804 7、-0.566 8 和1.445 2。所得Ke计算式如式（21）所示

最后，根据式（20）及式（21）建立的乘法模型式如式（22）所示

图2 为训练样本根据式（22）所得作物系数计算值及实际值变化过程图。两者平均值分别为0.399 和0.405，相差小于0.01，MAE 为0.11，MSE 为0.02，RMSE 为0.14，对于作物系数估算准确率P0.2为80%、P0.3为98%、P0.4为100%，误差较小，说明模型精度较高。

图2 训练样本计算作物系数Kc 与实际值比较 Fig. 2 Calculated crop coefficient Kc compared with actual value in training sample

3.2 模型检验

图3 为检验样本根据式（22）所建模型计算得到的作物系数Kc及实际值变化过程图。两者平均值分别为0.40 和0.44，仅相差0.04，MAE 为0.12，MSE 为0.02，RMSE 为0.15，对于作物系数估算准确率P0.2为81%、P0.3为96%、P0.4为100%，误差较小，模型预测精度较高。

图3 检验样本计算作物系数Kc 与实际值比较 Fig.3 Calculated crop coefficient Kc compared with actual value in testing sample

全生育期夏玉米实际作物系数与计算作物系数均值分别为0.42 和0.40，仅相差0.02，MAE 为0.12，RMSE为0.15，r 为0.91。对于作物系数估算准确率P0.2为81%、P0.3为97%、P0.4为100%。

表2 为以实际作物系数为依据，不同样本范围模型的计算精度评价指标值。结果表明本研究构建的作物系数模型能以较高精度计算夏玉米生长期内逐日作物系数，再利用联合国粮农组织（FAO）所推荐的估计作物实际蒸散量方法，可以进一步计算作物逐日蒸散量ET。

表2 不同样本作物系数计算精度评价指标值 Table 2 Evaluation index value of crop coefficient calculated accuracy in different samples

3.3 蒸散量计算及精度评价

图4 为依据式（2）计算所得训练集蒸散值与实际值对比图。两者平均值分别为3.88 和3.86 mm/d，相差小于0.1 mm/d，MAE 为1.0 mm/d，RMSE 为1.4 mm/d，对于蒸散量估算准确率P2为80%、P3为95%、P4为98%，误差较小，预测精度较高。

图5 为依据式（2）计算所得验证集蒸散值与实际值对比图。两者平均值分别为4.01 和4.16 mm/d，仅相差0.15 mm/d，MAE 为1.0 mm/d，RMSE 为1.3 mm/d，对于蒸散量估算准确率P2为81%、P3为96%、P4为100%，误差较小，模型预测精度较高。

图5 检验样本计算蒸散量ET 与实际值比较 Fig.5 Calculated evapotranspiration ET compared with actual value in testing sample

表3 为以实际作物系数Kc为依据，不同样本范围模型的计算精度评价指标值。结果表明本研究构建的作物系数模型能以较高精度计算夏玉米生长期内逐日作物系数，并利用 FAO 所推荐的估计作物实际蒸散量（Evapotranspiration，ET）方法所计算的作物逐日蒸散量ET 精度也较高。

为防止出现过拟合问题，本研究选取夏玉米6 月22日至7 月22 日连续时间作模型验证，所得蒸散量ET 实际值与估算值如图6 所示。两者平均值分别为2.33 和2.54 mm/d，仅相差0.21 mm/d，MAE 为0.83 mm/d，RMSE为2.82 mm/d，r 为0.85，精度较高。

表3 不同样本蒸散量计算精度评价指标值 Table 3 Evaluation index value of evapotranspiration calculated accuracy in different samples

图6 连续时段计算蒸散量ET 与实际值比较 Fig.6 Calculated evapotranspiration ET compared with actual value in continuous period

3.4 模型比较与评价

本研究采用相同数据并根据FAO 推荐作物系数及P-M 模型计算了全生长季的玉米蒸散量。FAO 建议将作物的生长期划分为生长初期、发育期、生长中期和生长后期这4 个阶段[27-30]。根据试验地区作物实际生长情况，夏玉米各生长阶段划分及对应的推荐作物系数如表4，其中发育期作物系数通过插值得到。

表4 夏玉米生长阶段划分及对应作物系数 Table 4 Summer maize growth stage division and corresponding crop coefficient

将计算结果与本研究作物系数模型计算的玉米蒸散量分别与其实际蒸散量进行了对比（图7），发现三者的趋势基本一致，但根据FAO 推荐作物系数计算的蒸散量ET 值整体偏高，计算精度明显低于本研究所建作物系数模型，具体计算结果如表5 所示。

基于本研究模型计算所得蒸散值与实际值均值仅相差0.07 mm/d，MAE 为1.00 mm/d，RMSE 为4.50 mm/d，r 为0.75，对于蒸散量估算准确率P2为81%、P3为95%、P4为99%。根据FAO 推荐作物系数计算所得蒸散值与实际值存在较大偏差，两者均值相差5.00 mm/d，MAE 为5.07 mm/d，RMSE 为6.09 mm/d，r 为0.75，对于蒸散量估算准确率P2为23%、P3为30%、P4为36%。本研究模型计算的蒸散值准确率相比FAO 推荐作物系数得到的计算值准确率（误差在2～3 mm/d 以内）提高了3 倍以上，且误差几乎都在4 mm/d 以内。因此本研究构建的气象-生理双函数模型，综合反应了气象和生物因子对玉米生长的共同影响，可以更高精度估算玉米作物系数和蒸散量。

表5 全生育期不同模型蒸散量计算精度评价指标值 Table 5 Evaluation index value of evapotranspiration calculated accuracy by different models during the whole growth period

图7 全生育期不同方法蒸散量ET 计算值与实际值比较 Fig.7 Calculated evapotranspiration ET compared with actual value based on different ways during the whole growth period

4 结 论

1）本研究基于地温和叶面积指数构建了气象-生理双函数乘法模型，可以反映作物逐日作物系数变化。该模型在训练和检验样本中计算精度均较高。在训练样本中，作物系数计算值与实际值均值相差0.01，平均绝对误差（Mean Absolute Deviation，MAE）为0.11，均方根误差（Root Mean Square Error，RMSE）为0.14，相关系数（Correlation Coefficient，r）为0.94；在检验样本中，两者均值相差0.04，MAE 为0.12，RMSE 为0.15，r 为0.87。该模型简化了夏玉米作物系数的计算，明确了地温和叶面积指数对作物系数的综合影响程度，提高了计算精度，可用于玉米作物系数的动态计算。

2）基于本研究作物系数模型计算的实际蒸散量精度较高，在训练样本中，蒸散量计算值与实际值均值相差0.1 mm/d，MAE 为1.0 mm/d，RMSE 为1.4 mm/d，r 为0.75；在检验样本中，两者均值相差0.15 mm/d，MAE 为1.0 mm/d，RMSE 为1.3 mm/d，r 为0.76。与联合国粮农组织（FAO）推荐作物系数计算结果相比，精度更高。

3）根据FAO 推荐作物系数计算蒸散值与实际值两者均值相差5.0 mm/d，MAE 为5.07 mm/d，RMSE 为6.09 mm/d，r 为0.75，本研究模型计算的蒸散值准确率相比FAO推荐作物系数得到的计算值准确率（误差在2～3 mm/d 以内）提高了3 倍以上。

本研究作物系数模型是基于五道沟水文实验站大型称重式蒸渗仪实测数据构建，本模型研究时间及空间尺度还有一定的局限性，通过蒸散发、潜水蒸发和包气带土壤水运移等多要素、多过程、面尺度、时尺度原型监测试验有待进一步完善，以大型蒸渗仪与大田地下水位同步变化的蒸散量试验有待进一步研究和实施完善。