CN115830446A

CN115830446A - 一种动态水体产品融合方法、装置、设备及可读存储介质

Info

Publication number: CN115830446A
Application number: CN202211491508.7A
Authority: CN
Inventors: 王义成; 饶品增; 王兴勇; 蒲立; 王燕晓; 赵月芬; 吕映; 甘醇; 焦醒; 范一为
Original assignee: China Institute of Water Resources and Hydropower Research
Current assignee: China Institute of Water Resources and Hydropower Research
Priority date: 2022-11-25
Filing date: 2022-11-25
Publication date: 2023-03-21
Anticipated expiration: 2042-11-25
Also published as: CN115830446B

Abstract

本发明公开了一种动态水体产品融合方法、装置、设备及可读存储介质。该方法包括获取多个动态MODIS产品数据集；获取水体和非水体样本，对各个动态MODIS产品样本数据集进行特征提取，训练水体分类模型；利用训练后的水体分类模型对动态MODIS产品数据集进行水体分类；对不同动态MODIS产品数据集的水体分类结果进行噪声处理；采用时空自适应融合方法对噪声处理后的不同动态MODIS产品数据集的水体分类结果进行融合，得到融合的水体分类结果。本发明采用基于同源数据的时空自适应融合方法对两组产品的分类结果进行融合，从而降低了分类误差和噪声等对结果的影响，也提高了分类结果的时空分辨率。

Description

一种动态水体产品融合方法、装置、设备及可读存储介质

技术领域

本发明涉及湖泊水体提取技术领域，具体涉及一种动态水体产品融合方法、装置、设备及可读存储介质。

背景技术

湖泊与人类生活和自然环境密切相关。准确认识湖泊的长期或动态变化至关重要。遥感卫星可以提供监测所必需的大量数据。陆地水体的长期动态已经被绘制出来。然而，对陆地湖泊水体进行高空间分辨率的短期连续监测仍然是一项艰巨的任务。适当的数据来源和可靠的方法都是长期监测水体的关键。

MODIS数据具有较高的时间和中等空间分辨率，是湖泊水体变化的一个很好的数据源。然而，高密度的光学遥感数据容易受到噪声的影响，这些噪声会对最终结果造成很大的误差。为了获得连续的地表水，修复被云、云影等噪声覆盖的水是至关重要的。

除此之外，一个合适的分类器对于长期的水提取是非常有用的。水提取方法比较成熟，大致可以分为采用一定指标的阈值法、有监督和无监督分类法以及其他一些方法。指数法(NDWI)简单可行，但时间序列水体提取难度较大，阈值需要在不同季节不断调整。监督分类方法是可行的，并且哪一种不需要确定阈值，只需要提供一些可靠的样本点。一些研究表明，随机森林(Random Forest,RF)分类器具有较好的鲁棒性，对土地覆盖和地表水具有高效的分类能力。

与单一数据源相比，联合使用高效的多个数据源有助于避免噪声，获得更好的结果。一些针对多数据源的数据同化和融合方法得到了很好的应用，如时空自适应反射率融合模型(STARFM)和时空反射率解混模型(STRUM)。然而，由于湖泊“水体”时间序列具有显著的动态特性，数据融合仍然是一个挑战。由于传感器本身和天顶角的不同，不同数据源的融合有显著差异。此外，由于误差重叠，使用融合产品可能会产生更大的误差。

发明内容

针对现有技术中的上述不足，本发明提供了一种动态水体产品融合方法、装置、设备及可读存储介质。

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

第一方面，本发明提出了一种动态水体产品融合方法，包括以下步骤：

S1、获取多个动态MODIS产品数据集；

S2、获取水体和非水体样本，对各个动态MODIS产品样本数据集进行特征提取，训练水体分类模型；

S3、利用训练后的水体分类模型对动态MODIS产品数据集进行水体分类；

S4、对不同动态MODIS产品数据集的水体分类结果进行噪声处理；

S5、采用时空自适应融合方法对噪声处理后的不同动态MODIS产品数据集的水体分类结果进行融合，得到融合的水体分类结果。

