基于立地混合效应的湖南丘陵平原区杉木多形立地指数模型研究

宋 争，黄 朗，胡 松，张雄清，朱光玉*

(1. 中国林业科学研究院林业科技信息研究所，北京 100091；2. 中南林业科技大学林学院，湖南 长沙 410004；3. 中国林业科学研究院林业研究所，北京 100091)

立地质量评价是对森林立地的宜林性或潜在生产力判断或预测的基础工作[1]，立地指数是一种最普遍的评价森林立地质量的方法[2-3]，分析某树种林分基准年龄时对应的林分优势高可用于分析森林立地生产潜力的基础[4]。除了根据优势高与年龄指标，还能够通过样地实测值和收获表蓄积量估算获得立地指数[5]。但由于立地条件的差异性或多样性，基于优势高-年龄的导向曲线模拟精度往往较低，同一树种在不同区域的立地指数估计结果随立地质量的不同而存在显著差异。因此，很多学者在评价立地质量时，优先分析环境要素（气候、植被、地形等）与立地指数的关系，得到包含环境因子的立地指数函数关系[6-7]。也有学者采用树高-年龄分级哑变量的方法，用样地林木平均高划分等级代替立地因子评价立地质量[8]。这些方法虽然提高了立地指数建模精度，但关于环境因子在模型中的表达仍需进一步探索。

目前，如何建立立地指数模型成为描述森林立地质量的重点和中心问题。研究指出有单形和多形两种关系，其中多形立地指数曲线结合优势木平均高和其他反映立地条件因子，共同评定立地质量[9-11]。最常用的2种拟合方法为固定基准年龄的静态方程和可变基准年龄的动态方程，前者依赖于选定的基准年龄，后者局限于严格的数据要求[12]。而混合效应模型是一种最常用的处理局部变异性的建模方法[13-14]。对于不同立地类型的优势高生长，混合模型不仅能够描述整体变化的平均趋势，还能够给出协方差和方差等相关数据以了解个体特点[15-16]。虽然在构建立地指数模型时，相关文献运用数量化方法I[17-18]，然而该方法前提是固定效应，如果考虑随机效应，该方法不适用。

鉴于此，本研究以湖南丘陵平原区杉木（Cunninghamia lanceolata(Lamb.) Hook.）人工林为研究对象，首先利用数量化方法I确定影响杉木优势高生长的显著性影响因子；然后以主导立地因子及其组合作随机效应，建立基于立地随机效应的杉木多形立地指数曲线模型，通过模型选优对不同立地因子水平的立地质量进行评价；最后以Kmeans聚类划分立地类型组，解决复杂立地类型的模型应用问题，以提高立地指数曲线模型的适用性和准确性。

1 材料与方法

1.1 研究区概况

湖南位于中亚热带，地貌大部分为丘陵，呈现大陆性季风湿润气候，全年平均气温17.5 ℃左右，平均降水量1200～1700 mm，日照时长1300～1800 h，雨水丰富且热量充沛。研究区乔木树种以杉木（Cunninghamia lanceolata(Lamb.)Hook）为主，其次为马尾松（Pinus massonianaLamb.）与拟赤杨（Alniphyllum fortunei(Hemsl.)Makino）等。主要小乔木与灌木有冬青（Ilex chinensisSims）、木姜子（Litsea euosmaW.W.Smith）、山茶（Camellia japonicaL.）、杜茎山（Maesa japonica(Thunb.) Moritzi.）等。主要草本植物有狗脊蕨（Woodwardia japonica(L. f.) Sm.）、铁芒萁（Dicranopteris linearis(Burm.) Underw.）等。

1.2 数据来源

以湖南丘陵平原区杉木人工纯林为研究对象，在湘乡、怀化、会同、株洲、永州、临湘、绥宁、桃源等26个县，设置样地（20 m × 30 m）共360块，记录相关立地因子，并对样地内活立木每木检尺。样地内每0.01 hm2选取最粗优势木1株，计算优势木平均高(H)。相关林分调查因子统计详见表1。

表1 杉木林分调查因子统计Table 1 Statistics of stand factors for Cunninghamia lanceolata

1.3 研究方法

1.3.1 立地分类方法 选取如表2所示的影响杉木优势高生长的6个立地因子，根据《中国森林立地》相关标准进行分级，将海拔每100 m划分为一级，分析不同海拔的差异。

表2 立地因子等级划分Table 2 The division of site factor grades

立地因子作自变量，各林分优势木平均高做因变量，运用数量化方法Ⅰ对杉木立地因子分别分析。基于方差分析中不同“Pr>F”大小，筛选显著因子并进一步得到主导因子。

1.3.2 非线性混合效应模型 为更准确地评价杉木立地质量，本研究考虑立地因子作随机效应对杉木立地指数模型的影响，随机效应包括了立地因子与其交互作用。非线性混合效应模型的形式为[19-20]：

