不同冠层阻力模型在夏玉米蒸散发计算中的优化应用

林馨贝，周岗，郑泽涛，赵璐,2*，梁川

（1.四川大学，成都610065；2.南方丘区节水农业研究四川省重点实验室，成都 610066）

0 引言

【研究意义】农田作物蒸散发（ET）包括土壤蒸发和植株蒸腾2 个部分[1]，是农田能量平衡和水分平衡的重要组成部分[2]。精确估算作物ET可为制定农田灌溉方案、评估水资源规划管理提供依据。【研究进展】如今ET测量方法较多，如可通过涡度相关方法、大型蒸渗仪法、水文学法等进行直接测定，也可通过遥感法、波文比能量平衡法、空气动力学法等间接测定[3-5]。随着测量精度不断提高，多种ET估算模型被相继提出，其中应用较多的是单源 Penman-Monteith 模型（ P-M ）、双源Shuttleworth-Wallance 模型（S-W）[6-9]、多源Clumping模型（C）[10-12]。

从理论上分析，多源复杂模型能更明确地反映ET过程和ET水分的来源，但多源复杂模型参数较多，容易出现“过度拟合”的情况导致参数拟合精度下降。同时，完整获得多源复杂模型所需的全部数据具有一定难度，数据中测量误差的累积也会影响模型精度，因此多源复杂模型的估算精度不一定高于单源模型[13-14]。在单源模型中，P-M 模型由于反映了各气候要素的综合影响，具有明确的物理意义，是ET研究中应用最为普遍的简单模型[15]。

P-M 模型估算农田ET的精度主要取决于冠层阻力（rc）的准确性。rc是P-M 模型中的重要参数，其受到叶片位置、土壤湿润状况及冠层内空气动力学特性等因素的综合影响，无法通过仪器直接测定[16]。目前rc计算的有效方法有3 种，第一种是利用实测ET值按P-M 或S-W 模型反推rc[17-18]；第二种是利用单叶气孔阻力，结合作物群体叶面积指数的空间垂直分布，计算得到群体的rc[19-21]；第三种是利用模型计算rc[22-23]。3 种方法中，利用模型法计算rc的应用较多，其中Jarvis 模型（JA 模型）全面考虑了太阳辐射、水汽压差、气温、土壤含水率及叶面积指数等因子与冠层阻力的综合关系，包含了各环境响应函数[24]，应用较为广泛。Srivastava 等[25]采用4 种方法估算rc，然后根据P-M 模型计算半湿润地区灌溉玉米的ET，得出JA 模型具有较高的可靠性；李彩霞等[26]基于JA 模型，引入冠层上、下部温差变量（ΔTc），提出了交替隔沟非充分供水下玉米rc的计算方法；周惠萍等[27]考虑氮元素对冠层导度的影响，建立了基于氮的JA 修正模型，提高了P-M 模型在不同施氮水平下估算番茄ET的准确性；李思恩等[28-29]基于JA 模型和流体运动的阻力定律，增加了土壤阻力，构建了耦合的表层阻力模型（CO 模型），CO 模型所求rc值可被看作水汽通过冠层和土壤的所克服的阻力平均值，该模型在作物生长季前期（LAI＜2）时模拟效果较好[30]。

【切入点】JA 模型和CO 模型中都采用了一些经验参数，由于不同区域环境因素和作物品种存在一定差别，因此在利用模型计算rc时，需先进行参数优化，以提高模拟精度。目前关于玉米ET研究中，大部分研究只选择多种rc的计算模型进行比较，或采用单一参数优化方法对rc的计算模型进行参数优化[25-26]，采用不同参数优化方法对同时对多种模型进行参数优化的研究较少，因此无法通过结果比较来筛选地区适用性最好的模型和参数优化方法的组合。

