传统图像获取技术的缺陷随着科技高速发展日益凸显，奈奎斯特釆样频率已经成为海量信息实时釆样和处理的桎梏，多数以巨大代价采样得到的数据最终在压缩过程中会被摒弃^[1]。由此采样得到的海量数据更成为了信号存储与传输的沉重负担，这给采样数据和传输的硬件设备带来了巨大挑战。

近年来由国内外知名学者Donoho^[2]和Candes等人^[3]正式提出了压缩感知理论，该理论指出一个可以压缩的信号或者在某个变换域是稀疏的信号，可以用一个与变换基毫不相关的测量矩阵将变换所得到的高维信号投影到一个低维空间上，然后通过求解一个类似欠定方程组的问题，选择优化算法从这些少量的投影中以高概率重构出来。

压缩感知主要包括三大核心内容：信号的稀疏表示，构造测量矩阵以及重构算法。而压缩感知理论中最关键的技术则是重构算法。压缩感知从数学意义上来讲就是求解最小L₀范数最优化问题，但是该问题属于NP难问题，需要大量的计算，因此把此问题转换为L₁范数予以解决，把非凸优化问题转换为凸优化问题。为了求解这一问题答案，科学家们近年来提出了许多优秀算法，常用算法包括，匹配追踪算法 (MP)、正交匹配追踪算法^[4](OMP) 以及类似的贪婪类算法^[5-7]。贪婪类算法拥有非常好的重建速度，算法计算很简单，但是在大多数情况下需要知道待重构信号的稀疏度，这是在工程实践当中很难得到的，并且精确重建理论保证性较弱，重建精度低。

从另一个方面探讨，科学研究中提出了可以将L₀范数转化为L₁范数问题求解，将非凸优化问题转换为凸优化问题，进而得到相似的解。近些年出现了一些比较好的凸优化算法，包括基追踪算法、匹配追踪算法、梯度投影算法等，但是这类算法在运行中需要大量的观测数据，因此，不适用于大规模实际应用问题。针对这类问题，文献[8-11]提出了基于L_p(0 < p<1) 范数实现算法，实验表明了在重构信号效果和可靠性方面，L_p( < p < 1) 范数明显强于L₀范数和L₁范数。

压缩感知理论是在远小于Nyquist采样率的条件下，用随机采样获取部分离散信号。由于压缩感知本质是非适应、非线性的可变换域内的信号重构算法，因此可以通过非线性重构算法重建出完美的重建信号。

压缩感知中信号重建问题，其实就是在欠定方程组中求出最优解的问题，本质是求解最小L₀范数的NP难问题。构造求解L₀问题模型为

$ \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)

式中，Φ为$ m \times m $测量矩阵，ψ为$ n \times n $稀疏基矩阵。

L₁范数优化依然存在计算复杂度高、需要大量观测样本并且传输数据间有较大冗余性等缺点。进一步改进重构算法性能，直接将这个问题的求解转化成s_p的L_p(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)

式中，0 < p < 1，$ f\left( s \right) = ||s|{|_p} = {\left( {\sum\limits_{i = 1}^\infty {|{s_i}{|^{^p}}} } \right)^{1/p}}$是目标函数。

式 (2) 的目的是在图 1中所示直线G(满足$ y = \mathit{\boldsymbol{ \boldsymbol{\varPhi} x}} $) 上找到一个点使得L_p球的半径最小，当0 < p < 1时，L_p球是内凸的，当球的半径不断增大时与G将相交于Y坐标轴上，这样的一个点是稀疏的。菱形L₁球是p取值的一个特殊情况，在一定的条件下也会有一个稀疏。但是当p>1时，L_p球外凸，当球逐渐膨胀时与图像切点不交于Y坐标轴上，此解是不稀疏的。L₂范数作为p>1时的一个特例，但也不能得到稀疏解。因为L₂-norm的最小化为$ x_1^2 + x_2^2 $, 因此找的时候是从原点开始逐渐增大半径$ r $，直到与$Ax = b $的直线相切为止过几何图形。而L_p-norm的最小化为$ x_1^{^p} + x_2^{^p} $，凹凸有致，根据p的取值不同区分球体是外凸还是内凹。如图 1所示。

