CN115346004B - 联合时空重建与cuda加速的遥感时间序列数据重建方法 - Google Patents

联合时空重建与cuda加速的遥感时间序列数据重建方法 Download PDF

Info

Publication number
CN115346004B
CN115346004B CN202211272949.8A CN202211272949A CN115346004B CN 115346004 B CN115346004 B CN 115346004B CN 202211272949 A CN202211272949 A CN 202211272949A CN 115346004 B CN115346004 B CN 115346004B
Authority
CN
China
Prior art keywords
time
space
reconstruction
reconstructed
pixel
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
CN202211272949.8A
Other languages
English (en)
Other versions
CN115346004A (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.)
Shenzhen Planning And Natural Resources Data Management Center Shenzhen Spatial Geographic Information Center
Original Assignee
Shenzhen Planning And Natural Resources Data Management Center Shenzhen Spatial Geographic Information Center
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 Shenzhen Planning And Natural Resources Data Management Center Shenzhen Spatial Geographic Information Center filed Critical Shenzhen Planning And Natural Resources Data Management Center Shenzhen Spatial Geographic Information Center
Priority to CN202211272949.8A priority Critical patent/CN115346004B/zh
Publication of CN115346004A publication Critical patent/CN115346004A/zh
Application granted granted Critical
Publication of CN115346004B publication Critical patent/CN115346004B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/50—Information retrieval; Database structures therefor; File system structures therefor of still image data
    • G06F16/58—Retrieval characterised by using metadata, e.g. metadata not derived from the content or metadata generated manually
    • G06F16/5866—Retrieval characterised by using metadata, e.g. metadata not derived from the content or metadata generated manually using information manually generated, e.g. tags, keywords, comments, manually generated location and time information
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/50—Information retrieval; Database structures therefor; File system structures therefor of still image data
    • G06F16/58—Retrieval characterised by using metadata, e.g. metadata not derived from the content or metadata generated manually
    • G06F16/587—Retrieval characterised by using metadata, e.g. metadata not derived from the content or metadata generated manually using geographical or spatial information, e.g. location
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T19/00—Manipulating 3D models or images for computer graphics
    • G06T19/20—Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Library & Information Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Architecture (AREA)
  • Computer Hardware Design (AREA)
  • Image Analysis (AREA)

Abstract

本发明涉及联合时空重建与CUDA加速的遥感时间序列数据重建方法，包括将遥感时间序列数据转换为时空三维向量，利用地物的时空相似性分别从时间和空间维度对缺失像元进行重建，然后利用CUDA来并行实现时间重建、空间重建和结果合并，形成最终重建数据集。本发明的有益效果是：本发明充分利用了遥感时间序列数据的时空相关信息，采用GPU并行计算，从质量和效率两个方面提升重建效果，具有重要的实际应用价值。

Description

