CN103400356B - 一种基于通用隐马尔可夫树模型的加权图像压缩感知方法 - Google Patents

Abstract

本发明提出了一种基于通用隐马尔可夫树模型的加权图像压缩感知方法，其将基于uHMT模型的迭代加权L1最小化算法应用于图像的压缩感知重建，根据自然图像的自身特性，采用uHMT模型代替HMT模型，直接给出模型参数，省去计算量较大的训练步骤，有效降低图像重建的计算复杂度。同时，由于权值的计算方法对重建精度有很大影响，通过本发明给出的改进的权值计算方法既将刻画小波系数的uHMT模型引入了图像重建，又降低了出现排斥模型现象的可能性，从而提高了图像的重建精度。

Description

一种基于通用隐马尔可夫树模型的加权图像压缩感知方法

技术领域

本发明涉及图像压缩领域，尤其涉及一种基于通用隐马尔可夫树模型的加权图像压缩感知方法。

背景技术

随着数字化和信息时代的来临，图像的数据量在不断扩大，给图像的传输处理带来了巨大压力，因此需要对图像进行有效的压缩和重构。近年来，Donoho等人提出了一种新的理论——压缩感知理论（CS）。该理论突破了传统奈奎斯特采样中关于采样速率的限制，采样速率不再决定于信号的带宽，而决定于信号中的结构和内容。其核心思想是将传统压缩过程中的采样和压缩过程合并进行，对信号进行随机投影得到少量观测值，然后通过重建算法从观测值中重建出原始信号。

目前CS重建算法包括贪婪算法和L1范数最小化算法两大类。标准的CS图像重建只利用图像稀疏性对其进行重建，并没有使用图像本身特殊的结构信息。但其实在图像的小波变换域中，存在很多先验结构，如图像小波系数幅值在尺度间呈指数衰减，与边缘纹理相对应的大系数呈树状分布等，这些先验特征引入图像CS重建中可提高图像重建质量。图像小波域的隐马尔可夫树（HMT）是刻画小波系数的分布特征的一种简洁而精确的树模型，自然图像的小波系数可由一种四叉树结构表示，如图1所示。对一个大小为M×N的图像，每一个位于尺度j＜log₂(N)的系数具有位于尺度j+1的4个子系数，每一个位于尺度j＞1的系数都有位于尺度j-1的唯一的父系数。HMT模型已广泛应用到图像复原、去噪中。

MarcoFD等2008年将小波HMT模型引入到分段光滑信号的压缩感知重建中，将HMT模型与迭代加权L1最小化（IRWL1）算法相结合，采用基于小波HMT模型的加权方式来重置权值，直接将HMT模型参数引入权值的修正中，再利用迭代加权L1最小化算法重建信号。

HMT模型及其参数估计

小波系数具有以下的两个统计特性：非高斯性和传递性。

在HMT模型中，非高斯性可以通过使用具有一个隐二值状态的两个零均值高斯分布的混合分布很好地刻画。对每个小波系数x_n（n∈[1,N],N为小波系数的长度）用变量S_n∈{1,2}表示不可观测的隐状态。S_n=1的高斯成分具有较小的方差，刻画与图像光滑部分对应的小系数分布特征；S_n=2的高斯成分具有较大方差，刻画与图像边缘纹理对应的大系数的分布特征。小波系数x_n的概率分布为

P(x_n)=P(S_n=1)P_x|S(x_n|S_n=1)+P(S_n=2)P_x|S(x_n|S_n=2)(1)

其中，P(S_n=k),(k=1,2)表示状态变量S_n的概率分布，

P(S_n=1)+P(S_n=2)=1(2)

P_x|S(x_n|S_n=k)为x_n的高斯条件概率密度函数，表示如下：

$P_{x|S}\left(x_n\right|S_n=k)=N(0,{\mathrm{\σ}}_{k;n}^2)---\left(3\right)$

这里

小波系数的传递性可用隐状态变量沿尺度的状态转移概率表征：任一结点小波系数的隐状态变量仅依赖于其父结点的隐状态。Markov状态转移性可由一个2×2维的状态转移矩阵A_n确定

$A_n=\left(\begin{array}{cc}{p}_{11}& p_{12}\\ p_{21}& p_{22}\end{array}\right)---\left(4\right)$

其中p_mn=P(S=m|S_p=n)表示在父系数的状态为n的条件下，其子系数的状态为m的概率。

