基于双重字典学的压缩感知图像超分辨率重建方法

摘要

本发明涉及基于双重字典学的压缩感知图像超分辨率重建方法，其包括如下步骤：冗余字典、编码字典参数训练，自回归模型权值参数训练，用训练好的冗余字典、编码字典、自回归模型权值参数对单帧低分辨率图像进行超分辨率重建，本算法具有重建效果好的特点，适用于医学成像、卫星遥感遥测、军事侦察与定位以及城市安防等诸多领域。

主权项

基于双重字典学习的压缩感知图像超分辨率重建方法，其特征在于如下步骤：a、冗余字典、编码字典参数训练：为冗余字典，Ψ＝[ψ₁，ψ₂，…，ψ_n]∈R^m×n为编码字典，m、n为正整数，其中m＝512，n＝49，双重字典指的是同时产生冗余字典和编码字典两个字典；第一步：读取图像库中超分辨率图像，将超分辨率图像转为灰度图像，然后分成大小为的小块样本，将得到的图像小块按照从左到右，从上到下，按列读取方式形成列向量，用s_i∈Rⁿ，i＝1，2，...，Q表示每个小块形成的列向量，Q为总的列向量的个数；第二步：计算s_i的方差Var(s_i)，只保留Var(s_i)大于阈值TH的向量，其中TH取值范围为：4.5～20，最终得到训练样本集合S＝[s₁，s₂，...s_M]，M大于120000；第三步：对公式(1)进行求解，采用迭代方法求解冗余字典Φ，编码字典Ψ，Θ为稀疏系数，λ、η为常数，η取等于1的数值，λ取值范围为0.05～0.2，表示求l₂范数，||₁表示求l₁范数：$\mathrm{fuction}=\{\Phi ,\Psi ,\Theta \}=\arg \underset{\Phi ,\Psi ,\Theta }{\min }\{{\left|\right|S-\mathrm{\Phi \Theta }\left|\right|}_2^2+\eta {\left|\right|\Theta -\mathrm{\Psi S}\left|\right|}_2^2+\lambda {\left|\Theta \right|}_1\}---\left(1\right)$(1)用高斯随机矩阵初始化冗余字典Φ，用单位矩阵初始化编码字典Ψ，用全零矩阵初始化稀疏系数Θ，迭代次数k＝0，最大迭代次数Max_Iter取800～1500，迭代收敛控制因子ε＝10^‑6；(2)定义(T_ζ[O])_i，j＝sign(O_i，j)max{|O_i，j|‑ζ，0}为阈值操作算子，ζ代表阈值操作变量，O代表阈值操作矩阵变量，O_i，j表示矩阵O中下标为(i，j)的元素，sign()为符号操作算子，取σ_Θ＝2||Φ^TΦ+ηI||_F，||||_F表示求Frobenius范数，I为单位矩阵，使用公式(2)更新当前的Θ值；${\Theta }^{k+1}=T_{\lambda /2{\sigma }_{\Theta }}[(1-\frac{\eta }{{\sigma }_{\Theta }}){\Theta }^k+\frac{1}{{\sigma }_{\Theta }}({\Phi }^T(S-\Phi {\Theta }^k)+\mathrm{\eta \Psi S})]---\left(2\right)$其中Θ^k+1、Θ^k分别代表迭代第k+1和第k步的Θ值，Φ^T代表Φ的转置；(3)定义操作π(d)＝d/max(1，||d||)，d为向量，该操作表示将向量投影到单位长度，定义σ_Φ＝2||ΘΘ^T||_F，使用公式(3)更新当前的Φ值：${\Phi }^{k+1}=\pi ({\Phi }^k+\frac{1}{{\sigma }_{\Phi }}(S-{\Phi }^k\Theta ){\Theta }^T)---\left(3\right)$此处π()表示对Φ的每一列进行单位长度投影，其中Φ^k+1、Φ^k分别代表迭代第k+1和第k步的Φ值，Θ^T代表Θ的转置；(4)计算σ_Ψ＝2||SS^T||_F，使用公式(4)更新当前的Ψ值：${\Psi }^{k+1}=\pi ({\Psi }^k+\frac{1}{{\sigma }_{\Psi }}(\Theta -{\Psi }^{k}S)S^T)---\left(4\right)$此处π()表示对Ψ的每一行进行单位长度投影，其中Ψ^k+1、Ψ^k分别代表迭代第k+1和第k步的Ψ值，S^T代表S的转置；(5)迭代次数k＝k+1；(6)将当前计算的Θ、Ψ、Φ值和前一次计算的值分别代入公式(1)，计算目标函数的值，判断和k≥Max_Iter条件是否满足，满足其中任意一条，停止迭代，输出Φ值和Ψ值，否则重复(2)到(5)，其中function^(k+1)、function^k分别代表第k+1和第k步计算的function值；b、自回归模型权值参数训练：第一步：读取图像库中超分辨率图像，将图像转为灰度图像，然后和一个大小为7×7，标准差为1.6的低频高斯卷积核进行卷积，得到低频图像，然后将原始图像和低频图像相减，差值反映图像的高频信息，我们将其命名为高频图像，将高频图像分成大小为的图像小块，n＝49，找到与超分辨率图像中s_i对应位置处的图像小块，并按照从左到右，从上到下，按列读取方式形成列向量，记为所有S＝[s₁，s₂，...s_M]对应的高频小块向量集合为第二步：使用K‑means分类算法将S^h分成K类{C₁，C₂，...C_K}，K取200，m_k表示C_k中向量的个数，使用公式(5)计算每一类的质心μ_k，k＝1，2，...，K，根据S^h的分类结果，将S＝[s₁，s₂，...s_M]也分成K类，表示为{S₁，S₂，...S_K}：${\mu }_k=\frac{1}{m_k}\underset{i=1}{\overset{m_k}{\Sigma }}{s}_i^h,s_i^h\in C_k---\left(5\right)$第三步：若s′_i表示s_i的中心像素值，q_i为s_i中s′_i的邻域像素值组成的向量，邻域大小取3×3，q_i为去掉中心像素后组成的向量，为8个元素的列向量，使用最小二乘法计算式中的a_k，a_k为8×1的向量，对所有分类进行相同的处理，得到自回归模型权值参数组合{a₁，a₂，...a_K}；第四步：输出自回归模型权值参数组合{a₁，a₂，...a_K}和每一类的质心{μ₁，μ₂，…μ_K}；c、图像超分辨率重建冗余字典Φ、编码字典Ψ、自回归模型权值参数组合{a₁，a₂，...a_K}、每一类的质心{μ₁，μ₂，...μ_K}为预先训练得到的，一次训练可一直使用；第一步：读入需要进行重建的低分辨率图像Y，若为灰度图像，使用双立方插值将Y插值到需要的大小，表示为X⁽⁰⁾；若图像为RGB三色图像，则将图像变换为YCbCr色彩空间，将Y分量插值到需要的大小，表示为X⁽⁰⁾，设则定义A、B为N″×N″维的系数矩阵；第二步：计算系数矩阵A和B，其分为以下几步：(1)将X⁽⁰⁾分成大小为的图像小块，记成为x_i，i＝1，2，...N，N表示图像分块的个数，N的值随输入图像的大小变化，分块的时候相邻的块之间有重叠，横向或纵向重叠4个像素点宽度；将X⁽⁰⁾和一个低频高斯卷积核进行卷积，得到低频图像，然后将X⁽⁰⁾和低频图像相减，差值反映X⁽⁰⁾的高频信息，在这里将其命名为高频图像，将高频图像分成大小为的图像小块，找到与x_i对应位置处的图像小块，并按照从左到右，从上到下，按列读取方式形成列向量，记为(2)计算和所有{μ₁，μ₂，...μ_K}之间的欧氏距离，找到距离最小的那个，它的下标记为k_i，找到{a₁，a₂，...a_K}中下标为k_i的权值参数，作为x_i的自回归模型权值参数(3)若x′_i表示x_i的中心像素值，χ_i为x_i中x′_i的邻域像素值组成的向量，使用公式(6)计算矩阵A；i、j是坐标变量，取正整数，取值范围为1～N″；(4)在X⁽⁰⁾图像中分成的小块中，寻找的L个相似小块，变量l＝1，2，...L，L表示相似小块的数量，为正整数，L取值范围为7～10，代表X⁽⁰⁾中和x_i相似的小块，计算中的值，其中表示归一化因子，h是常数，取值范围为65～70，设为权值向量，表示所有相似小块的中心像素集合，则使用公式(7)计算矩阵B：i、l是坐标变量，取正整数，取值范围为1～N″；第三步：预设定：γ₁取值范围0.008～0.01、γ₂取值范围0.04～0.1、γ₃取值6.5、P＝20、e＝10^‑6、Mid_Iter＝100以及最大的迭代次数Max_Iter＝150，设定常数矩阵τ＝0，初始化迭代次数k＝0；第四步：设I表示单位矩阵，I矩阵大小和矩阵A、B一样，D表示下采样矩阵，D根据重建倍数设定，H为高斯模糊矩阵，它是高斯卷积核，在重建因子为3时，高斯卷积核大小为7×7，标准差为1.6，重建因子为2时，高斯卷积核大小为5×5，标准差为0.9～1.1，在重建因子是4时，高斯卷积核大小为7×7，标准差为1.7～1.8，根据高斯卷积核设定，重建矩阵是一循环矩阵，计算公式(8)：X^(k+1/2)＝X^(k)+γ₃[(DH)^TY‑(DH)^TDHX^(k)]‑γ₁(I‑A)^T(I‑A)X^(k)‑γ₂(I‑B)^T(I‑B)X^(k)      (8)X^(k+1/2)、X^(k)分别代表迭代第k+1/2和第k步的重建结果；第五步：若R_i表示将x_i从X中裁剪出来，即：x_i＝R_iX，如果迭代次数k小于Mid_Iter，使用公式(9)计算稀疏系数Θ^(k+1/2)，Θ^(k+1/2)＝[α₁，α₂，...α_N]；否则，使用公式(10)计算α_i：Θ^(k+1/2)＝[ΨR₁X^(k+1/2)，ΨR₂X^(k+1/2)，…ΨR_NX^(k+1/2)]       (9)${\alpha }_i=\underset{\alpha }{\arg \min }\{{\left|\right|x_i-\mathrm{\Phi \alpha }\left|\right|}^2+{\gamma }_4{\left|\alpha \right|}_1\}---\left(10\right)$其中，γ₄为常数，取值范围为0.1～0.2，公式(10)采用特征符号寻找算法求解，具体过程如下：(a)定义向量θ∈R^m×1，θ_j代表向量θ中的第j个元素，θ_j∈{‑1，0，1}，初始化定义活动集合β＝{}，并将其初始化为空集；(b)对于α中为0的元素，计算找出j值，α_j代表α的第j各元素，如果则θ_j＝‑1，β＝β∪{j}，如果：则θ_j＝1，β＝β∪{j}；(c)选出Φ中下标为β的列向量组成选出α和θ中下标为β的元素分别组成和计算然后逐个检查和相应位置元素，看哪些元素改变了符号，符号改变指由正变负或由负变正，设有Num个位置的值改变符号，将中变符号位置的元素置0，每次只对一个位置的值操作，其它位置的值保持不变，用表示值改变了的新向量，则有Num种取值，将的Num种取值分别代入中，找出使值最小的那个的取值，把它赋值给α中下标为β的元素与中相应位置元素做相同的改变，将和α中元素变为0的下标从β中移除，更新θ＝sign(α)；(d)判断α中不为0的元素是否满足：若不满足，则执行(c)步，否则判断α中为0的元素是否满足：若不满足则执行(b)步，否则，返回α的值；第六步：用公式(11)计算Θ^(k+1)，定义操作(T_τ[Z])_i，j＝sign(Z_i，j)max{|Z_i，j|‑τ_i，j，0}，Z是该操作的对象：Θ^(k+1)＝T_τ[Θ^(k+1/2)]              (11)第七步：使用公式(12)计算X^(k+1)：$X^{(k+1)}={\left(\underset{i=1}{\overset{N}{\Sigma }}{R}_i^T{R}_i\right)}^{-1}\underset{i=1}{\overset{N}{\Sigma }}{R}_i^T\Phi {\alpha }_i---\left(12\right)$第八步：如果mod(k，P)＝＝0，且k≥Mid_Iter，将X^(k+1)取代第二步中的X⁽⁰⁾，重新计算A和B，并使用公式(13)计算τ_i，j：${\tau }_{i,j}=\frac{c{\sigma }_n^2}{{\sigma }_{i,j}+\delta }---\left(13\right)$c为一常数，σ_n是图像噪声的标准差，cσ_n取值范围为0.1～3.6，常规图像取0.1～0.6；δ为一比较小的常数，取0.35，σ_i，j计算如下：使用X^(k+1)，提取小块，并求与x_i相似的小块，对所有与x_i相似的小块向量计算然后提出的第j个数，计算这些数的标准差即为σ_i，j；第九步：迭代次数k＝k+1；第十步：判断和k≥Max_Iter，其中有一个条件成立，则停止迭代，返回超分辨率图像X；否则重复第四步到第十步；第十一步：若输入图像为灰度图像，直接输出X，若为彩色图像，则将CbCr分量插值到和X相同的大小，然后将亮度X和色彩CbCr转化为RGB空间，输出重建后的图像。