一种基于多源遥感的森林年际物候监测方法

摘要

本发明公开一种基于多源遥感的森林年际物候监测方法，首先，收集所有可用的云量低于80％的卫星遥感影像，然后对不同卫星遥感影像的整合方法进行修正，提高不同传感器的空间和光谱匹配度；接着运用改进的连续变化检测与分类模型，生成每日植被指数曲线；最后，基于每日合成影像采用逻辑回归模型检验增强型植被指数、归一化植被指数和地表水体指数提取最优森林年际SOS。本发明改进不同卫星数据的整合方法，增加观测频率；提出MCCDC模型，将辐射差异纳入考虑范围并优化模型算法，在保证精度的同时缩短计算时间，最终生成每日清晰无云遥感影像；引入3种植被指数估测森林年际SOS，评估不同指数在评估森林SOS的差异。

说明书

技术领域

本发明属于森林物候监测技术领域，具体涉及一种基于多源遥感的森林年际物候监测方法。

背景技术

近几十年来，针对地面物候观测网难以大规模监测与制图的局限性，遥感逐渐成为监测森林物候动态的有效手段。其中，一些低空间分辨率传感器由于具备高时间分辨率而被广泛应用于物候信息提取，如高级甚高分辨率辐射计(AVHRR)和中分辨率成像光谱仪(MODIS)具有高时间分辨率，并能够获取500米至1100米空间分辨率的遥感数据，具有大范围观测优势，国内外学者在此基础上提出众多模型算法估算森林生长季开始日期(SOS)。然而，低空间分辨率传感器在面向区域尺度观测时，在地表异质性强和破碎度较高的区域，图像中会包含大量的混合像元，从而造成严重的估计误差；中等空间分辨率传感器如Landsat(美国NASA的陆地卫星)具有30米的空间分辨率和长时间跨度，能够克服上述的限制，然而，其16天的重访周期和常见的云干扰使Landsat难以准确估测森林年SOS。

近年，欧洲航天局分别发射Sentinel-2的两颗卫星(Sentinel-2A和Sentinel-2B)，其与Landsat光谱分辨率相似，这为提高植被物候观察的时间分辨率提供了可行的方案。Landsat和Sentinel-2共同组合提供2.9天的平均观测间隔，大大增加了在区域尺度识别森林SOS的能力，目前NASA科学家提供了一套Landsat和Sentinel-2的整合产品，但是有研究表明整合产品并不适用于所有地区，并且该整合后的产品在中国的绝大多数地区是不可共享的，因此需要重新修改整合方法生成新的整合产品。另外对于持续多云的地区，这些观测远不足以准确评估森林SOS。

目前大多数遥感物候研究是基于高时间频率和粗空间分辨率的传感器进行的，不适用于区域尺度的森林物候变量的精准刻画。中等空间分辨率的传感器如Landsat往往用于刻画长时间序列的森林平均SOS，提取年SOS成为一个大的挑战。增加观察频率或合成无噪图像成为中等空间分辨率遥感监测植被年森林SOS的有效方法。

Melaas等人为增加观测频率检测多个Landsat轨道数字重叠的区域的植被物候，该方法通常不适用于大规模物候估计。Bolton等人研究发现，与基于MODIS的监测评估相比，整合Landsat 8和Sentinel-2时间序列提高了景观尺度评估植被物候模式的精确度。

合成无噪图像解决了多云雨地区出现的云和阴影干扰问题，是提高森林年SOS估测准确度可行策略之一。此外有学者尝试利用图像融合技术，利用现有的Landsat影像，从MODIS中创建高时间分辨率影像，例如Bhandari等人利用增强时空自适应反射率融合模型(ESTARFM)算法合成未获得或不合适的Landsat影像，最终生成间隔8天的植被指数时间序列来评价森林物候。然而，该算法在图像融合过程中不可避免会产生不确定性，并且这些不确定性在高度异质的景观中尤为常见。

综上所述，对于森林年际物候监测存在以下问题：

(1)目前美国宇航局科学家研发的Landsat和Sentinel-2集成产品在中国并不完全提供，因此改进Landsat和Sentinel-2的整合方法成为目前物候监测待解决问题之一。

(2)CCDC模型进行精确估测的前提是需要大量且清晰的Landsat观测值，若出现持续阴雨多云天气，则可能出现模型过拟合或欠拟合问题，不具有普适性，并且CCDC模型并未考虑不同传感器之间的辐射差异。

