森林作为陆地生态系统的重要组成部分，其净初级生产力(NPP)反映了森林固定大气CO₂的能力，是研究森林生态系统碳循环的关键，也是研究全球气候变化的基础。因此，准确估算森林生态系统NPP是很多学者关注的问题。目前，估算NPP的模型有多种，主要分为统计模型、光利用率模型和过程模型：统计模型以植被净第一生产力与气候因子之间的关系为基础，应用统计学方法建立回归模型(Lieth，1975)，较为典型的统计模型有迈阿密(Miami)模型(Lieth，1975)和筑后(Chikugo)模型(Uchijima et al.，1985)；光利用率模型是利用植被指数计算吸收光合有效辐射，进而计算森林NPP；过程模型研究森林植被的生理生态学特征及区域尺度上碳循环的过程，包括BIOME模型(Neilson et al.，1992)、BEPS模型(Liu et al.，1999)和TEM模型(Raich et al.，1991)等。这些模型在计算NPP时大部分只考虑了气候因素，忽略了外界干扰(火灾、病虫害、砍伐等)对NPP的影响。如陈国南(1987)利用迈阿密模型估算了我国植被的可能生产力； 何勇等(2007)利用植被与大气的相互作用模型AVIM及气象数据模拟了1971—2000年陆地生态系统NPP的变化特征； 赵国帅等(2011)利用GLOPEM-CEVSA模型估算了我国东北森林2000—2008年净初级生产力的时空分布格局； 王培娟等(2008)采用BEPS模型模拟了长白山森林植被NPP，并分析其主要影响因子及其敏感性。这些研究结果表明，东北森林NPP值较高，且温度与降水在一定范围内对森林NPP产生正反馈影响。但是这些研究估算的森林NPP精度只能达到70%左右，主要是由于模型使用的决定NPP大小的因子十分简单，导致模型忽略其他环境因子的作用以及异常条件对NPP的影响，也没有考虑到森林面积变化、森林类型变化对森林生态系统NPP的影响。

很多研究表明，干扰对森林NPP的影响远大于气象等其他环境因子，干扰后森林的再生长也是碳固定研究关注的问题。干扰的出现改变了森林的年龄结构，森林的NPP空间分布也发生了相应变化。因此，在长时间尺度上森林NPP的估算要考虑林分年龄的影响，这对进一步用碳循环模型准确估算森林生态系统碳源/汇分布至关重要。陆地生态系统碳循环模型InTEC模型考虑了气候及干扰对森林碳循环的双重影响，可模拟几百年来在火灾、病虫害及收获干扰下森林生态系统NPP的变化(Chen et al.，2003)。

本研究利用黑龙江省森林收获表、生物量相容性相对生长方程及叶片与细根周转率的数据，建立森林NPP与林分年龄之间的关系式，然后采用InTEC模型估算1901—2008年间小兴安岭森林NPP的时空变化，为进一步估算森林生态系统碳源/汇分布提供重要数据，也为森林管理与决策提供科学依据。

小兴安岭位于黑龙江省东北部(127°42′—130°14′ E，46°28′—49°21′N)，属低山丘陵地形，北部多台地、宽谷，中部为山势和缓的低山丘陵，南部山势较陡。该区属于温带大陆季风气候，冬季寒冷漫长，夏季温热短暂，年均气温-1～1 ℃，年降水量550～670 mm。区内森林茂密，树种较多，森林覆盖率72.6%，森林类型是以红松(Pinus koraiensis)为主的针阔混交林，主要树种为红松、兴安落叶松(Larix gmelinii)、樟子松(Pinus sylvestris var. mongolica)、山杨(Populus davidiana)、蒙古栎(Quercus mongolica)、紫椴(Tilia amurensis)、黑桦(Betula dahurica)和白桦(Betula platyphylla)等。主要土壤类型为暗棕壤，还有草甸土和沼泽土，部分高山地带为林地棕色针叶林土。

利用2005年小兴安岭86个森林清查固定样地的数据(树种组成、林分年龄、林分密度、立地指数、树高和胸径等)及收获表计算不同年龄森林的NPP。收获表由黑龙江省清查规划设计院2010年根据1990—2005年森林清查数据编制而成，包含估算林木蓄积与生物量的必需信息，如年龄、胸径、树高、活立木密度及蓄积增长量。根据黑龙江省主要树种非线性度量误差相容性生物量模型(董利虎等，2011)计算东北林区不同树种生物量，样地的林分生物量数据根据单木生物量与样地林木株树乘积获得。然后，采用林木生长方程模拟不同年龄下林木生物量，为计算林分NPP提供每年生物量增长量数据。对于混交林，要根据样地树种组成加权平均计算生物量及NPP。

