基于温室内气象数据的适用于冬季日光温室番茄蒸腾量各子模型对比研究

王 淼，罗新兰，尹佳琪，苏 慧,常姝婷

(1.锦州市气象局,辽宁 锦州 121000; 2.沈阳农业大学农学院，辽宁 沈阳 110866;3. 辽阳市气象局，辽宁 辽阳，111000； 4. 葫芦岛市气象局， 辽宁 葫芦岛，125000)

随着水资源的匮乏及科学有效灌溉的需求，节水灌溉及作物需水量的准确模拟计算成为热门课题。例如在温室或大棚等果蔬生产环境下，对节约优质适量灌溉的需求更加迫切。

北方寒区日光温室冬季仅在中午前后通风，其生产环境相对封闭，蒸腾作用消耗潜热降低植株温度，释放水汽增加温室内湿度，室内湿度的增加又会一定程度地抑制植株蒸腾作用，增加植株病虫害发生频率。Penman公式是世界粮农组织推荐的一种计算作物需水量的方法[1]，Monteith在1963年对此改进，提出以能量平衡和水汽扩散为基础的Penman-Monteith方程(以下简称P-M方程)，适用于多种环境下作物蒸发蒸腾量的直接计算[2-3]。Moritlle等[4]试验验证了P-M方程可以用于温室环境下作物蒸腾量的计算。然而，由于影响蒸腾的因素较多且很多参数很难精准测定，尤其是不同环境下边界层空气动力学阻力(ra)、气孔平均阻力(rc)，目前还缺乏通用、精准、可用于不同类型温室的子模型计算式。大多数学者认为ra与植株叶片和温室内空气之间对流热交换形式密切相关，不同风速下应采用不同的对流热交换模型[1]。另外，温室内风速是影响ra的主要因素，大部分研究应用P-M方程计算作物蒸腾量时，测量微风速后直接带入半经验公式即求出ra值[5-19]。冬季室内空气流通不均匀，微风速测量难度很大时，也可采用不含风速项的模拟式来计算ra值[20-22]。温室内环境(室内总辐射Rg、饱和水汽压差VPD、室内空气温度Ta)是影响子模型rc大小的重要因素，常采用3种影响因素的半经验式或植株整个生长周期的特征值来模拟rc值[5,16-18，20]，鲜有利用P-M反演式来计算。国内外研究中常将ra、rc子模型直接应用于蒸腾量试验中，而少有对各子模型对比分析择优后带入P-M方程对植株蒸腾速率模拟计算。

本文分别选出4组适用于北方寒区冬季日光温室内的ra、rc子模型，试验证明不同子模型对P-M方程计算植株蒸腾速率有明显影响，通过初步筛选、检验各项指标择优的方式选出基于此环境下模拟植株蒸腾量准确率最高的ra、rc子模型。北方寒区日光温室内具有弱光、高湿、低温及低风速(甚至为0 m·s-1)的环境特点，同时，本文还对多株植株同时进行试验，这为应用P-M反演式计算rc值提供了条件，子模型ra、rc的选取对P-M方程的准确模拟及冬季日光温室园艺作物灌溉需水量的确定具有重要参考意义。

1 材料与方法

1.1 试验概况

试验于2016年11月—2019年2月在沈阳农业大学北山科研基地日光温室(123°5′E，41°8′N)内进行。温室方位南偏西7°。东西长度60 m，跨度7.5 m，脊高3.5 m，北墙高2.3 m，墙体厚45 cm，北墙外侧设有冷棚，试验期间揭帘时间8∶30，盖帘时间14∶30。

供试番茄品种为适于北方日光温室长季节栽培的圣罗兰。植株采用南北成行基质袋培，每行铺有一条滴灌管，滴头与植株相对应，滴头出水口插于栽培袋中间位置，晴天的9∶00、12∶00、14∶00 分别灌溉100 ml·株-1·次-1,基质水分控制在田间持水量80%以上，地面覆盖塑料薄膜(无裸露土壤)，行距1 m，株距20 cm，共定植920株，定植密度为3.38株·m-2。选择长势良好的3株番茄作为测试样本。温室顶部通风口可人工开启，晴天中午室内温度超过30℃时通风口打开10 cm进行通风，平均通风时间为15 min。其余时间处于密闭状态，试验期间植株处于成熟采摘期，茎粗1.41～1.48 cm、株高188.3～193.6 cm。

