基于Landsat时序影像和LandTrendr算法的森林保护区植被扰动研究——以陕西柴松和太白山保护区为例

殷崎栋,柳彩霞,田 野

1 陕西生态环境规划设计院有限公司, 西安 710075 2 中国科学院空天信息创新研究院遥感科学国家重点实验室, 北京 100101 3 中国科学院城市环境研究所城市环境与健康重点实验室, 厦门 361021

加强和发展自然保护区事业是当前我国生态环境保护工作的一项重要而紧迫的任务。但随着我国经济的快速发展和对自然资源的需求日益增加,保护与开发的矛盾日益增加,“人工化”、“生境退化”的趋势日益严重,自然保护区面临着前所未有的压力[1]。

遥感技术的发展为自然保护区保护提供了切实可行的监测和保护手段,特别是陆地卫星Landsat提供30m空间分辨率的定期图像采集,平均重访周期为16d,可免费获得并追溯到1980年代早期,可支持植被和土地覆盖变化的长期观测和监测[2]。利用经过大气校正和地形校正的地表反射率产品,确保多个时间序列的图像具有可比性[3],使Landsat图像在长时间序列地表覆盖变化分析方面具有优越性。在长时间序列遥感数据的支持下,可以在不同的应用领域进行生态环境和土地变化的评估、追溯和历史重构,如植被监测[4-5]、土地利用变化[6-7]、环境扰动[8- 10]等。在国内自然保护区遥感监测方面,学者们已就崂山[11]、可可西里[12]、十八里峡[13]、佛坪[14]等自然保护区的人类活动进行了遥感监测。然而,先前的研究大部分集中在时间序列图像的离散信息提取,年际变化研究存在年份的空白,无法捕捉连续的地表变化特征,从而限制了数据的长期监测和生态环境稳定性评估能力。

Landsat时间序列数据可以监测森林扰动,但是需要合理处理因大气和几何校正、植被物候、太阳高度变化和传感器退化引起的数据噪声。Kennedy等[15]发现,在美国俄勒冈州西部,不同的光谱指数在准确检测细微干扰信号方面的能力各不相同。归一化植被指数(NDVI)的效率低于归一化燃烧比(NBR)或缨帽变换湿度(TCW),在森林监测领域,很多算法使用有限的光谱波段或指数来检测干扰,例如NBR[15]、森林指数(FI)[16]、NDVI和SWIR-NIR比值[17]或缨帽变换角度[18],因此,在检测森林扰动时仔细考虑比较光谱指标是非常重要的步骤。

目前,已经有多种基于Landsat时间序列的变化检测技术。例如LandTrendr[15]、BFAST[19-20]、CCDC[21-22]。由于时间序列算法的数据预处理和算法本身都非常耗时,为了达到更广泛的研究应用群体,LandTrendr算法得以在Google Earth Engine (GEE)平台上实施[23],从而提高了算法计算速度。与以前在本地服务器上运行变化检测算法相比,GEE中的LandTrendr利用云计算的优势更有效地检测森林的变化,为森林自然保护区的保护提供了大规模运用的机会。

了解森林保护区生态系统的组成、结构和功能的现状和长期趋势发展,可以为保护区管理提供科学可信的本底资料。陕西省地处黄河中游、长江上中游地区,属于国家生态环境建设重点区域。全省共有45个省级以上森林自然保护区和湿地自然保护区,其中国家级19个,省级26个。为了检验LandTrendr算法对不同生态气候分区的森林植被类型的敏感性是否存在差异,我们选取陕西省柴松和太白山自然保护区,利用GEE云计算平台中的LandTrendr算法进行森林变化检测分析,追踪近40年来陕西省柴松和太白山自然保护区的森林扰动。主要回答以下问题：(1)柴松省级自然保护区和太白山国家级自然保护区自保护区成立以来植被有明显的变化趋势吗？(2)如果保护区内植被发生了明显变化,变化的空间分布如何？(3)发生变化的原因是什么？

1 研究区概况

