毛竹Phyllostachys edulis林是中国南方一种重要的森林类型，具有生长快、周期短、产量高、用途广等特点。毛竹生物量的研究历来是竹类研究的一项重要内容，也是研究竹林生态系统物质循环的基础，对于评价毛竹的生产力及提高营林水平和综合利用其产品都有着重要的意义^[1]。毛竹林生物量不仅是反映林地生产力的指标，也是度量生态功能如碳储量的重要因子^[2]。立地是森林生长发育的重要物质基础^[3]。范叶青等^[4]研究发现：坡向、坡位对毛竹林分生物量有显著影响。王晨^[5]研究表明：不同海拔高度毛竹林的地上生物量存在显著差异。林分结构是林分功能的基础和表现^[6]。在一定立地和经营管理条件下，毛竹林生物量取决于竹林结构，合理的竹林结构能够提高竹林的产量^[7]。根据是否与毛竹的空间位置有关，毛竹林分结构可分为空间结构和非空间结构。刘恩斌等^[8]研究发现：影响毛竹林分生物量的主要非空间结构因子依次是毛竹株数与林分平均胸径。汤孟平等^[9]研究发现：近自然毛竹林的空间结构与生物量之间存在不可忽视的关系。以上研究表明立地因子和林分结构与毛竹林生物量有密切的关系。但是，现有研究仅从立地因子、林分非空间结构或空间结构单一方面分析毛竹林生物量的主要影响因子，而关于立地因子、林分空间结构和非空间结构对毛竹林生物量的综合影响研究尚少见报道。分析多因子对因变量的综合影响，应当包括确定主导因子以及分析主导因子的直接和间接影响。目前，主要采用相关分析、回归分析、方差分析等方法，确定影响毛竹林生物量的主导因子。方差分析可以判断各因子对毛竹林生物量是否具有显著影响，但无法评估影响大小^[4]。相关分析能够描述各因子与毛竹林生物量之间关系的密切程度，但没有考虑各因子之间的相互关系对生物量的影响^{[8, 10]}。回归分析可以分析因变量与多个自变量的关系，通过回归系数和显著性统计可以揭示各因子对毛竹林生物量的影响大小，但对多重共线性敏感，且理论上的假设限制也较苛刻，如变量独立性、正态性、线性等^[10]。随机森林(random forest，简称RF)是由BREIMAN^[11]提出的一种基于分类树的算法^[12]，可以评估解释变量的重要性^[13]。由于随机森林在理论上的假设限制较回归分析宽松，且具有不容易陷入过拟合、对变量的多元共线性不敏感、对数据集的适应性强等优点^[14]，因此，被广泛应用于评估变量的重要性和变量的筛选^[15-17]。但将随机森林应用于筛选出立地因子、林分空间结构和非空间结构中影响毛竹林生物量的主导因子研究尚少见报道。随机森林虽然考虑了解释变量之间的交互作用^[18]，但随机森林依然无法解析解释变量对因变量的直接、间接和总影响。结构方程模型(structural equation model，简称SEM)整合了路径分析、验证性因子分析与一般统计检验方法，可同时分析一组具有因果关系的方程式^[19]，可以用来揭示解释变量间的复杂关系及其对因变量的直接、间接和总影响^[20-21]。近年来，结构方程模型已广泛应用于分析气候、土壤、林分结构等因素对森林生物量的直接、间接和总影响^[22-23]。因此，将结构方程模型应用于分析影响毛竹林生物量的主导因子对毛竹林生物量的直接、间接和总影响的研究值得探讨。本研究以浙江省少受人为干扰的毛竹林为对象，通过随机森林从立地因子、林分非空间结构因子和空间结构因子中筛选出影响毛竹林生物量的主导因子。在此基础上，建立毛竹林生物量与主导因子的结构方程模型，揭示主导因子对毛竹林生物量的直接、间接和总影响，旨在为毛竹林可持续经营提供理论依据。

1.   研究区与研究方法

1.1   研究区概况

