高级检索

留言板

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

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

基于贝叶斯模型平均法的森林火灾预测模型构建研究

白海峰 刘晓东 牛树奎 何亚东

白海峰, 刘晓东, 牛树奎, 何亚东. 基于贝叶斯模型平均法的森林火灾预测模型构建研究[J]. 北京林业大学学报. doi: 10.12171/j.1000-1522.20200173
引用本文: 白海峰, 刘晓东, 牛树奎, 何亚东. 基于贝叶斯模型平均法的森林火灾预测模型构建研究[J]. 北京林业大学学报. doi: 10.12171/j.1000-1522.20200173
Bai Haifeng, Liu Xiaodong, Niu Shukui, He Yadong. Study on Construction of Forest Fire Prediction Model Based on Bayesian Model Averaging Method——Take Dali Prefecture, Yunnan Province as an example[J]. Journal of Beijing Forestry University. doi: 10.12171/j.1000-1522.20200173
Citation: Bai Haifeng, Liu Xiaodong, Niu Shukui, He Yadong. Study on Construction of Forest Fire Prediction Model Based on Bayesian Model Averaging Method——Take Dali Prefecture, Yunnan Province as an example[J]. Journal of Beijing Forestry University. doi: 10.12171/j.1000-1522.20200173

基于贝叶斯模型平均法的森林火灾预测模型构建研究

doi: 10.12171/j.1000-1522.20200173
基金项目: 国家自然科学基金项目(31770696)
详细信息
    作者简介:

    白海峰。主要研究方向:林火生态。Email:haifengbai@sina.com 地址:100083 北京市海淀区清华东路35号北京林业大学生态与自然保护学院

    通讯作者:

    刘晓东,博士,教授。主要研究方向:林火生态。Email:xd_liu@bjfu.edu.cn 地址:同上

  • 中图分类号: S762.2

