高级检索

留言板

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

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

塞罕坝华北落叶松人工林断面积预测模型

王冬至 张冬燕 张志东 黄选瑞

王冬至, 张冬燕, 张志东, 黄选瑞. 塞罕坝华北落叶松人工林断面积预测模型[J]. 北京林业大学学报, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072
引用本文: 王冬至, 张冬燕, 张志东, 黄选瑞. 塞罕坝华北落叶松人工林断面积预测模型[J]. 北京林业大学学报, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072
WANG Dong-zhi, ZHANG Dong-yan, ZHANG Zhi-dong, HUANG Xuan-rui. Prediction model for basal area of Larix principis-rupprechtii plantation in Saihanba of Hebei Province, northern China[J]. Journal of Beijing Forestry University, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072
Citation: WANG Dong-zhi, ZHANG Dong-yan, ZHANG Zhi-dong, HUANG Xuan-rui. Prediction model for basal area of Larix principis-rupprechtii plantation in Saihanba of Hebei Province, northern China[J]. Journal of Beijing Forestry University, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072

塞罕坝华北落叶松人工林断面积预测模型

doi: 10.13332/j.1000-1522.20170072
基金项目: 

“十二五”国家科技支撑计划项目 2015BAD09B01

林业公益性行业科研专项 20150430304

国家自然科学基金项目 31370636

详细信息
    作者简介:

    王冬至,博士,讲师。主要研究方向:森林可持续经营。Email: wangdz@126.com  地址: 071000  河北省保定市南市区乐凯南大街2596号河北农业大学西校区林学院

    通讯作者:

    黄选瑞,教授,博士生导师。主要研究方向:森林可持续经营。Email: hxr1962@163.com  地址:同上

  • 中图分类号: S758.5+7

