发布时间: 2017-04-16 摘要点击次数: 全文下载次数: DOI: 10.11834/jig.20170403 2017 | Volume 22 | Number 4 图像处理和编码

Lp范数压缩感知图像重建优化算法

1. 南昌航空大学无损检测技术教育部重点实验室, 南昌 330063;
2. 南昌航空大学信息工程学院, 南昌 330063
 收稿日期: 2016-11-18; 修回日期: 2016-12-28 基金项目: 国家自然科学基金项目（61663030，61663032）；江西省自然科学基金项目（20142BAB207021）；江西省教育厅科技项目（GJJ150753）；无损检测技术教育部重点实验室（南昌航空大学）开放基金项目（ZD29529005） 第一作者简介: 蒋沅 (1982-), 男, 副教授, 硕士生导师, 2010年于山东大学获控制理论与控制工程专业博士学位, 主要研究方向为控制理论与应用的研究、压缩感知与图像处理。E-mail:jiangyuan@nchu.edu.cn 中图法分类号: TP301.6 文献标识码: A 文章编号: 1006-8961(2017)04-0435-08

# 关键词

Improved search algorithm for compressive sensing image recovery based on Lp norm
Jiang Yuan1,2, Miao Shengwei1,2, Luo Huazhu1,2, Shen Pei1,2
1. Key Laboratory of Nondestructive Testing (Nanchang Hangkong University), Ministry of Education, Nanchang 330063, China;
2. Jiangxi province key laboratory of image processing and pattern recognition, Nanchang 330063, China
Supported by: National Natural Science Foundation of China (61663030, 61663032);Natural Science Foundation of Jiangxi Province, China (20142BAB207021)

# Abstract

Objective As one of the key technologies in compression sensing, a reconstruction algorithm plays a key role in scientific research.Commonly used reconstruction algorithms include the non-convex optimization algorithm based on L0 norm and the convex optimization algorithm based on L1 norm. However, these algorithms exhibit shortcomings, such as low precision reconstruction and extremely long operation time. To overcome these limitations and improve the precision and efficiency of the existing compressed sensing image reconstruction algorithm based on Lp norm reconstruction, an improved algorithm is presented in this study. Method Sequential quadratic programming is one of the most effective methods for solving constrained nonlinear programming problems in recent years. Compared with other algorithms, it demonstrates the most apparent advantages of good convergence, high computational efficiency, and edge search capability. The use of the Lagrange function series quadratic programming (SQP) method for the heather (Hesse) does not provide a positive definite matrix calculation for many problems; therefore, a value function should be introduced, the Hesse matrix sequence quadratic programming method should be corrected, and an image block compressed sensing technology should be integrated. This study proposes an image reconstruction algorithm based on Lp norm compression awareness. Result Under the same sampling rate of 40%, the signal-to-noise ratio (SNR) of the proposed algorithm is 34.28 dB, which is better than that of the block orthogonal matching pursuit (BOMP) algorithmThe BOMP algorithm's SNR is 33.76dB) and the algorithm presented when the penalty function is used as the correction method (30.23 dB). The time of the proposed algorithm is 190.55 s, which is faster than that of the BOMP algorithm (302.14 s) and that when the penalty function is used as the correction method (586.15 s). When the sampling rate is 50%, the SNR of the proposed algorithm is 35.42 dB, which is better than that of the BOMP algorithm (34.56 dB) and that when the penalty function is used as the correction method (31.38 dB). The computation time of the proposed algorithm is 196.67 s, which is faster than that of the BOMP algorithm (617.62 s) and that when the penalty function is used as the correction method (1 071.15 s). When the sampling rate is 60%, the SNR of the proposed algorithm is 36.33 dB, which is better than that of the BOMP algorithm (35.18 dB) and that when the penalty function is used as the correction method (33.57 dB). The computation time of the proposed algorithm is 201.72 s, which is faster than that of the BOMP algorithm (1 136.29 s) and that when the penalty function is used as the correction method (1 505.35 s). At a sampling rate of 70%, the SNR of the proposed algorithm is 38.62 dB, which is better than that of the BOMP algorithm (37.65 dB) and that when the penalty function is used as the correction method (35.17 dB). The calculation time of the proposed algorithm is 214.68 s, which is faster than that of the BOMP algorithm (1 802.42 s) and that when the penalty function is used as the correction method (2 415.81 s). Experimental results show that the improved algorithm is superior to the BOMP algorithm and algorithm when the penalty function is used as the correction method in terms of reconstruction precision and calculation time under the same sampling rate. When the sampling rate is higher, the precision of the reconstructed image will be higher and reconstruction time will be shorter. Conclusion The convex optimization algorithm for image reconstruction based on compression perception theory is superior to the non-convex optimization algorithm for image reconstruction.However, reconstruction precision emains low and reconstruction rate emains slow, and thus, the problem remains serious. This study proposes an image reconstruction algorithm based on Lp norm compression awareness. Experiment results indicate that the proposed algorithm generally improves image reconstruction precision and computation time. Further study should be conducted to establish an algorithm with higher precision and faster image reconstruction.

