毛竹(Phyllostachys edulis)是我国南方最重要、最典型的竹林资源类型, 占全国竹林面积的70%左右(漆良华等，2009)。浙江省是我国毛竹分布最丰富的省份之一，2011年全省毛竹林面积为73.04万hm²，毛竹总株数211 282万株。据研究，毛竹林乔木层年固碳量为5.097 t ·hm^-2(周国模等，2004)，是速生阶段杉木(Cunninghamia lanceolata)林的1.46倍(方晰等，2002)、热带山地雨林的1.33倍(李意德等, 1998)、苏南27年生杉木林的2.16倍(阮宏华等, 1997)，故毛竹林具有很强的碳汇功能。毛竹林具有独特的林分特征与经营措施：1)所有毛竹林分都是异龄、速生的；2)采伐和挖笋是经营毛竹林的必要手段。由于毛竹林分特有的生物学特性与经营措施，使毛竹林的固碳量达到一定程度就不再增加，而其他类型森林的固碳量一般随着林分平均年龄的变大而增多。目前对森林固碳潜力的研究大多是基于国家尺度对全国森林碳储量进行估算(Zhang et al., 2003；吴庆标等，2008), 或是对单个省的森林固碳潜力进行研究(张骏等，2010)，对林分尺度森林固碳潜力及具有巨大碳汇功能、特有经营措施毛竹林的固碳潜力的报道都比较少。因此在全球气候变化与全球碳循环受到国内外学者广泛关注的大背景下，研究毛竹林固碳潜力有着非常重要的意义。

合理的毛竹林分结构是充分发挥毛竹林碳汇功能的基础(Law et al., 2003；张建国等, 2004)。研究表明毛竹林生态系统在未遭受严重干扰的情况下, 都具有特定的林分结构规律(孟宪宇, 1996)。但近年来, 由于片面追求经济效益, 许多毛竹林由天然粗放经营转变为集约经营(徐秋芳等, 2003), 集约经营毛竹林在获得更高经济收益的同时, 破坏了毛竹林分固有的结构规律，从而严重影响毛竹林的可持续经营(陈双林等, 2004)与毛竹林生态效益的充分发挥，研究毛竹林分结构规律对提高毛竹林碳汇功能及毛竹林的可持续经营具有非常重要的理论与实践意义。毛竹林分结构的调控是提高毛竹林碳汇功能必要手段。到目前为止，已有林分结构调整的研究其目的大多是合理的采伐与经营森林(陈双林等，2003；Venema et al., 2005；Luis et al., 2003)，而有关以生态效益与固碳潜力为目的的毛竹林分结构调整研究少见报道。本研究以林分非空间结构为基础，对毛竹林固碳潜力进行估测与调控，其目的是：1)揭示毛竹林碳储量达到最大时所具有的非空间结构；2)以现有毛竹林为基础，提出一种以固碳潜力为目标的毛竹林经营模式，为充分发挥毛竹林的固碳功能提供技术支持。

1 基于非空间结构的毛竹林分固碳潜力模型构建 1.1 毛竹林分特征分析

在我国，毛竹是一种特殊的地域性植物，与其他植物不同的是：1)毛竹的胸径生长与高生长在第1年内完成；2)毛竹林多有隔年出笋成竹即大小年习性, 因此, 毛竹林林分年龄结构按2年划分1个龄级；3)对单株毛竹而言，当胸径生长与高生长完成后，胸径与竹高不再发生变化，但其生物量还在变化，这种变化主要是由年龄引起的；4)毛竹林是一种异龄林，故除胸径外，年龄也是其林分结构的重要特征；5)毛竹林分密度是影响毛竹林生长的重要因子之一, 在植株各年龄、各胸径阶段都制约着林分地上和地下部分的生长, 造成植株间对资源的竞争或相互干扰(Turkington et al., 2005), 且林分密度控制又是营林措施的一个有效手段, 因此林分密度的调整对毛竹林的固碳量影响很大。基于以上对毛竹林分特征的分析，本研究选择林分密度、林分年龄和林分胸径作为影响毛竹林碳储量的基础指标。在本研究中，用林分平均胸径表示毛竹林分胸径，用立竹度表示林分密度，考虑到平均年龄不能准确描述毛竹林分年龄，故用林分中各龄级所占比例来表示毛竹林分年龄。

