森林生态系统碳循环模型参数优化1)

陈晨　　　　　　　沃文伟　　　　　范文义

(东北林业大学，哈尔滨，150040)　　　　(黑龙江省水文水资源勘测局)　　　　(东北林业大学)



森林生态系统碳循环模型参数优化1)

陈晨沃文伟范文义

(东北林业大学，哈尔滨，150040)(黑龙江省水文水资源勘测局)(东北林业大学)

摘要以2011年帽儿山生态站通量塔观测的二氧化碳通量(总初级生产力和显热通量)数据为基础，使用集合卡尔曼滤波的顺序同化技术，对北方生态系统生产力(BEPS)模型的关键参数进行优化。结果表明：参数在季节尺度上变化显著，通常在展叶期迅速增大，夏季达到稳定，秋季落叶期降低。根据优化的参数，模型的总初级生产力、生态系统呼吸的模拟值精度显著提高，精度分别达到91%、96%，比优化之前的模拟值精度提高了8%和11%。说明集合卡尔曼滤波的参数优化可以明显改善模型模拟碳水通量的能力。

关键词BEPS模型；通量数据；集合卡尔曼滤波；数据同化；季节变化

分类号S715.3

Optimization of Ecosystem Carbon Cycle Model Parameters//

Chen Chen

(Northeast Forestry University, Harbin 150040, P. R. China); Wo Wenwei(Hydrology and Water Resources Survey Bureau of Heilongjiang Province); Fan Wenyi(Northeast Forestry University)//

Journal of Northeast Forestry University，2016,44(5)：15-19.

With CO2flux observation data (gross primary productivity (GPP) and sensible heat flux (LE)) from Maoershan Ecosystem Station flux tower in 2011, we used sequential data assimilation and ensemble Kalman filter technique to optimize some of the key parameters of the Boreal Ecosystem Productivity Simulator (BEPS) model by taking into account the errors in inputs, parameters, and observations. We optimized the parameters by data assimilation including maximum photosynthetic carboxylation rate (Vcmax), and the slope of stomatal conductance and net photosynthetic rate (M). Parameters were optimized in daily steps. The parameters varied significantly at seasonal scales, was with usually rapid increase in leaf expansion period, and in summer would reach a steady and decline in senescence of leaves. According to optimized parameters, model simulated value of GPP and RE flux were significantly increased. GPP and RE simulated accuracy reached 91% and 96%, and the precision values of the simulated values before optimization were increased by 8% and 11%, respectively. Ensemble Kalman filter parameter optimization could significantly improve model capacity in simulating the carbon and water fluxes.

KeywordsBEPS model； Flux data； Ensemble Kalman filter； Data assimilation； Seasonal variability

森林生态系统是陆地生态系统的主体，在维持全球碳平衡、减缓全球气候变化等方面起到重要的作用[1-2]，而生态系统和大气之间的碳、水和能量交换对地球气候影响显著[3]。关于森林生态系统的碳水能量循环及其驱动机制的研究已经成为当今全球气候变化研究的重点问题[4]。

陆地生态系统模型已经广泛应用于模拟生态系统生产力、水分消耗以及温室气体的排放。由于生态过程的复杂性及其时空变异的特征，模型参数值通常存在空间和时间的异质性，在应用到大区域的时候，这些参数通常是不能直接测量的(如冠层级别的参数)，但又显著地影响模型模拟碳、水和能量通量的能力。模型参数需要校准以确保模型模拟值和观测值一致，因此，模型参数的准确估计对于地表通量的模拟是十分重要的[5]。

数据同化技术是优化模型结构、状态和参数的新方法，已经在气象和海洋领域进行了广泛应用[6]。传统的模型参数估计方法通常假设参数是不随时间变化的，忽略了模型参数随时间变化的因素。数据同化技术结合了基于过程的模型、参数、观测值和状态变量的先验估计，对模型进行预测更新，并考虑参数误差和输入输出的关系以及模型结构误差所引起的不确定性。数据同化的主要目的是利用观测值寻找模型的状态变量的最佳估计。顺序同化技术在气候、水文和生态模型已经被广泛应用于递推估计的状态、时变参数和预测值的不确定性。集合卡尔曼滤波是传统的卡尔曼滤波的扩展，它基于蒙特卡洛方法和数据递归处理，使非线性动态系统具有高斯概率的最优估计，它可以递归地更新模型的参数和状态，提高对模型不确定性因素的处理。其基本思想是根据集合预报的结果来估计状态变量与观测变量之间的协方差，再利用观测资料和协方差来更新分析，得到分析集合，继续向前预报。