# Key words

compressive sensing; image reconstruction; norm; sampling; ratevalue function; sequence quadratic programming

# 1.1 Lp范数的信号重构理论框架

 $\mathop {{\rm{min}}}\limits_{x \in {\boldsymbol{\rm{R}}^{N \times 1}}} ||\mathit{\boldsymbol{ \boldsymbol{\varPsi} x|}}{\mathit{\boldsymbol{|}}_0}\;\;\;{\rm{s}}{\rm{.t}}{\rm{. }}\;\;y = \mathit{\boldsymbol{ \boldsymbol{\varPhi} \boldsymbol{\varPsi} x}}$ (1)

L1范数优化依然存在计算复杂度高、需要大量观测样本并且传输数据间有较大冗余性等缺点。进一步改进重构算法性能，直接将这个问题的求解转化成sp的Lp(0 < p < 1) 范数最优化问题，即

 $\mathop {{\rm{min}}}\limits_{x \in {\boldsymbol{\rm{R}}^{N \times 1}}} ||s|{|_p}\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;y = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}s$ (2)

# 1.2 基于拉格朗日函数Hesse矩阵和拉格朗日函数序列二次规划 (SQP) 方法的Lp范数优化算法

 $f\left( s \right)\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;c\left( s \right) = y-\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}s$ (3)

 $L\left( {s, \mathit{\pmb{\lambda}}} \right) = f\left( s \right)-{\mathit{\pmb{\lambda}} ^{\rm T}}{c_E}(s)$ (4)

 ${\mathit{\boldsymbol{J}}_E}\left( s \right) = \nabla {c_E}\left( s \right) = {[\nabla {c_1}(s)\;\;\;\;\nabla {c_2}\left( s \right)\;\;\; \cdots \;\;\;\nabla {c_{m1}}\left( s \right)]^{\rm T}}$ (5)

 $\nabla _{ss}^2L\left( {s, \mathit{\pmb{\lambda}}} \right) = \mathit{\boldsymbol{H}}\left( \mathit{\boldsymbol{s}} \right)-\sum\limits_{i = 1}^{{m_i}} {{\mathit{\pmb{\lambda}}_i}} {\nabla ^2}{c_i}(s)$ (6)

 $\nabla _{_{s\lambda }}^{^2}L\left( {s, \mathit{\pmb{\lambda}}} \right) =-{\mathit{\boldsymbol{J}}_E}{\left( x \right)^{\rm{T}}}$ (7)

 $\nabla _{_{\lambda s}}^{^2}L\left( {s, \mathit{\pmb{\lambda}}} \right) =-{\mathit{\boldsymbol{J}}_E}\left( x \right)$ (8)

 $\nabla _{_{\lambda \lambda }}^{^2}L\left( {s, \mathit{\pmb{\lambda}}} \right) = 0$ (9)

 ${\mathit{\boldsymbol{Z}}_k} = {\mathit{\boldsymbol{\theta }}_k}{\mathit{\boldsymbol{y}}_k} + \left( {1-{\mathit{\boldsymbol{\theta }}_k}} \right){\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{p}}_k}$ (19)

 ${\mathit{\boldsymbol{B}}_{k + 1}} = \frac{{{\mathit{\boldsymbol{B}}_k}-{\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{p}}_k}\mathit{\boldsymbol{p}}_{_k}^{\rm{T}}{\mathit{\boldsymbol{B}}_k}}}{{\mathit{\boldsymbol{p}}_{_k}^{\rm{T}}{\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{p}}_k}}} + \frac{{{\mathit{\boldsymbol{z}}_k}\mathit{\boldsymbol{z}}_{_k}^{\rm{T}}}}{{\mathit{\boldsymbol{p}}_{_k}^{\rm{T}}{\mathit{\boldsymbol{z}}_k}}}$ (20)

 ${\theta _k} = \left\{ \begin{array}{l} 1\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{p}}_{_k}^{\rm{T}}{\mathit{\boldsymbol{y}}_k} \ge 0.2\mathit{\boldsymbol{p}}_{_k}^{\rm{T}}{\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{p}}_k}\\ \frac{{0.8\mathit{\boldsymbol{p}}_{_k}^{\rm{T}}{\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{p}}_k}}}{{\mathit{\boldsymbol{p}}_{_k}^{\rm{T}}{\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{p}}_k}-\mathit{\boldsymbol{p}}_{_k}^{\rm{T}}{\mathit{\boldsymbol{y}}_k}}}\;\;\mathit{\boldsymbol{p}}_{_k}^T{\mathit{\boldsymbol{y}}_k} < 0.2\mathit{\boldsymbol{p}}_{_k}^{\rm{T}}{\mathit{\boldsymbol{B}}_k}{\mathit{\boldsymbol{p}}_k} \end{array} \right.$ (21)

 $\phi F\left( {s, \sigma } \right) = f\left( s \right)-\mathit{\boldsymbol{\mu }}{\left( s \right)^{\rm{T}}}{\boldsymbol{c}}\left( s \right) + \frac{1}{2}||\mathit{\boldsymbol{\sigma }}{\boldsymbol{c}}\left( s \right)|{|^2}$ (22)

 $\mu \left( s \right) = {[{\mathit{\boldsymbol{J}}_E}\left( \mathit{\boldsymbol{s}} \right){\mathit{\boldsymbol{J}}_E}{\left( s \right)^{\rm{T}}}]^{ -1}}{\mathit{\boldsymbol{J}}_E}(\mathit{\boldsymbol{s}})\nabla f\left( s \right)$ (23)

