基于时间序列卫星影像的桉树人工林月尺度林龄反演方法

摘要

本发明提供了基于时间序列卫星影像的桉树人工林月尺度林龄反演方法，涉及森林资源调查与林业定量遥感技术领域。该方法能够更加快速准确地获取桉树人工林林龄数据，弥补了传统林龄获取方法耗时耗力、成本高昂、数据获取周期长的不足；采用长时间序列遥感影像植被指数在桉树砍伐前后的变化进行林龄反演能够弥补基于林龄和遥感光谱反射之间的关系以建立林龄反演模型反演林龄的不足，实现时间序列林龄的精确提取；通过使用时空融合算法将高空间分辨率卫星遥感影像与低空间分辨率卫星遥感影像进行时空融合，重构月份因云雨影响而不能使用所缺失的影像数据，构建完整的时间序列影像，进而将桉树人工林林龄反演精确到月尺度。

说明书

技术领域

本发明属于森林资源调查与林业定量遥感技术领域，具体涉及一种基于时间序列卫星影像的桉树人工林月尺度林龄反演方法。

背景技术

桉树不仅生长快速，具有较好的木质特性，是纸浆和造纸生产中最重要的短纤维来源，而且在生态安全、调节全球气候变化及保障国家木材安全等方面也发挥着重要作用。桉树作为三大速生树种之一，成材周期短，通常以木材生产为主，因此，现阶段桉树人工林经营多以短轮伐期为主。

林龄作为森林的基本结构参数之一，对森林可持续经营、森林资源定量研究等方面有着不可或缺的作用。在基于时间序列光学影像的桉树人工林林分蓄积量的动态变化监测中，由于光学影像在高生物量或高郁闭度情况下容易出现信号饱和，导致通过基于时间序列光学影像的桉树人工林林分蓄积量的动态变化监测结果不准确，而将林龄信息作为辅助数据引入桉树人工林林分蓄积量反演模型则能够提高其动态变化监测结果的精度。林龄信息通常多来自森林资源清查，迫切需要一种先进的能够快速便捷准确地获取桉树人工林时间序列林龄数据的方法。

随着遥感技术的发展，一些学者开始利用卫星遥感影像反演林龄，如研究学者基于卫星遥感影像提取光谱及植被指数，建立林龄反演模型，如Jensen等基于TM影像，利用人工神经网络构建火炬松林龄反演模型；唐少飞等基于Sentinel-2影像，利用随机森林、支持向量机、神经网络等多种算法构建落叶松林龄估测模型。虽然先前研究结果表明基于卫星遥感影像能够实现林龄估测研究，但是，通过利用遥感光谱反射和林龄之间的关系建立林龄反演模型的方法限制因素太多，如林龄相似的近成熟林和成熟林光谱反射差异较小，难以精确区分；同时，林龄反演精度受物候条件影像较大，如展叶期与落叶期的林龄反演模型的精度差距很大，上述多种因素导致现有基于光谱反射与林龄相关关系的反演方法难以实现桉树人工林时间序列林龄的精确反演。因此，如何基于时间序列卫星遥感影像准确提取并及时更新桉树人工林林龄信息同时满足桉树人工林短轮伐期经营林龄数据的需求是短轮伐期桉树人工林资源经营管理拟解决的关键问题。

综上所述，桉树人工林林龄反演的现有研究还存在以下问题：

1、现阶段的桉树人工林林龄数据大多来自传统调查，而传统调查通常通过砍伐或者在树干上打孔，根据年轮测定树龄，这种方式不仅严重影响树木的生长而且耗时耗力成本高昂。此外，森林资源清查通常每五年进行一次，导致桉树人工林林龄数据更新频率较长，难以满足短轮伐期经营桉树人工林林龄数据的需求。

2、光学卫星影像因受云雨影响较大，导致难以获取高质量连续的时间序列影像数据，这在一定程度上限制了基于遥感影像的林龄反演结果多集中于年尺度，而难以进一步将林龄反演结果精确到月尺度，进而阻碍区域尺度桉树人工林蓄积量时间差异校正，导致区域蓄积量反演结果不准，给桉树人工林资源精准管理带来极大的困难。

3、利用遥感光谱反射和林龄之间关系建立林龄反演模型的方法限制因素太多，如林龄相似的近成熟林和成熟林光谱反射差异较小，难以精确区分；同时，林龄反演精度受物候条件影像较大，如展叶期与落叶期的林龄反演模型的精度差距很大，上述多种因素导致现有基于光谱反射与林龄相关关系的反演方法难以实现桉树人工林时间序列林龄的精确反演。

发明内容

针对现有技术存在的问题，本发明提供了一种基于时间序列卫星影像的桉树人工林月尺度林龄反演方法。旨在将桉树人工林林龄的反演精确到月尺度，提供更为准确的桉树人工林林龄数据，以支撑基于短轮伐期的桉树人工林经营需求。