陕西延安柴松省级自然保护区(东经108°39′—108°50′,北纬35°56′—36°07′)位于延安市富县西部(图1)。保护区成立于2004年。境内南北长约18 km,东西宽约16.5 km,海拔1000—1582 m之间,总面积17640 hm2(表1)。保护区植被属华北暖温带半湿润落叶阔叶林地带的北部落叶阔叶地带的西缘。森林植被均为天然次生林,成为黄土高原上的植物种类独特、多样、丰富的天然绿色“基因库”,其中价值较高的有柴松、油松、侧柏、栎类和山杨等86种。主要保护对象为柴松林及其森林生态系统。

陕西太白山国家级自然保护区(东经107°22′25″—107°51′30″和北纬33°49′30″—34°05′35″)位于陕西省太白、眉县、周至三县交界处(图1),保护区成立于1965年,1986年7月晋升为国家级自然保护区,占地面积56325 hm2(表1),主要保护对象为森林生态系统和自然历史遗迹。保护区地处秦岭山脉中段,太白山顶海拔3767 m,为秦岭的主峰。秦岭山脉为中国南北气候和植物的分界线,太白山则处于华北、华中和青藏高原三区生物交汇过渡地带,区内动植物资源丰富,植被垂直分带明显。高等植物有2000余种,国家重点保护植物有连香树、水青树、星叶草、太白红杉等21种。高等动物有270多种,国家保护动物有大熊猫、羚牛、豹等20多种。在太白山海拔3000 m以上的高山区,还保存着比较完整的第四纪冰川遗迹。该保护区是开展科学研究难得的“天然实验室”。

图1 研究区位置意图Fig.1 Location map of the study area

表1 研究区概况

2 数据和方法

2.1 数据来源与预处理

利用GEE平台自1984—2018年期间的全部存档Landsat TM/ETM/OLI的地表反射率数据构建时间序列数据集。由于蓝光波段对大气状况敏感,只选取了5个波段：绿(green)、红(red)、近红外(nir)和2个短波红外波段(swir1和swir2)。为了尽量减少因物候或太阳高度角变化引起的反射率变化,选取每年6月10日到8月30日的影像。由于Landsat OLI比TM/ETM+具有更高的辐射分辨率(12比特),我们利用了这两个传感器光谱值之间的统计协调函数对反射率进行归一化[24],接着进行提取年集合的中位反射率值,构建了具有最小云量图像的35年合成图,进而计算归一化植被指数NDVI[25](公式(1))和归一化燃烧比值NBR[26](公式(2))。通过比较不同比值指数和波段对森林扰动的敏感性(图2),NDVI和NBR对柴松保护区的森林扰动识别效果更好,单一波段受数据噪声的影响较大。综合考虑以往研究经验,本文利用NBR指数探测森林扰动。

(1)

(2)

图2 不同植被指数和波段对森林扰动的敏感性Fig.2 Sensitive of vegetation index and bands for forest disturbance NBR:归一化燃烧比值, Normalized Burn Ratio; NDVI:归一化植被指数, Normalized Difference Vegetation Index; swir2:短波红外波段2；swir1:短波红外波段1

2.2 LandTrendr算法探测森林扰动

利用Kennedy, Yang[15]开发的LandTrendr算法来检测保护区内森林的变化。LandTrendr算法的核心是一种时间分割算法,用于捕捉时间序列中长期、渐进或者短期的剧烈变化,可以分析每个像元的时间序列谱线来监测其是否发生变化。每个像素的输入是一个光谱指数或者波段的时间序列。通过设置去除噪声引起的尖峰(spikeThreshold)、识别潜在的断点(vertexCountOvershoot)、拟合轨迹和最佳分段数(pvalThreshold),寻找最佳模型(bestModelProportion)。因此,LandTrendr要求设定控制参数以保证变化检测的质量,通过实验分析,本研究对LandTrendr的参数如表2所示。在这里需要说明的是,基于保护区建立时间的考虑,对于柴松保护区,研究集中在2004—2018年,总计15a。对于太白山保护区,研究阶段为1984—2014年(总计31a)。因此对柴松保护区而言,研究时段相对较短,在LandTrendr的参数设置中,分割单元的数目可以稍少,maxSegments设置为4。根据以往研究经验,2个保护区的其他参数设置一致。