式中：Hij为 第i个立地类型下第j个样地的平均优势高，Tij为第i个 立地类型下第j个 样地的林龄，φij为参数向量，εij是误差项，M是立地类型数量，ni为第i个立地类型的样地数，β为（p×1）维固定效应向量，bi为（q×1） 维随机效应向量，Aij,Bij为设计矩阵，D为随机效应的方差—协方差矩阵，Ri为立地类型内的方差—协方差结构，σ2为方差，Gi为描述方差异质性的对角矩阵，Γi为描述随机效应自相关性的方差矩阵。

构建非线性混合效应模型，需要对固定与随机效应参数进行构造。通常情况下对所有可能的参数考虑随机效应，但这会造成模型不收敛[21-22]。为避免模型不收敛问题，本研究选择在1个参数上考虑立地的随机效应。其次，由于本研究中林分优势高与年龄不是重复调查数据，不存在样地内测量值之间的相关性。因此，假设误差项εij服从正态分布。

最后，本模型只考虑单个随机效应参数来确定随机效应参数方差，因此假定随机参数构造类型为无结构，并采用AIC、BIC、对数似然比（Loglikelihood）和确定系数R2选择最优混合效应模型。AIC、BIC值越小，Log-likelihood和R2越大，模型的拟合效果越好。

1.4 候选基础模型

基础模型不仅要能够很好描述林木优势高生长规律，还要具有生物学的可解释性[23-24]。其中立地指数曲线基础模型常为“S”形，本研究选择常用的立地指数曲线方程（M1～M8），拟合杉木林分优势高（见表3）。

表3 候选基础模型Table 3 Candidate basic models

1.5 模型精度评价

利用确定系数（R2） 、预估精度（P）、平均绝对误差（MAE）和均方根误差（RMSE）进行模型评价。

式中：yi为 第i个 样地林分优势高实测值，为第i个样地林分优势高估计值，n为样本数，tα为置信水平α=0.05时的t分布值，p为模型参数个数，为平均预估值。

2 结果与分析

2.1 因子显著性分析

林分优势木平均高受多种立地因子综合作用影响。根据优势木平均高与不同立地因子的相互关系，采用数量化方法Ⅰ对其进行显著性检验，结果（表4）表明：对优势木平均高具有显著影响的立地因子有海拔、坡度、坡向与土壤类型，其显著性顺序为土壤类型＞海拔＞坡向＞坡度。

表4 立地因子的显著性检验Table 4 Significance testing of site factors

2.2 基础模型拟合与选择

利用 Forstat2.1软件对候选基础模型M1～M8进行估计，相关基础模型参数估计值、确定系数（R2） 、预估精度（P）、平均绝对误差（MAE）及均方根误差（RMSE）如表5所示。

表5 候选基础模型拟合结果Table 5 The fitting results of candidate basic models

从8个候选基础模型（表5）拟合结果可知，这8个模型的拟合精度均较低（R2=0.4243～0.5644），也证明了区域性立地指数建模精度低的问题。其中，模型M4拟合效果最好，其确定系数（R2= 0.5644）与预估精度（P=0.9750）最大，平均绝对误差（MAE=6.7626）及均方根误差（RMSE=2.4100）最小。

2.3 非线性混合效应模型模拟

以主导立地因子作随机效应，考虑立地单因素与多因素交互作用，构建立地指数混合模型。结合基础模型参数个数和随机效应因子组合类型，拟合所有参数及其组合形式，最优随机效应立地指数模型通过AIC和BIC等评价指标获取（表6）。

计算结果显示，除HB、PX、TL等单因子在其它参数组合的拟合效果较优外，受水平数影响，其它随机效应组合均只在渐进参数a上时，模型收敛。因此，含立地随机效应的杉木多形立地指数模型表达式为：

式中：a，b，c为模型固定效应参数；μi为随机效应参数；H为林分平均优势木高；Age为林分年龄。

拟合多种随机效应因子组合类型后的结果如表6所示，AIC、BIC相对于基础模型M4都显著下降，Log-likelihood和R2均大幅提高，这说明基于立地随机效应的混合模型明显优于基础模型。此外，混合模型拟合效果的好坏，取决于立地因子显著性的高低，立地因子越显著，相应随机效应组合类型的混合模型拟合结果越好。其中，海拔、坡度、坡向、土壤类型4个主导立地因子交互作用的混合模型AIC、BIC最低，Log-likelihood及R2最高。根据立地分类的原则与方法，4个主导立地因子排列组合可得到立地类型，综上所述含立地类型的混合模型为最优模型。

表6 基于立地效应的混合模型拟合结果Table 6 The fitting results of mixed effects models considering site effect

2.4 模型应用

依据海拔、坡度、坡向、土壤类型4个主导因子，将研究区划分为78个立地类型，那么相应的混合模型随机效应参数有78个，势必会影响立地指数模型的应用。为了简化立地类型，将初始立地类型（ST）进行K-means聚类。参考森林立地分类的标准，确定系数≥0.99为聚类精度标准，合并得分值相近的立地类型为立地类型组（GST），使得本研究区的立地类型聚类为11种（表7）。

