光合作用利用光能进行水的光解, 释放O₂, 同化大气中的CO₂, 形成有机物, 是生态系统生产力形成的基础, 也是能量流动和物质循环的起点。准确模拟光合作用具有重要的科学意义和应用价值。在20世纪80年代初期, 随着Farquhar光合作用生化过程机理模型的提出, 这项工作有了长足的进展(Farquhar et al., 1980)。现今, Farquhar光合模型有了非常广泛的应用, 已成为叶片尺度光合作用模拟的主流机理模型。驱动模型运行的一些气象因子, 比如辐射、温度和湿度等, 在具体应用中容易获得。模型应用的一个难点是生理参数的确定, 特别是最大羧化速率(V_cmax)和最大电子传输速率(J_max)的确定。众所周知, 生理参数V_cmax、J_max在不同物种间差异很大, 同一个物种在不同季节间的变化也很显著(Xu & Baldocchi, 2003), 想要获得这些参数并将其应用到区域或者全球尺度的研究中, 难度很大。模型应用的另一难点在于尺度的扩展, Farquhar光合模型是叶片尺度的模型, 若想将该模型运用到冠层或者更大尺度的研究中, 需要进行合理的尺度转换(Jarvis, 1995), 而这种尺度转换是生态学研究中尚未解决的问题(Amthor, 1994)。以热带森林为例, 测定数以百计的树种的V_cmax和J_max, 即使不考虑人力物力上的巨大投入, 如何将测定结果组织起来, 以便可以代表整个冠层, 目前还未找到合适的方案(Sprintsin et al., 2012; Cubasch et al., 2013)。因此, 亟需寻找到一种合理的尺度扩展方案, 作为独立的、可靠的验证标准, 冠层尺度的V_cmax和J_max直接实测值是必不可少的。

涡度相关技术(eddy covariance technique)作为目前直接测定地表-大气间CO₂和水热通量的标准方法, 为原位、无破坏监测生态系统的光合作用(实质上植被大气之间的CO₂交换)提供了新的选择(de Pury & Farquhar, 1997; Wang & Leuning, 1998; Dai et al., 2004; Groenendijk et al., 2011; Kosugi et al., 2013)。在过去的30余年里, 该方法有了快速的发展, 得到了广泛的应用, 目前已成为国际通量观测网络(FLUXNET)的主要技术手段。据不完全统计, 目前全球安装涡度相关监测系统的研究站已经超过5 000个, 仅中国就超过了300个。考虑到涡度相关法可以提供冠层光合作用的实测值, 而叶片尺度的光合模型也很明确, 笔者认为可以在叶片光合模型的基础上提出一种算法, 直接反演出V_cmax和J_max这两个重要的生理参数在冠层上的数值。如果这一想法能够实现, 我们就有可能从生理学角度切入, 对大量冠层尺度的监测数据进行更深入的解读和挖掘。作为一个初步的尝试, 为了简化模型的数据和计算, 本文通过将冠层抽象为一片大叶, 提出了一种相对简洁却行之有效的方法来反演冠层的生理参数。下文中, 笔者将对该方法的技术细节进行分解, 希望该方法能引起读者的兴趣, 在读者的批评中获得进一步提高和完善。

1 方法

1.1 贝叶斯计算法

统计分析经常面对的一个挑战是用观测数据估计多元非线性模型的参数。

y=f(x,θ)+ε

在上式中, x为自变量, θ为待定参数向量, ε是模拟对象y变量的观测误差。基于频率统计学(frequentist statistics)的参数估计方法, 一般是在假定ε服从某一统计分布的条件下, 通过优化求出固定的θ值, 该值的可信度, 或在某一可信度下包含θ真值的范围。而基于贝叶斯统计学(Bayesian statistics)的方法, 则是假定θ服从某一先验分布π(θ) (prior distribution), 依据贝叶斯概率法则, 通过优化估计θ的后验分布(posterior distribution):

$$ p(\theta | y)=\frac{p(y| \theta)\pi (\theta)}{p(y)}$$

$$p(y)=\int p(y|\theta)\pi (\theta)d\theta$$

