【毕业学位论文】（Word原稿）利用CO2浓度观测提高地表碳通量模型估算精度的方法-统计教育学

利用 北京师范大学 摘要 在全球变化科学研究中，温室效应气体 地球系统中的循环成为了人们关注的焦点问题。 由于生态系统的不稳定性， 使得碳通量的估计具有巨大的复杂性和不确定性，提高碳通量的估计精度是全球变化研究中亟待解决的问题。本文通过 对 碳通量产生机制的分析， 对不同植被类型区域 的 陆地、海洋 经验碳通量采用不同参数调整的 优化方案 ， 基于极大似然方法思想 、矩阵分块理论和最小二乘思想，提出一种优化已有碳通量估算结果的方法，并通过随机模拟验证了该方法的有效性。进一 步，我们将该方法应用于实际数据 (005数据集，以及 007 文章中反演的大气传输算子、模型和观测误差方差矩阵 )，改进了碳通量的估计效果。 关键字 碳通量 大气反演 最小二乘估计 1. 引言 碳循环是地球上最大、最重要的生物地球化学循环。 在工业革命之前，全球生物圈与 大气系统处于平衡之中，碳循环也处于动态平衡状态。但是，工业革命 开始 ，特别是近几十年以来，人类对于物质资料和能源的需求剧增，包括土地利用的改变，植被和林木的破坏，化石燃料的大量使用 使得大气中 度显著增加，打破了地球生物圈和大气圈之间的动态平衡，引发了全球气候变暖，并带来了冰川融化，海平面上升等一系列的重大环境问题， 对人类的生存环境造成了极大的影响。 在非传统安全领域，气候变化问题无疑是近几年国际舞台的一大焦点。 陆地碳汇因其对大气 度的重要调节作用，以列入京都议定书，“巴厘岛线路图”及“哥本哈根世界气候大会”等国际气候谈判的重要内容。中国独特的地理形式，使其 对全球气候的影响不言而喻， 同时我国还面临着巨大的温室气体减排限排国际压力， 因此，中国不仅 需对自己的碳排放有清醒的认 识，而且对全球碳排放必须有足够深入的了解，这不仅是自然科学发展的必然过程，且更具有更 深远的政治意义。 限于 自然变化和人为活动对碳循环的双重驱动机制，以及地气系统碳交换和大气传输的复杂性，使得人们在人类活动与自然因素对碳源汇影响的科学认识上还有局限性。因此近几十年来，国内外碳循环研究已经引起了高度的重视。 目前用来估计区域和全球碳通量的主要方法有 3 种：直接测量 12，生态系统模拟 34567, 大气反演模型 89101112， 这些方法各具优势和不足。 直接测量是指利 用地面 多站点或单 站 点多塔联网的碳水通量观测 ，通过插值获取全球碳通量的分布情况的方法 。 虽然全球已有 500 多个通量观测站，但是由于生态系统特征的空间变异性，将点的碳通量观测数 据简单 外推 存在很大的误差 。 生态系统模型 是利用站点碳通量观测和陆地生态系统的光合作用机制建模，可以用来模拟陆地生态系统碳循环过程的变化特征以及预测未来趋势， 这种模型能很好地拟合观测站点局部的碳循环规律，相对于直接测量法已经有所改进，但应用到其他区域会存在很大的误差。 模型结果的可靠性受模型结构、参数取值和输入数据质量的影响，模型参数的误差可能 带来模型结果的系统性偏差。 大气反演是一种利用大气 度观测提取陆地以及海洋碳通量的方法。它充分利用了大气 度变化对生态系统碳循环变化响应的信息，可以较精确 地追踪出 大尺度 区域陆地生态系统碳汇的年际变化特征 13。早期的大气反演只将全球分为 11 个陆地区和 11 个海洋区，时间步长年或月 1415。随后出现了以地理统计学为基础利用多元地表参数来直接限制大气反演结果的方法 16。 从 减少设置侧边界条件带来的不确定性 角度考虑 ， 2007 年 17人发展了嵌套式的大气反演方法， 以生态 系统模型的输出结果作为经验碳通量。 于此同时，高时空分辨率的区域大气反演方法开始出现 18。 但大气反演结果的可靠性始终取决于大气 度观测数据的密度和大气输送模型的精度 ， 全球观测站点的稀疏在一定程度上也限制着大气反演方法的精度。 由于生态系统模型存在着很大的误差，使得生态系统模型的输出结果精度较低，从而使得由它作为经验通量得到的贝叶斯估计结果精度较低。因此我们需要对经验碳通量进行参数调整，利用碳浓度观测数据优化模型参数，以提高估计精度。同时， 在参数优化过程中，通常需要考虑不同植被功能类型 (如常绿针叶林、落叶针叶林、常绿阔叶林、落叶阔叶林、灌木、草地、农田和湿地等 ),对不同植被功能类型进行不同参数的调整 ,以减少模型的系统偏差。 本文利用 度观测数据优化 地表 量生态系统模型 参数，基于007， 2009)1920中的碳通量模型，提 出了一种新 的参数优化估计方法，数据模拟表明该方法有很高的估计精度。 2. 碳通量相关模型和结论 陆气同化优化碳通量估计模型如下： 0 , ) (其中模型分别表示地表生态系统碳通量回归模型和大气 反演碳通量 模型 ； s 为碳通量向量，它的维数与地表分割的区域数相同，各分量代表了特定 区域的碳通量；()f x 为地表生态系统碳通量的回归函数； x 为影响碳通量的协变量； 大气二氧化碳浓度观测向量（三维空间），维数与观测站点数相同； G 为大气 反演碳通量 模 型 中 的 响 应 矩 阵 ； 模 型 误 差 与 相 互 独 立 且c o v ( , ) 0 , E = 0 , E 0 。 在地表生态 系统 模型中， 影响 地表碳通量 的因素 包括陆地碳通量洋碳通量s、石化燃料碳排放s。其中R 考虑到陆面生态系统模型中 存在着误差， 本文基于 007， 2009)1920建立 的碳 通量 回归 模型 ： E R G P P O C E F F F I R E( , ) ( + )x s s s s s f(该模型中 , 是一个向量，维数与 ”表示向量内积。 例如当12( , , , ) ,ku u uu 12( , , , )kv v v 1 2 2( , , , )v u v u v里石化燃料碳排放和火灾碳排放是固定的，基于自下而上的估计方法得到，而对陆地和海洋碳通量进行参数调整。这一点上，一方面体现了对石化燃料燃烧碳排放数据的信心，另一方面也表明对陆地海洋碳通量有修正作用的大气观测数据的不足。 在模型估计中，我们考虑 5 种 植被功能类型 ，得到如下的模型： () E R G P P O C E F F F I R E 1 , 2 , 3 , 4 , 5i i i i s s s s s s (其中i 种植被类型 地块碳通量， 1, 2, 3, 4, 5i ； T T T T T 3 4 5( , , , )s s s s s s为各个各个i连接的 ( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) 为 待 估模型 参数 ， 模 型 误 差 i ( 1, 2 , 3, 4 , 5)i 与 观 测 误 差 相 互 独 立 且c o v ( , ) 0 , E = 0 , E 0 。 根据下述定理，得到参数估计结果。 定理 1 模型 ( 如果模型误差 服从正态分布，则 参数 的 最小二乘估计 是 ( 1 )( 2 )110( 5 ) ( ( ) ) ( ) X R G Q G X X R G Q G c（ X 的具体形式及定理证明 见附录 ) 3. 应用案例 据简介 建模 数据来自 007)17的结果 a) 观测0439m ,来自 005 数据集；碳通量 439n (5(年 )*12(月 )*50(区域 )， 其中 50 个区域的划分参考007)17 (50个区域的划分如图 1)。 碳通量 3000 维的数据排列顺序是第 1 个月 1区域 ， 第 2 个月 1区域 ， 第 3 个月 1区域 ， ， 第 60 个月 1区域 ； 图 1：北美 30 个区域划分和其余的全球 20 个区域划分 b) 三个矩阵数据 大气传输算子 5 4 3 9 3 0 0 0*G 是五年 (19990个区域的月度传输矩阵 ,通过结合 1模型 (这是一个由来自 计算得到 ； 观测误差方差矩阵 矩阵 ( 5 4 3 9 * 5 4 3 9 )R ， 其中 第 i 个月的误差标准差定义如下： 220 1 7 5. p p m v G V s R, 示残差分布的标准差，它是通过对 005中的平均月变差数据逐月进行计算得到 ； 0 175. 系统误差，所有站点数据均采用这个值。 模型误差方差矩阵 5 4 3 9 3 0 0 0*G 是以对角矩阵的形式构造，其中 对于20 个大区域，使用的是 年际差异版本中的先验通量不确定性，对于 30 个小区域，则采用 季度反演中估算大区域不确定性的方法 004)7，分别重新计算它们的不确定性。 c) , 分别服从正态分布 N ( , ) , N ( , ) ;0 R 0 Q 实验结果 基于上面的 数据和本文提出的优化方法，我们对经验碳通量进行参数调整， 得到的参数估计结果： 表 1：参数估计结果 估计结果 (1) (2) (3) (4) (5) 估计值 表 1 的参数估计结果看出，除 (4) 估计值比较接近 1，其余参数的估计值和 1 相差较大。因此，如果不进行参数调整（一般模型处理中认为()1,i 1, ,5i ），观测数据并不能使得目标函数 ( , )极小化。从统计似然观点看，参数均为 1 的模型并不能使得观测出现的可能性最大，然而气象学中的观测资料具有较好的公信度，所以认为模型（ ）的参数不全为 1。 同时，分别通过由参数调整和不进行参数调整两种模型得到的后验碳通量计算出 度，将它们与原始 测数据进行比较，得到参数调整的 度残差均方根误差 (于未进行参数调整的残差均方根误差 (从这个角度也说明对经验碳通量进行参数调整对于提高碳通量估计准确度具有很大的作用。 4. 估计方法模拟验证和精度估计 我们需要 验证这种估计方法的有效性以及研究它的 估计精度，根据假设检验的思想以及基于中心极限定理保证，采用统计模拟的方法为模型估计精度。模拟实验前首先需要取定一组参数（也就是 的真值），这组参数在之后的模拟过程中都是固定不变的。由于在优化过程中，需要通过观测E R G P P O C E F F F I R s s s s 此每次模拟都需要分别模拟真实的碳通量和二氧化碳浓度观测值。因此模拟过程总体上分为 2 个部分，他们之间的关系如下图所示 （矩阵 、 数据参考第 3 节 ） C F F I R *E R G P P O C E F F F I R ()ps s s s s s E R G P P O C E F F F I R E()s s s s s s G* T T 1 *1( ) ( ) ( ) ( )2p o p oJ s c s c G R + G Q G G P P O C E F F F I R E, , , ,s s s s E /3 /3 /3 ：碳通量参数优化模拟数据流程图 真实通量、模型通量、观测浓度的模拟方法 真实通量模拟方法：生态系统呼吸碳通量 生于 1,5均匀分布随机数； 生态系统总初级生产力碳通量 生于 1,4均匀分布随机数 . 海洋碳通量生于 1,6均匀分布随机数；真实火灾和人 为排放碳通量 3000维的 0向量，在模拟中，我们暂时不考虑这两项的值，于是先假设为 0，在后面的叙述中，如无特别说明，都默认这两项为 0。 模型 真实 通量 s 通过碳通量回归函数确定 : E R G P P O C E F F F I R E()s s s s s s 二氧化碳浓度观测值是通过 大气反演碳通量模型 o 拟得到。 参数及均方根误差的估计方法 对于调整参数的估计 1) 先对真实通量加以 Q 为方差的扰动误差，模拟得到经验碳通量E R G P P O C ( , , )s s s ,假设模型参数 未知，则碳通量估计 * E R G P P O C E F F F I R ()ps s s s s s 2) 结合第一步得到的二氧化碳浓度观测数据 经验碳通量，通过极小化目标函数 错误 !未找到引用源。 ， 得到 的估计值 。 重复模拟以上步骤 k 次，得到估计值 , ()算 估计量的均方根误差 ( ) 211 R M S E = ( )k 模拟结果 模拟中我们分别选取两组真值 = ( 1 , 2 , 3 , 4 , 5 ) ( 0 . 8 , 1 . 5 , 1 , 0 . 2 , 3 )和 ， 模拟结果如下： 表 2: 估计结果 (模拟次数 1000k = ) 估计结果 (1) (2) (3) (4) (5) 真值 1 2 3 4 5 估计均值 计 计结果 (1) (2) (3) (4) (5) 真值 估计均值 计 从表 2的结果可以看出统计模拟的估计值非常接近真值，说明估计相当准确；同时估计的 体现出 估计结果精度 是 很高 的。因此，从统计模拟结果可以看出：文章所提的新的参数估计方案可靠。 5. 结论与讨论 在碳通量研究发展过程中 ，提出了一些碳通量估计方法 (如直接测量、生态系统模型估计、大气反演等 )，但是碳通量的估计准确性，仍然是一个问题。本文通过对 碳通量模型中的 先验碳通量进行参数调整，利用一种新的估计方法 ,并将其应用于实际模型 (007)3)中 ,结果表明对经验碳通量进行参数调整是十分必要的。 同时利用随机模拟对估计方法的估计精度进行验证，得到了很高的估计精度。这样今后将这种方法融合到同化过程中，利用经过优化的模型进行数据同化，不论是在对历史碳汇的模拟还是对在未来碳汇的预测，都将得到更为可信的结果。 附录 定理 1 证明 ： 先验碳通量具有形式 *E R G P P S O E F I R + ) +ps s s s s s ，其中 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 01 0 1 0 1 0( ) 1 1 1 1 1 1 1 1 1111 由于误差服从正态分布，因此最优的 是使得目标函数 T 1 T 111( , ) ( ) R ( ) ( ) ( )22o o p pJ s s c c s s s G G s Q s 达到最小的参数值，这 里的目标函数 ( , )是关于 s 和 两个未知量的函数，因此先固定 ， 上述目标函数关于 s 的最小值解，即优化的生态系统碳通量为 - 1 - 1 1 T 1 * *0T T 1 * *0 ( ) ( + ) ( )( ) ( ) G T R Q G R R G Q G c s sc s s 从而 * - 1 - 1 1 T 1 *0 ( ) ( + ) ( ) G T R Q G R s c s T T 1 * *0 0 0T 1 *0 ( ) ( ) ( )( ) ( ) G G Q G R G Q G G G Q c c s s cG s T 1 * T 1 *00* T T - 1 1 T 1 *00* T T - 1 1 T T 1 *00* T T - 1 T 1 *011 ( ) ( ( ) ) R ( G ( ) ) ( ( ) ) ( ( ) )221( ) ( ) ( ) ( )21( ) ( ) ( ) ( )21( ) ( ) ( ) (2s c s c s s s ss c s cs c s cs c s G Q G R R R R G Q G G Q G G Q Q Q G R G Q G G Q G R R G Q G T T - 1 T T 1 *00* T T - 1 T 1 *00T T 1 *0* T T - 1 T T 1 *00* T T 1 *00)1( ) ( ) ( ) ( )21( ) ( ) ( ( ) ( )2( ) ( )1( ) ( ) ( ) ( ) ( )21( ) ( ) ( )2c s cs c s c s cs c s c G R G Q G G Q G R G Q G G Q G R R G Q G G R G Q G G Q G R G Q G R G Q G G Q G 们将矩阵、向量按照时间、类型进行分块： ( 1 ) ( 2 ) ( 5 ) ( 1 ) ( 2 ) ( 5 ) ( 1 ) ( 2 ) ( 5 )1 1 1 2 2 2 6 0 6 0 6 0()G G G G G G G G G 2 2 6 0 6 01 1 2 2 6 0 6 01 1 2 2 6T ( 1 ) T ( 5 ) T ( 1 ) T ( 5 ) T ( 1 ) T ( 5 ) E R E R E R E R E R E 1 ) T ( 5 ) T ( 1 ) T ( 5 ) T ( 1 ) T ( 5 ) P G P P G P P G P P G P P G P P G P 1 ) T ( 5 ) T ( 1 ) T ( 5 ) E O C E O C E O C E O C E O C () (s s s s s s ss s s s s s ss s s s s s0 6 0( 1 ) T ( 5 ) s令 ( ) ( ) ( )1 2 6 0G ( G G G ) 1 , 2 , 3 , 4 , 5i i ii 6 01 2 6 01 2 6 0T ( ) T ( ) T ( ) E R E R E ) T ( ) T ( ) P G P P G P P G P ) T ( ) T ( ) E O C E O C E O C ( ) 1 , 2 , 3 , 4 , 5 ( ) 1 , 2 , 3 , 4 , 5 ( ) 1 , 2 , 3 , 4 , 5i ii i ii i is s s s is s s s is s s s i取 E R G P P O C 1 , 2 , 3 , 4 , 5 X G G Gi i ii i i is s s * T T 1 *00T( 1 ) ( 1 )( 2 ) ( 2 ) 5 0 1 2 5 0( 5 ) ( 5 )1( ) ( ) ( ) ( )21()2 G R G Q G Q s c s X c X X X c令1 2 5()X X X X，则 ()J 可以看成是线性模型 ) 0 , D ( )y e e e X R + G Q G 的广义最小二乘估计目标函数的 此得到最优估计是 ( 1 )( 2 )110( 5 ) ( ( ) ) ( ) X R G Q G X X R G Q G c参考文献 1 S.