高级检索

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于混合效应的杂种落叶松人工幼龄林单木枯损模型

王涛 董利虎 李凤日

引用本文:
Citation:

基于混合效应的杂种落叶松人工幼龄林单木枯损模型

    作者简介: 王涛,博士生。主要研究方向:林分生长模型。Email:568716463@qq.com.cn 地址:150040 黑龙江省哈尔滨市香坊区和兴路26号东北林业大学林学院.
    通讯作者: 李凤日,教授,博士生导师。主要研究方向:林分生长模型。Email:fengrili@126.com 地址:同上
  • 基金项目:

    “十二五”国家科技支撑计划课题 2015BAD09B01

  • 中图分类号: S758.1

Individual tree mortality model for hybrid larch young plantations based on mixed effects

  • 摘要: 目的 利用固定间隔期复测数据,运用不同方法建立杂种落叶松人工幼龄林单木枯损模型,为确定杂种落叶松合理的经营措施和推广应用提供依据。方法 基于2003—2015年黑龙江省江山娇实验林场48块样地的复测数据,通过Logistic模型,利用全子集法和最大似然估计构建杂种落叶松单木枯损模型。使用列联表分析和分类率-阈值散点图,确定枯损模型预估时的最佳阈值。引入随机参数,构建样地水平广义线性混合模型。模型估计方法为自适应积分最大似然估计,模型筛选指标为Akaike信息标准(AIC)、贝叶斯信息标准(BIC)以及-2倍对数似然值。通过计算绝对平均偏差(Bias),绘制ROC曲线以及模型预估枯损率与实际枯损率直方图对两种模型的预测结果进行评价比较。结果 包含单木(林木胸径,DBH;胸径平方,DBH2)、林分(林分断面积,BA)、竞争(大于对象木树木断面积之和变形,BALD)3个水平变量组合的单木枯损模型拟合效果最佳。杂种落叶松枯损主要发生在小径阶且相对竞争较大时。单木枯损概率随DBH增加逐渐减小,随BALD、BA增加而逐渐增加。最佳阈值有效提高了模型预估效果,方差-协方差结构为无结构矩阵(UN)时,四参数混合模型的拟合结果最佳,其预估的林分枯损率更接近实际林分枯损率。结论 混合模型能够更有效地描述和预估杂种落叶松的单木枯损。阈值分析是提高二分类模型预测准确性的有效方法。杂种落叶松作为速生树种,幼龄时期应适时进行抚育间伐以减少枯损发生的概率。
  • 图 1  杂种落叶松人工幼龄林不同径阶枯损株数分布

    Figure 1.  Distribution for the number of mortality trees in different DBH classes of hybrid larch young plantations

    图 2  广义线性模型与广义线性混合模型ROC曲线比较

    Figure 2.  Comparison of the ROC curve between generalized linear model and generalized linear mixed model

    图 3  广义线性模型、广义线性混合模型预测枯损率与实际枯损率比较

    Figure 3.  Comparison of the actual mortality rate in the result of generalized linear model and generalized linear mixed model predicted

    图 4  杂种落叶松人工幼龄林枯损模型阈值点与分类率关系

    Figure 4.  Relationship of classification rate and threshold for the mortality model of hybrid larch young plantations

    表 1  杂种落叶松人工幼龄林基本因子统计表

    Table 1.  Statistics of basic characteristics about hybrid larch young plantations

    数据类别
    Data class
    样地数
    Sample plot number(N)
    林分平均胸径
    Mean DBH(Dg)/cm
    林分平均高
    Mean stand height(H)/m
    林分断面积/
    (m2·hm-2)
    Stand basal area (BA)/
    (m2·ha-1)
    每公顷存活株数/
    (株·hm-2)
    Living tree number per hectare(NHAL)/
    (plant·ha-1)
    存活树木蓄积/
    (m3·hm-2)
    Living tree volume per hectare(VOLL)/
    (m3·ha-1)
    每公顷枯死株数/
    (株·hm-2)
    Dead tree number per hectare(NHAD)/
    (plant·ha-1)
    枯死树木蓄积/
    (m3·hm-2)
    Dead tree volume per hectare(VOLD)/
    (m3·ha-1)
    建模数据
    Modeling data
    40 3.07~15.12
    (13.29±2.67)
    2.41~14.12
    (8.95±2.41)
    2.05~38.74
    (19.49±8.39)
    1775~4875
    (2707±758)
    6.19~200.82
    (95.60±4.65)
    0~1100
    (71±119)
    0~14.19
    (0.89±2.08)
    检验数据
    Testing data
    8 4.33~14.56
    (13.16±2.67)
    2.44~14.08
    (8.81±2.44)
    3.92~34.66
    (19.90±8.59)
    1642~4275
    (2633±627)
    10.27~184.18
    (98.10±5.18)
    0~475
    (50±88)
    0~10.99
    (0.63±1.96)
    注:括号内为平均值±标准差。Note: that in bracket is the mean value ± standard deviation.
    下载: 导出CSV

    表 2  阈值检验指标计算公式描述

    Table 2.  Formula description for the threshold test index

    检验指标Testing index 缩写Abbreviation 公式Formula
    真阳性率/灵敏度True positive rate/sensitivity TPR TPR=TP/(TP+FN)
    假阴性率False negative rate FNR FNR=FN/(TP+FN)=1-TPR
    真阴性率/特异度True negative rate/specificity TNR TNR=TN/(FP+TN)
    假阳性率False positive rate FPR FPR=FP/(FP+TN)=1-TNR
    错误分类率Mistake classification rate MCR MCR=FPR+FNR
    正确分类率Accuracy classification rate ACR ACR=(TP+TN)/(TP+TN+FP+FN)
    下载: 导出CSV

    表 3  杂种落叶松人工幼龄林不同间隔期建模结果

    Table 3.  Modeling results with different intervals for hybrid larch young plantations

    间隔期/a Interval/year AUC
    1 0.856
    2 0.861
    3 0.870
    4 0.836
    5 0.845
    下载: 导出CSV

    表 4  杂种落叶松人工幼龄林全子集模型筛选结果

    Table 4.  Results of model selection with All-sets for hybrid larch young plantations

    模型编号
    Model No.
    变量Variable AIC BIC -2倍对数似然值
    -2 log likelihood
    AUC VIF
    1 BALD 7170.6 7186.2 7166.6 0.795 2.71
    2 NHAL BALD 7053.3 7061.5 7033.2 0.806 2.86
    3 DBH DBH2 BALD 6903.7 6934.8 6895.7 0.826 3.16
    4 DBH DBH2 BALD BA 6804.3 6843.3 6794.3 0.870 4.43
    5 DBH DBH2 NHAL BALD RD 6932.9 6923.7 6918.2 0.865 4.60
    6 DBH DBH2 BA NHAL BALD DBA 6931.3 6939.7 6932.5 0.852 4.87
    7 DBH DBH2 BA NHAL BALD DBA RD 6921.3 6918.8 6912.3 0.852 5.26
    8 DBH DBH2 BA NHAL BAL BALD RD DBA 6885.1 6879.2 6875.3 0.854 6.02
    注:DBH为林木胸径,BAL指所有大于对象木的树木断面积之和,BALD为BAL与DBH的变形,RD为相对直径, DBA为胸径和林分断面积之比,BA为林分断面积,NHAL为每公顷存活株数。下同。Notes: DBH, tree diameter at breast height; BAL, sum of stand basal area of all greater than object wood; BALD, transformation of BAL and DBH; RD, relative diameter; DBA, ratio of DBH to stand basal area; BA, stand basal area; NHAL, NHAL is number of living trees per hectare. The same below.
    下载: 导出CSV

    表 5  杂种落叶松人工幼龄林Logistic广义线性模型参数估计

    Table 5.  Parameter estimation for Logistic generalized linear model of hybrid larch young plantations

    变量Variable 参数Parameter 自由度Freedom 估计值
    Estimated value
    标准误差SD Wald χ2 P
    截距Intercept b0 1 -2.5321 0.2282 123.0760 < 0.0001
    DBH b1 1 -0.6185 0.0648 90.9819 < 0.0001
    DBH2 b2 1 0.0260 0.0030 74.1734 < 0.0001
    BALD b3 1 0.3411 0.0338 101.9940 < 0.0001
    BA b4 1 0.0238 0.0205 15.3487 < 0.0001
    下载: 导出CSV

    表 6  杂种落叶松人工幼龄林广义线性混合模型模拟

    Table 6.  Simulation of the Logistic generalized linear mixed model of hybrid larch young plantations

    模型编号
    Model No.
    参数个数
    Number of parameters
    随机参数
    Random parameter
    筛选指标
    Screening criteria
    b0 b1 b2 b3 b4 AIC BIC -2倍对数似然值
    -2 log likelihood
    1 5 6804.30 6843.30 6794.30
    2 6 6313.60 6323.74 6301.60
    3 6 6282.87 6293.00 6270.87
    4 6 6420.65 6430.78 6408.65
    5 6 6523.29 6533.43 6511.29
    6 6 6392.90 6403.03 6380.90
    7 8 6240.89 6254.40 6224.89
    8 8 6247.94 6261.45 6231.94
    9 8 6259.96 6273.47 6243.96
    10 8 6235.31 6248.83 6219.31
    11 8 6259.53 6273.04 6243.53
    12 8 6281.58 6295.09 6265.58
    13 8 6215.97 6212.46 6210.46
    14 8 6349.88 6363.40 6333.88
    15 11 6284.48 6292.93 6274.48
    16 11 6203.66 6212.10 6193.66
    17 11 6205.95 6214.40 6195.95
    18 11 6208.55 6217.00 6198.55
    19 11 6214.05 6222.50 6204.05
    20 11 6225.70 6234.14 6215.70
    21 11 6271.29 6279.73 6261.29
    22 11 6242.72 6251.17 6232.72
    23 15 6203.10 6211.54 6193.10
    23 15 6204.59 6213.03 6194.59
    24 15 6207.70 6216.14 6197.70
    25 15 6186.54 6194.98 6176.54
    26 20 6236.21 6244.32 6226.13
    注:▲表示包括这个随机参数;b0b1b2b3b4分别为截矩、DBH、DBH2、BALD、BA的模型参数。Notes:▲ means it is included in the random parameters; b0, b1, b2, b3, b4 are model parameters of intercept, DBH, DBH2, BALD, BA.
    下载: 导出CSV

    表 7  基于不同方差-协方差结构的模型比较

    Table 7.  Comparison of model based on different variance-covariance structures

    模型编号
    Model No.
    方差-协方差结构
    Variance-covariance structure
    参数个数
    Number of parameters
    AIC BIC -2倍对数似然值
    -2 log likelihood
    Pearson χ2
    25 无结构矩阵Unstructured matrix(UN) 15 6186.54 6194.98 6176.54 0.93
    对角矩阵Diagonal matrix(VC) 9 6270.39 6285.59 6252.39 0.88
    复合矩阵Compound symmetry(CS) 7 6362.29 6374.11 6348.29 0.91
    一阶自回归结构First-order autoregressive(AR(1)) 16 6221.74 6225.09 6211.44 0.92
    一阶移动平移结构First-order moving-average(MA(1))
    一阶自回归移动平均结构
    First-order autoregressive moving-average(ARMA(1, 1))
    下载: 导出CSV

    表 8  杂种落叶松人工幼龄林枯损模型拟合与检验结果

    Table 8.  Fitting and testing results of mortality model of hybrid larch young plantations

    固定参数Fixed parameter 基础模型Basic model 混合模型Mixed model
    b0 -2.5321 -2.6368
    b1 -0.6185 -1.4545
    b2 0.0260 0.0392
    b3 0.3411 0.3456
    b4 0.0238 0.0522
    随机效应方差-协方差结构Random effect variance-covariance structure(G) $\left[\begin{array}{cccc}{0.2297} & {-0.0034} & {0.0851} & {-0.0945} \\ {-0.0034} & {0.0001} & {-0.0001} & {0.0012} \\ {0.0851} & {-0.0001} & {0.0357} & {-0.0373} \\ {-0.0945} & {0.0012} & {-0.0373} & {0.0405}\end{array}\right]$
    Pearson χ2 1.36 0.93
    AUC 0.8701 0.9178
    Wald χ2 < 0.0001 < 0.0001
    Bias 5.375 1.375
    下载: 导出CSV

    表 9  杂种落叶松人工幼龄林枯损模型阈值预测列联表分析

    Table 9.  Confusion matrix of the hybrid larch mortality model for four cut points

    列联表分析
    Confusion matrix
    实际变量分类Actual variety classification
    Positive 0=枯损Dead Negative 1=存活Living
    预测变量分类
    Classification of predicted variety
    Positive
    0=枯损Dead
    TPR 83.6%(A=0.08) FPR 17.4%(A=0.08)
    82.3%(B=0.06) 17.9%(B=0.06)
    32.4%(C=0.5) 0.8%(C=0.5)
    86.1%(D=0.1) 10%(D=0.1)
    Negative
    1=存活Living
    FNR 16.4%(A=0.08) TNR 82.6%(A=0.08)
    17.7%(B=0.06) 82.1%(B=0.06)
    67.6%(C=0.5) 99.2%(C=0.5)
    13.9%(D=0.1) 90.0%(D=0.1)
    注:ABCD为阈值点。Note:A, B, C, D are threshold points.
    下载: 导出CSV

    表 10  杂种落叶松人工幼龄林枯损模型不同阈值预测分类率

    Table 10.  Classification of the hybrid larch young plantation mortality model for four cut points

    阈值Threshold ACR MCR
    A=0.08 82.9% 33.7%
    B=0.06 82.1% 35.7%
    C=0.5 81.2% 68.4%
    D=0.1 88.9% 23.9%
    注:ACR(正确分类率):正确预测枯损与存活株数占总株数比率; MCR(错误分类率):假阳性率与假阴性率之和。Notes:ACR(accurate classification rate): correctly predicting the ratio of mortality and living trees to total trees; MCR(misclassification rate): the sum of false positive rate and false negative rate.
    下载: 导出CSV
  • [1] Lee Y. Predicting mortality for even-aged stands of lodgepole pine[J]. Forestry Chronicle, 1971, 47(1): 29-32. doi: 10.5558/tfc47029-1
    [2] Monserud R A. Simulation of forest tree mortality[J]. Forest Science, 1976, 22(4): 438-444.
    [3] Dobbertin M, Biging G S. Using the non-parametric classifier CART to model forest tree mortality[J]. Forest Science, 1998, 44(4): 507-516.
    [4] Lutz J A, Halpern C B.Tree mortality during early forest development: a long-term study of rates, causes, and consequences[J]. Ecological Monographs, 2006, 76(2): 257-275. doi: 10.1890/0012-9615(2006)076[0257:TMDEFD]2.0.CO;2
    [5] Yang Y, Titus S J, Huang S. Modeling individual tree mortality for white spruce in Alberta[J]. Ecological Modelling, 2003, 163(3): 209-222. doi: 10.1016/S0304-3800(03)00008-5
    [6] Monserud R A, Sterba H. Modeling individual tree mortality for Austrian forest species[J]. Forest Ecology & Management, 1999, 113(2-3): 109-123.
    [7] 杜纪山.落叶松林木枯损模型[J].林业科学, 1999, 35(2): 45-49. doi: 10.3321/j.issn:1001-7488.1999.02.008Du J S. Tree mortality model of Larix[J]. Scientia Silvae Sinicae, 1999, 35(2): 45-49. doi: 10.3321/j.issn:1001-7488.1999.02.008
    [8] Weiskittel A R, Hann D W, Kershaw J A J, et al. Forest growth and yield modeling[J]. Forest Growth & Yield Modeling, 2002, 7(2): 223-233.
    [9] Das A J. Quantifying tree mortality risk and spatial pattern in a temperate conifer forest[D]. Reston: Dissertations & Theses-Gradworks, 2007.
    [10] Yao X H, Titus S J, Macdonald S E. A generalized logistic model of individual tree mortality for aspen, white spruce, and lodgepole pine in Alberta mixedwood forests[J]. Canadian Journal of Forest Research, 2001, 31(2): 283-291.
    [11] Eid T, Tuhus E. Models for individual tree mortality in Norway[J]. Forest Ecology & Management, 2001, 154(1-2): 69-84.
    [12] Temesgen H, Mitchell S J. An individual-tree mortality model for complex stands of southeastern British Columbia[J]. Western Journal of Applied Forestry, 2005, 20(2): 101-109.
    [13] Adame P, Río M D, Cañellas I. Modeling individual-tree mortality in Pyrenean oak (Quercus pyrenaica Willd.) stands[J]. Annals of Forest Science, 2010, 67(8): 810. doi: 10.1051/forest/2010046
    [14] Holzwarth F, Kahl A, Bauhus J, et al. Many ways to die-partitioning tree mortality dynamics in a near-natural mixed deciduous forest[J]. Journal of Ecology, 2013, 101(1): 220-230. doi: 10.1111/jec.2012.101.issue-1
    [15] 向玮, 雷相东, 刘刚, 等.近天然落叶松云冷杉林单木枯损模型研究[J].北京林业大学学报, 2008, 30(6): 90-98. doi: 10.3321/j.issn:1000-1522.2008.06.014Xiang W, Lei X D, Liu G, et al. Individual tree mortality models for semi-natural larch-spruce-fir forests in Jilin Province, northeastern China[J]. Journal of Beijing Forestry University, 2008, 30(6): 90-98. doi: 10.3321/j.issn:1000-1522.2008.06.014
    [16] Hamilton D A. Extending the range of applicability of an individual tree mortality model[J]. Canadian Journal of Forest Research, 1990, 20(8): 1212-1218. doi: 10.1139/x90-160
    [17] Misir M, Misir N, Yavuz H. Modeling individual tree mortality for crimean pine plantations[J]. Journal of Environmental Biology, 2007, 28(2): 167-172.
    [18] Hamilton D A, Edwards B M. Modeling the probability of individual tree mortality[Stand development][D]. Ogden: Usda Forest Service Research Paper Int, 1976.
    [19] 王济川, 郭志刚.Logistic回归模型:方法与应用[M].北京:高等教育出版社, 2001.Wang J C, Guo Z G. Logistic regression models: methods and application[M]. Beijing:Higher Education Press, 2001.
    [20] King S L, Bennett K P, List S. Modeling noncatastrophic individual tree mortality using logistic regression, neural networks, and support vector methods[J]. Computers & Electronics in Agriculture, 2000, 27(1-3): 401-406.
    [21] Ma Z H, Peng C, Li W Z, et al. Modeling individual tree mortality rates using marginal and margianl and random effects regression modes[J]. Natural Resource Modeling, 2013, 26(2): 131-153. doi: 10.1111/j.1939-7445.2012.00124.x
    [22] Schabenberger O, Pierce F J.Contemporary statistical models for the plant and soil science[M]. Boca Raton: CRC Press, 2001.
    [23] Mailly D, Gaudreault M, Picher G, et al. A comparison of mortality rates between top height trees and average site trees[J]. Annals of Forest Science, 2009, 66(2): 202. doi: 10.1051/forest/2008084
    [24] 姜立春, 李凤日.混合效应模型在林业建模中的应用[M].北京:科学出版社, 2014.Jiang L C, Li F R. Application of mixed effects models in forestry modeling[M]. Beijing:Science Press, 2014.
    [25] Groom J D, Hann D W, Temesgen H. Evaluation of mixed-effects models for predicting Douglas-fir mortality[J]. Forest Ecology & Management, 2012, 276(4): 139-145.
    [26] Timilsina N, Staudhammer C L. Individual tree mortality model for slash pine in Florida: a mixed modeling approach[J]. Southern Journal of Applied Forestry, 2012, 36(4): 211-219. doi: 10.5849/sjaf.11-026
    [27] 王涛, 董利虎, 李凤日.杂种落叶松人工幼龄林林分枯损规律及枯损模型[J].东北林业大学学报, 2017, 45(5): 39-43. doi: 10.3969/j.issn.1000-5382.2017.05.008Wang T, Dong L H, Li F R. Mortality of stand trees for hybrid larch young plantation in Heilongjiang[J]. Journal of Northeast Forestry University, 2017, 45(5): 39-43. doi: 10.3969/j.issn.1000-5382.2017.05.008
    [28] Hein S, Weiskittel A R. Cutpoint analysis for models with binary outcomes: a case study on branch mortality[J]. European Journal of Forest Research, 2010, 129(4): 585-590. doi: 10.1007/s10342-010-0358-3
    [29] 祖笑锋, 倪成才, Gorden Nigh, 等.基于混合效应模型及EBLUP预测美国黄松林分优势木树高生长过程[J].林业科学, 2015, 51(3): 25-33.Zu X F, Ni C C, Gorden N, et al. Based on mixed-effects model and empirical best linear unbiased predictor to predict growth profile of dominant height[J]. Scientia Silvae Sinicae, 2015, 51(3): 25-33.
    [30] Müller J, Hothorn T. Maximally selected two-sample statistics as a new tool for the identification and assessment of habitat factors with an application to breeding-bird communities in oak forests[J]. European Journal of Forest Research, 2004, 123(3): 219-228. doi: 10.1007/s10342-004-0035-5
    [31] Hosmer D W, Lemeshow S. Applied logistic regression[M]. New York:John Wiley & Sons, 2000.
    [32] Bravooviedo A, Sterba H, Mdel R, et al. Competition-induced mortality for Mediterranean Pinus pinaster Ait. and P. sylvestris L.[J]. Forest Ecology & Management, 2006, 222(1): 88-98.
    [33] Chen H, Fu S, Roberta M, et al. Relative size and stand age determine Pinus banksiana mortality[J]. Forest Ecology & Management, 2008, 255(12): 3980-3984.
    [34] Crecente-campo F, Marshall P, Rodríguezsoalleiro R. Modeling non-catastrophic individual-tree mortality for Pinus radiata plantations in northwestern Spain[J]. Forest Ecology & Management, 2009, 257(6): 1542-1550.
    [35] Bolker B M, Brooks M E, Clark C J, et al. Generalized linear mixed models: a practical guide for ecology and evolution[J]. Trends in Ecology & Evolution, 2009, 24(3): 127-135.
    [36] 张含国, 张成林, 兰士波, 等.落叶松杂种优势分析及家系选择[J].南京林业大学学报(自然科学版), 2005, 29(3): 69-72. doi: 10.3969/j.issn.1000-2006.2005.03.017Zhang H G, Zhang C L, Lan S B, et al. The analysis and selection for families of hybrid larch[J]. Journal of Nanjing Forestry University(Natural Sciences Edition), 2005, 29(3): 69-72. doi: 10.3969/j.issn.1000-2006.2005.03.017
    [37] 杨书文, 鞠永贵, 张世英, 等.落叶松杂种优势的研究[J].东北林业大学学报, 1985, 13(1): 30-36.Yang S W, Jin Y G, Zhang S Y, et al. Research on larch hybrid vigor[J]. Journal of Northeast Forestry University, 1985, 13(1): 30-36.
  • [1] 向玮雷相东刘刚徐光陈光法 . 近天然落叶松云冷杉林单木枯损模型研究. 北京林业大学学报, 2008, 5(6): 90-98.
    [2] 牛亦龙董利虎李凤日 . 基于广义代数差分法的长白落叶松人工林地位指数模型. 北京林业大学学报, 2020, 42(2): 9-18. doi: 10.12171/j.1000-1522.20190036
    [3] 王冬至张冬燕张志东黄选瑞 . 塞罕坝华北落叶松人工林断面积预测模型. 北京林业大学学报, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072
    [4] 李耀翔姜立春 . 基于非线性混合模型的落叶松木材管胞长度模拟. 北京林业大学学报, 2013, 10(3): 18-23.
    [5] . 基于FORECAST模型的长白落叶松人工林经营措施对长期生产力的影响. 北京林业大学学报, 2012, 9(6): 1-6.
    [6] 徐奇刚雷相东国红李海奎李玉堂 . 基于多层感知机的长白落叶松人工林林分生物量模型. 北京林业大学学报, 2019, 41(5): 97-107. doi: 10.13332/j.1000-1522.20190035
    [7] 黄俊吴普特赵西宁 , . 多参数非线性降雨产流阈值模型试验研究. 北京林业大学学报, 2011, 8(1): 84-89.
    [8] 余黎雷相东王雅志杨英军王全军 . 基于广义可加模型的气候对单木胸径生长的影响研究. 北京林业大学学报, 2014, 11(5): 22-32. doi: 10.13332/j.cnki.jbfu.2014.05.007
    [9] 臧颢雷相东张会儒李春明卢军 . 红松树高-胸径的非线性混合效应模型研究. 北京林业大学学报, 2016, 13(6): 8-9. doi: 10.13332/j.1000-1522.20160008
    [10] 李春明李利学 . 基于非线性混合模型的栓皮栎树高与胸径关系研究. 北京林业大学学报, 2009, 6(4): 7-12.
    [11] 李春明 . 利用非线性混合模型进行杉木林分断面积生长模拟研究. 北京林业大学学报, 2009, 6(1): 44-49.
    [12] 毛斌彭立群李乐徐程扬 . 侧柏风景林美景度的林内色彩斑块非线性模型研究. 北京林业大学学报, 2015, 12(7): 68-75. doi: 10.13332/j.1000-1522.20140481
    [13] 陈瑜徐程扬李乐蔡丽丽 . 阔叶红松风景林单木景观质量评价与模型研究. 北京林业大学学报, 2014, 11(5): 87-93. doi: 10.13332/j.cnki.jbfu.2014.05.006
    [14] 姚丹丹徐奇刚闫晓旺李玉堂 . 基于贝叶斯方法的蒙古栎林单木枯死模型. 北京林业大学学报, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260
    [15] 周永学毛俊娟魏潇潇胡胜华刘杏娥高黎张洪江吴彩燕袁怀文颜绍馗殷亚方杨平白岗栓黄荣凤郑小贤王费新胡万良何亚平张莉俊李瑞王芳张璧光邓小文秦爱光戴思兰王胜华王兆印樊军锋费世民崔赛华李猛王正汪思龙谭学仁王小青孙向阳乔建平张克斌罗晓芳NagaoHirofumi赵天忠常旭杜社妮张岩刘燕王晓欢KatoHideo李昀高荣孚张双保陈放龚月桦范冰江玉林王海燕张占雄韩士杰李华徐嘉张旭刘云芳江泽慧孔祥文刘秀英侯喜录常亮陈宗伟陈秀明杨培华李媛良丁磊李晓峰郭树花IdoHirofumi任海青薛岩高建社张代贵陈学平张桂兰李考学徐庆祥蒋俊明费本华金鑫刘永红李雪峰续九如王晓东涂代伦丁国权张红丽 , . 非线性植被-侵蚀动力学模型初探. 北京林业大学学报, 2007, 4(6): 123-128.
    [16] 王岩孙宇瑞冶民生谢响明何磊蒋佳荔李绍才张学俭罗菊春侯旭柳新伟张文娟张金凤李云成朱妍高鹏盖颖贺庆棠王盛萍李永慈李吉跃吕建雄申卫军何静关文彬张华丽崔保山孙海龙廖学品唐守正王文棋昌明成仿云冯仲科张志强康向阳陆佩玲吴玉英马道坤李小飞于晓南石碧杨志荣王军辉张桂莲蒋湘宁关毓秀吴斌静洁路婷张平冬史剑波何权孙阁赵广杰陈永国王尚德蒲俊文张满良孙晓霞马克明彭少麟汪燕赵燕东胡文忠余新晓刘国华林威汪西林 , . 马尾松人工林直径分布神经网络模型研究. 北京林业大学学报, 2006, 3(1): 28-31.
    [17] 王曼霖董利虎李凤日 . 基于Possion回归混合效应模型的长白落叶松一级枝数量模拟. 北京林业大学学报, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204
    [18] 娄明华张会儒雷相东卢军 . 天然云冷杉针阔混交林单木胸径树高空间自回归模型研究. 北京林业大学学报, 2016, 13(8): 1-9. doi: 10.13332/j.1000-1522.20150491
    [19] 张西贾黎明张瑜郑聪慧 . 基于FVS的秦岭地区栓皮栎天然次生林单木模型构建. 北京林业大学学报, 2015, 12(5): 19-29. doi: 10.13332/j.1000-1522.20140356
    [20] 曹梦潘萍欧阳勋志臧颢吴自荣杨阳占常燕 . 基于哑变量的闽楠天然次生林单木胸径和树高生长模型研究. 北京林业大学学报, 2019, 41(5): 88-96. doi: 10.13332/j.1000-1522.20190026
  • 加载中