μ(s) 是式 (24) 最小二乘问题的解。

 $\mathop {{\rm{min}}}\limits_{\mu \in {\Re ^i}} ||\nabla f\left( s \right)-{\mathit{\boldsymbol{J}}_E}{\left( \mathit{\boldsymbol{s}} \right)^{\rm{T}}}\mu ||$ (24)

Fletcher价值函数在${\mathit{\boldsymbol{s}}_k}$处的梯度为

 $\begin{array}{l} \nabla \phi F\left( {{s_k}, \sigma } \right) = \nabla f\left( {{s_k}} \right)-\mathit{\boldsymbol{J}}_{_{{E_k}}}^{^{\rm{T}}}\mu \left( {{s_k}} \right)-\\ \nabla \mu \left( {{s_k}} \right)\mathit{\boldsymbol{c}}\left( {{s_k}} \right) + {\sigma ^{-1}}\mathit{\boldsymbol{J}}_{_{{E_k}}}^T\mathit{\boldsymbol{c}}({s_k}) \end{array}$ (25)

${\mathit{\boldsymbol{B}}_k}$替代迭代点处的Hesse矩阵所得的第$k$次迭代过程中搜索方向${d_k}$

 ${d_k} = \mathit{\boldsymbol{J}}_{_{{E_k}}}^{\rm{T}}d_{_k}^{^y} + {z_k}d_{_k}^{^z}$ (26)

 $d_{_k}^{^y} =-{({\mathit{\boldsymbol{J}}_{{E_k}}}\mathit{\boldsymbol{J}}_{_{{E_k}}}^{^{\rm{T}}})^{-1}}\mathit{\boldsymbol{c}}({s_k})$ (27)

 $\nabla f{\left( {{s_k}} \right)^{\rm{T}}}\mathit{\boldsymbol{J}}_{_{{E_k}}}^{\rm{T}}d_{_k}^{^y} =-\mu {\left( {{s_k}} \right)^{\rm{T}}}\mathit{\boldsymbol{c}}({s_k})$ (28)

 $\begin{array}{l} \nabla \phi \mathit{\boldsymbol{F}}{\left( {{s_k}, \sigma } \right)^{\rm{T}}}{d_k} = \nabla f{\left( {{s_k}} \right)^{\rm{T}}}{\boldsymbol{z}_k}\boldsymbol{d}_{_k}^{^z}-\\ \mathit{\boldsymbol{\mu }}{\left( {{s_k}} \right)^{\rm{T}}}\nabla \mathit{\boldsymbol{\mu }}{\left( {{s_k}} \right)^{\rm{T}}}{\mathit{\boldsymbol{d}}_k}-{\sigma ^{-1}}||\mathit{\boldsymbol{c}}\left( {{s_k}} \right)|{|^2} \end{array}$ (29)

 ${\mathit{\boldsymbol{W}}_k}{\mathit{\boldsymbol{d}}_k} + \nabla f\left( {{s_k}} \right)-{\mathit{\boldsymbol{J}}_{{E_k}}}{\mu _k} = 0$ (30)

 ${z^{\rm{T}}}{\mathit{\boldsymbol{W}}_k}{\mathit{\boldsymbol{d}}_k} =-z_{_k}^{\rm{T}}\nabla f\left( {{s_k}} \right)$ (31)

 $\begin{array}{l} \nabla \phi \mathit{\boldsymbol{F}}{\left( {{s_k}, \sigma } \right)^{\rm{T}}}{\mathit{\boldsymbol{d}}_k} =- {\left( {\mathit{\boldsymbol{d}}_{_k}^{^z}} \right)^{\rm{T}}}[z_{_k}^{\rm{T}}{\mathit{\boldsymbol{W}}_k}{\mathit{\boldsymbol{z}}_k}]\mathit{\boldsymbol{d}}_{_k}^{^z} -\\ {\left( {\mathit{\boldsymbol{d}}_{_k}^{^y}} \right)^{\rm{T}}}{\mathit{\boldsymbol{J}}_{{E_k}}}{\mathit{\boldsymbol{W}}_k}{\mathit{\boldsymbol{z}}_k}\mathit{\boldsymbol{d}}_{_k}^{^z} -\mathit{\boldsymbol{c}}{\left( {{s_k}} \right)^{\rm{T}}}\nabla \mu {\left( {{s_k}} \right)^{\rm{T}}}{\mathit{\boldsymbol{d}}_k} -{\sigma ^{ - 1}}||\mathit{\boldsymbol{c}}\left( {{s_k}} \right)|{|^2} \end{array}$ (32)

 ${\sigma ^{- 1}} \ge \left[{\frac{{{{\left( {\mathit{\boldsymbol{d}}_{_k}^{^y}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{J}}_{{E_k}}}{\mathit{\boldsymbol{W}}_k}{\mathit{\boldsymbol{z}}_k}\mathit{\boldsymbol{d}}_{_k}^{^z}-\mathit{\boldsymbol{c}}{{\left( {{s_k}} \right)}^{\rm{T}}}\nabla \mu {{\left( {{s_k}} \right)}^{\rm{T}}}{\mathit{\boldsymbol{d}}_k}}}{{||\mathit{\boldsymbol{c}}\left( {{s_k}} \right)|{|^2}}}} \right] + \bar \sigma$ (33)

1) 选择初始稀疏${s_0} \in {\boldsymbol{\rm{R}}^n}$，初始对称正定矩阵B0为1的单位矩阵, 初始价值参数σ0, 计算${z_0}$, E>0, 令K=0;

2) 首先将大小为N×N的原始图像进行分块，对每一小块图像用变换矩阵ψ进行稀疏变换，使用手段是小波变换，再选择高斯随机矩阵Φ作为测量矩阵，从而计算出该测量矩阵下的测量值$y$;