可选地，步骤S4具体包括：

对不同动态MODIS产品数据集的水体分类结果，剔除含云部分的分类结果。

可选地，步骤S5具体包括：

S51、根据第一动态MODIS产品数据集的第一时间间隔差和第二动态MODIS产品数据集的第二时间间隔差确定融合周期；

S52、将第一动态MODIS产品数据集在融合周期内第一空间分辨率的第一水体分类结果和第二水体分类结果进行像元分割，得到第二空间分辨率的第一像元分割结果和第二像元分割结果；

S53、根据第二动态MODIS产品数据集在融合周期内第二空间分辨率的第二水体分类结果获取对应的生产日期像元结果；

S54、依次遍历生产日期像元结果中的每个像元，根据该像元值所处值域区间以及第二水体分类结果中该像元值对应位置的像元值与第一像元分割结果中该像元值对应位置的像元值大小关系，确定第一时间间隔差下第二空间分辨率的第一水体融合分类结果；

S55、依次遍历生产日期像元结果中的每个像元，根据该像元值所处值域区间以及第二水体分类结果中该像元值对应位置的像元值与第二像元分割结果中该像元值对应位置的像元值大小关系，确定第一时间间隔差下第二空间分辨率的第二水体融合分类结果。

可选地，步骤S52具体包括：

将第一动态MODIS产品数据集在融合周期内第一空间分辨率的第一水体分类结果中每个像元分割为四个像元，并将分割后的四个像元值赋值为原像元值，得到第二空间分辨率的第一像元分割结果；

将第一动态MODIS产品数据集在融合周期内第一空间分辨率的第二水体分类结果中每个像元分割为四个像元，并将分割后的四个像元值赋值为原像元值，得到第二空间分辨率的第二像元分割结果。

可选地，步骤S54具体包括以下分步骤：

S541、依次遍历生产日期像元结果中的每个像元，判断该像元值所处值域区间处于第一值域区间还是第二值域区间；

若该像元值处于第一值域区间，则将第二水体分类结果中对应像元值作为第一水体融合分类结果的对应像元值；

若该像元值处于第二值域区间，则执行步骤S542；

S542、判断第二水体分类结果中该像元值对应位置的像元值是否等于第一像元分割结果中该像元值对应位置的像元值；若是，则执行步骤S543；否则执行步骤S544；

S543、将第二水体分类结果中该像元值对应位置的像元值或第一像元分割结果中该像元值对应位置的像元值作为第一水体融合分类结果的对应像元值；

S544、分别计算第二水体分类结果中对应位置像元的梯度值和第一像元分割结果中对应位置像元的梯度值；

S545、判断第一像元分割结果中对应位置像元的梯度值是否小于第二水体分类结果中对应位置像元的梯度值；若是，则将第一像元分割结果中该像元值对应位置的像元值作为第一水体融合分类结果的对应像元值；否则将第二水体分类结果中该像元值对应位置的像元值作为第一水体融合分类结果的对应像元值。

可选地，步骤S55具体包括以下分步骤：

S551、依次遍历生产日期像元结果中的每个像元，判断该像元值所处值域区间处于第一值域区间还是第二值域区间；

若该像元值处于第二值域区间，则将第二水体分类结果中对应像元值作为第二水体融合分类结果的对应像元值；

若该像元值处于第一值域区间，则执行步骤S552；

S552、判断第二水体分类结果中该像元值对应位置的像元值是否等于第二像元分割结果中该像元值对应位置的像元值；若是，则执行步骤S543；否则执行步骤S544；

S543、将第二水体分类结果中该像元值对应位置的像元值或第二像元分割结果中该像元值对应位置的像元值作为第二水体融合分类结果的对应像元值；

S544、分别计算第二水体分类结果中对应位置像元的梯度值和第二像元分割结果中对应位置像元的梯度值；

S545、判断第二像元分割结果中对应位置像元的梯度值是否小于第二水体分类结果中对应位置像元的梯度值；若是，则将第二像元分割结果中该像元值对应位置的像元值作为第二水体融合分类结果的对应像元值；否则将第二水体分类结果中该像元值对应位置的像元值作为第二水体融合分类结果的对应像元值。

