超拉普拉斯重叠组稀疏先验的稀疏角度CT 重建

齐子文，孔慧华*，李佳欣，潘晋孝

1 中北大学数学学院，山西 太原 030051；

2 信息探测与处理山西省重点实验室，山西 太原 030051

1 引言

X 射线计算机断层扫描(computed tomography,CT)成像技术可以获取人体和物体的内部结构信息，已广泛应用于工业无损检测、医学临床诊断等领域[1]。然而，在临床检查时额外的X 射线辐射可能导致癌症和其他遗传病变[2]。因此，降低辐射剂量和重建高质量图像是目前的研究重点。降低X 射线辐射量的方法主要分两方面：第一种方法是减少每个投影中的X 射线曝光，第二种方法是减少投影采样数，即稀疏投影[3]。第一种方法会产生含有噪声的投影，重建图像含有大量的噪声和伪影[4]。第二种方法是应用稀疏角度CT，但由于投影数据的采样不足，重建图像会产生严重的伪影[5]。

一些传统的CT 图像重建算法比如滤波反投影(filtered back projection,FBP)算法和联合代数重建(simultaneous algebraic reconstruction technique,SART)算法，这两个基础算法无法对不完整的投影数据重建出高质量的图像。压缩感知(Compressed sensing,CS)理论的出现使得图像处理中的不适定问题可以有效地解决[6]。基于CS 的图像重建算法主要利用图像的先验知识，特别是图像的稀疏性，如全变分[7]、字典学习[8]、小波变换等[9]。近年来，全变分(total variational,TV)算法被广泛应用于不同场景下的CT 重建，该算法可以很好地平滑图像以达到去噪的目的[10]。Sidky 等将TV 最小化引入到不完全投影数据重建中，提出了一种基于投影到凸集(projections onto convex sets,POCS)的高效迭代重建算法TV-POCS[11]。虽然TV 正则项能够很好地恢复图像的边缘信息，但是过渡的平滑图像会导致图像产生阶梯伪影。为了克服这一现象，Zhang 等将非零像素稀疏性和图像连续性假设两个正则项来代替TV 正则项，提出了一种基于大规模线性规划的稀疏角度层析重建算法[12]。连祥媛等提出一种多通道联合的广义全变分 (Total Generalized Variational,TGV)的能谱CT 迭代重建算法，利用通道之间的相关性来消除重建过程中产生的噪声和阶梯伪影[13]。Peng 等将组稀疏(group sparse representation,GSR)正则化引入稀疏角度CT 重建，利用图像中相似补丁构成的图像组来表示GSR 的基本单元，并且成功将其应用于稀疏角度CT 重建中[14]。Jon 等将重叠组稀疏(overlapping group sparsity,OGS)和超拉普拉斯先验(hyper-Laplacian,HL)相结合，提出了一种基于重叠组稀疏超拉普拉斯先验(overlapping group sparsity on hyper-Laplacian,OGSHL)的图像去噪算法[15]。

本文将OGS 正则项和HL 正则项应用于稀疏角度CT 图像重建。相比于传统的TV 正则项，OGS 正则项改进图像梯度的稀疏性，由于重叠组采用图像梯度的结构信息作为图像梯度稀疏性的度量标准，因此该算法能够在恢复图像边缘的同时消除阶梯伪影，基于HL 先验可以很好地近似图像梯度的重尾分布，对于重建图像细节的保留有着更好的效果。

2 理论推导

2.1 基于TV 的稀疏角度CT 重建算法

稀疏角度CT 图像重建与传统CT 图像重建的过程基本相同，设稀疏角度CT 重建中需要重建的图像为f，从探测器上获得的投影数据为p，数据采集的投影矩阵用A表示，稀疏角度CT 重建问题可表示为

在稀疏角度CT 重建的情况下，由于投影采样不足会使重建图像产生噪声和伪影。为了重建高质量图像，使用CS 理论将图像进行稀疏变换，例如对重建图像进行离散梯度变换(discrete gradient transform,DGT)，则任意图像的DGT 定义为

则对于任意图像f，其TV 表达式为

其中∇x，∇y分别表示水平方向和竖直方向上的梯度算子，其表达式分别为

故基于TV 正则项的稀疏角度CT 重建模型为