序列二次规划方法是近几年来公认的用来求解约束非线性规划问题的有效方法之一。与其他算法相比，它最突出的优点是收敛性好、计算效率高、边界搜索能力强。其基本思想是：通过每一步迭代求解一个二次规划子问题来确定一个下降方向，从而减少价值函数步长，重复此步骤直到求出原问题的最优解。

令式 (2) 转化为

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

(3)

式中，$ f(s) = \mathop {{\rm{min}}}\limits_{x \in {{{\boldsymbol{\rm{R}}}}^{N \times 1}}} ||s|{|_{\rm{p}}}$，Φ为$ m \times n $测量矩阵，式 (3) 的拉格朗日函数为

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

(4)

式中, $ \mathit{\pmb{\lambda}} = {({\lambda _1}, {\lambda _2}, \cdots {\lambda _n})^{\rm T}} $为拉格朗日乘子向量，$ {c_E}\left( s \right) $的Jacobi矩阵记为

$ {\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)

则拉格朗日函数关于s和λ的二阶梯度向量信息，得

$ \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{H}}\left( s \right) \in {\boldsymbol{\rm{R}}^{n \times n}} $为$f(s) $的海瑟 (Hesse) 矩阵。结合式 (6)—式 (9) 可得

$ \begin{array}{l} {\nabla ^2}L\left( {s, \mathit{\pmb{\lambda}}} \right){\rm{ }} = \left[{\begin{array}{*{20}{c}} {\nabla _{ss}^2L\left( {s, \mathit{\pmb{\lambda}}} \right)} & {\nabla _{_{x\lambda }}^{^2}L(s, \mathit{\pmb{\lambda}})}\\ {\nabla _{_{\lambda s}}^{^2}L\left( {s, \mathit{\pmb{\lambda}}} \right)} & {\nabla _{_{\lambda \lambda }}^{^2}L\left( {s, \mathit{\pmb{\lambda}}} \right)} \end{array}} \right] = \\ \left[{\begin{array}{*{20}{c}} {H\left( s \right)-\sum\limits_{i = 1}^{{m_1}} {{\mathit{\pmb{\lambda}}_i}} {\nabla ^2}{c_i}\left( x \right)} & {-{\mathit{\boldsymbol{J}}_E}{{\left( s \right)}^{\rm T}}}\\ {-{\mathit{\boldsymbol{J}}_E}\left( s \right)} & {\bf{0}} \end{array}} \right] \end{array} $

(10)

再将拉格朗日函数关于s的Hesse矩阵$ \nabla _{ss}^2L\left( {s, \mathit{\pmb{\lambda}}} \right) $记为$ W(s, \mathit{\pmb{\lambda}}) $，简化为

$ \begin{array}{l} {\nabla ^2}L\left( {s, \mathit{\pmb{\lambda}}} \right) = \left[{\begin{array}{*{20}{c}} {\nabla _{ss}^2L\left( {s, \mathit{\pmb{\lambda}}} \right)} & {\nabla _{_{x\lambda }}^{^2}L(s, \mathit{\pmb{\lambda}})}\\ {\nabla _{_{\lambda s}}^{^2}L\left( {s, \mathit{\pmb{\lambda}}} \right)} & {\nabla _{_{\lambda \lambda }}^{^2}L\left( {s, \mathit{\pmb{\lambda}}} \right)} \end{array}} \right] = \\ \left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{W}}\left( {s, \mathit{\pmb{\lambda}}} \right)} & {-{\mathit{\boldsymbol{J}}_E}{{\left( s \right)}^{\rm T}}}\\ {-{\mathit{\boldsymbol{J}}_E}\left( s \right)} & {\bf{0}} \end{array}} \right] \end{array} $

(11)

对于给定点$ {\mathit{\boldsymbol{T}}_k} = ({s_{k, }}{\mathit{\pmb{\lambda}}_k}) $，牛顿法迭代格式为

