河北省北部华北落叶松人工林立地指数空间分布预测*

李文博 吕振刚 黄选瑞 张志东

1.河北农业大学林学院 河北省林木种质资源与森林保护重点实验室 保定 071000； 2.河南省林业调查规划院 郑州 451200； 3.华中农业大学资源与环境学院 武汉 430070)

华北落叶松(Larixprincipis-rupprechtii)是河北省重要用材树种之一，主要分布于河北省北部地区，在维持该区木材生产、水土保持和森林游憩等方面发挥着重大作用(高敏等， 2017)。立地条件是林木生长发育和森林发展的重要基础，研究发现，气候、地形和土壤等立地条件会影响华北落叶松人工林生长及其适应性，进而影响立地生产力(王云霓等， 2013； 赵匡记， 2015)。当前，在林分水平关于华北落叶松人工林生产力的研究较多(王冬至等， 2015； 赵匡记， 2015)，但对于河北省北部地区华北落叶松人工林基于较大空间尺度立地生产力的确定及其影响因子的研究较少，不利于该地区华北落叶松人工林可持续经营(Bueisetal.， 2017)。

立地指数(site index, SI)指林分在基准年龄时的优势木高，其实质是将任意林分的树高最高值(林分优势木高)标准化至同一年龄的生长曲线(Littkeetal.， 2016)，因优势木高生长与材积生长密切相关且受林分结构和森林经营活动影响很小，因此SI被广泛用作同龄纯林的立地生产力评估指标。SI能够反映出区域土壤、气候等固有立地特征差异(McKenneyetal.， 2003； Skovsgaardetal.， 2008； Bueisetal.， 2017)，而气候、土壤和地形等因子的一些不利变化也会导致森林立地生产力下降并对森林生态系统造成较大威胁(Bravo-Oviedoetal.， 2011； Littkeetal.， 2016； Sabatiaetal.， 2014)，因此一些学者(Skovsgaardetal.， 2008； Noordermeeretal.， 2018)基于气候、地形和土壤等环境因子间接预测森林立地生产力。Nigh等(2004)建立不同气候变量组与SI的多元回归模型，得出不同树种SI对气候变量的响应程度不同； Bravo-Oviedo等(2011)通过气候、土壤因子与地中海油松(Pinuspinaster)SI的多元线性回归关系分析发现，在大的空间尺度，气候是该树种生产力的主要驱动力，但在同一气候区内，SI差异与土壤密切相关； Jiang等(2015)构建树种SI与气候和土壤不同组合的随机森林(random forest，RF)模型，认为基于气候和土壤的SI模型预测精度稍高于仅基于气候或土壤因子的SI模型。由此可知，SI预测精度受树种类型、研究区域和尺度、模型预测变量等因素影响，准确预测华北落叶松人工林SI及其空间变异规律，需综合分析SI与气候、地形和土壤等环境因子的关系。

近年来，随着林地需求的不断增加，预估森林立地生产力空间变异对森林经营和保护规划愈发重要(McKenneyetal.， 2003)。3S技术提高了较大区域尺度环境因子数据的精度和时效性，解决了传统调查成本高、样本获取量少的问题(Bravo-Oviedoetal.， 2011)，并能够绘制出信息明晰、分辨率较高的SI空间分布图，可为适地适树和制定精准提高森林质量的营、造林方案提供保证(Falkowskietal.， 2009； Parresoletal.， 2017)，因此很多研究基于3S多源空间环境因子数据构建了多元线性回归(multiple linear regression，MLR)(Huangetal.， 2017)、随机森林(RF)(Jiangetal.， 2015)、神经网络回归(Hlasnyetal.， 2017)等模型预测区域SI及其与环境因子的关系，并绘制了SI空间分布图。地统计学方法同样被用于森林立地生产力预测和绘图，如Fayad等(2016)基于机载和星载激光雷达提取的树高数据，构建与环境因子的回归克里格(regression Kriging, RK)模型预测法属圭亚那热带雨林高度并绘制了其空间分布图； Benitez等(2016)基于植被指数和土壤数据采用地理加权回归克里格(geographical weighted regression Kriging, GWRK)模型预测并绘制了亚马逊热带雨林地上生物量空间分布图。以往研究基于不同环境因子采用多种模型预测森林立地生产力的空间变异并对其结果进行可视化，但因林分类型、研究区域和尺度、模型预测变量等因素不同，模型预测精度存在较大差异。

