一种基于多源遥感影像的多云雨地区果树树龄测算方法

技术领域

本发明属于遥感技术领域，具体涉及一种基于多源遥感影像的多云雨地区果树树龄测算方法。

背景技术

树龄是评估果树生长状态的重要指标，果树的树龄测算对果树种植管理、产量品质预估、区域种植规划等具有重要作用。传统的树龄测算主要采取实地调查，费时费力，成本较高，且存在以“点”代“面”的问题，难以开展区域尺度的果园树龄统计。

近年来，卫星遥感技术为区域果树树龄测算提供了新的思路，现有技术中也出现了一些通过遥感技术对树龄进行测算的方法。例如，申请号为202010266245.4的中国专利《一种果树树龄识别方法、装置、设备及存储介质》以及申请号为202011610466.5的中国专利《一种基于密集时间遥感数据的新生红树林树龄评估方法》。

上述专利提出的方法均可实现基于卫星遥感影像进行树龄测算，但是要求采用多时相的中分辨率密集时序可见光影像(Lansat-TM或Sentinel-2A/2B)。然而在南方多云雨地区，受天气影响，可见光遥感影像缺失严重，难以获取密集时序的中分辨率遥感影像。现有技术中，能够实现密集时序观测的卫星通常空间分辨率较低，例如：Terra和Aqua卫星搭载的MODIS传感器每1-2天即可重复观测整个地球表面，可提供连续的植被指数观测数据，并通过连续16天的数据合成减轻了云雨天气的影响，可形成南方多云雨地区密集时序的植被指数数据；但由于其空间分辨率最高仅为250米，对于南方多云雨地区地块破碎的立地条件，空间分辨率太低，难以准确反映各个果园地块的植被指数。因此，上述树龄测算方法难以在多云雨地区开展应用。

此外，申请号为201910759290.0的中国专利《一种基于时间序列遥感数据的油棕树龄测算方法》提出了一种基于时间序列遥感数据的油棕树龄测算方法，利用遥感影像中时间序列遥感卫星数据，结合用于检测森林扰动的LandTrendr模型，得到油棕特许经营权范围内的森林扰动年份，进而推算出油棕的树龄。该方法利用年际间植被指数的变动进行树龄测算，对时序影像密集程度要求相对降低。但在南方多云雨地区，植被生长速度较快，且果园与农作物争地的现象普通存在。果园改种其他作物后，所在地块的植被指数仍能迅速增长。若以年为时间单位进行树龄测算，时间跨度较长，对于南方果树的树龄测算存在较大的不确定性。由此可见，该树龄测算方法同样难以在多云雨地区开展应用。

发明内容

本发明的目的是提供一种基于多源遥感影像的多云雨地区果树树龄测算方法，以解决现有基于遥感技术的树龄测算方法难以在多云雨地区开展应用的问题。

为达到上述目的，本发明所采用的技术方案是：

一种基于多源遥感影像的多云雨地区果树树龄测算方法，包括以下步骤：

步骤1：获取目标区域最新时相的多源高分辨率遥感影像进行处理，得到无云影的高分辨率影像产品；

步骤2：对步骤1得到的高分辨率影像产品进行地块边界提取和果园地块识别，得到目标区域地块边界矢量图和果园地块边界矢量图；

步骤3：获取目标区域密集时序的历年中分辨率植被指数，获取方法包括：

方法一：获取目标区域密集时序的历年中分辨率遥感影像进行处理，得到密集时序的历年中分辨率植被指数；

方法二：包括以下步骤，

步骤31：获取目标区域时序间断的历年中分辨率遥感影像进行处理，并分年度合成中分辨率植被指数，得到历年年度中分辨率植被指数；

步骤32：获取目标区域密集时序的历年低分辨率遥感影像进行处理，得到密集时序的历年低分辨率植被指数；

步骤33：以步骤31得到的历年年度中分辨率植被指数为基准，对步骤32得到的历年低分辨率植被指数进行分辨率提升，得到密集时序的历年中分辨率植被指数；

