基于云平台和BFAST算法的地表变化检测方法

周 旭 陈元鹏 刘岩涛 周 妍 李少帅 王 力

(1.自然资源部国土整治中心， 北京 100035； 2.中国地质大学(北京)土地科学技术学院， 北京 100083；3.中国科学院空天信息创新研究院遥感科学国家重点实验室， 北京 100101)

0 引言

受自然变迁和人类活动影响，地表自然环境和生态系统不断发生变化，分析这些变化信息是地理科学的一个重要研究方向[1-3]。遥感技术的发展增强了对地观测的能力，提供了重要的支撑数据，更便于准确掌握地表变化信息。迄今为止，行业领域已积累了大量遥感数据(如Landsat、MODIS、AVHRR等)，这些遥感数据以时间序列的方式记录了地表的动态变化，通过对变化情况的研究，能够在时间和空间上检测到地表的扰动，探究地表自然环境和生态系统发展演变的规律，支撑相关的科研与行政管理工作[4-7]。

基于遥感时序数据的地表变化检测属于时序变化分析研究，从遥感时序数据中开展变化分析存在一定困难，因为遥感时序数据中不仅包含季相性变化、趋势性变化和突发性变化，还存在数据处理过程中产生的大气散射、云影效应和几何误差等。因此，在遥感时序变化分析中，如何在季相性变化中即不受到误差干扰又能够有效判定和识别变化是核心问题。

目前，在遥感时序变化分析中，已存在多类方法，如阈值、差分、轨迹分割、分类、统计和回归[8-10]，每类方法均包含6种特征，分别为时序数据频率(影像数量)、光谱指数类型、变量维度(单变量/多变量)、实时性(非实时/近实时)、变点类型(突变/渐变)、检测单元(像元/亚像元)。但各类方法均存在一定弊端[11-12]。

从具体算法方面看，现阶段应用较为广泛的算法主要包括，如LandTrendr(Landsat-based detection of Trends)是轨迹分割类算法，采用单变量光谱指数NBR，检测地表变化，属于非实时检测；VCT(Vegetation change tracker)是阈值类算法，采用单变量光谱指数IFZ，识别地表变化，属于非实时检测；CCDC(Continuous change detection and classification)是统计和分类结合的算法，采用多变量光谱和指数，识别土地覆被类型的变化，属于近实时检测。但以上方法各有优劣，如LandTrendr、VCT在时序数据频率使用上相对较低(如年度均值)，一定程度上弱化了对地表连续变化的检测能力，并且两种算法属于非实时的检测，这在应对稀疏数据时会产生较大的结果延迟性；CCDC采用了多变量输入，而多变量的输入方法本身就会对分类结果带来误差，加之时序的影响因素，将进一步影响检测结果的准确性。相比之下，BFAST(Breaks for additive season and trend)采用单变量光谱指数能够相对减小检测结果的误差，采用高时序数据频率能够增强对地表连续变化的检测能力和响应速度，利于更精准、近实时地检测地表扰动变化，因此，在遥感时序变化分析方面更具有优势。

此外，除模型算法方面存在的问题外，还存在其他因素会影响基于遥感时序数据的地表变化检测，如受天气、卫星重访周期等因素的影响，遥感时间序列中往往存在大量缺失数据，而缺失数据对地表变化结果的准确性将产生较大影响[13-14]。更关键的是，大范围、长时间、高空间分辨率的遥感变化检测要应用大量遥感影像数据，需下载的数据量往往达TB级，庞大的数据量大大增加了数据处理和模型解算的成本，因此，节约成本、提高效率也是遥感时序地表变化检测的关键问题[15]。

为更好地改善以上问题，本文拟开展基于谷歌地球引擎(Google Earth Engine，GEE)云平台和BFAST算法的地表变化检测方法研究。首先利用GEE云平台，批量处理遥感时序数据，生成遥感时间序列数据集，以此提升数据处理和运算效率，降低因数据处理环节过多而引起的不确定性和误差；其次通过BFAST算法改善模型对遥感时间序列中确定阈值、对抗缺失数据等问题，增强模型对遥感时序数据中变化的分析能力，提升模型对遥感时序数据中的断点数据的识别精度；最后以河南某生态保护修复工程为实验区，开展实验验证。

1 研究区与数据源

1.1 研究区

