CN102988026B - 基于乘子法的自发荧光断层成像重建方法 - Google Patents

Abstract

一种基于乘子法的自发荧光断层成像重建方法，使用有限元分析方法扩散方程进行离散化，基于L1范数的惩罚项建立无约束条件最优化问题模型；得到所述无约束条件最优化问题模型的对偶模型；建立所述对偶模型的增广拉格朗日函数；简化增广拉格朗日函数的最大值函数；使用截断牛顿法求解增广拉格朗日函数的最大值；将增广拉格朗日函数的梯度作为目标向量的最快下降方向对目标向量进行更新；更新惩罚向量；计算目标函数值J(w)，如果||(J(w)^k-J(w)^k-1)||/||Φ^m||≥tol为真，计算k＝k+1并跳至步骤S4，否则，结束计算，其中，tol为目标函数的收敛效率阈值。本发明能够快速地在较大成像区域内得到准确可靠的光源分布信息，除正则化参数的以外其他参数都可以实现自适应调整，从而提高了成像的鲁棒性。

Description

基于乘子法的自发荧光断层成像重建方法

技术领域

本发明属于光学分子影像领域，涉及自发荧光断层成像重建方法，尤其是一种基于乘子法的自发荧光断层成像重建方法。

背景技术

光学分子影像技术是在分子细胞水平上实现生物体生理、病理变化的实时、无创、在体成像。荧光断层成像作为一种重要的光学分子影像成像模态，具有高灵敏度、低费用、无放射性等诸多优点，已经在新型药物研发，在体观测肿瘤发生、发展、转移机理等多个领域获得了广泛的应用。由于在近红外和可见光波长范围里，光子在生物体组织中传播时会发生复杂的散射，因而二维平面自发荧光成像往往不能提供自发荧光光源的准确位置，而且具有无法获得荧光光源深度信息的固有缺陷。自发荧光断层成像重建技术通过测量边界表面的荧光光强分布，实现了荧光光源的三维精确定位。

自发荧光断层成像重建是一个典型的病态问题。而且通常，表面荧光光强非常微弱，导致CCD探测器所采集到的信号混有大量干扰噪声，从而进一步增加了荧光光源重建的复杂性，限制了自发荧光成像技术在实际中的应用。研发一种能够准确、快速、鲁棒的荧光光源重建技术一直是国际上一个研究热点，但是到目前为止，还没有一种公认的技术能够实现这一目标。

现有的重建方法大都是将荧光光源重建看作是求解最优化问题。正则化方法是一个很通用的方法。其中，最常见的方法是使用基于L2范数的正则化方法，但是这类方法得到的重建结果往往过于平滑，精度也不够理想。考虑到在实际应用中，如肿瘤的早期检测和观察，荧光光源分布具有稀疏性分布的特点，因此，将稀疏性特点作为先验信息融入到优化问题模型中可以增强荧光光源重建的质量。因此基于L1范数的正则化方法能够很好地体现光源的稀疏性特点，但这会造成最优化问题目标函数不可微分，这就使得现有的很多用于L2范数正则化方法的光源重建方法不能直接用来求解基于L1范数的光源重建问题。除此之外，尽管现在已经存在一些可以被用来求解基于L1范数的正则化问题的数学方法，如迭代收缩方法，内点法等等，但是，自发荧光断层成像问题有病态性强，以及系统矩阵为实数稠密矩阵的特点，这使得以上的方法存在成像效率低，参数选择困难，甚至只能得到局部最优解而不能够得到准确的重建结果。

发明内容

本发明的目的是提供一种基于乘子法的自发荧光断层成像重建方法。

为实现上述目的，一种基于乘子法的自发荧光断层成像重建方法，包括步骤：

S1使用有限元分析方法扩散方程进行离散化，基于L1范数的惩罚项建立无约束条件最优化问题模型；

S2得到所述无约束条件最优化问题模型的对偶模型；