步骤4：通过步骤2得到的果园地块边界矢量图和步骤3得到的历年中分辨率植被指数，获得密集时序的果园地块历年植被指数；

步骤5：根据果园的植被指数变化特征，对步骤4获得的果园地块历年植被指数进行断点检测，获取断点位置；

步骤6：根据步骤5获取的断点位置，得到果园地块历年植被指数出现断点的最新时间，通过计算当前时间与出现断点最新时间的年度差距，得到果树的树龄。

进一步地，步骤1中，对多源高分辨率遥感影像进行处理包括几何纠正、辐射校正、大气纠正、云影检测、影像融合、镶嵌、裁剪。

进一步地，步骤2中，地块边界提取和果园地块识别采用深度学习算法。

进一步地，步骤31中，对历年中分辨率遥感影像进行处理，包括：

步骤311：对历年中分辨率遥感影像进行预处理，获得预处理影像，其中，预处理包括几何纠正、辐射校正、大气纠正处理；

步骤312：利用步骤311获得的预处理影像中的波段计算植被指数；

步骤313：将同一年度的中分辨率遥感影像进行植被指数合成，得到历年年度中分辨率植被指数。

进一步地，步骤312中，植被指数为归一化植被指数NDVI，步骤313中，植被指数合成采用最大值合成方法。

进一步地，步骤33中，分辨率提升的方法包括以下步骤：

步骤331：将历年低分辨率植被指数的影像插值为中分辨率影像；

步骤332：将步骤331得到的中分辨率影像与步骤31得到的历年年度中分辨率植被指数的影像进行融合。

进一步地，步骤331中，插值方法包括双线性插值法、cubic插值法。

进一步地，步骤332中，融合方法为基于像素的直接加权平均方法，其设置一定的融合窗口大小，并对每个像素在窗口大小范围内进行加权处理；其中，权重根据历年年度中分辨率植被指数中的NDVI值进行分配，具体算法如下：

上述算法中，n为奇数，是设定的融合窗口大小；L(i,j)是同年的年度中分辨率植被指数的影像中第i行j列的NDVI值，M(i,j)是经过插值的中分辨率影像中第i行j列的NDVI值，P(i,j)是第i行j列的权重，N(i,j)是经过融合计算得到的第i行j列的新值。

进一步地，步骤5中，在进行断点检测之前，可以对数据进行以下预处理步骤：

数据清洗：检查数据中是否存在错误、异常值或无效数据；

数据平滑：如果数据存在噪声或震荡，可以应用平滑技术来减少这些干扰，平滑方法包括移动平均、指数平滑、Loess平滑；

数据插值：如果数据中存在缺失值，可以使用插值技术填补这些缺失值，以使数据连续，插值方法包括线性插值、多项式插值、样条插值。

进一步地，步骤5中，断点检测采用滑动窗口方法。

由于采用上述技术方案，本发明具有以下有益效果：

本发明综合利用不同分辨率遥感影像的优势，解决了南方多云雨地区果树树龄测算难题。首先，本发明通过对目标区域最新时相的多源高分辨率遥感影像进行处理，得到无云影的高分辨率影像产品，避免云雨天气对测算结果造成影响。其次，本发明利用时序间断的中分辨率影像和密集时序的低分辨率影像，形成密集时序的历年中分辨率植被指数，并进一步形成密集时序的果园地块历年植被指数，通过对分辨率的提升，可以准确反映各个果园地块的植被指数，进而提高测算结果的准确性。最后，本发明通过对果园地块的历年植被指数进行断点检测，由此测算果树树龄，减少测算过程中的不确定性，进一步提高测算结果的准确性。

附图说明

图1是本发明提供的多云雨地区果树树龄测算方法的流程图；

图2是本发明提供的多云雨地区果树树龄测算方法的技术路线图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，基于本发明中的实施例，本领域技术人员在没有做出创造性劳动前提下所获得的其他实施例，均属于本发明保护的范围。

实施例

请参阅图1-2，一种基于多源遥感影像的多云雨地区果树树龄测算方法，包括步骤1-6。

