一种基于伪牛顿法的四元数域彩色图像压缩感知恢复方法

技术领域

本发明涉及一种基于伪牛顿法的四元数域彩色图像压缩感知恢复方法，属于数字图像处理技术领域。

背景技术

压缩感知(Compressed Sensing:CS)理论指出：只要信号是可压缩的或在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原始信号。其理论框架如图1所示，图中CS信息矩阵A=ΦΨ*，其中Φ为测量矩阵，Ψ为稀疏矩阵，上标“*”表示共轭转置。由图可知：

(1)CS的最终目的是用尽可能少的测量数据来恢复原始信号。

(2)CS理论框架主要分为三步：信号的稀疏表示、信号测量和信号恢复。其中稀疏矩阵Ψ可以根据信号自身的特点在正交基或字典中灵活选取。因此关于CS的研究主要集中在后两个步骤：测量矩阵Φ的设计和快速鲁棒的信号恢复算法。

对于图1中第三步，信号的恢复算法可以等价于求解如下优化问题：

$>\hat{x}=\arg \min J\left(x\right)$>$>s.t.J\left(x\right)={\left|\right|y-\mathrm{Ax}\left|\right|}_2^2+{\lambda }_1{\left|\right|{\psi }^{*}x\left|\right|}_0,---\left(1\right)$>

上式求解是个NP-Hard问题，由于l₁范数在一定条件下和l₀范数具有等价性，可得到相同的解。那么上式可以转化为l₁范数下的最优化问题：

$>\hat{x}=\arg \min J\left(x\right)$>$>s.t.J\left(x\right)={\left|\right|y-\mathrm{Ax}\left|\right|}_2^2+{\lambda }_1{\left|\right|{\psi }^{*}x\left|\right|}_1,---\left(2\right)$>

l₁范数下最优化问题又称为基追踪(Basis Pursuit:BP)，其常用实现算法有：伪牛顿法、内点法和梯度投影法。

四元数q是一种特殊的超复数，是复数的扩展形式。四元数q由1个实部和3个虚部组成：

q＝S(q)+X(q)i+Y(q)j+Z(q)k,    (3)

其中S(q),X(q),Y(q),表示实数域，i,j,k是三个虚数单位，满足如下性质：

i²＝j²＝k²＝ijk＝-1,    (4)

ij＝-ji＝k,jk＝-kj＝i,ki＝-ik＝j.    (5)

四元数的共轭和模分别定义为：

$>\bar{q}=S\left(q\right)-X\left(q\right)i-Y\left(q\right)j-Z\left(q\right)k,---\left(6\right)$>

$>\left|q\right|=\sqrt{S^2\left(q\right)+X^2\left(q\right){+Y}^2\left(q\right)+Z^2\left(q\right)}.---\left(7\right)$>

四元数q也可以用欧拉形式表达为幅度和相位的形式：

$>q=S\left(q\right)+X\left(q\right)i+Y\left(q\right)j+Z\left(q\right)k=\left|q\right|e^{{\eta }^{\theta }},$>

$>\left|q\right|=\sqrt{S^2\left(q\right)+X^2\left(q\right)+Y^2\left(q\right)+Z^2\left(q\right)}$>

$>\eta =\frac{X\left(q\right)i+Y\left(q\right)j+Z\left(q\right)k}{\sqrt{X^2\left(q\right)+Y^2\left(q\right)+Z^2\left(q\right)}},---\left(8\right)$>

θ＝arccos(S(q)/|q|).

η表示归一化的纯四元数，θ表示四元数相位角。

类似地，(3)和(8)的定义形式也适用于四元数矩阵x。

四元数矩阵x的l_p范数为：

$>{\left|\right|x\left|\right|}_p={\left(\underset{r=1}{\overset{m}{\Sigma }}{\left|x_r\right|}^p\right)}^{1/p}0<p<2---\left(9\right)$>

目前压缩感知对彩色图像的处理均采用传统的RGB三个分量单独处理的方式，而没有考虑三个通道数据之间的相关性，使得恢复算法的约束信息不够完善，恢复出来的图像效果没有达到理想的效果。