其中, p(y|θ)为似然函数(likelihood function), p(y)为边缘似然函数(marginal likelihood function)。当似然函数和边缘似然函数易得, 或可以通过某种近似法, 诸如马尔可夫链蒙特卡洛法(Markov chain Monte Carlo, Andrieu et al., 2010)、重点取样法(importance sampling, MacEachern et al., 1999)以及序贯蒙特卡洛法(sequential Monte Carlo, Gao & Zhang, 2012)等算出时, 这些方法便可用来求解待定参数的后验分布。然而对复杂的非线性多参数模型而言, 似然函数往往不可知, 或者计算边缘似然函数有困难。在此背景下, 近似贝叶斯计算法(Approximate Bayesian Computation)应运而生(Toni et al., 2009; Beaumont, 2010; Marin et al., 2012), 并被广泛应用(Csilléry et al., 2010; Vrugt & Sadegh, 2013; Hartig et al., 2014)。在此, 我们介绍如何应用Lenormand等(2013) 提出的自调近似贝叶斯计算法(APMC)来估计植物光合模型的几个重要参数。

1.2 光合作用模型

基于将冠层视为一片大叶的思维抽象, 叶片尺度的光合模型便得以直接应用于冠层之上, 简称大叶模型。本文采用了Farquhar模型来模拟冠层尺度的光合作用, 同时采用Ball-Berry模型来模拟气孔导度关系(Ball et al., 1987)。具体使用到的关系式见表1; 各缩写和符号代表的含义和参数初值见表2。模型细节请参考Medlyn等(2002a), von Caemmerer等(2009)和Damour等(2010)的文献。

1.3 APMC算法的具体运行

首先我们假设大叶模型的叶面CO₂分压(C_s)、 胞内氧分压(O_i)、 叶面相对湿度(h_s)可以用环境空气的实测值近似代替(附录I)。之所以做出如此假设, 是因为: 一方面, 我们无法知道哪些叶片的C_s、 O_i、 h_s值能代表整个冠层; 另一方面, 从大多数过程模型使用的边界层导度模块来看, 这些参数值在叶片表层与大气之间相差不会太大。但是, 胞间CO₂浓度(C_i)是个例外, 它是整个模型运行的一个关键, 不能直接假定其与环境CO₂浓度相等。用来估计C_i的方法主要有两种, 一种方法是通过迭代得到, 另一种方法是用Baldocchi (1994)提出的3次方程解。虽然Baldocchi (1994)、Yin和Struik (2009)等都认为3次方程解适用于正向模拟的绝大多数情况, 但我们的实践表明它不适用于逆向参数估计, 因为逆向模拟时,参数值变化范围可能很大, 因此本文最终采用了迭代法。我们在实践中还发现, 大多数情况下只需要低于10次的迭代就可得到稳态的C_i解。

表1中的方程(10)是控制光合速率的关键方程。它实际上是在Arrhenius方程(9)基础上引入了

$$f(S,H)=\frac{1+e^{(298S-H)/(298R)}}{ 1+e^{(TK\times S-H)/(TK\times R)}}$$

进行温度修正, 使得V_cmax、J_max在高温下减小。正如图1所示, 修正曲线对某些范围内的参数S和H的变化很敏感。作为实例, 我们仅把von Caemmerer等(2009)温度修正模型中的S值增加约10%, 就使修正曲线系数的变化区间由(0.0, 1.1)扩大到(0.0, 420), 并在0 ℃附近开始急剧下降。这种敏感性对于正向模拟不构成问题, 因为参数调试时可以根据经验, 不断缩小范围逐步进行调整。而在逆向模拟参数反演时, 参数调试依据的是规则而非经验, 这种敏感性会导致非优化解。针对这一问题, Medlyn等(2002b)的处理方法是固定H的值, 通过详细的分析和探索, 我们建立了一个温度修正方程, 可以避免参数反演时数值不稳定性, 同时得到优化的反演参数。这个新的温度修正曲线表达如下:

$$f{C_m,T_m}=\frac{1+e^{-C_mT_m}}{1+e^{C_m\times (TC-T_m)}}$$

