不同郁闭度下胸高直径对杉木冠幅特征因子的影响

符利勇，孙 华，张会儒，雷相东，雷渊才，唐守正

(1.中国林业科学研究院资源信息研究所，北京 100091;2.中南林业科技大学林业遥感信息工程研究中心，长沙 410004)

树冠是树木进行光合作用的重要场所，它决定树木的生活力和生产力，同时在树木生长过程中也是反映树木的长期竞争水平的重要指标［1］。在树冠结构中，冠幅和冠长率是树冠重要的两个特征因子［2-6］。在单木生长模型中经常用到冠幅和冠长率作为协变量预测树高或胸径生长量、生物量和树木枯损等［4，6］。同时利用冠幅和冠长率计算林木的竞争指数［7-10］。此外，冠幅也是可视化的重要参数［11］。

国内外对冠幅和冠长率的研究主要集中在定性和图表研究［4，7］，部分学者利用传统的回归方法建立冠幅或冠长率与一些林分因子，例如胸高直径、树高、胸高断面积、林分密度等因子的线性关系［12-15］，从而进一步分析这些因子对冠幅的影响。实际应用中，由于经营措施、立地条件，人为因素以及气候等条件的随机干扰，冠幅和冠长率与林分调查因子间可能呈现较复杂的非线性关系。与此同时，所调查的数据常为重复调查数据或多水平数据，如不同等级郁闭度中对树木冠幅或冠长率重复观察，这些数据的调查对象间可能存在有明显的自相关和异方差等［16］。而回归分析方法是假定数据间相互独立且非异质性［17-19］，是反应林分总体变化情况，对于研究对象(林分类型)个体间的差异程度，例如不同郁闭度等级对冠幅和冠长率的随机影响，传统的回归方法无法解决。为此，本论文以湖南省黄丰桥国有林场103块样地共2461株杉木为例，详细介绍利用混合模型方法分析不同郁闭度等级对冠幅和冠长率的影响，同时还得到冠幅和冠长率的预测模型。文中所有计算是在S-Plus软件nlme模块上实现。

1 材料与方法

1.1 研究区自然概况

研究区位于湖南省黄丰桥国有林场，该林场呈带状分布横跨株洲市攸县东西部，介于东经113°04'至113°43'，北纬27°06'至27°24'之间。东北部与江西的莲花、萍乡交界，东南与茶陵县接壤，西北部与株洲、醴陵毗邻。全场地貌以中低山为主，境内最高海拔1270 m，最低海拔115 m，坡度介于20—35°度之间。林场地处中亚热带季风湿润气候区，年均气温17.8℃，平均无霜期为292d，年均降水量1410.8 mm。全场现有林地面积10122.6 hm2，活立木蓄积879705 m3，其中有林地蓄积879688 m3，占活立木蓄积99.99%，四旁树蓄积17 m3，仅占0.01%。林场的森林覆盖率为86.24%。树种主要以杉木(Cunninghamia lanceolata(Lamb.)Hook)、松类(Pinus spp)为主，其中杉木面积3197.6 hm2，占用材林面积89.9%，蓄积593738 m3，占96.56%，全部为人工造林。

1.2 研究数据与方法

1.2.1 研究数据

按林分郁闭度(CD)大小将黄丰桥国有林场划分为5种不同林分类型:类型1 CD≥0.8;类型2 0.7≤CD＜0.8;类型3 0.6≤CD＜0.7;类型4 0.5≤CD＜0.6;类型5 CD＜0.5。在这5 种林分类型中分别随机抽取6、13、14、34和36块共103块样。样地最小面积为0.0066 hm2，最大面积为0.1052 hm2。每个样地调查内容为胸径5 cm以上的活立木胸径、树高、枝下高以及东、南、西、北四个方向的冠幅长度，共调查2461株。

1.2.2 指标测定

(1)冠幅(CW)

研究所定义冠幅为东、南、西和北四个方向的平均值，计算公式为［11］:

式中，CW东、CW西、CW南和CW北分别代表东、西、南、北冠幅长度(m)。

(2)冠长(CL)

冠长即树冠的长度，指树干第一个活枝到树梢的长度。

式中，H为树高(m)，˜H为树干第一个活枝高度(m)，简称枝下高。

(3)冠长率(CR)

冠长率为冠长与树高的比值［20］。

CR的值在0到1之间，当CR=0时表明树干没有冠幅，CR=1时为全树冠，这两种情形非常少见，将冠长率处于0＜CR＜1外的样木将全部剔除。

(4)大于对象木树木平均直径(MDL)