本发明的技术方案是：

基于时间序列卫星影像的桉树人工林月尺度林龄反演方法，其特征在于，该方法包括以下步骤：

步骤1：确定桉树人工林目标研究区，利用Google Earth Pro软件通过目视解译的方法在N-m年至N年各年的谷歌历史卫星影像中：1)选择多个分布在研究区的桉树所在像元即桉树样本点，并记录其坐标；2)选择多个分布在研究区的已被砍伐桉树所对应像元并记录其坐标和砍伐更新时间，并将桉树的砍伐更新时间定义为砍伐更新点；其中m为大于等于研究区桉树轮伐期的正整数；

步骤2：选择与N-m至N年各年分别对应的经过正射校正和经过地理配准的覆盖研究区的Landsat影像数据，在谷歌地球引擎GEE平台上调用这些影像数据并对这些影像数据进行影像的拼接、掩膜提取和辐射归一化预处理，获得N-m至N年各年研究区影像数据；

步骤3：对N-m至N年各年研究区影像数据进行参数提取，所述参数包括光谱反射率、包含归一化燃烧指数NBR的植被指数以及纹理信息；根据桉树样本点的坐标，在N-m至N年各年研究区影像数据中找到对应的像元，由这些像元及其对应的参数构建数据集，从该数据集划分出训练集，利用随机森林算法进行桉树的分类训练，构建研究区各年桉树人工林空间分布提取模型；将从各年研究区影像数据提取的参数分别对应输入研究区各年桉树人工林空间分布提取模型，获得N-m至N年研究区各年桉树人工林空间分布；

步骤4：分别利用N-m至N年研究区各年桉树人工林空间分布对对应的分别从各年研究区影像数据获得的归一化燃烧指数NBR空间分布进行掩膜提取，分别对应得到各年的桉树人工林NBR空间分布，从而获得N-m至N年研究区桉树人工林NBR空间分布；根据各砍伐更新点在桉树人工林NBR空间分布中对应像元的NBR值的变化规律确定NBR骤降幅度阈值；

步骤5：根据N-m至N年研究区桉树人工林NBR空间分布、NBR骤降幅度阈值，确定N年研究区桉树人工林空间分布中的真实桉树砍伐更新点及对应的砍伐更新年；

步骤6：基于N年研究区桉树人工林空间分布中各像元对应的砍伐更新年，计算出N年研究区桉树人工林年林龄,进而获得N年研究区桉树人工林年林龄空间分布；

步骤7：遍历N年研究区桉树人工林空间分布中各个像元对应的距离N年最近的砍伐更新年内所有月份的Landsat影像数据，遍历完成转至步骤11；在遍历过程中，若当前遍历月份的Landsat影像数据不能完全覆盖研究区，存在缺失影像数据，则执行步骤8；若当前遍历月份的Landsat影像数据能完全覆盖研究区，不存在缺失影像数据则转至步骤9；

步骤8：引入当前遍历月份及其前后月份的MODIS研究区影像数据，通过STARFM时空融合算法获得与当前遍历月份对应的Landsat影像数据空间分辨率相同的所述缺失影像数据；

步骤9：在GEE平台上针对当前遍历月份的研究区完整影像数据进行影像的拼接、掩膜提取和辐射归一化预处理，得到预处理后的当前遍历月份影像数据；

步骤10：在GEE平台上直接提取当前遍历月份研究区影像数据的NBR，利用各像元对应的距离N年最近的砍伐更新年研究区桉树人工林空间分布对当前遍历月份的研究区NBR空间分布进行掩膜提取，得到当前遍历月份的研究区桉树人工林NBR空间分布，返回步骤7；

步骤11：根据N年研究区桉树人工林空间分布中各像元在其对应的距离N年最近的砍伐更新年内各月的研究区桉树人工林归一化燃烧指数NBR随时间的变化规律，定义在A年B月刚被砍伐的桉树人工林为0年又0月岁，结合N年研究区桉树人工林年林龄空间分布，可确定在N年研究区桉树人工林空间分布中距离N年最近的砍伐更新年A年内的真实砍伐更新点所对应像元中的被砍伐桉树在N年C月的林龄为Q年L月，其中Q为N与A之间的差值，L为C与B之间的差值。

进一步地，根据所述的基于时间序列卫星影像的桉树人工林月尺度林龄反演方法，其特征在于，所述NBR骤降幅度阈值为0.4。