在河南某生态保护修复工程实施范围内选择局部区域作为研究区。该生态保护修复工程涉及安阳市、鹤壁市、新乡市、焦作市和济源市5市，总面积约1.4×104km2，实施期限为2018—2020年。选择的研究区位于河南省鹤壁市的淇滨区和山城区，地理坐标为35°47′50″～35°52′28″N、114°04′10″～114°11′01″E，区域内海拔120～550 m，暖温带大陆性季风气候，年平均气温14.6℃、平均降雨量548.8～863.4 mm[16]。研究区域东西宽10.41 km，南北长8.61 km，总面积8 963 hm2，区域内海拔由西南向东北逐渐降低，地形地貌从山地到平原呈阶梯式变化，层次较为分明。土地覆被类型多样，适合从不同土地覆被类型角度开展地表变化检测实验，验证实验结果。研究区地理区位及遥感影像数据如图1所示。

图1 研究项目区地理位置和遥感影像(真彩色、NDVI)Fig.1 Location of study area and remote sensing image

1.2 实验数据

Landsat系列卫星传感器数据的波段范围、时空分辨率较为一致[17-21]，因此，可以很好地重构长时间序列数据[22]。为最小化数据处理过程中的误差，选择Landsat8/OLI遥感影像作为数据源，经范围和云量等条件设置，最终筛选了98景遥感影像参与实验，时间跨度为2013年3月至2020年12月。研究区遥感影像数据统计量如图2所示。

图2 研究区遥感影像数据量统计Fig.2 Data volume statistics of remote sensing images

1.3 GEE云计算平台

GEE是当前世界上较为先进的处理卫星影像等地理空间观测数据的云计算平台。GEE云端数据库中集成了近40年的Landsat系列卫星的历史存档数据，给个人用户提供了强大的算力和云存储空间，同时提供了方便快捷的JavaScript语言API接口进行数据处理、算法实现以及结果分析[23-26]。本研究应用GEE云计算平台进行数据处理，减少了数据准备的前期工作，降低了数据处理与算法实现过程中对本地硬件设备的依赖。

2 研究方法

首先基于GEE的遥感时序数据构建，包括Landsat8/OLI地表反射率(Landsat surface reflectance)数据集的调用、条件筛选、镶嵌、裁剪，基于CFMask算法的遥感数据集进行云影掩膜，在数据预处理基础上开展光谱指数(NDVI)计算和时间序列数据集构建等。其次基于时序数据集与BFAST算法进行地表变化检测[27]。最后，基于Google Earth开展精度评价。具体技术路线如图3所示。

图3 技术路线Fig.3 Workflow of this study

2.1 数据处理和时序构建

数据处理和时序构建过程包括：①结合实际研究需求，基于GEE平台确定并生成研究区的矢量边界数据。②采用ImageCollection相关函数调用GEE平台内嵌的Landsat SR数据，数据主要为LC08/C01/T1_SR，以此数据作为实验的数据源。③采用Filter相关函数对数据进行条件筛选、影像剪裁，将质量较差、不符合实验标准的数据进行剔除，并统一剪裁至研究区矢量形状。④针对数据集，构建掩膜函数，开展基于CFMask算法的云影掩膜[28]。⑤计算NDVI数据，对步骤①～④精选的LC08/C01/T1_SR数据统一批量进行NDVI计算，并采用Reducer.percentile函数对NDVI数据进行噪声去除，符合实验标准NDVI数据集。⑥构建并导出NDVI时间序列数据集合(98景NDVI数据影像)，并生成时间序列数据立方体。

2.2 BFAST时序分解模型构建

基于时序数据的地表变化检测，是通过时间序列分解模型来探究光谱变量(如NDVI)和时间变量之间的关系，时序分解模型通常为加性的线性回归模型[29-31]，计算式为

Yt=Tt+St+et(t=1,2,…,n)

(1)

其中Tt=ai+bit(τi-1≤t<τi,i=1,2,…,m)

(2)

(3)

式中Yt——在时间t范围的观测数据

Tt——趋势项St——季节项

ai、bi——趋势项系数

et——残差项γj——振幅

δj——分段数量f——频率

其中振幅、分段数量为未知项，而频率为已知项。

联合式(1)～(3)可将式(1)改写为一个由基函数组成的广义线性模型，即

(4)

其中

sin(2πkt/f),cos(2πkt/f))
β=(ai,bi,γ1cosδ1,γ1sinδ1,…,
γkcosδk,γksinδk)T

式(4)中的β是一个未知参数集，通过最小二乘法估计。式(4)即为时序分解模型。该模型通过对未知参数集的估计实现对已知数据的有效拟合，进而开展时序变化检测。

2.3 时序变化检测

通过MOSUM(Moving sums of the residuals)对时序的结构变化进行检测，计算式为

(5)

h——MOSUM的带宽

n——历史阶段观测样本数量