其中：第一项为数据保真项，A为投影矩阵，p为投影数据，µ为保真项系数。第二项为TV 正则项，λ为正则项系数。

2.2 基于OGS-TV 的稀疏角度CT 重建算法

TV 正则项已被广泛应用，其虽然能够保留图像的锐利边缘，但是过渡的平滑图像会导致图像产生阶梯伪影。为了进一步提升图像的质量，消除重建中产生的阶梯伪影，Liu 等将一维信号的OGS 推广到二维图像中并成功应用于图像去噪[16]。

为了描述图像梯度的结构稀疏性，首先定义一个二维图像的像素组，形式如

其中：d(i,j),K是矩阵中所有元素组成的列向量，即：d(i,j),K=。根据式(8)可知当K=1 时，该正则项就转变为传统的TV 正则项，则OGS-TV 正则项可定义为

由于OGS 正则项考虑了每个像素点梯度的K×K方形领域内的信息，是一个非局部的稀疏先验，因此OGS-TV 正则项不仅拥有和TV 正则项一样的图像边缘保留能力，还能克服TV 重建中图像边缘产生的阶梯伪影。故基于OGS-TV 的稀疏角度CT 重建算法为

其中：第一项为数据保真项，第二项为OGS-TV 正则项，λ为正则项系数。

2.3 基于OGS-HL 的稀疏角度CT 重建算法

OGS-TV 虽然可以克服图像的阶梯伪影，但它难以恢复图像中较为复杂的纹理和细节。为进一步扩展OGS-TV，Kong 等强调图像的梯度服从重尾分布[17]，并且超拉普拉斯算子先验将会使得重尾分布获得更好的近似，提出的HL 先验的模型为

其中：‖·‖q表示拟范数lq，0＜q＜1。q为超拉普拉斯参数，通常参数q的范围为 0.5≤q≤1。根据式(11)，，将HL 先验代入OGS正则项得

根据式(12)，当超拉普拉斯参数q=1时，该正则项就转变为OGS-TV 正则项。OGS-HL 的正则项可表示为

故基于OGS-HL 的稀疏角度CT 重建算法，算法模型为

其中：第一项为数据保真项，第二项为OGS-HL 正则项，λ为正则项系数。

3 模型求解

OGS-HL 算法的模型是一个综合算法模型，TV，OGS-TV，OGS-HL 算法可以通过对模型设置不同的参数值来分别表示，如表1 所示。

表1 不同算法的参数设置Table 1 Parameter settings for different algorithms

由于算法模型包含的矩阵数据量大并且复杂，因此本文采用交替方向乘子 (alternating direction method of multipliers,ADMM)算法[18]，将模型转化为多个子问题来求解，式(14)引入变量u，将模型转化为

上述优化函数可分解为2 个子问题来分别更新f和u：

其中：式(17)采用梯度下降法来更新f：

式(18)进一步采用ADMM 算法来求解，引入辅助变量x1和x2，引入拉格朗日乘子w1，w2和惩罚参数δ，转化为以下包含多个变量的增广拉格朗日函数来更新：

Step1：更新x1和x2，根据式(20)，x1和x2的子问题可以由式(21)和式(22)来更新：

该子问题的结构比较复杂，故采用主分量最小化法(majorization minimization,MM)来求解上述最小化问题[19]。由于2 个子问题是相同的结构，下面只阐述式(21)的更新步骤。

为了方便后续计算，式(21)可简化为

首先，需找到φOH(x1)的优化目标函数，根据均值不等式，可以得到：

其中：a(i,j),K是 φOH(x1)中的任意一项，对于任意的a(i,j),K，x1(i,j),K≠0，并且a(i,j),K=x1(i,j),K时，对每一个稀疏组，代入式(24)可以得到 φOH(x1)的优化目标函数：

根据式 (24)能够证明 ψ(x1,a)≥φ(x1)，并且ψ(a,a)=φ(a)。通过简单的计算，ψ(x1,a)可以写为

其中：C 是一个常数，Λ(a)是对角矩阵，其对角的元素为

其中：r,t=1,2,...,n,l=1,2,...,n2，Λ(a)是通过二维卷积运算来求解，则式(23)中R(x1)的优化目标函数如：