事实上, Seller等(1996)也用了类似的方程来修正Rubisco的活化能。新的温度修正曲线与June等(2004)提出的简易温度修正的不同之处是, 新修正曲线可以很好地拟合f(S, H), 且修正系数的变化范围小, 避免了一些数值计算问题(图1)。并且, 我们可以比较直观地选择C_m和T_m的可变范围。实际上, 对任意正的C_m、T_m和TC, f(C_m, T_m)的变化只限于(0, 1)区间。

显示原图| 下载原图ZIP| 生成PPT

图1   温度修正模型比较。实线: von Caemmerer等(2009)的曲线(S = 710, H = 220 000), 点线: 增加S约10%的曲线(S = 790, H = 220 000), 虚线: 本文用的曲线(C_m = 0.3, T_m = 37)。

Fig. 1   Comparison of the temperature correction models. Solid line is the curve of von Caemmerer et al. (2009) when S = 710 and H = 220 000; dotted line is the curve when S is increased by 10% (S=790, H=220 000); dashed line is the curve calculated according to our new response curve (C_m = 0.3, T_m = 37).

APMC算法的具体实现参见附录I。必要的输入数据包括NEE、PPFD、TK、RH、WS、P、LeafW、C_a。选定了模型、待定参数、输入变量, 就确定了APMC循环计算中y = f(x,θ)的计算。

1.4 验证数据的获取

验证用的数据获取自日本的苫小牧观测站(42.73° N, 141.52° E), 由AsiaFlux提供。其中净生态系统CO₂交换(NEE)白天的数值可以视为生态系统昼间的净光合速率。NEE是涡度碳通量和碳储量通量之和。涡度碳通量是三维超声风速仪所测得的垂直风速和红外气体分析仪测定的高频CO₂浓度值, 是在时间尺度上求协方差得到的。而储量通量则是通过CO₂浓度的垂直廓线在时间的变化计算得到。有关的观测方法、周边环境和数据处理细节等, 在Hirano等(2003)、Wang等(2004)和Saigusa等(2008)的文中已有说明, 此处不再赘述。研究选择了2003年6月至9月期间的数据。因为此这期间, 归一化植被指数(NDVI)显示树叶的总量变化比较小, 使我们可以不考虑叶面积指数(LAI)或树叶的年龄因子的影响。把有缺损和有效光量子密度(PPFD)小于100的数据点消除后, 一共得到1β948个完整的、 30 min时间间隔的有效数据点。

2 结果

我们用了一万个蒙特卡洛粒子做模拟实验, 并设N_a = 100, 即每次循环计算都接受10β000粒子中最好的100个。从验证方法的角度出发, 我们不打算讨论这100粒子的参数值的统计分布, 而只比较由最佳粒子的参数值产生的NEE与观测值之间的差距。

拟合得到的方程参数如表3所示。用APMC方法估计的参数值全部落在参数值的上下限范围内, 这表明从植物生理学含义的角度看, 这个模拟结果是可以接受的。需要注意的是, 因为我们用了不同的温度修正方程, 得到的表1中方程(11)的V_opt和J_opt的意义与表1中方程(10)的V₂₅和J₂₅的意义是不一样的。仿效Medlyn等(2002b), 我们把方程(11)叫做速率峰值方程。

通过图2, 可以初步看出模拟的效果。图中数据点一共1β948个, 考虑到涡度相关法测定时本身带有一定的随机误差, 决定系数达到0.75体现了很好的模拟效果。从斜率来看, 1.04与1非常接近。以自由度为1β946的t检验衡量, 达到极显著。NEE的模拟与观测平均值分别为15.66和15.32 μmol·m^-2·s^-1, 标准方差值分别为9.08和11.12 μmol·m^-2·s^-1, 二者之差的平均值和标准方差分别为0.34和5.63 μmol·m^-2·s^-1。这表明最佳粒子的参数给出了无偏差的NEE模拟。

显示原图| 下载原图ZIP| 生成PPT

图2   净生态系统CO₂交换(NEE)模拟结果(实线)与观测值(符号)的比较。图中实线为1:1线。

Fig. 2   Relationship between modeled and observed net ecosystem exchange (NEE). The solid line indicates that model outputs equal observations.

3 讨论和结论

虽然APMC方法在很多学科领域已有应用, 但是在植物生态学上的应用比较少见。本文将重点放在了介绍用APMC反演光合参数的技术细节, 而非具体的植物生态学问题。笔者希望通过应用实例, 打开APMC方法在植物生理生态学研究中应用的大门, 使该方法得到更广泛的应用。