S3建立所述对偶模型的增广拉格朗日函数；

S4：简化增广拉格朗日函数的最大值函数；

S5使用截断牛顿法求解增广拉格朗日函数的最大值；

S6将增广拉格朗日函数的梯度作为目标向量的最快下降方向对目标向量进行更新；

S7更新惩罚向量；

S8计算目标函数值J(w)，如果||(J(w)^k-J(w)^k-1)||/||Φ^m||≥tol为真，计算k＝k+1并跳至步骤S4，否则，结束计算，其中，tol为目标函数的收敛效率阈值。

本发明能够快速地在较大成像区域内得到准确可靠的光源分布信息。并且，本发明除正则化参数的以外其他参数都可以实现自适应调整，从而提高了成像的鲁棒性。

附图说明

图1，本发明方法的流程图。

图2，小鼠躯干的非匀质模型示意图以及对应的四面体网格。

图3，非匀质模型体内光源位置示意图以及表面测量值。

图4，两种方法的成像结果示意图：(a)(b)(c)依次为Tikhonov方法(全域成像)、Tikhonov方法(局部成像)以及本发明所得到的成像结果。左侧一列为使用iso-surface的三维视图，右侧一列为在最大值点处进行切片得到的荧光光强分布图。

图5，两个光源成像结果示意图：左侧为使用iso-surface的三维视图，右侧为在最大值点处进行切片得到的荧光光强分布图。

图6，三个光源成像结果示意图：左侧为使用iso-surface的三维视图，右侧为在最大值点处进行切片得到的荧光光强分布图。

具体实施方式

本发明所提出的方法充分考虑了自发荧光断层成像病态性强，系统矩阵为稠密矩阵的特点，在构建和求解拉格朗日乘子函数的过程中运用了对偶技术、自适应的迭代收缩策略、以及预算子共轭梯度法等技术。因此，本发明不但能够实现基于L1范数的光源重建，更重要的是降低了对参数选择的敏感度和提高了重建的效率。

下面将结合附图详细描述本发明的重建方法，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

图1是本发明的重建方法的总体流程。

下面对本发明的重建方法所涉及的关键步骤进行逐一详细说明，具体形式如下面所述：

步骤101：使用有限元分析方法将扩散方程进行离散化，结合成像区域的四面体网格剖分，得到光源分布和边界测量数据之间的线性关系，然后使用基于L1范数的惩罚项建立起无约束条件最优化问题模型：

$\min J\left(w\right)=\frac{1}{2}{\left|\right|\mathrm{Aw}-{\mathrm{\Φ}}^m\left|\right|}_2^2+\mathrm{\λ}{\left|\right|w\left|\right|}_1^1$

其中A是系统矩阵，Φ^m为通过CCD相机所得到的被观测物体表面的光强值，λ是正则化参数，J(w)为最小二乘最优化问题的目标函数，w表示光源的分布信息，其初始值设定为零；

尽管辐射传输方程能够对光子在生物组织中的传输进行精确地描述，但是该方程是一个包含多个变量的微分-积分方程，直接对该方程进行求解非常复杂。巨大的计算量限制了该方程在实际情况中的应用，所以在实际应用大都是使用该方程的简化和近似模型。由于在近红外和可见光波长范围内，大多数生物体组织都具有高散射体吸收的特性，即生物体对光子的散射作用远远大于吸收作用，所以可以使用扩散方程代替辐射传输方程来描述光子在生物体组织中的传输：

$-\▿\·(D(x,\mathrm{\λ})\▿\mathrm{\Φ}(x,\mathrm{\λ}))+{\mathrm{\μ}}_a(x.\mathrm{\λ})\mathrm{\Φ}(x,\mathrm{\λ})=S(x,\mathrm{\λ})(x\∈\mathrm{\Ω})$

其中Ω表示整个成像区域，Φ(r)表光流密度，S(r)表示光源强度，μ_a(r)是组织吸收系数，μ’_s(r)表示约化组织反射系数，D(r)＝1/3(μ_a(r)+μ’_s(r))表示光学扩算系数。