联合时空重建与CUDA加速的遥感时间序列数据重建方法
技术领域
本发明涉及光学遥感图像处理技术领域，更确切地说，它涉及联合时空重建与CUDA加速的遥感时间序列数据重建方法。
背景技术
遥感时间序列数据可以提供地物变化的细节信息，被广泛应用在地区或全球环境变化研究中。然而，由于受到云、雾、冰雪和传感器本身故障影响，影像数据往往会出现信息缺失导致时空不连续现象，因此研究遥感影像重建方法，构建完善的遥感时间序列数据，是后续进行长时序遥感应用的重要基础，也是一项意义重大的课题。
目前，遥感数据重建算法主要包括空域重建方法和时域重建法。空域重建方法适用于小面积的重建问题，对于大面积的空间缺失，重建结果会出现过度平滑、模糊或重建不完整等现象。时域重建方法是利用多时相数据作为互补信息源，通过建立多时相之间的数学关系，实现特定区域缺失影像的信息重建，但在地物变化较快的地区应用效果不佳。此外，大多研究数据主要为低分辨率遥感数据，如MODIS，中高分辨率的数据研究较少。中高分辨率遥感数据的体量是低分辨率遥感数据的上百倍甚至上千倍，一般的重建方法应用在中高分辨率遥感数据上效率低下，方法实用性受限。
发明内容
本发明的目的是克服现有遥感时序数据重建方法在时空互补信息利用不足，且运行效率低的缺陷，提供了联合时空重建与CUDA加速的遥感时间序列数据重建方法。
第一方面，提供了联合时空重建与CUDA加速的遥感时间序列数据重建方法，包括：
步骤1、由CPU读取研究区的遥感时间序列数据，并保存在主机（Host）内存中，同时对所述遥感时间序列数据进行预处理，选择要重建的波段；
步骤2、以向量的形式将遥感时间序列数据表示为时空三维向量
Figure 289868DEST_PATH_IMAGE001
，其中I，J和Z分别表示空间经度、空间纬度和时间；
Figure 613533DEST_PATH_IMAGE002
代表一个具体像元的值，缺失像元的值为0；
步骤3、根据GPU显存大小与时空三维向量
Figure 924428DEST_PATH_IMAGE003
，获取待重建时空三维向量
Figure 726162DEST_PATH_IMAGE004
，并将所述待重建时空三维向量
Figure 71693DEST_PATH_IMAGE004
从CPU内存复制到GPU显存中；
步骤4、利用CUDA构建一种空间重建函数，沿着待重建时空三维向量
Figure 300680DEST_PATH_IMAGE004
的空间维度，对每个待重建像元进行并行的空间重建；
步骤5、利用CUDA构建一种时间重建函数，沿着待重建时空三维向量
Figure 98872DEST_PATH_IMAGE004
的时间维度，对每个待重建像元进行并行的时间重建；
步骤6、利用CUDA构建一种合并函数，并行的将待重建时空三维向量V上每个待重建像元的时间重建结果和空间重建结果加权平均，得到待重建时空三维向量
Figure 721875DEST_PATH_IMAGE004
的最终重建结果
Figure 921912DEST_PATH_IMAGE005
，并将
Figure 321800DEST_PATH_IMAGE005
从GPU显存复制到CPU内存，释放显存；
步骤7、根据待重建时空三维向量
Figure 607288DEST_PATH_IMAGE004
的最终重建结果
Figure 750825DEST_PATH_IMAGE005
，获取时空三维向量
Figure 70948DEST_PATH_IMAGE003
的最终重建结果
Figure 641737DEST_PATH_IMAGE006
；
步骤8、检查重建结果，若仍存在缺失，重复步骤3-7，直至完成达到重建要求。
作为优选，步骤1中，所述预处理包括几何校正、大气校正、影像配准和云检测。
作为优选，步骤3包括：
步骤3.1、根据GPU显存大小计算理论上最大处理面积
Figure 414521DEST_PATH_IMAGE007
，表示为：
Figure 96169DEST_PATH_IMAGE008
其中，
Figure 270799DEST_PATH_IMAGE009
表示显存大小，
Figure 543648DEST_PATH_IMAGE010
表示遥感数据每个像元所占内存大小，
Figure 679095DEST_PATH_IMAGE011
表示每个像元的长度；
步骤3.2、将所述最大处理面积
Figure 289068DEST_PATH_IMAGE007
与所述时空三维向量
Figure 193570DEST_PATH_IMAGE003
的空间面积进行比较，若所述最大处理面积
Figure 496375DEST_PATH_IMAGE007
大于所述时空三维向量
Figure 853538DEST_PATH_IMAGE003
的空间面积，则将所述时空三维向量
Figure 267202DEST_PATH_IMAGE003
作为待重建时空三维向量
Figure 26210DEST_PATH_IMAGE004
；若所述最大处理面积
Figure 499917DEST_PATH_IMAGE007
小于所述时空三维向量
Figure 344376DEST_PATH_IMAGE003
的空间面积，则采用格网切分的方式对所述时空三维向量
Figure 296152DEST_PATH_IMAGE003
进行空间维度分割，形成多个子区域的待重建时空三维向量
Figure 175246DEST_PATH_IMAGE004
，所述子区域的面积小于所述最大处理面积
Figure 819854DEST_PATH_IMAGE007
。
作为优选，步骤7中，若所述时空三维向量
Figure 151609DEST_PATH_IMAGE003
未经过空间维度分割，将待重建时空三维向量
Figure 907076DEST_PATH_IMAGE004
的最终重建结果
Figure 375097DEST_PATH_IMAGE005
作为时空三维向量
Figure 190607DEST_PATH_IMAGE003
的最终重建结果
Figure 9658DEST_PATH_IMAGE006
；若所述时空三维向量
Figure 834395DEST_PATH_IMAGE003
经过空间维度分割，将多个子区域的待重建时空三维向量
Figure 422502DEST_PATH_IMAGE004
的最终重建结果
Figure 143333DEST_PATH_IMAGE005
进行合并，得到时空三维向量
Figure 443821DEST_PATH_IMAGE003
的最终重建结果
Figure 541090DEST_PATH_IMAGE012
。
作为优选，步骤4包括：
步骤4.1、设置空间搜索半径，搜集待重建像元的空间邻居像元集合；
步骤4.2、考虑到光谱时序数据也存在缺失值，计算每个待重建像元与空间邻居像元的有效时序的交集，并将交集个数小于有效时序个数阈值的空间邻居像元删除；
步骤4.3、计算剩余空间邻居像元与待重建像元的时序相关性，生成相关性集合；相关系数的计算方式表示为：
Figure 249283DEST_PATH_IMAGE013
其中，
Figure 141016DEST_PATH_IMAGE014
为时序个数，
Figure 669081DEST_PATH_IMAGE015
为某个时序点，
Figure 570040DEST_PATH_IMAGE016
和
Figure 867161DEST_PATH_IMAGE017
分别为待重建像元和空间邻居像元在时序
Figure 195374DEST_PATH_IMAGE015
上的值，
Figure 210734DEST_PATH_IMAGE018
和
Figure 915385DEST_PATH_IMAGE019
分别为待重建像元和空间邻居像元的光谱时序数据的平均值；
步骤4.4、进一步筛选空间邻居像元，删除相关系数低于相关性阈值的空间邻居像元；如果筛选后的空间邻居像元个数为0，则无法重建，保持空值；否则，将相关系数作为相关权重和距离权重相乘，重新计算空间邻居像元的权重，通过加权平均最终得到待重建像元的空间维度估计值。
作为优选，步骤5包括：
步骤5.1、设置时间搜索半径，搜集待重建像元的时间邻居像元集合；
步骤5.2、根据步骤5.1所搜集的时间邻居像元集合计算待重建像元的时间维度估计值，表示为：
Figure 332591DEST_PATH_IMAGE020
其中，
Figure 831706DEST_PATH_IMAGE021
为待重建像元的时间维度估计值，m为时间邻居像元的个数，
Figure 334362DEST_PATH_IMAGE022
为第l个时间邻居像元的观测值，
Figure 577125DEST_PATH_IMAGE023
为当前时间邻居像元和待重建像元的时间间隔，
Figure 848837DEST_PATH_IMAGE024
为平衡参数，取值范围为0到1；当不存在时间邻居像元时，保留空值不变。
作为优选，步骤4和步骤5中，通过构建时空三维向量索引，利用索引的计算提取时间邻居像元与空间邻居像元；待重建像元
Figure 518853DEST_PATH_IMAGE025
的时空邻居集合分别通过以下公式算出：
Figure 508806DEST_PATH_IMAGE026
}
Figure 289680DEST_PATH_IMAGE027
}
上式中，
Figure 681478DEST_PATH_IMAGE028
和
Figure 522395DEST_PATH_IMAGE029
分别为空间邻居集合和时间邻居集合，SR和TR分别为空间搜索半径和时间搜索半径。
第二方面，提供了联合时空重建与CUDA加速的遥感时间序列数据重建装置，用于执行第一方面任一所述遥感时间序列数据重建方法，包括：
预处理模块，用于将研究区的遥感时间序列数据进行预处理，选择要重建的波段；
表示模块，用于以向量的形式将遥感时间序列数据表示为时空三维向量
Figure 999644DEST_PATH_IMAGE030
，其中I，J和Z分别表示空间经度、空间纬度和时间；
Figure 849788DEST_PATH_IMAGE031
代表一个具体像元的值，缺失像元的值为0；
第一获取模块，用于根据GPU显存大小与时空三维向量
Figure 830514DEST_PATH_IMAGE032
，获取待重建时空三维向量
Figure 842332DEST_PATH_IMAGE033
，并将所述待重建时空三维向量
Figure 541298DEST_PATH_IMAGE033
读入GPU显存中；
空间重建模块，用于利用CUDA构建一种空间重建函数，沿着待重建时空三维向量
Figure 195133DEST_PATH_IMAGE033
的空间维度，对每个待重建像元进行并行的空间重建；
时间重建模块，用于利用CUDA构建一种时间重建函数，沿着待重建时空三维向量
Figure 295944DEST_PATH_IMAGE033
的时间维度，对每个待重建像元进行并行的时间重建；
相加模块，用于利用CUDA构建一个合并函数，并行的将待重建时空三维向量V上每个待重建像元的时间重建结果和空间重建结果相加，得到待重建时空三维向量
Figure 478664DEST_PATH_IMAGE033
的最终重建结果
Figure 930505DEST_PATH_IMAGE034
，并将
Figure 856873DEST_PATH_IMAGE034
从显存导出；
第二获取模块，用于根据待重建时空三维向量
Figure 812190DEST_PATH_IMAGE033
的最终重建结果
Figure 634653DEST_PATH_IMAGE034
，获取时空三维向量
Figure 432845DEST_PATH_IMAGE032
的最终重建结果
Figure 32410DEST_PATH_IMAGE035
；
检查模块，用于检查重建结果，若仍存在缺失，第一获取模块、空间重建模块、时间重建模块、相加模块和第二获取模块重复运行，直至完成达到重建要求。
第三方面，提供了一种计算机存储介质，所述计算机存储介质内存储有计算机程序；所述计算机程序在计算机上运行时，使得计算机执行第一方面任一所述遥感时间序列数据重建方法。
本发明的有益效果是：本发明首先将遥感时间序列数据转换为时空三维向量，利用地物的时空相似性分别从时间和空间维度对缺失像元进行重建，然后将重建结果合并形成最终重建数据集。时间维度，采用改进的简单指数平滑函数进行逐像元重建，改善了传统方法因时间间隔大造成的误差。空间维度，采用改进的反距离加权函数进行逐像元重建，改善了传统方法因空间地物类型不同造成的误差。同时，利用CUDA来并行实现时间重建、空间重建和结果合并，提出空间分割策略和时空邻居搜寻方法，大大提高了重建效率。与现有单一时间重建方法和单一空间重建方法相比，本发明充分利用了遥感时间序列数据的时空相关信息，采用GPU并行计算，从质量和效率两个方面提升重建效果，具有重要的实际应用价值。
附图说明
图1为一种联合时空重建与CUDA加速的遥感时间序列数据重建方法的流程示意图；
图2为另一种联合时空重建与CUDA加速的遥感时间序列数据重建方法的流程示意图；
图3为本发明提供的遥感时间序列数据重建方法与对比方法的重建效果图；
图4为联合时空重建与CUDA加速的遥感时间序列数据重建装置的结构示意图。
具体实施方式
下面结合实施例对本发明做进一步描述。下述实施例的说明只是用于帮助理解本发明。应当指出，对于本技术领域的普通人员来说，在不脱离本发明原理的前提下，还可以对本发明进行若干修饰，这些改进和修饰也落入本发明权利要求的保护范围内。
实施例1
CUDA（compute unified device architecture）是英伟达公司提出的统一计算设备架构，可以结合CPU和GPU的优点进行异构并行计算，使得高并发的计算效率大幅提升。因此，结合空间重建和时间重建的各自优势，充分利用CUDA提升算法效率，是实现大区域遥感时间序列数据高效重建的有效途径。
在此基础上，本发明提供了一种联合时空重建与CUDA加速的遥感时间序列数据重建方法，充分利用了地物的空间相关性和时间相关性，分别从时间和空间两个维度设计算法进行影像重建，然后将时空重建结果进行合并生成最终重建结果，同时利用CUDA来优化重建算法，设计空间分割策略和时空邻居搜寻方法，提升大区域影像重建效率。该方法可实现中高分辨率遥感序列数据高质量重建，计算效率高、实用性强。
具体地，本发明提供的联合时空重建与CUDA加速的遥感时间序列数据重建方法，如图1所示，包括：
步骤1、由CPU读取研究区的遥感时间序列数据，并保存在主机内存中，同时对所述遥感时间序列数据进行预处理，选择要重建的波段。
步骤1中，预处理包括几何校正、大气校正、影像配准和云检测，该步骤可提高影像数据质量，并对云缺失区域进行标记。
步骤2、以向量的形式将遥感时间序列数据表示为时空三维向量
Figure 966868DEST_PATH_IMAGE036
，其中I，J和Z分别表示空间经度、空间纬度和时间；
Figure 632336DEST_PATH_IMAGE031
代表一个具体像元的值，缺失像元的值为0。
步骤3、根据GPU显存大小与时空三维向量
Figure 917824DEST_PATH_IMAGE032
，获取待重建时空三维向量
Figure 61360DEST_PATH_IMAGE033
，并将待重建时空三维向量
Figure 115904DEST_PATH_IMAGE033
从CPU内存复制到GPU显存中。
步骤3包括：
步骤3.1、根据GPU显存大小计算理论上最大处理面积
Figure 952273DEST_PATH_IMAGE037
，表示为：
Figure 459477DEST_PATH_IMAGE038
其中，
Figure 406705DEST_PATH_IMAGE039
表示显存大小，
Figure 581334DEST_PATH_IMAGE040
表示遥感数据每个像元所占内存大小，
Figure 323025DEST_PATH_IMAGE041
表示每个像元的长度；除以2是为了让空间重建与时间重建并行，因此占两倍显存；
步骤3.2、将最大处理面积
Figure 583105DEST_PATH_IMAGE037
与时空三维向量
Figure 68444DEST_PATH_IMAGE032
的空间面积进行比较，若最大处理面积
Figure 363160DEST_PATH_IMAGE037
大于时空三维向量
Figure 275752DEST_PATH_IMAGE032
的空间面积，则将时空三维向量
Figure 23128DEST_PATH_IMAGE032
作为待重建时空三维向量
Figure 312158DEST_PATH_IMAGE033
；若最大处理面积
Figure 195800DEST_PATH_IMAGE037
小于时空三维向量
Figure 544873DEST_PATH_IMAGE032
的空间面积，则采用格网切分的方式对时空三维向量
Figure 717229DEST_PATH_IMAGE032
进行空间维度分割，形成多个子区域的待重建时空三维向量
Figure 669004DEST_PATH_IMAGE033
，子区域的面积小于最大处理面积
Figure 548099DEST_PATH_IMAGE037
。本实施例中，假设
Figure 927127DEST_PATH_IMAGE037
大于时空三维向量
Figure 258883DEST_PATH_IMAGE032
的空间面积，因此，不对时空三维向量
Figure 14349DEST_PATH_IMAGE032
进行切割。
步骤4、利用CUDA构建一种空间重建函数，沿着待重建时空三维向量
Figure 13529DEST_PATH_IMAGE033
的空间维度，对每个待重建像元进行并行的空间重建。
步骤4中的空间重建函数是一种改进的反距离加权函数。原始反距离加权函数是利用空间领域内的已知观测数据来对未知数据进行估计，距未知点越近的空间邻居点被赋予较高的空间贡献权重。但反距离加权函数仅考虑两点之间的距离，忽略了地物类型是否一致。为了保证搜集到的空间邻居像元的有效性，本发明假设距离相近的同类地物的光谱变化曲线相似，利用待重建像元的时序光谱特征与邻近空间像元的时序光谱特征的相似度来筛选领域范围内空间邻居像元，从而改进反距离加权算法，具体实现如下：
Figure 829038DEST_PATH_IMAGE043
上述算法中的步骤可以描述为：
步骤4.1、设置空间搜索半径，搜集待重建像元的空间邻居像元集合（neighborS）；
步骤4.2、考虑到光谱时序数据也存在缺失值，计算每个待重建像元与空间邻居像元的有效时序的交集，并将交集个数小于有效时序个数阈值的空间邻居像元删除（算法第5-7行）；
步骤4.3、计算剩余空间邻居像元与待重建像元的时序相关性，生成相关性集合（neighborS_cor）；算法中的相似度采用皮尔森相关系数，假设待重建和空间邻居像元的光谱时序数据分别为
Figure 382510DEST_PATH_IMAGE044
和
Figure 348192DEST_PATH_IMAGE045
，则两者之间的相关系数的计算方式表示为：
Figure 60934DEST_PATH_IMAGE013
其中，
Figure 250606DEST_PATH_IMAGE014
为时序个数，
Figure 556954DEST_PATH_IMAGE015
为某个时序点，
Figure 654223DEST_PATH_IMAGE046
和
Figure 90977DEST_PATH_IMAGE047
分别为待重建像元和空间邻居像元在时序
Figure 248289DEST_PATH_IMAGE015
上的值，
Figure 41933DEST_PATH_IMAGE048
和
Figure 677314DEST_PATH_IMAGE019
分别为待重建像元和空间邻居像元的光谱时序数据的平均值；
步骤4.4、为了保证空间邻居像元的可信度，进一步筛选空间邻居像元，删除相关系数低于相关性阈值（
Figure 240013DEST_PATH_IMAGE049
）的空间邻居像元（第10-14行）；如果筛选后的空间邻居像元个数为0，则无法重建，保持空值；否则，将相关系数作为相关权重和距离权重（D_weight）相乘，重新计算空间邻居像元的权重（neighborS_weight），通过加权平均最终得到待重建像元的空间维度估计值（第15-19行）。
步骤5、利用CUDA构建一种时间重建函数，沿着待重建时空三维向量
Figure 568226DEST_PATH_IMAGE004
的时间维度，对每个待重建像元进行并行的时间重建。
步骤5中提出的时间重建函数是一种改进的简单指数平滑函数。原始简单指数平滑函数是假设时间邻居像元与待重建像元时间距离越近，则其贡献的权重越大，一般采用待重建像元之前的所有有效像元作为时间邻居像元进行建模。但是，采用所有时间邻居像元带来的误差较大，对此，本申请的步骤5包括：
步骤5.1、设置时间搜索半径，搜集待重建像元的时间邻居像元集合；
示例地，在遥感影像重返周期为r时，以前后2r作为时间窗口搜集待重建像元的时间邻居像元；这样在保证地物光谱变化不大的前提下，理论上有4个时间邻居像元，避免采用所有时间邻居像元带来的误差；
步骤5.2、根据步骤5.1所搜集的时间邻居像元计算待重建像元的时间维度估计值，表示为：
Figure 583587DEST_PATH_IMAGE020
其中，
Figure 288238DEST_PATH_IMAGE021
为待重建像元的时间维度估计值，m为时间邻居像元的个数，
Figure 705444DEST_PATH_IMAGE050
为第l个时间邻居像元的观测值，
Figure 204558DEST_PATH_IMAGE051
为当前时间邻居像元和待重建像元的时间间隔，
Figure 707215DEST_PATH_IMAGE024
为平衡参数，取值范围为0到1；当不存在时间邻居像元时，保留空值不变。
此外，步骤4和步骤5中均需要搜寻待重建像元的时间和空间邻居像元，但通过全局搜索和距离比较的方式极为耗时，本发明考虑到时空向量的排列顺序已经隐含了时空邻近关系，因此通过构建时空三维向量索引，可以利用索引的计算快速提取时间邻居像元与空间邻居像元。具体来说，对于一个待重建像元
Figure 684398DEST_PATH_IMAGE025
，
Figure 221690DEST_PATH_IMAGE052
分别为经度、纬度和时间维度上的索引，待重建像元
Figure 891705DEST_PATH_IMAGE025
的时空邻居集合分别通过以下公式算出：
Figure 147237DEST_PATH_IMAGE053
}
Figure 928112DEST_PATH_IMAGE054
}
上式中，
Figure 54331DEST_PATH_IMAGE055
和
Figure DEST_PATH_IMAGE056
分别为空间邻居集合和时间邻居集合，SR和TR分别为空间搜索半径和时间搜索半径。
步骤6、利用CUDA构建一种合并函数，并行的将待重建时空三维向量V上每个待重建像元的时间重建结果和空间重建结果加权平均，得到待重建时空三维向量
Figure 301772DEST_PATH_IMAGE004
的最终重建结果
Figure 903655DEST_PATH_IMAGE005
，并将
Figure 363586DEST_PATH_IMAGE005
从GPU显存复制到CPU内存，释放显存。
步骤7、根据待重建时空三维向量
Figure 104DEST_PATH_IMAGE004
的最终重建结果
Figure 887288DEST_PATH_IMAGE005
，获取时空三维向量
Figure 710888DEST_PATH_IMAGE003
的最终重建结果
Figure 974510DEST_PATH_IMAGE012
。
由于本实施例在步骤3.2中不对时空三维向量
Figure 199955DEST_PATH_IMAGE003
进行切割，因此，待重建时空三维向量
Figure 258041DEST_PATH_IMAGE004
的最终重建结果
Figure 834516DEST_PATH_IMAGE005
可以直接作为时空三维向量
Figure 964146DEST_PATH_IMAGE003
的最终重建结果
Figure 919464DEST_PATH_IMAGE006
。
步骤8、检查重建结果，若仍存在缺失，重复步骤3-7，直至完成达到重建要求。
步骤8中，本申请不对具体的重建要求进行限定，示例地，重建要求为缺失率小于10%。又示例地，重建要求也可以为完成所有缺失像元重建。
实施例2
如图2所示，本申请还提供了另一种联合时空重建与CUDA加速的遥感时间序列数据重建方法。在本实施例中，
Figure 538664DEST_PATH_IMAGE007
小于时空三维向量
Figure 212222DEST_PATH_IMAGE003
的空间面积，因此采用格网切分的方式对时空三维向量
Figure 942280DEST_PATH_IMAGE003
进行空间维度分割，形成多个子区域的待重建时空三维向量
Figure 35262DEST_PATH_IMAGE004
，子区域的面积小于最大处理面积
Figure 825363DEST_PATH_IMAGE007
。需要说明的是，若不对S进行切分，则超过GPU处理范围，无法计算。
此外，在得到多个子区域的待重建时空三维向量
Figure 720638DEST_PATH_IMAGE004
的最终重建结果
Figure 254388DEST_PATH_IMAGE005
后，将其进行合并，得到时空三维向量
Figure 184298DEST_PATH_IMAGE003
的最终重建结果
Figure 879721DEST_PATH_IMAGE012
。
实施例3
在本实施例中结合模拟实验结果对本发明的效果进一步分析：
本模拟实验数据为美国Kingsbury县2018年全年的Landsat-7 EMT+遥感影像，对比实验包括三种常用的重建方法：基于反距离加权法的空间重建法（IDW）、基于简单指数平滑的时间重建法（SES）和基于谐波分析的时间重建方法（HANTS）。图3为Kingsbury县在2018年第151天的原始影像和不同方法的重建影像。可以看出，IDW没有完全重建缺失区域，SES和HANTS方法重建区域与附近地物有明显的“边界”现象（水体区域尤为明显），而本发明很好的重建了原始影像原本的条带缺失和模拟缺失。另外，计算模拟缺失区域的重建像元与真值间的R2，也可以得出本发明的效果最好。
为了验证本发明采用CUDA提高遥感时序数据的重建效率的可靠性，本发明模拟构建5组遥感时间序列数据集，分别计算采用纯CPU和CUDA优化下的算法重建时间。5组遥感时间序列数据集的时间序列相同，均为365，空间大小分别为： 100×100像元、500×500像元、1000×1000像元、2000×2000像元和3000×3000像元。另外，本实验基于Google提供的免费云计算环境Google Colaboratory，具体硬件配置：2.2GHz Intel Xeon(R)处理器，25G内存，NVIDIA Tesla V100显卡，16G显存。对比结果如表1所示，可以看出利用CUDA优化算法来进行遥感时序数据重建相比较纯CPU有明显提升，加速比可达到700以上，效率提升巨大。
表1 纯CPU和CUDA优化下算法耗时比较
实验图像大小 CPU（s） CUDA（s） 加速比
100×100×365 72.12 0.24 300.5
500×500×365 1903.45 2.76 689.66
1000×1000×365 7349.8 10.05 731.32
2000×2000×365 29841.24 41.38 721.15
3000×3000×365 66908.32 94.42 708.62
综上，本发明利用地物在时间和空间的相关性，充分结合遥感时序数据重建中空间重建和时间重建的各自优势，分别在空间重建中提出改进的反距离加权方法，在时间重建中提出改进的简单指数平滑方法，然后将两者重建结果合并作为最终的重建结果，解决了传统重建方法存在的重建不完全和重建结果与现实割裂的缺点。同时，采用CUDA加速重建算法，提出空间分割策略和时空邻居搜寻方法，极大地提高了重建效率，实现了高效的大范围遥感时序数据重建，具有重要的应用价值。

