[发明专利]基于双重字典学习的压缩感知图像超分辨率重建方法

摘要

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

主权项

1.基于双重字典学习的压缩感知图像超分辨率重建方法，其特征在于如下步骤：1、冗余字典、编码字典参数训练：为冗余字典，Ψ＝[ψ₁，ψ₂，…，ψ_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₁范数：fuction={Φ,Ψ,Θ}=argminΦ,Ψ,Θ{||S-ΦΘ||22+η||Θ-ΨS||22+λ|Θ|1}---(1)]]>(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)更新当前的Θ值；Θk+1=Tλ/2σΘ[(1-ησΘ)Θk+1σΘ(ΦT(S-ΦΘk)+ηΨS]---(2)]]>其中Θ^k+1、Θ^k分别代表迭代第k+1和第k步的Θ值，Φ^T代表Φ的转置。(3)定义操作π(d)＝d/max(1，||d||)，d为向量，该操作表示将向量投影到单位长度，定义σ_Φ＝2||ΘΘ^T||_F，使用公式(3)更新当前的Φ值：Φk+1=π(Φk+1σΦ(S-ΦkΘ)ΘT)---(3)]]>此处π()表示对Φ的每一列进行单位长度投影，其中Φ^k+1、Φ^k分别代表迭代第k+1和第k步的Φ值，Θ^T代表Θ的转置。(4)计算σ_Ψ＝2||SS^T||_F，使用公式(4)更新当前的Ψ值：Ψk+1=π(Ψk+1σΨ(Θ-ΨkS)ST)---(4)]]>此处π()表示对Ψ的每一行进行单位长度投影，其中Ψ^k+1、Ψ^k分别代表迭代第k+1和第k步的Ψ值，S^T代表S的转置。(5)迭代次数k＝k+1；(6)将当前计算的Θ、Ψ、Φ值和前一次计算的值分别代入公式(1)，计算目标函数的值，判断||function(k+1)-function(k)||22/||function(k)||22<ϵ]]>(function^(k+1)、function^k分别代表第k+1和第k步计算的function值)和k≥Max_Iter条件是否满足，满足其中任意一条，停止迭代，输出Φ值和Ψ值；否则重复(2)到(5)。2、自回归模型权值参数训练：第一步：读取图像库中超分辨率图像，将图像转为灰度图像，然后和一个低频高斯卷积核进行卷积(高斯卷积核大小为7×7，标准差为1.6)，得到低频图像，然后将原始图像和低频图像相减，差值反映图像的高频信息(这里将其命名为高频图像)，将高频图像分成大小为的图像小块，n＝49，找到与超分辨率图像中s_i对应位置处的图像小块，并依从左到右，从上到下，按列读取方式形成列向量，记为所有S＝[s₁，s₂，...s_M]对应的高频小块向量集合为Sh=[s1h,s2h,...sMh].]]>第二步：使用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}：μk=1mkΣi=1mksih]]>sih∈Ck---(5)]]>第三步：若s′_i表示s_i的中心像素值，q_i为s_i中s′_i的邻域像素值组成的向量，邻域大小取3×3(包括中心像素，q_i为去掉中心像素后组成的向量，为8个元素的列向量)，使用最小二次法计算式中的α_k，α_k为8×1的向量，对所有分类进行相同的处理，得到自回归模型权值参数组合{α₁，α₂，...α_K}；第四步：输出自回归模型权值参数组合{α₁，α₂，...α_K}和每一类的质心{μ₁，μ₂，...μ_K}；3、图像超分辨率重建冗余字典Φ、编码字典Ψ、自回归模型权值参数组合{α₁，α₂，...α_K}、每一类的质心{μ₁，μ₂，...μ_K}为预先训练得到的，一次训练可一直使用。第一步：读入需要进行重建的低分辨率图像Y，若为灰度图像，使用双立方插值将Y插值到需要的大小，表示为X⁽⁰⁾；若图像为RGB三色图像，则将图像变换为YCbCr色彩空间，将Y分量插值到需要的大小，表示为X⁽⁰⁾，设X⁽⁰⁾∈R^N″×1，则定义A、B为N″×N″维的系数矩阵；；第二步：计算系数矩阵A和B，其分为以下几步：(1)将X⁽⁰⁾分成大小为的图像小块，记为成x_i，i＝1，2，...N，N表示图像分块的个数(随输入图像的大小变化)，分块的时候相邻的块之间有重叠(横向或纵向重叠4个像素点宽度)；将X⁽⁰⁾和一个低频高斯卷积核进行卷积，得到低频图像，然后将X⁽⁰⁾和低频图像相减，差值反映X⁽⁰⁾的高频信息(这里将其命名为高频图像)，将高频图像分成大小为的图像小块，找到与x_i对应位置处的图像小块，并依从左到右，从上到下，按列读取方式形成列向量，记为(2)计算和所有{μ₁，μ₂，...μ_K}之间的欧氏距离，找到距离最小的那个，它的下标记为k_i，找到{α₁，α₂，...α_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)αi=argminα{||xi-Φα||2+γ4|α|1}---(10)]]>其中，γ₄为常数，取值范围为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)步，否则，返回α的值(即α_i＝α)；第六步：用公式(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)=(Σi=1NRiTRi)-1Σi=1NRiTΦαi---(12)]]>第八步：如果mod(k，P)＝＝0，且k≥Mid_Iter，将X^(k+1)取代第二步中的X⁽⁰⁾，重新计算A和B，并使用公式(13)计算τ_i，j：τi,j=cσn2σi,j+δ---(13)]]>c为一常数，σ_n是图像噪声的标准差，cσ_n取值范围为0.1～3.6，常规图像取0.1～0.6；δ为一比较小的常数，取0.35，σ_i，j计算如下：使用X^(k+1)，提取小块，并求与x_i相似的小块，对所有与x_i相似的小块向量计算然后提出l＝1，2，...L的第j个数，计算这些数的标准差即为σ_i，j。第九步：迭代次数k＝k+1；第十步：判断和k≥Max_Iter，其中有一个条件成立，则停止迭代，返回超分辨率图像X；否则重复第四步到第十步；第十一步：若输入图像为灰度图像，直接输出X，若为彩色图像，则将CbCr分量插值到和X相同的大小，然后将亮度X和色彩CbCr转化为RGB空间，输出重建后的图像。