-
单木生长模型以单木为基本单位,通过数学模型来模拟林木的生长过程,揭示林木生长规律,并且能够提供较详细的林分结构特征及动态变化信息,对于森林经营管理和保护等方面有十分重要的实践价值[1-3]。例如,李慧婷[2]采用5种理论生长方程建立了青海海东地区华北落叶松的树干材积生长模型。与传统的非线性回归方法相比,混合效应模型能够有效提高模型的估计精度[3-8]。李春明和唐守正等[4]建立了考虑样地的随机效应、观测数据的时间序列相关性及间伐强度的落叶松云冷杉林断面积生长混合模型,结果表明混合效应模型显著地提高了林分断面积的预估精度。
单木生长模型在区域水平上的应用具有一定的局限性,其影响因素主要有树木立地条件、人为经营措施、空间结构多样性等[9-11]。张雷等[12]研究发现,随海拔的升高,祁连山北坡青海云杉林的结构出现明显变化,平均胸径增加,平均树高呈“单峰”变化。张中惠等[13]探究宁夏六盘山华北落叶松单木树高生长对立地因子的响应,发现影响最大的立地因子是海拔,其次是坡向、坡度。海拔、坡向与坡位在地形上改变了温度、湿度、光照等资源的分配,进而影响植被分布与林木生长发育[14]。因此,研究林木在不同环境下的生长特点,不仅对科学描述林木生长动态十分重要,还能为森林多种功能与效益的发挥提供基础参考。
青海云杉(Picea crassifolia Kom.)是青海地区主要的森林更新树种和造林树种,在涵养水源、水土保持、改善环境等方面的作用显著。国内外学者对该树种的生理生态特征、树轮生长特性等方面进行了广泛研究[15-17]。但是,关于青海云杉单木生长模型的相关研究较少[2]。因此,本研究以青海云杉为研究对象,分析不同起源、不同坡位下单木胸径生长的规律,构建单木胸径生长模型和考虑起源和坡位影响的混合效应模型,旨在为今后有效保护、合理经营青海云杉林提供经验模型和有益参考。
-
采样点位于青海省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。
采样点
Sample region起源
Origin坡位
Slope position样本量
Sample size树轮宽度
Tree-ring width/mm时段
Period尖扎县 天然 上坡 16 5.25/1.86/0.25 1933—2022 中坡 3 9.57/1.53/0.21 1927—2022 下坡 13 3.82/1.72/0.17 1941—2022 人工 上坡 15 5.65/1.60/0.14 1969—2022 下坡 13 5.00/1.90/0.28 1969—2022 祁连县 人工 上坡 14 2.78/0.90/0.13 1747—2022 中坡 15 3.22/0.80/0.10 1833—2022 下坡 15 3.36/1.18/0.21 1824—2022 大通县 天然 上坡 3 3.39/1.76/0.16 1947—2022 中坡 6 5.24/1.90/0.16 1950—2022 下坡 9 5.54/2.13/0.37 1950—2022 平坡 24 5.93/1.90/0.21 1941—2022 人工 平坡 16 7.13/3.11/0.49 1965—2022 注:树轮宽度表示为最大轮宽/平均轮宽/最小轮宽。
Note: Tree-ring width present as maximum tree-ring width/ mean tree-ring width/ minimum tree-ring width.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。树龄为树芯的年轮数加上林木生长到胸径位置时的年限。
-
结合树木的生长规律及前人的研究成果,本研究选取理查德(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包进行混合效应模型的参数估计。 -
为了对模型的拟合结果进行评价和比较,采用的指标有决定系数(R2)和均方根误差(RMSE)。采用全部数据计算模型估计值的精度指标,选择总体相对误差(TRE)、平均系统误差(MSE)、平均预估误差(MPE)和平均百分标准误差(MPSE)4项误差指标对模型进行综合评价。根据R2较大, 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。
天然林不同坡位中,青海云杉在上坡的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模型的拟合结果表现最好,R2为0.915;对于人工林,Korf模型的拟合效果最优,R2为0.946。2个模型的TRE在 ± 1%以内,MSE在 ± 2%以内,MPE均在3%以内,MPSE在35%以内。
起源
Origin模型
Model参数预估值 Parameter estimates 评价指标 Evaluation indices a b c R2 RMSE TRE/% MSE/% MPE/% MPSE/% 天然 Natural Gompertz 24.25 2.888 8 0.056 3 0.915 2.32 −0.09 −0.10 2.39 29.72 人工 Planted Korf 228.75 6.374 4 0.202 6 0.946 1.81 0.48 1.22 1.92 34.76 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%范围内。
起源
Origin坡位
Slope position模型
Model参数预估值 Parameter estimates 评价指标 Evaluation indices a b c R2 RMSE TRE/% MSE/% MPE/% MPSE/% 天然 Natural 上坡 Upper slope Richards 56.54 1.018 5 0.008 3 0.982 1.18 0.19 −0.10 4.12 26.84 中坡 Middle slope Logistic 20.28 8.330 2 0.108 6 0.915 1.84 1.16 0.56 5.76 27.12 下坡 Lower slope Gompertz 24.07 3.104 2 0.054 8 0.913 2.43 0.51 0.32 5.09 35.28 平坡 Flat slope Korf 938.03 7.604 9 0.179 5 0.996 0.48 0.58 0.31 3.75 23.19 人工 Planted 上坡 Upper slope Korf 3 020.36 8.385 8 0.104 4 0.940 1.88 1.18 1.95 2.85 33.90 中坡 Middle slope Richards 27.40 0.957 6 0.008 6 0.973 1.05 0.29 0.34 2.90 25.05 下坡 Lower slope Korf 15 768.94 10.470 7 0.103 6 0.980 1.51 0.90 1.72 2.55 28.41 平坡 Flat slope Logistic 16.42 17.207 7 0.233 5 0.955 1.18 −0.41 −0.28 4.40 19.90 Table 3. The parameter estimates and evaluation indices of optimal individual-tree DBH growth models 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。
模型
Model混合效应参数
Mixed-effects parameters异方差函数
Heteroscedasticity functionAIC BIC LogLik Logistic 无 无 13 742.680 13 765.520 −6 867.341 b、c 无 12 798.950 12 856.050 −6 389.477 Korf 无 无 13 657.260 13 680.090 −6 824.628 b 无 12 675.170 12 709.420 −6 331.584 Gompertz 无 无 13 709.540 13 732.380 −6 850.770 a、c 无 12 700.520 12 757.610 −6 340.259 常数加幂函数 12 396.490 12 465.000 −6 186.246 幂函数 12 394.490 12 457.300 −6 186.247 指数函数 不收敛 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%。
项目
Items参数
Parameter基础模型
Based model混合效应模型
Mixed effect model固定参数 a 20.24(p<0.000 1) 23.45(p<0.000 1) b 2.194 4(p<0.000 1) 2.558 4(p<0.000 1) c 0.044 8(p<0.000 1) 0.046 7(p<0.000 1) 异方差函数 幂函数 0.550 1 拟合统计量 R2 0.529 0.702 RMSE 5.23 4.14 检验指标 TRE/% −0.28 0.03 MSE/% −1.82 −0.30 MPE/% 5.33 4.23 MPSE/% 34.79 29.54 Table 5. Parameter estimation and Goodness-of-fit of based model and mixed effect model
-
利用最优模型对青海云杉胸径的生长拟合见图5。由图可以看出,青海云杉胸径总生长量随年龄逐渐增大,其中天然林上坡和人工林下坡的单木胸径生长趋势较好。天然林整体、天然林中坡和天然林下坡的单木连年生长量CAI先增大后减小,平均生长量MAI逐年下降;天然林上坡与平坡的单木CAI和MAI均呈现逐年下降趋势。除了人工林平坡,人工林的单木CAI和MAI均逐年下降。
Diameter Growth Models of Picea crassifolia from Different Habitats in Qinghai Province
- Received Date: 2023-03-13
- Accepted Date: 2023-07-18
- Available Online: 2024-02-20
Abstract: