青海省不同生境下青海云杉胸径生长模型研究

马 浩，陈科屹，徐干君，党禹杰，何友均，王建军＊

(1. 国家林业和草原局西北调查规划院，陕西西安 710048；2. 中国林业科学研究院林业科技信息研究所 北京 100091)

单木生长模型以单木为基本单位，通过数学模型来模拟林木的生长过程，揭示林木生长规律，并且能够提供较详细的林分结构特征及动态变化信息，对于森林经营管理和保护等方面有十分重要的实践价值[1-3]。例如，李慧婷[2]采用5 种理论生长方程建立了青海海东地区华北落叶松的树干材积生长模型。与传统的非线性回归方法相比，混合效应模型能够有效提高模型的估计精度[3-8]。李春明和唐守正等[4]建立了考虑样地的随机效应、观测数据的时间序列相关性及间伐强度的落叶松云冷杉林断面积生长混合模型，结果表明混合效应模型显著地提高了林分断面积的预估精度。

单木生长模型在区域水平上的应用具有一定的局限性，其影响因素主要有树木立地条件、人为经营措施、空间结构多样性等[9-11]。张雷等[12]研究发现，随海拔的升高，祁连山北坡青海云杉林的结构出现明显变化，平均胸径增加，平均树高呈“单峰”变化。张中惠等[13]探究宁夏六盘山华北落叶松单木树高生长对立地因子的响应，发现影响最大的立地因子是海拔，其次是坡向、坡度。海拔、坡向与坡位在地形上改变了温度、湿度、光照等资源的分配，进而影响植被分布与林木生长发育[14]。因此，研究林木在不同环境下的生长特点，不仅对科学描述林木生长动态十分重要，还能为森林多种功能与效益的发挥提供基础参考。

青海云杉（Picea crassifoliaKom.）是青海地区主要的森林更新树种和造林树种，在涵养水源、水土保持、改善环境等方面的作用显著。国内外学者对该树种的生理生态特征、树轮生长特性等方面进行了广泛研究[15-17]。但是，关于青海云杉单木生长模型的相关研究较少[2]。因此，本研究以青海云杉为研究对象，分析不同起源、不同坡位下单木胸径生长的规律，构建单木胸径生长模型和考虑起源和坡位影响的混合效应模型，旨在为今后有效保护、合理经营青海云杉林提供经验模型和有益参考。

1 研究区概况

青海位于我国西北部，地理坐标为31°36'～39°19' N，89°35'～103°04' E，属于高原大陆性气候，冬季寒冷干燥，夏季短暂且凉爽湿润，日照时间长，太阳辐射强。年平均气温-5.1～9.0℃，年均降水量15～750 mm，降水分布由东南向西北逐渐减少[18]。该地区的地形、地貌复杂多样，地势总体呈西北高东南低，东部多山，西部为高原和盆地，河流湖泊分布多，平均海拔3 000 m 以上[18-19]。主要的植被类型有常绿针叶林、落叶阔叶林、高寒灌丛、高寒草甸等，主要建群种有青海云杉、祁连圆柏（Juniperus przewalskiiKom.）、青杨（Populus cathayanaRehd.）、油松（Pinus tabuliformisCarr.） 等， 主要灌木有沙棘（Hippophae rhamnoidesLinn.）、柠条（Caragana korshinskiiKom.）、锦鸡儿属（CaraganaFabr.）等。

2 材料与方法

2.1 数据收集与处理

采样点位于青海省3 个县：尖扎县、祁连县、大通回族土族自治县（简称“大通县”）。采样点选取受人为干扰少、远离道路及防火道、环境较为稳定的地带。在不同起源、不同立地条件下的典型林分中选取标准木作为样木。其中，以坡位划分，在不同海拔、坡向、坡度内进行采样，保证数据的区域代表性和均匀性。于2022 年7 月采集青海云杉样木共172 株。其中，天然起源的青海云杉74 株（上坡19 株，中坡9 株，下坡22 株，平坡24株），人工起源的青海云杉88 株（上坡29 株，中坡15 株，下坡28 株，平坡16 株）。在每1 株样木的胸径高（1.3 m）的位置，用生长锥分别在南、北2 个方向上钻取树芯，将树芯样本保存并分类编码，记录样木胸径（DBH）、树种、起源、坡位等信息。