取消药品加成前后住院病人性别、年龄构成、医保身份构成，平均住院天数无差异。2016年 26 550人次 ，男性10 832人次，女性15 718人次；平均年龄为56.25岁；医保身份城镇职工20 433人次，城镇居民6 117人次；平均住院天数10.10天。2017年 27 304人次，男性10 796人次，女性16 508人次；平均年龄为56.45；医保身份城镇职工21 110人次，城镇居民 6 194人次；平均住院天数10.07天。

通常根据观测样本的数量选取，如h=n或h=n/4[32-33]，当时间序列处于稳定状态时，|MOt|接近于0且随机波动，而当时间序列出现较大变化时，|MOt|将系统性的偏离0，而当|MOt|偏离0超过95%的显著性边界时，即判定时间序列结构出现断点。

在时序变化检测过程中，需要在时序中设置并明确两个阶段：①历史阶段，此阶段基于观测数据进行模型训练。②检测阶段，此阶段基于训练好的模型对数据进行拟合，进而判别检测阶段数据的结构变化以及断点。

确定历史阶段的方法通常为：①基于先验知识。②基于Reversed-Ordered-Cumulative sum(ROC)。③基于Bai Perron(BP)断点估计方法[34]。而在本实验中，因生态保护修复工程实施周期为2018—2020年，因此基于先验知识，拟定检测阶段为2018—2020年，历史阶段为2013—2018年。

3 结果与分析

3.1 数据处理和时序构建

基于GEE云平台，实验得到符合要求的NDVI时间序列数据，数据共计98景，从统计描述看，98景NDVI影像数据的均值在0.037～0.690范围内波动，波动周期符合年度季相特征。标准差在0.03～0.16范围内波动，稳定性较好，数据质量较好，时序NDVI数据构建结果符合预期。研究区数据统计结果见图4。

图4 时序NDVI数据的均值和标准差Fig.4 Mean and variance of time series NDVI data

3.2 地表变化检测结果

地表变化检测结果主要包括：①检测的扰动(断点)发生时间，即扰动(断点)发生的年份和月份。②检测的扰动(断点)变化强度，即模型预测值和实际观测值之间的残差中值。对98景NDVI数据立方体执行算法后，结果如图5所示。

图5a为扰动变化强度，值域为-0.47～0.36，均值为-0.052，标准差为0.065，从统计结果并结合图5a视觉效果，可知研究区内2018—2020年的扰动变化以负向为主。

图5b～5d分别为2018—2020年3个年度的扰动时间(以月度计量)，其中2018年发生扰动的像元5 184个，扰动主要发生在8月之后，为4 414个；2019年发生扰动的像元14 009个，其中发生在8月之前的扰动为7 240个，8月之后的为6 769个；2020年发生扰动的像元3 435个，其中发生在8月之前的扰动为1 396个，8月之后的扰动为2 039个。从扰动时间统计结果并结合图5b～5d视觉效果，可知2018—2020年间，扰动发生的时间主要集中在下半年，其中2019年扰动发生的像元数量最大、面积最广。

图5 研究区扰动强度与发生时间Fig.5 Disturbance intensity and time in study area

3.3 结果验证与精度评价

基于上述方法能够有效、准确检测到地表变化，采用随机抽样方法分别从扰动变化强度、扰动发生时间的检测结果中抽取像元，组成参考数据集，同时基于Google Earth高分辨率影像数据，进行结果验证和精度评价。

首先，从检测结果中选取扰动较为集中连片的区域(2020年修建的道路)开展结果验证，如图6所示。图6a为发生在2020年的扰动识别结果。图6b和图6c为Google Earth高分辨率遥感影像，图6b为2019年影像，图6c为2020年影像，图6d为2020年影像标记点细节，图6a～6d中的标记点(红色十字和黄色图钉)为同一点，是识别结果中的第23 733号像元，图6e为标记点的时间序列数据，和扰动识别结果发生的时间节点。

图6 2020年扰动识别结果Fig.6 Disturbance identification results in 2020

对比图6a和图6c可看出，图6a的识别结果准确，虽然图6a的空间分辨率较低，但参照图6c高分辨率影像可看出，图6a对道路轮廓和形态的识别结果清晰；对比图7b和图7c可看出，该道路在2019年确未修建，而是在2020年修建完成，验证了图6a识别结果时间节点的准确性。

图7 2019年扰动识别结果Fig.7 Disturbance identification results in 2019

从图6e时间序列数据可看出，在2020年之前，预测变量对观测变量进行了良好的拟合，直到2020年第2季度，标记点所在的空间位置地表发生了连续剧烈的扰动(即道路的修建)，进而导致观测变量大幅偏移了预测变量。通过图6e可以进一步佐证图6a识别结果的准确性，也说明基于本文提出的方法可有效检测和识别地表变化过程。

