玉米蒸散发变化的归因—基于贝叶斯原理参数优化的双源模型

2022-11-30赵子敬魏国孝刘红娟田强龙邱中齐王福兵杨泽伟

节水灌溉 2022年11期

赵子敬，魏国孝，刘红娟，田强龙，邱中齐，王福兵，杨泽伟

（兰州大学资源环境学院，兰州 730000）

0 引言

几乎所有的陆地-大气耦合过程都受到陆地水储存和通量的显着调节，这些储存和通量也会影响碳、养分、生态系统结构和功能，并限制社会的淡水供应[1]。蒸散发（ET）包括土壤蒸发（ETs）和植物散发（ETc）是气候系统的中心过程，是水、能量、循环的纽带[2]。ET的准确评估对有效灌溉和实现水资源精细化管理有重要作用[3]。ET估算的准确性与模型的选择和参数有关，并且取决于对气候和生物因素的准确评估。因此，了解气候和植被对ET的相对影响对于预测水循环将如何响应未来的气候和生物变化具有重要意义[4]。

ET受多种因素控制，有研究表明，ET的变化是由太阳辐射、空气温度、水汽压、风速等变量驱动的[5-7]。不同气候区域对ET影响最大的气候因素不同，影响ET过程的不同因素之间复杂的相互依赖关系使得准确量化ET具有挑战性[8]。有些研究对蒸散发变化趋势主要集中在气候因素上，而没有考虑植被的作用。然而，植被在蒸散发中起着至关重要的作用，尤其是通过植物控制ETc的过程。目前有很多方法可以分离气候和植物驱动因素对ET的影响[4,9,10]。Stoy等[4]将涡动相关方差测量与线性扰动分析方法相结合，分离出美国东南部3个生态系统中物理和生物因素对ET的相对贡献。YANG 等[10]对PenPan 蒸发模型进行了修正，推导出一阶偏微分形式的蒸发皿蒸发方程，以量化气候因子的贡献。LI 等[9]利用Penman-Monteith模型、修正的作物系数法、Priestly-Taylor 模型、线性回归模型解释了2007-2013年中国西北玉米地ET年际变化。在该研究中，他们证实了植被驱动因子在调控ET过程中起着重要作用，强调准确估算冠层导度在ET建模和预测中的重要性，并为区分气候和植被对ET变化的贡献提供了一种新方法。

为了明确考虑通量相互作用，Shuttleworth 和Wallace[11]提出了一个双源模型(以下称为“SW 模型”)，该模型结合了两个单源模型:一个是土壤蒸发(ETs)，另一个是植物散发(ETc)。SW 模型旨在捕捉稀疏植被条件下叶片下的裸土状况，考虑冠层和土壤对总能量的贡献。在模型中，植物冠层阻力(rcs)和土壤阻力(rss)和冠层高度与参考水平之间的空气动力阻力(raa)以及叶片到冠层高度和土壤表面到冠层高度之间的空气动力阻力(rca和rsa)分别通过控制地表与大气之间的能量传递来调节植物的能量传递[12]。然而，SW 模型中复杂参数和模型结构的不确定使得SW 模型模拟ET具有挑战性。近年来，在开发用于驱动和处理参数和模型结构不确定性的系统框架方面取得了重大进展[13,14]。特别是，一种新的自适应马尔可夫链蒙特卡罗(MCMC)算法，又称为差分进化自适应都市(Differential Evolution Adaptive Metropolis )DREAM 算法，被用来在贝叶斯框架内有效估计驱动数据和参数不确定性[15,16]。所以根据贝叶斯原理利用DREAM 算法优化SW 模型参数，这不仅简化了参数的获取，还可以降低参数的不确定性。偏导数方法可以量化驱动因子对ET变化的相对贡献，所以本研究将SW 模型中的ETs和ETc分开，定量评估了冠层上方和土壤表面的可用能量(A和As)、饱和水汽压差(es-ea)和阻力raa、rsa、rca、rcs、rss在ET变化中的贡献，从而确定ET变化的主要驱动因子。