因此，HMT模型可由这些参数确定：混合方差转移矩阵A_n；根结点状态{P(S₁=1),P(S₁=2)}。为简化模型，同一尺度上的所有系数的模型参数可视为相等的，并且小波系数的幅值在尺度间呈指数形式下降，即HMT模型的在尺度j的高斯成分方差有以下关系：

${\mathrm{\σ}}_{k;j}^2=C_{{\mathrm{\σ}}_k}\·2^{-j{\mathrm{\α}}_k},k=\mathrm{1,2}---\left(5\right)$

从而可得到指数形式的状态转移概率矩阵

$A_j=\left(\begin{array}{cc}1-C_{11}{2}^{-{\mathrm{\γ}}_{1}j}& C_{11}-2^{-{\mathrm{\γ}}_{1}j}\\ \frac{1}{2}-C_{22}{2}^{-{\mathrm{\γ}}_{2}j}& \frac{1}{2}+C_{22}{2}^{-{\mathrm{\γ}}_{2}j}\end{array}\right)---\left(6\right)$

因此，HMT模型可由如下的参数向量确定：

$\mathrm{\&Theta;}=\{{\mathrm{\&alpha;}}_1,C_{{\mathrm{\&alpha;}}_1},{\mathrm{\&alpha;}}_2,C_{{\mathrm{\&alpha;}}_2},{\mathrm{\&gamma;}}_1,C_{11},{\mathrm{\&gamma;}}_2,C_{22},p_{j_0}\left(1\right),p_{j_0}\left(2\right)\},1<j<J---\left(7\right)$

其中k=1,2是根系数的状态概率。

模型参数可以通过EM算法训练得到，已知了小波系数x和模型参数Θ，可通过Upward-Downward算法获得小波系数的后验概率P(S=k|x,Θ)。

(2)基于HMT的迭代加权L1最小化算法

迭代加权L1方法（IRWL1）是通过求解加权L1的最小化问题以增强重构信号的稀疏性。在第i+1次迭代中，迭代加权L1最小化要求解如下形式的解：