式(28)可以得出 ξ(x1,a)≥R(x1)，ξ(a,a)=R(a)。故 ξ(x1,a)也满足MM 算法的前提条件，为了极小化R(x1)，使用MM 算法迭代，初始化，反复极小化替代函数，得到

根据欧拉-拉格朗日方程，式(29)可写为

其中：E是一个全1 的单位矩阵，S(x1)=diag(|x1|2q-2)。

Step2：更新w1和w2，拉格朗日乘子w1和w2的更新公式为

Step3：更新u，u的子问题的更新公式为

上述u子问题的优化目标是二次泛函，故可以求解下述等效正态方程：

显然，式(34)含有二维卷积运算，因此不能直接求出u子问题的最优解，空域中的卷积运算在频域中会转化为普通的乘法运算，故采用二维快速傅里叶变换与逆变换来求解：

其中：F和F-1分别是快速傅里叶变换和快速傅里叶逆变换。

基于OGS-HL 的稀疏角度CT 重建算法步骤：

4 实验结果

4.1 实验设置

为验证OGS-HL 算法在稀疏角度CT 重建下的效果，本文选取其他经典的算法与OGS-HL 算法作对比。实验采用仿真模型和真实模型来验证算法。实验配置为：11th Gen Intel (R) Core (TM) i5-11260H @ 2.60 GHz 的CPU，Nvidia RTX3050ti (4 GB)的GPU，16 GB 的内存，实验软件为MATLAB R2019b。本文算法与选取的对比算法均在MATLAB 和使用MEX 函数编译的C++的混合模式中实现，接口在MATLAB中实现。

4.2 仿真小鼠实验

仿真实验采用MOBY 软件生成的小鼠胸腔截面作为测试模型，在血液中注入1.2%的碘造影剂，胸腔模型尺寸为20 mm×20 mm，分辨率为512×512，如图1 所示。

图1 仿真小鼠胸腔模型和碘造影剂Fig.1 The simulated mouse thorax phantom and the iodine contrast agent

投影数据由等距扇形束几何采集，X 射线源到旋转中心的距离为100 mm，物体半径为10 mm，每个探测器长度为20 mm，探测器有320 个单元，每个单元的长度为0.0625 mm。每条X 射线发射的光子数为5×104，计算沿每条X 射线路径的预期光子数，为了模拟数据噪声，根据泊松分布生成随机数，其中方差为上述期望的光子数，然后进行对数运算，得到无噪声和有噪声的投影数据。使用分裂布雷格曼(Split-Bregman,SB)算法的重建的无噪声图像作为参考图像[20]。在360°范围内进行等间隔采样，设置的采样间隔角度分别为6°、4°、3°和2°，得到的稀疏投影角度个数为60、90、120 和180。本节实验采用归一化均方根误差(NRMSE)、峰值信噪比(PSNR)和结构相似度指数(SSIM)来评价噪声投影数据中重建图像的质量。

在OGS-TV 算法中，经过多次实验，参数为µ=50，λ=2×10-4，δ=1，q=1，K=3。在OGS-HL算法中，参数为 µ=50，λ=2×10-4，δ=1，q=0.8，K=3。

图2 展示了6 种不同算法下不同角度的仿真小鼠重建图像，图3 展现了仿真小鼠骨的部分ROI 区域。从图2 可以看出，FBP 和SART 两种算法在四个不同的角度下重建的图像中含有大量的伪影和噪声。从图3 的ROI 区域中仿真小鼠的骨切面可以看出TV 算法重建的图像出现阶梯伪影，使得骨切面的边缘变得不清晰。TGV 算法虽然能够克服TV 算法出现的问题，但对60、90 和120 角度下骨切面周围产生的条纹伪影去除效果较差。OGS-TV 算法在稀疏角度下重建的图像都有较大的改进，但对噪声的抑制效果一般。OGS-HL 算法不仅能够消除稀疏角度下产生的伪影，并且能够很好地抑制噪声，对小鼠骨切面周围的噪声消除以及骨切面的边缘保护有着良好的效果。

图2 仿真小鼠重建结果: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL。从上到下依次为投影幅数为60、90、120 和180 的重建图像Fig.2 The reconstruction results of the simulated mouse model: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL.From top to bottom,the reconstructed images are from 60,90,120 and 180

