Identification and prediction models of driving factors for forest fires in Daxing’an Mountains of Inner Mongolia, northern China
-
摘要:目的
选择和验证适合研究区的林火预测模型，明确火灾发生的关键驱动因素并绘制火险区划图，为森林火灾预防和管理工作提供科学依据和决策支持。
方法基于1981—2020年历史火灾数据，融合多源（气象条件、地形、植被、人类活动和社会经济）数据，比较4种机器学习方法在内蒙古大兴安岭林火预测中的适用性，同时，基于显著影响火灾发生的驱动因素绘制火灾发生可能性地图和火险区划图。
结果（1）增强回归树模型接受者操作特性曲线下的面积值为0.967，随机森林模型的AUC为0.947，均表现出优异的预测性能。Logistic回归模型和Gompit回归模型的预测准确率较上两种略低，AUC分别为0.852、0.851，也满足研究区的基本预测要求。（2）气象因素气温日较差、日最小相对湿度是影响内蒙古大兴安岭林火发生的主导因素；海拔在驱动因素的相对重要性排序中位居前列；人类活动和社会经济因素（如距公路的距离、距火灾瞭望塔的距离、人均GDP等）对林火发生也有一定影响。（3）内蒙古大兴安岭东部和东南部存在大面积火灾中、高风险区，北部中俄边境和西南部中蒙边境也有较高的火灾风险。火灾发生前一年秋季防火期的平均气温、平均地表温度等因素会影响第2年森林火灾的发生。
结论与其他3种模型相比，增强回归树模型是最适合内蒙古大兴安岭林火发生的预测模型。气象因子、海拔显著影响内蒙古大兴安岭林火发生，人类活动和社会经济因素对火灾发生也有一定的影响。研究区的中高火险区域主要集中在东部和东南部，北部和西南部也有一定的火灾风险。
Abstract:ObjectiveThis paper aims to select and validate suitable forest fire prediction models for the study area, identify key driving factors of fire occurrence, and map fire risk zoning, then providing scientific basis and decision support for forest fire prevention and management.
MethodUsing historical fire data from 1981 to 2020 and integrating multi-source data (meteorological conditions, topography, vegetation, human activities, and socio-economic factors), the applicability of four machine learning methods in predicting forest fires in the Daxing’an Mountains of Inner Mongolia of northern China was compared. Based on the significant factors influencing fire occurrence, maps of fire occurrence probability and fire risk zoning were generated.
Result(1) The boosted regression tree model (BRT) showed an area under the curve (AUC) value of 0.967, and the random forest model (RF) achieved an AUC of 0.947, both demonstrating excellent predictive performance. The predictive accuracy of the Logistic regression model (LR) and the Gompit regression model (GR) was slightly lower than former two models, but still met the basic predictive requirements for the study area, with AUC values of 0.852 and 0.851, respectively. (2) Meteorological factors, such as diurnal temperature range and daily minimum relative humidity, were the dominant factors influencing forest fires in the Daxing’an Mountains of Inner Mongolia. Elevation also ranked high in the relative importance of driving factors. Human activities and socio-economic factors, such as distance to roads, distance to fire lookout towers, and per capita GDP, also had some influence on fire occurrence. (3) Large areas of medium to high fire risk were present in the eastern and southeastern parts of the Daxing’an Mountains of Inner Mongolia, while the northern China-Russia border and the southwestern China-Mongolia border also exhibited elevated fire risk. Factors such as average temperature and average surface temperature during fire prevention period in autumn of previous year influenced forest fire occurrences in the following year.
ConclusionAmong the four models compared, the BRT was identified as the most suitable one for predicting forest fire occurrence in the Daxing’an Mountains of Inner Mongolia. Meteorological factors and elevation significantly influence fire occurrence, while human activities and socio-economic factors also have a certain impact on the occurrence of fires. The high and medium fire risk areas are primarily concentrated in the eastern and southeastern parts of the study area, with some fire risks presenting in the northern and southwestern regions.
-
氮（N）、磷（P）、钾（K）元素是植物生长发育的三要素，科学合理的施肥能促进植物开花、结实，显著提高产量和品质，而盲目不合理甚至过量施肥不仅限制作物产量的提高，增加肥料投入成本，降低肥料的吸收利用效率，造成不必要的浪费，还会减少土壤微生物和动物数量，破坏土壤结构，造成土壤养分富集和环境污染[1-2]。林木种子园由于经营周期长，长期施肥方案的确定尤为重要。结合物种自身需肥规律、生长土壤供肥特性、肥料效应及养分间的相互作用规律，确定氮、磷、钾最优配比和适宜施用量的配方施肥方式是农林综合运用现代科技成果的一项高效技术措施。目前配方施肥研究在农作物和园艺方面卓有成效[3-5]，林业配施方面的研究多针对苗期和林分生长[6-11]，种子园配施研究在白桦（Betula platyphylla）[12]、杉木（Cunninghamia lanceolata）[13-15]、油松（Pinus tabuliformis）[16]、樟子松（P. syluestris var. mongolica）[17]、云南松（P. yunnanensis）[18]、红松（P. koraiensis）[19]、华山松（P. armandi）[20]、马尾松（P. massoniana）[21]等树种中也有开展，但长期连续多元素配比施肥对落叶松种子园种子产量性状影响方面的研究鲜见报道。
日本落叶松（Larix kaempferi）是我国主要速生用材树种，良种主要来源于种子园。其结实大、小年间隔期长，导致其产量极不稳定，难以满足社会对林木良种的需求。自然未改良的种子仍是造林用种子的主体，而良种苗木的匮乏限制了人工林质量的提升[22]。近年来，针对日本落叶松无性系种子园施肥方面的研究发现，单独施加N素、P素或者N、P、K不同配比施肥都可影响种实性状、产量和品质，并提高抗逆性[23-27]。另外，不同肥料间具有耦合、协同和激励作用[26]。因此，为促进肥料利用率，有效提高种子园产量和品质，在前人研究的基础上，本研究以日本落叶松无性系初级种子园中无性系为试验材料，通过不同施肥处理探明提高种子园球果产量的最佳施肥配比和最适用量，旨在为提出节约环保高效的落叶松种子园经营技术提供依据。
1. 试验材料与方法
1.1 试验地自然概况
试验地设在甘肃小陇山沙坝国家落叶松/云杉良种基地的日本落叶松无性系初级种子园，位于小陇山林区腹地，处在山麓坡脚、沟谷地带，地势平坦，紧邻河流，灌溉方便。地理位置为105°51′27″ ~ 105°54′51″N、34°07′28″ ~ 34°10′22″E，海拔1 560 ~ 2 019 m，相对高差700 ~ 1 000 m，坡度一般为28° ~ 40°。地貌上属秦岭西段，与青藏高原相邻，受西北干冷气流影响大。试验地属我国暖温带南缘与北亚热带的过渡地带，气候温暖湿润，为温带湿润区。气候特点是冬无严寒夏无酷暑，春季气温回升快，常有倒春寒和春旱发生。年均气温7.2 ℃，极端最高气温32.0 ℃，极端最低气温–27.0 ℃。年均降水量460.0 ~ 1 012.2 mm，降雨集中分布在7、8、9月，一般春秋雨量占年降雨量的70% ~ 80%，降雨年际变化较大，枯、丰年降水量相差近1倍，年蒸发量989 ~ 1 658 mm，平均相对湿度68% ~ 78%；干燥度1.25 ~ 1.50，径流深度100 mm左右。平均年日照时数1 520 ~ 2 313 h，≥10 ℃的有效积温2 480 ℃，无霜期154 ~ 185 d。土壤以山地褐土、山地棕壤为主，土层厚度30 ~ 60 cm，土壤有机质含量1.71% ~ 5.35%，pH值5.72 ~ 7.90，显弱酸性或弱碱性，速效N、P缺乏，且土壤熟化程度高[28]。
施肥前（2008年4月中旬）测定了种子园土壤肥力本底，取样和分析方法参见文献[29]，测定结果见表1。
表 1 日本落叶松无性系种子园施肥前土壤N、P、K养分含量测定结果Table 1. Determination results of soil nitrogen (N), phosphorus (P) and potassium (K) nutrients before fertilization in Larix kaempferi clonal seed orchard土层深度
Depth of soil layer/cmpH值
pH value全氮
Total N/%全磷
Total P/%全钾
Total K/%速效氮
Available N/(mg·kg−1)速效磷
Available P/(mg·kg−1)速效钾
Available K/(mg·kg−1)0 ~ 30 6.170 ± 0.360 0.310 ± 0.068 0.049 ± 0.026 1.570 ± 0.180 37.490 ± 7.680 3.760 ± 0.650 121.450 ± 21.360 1.2 试验材料
试验材料为日本落叶松嫁接无性系初级种子园内正常开花结实的12号无性系分株。建园材料为2年生嫁接苗（嫁接接穗来自于辽宁日本落叶松优树），1989年依据无性系错位配列设计营建种子园，株行距6 m × 6 m，整地规格80 cm × 80 cm × 60 cm。为消除无性系效应，施肥用母株为同一无性系分株[30]，且生长势基本一致，树冠完整，树体没有经过人为修剪、环剥等经营措施，树体平均胸径29.8 cm，平均树高15.6 m，平均冠幅8.9 m。
1.3 试验设计与方法
试验以尿素（含有效N 46.4%）、过磷酸钙（含有效P2O5 ≥ 12.0%）、硫酸钾（含有效K2O ≥ 52.0%）3种肥料作为氮（N）、磷（P）、钾（K）3种营养元素来源。设计了N1P1K1 、N1P2K1和N1P1K2 3种配方肥，分别对应的尿素、过磷酸钙、硫酸钾的质量比为1∶1∶1、1∶2∶1和1∶1∶2，即N P K纯养分含量比为1∶0.26∶0.79、1∶0.52∶0.79和1∶0.26∶1.58。
施肥方式采用环状沟施肥，即：以试验母株树干基部为圆心，树冠正投影1/2处为半径（落叶松细根≤ 5 mm）生物量在树干和树冠投影边缘的中点处最大，且随着年龄的增加呈现出细根向深层土壤增生的趋势）[31-32]，用铁锹或䦆头开挖宽约30 cm、深30 ~ 40 cm的环状沟渠，按照试验设计将肥料均匀的洒在沟渠内，然后覆土。沿试验母株树冠正投影外1.5 m左右垂直埋入深60 cm的一圈塑料膜与其他母树隔离，并将同一试验区组布置在种子园相同小区内。
（1）球果产量对氮磷钾配比水平响应试验：2008—2010年春季（3月中旬前后）进行同一施肥量（2.5 kg/株）3种氮磷钾（N∶P∶K）营养元素配方施肥试验（表2），以不施肥作为对照（CK）。试验1设在种子园15大区，随机区组设计，每个区组施肥3株，共3个区组，4个处理共36株（含CK）。
表 2 氮磷钾配方施肥试验 kg/株Table 2. Experiment design of N, P and K formula fertilizationkg/plant 处理 Treatment 水平Ⅰ Level Ⅰ 水平Ⅱ Level Ⅱ 水平Ⅲ Level Ⅲ 全氮 Total N 全磷 Total P 全钾 Total K 全氮 Total N 全磷 Total P 全钾 Total K 全氮 Total N 全磷 Total P 全钾 Total K CK 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 N1P1K1 1.22 0.32 0.96 2.44 0.63 1.93 3.66 0.95 2.89 N1P2K 1.08 0.56 0.85 2.16 1.13 1.71 3.25 1.69 2.56 N1P1K2 0.88 0.23 1.39 1.76 0.46 2.78 2.64 0.69 4.17 注：水平Ⅰ. 2.5 kg/株；水平Ⅱ. 5.0 kg/株；水平Ⅲ. 7.5 kg/株。 Notes: level Ⅰ, 2.5 kg/plant; level Ⅱ, 5.0 kg/plant; level Ⅲ, 7.5 kg/plant. （2）球果产量对氮磷钾不同施肥水平连年配方施肥的响应试验：确定影响球果产量的最佳氮磷钾（N∶P∶K）营养元素配比后，2011—2018年3月（散粉前10 d左右）以水平Ⅰ（2.5 kg/株）、水平Ⅱ（5.0 kg/株）、水平Ⅲ（7.5 kg/株）开展连年施肥试验，各配方3水平中氮磷钾含量见表2，以不施肥作为CK。试验2设在种子园14大区，3个处理和CK按照随机区组设计，每个区组3株树，共3个区组，共36个分株（含CK）。
1.4 试验数据收集与分析
种子成熟后，单株采收球果，用电子天平称取球果鲜质量（精度0.1 ㎏）。
数据采用SPSS16.0进行方差分析（ANOVA、GLM）和多重比较，方差分析模型为：
xij=μ+αi+eij (1) xijk=μ+αi+βj(i)+eijk (2) 式中：
μ 为群体均值，xij 、xijk 为第i（i=1，2，3，⋯ ，n）个处理的第j（j=1，2，3，⋯ ，m）个k（k=1，2，3）重复的观测值，αi 表示A因子第i个处理的效应值，βj(i) 为B因子j处理i水平的效应值，eij 和eijk 为随机误差。2. 结果与分析
2.1 球果产量对氮磷钾配比水平的响应
按照CK处理的球果产量，将结实年份划分为3类：Ⅰ类为结实小年，单株平均球果产量 ≤ 12.50 kg/株；Ⅱ类为结实平年，12.50 kg/株 < 单株平均产量球果 ≤ 45.00 kg/株；Ⅲ类为结实大年，单株平均产量球果 ≥ 45.00 kg/株。
2008年为结实大年，4个处理下球果产量从大到小依次为：N1P2K1 > N1P1K2 > N1P1K1 > CK，分别达到（97.63 ± 0.87） kg/株、（95.12 ± 0.81） kg/株、（94.61 ± 0.83） kg/株以及（82.60 ± 0.89） kg/株，N1P2K1分别比CK、N1P1K1、N1P1K2提高了18.15%、3.17%和2.63%。多重比较显示，N1P2K1与CK、N1P1K1、 N1P1K2均存在显著差异，而N1P1K1 和N1P1K2差异不显著，但均与CK具有显著差异。2009年和2010年分别为不结实和或极少量结实年份，不纳入分析。
因此，对速效N、速效P缺乏的土壤进行配比施肥，N、P、K 3种配比皆可显著提高球果产量；但在N、P、K配比肥料中，在促进球果产量方面，提高P素含量比提高K素含量效果显著，且K素昂贵，增加K素含量会造成资源浪费，增加成本。
2.2 球果产量对氮磷钾N1P2K1配方肥的连年响应
2011年为Ⅲ类结实大年，方差分析和多重比较结果显示，在3个不同水平（2.5、5.0、7.5 kg/株）配施N1P2K1肥，球果平均产量与对照显著差异（F球果年产量 = 28.751 > F（0.05），P < 0.05），不同施肥水平均有效提高了球果年产量，其中5.0 kg/株处理效果最好，球果产量达到（75.56 ± 2.35） kg/株，与其他处理相比差异显著，比CK、2.5 kg/株、7.5 kg/株处理分别提高了20.17%、5.12%、4.34%。2.5 kg/株和7.5 kg/株处理间差异不显著，但均与CK具有显著差异，球果产量显著高于CK（表3）。
表 3 日本落叶松结实母树球果产量对N、P、K配比连年施肥响应的多重比较Table 3. Multiple comparison of response of cone yield of L. kaempferi seed-bearing mother plant toN, P, K ratio and successive fertilization施肥水平/（kg·株−1）
Fertilization
level/(kg·plant−1)结实母株球果平均年产量/（kg·株−1）
Average annual production of cones from fruit-bearing mother plant/(kg·plant−1)年度累积球果
产量/（kg·株−1）
Annual cumulative
cone yield/(kg·plant−1)2011 2012 2013 2014 2015 2016 2017 2018 0（CK） 62.88 ± 3.08c 0 0.00 ± 0.00c 55.33 ± 1.50c 0.00 ± 0.00c 0.43 ± 0.65c 56.22 ± 2.44b 0 174.87 ± 4.83d 2.5 71.88 ± 2.52b 0 3.71 ± 1.90a 72.26 ± 3.96a 24.59 ± 1.94a 34.22 ± 2.07a 60.78 ± 2.22a 0 267.43 ± 5.15a 5.0 75.56 ± 2.35a 0 3.26 ± 0.27b 70.51 ± 3.23a 22.30 ± 1.75b 29.62 ± 1.14b 58.24 ± 2.05b 0 259.49 ± 6.44b 7.5 72.42 ± 3.75b 0 3.12 ± 0.23b 66.24 ± 2.07b 20.93 ± 1.55b 33.17 ± 1.39a 57.12 ± 2.03b 0 253.01 ± 6.28c F值F value 28.751 590.088 59.997 538.117 1095.422 6.886 498.240 注：小写字母表示不同处理间在P < 0.05水平上的的差异显著性。数值为均值 ± 标准差。下同。Notes: different lowercase letters indicate the difference significance at P < 0.05 level between different treaments
(P < 0.05). The values are mean ± SD. Same as below.2012年为I类结实小年，施肥处理和CK均未结实，球果产量为0。2013年为I类结实小年，施肥处理出现少量结实，且与CK差异显著；2.5 kg/株处理球果产量达到（3.71 ± 1.90） kg/株，显著高于其他处理（P < 0.05）。2014年为III类结实大年，球果产量范围为55.33 ~ 72.26 kg/株；多重比较显示，2.5 kg/株处理与5.0 kg/株处理球果产量最高且二者差异不显著，但与CK和7.5 kg/株处理相比具有显著差异，其中效果最为显著的是2.5 kg/株处理，比CK提高了20%以上（表3）。2015年为I类结实小年，但施肥处理后球果产量达到II类平年水平；2.5 kg/株处理结实量最大，且与5.0 kg/株处理、7.5 kg/株处理具有显著差异（P < 0.05）。2016年为I类结实小年，但施肥处理后球果产量达到II类平年水平，本年度CK出现极少量结实；2.5 kg/株处理结实量最高，其次就是7.5 kg/株处理，二者球果产量差异不显著，但均与5.0 kg/株处理差异显著（P < 0.05），施肥水平对球果产量影响显著。2017年为III类结实大年，球果产量范围为56.22 ~ 60.78 kg/株；2.5 kg/株处理产量最高，达到（60.78 ± 2.22） kg/株，与CK和其他施肥水平具有显著差异（P < 0.05），而其他3个处理差异不显著。2018年出现极端灾害天气，在4月1—3日、4月6—7日连续2次降雪并发生冰冻，沙坝良种基地气象资料显示6—7日最低气温达到−12.7 ℃，结实母株受到严重冻害，不纳入分析。
分析以上情况可知，按照N1P2K1配比，当施肥量达到5.0 kg/株时，球果产量增加至最高水平；当施肥量达到7.5 kg/株时，球果产量没有随着施肥量的增加提高，反而略有降低趋势，与2.5 kg/株施肥水平下的球果产量相当。这表明5.0 kg/株的施肥量当年能够满足母株生长、结实对营养元素的需求，施肥不足或者过多施肥均会限制球果产量的增加。
连续施肥3年后，不同水平施肥下的母株出现了连年结实的现象，且2.5 kg/株施肥量下的母株结实量最高。随着施肥年限的延长，当施肥量 > 2.5 kg/株时，母株球果产量不显著增加，且随着施肥量的增加，球果产量有降低的趋势。N1P2K1不同水平连续施肥第5年和第6年后，母树球果产量虽然不及大年多，但是随着施肥年限的增加，球果产量呈增多的趋势。同为结实小年，施肥第5和第6年的球果产量明显高于第3年，且结实大年第7年的球果产量低于前2次结实大年，这可能与母树没有间歇、连续结实营养消耗有关。因此，N1P2K1配方肥能有效提高球果产量，且长期连年施肥效果显著。
2.3 不同年度间球果平均产量分析
2.3.1 结实大年球果产量分析
自然状态下日本落叶松的球果产量有明显的“大、小年”现象，间隔期一般为3 ~ 5年，这是因为结实需要消耗大量营养，树体需要通过自我调节来恢复树势。对N1P2K1配施8年中出现的大、小年分别进行施肥响应分析发现，连续施肥8年间，出现了3次结实大年，分别在施肥第1年（当年）、第4年和第7年。方差分析显示，结实大年中试验母树的球果产量对不同施肥水平的响应均具有显著差异（F球果年产量=45.460 > F（0.05），P < 0.05）。计算3次大年中的球果年平均产量，发现结实母株球果产量对不同施肥水平响应的大小依次为：2.5 kg/株 > 5.0 kg/株> 7.5 kg/株 > CK（表4）。4个处理可分为A、B两组，CK和7.5 kg/株处理为A组，2.5 kg/株处理和5.0 kg/株处理为B组。A、B两组内各处理差异不显著，A、B两组之间差异显著。随着施肥年限的延长，结实母株的球果年产量并没有随着施肥量的增加而增加，反而呈降低趋势，表明按照N1P2K1配比连年施肥，适宜的施肥量范围为2.5 ~ 5.0 kg/株。综合考虑成本，认为2.5 ~ 5.0 kg/株为结实大年的最适宜施肥量。
表 4 结实母株大小年球果产量对N1P2K1不同水平施肥响应的多重比较 kg/株Table 4. Multiple comparison of cone yield of seed-bearing mother plants in different levels of N1P2K1kg/plant 施肥水平
Fertilization level结实大年球果产量
Cone yield in fruit-bearing large year结实小年球果产量
Cone yield in fruit-bearing small year0（CK） 58.14 ± 4.15b 0.00 ± 0.00 b 2.5 68.31 ± 6.14a 12.50 ± 14.41a 5.0 68.10 ± 7.81a 11.04 ± 12.63a 7.5 65.26 ± 6.92b 11.44 ± 13.55a 方差分析表明，3次结实大年间母株球果产量具有显著差异（F球果年产量=109.199 > F（0.05），P < 0.05）但随着施肥年限延长，球果产量并没有显著增加（表5）。
表 5 结实母株大年间球果产量对N1P2K1不同水平连年施肥响应的多重比较Table 5. Multiple comparison in response of annual cone yield of seed-bearing mother plant to different levels of successive N1P2K1 fertilization in fruit-bearing large year结实年份
Year of fruitfulness球果产量/（kg·株−1）
Cone yield/(kg·plant−1)2017 58.09 ± 2.72b 2014 66.09 ± 7.21b 2011 70.68 ± 5.57a 2.3.2 结实小年球果产量分析
连续施肥8年间，出现了5次结实小年，分别是2012、2013、2015、2016和2018年，不同施肥水平下球果产量的大小依次为2.5 kg/株 > 7.5 kg/株 > 5.0 kg/株 > CK。方差分析显示，N1P2K1配方肥3个施肥水平下结实母株的球果产量均显著高于CK（F球果年产量 = 88.450 > F（0.05），P < 0.05），且不同施肥水平间差异不显著。不同施肥水平都能促进母树提高球果产量，其中2.5 kg/株处理的施肥量最低，表明在N1P2K1配比下，2.5 kg/株的施肥量能满足母树开花结实消耗的养分供给，起到调节树势的效果。
方差分析显示，5个结实小年年份之间，结实母株的球果产量对施肥年度的响应具有显著差异（F球果年产量=114.447 > F（0.05），P < 0.05）。除特殊年份外（譬如冰冻雪灾、病虫危害等），小年中结实母株的球果产量变化趋势与大年相反，在N1P2K1配比水平连年施肥下，随着施肥年限的增加，球果产量逐年增加（表6）。
表 6 结实母株小年间球果产量对N1P2K1不同水平连年施肥响应的多重比较Table 6. Multiple comparison in response of annual cone yield of seed-bearing mother plant to different levels of successive N1P2K1 fertilization in fruit-bearing small year结实年份
Year of fruitfulness球果产量/（kg·株−1）
Cone yield/(kg·plant−1)2012 0.00 ± 0.00c 2013 2.52 ± 1.51c 2015 17.38 ± 14.50b 2016 24.01 ± 10.18a 2018 0.00 ± 0.00c 2.4 结实母树球果累积产量对氮磷钾配比连年施肥的响应
连续施肥8年结实母株球果累积产量的方差分析表明，不同水平连年施肥间具有显著差异（F球果年累积产量=498.240 > F（0.05），P < 0.05）。多重比较显示，4个处理母树8年球果累积产量均具有显著差异（表3），2.5 kg/株、5.0 kg/株、7.5 kg/株处理对应的球果累积产量依次高出CK52.93%、48.39%和44.68%。其中，2.5 kg/株处理的球果累积产量最高，为（267.43 ± 5.15） kg，比5.0 kg/株处理和7.5 kg/株处理提高了3.06%和5.70%。不同施肥水平都显著增加了球果累积产量，但是当连年施肥量 > 2.5 kg/株的时候，球果累积产量随着施肥水平提高而呈下降趋势。
N1P2K1配比下不同水平连年施肥对应的结实母株的球果产量均显著高于CK。2.5 kg/株处理在7个年份中对应的年球果产量均为最高，且累计产量高于其他3个处理。
研究认为，综合考虑经济成本，日本落叶松无性系种子园适宜的施肥方式为以2.5 kg/株施加N1P2K1配方肥。
3. 讨 论
3.1 球果产量对氮磷钾不同配比水平的响应
N、P、K是植物生长发育必不可少的矿质营养元素，其中，N是合成叶绿素、蛋白质、核酸、部分激素、酶、氨基酸、多酚类等不可或缺的组成部分[33]，伴随植物生长发育生命过程的始终；P以多种途径参与植物体内的代谢过程，是植物体中核酸、蛋白质、磷脂等许多重要化合物的组成成分，占植物干质量的0.05% ~ 0.50%，对保障植物生长及产量起着不可替代的作用[26]；K是植物吸收的数量仅次于N素的矿物质元素，在植物体内不形成任何分子中的稳定结构部分，但具有维持细胞渗透压、调节气孔运动、保障酶活性、增加光合作用、促进体内糖分的运输、提高产量和品质等重要作用[34]。N素缺乏会导致叶绿素合成不足、光合作用减弱，植物的生长发育和结实就会受到限制；同样P素缺乏，会对N素的吸收、植物的光合作用、营养生长、生殖生长、产量和品质、生物量的分配、根系统结构、细根的形成和寿命等过程造成严重影响[26]。N、P、K 3元素之间具有耦合、协同和激励作用，单独施入N素或者P素虽然能改善种实性状、提高球果产量、增加结实母树抗逆性，但是单一元素的施用有很多弊端，会导致其他元素利用效率下降、甚至肥效丧失，从而造成资源浪费[1]和点源污染[35-37]。据统计，我国化肥利用效率数据中，氮肥为30% ~ 35%，磷肥为10% ~ 25%，钾肥为35% ~ 50%，肥料利用效率低[38]。
因此，合理配比施肥能有效提高肥料的吸收利用和转化效率[27]，葛艺早等[14]在杉木种子园通过氮磷钾合理配比精准施肥，使杉木种子千粒质量、种子内可溶性糖含量、蛋白质含量显著增加，有效提高了种子产量和品质。通过N、P、K合理配比施肥促进种子园母树开花结实、提高种子产量和品质在我国马尾松、樟子松、红松种子园，国外南方松（Pinus palustris）、火炬松（P. taeda）、湿地松（P. elliottii）、欧洲赤松（P. sylvestris）等树种中也有类似报道[17,21,38-39]。N素是土壤中最易耗竭和限制植物生长的元素之一[40]，自然土壤中有效P含量低，难以满足植物最佳生长的需要。施肥作为一种养分调控手段，能合理及时补充母树生长、开花、结实所需要的营养物质，提高种子园母树的种子产量和品质。
日本落叶松种子园速效N、P缺乏，单独施P素并不能显著改善日本落叶松结实母株种实性状，且施肥量不宜超过0.18 kg/株，否则P素将在土壤中累积[26]。本研究表明，N1P1K1、N1P2K1、N1P1K2 3种不同配比水平施肥都能有效提高结实母树的球果产量，但是N1P2K1配比下结实母树的球果产量最高。配施中增加P素含量能有效提高结实母树的球果产量，其既满足了结实母树对P素的最佳需求，还促进了母树对N素和其他营养元素的吸收、转化和利用[26]，而N素利用效率的提高能促进种子蛋白质、可溶性糖含量的提高，从而有效提高种子品质[14]。K素虽在植物生长开花结实过程中起着催化剂和酶的促进作用，但是配施中增加K素含量，球果产量提高不显著。K肥是极易淋溶的元素，含量受土壤母质的影响[29]，因此，N、P、K配比施肥中，提高K素含量不仅造成肥料的浪费还增加经济成本。所以日本落叶松种子园进行N、P、K配比施肥，不建议提高K素含量。综合以上分析，日本落叶松种子园提高结实母树球果产量适宜的配比水平为N1P2K1，单株施肥量为2.5 kg/株，N、P、K纯养分含量依次为1.08、0.85、0.56 kg/株。
N、P、K配比水平施肥能促进结实，提高肥料间的利用效率，这除了与配比水平有关外，还与不同树种、以及同一树种不同发育时（年龄）期、施肥时间、施肥次数、施肥方式等因素有关[40-42]，比如欧洲赤松7 ~ 10年生无性系种子园中，每年早春连施3年N（150 ~ 250 kg/hm2）、P（250 kg/hm2）、K（200 kg/hm2）配比效果最好[38]；油松中施肥促进开花结实试验表明，结合灌溉施肥显著高于单独施肥[43]；火炬松二代种子园夏季施N肥促进开花效果优于8月施肥[38]。同时，P素在土壤中移动性差、扩散速率低，且施入土壤中的至少70% ~ 90%的P素与Fe、Al和Ca结合固定成难溶性磷酸盐，很大程度上降低了施肥的有效性[44]。因此，本研究采用的施肥方式为：在每年春季散粉前6 ~ 10 d，在树冠正投影（细根主要分布区域）的1/2处，深度30 ~ 40 ㎝，采用环状沟施肥，缩短P在土壤中的移动距离，增加P素的当季利用效率，从而提高结实母株的球果年产量。
3.2 球果产量对氮磷钾（N1P2K1）连年不同水平施肥的响应
结实要消耗树体大量营养物质，日本落叶松在自然状态下大小年明显，在一次大量结实后，一般间隔3 ~ 5年，且大小年间种子产量相差极为悬殊[30]。本研究表明，N1P2K1配比下，不同水平环状沟连年施肥都能有效提高结实母株的球果年产量，连续施肥3年后不同施肥水平母株出现了连年结实的现象。不论是结实大年还是结实小年，2.5 kg/株施肥水平对母株年球果产量的提高效果最显著，结实大年中5.0 kg/株水平施肥次之。随N1P2K1配方施肥量提高，球果年产量在结实大年没有显著增加，反而呈减小的趋势，但在结实小年略有增加但差异不显著。连续8年施肥后，结实母株的球果累积产量不同水平施肥间均具有显著差异，且在2.5 kg/株施肥水平下结实母株的球果累积产量最高，分别比CK、5.0 kg/株、7.5 kg/株施肥水平下的产量增加了52.93%、48.39%和44.68%。因此，20年以上的落叶松结实母株，适宜的施肥水平应该为2.5 kg/株，其能够满足结实母株生长、开花结实对养分的需求，使母株结实消耗的大量营养物质得到及时补充，能有效调整结实周期。
3.3 球果产量对氮磷钾（N1P2K1）配方施肥年度间的响应
日本落叶松结实母株球果年产量年度间的差异是多重因素综合影响的结果。连续施肥第2年，不同水平施肥和CK对应的母树都没有结实，这主要是因为施肥当年为结实大年，虽然进行了施肥，但是结实消耗的大量营养物质还未得到有效补充；连续施肥第3年，除CK外，不同水平施肥的母树都开始有少量结实，表明连续施肥后，树体结实消耗的营养得到了补充，树势开始恢复，出现了少量结实。随着施肥年限的延长，除极端灾害性气候影响外，小年间的球果产量呈递增的趋势，达到平年结实水平，且年度间球果产量具有显著差异；而结实大年的球果产量与结实小年相反，随施肥年限的延长，球果产量递减。结实平年中对营养的消耗较低，连年施肥使营养得到持续补充，结实不再出现间隔期；而在结实大年中，结实短期内消耗母株大量营养，需要较长时间的恢复。
母树结实量除了受施肥效应影响外，还与本身遗传因素（结母树年龄）、气象因子（温度、光照、湿度、风速、降水等）、花粉、病虫害等因子有关，相关因素对种子产量和品质的影响在杉木、红松、华山松、马尾松等种子园中都有相似报道[45-46]。本试验中，2018年4月初受到了连续低温冰冻灾害性气候的影响，CK和不同水平施肥下的母树均没有结实。
3.4 问题与今后的研究方向
种子园内结实母株的球果产量除了受N、P、K配比水平、施肥量影响外，还与施肥时间、施肥次数、施肥方式等有关[41-42,47]，是多重因素综合影响的结果。本研究只进行了N、P、K配比3水平、N1P2K1最佳配比4个施肥水平下环状沟春季施肥的研究，且由于受试验株数的限制，N、P、K配比水平和N1P2K1配方施肥量水平数设计较少。因此，今后应在施肥时间、肥料配方和施肥方式等方面研究施肥对球果产量和种子品质影响[45,48]，同时开展落叶松种子园专用肥配方和适应低氮和低磷无性系良种选择研究[26]，有效提高日本落叶松种子园种子产量和品质，满足社会对良种的需求。
4. 结 论
（1）日本落叶松种子园速效N、速效P缺乏，合理的N、P、K含量配比能有效提高日本落叶松结实母株球果产量，20年以上的日本落叶松无性系种子园N、P、K适宜的配比水平为N1P2K1（1∶0.52∶0.79）。
（2）N1P2K1配比不同水平连年施肥，除施肥当年5.0 kg/株结实母株球果产量最高外，不论是大年还是小年均以2.5 kg/株结实母树球果产量最高，且连续8年结实母株球果累计产量最高。日本落叶松种子园N1P2K1配比水平适宜的施肥量宜选择2.5 kg/株，N、P、K纯养分含量分别是1.08、0.85和0.56 kg/株。
（3）日本落叶松结实母树球果产量受多重因素综合影响，除了施肥效应外，本身遗因素、气象因子、病虫危害等都将影响结实母株球果产量，今后重视施肥管理的同时应加强病虫、气候等灾害的预防和低N低P下N、P高效利用无性系的选择。
-
![]()
图 2 历史火点位置（a）与随机点的位置（b）
Figure 2. Historical fire point location (a) and random point location (b)
![]()
图 7 基于BRT模型全样本结果的边际效应
Figure 7. Marginal effects based on full sample results of BRT model
表 1 变量的简写及其单位
Table 1 Abbreviation of variables and their units
因素
Factor变量
Variable简写
Abbreviation单位
Unit气象因素
Meteorological factor日平均气温 Average daily temperature Temp ℃ 日最高气温 Maximum daily temperature Max_temp ℃ 气温日较差 Daily temperature range Temp_diff ℃ 日均相对湿度 Daily average relative humidity Hum % 日最小相对湿度 Daily minimum relative humidity Minhum % 日均地表温度 Daily average surface temperature G_temp ℃ 日最高地表温度 Daily maximum surface temperature Maxg_temp ℃ 日均降水量 Daily mean precipitation Prec mm 日均风速 Daily mean wind speed Win m/s 日照时数 Sunshine hours Sun h 月均气温 Monthly mean temperature Meantemp ℃ 月均降水量 Monthly mean precipitation Meanprec mm 月均日照时数 Monthly mean sunshine hours Meansun h 月均湿度 Monthly mean humidity Meanhum % 月均地表温度 Monthly mean surface temperature Mg_temp ℃ 春季防火期平均气温 Average temperature during spring fire season TempSpr ℃ 春季防火期平均湿度 Average humidity during spring fire season HumSpr % 春季防火期平均地表温度
Average surface temperature during spring fire seasonG_tempSpr ℃ 春季防火期平均降水量 Average precipitation during spring fire season PrecSpr mm 春季防火期平均日照时数
Average sunshine hours during spring fire seasonSunSpr h 秋季防火期平均气温 Average temperature during fall fire season TempAut ℃ 秋季防火期平均地表温度
Average surface temperature during fall fire seasonG_tempAut ℃ 秋季防火期平均湿度 Average humidity during fall fire season HumAut % 秋季防火期平均降水量 Average precipitation during fall fire season PrecAut mm 秋季防火期平均日照时数
Average sunshine hours during fall fire seasonSunAut h 人为因素
Human factor与最近居民点的距离 Distance to the nearest settlement Dis_res km 与最近道路的距离 Distance to the nearest road Dis_road km 与最近铁路的距离 Distance to the nearest railroad Dis_rail km 与最近瞭望塔的距离 Distance to the nearest watchtower Dis_watch km 植被状况
Vegetation status植被类型 Vegetation type Vegetation type 归一化植被指数 Normalized difference vegetation index NDVI 地形信息
Topographic information海拔 Altitude Dem m 坡向 Aspect Aspect 坡度 Slope Slope ° 社会经济因素
Socio-economic factor人均 GDP GDP per capita GDP 万元/km2
104 CNY/km2人口密度 Population density Pop 人/km2
person/km2注：春季、秋季防火期相关因子指火灾发生前一年的数据。Note: spring and fall fire season weather factors refer to data from the year prior to the fire. 
下载: 导出CSV
表 2 4种模型预测准确性和拟合度的比较
Table 2 Comparison of prediction accuracy and fit of four models
样本
Sample模型
Model最佳临界值
Optimal critical valueAUC值
AUC value预测准确率
Prediction accuracy训练数据集
Training data set测试数据集
Test data set样本1
Sample 1LR 0.378 0.855 0.776 0.842 GR 0.307 0.847 0.772 0.839 RF 0.338 0.949 0.874 0.876 BRT 0.426 0.959 0.895 0.897 样本2
Sample 2LR 0.318 0.854 0.758 0.854 GR 0.352 0.849 0.771 0.852 RF 0.320 0.943 0.867 0.894 BRT 0.431 0.969 0.904 0.921 样本3
Sample 3LR 0.345 0.843 0.760 0.863 GR 0.313 0.850 0.758 0.837 RF 0.320 0.943 0.867 0.876 BRT 0.358 0.970 0.901 0.911 样本4
Sample 4LR 0.361 0.856 0.770 0.857 GR 0.347 0.857 0.774 0.841 RF 0.292 0.946 0.867 0.892 BRT 0.405 0.967 0.903 0.905 样本5
Sample 5LR 0.346 0.858 0.773 0.872 GR 0.350 0.855 0.777 0.857 RF 0.300 0.939 0.866 0.887 BRT 0.388 0.969 0.900 0.900 全样本
Full samplesLR 0.388 0.852 0.755 GR 0.400 0.851 0.769 RF 0.281 0.947 0.863 BRT 0.349 0.967 0.899
下载: 导出CSV
-
[1] Thom D, Seidl R. Natural disturbance impacts on ecosystem services and biodiversity in temperate and boreal forests[J]. Biological Reviews, 2016, 91: 760−781. doi: 10.1111/brv.12193
[2] Pang Y Q, Li Y D, Feng Z K, et al. Forest fire occurrence prediction in China based on machine learning methods[J/OL]. Remote Sensing, 2022, 14: 5546[2024−01−21]. https://doi.org/10.3390/rs14215546.
[3] Mundo I A, Wiegand T, Kanagaraj R, et al. Environmental drivers and spatial dependency in wildfire ignition patterns of northwestern Patagonia[J]. Journal of Environmental Management, 2013, 123: 77−87.
[4] Phelps N, Woolford D G. Comparing calibrated statistical and machine learning methods for wildland fire occurrence prediction: a case study of human-caused fires in Lac La Biche, Alberta, Canada[J]. International Journal of Wildland Fire, 2021, 30: 850−870. doi: 10.1071/WF20139
[5] Nguyen N T, Dang B T N, Pham X C, et al. Spatial pattern assessment of tropical forest fire danger at Thuan Chau area (Vietnam) using GIS-based advanced machine learning algorithms: a comparative study[J]. Ecological Informatics, 2018, 46: 74−85. doi: 10.1016/j.ecoinf.2018.05.009
[6] 曾爱聪, 蔡奇均, 苏漳文, 等. 基于MODIS卫星火点的浙江省林火季节变化及驱动因子[J]. 应用生态学报, 2020, 31(2): 399−406. Zeng A C, Cai Q J, Su Z W, et al. Seasonal variation and driving factors of forest fire in Zhejiang Province, China, based on MODIS satellite hot spots, Chinese[J]. Journal of Applied Ecology, 2020, 31(2): 399−406.
[7] Guris S, Akya E C, Un T. Estimating of probability of home-ownership in rural and urban areas: Logit, Probit and Gompit model[J]. European Journal of Social Sciences, 2011, 21: 405−411.
[8] 苏漳文, 曾爱聪, 蔡奇均, 等. 基于Gompit回归模型的大兴安岭林火预测模型及驱动因子研究[J]. 林业工程学报, 2019, 4(4): 135−142. Su Z W, Zeng A C, Cai Q J, et al. Study on prediction model and driving factors of forest fire in Da Hinggan Mountains using Gompit regression method[J]. Journal of Forestry Engineering, 2019, 4(4): 135−142.
[9] 高超, 林红蕾, 胡海清, 等. 我国林火发生预测模型研究进展[J]. 应用生态学报, 2020, 31(9): 3227−3240. Gao C, Lin H L, Hu H Q, et al. A review of models of forest fire occurrence prediction in China[J]. Chinese Journal of Applied Ecology, 2020, 31(9): 3227−3240.
[10] 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: 117−129. doi: 10.1016/j.foreco.2012.03.003
[11] Eskandari S, Amiri M, Sadhasivam N, et al. Comparison of new individual and hybrid machine learning algorithms for modeling and mapping fire hazard: a supplementary analysis of fire hazard in different counties of Golestan Province in Iran[J]. Natural Hazards, 2020, 104: 305−327. doi: 10.1007/s11069-020-04169-4
[12] Chang C, Chang Y, Xiong Z P, et al. Predicting grassland fire-occurrence probability in Inner Mongolia Autonomous Region, China[J/OL]. Remote Sensing, 2023, 15(12): 2999[2024−01−23]. https://doi.org/10.3390/rs15122999.
[13] Ruffault J, Mouillot F. Contribution of human and biophysical factors to the spatial distribution of forest fire ignitions and large wildfires in a French Mediterranean region[J]. International Journal of Wildland Fire, 2017, 26: 498−508. doi: 10.1071/WF16181
[14] Wu Z W, He H S, Keane R E, et al. Current and future patterns of forest fire occurrence in China[J]. International Journal of Wildland Fire, 2020, 29: 104−119. doi: 10.1071/WF19039
[15] Massada A B, Syphard A D, Stewart S I, et al. Wildfire ignition-distribution modelling: a comparative study in the Huron-Manistee National Forest, Michigan, USA[J]. International Journal of Wildland Fire, 2013, 22: 174−183. doi: 10.1071/WF11178
[16] 刘璞. 内蒙古大兴安岭典型森林积雪化学特征研究[D]. 哈尔滨: 东北林业大学, 2022. Liu P. Study up on chemical characteristics of snow cover in typical forests of Daxing’an Mountains, Inner Mongolia[D]. Harbin: Northeast Forestry University, 2022.
[17] 赵鹏武, 葛晨赫, 王嘉夫, 等. 内蒙古大兴安岭森林火险等级评价与区划研究[J]. 西南林业大学学报(自然科学), 2021, 41(2): 151−158. Zhao P W, Ge C H, Wang J F, et al. Study on rating and classification of forest fire risk in Daxing’an Mountains, Inner Mongolia[J]. Journal of Southwest Forestry University (Natural Sciences), 2021, 41(2): 151−158.
[18] 徐化成. 中国大兴安岭森林[M]. 北京: 科学出版社, 1998: 40−41. Xu H C. Da Hinggan Ling Mountains forests in China[M]. Beijing: Science Press, 1998: 40−41.
[19] 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: 1989−2004. doi: 10.1007/s10980-013-9935-4
[20] Littell J S, McKenzie D, Peterson D L, et al. Climate and wildfire area burned in western U. S. ecoprovinces, 1916−2003[J]. Ecological Applications, 2009, 19: 1003−1021. doi: 10.1890/07-1183.1
[21] Zhang H J, Qi P C, Guo G M. Improvement of fire danger modelling with geographically weighted logistic model[J]. International Journal of Wildland Fire, 2014, 23: 1130−1146. doi: 10.1071/WF13195
[22] 内蒙古自治区统计局. 改革开放30年的内蒙古[M]. 北京: 中国统计出版社, 2008. Inner Mongolia Autonomous Region Bureau of Statistics. Inner Mongolia in the 30 years of reform and opening up[M]. Beijing: China Statistics Press, 2008.
[23] 白海峰, 刘晓东, 牛树奎, 等. 基于贝叶斯模型平均法的森林火灾预测模型构建研究以云南省大理州为例[J]. 北京林业大学学报, 2021, 43(5): 44−52. Bai H F, Liu X D, Niu S K, et al. Construction of forest fire prediction model based on Bayesian model averaging method: taking Dali Prefecture, Yunnan Province of southwestern China as an example[J]. Journal of Beijing Forestry University, 2021, 43(5): 44−52.
[24] Frank E H. Regression modeling strategies-with applications to linear models, logistic regression, and survival analysis[M]. New York: Springer, 2001.
[25] Guo F T, Selvalakshmi S, Lin F F, et al. Geospatial information on geographical and human factors improved anthropogenic fire occurrence modeling in the Chinese boreal forest[J]. Canadian Journal of Forest Research, 2016, 46: 582−594. doi: 10.1139/cjfr-2015-0373
[26] 陈锋, 林向东, 牛树奎, 等. 气候变化对云南省森林火灾的影响[J]. 北京林业大学学报, 2012, 34(6): 7−15. Chen F, Lin X D, Niu S K, et al. Influence of climate change on forest fire in Yunnan Province, southwestern China[J]. Journal of Beijing Forestry University, 2012, 34(6): 7−15.
[27] Wu Z W, He H S, Yang J, et al. Relative effects of climatic and local factors on fire occurrence in boreal forest landscapes of northeastern China[J]. Science of the Total Environment, 2014, 493: 472−480. doi: 10.1016/j.scitotenv.2014.06.011
[28] 舒展. 气候变化对大兴安岭塔河林业局森林火灾的影响研究[D]. 哈尔滨: 东北林业大学, 2011. Shu Z. Impacts research of climate change on forest fires in Tahe Forestry Bureau in Great Xing’an Mountain Region[D]. Harbin: Northeast Forestry University, 2011.
[29] 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.
[30] Guo F T, Zhang L J, Jin S, et al. Modeling anthropogenic fire occurrence in the boreal forest of China using logistic regression and random forests[J/OL]. Forests, 2016, 7: 250[2023−01−23]. https://www.mdpi.com/1999-4907/7/11/250.
[31] Su Z W, Hu H Q, Tigabu M, et al. Geographically weighted negative binomial regression model predicts wildfire occurrence in the Great Xing’an Mountains better than negative binomial model[J/OL]. Forests, 2019, 10: 377[2023−01−23]. http://dx.doi.org/10.3390/f10050377.
[32] 周粉粉, 郭蒙, 钟超, 等. 呼伦贝尔草原火时空格局及特征分析[J]. 地理科学, 2022, 42(10): 1838−1847. Zhou F F, Guo M, Zhong C, et al. Spatial-temporal patterns and characteristics of grassland fire in the Hulunbuir Grassland[J]. Scientia Geographica Sinica, 2022, 42(10): 1838−1847.
[33] 杨少斌, 曹萌, 祝鑫海, 等. 2001—2019年内蒙古大兴安岭北部原始林区森林火灾发生规律研究[J]. 灾害学, 2022, 37(3): 122−128. Yang S B, Cao M, Zhu X H, et al. Research on occurrence law of forest fire in the northern primitive forest area of Greater Khingan Range, Inner Mongolia from 2001 to 2019[J]. Journal of Catastrophology, 2022, 37(3): 122−128.
-
期刊类型引用(5)
1. 吕焕哲，彭渤，孙旭，陈希亚，胡雨萱. 国土空间规划视域下国家公园适宜范围识别研究——以拟建洞庭湖国家公园为例. 湖南师范大学自然科学学报. 2024(03): 37-44 . 
百度学术
2. 邱胜荣，唐小平，刘洋，缐海彦，黄璐，高琦，樊远征. 国家公园范围确定方法及在秦岭国家公园候选区的实践. 国家公园(中英文). 2024(10): 620-631 . 百度学术
3. 曾丽婷 ，戴冰武 ，曾毅 ，吕焕哲 ，谭畅 ，王一淳 . 基于生态安全格局的国家公园边界划定研究——以拟建洞庭湖国家公园为例. 国土资源导刊. 2023(01): 6-13 . 百度学术
4. 孙乔昀，李娜，张玉钧. 面向国土景观风貌管控的中国国家公园空间布局研究. 中国园林. 2022(04): 26-31 . 百度学术
5. 黄晗雯，刘洋，陈君帜，李晓民. 拟建松嫩鹤乡国家公园范围划定方法探讨. 风景园林. 2022(07): 26-29 . 百度学术
其他类型引用(3)
计量
- 文章访问数: 218
- HTML全文浏览量: 44
- PDF下载量: 52
- 被引次数: 8