Claims (9)

1.联合时空重建与CUDA加速的遥感时间序列数据重建方法，其特征在于，包括：
步骤1、由CPU读取研究区的遥感时间序列数据，并保存在主机内存中，同时对所述遥感时间序列数据进行预处理，选择要重建的波段；
步骤2、以向量的形式将遥感时间序列数据表示为时空三维向量
Figure DEST_PATH_IMAGE001
，其中I，J和Z分别表示空间经度、空间纬度和时间；
Figure 749615DEST_PATH_IMAGE002
代表一个具体像元的值，缺失像元的值为0；
步骤3、根据GPU显存大小与时空三维向量
Figure DEST_PATH_IMAGE003
，获取待重建时空三维向量
Figure 260231DEST_PATH_IMAGE004
，并将所述待重建时空三维向量
Figure 321859DEST_PATH_IMAGE004
从CPU内存复制到GPU显存中；
步骤4、利用CUDA构建一种空间重建函数，沿着待重建时空三维向量
Figure 44964DEST_PATH_IMAGE004
的空间维度，对每个待重建像元进行并行的空间重建；
步骤5、利用CUDA构建一种时间重建函数，沿着待重建时空三维向量
Figure 390495DEST_PATH_IMAGE004
的时间维度，对每个待重建像元进行并行的时间重建；
步骤6、利用CUDA构建一个合并函数，并行的将待重建时空三维向量V上每个待重建像元的时间重建结果和空间重建结果加权平均，得到待重建时空三维向量
Figure 993383DEST_PATH_IMAGE004
的最终重建结果
Figure DEST_PATH_IMAGE005
，并将
Figure 322733DEST_PATH_IMAGE005
从GPU显存复制到CPU内存，释放显存；
步骤7、根据待重建时空三维向量
Figure 52792DEST_PATH_IMAGE004
的最终重建结果
Figure 3562DEST_PATH_IMAGE005
，获取时空三维向量的最终重建结果
Figure 731346DEST_PATH_IMAGE006
；
步骤8、检查重建结果，若仍存在缺失，重复步骤3-7，直至完成达到重建要求。
2.根据权利要求1所述的联合时空重建与CUDA加速的遥感时间序列数据重建方法，其特征在于，步骤1中，所述预处理包括几何校正、大气校正、影像配准和云检测。
3.根据权利要求1所述的联合时空重建与CUDA加速的遥感时间序列数据重建方法，其特征在于，步骤3包括：
步骤3.1、根据GPU显存大小计算理论上最大处理面积
Figure DEST_PATH_IMAGE007
，表示为：
Figure DEST_PATH_IMAGE009
其中，
Figure 79151DEST_PATH_IMAGE010
表示显存大小，
Figure DEST_PATH_IMAGE011
表示遥感数据每个像元所占内存大小，
Figure 130677DEST_PATH_IMAGE012
表示每个像元的长度；
步骤3.2、将所述最大处理面积
Figure 185221DEST_PATH_IMAGE007
与所述时空三维向量
Figure 411803DEST_PATH_IMAGE003
的空间面积进行比较，若所述最大处理面积
Figure 935319DEST_PATH_IMAGE007
大于所述时空三维向量
Figure 741601DEST_PATH_IMAGE003
的空间面积，则将所述时空三维向量
Figure 119493DEST_PATH_IMAGE003
作为待重建时空三维向量
Figure 985818DEST_PATH_IMAGE004
；若所述最大处理面积
Figure 245898DEST_PATH_IMAGE007
小于所述时空三维向量
Figure 370717DEST_PATH_IMAGE003
的空间面积，则采用格网切分的方式对所述时空三维向量
Figure 399853DEST_PATH_IMAGE003
进行空间维度分割，形成多个子区域的待重建时空三维向量
Figure 905921DEST_PATH_IMAGE004
，所述子区域的面积小于所述最大处理面积
Figure 387718DEST_PATH_IMAGE007
。
4.根据权利要求3所述的联合时空重建与CUDA加速的遥感时间序列数据重建方法，其特征在于，步骤7中，若所述时空三维向量
Figure 801382DEST_PATH_IMAGE003
未经过空间维度分割，将待重建时空三维向量
Figure 435756DEST_PATH_IMAGE004
的最终重建结果
Figure 175042DEST_PATH_IMAGE005
作为时空三维向量
Figure 144135DEST_PATH_IMAGE003
的最终重建结果
Figure 95911DEST_PATH_IMAGE006
；若所述时空三维向量
Figure 302901DEST_PATH_IMAGE003
经过空间维度分割，将多个子区域的待重建时空三维向量
Figure 199706DEST_PATH_IMAGE004
的最终重建结果
Figure 656096DEST_PATH_IMAGE005
进行合并，得到时空三维向量
Figure 411562DEST_PATH_IMAGE003
的最终重建结果
Figure 269797DEST_PATH_IMAGE006
。
5.根据权利要求1所述的联合时空重建与CUDA加速的遥感时间序列数据重建方法，其特征在于，步骤4包括：
步骤4.1、设置空间搜索半径，搜集待重建像元的空间邻居像元集合；
步骤4.2、考虑到光谱时序数据也存在缺失值，计算每个待重建像元与空间邻居像元的有效时序的交集，并将交集个数小于有效时序个数阈值的空间邻居像元删除；
步骤4.3、计算剩余空间邻居像元与待重建像元的时序相关性，生成相关性集合；相关系数的计算方式表示为：
Figure DEST_PATH_IMAGE013
其中，
Figure 367197DEST_PATH_IMAGE014
为时序个数，
Figure DEST_PATH_IMAGE015
为某个时序点，
Figure 842040DEST_PATH_IMAGE016
和
Figure DEST_PATH_IMAGE017
分别为待重建像元和空间邻居像元在时序
Figure 181624DEST_PATH_IMAGE015
上的值，
Figure 894365DEST_PATH_IMAGE018
和
Figure DEST_PATH_IMAGE019
分别为待重建像元和空间邻居像元的光谱时序数据的平均值；
步骤4.4、进一步筛选空间邻居像元，删除相关系数低于相关性阈值的空间邻居像元；如果筛选后的空间邻居像元个数为0，则无法重建，保持空值；否则，将相关系数作为相关权重和距离权重相乘，重新计算空间邻居像元的权重，通过加权平均最终得到待重建像元的空间维度估计值。
6.根据权利要求1所述的联合时空重建与CUDA加速的遥感时间序列数据重建方法，其特征在于，步骤5包括：
步骤5.1、设置时间搜索半径，搜集待重建像元的时间邻居像元集合；
步骤5.2、根据步骤5.1所搜集的时间邻居像元集合计算待重建像元的时间维度估计值，表示为：
Figure DEST_PATH_IMAGE021
其中，
Figure 428245DEST_PATH_IMAGE022
为待重建像元的时间维度估计值，m为时间邻居像元的个数，
Figure DEST_PATH_IMAGE023
为第l个时间邻居像元的观测值，
Figure 390385DEST_PATH_IMAGE024
为当前时间邻居像元和待重建像元的时间间隔，
Figure DEST_PATH_IMAGE025
为平衡参数，取值范围为0到1；当不存在时间邻居像元时，保留空值不变。
7.根据权利要求1所述的联合时空重建与CUDA加速的遥感时间序列数据重建方法，其特征在于，步骤4和步骤5中，通过构建时空三维向量索引，利用索引的计算提取时间邻居像元与空间邻居像元；待重建像元
Figure 5431DEST_PATH_IMAGE026
的时空邻居集合分别通过以下公式算出：
Figure DEST_PATH_IMAGE027
}
Figure 369416DEST_PATH_IMAGE028
}
上式中，
Figure DEST_PATH_IMAGE029
和
Figure 543039DEST_PATH_IMAGE030
分别为空间邻居集合和时间邻居集合，SR和TR分别为空间搜索半径和时间搜索半径。
8.联合时空重建与CUDA加速的遥感时间序列数据重建装置，其特征在于，用于执行权利要求1至7任一所述遥感时间序列数据重建方法，包括：
预处理模块，用于将研究区的遥感时间序列数据进行预处理，选择要重建的波段；
表示模块，用于以向量的形式将遥感时间序列数据表示为时空三维向量
Figure 461317DEST_PATH_IMAGE001
，其中I，J和Z分别表示空间经度、空间纬度和时间；
Figure 362277DEST_PATH_IMAGE002
代表一个具体像元的值，缺失像元的值为0；
第一获取模块，用于根据GPU显存大小与时空三维向量
Figure 298878DEST_PATH_IMAGE003
，获取待重建时空三维向量
Figure 627091DEST_PATH_IMAGE004
，并将所述待重建时空三维向量
Figure 970347DEST_PATH_IMAGE004
读入GPU显存中；
空间重建模块，用于利用CUDA构建一种空间重建函数，沿着待重建时空三维向量
Figure 674998DEST_PATH_IMAGE004
的空间维度，对每个待重建像元进行并行的空间重建；
时间重建模块，用于利用CUDA构建一种时间重建函数，沿着待重建时空三维向量
Figure 216838DEST_PATH_IMAGE004
的时间维度，对每个待重建像元进行并行的时间重建；
相加模块，用于利用CUDA构建一个合并函数，并行的将待重建时空三维向量V上每个待重建像元的时间重建结果和空间重建结果相加，得到待重建时空三维向量
Figure 732264DEST_PATH_IMAGE004
的最终重建结果
Figure 359555DEST_PATH_IMAGE005
，并将
Figure 602317DEST_PATH_IMAGE005
从显存导出；
第二获取模块，用于根据待重建时空三维向量
Figure 201926DEST_PATH_IMAGE004
的最终重建结果
Figure 871941DEST_PATH_IMAGE005
，获取时空三维向量
Figure 250444DEST_PATH_IMAGE003
的最终重建结果
Figure 31318DEST_PATH_IMAGE006
；
检查模块，用于检查重建结果，若仍存在缺失，第一获取模块、空间重建模块、时间重建模块、相加模块和第二获取模块重复运行，直至完成达到重建要求。
9.一种计算机存储介质，其特征在于，所述计算机存储介质内存储有计算机程序；所述计算机程序在计算机上运行时，使得计算机执行权利要求1至7任一所述遥感时间序列数据重建方法。
CN202211272949.8A 2022-10-18 2022-10-18 联合时空重建与cuda加速的遥感时间序列数据重建方法 Active CN115346004B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211272949.8A CN115346004B (zh) 2022-10-18 2022-10-18 联合时空重建与cuda加速的遥感时间序列数据重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211272949.8A CN115346004B (zh) 2022-10-18 2022-10-18 联合时空重建与cuda加速的遥感时间序列数据重建方法

