基于Bayesian的Penman-Monteith模型简易大棚蒸散模拟

李 群 汪 星 汪有科，3 纳文娟 唐 燕 吴久江

(1.西北农林科技大学水利与建筑工程学院， 陕西杨凌 712100； 2.宁夏大学农学院， 银川 752201；3.西北农林科技大学水土保持研究所， 陕西杨凌 712100)

0 引言

大量研究表明Penman-Monteith(PM)模型可以应用在温室环境蒸散估算[1-2]，并且目前是温室蒸散模拟应用最广泛的机理模型[3-4]。当前国内外学者对温室PM蒸散模拟的研究多以智能自控温室为载体，简易大棚应用PM模型估算蒸散量少有报道。我国90%以上的温室属于简易大棚[5-6]，如小跨度塑料大棚的蒸散模拟鲜有报道[7-9]。由于温室环境与大田环境相比风速明显更低，PM模型的关键参数空气动力学阻力ra使用常用的Perrier公式往往会被过高估计，进而导致PM模型的模拟效果较差[10-13]。为解决这一问题，BOULARD等[14]提出一种基于温室外部环境估算温室蒸散的方法，李军等[7]在我国长江下游单栋塑料大棚中使用此方法模拟塑料大棚小白菜蒸散结果较好。但这种方法要求温室空气与室外环境耦合紧密，即需长期保持温室自然通风，不适合在我国北方地区冬季使用。陈新明等[15]与王健等[12]在日光温室环境中将风速取0并使用THOM等[16]的公式计算ra，这种方法有效解决了温室低风速条件下ra被过高估计的问题，但对于大棚有通风情况下的蒸散估算仍未涉及。孙怀卫等[8]建立了华东地区塑料大棚内部风速与外部风速的线性关系，并基于能量平衡和紊流扩散理论建立了韭菜蒸腾模型。这种方法较为繁琐，需要使用计算流体动力学(Computational fluid dynamics，CFD)分析棚内气流场特征。此外，上述研究使用的PM模型其约束函数形式均是唯一的，这会导致模型模拟精度缺乏比较，也没有考虑模型的不确定性。

PM模型结构较为复杂，子模型涉及的待率定参数较多，模型的不确定性可能较强，基于上述考虑选择使用Bayesian方法作为模型参数率定方法。Bayesian方法将参数的概率分布、实测数据、参数波动范围和不确定性的相关假设相结合，既考虑了模型参数、输入变量与输出结果的不确定性，也考虑到了已知信息以及参数之间的互相影响[17]，是一种适合高维模型参数率定的优秀方法[18-19]。近20年来，Bayesian方法在大田环境的蒸腾、蒸散模拟及模型不确定性分析中都得到了较好的应用[18-21]，但在温室蒸散模拟研究中很少使用Bayesian方法，对温室蒸散模型不确定性的分析更是鲜有报道。

本研究在关中地区简易大棚中使用温室中应用最为广泛的PM蒸散机理模型，进行简易大棚日尺度下章姬草莓蒸散模拟研究，以期为简易大棚蒸散模拟提供参考。

1 材料与方法

1.1 研究区概况

试验区位于陕西省渭南市白水县雷牙镇刘家卓村(35°12′N，109°36′E，海拔809 m)，年平均太阳总辐射5 363.52 MJ/m2，年平均气温11.4℃，年平均降水量577.8 mm。试验地0～0.3 m土层土壤质地均一，土壤渗透性良好，砂砾、粉粒和黏粒3种粒级质量分数分别为70.1%、21.8%、8.1%，属于沙质土。0～0.3 m土层土壤容重为1.39 g/cm3，环刀法测得田间持水率为35%，计划湿润层深度取章姬草莓根系平均深度0.3 m。定植时间为2018年8月31日傍晚，定植期共12 d。使用宽窄行起垄方式栽植，9垄，下宽0.45 m，上宽0.3 m，垄高0.35 m，间距0.25 m，每垄660株，栽植2行，行距0.3 m，株距0.2 m，每公顷约栽90 000株。栽植时同一行植株的花序朝同一方向，使草莓苗弓背朝花序预定生长方向，苗心露出畦面，根系平展埋入疏松土层。草莓苗栽植后及时灌溉定植水，灌至0.3 m深度土层使土壤含水率达到田间持水率，定植水每2 d灌溉一次，持续13 d。使用国产迷宫式滴灌带，滴灌带内径16 mm，额定流量2 L/h，滴孔间距200 mm，湿润比取0.35，灌溉水利用效率取0.9。使用黑色PVE地膜以保温保墒，定植后全地膜覆盖栽培。