在实际应用过程中，边界上光学信号的采集工作都是在暗箱中进行的，所以假定没有光子通过边界进入成像区域Ω，这符合使用Robin边界条件的要求：

$\mathrm{\Φ}(x,\mathrm{\λ})+2A(x;n,n^{\′})D(x,\mathrm{\λ})(v\left(x\right)\·\▿\mathrm{\Φ}(x,\mathrm{\λ}))=0(x\∈\∂\mathrm{\Ω})$

结合有限元分析方法，我们可以得到扩散方程及其边界条件的弱解形式：

${\∫}_{\mathrm{\Ω}}(D(x,\mathrm{\λ})(\▿\mathrm{\Φ}(x,\mathrm{\λ}))(\▿\mathrm{\ψ}(x,\mathrm{\λ}))+{\mathrm{\μ}}_a(x,\mathrm{\λ})\mathrm{\Φ}(x,\mathrm{\λ})\mathrm{\ψ}(x,\mathrm{\λ}))\mathrm{dx}$

$+{\∫}_{\∂\mathrm{\Ω}}\frac{1}{2A_n\left(x\right)}\mathrm{\Φ}(x,\mathrm{\λ})\mathrm{\ψ}(x,\mathrm{\λ})\mathrm{dx}={\∫}_{\mathrm{\Ω}}S(x,\mathrm{\λ})\mathrm{\ψ}(x,\mathrm{\λ})\mathrm{dx}$

其中，Ψ(x，λ)为基函数。

根据对成像区域进行四面体剖分得到的网格，结合对应的基函数，可以将上述的弱解形式离散化，得到下面的等式：

PΦ＝FS

其中，P是正定矩阵，Φ表示边界结点上光强，F表示光源加权矩阵，S表示成像区域内的光强分布。

由于Φ中包含不可测量结点，所以去除掉P^-1F中与之对应的行向量，并得到：

AS＝Φ^m

其中A∈R^M×N表示系统矩阵，Φ^m表示测量值。

由于荧光断层成像是一个有严重病态性的逆问题，现有技术都是需要融入一个惩罚项函数作为正则化项使得求解过程变得更加稳定，本发明方法采用基于L1范数的惩罚项将光源的稀疏性特点作为先验知识融入到光源的重建中，使得荧光断层成像重建问题转化为具有如下形式的优化问题模型：

$\min J\left(w\right)=\frac{1}{2}{\left|\right|\mathrm{Aw}-{\mathrm{\Φ}}^m\left|\right|}_2^2+\mathrm{\λ}{\left|\right|w\left|\right|}_1^1$

其中，w表示光源的分布信息，λ表示正则化参数。

步骤102：使用对偶变换将步骤1中所述的无约束条件最优化问题模型进行转化，得到相应的对偶模型：

$\max E(\mathrm{\α},v)=-\frac{1}{2}{\left|\right|\mathrm{\α}-{\mathrm{\Φ}}^m\left|\right|}^2+\frac{1}{2}{\left|\right|{\mathrm{\Φ}}^m\left|\right|}^2-{\mathrm{\δ}}_{\mathrm{\λ}}^{\∞}\left(v\right)$ subject to  v-A^Tα＝0

该对偶模型是一个约束条件最优化问题模型，目标解向量α的维数要远小于w的维数，这将会显著的提高成像的效率。

由于步骤101中所提出的最优化问题目标函数明显是不可微分的，这使得现有很多基于L2范数的荧光光源重建方法不能直接用来求解该模型。针对上述问题，本发明使用对偶变换，将步骤101中的非约束条件最优化问题模型进行转化，并得到与之等价的约束条件最优化问题模型。原模型中的目标变量w和对偶模型中的目标变量α的维数分别是M和N，M和N分别对应于系统矩阵的维数，还分别对应了表面荧光可测量结点的维数和成像区域所有结点的维数。显然，运用对偶变换将明显降低模型目标变量的维数。同时，转化所得到的约束条件最优化问题模型也为下一步使用增广拉格朗日函数技术奠定了基础。