可选地，像元的梯度值的计算方法为：

其中，gradient_i为像元i的梯度值，T为像元i周围像元个数，dx_it为像元i与周围像元t的距离，dy_it为像元i与周围像元t的差值。

第二方面，本发明提出了一种动态水体产品融合装置，包括：

数据获取模块，用于获取多个动态MODIS产品数据集；

模型训练模块，用于获取水体和非水体样本，对各个动态MODIS产品样本数据集进行特征提取，训练水体分类模型；

水体分类模块，用于利用训练后的水体分类模型对动态MODIS产品数据集进行水体分类；

噪声处理模块，用于对不同动态MODIS产品数据集的水体分类结果进行噪声处理；

数据融合模块，用于采用时空自适应融合方法对噪声处理后的不同动态MODIS产品数据集的水体分类结果进行融合，得到融合的水体分类结果。

第三方面，本发明提出了一种动态水体产品融合设备，包括：存储器和处理器；

所述存储器，用于存储程序；

所述处理器，用于执行所述程序，实现如上述的动态水体产品融合方法的各个步骤。

第四方面，本发明提出了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时，实现如上述的动态水体产品融合方法的各个步骤。

本发明具有以下有益效果：

由于云、云影等噪声的广泛存在，利用遥感数据进行持续、高时空分辨率的湖泊水位动态监测仍有局限性，因此本发明提出了一种融合方法来绘制高精度连续时间序列的湖泊水域。为了减少噪声的影响，保证最终结果的高时空分辨率，采用MODIS两套产品(包括MOD09A1和MOD13Q1)进行水体提取。分类器采用随机森林(Random Forest,RF)方法，获得初步的分类结果。采用基于同源数据的时空自适应融合方法(Homologous Data-BasedSpatial and Temporal Adaptive Fusion Method,HDSTAFM)对两组产品的分类结果进行融合，从而降低了分类误差和噪声等对结果的影响，也提高了分类结果的时空分辨率。

附图说明

图1为本发明实施例1中一种动态水体产品融合方法的流程示意图；

图2为本发明实施例1中动态水体产品融合过程示意图；

图3为本发明实施例1中网格像元融合示意图；

图4为本发明实施例1中某湖泊及周边地区的分类结果，并和卫星测高数据进行对比示意图；

图5为本发明实施例1中某湖泊每8天的分类结果对应的水体面积和湖口站断面的实测水位的相关关系图；

图6为本发明实施例2中一种动态水体产品融合装置的结构示意图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

实施例1

如图1至图3所示，本发明实施例提供了一种动态水体产品融合方法，包括以下步骤S1至S5：

S1、获取多个动态MODIS产品数据集；

在本发明的一个可选实施例中，本实施例选取数据集样本主要参考MOD44W产品和对应的真假彩色影像。为了选择准确的样本，参考湖泊区域的MOD44W产品和对应的MODIS真假彩色产品(MOD09A1或MOD13Q1)，去除有云的像元，肉眼判断水体和非水体，分别标记为1和0，并记录对应的日期。

MOD09A1是一组地表反射率产品，包含1-7波段，分辨率约为500m，为8天网格产品。每一个像元都与8天内的最佳观测有关，这些观测基于高观测覆盖率、低视角、无云或云影和气溶胶装载。它还提供包含MOD35云/雪/冰标志的质量波段(QA)数据，以及所有像元对应的生产日期(DOY)。

MOD13Q1主要用于植被监测。它包含NDVI和EVI，以及红色、近红外、蓝色和中红外波段，分别对应波段1、波段2、波段3、波段7。该产品每16天提供一次，空间分辨率约为250米。此外，它还提供质量波段(QA)和DOY相关数据。