(3)CCDC模型利用前12个清晰观测值建立初始时间序列来模拟SR运动轨迹，跟踪SR变化周期，并对比模型预测值与实际观测值，通过设定阈值，判断土地覆盖是否发生变化，迭代更新模型参数，预测任意期望日期的Landsat地表反射率(SR)。因此，CCDC模型对每个观测值逐像元进行判断，并迭代更新模型参数，计算量大，工作效率低，针对Landsat场景中地表覆盖类型不变的地块，逐像元比较实际观测值和模型预测值十分耗时。

发明内容

本发明解决的技术问题：本发明提供一种具有普适性、计算效率高、预测结果精准的基于多源遥感的森林年际物候监测方法。

技术方案：为了解决上述技术问题，本发明采用的技术方案如下：

一种基于多源遥感的森林年际物候监测方法，首先，收集所有可用的云量低于80％的卫星遥感影像，然后对不同卫星遥感影像的整合方法进行修正，提高不同传感器的空间和光谱匹配度；接着运用改进的连续变化检测与分类模型，生成每日植被指数曲线；最后，基于每日合成影像采用逻辑回归模型检验增强型植被指数、归一化植被指数和地表水体指数提取最优森林年际SOS。具体步骤包括：

S1:获取数据，并对数据进行预处理；

S2:建立连续变化检测与分类模型，利用连续变化检测与分类模型生成每日植被指数曲线；

S3：将每个像元的建立连续变化检测与分类模型拟合系数输入支持向量机分类器进行森林覆盖制图；

S4：进行年际森林生长季开始日期信息提取。

进一步地，步骤S1中，数据包括Landsat ETM+/OLI的陆地卫星影像、Sentinel-2影像、森林资源二类调查数据和物候实测数据；

进一步地，步骤S1中，数据的预处理包括大气校正、云和阴影掩膜、图像配准和图像光谱归一化处理。

进一步地，步骤S2中，连续变化检测与分类模型考虑像元的不同变化类型，在获取足够清晰像元基础上拟合模型，给定判断标准来确定观测期内是否有土地覆盖变化，如果有，则逐一判断每个观测值，并识别土地覆盖变化发生的具体日期，最后通过给定的模型生成每日遥感影像。

进一步地，步骤S3中，选用径向基函数作为核函数，定义成本参数C，参数γ设定为输入特征个数的倒数，得到若干系数作为支持向量机分类的输入特征。

进一步地，步骤S4中，针对云和阴影区域生成逐日地表反射率数据，并与清晰的Landsat和Sentinel-2观测值融合，计算归一化植被指数、增强植被指数和地表水体，采用逻辑斯蒂回归模型分析3种植被指数在整个观测期间的逐像元波动模式，将拟合模型上升阶段的变化率的最大值对应的日期作为森林年际物候。

有益效果：与现有技术相比，本发明具有以下优点：

(1)本发明通过采用线性模型(斜率和斜距)来归一化和最小化多源数据的辐射差异，将Sentinel-2数据与Landsat系列影像数据整合，生成具有较高时间分辨率的Landsat影像，周期约为2.9日，比目前常用的ESTARFM算法合成Landsat图像更加准确可信。

(2)改良了CCDC模型，在CCDC模型基础上，引入多源数据增加观测频率，并提出一个合适的判断标准来确定观察期内土地覆盖是否发生变化。如果有，再逐一判断每个观测，并定位土地覆盖变化发生日期，在确保精度的同时缩短计算时间，提高计算效率。

(3)通过Landsat和Sentinel-2的整合以及合成清晰无云影像消除了大部分的云、阴影噪声，使得物候监测不再受困于多云雨天气地区森林物候信息提取，在大部分地区更具普适性。此外本方法能将区域尺度森林SOS制图产品扩展至30米空间尺度，并有扩展至10米观测尺度的直接潜力。

(4)采用归一化植被指数(NDVI)、增强型植被指数(EVI)和地表水体指数(LSWI)三种植被指数提取森林年际SOS，通过择优比对获得更加精确的物候信息。

附图说明

图1是基于多源遥感的森林年际物候监测方法流程图；

图2是MCCDC模型示意图；

图3是合成图像(右)与真实图像(左)的目视对比；

图4中(a)是2013年至2019年森林覆盖制图；(b)2013年至2019年森林面积变化；

图5是三种植被指数估测SOS的验证结果；

图6中(a)是基于3种植被指数提取的SOS空间分布，(b)与(c)2013-2019年2个不同地点年际SOS动态变化。

具体实施方式

下面结合具体实施例，进一步阐明本发明，实施例在以本发明技术方案为前提下进行实施，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围。

