利用CO2浓度观测提高地表碳通量模型估算精度的方法.doc

利用CO2浓度观测提高地表碳通量模型估算精度的方法北京师范大学 郑珩、黄文贤、黄勍摘要 在全球变化科学研究中，温室效应气体CO2在地球系统中的循环成为了人们关注的焦点问题。由于生态系统的不稳定性，使得碳通量的估计具有巨大的复杂性和不确定性，提高碳通量的估计精度是全球变化研究中亟待解决的问题。本文通过对碳通量产生机制的分析，对不同植被类型区域的陆地、海洋经验碳通量采用不同参数调整的优化方案，基于极大似然方法思想、矩阵分块理论和最小二乘思想，提出一种优化已有碳通量估算结果的方法，并通过随机模拟验证了该方法的有效性。进一步，我们将该方法应用于实际数据(GLOBALVIEW- CO2,2005数据集，以及Feng Deng 2007文章中反演的大气传输算子、模型和观测误差方差矩阵)，改进了碳通量的估计效果。关键字 碳通量 大气反演 最小二乘估计1. 引言碳循环是地球上最大、最重要的生物地球化学循环。在工业革命之前，全球生物圈与大气系统处于平衡之中，碳循环也处于动态平衡状态。但是，工业革命 开始，特别是近几十年以来，人类对于物质资料和能源的需求剧增，包括土地利用的改变，植被和林木的破坏，化石燃料的大量使用使得大气中CO2浓度显著增加，打破了地球生物圈和大气圈之间的动态平衡，引发了全球气候变暖，并带来了冰川融化，海平面上升等一系列的重大环境问题，对人类的生存环境造成了极大的影响。在非传统安全领域，气候变化问题无疑是近几年国际舞台的一大焦点。陆地碳汇因其对大气CO2浓度的重要调节作用，以列入京都议定书，“巴厘岛线路图”及“哥本哈根世界气候大会”等国际气候谈判的重要内容。中国独特的地理形式，使其对全球气候的影响不言而喻，同时我国还面临着巨大的温室气体减排限排国际压力，因此，中国不仅需对自己的碳排放有清醒的认识，而且对全球碳排放必须有足够深入的了解，这不仅是自然科学发展的必然过程，且更具有更深远的政治意义。限于自然变化和人为活动对碳循环的双重驱动机制，以及地气系统碳交换和大气传输的复杂性，使得人们在人类活动与自然因素对碳源汇影响的科学认识上还有局限性。因此近几十年来，国内外碳循环研究已经引起了高度的重视。 目前用来估计区域和全球碳通量的主要方法有3种：直接测量12，生态系统模拟34567, 大气反演模型89101112，这些方法各具优势和不足。直接测量是指利用地面多站点或单站点多塔联网的碳水通量观测，通过插值获取全球碳通量的分布情况的方法。虽然全球已有500多个通量观测站，但是由于生态系统特征的空间变异性，将点的碳通量观测数据简单外推存在很大的误差。生态系统模型是利用站点碳通量观测和陆地生态系统的光合作用机制建模，可以用来模拟陆地生态系统碳循环过程的变化特征以及预测未来趋势，这种模型能很好地拟合观测站点局部的碳循环规律，相对于直接测量法已经有所改进，但应用到其他区域会存在很大的误差。模型结果的可靠性受模型结构、参数取值和输入数据质量的影响，模型参数的误差可能带来模型结果的系统性偏差。 大气反演是一种利用大气CO2浓度观测提取陆地以及海洋碳通量的方法。它充分利用了大气CO2浓度变化对生态系统碳循环变化响应的信息，可以较精确地追踪出大尺度区域陆地生态系统碳汇的年际变化特征13。早期的大气反演只将全球分为11个陆地区和11个海洋区，时间步长年或月1415。随后出现了以地理统计学为基础利用多元地表参数来直接限制大气反演结果的方法16。从减少设置侧边界条件带来的不确定性角度考虑，2007年Deng等 17人发展了嵌套式的大气反演方法，以生态系统模型的输出结果作为经验碳通量。于此同时，高时空分辨率的区域大气反演方法开始出现18。但大气反演结果的可靠性始终取决于大气CO2浓度观测数据的密度和大气输送模型的精度，全球观测站点的稀疏在一定程度上也限制着大气反演方法的精度。由于生态系统模型存在着很大的误差，使得生态系统模型的输出结果精度较低，从而使得由它作为经验通量得到的贝叶斯估计结果精度较低。因此我们需要对经验碳通量进行参数调整，利用碳浓度观测数据优化模型参数，以提高估计精度。同时，在参数优化过程中，通常需要考虑不同植被功能类型(如常绿针叶林、落叶针叶林、常绿阔叶林、落叶阔叶林、灌木、草地、农田和湿地等),对不同植被功能类型进行不同参数的调整,以减少模型的系统偏差。 本文利用CO2浓度观测数据优化地表CO2通量生态系统模型参数，基于Peters(2007，2009)1920中的碳通量模型，提出了一种新的参数优化估计方法，数据模拟表明该方法有很高的估计精度。2. 碳通量相关模型和结论陆气同化优化碳通量估计模型如下：其中模型分别表示地表生态系统碳通量回归模型和大气反演碳通量模型；为碳通量向量，它的维数与地表分割的区域数相同，各分量代表了特定区域的碳通量；为地表生态系统碳通量的回归函数；为影响碳通量的协变量；为大气二氧化碳浓度观测向量（三维空间），维数与观测站点数相同；为大气反演碳通量模型中的响应矩阵；模型误差与相互独立且。在地表生态系统模型中，影响地表碳通量的因素包括陆地碳通量、海洋碳通量、石化燃料碳排放和火灾碳排放。其中是从海洋生态系统模型中得到，是从陆地生态系统模型中得到。 考虑到陆面生态系统模型中存在着误差，本文基于Peters(2007，2009)1920建立的碳通量回归模型：该模型中,是一个向量，维数与、相同；“”表示向量内积。例如当时，。这里石化燃料碳排放和火灾碳排放是固定的，基于自下而上的估计方法得到，而对陆地和海洋碳通量进行参数调整。这一点上，一方面体现了对石化燃料燃烧碳排放数据的信心，另一方面也表明对陆地海洋碳通量有修正作用的大气观测数据的不足。 在模型估计中，我们考虑种植被功能类型，得到如下的模型： 其中为第i种植被类型地块碳通量，；为各个拉直后连接成的碳通量数据；而各个连接的为待估模型参数，模型误差与观测误差相互独立且。根据下述定理，得到参数估计结果。定理1 模型中如果模型误差和观测误差服从正态分布，则参数的最小二乘估计是（的具体形式及定理证明见附录)3. 应用案例3.1 数据简介 建模数据来自Deng(2007)17的结果a) 观测维数,来自GLOBALVIEW-CO2,2005数据集；碳通量的维数(5(年)*12(月)*50(区域)，其中50个区域的划分参考Deng(2007)17 (50个区域的划分如图1)。碳通量3000维的数据排列顺序是第1个月1-50个区域，第2个月1-50个区域，第3个月1-50个区域，第60个月1-50个区域；图1：北美30个区域划分和其余的全球20个区域划分b) 三个矩阵数据 大气传输算子是五年(1999-2003)50个区域的月度传输矩阵,通过结合Biome-BGC21模型(这是一个由来自VCEP/NCAR分析数据的每天气象数据驱动的生态系统模型)计算得到； 观测误差方差矩阵矩阵，其中第个月的误差标准差定义如下：, 表示残差分布的标准差，它是通过对GLOBALVIEW-CO2 2005中的平均月变差数据逐月进行计算得到；为系统误差，所有站点数据均采用这个值。 模型误差方差矩阵是以对角矩阵的形式构造，其中对于20个大区域，使用的是TransCom 3 年际差异版本中的先验通量不确定性，对于30个小区域，则采用TransCom 3 季度反演中估算大区域不确定性的方法Gurney(2004)7，分别重新计算它们的不确定性。c) 分别服从正态分布3.2 实验结果基于上面的数据和本文提出的优化方法，我们对经验碳通量进行参数调整，得到的参数估计结果：表1：参数估计结果估计结果估计值1.79520.84121.34201.09910.8206 从表1的参数估计结果看出，除估计值比较接近1，其余参数的估计值和1相差较大。因此，如果不进行参数调整（一般模型处理中认为），观测数据并不能使得目标函数极小化。从统计似然观点看，参数均为1的模型并不能使得观测出现的可能性最大，然而气象学中的观测资料具有较好的公信度，所以认为模型（2.1式）的参数不全为1。 同时，分别通过由参数调整和不进行参数调整两种模型得到的后验碳通量计算出CO2浓度，将它们与原始CO2观测数据进行比较，得到参数调整的CO2浓度残差均方根误差(1.4030)小于未进行参数调整的残差均方根误差(1.4104).从这个角度也说明对经验碳通量进行参数调整对于提高碳通量估计准确度具有很大的作用。4. 估计方法模拟验证和精度估计我们需要验证这种估计方法的有效性以及研究它的估计精度，根据假设检验的思想以及基于中心极限定理保证，采用统计模拟的方法为模型估计精度。模拟实验前首先需要取定一组参数（也就是的真值），这组参数在之后的模拟过程中都是固定不变的。由于在优化过程中，需要通过观测和经验碳通量来估计参数。因此每次模拟都需要分别模拟真实的碳通量和二氧化碳浓度观测值。因此模拟过程总体上分为2个部分，他们之间的关系如下图所示（矩阵数据参考第3节）图2：碳通量参数优化模拟数据流程图4.1. 真实通量、模型通量、观测浓度的模拟方法 真实通量模拟方法：生态系统呼吸碳通量产生于1,5均匀分布随机数；生态系统总初级生产力碳通量产生于1,4均匀分布随机数. 海洋碳通量产生于1,6均匀分布随机数；真实火灾和人为排放碳通量为3000维的0向量，在模拟中，我们暂时不考虑这两项的值，于是先假设为0，在后面的叙述中，如无特别说明，都默认这两项为0。模型真实通量通过碳通量回归函数确定:二氧化碳浓度观测值是通过大气反演碳通量模型模拟得到。4.2. 参数及均方根误差的估计方法对于调整参数的估计1) 先对真实通量加以为方差的扰动误差，模拟得到经验碳通量,假设模型参数未知，则碳通量估计2) 结合第一步得到的二氧化碳浓度观测数据和经验碳通量，通过极小化目标函数，得到的估计值。重复模拟以上步骤次，得到估计值,计算估计量的均方根误差4.3. 模拟结果模拟中我们分别选取两组真值=，模拟结果如下：表2:估计结果(模拟次数) 估计结果真值12345估计均值1.0041.99992.99983.98964.9903估计RMSE0.00150.00270.00170.01590.0153估计结果真值0.81.510.23估计均值0.80021.49981.00000.20102.9931估计RMSE0.00100.00180.00080.00280.0099 从表2的结果可以看出统计模拟的估计值非常接近真值，说明估计相当准确；同时估计的RMSE大部分在0.003以下，也体现出估计结果精度是很高的。因此，从统计模拟结果可以看出：文章所提的新的参数估计方案可靠。 5. 结论与讨论 在碳通量研究发展过程中，提出了一些碳通量估计方法(如直接测量、生态系统模型估计、大气反演等)，但是碳通量的估计准确性，仍然是一个问题。本文通过对碳通量模型中的先验碳通量进行参数调整，利用一种新的估计方法,并将其应用于实际模型(Deng(2007)3)中,结果表明对经验碳通量进行参数调整是十分必要的。同时利用随机模拟对估计方法的估计精度进行验证，得到了很高的估计精度。这样今后将这种方法融合到同化过程中，利用经过优化的模型进行数据同化，不论是在对历史碳汇的模拟还是对在未来碳汇的预测，都将得到更为可信的结果。附录定理1证明：先验碳通量具有形式，其中由于误差服从正态分布，因此最优的是使得目标函数达到最小的参数值，这里的目标函数是关于和两个未知量的函数，因此先固定，上述目标函数关于的最小值解，即优化的生态系统碳通量为从而因此我们需要极小化目标函数为 为了得到参数估计的显示解，我们将矩阵、向量按照时间、类型进行分块：令取因此极小化的目标函数具有下面形式令，则可以看成是线性模型的广义最小二乘估计目标函数的0.5倍，因此得到最优估计是参考文献1 Wofsy, S. C., Goulden, M. L., Munger, J. W., Fan, S. M., Bakwin. P. S., and co-authors. 1993. Net exchange of CO2 in a midlatitude forest.Science 260,1314.2 Takahashi, T., Wanninkhof, R., Feely, R. A., Weiss, R. F., Chipman, D. W., and co-authors. 1999. Net sea-air CO2 flux over the global oceans:An improved estimate based on the sea-air pCO2 difference.In Proceedings of the 2nd International Symposium on CO2 in the Oceans,the 12th global environment Tsukuba,18-22 January 1999.Tsukuba,Japan,18-01,4 pp.3 Sellers, P. J., Mintz, Y., Sud. Y. C., and Dalcher, A. 1986. A simple biosphere model(SiB)for use within general circulation models. J. Atmos.Sci. 43, 505-531.4 Sellers, P. J., Randall. D. A., Collatz, G. J., Berry, J. A., Field, C. B., and coauthors. 1996. A revised land surface parameterization(SiB2) for atmospheric GCMs .l. Model formulation.J.Climate 9, 706-737.5 Running, S. W. and Gower, S. T. 1991. FOREST-BGC, a general model of forest ecosystem processes for regional applications,II. Dynamic carbon allocation and nitrogen budgets. Tree Physiol. 9, 147-160.6 Running, S. W. and Coughlan, J. C. 1998. A general model of forest ecosystem processes for regional applications, I. hydrologic balance, canopy gas exchange and primary production processes. Ecol. Model. 42, 125-154.7 Liu, J., Chen, J.M., and Chen, W. 1999. Net primary productivity distribution in BOREAS region from a process model driven by statellite and surface data. J.Geophy.Res. 104, 27735-27754.118 Kaminski, T., Heimann, M. and Giering, R.1999. A coarse grid threedimensional global inverse model of the atmospheric transport:1.Adjoint model and Jacobian matrix.J.Geophys. Res. 104(D15), 18535-18554.9 Rdenback, C., houweling, S., Gloor, M.and Heimann, M. 2003. CO2 flux history 1982-2001 inferred from atmospheric data using a global inversion of atmospheric transport. Atmos. Chem.Phys. 3. 1919-1964.10 Gurney, K. R., Law, R. M., Denning, A. S., Rayner, P. J., Baker, D., and co-authors. 2003. TransCom3 CO2 inversion intercomparison:1.Annual mean control results and sensitivity to transport and prior flux information. Tellus 55B, 555-579.11 Gurney, K. R., Law, R. M., Denning, A. S., Rayner, P. J.,Baker, D., and coauthors. 2004. TransCom3 CO2 inversion intercomparison:Mcdel mean results for the estimation of seasonal carbon sources and sinks.Global Biogeochem. Cycles 18(1),GB1010.12 Baker, D. F., Law, R, M., Gurney, K., R., Rayner, P., Peylin, and co-author. 2006. TranCom3 inversion intercomparison: Impact of transport model errors on the interannual variability of regional CO2 fluxes, 1998-2003. Global Biogeochem. Cycles:20(1),GB1002.13 Denman, K.L., G. Brasseur,A. Chidthaisong et al.,2007.Couplings Between Changes in the Climate System and Biogeochemistry.In:Climate Change 2007:The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change Solomon, S., D. Qin, M. Manning, Z. Chen, M. Marquis, K. B. Averyt, M. Tignor and H.L.Miller(eds.).Cambridge University Press,Cambridge,United Kingdom and New York,NY,USA.14 Denning AS, M Holzer; KR Gurney et al,1999. Three-dimensional transport and concentration of SF6-A model intercomparison study(TransCom2).Tellus, 52B, 266-297.15 Guiney KR, RM Law,AS Denning et al.,2002.Towards robust regional estimates of CO2 sources and sinkes using atmospheric models.Nature, 415(6872),626-630.16 Michalak A,Bruhwiler, T, and Tans P, 2004. A geostatistical approach to surface flux estimation of atmospheric trac gases. Journal of Geophysical Research, 109(D14),doi:10.1029/2003JD004422.17 Deng F, JM Chen,and Ishizawa, M et al.,2007.Global monthly CO2 flux inversion with a focus over North America .Tellus 59B, 179-190.18 Schuh AE, AS Denning, and M Uliasz et al. 2009. Seeing the forest through the trees:Recovering large-scale carbon flux biases in the midst of small-scale variability. Journal of Geophysical Research, 114,G3,doi:10.1029/2008JG000842, 2009.19 Peter W, AR Jackson, and C Sweeney,2007. An atmospheric perspective on North American carbon dioxide exchange:CarbonTracker.PNAS,104(48):18925-18930.20 Peter W, MC Krol,and van der Werf G.R et al.,2009.Seven years of recent Europenan net terrestrial carbon dioxide exchange constrained by atmospheric observations. Global Change Biology,16(4), 1317-1337. A Method to Improve The Accuracy of Surface Carbon Flux Estimate Using CO2 Concentration ObservationAbstract In nowadays study of global environment change,