通量数据常用于校准和验证碳循环生态系统模型，连续的涡度相关通量数据，可以提供研究碳循环生态系统模型的一些关键生理生化参数的量值和季节性变化的条件。北方生态系统生产力(BEPS)模型经过多次的修改和更新，详细解释了光合作用、呼吸作用、能量分配、水文和土壤生物化学过程，以及模拟净生态系统生产力(NEP)。一些研究已经应用集合卡尔曼滤波对BEPS模型的参数进行了不同程度的优化，Mo et al[7]使用集合卡尔曼滤波和加拿大FRCN数据，对BEPS模型中的5个参数进行了优化，证实了通过涡度相关数据获得模型参数的季节变化和年际变化的可行性，提高了模型模拟通量的估计；Lin et al[8]在Mo et al的基础上基于两阶段数据同化方案TSDA，使用MODIS遥感数据获得土壤水分，对BEPS模型的参数和状态变量(土壤湿度)同时进行了同化；Ju et al[9]使用集合卡尔曼滤波和千烟洲涡度相关数据，对BEPS模型进行了参数优化，比较了参数与土壤含水率之间的关系，模拟了干旱对碳水通量的影响。

帽儿山地区是我国典型的天然次生林和人工林分布区，在东北林区具有代表性，相关学者已经利用碳循环模型对东北林区碳水通量模拟进行了研究，但模型参数大多是引用文献值或使用经验值，具有一定的主观性，缺少对优化模型参数的研究，这样难以深入地了解碳水循环过程，模拟结果也存在较大的不确定性。因此，本研究将帽儿山通量塔获取的通量数据使用集合卡尔曼滤波的数据同化方法同化到BEPS模型，目的是对BEPS模型中，与光合作用、蒸腾作用相关的两个重要参数最大羧化速率(Vcmax)和BWB方程斜率(M)进行了优化，从而获得这些参数的季节变化，提高模型的模拟精度。

1研究区概况

帽儿山森林生态站位于黑龙江省尚志市帽儿山东北林业大学实验林场，平均海拔400 m，平均坡度15°，土壤类型主要为暗棕壤。帽儿山地区属于长白山山脉的支脉，张广才岭西北部的余脉。年平均气温3.1 ℃，年降水量629 mm，年平均蒸散量864 mm，四季分明，夏季温暖湿润，冬季寒冷干燥，属于典型的大陆性季风气候，5—10月为生长季。

现有植被主要是阔叶红松林经过多次不同程度的人为干扰后形成的东北北部山区典型的天然次生林和人工林。主要树种有胡桃楸(Juglansmandshurica)、水曲柳(Fraxinusmandshurica)、五角槭(Acermono)、春榆(Ulmusjaponica)、白桦(Betulaplatyphylla)、山杨(Populousdavidiana)、大青杨(Populusussuriensis)、枫桦(Betulacostata)和紫椴(Tiliaamurensis)等，还包括少量红松(Pinuskoraiensis)和落叶松(Larixgmmelinii)等。

帽儿山通量观测塔建于2007年，塔高36 m。通量塔周围植被为天然阔叶落叶林，林龄约为60 a，林层高约18 m，林分复层结构，下层灌木发育良好。安装塔上的通量系统一直运转，记录每半小时生态系统与大气之间的碳、水和能量交换。装配在塔上的三维声波风速计和温度计(CSAT3，Campbell Scientific Inc.,USA)分别测定风速和温度的变化，利用闭路红外气体分析仪(LI-7000，Licor Inc，USA)和三维超声波风速仪分别测定CO2和水蒸气的变化[10-11]。

2研究方法

2.1数据收集和处理