为准确预测华北落叶松人工林SI及其空间变异规律，本研究基于河北省北部地区多源环境因子数据，筛选主要环境因子，应用MLR、RF、RK和GWRK模型预测研究区华北落叶松人工林SI并绘制其空间分布图，旨在解决3个问题： 1) 哪些环境因子是影响研究区华北落叶松人工林SI的主要因子； 2) 如何在MLR、RF、RK和GWRK 4种预测模型中选择最优预测模型，模型精度表现如何； 3) 华北落叶松人工林SI的空间分布格局及其对主要环境因子的响应机制是什么。

1 研究区概况

研究区位于河北省承德市围场满族蒙古族自治县(41°31′—42°37′N，116°41′—118°20′E)(图1)，北接内蒙古，地处内蒙古高原和冀北山地的过渡带，海拔663～2 042 m，分坝上和坝下2种自然地貌类型。坝上气候寒冷，年均温-1.2～4.9 ℃，年均降水量445 mm； 坝下四季分明，年均温-0.4～6.9 ℃，年均降水量441 mm，区域降水少且多集中于6—9月。受地形和气候影响，形成了褐土、棕壤、灰色森林土、风沙土和草甸土等土壤类型。研究区乔木主要有华北落叶松、樟子松(Pinussylvestrisvar.mongolica)、云杉(Piceaasperata)、蒙古栎(Quercusmongolica)等； 灌木主要有沙棘(Hippophaerhamnoides)、山刺玫(Rosadaverica)、稠李(Prunuspadus)等； 草本植物主要有蒲公英(Herbataraxaci)、苔草(Carextristachya)、老鹳草(Geaniumdaharicumvar.alpinum)等。

图1 研究区地理位置及SI样地点

2 数据收集与研究方法

2.1 数据收集

2.1.1 立地指数 采用塞罕坝机械林场和木兰国有林场管理局森林资源二类调查数据及近6年样地解析木数据。调查样木健康且均无损害、无虫害和人为活动干扰，其中，森林资源二类调查数据包含优势树种、林龄、胸径和树高等测树因子，研究主要选取优势树种为华北落叶松人工纯林的小班样地； 样地解析木采取5株优势木和5株平均木的方法，分别确定能够代表每块样地林分水平的1株优势木和1株平均木进行解析，按5年一龄阶、2 m一分段整化。由于不同立地林地生产力评价结果会因基准年龄误取而有所偏差(吴恒等， 2015)，故选取研究区内不同立地条件下20块华北落叶松人工近熟林(≥ 30年且年龄为5的倍数)样地的优势木解析资料分析树高生长过程。由图2可知，优势高生长在林龄20年后逐渐趋于稳定，30年时对应的样本量最大，故确定基准年龄为30年。由于森林资源调查树高为样地平均木高，因此筛选研究区内林龄≥ 20年的92块不同立地条件样地解析木数据建立华北落叶松人工林样地优势木-平均木树高转换方程(表1、图3)，将研究区内基准年龄为30年的华北落叶松人工林样地平均木高转换为立地指数，剔除±3倍标准差外的样地数据，最终对1 179块样地数据进行描述性统计分析(表1)。

2.1.2 环境因子 气候数据来源于世界气候数据库，空间分辨率30弧秒(约为1 km空间分辨率)，与研究区基础地理数据配准、裁剪后产生包括气温、降水等19个生物气候变量数据(Hijmansetal.， 2005； http:∥www.worldclim.org/)； 土壤数据主要来源于北京师范大学中国土壤数据集(Weietal.， 2013)，空间分辨率30弧秒。为探究土壤化学性质与立地指数的相关关系，采用研究区内土壤全氮、全磷等9个土壤因子栅格数据集； 从地理空间数据云(http:∥www.gscloud.cn/)下载研究区30 m分辨率DEM影像，影像预处理后再基于DEM生成坡向、坡度和海拔3个地形因子。