【拟解决的关键问题】怀来位于半干旱气候区，光热充足，适合夏玉米的生长[31]。夏玉米生长过程中耗水较多[32]，在干旱条件下易因缺水而减产[33-34]，因此需精确计算夏玉米ET以通过灌溉来保证玉米的正常生长发育。为此，本文基于怀来站点2013年的涡度相关数据和气象数据，以P-M 模型反推的实际rc为标准，分别采用最小二乘法和蚁群算法率定JA 模型和CO 模型中的相关参数，进一步估算2014年玉米的rc和ET，筛选出最优模型和参数优化方法的组合，为提高rc的计算模型在怀来地区的适用性和计算精度提供理论依据。

1 材料与方法

1.1 研究区概况

中国科学院怀来遥感综合试验站位于河北省怀来县东花园镇（115.7880°E，40.3491°N，海拔480 m）。站点所在区域处于半干旱地区，年均气温为10.1℃，最高气温39℃，最低气温达-20℃以下，年均降水总量约为370 mm，年均风速为3.4 m/s，试验站所在区域以沙质冲积土为主，主要种植夏玉米，品种为郑单958，每年5月上旬播种，9月中旬成熟[35-36]。

怀来站有涡度相关系统和自动气象站各2 套和1套土壤水分无线传感器网络，可对地表通量和气象要素进行实时观测[37]。

1.2 数据处理

怀来站点2013年和2014年的气象数据和涡度相关数据来自国家青藏高原科学数据中心提供的海河流域多尺度地表通量与气象要素观测数据集（http://data.tpdc.ac.cn/zh-hans/），内容包括净辐射Rn、土壤热通量G、显热通量H、潜热通量LE、风速W、降水量P、土壤含水率θ、空气湿度RH和气温T。仪器由于故障维修等因素，部分时间段出现数据缺失，日尺度数据由当日有效数据取均值进行插补[38]，连续多日缺失的热通量数据可以由非线性回归方法进行插补[39]。

在裸露地面或者植被较低矮处，由热力学第一定律可推导出地表能量平衡方程[40]：

式中：Rn为净辐射（W/m2）；G为土壤热通量（W/m2）；LE为潜热通量（W/m2）；H为显热通量（W/m2）。式（1）左右二端相等时，代表能量闭合，反之则不闭合。能量不闭合的现象在测量系统中普遍存在，本文采用波文比能量强制闭合法使能量闭合[41]。

式中：D为能量不闭合差（W/m2）；β为波文比；ΔLE和ΔH为修正值，分别将其添加到对应项即可得到LE和H修正后的数值。分析能量闭合度能够有效提高测量精度[42-43]。

本文在进行参数拟合时，在观测数据的基础上剔除了降水量大于0 的观测日，以提高后续参数优化精度。

1.3 计算方法

1.3.1 Penman-Monteith 模型

1965年，Monteith 从能量平衡理论和大气扩散理论出发，在 Penman 等的研究基础上得到了Penman-Monteith（P-M）模型[7,38]，计算式为：

式中：λ为汽化潜热（MJ/kg）；ET为实际蒸散量（mm）；Δ代表饱和水汽压与温度关系曲线的斜率；Rn为净辐射（W/m2）；G为土壤热通量（W/m2）；ρ为空气密度（kg/m3）；Cp为空气定压比热（J/（kg·℃））；VPD为饱和水汽压差（kPa）；ra为空气动力学阻力（s/m）；rc为冠层阻力（s/m）；γ为干湿表常数（kPa/℃）。

空气动力学阻力ra可以通过计算[40]为：

式中：zm、zh分别为测风速和测湿度的高度（zm=2.5m，zh=5m）；d为0 平面位移高度（m）；hc为作物平均高度（m）；zom为控制动量传递的粗糙度长度（m）；zoh为控制热通量和水汽传输的粗糙长度。K为卡曼常数，取0.41，uz为参考高度的风速（m/s）。

1.3.2 Jarvis 冠层阻力模型

Jarvis 模型（JA 模型）中包含了环境因素的响应函数，考虑了辐射、温度、湿度和土壤含水率对rc的影响。模型中各函数相互独立，该模型是一个多变量的综合性模型[24,44]：