进一步地，根据所述的基于时间序列卫星影像的桉树人工林月尺度林龄反演方法，其特征在于，所述步骤5中所述的确定N年研究区桉树人工林空间分布中的真实桉树砍伐更新点及对应的砍伐更新年的方法为：将N-m至N年各年研究区桉树人工林空间分布中所在像元发生NBR骤降对应的时间点定义为桉树伪砍伐更新点；然后以轮伐期为时间步长，在时间步长内：若像元的NBR值只骤降1次，且骤降的幅度大于等于NBR骤降幅度阈值，则该次NBR值骤降对应的伪砍伐更新点则为真实砍伐更新点；若像元的NBR值骤降大于1次，且骤降的幅度均大于等于NBR骤降幅度阈值，则取最后1次NBR值骤降对应的伪砍伐更新点为真实砍伐更新点；各真实砍伐更新点分别对应的当前年定义为各真实砍伐更新点的砍伐更新年。

进一步地，根据所述的基于时间序列卫星影像的桉树人工林月尺度林龄反演方法，其特征在于，基于N年研究区桉树人工林空间分布中各像元对应的砍伐更新年，计算N年研究区桉树人工林年林龄的方法为：定义各像元内的桉树在砍伐更新年为0岁，利用N年分别减去各像元对应的距离N年最近的砍伐更新年，则得到各像元内桉树在N年的年林龄；若N年研究区桉树人工林中部分像元在N-m至N年均无砍伐更新年，则视该部分像元内桉树在N年的年林龄为大于m岁。

与现有技术相比，本发明提出的技术方案具有以下有益效果：

(1)提供了一种能够更加快速准确地获取桉树人工林林龄数据的更为先进的方法，弥补了传统林龄获取方法耗时耗力、成本高昂、数据获取周期长的不足。

(2)通过使用时空融合算法将高空间分辨率卫星遥感影像与低空间分辨率卫星遥感影像进行时空融合，重构因云雨影响而不能使用所缺失月份的影像数据，从而构建完整的时间序列影像，进而将桉树人工林林龄反演结果进一步精确到月尺度，为区域桉树人工林蓄积量时间差异校正提供基础数据。

(3)采用长时间序列遥感影像植被指数在桉树砍伐前后的变化进行林龄反演能够弥补基于林龄和遥感光谱反射之间的关系以建立林龄反演模型反演林龄的不足，实现时间序列林龄的精确提取。

附图说明

为了更清楚地说明本发明实施例中的具体方式，下面将对实施例中涉及的相关附图做简单说明，下面的附图仅仅是本发明的优选实施例，对于本领域普通技术人员来说，在没有创造性改变的前提下，可以根据这些附图获得其它的附图。

图1为本实施方式基于时间序列卫星影像的桉树人工林月尺度林龄反演方法的流程示意图；

图2为本实施方式2006-2020年研究区各年桉树人工林空间分布图，其中(a)为2006年研究区桉树人工林空间分布图；(b)为2007年研究区桉树人工林空间分布图；(c)为2008年研究区桉树人工林空间分布图；(d)为2009年研究区桉树人工林空间分布图；(e)为2010年研究区桉树人工林空间分布图；(f)为2011年研究区桉树人工林空间分布图；(g)为2012年研究区桉树人工林空间分布图；(h)为2013年研究区桉树人工林空间分布图；(i)为2014年研究区桉树人工林空间分布图；(j)为2015年研究区桉树人工林空间分布图；(k)为2016年研究区桉树人工林空间分布图；(l)为2017年研究区桉树人工林空间分布图；(m)为2018年研究区桉树人工林空间分布图；(n)为2019年研究区桉树人工林空间分布图；(o)为2020年研究区桉树人工林空间分布图；

图3(a)为本实施方式某砍伐更新点对应的由黑色正方形所代表的像元示意图；(b)为(a)图中黑色正方形代表的像元对应的砍伐更新点的砍伐更新年及前后年的归一化燃烧指数NBR的变化曲线图；

图4为本实施方式2020年研究区桉树人工林年林龄空间分布图；

图5为本实施方式2020年研究区桉树人工林月林龄空间分布图。

具体实施方式

为了便于理解本申请，下面将参照相关附图对本申请进行更全面的描述。附图中给出了本申请的较佳实施方式。但是，本申请可以以许多不同的形式来实现，并不限于本文所描述的实施方式。相反地，提供这些实施方式的目的是使对本申请的公开内容理解的更加透彻全面。

本实施方式以基于2006-2020年的Landsat系列影像数据反演2020年广西壮族自治区桉树人工林月尺度林龄为例对本发明的基于时间序列卫星影像的桉树人工林月尺度林龄反演方法进行详细描述。图1是本实施方式基于时间序列卫星影像的桉树人工林月尺度林龄反演方法的流程示意图，如图1所示，该方法包括以下步骤：

步骤1：确定桉树人工林目标研究区，利用Google Earth Pro软件通过目视解译的方法在N-m年至N年各年的谷歌历史卫星影像中：1)选择多个分布在研究区的桉树所在像元即桉树样本点，并记录其坐标；2)选择多个分布在研究区的已被砍伐桉树所对应像元并记录其坐标和砍伐更新时间，并将桉树的砍伐更新时间定义为砍伐更新点；其中m为大于等于研究区桉树轮伐期的正整数；

