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

摘要

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

主权项

一种基于稀疏表示的遥感图像去云方法，其特征在于，包括以下几个步骤：步骤一、提取云掩膜矩阵，得到新的图像；将厚云提取为一个与原始图像大小相同的掩膜矩阵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)式的限制下，求图像的最稀疏表示变为求解：$\begin{array}{cc}\{{{\alpha }_n}^{\mathrm{opt}},{{\alpha }_t}^{\mathrm{opt}}\}=\arg \underset{\{{\alpha }_n,{\alpha }_t\}}{\min }{\left|\right|{\alpha }_t\left|\right|}_0+{\left|\right|{\alpha }_t\left|\right|}_0& s.t.X=T_t{\alpha }_t+T_n{\alpha }_n---\left(2\right)\end{array}$其中，α_n^opt表示光滑部分所求的表示系数，α_t^opt表示纹理部分所求的表示系数，||·||₀表示l₀范数，表示计算非零项的个数；根据基追踪算法将(1)式转化为求解l₁范数问题，同时考虑图像含有噪声以及对光滑部分进行全变差调整，(2)式的求解变为：$\{{{\alpha }_n}^{\mathrm{opt}},{{\alpha }_t}^{\mathrm{opt}}\}=\arg \underset{\{{\alpha }_n,{\alpha }_t\}}{\min }{\left|\right|{\alpha }_t\left|\right|}_1+{\left|\right|{\alpha }_t\left|\right|}_1+\lambda {{\left|\right|X-T_t{\alpha }_t-T_n{\alpha }_n\left|\right|}_2}^2+{\gamma }^{\mathrm{TV}}\left\{T_n{\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＝T_t⁺X_t，其中T_n⁺为Curvelet反变换，T_t⁺为LDCT反变换，则(3)式变为：$\{{X_n}^{\mathrm{opt}},{X_t}^{\mathrm{opt}}\}=\arg \underset{\{{\alpha }_n,{\alpha }_t\}}{\min }{\left|\right|{{T_n}^{+}X}_n\left|\right|}_1+{\left|\right|{T_t}^{+}{X}_t\left|\right|}_1+\lambda {{\left|\right|X-X_t-x_n\left|\right|}_2}^2+{\gamma }^{\mathrm{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变换系数α_n＝T_t⁺(X_t+R)；使用阈值为δ的软阈值函数处理系数α_t，得到α_t，再利用公式X_t＝T_tα_t重构X_n；其中T_t⁺表示LDCT反变换；所述的软阈值函数为：$\alpha =\left\{\begin{array}{cc}\mathrm{sign}\left(\alpha \right)\left(\right|\alpha |-\delta )& \left|\alpha \right|\ge \delta \\ 0& \left|\alpha \right|<\delta \end{array}\right.$(c)对由(a)得到的图像的光滑部分X_n进行全变差调整，全变差形式由下式给出：$X_n=X_n-\mu \frac{\partial \mathrm{TV}\left\{X_n\right\}}{{\partial X}_n}=X_n-\mu ▿\left(\frac{{▿X}_n}{\left|{▿X}_n\right|}\right)$其中，为梯度算子，为求偏导数符号，在前1000次迭代的过程中，令μ＝2，其他情况，令μ＝1；(d)更新迭代阈值：δ＝δ‑λ(e)判断迭代次数是否小于N，如果小于N，判断迭代阈值和步长之间的关系，如果δ＞λ，返回(a)，对图像继续分解；如果δ＜λ，则对图像分解结束，停止迭代，如果迭代次数等于N，停止迭代，得到图像在字典T_n和字典T_t表示下最终的稀疏系数α_t，α_n；α_t为X在字典T_n表示下的最终的稀疏系数，α_n为X在字典T_t表示下的最终的稀疏系数；步骤三、根据步骤二得到的稀疏系数α_t，α_n，使用公式X＝T_tα_t+T_nα_n重构出修复后的去云图像$\bar{X}=T_t{\alpha }_t+T_n{\alpha }_n$。