InTEC模型估算森林生态系统NPP所需的气象数据包括月平均气温、月总降水量、月平均太阳辐射和年大气CO₂浓度。1901—2008年月平均气温及降水量采用UK气候研究中心全球0.5度分辨率数据集(CRU3.0，UK Climate Research Unit 3.0)。1948—2008年月平均太阳辐射数据采用NCEP/NCAR(national center for atmosphere research/national centers for environmental prediction)的T62高斯再分析数据。1901—1948年月平均太阳辐射数据采用Thornton等(1999)的估算方法，并根据1950—1959年10年的月平均太阳辐射校正1948年前计算的月平均太阳辐射。1958—2009年大气CO₂浓度来自于美国夏威夷的冒纳罗亚观测站(20° N，156° W)观测数据集。1901—1957年大气CO₂浓度数据根据CGCM2天气模型(coupled global changed model 2)的模拟数据插值计算得到，前提是假设研究区域大气CO₂浓度没有空间变化。所有的数据均采用UTM投影，WGS84坐标系，数据分辨率为1 000 m。

植被参数数据也是模型估算NPP的重要输入参数，包括叶面积指数LAI、植被类型、参考年NPP和林分年龄。小兴安岭2003年的叶面积指数数据采用8天250 m分辨率无云MODIS影像上定量估算的结果。选取每月平均10天叶面积指数中最大值为全年的LAI最大值，具体的MODIS影像处理方法及叶面积指数的定量反演方法参见文献Deng等(2006)。森林植被类型数据集采用2006年东北地理与农业生态研究所出版的1∶1000 000东北植被分布图。东北林区8大主要树种包括红松、兴安落叶松、樟子松、白桦、紫椴、蒙古栎、山杨和黑桦，主要森林类型有针叶混交林、阔叶混交林和针阔混交林。对小兴安岭二类调查小班数据进行泰森多边形插值获得林分年龄数据。本研究采用2003和2005年的MODIS 陆地4级标准NPP产品数据作为参考与验证的标准。

土壤田间持水力和萎蔫点数据采用国际地圈生物圈计划发布的全球选择性土壤特性格网数据(international geosphere-biosphere program，global gridded surfaces of selected soil characteristics)。土壤深度数据来自于美国田纳西州橡树岭国家实验室分布式数据中心(Oak Ridge National Laboratory Distributed Active Archive Center，Tennesse，USA)提供的全球土壤质地和持水力数据集。土壤沙土、黏土及淤泥含量数据来源于联合国粮农组织(FAO)和维也纳国际应用系统研究所(IIASA)构建的世界和谐土壤数据库(harmonized world soil database，HWSD)。

林分NPP为生物量(W)的增长量、凋落物量、细根更新的量及树木死亡量之和。树木在达到过熟林之前每年的死亡量很小，这里忽略不计； 树木达到过熟林之后死亡量增加，在InTEC模型中按照发生干扰后死亡进行处理。得到NPP的计算公式如下:

式中：dW/dt为在特定的年龄(t)下每年树木生物量中干物质总量的变化量，根据收获表计算获得； L_f和L_fr为每年叶片和细根的更新量。

生物量的计算分为4部分: 树干生物量W_s、枝生物量W_b、叶生物量W_f和粗根生物量W_r。根据董利虎等(2011)建立的黑龙江省主要树种非线性度量误差相容性生物量模型，计算东北林区不同树种生物量。总生物量及各部分生物量计算公式如下:

总生物量 W_t:

地上生物量W_a:

根生物量W_r:

茎生物量 W_s:

冠生物量 W_c:

枝生物量 W_b:

叶生物量W_f:

式中: a₁，b₁，c₁，r₁，r₂，r₃，r₄，r₅，r₆，r₇，r₈为利用实测数据拟合相应生物量方程得到的模型系数值；D为胸径；W为树高。

这一方法解决了4部分生物量分别计算加和不等于总生物量的问题。样地的生物量数据根据单木生物量与样地林木株树乘积获得。东北林区8个树种的生物量模型及参数见表 1。为了避免样地调查数据林分年龄范围的限制导致的叶生物量模型估算结果没有上限，这里假设叶生物量的最大值出现在生物量增长速度最大时。

表 1 生物量相容性模型及相应参数 Tab.1 Compatible model and its corresponding coefficients for biomass estimations