本实施方式确定的桉树人工林目标研究区是广西壮族自治区，桉树的轮伐期通常为3-5年。本实施方式利用Google Earth Pro软件通过目视解译的方法在2006-2020年(年跨度比桉树的轮伐期足够长，才能够反演出林龄)各年的谷歌历史卫星影像中选择1000个分布在研究区的桉树样本点并记录其坐标，通过目视解译的方法在2006-2020年的谷歌历史卫星影像中任意选择600个桉树砍伐更新点并记录对应的像元坐标和砍伐更时间。

在该步骤的实施过程中，在Google Earth Pro软件中打开2006-2020年各年谷歌历史卫星影像，在谷歌历史卫星影像中找到广西壮族自治区的位置，由于桉树多为人工林，树与树之间有明显间隔，从谷歌历史卫星影像看，桉树人工林纹理特征明显，通过目视解译的方法在每年的谷歌历史卫星影像上选择1000个均匀分布在广西壮族自治区内的桉树样本点，共15000个桉树样本点，并记录各桉树样本点的坐标；通过目视解译的方法在2006-2020年谷歌历史卫星影像中选择分布在广西壮族自治区内的已被砍伐桉树对应的600个像元，其中，先在2006-2020年谷歌历史卫星影像中找到桉树集中分布的区域，通过目视解译的方法找到600个已被砍伐桉树对应的像元，2006-2020年每年找一部分，共600个像元并记录各像元的坐标，一个像元对应一个桉树砍伐更新点，对比已选择的600个像元前后几年的历史卫星影像，确定各像元内桉树的砍伐更新时间并记录。定义各像元内桉树的砍伐更新时间为各像元对应的桉树砍伐更新点。

步骤2：选择与N-m至N年各年分别对应的经过正射校正和经过地理配准的覆盖研究区的Landsat影像数据，在谷歌地球引擎GEE平台上调用这些影像数据并对这些影像数据进行影像的拼接、掩膜提取和辐射归一化预处理，获得N-m至N年各年研究区影像数据。

本实施方式选取与2006-2020年各年分别对应的经过正射校正和经过地理配准的覆盖研究区的最新Landsat影像数据，在谷歌地球引擎GEE(Google Earth Engine)平台上调用这些影像数据，并对这些影像数据进行预处理，包括影像的拼接、掩膜提取和辐射归一化，获得预处理后的2006-2020年各年的研究区影像数据。

在本实施方式的具体实施过程中，综合考虑Landsat系列卫星的发射时间，选取2006-2011年的Landsat-5影像数据，2012年的Landsat-7影像数据，2013-2020年的Landsat-8影像数据。在GEE平台上直接调用些影像数据。某年的研究区影像数据是使用某年的来自不同月份的多幅Landsat影像进行拼接得到的。例如使用能完全覆盖研究区的2006年不同月份的Landsat-5影像数据，通过将2006年不同月份的Landsat-5影像数据进行拼接，研究区矢量边界图可从网上下载，再利用研究区矢量边界图进行掩膜提取，将拼接后的影像中研究区内的部分提取出来，研究区外的部分则舍去，获得2006年研究区影像数据。以此类推，再获取2007-2020年各年的研究区影像数据。由于Landsat-5携带的TM传感器、Landsat-7携带的ETM+传感器和Landsat-8携带的OLI传感器所采集的影像之间存在一定的地表反射差异，因此需要进行辐射归一化处理，使用GEE自带的辐射归一化算法进行处理，2006-2020年研究区影像数据直接输入辐射归一化算法中即可得到辐射归一化后的2006-2020年各年的研究区影像数据。

步骤3：对各年研究区影像数据进行参数提取，所述参数包括光谱反射率、包含归一化燃烧指数NBR的植被指数以及纹理信息；根据桉树样本点的坐标，在N-m至N年各年研究区影像数据中找到对应的像元，由这些像元及其对应的参数构建数据集，从该数据集划分出训练集，利用随机森林算法进行桉树的分类训练，构建研究区各年桉树人工林空间分布提取模型；将从各年研究区影像数据提取的参数分别对应输入研究区各年桉树人工林空间分布提取模型，获得N-m至N年研究区各年桉树人工林空间分布；

本实施方式对预处理后的2006-2020年研究区影像数据进行参数提取，包括光谱反射率、植被指数以及纹理信息。利用步骤1所采集的桉树样本数据，通过ArcGIS软件构建数据集，采用随机森林算法进行分类，构建桉树人工林空间分布提取模型，基于预处理后的2006-2020年研究区影像数据，获得2006-2020年各年研究区桉树人工林空间分布。

