宁夏南部的六盘山、中部的罗山、西北部的贺兰山3个国家级自然保护区(简称“三山”)，作为绿色生态屏障和水源涵养林区，具有维护生物多样性、保持水土、防风固沙等生态功能，是我国西北重要的生态安全屏障，也是黄河流域生态保护和高质量发展国家战略的重要组成单元。森林的动态变化一定程度上影响着森林生态系统的稳定和生态功能的发挥，而森林的动态变化包括森林突发性变化和植被时序趋势变化^[1]。森林干扰监测属于研究森林突发性变化，主要表现为森林火灾、森林采伐以及因受雨雪和病虫灾害等所引起的森林改变^[2]。近30 a来，随着宁夏经济社会快速发展，森林干扰事件也逐渐增加。因此，开展“三山”长时间序列的森林干扰监测，对于守好“三山”自然保护区，维护西北乃至黄河流域生态安全具有重要的理论和实践意义。

传统的森林干扰监测主要采用实地测量和遥感影像人工识别的方法，虽然监测精度较高，但是工作量大、成本高且时效性差；因而发展高精度、自动化的基于长时序、高分辨率遥感影像的森林干扰监测模型显得极为迫切^[3]。近年来，随着美国地质调查局(United States Geological Survey, USGS)免费开放Landsat堆栈数据，基于Landsat影像的森林干扰监测得到长足的发展。Landsat提供30 m空间分辨率、16 d时间分辨率的多波段影像集，并可免费追溯到1984年^[4]，使其成为主要的长时序森林干扰监测的遥感数据源。基于Landsat时间序列堆栈的森林干扰监测算法有很多，如VCT(vegetation change tracker)、CCDC(continuous change detection and classification)、BFAST(breaks for additive and seasonal trend)、Vogelmann、Brooks、LandTrendr等^[5]，目前比较流行的是BFAST和LandTrendr算法。BFAST算法的优点是不需要设置阈值或变化轨迹，而是将时间序列分成季节、趋势和残差3个组分，可减少背景噪声和季节变化对监测结果的影响，可用于多源遥感数据；缺点是该算法无法检测出再次干扰过程，同时模型也较为复杂^[6]。LandTrendr算法的优点是可同时监测出森林干扰和恢复发生的时间、地点、面积和强度，可通过设置光谱时间轨迹的参数减少时间分割过程中出现的过拟合问题和背景噪声，适用于同源遥感数据；缺点是参数阈值设置较为复杂，数据运算量大^[7-8]。2018年Kennedy等^[9]在GEE(Google Earth Engine)平台上实现LandTrendr算法的应用，有效解决该算法需要大量的数据管理、预处理以及计算资源的工作要求。GEE是目前世界上最先进、处理能力最强的遥感云平台。GEE结合LandTrendr已成为大区域、长时序森林干扰监测的最有利工具^[10]。此外，通过查阅国内外文献，未见关于宁夏“三山”森林干扰研究的报道。

因此，笔者将基于GEE遥感平台中的Landsat 5—8时间序列堆栈(1990—2020年)，通过LandTrendr算法监测近30 a来宁夏“三山”的森林干扰，包括干扰年份、干扰面积、干扰持续时间等，并借谷歌(Google)现有影像以及历史数据，对监测结果进行精度评价，最后统计分析“三山”森林干扰的时空分布特征，以期探讨宁夏“三山”的森林干扰程度和趋势是否存在差异，以及造成“三山”森林干扰的主要原因。

