基于最优性原理的普适性碳水通量耦合估算方法研究

谭 深,王 焓,*

1 清华大学 地球系统科学系,北京 100084

2 清华大学 全球变化研究院,北京 100875

植物光合作用与地表蒸散是全球碳、水循环的关键环节,也是自然生态系统生产和服务的基础。精准计算光合及蒸散过程产生的碳水通量既是地球系统科学的重要研究内容,也是实现资源高效管理、利用的前提[1]。然而在地球系统动态预测[2]、基于遥感(RS)技术的资源监测研究中,碳、水通量计算方法表现出无法忽视的不确定性亟待解决[3]。部分研究对机理认识不够深入,无法准确、完整地表达地表过程对环境变化的响应是其重要原因[1,4—6]。当前研究常通过结合土地覆被分类数据产品,将属性相同或相近的区域划分为同一类别(或植被功能型)作为基本单元,根据经验或实测数据为每个类别赋予“特征参数”。样本代表性不足和机理缺失导致站点尺度构建的模型在时空外推过程中产生较强的不确定性[7]。另外,部分研究剥离耦合的光合和蒸腾过程,忽略二者联系与相互约束,可能导致结果误差[8]。因此,深入理解植物光合过程和地表蒸散过程,厘清二者联系,同时避免采用特征参数将地表属性离散化,从而提升碳水通量计算方法的稳定性和应用的可靠性,具有重要研究意义[9]。

量化植物行为对环境的响应规律、简化庞杂的参数化方案,是避免特征参数对模型的限制、提升模型稳定性和普适性的有效途径[10]。生态水文最优性原理(以下简称最优性原理)通过对大量观测数据的分析和归纳后认为,植被能够适应环境的变化,调整自身理化性质达到资源最优利用,实现净碳收益最大化[4,11—13]。Wang等[14]整合最优性原理和光合作用研究成果,构建了不包含特征参数的C3植物普适性光合模型。该模型在全球范围的验证研究中取得了良好的表现[15],作为新一代植被模型的核心,在全球总初级生产力监测[15]、作物产量模拟和预测[16]、全球变化等研究中得到了广泛的应用[17],能够与彭曼 (PM,Penman-Monteith)公式结合估算蒸腾[18]。

我国幅员辽阔,下垫面复杂多样,碳、水通量的估算往往需要借助参数化方法和特征参数实现。受大尺度土地覆被数据产品准确性和模型参数化方案代表性不足的影响,不同研究估算的全国尺度碳水通量结果差异显著。本研究分别从站点尺度和全国尺度,论证基于最优性原理构建的普适性光合模型P model耦合估算GPP与ET的可行性。研究目标的实现标志着无需率定模型参数即可准确估算光合和蒸散速率。从而为农业管理、区域资源评估、陆面模式发展等研究提供方法借鉴和数据支持。

1 方法与数据

1.1 普适性生产力模型P model

最优性原理认为植物能够通过对环境的适应,调整自身理化性质,实现对能量、水分、养分资源的优化利用[4,13]。光合作用生理生化模型中,瞬时光合速率由光反应和暗反应过程共同限制[19],二者均受由气孔开闭程度调节的胞间二氧化碳浓度(ci,μmol/mol)影响。针对气孔行为的最优性理论研究从水分的角度,认为植物在生长过程中,能够权衡气孔开放所产生的光合收益和蒸腾损耗,最大化碳同化效率,实现环境适应[12,20]。Wright等[21]和Prentice等[22]在其基础上,将植物维持光合能力产生的损耗纳入权衡网络,提出了最低消耗假说：即植物的气孔调节倾向于最小化用于维持羧化和蒸腾速率的碳消耗,并与经典光合作用模型整合,得出ci的理论模型：

(1)

(2)

式中,ca环境二氧化碳浓度(μmol/mol)；Γ*二氧化碳补偿点(Pa),表达为对温度Ta的函数：Γ*=4.08×exp [(27055.67/8.3145)(1/298.15-1/Ta)]；η*为水的粘滞阻力,表达为对温度的函数：η*=exp {580[1/(Ta-138)]-1/160}；VPD为饱和水汽压差 (Vapor Pressure Deficit, Pa)；K为Rubisco酶的米氏系数(Pa),表达为：K=Kc(1+209460/KO),其中,Kc和KO分别为Ta温度下羧化和氧合反应的米氏系数,分别通过Kc=40.41×exp [(64805.5/8.3145)(1/298.15-1/Ta)]和KO=27480×exp [(36164/8.3145)(1/298.15-1/Ta)]定量计算。Wang等[14]在全球尺度,利用叶片稳定碳同位素观测数据,律定公式(1)中β的取值为146,并对该理论模型进行了验证。