31个环境因子栅格数据均采用Asia_Lambert_Conformal_Conic投影坐标系，空间分辨率统一为30 m。为筛选出对SI影响大且因子间相关性低的变量，首先对31个环境因子和SI进行Pearson相关分析，得到显著影响(P< 0.05)SI的15个环境因子(表2)； 然后计算环境因子的方差膨胀系数VIF，剔除VIF > 15的环境因子，筛选得到海拔(DEM)、最干月降水(BIO14)、土壤有机质(SOM)、土壤全氮(TN)和土壤全磷(TP)5个环境因子； 最后建立SI与5个环境因子的逐步回归方程，进一步筛选出与SI显著相关的DEM、BIO14、TN和TP共4个环境因子，VIF分别为3.01、2.64、1.48和1.29，不存在多重共线性，4个环境因子的空间分布如图4所示。

图2 优势高生长变异系数

图3 样地优势木-平均木树高转换方程

表1 森林调查数据描述性统计结果

表2 环境因子描述性统计①

2.2 研究方法

2.2.1 随机森林(RF) RF是一种非参数决策树学习算法，每株决策树依赖于独立同分布的随机向量，最初是为分类和回归问题设计的，在机器学习和数据挖掘方面应用广泛(Breiman， 2001)。设响应变量(SI)为一维向量Y，环境因子为p维向量X，其中，Y可为连续、二进制、分类或生存变量集合，X可为连续或离散变量集合。本研究SI为连续变量，故实际上是在非参数回归框架下求X=x(X∈x)时的响应函数h(x)，此估算基于训练数据λ={(X1，Y1)，…，(Xn，Yn)}，其中，(Xi，Yi)与(X，Y)独立同分布，h(x)可表示为：

h(x)=E[Y∣X=x]，

(1)

即RF回归模型输出为所有决策树结果的平均值。模型基于R3.5.1调用“randomForest”包实现，决策树株数(ntree)和树节点(mtry)是模型的2个重要设置参数，针对RF解决回归问题，Breiman(2001)建议mtry默认值为全部响应变量的1/3(取整)，ntree也代表重采样参数，一般选取整体误差率趋于稳定时对应的值(欧强新等， 2018)。

2.2.2 回归克里格(RK) RK为通过建立响应变量与辅助变量之间的线性模型并对回归残差进行OK(ordinary Kriging)插值从而预测响应变量空间特征的统计模型(Hengletal.， 2007)。设样地为S，SI为响应变量，则：

SIRK(S)=SIMLR(S)+εRK(S)；

(2)

(3)

式中： SIMLR(S)和εRK(S)分别为模型趋势项和残差项，趋势项主要通过建立MLR模型获得，将MLR模型残差基于ArcGIS进行OK插值获得残差项；k为环境因子总数，k=4(图4)；bi为第i个环境因子Fi的模型系数；IMLR为截距。

2.2.3 地理加权回归克里格(GWRK) GWRK与RK模型的最大区别在于模型趋势项不同，GWRK首先采用局部回归方法构建模型趋势项，然后将局部回归残差进行OK插值(杨顺华等， 2015)，其表达式为：

SIGWRK(S)=SIGWR(S)+εGWRK(S)。

(4)

式中： SIGWR(S)为响应变量在样地S处的地理加权回归(geographical weighted regression，GWR)预测值。

构建GWR的目的是探究响应变量与辅助变量之间是否存在空间非平稳性关系，可由下式表示：

(5)

式中：Fi(S)为第i个环境因子在样地S处的实测值；βi(S)为第i个环境因子在样地S处的拟合系数；IGWR(S)为截距。

图4 环境因子空间分布

