CN104156923A - 一种基于稀疏表示的多光谱遥感图像去云方法 - Google Patents

Abstract

本发明涉及一种基于稀疏表示的多光谱遥感图像去云方法，利用多时像方法的同时，利用SAR影像和不同时相的低分辨多光谱影像分别和高分辨多光谱影像的对应图像块建立稀疏字典对，然后利用稀疏编码重建云区部分信息。其中，低分辨光学影像字典提供主要的低频信息，SAR图像字典提供主要的高频信息。本发明相对与传统的仅利用SAR图像去云的方法能得到更接近真实场景的重建图像，在确保低频信息完整的同时更好的保留了图像的高频信息。该方法同时适用于小区域云和大区域云的去除，均能得到较理想的结果。

Description

一种基于稀疏表示的多光谱遥感图像去云方法

技术领域

本发明属于计算机图像的处理，涉及一种多光谱遥感图像的处理，具体涉及一种基于稀疏表示的多光谱遥感图像去云方法。

背景技术

光学影像作为重要的遥感数据源，在资源、气象、环境监测等领域具有重要的应用价值。但光学影像容易受到天气尤其是云雾的影响，根据数据统计每年大概2/3的光学影像都受到云雾的干扰，这使得光学影像丢失了中大量的重要信息，并且给后续图像的判读和分析带来巨大的困难，例如图像识别和分类。因此，运用图像处理技术，去除或削弱云雾的影响，提高遥感数据的利用率具有非常重要的研究背景和实用意义。

传统的遥感图像去云方法主要有基于图像修复、基于多光谱和基于多时相三种方法。基于图像修复的方法是利用含云图像本身的信息来重建云区信息，由于缺少辅助图像的信息使得所得结果不够准确。基于多光谱的方法大多数都是利用多光谱图像中不同波段的信息来重构云区部分。基于多时相的方法主要是利用不同时相的无云遥感图像来预测当前图像的含云区域。

文献“Cloud Removal of Optical Image Using SAR Data for ALOS Application.Experimenting on Simulated ALOS Data”中提出了利用差值技术和SAR图像信息来去除多光谱遥感影像中的云雾区域。但是该方法不具有普适性，而且对于某些特定的目标去云效果较差。

综上所述，传统的利用SAR图像来去除云雾影像的方法具有不稳定性等不足。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于稀疏表示的多光谱遥感图像去云方法，将SAR图像和多时像方法结合，通过基于学习的方法重建含云图像中的云区信息。

技术方案

一种基于稀疏表示的多光谱遥感图像去云方法，其特征在于步骤如下：

步骤1：根据Landsat-5TM影像的特点，利用基于阈值的方法检测给定目标图像I中的像素点，确定像素点的像素值大于阈值TH＝230的像素点为云区的候选点；然后将像素值低于50，且按照光照方向在云区附近大小为原图像尺寸25％的区域内的点为云的阴影点；标记云区的候选点和阴影点为1，其它为0，得到云区标记图像F；

步骤2：提取多个无云的高分辨率多光谱影像HRI、低分辨率多光谱影像LRI和SAR图像的特征得到对应的特征图像，分别得到HRI特征块H＝{x₁,x₂,…,x_N}、LRI特征快L＝{y₁,y₂,…,y_N}和SAR特征块S＝{z₁,z₂,…,z_N}；其中，LRI通过双线性差值的方法放大到与HRI具有相同大小，N表示图像块的数量，x_i、y_i和z_i，分别表示三种图像中大小的图像块拉直后的列向量，i＝1,2,…,N；

利用HRI和LRI特征块学习得到稀疏字典D_HL

${D_{\mathrm{HL}}}^{*}=\underset{D_{\mathrm{HL}},}{\arg }\underset{A_{\mathrm{HL}}}{\min }\{{\left|\right|X_{\mathrm{HL}}-D_{\mathrm{HL}}{A}_{\mathrm{HL}}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|A_{\mathrm{HL}}\left|\right|}_1\}$

利用HRI和SAR特征块学习得到稀疏字典D_HS

${D_{\mathrm{HS}}}^{*}=\underset{D_{\mathrm{HS}},}{\arg }\underset{A_{\mathrm{HS}}}{\min }\{{\left|\right|X_{\mathrm{HS}}-D_{\mathrm{HS}}{A}_{\mathrm{HS}}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|A_{\mathrm{HS}}\left|\right|}_1\}$

其中：X_HL＝[H；L]，D_HL＝[D_H1；D_L]，A_HL＝[A_H；A_L]是X_HL的稀疏表示系数，X_HS＝[H；S]，D_HS＝[D_H2；D_S]，A_HS＝[A_H；A_S]是X_HS的稀疏表示系数；