发明内容

发明目的：针对传统彩色图像压缩感知恢复算法中的将RGB三个通道数据之间的相关性割裂的问题与不足，本发明结合压缩感知理论和四元数信号处理理论，提供一种基于伪牛顿法的四元数域彩色图像压缩感知恢复方法。将彩色图像转化为四元数矩阵信号，采用四元数矩阵信号的欧拉表示形式，在压缩感知恢复算法中加入了四元数矩阵信号的幅度和相位作为新的约束项，使恢复的图像更为平滑。

为了将彩色图像三个通道数据之间的相关性相结合，将彩色图像的RGB三个分量分别放在四元数的一个实部和两个虚部，这样将彩色图像中的信息转化为四元数矩阵。同时我们将四元数信号表示成欧拉形式，即幅度和相位的形式，这样我们就保留了三个通道数据之间的相关性，从而能更好的恢复彩色图像。

技术方案：一种基于伪牛顿法的四元数域彩色图像压缩感知恢复方法，包括编码端对彩色图像进行稀疏表示及测量编码，解码端对接收的数据进行解码重构，得到恢复图像。具体包括以下步骤：

步骤A、将大小为N×N的彩色图像T的RGB三个通道的数据分别赋值给N×N的四元数矩阵Q的一个实部和两个虚部，并将第三个虚部置零，即

Q=R(T)+G(T)i+B(T)j+0k；i,j,k是三个虚数单位，R(T)，G(T)和B(T)分别表示彩色图像R分量、G分量和B分量的数据值；

步骤B、用大小为N×N的稀疏矩阵W对四元数矩阵Q进行稀疏，得到稀疏后的四元信号矩阵Q₀，稀疏矩阵W为离散小波变换(DWT)矩阵；

步骤C、用大小为M×N的观测矩阵Φ对稀疏后的四元信号矩阵Q₀进行观测，得到M×N的观测值Y₁，观测矩阵Φ为随机高斯矩阵；(注：对四元信号矩阵Q₀进行观测相当于按每列进行观测处理)其中，M＜＜N；

步骤D、用观测矩阵Φ对稀疏后的四元信号矩阵Q₀的转置进行观测，得到M×N的观测值Y₂；(注：对四元信号矩阵Q₀的转置进行观测相当于对四元信号矩阵Q₀按每行进行观测处理)

步骤E、设置压缩感知优化方程迭代的初始值Q_t，Q_t＝Φ*×Y(这里的Y分别指Y₁和Y₂)，上标“*”表示共轭转置，并且将得到的四元数信号矩阵Q_t改写为幅度和相位的形式：η表示归一化的纯四元数，θ_t表示四元数相位角；表示两个矩阵的点乘积，即矩阵对应位置的元素相乘；

步骤F、用观测值Y₁和Y₂，以及稀疏矩阵W和观测矩阵Φ，结合四元数信号压缩感知恢复算法，从观测值中高概率重构出四元数信号和其中是通过观测值Y₁重构得到的信号，是通过观测值Y₂重构得到的信号；

步骤G、将得到的四元数信号和相加求平均，得到四元数矩阵将的实部赋值给彩色图像的R通道，将的两个虚部分别赋值给彩色图像的G通道和B通道，得到恢复的彩色图像。

将彩色图像T的RGB三个通道的数据分别赋值给大小为N×N四元数矩阵Q的一个实部和两个虚部过程中，两个虚部的选择可以是任意的，但一定要赋值给实部，保证相位有意义。

所述从观测值中高概率重构出四元数信号和的步骤为：

步骤1)、将压缩感知恢复算法中经典的l₁范数最优化问题

$>\hat{x}=\arg \min J\left(x\right)$>$>s.t.J\left(x\right)={\left|\right|y-\mathrm{Ax}\left|\right|}_2^2+{\lambda }_1{\left|\right|{\psi }^{*}x\left|\right|}_1$>

改写为l_p(其中0<p<2)范数下的最优化问题，同时用幅度和相位作为压缩感知优化问题新的约束项，得到目标函数J(Q)；

$>\hat{Q}=\arg \min J\left(Q\right)$>