1.2 毛竹林(样地)生物量模型

单株毛竹生物量采用公式(1)计算，该模型是根据常用生物量模型与毛竹的生长规律，结合浙江省毛竹生物量调查数据建立的, 模型的相关系数R²=0.937, 在0.05置信水平下的预估精度为96.43%, 总系统误差为-0.021%, 符合生物量估算精度要求(周国模，2006)。

(1)

式中：M为单株毛竹生物量(kg)，A为单株毛竹年龄(度)，D为单株毛竹胸径(cm)。在计算单株毛竹生物量的基础上，通过累加样地中所有毛竹的生物量即可得毛竹林分生物量。

毛竹林分平均胸径D_g可用下式计算：

(2)

式中：N为林分立竹度，d_i为第i株毛竹的胸径。

令λ为4度及以上竹的度数，a₀, a₁, a₂, a₃, a₄, a₅为参数，x₁, x₂, x₃, x₄分别为1度竹、2度竹、3度竹和≥4度竹所占毛竹样地总株数的百分比，则由(1)式推导得毛竹林分(样地)生物量M_Total模型为：

(3)

对(3)式进行改进得(4)式：

(4)

式中：c₁, c₂, c₃, c₄均为≥0的参数，c₅, c₆均为>1的参数。公式(4)的生物学意义在于：1)c₁x₁+c₂x₂+c₃x₃+c₄x₄，D_g和N是影响毛竹林分生物量的3个基本因素，分别表示毛竹林分年龄、林分平均胸径和林分立竹度，且公式(4)所示结构与常用生物量模型(曾慧卿等，2007；唐守正等，2006)也相匹配; 2)在毛竹林中存在挖笋、采伐、自然死亡等现象，这样就会使毛竹林分生物量、林分年龄、林分平均胸径D_g和立竹度N发生改变，而公式(4)所示结构能反映毛竹林分的相互变化关系; 3)c₁x₁+c₂x₂+c₃x₃+c₄x₄反映出毛竹林是异龄林; 4)当c₁, c₂, c₃, c₄中的1个或2个为零时表明毛竹林中没有零参数所对应度数的毛竹; 5)毛竹林具有碳积累速度与生长速度快的特点(刘恩斌等，2009)，而D_g与N的幂c₅, c₆能反映毛竹林的这一特点; 6)由于D_g与N对毛竹林分生物量的影响比x₁, x₂, x₃, x₄大，所以模型中的c₅和c₆大于1而x₁, x₂, x₃, x₄的幂为1。可见公式(4)具有明确的生物学意义，经检验公式(4)的拟合精度(R²=0.994)高于公式(3)(R²=0.988)。

毛竹林分生物量M_Total乘以毛竹碳转化系数0.504 2(周国模等，2004)后，就可得毛竹林固碳量。

1.3 毛竹林固碳潜力模型

要构建毛竹林固碳潜力模型, 首先要建立毛竹林分立竹度、胸径的约束条件, 即建立株数、胸径之间的对应关系。据研究，异龄林株数按径级的分布可用负指数分布表示(Meyer, 1952；汤孟平，2004)，毛竹林分立竹度与胸径的约束条件为：

(5)

式中：k, a为参数, 为了使公式(5)适应范围更广，现对模型加以改进，改进后的模型为：

(6)

式中：b为参数。

其次要建立年龄约束条件。毛竹林分的年龄具有如下特点：1)毛竹林是异龄林；2)林分中低龄级毛竹的株数比高龄级的要多；3)从生物量的角度来考虑，对某一毛竹林而言，林分的年龄越大其吸收固定CO₂的量就越多。综合考虑这3方面的因素，即可建立毛竹林分年龄的约束条件，即公式(8)的约束条件。

结合以上的分析，可得基于非空间结构的毛竹林分最大生物量maxM_Total模型：

(7)