目前，利用SW 模型估算ET时很少考虑阻力的变化，将ET的年际变化归因于这些阻力的研究也很少受到关注。本研究旨在探讨大满站玉米生态系统蒸散发变化的潜在驱动因素，利用偏导数方法量化气候因素和植被因素对ETs和ETc变化的相对贡献。而阻力因子主要受到气候和植物因素的影响。所以考虑了阻力驱动因子的变化，并对SW 模拟的ETs和ETc的变化进行归因分析。这一结果为了解玉米ET的年际变化提供了有用的见解，并强调了阻力在调控ET中的重要作用。

1 材料与方法

1.1 研究区概况

为了评估不同驱动因子对ETs和ETc年际变化的影响，本研究选择了2016-2018年7-9月份大满通量站的数据。玉米生长试验是在中国西北部甘肃省张掖市的大满站进行的。大满绿洲位于黑河流域中游，是中国西北干旱区第二大内陆河流域。黑河流域中游地区以绿洲灌溉农业为特征，是生活用水和农业用水量大的地区。它的年平均降水量和气温分别为125 mm 和7.2 °C（1960-2000年），年平均潜在蒸发量约为2 290 mm，年平均日照时数3 106 h，无霜期148 d。土壤类型以粉质黏土为主，冻土层深度约143 mm。研究区为典型的灌溉农业区，主要水源为祁连山融雪。玉米和春小麦是该地区种植的主要农作物。玉米一般在4月下旬播种，9月中旬收获，种植行距40 cm，株距30 cm，研究区作物密度约为6.6 万株/hm2。大满站的实验位置和研究地点见图1。

图1 大满站的实验位置和研究地点Fig.1 Experiment locations and study sites at Daman Station

1.2 测量与处理

大满站的观测数据是从黑河流域联合遥测实验研究（HiWATER）项目的野外观测系统中收集的，详细描述参考李等[17]。涡度（EC）相关系统的架设高度为4.5 m。涡度数据中包含潜热通量和感热通量，均为半小时尺度，数据处理包括：野值点剔除，延迟时间校正，三维旋转[18]，频率响应修正，超声虚温修正和密度校正等，频率响应校正等[19]。并删除由于降水、水汽凝结、系统故障导致的虚假数据。站点仪器测量的数据主要有降雨量、气温、风速、大气压强、土壤水分、土壤热通量、净辐射等数据。在生长季，叶面积指数(LAI)大约每10 d测量一次。

1.3 模型与方法

1.3.1 Shuttleworth-Wallace(SW)模型

SW[13]模型包括植物蒸腾的一维模型和土壤蒸发的一维模型。在模型中，地表和冠层阻力调节土壤和植物表面之间的质量和热量传递，而空气动力阻力调节两表面和大气边界层之间的阻力[13,14]。土壤蒸发和植物蒸腾用以下公式[13]计算：

式中：λET为来自植物（λT）和土壤（λE）的潜热通量之和，W/m2；Δ 为饱和水汽压曲线斜率，kPa/K；γ为湿度常数，kPa/K；Cp为干燥空气的比热容，1 013 J/(kg•K-1)；es和ea分别为饱和水汽压和实际水汽压，kPa；raa为植物冠层高度与参考高度之间的空气动力学阻力，s/m；rca为冠层内边界阻力，s/m；rsa为地面与冠层高度间的空气动力阻力，s/m；rcs为植物群体冠层阻力，s/m；rss为下垫面土壤表面阻力，s/m。

A和As(W/m2)分别为冠层和土壤表面以上的可用能量，定义为：

式中：Rn和Rns分别为冠层上方和地面接收的太阳净辐射，W/m2；G为土壤热通量，W/m2；Ka为消光系数；LAI为叶面积指数。阻力模型的计算参考魏等[1]。

1.3.2 贝叶斯推理和DREAM 算法

参数的后验概率分布由贝叶斯定理计算：

式中：π(θ|M)表示模型M下的θ的先验密度；p(D|θ,M)是模型M及参数θ的联合似然性。边际似然或贝叶斯模型证据是：