$>s.t.J\left(Q\right)={\left|\right|Y-\Phi {\mathrm{WQW}}^{\text{'}}|}_2^2+{\lambda }_1\underset{i=1}{\overset{N}{\Sigma }}{{({\left({\mathrm{WQW}}^{\text{'}}\right)}_i^2+ϵ)}^{p/2}}^{}+{\lambda }_2\underset{i=1}{\overset{N}{\Sigma }}{({\left(D\right|Q\left|\right)}_i^2+ϵ)}^{p/2}+{\lambda }_3\underset{i=1}{\overset{N}{\Sigma }}{({\left(\mathrm{D\theta }\right)}_i^2+ϵ)}^{p/2}$>

其中ε是一个很小的常数，常取1.0×10^-6，D代表梯度算子，W′表示W的转置；

步骤2)、对目标函数J(Q)关于四元数信号Q求导，得到目标函数的梯度ΔJ(Q)，

通过测量梯度的变化构造一个近似的海森矩阵H(Q)；

步骤3)、设置迭代的步长因子γ，γ取值为0.7；

步骤4)、用梯度ΔJ(Q_t)，海森矩阵H(Q_t)及步长因子γ，通过公式Q₁=Q_t-γ[H(Q_t)]^-1ΔJ(Q_t)求得经过第一次迭代后的值Q₁；

步骤5)、将Q_t更新成Q₁，重新计算梯度ΔJ(Q_t)和海森矩阵H(Q_t)，重复步骤I和J，直到迭代过程收敛；

步骤6)、如果观测值为Y₁时，收敛后的Q₁不再赋给Q_t，把Q₁赋给作为Y₁的处理结果；

如果观测值为Y₂时，收敛后的Q₁不再赋给Q_t；把Q₁赋给作为Y₂的处理结果。

有益效果：本发明的彩色图像压缩感知恢复算法，在编码时，将彩色图像的RGB三个分量分别赋值给四元数信号的一个实部和两个虚部，将彩色图像的压缩感知恢复问题转换到四元数域的压缩感知恢复问题，对稀疏后四元数矩阵信号Q₀按行进行观测得到观测值Y₁，对Q₀按列进行观测得到观测值Y₂；进行解码时，通过观测值Y₁和Y₂，以及稀疏矩阵W和观测矩阵Φ，得到解码重构恢复的四元数矩阵和将四元数矩阵和相加求平均，得到结果四元数矩阵将四元数矩阵的实部及两个虚部，分别赋值给彩色图像的RGB三个通道得到结果图像。本发明的方法与现有的将彩色图像RGB三个分量数据单独处理的方法相比，只需进行一次压缩感知就能恢复彩色图像，比传统的进行三次压缩感知更节约时间。同时本发明利用了四元数矩阵的幅度和相位信息，结合了三个通道之间的相关性，使得彩色图像的恢复更为精确，这是传统的彩色图像处理技术难以企及的。

附图说明

图1为压缩感知的理论框架图；

图2为本发明方法与现有方法的比较流程图；

图3为采样率从60%到85%时，按列进行处理时lena图像的恢复效果比较图；

图4为采样率从60%到85%时，按行进行处理时lena图像的恢复效果比较图；

图5为采样率从60%到85%时，按行处理与按列处理相结合时lena图像的恢复效果比较图；

图6为采样率从60%到85%时，单独按行处理、单独按列处理、按行处理与按列处理相结合的方法所对应的图像的不同恢复效果比较图。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图2所示，基于伪牛顿法的四元数域彩色图像压缩感知恢复方法，按照前述步骤进行：

步骤1、将彩色图像表示成四元数矩阵的形式，一幅64×64彩色lena图像可以用一个64×64的四元数矩阵来存储其三个通道的值。

例如，彩色图像右上角3×3的像素矩阵是：

$>\left(\begin{array}{ccc}\left(\mathrm{151,111,108}\right)& \left(\mathrm{84,37,34}\right)& \left(\mathrm{105,65,66}\right)\\ \left(\mathrm{130,102,98}\right)& \left(\mathrm{59,28,24}\right)& \left(\mathrm{80,62,57}\right)\\ (\mathrm{93,99},161)& \left(\mathrm{20,29,65}\right)& \left(\mathrm{55,62,78}\right)\end{array}\right)$>