大于对象木(被分析的树木)树木平均直径是指在样地中大于对象木所有树木直径的平均值，反映树木在林分中竞争强度大小［16］。

式中，MDLi为大于第i株对象木所有树木的平均直径，ni为样地中大于对象木直径的树木株树。

1.2.3 数据分析

(1)基础模型

论文选用Logistic模型分析冠幅和冠长率与林木胸高直径(D)关系［10］。为降低树木间由竞争所产生的随机效应对冠幅和冠长率影响，模型中将增加因变量MDL，表达式如下:

式中，y为冠幅(CW)或冠长率(CR);φ1—φ4为待估参数。

(2)混合模型

Pinheiro和Bates［23］给出的非线性混合效应模型(称NLMEMs)表达式为:

式中，M为郁闭度等级数;ni为第i个郁闭度等级重复观测次数;yij表示第i个郁闭度等级中第j次观测的因变量值(此次研究为冠幅或冠长率)，f是关于参数向量φij和连续变量υij的非线性函数;β为p×1维固定效应参数;ui为q×1维的随机效应参数，假定服从期望为0方差为ψ的正态分布;φij为形式参数(简称形参)，它与β和 ui呈线性函数关系;Aij与 Bij分别为β和 ui的设计矩阵;εij是随机误差项，假定对于所有的i和j，εij都服从期望为0方差为R的正态分布，并假定随机效应参数 ui与误差项εij之间相互独立。关于混合模型的详细介绍，见符利勇等［22］，符利勇和唐守正［23］等。

(3)模型评价

利用模型(6)分析树冠和冠长率时，在103块样地中随机抽取69块样地作为建模数据，剩下34块样地作为检验数据。评价指标为平均残差(¯e)、残差方差(η)和均方误差(δ)，相应的计算公式如下:

2 结果与分析

2.1 冠幅和冠长率结构

2.1.1 各林分类型基本信息

种林分类型年龄、密度、胸径、树高和枝下高的统计信息见表1。从表中得知，所有调查样地中，林分年龄分布在5—35a间，其中各类型林分平均年龄非常接近(类型1，2，4为19a，类型3和5为20a)。对于密度，各调查样地间相差较大，分布在237.62—1719.272株/hm2间，而林分类型间平均密度相差较小。对于胸径，各调查样地主要分布在5—39.5 cm间，而林分类型间平均胸高相差较小。对于树高，各调查样地主要分布在4.3—25.8 m间，而林分类型间相差不大。对于枝下高，各调查样地主要分布在1—18.2 m间，同样各林分类型间相差不大。

表1 5种林分类型统计信息Table1 Summary statistics for 5 stand types

2.1.2 冠幅结构

各类型冠幅大小基本信息见表1。从表1中得知5种类型的冠幅大小在区间(0.65—7 m)之间，其中，类型5的平均冠幅最大，为3.85 m，类型2最小为3.07 m，5种类型的总平均冠幅为3.37 m。5种类型冠幅变异系数分别为29.13、21.82、25、24.93和25.45。为清楚地描述冠幅的结构分布，以0.5 m为级距对冠幅进行分级(尝试多种级距，当级距很小时，多个级距间冠幅分布为0，很大时冠幅集中分布到几个级距间，综合考虑级距为0.5时能较好反映冠幅结构)，其中大于5 m为上界，图1为各类型中冠幅等级占观测数的百分比分布。显然5种类型的冠幅大小主要分布在区间(2.5—4.0 m)之间，分别占 65.82%、69.56%、70.79%、58.15%和53.21% ，小于2.0 m和大于5.0 m分布较少。

图1 各类型中冠幅等级占观测数的百分比分布Fig.1 Crown width classes accounted for the percentage distribution of observations for each type

2.1.3 冠长率结构

各类型冠长率基本统计量见表2。从表2中得知，5种类型的冠长率分布在区间(0.02—0.89)之间，总平均冠长率为0.47，其中类型2的冠长率最大为0.52，类型3最小为0.44。对于变异系数，类型1最大，为45.32，类型4最小为29.98。与冠幅一样，对冠长率按0.1级距进行分级(级距选择方法与冠幅类同)，图2为各类型中冠长率等级占观测数的百分比分布。从图中得知5种类型冠长率主要分布在区间(0.3—0.7)之间，分别占55.69% 、71.93% 、67.01% 、82.22%和79.28%。在区间小于0.2 和大于 0.8 冠长率分布较少，大于0.9的冠长率几乎没有。

表2 5种类型各自冠长率统计量Table2 Summary of live crown ratio for 5 types

2.2 各类型冠幅和冠长率与胸径关系

2.2.1 各类型冠幅与胸径关系