步骤103：使用增广拉格朗日函数函数技术，建立步骤2所述对偶模型的增广拉格朗日函数：

$L(\mathrm{\α},v,w,\mathrm{\μ})=E(\mathrm{\α},v)-w^T(A^T\mathrm{\α}-v)-\frac{\mathrm{\μ}}{2}{\left|\right|A^T\mathrm{\α}-v\left|\right|}_2^2$

其中μ为惩罚向量，w为拉格朗日乘子。

在步骤102所得出的对偶最优化目标函数maxE(α，v)的基础上添加增广拉格朗日惩罚项 $-w^T(A^T\mathrm{\α}-v)-\frac{\mathrm{\μ}}{2}{\left|\right|A^T\mathrm{\α}-v\left|\right|}_2^2$ 就可以得到模型的增广拉格朗日函数；当μ＝0时，得到的是模型的拉格朗日函数。此时，表示荧光光源分布信息的w成为了拉格朗日乘子。

步骤104：将v表示为α的削顶函数增广拉格朗日函数的最大值函数可以化简为：

$\max L(\mathrm{\α},w,\mathrm{\μ})=\max -\frac{1}{2}{\left|\right|\mathrm{\α}-{\mathrm{\Φ}}^m\left|\right|}_2^2-\frac{\mathrm{\η}}{2}{\left|\right|{\mathrm{Shrink}}_{\mathrm{\λ}}(A^T\mathrm{\α}+\frac{w}{\mathrm{\η}})\left|\right|}_2^2$

从而降低了模型求解的复杂度。

将步骤104中的增广拉格朗日乘子函数看作是v的函数并展开为向量的形式，由于最大化函数中的每一项都能够被展开为N维，并且都包含v_j，由此可以得到：

$L(\mathrm{\α},v,w,\mathrm{\μ})=-\frac{1}{2}{\left|\right|\mathrm{\α}-{\mathrm{\Φ}}^m\left|\right|}^2+\frac{1}{2}{\left|\right|{\mathrm{\Φ}}^m\left|\right|}^2-\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\left(\frac{\mathrm{\μ}}{2}{(v_j-{(\frac{w}{\mathrm{\μ}}+A^T\mathrm{\α})}_j)}^2+{\mathrm{\δ}}_{\mathrm{\λ}}^{\∞}\left(v_j\right)\right)$

其中v_j表示向量的第j个元素。结合指示函数的性质，v(α)的最大值可以用包含α的削顶函数表示，并得到新的增广拉格朗日函数：

$\max L(\mathrm{\α},w,\mathrm{\μ})=\max -\frac{1}{2}{\left|\right|\mathrm{\α}-{\mathrm{\Φ}}^m\left|\right|}_2^2-\frac{\mathrm{\η}}{2}{\left|\right|{\mathrm{Shrink}}_{\mathrm{\λ}}(A^T\mathrm{\α}+\frac{w}{\mathrm{\η}})\left|\right|}_2^2$

从而降低了模型求解的复杂度。其中，Shirnk_λ(w)是软阈值函数，其定义为：

${\mathrm{Shirnk}}_{\mathrm{\λ}}\left(w\right)={\left(\max \left(\right|w_j|-\mathrm{\λ},0)\frac{w_j}{\left|w_j\right|}\right)}_j(j=1,...,n),$

削顶函数的定义为：

${\mathrm{CL}}_{\mathrm{\λ}}\left(y\right):={\left(\min \right|y_i|,\mathrm{\λ})\frac{y_i}{\left|y_i\right|})}_{j=1}^n$

步骤105：使用融合积极集策略的截断牛顿法求解增广拉格朗日函数的最大值maxL_α(α，w，μ)。

