CN103020912B - 一种结合波段聚类和稀疏表达的遥感影像复原方法 - Google Patents

Abstract

本发明提供一种结合波段聚类和稀疏表达的遥感影像复原方法，以提高高光谱遥感影像的空间分辨率为目的，针对高光谱影像光谱维信息丰富以及不同波段噪声强度不同的特点，构建多波段的影像复原模型，利用波段的间高相似性及冗余信息进行相互约束和互补，最终得到高质量的高光谱影像。首先对高光谱影像进行波段聚类，将大量波段分为相关信息差别较大的少量类别；然后采用压缩感知理论对相同类别的一簇波段构建整体的变分训练多波段字典，通过该字典完成影像复原工作。本发明充分利用了多个波段间的相关性对目标影像进行恢复，并保持其光谱特性，复原结果具有更强的空间信息和光谱信息保持性能。

Description

一种结合波段聚类和稀疏表达的遥感影像复原方法

技术领域

本发明属于遥感图像处理技术领域，涉及一种基于波段聚类和稀疏表达的遥感影像复原方法，对于遥感影像中存在着较高相关性的波段，在此基础上利用高相关性波段组成的类中各波段之间的约束和互补关系，构建多通道的影像复原模型获得高质影像。

背景技术

随着信息处理技术的发展，人们对影像的要求也越来越高，但由于卫星、传感器或者大气等多方面的影响，最终获取的影像上可能会存在一些噪声或者光学模糊现象，所得到的影像往往分辨率很低而且会有退化变形和受噪声污染，导致遥感影像的高空间分辨率和高光谱分辨率无法同时获得，难以满足人们的视觉需要。在获取高质量的影像比较困难或者代价比较昂贵的时候，使用影像复原算法是很有必要的。影像复原主要可分为基于单幅影像复原和多幅影像复原的技术，其主要解决影像去噪，去模糊和超分重建问题。基于多幅影像的复原方法通常是指通过获取多幅不同时相或者不同位移的遥感影像，对多幅具有互补信息的低分辨率影像进行处理，重构出一幅高分辨率影像的技术，其方法有非均匀内插方法，基于多幅影像的Huber-Markov先验模型等。基于单幅影像的复原方法有基于邻域嵌套的超分辨率重建（LLE），基于样本学习的影像重建等。卫星遥感影像具有成像周期较长、地表情况复杂多变、大气辐射干扰严重等原因，使得同一场景的多时相遥感影像的获取较为困难，同时涉及到图像配准问题，而目前的配准算法其精度仍有待提高。可见发展遥感影像单幅影像复原的方法是重要的，但是现有针对多光谱遥感影像的复原算法大多将遥感影像视为普通图像进行复原，并没有考虑到遥感影像的光谱特性。因此，充分利用遥感影像丰富的光谱信息，研究针对单平台单时相遥感影像的复原具有较高的经济价值和现实意义。

发明内容

本发明的目的在于，针对现有技术的上述不足，提供一种结合波段聚类和稀疏表达的遥感影像复原方法，在利用波段间的高度相关性来进行复原的基础上，加入对噪声稳健的稀疏先验，对同一类具有较高相关性的波段同时训练字典并进行影像复原，这样构建的字典可以避免但波段训练字典时信息不足的情况，能够加入其他波段的信息构建一种多通道的稀疏模型，利用多波段直接的相关性，纠正单波段复原中字典训练不当导致复原不佳的情况。

本发明的技术方案为一种结合波段聚类和稀疏表达的遥感影像复原方法，包括以下步骤：

步骤1，将待处理影像根据相关系数指标，划分为具有相关性差异的多个波段类别；

步骤2，对每一波段类别的波段求取最佳字典，包括以下子步骤，

步骤2.1，设第m波段类X^m中包含L个波段，在每个波段上以像素点(x_i,y_j)为中心，设一个n×n的窗口，将L个同一光谱的n×n的窗口块提取出来组合为维度为nn×L的列向量nn=n×n；同时利用列向量组成第m波段类的样本库；

步骤2.2，稀疏系数向量求解，第一次迭代时利用步骤2,1所得样本库中部分列向量构建该m类的初始字典，基于初始字典采用OMP算法求解列向量的稀疏系数向量；后续迭代利用上一次执行步骤2.3更新的新字典采用OMP算法求解列向量的稀疏系数向量；

步骤2.3，字典更新，包括采用K-SVD算法更新字典，更新方式如下，

