基于多元统计分析的马尾松人工林健康评价研究——以广西热带林业实验中心为例

赵勇钧，谢阳生，王建军，李 晗，葛方兴，孟京辉

（1.北京林业大学 森林资源和环境管理国家林业局重点开放性实验室，北京 100083；2.中国林业科学研究院资源信息研究所，北京 100091）

森林的健康状况直接关系到全球生态环境安全和人类社会的可持续发展[1]，如何开展森林健康经营，实现森林生态系统的恢复和森林多功能的发挥，已成为世界关注和研究的焦点领域[2]。森林健康评价是进行森林健康经营的基础[3]，健康的森林不能简单地认为是没有病虫灾害，而是将其维持到较低水平，更好维护健康森林中的食物链和生物多样性，保持森林结构的稳定性[4-5]。

评价指标的选取和评价方法是森林健康研究的核心内容。迄今，国内外对此进行了大量研究，如Costanza[6]认为健康的生态系统应该体现出完整性、稳定性和功能性，提出系统活力、系统组织和系统恢复力3 个指标；美国林业局与环境保护署提出初级生产力、龄级分布、土壤养分含量、树冠情况、林下植被多样性、病虫害等指标对全国森林进行健康评价[7]；肖风劲[8]提出森林生态系统的生态要素、生理要素、胁迫要素、环境要素和气象要素5 个方面共19 个指标；甘敬[9]基于BP 神经网络确立林分层次结构、病虫害程度和土壤厚度3 个指标为森林健康快速评价指标。目前在健康指标选取过程中，由于研究尺度和研究方法以及对森林健康的标准存在差异，不同学者所考虑的指标也大多不一样，存在有些指标过于理论化难以量化，指标之间的信息重叠，指标过少很难有效反映森林生态系统健康的真实状况，指标过多导致收集应用困难等情况，且不同区域、不同森林类型都会导致选取关键指标和标准的差异。由于森林自身的复杂性，使得森林健康的评价难度较大，迄今尚未存在统一的标准方法，在现有的森林健康评价方法中，应用最多的是层次分析法、主成分分析法、综合指数法、模糊综合评判法、聚类分析法、指示物种法、人工神经网络法、健康距离法等[10]，对于同一块森林生态系统的健康程度而言，不同的评价方法对森林健康评价的结果也不同，且评价方法各有其优缺点，关于权重的确定方面，有些方法比较主观，如健康距离法、综合指数评价法；在定性定量方面，有的方法是定量分析，如主成分分析法，有的是将两者结合，如层次分析法、模糊综合评判法。目前森林健康评价越来越呈现出多种方法相结合的趋势。

因此，评价方法与准则的确定以及实用性是目前森林健康评价中存在的主要问题。本研究针对我国南方主要造林树种马尾松Pinus massoniana，以广西热带林业实验中心的马尾松人工林为例，在科学选择指标的基础上，建立因子分析、聚类分析和判别分析综合量化评价体系，以期为马尾松人工林的健康评价提供方法和技术支撑，为该地区马尾松人工林健康经营提供科学依据。

1 研究区概况与数据来源

研究区位于广西壮族自治区凭祥市的中国林业科学研究院热带林业实验中心（以下简称“ 热 林 中 心”，21°57′47″～22°19′27″N，106°39′50″～ 106°59′30″E）。气候属于南亚热带半湿润- 湿润季风性气候，日照充足，雨量充沛，年降水量1 062～1 772 mm，年蒸发量1 261～1 388 mm，年均气温21～23℃，极端高温40.3℃，极端低温-1.5℃。区内地形主要以丘陵为主，平均海拔500～800 m，土壤以砖红壤、红壤为主。区内物种资源非常丰富，乔木树种以马尾松Pinus massoniana、杉木Cunninghamia lanceolata、 米 老 排Mytilaria laosensis、 红 椎Castanopsis hystrix、灰 木 莲Manglietia glauca等为主，灌木主要是山苍子Litsea cubeba、盐肤木Rhus chinensis、三叉苦Evodia lepta等，草本主要是蔓生绣竹Microstegium gratun、铁线蕨Adiantum capillas-veneris、五节芒Miscanthus floridulu等。

