CN102800058A - 基于稀疏表示的遥感图像去云方法 - Google Patents

Abstract

本发明公开了一种基稀疏表示的遥感图像去云方法，具体涉及一种基于稀疏表示的遥感图像去云方法，包括：提取云掩膜矩阵得到新的图像、初始化迭代参数、固定图像的纹理部分更新图像的光滑部分、固定图像的光滑部分更新图像的纹理部分、使用全变差函数对光滑部分进行调整、更新迭代阈值、判断分解是否完成几个步骤。本发明提出一种基于稀疏表示的遥感图像去云方法，对图像的全变差调整函数的参数进行调整，即在使用块坐标松弛算法来分解图像时，通过改变全变差调整函数的参数来调节算法的收敛速度，提高算法的效率和对图像的分解效果。在尽可能不破坏图像原始信息的前提下，该方法去除厚云的效果比较明显。

Description

基于稀疏表示的遥感图像去云方法

技术领域

本发明属于遥感图像去云技术领域，具体涉及一种基于稀疏表示的遥感图像去云方法。

背景技术

遥感图像已在军事侦察、地质解译、石油勘探、天气预报、农作物长势检测等诸多领域得到了广泛的应用。目前通过遥感获得的图像大部分是光学图像，而光学图像的质量极易受到气候因素的影响，云覆盖就是其中一种。从遥感物理的角度来说，云属于大气气溶胶的范畴，是地球大气中具有一定稳定性的，沉降速度小的，尺度范围在10^-3μm~10μm之间的液态离子或固态离子所组成的混合物集合体。云覆盖不仅影响遥感图像的解译精度，而且也是造成数据缺乏的一个重要因素。尤其是在多云多雾地区，云覆盖使得遥感图像模糊、重要地物信息被覆盖，有用的遥感信息变得匮乏。因此如何有效地去除或减少云的影响，不仅是对遥感图像进行准确解译的基础，也是增加遥感数据有效性的一个重要途径，具有十分重要的实际意义。

去云处理是图像预处理中的一个很重要的环节。在去云方面，除合成孔径雷达传感器能够穿透云层获取信息外，其它传感器均未能彻底解决影像数据的云覆盖问题。因此，从九十年代开始，遥感图像的去云技术开始成为一个活跃的研究领域，在这一发展过程中已经提出了许多算法和思路，其中具有代表性的有：多光谱图像去云法、多传感器图像融合去云法、同态滤波去云法等。

多光谱图像去云法是采用一种对云较敏感的传感器，专门探测云的信息，然后在原始图上减去云图，得到去云以后的图像，或者是不添加任何专门的传感器，而是利用多光谱影像中某些波段对云的较强敏感性来提取云的信息。但是该方法要求的前提很高，需要有冗余的波段，对其应用场合有诸多限制。

多传感器图像融合去云法的原理是在现有条件下，利用不同传感器在不同时间获取的数据，对有云层覆盖地区的影像进行替换，以消除云覆盖的影响。但是利用该方法去云需要解决以下几个问题：首先要解决云覆盖区域多传感器的配准问题，即为了使替补图像能精确的替补有云的原始图像，必须对两幅图像进行几何纠正；其次要解决可能存在的辐射差异问题；而且要保证两幅图像在相同的区域不能同时有云。

同态滤波去云法认为遥感图像中云覆盖所引起的噪声提高了图像的低频信息，而削弱了图像的高频信息，因此可以通过使用一个适当的滤波器来将云从图像中去除。虽然这样的处理能够去除一定的云覆盖，但是对没有云覆盖的清晰区域也会造成一定的影响，而且该方法对厚云的处理效果不理想。

以上算法大部分都只能对比较薄的云进行去除，对厚云的去除效果很不理想，而实际中云的薄厚具有随机性，所以对厚云的去除是一种亟待解决的问题。

由于通过变换来表示信号(图像)在许多应用中具有简单性和有效性而一直受到研究人员的重视，传统的信号(图像)表示理论大多基于非冗余的正交基函数的变换，如傅里叶变换等。最近几年，在正交小波变换的基础上，又提出了许多新的变换，如曲波(Curvelet)、轮廓波(Contourlet)等变换。基于这些变换，普遍采用的是超完备(overcomplete)冗余表示，其基本思想是基函数用称之为字典的超完备的冗余函数来取代，字典中的元素称为原子，信号(图像)由原子的线性组合来表示。其中原子的个数比信号(图像)的维数大得多，由此产生了系数的冗余。稀疏表示就是用最小的系数来表示信号(图像)。比如申请号为201010522273.4的发明专利介绍了利用Curvelet冗余字典稀疏表示图像的方法，但是这种方法只使用了一种冗余字典，没有很好的利用不同字典的对图像特定部分的表示特点。