对应于字典中的某列原子，找到利用该列原子的影像列以及相应稀疏系数向量，通过对每列原子所得的残差做奇异值分解来更新该原子，使通过字典线性表达所得到的值与原始影像的值之差满足能量最小；

步骤2.4，判断是否达到迭代停止条件，是则确定第m波段类的当前字典为最佳字典，否则返回执行步骤2.2直到达到迭代停止条件；

步骤3，根据步骤2所得每一波段类别的最佳字典进行最小二乘求解得到高质影像。

而且，步骤1实现方式为，计算相关系数矩阵，以相关系数为指标，利用K-means算法对所有波段进行非监督聚类分为k类，每类中的波段相关性达到最大。

本发明的特色在于，不同波段之间的相关性存在差异，相邻波段的相关性一般较强，但是除了相邻波段以外，其它波段也可能与该波段有较强的相关性，因此可通过波段聚类的方法对高光谱影像波段分类，使高相关性波段聚为一类，改善单波段处理信息不足的问题；对于高光谱影像中处于同类波段构建多波段的稀疏先验，在稀疏先验本身对噪声稳健的特点下，在训练字典和用变分模型求解时，其它波段的信息加入在迭代过程中也起到重要作用，其相互之间形成约束和信息的互补，有效控制了复原精度。总之，本发明提出的方法可有效的用于遥感影像质量改善，获得更加准确的复原结果。

附图说明

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

具体实施方式

本发明技术方案可采用计算机软件技术实现自动运行流程。以下结合附图和实施例详细说明本发明技术方案。

如图1所示，本发明实施例包括以下步骤：

步骤1，波段选择及分类：将待处理影像根据相关系数指标，划分为具有相关性差异的多个波段类别。包括以下子步骤：

步骤1.1，对待处理高光谱影像计算相关系数矩阵，相关系数矩阵中两波段影像xⁱ,x^j（i,j＝1,...,B，i≠j，B表示该高光谱影像波段总数）之间的相关系数表示为其中分别表示影像xⁱ,x^j的平均值，N表示单波段的像元总数，像元标号a=1,...，N：

$r_{x^i{x}^j}=\frac{\underset{a=1}{\overset{N}{\mathrm{\Σ}}}(x_n^i-{\stackrel{\‾}{x}}^i)(x_n^j-{\stackrel{\‾}{x}}^j)}{\sqrt{\underset{a=1}{\overset{N}{\mathrm{\Σ}}}{(x_n^i-{\stackrel{\‾}{x}}^i)}^2}\sqrt{\underset{a=1}{\overset{N}{\mathrm{\Σ}}}{(x_n^j-{\stackrel{\‾}{x}}^j)}^2}}---\left(1\right)$

步骤1.2，以相关系数为准则，用非监督分类（如K-means）对波段聚类，可用如下优化公式表示：

$J=\underset{C}{\arg \min }\underset{b=1}{\overset{K}{\mathrm{\Σ}}}\underset{x^i\∈c_b}{\mathrm{\Σ}}R(x^i,{\mathrm{\μ}}_b)---\left(2\right)$

其中R(·,·)表示相关系数测度，得到的K个聚类记为c₁，c₂…c_K，μ_b是第c_b类的聚类中心，b＝1,2...，K，xⁱ表示第i个波段影像。如果两个波段影像xⁱ,x^j相关系数越高代表这两个波段之间越接近，选择相关系数大的波段影像分为一类，更能准确复原信息。

步骤2，对每一波段类别的波段求取最佳字典。本步骤对处于相同类别的波段构建多通道的正则化模型，建立基于多波段的稀疏先验影像复原模型。

步骤2.1，设第m波段类X^m中包含L个波段，在每个波段上以像素点(x_i,y_j)为中心，设一个n×n的窗口，将L个同一光谱的n×n的窗口块提取出来组合为维度为nn×L(nn=n×n)的列向量同时利用列向量组成第m类的样本库。

依据压缩感知（稀疏表达）理论，n×n图像块可表示为一组原子[φ₁,...,φ_h,...,φ_k]_h∈k∈Φ的线性组合x＝Φα，其对应字典（其维度满足条件：k>n）和稀疏向量α。由于高光谱影像波段数多，无法同时对整幅影像的多波段类构建过完备字典，因此可利用多通道全局先验（Cross-channelGlobalPrior），在高光谱影像中，对于第m波段类X^m，假设其包含L个波段，xⁱ(i＝1,...,L)表示第i波段影像，将该类影像中处于同一光谱的每个像元，在该像元所在波段中选取大小为n×n的窗口并将其表示为列向量，然后按照波段先后顺序依次排列该光谱中L个n×n的窗口，形成新的列向量其中m为类别号，i表示列向量号，N为单波段影像具有总像元数。同时利用列向量组成第m类的样本库。如公式3表示提取其中R_i为从第m波段类X^m提取列向量的矩阵，

