基于最大胸径生长率的浙江省杉木人工林立地质量评价

杨海宾， 张茂震， 丁丽霞， 于晓辉

（ 1. 浙江农林大学 省部共建亚热带森林培育国家重点实验室， 浙江 杭州311300； 2. 浙江农林大学 环境与资源学院， 浙江 杭州311300）

森林立地质量评价是实现科学造林、 合理高效利用林地的重要保证。 掌握森林立地质量可为森林的可持续发展制定更科学的经营措施及提高森林的经济、 生态、 社会效益提供理论基础［1］。 传统的森林立地质量评价采用的是地位指数或地位级作为重要指标， 即根据林分优势木的平均高和林分平均高2 种直接评价方法评定立地质量［1-3］。 虽然这2 种方法在同龄林的质量评价中精度较高， 但需要进行大量的实地调查、 采伐并制作解析木， 导致经济成本和生态成本较高， 永久性样地和代表性试验样地缺乏， 因而传统森林立地质量评价所需的大尺度数据相对缺乏［3-4］， 实用性较差。 随着近年来遥感技术应用于立地质量的评价研究［5-6］， 大尺度的立地分类逐渐克服了传统森林立地质量评价中生境制图时的耗时、 耗人力的缺点。 此外， 统计学方法用于森林立地质量评价和分类的实践也取得了良好效果， 如聚类分析、 多元回归分析、 数量化理论［7-12］等各种统计方法的应用使得定量研究在立地分类和立地质量评价综合系统的构建中有了更广泛的空间［9-14］。 然而， 尽管这些方法可以增加森林立地质量评价和分类结果的经济产出， 并在一定条件下提高精度， 但他们基本上仍沿用传统的评价指标， 在成本与效益之间没有突出改善， 难以适应大范围森林立地质量评价［15-17］。 国家森林资源连续清查数据是覆盖范围广、 精度可靠的数据， 但在较大范围的森林立地评价研究应用却非常少， 相关案例主要是结合遥感影像以样地总生物量为评价指标。 而样地总生物量会受年龄、 人为活动等多种干扰因素的影响， 难以准确评价森林立地质量。本研究利用近期5 期浙江省国家森林资源清查（national forest inventory， NFI）样地数据， 提取杉木Cunninghamia lanceolata林样地优势木的最大胸径生长率， 作为备选立地质量评价指标， 结合数字高程模型（DEM）和土壤数据， 运用数量化理论Ⅰ方法建立浙江省杉木人工林立地质量评价数量化理论模型。

1 研究区概况与研究数据

1.1 研究区概况

浙江省地处中国东南沿海长江三角洲南翼， 属亚热带季风气候， 季风显著， 四季分明， 年气温适中， 光照较多， 雨量丰沛， 空气湿润， 雨热季节变化同步， 气候资源类型多样， 气象灾害繁多。 年平均气温15.0～18.0 ℃， 年平均降水量为980～2 000 mm， 年平均日照时数1 710～2 100 h。

全省陆域面积为1 055×104hm2， 其中山地和丘陵占74.6%。 全省林地面积667.97×104hm2， 森林面积584.42×104hm2（其中85.95×104hm2为国家特别规定灌木林）； 杉木人工林纯林面积96.74×104hm2， 占森林面积的19.43%。 全省大部分地区适宜杉木生长， 南部属于杉木中心产区。

1.2 研究数据

1.2.1 立地因子数据 立地因子（包括地形地貌、 土壤、 气候、 植被、 生物等）对树木的生长有重要影响。 坡度、 坡向、 海拔、 土壤、 坡位、 湿度、 太阳辐射等地形及环境因素会不同程度影响杉木的生长［18-19］，其中地形因子主要包括海拔、 坡位、 坡向、 坡度等。 湿度及太阳辐射由于数据较难获取［20］， 本研究考虑其与地区所在的纬度关系较大， 因此将纬度列为立地因子之一。 地形因子可由DEM 处理得到， 浙江省DEM 由地理空间数据云网站上下载得到。 DEM 提取纬度、 海拔、 坡位、 坡向、 坡度等信息在Arc GIS软件中完成。