在ArcGIS中建立GWR模型，对其局部回归残差进行OK插值，得到εGWRK(S)，二者相加即为GWRK模型最终预测结果。

RK和GWRK模型残差空间相关性通过半变异函数模型进行分析，采用最优半变异函数模型的插值参数进行OK插值。

2.2.4 模型精度评价 模型验证基于10折交叉验证法，即每次随机选取1 179块样地的90%作为模型训练集，10%作为测试集。采用模型决定系数(R2)、均方根误差(root mean square error，RMSE)、绝对平均误差(mean absolute error，MAE)和赤池信息准则(Akaike information criterion，AIC)进行评价和选择，计算公式如下：

(6)

(7)

(8)

AIC=-2lgL+2p。

(9)

2.2.5 统计分析 因子筛选、MLR和RF模型建立、RK和GWRK预测模型残差半变异函数分析以及模型精度指标计算均基于R3.5.1 进行，地统计学的RK和GWRK预测模型实现及SI空间预测主要基于ArcGIS 10.2进行，SI与主要环境因子的相关和偏相关分析主要在R3.5.1软件中完成。

3 结果与分析

3.1 模型精度评价

由MLR、RF、RK和GWRK模型预测值与实测值残差平均分布(表3)可知，整体上，4种模型残差均为正态分布(P>0.05)； MLR和RF模型的空间自相关性明显(P<0.05)，RK和GWRK模型残差空间分布随机(P>0.05)。由表4可知，MLR、RF、RK和GWRK模型预测精度交叉验证结果相差较大： MLR模型R2最低，RMSE、MAE和AIC最大，模型精度最低； GWRK模型精度最高； RF模型精度一般； RK与GWRK模型精度相近。相比MLR模型，RF、RK和GWRK模型的平均预测精度均有较大提升，R2分别提高0.45、1.00和1.07倍，RMSE、MAE和AIC依次呈下降趋势； 地统计学模型显著优于RF模型，基于空间局部回归的GWRK模型对SI空间变异的预测能力优于基于空间全局回归的RK模型(R2提高0.4%，RMSE、MAE和AIC分别降低0.5%、1.0%和1.6%)。分别选取10折交叉验证中精度最高的MLR、RF、RK和GWRK模型预测研究区SI空间分布，虽然4种模型预测SI分布趋势大致相同，但MLR与RF模型预测结果“锯齿化”明显，而RK和GWRK模型预测结果高低值边缘曲线更加细腻，GWRK模型不同SI的小区域明显增多(图5)。综上所述，选择GWRK作为研究区SI的最终预测模型。

表3 模型残差描述性统计和正态性检验

3.2 SI与主要环境因子的空间分布格局

由表5可知，Moran’I的分布范围为0.28～0.84，表明模型系数具有显著空间自相关性(Z>2.58，P<0.01)，环境因子对SI的影响随地理位置改变而变化。由GWRK模型预测的SI空间分布(图5D)可知，立地指数整体上西北高、东南低，与DEM和BIO14的空间分布格局一致(图4A、B)，表明DEM和BIO14共同主导SI的东西差异，但DEM空间局部“聚类”较BIO14更显著(Moran’I: 0.65 > 0.38)，DEM对SI东西变异的主导性更强； TP和TN的空间分布格局较复杂，除西部边缘及中南部较高外，均呈南低北高的空间分布格局(图4C、D)，二者共同主导SI的南北差异。

表4 预测模型交叉验证统计量

将SI自动划分为Ⅰ～Ⅴ级，分别代表高(14.66～20.43 m)、较高(12.53～14.66 m)、中(10.39～12.53 m)、较低(8.40～10.39 m)和低(< 8.40 m)立地生产力(图5D)。第Ⅰ等级占研究区总面积的10.7%，主要分布于西北部坝上地区； 第Ⅱ等级占研究区总面积的15.3%，绝大部分分布于西部，北部分布适中，极少数(约占第Ⅱ等级的1.1%)分布于东部； 第Ⅲ等级占研究区总面积的25.7%，主要分布于第Ⅰ、Ⅱ等级和Ⅳ、Ⅴ等级边缘区域及东部； 第Ⅳ、Ⅴ等级占研究区总面积比例较大(48.4%)，集中于东部和南部的坝下地区。研究区华北落叶松人工林立地生产力区域两极分化较为明显，西北坝上和东南坝下地区SI差距较大，为7.99～15.43 m。