$X_i^m=R_i{X}^m---\left(3\right)$

步骤2.2，稀疏系数向量求解，第一次迭代时利用步骤2,1所得样本库中部分列向量构建该m类的初始字典，基于初始字典求解列向量的稀疏系数向量；后续迭代利用上一次执行步骤2.3更新的新字典求解列向量的稀疏系数向量。本步骤用于估计稀疏系数向量：

影像复原是为了将低质影像Y恢复为高质量的影像X。其能够被表示为

Y＝DHX+v（4）

其中H和D是退化算子，v是噪声。当H和D为单位阵时，影像恢复过程在于去噪；当D为单位阵而H是一个模糊算子时，影像恢复过程在于去模糊；当D是降采样算子而H是模糊算子时，影像恢复就是超分辨率重建。

根据影像恢复基本公式，对应于第m类其恢复问题在加入全局先验的稀疏表达中可以表示为：

$\{{\widehat{\mathrm{\α}}}_i^m,{\hat{X}}^m\}=\underset{{\mathrm{\α}}_i^m,X^m}{\arg \min }{\mathrm{\λ}}^m{\left|\right|{\mathrm{DHX}}^m-Y^m\left|\right|}_2^2+\underset{i}{\mathrm{\Σ}}{\mathrm{\μ}}_i^m{\left|\right|{\mathrm{\α}}_i^m\left|\right|}_0+\underset{i}{\mathrm{\Σ}}{\left|\right|\mathrm{\Φ}{\mathrm{\α}}_i^m-R_i{X}^m\left|\right|}_2^2---\left(5\right)$

该公式的意思是求解满足该代价函数能量最小时的和X^m，其中表示nn×L（nn=n×n）的列向量的稀疏向量、Y^m表示输入的低质影像Y的第m波段类、X^m表示待恢复的高质量影像X的第m波段类,分别为和X^m最优估计值，和λ^m为第m波段类的正则化参数，通过调节和λ^m控制各正则项之间的比例关系达到影像最佳复原程度。公式（5）中第一项为数据一致项，第二项为多通道稀疏先验，第三项为多通道全局先验。这种多通道的先验模型能够有效保持光谱上的稀疏一致性。

求解稀疏系数向量α时，每次迭代过程中固定H、D以及Ф，初次迭代需要构建初始字典Ф。第m类中的的稀疏系数向量为满足，

${\widehat{\mathrm{\α}}}_i^m=\arg \min {\mathrm{\μ}}_i^m{\left|\right|{\mathrm{\α}}_i^m\left|\right|}_0+\left|\right|D{\mathrm{\α}}_i^m-X_i^m\left|\right|---\left(6\right)$

当更新稀疏系数向量时，可将H、D以及Ф看为整体。由于在每次迭代中固定H、D以及Ф，该问题可由正交追踪匹配算法（OMP）求解。

步骤2.3，字典更新，更新方式如下，

对应于字典中的某列原子，找到利用该列原子的影像列以及相应稀疏系数向量，通过对每列原子所得的残差做奇异值分解来更新该原子，使通过字典线性表达所得到的值与原始影像的值之差满足能量最小。

在该步骤中，固定H、D以及α更新字典。逐次更新字典Ф中每列φ_h(h＝1,2,...,k)。首先找出利用第h列原子的索引ω_h＝{d|α_d(h)≠0}，其中α_d(h)表示第h行值不为0的d个稀疏向量：

${\left|\right|R_d{X}^m-\underset{h=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\φ}}_h{\left({\mathrm{\α}}_d^m\right)}_h^T\left|\right|}_F^2={\left|\right|(R_d{X}^m-\underset{h\≠h_0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\φ}}_h{\left({\mathrm{\α}}_d^m\right)}_h^T)-{\mathrm{\φ}}_{h_0}{\left({\mathrm{\α}}_d^m\right)}_{h_0}^T\left|\right|}_F^2---\left(7\right)$

其中R_d表示从X^m中提取出d个列向量，表示索引中d个稀疏向量组的h行向量，φ_h表示第h列原子，h₀表示需要更新原子列序号，表示第m类稀疏向量中d个稀疏向量。