在该步骤的实施过程中，在GEE平台上可以直接提取预处理后的2006-2020年研究区影像数据的光谱反射率、植被指数和纹理信息。其中，光谱反射率是卫星搭载的传感器不同的波段被不同地物反射的光通量与入射到物体上的光通量之比；植被指数包括：增强植被指数EVI、归一化植被指数NDVI、归一化植被指数NDVIA、陆地叶绿素指数MTCI、新型倒红边叶绿素指数IRECI、植被衰减指数PSRI、转化叶绿素吸收反射率指数TCARI、归一化差值水体指数NDWI、修正的叶绿素吸收反射率指数MCARI、重归一化植被指数RDVI、三角植被指数TVI、土壤调节植被指数SAVI、水分胁迫指数MSI、地表水分指数LSWI、归一化燃烧指数NBR、增强归一化红边植被指数mNdvired_edge、改进红边比值植被指数MSRred_edge和叶绿素红边指数CIred_edge；纹理信息包括：均差、方差、同质性、对比度、相异性、信息熵、二阶矩和相关性。

下面以2006年研究区桉树人工林空间分布的获取过程为例说明各年研究区桉树人工林空间分布的获取过程：将预处理后的2006年研究区影像数据在ArcGIS中打开，将由步骤1获得的2006年1000个桉树样本点的坐标输入ArcGIS，则能在预处理后的2006年研究区影像数据中找到对应的点。由于预处理后的2006-2020年研究区影像数据均为栅格图像，因此，是将1000个桉树样本点对应到2006年研究区影像数据中。利用ArcGIS中的栅格图像多值提取到点功能，将从2006年研究区影像数据所提取的对应像元的光谱反射率、植被指数和纹理信息提取到1000个桉树样本点中，则1000个桉树样本点各点的属性表都有其对应的光谱反射率、植被指数和纹理信息。以1000个桉树样本点的分类结果为桉树，1000个桉树样本点各点的属性表对应的光谱反射率、植被指数和纹理信息为特征，构建2006年研究区桉树人工林空间分布提取模型的数据集，取数据集中75％的数据作为训练集，剩余25％为验证集，利用随机森林算法进行桉树的分类训练，即可构建2006年研究区桉树人工林空间分布提取模型。再将2006年研究区影像数据输入上述构建的2006年研究区桉树人工林空间分布提取模型中，随机森林算法根据上述提取出的2006年研究区影像数据特征(光谱反射率、植被指数和纹理信息)，即可获得2006年研究区桉树人工林空间分布。2006-2020年各年研究区桉树人工林空间分布均是栅格图像，桉树以像元的形式组成桉树人工林空间分布。

以此类推，对2007-2020年研究区各年影像数据进行上述操作，即可获得2007-2020年研究区各年桉树人工林空间分布提取模型，进而得到2006-2020年研究区各年桉树人工林空间分布，如图2所示为2006-2020年研究区各年桉树人工林空间分布图，其中(a)为2006年研究区桉树人工林空间分布图，(b)为2007年研究区桉树人工林空间分布图，(c)为2008年研究区桉树人工林空间分布图，(d)为2009年研究区桉树人工林空间分布图，(e)为2010年研究区桉树人工林空间分布图，(f)为2011年研究区桉树人工林空间分布图，(g)为2012年研究区桉树人工林空间分布图，(h)为2013年研究区桉树人工林空间分布图，(i)为2014年研究区桉树人工林空间分布图，(j)为2015年研究区桉树人工林空间分布图，(k)为2016年研究区桉树人工林空间分布图，(l)为2017年研究区桉树人工林空间分布图，(m)为2018年研究区桉树人工林空间分布图，(n)为2019年研究区桉树人工林空间分布图，(o)为2020年研究区桉树人工林空间分布图。

步骤4：分别利用N-m至N年研究区各年桉树人工林空间分布对对应的分别从各年研究区影像数据获得的归一化燃烧指数NBR空间分布进行掩膜提取，分别对应得到各年的桉树人工林NBR空间分布，从而获得N-m至N年研究区桉树人工林NBR空间分布；根据各砍伐更新点在桉树人工林NBR空间分布中对应像元的NBR值的变化规律确定NBR骤降幅度阈值。

