CN110555818B - 一种卫星图像序列云区修补方法和装置 - Google Patents
一种卫星图像序列云区修补方法和装置 Download PDFInfo
- Publication number
- CN110555818B CN110555818B CN201910857550.8A CN201910857550A CN110555818B CN 110555818 B CN110555818 B CN 110555818B CN 201910857550 A CN201910857550 A CN 201910857550A CN 110555818 B CN110555818 B CN 110555818B
- Authority
- CN
- China
- Prior art keywords
- cloud
- image
- area
- images
- multispectral
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 57
- 230000008439 repair process Effects 0.000 claims abstract description 85
- 230000004927 fusion Effects 0.000 claims abstract description 44
- 238000004364 calculation method Methods 0.000 claims description 25
- 238000000605 extraction Methods 0.000 claims description 13
- 238000005429 filling process Methods 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 description 23
- 238000012545 processing Methods 0.000 description 20
- 238000002310 reflectometry Methods 0.000 description 16
- 238000001514 detection method Methods 0.000 description 13
- 238000010586 diagram Methods 0.000 description 13
- 230000008569 process Effects 0.000 description 8
- 230000006870 function Effects 0.000 description 6
- 230000008859 change Effects 0.000 description 5
- 238000001914 filtration Methods 0.000 description 5
- 238000012163 sequencing technique Methods 0.000 description 5
- 238000004891 communication Methods 0.000 description 4
- 238000004590 computer program Methods 0.000 description 4
- 238000012937 correction Methods 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000003384 imaging method Methods 0.000 description 3
- 230000003993 interaction Effects 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000000877 morphologic effect Effects 0.000 description 3
- 238000007500 overflow downdraw method Methods 0.000 description 3
- 230000005855 radiation Effects 0.000 description 3
- 230000009467 reduction Effects 0.000 description 3
- 230000003595 spectral effect Effects 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 230000000007 visual effect Effects 0.000 description 3
- 238000012952 Resampling Methods 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- VUFNLQXQSDUXKB-DOFZRALJSA-N 2-[4-[4-[bis(2-chloroethyl)amino]phenyl]butanoyloxy]ethyl (5z,8z,11z,14z)-icosa-5,8,11,14-tetraenoate Chemical compound CCCCC\C=C/C\C=C/C\C=C/C\C=C/CCCC(=O)OCCOC(=O)CCCC1=CC=C(N(CCCl)CCCl)C=C1 VUFNLQXQSDUXKB-DOFZRALJSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 230000006837 decompression Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000007730 finishing process Methods 0.000 description 1
- 239000012634 fragment Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000011065 in-situ storage Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000001151 other effect Effects 0.000 description 1
- 230000001953 sensory effect Effects 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/22—Matching criteria, e.g. proximity measures
-
- 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
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4007—Scaling of whole images or parts thereof, e.g. expanding or contracting based on interpolation, e.g. bilinear interpolation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/90—Dynamic range modification of images or parts thereof
- G06T5/94—Dynamic range modification of images or parts thereof based on local image properties, e.g. for local contrast enhancement
-
- 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/12—Edge-based segmentation
-
- 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/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/40—Extraction of image or video features
- G06V10/58—Extraction of image or video features relating to hyperspectral data
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
- G06V10/761—Proximity, similarity or dissimilarity measures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/13—Satellite images
-
- 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/10016—Video; Image sequence
-
- 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
- G06T2207/10036—Multispectral image; Hyperspectral image
-
- 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/20021—Dividing image into blocks, subimages or windows
-
- 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/20212—Image combination
- G06T2207/20221—Image fusion; Image merging
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Multimedia (AREA)
- Evolutionary Computation (AREA)
- Astronomy & Astrophysics (AREA)
- Artificial Intelligence (AREA)
- Remote Sensing (AREA)
- Data Mining & Analysis (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Computing Systems (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- General Engineering & Computer Science (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Image Processing (AREA)
Abstract
本申请实施例提出一种卫星图像序列云区修补方法和装置。实现方案为：从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像；根据各单波段图像之间的相似度，计算修补多光谱图像所需要的参考图像及修补顺序；根据单波段图像对应的云掩膜图像，确定多光谱图像的修补区域；根据修补顺序，将多光谱图像的修补区域的像素值，替换为参考图像的非云区域的像素值，得到替换修补图像；记录修补区域中每张多光谱图像对应的各参考图像的编号值，得到替换参考掩膜图像；通过泊松融合对目标区域中的各替换修补图像和各替换参考掩膜图像进行云修补。本申请提高了卫星图像序列云区修补的效率和准确度。
Description
技术领域
本申请涉及遥感图像处理技术，尤其涉及一种卫星图像序列云区修补方法和装置。
背景技术
云在遥感卫星图像中具有较高的辐射亮度。卫星影像在成像过程中由于受到云层遮挡，导致原地物光谱失真，影响遥感数据产品的生产与应用。云区的检测与修补一直是遥感图像处理中的重要课题。早期云检测与修补技术由于受遥感数据来源限制而多基于单幅或少量图像，云检测与云区修补的精度有限。随着遥感技术的发展与投入运行的遥感卫星的增加，目前可以获取的卫星数据已有海量的规模。利用现有海量历史遥感数据，基于定量遥感理论反演云与云下阴影的分布，基于同区域海量历史影像插补云与云下阴影区域，获取高质量无云遥感图像成为目前遥感数据预处理的研究热点。
发明内容
本申请实施例提供一种卫星图像序列云区修补方法和装置，以解决现有技术中的一个或多个技术问题。
第一方面，本申请实施例提供了一种卫星图像序列云区修补方法，包括：
从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像；
根据各单波段图像之间的相似度，计算修补多光谱图像所需要的至少一个参考图像，以及各参考图像的修补顺序，参考图像选自于多张多光谱图像中一张或多张；
根据单波段图像对应的云掩膜图像，确定多光谱图像的修补区域；
根据修补顺序，将多光谱图像的修补区域的像素值，替换为参考图像的非云区域的像素值，得到替换修补图像；
记录修补区域中每张多光谱图像对应的各参考图像的编号值，得到替换参考掩膜图像；
通过泊松融合对目标区域中的各替换修补图像和各替换参考掩膜图像进行云修补。
在一种实施方式中，从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像，包括：
从多光谱图像中识别云区域和云下阴影区域，生成待选云掩膜图像；
将待选云掩膜图像中小于第一预设像素值的云区域和云下阴影区域重新标记为清晰地表，将小于第二预设像素值的清晰地表重新标记为云区域；
将重新标记的云区域以半径不超过第一扩展像素的尺度膨胀，并将重新标记的云下阴影区域以半径不超过第二扩展像素的尺度膨胀，获得云掩膜图像。
在一种实施方式中，根据各单波段图像之间的相似度，计算修补多光谱图像所需要的至少一个参考图像，以及各参考图像的修补顺序，包括：
对多张多光谱图像和多张云掩膜图像，按照预设波段采用最近邻插值法，生成各单波段图像的缩略图和各云掩膜图像的缩略图；
根据各单波段图像的缩略图，计算各单波段图像之间的相似度；
根据各参考图像与单波段图像的相似度排序，依次使用参考图像的缩略图中的非云区域像素，填充单波段图像对应的多光谱图像中的标记区域，标记区域由单波段图像对应的云掩膜图像的缩略图标记；
在标记区域被填充完毕或填充区域不再减少的情况下，记录填充过程中用到的参考图像和对应的修补顺序。
在一种实施方式中，通过泊松融合对目标区域中的各替换修补图像和各替换参考掩膜图像进行云修补，包括：
将多光谱图像修补区域外的引导向量值为0；
将多光谱图像的边界划分为修补区域边界和修补区域内部边界；
对于修补区域内部边界位置的引导向量值，采用来自相同参考图像的邻近像素值。
在一种实施方式中，根据各单波段图像的缩略图，计算各单波段图像之间的相似度，包括：
通过结构相似性指数、图像获取时间、云盖量计算相似度，计算公式如下：
其中，a、b、c为加权系数，Sij为单波段图像的缩略图Ii与Ij的相似度，SSIM(Ii,Ij)为缩略图Ii与Ij的结构相似性指数，Tij为缩略图Ii与Ij的获取时间间隔，Ci与Cj分别为缩略图Ii与Ij的云盖量，Cij为缩略图Ii与Ij的公共云区域与云下阴影区域。
在一种实施方式中，结构相似性指数的计算公式为:
其中，SSIM(Ii,Ij)为缩略图Ii与Ij的结构相似性指数，μi与μj为缩略图Ii与Ij的对应云掩膜区域值都为1的像素位置的均值，σi与σj为缩略图Ii与Ij的对应云掩膜区域值都为1的像素位置的标准差，σij为缩略图Ii与Ij的对应云掩膜区域值都为1的像素位置的协方差，云掩膜区域值包括云区域值和云下阴影区域值。
第二方面，本申请实施例提供了一种卫星图像序列云区修补装置，包括：
云掩膜提取整饰模块，用于从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像；
相似度计算模块，用于根据各单波段图像之间的相似度，计算修补多光谱图像所需要的至少一个参考图像，以及各参考图像的修补顺序，参考图像选自于多张多光谱图像中一张或多张；
参考数据选定模块，用于根据单波段图像对应的云掩膜图像，确定多光谱图像的修补区域；
直接替换修补模块，用于根据修补顺序，将多光谱图像的修补区域的像素值，替换为参考图像的非云区域的像素值，得到替换修补图像；
替换记录模块，用于记录修补区域中每张多光谱图像对应的各参考图像的编号值，得到替换参考掩膜图像；
泊松融合修补模块，用于通过泊松融合对替换修补图像和替换参考掩膜图像进行云修补。
在一种实施方式中，云掩膜提取整饰模块，包括：
初始单元，用于从多光谱图像中识别云区域和云下阴影区域，生成待选云掩膜图像；
标记单元，用于将待选云掩膜图像中小于第一预设像素值的云区域和云下阴影区域重新标记为清晰地表，将小于第二预设像素值的清晰地表重新标记为云区域；
膨胀单元，用于将重新标记的云区域以半径不超过第一扩展像素的尺度膨胀，并将重新标记的云下阴影区域以半径不超过第二扩展像素的尺度膨胀，获得云掩膜图像。
在一种实施方式中，相似度计算模块，包括：
插值单元，用于对多张多光谱图像和多张云掩膜图像，按照预设波段采用最近邻插值法，生成各单波段图像的缩略图和各云掩膜图像的缩略图；
计算单元，用于根据各单波段图像的缩略图，计算各单波段图像之间的相似度；
填充单元，用根据各参考图像与单波段图像的相似度排序，依次使用参考图像的缩略图中的非云区域像素，填充单波段图像对应的多光谱图像中的标记区域，标记区域由单波段图像对应的云掩膜图像的缩略图标记；
记录单元，用于在标记区域被填充完毕或填充区域不再减少的情况下，记录填充过程中用到的参考图像和对应的修补顺序。
在一种实施方式中，泊松融合修补模块，包括：
设置单元，用于将多光谱图像修补区域外的引导向量值为0。
划分单元，用于将多光谱图像的边界划分为修补区域边界和修补区域内部边界；
融合单元，用于对于修补区域内部边界位置的引导向量值，采用来自相同参考图像的邻近像素值。
上述申请中的一个实施例具有如下优点或有益效果：一方面，提高了云检测与修补的精度与适用性，另一方面，在结果精度与运算复杂度上更能满足海量数据自动处理的需求。
上述可选方式所具有的其他效果将在下文中结合具体实施例加以说明。
附图说明
附图用于更好地理解本方案，不构成对本申请的限定。其中：
图1示出根据本申请一实施例的卫星图像序列云区修补方法流程图。
图2示出根据本申请一实施例的遥感卫星多光谱图像自动云修补流程图。
图3示出根据本申请一实施例的云掩膜整饰前后比较示例图。
图4示出根据本申请一实施例的泊松融合法云修补的算法原理示意图。
图5示出根据本申请一实施例的直接替换法云修补与泊松融合法云修补结果示例图。
图6示出根据本申请一实施例的同区域多图像云修补前后及替换修补参考掩膜示例图。
图7示出根据本申请一实施例的卫星图像序列云区修补装置结构框图。
图8示出根据本申请一实施例的云掩膜提取整饰模块结构框图。
图9示出根据本申请一实施例的相似度计算模块结构框图。
图10示出根据本申请一实施例的泊松融合修补模块结构框图。
具体实施方式
以下结合附图对本申请的示范性实施例做出说明，其中包括本申请实施例的各种细节以助于理解，应当将它们认为仅仅是示范性的。因此，本领域普通技术人员应当认识到，可以对这里描述的实施例做出各种改变和修改，而不会背离本申请的范围和精神。同样，为了清楚和简明，以下的描述中省略了对公知功能和结构的描述。
本申请的主要思路为：对于遥感卫星在同一地理区域(目标区域)获取的不同时间的多幅多光谱图像数据，假设地表覆盖类型没有发生大的变化，而导致变化的主要因素是运动的云与云下阴影。利用地表反射率产品包含的质量波段实现云与云下阴影掩膜的自动生成，采用缩略图的方式快速获取图像之间的相似性排序，进而确定修补每幅多光谱图像所采用的多幅参考图像，然后对每幅多光谱图像中云掩膜所标记的云区域与云下阴影区域，按照修补排序将参考图像对应像素位置直接替换到待修补的多光谱图像中，得到每幅多光谱图像的直接替换修补图像，同时记录每个替换像素对应的参考图像的编号，得到替换参考掩膜图像，最后对直接替换修补图像与替换参考掩膜图像使用泊松融合算法消除修补区域的光谱差异，得到每幅多光谱图像无缝的云与云下阴影修补结果。
其中的遥感卫星多光谱图像数据，其典型代表为Landsat卫星，包含Landsat4/5/7/8，传感器有TM、ETM+、OLI，空间分辨率不低于30米的全色、多光谱图像，不包含热红外图像。本申请数据源不限定特指Landsat数据源，可适用目前多数遥感卫星多光谱图像产品数据，要求数据经过系统几何校正且数据之间平均配准误差不超过1个像素，要求包含用于提取云掩膜的质量波段，数据数据处理级别使用地表反射率数据，也可以使用原始DN值数据，处理时要求同一区域至少有两幅图像。
图1示出根据本申请一实施例的卫星图像序列云区修补流程图。如图1所示，该数据存储方法包括：
步骤S11、从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像。
在一种实施方式中，从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像，包括：
从多光谱图像中识别云区域和云下阴影区域，生成待选云掩膜图像。
在一种示例中，该步骤主要进行云掩膜提取，对所有待处理的多光谱图像，按照区域归类，逐个区域处理，对于每个待处理图像，根据多光谱图像数据产品包含的质量波段中标记的云与云下阴影区域，生成单波段云掩膜图像。
其中，按照区域归类，是将相同的航带号的所有图像划归到同一区域，比如Landsat卫星航带编号采用WRS(world wide reference system)参考系，结合星下点成像特性而形成的固定地面参考网格，WRS网格的二维坐标采用Path和Row进行标识，具有相同Path与Row的待处理数据被划归为同一区域。中国的遥感卫星多光谱图像数据命名中同样包含与Path和Row标识类似的命名规则。所述质量波段，为多光谱图像产品数据自带数据，通常以不同的标记记录云、云下阴影、雪、水体等标记信息。所述单波段云掩膜图像，是从将质量波段数据重新简化标记，具体是将背景标记0值，清晰地表标记1值，厚云与卷云-高标记2值，云下阴影标记3值，使用8位记录。
步骤S12、根据各单波段图像之间的相似度，计算修补多光谱图像所需要的至少一个参考图像，以及各参考图像的修补顺序，参考图像选自于多张多光谱图像中一张或多张。
在一种实施方式中，将待选云掩膜图像中小于第一预设像素值的云区域和云下阴影区域重新标记为清晰地表，将小于第二预设像素值的清晰地表重新标记为云区域；
将重新标记的云区域以半径不超过第一扩展像素的尺度膨胀，并将重新标记的云下阴影区域以半径不超过第二扩展像素的尺度膨胀，获得云掩膜图像。
在一种示例中，该步骤主要进行云掩膜整饰，利用形态学算法扩大云与云下阴影区域的范围。
利用形态学算法扩大云与云下阴影区域的范围，包含两步处理，第一步减小离散化水平，具体方法是将云掩膜中小于4个像素(第一预设像素值)的云与云下阴影重标记为清晰地表值1，将小于4个像素(第二预设像素值)的清晰地表重标记为云值2，第二步是膨胀云与云下阴影区域，将标记2值的云区域以半径不超过5个像素(第一扩展像素)的尺度膨胀，标记3值的云下阴影区域以半径不超过10个像素(第二扩展像素)的尺度膨胀。
步骤S13、根据单波段图像对应的云掩膜图像，确定多光谱图像的修补区域。
在一种实施方式中，根据各单波段图像之间的相似度，计算修补多光谱图像所需要的至少一个参考图像，以及各参考图像的修补顺序，包括：
对多张多光谱图像和多张云掩膜图像，按照预设波段采用最近邻插值法，生成各单波段图像的缩略图和各云掩膜图像的缩略图；
根据各单波段图像的缩略图，计算各单波段图像之间的相似度。
在一种实施方式中，通过结构相似性指数、图像获取时间、云盖量计算该相似度，计算公式如下：
其中，a、b、c为加权系数，Sij为该单波段图像的缩略图Ii与Ij的相似度，SSIM(Ii,Ij)为缩略图Ii与Ij的结构相似性指数，Tij为缩略图Ii与Ij的获取时间间隔，Ci与Cj分别为缩略图Ii与Ij的云盖量，Cij为缩略图Ii与Ij的公共云区域与云下阴影区域。
在一种实施方式中，结构相似性指数的计算公式为:
其中，SSIM(Ii,Ij)为缩略图Ii与Ij的结构相似性指数，μi与μj为缩略图Ii与Ij的对应云掩膜区域值都为1的像素位置的均值，σi与σj为缩略图Ii与Ij的对应云掩膜区域值都为1的像素位置的标准差，σij为缩略图Ii与Ij的对应云掩膜区域值都为1的像素位置的协方差，该云掩膜区域值包括云区域值和云下阴影区域值。
在一种示例中，该步骤主要通过计算相似性来确定修补区域，对同区域的所有多光谱图像，以及对应的云掩膜图像，生成缩略图，利用缩略图快速计算所有多光谱图像之间的相似性。
在缩略图的生成过程中，对同区域的所有多光谱图像，仅选取受大气影响明显的蓝光(0.450–0.515μm)单个波段按照最近邻插值方法生成缩略图，对应的云掩膜图像按照最近邻插值方法生成缩略图，缩小比例为原始尺寸的四分之一。所述利用缩略图快速计算所有多光谱图像之间的相似性，是使用结构相似性指数SSIM(Structural SimilarityIndex，结构相似性)，统计两幅缩略图中公共清晰地表区域，即对应云掩膜区域值都为1的像素位置。所述的计算相似性，是对于同区域的每一幅多光谱图像，依据SSIM、图像获取时间、云盖量获取其它图像与该图像的相似性。
步骤S14、根据修补顺序，将多光谱图像的修补区域的像素值，替换为参考图像的非云区域的像素值，得到替换修补图像。
在一种实施方式中，根据各参考图像与单波段图像的相似度排序，依次使用参考图像的缩略图中的非云区域像素，填充单波段图像对应的多光谱图像中的标记区域，标记区域由单波段图像对应的云掩膜图像的缩略图标记；
在标记区域被填充完毕或填充区域不再减少的情况下，记录填充过程中用到的参考图像和对应的修补顺序。
在一种示例中，该步骤主要确定云修补参考数据，对同区域的每一个多光谱图像，利用图像缩略图、云掩膜缩略图及相似性排序的结果，计算修补每一幅图像所需要的参考图像数量与修补顺序，得到云修补参考查找表。
确定云修补参考数据，是对同区域的每一个多光谱图像的缩略图、对其云掩膜缩略图所标记的云与云下阴影区域，依次使用相似性高的图像缩略图非云区域像素填充，直到填充完毕或填充区域不再减少为止，记录用到的参考图像数量、名称、修补顺序。所述的云修补参考查找表，记录同区域的每一个多光谱图像云修补所需要的参考图像数量、名称、修补顺序，保存为外部文本文件。
步骤S15、记录修补区域中每张多光谱图像对应的各参考图像的编号值，得到替换参考掩膜图像。
在一种示例中，该步骤主要进行直接替换法云修补，对同区域的每一个多光谱图像，对于其云掩膜所标记的云与云下阴影待修补区域，根据云修补参考查找表所确定的参考图像修补顺序，将参考图像非云区域像素值逐波段替换到待修补区域，得到替换修补图像，同时记录修补区域每个像素对应参考图像的编号值，得到替换参考掩膜图像。
直接替换法云修补，具体是对同区域的每一个多光谱图像，对于其云掩膜所标记的云与云下阴影待修补区域，根据云修补参考查找表所确定的参考图像修补顺序，将参考图像非云区域像素值逐波段替换到待修补区域，得到替换修补图像。所述替换参考掩膜图像，为单波段8位图像，其中像素0值对应云掩膜图像标记的背景值与清晰地表值，非0值对应云掩膜图像中的云与云下阴影区域，非0值取值于修补云区所用参考图像的编号值，比如i值表示该位置像素云修补时来自第i幅参考图像，i∈[1,n]，n为参考图像数量。
步骤S16、通过泊松融合对目标区域中的各替换修补图像和各替换参考掩膜图像进行云修补。
在一种示例中，该步骤主要进行泊松融合法云修补，对同区域的每一个多光谱替换修补图像与替换参考掩膜图像，利用泊松融合完成云修补，得到最终多光谱云修补图像。
泊松融合法云修补，与目前无人机图像无缝拼接所采用的泊松融合在求解过程是一致的，差异在于泊松融合求解所使用的引导向量域的计算，修补区域外的引导向量值为0，修补区域内引导向量值按照一般泊松融合的引导向量运算，将边界划分为修补区域边界与修补区域内部边界两种，修补区域边界的引导向量值按照一般泊松融合边界条件运算，对于修补区域内部边界位置的引导向量值，计算仅采用来自相同参考图像的临近像素值，实现时采用4邻域或8邻域窗口滤波的方式。
与现有技术相比，本申请提供了一种针对遥感卫星多光谱图像数据自动云修补解决方案。利用多光谱数据产品包含的质量波段提取的高质量云与云下阴影像元，配合整饰算法自动生成高质量云掩膜。利用生成缩略图的方式实现快速计算SSIM相似性，并依据SSIM、图像获取时间、云盖量综合确定每幅图像对其它图像的相似性排序，进而快速确定修补每幅图像所需参考图像的数量与排序。最后在简单直接替换修补的基础上，利用泊松融合完成高质量云修补。算法完全自动化，整个云检测与云修补过程无需人机交互，用户仅需对最终的检测结果进行简单的检查。涉及的关键步骤采用成熟的算法实现，具有较高的稳定性与适用性。对于以Landsat系列卫星为代表的多光谱图像海量数据自动生产高质量云修补数据产品，提供了关键的技术支撑。本技术是基于时间序列的地表覆盖分类研究中必要的数据预处理过程，用于消除云与云下阴影对序列光谱的影响。相较于之前的基于简单滤波的方法明显具有更高的精度。
以下提供本申请的几个实施例，以阐述本申请的具体实施过程：
实施例1
实施例1中，实现遥感卫星多光谱图像自动云修补流程如图2所示：
步骤S201、多光谱图像质量波段(输入数据)：质量波段通常为多光谱图像产品数据自带数据，比如Landsat产品数据为*.tar.gz压缩包格式，解压后得到的质量波段数据为单波段16位整型tif图像，以不同的位数记录云、云下阴影、雪、水体等标记信息，其中原始DN值数据产品记录在*_bqa.tif文件中，地表反射率数据产品记录在*_pixel_qa.tif文件中。质量波段与对应产品的原始DN值数据、或者地表反射率数据具有相同的图像尺寸、地图投影及坐标，所以质量波段的像素位置与图像数据是对应的。数据版本方面，2018年之前的PreCollection版本云与云下阴影检测采用CAAC算法，2018年以后新版的Collection 1版本采用CFMask算法具有更好的云与云下阴影检测精度，建议使用Collection 1新版数据。对于中国的遥感卫星多光谱图像产品数据，质量波段的标记与Landsat产品类似，具体数据格式与标记规则存在很多差异。
实际处理数据时，有些多光谱数据的云盖量很大，一致达到100％，对于这类数据的云修补意义不大以致无法完成，所以在选择输入数据时，我们将辅助文件中标称云盖量大于80％的数据排除在处理之外。
步骤S202、云掩膜提取：多光谱图像产品数据多为*.tar.gz压缩包格式，假定有m个待处理的压缩包数据，m取值可以为千万数量级的海量数据，算法首先通过文件查询或者读取列表的方式获取所有m个数据的文件路径，并记录到内存中，然后将所有数据按照区域归类，即将相同的航带号的所有图像划归到同一区域。Landsat卫星航带编号采用WRS参考系，结合星下点成像特性而形成的固定地面参考网格，WRS网格的二维坐标采用Path和Row进行标识。Path与Row编号通过压缩包名称即可获取，例如文件名“LC081190292016071101T1-SC20180904092105.tar.gz”中的“LC08”后面的“119”为Path值，“029”为Row值。将m个数据中具有相同Path与Row的待处理数据被划归为同一区域。中国遥感卫星多光谱图像数据产品命名方式类似，通过文件名称同样可以按照区域划分数据。算法程序一次完成的处理流程完成对同一个区域数据的处理，直至所有区域处理完成。
假定某一区域一共包含n个待处理数据，算法首先自动将所有数据采用并行的方式快速解压到硬盘，并记录解压后的质量数据与地表反射率数据路径名称。云掩膜提取是将n个质量数据由16位图像简化为只标记云与云下阴影的8位图像。具体是将背景标记0值，清晰地表标记1值，云(Cloud、High ConfidenceCirrus)标记2值，云下阴影(Cloud Shadow)标记3值。云掩膜提取后得到n个云掩膜tif格式单波段8位Byte型图像，保存到硬盘。
步骤S203、云掩膜整饰：，直接利用质量波段提取的云掩膜图像不建议直接用于后续云修补，实验发现质量波段提取的掩膜图像主要存在两个问题：1)掩膜像素在云边界及薄云等处经常出现过于离散的情况，这主要是由于云检测算法目前是基于阈值的逐像素判断，并没有考虑相邻像素的位置关系，过于离散的掩膜会降低后续云修补的质量与稳定性，增加运算复杂度；2)掩膜像素在云边界处的精度不高，对于云下阴影的检测精度更低，在边界处经常出现云与云下阴影的漏检测。云掩膜整饰的目的是解决上述两个问题。
云掩膜整饰主要包含两步处理：
①减小离散化水平，具体方法是将云掩膜中小于4个像素的云与云下阴影重标记为清晰地表值1，将小于4个像素的清晰地表重标记为云值2；
②膨胀云与云下阴影区域，将标记2值的云区域以半径不超过5个像素的尺度膨胀，标记3值的云下阴影区域以半径不超过10个像素的尺度膨胀，扩充边界能在很大程度上解决漏检测的问题。图3是云掩膜整饰前后比较示例，其中，图3a为原始图像，图3b为QA波段提取云掩膜，图3c为云掩膜整饰结果。从图3a、图3b、图3c的云掩膜整饰前后对比可以看出，一些小像素碎块被去除，云下阴影区域相比较云区域尺寸膨胀的更多。
对于硬盘中的n个云掩膜tif图像，云掩膜整饰过程利用形态学算法优化，结合并行处理策略实现对云掩膜图像的快速自动更新。
步骤S204、云掩膜图像(中间结果)：云掩膜图像是与每个原始图像像素尺寸相同的单波段8位Byte型数据类型图像，像素值0代表背景，像素值1代表清晰地表，像素值2代表云，像素值3代表云下阴影，云掩膜图像保存到硬盘。
步骤S206、云掩膜重采样：为了提高后续相似性排序与确定云修补参考数据处理步骤的处理速度，对云掩膜进行重采样，按照最近邻插值方法生成缩略图，缩小比例为原始尺寸的四分之一。对于硬盘中的n个云掩膜tif图像，使用并行处理快速生成n个云掩膜缩略图保存到硬盘。
步骤S208、云掩膜缩略图(中间结果)：云掩膜缩略图是原始图像四分之一尺寸的单波段8位Byte型数据图像，像素值0代表背景，像素值1代表清晰地表，像素值2代表云，像素值3代表云下阴影，云掩膜缩略图保存到硬盘。
步骤205、多光谱地表反射率(输入数据)：多光谱图像产品数据通常为*.tar.gz压缩包格式，解压后得到的地表反射率数据为多个单波段16位有符号整型tif图像。每个tif图像记录由可见光到中波红外范围的不同光谱范围的波段。
如果数据为原始DN值，则通常为多个单波段16位无符号整型tif图像。
步骤S207、地表反射率重采样：为了提高后续相似性排序与确定云修补参考数据处理步骤的处理速度，对受大气影响明显的蓝光(0.450–0.515μm)单个波段的地表反射率图像进行重采样，按照最近邻插值方法生成缩略图，缩小比例为原始尺寸的四分之一。对于硬盘中的n个地表反射率tif图像，使用并行处理快速生成n个地表反射率缩略图保存到硬盘。
步骤S209、地表反射率缩略图(中间结果)：地表反射率缩略图是蓝光(0.450–0.515μm)波段的地表反射率图像四分之一尺寸的单波段16位整型数据图像，地表反射率缩略图保存到硬盘。
步骤S210、计算相似性：计算相似性是对于同区域的每一幅多光谱图像，依据SSIM、图像获取时间、云盖量计算其它图像与该图像的相似性。
SSIM的计算使用缩略图减小运算量，提高效率。对于硬盘中的n个多光谱地表反射率缩略图，以及对应的云掩膜缩略图，任意两幅缩略图Ii与Ij，其中i,j∈n，统计Ii与Ij中公共的清晰地表区域，即对应云掩膜区域值都为1的像素位置，的均值μi与μj，标准差σi与σj，协方差σij。则Ii与Ij的SSIM计算公式如下：
其中C为常数，比如C＝2，目的是使得SSIM的取值范围介于-1到1之间。SSIM值越大代表两幅图像相似性越高。
图像获取时间记录在多光谱地表反射率产品数据的辅助文件中*_MTL.txt，记两幅缩略图Ii与Ij的获取时间间隔天数为Tij，则Tij值越小代表两幅图像相似性越高。
云盖量通过统计缩略图Ii与Ij对应的云掩膜缩略图Mi与Mj中的云与云下阴影区域像素数得到，记Ii与Ij的云盖量分别Ci与Cj，记公共云与云下阴影区域为Cij，则对于图像Ii，当图像Ij的Cj与Cij的取值越小越有利于后续云区的修补。
最后，Ii对Ij的相似性Sij计算公式如下：
其中Mij为云掩膜缩略图Mi与Mj中公共的非0像素数量，包含了非背景区域的清晰地表与云区，a、b、c为经验常数，例如a＝1、b＝1、c＝1，这样Sij的值越大代表Ii对Ij的相似性越高。
步骤S211、确定云修补参考数据：对同区域的每一个多光谱图像确定对其进行云修补所需的参考图像数量、文件名及修补顺序。
确定云修补参考数据的计算使用缩略图减小运算量，提高效率。对于硬盘中的n个多光谱地表反射率图像，记当前处理到第i幅图像，则对云掩膜缩略图Mi，需在其余图像中查找与i图相似性值的最大的j图作为第一个参考图像，即相似性
虚拟云修补则是使用该参考图像的云掩膜缩略图Mj，将在Mj中的清晰地表填充到对应地理位置下Mi的云与云下阴影区域中，即对于Mi所谓像素点，如果值大于1，而该位置在Mj中值等于1，则将Mi的该位置值设为1。依次类推，直到云掩膜Mi中的云区修补完毕，或者用尽所有参考图像。记录修补第i幅图所需的参考图像数量、文件名及修补顺序依次类推，直到所有n幅图像全部确定了云修补所用的参考图像。
需要注意的是，实际处理数据过程中，由于云检测结果的复杂性，真实数据往往很少有能将云掩膜中的所有云区域都填充完毕的情况，更多的是处理到填充区域不再减少为止，此时剩余区域通常是一些高亮的地表，在所有参考图像的云掩膜中都被误检测为云。一般情况下修补一幅图像所需参考图像的数量在5幅左右。
步骤S212、云修补参考查找表(中间结果)：记录同区域的每一个多光谱图像云修补所需要的参考图像数量、文件名、修补顺序，保存为外部文本文件。
步骤S213、直接替换法云修补，对同区域的每一个多光谱图像，对于其云掩膜所标记的云与云下阴影待修补区域，根据云修补参考查找表所确定的参考图像修补顺序，将参考图像非云区域像素值逐波段替换到待修补区域，得到替换修补图像。同时为了标示图像需补区域每个像素的修补来源，记录替换参考掩膜。
记当前修补图像为T，对应云掩膜MT，修补参考图像为Ri，参考图像云掩膜为
i∈[1,m]，m为参考图像的数量，记Tj(x,y)表示修补图像第j个波段像素位置为(x,y)的像素,j∈[1,b]，b为波段数量。记替换参考掩膜为
则直接替换法云修补的处理过程为:
图5a是原始图像，图5b是替换参考掩膜示例，修补区域内不同灰度级表示来自不同参考图像。图5c是直接替换法云修补结果示例，可见直接替换结果在修补区域边界及修补区域内部不同参考图像边界存在明显目视辐射差异，所以需要后续基于泊松融合的云修补消除辐射差异。
步骤S214、替换参考掩膜(中间结果)：用于标示图像需补区域每个像素的修补来源，为与云掩膜尺寸相同的单波段8位图像，其中像素0值对应云掩膜图像标记的背景值与清晰地表值，非0值对应云掩膜图像中的云与云下阴影区域，非0值取值于修补云区所用参考图像的编号值，比如i值表示该位置像素云修补时来自第i幅参考图像，i∈[1,m]，m为参考图像数量。
每幅图像的替换参考掩膜与使用的参考图像文件名称列表作为最终云修补产品数据的辅助文件提供给用户，方便用户追溯每幅图像的修补区域与每个像素的修补来源，这对于后续可能的地表覆盖变化监测研究是必要的。
步骤S215、替换修补数据(中间结果)：与原始多光谱地表反射率数据格式相同，差异在于其云掩膜所标记的云与云下阴影待修补区域像素值被替换为来自不同参考图像的像素值。
步骤S216、泊松融合法云修补：对同区域的每一个多光谱替换修补图像与替换参考掩膜图像，利用泊松融合完成云修补，得到最终多光谱云修补图像。
泊松融合目前在无人机图像无缝镶嵌领域有着广泛的应用。无人机图像无缝拼接所解决的问题是所有图形之间匀色与消除拼接缝，运算完毕后结果图像的所有像素值都会发生变化，整幅图像达到均匀增强的色彩表现，而与之不同的是，本算法使用的泊松融合要求仅对修补区域的像素值进行处理，修补区域外的像素值不变，并且需要消除修补区域内不同参考图像像素区域之间的边界。
图4是本算法使用的泊松融合法云修补的算法原理示意图，下面结合图4中的定义描述泊松融合法云修补的算法实施原理。记S为云修补的目标图像，Ω为待修补的云与云下阴影区域，改区域的边界像素记
记f*为S减去Ω范围的标量函数，记f为Ω范围的未知标量函数，g为直接替换法得到的修补区域像素值标量函数，其中g1与g2表示来自不同的参考图像，所以g表示为分段函数的形式，g1与g2之间具有复杂的边界，由g计算得到的引导向量域记为V。
泊松融合的目的是求解f，使得修补区域与S具有最佳一致性的过渡，运算上使得Ω内部梯度变化最小，记
为梯度运算，则有
且在边界处有
记
为拉普拉斯算子，当梯度的二阶偏导为0时取得极小值，由拉格朗日方程，有Δf|Ω＝0，表示为引导向量域最小化问题，有
其解是狄利克雷边界条件下泊松方程的唯一解，即Δf|Ω＝divV，同样在边界处有
其中
是V＝(u,v)的散度。
泊松融合的目的是求解Ω内的修正函数
使得
因此在Ω内需要的修正
便通过边界上的
与误差f*-g作为边界条件，向Ω内插值求解。即
且在边界处有
本算法采用的泊松融合求解过程仍然可以使用无人机图像无缝拼接的编程实现，这在目前很多开源项目中都有涉及，并且包括多种加速优化版本可供选择。不同的是引导向量域V的求解，本算法计算V的方式略复杂些，将无人机图像泊松融合的统一边界类型细化，分成了修补区域边界与修补区域内部边界两种。
修补区域边界，即
修补区域外V的值为0，修补区域边界的V值按照一般泊松融合边界条件运算，修补区域内V的值按照一般泊松融合的引导向量域运算，具体编程实现时采用4邻域或8邻域窗口滤波的方式。例如4邻域修补区域内V的值计算公式为：
V(x,y)＝g(x+1,y)+g(x,y+1)+g(x-1,y)+g(x,y-1)-4×g(x,y)
而使用8邻域修补区域内V的值计算公式为：
V(x,y)＝g(x+1,y)+g(x,y+1)+g(x-1,y)+g(x,y-1)+g(x+1,y+1)+g(x-1,y+1)+g(x+1,y-1)+g(x-1,y-1)-8×g(x,y)
修补区域内部边界，即表示来自不同的参考图像g1与g2之间的边界，见图4。对于修补区域内部边界位置的像素V值，计算仅采用来自相同参考图像的临近像素值。例如使用4邻域窗口滤波时，图4中g1区域的像素点c的V值与g2区域像素点4的V值计算公式分别为：
V(c)＝g1(b)+g1(e)-2×g1(c)
V(4)＝g2(2)+g2(5)-2×g2(4)
而使用8邻域窗口滤波时，图4中g1区域的像素点c的V值与g2区域像素点4的V值计算公式分别为：
V(c)＝g1(a)+g1(b)+g1(d)+g1(e)+g1(f)-5×g1(c)
V(4)＝g2(1)+g2(2)+g2(3)+g2(5)+g2(6)-5×g2(4)
通过对引导向量域的定制，本算法可以很好地保持修补区域的光谱特性，并且不会改变未修补区域的像素值，这对于数据后续的遥感分类研究非常重要。
图5d是泊松融合云修补结果示例，可见基于泊松融合的云修补有效消除了直接替换结果在修补区域边界及修补区域内部不同参考图像边界存在的明显目视辐射差异。
步骤S217、多光谱图像云修补产品(输出数据)：最终云修补产品数据，与原始多光谱地表反射率数据格式相同，差异在于其云掩膜所标记的云与云下阴影待修补区域像素值被替换为来自不同参考图像的像素值，并通过泊松融合完成光谱修正。每幅云修补产品数据按照区域放置在不同的以Path、Row及原始文件名命名的文件夹下，例如：../119/029/LC081190292013041401T1-SC20180831224406/。目前每幅云修补产品包含7个文件，包括：
*sr_cr.tif为云修补图像，包含所有波段，采用LZW压缩存储，*为原始数据文件名
*sr_crm.tif为替换参考掩膜，单波段byte型LZW压缩存储，像素值对应参考图像编号
*sr_cr.txt文本文件记录云盖量统计与当前图像所采用的修补图像编号与名称
*.jpg为修补图像与掩膜图像的jpg缩略图，
*MTL/ANG.txt为C1产品原始辅助文件拷贝
本申请的C++算法实例已经在PC平台上实现，前期通过大量实验数据验证了算法的有效性与鲁棒性，目前已经用于海量数据的云修补产品生产，并完成了第一个版本的中国区域2013年至2017年连续5年Landsat Collection 1版本地表反射率数据的云修补，完整覆盖中国区域包含528个航带编号，数据总量超过5万幅，约30TB。图6是中国东北部(path:119/row:029)区域部分图像的云修补前后及替换修补参考掩膜示例图，其中列61对应的是原始图像，列62对应的是替换修补参考掩膜，列63对应的是修补结果图像。
图7示出根据本申请一实施例的卫星图像序列云区修补装置结构框图。
如图7所示，该卫星图像序列云区修补装置包括：
云掩膜提取整饰模块71，用于从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像；
相似度计算模块72，用于根据各所述单波段图像之间的相似度，计算修补所述多光谱图像所需要的至少一个参考图像，以及各所述参考图像的修补顺序，所述参考图像选自于多张所述多光谱图像中一张或多张；
参考数据选定模块73，用于根据所述单波段图像对应的云掩膜图像，确定所述多光谱图像的修补区域；
直接替换修补模块74，用于根据所述修补顺序，将所述多光谱图像的修补区域的像素值，替换为所述参考图像的非云区域的像素值，得到替换修补图像；
替换记录模块75，用于记录所述修补区域中每张多光谱图像对应的各参考图像的编号值，得到替换参考掩膜图像；
泊松融合修补模块76，用于通过泊松融合对所述替换修补图像和所述替换参考掩膜图像进行云修补。
图8示出根据本申请一实施例的云掩膜提取整饰模块结构框图。
如图8所示，该云掩膜提取整饰模块71包括：
初始单元711，用于从所述多光谱图像中识别云区域和云下阴影区域，生成待选云掩膜图像；
标记单元712，用于将所述待选云掩膜图像中小于第一预设像素值的云区域和云下阴影区域重新标记为清晰地表，将小于第二预设像素值的清晰地表重新标记为云区域；
膨胀单元713，用于将重新标记的云区域以半径不超过第一扩展像素的尺度膨胀，并将重新标记的云下阴影区域以半径不超过第二扩展像素的尺度膨胀，获得所述云掩膜图像。
图9示出根据本申请一实施例的相似度计算模块结构框图。
如图9所示，该相似度计算模块72包括：
插值单元721，用于对多张所述多光谱图像和多张所述云掩膜图像，按照预设波段采用最近邻插值法，生成各所述单波段图像的缩略图和各所述云掩膜图像的缩略图；
计算单元722，用于根据各所述单波段图像的缩略图，计算各所述单波段图像之间的相似度；
填充单元723，用根据各所述参考图像与所述单波段图像的相似度排序，依次使用所述参考图像的缩略图中的非云区域像素，填充所述单波段图像对应的多光谱图像中的标记区域，所述标记区域由所述单波段图像对应的云掩膜图像的缩略图标记；
记录单元724，用于在所述标记区域被填充完毕或填充区域不再减少的情况下，记录填充过程中用到的参考图像和对应的修补顺序。
图10示出根据本申请一实施例的泊松融合修补模块结构框图。
如图10所示，该泊松融合修补模块76包括：
设置单元761，用于将所述多光谱图像修补区域外的引导向量值为0。
划分单元762，用于将所述多光谱图像的边界划分为修补区域边界和修补区域内部边界；
融合单元763，用于对于所述修补区域内部边界位置的引导向量值，采用来自相同参考图像的邻近像素值。
此处描述的系统和技术的各种实施方式可以在数字电子电路系统、集成电路系统、专用集成电路(Application Specific Integrated Circuits，ASIC)、计算机硬件、固件、软件、和/或它们的组合中实现。这些各种实施方式可以包括：实施在一个或者多个计算机程序中，该一个或者多个计算机程序可在包括至少一个可编程处理器的可编程系统上执行和/或解释，该可编程处理器可以是专用或者通用可编程处理器，可以从存储系统、至少一个输入装置、和至少一个输出装置接收数据和指令，并且将数据和指令传输至该存储系统、该至少一个输入装置、和该至少一个输出装置。
这些计算程序(也称作程序、软件、软件应用、或者代码)包括可编程处理器的机器指令，并且可以利用高级过程和/或面向对象的编程语言、和/或汇编/机器语言来实施这些计算程序。如本文使用的，术语“机器可读介质”和“计算机可读介质”指的是用于将机器指令和/或数据提供给可编程处理器的任何计算机程序产品、设备、和/或装置(例如，磁盘、光盘、存储器、可编程逻辑装置(programmable logic device，PLD))，包括，接收作为机器可读信号的机器指令的机器可读介质。术语“机器可读信号”指的是用于将机器指令和/或数据提供给可编程处理器的任何信号。
为了提供与用户的交互，可以在计算机上实施此处描述的系统和技术，该计算机具有：用于向用户显示信息的显示装置(例如，CRT(Cathode RayTube，阴极射线管)或者LCD(液晶显示器)监视器)；以及键盘和指向装置(例如，鼠标或者轨迹球)，用户可以通过该键盘和该指向装置来将输入提供给计算机。其它种类的装置还可以用于提供与用户的交互；例如，提供给用户的反馈可以是任何形式的传感反馈(例如，视觉反馈、听觉反馈、或者触觉反馈)；并且可以用任何形式(包括声输入、语音输入或者、触觉输入)来接收来自用户的输入。
可以将此处描述的系统和技术实施在包括后台部件的计算系统(例如，作为数据服务器)、或者包括中间件部件的计算系统(例如，应用服务器)、或者包括前端部件的计算系统(例如，具有图形用户界面或者网络浏览器的用户计算机，用户可以通过该图形用户界面或者该网络浏览器来与此处描述的系统和技术的实施方式交互)、或者包括这种后台部件、中间件部件、或者前端部件的任何组合的计算系统中。可以通过任何形式或者介质的数字数据通信(例如，通信网络)来将系统的部件相互连接。通信网络的示例包括：局域网(Local Area Network，LAN)、广域网(Wide Area Network，WAN)和互联网。
计算机系统可以包括客户端和服务器。客户端和服务器一般远离彼此并且通常通过通信网络进行交互。通过在相应的计算机上运行并且彼此具有客户端-服务器关系的计算机程序来产生客户端和服务器的关系。
应该理解，可以使用上面所示的各种形式的流程，重新排序、增加或删除步骤。例如，本发申请中记载的各步骤可以并行地执行也可以顺序地执行也可以不同的次序执行，只要能够实现本申请公开的技术方案所期望的结果，本文在此不进行限制。
上述具体实施方式，并不构成对本申请保护范围的限制。本领域技术人员应该明白的是，根据设计要求和其他因素，可以进行各种修改、组合、子组合和替代。任何在本申请的精神和原则之内所作的修改、等同替换和改进等，均应包含在本申请保护范围之内。
Claims (10)
1.一种卫星图像序列云区修补方法，其特征在于，包括：
从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像；
对多张所述多光谱图像和所述多张单波段图像对应的多张云掩膜图像，按照预设波段采用最近邻插值法，生成各所述单波段图像的缩略图和各所述云掩膜图像的缩略图；
根据各所述单波段图像的缩略图，计算结构相似性指数和云盖量，通过所述结构相似性指数、所述云盖量、缩略图的获取时间间隔，计算各所述单波段图像之间的相似度；
根据各所述单波段图像之间的相似度，计算修补所述多光谱图像所需要的至少一个参考图像，以及各所述参考图像的修补顺序，所述参考图像选自于多张所述多光谱图像中一张或多张；
根据所述单波段图像对应的云掩膜图像，确定所述多光谱图像的修补区域；
根据所述修补顺序，将所述多光谱图像的修补区域的像素值，替换为所述参考图像的非云区域的像素值，得到替换修补图像；
记录所述修补区域中每张多光谱图像对应的各参考图像的编号值，得到替换参考掩膜图像；
通过泊松融合对所述目标区域中的各所述替换修补图像和各所述替换参考掩膜图像进行云修补。
2.根据权利要求1所述的方法，其特征在于，从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像，包括：
从所述多光谱图像中识别云区域和云下阴影区域，生成待选云掩膜图像；
将所述待选云掩膜图像中小于第一预设像素值的云区域和云下阴影区域重新标记为清晰地表，将小于第二预设像素值的清晰地表重新标记为云区域；
将重新标记的云区域以半径不超过第一扩展像素的尺度膨胀，并将重新标记的云下阴影区域以半径不超过第二扩展像素的尺度膨胀，获得所述云掩膜图像。
3.根据权利要求1所述的方法，其特征在于，根据各所述单波段图像之间的相似度，计算修补所述多光谱图像所需要的至少一个参考图像，以及各所述参考图像的修补顺序，包括：
根据各所述参考图像与所述单波段图像的相似度排序，依次使用所述参考图像的缩略图中的非云区域像素，填充所述单波段图像对应的多光谱图像中的标记区域，所述标记区域由所述单波段图像对应的云掩膜图像的缩略图标记；
在所述标记区域被填充完毕或填充区域不再减少的情况下，记录填充过程中用到的参考图像和对应的修补顺序。
4.根据权利要求1所述的方法，其特征在于，通过泊松融合对所述目标区域中的各所述替换修补图像和各所述替换参考掩膜图像进行云修补，包括：
所述多光谱图像修补区域外的引导向量值为0；
将所述多光谱图像的边界划分为修补区域边界和修补区域内部边界；
对于所述修补区域内部边界位置的引导向量值，采用来自相同参考图像的邻近像素值。
5.根据权利要求3所述的方法，其特征在于，通过所述结构相似性指数、所述云盖量、缩略图的获取时间间隔，计算各所述单波段图像之间的相似度，包括：
通过结构相似性指数、图像获取时间、云盖量计算所述相似度，计算公式如下：
其中，a、b、c为加权系数，Sij为所述单波段图像的缩略图Ii与Ij的相似度，SSIM(Ii,Ij)为所述缩略图Ii与Ij的结构相似性指数，Tij为所述缩略图Ii与Ij的获取时间间隔，Cj为所述缩略图Ij的云盖量，Cij为所述缩略图Ii与Ij的公共云区域与云下阴影区域，Mij为云掩膜缩略图Mi与Mj中公共的非0像素数量。
6.根据权利要求5所述的方法，其特征在于，所述结构相似性指数的计算公式为:
其中，SSIM(Ii,Ij)为所述缩略图Ii与Ij的结构相似性指数，μi与μj为所述缩略图Ii与Ij的对应云掩膜区域值都为1的像素位置的均值，σi与σj为所述缩略图Ii与Ij的对应云掩膜区域值都为1的像素位置的标准差，σij为所述缩略图Ii与Ij的对应云掩膜区域值都为1的像素位置的协方差，所述云掩膜区域值包括云区域值和云下阴影区域值，C为常数。
7.一种卫星图像序列云区修补装置，其特征在于，包括：
云掩膜提取整饰模块，用于从目标区域的多张多光谱图像中分别识别并扩大云区域和云下阴影区域，生成多张单波段图像；
相似度计算模块，用于对多张所述多光谱图像和所述多张单波段图像对应的多张云掩膜图像，按照预设波段采用最近邻插值法，生成各所述单波段图像的缩略图和各所述云掩膜图像的缩略图；根据各所述单波段图像的缩略图，计算结构相似性指数和云盖量，通过所述结构相似性指数、所述云盖量、缩略图的获取时间间隔，计算各所述单波段图像之间的相似度；根据各所述单波段图像之间的相似度，计算修补所述多光谱图像所需要的至少一个参考图像，以及各所述参考图像的修补顺序，所述参考图像选自于多张所述多光谱图像中一张或多张；
参考数据选定模块，用于根据所述单波段图像对应的云掩膜图像，确定所述多光谱图像的修补区域；
直接替换修补模块，用于根据所述修补顺序，将所述多光谱图像的修补区域的像素值，替换为所述参考图像的非云区域的像素值，得到替换修补图像；
替换记录模块，用于记录所述修补区域中每张多光谱图像对应的各参考图像的编号值，得到替换参考掩膜图像；
泊松融合修补模块，用于通过泊松融合对所述替换修补图像和所述替换参考掩膜图像进行云修补。
8.根据权利要求7所述的装置，其特征在于，所述云掩膜提取整饰模块，包括：
初始单元，用于从所述多光谱图像中识别云区域和云下阴影区域，生成待选云掩膜图像；
标记单元，用于将所述待选云掩膜图像中小于第一预设像素值的云区域和云下阴影区域重新标记为清晰地表，将小于第二预设像素值的清晰地表重新标记为云区域；
膨胀单元，用于将重新标记的云区域以半径不超过第一扩展像素的尺度膨胀，并将重新标记的云下阴影区域以半径不超过第二扩展像素的尺度膨胀，获得所述云掩膜图像。
9.根据权利要求7所述的装置，其特征在于，所述相似度计算模块，包括：
填充单元，用于根据各所述参考图像与所述单波段图像的相似度排序，依次使用所述参考图像的缩略图中的非云区域像素，填充所述单波段图像对应的多光谱图像中的标记区域，所述标记区域由所述单波段图像对应的云掩膜图像的缩略图标记；
记录单元，用于在所述标记区域被填充完毕或填充区域不再减少的情况下，记录填充过程中用到的参考图像和对应的修补顺序。
10.根据权利要求7所述的装置，其特征在于，所述泊松融合修补模块，包括：
设置单元，用于将所述多光谱图像修补区域外的引导向量值为0；
划分单元，用于将所述多光谱图像的边界划分为修补区域边界和修补区域内部边界；
融合单元，用于对于所述修补区域内部边界位置的引导向量值，采用来自相同参考图像的邻近像素值。
Priority Applications (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910857550.8A CN110555818B (zh) | 2019-09-09 | 2019-09-09 | 一种卫星图像序列云区修补方法和装置 |
| US17/641,571 US20230094679A1 (en) | 2019-09-09 | 2019-10-17 | Satellite Image Sequence Cloud Region Repairing Method and Apparatus |
| PCT/CN2019/111581 WO2021046965A1 (zh) | 2019-09-09 | 2019-10-17 | 一种卫星图像序列云区修补方法和装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910857550.8A CN110555818B (zh) | 2019-09-09 | 2019-09-09 | 一种卫星图像序列云区修补方法和装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110555818A CN110555818A (zh) | 2019-12-10 |
| CN110555818B true CN110555818B (zh) | 2022-02-18 |
Family
ID=68739953
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910857550.8A Active CN110555818B (zh) | 2019-09-09 | 2019-09-09 | 一种卫星图像序列云区修补方法和装置 |
Country Status (3)
| Country | Link |
|---|---|
| US (1) | US20230094679A1 (zh) |
| CN (1) | CN110555818B (zh) |
| WO (1) | WO2021046965A1 (zh) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20230245445A1 (en) * | 2020-07-13 | 2023-08-03 | Esen Sistem Entegrasyon Ve Muhendislik Hizm. San Ve Tic Ltd. Sti. | An object detection method |
| CN112102180B (zh) * | 2020-08-21 | 2022-10-11 | 电子科技大学 | 一种基于Landsat光学遥感图像的云识别方法 |
| CN112465846B (zh) * | 2020-11-26 | 2023-04-07 | 西安电子科技大学 | 基于填充策略的含云遥感图像压缩方法 |
| CN113343871B (zh) * | 2021-06-17 | 2022-03-18 | 哈尔滨工业大学 | 基于泊松融合+直方图匹配的高分四号多光谱遥感影像中舰船目标的仿真方法及系统 |
| CN113656419B (zh) * | 2021-07-30 | 2023-06-13 | 北京市遥感信息研究所 | 全球地表反射率数据集构建及更新方法及装置 |
| CN113790818B (zh) * | 2021-08-10 | 2023-09-26 | 中国电子科技集团公司第十三研究所 | 可见光热反射测温方法及测温设备 |
| CN113838065B (zh) * | 2021-09-23 | 2022-06-21 | 江苏天汇空间信息研究院有限公司 | 一种基于影像标记的自动化去云方法 |
| CN114821136B (zh) * | 2022-04-11 | 2023-04-21 | 成都信息工程大学 | 一种自适应的云微粒子图像数据处理方法 |
| CN116503742B (zh) * | 2023-06-26 | 2023-09-08 | 自然资源部第二海洋研究所 | 一种基于超像素和图割算法的遥感图像云替补方法 |
| CN117575970B (zh) * | 2024-01-15 | 2024-04-16 | 航天宏图信息技术股份有限公司 | 基于分类的卫星影像自动处理方法、装置、设备及介质 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104008543A (zh) * | 2014-05-12 | 2014-08-27 | 河海大学 | 一种图像融合质量评价方法 |
| CN105046708A (zh) * | 2015-07-14 | 2015-11-11 | 福州大学 | 一种与主观感知相一致的颜色校正客观评估方法 |
| CN106485740A (zh) * | 2016-10-12 | 2017-03-08 | 武汉大学 | 一种结合稳定点和特征点的多时相sar图像配准方法 |
| CN106940887A (zh) * | 2017-03-09 | 2017-07-11 | 中国科学院遥感与数字地球研究所 | 一种gf‑4卫星序列图像云与云下阴影检测方法 |
| CN107103295A (zh) * | 2017-04-20 | 2017-08-29 | 苏州中科天启遥感科技有限公司 | 光学遥感影像云检测方法 |
| CN108230376A (zh) * | 2016-12-30 | 2018-06-29 | 北京市商汤科技开发有限公司 | 遥感图像处理方法、装置和电子设备 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20070119879A (ko) * | 2006-06-16 | 2007-12-21 | 삼성전자주식회사 | 영상의 해상도 변환 방법 및 그 장치 |
| CN105678777B (zh) * | 2016-01-12 | 2018-03-13 | 武汉大学 | 一种多特征联合的光学卫星影像云与云阴影检测方法 |
| CN107230195B (zh) * | 2017-07-12 | 2020-09-18 | 中国科学院遥感与数字地球研究所 | 一种影像处理方法和装置 |
| CN108109126B (zh) * | 2018-01-12 | 2020-10-27 | 适普远景遥感信息技术（北京）有限公司 | 一种基于卫星遥感影像的目标区域填充与融合处理方法 |
-
2019
- 2019-09-09 CN CN201910857550.8A patent/CN110555818B/zh active Active
- 2019-10-17 WO PCT/CN2019/111581 patent/WO2021046965A1/zh active Application Filing
- 2019-10-17 US US17/641,571 patent/US20230094679A1/en active Pending
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104008543A (zh) * | 2014-05-12 | 2014-08-27 | 河海大学 | 一种图像融合质量评价方法 |
| CN105046708A (zh) * | 2015-07-14 | 2015-11-11 | 福州大学 | 一种与主观感知相一致的颜色校正客观评估方法 |
| CN106485740A (zh) * | 2016-10-12 | 2017-03-08 | 武汉大学 | 一种结合稳定点和特征点的多时相sar图像配准方法 |
| CN108230376A (zh) * | 2016-12-30 | 2018-06-29 | 北京市商汤科技开发有限公司 | 遥感图像处理方法、装置和电子设备 |
| CN106940887A (zh) * | 2017-03-09 | 2017-07-11 | 中国科学院遥感与数字地球研究所 | 一种gf‑4卫星序列图像云与云下阴影检测方法 |
| CN107103295A (zh) * | 2017-04-20 | 2017-08-29 | 苏州中科天启遥感科技有限公司 | 光学遥感影像云检测方法 |
Non-Patent Citations (3)
| Title |
|---|
| Automatic Cloud Removal from Multi-Temporal Landsat Collection 1 Data Using Poisson Blending;Changmiao Hu et al.;《2019 IEEE International Geoscience and Remote Sensing Symposium》;20190802;正文第2章 * |
| Cloud Removal From Multitemporal Satellite Images Using Information Cloning;Chao-Hung Lin et al.;《IEEE Transactions on Geoscience and Remote Sensing》;20120608;第51卷(第1期);全文 * |
| 多时相遥感影像厚云去除;黄微 等;《遥感信息》;20160315;第31卷(第6期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| WO2021046965A1 (zh) | 2021-03-18 |
| CN110555818A (zh) | 2019-12-10 |
| US20230094679A1 (en) | 2023-03-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110555818B (zh) | 一种卫星图像序列云区修补方法和装置 | |
| Zhu et al. | A modified neighborhood similar pixel interpolator approach for removing thick clouds in Landsat images | |
| Makarau et al. | Haze detection and removal in remotely sensed multispectral imagery | |
| Mostafa | A review on various shadow detection and compensation techniques in remote sensing images | |
| Soille | Morphological image compositing | |
| CN109872278B (zh) | 基于u形网络和生成对抗网络的图像云层移除方法 | |
| CN107067405B (zh) | 基于尺度优选的遥感影像分割方法 | |
| US8503761B2 (en) | Geospatial modeling system for classifying building and vegetation in a DSM and related methods | |
| CN109308688B (zh) | 一种可见光和近红外波段厚云及阴影去除方法 | |
| Chen et al. | Structure-aware image inpainting using patch scale optimization | |
| Grigillo et al. | Automated building extraction from IKONOS images in suburban areas | |
| CN115951350A (zh) | 永久散射体点提取方法、装置、设备及介质 | |
| Ju et al. | Remote sensing image haze removal using gamma-correction-based dehazing model | |
| CN111310771A (zh) | 遥感影像的道路图像提取方法、装置、设备及存储介质 | |
| Su et al. | Superpixel-based weighted collaborative sparse regression and reweighted low-rank representation for hyperspectral image unmixing | |
| CN114092326A (zh) | 一种Web地图瓦片更新处理方法及装置 | |
| Albanwan et al. | 3D iterative spatiotemporal filtering for classification of multitemporal satellite data sets | |
| CN112967305A (zh) | 一种复杂天空场景下的图像云背景检测方法 | |
| Yu et al. | Colour balancing of satellite imagery based on a colour reference library | |
| CN109978982B (zh) | 一种基于倾斜影像的点云快速上色方法 | |
| Amirkolaee et al. | Convolutional neural network architecture for digital surface model estimation from single remote sensing image | |
| CN115526795A (zh) | 一种基于区域匹配和颜色迁移的无人机图像阴影补偿方法 | |
| CN113516059B (zh) | 固体废弃物的识别方法、装置、电子设备及存储介质 | |
| CN115456886A (zh) | 一种基于深度学习和光照模型的航空遥感影像阴影去除方法 | |
| CN115082533A (zh) | 一种基于自监督的临近空间遥感图像配准方法 |
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 |