本研究收集了2011年1月1日至2011年12月31日帽儿山森林生态站的气象数据和通量数据。入射太阳辐射、空气温度、相对湿度、风速和降水量的半小时气象观测数据等用于驱动BEPS模型；总初级生产力(GPP)、生态系统呼吸(RE)、潜热通量(LE)、显热通量(SH)等涡度相关通量测量值用于参数优化。对通量数据进行质量控制和选择，对于缺省数据采用间隙填充的方法进行插补。

帽儿山生态站的森林主要是天然阔叶落叶林，叶面积指数(LAI)的季节性变化相对较大。叶面积指数通过LAI-2000植物冠层分析仪和TRAC冠层分析仪测量，并用测量值经过3次样条插值法估算出LAI每日的值(见图1)。

图1　2011年帽儿山生态站叶面积指数的变化

2.2模型描述

本研究应用BEPS生态系统模型[12-15]模拟帽儿山生态站的森林生态系统碳水交换的过程，BEPS生态系统模型是一个基于过程的模型，它包括能量分配、光合作用、自养呼吸、土壤有机质分解、水文过程和土壤热交换模块。

光合速率的计算：

A=min(Ac，Aj)-Rd；

(1)

Ac=VcmaxfV(T)(Ci-Γ)/(Ci+Kc(1+Oi/Ko))；

(2)

Aj=JmaxfJ(T)(Ci-Γ)/4(Ci-2Γ)。

(3)

式中：A、Ac、Aj分别表示净光合速率、Rubisco酶总限制光合速率、限制总光合速率；Rd是白天叶片暗呼吸；Vcmax是在25 ℃的最大羧化速率；Jmax是在25 ℃的电子传递速率；Ci、Oi分别表示细胞间CO2和氧气浓度；Γ是不包含暗呼吸的CO2补偿点；Kc、Ko分别表示CO2和O2的米氏常数；fV(T)、fJ(T)分别表示空气温度(T)对Vcmax、Jmax的响应函数，Jmax=2.39×Vcmax-14.2，所以只需得到Vcmax的值，通过计算就可以得到Jmax的值。

气孔导度的计算：气孔导度gs是根据修改后的Ball-Woodrow-Berry经验模型计算得出，它与光合速率(A)线性相关。

gs=m(Ahr/Cs)+go。

(4)式中：gs为主气孔导度；go为残余导度；hr、Cs分别表示叶面相对湿度和CO2浓度；m是一个经验系数。

gs=fw(m(Ahr/Cs))+go。

(5)式中：fw是一个土壤水分胁迫因子，fw和m组合在一起重写为M，M为Ball-Woodrow-Berry方程的斜率。

BEPS模型运行：模型驱动数据包括气象变量(太阳辐射、长波辐射、温度、相对湿度、风速、降水量)和叶面积指数。模型输出数据包括GPP、NPP、NEP、自养呼吸和异养呼吸、净辐射、LE和SH、土壤湿度和温度。BEPS模型模拟的时间步长设置为30 min。本研究中使用模型输出的每日GPP作为观测值。

2.3参数优化方法

BEPS模型的参数优化是通过基于集合卡尔曼滤波的参数估计方法实现的。集合卡尔曼滤波是由Evensen[16]在1994年基于传统卡尔曼滤波开发的顺序数据同化方法。在线性情况下，卡尔曼滤波是最佳的递归处理算法[17]。在非线性情况下，扩展卡尔曼滤波将非线性线性化[18-19]。集合卡尔曼滤波作为扩展卡尔曼滤波的替代，它基于蒙特卡洛方法生成集合，并根据集合预报的结果来估计状态变量和观测变量之间的协方差。

集合卡尔曼滤波数据同化过程中有两组动态集合，分别是参数集合和通量预测集合。他们分别是通过对初始值添加服从高斯概率分布的扰动噪声得到一组集合。

(6)

式中：下标t表示时间，上标i+表示分析值、i-表示预测值。

参数成员集合更新为：

(7)

卡尔曼增益矩阵的计算：

(8)