Prediction model for basal area of Larix principis-rupprechtii plantation in Saihanba of Hebei Province, northern China

  • 摘要: 如何实现林分水平和单木水平预测林分断面积的兼容性,并提高不同水平断面积模型的预估精度,在森林经营过程中是一个亟待解决的科学问题。本文以华北暖温带华北落叶松人工林为研究对象,运用105块连续观测的固定样地数据,首先采用Gauss-Newton算法建立林分水平和单木水平断面积生长模型;其次分别采用不同形式的逻辑斯蒂方程对单木生存概率方程进行拟合;最后将不同水平最优预测模型进行组合并建立组合方程,采用最小二乘法估计组合方程参数,以提高对不同水平断面积的预测精度。结果表明:在约束参数法中,林分密度模型、林分断面积预测模型和单木断面积模型均具有较好的预测效果,并均能解释90%以上的变异;在分解法中采用逻辑斯蒂方程来预测单木生存概率和林分密度,经检验获得的ROC曲线下面积为0.906,表明该方程可以较好地预测林木生存概率;在组合预测法中,采用不同水平的最优模型进行组合后的预测效果最佳。在预测林分密度和断面积时,组合预测方程预测精度最高,林分水平模型预测精度次之,单木水平模型预测精度最低。组合预测法能够预测不同水平下的林分密度、立木生存概率、林分断面积及单木断面积,提高了模型预测精度,为预测林分生长动态、空间结构变化及经营效果评价等提供参考依据。
  • 图  1  不同水平最优模型残差分布

    M4、M7、M9分别代表林分密度、林分断面积和单木断面积预测模型,即对应文中的公式(4)、(7)和(9)。

    Figure  1.  Residual distribution of different level optimal models

    M4, M7, M9 represent prediction models of stand density, basal area for stand, basal area models for individual trees, which corresponding to formula (4), (7) and (9) of the paper, respectively.

    图  2  单木生存概率ROC曲线

    Figure  2.  ROC curve for the individual tree mortality model fitted

    表  1  建模数据与检验数据统计

    Table  1.   Summary statistics for modeling and validation data sets

    数据
    Data
    变量
    Variable
    2011年调查Inventory in 20112016年调查Inventory in 2016
    最小值
    Min.
    最大值
    Max.
    均值
    Mean
    标准差
    SD
    最小值
    Min.
    最大值
    Max.
    均值
    Mean
    标准差
    SD
    建模数据Modelling data(n=70)林分年龄Stand age134231.510.22184737.09.21
    胸径Diameter at breast height(DBH)/cm6.2728.0217.175.619.0631.5219.456.13
    平均高Mean height/m4.7116.3411.042.336.5419.0212.842.97
    林分密度/(株·hm-2)Stand density/(tree·ha-1)900.03 353.01 927.7958.4500.02 282.01 295.6611.2
    林分断面积/(m2·hm-2)Stand basal area/(m2·ha-1)8.8267.4336.4512.0711.5190.3142.7117.38
    立地质量Site index/m4.6616.8112.362.917.0618.5614.143.43
    生存概率Survival probability/%0.351.000.770.17
    检验数据Validation data(n=35)林分年龄Stand age134033.510.43184738.09.46
    胸径Diameter at breast height(DBH)/cm6.8729.1117.585.079.1232.4220.245.97
    平均高Mean height/m5.1117.2311.802.047.0518.9213.273.09
    林分密度/(株·hm-2)Stand density/(tree·ha-1)900.03 178.01 802.7907.8500.02 077.01 226.7578.4
    林分断面积/(m2·hm-2)Stand basal area/(m2·ha-1)9.0270.5837.0912.5410.9785.6645.0715.06
    立地质量Site index/m5.2416.3313.113.136.9819.2415.373.55
    生存概率Survival probability/%0.411.000.740.19
    注:n为样本数。Note:n is sample number.
    下载: 导出CSV

    表  2  不同水平模型参数估计与统计检验

    Table  2.   Parameter estimating and goodness-of-fit statistics in different levels

    模型Model参数Parameter估计值EstimateSER2BiasRMSE
    M4a10.005 10.003 50.900 21.952 52.186 3
    a20.406 50.399 5
    a31.271 30.526 7
    M5a1-0.405 80.095 40.804 32.005 62.334 2
    a2-0.738 70.415 8
    a3-0.045 20.009 0
    M6b1-0.000 80.000 30.812 31.280 93.155 0
    b27.143 82.056 7
    b33.222 40.684 5
    M7b124.978 61.529 10.914 40.899 71.833 4
    b21.215 30.191 7
    b30.799 20.040 5
    M8b10.999 30.000 80.867 11.008 42.304 1
    b20.000 50.000 3
    b3-2.853 60.677 6
    M9c00.011 20.002 30.913 20.097 50.918 7
    c1-4.284 30.163 3
    c24.396 30.151 7
    c3-24.096 12.385 7
    M10c00.000 20.000 10.908 80.105 81.007 6
    c10.756 70.111 3
    c2-0.528 10.262 4
    c31.093 00.081 3
    c4-0.042 20.012 1
    注:M4~M10对应文中公式(4)~(10)。Note: M4-M10 correspond to formula (4)-(10) of the paper, respectively.
    下载: 导出CSV

    表  3  不同水平生存概率模型参数估计与统计检验

    Table  3.   Parameter estimating and goodness-of-fit statistics in different levels of survival probability

    模型Model参数Parameter估计值EstimateSER2BiasRMSE
    M11f0-3.294 00.780 40.910 41.022 42.077 1
    f10.000 50.000 3
    f20.009 20.002 9
    f324.111 617.138 5
    M12f0-3.663 00.616 00.904 61.075 82.332 4
    f111.015 14.372 4
    f20.036 00.007 8
    f30.007 90.002 3
    注:M11、M12对应文中公式(11)、(12)。Note: M11,M12 correspond to formula (11),(12) of the paper, respectively.
    下载: 导出CSV

    表  4  林分、单株生长模型与组合方程预测林分变量评价

    Table  4.   Evaluation statistics for predicted stand variables from the stand and individual-tree growth models and forecast combination

    变量
    Variable
    每公顷株数预测Predicted tree number per ha林分断面积预测Predicted stand basal area
    RMSERRMSEMADRMSERRMSEMAD
    林分生长Stand growth305.739 47.151.933 82.054 78.440.880 1
    单木生长Individual tree growth364.008 713.152.406 12.931 615.660.973 4
    组合预测Combination forecast284.446 35.421.603 21.664 2-5.610.791 6
    下载: 导出CSV
  • [1] ZHANG S, AMATEIS R L, BURKHART H E. Constraining individual tree diameter increment and survival models for loblolly pine plantations[J]. Forest Science, 1997, 43(6): 414-423. http://europepmc.org/abstract/AGR/IND21235335
    [2] CIESZEWSKI C J. Comparing fixed-and variable-base-age site equations having single versus multiple asymptotes[J]. Forest Science, 2002, 48(1): 7-23. http://d.old.wanfangdata.com.cn/NSTLQK/NSTL_QKJJ027400739/
    [3] TEWARI V P, ÁLVAREZ-GONZALEZ J G, GARCÍA O. Developing a dynamic growth model for teak plantations in India[J]. Forest Ecosystems, 2014, 1(1): 9-17. doi:  10.1186/2197-5620-1-9
    [4] RITCHIE M W, HANN D W. Implications of disaggregation in forest growth and yield modeling[J]. Forest Science, 1997, 43 (2): 223-233. http://europepmc.org/abstract/AGR/IND21235281
    [5] QIN J, CAO Q V. Using disaggregation to link individual-tree and whole-stand growth models[J]. Canadian Journal of Forest Research, 2006, 36: 953-960. doi:  10.1139/x05-284
    [6] GONZALEZ J G, ZINGG A, GADOW K V. Estimating growth in beech forests: a study based on longterm experiments in Switzerland[J]. Annals of Forest Science, 2009, 307(9): 1-13. doi:  10.1051%2Fforest%2F2009113
    [7] RITCHIE M W, HANN D W. Implications of disaggregation in forest growth and yield modeling[J]. Forest Science, 1997, 43 (2): 223-233. http://europepmc.org/abstract/AGR/IND21235281
    [8] BURKHART H E, TOME M. Modeling forest trees and stands[M]. Berlin: Springer, 2012.
    [9] 陈清, 张令峰, 傅松玲.树木年龄和断面积对加拿大北方林树木死亡率的影响[J].应用生态学报, 2011, 22(9): 2477-2481. http://d.old.wanfangdata.com.cn/Periodical/yystxb201109038

    CHEN Q, ZHANG L F, FU S L. Effects of tree age and basal area on boreal forest tree mortality in Canada[J]. Chinese Journal of Applied Ecology, 2011, 22(9): 2477-2481. http://d.old.wanfangdata.com.cn/Periodical/yystxb201109038
    [10] 雷相东, 李永慈, 向玮.基于混合模型的单木断面积生长模型[J].林业科学, 2009, 45(1): 74-80. doi:  10.3321/j.issn:1001-7488.2009.01.014

    LEI X D, LI Y C, XIANG W. Individual basal area growth model using multi-level linear mixed model with repeated measures[J]. Scientia Silvae Sinicae, 2009, 45(1): 74-80. doi:  10.3321/j.issn:1001-7488.2009.01.014
    [11] 符利勇, 唐守正, 张会儒, 等.基于多水平非线性混合效应蒙古栎林单木断面积模型[J].林业科学研究, 2015, 28(1): 23-31. http://d.old.wanfangdata.com.cn/Periodical/lykxyj201501004

    FU L Y, TANG S Z, ZHANG H R, et al. Multilevel nonlinear mixed-effects basal area models for individual trees of Quercus mongolica[J]. Forest Research, 2015, 28(1): 23-31. http://d.old.wanfangdata.com.cn/Periodical/lykxyj201501004
    [12] 倪成才, 王庆丰.火炬松人工林胸高断面积差分模型的拟合与筛选[J].北京林业大学学报, 2011, 33(3): 1-7. doi:  10.3969/j.issn.1671-6116.2011.03.001

    NI C C, WANG Q F. Model selection and fit of algebraic difference models for basal area of loblolly pine plantations[J]. Journal of Beijing Forestry University, 2011, 33(3): 1-7. doi:  10.3969/j.issn.1671-6116.2011.03.001
    [13] 李春明, 唐守正.基于非线性混合模型的落叶松云冷杉林分断面积模型[J].林业科学, 2010, 46(7): 106-113. http://d.old.wanfangdata.com.cn/Periodical/lykx201007016

    LI C M, TANG S Z. The basal area model of mixed stands of Larix olgensis, Abies nephrolepis and Picea jezoensis based on nonlinear mixed model[J]. Scientia Silvae Sinicae, 2010, 46(7): 106-113. http://d.old.wanfangdata.com.cn/Periodical/lykx201007016
    [14] 张雄清, 张建国, 段爱国.杉木人工林林分断面积生长模型的贝叶斯法估计[J].林业科学研究, 2015, 28(4): 538-542. doi:  10.3969/j.issn.1001-1498.2015.04.013

    ZHANG X Q, ZHANG J G, DUAN A G. Application of bayesian method in stand basal area prediction of Chinese fir plantation[J]. Forest Research, 2015, 28(4): 538-542. doi:  10.3969/j.issn.1001-1498.2015.04.013
    [15] 张雄清, 雷渊才, 陈新美.林分断面积组合预测模型权重确定的比较[J].林业科学, 2011, 47(7): 36-41. http://d.old.wanfangdata.com.cn/Periodical/lykx201107006

    ZHANG X Q, LEI Y C, CHEN X M. Comparison of weight computation in stand basal area combined model[J]. Scientia Silvae Sinicae, 2011, 47(7): 36-41. http://d.old.wanfangdata.com.cn/Periodical/lykx201107006
    [16] VALBUENA P, DELPESO C, BRAVO F. Stand density management diagrams for two mediterranean pine species in eastern spain[J]. Investigación Agraria: Sistemas Recursos Forestales, 2008, 17(2): 97-104. doi:  10.5424/srf/2008172-01026
    [17] ZHANG X, LEI Y, CAO Q V, et al. Improving tree survival prediction with forecast combination and disaggregation[J]. Canadian Journal of Forest Research, 2011, 41: 1928-1935. doi:  10.1139/x11-109
    [18] ANDREA H, CAO Q V, ALVAREZ J G, et al. Compatibility of whole-stand and individual-tree models using composite estimators and disaggregation[J]. Forest Ecology and Management, 2015, 348(11): 46-56 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=38160e55b539dc3a70a82cf617f9ea25
    [19] ARANDA D U, GRANDAS J A, ALVAREZ J G, et al. Site quality curves for birch stands in north-western Spain[J]. Silva Fennica, 2006, 40 (4): 631-644. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=e6b7779e4c77cdd32c556b781c213438
    [20] 王冬至, 张冬燕, 王方, 等.塞罕坝主要立地类型针阔混交林树高曲线构建[J].北京林业大学学报, 2016, 38(10): 7-14. doi:  10.13332/j.1000-1522.20150359

    WANG D Z, ZHANG D Y, WANG F, et al. Height curve construction of needle and broadleaved mixed forest under main site types in Saihanba, Hebei of northern China[J]. Journal of Beijing Forestry University, 2016, 38(10): 7-14. doi:  10.13332/j.1000-1522.20150359
    [21] GARCIA O. A parsimonious dynamic stand model for interior spruce in British Columbia[J]. Forest Science, 2011, 57 (4): 265-280. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=47de0ae6bf440ef1a8dd3cec62145e98
    [22] GARCIA O. Building a dynamic growth model for trembling aspen in western Canada without age data[J]. Canadian Journal of Forest Research, 2013, 43 (3): 256-265. doi:  10.1139/cjfr-2012-0366
    [23] GARCIA O, BURKHART H E, AMATEIS R L. A biologically-consistent stand growth model for loblolly pine in the Piedmont physiographic region, USA[J]. Forest Ecology and Management, 2011, 262 (11): 2035-2041. doi:  10.1016/j.foreco.2011.08.047
    [24] CAO Q V. Linking individual-tree and whole-stand models for forest growth and yield prediction[J]. Forest Ecosystems, 2014, 1: 1-8. doi:  10.1186/2197-5620-1-1
    [25] JUMA R, PUKKALA T, DEMIGUEL S, et al. Evaluation of different approaches to individual tree growth and survival modelling using data collected at irregular intervals-a case study for Pinus patula in Kenya[J]. Forest Ecosystems, 2014, 8(1): 1-14. doi:  10.1186/s40663-014-0014-3
    [26] BATES J M, GRANGER C W J. The combination of forecasts[J]. A Quarterly Journal of Operations Research, 1969, 20 (4): 451-468. doi:  10.1057/jors.1969.103
    [27] VANCLAY J K. Modelling forest growth and yield: application to mixed tropical forests[M]. Wallingford: CAB International, 1994.
    [28] ZHANG X, LEI Y. A linkage among whole-stand model, individual-tree model and diameter-distribution model[J]. Journal of Forest Science, 2010, 56: 600-608. doi:  10.17221/102/2009-JFS
    [29] CRECENTE-CAMPO F, SOARES P, TOME M, et al. Modelling annual individual-tree growth and mortality of Scots pine with data obtained at irregular measurement intervals and containing missing observations[J]. Forest Ecology and Management, 2010, 260: 1965-1974. doi:  10.1016/j.foreco.2010.08.044
    [30] CAO Q V. Prediction of annual diameter growth and survival for individual trees from periodic measurements[J]. Forest Science, 2000, 46: 127-131. http://europepmc.org/abstract/AGR/IND22301989
    [31] NORD-LARSEN T. Modeling individual-tree growth from data with highly irregular measurement intervals[J]. Forest Science, 2006, 52: 198-208.
    [32] CAO Q V, STRUB M. Evaluation of four methods to estimate parameters of an annual tree survival and diameter growth model[J]. Forest Science, 2008, 54 (6): 617-624. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=456a6d7c29b0fedca89786f7af26405a
    [33] WYKOFF W R. A basal area increment model for individual conifers in the northern Rocky Mountains[J]. Forest Science, 1990, 36 (4): 1077-1104. http://europepmc.org/abstract/AGR/IND91008565
    [34] MONSERUD R A, STERBA H. A basal area increment model for individual trees growing in even-and uneven-aged forest stands in Austria[J]. Forest Ecology and Management, 1996, 80 (3): 57-80. https://www.sciencedirect.com/science/article/pii/0378112795036385
  • [1] 张健飞, 王淳, 徐雯雯, 黄选瑞, 张志东.  华北落叶松不同代际人工林土壤养分及细菌群落变化特征 . 北京林业大学学报, 2020, 42(3): 36-45. doi: 10.12171/j.1000-1522.20190256
    [2] 马静怡, 黄华国, 黄侃, 邢路.  基于16线阵TLS数据的单木识别及林分断面积估测研究 . 北京林业大学学报, 2018, 40(8): 23-32. doi: 10.13332/j.1000-1522.20180016
    [3] 张文一, 景天忠, 严善春.  基于机器学习的落叶松毛虫发生面积预测模型 . 北京林业大学学报, 2017, 39(1): 85-93. doi: 10.13332/j.1000-1522.20160205
    [4] 赵天梁.  山西华北落叶松群落物种多样性 . 北京林业大学学报, 2017, 39(6): 45-50. doi: 10.13332/j.1000-1522.20170067
    [5] 张素芳, 张磊, 赵佳丽, 张莉, 张含国.  长白落叶松小RNA测序和其靶基因预测 . 北京林业大学学报, 2016, 38(12): 64-72. doi: 10.13332/j.1000-1522.20150404
    [6] 王冬至, 张冬燕, 王方, 张志东, 黄选瑞.  塞罕坝主要立地类型针阔混交林树高曲线构建 . 北京林业大学学报, 2016, 38(10): 7-14. doi: 10.13332/j.1000-1522.20150359
    [7] 刘宇, 郭建斌, 王彦辉, 刘泽彬, 邓秀秀, 张桐, 熊伟, 左海军.  宁夏六盘山不同密度华北落叶松人工林枯落物水文效应 . 北京林业大学学报, 2016, 38(8): 36-44. doi: 10.13332/j.1000-1522.20160007
    [8] 姚丹丹, 雷相东, 张则路.  基于贝叶斯法的长白落叶松林分优势高生长模型研究 . 北京林业大学学报, 2015, 37(3): 94-100. doi: 10.13332/j.1000-1522.20140221
    [9] 王冬至, 张志东, 牟洪香, 李永宁, 黄选瑞.  结构方程模型在落叶松林经营中的应用 . 北京林业大学学报, 2015, 37(3): 69-75. doi: 10.13332/j.1000-1522.20140326
    [10] 赵匡记, 王利东, 王立军, 贾忠奎, 马履一.  华北落叶松蓄积量及生产力研究 . 北京林业大学学报, 2015, 37(2): 24-31. doi: 10.13332/j.cnki.jbfu.2015.02.011
    [11] 李丹, 戴巍, 闫志刚, 王霓虹.  基于模糊层次分析法的落叶松人工林生境评价系统 . 北京林业大学学报, 2014, 36(4): 75-81. doi: 10.13332/j.cnki.jbfu.2014.04.015
    [12] 王巍伟, 吕瑞恒, 刘勇, 陈晓, 李国雷, .  不同氮含量华北落叶松叶凋落物在不同间伐强度林内氮释放规律研究 . 北京林业大学学报, 2014, 36(3): 63-68. doi: 10.13332/j.cnki.jbfu.2014.03.009
    [13] 倪成才, 王庆丰.  火炬松人工林胸高断面积差分模型的拟合与筛选 . 北京林业大学学报, 2011, 33(3): 1-7.
    [14] 李庆梅, 刘艳, 陈怀梁, 马风云, 侯龙鱼.  华北落叶松和日本落叶松种子超干保存活力和生理特性研究 . 北京林业大学学报, 2011, 33(2): 25-30.
    [15] 张雄清, 雷渊才, 陈新美, 王金增.  组合预测法在林分断面积生长预估中的应用 . 北京林业大学学报, 2010, 32(4): 6-11.
    [16] 李国雷, 刘勇, 吕瑞恒, 于海群, 李瑞生.  华北落叶松人工林密度调控对林下植被发育的作用过程 . 北京林业大学学报, 2009, 31(1): 19-24.
    [17] 罗云建, 张小全, 王效科, 朱建华, 张治军, 孙贵生, 高峰, .  华北落叶松人工林生物量及其分配模式 . 北京林业大学学报, 2009, 31(1): 13-18.
    [18] 刘勇, 李国雷, 林平, 姜辉, 于海群, 吕瑞恒.  华北落叶松人工幼、中龄林土壤肥力变化 . 北京林业大学学报, 2009, 31(3): 17-23.
    [19] 刘春延, 谷建才, 李吉跃, 陈平, 陆贵巧, 田国恒, .  塞罕坝华北落叶松生长与气候因子的相关分析 . 北京林业大学学报, 2009, 31(4): 102-105.
    [20] 岳永杰, 余新晓, 刘彦, 樊登星, 史宇, 王小平, 陈俊崎, .  华北落叶松林不同发育阶段种群分布格局研究 . 北京林业大学学报, 2008, 30(supp.2): 171-176.
  • 加载中