该步骤的实施过程中，在ArcGIS软件中打开基于步骤3所获得的2006-2020年研究区各年桉树人工林空间分布和从2006-2020年各年研究区影像数据提取的归一化燃烧指数NBR，其中，所提取的2006-2020年各年研究区归一化燃烧指数NBR以空间分布的形式存在，也即从步骤3能获得2006-2020年各年研究区归一化燃烧指数NBR空间分布，利用ArcGIS软件中的掩膜提取功能，利用2006-2020年研究区各年桉树人工林空间分布对从2006-2020年各年研究区影像数据提取的归一化燃烧指数NBR空间分布进行掩膜提取，分别对应得到各年的桉树人工林归一化燃烧指数NBR空间分布，从而得到2006-2020年研究区桉树人工林的归一化燃烧指数NBR空间分布。在ArcGIS中将由步骤1所获得的600个砍伐更新点展开，每个砍伐更新点对应于2006-2020年研究区桉树人工林的归一化燃烧指数NBR空间分布中对应位置的像元值，观察600个砍伐更新点对应于2006-2020年研究区桉树人工林的归一化燃烧指数NBR空间分布中对应像元的NBR值的变化，发现各桉树砍伐更新点对应像元的归一化燃烧指数NBR值呈“平稳-骤降-上升-平稳”的变化，与桉树轮伐期内其他年相比，砍伐更新点对应的当前年的NBR值迅速骤降，如图3中(a)图示出的像元的NBR值的变化如图(b)所示，砍伐更新点对应的2016年的NBR值迅速骤降。本实施方式中多数砍伐更新点对应像元的归一化燃烧指数NBR的骤降幅度在0.45左右，取骤降幅度最小值为0.4，因此，本实施方式将NBR骤降幅度阈值确定为0.4。

步骤5：根据N-m至N年研究区桉树人工林NBR空间分布、NBR骤降幅度阈值，确定N年研究区桉树人工林空间分布中的真实桉树砍伐更新点及对应的砍伐更新年；

方法为：将N-m至N年各年研究区桉树人工林空间分布中所在像元发生NBR骤降对应的时间点定义为桉树伪砍伐更新点；然后以轮伐期为时间步长，在时间步长内：若像元的NBR值只骤降1次，且骤降的幅度大于等于NBR骤降幅度阈值，则该次NBR值骤降对应的伪砍伐更新点则为真实砍伐更新点；若像元的NBR值骤降大于1次，且骤降的幅度均大于等于NBR骤降幅度阈值，则取最后1次NBR值骤降对应的伪砍伐更新点为真实砍伐更新点；各真实砍伐更新点分别对应的当前年定义为各真实砍伐更新点的砍伐更新年。

基于2006-2020年研究区桉树人工林的归一化燃烧指数NBR空间分布，将2006-2020年各年研究区中桉树人工林空间分布中所在像元发生NBR骤降对应的时间点定义为桉树伪砍伐更新点，根据预设时间步长和归一化燃烧指数NBR骤降幅度阈值作去伪处理，确定2020年研究区桉树人工林空间分布中的真实桉树砍伐更新点及对应的砍伐更新年；

在该步骤的实施过程中，将2006-2020年各年研究区桉树人工林空间分布中所在像元发生NBR骤降对应的时间点定义为伪桉树砍伐更新点，各伪桉树砍伐更新点对应的当前年定义为伪桉树砍伐更新点的砍伐更新年。实际上当像元所在区域未发生砍伐事件时，归一化燃烧指数NBR也会下降或者上升，例如由于桉树人工林每年的生长条件不同，雨水光照充足时归一化燃烧指数NBR就会上升，发生干旱或病虫害时归一化燃烧指数NBR就会下降。本发明将像元发生归一化燃烧指数NBR下降时所对应的时间点统一定义为伪砍伐更新点，然后做去伪处理，将砍伐情况以外因素引起NBR下降的错误砍伐更新点去掉，得到真实的砍伐更新点。综合考虑桉树人工林砍伐的具体情况，由于桉树一般在5年砍伐一次，因此本实施方式预设时间步长为5年。在前后共5年内若像元的归一化燃烧指数NBR只骤降1次，且骤降的幅度大于等于归一化燃烧指数NBR的骤降幅度阈值，则该像元对应的伪砍伐更新点则为真实砍伐更新点，骤降发生时对应的当前年为该真实砍伐更新点的砍伐更新年；若在前后共5年内像元的归一化燃烧指数NBR骤降大于1次，且骤降的幅度均大于等于归一化燃烧指数NBR的骤降幅度阈值，则取最后1次NBR骤降对应的伪砍伐更新点为真实砍伐更新点，对应的当前年为真实砍伐更新点的砍伐更新年，其他伪砍伐更新点为错误砍伐更新点。

综上，基于2020年研究区桉树人工林空间分布，将伪砍伐更新点减去错误砍伐更新点即可得到真实砍伐更新点，每个真实砍伐更新点对应一个砍伐更新年，统计2020年研究区桉树人工林空间分布各像元对应的砍伐更新年。

步骤6：基于N年研究区桉树人工林空间分布中各像元对应的砍伐更新年，计算出N年研究区桉树人工林年林龄,进而获得N年研究区桉树人工林年林龄空间分布；