通过对模型中的参数进行敏感性分析，选择最敏感的两个参数Vcmax和M进行优化，这两个参数对于模拟GPP有着显著的影响。Vcmax的初始值、标准差和范围分别为10、4和5～80，M的初始值、标准差和范围分别为4、0.4和4～14。这两个参数通过同化通量塔站点的总初级生产力(GPP)和潜热(LE)通量的观测值进行优化。参数优化是以每天的步长进行。通量数据的观测误差是与通量大小相关，LE和GPP的观测误差根据试验被设置为15%[20]，GPP和LE通量误差的相互相关被假定为可忽略不计。集合的大小应足够大以确保误差协方差的正确估计，而太大的集合会带来沉重的计算负担，经过试验集合大小设置为100。当测量的GPP和LE可用时，对参数进行优化并进行更新，更新的参数用来模拟接下来一天的通量。

3结果与分析

3.1优化参数的季节性变化

在集合卡尔曼滤波数据同化过程中，参数加入校正值后每日更新，然后将更新后的参数应用到下一天的模型计算中。两个参数与光合作用和呼吸作用相关，所以只在生长季进行优化，优化的参数在每天的步长都有一定的波动。

Vcmax表示在参考温度25 ℃时叶片光合作用的能力，它主要与叶片含氮量线性相关，其次也受光合作用酶、叶龄以及土壤水分胁迫等其他因素的影响。图2表现了Vcmax在2011年的季节性变化。2011年，该参数在10～65 μmol·m-2·s-1的范围变化，全年Vcmax的平均值约为45 μmol·m-2·s-1。光合作用随着生长季的变化而变化，Vcmax初始值设为10，参数在5—6月份迅速增大、7月中下旬逐渐达到峰值(65 μmol·m-2·s-1)，9月份叶片衰老，参数值降低。

M表示在ball-wood-berry模型中光合作用和气孔导度关系的斜率，它也表现出明显的季节性变化。图3显示了M在2011年的季节性变化。2011年，该参数在4～11的范围变化，全年M的平均值约为8。在2011年，M的初始值设为4，在早期生长季约为6，然后在6月底迅速增加到10，在7月略有下降且有很大的波动，在8月中旬达到峰值11，9月份又下降到8，这和9月份的干旱有关，最终在生长季结束前达到相对稳定的值为9。2011年的年降水量约为620 mm，是一个相对湿润的年份。

图2　参数Vcmax在2011年的季节性变化

图3　参数M在2011年的季节性变化

3.2优化后BEPS模型在通量模拟中的表现

模型使用优化的参数对通量进行模拟，通过比较通量的观测值和模拟值，评价模型参数优化的效果。由图4可知，参数优化后，模拟GPP达到观测GPP的91%，这一结果与使用校准的参数集(Vcmax=57.7，M=8)相比，模型模拟的效果有显著的改善，精度提高了8%。由图5可知，参数优化后，模拟RE达到观测RE的96%，与使用校准的参数集相比，模拟的精度提高了11%。

图4　2011年预测和观测每日GPP之间的比较

图5　2011年预测和观测每日RE之间的比较

4结论与讨论

结合涡度相关通量观测数据和BEPS生态系统模型过程模型，探讨了使用集合卡尔曼滤波数据同化方法优化模型关键参数。在25 ℃时，Vcmax、M这两个参数分别影响光合作用和蒸腾作用过程，提高了生态系统模型所使用参数的准确性，减少了直接使用模型参数模拟通量的误差。参数优化结果表明，这些参数表现出明显的季节性变化，光合作用呈现突出的季节性变化，通常在生长季节初期迅速增加并达到一个稳定的水平，然后在生长季节的末期突然下降。以往参数设置为稳定值，这会导致在不正常(如干旱等)的气候条件下，模型预测值的明显偏差。

通过集合卡尔曼滤波的参数优化，模型模拟的GPP、RE和水汽通量的估计值的精度在半小时和每日的时间尺度可以显著提高。使用集合卡尔曼滤波的数据同化可以获得时变参数，优化的参数可以显著改善生态模型的准确性。

但是需要明确以下几点：(1)本研究所使用的通量数据缺少小时观测GPP，所以使用每日GPP作为集合卡尔曼滤波中的观测值，并与模型输出日GPP进行比较，若使用小时观测GPP，参数优化结果的精度会有所提高。(2)本文优化的参数在实地测量中很难获得，Vcmax虽然可以使用LI-6400光合仪测定，但难度较大，如能获取一部分参数实测值与模拟的参数进行比较能提高参数优化的精度。(3)模型参数优化后，通量观测值和模拟值仍存在一定差异，这可能是由于模型结构方面的不确定性引起的，以后的研究重点应该放在模型机理结构上。