则相对应的四元数矩阵表示为：

$>\left(\begin{array}{ccc}(151+111i+108j+0k)& (84+37i+34j+0k)& (105+65i+66j+0k)\\ (130+102i+98j+0k)& (59+28i+24j+0k)& (80+62i+57j+0k)\\ (93+99i+161j+0k)& (20+29i+65j+0k)& (55+62i+78j+0k)\end{array}\right)$>

因此，我们可以得到大小为64×64的四元数矩阵Q。

步骤2、将四元数矩阵Q在冗余字典下进行稀疏分解，实验中我们取的是大小为64×64离散小波变换(DWT)矩阵W，根据以下公式进行稀疏分解：

Q₀＝W×Q×W′,W′表示矩阵W的转置.    (10)

例如，彩色图像右上角3×3的像素矩阵稀疏后得到：

$>\left(\begin{array}{ccc}(4-2i-j+0k)& (6-i-6j+0k)& (4-i-3j+0k)\\ (-2+0i-j+0k)& (-6+i+0j+0k)& (-3-i+2j+0k)\\ (-2-i-j+0k)& (5-i+0j+0k)& (2-i-j+0k)\end{array}\right)$>

对得到的3×3四元数矩阵中四元数幅度的值接近0的元素，我们将其置零，得到如下更为稀疏的四元数矩阵：

$>\left(\begin{array}{ccc}(4-2i-j+0k)& (6-i-6j+0k)& (4-i-3j+0k)\\ (0+0i+0j+0k)& (-6+i+0j+0k)& (0+0i+0j+0k)\\ (0+0i+0j+0k)& (5-i-0j+0k)& (0+0i+0j+0k)\end{array}\right)$>

因此，我们可以得到大小为64×64的稀疏四元数矩阵Q₀，并且Q₀中的非零元素远远少于N×N。

步骤3、用大小为M×N的随机高斯矩阵Φ对稀疏后的四元数信号Q₀的按行进行观测，得到M×N的观测值Y₁，接着对Q₀按列进行观测，得到M×N的观测值Y₂。

步骤401、对原始四元数矩阵Q进行压缩感知编码以后，我们需要通过解决如下优化问题来解码恢复彩色图像：

$>\hat{Q}=\arg \min J\left(Q\right)$>

$>s.t.J\left(Q\right)={\left|\right|Y-\Phi {\mathrm{WQW}}^{\text{'}}\left|\right|}_2^2+{\lambda }_1{\left|\right|{\mathrm{WQW}}^{\text{'}}\left|\right|}_p^p+{\lambda }_2{\left|\right|D\left|Q\right||}_p^p+{\lambda }_3{\left|\right|\mathrm{D\theta }\left|\right|}_p^p---\left(11\right)$>

其中D代表梯度算子，形式如下：

$>D=\left(\begin{array}{ccccccc}-1& 1& 0& 0& .& .& 0\\ 0& -1& 1& 0& .& .& 0\\ .& .& .& .& .& .& .\\ 0& 0& 0& 0& 0& -1& 1\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right)---\left(12\right)$>

步骤402、为了避免l_p范数在初始点的不可微，我们将l_p范数改写为如下形式：

$>{\left|\right|f\left|\right|}_p^p\approx \underset{i=1}{\overset{N}{\Sigma }}{({\left|{\left(f\right)}_i\right|}^2+\xi )}^{p/2}---\left(13\right)$>

其中ξ是一个很小的常数，常取1.0×10^-6。这样(11)式就可以改写为：

$>\hat{Q}=\arg \min J\left(Q\right)$>

