CN103810755B - 基于结构聚类稀疏表示的压缩感知光谱图像重建方法 - Google Patents

Abstract

本发明公开了一种基于结构聚类稀疏表示的压缩感知光谱图像重建方法，解决了现有光谱图像重建方法没有充分利用空间和谱间相关性、难以精确恢复光谱图像局部结构的难题。其实现步骤是：1.对输入的光谱图像编码感知数据进行反投影，得到初始重建光谱图像；2.对重建光谱图像进行分块，得到一系列交叠的三维光谱图像块；3.利用基于结构聚类的稀疏表示方法对三维光谱图像块进行去噪；4.利用去噪后的光谱图像块恢复整个光谱图像；5.利用反投影技术对光谱图像进行更新；对步骤2‑5进行迭代，得到最终重建结果。实验结果表明本发明能重建出更加精细的光谱图像结构，重建的光谱图像具有更高的信噪比。

Description

基于结构聚类稀疏表示的压缩感知光谱图像重建方法

技术领域

本发明涉及图像处理技术领域，具体涉及一种压缩感知光谱图像的重建方法，主要用于对光谱图像的高质量恢复。

背景技术

光谱图像谱图合一的成像特点使其广泛应用于植被生物监测、土地动态监测、矿物分布调查以及目标侦查、近海环境监测等科研领域和军事领域。为了获得高分辨率光谱图像，众多学者分别从改善光学器件效能、研制高性能图像传感器等方面做了大量的研究工作，但是已取得的成果仍难以满足人们对高分辨率光谱图像的需求。其主要原因在于，传统的光谱成像系统的内在成像机理很难做到空间、时间、谱间同时高分辨。

为了克服传统光谱图像成像系统的内在缺陷，近年来国内外一些研究者，如美国杜克大学研究者于2007年在国际权威期刊Optics Express上发表了题为“Single-shotcompressive spectral imaging with a dual-disperser architecture”(Vol.15,No.21,pp.14013-14027,2007)的论文，提出了基于压缩感知理论的高分辨率计算光谱成像新方法。首先利用感知编码获得混叠光谱图像信息，然后通过计算重建的方法重建出原始光谱图像。相比传统的光谱成像技术，基于编码感知的光谱成像技术能获得更高的空间、时间和谱间分辨率。目前，基于压缩感知理论，计算光谱成像的感知编码技术相对成熟，而光谱图像的重建理论和方法研究相对滞后。现有多数光谱图像重建方法仍然采用传统的全变差、小波变换等灰度图像表示方法，无法有效利用光谱图像的空间、谱间相关性，难以重建出高质量的光谱图像。较差的光谱图像重建性能已成为制约计算光谱成像方法发展与应用的一个重要瓶颈问题。

光谱图像的观测过程可以表示为y＝Hx+n，其中x是高质量光谱图像，H表示观测矩阵，n代表观测过程中引入的噪声，y是得到的观测结果。光谱重建的目的就是根据观测结果y和观测矩阵H利用光谱重建方法获得高质量光谱图像x。光谱图像不同波段和局部区域包含大量的相关性，如何挖掘和利用光谱图像的谱间和空间相关性是重建高质量光谱图像的一个关键因素。

现有的光谱重建方法，如美国杜克大学的学者A.Wagadarikar等人于2008年在Applied Optics期刊上发表的“Single disperser design for coded aperturesnapshot spectral imaging”(Vol.47,No.10,pp.B44-B51,2008)一文中使用的重建方法仍然利用传统的小波域稀疏表示模型对各个谱段图像分别进行重建，这种重建方法只是在小波基下对每个谱段图像进行稀疏约束，没有充分利用到光谱图像的空间相关性和谱间相关性，难以重建出高质量光谱图像。另外，通过稀疏系数之间的结构相关性进行约束可以有效提升图像稀疏分解精度，这种结构稀疏模型已被证明是一个更加鲁棒，更加精确的稀疏模型。但是，现有的光谱重建技术主要是利用谱段图像的稀疏性进行约束，并没有利用到稀疏系数之间的相关性信息，因而难以精确重建出高质量的光谱图像。

发明内容