1.2 数据采集与处理

试验区简易大棚中央距离地面1 m处设有小型气象站(RR-9100型)自动监测记录气象数据(步长1 h)，测量指标包括空气温度、空气相对湿度、风速、太阳辐射，试验数据采集时间段为2018年8月31日至2020年3月15日，包含章姬草莓两个完整生育阶段。

章姬草莓冠层高度和叶面积指数每5 d测量一次，待测样本植株每次随机选取50株取平均值。测量所选取植株所有有效长度和有效宽度，根据实测折算系数计算每株草莓的叶面积并取平均值，最后根据种植密度得到叶面积指数。实测蒸散量根据水量平衡原理，采用称量法测量草莓植株蒸散量。使用电子天平对桶栽草莓进行连续称量，试验桶内径0.5 m，高0.5 m，表面覆盖和棚内草莓相同的地膜。桶栽草莓共6桶，均匀分布在一个简易大棚中。实测蒸散量为6桶草莓蒸散量的平均值。

模型精度评价使用平均相对误差(MAE)、决定系数(R2)和威尔莫特一致性指数(D)，模型精度可靠的标准为R2和D均大于0.8且MAE小于20%。使用Bayesian方法估算PM模型中参数的后验分布，并分析参数的不确定性。

1.3 PM模型

PM模型是一种基于能量平衡原理的蒸散机理模型，是目前温室蒸散模拟应用最广泛的机理模型[3-4]，PM模型一般形式为

(1)

式中λET——潜热通量，W/m2

Δ——饱和水汽压-温度关系曲线在平均气温T处的切线斜率，kPa/K

Rn——净辐射通量，W/m2

G——土壤热通量，W/m2

ρa——空气密度，取1.29 kg/m3

Cp——空气定压比热容，取1 004 J/(kg·K)

es——饱和水汽压，kPa

ea——实际水汽压，kPa

rc——冠层边界层阻力，s/m

γ——湿度计常数，取0.067 kPa/K

ra——空气动力学阻力，s/m

其中，ra和rc是两个关键参数。

1.3.1空气动力学阻力模型

1.3.1.1Perrier公式法

ra模型包括PERRIER[22]提出的Perrier公式，此公式已在通风条件较好的大田环境下得到普遍认可[18,23]。多数研究认为该公式在温室低风速环境中并不适用[10-13]。最终Perrier公式可以表达为仅与冠层高度hc和参考高度处风速uz相关，计算式为

(2)

1.3.1.2Thom风速零值法

(3)

1.3.1.3热传输系数法

温室内水汽以涡流方式扩散，ra可以通过热传输系数计算[26]。温室的对流类型可分为自由对流、强迫对流和混合对流，计算热传输系数hs时先判定出对流类型，之后根据不同对流类型的公式计算对应的hs[27-28]。表1列出了热传输系数法中3种对流类型出现的比例，以及各类对流类型条件下与hs相关的因子。由表1可知，温室通风时绝大部分时间温室内的对流都属于强迫对流，强迫对流条件下冠层高度风速u与hs正相关，叶片特征长度dc与hs负相关；不通风时绝大部分时间属于混合对流，冠层温度与空气温度的差值Δt、u分别与hs正相关。

表1 3种对流出现比例以及与hs相关的因子Tab.1 Proportion of three kinds of convection and hs related factors

将3类对流条件下计算出的热传输系数hs代入可计算得到ra[26]，即

(4)

式中LAI——叶面积指数

1.3.2冠层边界层阻力模型

1.3.2.1基于叶片气孔阻力的冠层边界层阻力模型

BAILEY等[26]通过分析叶片气孔阻力rs与气象因子的响应，构建了基于有效叶面积指数LAIe的冠层边界层阻力rc模型，计算式为

(5)

(6)

rs是确定rc的关键参数，由气象因子、水分条件所决定[29-30]。本研究简易大棚章姬草莓实施充分灌溉，对rs与简易大棚内气象因子的相关关系分析结果显示，太阳辐射Rs为影响rs的主要气象因子，这一结果与文献[31-32,10]对温室番茄和甜椒的研究结果相似。图1为最终通过非线性拟合得出的Rs与rs的关系，其关系式为