式中：θF为田间持水率（m3/m3）；θW为凋萎含水率（m3/m3）。a1、a2、a3是经验参数，3 个经验参数和rcmin通过实测数据拟合得出。各参数拟合初值取a1=100[45]，a2=-0.12[24]，a3=0.0016[45]，rcmin=30 s/m[45]。

1.3.3 耦合冠层阻力模型

耦合冠层阻力模型（CO 模型）同时考虑了植物蒸腾和土壤蒸发的过程，所求rc指水汽通过冠层和土壤时克服的平均阻力，计算式为[28-31]：

1.3.4 利用Penman-Monteith 模型反算冠层阻力

利用站点观测的ET，通过P-M 模型反算冠层阻力：

式中：λETEC为修正后的潜热通量（W/m2），其余参数含义同式（1）。

1.4 敏感度分析

采用BP 神经网络模型分析rc与各影响因子的敏感度。BP 神经网络模型是一种按照误差反向传播的多层神经网络模型。本文中，模型输入值为各影响因子（Rn、T、LAI、VPD和θ），输出值为rc。使用BP神经网络模型控制输入值，排除某一影响因子，进而改变模型的输出值，通过输出值的变化比较rc对各影响因子的敏感程度。排除某一因子后输出值变化越大，则rc对被排除的因子敏感程度越高，反之越低[46-47]。

1.5 参数优化算法

1.5.1 最小二乘法

最小二乘法（LSM）是一种进行曲线拟合的方法，在处理大量数据和分析误差中具有广泛的应用。其基本原理[48]为：求一拟合目标函数φ(x)，使其在各点的偏差平方和达到最小：

式中：δ为拟合函数与原始数据差值；φ(xi)为拟合目标函数；yi为原始数据值。

1.5.2 蚁群算法

蚁群算法（ACO）通过模拟蚂蚁觅食过程中个体间依靠信息素进行交流的过程，进而得到从蚁群到食物源的最短路径的过程来进行优化[49-50]。采用蚁群算法求优化解的流程为：①预估模型参数范围，初始化蚁群A(k)。预设蚁群规模、信息素蒸发系数、信息素强度系数和蚂蚁爬行速度参数。②为每一只蚂蚁选择下一个节点。遍历所有节点后计算最佳路径。③更新信息素矩阵，目标至越优的方格信息素越强。④检查终止条件，看是否达到最大迭代次数或者最优蚂蚁是否达到要求精度。满足条件则进行下一步，不满足则返回步骤②，重复执行步骤②、③、④。⑤输出最优值。

1.6 模型精度评价参数

本文采用决定系数R2，均方根误差RMSE和一致性指数d来评价模型拟合的精度。

式中：Emi和EOi分别为模拟值和实测值；是实测值的平均值；n为样本容量。本文中，利用λETEC带入P-M 模型中反算所得的为实测值，利用模型拟合出的rc为模拟值；用相关仪器测得ETEC的为实测值，利用模型得到的ET为模拟值。

2 结果与分析

2.1 气象数据的逐日变化

2013、2014年夏玉米全生育期的净辐射Rn，风速W，温度T，饱和水汽压差VPD，降水量P和土壤含水率θ的逐日变化如图1所示。2013年和2014年玉米全生育期的Rn基本在50～250 W/m2间变化，2013年Rn的最大值为245.06 W/m2，最小值为23.04 W/m2，2014年Rn的最大值为218.95 W/m2，最小值为25.04 W/m2，6月末至7月初的Rn值在一个较高水平，之后随时间呈下降趋势；风速W的值在2014年波动较大，2013年的最大值出现在5月19日，为5.28 m/s，2014年的最大值出现在5月25日，为5.97 m/s，随后W呈下降趋势；在夏玉米的生育期内，温度T的变化范围主要在15～30℃之间，趋势大致呈抛物线形式，在7、8月达到峰值，2013年的峰值为27.36℃，2014年的峰值为28.96℃；饱和水汽压差VPD整体随时间呈下降趋势，2013年的最大值和最小值分别为2.31kPa 和0.29 kPa，在2014年波动较大，最大值和最小值分别为3.30 kPa 和0.15 kPa；2013年生育期内总降水量P为374.5 mm，2014年生育期内总P为293.4mm。P在时间上有明显的分布不均现象，主要集中于6月下旬和7月下旬，有时会出现较大的日降水量。2013年的日P峰值出现在7月15日，为34.8 mm，2014年的日降水量峰值出现在8月28日，为24.6 mm；土壤含水率θ的变化趋势明显受到T和P的影响。每年7月T的水平较高，θ的值保持在一个较低的水平，在有降水时变化不明显。进入8月后T开始下降，此时当有降水后，θ的值会有明显上升。