步骤3：将云区标识图F中各个云区以最小矩形标识，对每个标识矩形从外向内按顺时针方向依次取相同大小的图像块并拉直后的图像为待求图像块x，同时在对应的LRI和SAR中取相应位置的图像特征块并拉直，记为y，z；

计算y，z分别在D_HL和D_HS下的稀疏表示系数：

${\mathrm{\α}}_1=\underset{{\mathrm{\α}}_1}{\min }{\left|\right|y-D_L{\mathrm{\α}}_1\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\α}}_1\left|\right|}_1$

${\mathrm{\α}}_2=\underset{{\mathrm{\α}}_2}{\min }{\left|\right|y-D_S{\mathrm{\α}}_2\left|\right|}_2^2+\mathrm{\λ}{\left|\right|{\mathrm{\α}}_2\left|\right|}_1$

利用LRI预测得到的图像块x₁＝D_H1·α₁

利用SAR预测得到的图像块x₂＝D_H2·α₂

继而得到待求图像块x^*＝x₁+x₂＝D_H1·α₁+D_H2·α₂；

步骤4：依次用重建后的图像块替换图像I中对应的图像块，得到去云后的图像I′。

在步骤3完成后，利用调整重建图像块的灰度值，然后再用调整灰度值后的图像块替换图像I中对应的图像块，得到去云后的图像I′；其中：/Ω表示图像I中的非云区域，为有效数据区域；Ω为含云区域。

：所述高分辨率多光谱影像HRI、低分辨率多光谱影像LRI和SAR图像块的大小为5×5。

所述两个稀疏字典对D_HL或D_HS的大小为512。

所述训练样本图像块数为10000。

有益效果

本发明提出的一种基于稀疏表示的多光谱遥感图像去云方法，利用多时像方法的同时，利用SAR影像和不同时相的低分辨多光谱影像分别和高分辨多光谱影像的对应图像块建立稀疏字典对，然后利用稀疏编码重建云区部分信息。其中，低分辨光学影像字典提供主要的低频信息，SAR图像字典提供主要的高频信息。

本发明相对与传统的仅利用SAR图像去云的方法能得到更接近真实场景的重建图像，在确保低频信息完整的同时更好的保留了图像的高频信息。该方法同时适用于小区域云和大区域云的去除，均能得到较理想的结果。

具体实施方式

现结合实施例对本发明作进一步描述：

合成孔径雷达(Synthetic Aperture Radar,SAR)是20世纪50年代末研制成果的一种主动微波传感器，也是微波传感器中发展最为迅速和最有成效的传感器之一。作为一种主动传感器，它能够不受光照和天气条件的限制实现全天时、全天候对地观测，还可以透过地表和植被获取地表以下信息。因此SAR影像是一种无云图像。

传统的稀疏表示理论

在稀疏表示中，矩阵D∈R^n×m(n<m)表示拥有m个原子的过完备字典，其中每一个列表示一个基本原子。对于一个给定的信号x，可以表示为D中若干个原子的一个稀疏的线性组合，也就是说有很少的几个非零系数。寻找x的最稀疏的表示问题可以表示为：

$\widehat{\mathrm{\&alpha;}}=\arg \min {\left|\right|\mathrm{\&alpha;}\left|\right|}_{0}s.t.{\left|\right|x-\mathrm{D\&alpha;}\left|\right|}_2^2<\mathrm{\&epsiv;}---\left(1\right)$

这是一个NP难的问题。但是当稀疏表示系数α足够的稀疏时，该问题可以转换为求解l₁范数的问题：

$\widehat{\mathrm{\&alpha;}}=\arg \min {\left|\right|\mathrm{\&alpha;}\left|\right|}_{1}s.t.{\left|\right|x-\mathrm{D\&alpha;}\left|\right|}_2^2<\mathrm{\&epsiv;}---\left(2\right)$

利用Lagrange正则化，上式可表示为：

$\widehat{\mathrm{\&alpha;}}=\arg \min {{\left|\right|x-\mathrm{D\&alpha;}\left|\right|}_2^2+\mathrm{\&lambda;}\left|\right|\mathrm{\&alpha;}\left|\right|}_1---\left(3\right)$

其中，λ是正则化参数用来平衡x的精度和稀疏度。

基于稀疏表示的图像去云

(1)根据Landsat-5TM影像的特点，对于给定目标图像I，利用基于阈值的方法确定候选云区，然后通过光照方向和云区位置信息判断云阴影区域。得到云区标记图像F。

(2)通过大量无云的高、低分辨率多光谱影像和SAR图像，提取低分辨多光谱影像和SAR影像的特征得到对应的特征图像。利用高分辨率多光谱影像图像块和对应的低分辨率多光谱影像特征块、高分辨率多光谱影像图像块和对应的SAR影像特征块分别通过学习的方法建立两个稀疏字典D_HL和D_HS。