N年研究区桉树人工林年林龄计算方法为：定义各像元内的桉树在砍伐更新年为0岁，利用N年分别减去各像元对应的距离N年最近的砍伐更新年，则得到各像元内桉树在N年的林龄；若N年研究区桉树人工林中部分像元在N-m至N年均无砍伐更新年，则视该部分像元内桉树在N年的林龄为大于m岁。

本实施方式中，基于由步骤5统计2020年研究区桉树人工林空间分布中各像元对应的砍伐更新年，计算出2020年研究区桉树人工林年林龄，进而获得2020年研究区桉树人工林年林龄空间分布。在该步骤的实施过程中，在GEE平台中通过编程计算2020年研究区桉树人工林年林龄，编程思路为：基于2020年研究区桉树人工林空间分布中所在像元(2020年研究区桉树人工林空间分布为栅格图像，因此桉树人工林是以像元的形式存在)，查询各像元对应的真实砍伐更新点的砍伐更新年，明确各像元内的桉树距离N年最近的砍伐更新年。由于真实砍伐更新点对应像元内的桉树在对应砍伐更新年被砍伐，定义各像元内的桉树在砍伐更新年为0岁，因此，在明确各像元对应的距离N年最近的砍伐更新年后，利用2020年分别减去各像元对应的距离N年最近的砍伐更新年，则得到各像元内的桉树的林龄。如：2020年研究区桉树人工林空间分布中某像元对应的距离2020年最近的砍伐更新年为2015年，因此，该像元中的桉树在2015年为0岁，在2020年为5岁。若2020年研究区桉树人工林中部分像元在2006年-2020年均无砍伐更新年，则在2020年视其为大于14岁。综上，可获得2020年研究区桉树人工林空间分布中的所有像元对应桉树的年林龄，进而得到2020年研究区桉树人工林年林龄空间分布。利用上述方法也可获得2007-2019各年研究区桉树人工林年林龄空间分布，本实施方式由于没有使用2006年之前的影像数据，因此，2006年研究区桉树人工林年林龄无法反演。如图4所示为2020年研究区桉树人工林年林龄空间分布。

步骤7：遍历N年研究区桉树人工林空间分布中各个像元对应的距离N年最近的砍伐更新年内所有月份的Landsat影像数据，遍历完成转至步骤11；在遍历过程中，若当前遍历月份的Landsat影像数据不能完全覆盖研究区，存在缺失影像数据，则执行步骤8；若当前遍历月份的Landsat影像数据能完全覆盖研究区，不存在缺失影像数据则转至步骤9；

在本实施例中，基于步骤5统计的2020年研究区桉树人工林空间分布中各像元对应的砍伐更新年，能够确定2020年桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年，遍历2020年研究区各像元对应的距离2020年最近的砍伐更新年内所有月份的Landsat-8影像数据，若当前遍历月份的Landsat-8影像数据能完全覆盖研究区，则执行步骤9；若当前遍历月份的Landsat-8影像数据不能完全覆盖研究区，则执行步骤8。

步骤8：引入当前被遍历的月份及其前后月份的MODIS研究区影像数据，通过STARFM时空融合算法获得与当前遍历月份对应的Landsat影像数据空间分辨率相同的所述缺失影像数据；

在本实施方式中，Landsat-8影像数据空间分辨率高，但时间分辨率低，重访周期为16天，由于卫星影像数据受天气影响，影像数据有时会涵盖大量的云导致不能使用，因此，在同一个月份同一区域Landsat-8只有两张卫星影像，且受云的影像并不一定能使用；MODIS影像数据时间分辨率高，重放周期为1天，同一月份同一区域每天都有一张卫星影像，但MODIS影像空间分辨率低，林龄反演效果相对较差，因此，利用MODIS影像数据与Landsat-8影像数据进行时空融合能够弥补当前遍历月份缺失的Landsat-8影像数据。例如，2020年研究区桉树人工林空间分布中某像元对应的距离2020年最近的砍伐更新年为2020年，而2020年4月份的Landsat-8影像不能够完全覆盖研究区，将2020年研究区3月31日的Landsat-8影像数据、2020年研究区3月31日的MODIS影像数据、2020年研究区5月2日的Landsat-8影像数据，2020年研究区5月2日的MODIS影像数据和2020年研究区4月15日的MODIS影像数据作为输入数据放到STARFM时空融合算法中运行，即可获得研究区2020年4月的具有与Landsat-8相同空间分辨率的能够完全覆盖研究区的影像数据。按上述方法，获得2020年每个月份均可完全覆盖研究区的2020年研究区完整影像数据，从而获得2020年研究区桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年内每个月均可完全覆盖研究区的完整影像数据。

步骤9：在GEE平台上针对当前遍历月份的研究区完整影像数据进行影像的拼接、掩膜提取和辐射归一化预处理，得到预处理后的当前遍历月份影像数据；