毛竹林分最大生物量maxM_Total乘以毛竹碳转化系数0.504 2(周国模等，2004)后，就可得毛竹林分固碳潜力。毛竹林分最大生物量maxM_Total可分为2部分。

第1部分为：

(8)

第2部分为：

(9)

采用Matlab提供的Linprog函数可对第1部分进行优化求解，具体优化过程按Linprog函数语法规则进行。根据函数极值求解原理，应用Matlab提供的Solve函数可对第2部分进行优化求解，具体优化步骤如下：1)把N=ke^-aD_g+b代入公式(9)中; 2) ; 3)用solve函数求方程(akD_g-k)exp(-aD_g)-b=0的根; 4)把求得的根代入N=ke^-aD_g+b，可得N的值; 5)由可得M₂的最大值。

2 模型应用 2.1 研究区概况

浙江省位于我国东南沿海(118°01′—123°10′E，27°06′—31°11′N)，气候温暖，降水充沛，地势呈自西南向东北倾斜趋势，西南部分为平均海拔800 m的山区，中部为丘陵、盆地，海拔100~500 m，土壤以黄壤和红壤为主，植被类型主要有针叶林、针阔叶混交林、常绿阔叶林、落叶阔叶林和竹林。浙江省是我国毛竹分布最丰富的省份之一，2011年全省毛竹林立竹量2 715株·hm^-2，当年生新竹占毛竹总株数的19.2%。

2.2 数据来源 2.2.1 全省毛竹连续清查样地

浙江省于1979年建立了森林资源连续清查体系，以5年为1个复查周期。共设置固定样地4 250个，样点格网为4 km×6 km，样地形状为正方形，边长28.28 m，本研究利用2004年的调查数据，选择基本为毛竹纯林的样地245块，每块样地毛竹胸径5~15 cm，毛竹株数18~416株。样地内5 cm以上的竹子均要调查记载，调查内容主要有胸径和年龄(当年生竹记为1度竹；2~3年生竹记为2度竹，4~5年生竹记为3度竹，依此类推)。

2.2.2 典型毛竹调查样地

由于本研究的目标是毛竹的固碳潜力，考虑到全省毛竹连续清查样地中包含的毛竹林分结构不够全面，故在临安市与安吉县不同立地条件的毛竹林中选择了与连续清查样地无重复的150块典型调查样地，其大小、调查方法与毛竹连续清查样地全完一致。

2.3 浙江省毛竹林分结构现状分析

将245块毛竹样地的胸径、年龄数据合并，并计算浙江省某径阶、某龄级的毛竹概率，再应用二元最大熵函数测量全省毛竹胸径、年龄联合分布信息(刘恩斌等，2010)，结果表明测量精度很高(回归离差平方和为9.976 8e-05，R²=0.996 0)，然后用Matlab软件绘制省域尺度毛竹胸径、年龄二元概率密度图(图 1)。从图 1可以看出: 1)1度毛竹在9径阶附近株数最多，2度毛竹在8径阶附近株数最多，3度毛竹在8径阶附近株数最多，≥4度毛竹在9径阶附近株数最多; 2) 1~2度、胸径8~9 cm的毛竹株数最多，随着年龄、胸径的增大毛竹株数在逐渐的减少; 3)对于同一径阶的毛竹，其株数的变化趋势是先增后减。

再对上述395块毛竹样地实测数据进行统计，得浙江省毛竹林分(样地)平均胸径为8.671 3 cm、林分平均株数为2 425株·hm^-2, 林分(样地)平均碳储量为13 002.379 7(kg·hm^-2)。

2.4 浙江省毛竹林固碳潜力模型参数拟合

把200块全省毛竹连续清查样地与100块毛竹典型调查样地数据作为建模数据，应用SPSS软件的非线性回归功能对模型(4)进行拟合，其参数最后5次迭代过程见表 1。从表 1可以看出，模型(4)的拟合参数值为：(c₁, c₂, c₃, c₄, c₅, c₆)=(0.240 4, 0.252 5, 0.240 5, 0.251 7, 1.704 4, 1.016 5)，此时R²=0.994, 再把建模数据以外的连续清查样地与典型调查样地作为检验样本，对模型(4)进行检验，此时R²=0.981，说明毛竹林分生物量模型的拟合精度与检验精度都非常高。

