一种X射线CT图像的增广拉格朗日迭代重建方法

摘要

一种X射线CT图像的增广拉格朗日迭代重建方法，其特征在于：依次包括如下步骤：（1）获取CT设备的系统参数和低剂量扫描协议下的投影数据；（2）对步骤（1）中的投影数据进行逐个数据点上的方差估计，并对步骤（1）中的投影数据进行滤波反投影得到初始图像；（3）以步骤(2)中得到的初始图像作为迭代的初始图像进行迭代重建，根据迭代公式、 进行循环迭代，获得最终的重建图像。本发明同时提出了对上述迭代公式的优化算法。本发明适用性宽，迭代次数少，成像质量高。

说明书

技术领域

本发明涉及一种医学影像的图像处理技术领域，具体涉及一种X射线CT图像的增广拉格朗日迭代重建方法。 

背景技术

X射线CT扫描已经广泛应用于临床医学影像诊断，但是X射线辐射剂量（Radiation Dose）会对人体造成辐射损伤、甚至诱发恶性肿瘤等不可预测的伤害。为了降低对使用者的损害,在保证重建图像质量的前提下，最大限度地降低X射线剂量已成为医学CT成像领域的迫切要求。 

为了降低X射线辐射剂量，可以通过各种硬件技术及软件技术降低CT扫描中的X射线使用剂量。常见的方法有降低管电流、降低X射线曝光时间以及减少投影数据的采集角度。但是，通过降低X射线辐射剂量进行低剂量（Low Dose）X射线CT图像重建时会因投影数据中存在大量噪声而导致重建图像质量严重退化，难以满足临床诊断需求。 

目前，CT厂商提出了许多低剂量成像新技术，如管电流调制TCM（Tube Current Modulation）技术，可以在保证图像质量的前提下，一定程度上减少辐射剂量。除了CT厂商采用的硬件技术和优化的扫描协议来降低剂量外，在不改变现有硬件设备的情况下，基于软件的低剂量重建技术受到广泛研究，整体可分为三种策略进行：策略一是直接对CT图像进行噪声滤波；策略二是先对CT投影数据（Projection Data）进行滤波，再采用传统的滤波反投影法进行重建；策略三是CT图像迭代重建。 

策略一未考虑CT图像噪声的来源，属于后处理技术，存在的主要问题是：CT图像噪声分布特性复杂，难以实现高性能的滤波器。策略二直接对投影数据进行恢复，可以将投影数据的统计特性纳入到滤波器的设计中，使得这类算法具有高效且重建图像均匀性好等优点。策略三充分考虑投影数据的统计特性和CT系统的物理模型，特别适合于低剂量CT图像迭代重建。 

通过策略三实现CT图像重建，需要对相应的目标函数（Objective Function）进行正则化（Regularization）设计，并给出相应的求解技术。 

现有CT图像迭代重建中的正则化分为两类:非光滑的正则化和光滑的正 则化。CT图像迭代重建中的非光滑的正则化项，比如全变分和l₁正则化，使得相应的目标函数无法用传统的梯度下降（Gradient Descend）方法进行求解。虽然可以使用光滑的正则化函数进行近似，但是相应算法的计算量大且收敛速度慢。对于CT图像迭代重建中的光滑正则化，虽然可以用传统的梯度下降方法求解，比如非线性共轭梯度方法（NCG），但是计算量也往往很大且收敛速度较慢。因此，针对现有技术不足，提供一种适用性宽、计算简单的X射线CT图像重建方法以解决现有技术不足甚为必要。 

发明内容

本发明提供一种X射线CT图像的增广拉格朗日迭代重建方法,该方法适用性宽、迭代次数少，成像质量高。 

本发明的上述目的通过如下技术手段实现。 

一种X射线CT图像的增广拉格朗日迭代重建方法，依次包括如下步骤： 

（1）、获取CT设备的系统参数和低剂量扫描协议下的投影数据y； 

（2）、对步骤（1）中的投影数据y进行逐个数据点上的方差估计，并对步骤（1）中的投影数据y进行滤波反投影得到初始图像； 