在本实施方式中，在GEE平台上对对应像元的砍伐更新年内当前遍历月份的研究区完整影像数据进行影像的拼接、掩膜提取和辐射归一化预处理，得到预处理后的2020年研究区桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年当前遍历月份预影像数据。在该步骤的实施过程中，在GEE平台打开对应像元的距离2020年最近的砍伐更新年内当前遍历月份的研究区完整影像数据进行影像拼接，利用研究区矢量边界图进行掩膜提取，将矢量边界外的部分去掉，再使用GEE自带的辐射归一化算法进行处理，得到预处理后的2020年研究区桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年内当前月份影像数据。

步骤10：在GEE平台上直接提取当前遍历月份研究区影像数据的NBR，利用各像元对应的距离N年最近的砍伐更新年研究区桉树人工林空间分布对当前遍历月份的研究区NBR空间分布进行掩膜提取，得到当前遍历月份的研究区桉树人工林NBR空间分布，返回步骤7；

在本实施方式中，在GEE平台上直接提取2020年研究区桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年当前遍历月份研究区影像数据的NBR。2020年研究区桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年当前遍历月份研究区影像数据的归一化燃烧指数NBR是以空间分布的形式存在，即2020年研究区桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年当前遍历月份研究区归一化燃烧指数NBR空间分布。利用2020年研究区桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年的研究区桉树人工林空间分布对对应像元的砍伐更新年当前遍历月份的研究区归一化燃烧指数NBR空间分布进行掩膜提取，得到2020年研究区桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年当前遍历月份的桉树人工林NBR空间分布。返回步骤7，直至获得2020研究区桉树人工林空间分布中各像元对应的距离2020年最近的砍伐更新年所有月份的研究区桉树人工林归一化燃烧指数NBR空间分布。例如，2020年某像元对应桉树砍伐更新年为2015年，且2015年某一或若干月份研究区Landsat影像数据不能完全覆盖研究区，通过步骤7至步骤10，则能够获得2015年研究区桉树人工林各个月份的研究区桉树人工林归一化燃烧指数NBR空间分布。以此类推，获得2020年研究区桉树人工林空间分布中各像元对应的真实砍伐更新点所对应的距离2020年最近的砍伐更新年每个月的研究区桉树人工林归一化燃烧指数NBR空间分布。

步骤11：根据N年研究区桉树人工林空间分布中各像元在其对应的距离N年最近的砍伐更新年内各月的研究区桉树人工林归一化燃烧指数NBR随时间的变化规律，定义在A年B月刚被砍伐的桉树人工林为0年又0月岁，结合N年研究区桉树人工林年林龄空间分布，可确定在N年研究区桉树人工林空间分布中距离N年最近的砍伐更新年A年内的真实砍伐更新点所对应像元中的被砍伐桉树在N年C月的林龄为Q年L月，其中Q为N与A之间的差值，L为C与B之间的差值。

在本实施方式中，通过观察2020年研究区桉树人工林空间分布中各像元在其对应的距离2020年最近的砍伐更新年内各月的研究区桉树人工林归一化燃烧指数NBR随时间的变化规律，桉树人工林在刚被砍伐时的NBR值最低，随着桉树的生长，其归一化燃烧指数NBR会逐渐上升后趋于平缓，因此，以刚被砍伐的桉树人工林为0年又0月岁。例如：2020年研究区桉树人工林空间分布中某像元对应的距离2020年最近的砍伐更新年为2015年，基于2020年研究区桉树人工林空间分布中各像元对应的真实砍伐更新点所对应的距离2020年最近的砍伐更新年研究区桉树人工林NBR空间分布，可以得到2015年各月研究区桉树人工林NBR空间分布，通过该像元在2015年12个月内归一化燃烧指数NBR值的变化规律，发现该像元在2015年归一化燃烧指数NBR值在5月最低，因此，该像元内的桉树人工林在2015年5月为0年又0月岁，以被砍伐月为基准，往后累计计算桉树人工林月林龄，则该像元内的桉树人工林在2015年12月为0年又7月岁，结合2020年研究区桉树人工林年林龄空间分布，进而计算出该像元内的桉树人工林在2020年12月为5年又7月岁，以此类推，可以将2020年研究区桉树人工林空间分布中所有像元内的桉树人工林的年林龄细化至月林龄，得到2020年研究区桉树人工林月林龄空间分布，进而将2020年桉树人工林年林龄反演精确到月尺度。在利用上述方法也可获得2007-2019各年研究区桉树人工林月林龄空间分布，本实施方式由于没有使用2006年之前的影像数据，因此，2006年研究区桉树人工林年林龄无法反演，以至于2006年研究区桉树人工林月林龄无法反演。如图5所示为2020年研究区桉树人工林月林龄空间分布，其中图例中x/x表示x年又x月岁，例如1/01表示1年又1月岁。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解；其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；因而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。