$ {\mathit{\boldsymbol{T}}_{k + 1}} = {\mathit{\boldsymbol{T}}_k} + {\mathit{\boldsymbol{p}}_k} $

(12)

这里$ {\mathit{\boldsymbol{p}}_k} = ({s_{k, }}{r_k}) $利用牛顿拉普森公式，满足线性方程组

$ \left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{W}}\left( {s, \mathit{\pmb{\lambda}}} \right)} & {-{\mathit{\boldsymbol{J}}_E}{{\left( s \right)}^{\rm T}}}\\ {-{\mathit{\boldsymbol{J}}_E}\left( s \right)} & {\bf{0}} \end{array}} \right]\left[{\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{d}}_k}}\\ {{\mathit{\boldsymbol{r}}_k}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} {-\nabla f\left( {{s_k}} \right) + {\mathit{\boldsymbol{J}}_E}{{\left( s \right)}^{\rm{T}}}{\mathit{{\lambda}}_k}}\\ {\mathit{\boldsymbol{c}}(s)} \end{array}} \right] $

(13)

不难发现，只要矩阵$ \mathit{\boldsymbol{d}} \ne {\boldsymbol{0}} $满秩且$\mathit{\boldsymbol{W}}(s, \mathit{\pmb{\lambda}}) $是正定的，式 (13) 就是非奇异的系数矩阵，并且该方程组有且仅有唯一解。由于每一次迭代在求解方程组式 (13) 中所获得的数值不太稳定，故考虑将它转化为一个严格凸二次规划问题。转化条件是式 (3) 的解点s_*处最优性二阶充分条件成立，即对满足$ {\mathit{\boldsymbol{J}}_E}{({s^*})^{\rm{T}}}\mathit{\boldsymbol{d}} = \mathit{\boldsymbol{0}} $的任一向量$\mathit{\boldsymbol{d}} \ne {\boldsymbol{0}} $成立。

这时，当α>0且无限接近于0时，有

$ \mathit{\boldsymbol{W}}\left( {{s^*}, {\mathit{\pmb{\lambda}}^*}} \right) + \frac{1}{{2\alpha }}{\mathit{\boldsymbol{J}}_E}{\left( {{s^*}} \right)^{\rm{T}}}{\mathit{\boldsymbol{J}}_E}({s^*}) $

(14)

通过分析可得，式 (14) 正定，则可以用一个正定矩阵$ \mathit{\boldsymbol{B}}\left( {s, \mathit{\pmb{\lambda}}} \right){\rm{ }} = {\rm{ }}W\left( {s, \mathit{\pmb{\lambda}}} \right) + \frac{1}{{2\alpha }}\mathit{\boldsymbol{J}}_{_E}^{\rm{T}}{\mathit{\boldsymbol{J}}_E}(s) $来替代，记为

$ \mathit{\boldsymbol{B}}\left( {s, \mathit{\pmb{\lambda}}} \right){\rm{ }} = \mathit{\boldsymbol{W}}\left( {s, \mathit{\pmb{\lambda}}} \right) + \frac{1}{{2\alpha }}\mathit{\boldsymbol{J}}_{_E}^{\rm{T}}{\mathit{\boldsymbol{J}}_E}(s) $

(15)

当$ \left( {s, \mathit{\pmb{\lambda}}} \right) \to ({s^*}, {\mathit{\pmb{\lambda}}^*}) $时，矩阵$ \mathit{\boldsymbol{B}}({s^*}, {\mathit{\pmb{\lambda}}^*}) $正定。式 (15) 通过变形可得

$ \mathit{\boldsymbol{B}}\left( {s, \mathit{\pmb{\lambda}}} \right){\rm{ }}{d_k}-{\mathit{\boldsymbol{J}}_E}{\left( s \right)^{\rm{T}}}\mathit{\boldsymbol{\bar \lambda }} =-\nabla f\left( {{s_k}} \right) $

(16)

因此式 (13) 可以等价为

