CN103020939B - 利用多时相数据去除光学遥感影像大面积厚云的方法 - Google Patents

Abstract

本发明公开了一种利用多时相数据去除光学遥感影像大面积厚云的方法：如果光学遥感影像中存在大面积的厚云，而该区域其他时相的影像存在无云数据，则可利用它们的互补信息对云区数据进行修复重建。首先利用所有时相的无云数据进行字典学习，自适应地顾及影像间的相关性，学习出一个过完备的字典及该影像的最佳稀疏表示系数，最后对厚云区的数据修复重建。本发明利用不同的多时相影像厚云区的互补信息，以各影像的相关性作为权重，借助新兴的稀疏表达理论填补影像厚云区数据，不仅取得了较高的精度，还拓展了大面积厚云去除的思路，具有重要的实际意义。

Description

利用多时相数据去除光学遥感影像大面积厚云的方法

技术领域

本发明提出一种利用多时相数据去除光学遥感影像大面积厚云的方法，以自适应光谱加权的字典学习方法有效利用其他时相数据与待处理影像厚云区的相关性，实现厚云区数据的填补，涉及光学遥感图像处理技术领域。

背景技术

在实际成像时，光学遥感影像通常受到天气条件的影响，其中最为常见的影响因素就是云。大多数光学传感器获得的遥感影像都会受到不同程度的云污染，云层掩盖了下层地物的真实信息。尤其是厚云，完全掩盖了地物辐射信息，造成影像信息的缺失，明显降低了影像质量，严重影响了后续的影像分析工作和进一步应用。因此，有效去除厚云，对于改善影像质量，有效提取地面目标物体信息，避免数据浪费，提高影像的可用性和利用率，具有重要的意义。

目前，国际上遥感影像厚云去除方法主要有两大类：

针对单幅影像的厚云，由于没有其他的参考信息，往往是对影像自身从补绘角度出发进行去除。发展比较成熟的主要有内插的方法，例如根据厚云区邻域无云区数据进行最邻近、双线性、双三次内插以及克里金插值（Kriging）等；偏微分方程（PDE）的方法，通过确定厚云区邻域扩散的方向和信息，实现厚云区边界外围信息自动向内扩散，达到云区填补的目的；总变分（TV）方法，由待修复厚云区边界点逐步向内进行各向异性的扩散，但边界处往往比较模糊；纹理合成的方法，能够较好地保持厚云区的结构和纹理信息。这类方法利用了厚云区与周围无云区的统计学特性、空间实体分布特性和自相关性，对较小的厚云区有比较好的去除效果，但对于大面积的厚云却乏而无力。

而针对多时相影像数据，人们则考虑利用该系列影像的互补信息进行厚云去除处理。一方面，对于影像中待处理厚云区域，可采用同区域、相近时相的影像数据进行替代，即完成厚云区的数据修复。该方法要求其他时相的影像在时间上比较靠近，地物无明显变化，且多幅影像的厚云区域不能相互重叠，这使得其实际应用价值大打折扣。另外即使找到了符合要求的多时相互补影像，在替换厚云区时，由于相近时相的影像之间存在色调差异，镶嵌的影像会出现明显的痕迹，很难达到无缝拼接的效果；虽然可以通过一定的技术调节拼接的可视化效果，但会影响到进一步的定量分析应用。另一方面，也有利用多时相影像融合达到插补去厚云的目的。这类方法同样要求厚云区域不能重叠，且厚云的去除不彻底，容易影响影像的判读。

上述方法对于去除小块厚云具有较好的效果，但对于通常存在的大面积厚云，往往不能有效利用多时相影像间的互补信息，故而不能彻底去除。

发明内容

针对现有去云方法的缺点，本发明提出一种利用多时相数据去除光学遥感影像大面积厚云的方法，充分挖掘多时相影像的互补信息，采用光谱自适应加权字典学习方法填充影像的大面积厚云区域。

本发明的技术方案为一种利用多时相数据去除光学遥感影像大面积厚云的方法，包括以下步骤：