图(2) / 表 (4)
计量
  • 文章访问数:  558
  • HTML全文浏览量:  107
  • PDF下载量:  22
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-03-10
  • 修回日期:  2017-03-27
  • 刊出日期:  2017-07-01

塞罕坝华北落叶松人工林断面积预测模型

doi: 10.13332/j.1000-1522.20170072
    基金项目:

    “十二五”国家科技支撑计划项目 2015BAD09B01

    林业公益性行业科研专项 20150430304

    国家自然科学基金项目 31370636

    作者简介:

    王冬至,博士,讲师。主要研究方向:森林可持续经营。Email: wangdz@126.com  地址: 071000  河北省保定市南市区乐凯南大街2596号河北农业大学西校区林学院

    通讯作者: 黄选瑞,教授,博士生导师。主要研究方向:森林可持续经营。Email: hxr1962@163.com  地址:同上
  • 中图分类号: S758.5+7

摘要: 如何实现林分水平和单木水平预测林分断面积的兼容性,并提高不同水平断面积模型的预估精度,在森林经营过程中是一个亟待解决的科学问题。本文以华北暖温带华北落叶松人工林为研究对象,运用105块连续观测的固定样地数据,首先采用Gauss-Newton算法建立林分水平和单木水平断面积生长模型;其次分别采用不同形式的逻辑斯蒂方程对单木生存概率方程进行拟合;最后将不同水平最优预测模型进行组合并建立组合方程,采用最小二乘法估计组合方程参数,以提高对不同水平断面积的预测精度。结果表明:在约束参数法中,林分密度模型、林分断面积预测模型和单木断面积模型均具有较好的预测效果,并均能解释90%以上的变异;在分解法中采用逻辑斯蒂方程来预测单木生存概率和林分密度,经检验获得的ROC曲线下面积为0.906,表明该方程可以较好地预测林木生存概率;在组合预测法中,采用不同水平的最优模型进行组合后的预测效果最佳。在预测林分密度和断面积时,组合预测方程预测精度最高,林分水平模型预测精度次之,单木水平模型预测精度最低。组合预测法能够预测不同水平下的林分密度、立木生存概率、林分断面积及单木断面积,提高了模型预测精度,为预测林分生长动态、空间结构变化及经营效果评价等提供参考依据。