$ \left[{\begin{array}{*{20}{c}} {\mathit{\boldsymbol{B}}\left( {s, \mathit{{\lambda}}} \right)} & {-{\mathit{\boldsymbol{J}}_E}{{\left( s \right)}^{\rm T}}}\\ {-{\mathit{\boldsymbol{J}}_E}\left( s \right)} & 0 \end{array}} \right]{\rm{ }}\left[{\begin{array}{*{20}{c}} {{\boldsymbol{d}_k}}\\ {{{\mathit{{\bar \lambda }}}_k}} \end{array}} \right] = - \left[{\begin{array}{*{20}{c}} {-\nabla f\left( {{s_k}} \right)}\\ {\mathit{\boldsymbol{c}}(s)} \end{array}} \right] $

(17)

式中，$ {{\bar \lambda }_k} = {\lambda _k} + {r_k} + \frac{1}{{2\alpha }}{\mathit{\boldsymbol{J}}_E}{\left( s \right)^{\rm{T}}}{\mathit{\boldsymbol{J}}_E}\left( s \right) $。之后，可以把方程组 (17) 转化为严格凸二次规划问题来求解。

在拉格朗日函数Hesse矩阵SQP方法中，因为在计算拉格朗日函数在迭代点$ {s_k} $处的Hesse矩阵$ W(s, \mathit{{\lambda}}) $时计算量巨大，为了弥补这一缺陷，本文采用修正Hesse矩阵SQP的方法。利用对称正定矩阵${\mathit{\boldsymbol{B}}_\mathit{\boldsymbol{k}}} $代替迭代点处的Hesse矩阵，可以在很大程度上减少计算量。

修正拉格朗日函数$ L\left( {s, \mathit{\pmb{\lambda}}} \right) $的Hesse矩阵，用对称正定矩阵${\mathit{\boldsymbol{B}}_\mathit{\boldsymbol{k}}} $替代Hesse矩阵。令

$ \begin{array}{l} {\mathit{\boldsymbol{p}}_k} = {\mathit{\boldsymbol{s}}_{k + 1}}-{\mathit{\boldsymbol{s}}_k}\\ {\mathit{\boldsymbol{y}}_k} = {\nabla _s}L\left( {{\mathit{\boldsymbol{s}}_{k + 1}}, {\mathit{\pmb{\lambda}}_{k + 1}}} \right)-{\nabla _s}L({\mathit{\boldsymbol{s}}_{k + 1}}, {\mathit{\pmb{\lambda}}_{k + 1}}) \end{array} $

(18)

式中，向量$ {\mathit{\boldsymbol{s}}_k} $和$ {\mathit{\boldsymbol{y}}_k} $必须满足$ \mathit{\boldsymbol{p}}_{_k}^{^{\rm{T}}}{\mathit{\boldsymbol{y}}_k} $>0的条件，但是式 (18) 确定的向量$ {\mathit{\boldsymbol{s}}_k} $和$ {\mathit{\boldsymbol{y}}_k} $不一定满足这一条件，为此，可以对向量$ {\mathit{\boldsymbol{y}}_k} $进行修正。令

$ {\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} $矩阵为

$ {\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)

由式 (18)—式 (21) 可以保持$ {\mathit{\boldsymbol{B}}_k} $正定。

在迭代过程中，为了保证SQP方法的全局收敛特性，常常借助于价值函数来确定搜索步长，$ {d_k} $作为$ {\mathit{\boldsymbol{B}}_k} $代替迭代点处的Hesse矩阵所得的第$k $次迭代过程中的搜索方向。这个搜索方向$ {d_k} $拥有一个很好的性质，它是很多价值函数的下降方向。引入Fletcher价值函数便得到算法的全局收敛。Fletcher价值函数为

$ \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)

式中，μ(s) 为乘子向量，σ>0为价值参数。若函数c(s) 的Jacobi矩阵$ {\mathit{\boldsymbol{J}}_E}\left( \mathit{\boldsymbol{s}} \right) = \nabla \mathit{\boldsymbol{c}}{(\mathit{\boldsymbol{s}})^{\rm{T}}} $是行满秩的，则乘子向量可取为

$ \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)