图3为各类型中树木冠幅大小与胸径的散点分布图，同时也给定了冠幅与胸径的线性回归方程和决定系数。图中看出冠幅大小与胸径总体呈正相关，其中，类型1的决定系数最大，类型5次之，最小是类型2。5种类型的树木冠幅与胸径有一定的线性相关性，但通过相关性检验得知在可靠性α=0.05时都相关不显著。类型1具有较大的相关系数ρ=0.74可能与观测数较少有关。因此说明冠幅大小与胸径可能存在较复杂的非线性关系。

图2 各类型中冠长率等级占观测数的百分比分布Fig.2 Live crown ratio classes accounted for the percentage distribution of observations for each type

图3 各类型冠幅大小与胸径散点分布Fig.3 Scatters distribution of crown width and diameter for each type

2.2.2 各类型冠长率与胸径关系

图4为5种类型的树木冠长率与胸径的散点分布图。从图中散点走势得知，冠长率与胸径大小总体而言呈反比，类型2和类型3最明显。同时从图中还可得知各类型的决定系数都非常小，尤其是类型2几乎为零，从而说明冠长率与胸径大小无线性关系，通过相关系数检验也进一步得到证实。因此对于各种类型，冠长率可能与胸径呈复杂的非线性关系。

2.3 冠幅和冠长率模型

2.3.1 冠幅模型

从上节得知冠幅和冠长率与胸径线性相关性较弱，因此论文选用模型(5)对它们进行拟合。由于考虑郁闭度对冠幅和冠长率的影响，故把郁闭度作为随机效应因子构造混合模型。对作用在模型(5)上的15种不同形式参数构造类型进行分析与比较(评价指标AIC和BIC越小越好，参见文献［21］)得知，当由郁闭度产生的随机效应同时作用在形式参数φ2和φ3上时评价指标最小，AIC=3770.23，BIC=3819.25，因此选择该模型分析郁闭度对冠幅的影响［26］，模型表达式为:

图4 各类型冠长率大小与胸径散点分布Fig.4 Scatters distribution of live crown ratio and diameter for each type

其中，M=5;CWij是郁闭度第i等级第j次观测树木的冠幅大小(m);β1—β4为固定效应参数;ui=(u2i，u3i)为郁闭度第i等级所产生的2×2维随机效应，ψ为ui的方差;εij为误差项，方差为σ2;对于任意的i和j，假定 ui和 εij相互独立;并假定 εi=(εi1，…，εini)T服从期望为0方差为 Ri的正态分布［25-26］。

Gi是用来描叙对象内方差异质性的ni×ni维对角矩阵;Γi是用来描叙对象内误差自相关性的ni×ni维矩阵。由于本实例中各样地内观测数据之间没有明显的相关性，故Γi为单位矩阵。为解释异方差性，采用指数函数作为残差方差模型［21］:

式中，γ为待估常数。

利用建模数据对模型(11)进行计算得到各参数估计值见表3。模型(11)中，给定D和MDL，冠幅与随机效应参数u2i和u3i呈正比，即u2i和u3i越大冠幅越大。从表3中得知，类型5对应的随机效应参数最大，类型3次之，最小的是类型1，从而说明在相同的胸高直径和竞争条件下，类型5的冠幅最大，而类型1的冠幅最小。利用模型(11)对检验数据进行预测时，得到的评价指标¯e=0.3124、η=1.0132、δ=1.0603，而模型(5)为¯e=0.5731、η=1.3201、δ=1.4391，因此说明模型(11)有较高的预测精度，同时也进一步证实混合模型比传统的回归模型精度要高［22-23］。

表3 模型(11)参数估计值Table3 Parameter estimates of model(11)

2.3.2 冠长率模型

与上节相似，同样利用基础模型(5)构建冠长率混合模型，随机效应因子为郁闭度。通过对15种不同形式参数构造类型进行分析与比较得知当随机效应同时作用在形参φ1和φ4上时对应的AIC和BIC最小，AIC=－1524，BIC=－1474.99，因此选择该模型分析郁闭度对冠长率的影响，模型表达式为:

式中，RCVij为郁闭度第i等级中第j次观测树木的冠长率，ui=(u1i，u4i)为郁闭度第i等级所产生的随机效应。模型(14)中除随机效应作用的形参不同外(作用在形参φ1和φ4上)，模型中参数结构和分布与模型(11)中定义完全相同，本节不再重复阐述。

