CN116245757B - 多模态数据的多场景通用性遥感影像云修复方法和系统 - Google Patents

多模态数据的多场景通用性遥感影像云修复方法和系统 Download PDF

Info

Publication number
CN116245757B
CN116245757B CN202310130461.XA CN202310130461A CN116245757B CN 116245757 B CN116245757 B CN 116245757B CN 202310130461 A CN202310130461 A CN 202310130461A CN 116245757 B CN116245757 B CN 116245757B
Authority
CN
China
Prior art keywords
pixel point
spectrum
image
target
reference image
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310130461.XA
Other languages
English (en)
Other versions
CN116245757A (zh
Inventor
梁治华
丁志平
魏海硕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Aisi Times Technology Co ltd
Original Assignee
Beijing Aisi Times Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Aisi Times Technology Co ltd filed Critical Beijing Aisi Times Technology Co ltd
Priority to CN202310130461.XA priority Critical patent/CN116245757B/zh
Publication of CN116245757A publication Critical patent/CN116245757A/zh
Application granted granted Critical
Publication of CN116245757B publication Critical patent/CN116245757B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/77—Retouching; Inpainting; Scratch removal
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/194—Segmentation; Edge detection involving foreground-background segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20081—Training; Learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20084—Artificial neural networks [ANN]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Processing (AREA)

Abstract

本申请涉及电数字数据处理技术领域，提供了一种多模态数据的多场景通用性遥感影像云修复方法和系统。该方法通过获取目标区域的目标影像和来自不同卫星的多模态遥感影像，基于多模态遥感影像的不同成像时间，分别确定最佳光学参考影像、最佳SAR参考影像，然后通过引入土地覆盖数据确定相似像素点集合(高质量相似像素点)，并通过对高质量相似像素点的合理权重分配与加权计算重建了缺失信息，实现对目标影像中云遮挡区域的修复，得到目标影像对应的云修复影像。通过该方法进行云修复，无论土地覆盖是否突变、发生何种程度的变化、土地覆盖突变是否分异等多种不同场景下，均能够取得更高精度、更佳效果的云修复结果，可以适用于不同的具体应用场景。

Description

多模态数据的多场景通用性遥感影像云修复方法和系统
技术领域
本申请涉及电数字数据处理技术领域，特别涉及一种多模态数据的多场景通用性遥感影像云修复方法和系统。
背景技术
光学影像在地球观测中发挥着重要作用，但由于云污染导致光学影像中大量有效数据缺失，降低了光学遥感影像的可用性，因此，去除光学影像中的云(也称云修复)一直是遥感领域的研究重点。
现有技术中，遥感影像云修复方法大致可分为：即基于空间的方法、基于多光谱的方法、基于时间的方法、混合方法以及基于深度学习的方法。然而，上述遥感影像云修复方法面临着应用场景单一、预测结果离散、所需辅助数据量大、计算成本高等问题。
因此，需要提供一种针对上述现有技术不足的改进技术方案。
发明内容
本申请的目的在于提供一种多模态数据的多场景通用性遥感影像云修复方法和系统，能够利用多个遥感卫星(比如Landsat-8、Landsat-9、Sentinel-1和Sentinel-2等)中不同卫星的成像时间与重访周期差异，确定最佳光学、SAR参考影像，进而通过权重分配和加权计算对云遮挡区域进行信息重建，以解决或缓解上述现有技术中存在的问题。
为了实现上述目的，本申请提供如下技术方案：
第一方面，本申请提供了一种多模态数据的多场景通用性遥感影像云修复方法，包括：
获取目标区域的多模态遥感影像和目标影像；所述多模态遥感影像包括来自不同卫星的光学遥感影像和来自不同卫星的合成孔径雷达SAR影像；所述目标影像为包含云遮挡区域的遥感影像；
基于所述多模态遥感影像的不同成像时间，分别从所述光学遥感影像和SAR影像中确定所述目标影像的最佳光学参考影像以及所述目标影像的最佳SAR参考影像；
基于所述目标影像、所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据，确定相似像素点集合；
依次对所述相似像素点集合中的各个像素点进行权重分配和加权计算，得到所述目标影像对应的云修复影像。
本申请中，通过获取来自不同卫星的光学遥感影像(如Landsat-8、Landsat-9)和来自不同卫星的合成孔径雷达SAR影像(如Sentinel-1、Sentinel-2)，组建虚拟星座，然后根据虚拟星座中不同成像时间和重访周期，分别确定目标影像的最佳光学参考影像和最佳SAR参考影像；通过引入土地覆盖数据以确定相似像素点集合，随后对相似像素点集合中的每个像素点进行权重分配和加权计算，以对目标影像中每个被云遮挡的目标像素点进行信息重建，最终得到目标影像的云修复影像。由于采用不同卫星、不同成像时间的光学影像与SAR影像作为参考影像，并引入土地覆盖数据对相似像素点进行筛选，在此基础上对相似像素点进行合理的权重分配和加权计算，使得相似像素点能够更综合地表达出目标像素点处地物景观的自然连贯性和演变情况，进而在多种地物变化场景下均能够有效的捕捉土地覆盖变化信息，并取得精度更高、效果更佳的云修复结果。
在一种可能的实现方式中，所述基于所述多模态遥感影像的不同成像时间，分别从所述光学遥感影像和SAR影像中确定所述目标影像的最佳光学参考影像以及所述目标影像的最佳SAR参考影像，具体为：
分别从所述光学遥感影像、SAR影像中选择成像时间在所述目标影像的成像时间之前，且与所述目标影像的成像时间间隔最短的无云影像；
将从所述光学遥感影像中选出的无云影像作为所述目标影像的最佳光学参考影像；将从所述SAR影像中选出的无云影像作为所述目标影像的最佳SAR参考影像。
本实施例中，由于多个卫星组成的虚拟星座中各卫星的时间分辨率不同，导致不同卫星的成像时间与重访周期也存在差异，当虚拟星座中任意一景目标影像中有云遮挡区域时，总能分别从不同来源的光学影像和不同来源的SAR影像中选出一景成像时间在目标影像的成像时间之前，且与目标影像的成像时间间隔最短的无云光学影像和无云SAR影像，作为对应的最佳光学参考影像和最佳SAR参考影像，在一定程度上降低了土地覆盖类型突变的可能性及其所带来的影响。
在一种可能的实现方式中，所述基于所述目标影像、所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据，确定相似像素点集合，具体为：
基于所述目标影像中云遮挡区域的中心点到云遮挡区域的边界之间的最短距离，构建相似像素的搜索区域；
根据所述搜索区域，对所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据进行求交集运算，以确定相似像素点集合。
在一种可能的实现方式中，所述相似像素点集合包括高质量相似像素点集合；
对应地，根据所述搜索区域，对所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据进行求交集运算，以确定相似像素点集合，具体为：
分别计算目标像素点与所述最佳光学参考影像在所述搜索区域内的所有像素点、所述目标像素点与所述最佳SAR参考影像在所述搜索区域内的所有像素点之间的光谱相似性，对应得到第一光谱相似性和第二光谱相似性；其中，所述目标像素点为所述目标影像中云遮挡区域内任一像素点；
根据预设相似性条件以及所述第一光谱相似性、所述第二光谱相似性，分别对所述最佳光学参考影像在所述搜索区域内的所有像素点、所述最佳SAR参考影像在所述搜索区域内的所有像素点进行分别筛，对应得到第一备选相似像素点集合、第二备选相似像素点集合；
在所述搜索区域内从所述土地覆盖数据中搜索与所述目标像素点具有相同土地覆盖类型的像素点，以得到第三备选相似像素点集合；
将所述第一备选相似像素点集合、所述第二备选相似像素点集合与所述第三备选相似像素点集合进行求交集运算，得到所述高质量相似像素点集合。
在一种可能的实现方式中，所述相似像素点集合还包括一般相似像素点集合；
当所述第一备选相似像素点集合中的像素个数为0或所述第三备选相似像素点集合中的像素个数为0时，将所述第二备选相似像素点集合作为一般相似像素点集合。
在一种可能的实现方式中，所述依次对所述相似像素点集合中的各个像素点进行权重分配和加权计算，得到所述目标影像对应的云修复影像，具体为：
基于所述相似像素点集合的光谱空间变化信息与空间位置，对所述目标像素点的光谱值进行预测，得到第一光谱预测值；
基于所述相似像素点集合的光谱时间变化信息，对所述目标像素点的光谱值进行预测，得到第二光谱预测值；
分别计算所述第一光谱预测值的权重、所述第二光谱预测值的权重，以确定所述目标像素点最终的光谱预测值，从而得到所述目标影像对应的云修复影像。
在一种可能的实现方式中，所述基于所述相似像素点集合的光谱空间变化信息与空间位置，对所述目标像素点的光谱值进行预测，得到第一光谱预测值，具体为：
根据所述最佳光学参考影像和所述最佳SAR参考影像分别在所述相似像素点集合的各像素点处的光谱值，以及所述最佳光学参考影像和所述最佳SAR参考影像分别在目标像素点位置处的光谱值，计算得到第一光谱差值；
根据所述相似像素点集合中各像素点的位置、所述目标像素点的位置，计算得到空间距离；
分别对所述第一光谱差值、所述空间距离进行归一化处理，并计算得到所述相似像素点集合中各像素点的第一权重；
根据所述目标影像在所述相似像素点集合的各像素点处的光谱值以及所述相似像素点集合中各像素点的第一权重，对所述目标像素点的光谱值进行预测，得到第一光谱预测值。
在一种可能的实现方式中，所述基于所述相似像素点集合的光谱时间变化信息，对所述目标像素点的光谱值进行预测，得到第二光谱预测值，具体为：
计算t1时刻在最佳光学参考影像中相似像素点集合中各像素点与目标像素点的光谱值之差，以得到第二光谱差值；其中，所述t1时刻为所述最佳光学参考影像的成像时间；
计算t2时刻相似像素点集合中各像素点与目标像素点在最佳SAR参考影像中的光谱值之差，以得到第三光谱差值；其中，所述t2时刻为所述最佳SAR参考影像的成像时间；
分别基于所述第二光谱差值、所述第三光谱差值进行计算，对应得到所述相似像素点集合中各像素点的第二权重、第三权重；
根据所述第二权重、所述第三权重，以及所述相似像素点集合在所述t1时刻、所述t2时刻、所述t3时刻的光谱值，对所述目标像素点的光谱值进行预测，得到第二光谱预测值。
在一种可能的实现方式中，所述分别计算所述第一光谱预测值的权重、所述第二光谱预测值的权重，以确定所述目标像素点最终的光谱预测值，从而得到所述目标影像对应的云修复影像，具体为：
根据所述相似像素点集合与所述目标像素点的平均光谱差值Rss，计算所述第一光谱预测值的权重；
根据所述相似像素点集合的综合平均光谱差值Rst，计算所述第二光谱预测值的权重；
根据所述第一光谱预测值、所述第一光谱预测值的权重、所述第二光谱预测值、所述第二光谱预测值的权重，计算所述目标像素点最终的光谱预测值，从而得到所述目标影像对应的云修复影像。
第二方面，本申请提供一种多模态数据的多场景通用性遥感影像云修复系统，包括：
获取单元，配置为获取目标区域的多模态遥感影像和目标影像；所述多模态遥感影像包括来自至少一个卫星的光学遥感影像和来自至少一个卫星的合成孔径雷达SAR影像；所述目标影像为包含云遮挡区域的遥感影像；
确定单元，配置为基于所述多模态遥感影像的不同成像时间，分别从所述光学遥感影像和SAR影像中确定所述目标影像的最佳光学参考影像以及所述目标影像的最佳SAR参考影像；
筛选单元，配置为基于所述目标影像、所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据，确定相似像素点集合；
修复单元，配置为依次对所述相似像素点集合中的各个像素点进行权重分配和加权计算，得到所述目标影像对应的云修复影像。
有益效果：
本申请的技术方案中，通过获取目标区域的多模态遥感影像和目标影像；多模态遥感影像包括来自不同卫星的光学遥感影像和来自不同卫星的合成孔径雷达SAR影像；目标影像为包含云遮挡区域的遥感影像；然后基于多模态遥感影像的不同成像时间，分别从光学遥感影像和SAR影像中确定目标影像的最佳光学参考影像以及目标影像的最佳SAR参考影像；接着基于目标影像、最佳光学参考影像、最佳SAR参考影像以及预先获取的目标区域的土地覆盖数据，确定相似像素点集合；依次对相似像素点集合中的各个像素点进行权重分配和加权计算，得到目标影像对应的云修复影像。由于采用不同卫星、不同成像时间的光学影像与SAR影像作为参考影像，并引入土地覆盖数据对相似像素点进行筛选，在此基础上对相似像素点进行权重分配和加权计算，使得相似像素点能够更综合地表达出目标像素点处地物景观的自然连贯性和演变情况。实验结果表明，无论目标区域的土地覆盖是否突变、发生何等程度的变化、土地覆盖突变是否分异等多种不同场景下，通过本申请的技术方案均可取得精度更高、效果更好的云修复效果。
由于SAR影像与光学影像相比具有普遍的离散特性，本申请的技术方案中，采用对相似像素点进行合理分配权重和加权计算，使得SAR影像作为参考影像进行云修复的结果并未继承SAR影像的离散特性，使修复后的目标像素点能够与周边的无云光学影像具有较好的自然过渡性，进而提高云修复的精度，同时也降低了对计算机算力的需求。
此外，采用本申请的技术方案，可以选择目前可广泛获取的、免费的多个卫星作为数据来源，组建虚拟星座以提供多模态数据，然后从多模态数据中选择极少量的遥感影像作为参考影像，并能够在现有遥感卫星数据条件下尽最大可能地对地物景观(土地覆盖)的变化进行有效捕捉，进一步提高云修复的精度，而且由于仅需要从现有不同来源的多模态数据中选择最佳光学参考影像和最佳SAR影像，使得本方法的实施不再受到需要大量高质量数据的限制，提高了其适用性。
附图说明
构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。其中：
图1为根据本申请的一些实施例提供的多模态数据的多场景通用性遥感影像云修复方法的流程示意图；
图2为根据本申请的一些实施例提供的相似像素点筛选的逻辑示意图；
图3为根据本申请的一些实施例提供的多模态数据的多场景通用性遥感影像云修复系统的结构示意图。
具体实施方式
为了便于理解本申请的技术方案，下面对相关术语进行说明。
遥感(Remote Sensing)是指非接触的、远距离的探测技术，一般指运用传感器/遥感器对物体的电磁波的辐射、反射特性的探测，获取其反射、辐射或散射的电磁波信息(如电场、磁场、电磁波、地震波等信息)，并进行提取、判定、加工处理、分析与应用的一门科学和技术。
人造卫星(简称卫星)是主要的遥感平台，是遥感过程中乘载传感器/遥感器的运载工具，它如同在地面摄影时安放照相机的三脚架，是在空中或空间安放传感器/遥感器的装置。
传感器/遥感器是远距离感测地物环境辐射或反射电磁波的仪器，传感器的种类有多种，除常见的可见光摄影机、红外摄影机、紫外摄影机外，还有红外扫描仪、多光谱扫描仪、微波辐射和散射计、侧视雷达、专题成像仪、成像光谱仪等，传感器正在向多光谱、多极化、微型化和高分辨率的方向发展。
陆地卫星(Landsat)是美国NASA发射的一系列太阳同步近极地圆形轨道卫星的总称，总共包括8颗卫星，其上搭载的传感器包括：多光谱扫描仪(Multispectral Scanner，MSS)、专题制图仪(Thematic Mapper，TM)、增强专题成像仪(Enhanced Thematic MapperPlus，ETM+)、陆地成像仪(Operational Land Imager，OLI)和热红外传感器(ThermalInfrared Sensor，TIRS)2个传感器等，能够获取大量免费的遥感影像，其空间分辨率为30～100米。自2021年9月Landsat-9成功发射，Landsat-8与其协同工作后，Landsat影像的时间分辨率降至8天。
哨兵系列卫星(Sentinel-1～Sentinel-6)是欧洲太空总署哥白尼计划的一部分，主要目标便是对地球的观测任务。其中，Sentinel-2在A、B两颗卫星的协同工作下，其遥感影像的时间分辨率达到5天，Sentinel-1同样在A、B两颗卫星的协同工作下，其合成孔径雷达(synthetic aperture radar，SAR)数据的时间分辨率为6天。同时，哨兵系列卫星携带的多光谱成像仪(Multispectral Instrument，MSI)可覆盖13个光谱波段，能够提供多种空间分辨率的光学遥感影像。
因遥感技术能够周期性地捕获地表、海洋、大气等信息，使得其在环境监测、资源调査、灾害应急等方面发挥着重要作用。然而，在任何时间点获取的遥感影像中，总有相当一部分地球表面被云层覆盖，云层的覆盖严重阻碍了光学频带电磁波的有效传输，导致光学遥感影像中的有效数据缺失，大大降低了光学遥感影像的可用性。因此，如何能够根据实际应用需求，有效重建云覆盖区域地表信息，获得无云的高质量遥感影像，始终是遥感领域急需解决的重要问题。过去的几十年以来，已经陆续的提出了许多光学影像云修复方法。总体来看，传统遥感影像云修复方法大致可分为：即基于空间的方法、基于多光谱的方法、基于时间的方法、混合方法以及基于深度学习的方法，下面对现有的云修复方法进行详细介绍。
基于空间的方法主要是根据遥感影像局部或非局部区域的空间自相关性，利用无云覆盖区域的信息来填充云覆盖区域，包括最近邻插值、克里金插值、传播扩散法、基于样本的方法等。基于空间的方法的重要前提假设是被云遮挡的区域所缺失数据和其他无云区域的数据具有相同的统计或几何结构，多适用于面积较小的云遮挡区域，当云遮挡区域逐渐增大，地物类型逐渐复杂时，由于基于空间的方法缺少可用于判定修复保真度的参考影像，往往会导致云修复结果过渡平滑，失真现象较为严重。
基于多光谱的方法主要是根据不同波段的遥感影像具有不同穿云透雾能力，通过模拟穿透性较强波段所获得的完整地物信息与穿透性较弱波段所获得的不完整地物信息的潜在相关关系，利用穿透性较强波段的完整地物信息来重建穿透性较弱波段的不完整地物信息，但基于多光谱的方法大多适用于薄云或雾、霾的修复重建，当云层变厚时，几乎所有的光学波段都无法有效穿透，也就导致几乎所有的光学波段遭受信息缺失。
基于时间的方法主要是利用同一区域不同时间的影像来获得云遮盖区域的互补信息，与上述两种方法相比，其往往表现出更好的云修复效果。基于时间的方法可进一步分类为时间替换和时间滤波，其中时间替换法又可进一步细分为直接替换法和间接替换法。直接替换法是将相同区域不同时间内的无云影像直接替换目标影像的云覆盖区，当时间序列较长，可选的无云影像较多时，可在其中选择最优值进行替换，但在实际情况下，往往有云覆盖的目标影像与相同区域的无云影像的时间间隔较长，且由于探测器角度、太阳角度、大气影响等因素，不同场景下相同位置的影像之间存在差异，进而也就导致大多数情况下直接替换法效果较差。间接替换法是通过一些列的变换、校正，尽可能的消除因参考影像与要修复的目标影像的时间差异所带来的像素值差异，如直方图匹配法、邻域相似像素插值法(NSPI)、多时间加权线性回归等。间接替换法与直接替换法相比，云修复效果明显提升，但受卫星重访周期与天气因素影响，要进行云修复的目标影像往往与无云参考影像的成像时间间隔较长，如果在此时间间隔内土地覆盖类型发生突变，进而会破坏像素之间的相关性，明显降低云修复效果。时间滤波法多使用相同区域的长时间序列的遥感影像对目标影像进行云修复，其通过采集长时间序列内同一像素在所有时间点的像素值并对其进行拟合处理。该方法将云污染区域看作时间序列中的噪声，并根据拟合的波动变化规律对噪声进行修正，如滑动窗口滤波器，基于函数的曲线拟合，以及频域方法。但时间滤波法侧重于反映地表地物的自然变化规律，若地表土地覆盖类型发生突变，该方法将无法有效捕捉变化信息，进而导致云修复结果存在明显误差。
混合方法通过综合空间、光谱、时间中的两种或三种类型的相关性，以期综合利用上述三种方法的优势，进而有效重建云遮挡区域信息，例如，基于光谱角度映射器的时空相似性(SAMSTS)算法，数据插值经验正交函数(DINEOF)方法，时空MRF模型指导下的相似像素替换法(STMRF)。混合方法的基本思想是：选择一定数量的无云影像作为参考影像，当无云参考影像的数量较多时，往往会取得较好的云修复效果，然而，要进行云修复的目标影像往往与作为参考的无云光学影像的成像时间间隔较长，在此时间间隔内，若土地覆盖类型发生变化，基于混合方法所捕捉到的各像素之间的相关关系将不再有效，也就无法得到合理有效的重建信息。
近十年来，随着深度学习的不断发展，一些基于深度的学习的多时相遥感影像云修复方法不断被提出，比如，利用统一时空谱深度卷积神经网络(STS-CNN)来重建多时相的Landsat数据；又比如通过建立Cloud-GAN网络学习到了多云图像和无云图像之间的映射，实现了多时相Sentinel-2影像中的薄云修复；还有些学者利用深度卷积生成对抗网络(DCGAN)修复了海表温度图像中因云遮挡而缺失的部分。然而基于深度学习的方法之所以具有强大的学习与数据生成能力的前提条件是，大量可用于训练的高质量多时相光学遥感数据，当输入数据质量较低(如大面积云覆盖)或输入数据数量较少时，极易面临训练和预测的不稳定性挑战。
基于上述分析可知，虽然现有的遥感影像云修复方法有多种，比如，基于空间的方法、基于多光谱的方法、基于时间的方法、混合方法以及基于深度学习的方法，但是，基于空间的方法和基于多光谱的方法因其本身的局限性，在实际应用中云修复效果相对较差，基于时间的多时相云修复方法、混合方法以及基于深度学习的云修复方法正逐渐成为主流，但基于时间的和混合方法，因其可获得参考影像往往与目标影像的时间间隔相对较长，往往忽略了在此期间内的土地覆盖类型变化情况，使其仅能适用于地表土地覆盖类型无变化或符合一般自然变化规律的场景，无法有效应对土地覆盖类型发生突变的场景，而基于深度学习的方法往往需要大量的训练数据，模型较为复杂，计算量与计算成本相对较大，还普遍面临着预测的不稳定性。
为此，本申请综合利用广泛可获取的卫星影像，如landsat系列(比如Landsat-8、Landsat-9)、哨兵卫星(Sentinel-1、Sentinel-2)等光学影像，以及不同卫星的不同成像时间、重访周期差异及穿云透雾等优势，提供了一种多模态数据的多场景通用性遥感影像云修复方法(Multi-scene universal cloud removal based on multi-modal data，简称MUCRMD)，该技术方案通过从目前可普遍获得的多模态遥感影像中选择极少量高质量的无云影像为参考影像，且参考影像包括光学和SAR两种不同的成像方式的影像，利用合理的权重分配和加权计算将光学参考影像和SAR影像相结合，在现有遥感卫星数据条件下尽最大可能的有效捕捉到了地物变化信息，形成了无论云遮挡区域地物类型是否变化、发生何种程度的变化，均可有效适用的多场景通用性概念框架，在一定程度上保证了预测的稳定性，降低了SAR数据的离散性继承，提高了云遮挡区域信息重建效率。
下面将参考附图并结合实施例来详细说明本申请。
各个示例通过本申请的解释的方式提供而非限制本申请。实际上，本领域的技术人员将清楚，在不脱离本申请的范围或精神的情况下，可在本申请中进行修改和变型。例如，示为或描述为一个实施例的一部分的特征可用于另一个实施例，以产生又一个实施例。因此，所期望的是，本申请包含归入所附权利要求及其等同物的范围内的此类修改和变型。
示例性方法
本申请实施例提供一种多模态数据的多场景通用性遥感影像云修复方法，如图1、图2所示，该方法包括：
步骤S101、获取目标区域的多模态遥感影像和目标影像；多模态遥感影像包括来自不同卫星的光学遥感影像和来自不同卫星的合成孔径雷达SAR影像；目标影像为包含云遮挡区域的遥感影像。
其中，目标影像是任意一幅包含云遮挡区域的光学遥感影像。需要说明的是，根据成像原理不同，遥感影像包括光学和微波两种成像方式，其中，光学成像与微波成像主要区别在于其接受电磁波的范围，光学传感器是接收可见光到红外区，微波传感器主要是微波区。光学成像卫星通常是指对目标在可见光、近红外和短波红外电磁谱段进行成像观测，获取和分析被观测对象的光学特征信息的卫星。光学成像卫星占遥感卫星总量的52.7％，是最为主要的遥感卫星。运用雷达成像被称为微波遥感卫星，典型的如合成孔径雷达(SAR)卫星。与光学遥感器相比，具有全天候、全天时及能穿透云层的成像特点。
多模态遥感影像指的是不同数据来源的遥感影像，特别地，可以是来自不同卫星所搭载的传感器获取的遥感影像，也可以是来自于同一个卫星所搭载的不同传感器所获得的遥感影像。本申请实施例中，将多个卫星(比如Landsat-8、Landsat-9、Sentinel-1和Sentinel-2卫星)组建虚拟星座，通过综合利用不同遥感卫星的成像时间与重访周期差异，以及光学、SAR影像的各自成像优势，能够有效实现中分辨率光学影像云去除。
步骤S102、基于多模态遥感影像的不同成像时间，分别从光学遥感影像和SAR影像中确定目标影像的最佳光学参考影像以及目标影像的最佳SAR参考影像。
卫星成像时间指的是所搭载的遥感器/传感器对地物进行采样的时间点。不同卫星的成像时间与该卫星的重访周期(也称时间分辨率、卫星重复周期)有关。重复周期/轨道重访周期指的是卫星相对地球而言，两次轨道完全重合所需时间间隔，即人造卫星第一次经过某个星下点(地球中心与人造卫星的连线在地球表面上的交点)和第二次经过同一个星下点的时间间隔。比如Landsat系列卫星的时间分辨率为8天，Sentinel系列卫星的时间分辨率为5天。重访周期越短，卫星会以更高的时间分辨率对地面覆盖区域进行采用，对较短时间内有可能发送变化的地物景观观测越有利。
本申请实施例中，为了克服单个卫星时间分辨率的限制，从多个卫星组成的虚拟星座中选择分别选择目标影像的最佳光学参考影像以及目标影像的最佳SAR参考影像，充分利用不同卫星重访时间的差异，使得参考影像与目标影像的时间间隔尽可能达到最短，以最小化地物景观变化对云修复带来的影响。
步骤S103、基于目标影像、最佳光学参考影像、最佳SAR参考影像以及预先获取的目标区域的土地覆盖数据，确定相似像素点集合。
土地覆盖数据是表征地球表面当前所具有的自然和人为影响所形成的覆盖物的数据，根据地球表面覆盖物的不同，土地覆盖可以有多种类型，如水体、森林、草场、农田、土壤、冰川、湖泊、沼泽湿地及道路等。根据数据组织的方式，土地覆盖数据可以是矢量数据，也可以是栅格数据。在矢量格式的土地覆盖数据中，每一个矢量图斑对应一种土地覆盖类型，为了便于运算，矢量格式的土地覆盖数据在与遥感影像进行运算之前，可以先进行栅格化，使其变为栅格格式的土地覆盖数据(简称栅格数据)。所得到的栅格数据中，像元的光谱值表示对应的土地覆盖类型。土地覆盖数据为栅格数据时，因遥感影像(包括光学遥感影像和SAR影像)均为由多个像素点组成的栅格数据，更方便后续的计算和处理。其中，栅格数据的每个像素点表示地球表面一定的实际面积大小，这种对应关系称为空间分辨率。
由于目标影像的成像时间、最佳光学参考影像的成像时间以及最佳SAR参考影像的成像时间均不同，本申请实施例中，通过引入土地覆盖数据以确定相似像素点集合，从而能够确定参考影像(包括最佳光学参考影像和最佳SAR参考影像)与目标影像在相似像素点集合中各相似像素点处的地物景观的变化情况，为后续对目标像素点进行高精度预测奠定了基础。
步骤S104、依次对相似像素点集合中的各个像素点进行权重分配和加权计算，得到目标影像对应的云修复影像。
本申请实施例中，根据相似像素点集合中的各个像素点在参考影像中的光谱值，通过合理的权重分配，并进行加权计算，巧妙地建立参考影像与目标影像之间的对应关系，最终得到目标影像对应的云修复影像。该方法不仅避免了现有的基于深度学习模型和SAR影像进行云修复带来的离散性继承的问题，而且大大降低了云修复所需的计算量，提高了云修复的效率。
本申请实施例中，通过获取来自不同卫星的光学遥感影像和来自不同卫星的合成孔径雷达SAR影像，进而组建虚拟星座，然后根据虚拟星座中不同成像时间和重访周期，分别确定目标影像的最佳光学参考影像和最佳SAR参考影像；然后通过引入土地覆盖数据确定相似像素点集合，随后对相似像素点集合中的每个像素点进行权重分配和加权计算，以对目标影像中每个被云遮挡的目标像素点进行信息重建，最终得到目标影像的云修复影像。由于采用不同卫星、不同成像时间的光学影像与SAR影像作为参考影像，并引入土地覆盖数据对相似像素点进行筛选，使得相似像素点能够更综合地表达出目标像素点处地物景观的自然连贯性和演变情况，进而使得多种地物变化场景下均能得到均可更有效的捕捉土地覆盖变化信息，然后进行合理的权重分配和加权计算，取得精度更高、效果更佳的云修复结果。
在本申请的一些实施例中，基于多模态遥感影像的不同成像时间，分别从光学遥感影像和SAR影像中确定目标影像的最佳光学参考影像以及目标影像的最佳SAR参考影像，具体为：分别从光学遥感影像、SAR影像中选择成像时间在目标影像的成像时间之前，且与目标影像的成像时间间隔最短的无云影像；将从光学遥感影像中选出的无云影像作为目标影像的最佳光学参考影像；将从SAR影像中选出的无云影像作为目标影像的最佳SAR参考影像。
为减少地物变化对云修复的影响，本申请实施例中，在确定需要修复的目标影像之后，基于目标影像的成像时间，从虚拟星座提供的多模态数据中分别确定目标影像的最佳光学参考影像和最佳SAR参考影像。
具体来说，目标影像确定后，由于虚拟星座包括多种不同来源的光学影像，可以从多种不同来源的光学影像确定一景成像时间在目标影像的成像时间之前，且与目标影像的成像时间间隔最短的无云光学影像，作为最佳光学参考影像，由此，在一定程度上极大提高了无云光学影像的选择空间，可以尽可能地选择距离目标影像最近的无云光学影像作为最佳光学参考影像，在一定程度上降低了土地覆盖类型突变的可能性及其对云修复所带来的误差影响。
同时，由于SAR影像具有穿透云雾的特性，故所有的SAR影像均无云遮挡，因此只需从不同来源的SAR影像中确定一景成像时间在目标影像的成像时间之前，且与目标影像的成像时间间隔最短的SAR影像，即可作为最佳SAR参考影像。
比如由Landsat-8、Landsat-9、Sentinel-2A和Sentinel-2B卫星组建的虚拟星座，可以提供4中不同来源的光学影像以及SAR影像，这样，对于确定的目标区域，理论上每月约有10景光学影像和5景SAR影像可以作为备选参考影像，且从成像时间来说，SAR影像穿插分布在光学影像之间，因此，目标影像的成像时间与SAR影像的成像时间之间最大时间间隔为5天，同时也存在时间间隔为4、3、2、1、0天等情景，进一步缩短了参考影像与目标影像的成像时间间隔，保证了地物景观之间的连续性。
此外，由于SAR影像具有穿透云雾的成像特性，基于SAR影像可以有效判断无云光学影像与目标影像的成像时间间隔内云遮挡区域的土地覆盖类型是否发生变化，并且能够尽可能地捕捉这些土地覆盖变化信息，避免因土地覆盖突变而导致一系列云修复的误差。
在一些实施例中，基于目标影像、最佳光学参考影像、最佳SAR参考影像以及预先获取的目标区域的土地覆盖数据，确定相似像素点集合，具体为：基于目标影像中云遮挡区域的中心点到云遮挡区域的边界之间的最短距离，构建相似像素的搜索区域；根据搜索区域，对最佳光学参考影像、最佳SAR参考影像以及预先获取的目标区域的土地覆盖数据进行求交集运算，以确定相似像素点集合。
本申请实施例中，将目标影像，成像时间在目标影像之前的最佳光学参考影像、最佳SAR参考影像，以及预先获取的目标区域的土地覆盖数据作为输入数据，首先根据目标影像中云遮挡区域的范围确定相似像素的搜索区域，然后根据该搜索区域以及一系列的相似性原则筛选出目标影像中云遮挡区域中个像素对应的相似像素点集合，为后续根据相似像素点集合预测云遮挡区域中每个像素点的光谱值奠定基础。
需要说明的是，由于地表地物的连贯性与自然过渡性，可以确定在云遮挡区域的相邻区域应当存在与云遮挡区域的光谱特征相似的像素点。基于此，本申请实施例中，通过设置相连区域的大小即可确定相似像素的搜索区域，然后根据搜索区域，对最佳光学参考影像、最佳SAR参考影像以及预先获取的目标区域的土地覆盖数据进行求交集运算，以确定相似像素点集合。
具体地，参见图2，本申请实施例的相似像素点集合由如下步骤确定：首先，计算云遮挡区域的中心点到云遮挡区域的边界之间的最短距离d；然后取该最短距离d的一半(即d/2)作为云遮挡区域向外扩张的缓冲半径，通过缓冲得到云遮挡区域的缓冲区，该缓冲区即为云遮挡区域的邻域，即相似像素的搜索区域，然后以光谱相似性为基础构建用于筛选相似像素点的相似性原则，并基于最佳光学参考影像、最佳SAR参考影像、土地覆盖数据进行求交集运算，从而得到最终的相似像素点集合。
在一些实施例中，相似像素点集合包括高质量相似像素点集合；对应地，根据搜索区域，对最佳光学参考影像、最佳SAR参考影像以及预先获取的目标区域的土地覆盖数据进行求交集运算，以确定相似像素点集合，具体为：分别计算目标像素点与最佳光学参考影像在搜索区域内的所有像素点、目标像素点与最佳SAR参考影像在搜索区域内的所有像素点之间的光谱相似性，对应得到第一光谱相似性和第二光谱相似性；根据预设相似性条件以及第一光谱相似性、第二光谱相似性，分别对最佳光学参考影像在搜索区域内的所有像素点、最佳SAR参考影像在搜索区域内的所有像素点进行分别筛，对应得到第一备选相似像素点集合、第二备选相似像素点集合；在搜索区域内从土地覆盖数据中搜索与目标像素点具有相同土地覆盖类型的像素点，以得到第三备选相似像素点集合；将第一备选相似像素点集合、第二备选相似像素点集合与第三备选相似像素点集合进行求交集运算，得到高质量相似像素点集合。
其中，目标像素点为目标影像中云遮挡区域内任一像素点。
本申请实施例中，参考影像在搜索区域内任一像素点与目标像素点之间的光谱相似性称为第一光谱相似性，可以采用平均绝对误差、均方误差或者均方根误差来表征。
本申请实施例采用均方根误差RMSD表示第一光谱相似性，其公式如下：
式中：RMSDoi表示第一光谱相似性；最佳光学参考影像的波段总数为n_oi个，n_oi个波段中的第b_oi个波段用b_oi表示；(x，y)表示目标像素点的坐标，(xi，yi)表示搜索区域内的第i个像素点的坐标，t1为最佳光学参考影像的成像时间，L(xi，yi，t1，b_oi)表示t1时刻的最佳光学参考影像中位于第b_oi个波段的(xi，yi)处的光谱值，L(x，y，t1，b_oi)表示t1时刻的最佳光学参考影像中位于b_oi波段的目标像素点(x，y)处的光谱值。
目标像素点与最佳SAR参考影像在搜索区域内的所有像素点之间的光谱相似性称为第二光谱相似性，本申请实施例中，采用均方根误差RMSD表示第二光谱相似性，其公式如下：
式中：RMSDsar表示第二光谱相似性；最佳SAR参考影像的波段数量为n_sar个，n_sar个波段中的第b_sar个波段用b_sar表示；(x，y)表示目标像素点的坐标，(xi，yi)表示搜索区域内的第i个像素点的坐标，t2为最佳SAR参考影像的成像时间，L(xi，yi，t2，b_sar)表示t2时刻的最佳SAR参考影像中位于b_sar波段的(xi，yi)处的光谱值，L(x，y，t2，b_sar)表示t2时刻的最佳SAR参考影像中位于b_sar波段的目标像素点(x，y)处的光谱值。
应理解，随着第一光谱相似性RMSDoi、第二光谱相似性RMSDsar数值的增大，对应的最佳光学参考影像、最佳SAR参考影像中的像素点与目标像素点的相似性也就越小。基于此，本申请实施例中，根据预设相似性条件以及第一光谱相似性、第二光谱相似性，分别对最佳光学参考影像在搜索区域内的所有像素点、最佳SAR参考影像在搜索区域内的所有像素点进行分别筛，对应得到第一备选相似像素点集合、第二备选相似像素点集合。
具体地，预设相似性条件可以是根据目标像素点的自身光谱值的百分比例来设置第一光谱相似性、第二光谱相似性可允许的光谱误差范围，进而根据每个目标像素点的光谱特性形成其独有的光谱相似性判断阈值。
其中，第一光谱相似性对应的相似性条件如下：
式中，abs(*)表示取绝对值，m_oi表示针对最佳光学参考影像所确定的目标像素点(x，y)处的光谱值的百分比。
根据公式(3)，即可从最佳光学参考影像在搜索区域内的所有像素点中筛选出第一光谱相似性满足公式(3)的像素点，即得到第一备选相似像素点集合。
第二光谱相似性对应的相似性条件如下：
式中，m_sar表示针对最佳SAR参考影像所确定的目标像素点(x，y)处的光谱值的百分比。
根据公式(4)，即可从最佳SAR参考影像在搜索区域内的所有像素点中筛选出第二光谱相似性满足公式(4)的像素点，即得到第二备选相似像素点集合。
这里，m_oi、m_sar可以是根据先验知识所确定的经验值，通常，可以将其设置为大于0％且小于100％的数值。优选地，当m_oi、m_sar的取值为2％～4％时，可以达到较佳的效果。当地球表面覆盖的地物类型比较复杂时，也可以适当提高m_oi、m_sar的取值。
在分别基于最佳光学参考影像、最佳SAR参考影像确定第一备选相似像素点集合、第二备选相似像素点集合之后，还需要根据土地覆盖数据，在搜索区域内确定与目标像素点具有相同土地覆盖类型的像素点组成的区域，以得到第三备选相似像素点集合。
最后，将从最佳光学参考影像中搜索到的第一备选相似像素点集合、从最佳SAR参考影像中搜索到的第二备选相似像素点集合以及从土地覆盖数据中搜索到的第三备选相似像素点集合进行求交集运算，得到高质量相似像素点集合，相应的表达式如下：
式中，land cover type(xi，yi)表示像素点(xi，yi)处的土地覆盖类型，landcover type(x，y)表示目标像素点处的土地覆盖类型。
实际应用中，在极少数情况下，第一备选相似像素点集合、第三备选相似像素点集合可能为空集，此时，为保证这种极少数情况下也能够对目标像素点顺利地进行云修复，在本申请的一些实施例中，相似像素点集合为一般相似像素点集合；当第一备选相似像素点集合中的像素个数为0或第三备选相似像素点集合中的像素个数为0时，将第二备选相似像素点集合作为一般相似像素点集合，也就是采用从最佳SAR参考影像中确定的备选相似像素点相似像素点集合，参与后续的运算。
此外，当基于公式(5)对目标像素点、最佳光学参考影像、最佳SAR参考影像、土地利用覆盖数据进行交集运算得到的结果是空集时，也就是说，在目标像素点处，从最佳光学参考影像、最佳SAR参考影像、土地利用覆盖数据没有找到高质量相似像素点，则将第二备选相似像素点集合作为一般相似像素点集合。
在一些实施例中，依次对相似像素点集合中的各个像素点进行权重分配和加权计算，得到目标影像对应的云修复影像，具体为：基于相似像素点集合的光谱空间变化信息与空间位置，对目标像素点的光谱值进行预测，得到第一光谱预测值；基于相似像素点集合的光谱时间变化信息，对目标像素点的光谱值进行预测，得到第二光谱预测值；分别计算第一光谱预测值的权重、第二光谱预测值的权重，以确定目标像素点最终的光谱预测值，从而得到目标影像对应的云修复影像。
在获取相似像素点集合后，本申请实施例中，根据目标像素点与相似像素点之间的光谱差异、空间距离差异，对相似像素点进行合理的权重分配，并根据相似像素点集合的光谱空间变化信息与空间位置完成基于光谱-空间信息的预测得到第一光谱预测值，根据相似像素点集合的光谱时间变化信息完成基于光谱-时间信息预测得到第二光谱预测值；最后对基于光谱-空间信息的预测和基于光谱-时间信息预测的结果进行综合加权平均，完成综合光谱-空间-时间信息的预测。
具体地，首先进行基于光谱-空间信息的预测。在一些实施例中，基于相似像素点集合的光谱空间变化信息与空间位置，对目标像素点的光谱值进行预测，得到第一光谱预测值，具体为：根据最佳光学参考影像和最佳SAR参考影像分别在相似像素点集合的各像素点处的光谱值，以及最佳光学参考影像和最佳SAR参考影像分别在目标像素点位置处的光谱值，计算得到第一光谱差值；根据相似像素点集合中各像素点的位置、目标像素点的位置，计算得到空间距离；分别对第一光谱差值、空间距离进行归一化处理，并计算得到相似像素点集合中各像素点的第一权重；根据目标影像在相似像素点集合的各像素点处的光谱值以及相似像素点集合中各像素点的第一权重，对目标像素点的光谱值进行预测，得到第一光谱预测值。
经过前述的相似像素点筛选步骤，使得相似像素点与目标像素点在光谱上具有较高的相似性，因此，在基于光谱-空间信息进行预测时可以直接利用得到的相似像素点在t3时刻的光谱信息对目标像素点的光谱值进行预测。
考虑到不同的相似像素点与目标像素点的光谱差值、空间距离不同，导致不同的相似像素点对目标像素点具有不同的贡献度。具体来说，从光谱差异来看，相似像素点与目标像素点的光谱差值越小，其在对目标像素点的预测过程中的贡献度就越大；从空间距离来看，距离目标像素点越近的相似像素点往往与目标像素点越相似，其在对目标像素点的预测过程中的贡献度就越大。此外，为了尽可能提高云遮挡区域与相邻区域的连贯性和自然过渡性，在计算不同相似像素点的权重时，综合考虑了光谱和距离的因素。
应当理解，在搜索得到高质量相似像素的情况下，优先使用高质量像素点对目标像素点进行预测，反之，在没有搜索到高质量像素点的情况下，使用一般相似像素点作为参与计算的相似像素点对目标像素点进行预测。
首先，根据最佳光学参考影像和最佳SAR参考影像分别在相似像素点集合的各像素点处的光谱值，以及最佳光学参考影像和最佳SAR参考影像分别在目标像素点位置处的光谱值，计算第一光谱差值。
本申请实施例中，第一光谱差值基于相似像素点与目标像素点在最佳光学参考影像中的光谱差值，以及，相似像素点与目标像素点在最佳SAR参考影像中的光谱差值进行计算，其中，第一光谱差值可以采用相似像素点与目标像素点相应类型的光谱值直接求差值得到，也可以通过计算相似像素点与目标像素点之间的RMSD，并用RMSD表示第一光谱差值。
采用相似像素点与目标像素点之间的RMSD表示第一光谱差值时，其计算公式如下：
式中，RMSDi表示相似像素点(xi，yi)与目标像素点(x，y)的第一光谱差值。
然后，根据相似像素点的空间位置、目标像素点的空间位置，计算空间距离，具体计算公式如下：
式中，Di表示相似像素点(xi，yi)与目标像素点(x，y)的空间距离。
当云遮挡区域的面积较大时，相似像素点与目标像素点之间的空间距离也比较大，有可能造成空间距离与第一光谱差值的数值不可比，即二者不在统一量纲标准下，因此，分别对第一光谱差值和空间距离进行归一化处理，以将第一光谱差值和空间距离缩放到统一量纲下。对第一光谱差值进行归一化处理的计算公式如下：
式中，表示归一化处理后的第一光谱差值，RMSDi_max、RMSDi_min分别表示第一光谱差值的最大值、最小值。
对空间距离进行归一化处理的计算公式如下：
式中，表示归一化处理后的空间距离，Di_max、Di_min分别空间距离的最大值、最小值。
接着，根据归一化处理后的第一光谱差值、归一化处理后的空间距离，计算相似像素点的第一权重，其将表达式如下：
式中，Wi表示相似像素点(xi，yi)的第一权重，N表示相似像素点的数量。
在搜索得到高质量相似像素的情况下，采用高质量相似像素点作为相似像素点带入公式(6)～(10)，即采用最佳光学参考影像与最佳SAR参考影像在相似像素点处的光谱值进行运算，若没有搜索到高质量像素点，只找到一般相似像素点，则将一般相似像素点作为相似像素点进行运行，此时只需将公式(6)中光学影像的光谱差值计算部分去除即可，其余公式保持不变，进而计算得到相似像素点的第一权重。
最后，根据目标影像在相似像素点集合的各像素点处的光谱值以及相似像素点集合中各像素点的第一权重，对其进行加权平均计算，即可得到基于光谱-空间预测的目标像素点的光谱预测值，也就是第一光谱预测值。其计算公式如下：
式中，Lss(x，y，t3，b_oi)表示第一光谱预测值，也就是t3时刻目标像素点(x，y)处的b_oi波段的光谱预测值，L(xi，yi，t3，b_oi)表示t3时刻目标影像中相似像素点(xi，yi)处的b_oi波段的光谱值。
然后进行基于光谱-时间信息的预测。在一些实施例中，基于相似像素点集合的光谱时间变化信息，对目标像素点的光谱值进行预测，得到第二光谱预测值，具体为：计算t1时刻在最佳光学参考影像中相似像素点集合中各像素点与目标像素点的光谱值之差，以得到第二光谱差值；计算t2时刻相似像素点集合中各像素点与目标像素点在最佳SAR参考影像中的光谱值之差，以得到第三光谱差值；分别基于第二光谱差值、第三光谱差值进行计算，对应得到相似像素点集合中各像素点的第二权重、第三权重；根据第二权重、第三权重，以及相似像素点集合在t1时刻、t2时刻、t3时刻的光谱值，对目标像素点的光谱值进行预测，得到第二光谱预测值。
其中，t1时刻为最佳光学参考影像的成像时间，t3时刻为目标影像的成像时间；t2时刻为最佳SAR参考影像的成像时间。
本申请实施例中，分别就相似像素点集合为高质量相似像素点以及相似像素点为一般相似像素点两种情况，对目标像素点的光谱值进行基于光谱-时间信息的预测。
下面详细说明相似像素点集合为高质量像素点时，基于光谱-时间信息对目标像素点进行预测的步骤。
为了方便描述，将最佳光学参考影像的成像时间t1到目标影像的成像时间t3的时间间隔称为第一时间段，将最佳SAR参考影像的成像时间t2到目标影像的成像时间t3的时间间隔称为第二时间段，由于高质量相似像素点在t1时刻到t3时刻始终与目标像素点保持较高的相似性，因此，高质量相似像素点在不同时间段内的变化也与目标相似像素点的变化基本保持相同，即高质量相似像素点分别在第一时间段、第二时间段内的光谱变化值与目标像素点在第一时间段、第二时间段内的光谱变化值将对应保持较高的相似性。基于此，本申请实施例中，分别根据第一时间段、第二时间段内高质量相似像素点的光谱变化值，推测相似像素点的光谱变化值，然后根据得到的光谱变化值分别推测t3时刻目标像素点的光谱值，最后对两次分别推算到的t3时刻目标像素点的光谱值进行加权平均计算得到基于光谱-时间信息的第二光谱预测值。
具体地，首先计算第二光谱差值，即t1时刻高质量相似像素点位置处的光谱值与目标像素点位置处的光谱值之差，也就是t1时刻高质量像素点与目标像素点之间的光谱相似性，这里仍然采用均方根误差对第二光谱差值进行表达，其计算公式如下：
式中，RMSDi_t1表示第二光谱差值，(xi，yi)表示第i个高质量相似像素点的坐标，(x，y)表示目标像素点的坐标。
然后，计算第三光谱差值，即t2时刻在最佳SAR参考影像中高质量相似像素点与目标像素点的光谱值之差，计算公式如下：
接着，对RMSDi_t1、RMSDi_t2分别进行归一化处理以消除量纲差异，其计算公式如下：
式中，表示对RMSDi_t1进行归一化处理后的结果，/>表示对RMSDi_t2进行归一化处理的结果，RMSDi_t1_min、RMSDi_t1_max分别表示RMSDi_t1的最小值、最大值，RMSDi_t2_min、RMSDi_t2_max分别表示RMSDi_t2的最小值、最大值。
随后，结合上述步骤得到的相似像素点与目标像素点之间的空间距离，分别基于归一化处理后的第二光谱差值、第三光谱差值，进行计算，对应得到高质量相似像素点的第二权重、第三权重，其计算公式如下：
式中，RMSDi_t1、RMSDi_t2分别代表在t1、t2时刻高质量相似像素点的光谱相似性，分别代表归一化处理后的光谱相似性，Wi_t1、Wi_t2分别代表在t1、t2时刻高质量相似像素点的权重，/>分别表示第i个、第j个高质量相似像素点与目标像素点的归一化处理后的空间距离。
接着，根据高质量相似像素点在t1时刻、t2时刻、t3时刻的光谱值，基于第二权重、第三权重分别推测目标像素点在t3时刻的光谱值，其计算公式如下：
式中，Lst_t1(x,y,t3,b_oi)、Lst_t2(x,y,t3,b_oi)分别表示基于第一时间段、第二时间段高质量相似像素点的光谱变化值预测得到的目标像素点的光谱值。
下一步，根据高质量相似像素点在t1、t2时刻与t3时刻的光谱差值，为Lst_t1(x,y,t3,b_oi)、Lst_t2(x,y,t3,b_oi)分别赋予不同的权重，具体计算公式如下：
式中，Rt1、Rt2分别代表t1、t2时刻高质量相似像素点与t3时刻高质量相似像素点的光谱差值(也称为相似性)。
需要说明的是，在公式(20)中，由于t1时刻的最佳光学参考影像、t3时刻的目标影像均为光学影像，其波段可一一对应，因此，在计算公式(20)时，直接使用相同波段的相同位置，即一个高质量相似像素点处的不同时刻(即t1、t3时刻)的光谱值进行相减即可得到Rt1。而公式(21)中，由于t2是最佳SAR参考影像的成像时间，t3是目标影像的成像时间，也就是说，t2、t3分别对应光学影像和SAR影像，而光学影像与SAR影像的波段无法一一对应，因此，在高质量相似像素点处，将t3时刻的目标影像的波段逐个与最佳SAR参考影像中的波段进行相减运算，以得到高质量相似像素点处t2时刻与t3时刻的光谱差值Rt2。
最后，根据Rt1、Rt2计算基于光谱-时间信息对目标像素点的最终光谱预测值，即第二光谱预测值，如公式(22)所示，公式(22)如下：
以上为在相似像素点集合为高质量相似像素点的情况下，基于光谱-时间信息对目标像素点进行预测的步骤。
当相似像素点集合为一般相似像素点时，也就是目标像素点没有高质量相似像素点的情况下，那么，只需将上述公式(12)～公式(22)中有关最佳光学参考影像的参数去掉即可，其余计算步骤与高质量相似像素点的步骤相同，即可基于光谱-时间信息对目标像素点进行预测以得到第二光谱预测值。
在对目标像素点分别进行基于光谱-空间信息的预测和基于光谱-时间信息的预测之后，在一些实施例中，分别计算第一光谱预测值的权重、第二光谱预测值的权重，以确定目标像素点最终的光谱预测值，从而得到目标影像对应的云修复影像，具体为：根据相似像素点集合与目标像素点的平均光谱差值Rss，计算第一光谱预测值的权重；根据相似像素点集合的综合平均光谱差值Rst，计算第二光谱预测值的权重；根据第一光谱预测值、第一光谱预测值的权重、第二光谱预测值、第二光谱预测值的权重，计算目标像素点最终的光谱预测值，从而得到目标影像对应的云修复影像。
本申请实施例中，基于光谱-空间信息的预测主要基于目标影像中的高质量相似点在不同空间距离内光谱变化情况进行预测，因此，其辐射测量更加一致，即云遮挡区域与其周围区域可以保持较好的连贯性。基于光谱-时间信息的预测通过利用高质量相似像素点在不同时间段内的光谱变化情况，成功捕捉到了地物的变化信息，使得预测结果更加贴近地物的演变进程。因此，为综合考虑地物景观的自然连贯性与演变情况，通过将这两个结果加权组合以获得精度更高的预测结果。
首先，根据相似像素点集合与目标像素点的平均光谱差值Rss，计算第一光谱预测值的权重，也就是说，利用高质量相似像素点与目标像素点的平均光谱差值Rss来表示地物的连贯程度，其计算公式如下：
然后，根据相似像素点集合的综合平均光谱差值Rst，计算第二光谱预测值的权重，也可以理解为，利用高质量相似像素点或者一般相似像素点在第一时间段和第二时间段的综合平均光谱差值来表示地物的变化程度，其计算公式如下：
式中，由于第一时间段和第二时间段是两个不同的时间段，也就会产生不同的平均光谱差值，因此需根据相似像素点在t1或t2时刻与t3时刻的光谱差值情况分别赋予不同的权重，经过平均加权计算后得到Rst。具体来说，基于第一时间段的光谱差值Rt1计算得到对应的平均光谱差值，即公式(24)等号右端第一项，基于第二时间段内的光谱差值Rt2计算得到对应的平均光谱差值，即公式(24)等号右端第二项，然后经平均加权计算后得到Rst。
接着，根据平均光谱差值Rss、第一时间段和第二时间段的综合平均光谱差值Rst，分别计算得到第一光谱预测值的权重和第二光谱预测值的权重。
最后，根据第一光谱预测值Lss、第一光谱预测值的权重、第二光谱预测值Lst、第二光谱预测值的权重，计算目标像素点最终的光谱预测值，从而得到目标影像对应的云修复影像，即根据地物的连贯程度和变化程度权重计算得到目标像素点的最终预测值Lt3，计算公式如下：
式中，表示第一光谱预测值的权重，/>表示第二光谱预测值的权重。
可以理解，在云遮挡区域范围较大的情况下，可能会出现极个别的目标像素点既无高质量相似像素点又无一般像素点的特殊情况，此时，本申请的一些实施方式中，使用局部线性直方图匹配来对该目标像素点进行预测。
具体地，首先，通过公共像素的平均值和标准偏差，计算获得增益(gain)和偏差(bias)，对应的公式如下：
bias＝μP-μF*gain (27)，
式中，μP表示目标影像中公共像素的平均值，μF表示t2时刻遥感影像中公共像素的平均值，公共像素指的是将目标影像与t2时刻遥感影像进行求交集运算得到的像素集合，σP表示目标影像中公共像素的标准偏差，σF表示t2时刻遥感影像中公共像素的标准偏差。
然后，将增益与t2时刻遥感影像中目标像素点处的光谱值进行相乘，然后加上偏差，最终获得目标像素点的预测值。对应的计算公式如下：
Lt3(x,y,t3,b_oi)＝gain*L(x,y,t2,b_oi)+bias (28)，
式中，L(x，y，t2，b_oi)表示t2时刻遥感影像中目标像素点处的光谱值。
在一个具体的示例中，为了验证本申请所提供的技术方案在多场景下的通用性和适应能力，分别选择了自然因素影响下大部分区域土地覆盖类型发生分异变化场景(以下简称场景1)、人为因素影响下大部分区域土地覆盖类型发生分异变化场景(以下简称场景2)、大部分区域土地覆盖类型发生相同变化场景(以下简称场景3)、大部分区域土地覆盖类型无变化场景(以下简称场景4)等进行实验。
以场景1为例，实验所在目标区域为我国J省某市，目标影像为Landsat 8影像，成像时间t3为2019年9月20日，云遮挡区域面积约为35km2。在目标影像成像时间之前且与其时间间隔最短的SAR影像(即最佳SAR参考影像)为Sentinel-1影像，成像时间t2为2019年9月18日。在目标影像成像时间之前且与其时间间隔最短的光学遥感影像(即最佳光学参考影像)为Sentinel-2影像，成像时间t1为2019年9月9日。土地覆盖数据来自Esri_Land_Cover数据集。
对上述遥感影像进行可视化显示后，可以看出：云遮挡区域以及相邻区域在成像时间t1到成像时间t3时间段内(即第一时间段)大部分区域的土地覆盖类型发生了不同程度的变化，具体来说，在第一时间段内，可视化窗口中，区遮挡区域及其相邻区域的左上部的土地覆盖类型变化表现为：滩涂、沙洲逐渐加高，陆续生长出一些植被；区遮挡区域及其相邻区域的右下部的土地覆盖由水体转变为滩涂、沙洲；而左上部与右下部之间的区域仍为水体，土地覆盖类型未发生明显变化，可见，在自然因素驱动下，云遮挡区域及其相邻区域的土地覆盖类型均发生了分异性变化。
基于上述数据，分别采用改进的邻域相似像素插值器(MNSPI)算法、MNSPI_SAR算法、以及本申请提供的MUCRMD对目标影像进行云修复，并采用RMSE(Root Mean SquareError)对云修复的结果进行精度评估，各波段的精度如表1所示，表1如下：
表1 MNSPI、MNSPI_SAR、MUCRMD的精度评估对比
其中，MNSPI_SAR算法为利用SAR影像寻找相似像素点，进而再利用SAR图像中相同位置的无云光学像素替换有云光学像素的方法对t3时刻云遮挡区域进行预测。由于该方法与MNSPI具有一定的相似性，只不过参考数据由光学影像变为了SAR影像，因此将该方法简称为MNSPI_SAR算法。从上表可以看出，在自然因素影响下大部分区域土地覆盖类型发生分异变化场景下，基于本申请提供的多模态数据的多场景通用性遥感影像云修复方法得到的各个波段预测结果的精度均有显著提高。
对场景2、场景3、场景4进行实验，均得到与场景1类似的结果，在此不再一一列举。
综上所述，本申请的技术方案中，通过获取目标区域的多模态遥感影像和目标影像；多模态遥感影像包括来自不同卫星的光学遥感影像和来自不同卫星的合成孔径雷达SAR影像；目标影像为包含云遮挡区域的遥感影像；然后基于多模态遥感影像的不同成像时间，分别从光学遥感影像和SAR影像中确定目标影像的最佳光学参考影像以及目标影像的最佳SAR参考影像；接着基于目标影像、最佳光学参考影像、最佳SAR参考影像以及预先获取的目标区域的土地覆盖数据，确定相似像素点集合；依次对相似像素点集合中的各个像素点进行权重分配和加权计算，得到目标影像对应的云修复影像。由于采用不同卫星、不同成像时间的光学影像与SAR影像作为参考影像，并引入土地覆盖数据对相似像素点进行筛选，使得相似像素点能够更综合地表达出目标像素点处地物景观的自然连贯性和演变情况。实验结果表明，无论目标区域的土地覆盖是否突变、发生何等程度的变化、土地覆盖突变是否分异等多种不同场景下，通过本申请的技术方案均可取得精度更高、效果更好的云修复效果。
本申请实施例提供的方法综合利用目前可广泛免费获得的中分辨率遥感卫星(如Landsat-8、Landsat-9、Sentinel-1和Sentinel-2)组建虚拟星座，借助虚拟星座中不同卫星的成像时间与重访周期差异，确定最佳光学、SAR参考影像，然后通过引入土地覆盖数据确定高质量相似像素点，并通过对高质量相似像素点的合理权重分配与加权计算重建了缺失信息。结果表明，无论土地覆盖是否突变，发生何种程度的变化，土地覆盖突变是否分异等多种不同场景下，与现有的云修复方法相比，本申请实施例所提供的方法均可取得精度更高、效果更佳的云修复结果，证明了该方法不再受传统遥感影像云修复方法的单一场景限制，可以适用于不同的具体应用场景。
示例性系统
本申请实施例提供一种多模态数据的多场景通用性遥感影像云修复系统，如图3所示，该系统包括：获取单元301、确定单元302、筛选单元303和修复单元304。其中：
获取单元301，配置为获取目标区域的多模态遥感影像和目标影像；多模态遥感影像包括来自至少一个卫星的光学遥感影像和来自至少一个卫星的合成孔径雷达SAR影像；目标影像为包含云遮挡区域的遥感影像。
确定单元302，配置为基于多模态遥感影像的不同成像时间，分别从光学遥感影像和SAR影像中确定目标影像的最佳光学参考影像以及目标影像的最佳SAR参考影像。
筛选单元303，配置为基于目标影像、最佳光学参考影像、最佳SAR参考影像以及预先获取的目标区域的土地覆盖数据，确定相似像素点集合。
修复单元304，配置为依次对相似像素点集合中的各个像素点进行权重分配和加权计算，得到目标影像对应的云修复影像。
本申请实施例提供的多模态数据的多场景通用性遥感影像云修复系统能够实现上述任一实施例提供的多模态数据的多场景通用性遥感影像云修复方法的流程、步骤，并达到相同的技术效果，在此不再一一赘述。
以上所述仅为本申请的优选实施例，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (3)

1.一种多模态数据的多场景通用性遥感影像云修复方法，其特征在于，包括：
获取目标区域的多模态遥感影像和目标影像；所述多模态遥感影像包括来自不同卫星的光学遥感影像和来自不同卫星的合成孔径雷达SAR影像；所述目标影像为包含云遮挡区域的遥感影像；
基于所述多模态遥感影像的不同成像时间，分别从所述光学遥感影像和SAR影像中确定所述目标影像的最佳光学参考影像以及所述目标影像的最佳SAR参考影像；
基于所述目标影像、所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据，确定相似像素点集合；
依次对所述相似像素点集合中的各个像素点进行权重分配和加权计算，得到所述目标影像对应的云修复影像；
所述基于所述目标影像、所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据，确定相似像素点集合，具体为：
基于所述目标影像中云遮挡区域的中心点到云遮挡区域的边界之间的最短距离，取该最短距离的一半作为云遮挡区域向外扩张的缓冲半径，通过缓冲得到云遮挡区域的缓冲区，将该缓冲区作为相似像素的搜索区域；
根据所述搜索区域，对所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据进行求交集运算，以确定相似像素点集合；
所述相似像素点集合包括高质量相似像素点集合；
对应地，根据所述搜索区域，对所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据进行求交集运算，以确定相似像素点集合，具体为：
分别计算目标像素点与所述最佳光学参考影像在所述搜索区域内的所有像素点、所述目标像素点与所述最佳SAR参考影像在所述搜索区域内的所有像素点之间的光谱相似性，对应得到第一光谱相似性和第二光谱相似性；其中，所述目标像素点为所述目标影像中云遮挡区域内任一像素点；
根据预设相似性条件以及所述第一光谱相似性、所述第二光谱相似性，分别对所述最佳光学参考影像在所述搜索区域内的所有像素点、所述最佳SAR参考影像在所述搜索区域内的所有像素点进行分别筛选，对应得到第一备选相似像素点集合、第二备选相似像素点集合；
在所述搜索区域内从所述土地覆盖数据中搜索与所述目标像素点具有相同土地覆盖类型的像素点，以得到第三备选相似像素点集合；
将所述第一备选相似像素点集合、所述第二备选相似像素点集合与所述第三备选相似像素点集合进行求交集运算，得到所述高质量相似像素点集合；
所述依次对所述相似像素点集合中的各个像素点进行权重分配和加权计算，得到所述目标影像对应的云修复影像，具体为：
基于所述相似像素点集合的光谱空间变化信息与空间位置，对所述目标像素点的光谱值进行预测，得到第一光谱预测值；
基于所述相似像素点集合的光谱时间变化信息，对所述目标像素点的光谱值进行预测，得到第二光谱预测值；
分别计算所述第一光谱预测值的权重、所述第二光谱预测值的权重，以确定所述目标像素点最终的光谱预测值，从而得到所述目标影像对应的云修复影像；
所述基于所述相似像素点集合的光谱空间变化信息与空间位置，对所述目标像素点的光谱值进行预测，得到第一光谱预测值，具体为：
根据所述最佳光学参考影像和所述最佳SAR参考影像分别在所述相似像素点集合的各像素点处的光谱值，以及所述最佳光学参考影像和所述最佳SAR参考影像分别在目标像素点位置处的光谱值，计算得到第一光谱差值；
采用相似像素点与目标像素点之间的RMSD表示第一光谱差值，其计算公式如下：
，
式中，最佳光学参考影像的波段总数为n_oi个，n_oi个波段中的第b_oi个波段用b_oi表示；t 1 为最佳光学参考影像的成像时间，L（x i ，y i ，t 1 ，b_oi）表示t 1 时刻的最佳光学参考影像中位于第b_oi个波段的（x i ，y i ）处的光谱值，L（x，y，t 1 ，b_oi）表示t 1 时刻的最佳光学参考影像中位于b_oi波段的目标像素点（x，y）处的光谱值；最佳SAR参考影像的波段数量为n_ sar个，n_sar个波段中的第b_sar个波段用b_sar表示；t 2 为最佳SAR参考影像的成像时间，L （x i ，y i ，t 2 ，b_sar）表示t 2 时刻的最佳SAR参考影像中位于b_sar波段的（x i ，y i ）处的光谱值，L （x，y，t 2 ，b_sar）表示t 2 时刻的最佳SAR参考影像中位于b_sar波段的目标像素点（x，y）处的光谱值；RMSD i 表示相似像素点（x i ，y i ）与目标像素点（x，y）的第一光谱差值；
根据所述相似像素点集合中各像素点的位置、所述目标像素点的位置，计算得到空间距离；其计算公式如下：
，
式中，D i 表示相似像素点（x i ，y i ）与目标像素点（x，y）的空间距离；
分别对所述第一光谱差值、所述空间距离进行归一化处理，并计算得到所述相似像素点集合中各像素点的第一权重；
对第一光谱差值进行归一化处理的计算公式如下：
，
式中，表示归一化处理后的第一光谱差值，RMSD i_max 、RMSD i_min 分别表示第一光谱差值的最大值、最小值；
对空间距离进行归一化处理的计算公式如下：
，
式中，表示归一化处理后的空间距离，D i_max 、D i_min 分别空间距离的最大值、最小值；
第一权重的表达式如下：
，
式中，W i 表示相似像素点（x i ，y i ）的第一权重，N表示相似像素点的数量；
根据所述目标影像在所述相似像素点集合的各像素点处的光谱值以及所述相似像素点集合中各像素点的第一权重，对所述目标像素点的光谱值进行预测，得到第一光谱预测值；
第一光谱预测值的计算公式如下：
，
式中，L ss （x，y，t 3 ，b_oi）表示第一光谱预测值，也就是t 3 时刻目标像素点（x，y）处的b_oi波段的光谱预测值，L（x i ，y i ，t 3 ，b_oi）表示t 3 时刻目标影像中相似像素点（x i ，y i ）处的b_oi波段的光谱值；
所述基于所述相似像素点集合的光谱时间变化信息，对所述目标像素点的光谱值进行预测，得到第二光谱预测值，具体为：
计算t 1 时刻在最佳光学参考影像中相似像素点集合中各像素点与目标像素点的光谱值之差，以得到第二光谱差值；第二光谱差值的计算公式如下：
，
式中，RMSD i_t1 表示第二光谱差值，（x i ，y i ）表示第i个高质量相似像素点的坐标，（x，y）表示目标像素点的坐标；
计算t 2 时刻在最佳SAR参考影像中相似像素点集合中各像素点与目标像素点的光谱值之差，以得到第三光谱差值；第三光谱差值的计算公式如下：
，
对、/>分别进行归一化处理，其计算公式如下：
，
，
式中，表示对/>进行归一化处理后的结果，/>表示对进行归一化处理的结果，/>、/>分别表示的最小值、最大值，/>、/>分别表示/>的最小值、最大值；
分别基于所述第二光谱差值、所述第三光谱差值进行计算，对应得到所述相似像素点集合中各像素点的第二权重、第三权重；其计算公式如下：
，
，
式中，RMSD i_t1 、RMSD i_t2 分别代表在t1、t2时刻高质量相似像素点的光谱相似性，、/>分别代表归一化处理后的光谱相似性，W i_t1 、W i_t2 分别代表在t1、t2时刻高质量相似像素点的权重，/>、/>分别表示第i个、第j个高质量相似像素点与目标像素点的归一化处理后的空间距离；
根据高质量相似像素点在t 1 时刻、t 2 时刻、t 3 时刻的光谱值，基于第二权重、第三权重分别推测目标像素点在t 3 时刻的光谱值，其计算公式如下：
，
，
式中，、/>分别表示基于第一时间段、第二时间段高质量相似像素点的光谱变化值预测得到的目标像素点的光谱值；
其中，所述t 1 时刻为所述最佳光学参考影像的成像时间；所述t 2 时刻为所述最佳SAR参考影像的成像时间；所述t 3 时刻为所述目标影像的成像时间；第一时间段为最佳光学参考影像的成像时间t 1 到目标影像的成像时间t 3 的时间间隔；第二时间段为最佳SAR参考影像的成像时间t 2 到目标影像的成像时间t 3 的时间间隔；
计算高质量相似像素点在t 1 、t 2 时刻与t 3 时刻的光谱差值，具体计算公式如下：
，
，
式中，R t1 、R t2 分别代表t 1 、t 2 时刻高质量相似像素点与t 3 时刻高质量相似像素点的光谱差值；
根据R t1 、R t2 计算第二光谱预测值；其计算公式如下：
，
所述分别计算所述第一光谱预测值的权重、所述第二光谱预测值的权重，以确定所述目标像素点最终的光谱预测值，从而得到所述目标影像对应的云修复影像，具体为：
根据所述相似像素点集合与所述目标像素点的平均光谱差值R ss ，计算所述第一光谱预测值的权重；
根据所述相似像素点集合的综合平均光谱差值R st ，计算所述第二光谱预测值的权重；
平均光谱差值R ss 的计算公式如下：
，
综合平均光谱差值R st 的计算公式如下：
，
第一光谱预测值的权重的表达式如下：
，
第二光谱预测值的权重的表达式如下：
，
根据所述第一光谱预测值、所述第一光谱预测值的权重、所述第二光谱预测值、所述第二光谱预测值的权重，计算所述目标像素点最终的光谱预测值，从而得到所述目标影像对应的云修复影像；
所述目标像素点最终的光谱预测值的计算公式如下：
，
式中，表示第一光谱预测值的权重，/>表示第二光谱预测值的权重。
2.根据权利要求1所述的多模态数据的多场景通用性遥感影像云修复方法，其特征在于，所述基于所述多模态遥感影像的不同成像时间，分别从所述光学遥感影像和SAR影像中确定所述目标影像的最佳光学参考影像以及所述目标影像的最佳SAR参考影像，具体为：
分别从所述光学遥感影像、SAR影像中选择成像时间在所述目标影像的成像时间之前，且与所述目标影像的成像时间间隔最短的无云影像；
将从所述光学遥感影像中选出的无云影像作为所述目标影像的最佳光学参考影像；将从所述SAR影像中选出的无云影像作为所述目标影像的最佳SAR参考影像。
3.一种多模态数据的多场景通用性遥感影像云修复系统，其特征在于，包括：
获取单元，配置为获取目标区域的多模态遥感影像和目标影像；所述多模态遥感影像包括来自至少一个卫星的光学遥感影像和来自至少一个卫星的合成孔径雷达SAR影像；所述目标影像为包含云遮挡区域的遥感影像；
确定单元，配置为基于所述多模态遥感影像的不同成像时间，分别从所述光学遥感影像和SAR影像中确定所述目标影像的最佳光学参考影像以及所述目标影像的最佳SAR参考影像；
筛选单元，配置为基于所述目标影像、所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据，确定相似像素点集合；
修复单元，配置为依次对所述相似像素点集合中的各个像素点进行权重分配和加权计算，得到所述目标影像对应的云修复影像；
所述基于所述目标影像、所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据，确定相似像素点集合，具体为：
基于所述目标影像中云遮挡区域的中心点到云遮挡区域的边界之间的最短距离，取该最短距离的一半作为云遮挡区域向外扩张的缓冲半径，通过缓冲得到云遮挡区域的缓冲区，将该缓冲区作为相似像素的搜索区域；
根据所述搜索区域，对所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据进行求交集运算，以确定相似像素点集合；
所述相似像素点集合包括高质量相似像素点集合；
对应地，根据所述搜索区域，对所述最佳光学参考影像、所述最佳SAR参考影像以及预先获取的所述目标区域的土地覆盖数据进行求交集运算，以确定相似像素点集合，具体为：
分别计算目标像素点与所述最佳光学参考影像在所述搜索区域内的所有像素点、所述目标像素点与所述最佳SAR参考影像在所述搜索区域内的所有像素点之间的光谱相似性，对应得到第一光谱相似性和第二光谱相似性；其中，所述目标像素点为所述目标影像中云遮挡区域内任一像素点；
根据预设相似性条件以及所述第一光谱相似性、所述第二光谱相似性，分别对所述最佳光学参考影像在所述搜索区域内的所有像素点、所述最佳SAR参考影像在所述搜索区域内的所有像素点进行分别筛选，对应得到第一备选相似像素点集合、第二备选相似像素点集合；
在所述搜索区域内从所述土地覆盖数据中搜索与所述目标像素点具有相同土地覆盖类型的像素点，以得到第三备选相似像素点集合；
将所述第一备选相似像素点集合、所述第二备选相似像素点集合与所述第三备选相似像素点集合进行求交集运算，得到所述高质量相似像素点集合；
所述依次对所述相似像素点集合中的各个像素点进行权重分配和加权计算，得到所述目标影像对应的云修复影像，具体为：
基于所述相似像素点集合的光谱空间变化信息与空间位置，对所述目标像素点的光谱值进行预测，得到第一光谱预测值；
基于所述相似像素点集合的光谱时间变化信息，对所述目标像素点的光谱值进行预测，得到第二光谱预测值；
分别计算所述第一光谱预测值的权重、所述第二光谱预测值的权重，以确定所述目标像素点最终的光谱预测值，从而得到所述目标影像对应的云修复影像；
所述基于所述相似像素点集合的光谱空间变化信息与空间位置，对所述目标像素点的光谱值进行预测，得到第一光谱预测值，具体为：
根据所述最佳光学参考影像和所述最佳SAR参考影像分别在所述相似像素点集合的各像素点处的光谱值，以及所述最佳光学参考影像和所述最佳SAR参考影像分别在目标像素点位置处的光谱值，计算得到第一光谱差值；
采用相似像素点与目标像素点之间的RMSD表示第一光谱差值，其计算公式如下：
，
式中，最佳光学参考影像的波段总数为n_oi个，n_oi个波段中的第b_oi个波段用b_oi表示；t 1 为最佳光学参考影像的成像时间，L（x i ，y i ，t 1 ，b_oi）表示t 1 时刻的最佳光学参考影像中位于第b_oi个波段的（x i ，y i ）处的光谱值，L（x，y，t 1 ，b_oi）表示t 1 时刻的最佳光学参考影像中位于b_oi波段的目标像素点（x，y）处的光谱值；最佳SAR参考影像的波段数量为n_ sar个，n_sar个波段中的第b_sar个波段用b_sar表示；t 2 为最佳SAR参考影像的成像时间，L （x i ，y i ，t 2 ，b_sar）表示t 2 时刻的最佳SAR参考影像中位于b_sar波段的（x i ，y i ）处的光谱值，L （x，y，t 2 ，b_sar）表示t 2 时刻的最佳SAR参考影像中位于b_sar波段的目标像素点（x，y）处的光谱值；RMSD i 表示相似像素点（x i ，y i ）与目标像素点（x，y）的第一光谱差值；
根据所述相似像素点集合中各像素点的位置、所述目标像素点的位置，计算得到空间距离；其计算公式如下：
，
式中，D i 表示相似像素点（x i ，y i ）与目标像素点（x，y）的空间距离；
分别对所述第一光谱差值、所述空间距离进行归一化处理，并计算得到所述相似像素点集合中各像素点的第一权重；
对第一光谱差值进行归一化处理的计算公式如下：
，
式中，表示归一化处理后的第一光谱差值，RMSD i_max 、RMSD i_min 分别表示第一光谱差值的最大值、最小值；
对空间距离进行归一化处理的计算公式如下：
，
式中，表示归一化处理后的空间距离，D i_max 、D i_min 分别空间距离的最大值、最小值；
第一权重的表达式如下：
，
式中，W i 表示相似像素点（x i ，y i ）的第一权重，N表示相似像素点的数量；
根据所述目标影像在所述相似像素点集合的各像素点处的光谱值以及所述相似像素点集合中各像素点的第一权重，对所述目标像素点的光谱值进行预测，得到第一光谱预测值；
第一光谱预测值的计算公式如下：
，
式中，L ss （x，y，t 3 ，b_oi）表示第一光谱预测值，也就是t 3 时刻目标像素点（x，y）处的b_oi波段的光谱预测值，L（x i ，y i ，t 3 ，b_oi）表示t 3 时刻目标影像中相似像素点（x i ，y i ）处的b_oi波段的光谱值；
所述基于所述相似像素点集合的光谱时间变化信息，对所述目标像素点的光谱值进行预测，得到第二光谱预测值，具体为：
计算t 1 时刻在最佳光学参考影像中相似像素点集合中各像素点与目标像素点的光谱值之差，以得到第二光谱差值；第二光谱差值的计算公式如下：
，
式中，RMSD i_t1 表示第二光谱差值，（x i ，y i ）表示第i个高质量相似像素点的坐标，（x，y）表示目标像素点的坐标；
计算t 2 时刻在最佳SAR参考影像中相似像素点集合中各像素点与目标像素点的光谱值之差，以得到第三光谱差值；第三光谱差值的计算公式如下：
，
对、/>分别进行归一化处理，其计算公式如下：
，
，
式中，表示对/>进行归一化处理后的结果，/>表示对进行归一化处理的结果，/>、/>分别表示的最小值、最大值，/>、/>分别表示的最小值、最大值；
分别基于所述第二光谱差值、所述第三光谱差值进行计算，对应得到所述相似像素点集合中各像素点的第二权重、第三权重；其计算公式如下：
，
，
式中，RMSD i_t1 、RMSD i_t2 分别代表在t1、t2时刻高质量相似像素点的光谱相似性，、/>分别代表归一化处理后的光谱相似性，W i_t1 、W i_t2 分别代表在t1、t2时刻高质量相似像素点的权重，/>、/>分别表示第i个、第j个高质量相似像素点与目标像素点的归一化处理后的空间距离；
根据高质量相似像素点在t 1 时刻、t 2 时刻、t 3 时刻的光谱值，基于第二权重、第三权重分别推测目标像素点在t 3 时刻的光谱值，其计算公式如下：
，
，
式中，、/>分别表示基于第一时间段、第二时间段高质量相似像素点的光谱变化值预测得到的目标像素点的光谱值；
其中，所述t 1 时刻为所述最佳光学参考影像的成像时间；所述t 2 时刻为所述最佳SAR参考影像的成像时间；所述t 3 时刻为所述目标影像的成像时间；第一时间段为最佳光学参考影像的成像时间t 1 到目标影像的成像时间t 3 的时间间隔；第二时间段为最佳SAR参考影像的成像时间t 2 到目标影像的成像时间t 3 的时间间隔；
计算高质量相似像素点在t 1 、t 2 时刻与t 3 时刻的光谱差值，具体计算公式如下：
，
，
式中，R t1 、R t2 分别代表t 1 、t 2 时刻高质量相似像素点与t 3 时刻高质量相似像素点的光谱差值；
根据R t1 、R t2 计算第二光谱预测值；其计算公式如下：
，
所述分别计算所述第一光谱预测值的权重、所述第二光谱预测值的权重，以确定所述目标像素点最终的光谱预测值，从而得到所述目标影像对应的云修复影像，具体为：
根据所述相似像素点集合与所述目标像素点的平均光谱差值R ss ，计算所述第一光谱预测值的权重；
根据所述相似像素点集合的综合平均光谱差值R st ，计算所述第二光谱预测值的权重；
平均光谱差值R ss 的计算公式如下：
，
综合平均光谱差值R st 的计算公式如下：
，
第一光谱预测值的权重的表达式如下：
，
第二光谱预测值的权重的表达式如下：
，
根据所述第一光谱预测值、所述第一光谱预测值的权重、所述第二光谱预测值、所述第二光谱预测值的权重，计算所述目标像素点最终的光谱预测值，从而得到所述目标影像对应的云修复影像；
所述目标像素点最终的光谱预测值的计算公式如下：
，
式中，表示第一光谱预测值的权重，/>表示第二光谱预测值的权重。/>
CN202310130461.XA 2023-02-08 2023-02-08 多模态数据的多场景通用性遥感影像云修复方法和系统 Active CN116245757B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310130461.XA CN116245757B (zh) 2023-02-08 2023-02-08 多模态数据的多场景通用性遥感影像云修复方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310130461.XA CN116245757B (zh) 2023-02-08 2023-02-08 多模态数据的多场景通用性遥感影像云修复方法和系统

Publications (2)

Publication Number Publication Date
CN116245757A CN116245757A (zh) 2023-06-09
CN116245757B true CN116245757B (zh) 2023-09-19

Family

ID=86625691

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310130461.XA Active CN116245757B (zh) 2023-02-08 2023-02-08 多模态数据的多场景通用性遥感影像云修复方法和系统

Country Status (1)

Country Link
CN (1) CN116245757B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117576236A (zh) * 2023-11-13 2024-02-20 宁波大学 联合多时相sar和光学信息的缺失光学图像重建方法

Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104156923A (zh) * 2014-08-12 2014-11-19 西北工业大学 一种基于稀疏表示的多光谱遥感图像去云方法
CN107590782A (zh) * 2017-08-21 2018-01-16 西北工业大学 一种基于全卷积网络的高分辨率光学图像厚云去除方法
CN108280810A (zh) * 2018-01-09 2018-07-13 北方工业大学 一种单时相光学遥感图像云覆盖区修复的自动处理方法
CN108917933A (zh) * 2018-04-10 2018-11-30 中国科学院上海技术物理研究所 一种基于地面高温点的宽覆盖相机几何定位方法
CN110287898A (zh) * 2019-06-27 2019-09-27 苏州中科天启遥感科技有限公司 一种光学卫星遥感影像云检测方法
CN111881816A (zh) * 2020-07-27 2020-11-03 安徽省（水利部淮河水利委员会）水利科学研究院（安徽省水利工程质量检测中心站） 一种长时序的河湖围埂养殖区域监测方法
CN112381815A (zh) * 2020-11-26 2021-02-19 浙江易智信息技术有限公司 多维融合遥感大数据去云方法
CN113850139A (zh) * 2021-08-25 2021-12-28 南京林业大学 一种基于多源遥感的森林年际物候监测方法
CN114119617A (zh) * 2021-11-12 2022-03-01 武汉大学 一种多光谱卫星遥感影像的内陆盐湖卤虫带提取方法
CN114140701A (zh) * 2021-12-10 2022-03-04 广西师范大学 一种多尺度多时相遥感影像去云的方法
CN114202698A (zh) * 2021-12-15 2022-03-18 南京大学 联合主被动遥感数据的水体精准提取方法
CN114418970A (zh) * 2021-12-31 2022-04-29 中国科学院空天信息创新研究院 基于卫星遥感的雾霾分布与气溶胶光学厚度检测方法和装置
CN114511786A (zh) * 2022-04-20 2022-05-17 中国石油大学(华东) 融合多时相信息和分通道密集卷积的遥感图像去云方法
CN115131674A (zh) * 2022-06-24 2022-09-30 武汉大学 基于深度低秩网络的多时相光学遥感影像云检测方法
CN115165773A (zh) * 2022-08-01 2022-10-11 上海兰桂骐技术发展股份有限公司 一种基于Google Earth Engine的水稻面积提取方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9721181B2 (en) * 2015-12-07 2017-08-01 The Climate Corporation Cloud detection on remote sensing imagery
WO2020160641A1 (en) * 2019-02-04 2020-08-13 Farmers Edge Inc. Shadow and cloud masking for remote sensing images in agriculture applications using multilayer perceptron
US20230026811A1 (en) * 2021-07-15 2023-01-26 Ping An Technology (Shenzhen) Co., Ltd. System and method for removing haze from remote sensing images

Patent Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104156923A (zh) * 2014-08-12 2014-11-19 西北工业大学 一种基于稀疏表示的多光谱遥感图像去云方法
CN107590782A (zh) * 2017-08-21 2018-01-16 西北工业大学 一种基于全卷积网络的高分辨率光学图像厚云去除方法
CN108280810A (zh) * 2018-01-09 2018-07-13 北方工业大学 一种单时相光学遥感图像云覆盖区修复的自动处理方法
CN108917933A (zh) * 2018-04-10 2018-11-30 中国科学院上海技术物理研究所 一种基于地面高温点的宽覆盖相机几何定位方法
CN110287898A (zh) * 2019-06-27 2019-09-27 苏州中科天启遥感科技有限公司 一种光学卫星遥感影像云检测方法
CN111881816A (zh) * 2020-07-27 2020-11-03 安徽省（水利部淮河水利委员会）水利科学研究院（安徽省水利工程质量检测中心站） 一种长时序的河湖围埂养殖区域监测方法
CN112381815A (zh) * 2020-11-26 2021-02-19 浙江易智信息技术有限公司 多维融合遥感大数据去云方法
CN113850139A (zh) * 2021-08-25 2021-12-28 南京林业大学 一种基于多源遥感的森林年际物候监测方法
CN114119617A (zh) * 2021-11-12 2022-03-01 武汉大学 一种多光谱卫星遥感影像的内陆盐湖卤虫带提取方法
CN114140701A (zh) * 2021-12-10 2022-03-04 广西师范大学 一种多尺度多时相遥感影像去云的方法
CN114202698A (zh) * 2021-12-15 2022-03-18 南京大学 联合主被动遥感数据的水体精准提取方法
CN114418970A (zh) * 2021-12-31 2022-04-29 中国科学院空天信息创新研究院 基于卫星遥感的雾霾分布与气溶胶光学厚度检测方法和装置
CN114511786A (zh) * 2022-04-20 2022-05-17 中国石油大学(华东) 融合多时相信息和分通道密集卷积的遥感图像去云方法
CN115131674A (zh) * 2022-06-24 2022-09-30 武汉大学 基于深度低秩网络的多时相光学遥感影像云检测方法
CN115165773A (zh) * 2022-08-01 2022-10-11 上海兰桂骐技术发展股份有限公司 一种基于Google Earth Engine的水稻面积提取方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Cloud and cloud shadow removal of landsat 8 images using Multitemporal Cloud Removal method;Danang Surya Candra等;《Danang Surya Candra》;1-5页 *
多云多雨地区光学与SAR影像融合增强方法研究;于彬;《中国优秀硕士学位论文全文数据库 (基础科学辑)》;第2019年卷(第11期);A008-144页 *

Also Published As

Publication number Publication date
CN116245757A (zh) 2023-06-09

Similar Documents

Publication Publication Date Title
Henriksen et al. Extracting accurate and precise topography from LROC narrow angle camera stereo observations
CN104299228B (zh) 一种基于精确点位预测模型的遥感影像密集匹配方法
CN112488924A (zh) 一种图像超分辨率模型训练方法、重建方法及装置
CN114092835B (zh) 基于不同时空分辨率的归一化植被指数数据时空融合方法
CN107507152A (zh) 一种基于多图像局部插值的遥感图像缺失数据修复方法
CN107966210B (zh) 基于高光谱图像的热红外融合重构方法
Zhu et al. Robust registration of aerial images and LiDAR data using spatial constraints and Gabor structural features
CN108932708A (zh) 基于超分辨增强的星载多光谱遥感图像融合方法
CN116245757B (zh) 多模态数据的多场景通用性遥感影像云修复方法和系统
CN116402942A (zh) 一种融合多尺度图像特征的大规模建筑物三维重建方法
CN114529830A (zh) 基于混合卷积网络遥感图像时空融合方法
CN115471749A (zh) 地外探测无监督学习的多视角多尺度目标识别方法及系统
Toutin et al. Fusion of Radarsat-2 polarimetric images for improved stereo-radargrammetric DEM
CN111368716B (zh) 一种基于多源时空数据的地质灾害灾毁耕地提取方法
Bertone et al. Highly Resolved Topography and Illumination at Mercury’s South Pole from MESSENGER MDIS NAC
US20230089827A1 (en) Method for selecting stereo pairs of aerial or satellite images to generate elevation data
CN107194334B (zh) 基于光流模型的视频卫星影像密集匹配方法及系统
Liu et al. High-spatial-resolution nighttime light dataset acquisition based on volunteered passenger aircraft remote sensing
Bailey et al. Determining large scale sandbar behaviour
Qin et al. A coarse elevation map-based registration method for super-resolution of three-line scanner images
CN114494039A (zh) 一种水下高光谱推扫图像几何校正的方法
CN113029332A (zh) 卫星云图的预测方法、装置和处理器
Sarabandi et al. Building inventory compilation for disaster management: Application of remote sensing and statistical modeling
Andaru et al. Multitemporal UAV photogrammetry for sandbank morphological change analysis: evaluations of camera calibration methods, co-registration strategies, and the reconstructed DSMs
Oshio et al. Generating DTM from DSM Using a Conditional GAN in Built-Up Areas

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant