六盘山华北落叶松单木树高对立地因子和林分特征的响应

2022-03-26张中惠王彦辉郭建斌

林业科学研究 2022年1期

张中惠，王彦辉，郭建斌＊，王晓

(1.北京林业大学水土保持学院，北京 100083；2.中国林业科学研究院森林生态环境与保护研究所，国家林业和草原局森林生态环境重点实验室，北京 100091)

开展森林精准经营[1]，是未来林业的发展方向，这需区分考虑林分内每株树的生长特征及受不同因素的影响，其中单木树高是最重要的林木特征之一。华北落叶松是我国北方暖温带湿润半湿润气候区山地的主要造林树种，具有速生及适应性强等优点[2]。华北落叶松单木树高生长过程同时受林分结构因子（如林龄、密度、郁闭度等）[3]及立地因子（气候、地形、土壤等）[4]的影响。对华北落叶松林的前期研究表明，林龄直接影响单木树高生长和其他林分结构特征变化[5]；林分密度和林分郁闭度等林分结构可影响树木对环境资源（养分、水分、光照等）的占有和利用，从而影响单木树高的生长、竞争、死亡和产量[6-7]。影响林木树高生长的主要立地因子是海拔、坡向、坡度和土层厚度等[8]，它们主要通过影响光照、温度、养分和水分而影响树高生长过程[9]。单木模型以树木个体自身生长特征为基础，从林木生长的竞争机制出发，模拟单木生长过程[10]。在以往研究中，采用解析木数据建立了兴安落叶松（Larix gmelinii(Rupr.)Kuzen）的单木树高模型，但并未考虑林分特征因子以及立地因子对单木树高生长的影响[11]；为反映多个因素对单木树高生长的综合影响，可用外包线法首先确定生长指标对单因子变化的响应函数，然后建立多因子耦合模型，如杨文娟[12]在祁连山研究青海云杉（Picea crassifoliaKom.）林分生长响应林龄、密度、海拔和坡向的模型，田奥建立了六盘山华北落叶松林分树高与林龄、密度和海拔的耦合模型[3]。在六盘山区华北落叶松林以往的研究中都只建立林龄-树高模型[13]，但还没进行过单木树高生长的多因素影响及耦合模型研究。因此，本研究在六盘山区选择华北落叶松人工林典型样地，开展每木调查及解析木分析，定量研究单木树高受立地条件和林龄及其他林分结构的影响，并建立能反映多因子影响的树高生长耦合模型，以便为选择适宜造林立地、预测单木树高生长和林分结构变化、开展森林精准经营等提供理论和技术支撑。

1 材料与方法

1.1 研究区概况

本研究在宁夏六盘山南段东坡的香水河小流域（106°12′10.6″～106°16′30.5″ E，35°27′22.5″～35°33′29.7″ N）进行，海拔变化在2 070～2 931 m，属温带半湿润气候，年均气温3.7 ℃，年均降水量671 mm[14]。流域内土壤以灰褐土为主。小流域森林覆盖率高达82.91%；以华山松（Pinus armandiiFranch.）、白桦（Betula platyphyllaSuk.）等天然次生林为主，占小流域面积的58.51%；人工林以华北落叶松纯林为主，占小流域面积的24%[15]；灌丛面积占12.01%，主要有西北栒子（Cotoneaster zabeliiSchneid.）、沙棘（Hippophae rhamnoidesLinn.）等；草地面积占16.25%，主要有狼针茅（Stipa baicalensisRoshev.）、早熟禾（Poa annuaL.）等；草甸面积占4.64%，主要有苔草（Carex tristachyaSpp.）、蕨（Pteridium aguilinumL.）等[15]。

1.2 研究方法

1.2.1 样地布设经全面踏査后，沿海拔梯度（2 000～2 200、2 200～2 400、2 400～2 600、2 600～2 800 和2 800～3 000 m）选择了23 块面积20 m ×20 m 的华北落叶松纯林样地（表1），记录样地经纬度及海拔、坡度、坡向等立地条件。调查中，将正北方向记为0 度，顺时针偏离正北180 度以内的坡向为正，逆时针偏离正北180 度以内的坡向为负。

表1 华北落叶松人工纯林样地基本情况Table 1 Basic information on sample plots of pure plantation of Larix principis-rupprechtii

1.2.2 树木调查调查林分的郁闭度、密度，然后对各样地实测全部林木的树高、胸径、冠幅及枝下高等。用相对树高表征树木优势度[16]。优势度计算公式为：

然后按Kraft 树冠优势度分级标准[17]，将各样地内的林木个体优势度分为3 级：I 代表优势木、Ⅱ代表平均木、Ⅲ代表被压木，每个样地至少选择优势木和平均木各2 株及被压木1 株，进行解析木调查（共116 株）。