发明内容

本发明的目的是是为了解决上述现有技术中存在的问题，提出一种基于稀疏表示的遥感图像去云方法，在尽可能不破坏图像原始信息的前提下，去除遥感图像中的厚云。

本发明的基于稀疏表示的遥感图像去云方法，主要包括以下步骤：

步骤一、提取云掩膜矩阵，得到新的图像；

步骤二、对步骤一得到的新的图像X进行优化，即用块坐标松弛算法将图像X稀疏分解为光滑部分T_n和纹理部分T_t，获取最终的稀疏系数α_t，α_n；

步骤三、根据步骤二得到的稀疏系数α_t，α_n，使用公式X＝T_tα_t+T_nα_n重构出修复后的去云图像

本发明的优点在于：

（1）本发明利用稀疏表示的方法，将遥感图像中被云覆盖的区域用字典中的原子进行填充，不局限于传统的只能处理薄云的方法；

（2）通过对全变差函数的参数进行调整，使图像分解的速度加快，提高了分解的质量。

附图说明

图1本发明的去云算法流程图；

图2本发明使用的带云的遥感图像；

图3本发明的去云效果图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明的基于稀疏表示的遥感图像去云方法，对图像的全变差调整函数的参数进行调整，即在使用块坐标松弛算法来分解图像时，通过改变全变差调整函数的参数来调节算法的收敛速度，提高算法的效率。

本发明是一种基于稀疏表示的遥感图像去云方法，流程如图1所示，包括以下几个步骤：

步骤一、提取云掩膜矩阵，得到新的图像。

将厚云提取为一个与原始图像大小相同的掩膜矩阵M，掩膜矩阵M的元素，在有厚云遮挡的地方置为“0”，在没有厚云遮挡的地方置为“1”。设原始的图像为

则将原始图像与掩膜矩阵M对应位置相乘，得到新的图像为X，即

从而直接将含云的区域作为丢失像素的区域来处理，然后对该区域用稀疏字典的原子进行修复，其中修复后的图像，也就是进行去云处理的图像。

步骤二、对步骤一得到的新的图像X进行优化，即用块坐标松弛算法将图像X稀疏分解为光滑部分T_n和纹理部分T_t，获取最终的稀疏系数α_t，α_n。

采用Curvelet变换字典和离散余弦变换(LDCT，Local Discrete Cosine Transform)字典，即分别由Curvelet变换和离散余弦变换所产生的两组过完备正交基来表示图像的纹理部分和光滑部分，其中每一个基中的向量，也称为一个原子。Curvelet变换理论是由Candes于1999年提出，它由脊波理论衍生而来。Curvelet变换克服了脊波变换固定尺度的缺陷，对曲线状奇异特征具有稀疏的表示。LDCT适合于对图像平滑部分或周期部分的稀疏表示。由于该变换主要采用重叠分块的方法，具有较高的冗余度，因此LDCT是一种超完备变换。

根据Meyer的观点，图像可以分为分片光滑部分(以下记为光滑部分)和纹理部分。利用上面所述的两种字典对掩膜后的图像稀疏分解为光滑部分和纹理部分，即：

X＝X_n+X_t＝T_tα_t+T_nα_n    (1)

其中：X_n表示光滑部分，X_t表示纹理部分，T_t为LDCT字典，T_n为Curvelet字典，α_n为X在字典T_n表示下的稀疏系数，α_t为X在字典T_t表示下的稀疏系数。

所以在(1)式的限制下，求图像的最稀疏表示变为求解：

$\{{{\mathrm{\α}}_n}^{\mathrm{opt}},{{\mathrm{\α}}_t}^{\mathrm{opt}}\}=\arg \underset{\{{\mathrm{\α}}_n,{\mathrm{\α}}_t\}}{\min }{\left|\right|{\mathrm{\α}}_t\left|\right|}_0+{\left|\right|{\mathrm{\α}}_t\left|\right|}_{0}s.t.X=T_t{\mathrm{\α}}_t+T_n{\mathrm{\α}}_n---\left(2\right)$

