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

采样点位于青海省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。

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

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

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

式中，$ {y}_{ijk} $是第$ i $个第一水平中的第$ j $个第二水平内的第$ k $次观察值，$ m $、$ {m}_{i} $分别是第一水平、第二水平的分组数量，$ {n}_{ij} $是第$ i $个第一水平中的第$ j $个第二水平内的观测次数，$ f $是含有参数向量$ {\varphi }_{ijk} $和协变量向量$ {v}_{ijk} $的非线性函数，$ {A}_{ijk} $是设计矩阵，$ \beta $是（p × 1）维固定效应向量，$ {B}_{i,jk}{b}_{i} $、$ {B}_{ijk}{b}_{ij} $分别是第一水平、第二水平的随机效应设计矩阵，$ {b}_{i} $、$ {b}_{ij} $分别是第一水平、第二水平的随机参数向量，$ {D}_{1} $、$ {D}_{2} $分别是第一水平、第二水平的随机参数方差-协方差矩阵，$ {b}_{i} $和$ {b}_{ij} $不相关，$ {\epsilon }_{ijk} $是服从正态分布的误差项，$ {\sigma }^{2} $是方差，$ {R}_{ij} $是第$ i $个第一水平中的第$ j $个第二水平内的方差-协方差矩阵。

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

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

以林木生长到胸高位置的年龄为起始年限（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  不同起源与坡位青海云杉单木胸径生长曲线

Figure 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个相交点。

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

图  2  不同起源单木胸径生长拟合曲线对比

Figure 2.  The fitting curves for Individual-tree DBH growth models at different origins

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

图  3  不同坡位青海云杉天然林单木胸径生长拟合曲线对比

Figure 3.  The fitting curves for Individual-tree DBH growth model for Picea crassifolia natural forests at different slope positions

图  4  不同坡位青海云杉人工林单木胸径生长拟合曲线对比

Figure 4.  The fitting curves for Individual-tree DBH growth model for Picea crassifolia plantations at different slope positions

根据生长模型拟合结果，本研究分别以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。

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

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

图  5  青海云杉单木胸径生长拟合曲线

Figure 5.  Individual-tree DBH growth curve for Picea crassifolia

总体来看，青海云杉生长到胸高位置后（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模型）。各最优胸径生长模型的R²均在0.913以上， TRE和MSE均在 ± 2%以内，MPE大多在5%以内，MPSE在35%以内（表2～3）。这与同一地区生长模型研究得到的结果接近。马克西等^[34]构建了青海省4个主要树种组的胸径生长率模型，其中模型的R²均在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]。因此，未来气候变化下的森林管理与保护应充分考虑更多的地形、气候差异对林木生长的影响。

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