本发明的特性决定了重建荧光光源所需的计算量主要集中在求解步骤104所述增广拉格朗日函数的最大值，所以在本步骤中，通过综合使用牛顿法和积极集策略提高最大值的求解效率。下面对本步骤所涉及的子步骤进行详细说明：

步骤201：计算积极集A₊，该矩阵是系统矩阵A的子矩阵，由J₊＝{j∈{1，2，...，n}：|x_k+μ_kA^Tα|＞|λμ|_k}得到。μ和w所包含元素的数值大小直接决定了积极集的大小，所以在迭代过程中，随着μ的增大和w的稀疏性逐渐增强，牛顿法的每次迭代也会变得更快。

步骤202：集合积极集计算出增广拉格朗日函数的海森矩阵和雅克比矩阵：

${{\▿}^2}_{\mathrm{\α}}L\left(\mathrm{\α}\right)=I^m+{\mathrm{\μ}}_k{A}_{+}{A}_{+}^T$

${\▿}_{\mathrm{\α}}L\left(\mathrm{\α}\right)=\mathrm{\α}-{\mathrm{\Φ}}^m+{\mathrm{AST}}_{\mathrm{\λ}{\mathrm{\μ}}_k}(w_k+{\mathrm{\μ}}_k{A}^T\mathrm{\α})$

通常情况下，牛顿法最耗费计算时间的就是求解海森矩阵和雅克比矩阵，这里使用积极集策略将会显著减少求解时间。本发明中，通过选取海森矩阵的对角元素得到预共轭梯度算子P：

$P=\mathrm{diag}\left({\▿}^{2}J\left(\mathrm{\α}\right)\right)=\mathrm{diag}(-I_m-{\mathrm{\μ}}_k{A}_{+}{A}_{+}^T);$

步骤203：通过使用预共轭梯度法计算得至Δα，使用标准后向检索技术计算步长γ，实现对α的更新：

α_k+1＝α_k+Δα_k·γ；

步骤106：将作为w的最快下降方向，并结合软阈值函数技术，使用下面的式子对拉格朗日乘子向量进行更新：

$w_{k+1}={\mathrm{Shink}}_{{\mathrm{\λ\μ}}_k}(w_k+{\mathrm{\μ}}_k{A}^T{\mathrm{\α}}_k);$

将 $\∂\max L(\mathrm{\α},v,w,\mathrm{\μ})/\∂(\mathrm{\α},v,\mathrm{\μ})=-A^T{\mathrm{\α}}_t+v_{t}w$ 的最快下降方向，使用软阈值函数技术得到：

$w_{k+1}={\mathrm{Shink}}_{\mathrm{\λ}{\mathrm{\μ}}_k}(w_k+{\mathrm{\μ}}_k{A}^T{\mathrm{\α}}_k)$

保证了求解得到的光源分布具有显著的稀疏性。

步骤107：使用保守策略更新惩罚向量μ：

μ_k+1＝2μ_k

步骤108：计算目标函数值J(w)，如果||(J(w)^k-J(w)^k-1)||/||Φ^m||≥tol为真，计算k＝k+1并跳至步骤104，否则，结束计算；其中，tol为目标函数的收敛效率阈值。

运行结果

为了验证本发明方法的性能，我们利用小鼠CT数据集进行荧光光源重建的仿真实验。

请参阅图2，使用Tecplot软件显示仿真实验中所使用的基于小鼠CT数据集的非匀质模型。截取小鼠胸部2.5厘米长的一段躯干作为成像区域，该区域内包含肌肉、骨骼、心脏、肺、肝脏和脾脏六种组织(见表1)，每种组织分别定义了对应的光学特性参数，如图(a)所示。将该模型划分成结点数为4093，四面体个数为20068的网格网格，如图(b)所示。

表1不同组织的光学系数