用于参数估计的似然函数p(D|θ,M)是根据观测误差的分布指定的。在时间t的每个观测值D（t）中的误差有e( )t表示为：

假设e(t)服从均值为零的高斯分布，似然函数可以表示为：

式中：n是观测次数；σ表示误差方差。

本文使用DREAM 算法来探索ET模型的参数空间，并优化SW 模型中的参数。DREAM 采样方案是对混洗复杂进化Metropolis (SCEM-UA)的全局优化算法的改编。Vrigt 等[15,16]更详细地描述了该算法。

1.4 评估模型性能的传统统计指标

用于评估模型性能的传统误差度量包括决定系数（R2）、一致性指数（IA）、相对误差度量（EF）、均方根误差（RMSE）、平均偏差误差（MBE）。

1.5 不同驱动因子对ET变化趋势的贡献

本文采用偏微分方程研究气候和植被因素对ET变化趋势的贡献。通过这种方法，可以分离出驱动生长季蒸散发变化的气候和植物因素的相对贡献。李等[9]使用的Penman-Monteith模型将ET表示为4个气候变量和一个植物变量的连续函数:

本研究把SW 模型中的ETs和ETc分别表示为两个不同的连续函数，并假设各个变量是相互独立的。用偏导数计算ETs和ETc的变化，并分析各驱动因子对ETs和ETc的相对贡献。并用作物潜热通量(λT)和土壤潜热通量(λE)(W/m2)来表示ETc和ETs。

λE可以表示为A、As、es-ea、raa、rsa、rss的连续函数：

同样地，λT可以表示为A、As、es-ea、raa、rca、rcs的连续函数：

λE和λT的相对变化可以用泰勒展开级数表示为一阶偏微分形式。

λE的相对变化：

类似的,λT的相对变化可以表达为：

λET是通过涡度观测系统测得的，并利用SW 模型和实测气象资料计算λE和λT驱动因子的偏导数。为方便表述，下文用LE表示潜热通量。

2 结果与分析

2.1 能量闭合分析

所有表面能量和质量交换模型都基于基本守恒定律；即能量守恒和质量守恒。能量守恒方程的主要组成部分，通常称为“能量平衡闭合”，可以描述为：

式中：H为显热湍流通量，LE为潜热湍流通量，Rn为净辐射，G为地表热通量。能量平衡比（EBR），定义为H和LE之和与Rn和G之差的比值。在大多数站点EBR为70%～90%，这称之为能量不闭合问题。

2016-2018年基于半小时白天湍流通量（H+LE）总和与可用能量（Rn-G）的线性回归见图2。由于晚间风速小，或者湍流不发生，以平流的方式传递通量，导致EC测得的通量严重不闭合。因此，本研究只选择白天的数据，并定义“Rn-G>0”为白天。结果表明，2016-2018年的EBR分别为0.75、0.89、0.84，决定系数（R2）也都在0.7 以上，说明拟合效果良好。

图2 2016-2018年能量平衡闭合Fig.2 Energy balance closure from 2016 to 2018

2.2 参数优化及模型性能

SW 模型中有7 个参数gmax、Q50、D50、Kq、Ka、b1和b2（表1）。每个参数的先验概率密度被指定为均匀分布。2016-2018年7-9月的数据用于估计每个参数的后验概率密度函数，并使用DREAM 算法生成了50 000N个样本来实现。在计算中，链数N等于相关模型中的参数数，对于SW 模型，N等于7。对于每个链，前10 000 个样本作为老化数据丢弃，其余40 000个样本用于设置后验密度函数。优化后的SW 模型参数值的最大似然估计见表1。

表1 模型参数优化后的最大似然估计Tab.1 Maximum likelihood estimation of model parameters after optimization

使用2016-2018年白天每半小时数据评估了优化后的SW模型性能。使用模型的校准参数，并运行SW 模型来估计每半小时的LE值。SW模型模拟与测量LE值的线性回归见图3。一般来说，SW 模型与测量的LE有较好的拟合效果(R2>0.77)。图3可以看出，虽然模型和观测值仍存在一定的差异，但SW 模型能够捕捉到蒸散发的动态变化。