图5 不同预测模型SI空间分布

表5 GWR模型描述性统计

3.3 SI与主要环境因子的关系

研究区华北落叶松人工林SI预测值与主要环境因子的相关分析(表6)表明，在大的空间尺度上，SI与DEM、BIO14、TN和TP呈显著正相关关系(P< 0.01)； 偏相关分析表明，控制其他环境因子不变，各立地生产力等级内SI均对BIO14的变化响应较弱； 第Ⅰ、Ⅴ等级SI主要受限于TN，第Ⅱ、Ⅲ等级SI与DEM的变化均呈显著正相关关系(P< 0.01)，而第Ⅱ、Ⅲ等级SI分别随TP、TN增加而增大(P< 0.01)。 结合图5、6可知，华北落叶松人工林喜海拔较高(约> 1 450 m)、土壤TP含量(0.07～0.08 g·kg-1)以及TN含量(0.2～0.3 g·kg-1)适中的区域，土壤TN、TP含量过高/低均会限制华北落叶松人工林生长。

表6 各立地生产力等级SI与环境因子的偏相关关系①

图6 SI与其敏感环境因子的关系

4 讨论

本研究探讨研究区SI随环境因子的变化规律，进一步验证了重要环境因子作为辅助变量预测立地生产力的可行性(Bravo-Oviedoetal.， 2011； Jiangetal.， 2015)； 但不同立地生产力因研究区域以及林分类型等差异，筛选出的重要环境因子不尽相同(Huangetal.， 2017； 雷相东等， 2018)。综合考虑SI样本数量和环境因子数量，本研究最终选取与立地指数相关性高(P< 0.01)、彼此之间相关性低(VIF < 4)的海拔、最干月降水、土壤全氮和全磷作为影响研究区立地指数的重要环境变量，最干月降水在数值上存在较大变异(CV=43.75%)，与立地指数呈显著正相关关系(P< 0.01)，在空间上也具有较强空间相关性(Moran’I=0.54，P< 0.01)。为进一步探讨其对立地指数的影响，本研究保留该解释变量，剔除多重共线性问题严重(VIF > 15)的生长季温度和降水等气候因子。

立地生产力模型能够应用于森林经营实践至少需满足2个条件： 解释立地指数变异至少50%； 辅助变量易于获取(Blythetal.， 1981)。本研究除MLR模型外，其他3种预测模型均满足该条件，能够解释立地指数变异的54.4%～78.0%，其中地统计学模型(RK和GWRK)的解释能力更强(表4)。RK和GWRK模型预测残差均采用属于线性克里格的普通克里格进行插值，其主要差异是趋势项建模方法不同，RK建立的是环境因子与立地指数的全局回归模型(MLR)，而GWRK建立的是局部回归的地理加权回归模型(GWR)。由于MLR模型仅为整个研究区数据的平均化呈现，丢弃了大量关系空间变异和模型性能的潜在有价值信息(Fotheringhametal., 1998)，而GWR模型基于连续空间权重函数估测除观测样地外其他位置参数，克服了全局回归模型中不同空间位置环境因子权重相同的问题(Georganosetal.， 2017)，因此GWRK模型精度更高(表4)； 但也有研究得出相反结论(Jinetal.， 2016)，可能是受数据完整性和空间尺度大小等因素的影响。本研究为保证所有模型输入变量一致，只选取4个主要环境因子进行建模，虽然RF模型适用于数据维度较高的情况，且不需要预先选择解释变量(Farrellyetal.， 2011)，但当解释变量为与立地指数显著相关(P< 0.01)的15个环境因子时，其预测精度(R2=0.671, RMSE=1.817, MAE=1.346, AIC=281.824)低于地统计学模型，并不是最优模型。相比RF和地统计学模型，MLR模型只能解释预测变量的部分方差变异，与以往研究结论(Zhangetal.， 2017)一致。因此本研究认为，当预测变量具有显著空间变异规律或预测变量与解释变量之间并不是简单的全局线性关系时，应避免使用MLR模型。