式中, $ {\mathit{\boldsymbol{z}}_k} $的列向量为零空间$ N({J_{{E_k}}}) $的一组基。利用式 (3) 可行性条件有

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

(27)

然后，利用最小二乘子$ \mu ({s_k}) $的式 (26) 可得

$ \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)

另外, 由KT条件有

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

(30)

式中，$ {\mu _k} $为拉格朗日乘子, 那么有

$ {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)

因此，当Hesse矩阵正定时，若取价值参数满足

$ {\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)

式中，$ \bar \sigma > 0 $，则$ {d_k} $是Fletcher价值函数的下降方向。

本文算法具体步骤如下：

1) 选择初始稀疏$ {s_0} \in {\boldsymbol{\rm{R}}^n} $，初始对称正定矩阵B₀为1的单位矩阵, 初始价值参数σ₀, 计算${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进行逆变换，获得每一块分块图像的重构图像，然后按照顺序依次重组获取原始重构图像。

为了实验证明本文改进的压缩感知图像重构算法性能，在实验中选取尺度为256×256像素的boat图像，观测矩阵则采用高斯随机矩阵。运用本文改进算法进行图像重构实验，重构结果如图 2所示。

针对命名为boat的原始图像，当采样率分别是18.2%、19.7%时，重构出来的图像所对应的信噪比分别为67.352 1 dB、73.108 6 dB, 从图 2中可以看出当采样率到达20.5%，信噪比高达85.154 3 dB时，几乎完全达到了图像的精确重构。随着采样率不断地提高，图像重构精度也进一步提高。

依然采用boat图像作为初始图像，在采样率为30%时，分别采用本文算法、文献[7]提出的BOMP (block orthogonal matching pursuit) 算法和文献[9]中算法，分别没有采用分块技术及采用分块技术进行图像重构，图像重建结果如图 3所示。

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

分块后，3种算法在采样率为M/N=0.5时，图像重构效果进行比较，选择尺度为256像素×256像素的Lena、peppers两幅图像进行实验对比，图像重构结果如图 4所示。

测试图像Lena在采样率为50%时，本文重构算法、BOMP重构算法、文献[9]重构算法所对应的信噪比分别为35.42 dB, 34.56 dB, 31.38 dB。同样当测试图像为peppers，采样率为50%时，本文重构算法、BOMP重构算法、文献[9]重构算法信噪比分别为34.53 dB, 33.72 dB, 31.26 dB。无论对于Lena还是peppers图像来说，本文重构信噪比都高于BOMP重构算法和文献[9]重构算法。

表 1和表 2给出了3种算法在不同采样率下两幅重构图像客观结果，图像重构信噪比以及计算时间。

从表 1、表 2中可以看出，本文重构算法在采样率为40%至70%时的图像重构信噪比都大于BOMP重构算法和文献[9]中提到的重构算法的信噪比，并且在算法计算时间上，本文重构算法比后面两种算法，计算时间明显更短。相比本文算法，BOMP重构算法和文献[9]重构算法随着采样率的不断提高算法计算时间也明显增大，而本文重构算法计算时间比较平稳。因此，实验充分验证了本文重构算法的有效性和可行性。

在研究现有L_p范数的压缩感知图像重构算法基础上，针对SQP方法中海瑟矩阵不正定导致计算量很大的问题，引入价值函数，修正Hesse矩阵的序列二次规划方法并结合图像分块压缩感知技术，提出了L_p范数压缩感知图像重构算法。通过实验对本文算法、BOMP重构算法、文献[9]重构算法在信噪比和算法计算时间进行对比，在不同采样率下，本文算法都明显优于其他两种算法，而且在采样率仅为20.5%时, 信噪比高达85.154 3 dB, 重构图像比较清晰。本文算法的最大优点在于采用了分块压缩感知技术，提高图像重构效率，降低了重构时间，缺点是在图像采样率比较低的情况下，存在图像干扰块效应。接下来研究方向是如何在采样率低的情况下，高精度地还原图片，消除图像干扰块效应。本文算法加入了滤子技术，在一定程度上解决了初值预估计更加接近可行域，使得算法在图像重构效率得到提高。