本发明的目的在于针对上述现有光谱重建技术的不足，提出一种基于结构聚类稀疏表示的压缩感知光谱图像重建方法，以提升光谱图像的重建质量。

实现本发明的技术思路是：构造基于三维字典的光谱图像联合结构稀疏表示模型，使得三维字典既包含了光谱图像的空间局部信息又包含了谱段之间的相关性信息，以实现光谱图像的空间和谱间信息的高精度建模；联合结构稀疏表示通过对稀疏系数之间的相关性进行约束来提升图像稀疏分解精度。具体步骤包括如下：

（1）输入光谱图像观测结果y∈R^h(w+n-1)m和观测矩阵H∈R^{h(w+n-1)m×hwn}，其中h,w分别表示每个谱段图像的高和宽，n表示光谱数，m表示观测次数，R代表实数域，x表示待求解的原始光谱图像，x∈R^hwn；

（2）设定反投影系数δ、稀疏系数阈值t₁、最大迭代步数P、更新步长L、光谱图像块类别数K、图像块相似度阈值τ和相似度权重参数h₀。初始恢复的光谱图像为x⁽⁰⁾=H^Ty，x⁽⁰⁾∈R^hwn，设当前迭代次数p＝0；

（3）对初始估计光谱图像x⁽⁰⁾进行分块，得到N个大小为s×s×n的三维光谱图像块，并对其进行向量化表示，得到图像块1≤i≤N。对每一个图像块x_i ⁽⁰⁾根据下列公式进行最近邻(K-NN)搜索，得到其相似块的下标集合和相似权重：

$G_i=\left\{j\right|{\left|\right|{x_i}^{\left(0\right)}-{x_j}^{\left(0\right)}\left|\right|}_2^2\≤\mathrm{\τ},1\≤j\≤N\},1\≤i\≤N,$

$w_{\mathrm{ij}}=\frac{1}{s_0}{e}^{-{\left|\right|{x_i}^{\left(0\right)}-{x_j}^{\left(0\right)}\left|\right|}_2^2/h_0},j\∈G_i,$

其中，G_i为图像块x_i ⁽⁰⁾的相似图像块的下标集，w_ij为刻画图像块x_i ⁽⁰⁾和x_j ⁽⁰⁾相似度的权重，s₀为归一化系数；

（4）设x_i ⁽⁰⁾对应的稀疏表示字典为，1≤k_i≤K，k_i表示第i个图像块x_i ⁽⁰⁾所属的类别号，并初始化字典为DCT变换基；

对以下步骤进行迭代重建光谱图像：

（5）对当前的光谱图像x^(p)按下式进行反投影重建，得到反投影结果x^(p+1/2)：

x^(p+1/2)＝x^(p)+δH^T(y-Hx^(p))；

（6）对反投影结果x^(p+1/2)进行分块，得到N个三维图像块1≤i≤N；

（7）利用相似图像块下标集G_i和权重w_ij按下式计算每个图像块x_i ^(p+1/2)所对应的k-NN聚类中心

${\mathrm{\μ}}_i=\underset{j\∈G_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{x_j}^{(p+1/2)},1\≤i\≤N;$

（8）利用三维字典和聚类中心μ_i对三维图像块按下式进行稀疏重建得到更新后的三维光谱图像块

$x_i^{(p+1)}=D_{k_i}\·\mathrm{soft}({D_{k_i}}^T(x_i^{(p+1/2)}-{\mathrm{\μ}}_i),t_1)+{\mathrm{\μ}}_i,1\≤i\≤N,$

其中，为软阈值函数，t₁为阈值；

（9）利用更新后的三维图像块1≤i≤N，按下式重建出优化后的光谱图像x^{(p +1)}：

$x^{(p+1)}={\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{R_i}^T{R}_i\right)}^{-1}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{R_i}^T{x}_i^{(p+1)}\right),$

其中，表示从光谱图像x^(p+1/2)得到光谱图像块对应的生成矩阵，即R_i ^T表示为R_i的转置；

（10）当mod(p,L)＝0时利用优化后的光谱图像x^(p+1)更新相似图像块下标集G_i和权重w_ij，其中mod(p,L)表示当前迭代步数p除以更新步长L后的余数；