根据索引确定利用h₀号原子的中d个稀疏向量组，要优化和使得（7）式最小，可计算其稀疏表达残差E_l（representationerror）,

$E_l=R_d{X}^m-\underset{h\≠l}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\φ}}_h{\left({\mathrm{\α}}_d^m\right)}_h^T---\left(8\right)$

对E_l构建进行奇异值分解（SVD）得E_l＝UΔV^T，其中U和V为各列正交的特征向量，Δ为特征值。选择U的第一列更新字典中第h₀列原子，V与Δ(1,1)的乘积更新稀疏系数向量{α_d(l)}(d∈ω_l)。

步骤2.4，判断是否达到迭代停止条件，是则确定第m波段类的当前字典为最佳字典，否则返回执行步骤2.2直到达到迭代停止条件。实施例设置迭代停止条件为达到迭代次数，即重复2,3步骤迭代J次后，得到优化的字典Ф。J为预设的迭代次数。对下式（9）求解

${\hat{X}}^m=\underset{X^m}{\arg \min }\mathrm{\λ}{\left|\right|{\mathrm{DHX}}^m-Y^m\left|\right|}_2^2+\underset{i}{\mathrm{\Σ}}{\left|\right|\mathrm{\Φ}{\mathrm{\α}}_i^m-R_i{X}^m\left|\right|}_2^2---\left(9\right)$

可得第m类的影像复原结果，

${\hat{X}}^m={(\mathrm{\λ}{\left(\mathrm{DH}\right)}^T\left(\mathrm{DH}\right)+\underset{i}{\mathrm{\Σ}}{R}_i^T{R}_i)}^{-1}(\mathrm{\λ}{\left(\mathrm{DH}\right)}^T{Y}^m+\underset{i}{\mathrm{\Σ}}{R}_i^{T}D{\widehat{\mathrm{\α}}}_i^m)---\left(10\right)$

表示矩阵R_i的转置。

根据步骤2.1-2.4，依次处理各波段类，得到每一波段类别的最佳字典。

步骤3，根据步骤2所得每一波段类别的最佳字典进行最小二乘求解得到高质影像。

步骤2中构建多通道稀疏先验和多通道全局先验，建立总体影像复原模型，见公式(5)；将步骤2所得每一类的最佳字典带入最终的影像复原模型中进行最小二乘求解即可得到最终的高光谱影像复原结果。

具体实施时，可设计程序为对某一待恢复波段类进行迭代更新稀疏系数向量和字典，实现基于多波段的稀疏先验联合求解，得到该波段类复原结果。如果待处理波段类不为空，继续处理下一波段类，直到待处理波段类为空得到高光谱影像复原结果。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (2)

1.一种结合波段聚类和稀疏表达的遥感影像复原方法，其特征在于，包括以下步骤：

步骤1，将待处理影像根据相关系数指标，划分为具有相关性差异的多个波段类别；

步骤2，对每一波段类别的波段求取最佳字典，包括以下子步骤，

步骤2.1，设第m波段类X^m中包含L个波段，在每个波段上以像素点(x,y)为中心，设一个n×n的窗口，将L个同一光谱的n×n的窗口块提取出来组合为长度为nn×L的列向量nn＝n×n；同时利用列向量组成第m波段类的样本库，i＝1,...,N，N为单波段影像具有总像元数；

步骤2.2，稀疏系数向量求解，第一次迭代时利用步骤2.1所得样本库中部分列向量构建第m波段类的初始字典，基于初始字典采用OMP算法求解列向量的稀疏系数向量；所述OMP算法为正交追踪匹配算法；

步骤2.3，字典更新，包括采用K-SVD算法更新字典，更新方式如下，

对应于字典中的某列原子，找到利用该列原子的影像列以及相应稀疏系数向量，通过对每列原子所得的残差做奇异值分解来更新该原子；

步骤2.4，判断是否达到迭代停止条件，是则确定第m波段类的当前字典为最佳字典，否则返回执行步骤2.2直到达到迭代停止条件，后续迭代执行步骤2.2时利用上一次执行步骤2.3更新的新字典采用OMP算法求解列向量的稀疏系数向量；所述迭代停止条件为达到预设的迭代次数；

步骤3，根据步骤2所得每一波段类别的最佳字典进行最小二乘求解得到高质影像。

2.根据权利要求1所述结合波段聚类和稀疏表达的遥感影像复原方法，其特征在于：步骤1实现方式为，计算相关系数矩阵，以相关系数为指标，利用K-means算法对所有波段进行非监督聚类分为k类，选择相关系数大的波段影像分为一类。