表 1 生物量相容性模型及相应参数 Tab.1 Compatible model and its corresponding coefficients for biomass estimations 树种Tree species ln a 1 b 1 c 1 r 1 r 2 r 3 r 4 r 5 r 6 r 7 r 8 樟子松 Pinus sylvestris var. mongolica －3.12 3.23 0.51 0.34 －0.49 0.35 0.18 1.22 －1.30 2.39 －0.40 兴安落叶松 Larix gmelinii －2.22 2.00 0.51 0.31 0.79 －0.86 5.6 0.09 －1.32 2.78 －0.77 红松 Pinus koraiensis －2.41 2.47 － 0.43 0.63 －0.94 0.13 1.10 －0.85 7.32 －0.86 黑桦 Betula dahurica －2.73 2.29 0.37 1.47 0.04 －0.74 0.06 1.88 －1.53 0.68 －0.35 紫椴 Tilia amurensis －3.17 1.88 0.91 1.88 －0.90 0.26 0.52 0.41 －0.86 5.60 －1.08 白桦 Betula platyphylla －2.89 2.08 0.69 1.05 0.48 －0.94 0.31 0.94 －1.01 1.24 －0.61 山杨 Populus davidiana －4.09 2.36 0.72 3.74 0.11 －1.14 0.01 1.59 －0.66 10.54 －1.26 蒙古栎 Quercus mongolica －2.92 2.21 0.63 1.62 0.51 －1.27 0.11 1.41 －1.22 1.70 －0.75

林木生长方程用来模拟生物量随年龄的变化情况。林木生长方程包括逻辑斯蒂方程、单分子式、坎派兹方程、考尔夫方程及理查德方程。这些生长方程均呈现出“S”形曲线，先慢后快，再慢最后达到平稳。这里采用理查德方程，如下:

式中: W_t为通过生物量方程估算的生物量，随着年龄(t)变化而变化； A，r和c分别为最大生物量、增长速率和林木同化能力，是生物量与年龄模拟关系的系数。

每年凋落物量等于叶生物量乘以叶周转率公式如下：

式中: T_f为叶周转率； C_f为叶片碳含量与生物量的转换系数，根据以往的研究结果(于颖等，2012)，这里取值0.44。对于每年凋落的阔叶林来说，叶周转率等于1，对于针叶林，其生长周期受温度、光照、营养、病虫害、大气污染等多因素影响，很难通过测量获得。本文采用White等(2000)的不同树种叶周转率数据。

细根更新量依赖于叶的更新量及NPP分配到细根与叶的比例，公式如下：

式中: e为NPP分配到细根与叶的比例系数(后文表述为e值)。

叶更新量与细根更新量二者关系的原理见文献Thornton(1998)。不同树种的e值与叶周转率从文献White等(2000)中获得，见表 2。

表 2 不同树种叶周转率和e值 Tab.2 Leaf turnover rates and e values for different species

表 2 不同树种叶周转率和e值 Tab.2 Leaf turnover rates and e values for different species 树种Tree species 叶周转率 Leaf turnover rate e 樟子松 Pinus sylvestris var mongolica 0.385 1.03 兴安落叶松 Larix gmelinii 1.000 1.4 红松 Pinus koraiensis 0.400 1.40 黑桦 Betula dahurica 1.000 1.26 紫椴 Tilia amurensis 1.000 1.20 白桦 Betula platyphylla 1.000 1.26 山杨 Populus davidiana 1.000 1.20 蒙古栎 Quercus mongolica 1.000 1.20

InTEC模型估算长时间尺度森林NPP时首先设定初始年NPP，即工业革命前期森林NPP，然后根据气候及林分年龄对NPP的影响推算出100多年后森林NPP的变化值，用参考年2003年森林NPP作为校正值，反复迭代调整初始年NPP，直到估算的NPP值与参考年NPP值相等或相差1%以内迭代停止，此时得到准确的工业前期森林NPP，再根据NPP与年龄关系计算历年来森林NPP变化情况。

模拟根据2.2方法计算的NPP与林分年龄(t)的关系，采用下式(Chen et al.，2003)，其中M，b，d，g是模型的系数。将不同年龄下的NPP除以NPP的最大值得到归一化NPP与年龄的关系，用于计算不同树种历年NPP(t)值：