（11）当mod(p,L)＝0时利用优化后的光谱图像x^(p+1)更新三维字典

（12）当mod(p,L)＝0时更新稀疏系数阈值t₁；

（13）p＝p+1，当p≤P时，返回步骤（5）；否则终止迭代，并输出最终光谱图像重建结果x＝x^(P)。

本发明与现有技术相比具有以下优点：

第一，本发明是基于三维图像块的光谱图像重建方法，该三维图像块既包含了光谱图像的空间维度信息又包含了光谱维度信息，能够充分利用到光谱图像的空间稀疏特性和谱间相关性；另外三维自适应字典学习方法的应用克服了传统字典难以表征图像中出现的各种局部结构的缺点，能够更好的刻画光谱图像的结构特征；

第二，本发明由于在图像重建中利用到了稀疏系数的结构相关性，能更加精确恢复光谱图像；

第三，本发明由于在参数设置中对阈值参数进行自适应的选取，使重构算法具有更好的鲁棒性。

附图说明

图1为本发明的总流程图；

图2为现有aviris光谱图像的第6个谱段图像；

图3(a)为用现有CASSI方法在观测8次时得到的重构光谱图像的第6个谱段图像，图3(b)为用本发明方法在观测8次时得到的重构光谱图像的第6个谱段图像；

图4(a)为用现有CASSI方法在观测15次时得到的重构光谱图像的第6个谱段图像，图4(b)为用本发明方法在观测15次时得到的重构光谱图像的第6个谱段图像；

图5(a)为用现有CASSI方法在观测20次时得到的重构光谱图像的第6个谱段图像，图5(b)为用本发明方法在观测20次时得到的重构光谱图像的第6个谱段图像。

具体实施方式

参照图1，本发明是基于结构聚类稀疏表示的压缩感知光谱图像重建方法，其实现步骤如下：

步骤1：输入观测数据。

输入光谱图像观测结果y∈R^h(w+n-1)m和观测矩阵H∈R^{h(w+n-1)m×hwn}，其中h,w分别表示每个谱段图像的高和宽，n表示光谱数，m表示观测次数，R代表实数域，x表示待求解的原始光谱图像，x∈R^hwn。

步骤2：设定反投影系数δ、稀疏系数阈值t₁、最大迭代步数P、更新步长L、光谱图像块类别数K、图像块相似度阈值τ和相似度权重参数h₀。初始恢复的光谱图像为x⁽⁰⁾=H^Ty，x⁽⁰⁾∈R^hwn，设当前迭代次数p＝0。

步骤3：对初始估计光谱图像x⁽⁰⁾进行分块，得到N个大小为s×s×n的三维光谱图像块，并对其进行向量化表示，得到图像块1≤i≤N。对每一个图像块x_i ⁽⁰⁾根据下列公式进行最近邻(K-NN)搜索，得到其相似块的下标集合和相似权重：

$G_i=\left\{j\right|{\left|\right|{x_i}^{\left(0\right)}-{x_j}^{\left(0\right)}\left|\right|}_2^2\≤\mathrm{\τ},1\≤j\≤N\},1\≤i\≤N,$

$w_{\mathrm{ij}}=\frac{1}{s_0}{e}^{-{\left|\right|{x_i}^{\left(0\right)}-{x_j}^{\left(0\right)}\left|\right|}_2^2/h_0},j\∈G_i,$

其中，G_i为图像块x_i ⁽⁰⁾的相似图像块的下标集，w_ij为刻画图像块x_i ⁽⁰⁾和x_j ⁽⁰⁾相似度的权重，s₀为归一化系数。

步骤4：设x_i ⁽⁰⁾对应的稀疏表示字典为1≤k_i≤K，k_i表示第i个图像块x_i ⁽⁰⁾所属的类别号，并初始化字典为DCT变换基。

对以下步骤进行迭代重建光谱图像：

步骤5：对当前的光谱图像x^(p)按下式进行反投影重建，得到反投影结果x^(p+1/2)：

x^(p+1/2)＝x^(p)+δH^T(y-Hx^(p))。

步骤6：对反投影结果x^(p+1/2)进行分块，得到N个三维图像块1≤i≤N。