1.2.2 样地数据 浙江省1989-2009 年共5 期NFI 样地数据， 包括样地和样木数据。 由于浙江省在1989 年后有部分样地位置移动， 1989 年清查数据仅作分析参考用， 实际应用为1994-2009 年共4 期数据。 4 期NFI 样地数据中， 有连续复位样木的杉木人工林样地529 个， 用于建模； 另有用于模型检验的临安地区2004-2009 年杉木林固定样地51 个。 所有样地数据包含样地编号、 样地位置、 行政区、 地理地貌信息、 土壤信息和样地内植被信息等。 样木数据包含样地编号、 样木编号、 树种、 胸径等。 由于样地是固定（永久）的， 样地号和样木号在不同调查期是固定的， 因而可实现每株样木生长的连续监测。

2 研究方法

2.1 基于数量化理论Ⅰ的立地质量评价模型

数量化理论Ⅰ是一种类似多元回归的分析方法［15-16］， 用于自变量都是定性变量、 基准变量是定量变量的因素分析与预测问题， 采用说明性多变量模拟线性表示式中基准变量的定量变化。 假定有如下的线性模型：

假设自变量中有h个是定量的， 它们在第i个样本中的数据为xi（u）（u=1, 2,…,n）， 有m个是定性的，即m个项目， 其中第j项目有rj个类目， 它们在第i个样本中的反应式基准变量的数据为yi（i=1,2, …,n）。

利用最小二乘法， 即Bu（u=1,2, …,h）和bjk（j=1,2, …,m;k=1,2, …,rj）的最小二乘估计构造正规方程：

求解正规方程组（2）得预测方程, 并对方程精度进行检验， 同时通过偏相关系数、 方差比和范围评价每个立地因子对因变量作用的大小。

数量化理论Ⅰ相对于一般回归分析， 其不同之处在于可把定性变量与定量变量同时纳于回归模型。在数量化理论模型中， 因变量和自变量被分别称为基准变量和说明变量。 说明变量中定性变量称为项目， 每个定性变量的划分称为类目， 定性变量的类目类似于定量变量的取值区间。

2.2 立地因子的量化

根据数据易采集、 便于立地质量评价模型推广的原则， 参考前人研究成果， 选择立地因子。 本研究选取影响杉木生长的地形因子与土壤类型作为立地因子。 地形因子包括地貌、 纬度、 海拔、 坡度、 坡向、 坡位等， 其中纬度、 海拔、 坡度是定量变量， 地貌、 坡向、 坡位与土壤类型是定性变量。 根据数量化理论Ⅰ对变量的约定， 结合样地立地因子数据， 对定性立地因子进行类目划分（表1）， 并对项目（因子）的类目利用示性函数（0, 1）数量化方法对定性因子进行定量化处理， 即：

将每个样地的定性因子代入式（3）中， 与定量因子一起编制立地因子数量化反应表（表2）， 用于立地质量模型建立与评价。

2.3 立地质量评价指标确立

反映森林立地质量最主要和直接的指标是生产力。 可观测的森林立地生产力指标主要有特定年龄的树高、 胸径和蓄积量， 或者其生长量， 应用最多的是特定年龄的优势木树高。 然而， 优势木树高和年龄的准确测量不仅成本很高， 而且难度特大。 胸径生长率是最能直观反映立地质量的森林立地质量评价因子。 现有的NFI 数据信息量丰富， 且数据可靠性好。 胸径是NFI 调查中最为重要的测树因子， 但在NFI中不调查单株样木的年龄。 从NFI 中提取的胸径生长率只能代表在某一清查期内某棵树胸径的年均生长率， 并不能直接反映所在地段的立地质量好坏。 为了充分利用NFI 数据开展森林立地质量评价， 本研究提出1 个新的立地质量评价指标： 最大胸径生长率指标。 以最大胸径生长率指标为立地评价指标， 与地形地貌、 环境、 土壤等易获取因子建立模型， 可以用于估计包含无林地在内的所有林地的潜在生产力。一般来说， 林木的生长规律是从幼龄林到成熟、 过熟林， 胸径和树高的生长速度都呈逐渐降低的趋势，即在树木的生长过程中生长率都会出现一个高峰。 因此， 在NFI 清查次数足够多的情况下,每个样地复位样木的生长率都会在某一时期出现最高值， 即最大胸径生长率， 它代表了该地块的最大生产潜力。 若以R表示最大胸径生长率， 则：