将样本带回实验室后，按树轮年代法的标准方法进行风干、固定、切割/打磨等预处理；使用LinTab 年轮分析仪对树轮宽度测定，测定精度为0.01 mm。年轮条测量获得的测量值为某一方向上的半径生长量，对测量值的平均值乘以2，得到该年龄下的胸径生长量。计算胸径平均生长量（MAI）和胸径连年生长量（CAI），龄阶间距为5 a。样本信息见表1。

表1 样本信息 Table 1 Basic information of samples

此外，在典型天然林中选取4 株标准木进行树干解析，记录树干长、胸径以及圆盘的直径与年龄，以确定林木地径生长到胸径位置时的年限。解析木的胸径范围为25.9～34.8 cm，树龄范围为69～135 a，树高范围为16.1～26.5 m，到达胸径位置的年龄范围为21～26 a。根据解析木数据分析，青海云杉地径生长到胸径位置时的年限为24 a。树龄为树芯的年轮数加上林木生长到胸径位置时的年限。

2.2 研究方法

2.2.1 生长模型建立 结合树木的生长规律及前人的研究成果，本研究选取理查德（Richards）方程、逻辑斯蒂（Logistic）方程、考尔夫（Korf）方程、坎派兹（Gompert）方程和单分子式（Mitscherlich）方程5 种树木生长理论方程对胸径生长量与年龄的关系曲线进行拟合[2,20]。利用ForStat 2.2 统计软件拟合生长模型。

2.2.2 非线性混合效应模型建立 非线性混合效应模型的建立考虑了回归函数依赖于固定和随机效应的非线性关系[21]。多水平非线性混合效应模型的形式（以两水平模型为例）如下：

式中，yijk是第i个第一水平中的第j个第二水平内的第k次观察值，m、mi分别是第一水平、第二水平的分组数量，nij是第i个第一水平中的第j个第二水平内的观测次数，f是含有参数向量φijk和协变量向量vijk的非线性函数，Aijk是设计矩阵，β是（p×1）维固定效应向量，Bi,jkbi、Bijkbij分别是第一水平、第二水平的随机效应设计矩阵，bi、bij分别是第一水平、第二水平的随机参数向量，D1、D2分别是第一水平、第二水平的随机参数方差-协方差矩阵，bi和bij不相关，εijk是服从正态分布的误差项，σ2是方差，Rij是第i个第一水平中的第j个第二水平内的方差-协方差矩阵。

对于参数效应的确定，本研究将所有不同随机效应参数组合的模型都进行拟合。模型的拟合优度指标为赤池信息准则（AIC）、贝叶斯信息准则（BIC）和对数似然值（LogLik），一般采用最小AIC和BIC，以及最大Loglik值的标准确定具有最优参数组合的模型。随机效应参数的方差-协方差矩阵D反映了随机效应在个体之间的差异性。采用常见的广义正定矩阵作为随机效应参数的方差-协方差矩阵。方差-协方差矩阵Rij主要用于解决数据中存在的自相关和异方差问题[22]。由于胸径生长量数据不是重复观测数据，因此不考虑数据的自相关问题。通过加权回归方法来消除异方差问题，常用的异方差函数有常数加幂函数、幂函数和指数函数。利用R 语言的nlme 包进行混合效应模型的参数估计。

2.2.3 模型评价 为了对模型的拟合结果进行评价和比较，采用的指标有决定系数（R2）和均方根误差（RMSE）。采用全部数据计算模型估计值的精度指标，选择总体相对误差（TRE）、平均系统误差（MSE）、平均预估误差（MPE）和平均百分标准误差（MPSE）4 项误差指标对模型进行综合评价。根据R2较大，RMSE、TRE、MSE、MPE和MPSE较小原则选取最优模型。

3 结果与分析

3.1 胸径生长过程

以林木生长到胸高位置的年龄为起始年限（24 a），5 a 为龄阶，分析不同起源、坡位生长下青海云杉的胸径生长过程（图1）。CAI与MAI均呈现“升—降—平缓变化”的趋势，MAI大于CAI。其中，天然林内青海云杉的快速生长期为29—44 a，CAI与MAI均在0.40 cm 以上；29 a时CAI与MAI相交且有最大值，分别为0.509 cm、0.515 cm。人工林内青海云杉的速生期在29—39 a，CAI与MAI均大于0.40 cm，34 a 时CAI与MAI有最大值，分别为0.473 cm、0.434 cm；CAI与MAI出现多个相交点，MAI在39 a 快速下降后变化不大，CAI在树龄164—259 a 波动较大，在244 a 出现第2 个峰值0.302 cm。