1.2.3 多因素影响的树高生长耦合模型建立上外包线法可以剥离出各单一因子影响并分析各单一因子影响以确定响应函数类型[18]。建立模型过程：（1）将林龄对应的各树高值（H）除以理论最大树高Hmax（33 m），得到相对树高[3]；把相对树高变化范围分成若干区段，在每段中选取大于本段数据平均值加一倍标准差的数据点[19]，或在一些数据偏少的区段中直接选择最大数据点，拟合外包线f(x1)；（2）同理，将步骤1 中消除了林龄因子影响的数据与下一个影响因子做散点图，并得到响应函数f(x2)；（3）用同样方法逐个消除其他因子的影响，直到确定相对树高对最后一个因子的响应函数f(xn)；（4）连乘响应函数，得到受多因素影响的树高生长耦合模型：

式中，H为单木树高，f(x1)、f(x2)、...等分别表示单木树高生长对各单一因素x1、x2、...等的响应函数，n为考虑的影响因子数。

把116 株解析木数据分成：92 株解析木实测数据拟合模型参数，剩余24 株实测解析木数据验证模型拟合优度，选择验证解析木数据时注意了保持解析木在林龄、密度、郁闭度、海拔、坡度、坡向、优势度范围内的均匀分布。评价模型拟合优度的指标包括决定系数R2、均方根误差RMSE、总相对误差TRE。R2越接近1，RMSE和TRE越小，说明模型拟合越好。计算公式为:

式中，Yi和为第i个单木树高的实测值和预测值，为所有实测树高的平均值，N为用于参数率定或检验的样本数。

1.2.4 数据处理用统计Excel 软件处理野外调查数据，用统计软件Origin 进行绘图及曲线拟合，并用SPSS25.0 软件进行典型相关分析等统计分析，用软件1st0pt 进行耦合模型的参数拟合。

2 结果与分析

2.1 单木树高对立地因子及林分特征的响应

2.1.1 单木树高与立地因子及林分特征的相关性单木树高与各影响因子的相关分析表明（表2），各因子的相关系数排序为：林龄＞海拔＞郁闭度＞优势度＞林分密度＞坡向＞坡度。

表2 单木树高与各因子的相关性Table 2 Correlation between tree height and influencing factors

2.1.2 单木树高对各单一因子的响应基于所有解析木数据分析得到的单木树高对林龄的响应见图1-A 中的上外包线，单木相对树高随林龄增加而逐渐升高，在21 a 以前增长较快，随后增速变缓，在45 a 时最高可达27.68 m。确定树高生长过程符合Richard 方程：f(age)=a×(1-exp(-b×age))∧c，R2=0.997 5。单木树高对海拔的响应见图1-B，树高在海拔2 000～2 200 m 随海拔升高逐渐变大，之后随海拔继续升高而逐渐下降。树高对海拔的响应呈三次多项式：f(ele)=d×ele3+e×ele2+f×ele+g，R2=0.856 7。单木树高对郁闭度的响应见图1-C 中的外包线，单木树高在郁闭度0.3～0.56 范围内随郁闭度增大而逐渐增大，在郁闭度＞ 0.56 后逐渐减小。树高对郁闭度的响应呈二次多项式：f(cd)=h×cd2+i×cd+j，R2=0.662 2。单木树高对优势度的响应特征见图1-D 中的外包线，可知单株树高随优势度增加先平缓增加，在优势度为-0.2 后增速加快，在优势度为0.08 后增速渐趋平缓。树高对优势度的响应符合S型曲线：f(dom)=k+(l-k)/[1+exp((dom-m)/n)]，R2=0.942 9。单木树高对林分密度的响应见图1-E 中的外包线，可知树高随密度增加在400～1 200 株·hm-2范围内逐渐增加，之后逐渐下降。单木树高对林分密度的响应呈二次多项式关系：f(den)=o×den2+p×den+q，R2=0.835 6。单木树高对坡向的响应见图1-F 中的外包线，可知最适坡向为阴坡和半阴坡，在坡向180°到0°的范围内，单木树高随坡向靠近正北方向逐渐增加。响应函数为二次多项式：f(asp)=r×asp2+s×asp+t，R2=0.522 4。单木树高对坡度的响应见图1-G 中外包线，单木树高在坡度达23.5°之前随坡度增大而增加，之后转而减小，最适坡度在20°～27°。响应函数为：f(slope)=u×slope2+v×slope+w，R2=0.649 7。

图1 单木树高对立地条件和林分特征的响应Fig.1 Response of single tree height to site conditions and stand structure characteristics

2.2 单木树高的多因子耦合预测模型建立与检验

