高级检索

留言板

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

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

基于贝叶斯方法的蒙古栎林单木枯死模型

姚丹丹 徐奇刚 闫晓旺 李玉堂

姚丹丹, 徐奇刚, 闫晓旺, 李玉堂. 基于贝叶斯方法的蒙古栎林单木枯死模型[J]. 北京林业大学学报, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260
引用本文: 姚丹丹, 徐奇刚, 闫晓旺, 李玉堂. 基于贝叶斯方法的蒙古栎林单木枯死模型[J]. 北京林业大学学报, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260
Yao Dandan, Xu Qigang, Yan Xiaowang, Li Yutang. Individual-tree mortality model of Mongolian oak forests based on Bayesian method[J]. Journal of Beijing Forestry University, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260
Citation: Yao Dandan, Xu Qigang, Yan Xiaowang, Li Yutang. Individual-tree mortality model of Mongolian oak forests based on Bayesian method[J]. Journal of Beijing Forestry University, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260

基于贝叶斯方法的蒙古栎林单木枯死模型

doi: 10.13332/j.1000-1522.20180260
基金项目: 林业行业公益性科研项目(201504303)
详细信息
    作者简介:

    姚丹丹。主要研究方向:森林经理学。Email:371352079@qq.com  地址:100091北京市海淀区香山路东小府1号中国林业科学研究院资源信息研究所

    通讯作者:

    闫晓旺,教授级高级工程师。主要研究方向:森林经理学。Email:ccyxw2005@163.com  地址:130022吉林省长春市人民大街4756号吉林省林业调查规划院

  • 中图分类号: S758.5