小兴安岭林区各树种组成林分的NPP及其各个组分碳贮量随年龄的变化见图 1。各树种的总NPP均随着年龄的增长迅速增加到最大值，即达到林木的成熟年龄，然后随着年龄的继续增加逐渐减少，直至达到稳定值。幼龄林和中龄林，碳主要固定在树干、枝及粗根上；成熟林后，林木的NPP开始迅速减小；老龄林，林木的NPP主要取决于叶和细根的周转更新，在林分年龄较大时，叶更新在达到最大值后不再变化，NPP基本稳定不变。以往的研究(He et al.，2012)发现，在老龄林中，森林NPP随着年龄的增加有较小的下降，这是老龄林叶子更新能力下降所引起的。

	图 1 NPP及其各个组分碳贮量随林分年龄的变化 Fig. 1 NPP and its components carbon storage change with stand age

不同林分NPP最大值出现的年龄因林分类型和树种组成的不同而不同(图 1)。对于纯林，人工林固碳能力高于天然林。在纯林中，针叶林的生长速度最慢，NPP达到最大值时晚于阔叶林，不同的针叶纯林一般在20～30年NPP达到最大，而阔叶林10年左右NPP就达到了最大值。对于混交林，针阔混交林的固碳能力大于针叶混交林及阔叶混交林，最大值为256 gC·m^-2a^-1。阔叶混交林生长速度最快，树龄17 年时NPP达到最大值，其次是针阔混交林，树龄23 年时NPP达到最大。针叶混交林生长速度最慢，43年树龄时NPP达到最大。

本研究统计了归一化NPP与年龄的关系(图 2)。天然针叶林的固碳能力大于天然阔叶林与针阔混交林。单一树种针叶林纯林在NPP达到最大之后随年龄的增加迅速减小，这一趋势比单一树种阔叶林纯林明显。这与He等(2012)的研究结果相同。生产力较大的树种，NPP随年龄的变化也较大，呈现出年轻时快速增长。年老时快速下降的状态。不同林分类型NPP与年龄关系方程的拟合系数与模型决定系数R²见表 3。

	图 2 不同林分类型归一化NPP与年龄关系 Fig. 2 Relationships between normalized NPP and age for different stands

表 3 不同林分类型NPP与年龄关系方程的拟合系数与模型决定系数R² Tab.3 Coefficients and R² of age-related NPP functions for different species for eq. (12)

表 3 不同林分类型NPP与年龄关系方程的拟合系数与模型决定系数R2 Tab.3 Coefficients and R2 of age-related NPP functions for different species for eq. (12) 林分类型Stand type M b g d R 2 樟子松林 Pinus sylvestris var. mongolica stand 172.477 0.628 6.852 3.437 0.996 1 针叶林Coniferous forest 166.656 1.712 35.041 0.495 0.993 8 针阔混交林Coniferous and broad-leaved mixed forest 144.989 2.888 20.944 0.465 0.998 0 蒙古栎林 Quercus mongolica stand 48.000 5.808 55.391 0.122 0.999 8 人工杨树林Planted poplar stand 120.076 4.789 9.420 0.458 0.999 8 山杨林 Populus davidiana stand 96.369 3.176 32.026 0.292 0.999 7 阔叶林Broad-leaved forest 84.991 4.048 20.142 0.388 0.997 5 兴安落叶松林 Larix gmelinii stand 397.404 0.803 7.186 2.513 0.996 9 红松林 Pinus koraiensis stand 181.997 0.053 3.147 5.228 0.998 2 黑桦林 Betula dahurica stand 32.278 4.720 43.153 0.147 0.999 9 紫椴林 Tilia amurensis stand 143.107 2.394 38.398 0.122 0.999 7 白桦林 Betula platyphylla stand 150.000 1.892 14.113 1.050 0.997 8

综上所述，从森林生态系统固碳能力的角度分析，小兴安岭人工林的固碳能力要大于天然林，因为天然林在年龄较大时成为中性碳库，而人工林仍然是个碳汇。

根据2005年小兴安岭86个森林清查固定样地数据(树种组成、林分年龄、林木密度、立地指数、树高、胸径等)以及式(1)计算2005年样地NPP，与InTEC模型模拟的2005年NPP值及MODIS的NPP产品进行比较，结果见图 3和4。InTEC模型模拟的森林NPP与实测数据有较好的一致性，精度为85.1%。MODIS的NPP产品比实测NPP数据偏高，主要是因为MODIS的NPP产品采用的是光利用率模型，没有区分阴阳叶分别计算光合作用，高估了阴叶的有效光合利用率，使得计算的森林NPP结果偏大。

	图 3 InTEC模型模拟的NPP与实测值比较 Fig. 3 NPP comparisons of InTEC modeling and measurements

	图 4 MODIS的NPP产品与实测值比较 Fig. 4 Comparisons of MODIS NPP and measurements