图3 2016-2018年SW 模型和实测LE线性拟合Fig.3 The SW model was fitted with the observed LE from 2016 to 2018

半小时蒸散量实测值与模拟值的R2、RMSE、IA、EF、MBE见表2。基于对这些传统误差度量的分析，这三年的IA和EF比较接近，而MBE分别为-9.61、13.82、-19.11。从总体上来看经过参数优化后的SW模型能够准确模拟ET。

表2 模型性能的传统统计指标Tab.2 Traditional statistical metrics of model performance

2.3 ET及驱动因子的年平均变化

2016-2018年SW 模拟的ET、ETs、ETc、A、As、es-ea、raa、rsa、rss、rca、rcs白天每半小时的平均值见表3。2016-2018年ET的实测值分别为195.62、232.17、213.66 W/ m2，从表3可

表3 2016-2108年白天每半小时的均值Tab.3 Mean values for every half hour during daytime from 2016 to 2108

以看出2016年和2017年SW 模型模拟的ET均值与实测ET相差不大，2018年则模拟的偏大。2018年ETs为138.07 W/m2，明显比前两年偏大，主要是因为As模拟值偏大，导致分配给冠层的可用能量偏小，所以ETc的值就偏小。阻力方面，其中raa、rsa、rca三年变化不明显。2017年rss值为317.13 s/m，明显比2016年134.37 s/m 和2018年85.17 s/m 大；2018年rcs值为217.73 s/m，与其余两年相比偏大。这些差异可能是导致玉米农田ET年际变化的潜在因素。

2.4 驱动ETs和ETc变化的因素

由A、As、es-ea、raa、rsa、rss对ETs变化的贡献率和A、As、es-ea、raa、rca、rcs对ETc变化的贡献率见图4。它由方程式（14）和（15）得到，贡献率的总和应该等于1，然而由于各种不确定的因素导致结果并不严格等于1。从图4可以看出ETs与A、As、es-ea、raa呈正相关与rsa和rss呈负相关。与2016年相比，ETs分别增加了2.56 W/m2、72.9 W/m2（表3），这主要是因为As的增加。2017年和2018年As分别增加了55.09 W/m2和156.93 W/m2。ETs的增加主要归因于As、es-ea、A的变化（图4），并且rsa和rss的变化抵消了部分对ETs的正贡献。其中，对ETs变化最大的正相关驱动因子为As，最大负相关驱动因子为rss，分别贡献了124%、95.4%和-99.1%、-26.9%。对于影响ETc变化的驱动因子对2017 和2018年的贡献有较大区别。其中，2017年ETc与A、es-ea呈正相关，与As呈负相关；而在2018年，ETc与A、es-ea呈负相关，与As呈正相关。rcs对于ETc的变化都是正贡献，并且占了很大比例，分别为68.8%和71.9%，间接说明植物因素在调节ETc中起着不可忽视的作用。

图4 2016-2018年不同驱动因子对ETS和ETC变化的相对贡献Fig.4 The relative contributions of the different driving factors to changes in ETs and ETc,from 2016 to 2018

3 讨论

3.1 能量闭合分析

涡度相关（EC）法是基于对垂直速度波动（w′）和标量浓度波动（c′）乘积的直接测量，并假设平均垂直速度可以忽略不计，从而直接估计H和LE；Rn和G分别是通过净辐射仪和土壤热通量仪测量的。表面能量收支不闭合的原因有很多，Mahrt[20]讨论了几个原因，其中包括：①在非常接近表面的地方测量的各种通量分量之间的源区域缺乏重合，例如叶子的蒸发和来自炎热干燥土壤表面的显热；②测量时间为典型的30 min可能会错过由非常低的频率波动引起的协方差；③仪器位置的平均垂直速度不可能系统的为零，从而产生垂直平流通量；④与传感器分离、频率响应、对准问题以及来自塔或仪器安装结构的干扰有关的测量误差等。而Foken[21]表明以前关于测量误差或存储项是未闭合能量平衡的原因的假设站不住脚，因为即使是来自校准传感器的湍流通量、净辐射和地面热通量也无法闭合能量平衡。作为一个假设，他将能量平衡闭合问题假设为一个尺度问题，并且闭合只可能在景观尺度上发生。由于能量不闭合的原因至今没有一个明显的定论，并且在极端情况下，整个能量不闭合可能归因于潜热通量[22]。研究结果表明，H+LE的真实表面通量可能比EC 测量值高20%～30%。能量不闭合可能会对SW 模型的不确定性产生影响，模型使用的H值比实际值高，从而导致模拟的LE偏高。