步骤1：获取目标区域最新时相的多源高分辨率遥感影像进行处理，得到无云影的高分辨率影像产品。其中，对多源高分辨率遥感影像进行处理包括几何纠正、辐射校正、大气纠正、云影检测、影像融合、镶嵌、裁剪。

南方地区云雨天气较多，卫星拍摄的影像上通常会有云影遮挡，造成部分地面无法被观测到，解决的方法是采用多景影像进行拼接，将各景影像中无云影的部分合并成一景，得到完全没有遮挡的地面观测影像，所用到的多景影像可以是同一颗卫星在不同时段拍摄的影像，也可以是不同卫星(多源)在不同时段拍摄的影像。如果能找到一景完全没有云影遮挡的影像，只用一景也可以。

具体案例中，为了提取广西壮族自治区武鸣县双桥镇2023年的柑橘种植地块，获取了高分一号卫星2023年7月、资源三号卫星2023年8月和北京二号卫星2023年9月各一景遥感影像，对获取的遥感影像进行预处理，包括几何纠正、辐射校正、大气纠正、云影检测、影像融合、镶嵌、裁剪等，最终形成一景目标区域无云影的2米分辨率影像产品。

步骤2：对步骤1得到的高分辨率影像产品进行地块边界提取和果园地块识别，得到目标区域地块边界矢量图和果园地块边界矢量图。其中，地块边界提取和果园地块识别采用深度学习算法。

具体案例中，基于上述2米分辨率影像产品，采用深度学习算法进行地块边界提取和柑橘地块识别，得到目标区域地块边界矢量图和柑橘种植地块边界矢量图。其中，地块边界提取和柑橘地块识别有多种方法，案例中采用的深度学习算法作为一种优选方法，通常通过人工目视解译也可以进行提取和识别，但比较费时费力。

步骤3：获取目标区域密集时序的历年中分辨率植被指数。其中，获取方法包括方法一和方法二。若方法一能够实现，优先选择方法一，若方法一无法实现，则采用方法二，在南方多云雨地区，获取密集时序的遥感影像可能相对困难，因此多数情况下会使用方法二。

方法一：获取目标区域密集时序的历年中分辨率遥感影像进行处理，得到密集时序的历年中分辨率植被指数。

方法二：包括步骤31-33。

步骤31：获取目标区域时序间断的历年中分辨率遥感影像进行处理，并分年度合成中分辨率植被指数，得到历年年度中分辨率植被指数。

其中，对历年中分辨率遥感影像进行处理，包括步骤311-313。

步骤311：对历年中分辨率遥感影像进行预处理，获得预处理影像，其中，预处理包括几何纠正、辐射校正、大气纠正处理。

步骤312：利用步骤311获得的预处理影像中的波段计算植被指数。本实施例中，植被指数为归一化植被指数NDVI，其计算公式如下：

NDVI＝(NIR-R)/(NIR+R)

其中，NIR为近红外波段数值，R为红光波段数值。

具体案例中，获取了2001-2023年覆盖广西壮族自治区武鸣县双桥镇的30米分辨率Lansat-TM影像总共79景，各个年度通常只有2-5景影像，其中2004年最多有8景。经过几何纠正、辐射校正、大气纠正等影像预处理后，利用影像中的近红外波段和红光波段计算各景的NDVI。

在其他实施例中，植被指数也可以是其他的指数，例如归一化燃烧比率NBR，归一化水体指数NDWI。

步骤313：将同一年度的中分辨率遥感影像进行植被指数合成，得到历年年度中分辨率植被指数。通过合成处理，可有效减少影像中云影的干扰。本实施例中，植被指数合成采用最大值合成方法。最大值合成，即某一位置的像元值采用同一年度中所有影像同一位置中的最大值。在其他实施例中，合成方法还可以是平均值合成、总和合成等方法。

具体案例中，2023年只有3景影像，拍摄时间分别是10月16日、11月17日和12月08日，将这3景影像计算的NDVI进行最大值合成得到2023年度Landsat30米分辨率NDVI合成影像。通过相同方法，得到2001-2023各年的年度Landsat30米分辨率NDVI合成影像。

