高级检索

留言板

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

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

基于零膨胀模型及混合效应模型相结合的蒙古栎林林木进界模拟研究

李春明 李利学

引用本文:
Citation:

基于零膨胀模型及混合效应模型相结合的蒙古栎林林木进界模拟研究

    作者简介: 李春明,博士,副研究员。主要研究方向:森林生长模型。Email:lichunm@ifrit.ac.cn 地址:100091中国林业科学研究院资源信息研究所.
  • 中图分类号: S758.5

Simulating study on tree recruitment of Quercus mongolica based on zero-inflated model and mixed effect model methods

  • 摘要: 目的 林木的进界是确保森林长期维持的基本条件,而进界模型能够预测森林的发展,是量化森林生态系统未来健康和生产力的基础。方法 以吉林省1995年设立的295块蒙古栎固定样地数据为例,构建基于林分因子、立地因子及气象因子的蒙古栎林林木进界模型。模型的基本形式包括泊松分布和负二项分布两种离散形式。考虑到样地中存在大量零值的问题,在这些基础模型上考虑加入零膨胀模型。为了解决模型存在的嵌套和纵向数据问题,在构建模型时把样地的随机效应考虑进去。最后利用验证数据来验证。结果 林分算数平均直径和林分公顷株数是影响林木进界概率和数量最重要的影响因子,并且均与林木进界概率和数量呈反比。立地和气象因子中的各项因子对进界均没有产生明显影响。负二项分布模型由于考虑了数据过度离散问题,模拟精度要高于泊松分布;在考虑样地的随机效应后,除了标准负二项分布模型外所有模型都明显提高了模型的模拟精度;同时考虑随机效应和零膨胀的负二项分布模型,其模型的模拟效果最好,验证结果也支持此结论。结论 为了确保进界的发生,在进行森林经营时,确定合理的初植和经营密度至关重要。
  • 表 1  蒙古栎样地各因子统计表

    Table 1.  Statistical table of factors in Quercus mongolica sample plots

    影响进界的因子
    Factor affecting recruitment
    调查年份
    Survey year
    指标
    Index
    平均值 (标准差)
    Mean (standard deviation)
    最大值
    Max.
    最小值
    Min.
    林分和单木因子
    Stand and single wood factor
    1999 胸径 DBH (d)/cm 12.7 (7.5) 82.7 5
    林分平均直径 Average stand diameter (Dg)/cm 15.3 (4.2) 30.9 6.3
    公顷断面积/(m2·hm− 2) Basal area/(m2·ha− 1) 22.9 (9.4) 57.7 3.0
    公顷株数/(株·hm− 2) Stand density/(tree·ha− 1) 1 343 (621) 3 317 250
    林分和单木因子
    Stand and single wood factor
    2004 胸径 DBH (d)/cm 12.9 (7.6) 83.6 5
    林分平均直径 Average stand diameter (Dg)/cm 15.7 (4.1) 31.3 6.3
    公顷断面积/(m2·hm− 2) Basal area/(m2·ha− 1) 24.1 (8.8) 59.7 3.8
    公顷株数/(株·hm− 2) Stand density/(tree·ha− 1) 1 370 (618) 3 250 350
    林分和单木因子
    Stand and single wood factor
    2009 胸径 DBH (d)/cm 13.5 (8.0) 84.9 5
    林分平均直径 Average stand diameter (Dg)/cm 16.5 (4.1) 32.6 6.6
    公顷断面积/(m2·hm− 2) Basal area/(m2·ha− 1) 26.2 (8.5) 61.2 7.3
    公顷株数/(株·hm− 2) Stand density/(tree·ha− 1) 1 347 (584) 3 550 367
    立地因子
    Site factor
    海拔 Elevation/m 596 (196) 1 280 100
    坡度 Slope degree/(°) 21 (8.5) 45 0
    坡向 Slope aspect 按方位角从0º ~ 360º,共分成9个坡向,每个坡向大概45º,用数字1 ~ 9表示
    According to the azimuth angle from 0−360 degrees, it is divided into 9 slope directions, each of which is about 45 degrees, represented by No. 1−9
    下载: 导出CSV

    表 2  主要气象变量统计表

    Table 2.  Satistical table of main climate variables

    气象因子
    Climate factor
    最大值
    Max.
    最小值
    Min.
    平均值 (标准差)
    Mean (standard deviation)
    年份
    Year
    年平均温度
    Mean annual temperature (MAT)/℃
    6.56 1.8 4.2 (0.8) 1999—2004
    6.96 2.3 4.6 (0.8) 2004—2009
    6 1.3 3.7 (0.8) 2009—2014
    最暖月平均气温
    Mean temperature of the warmest month (MWMT)/℃
    23.58 18.6 21.1 (0.9) 1999—2004
    22.68 18.4 20.8 (0.8) 2004—2009
    22.78 18.3 20.8 (0.8) 2009—2014
    最冷月平均气温
    Mean temperature of the coldest month (MCMT)/℃
    − 11.5 − 18.4 − 15.7 (1.1) 1999—2004
    − 10.1 − 16.9 − 14.4 (1.1) 2004—2009
    − 12.1 − 19 − 16.5 (1.1) 2009—2014
    年平均降水量
    Mean annual precipitation (MAP)/mm
    920.6 498.6 646.3 (75.3) 1999—2004
    1 038.8 475.2 640.9 (113.3) 2004—2009
    1 139.2 518.2 705.3 (126.4) 2009—2014
    年平均夏季 (5月至9月)降水量
    Mean annual summer (from May to September) precipitation (MSP)/mm
    731.4 399.2 514.9 (53.7) 1999—2004
    827.6 379.2 512.9 (88.9) 2004—2009
    884.2 389.8 539.4 (99.0) 2009—2014
    无霜期天数
    Number of frost-free days (NFFD)
    196.8 156.6 176.6 (7.9) 1999—2004
    200.2 158.6 178.9 (8.1) 2004—2009
    194.6 154.8 174.0 (7.7) 2009—2014
    上一年8月至当年7月的降雪量
    Snowfall between August in previous year and July in current year (PAS)/ mm
    115.8 35.2 54.9 (16.3) 1999—2004
    123.2 35 60.9 (14.8) 2004—2009
    154.6 52.4 80.5 (19.7) 2009—2014
    下载: 导出CSV

    表 3  林分进界模型模拟结果

    Table 3.  Simulation results of stand-level recruitment models

    参数
    Parameter
    基础模型 Basic model
    M1M2M3M4
    不考虑随机效应
    No random effect
    模型的零部分
    Zero component of the model
    $ {\alpha _0}$ 5.730 8 (0.990 3)*** 5.726 6 (1.000 1)***
    $ {\alpha _1}$ − 0.203 8 (0.046 4)*** − 0.202 4 (0.046 8)***
    $ {\alpha _2}$ − 1.110 0 (0.265 1)*** − 1.110 3 (0.267 7)***
    计数的部分
    Positive count component of the model
    $\; {\beta _0}$ 9.713 9 (0.026 7)*** 9.237 8 (0.027 6)*** 8.446 8 (0.469 7)*** 8.065 1 (0.295 0)***
    $\; {\beta _1}$ − 0.313 0 (0.001 7)*** − 0.274 4 (0.001 7)*** − 0.214 1 (0.021 4)*** − 0.187 4 (0.013 8)***
    $\; {\beta _2}$ − 0.756 0 (0.008 3)*** − 0.621 9 (0.008 5)*** − 0.771 5 (0.149 9)*** − 0.583 3 (0.095 2)***
    $ k$ 2.437 2 0.753 2
    AIC 55 314 39 871 4 877.6 4 671.8
    BIC 55 327 39 896 4 894.2 4 700.9
    − 2logL 55 308 39 859 4 869.6 4 657.8
    参数
    Parameter
    考虑混合效应模型 With mixed effect model
    M5 M6 M7 M8
    考虑随机效应部分
    With random effect
    模型的零部分
    Zero component of the model
    $ {\alpha _0}$ 5.729 5 (0.990 3)*** 5.716 6 (0.993 3)***
    $ {\alpha _1}$ − 0.203 8 (0.046 4)*** − 0.202 6 (0.046 5)***
    $ {\alpha _2}$ − 1.109 8 (0.265 2)*** − 1.108 0 (0.265 9)***
    计数的部分
    Positive count component of the model
    $\; {\beta _0}$ 10.453 3 (0.150 7)*** 9.587 2 (0.113 4)*** 8.995 2 (0.552 7)*** 8.263 9 (0.354 6)***
    $\; {\beta _1}$ − 0.391 3 (0.006 4)*** − 0.259 9 (0.006 3)*** − 0.240 6 (0.026 7)*** − 0.202 0 (0.017 8)***
    $\; {\beta _2}$ − 0.992 8 (0.017 5)*** − 1.220 0 (0.020 8)*** − 0.926 2 (0.169 5)*** − 0.687 1 (0.107 8)***
    $ k$ 2.292 2 0.525 7
    AIC 21 504 15 054 4 882.4 4 645.0
    BIC 21 518 15 078 4 899.8 4 672.7
    − 2logL 21 496 15 040 4 872.4 4 629.0
    随机效应方差协方差矩阵 Covariance matrix of random effect variance 2.976 0 0.673 0 0.000 1 0.254 7
    注:k是负二项分布的待估参数;M1为标准泊松分布进界模型,即公式 (1);M2为零膨胀泊松分布进界模型,即公式 (2);M3为标准负二项分布进界模型,即公式 (3);M4为零膨胀负二项分布进界模型,即公式 (4);M5为考虑随机效应的标准泊松分布进界模型;M6为考虑随机效应的零膨胀泊松分布进界模型;M7为考虑随机效应的标准负二项分布进界模型;M8为考虑随机效应的零膨胀负二项分布进界模型。表中括号内的值为标准差;*为0.01 < P < 0.05,**为0.001 < P < 0.01,***P < 0.001;$ {\alpha _0}$$\; {\beta _0}$为截距参数,$ {\alpha _1}$$\; {\beta _1}$为林分算数平均直径参数值,$ {\alpha _2}$$\; {\beta _2}$为林分公顷株数参数值。Notes: k is a parameter to be estimated in negative binomial distribution; M1 is Poisson distribution recruitment model, i.e. formula (1); M2 is zero-inflated Poisson distribution recruitment model, i.e. formula (2); M3 is negative binomial distribution recruitment model, i.e. formula (3); M4 is zero-inflated negative binomial distribution recruitment model, i.e. formula (4); M5 is Poisson distribution recruitment model based on mixed effect model; M6 is zero-inflated Poisson distribution recruitment model based on mixed effect model; M7 is negative binomial distribution recruitment model based on mixed effect model; M8 is zero-inflated negative binomial distribution recruitment model based on mixed effect model. Values in brackets are standard deviation; * is 0.01 < P < 0.05,** is 0.001 < P < 0.01,*** is P < 0.001; $ {\alpha _0}$ and $\; {\beta _0}$ are intercept parameters, $ {\alpha _1}$ and $\; {\beta _1}$ are stand arithmetic average diameter parameters, $ {\alpha _2}$ and $\; {\beta _2}$ are stand number of per hectare parameters.
    下载: 导出CSV

    表 4  利用LRT指标对模型进行比较的结果

    Table 4.  Comparing results of the model based on LRT index

    模型 ModelLRT值 LRT valueP模型 ModelLRT 值 LRT valueP模型 ModelLRT值 LRT valueP
    M1/M2 15 449 < 0.000 1 M1/M5 33 812 < 0.000 1 M3/M4 211.8 < 0.000 1
    M2/M6 24 819 < 0.000 1 M3/M7 2.8 > 0.05 M4/M8 28.8 < 0.000 1
    下载: 导出CSV

    表 5  模型验证结果

    Table 5.  Validation results of models

    评价指标
    Evaluating index
    不考虑随机效应 No random effect考虑随机效应 With random effect
    M1M2M3M4M5M6M7M8
    $ {R^2}$ 0.40 0.55 0.50 0.65 0.75 0.91 0.54 0.84
    RMSE/株 RMSE/tree
    163.8 146.2 158.7 139.2 125.3 161.8 153.8 116.5
    $ |\bar E|$/株 $ |\bar E|$/tree 98.8 95.9 83.9 81.7 65 70.2 85 61.7
    下载: 导出CSV
  • [1] Kangas A, Maltamo M. Forest inventory[M]. Dordrecht: Springer, 2006.
    [2] Li R, Weiskittel A R, Kershaw J A, Jr. Modeling annualized occurrence, frequency, and composition of ingrowth using mixed-effects zero-inflated models and permanent plots in the Acadian forest region of North America[J]. Canadian Journal of Forest Research, 2011, 41: 2077−2089. doi: 10.1139/x11-117
    [3] Gauthier M M, Guillemette F, Bedard S. On the relationship between saplings and ingrowth in northern hardwood stands[J]. Forest Ecology and Management, 2015, 358: 261−271. doi: 10.1016/j.foreco.2015.09.020
    [4] Russell M B, Westfall J A, Woodall C W. Modeling browse impacts on sapling and tree recruitment across forests in the northern United States[J]. Canadian Journal of Forest Research, 2017, 47: 1474−1481. doi: 10.1139/cjfr-2017-0155
    [5] Shifley S R, Ek A R, Burk T E. A generalized methodology for estimating forest ingrowth at multiple threshold diameters[J]. Forest Science, 1993, 39: 776−798.
    [6] Zhang X Q, Lei Y C, Cai D X, et al. Predicting tree recruitment with negative binomial mixture models[J]. Forest Ecology and Management, 2012, 270: 209−215. doi: 10.1016/j.foreco.2012.01.028
    [7] Weiskittel A R, Hann D W, Kershaw J A, et al. Forest growth and yield modeling[M]. New Jersey: Wiley-Blackwell, 2011.
    [8] Vanclay J K. Modelling regeneration and recruitment in a tropical rain forest[J]. Canadian Journal of Forest Research, 1992, 22(9): 1235−1248. doi: 10.1139/x92-165
    [9] Trasobares A, Pukkala T, Miina J. Growth and yield model for uneven-aged mixtures of Pinus sylvestris L. and Pinus nigra Arn. in Catalonia, northeast Spain[J]. Annual of Forest Science, 2004, 61: 9−24. doi: 10.1051/forest:2003080
    [10] Liang J, Buongiorno J, Monserud R A. Growth and yield of all-aged Douglas fir western hemlock forest stands: a matrix model with stand diversity effects[J]. Canadian Journal of Forest Research, 2005, 35: 2368−2381. doi: 10.1139/x05-137
    [11] Hao Q, Meng F, Zhou Y, et al. A transition matrix growth model for uneven-aged mixed-species forests in the Changbai Mountains, northeastern China[J]. New Forest, 2005, 29: 221−231. doi: 10.1007/s11056-005-5657-z
    [12] Xiang W, Lei X D, Zhang X Q. Modelling tree recruitment in relation to climate and competition in semi-natural Larix-Picea-Abies forests in northeast China[J]. Forest Ecology and Management, 2016, 382: 100−109. doi: 10.1016/j.foreco.2016.09.050
    [13] Adame P, del Rìo M, Cañellas I. Ingrowth model for pyrenean oak stands in north-western Spain using continuous forest inventory data[J]. European Journal of Forest Research, 2010, 129: 669−678. doi: 10.1007/s10342-010-0368-1
    [14] Lexerød N L. Recruitment models for different tree species in Norway[J]. Forest Ecology and Management, 2005, 206: 91−108. doi: 10.1016/j.foreco.2004.11.001
    [15] Bravo F, Pando V, Ordóñez C, et al. Modelling ingrowth in mediterranean pine forests: a case study from scots pine (Pinus sylvestris L.) and mediterranean maritime pine (Pinus pinaster Ait.) stands in Spain[J]. Forest System, 2008, 17: 250−260.
    [16] Klopcic M, Poljanec A, Boncina A. Modelling natural recruitment of European beech (Fagus sylvatica L.)[J]. Forest Ecology and Management, 2012, 284: 142−151. doi: 10.1016/j.foreco.2012.07.049
    [17] Yang Y, Huang S. Two-stage ingrowth models for four major tree species in Alberta[J]. European Journal of Forest Research, 2015, 134: 991−1004. doi: 10.1007/s10342-015-0904-0
    [18] Lambert D. Zero-inflated Poisson regression, with an application to defects in manufacturing[J]. Technometrics, 1992, 34: 1−14. doi: 10.2307/1269547
    [19] Hall D B. Zero-inflated Poisson and binomial regression with random effects: a case study[J]. Biometrics, 2000, 56: 1030−1039. doi: 10.1111/j.0006-341X.2000.01030.x
    [20] Ripley T, Scrimgeour G, Boyce M S. Bull trout (Salvelinus confluentus) occurrence and abundance influenced by cumulative industrial developments in a Canadian boreal forest watershed[J]. Canadian Journal of Fisheries and Aquatic Sciences, 2005, 62: 2431−2442. doi: 10.1139/f05-150
    [21] Cunningham R B, Lindenmayer D B. Modeling count data of rare species: some statistical issues[J]. Ecology, 2005, 86: 1135−1142. doi: 10.1890/04-0589
    [22] Fortin M, DeBlois J. Modeling tree recruitment with zero-inflated models: the example of hardwood stands in southern Québec, Canada[J]. Forest Science, 2007, 53(4): 529−539.
    [23] Zell J, Rohner B, Thürig E, et al. Modeling ingrowth for empirical forest prediction systems[J]. Forest Ecology and Management, 2019, 433: 771−779. doi: 10.1016/j.foreco.2018.11.052
    [24] Fang Z, Bailey R L. Nonlinear mixed effects modeling for slash pine dominant height growth following intensive silvicultural treatments[J]. Forest Science, 2001, 47: 287−300.
    [25] Wang Y, Lemay V M, Baker T G. Modeling and prediction of dominant height and site index of Eucalyptus globules plantations using a nonlinear mixed-effects model approach[J]. Canadian Journal of Forest Research, 2007, 37(8): 1390−1403. doi: 10.1139/X06-282
    [26] Zhao D, Wilson M, Borders B E. Modeling response curves and testing treatment effects in repeated measures experiments: a multilevel nonlinear mixed-effects models approach[J]. Canadian Journal of Forest Research, 2005, 35(1): 122−132. doi: 10.1139/x04-163
    [27] Yoshida T, Noguchi M, Akibayashi Y, et al. Twenty years of community dynamics in a mixed conifer-broad-leaved forest under a selection system in northern Japan[J]. Canadian Journal of Forest Research, 2006, 36: 1363−1375. doi: 10.1139/x06-041
    [28] Alves L F, Vieira S A, Scaranello M A, et al. Forest structure and live aboveground biomass variation along an elevational gradient of tropical Atlantic moist forest (Brazil)[J]. Forest Ecology and Management, 2010, 260: 679−691. doi: 10.1016/j.foreco.2010.05.023
    [29] Batllori E, Camarero J J, Ninot J M, et al. Seedling recruitment, survival and facilitation in alpine Pinus uncinata tree line ecotones. Implications and potential responses to climate warming[J]. Global Ecology and Biogeography, 2009, 18(4): 460−472. doi: 10.1111/j.1466-8238.2009.00464.x
    [30] He Z B, Zhao W Z, Zhang L J, et al. Response of tree recruitment to climatic variability in the alpine treeline ecotone of the Qilian Mountains, northwestern China[J]. Forest Science, 2013, 59(1): 118−126. doi: 10.5849/forsci.11-044
    [31] Wang T, Wang G, Innes J L, et al. Climate AP: an application for dynamic local downscaling of historical and future climate data in Asia Pacific[J]. Frontiers of Agricultural Science and Engineering, 2017, 4(4): 448−458. doi: 10.15302/J-FASE-2017172
    [32] Zuur A F, Ieno E N, Walker N J, et al. Mixed effects models and extensions in ecology with R[M/OL]. [2019−01−10]. https://link.springer.com/content/pdf/10.1007%2F978-0-387-87458-6.pdf.
    [33] Akaike H. A new look at the statistical model identification[J]. IEEE Transactions on Automatic Control, 1974, 19(6): 716−723. doi: 10.1109/TAC.1974.1100705
    [34] Vonesh E F, Chinchilli V M. Linear and nonlinear models for the analysis of repeated measurements[M]. New York: Marcel Dekker, 1997.
    [35] Calama R, Montero G. Interregional nonlinear height-diameter model with random coefficients for stone pine in Spain[J]. Canadian Journal of Forest Research, 2004, 34: 150−163. doi: 10.1139/x03-199
    [36] Bose A K, Harvey B D, Brais S. Sapling recruitment and mortality dynamics following partial harvesting in aspen-dominated mixedwoods in eastern Canada[J]. Forest Ecology and Management, 2014, 329: 37−48. doi: 10.1016/j.foreco.2014.06.004
    [37] De Avila A L, Schwartz G, Ruschel A R, et al. Recruitment, growth and recovery of commercial tree species over 30 years following logging and thinning in a tropical rain forest[J]. Forest Ecology and Management, 2017, 385: 225−235. doi: 10.1016/j.foreco.2016.11.039
  • [1] 张琴林天喜王贵春孙国文范秀华 . 红松、蒙古栎和色木槭凋落物混合分解研究. 北京林业大学学报, 2014, 36(6): 106-111. doi: 10.13332/j.cnki.jbfu.2014.06.020
    [2] 张怡卓苏耀文李超门洪生 . 蒙古栎木材MOR与MOE的近红外光谱预测模型分析. 北京林业大学学报, 2016, 38(8): 99-105. doi: 10.13332/j.1000-1522.20150505
    [3] 张杰李健康段安安孙宇涵刘洪山武素然李贵芬李云 . 不同质量浓度NAA、IBA对栓皮栎、蒙古栎黄化嫩枝扦插生根的影响. 北京林业大学学报, 2019, 41(7): 128-138. doi: 10.13332/j.1000-1522.20190103
    [4] 张晓红张会儒 . 蒙古栎次生林垂直结构特征对目标树经营的响应. 北京林业大学学报, 2019, 41(5): 56-65. doi: 10.13332/j.1000-1522.20190046
    [5] 李艳丽杨华邓华锋 . 蒙古栎−糠椴天然混交林空间格局研究. 北京林业大学学报, 2019, 41(3): 33-41. doi: 10.13332/j.1000-1552.20180236
    [6] 樊莹乔雪涛赵秀海 . 长白山自然保护区蒙古栎幼树生理生长特性随海拔梯度的变化. 北京林业大学学报, 2019, 41(11): 1-10. doi: 10.13332/j.1000-1522.20190095
    [7] 郭斌 . 栎属近缘种指纹图谱构建及遗传结构. 北京林业大学学报, 2018, 40(5): 10-18. doi: 10.13332/j.1000-1522.20170444
    [8] 陈彬杭温晓示张树斌柴世品孙晗王襄平 . 吉林北部山区长白落叶松林径向生长对气候干暖化的响应. 北京林业大学学报, 2018, 40(12): 18-26. doi: 10.13332/j.1000-1522.20180333
    [9] 欧强新雷相东沈琛琛宋国涛 . 基于随机森林算法的落叶松−云冷杉混交林单木胸径生长预测. 北京林业大学学报, 2019, 41(9): 9-19. doi: 10.13332/j.1000-1522.20180266
    [10] 周超凡张会儒徐奇刚雷相东 . 基于相邻木关系的林层间结构解析. 北京林业大学学报, 2019, 41(5): 66-75. doi: 10.13332/j.1000-1522.20190051
    [11] 欧光龙胥辉王俊峰肖义发陈科屹郑海妹 . 思茅松天然林林分生物量混合效应模型构建. 北京林业大学学报, 2015, 37(3): 101-110. doi: 10.13332/j.1000-1522.20140316
    [12] 姚丹丹徐奇刚闫晓旺李玉堂 . 基于贝叶斯方法的蒙古栎林单木枯死模型. 北京林业大学学报, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260
    [13] 陈国栋丁佩燕尹忠东 . 基于混合效应模型的新疆天山云杉单木胸径预测模型构建. 北京林业大学学报, 2020, 35(): 1-11. doi: 10.12171/j.1000-1522.20190236
    [14] 臧颢雷相东张会儒李春明卢军 . 红松树高-胸径的非线性混合效应模型研究. 北京林业大学学报, 2016, 38(6): 8-9. doi: 10.13332/j.1000-1522.20160008
    [15] 王涛董利虎李凤日 . 基于混合效应的杂种落叶松人工幼龄林单木枯损模型. 北京林业大学学报, 2018, 40(10): 1-10. doi: 10.13332/j.1000-1522.20170437
    [16] 王曼霖董利虎李凤日 . 基于Possion回归混合效应模型的长白落叶松一级枝数量模拟. 北京林业大学学报, 2017, 39(11): 45-55. doi: 10.13332/j.1000-1522.20170204
    [17] 李春明 . 利用非线性混合模型进行杉木林分断面积生长模拟研究. 北京林业大学学报, 2009, 31(1): 44-49.
    [18] 王树力周健平 , . 基于结构方程模型的林分生长与影响因子耦合关系分析. 北京林业大学学报, 2014, 36(5): 7-12. doi: 10.13332/j.cnki.jbfu.2014.05.011
    [19] 颜伟段光爽王一涵孙钊周桃龙符利勇 . 河南省栎类和杨树林分断面积和蓄积生长模型构建. 北京林业大学学报, 2019, 41(6): 55-61. doi: 10.13332/j.1000-1522.20180311
    [20] 李杨亢新刚 . 长白山云冷杉针阔混交林林木空间利用率混合模型. 北京林业大学学报, 2020, 42(5): 71-79. doi: 10.12171/j.1000-1522.20190112
  • 加载中