式（4）中：k为清查期，rk为第k期的胸径生长率，n为森林资源清查总期数。

表1 定性立地因子类目划分与量化编号Table 1 Classification and quantitative numbering of qualitative site factors

表2 立地因子反应表Table 2 Response table of site factors

NFI 清查间隔期为5 a， 杉木的龄级期也是5 a， 5 期NFI 清查数据即可覆盖5 个龄级期的生长率。假定杉木轮伐期为25 a， 近熟林进入成熟林时全部采伐， 则5 次清查能包含全部最大胸径生长率。 因此， NFI 与龄级之间的关系如图1 所示。 图1 中， 样地样木的龄级为Ⅰ～Ⅴ， 成熟林全部采伐（由于现实森林的林龄不同， 因而有A、 B、 C、 D、 E 等5 种情形）。

样木胸径生长率的计算方法如下：

式（5）中：rijk为第i个样地、 第j株复位样木、 第k期清查的胸径生长率；Dij（k+1）为k+1 期样木胸径；Dijk为第k期样木胸径。

样地胸径生长率计算建立在样木生长率的基础上， 设第i个样地、 第k期清查的胸径生长率为rik， 则：

样地优势木生长率比样地生长率更能代表样地所在位置的立地生产潜力， 样地优势木胸径生长率计算公式如下：

式（7）中： rikD即第i样地第k期3 株优势木生长率的均值。

根据样地的胸径生长率， 从每个样地多期调查数据所生成的胸径生长率时间系列中选取最大者：

式（8）中：riD相当于样地i位置幼龄林优势木的平均胸径生长率。

图1 清查期次与样木龄级分布Figure 1 Inventory period and age class distribution of sample trees

2.4 立地质量评价精度检验

2.4.1 模型显著性检验 复相关系数是估计回归精度的重要指标之一， 复相关系数越大， 则估测相关越密切， 说明模型对优势木最大胸径生长率均值的估计效果越好， 所得的数量化得分表就越可靠［21］。 根据复相关系数检验线性关系是否显著进行F检验， 可确定坡位、 坡向、 土壤类型、 纬度位置、 海拔、 坡度等变量的线性关系是否显著［22］。

2.4.2 模型预测精度检验 模型预测精度检验采用均方根误差（RMSE）作为模型的评价指标。 均方根误差是用来衡量实际测量值同模型预测值之间的偏差， 它可作为衡量测量精度的一种数值指标［23］。 在模型建立后， 为对模型拟合性能进行检验， 运用浙江省杭州市临安区的实际测量值代入模型计算最大胸径生长率进行实地检验， 采用平均绝对误差、 RMSE 值进行评估， 反映模型估测的精度。

2.4.3 方差比检验模型检验 采用方差比方法进行。 该检验方法具有精确的样本分布， 而不依赖于大样本渐近极限分布； 在通常的时间序列数据非正态分布的情况下， 这种非参数方差比检验具有更高的检验功效。 方差比检验对于一段时间内某一时间间隔的方差比例效果较好。 方差比检验是指各因子（项目）方差占模型方差的百分比， 表示立地因子（坡位、 坡向、 土壤类型、 纬度位置、 海拔、 坡度）对因变量（最大胸径生长率）的作用是否显著。

3 结果与分析

图2 杉木单株生长率比较Figure 2 Comparison of growth rate of C. lanceolata

3.1 基于NFI 的立地质量评价指标