浙江省地处中国东南沿海(27°06′~31°11′N，118°01′~123°10′E)，东临东海，南接福建，西与江西、安徽相连，北与上海、江苏接壤。陆地总面积10.18 万km²，约占全国的1.06%，其中：山地和丘陵占70.4%，平原和盆地占23.2%，河流和湖泊占6.4%，故有“七山一水二分田”之说。在全国植被分区上属于中亚热带常绿阔叶林北部亚地带和南部亚地带，主要植被为针叶林、针阔混交林、阔叶林、灌木林和竹林。地带性植被为中亚热带常绿阔叶林。

根据2017年《浙江省森林资源及其生态功能价值公告》：全省森林605.68 万hm²，森林覆盖率为60.96%。竹林面积91.02 万hm²，占森林面积的15.03%。毛竹林面积79.75 万hm²，占竹林面积的87.62%。

1.2   数据来源

2017−2019年，根据浙江省森林资源一类调查的系统抽样样地中的毛竹林样地分布较多的10个县(市、区)为调查区域(表1)，包括：浙江南部的丽水市庆元县和温州市泰顺县，西部的衢州市常山县，中部的金华市武义县和绍兴市诸暨市，东部的宁波市宁海县和台州市黄岩区以及北部的湖州市安吉县、杭州市临安区和宁波市余姚市。在调查区域内，选择少受人为干扰的近自然毛竹林，设置10 m×10 m的样地共52个，对每个样地进行每木调查，测定胸径、年龄、胸高节长等因子。根据样地每木调查结果，按照胸径分布范围，在大、中、小径级各随机选取1株梢头完整且无病虫害的样竹，共156株。用锯子在样竹根颈处将其伐倒，用皮尺测量竹高。自竹秆基部向上，将竹节从小到大依次编号，称之为节号。用围尺测量每一个竹节的竹节长和竹节中央直径。用游标卡尺测量基部、胸高处和1/2高处东、南、西、北的壁厚，取4个方位壁厚的平均值作为对应位置的壁厚。取3株样竹秆形结构因子的平均值作为样地秆形结构指标。在样地中心挖取1个土壤剖面，用围尺量取土层厚度。取表层土样(0~20 cm)作为样地土壤样品，用自封袋保存好土样，将样品带回实验室进行土壤理化性质测定。pH采用电位法测定；土壤有机质采用浓硫酸-重铬酸钾外加热法测定；碱解氮采用碱解扩散法测定；有效磷采用盐酸-氟化铵法测定；速效钾采用火焰分光光度法测定。

表 1  研究地区概况

Table 1  Overview of the study area

地名	纬度(N)	经度(E)	年均气温/℃	年均降水量/mm	样地数量/个	样竹数量/株
余姚	30°03′	121°09′	16.2	1 361	8	24
临安	30°23′	118°51′	16.4	1 628	4	12
诸暨	29°43′	120°32′	16.3	1 373	6	18
安吉	30°38′	119°40′	15.8	1 420	4	12
宁海	29°29′	121°25′	16.4	1 480	6	18
泰顺	27°30′	119°42′	17.9	1 670	4	12
黄岩	28°38′	121°17′	17.0	1 676	4	12
武义	28°54′	119°48′	17.9	1 546	6	18
常山	28°51′	118°30′	16.3	1 700	6	18
庆元	27°27′	119°30′	17.4	1 760	4	12

下载: 导出CSV 

| 显示表格

1.3   林分结构指数的计算

根据是否与毛竹空间位置有关，毛竹林分结构包括空间结构和非空间结构。本研究选择立竹度(样地株数)、林分平均胸径、郁闭度、年龄均匀度、秆形结构因子作为林分非空间结构指标，以聚集指数、竞争指数和年龄隔离度作为空间结构的基础指标。统计特征见表2。

表 2  毛竹林生物量、结构因子和立地因子描述性统计特征

Table 2  Descriptive statistical characteristics of structure factors of bamboo forests