1.2 观测项目

室外及室内中部安装HOBO自动气象站以获取试验所需基本气象数据(辐射、温度、相对湿度、土壤含水量等)；番茄群体中间安装北京雨根公司的环境及茎流监测仪，测定项目包括：植株上方的太阳净辐射、植株的茎流值、空气温湿度、叶温(用与叶片相同颜色的绿色别针将热电偶与叶片固定)；采用热球风速仪测量植株上方微风速。试验期间24 h连续测试，取样间隔10 min，蒸腾速率估算中用到的相关气象要素观测项目及仪器位置如表1所示。基质含水量探针置于栽培袋中间位置。

表1 冬季日光温室内气象要素观测项目及仪器位置Table 1 Meteorological elements tested items and instrument position in winter solar greenhouse

单株植株实测蒸腾速率[4,9]计算公式：

(1)

式中，Tf为单株植株蒸腾速率(mm·h-1)；La为植株监测探针上方的叶面积(m2)；LAI为单株叶面积指数(植株叶面积/植株占地面积)；Fs为质量形式茎秆液流速率(g·h-1)，根据茎流计流速公式计算。

1.3 分析方法

1.3.1 叶面积指数(LAI) 叶面积采用主干上主叶长宽系数法计算，测量主叶的叶长(包括叶柄长度)和最大叶宽,计算公式[23]：

(2)

式中，LA为叶面积(cm2)；a为叶片长度(cm)；b为叶片宽度(cm)。LAI=LA/单株植株占地面积。

1.3.2 消光系数(K) 根据比尔-朗伯特(Beer-Lambert)定律得到：

(3)

式中，LAI为单株叶面积指数；I，I0分别为单株尺度的地面和植株上方的总辐射(W·m-2)。计算每日9∶00—14∶00每10 min的消光系数，根据冯新灵等[24]研究建议剔除不合理的消光系数数据(即K>1的数据)，再以日为单位将合理数据均值后作为每日的消光系数。

1.3.3 采用Penman-Monteith方程模拟日光温室基质袋培单株植株蒸腾量 在考虑土壤热通量、并结合适用于本环境的各边界层空气动力学阻力及气孔平均阻力计算式，Penman-Monteith方程计算温室内单株植株蒸腾量的基本式可转换[14,25-26]为:

(4)

联合国粮食及农业组织(简称FAO)推荐土壤热通量(G)白天采用0.5Rn，夜间采用0.1Rn[4]，本文根据北方寒区日光温室揭、盖帘时间对应的温室内日照时数对G进行修正，白天采用0.35Rn，夜间采用0.13Rn，已在前文对取值做了敏感性检验[11],本文不再赘述。

1.3.4 边界层空气动力学阻力模拟 于冬季日光温室内进行，仅在晴天中午室温超过30℃时通风(通风口开度10 cm，平均通风15 min)，适用于此环境下ra的4种计算方法如下:

(1)根据Campbell[22]研究结果，由叶片温度与周围空气温度差值引起的自由对流对ra起决定作用时，计算公式：

(5)

式中，ra1为单一叶片空气动力学阻力(s·m-1)；d为叶片特征长度(m)(本试验取值0.15 m)；Tl为叶片温度(℃)；Ta为周围空气温度(℃)；ra为单株植株边界层空气动力学阻力(s·m-1)；LAI为单株叶面积指数。

(2)根据康绍忠等[14]和刘贤赵等[19]研究结果，边界层空气动力学阻力是不同作物风速廓线的函数，由湍流传输阻力与叶片层流边界层阻力组成，可表示为：

ra=rg+rb

(6)