本实施例公开一种基于多源遥感的森林年际物候监测方法，方法流程图如图1所示：首先，收集所有可用的云量低于80％的Landsat和Sentinel-2影像，然后对Landsat和Sentinel-2的整合方法进行修正，提高不同传感器的空间和光谱匹配度。接着运用改进的连续变化检测与分类(MCCDC)模型，生成30m空间分辨率的每日植被指数曲线。最后，基于每日合成影像采用逻辑回归模型检验增强型植被指数(EVI)、归一化植被指数(NDVI)和地表水体指数(LSWI)提取最优森林年际SOS。具体包括以下步骤：

S1:获取卫星遥感数据，并对数据进行预处理；

首先，获取数据，本发明的数据包括Landsat ETM+/OLI的陆地卫星影像、Sentinel-2影像、森林资源二类调查数据、物候实测数据。

(1)从美国地质调查局(USGS)网站下载2013年1月1日至2019年12月30日云量小于80％的Landsat ETM+/OLI的陆地卫星影像，下载的数据轨道号为131/42，共计188景影像。ETM+影像为89景，OLI影像为99景。所选Landsat图像的基本信息如表1所示，其中LandsatETM+/OLI均为地表反射率产品；

(2)从欧洲航空局(ESA)网站下载2017-2019年3-11月云量小于80％的Sentinel-2大气表观反射率(TOA)辐射率数据，所用光谱波段对应于Landsat(轨道号T47RPJ，Sentinel-2A58景和Sentinel-2B38景)，基本信息如表1所示；

(3)由林业部门提供2016年的研究区森林资源二类调查数据，用于模型精度验证；

(4)40个物候观测样地，其中28个观测数据来源于场地工作人员记录，12个观测数据由作者实地调查获取，用于模型精度验证。

表1获取的遥感影像基本信息

然后，对获取的数据进行预处理，数据预处理包括以下步骤：

(1)大气校正：Landsat的SR数据直接由USGS提供，Landsat ETM+使用LEDAPS算法，OLI影像使用LaSRC算法进行大气校正。Sentinel-2数据由ESA提供，TOA使用Sen2cor辐射工具校正为SR产品。

(2)云和阴影掩膜：云和雪在可见光波段具有高亮度特征，而云影则为低亮度，容易与其他土地覆盖类型混淆，造成干扰。因此，本发明使用Fmask 4.0软件对Landsat和Sentinel-2影像进行云及对应阴影干扰源去除。

(3)图像配准：Landsat和Sentinel-2数据空间位置不匹配，需要进行图像配准。以Sentinel-2为参考影像，在Landsat与Sentinel-2重叠区域手工选取多个地面控制点，拟合一阶多项式进行批量配准，配准的几何误差小于0.5像素。

(4)Landsat和Sentinel-2的图像光谱归一化处理：Landsat和Sentinel-2传感器在特定波长下具有不同的光谱响应，因此需要进行逐波段调整以消除传感器差异引起的辐射差异。

本发明采用通用系数和c因子法进行BRDF校正，然后在采集日期相同或相邻的影像对上识别多个像元样本(确保影像对没有发生土地覆盖变化)，然后拟合线性模型(斜率和截距)进行逐波段辐射归一化操作，并采用最邻近点重采样方法将Sentinel-2影像重采样到30m。

本发明通过修改Sentinel-2数据与Landsat数据的整合方法，采用线性模型(斜率和斜距)来归一化和最小化多源数据的辐射差异，以生成更一致的合成图像。

S2:建立连续变化检测与分类模型，利用连续变化检测与分类模型生成每日植被指数曲线；

MCCDC模型考虑像元的不同变化类型，在获取足够清晰像元基础上拟合模型(式1)，给定判断标准来确定观测期内是否有土地覆盖变化。如果有，则逐一判断每个观测值，并识别土地覆盖变化发生的具体日期，最后通过给定的模型生成每日遥感影像。

首先将每个波段同一位置的观察期内所有清晰的像元按照图像采集时间排序，基于式(1)拟合出完整时间序列模型，然后计算模型R

根据式(2)确定一个土地覆盖变化事件后，需要准确提取其在某一具体变化日期的发生情况。对于发生土地覆盖变化的像元，本发明基于CCDC模型的基本原理，将前24个相同位置的清晰像素拟合一个初始化的时间序列模型，然后判断该像元后续的4个连续观测值与预测值的差异是否超出模型均方根误差(RMSE)的3倍，如果超出范围，则将第1次出现的时间定义为土地覆盖变化的发生时间。

式中：x是儒略日期，T是一年中的天数，a

图2给出了不同位置2个像元的NIR和SWIR1波段拟合曲线，可以看见MCCDC在考虑土地覆盖变化发生时能很好地拟合SR的规律性变化，其中(a)和(b)中的拟合曲线表明第一个像元的土地覆盖类型未发生变化，(c)和(d)中的拟合曲线表明第二个像元在观测期间土地覆盖发生变化。