利用建模数据对模型(14)进行计算得到参数估计值(表4)。从表4中得知，类型1对应的随机效应参数最大，类型2次之，逐渐递减，最小是类型5，从而说明在相同的胸高直径和竞争条件下，类型1的冠长率最大，而类型5的冠长率最小。同样利用模型(14)对检验数据进行预测时，得到的结论与上节完全相同，即相对回归模型，混合模型(14)有较高的预测精度，在实际应用中可以利用该模型预测冠长率。

表4 模型(14)参数估计值Table4 Parameter estimates of model(14)

3 结论与讨论

(1)从表1中得知各类型中杉木的冠幅主要分布在(2.5—4.0 m)之间，对于小于1.5 m和大于5.0 m分布较少，这与所调查的林分特性及杉木本身有关。类型1的变异系数比其它类型要大，可能存在以下几个原因:

1)在郁闭度较小的林分生长空间中，更新的小树较多，相比与优势木的冠幅，进界幼树的冠幅小得很多;

2)从表1中看出类型1的最大冠幅为5.05而最小冠幅为0.65，因此林分中郁闭度较小可能是由于林分中优势树种呈团状分布，从而使得林分中林隙较多，进界的幼树较多;

3)可能是由于类型1中观测数较少产生。对于冠长率类型1的变异系数最大而类型5变异系数较少，这可能是由于在郁闭度小的林分中，存在一定数量的幼树和优势木，它们对应的冠长率根据树种特性和立地条件差异较大，而在郁闭度较大的林分中，对于一些中庸木和亚优势木可能由于竞争不过优势木而死亡，从而存活下来的主要是优势树种，其冠长率差异较小(表2)。

(2)通常冠幅大小与直径有着显著的相关性［11］，从图1中看出类型1的冠幅与胸径相关系数较高(R=0.7430)，而其它类型相对较小，这可能是由于在郁闭度较小的林分空间中，林木竞争关系相对较弱，树木可以无约束的生长，因此冠幅与胸径有明显的关系。对于类型2，冠幅与胸径相关系数最小(R=0.5292)，这可能是由于在该郁闭度等级中，树种存在一定的竞争关系，同时局部的立地条件对树木的生长影响较大。因此对于胸径较大的树种，如果周围存在较多的竞争树种，其冠幅可能由于自身竞争能力弱而影响生长，相反对于幼树，如果立地条件好周边又没有竞争树种，其冠幅也可能较大。对于类型5，由于郁闭度较高，因此成活的树种一般为优势树种，此时树种的冠幅与胸径的影响不会很大(图1)。对于冠长率与胸径关系，与冠幅完全相似。

如果把所有调查样地按照相同林分年龄分类并分析不同年龄大小间冠幅和冠长率与直径间的关系，结果得知，它们之间有较大的相关性。因此说明，冠幅和冠长率与胸径线性相关较弱的一部分原因可能还与同一林分类型中年龄之间差异较大引起。除此之外，冠幅和冠长率的生长还受经营措施影响，例如间伐次数和间伐强度等。历史资料显示，在本研究区域中，根据林场生产实际和林分生长状况，共抚育间伐分两次。第1次为生长伐，间伐的林龄为9—11a;第2次为定株伐，间伐的林龄为16—18a，两次间伐强度为林分株数的25%左右。冠幅和冠长率与直径不线性相关可能与间伐有一定关系，目前作者正在收集研究区域林分经营的历史资料，尚需进一步研究。同时关于其它林分因子对冠幅和冠长率的具体影响，将做进一步研究。

(3)除用定性的方法描述郁闭度对冠幅和冠长率影响外，此次研究还应用的非线性混合效应模型方法。该方法为常见的一种用来处理重复调查数据或纵向数据的工具［21］。通过对冠幅和冠长率分析得到，在相同的胸高直径和竞争条件下，类型5的冠幅最大，而类型1的冠幅最小。对于冠长率，恰好相反，即类型1的冠长率最大，而类型5的冠长率最小。因此利用此次研究结论可以对林分生长收获预测和经营决策提供依据。冠幅和冠长率通常与多个林分因子有显著影响，例如林分密度、林木竞争指数、枝下高等［10-11］，此次研究是在不同等级的郁闭度下分析冠幅和冠长率与胸高直径以及大于对象木树木平均胸高直径之间关系，通过检验数据验证，论文给出的冠幅模型和冠长率模型预测精度较高，在实际应用中可以利用它们对冠幅和冠长率进行预测。

［1］ Biging G S，Dobbertin M.A comparison of distance-dependent competition measures for height and basal area growth of individual conifer trees.Forest Science，1992，38(3):695-720.

［2］ Assmann E.The Principles of Forest Yield Studies.Oxford:Pergamon Press，1970:506-506.

［3］ Spurr SH，Barnes B V.Forest Ecology.3rd ed.New York:John Wiley，1980:527-527.