表7 立地类型聚类结果Table 7 Classification of cluster with site types

将聚类后的立地类型组（GST）作为随机效应构建立地指数非线性混合模型，利用AIC、BIC、Log-likelihood及R2等4个指标进行模型选择（表8），并与基础模型（None）、初始立地类型（ST）模拟结果进行比较。如表8所示，考虑立地的随机效应后，模型确定系数从0.5644提高到0.8089～0.8117，其中基于立地类型组的混合效应模型模拟精度最高（R2=0.8117，Loglikelihood=-704.9879）。因此，基于立地随机效应的立地指数曲线模型可以极大提高模型精度，同时利用K-means聚类划分立地类型组的方法，不仅保证了预估精度，而且还简化了立地类型，提高立地指数模型的实用性。

表8 模型精度评价Table 8 Model evaluation

利用该模型绘制湖南丘陵平原区杉木多形立地指数曲线（图1）。图中展示了不同立地类型组的立地指数曲线图，各立地类型组的模型模拟效果均较好。

图1 多形立地指数曲线Fig. 1 Polymorphic site index curve

3 讨论

林分优势高生长受林分密度影响很小，是同龄林立地质量评价采用最广泛的指标，但关于林分优势高的定义及测定标准并不一致[25]。赵美丽等[26]对樟子松（Pinus elliottiiEngelm.）同龄林最高、最粗木平均优势高的相关性、稳定性研究后认为，最粗6株木法能科学得到林分优势高，该方法具有简易、稳定的优点，所以本研究以每100 m2选择1株最粗木的方法计算平均优势高。

在森林立地分类与评价的因子中，地形属间接因素，但综合反映着气候和土壤特性。本研究在进行显著性研究中用数量化方法I，以分析不同立地因子对杉木优势高生长的影响，结果表明，海拔、土壤类型、坡向、坡度与杉木优势高生长有强相关性。不同立地因子对优势高生长的影响也不同。产生显著性差异的原因有以下几点：从海拔来说，南方丘陵平原地区热量和降水量的垂直分异总的来说决定了研究对象的生长与分布，海拔为一个主要因素。从土壤类型与坡度来说，杉木对土壤湿度、养分、通气性等要求较高，杉木生长由土壤类型的差异直接影响，坡度间接表现土壤特性。从坡向来说，不同坡向热量光照差异，对杉木不足以引起其生长的很大差异，故坡向相比较而言，显著性影响要低于土壤类型、海拔，但高于坡度的影响作用。王冬至等[27]在研究塞罕坝华北落叶松（Larix principisrupprechtiiMayr）人工林地位指数模型时，将林分优势高与不同立地因子进行相关性分析，发现林分优势高与立地因子显著相关。

基于林分优势高-年龄关系的立地指数，已成为表征森林立地生产力的一个重要参数，应用该方法的前提在于导向曲线具有较高的模拟精度。然而对于某一特定的树种，导向曲线的形状会因气候、土壤、地形等的不同而有所变化，使用同一导向曲线模拟树高生长必然会造成较大误差[28]。而混合模型在估计总体的平均效应和个体差异上具有独特优势，可以通过多参数方法构建立地指数模型[29]，有利于增加模型的评估精度。Wang等[30]利用哑变量与混合模型法分别模拟火炬松（Pinus taedaL.）树高生长，最终表明混合模型方法准确性更好。因此，为了解释不同立地水平的优势高生长差异，本研究采用单水平的非线性混合效应模型进行立地随机效应模拟，确定系数（R2）从基础模型拟合效果的0.4243～0.5644提高到0.5655～0.8089，以立地因子及其组合作随机效应的混合模型对比基础模型优势显著。

本研究以K-means聚类为基础，将复杂的立地类型（ST）聚类划分成立地类型组（GST），并在立地指数模型中考虑GST为随机效应，构建多形立地指数模型（R2=0.8117）。此方法较原来基础模型极大提高了模型预估精度，并且验证了立地类型聚类构造随机效应的方法具可行性。

4 结论

本研究采用非线性混合效应模型方法，构建了含立地随机效应的湖南丘陵平原区杉木多形立地指数曲线模型，利用Forstat求解模型参数，得出如下结论：海拔、坡度、坡向与土壤类型是影响湖南丘陵平原地区杉木优势高生长的显著性因子，以显著性因子及其组合作随机效应可以显著提高立地指数建模精度，且包含立地类型的混合模型拟合精度最高，拟合精度的高低与主导立地因子的显著性紧密相关；以K-means聚类将杉木78个初始立地类型划分为11个立地类型组，包含立地类型组的混合模型在保留预估精度的前提下提高模型实用性。因此，基于立地随机效应的多形立地指数模型可以显著提高建模精度，也为区域性复杂立地类型的立地质量评价提供了思路与方法。