统计 指标	年龄 均匀度	林分平均 胸径/cm	立竹度/ (株·hm−2)	郁闭度	竞争 指数	年龄 隔离度	聚集 指数	平均高/m	平均 胸高竹 节长/cm	平均 基部 壁厚/mm	平均 胸高 壁厚/mm	平均 1/2高 壁厚/mm	平均 1/2高 节号
均值	0.86	10.51	4 000.00	0.80	5.47	0.59	0.94	14.60	22.70	16.74	12.07	6.57	28.00
标准差	0.10	1.28	1 367.97	0.09	1.38	0.11	0.12	1.68	3.38	1.99	12.37	0.67	2.70
极小值	0.66	7.08	1 600.00	0.60	2.88	0.30	0.73	10.80	16.40	12.49	8.52	5.25	24.00
极大值	1.00	12.94	8 100.00	0.95	8.81	0.84	1.24	18.40	32.00	21.97	99.22	8.33	39.00
变异系数	0.12	0.12	0.34	0.11	0.25	0.18	0.13	0.11	0.15	0.12	1.02	0.10	0.09

下载: 导出CSV 

| 显示表格

统计 指标	平均 竹节数	平均1/2 高竹节 长/cm	平均最 长竹节 长/cm	海拔 高度/m	坡度/ (°)	土层厚 度/cm	pH	土壤 有机质/ (g·kg−1)	碱解氮/ (mg·kg−1)	有效磷/ (mg·kg−1)	速效钾/ (mg·kg−1)	毛竹林 生物量/ (t·hm−2)
均值	60.00	34.70	35.40	35.40	33.60	37.00	5.54	30.90	242.24	21.40	94.65	64.06
标准差	4.44	3.66	4.20	267.97	6.20	14.45	0.36	17.80	109.05	17.81	56.14	23.47
极小值	51.00	25.90	19.00	71.40	20.00	19.50	4.53	14.00	70.00	3.35	34.00	23.26
极大值	72.00	42.10	43.10	1 243.00	46.00	88.00	6.25	88.90	530.88	76.70	327.00	124.60
变异系数	0.07	0.11	0.12	0.81	0.18	0.39	0.07	0.58	0.45	0.83	0.59	0.37

下载: 导出CSV 

| 显示表格

1.3.1   毛竹林分平均胸径

毛竹林分平均胸径计算公式为：

式(1)中：N为样地中毛竹总株数；d_i为第i株毛竹的胸径(cm)；D_g为毛竹林分平均胸径(cm)。

1.3.2   年龄均匀度

毛竹是异龄纯林，考虑到平均年龄无法准确描述毛竹林年龄特征，因此，采用不同龄级毛竹株数的分布均匀程度来代表毛竹林的年龄结构^[24]。毛竹林各龄级株数分布均匀程度可用Shannon均匀度指数描述^[24]。计算公式为：

式(2)中：S为龄级数；P_i是第i龄级的株数比例；A是年龄均匀度。

1.3.3   聚集指数

聚集指数是检验种群空间分布格局的常用指数。计算公式如下所示：

式(3)中：r_i表示第i株毛竹到其最近邻毛竹的距离；F表示样地面积；N表示样地内毛竹总株数；R是林分聚集指数。

1.3.4   竞争指数

采用基于Voronoi图的Hegyi竞争指数^[25]。计算公式如下所示：

式(4)中：d_i是第i株对象竹的胸径；d_j是第j株相邻竹的胸径；L_ij是第i株对象竹到第j株相邻竹的距离；n_i是基于Voronoi图的所有相邻竹株数；I_Ci是第i株对象竹的竞争指数。

取样地内全部竹子竞争指数的平均值作为林分竞争指数。计算公式为：

式(5)中：I_C表示林分竞争指数；N表示样地内毛竹总株数；I_Ci表示第i空间结构单元中对象竹的竞争指数。

1.3.5   年龄隔离度

本研究采用全混交度描述毛竹林的空间隔离度^[26]。全混交度的计算公式为:

式(6)中：Mc_i表示第i空间结构单元中对象竹的全混交度；D_i表示空间结构单位的Simpson指数，$D_{i}=1-\displaystyle\sum\limits_{j=1}^{s_{i}} p_{j}^{2} $，D_i的取值范围为[0, 1]；c_i表示对象竹的最近邻竹中成对相邻竹不是同一竹龄的个数；n_i表示最近邻竹的株数；$\dfrac{{{c_i}}}{{{n_i}}}$表示最近邻竹的年龄隔离度；M_i表示简单混交度，${M_i} = \dfrac{\,1\,}{\,n\,}\displaystyle\sum\limits_{j = 1}^{{n_i}} {{v_{i,j}}} $，如果第i对象竹和第j最近邻竹年龄相同，${{v_{i,j}}}=$0，否则等于1。