S3：进行森林覆盖制图

本发明将每个像元的MCCDC模型拟合系数输入支持向量机(SVM)分类器进行森林覆盖制图，选用径向基函数(RBF)作为核函数，成本参数C定义为10，参数γ通常设定为输入特征个数的倒数，共48个系数(每个波段8个系数，共6个波段)作为SVM分类的输入特征。

S4:年际森林生长季开始日期(SOS)信息提取

针对云和阴影区域生成逐日地表反射率数据，并与清晰的Landsat和Sentinel-2观测值融合，计算归一化植被指数(NDVI)、增强植被指数(EVI)和地表水体(LSWI)，采用逻辑斯蒂回归模型分析3种植被指数在整个观测期间的逐像元波动模式，将拟合模型上升阶段的变化率(ROC)的最大值对应的日期作为森林SOS。

其中ρ

S5:对精度评价

通过计算R

该方法有以下优势：(1)改进Landsat和Sentinel-2的整合方法，增加观测频率；(2)提出MCCDC模型，将Landsat 7ETM+与Landsat 8OLI的辐射差异纳入考虑范围并优化模型算法，在保证精度的同时缩短计算时间，最终生成每日清晰无云遥感影像；(3)引入3种植被指数估测森林年际SOS，评估不同指数在评估森林SOS的差异。

对本实施例的森林年际物候模型进行结果分析

(1)多源数据的整合和合成Landsat图像评价

表2展示了Landsat 8与Landsat 7和Sentinel-2的SR匹配的斜率、截距和R

表2 Landsat和Sentinel-2的逐波段光谱归一化变换系数

表3显示了Landsat 7/8的合成波段与其对应的真实波段的R

表3合成图像与真实图像各波段R

通过目视比较2013年冬季和2019年春季的原始图像与合成图像(图3)，对于不同的土地覆盖类型，合成图像与原始图像十分相似，更重要的是，合成图像完全克服了噪声的影响，如云、阴影和Landsat 7的SLC-off图像缝隙。对于一些特征地物，如湖泊、道路、森林、建筑物均能够较好地刻画出来。

(2)森林覆盖动态变化

通过对比2016年森林资源二类调查数据与2016年森林覆盖分类图进行对比，结果显示森林的用户精度为89％，生产者精度为88％，总体精度达到88％。在此基础上，利用Landsat影像数据对2019年森林覆盖分类图进行验证，结果表明：森林与非森林的用户精度、生产者精度和总体精度均大于88％。总体而言，这些精度表明现有方法获得的森林覆盖图可信(表4)。森林面积从2013年的4015.69平方公里下降到2019年的4010.17平方公里(图4)。

表4森林覆盖制图精度评价

(3)森林SOS提取

本发明从2013年、2018年和2019年的40个样本中采集森林SOS信息，并将基于遥感提取的SOS与野外观测数据进行对比，获取不同植被指数估测的森林SOS的R

从图5上看，EVI预测的SOS与实测SOS的R

图6(a)显示了2013年和2019年从三个指数中提取的SOS分布。图6(b)和(c)显示了2013年至2019年两个不同站点的年度SOS动态。NDVI和LSWI提取的SOS的时空分布相似，大多数SOS主要来自80-150天。具体而言，从NDVI和LSWI提取的研究区西部SOS以80-120天为主，底部和东南角的SOS以120-150天为主。北部SOS晚于150天(图6(a))。总体而言，南部SOS早于北部，西部SOS早于东部SOS。对于一些研究区域，三个指数的最大SOS差异不超过一周(图6(b))。然而，在研究区的其他小部分，不同指标的SOS差异极大。当基于EVI的SOS小于40天时，来自NDVI和LSWI的SOS晚于90天，差异达到两个月(图6(c))。

本发明改进后的模型比现有模型要复杂很多(参数更多)，既考虑了趋势项，又考虑了年际变动；本发明引入Sentinel-2数据，可以明显增加观测点的密度，弥补了CCDC模型在对于常年多云和多雨的地区可能存在精度偏低的现象，使该模型更加具有普遍性；本发明提出的MCCDC模型在模型精度和计算效率之间做了个取舍，在保证一定精度的同时，减少了计算时间。传统的CCDC模型计算可能需要几天的时间，但是MCCDC模型只需要十几个小时即可完成；CCDC模型的建立未考虑Landsat5，7和8的传感器差异。事实上不同地物的反射率会有所差异。本文通过线性回归消除它们之间的反射率差异，生成更为准确的合成图像。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。