Study on Construction of Forest Fire Prediction Model Based on Bayesian Model Averaging Method——Take Dali Prefecture, Yunnan Province as an example

  • 摘要:   目的  本文基于贝叶斯模型平均方法,结合二项逻辑斯蒂回归模型,构建云南省大理州森林火灾发生预测模型,以期提高林火预测精度,为研究地区林火管理提供技术支持。  方法  利用2000—2013年大理州林火数据及对应的气象数据,分别运用二项逻辑斯蒂回归模型和贝叶斯模型平均法,对该地区森林火灾对气象因子的响应进行实证分析。二项逻辑斯蒂回归模型为单一模型,建模前通过对各解释变量进行多重共线性检验,剔除有显著共线性的解释变量,然后通过逐步回归法,筛选最终变量并进行参数拟合。贝叶斯平均模型为组合模型,基于贝叶斯模型平均法建模时,采用奥卡姆窗的方法来适当调整模型空间,并以5个最优模型的后验概率作为权重进行加权建模。将全样本数据随机分成80%的训练样本和20%的测试样本,基于训练样本建立模型,对测试样本进行预测,通过对比观测值和预测值计算模型的准确率。  结果  通过二项逻辑斯蒂模型拟合,优度为0.783,预测精度为0.718。通过贝叶斯平均模型拟合,优度为0.868,预测精度为0.807。2个模型预测结果对比显示,在训练集中,贝叶斯平均模型的预测准确率比二项逻辑斯蒂回归模型高9.3%;在测试集中,贝叶斯平均模型的预测准确率比二项逻辑斯蒂回归模型高8.9%。  结论  在基于气象因子的大理州林火发生预测模型构建研究中,贝叶斯平均模型的拟合优度和预测精度均高于二项逻辑斯蒂模型,表明贝叶斯模型平均法具有一定的现实应用意义,可用于提高研究地区林火预测精度,有利于森林火灾的决策管理。
  • 图  1  研究区2000—2013年火点分布

    Figure  1.  The distribution of fire points in the study area from 2000 to 2013

    图  2  Step_LR模型的ROC曲线图

    AUC为曲线下面积;下同。AUC, area under curve;the same below.

    Figure  2.  ROC curves of Step_LR model

    图  3  BMA模型可视化

    模型编号Model No. 注:图示根据奥卡姆窗被选中的98个模型以及每个模型各自选中的变量。横轴为模型编号,宽度表示该模型的后验概率大小,纵轴为解释变量代码。红色表示该变量与被解释变量存在正相关关系,蓝色表示存在负相关关系,无颜色即表示该变量没有被选入该模型。Note: The figure shows the 98 models selected according to the Occam’s window and the variables selected by each model. The x-axis refers to the model number, the width represents the posterior probability of the model, and the y-axis is equidistant, showing the code of each explanatory variable. The red indicates that the variable has a positive correlation with the explained variable, the blue indicates that there is a negative correlation, and the variable without colour is not selected into the model.

    Figure  3.  BMA model visualization

    图  4  BMA_LR模型的ROC曲线图

    Figure  4.  ROC curves of BMA_LR model

    表  1  模型变量的基本统计描述

    Table  1.   The basic statistical description of model variables

    模型变量 Model variables变量代码 Variable code最小值 Min.最大值 Max.均值 Mean标准差 SD
    日平均风速 Daily average wind speed/(m·s−1) WIN_avg 0.80 10.80 3.64 1.38
    日最大风速 Daily maximum wind speed/(m·s−1) WIN_max 3.10 20.60 9.21 2.35
    日照时数 Sunshine hours/h SSD 2.20 12.20 9.24 1.80
    日平均本站气压 Daily average pressure/kPa PRS_avg 79.33 80.75 80.05 0.22
    日平均气温 Daily average temperature/℃ Tavg 4.20 24.10 15.88 3.34
    日最高气温 Daily maximum temperature/℃ Tmax 12.10 31.00 23.34 3.17
    日最低气温 Daily minimum temperature/℃ Tmin −0.80 18.20 8.59 3.92
    日平均水汽压 Daily average water vapor pressure/kPa VP_avg 0.27 1.68 0.71 0.21
    日平均相对湿度 Daily average relative humidity/% RH_avg 21.00 72.00 41.46 8.27
    日最小相对湿度 Daily minimum relative humidity/% RH_min 6.00 46.00 18.78 5.71
    20—20 时降水量 20:00—20:00 precipitation/mm Pre 0 3.00 0.03 0.23
    细小可燃物湿度码 Fine fuel moisture code FFMC 79.33 97.56 94.6 1.70
    粗腐殖质湿度码 Duff moisture code DMC 18.18 342.68 113.29 52.55
    干旱码 Drought code DC 61.91 660.71 373.32 98.28
    初始蔓延指数 Initial spread index ISI 1.44 15.63 10.01 1.99
    累积指数 Build up index BUI 23.92 339.76 128.54 48.63
    火险天气指数 Fire weather index FWI 5.35 49.05 33.47 7.14
    火点 Fire point Fire 0 1 0.50 0.50
    注:各模型变量样本数为1 102。Note: The sample number of each model variable is 1 102.
    下载: 导出CSV

    表  2  变量的多重共线性检验

    Table  2.   Multicollinearity test of variables

    变量 VariablesWIN_avgWIN_maxSSDPRS_avgTmaxTminRH_avgRH_minFFMCISI
    VIF值 VIF value8.671.409.137.686.678.128.113.791.959.35
    下载: 导出CSV

    表  3  Step_LR模型参数拟合

    Table  3.   The parameters estimation of Step_LR model

    变量
    Variables
    估计系数
    Estimate
    标准误差
    Std.error
    Z
    Z value
    P
    P value
    (Intercept) −13.006 2.225 −2.923 0.003
    WIN_max 0.013 0.003 3.932 0.000
    Tmax 0.049 0.004 11.281 0.000
    Tmin −0.030 0.004 −8.448 0.000
    RH_avg −0.129 0.013 −9.903 0.000
    RH_min 0.055 0.020 2.733 0.006
    FFMC 0.077 0.023 2.091 0.037
    下载: 导出CSV

    表  4  基于贝叶斯后验概率的模型平均

    Table  4.   Model average based on Bayesian posterior probability

    Variablesp! = 0SDmodel 1model 2model 3model 4model 5
    Intercept 100 45.070 −6.028 82.170 −5.490 86.100 −3.648
    WIN_avg 1.3 0.002
    WIN_max 11.5 0.003
    SSD 5.9 0.004
    PRS_avg 23.4 0.005 −0.011 −0.012
    Tavg 26.8 0.021 0.048
    Tmax 98.0 0.015 0.054 0.054 0.028 0.063 0.046
    Tmin 26.6 0.011 −0.024
    VP_avg 95.6 0.023 −0.081 −0.082 −0.077 −0.102 −0.064
    RH_avg 19.7 0.036 −0.038
    RH_min 33.9 0.033 0.061
    Pre 0.0 0
    FFMC 0.3 0.004
    DMC 16.6 0.005
    DC 6.0 0.001
    ISI 0.0 0
    BUI 70.0 0.006 0.008 0.007 0.007 0.007 0.007
    FWI 8.3 0.013
    nVar 3 4 5 5 4
    BIC −7 039 −7 038 −7 038 −7 038 −7 037
    post prob 0.062 0.057 0.056 0.039 0.035
    下载: 导出CSV

    表  5  Step_LR模型和BMA_LR模型中最终指标体系及预测准确率

    Table  5.   The final indicator system and prediction accuracy in the Step_LR & BMA_LR model

    模型 Model模型指标体系 The indicator system预测准确率 Prediction accuracy/%
    训练集 The training sample (80%)测试集 The test sample (20%)
    Step_LR modelWIN_max, Tmax, Tmin, RH_min, RH_avg, FFMC73.371.8
    BMA_LR modelPRS_avg, Tavg, Tmax, Tmin, VP_avg, RH_avg, RH_min, BUI82.680.7
    下载: 导出CSV
  • [1] Rigo D D, Giorgio L, Durrant T H, et al. Forest fire danger extremes in Europe under climate change: variability and uncertainty[M]. Luxembourg: Publications Office of the European Union, 2017.
    [2] 田晓瑞, 宗学政, 舒立福, 等. ENSO事件对中国森林火险天气的影响[J]. 应用生态学报, 2020, 31(5):65−73.

    Tian X R, Zong X Z, Shu L F, et al. Impacts of ENSO events on forest fire weather of China[J]. Chinese Journal of Applied Ecology, 2020, 31(5): 65−73.
    [3] 白夜, 武英达, 贾宜松, 等. 2019—2020年澳大利亚气候异常与山火爆发的关系分析及应对策略[J]. 中国应急救援, 2020(2):23−27. doi:  10.3969/j.issn.1673-5579.2020.02.006

    Bai Y, Wu Y D, Jia Y S, et al. Link between climate anomaly and Australia bushfires in 2019-2020[J]. China Emergency Rescue, 2020(2): 23−27. doi:  10.3969/j.issn.1673-5579.2020.02.006
    [4] 赵凤君, 舒立福. 森林草原火灾扑救安全学[M]. 北京: 中国林业出版社, 2015.

    Zhao F J, Shu L F. Forest and grassland fire fighting safety[M]. Beijing: China Forestry Publishing House, 2015.
    [5] 岳超, 罗彩访, 舒立福, 等. 全球变化背景下野火研究进展[J]. 生态学报, 2020, 40(2):385−401.

    Yue C, Luo C F, Shu L F, et al. A review on wildfire studies in the context of global change[J]. Acta Ecologica Sinica, 2020, 40(2): 385−401.
    [6] Marlon J R, Bartlein P J, Gavin D G, et al. Long-term perspective on wildfires in the western USA[J]. Proceedings of the National Academy of Sciences of the United States of America, 2012, 109(9): 3203−3204.
    [7] Westerling A L. Increasing western US forest wildfire activity: sensitivity to changes in the timing of spring[J]. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 2016, 371: 1−10.
    [8] 潘登, 郁培义, 吴强. 基于气象因子的随机森林算法在湘中丘陵区林火预测中的应用[J]. 西北林学院学报, 2018, 33(3):175−183.

    Pan D, Yu P Y, Wu Q. Application of random forest algorithm on the forest fire prediction based on meteorological factors in the hilly area, central Hunan Province[J]. Journal of Northwest Forestry University, 2018, 33(3): 175−183.
    [9] North M P., Stephens S L, Collins B M, et al. Reform forest fire management[J]. Science, 2015, 349: 1280−1281. doi:  10.1126/science.aab2356
    [10] Fischer A P, Spies T A, Steelman T A, et al. Wildfire risk as a socioecological pathology[J]. Frontiers in Ecology and the Environment, 2016, 14(5): 276−284. doi:  10.1002/fee.1283
    [11] Zhang G, Wang M, Liu K. Forest fire susceptibility modeling using a convolutional neural network for Yunnan Province of China[J]. International Journal of Disaster Risk Science, 2019, 10(3): 386−403. doi:  10.1007/s13753-019-00233-1
    [12] Murphy T E, Tsang S W, Leo L S, et al. Bayesian model averaging for selection of a risk prediction model for death within thirty days of discharge: The Silver-Ami study.[J]. International Journal of Statistics in Medical Research, 2019, 8: 1−7. doi:  10.6000/1929-6029.2019.08.01
    [13] Huang H, Liang Z, Li B, et al. Combination of multiple data-driven models for long-term monthly runoff predictions based on Bayesian model averaging[J]. Water Resources Management, 2019, 33(9): 3321−3338. doi:  10.1007/s11269-019-02305-9
    [14] 王倩, 师鹏飞, 宋培兵, 等. 基于贝叶斯模型平均法的洪水集合概率预报[J]. 水电能源科学, 2016(6):64−66.

    Wang Q, Shi P F, Song P B, et al. Multi-model ensemble flood probability forecasting based on BMA[J]. Water Resources and Power, 2016(6): 64−66.
    [15] 张畅, 陈新军. 海洋环境因子对澳洲鲐亲体补充量关系的影响: 基于贝叶斯模型平均法的研究[J]. 海洋学报, 2019, 41(2):99−106.

    Zhang C, Chen X J. The effect of environmental factors on stock-recruitment relationship of spotted mackerel-based on Bayesian model averaging method[J]. Haiyang Xuebao, 2019, 41(2): 99−106.
    [16] 李丽琴. 云南省森林火灾发生与气象因子之间的关系研究[D]. 北京: 北京林业大学, 2010.

    Li L Q. Study on the relationship between forest fires and the meteorological factors in Yunnan[D]. Beijing: Beijing Forestry University, 2010.
    [17] 周明昆, 王永平, 高月忠. 气象因子对云南大理森林火灾的影响[J]. 四川林业科技, 2012, 33(6):96−99. doi:  10.3969/j.issn.1003-5508.2012.06.022

    Zhou M K, Wang Y P, Gao Y Z. Effects of meteorological factors on forest fires in Dali, Yunnan[J]. Journal of Sichuan Forestry Science and Technology, 2012, 33(6): 96−99. doi:  10.3969/j.issn.1003-5508.2012.06.022
    [18] Martell D L, Otukol S, Stocks B J. A logistic model for predicting daily people-caused forest fire occurrence in Ontario[J]. Canadian Journal of Forest Research, 1987, 17(5): 394−401. doi:  10.1139/x87-068
    [19] 苏漳文, 刘爱琴, 郭福涛, 等. 福建林火发生的驱动因子及空间格局分析[J]. 自然灾害学报, 2016, 25(2):110−119.

    Su Z W, Liu A Q, Guo F T, et al. Driving factors and spatial distribution pattern of forest fire in Fujian Province[J]. Journal of Natural Disasters, 2016, 25(2): 110−119.
    [20] 于建龙, 刘乃安. 我国大兴安岭地区森林雷击火发生的火险天气等级研究[J]. 火灾科学, 2010, 19(3):131−137. doi:  10.3969/j.issn.1004-5309.2010.03.004

    Yu J L, Liu N A. Lightning - caused wildland fire weather danger rating in Daxing’anling region[J]. Fire Safety Science, 2010, 19(3): 131−137. doi:  10.3969/j.issn.1004-5309.2010.03.004
    [21] Bisquert M, Caselles E, Sánchez J M, et al. Application of artificial neural networks and logistic regression to the prediction of forest fire danger in Galicia using MODIS data[J]. International Journal of Wildland Fire, 2012, 21(8): 1025−1029. doi:  10.1071/WF11105
    [22] Oliveira S, Oehler F, San-Miguel-Ayanz J, et al. Modeling spatial patterns of fire occurrence in Mediterranean Europe using multiple regression and random forest[J]. Forest Ecology and Management, 2012, 275(4): 117−129.
    [23] 陈岱. 基于Logistic回归模型的大兴安岭林火预测研究[J]. 林业资源理, 2019(1):116−122.

    Chen D. Prediction of forest fire occurrence in Daxing’an Mountains based on logistic regression model[J]. Forest Resources Management, 2019(1): 116−122.
    [24] Raftery A E, Gneiting T, Balabdaoui F, et al. Using Bayesian model averaging to calibrate forecast ensembles[J]. Monthly Weather Review, 2005, 133(5): 1155−1174. doi:  10.1175/MWR2906.1
    [25] 梁慧玲, 林玉蕊, 杨光, 等. 基于气象因子的随机森林算法在塔河地区林火预测中的应用[J]. 林业科学, 2016, 52(1):89−98.

    Liang H L, Lin Y R, Yang G, et al. Application of random forest algorithm on the forest fire prediction in Tahe Area based on meteorological factors[J]. Scientia Silvae Sinicae, 2016, 52(1): 89−98.
    [26] 顾先丽, 吴志伟, 张宇婧, 等. 气候变化背景下江西省林火空间预测[J]. 生态学报, 2020, 40(2):667−677.

    Gu X L, Wu Z W, Zhang Y J, et al. Prediction research of the forest fire in Jiangxi Province in the background of climate change[J]. Ecological Sinica, 2020, 40(2): 667−677.
    [27] Chang Y, Zhu Z L, Bu R C, et al. Predicting fire occurrence patterns with logistic regression in Heilongjiang Province, China[J]. Landscape Ecology, 2013, 28(10): 1989−2004. doi:  10.1007/s10980-013-9935-4
    [28] Guo F T, Su Z W, Wang G Y, et al. Understanding fire drivers and relative impacts in different Chinese forest ecosystems[J]. Science of the Total Environment, 2017, 605: 411−425.
    [29] Flannigan M D, Krawchuk M A, Groot W J D, et al. Implications of changing climate for global wildland fire[J]. International Journal of Wildland Fire, 2009, 18(5): 483−507. doi:  10.1071/WF08187
    [30] Loepfe L, Rodrigo A, Lloret F. Two thresholds determine climatic control of forest fire size in Europe and northern Africa[J]. Regional Environmental Change, 2014, 14(4): 1395−1404. doi:  10.1007/s10113-013-0583-7
    [31] 蔡奇均, 曾爱聪, 苏漳文, 等. 基于Logistic回归模型的浙江省林火发生驱动因子分析[J]. 西北农林科技大学学报, 2020, 48(2):108−115.

    Cai Q J, Zeng A C, Su Z W, et al. Driving factors of forest fire in Zhejiang Province based on logistic regression model[J]. Journal of Northwest A & F University, 2020, 48(2): 108−115.
  • [1] 袁硕, 李超, 陈昊, MuhammadAmir, 刘晓东.  福建省将乐县生物防火林带阻隔网空间布局与规划 . 北京林业大学学报, doi: 10.12171/j.1000-1522.20190365
    [2] 曾爱聪, 郭新彬, 郑文霞, 魏帽, 靳全锋, 郭福涛.  基于MODIS卫星火点数据的浙江省林火时空动态变化特征 . 北京林业大学学报, doi: 10.12171/j.1000-1522.20190483
    [3] 周浪, 樊坤, 瞿华, 吕媛媛, 张正宜.  基于Sparse-DenseNet模型的森林火灾识别研究 . 北京林业大学学报, doi: 10.12171/j.1000-1522.20190371
    [4] 吴燕, 胡琦鹏, 白永超, 陈露, 唐仲秋, 侯智霞.  笃斯越桔伴生植物矿质元素与土壤肥力因子的多元分析 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20180420
    [5] 姚丹丹, 徐奇刚, 闫晓旺, 李玉堂.  基于贝叶斯方法的蒙古栎林单木枯死模型 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20180260
    [6] 邵艳莹, 吴秀芹, 张宇清, 秦树高, 吴斌.  内蒙古地区植被覆盖变化及其对水热条件的响应 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20170414
    [7] 王坚, 刘晓东.  福建三明2000—2011年森林火灾PM2.5排放量的估算 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20180075
    [8] 朱娅坤, 秦树高, 张宇清, 张举涛, 邵艳莹, 高岩.  毛乌素沙地植被物候动态及其对气象因子变化的响应 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20180020
    [9] 宫大鹏, 康峰峰, 刘晓东.  新巴尔虎草原火时空分布特征及对气象因子响应 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20170402
    [10] 段文军, 王成, 张昶, 宋阳, 郝泽周, 徐心慧, 金一博, 王子研.  夏季3种生境森林内空气颗粒物变化特征 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20160358
    [11] 杨路明, 秦树高, 刘振, 朱林峰, 刘峰.  毛乌素沙地不同下垫面地表凝结水形成过程 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20150375
    [12] 张璇, 张会兰, 王玉杰, 王云琦, 刘春霞, 杨坪坪, 潘声雷.  缙云山典型树种树干液流日际变化特征及与气象因子关系 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20150389
    [13] 姚丹丹, 雷相东, 张则路.  基于贝叶斯法的长白落叶松林分优势高生长模型研究 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20140221
    [14] 李明泽, 王雪, 高元科, 付瑜, 范文义.  大兴安岭植被指数年际变化及影响因子分析 . 北京林业大学学报, doi: 10.13332/j.1000-1522.20140220
    [15] 曹林, 代劲松, 徐建新2, 许子乾, 佘光辉.  基于机载小光斑LiDAR 技术的亚热带森林参数信息优化提取 . 北京林业大学学报, doi: 10.13332/j.cnki.jbfu.2014.05.009
    [16]
    陈锋, 林向东, 牛树奎, 王叁, 李德, 
    气候变化对云南省森林火灾的影响 . 北京林业大学学报,
    [17] 徐爱俊, 方陆明, 楼雄伟.  基于可见光视频的森林火灾识别算法 . 北京林业大学学报,
    [18] 张彦林, 冯仲科, 姚山, 董斌, .  城市森林火灾时空蔓延模型构建 . 北京林业大学学报,
    [19] 张军国, 李文彬, 韩宁, 阚江明.  基于ZigBee无线传感器网络的森林火灾监测系统的研究 . 北京林业大学学报,
    [20] 田勇臣, 刘少刚, 赵刚, 胡健, 李文彬.  森林火灾蔓延多模型预测系统研究 . 北京林业大学学报,
  • 加载中