宁夏“三山”包括六盘山、罗山、贺兰山3个国家级自然保护区(图 1)，分别成立于1988、2002和1988年。六盘山(E 106°09′~106°30′，N 35°15′~35°41′)，总面积为6.94万hm²，主要森林植被为华北落叶松(Larix principis-rupprechtii)、油松(Pinus tabulaeformis)、华山松(Pinus armandii)、山杨(Populus davidiana)、白桦(Betula platyphylla)、辽东栎(Quercus liaotungensis)、槭树(Acer miyabei)和山柳(Salix pseudotangii)等^[11]。罗山(E 106°04′~106°24′，N 37°11′~37°25′)，总面积为3.44万hm²，主要森林植被为青海云杉(Picea crassifolia)、油松、山杨、白桦等^[12]。贺兰山(E 105°49′~106°41′，N 38°19′~39°22′)，总面积为19.24万hm²，主要森林植被为青海云杉、灰榆(Ulmus glaucescens)、油松、山杨、白桦等^[13]。

图 1(Fig. 1)

图 1 宁夏“三山”位置 Fig. 1 Location of "Three Mountains" of Ningxia

利用GEE遥感云平台(https://code.earthengine.google.com) 提供的1990—2020年期间的T1级别(质量最高)的Landsat 5(TM)、Landsat 7(ETM+) 和Landsat 8(OLI/TIRS)地表反射率数据产品，构建时间序列堆栈(LTSS)。为尽量减少因物候或太阳高度角变化引起的反射率变化，笔者选取每年6—9月的影像，并对时间序列堆栈中的水体、云和云阴影进行掩膜处理，最后采用中值法重构目标年份最小云量的季节合成影像。

笔者选用归一化燃烧率(normalized burn ratio，NBR)构建时间序列数据进行森林干扰监测。研究^{[1, 14]}表明，NBR对多种类型森林干扰事件具有最大的敏感性，比采用归一化植被指数(normalized difference vegetation index, NDVI)、缨帽变换湿度指数(tasseled cap transformation wetness, TCW)，更能降低噪声的干扰，错误率也最低。

式中：N为归一化燃烧率，量纲为1；P₁和P₂分别为近红外(850~880 nm)和短波红外(1 570~1 650 nm)地表反射率。

LandTrendr时间序列分割算法(简称LandTrendr算法)^{[7, 9]}有6个步骤：1)提取时间序列轨迹；2)识别潜在分割点；3)用角度变化删除多余分割点；4)基于点与点连接或回归的方法拟合分割点；5)创建连续简化的拟合轨迹；6)利用p值选择最佳拟合轨迹。表 1是笔者所用到的参数及取值。

表 1(Tab. 1)

表 1 LandTrendr算法参数设置 Tab. 1 Parameters used in LandTrendr algorithm 参数名称Parameter 笔者取值Value of this article 数据类型Data type 顶点阈值Spike threshold 0.9 浮点型Float 最大分段数Max segments 6 整型Integer 可超过的顶点数量Vertex count overshoot 3 整型Integer 拟合的最大P值Fitted max P 0.05 浮点型Float 最佳模型比Optimal model proportion 0.75 浮点型Float 恢复阈值Recovery threshold 0.25 浮点型Float 防止一年恢复Prevent one year recovery TRUE 布尔型Boolean 需要的最小观测数Min observations needed 6 整型Integer

表 1 LandTrendr算法参数设置 Tab. 1 Parameters used in LandTrendr algorithm

森林干扰制图前需先对森林分布范围进行掩膜，以排除水体、裸地等非森林地物。笔者采用计算像元的缨帽变换角(tasseled cap transformation angle, TCA)，并通过设定阈值提取森林像元。提取的像元在整个时序中被识别为持续森林^[1]，计算方法如下：

式中：T为缨帽变换角；T₁为缨帽变换湿度分量；T₂为缨帽变换绿度分量。

为准确评估森林干扰的准确性，在研究时段和研究区内选择森林干扰监测变化像元和不变像元各90个，利用Google Earth Pro的高时空分辨率影像对这180个像元进行人工目视解译(记录像元是否发生变化和变化时间)，并将目视解译结果与采用LandTrendr算法干扰监测结果相对比，最终得到监测精度的混淆矩阵，然后根据混淆矩阵再计算出总体精度、生产者精度、用户精度、Kappa系数等^[5]。