图1 不同起源与坡位青海云杉单木胸径生长曲线Fig. 1 Individual-tree DBH growth for Picea crassifolia at different origins and slope positions

天然林不同坡位中，青海云杉在上坡的CAI与MAI在29—84 a 均保持0.35 cm 以上的较高值，为生长速生期，并且出现2 个相交点；82 a之后CAI先降后升。在中坡的青海云杉速生期为29—49 a，其中CAI与MAI在29—39 a 均大于0.50 cm，29—94 a 之间波动幅度较大，出现3 处相交点。在下坡的青海云杉的CAI与MAI变化平缓，速生期为29—64 a，24—104 a 之间两者出现3 个相交点。在平坡的青海云杉CAI与MAI逐年下降，速生期为29—54 a，其中CAI与MAI在29—44 a 均在0.50 cm 以上。

人工林不同坡位中，青海云杉MAI与CAI多次相交。在上坡的青海云杉CAI与MAI先降再平缓变化，生长速生期为29—39 a，CAI与MAI均在0.35 cm 以上；229—289 a 之间CAI的波动幅度较大。在中坡的青海云杉CAI与MAI变化平缓，但两者均小于0.30 cm。在下坡的青海云杉CAI与MAI逐年下降后平缓变化；速生期为29—39 a；CAI波动幅度较大，159—204 a 时期CAI先升后降。在平坡的青海云杉呈现先升后降的变化趋势，速生期在29—54 a，出现1 个相交点。

3.2 胸径生长模型的拟合与评价

3.2.1 胸径生长模型的拟合与评价 采用5 种生长模型拟合不同起源青海云杉单木胸径生长量与年龄关系曲线，模型拟合效果良好（表2，图2）。天然林单木生长模型中，Gompertz 模型的拟合结果表现最好，R2为0.915；对于人工林，Korf 模型的拟合效果最优，R2为0.946。2 个模型的TRE在 ± 1%以内，MSE在 ± 2%以内，MPE均在3%以内，MPSE在35%以内。

图2 不同起源单木胸径生长拟合曲线对比Fig. 2 The fitting curves for Individual-tree DBH growth models at different origins

表2 不同起源最优单木胸径生长模型参数及评价指标Table 2 The parameter estimates and evaluation indices of optimal individual-tree DBH growth models at different origins

不同坡位青海云杉单木胸径生长模型的拟合结果如表3、图3～4 所示。上坡青海云杉单木胸径生长拟合效果最优的生长模型是Richards 模型，中坡的最优生长模型是Logistic 模型，下坡的最优生长模型是Gompertz 模型，平坡的最优生长模型是Korf 模型。其R2分别为0.982、0.915、0.913、0.996，TRE和MSE在 ± 2%以内，MPE在6%范围内，MPSE大部分在35%范围内。人工林不同坡位下，上坡、中坡、下坡和平坡的最优模型分别是Korf 模型、Richards 模型、Korf 模型、Logistic模型。其R2分别为0.940、0.973、0.980、0.955，TRE和MSE在 ± 2%以内，MPE在5%范围内，MPSE在35%范围内。

图3 不同坡位青海云杉天然林单木胸径生长拟合曲线对比Fig. 3 The fitting curves for Individual-tree DBH growth model for Picea crassifolia natural forests at different slope positions

表3 不同坡位青海云杉最优单木胸径生长模型参数及评价指标 Table 3 The parameter estimates and evaluation indices of optimal individual-tree DBH growth models at different slope positions

3.2.2 混合效应模型的拟合与评价 根据生长模型拟合结果，本研究分别以Gompertz 模型、Korf 模型和Logistic 模型作为基础模型，采用全部数据，在基础模型中考虑起源、坡位的混合效应，对所有不同随机参数组合进行拟合和检验，各基础模型的最优模型结果见表4。由表可知，混合效应模型的检验结果均优于基础模型，考虑混合效应显著提高了模型的拟合效果。Gompertz 模型、Korf模型、Logistic 模型的最优随机参数组合分别为(a、c)、(b)、(b、c)。以Korf 模型为基础模型的混合效应模型具有最小的AIC、BIC和最大的LogLik，但是Korf 模型在考虑异方差函数过程中不收敛。因此，选择拟合结果其次的Gompertz 模型（AIC=12 700.520，BIC= 12 757.610，LogLik= -6 340.259），并进行异方差校正。由表4 可知，3 种异方差函数均能改善模型的拟合效果。其中，添加幂函数的混合效应模型的拟合结果最优，AIC为12 394.490，BIC为12 457.300，LogLik为-6 186.247。