3.1.1 单株胸径生长率与胸径的关系 根据NFI 数据分析， 杉木生长率从幼林到成熟林呈下降趋势。 图2A 和图2B 显示了1994-1999 年和2004-2009 年调查间隔期内杉木单株生长率按胸径分布的变化。 首先， 胸径生长率随胸径增大而降低， 2 个时期趋势相同。 其次， 其趋势为单调下降。 再次， 曲线在纵轴上分布较宽， 特别是当胸径大约为20～30 cm 时， 生长率出现较大的变化范围。 图2B（2004-2009 年）最为明显。 它们代表了多条单调下降的曲线， 其本质是不同立地质量的生长率差异。 这一特点图2B 较图2A 突出， 原因是2004 年这有部分位置移动后的新样地加入。 2 幅图同时显示， 随着胸径增大， 其生长率分化明显， 而在较小胸径阶段生长率随胸径增大而较快下降， 且在纵轴上分布区间较窄。

3.1.2 单株胸径生长率随时间变化趋势 单株胸径生长率连续变化趋势与年龄有关。 在无年龄的情况下， 这种变化在连续多期之间呈单调曲线形式上升或下降。 图3 为单株木胸径生长率随胸径增长而的变化情况。 为了图示清晰，将每个胸径值所对应的多株样木胸径生长率取均值， 每条曲线代表1 个胸径级。 图3 中D8 代表第1 次清查时胸径为8 cm， D12 代表第1 次清查时胸径为12 cm， ……。 余类推。 从图3 可以发现： 小径级样木胸径生长率大， 大径级样木胸径生长率小。 以后一直保持这一规律。 对于同一株样木来说， 清查次数增加就是年龄的增加， 随着年龄的增长， 样木胸径生长率下降， 最后趋于平缓。 结果表明： 在杉木生长过程中， 单株胸径生长率存在最大值， 这个最大值代表所在地段的森林立地生产力。 因此， 最大胸径生长率可以被用作森林立地质量评价的直接指标。

3.1.3 样地最大胸径生长率 通过多期复位样木分析， 提取每个样地的最大胸径生长率。 由于样地采伐、 灾害等因素影响， 用样地全部样木胸径生长率的平均值代表样地生产潜力的状态不适合所有样地。 表3 显示： 用3 株最大胸径生长率取平均与样地的立地因子有更好的相关性。 因此， 为了避免偶然因素的干扰， 每个样地的最大胸径生长率取3 棵优势木胸径生长率的均值。共提取529 个样地的最大胸径生长率。

图3 杉木不同胸径径级平均生长率随时间变化趋势Figure 3 Trend of average growth rate of different DBH grades of C. lanceolata with time

表3 最大胸径生长率与立地因子之间皮尔森相关系数Table 3 Pearson correlation coefficient between the maximum DBH growth rate and the site factor

3.2 杉木林立地质量评价模型

以优势木最大胸径生长率均值作为因变量， 坡位、 坡向、 土壤类型、 纬度位置、 海拔、 坡度作为自变量， 采用数量化理论Ⅰ建立杉木林立地质量评价模型。 4 期529 个NFI 的样地数据参与数量化拟合，运用软件SPSS 19.0 完成线性回归分析模型建立与评价。 具体结果见表4。

根据表4 数据， 本研究最终建立浙江省杉木林立地质量评价预测模型公式为：

式（9）中： δ（1）为样地纬度位置， δ（2）为样地海拔， δ（3）为样地坡度， δ（4,i）为样地坡位第i个项目， δ（5,i）为样地坡向第i个项目，δ（6,i）为样地土壤类型第i个项目。 表4 显示模型统计检验F=14.723， 达极显著水平， 模型复相关系数R为0.606。 方差分析（表5）表明： 各项目对模型贡献均达显著水平。

3.3 模型检验