用395块毛竹样地与SPSS软件的非线性回归功能，对模型(6)进行拟合，其参数值为：(k, a, b)=(1 173.839 1, 0.039 5, -376.953)，此时R²=0.958 7。

2.5 浙江省毛竹分固碳潜力

在模型参数拟合结果的基础上, 由Matlab提供的Linprog函数对第1部分进行求解，其结果为：当x₁，x₂，x₃和x₄均为0.25时，目标函数值M₁达到最大。应用Matlab提供的Solve函数对第2部分进行求解，其结果为：当N为4 363株·hm^-2，D_g为12.169 1 cm时，目标函数值M₂达到最大。因此当x₁，x₂，x₃和x₄均为0.25，N为4 363株·hm^-2及D_g为12.169 1 cm时，M_Total达到最大，maxM_Total为83 737.038 7 kg·hm^-2，转换为碳储量为42 220.214 9 kg·hm^-2。

根据模型参数拟合结果与模型第1部分求解结果，用Matlab软件作函数M_Total=(c₁x₁+c₂x₂+c₃x₃+c₄x₄)D_g^c₅N^c₆的图像(图 2)。从图 2可以看出: 1)毛竹林分碳储量随着林分平均胸径的增大与样地株数的增加迅速增加; 2)林分平均胸径对毛竹林分碳储量的影响比样地株数要小; 3)林分平均胸径大的曲线斜率比林分平均胸径小的曲线斜率要大，株数大的曲线斜率比株数小的曲线斜率要大。但由于各种条件的限制，在实际中不存在图 2所示毛竹林分碳储量的最大值。

2004年毛竹林分(样地)平均碳储量为13 002.379 7 kg·hm^-2，而毛竹林分(样地)的最大碳储量可达42 220.214 9 kg·hm^-2, 毛竹林2004年的固碳量不足潜在固碳量的1/3，可见浙江毛竹的固碳潜力巨大。

2.6 浙江省毛竹林固碳潜力调控措施

本研究表明，当x₁，x₂，x₃和x₄均为0.25，N为4 363株·hm^-2，D_g为12.169 1 cm时，毛竹林分的固碳量达到最大，结合2004年毛竹林分结构现状，要使毛竹林分固碳量达到最大，应采取如下调控措施：1)不挖笋不采伐、增加毛竹株数是毛竹林初始的经营模式，但当样地中毛竹株数超过4 363株·hm^-2时，应采用挖笋与择伐等手段减少毛竹株数，以增大毛竹林分平均胸径；2)从c₅, c₆参数值及图 2看出，林分株数对毛竹林分固碳的影响最大，林分平均胸径次之，从c₀, c₁, c₂, c₃参数值及第1部分的优化结果可以看出，各度竹所占比例对毛竹林分固碳的影响最小，因此要提高毛竹林固碳量，首先要保证林分中有一定数量的毛竹，其次要保留林分中龄级较小且胸径较大的毛竹；3)在实际毛竹林分中，1度竹的株数最多，从各度竹株数所占比例相等可知，挖笋是减少毛竹林分株数的最主要手段，且地径较小的笋是选择的主要对象；4)龄级较大、胸径较小的毛竹是择伐的主要对象。

3 结论与讨论

综合以上分析，可以得出如下结论：1)本文构建了具有生物学意义的毛竹林分生物量模型，从拟合结果与检验结果可以看出本研究构建的模型精度很高，在此基础上，结合毛竹林分特征，构建了毛竹林分固碳潜力模型；2)当x₁，x₂，x₃和x₄均为0.25，N为4 363株·hm^-2，D_g为12.169 1 cm时，毛竹林分的固碳量达到最大，其值为42 220.214 9 kg·hm^-2。

但由于所用数据主要来源于森林资源清查资料和临安、安吉的典型调查资料，不能够涵盖全省所有毛竹林分的立地状况，这样就降低了回归关系式的使用效果，从而使固碳潜力估算与林分结构调整结果在某些地区偏差较大。