取样地内全部竹子年龄隔离度平均值作为林分年龄隔离度。计算公式为：

式(7)中：M表示林分年龄隔离度；N表示样地内毛竹总株数；M_Ci表示第i空间结构单元中对象竹的年龄隔离度。

1.3.6   秆形结构因子

竹秆是毛竹的主体，秆形结构是毛竹研究的重要内容之一。本研究选择8个具有代表性秆形结构指标作为毛竹林分结构指标。9个指标分别为平均高、平均竹节数、平均胸高节长、平均1/2高节长、平均1/2高节号、平均最长竹节长、平均基部壁厚、平均胸高壁厚和平均1/2高壁厚。

1.4   毛竹林生物量

单株毛竹生物量采用以下公式计算^[27]：

式(8)中：W为单株毛竹生物量(kg·株⁻¹)；A为单株毛竹年龄(度)；D为毛竹胸径(cm)。毛竹林生物量指标是在计算单株毛竹生物量的基础上，累加样地中所有毛竹的生物量求得。

1.5   立地因子数量化

以坡度、坡向、坡位、土层厚度、海拔高度、土壤有机质、碱解氮、有效磷、速效钾和土壤pH作为样地立地指标。坡向与坡位为描述性指标，对其进行数量化。坡位分为3个等级：上坡为1，中坡为2，下坡为3^[28]。坡向分为3个等级：阳坡为1，阴坡为2，半阴半阳坡为3^[29]。立地因子描述性统计见表2。

1.6   统计分析方法

采用随机森林和结构方程模型相结合的方法，分析毛竹林生物量与立地因子和结构因子之间的关系。首先，采用随机森林，从立地因子和结构因子中筛选出影响毛竹林生物量的主导因子，然后，利用结构方程模型分析主导因子对毛竹林生物量的直接、间接和总影响。

1.6.1   随机森林

随机森林通过自助法重采样(bootstrap)从n个原始样本集中有放回地重复随机抽取n个样本生成袋内数据，未被采样的数据作为袋外数据(out-of-bag data)。袋外数据作为测试集用于计算袋外误差(out-of-bag error)。袋内数据作为决策树的根节点，在决策树的每个节点处随机选择m个特征(mtry)，并在m个特征中选择1个最具有分类能力的特征进行分裂，决策树最大限度地生长，不做任何修剪，所有决策树组合成随机森林。随机森林可以用于分类和回归，当用于回归时，将所有决策树预测值的平均值作为最终值。本研究通过置换检验，根据袋外观测值和模型预测值之间均方误差的增加量(IncMSE)来评估特征变量的重要性^[13,30]。根据特征变量的重要性及影响显著性判断影响毛竹林生物量的主导因子。随机森林分析均通过R 3.5.3进行，随机森林的建立基于randomForest包完成，变量的显著性和重要性基于rfPermute包完成。

1.6.2   结构方程模型

结构方程模型可以量化变量之间的直接和间接影响，检验整体模型是否具有统计学意义，并提供对数据的准确描述。结构方程模型分析过程主要包括以下3个步骤。①模型设定：基于已有知识和理论分析，预设因子之间的关系，构造初始模型。②模型估计：在估计拟合指标时采用卡方检定，估计参数时采用极大似然估计法。为更好满足结构方程模型的正态性和线性假设，结构方程模型分析所用数据均进行自然对数转换^[23]。③模型评价与修正：根据卡方统计量相对应的P值、拟合优度指数(goodness of fit，I_GFI)、相对拟合指数(comparative fit index，I_CFI)和标准化残差均方根(standardized root mean square residual，R_SRMR)来评价模型，具体的准则为P≥0.05，I_GFI＞0.95，I_CFI＞0.95，R_SRMR＜0.08^[31-32]。若模型拟合效果不佳，根据路径修正系数和理论分析对模型进行修正。结构方程模型的建立及运算均通过软件AMOS 21.0完成。

2.   结果与分析