3.2 不确定性和误差分析

关于DREAM 算法的效率，SW 模型的接受率远高于之前研究中使用的一些MCMC 算法获得的接受率[23]。DREAM 算法可以有效地处理涉及高维、多模态、非线性的问题。使用DREAM 算法优化的参数Ka，可以明确地将总可用能量划分为SW 模型中冠层和土壤吸收的能量。方程（5）的分析发现Ka的变化不仅可以解释消光效应，还可以纠正能量驱动数据的错误。这也意味着使用标定数据估算的Ka值实际上不仅仅是真正的消光系数，还包括SW 模型中的能量不平衡校正。从这个分析中，可以看到，Ka不仅涉及到冠层和土壤表面之间的能量分布，还涉及到能量的不平衡。因此，参数Ka对SW 模型的性能有很大的影响。此外，由于Ka影响冠层和土壤的可用能量分配，Ka值较低导致分配给土壤表面可用能量偏大，从而导致模拟的ETs偏大。

试验结果可能存在不确定因素。首先，对ETs和ETc估计的误差可能归因于在泰勒展开式中考虑的不同驱动因子独立性的假设。事实上，这些变量并不是完全独立的，它们确实在一定程度上相互影响，这导致试验分析存在一定的不确定性。其次，这些变量并不是固定的值，而是随着环境的变化而变化，同时也会受到测量误差的影响。再次，泰勒偏微分方程适用于两个比较相近的值，由于所考虑的变量的年变异性造成的差异导致了较大的误差。在本研究中，没有考虑人为因素和下垫面条件的影响，因此对每个变量的分析引入了一些不确定性因素。除了测量和计算误差外，模型的选择还存在一定的不确定性。以往的研究表明，冠层阻力和空气动力阻力模型的不同会导致不同的模拟结果[24]。此外，模型的结构也是不确定性的进一步来源。本研究只是用DREAM 算法优化SW 模型参数使模型达到良好的拟合效果，并没有土壤蒸发和植物散发的实测值，这也是导致结果不确定的重要因素。

3.3 归因分析

本文利用泰勒展开法求解偏微分方程，分析了气候和植被变量对SW 模型中ETs和ETc的相对贡献。结果表明，土壤表面可用能量(As)是影响ETs变化最主要的外因，这或许与分配到土壤表面的净辐射有关。这与ZHANG 等[25]揭示净辐射下降是影响青藏高原东部蒸发皿蒸发的主要因素相符。而LIU 等[5]认为风速是影响是青藏高原蒸发皿蒸发的主要因素，这可能与风速影响了空气动力阻力和饱和水气压差有关；并且李等[26]发现空气动力学项是造成东北气象站参考作物蒸散发量下降趋势的主要原因。由此可见，ETs年际变化主要受气候变化的影响，尤其是As的影响。对于ETc而言，rcs变化是导致ETc变化的主要植物变量。Wilson 等[27]研究了北美温带阔叶林ET的年际变化，结果表明不同年份ET的变化主要是由于冠层导度造成的，其次是净辐射的变化。LI等[28]表明冠层导度主导了能量分配过程，然后发现ET的年际变化主要归因于冠层导度的变化而不是气候因素的变化[9]。主要是因为冠层导度是调节ETc的重要因素。这些结果表明，植物驱动因素对ET有重要影响，在未来气候条件下模拟蒸散发应更加重视冠层导度的作用。