表2 LandTrendr参数

对森林扰动发生的空间和时间分布检测对森林保护和政策管理非常重要。相比传统的土地覆盖分类和变化检测技术,LandTrendr可以成功捕获森林扰动信号进行空间和时间制图。Landsat系列卫星能提供长达近40年的卫星观测数据,使对森林保护区的长期观测和监测成为可能。尽管已有很多利用LandTrendr或者相似时间序列检测算法在本地服务器的使用案例[15, 27- 29],算法实施带来的数据准备和运行带来巨大的时间成本和硬件成本。但是,在GEE支持下的LandTrendr算法极大的提高了算法效率。用户不再需要从美国地质调查局或者其他卫星地面站下载大量的原始有云卫星图像,而是基于Google云平台合成用户感兴趣区域的无云或者少云高质量拼接图像,进而利用LandTrendr的主要作者和GEE合作开发的算法进行森林扰动探测。用户可以通过Java script编程加载感兴趣结果图层和时间序列卫星图像,在GEE平台中对算法结果进一步目视判读和参数调整,最后只需要导出结果图层,如扰动年份等。

2.3 精度评价

为了准确评估森林扰动变化的准确性,在研究时间内,在柴松保护区内选择了100个像元表示变化,另外100个像元代表不变区域。对每个像元,通过人工目视判读合成的15幅图像,并使用高空间分辨率的Google Earth(GE)图像对土地覆盖进一步确认,记录该像元是否发生了变化,何时(发生年份)发生了变化,从而组成200个像元的验证样本集,生成混淆矩阵,计算用户精度、生产者精度和总体精度。

3 结果

3.1 LandTrendr算法的精度评价