请参阅图3。在本试验中，将一个光强设置为2nW/mm³，直径为0.6mm的球形点光源放置在小鼠体内，球心位置为(18.24mm，31.58mm，47.29mm)。其中光源位置和小鼠躯干的非匀质模型如图(a)所示，仿真得到的表面测量值如图(b)所示。

请参阅图4。为了更好地说明算法的性能，我们将本发明和现有最为流行的，被作为衡量重建性能的标准的基于Tikhonov正则化的成像方法进行对比。本发明使用了两种不同的成像条件突出其性能，一种是将小鼠全身作为成像区域的全域成像，还有一个是使用了可行光源区域策略(PSR)的局部成像。这里的局部成像区域定位为{(x，y，z)|46mm＜z＜50mm}。PSR的基本思想是人为地将光源限定在某个较小的区域内来减少成像的复杂性。该方法虽然能提高成像的质量，但是在实际应用如何选定这样一个较小的区域中往往是很困难的，会极大地限制荧光断层成像应用的范围。图4直观地显示了两种重建方法的性能。图中左侧一列为使用iso-surface的三维视图，右侧为在最大值点处进行切片得到的荧光光强分布图，黑色圆圈表示光源的正确位置。(a)图显示了使用Tikhonov正则化方法在全域范围内进行成像得到的重建结果：在定位精度方面，光源的中心位置偏离了正确位置(圆圈表示的区域)5.6681mm，在重建光源的最大能量方面，最大值只有0.00104nW/mm³，所需的成像时间是1301.2s；(b)图显示了使用Tikhonov正则化方法在局部范围内进行成像得到的重建结果：在定位精度方面，光源的中心位置偏离了正确位置有2.8132mm，在重建光源的最大能量方面，最大值只有0.00617nW/mm³，成像所需的时间是为216s；无论是区域成像还是局部成像Tikhonov正则化方法所得到的结果都显得过于平滑；(c)图显示了使用本发明在全域范围内进行成像得到的重建结果：在定位精度方面，光源的中心位置偏离了正确位置0.3995mm，在重建光源的最大能量方面，最大值为0.0143nW/mm³，成像所需时间为4.051s。本发明实现了亚毫米以下的定位精度，而且在成像效率方面相比局部的Tikhonov正则化方法提高了40倍以上。

请参阅图5和图6，在单光源的对比试验的基础上，进行双光源和三个光源的仿真实验，从而进一步说明本发明的多光源重建性能。在单光源实验的基础上添加了两个球形光源，大小均为0.6mm，对应区域的光强分别为1.6nW/mm³和2nW/mm³。通过定量分析可以得到，在双光源和三个光源试验中，本发明都获得了定位精度＜1.0毫米，重建时间＜5.0s的显著结果。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (10)

1.一种基于乘子法的自发荧光断层成像重建方法，包括步骤：

S1使用有限元分析方法扩散方程进行离散化，基于L1范数的惩罚项建立无约束条件最优化问题模型；

S2得到所述无约束条件最优化问题模型的对偶模型；

S3建立所述对偶模型的增广拉格朗日函数；

S4简化增广拉格朗日函数的最大值函数；

S5使用截断牛顿法求解增广拉格朗日函数的最大值；

S6将增广拉格朗日函数的梯度作为目标向量的最快下降方向对目标向量进行更新；

S7更新惩罚向量；

S8计算目标函数值J(w)，如果||(J(w)^k-J(w)^k-1)||/||Φ^m||≥tol为真，计算k＝k+1并跳至步骤S4，否则，结束计算，其中，tol为目标函数的收敛效率阈值，Φ^m为通过CCD相机所得到的被观测物体表面的光强值。

2.根据权利要求1所述的方法，其特征在于所述无约束条件最优化问题模型由下式表达：

$\min J\left(w\right)=\frac{1}{2}{\left|\right|\mathrm{Aw}-{\mathrm{\Φ}}^m\left|\right|}_2^2+\mathrm{\λ}{\left|\right|w\left|\right|}_1^1$