如表 2所示，LandTrendr算法的用户精度、生产者精度、总体精度均能＞90%，Kappa系数达到0.87，麦克尼马尔检验结果表明，人工目视解译与LandTrendr森林干扰监测结果无显著差异(P=0.39＞0.05)。

表 2(Tab. 2)

表 2 LandTrendr算法精度评价 Tab. 2 Accuracy assessment for LandTrendr algorithm 类别 Class 变化像元 Changed pixels 稳定像元 Stable pixels 像元总数 Number of total pixels 用户精度 User accuracy/% 变化像元Changed pixels 82 8 90 91.11 稳定像元Stable pixels 4 86 90 95.56 像元总数Total number of pixels 86 94 180 生产者精度Producer accuracy/% 95.35 91.49 总体精度Overall accuracy/% 93.33 Kappa系数Kappa coefficient 0.87 麦克尼马尔检验McNemar's test 0.39＞0.05

表 2 LandTrendr算法精度评价 Tab. 2 Accuracy assessment for LandTrendr algorithm

LandTrendr算法的森林干扰监测结果见表 3。按照森林干扰的面积排序，1991—2020年“三山”森林干扰总面积：六盘山(1 231.40 hm²)＞贺兰山(90.24 hm²)＞罗山(23.87 hm²)，六盘山是贺兰山13.65倍，是罗山的51.59倍，贺兰山是罗山的3.78倍。按照森林干扰的总面积与保护区的面积比排序，1991—2020年“三山”森林干扰总面积比：六盘山(1.77%)＞罗山(0.07%)＞贺兰山(0.05%)，六盘山是罗山25.29倍，是贺兰山的35.40倍，罗山是贺兰山的1.40倍。如表 3所示，六盘山森林干扰面积最大值(511.72 hm²)发生在1997年；罗山森林干扰面积最大值(7.65 hm²)发生在1999年；贺兰山森林干扰面积最大值(32.29 hm²)发生在1991年。森林干扰面积标准差为：六盘山(97.75 hm²)＞贺兰山(7.45 hm²)＞罗山(1.98 hm²)，表明六盘山森林干扰面积变化波动较大，其次为贺兰山和罗山。

表 3(Tab. 3)

表 3 “三山”干扰面积统计 Tab. 3 Disturbance area statistics of "Three Mountains"  hm2 六盘山(总面积69 419.76 hm2) Liupan Mountain (total area 69 419.76 hm2) 罗山(总面积34 395.49 hm2) Luo Mountain (total area 34 395.49 hm2) 贺兰山(总面积192 405.45 hm2) Helan Mountain (total area 192 405.45 hm2) 年份 Year 干扰面积 Disturbance area 年份 Year 干扰面积 Disturbance area 年份 Year 干扰面积 Disturbance area 1991 137.83 1991 0 1991 32.29 1992 6.58 1992 0 1992 0 1993 0 1993 0 1993 26.38 1994 31.36 1994 0 1994 0.77 1995 171.91 1995 0 1995 1.33 1996 92.38 1996 0 1996 0 1997 511.72 1997 0.93 1997 0 1998 0.88 1998 0 1998 0.45 1999 1.68 1999 7.65 1999 6.86 2000 34.90 2000 2.86 2000 0 2001 3.24 2001 0 2001 0 2002 10.25 2002 0 2002 0 2003 2.04 2003 0 2003 0 2004 7.67 2004 2.07 2004 1.47 2005 32.04 2005 7.00 2005 0 2006 12.89 2006 0 2006 0 2007 12.76 2007 0 2007 0 2008 3.58 2008 0 2008 4.05 2009 11.34 2009 0 2009 0 2010 37.03 2010 0 2010 0.27 2011 4.46 2011 0 2011 0.77 2012 3.67 2012 0 2012 0 2013 6.06 2013 0 2013 0 2014 1.24 2014 0 2014 6.60 2015 14.69 2015 0 2015 1.50 2016 6.41 2016 0 2016 0 2017 0.05 2017 0 2017 0 2018 16.47 2018 0 2018 3.20 2019 3.44 2019 0 2019 0 2020 52.83 2020 3.36 2020 4.30 合计Total 1 231.40 合计Total 23.87 合计Total 90.24 平均Mean 41.05 平均Mean 0.80 平均Mean 3.01 标准差 Standard deviation 97.75 标准差 Standard deviation 1.98 标准差 Standard deviation 7.45