图(4)表(10)
计量
  • 文章访问数:  784
  • HTML全文浏览量:  89
  • PDF下载量:  12
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-12-06
  • 录用日期:  2018-07-23
  • 刊出日期:  2018-10-01

基于混合效应的杂种落叶松人工幼龄林单木枯损模型

    通讯作者: 李凤日, fengrili@126.com
    作者简介: 王涛,博士生。主要研究方向:林分生长模型。Email:568716463@qq.com.cn 地址:150040 黑龙江省哈尔滨市香坊区和兴路26号东北林业大学林学院
  • 东北林业大学林学院,黑龙江 哈尔滨 150040
基金项目:  “十二五”国家科技支撑计划课题 2015BAD09B01

摘要: 目的利用固定间隔期复测数据,运用不同方法建立杂种落叶松人工幼龄林单木枯损模型,为确定杂种落叶松合理的经营措施和推广应用提供依据。方法基于2003—2015年黑龙江省江山娇实验林场48块样地的复测数据,通过Logistic模型,利用全子集法和最大似然估计构建杂种落叶松单木枯损模型。使用列联表分析和分类率-阈值散点图,确定枯损模型预估时的最佳阈值。引入随机参数,构建样地水平广义线性混合模型。模型估计方法为自适应积分最大似然估计,模型筛选指标为Akaike信息标准(AIC)、贝叶斯信息标准(BIC)以及-2倍对数似然值。通过计算绝对平均偏差(Bias),绘制ROC曲线以及模型预估枯损率与实际枯损率直方图对两种模型的预测结果进行评价比较。结果包含单木(林木胸径,DBH;胸径平方,DBH2)、林分(林分断面积,BA)、竞争(大于对象木树木断面积之和变形,BALD)3个水平变量组合的单木枯损模型拟合效果最佳。杂种落叶松枯损主要发生在小径阶且相对竞争较大时。单木枯损概率随DBH增加逐渐减小,随BALD、BA增加而逐渐增加。最佳阈值有效提高了模型预估效果,方差-协方差结构为无结构矩阵(UN)时,四参数混合模型的拟合结果最佳,其预估的林分枯损率更接近实际林分枯损率。结论混合模型能够更有效地描述和预估杂种落叶松的单木枯损。阈值分析是提高二分类模型预测准确性的有效方法。杂种落叶松作为速生树种,幼龄时期应适时进行抚育间伐以减少枯损发生的概率。