图1 简易大棚章姬草莓rs与Rs的关系Fig.1 Relationship between rs and Rs of Akihime strawberry in basic greenhouse

(7)

综合式(5)～(7)，最终构建出包含气象因子Rs和作物LAIe的双因子rc模型。

1.3.2.2基于冠层导度的冠层边界层阻力模型

冠层边界层阻力rc是冠层导度gc的倒数，通过构建气象因子、土壤含水率与冠层导度gc的关系

(8)

式中gmax——冠层顶部叶片最大气孔导度，m/s

KQ——短波辐射消光系数

Qh——冠层顶部可见光通量密度，W/m2

Q50——冠层导度达到最大值一半时的可见光通量，W/m2

VPD——空气饱和水汽压差，kPa

D50——冠层导度达到最大值一半时的饱和水汽压差，kPa

F2——水分胁迫系数

现场观摩结束后，新洋丰农艺师在室内会议室，为参会人员详细地讲解茄子高产高效栽培技术要点及洋丰百倍邦套餐肥的优势与特点，农户们听后激动不已。海南荆岛公司张斌经理发布钜惠订货政策后，农户们一算账，用洋丰百倍帮套餐肥不仅茄子长势比用进口肥更好，而且按活动政策现场订肥，一亩地投入还要省几十甚至上百元。

求倒数后即可确立冠层边界层阻力rc2模型[19,33]。本研究中不存在水分胁迫，因此F2取值恒为1。此模型共4个待定参数，分别为gmax、KQ、Q50、D50。

1.3.2.3基于叶片气孔导度尺度提升的冠层边界层阻力模型

冠层边界层阻力rc与冠层导度gc是倒数关系，因此构建基于气孔导度gs的rc模型的一个基本思路是将gs尺度转换到gc，两种常见的gs模型是Jarvis模型[34]和Ball模型[35]。本研究参照张宝忠等[36]的方法，基于gs实测值建立gs与气象因子(光合有效辐射PAR、VPD和空气温度Ta[37-43])的Jarvis模型，再以光合有效辐射PAR为尺度转换因子，构建出基于气象因子PAR、VPD和Ta的冠层导度gc模型。

gc模型可表示为gs在叶面积指数上的积分，即

(9)

通过PARa作为尺度转换因子得到gc模型为

(10)

式中a1、a2、a3、a4——待率定经验参数

PARa——叶片截获的光合有效辐射，μmol/(m2·s)

K——消光系数

PARh——冠层顶部高度的光合有效辐射，μmol/(m2·s)

最后求倒即可得到rc3模型，此模型共5个待定参数，分别为K、a1、a2、a3、a4。

1.3.3PM模型汇总

将上述PM模型中2个关键参数ra和rc的模型两两组合，共得到9种模型(表2)。其中涉及rc2和rc3的6种模型(*标记)使用Bayesian方法进行参数估计。

表2 PM模型汇总Tab.2 Penman-Monteith models summary

1.4 Bayesian参数估计方法

PM模型结构较为复杂，子模型rc2和rc3涉及的待率定参数较多(分别为4个和5个)，模型的不确定性较强，基于上述考虑选择使用Bayesian方法作为模型参数率定方法。

子模型rc2中的参数gmax、KQ、Q50、D50和子模型rc3中的参数K、a1、a2、a3、a4作为参数向量β的组成部分，先验信息见表3。

表3 PM模型参数先验分布Tab.3 Prior distribution of parameters in Penman-Monteith models

假设模型模拟值S与实测值O之间的差异即模型误差独立且服从方差σ2为未知常数的正态分布，模型误差均值为0[47-48]，则包含n个观测值的实测数列似然函数为[19]

(11)

式中Oi——第i个蒸散实测值

xi——模型输入向量

S(xi,β)——参数向量β下模型输出值

σ——模型误差的标准差

后验分布抽样使用基于Gibbs抽样的马尔科夫-蒙特卡罗(MCMC)方法，使用WinBUGS软件实现，运行3条平行MCMC链，迭代次数30万次并酌情增加。

2 结果与讨论

2.1 简易大棚日尺度空气动力学阻力ra模型筛选