式中，ra为边界层空气动力学阻力(s·m-1)，rg为湍流传输阻力(s·m-1)，rb为叶片层流边界层阻力(s·m-1)。

根据Thom等[22]提出的适用于温室内计算湍流传输阻力的公式，当Uz=0时，公式简化为：

式中，Z为风速测量高度(m)，Z0为地面粗糙度(m)，d为零平面位移长度(m)。

根据康绍忠等[14]研究，对于一般作物：

Z0=0.13hc，d=0.64hc

式中，hc为番茄冠层高度(m)，其余符号意义同前。

根据此公式可避免由于温室内风速过低，出现湍流传输阻力无穷大的情况。

叶片层流边界层阻力rb与叶片特征宽度、风速有关, 当温室内风速不为零时，叶片层流边界层阻力不能忽略，此时rb可表示为：

式中，c为经验常数(0.01 m·s-10.5)，u为冠层上方风速(m·s-1)，w为叶片特征宽度(m)。

(3)根据Bailey等[1]关于温室内作物叶片的研究表明，边界层空气动力学阻力(ra)与叶面和空气之间的对流热交换系数(hp,a)的关系式为：

(7)

不同对流情况下，hp,a的计算式为：

其中

式中，Kc为空气的热导率(0.0267 W·m-1·K-1)，Gr为格拉晓夫数，Pr为普朗特数，β为空气的热膨胀系数(0.00375 K-1)，g为重力加速度(9.8 m·s-2)，Tl为叶面温度(℃)，Ta为室内空气温度(℃)，Re为雷诺数，v为空气运动粘性系数(1.6×10-5m2·s-1)，a为空气热扩散率(2.29×10-5m2·s-1)，u为室内空气流速(m·s-1)，l为作物叶片特征长度，番茄取0.15 m。

(4)Stanghellini[27]应用叶片与空气之间温度差值及温室内空气流速的经验公式计算ra值，公式可总结为：

ra=1174L0.5(L|Tl-Ta|+207u2)-0.25

(8)

式中，L为叶片特征长度(m)，Tl为叶片温度(℃)，Ta为室内空气温度(℃)，u为作物冠层上方微风速(m·s-1)。

1.3.5 气孔平均阻力模拟 作物通过叶片气孔与周围环境进行显热交换和水分传输，交换和传输过程中受气孔阻力控制，气孔阻力又受诸多气象因子影响，主要包括温室内植株上方总辐射(Rg)、温室内饱和水汽压差(VPD)、温室内空气温度(Ta)。本文对处于成熟采摘期的番茄植株假定其在大叶模式下进行气孔平均阻力的模拟，根据Boulard等[28]的研究表明，作物气孔平均阻力可由以下各统计方程来计算：

(1)根据Pollet等[15]的研究结果此方程可总结为：

×[1+0.016(Ta-16.4)2]

(9)

式中，Rg为温室内植株上方总辐射(W·m-2)，Ta为温室内空气温度(℃)，D为温室内饱和水汽压差(kPa)。

(2)北方寒区日光温室冬季VPD变化范围很小，根据Pollet等[20]及Boulard等[28]研究，表达式可总结为Rg、Ta的方程，即：

(10)

(3)Demrati等[29]认为冬季日光温室处于密闭状态下，可直接将室内总辐射带入下式计算出作物冠层对水汽的阻力(即气孔平均阻力)：

(11)

(4)由于影响温室内气孔平均阻力的主要因素包括温室内总辐射、温度、饱和水汽压差等气象因子，因此可以通过反演的方式计算出植株样本a的rc值作为植株样本b、c的rc值并带入P-M方程计算另两株植株的蒸腾速率模拟值，rc计算表达式如下：

(12)

式中，λE为植株a的实测蒸腾潜热值(W·m-2)，ra为植株a的边界层空气动力学阻力(s·m-1)，其余物理量为植株a单株尺度对应的数值且意义同前。

2 结果与分析

2.1 不同子模型模拟日光温室内各边界层空气动力学阻力(ra)与气孔平均阻力(rc)对比分析