［4］ Monserud R A，Sterba H.A basal area increment model for individual trees growing in even-and uneven-aged forest stands in Austria.Forest Ecology and Management，1996，80(1/3):57-80.

［5］ Kalliovirta J，Tokola T.Functions for estimating stem diameter and tree age using tree height，crown width and existing stand database information.Silva Fennica，2005，39(2):227-248.

［6］ Russell M B，Weiskittel A R.Maximum and largest crown width equations for 15 tree species in Maine.Northern Journal of Applied Forestry，2011，28(2):84-91.

［7］ Wykoff W R.A basal area increment model for individual conifers in the Northern Rocky Mountains.Forest Science，1990，36(4):1077-1104.

［8］ Biging G S，Dobbertin M.Evaluation of competition indices in individual tree growth models.Forest Science，1995，41(2):360-377.

［9］ Monleon V J，Azuma D，Gedney D.Equations for predicting uncompacted crown ratio based on compacted crown ratio and tree attributes.Western Journal of Applied Forestry，2004，19(4):260-267.

［10］ Toney C，Reeves M C.Equations to convert compacted crown ratio to uncompacted crown ratio for trees in the interior west.Western Journal of Applied Forestry，2009，24(2):76-82.

［11］ Lei X D，Zhang Z L，Cheng X G.Crown-width prediction models for several tree species including Larix olgensis in northeastern China.Journal of Beijing Forestry University 2006，28(6):75-79.

［12］ Gill SJ，Biging G S，Murphy E C.Modeling conifer tree crown radius and estimating canopy cover.Forest Ecology and Management，2000，126(3):405-416.

［13］ Brabg D C.A local basal area adjustment for crown width prediction.Northern Journal of Applied Forestry，2001，18(1):22-28.

［14］ Bechtold W A.Crown-diameter prediction models for 87 species of stand-grown trees in the eastern United States.Southern Journal of Applied Forestry，2003，27(4):269-278.

［15］ Bechtold W A.Largest-crown-width prediction models for 53 species in the western United States.Western Journal of Applied Forestry，2004，19(4):245-251.

［16］ Uzoh F C C，Oliver WW.Individual tree diameter increment model for managed even-aged stands of ponderosa pine throughout the western United States using a multilevel linear mixed effects model.Forest Ecology and Management，2008，256(3):438-445.

［17］ Biging G S.Improved estimates of site index curves using a varying-parameter model.Forest Science，1985，31(1):248-259.

［18］ Keselman H J，Algina J，Kowalchuk R K，Wolfinger R D.A comparison of recent approaches to the analysis of repeated measurements.British Journal of Mathematical and Statistical Psychology，1999，52(1):63-78.

［19］ Garrett M F，Laird N M，Ware J H.Applied Longitudinal Analysis.New Jersey:John Wiley and Sons，Inc.，2004.

［20］ Monserud R A.Methodology for Simulating Wisconsin Northern Hardwood Stand Dynamics［D］.Madison:University of Wisconsin，1975:156-156.

［21］ Pinherio JC，Bates D M.Mixed-Effects Models in Sand S-PLUS.New York:Spring-Verlag，2000.

［22］ Fu L Y，Tang SZ.Dominant height for Chinese fir plantation using nonlinear mixed effects model based on linearization algorithm.Scientia Silvae Sinicae，2012，48(7):66-71.

［23］ Fu L Y，Li Y C，Li C M，Tang SZ.Analysis of the basal area for Chinese fir plantation using two kinds of nonlinear mixed effects model(two levels).Scientia Silvae Sinicae，2012，48(5):36-43.

［24］ Pinheiro J C，Bates D M.Approximations to the log-likelihood function in the nonlinear mixed effects model.Journal of Computational and Graphical Statistics，1995，4(1):12-35.

［25］ Davidian M，Giltinan D M.Nonlinear Models for Repeated Measurement Data.New York:Chapman＆ Hall，1995.

［26］ Meng SX，Huang SM.Improved calibration of nonlinear mixed-effects models demonstrated on aheight growth function.Forest Science，2009，55(3):239-248.

参考文献:

［11］ 雷相东，张则路，陈晓光.长白落叶松等几个树种冠幅预测模型的研究.北京林业大学学报，2006，28(6):75-79.

［22］ 符利勇，唐守正.基于非线性混合模型的杉木优势木平均高.林业科学，2012，48(7):66-71.

［23］ 符利勇，李永慈，李春明，唐守正.利用2种非线性混合效应模型(2水平)对杉木林胸径生长量的分析.林业科学，2012，48(5):36-43.