（3）、以步骤(2)中得到的初始图像作为迭代的初始图像进行迭代重建，获得最终的重建图像。 

上述步骤（3）中迭代重建采用如下方法进行， 

所述迭代重建针对如下X射线CT图像重建模型： 

$\left(\begin{array}{cc}\underset{x}{\min }& \Psi \left(\mathrm{D\mu }\right)\\ s.t.& y=\mathrm{A\mu }\end{array}\right)---\left(1\right);$

其中，y=(y₁,y₂,…,y_M)^T为步骤（1）得到的投影数据，μ=(μ₁,μ₂,…,μ_N)^T为待重建图像，A为M×N维的系统矩阵，T表示矩阵转置，函数Ψ(Dμ)定义如下： 

$\Psi \left(\mathrm{D\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{\Phi }_r(\underset{p=1}{\overset{P}{\Sigma }}{\left|{[D}_p\mu {]}_r\right|}^m)---\left(2\right);$

其中κ_r>0控制图像的空间分辨率,D_p为L×N维的矩阵，p=1,2,...,P, 为实数集，Φ_r为保持边界的势函数； 

引入新的约束变量，（1）式写成另外一个等价形式： 

$\left(\begin{array}{cc}\underset{x,v}{\min }& \Psi \left(x\right)\\ s.t.& y=\mathrm{A\mu },x=\mathrm{D\mu }\end{array}\right)---\left(3\right);$

（3）式的增广拉格朗日函数为： 

$L(x,\mu ,\lambda ,\gamma )=\Psi \left(x\right)-{\lambda }^T(x-\mathrm{D\mu })+\frac{\beta }{2}{\left|\right|x-\mathrm{D\mu }\left|\right|}^2-{\gamma }^T(y-\mathrm{A\mu })+\frac{\eta }{2}{\left|\right|y-\mathrm{A\mu }\left|\right|}_W^2---\left(4\right);$

其中λ，γ分别为对应约束的拉格朗日乘子，β>0，η>0为惩罚参数，W=diag{w_i}为M×M对角矩阵，{w_i}为步骤（2）数据方差的倒数， ${\left|\right|x\left|\right|}_W^2=x^T\mathrm{Wx};$

对（4）式采用交替最小化方法进行求解，设在第k步迭代点为（x_k,μ_k)，对应的拉格朗日乘子为λ_k,γ_k，则下一个迭代点由下式得到： 

$x_{k+1}=\arg \underset{x}{\min }\Psi \left(x\right)+\frac{\beta }{2}{\left|\right|x-({\mathrm{D\mu }}_k+{\lambda }_k/\beta )\left|\right|}^2---\left(5\right);$

${\mu }_{k+1}=\arg \underset{\mu }{\min }-{\lambda }_k^T(x_{k+1}-\mathrm{D\mu })+\frac{\beta }{2}{\left|\right|x_{k+1}-\mathrm{D\mu }\left|\right|}^2-{\gamma }_k^T(y-\mathrm{A\mu })+\frac{\eta }{2}{\left|\right|y-\mathrm{A\mu }\left|\right|}_W^2---\left(6\right);$

更新拉格朗日乘子λ_k+1=λ_k-β(x_k+1-Dμ_k+1),γ_k+1=γ_k-η(y-Aμ_k+1)，循环执行公式（5）和公式（6），当循环次数达到预设的次数时停止迭代运算，得到最终的CT图像μ。 

上述对(5)式进行求解具体包括， 

先将(5)式写成多个一维的最小化形式： 

$x_{k+1}^r=\arg \underset{x^r}{\min }\Psi \left(x^r\right)\frac{{\mathrm{\beta \lambda }}_k^r}{2}{(x^r-{\rho }_k^r)}^2---\left(7\right);$

其中是ρ_k=Dμ_k+λ_k/β的第r个分量； 

再根据（7）式及Shrinkage法则对(7)式求解，得到的具体解。 

具体的，当m=1，Φ_r(x)=x/δ-log(1+x/δ),δ>0， $\Psi \left(\mathrm{D\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{\Phi }_r\left({\left|\right|{\left[\mathrm{D\mu }\right]}_r\left|\right|}_1\right)$时， 

(7)式的解为： 

$x_{k+1}^r=\frac{{\xi }_k^r+\sqrt{{\left({\xi }_k^r\right)}^2+4\delta {\left|\right|{\rho }_k^r\left|\right|}_1}}{2}\mathrm{sgn}\left({\rho }_k^r\right);$

其中${\xi }_k^r={\left|\right|{\rho }_k^r\left|\right|}_1-\delta -{\kappa }_r/\left(\mathrm{\delta \beta }\right).$

当$m=2,{\Phi }_r\left(x\right)=\sqrt{x},$$\Psi \left(\mathrm{R\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{\left|\right|\left[R\mu {]}_r\right||}_2,$时，(7)式的解为： 

$x_{k+1}^r=\max \{{\left|\right|{\rho }_k^r\left|\right|}_2-\frac{{\kappa }_r}{\xi },0\}\frac{{\rho }_k^r}{{\left|\right|{\rho }_k^r\left|\right|}_2};$

当m=1,Φ,(x)=x,$\Psi \left(\mathrm{R\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{\left|\right|{\left[R\mu \right]}_r\left|\right|}_1,$时, 

（7）式的解为： 

$x_{k+1}^r=\max \{{\left|\right|{\rho }_k^r\left|\right|}_1-\frac{{\kappa }_r}{\xi },0\}\mathrm{sign}\left({\rho }_k^r\right).$

上述（6）式具体采用自适应的非单调谱梯度算法进行计算，其迭代格式如下： 

μ_k+1=μ_k-α_kg_k(μ_k)          (9)； 

设置步长的自适应选取公式： 

其中τ<1为充分靠近1的常数,${\alpha }_k^{\mathrm{BB}1}=s_{k-1}^T{s}_{k-1}/{y_{k-1}}^T{s}_{k-1},$${\alpha }_k^{\mathrm{BB}2}={z_{k-1}}^T{s}_{k-1}/{z_{k-1}}^T{z}_{k-1},$Q_k(u)的梯度g_k=βD^T(Dμ-x_k+1)+D^Tλ_k+ηA^TW(y-Aμ)+A^Tγ_k,s_k-1=μ_k-μ_k-1,z_k-1=g_k-g_k-1，且步长α_k满足自适应的非单调线搜索条件，自适应非单调线搜索形式为： 

$Q_k({\mu }_k-{\mathrm{\alpha g}}_k)\le f_r-\mathrm{\alpha \theta }{g}_k^T{g}_k---\left(11\right);$

其中0<θ<1为常数,α>0为初始步长,f_r为参考函数值； 

循环执行（9）式，当满足终止条件时迭代终止，得到（6）式的解。 

上述步骤（1）中获取的CT设备的系统参数包括X射线入射光子强度I₀ 和系统电子噪声的方差

上述步骤（2）对步骤（1）中的投影数据y进行滤波反投影得到初始图像具体是采用FBP方法对投影数据y进行重建得到初始图像。 

本发明的X射线CT图像的增广拉格朗日迭代重建方法，适用于一般的正则化技术，包括光滑和非光滑的正则化技术都适用，迭代次数少，可实现快速X射线CT图像高质量迭代重建。 

附图说明

利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限制。 

图1为本发明方法的流程示意图； 

图2为本发明实施例所采用的XCAT体模图像。 

图3为通过滤波反投影法重建的初始图像； 

图4为通过本发明的方法重建的图像； 

图5为通过非线性共轭梯度方法重建的图像； 

图6为通过预条件共轭梯度方法重建的图像； 

图7为三种算法的信噪比曲线。 

具体实施方式

结合以下实施例对本发明作进一步描述。 

实施例1。 

一种X射线CT图像的增广拉格朗日迭代重建方法，依次包括如下步骤： 

（1）、获取CT设备的系统参数和低剂量扫描协议下的投影数据y；获取的CT设备的系统参数主要包括X射线入射光子强度I₀和系统电子噪声的方差等。 

（2）、对步骤（1）中的投影数据y进行逐个数据点上的方差估计，并采用FBP方法对步骤（1）中的投影数据y进行重建得到初始图像。 

（3）、以步骤(2)中得到的初始图像作为迭代的初始图像进行迭代重建，获得最终的重建图像。 

步骤（3）中迭代重建采用如下方法进行， 

迭代重建针对如下X射线CT图像重建模型： 

$\left(\begin{array}{cc}\underset{x}{\min }& \Psi \left(\mathrm{D\mu }\right)\\ s.t.& y=\mathrm{A\mu }\end{array}\right)---\left(1\right);$

其中，y=(y₁,y₂,…,y_M)^T为步骤（1）得到的投影数据，μ=(μ₁,μ₂,…,μ_N)^T为待重建图像，A为M×N维的系统矩阵，T表示矩阵转置，函数Ψ(Dμ)定义如下： 

$\Psi \left(\mathrm{D\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{\Phi }_r(\underset{p=1}{\overset{P}{\Sigma }}{\left|{[D}_p\mu {]}_r\right|}^m)---\left(2\right);$

其中κ_r>0控制图像的空间分辨率,D_p为L×N维的矩阵，p=1,2,...,P, 为实数集，Φ_r为保持边界的势函数； 

引入新的约束变量，（1）式写成另外一个等价形式： 

$\left(\begin{array}{cc}\underset{x,v}{\min }& \Psi \left(x\right)\\ s.t.& y=\mathrm{A\mu },x=\mathrm{D\mu }\end{array}\right)---\left(3\right);$

（3）式的增广拉格朗日函数为： 

$L(x,\mu ,\lambda ,\gamma )=\Psi \left(x\right)-{\lambda }^T(x-\mathrm{D\mu })+\frac{\beta }{2}{\left|\right|x-\mathrm{D\mu }\left|\right|}^2-{\gamma }^T(y-\mathrm{A\mu })+\frac{\eta }{2}{\left|\right|y-\mathrm{A\mu }\left|\right|}_W^2---\left(4\right);$

其中λ，γ分别为对应约束的拉格朗日乘子，β>0，η>0为惩罚参数，W=diag{w_i}为M×M对角矩阵，{w_i}为步骤（2）数据方差的倒数， ${\left|\right|x\left|\right|}_W^2=x^T\mathrm{Wx};$

对（4）式采用交替最小化方法进行求解，设在第k步迭代点为(x_k,μ_k)，对应的拉格朗日乘子为λ_k，γ_k，则下一个迭代点由下式得到： 

$x_{k+1}=\arg \underset{\text{x}}{\min }\Psi \left(x\right)+\frac{\beta }{2}{\left|\right|x-({\mathrm{D\mu }}_k+{\lambda }_k/\beta )\left|\right|}^2---\left(5\right);$

${\mu }_{k+1}=\arg \underset{\mu }{\min }{-\lambda }_k^T(x_{k+1}-\mathrm{D\mu })+\frac{\beta }{2}{\left|\right|x_{k+1}-\mathrm{D\mu }\left|\right|}^2-{\gamma }_k^T(y-\mathrm{A\mu })+\frac{\eta }{2}{\left|\right|y-\mathrm{A\mu }\left|\right|}_W^2---\left(6\right);$

更新拉格朗日乘子λ_k+1=λ_k-β(x_k+1-Dμ_k+1),γ_k+1=γ_k-η(y-Aμ_k+1)，循环执行公式（5）和公式（6），当循环次数达到预设的次数时停止迭代运算，得到最终 的CT图像μ。 

本发明的方法给出的CT图像重建模型包含了一般的正则化项，对于包含一般正则化形式的CT图像重建模型都可以用本发明的方法进行求解。故重建模型的目标函数不要求光滑，对于光滑和非光滑的目标函数本发明的方法都可以求解，因此本发明适用性宽。 

综上所述，本发明的X射线CT图像的增广拉格朗日迭代重建方法，适用于一般的正则化技术，包括光滑和非光滑的正则化技术都适用，迭代次数少，可实现快速X射线CT图像高质量迭代重建。 

实施例2。 

一种X射线CT图像的增广拉格朗日迭代重建方法，依次包括如下步骤：（1）、获取CT设备的系统参数和低剂量扫描协议下的投影数据y；获取的CT设备的系统参数主要包括X射线入射光子强度I₀和系统电子噪声的方差等。 

（2）、对步骤（1）中的投影数据y进行逐个数据点上的方差估计，并采用FBP方法对步骤（1）中的投影数据y进行重建得到初始图像。 

（3）、以步骤（2）中得到的初始图像作为迭代的初始图像进行迭代重建，获得最终的重建图像。 

步骤（3）中迭代重建采用如下方法进行， 

迭代重建针对如下X射线CT图像重建模型： 

$\left(\begin{array}{cc}\underset{x}{\min }& \Psi \left(\mathrm{D\mu }\right)\\ s.t.& y=\mathrm{A\mu }\end{array}\right)---\left(1\right);$

其中，y=(y₁,y₂,…,y_M)^T为步骤（1）得到的投影数据，μ=(μ₁,μ₂,…,μ_N)^T为待重建图像，A为M×N维的系统矩阵，T表示矩阵转置，函数Ψ（Dμ）定义如下： 

$\Psi \left(\mathrm{D\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{\Phi }_r(\underset{p=1}{\overset{P}{\Sigma }}{\left|{[D}_p\mu {]}_r\right|}^m)---\left(2\right);$

其中κ_r>0控制图像的空间分辨率,D_p为L×N维的矩阵，p=1,2,...,P, 为实数集，Φ_r为保持边界的势函数； 

引入新的约束变量，（1）式写成另外一个等价形式： 

$\left(\begin{array}{cc}\underset{x,v}{\min }& \Psi \left(x\right)\\ s.t.& y=\mathrm{A\mu },x=\mathrm{D\mu }\end{array}\right)---\left(3\right);$

（3）式的增广拉格朗日函数为： 

$L(x,\mu ,\lambda ,\gamma )=\Psi \left(x\right)-{\lambda }^T(x-\mathrm{D\mu })+\frac{\beta }{2}{\left|\right|x-\mathrm{D\mu }\left|\right|}^2-{\gamma }^T(y-\mathrm{A\mu })+\frac{\eta }{2}{\left|\right|y-\mathrm{A\mu }\left|\right|}_W^2---\left(4\right);$

其中λ，γ分别为对应约束的拉格朗日乘子，β>0，η>0为惩罚参数，W=diag{w_i}为M×M对角矩阵，{w_i}为步骤（2）数据方差的倒数， ${\left|\right|x\left|\right|}_W^2=x^T\mathrm{Wx};$

对（4）式采用交替最小化方法进行求解，设在第k步迭代点为(x_k,μ_k)，对应的拉格朗日乘子为λ_k，γ_k，则下一个迭代点由下式得到： 

$x_{k+1}=\arg \underset{x}{\min }\Psi \left(x\right)+\frac{\beta }{2}{\left|\right|x-({\mathrm{D\mu }}_k+{\lambda }_k/\beta )\left|\right|}^2---\left(5\right);$

${\mu }_{k+1}=\arg \underset{\mu }{\min }{-\lambda }_k^T(x_{k+1}-\mathrm{D\mu })+\frac{\beta }{2}{\left|\right|x_{k+1}-\mathrm{D\mu }\left|\right|}^2-{\gamma }_k^T(y-\mathrm{A\mu })+\frac{\eta }{2}{\left|\right|y-\mathrm{A\mu }\left|\right|}_W^2---\left(6\right);$

更新拉格朗日乘子λ_k+1=λ_k-β(x_k+1-Dμ_k+1),γ_k+1=γ_k-η(y-Aμ_k+1)，循环执行公式（5）和公式（6），当循环次数达到预设的次数时停止迭代运算，得到最终的CT图像μ。 

对(5)式进行精确求解具体包括， 

先将(5)式写成多个一维的最小化形式： 

$x_{k+1}^r=\arg \underset{x^r}{\min }\Psi \left(x^r\right)+\frac{{\mathrm{\beta \lambda }}_k^r}{2}{(x^r-{\rho }_k^r)}^2---\left(7\right);$

其中是ρ_k=Dμ_k+λ_k/β的第r个分量； 

再根据（7）式及Shrinkage法则对(7)式求解，得到的具体解。 

以m=1，Φ_r(x)=x/δ-log(1+x/δ),δ>0， $\Psi \left(\mathrm{D\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{{\Phi }_r\left(\right|\left|{\left[\mathrm{D\mu }\right]}_r\right||}_1)$为例，此例为光滑的目标函数。 

(7)式的解为： 

$x_{k+1}^r=\frac{{\xi }_k^r+\sqrt{{\left({\xi }_k^r\right)}^2+4\delta {\left|\right|{\rho }_k^r\left|\right|}_1}}{2}\mathrm{sgn}\left({\rho }_k^r\right)---\left(8\right);$

其中${\xi }_k^r={\left|\right|{\rho }_k^r\left|\right|}_1-\delta -{\kappa }_r/\left(\mathrm{\delta \beta }\right).$

（6）式具体采用自适应的非单调谱梯度算法进行计算，其迭代格式如下： 

μ_k+1=μ_k-α_kg_k（μ_k）          (9)； 

设置步长的自适应选取公式： 

其中τ<1为充分靠近1的常数,${\alpha }_k^{\mathrm{BB}1}=s_{k-1}^T{s}_{k-1}/{y_{k-1}}^T{s}_{k-1},$${\alpha }_k^{\mathrm{BB}2}={z_{k-1}}^T{s}_{k-1}/{z_{k-1}}^T{z}_{k-1},$Q_k(u)的梯度g_k=βD^T(Dμ-x_k+1)+D^Tλ_k+ηA^TW(y-Aμ)=A^Tγ_k,s_k-1=μ_k-μ_k-1,z_k-1=g_k-g_k-1，且步长α_k满足自适应的非单调线搜索条件，自适应非单调线搜索形式为： 

$Q_k({\mu }_k-{\mathrm{\alpha g}}_k)\le f_r-\mathrm{\alpha \theta }{g}_k^T{g}_k---\left(11\right);$

其中0<θ<1为常数,α>0为初始步长,f_r为参考函数值； 

循环执行（9）式，当满足终止条件时迭代终止，得到（6）式的解。 

本发明的方法给出的CT图像重建模型包含了一般的正则化项，对于包含一般正则化形式的CT图像重建模型都可以用本发明的方法进行求解。故重建模型的目标函数不要求光滑，对于光滑和非光滑的目标函数本发明的方法都可以求解，因此本发明适用性宽。 

具体提出了自适应步长选取公式，并将自适应的非单调线搜索与提出的步长公式结合起来给出的求解无约束优化问题的新方法。进一步，将求解无约束优化问题的方法与增广拉格朗日方法结合起来用于CT图像重建，大大减少了计算的复杂程度，重建图像质量高。 

综上所述，本发明的X射线CT图像迭代重建方法，适用于一般的正则化技术，包括光滑和非光滑的正则化技术都适用，迭代次数少，可实现快速X 射线CT图像高质量迭代重建。 

实施例3 

一种X射线CT图像的增广拉格朗日迭代重建方法，其他内容与实施例2相同，不同之处在于：针对如下非光滑的目标函数进行(5)式求解，具体如下： 

对于（5）式，我们可以对其精确求解。(5)式可以写成多个一维的最小化问题： 

$x_{k+1}^r=\arg \underset{x^r}{\min }\Psi \left(x^r\right)+\frac{{\mathrm{\beta \lambda }}_k^r}{2}{(x^r-{\rho }_k^r)}^2---\left(7\right);$

其中是ρ_k=Dμ_k+λ_k/β的第r个分量。 

在本实例中，设置非光滑的目标函数$\Psi \left(\mathrm{R\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{\left|\right|{\left[R\mu \right]}_r\left|\right|}_2,$时， 

(7)式的解为： 

$x_{k+1}^r=\max \{{\left|\right|{\rho }_k^r\left|\right|}_2-\frac{{\kappa }_r}{\xi },0\}\frac{{\rho }_k^r}{{\left|\right|{\rho }_k^r\left|\right|}_2}---\left(8\right);$

可见，本发明的X射线CT图像的增广拉格朗日迭代重建方法，适用于一般的正则化技术，包括光滑和非光滑的正则化技术都适用。 

实施例4 

一种X射线CT图像的增广拉格朗日迭代重建方法，其他内容与实施例2相同，不同之处在于：针对如下非光滑的目标函数进行(5)式求解，具体如下： 

对于（5）式，我们可以对其精确求解。(5)式可以写成多个一维的最小化问题： 

$x_{k+1}^r=\arg \underset{x^r}{\min }\Psi \left(x^r\right)+\frac{{\mathrm{\beta \lambda }}_k^r}{2}{(x^r-{\rho }_k^r)}^2---\left(7\right);$

其中是ρ_k=Dμ_k+λ_k/β的第r个分量。 

在本实例中，设置非光滑的目标函数：m=1,Φ,(μ)=μ, $\Psi \left(\mathrm{R\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{\left|\right|{\left[R\mu \right]}_r\left|\right|}_1,$时, 

（7）式的解为：$x_{k+1}^r=\max \{{\left|\right|{\rho }_k^r\left|\right|}_1-\frac{{\kappa }_r}{\xi },0\}\mathrm{sign}\left({\rho }_k^r\right)---\left(8\right);$

可见，本发明的X射线CT图像的增广拉格朗日迭代重建方法，适用于一般的正 则化技术，包括光滑和非光滑的正则化技术都适用。 

实施例5。 

本实施例以计算机仿真的数字体模数据为例来描述本发明所述方法的具体实施过程，如图1所示，本实施例的实施过程如下。 

（1）采用图2所示的XCAT数字体模图像作为本发明的计算机仿真实验对象。体模图像的大小设为512×512，模拟CT机的X射线源到旋转中心和探测器的距离分别为570mm和1040mm，旋转角在[0,2π]间采样值为1160，每个采样角对应672个探测器单元，探测器单元的大小为1.40mm。通过CT系统仿真生成大小为1160×672的投影数据y。 

其中X射线的入射光子强度I₀为2.5×10⁴，系统电子噪声的方差为10.0（在实际的CT数据采集中，投影数据和系统参数即入射光子强度I₀和系统电子噪声的方差均可以直接获取）。 

（2）对步骤1中模拟生成的CT投影数据y和系统参数I₀和对投影数据y使用传统扇形束滤波反投影FBP（Filtered Back-Projection）算法进行重建得到迭代的初始图像，传统扇形束滤波反投影算法所采用的截止频率为奈奎斯特频率。 

设置μ₀为初始图像,x₀=Dμ₀,β=2¹⁰,η=2¹³,λ₀=γ₀=0,采用（5）式和（6）式进行迭代计算得到。 

$x_{k+1}=\arg \underset{x}{\min }\Psi \left(x\right)+\frac{\beta }{2}{\left|\right|x-({\mathrm{D\mu }}_k+{\lambda }_k/\beta )\left|\right|}^2---\left(5\right);$

${\mu }_{k+1}=\arg \underset{\mu }{\min }{-\lambda }_k^T(x_{k+1}-\mathrm{D\mu })+\frac{\beta }{2}{\left|\right|x_{k+1}-\mathrm{D\mu }\left|\right|}^2-{\gamma }_k^T(y-\mathrm{A\mu })+\frac{\eta }{2}{\left|\right|y-\mathrm{A\mu }\left|\right|}_W^2---\left(6\right);$

对于（5）式，我们可以对其精确求解。(5)式可以写成多个一维的最小化问题： 

$x_{k+1}^r=\arg \underset{x^r}{\min }\Psi \left(x^r\right)+\frac{{\mathrm{\beta \lambda }}_k^r}{2}{(x^r-{\rho }_k^r)}^2---\left(7\right);$

其中是ρ_k=Dμ_k+λ_k/β的第r个分量。在本实例中，m=1， Φ_r(x)=x/δ-log(1+x/δ),δ>0，$\Psi \left(\mathrm{D\mu }\right)=\underset{r}{\Sigma }{\kappa }_r{\Phi }_r\left({\left|\right|{\left[\mathrm{D\mu }\right]}_r\left|\right|}_1\right),$

由上式以及Shrinkage法则可以得到(7)式的解为： 

$x_{k+1}^r=\frac{{\xi }_k^r+\sqrt{{\left({\xi }_k^r\right)}^2+4\delta {\left|\right|{\rho }_k^r\left|\right|}_1}}{2}\mathrm{sgn}\left({\rho }_k^r\right)---\left(8\right);$

其中${\xi }_k^r={\left|\right|{\rho }_k^r\left|\right|}_1-\delta -{\kappa }_r/\left(\mathrm{\delta \beta }\right).$

对于（6）式，采用自适应的非单调谱梯度算法，其迭代格式如下： 

μ_k+1=μ_k-α_kg_k(μ_k)          (9)； 

其中，设置一个步长的自适应选取公式： 

其中τ<1为充分靠近1的常数,${\alpha }_k^{\mathrm{BB}1}=s_{k-1}^T{s}_{k-1}/{y_{k-1}}^T{s}_{k-1},$${\alpha }_k^{\mathrm{BB}2}={z_{k-1}}^T{s}_{k-1}/{z_{k-1}}^T{z}_{k-1},$Q_k(u)的梯度g_k=βD^T(Dμ-x_k+1)+D^Tλ_k+ηA^TW(y-Aμ)+A^Tγ_k,s_k-1=μ_k-μ_k-1,z_k-1=g_k-g_k-1,且步长α_k满足自适应的非单调线搜索条件。循环执行（9）式满足一定的终止条件迭代终止，可以得到（6）式的一个解。 

（3）更新拉格朗日乘子λ_k+1=λ_k-β(x_k+1-Dμ_k+1),γ_k+1=γ_k-η(y-Aμ_k+1)。 

（4）循环执行，当循环次数达到预设的次数时即停止迭代运算，得到最终的CT图像。 

本发明的方法给出的CT图像重建模型包含了一般的正则化项，对于包含一般正则化形式的CT图像重建模型都可以用本发明的方法进行求解。故重建模型的目标函数不要求光滑，对于光滑和非光滑的目标函数本发明的方法都可以求解，因此本发明适用性宽。 

具体计算提出了自适应步长选取公式，并将自适应的非单调线搜索与提出的步长公式结合起来给出的求解无约束优化问题的新方法。进一步将求解无约束优化问题的方法与增广拉格朗日方法结合起来用于CT图像重建，大大减少了计算的复杂程度，重建图像质量高。 

为了验证本发明所示方法的效果，对模拟生成的同组CT投影数据直接采用滤波反投影的方法进行重建得到迭代初始图像，如图3所示。将迭代初始图像分别使用本发明公开的重建方法以及经典的非线性共轭梯度方法和预条件共轭梯度方法进行重建，分别得到重建结果为图4，图5，图6所示。 

比较图4、图5和图6可以看出，图4相对于图5和图6图像质量大为改善，可见本发明公开的方法重建的图像质量高。 

图7给出了三种方法的信噪比曲线，可以看出本发明公开的重建方法的收敛速度比非线性共轭梯度方法快，而且与非线性共轭梯度方法和预条件共轭梯度方法相比本发明提出的方法重建的图像质量明显要好。 

综上所述，本发明的X射线CT图像的增广拉格朗日迭代重建方法，适用于一般的正则化技术，包括光滑和非光滑的正则化技术都适用，迭代次数少，可实现快速X射线CT图像高质量迭代重建。 

最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。