基于模型模拟蒸腾速率的准确程度与ra、rc的取值有直接关系[2]，目前对于不同类型温室中ra、rc计算公式的选取仍没有一致的结论。为了验证本文提出的几种计算ra、rc的方法对P-M模型模拟精确程度的影响，将12月20日及12月21日的试验数据计算的各ra、rc值带入P-M方程计算逐时蒸腾速率，以对比各子模型在不同时刻的模拟效果。以1.3.4、1.3.5节中的公式为基础，计算单株植株蒸腾速率的相对变化率，白天各时刻蒸腾速率计算结果如图1所示。

由图1(a)可见不同的ra计算式对白天逐时蒸腾速率影响最大的时间是11∶00—13∶00，蒸腾速率平均变化量为8.01%，影响最小的时间为9∶00及15∶00左右，蒸腾速率平均变化量为0.76%，其余时间蒸腾速率平均变化量为5.11%；由图1(b)可见不同的rc计算式对白天逐时蒸腾速率影响最大的时间是11∶00—13∶00，蒸腾速率平均变化量为42.21%，影响最小的时间段为9∶00—10∶00，蒸腾速率平均变化量为27.48%。

由于北方寒区冬季严寒，日光温室于夜间覆盖保温被以保证室内温度满足植株正常生长需要，而此时蒸腾量很小。为了进一步明确子模型选取对植株蒸腾作用的显著影响，本试验对夜间番茄蒸腾速率作简要计算说明，对植株夜间的蒸腾速率值比较分析，给出不同ra、rc组合形式对夜间植株蒸腾量的影响情况，结果如表2所示，不同的ra计算式对夜间蒸腾速率影响较大，且影响程度相近，平均变化量为46.13%；各rc计算式对夜间蒸腾速率影响程度较一致，平均变化量为30.25%。

表2 基于不同子模型的冬季日光温室内夜间番茄各边界层空气动力学阻力与气孔平均阻力模拟的蒸腾速率/(mm·h-1)Table 2 Transpiration rate of different aerodynamic resistance (ra)、average stomatal resistance (rc)on different submodel in solar greenhouse during night winter

总体上看，不同的ra、rc对P-M方程模拟值的影响程度均是夜间大于白天，白天rc对P-M方程模拟值的影响程度明显大于ra，夜间则相反，因此在应用P-M方程计算植株蒸腾速率时子模型ra、rc的选用是非常重要的。

2.2 日光温室内单株尺度子模型运行结果

选取试验期间有典型天气的2017-12-13—2017-12-25(其中13日与17日为雾霾天，14日与24日为阴天，其余时间为晴天)连续13 d的测试值。根据1.3.4与1.3.5节各公式计算边界层空气动力学阻力(ra)与气孔平均阻力(rc)，各公式模拟值的对比见图2。

由图2(a)可以看出，基于1.3.4节此环境下公式Campbell(1977)与公式Stanghellini(1995)的变化趋势相同，波动范围不同，公式Campbell(1977)计算的ra平均日波动范围在晴天为185～489 s·m-1，阴天为249～357 s·m-1；公式Stanghellini(1995)计算的ra平均日波动范围在晴天为450～1 403 s·m-1，阴天为649～1 011 s·m-1。而公式刘贤赵(2001)、公式Bailey(1993)计算的ra值波动振幅过于平稳。

由图2(b)可以看出，冬季日光温室除午间通风时段，其余时间处于密闭状态，室内湿度一直维持较高值，VPD变化范围很小，考虑VPD时，rc的波动范围反而不那么明显，所以如公式Pollet(2000)所示仅考虑Rg与Ta时rc的计算值在全天均呈现明显的波动性。鉴于此，1.3.5节中公式Pollet(2000)及反演式计算的rc平均日波动范围在晴天分别为261～1 359、253～1 356 s·m-1；阴天分别为250～1 291、235～1 260 s·m-1。此环境下考虑影响rc的3个因素(植株上方总辐射,Rg;饱和水汽压差,VPD;室内平均空气温度,Ta)及一个因素(Rg)时在白天时段rc模拟值基本相同，这与实际情况不符。

