Development of climate-sensitive nonlinear mixed-effects tree height-DBH model for Cunninghamia lanceolata
-
摘要:目的 建立基于林分优势高和气候因子的杉木树高−胸径非线性混合效应模型,为杉木生长研究和经营管理提供理论依据。方法 基于2020年国家森林资源年度监测评价广西壮族自治区试点25块杉木样地的每木胸径和树高实测数据及样地位置对应的气候数据,选择7个常用的树高−胸径模型,筛选出模拟精度最高的模型作为基础模型,再引入代表林分竞争、立地条件和气候因子的变量构建广义非线性模型,并在此基础上,加入样地效应构建杉木非线性混合效应模型。最后,运用十折交叉验证法对3种模型进行检验。结果 Chapman-Richards模型为最佳杉木树高−胸径关系基础模型,林分优势高、林分断面积和年平均降水量与树高生长显著相关,用于构建广义非线性模型,对比分析确定随机参数为3个的组合构造非线性混合效应模型。基础模型、广义非线性模型、非线性混合效应模型的调整决定系数分别为0.674 2、0.797 3和0.857 3,平均绝对误差分别为1.607 5、1.270 1和1.010 6 m,均方根误差分别为2.032 1、1.632 1和1.338 4 m,相对均方根误差分别为20.796 4%、16.703 3%和13.697 3%,混合效应模型呈现出更好的拟合效果。结论 引入林分优势高和气候因子的杉木树高−胸径非线性混合效应模型可以较好地描述杉木树高胸径曲线,适用于大范围的树高预测。Abstract:Objective The nonlinear mixed-effects tree height-DBH model of Cunninghamia lanceolata based on stand dominant height and climate factors is established, which provides theoretical basis for the research on growth and forest management.Method Based on the annual monitoring and evaluation of national forest resources of Guangxi Zhuang Autonomous Region, southern China in 2020, this study used the data of DBH and height of each tree, climate data of 25 Cunninghamia lanceolata sample plots, chose the basic model with the highest simulation accuracy among seven common height-DBH models. On this basis, stand competition, site condition and climatic factors were used to build generalized nonlinear model, then used the sample plot effect to build nonlinear mixed-effects model. The 10-fold cross-validation method was applied to the test of three models.Result Chapman-Richards model was the basic height-DBH model with the highest accuracy. The stand dominant height, basal area of forest stands and the mean annual precipitation were significantly related to the tree height growth, which were used to build generalized nonlinear model. Through comparative analysis, the study selected three random parameters to build nonlinear mixed-effects model. The adjustment determination coefficient of basic model, generalized nonlinear model and nonlinear mixed-effects model were 0.674 2, 0.797 3 and 0.857 3, respectively, the mean absolute errors were 1.607 5, 1.270 1 and 1.010 6 m, the root-mean-square errors were 2.032 1, 1.632 1 and 1.338 4 m, and the relative root mean square errors were 20.796 4%, 16.703 3% and 13.697 3%, respectively. The nonlinear mixed-effects model showed the best fitting effect.Conclusion Using nonlinear mixed-effects tree height-DBH model based on stand dominant height and climatic factors can better describe the height-DBH curve of Cunninghamia lanceolata, which is suitable for the prediction of tree height on a large scale.
-
胸径和树高作为森林资源调查监测的最基本测量因子,是计算森林蓄积量、生物量和碳储量,衡量林分生长潜力和立地质量的重要因子[1]。胸径能快速准确获取,而树高测定相对困难且不太精确。生产实践中,往往通过样木胸径数据与抽样样地样木树高全测数据,或者利用几株林分平均木的树高数据,构建树高−胸径回归模型,再推算林分总体其他样木树高。传统的Richards、Logistic、Korf和Gompertz等回归模型[2−6],能很好地描述树高−胸径生长的生物学特征,但林木生长处于林木竞争、立地等复杂环境中,不同样地同一树种的树高与胸径的关系不尽相同,单一的生长曲线难以描述所有可能的树高与胸径的关系[7]。对此,众多专家学者逐渐将反映立地条件、林分密度、森林类型等因子纳入到树高与胸径关系研究中[8−11],BP神经网络、分位数回归、哑变量模型和混合效应模型[12−14],以及树高分级方法[15−17]不断创新,符利勇[18]通过建立非线性混合效应模型,解决了树高−胸径模型在个体差异与样地方面的问题。董云飞[19]研究发现传统非线性模型、BP神经网络模型、非线性混合模型3种标准树高曲线建立方法,拟合精度和预测效果逐步有所提高。基于上述方法所做的研究大多针对较小尺度,模型精度得到了一定的提高,然而推广应用到省级乃至全国大区域,由于树种的立地与土壤、林分特征、气候差异等因子难以全面准确掌握,部分方法的拟合效果就不甚理想。
林木生长不仅受自身条件、竞争和立地等因子影响,还与光照、温度和降水等气候因子密切相关[9,20−21],其径向生长对气候因子的响应十分敏感[22−24]。然而,传统的树高胸径模型未考虑气候因子,没有客观全面反映气候环境对林木生长影响。近年来国内外学者开始聚焦于气候因子对树高−胸径曲线的影响[25−28],但极少在大尺度上分析气候因素的影响,本研究将省域范围的气候因子作为自变量引入树高胸径模型。
杉木(Cunninghamia lanceolata)是我国南方种植广泛的主要用材树种之一,《2021中国林草资源及生态状况》结果显示[29],在乔木林主要优势树种(组)中,杉木林1 240.27 × 104 hm2,占全国乔木林面积的6.21%,蓄积113 667.90 × 104 m3,占全国乔木林蓄积的5.99%,在我国的木材供给、水土保持、固碳增汇等方面起到了十分重要的作用。杉木为亚热带树种,广西壮族自治区是其主产区,杉木蓄积居全国31个省级单位第一。
本研究基于国家森林资源年度监测评价广西壮族自治区试点杉木样地的每木胸径和树高实测数据,从7种常用的树高胸径关系理论方程筛选出最优模型作为基础模型,加入林分优势高、林分每公顷断面积和气候因子等协变量建立广义非线性模型,并在此基础上引入样地随机效应构建了杉木树高−胸径非线性混合效应模型,比较3种不同模型拟合杉木树高胸径关系的效果,以期为该区域杉木生长研究和经营管理提供理论依据和数据支持。
1. 研究区概况
研究区地处中国南部,涉及广西壮族自治区全区(104°28′ ~ 112°04′E,20°54′ ~ 26°23′N),北回归线横贯全区中部,位于云贵高原向东南沿海丘陵过渡地带,地势总貌为西北高东南低,由西北向东南倾斜,四周山地环绕。属于亚热带季风气候区,年平均气温20.7 ℃,最热月为7月,月均气温23 ~ 29 ℃,年日照时数1 519.1 h,年均降雨量1 542.5 mm。水热条件良好,林木生长量是全国平均水平的2 ~ 3倍。全区人工林主要树种有桉树(Eucalyptus robusta)、杉木和马尾松(Pinus massoniana)等,面积869.36万 hm2,居全国第一;天然林主要以常绿落叶阔叶混交林为主。
2. 材料与方法
2.1 数据来源
研究数据来源于2020年国家森林资源年度监测评价广西壮族自治区试点项目,在全区森林资源连续清查杉木固定样地中随机筛选25块纯林样地,样地面积为0.066 7 hm2,样地分布如图1所示。于2020年11—12月对每个样地进行每木检尺并测量树高,共有杉木3 098株,记录样地坐标、海拔、坡度、坡向、土层厚度、腐殖质厚度等立地因子,树种、胸径、平均年龄等林木林分因子。以样地调查数据为基础,计算林分每公顷断面积、林分株数密度以及优势树种组里的最高6株样木的算术平均值作为林分优势高。
气候因子数据来源于Wang等[30]研发的ClimateAP软件。运用该软件逐个输入25块样地的坐标、海拔,根据林分平均年龄数据来提取对应时间段的年平均气温、最热月温度、年平均降水量等气候数据,再计算得出平均气候因子。如某一杉木样地平均年龄为28年,则提取此样地所在位置近28年的气候因子,取平均值即为该样地的气候因子,其他样地依次类推,从而获取所有样地的气候因子数据。样地调查数据和气候数据统计情况见表1。
表 1 样地因子和气候因子统计Table 1. Statistical table of survey factors and climatic variables in sample plots变量 Variable 变量符号
Variable symbol最小值
Min. value最大值
Max. value平均值
Mean标准差
SD胸径
Diameter at breast height/cmD 5.0 40.9 11.7 5.1 树高
Tree height/mH 3.0 21.8 9.8 3.6 林分断面积/(m2·hm−2)
Stand basal area/(m2·ha−1)SBA 5.61 52.88 28.11 12.49 林分密度/(株·hm−2)
Stand density/(tree·ha−1)NT 375 5 640 2 521 1 354 平均年龄/a
Average age/yearA 4 43 14 8 林分优势高
Stand dominant height/mHdom 9.8 21.0 14.7 2.4 年平均气温
Annual mean temperature/℃MAT 13.50 23.50 18.28 2.25 最热月平均气温
Mean temperature of the warmest month/℃MWMT 23.50 29.50 27.06 1.66 最冷月平均气温
Mean temperature of the coldest month/℃MCMT −1.10 17.40 5.34 4.14 平均温差
Mean temperature difference/℃TD 8.80 21.50 18.26 3.20 年平均降水量
Annual mean precipitation/mmPMA 486.00 2 083.00 1 592.84 374.27 湿热指数
HumidexAHM 11.40 62.00 20.43 11.02 0 ℃以下天数
Days below 0 ℃DD_0 0 25 5 8 5 ℃以上天数
Days above 5 ℃DD5 3 271 6 548 4 840 752 18 ℃以下天数
Days below 18 ℃DD_18 116 2 131 1 126 489 18 ℃以上天数
Days above 18 ℃DD18 482 1 990 1 221 358 2.2 模型构建
2.2.1 基础模型的构建
根据以往研究[31−32],非线性模型较线性模型能更好描述树高胸径关系。本文选取生态学意义较好的7个非线性树高−胸径模型,包括经典的Richards、Logistic、Korf和Gompertz模型(表2),从中筛选预测精度最高的作为杉木树高−胸径关系的基础模型。
表 2 树高−胸径候选模型Table 2. Candidate models of tree height-DBH模型名称
Model name模型公式
Model formula公式序号
Equation No.Curtis H=1.3+a1[D/(1+D)]a2+ε (1) Meyer H=1.3+a1[1−exp(−a2D)]+ε (2) Wykoff H=1.3+exp[a1+a2/(D+1)]+ε (3) Chapman-Richards H=1.3+a1[1−exp(−a2D)]a3+ε (4) Logistic H=1.3+a1/[1+a2exp(−a3D)]+ε (5) Weibull H=1.3+a1[1−exp(−a2Da3)]+ε (6) Gompertz H=1.3+a1exp[−a2exp(−a3D)]+ε (7) 注:a1、a2、a3为模型系数,ε为误差项。Notes: a1, a2, a3 are the model coefficients, ε is the error term. 2.2.2 广义非线性模型的构建
树高胸径的关系受树种自身条件、竞争、立地质量和气候等多重因子综合影响,为进一步提高模型预测精度,将林分平均年龄(A),包括林分断面积(SBA)、林分密度(NT)在内的竞争因子,以林分优势高(Hdom)为代表的立地因子,以及年平均气温(TMA)、平均温差(TD)和年平均降水量(PMA)等在内的气候因子与树高进行相关性分析,选择与树高相关性强的变量以不同的组合形式引入筛选得到的最优基础模型进行回归系数显著性检验和共线性检验确定最终变量构建广义非线性模型。
2.2.3 混合效应模型的构建
采用非线性混合效应模型方法构建杉木树高−胸径模型,非线性混合效应模型形式如下所示。
\left\{ {\begin{array}{*{20}{c}} {{H_{ij}} = f\left( {\beta ,{b_i},{D_{ij}}} \right) + {{\boldsymbol{\varepsilon}} _{ij}}}\\ {{b_i} \sim N(0,\;{\boldsymbol{G}})}\\ {{{\boldsymbol{\varepsilon}} _{ij}} \sim N\left( {0,{{\boldsymbol{R}}_{ij}}} \right)} \end{array}} \right. 式中:Hij和Dij分别为第i个样地第j株样木的树高和胸径,f(β,bi,Dij)描述胸径与树高之间的函数关系,β为固定效应参数估计值,bi为样地随机效应参数估计值,G为随机效应方差−协方差矩阵;εij为误差向量;Rij为样地内方差−协方差矩阵。
在构建混合效应模型时,将基础模型所有参数作为随机效应参数,通过比较不同随机参数组合模型的拟合精度指标确定最优模型。选取对角矩阵(DM)、复合对称矩阵(CS)、广义正定矩阵(general positive-definite matrix) 3种林业常用结构来确定最优随机效应方差−协方差结构。
2.3 模型的选择与检验
根据平均绝对误差(mean absolute error,MAE)、均方根误差(root mean square error,RMSE)、相对均方根误差(relative root mean square error,rRMSE)和调整决定系数(adjusted coefficient of determination,R2 adj)来筛选最优基础模型,将赤池信息准则(AIC)、贝叶斯准则(BIC)和对数似然值(Loglik)作为表示广义非线性模型、非线性混合效应模型的固定效应部分的拟合优度的指标。统计指标表达式如下所示。
{\rm{MAE = }}\frac{{\displaystyle \sum {\left| {{H_{ij}} - {{\hat H}_{ij}}} \right|} }}{n} (9) {\rm{RMSE}} = \sqrt {\frac{{\displaystyle \sum {{{\left( {{H_{ij}} - {{\hat H}_{ij}}} \right)}^2}} }}{n}} (10) {\rm{rRMSE}} = \frac{{{\rm{RMSE}}}}{{\overline {{H_{ij}}} }} \times 100\% (11) {R^2} = 1 - \frac{{\displaystyle \sum {{{\left( {{H_{ij}} - {{\hat H}_{ij}}} \right)}^2}} }}{{\displaystyle \sum {{{\left( {{H_{ij}} - \overline {{H_{ij}}} } \right)}^2}} }} (12) R_{{\rm{adj}}}^2 = 1 - \frac{{\left( {1 - {R^2}} \right)(n - 1)}}{{(n - k - 1)}} (13) {\rm{AIC}} = - 2\ln L + 2{{k}} (14) {\rm{BIC}} = - 2\ln L + {{k}}\ln n (15) {\rm{LRT}} = 2\left( {\ln {{L_1}} - \ln {{L_2}} } \right) (16) 式中: H_{ij} 为观察值, {\hat H_{ij}} 为估计值,\overline {{H_{ij}}} 为观察值的均值,n为样木株数,k为自变量的个数,L为模型的对数似然值,L1和L2分别为模型1和模型2的对数似然值。
混合效应模型的随机效应部分参数值\;{\mu _i} 检验需通过检验样本已知信息求得,使用Vonesh等提出的方法[33](式17)计算。采用指数函数(Ve)、幂函数(Vp)和常数加幂函数(Vc)形式来消除数据的异方差性。
{\mu _i} \approx {\boldsymbol{G}}{\boldsymbol{Z}}_i^{\rm{T}}{\left( {{{\boldsymbol{Z}}_i}{\boldsymbol{G}}{\boldsymbol{Z}}_i^{\rm{T}} + {{\boldsymbol{R}}_i}} \right)^{ - 1}}{{\boldsymbol{e}}_i} (17) V_{\rm{e}} = {\sigma ^2}\exp \left( {2\alpha {u_i}} \right) (18) V_{\rm{p}} = {\sigma ^2}\exp \left( {u_i^{2\alpha }} \right) (19) V_{\rm{c}} = {\sigma ^2}{(\alpha + u_i^\beta )^2} (20) 式中:Zi为已知设计矩阵,Ri为样地方差−协方差矩阵,ei为树高实测值与固定参数计算树高值之间的误差向量,{\sigma ^2} 为误差扩散的比例因子,α和β为函数的待估参数,u_i 基于固定效应参数计算得到的预测值。
3. 结果与分析
3.1 基础模型结果
对Curtis、Meyer、Wykoff和Chapman-Richards等7个树高−胸径模型进行拟合,比较各个模型的MAE、RMSE、rRMSE、R2 adj 3个指标值。结果表明,各模型的参数均表现为显著,7个模型都呈现较好的适应性,MAE值处于1.607 5 ~ 1.715 5 m,RMSE值处于2.032 1 ~ 2.102 9 m,rRMSE值处于20.796 4% ~ 21.521 4%,R2 adj处于0.573 0 ~ 0.674 2之间(表3)。根据MAE值、RMSE值和rRMSE值越小,R2 adj值越大的原则筛选最优基础模型,Chapman-Richards模型体现出最好的适应性,MAE值为1.607 5 m,RMSE值为2.032 1 m,rRMSE值为20.796 4%,R2 adj值0.674 2,选用Chapman-Richards模型作为杉木树高−胸径关系基础模型,模型具体形式如下所示。
表 3 基础模型拟合结果Table 3. Fitting results of basic models模型名称
Model name参数估计值 Parameter estimate MAE/m RMSE/m rRMSE/% R2 adj a1 a2 a3 Curtis 21.552 3*** 10.449 5*** 1.632 2 2.051 7 20.996 1 0.668 3 Meyer 22.296 5*** 0.043 3*** 1.715 5 2.102 9 21.521 4 0.573 0 Wykoff 3.106 3*** −11.399 5*** 1.634 6 2.052 2 21.002 4 0.664 3 Richards 15.478 3*** 0.124 7*** 1.988 9*** 1.607 5 2.032 1 20.796 4 0.674 2 Logistic 14.274 7*** 8.671 5*** 0.231 9*** 1.608 4 2.039 2 20.868 8 0.668 1 Weibull 14.798 3*** 0.022 9*** 1.516 9*** 1.612 9 2.039 7 20.874 3 0.673 1 Gompertz 14.954 7*** 3.171 4*** 0.158 7*** 1.612 4 2.039 2 20.869 3 0.671 3 注:***表示显著(P < 0.001)。MAE.平均绝对误差;RMSE.均方根误差;rRMSE.相对均方根误差;R2 adj.调整决定系数。Notes: *** means significant at the P < 0.001 level. MAE, mean absolute error; RMSE, root mean square error; rRMSE, relative root mean square error; R2 adj, adjusted coefficient of determination. H = 1.3 + 15.478\;3\left(1 - \exp {( - 0.124\;7D)}\right)^{1.988\;9} + \varepsilon (21) 3.2 广义非线性模型结果
依次将林分平均年龄、林分断面积、林分密度、林分优势高、年平均气温、平均温差和年平均降水量等变量与树高进行相关性分析。考虑到竞争、立地质量和气候因子内部的共线性,以及构建广义非线性模型的不宜过于复杂,分别在代表竞争、立地质量和气候条件的因子中选择与树高相关性最高的因子引入基础模型(PMA),最终确定林分断面积(SBA)、林分优势高(Hdom)、年平均降水量3个变量进入模型,3者与树高相关系数分别为0.415 2、0.481 0和−0.132 6,其中林分断面积代表竞争情况,林分优势高代表立地质量条件,年平均降水量代表气候条件。将3个变量以不同的组合方式代入最优基础模型Chapman-Richards模型,比较不同组合形式模型的MAE、RMSE和R2 adj值,确定最终的广义非线性模型,具体形式如下所示。
\begin{split} H = &1.3 + \left( {{a_1} + {a_4}{H_{{\rm{dom}}}}} \right)\times \\&{\left\{ {1 - \exp \left[ { - \left( {{a_2} + {a_5}S_{\rm{ BA}}} \right)D} \right]} \right\}^{\left( {{a_3} + {a_6}P_{\rm{MA }}} \right)}} + \varepsilon \end{split} (22) 式中:a1 ~ a6为模型系数。
广义非线性模型的AIC = 11 841.2,BIC = 11 883.5,Loglik = −5 913.63,MAE = 1.270 1 m,RMSE = 1.632 1 m,rRMSE = 16.703 3%,R2 adj = 0.797 3,从表4可以看出,各变量都显著影响树高的生长。
表 4 广义非线性模型参数估计结果Table 4. Results of parameter estimation of generalized nonlinear models参数
Parameter估计值
Estimate标准差
SDt P a1 6.226 0 0.452 2 13.766 < 0.000 1 a2 0.085 6 0.005 2 16.372 < 0.000 1 a3 1.507 0 0.086 0 17.514 < 0.000 1 a4 0.570 5 0.025 8 22.152 < 0.000 1 a5 0.001 1 0.000 1 13.432 < 0.000 1 a6 0.000 1 0.000 0 2.803 0.005 1 3.3 非线性混合效应模型结果
在广义非线性模型的基础上引入样地随机效应,对63个不同组合形式进行模拟,其中36个拟合收敛。当考虑一个随机参数时,分别添加在6个参数上都拟合收敛,在a2上显示拟合精度最高;2个随机参数添加时,15种组合都显示拟合收敛,在a2、a3上拟合精度最高,AIC和BIC分别为10 814.9和10 875.3;添加3个随机参数的各种组合中,20种组合有12个组合收敛,在a2、a3和a6上的拟合精度最高,AIC和BIC分别为10 811.8和10 890.3;添加3个以上随机参数时,22种组合仅3个收敛。不同随机参数个数的最优拟合精度AIC、BIC、Loglik结果显示如表5。添加了一定数量随机参数均比广义非线性基础模型的拟合效果好,随着添加随机参数数量增加,呈现AIC值逐步减少,Loglik值逐步增大的趋势。但随着参数数量的增加,拟合效果并没有显著提升,而且拟合不稳定性也增加,为防止过度拟合同时兼顾稳定性,本文选择在a2、a3和a6 3个参数上考虑随机效应。确定最后的非线性混合效应模型如下。
表 5 混合效应模型随机参数组合的部分拟合结果Table 5. Partial results of combinations of random parameters for mixed-effects model模型
Model随机参数
Random parameterAIC BIC Loglik LRT P 广义非线性模型
Generalized nonlinear model无 None 11 841.2 11 883.5 −5 913.63 模型1 Model 1 a2 10 904.4 10 952.7 −5 444.20 模型2 Model 2 a2、a3 10 814.9 10 875.3 −5 397.49 93.41 < 0.000 1 模型3 Model 3 a2、a3、a6 10 811.8 10 890.3 −5 392.92 9.14 0.027 5 注:模型2的LRT和P值为模型2与1相比较而得;模型3的LRT和P值为模型3与2相比较而得。Notes: the LRT and P in the line of model 2 are caululated by model 2 and 1. The LRT and P in the line of model 3 are caululated by model 3 and 2. \begin{split} H =& 1.3 + \left( {{a_1} + {a_4}{H_{{\rm{dom}}}}} \right)\{ 1 - \exp [ - \left( {{a_2} + {b_1}} \right) - \\ &{a_5}S_{\rm{ BA}}]D{\} ^{\left[ {\left( {{a_3} + {b_2}} \right) + \left( {{a_6} + {b_3}} \right)P_{\rm{MA}}} \right]}} + \varepsilon \end{split} (23) 式中:a1 ~ a6为固定参数,b1 ~ b3为随机参数。
在模型23上先后引入对角矩阵、复合对称矩阵、广义正定矩阵3种随机参数协方差结构,确定广义正定矩阵为最优随机效应方差−协方差结构后,再分别引入常数加幂函数、幂函数和指数函数3种异方差函数消除方差异质性,模型拟合结果如表6所示。根据结果可以看出加入幂函数后模型的AIC和BIC值最小,Loglik值最大,故将其选定为最终的异方差函数。最终得到的杉木树高−胸径非线性混合模型的表达式如下。
表 6 基于不同异方差函数的非线性混合效应模型比较Table 6. Comparison of nonlinear mixed-effects models based on different variance functions模型
Model异方差函数 Heteroscedasticity function AIC BIC Loglik LRT P 模型3
Model 3无 None 10 811.8 10 890.3 −5 392.92 模型3.1
Model 3.1常数加幂函数 Constant plus power function 10 760.2 10 850.8 −5 365.12 55.60 < 0.000 1 模型3.2
Model 3.2幂函数 Power function 10 758.2 10 842.7 −5 365.12 55.60 < 0.000 1 模型3.3
Model 3.3指数函数 Exponent function 10 782.0 10 866.5 −5 377.00 31.84 < 0.000 1 注:LRT和P值为模型3.1、3.2、3.3与模型3相比较而得。Note: LRT and P are caululated by comparison of model 3.1, 3.2, 3.3 and model 3. \begin{array}{l} H = 1.3 + \left( {5.355\;1 + 0.724\;0{H_{{\rm{dom}}}}} \right)\{ 1 - \exp [ - (0.065\;4 + \\ {b_1}) - 0.000\;9S_{\rm{BA}}]D{\} ^{\left[ {\left( {0.971\;3 + {b_2}} \right) + \left( {0.000\;3 + {b_3}} \right)P_{\rm{MA}}} \right]}} + \varepsilon \end{array} \begin{array}{l} {{\boldsymbol{b}}_i} = \left[ {\begin{array}{*{20}{l}} {{b_1}}\\ {{b_2}}\\ {{b_3}} \end{array}} \right]\sim {\boldsymbol{N}}\left\{ {\left[ {\begin{array}{*{20}{l}} 0\\ 0\\ 0 \end{array}} \right]} \right.,\;{\boldsymbol{G}} =\\ \left. { \left( {\begin{array}{*{20}{c}} {9.263\;6 \times {{10}^{ - 3}}} & {2.899\;7 \times {{10}^{ - 2}}} & { - 8.997\;7 \times {{10}^{ - 6}}}\\ {2.899\;7 \times {{10}^{ - 2}}} & {9.242\;5 \times {{10}^{ - 1}}} & { - 3.077\;7 \times {{10}^{ - 4}}}\\ { - 8.997\;7 \times {{10}^{ - 6}}} & { - 3.077\;7 \times {{10}^{ - 4}}} & {1.299\;7 \times {{10}^{ - 7}}} \end{array}} \right)} \right\}, \end{array} \begin{split} &{{\boldsymbol{e}}_i} \sim N\left( {0,{{\boldsymbol{R}}_i} = 0.397\;0{\boldsymbol{G}}_i^{0.5}{{\boldsymbol{\varGamma}} _i}{\boldsymbol{G}}_i^{0.5}} \right)\\&\qquad{\mathop{\rm var}} \left( {{{\boldsymbol{e}}_i}} \right) = 0.397\;0\hat y_{ij}^{0.675\;1} \end{split} (24) 式中:{\boldsymbol{\varGamma}} 为误差自相关矩阵。
3.4 模型比较
从树高−胸径基础模型,到广义非线性模型、非线性混合效应模型,随着竞争、立地和气候因子的加入,以及考虑样地随机效应后,拟合精度逐步提高,各模型拟合情况如表7所示。AIC值从基础模型的13 223.7减少到非线性混合效应的模型的10 758.2,BIC值从13 247.8减少到10 842.7,MAE、RMSE和rRMSE值也不同程度减少,R2 adj值从0.674 2增加到0.857 3。从不同模型的残差图(图2)可以看出,引入竞争、立地和气候因子的广义非线性模型,以及非线性混合效应模型,残差在同质性和正态性方面得到显著改善。
表 7 不同模型比较Table 7. Comparison of different models模型 Model AIC BIC Loglik MAE/m RMSE/m rRMSE/% R2 adj 基础模型
Basic model13 223.7 13 247.8 −6 607.84 1.607 5 2.032 1 20.796 4 0.674 2 广义非线性模型 Generalized nonlinear model 11 841.2 11 883.5 −5 913.63 1.270 1 1.632 1 16.703 3 0.797 3 混合效应模型
Mixed-effects model1 0758.2 10 842.7 −5365.12 1.010 6 1.338 4 13.697 3 0.857 3 3.5 模型检验
采用十折交叉验证的方法对3种模型精度进行验证,基础模型的MAE、RMSE、rRMSE和R2 adj平均值和标准差分别为(1.611 1 ± 0.046 3)m、(2.036 8 ± 0.043 5) m、(20.830 5 ± 0.529 7)%和(0.674 5 ± 0.021 4),广义非线性模型的MAE、RMSE、rRMSE和R2 adj平均值和标准差分别为(1.270 1 ± 0.047 6) m、(1.627 2 ± 0.040 9) m、(16.640 9 ± 0.511 6)%和(0.787 9 ± 0.016 7),非线性混合效应模型的MAE、RMSE、rRMSE和R2 adj平均值和标准差分别为(1.010 5 ± 0.039 0) m、(1.328 9 ± 0.047 5) m、(13.587 8 ± 0.470 8)%和(0.858 3 ± 0.015 8),3种模型的拟合精度从高到低依次为非线性混合效应模型、广义非线性模型、基础模型,与建模数据研究结果一致。
4. 讨 论
本研究利用广西壮族自治区的杉木林样地数据,引入林分优势高和气候因子,建立了非线性混合效应模型。通过筛选常见7个预测能力强、生物学意义较好的树高−胸径模型,得到Chapman-Richards为最优基础模型。通过分析影响树高生长的竞争、立地、气候等因子,确定林分优势高、林分断面积和年平均降水量3个变量进入模型。林分优势高相对于海拔、坡度、坡向、土壤厚度、腐殖质厚度等立地条件因子,其具备较强代表性,在同一树种不同地域更能直接反应立地好坏,且在生产实践中更易获取,近年来受到了很多研究林分生长方向的学者青睐。Sharma等[12]以林分优势高、林分断面积和林分密度代表立地特征构建加拿大安大略省冷杉(Abies fabri)、红松(Pinus koraiensis)、白桦(Betula platyphylla)等树种的树高−胸径混合效应模型,相对于传统Chapman-Richards模型取得了更好的拟合效果;郭嘉等[17]在研究秦岭林区松栎林引入优势木树高作为立地变量建立树高−胸径模型,有效提高了模型拟合精度,本研究得出了类似的结论。
气候变化深刻影响树木生长[24,34−35],在树木年轮气候学上大量研究表明树木径向生长与气候因子存在一定的规律。本研究经相关性分析得到杉木树高与年平均气温、湿热指数等因子正相关但不显著,与年平均降水量、平均温差等因子负相关,基础模型中引入年平均降水量有利于提高拟合精度。杉木生长习性为喜温暖湿润,不耐严寒及湿热,水湿条件的影响大于温度条件。杉木研究样地地域广阔但总体位于北回归线附近,属亚热带季风气候区,是杉木的温度适宜区,此区域温度不是其树高生长主要影响因素,但沿海较强的降雨易形成湿热天气,不利于养分积累而影响杉木生长,表现与树高呈现负相关的情况。Wang等[36]在研究基于气象变量的红松人工林单木胸径生长模型时发现,生长季最大降水量明显影响红松胸径生长,与其负相关,认为过于湿润的气候条件将抑制林木生长,与本文得到一致的结论。但是,树种的遗传特征主要决定了树木生长对温度和降水存在不同的响应[37],适地适树地分析气候因子对于树木生长很有必要。Tian等[38]研究发现年平均气温大小对于五角槭(Acer mono)的树高生长影响微不足道,但对冷杉和胡桃楸(Juglans mandshurica)的树高生长明显影响,年平均降水量显著影响蒙古栎(Quercus mongolica)的树高−胸径异速生长,对于五角槭的异速生长影响却很小,树种的特异性导致对环境因子表现出不同程度的敏感性反应。
混合效应模型相对于传统的基础模型和广义非线性模型,具有更强灵活性,拟合结果来看,非线性混合效应模型优于广义非线性模型和基础模型,平均绝对误差MAE和均方根误差RMSE值较基础模型分别减少了0.667 6 m和0.779 6 m,R2 adj增加了0.208 4。混合效应能兼顾样地的固定效应和随机效应,采用多种异方差函数和自相关结构可以有效改善模型误差的方差异质性和数据自相关性,适用于区域大尺度重复观测样地的林业数据。但要深入理解各变量的生物学意义和各变量间的共线性,准确把握基础模型形式参数的构造[39−41]。本研究计算所有可能的63个随机效应组合后,也尝试用5个随机变量进入模型,但形式参数个数较多且同时考虑多个水平逐级嵌套随机效应使得模型计算收敛困难,排除未收敛的组合,对比筛选出随机变量数量为3的最佳形式参数组合。
5. 结 论
本研究将经典Chapman-Richards模型、广义非线性模型、非线性混合模型应用于杉木树高−胸径曲线中,先后建立起3种树高−胸径模型,结果显示3种模型的拟合精度逐步提高,非线性混合效应模型效果最好。近年来国家启动林草生态综合监测工作,以国家森林资源连续清查固定样地框架为基础,覆盖全国各省的样地坐标、海拔信息能提取丰富的气候因子数据,并且样地调查将林分平均优势高纳入年度调查因子,为构建国家和地方尺度的树高−胸径模型提供了数据基础,未来可尝试运用非线性混合效应模型开展宽区域、多树种、更精确的数表研建。
-
表 1 样地因子和气候因子统计
Table 1 Statistical table of survey factors and climatic variables in sample plots
变量 Variable 变量符号
Variable symbol最小值
Min. value最大值
Max. value平均值
Mean标准差
SD胸径
Diameter at breast height/cmD 5.0 40.9 11.7 5.1 树高
Tree height/mH 3.0 21.8 9.8 3.6 林分断面积/(m2·hm−2)
Stand basal area/(m2·ha−1)SBA 5.61 52.88 28.11 12.49 林分密度/(株·hm−2)
Stand density/(tree·ha−1)NT 375 5 640 2 521 1 354 平均年龄/a
Average age/yearA 4 43 14 8 林分优势高
Stand dominant height/mHdom 9.8 21.0 14.7 2.4 年平均气温
Annual mean temperature/℃MAT 13.50 23.50 18.28 2.25 最热月平均气温
Mean temperature of the warmest month/℃MWMT 23.50 29.50 27.06 1.66 最冷月平均气温
Mean temperature of the coldest month/℃MCMT −1.10 17.40 5.34 4.14 平均温差
Mean temperature difference/℃TD 8.80 21.50 18.26 3.20 年平均降水量
Annual mean precipitation/mmPMA 486.00 2 083.00 1 592.84 374.27 湿热指数
HumidexAHM 11.40 62.00 20.43 11.02 0 ℃以下天数
Days below 0 ℃DD_0 0 25 5 8 5 ℃以上天数
Days above 5 ℃DD5 3 271 6 548 4 840 752 18 ℃以下天数
Days below 18 ℃DD_18 116 2 131 1 126 489 18 ℃以上天数
Days above 18 ℃DD18 482 1 990 1 221 358 表 2 树高−胸径候选模型
Table 2 Candidate models of tree height-DBH
模型名称
Model name模型公式
Model formula公式序号
Equation No.Curtis H = 1.3 + {a_1}{\left[ {D/\left( {1 + D} \right)} \right]^{{a_2}}} + \varepsilon (1) Meyer H = 1.3 + {a_1}\left[ {1 - \exp \left( { - {a_2}D} \right)} \right] + \varepsilon (2) Wykoff H = 1.3 + \exp \left[ {{a_1} + {a_2}/(D + 1)} \right] + \varepsilon (3) Chapman-Richards H = 1.3 + {a_1}{\left[ {1 - \exp \left( { - {a_2}D} \right)} \right]^{{a_3}}} + \varepsilon (4) Logistic H = 1.3 + {a_1}/\left[ {1 + {a_2}\exp \left( { - {a_3}D} \right)} \right] + \varepsilon (5) Weibull H = 1.3 + {a_1}\left[ {1 - \exp \left( { - {a_2}{D^{{a_3}}}} \right)} \right] + \varepsilon (6) Gompertz H = 1.3 + {a_1}\exp \left[ { - {a_2}\exp \left( { - {a_3}D} \right)} \right] + \varepsilon (7) 注:a1、a2、a3为模型系数,ε为误差项。Notes: a1, a2, a3 are the model coefficients, ε is the error term. 表 3 基础模型拟合结果
Table 3 Fitting results of basic models
模型名称
Model name参数估计值 Parameter estimate MAE/m RMSE/m rRMSE/% R2 adj a1 a2 a3 Curtis 21.552 3*** 10.449 5*** 1.632 2 2.051 7 20.996 1 0.668 3 Meyer 22.296 5*** 0.043 3*** 1.715 5 2.102 9 21.521 4 0.573 0 Wykoff 3.106 3*** −11.399 5*** 1.634 6 2.052 2 21.002 4 0.664 3 Richards 15.478 3*** 0.124 7*** 1.988 9*** 1.607 5 2.032 1 20.796 4 0.674 2 Logistic 14.274 7*** 8.671 5*** 0.231 9*** 1.608 4 2.039 2 20.868 8 0.668 1 Weibull 14.798 3*** 0.022 9*** 1.516 9*** 1.612 9 2.039 7 20.874 3 0.673 1 Gompertz 14.954 7*** 3.171 4*** 0.158 7*** 1.612 4 2.039 2 20.869 3 0.671 3 注:***表示显著(P < 0.001)。MAE.平均绝对误差;RMSE.均方根误差;rRMSE.相对均方根误差;R2 adj.调整决定系数。Notes: *** means significant at the P < 0.001 level. MAE, mean absolute error; RMSE, root mean square error; rRMSE, relative root mean square error; R2 adj, adjusted coefficient of determination. 表 4 广义非线性模型参数估计结果
Table 4 Results of parameter estimation of generalized nonlinear models
参数
Parameter估计值
Estimate标准差
SDt P a1 6.226 0 0.452 2 13.766 < 0.000 1 a2 0.085 6 0.005 2 16.372 < 0.000 1 a3 1.507 0 0.086 0 17.514 < 0.000 1 a4 0.570 5 0.025 8 22.152 < 0.000 1 a5 0.001 1 0.000 1 13.432 < 0.000 1 a6 0.000 1 0.000 0 2.803 0.005 1 表 5 混合效应模型随机参数组合的部分拟合结果
Table 5 Partial results of combinations of random parameters for mixed-effects model
模型
Model随机参数
Random parameterAIC BIC Loglik LRT P 广义非线性模型
Generalized nonlinear model无 None 11 841.2 11 883.5 −5 913.63 模型1 Model 1 a2 10 904.4 10 952.7 −5 444.20 模型2 Model 2 a2、a3 10 814.9 10 875.3 −5 397.49 93.41 < 0.000 1 模型3 Model 3 a2、a3、a6 10 811.8 10 890.3 −5 392.92 9.14 0.027 5 注:模型2的LRT和P值为模型2与1相比较而得;模型3的LRT和P值为模型3与2相比较而得。Notes: the LRT and P in the line of model 2 are caululated by model 2 and 1. The LRT and P in the line of model 3 are caululated by model 3 and 2. 表 6 基于不同异方差函数的非线性混合效应模型比较
Table 6 Comparison of nonlinear mixed-effects models based on different variance functions
模型
Model异方差函数 Heteroscedasticity function AIC BIC Loglik LRT P 模型3
Model 3无 None 10 811.8 10 890.3 −5 392.92 模型3.1
Model 3.1常数加幂函数 Constant plus power function 10 760.2 10 850.8 −5 365.12 55.60 < 0.000 1 模型3.2
Model 3.2幂函数 Power function 10 758.2 10 842.7 −5 365.12 55.60 < 0.000 1 模型3.3
Model 3.3指数函数 Exponent function 10 782.0 10 866.5 −5 377.00 31.84 < 0.000 1 注:LRT和P值为模型3.1、3.2、3.3与模型3相比较而得。Note: LRT and P are caululated by comparison of model 3.1, 3.2, 3.3 and model 3. 表 7 不同模型比较
Table 7 Comparison of different models
模型 Model AIC BIC Loglik MAE/m RMSE/m rRMSE/% R2 adj 基础模型
Basic model13 223.7 13 247.8 −6 607.84 1.607 5 2.032 1 20.796 4 0.674 2 广义非线性模型 Generalized nonlinear model 11 841.2 11 883.5 −5 913.63 1.270 1 1.632 1 16.703 3 0.797 3 混合效应模型
Mixed-effects model1 0758.2 10 842.7 −5365.12 1.010 6 1.338 4 13.697 3 0.857 3 -
[1] 邓祥鹏, 许芳泽, 赵善超, 等. 基于贝叶斯法的新疆天山云杉树高−胸径模型研究[J]. 北京林业大学学报, 2023, 45(1): 11−20. Deng X P, Xu F Z, Zhao S C, et al. Tree height-DBH model for Picea schrenkiana in Tianshan Mountain, Xinjiang of northwestern China based on Bayesian method[J]. Journal of Beijing Forestry University, 2023, 45(1): 11−20.
[2] Curtis R O. Height-diameter and height-diameter-age equations for second-growth Douglas-fir[J]. Forestry Science, 1967, 13(4): 365−375.
[3] Fang Z X, Bailey R L. Height-diameter models for tropical forests on Hainan Island in southern China[J]. Forest Ecology and Management, 1998, 110: 315−327. doi: 10.1016/S0378-1127(98)00297-7
[4] 王明亮, 唐守正. 标准树高曲线的研制[J]. 林业科学研究, 1997, 10(4): 259−264. Wang M L, Tang S Z. Research on universal height-diameter curves[J]. Forest Research, 1997, 10(4): 259−264.
[5] 骆期邦, 曾伟生, 彭长清. 可变参数相对树高曲线模型及其应用研究[J]. 林业科学, 1997, 33(3): 202−210. Luo Q B, Zeng W S, Peng C Q. Variable relative tree height curve model and its application in tree volume estimation[J]. Scientia Silvae Sinicae, 1997, 33(3): 202−210.
[6] 赵俊卉, 亢新刚, 刘燕. 长白山主要针叶树种最优树高曲线研究[J]. 北京林业大学学报, 2009, 31(4): 13−18. Zhao J H, Kang X G, Liu Y. Optimal height-diameter models for dominant coniferous species in Changbai Mountain, northeastern China[J]. Journal of Beijing Forestry University, 2009, 31(4): 13−18.
[7] Temesgen H, Gadow K V. Generalized height-diameter models: an application for major tree species in complex stands of interior British Columbia[J]. European Journal of Forest Research, 2004, 123: 45−51. doi: 10.1007/s10342-004-0020-z
[8] 赵俊卉, 亢新刚, 张慧东, 等. 长白山3个主要针叶树种的标准树高曲线[J]. 林业科学, 2010, 46(10): 191−194. Zhao J H, Kang X G, Zhang H D, et al. Generalized height-diameter models for three main coniferous trees species in Changbai Mountain[J]. Scientia Silvae Sinicae, 2010, 46(10): 191−194.
[9] 赵俊卉, 亢新刚, 张慧东, 等. 长白山主要针叶树种胸径和树高变异系数与竞争因子的关系[J]. 应用生态学报, 2009, 20(8): 1832−1837. Zhao J H, Kang X G, Zhang H D, et al. Relationships between coefficient of variation of diameter and height competition index of main coniferous trees in Changbai Mountains[J]. Chinese Journal of Applied Ecology, 2009, 20(8): 1832−1837.
[10] 蒋益, 邓华锋, 高东启, 等. 用度量误差模型方法建立油松树高曲线方程组[J]. 东北林业大学学报, 2015, 43(5): 126−129. Jiang Y, Deng H F, Gao D Q, et al. Constructing height-diameter curve equations with measurement error models for Chinese pine stands[J]. Journal of Northeast Forestry University, 2015, 43(5): 126−129.
[11] 娄明华, 张会儒, 雷相东, 等. 基于空间自相关的天然蒙古栎阔叶混交林林木胸径−树高模型[J]. 林业科学, 2017, 53(6): 67−76. Lou M H, Zhang H R, Lei X D, et al. Individual diameter-height models for mixed Quercus mongolica broadleaved natural stands based on spatial autocorrelation[J]. Scientia Silvae Sinice, 2017, 53(6): 67−76.
[12] Sharma M, Parton J. Height-diameter equations for boreal tree species in Ontario using a mixed-effects modeling approach[J]. Forest Ecology and Management, 2007, 249(3): 187−198. doi: 10.1016/j.foreco.2007.05.006
[13] Felipe C, Margarida T, Paula S, et al. A generalized nonlinear mixed-effects height-diameter model for Eucalyptus globulus L. in northwestern Spain[J]. Forest Ecology and Management, 2009, 259: 943−952.
[14] Karol B, Lauri M, Paula S. Mixed-effects generalized height-diameter model for young silver birch stands on post-agricultural lands[J]. Forest Ecology and Management, 2020, 460: 117901. doi: 10.1016/j.foreco.2020.117901
[15] 李杰. 基于分级的福建将乐地区栲树树高曲线模型研究[J]. 西北农林科技大学学报, 2019, 47(11): 34−42. Li J. Classification based height-diameter model for Castanopsis fargesii in Jiangle, Fujian[J]. Journal of Northwest A&F University, 2019, 47(11): 34−42.
[16] 李海奎, 法蕾. 基于分级的全国主要树种树高−胸径曲线模型[J]. 林业科学, 2011, 47(10): 83−90. Li H K, Fa L. Height-diameter model for major tree species in China using the classified height method[J]. Scientia Silvae Sinicae, 2011, 47(10): 83−90.
[17] 郭嘉, 孙帅超, 田相林, 等. 引入优势木树高建立的秦岭林区松栎林树高−胸径模型[J]. 东北林业大学学报, 2019, 47(11): 66−72. Guo J, Sun S C, Tian X L, et al. Predicting tree height from tree diameter and dominant height for pine-oak forests in Qinling Mountains[J]. Journal of Northeast Forestry University, 2019, 47(11): 66−72.
[18] 符利勇. 非线性混合效应模型及其在林业上应用[D]. 北京: 中国林业科学研究院, 2012. Fu L Y. Nonlinear mixed effects model and its application in forestry[D]. Beijing: Chinese Academy of Forestry, 2012.
[19] 董云飞. 福建杉木人工林单木主长模型的研究[D]. 北京: 北京林业大学, 2015. Dong Y F. Study on individual tree growth model for Chinese fir plantation in Fujian[D]. Beijing: Beijing Forestry University, 2015.
[20] 张海平, 李凤日, 董利虎, 等. 基于气象因子的白桦天然林单木直径生长模型[J]. 应用生态学报, 2017, 28(6): 1851−1859. Zhang H P, Li F R, Dong L H, et al. Individual tree diameter increment model for natural Betula platyphylla forests based on meteorological factors[J]. Chinese Journal of Applied Ecology, 2017, 28(6): 1851−1859.
[21] 杨鑫, 王建军, 杜志, 等. 基于气候因子的兴安落叶松天然林单木直径生长模型构建[J]. 北京林业大学学报, 2022, 44(8): 1−11. Yang X, Wang J J, Du Z, et al. Development of individual-tree diameter increment model for natural Larix gmelinii forests based on climatic factors[J]. Journal of Beijing Forestry University, 2022, 44(8): 1−11.
[22] 王淼, 白淑菊, 陶大立, 等. 大气增温对长白山林木直径生长的影响[J]. 应用生态学报, 1995, 6(2): 128−132. Wang M, Bai S J, Tao D L, et al. Effect of rise in air-temperature on tree ring growth of forest on Changbai Mountain[J]. Chinese Journal of Applied Ecology, 1995, 6(2): 128−132.
[23] 于健, 徐倩倩, 刘文慧, 等. 长白山东坡不同海拔长白落叶松径向生长对气候变化的响应[J]. 植物生态学报, 2016, 40(1): 24−35. doi: 10.17521/cjpe.2015.0216 Yu J, Xu Q Q, Liu W H, et al. Response of radial growth to climate change for Larix olgensis along an altitudinal gradient on the eastern slope of Changbai Mountain, Northeast China[J]. Chinese Journal of Plant Ecology, 2016, 40(1): 24−35. doi: 10.17521/cjpe.2015.0216
[24] 韩艳刚, 周旺明, 齐麟, 等. 长白山树木径向生长对气候因子的响应[J]. 应用生态学报, 2019, 30(5): 1513−1520. Han Y G, Zhou W M, Qi L, et al. Tree radial growth-climate relationship in Changbai Mountain, Northeast China[J]. Chinese Journal of Applied Ecology, 2019, 30(5): 1513−1520.
[25] Wang X, Fang J, Tang Z, et al. Climatic control primary forest structure and DBH-height allometry in northeast China[J]. Forest Ecology and Management, 2006, 234(1): 264−274.
[26] Leites L P, Robinson A P, Rehfeldt G E, et al. Height-growth response to climatic changes differs among populations of Douglas-fir: a novel analysis of historic data[J]. Ecological Applications, 2012, 22(1): 154−165. doi: 10.1890/11-0150.1
[27] Yang Y, Huang S. Effects of competition and climate variables on modelling height to live crown for three boreal tree species in Alberta, Canada[J]. European Journal of Forest Research, 2018, 137(2): 153−167. doi: 10.1007/s10342-017-1095-7
[28] Fortin M, van Couwenberghe R, Perez V, et al. Evidence of climate effects on the height-diameter relationships of tree species [J/OL]. Annals of Forest Science, 2019, 76(1): 1[2022−12−19]. https://link.springer.com/article/10.1007/s13595-018-0784-9.
[29] 国家林业和草原局. 2021中国林草资源及生态状况[M]. 北京: 中国林业出版社, 2022. National Forestry and Grassland Administration. Forest and grassland resources and ecological status in China in 2021 [M]. Beijing: China Forestry Publishing House, 2022.
[30] Wang T, Hamann A, Spittlehouse D L, et al. Climatewna-high-resolution spatial climate data for western north America[J]. Journal of Applied Meteorology and Climatology, 2012, 51(1): 16−29. doi: 10.1175/JAMC-D-11-043.1
[31] Sharma M. Comparing height-diameter relationships of boreal tree species grown in plantations and natural stands[J]. Forest Science, 2016, 62(1): 70−77.
[32] Han Y G, Lei Z Y, Ciceu A, et al. Determining an accurate and cost-effective individual height-diameter model for Mongolian pine on sandy land [J/OL]. Forests, 2021, 12: 1144[2023−01−19]. https://doi.org/10.3390/f12091144.
[33] Vonesh E F, Chinchilli V M. Linear and nonlinear models for the analysis of repeated measurements [M]. New York: Marcel Dekker Inc., 1997.
[34] Neumann M, Mues V, Moreno A, et al. Climate variability drives recent tree mortality in Europe[J]. Global Change Biology, 2017, 23: 4788−4797. doi: 10.1111/gcb.13724
[35] 刘敏, 毛子军, 厉悦, 等. 不同纬度阔叶红松林红松径向生长对气候因子的响应[J]. 应用生态学报, 2016, 27(5): 1341−1352. Liu M, Mao Z J, Li Y, et al. Response of radial growth of Pinus koraiensis in broad-leaved Korean pine forests with different latitudes to climatical factors[J]. Chinese Journal of Applied Ecology, 2016, 27(5): 1341−1352.
[36] Wang M, Zhao Y H, Zhen Z, et al. Individual-tree diameter growth model for Korean pine plantations based on optimized interpolation of meteorological variables[J]. Journal of Forestry Research, 2021, 32(4): 1535−1552.
[37] Franceschini T, Schneider R. Influence of shade tolerance and development stage on the allometry of ten temperate tree species[J]. Physiological Ecology, 2014, 176(3): 739−749.
[38] Tian D Y, Jiang L C, Shahzad M K, et al. Climate-sensitive tree height-diameter models for mixed forests in Northeastern China[J]. Agricultural and Forest Meteorology, 2022, 326: 1−12.
[39] Fang Z, Bailey R L. Nonlinear mixed effects modeling for slash pine dominant height growth following intensive silvicultural treatments[J]. Forest Science, 2001, 47(3): 287−300.
[40] Calama R, Montero G. Interregional nonlinear height-diameter model with random coefficients for stone pine in Spane[J]. Canadian Journal of Forest Research, 2004, 34: 150−163. doi: 10.1139/x03-199
[41] 符利勇, 唐守正, 张会儒, 等. 基于多水平非线性混合效应蒙古栎林单木断面积模型[J]. 林业科学研究, 2015, 28(1): 23−31. Fu L Y, Tang S Z, Zhang H R, et al. Multilevel nonlinear mixed-effects basal area models for individual trees of Quercus mongolica[J]. Forest Research, 2015, 28(1): 23−31.