表 3 “三山”干扰面积统计 Tab. 3 Disturbance area statistics of "Three Mountains"

图 2a~c反映的是“三山”森林年干扰面积与保护区的面积比(年干扰面积比)随时间的变化特征，总体上都呈现出不规律的波动。但是，“三山”的年干扰面积比都有减少的趋势(斜率)，六盘山减少较为明显，平均每年减少52.60×10^-4%(3.65 hm²)；其次为贺兰山，平均每年减少1.58×10^-4%(0.30 hm²)；最少的为罗山，平均每年减少0.40×10^-4%(0.014 hm²)。1991—2020年，六盘山森林干扰面积减少15.78×10^-2%(109.54 hm²)，贺兰山减少0.47×10^-2%(9.12 hm²)，罗山减少0.12×10^-2%(0.42 hm²)。六盘山是贺兰山33.57倍，是罗山的131.50倍，贺兰山是罗山的3.92倍。图 2d是“三山”年均降水量随时间的变化特征，按照变异系数排序：贺兰山(27.4%)＞罗山(25.0%)＞六盘山(21.1%)。

图 2(Fig. 2)

图 2 年干扰面积比和年均降水量随时间变化特征 Fig. 2 Change characteristics of annual disturbance area ratio and average annual precipitation with time

图 3反映1991—2020年“三山”森林干扰的空间分布情况和干扰的持续时间与年份。“三山”中，六盘山的森林干扰主要发生在六盘山的北部和东侧2块“飞地”，1991—2020年间，六盘山的森林干扰总面积为1 231.40 hm²，2块“飞地”所占比例就达到63.63%(783.55 hm²)。罗山的森林干扰主要发生在罗山保护区的北部，俗称“大罗山”，由于历史原因大罗山建有省级电视调频转播台和地震监测台，设施的扩建和人为活动较为频繁；而“小罗山”(罗山保护区的南部)和“飞地”几乎没有森林干扰。贺兰山的森林干扰主要发生贺兰山的中部。这一地区连续分布有贺兰山岩画、苏峪口国家森林公园、拜寺口双塔、滚钟口、西夏王陵等重要的风景名胜区，旅游开发强度较大，人类干扰频繁。

图 3(Fig. 3)

图 3 “三山”森林干扰年份和持续时间空间格局 Fig. 3 Spatial pattern of forest-disturbed year and duration of "Three Mountains"

如图 3所示，“三山”的森林干扰持续时间都在1~3 a之间，且多集中在1 a。六盘山92.62%的森林干扰持续时间为1 a，5.54%的干扰持续时间为2 a，1.84%的干扰持续时间为3 a；罗山90.71%的干扰持续时间为1 a，8.99%的干扰持续时间为2 a，0.3%的干扰持续时间为3 a；贺兰山57.88%的干扰持续时间为1 a，5.71%的干扰持续时间为2 a，36.41%的干扰持续时间为3 a。