Individual-tree mortality model of Mongolian oak forests based on Bayesian method

  • 摘要: 目的贝叶斯统计法在提高模型参数稳定性上有较大的优势,研究贝叶斯方法在单木枯死模型中的应用,改进模型参数的估计方法,为蒙古栎天然林林分生长收获与经营管理提供参考。方法以蒙古栎天然异龄林为对象,基于202块固定样地数据,利用二分类Logistic模型构建基于经典概率统计法、贝叶斯法和分层贝叶斯法的蒙古栎单木枯死模型。随机抽取80%的数据用于建立模型,剩下的20%用于检验模型,利用经典概率统计法(非线性最小二乘法)、有先验信息的贝叶斯统计法和无先验信息的分层贝叶斯统计法进行参数估计,分析模型的表现和参数分布。模型的拟合效果通过计算ROC曲线下的面积AUC(Under Curve)来判断,并利用Pearson-χ2检验来检验模型的拟合优度。结果(1)贝叶斯法与传统极大似然法的估计值相近,且其估计参数的标准差小于传统方法。(2)贝叶斯法估计参数的可信区间最小,比传统极大似然法的置信区间小6.0% ~ 31.8%。层次贝叶斯法估计参数的可信区间最大,比传统极大似然法的置信区间大11.2% ~ 185.0%。(3)拟合效果最好的是层次贝叶斯法,其模型AUC值为0.83,贝叶斯法与传统极大似然法模型的AUC值均为0.73。结论层次贝叶斯法在拟合枯死模型方面具有明显的优势,拟合效果最好,模型预估精度最高。
  • 图  1  单木枯死模型参数的贝叶斯估计的95%可信区间

    C:极大似然法;B:贝叶斯法;HB:层次贝叶斯法。下同。C represents estimated results of maximum likelihood method; B represents estimated results of Bayesian method; HB represents estimated results of Hierarchical Bayesian method. Same as below.

    Figure  1.  Bayesian 95% credible and classical 95% confidence intervals for parameters in mortality mode

    图  2  枯死模型的贝叶斯、层次贝叶斯参数的后验分布和极大似然估计参数的似然分布

    Figure  2.  Posterior distribution of Bayesian method and likelihood distribution of maximum likelihood method for morality model

    图  3  3种方法拟合枯死模型的ROC曲线

    Figure  3.  ROC curves of mortality model based on three methods

    表  1  建模数据和检验数据样地基本因子统计量

    Table  1.   Summary statistics of sample plots for calibration and validation data

    项目
    Item
    株数
    Number of trees
    自变量
    Independent variable
    平均值
    Average
    最大值
    Maximum
    最小值
    Minimum
    标准差
    Standard deviation
    建模数据 Calibration data 11 118 D 12.386 5 86.806 5 5.000 0 7.456 7
    BAL 0.846 3 2.040 9 0.000 0 0.397 1
    N 1 599 3 950 467 675
    Nm 135 733 17 134
    检验数据 Validation data 2 780 D 12.155 9 73.690 9 5.003 6 7.440 2
    BAL 0.841 0 1.881 1 0.000 0 0.390 6
    N 1 482 3 950 217 698
    Nm 202 700 17 135
    注:D为林木胸径(cm),BAL为林木胸高断面积(m2/hm2),N为林分密度(株/hm2),Nm为林分枯死株数(株/hm2)。Notes: D represents DBH (cm); BAL represents basal area of breast height of forest trees (m2/ha); N represents stand density (tree/ha); Nm represents number of dead trees (tree/ha).
    下载: 导出CSV

    表  2  参数的先验信息

    Table  2.   Prior distribution of each parameter for base model and mixed-effects model

    模型 Model参数 Parameter先验分布 Prior distribution
    基础模型
    Base model
    a0 a0 ~ N (−2.148 3, 0.149 52)
    a1 a1 ~ N (−0.126, 0.010 52)
    a2 a2 ~ N (1.294 3, 0.093 92)
    a3 a3 ~ N (0.176 3, 0.007 32)
    a4 a4 ~ N (−1.893 9, 0.444 82)
    混合模型
    Model with
    random effects
    a0 a0 ~ N (0, 10 002)
    a1 a1 ~ N (0, 10 002)
    a2 a2 ~ N (0, 10 002)
    a3 a3 ~ N (0, 10 002)
    a4 a4 ~ N (0, 10 002)
    σa0 (1/σa0)2 ~ Gamma (0.001, 0.001)
    σ (1/σ)2 ~ Gamma (0.001, 0.001)
    下载: 导出CSV

    表  3  观测值和预测值的列联表

    Table  3.   Contingency table of the observed values and predicted one

    项目 Item预测值 Predicted value总和 Total
    枯死 Dead存活 Alive
    观测值 Observed value 枯死 Dead a11 a12 A1m
    存活 Alive a21 a22 A2m
    和 Total An1 An2 Anm
    敏感度 Sensitivity=a11/A1m, 特异度 Specificity=a22/A2m
    注:a11表示正确预测枯死林木的株数,a12表示错误预测枯死林木的株数,a21表示错误预测林木存活的株数,a22表示正确预测存活林木的株数,A1m表示枯死实测值,A2m表示存活实测值,An1表示枯死预测值,An2表示存活预测值。Notes: a11 represents the number of dead trees which is correctly predicted, a12 represents the number of dead trees which is not correctly predicted, a21 represents the number of alive trees which is not correctly predicted, a22 represents the number of alive trees which is correctly predicted, A1m represents the measured value of dead trees, A2m represents the measured value of alive trees, An1 represents the predicted value of dead trees, An2 represents the predicted value of alive trees.
    下载: 导出CSV

    表  4  基于建模数据的参数估计结果

    Table  4.   Parameter estimates of classical and Bayesian methods using calibration

    方法
    Method
    参数
    Parameter
    估计值
    Estimated value
    标准差
    Standard deviation
    置信区间 Credible interval
    2.50%97.50%
    最大似然 Maximum likelihood a0 −2.148 3 0.149 5 −2.447 2 −1.849 4
    a1 −0.126 0 0.010 5 −0.147 1 −0.104 9
    a2 1.294 3 0.093 9 1.106 6 1.482 1
    a3 0.176 3 0.073 2 0.029 8 0.322 8
    a4 −1.893 9 0.444 8 −2.783 5 −1.004 3
    贝叶斯统计 Bayesian method a0 −2.136 0 0.138 8 −2.413 6 −1.858 4
    a1 −0.123 1 0.007 2 −0.137 5 −0.108 7
    a2 1.271 0 0.078 8 1.113 3 1.428 7
    a3 0.136 1 0.062 3 0.011 6 0.260 6
    a4 −1.886 0 0.415 4 −2.716 8 −1.055 2
    层次贝叶斯 Hierarchical bayesian method a0 −2.641 0 0.289 7 −3.248 0 −2.176 0
    a1 −0.145 4 0.013 0 −0.168 6 −0.118 0
    a2 1.595 0 0.176 0 1.311 0 1.963 0
    a3 0.274 8 0.084 9 0.108 7 0.444 0
    a4 −3.082 0 1.337 0 −5.523 0 −0.452 2
    σa0 1.030 0 0.158 9 0.742 8 1.360 0
    下载: 导出CSV

    表  5  3种方法拟合的枯死模型误差统计量

    Table  5.   Error statistics of morality model by three methods

    方法 Method建模数据 Calibration data验证数据 Validation data
    枯死数
    Number of
    dead trees
    预测枯死数
    Predicted number
    of dead trees
    χ2AUCAIC (DIC)枯死数
    Number of
    dead trees
    预测枯死数
    Predicted number
    of dead trees
    χ2AUC (DIC)
    传统方法
    Classical method
    1 338 1 349 0.09 0.734 7 407.4 375 404 2.08 0.726
    贝叶斯统计
    Bayesian method
    1 338 1 336 0.01 0.734 7 407.21 375 394 0.92 0.727
    层次贝叶斯统计
    Hierarchical Bayesian method
    1 338 1 335 0.01 0.825 6 771.26 375 389 0.50 0.800
    下载: 导出CSV
  • [1] 刘平, 马履一, 贾黎明, 等. 油松林木枯损率模型研究[J]. 林业资源管理, 2008(2):51−56. doi:  10.3969/j.issn.1002-6622.2008.02.012

    Liu P, Ma L Y, Jia L M, et al. Study on tree mortality model for Pinus tabulaeformis plantation[J]. Forest Resources Management, 2008(2): 51−56. doi:  10.3969/j.issn.1002-6622.2008.02.012
    [2] Fahey T J, Battles J J, Wilson G F. Responses of early successional northern hardwood forests to changes in nutrient availability[J]. Ecological Monographs, 1998, 68: 183−212. doi:  10.1890/0012-9615(1998)068[0183:ROESNH]2.0.CO;2
    [3] Kobe R K. Intraspecific variation in sapling mortality and growth predicts geographic variation in forest composition[J]. Ecological Monograph, 1996, 66: 181−201. doi:  10.2307/2963474
    [4] Hamilton D. A logistic model of mortality in thinned andunthinned mixed conifer stands of northern Idaho[J]. Forest Science, 1986, 32: 989−1000.
    [5] Wyckoff P H, Clark J S. Predicting tree mortality from diameter growth: a comparison of maximum likelihood and Bayesian approaches[J]. Canadian Journal of Forest Research, 2000, 30(1): 156−167. doi:  10.1139/x99-198
    [6] Metcalf C J E, McMahon S M, Clark J S. Overcoming data sparseness and parametric constraints in modeling of tree mortality: a new nonparametric Bayesian model[J]. Canadian Journal of Forest Research, 2009, 39(9): 1677−1687. doi:  10.1139/X09-083
    [7] Lu L, Wang H, Chhin S, et al. A Bayesian model averaging approach for modelling tree mortality in relation to site, competition and climatic factors for Chinese fir plantations[J]. Forest Ecology and Management, 2019, 440: 169−177.
    [8] Green E J, Roesch F A, Smith A F M, et al. Bayesian estimation for the three-parameter Weibull distribution with tree diameter data[M]. Washington: International Biometric Society, 1994.
    [9] Bullock B P, Boone E L. Deriving tree diameter distributions using Bayesian model averaging[J]. Forest Ecology and Management, 2007, 242(2−3): 127−132. doi:  10.1016/j.foreco.2007.01.024
    [10] Li R X, Stewart B, Weiskittel A. A Bayesian approach for modelling non-linear longitudinal/hierarchical data with random effects in forestry[J]. Forestry, 2012, 85(1): 17−25. doi:  10.1093/forestry/cpr050
    [11] Green E J, Strawderman W E. A Bayesian growth and yield model for slash pine plantations[J]. Journal of Applied Statistics, 1996, 23(2−3): 285−300. doi:  10.1080/02664769624251
    [12] Nystrom K, Stahl G. Forecasting probability distributions of forest yield allowing for a Bayesian approach to management planning[J]. Silva Fennica, 2001, 35(2): 185−201.
    [13] 张雄清, 张建国, 段爱国. 基于贝叶斯法估计杉木人工林树高生长模型[J]. 林业科学, 2014, 50(3):69−75.

    Zhang X Q, Zhang J G, Duan A G. Tree-height growth model for Chinese fir plantation based on Bayesian method[J]. Scientia Silvae Sinicae, 2014, 50(3): 69−75.
    [14] Zhang X Q, Duan A G, Zhang J G. Estimating tree height-diameter models with the Bayesian method[J/OL]. The Scientific World Journal, 2014: 1−9 [2018−01−06]. https://www.ncbi.nlm.nih.gov/pubmed/24711733.
    [15] Tatsumi S, Owari T. Bayesian modeling of neighborhood competition in uneven-aged mixed-species stands[J]. Formath, 2013, 12: 191−209. doi:  10.15684/formath.12.191
    [16] 马武. 蒙古栎林单木生长模型系研究[D]. 北京: 中国林业科学研究院, 2012.

    Ma W. Growth model for individual-tree in natural Quercus mongolica forests[D]. Beijing: Chinese Academy of Forestry, 2012.
    [17] Fang Z, Bailey R L. Nonlinear mixed effects modeling for slash pine dominant height growth following intensive silvicultural treatments[J]. Forest Science, 2001, 47(3): 287−300.
    [18] Calegario N, Daniels R F, Maestri R, et al. Modeling dominant height growth based on nonlinear mixed-effects model: a clonal eucalyptus plantation case study[J]. Forest Ecology and Management, 2005, 204(1): 11−21. doi:  10.1016/j.foreco.2004.07.051
    [19] Patenaud G, Milne R, Van Oijen M, et al. Integrating remote sensing datasets into ecological modelling: a Bayesian approach[J]. International Journal of Remote Sensing, 2008, 29(5): 1295−1315. doi:  10.1080/01431160701736414
    [20] Klemedtsson L, Jansson P E, Gustafsson D, et al. Bayesian calibration method used to elucidate carbon turnover in forest on drained organic soil[J]. Biogeochemistry, 2008, 89(1): 61−79. doi:  10.1007/s10533-007-9169-0
    [21] Russell M B. Influence to prior distributions and random effects on count regression models: implications for estimating standing dead tree abundance[J]. Environment and Ecological Statistics, 2015, 22: 145−160. doi:  10.1007/s10651-014-0290-7
  • [1] 陈国栋, 杜研, 丁佩燕, 郭珂歆, 尹忠东.  基于混合效应模型的新疆天山云杉单木胸径预测模型构建 . 北京林业大学学报, 2020, 42(7): 12-22. doi: 10.12171/j.1000-1522.20190236
    [2] 娄明华, 张会儒, 雷相东, 白超, 杨同辉.  天然栎类阔叶混交林林分平均高与平均胸径关系模型 . 北京林业大学学报, 2020, 42(9): 37-50. doi: 10.12171/j.1000-1522.20190463
    [3] 李春明, 李利学.  基于零膨胀模型及混合效应模型相结合的蒙古栎林林木进界模拟研究 . 北京林业大学学报, 2020, 42(6): 59-67. doi: 10.12171/j.1000-1522.20190216
    [4] 胡雪凡, 张会儒, 周超凡, 张晓红.  不同抚育间伐方式对蒙古栎次生林空间结构的影响 . 北京林业大学学报, 2019, 41(5): 137-147. doi: 10.13332/j.1000-1522.20190037
    [5] 欧强新, 雷相东, 沈琛琛, 宋国涛.  基于随机森林算法的落叶松−云冷杉混交林单木胸径生长预测 . 北京林业大学学报, 2019, 41(9): 9-19. doi: 10.13332/j.1000-1522.20180266
    [6] 李艳丽, 杨华, 邓华锋.  蒙古栎−糠椴天然混交林空间格局研究 . 北京林业大学学报, 2019, 41(3): 33-41. doi: 10.13332/j.1000-1552.20180236
    [7] 曹梦, 潘萍, 欧阳勋志, 臧颢, 吴自荣, 杨阳, 占常燕.  基于哑变量的闽楠天然次生林单木胸径和树高生长模型研究 . 北京林业大学学报, 2019, 41(5): 88-96. doi: 10.13332/j.1000-1522.20190026
    [8] 王涛, 董利虎, 李凤日.  基于混合效应的杂种落叶松人工幼龄林单木枯损模型 . 北京林业大学学报, 2018, 40(10): 1-10. doi: 10.13332/j.1000-1522.20170437
    [9] 娄明华, 张会儒, 雷相东, 卢军.  天然云冷杉针阔混交林单木胸径树高空间自回归模型研究 . 北京林业大学学报, 2016, 38(8): 1-9. doi: 10.13332/j.1000-1522.20150491
    [10] 邱志明, 苏声欣, 钟智昕, 唐盛林, 林谦佑.  柳杉人工林行列疏伐异龄混交林经营研究 . 北京林业大学学报, 2015, 37(3): 44-54. doi: 10.13332/j.1000-1522.20140317
    [11] 张西, 贾黎明, 张瑜, 郑聪慧.  基于FVS的秦岭地区栓皮栎天然次生林单木模型构建 . 北京林业大学学报, 2015, 37(5): 19-29. doi: 10.13332/j.1000-1522.20140356
    [12] 李响, 甄贞, 赵颖慧.  基于局域最大值法单木位置探测的适宜模型研究 . 北京林业大学学报, 2015, 37(3): 27-33. doi: 10.13332/j.1000-1522.20140313
    [13] 姚丹丹, 雷相东, 张则路.  基于贝叶斯法的长白落叶松林分优势高生长模型研究 . 北京林业大学学报, 2015, 37(3): 94-100. doi: 10.13332/j.1000-1522.20140221
    [14] 余黎, 雷相东, 王雅志, 杨英军, 王全军.  基于广义可加模型的气候对单木胸径生长的影响研究 . 北京林业大学学报, 2014, 36(5): 22-32. doi: 10.13332/j.cnki.jbfu.2014.05.007
    [15] 张琴, 林天喜, 王贵春, 孙国文, 范秀华.  红松、蒙古栎和色木槭凋落物混合分解研究 . 北京林业大学学报, 2014, 36(6): 106-111. doi: 10.13332/j.cnki.jbfu.2014.06.020
    [16] 李丹, 戴巍, 闫志刚, 王霓虹.  基于模糊层次分析法的落叶松人工林生境评价系统 . 北京林业大学学报, 2014, 36(4): 75-81. doi: 10.13332/j.cnki.jbfu.2014.04.015
    [17] 陈瑜, 徐程扬, 李乐, 蔡丽丽.  阔叶红松风景林单木景观质量评价与模型研究 . 北京林业大学学报, 2014, 36(5): 87-93. doi: 10.13332/j.cnki.jbfu.2014.05.006
    [18] 董灵波, 刘兆刚.  基于形态结构特征参数的樟子松人工林单木可视化研究 . 北京林业大学学报, 2011, 33(5): 20-27.
    [19] 李春明.  随机截距效应在模拟杉木人工林单木胸径生长量中的应用 . 北京林业大学学报, 2011, 33(4): 7-12.
    [20] 向玮, 雷相东, 刘刚, 徐光, 陈光法.  近天然落叶松云冷杉林单木枯损模型研究 . 北京林业大学学报, 2008, 30(6): 90-98.
  • 加载中