3.3.1 方差比检验 方差比是各因子（项目）方差占模型方差的百分比， 表示因子对因变量（最大胸径生长率）的作用是否显著。 按公式计算各影响因素的方差， 结果见表5。 结果表明： 对胸径生长率的贡献程度从大到小依次是坡向、 坡位、 土壤类型、 坡度、 纬度、 海拔。 其中坡向对杉木林胸径生长率的贡献率最大， 坡向会影响太阳光照量和太阳辐射强度， 造成不同坡向的日照强度和热量不同， 从而对杉木林的生长产生较大影响。 坡位对杉木林的胸径生长率的贡献率较大， 不同坡位具有相应的小气候环境， 影响着湿度、 水分和温度等条件。 土壤类型对胸径生长率的贡献同样较大的原因是由于土壤类型不同， 杉木生长的土壤营养环境不同， 从而对杉木林的生长产生一定影响。 从表4 和表5 可以看出： 纬度和海拔的影响较小， 主要是因为本研究所选择的研究区为浙江省， 纬度变化范围较小； 同时浙江省内山地平均海拔一般在200～1 000 m， 而本研究杉木人工林样地分布在400～800 m 较窄的范围内， 两者的差距对于杉木林来说不能形成较大区分度， 因此纬度和海拔的对于胸径生长率的贡献率较小。

3.3.2 基于临安实地数据的模型拟合精度分析 在模型建立后， 为对模型拟合性能进行检验， 运用浙江省杭州市临安区的实际测量值代入模型计算最大胸径生长率进行实地检验， 采用平均绝对误差、 RMSE值进行评估， 对全部样本和临安地区样本的观测值和预测值之间的误差进行评价， 反映模型估测的精度。 从表6 可以看出： 临安地区样本和全部样地作为模型估计结果检验的平均绝对误差分别为3.407 和2.125， 整体平均绝对误差相对较小， 而且随着样本的增加， 模型的拟合精度提高， 表明利用最大生长率进行数量化理论拟合是精度相对较高的方法。 均方根误差能够反映预测值和测量值之间的偏差程度，数据显示临安区样本和全部样本检验的均方根误差为0.518 和1.421， 而计算两者的相对均方根误差的值分别为0.073 和0.062， 相对均方根误差越小， 相对均值的离散程度越小， 模型的精度也就越高。

表4 杉木林立地因子数量化理论Ⅰ得分及t 检验表Table 4 Score and t-test table of the quantification theory Ⅰmodel of C. lanceolata

表5 各类立地因子影响因素方差比Table 5 Variance ratio of the influencing factors of various site factors

表6 模型拟合误差分析表

4 结论与讨论

基于NFI 的最大胸径生长率是杉木人工林立地质量评价的可靠指标， 在经济上节约、 高效。 从建模因子的贡献率来看， 对胸径生长率的贡献程度从大到小依次是坡向、 坡位、 土壤类型、 坡度、 纬度、 海拔。 利用最大胸径生长率， 结合DEM 数据建立数量化理论Ⅰ模型， 模型拟合效果和检验结果均比较理想。 在缺乏树龄和树高数据的情况下， 可以利用基于多期NFI 数据的样地最大胸径生长率作为评价指标进行森林立地质量评价。 不仅可以用于现有杉木人工林立地质量评价， 还可以用于无林地立地生产潜力预测。

本研究提出基于NFI 的最大胸径生长率作为杉木人工林立地质量评价指标， 结合地理信息建立立地质量评价模型， 结果符合研究区内杉木林生长立地条件现状。 但限于研究区NFI 数据积累的时间跨度有限， 模型精度还有进一步提高的空间。 由于NFI 每5 a 清查1 次， 其数据的时间分辨率为5 a。 因此，NFI 清查期数越多， 在同一个样地所测样木的年龄阶段越多， 从中提取的最大胸径生长率对相应立地生产力的代表性越接近真实。 随着NFI 数据的不断积累， 基于NFI 的森林立地质量评价精度会不断提高。同时， 基于数量化理论Ⅰ的杉木最大胸径生长率估计模型所解决的问题是得到所有森林立地上的假设最大胸径生长率， 只是一个以生长率代表立地质量的初步展现， 在此基础上还可以有进一步的分类、 聚类及其检验等分析。