MOD44W产品是全球地表水数据集，空间分辨率为250m。该产品提供了相对可靠的地表水范围，可用于训练分类模型。值得注意的是，这些MODIS产品是从一级大气档案和分布(LAADS,https://ladsweb.modaps.eosdis.nasa.gov)下载的。

在本发明的一个可选实施例中，本实施例结合遥感影像，人工获取湖区的水体和非水体样本点，来对动态MODIS产品的样本数据集进行特征提取。

MOD09A1包含7个波段，为了使得分类结果更加可靠，本实施例计算三个水体相关的指数：NDVI、NDWI和MNDWI。总共组成10个特征用于湖泊周边水体和非水体分类。MOA13Q1包含4个波段和两个植被指数(NDVI和EVI)，均被用做分类特征。

根据提取的分类特征和样本位置，记录每个样本对应日期的所有特征值，构建成一套水体—非水体分类样本，并按照7：3的比例区分训练样本和测试样本。采用随机森林算法和训练样本构建分类模型，并进行分类。

在本发明的一个可选实施例中，本实施例利用训练后的水体分类模型和所有特征变量分类整个湖泊区域的所有像元，得到湖泊及周边的水体—非水体分类结果。

在本发明的一个可选实施例中，步骤S4具体包括：

对不同动态水体产品数据集的水体分类结果，剔除含云部分的分类结果，标记为-1。

本实施例针对两套MODIS产品(MOD09A1和MOD13Q1)去噪后的结果为两套不同时空分辨率的水体产品(8天500m和16天250米)。

在本发明的一个可选实施例中，步骤S5具体包括：

S52、将第一动态MODIS产品数据集在融合周期内第一空间分辨率的第一水体分类结果和第二水体分类结果进行像元分割，得到第二空间分辨率的第一像元分割结果和第二像元分割结果；具体包括：

S54、依次遍历生产日期像元结果中的每个像元，根据该像元值所处值域区间以及第二水体分类结果中该像元值对应位置的像元值与第一像元分割结果中该像元值对应位置的像元值大小关系，确定第一时间间隔差下第二空间分辨率的第一水体融合分类结果；具体包括以下分步骤：

若该像元值处于第二值域区间，则执行步骤S542；

S55、依次遍历生产日期像元结果中的每个像元，根据该像元值所处值域区间以及第二水体分类结果中该像元值对应位置的像元值与第二像元分割结果中该像元值对应位置的像元值大小关系，确定第一时间间隔差下第二空间分辨率的第二水体融合分类结果；具体包括以下分步骤：

若该像元值处于第一值域区间，则执行步骤S552；

具体而言，本实施例中根据步骤S1至S4可以基于MOD09A1产品获取8天500米的地表水体分类结果；基于MOD13Q1产品获取16天250米的地表水体分类结果。在这种情况下，为了获得更高时空分辨率的地表水体数据，本实施例采用时空自适应融合方法来融合这两套产品结果。其主要思路为：考虑同一个时间段内(16天)的分类结果进行融合。在16天内，MOD09A1有两期分类产品(假设为2010001和2010009天的分类结果)，MOD13Q1有一套分类产品(对应的产品日期为2010001)。对于前者来说，2010001天的分类结果是基于2010001-2010008天的最佳像元的分类结果，而201009天的分类结果为2010009-2010016天的最佳像元的分类结果；对于后者来说，2010001天的分类结果为2010001-2010016天的最佳像元的分类结果。在这种情况下，可以考虑将以16天为一个周期来融合两套分类结果。

首先，将2010001天和2010009天的MOD09A1分类结果数据(500m空间分辨率)直接分割成250m，即将一个像元分割为四个像元。像元值(0或1)没有改变；也就是说，如果初始像元值为1，则4个像元(通过除法得到)的值也等于1。他们分割后的像元数量和MOD13Q1在2010001天的分类结果像元数量完全相同。图3展示了部分示例像元的水体和非水体数值。A1是分割后的MOD09A12010001天的任一像元。对应的位置为分割后MOD09A1 2010009天的A2像元和MOD13Q1 2010001天的B像元。C像元为B像元在2010001天对应的生产日期(Theday of The year,DOY)，取值范围为1～16。D1和D2分别为需要判断的2010001天和2010009天的融合像元值。

然后比较MOD09A1和MOD13Q1在2010001天的结果，存在以下几种情况：

如果1≤C<9，则说明B是准确的，D1的数值直接赋值给B，即MOD13Q1的结果(即16天250m分辨率的最优值)比MOD09A1的结果(即8天500m分辨率的最优值)更加准确。

如果9≤C<17，且A1＝B，则说明两个分类结果完全一致，D1的值可以分配给A1或B。

如果9≤C<17，且A1≠B，则需要比较它们与周围像元的相关性，以确定最终结果。本实施例采用一个梯度方程来计算他们的数值和周边像元相关性。对比A1和B的梯度值，如果gradient_A1<gradient_B，则说明A1与窗口内周围像元的相关性比B更合理，则D1＝A1；否则D1＝B。

像元的梯度值的计算方法为：

其中，gradient_i为像元i的梯度值，T为像元i周围像元个数，dx_it为像元i与周围像元t的距离，dy_it为像元i与周围像元t的差值。如果像元窗口设置为3×3，则T＝8，此时dx_it＝1或

和dy_it＝0或±1。

本实施例中D1和D2的融合规则可以表示为：

通过上述分类结果的时空自适应融合，可以计算得2010001天和2010009天的所有对应的像元值(D1和D2)。该结果具有250m的空间分辨率，时间间隔差为8天。根据以上规则，我们可以获取得到每8天250米的水体分类结果。

如图4所示，为某湖泊及周边地区的分类结果，并和卫星测高数据进行对比示意图。如图5所示，为某湖泊每8天的分类结果对应的水体面积和湖口站断面的实测水位的相关关系图。

实施例2

如图6所示，本发明实施例基于实施例1所描述的一种动态水体产品融合方法，提供了一种动态水体产品融合装置，包括：

数据获取模块，用于获取多个动态MODIS产品数据集；

本发明实施例2提供的一种动态水体产品融合装置具有上述实施例1中一种动态水体产品融合方法的有益效果。

实施例3

本发明实施例基于实施例1所描述的一种动态水体产品融合方法，提供了一种动态水体产品融合设备，包括：存储器和处理器；

所述存储器，用于存储程序；

本发明实施例3提供的一种动态水体产品融合设备具有上述实施例1中一种动态水体产品融合方法的有益效果。

实施例4

本发明实施例基于实施例1所描述的一种动态水体产品融合方法，提供了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时，实现如上述的动态水体产品融合方法的各个步骤。

本发明实施例4提供的一种计算机可读存储介质具有上述实施例1中一种动态水体产品融合方法的有益效果。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims

1.一种动态水体产品融合方法，其特征在于，包括以下步骤：

S1、获取多个动态MODIS产品数据集；

2.根据权利要求1所述的一种动态水体产品融合方法，其特征在于，步骤S4具体包括：

3.根据权利要求1所述的一种动态水体产品融合方法，其特征在于，步骤S5具体包括：

4.根据权利要求3所述的一种动态水体产品融合方法，其特征在于，步骤S52具体包括：

5.根据权利要求4所述的一种动态水体产品融合方法，其特征在于，步骤S54具体包括以下分步骤：

若该像元值处于第二值域区间，则执行步骤S542；

6.根据权利要求4所述的一种动态水体产品融合方法，其特征在于，步骤S55具体包括以下分步骤：

若该像元值处于第一值域区间，则执行步骤S552；

7.根据权利要求5或6所述的一种动态水体产品融合方法，其特征在于，像元的梯度值的计算方法为：

8.一种动态水体产品融合装置，其特征在于，包括：

数据获取模块，用于获取多个动态MODIS产品数据集；

9.一种动态水体产品融合设备，其特征在于，包括：存储器和处理器；

所述存储器，用于存储程序；

所述处理器，用于执行所述程序，实现如权利要求1-7中任一项所述的动态水体产品融合方法的各个步骤。

10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时，实现如权利要求1-7中任一项所述的动态水体产品融合方法的各个步骤。