步骤1，将所需处理的多时相影像序列进行几何校正，得到同一区域不同时相的影像；

步骤2，对同一区域不同时相的影像进行厚云区检测，得到各时相影像的云掩膜，并计算各时相影像非厚云区的相关系数；

步骤3，将同一区域不同时相的影像按照时间序列组成多维的空间光谱影像，对空间光谱影像划分影像子块，将所有的影像子块重组为二维的矩阵；将步骤2所得各时相影像的云掩膜按照时间序列组成空间光谱云掩膜，对空间光谱云掩膜划分掩膜子块，将所有的掩膜子块重组为二维的矩阵；

步骤4，根据以上的影像子块和掩膜子块的矩阵进行字典学习，重建影像厚云区域。

而且，步骤2进行厚云区检测时，将影像中的异常值也标记为厚云区。

而且，步骤3中，对空间光谱影像采用尺寸为n×n的窗口滑动划分成大小为n×n×k的影像子块，其中，k是影像的光谱维数，即同一区域不同时相的影像的个数，然后将每一个影像子块重组为一个列向量，并将所有的列向量按滑动顺序组成一个二维的矩阵；并且，对空间光谱云掩膜采用尺寸为n×n的窗口滑动划分成大小为n×n×k的掩膜子块，然后将每一个掩膜子块重组为一个列向量，并将所有的列向量按滑动顺序组成一个二维的矩阵。

而且，步骤4中设字典的大小为kn²×m，m≥256且m＞kn²；根据步骤2所得各时相影像非厚云区的相关系数，自适应地确定多时相数据在字典学习过程中的权重。

本发明提出的利用多时相数据去除光学遥感影像大面积厚云方法，它是针对同一地区、不同时相所获得的一组时间序列影像进行字典学习、稀疏重建云去除处理的方法。该方法不必要求不同时相的影像数据在时间上比较邻近，也不必要求厚云区域不能重叠，只要多幅不同时相的影像存在足够的互补信息，就能够较好地对厚云区域数据进行修复重建。且在厚云区修复重建的过程中，对多时相影像进行光谱重组，按照多/高光谱影像的思路对之进行处理，以影像间的相关系数衡量其光谱相关性，自适应地根据不同时相影像的相关性考虑其在修复中所起作用的大小，利用稀疏表达理论依据学习的过完备字典最终实现对各个时相厚云区数据的同时填补。本发明提出的利用多时相数据去除光学遥感影像大面积厚云的方法，可以有效利用多时相影像的互补信息，实现大面积厚云的彻底去除，明显提高影像的辐射质量，大幅提升影像的综合应用潜力，在目标识别、地物分类、矿物勘察、军事侦察等方面具有重要的应用价值。因此，光学遥感影像大面积厚云去除方法不仅具有非常重要的学术价值而且具有重要的现实意义。

附图说明

图1是本发明的实施例流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对根据本发明一个实施例的多时相光学遥感影像大面积厚云去除方法进一步详细说明。应当理解，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。

本发明技术方案可采用计算机软件技术实现自动运行流程。以下结合图1详述实施例去除大面积厚云的具体步骤。

步骤1，将所需处理的多时相影像序列进行几何校正，得到同一区域不同时相的影像。

本发明提出的利用多时相数据去除光学遥感影像大面积厚云方法，在对影像处理前，需要进行几何校正，以便得到多时相影像同一区域的数据，具体几何校正为现有技术。实施例的同区域多时相影像是经过几何校正后的400×400的7幅不同时相的MODIS影像。

步骤2，对同一区域不同时相的影像进行厚云区检测，得到云掩膜数据，并计算各时相影像非厚云区的相关系数。