图(4) / 表 (5)
计量
  • 文章访问数:  37
  • HTML全文浏览量:  13
  • PDF下载量:  4
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-06-19
  • 修回日期:  2021-01-07

基于贝叶斯模型平均法的森林火灾预测模型构建研究

doi: 10.12171/j.1000-1522.20200173
    基金项目:  国家自然科学基金项目(31770696)
    作者简介:

    白海峰。主要研究方向:林火生态。Email:haifengbai@sina.com 地址:100083 北京市海淀区清华东路35号北京林业大学生态与自然保护学院

    通讯作者: 刘晓东,博士,教授。主要研究方向:林火生态。Email:xd_liu@bjfu.edu.cn 地址:同上
  • 中图分类号: S762.2

摘要:   目的  本文基于贝叶斯模型平均方法,结合二项逻辑斯蒂回归模型,构建云南省大理州森林火灾发生预测模型,以期提高林火预测精度,为研究地区林火管理提供技术支持。  方法  利用2000—2013年大理州林火数据及对应的气象数据,分别运用二项逻辑斯蒂回归模型和贝叶斯模型平均法,对该地区森林火灾对气象因子的响应进行实证分析。二项逻辑斯蒂回归模型为单一模型,建模前通过对各解释变量进行多重共线性检验,剔除有显著共线性的解释变量,然后通过逐步回归法,筛选最终变量并进行参数拟合。贝叶斯平均模型为组合模型,基于贝叶斯模型平均法建模时,采用奥卡姆窗的方法来适当调整模型空间,并以5个最优模型的后验概率作为权重进行加权建模。将全样本数据随机分成80%的训练样本和20%的测试样本,基于训练样本建立模型,对测试样本进行预测,通过对比观测值和预测值计算模型的准确率。  结果  通过二项逻辑斯蒂模型拟合,优度为0.783,预测精度为0.718。通过贝叶斯平均模型拟合,优度为0.868,预测精度为0.807。2个模型预测结果对比显示,在训练集中,贝叶斯平均模型的预测准确率比二项逻辑斯蒂回归模型高9.3%;在测试集中,贝叶斯平均模型的预测准确率比二项逻辑斯蒂回归模型高8.9%。  结论  在基于气象因子的大理州林火发生预测模型构建研究中,贝叶斯平均模型的拟合优度和预测精度均高于二项逻辑斯蒂模型,表明贝叶斯模型平均法具有一定的现实应用意义,可用于提高研究地区林火预测精度,有利于森林火灾的决策管理。