其中，α_n  ^opt表示光滑部分所求的表示系数，α_t  ^opt表示纹理部分所求的表示系数，||·||₀表示l₀范数（0范数），表示计算非零项的个数。(2)式所描述的问题是一个非凸和高度非光滑的，即对(2)式的求解在优化问题中是一个NP难问题，很难进行数值求解，因此根据基追踪算法将其转化为求解l₁范数问题。而对于l₁范数（1范数）问题，很容易通过数值计算的方法来解决。同时考虑图像含有噪声以及对光滑部分进行全变差调整，(2)式的求解变为：

$\{{{\mathrm{\α}}_n}^{\mathrm{opt}},{{\mathrm{\α}}_t}^{\mathrm{opt}}\}=\arg \underset{\{{\mathrm{\α}}_n,{\mathrm{\α}}_t\}}{\min }{\left|\right|{\mathrm{\α}}_t\left|\right|}_1+{\left|\right|{\mathrm{\α}}_t\left|\right|}_1+\mathrm{\λ}{{\left|\right|X-T_t{\mathrm{\α}}_t-T_n{\mathrm{\α}}_n\left|\right|}_2}^2+\mathrm{\γTV}\left\{T_n{\mathrm{\α}}_n\right\}---\left(3\right)$

其中，λ为拉格朗日乘子，γ为正则化参数，TV{T_nα_n}表示对图像的光滑部分进行全变差(TV，Total Variation)调整的部分；

由于X_n＝T_nα_n，X_t＝T_tα_t，所以α_n＝T_n ⁺X_n，

其中T_n ⁺为Curvelet反变换，

为LDCT反变换，则(3)式变为：

$\{{X_n}^{\mathrm{opt}},{X_t}^{\mathrm{opt}}\}=\arg \underset{\{{\mathrm{\α}}_n,{\mathrm{\α}}_t\}}{\min }{\left|\right|T_n^{+}+X_n\left|\right|}_1+{\left|\right|T_t^{+}{X}_t\left|\right|}_1+\mathrm{\λ}{{\left|\right|X-X_t-X_n\left|\right|}_2}^2+\mathrm{\γTV}\left\{X_n\right\}---\left(4\right)$

传统的TV图像修复模型通过最小化能量泛函完成对缺损区域的修复，它可以保持尖锐的边缘且数值实现非常简单，但是该算法对参数的选取较敏感，且运算量大。

对(4)式，采用块坐标松弛算法(BCR，block-coordinate relaxation)求解，其主要思路是对变换系数按块进行更新，以达到对图像的分解，具体包括：

迭代次数N的取值范围可以为[300,3000]，迭代阈值一般固定为δ＝255，步长λ为迭代次数N与迭代阈值δ的函数，服从公式λ＝δ/(N-1)，正则化参数γ∈[0.5,2]。本发明设置迭代次数N＝1500，步长λ＝δ/(N-1)＝255/1499，正则参数γ＝0.5，并初始化光滑部分X_n＝X，纹理部分X_t＝0，并按以下步骤迭代N次；

(a)固定纹理部分X_t不变，更新光滑部分X_n。

先计算残差R＝X-X_n-X_t，然后计算X_n+R的curvelet变换系数α_n＝T_n ⁺(X_n+R)。使用阈值为δ的软阈值函数处理系数α_n，得到α_n，再利用公式X_n＝T_nα_n重构X_n。其中T_n ⁺表示curvelet反变换；

(b)固定由(a)得到的光滑部分X_n不变，更新纹理部分X_t。

先计算残差R＝X-X_n-X_t，然后计算X_t+R的LDCT变换系数

使用阈值为δ的软阈值函数处理系数α_t，得到α_t，再利用公式X_t＝T_tα_t重构X_n。其中

表示LDCT反变换。

所述的软阈值函数为：