表4 混合效应模型的模拟结果 Table 4 Fitting results of mixed-effect model

青海云杉胸径生长基础模型和混合效应模型的参数估计以及拟合结果如表5 所示。从表中可以看出，最优混合效应模型具有更好的拟合精度，其R2为0.702，较基础模型提高了32.7%；RMSE为4.14，较基础模型下降了20.8%。十折交叉验证法检验结果显示，混合效应模型的检验指标均优于基础模型。混合效应模型TRE、MSE、MPE、MPSE的降幅分别为89.3%、83.5%、20.6%、15.1%。

3.2.3 模型的预测 利用最优模型对青海云杉胸径的生长拟合见图5。由图可以看出，青海云杉胸径总生长量随年龄逐渐增大，其中天然林上坡和人工林下坡的单木胸径生长趋势较好。天然林整体、天然林中坡和天然林下坡的单木连年生长量CAI先增大后减小，平均生长量MAI逐年下降；天然林上坡与平坡的单木CAI和MAI均呈现逐年下降趋势。除了人工林平坡，人工林的单木CAI和MAI均逐年下降。

图5 青海云杉单木胸径生长拟合曲线Fig. 5 Individual-tree DBH growth curve for Picea crassifolia

4 讨论

总体来看，青海云杉生长到胸高位置后（24 a），单木胸径生长量随着年龄的增加而逐渐减少（图1）。青海云杉生长缓慢，树龄20 a 左右的林木树高约1.4 m[23-24]。温度、降水、光、空间等是该地区林木生长的限制因素[25-27]。随着胸径生长量的增大，林木个体在林分中能够获取的可利用资源变少，养分吸收的有效性降低，使得林木在后期生长缓慢。不同起源对比，天然起源的青海云杉在树龄29—44 a 保持相对较高的单木胸径生长量，CAI和MAI均在0.40 cm 以上，比人工林的持续时间长（29—39 a）；随后CAI和MAI的变化趋势平缓（图1）。在西北地区已有类似的研究结果[28-30]。王学福[29]研究发现祁连山青海云杉胸径连年生长量高峰期在40—50 a，平均生长量高峰期在90—110 a。杨文娟[30]指出青海云杉单木胸径的最大连年生长量出现在40 a 时，最大平均生长量出现在70 a 时，分别为0.30 cm、0.22 cm。不同坡位条件下，胸径的生长规律具有差异（图1）。坡位上的差异主要表现在土壤肥力和土壤水分方面，其中低坡位能够最有效地提供林木生长所需的养分，进而促进林木的径向生长[31-33]。苏妮尔等[32]指出下坡位的土壤养分含量最高，是大径材红皮云杉培育的最优坡位。

单木胸径生长模型的拟合结果显示，不同生境压力下林木生长过程具有不同的理论生长方程，但是大部分都符合传统的“S”形生长曲线（除了Mitscherlich 模型）。各最优胸径生长模型的R2均在0.913 以上，TRE和MSE均在 ± 2%以内，MPE大多在5%以内，MPSE在35%以内（表2～3）。这与同一地区生长模型研究得到的结果接近。马克西等[34]构建了青海省4 个主要树种组的胸径生长率模型，其中模型的R2均在0.93 以上，MPE均在0.3%之内；胸径生长量预估的MPE均在2%之内，MPSE均在60%左右。

在基础模型上考虑起源和坡位对胸径生长量的影响，建立了青海云杉单木胸径生长量混合效应模型，结果显示混合模型的拟合效果均优于基础模型（表4～5）。在模型检验中，混合模型仍然呈现较高的拟合精度，其TRE、MSE、MPE和MPSE较基础模型分别下降89.3%、83.5%、20.6%、15.1%。这与大多数学者的研究结论一致[4,6,8,21,26]。此外，模型拟合的精度还可能受到海拔、温度、降水、干旱等因素的影响[12,17,35]。因此，未来气候变化下的森林管理与保护应充分考虑更多的地形、气候差异对林木生长的影响。

5 结论

不同生境压力下的青海云杉胸径生长有不同的变化规律。单木胸径生长量随着年龄的增加而下降，其中天然林的速生期较人工林的持续时间长，不同坡位条件下胸径生长规律具有差异。青海云杉单木胸径混合效应模型的拟合效果优于基础模型。因此，在未来的森林经营及保护中，这些模型适用于青海省不同生境压力下青海云杉的林木胸径生长量估算。