表(5)
计量
  • 文章访问数:  456
  • HTML全文浏览量:  90
  • PDF下载量:  28
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-05-09
  • 录用日期:  2019-09-05
  • 网络出版日期:  2020-05-27
  • 刊出日期:  2020-06-01

基于零膨胀模型及混合效应模型相结合的蒙古栎林林木进界模拟研究

    作者简介: 李春明,博士,副研究员。主要研究方向:森林生长模型。Email:lichunm@ifrit.ac.cn 地址:100091中国林业科学研究院资源信息研究所
  • 1. 中国林业科学研究院资源信息研究所,北京 100091
  • 2. 河北省承德县五道河林场,河北 承德 067407

摘要: 目的林木的进界是确保森林长期维持的基本条件,而进界模型能够预测森林的发展,是量化森林生态系统未来健康和生产力的基础。方法以吉林省1995年设立的295块蒙古栎固定样地数据为例,构建基于林分因子、立地因子及气象因子的蒙古栎林林木进界模型。模型的基本形式包括泊松分布和负二项分布两种离散形式。考虑到样地中存在大量零值的问题,在这些基础模型上考虑加入零膨胀模型。为了解决模型存在的嵌套和纵向数据问题,在构建模型时把样地的随机效应考虑进去。最后利用验证数据来验证。结果林分算数平均直径和林分公顷株数是影响林木进界概率和数量最重要的影响因子,并且均与林木进界概率和数量呈反比。立地和气象因子中的各项因子对进界均没有产生明显影响。负二项分布模型由于考虑了数据过度离散问题,模拟精度要高于泊松分布;在考虑样地的随机效应后,除了标准负二项分布模型外所有模型都明显提高了模型的模拟精度;同时考虑随机效应和零膨胀的负二项分布模型,其模型的模拟效果最好,验证结果也支持此结论。结论为了确保进界的发生,在进行森林经营时,确定合理的初植和经营密度至关重要。