$>s.t.J\left(Q\right)={\left|\right|Y-\Phi {\mathrm{WQW}}^{\text{'}}|}_2^2+{\lambda }_1\underset{i=1}{\overset{N}{\Sigma }}{{({\left({\mathrm{WQW}}^{\text{'}}\right)}_i^2+ϵ)}^{p/2}}^{}+{\lambda }_2\underset{i=1}{\overset{N}{\Sigma }}{({\left(D\right|Q\left|\right)}_i^2+ϵ)}^{p/2}+{\lambda }_3\underset{i=1}{\overset{N}{\Sigma }}{({\left(\mathrm{D\theta }\right)}_i^2+ϵ)}^{p/2}---\left(14\right)$>

其中ε是一个很小的常数，常取1.0×10^-6，D代表梯度算子，W′表示W的转置；

步骤403、对目标函数J(Q)求导，得到目标函数的梯度ΔJ(Q)，通过测量梯度的变化构造一个近似的海森矩阵H(Q)。

步骤404、设置迭代的初始值Q_t＝Φ*×Y(这里的Y分别指Y₁和Y₂)，通过如下公式计算Q_t的幅度|Q_t|和相位θ_t：

$>Q=S\left(Q\right)+X\left(Q\right)i+Y\left(Q\right)j+0k=\left|Q\right|e^{{\eta }^{\theta }}$>

$>\left|Q\right|=\sqrt{S^2\left(Q\right)+X^2\left(Q\right)+Y^2\left(Q\right)}$>

$>\eta =\frac{I\left(Q\right)i+J\left(Q\right)j}{\sqrt{I^2\left(Q\right)+J^2\left(Q\right)}}---\left(15\right)$>

θ=arccos(S(Q)/|Q|)

例如将得到的64×64的四元数矩阵Q_t的右上角3×3的元素改写为幅度和相位的形式：

$>\left|Q_t\right|=\left(\begin{array}{ccc}202.1930& 133.8469& 131.0572\\ 163.6490& 122.6051& 115.8836\\ 109.8818& 120.3578& 190.3418\end{array}\right)$>

$>{\theta }_t=\left(\begin{array}{ccc}0.7275& 0.5929& 0.6022\\ 0.6528& 0.5882& 0.5630\\ 0.5617& 0.6049& 0.5626\end{array}\right)$>

步骤405、设置迭代的步长因子γ，γ取值常设为(0,1]之间，本发明γ取值为0.7。

步骤406、结合梯度ΔJ(Q)和海森矩阵H(Q)用如下伪牛顿迭代式求得经过第一次迭代的值Q₁。

Q₁=Q_t—γ[H(Q_t)]^-1ΔJ(Q_t)    (16)

步骤407、将Q_t更新成Q₁，重新计算梯度ΔJ(Q)和海森矩阵H(Q)，重复步骤406到步骤407，直到迭代过程收敛。

步骤5、将通过观测值Y₁和Y₂重构得到的四元数矩阵和相加求平均，得到四元数矩阵将四元数矩阵的实部及两个虚部，分别赋值给彩色图像的RGB三个通道得到结果图像T。

为了验证本发明方法的效果，进行了以下实验：

1、实验条件：

在一台计算机上进行验证实验，该计算机配置为AMD单核处理器(3200兆赫兹)和4096兆字节随机存取存储器(RAM)，编程语言用的是Matlab(8.0版本)。

2、实验方法：

使用大小为64×64的彩色图像lena作为实验对象，采用以下步骤完成lena图像的压缩感知恢复。

1)采用具体实施方式中的步骤1，将彩色图像二维矩阵转化为四元数二维矩阵Q；

2)采用具体实施方式中的步骤2，对四元数矩阵Q进行稀疏，得到稀疏后的四元数矩阵Q₀；

3)设置不同的采样率，对四元数矩阵Q₀的按行进行观测，得到观测值Y₁。将公式(14)中的三个约束项的权重设置为

λ₁=0.01，λ₂=0.0114，λ₃=0.0028，p＝0.9，迭代的步长因子r=0.7，用具体实施方式中的步骤401到步骤406，得到恢复矩阵

4)设置不同的采样率，对四元数矩阵Q₀的按列进行观测，得到观测值Y₂。将公式(14)中的三个约束项的权重设置为

λ₁=0.005，λ₂=0.0124，λ₃=0.0048，p＝0.9，迭代的步长因子r=0.7，用具体实施方式中的步骤401到步骤406，得到恢复矩阵