图3 仿真小鼠重建图像ROI 区域: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL。从上到下依次为投影幅数为60、90、120 和180 的ROI 区域Fig.3 The reconstruction results of the simulated mouse model: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL.From top to bottom,the ROI regions are 60,90,120 and 180

图4 绘制了五种算法在重建图像过程中在四个不同角度下NRMSE、PSNR 和SSIM 评价指标的收敛曲线图。从图4 可以看出，OGS-HL 算法对于图像重建的效果对比其他算法有所提升，其在180 角度下的提升较小，在120 角度下重建的效果提升，在60 和90 角度下的重建效果提升较为明显。仿真试验结果表明，OGS-HL 算法在稀疏角度重建的条件下会保持良好的优势，能够提升重建图像的质量。

4.3 临床小鼠实验

为了进一步说明OGS-HL 算法的有效性，本节实验使用MARS (Medipix All Resolution System)微型CT 上采集的来自真实临床前小鼠的投影数据。重建的CT 图像覆盖面积为18.41 mm×18.41 mm，分辨率为512×512。投影数据采用等距扇形束几何采集，X射线源到旋转中心的距离为158 mm，到探测器的距离为255 mm。在360°范围内进行等间隔采样，设置的采样间隔角度分别为6°、4°、3°和2°，得到的稀疏投影角度个数为60、90、120 和180。

在OGS-TV 算法中，经过多次实验，参数为µ=50，λ=2×10-4δ=0.1，q=1，K=3。在OGS-HL算法中，参数为 µ=50，λ=2×10-4，δ=0.1，q=0.8，K=3。

图5 展示了6 种不同算法下不同角度的临床真实小鼠重建图像，图6 展现了临床真实小鼠骨的部分ROI 区域。从图5 可以看出，每个算法重建的效果与以上仿真小鼠重建实验的效果基本相同，在四个不同的角度下FBP 算法和SART 算法重建过程产生了大量的伪影和噪声。从图6 的ROI 区域看出，TV 算法下小鼠骨切面的边缘出现伪影，且含有少量的噪声，TGV 算法在60、90 角度和120 角度下的噪声去除并不明显，OGS-TV 算法也含有少量的噪声，OGS-HL对于骨切面周围噪声的消除和边缘的保护有良好的效果。从实验结果和ROI 区域可以看出，OGS-HL 算法在稀疏角度条件下的重建效果优于其他算法。

图5 临床小鼠模型重建结果: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL。从上到下依次为投影幅数为60、90、120 和180 的重建图像Fig.5 The reconstruction results of the clinical mouse model: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL.From top to bottom,the reconstructed images are from 60,90,120 and 180

图6 临床小鼠模型重建图像ROI 区域: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL。从上到下依次为投影幅数为60、90、120 和180 的ROI 区域Fig.6 The reconstruction results of the clinical mouse model: (a) FBP;(b) SART;(c) TV;(d) TGV;(e) OGS-TV;(f) OGS-HL.From top to bottom,the ROI regions are 60,90,120,and 180

5 结论

为了进一步提升稀疏角度CT 图像的重建质量，本文将OGS 和HL 应用于稀疏角度CT 图重建，得到OGS-HL 稀疏角度CT 迭代重建算法。相比于传统的TV 正则项，OGS 正则项考虑每个像素点梯度重叠组的信息，度量图像梯度稀疏性，克服了图像重建过程中产生的阶梯伪影，取得了不错的重建效果。图像的梯度服从重尾分布，HL 先验可以很好地近似图像梯度的重尾分布，达到抑制噪声的效果。通过以上仿真小鼠和临床小鼠两个实验结果表明，OGS-HL 算法对于稀疏角度CT 图像重建的图像质量提升明显，具有较好的鲁棒性。本文是基于模型的稀疏CT 重建方法，侧重于物理模型和数学算法，而基于深度学习的稀疏CT 重建方法则注重于数据驱动和自主学习能力，更侧重于大量的CT 数据作为训练样本，通过学习数据中的信息来实现重建。由于本文进行的实验没有大量的训练数据集，故采用基于模型的算法来实现重建。在参数选择上，正则项参数 λ和惩罚参数 δ的选择对重建图像的效果极为重要，但是，本文这两个参数的选取为经验选取，导致进行实验时耗费的时间较长，之后的研究会进一步改进模型，使算法模型更加灵活简便。