基于最优性原理,关于植物光合作用生化过程的协同限制假说认为,在长期的适应下,植物能够在周到月的尺度调节叶片内Rubisco酶活性,使得羧化反应速率和光反应速率趋于一致,以实现对资源的充分利用[23]。Wang等[14]将协同限制假说与ci的模型,以及光合生理生化模型进一步整合,构建了C3植物普适性光合模型P model,并进行了全球尺度的验证。Stocker等[15]进一步考虑土壤水分对光合过程的胁迫[24],给出最新版本P model表达式和模型开源代码(https://github.com/stineb/rpmodel)：

(3)

其中,GPP为总初级生产力,是全部叶片光合作用在冠层尺度上的体现。Iobs为冠层截获的光合有效辐射 (mol photon/(m2·s)),表达为光合有效辐射(PAR, Photosynthetically Active Radiation)与冠层截获能力fAPAR(fraction of Absorbed Photosynthetically Active Radiation)的乘积,即：Iobs=fAPAR×PAR,c*为固定值0.41,β(θ) 为由土壤含水量SWC (Soil Water Content, m3/m3)驱动的水分胁迫项。φo为内禀光量子效率(mol CO2/mol),可表达为温度Ta(℃)的函数[25]:

(4)

式中的m项体现了CO2对光合作用的限制,其表达式为：

m=(ci-Γ*)/(ci-2Γ*)

(5)

关于P model的详细推导过程可参考Wang等[14]和Stocker等[15]。

基于P model的C3植物的光合速率计算方案可拓展至C4植物。C4植物具有比C3植物更强的CO2亲和力[26],所以我们假定C4植物的光合过程不受到环境CO2浓度的限制,即令m=1[18,27]。另外,C4植物的内禀光量子效率响应函数采用Kubien等[28]提出的计算方案

(6)

1.2 蒸散估算

水蒸气分子经由叶片气孔向外扩散的速率通过扩散方程(the Fick′s law)表达。由于GPP代表冠层所有叶片的光合贡献,冠层尺度的气孔导度(Gs)可表达为：

(7)

其中,1.6表达CO2分子与水分子扩散速率的差异；χ为气孔内外二氧化碳浓度分压ci与ca的比值,即ci/ca。对于C3植物,其取值可通过P model定量预测(公式1)；而对于C4植物,该比值变化相对保守。本研究采取Collatz等[29]的方案,将C4植物χ取为固定值0.45。

地表蒸散包括冠层蒸腾(Ec)和非生物蒸发(Es,包括土壤和冠层截留蒸发)两个主要部分,对应的潜热通量表达为：

λE=λEc+λEs

(8)

式中,λ为水的汽化潜热(MJ/kg)。蒸腾对应的潜热通量λEc可通过PM公式计算[30]：

(9)

其中,Δ为饱和水汽压随空气温度变化的斜率(kPa/K),γ为干湿表常数(kPa/℃),Qn,c为冠层吸收的辐射能量(W/m),表达为地表可利用能量(净辐射Rn与土壤热通量G之差)与冠层截获能力的乘积[30—31]。ρ为空气密度(kg/m3),cp为干空气定压比热(J kg-1K-1),ga为空气动力学阻抗。由于摩擦风速观测的缺失,本研究利用风速的函数计算ga[32]：

1/ga=208/u

(10)

研究表明,在没有外界干扰的情况下,蒸腾(T)与蒸散的比值(T/ET)受到能量和供水条件的协同限制,会在一定区间内((70±9)%)变化[33]。Tan等[18]基于这一原理建立了T/ET响应环境因子(Rn、Ta、SWC、fAPAR)的经验关系,并认为该关系在时间序列上保持不变。本研究沿用这一方案通过蒸腾计算蒸散。考虑到植物对环境变化的适应周期,本文以周为步长计算站点尺度GPP与ET；而对于全国尺度的通量模拟实验,结果时间分辨率为8 d,空间分辨率为500m,与MODIS产品保持一致。

1.3 数据

本研究分别从站点尺度和全国尺度论证碳水通量模拟结果的合理性。站点尺度的模型验证基于第二批发布的中国典型生态系统碳水通量数据集ChinaFLUX实现[34—36],共使用来自7个站点(3个森林站点、3个草地站点和1个农田站点)54站年观测资料(表1)。其中,禹城站下垫面为冬小麦和夏玉米轮作农田,玉米生长季内观测资料用于C4模型的论证。原始10Hz通量记录经过质量控制与缺失数据插补,汇总至30分钟数据集；再基于夜间的观测数据,建立呼吸速率与温度的关系,将观测碳通量拆分成GPP与呼吸[34]。本研究选用MODIS产品 (MYD15A2H)作为遥感fAPAR输入：首先利用该产品质量控制文件排除低质量数据,再经过Savitzky-Golay滤波器去除高频噪声后[37],构建完整时间序列的fAPAR数据集作为输入。

表1 本文涉及ChinaFLUX站点信息

全国尺度碳、水通量估算基于谷歌地球引擎云计算平台实现[38]。针对数据池中的MYD15A2H产品进行了异常值剔除和质量控制,同样利用Savitzky-Golay滤波器平滑异常值并构建时间序列完整的数据集,结合汇总、筛选后的GLDAS(V2.1)气象驱动产品、ERA5土壤水分再分析产品共同输入P model。为了体现C3与C4植物的通量贡献差异,研究采用全球尺度的C4植物面积比例产品[39],以每个格点内两类植物通量贡献的面积加权平均作为总地表通量。

GPP估算结果通过与2007—2015年GOME-2(Global Ozone Monitoring Experiment-2, V27)太阳诱导叶绿素荧光(SIF, Sun-Induced Chlorophyll Fluorescence)产品对比,评估结果空间特征的合理性[40—41]。逐月SIF产品经过异常值剔除,计算多年平均作为对照。基于经验关系估算的T/ET结果与Niu等[42]基于参数化方法估算的我国T/ET产品(2003—2015年)对比[43]；2003—2018年ET估算结果则与PML-V2[31]、GLEAM-V3.3[44—45]和MODIS(MOD16A2)[46]3种基于参数化方法发展的全球ET产品平均值对比(PML截止于2017年)。

2 结果与分析

2.1 GPP估算结果

总体而言,站点尺度GPP的估算结果与实测数据表现出良好的一致性(图1)。综合全部样本,P model对GPP的预测精度相关系数R2= 0.61,均方根误差 (RMSE, Root Mean Square Error)= 2.1gC/d,拟合斜率为0.96。模型在具有显著物候特征的站点(如长白山站和海北站)、以及植被覆盖度更高的站点(千烟洲站)具有更好的预测效果。模型估算结果与实测数据鼎湖山站表现出很好的一致性,但是存在高估(拟合斜率约为1.48)。该站点地处广东,云雨天气导致遥感fAPAR相较于其他站点具有更少的有效观测和和更低的质量,在3-6月甚至存在存在长时间(超过6周)缺少有效数据的情况,导致时间重建后的fAPAR存在高估,影响估算效果。模型在植被稀疏的当雄站表现出较低的相关性。由于当雄站地处高寒草甸,遥感观测信号包含大量非植被信息,影响模型估算结果。这也为模型的应用提供了重要经验,即选择质量更好、精度更高的fAPAR产品作为输入能够提升模型表现,后续卫星精细谱段设计工作也应当提升植被指数相关波段对光合过程的敏感程度[47]。全国范围内2007—2015年平均GPP与同时间段的平均SIF信号空间分布一致(图2),但由于值域分布和SIF产品空间分辨率较低的原因,模拟GPP在我国南方地区表现出更加丰富的空间细节和更大的值域变化范围。证明P model 能够在无需特征参数标定的情况下,提供稳定的GPP估算结果。GPP估算模型的稳定表现也是后续准确估算ET重要保证。

图1 GPP预测结果与通量观测对比

图2 基于P model估算我国2007—2015年平均总初级生产力(gC/a)和 2007—2015年平均GOME-2观测SIF信号/(mW/sr)

2.2 ET估算结果

本研究沿用Tan等[18]提出的全球尺度T/ET估算方法,与ChinaFLUX的站点观测表现出较好的一致性(R2= 0.90,图3)。该经验关系主要由fAPAR驱动,表明在植被茂盛、环境湿润的地区,冠层对辐射的截获能力更强,蒸腾占ET的比例更高。由于该方法缺少对冠层截留蒸发过程的独立刻画,所得结果在植被覆盖程度较高、叶面积较大的站点略有低估(图3)。另外,基于该经验关系法的T/ET估算结果与Niu等[42]在中国发布的基于参数化方案构建的T/ET产品空间趋势接近,但变化范围较小(图4)在植被覆盖程度较低的西藏阿里、日喀则等地具有更高的蒸腾贡献比例。总体而言,本研究采用的经验方法所得结果更符合已有研究提出的(70±9)%变化区间[33]

图3 T/ET预测结果与站点观测对比

图4 经验方法与参数化方法T/ET计算结果

与站点ET观测的对比结果表明,基于P model的普适性蒸散模型表现稳定(图5)：相关系数R2= 0.66,RMSE = 0.85mm/d,拟合斜率为1.04。模型在内蒙古和禹城站表现稍差,主要由于内蒙站地处草地,通量贡献源区范围内可能存在部分C4草本植物,影响通量观测的代表性。另外,本研究将C4植物的χ简化为常数(0.45±0.16)。根据对C3植物的环境响应规律,夏季χ值通常较冬季更高,由于C4植物缺少对环境的动态响应,可能高估禹城站玉米种植期间的导度,进而增加蒸散估算的不确定性。

图5 ET估算结果与通量观测对比

总体而言,全国尺度的ET估算结果与其他方法结果具有较好的空间一致性,但值域偏低(图6)。部分南方地区,如四川盆地和贵州地区,P model估算的结果要显著低于参数化方法的平均结果。本研究的方法将非生物蒸发(包括土壤蒸发和冠层截留蒸发)作为整体通过T/ET刻画。南方地区由于冠层结构复杂、雨水充沛,可能导致低估截留蒸发。另有研究表明,部分参数化产品在我国较湿润的流域存在高估现象[48]。因此,本研究方法所得结果仍然需要更加系统的水平衡验证和时间序列变化趋势验证,以论证其稳定性。

图6 基于P model和参数化模型估算我国2003—2018年平均ET(mm/a)空间分布

3 讨论

本研究基于生态水文最优性原理,根据植物对环境的适应规律,拓展已有普适性C3植物光合模型P model至C4植物；并以气孔导度作为结点,耦合估算我国碳、水通量。与地面观测、遥感SIF产品和同类产品的对比证明了本文采用的普适性模型的应用潜力。传统参数化方法以地物类型作为基本研究单元的处理思路能够在充分率定参数和输入数据完整的研究中取得良好表现,但可能导致模型在未来情景或其他区域拓展的过程中产生问题,其原因包括以下几个方面。首先,模型选择、参数标定的过程受到研究区域特点、观测样本量、模型结构、拟合方法等因素影响。甚至相同站点的不同观测方法可能导致参数拟合结果的差异[49—50]。严重增加了模型标定的工作量,为模型的应用造成不便,也限制了模型的外推效果。

图7 标准状态(鼎湖山站平均状态)下不同导度模型模拟结果对比

另外,土地覆被产品的精确性和准确性也会影响通量计算精度。常见遥感分类产品空间分辨率介于10—103m之间,低分辨率分类产品在地物复杂的地区会产生混合像元,影响分类效果。对于参数化模型而言,不同地类的特征参数具有显著差异,甚至于归属统一地类的不同物种间的参数都会存在很大差异[49]。如图7示例为假设的标准状态下(接近鼎湖山站的平均状态,温度=20 ℃,VPD=10 hPa,fAPAR=0.8,ca=380 μmol/mol,SWC=0.2 m3/m3,PPFD=275 μmol photon/(m2s))常见气孔导度模型(Jarvis模型[51]、Ball-Berry模型[52]、Leuning模型[53]和Medlyn模型[20],模型参数参考已有研究[20,46,50])在使用不同特征参数情况下模拟结果的差异。能够看出,在土地覆被产品提供错误分类信息的情况下,模拟结果将会产生显著偏差,普适性模型则不会由于地类混淆产生问题。然而,本研究利用P model在进行大尺度通量模拟的过程中采用了较低分辨率的C4植物面积比例数据,通量模型复杂的非线性结构会对估算结果造成一定影响。

由于建模过程对客观规律的认识不够深入,部分模型假定的环境响应规律与真实状况不完全一致,但参数化过程和不充分的观测样本可能掩盖这一问题。图8展示了常用的导度模型在面对升温和CO2增浓过程中的敏感性对比,除可以通过参数率定消除的系统偏差外,不同模型对环境变化的响应趋势有所差异。如Jarvis类导度模型常假设实际导度是最大导度与环境胁迫项的乘积,并假设温度增加会降低环境对导度的胁迫[46],结果与水碳耦合类模型的估算结果产生相悖的趋势,该模型也忽视了环境CO2浓度增加对导度的影响(图8)。

图8 不同导度模型响应环境因子变化示意

为了满足应用的要求、提升模型的精度,通量估算过程中的参数化过程逐渐细化。增加模型中的参数数量和模型的复杂程度不仅对其运行的稳定性和普适性产生严重阻碍,更可能导致其稳定性下降。因此,避免采用过为复杂的模型结构、减少参数数量、从而避免可能产生的过拟合现象是未来模型发展需要考虑的重点[5]。

本文采用的普适性碳水通量估算方案仍然需要发展和改进。前文提到,当前C4植物光合和蒸散计算方案基于C3植物的计算方案,考虑其具有更好的CO2亲和力而简化得到,并且χ的估算方法也缺少对环境的动态响应。未来计划通过更多的地面实验,探索C4植物气孔的环境适应规律,改进C4植物蒸腾、蒸散计算方案。

本研究的GPP、ET估算结果在四川盆地略有低估(图2、图6)。造成这一现象的主要原因是云层对光学遥感观测的污染。四川盆地地处亚热带,特殊的地形条件导致其内部湿热、多云,对遥感观测和地表参数反演造成了严重影响。本研究采用的光能利用率模型和蒸散模型对遥感fAPAR产品具有较强的敏感性,在该地区易受到影响导致结果产生偏差。已有研究证明,融合多种数据源的遥感通量观测产品能够在局部地区取得较好的观测精度[54]。计划在后续的研究中同化多源遥感数据,提升监测方案的稳定性。这一不足也表明,非参数化模型对输入遥感参数精度具有较高的要求。当前使用的MODIS产品在反演过程中缺少对散射辐射的考虑,在阴天可能导致误差[55—56]。

另外,蒸散计算过程中的重要系数T/ET的相关研究近年来成为了研究的热点[57],但不同方法之间仍然存在巨大争议[58]。本文采用了简单的经验拟合方法,以期在保持模型普适性的前提下便于后续研究改进。根据本文结果,该方案能够与参数化方法取得较好的一致性,但无法反应非生物蒸发在时间序列上的变化趋势。已有研究表明,T/ET取值会在(70±9)%产生变化[33],因此该经验方案在地表覆被或气候类型变化剧烈的地区可能产生一定的误差。后续计划通过与其他蒸散分离方法[59—60]、非参数化蒸散产品对比[61],构建更加严谨、准确的蒸发计算方案,提升蒸散估算精度。

4 结论

本研究基于生态水文最优性原理,将已有的普适性C3植物光合模型P model拓展至C4植物,并以气孔导度作为结点,耦合估算我国碳、水通量。站点尺度通量验证基于ChinaFLUX数据集7个站点共54站年的地面观测实现。结果显示,P model对GPP的估算结果相关系数R2= 0.61,RMSE = 2.1gC/d,拟合斜率为0.96；ET的估算结果相关系数R2= 0.66,RMSE = 0.85mm/d,拟合斜率为 1.04。全国尺度通量计算基于Google Earth Engine云平台实现。GPP模拟结果与遥感SIF观测具有较好的一致性；ET模拟结果与其他产品相比空间趋势较为一致。证明基于最优性原理构建的普适性碳水通量耦合估算模型稳定可靠,能够在无法获取高精度土地覆被信息以及模型训练样本不足的情况下取得稳定的计算结果。