-
单木模型通常以树木个体的生长信息为基础，从林木生长的竞争机制出发，模拟林分内不同大小树木的生长过程[1-2]。与林分和径阶水平模型相比，单木模型不仅能够提供较为详细的林分结构特征及动态变化信息，而且还是研究林木生长量、收获量以及生物量的重要基础，因此单木生长模型是定量研究林木自身生长过程及其对各种干扰过程响应机制的有效手段。根据模型中是否包含距离信息，单木生长可分为与距离有关和与距离无关模型[3]，其中由于前者需要详细的距离信息，这极大程度上限制该类模型在林业生产中的应用，因此与距离无关的单木生长模型一直是国内外林业研究的重要内容。
传统研究认为林木直径生长主要受树木自身大小、竞争因子和立地条件的影响[1, 4-5]。但近些年来，国内外学者逐渐关注到物种多样性和林分空间结构特征对树木生长的影响。在物种多样性方面，Vannoppen等[6]发现林分物种多样性对北欧地区山毛榉（Fagus sylvatica Linn.）和夏栎（Quercus robur Linn.）直径生长量的贡献量分别为2%和21%；Wang等[7]研究表明中国东北地区近熟林和过熟林内分别有66.7%和64.3%的样地内生物量和木本植物多样性间存在显著相关性；Liang等[8]基于全球44个国家约77万个固定样地数据的研究结果表明：生物多样性与生产力间存在显著的正相关；谭凌照等[9]发现：在蛟河阔叶红松（Pinus koraiensis Sieb. et Zucc.）林中物种多样性均匀度指数与生产力显著负相关。在林分空间结构特征方面，汤孟平等[10]研究表明：较高聚集程度、年龄隔离度以及较低的最邻近竹株数和竞争指数均有利于增加毛竹（Phyllostachys edulis (Carr.) H.de Lehaie）林单位面积生物量；吕延杰等[11]研究表明：当红皮云杉（Picea koraiensis Nakai）和冷杉（Abies fabri (Mast.) Craib）处于随机分布时胸径生长量最大，且其随着混交度的增加而增加。总的来看，物种多样性与林分空间结构特征对单木直径生长的作用已得到广泛认同，但现阶段还鲜有将其整合到单木生长模型中的研究。
现阶段单木生长模型数据源主要包括2类，即固定样地复测数据[1-2, 4-5]和单木解析数据[12]。从统计学角度来看，2种数据均存在明显的嵌套效应（即区域→样地→单木）、空间自相关性以及时间自相关性等问题，这些均严重违背传统统计学模型中对数据正态性、独立性和方差齐性的要求[13]，致使所建模型不能很好满足林业生产需要。与传统回归模型不同，混合模型只要求误差项为正态分布，而对模型误差项之间是否相互独立和是否等方差不做要求，所以混合模型能够处理更多类型的数据。目前，国内外学者采用混合模型对单木生长模型也进行了一些探索，如姜立春等[12]采用解析木数据建立了兴安落叶松（Larix gmelinii (Rupr.)Kuzen）单木直径和树高模型，但并未考虑林分变量、物种多样性和结构参数对生长的影响；彭娓等[4]基于固定样地复测数据建立了长白落叶松（Larix olgensis Henry）人工林的单木直径生长模型；王蒙等[5]则进一步考虑了抚育间伐对落叶松人工林单木直径生长的影响；祖笑锋等[14]采用解析木数据建立了美国黄松（Pinus ponderosa Dougl. ex Laws.）林分优势高生长模型。这些研究均表明，混合模型不仅能够大幅度提升模型的拟合效果，同时还能有效反应区域环境、林分、经营等变量对林木生长过程的潜在影响。
由于长期以来过度的人为采伐和严重的火灾干扰，当地顶级群落兴安落叶松林已逐渐退化为以兴安落叶松和白桦（Betula platyphylla Suk.）为主要优势树种的天然次生林，其结构、功能和生产力均与原始顶级群落存在显著差异。因此，通过建立该地区次生林中兴安落叶松单木直径生长模型，对揭示林木生长与林分因子、单木因子、物种多样性以及空间结构特征间的关系具有重要作用。为此，本研究以兴安落叶松天然次生林为研究对象，利用混合模型技术在传统单木生长模型基础上逐渐引入林分因子、单木因子、物种多样性以及林分空间结构参数，进而建立该地区兴安落叶松单木直径生长模型，为该地区次生林的有效恢复提供理论依据和技术支撑。
-
盘古林场（52°16′38″～52°47′04″ N、123°20′02″～124°21′40″ E）属于大兴安岭地区塔河林业局管辖，总经营面积15.2万hm2。该区属于大兴安岭中部伊勒呼里山北麓，年均气温−3℃，极端最低温和最高温分别为−53℃和36℃，平均生长季约100 d；地带性土壤为棕色针叶林土，其中土壤有机质、全氮和全磷平均含量分别为330.48 g·kg−1、7.01 mg·kg−1和81.24 g·kg−1[15]；地带性植被为兴安落叶松，但现阶段主要为以兴安落叶松和白桦为主要树种的天然次生林，其伴生树种主要包括樟子松（Pinus sylvestris var. mongolica Litv.）、红皮云杉、山杨（Populus davidiana Dode）等。
-
在全面踏查基础上，于2011年7—8月份在盘古林场选择不同年龄、不同密度、不同立地条件的林分设置固定样地18块。各样地面积20 m × 30 m，采用相邻格网调查法，实测样地内直径 ≥ 5 cm树木的树种、胸径、树高、枝下高、冠幅、坐标以及状态等指标。每块样地按等断面积径级标准木法，将林木大小划分为5个等级，将各等级林木的平均胸径和树高作为选择解析木的标准，并进一步参考该样地树种组成在样地附近选伐样木。将解析木树干按1 m区分段进行区分，分别在胸高处（1.3 m）和各区分处截取3 cm厚圆盘，并在非工作面上标出南、北方向。在室内将圆盘的工作面刨光，并通过髓心分别画出东西、南北2条直线，进而采用200 dpi分辨率扫描成jpg图像。将图像导入WinDENDRO软件中，根据图像中早晚材色泽差异识别出各年轮具体位置，且要求各方向中识别出的年轮数必须严格相等，以各方向各年轮宽度的平均值作为该年龄下的平均直径。从18块样地中共收集兴安落叶松解析样木50株（其余40株均为其他树种）。林分和单木统计数据如表1所示。
表 1 大兴安岭天然落叶松林分和单木调查因子统计
Table 1. Descriptive statistics of stand- and tree-level of natural Lari xgmelinii stand in Daxing′anling Mountains
组别
Group变量
Variable最小值
Min.最大值
Max.平均值
Mean标准差
Std.变异系数
C.V./%林分 Stand (n = 18) 平均胸径 Mean DBH/cm 7.00 14.80 10.96 2.04 18.64 平均树高 Mean HT/m 6.66 19.55 12.19 2.73 22.44 株数 Number/(trees·hm−2) 567.00 3850.00 1887.69 777.63 41.19 蓄积 Volume/(m3·hm−2) 144.66 646.67 366.71 136.06 37.10 郁闭度 Crown density 0.40 0.85 0.67 0.08 11.94 Margalef指数 Ma 0.19 0.87 0.54 0.20 36.82 Simpson指数 Sp 0.05 0.73 0.50 0.19 37.66 Shannon-Wiener指数 Sw 0.21 0.57 0.37 0.11 28.46 Pielow指数 Pw 0.25 1.41 1.02 0.30 29.14 混交度 M 0.03 0.48 0.25 0.11 45.36 大小比 U 0.42 0.53 0.50 0.02 3.73 角尺度 W 0.49 0.59 0.54 0.03 5.20 竞争指数 H 0.24 0.34 0.29 0.02 7.48 单木 Tree (n = 50) 年龄 Age/a 20.00 234.00 70.60 42.88 60.74 胸径 DBH/cm 6.50 38.80 16.44 7.34 44.66 树高 HT/m 6.10 26.00 14.64 4.16 28.42 冠幅 CW/m 0.70 2.88 1.69 0.51 30.08 冠长 CL/m 2.60 21.00 6.94 3.03 43.66 Note: DBH: diameter at breast height; HT: the height of tree; CW: crown width; CL: crown length; Ma, Sp, Sw and Pw are the diversity indexes of Margalef, Simpson, Shannon-Wiener and Pielow, respectively;U, M, W and H are respectively the indexes of neighborhood comparison, species mingling, uniform angle and competition. -
由于树木直径生长通常具有典型的非线性特征，因此选用一些常用的非线性曲线作为备选模型，包括Mitscherlich方程、Logistic方程、Weibull方程、Gompertz方程、Schumacher方程、Korf方程、Hosfled方程以及Richards方程，各方程形式详见孟宪宇[3]，此处不再赘述。根据模型参数估计值、拟合效果以及残差分布图等特征，确定适于描述天然落叶松单木直径生长过程的基础模型。考虑到树木所处环境会对林木直径生长产生影响，因此除将常规林分因子和单木因子加入到模型中外，本研究还尝试引入林分物种多样性、林分空间结构参数指标，对基础模型的参数进行再参数化处理。林分因子主要包括林分平均胸径（Dg）、平均树高（Hg）、株数密度（N）、郁闭度（Dc）等，单木因子包括胸径（DBH）、树高（HT）、冠幅（CW）、冠长（CL）等，物种多样性指标包括Margalef丰富度指数（Ma）、Simpson多样性指数（Sp）、Shannon-Wiener指数（Sw）、Pielow均匀度指数（Pw），林分空间结构指标包括混交度（M）、大小比（U）、角尺度（W）和Heygi竞争指数（H）。此外，还进一步考虑上述各变量间的组合及其变换形式（如对数、平方等）。各物种多样性指数和结构参数具体定义详见谭凌照等[9]和吕延杰等[11]，此处不再赘述。
-
采用混合模型建立兴安落叶松单木直径生长模型，一是为通过引入随机因素而提高模拟精度，二是为了有效处理解析木数据中存在的异方差和时间自相关问题。根据Pinheiro等定义[13]，非线性混合效应模型的基本形式为[13]：
$ \left\{ {\begin{array}{*{20}{l}} {{y_{ij}} = f\left({{\varphi _{ij}},{v_{ij}}} \right) + {e_{ij}}}\\ {{\varphi _{ij}} = {A_{ij}}\beta + {B_{ij}}{b_i}} \end{array}} \right.$
(1) 式中：
$ {y}_{ij} $ 为第i棵树的第j次观测值；f()为含有参数向量$ {\varphi }_{ij} $ 和协变量$ {v}_{ij} $ 的函数；$ {A}_{ij} $ 为已知固定效应的设计矩阵；$ \beta $ 为未知的固定参数向量；$ {B}_{ij} $ 为已知随机效应的设计矩阵；$ {b}_{i} $ 为带有方差协方差结构矩阵的未知随机效应参数向量，其服从期望值为0、方差为D的正态分布；D为随机效应的协方差矩阵；$ {e}_{ij} $ 为服从正态分布（期望为0、方差为$ {\sigma }^{2}{R}_{i} $ ）的误差项；Ri为第i株树木的方差协方差矩阵。因各样地解析木数量相对较少，因此本研究仅考虑单木水平的随机因子。对已确定的基础模型，构建混合模型时还需确定随机参数、方差协方差结构、时间序列结构以及异方差校正函数共4类参数[12-13]。对随机参数，通常认为所有参数都有成为随机效应参数的可能，因此可对不同随机参数的组合模型进行拟合，通过多项统计指标来评价模型拟合效果。对混合模型中的方差协方差结构，本研究检验林业上3种常用结构，即复合对称矩阵CS、对角矩阵UN(1) 和广义正定矩阵UN；为进一步提高模型拟合效果，采用复合对称结构CS、一阶自回归结构AR(1)、一阶自回归与滑动平均结构ARMA(1,1)以及高斯结构corGaussian来描述解析木数据中的时间序列相关性。采用指数函数varExp、幂函数varPower和常数加幂函数varConstPower来校正解析木数据中的异方差现象。各模型拟合优度通过赤池信息准则（AIC）、贝叶斯准则（BIC）和负对数似然值（−LL）来确定。此外，对不同参数个数的模型，还需进行似然比检验（LRT）以最终确定最优模型。
-
混合模型中包含固定参数β和随机参数bi两部分，其中固定效应参数与普通模型相同，而随机参数则需要根据模型拟合结果和检验数据重新计算，其具体公式为[13]：
$ {\widehat{b}}_{k}={\widehat{D}{\widehat{\mathrm{Z}}}_{k}({\widehat{\mathrm{Z}}}_{k}\widehat{D}{{\widehat{\mathrm{Z}}}_{k}}'+{\widehat{\mathrm{R}}}_{k})}^{-1}{\widehat{\mathrm{e}}}_{k} $
(2) 式中：
$ \widehat{D} $ 为随机参数的方差协方差结构；$ {\widehat{Z}}_{k} $ 为第k棵树的设计矩阵；$ {\widehat{R}}_{k} $ 为误差项的方差协方差矩阵；$ {\widehat{e}}_{k} $ 为基础模型的误差项。各基础模型和混合效应模型的精度评价和检验均采用五折交叉检验法进行，即建模过程需采用全部样本数据，而检验过程则将数据随机分为5份，每份单独进行评价，以其平均值作为最终结果。各模型采用调整系数（Ra2）、平均绝对误差（Bias）和均方根误差（RMSE）3项指标进行评价和精度检验[4-5]，其中R2越大越好，Bias和RMSE均越接近0越好。
-
分别采用S-Plus软件中的NLS模块拟合各基础模型，采用NLME模块拟合各混合效应模型；其余数据统计和图表绘制均采用EXCEL完成。
-
在8种候选模型中，仅Mitscherlich、Schumacher、Logistic和Gompertz方程能够收敛，且各模型参数均能通过α = 0.05水平的显著性检验（表2）。虽然各项指标均显示Gompertz方程的拟合效果最好，但方差分析结果表明Mitscherlich方程和Gompertz方程间的差异不显著（P = 0.09），且在后续再参数化过程中也较难收敛，因此在综合评估各模型拟合优度、模型参数个数以及残差分布图后，选择Mitscherlich方程为大兴安岭地区兴安落叶松单木直径生长的最优基础模型。
表 2 基础模型参数估计及拟合优度
Table 2. The parameter estimation and goodness-of-fit for the base models
编号
NO.模型
Model公式
Formulation参数估计值
Parameter estimation拟合优度
Goodness-of-fita b c 调整系数
Ra2平均绝对误差
Bias/cm均方根误差
RMSE/cmB1 Mitscherlich $y=a(1-\mathrm{E}\mathrm{x}\mathrm{p}(-b{t}\left)\right)$ 23.41 0.01314 0.6285 2.47 3.17 B2 Schumacher $y=a\mathrm{E}\mathrm{x}\mathrm{p}(-b/{t})$ 18.70 23.93 0.5897 2.77 3.38 B3 Logistic $y=a/(1+b\mathrm{E}\mathrm{x}\mathrm{p}(-c{t}\left)\right)$ 18.75 5.755 0.04191 0.6243 2.42 3.17 B4 Gompertz $y=aE\mathrm{x}\mathrm{p}(-b\mathrm{E}\mathrm{x}\mathrm{p}(-c{t}\left)\right)$ 21.92 2.232 0.02345 0.6302 2.39 3.14 注：y表示胸径，t表示年龄，a、b和c为待估参数;
Note:y is DBH, t is age, and a、b、c are the estimated parameters.对不同单木进行Mitscherlich式拟合，结果表明模型收敛率达到82.86%，可以有效反映所有建模样本的整体特征。Mitscherlich模型的参数a和b与林分变量、单木变量、物种多样性、结构参数的相关性分析（图1）。可以看出，参数a与林分变量中的平均树高（r = 0.5289）、单木变量中的胸径（r = 0.7871）、树高（r = 0.7962）和冠幅（r = 0.5264）均显著正相关，而参数b则仅与林分变量中的平均树高（r = −0.5042）显著负相关。为进一步揭示物种多样性和林分结构参数对直径生长的影响，结合相关性分析、变量共线性以及参数可获得性等问题，综合确定参数a的候选变量包括林分平均树高、单木胸径以及混交度和角尺度，而参数b的候选变量包括林分平均胸径、Pielow均匀度指数（与Sw指数的r = 0.7704）以及Heygi竞争指数，其中平均树高和单木胸径为参数a的固定选项，平均胸径为参数b的固定选项。
图 1 Mitscherlich模型参数a和b与林分变量、单木变量、物种多样性以及林分结构参数的相关性
Figure 1. Relations between the parameters of Mitscherlich and the indices of stand-, tree-level, diversity and spatial structures
因候选变量个数较多，因此本研究采用相乘思想构造函数对Mitscherlich模型中参数a和b进行预测，具体形式如下：
$ a={a}_{0}{X}_{1}^{{a}_{1}}{X}_{2}^{{a}_{2}}\cdots {X}_{N}^{{a}_{N}} $
(3) $ b={b}_{0}{X}_{1}^{{b}_{1}}{X}_{2}^{{b}_{2}}\cdots {X}_{N}^{{b}_{N}} $
(4) 式中：a0, a1, ···,aN为Mitscherlich模型参数a的待估系数，b0,b1,···, bN为Mitscherlich模型参数b的待估系数，N为候选变量个数，
$ {X}_{1}{，X}_{2}，\cdots {，X}_{N} $ 分别为各候选变量。当加入林分和单木变量后，模型M0的拟合优度较基础模型B1显著增加约18%（以Ra2值为例，下同；表3）。对混交度、角尺度、Pielow均匀度以及竞争指数4个变量而言，当仅考虑1个变量时，Heygi竞争指数对模型的提升效果最为明显，其精度较模型M0显著增加约7.48%；当考虑2个指标时，混交度和Heygi竞争指数的组合使模型精度提升约10.95%；当考虑3个指标时，混交度&角尺度&Heygi组合的效果最佳（11.21%）；当考虑4个指标时，模型拟合效果较M0增加约11.22%。方差分析结果表明：模型M0与B1、M4与M1、M6与M4、M14与M6均存在极显著差异（P < 0.0001），而模型M14与M15间的差异则不显著（P = 0.4157）。因此，在综合考虑各模型参数估计值、残差分布以及拟合优度等指标基础上，综合确定M14为兴安落叶松单木直径生长的最优广义模型，其具体形式为：
表 3 基于林分、单木、物种多样性和结构参数的单木直径生长模型拟合优度
Table 3. Goodness-of-fit for diameter growth model with stand-level, tree-level, diversity and structure variables
模型
Model参数a
Parameter a参数b
Parameter b拟合优度
Goodness-of-fit方差分析
ANOVAHg DBH M W Dg Pw H 调整系数
Ra2平均绝对误差
Bias/cm均方根误差
RMSE/cmB1 0.6285 2.47 3.17 M0 √ √ √ 0.7431 2.11 2.69 < 0.0001 M1 √ √ √ √ 0.7711 1.96 2.52 M2 √ √ √ √ 0.7453 2.11 2.68 M3 √ √ √ √ 0.7456 2.04 2.66 M4 √ √ √ √ 0.7987 1.92 2.43 < 0.0001 M5 √ √ √ √ √ 0.7984 1.84 2.39 M6 √ √ √ √ √ 0.8245 1.67 2.23 < 0.0001 M7 √ √ √ √ √ 0.7603 2.01 2.58 M8 √ √ √ √ √ 0.802 0 1.88 2.41 M9 √ √ √ √ √ 0.7912 1.83 2.40 M10 √ √ √ √ √ 0.8163 1.73 2.27 M11 √ √ √ √ √ √ 0.8257 1.68 2.22 M12 √ √ √ √ √ √ 0.8186 1.72 2.28 M13 √ √ √ √ √ √ 0.7984 1.82 2.38 M14 √ √ √ √ √ √ 0.8264 1.68 2.21 < 0.0001 M15 √ √ √ √ √ √ √ 0.8265 1.68 2.21 0.4157 注：√表示包含该参数；
Note: √ indicates that the variable was included.$ y\left(t\right)={a}_{0}{\cdotHg}^{{a}_{1}}\cdot{DBH}^{{a}_{2}}\cdot {M}^{{a}_{3}}\cdot{W}^{{a}_{4}}\cdot(1-\mathrm{E}\mathrm{x}\mathrm{p}(-{b}_{0}{\cdotDg}^{{b}_{1}}{\cdotH}^{{b}_{2}}\cdot{t}\left)\right) $
(5) 式中：y(t)为第t年时树木的直径；DBH为解析木的胸径（伐倒时）；其余变量如前所述。
-
在式5基础上，引入随机参数、方差协方差结构、异方差校正函数和时间序列结构等构造兴安落叶松单木直径生长的混合效应模型（表4）。对随机参数而言，仅有8种参数组合能够收敛；当考虑1个随机参数时，模型E2的AIC值最小（5628.26）；当考虑2个随机参数时，模型E8的AIC值最小（5630.02），但较E2略有增加（△AIC = 1.76），且似然比检验结果表明模型E8与E2间的差异也不显著（P = 0.6276）。因此，将a1作为兴安落叶松单木直径生长混合效应模型的最优随机参数组合。
表 4 兴安落叶松单木直径生长混合效应模型筛选结果
Table 4. Selection results of diameter growth model of natural Larix gmelinii stand when using mixed-effects regression techniques
项目
Items模型
Model参数
Parameter参数个数
Parameter number赤池信息准则
AIC贝叶斯信息准则
BIC对数似然值
Log-likelihood似然比检验
LRTP值
P−value随机参数
Random parameterE1 None 9 8773.91 8824.23 −4377.96 E2 a1 10 5628.26 5684.17 −2804.79 3158.93 < 0.0001 E3 a2 10 5629.59 5685.50 −2804.79 E4 a3 10 5644.69 5700.61 −2812.35 E5 a4 10 5631.07 5686.99 −2805.54 E6 b1 10 5321.24 5377.15 −2650.62 E7 a1,a2 11 5630.26 5691.76 −2804.13 E8 a1,a4 11 5630.02 5691.53 −2804.01 0.2353 0.6276 E9 a3,a4 11 5632.58 5694.08 −2805.29 方差函数
Variance functionE10 幂函数 Power 11 5611.80 5673.31 −2794.90 18.4562 < 0.0001 E11 指数函数 Exp 11 5625.34 5686.85 −2801.67 E12 常数加幂函数
ConstPower12 5583.78 5650.88 −2778.89 30.0214 < 0.0001 时间序列结构
Autocorrelation structureE13 一阶自回归与滑动平均
ARMA(1,1)12 3108.50 3175.59 −1542.25 2505.31 < 0.0001 E14 高斯结构
corGaussian13 −426.16 −353.47 226.08 3536.65 < 0.0001 鉴于模型E2中仅包含1个随机参数，本研究不再进行混合模型方差协方差结构的选择。3种异方差校正函数检测结果表明，常数加幂函数（模型E12）对异方差的校正效果显著优于幂函数（模型E10；P < 0.01）和指数函数（模型E11；P < 0.01）。但以常数加幂函数（模型E12）为基础构造的包含时间序列结构的混合效应模型均不能收敛，因此本研究选择次优的幂函数（即E10）为兴安落叶松单木直径生长混合效应模型的最适异方差校正函数。在模型E10基础上，进一步测试4种时间序列结构对模型拟合效果的影响，结果表明仅一阶自回归与滑动平均结构ARMA(1,1)和高斯结构corGaussian能够收敛，且corGaussian结构显著优于ARMA(1,1)结构（P < 0.01）；模型E13与E14间的差异同样达到极显著水平（P < 0.01）。因此，将模型E14作为兴安落叶松单木直径生长的最优混合效应模型。
兴安落叶松单木直径生长基础模型、广义模型和混合模型的参数估计以及拟合结果如表5所示。除混合模型参数a0和b2
外（P = 0.076 4和0.0953），其余各参数均达到显著水平（P < 0.05）。3个模型统计结果表明：混合模型的Ra2分别较基础模型和广义模型增加约54%和17%，而平均绝对误差和均方根误差分别下降约72%、70%和59%、57%；混合模型不仅显著缩小了残差的范围，而且异方差现象也得到了显著削弱（图2），表明模型E14能够满足兴安落叶松单木直径预测的需求。 表 5 兴安落叶松单木直径生长基础模型、广义模型和混合模型的参数估计和拟合结果
Table 5. Parameter estimation and fitting results of base-, general- and mixed-model for natural Larix gmelinii stand
项目
Items参数
Parameter基础模型
Base model (B1)广义模型
General model (M14)混合模型
Mixed model (E14)固定参数
Fixed parametera or a0 23.41 (< 0.0001) 178.95 (< 0.0001) 11.70 (0.0478) a1 −0.001994 (< 0.0001) −0.5389 (0.0284) a2 −0.5328 (< 0.0001) 1.466 (0.0004) a3 0.7033 (< 0.0001) 0.4928 (0.0764) a4 0.1614 (< 0.0001) 0.3669 (0.0005) b or b0 0.01314 (< 0.0001) 0.005387 (< 0.0001) 0.05265 (0.0473) b1 1.1188 (< 0.0001) −1.068 (0.0003) b2 2.1024 (< 0.0001) −0.05710 (0.0953) 模型误差
Error variance$ {\sigma }^{2} $ 10.06 4.865 0.6647 随机参数
Random variance$ {\sigma }_{{\mathrm{a}}_{1}}^{2} $ 0.1854 方差函数
Variance function幂函数 Power −0.01106 时间序列结构
Autocorrelation structure范围 Range 2.143 块金值 Nugget 0.1000 拟合统计量
Goodness-of-fit确定系数 Ra2 0.6285 0.8264 0.9710 平均绝对误差 Bias/cm 2.47 1.68 0.69 均方根误差 RMSE/cm 3.17 2.21 0.94 
图 2 兴安落叶松单木直径生长基础模型、广义模型和混合模型的残差分布
Figure 2. Residual distribution of base-, general- and mixed-model for natural Larix gmelinii stand in Daxing′anling Mountains
-
五折交叉检验结果表明，兴安落叶松单木直径生长混合模型的精度同样显著高于基础模型和广义模型（表6），能够满足林业生产的需要。结合表1，本研究假设林分平均胸径Dg = 10 cm, 林分平均高Hg = 12 m, 单木胸径DBH = 15 cm，角尺度W = 0.5，混交度M = 0.5以及林分平均竞争指数H = 0.3为一组基础的模拟情景，即在每次模拟中仅改变1个变量的值而保持其他不变。以混合效应模型E14的固定效应参数为基础，模拟结果表明单木直径生长量随林分平均胸径、平均树高以及竞争指数的增大而减小，但随着单木胸径、角尺度和混交度的增加而增加（图3）。整体来看，各变量对单木直径大小的贡献量（定义为各变量每增加1个单位时对应的胸径增长或减少量）依次表现为：混交度（+4.71 cm） > 角尺度（+3.62 cm） > 竞争指数（−0.73 cm） > 单木胸径（+0.45 cm） > 林分平均胸径（−0.34 cm） > 林分平均树高（−0.21 cm）。
表 6 兴安落叶松单木直径生长基础模型、广义模型和混合模型的独立样本检验
Table 6. Validation of base-, general- and mixed-model for natural Larix gmelinii stand in Daxing'anling Mountains
模型
Model调整系数
Ra2平均绝对误差
Bias/cm均方根误差
RMSE/cm基础模型
Base model [B1]0.6273 2.52 3.25 广义模型
General model [M14]0.8577 1.92 2.14 混合模型
Mixed model [E14]0.9838 0.56 0.72 
图 3 不同变量对兴安落叶松林单木直径生长的影响
Figure 3. Effects of different variables on the diameter growth of natural Larix gmelinii Stand in Daxing'anling Mountains
-
在天然异龄林中，光照是限制树木生长的重要因素之一。显然，大树较小树具有显著的竞争优势，因此初始胸径的大小可充分代表树木在林分中的竞争态势，即一定范围内树木胸径越大，其直径生长量也越大。根据模型E14测算结果，树木初始胸径每增加1 cm时，单木直径生长量可平均增加4.50 cm左右（生长至150 a 时，下同）；测算结果表明，林分平均胸径每增加1 cm，单木直径生长量将会下降约3.40 cm，而林分平均树高每增加1 m，单木直径生长量则会下降2.11 cm。综上所述，可通过适当的采伐来降低林分整体的竞争水平，这对促进林木直径生长具有重要作用。
除了上述基本林分和单木因子外，本研究还进一步测试了4种物种多样性指数和4种林分空间结构参数对模型拟合精度的影响（图1）。对天然林而言，林分内不同树种的搭配是经过长期的自然选择形成的，因此各树种间存在较强的优势互补关系，这是进行物种多样性与林分生产力关系研究的理论基础。结果表明，Mitscherlich模型参数a与所有物种多样性指数均不相关，而参数b则与Sw呈负相关关系、与Pw呈正相关关系，这与谭凌照等[9]研究结果相似。鉴于Sw与Pw间显著相关（r = 0.77），因此本研究仅尝试将Pw引入到模型中，但结果表明模型（M15）拟合效果并不理想，其与模型M14间的差异不显著（P = 0.4157）。因此，本研究所建模型中未包含任何生物多样性指标，但这可能与研究区域物种多样性水平整体偏低有关。在物种多样性丰富的地区，该结论是否成立仍有待于进一步检验。此外，本研究仅考虑了林分物种多样性指标，而与林木生长密切相关的结构多样性并未充分考虑[16]，这也有待于进步研究。
林分空间结构参数是近年来国内外森林经营领域研究的热点。在此基础上，惠刚盈及其团队提出的结构化森林经营理论被认为是实现森林可持续经营的重要途径之一[17]，但现阶段还很少有将其整合到传统林木生长模型的研究中。汤孟平等[10]、吕延杰等[11]虽然揭示了林分空间结构参数与林木生长的关系，但并未给出实用可行的预测模型。本研究在他们的基础上，进一步研究发现Mitscherlich模型参数a与混交度（r = −0.23）、大小比（r = −0.19）和角尺度（r = −0.26）间均存在较强的相关性，而参数b则仅与Heygi竞争指数相关（r = −0.12）。经过综合评估，本研究试图将混交度和角尺度引入到参数a中，而将Heygi竞争指数引入到参数b中，结果表明各模型Ra2值分别+3.77%、+0.30%和+7.48%。因此，最终将混交度、角尺度和Heygi竞争指数引入到模型中，模型模拟结果表明林木混交程度越高、水平分布越随机，其直径生长量越大；林木竞争越激烈，其直径生长量越小，这些均与前人研究结果一致[10-11]。模型E14模拟结果表明，混交度、角尺度和Heygi竞争指数的数值每增加10%，其对应的直径生长量平均分别为+4.71 cm、+3.62 cm和−0.73 cm，但因不同的林分特征和生长阶段而异（图3）。
与传统统计学模型相比，混合模型能够通过固定参数反应总体的变化趋势，而通过随机参数来反映树木个体生长的差异。此外，混合模型还能提供其他一些丰富的信息，如方差协方差结构能够反应不同林分条件下树木直径生长过程的差异，这对复杂的异龄混交林而言具有明显优势；模拟结果表明，兴安落叶松单木直径生长的混合模型仅包含1个随机参数，其方差
$ {\sigma }_{{\mathrm{a}}_{1}}^{2} =0.1\;854$ 。异方差函数则能有效解决解析木数据中的时间自相关性，进而确保得到的参数估计值是无偏的；模拟结果表明兴安落叶松单木直径生长混合模型的最适异方差校正函数为幂函数。模型模拟和检验结果表明，相较于广义模型（M14），混合模型E14的各项指标和残差分布图均得到显著改善，这与其他学者的结果一致[12-13]。但需要强调的是，混合模型的检验和应用过程较传统模型要复杂的多，必须借助少量样本、模型参数估计结果（表5）和式（2）实现不同林分条件下树木个体直径生长过程中随机参数的确定，因此开发具有针对性的应用软件是当务之急。此外，由于样本量不足的问题，本研究仅考虑了单木水平的随机因子，而没有同时考虑区域、样地和单木等多水平的随机效应，随着数据的积累，可在这方面继续深入研究。 -
林分结构显著影响兴安落叶松单木的直径生长过程。当引入林分角尺度和竞争指数后，广义基础模型M0的Ra2值较基础模型B1增加约18%；当进一步引入混交度、角尺度和竞争指数3个结构参数后，广义最优模型M14的（Ra2）值进一步增加约11%；而混合模型（E14）则使模型精度进一步提升约17%，且模型异方差得到显著削弱。模拟结果表明：兴安落叶松单木直径生长随林分平均胸径、平均树高和平均竞争指数的增加而减小，但随单木胸径、林分角尺度（接近随机分布，即W = 0.5）和林分混交度的增加而增加。各变量对单木直径大小的贡献量（各指标每增加10%）依次表现为：混交度（+4.71 cm） > 角尺度（+3.62 cm） > 竞争指数（−0.73 cm） > 单木胸径（+0.45 cm） > 林分平均胸径（−0.34 cm） > 林分平均树高（−0.21 cm）。因此，在后续森林经营中，应通过采伐逐渐降低林分的竞争水平，注重引入阔叶树种增加林木的混交程度，并逐步调整林木的水平格局为随机分布模式。
基于空间结构参数的大兴安岭天然落叶松单木直径生长模型
Diameter Growth Model Using Spatial Structure Parameters of Natural Larix gmelinii Stand in Daxing'anling Mountains, Northeast China
-
摘要:
目的 研究林分空间结构对兴安落叶松单木直径生长的影响，为大兴安岭天然落叶松次生林的有效恢复提供理论依据和技术支撑。 方法 以大兴安岭天然落叶松次生林为研究对象，基于塔河林业局盘古林场50株兴安落叶松解析木数据，采用S-Plus软件中的NLS函数对8种常用的林木生长模型进行拟合优度评价，进而确定最优基础模型；通过再参数化方法，考虑常用林分因子、单木因子、物种多样性以及林分空间结构等参数对林木直径生长的影响；最终，采用NLME函数构建兴安落叶松单木直径生长的最优混合效应模型，并系统评价不同模型的拟合优度。 结果 兴安落叶松单木直径生长的最优基础模型为Mitscherlich方程（调整系数Ra2 = 0.6285），其参数a分别与林分平均树高、单木胸径、混交度和角尺度显著相关，而参数b则与林分平均胸径和竞争指数显著相关；据此，建立的兴安落叶松广义直径生长模型的Ra2值相较基础模型提升约31%（Ra2 = 0.8264）；通过引入随机参数（参数a1）、异方差函数（幂函数varPower）、时间序列函数（高斯结构corGaussian），单木水平的混合效应模型进一步将Ra2值提升约17%（Ra2 = 0.9710）；五折交叉检验结果表明，混合模型的调整系数、平均绝对误差和均方根误差分别为0.9838、0.56 cm和0.72 cm，能够满足单木直径生长预测的需要。模型模拟结果表明，单木直径生长量随着单木胸径、角尺度和混交度的增加而增加，但随着林分平均胸径、平均树高和竞争指数的增加而减小；当各变量数值每增加10%时，其对应的胸径生长量依次为：+0.45 cm、+3.62 cm、+4.71 cm、−0.34 cm、−0.21 cm和−0.73 cm，但取决于不同的林分特征和生长阶段。 结论 林分角尺度、混交度和竞争水平对兴安落叶松单木直径生长具有显著作用，在后续森林经营中应通过合理的采伐来降低林木的竞争水平，并注重创建和维护林木的随机分布格局和混交程度。 Abstract:Objective To study the effects of stand spatial structure on individual diameter growth of natural Larix gmelinii stand so as to provide reference and support for the effective recovery of secondary forests in the region. Method Based on 50 sample trees of L. gmelinii collected from the secondary forests in Pangu Forest Farm of Tahe Forestry Bureau, eight commonly used forest growth models were compared and evaluated using NLS function in S-Plus software, and the best basic model was determined. Then, the effects of commonly used stand factors, individual tree factors, species diversity, and stand spatial structure on the diameter growth were quantified using the parameterized method. Finally, the optimal mixed-effect model of growth for L. gmelinii trees were developed using NLME function, and the goodness-of-fit of different models were systematically evaluated. Result Mitscherlich equation was proved to be the optimal basic model of diameter growth for natural L. gmelinii stand, in which the coefficient of determination (Ra2) of the function was 0.6285. The parameter a highly related with the DBH of individual trees, and the mean tree height, the mingling and the uniform of the stand, while the parameter b mainly related with the mean DBH and the competition index of stand. The generalized diameter growth model could improve the value of Ra2 by about 31% (Ra2 = 0.8264) when comparing with the basic model. However, the technique of mixed-effect model could further improve the value of Ra2 by about 17% (Ra2 = 0.9710) when the random parameters (a1), heteroscedasticity function (varPower), and autocorrelation structure (corGaussian) were introduced. The results of 5-fold cross-validation showed that the coefficient of determination, average absolute error and root mean square of the estimated mixed-effect model were 0.9838, 0.56 cm and 0.72 cm respectively, indicating the outstanding prediction accuracy of the model. The simulation results showed that the amount of diameter growth increased significantly with the increase of the DBH of individual tree, the mingling and the uniform values (the closer to 0.5, the better), but decreased dramatically with the increases of the average DBH, the average tree height, and the Heygi's competition index. The corresponding amount of diameter growth would be: +0.45 cm, +3.62 cm, +4.71 cm, −0.34 cm, −0.21 cm and −0.73 cm, respectively, if the value of each variable increased by ten percent, however it highly depended on the stand characteristics and growth stages. Conclusion The levels of stand uniform and competition will affect the diameter growth of natural Larix gmelinii stand significantly, thus decreasing the competition levels of forest and creating (or maintaining) the random distributions and higher mingling degrees of the trees by appropriately cutting are highly recommended for the forest management in the near future. -
Key words:
- Larix gmelinii
- / spatial structure
- / diameter growth
- / biodiversity
- / mixed model
-
图 1 Mitscherlich模型参数a和b与林分变量、单木变量、物种多样性以及林分结构参数的相关性
Figure 1. Relations between the parameters of Mitscherlich and the indices of stand-, tree-level, diversity and spatial structures
图 2 兴安落叶松单木直径生长基础模型、广义模型和混合模型的残差分布
Figure 2. Residual distribution of base-, general- and mixed-model for natural Larix gmelinii stand in Daxing′anling Mountains
图 3 不同变量对兴安落叶松林单木直径生长的影响
Figure 3. Effects of different variables on the diameter growth of natural Larix gmelinii Stand in Daxing'anling Mountains
表 1 大兴安岭天然落叶松林分和单木调查因子统计
Table 1. Descriptive statistics of stand- and tree-level of natural Lari xgmelinii stand in Daxing′anling Mountains
组别
Group变量
Variable最小值
Min.最大值
Max.平均值
Mean标准差
Std.变异系数
C.V./%林分 Stand (n = 18) 平均胸径 Mean DBH/cm 7.00 14.80 10.96 2.04 18.64 平均树高 Mean HT/m 6.66 19.55 12.19 2.73 22.44 株数 Number/(trees·hm−2) 567.00 3850.00 1887.69 777.63 41.19 蓄积 Volume/(m3·hm−2) 144.66 646.67 366.71 136.06 37.10 郁闭度 Crown density 0.40 0.85 0.67 0.08 11.94 Margalef指数 Ma 0.19 0.87 0.54 0.20 36.82 Simpson指数 Sp 0.05 0.73 0.50 0.19 37.66 Shannon-Wiener指数 Sw 0.21 0.57 0.37 0.11 28.46 Pielow指数 Pw 0.25 1.41 1.02 0.30 29.14 混交度 M 0.03 0.48 0.25 0.11 45.36 大小比 U 0.42 0.53 0.50 0.02 3.73 角尺度 W 0.49 0.59 0.54 0.03 5.20 竞争指数 H 0.24 0.34 0.29 0.02 7.48 单木 Tree (n = 50) 年龄 Age/a 20.00 234.00 70.60 42.88 60.74 胸径 DBH/cm 6.50 38.80 16.44 7.34 44.66 树高 HT/m 6.10 26.00 14.64 4.16 28.42 冠幅 CW/m 0.70 2.88 1.69 0.51 30.08 冠长 CL/m 2.60 21.00 6.94 3.03 43.66 Note: DBH: diameter at breast height; HT: the height of tree; CW: crown width; CL: crown length; Ma, Sp, Sw and Pw are the diversity indexes of Margalef, Simpson, Shannon-Wiener and Pielow, respectively;U, M, W and H are respectively the indexes of neighborhood comparison, species mingling, uniform angle and competition. 
下载: 导出CSV
表 2 基础模型参数估计及拟合优度
Table 2. The parameter estimation and goodness-of-fit for the base models
编号
NO.模型
Model公式
Formulation参数估计值
Parameter estimation拟合优度
Goodness-of-fita b c 调整系数
Ra2平均绝对误差
Bias/cm均方根误差
RMSE/cmB1 Mitscherlich $y=a(1-\mathrm{E}\mathrm{x}\mathrm{p}(-b{t}\left)\right)$ 23.41 0.01314 0.6285 2.47 3.17 B2 Schumacher $y=a\mathrm{E}\mathrm{x}\mathrm{p}(-b/{t})$ 18.70 23.93 0.5897 2.77 3.38 B3 Logistic $y=a/(1+b\mathrm{E}\mathrm{x}\mathrm{p}(-c{t}\left)\right)$ 18.75 5.755 0.04191 0.6243 2.42 3.17 B4 Gompertz $y=aE\mathrm{x}\mathrm{p}(-b\mathrm{E}\mathrm{x}\mathrm{p}(-c{t}\left)\right)$ 21.92 2.232 0.02345 0.6302 2.39 3.14 注：y表示胸径，t表示年龄，a、b和c为待估参数;
Note:y is DBH, t is age, and a、b、c are the estimated parameters.
下载: 导出CSV
表 3 基于林分、单木、物种多样性和结构参数的单木直径生长模型拟合优度
Table 3. Goodness-of-fit for diameter growth model with stand-level, tree-level, diversity and structure variables
模型
Model参数a
Parameter a参数b
Parameter b拟合优度
Goodness-of-fit方差分析
ANOVAHg DBH M W Dg Pw H 调整系数
Ra2平均绝对误差
Bias/cm均方根误差
RMSE/cmB1 0.6285 2.47 3.17 M0 √ √ √ 0.7431 2.11 2.69 < 0.0001 M1 √ √ √ √ 0.7711 1.96 2.52 M2 √ √ √ √ 0.7453 2.11 2.68 M3 √ √ √ √ 0.7456 2.04 2.66 M4 √ √ √ √ 0.7987 1.92 2.43 < 0.0001 M5 √ √ √ √ √ 0.7984 1.84 2.39 M6 √ √ √ √ √ 0.8245 1.67 2.23 < 0.0001 M7 √ √ √ √ √ 0.7603 2.01 2.58 M8 √ √ √ √ √ 0.802 0 1.88 2.41 M9 √ √ √ √ √ 0.7912 1.83 2.40 M10 √ √ √ √ √ 0.8163 1.73 2.27 M11 √ √ √ √ √ √ 0.8257 1.68 2.22 M12 √ √ √ √ √ √ 0.8186 1.72 2.28 M13 √ √ √ √ √ √ 0.7984 1.82 2.38 M14 √ √ √ √ √ √ 0.8264 1.68 2.21 < 0.0001 M15 √ √ √ √ √ √ √ 0.8265 1.68 2.21 0.4157 注：√表示包含该参数；
Note: √ indicates that the variable was included.
下载: 导出CSV
表 4 兴安落叶松单木直径生长混合效应模型筛选结果
Table 4. Selection results of diameter growth model of natural Larix gmelinii stand when using mixed-effects regression techniques
项目
Items模型
Model参数
Parameter参数个数
Parameter number赤池信息准则
AIC贝叶斯信息准则
BIC对数似然值
Log-likelihood似然比检验
LRTP值
P−value随机参数
Random parameterE1 None 9 8773.91 8824.23 −4377.96 E2 a1 10 5628.26 5684.17 −2804.79 3158.93 < 0.0001 E3 a2 10 5629.59 5685.50 −2804.79 E4 a3 10 5644.69 5700.61 −2812.35 E5 a4 10 5631.07 5686.99 −2805.54 E6 b1 10 5321.24 5377.15 −2650.62 E7 a1,a2 11 5630.26 5691.76 −2804.13 E8 a1,a4 11 5630.02 5691.53 −2804.01 0.2353 0.6276 E9 a3,a4 11 5632.58 5694.08 −2805.29 方差函数
Variance functionE10 幂函数 Power 11 5611.80 5673.31 −2794.90 18.4562 < 0.0001 E11 指数函数 Exp 11 5625.34 5686.85 −2801.67 E12 常数加幂函数
ConstPower12 5583.78 5650.88 −2778.89 30.0214 < 0.0001 时间序列结构
Autocorrelation structureE13 一阶自回归与滑动平均
ARMA(1,1)12 3108.50 3175.59 −1542.25 2505.31 < 0.0001 E14 高斯结构
corGaussian13 −426.16 −353.47 226.08 3536.65 < 0.0001
下载: 导出CSV
表 5 兴安落叶松单木直径生长基础模型、广义模型和混合模型的参数估计和拟合结果
Table 5. Parameter estimation and fitting results of base-, general- and mixed-model for natural Larix gmelinii stand
项目
Items参数
Parameter基础模型
Base model (B1)广义模型
General model (M14)混合模型
Mixed model (E14)固定参数
Fixed parametera or a0 23.41 (< 0.0001) 178.95 (< 0.0001) 11.70 (0.0478) a1 −0.001994 (< 0.0001) −0.5389 (0.0284) a2 −0.5328 (< 0.0001) 1.466 (0.0004) a3 0.7033 (< 0.0001) 0.4928 (0.0764) a4 0.1614 (< 0.0001) 0.3669 (0.0005) b or b0 0.01314 (< 0.0001) 0.005387 (< 0.0001) 0.05265 (0.0473) b1 1.1188 (< 0.0001) −1.068 (0.0003) b2 2.1024 (< 0.0001) −0.05710 (0.0953) 模型误差
Error variance$ {\sigma }^{2} $ 10.06 4.865 0.6647 随机参数
Random variance$ {\sigma }_{{\mathrm{a}}_{1}}^{2} $ 0.1854 方差函数
Variance function幂函数 Power −0.01106 时间序列结构
Autocorrelation structure范围 Range 2.143 块金值 Nugget 0.1000 拟合统计量
Goodness-of-fit确定系数 Ra2 0.6285 0.8264 0.9710 平均绝对误差 Bias/cm 2.47 1.68 0.69 均方根误差 RMSE/cm 3.17 2.21 0.94
下载: 导出CSV
表 6 兴安落叶松单木直径生长基础模型、广义模型和混合模型的独立样本检验
Table 6. Validation of base-, general- and mixed-model for natural Larix gmelinii stand in Daxing'anling Mountains
模型
Model调整系数
Ra2平均绝对误差
Bias/cm均方根误差
RMSE/cm基础模型
Base model [B1]0.6273 2.52 3.25 广义模型
General model [M14]0.8577 1.92 2.14 混合模型
Mixed model [E14]0.9838 0.56 0.72
下载: 导出CSV
-
[1] Lhotka J M, Loewenstein E F. An individual-tree diameter growth model for managed uneven-aged oak-shortleaf pine stands in the Ozark Highlands of Missouri, USA[J]. Forest Ecology and Management, 2011, 261(3): 770-778. doi: 10.1016/j.foreco.2010.12.008 [2] Uzoh F C C, Oliver W W. Individual tree diameter increment model for managed even-aged stands of ponderosa pine throughout the western United States using a multilevel linear mixed effects model[J]. Forest Ecology and Management, 2008, 256(3): 438-445. doi: 10.1016/j.foreco.2008.04.046 [3] 孟宪宇. 测树学(第三版)[M]. 北京: 中国林业出版社, 2006. [4] 彭 娓, 李凤日, 董利虎. 黑龙江省长白落叶松人工林单木生长模型[J]. 南京林业大学学报: 自然科学版, 2018, 42(3):19-27. [5] 王 蒙, 李凤日. 基于抚育间伐效应的长白落叶松人工林单木直径生长模型[J]. 南京林业大学学报: 自然科学版, 2018, 42(3):28-36. [6] Vannoppen A, Kint V, Ponette Q, et al. Tree species diversity impacts average radial growth of beech and oak trees in Belgium, not their long-term growth trend[J]. Forest Ecosystems, 2019, 6: article id 10. doi: 10.1186/s40663-019-0169-z [7] Wang J, Cheng Y X, Zhang C Y, et al. Relationships between tree biomass productivity and local species diversity[J]. Ecosphere, 2016, 7(11): e01562. [8] Liang J J, Crowther T W, Picard N, et al. Positive biodiversity-productivity relationship predominant in global forests[J]. Science, 2016, 354(6309): article id aaf8957. doi: 10.1126/science.aaf8957 [9] 谭凌照, 范春雨, 范秀华. 吉林蛟河阔叶红松林木本植物物种多样性及群落结构与生产力的关系[J]. 植物生态学报, 2017, 41(11):29-36. [10] 汤孟平, 徐文兵, 陈永刚, 等. 天目山近自然毛竹林空间结构与生物量的关系[J]. 林业科学, 2011, 47(8):1-6. doi: 10.11707/j.1001-7488.20110801 [11] 吕延杰, 杨 华, 张 青, 等. 云冷杉天然林林分空间结构对胸径生长量的影响[J]. 北京林业大学学报, 2017, 39(9):41-47. [12] 姜立春, 杜书立. 基于非线性混合模型的东北兴安落叶松树高和直径生长模拟[J]. 林业科学研究, 2012, 25(1):11-16. doi: 10.3969/j.issn.1001-1498.2012.01.003 [13] Pinheiro J C, Bates D M. Mixed-effects Models in S and S-PLUS[M]. New York: Springer, 2001. [14] 祖笑锋, 倪成才, Nigh G, 等. 基于混合效应模型及EBLUP预测美国黄松林分优势木树高生长过程[J]. 林业科学, 2015, 51(3):28-36. [15] 张会儒, 李凤日, 赵秀海, 等. 东北过伐林可持续经营技术[M]. 北京: 中国林业出版社. 2016. [16] Dãnescu A, Albrecht A T, Bauhus J. Structural diversity promotes productivity of mixed, uneven-aged forests in southwestern Germany[J]. Oecologia, 2016, 182(2): 319-333. doi: 10.1007/s00442-016-3623-4 [17] 赵中华, 惠刚盈. 21世纪以来我国首创的森林经营方法[J]. 北京林业大学学报, 2019, 41(12):50-57. doi: 10.12171/j.1000-1522.20190325 -
点击查看大图
图(3) / 表(6)
计量
- 文章访问数: 3835
- HTML全文浏览量: 2516
- PDF下载量: 76
- 被引次数: 0