English Abstract

王冬至, 张冬燕, 张志东, 黄选瑞. 塞罕坝华北落叶松人工林断面积预测模型[J]. 北京林业大学学报, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072
引用本文: 王冬至, 张冬燕, 张志东, 黄选瑞. 塞罕坝华北落叶松人工林断面积预测模型[J]. 北京林业大学学报, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072
WANG Dong-zhi, ZHANG Dong-yan, ZHANG Zhi-dong, HUANG Xuan-rui. Prediction model for basal area of Larix principis-rupprechtii plantation in Saihanba of Hebei Province, northern China[J]. Journal of Beijing Forestry University, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072
Citation: WANG Dong-zhi, ZHANG Dong-yan, ZHANG Zhi-dong, HUANG Xuan-rui. Prediction model for basal area of Larix principis-rupprechtii plantation in Saihanba of Hebei Province, northern China[J]. Journal of Beijing Forestry University, 2017, 39(7): 10-17. doi: 10.13332/j.1000-1522.20170072
  • 单木水平和全林分水平模型是预测森林生长和产量预测模型中的重要组成部分,在森林可持续经营管理工作中具有重要作用[1]。单木模型以单株立木为基本单位,能够详尽描述立木生存和生长过程等信息[1-2];全林分模型是以林龄、立地及林分密度等因子模拟林分水平生长和收获的预测模型[3]。不同水平预测模型具有各自的优点与不足,全林分模型简单、粗放、便于应用,能够直接预测林分水平的变量,但无法反映单木水平的详细信息[4-6];而将单木水平模型应用到林分水平进行预测时,由于误差积累会产生较大的预测误差,因此预测结果精确度低。在模型应用过程中,单木生长模型在林分生长与产量的预测中往往会产生较大的预测误差,如果关注林分水平预测与经营,就会促使森林经营管理者采用全林分模型[5],而全林分模型在预测单木生长过程中又不可靠[7]。因此,建立既能获得良好林分水平预测结果,同时还能了解详细单木水平信息[8]的相容性模型,在森林经营管理中具有重要意义。

    目前在构建单木水平断面积预测模型时,主要采用最小二乘法、线性回归分析法[9]、混合模型法[10-11]等方法建立单木断面积生长模型;在构建林分水平断面积预测模型时,主要采用约束型极大似然估计(Restricted maximum likelihood,REML)法、极大似然估计(Maximum likelihood,ML)法[12]、混合模型法[13]及贝叶斯法[14]等方法来建立林分断面积生长预测模型;此外,在预测林分水平断面积过程中,张雄清等[15]基于误差平方和法、方差协方差法和最优加权法,建立了不同水平的断面积生长组合预测模型。由于单木水平变量能够预测林分水平变量,模型预测精度取决于林分水平模型的预测精度[16-17]。目前对于人工林单水平[10-11, 13](单木水平或林分水平)断面积的预测研究较多,而有关同时预测单木水平和林分水平断面积的研究较少,因此如何构建人工林不同水平的相容性断面积预测模型是一个亟待解决的科学问题。

    本研究以华北暖温带华北落叶松(Larix principis-rupprechtii)人工林为研究对象,对不同水平断面积预测模型通过约束参数法、分解法[18]和组合法3种方法,将单木水平和林分水平模型组合在一起。首先,采用约束参数法,建立林分水平和单木水平断面积预测模型,优化林分径级分布及单木水平预测;其次,在分解方法中,通过预测林分生长和死亡率来调整单木断面积生长模型,对林分断面积进行预测;最后将不同水平断面积预测模型进行组合,确定最优模型的组合并进行预测。

    • 塞罕坝机械林场总场(41°22′~42°58′ N,116°53′~118°31′ E)位于河北省最北部,由内蒙古克什克腾旗和多伦县的部分地区组成,属阴山山脉与大兴安岭余脉的交接地带。林区属华北暖温带立地类型区域,海拔1 010~1 940 m,年均气温-1.2 ℃,年均最高气温33.4 ℃,年均最低气温-43.3 ℃;年均日照时数2 548.7 h;年均降水量约452.2 mm。研究区主要乔木树种有华北落叶松(Larix principis-rupprechtii)、白桦(Betula platyphylla)、云杉(Picea asperata)、樟子松(Pinus sylvestris var. mongolica)、山杨(Populus davidiana)、油松(Pinus tabuliformis)、蒙古栎(Quercus mongolica)等,主要灌木树种有绣线菊(Spiraea salicifolia)、沙棘(Hippophae rhamnoides)、胡枝子(Lespedeza bicolor)、山刺玫(Rosa daverica)、稠李(Prunus padus)、华北忍冬(Lonicera tatarinowii)、大叶小檗(Berberis ferdinandi-coburgii)、栓翅卫矛(Euonymus phellomanus)等,主要草本植物有蒲公英(Herba taraxaci)、野蔷薇(Rosa multiflora)、曼陀罗(Datura stramonium)、苔草(Carex tristachya)、草地老鹳草(Geanium daharicum var. alpinum)、藜芦(Veratrum nigrum)、地榆(Sanguisorba officinalis)、唐松草(Thalictrum aquilegifolium)等.

    • 2011年7—8月份在河北省塞罕坝机械林场总场下属的北曼甸、大唤起、阴河、千层板、第三乡5个分林场,共设置华北落叶松人工林固定样地105块,样地大小均为30 m×30 m,对样地内胸径大于5 cm的立木进行编号,观测立木相对坐标并进行每木检尺,主要观测因子包括立地因子(海拔、坡度、坡位、坡向及土层厚度等)和林分因子(胸径、树高、冠幅、枝下高、密度及郁闭度等)。所有固定样地均在2011年11月份进行了传统经营(原有株数的30%),并于2016年7—8月份对所有固定样地进行了复测,复测因子与样地设置观测因子相同。研究中将固定样地连续观测数据分为建模数据和检验数据如表 1所示。

      表 1  建模数据与检验数据统计

      Table 1.  Summary statistics for modeling and validation data sets

      数据
      Data
      变量
      Variable
      2011年调查Inventory in 20112016年调查Inventory in 2016
      最小值
      Min.
      最大值
      Max.
      均值
      Mean
      标准差
      SD
      最小值
      Min.
      最大值
      Max.
      均值
      Mean
      标准差
      SD
      建模数据Modelling data(n=70)林分年龄Stand age134231.510.22184737.09.21
      胸径Diameter at breast height(DBH)/cm6.2728.0217.175.619.0631.5219.456.13
      平均高Mean height/m4.7116.3411.042.336.5419.0212.842.97
      林分密度/(株·hm-2)Stand density/(tree·ha-1)900.03 353.01 927.7958.4500.02 282.01 295.6611.2
      林分断面积/(m2·hm-2)Stand basal area/(m2·ha-1)8.8267.4336.4512.0711.5190.3142.7117.38
      立地质量Site index/m4.6616.8112.362.917.0618.5614.143.43
      生存概率Survival probability/%0.351.000.770.17
      检验数据Validation data(n=35)林分年龄Stand age134033.510.43184738.09.46
      胸径Diameter at breast height(DBH)/cm6.8729.1117.585.079.1232.4220.245.97
      平均高Mean height/m5.1117.2311.802.047.0518.9213.273.09
      林分密度/(株·hm-2)Stand density/(tree·ha-1)900.03 178.01 802.7907.8500.02 077.01 226.7578.4
      林分断面积/(m2·hm-2)Stand basal area/(m2·ha-1)9.0270.5837.0912.5410.9785.6645.0715.06
      立地质量Site index/m5.2416.3313.113.136.9819.2415.373.55
      生存概率Survival probability/%0.411.000.740.19
      注:n为样本数。Note:n is sample number.
    • 不同水平断面积生长均是与年龄相关的动态变化过程,研究中选择立地质量、林分密度作为研究不同水平断面积生长的动态约束指标,并用来估计不同水平断面积动态变化过程。

    • 地位指数是评价森林立地质量的重要指标,动态地位指数模型是随时间变化反应任意年龄立木生长的真实状态[19],同时该模型还具有多态性与多渐近线两个特点。王冬至等[20]基于广义代数差分法(GADA)建立了华北落叶松人工林动态地位指数模型,模型表达式为:

      $$ {H_2} = \exp \left( {{X_0}} \right){\left( {1 - \exp \left( { - 0.086{t_2}} \right)} \right)^{1.981 + 1/{X_0}}} $$ (1)
      $$ {L_0} = \ln \left( {1 - \exp \left( { - 0.086{t_1}} \right)} \right) $$ (2)
      $$ {X_0} = 0.5\left( {\left( {\ln {H_1} - 1.981{L_0}} \right) + \sqrt {{{\left( {\ln {H_1} - 1.981{L_0}} \right)}^2} - 4{L_0}} } \right) $$ (3)

      式中:H1H2分别为林分年龄在t1t2时的林分优势高,m;X0L0为与立地相关的参数。

    • 树木死亡和生存是一个复杂多变的过程,因此预测立木生存力比较困难。在林分密度模型中,假设株数为N,优势高为H,当前林分株数和优势高关系可表述为:dN/dH=a1Ha2Na3a1, a2, a3为模型估计参数值。此假设在林分生长模型中被广泛使用[3, 21-23],对t1t2两个时间内的方程两边进行微分,可得到每公顷株数的转换方程:

      $$ {N_2} = {\left[ {N_1^{1 - {a_3}} + {a_1}\frac{{{a_3} - 1}}{{{a_2} + 1}}\left( {H_2^{{a_2} + 1} - H_1^{{a_2} + 1}} \right)} \right]^{1/\left( {1 - {a_3}} \right)}} $$ (4)
      $$ {N_2} = 1\;000{\left[ {{{\left( {\frac{{{N_1}}}{{1\;000}}} \right)}^{1 - {a_3}}} + {a_1}\frac{{{a_3} - 1}}{{{a_2} + 1}}\left( {{{\left( {\frac{{{H_2}}}{{10}}} \right)}^{{a_2} + 1}} - {{\left( {\frac{{{H_1}}}{{10}}} \right)}^{{a_2} + 1}}} \right)} \right]^{1/\left( {1 - {a_3}} \right)}} $$ (5)

      式中:N1N2为调查期初和期末每公顷株数,株/hm2H为林分优势高,m;a1a2a3为模型参数。

    • 林分生产力W=GHG为林分断面积,m2/hm2H为林分优势高,m。林分生产力的变化可以反映林分增长量和死亡率之间的差异,本研究基于林分生产力模型来构建林分断面积生长模型。在纯林中,增长量可以由b1Hb2Nb3表示,死亡率可表示为-k(W/N)(dN/dH)=-kW(dlnN/dH),b1b2b3为模型参数,k为死亡木与活立木大小的比值,且为常数。因此假定林分生产力模型[3, 22-24]为:dW/dH = b1Hb2Nb3+kW(dlnN/dH)。在这种情况下,断面积转移方程可描述为:

      $$ {G_2} = N_2^{{b_3}}\left( {{G_1}{H_1}N_1^{ - {b_3}} + {b_1}\left( {H_2^{{b_2} + 1} - H_1^{{b_2} + 1}} \right)/\left( {{b_2} + 1} \right)} \right)/{H_2} $$ (6)
      $$ {G_2} = {\left( {\frac{{{N_2}}}{{1\;000}}} \right)^{{b_3}}}\left( {{G_1}\frac{{{H_1}}}{{10}}{{\left( {\frac{{{N_1}}}{{1\;000}}} \right)}^{ - {b_3}}} + {b_1}\left( {{{\left( {\frac{{{H_2}}}{{10}}} \right)}^{{b_2} + 1}} - {{\left( {\frac{{{H_1}}}{{10}}} \right)}^{{b_2} + 1}}} \right)/\left( {{b_2} + 1} \right)} \right)/ \left( {\frac{{{H_2}}}{{10}}} \right) $$ (7)
      $$ {G_2} = {G_1}N_2^{1 - {b_1}H_2^{{b_2}}}N_1^{{b_1}H_2^{{b_2}} - 1}{\left( {\frac{{{H_2}}}{{{H_1}}}} \right)^{{b_3}}} $$ (8)

      式中:G1G2为林分年龄在t1t2时的林分断面积,m2/hm2

    • 单木生长模型通常被用来预测单木水平胸径和断面积生长过程,本研究中选用Cao[24]和Andrea等[18]构建的单木生长模型作为华北落叶松人工林单木生长模型的基础模型,其表达式为:

      $$ {g_2} = {g_1} + {c_0}g_1^{{c_1}}G_1^{{c_2}}\exp \left( {{c_3}\frac{{{d_1}}}{{d_1^2}}} \right) $$ (9)
      $$ {g_2} = {g_1} + {c_0}{\left( {\frac{{{t_2}}}{{{t_1}}}} \right)^{{c_1}}}H_1^{{c_2}}G_1^{{c_3}}d_1^{{c_4}} $$ (10)

      式中:g1g2为年龄在t1t2时的单木断面积,m2d1为年龄在t1时的胸径,cm;c0, c1c2c3c4为模型参数估计值。

    • 研究中采用分解法,通过调整林分株数与断面积参数,将单木水平和林分水平断面积生长模型进行融合。为了准确预测立木死亡率,基于两种不同形式的逻辑斯蒂方程来预测单木生存概率和林分密度[25]。假设确定一个生存概率与单位面积株数,二者相乘就可以得到林分水平生存概率。在林分内,当P=0时林木就被认为是死亡,当P=1时林木就被认为是生存,每个样地林木生存概率由以下公式进行计算。

      $$ P = \frac{1}{{1 + \exp \left( {{f_0} + {f_1}{d_1} + {f_2}{G_1} + {f_3}\frac{{{d_1}}}{{d_1^2}}} \right)}} $$ (11)
      $$ P = \frac{1}{{1 + \exp \left( {{f_0} + {f_1}{N_1} + {f_2}{G_1} + {f_3}{d_1}} \right)}} $$ (12)

      式中:P为生存概率;f0, f1, f2, f3为模型参数估计值。

    • 为了分析不同水平对模型预测精度产生的影响,采用交叉验证法[26-27]从组合方程中得出林分变量,运用几何平均回归算法[20]对模型参数进行估计。这种方法有效地结合了不同模型信息,可以降低模型产生的误差,进而提高模型预测精度[28]。林分密度和林分断面积公式可表述为:

      $$ N_2^{\rm{C}} = {W^{\rm{N}}}N_2^{\rm{T}} + \left( {1 - {W^{\rm{N}}}} \right)N_2^{\rm{S}} $$ (13)
      $$ G_2^{\rm{C}} = {W^{\rm{G}}}G_2^{\rm{T}} + \left( {1 - {W^{\rm{G}}}} \right)G_2^{\rm{S}} $$ (14)

      式中:N2CG2C为林分密度和林分断面积预测值,N2TG2T为单木水平模型估计值,N2SG2S为林分水平模型估计值,WNWG为权重。

    • 基于常见的统计量对不同模型拟合精度和预测精度进行比较,即用平均绝对误差(MAD)、均方根误差(RMSE)、相对均方根误差(RRMSE)、确定系数(R2)及绝对误差(Bias)等指标对模型精度进行检验。不同统计量计算公式如式(15)~(19)所示。

      $$ \text{MAD=}\sum\limits_{i=1}^{n}{\left| {{y}_{i}}-{{\widehat{y}}_{i}} \right|}/n $$ (15)
      $$ \text{RMSE=}\sqrt{{{\sum\limits_{i=1}^{n}{\left( {{y}_{i}}-{{\widehat{y}}_{i}} \right)}}^{2}}/n-p} $$ (16)
      $$ \text{Bias=}\sum\limits_{i=1}^{n}{\left( {{y}_{i}}-{{\widehat{y}}_{i}} \right)}/n $$ (17)
      $$ \text{RRMSE=}\sqrt{\sum\limits_{i=1}^{n}{{{\left( {{y}_{i}}-{{\widehat{y}}_{i}} \right)}^{2}}}/\left( n-1 \right)}/\sum\limits_{i=1}^{n}{{{{\hat{y}}}_{i}}/n} $$ (18)
      $$ {{R}^{2}}\text{=1-}\sum\limits_{i=1}^{n}{{{\left( {{y}_{i}}-{{\widehat{y}}_{i}} \right)}^{2}}}/\sum\limits_{i=1}^{n}{{{\left( {{y}_{i}}-{{{\bar{y}}}_{i}} \right)}^{2}}} $$ (19)

      式中:yi、${{\widehat{y}}_{i}}$、y分别为观测值、预测值和平均值,n为样本数,p为模型参数个数。

    • 华北落叶松林分水平和单木水平断面积生长模型及林分密度预测方程的参数估计值、标准误(SE)、参数95%置信区间、确定系数(R2)、绝对误差(Bias)及均方根误差(RMSE)如表 2所示。在林分密度预测模型即模型(4)、(5)中,模型(4)的绝对误差(Bias)及均方根误差(RMSE)小于模型(5),且模型(4)的确定系数(R2)较大,并解释了林分密度90%以上的变异;在林分断面积预测模型即模型(6)、(7)和(8)中,模型(7)的绝对误差(Bias)及均方根误差(RMSE)均小于模型(6)和模型(8),且模型(7)的确定系数(R2)最大,并解释了林分断面积91%以上的变异;因此,确定模型(4)和模型(7)为华北落叶松人工林最优的林分密度模型和林分断面积预测模型。在单木断面积生长预测模型即模型(9)和(10)中,模型确定系数(R2)、绝对误差(Bias)及均方根误差(RMSE)均较为接近,但模型(9)略优于模型(10),二者均能解释单木水平断面积生长90%以上的变异,因此在预测单木水平断面积生长时可选择其中任一模型。运用拟合最优模型对林分水平与单木水平断面积及株数进行了预测,模型残差分布如图 1所示。

      表 2  不同水平模型参数估计与统计检验

      Table 2.  Parameter estimating and goodness-of-fit statistics in different levels

      模型Model参数Parameter估计值EstimateSER2BiasRMSE
      M4a10.005 10.003 50.900 21.952 52.186 3
      a20.406 50.399 5
      a31.271 30.526 7
      M5a1-0.405 80.095 40.804 32.005 62.334 2
      a2-0.738 70.415 8
      a3-0.045 20.009 0
      M6b1-0.000 80.000 30.812 31.280 93.155 0
      b27.143 82.056 7
      b33.222 40.684 5
      M7b124.978 61.529 10.914 40.899 71.833 4
      b21.215 30.191 7
      b30.799 20.040 5
      M8b10.999 30.000 80.867 11.008 42.304 1
      b20.000 50.000 3
      b3-2.853 60.677 6
      M9c00.011 20.002 30.913 20.097 50.918 7
      c1-4.284 30.163 3
      c24.396 30.151 7
      c3-24.096 12.385 7
      M10c00.000 20.000 10.908 80.105 81.007 6
      c10.756 70.111 3
      c2-0.528 10.262 4
      c31.093 00.081 3
      c4-0.042 20.012 1
      注:M4~M10对应文中公式(4)~(10)。Note: M4-M10 correspond to formula (4)-(10) of the paper, respectively.

      图  1  不同水平最优模型残差分布

      Figure 1.  Residual distribution of different level optimal models

    • 应用两种不同形式的逻辑斯蒂方程来预测单木生存概率,模型参数估计值、标准误(SE)、参数95%置信区间、确定系数(R2)、绝对误差(Bias)及均方根误差(RMSE)如表 2所示,其中模型(11)的绝对误差(Bias)及均方根误差(RMSE)小于模型(12),且模型(11)的确定系数(R2)较大,并解释了生存概率91%以上的变异;采用受试者工作特性曲线(ROC)对模型预测效果进行判别,其中Hosmer-Lemeshow统计量的显著性水平为P<0.001,获得的ROC曲线下面积为0.906,在左侧曲线上升较快,至顶后接近水平向右延伸,可以较好地预测林木生存概率,表明模型的判别效果较好。

      表 3  不同水平生存概率模型参数估计与统计检验

      Table 3.  Parameter estimating and goodness-of-fit statistics in different levels of survival probability

      模型Model参数Parameter估计值EstimateSER2BiasRMSE
      M11f0-3.294 00.780 40.910 41.022 42.077 1
      f10.000 50.000 3
      f20.009 20.002 9
      f324.111 617.138 5
      M12f0-3.663 00.616 00.904 61.075 82.332 4
      f111.015 14.372 4
      f20.036 00.007 8
      f30.007 90.002 3
      注:M11、M12对应文中公式(11)、(12)。Note: M11,M12 correspond to formula (11),(12) of the paper, respectively.

      图  2  单木生存概率ROC曲线

      Figure 2.  ROC curve for the individual tree mortality model fitted

    • 用林分生长模型得出的每公顷株数和断面积(N2SG2S)与单木生长模型得到的(N2TG2T)合并,得出方程(13)和(14)的参数估计值即WN=0.276 7、WG=0.168 5。研究中采用逻辑斯蒂模型对林木生存概率进行了预测,根据林分水平、单木水平断面积预测模型及组合预测方程,对每公顷株数和林分断面积预测误差进行了比较(如表 4所示),结果表明在预测林分株数和林分断面积时,组合预测方程预测效果最好,其预测值的平均绝对误差(MAD)、均方根误差(RMSE)和相对均方根误差(RRMSE)最小。

      表 4  林分、单株生长模型与组合方程预测林分变量评价

      Table 4.  Evaluation statistics for predicted stand variables from the stand and individual-tree growth models and forecast combination

      变量
      Variable
      每公顷株数预测Predicted tree number per ha林分断面积预测Predicted stand basal area
      RMSERRMSEMADRMSERRMSEMAD
      林分生长Stand growth305.739 47.151.933 82.054 78.440.880 1
      单木生长Individual tree growth364.008 713.152.406 12.931 615.660.973 4
      组合预测Combination forecast284.446 35.421.603 21.664 2-5.610.791 6
    • 为了提高对不同水平断面积模型的预测精度,本研究通过约束参数法、分解法分别建立了不同水平的断面积预测模型,并采用组合预测法将不同水平最优预测模型进行组合,以提高模型预测精度及不同水平断面积预测模型的兼容性。

      通过约束参数法建立了不同水平断面积预测模型,其中单木断面积生长模型的参数具有生物学意义,林分初始密度越大,单木断面积增长越慢,当林分密度达到最小值时,单木断面积增长迅速。竞争是影响单木断面积生长的重要因素,而立地指数在断面积生长模型中对立木胸径生长影响较小,可能是模型中除了立地质量参数外,还包含了反应立木大小的参数(胸径与胸径平方的比),Crecente-Campo等[29]也得到了同样的结果。

      由于单木分解变量能够预测林分水平变量,单木水平模型预测精度取决于林分水平模型的预测精度[5, 16-17]。为了分析不同水平对模型精度产生的影响,研究中采用逻辑斯蒂方程来预测立木生存概率,结合生长率变量[30]拟合单木断面积模型和生存概率模型参数,许多学者[5, 29, 31-32]曾做过类似的研究。在单木生长模型中,使用立木大小变量和竞争变量来构建模型,Wykoff[33]和Monserud等[34]做过同样的研究,并将断面积生长模型在混交林中做了发展与应用,解释了方差是由于立木大小和竞争变量引起的。

      采用组合预测法对不同水平断面积进行预测,并认为不同水平断面积模型的预测结果是无偏的,且预测模型的协方差是已知的。然而,如果协方差是未知的,就可能被忽略,每个模型的预测可以通过其均方误差逆加权进行计算。从组合方法中得出林分变量[17],这种方法有效地结合了不同模型信息,降低模型产生的误差,进而提高了模型预测精度。组合预测法能够准确预测不同水平林分密度、生存概率、单木断面积、林分断面积[28],因此,在结构复杂的混交林或天然次生林中如何应用还需进一步研究。

参考文献 (34)

目录

    /

    返回文章
    返回