实施例对同一区域7幅不同时相的MODIS影像进行厚云检测。实施例采用可见光第2通道（0.841-0.876μm）和第1通道（0.620-0.670μm）的比值R₂/R₁作为云指数以进行厚云检测，其原理是云在这两个通道上的反射率很接近，而水体和植被在这两个通道上有较大差异，当该比值R₂/R₁介于0.9和1.1之间即认为是厚云，但此法不适用于沙漠地区。具体实施时，本领域技术人员可以根据不同传感器的特点，查找相关资料及凭借经验知识确定其他影像的云指数计算公式以及厚云区云指数的范围。另外，影像中也会存在一些异常值或无效值，处理中可同时将这类值标记为厚云区，以便后续一并进行修复重建。在厚云检测后，得到各幅影像的云掩膜，据此计算各幅影像非厚云区的相关系数（CC），用以衡量各影像的相关性，自适应地确定其在去云修复过程中所起的作用。其中，影像的相关系数是通过公式（1）计算的。

$\mathrm{CC}(x,y)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(x_i-{\mathrm{\μ}}_x)(y_u-{\mathrm{\μ}}_y)}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-{\mathrm{\μ}}_x)}^2\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(y_i-{\mathrm{\μ}}_y)}^2}}}---\left(1\right)$

其中x，y是两幅影像（这里把影像转换为一位数组处理），N是影像的像素个数，x_i、y_i是影像x，y中的第i个像素，i的取值为1,2,…N，μ_x、μ_y分别是其对应的影像均值。

步骤3，将同一区域不同时相的影像按照时间序列组成多维的空间光谱影像，对空间光谱影像划分影像子块，将所有的影像子块重组为二维的矩阵；将步骤2所得各时相影像的云掩膜按照时间序列组成空间光谱云掩膜，对空间光谱云掩膜划分掩膜子块，将所有的掩膜子块重组为二维的矩阵。

实施例中，把7幅影像按照时间序列光谱重组，组成一个新的400×400×7的空间光谱影像，可按从左到右、从上到下每次移动一个像素的滑动方式将其分为157609个4×4×7的影像子块，然后按照列向把影像子块展开成112×1的列向量，最后把所有的列向量以影像滑块的先后顺序组成一个112×157609的二维矩阵，这个由影像子块向量构成的矩阵即是后面需要处理的矩阵。同时，对云掩膜也进行类似的处理，最终得到一个112×157609的掩膜矩阵，即掩膜子块向量构成的矩阵。

步骤4，根据以上的影像子块和掩膜子块的矩阵进行字典学习，重建影像厚云区域。

在上述步骤的基础上，即可进行字典学习。实施例采用KSVD字典学习的方式。KSVD的目的即是在样本Y的基础上，学习出一个符合稀疏标准T₀的字典D，同时得到字典的表示系数X，如公式（2）所示。其中，x_i表示X的第i列元素，||x_i||₀表示0范数，即x_i中非零元素的个数。

$\underset{D,X}{\min }\left\{{\left|\right|Y-\mathrm{DX}\left|\right|}_F^2\right\}s.t.\∀i,{\left|\right|x_i\left|\right|}_0\≤T_0---\left(2\right)$

步骤4中设字典的大小为kn²×m，m≥256且m＞kn²；根据步骤2所得各时相影像非厚云区的相关系数，自适应地确定多时相数据在字典学习过程中的权重。

实施例中，设定字典的大小为112×256，初值为列向归一化的离散余弦变换（DCT）基。目前字典学习的过程主要包括下列两个过程：

①在给定字典的情况下，求解字典的表示系数。求解字典的表示系数，目前通用的是追踪算法，主要有匹配追踪（MP）、正交匹配追踪（OMP）、正则化正交匹配追踪（ROMP）、压缩感知追踪（CoSaMP）、子空间追踪（SP）、基追踪（BP）以及树型追踪（TMP）等，综合算法的精度、复杂度和效率，建议采用正交匹配追踪（OMP）即可。

②在已知字典表示系数的情况下，更新字典。更新字典是按原子（字典的每一列就是一个原子）进行的，即每一次只更新一个原子。如公式（3）更新第k（k＝1,2,…,112）个原子时，首先计算不含该原子时影像表示的残差E_k

$E_k=Y-\underset{j\≠k}{\mathrm{\Σ}}{d}_j{x}^j---\left(3\right)$