$\mathrm{\&alpha;}=\left\{\begin{array}{cc}\mathrm{sign}\left(\mathrm{\&alpha;}\right)\left(\right|\mathrm{\&alpha;}|-\mathrm{\&delta;})& \left|\mathrm{\&alpha;}\right|\&GreaterEqual;\mathrm{\&delta;}\\ 0& \left|\mathrm{\&alpha;}\right|<\mathrm{\&delta;}\end{array}\right.$

(c)对由(a)得到的图像的光滑部分X_n进行全变差调整，对应式（4）中的γTV{X_n}项，确保X_n尽可能光滑，以保证由(4)式得到的稀疏系数尽可能最稀疏。全变差形式由下式给出：

$X_n=X_n-\mathrm{\&mu;}\frac{\&PartialD;\mathrm{TV}\left\{X_n\right\}}{{\&PartialD;X}_n}=X_n-\mathrm{\&mu;}\&dtri;\left(\frac{{\&dtri;X}_n}{\left|{\&dtri;X}_n\right|}\right)$

其中，

为梯度算子，

为求偏导数符号，在前1000次迭代的过程中，令μ=2，其他情况，令μ=1；

(d)更新迭代阈值

δ＝δ-λ

(e)判断迭代次数是否小于N，如果小于N，判断迭代阈值和步长之间的关系，如果δ＞λ，返回(a)，对图像继续分解；如果δ＜λ，则对图像分解结束，停止迭代，如果迭代次数等于N，停止迭代，得到图像在上述两种字典表示下最终的稀疏系数α_t，α_n；α_t为X在字典T_n表示下的最终的稀疏系数，α_n为X在字典T_t表示下的最终的稀疏系数。步骤三、根据步骤二得到的稀疏系数α_t，α_n，，使用公式X＝T_tα_t+T_nα_n重构出修复后的去云图像

$\stackrel{\&OverBar;}{X}=T_t{\mathrm{\&alpha;}}_t+T_n{\mathrm{\&alpha;}}_n$

带云的遥感图像如图2所示，应用本发明的去云方法，得到的去云效果如图3所示，从图中可以看出，去云效果很好。

Claims (4)

1.一种基于稀疏表示的遥感图像去云方法，其特征在于，包括以下几个步骤：

步骤一、提取云掩膜矩阵，得到新的图像；

将厚云提取为一个与原始图像大小相同的掩膜矩阵M，掩膜矩阵M的元素，在有厚云遮挡的地方置为“0”，在没有厚云遮挡的地方置为“1”；设原始的图像为

则将原始图像与掩膜矩阵M对应位置相乘，得到新的图像为X，即

步骤二、对步骤一得到的新的图像X进行优化，即用块坐标松弛算法将图像X稀疏分解为光滑部分T_n和纹理部分T_t，获取最终的稀疏系数α_t，α_n；

利用采用Curvelet变换字典和离散余弦变换字典对掩膜后的图像稀疏分解为光滑部分和纹理部分，即：

X＝X_n+X_t＝T_tα_t+T_nα_n    (1)

其中：X_n表示光滑部分，X_t表示纹理部分，T_t为LDCT字典，T_n为Curvelet字典，α_n为X在字典T_n表示下的稀疏系数，α_t为X在字典T_t表示下的稀疏系数；

在(1)式的限制下，求图像的最稀疏表示变为求解：

$\{{{\mathrm{\&alpha;}}_n}^{\mathrm{opt}},{{\mathrm{\&alpha;}}_t}^{\mathrm{opt}}\}=\arg \underset{\{{\mathrm{\&alpha;}}_n,{\mathrm{\&alpha;}}_t\}}{\min }{\left|\right|{\mathrm{\&alpha;}}_t\left|\right|}_0+{\left|\right|{\mathrm{\&alpha;}}_t\left|\right|}_{0}s.t.X=T_t{\mathrm{\&alpha;}}_t+T_n{\mathrm{\&alpha;}}_n---\left(2\right)$

其中，α_n ^opt表示光滑部分所求的表示系数，α_t ^opt表示纹理部分所求的表示系数，||·||₀表示l0范数，表示计算非零项的个数；根据基追踪算法将(1)式转化为求解l1范数问题，同时考虑图像含有噪声以及对光滑部分进行全变差调整，(2)式的求解变为：

$\{{{\mathrm{\&alpha;}}_n}^{\mathrm{opt}},{{\mathrm{\&alpha;}}_t}^{\mathrm{opt}}\}=\arg \underset{\{{\mathrm{\&alpha;}}_n,{\mathrm{\&alpha;}}_t\}}{\min }{\left|\right|{\mathrm{\&alpha;}}_t\left|\right|}_1+{\left|\right|{\mathrm{\&alpha;}}_t\left|\right|}_1+\mathrm{\&lambda;}{{\left|\right|X-T_t{\mathrm{\&alpha;}}_t-T_n{\mathrm{\&alpha;}}_n\left|\right|}_2}^2+\mathrm{\&gamma;TV}\left\{T_n{\mathrm{\&alpha;}}_n\right\}---\left(3\right)$

其中，λ为拉格朗日乘子，γ为正则化参数，TV{T_nα_n}表示对图像的光滑部分进行全变差调整的部分；

由于X_n＝T_nα_n，X_t＝T_tα_t，所以α_n＝T_n ⁺X_n，

其中T_n ⁺为Curvelet反变换，

为LDCT反变换，则(3)式变为：

$\{{X_n}^{\mathrm{opt}},{X_t}^{\mathrm{opt}}\}=\arg \underset{\{{\mathrm{\&alpha;}}_n,{\mathrm{\&alpha;}}_t\}}{\min }{\left|\right|T_n^{+}+X_n\left|\right|}_1+{\left|\right|T_t^{+}{X}_t\left|\right|}_1+\mathrm{\&lambda;}{{\left|\right|X-X_t-X_n\left|\right|}_2}^2+\mathrm{\&gamma;TV}\left\{X_n\right\}---\left(4\right)$

对(4)式，采用块坐标松弛算法求解，具体包括：

设置迭代次数N、迭代阈值δ，步长λ，λ＝δ/(N-1)，正则化参数γ，初始化光滑部分X_n＝X，纹理部分X_t＝0，并按以下步骤进行迭代；

(a)固定纹理部分X_t不变，更新光滑部分X_n；

先计算残差R＝X-X_n-X_t，然后计算X_n+R的curvelet变换系数α_n＝T_n ⁺(X_n+R)；使用阈值为δ的软阈值函数处理系数α_n，得到α_n，再利用公式X_n＝T_nα_n重构X_n；其中T_n ⁺表示curvelet反变换；

(b)固定由(a)得到的光滑部分X_n不变，更新纹理部分X_t；

先计算残差R＝X-X_n-X_t，然后计算X_t+R的LDCT变换系数

使用阈值为δ的软阈值函数处理系数α_t，得到α_t，再利用公式X_t＝T_tα_t重构X_n；其中

表示LDCT反变换；

所述的软阈值函数为：

$\mathrm{\&alpha;}=\left\{\begin{array}{cc}\mathrm{sign}\left(\mathrm{\&alpha;}\right)\left(\right|\mathrm{\&alpha;}|-\mathrm{\&delta;})& \left|\mathrm{\&alpha;}\right|\&GreaterEqual;\mathrm{\&delta;}\\ 0& \left|\mathrm{\&alpha;}\right|<\mathrm{\&delta;}\end{array}\right.$

(c)对由(a)得到的图像的光滑部分X_n进行全变差调整，全变差形式由下式给出：

$X_n=X_n-\mathrm{\&mu;}\frac{\&PartialD;\mathrm{TV}\left\{X_n\right\}}{{\&PartialD;X}_n}=X_n-\mathrm{\&mu;}\&dtri;\left(\frac{{\&dtri;X}_n}{\left|{\&dtri;X}_n\right|}\right)$

其中，

为梯度算子，

为求偏导数符号，在前1000次迭代的过程中，令μ=2，其他情况，令μ=1；

(d)更新迭代阈值：

δ＝δ-λ

(e)判断迭代次数是否小于N，如果小于N，判断迭代阈值和步长之间的关系，如果δ＞λ，返回(a)，对图像继续分解；如果δ＜λ，则对图像分解结束，停止迭代，如果迭代次数等于N，停止迭代，得到图像在上述两种字典表示下最终的稀疏系数α_t，α_n；

α_t为X在字典T_n表示下的最终的稀疏系数，α_n为X在字典T_t表示下的最终的稀疏系数；

步骤三、根据步骤二得到的稀疏系数α_t，α_n，使用公式X＝T_tα_t+T_nα_n重构出修复后的去云图像

$\stackrel{\&OverBar;}{X}=T_t{\mathrm{\&alpha;}}_t+T_n{\mathrm{\&alpha;}}_n.$

2.根据权利要求1所述的一种基于稀疏表示的遥感图像去云方法，其特征在于，所述的步骤二中，迭代次数N的取值范围为[300,3000]。

3.根据权利要求1所述的一种基于稀疏表示的遥感图像去云方法，其特征在于，所述的步骤二中，迭代阈值δ＝255。

4.根据权利要求1所述的一种基于稀疏表示的遥感图像去云方法，其特征在于，所述的步骤二中，正则化参数γ∈[0.5,2]。

Images