1、利用 CO2 浓度观测提高地表碳通量模型估算精度的方法北京师范大学 郑珩、黄文贤、黄勍摘要 在全球变化科学研究中，温室效应气体 CO2 在地球系统中的循环成为了人们关注的焦点问题。由于生态系统的不稳定性，使得碳通量的估计具有巨大的复杂性和不确定性，提高碳通量的估计精度是全球变化研究中亟待解决的问题。本文通过对碳通量产生机制的分析，对不同植被类型区域的陆地、海洋经验碳通量采用不同参数调整的优化方案，基于极大似然方法思想、矩阵分块理论和最小二乘思想，提出一种优化已有碳通量估算结果的方法，并通过随机模拟验证了该方法的有效性。进一步，我们将该方法应用于实际数据(GLOBALVIEW- CO2,2005
2、 数据集，以及 Feng Deng 2007 文章中反演的大气传输算子、模型和观测误差方差矩阵)，改进了碳通量的估计效果。关键字 碳通量 大气反演 最小二乘估计1. 引言碳循环是地球上最大、最重要的生物地球化学循环。在工业革命之前，全球生物圈与大气系统处于平衡之中，碳循环也处于动态平衡状态。但是，工业革命 开始，特别是近几十年以来，人类对于物质资料和能源的需求剧增，包括土地利用的改变，植被和林木的破坏，化石燃料的大量使用使得大气中 CO2 浓度显著增加，打破了地球生物圈和大气圈之间的动态平衡，引发了全球气候变暖，并带来了冰川融化，海平面上升等一系列的重大环境问题，对人类的生存环境造成了极大的影
3、响。在非传统安全领域，气候变化问题无疑是近几年国际舞台的一大焦点。陆地碳汇因其对大气 CO2 浓度的重要调节作用，以列入京都议定书 ， “巴厘岛线路图”及“哥本哈根世界气候大会”等国际气候谈判的重要内容。中国独特的地理形式，使其对全球气候的影响不言而喻，同时我国还面临着巨大的温室气体减排限排国际压力，因此，中国不仅需对自己的碳排放有清醒的认识，而且对全球碳排放必须有足够深入的了解，这不仅是自然科学发展的必然过程，且更具有更深远的政治意义。限于自然变化和人为活动对碳循环的双重驱动机制，以及地气系统碳交换和大气传输的复杂性，使得人们在人类活动与自然因素对碳源汇影响的科学认识上还有局限性。因此近几十
4、年来，国内外碳循环研究已经引起了高度的重视。目前用来估计区域和全球碳通量的主要方法有 3 种：直接测量 12，生态系统模拟 34567, 大气反演模型 89101112，这些方法各具优势和不足。直接测量是指利用地面多站点或单站点多塔联网的碳水通量观测，通过插值获取全球碳通量的分布情况的方法。虽然全球已有 500 多个通量观测站，但是由于生态系统特征的空间变异性，将点的碳通量观测数据简单外推存在很大的误差。生态系统模型是利用站点碳通量观测和陆地生态系统的光合作用机制建模，可以用来模拟陆地生态系统碳循环过程的变化特征以及预测未来趋势，这种模型能很好地拟合观测站点局部的碳循环规律，相对于直接测量法已
5、经有所改进，但应用到其他区域会存在很大的误差。模型结果的可靠性受模型结构、参数取值和输入数据质量的影响，模型参数的误差可能带来模型结果的系统性偏差。大气反演是一种利用大气 CO2 浓度观测提取陆地以及海洋碳通量的方法。它充分利用了大气 CO2 浓度变化对生态系统碳循环变化响应的信息，可以较精确地追踪出大尺度区域陆地生态系统碳汇的年际变化特征 13。早期的大气反演只将全球分为 11 个陆地区和 11 个海洋区，时间步长年或月 1415。随后出现了以地理统计学为基础利用多元地表参数来直接限制大气反演结果的方法 16。从减少设置侧边界条件带来的不确定性角度考虑，2007 年 Deng 等 17人发展
6、了嵌套式的大气反演方法，以生态系统模型的输出结果作为经验碳通量。于此同时，高时空分辨率的区域大气反演方法开始出现 18。但大气反演结果的可靠性始终取决于大气 CO2 浓度观测数据的密度和大气输送模型的精度，全球观测站点的稀疏在一定程度上也限制着大气反演方法的精度。由于生态系统模型存在着很大的误差，使得生态系统模型的输出结果精度较低，从而使得由它作为经验通量得到的贝叶斯估计结果精度较低。因此我们需要对经验碳通量进行参数调整，利用碳浓度观测数据优化模型参数，以提高估计精度。同时，在参数优化过程中，通常需要考虑不同植被功能类型(如常绿针叶林、落叶针叶林、常绿阔叶林、落叶阔叶林、灌木、草地、农田和湿地
7、等),对不同植被功能类型进行不同参数的调整,以减少模型的系统偏差。本文利用 CO2 浓度观测数据优化地表 CO2 通量生态系统模型参数，基于Peters(2007，2009) 1920中的碳通量模型，提出了一种新的参数优化估计方法，数据模拟表明该方法有很高的估计精度。2. 碳通量相关模型和结论陆气同化优化碳通量估计模型如下：* MERGEFORMAT (2.1)0sxcGf(,)其中模型分别表示地表生态系统碳通量回归模型和大气反演碳通量模型； 为s碳通量向量，它的维数与地表分割的区域数相同，各分量代表了特定区域的碳通量； 为地表生态系统碳通量的回归函数； 为影响碳通量的协变量；()fx x为大
8、气二氧化碳浓度观测向量（三维空间） ，维数与观测站点数相同； 为大oc G气反演碳通量模型中的响应矩阵；模型误差 与 相互独立且。v(,)0,E=,0在地表生态系统模型中，影响地表碳通量的因素包括陆地碳通量、海洋碳通量 、石化燃料碳排放 和火灾碳排放 。其ERGP()sOCE()sF()sFIRE()s中 是从海洋生态系统模型中得到， 是从陆地生态系统模型中得到。OC ERGPs考虑到陆面生态系统模型中存在着误差，本文基于 Peters(2007，2009) 1920建立的碳通量回归模型：* MERGEFORMAT ERGPOCEFIRE(,)(+)xsssAf(2.2)该模型中, 是一个向量
9、，维数与 、 相同；“ ”表示向量内积。例如当ERsGPA时， 。这里石化燃料碳12(,),ku 12(,)kv 12(,)kuvvA排放和火灾碳排放是固定的，基于自下而上的估计方法得到，而对陆地和海洋碳通量进行参数调整。这一点上，一方面体现了对石化燃料燃烧碳排放数据的信心，另一方面也表明对陆地海洋碳通量有修正作用的大气观测数据的不足。在模型估计中，我们考虑 种植被功能类型，得到如下的模型：5* ()ERGPOCEFIRE1,2345iiiiiiiosssic MERGEFORMAT (2.3)其中 为第 i 种植被类型地块碳通量， ； 为is ,iTT12345(,)ss各个 拉直后连接成的
10、碳通量数据；而各个 连接的 为i i()()()待估模型参数，模型误差 与观测误差 相互独立且i(1,2345)。cov(,)0,E=,0i i根据下述定理，得到参数估计结果。定理 1 模型* MERGEFORMAT (2.3)中如果模型误差 和观测误差 服从正态i 分布，则参数 的最小二乘估计是(1)2110(5)()()XRGQXRGQ c（ 的具体形式及定理证明见附录)X3. 应用案例3.1数据简介 建模数据来自 Deng(2007)17的结果a) 观测 维数 ,来自 GLOBALVIEW-CO2,2005 数据集；碳通量0c5439m的维数 (5(年)*12(月)*50( 区域)，其中
11、 50 个区域的划分参考snDeng(2007)17 (50 个区域的划分如图 1)。碳通量 3000 维的数据排列顺序是第 1 个月 1-50 个区域，第 2 个月 1-50 个区域，第 3 个月 1-50 个区域，第 60 个月 1-50 个区域；图 1：北美 30 个区域划分和其余的全球 20 个区域划分b) 三个矩阵数据 大气传输算子 是五年(1999-2003)50个区域的月度传输矩54390*G阵,通过结合Biome-BGC 21模型(这是一个由来自VCEP/NCAR分析数据的每天气象数据驱动的生态系统模型)计算得到； 观测误差方差矩阵矩阵 ，其中第 个月的误差标准差定义(5439
12、*)Ri如下： , 表示残差分布的标准差，它22017.pmvGVsdis是通过对GLOBALVIEW-CO 2 2005中的平均月变差数据逐月进行计算得到； 为系统误差，所有站点数据均采用这个值。0175.pmv 模型误差方差矩阵 是以对角矩阵的形式构造，其中对于54390*G20 个大区域，使用的是 TransCom 3 年际差异版本中的先验通量不确定性，对于 30 个小区域，则采用 TransCom 3 季度反演中估算大区域不确定性的方法 Gurney(2004)7，分别重新计算它们的不确定性。c) 分别服从正态分布,N(,),;0RQ3.2 实验结果基于上面的数据和本文提出的优化方法，
13、我们对经验碳通量进行参数调整，得到的参数估计结果：表 1：参数估计结果估计结果(1)(2)(3)(4)(5)估计值 1.7952 0.8412 1.3420 1.0991 0.8206从表 1 的参数估计结果看出，除 估计值比较接近 1，其余参数的估计值(4)和 1 相差较大。因此，如果不进行参数调整（一般模型处理中认为 ()1,i） ，观测数据并不能使得目标函数 极小化。从统计似然观点看，,5i (,)Js参数均为 1 的模型并不能使得观测出现的可能性最大，然而气象学中的观测资料具有较好的公信度，所以认为模型（2.1 式）的参数不全为 1。同时，分别通过由参数调整和不进行参数调整两种模型得到
14、的后验碳通量计算出 CO2 浓度，将它们与原始 CO2 观测数据进行比较，得到参数调整的 CO2浓度残差均方根误差(1.4030)小于未进行参数调整的残差均方根误差(1.4104).从这个角度也说明对经验碳通量进行参数调整对于提高碳通量估计准确度具有很大的作用。4. 估计方法模拟验证和精度估计我们需要验证这种估计方法的有效性以及研究它的估计精度，根据假设检验的思想以及基于中心极限定理保证，采用统计模拟的方法为模型估计精度。模拟实验前首先需要取定一组参数（也就是 的真值） ，这组参数在之后的模拟过程中都是固定不变的。由于在优化过程中，需要通过观测 和经验碳通量oc来估计参数。因此每次模拟都需要分
15、别模拟真实的ERGPOCEFIREsss碳通量和二氧化碳浓度观测值。因此模拟过程总体上分为 2 个部分，他们之间的关系如下图所示（矩阵 数据参考第 3 节）Q、ERsAGPsOCEsAFIREs*ERGPOCEFIRE()psssAERGPOCEFIRE()sssA*TT1*1()()()2popoJscsc+QERGPOCEFIRE,sssERGPOC/3/so图 2：碳通量参数优化模拟数据流程图4.1. 真实通量、模型通量、观测浓度的模拟方法真实通量模拟方法：生态系统呼吸碳通量 产生于1,5均匀分布随机数；ERs生态系统总初级生产力碳通量 产生于1,4均匀分布随机数 . 海洋碳通量GPs产
16、生于1,6 均匀分布随机数；真实火灾和人为排放碳通量 为3000OCEs FIRE,s维的0向量，在模拟中，我们暂时不考虑这两项的值，于是先假设为0，在后面的叙述中，如无特别说明，都默认这两项为0。模型真实通量 通过碳通量回归函数确定:sERGPOCEFIRE()sssA二氧化碳浓度观测值是通过大气反演碳通量模型 模拟得到。oc4.2. 参数及均方根误差的估计方法对于调整参数的估计1) 先对真实通量加以 为方差的扰动误差，模拟得到经验碳通量Q,假设模型参数 未知，则碳通量估计ERGPOCE(,)s*RGPOCEFIRE()psssA2) 结合第一步得到的二氧化碳浓度观测数据 和经验碳通量，通过
17、极小oc化目标函数* MERGEFORMAT 错误！未找到引用源。，得到 的估计值 。重复模拟以上步骤 次，得到估计值, 计算估计量的均方根误差k()k()21RMSE=i4.3. 模拟结果模拟中我们分别选取两组真值 = ， 模拟结果(,2345)(0.8,1.2,3)和如下：表2: 估计结果(模拟次数 )k估计结果(1)(2)(3)(4)(5)真值 1 2 3 4 5估计均值 1.004 1.9999 2.9998 3.9896 4.9903估计 RMSE 0.0015 0.0027 0.0017 0.0159 0.0153估计结果(1)(2)(3)(4)(5)真值 0.8 1.5 1 0.
18、2 3估计均值 0.8002 1.4998 1.0000 0.2010 2.9931估计 RMSE 0.0010 0.0018 0.0008 0.0028 0.0099从表2的结果可以看出统计模拟的估计值非常接近真值，说明估计相当准确；同时估计的RMSE大部分在0.003以下，也体现出估计结果精度是很高的。因此，从统计模拟结果可以看出：文章所提的新的参数估计方案可靠。 5. 结论与讨论在碳通量研究发展过程中，提出了一些碳通量估计方法(如直接测量、生态系统模型估计、大气反演等)，但是碳通量的估计准确性，仍然是一个问题。本文通过对碳通量模型中的先验碳通量进行参数调整，利用一种新的估计方法,并将其应
19、用于实际模型(Deng(2007) 3)中,结果表明对经验碳通量进行参数调整是十分必要的。同时利用随机模拟对估计方法的估计精度进行验证，得到了很高的估计精度。这样今后将这种方法融合到同化过程中，利用经过优化的模型进行数据同化，不论是在对历史碳汇的模拟还是对在未来碳汇的预测，都将得到更为可信的结果。附录定理 1 证明：先验碳通量具有形式 ， 其中*ERGPSOEFIRE+)psss（010101010101010101( ) 由于误差服从正态分布，因此最优的 是使得目标函数T1T1(,)()R()()()22ooppJssccssGsQ达到最小的参数值，这里的目标函数 是关于 和 两个未知量的函
20、数，,J因此先固定 ，上述目标函数关于 的最小值解，即优化的生态系统碳通量为s-1-1T*0T()+)()(QGGRpscs从而*-1-1T*0()(+)()GRQGp ps csT0 01*0()()pc cs因此我们需要极小化目标函数为 T1 *T1*00*T-1 0T-1TT1*0*T-1*1()()R(G)()()22()()()()1()2 ppp pp pJscscsscsc ss QQRGR 0T-TT10*T-1*0T*0*T-1TT1*0 0T*0()()()()()1)()()2()p ppp ppcsc scscsc ssc GRQQGRG为了得到参数估计的显示解，我们将矩阵、向量按照时间、类型进行分块： (1)2)(5)(12)(5)(1)2)(5)6060 GG1122606011226T()()T()()T()()TERERERER()(5)()(5)(1)(5)GPPGPGPPGPT()()T()()TOCEEOCEOCEsssssss 060()T()OCEs 令 ()()()12601,34,5iiii12601260T()()T()ERER()()()GPPGPT()()T()OCEEOCE ,2, 1,34,5iiiiiiiiiiiississi