步骤32：获取目标区域密集时序的历年低分辨率遥感影像进行处理，得到密集时序的历年低分辨率植被指数。

具体案例中，获取了2001-2023年每隔16天一期的MODIS MOD13Q1影像产品总共529景，其中每年都有23景。该影像产品包含有16天最大值合成的NDVI数据，空间分辨率为250米。通过ENVI软件读取其中的NDVI数据，并将其进行裁剪，得到目标区域2001-2023年每隔16天一期的密集时序的NDVI数据。

步骤33：以步骤31得到的历年年度中分辨率植被指数为基准，对步骤32得到的历年低分辨率植被指数进行分辨率提升，得到密集时序的历年中分辨率植被指数。

其中，分辨率提升的方法包括步骤331和步骤332。

步骤331：将历年低分辨率植被指数的影像插值为中分辨率影像。插值方法包括双线性插值法、cubic插值法，除此之外还可以采用其他方法。本实施例中，采用双线性插值法，将MODIS250米分辨率的NDVI影像插值为30米分辨率，形成MODIS 30米分辨率NDVI影像(即中分辨率影像)。经过插值后形成的30米分辨率的NDVI影像中，像元的NDVI值继承了250米分辨率图像中相应位置像元的NDVI值，仍难以反映各地块间的植被指数差异，因此需要进行下一步的融合处理。

步骤332：将步骤331得到的中分辨率影像与步骤31得到的历年年度中分辨率植被指数的影像进行融合。融合方法为基于像素的直接加权平均方法，其设置一定的融合窗口大小，并对每个像素在窗口大小范围内进行加权处理；其中，权重根据历年年度中分辨率植被指数中的NDVI值进行分配，具体算法如下：

上述算法中，n为奇数，是设定的融合窗口大小；L(i,j)是同年的年度中分辨率植被指数的影像中第i行j列的NDVI值，M(i,j)是经过插值的中分辨率影像中第i行j列的NDVI值，P(i,j)是第i行j列的权重，N(i,j)是经过融合计算得到的第i行j列的新值。

本实施例中，是将插值后的MODIS 30米分辨率NDVI影像与步骤31形成的历年年度Landsat30米分辨率NDVI合成影像进行融合，其中权重根据历年年度Landsat30米分辨率NDVI合成影像中的NDVI值进行分配。具体算法中，L(i,j)是同年的年度Landsat30米分辨率NDVI合成影像中第i行j列的值，M(i,j)是经过插值的MODIS 30米分辨率NDVI影像中第i行j列的值。

具体案例：如表1所示，经过插值的MODIS 30米分辨率NDVI影像中，第5行第5列像元的值为M55。

表1.MODIS 30米像元

如表2所示，在年度Landsat30米分辨率NDVI合成影像中，相同位置及其周围相邻位置像元值分别为L11、L12、…、L99。

表2.Landsat 30米像元

选定插值的窗口大小为9，即窗口是9*9个像元范围，则第5行第5列像元的融合算法为：

P55＝L55/(L11+L12+…+L99)

N55＝P55*(M11+M12+…+M99)

其中P55是第5行第5列像元的权重，N55是第5行第5列像元的新值。

通过这种方法，将MODIS250米分辨率的NDVI影像提升至30米分辨率，并尽可能保留了年度Landsat30米分辨率NDVI合成影像的纹理信息。此外，采用此方法，以各年的年度Landsat30米分辨率NDVI合成影像为权重的基准，对各年度中所有的MODIS250米分辨率的NDVI图进行分辨率提升(2001-2023年每年有23景)，得到密集时序的30米分辨率NDVI图(即密集时序的历年中分辨率植被指数)。

步骤4：通过步骤2得到的果园地块边界矢量图和步骤3得到的历年中分辨率植被指数，获得密集时序的果园地块历年植被指数。