(3)针对目标图像I，利用云区标识图F，将各个云区以最小矩形标识。对每个标识矩形从外向内按顺时针方向依次取相同大小的图像块并拉直，记为x(x是待求图像块)，同时在对应的LRI和SAR中取相应位置的图像特征块并拉直，记为y，z。计算y，z分别在D_HL和D_HS下的稀疏表示系数，利用稀疏编码重建含云区域信息。

具体实施例：

(1)云区和阴影检测

由于在Landsat-5TM影像中云区部分具有较高的亮度值，因此可以利用基于阈值的方法来检测云区。对于给定含云图像I，逐个判断I中的像素点，若当前点的像素值大于给定阈值TH＝230时，认为该点是云区的候选点；否则认为不是云区。但是由于其他地物的干扰，云区候选点中存在许多孤立的像素点，因此需要去除这些孤立的干扰点以确定云区。当云区像素点被确定后，云的阴影部分可以通过云的位置和光照方向来确定，将像素值低于50，且按照光照方向在云区附近大小为原图像尺寸25％的区域内的点为云的阴影点；以此方法得到云区和阴影的标记图F，其中云和阴影像素点标识值为1，其它为0。

(2)字典对学习

定义高、低分辨多光谱无云影像为HRI和LRI，SAR影像为SAR，HRI图像块、LRI特征快和SAR特征块分别表示为H＝{x₁,x₂,…,x_N}，L＝{y₁,y₂,…,y_N}和S＝{z₁,z₂,…,z_N}。其中，LRI通过双线性差值的方法放大到与HRI具有相同大小，N表示图像块的数量，x_i、y_i和z_i(i＝1,2,…,N)分别表示三种图像中大小的图像块拉直后的列向量。

我们利用HRI和LRI特征块，HRI和SAR特征块分别学习两个稀疏字典对D_HL和D_HS。根据公式(3)可利用下式学习得到D_HL和D_HS：

${D_{\mathrm{HL}}}^{*}=\underset{D_{\mathrm{HL}},}{\arg }\underset{A_{\mathrm{HL}}}{\min }\{{\left|\right|X_{\mathrm{HL}}-D_{\mathrm{HL}}{A}_{\mathrm{HL}}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|A_{\mathrm{HL}}\left|\right|}_1\}---\left(4\right)$

${D_{\mathrm{HS}}}^{*}=\underset{D_{\mathrm{HS}},}{\arg }\underset{A_{\mathrm{HS}}}{\min }\{{\left|\right|X_{\mathrm{HS}}-D_{\mathrm{HS}}{A}_{\mathrm{HS}}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|A_{\mathrm{HS}}\left|\right|}_1\}---\left(5\right)$

其中，X_HL＝[H；L]，D_HL＝[D_H1；D_L]，A_HL＝[A_H；A_L]是X_HL的稀疏表示系数，X_HS＝[H；S]，D_HS＝[D_H2；D_S]，A_HS＝[A_H；A_S]是X_HS的稀疏表示系数。此处，三种图像块的大小均为5×5，字典大小为512，训练样本图像块数为10000。

为了使训练简单有效，采用下述4个算子来提取LRI和SAR的特征：

f₁＝[-1，0，1]，f₂＝f₁ ^T  (6)

f₃＝[1，0，-2，0，1]，f₄＝f₃ ^T

(3)云区信息重建

对于给定的含云图像I，利用云区标识图F，将各个云区以最小矩形标识。对每个标识矩形从外向内按顺时针方向依次取大小的图像块并拉直，记为x(x是待求图像块)，同时在对应的LRI和SAR中取相应位置的图像特征块并拉直，记为y，z。

对当前y，z，利用上步训练得到的字典对D_HS和D_HL，计算稀疏表示系数α₁和α₂：

${\mathrm{\&alpha;}}_1=\underset{{\mathrm{\&alpha;}}_1}{\min }{\left|\right|y-D_L{\mathrm{\&alpha;}}_1\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|{\mathrm{\&alpha;}}_1\left|\right|}_1---\left(7\right)$

${\mathrm{\&alpha;}}_2=\underset{{\mathrm{\&alpha;}}_2}{\min }{\left|\right|y-D_S{\mathrm{\&alpha;}}_2\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|{\mathrm{\&alpha;}}_2\left|\right|}_1---\left(8\right)$

然后利用下式计算利用LRI预测得到的图像块x₁和用SAR预测得到的图像块x₂：

x₁＝D_H1·α₁  (9)

x₂＝D_H2·α₂  (10)

根据公式(9)和公式(10)得到预测图像块x^*：

x^*＝x₁+x₂＝D_H1·α₁+D_H2·α₂  (11)

由于重建后的图像块和原始图像块有一定的灰度差异，我们定义图像I中的有效数据区域(非云区域)为/Ω，含云区域为Ω。利用下式调整重建图像块的灰度值：