Publications (2)

Publication Number Publication Date
CN115346004A CN115346004A (zh) 2022-11-15
CN115346004B true CN115346004B (zh) 2023-01-31

Family

ID=83957621

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211272949.8A Active CN115346004B (zh) 2022-10-18 2022-10-18 联合时空重建与cuda加速的遥感时间序列数据重建方法

Country Status (1)

Country Link
CN (1) CN115346004B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103902802A (zh) * 2012-12-29 2014-07-02 中国科学院深圳先进技术研究院 一种顾及空间信息的植被指数时间序列数据重建方法
CN109902259A (zh) * 2019-02-25 2019-06-18 中国科学院地理科学与资源研究所 一种轻量级的缺失时空数据的重构方法
CN111932457A (zh) * 2020-08-06 2020-11-13 北方工业大学 遥感影像高时空融合处理算法及装置

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8379947B2 (en) * 2010-05-28 2013-02-19 International Business Machines Corporation Spatio-temporal image reconstruction using sparse regression and secondary information
JP6032467B2 (ja) * 2012-06-18 2016-11-30 株式会社日立製作所 時空間データ管理システム、時空間データ管理方法、及びそのプログラム
CN106934783B (zh) * 2017-03-02 2020-01-31 宁波大学 一种高频次遥感时间序列数据的时域重建方法
CN113112591B (zh) * 2021-04-15 2022-08-26 宁波甬矩空间信息技术有限公司 基于耦合稀疏张量分解的多时相遥感影像时空谱融合方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103902802A (zh) * 2012-12-29 2014-07-02 中国科学院深圳先进技术研究院 一种顾及空间信息的植被指数时间序列数据重建方法
CN109902259A (zh) * 2019-02-25 2019-06-18 中国科学院地理科学与资源研究所 一种轻量级的缺失时空数据的重构方法
CN111932457A (zh) * 2020-08-06 2020-11-13 北方工业大学 遥感影像高时空融合处理算法及装置

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Missing Data Reconstruction in Remote Sensing Image With a Unified Spatial–Temporal–Spectral Deep Convolutional Neural Network;Qiang Zhang 等;《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》;20181231;第56卷(第8期);第4274-4288页 *