由于ra、rc变化具有波动性，由图2分析结果初步确定1.3.4节中由空气温度、叶温及植株生理指标(叶片特征长度、LAI)代入公式Campbell(1977)、公式Stanghellini(1995)计算的ra子模型、1.3.5节中根据室内植株上方总辐射(Rg)、空气温度(Ta)代入公式Pollet(2000)及P-M方程反演式计算出的rc模拟值符合其随时间波动的特征。

2.3 不同子模型模拟日光温室内边界层空气动力学阻力(ra)与气孔平均阻力(rc)的验证与筛选

为了进一步说明子模型对P-M方程模拟植株蒸腾速率的影响程度以及模型的通用性，分别将ra于12月18日—1月1日及rc于12月11日—12月25日的子模型计算值代入P-M方程并求得蒸腾速率模拟值与实测值的Pearson相关系数、平均相对误差(MRE)、平均绝对误差(MAE)、均方根误差(RMSE)及模型效率指数(EF)[30-31]，结果见表3。由表3分析得出，北方寒区冬季日光温室生产环境中各ra子模型的模拟值与实测值相关性最强的是Campbell与Bailey模拟式，其与实测值对比的相关系数均为0.826(刘贤赵0.813；Stanghellini 0.822)；平均绝对误差(MAE)Campbell模型为0.02，其余3个子模型均为0.03；Campbell模型的均方根误差(RMSE)最小；模型效率指数(EF)Campbell>Bailey=Stanghellini>刘贤赵。因为本试验测试时间为冬季，温室中午前后开窗通风时室内风速垂直流通且并不均匀，其余时间处于密闭状态，微风速的测量偏差较大。通过各检验系数结合子模型整体波动性，综合分析得出选用包含空气温度、植株叶片温度及植株生理指标(叶片特征长度、单株植株叶面积指数LAI)等参数模拟单株植株边界层空气动力学阻力值；各rc子模型计算得出的模拟值中反演式与实测值对比的相关系数为0.955(Pollet 0.861；Boulard 0.863，Demrati 0.860)，反演式的平均相对误差(MRE)、平均绝对误差(MAE)、均方根误差(RMSE)均低于其余三者，反演式的模型效率指数(EF)为0.95，其准确率高于其余三者，所以本文选用反演式模拟单株植株气孔平均阻力值。

表3 各子模型模拟的植株蒸腾速率模拟值与实测值的分析检验Table 3 Inspection analysis between simulated on different submodel and measured of transpiration rate

3 讨论与结论

3.1 讨论

基于模型模拟蒸腾速率的准确程度与边界层空气动力学阻力、气孔平均阻力的取值有直接关系[2]，而二者数值的精确获得很难。冠层内外风速及温度的变化会产生不同的气流模式类型[1],气孔平均阻力的取值不仅与太阳辐射、VPD、空气温度、作物水分等有关[32],且还受冠层结构、冠层上方气流及土壤条件等影响[12]。已有部分研究于空气流通较均匀的塑料大棚[5]、玻璃温室[6]、防虫遮阳网温室[7]、日光温室[8-13,24]中直接采用运用微风速计算的ra值、半经验公式总结的rc值或总结出的特征值模拟整个试验周期的子模型运行结果，本试验验证的4组ra、rc计算式对P-M方程模拟植株蒸腾速率有明显影响，验证各组子模型准确率后再运用其计算植株蒸腾量值是非常必要的。