图1 2013年和2014年玉米生育期内气象数据的逐日变化Fig.1 Daily variation of meteorological factors during growing season of maize in 2013 and 2014

2.2 实际冠层阻力的逐日变化及敏感程度分析

2013年、2014年的夏玉米全生育期的实际rc的逐日变化如图2所示。实际rc在夏玉米生育前期时最高，接近5000 s/m，之后明显减少，到6月下旬低于1000 s/m。结合图1 与图2，各气象因子和LAI的变化会导致rc的值变化，如6、7月的T、Rn在一个较高水平，此时rc值下降较快。每年6月下旬有较为集中的降水，此时rc出现明显的上升。

图2 实际冠层阻力的逐日变化Fig.2 Daily variation of measured canopy resistance

利用BP 神经网络模型分析rc对各影响因子的敏感程度，结果如表1所示。在保留全部影响因子的条件下，输出值与实测值间的相关性较高，R2=0.92，RMSE=346.61 s/m，d=0.90。当控制Rn输入后，R2=0.23，RMSE=1108.20 s/m，d=0.62，说明输出值和实际值的相关关系明显下降，rc对Rn的敏感程度较高；当控制T输入后，R2=0.92，RMSE=361.33 s/m，d=0.89；控制VPD输入后，R2=0.92，RMSE=350.40 s/m，d=0.90；控制θ输入后，R2=0.91，RMSE=374.41 s/m，d=0.89。可见rc对T、VPD和θ的敏感性程度都较低。当控制LAI输入后，R2=0.42，RMSE=959.76 s/m，d=0.63，rc对LAI的变化较为敏感。综合可知rc对这5 个因子的敏感程度为：Rn＞LAI＞θ＞T＞VPD。

表1 敏感程度分析结果Table 1 Results of sensitivity analysis

2.3 模型参数拟合结果

以2013年实测的涡度相关数据通过P-M 模型反算得到的实际rc为标准，分别采用最小二乘法和蚁群算法对JA 模型和CO 模型中的参数进行拟合，拟合结果见表2。2 种优化算法对参数拟合的效果较好，模型估计值和实测值的变化趋势相同，R2的值均超过0.8。总体上CO 模型比JA 模型的拟合效果更优，其中使用蚁群算法优化的 CO 模型拟合结果最好（R2=0.89，RMSE=410.90 s/m，d=0.88），参数拟合结果为b1=2.20，b2=0.18，b3=0.01，b4=-0.18，b5=-1.39，b6=-0.45。用蚁群算法优化的CO 模型在6月中旬以前的估计值比较接近实测值，之后的估计值和实测值之间有一定差距。用蚁群算法优化的JA 模型估计rc的效果略好（R2=0.83，RMSE=515.52 s/m，d=0.85），参数拟合值分别为rcmin=7.43 s/m，a1=360.80，a2=0.25，a3=-0.003。