2.1   基于随机森林筛选影响毛竹林生物量的主导因子

采用随机森林，从立地因子和结构因子中分别筛选出影响毛竹林生物量的主导因子。各特征变量的重要性及影响显著性如图1和图2所示。由图1可见：结构因子中立竹度、竞争指数和林分平均胸径对毛竹林生物量的影响达极显著水平(P＜0.01)，3个特征变量的均方误差增加量分别为16.89%、14.25%和10.68%。由图2可见：立地因子中海拔高度和土层厚度对毛竹林生物量的影响达显著(P＜0.05)或极显著水平(P＜0.01)，2个特征变量的均方误差增加量分别为3.52%和9.33%。这表明立竹度、林分平均胸径、竞争指数、土层厚度和海拔高度是立地因子和结构因子中影响毛竹林生物量的主导因子。随机森林分析结果表明：毛竹秆形结构因子对毛竹林生物量影响较小，可能的原因是受物种遗传特性影响，毛竹林整体秆形结构因子差异较小^[33-34]，因此对毛竹林生物量影响较小。

图 1  结构因子重要性排序

Figure 1  Ranking importance of structure factor

下载: 全尺寸图片 幻灯片

图 2  立地因子重要性排序

Figure 2  Ranking importance of site factor

下载: 全尺寸图片 幻灯片

2.2   基于结构方程模型分析主导因子对毛竹林生物量的影响

立地因子对林分结构有一定的影响^{[1, 35]}，林分结构之间也存在一定的关系^[36-37]。在此基础上，预设5个主导因子之间的影响路径，建立主导因子与生物量关系的结构方程模型，初始模型如图3A所示，模型运行结果见图3B。根据模型评价及检验指标(${\chi ^2}$=4.888，P=0.180，R_SRMR=0.083 8，I_CFI=0.994 0，I_GFI=0.969 0)。可见，初始模型中所预设的路径能够被接受，也能较好体现所采集的数据。由图3B可知：立竹度、林分平均胸径和海拔高度对毛竹林生物量有显著的直接影响，它们之间的标准化通径系数分别为0.914(P＜0.001)、0.739(P＜0.001)和−0.069(P＜0.001)。土层厚度和竞争指数对毛竹林生物量没有显著的直接影响，但通过显著影响其他因子间接影响毛竹林生物量。土层厚度通过显著影响立竹度和林分平均胸径间接影响毛竹林生物量。竞争指数通过显著影响林分平均胸径间接影响毛竹林生物量。竞争指数和立竹度之间的路径系数为0.871(P＜0.001)，说明通过挖笋、伐竹等手段是调节林分竞争水平的有效方式之一。

图 3  毛竹林生物量和主导因子的初始模型及运行结果

Figure 3  Intial model and result between stand biomass and predominant factors

下载: 全尺寸图片 幻灯片

由表3可见：林分平均胸径对毛竹林生物量的总影响最大，达0.739；立竹度对毛竹林生物量的直接影响最大，达0.914；土层厚度对毛竹林生物量的间接影响最大，达0.492；海拔高度、土层厚度和竞争指数对毛竹林生物量的间接影响大于直接影响；海拔高度对毛竹林生物量的总影响为负，即随着海拔的增大，毛竹林生物量逐渐减小，这一结论与范叶青等^[28]研究结果一致；土层厚度对毛竹林生物量的总影响为正，即随着土层厚度的增大，毛竹林生物量逐渐增大，这与许多研究结论一致^{[1, 35]}。土层厚度越大，土壤含水量越高^[38]，有机质含量越高^[39]，毛竹生长发育所需的养分供给越充分，故毛竹林生物量越大；立竹度和林分平均胸径对毛竹林生物量的总影响为正。说明立竹度和林分平均胸径越大，毛竹林生物量越大，这一结论与刘恩斌等^[8]和范叶青等^[4]研究结论一致；立竹度通过竞争指数、林分平均胸径对毛竹林生物量产生负间接影响，间接影响为−0.262。说明随着立竹度的增加，林木内争夺阳光、养分等变得激烈，林分平均胸径呈下降趋势，这一结论与林建忠等^[37]研究结论一致。立地因子对毛竹林生物量的总影响弱于结构因子(表3)，可能是浙江省优越的水热条件削弱了毛竹林生物量在立地因子间的差异^[40]。