$\hat{x}=x^{*}\&times;\frac{{\left|\right|x^{/\mathrm{\&Omega;}}\left|\right|}^2}{{\left|\right|x^{*\mathrm{\&Omega;}}\left|\right|}^2}---\left(12\right)$

依次用重建后的图像块替换图像I中对应的图像块，即可得到去云后的图像I′。

Claims (5)

1.一种基于稀疏表示的多光谱遥感图像去云方法，其特征在于步骤如下：

步骤1：根据Landsat-5TM影像的特点，利用基于阈值的方法检测给定目标图像I中的像素点，确定像素点的像素值大于阈值TH＝230的像素点为云区的候选点；然后将像素值低于50，且按照光照方向在云区附近大小为原图像尺寸25％的区域内的点为云的阴影点；标记云区的候选点和阴影点为1，其它为0，得到云区标记图像F；

步骤2：提取多个无云的高分辨率多光谱影像HRI、低分辨率多光谱影像LRI和SAR图像的特征得到对应的特征图像，分别得到HRI特征块H＝{x₁,x₂,…，x_N}、LRI特征快L＝{y₁,y₂,…,y_N}和SAR特征块S＝{z₁,z₂,…,z_N}；其中，LRI通过双线性差值的方法放大到与HRI具有相同大小，N表示图像块的数量，x_i、y_i和z_i，分别表示三种图像中大小的图像块拉直后的列向量，i＝1,2,…,N；

利用HRI和LRI特征块学习得到稀疏字典D_HL

${D_{\mathrm{HL}}}^{*}=\underset{D_{\mathrm{HL}},}{\arg }\underset{A_{\mathrm{HL}}}{\min }\{{\left|\right|X_{\mathrm{HL}}-D_{\mathrm{HL}}{A}_{\mathrm{HL}}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|A_{\mathrm{HL}}\left|\right|}_1\}$

利用HRI和SAR特征块学习得到稀疏字典D_HS

${D_{\mathrm{HS}}}^{*}=\underset{D_{\mathrm{HS}},}{\arg }\underset{A_{\mathrm{HS}}}{\min }\{{\left|\right|X_{\mathrm{HS}}-D_{\mathrm{HS}}{A}_{\mathrm{HS}}\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|A_{\mathrm{HS}}\left|\right|}_1\}$

其中：X_HL＝[H；L]，D_HL＝[D_H1；D_L]，A_HL＝[A_H；A_L]是X_HL的稀疏表示系数，X_HS＝[H；S]，D_HS＝[D_H2；D_S]，A_HS＝[A_H；A_S]是X_HS的稀疏表示系数；

步骤3：将云区标识图F中各个云区以最小矩形标识，对每个标识矩形从外向内按顺时针方向依次取相同大小的图像块并拉直后的图像为待求图像块x，同时在对应的LRI和SAR中取相应位置的图像特征块并拉直，记为y，z；

计算y，z分别在D_HL和D_HS下的稀疏表示系数：

${\mathrm{\&alpha;}}_1=\underset{{\mathrm{\&alpha;}}_1}{\min }{\left|\right|y-D_L{\mathrm{\&alpha;}}_1\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|{\mathrm{\&alpha;}}_1\left|\right|}_1$

${\mathrm{\&alpha;}}_2=\underset{{\mathrm{\&alpha;}}_2}{\min }{\left|\right|y-D_S{\mathrm{\&alpha;}}_2\left|\right|}_2^2+\mathrm{\&lambda;}{\left|\right|{\mathrm{\&alpha;}}_2\left|\right|}_1$

利用LRI预测得到的图像块x₁＝D_H1·α₁

利用SAR预测得到的图像块x₂＝D_H2·α₂

继而得到待求图像块x^*＝x₁+x₂＝D_H1·α₁+D_H2·α₂；

步骤4：依次用重建后的图像块替换图像I中对应的图像块，得到去云后的图像I′。

2.根据权利要求1所述基于稀疏表示的多光谱遥感图像去云方法，其特征在于：在步骤3完成后，利用调整重建图像块的灰度值，然后再用调整灰度值后的图像块替换图像I中对应的图像块，得到去云后的图像I′；其中：/Ω表示图像I中的非云区域，为有效数据区域；Ω为含云区域。

3.根据权利要求1所述基于稀疏表示的多光谱遥感图像去云方法，其特征在于：所述高分辨率多光谱影像HRI、低分辨率多光谱影像LRI和SAR图像块的大小为5×5。

4.根据权利要求1所述基于稀疏表示的多光谱遥感图像去云方法，其特征在于：所述两个稀疏字典对D_HL或D_HS的大小为512。

5.根据权利要求1所述基于稀疏表示的多光谱遥感图像去云方法，其特征在于：所述训练样本图像块数为10000。