按相关系数的大小排序，从单木树高响应林龄的模型开始，逐个增加耦合其他因子，得出模型（式6～式12）。为便于生产应用，本研究同时建立包括林龄、林分密度、海拔、坡向、坡度这些独立变量的模型（式13），然后在此基础上增加优势度建立了模型式（14）。从表3 可看出，从式（6）到式（10），伴随R2升高，RMSE和TRE降低；式（11）为加入坡向后R2降低，RMSE和TRE均上升；式（12）综合考虑了所有因子，其R2升高；式（13）只考虑了独立因子，所以R2下降；式（14）在式（13）的基础上加入了优势度，导致R2升高到和式（12）相近的水平。式（6）～（14）的TRE不超过 ± 0.3%，式（10）以及式（12）的TRE最小，这3 个评价指标均证明，式（12）是最优的多因子耦合模型，其次为式（14）。

表3 单木树高生长模型及其拟合优度Table 3 Forecasting and calculating models of tree height growth per plant

用24 株解析木的实测数据对模型（6～14）进行检验（如图2）：式（10）、式（12）、式（14）比其他模型拟合的要好，预测点最靠近并均匀分布在45°线两侧。

图2 单木树高生长模型（式6～14）的检验Fig.2 Test of tree heights of models of Eq.6-14

3 讨论

3.1 林分特征因子和立地因子对单木树高的影响

林龄是影响树高生长的最重要因子。华北落叶松树高随林龄增加的生长过程表现为“S”型曲线，可划分为生长初期、快速生长期、平缓生长期3 个阶段[20]。如图1-A 中下外包线表示的树高生长过程S 型曲线特征比上外包线更明显，本研究表明采用Richard 方程能很好表示树高生长过程[21]。林冠郁闭度直接或间接地影响林木生长及众多服务功能[22]，本研究中，华北落叶松单木树高在郁闭度为0.49～0.64 时最好，相比多功能森林经营要求的合理郁闭度范围0.6～0.8[23]有些偏低但也差异不大。优势度能直接反映林木对自然资源的利用能力[24]，单木树高随优势度增大过程呈“S”型曲线特征，在优势度超过0.08 后，树高增大速率趋于平缓，说明优势木和亚优势木的树高可能会接近或达到立地质量决定的潜在树高，是获得木材生产经济效益的主要贡献者。林分密度会影响单木的光照、水分、养分等资源占有量[6]，华北落叶松单木树高生长的最适密度平均为1 100～1 300 株·hm-2[25]。海拔虽不直接影响树木生长，但会通过影响降水、温度、蒸散等环境因子间接影响树木生长[26-27]。本研究表明六盘山区最适合华北落叶松生长的海拔范围为2 000～2 400 m，降水不足是六盘山区树木生长的主要限制因子，在海拔超过2 200 m 后，海拔升高导致的气温降低逐渐成为限制树木生长的主要因素，使得树高随海拔增加而逐渐降低[3]。坡向与坡度会影响地面得到的太阳辐射，阳坡接受的太阳照射比阴坡更多，坡度大时接受太阳垂直照射的程度越大[28]，从而影响潜在蒸散和水分条件[3]。本研究表明，单木树高在阴坡和半阴坡最大，这是因光照较弱、蒸散较小导致植物可用水分较多和受干旱胁迫较轻；单木树高随坡度增大表现为先增后减，最适坡度范围为20°～25°，类似于河北塞罕坝林场的华北落叶松生长研究结果[27]，因地势平缓时导致的林木相互遮光明显以及坡度过陡导致的土壤偏薄和水土流失均不利于林木生长。

3.2 多因素耦合树高模型的优点及其可能应用

在自然环境中，林木生长同时受多因子影响，因而时空差异很大，所以综合考虑多因子影响的模型预测精度一般较高。然而，由于缺少调查资料或追求简单易用等原因，以往的华北落叶松生长模型常仅考虑林龄影响[8]，忽略其他因子影响，降低了预测精度和限制了应用范围。相比之下，建立多因素耦合模型，便于准确预测在不同立地及林分结构条件下的生长变化。在以往的六盘山华北落叶松林分生长指标模型中考虑了林龄、密度和海拔这3 个主要因素[3]，在祁连山青海云杉林分生长指标模型中考虑了林龄、密度、海拔和坡向[10]。在本研究中，经过检验比较模型（6～14）得出，六盘山华北落叶松单木树高生长模型在综合考虑因子的影响时精度最佳，优点是能考虑所有立地条件和林分特征因子对树高生长的影响，会更准确地预测树高。模型在考虑林龄时能准确预测立地条件和林分结构平均情况下的树木生长；考虑海拔的影响可反映树高生长随海拔变化的空间差异；继续增加考虑了郁闭度和林分密度的影响，能在一定程度上反映林分结构对树木生长的影响，利于以林分为空间单元预测单木树高的时空变化；增加优势度的影响能精确到单株精准经营；继续增加考虑坡向和坡度，因为本研究区位于半湿润区，树高的坡向响应较弱，且本研究选取不同坡向的样地较少，无法准确体现坡向对树高生长的影响，在未来研究中应增加所缺坡向的样地，进一步研究不同坡向下的树高生长差异，使模型更精准。综上，可根据研究目的或生产需求，从这里建立的单木树高耦合模型（式6～14）中选择应用。