步骤7：利用相似图像块下标集G_i和权重w_ij按下式计算每个图像块所对应的k-NN聚类中心

${\mathrm{\μ}}_i=\underset{j\∈G_i}{\mathrm{\Σ}}{w}_{\mathrm{ij}}{x_j}^{(p+1/2)},1\≤i\≤N.$

步骤8：利用三维字典和聚类中心μ_i对三维图像块按下式进行稀疏重建得到更新后的三维光谱图像块

其中，为软阈值函数，t₁为阈值。

步骤9：利用更新后的三维图像块1≤i≤N，按下式重建出优化后的光谱图像x^(p+1)：

$x^{(p+1)}={\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{R_i}^T{R}_i\right)}^{-1}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{R_i}^T{x}_i^{(p+1)}\right),$

其中，表示从光谱图像x^(p+1/2)得到光谱图像块对应的生成矩阵，即R_i ^T表示为R_i的转置。

步骤10：当mod(p,L)＝0时利用优化后的光谱图像x^(p+1)更新相似图像块下标集G_i和权重w_ij，其中mod(p,L)表示当前迭代步数p除以更新步长L后的余数，具体更新步骤如下：

（10a）对光谱图像x^(p+1)进行分块，得到N个大小为s×s×n的三维光谱图像块，并对其进行向量化表示，得到图像块1≤i≤N；

（10b）对每个图像块x_i ^(p+1)按如下公式更新下标集G_i和权重w_ij：

$G_i=\left\{j\right|{\left|\right|{x_i}^{(p+1)}-{x_j}^{(p+1)}\left|\right|}_2^2\≤\mathrm{\τ},1\≤j\≤N\},1\≤i\≤N,$

$w_{\mathrm{ij}}=\frac{1}{s_0}{e}^{-{\left|\right|{x_i}^{(p+1)}-{x_j}^{(p+1)}\left|\right|}_2^2/h_0},j\∈G_i,$

其中，s₀为归一化系数。

步骤11：当mod(p,L)＝0时利用优化后的光谱图像x^(p+1)更新三维字典具体更新步骤如下：

（11a）对光谱图像x^(p+1)进行分块，得到N个大小为s×s×n的三维光谱图像块，并对其进行向量化表示，得到1≤i≤N；

（11b）对N个三维光谱图像块利用K均值方法划分成K类，为每个三维光谱图像块x_i ^(p+1)找到其所属类的类别号k_i，1≤k_i≤K，对每一类三维图像块利用主成分分析法学习一个字典，得到三维字典集D＝{D₁,D₂,...,D_k,...,D_K}，类别号k_i所对应的字典即为更新后的三维字典，对于三维光谱图像块x_i ^(p+1)，将采用三维字典进行稀疏编码。

步骤12：当mod(p,L)＝0时根据贝叶斯公式对稀疏系数阈值t₁进行自适应更新：

$t_1=\frac{2\sqrt{2}{{\mathrm{\σ}}_n}^2}{\mathrm{\σ}+\mathrm{\ϵ}},$

其中，σ_n为噪声标准差，ε的值为2.2×10^-16，L_i＝|G_i|，|G_i|表示下标集G_i中元素的个数，σ为稀疏系数对应的标准差，α_j为稀疏系数。

步骤13：p＝p+1，当p≤P时，返回步骤（5）；否则终止迭代，并输出最终光谱图像重建结果x＝x^(P)。

本发明的效果可以通过如下仿真实验进一步说明：

1.仿真条件：

本发明的仿真运行系统为Intel(R)Core(TM)i5-2400 CPU @3.10GHz，64位Windows操作系统，仿真软件采用MATLAB(R2013a)，实验数据采用现有90谱段aviris光谱图像数据。

观测方法采用美国杜克大学的学者M.E.Gehm，R.Johm，D.J.Brady，R.M.Willet和T.J.Schualz在论文“Single-shot compressive spectral imaging with a dual-disperser architecture”OPTICS EXPRESS,Vol.15,No.21,pp.14013-14027,2007中提出的对光谱图像各个谱段先编码后错位叠加的观测方法。