LandTrendr的时间分割和拟合可以捕捉森林向非森林转变的光谱特征。表3所示的混淆矩阵说明,生产者精度和用户精度均达到85%以上,总体精度达到93%。对于检测出森林变化的像元,进一步评估了变化年份。表4给出了变化年份的精度,2007年之后的变化年份生产精度都达到了80%以上。用户精度除2007年和2005年稍低以外(分别是70%和75%, 其他年份的用户精度达100%,变化年份检测的总体精度为89%。

表4 森林扰动年份精度评估

3.2 森林扰动的历史空间分布

图3所示的扰动图显示了2004—2018年间柴松保护区森林扰动的空间分布和变化年份。森林扰动发生的区域主要集中在柴松保护区的实验区和缓冲区,但是核心区也存在部分森林扰动。扰动类型主要为人类活动所致。总体上,2004—2018年间,共有42.74 hm2的森林发生了扰动(表5),转化为非森林,共计251个斑块,年变化最大面积值和斑块数均发生在2008年,2007、2015年和2005年次之,森林向非森林转化发生的斑块数依次为57、24和20块(图3)。

图3 柴松保护区扰动空间和历史分布Fig.3 Disturbance history of Chaisong Reserve(a. 扰动年分布； b. 实验区扰动； c. 缓冲区和核心区扰动; d-e. Landsat图像(NIR/R/G); f-h. GE高清卫星影像；P1位置[108.723276, 36.082742]; P2位置[108.697087, 36.024953])

表5 柴松保护区森林扰动图斑统计

图4所示的扰动图显示了1984—2018年间太白山保护区森林扰动的空间分布和变化年份。森林扰动发生的区域主要集中在太白山保护区核心区。扰动主要由自然因素引起,如第四纪冰川融化形成的侵蚀沟和其他自然因素造成的森林退化。1984—2018年间,共有23.68 hm2的森林发生了扰动(表6),转化为非森林,共计57个斑块数,年变化最大面积值和斑块数均发生在2013年,1994年次之(7块,1.87 hm2)(图4)。

4 讨论

4.1 森林扰动原因

本文结果表明,LandTrendr可以捕获不同气候区森林扰动的空间分布特征。在降雨量偏少的柴松保护区,扰动最为频繁的年份集中在2007—2008年间,这是由于保护区内违法建设了一批以石油开采为目的的钻井平台及附属设施(图3),从而造成了森林保护区森林覆盖的破碎化。自2016年之后森林扰动有减少趋势。以环保部为执行主体的自然保护区“绿盾”行动自2017年在国家级保护区实施以来,对省市级自然保护区的保护和监管起到了督促作用,在保护区内的部分违法违规设施正在逐步清退。

图4 太白山保护区扰动空间和历史分布Fig.4 Disturbance history of Taibaishan Reserve Disturbance yeara. 扰动年分布； b. 冰川消融扰动； c. 其他自然扰动; d-f. Landsat图像(NIR/R/G); g-i. GE高清卫星影像；P1位置[107.79248, 33.94858]; P2位置[107.687327, 33.966819]

表6 太白山保护区森林扰动图斑统计

在降雨量偏多的太白山保护区,森林扰动活动远远少于柴松保护区,扰动最为频繁的年份集中在2013年,且扰动活动与气候变化和冰川融化造成的地形地貌变化相关(图4)。长远来看,冰川活动的变化,特别是冰川融化导致的裸露冰碛石移动,会影响到冰川末端的高山森林带,可能会改变高寒地区森林的分布[30-31]。因此,气候变化对自然保护区的生态环境保护提出了更多挑战[32],进行长时间序列的遥感监测能为处于高寒地带人类难以到达的森林生态系统类保护区提供可参考依据。

4.2 未来工作

进一步确定自然保护区内森林退化的诱导因素,哪些是人类活动,特别是非法人类活动造成的,如违章建筑、违章工矿设施。哪些是自然因素造成的,如森林病虫害、气候变化等。本研究涉及的柴松自然保护区的关键扰动因素是因发展经济驱动的人类活动,有少量活动甚至是在保护区的核心区。虽然造成森林扰动的土地面积与整个保护区的面积相比,占比并不大,但是森林扰动破坏了保护区动物栖息地的连通性,造成了自然景观的破碎化,开采活动造成的噪声、粉尘及其他行为对动植物的影响还需进一步评估。根据我国《自然保护区条例》规定：自然保护区内保存完好的天然状态的生态系统以及珍稀、濒危动植物的集中分布地,应当划为核心区,禁止任何单位和个人进入；除依照本条例第二十七条的规定经批准外,也不允许进入从事科学研究活动；核心区外围可以划定一定面积的缓冲区,只准进入从事科学研究观测活动。缓冲区外围划为实验区,可以进入从事科学试验、教学实习、参观考察、旅游以及驯化、繁殖珍稀、濒危野生动植物等活动。未来我们需要确定核心区、缓冲区、实验区的扰动范围和时间,结合动植物多样性观测数据(如红外相机捕捉的野生动物活动区域),开展保护区森林扰动对动物栖息地的影响研究[33-34]。

5 结论

本文基于光谱轨迹变化检测算法对森林保护区森林扰动进行了评估。利用GEE云计算平台的LandTrendr算法和Landsat影像,对陕西省柴松自然保护区和太白山自然保护区的森林向非森林转化进行了检测。对变化像元和稳定像元的总体检测精度达到93%,对扰动年份的总体检测精度达到89%。柴松保护区的扰动年份主要发生在2008年,由人类活动所致。太白山自然保护区的扰动年份主要发生在2013年,由自然因素所致。柴松保护区扰动面积为42.74 hm2,太白山保护区的扰动面积为 23.68 hm2,柴松保护区的扰动斑块数目也明显大于太白山保护区。本文研究结果可以帮助科研人员和管理人员了解当前两个保护区的森林生长状态,同时为自然保护区生态环境保护评估提供基线信息。