其中，Y表示顾及了多时相数据在字典学习中的权重的所有影像子块组成的二维矩阵；d_j表示字典的第j个原子，j=1,2,…112；x^j表示X的第j列元素。

当加入稀疏约束，也就是云掩膜子块矩阵Ω_k的时候，残差变为将其按照公式（4）进行SVD分解，更新后的第k个原子d_k为左酉矩阵U的第一列，其对应的表示系数为右酉矩阵V与奇异值矩阵Λ第一个元素的乘积。

$E_k^R={\mathrm{U\ΛV}}^T---\left(4\right)$

更新字典按照原子的先后顺序依次进行，当所有的原子都遍历更新一次，即完成字典的一次更新。如果要取得更好的实验结果，则需要进行更多次数的字典更新。在①②过程结束后，便可得到7幅影像组成的空间光谱影像在学习的过完备字典D上的稀疏表示。最后，利用该字典及其表示系数的乘积就可以重建厚云区。

本领域普通技术人员可以理解，利用本发明不仅对多时相影像数据的要求不像常规方法那样苛刻，还可以合理地利用多时相影像之间的相关性去除遥感影像中大面积厚云。

本发明提出的上述方法对于去除光学遥感影像中的大面积厚云具有广泛的通用性，受客观因素的限制较少。由模拟实验实际测试结果表明，该方法具有较高的精度，修复后的影像与完好影像之间的相关系数可以达到0.99以上。

应该注意到并理解，在不脱离后附的权利要求所要求的本发明的精神和范围的情况下，能够对上述详细描述的本发明做出各种修改和改进。因此，要求保护的技术方案的范围不受所给出的任何特定示范教导的限制。

Claims (3)

1.一种利用多时相数据去除光学遥感影像大面积厚云的方法，其特征在于，包括以下步骤：步骤1，将所需处理的多时相影像序列进行几何校正，得到同一区域不同时相的影像；

步骤2，对同一区域不同时相的影像进行厚云区检测，得到各时相影像的云掩膜，并计算各时相影像非厚云区的相关系数如下，

$\mathrm{CC}(x,y)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(x_i-{\mathrm{\μ}}_x)(y_i-{\mathrm{\μ}}_y)}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-{\mathrm{\μ}}_x)}^2}\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(y_i-{\mathrm{\μ}}_y)}^2}}$

其中，x、y是两幅影像，N是影像的像素个数，x_i、y_i是影像x、y中的第i个像素，i的取值为1,2,…N，μ_x、μ_y分别是两幅影像对应的影像均值；

步骤3，将同一区域不同时相的影像按照时间序列组成多维的空间光谱影像，对空间光谱影像划分影像子块，将所有的影像子块重组为二维的矩阵；将步骤2所得各时相影像的云掩膜按照时间序列组成空间光谱云掩膜，对空间光谱云掩膜划分掩膜子块，将所有的掩膜子块重组为二维的矩阵；实现方式如下，

对空间光谱影像采用尺寸为n×n的窗口滑动划分成大小为n×n×k的影像子块，其中，k是影像的光谱维数，即同一区域不同时相的影像的个数，然后将每一个影像子块重组为一个列向量，并将所有的列向量按滑动顺序组成一个二维的矩阵；并且，对空间光谱云掩膜采用尺寸为n×n的窗口滑动划分成大小为n×n×k的掩膜子块，然后将每一个掩膜子块重组为一个列向量，并将所有的列向量按滑动顺序组成一个二维的矩阵；

步骤4，根据以上的影像子块和掩膜子块的矩阵进行字典学习，重建影像厚云区域。

2.根据权利要求1所述利用多时相数据去除光学遥感影像大面积厚云的方法，其特征在于：步骤2进行厚云区检测时，将影像中的异常值也标记为厚云区。

3.根据权利要求2所述利用多时相数据去除光学遥感影像大面积厚云的方法，其特征在于：步骤4中设字典的大小为kn²×m，m≥256且m＞kn²；根据步骤2所得各时相影像非厚云区的相关系数，自适应地确定多时相数据在字典学习过程中的权重。