仿真实验中采用峰值信噪比PSNR指标来评价光谱重建效果，其定义为：

$\mathrm{PSNR}=10\·{\log }_{10}\left(\frac{{\mathrm{MAX}}^2}{\mathrm{MSE}}\right),$

$\mathrm{MSE}=\frac{{\left|\right|x-\hat{x}\left|\right|}_F^2}{h\×w\×n},$

其中MAX为光谱图像像素值的峰值，本实验中峰值MAX取1，MSE为原光谱图像x和重构图像之间的均方误差，x,峰值信噪比PSNR越大表示重建结果越好。

2.仿真内容：

实验一，对现有aviris光谱图像在观测8次时分别利用现有CASSI方法和本发明方法进行重建，结果如图3，其中图3(a)为现有CASSI方法得到的重建结果中的第6个谱段图像，图3(b)为本发明方法得到的重建结果中的第6个谱段图像。

实验二，对现有aviris光谱图像在观测15次时分别利用现有CASSI方法和本发明方法进行重建，结果如图4，其中图4(a)为现有CASSI方法得到的重建结果中的第6个谱段图像，图4(b)为本发明方法得到的重建结果中的第6个谱段图像。

实验三，对现有aviris光谱图像在观测20次时分别利用现有CASSI方法和本发明方法进行重建，结果如图5，其中图5(a)为现有CASSI方法得到的重建结果中的第6个谱段图像，图5(b)为本发明方法得到的重建结果中的第6个谱段图像。

从图3-图5可以看出，本发明基于结构聚类稀疏表示的压缩感知重构方法恢复出的结果具有更好的边缘和纹理结构，能够更佳的逼近原始图像，视觉效果更好。

表1给出了8次观测、10次观测、15次观测、20次观测和30次观测五种情况下用现有CASSI方法和本发明方法对现有aviris光谱图像观测结果进行重建得到的重建光谱图像的峰值信噪比(PSNR)。

表1aviris光谱图像重建结果峰值信噪比（单位：dB）

观测次数

8次

10次

15次

20次

30次

CASSI方法	42.53	43.00	44.46	45.78	48.88
						本发明方法	43.13	45.16	46.03	47.94	50.45

从表1中我们可以看出，在这5种观测情况下本发明方法对aviris光谱图像的重建结果的PSNR值都要高于CASSI方法，即本发明方法具有更好的光谱图像重建效果。

Claims (1)

1.一种基于结构聚类稀疏表示的压缩感知光谱图像重建方法，包括如下步骤：

(1)输入光谱图像观测结果y∈R^h(w+n-1)m和观测矩阵H∈R^{h(w+n-1)m×hwn}，其中h,w分别表示每个谱段图像的高和宽，n表示光谱数，m表示观测次数，R代表实数域，x表示待求解的原始光谱图像，x∈R^hwn；

(2)设定反投影系数δ、稀疏系数阈值t₁、最大迭代步数P、更新步长L、光谱图像块类别数K、图像块相似度阈值τ和相似度权重参数h₀，初始恢复的光谱图像为x⁽⁰⁾＝H^Ty，x⁽⁰⁾∈R^hwn，设当前迭代次数p＝0；

(3)对初始估计光谱图像x⁽⁰⁾进行分块，得到N个大小为s×s×n的三维光谱图像块，并对其进行向量化表示，得到图像块1≤i≤N，对每一个图像块x_i ⁽⁰⁾根据下列公式进行最近邻(K‐NN)搜索，得到其相似块的下标集和相似权重：

$G_i=\left\{j\right|\left|\right|{x_i}^{\left(0\right)}-{x_j}^{\left(0\right)}|{|}_2^2\≤\mathrm{\τ},1\≤j\≤N\},1\≤i\≤N,$

$w_{ij}=\frac{1}{s_0}{e}^{-\left|\right|{x_i}^{\left(0\right)}-{x_j}^{\left(0\right)}|{|}_2^2/h_0},j\∈G_i,$

其中，G_i为图像块x_i ⁽⁰⁾的相似图像块的下标集，w_ij为刻画图像块x_i ⁽⁰⁾和x_j ⁽⁰⁾相似度的权重，s₀为归一化系数；