参考文献

[1]DIXON R K, SOLOMON A M, BROWN S, et al. Carbon pools and flux of global forest ecosystems[J]. Science,1994,263:185-190.

[2]YUDE P, BIRDSEY R A, JINGYUN F, et al. A large and persistent carbon sink in the world’s forests[J]. Science,2011,333:988-993.

[3]BONAN G B. Forests and Climate Change: Forcings, Feedbacks, and the Climate Benefits of Forests[J]. Science,2008,320:1444-1449.

[4]CAO M, WOODWARD F I. Dynamic responses of terrestrial ecosystem carbon cycling to global climate change[J]. Nature,1998,393:249-252.

[5]LUO YIQI, WENG ENSHENG, WU XIAOWEN, et al. Parameter identifiability, constraint, and equifinality in data assimilation with ecosystem models [J]. Ecological Applications, 2009, 19(3):571-574.

[6]MORADKHANI H, SOROOSHIAN S, GUPTA H V, et al. Dual state-parameter estimation of hydrological models using ensemble Kalman filter[J]. Advances in Water Resources,2005,28(2):135-147.

[7]MO X, CHEN J M, JU W, et al. Optimization of ecosystem model parameters through assimilating eddy covariance flux data with an ensemble Kalman filter[J]. Ecological Modelling,2008,217(1):157-173.

[8]JU W, WANG S, YU G, et al. Modeling the impact of drought on canopy carbon and water fluxes for a subtropical evergreen coniferous plantation in southern China through parameter optimization using an ensemble Kalman filter[J]. Biogeosciences,2010,7(3):845-857.

[9]ZHU L, CHEN J M, QIN Q, et al. Optimization of ecosystem model parameters using spatio-temporal soil moisture information[J]. Ecological Modelling,2009,220(18):2121-2136.

[10]王静,王兴昌,王传宽.基于不同浓度变量的温带落叶阔叶林CO2储存通量的误差分析[J]. 应用生态学报,2013,24(4):975-982.

[11]焦振,王传宽,王兴昌.温带落叶阔叶林冠层CO2浓度的时空变异[J].植物生态学报,2011,35(5):512-522.

[12]LIU J, CHEN J M, CIHLAR J, et al. A process-based boreal ecosystem productivity simulator using remote sensing inputs[J]. Remote Sensing of Environment,1997,62(97):158-175.

[13]CHEN J M, LIU J, CIHLAR J, et al. Daily canopy photosynthesis model through temporal and spatial scaling for remote sensing applications[J]. Ecological Modelling,1999,124(2/3):99-119.

[14]JU W, CHEN J M, BLACK T A, et al. Modelling multi-year coupled carbon and water fluxes in a boreal aspen forest[J]. Agricultural & Forest Meteorology,2006,140(1):136-151.

[15]CHEN B, CHEN J M, JU W. Remote sensing-based ecosystem-atmosphere simulation scheme (EASS)-Model formulation and test with multiple-year data[J]. Ecological Modelling,2007,209(2/3/4):277-300.

[16]EVENSEN G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[J]. Journal of Geophysical Research,1994,99(C5):10143-10162.

[17]KALMAN R. New approach to linear filtering and prediction problems[J]. Trans AMSE J Basin Eng，1960,82(D):35-45.

[18]JAZWINSKI A H. Stochastic Processes and Filtering Theory[M]. New York: Academic Press,1970.

[19]EVENSEN G. The Ensemble Kalman Filter: theoretical formulation and practical implementation[J]. Ocean Dynamics,2003,53(4):343-367.

[20]LASSLOP G, REICHSTEIN M, KATTGE J, et al. Influences of observation errors in eddy flux data on inverse model parameter estimation[J]. Biogeosciences Discussions,2008,5(1):751-785.

收稿日期：2016年1月17日。

第一作者简介：陈晨，男，1990年2月生，东北林业大学林学院，硕士研究生。E-mail：424520081@qq.com。通信作者：范文义，东北林业大学林学院，教授。E-mail：fanwy@163.com。

1)“十二五”国家科技支撑项目(2011BAD37B01)；国家自然科学基金项目(31300533)。

责任编辑：王广建。