2011年在热林中心区域面积内按照系统抽样方法在公里网（1 km×1 km）交叉点布设监测样地，共设置238 个星状圆形样地（图1）。样地设计如图2所示，样地面积为400 m2，在每个样地内设置3 个半径R=6.51 m 的圆形子样（A、B、C），子样圆内分别设一个4 m×4 m 幼树灌木样方（a、b、c）和一个1 m×1 m 的幼苗草本样方（1、2、3）。在样地中心（O）记录海拔、坡度、坡向等立地因子，在设置的3 个圆形样地中对所有胸径≥5 cm 的乔木进行每木检尺，记录树种、胸径、树高、郁闭度、起源、干型质量等，对幼树灌木层调查记录灌木平均高、株数及盖度，对幼苗草本样方调查记录更新幼苗数和草本植物种类、数量、高度、盖度。在每个乔木样地中心点附近8.49 m2范围内选取土壤剖面，按照土壤发生层进行取样，采取O 层、A 层和B 层土壤样品，采样时小心除去土壤表层枯落物，在不同土层的不同部位用环刀取土，每层取3 个土样，然后将土样混合均匀装入采样袋内。本研究从238 个系统抽样样地中选取马尾松为优势树种的84 个人工林样地，其中幼龄林12 个，中龄林34 个、近熟林27 个、成熟林11 个，马尾松人工林样地重新编号为M1～M84。

2 研究方法

2.1 评价指标选取及量化

图1 热带林业实验中心样地分布Fig.1 Distribution of sampling plots of tropical forestry experimental center

图2 热带林业实验中心系统抽样样地设计Fig.2 System sample plot design of tropical forestry experimental center

本研究依据森林生态系统的健康特征[11]，遵循科学性、客观性、可操作性的原则[12]，选取可获取性、可量化性、有代表性的评价指标13 个，包括林分平均胸径、平均树高、每公顷蓄积、地位指数、密度、郁闭度、树种组成指数、灌草盖度、腐质层厚度、土壤厚度、土壤有机质含量、灾害程度、天然更新等级。其中，灾害程度和天然更新等级两个定性指标根据国家森林资源连续清查调查技术规程[13]，同时参考甘敬、施明辉等[14-15]对定性指标的量化标准将其转化为数值参与评价（表1）。剩余的11 个定量指标则直接参与评价，部分指标计算公式如下：

表1 定性指标量化标准Table1 Qualitative indicator quantitative standard

树种组成指数用样地内乔木树种Simpson 多样性指数表示[16]：

式中：S为调查样地林分中树种Simpson 多样性指数，即树种组成指数，p为样地内树种个数；Ni为第i个树种占样地中所有树种的比例。

地位指数的计算采用广西马尾松人工林多形曲线地位指数模型[17]：

式中：DH为林分优势高，T为林分年龄，T0=20 a。

2.2 健康评价指标体系的构建

本研究基于多元统计分析中的因子分析法、聚类分析法和判别分析法[18]对热林中心的马尾松人工林进行健康评价研究。通过对各样地的指标数据进行因子分析，提取共性因子，分别计算各因子得分值和各样地健康指数；对每个样地的健康指数值进行聚类分析划分健康等级，得到各样地的健康状况分类；同时基于每块样地各因子得分值建立判别函数，利用判别模型对各样地进行健康状况归类，通过归类对聚类分析结果进行回判验证。

具体步骤：1）对指标数据进行KMO（Kaiser-Meyer-Olkin）和Bartlett 球形度检验[19-20]，以验证数据是否适合因子分析。

2）查看因子分析结果以及解释的总方差，按照特征值（主成分对原变量的解释能力）大于1的原则提取主成分，得到因子得分系数矩阵。

3）对原始数据标准化处理。评价指标之间由于量纲不统一，缺乏可比性，因此需要将原始数据进行标准化处理。评价指标分为3 类：①正向型指标。数值越大，该指标森林健康程度越高，如平均胸径、平均树高、地位指数等。②逆向型指标。数值越小，该指标森林健康程度越高，如灾害程度，本研究在定性指标量化过程中已将逆向指标正向化（表1）。③极值型指标。数值在某一区间内时，该指标森林健康程度最佳，离此值距离越远，健康程度越低，如密度、郁闭度。由于极值型指标在本试验中影响不大，故指标都采用同样的标准化处理方法。本研究采用“中心化”处理法对评价指标进行标准化处理，计算公式如下：

式中：Uij为第i个指标的第j个样地观测值的标准化数据，Xij为第i个指标第j个样地观测值，、Si分别为第i个指标的平均值和标准差。

4）计算各因子得分值和各样地健康指数。计算公式如下：

式中：FHI 为各样地健康指数，Fjm为第j个样地第m个因子的得分，Wm为第m个因子的权重系数，Dmi为第i个指标的第m个因子得分系数，Cm为因子分析后得到各因子的贡献率。

5）采用Ward 聚类法对各个样地健康指数值进行聚类分析划分健康等级。根据热林中心的特点和实际情况，参考国家森林资源连续清查调查技术规程[13]中关于森林健康等级的划分标准以及现有研究成果[21-22]，将热林中心马尾松人工林健康划分为健康（Ⅰ）、亚健康（Ⅱ）、不健康（Ⅲ）3 个等级。