5)采用具体实施方式中的步骤5，得到结果图像T。

3、实验结果的评价指标：

实验中采用如下公式计算彩色图像的峰值信噪比：

两个M×N的彩色图像I和K，它们RGB各个通道的均方差定义为：

$>{\mathrm{MSE}}_r=\frac{1}{\mathrm{MN}}\underset{i=0}{\overset{M-1}{\Sigma }}\underset{j=0}{\overset{N-1}{\Sigma }}{\left|\right|I(i,j,1)-K(i,j,1)\left|\right|}^2$>

$>{\mathrm{MSE}}_g=\frac{1}{\mathrm{MN}}\underset{i=0}{\overset{M-1}{\Sigma }}\underset{j=0}{\overset{N-1}{\Sigma }}{\left|\right|(i,j,2)-K(i,j,2)\left|\right|}^2---\left(17\right)$>

$>{\mathrm{MSE}}_b=\frac{1}{\mathrm{MN}}\underset{i=0}{\overset{M-1}{\Sigma }}\underset{j=0}{\overset{N-1}{\Sigma }}{\left|\right|(i,j,3)-K(i,j,3)\left|\right|}^2$>

RGB各个通道的峰值信噪比定义为：

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

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

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

其中，MAX_r,MAX_g,MAX_b表示图像I和K的RGB各个通道图像点颜色的最大值，如果每个采样点用8位表示，那么它们就是255。

彩色图像I和K的峰值信噪比定义为：

PSNR＝(PSNR_r+PSNR_g+PSNR_b)/3    (19)

图3显示了采样率从60%到85%时，按列进行处理时lena图像的恢复效果。表1为采样率与峰值信噪比。

表1采样率与峰值信噪比

单独按列进行处理时，原始图像的行和行之间的相关性被割裂了。从图3的恢复结果可以看到，恢复图像的列与列之间存在明显的裂痕。当采样率提高时，这种裂痕的数量变少了，说明在提高采样率时，每列的恢复效果更为精确，可以近似地保持行与行之间相关性，但裂痕仍然会存在。

从表1的结果可以看到，当采样率提高时，恢复图像与原始图像的峰值信噪比提高了，这与压缩感知的理论知识是一致的。

图4显示了采样率从60%到85%时，按行处理时lena图像的恢复效果。表2为采样率与峰值信噪比。

表2采样率与峰值信噪比

单独按行进行处理时，原始图像的列和列之间的相关性被割裂了，从图5的恢复结果可以看到，恢复图像的行与行之间存在明显的裂痕。同单独按行处理一样，当采样率提高时，这种裂痕的数量变少了，说明在提高采样率时，每行的恢复效果更为精确，可以近似地保持列与列之间相关性，但裂痕仍然会存在。

从表2的结果可以看到，当采样率提高时，恢复图像与原始图像的峰值信噪比提高了，这与压缩感知的理论知识是一致的。对比表1和表2我们发现，单独按行进行处理比单独按列进行处理得到的恢复图像与原始图像的峰值信噪比要高，这说明按行进行处理时，每行的稀疏度比每列的稀疏度要大，这一点与压缩感知的理论知识也是一致的。

图5显示了采样率从60%到85%时，按行处理与按列处理相结合图像的恢复效果。表3为采样率与峰值信噪比。

使用按行处理与按列处理相结合的方法，从图5的恢复结果可以看到，恢复图像行与行之间的裂痕、恢复图像列与列之间的裂痕得到了明显的改善，在采样率提高时，恢复图像更为平滑，取得了较理想的结果。

从表3的结果可以看到，当采样率提高时，恢复图像与原始图像的峰值信噪比提高了，这与压缩感知的理论知识是一致的。对比表2和表3我们发现，当采用按行处理与按列处理相结合这种方法时，恢复图像与原始图像的峰值信噪比得到了小幅度的改善，这与我们预计的结果是一致的。

为了更好的对比不同方法的恢复结果，图6显示了采样率从60%到85%时，单独按行处理、单独按列处理、按行处理与按列处理相结合的方法所对应的图像的不同恢复效果。