English Abstract

  • 进界被定义为在某一段时间内,抽样样地上树木的生长(通常是指树木的树高和胸径)达到规定的临界尺寸[1-2]。进界是确保森林长期维持的基本条件[3],而进界模型作为森林生长及收获预估系统中重要的组成部分,能够预测森林的发展,是量化森林生态系统未来健康和生产力的基础[4]。作为森林发展过程中重要的组成部分,进界并没有像生存生长、枯损和采伐那样得到应有的重视[5]。忽略进界将会导致森林生长及收获预估模型产生大的偏差[6]

    目前模拟林分的进界主要包括静态的和动态的两种方法[7]。静态的进界模型不考虑林分特征,假定进界的平均量在某一段时间内是不变的,这种方法通常被用到林分收获表里[8]。然而,当初始林分条件变化非常大或者准确的生长过程表缺乏时,这个简单假设通常是不满足要求的。动态的进界模型通过统计回归模型,把立地和林分条件作为解释变量,来预测未来林分的进界,能够提供相对准确的预测。最开始的动态模拟方法是采用线性和非线性进界模型来预测某一时期内进界数量,基于最小二乘方法[5, 9-11]。然而,这些模型总是预测进界会发生,事实上并不完全是这样,进界数据中经常会包括大量的零值问题[2, 12]。后来基于logistic模型的两步法被推荐来解决此问题,并被认为是一个很好的方法[8, 13]。在两步法模型中,第1个方程基于一系列协变量,来判断林分是否发生进界(假设因变量服从Bernoulli分布),第2个方程是在给定进界发生的情况下,基于相同或不同的协变量,建立进界株数与协变量的线性或非线性方程来估计进界的数量[14-17]。上述这些研究都把进界株数看成服从正态分布,但林分水平进界株数是典型的计数数据,是离散随机变量,服从二项分布、泊松分布或负二项分布,不服从正态分布。如果按正态分布就会造成大的误差。另外,在进行大量野外调查时,通常包括大量的样地数据,而这些样地通常会出现过量的零值,即大量样地没有进界情况。零膨胀模型(zero-inflated model)包括两个分布的联合概率,可解释存在过量零值的问题,并把主要影响因子作为协变量引入到模型中去[18]。很多领域都用到了零膨胀模型,例如制造业、园艺研究、野生动物监测等[19-21]。与两步法相类似,零膨胀模型的第1步是判断样地内进界发生与否,采用服从Bernoulli分布的二分量logistic模型来模拟。第2步是模拟协变量影响事件的丰度问题,选择零膨胀泊松分布或零膨胀负二项分布来模拟进界株数。Fortin等[22]最早选择零膨胀泊松模型和零膨胀离散Weibull分布来模拟硬木林的进界。由于泊松分布要求均值等于方差,而Weibull分布很难解释可加参数。零膨胀负二项分布模型逐渐被发现能够比零膨胀泊松模型模拟结果更好[2, 4, 6, 12]

    在时空尺度上,林分进界是一个复杂的过程,具有很高的变化性、复杂性和随机性,很难准确预测。受林分特征、气候和地理特征等多个因子及其交互作用的影响,估计这些变量的交互作用对进界的影响是一个挑战。在构建进界模型时,很多数据都是基于永久性固定样地的多次定期观测,每个林分包括多个样地,这些数据存在着多层嵌套结构以及多次观测间的自相关性及异方差,存在着数据和统计上的挑战。通过使用混合效应模型来模拟进界的发生和数量,可以成功地量化进界的随机性[2, 23]。虽然混合效应模型方法已经被广泛应用到了森林生长模拟研究中[24-26],但是在模拟进界中考虑随机效应的研究并不多见。只有Li等[2]和Zell等[23]考虑了样地的随机效应,并且表明与不考虑随机效应相比,模型精度明显提高,而其他几位学者考虑随机效应后模型不能够收敛[4, 6, 12]

    蒙古栎(Quercus mongolica)是中国东北林区中主要的次生林树种,多成大面积纯林或与红松(Pinus koraiensis)、云杉(Picea jezoensis)、冷杉(Abies nephrolepis)、水曲柳(Fraxinus mandshurica)及白桦(Betula platyphylla)等树种混生。以吉林省多次观测的蒙古栎林固定样地数据为例,把数据分成两部分,其中1999年和2009年两期共236块样地为模拟数据(80%),2004年选择了59块样地为验证数据(20%)。基于模拟数据,选择林分因子、立地因子和气象因子为协变量,采用标准泊松分布模型(Poisson模型)、标准负二项分布模型(NB)、零膨胀泊松分布模型(ZIP)和零膨胀负二项分布模型(ZINB)构建林分进界模型,并在这些模型的基础上考虑样地的随机效应,利用AIC、BIC和− 2loglikelihood(− 2logL)等模型拟合效果指标找出模拟精度最高的模型,并选择验证数据进行精度验证。

    • 本研究所用数据来源于吉林省295块从1999—2014年的蒙古栎林固定样地调查数据,样地大小为0.06 hm2,每5年调查一次,共选择了3次作为研究对象。在模型构建过程中,通常把研究数据分成模拟数据(一般占80%左右)和验证数据(一般占20%)两部分,1999年和2009年两期共236块样地数据为模拟数据(80%),2004年选择的59块样地为验证数据(20%)。样地中蒙古栎树种占7成以上,其他树种包括白桦、红松、水曲柳、黄菠萝(Phellodendron amurense)、云杉和冷杉等。调查的样地因子主要包括海拔、坡度、坡向、坡位、土壤厚度、优势树种、起源、平均林龄、平均树高、郁闭度等因子。调查的样木因子包括对胸径大于5 cm的树木每木胸径进行检尺,对上期胸径没有达到5 cm,而本期达到5 cm的树木进行挂牌并统计每个样地的株数。外业调查结束后进行每公顷进界株数、每公顷株数、每公顷断面积等指标的计算。具体样地因子见表1

      表 1  蒙古栎样地各因子统计表

      Table 1.  Statistical table of factors in Quercus mongolica sample plots

      影响进界的因子
      Factor affecting recruitment
      调查年份
      Survey year
      指标
      Index
      平均值 (标准差)
      Mean (standard deviation)
      最大值
      Max.
      最小值
      Min.
      林分和单木因子
      Stand and single wood factor
      1999 胸径 DBH (d)/cm 12.7 (7.5) 82.7 5
      林分平均直径 Average stand diameter (Dg)/cm 15.3 (4.2) 30.9 6.3
      公顷断面积/(m2·hm− 2) Basal area/(m2·ha− 1) 22.9 (9.4) 57.7 3.0
      公顷株数/(株·hm− 2) Stand density/(tree·ha− 1) 1 343 (621) 3 317 250
      林分和单木因子
      Stand and single wood factor
      2004 胸径 DBH (d)/cm 12.9 (7.6) 83.6 5
      林分平均直径 Average stand diameter (Dg)/cm 15.7 (4.1) 31.3 6.3
      公顷断面积/(m2·hm− 2) Basal area/(m2·ha− 1) 24.1 (8.8) 59.7 3.8
      公顷株数/(株·hm− 2) Stand density/(tree·ha− 1) 1 370 (618) 3 250 350
      林分和单木因子
      Stand and single wood factor
      2009 胸径 DBH (d)/cm 13.5 (8.0) 84.9 5
      林分平均直径 Average stand diameter (Dg)/cm 16.5 (4.1) 32.6 6.6
      公顷断面积/(m2·hm− 2) Basal area/(m2·ha− 1) 26.2 (8.5) 61.2 7.3
      公顷株数/(株·hm− 2) Stand density/(tree·ha− 1) 1 347 (584) 3 550 367
      立地因子
      Site factor
      海拔 Elevation/m 596 (196) 1 280 100
      坡度 Slope degree/(°) 21 (8.5) 45 0
      坡向 Slope aspect 按方位角从0º ~ 360º,共分成9个坡向,每个坡向大概45º,用数字1 ~ 9表示
      According to the azimuth angle from 0−360 degrees, it is divided into 9 slope directions, each of which is about 45 degrees, represented by No. 1−9
    • 在构建进界模型时,首先考虑的就是林分因子的影响。作为林分密度指标,林分公顷断面积[2, 4, 6, 22-23]、林分公顷株数[2, 16]、林分算数平均直径[6]以及林分平方平均直径[16]通常情况下会被选择一个或其中几个的影响。除了林分密度和平均胸径外,树种组成和竞争对进界也具有重要的影响[14, 17, 27]

    • 立地因子对进界的影响,主要包括坡度(slope degree)、坡向(slope aspect)、海拔(elevation)及其交互作用的影响[6, 16, 28]。本研究选择了这3个指标值来分析其对进界株数的影响。

    • 气候可能会增加树木进界的可变性和不确定性[12],并引起一些学者的注意[29]。气象因子对进界的影响方面,主要考虑温度、降水量、霜冻及干旱等的影响[16, 23, 30]。本研究利用ClimateAP v2.10软件[31],通过输入各样地的空间位置坐标和海拔,然后获取各样地每年的各类气象因子变量,最后对5年调查间隔期的变量进行平均计算。主要气象因子见表2

      表 2  主要气象变量统计表

      Table 2.  Satistical table of main climate variables

      气象因子
      Climate factor
      最大值
      Max.
      最小值
      Min.
      平均值 (标准差)
      Mean (standard deviation)
      年份
      Year
      年平均温度
      Mean annual temperature (MAT)/℃
      6.56 1.8 4.2 (0.8) 1999—2004
      6.96 2.3 4.6 (0.8) 2004—2009
      6 1.3 3.7 (0.8) 2009—2014
      最暖月平均气温
      Mean temperature of the warmest month (MWMT)/℃
      23.58 18.6 21.1 (0.9) 1999—2004
      22.68 18.4 20.8 (0.8) 2004—2009
      22.78 18.3 20.8 (0.8) 2009—2014
      最冷月平均气温
      Mean temperature of the coldest month (MCMT)/℃
      − 11.5 − 18.4 − 15.7 (1.1) 1999—2004
      − 10.1 − 16.9 − 14.4 (1.1) 2004—2009
      − 12.1 − 19 − 16.5 (1.1) 2009—2014
      年平均降水量
      Mean annual precipitation (MAP)/mm
      920.6 498.6 646.3 (75.3) 1999—2004
      1 038.8 475.2 640.9 (113.3) 2004—2009
      1 139.2 518.2 705.3 (126.4) 2009—2014
      年平均夏季 (5月至9月)降水量
      Mean annual summer (from May to September) precipitation (MSP)/mm
      731.4 399.2 514.9 (53.7) 1999—2004
      827.6 379.2 512.9 (88.9) 2004—2009
      884.2 389.8 539.4 (99.0) 2009—2014
      无霜期天数
      Number of frost-free days (NFFD)
      196.8 156.6 176.6 (7.9) 1999—2004
      200.2 158.6 178.9 (8.1) 2004—2009
      194.6 154.8 174.0 (7.7) 2009—2014
      上一年8月至当年7月的降雪量
      Snowfall between August in previous year and July in current year (PAS)/ mm
      115.8 35.2 54.9 (16.3) 1999—2004
      123.2 35 60.9 (14.8) 2004—2009
      154.6 52.4 80.5 (19.7) 2009—2014
    • 标准的泊松概率密度函数为:

      $ {f_{\rm{P}}}(y) = \frac{{{\lambda ^y}{\rm{e}^{ - \lambda }}}}{{y!}} $

      (1)

      式中:$y$是进界数量的随机变量,$\lambda $$y$的数学期望和方差。

      应用一个标准的泊松分布来模拟进界的数量没有充分考虑数据中存在的大量零值问题。

      零膨胀泊松分布模型(ZIP)经常被用作考虑过量零值问题的模拟方程。在零膨胀泊松分布模型中,零的数量被分成了两个部分:第1部分是伯努利分布原理引起的零计数;第2部分是由泊松过程引起的零计数。在零膨胀泊松模型中,概率为1 − $p$的非零计数假定为泊松分布。具体的概率密度函数如下:

      $ {f_{\rm{ZIP}}}(y) = \left\{ {\begin{array}{*{20}{l}} {p + (1 - p){\rm{e}^{ - \lambda }},\;\;{\rm{ }}y = 0}\\ {(1 - p)\dfrac{{{\lambda ^y}{\rm{e}^{ - \lambda }}}}{{y!}},\;\;{\rm{ }}y > 0} \end{array}} \right. $

      (2)
    • 虽然ZIP模型足以模拟包含大量零值的数据,但是当数据表现出过度零散时,就不能很好的适应了,因为泊松分布需要严格的均值和方差相等。不同于泊松模型,负二项模型包含一个可加参数来解释过度离散。负二项分布最简单的概率密度函数:

      $ {f_{\rm{NB}}}(y) = \frac{{\Gamma \left(y + \dfrac{1}{k}\right)}}{{\Gamma \left(\dfrac{1}{k}\right)y!}}{\left(\frac{1}{{\mu k + 1}}\right)^{1/k}}{\left(\frac{{\mu k}}{{\mu k + 1}}\right)^y} $

      (3)

      式中:$\Gamma $为伽马函数,$k$$\;\mu$是待估参数[32]

      负二项分布的方差是${\rm{Var}}(y) = \mu + k{\mu ^2}$。当$k$趋于0时,NB分布就变成了泊松分布。与ZIP相类似,零膨胀负二项分布(ZINB)的概率密度函数如下:

      $ {f_{\rm{ZINB}}}(y) = \left\{ {\begin{array}{*{20}{l}} {p + \left(1 - p\right){{\left(\dfrac{1}{{\mu k + 1}}\right)}^{1/k}},\;\;\;{\rm{ }}y = 0}\\ {(1 - p){f_{\rm{NB}}}(y),\;\;\;y > 0} \end{array}} \right. $

      (4)
    • 在广义的线性模型(GLM)中,包括基本的3个部分:(1)因变量y的分布;(2)连接函数;(3)系统的线性预测[32]。在本研究中,因变量的分布,对于计数部分,采取泊松分布或者是负二项分布;对于零部分,采取伯努利分布或者是伯努利分布与泊松分布或者是负二项分布的结合。连接函数是一个把线性预测值与因变量的希望值结合起来。

      影响进界数量的关键因素是与林分因子、立地因子及气象因子相关的。标准的泊松分布和零膨胀泊松分布的线性预测为:

      $\ln \lambda = X\beta = {\beta _0} + {\beta _1}{x_1} + {\beta _2}{x_2} + \cdot \cdot \cdot+ {\beta _n}{x_n}$

      (5)

      对于标准的负二项分布和零膨胀负二项分布模型,对数连接函数如下:

      $\ln \mu = X\beta = {\beta _0} + {\beta _1}{x_1} + {\beta _2}{x_2} + \cdot \cdot \cdot + {\beta _n}{x_n}$

      (6)

      在伯努利分布中的参数$p$,利用logistic模型来估计:

      $ p{\rm{ = }}\left( {\dfrac{1}{{1{\rm{ + }}{{\rm{e}}^{ - ({\alpha _0} + {\alpha _1}{x_1} + {\alpha _2}{x_2} + \cdot \cdot \cdot + {\alpha _n}{x_n})}}}}} \right) $

      (7)

      式(5)、(6)和(7)3个模型中的自变量来自于林分因子、立地因子及气象因子经过逐步回归后具有显著影响的因子。在模拟过程中,由于存在着嵌套结构,以及时间序列相关性,需要考虑样地的随机效应。零膨胀模型包括两部分随机效应,一部分是在伯努利分布上,另一部分是在服从泊松或负二项分布的计数部分。本研究在两部分均考虑了随机截距效应,服从均值为0的正态分布。

      本研究选择了林分因子、立地因子和气象因子作为协变量。具体的林分因子变量包括公顷株数、公顷断面积、林分算数平均直径、林分平方平均直径、林龄以及郁闭度。立地因子包括坡度、坡向、海拔及变换形式SIC和CE等变量。气象因子选择年平均温度(MAT)、最暖月平均气温(MWMT)、最冷月平均气温(MCMT)、年平均降水量(MAP)、年平均夏季(5月至9月)降水量(MSP)、无霜期天数(NFFD)和上一年8月至当年7月的降雪量(PAS)等因子。利用逐步回归方法选择对进界有显著影响的因子($\alpha = 0.05$)。整个计算和模拟过程利用SAS9.4软件的SAS NLMIXED模块。

    • 采用AIC(akaike information criterion)、BIC(Bayesian information criterion)和− 2×对数似然值(− 2loglikelihood,− 2logL)这3个指标来比较不同模型间的模拟效果[33]。3个值越小,表明模拟效果越好。对于两个模型间模拟效果差异是否显著,采用似然比卡方检验(LRT)来评价($\alpha = 0.05$)。在模型模拟结束后,要利用验证数据对模型进行精度验证,由于考虑了样地的随机效应,则首先要计算出各样地的随机效应参数值。随机效应参数计算方法可参考文献[34],具体的计算公式:

      ${\hat { b}_i} = \hat { D}\hat { Z}_i^{\rm{T}}{\left({\hat { R}_i} + {\hat { Z}_i}\hat { D}\hat { Z}_i^{\rm{T}}\right)^{ - 1}}{ {\hat{ \varepsilon} _i}}$

      (8)

      式中:${\hat{ b}_i}$是带有方差协方差矩阵的随机效应向量,$\hat { D}$为样地间随机效应方差协方差矩阵,${\hat { R}_i}$为样地$i$的误差效应方差协方差矩阵,${\hat{ \varepsilon} _i}$为样地$i$残差向量,${\hat { Z}_i}$为样地$i$设计矩阵。

      根据验证数据和模拟结果,计算出每个样地的进界株数,然后选择确定系数(${R^2}$)、均方根根误差($\rm{RMSE}$)和平均绝对残差($|\bar E|$)3个模型精度评价指标对验证结果进行效果评价[35]

      ${R^2}{\rm{ = 1 - }}\frac{{\displaystyle\sum\limits_{i = 1}^m {{{(\rm{est} - \rm{mea})}^2}} }}{{\displaystyle\sum\limits_{i = 1}^m {{{(\rm{mea} - \overline{mea})}^2}} }}$

      (9)

      $ \begin{array}{c} \;\\ \rm{RMSE} = \sqrt {\dfrac{{\displaystyle\sum\limits_{i = 1}^m {{{(\rm{est} - \rm{mea})}^2}} }}{{m - 1}}} \end{array}$

      (10)

      $\left|\bar E\right| = \sum\limits_{i = 1}^m {|\rm{est} - \rm{mea}|} /m$

      (11)

      式中:$\rm{est}$为估计值,$\rm{mea}$为观测值,$\rm{\overline{mea} }$为观测平均值,$m$为样地数量。

    • 分别选择泊松分布和负二项分布进行模拟,结果见表3。研究表明:对进界有明显影响的因子包括林分算数平均直径和林分公顷株数。而所有的立地因子和气象因子均对进界没有显著影响。随着林分算数平均直径和公顷株数的增加,进界概率明显降低。与标准泊松分布相比,标准负二项分布的AIC值、BIC值和− 2logL值明显降低,说明数据存在着过度离散问题,采取标准负二项分布效果更好。在所有参与建模的236块样地中,没有发生进界的样地占比接近22%。这意味着数据存在着过量零值的可能。因此,进一步选择零膨胀泊松分布模型和零膨胀负二项分布模型进行模拟。结果表明,无论是泊松分布还是负二项分布,加入零膨胀模型后,模型的模拟精度明显提高(AIC值、BIC值和− 2logL值越来越小)。

      表 3  林分进界模型模拟结果

      Table 3.  Simulation results of stand-level recruitment models

      参数
      Parameter
      基础模型 Basic model
      M1M2M3M4
      不考虑随机效应
      No random effect
      模型的零部分
      Zero component of the model
      $ {\alpha _0}$ 5.730 8 (0.990 3)*** 5.726 6 (1.000 1)***
      $ {\alpha _1}$ − 0.203 8 (0.046 4)*** − 0.202 4 (0.046 8)***
      $ {\alpha _2}$ − 1.110 0 (0.265 1)*** − 1.110 3 (0.267 7)***
      计数的部分
      Positive count component of the model
      $\; {\beta _0}$ 9.713 9 (0.026 7)*** 9.237 8 (0.027 6)*** 8.446 8 (0.469 7)*** 8.065 1 (0.295 0)***
      $\; {\beta _1}$ − 0.313 0 (0.001 7)*** − 0.274 4 (0.001 7)*** − 0.214 1 (0.021 4)*** − 0.187 4 (0.013 8)***
      $\; {\beta _2}$ − 0.756 0 (0.008 3)*** − 0.621 9 (0.008 5)*** − 0.771 5 (0.149 9)*** − 0.583 3 (0.095 2)***
      $ k$ 2.437 2 0.753 2
      AIC 55 314 39 871 4 877.6 4 671.8
      BIC 55 327 39 896 4 894.2 4 700.9
      − 2logL 55 308 39 859 4 869.6 4 657.8
      参数
      Parameter
      考虑混合效应模型 With mixed effect model
      M5 M6 M7 M8
      考虑随机效应部分
      With random effect
      模型的零部分
      Zero component of the model
      $ {\alpha _0}$ 5.729 5 (0.990 3)*** 5.716 6 (0.993 3)***
      $ {\alpha _1}$ − 0.203 8 (0.046 4)*** − 0.202 6 (0.046 5)***
      $ {\alpha _2}$ − 1.109 8 (0.265 2)*** − 1.108 0 (0.265 9)***
      计数的部分
      Positive count component of the model
      $\; {\beta _0}$ 10.453 3 (0.150 7)*** 9.587 2 (0.113 4)*** 8.995 2 (0.552 7)*** 8.263 9 (0.354 6)***
      $\; {\beta _1}$ − 0.391 3 (0.006 4)*** − 0.259 9 (0.006 3)*** − 0.240 6 (0.026 7)*** − 0.202 0 (0.017 8)***
      $\; {\beta _2}$ − 0.992 8 (0.017 5)*** − 1.220 0 (0.020 8)*** − 0.926 2 (0.169 5)*** − 0.687 1 (0.107 8)***
      $ k$ 2.292 2 0.525 7
      AIC 21 504 15 054 4 882.4 4 645.0
      BIC 21 518 15 078 4 899.8 4 672.7
      − 2logL 21 496 15 040 4 872.4 4 629.0
      随机效应方差协方差矩阵 Covariance matrix of random effect variance 2.976 0 0.673 0 0.000 1 0.254 7
      注:k是负二项分布的待估参数;M1为标准泊松分布进界模型,即公式 (1);M2为零膨胀泊松分布进界模型,即公式 (2);M3为标准负二项分布进界模型,即公式 (3);M4为零膨胀负二项分布进界模型,即公式 (4);M5为考虑随机效应的标准泊松分布进界模型;M6为考虑随机效应的零膨胀泊松分布进界模型;M7为考虑随机效应的标准负二项分布进界模型;M8为考虑随机效应的零膨胀负二项分布进界模型。表中括号内的值为标准差;*为0.01 < P < 0.05,**为0.001 < P < 0.01,***P < 0.001;$ {\alpha _0}$和$\; {\beta _0}$为截距参数,$ {\alpha _1}$和$\; {\beta _1}$为林分算数平均直径参数值,$ {\alpha _2}$和$\; {\beta _2}$为林分公顷株数参数值。Notes: k is a parameter to be estimated in negative binomial distribution; M1 is Poisson distribution recruitment model, i.e. formula (1); M2 is zero-inflated Poisson distribution recruitment model, i.e. formula (2); M3 is negative binomial distribution recruitment model, i.e. formula (3); M4 is zero-inflated negative binomial distribution recruitment model, i.e. formula (4); M5 is Poisson distribution recruitment model based on mixed effect model; M6 is zero-inflated Poisson distribution recruitment model based on mixed effect model; M7 is negative binomial distribution recruitment model based on mixed effect model; M8 is zero-inflated negative binomial distribution recruitment model based on mixed effect model. Values in brackets are standard deviation; * is 0.01 < P < 0.05,** is 0.001 < P < 0.01,*** is P < 0.001; $ {\alpha _0}$ and $\; {\beta _0}$ are intercept parameters, $ {\alpha _1}$ and $\; {\beta _1}$ are stand arithmetic average diameter parameters, $ {\alpha _2}$ and $\; {\beta _2}$ are stand number of per hectare parameters.
    • 选择标准泊松分布、标准负二项分布、零膨胀泊松分布模型和零膨胀负二项分布模型作为基础模型,在此基础上考虑样地水平的随机截距效应。由于零膨胀模型在伯努利分布的logistic方程上随机效应不收敛,因此只考虑了计数部分的随机效应。在计数部分只有截距效应收敛。研究结果(表3)表明,考虑了样地的随机效应后,除了标准负二项分布模拟结果外,标准泊松分布模型、零膨胀泊松分布模型和零膨胀负二项分布模型的模拟效果明显提高(AIC、BIC和− 2logL值都明显降低)。进一步利用LRT对模型之间的差异进行比较评价,结果见表4

      表 4  利用LRT指标对模型进行比较的结果

      Table 4.  Comparing results of the model based on LRT index

      模型 ModelLRT值 LRT valueP模型 ModelLRT 值 LRT valueP模型 ModelLRT值 LRT valueP
      M1/M2 15 449 < 0.000 1 M1/M5 33 812 < 0.000 1 M3/M4 211.8 < 0.000 1
      M2/M6 24 819 < 0.000 1 M3/M7 2.8 > 0.05 M4/M8 28.8 < 0.000 1

      表4中可以看出,除了标准负二项分布外,标准泊松分布模型、零膨胀泊松分布模型和零膨胀负二项分布模型的模拟效果都比不考虑随机效应效果要好,并且均达到显著水平($\alpha = 0.05$P < 0.000 1)。其中考虑样地随机效应的零膨胀负二项分布模型的模拟精度最高。

    • 利用确定系数(${R^2}$)、均方根误差(RMSE)和平均绝对残差($|\bar E|$)3个模型精度评价指标对验证数据进行计算。其中不考虑随机效应的模型(M1、M2、M3、M4)直接利用表3的模拟结果。而考虑随机效应的模型(M5、M6、M7、M8),首先利用式(8)计算出每个样地的随机效应值,然后再计算3个评价指标。具体计算结果见表5

      表 5  模型验证结果

      Table 5.  Validation results of models

      评价指标
      Evaluating index
      不考虑随机效应 No random effect考虑随机效应 With random effect
      M1M2M3M4M5M6M7M8
      $ {R^2}$ 0.40 0.55 0.50 0.65 0.75 0.91 0.54 0.84
      RMSE/株 RMSE/tree
      163.8 146.2 158.7 139.2 125.3 161.8 153.8 116.5
      $ |\bar E|$/株 $ |\bar E|$/tree 98.8 95.9 83.9 81.7 65 70.2 85 61.7

      表5中可以看出,利用验证数据对模型的拟合效果进行评价,虽然与模拟数据的规律不完全一致,大体上随着模型的模拟效果提高,模型的相关系数提高,而均方根误差和平均绝对残差值降低。在考虑随机效应后,与不考虑随机效应模型相比,精度都有所提高。模型验证精度最高的是考虑样地随机效应的零膨胀负二项分布模型,这与表3中的模拟结果一致。因此,最后选择的最佳林分进界模型形式及概率分布形式如下。

      模型形式:

      $\ln \mu = X\beta = 8.263 \; 9 - 0.202 \; 0{x_1} - 0.687 \; 1{x_2}$

      分布密度函数服从零膨胀负二项分布形式(式(4)),其中具体参数值如下:

      $ \begin{split} &p = \dfrac{1}{{1 + \exp ( - (5.716 \; 6 - 0.202 \; 6{x_1} - 1.108 \; 0{x_2}))}}, \\ &k = 0.525 \; 7 \end{split} $

      式中:${x_1}$为林分算数平均直径,${x_2}$为林分公顷株数。

    • (1)作为林分密度指标,林分公顷株数与林木的进界呈反比,随着林分公顷株数的增加,林分内林木的进界概率和进界株数都明显降低。林分公顷株数过大,对于光照、营养和空间的竞争加剧,更小的树木容易枯损,很难进界。林分算数平均直径与林木的进界呈反比,随着林分算数平均直径的增加,林分内林木的进界概率和进界株数明显降低。林分算数平均直径越大,说明上层林分占有空间能力越大,也不利于更小的树木发生进界。因此,为了促进树木的进界,就要适时采取科学的经营措施。本研究中,立地因子和气象因子中的各类因子对进界的影响均不显著,可能的原因是本研究中选择的样地主要分布在吉林省东南部,样地之间的立地和气象等因子差异小,不足以产生大的影响。这两个因子对林分进界虽然没有显著影响,但是并不代表不重要。在今后的研究中,如果数据量足够大,或者是样地范围广,则构建进界模型时还需要考虑气象因子和立地因子,再进行模拟。

      (2)负二项分布模型的模拟和验证精度都要好于泊松分布,说明数据中存在过度离散问题,在实际中是不能忽视的。无论是泊松分布还是负二项分布模型,在考虑数据中存在过量零值问题后,即考虑加入零膨胀模型后,模型的模拟精度都显著提高。这说明样地中存在着大量的零值包括两种情况,一种是分布中的零值问题,一种是固有的零值问题。除了负二项分布模型外,与不考虑随机效应模型相比,标准泊松分布、零膨胀泊松分布和零膨胀负二项分布在考虑随机效应后,模拟精度明显提高。同时,利用验证数据进行验证也同样印证了此结论。

      (3)与一些学者结论相一致,在不考虑随机效应时,零膨胀负二项分布模型的模拟效果最好[6, 12],当考虑随机效应后,零膨胀负二项分布模型的模拟效果最好[2],利用验证数据进行验证也支持此结论。但是,Zell等[23]在模拟进界时,认为零膨胀泊松分布在考虑随机效应后的模拟精度最高。造成两种不同结论的原因可能是数据本身离散程度的问题。零改变模型(zero-alter model)也经常被用来模拟林木的进界[2, 6, 12]。本研究也采取了零改变模型进行模拟,但是和零膨胀模型的模拟精度几乎相同,因此没有列入结果中。

      (4)很多因子对林木的进界有影响,并不仅仅限于林分因子、立地因子和气象因子的影响。采伐对进界也具有重大的影响,影响的效果取决于采伐的方式、采伐时间以及不同树种对采伐的响应[12, 36-37]。随着以后研究数据的逐渐完善,研究区域逐渐扩大,还需考虑更多因子的影响,以及这些因子之间的交互作用。

    • 本研究基于零膨胀模型及混合效应模型相结合的方法,采用泊松分布和负二项分布两种形式,利用逐步回归参数选择方法,构建了基于林分因子、立地因子及气象因子的吉林省东南部蒙古栎林林分水平进界模型。模拟结果表明,林分算数平均直径和林分公顷株数是影响林木进界概率和数量最重要的影响因子,并且均与林木进界概率和进界的数量呈反比。立地和气象因子中的所有指标均对进界没有产生明显影响。负二项分布模型由于考虑了数据过度离散问题,模拟精度要高于泊松分布;在考虑样地的随机效应后,除了标准负二项分布模型外所有模型都明显提高了模型的模拟精度;同时考虑随机效应和零膨胀的负二项分布模型,其模型的模拟效果最好,精度最高。

参考文献 (37)

目录

    /

    返回文章
    返回