利用APMC方法进行冠层光合参数反演的研究还未见报道。过去光合参数的获得(这里指V_cmax和J_max), 大多是通过测定光合速率对胞间CO₂浓度上升的响应曲线, 再结合Farquhar模型计算得到(详见Gu等(2010)及其参考文献)。正因如此, 笔者无法通过对比本研究与其他方法所得的冠层参数值来说明本方法的可靠性, 但是, 从模型在验证站点运行的情况来看, 还是非常理想的。本文研究结果表明, 实测值和根据反演参数模拟出来的数值具有很好的相关性, 并且斜率接近1。考虑到涡度相关法本身就具有较大的随机误差,该方法在具体应用上具有可行性。实际上, 如果设定的参数初始值范围适当, APMC方法几乎总是可以得出与观测NEE有良好线性关系的结果, 做出非常好的模拟。我们同时还将该方法应用在了叶片尺度以及其他通量站点, 也得到了很好的拟合效果(未发表数据)。

但是, 使用APMC方法进行冠层光合参数反演, 并非完美。比如, 如果对参数初始值的随机取样范围不加限制, 得到的参数值可能超出常识范围。这不是参数估计方法的问题, 而是因为模型方程的过度参数化(Medlyn et al., 2002b)。实际上, 因为植物的光合速率取决于A_c和A_j中小的一个, 直接影响最佳光合温度的参数就有10个。为了使曲线在拐点处有很好的表现, 所选择的驱动变量,要有一定的梯度, 当梯度不足(比如高温高湿的热带雨林)时, 参数初始值的选择就显得尤为重要。

本文将一种新的方法应用到已经成熟的光合作用模型中, 完成参数反演的目的。从植物生态学的角度看, 这个参数反演工作的重要性可以概括为以下3点:

(1) 该方法有望成为叶片生理参数反演的一种新选择。长期以来, V_cmax和J_max都是通过测定CO₂浓度响应曲线, 结合Farquhar模型计算出来的。进行CO₂浓度响应曲线的测定, 操作比较复杂, 并且基于这些数据来计算参数有较大的不确定性(Gu et al., 2010)。而我们的方法可以通过日常条件下的一系列监测值来实现反推生理参数的目的, 可以作为叶片生理参数反演的一种新选择。

(2) 使直接获得冠层整体生理参数成为可能。对森林植被而言, 想要对整个冠层进行控制, 进行CO₂响应分析, 在目前的技术条件下几乎没有可行性。然而, 我们这个方法提供了这种可能性, 利用涡度相关技术所测定的生态系统净CO₂交换量即可反推出整个冠层的光合生理参数。这是一个全新的视角, 对深入挖掘和利用涡度通量数据, 以及获得冠层生理参数都有重要价值。

(3) 有益于叶片到冠层尺度扩展问题的解决。如何将叶片尺度的机理模型应用到由众多树种组成的森林冠层上, 一直都没有很好的解决办法。一个比较传统的方法是通过选择典型树种向阳的几片代表性的叶片, 用光合叶室测定CO₂浓度响应曲线来推导出这几个种的光合生理参数值, 然后将平均值代入叶片尺度模型, 进行计算模拟。不论是多层模型还是大叶模型, 基本上都采用这种思路。对于农田生态系统而言, 物种单一, 冠层低矮, 在实际的应用中, 上述方法尚有一定的可行性。对于森林冠层, 特别是物种多样的热带森林, 该方法的应用则存在诸多困难。一方面, 要进入高达80 m的林冠层进行CO₂浓度响应曲线的测量是非常困难的; 另一方面, 热带森林物种多样, 很难确定到底选哪些种、取多少种能代表这个生态系统。本研究所提出的方法, 为这个问题的解决提供了一个新的选择。APMC模型结合方法机理明确, 模拟效果好。这个方法的提出和应用对光合模型参数的获取或许有积极的意义, 为解决长期困扰我们的尺度扩展问题提供了新的视角。

致谢 感谢日本的苫小牧观测站提供观测数据。

The authors have declared that no competing interests exist.

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子