将4 种方法所得参数分别带入对应模型，利用2014年的气象数据对夏玉米全生育期ET进行估计，并用涡度相关系统实测的ET值进行验证，结果见表3。在JA 模型中，采用最小二乘法得到的R2=0.65，RMSE=1.20 mm，d=0.76，采用蚁群算法得到的R2=0.67，RMSE=1.16 mm，d=0.76，二者对JA 模型的拟合效果均较好；在CO 模型中，采用最小二乘法得到的R2=0.70，RMSE=1.11 mm，d=0.75，采用蚁群算法得到的R2=0.72，RMSE=1.07 mm，d=0.75，因此在计算ET时，误差最小的是使用蚁群算法优化后的CO 模型，总体上CO 模型的计算结果要优于JA 模型，且在2 种模型中采用蚁群算法得到的结果精度更高。

3 讨论

研究表明，rc在夏玉米生育前期达到最大值，之后随着夏玉米生长逐渐下降，由于夏玉米生育期是5—9月，进入5月后，温度开始上升，日照时间延长，加速了水分的蒸发，同时玉米叶面积指数变大，气候变化和夏玉米生长使气孔开度变大，最终导致了rc值变小。rc对Rn的敏感程度最高，因为Rn是地表能量的重要组成部分，在地表能量平衡和水分平衡中有直接参与[50]；LAI对rc的影响来自于作物本身的生理参数变化；其他气候因素（θ，T和VPD）的变化对rc有影响，但是影响程度远不及Rn和LAI明显。这与Li 等[46]、Ding 等[52]和Zhang 等[53]研究结果一致。

优化后CO 模型在估算rc和ET时的精度优于JA模型，JA 模型对rc估算值偏小。CO 模型中不仅考虑了水分通过植物冠层时克服的阻力，还考虑了水分从土壤表面蒸发时克服的阻力[51]，无论对于稀疏还是密闭的植被都具有一定应用性。夏玉米生育前期，由于植被稀疏，土壤表面的水分蒸发量占夏玉米ET的绝大部分，因此忽略土壤阻力的作用将会带来一定的误差。JA 模型没有考虑土壤表面的阻力，因此估计的rc值偏小，可能导致估算的ET值偏大[29-30]。2 种优化方法中，蚁群算法的鲁棒性好，全局搜索能力强，能够收敛于全局最优解，在求解过程中能够不依赖人工调整，适用于优化求解问题[54]。而最小二乘法是线性估计，当变量间为明显的非线性关系时，最小二乘法得到的结果误差较大，因此最小二乘法在使用上局限性较大。不论是对参数进行优化，还是对2014年夏玉米的ET进行估算，蚁群算法的精度均优于最小二乘法，可成为参数优化的推荐算法。

由于数据质量的影响，本文仅选取了夏玉米2a生长期涡度相关数据和气象数据。在以后的研究中应选择时间跨度更长的数据，以提高夏玉米rc的参数拟合精度，进一步提高ET估计值的可靠度，为怀来地区的夏玉米作物需水量估计和农业水管理提供更为合理的科学参考。

4 结论

1）rc在玉米的生育期呈下降趋势，在生育前期rc值较大，接近5000 s/m。rc与各气象因子和LAI存在明显的非线性关系，采用BP 神经网络模型进行敏感分析，得到rc对各影响因子敏感程度为：Rn＞LAI＞θ＞T＞VPD。

2）最小二乘法和蚁群算法对rc的计算模型中参数拟合的效果较好，模型估计值和实测值的变化趋势相同。总体上CO 模型比JA 模型的拟合效果更优，其中使用蚁群算法优化的CO 模型拟合结果最好，R2=0.89，RMSE=410.90 s/m，d=0.88，且6月中旬以前的估计值比较接近实测值，之后的估计值和实测值之间有一定差距。

3）利用2014年气象数据，采用优化后的rc的计算模型计算2014年的玉米ET，CO 模型的计算ET值精度高于JA 模型，其中采用蚁群算法优化后的CO模型计算ET精度最高，R2=0.70，RMSE=31.44 mm，d=0.75，采用最小二乘法优化后的JA 模型计算ET的精度较低，R2=0.70，RMSE=31.44 mm，d=0.76。