图(3) / 表 (5)
计量
  • 文章访问数:  800
  • HTML全文浏览量:  383
  • PDF下载量:  51
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-08-07
  • 修回日期:  2019-01-25
  • 网络出版日期:  2019-09-09
  • 刊出日期:  2019-09-01

基于贝叶斯方法的蒙古栎林单木枯死模型

doi: 10.13332/j.1000-1522.20180260
    基金项目:  林业行业公益性科研项目(201504303)
    作者简介:

    姚丹丹。主要研究方向:森林经理学。Email:371352079@qq.com  地址:100091北京市海淀区香山路东小府1号中国林业科学研究院资源信息研究所

    通讯作者: 闫晓旺,教授级高级工程师。主要研究方向:森林经理学。Email:ccyxw2005@163.com  地址:130022吉林省长春市人民大街4756号吉林省林业调查规划院
  • 中图分类号: S758.5

摘要: 目的贝叶斯统计法在提高模型参数稳定性上有较大的优势,研究贝叶斯方法在单木枯死模型中的应用,改进模型参数的估计方法,为蒙古栎天然林林分生长收获与经营管理提供参考。方法以蒙古栎天然异龄林为对象,基于202块固定样地数据,利用二分类Logistic模型构建基于经典概率统计法、贝叶斯法和分层贝叶斯法的蒙古栎单木枯死模型。随机抽取80%的数据用于建立模型,剩下的20%用于检验模型,利用经典概率统计法(非线性最小二乘法)、有先验信息的贝叶斯统计法和无先验信息的分层贝叶斯统计法进行参数估计,分析模型的表现和参数分布。模型的拟合效果通过计算ROC曲线下的面积AUC(Under Curve)来判断,并利用Pearson-χ2检验来检验模型的拟合优度。结果(1)贝叶斯法与传统极大似然法的估计值相近,且其估计参数的标准差小于传统方法。(2)贝叶斯法估计参数的可信区间最小,比传统极大似然法的置信区间小6.0% ~ 31.8%。层次贝叶斯法估计参数的可信区间最大,比传统极大似然法的置信区间大11.2% ~ 185.0%。(3)拟合效果最好的是层次贝叶斯法,其模型AUC值为0.83,贝叶斯法与传统极大似然法模型的AUC值均为0.73。结论层次贝叶斯法在拟合枯死模型方面具有明显的优势,拟合效果最好,模型预估精度最高。