(4)设x_i ⁽⁰⁾对应的稀疏表示字典为1≤k_i≤K，k_i表示第i个图像块x_i ⁽⁰⁾所属的类别号，并初始化字典为DCT变换基；

对以下步骤进行迭代重建光谱图像：

(5)对当前的光谱图像x^(p)按下式进行反投影重建，得到反投影结果x^(p+1/2)：

x^(p+1/2)＝x^(p)+δH^T(y-Hx^(p))；

(6)对反投影结果x^(p+1/2)进行分块，得到N个三维图像块1≤i≤N；

(7)利用相似图像块下标集G_i和权重w_ij按下式计算每个图像块所对应的k‐NN聚类中心

${\mathrm{\μ}}_i=\underset{j\∈G_i}{\Σ}{w}_{ij}{x_j}^{(p+1/2)},1\≤i\≤N;$

(8)利用三维字典和聚类中心μ_i对三维图像块按下式进行稀疏重建得到更新后的三维光谱图像块

$x_i^{(p+1)}=D_{k_i}\·soft\left({D_{k_i}}^T\right(x_i^{(p+1/2)}-{\mathrm{\μ}}_i),t_1)+{\mathrm{\μ}}_i,1\≤i\≤N,$

其中，为软阈值函数，t₁为阈值；

(9)利用更新后的三维图像块1≤i≤N，按下式重建出优化后的光谱图像x^(p+1)：

$x^{(p+1)}={\left(\underset{i=1}{\overset{N}{\Σ}}{R_i}^T{R}_i\right)}^{-1}\left(\underset{i=1}{\overset{N}{\Σ}}{R_i}^T{x}_i^{(p+1)}\right),$

其中，表示从光谱图像x^(p+1/2)得到光谱图像块对应的生成矩阵，即R_i ^T表示为R_i的转置；

(10)当mod(p,L)＝0时利用优化后的光谱图像x^(p+1)更新相似图像块下标集G_i和权重w_ij，其中mod(p,L)表示当前迭代步数p除以更新步长L后的余数：

(10a)对光谱图像x^(p+1)进行分块，得到N个大小为s×s×n的三维光谱图像块，并对其进行向量化表示，得到图像块1≤i≤N；

(10b)对每个图像块x_i ^(p+1)按如下公式更新下标集G_i和权重w_ij：

$G_i=\left\{j\right|\left|\right|{x_i}^{(p+1)}-{x_j}^{(p+1)}|{|}_2^2\≤\mathrm{\τ},1\≤j\≤N\},1\≤i\≤N,$

$w_{ij}=\frac{1}{s_0}{e}^{-\left|\right|{x_i}^{(p+1)}-{x_j}^{(p+1)}|{|}_2^2/h_0},j\∈G_i,$

其中，s₀为归一化系数；

(11)当mod(p,L)＝0时利用优化后的光谱图像x^(p+1)更新三维字典

(11a)对光谱图像x^(p+1)进行分块，得到N个大小为s×s×n的三维光谱图像块，并对其进行向量化表示，得到图像块1≤i≤N；

(11b)对N个三维光谱图像块利用K均值方法划分成K类，为每个三维光谱图像块x_i ^(p+1)找到其所属类的类别号k_i，1≤k_i≤K，对每一类三维图像块利用主成分分析法学习一个字典，得到三维字典集D＝{D₁,D₂,...,D_k,...,D_K}，类别号k_i所对应的字典即为更新后的三维字典对于三维光谱图像块x_i ^(p+1)，将采用三维字典进行稀疏编码；

(12)当mod(p,L)＝0时利用优化后的光谱图像x^(p+1)更新稀疏系数阈值t₁：

$t_1=\frac{2\sqrt{2}{{\mathrm{\σ}}_n}^2}{\mathrm{\σ}+\mathrm{\ϵ}},$

其中，σ_n为噪声标准差，ε的值为2.2×10^-16，L_i＝|G_i|，|G_i|表示下标集G_i中元素的个数，σ为稀疏系数对应的标准差，α_j为稀疏系数；

(13)p＝p+1，当p≤P时，返回步骤(5)；否则终止迭代，并输出最终光谱图像重建结果x＝x^(P)。