根据1991—2020年“三山”森林干扰总面积分析结果，六盘山的干扰总面积最高，其次为贺兰山和罗山(表 3)，考虑到“三山”占地面积大小不同，该文计算了“三山”森林干扰总面积比，其排序结果为：六盘山(1.77%)＞罗山(0.07%)＞贺兰山(0.05%)，可见六盘山的干扰程度要远高于罗山(25.29倍)和贺兰山(35.40倍)。根据文献^[15]记载，1991—2020年间，六盘山仅在2007年5月16日的小南川林区发生过火灾，过火面积9.59 hm²，该文也监测到了这次森林火灾(表 3)。因此，森林火灾并不是造成六盘山干扰程度较高的主要原因。据统计，长期生活在六盘山自然保护区范围内的农户约3 500人，分布在泾源、隆德2县6个乡(镇)约26个行政村，开荒垦殖、伐木、放牧、挖药、薪炭等多种人为活动的影响，都有可能造成森林干扰事件的发生^[16]。从“三山”森林干扰的空间分布(图 3)可见，六盘山的森林干扰主要发生在北部和东侧2块“飞地”，而这些区域也是六盘山人类活动较为集中的区域^[17]。此外，由图 3分析可知，六盘山森林干扰持续时间为1 a的面积比例最大，达92.62%，扰持续时间为2和3 a的分别只有5.54%和1.84%，结合Google影像目视解译，六盘山森林干扰主要是人类活动造成的急剧干扰事件，而不是森林病虫害，因为病虫害对森林的干扰持续时间会较长^[18]。

1991—2020年，罗山的森林干扰总面积比只有0.07%，且主要发生在罗山保护区北部的“大罗山”(图 3)，由于历史原因，大罗山建有省级电视调频转播台和地震监测台等设施，山区道路的修建和建筑设施的扩建以及人为活动都较为频繁^[19]。此外，罗山90.71%的干扰持续时间为1 a，因此，可以推断罗山森林干扰的主要原因也是人类活动。与六盘山和罗山相比，虽然贺兰山干扰持续时间为1 a的面积比也超过50%(57.88%)，但是持续时间为3 a的面积比却达到36.41%，可见贺兰山森林干扰的原因与六盘山和罗山不完全相同，应是人为干扰和自然干扰共同作用的结果。人为干扰的原因是贺兰山分布有众多著名的风景名胜区，旅游开发强度较大，人类干扰频繁^[20]。自然干扰的原因可能和贺兰山所处地理位置的气候条件有关。贺兰山地处我国温带草原区与荒漠区的过渡地带，气候干旱少雨，植被盖度低，当遇到降水减少、气温升高、湿度下降等不利于森林植被正常生长的年份，就有可能造成森林干扰事件的发生^[21]。根据对“三山”年均降水量的分析，贺兰山的降水量明显低于六盘山和罗山，且贺兰山降水变异还要高于六盘山和罗山(图 2d)，这就导致贺兰山的森林干扰受自然因素的影响较大^[21]；由图 2d中贺兰山的年均降水量变化可见，1991—1993年贺兰山降水量持续降低，并达到30 a来的最低值，与之对应，干扰持续时间为3 a的区域也达到最大；此外，通过查阅相关文献，并未发现贺兰山近30 a有大规模病虫害的发生。综合分析，贺兰山的森林干扰是人为活动和自然条件共同作用的结果。

1) 1991—2020年，宁夏“三山”森林干扰主要发生在六盘山的北部和东侧2块“飞地”、罗山的北部即“大罗山”和贺兰山的中部。“三山”的森林干扰持续时间都在1~3 a内，且多集中在1 a。

2) 宁夏“三山”中，六盘山的干扰程度(1.77%)要远高于罗山(0.07%)和贺兰山(0.05%)，“三山”森林干扰的原因不完全相同，六盘山和罗山森林干扰主要是人为活动造成，而贺兰山森林干扰是人为活动和自然条件共同作用的结果。

3) 宁夏“三山”的森林干扰面积比每年都有减少的趋势，六盘山减少最为明显(52.60×10^-4%/a)，其次为贺兰山(1.58×10^-4%/a)和罗山(0.40×10^-4%/a)。基于GEE平台的LandTrendr时间序列分割算法，可及时、准确地获取森林干扰信息。