English Abstract

姚丹丹, 徐奇刚, 闫晓旺, 李玉堂. 基于贝叶斯方法的蒙古栎林单木枯死模型[J]. 北京林业大学学报, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260
引用本文: 姚丹丹, 徐奇刚, 闫晓旺, 李玉堂. 基于贝叶斯方法的蒙古栎林单木枯死模型[J]. 北京林业大学学报, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260
Yao Dandan, Xu Qigang, Yan Xiaowang, Li Yutang. Individual-tree mortality model of Mongolian oak forests based on Bayesian method[J]. Journal of Beijing Forestry University, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260
Citation: Yao Dandan, Xu Qigang, Yan Xiaowang, Li Yutang. Individual-tree mortality model of Mongolian oak forests based on Bayesian method[J]. Journal of Beijing Forestry University, 2019, 41(9): 1-8. doi: 10.13332/j.1000-1522.20180260
  • 林木枯死和种群动态变化存在密切关系,也会对林分的结构与功能等方面产生影响,林分枯死模型可以在模型的基础上对森林演变规律进行分析,同时明确林分生长的变化情况以及影响因素[1]。目前林业领域关于林分的生长−死亡关系[2-3]的研究已经有很多,不过在具体实践中,需要基于固定样地中采集的数据确定出死亡率[4],而与此相关的数据采集难度较大。传统的林分枯死模型有一定局限性,无法对小生长率条件下的死亡概率进行分析描述,应用非参数模型是一种可行的解决方法。贝叶斯模型作为一种非参数模型,其特征表现为应用先验信息和样本信息进行统计分析,且具有不要求参数服从正态分布的优点。

    使用贝叶斯方法估计森林生长模型的参数近年来正逐渐成为林业学者们研究的热点,主要包括单木枯死模型[5-7]、直径分布[8-9]、树高曲线模型[10-11]、断面积分布模型[12]等。张雄清等[13]使用贝叶斯法估计了江西杉木(Cunninghamia lanceolata)人工林的树高生长模型参数,研究结果表明,利用贝叶斯法可以提高模型的可靠性和预测精度。同时又由于贝叶斯法能够利用先验信息与样本信息进行统计推断,在估计模型前可加入历史文献资料或专家主观经验,这一特点十分适合林业建模工作。Zhang等[14]分别运用传统法、无信息先验和有信息先验的贝叶斯方法对树高曲线模型的拟合效果进行比较分析,发现传统法、无先验信息和有先验信息贝叶斯方法拟合的模型检验结果相近,但是贝叶斯估计参数的可信区间比传统法更为集中[14]。另外,由于有先验信息的加入,贝叶斯法可以克服样本量小的问题,Wyckoff等[5]使用贝叶斯统计解决了林木死亡模型中数据稀少的问题。贝叶斯统计在国内生长模型方面的研究渐成趋势[12-13],但由于生长数据的相关和层次性,层次贝叶斯模型的应用不多见[13,15]

    本文以蒙古栎(Quercus mongolica)天然异龄林为研究对象,基于202块固定样地数据,利用二分类Logistic模型构建基于经典概率统计法、贝叶斯法和分层贝叶斯法的蒙古栎单木枯死模型,探讨贝叶斯统计法与分层贝叶斯法在单木枯死模型中的应用。

    • 本文的研究区域处于吉林省延边东北部,地理坐标为129°47′ ~ 131°57′ E、43°00′ ~ 44°00′ N。海拔在360 ~ 1 478 m之间,坡度较为平缓,一般在25°以内。土质肥沃,属于玄武岩中低山灰化土灰棕壤类型。属温带大陆性季风气候,降水充足,年均气温一般不超过4.0 ℃。区域的极端温差在− 38.0 ~ 38.0 ℃区间内。该区植物种类较多,森林资源丰富。以针阔混交林为主,主要树种包括樟子松(Pinus sylvestris var. mongolica)、红豆杉(Taxus wallichiana var. chinensis)、云杉(Picea asperata)、冷杉(Abies fabri)、红松(Pinus koraiensis)、长白落叶松(Larix olgensis)等。

    • 研究数据来源于吉林省汪清林业局复位调查样地数据。每块固定样地的面积为0.06 hm2,共207块,对其中的异常数据进行剔除处理,确定出满足要求的标准地有202块。林分调查中选择的林分因子主要包括优势树种、平均高、平均年龄等测树因子以及坡度、坡向、海拔等立地因子。对所得结果进行统计处理时,选择了间距10年的两次测量值作为一个独立记录,确定出13 897个记录。选择其中的80%的数据进行建模,其余20%作为测试集进行模型检验。表1是样地调查数据的统计信息。

      表 1  建模数据和检验数据样地基本因子统计量

      Table 1.  Summary statistics of sample plots for calibration and validation data

      项目
      Item
      株数
      Number of trees
      自变量
      Independent variable
      平均值
      Average
      最大值
      Maximum
      最小值
      Minimum
      标准差
      Standard deviation
      建模数据 Calibration data 11 118 D 12.386 5 86.806 5 5.000 0 7.456 7
      BAL 0.846 3 2.040 9 0.000 0 0.397 1
      N 1 599 3 950 467 675
      Nm 135 733 17 134
      检验数据 Validation data 2 780 D 12.155 9 73.690 9 5.003 6 7.440 2
      BAL 0.841 0 1.881 1 0.000 0 0.390 6
      N 1 482 3 950 217 698
      Nm 202 700 17 135
      注:D为林木胸径(cm),BAL为林木胸高断面积(m2/hm2),N为林分密度(株/hm2),Nm为林分枯死株数(株/hm2)。Notes: D represents DBH (cm); BAL represents basal area of breast height of forest trees (m2/ha); N represents stand density (tree/ha); Nm represents number of dead trees (tree/ha).
    • 由于本次研究使用的数据与马武[16]的数据相同,因此基础模型直接选择了马武[16]建立的模型。具体表达式为:

      $$ \left\{ \begin{array}{l} \!\!\!\! P = \dfrac{{{{\rm{e}}^{{a_0} + {a_1}D + {a_2}{\rm{BAL}} + {a_3}{\rm{BAP}} + {a_4}{\rm{ASE}}}}}}{{1 + {{\rm{e}}^{{a_0} + {a_1}D + {a_2}{\rm{BAL}} + {a_3}{\rm{BAP}} + {a_4}{\rm{ASE}}}}}} + \varepsilon \\ \!\!\!\!{\rm{BAP}}\;{\rm{ = }}\;\left( {g/G} \right)\; \times 10\;000\\ \!\!\!\!{\rm{ASE}}\;{\rm{ = }}\;{\cos}\;{A} \times \left( {\ln \;E/1\;000} \right)\; \times \;\tan \;S \end{array} \right. $$ (1)

      式中:P为林木10年间的枯死概率,只有0和1两个值,P = 1时,代表林木枯死,P = 0时,代表林木存活;D、BAL、BAP、ASE为自变量,其中,D 为林木胸径,BAL为林木胸高断面积,BAP是单木胸高断面积与样地胸高断面积的比,ASE是样地坡向、坡度和海拔的综合影响因子;g 为林木胸高断面积; G 为样地胸高断面积;A 为样地坡向,正北方为 0° ,正西方为 90° ,依次类推; E 为样地海拔; S 为样地坡度;a0a1a2a3a4为待估参数,$\varepsilon $是误差项。

    • 在公式(1)的基础上,设置样地为随机效应因子。以AIC值(赤池信息准则)越小越好的原则进行筛选处理,从而确定出最合适的混合效应模型。对应的模型表达式如下:

      $$ \left\{ {\begin{array}{*{20}{l}} \!\!\!\!\!{{{P_{ij}} = \dfrac{{{{\rm{e}}^{\left( {\left( {{a_{0i}} + a{r_i}} \right) + {a_1}{d_{ij}} + {a_2}{\rm{BA}}{{\rm{L}}_{ij}} + {a_3}{\rm{BA}}{{\rm{P}}_{ij}} + {a_4}{\rm{AS}}{{\rm{E}}_i}} \right)}}}}{{1 + {{\rm{e}}^{\left( {\left( {{a_{0i}} + a{r_i}} \right) + {a_1}{d_{ij}} + {a_2}{\rm{BA}}{{\rm{L}}_{ij}} + {a_3}{\rm{BA}}{{\rm{P}}_{ij}} + {a_4}{\rm{AS}}{{\rm{E}}_i}} \right)}}}}}}\\ \!\!\!\!\!{{{\beta _i} \! = \! \left( {{a_{0i}},{a_{1i}},{a_{2i}},{a_{4i}}} \right) = f\left( {\lambda ,{b_i}} \right) = \left[ {\begin{array}{*{20}{c}} \!\!\!\!{{a_0} + {a_{r0i}}}\!\!\!\!\!\!\\ \!\!\!\!{{a_1} + 0}\!\!\!\!\\ \!\!\!\!{{a_2} + 0}\!\!\!\!\\ \!\!\!\!{{a_3} + 0}\!\!\!\!\\ \!\!\!\!{{a_4} + 0}\!\!\!\! \end{array}} \right], {\text{为参数向量}}}}\\ \!\!\!\!\!{{\lambda = \left( {{a_{0i}},{a_{1i}},{a_{2i}},{a_{3i}},{a_{4i}}} \right),{\text{是固定效应向量}}}}\\ \!\!\!\!\!{{{b_i} = {a_{r0i}},{\text{且}}{b_i} \sim N\left( {0,\sigma _{{a_0}}^2} \right),{\text{是随机效应向量}}}}\\ \!\!\!\!\!{{{\varepsilon _{ij}} \sim N\left( {0,{\sigma ^2}} \right),{\text{是残差向量}}}} \end{array}} \right. $$ (2)

      式中:Pij是第i块样地第j株林木的枯死概率,dij是第i块样地第j株林木的胸径信息,BALij是第i块样地中大于对象木j的所有林木胸高断面积和;ASEi具体表示相应的第i块样地坡向综合影响因子;a0a1a2a3a4是固定效应参数,ar0i是随机效应参数,$\varepsilon_{ij} $是误差项。

    • 贝叶斯法是将参数θ视作连续变量,根据参数的先验信息与样本信息由全概率公式计算后验概率得到参数θ的条件概率分布的统计方法。即令模型形式为:

      $$ y = f\left( {x|\theta } \right) $$ (3)

      则根据先验信息πθ)、样本信息py|θ),由全概率公式得到给定样本数据y的参数θ的条件分布πθ|y):

      $$ \pi (\theta |y) = \frac{{p(y|\theta )\pi \left( \theta \right)}}{{p\left( y \right)}} $$ (4)

      本研究中,使用MCMC算法中的GIBBS抽样法估算参数。

      相较于传统统计方法,贝叶斯法增加了对先验信息的考虑,把样本与参数都视为随机变量,不要求服从正态分布,尤其在小样本量的估算时,贝叶斯法估计的参数的可信区间相较于传统极大似然估计的置信区间更小,模型精度也更高。

      本研究先基于公式(1)使用传统的极大似然法得到模型参数的估计结果,并将该结果作为贝叶斯法的估计参数的先验信息,结合样本信息使用贝叶斯法估算参数的后验概率分布,得到模型的估计结果。

    • 混合模型近年来在林业领域的应用越来越广[17-18],其特征表现为存在显著清晰的层次结构,包含了固定效应与随机效应,可代表总体到个体的多个层次。这种模型可对总体的一般变化趋势进行描述分析,同时也可以反映出个体间的差异。混合效应模型在贝叶斯框架中,一般基于层次贝叶斯法进行参数估计。混合效应模型在参数估计时传统的可选择的方法一般为最大似然或最小二乘方法,而目前的研究中使用层次贝叶斯法进行参数估计也渐成趋势。

      层次贝叶斯法较之贝叶斯法的区别就是在先验参数a中设置一个层次。

    • 贝叶斯统计进行分析过程中,需要对这些模型选择满足要求的先验信息,与此相关的信息主要来源于最大似然所得结果,一般情况下都选择正态分布:$\lambda \sim {{N}}\left( {{\lambda ^*},M} \right)$$\lambda ^*$是相关的估计值,可以据此进行统计回归确定出目标参数。

      层次贝叶斯统计法进行分布分析时,可选择无先验信息条件下的分布进行处理。固定λ条件下,其对应的分布可满足正态分布$\lambda \sim {{N}}( {0,{{1\; 000}^2}} )$,同时选择$\sigma _{{a_{_0}}}^{ - 1}$服从Gamma分布。表2显示了对应的先验信息相关结果。

      表 2  参数的先验信息

      Table 2.  Prior distribution of each parameter for base model and mixed-effects model

      模型 Model参数 Parameter先验分布 Prior distribution
      基础模型
      Base model
      a0 a0 ~ N (−2.148 3, 0.149 52)
      a1 a1 ~ N (−0.126, 0.010 52)
      a2 a2 ~ N (1.294 3, 0.093 92)
      a3 a3 ~ N (0.176 3, 0.007 32)
      a4 a4 ~ N (−1.893 9, 0.444 82)
      混合模型
      Model with
      random effects
      a0 a0 ~ N (0, 10 002)
      a1 a1 ~ N (0, 10 002)
      a2 a2 ~ N (0, 10 002)
      a3 a3 ~ N (0, 10 002)
      a4 a4 ~ N (0, 10 002)
      σa0 (1/σa0)2 ~ Gamma (0.001, 0.001)
      σ (1/σ)2 ~ Gamma (0.001, 0.001)
    • 二分类逻辑回归模型需要设定临界概率以确定目标值为1还是0,即该株树是否枯死。本研究临界概率的确定依据MST原则进行。MST原则确定临界概率的具体流程:(1)首先分别从0到1以0.01为间隔,确定出100个不同的临界概率的敏感度和特异度(见表3),其中敏感度是正确预测林木的枯死数与枯死观测值的比值,特异度是正确预测树木存活株数与存活观测株数的比值;(2)以敏感度与特异度相加之和所得的最大概率阈值为临界概率,本文通过分析计算确定的对应临界概率为0.22。

      表 3  观测值和预测值的列联表

      Table 3.  Contingency table of the observed values and predicted one

      项目 Item预测值 Predicted value总和 Total
      枯死 Dead存活 Alive
      观测值 Observed value 枯死 Dead a11 a12 A1m
      存活 Alive a21 a22 A2m
      和 Total An1 An2 Anm
      敏感度 Sensitivity=a11/A1m, 特异度 Specificity=a22/A2m
      注:a11表示正确预测枯死林木的株数,a12表示错误预测枯死林木的株数,a21表示错误预测林木存活的株数,a22表示正确预测存活林木的株数,A1m表示枯死实测值,A2m表示存活实测值,An1表示枯死预测值,An2表示存活预测值。Notes: a11 represents the number of dead trees which is correctly predicted, a12 represents the number of dead trees which is not correctly predicted, a21 represents the number of alive trees which is not correctly predicted, a22 represents the number of alive trees which is correctly predicted, A1m represents the measured value of dead trees, A2m represents the measured value of alive trees, An1 represents the predicted value of dead trees, An2 represents the predicted value of alive trees.
    • 模型的精度检验主要从模型拟合效果和拟合优度检验两个方面进行。

      本研究使用ROC曲线下面积AUC值与Pearson-χ2检验来判断模型的拟合效果与拟合优度。ROC曲线是以不同阈值下特异度与敏感度所绘制的曲线图,其曲线下面积AUC值越大越好,一般来说,AUC值为0.5 ~ 0.7时模型判断准确度较低,0.7 ~ 0.9时判断准确度较高,0.9以上则为优秀水平。

      由于本研究使用的基础模型是二分类Logistic模型,因此选用Pearson-χ2检验来检测模型预测的枯死树木株数与样地实测枯死株数有无显著差异。其统计量公式如下:

      $$ {\chi ^2} = \frac{{{{({\rm{N}}{{\rm{d}}_{{\rm{obs}}}} - {\rm{N}}{{\rm{d}}_{{\rm{pred}}}})}^2}}}{{{\rm{N}}{{\rm{d}}_{{\rm{pred}}}}}} $$ (5)

      式中:χ2是自由度为nχ2统计值,Ndobs和Ndpred是不同分组(树种、胸径等)的死亡林木株数的实测值和预测值。

    • 拟合结果如表4所示。对比分析结果可知贝叶斯法与传统极大似然法两种方法拟合的参数估计值相近,但贝叶斯统计标准差较小,而层次贝叶斯法的估计结果与前两种方法略有不同,且估计参数的标准差均较大。

      表 4  基于建模数据的参数估计结果

      Table 4.  Parameter estimates of classical and Bayesian methods using calibration

      方法
      Method
      参数
      Parameter
      估计值
      Estimated value
      标准差
      Standard deviation
      置信区间 Credible interval
      2.50%97.50%
      最大似然 Maximum likelihood a0 −2.148 3 0.149 5 −2.447 2 −1.849 4
      a1 −0.126 0 0.010 5 −0.147 1 −0.104 9
      a2 1.294 3 0.093 9 1.106 6 1.482 1
      a3 0.176 3 0.073 2 0.029 8 0.322 8
      a4 −1.893 9 0.444 8 −2.783 5 −1.004 3
      贝叶斯统计 Bayesian method a0 −2.136 0 0.138 8 −2.413 6 −1.858 4
      a1 −0.123 1 0.007 2 −0.137 5 −0.108 7
      a2 1.271 0 0.078 8 1.113 3 1.428 7
      a3 0.136 1 0.062 3 0.011 6 0.260 6
      a4 −1.886 0 0.415 4 −2.716 8 −1.055 2
      层次贝叶斯 Hierarchical bayesian method a0 −2.641 0 0.289 7 −3.248 0 −2.176 0
      a1 −0.145 4 0.013 0 −0.168 6 −0.118 0
      a2 1.595 0 0.176 0 1.311 0 1.963 0
      a3 0.274 8 0.084 9 0.108 7 0.444 0
      a4 −3.082 0 1.337 0 −5.523 0 −0.452 2
      σa0 1.030 0 0.158 9 0.742 8 1.360 0

      分析图1所得结果可知:这3种方法置信区间存在明显的重合,其中贝叶斯法估计的模型参数的可信区间最小,明显小于最大似然法估计参数的置信区间;层次贝叶斯统计所得结果的置信区间较为分散,大于最大似然的置信区间,最大可达到185.00%。

      图  1  单木枯死模型参数的贝叶斯估计的95%可信区间

      Figure 1.  Bayesian 95% credible and classical 95% confidence intervals for parameters in mortality mode

      图2给出了用3种方法确定出的参数分布情况。对所得结果进行对比分析可知:贝叶斯法(B)统计对应的分布最为集中,且与传统极大似然法(C)估计出的参数分布出现大面积的重合;层次贝叶斯法(HB)统计的参数分布最为分散,这是因为随机效应的加入,增加了参数的不稳定性,包含了更多的样地差异信息。

      图  2  枯死模型的贝叶斯、层次贝叶斯参数的后验分布和极大似然估计参数的似然分布

      Figure 2.  Posterior distribution of Bayesian method and likelihood distribution of maximum likelihood method for morality model

    • 3种方法基于建模集与验证集的模型检验结果如表5所示。

      表 5  3种方法拟合的枯死模型误差统计量

      Table 5.  Error statistics of morality model by three methods

      方法 Method建模数据 Calibration data验证数据 Validation data
      枯死数
      Number of
      dead trees
      预测枯死数
      Predicted number
      of dead trees
      χ2AUCAIC (DIC)枯死数
      Number of
      dead trees
      预测枯死数
      Predicted number
      of dead trees
      χ2AUC (DIC)
      传统方法
      Classical method
      1 338 1 349 0.09 0.734 7 407.4 375 404 2.08 0.726
      贝叶斯统计
      Bayesian method
      1 338 1 336 0.01 0.734 7 407.21 375 394 0.92 0.727
      层次贝叶斯统计
      Hierarchical Bayesian method
      1 338 1 335 0.01 0.825 6 771.26 375 389 0.50 0.800

      表5可知:3种方法基于建模数据的AUC值在0.734 ~ 0.825范围内,模型的诊断表现良好,其中,极大似然估计和贝叶斯估计的模型AUC值均为0.734;层次贝叶斯法估计的模型AUC值为0.825,层次贝叶斯法精度最高。对应的ROC曲线如图3所示。在AIC(DIC)表现上,层次贝叶斯法估计的模型表现同样最好,DIC最小,为6 771.26。3种方法基于验证集的结果显示:层次贝叶斯法同样表现最好,AUC值最大,为0.800;贝叶斯法的AUC与传统极大似然估计的AUC相近,分别为0.727与0.726。

      图  3  3种方法拟合枯死模型的ROC曲线

      Figure 3.  ROC curves of mortality model based on three methods

      χ2检验的结果表现上,除传统极大似然方法估计的模型χ2值为0.09以外,贝叶斯法与层次贝叶斯法估计的模型χ2值均为0.01,说明3种方法的模型预测枯死林木株数与样地实测枯死株数并无差异,均通过χ2检验,具有统计意义。

    • 本研究建立了基于贝叶斯方法的单木枯死模型,发现贝叶斯统计的参数后验分布与传统极大似然法的参数似然分布的重合度高,且更为集中,说明贝叶斯法估计参数的后验结果会受到先验信息不小的影响。张雄清等[13]对树高曲线的研究也得出了类似的研究结果。由于考虑了随机效应,层次贝叶斯法的参数后验分布对比贝叶斯法的参数后验分布则较为扩散,这与Patenaud等[19]使用层次贝叶斯统计在生物学区域模型中的研究结果类似。

      对比枯损模型中使用贝叶斯法估计参数的标准差和使用传统极大似然估计的参数标准差,在枯损模型中,有先验信息贝叶斯统计的参数标准差比最大似然法都要低,范围在6.00% ~ 31.74%。这可能是由于贝叶斯法本身就引入了来自于传统极大似然法的估计结果作为先验信息,使参数估计的结果更加稳定[20]。层次贝叶斯法估计参数的标准差则远大于其余两种方法估计参数的标准差,Russell[21]指出这是由于随机效应的加入,模型解释了更多的样地间差异。

      从模型的拟合效果来看,3种方法均通过χ2检验,但其中表现最好的是层次贝叶斯法的模型结果,其在检验样本中表现的AUC值为0.800,高于贝叶斯法和传统极大似然法估计的0.727和0.726,AIC(DIC)指标也说明了相同的结论。这与Russell[21]在枯损模型中得出的结论类似,即加入了随机效应的层次贝叶斯模型的拟合精度高于只含有固定效应的贝叶斯模型。因此,层次贝叶斯模型在建模过程中既考虑到了先验信息,又包含有随机效应,在估计单木枯死模型参数上有较大优势。

参考文献 (21)

目录

    /

    返回文章
    返回