对比分析InTEC模型模拟的小兴安岭森林NPP值与2005年对应区域MODIS的NPP数据产品，结果见图 5和表 4。本研究采用InTEC模型估算的2005年小兴安岭地区NPP值为96～642 gC·m^-2a^-1，平均值为313 gC·m^-2a^-1。2005年小兴安岭地区MODIS 数据估算的NPP值为156～696 gC·m^-2 a^-1，平均值为368 gC·m^-2 a^-1，从NPP值分布范围比较来看，MODIS数据产品高估了NPP。小兴安岭中南部NPP较高，北部及最南部区域NPP值较低，因为中南部海拔为300～600 m，低山丘陵局地生态环境有利于植被的生长，是森林分布面积较大的区域，以红松、樟子松等针叶林、针阔混交林为主，生产力较高。而北部地势平缓，最南部地势陡峭，均不利于植被的生长，森林面积分布较小，所以NPP值较低。小兴安岭与大兴安岭接壤处NPP值较中部低，主要是气温与中部相比偏低，且树种分布以白桦、蒙古栎为主，生产力较低。

	图 5 InTEC模型估算的NPP与MODIS的NPP Fig. 5 NPP estimated from InTEC model and MODIS NPP product

表 4 InTEC模拟NPP与MODIS的NPP统计值 Tab.4 NPP statistics based on InTEC and MODIS

表 4 InTEC模拟NPP与MODIS的NPP统计值 Tab.4 NPP statistics based on InTEC and MODIS gC·m-2 a-1 项目Item 最大值Max. 最小值Min. 平均值Mean InTEC NPP 642 96 313 MODIS NPP 696 156 368

小兴安岭森林NPP从工业前期的239 gC·m^-2 a^-1增长到2008年的314 gC·m^-2 a^-1(图 6)，增加了30%。1900—1950年期间NPP逐年增加，到1960年NPP开始下降，主要是由于19世纪60年代外界的干扰持续增加，包括砍伐、病虫害和火灾的影响，使得NPP减少，这种趋势持续到1980年。1980—1990年小兴安岭森林NPP又开始呈现增加趋势，是由于气候变化、CO₂和氮沉降对NPP增加的综合影响弥补了外界干扰引起的NPP消耗。1990年 NPP开始下降，是由于近些年火灾病虫害的发生使得森林生长受到影响，干扰后森林年龄变小，多以幼龄林为主，NPP减少。

	图 6 InTEC模型估算小兴安岭森林NPP变化 Fig. 6 Change of forest NPP in Xiaoxing’an Mountains based on InTEC model

本研究发现小兴安岭幼龄林NPP随着年龄的变化迅速增加，在达到最大值后开始下降，最后达到稳定。基于InTEC模型估算的森林NPP与样地数据计算的森林NPP有较好的一致性，而MODIS的NPP产品却比模型估算的NPP偏高，这进一步说明MODIS的NPP产品需要验证调整后方可使用。

InTEC模型模拟结果表明，小兴安岭中南部NPP较高，北部及最南部区域NPP值较低，说明NPP受海拔及森林类型的影响，在海拔适中的地区NPP较高，对于平坦及陡峭山区，NPP值偏低。1901—2008年间小兴安岭森林NPP总体增加了30%，主要受气候变暖、大气N沉降、大气CO₂浓度增加的影响，在1920，1960和1990年由于外界的干扰使得森林NPP有下降趋势，但是干扰后森林恢复生长，使得NPP再次增加，所以适当的干扰也会促进森林NPP增加。

研究采用的总生物量模型估算精度较高，达到90%，茎生物量估算精度为87.5%。生物量模型对粗根和叶的估算精度相对较低，但是也超过了80%。本研究估算的NPP较其他研究低，主要是由于没有考虑下层林木NPP。

除了生物量模型的误差外，还有其他误差来源: 1)平均胸径和树高的使用平滑了生物量随年龄的变化； 2)从文献中查找的不同树种叶与细根的周转率、e值并不能完全代表本研究的树种； 3)假设叶子更新达到最大值后保持不变，但是有时在老龄林中叶子更新量会出现降低的趋势。由以上因素引起的NPP计算的不确定性误差在25%左右，并且随着年龄的变化有增加趋势。

本研究模拟了东北森林主要树种NPP与年龄的关系并分析其时空变化情况，为碳循环模型在东北森林中应用提供了重要的基础数据，也为森林决策者经营管理东北森林提供了重要参考依据。另外，NPP随年龄变化的同时也受森林立地环境的影响，不同的立地环境下NPP与年龄的关系不同，这是需要进一步研究的方向。