图7a为扰动识别结果(2019年修建的地表建筑物)，图7b和图7c为Google Earth高分辨率遥感影像，图7b为2018年影像，图7c为2019年影像，图7a～7c中的标记点(红色十字和黄色图钉)为同一点，是识别结果中第28182号像元，图7d为标记点的时间序列数据，和扰动识别结果发生的时间节点。对比图7a～7c可得出，与图6结论一致，图7a对地表建筑物的轮廓和形态识别清晰，识别结果时间节点准确，而通过图7d进一步验证了识别结果的可靠性和准确性。

从图6和图7时序数据上看，属于规律变化后的突变，即在地表变化发生前，地表覆被基本呈规律性变化，而图8a～8e则验证了本文提出的方法在地表覆被不规律变化中对扰动的检测能力，其中图8a～8d为Google Earth高分辨率遥感影像，图8e为标记点(黄色图钉)的时间序列数据和扰动识别结果发生的时间节点。如图8a～8d所示，该区域地表覆被变化复杂，类型包括植被、裸土、水体以及三者的混合存在，2015年5月区域内以裸土为主略带植被，而2016年3月区域内则出现大面积水体，而2018年7月区域内水体大幅减少，受季节性因素影响区域内同时存在大量植被，而至2020年11月区域内水体基本转化为了裸土和植被。

图8 地表覆被不规律变化中的扰动检测能力Fig.8 Disturbance detecting capability in irregular changes of land surface cover

从图8e时间序列数据看，标记点2015年4季度NDVI存在明显的下降且由正转负，说明该阶段其空间位置的地表覆被已由植被和裸土转换为水体；2016—2018年，NDVI在-0.4～0.4之间波动，说明该阶段其空间位置的地表覆被受季节性因素影响，为植被、裸土、水体交替存在；从2018年开始，NDVI逐步在0以上形成波动，说明该阶段其空间位置的地表覆被中的水体逐渐消失。由此可见，时间序列数据与Google Earth高分辨率遥感影像呈现的趋势一致。虽然，从图8e的时间序列数据看，NDVI总体呈先下降而后升高趋势，但本文提出的方法对该趋势仍进行了有效拟合，并在检测阶段(2018—2020年)准确判定了扰动(断点)发生的年份和月份。进一步说明了，本文方法在地表NDVI复杂变化的条件下，对地表变化识别仍具有准确性。

3.3.2精度评价

通过图6～8从局部验证本文提出的方法对于检测地表变化的有效性和准确性，但还需通过样点验证进一步开展精度评价，扩大验证范围，提升结果验证和精度评价的稳健性。地表年度变化扰动的精度评价样点选取是较为困难的，尤其是历史已经发生过的变化扰动难以追溯。 因此，本文选择逆向精度评价，即在年度检测结果中抽样，将样点叠置在Google Earth高分辨率影像中开展精度评价验证。

分别从2018—2020年扰动检测结果中随机选取样点，2018年选取样点像元625个，2019年选取样点像元1 260个，2020年选取样点像元412个。其中，并非选取的所有样点都可用于精度评价，如选取的年度样点其空间位置所在的Google Earth中没有本年度和邻近年度高分辨率遥感影像数据，则无法开展精度评价，则该样点即为无效样点。通过对样点数据的筛查， 2018—2020年选取样点中的有效样点占比分别为85.1%、77.5%、89.1%，数量分别为532、976、367个像元，有效样点数量可满足精度评价条件。将有效样点叠置在Google Earth高分辨率影像中，通过目视解译，开展精度评价，2018—2020年的扰动检测结果正确样点率分别为86.5%、80.7%、87.7%，总体正确样点率为83.7%，正确率较高，如表1所示。

表1 精度评价Tab.1 Accuracy evaluation results

4 结束语

基于GEE和BFAST算法的地表变化检测方法，对长时间序列卫星遥感影像的变化信息进行了提取。实验结果分析表明，通过GEE能够高效获取并生成符合实验标准和要求的NDVI长时间序列数据，可以大幅度减少数据准备的前期工作，极大地降低了数据处理与算法实现过程中对本地硬件设备的依赖，节约了大量人力物力成本。通过BFAST算法对NDVI长时间序列数据进行变化检测，能够对时间序列的季相趋势进行充分拟合，分析趋势性变化，也能够准确检测到时间序列中的突发性变化，进而准确识别并检测到地表的变化扰动信息。针对本文实验区，地表变化检测结果的总体精度为83.7%，2018—2020年分年度检测结果精度分别为86.5%、80.7%、87.7%，精度较高。