6）根据每块样地的各个因子得分值建立判别函数，利用Fisher 判别法，建立马尾松人工林健康判别模型，分别用自身验证法和交互验证法进行回判验证。

数据分析在Excel 2013 和SPSS 20.0 中进行。

3 结果与分析

3.1 因子分析

由表2可知，选取的评价指标经过检验之后，Bartlett 球形度统计量为536.027，相应的显著水平为0.000 ＜0.05，说明相关系数矩阵与单位阵有显著差异。同时，KMO 系数为0.603 ＞0.5，根据KMO 度量标准[20]可知实验数据可以用于因子分析。

表2 KMO取样足够量数和Bartlett球形度检验Table2 KMO sampling sufficient quantity and Bartlett sphericity test

本研究共有13 项指标参与分析，按照特征值（主成分对原变量的解释能力）大于1 的原则，共提取5 个主成分，如表3所示，旋转后方差累积方差贡献率为72.68%，总体上反映了原始数据的大部分信息，因子分析效果较为理想。

旋转成分矩阵中的各元素被称为因子载荷，代表了指标与主成分之间的相关系数，其绝对值越大，说明相关程度越高，由表4可以看出，平均胸径、平均树高、每公顷蓄积、地位指数在第一主成分的载荷明显大于其他指标，说明上述4 个指标之间的相关性较高，因此可视为一组，这些指标主要体现了森林生产力信息，可以将其理解为森林生产力指标；同理，第二主成分中的腐殖质厚度、土壤厚度、土壤有机质含量构成一组，体现了林地土壤信息，可以将其视为土壤指标；第三主成分中的密度和郁闭度构成一组，第四主成分中的灾害危害情况、天然更新、灌草盖度构成一组，第五主成分中的树种组成指数构成一组，这3 个主成分都同样体现了生态系统的稳定性，可以解释为生态系统稳定性指标。本研究构建的马尾松人工林健康评价指标体系如表5所示。

表3 各成分的特征值†Table3 Characteristic values of each component

表4 旋转成分矩阵Table4 Rotation component matrix

根据式（3）、式（4）、式（5）、式（6）和表4，可以计算出各样地的5 个因子得分值和综合得分值，如表6所示。

3.2 聚类分析与综合评价

基于各样地因子综合得分FHI，采用Ward 聚类，把84 个样地划分为3 个聚类组（图3）。Ⅰ组包括27 个样地，FHI 值0.249～1.128，为健康状态的样地，占调查样地的32.1%；Ⅱ组包括30块样地，FHI 值-0.270～0.213，为亚健康状态的样地，占调查样地的35.8%；Ⅲ组包括27 块样地，FHI 值-0.287～-1.107，为不健康的样地，占调查样地的32.1%（表5和图3）。表明热林中心马尾松人工林健康状况总体上分布较为均匀，但非健康林分仍占67.9%，健康状况不容乐观，需要进行健康抚育经营。

表5 马尾松人工林健康评价指标体系Table5 Health evaluation index system of Pinus massoniana plantation

表6 各样地因子得分Table6 Scores of various plot factors

图3 84 个马尾松人工林样地Ward 聚类图像Fig.3 Ward clustering map of 84 Pinus massoniana plantations

根据龄组数据，将热林中心马尾松人工林按幼龄林、中龄林、近熟林、成熟林4 个龄组统计其健康情况，如表7所示。为了更直观地表示分布，将数表用柱状图表示，从图4可以看出，幼龄林中无健康林分，亚健康林分占16.7%，不健康林分占83.3%；中龄林中健康林分占26.5%，亚健康林分占44.1%，不健康林分占29.4%；近熟林中健康林分占59.3%，亚健康林分占29.6%，不健康林分占11.1%；成熟林中健康林分占18.2%，亚健康林分占45.5%，不健康林分占36.3%。整体上不同龄组马尾松人工林健康状况从优到劣的排序依次为：近熟林＞中龄林＞成熟林＞幼龄林。

表7 不同龄级马尾松人工林健康状况评价统计Table7 Statistics on health assessment of different ages of Pinus massoniana plantations

图4 不同龄级马尾松人工林健康评价结果Fig.4 Health evaluation results of different ages of Pinus massoniana plantations

3.3 判别分析

根据各样地5 个主因子得分（表6）和聚类分析结果（图3），采用Fisher 线性判别以检验聚类分析结果的准确性。从表8可以看出，第一判别函数解释了96.4%的方差，第二判别函数解释了3.6%的方差，两个判别函数解释了全部方差。第一判别函数与因变量具有很强的相关性，相关系数（r）为0.909，说明第一判别函数是重要的分类维度，第二判别式与因变量的相关系数（r）为0.390，利用两个判别式可以较好地区分各类，如图5所示。表9是对两个判别函数的显著性检验，由Wilks’Lambda 检验结果，认为两个判别函数在0.05 的显著水平上是显著的。