$\left\{\begin{array}{cc}{\tilde{x}}^{(i+1)}\arg & \underset{x}{m}\mathrm{in}{\left|\right|W^{\left(i\right)}x\left|\right|}_1\\ \mathrm{subjectto}& \mathrm{\&Phi;\&Psi;x}=y\end{array},---\left(8\right)\right.$

其中，W⁽ⁱ⁾是对角权值矩阵，其对角元素

$w_n^{\left(i\right)}={\left(\right|{\tilde{x}}_n^{(i-1)}+\mathrm{\&epsiv;}\left|\right)}^{-1}---\left(9\right)$

其中ε是正则化常量。

由于迭代加权L1方法没有考虑小波系数中存在的结构信息，因此在基于HMT模型的迭代加权L1最小化算法中，将HMT模型参数引入权值的修正中，修正的权值为：

$w_n^{\left(i\right)}={(P(S_n=2{|\tilde{x}}^{(i-1)},\mathrm{\&Theta;})+\mathrm{\&delta;})}^{-p}---\left(10\right)$

其中，δ是正则化参数，p用来调节对小系数惩罚强度的参数。用该修正的权值代替IRWL1中的权值，加强了被重建的系数中存在的尺度间的树结构相关性，获得了优于标准IRWL1的重建性能。

HMT+IRWL1算法首先将初始权值W⁽⁰⁾=1的公式（8）的解做为初始解，然后利用EM算法对HMT的参数Θ进行估计，接着，执行权值更新和求解加权L1范数最小化重建。在权值更新步中，通过Upward-Downward算法获得小波系数的后验概率P(S=k|x,Θ)，再利用公式（10）对权值进行更新；在重建步中，将更新的权值用于公式（8），对其进行求解，获得系数的新的更新解当满足算法迭代停止准则时，算法停止。

上述重建算法结合了小波系数的先验结构和加权L1最小化算法，在重建精度上有一定优势，但是仍存在这样的缺点：

在对HMT模型参数Θ进行估计时，利用EM算法需经过大量的训练才能得到具体的参数值，计算量较大，算法计算开销大。

采样率接近50%时，重建精度略低于迭代加权L1最小化法，也就是说尽管在大多数情况下，HMT模型能精确刻画小波系数分布，但在个别时刻，小波系数存在排斥隐马尔可夫树模型的情况。

发明内容

为了解决HMT模型参数估计计算复杂度高的问题，结合自然图像自身特性，本发明采用uHMT模型代替HMT，该uHMT模型直接给出模型参数，不需要通过EM算法训练处模型参数，也即uHMT模型省去了计算量较大的训练步骤，有效降低了计算复杂度；同时，针对现有技术的第2个缺点，本发明给出的改进的权值计算方式，将HMT模型引入重建算法的同时又降低了小波系数偶尔出现排斥模型的可能，从而提高重建值的精确度，提高图像的重建质量。

本发明提出了一种基于通用隐马尔可夫树模型的加权图像压缩感知方法，其包括以下几个步骤：

步骤1、输入图像A,对其进行小波变换，得到A的小波系数矩阵X；

步骤2、通过高斯随机矩阵对所述小波系数矩阵X进行测量，得到测量值矩阵Y；

步骤3、对测量矩阵Y按列进行分块，每一块是一个一维列向量，记第k列为Y_k，分别对每一列进行如下处理；

步骤4、设置迭代次数为i,通过求解P1问题对Y_k进行重建，得到重建初始值

步骤5、根据第i-1次迭代得到的解和uHMT模型参数Θ计算第i次更新的权值

步骤6、将更新的权值用于WP1问题，求解获得第i次的更新解

步骤7、当迭代次数达到给定的最大值时，迭代停止，得到第k列的重建稀疏值否则，重复步骤5和6，直至满足迭代停止准则；

步骤8、k=k+1,重复步骤4至7，直至矩阵Y的所有列都得以重建。将所有列的重建稀疏值重组成稀疏值矩阵

步骤9、对稀疏值进行小波逆变换，得到图像重建值。

其中，所述步骤5为通过下式计算更新的权值

其中，N为重建向量的长度。

对n∈[1,N]，中的对角元素计算如下，

$w_{n,n}^{\left(i\right)}=\left\{\begin{array}{cc}{\left(\right|{\tilde{x}}_{k,n}^{(i-1)}|+\mathrm{\&epsiv;})}^{-1}& ,i=2L-1\\ {(P(S_n=2|{\tilde{X}}_k^{(i-1)},\mathrm{\&Theta;})+\mathrm{\&delta;})}^{-q}& ,i=2L\end{array}\right.$

其中，L为正整数，表示的第n个元素；ε和δ均为正则化参数，ε∈[5,20]，δ∈[10^-1,10^-10]；q是用来调节对小系数惩罚强度的参数，q∈[0.1,1]；状态概率可用Viterbi算法求解得到。

其中，步骤5所述uHMT模型参数为

$\mathrm{\&Theta;}=\{{\mathrm{\&alpha;}}_1=3.1,C_{{\mathrm{\&alpha;}}_1}=2^{11},{\mathrm{\&alpha;}}_2=2.25,C_{{\mathrm{\&alpha;}}_2}=2^{11},{\mathrm{\&gamma;}}_1=1,C_{11}=2^{2.3},{\mathrm{\&gamma;}}_2=0.4,C_{22}=2^{0.5},p_{j_0}\left(1\right)=0.5,p_{j_0}\left(2\right)=0.5\}$

其中，α_k,（k=1,2）是用于决定模型中位于尺度j的小波系数混合方差的参数，γ_k，C_kk，（k=1,2）是决定状态转移矩阵的参数，(k=1,2)表示uHMT模型中根系数的状态概率。

其中所述步骤4中的P1问题表示如下：

$\left\{\begin{array}{cc}\tilde{x}=\mathrm{ar}\underset{x}{g}& \min {\left|\right|x\left|\right|}_1\\ \mathrm{subjectto}& \mathrm{\&Phi;\&Psi;x}=y\end{array}\right.$

其中，Φ是测量矩阵，Ψ是小波变换矩阵，y是测量值。argmin表示使得||x||₁最小的x的值，||x||₁表示x的1-范数，subjectto表示x满足这样一个约束条件。

其中步骤6所述WP1问题表示如下：

$\left\{\begin{array}{cc}\tilde{x}=\mathrm{ar}\underset{x}{g}& \min {\left|\right|W_x\left|\right|}_1\\ \mathrm{subjectto}& \mathrm{\&Phi;\&Psi;x}=y\end{array}\right.$

其中，W是权值矩阵，Φ是测量矩阵，Ψ是小波变换矩阵，y是测量值，argmin表示使得||x||₁最小的x的值，||x||₁表示x的1-范数，subjectto表示x满足这样一个约束条件。

本发明将基于uHMT模型的迭代加权L1最小化算法应用于图像的压缩感知重建。本发明根据自然图像的自身特性，采用uHMT模型代替HMT模型，直接给出模型参数，能够省去计算量较大的训练步骤，有效降低图像重建的计算复杂度。同时，由于权值的计算方法对重建精度有很大影响，通过本发明给出的改进的权值计算方法既将刻画小波系数的uHMT模型引入了图像重建，又降低了出现排斥模型现象的可能性，从而提高了图像的重建精度。

附图说明

图1为现有技术中图像的小波四叉树结构。

图2为与本发明实施例一致的方法流程图。

图3为当采样率为50%时，通过与本发明实施例一致的方法与GPSR，HMT+IRWL1重建出来的Lena图的仿真对比图。

其中，图（a）为原始图像；（b）为GPSR重建图像，（c）为HMT+IRWL1重建图像；（d）为通过本发明的方法重建的图像。

图4为采用本发明的方法与GPSR，HMT+IRWL1重建出来的Lena图的峰值信噪比PSNR和重构时间t随采样率变化的对比表。

具体实施例

本发明技术适用于自然图像的压缩感知重建。将uHMT模型引入迭代加权L1方法中，改进权值计算方式，克服了现有技术的缺点，重建出原始图像。如图2所示，本发明提出一种基于通用隐马尔可夫树模型的加权图像压缩感知方法，其包括以下几个步骤组成：

步骤1、输入图像A,对其进行小波变换，得到A的小波系数矩阵X;

步骤2、通过高斯随机矩阵对X进行测量得到测量值矩阵Y；

步骤3、对测量矩阵Y按列进行分块，每一块是一个一维列向量，记第k列为Y_k，分别对每一列进行如下处理；

步骤4、设置迭代次数为i,通过求解P1问题对Y_k进行重建得到初始值 ${\tilde{X}}_k^{\left(0\right)};$

步骤5、根据第i-1次迭代得到的解和uHMT模型参数Θ计算更新的权值W⁽ⁱ⁾；

步骤6、将更新的权值用于WP1问题，求解获得第i次更新解

步骤7、当迭代次数达到给定的最大值时，迭代停止，得到重建稀疏值否则，重复步骤5和6，直至满足迭代停止准则；

步骤8、k=k+1,重复步骤4至7，直至矩阵Y的所有列都得以重建。将所有列的重建稀疏值重组成稀疏值矩阵

步骤9、对稀疏值进行小波逆变换，得到图像重建值。

2.2.1uHMT模型参数

现有技术主要针对的是分段光滑的实信号，模型参数需通过EM算法经过大量的训练才能得到具体的参数值，而本发明技术针对的是自然图像信号，，大多数自然图像信号的归一化小波系数的方差衰减以及传递性衰减范围非常相似，根据这一特性，本发明用uHMT模型代替HMT模型，该模型直接给出模型参数，能省去HMT模型求参数时较大的训练步骤。uHMT模型参数为

$\mathrm{\&Theta;}=\{{\mathrm{\&alpha;}}_1=3.1,C_{{\mathrm{\&alpha;}}_1}=2^{11},{\mathrm{\&alpha;}}_2=2.25,C_{{\mathrm{\&alpha;}}_2}=2^{11},{\mathrm{\&gamma;}}_1=1,C_{11}=2^{2.3},{\mathrm{\&gamma;}}_2=0.4,C_{22}=2^{0.5},p_{j_0}\left(1\right)=0.5,p_{j_0}\left(2\right)=0.5\}---\left(11\right)$

其中，α_k,（k=1,2）是用于决定模型中位于尺度j的小波系数混合方差的参数，γ_k，C_kk，（k=1,2）是决定状态转移矩阵的参数，(k=1,2)表示uHMT模型中根系数的状态概率。

该模型参数经验证对于大多数自然图像都是适用的。

2.2.2更新权值的计算

本发明定义了如下的权值更新计算公式，权值矩阵是一个N×N的对角矩阵，其对角元素计算公式如下：

$w_{n,n}^{\left(i\right)}=\left\{\begin{array}{cc}{\left(\right|{\tilde{x}}_{k,n}^{(i-1)}|+\mathrm{\&epsiv;})}^{-1}& ,i=2L-1\\ {(P(S_n=2|{\tilde{X}}_k^{(i-1)},\mathrm{\&Theta;})+\mathrm{\&delta;})}^{-q}& ,i=2L\end{array}\right.$

其中，n∈[1,N]，L为正整数，表示的第n个元素；ε和δ均为正则化参数，ε∈[5,20]，δ∈[10^-1,10^-10]，q是用来调节对小系数惩罚强度的参数，q∈[0.1,1]；状态概率可用Viterbi算法求解得到。

新的权值既将刻画小波系数的uHMT模型引入了图像重建，又降低了出现排斥模型现象的可能性。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

1.一种基于通用隐马尔可夫树模型的加权图像压缩感知方法，其包括以下步骤：

步骤1、输入图像A,对其进行小波变换，得到A的小波系数矩阵X；

步骤2、通过高斯随机矩阵对所述小波系数矩阵X进行测量，得到测量值矩阵Y；

步骤3、对测量值矩阵Y按列进行分块，每一块是一个一维列向量，记第k列为Y_k，分别对每一列进行如下处理；

步骤4、设置迭代次数为i,通过求解P1问题对Y_k进行重建，得到重建初始值

步骤5、根据第i-1次迭代得到的解和uHMT模型参数Θ计算第i次更新的权值

所述步骤5中更新的权值的计算公式如下：

其中，N为重建向量的长度；

对n∈[1,N]，中的对角元素计算如下，

$w_{n,n}^{\left(i\right)}=\left\{\begin{array}{cc}{\left(\right|{\tilde{x}}_{k,n}^{(i-1)}|+\mathrm{\&epsiv;})}^{-1},& i=2L-1\\ {\left(P\right(S_n=2|{\tilde{X}}_k^{(i-1)},\mathrm{\&Theta;})+\mathrm{\&delta;})}^{-q},& i=2L\end{array}\right.$

其中，L为正整数，表示的第n个元素；ε和δ均为正则化参数，ε∈[5,20]，δ∈[10^-1,10^-10]，q是用来调节对小系数惩罚强度的参数，q∈[0.1,1]；状态概率可用Viterbi算法求解得到；

步骤6、将更新的权值用于WP1问题，求解获得第i次更新解

步骤7、当迭代次数达到给定的最大迭代次数时，迭代停止，得到第k列的重建稀疏值否则，重复步骤5和6，直至满足迭代停止准则；

步骤8、k＝k+1,重复步骤4至7，直至矩阵Y的所有列都得以重建；将所有列的重建稀疏值重组成稀疏值矩阵

步骤9、对稀疏值矩阵进行小波逆变换，得到图像重建值。

2.如权利要求1所述基于通用隐马尔可夫树模型的加权图像压缩感知方法，

其中所述步骤4中的P1问题表示如下：

$\left\{\begin{array}{c}\tilde{x}=\underset{x}{\arg }\min {||x||}_1\\ \begin{array}{cc}subjectto& \mathrm{\&Phi;}\mathrm{\&Psi;}x=y\end{array}\end{array}\right.$

其中，Φ是测量矩阵，Ψ是小波变换矩阵，y是测量值，argmin表示使得||x||₁最小的x的值，||x||₁表示x的1范数，subjectto表示x满足这样一个约束条件。

3.如权利要求1所述基于通用隐马尔可夫树模型的加权图像压缩感知方法，其中步骤5所述uHMT模型参数为：

其中，α_k,k＝1,2是用于决定模型中位于尺度j的小波系数混合方差的参数；γ_k，C_kk，k＝1,2是决定状态转移矩阵的参数；，k＝1,2表示uHMT模型中根系数的状态概率。

4.如权利要求1所述基于通用隐马尔可夫树模型的加权图像压缩感知方法，其中步骤6所述WP1问题表示如下：

$\left\{\begin{array}{c}\tilde{x}=\underset{x}{\arg }\min {||Wx||}_1\\ \begin{array}{cc}subjectto& \mathrm{\&Phi;}\mathrm{\&Psi;}x=y\end{array}\end{array}\right.$

其中，W是权值矩阵，Φ是测量矩阵，Ψ是小波变换矩阵，y是测量值，argmin表示使得||x||₁最小的x的值，||x||₁表示x的1-范数，subjectto表示x满足这样一个约束条件。