3) 若$|z_{_k}^{\rm{T}}f\left( {{s_k}} \right)|| + ||\mathit{\boldsymbol{c}}\left( {{s_k}} \right)|| \le E$，停止，解得${\mathit{\boldsymbol{s}}_k}$；否则，转步骤4);

4) 计算式 (26) 得到迭代点处的搜索方向${\mathit{\boldsymbol{d}}_k}$，并修正乘子μ(s), 在${\mathit{\boldsymbol{s}}_k} + {\mathit{\boldsymbol{d}}_k}$处用最小二乘法估计式 (26)，计算$\mu ({\mathit{\boldsymbol{s}}_{k + 1}})$以及$\Delta \mu \left( {{\mathit{\boldsymbol{s}}_k}} \right) = \mu \left( {{\mathit{\boldsymbol{s}}_{k + 1}}} \right)-\mu ({\mathit{\boldsymbol{s}}_k})$;

5) 通过式 (33) 估计价值参数${\sigma _k}$, 修正Hesse矩阵求解式 (20) 中${\mathit{\boldsymbol{B}}_{k + 1}}$;

6) 令K:=K+1，${\mathit{\boldsymbol{s}}_{k + 1}} = {\mathit{\boldsymbol{s}}_k} + \mathit{\boldsymbol{\mu }}\left( s \right){\mathit{\boldsymbol{d}}_k}$, 返回步骤2);

7) 采用步骤1) 中相同的稀疏矩阵ψ对所得的稀疏系数s进行逆变换，获得每一块分块图像的重构图像，然后按照顺序依次重组获取原始重构图像。