由表4可知，ra模型中ra1表现最差，ra2与ra3表现相当；分析ra1-rc3、ra2-rc3、ra3-rc3，可得出相同的结论。ra1表现差的原因是ra1模型中风速在分母，温室环境中日均风速较低，因此会高估ra从而影响最终模型的估计效果，因此温室使用ra1结果较好的报道很少[49]，一些研究甚至建议不在温室中使用这种方法[10-13]。但本研究中ra1-rc3模型完全符合3个精度指标，ra1-rc2模型仅率定年D指标稍微超出要求。因此，简易大棚日尺度条件下虽然ra1较ra2和ra3表现差，但并非不可用，这一结论有别于现有温室环境的报道。Perrier公式法能在研究区使用的原因是棚栽草莓通风时间较长，另外简易大棚的密闭性较高端温室差，存在漏风情况，这些原因导致研究区简易大棚内部日均风速较一般温室高，进而使Perrier公式计算结果不会被过高估计。本研究统计的棚内平均风速为0.36 m/s(2018率定年与2019检验年2年平均值)，接近已有文献中日光温室风速日均值上限(0～0.36 m/s)[50]，使用Perrier方法计算出的ra1为253.41，与已有的ra报道差别不大[29]。但由于ra1模型与ra2模型的复杂程度相当且远小于ra3，ra2模型与ra3模型精度相当，因此从模型精度与复杂程度两个角度考虑，研究区简易大棚日尺度条件下推荐使用ra2模型。

表4 PM模型蒸散模拟精度评价Tab.4 Accuracy evaluation of evapotranspiration simulation in Penman-Monteith models

2.2 简易大棚日尺度冠层边界层阻力rc模型筛选

简易大棚日尺度蒸散量模拟值与实测值相关关系如图2所示，可直观看出除ra1-rc1和ra2-rc1拟合效果较差以外，其余各PM蒸散模型的日蒸散量模拟值与实测值能较好地沿1∶1直线分布。ra3-rc1虽然能较好沿直线1∶1分布，但散点明显更分散，决定系数R2也较低(0.75)。分析PM蒸散模型精度评价(表4)，发现包含rc1子模型的PM蒸散模型表现最差。可能的原因是rc1子模型考虑的作物-环境因子不全面，即模型缺少关键的作物-环境因子，进而导致一些关键的生理-物理过程被忽略，最终影响模型精度。rc1涉及的因子仅有太阳有效辐射和叶面积指数，rc2较rc1增加了描述冠层结构的参数项并增加了棚内水汽压亏缺这一环境因子，rc3较rc1增加了温度、饱和水汽压亏缺等环境因子。rc1较rc2、rc3缺少的这些因子都被证明是估算冠层导度的关键因子[51]。

图2 简易大棚日尺度蒸散量模拟值与实测值相关关系Fig.2 Correlations between simulated values of daily evapotranspiration and measured values in basic greenhouse

分析表4中rc子模型rc2和rc3，发现rc3较rc2性能更好(ra2-rc3优于ra2-rc2，ra3-rc3优于ra3-rc2)。对比rc3与rc2两个子模型，发现rc3较rc2多考虑了Ta这一环境因子，但子模型rc3优于rc2的原因并不一定是因为Javrvis模型中rc3多考虑了Ta因子。CHEN等[52]通过对包含不同数量关键环境因子(1～4个)、不同约束函数形式共计191个Javrvis模型进行了蒸腾模拟研究，发现Javrvis模型中嵌套增加关键环境因子并不总是能提高模型性能，模型性能还取决于Javrvis模型的约束函数形式。因此，rc3精度更高的一个可能原因是本文在组建Javrvis模型时所选择的关键环境因子以及其约束函数形式综合起来更加适合研究区简易大棚的环境。此外，rc3子模型较rc2子模型包含有更多的待率定经验参数，与之对应的是包含有更多的先验分布信息，这意味着Bayesian分析方法在参数率定的过程中可以利用更多的先验分布信息，这对减小模型不确定性、提高模型性能具有积极意义。因此，rc3精度更高的另一个可能原因是Bayesian参数率定过程中先验信息更加充分。

rc3子模型中以光合有效辐射PARa为尺度转换因子，以光合有效辐射PARa和饱和水汽压差VPD构建的冠层导度估算模型被证明能较好地实现叶片气孔导度向冠层导度的尺度转化提升[38]，但此模型中包含有难以测量的因子——消光系数K，K在测量过程中容易产生测量误差，误差累积到达一定数量必然会影响模型精度。因此，当模型中包含难以实测的因子时，在先验信息充分的条件下建议尝试将此因子当作未知参数，结合Bayesian方法计算出这些难以实测因子的估计值。