具体案例中，利用步骤2获取的果园地块边界矢量图和步骤3获取的密集时序30米分辨率NDVI图，统计各个时段各个果园地块边界内的NDVI平均值，得到各个果园地块的密集时序NDVI值(即密集时序的果园地块历年植被指数)。其中，数值统计在ARCMAP软件中进行。

步骤5：根据果园的植被指数变化特征，对步骤4获得的果园地块历年植被指数进行断点检测，获取断点位置。

在进行断点检测之前，可以对数据进行以下预处理步骤：

数据清洗：检查数据中是否存在错误、异常值或无效数据；

数据平滑：如果数据存在噪声或震荡，可以应用平滑技术来减少这些干扰，平滑方法包括移动平均、指数平滑、Loess平滑，除此之外，还可以采用其他方法；

数据插值：如果数据中存在缺失值，可以使用插值技术填补这些缺失值，以使数据连续，插值方法包括线性插值、多项式插值、样条插值，除此之外，还可以采用其他方法。

例如，进行断点检测前，采用移动平均平滑方法对数据进行平滑处理，以进一步减轻云影对数据的干扰。移动平均平滑是一种常见的数据平滑方法，用于减少随机噪声对时序数据的干扰。该方法的主要思想是利用一定窗口大小内的数据点的平均值来替代原始数据点。本实施例中，窗口大小设置为5。

断点检测采用滑动窗口方法，滑动窗口方法是将时序数据分割成固定长度的窗口，然后通过计算窗口内的统计特征(如均值、方差等)来判断窗口内是否存在断点。常见的滑动窗口方法有滑动平均法和滑动方差法。本实施例采用滑动平均法。

本实施例中，使用IDL语言编写了detect_breakpoints函数来实现时序断点检测的滑动平均算法。函数接收三个参数：data表示输入数据序列，window_size表示滑动窗口的大小，threshold表示判断断点的阈值。根据MODIS数据序列的特点，设置滑动窗口的大小为23，阈值为0.5。

函数首先通过遍历数据序列，每次取出一个大小为window_size的滑动窗口，然后计算窗口内数据的平均值；接着，函数判断窗口后的一个数据是否与窗口内的平均值之差超过了给定的阈值threshold。如果超过了阈值，说明该数据点为断点，将断点位置添加到breakpoints数组中。最后，函数返回所有检测到的断点位置。

具体代码如下：

FUNCTION detect_breakpoints,data,window_size,threshold

breakpoints＝[]

FOR i＝0,N_ELEMENTS(data)-window_size-1DO BEGIN

window＝data[i:i+window_size-1]

avg＝TOTAL(window)/window_size

IF ABS(data[i+window_size]-avg)GT threshold THEN BEGIN

breakpoints＝[breakpoints,i+window_size]

ENDIF

ENDFOR

RETURN,breakpoints

END

步骤6：根据步骤5获取的断点位置，得到果园地块历年植被指数出现断点的最新时间，通过计算当前时间与出现断点最新时间的年度差距，得到果树的树龄。

具体案例中，某像元的断点时间出现在2002、2014和2019年，最近的断点时间为2019年，影像的最新日期为2023年，则树龄为2023-2019＝4年。

本发明综合利用不同分辨率遥感影像的优势，解决了南方多云雨地区果树树龄测算难题。首先，本发明通过对目标区域最新时相的多源高分辨率遥感影像进行处理，得到无云影的高分辨率影像产品，避免云雨天气对测算结果造成影响。其次，本发明利用时序间断的中分辨率影像和密集时序的低分辨率影像，形成密集时序的历年中分辨率植被指数，并进一步形成密集时序的果园地块历年植被指数，通过对分辨率的提升，可以准确反映各个果园地块的植被指数，进而提高测算结果的准确性。最后，本发明通过对果园地块的历年植被指数进行断点检测，由此测算果树树龄，减少测算过程中的不确定性，进一步提高测算结果的准确性。

上述说明是针对本发明较佳可行实施例的详细说明，但实施例并非用以限定本发明的专利申请范围，凡本发明所提示的技术精神下所完成的同等变化或修饰变更，均应属于本发明所涵盖专利范围。