English Abstract

白海峰, 刘晓东, 牛树奎, 何亚东. 基于贝叶斯模型平均法的森林火灾预测模型构建研究[J]. 北京林业大学学报. doi: 10.12171/j.1000-1522.20200173
引用本文: 白海峰, 刘晓东, 牛树奎, 何亚东. 基于贝叶斯模型平均法的森林火灾预测模型构建研究[J]. 北京林业大学学报. doi: 10.12171/j.1000-1522.20200173
Bai Haifeng, Liu Xiaodong, Niu Shukui, He Yadong. Study on Construction of Forest Fire Prediction Model Based on Bayesian Model Averaging Method——Take Dali Prefecture, Yunnan Province as an example[J]. Journal of Beijing Forestry University. doi: 10.12171/j.1000-1522.20200173
Citation: Bai Haifeng, Liu Xiaodong, Niu Shukui, He Yadong. Study on Construction of Forest Fire Prediction Model Based on Bayesian Model Averaging Method——Take Dali Prefecture, Yunnan Province as an example[J]. Journal of Beijing Forestry University. doi: 10.12171/j.1000-1522.20200173
  • 近年来,全球范围内森林火灾的发生频次不断增大,不仅威胁着人们的生命财产安全,同时还造成巨大的资源损失和环境破坏。2017年为欧洲有史以来遭受森林火灾侵袭最为严重的年份之一,火灾对西班牙、葡萄牙和意大利等国均造成了灾难性的事件[1],2019—2020年火险期澳大利亚的野火持续了5个月,引起了全球性持续关注[2-3]。我国也是世界上森林火灾发生较为严重的国家之一,2000—2015年发生森林火灾的次数平均为7 632次/年,火场总面积平均为230 622 hm2/年,受害森林面积平均为94 864 hm2/年,人员伤亡平均为111人/年[4]。在全球气候变暖的背景下,火活动将增加,火险期延长,野火发生概率升高[5-8],准确估计火灾概率对于减少林火的负面影响起着至关重要的作用[9-11]。因此,准确的森林火灾发生预测模型的构建是林火发生预测的重要手段,对防控森林火灾具有十分重要的意义。

    随着对森林火灾发生的认识不断深入,从早期的线性模型到计数模型,森林火灾发生与气象因子的关系模型结构日趋复杂,然而这种模型结构的复杂化不能降低林火发生与气象因子关系模型的不确定性。而贝叶斯模型平均法则是近年来文献报告中处理模型不确定性一个很好的方法,虽目前尚未将其应用于森林火灾预测,但在其他方面如医疗、水文、渔业的应用[12-15]表明,贝叶斯模型平均算法具有较高的预测效果和模型稳定性。云南省地处我国西南,是我国三大林区之一,但同时也是我国的林火高发区。由于受地理位置、地形地势、气候以及森林资源分布和人为活动等影响,云南的森林火灾较为严重[16]。因此,本文基于云南省大理州林火数据和气象数据,运用R统计软件,分别应用逻辑斯蒂回归模型和贝叶斯模型平均法建立林火发生的单一模型和组合模型,通过不同模型拟合结果的对比分析,判断贝叶斯模型平均法在构建区域林火预测模型中的适用性。

    • 大理州地处云南省中部偏西,地理位置98°52′ ~ 101°03′ E、24°41′ ~ 26°42′ N,总面积2. 945 9万km2,其中山区面积占总面积的93.4%。地势西北高,东南低,平均海拔2 090 m,地貌复杂多样,土壤类型以紫色土和红壤土为主。气候属于低纬高原季风气候,立体气候特点显著,干湿季节分明,年温差小、日温差大。年均气温15.7 ℃,年均降水量836 mm,年日照时数2 072 ~ 2 693 h,无霜期225 ~ 345 d。大理州森林面积173.37万hm2,森林覆盖率61.22%,森林蓄积量1.1376亿m3,主要优势树种有云南松(Pinus yunnanensis)、华山松(P. armandii)、铁杉(Tsuga chinensis)、冷杉(Abies fabri)、马尾杉(Phlegmariurus phlegmaria)、思茅松(Pinus kesiya)等。

      大理州属于森林火灾高发区,该地区冬春季气温高、降水少、风力大,极易发生森林火灾[17]。据统计,2000—2013年间,大理州共计发生森林火灾552次(图1),过火面积高达11 027 hm2,其中有林地面积为7 002 hm2。由于大理州森林火灾集中发生于1—6月份(551次),因此本文以该时段为研究对象,选取了该时段内的区域林火数据和气象数据。

      图  1  研究区2000—2013年火点分布

      Figure 1.  The distribution of fire points in the study area from 2000 to 2013

    • 本文森林火灾数据来自云南省大理州森林防火指挥部办公室2000—2013年林火发生情况数据,包括火灾发生时间、起火地点、起火原因、过火面积等。

      气象数据来自国家气象科学数据中心基本气象数据,为大理气象站(区站号:56751)2000—2013年1—6月的逐日气象数据。本文通过对气象数据进行预处理,剔除掉缺测数据过多的气象因子,保留11个气象因子:日平均风速(m/s)、日最大风速(m/s)、日照时数(h)、日平均本站气压(kPa)、日平均气温(℃)、日最高气温(℃)、日最低气温(℃)、日平均水汽压(kPa)、日平均相对湿度(%)、日最小相对湿度(%)、20:00—次日20:00降水量(mm)。

      另外,考虑到森林火灾的发生不仅与当时的气象条件有关,而且还与前期气象条件密切相关,气象因子对林火发生的影响具有一定滞后性,因此引入加拿大森林火险天气指标系统(fire weather index,FWI)。FWI指标体系以时滞——平衡含水率理论为基础,将气象条件和可燃物含水率有机地联系起来,通过天气条件的变化计算可燃物含水率的变化,进而确定潜在火险等级,这在一定程度上代表了当前气象条件与前期气象条件的综合效应,因此本文将FWI指标体系看成气象因子,与11个气象因子共同设为自变量,进行森林火灾发生概率模型构建。

    • 本文以火灾发生日气象数据为自变量,因变量Y = 1,同时在构建判别模型时,需要创建一定比例的非火点(Y = 0)作为对照。本文按1∶1比例选取对照点[18-19],对照非火点创建过程中遵循时间和空间上的完全随机[20]

      本文将全样本数据随机分成80%的训练样本和20%的测试样本,基于训练样本建立模型,对测试样本进行预测,通过对比观测值和预测值计算出模型的准确率。本文运用R软件中的glm函数进行二项逻辑斯蒂回归的计算,逐步回归是通过R软件中的step.glm函数实现的,贝叶斯模型平均算法是通过R软件中的BMA,MASS程序包实现的。

    • 逻辑斯蒂回归模型(Logistic regression model, LR)是目前国内外普遍运用的林火预测模型[18, 21-23]。将是否发生火灾的概率P设定为因变量YY = 1为发生林火,Y = 0为未发生林火),各气象因子设定为自变量X,建立Logistic回归模型,模型表达式如下:

      $$P = \frac{1}{{1 + {e^{ - \left( {{\beta _0} + {\beta _1}{X_1} + {\beta _2}{X_2} + \cdots + {\beta _n}{X_n}} \right)}}}}$$ (1)

      式中:P为林火发生的概率;βi(i = 1, 2, …, n)为模型中影响林火发生的自变量Xii = 1, 2, …, n)的逻辑斯蒂回归系数,n为模型自变量的个数。

    • 贝叶斯模型平均法(Bayesian model averaging, BMA)是由Rafery等[24]提出的一种利用多模型组合进行概率预报的统计方法。BMA采用贝叶斯公式,将带有先验信息的未知参数分布与描述对象新信息的似然函数结合,获得对象的后验分布,基于后验分布对对象进行统计和推断,从而不断完善对对象的认识[14]

      BMA是通过估算潜在变量X的所有可能组合的模型并在需要的组合上构建加权平均来解决问题的。通常情况下,根据已有的主观信息指定先验分布,以模型后验概率(posterior model probabilities,PMP)为权重对可能的单项模型进行加权平均,以解释变量的后验包含概率(posterior inclusion probabilities,PIP)大小作为选择解释变量的客观标准[15]

      在Logistic回归模型中,基于BMA方法可得到的每个解释变量系数βi的后验概率,即为解释变量的后验包含概率PIP,并对参数向量$\theta \left( {{\beta _0}, \;{\beta _1}, \;\cdots,\;{\beta _n}} \right)$进行相应的统计推断,可以构建出BMA-Logistic回归模型,具体情况如下:

      单个模型的后验概率表达式为:

      $${\rm{P}}\left( {{M_r}{\rm{|}}D} \right) = \frac{{P\left( {D{\rm{|}}{M_r}} \right)P\left( {{M_r}} \right)}}{{\displaystyle\mathop \sum \nolimits_{j = 1}^k P\left( {D{\rm{|}}{M_j}} \right)P\left( {{M_j}} \right)}}$$ (2)
      $$P\left( {D{\rm{|}}{M_r}} \right) = \mathop \int \nolimits P\left( {D{\rm{|}}{\theta _r},{M_r}} \right)P\left( {{\theta _r}{\rm{|}}{M_r}} \right)d{\theta _r}$$ (3)

      式中:${M_r}$代表模型空间D中的第r个模型,k表示模型的个数,$P\left( {{M_r}} \right)$为模型${M_r}$的先验概率分布,$P\left( {D{\rm{|}}{M_r}} \right)$为模型${M_r}$对应的似然函数积分,${\theta _r}$表示模型${M_r}$的参数向量,$P\left( {{\theta _r}{\rm{|}}{M_r}} \right)$表示模型${M_r}$所对应的参数先验概率分布,$P\left( {D{\rm{|}}{\theta _r},{M_r}} \right)$表示模型${M_r}$所对应的似然函数。

      根据贝叶斯公式,参数向量$\theta $的后验密度分布是模型空间条件下参数$\theta $后验密度分布的加权平均值,权重为模型后验概率${\rm{P}}\left( {{M_r}{\rm{|}}D} \right)$。BMA还需要指定模型先验概率和参数的先验信息,即设定$P\left( {{M_r}} \right)$。一般情况,可以设定相等的模型先验概率,即设定均匀分布的模型空间,而对于参数的先验信息,可以设定单位信息先验(Unit information prior)。

    • 本文采用ROC(receiver operating characteristic)曲线和AUC(area under curve)对预测效果进行检验。AUC作为数值可以直观的评价分类器的好坏,值越大分类效果越好。一般认为,AUC值等于0.5时相当于一个完全的随机预测,在(0.5, 0.7]区间,准确性较低,在(0.7, 0.8]区间,准确性中等,在(0.8, 0.9]区间,准确性较好,在(0.9, 1]区间,具有高准确性[25-26]

      通过ROC曲线分析法可以得到模型的敏感度值和特异性值,通过约登指数(Youden index)公式(约登指数 = 敏感度值 + 特异性值 − 1)可以计算分类阈值(cut-off point),进而对预测概率进行分类。如果计算得到的预测概率值大于该阈值则认为会发生林火,小于该阈值则认为不会发生林火[25, 27]

    • 根据大理州2000—2013年1—6月火点数据,按1∶1比例选取对照点,应用SPSS软件对火点和对照点组成的全样本数据进行基本统计(表1)。

      表 1  模型变量的基本统计描述

      Table 1.  The basic statistical description of model variables

      模型变量 Model variables变量代码 Variable code最小值 Min.最大值 Max.均值 Mean标准差 SD
      日平均风速 Daily average wind speed/(m·s−1) WIN_avg 0.80 10.80 3.64 1.38
      日最大风速 Daily maximum wind speed/(m·s−1) WIN_max 3.10 20.60 9.21 2.35
      日照时数 Sunshine hours/h SSD 2.20 12.20 9.24 1.80
      日平均本站气压 Daily average pressure/kPa PRS_avg 79.33 80.75 80.05 0.22
      日平均气温 Daily average temperature/℃ Tavg 4.20 24.10 15.88 3.34
      日最高气温 Daily maximum temperature/℃ Tmax 12.10 31.00 23.34 3.17
      日最低气温 Daily minimum temperature/℃ Tmin −0.80 18.20 8.59 3.92
      日平均水汽压 Daily average water vapor pressure/kPa VP_avg 0.27 1.68 0.71 0.21
      日平均相对湿度 Daily average relative humidity/% RH_avg 21.00 72.00 41.46 8.27
      日最小相对湿度 Daily minimum relative humidity/% RH_min 6.00 46.00 18.78 5.71
      20—20 时降水量 20:00—20:00 precipitation/mm Pre 0 3.00 0.03 0.23
      细小可燃物湿度码 Fine fuel moisture code FFMC 79.33 97.56 94.6 1.70
      粗腐殖质湿度码 Duff moisture code DMC 18.18 342.68 113.29 52.55
      干旱码 Drought code DC 61.91 660.71 373.32 98.28
      初始蔓延指数 Initial spread index ISI 1.44 15.63 10.01 1.99
      累积指数 Build up index BUI 23.92 339.76 128.54 48.63
      火险天气指数 Fire weather index FWI 5.35 49.05 33.47 7.14
      火点 Fire point Fire 0 1 0.50 0.50
      注:各模型变量样本数为1 102。Note: The sample number of each model variable is 1 102.
    • 多重共线性(Multicollinearity)是指线性回归模型中解释变量之间存在某种密切相关关系。多重共线性是普遍存在的,通常情况下,适度的共线性不成问题,但严重的共线性会导致解释变量的显著性检验失去意义及模型估计产生一定偏差甚至无效。因此在涉及多个解释变量时,应首先对其进行多重共线性检验。本文运用方差膨胀因子(Variance inflation factor, VIF)诊断法对解释变量进行多重共线性检验。通常,VIF值越大,说明多重共线性就越显著,一般认为VIF大于10时,解释变量之间具有显著的共线性[28]。经检验,WIN_avg、WIN_max、SSD、PRS_avg、Tmax、Tmin、RH_avg、RH_min、FFMC、ISI共计10个气象因子的VIF小于10(表2),进入模型拟合阶段。

      表 2  变量的多重共线性检验

      Table 2.  Multicollinearity test of variables

      变量 VariablesWIN_avgWIN_maxSSDPRS_avgTmaxTminRH_avgRH_minFFMCISI
      VIF值 VIF value8.671.409.137.686.678.128.113.791.959.35
    • 基于训练样本,运用R软件中的glm函数对大理州火点和对照点以及WIN_avg、WIN_max、SSD、PRS_avg、Tmax、Tmin、RH_avg、RH_min、FFMC、ISI共计10个气象因子数据进行Logistic模型拟合,通过逐步回归法,逐步剔除模型中不显著的变量,得到的最终变量为WIN_max、Tmax、Tmin、RH_avg、RH_min、FFMC共计6个变量,选择最终变量进行模型拟合,来构建Step_Logistic(Step_LR)模型,拟合参数如表3所示。经逐步回归筛选出的最终变量与该地区林火发生均有显著相关性,其中WIN_max、Tmax、Tmin、RH_avg、RH_min在P < 0.01水平上显著相关,FFMC在P < 0.05水平上显著相关。

      表 3  Step_LR模型参数拟合

      Table 3.  The parameters estimation of Step_LR model

      变量
      Variables
      估计系数
      Estimate
      标准误差
      Std.error
      Z
      Z value
      P
      P value
      (Intercept) −13.006 2.225 −2.923 0.003
      WIN_max 0.013 0.003 3.932 0.000
      Tmax 0.049 0.004 11.281 0.000
      Tmin −0.030 0.004 −8.448 0.000
      RH_avg −0.129 0.013 −9.903 0.000
      RH_min 0.055 0.020 2.733 0.006
      FFMC 0.077 0.023 2.091 0.037
    • 应用ROC曲线分析法对Step_LR模型的预测能力进行拟合优度检验,并计算林火发生的分类阈值。图2为预测模型的ROC曲线图,ROC曲线下的面积(AUC值)为0.783,显著性水平P < 0.001,说明模型的拟合优度为中等水平。根据约登指数公式可得林火发生的分类阈值为0.498,以该值为分界点,林火发生的预测概率值(P)大于0.498视为有林火发生,小于0.498则视为无林火发生,进一步计算得到20%测试数据集对林火发生的预测概率为0.718。

      图  2  Step_LR模型的ROC曲线图

      Figure 2.  ROC curves of Step_LR model

    • 本文共有17个解释变量,利用各解释变量建立林火发生的Logistic模型,可能存在的模型个数高达217个,贝叶斯模型平均法可遍历模型空间中的每一个模型,根据各模型的后验概率来衡量其对林火发生的相对重要性,通过模型空间调整以确定较优模型。将Fire(1和0)设定为因变量Y,各气象因子设定为自变量X,建立BMA模型,使用R程序的“BMA package”进行计算。

      本文设定相等的模型先验概率即均匀分布的模型空间,对于参数的先验信息,设定为单位信息先验;同时,采用了奥卡姆窗(Occam’s window)的方法来适当调整模型空间,即减少一定的模型数量,本文设定当一个模型的后验概率(PMP)小于最佳模型后验概率的5%时,则从模型空间中被剔除。程序结果如图3所示,奥卡姆窗筛选了98个较优模型,每个模型包含了部分变量。

      图  3  BMA模型可视化

      Figure 3.  BMA model visualization

    • 本文基于贝叶斯估计的后验概率对98个较优模型做了贝叶斯模型平均,程序结果见表4

      表 4  基于贝叶斯后验概率的模型平均

      Table 4.  Model average based on Bayesian posterior probability

      Variablesp! = 0SDmodel 1model 2model 3model 4model 5
      Intercept 100 45.070 −6.028 82.170 −5.490 86.100 −3.648
      WIN_avg 1.3 0.002
      WIN_max 11.5 0.003
      SSD 5.9 0.004
      PRS_avg 23.4 0.005 −0.011 −0.012
      Tavg 26.8 0.021 0.048
      Tmax 98.0 0.015 0.054 0.054 0.028 0.063 0.046
      Tmin 26.6 0.011 −0.024
      VP_avg 95.6 0.023 −0.081 −0.082 −0.077 −0.102 −0.064
      RH_avg 19.7 0.036 −0.038
      RH_min 33.9 0.033 0.061
      Pre 0.0 0
      FFMC 0.3 0.004
      DMC 16.6 0.005
      DC 6.0 0.001
      ISI 0.0 0
      BUI 70.0 0.006 0.008 0.007 0.007 0.007 0.007
      FWI 8.3 0.013
      nVar 3 4 5 5 4
      BIC −7 039 −7 038 −7 038 −7 038 −7 037
      post prob 0.062 0.057 0.056 0.039 0.035

      表4可知,Tmax的后验包含概率为98.0%,VP_avg的后验包含概率为95.6%,BUI的后验包含概率为70.0%,这3个变量的后验包含概率较大,意味着这3个变量对林火发生有较大影响。其他解释变量的后验包含概率相对较小,说明这些解释变量对林火发生的影响力相对减弱。而从模型的稳定性角度看,最佳模型(model 1)的后验概率仅为0.062,前5个模型的累计后验概率为0.249,由此可见,模型的不确定性在该数据集中是相当大的。

      本文以5个最优模型的后验概率作为权重进行加权,来构建BMA_Logistic(BMA_LR)组合模型。加权整合后的平均模型为:

      $$P = \frac{1}{{1 + {{\rm{e}}^{ - \left( {7.233 - 0.001{\rm{PRS}}\_{\rm{avg}} + 0.003{\rm{Tavg}} + 0.012{\rm{Tmax}} - 0.001{\rm{Tmin}} - 0.020{\rm{VP}}\_{\rm{avg}} - 0.001{\rm{RH}}\_{\rm{avg}} + 0.002{\rm{RH}}\_{\rm{min}} + 0.002{\rm{BUI}}} \right)}}}}$$
    • 应用ROC曲线分析法对BMA_LR模型的预测能力进行拟合优度检验,并计算林火发生的分类阈值。图4为预测模型的ROC曲线图,且ROC曲线下的面积(AUC值)为0.868,显著性水平P < 0.001,说明模型的拟合优度较高。根据约登指数公式可得林火发生的分类阈值为0.562,以该值为分界点,林火发生的预测概率值大于0.562视为有林火发生,小于0.562则视为无林火发生,进一步计算得到20%测试数据集对林火发生的预测概率为0.807,结果显示模型具有较高的预测能力,可用于大理州林火发生的预测预报。

      图  4  BMA_LR模型的ROC曲线图

      Figure 4.  ROC curves of BMA_LR model

    • Step_LR模型和BMA_LR模型的最终指标体系对比如表5所示。在BMA_LR模型中,Tmax、VP_avg和BUI 3个变量的后验包含概率最大,这3个变量只有Tmax进入了Step_LR模型的指标体系。

      表 5  Step_LR模型和BMA_LR模型中最终指标体系及预测准确率

      Table 5.  The final indicator system and prediction accuracy in the Step_LR & BMA_LR model

      模型 Model模型指标体系 The indicator system预测准确率 Prediction accuracy/%
      训练集 The training sample (80%)测试集 The test sample (20%)
      Step_LR modelWIN_max, Tmax, Tmin, RH_min, RH_avg, FFMC73.371.8
      BMA_LR modelPRS_avg, Tavg, Tmax, Tmin, VP_avg, RH_avg, RH_min, BUI82.680.7
    • 根据模型变量选择结果,分别对Step_LR模型和BMA_LR模型预测准确率进行计算。在训练集中,BMA_LR模型的预测准确率比Step_LR模型高9.3%,在测试集中,BMA_LR模型的预测准确率比Step_LR模型高8.9%,结果显示BMA_LR模型具有较高的预测能力,可适用于云南省大理州林火预测。

    • 本文分别应用二项逻辑斯蒂回归模型和贝叶斯模型平均法建立了云南省大理州林火发生的单一模型和组合模型,比较模型拟合结果表明:贝叶斯平均模型的预测准确率比逐步逻辑斯蒂回归模型高8.9%,说明基于贝叶斯模型平均法的组合模型林火预测效果优于逐步逻辑斯蒂回归模型,可用于该地区的林火发生预测预报。

      本文构建的贝叶斯平均模型包含日平均本站气压、日平均气温、日最高气温、日最低气温、日平均水汽压、日平均相对湿度、日最小相对湿度、BUI共计8个气象因子,表明这8个气象因素是影响该地区林火发生的重要驱动气象因子,特别是日平均水汽压、日最高气温和BUI后验包含概率均在70%以上,表明这3个气象因子是影响该地区林火发生的主要驱动气象因子。

      在构建逻辑斯蒂回归模型时,首先进行了多重共线性检验,剔除了日平均气温、日平均水汽压、20:00—次日20:00时降水量、DC、DMC、BUI、FWI共计7个气象因子。在这7个被剔除的气象因子中,日平均气温、日平均水汽压、BUI这3个气象因子则进入了最终贝叶斯平均模型,可见进行多重共线性检验存在将重要驱动因子提前去除的风险。这是因为进行多重共线性检验考虑的是解释变量之间的相关性,而并未考虑被剔除的变量对被解释变量的影响,因此,提前对解释变量进行多重共线性检验,并根据检验结果对存在多重共线性的变量进行剔除,则有可能剔除了对林火发生有显著影响的气象因子,这就会造成预测准确率的降低。

      本文通过分析森林火灾对气象因子的响应讨论了贝叶斯模型平均法对林火发生气象因子的选择及林火预测的适用性,实际上,自然环境中影响林火发生的生态因子和非生态因子很多。已有研究表明,森林火灾发生还与植被类型、林型、地形、人为活动、社会经济[23, 28-31]等多个因素息息相关,因此若想更加及时有效的预测森林火灾的发生,还需进一步分析上述因素中的重要驱动因子及其对林火发生的影响。今后研究中应纳入以上变量,筛选关键驱动因子,建立组合模型,进一步提高预测效果。

参考文献 (31)

目录

    /

    返回文章
    返回