表 3  主导因子对毛竹林生物量的标准化直接、间接和总影响

Table 3  Standardized direct, indirect and total effect of predominant factors on stand biomass

潜变量	显变量	直接影响	间接影响	总影响
毛竹林生物量	土层厚度	−0.014	0.492	0.478
	海拔	−0.069	−0.215	−0.284
	立竹度	0.914	−0.262	0.652
	林分平均胸径	0.739	0	0.739
	竞争指数	−0.012	0.325	0.312

下载: 导出CSV 

| 显示表格

3.   讨论

立地条件是影响林木生长的重要因素。为揭示立地条件与毛竹林生物量的关系，以往的研究采用方差分析^[4-5]、偏相关分析^[4]等方法分析立地因子与毛竹林生物量的直接关系，忽视了立地因子对毛竹林生物量的间接影响。本研究通过结构方程模型解析立地因子对毛竹林生物量的影响路径，表明立地因子主要通过影响毛竹林生物量的构件因子(立竹度、林分平均胸径)间接影响毛竹林生物量。

林分结构因子及其相互关系是竹林生产经营过程结构调整的重要依据。刘恩斌等^[8]采用相关分析对立竹度和毛竹林生物量的关系进行了研究，认为立竹度和毛竹林生物量密切相关。本研究中，立竹度与毛竹林生物量之间的路径系数达0.914，表明立竹度对毛竹林生物量的直接影响最大，此结论与刘恩斌等^[8]研究结论一致，但与刘恩斌等^[8]研究方法不同。本研究考虑了林分结构之间的相互影响，并应用结构方程模型量化了立竹度对毛竹林生物量的间接影响，立竹度通过林分平均胸径和竞争指数对毛竹林生物量产生负的间接影响(图3B)，因而总影响(0.652)弱于直接影响(0.914)。说明仅研究立竹度与毛竹林生物量的直接关系，在一定程度上会高估立竹度对毛竹林生物量的正向影响。竞争指数对毛竹林生物量没有显著的直接影响，但通过林分平均胸径对毛竹林生物量产生有正向的间接影响，即毛竹林生物量随着竞争强度增大而增大。可能的原因是，缺乏人工管理的近自然毛竹林其竞争水平普遍较低，而过低的竞争压力毛竹林无法达到高产。维持在适当强度竞争状态，毛竹林生物量才有可能高产^[9]。

毛竹林具有大小年周期的动态变化特征。随着时间的推移，影响毛竹林生物量的主导因子可能发生变化。在进一步研究中，有必要设置固定标准地进行长期监测，建立多个结构方程模型，准确解析主导因子动态变化及其对毛竹林生物量的影响。

4.   结论

本研究以少受人为干扰的毛竹林为研究对象，通过随机森林筛选立地因子和结构因子中影响毛竹林生物量的主导因子。在此基础上，通过构建毛竹林生物量和主导因子的结构方程模型，进一步解析了主导因子对毛竹林生物量的直接、间接和总影响，得出以下主要结论：①毛竹林生物量的形成是一个复杂的过程，受立地、非空间结构和空间结构多因素共同影响，各因素之间并非相互独立，而是存在一定内在联系，共同影响毛竹林的生物量。②采用随机森林，揭示了立地因子和结构因子中影响毛竹林生物量的主导因子，包括：立竹度、林分平均胸径、竞争指数、土层厚度和海拔高度。③结构方程模型表明，立地因子是通过影响林分结构因子对毛竹林生物量产生间接影响，林分结构因子之间的互相影响也对毛竹林生物量有间接影响。④非空间结构因子立竹度和林分平均胸径作为毛竹林生物量的主要构件因子与毛竹林生物量的关系最为密切，立地因子海拔高度、土层厚度和空间结构因子竞争指数主要通过影响毛竹林生物量的构件因子而间接影响毛竹林生物量。⑤在毛竹林经营中，应当在考虑立地因子对毛竹林生物量影响的前提下，调整并保持合理立竹度，调节林分竞争水平，以提高毛竹林生物量。