Bailey等[1]在英国(5月)及西班牙(7月)空气流通较均匀的温室内利用叶片边界层空气流通模式模拟ra值。结果表明：开窗通风时，强迫对流模拟式更接近实测值，两地区的模拟值与实测值平均差值分别为2.8%、3.2%；康绍忠等[14]和刘贤赵等[19]在Thom等[22]提出的温室内湍流传输阻力计算式基础上提出叶片层流边界层阻力值，提高了冬季大棚作物需水量计算的准确率；Stanghellini[27]利用叶片与空气温差及室内空气微风速的半经验公式也得到了较好的模拟效果；Campbell等[25]总结出基于叶温、空气温度、植株生理指标(LAI、叶片特征长度)为基础的模拟式，避免了由于冬季温室内空气流通不均匀时误差较大的问题，经检验在冬季日光温室空气流通不均匀时采用此模拟式可提高试验准确率。经对比分析，本文选用Campbell等[25]总结出的ra模拟式，并将LAI换算为单株尺度，当室内微风速测量困难时可采用不含有风速项的ra模拟式。

国内外相关研究中多是利用温室内影响气孔平均阻力的气象因子(Rg、VPD、Ta)总结的半经验公式来模拟rc值[20,28-29]，本文利用反演式求得rc值，经验证利用反演式模拟植株蒸腾速率的误差较小，且Pearson相关系数达0.955，准确率明显提高。

本试验在北方寒区很少通风的日光温室内进行，此条件下空气流速不均匀，微风速测量难度较大。本文仅对冬季日光温室内ra、rc筛选择优，并未在夏季日光温室内验证。在今后试验中将增加夏季空气流通较均匀时各子模型的筛选验证。试验番茄为基质袋培，地面覆盖塑料薄膜，无裸露土壤。番茄根部水热状况与盆栽或土壤栽培有区别，各组子模型的适用性对蒸腾结果的影响有待进一步研究及优化。

总体来看，本文主体采用应用较广泛的P-M方程模拟温室单株植株蒸腾量，对边界层空气动力学阻力及气孔平均阻力不同子模型模拟式作对比分析；本试验选取的ra、rc子模型模拟植株蒸腾速率的准确率较高，特别是运用反演式计算的rc值达到了较优的模拟效果。试验对比分析结果可为北方寒区冬季日光温室空气流通并不均匀(微风速测量难度较大)、选取多株植株同时进行试验时子模型的选取提供参考。

3.2 结 论

(1)总体来看，ra、rc随时间波动变化，利用不同ra、rc模拟式计算北方寒区日光温室蒸腾速率的相对变化率有差别。不同子模型模拟的ra在白天、夜间平均变化量分别为5.73%、46.13%；不同子模型模拟的rc在白天、夜间平均变化量分别为37.10%、30.25%，rc模拟式对白天蒸腾速率的影响程度明显大于ra的选择，虽然夜间不同的ra计算式对蒸腾速率模拟值影响程度大于rc约15.88%，但由于北方寒区日光温室冬季番茄蒸腾作用主要集中于白天，因此模拟植株蒸腾速率时子模型rc的选择尤为关键。

(2)本试验条件下，选取的子模型计算得出的晴天、阴天番茄植株边界层空气动力学阻力(ra)的变化范围分别为185～489、249～357 s·m-1，气孔平均阻力rc变化范围分别为253～1 356、235～1 260 s·m-1，应用P-M方程模拟日光温室作物蒸腾量时应根据温室内特定的气象环境选择合适的ra、rc模拟式。

(3)本研究验证了ra、rc(各4组)子模型组合式分别代入P-M方程模拟植株蒸腾速率并与实测值对比分析，结果表明：Campbell模型的Pearson相关系数、MRE、MAE、RMSE、EF指标均优于其余3种模型，又结合子模型的波动特征，综合考虑本试验的ra值选取Campbell模型模拟；反演式的各指标均优于其他模拟式且其Pearson相关系数为0.955、MRE为7.16%、MAE为0.01、RMSE为8.56、EF为0.95，因此，本文选取反演式计算rc值。建议在温室内微风速测量难度较大时，可采用测试叶温、空气温度、LAI及叶片特征长度的方式来计算ra值；由于rc值不受植株自身条件影响，仅与室内环境要素有关，因此当试验同时选取多株植株时，可以通过其中某一植株的反演式计算出的rc值同样作为其余植株的rc值代入P-M方程计算相应单株植株的蒸腾量可达到较优的模拟效果。