Also Published As

Publication number Publication date
CN115346004A (zh) 2022-11-15

Similar Documents

Publication Publication Date Title
Uhrig et al. Sparsity invariant cnns
Terriberry et al. GPU accelerating speeded-up robust features
CN101976435B (zh) 基于对偶约束的联合学习超分辨方法
CN103279939B (zh) 一种图像拼接处理系统
CN103034982B (zh) 一种基于变焦视频序列的图像超分辨率重建方法
CN102831598B (zh) 多分辨率NMF和Treelet融合的遥感图像变化检测方法
CN103871039B (zh) 一种sar图像变化检测差异图生成方法
Gou et al. Remote sensing image super-resolution reconstruction based on nonlocal pairwise dictionaries and double regularization
CN107516322A (zh) 一种基于对数极空间的图像物体大小和旋转估计计算方法
Wang et al. Dilated projection correction network based on autoencoder for hyperspectral image super-resolution
Dong et al. Real-world remote sensing image super-resolution via a practical degradation model and a kernel-aware network
Li et al. HS2P: Hierarchical spectral and structure-preserving fusion network for multimodal remote sensing image cloud and shadow removal
Donate et al. Improved reconstruction of images distorted by water waves
Fu et al. Three-dimensional singular spectrum analysis for precise land cover classification from UAV-borne hyperspectral benchmark datasets
Dong et al. A unified approach of multitemporal SAR data filtering through adaptive estimation of complex covariance matrix
Guo et al. A flexible object-level processing strategy to enhance the weight function-based spatiotemporal fusion method
Wang et al. Group shuffle and spectral-spatial fusion for hyperspectral image super-resolution
Liu et al. SACTNet: Spatial attention context transformation network for cloud removal
CN115346004B (zh) 联合时空重建与cuda加速的遥感时间序列数据重建方法
CN111696167A (zh) 自范例学习引导的单张影像超分辨率重构方法
Guo et al. Blind single-image-based thin cloud removal using a cloud perception integrated fast Fourier convolutional network
Jung et al. Intensity-guided edge-preserving depth upsampling through weighted L0 gradient minimization
CN106056551A (zh) 基于局部相似样例学习的稀疏去噪方法
Jenicka et al. Distributed texture-based land cover classification algorithm using hidden Markov model for multispectral data
Luo et al. Real-time pedestrian detection method based on improved YOLOv3

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