立地生产力整体分布趋势与海拔和最干月降水分布规律一致，均由西北向东南递减(图5D)。研究区地处华北针叶林带，森林覆盖面积大，而东西海拔变化再分配影响林分生长的光照和地表土壤水热条件，进而对土壤养分及酶活性、土壤微生物、凋落物质量和根系产生影响(Fiereretal.， 2011； Wangetal.， 2009)，导致华北落叶松人工林生长差异。受地形影响，研究区最干季(一般为冬季)均温-17.88～-7.07 ℃，西北低、东南高。温带落叶乔木光合作用最低温度为-3～-1 ℃，而冬季木本植物的越冬器官-25 ℃仍未停止呼吸(潘瑞炽， 2008)，因此，在最干月华北落叶松呼吸作用会消耗较多的储存养分以维持基本生命活动，且在低温胁迫条件下可能会造成细胞缺水(Bray， 1997)。水是植物细胞中最重要的分子，降水会影响整个土壤-植物-大气连续体(soil-plant-atmosphere continuum，SPAC)的物质循环(Pengetal.， 2007； Philip， 1996)，而最干月降水量大的地区能够增强华北落叶松人工林对低温干旱条件的耐受能力(Bray， 1997)，因此立地指数高，但主要影响宏观区域尺度(表6)。土壤氮磷是植物生长发育的两大限制因子和必需的矿质元素(卢同平等， 2017)，土壤全氮和全磷含量增加在一定程度上可促进植物吸收氮磷营养元素，而植物氮磷含量反映植物养分的地理分布格局和生物地球化学循环的区域特征(樊江文等， 2014)，是对区域立地因子、生物群落以及人类活动等生物和非生物因子长期适应的结果(卢同平等， 2017)，如西北边缘地区的华北落叶松人工林生长受海拔和土壤全磷含量的共同限制(图5、表5)。温带华北落叶松人工林土壤全氮、全磷含量过高或过低均会导致生态系统中氮磷的不平衡输入，可能诱发植物从氮限制向磷限制转变(Dengetal.， 2016)，不利于华北落叶松生长发育，本研究进一步验证了该结果(表6、图6)。总体来看，华北落叶松人工林立地生产力等级高的地区海拔约高于1 400 m(图6)，且在土壤全氮和全磷含量适中的地区竞争优势更为突出。在森林经营实践中，管理者应充分认识到立地指数是林木生长特性与多个环境因子综合作用的结果(Littkeetal.， 2016)，从而采取有效的经营措施提高立地质量，促进树木生长。

5 结论

空间明晰的森林立地生产力图的绘制以及立地生产力与重要环境因子之间关系的确定一直是森林经营实践中的热点。本研究建立的立地指数与主要环境因子之间的地理加权回归克里格模型，可用于绘制立地指数空间分布图，探讨研究区华北落叶松人工林立地生产力对主要环境因子的响应机制。在景观尺度上，华北落叶松人工林立地指数对主要环境因子响应较强，海拔、最干月降水较高以及土壤氮磷含量适中的研究区西北部，是华北落叶松人工林适生区，而海拔、最干月降水较低以及土壤氮磷含量不平衡的东、南部，华北落叶松人工林的适应性较差。因此，建议在东、南部地区的森林经营实践中，可通过控制林分密度(约2 170株·hm-2)(任丽娜等， 2012)、营造混交林来缓解气候变化和海拔对树木生长的限制，或向土壤中添加适量的氮磷肥以平衡土壤氮磷比，从而促进华北落叶松人工林生长以及生产力提高。