2.3 简易大棚日尺度不同生育阶段的参数不确定性分析

同时满足3个精度指标的PM模型共有5个(ra1-rc3、ra2-rc2、ra2-rc3、ra3-rc2和ra3-rc3)，图3所示5个模型中涉及的大部分参数的95%置信度后验分布区间在不同生育阶段与先验分布区间相比明显缩小，这说明Bayesian方法在参数率定过程中减小了模型参数的不确定性，参数率定效果较好。某些参数仅在特定的生育阶段率定效果较好，在膨果期中参数D50、KQ、Q50和gmax率定效果均较好。在模型ra2-rc2和ra3-rc2中，参数D50的后验分布区间较先验分布区间分别缩小了96.44%和93.98%，参数KQ的后验分布区间较先验分布区间分别缩小了56.08%和24.15%，参数Q50的后验分布区间较先验分布区间分别缩小了97.78%和99.50%，参数gmax的后验分布区间较先验分布区间分别缩小了99.75%和99.68%。仅特定生育阶段参数率定效果好的原因是先验信息不够充分[53]。本研究参照已有文献中参数的先验信息均没有区分生育阶段(表3)，这是造成先验信息不充分的一个原因。例如叶面积LAI在草莓生长初期随着植株的快速生长叶面积增大较快，花期时叶面积才趋于稳定，如果使用全生育期LAI均值反推出的参数作为先验信息，必然会影响到Bayesian参数率定结果。本研究中参数gmax随生育期变化明显，在草莓生长初期参数值较大，在膨果期骤降至极小值，这一变化规律与青藏高原高山草甸[19]、西北地区胡杨林[54]和美国威斯康星州草地[55]相同。因此，D50、KQ、Q50和gmax这类参数与植株生理特性和环境辐射关系密切，对作物生育阶段的变化比较敏感，使用区分生育阶段的先验信息进行参数率定效果会更好。

图3 2018年PM蒸散模型日尺度不同生育阶段参数最优估计值和95%置信区间Fig.3 Optimal estimation and 95% confidence interval of parameters in different growth stages of PM models in 2018

参数a2率定效果最好，在草莓全生育阶段中参数a2后验分布区间都较先验分布区间明显缩小。ra1-rc3模型中参数a2的后验分布区间较先验分布区间在4个生育阶段分别缩小了72.32%、88.59%、56.10%和86.32%，ra2-rc3模型分别缩小了97.65%、92.38%、93.31%和98.24%，ra3-rc3模型分别缩小了97.65%、92.38%、93.31%和98.24%。参照对D50、KQ、Q50和gmax这类参数有关生育阶段敏感性的分析，同样使用不够充分的先验信息却得到了在全生育阶段都较好的参数率定结果，说明参数a2对作物生育阶段的变化不敏感。此外，还可以从rc3模型分析得出相同的结论。本研究使用的rc3模型中参数a2仅与VPD相关，VPD反映了空气的干燥程度，很明显这一环境因子在简易大棚环境中不随生育阶段发生显著变化。

综上，结合模型模拟精度评价与模型参数不确定性分析，筛选出ra2-rc3和ra3-rc3两种适宜研究区简易大棚日尺度下的草莓PM蒸散模型，由于ra3-rc3模型涉及的热传输系数法计算繁琐，故推荐优先使用ra2-rc3模型。

3 结论

(1)Bayesian参数估计方法使PM模型在简易大棚日尺度蒸散估算精度较高，应用Bayesian参数估计方法的6种PM模型中有5种在模型率定年和模型检验年同时满足3个精度指标。

(2)Bayesian参数估计方法能有效减小简易大棚日尺度PM模型中部分参数的不确定性，在模型ra2-rc3中参数a2的后验分布区间较先验分布区间在4个生育阶段分别缩小了97.65%、92.38%、93.31%和98.24%，在模型ra2-rc2中参数D50、KQ、Q50和gmax在膨果期后验分布区间较先验分布区间分别缩小了96.44%、56.08%、97.78%和99.75%。

(3)筛选出ra2-rc3和ra3-rc32种适宜简易大棚日尺度蒸散模拟的PM模型，其中最优模型为ra2-rc3。