表8 Fisher判别函数特征值Table8 Fisher discriminant function eigenvalues

图5 基于Fisher 判别函数的各类散点图像Fig.5 Various scatter plots based on Fisher discriminant function

表9 Wilks’ Lambda检验Table9 Test of Wilks’ Lambda

经过判别分析得到标准化的Fisher 判别函数系数（表10），由此得到判别函数的表达式为：

将各样地的主成分得分f1、f2、f3、f4、f5分别代入判别函数表达式（7）（8）中计算函数值，比较它们距各类中心的距离，然后待判样本被划分到距离最小的那一类，可以得到各个样地的判别分类结果，如表11所示，自身验证法对Ⅰ～Ⅲ健康等级的判断正确率分别为100%、93.5%、100%，平均为97.8%，交互验证对Ⅰ～Ⅲ健康等级的判断正确率分别为100%、90.3%、92.6%，平均为94.3%。综合两类验证方法，平均判断正确率为96.1%，说明所建立的判别函数是有效的，聚类分析结果是比较准确的，研究所采用的模型是可行的。

表10 标准化的Fisher线性判别式函数系数Table10 Standardized Fisher linear discriminant function coefficients

表11 不同健康等级马尾松人工林的验证结果Table11 Verification results of different health grades of Pinus massoniana plantation

4 结论与讨论

森林健康评价是我们了解森林健康状态的一种重要技术手段，做好森林健康评价可为森林健康经营及林业发展提供理论依据。在评价过程中建立一个准确的指标体系是评价森林健康的首要任务，而指标的选取和评价方法的选择是关键。由于森林生态系统的复杂性，目前对森林健康评价指标体系的构建还没有明确统一的标准。例如，黄国胜[12]等从林分活力、林分结构、森林恢复力抵抗力3 个方面选取树种多样性、群落结构、自然度、林木枯损率、郁闭度、灾害等级等9 个指标，利用层次分析法在林分尺度上对三峡库区森林健康进行定量评价；楚春晖[23]等选取平均胸径、平均树高、群落结构、郁闭度、每亩蓄积量、林分起源等12 个指标，利用SOM 神经网络对五指山市的森林健康状况进行评价。在上述健康评价中，存在没有考虑森林经营目的的情况，不同林种的森林不能用同一个健康状态标准进行评价，应该具体问题具体分析，且关于森林健康评价的研究大多只是评价出健康结果，并未对评价结果的正确性和合理性进行验证。尽管目前有很多种数学方法[10,24-27]应用于森林健康评价中，但是这些方法都或多或少存在缺陷，例如评价指标之间的相关性和指标权重的主观性，因此将多种数学方法结合起来对森林健康进行定量研究不失为一种有意义的尝试。

本研究在热林中心系统抽样样地中选取84 块马尾松人工林样地作为研究对象，使用因子分析法构建了以林分为尺度的马尾松人工林健康评价指标体系，选取林分平均胸径、平均树高、每公顷蓄积、地位指数、密度、郁闭度、树种组成指数、灌草盖度、腐质层厚度、土壤厚度、土壤有机质含量、灾害程度、天然更新等级13 个代表性评价指标，涵盖了森林生产力、土壤条件、生态系统稳定性3 个评价方面，较全面地反映了森林生态系统的真实情况。通过聚类分析将马尾松人工林分为3 个健康等级，聚类结果表明，研究区马尾松人工林中处于健康状态的林分占32.1%，亚健康占35.8%，不健康占32.1%，且不同龄组马尾松人工林健康状况从优到劣的排序依次为：近熟林＞中龄林＞成熟林＞幼龄林。根据聚类分析结果采用Fisher 线性判别进行判别分析，并进行自身验证和交叉验证，平均判断正确率为96.1%，具有较高的判别精度，说明所建立的判别函数是有效的，聚类分析结果是比较准确的，研究所采用的模型是可行的。

通过研究结果可以看出研究区大部分马尾松林分还处于非健康状态，近熟林健康状况优于中龄林，中龄林优于成熟林，成熟林优于幼龄林。所以建议对近熟林和中龄林进行疏伐、间伐和卫生伐，伐掉林内的病弱植株，调整林分密度，促进林木生长。对于成熟林应进行择伐或主伐，实现森林更新。对幼龄林要加大抚育力度，一方面对林木适时采取修枝、割灌等人为干预措施，促进天然更新与生长；另一方面，实施松土、灌溉、施肥等林地抚育措施，以提高林地生产力。