其中，其中A是系统矩阵，λ是正则化参数，J(w)为最小二乘最优化问题的目标函数，w表示光源的分布信息，其初始值设定为零。

3.根据权利要求2所述的方法，其特征在于所述对偶模型由下式表达：

$\begin{array}{cccc}\max E(\mathrm{\α},v)=-\frac{1}{2}{\left|\right|\mathrm{\α}-{\mathrm{\Φ}}^m\left|\right|}^2+\frac{1}{2}{\left|\right|{\mathrm{\Φ}}^m\left|\right|}^2-{\mathrm{\δ}}_{\mathrm{\λ}}^{\∞}\left(v\right)& \mathrm{subject}& \mathrm{to}& v-A^T\mathrm{\α}=0\end{array}$

该对偶模型是一个约束条件最优化问题模型，E(α，v)是新的目标函数，为指示函数，目标解向量α的维数要远小于w的维数。

4.根据权利要求3所述的方法，其特征在于所述增广拉格朗日函数由下式表达：

$L(\mathrm{\α},v,w,\mathrm{\μ})=E(\mathrm{\α},v)-w^T(A^T\mathrm{\α}-v)-\frac{\mathrm{\μ}}{2}{\left|\right|A^T\mathrm{\α}-v\left|\right|}_2^2$

其中μ为惩罚向量，w为拉格朗日乘子。

5.根据权利要求4所述的方法，其特征在于：通过削顶函数技术将v表示成α的函数：

$v\left(\mathrm{\α}\right)={\mathrm{CL}}_{\mathrm{\λ}}(\frac{w}{\mathrm{\μ}}+A^T\mathrm{\α})$

其中削顶函数的表达式为：

${\mathrm{CL}}_{\mathrm{\λ}}\left(y\right):=\left(\min \right|y_i|,\mathrm{\λ})\frac{y_i}{\left|y_i\right|}{)}_{j=1}^n.$

6.根据权利要求5所述的方法，其特征在于所述使用截断牛顿法求解增广拉格朗日函数的最大值包括：

计算积极集A₊；

集合积极集计算出增广拉格朗日函数的海森矩阵和雅克比矩阵：

通过使用预共轭梯度法计算得到Δα，使用标准后向检索技术计算步长γ，实现对α的更新。

7.根据权利要求6所述的方法，其特征在于所述雅克比矩阵和海森矩阵由下式计算得到：

${\▿}_{\mathrm{\α}}^{2}L\left(\mathrm{\α}\right)=I^m+{\mathrm{\μ}}_k{A}_{+}{A}_{+}^T$

${\▿}_{\mathrm{\α}}L\left(\mathrm{\α}\right)=\mathrm{\α}-{\mathrm{\Φ}}^m+{\mathrm{AST}}_{{\mathrm{\λ\μ}}_k}(w_k+{\mathrm{\μ}}_k{A}^T\mathrm{\α})$

其中Shirnk_λ(w)为标准软阈值函数，其表达式为：

${\mathrm{Shirnk}}_{\mathrm{\λ}}\left(w\right)={\left(\max \left(\right|w_j|-\mathrm{\λ},0)\frac{w_j}{\left|w_j\right|}\right)}_j(j=1,...,n).$

8.根据权利要求7所述的方法，其特征在于选取海森矩阵的对角元素得到预共轭梯度算子P：

$P=\mathrm{diag}\left({\▿}^{2}J\left(\mathrm{\α}\right)\right)=\mathrm{diag}({-I}_m-{\mathrm{\μ}}_k{A}_{+}{A}_{+}^T).$

9.根据权利要求6所述的方法，其特征在于：积极集A₊是系统矩阵A的子矩阵，包含了由J₊＝{j∈{1，2，...，n}：|x_k+μ_kA^Tα|＞|λμ|_k}所得到的积极集。

10.根据权利要求1所述的方法，其特征在于使用软阈值函数更新目标向量。