English Abstract

  • 树木枯损指在环境干扰和遗传特性共同作用下发生的树木生命力逐渐减弱直至死亡的过程[1-2],它受多种因素协同作用,且具有一定的物种特异性,是一个高度随机的过程[3]。根据树木死亡原因不同,枯损可分为非自然枯损和自然枯损两种[4-5]。非自然枯损包括由虫害、火灾、干旱等一系列突发原因造成的树木死亡;自然枯损包括对光、湿度和养分的竞争导致的死亡,自然枯损也称为竞争枯损[6],且相对而言更容易预测。单木枯损或存活模型是林分生长和收获预估模型系统的重要组成部分[7-8],对描述林分动态有着重要的意义[9]。树木的寿命决定了碳存在的时间以及森林生物量中碳汇量的大小,因此理解预测树木枯损是森林生态、生物地球化学、林分动态分析中不可或缺的过程。随着森林集约化经营管理水平的提高,特别是人工商品林的生产和经营管理,使得人们更加关注林木在未来时刻的生长和存活状态,枯损模型的研究为造林初值密度设计、间伐、主伐等营林生产活动提供了有力的依据[7]

    在枯损研究中应用最广的是单木枯损模型[10-15],通过以单木为研究对象计算其存活或枯损的概率,同径阶和林分枯损模型相比较,具有更高的预估精度。单木枯损模型研究发现单木枯损同树木活力、单木竞争以及林分特性有着重要的关系。早期枯损研究中,常用Logistic回归构建单木枯损模型[16-17]。Logistic回归有着极佳的统计学特性,其预测值介于0~1之间,与单木枯损概率的分布范围相契合。枯损本质上是一个非高斯二项分布,所以Logistic回归模型采用最大似然估计(MLE)对模型进行估计,其模型评价、回归系数解释以及统计推断有其特定的形式[18-19]。随着数据收集以及计算机方法的提高,更多的方法被用来进行单木枯损模型的研究,包括神经网络、边际效应模型等[20-21]。传统单木枯损研究的主要问题之一是忽略了数据中的固有多层次结构,树木嵌套在样地当中,样地中的有效单木数据往往都是重复测量得到的,不同水平的数据具有嵌套性,数据之间存在着一定的交互作用[22],同一样地的树木会受到竞争和环境因素的协同作用,利用多层次数据进行建模,模型误差结果偏大。另一个主要问题是模型预估时未能给出合理的阈值,传统模型多选择阈值0.5[4-5]。如果一个案例的预测事件概率大于0.5,就将其界定为预测事件发生,否则认定该事件不发生。实际上0.5作为枯损发生与否的判别依据并不准确,该值往往不能对枯损事件进行精准的预估,使得模型预估精度受到限制[13]。利用混合模型解决单木枯损模型多水平数据自相关是单木枯损研究的一次新尝试,通过将固定参数和随机参数相结合,满足独立同分布的假设的同时[23],得到渐进无偏参数估计,并利用方差-协方差结构来反映数据间的相关性和异质性,大大减小了模型误差[24-26]。除了利用混合模型外,在模型预估过程中,通过ROC曲线和分类表相结合,确定更加合理的阈值,能够有效提高模型的预估精度,方便模型的应用。

    杂种落叶松(Larix kaempferi × Larix olgensis)是以当地种源的日本落叶松(Larix kaempferi)为亲本之一,长白落叶松(Larix olgensis)为另一亲本的属内种间杂交品种,其生长速度具有明显优势,为喜光的非耐阴树种[27],国内暂没有对于该树种的枯损研究。本文利用广义线性混合模型(GLMM),建立样地水平杂种落叶松枯损模型,通过阈值分析与列联表分析和传统单木枯损模型进行比较[28],有助于了解杂种落叶松枯损的主要规律,为杂种落叶松人工林动态预测和制定合理的经营措施提供基础,同时对东北地区推广杂种落叶松林具有重要现实意义。

    • 数据来自于1998年黑龙江省林业科学院江山娇实验林场营造的20hm2杂种落叶松人工林,造林地为蒙古栎(Quercus mongolica)天然林皆伐迹地。2003年黑龙江省林业科学研究所在实验林内设置48块固定标准地,造林密度从2000株/hm2到6000株/hm2,样地面积从0.03hm2到0.0621hm2(平均0.04hm2)。各标准地集中连片,立地条件基本相同。2015年时林龄为19年。自2003年设置杂种落叶松人工林固定标准地以来,每年对胸径、树高(分平均木和优势木测量5~7株单木)以及树木状态(枯损,设为1;存活,设为0)进行观测,截止2015年已持续12期观测,标准地至今未经过抚育间伐,数据保存完整。本文按照样地进行随机抽选,共抽选出48块样地。利用40块样地的数据进行模型构建,8块样地的数据进行模型检验。样地基本因子的统计信息见表 1。数据中枯损株数按径阶的分布信息见图 1

      表 1  杂种落叶松人工幼龄林基本因子统计表

      Table 1.  Statistics of basic characteristics about hybrid larch young plantations

      数据类别
      Data class
      样地数
      Sample plot number(N)
      林分平均胸径
      Mean DBH(Dg)/cm
      林分平均高
      Mean stand height(H)/m
      林分断面积/
      (m2·hm-2)
      Stand basal area (BA)/
      (m2·ha-1)
      每公顷存活株数/
      (株·hm-2)
      Living tree number per hectare(NHAL)/
      (plant·ha-1)
      存活树木蓄积/
      (m3·hm-2)
      Living tree volume per hectare(VOLL)/
      (m3·ha-1)
      每公顷枯死株数/
      (株·hm-2)
      Dead tree number per hectare(NHAD)/
      (plant·ha-1)
      枯死树木蓄积/
      (m3·hm-2)
      Dead tree volume per hectare(VOLD)/
      (m3·ha-1)
      建模数据
      Modeling data
      40 3.07~15.12
      (13.29±2.67)
      2.41~14.12
      (8.95±2.41)
      2.05~38.74
      (19.49±8.39)
      1775~4875
      (2707±758)
      6.19~200.82
      (95.60±4.65)
      0~1100
      (71±119)
      0~14.19
      (0.89±2.08)
      检验数据
      Testing data
      8 4.33~14.56
      (13.16±2.67)
      2.44~14.08
      (8.81±2.44)
      3.92~34.66
      (19.90±8.59)
      1642~4275
      (2633±627)
      10.27~184.18
      (98.10±5.18)
      0~475
      (50±88)
      0~10.99
      (0.63±1.96)
      注:括号内为平均值±标准差。Note: that in bracket is the mean value ± standard deviation.

      图  1  杂种落叶松人工幼龄林不同径阶枯损株数分布

      Figure 1.  Distribution for the number of mortality trees in different DBH classes of hybrid larch young plantations

    • 为建立杂种落叶松单木枯损混合效应模型,选取Logistic模型作为枯损的基础模型。Logistic模型是用来处理二分类数据最为有效的广义线性模型(GLM),其基本模型形式为:

      $ P=\frac{1}{1+\mathrm{e}^{-\left(\alpha+\beta X_{i}\right)}}=\frac{\mathrm{e}^{\alpha+\beta X_{i}}}{1+\mathrm{e}^{\alpha+\beta X_{i}}} $

      (1)

      式中:P为事件发生的概率,在本文中指的是树木3年间隔期内发生枯损的概率;Xi为预测变量矩阵;α为常数项;β为预测变量的系数向量。

      可以通过“logit”链接函数转换为线性形式:

      $ \ln \left(\frac{P}{1-P}\right)=\alpha+\boldsymbol{\beta} \boldsymbol{X}_{i} $

      (2)

      本文在此基础上,引入样地水平随机效应,构建基于Logistic回归的单水平广义线性混合模型(GLMM),模型最终形式可以表达为:

      $ $$\left\{ {\matrix{ {{\mathop{\rm logit}\nolimits} (P) = \alpha + {\boldsymbol{\beta }}{\boldsymbol{X}_i} + {\boldsymbol{Z}_i}{\boldsymbol{b}_i} + {\boldsymbol{e}_i}} \hfill \cr \begin{equation} \boldsymbol{b}_{i} \sim N(0, \boldsymbol{G}) \end{equation} \hfill \cr {{\boldsymbol{e}_i}\sim N\left( {0, {\sigma ^2}{\boldsymbol{R}_i}} \right)} \hfill \cr } } \right.$$ $

      (3)

      式中:α为常数项;βp×1为预测变量的系数向量;Xi(ni×p)和Zi(ni×q)分别代表对应于p个固定效应和q个随机效应的设计矩阵固定效应参数向量;bi代表对应于第i个观测对象的随机效应参数;eini×1维误差向量;Gq×q维随机效应协方差矩阵;Rini×ni维组内协方差矩阵;σ2为方差;ni为研究水平样本数。

    • 枯损为小概率事件,连年观测的数据可能无法准确反映林分的枯损规律,而国内外通常采用的5年间隔期对于速生性极强的杂种落叶松可能也不适用。本文通过将建模数据进行整理,按照1~5年间隔期,整理出5个不同的数据集并分别建模,为方便单木枯损模型的建立,将数据中的变量分为3类:(1)单木变量(I),包括DBH、DBH2;(2)林分变量(ST),包括BA、NHAL;(3)竞争变量(C),包括BAL、BALD、RD、DBA。利用全子集法(All sub-set regression),将8个备选变量进行组合,分别建立模型。其中BAL指所有大于对象木树木断面积之和,BALD为BAL与DBH的变形,RD为相对直径,DBA为胸径和林分断面积之比,具体见公式(4)~(7)。在模型选择时,应该考虑变量之间的多重共线性问题。多重共线性会使得模型的统计推断严重偏离实际。选择方差膨胀因子(VIF)作为多重共线性的检验指标,当VIF>10时,模型存在多重共线性。本文选择保留0 < VIF≤10的模型,最后通过模型拟合优度比较选出杂种落叶松单木枯损模型最佳间隔期和最佳基础模型形式。

      $ \begin{gathered} \mathrm{BAL}_{i}=\sum_{j=i}^{n} \pi\left(\frac{\mathrm{DBH}_{j}}{2}\right)^{2},\left(\mathrm{DBH}_{j} \in P\right), \hfill \\ P = \left\{ {{\text{DB}}{{\text{H}}_j}, j = 1, 2, \cdot \cdot \cdot , n|{\text{DB}}{{\text{H}}_j}{\text{ > DB}}{{\text{H}}_i}} \right\} \hfill \\ \end{gathered} $

      (4)

      $ \operatorname{DALD}_{i}=\frac{\mathrm{BAL}_{i}}{\ln \left(\mathrm{DBH}_{i}+1\right)} $

      (5)

      $ \mathrm{RD}_{i}=\frac{\mathrm{DBH}_{i}}{D_{\mathrm{g}}} $

      (6)

      $ \mathrm{DBA}_{i}=\frac{\mathrm{DBH}_{i}}{\mathrm{BA}} $

      (7)

      式中:DBHj为第j株树木的胸径,且大于对象木iDg为对象木i所在林分平均直径;BA表示对象木i所在林的林分断面积。

    • 单木枯损数据中存在着调查时间上的自相关性和不同样地水平上的异方差性,通过混合模型中的组内方差-协方差结构以及随机效应方差-协方差结构可以对上述两个问题进行修正。常采用公式(8)来描述时间相关性的组内方差-协方差结构(Ri)。本文利用3种组内方差-协方差结构来解决时间自相关问题,包括一阶自回归结构(AR(1))、一阶移动平移结构(MA(1))和一阶自回归移动平均结构(ARMA(1, 1))。利用3种随机效应参数方差-协方差结构对样地水平上的异质性进行描述,包括无结构矩阵(UN)、复合对称矩阵(CS)、对角矩阵(DM)等[24]。具体操作过程为:

      (1) 确定单木枯损的基础模型后,假定随机效应的方差-协方差结构为无结构矩阵(UN)。然后,利用SAS 9.4的PROC GLIMMIX模块拟合包含不同随机参数个数组合的模型,从中选出模拟精度最高的模型作为最佳的随机参数组合形式。

      (2) 通过比较拟合优度指标(AIC、BIC、-2倍对数似然值、Pearson χ2),确定步骤(1)中所得到的最佳模型的方差-协方差结构,包括预先假定的无结构矩阵(UN)、复合对称矩阵(CS)和对角矩阵(DM)。

      (3) 测试一阶自回归结构(AR(1))、一阶移动平移结构(MA(1))和一阶自回归移动平均结构(ARMA(1, 1))这3种组内方差-协方差结构对于混合模型的影响。

      $ \boldsymbol{R}_{i}={\sigma}^{2} \boldsymbol{G}_{i}^{0.5} \boldsymbol{T}_{i} \boldsymbol{G}_{i}^{0.5} $

      (8)

      式中:σ2为模型的误差方差值;Ti为组内误差相关性结构;Gi为描述方差异质性的对角矩阵。

    • Logistic回归模型通过SAS 9.4中PROC LOGISTIC过程进行拟合,广义线性混合模型通过SAS 9.4中PROC GLIMMIX过程进行拟合。PROC GLIMMIX过程默认的参数估计方法为限制性残差伪似然估计(RSPL),该方法不能提供真实的似然值。本文选用自适应积分最大似然估计(QUAD)作为参数估计方法,该方法能提供真实的似然值。通过计算Akaike信息标准(AIC)、Bayesian信息标准(BIC)以及-2倍对数似然值进行模型筛选。选择模型的标准是其值达到最小的自变量子集为最优子集,具体公式如下:

      $ - 2倍对数似然值 = - 2{\text{In}}(L) $

      (9)

      $ \mathrm{AIC}=2 p-2 \ln (L)=n \ln (\mathrm{SSE})-n \ln (n)+2 p $

      (10)

      $ \begin{gathered} {\text{BIC = - 2In(}}{L) + k}{\text{In(}}{n) = n}{\text{In(SSE) - }}{n}{\text{In(}}{n) + } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{In(}}{n)p} \hfill \\ \end{gathered} $

      (11)

      式中:L为估计模型的似然函数最大值,n为样本数,p为参数个数,SSE为模型残差平方和。

      本文采用的拟合优度检验方法包括模型参数显著性检验、Wald χ2检验、Pearson χ2检验以及受试者工作特征曲线(Receiver operating characteristic,简称ROC曲线)。Pearson χ2检验常用来检验数据是否存在过离散现象,χ2值越接近1,过离散程度越小。ROC曲线是以假阳性率为横坐标,真阳性率为纵坐标构建的曲线图,曲线下方面积为AUC值,其值可近似等价于传统线性回归中的R2,AUC值大于0.7时,说明模型符合精度要求[19]

      混合模型中固定效应部分的检验与传统的检验方法相同。然而,随机效应部分的检验需要二次抽样来计算随机参数值。随机效应值是采用EBLUP-经验线性无偏最优预测法(empirical best linear unbiased predictor)进行计算(公式(12)),其中$\hat{b}_{i}$为随机向量bi的条件分布数学期望;为随机向量bi协方差矩阵的估计值;$\boldsymbol{\hat{Z}}_{i}$和$\boldsymbol{\hat{Z}}_{i}^{\mathrm{T}}$分别表示随机效应设计矩阵(Zi)及其转置(ZiT)的估计值;$\boldsymbol{\tilde{y}}_{i}$表示观测值向量;$\hat{\boldsymbol{\beta}}_{i}$为固定效应参数向量估计值,A为其关联矩阵;Xi为预测变量矩阵;f($\boldsymbol{A} \hat{\boldsymbol{\beta}}_{i}$, Xi)为固定效应参数计算的预测值[29]

      模型检验指标为绝对平均偏差(Bias)(公式(13)),其中Yi为第i个样地的实际枯损株数,πij为模型第i个样地模型预测枯损株数,N为检验数据样地个数,具体定义参见文献[25]。通过绘制广义线性模型、广义线性混合模型预测枯损率与实际枯损率比较图,比较基础单木枯损混合模型对不同径阶数据在预测时的差异。

      $ \hat{b_{i}}=\hat{\boldsymbol{D}} \hat{\boldsymbol{Z}}_{i}^{\mathrm{T}}\left(\hat{\boldsymbol{Z}}_{i} \hat{\boldsymbol{D}} \hat{\boldsymbol{Z}}_{i}^{\mathrm{T}}\right)^{-1}\left[\tilde{\boldsymbol{y}}_{i}-f\left(\boldsymbol{A} \hat{\boldsymbol{\beta}}_{i}, \boldsymbol{X}_{i}\right)\right] $

      (12)

      $ {Bias}=\frac{\boldsymbol{\sum}\left|\boldsymbol{Y}_{i}-\boldsymbol{\pi}_{i j}\right|}{N} $

      (13)
    • 常利用阈值(threshold)或概率界限(cut point)将枯损概率转化为真实事件(枯损或存活),来决定枯损是否发生[30]。传统枯损模型多选择阈值0.5[4-5],如果单木枯损预测概率大于0.5,就将其界定为枯损,否则记为存活。理论上,该值仅适用于林分内存活木株数与枯损木株数相近的情况,实际上枯损是小概率事件,枯损和存活二者相等的情况几乎没有。阈值分析是提高模型预估的有效手段,目前较为流行的阈值确定方法包括固定阈值[6, 9, 25]、随机阈值[12]和割点阈值[6]。固定阈值通常认为正确分类率最大时的阈值为最佳阈值点;随机阈值常选用灵敏度与特异度曲线交点时的阈值为最佳阈值;割点阈值选取树种的林分平均观测枯损率作为最佳阈值。

      阈值分析的基本原理是利用二分类列联表,对实际枯损、存活株数和预测枯损、存活株数之间不同关系进行计算,利用相关检验标准进行阈值选择。二分类列联表也称为混淆矩阵,它给出了两种不同类型的错误分类和两种不同类型的正确分类。数学上定义事件发生为阳性(P),而未发生为阴性(N)。本文中枯损为阳性,存活为阴性。当预测结果与观测结果相同时记为真(T),预测结果与观测结果不同时记为假(F)。TP代表预测结果与实际观测结果均为枯损,称为真阳性事件;TN代表预测结果与实际观测结果均为存活,称为真阴性事件;FP代表预测结果为枯损而实际观测结果为存活,称为假阳性事件;FN代表预测结果为存活而观测结果为枯损,称为假阴性事件。具体阈值检验指标计算公式见表 2,其中TP与FN之和恰好是观测枯损总株数,FP与TN之和为观测存活总株数。本文通过SAS 9.4软件包随机产生1000个0~1的随机阈值,计算该阈值点下基础模型与混合模型预测的分类率。构建列联表以及分类率-阈值散点图,以灵敏度、特异度和总分类正确率作为评判标准,探究不同阈值点对枯损模型的意义,包括:(1)传统阈值0.5;(2)灵敏度同特异度相等时的阈值点[31];(3)错误分类率最低时的阈值点[28]

      表 2  阈值检验指标计算公式描述

      Table 2.  Formula description for the threshold test index

      检验指标Testing index 缩写Abbreviation 公式Formula
      真阳性率/灵敏度True positive rate/sensitivity TPR TPR=TP/(TP+FN)
      假阴性率False negative rate FNR FNR=FN/(TP+FN)=1-TPR
      真阴性率/特异度True negative rate/specificity TNR TNR=TN/(FP+TN)
      假阳性率False positive rate FPR FPR=FP/(FP+TN)=1-TNR
      错误分类率Mistake classification rate MCR MCR=FPR+FNR
      正确分类率Accuracy classification rate ACR ACR=(TP+TN)/(TP+TN+FP+FN)
    • 利用不同间隔期数据集建立模型,模型拟合结果见表 3。由表 3可以看出,使用3年间隔期建模,模型拟合AUC值最高。表 4列出了全子集法建模结果,自变量个数相同的组合仅列出AIC值最小、AUC值最大并且方差膨胀因子满足约束条件时的模型形式。从表 4可以看出,当自变量个数为1~4时,模型拟合AUC值逐渐增大,当参数个数继续增加时,AUC值逐渐减小。本文选用表 4中的模型4作为最优基础模型,此模型同时包含了单木、竞争和林分3个水平的变量,预估AUC值最大,且无多重共线性。具体模型见公式(14)。表 5为该模型的参数估计,从中可以看出各变量的χ2检验均显著。

      $ \begin{array}{c}{P=1 /\left(1+\exp \left[-\left(b_{0}+b_{1} \mathrm{DBH}+b_{2} \mathrm{DBH}^{2}+\right.\right.\right.} \\ {b_{3} \mathrm{BALD}+b_{4} \mathrm{BA} ) ] )}\end{array} $

      (14)

      表 3  杂种落叶松人工幼龄林不同间隔期建模结果

      Table 3.  Modeling results with different intervals for hybrid larch young plantations

      间隔期/a Interval/year AUC
      1 0.856
      2 0.861
      3 0.870
      4 0.836
      5 0.845

      表 4  杂种落叶松人工幼龄林全子集模型筛选结果

      Table 4.  Results of model selection with All-sets for hybrid larch young plantations

      模型编号
      Model No.
      变量Variable AIC BIC -2倍对数似然值
      -2 log likelihood
      AUC VIF
      1 BALD 7170.6 7186.2 7166.6 0.795 2.71
      2 NHAL BALD 7053.3 7061.5 7033.2 0.806 2.86
      3 DBH DBH2 BALD 6903.7 6934.8 6895.7 0.826 3.16
      4 DBH DBH2 BALD BA 6804.3 6843.3 6794.3 0.870 4.43
      5 DBH DBH2 NHAL BALD RD 6932.9 6923.7 6918.2 0.865 4.60
      6 DBH DBH2 BA NHAL BALD DBA 6931.3 6939.7 6932.5 0.852 4.87
      7 DBH DBH2 BA NHAL BALD DBA RD 6921.3 6918.8 6912.3 0.852 5.26
      8 DBH DBH2 BA NHAL BAL BALD RD DBA 6885.1 6879.2 6875.3 0.854 6.02
      注:DBH为林木胸径,BAL指所有大于对象木的树木断面积之和,BALD为BAL与DBH的变形,RD为相对直径, DBA为胸径和林分断面积之比,BA为林分断面积,NHAL为每公顷存活株数。下同。Notes: DBH, tree diameter at breast height; BAL, sum of stand basal area of all greater than object wood; BALD, transformation of BAL and DBH; RD, relative diameter; DBA, ratio of DBH to stand basal area; BA, stand basal area; NHAL, NHAL is number of living trees per hectare. The same below.

      表 5  杂种落叶松人工幼龄林Logistic广义线性模型参数估计

      Table 5.  Parameter estimation for Logistic generalized linear model of hybrid larch young plantations

      变量Variable 参数Parameter 自由度Freedom 估计值
      Estimated value
      标准误差SD Wald χ2 P
      截距Intercept b0 1 -2.5321 0.2282 123.0760 < 0.0001
      DBH b1 1 -0.6185 0.0648 90.9819 < 0.0001
      DBH2 b2 1 0.0260 0.0030 74.1734 < 0.0001
      BALD b3 1 0.3411 0.0338 101.9940 < 0.0001
      BA b4 1 0.0238 0.0205 15.3487 < 0.0001
    • 在假定随机效应的方差-协方差结构为无结构矩阵(UN)的基础上,拟合包含不同随机参数个数组合的模型,从中选出模拟精度最高的模型作为最佳的随机参数组合形式。表 6列出了不同随机参数组合的收敛拟合结果,共计26个模型。利用-2倍对数似然值、AIC、BIC值对表 6中模型进行比较,其值越小说明模型的拟合效果越好。所有混合模型的拟合效果均优于固定模型(表 6中模型1),此时的模型1即为不包含随机效应的基础模型。当仅考虑一个随机效应参数,模型3的AIC值最小,说明增加一个随机参数提高了模型拟合精度,此时随机参数为b1。加入两个随机参数组合时,模型13的拟合效果最佳,此时随机参数为b1b4。加入3个随机参数组合时,模型16拟合效果最佳,随机参数为b0b1b3。加入4个随机参数组合时,模型25拟合效果最佳,随机参数为b1b2b3b4。随着随机参数的增加,拟合效果逐渐增强,当随机参数增加为5个时,拟合效果开始下降,最终确定最佳随机参数组合为包含4个随机参数的模型25。

      表 6  杂种落叶松人工幼龄林广义线性混合模型模拟

      Table 6.  Simulation of the Logistic generalized linear mixed model of hybrid larch young plantations

      模型编号
      Model No.
      参数个数
      Number of parameters
      随机参数
      Random parameter
      筛选指标
      Screening criteria
      b0 b1 b2 b3 b4 AIC BIC -2倍对数似然值
      -2 log likelihood
      1 5 6804.30 6843.30 6794.30
      2 6 6313.60 6323.74 6301.60
      3 6 6282.87 6293.00 6270.87
      4 6 6420.65 6430.78 6408.65
      5 6 6523.29 6533.43 6511.29
      6 6 6392.90 6403.03 6380.90
      7 8 6240.89 6254.40 6224.89
      8 8 6247.94 6261.45 6231.94
      9 8 6259.96 6273.47 6243.96
      10 8 6235.31 6248.83 6219.31
      11 8 6259.53 6273.04 6243.53
      12 8 6281.58 6295.09 6265.58
      13 8 6215.97 6212.46 6210.46
      14 8 6349.88 6363.40 6333.88
      15 11 6284.48 6292.93 6274.48
      16 11 6203.66 6212.10 6193.66
      17 11 6205.95 6214.40 6195.95
      18 11 6208.55 6217.00 6198.55
      19 11 6214.05 6222.50 6204.05
      20 11 6225.70 6234.14 6215.70
      21 11 6271.29 6279.73 6261.29
      22 11 6242.72 6251.17 6232.72
      23 15 6203.10 6211.54 6193.10
      23 15 6204.59 6213.03 6194.59
      24 15 6207.70 6216.14 6197.70
      25 15 6186.54 6194.98 6176.54
      26 20 6236.21 6244.32 6226.13
      注:▲表示包括这个随机参数;b0b1b2b3b4分别为截矩、DBH、DBH2、BALD、BA的模型参数。Notes:▲ means it is included in the random parameters; b0, b1, b2, b3, b4 are model parameters of intercept, DBH, DBH2, BALD, BA.

      表 7列出了模型25在考虑3种不同方差-协方差结构时的结果比较,无结构矩阵(UN)作为方差-协方差结构时具有最小的AIC、BIC、-2倍对数似然值,且Pearson χ2更接近1,说明UN结构具有更好的拟合优度,并且所构建的模型不存在过离散情况。在尝试加入解决时间自相关问题的组内方差-协方差结构时,包括一阶自回归结构(AR(1))、一阶移动平移结构(MA(1))和一阶自回归移动平均结构(ARMA(1, 1)),混合模型拟合无法收敛或精度没有提高,考虑模型的简洁性,本文仅考虑随机效应部分的方差-协方差结构。

      表 7  基于不同方差-协方差结构的模型比较

      Table 7.  Comparison of model based on different variance-covariance structures

      模型编号
      Model No.
      方差-协方差结构
      Variance-covariance structure
      参数个数
      Number of parameters
      AIC BIC -2倍对数似然值
      -2 log likelihood
      Pearson χ2
      25 无结构矩阵Unstructured matrix(UN) 15 6186.54 6194.98 6176.54 0.93
      对角矩阵Diagonal matrix(VC) 9 6270.39 6285.59 6252.39 0.88
      复合矩阵Compound symmetry(CS) 7 6362.29 6374.11 6348.29 0.91
      一阶自回归结构First-order autoregressive(AR(1)) 16 6221.74 6225.09 6211.44 0.92
      一阶移动平移结构First-order moving-average(MA(1))
      一阶自回归移动平均结构
      First-order autoregressive moving-average(ARMA(1, 1))

      表 8列出了基本模型(表 6中模型)和混合模型(表 6中模型25)的拟合结果和检验结果,包括固定参数值、随机参数方差-协方差组成、拟合指标(Pearson χ2、Wald χ2、AUC)和检验指标(Bias)。结果表明:基于样地效应的混合模型优于基本模型,AUC值由0.8701提高到0.9178。χ2检验 < 0.0001说明模型有显著意义,Pearson χ2约等于1,模型过离散情况得到修正。模型检验偏差值Bias显示混合模型检验精度高于基础模型。图 2为广义线性模型与广义线性混合模型ROC曲线图。

      表 8  杂种落叶松人工幼龄林枯损模型拟合与检验结果

      Table 8.  Fitting and testing results of mortality model of hybrid larch young plantations

      固定参数Fixed parameter 基础模型Basic model 混合模型Mixed model
      b0 -2.5321 -2.6368
      b1 -0.6185 -1.4545
      b2 0.0260 0.0392
      b3 0.3411 0.3456
      b4 0.0238 0.0522
      随机效应方差-协方差结构Random effect variance-covariance structure(G) $\left[{array}{cccc}{0.2297} & {-0.0034} & {0.0851} & {-0.0945} \\ {-0.0034} & {0.0001} & {-0.0001} & {0.0012} \\ {0.0851} & {-0.0001} & {0.0357} & {-0.0373} \\ {-0.0945} & {0.0012} & {-0.0373} & {0.0405}{array}\right]$
      Pearson χ2 1.36 0.93
      AUC 0.8701 0.9178
      Wald χ2 < 0.0001 < 0.0001
      Bias 5.375 1.375

      图  2  广义线性模型与广义线性混合模型ROC曲线比较

      Figure 2.  Comparison of the ROC curve between generalized linear model and generalized linear mixed model

      图 3为两个模型在对不同径阶枯损率预估时,与实际枯损率之间的差异比较图。从图 3可以看出,枯损主要集中在2~10cm径阶,大于10cm径阶时枯损率逐渐下降并趋于稳定。截止到调查年份为止,该19年生杂种落叶松人工林的最大枯损径阶为20cm。两个模型的预估的枯损率计算结果显示,二者均在一定程度上低估了单木发生枯损的可能。从总体上看,混合模型计算的枯损率更加接近实际枯损率,综上所述,混合模型比基础Logistics模型具有更高的预测精度。

      图  3  广义线性模型、广义线性混合模型预测枯损率与实际枯损率比较

      Figure 3.  Comparison of the actual mortality rate in the result of generalized linear model and generalized linear mixed model predicted

    • 图 4列出了4个潜在阈值点ABCD,点A所在曲线为二分类事件总体的错误分类率随临界值的变化趋势线,点A表示错误分类率最小的点,此时的阈值为0.08。B点所在两条曲线分别是发生事件错误分类率以及未发生事件错误分类率随临界值变化趋势的交点,此时阈值为0.06。C点为0.5,是理论上枯损存活发生概率的临界值。D点为混合模型计算的错误分类率最小的点,此时的阈值为0.1。

      图  4  杂种落叶松人工幼龄林枯损模型阈值点与分类率关系

      Figure 4.  Relationship of classification rate and threshold for the mortality model of hybrid larch young plantations

      表 9为4个阈值点列联表分析具体结果:C=0.5时,灵敏度(TPR)为32.4%,特异度(TNR)为99.2%,此时表 10C点的预测存活树木分类正确率(ACR)为81.2%,仅有32.4%的枯损被正确分类,67.6%的枯损树木被预估成存活树木,阈值0.5明显不适用。阈值调整后,假阴性率(FNR)显著降低,二者相交于B点,此时阈值为0.06。由表 9可以看出:枯损正确分类率(TPR)明显提高,82.3%的枯损树木被正确的分类,存活树木的分类率与之近似相等,此时正确分类率(ACR)为82.1%,B点优于C点的预估判断。错误分类率(MCR)在阈值为0.08时达到最小,此阈值点下模型的特异度(TNR)为83.6%,灵敏度(TPR)为82.6%,正确分类率(ACR)为82.9%,均高于阈值0.06,说明在该阈值点下枯损树木和存活树木都能被更好地进行预估分类,确定错误分类率最低时的点A(0.08)作为基础模型最佳阈值点。同理,对混合模型进行阈值分析,求得错误分类率最低点D=0.1,此时表 10中错误分类率(MCR)为23.9%,正确分类率为88.9%,从中可以看出混合模型能够提高模型的预估精度。

      表 9  杂种落叶松人工幼龄林枯损模型阈值预测列联表分析

      Table 9.  Confusion matrix of the hybrid larch mortality model for four cut points

      列联表分析
      Confusion matrix
      实际变量分类Actual variety classification
      Positive 0=枯损Dead Negative 1=存活Living
      预测变量分类
      Classification of predicted variety
      Positive
      0=枯损Dead
      TPR 83.6%(A=0.08) FPR 17.4%(A=0.08)
      82.3%(B=0.06) 17.9%(B=0.06)
      32.4%(C=0.5) 0.8%(C=0.5)
      86.1%(D=0.1) 10%(D=0.1)
      Negative
      1=存活Living
      FNR 16.4%(A=0.08) TNR 82.6%(A=0.08)
      17.7%(B=0.06) 82.1%(B=0.06)
      67.6%(C=0.5) 99.2%(C=0.5)
      13.9%(D=0.1) 90.0%(D=0.1)
      注:ABCD为阈值点。Note:A, B, C, D are threshold points.

      表 10  杂种落叶松人工幼龄林枯损模型不同阈值预测分类率

      Table 10.  Classification of the hybrid larch young plantation mortality model for four cut points

      阈值Threshold ACR MCR
      A=0.08 82.9% 33.7%
      B=0.06 82.1% 35.7%
      C=0.5 81.2% 68.4%
      D=0.1 88.9% 23.9%
      注:ACR(正确分类率):正确预测枯损与存活株数占总株数比率; MCR(错误分类率):假阳性率与假阴性率之和。Notes:ACR(accurate classification rate): correctly predicting the ratio of mortality and living trees to total trees; MCR(misclassification rate): the sum of false positive rate and false negative rate.
    • 杂种落叶松广义线性模型最终由变量DBH、DBH2、BALD、BA组成,变量均通过χ2显著性检验,显著性水平α=0.05。Weiskittel等[8]建议在构建枯损模型时应当考虑单木、林分和竞争3个水平的变量,这与本文的研究结果相一致。其中胸径(DBH)是用来描绘树木活力的重要单木水平变量[12],Monserud等[6]利用澳大利亚地区森林资源清查数据研究单木枯损时,提出胸径及其变形对于多种针叶树种的枯损研究都是最为基础的研究变量,建议除单木胸径外,还需要考虑DBH2和DBH-1两个变量。本文选择加入DBH2变量的目的是为了对未来大径阶杂种落叶松枯损进行模拟,更好地描述树木生长的U型曲线。简单地说,树木在小径阶时由于竞争激烈,枯损发生的概率较高。随着树木的生长,死亡树木遭到淘汰,竞争逐渐减弱,枯损概率降低。当径阶增加到一定值时,树木进入老龄化,枯损概率又开始增加。很多学者使用DBH-1这一变量构造U型曲线[13, 32],本文在使用这一变量时,模型产生了严重的错误偏差。正如Groom等[25]在研究中指出的,当数据中包含较多小径阶树木时(例如,胸径 < 5.0cm),该变量会使得模型过高估计枯损概率发生的概率。DBH系数为负说明杂种落叶松枯损的概率随着树木胸径的增加而逐渐减小,径阶越小,对外界干扰的抵抗能力越小,发生枯损的概率越高。这是因为树木生理学以及环境压力导致的结果[14]。林分断面积(BA)是用来描述林分每公顷株数以及单木直径的林分水平变量,体现了林分内树木的拥挤程度,也体现了杂种落叶松在垂直方向对光照,水平方向对空间、水分、养分的需求,在多个单木枯损模型中出现[8, 10]。模型系数为正,表明与枯损概率之间呈正比,随着BA的增加,树木发生枯损的概率逐渐增加,Temesgen等[12]的研究表明林分断面积是枯损模型中高枯损率发生过程中最为重要的变量之一。早期研究常使用胸径和林分平均直径之比表述单木之间的竞争关系[18],随后Monserud[2]的研究表明用树冠竞争因子描述竞争关系更为有效,但是树冠竞争因子有着很大的变异性,且在实际林业外业中测量较少,所以未被推广使用。BALD被证明是近些年单木枯损模型中表现最好的单边竞争变量[32],本文中该变量系数为正,表明变量BALD与枯损概率之间呈正比关系,随着BALD的增加,对象木在林分中越处于劣势状态,受到的来自其他单木的竞争更强,单木枯损发生的概率逐渐增加,需要强调指出的是该变量不适合在已经间伐过的林分中使用[18]

      在杂种落叶松单木枯损模型的基础上加入随机参数,选择UN矩阵作为方差-协方差结构,构建样地水平广义线性混合模型。从图 2可以看出,混合模型有着更高的模型精度。图 3表明通过混合模型预估的林分不同径阶枯损率更接近实际枯损率。Timilsina等[26]通过混合模型对美国弗洛利达地区湿地松(Pinus elliottii)枯损模型进行研究,提出混合模型能够大大增加枯损的预估效果。混合模型之所以能够提高模型的预估精度,主要是因为它能够将林业中这种多变量多层次数据中一些潜在或未观测的记录,包括坡度、土壤特性、排水等因素包含在模型中,极大程度地考虑了样地与样地之间的差异性对模型的影响[25]。本文选用自适应积分最大似然估计对广义线性混合模型进行拟合,比限制伪似然估计和拉普拉斯估计能提供更加精准的函数逼近[35]。在林业混合模型构建中,除了考虑数据的多层次结构(如样地水平、单木水平)带来的异质性影响外,还需要考虑时间自相关性造成的模型偏差。本文考虑了利用自回归结构或移动平移结构对时间自相关进行处理,但是没有取得预计的效果。Timilsina等[26]解释了发生这种可能的原因:枯损是高度随机的过程,主要分为竞争枯损也称为自然枯损,以及突发枯损也称为非自然枯损两个种类,其中自然枯损受时间影响较小,而非自然枯损受时间的影响极大。在研究病虫害、火灾、风暴等气候因素造成的森林突发性枯损时,时间自相关处理能够显著提高模型的效果。本文研究对象多为竞争造成的枯损,受时间自相关影响程度较小。随着全球气候的不断变化,干旱、洪水等极端天气对森林造成的破坏程度与日俱增,结合气候因子对树木枯损进行分析,是未来枯损研究的一个发展方向。

      Chen等[33]对北美短叶松(Pinus banksiana)进行枯损研究的数据中存活树木的数量n=14797、枯损树木的数量n=1235、林分总枯损率=7.7%,此时最佳阈值点是灵敏度与特异度的交点,阈值为0.0855。Hein等[28]在研究枝条枯损时确定错误分类率最小值时的阈值点为最佳阈值点,枝条的存活率为31.7%,其最佳阈值为0.328。Crecente-campo等[34]确定选用的固定阈值方法更适宜西班牙地区辐射松(Pinups radiata)单木枯损模型,此时的固定阈值接近林分实际枯损发生率(0.71),而在另一枯损模型中又选择了灵敏度与特异度曲线交点作为最佳阈值,说明对于不同数据,阈值选择可能不同。结合上述研究和本文的阈值分析得出以下结论:(1)最佳阈值的确定是对一定区域内特定树种或林分的大数据平均处理,本文得出的错误分类率最小值时的最佳阈值是对于牡丹江地区杂种落叶松枯损而言,建模者在进行二分类数据处理时应该选择多种方法进行比较,选择最适宜研究数据的方法。(2)在研究枯损模型时不能为了追求较高的枯损判断正确率,而牺牲存活树木的判断正确率,否则会有很多实际存活的树木被错误地判断为死亡树木,反之亦然。假阳性率(FPR)为枯损树木被错误地判断为存活树木的比率,该结果会增加林分现存株数密度,从而影响后续经营方案设计。假阴性率(FPR)为存活树木被错误地判断为枯损树木的比率,该结果会降低林分蓄积,影响该树种生长模型的研究。(3)阈值与林分枯损率之间有着明显的关系,针对枯损事件发生记为1时,相对而言林分枯损率越大,得到的最佳阈值越大;林分枯损率越小得到的最佳阈值越小。阈值分析是枯损模型研究中非常重要的环节,提供了将枯损概率转换为实际枯损事件的判断值,有效地提高了模型的判别精度,在枯损模型建立过程中应对其进行详细的计算。

      杂种落叶松作为黑龙江省东南部地区大力推广的速生性树种,对其幼龄时期的枯损进行研究有助于准确预测杂种落叶松林的生长动态,为确定杂种落叶松合理的初值密度、抚育间伐时间、强度和间隔期以及人工整枝等经营措施提供依据。杂种落叶松在胸径、树高等生长特性方面较亲本具有明显优势,张含国等[36]通过在牡丹江林口县种子园的杂种落叶松研究显示,杂种落叶松在胸径生长方面超过亲本13.4%。杨书文等[37]通过对杂种落叶松进行研究表明,杂种落叶松在抗寒性、耐阴性、抗病性等方面具有显著提高。杂种落叶松作为强阳性树种,竞争是其枯损的主要影响因子,竞争包括对养分、水分的需求,这与其速生性有着紧密联系。幼龄时期生长较快对于养分和水分的需求很大,加剧了林木间竞争的发生,在早期营林过程中应采取一定措施加强杂种落叶松的养分供应。枯损发生为非瞬时的,一旦处于竞争劣势,杂种落叶松会处于一个停止生长的状态,随着时间推移而逐渐枯死,在森林经营时可通过抚育措施控制树木间距,确保单木存活空间,从而减少枯损发生的可能。林分密度和树种生长速度会加剧林分内树木之间的竞争,加快树木枯损,这主要体现在自然枯损方面,自然枯损能够通过模型形式体现。而对于树种抗寒性、耐阴性、抗病性主要对非自然枯损造成影响,杂种落叶松在这些特性上均相对于亲本有着显著的提高,理论上杂种落叶松发生枯损的可能性应该较亲本低,具体的模型构建属于枯损机理模型,目前由于缺乏相应数据,无法进行研究,以后研究时可将此作为重点。

    • 杂种落叶松人工幼龄林单木枯损模型研究表明:与传统的广义线性模型(GLM)相比,广义线性混合模型(GLMM)包含了更多的信息,能够更有效地描述杂种落叶松的单木枯损。阈值分析能够同时判断树木枯损和存活,是提高二分类模型预测准确性的有效方法。杂种落叶松作为速生性极强的树种,为了更好地推广应用,降低单木枯损的发生概率,幼龄时期应及时开展抚育间伐、人工整枝等经营措施调整林木的竞争状态,确保存活木的生长空间。

参考文献 (37)

目录

    /

    返回文章
    返回