# 2 图像重构算法结果及分析

3种算法都采用图像分块技术，从图 3中可以看出与未分块算法图像重构质量有明显地提高，但是重构图像边界处清晰度不足。将图像分成尺度相同的块，而每个图像块具有不同图像信息，稀疏度小的块，重构效果更佳。因此，使用图像分块技术导致图像边界不连续，产生块效应。

Table 1 Performance comparison of three reconstruction algorithms of image Lena

 采样率/% 信噪比/dB 计算时间/s 本文算法 BOMP算法 文献[9] 本文算法 BOMP算法 文献[9] 40 34.28 33.60 30.23 190.55 302.14 586.15 50 35.42 34.56 31.38 196.67 617.62 1 071.15 60 36.33 35.18 33.57 201.72 1 136.29 1 505.35 70 38.62 37.65 35.17 214.68 1 802.42 2 415.81

Table 2 Performance comparison of three reconstruction algorithms of image peppers

 采样率/% 信噪比/dB 计算时间/s 本文算法 BOMP算法 文献[9] 本文算法 BOMP算法 文献[9] 40 33.25 32.83 30.23 186.52 282.87 614.92 50 34.53 33.72 31.26 212.35 567.23 968.41 60 35.37 34.62 33.43 218.46 981.74 1 602.04 70 37.46 36.35 35.67 220.67 1 729.65 2 315.26

# 参考文献

• [1] Ren Y M, Zhang Y N, Li Y. Progress and prospect of compressed sensing and its application in image processing[J]. Journal of Automation, 2014, 40(8): 1563–1575. [任越美, 张艳宁, 李映. 压缩感知及其图像处理应用研究进展与展望[J]. 自动化学报, 2014, 40(8): 1563–1575. ] [DOI:10.3724/SP.J.1004.2014.01563]
• [2] Donoho D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289–1306. [DOI:10.1109/TIT.2006.871582]
• [3] Candès E J, Romberg J, Tao T. Robust uncertainty principles:exact signal reconstruction from highly incomplete frequency information[J]. IEEE Transactions on Information Theory, 2006, 52(2): 489–509. [DOI:10.1109/TIT.2005.862083]
• [4] Tropp J A, Gilbert A C. Signal recovery from random measurements via orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655–4666. [DOI:10.1109/TIT.2007.909108]
• [5] Wu R, Huang W, Chen D R. The exact support recovery of sparse signals with noise via orthogonal matching pursuit[J]. IEEE Signal Processing Letters, 2013, 20(4): 403–406. [DOI:10.1109/LSP.2012.2233734]
• [6] Donoho D L, Tsaig Y, Drori I, et al. Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2012, 58(2): 1094–1121. [DOI:10.1109/TIT.2011.2173241]
• [7] Li S D, Pei W J, Hu G Q, et al. OMP reconstruction algorithm via Bayesian model and its application[J]. Systems Engineering and Electronics, 2015, 37(2): 246–252. [李少东, 裴文炯, 胡国旗, 等. 贝叶斯模型下的OMP重构算法及应用[J]. 系统工程与电子技术, 2015, 37(2): 246–252. ] [DOI:10.3969/j.issn.1001-506X.2015.02.03]
• [8] Ning F L, He B J, Wei J. An algorithm for image reconstruction based on lp norm[J]. Acta Physica Sinica, 2013, 62(17): 282–289. [宁方立, 何碧静, 韦娟. 基于lp范数的压缩感知图像重建算法研究[J]. 物理学报, 2013, 62(17): 282–289. ] [DOI:10.7498/aps.62.174212]
• [9] Ma C, Liu Z, Zhen X X. Combining penalty function with sequential quadratic programming method for lp norm minimization[J]. Computer Engineering and Applications, 2013, 49(18): 212–216. [马聪, 刘哲, 甄小仙. 结合罚函数与序列二次规划的lp范数优化方法[J]. 计算机工程与应用, 2013, 49(18): 212–216. ] [DOI:10.3778/j.issn.1002-8331.1112-0516]
• [10] Chartrand R, Yin W. Iteratively reweighted algorithms for compressive sensing[C]//Proceedings of 2008 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). Las Vegas, NV, USA:IEEE, 2008:3869-3872.[DOI:10.1109/ICASSP.2008.4518498]
• [11] Li Z C, Tang J H. Unsupervised feature selection via nonnegative spectral analysis and redundancy control[J]. IEEE Transactions on Image Processing, 2015, 24(12): 5343–5355. [DOI:10.1109/TIP.2015.2479560]