通过双重低秩张量分解结合自适应变换，实现低剂量光谱CT的低复杂度重建

《Medical Image Analysis》：Low-complexity reconstruction of low-dose spectral CT via double low-rank tensor factorization with adaptive transforms

　　刘春燕|吴同乐|王宏|胡殿林|张斌|王建军
中国重庆西南大学数学与统计学院，邮编400715

**摘要**
光谱计算机断层扫描（CT）是一种医学成像技术，通过同步收集多能量X射线数据并分析不同能量范围内物质的吸收特性，为医生提供更详细和准确的成像信息。然而，每个能量区间内的光子计数不足可能会导致重建图像中出现显著噪声，尤其是在低剂量条件下。本文提出了一种基于空间因素去噪的稀疏视图光谱CT重建方法，该方法采用双重低秩张量分解框架。该方法通过充分利用全局光谱低秩性、群体稀疏性和非局部自相似性，有效抑制混合噪声并降低计算复杂度。具体而言，本文使用双重低秩张量分解来模拟光谱CT图像中固有的全局低秩结构。为了充分利用空间因子中的非局部自相似性，我们在空间和光谱维度上构建了三阶张量块，并通过引入加权张量核范数约束来增强其低秩性。同时，为了促进结构信息的共享，我们对空间因子的空间梯度施加了加权群体稀疏性约束。因此，所提出的模型能够通过巧妙结合全局、非局部和局部先验信息，真实地表征光谱CT图像的潜在结构。此外，我们使用近端交替最小化算法来求解所提出的模型，并为其收敛性提供了理论保证。广泛的数值模拟和临床患者数据实验表明，与现有的流行算法相比，所提出的方法在降低计算复杂性和提高光谱CT图像质量方面具有显著优势。

**引言**
光谱计算机断层扫描（CT）利用光子计数探测器（PCD）从狭窄的能量区间获取多通道投影，从而生成由多个能量区间组成的光谱CT图像。该技术利用不同能量X射线对物质的吸收特性差异，实现比传统CT技术更精细的组织表征（Boll等人，2010年）、材料分解（Johnson等人，2007年）和病变检测（Lv等人，2012年）。然而，将X射线谱划分为多个能量区间会显著减少每个能量区间内的光子数量，导致该区间仅捕获总入射光子的一小部分。这种光子数量不足直接导致单个能量区间的噪声显著增加。此外，当前的CT成像技术通常涉及高辐射剂量，而这种高剂量暴露可能会增加患者患癌症的风险（Brenner和Hall，2007年；HalL和Brenner，2014年）。因此，开发一种高效的低剂量重建方法以辅助疾病诊断在实际临床应用中具有重要意义。

**经典迭代重建方法**，如代数重建技术（ART）（Gordon等人，1970年）和同时代数重建技术（SART）（Andersen和Kak，1984年），可以在一定程度上减少图像噪声和伪影。然而，它们在没有先验信息的情况下难以获得满意的CT图像重建效果。利用图像先验信息是去噪光谱CT图像的关键。为此，研究人员提出了一系列统计迭代重建（SIR）（Elbakri和Fessler，2002年）算法，这些算法整合了先验知识以改善图像重建并减少辐射剂量。SIR方法的核心优势在于它可以充分利用图像中的固有先验信息，从而有效指导重建过程。目前，流行的先验可以分为空间相关性（Tian等人，2011年；Chang等人，2017年）和光谱相关性（He等人，2022年；Wu等人，2024年）。经典的空间相关性正则化方法，如总变分（TV）（Peng等人，2020年；Chu等人，2013年）正则化，是维护图像处理中局部空间分段平滑性的常用工具。因此，许多研究人员引入了各种类型的TV正则化方法来更好地抑制噪声并保持边缘信息。例如，Rigie和La Rivire（2015年）引入了总核变分（TNV）方法来改善光谱CT图像的重建效果。此外，Zeng等人（2016年）提出了一种基于结构张量总变分的加权最小二乘法（PWLS-STV），通过引入STV正则化来有效消除TV中出现的斑块伪影。为了提高收敛速度并降低计算复杂度，Fang等人（2015年）提出了一种基于张量TV（TTV）模型的低剂量CT灌注（CTP）方法，该方法结合了血管结构的空间相关性和血液流动的时间连续性。在Li等人（2019年）的研究中，Li等人在低秩张量Tucker分解框架中引入了时空总变分正则化（LRTD-STTV）来改善脑灌注CT的重建效果。

尽管基于低秩矩阵/张量建模的现有TV正则化方法已经取得了相对满意的恢复效果，但它们在表示光谱CT图像的空间相关性方面能力不足，难以实现最佳去噪性能。相比之下，基于非局部自相似性（NSS）（Wu等人，2019b）的正则化方法可以更有效地表示图像的空间相关性。与传统的局部空间相关性方法（如TV（Yuan，2016年）不同，NSS正则化方法可以通过利用不同空间位置相似块的信息，更有效地捕捉复杂结构和纹理特征，从而在图像重建和去噪任务中表现出更好的性能。因此，研究人员使用张量稀疏编码或4D协作过滤技术从空间-光谱域中提取和聚合相关特征，并提出了TDL（Zhang等人，2017年）、l0TDL（Wu等人，2018a）和SSCMF（Wu等人，2018b）等重建算法来有效表示NSS先验。然而，这些方法未能对非局部相似块进行聚类、向量化及堆叠，导致重建图像的清晰度和细节保留不足（Xie等人，2016年；Wu等人，2019a）。为此，研究人员通过聚类、向量化和堆叠相似块构建了基本单元，并采用低秩矩阵/张量建模方法更精确地捕捉图像结构，从而提高重建质量。例如，基于块的低秩惩罚（Kim等人，2015年）、非局部低秩和稀疏矩阵分解（Niu等人，2018年）以及非局部张量分解（Xia等人，2019年；Wang等人，2021年；X. Chen等人，2022年；Liu等人，2024年）。此外，为了充分利用内在的低秩性，Lin等人（2025年）提出了一种低秩四元数重建算法，该算法将非局部相似块构建为四元数矩阵，并引入了调整后的加权Schatten-p范数来增强低秩表示。然而，现有方法要么忽略了能量区间间的群体稀疏性，要么在非局部张量构建中产生高计算成本，要么缺乏统一的空间-光谱先验，导致在低剂量条件下细节模糊和伪影出现。因此，上述方法主要存在以下两个局限性：

1. **稀疏性表示不足**：光谱CT的不同能量区间捕获来自同一横截面的信息，因此每个能量区间图像中的平滑区域通常位于相同位置。这意味着当空间差异图像沿光谱维度展开成矩阵时，其大多数行向量的?2范数为零或接近零，这被称为群体稀疏性。图1(b)展示了沿光谱维度展开的胸部图像的空间差异图像及其对应的直方图分布。具体来说，从图1(b1)和(b2)可以看出，每行的元素几乎呈现全零或全非零的分布模式，非零行的数量远少于零行的数量，这表明这些元素具有共享的群体稀疏结构。图1(b3)和(b4)分别展示了图1(b1)和(b2)中每行的?2范数值的直方图分布。显然，大多数?2范数值为零或接近零。然而，传统的TV方法仅通过?1范数来约束空间差异图像的稀疏性，未能充分利用差异图像中的群体稀疏结构，从而无法充分探索不同能量区间之间的共同特征。
2. **复杂性增加和缺乏空间相关性**：随着能量区间数量的增加，构建的非局部三阶张量的第三维度也随之增加。所有非局部三阶张量的第三维度的增加显著提高了后续低秩张量近似的计算复杂度。现有的低秩张量分解（LRTF）框架可以有效降低计算复杂度，因此在图像噪声去除（Zheng等人，2020年；Zhang等人，2022年；Y. Chen等人，2022年；Zhang等人，2024年）和图像重建（Fan等人，2015年；B. Liu等人，2023年）等应用中得到广泛应用。然而，上述方法要么仅采用单一的空间-光谱分解结构，要么未能充分考虑空间相关性（例如NSS），因此难以有效恢复重建图像中的细微细节和结构信息。

为了解决上述问题，本文提出了一种非局部双重LRTF（NDLRTF）方法，通过光谱-空间协作正则化来实现混合噪声抑制（见图2）。如图1所示，我们观察到胸部图像的三个关键特性：低秩性、群体稀疏性和NSS。所提出的方法将这三个特性整合到一个统一框架中。该方法不仅通过张量分解将高维图像投影到低维光谱子空间中以减少数据冗余，还利用空间先验约束来保留结构细节和光谱信息，从而实现光谱CT图像特性的协作优化。具体而言：
(1) 首先，基于滤波反投影（FBP）获得初始图像，并在其DLRTF框架内对其全局光谱低秩性进行建模。即，光谱CT图像X通过模式-3张量-矩阵乘积分解为空间因子和低维光谱基（X=B×3L1和X=C×3L2，×3表示模式-3张量-矩阵乘积，即X(i1,i2,j)=∑i3B(i1,i2,i3)L1(j,i3)）。对于空间因子B，提取图像块，并将不同能量区间内空间位置相似的块堆叠起来构建非局部三阶张量，其低秩性受到加权张量核范数（WTNN）的约束。
(2) 其次，我们对空间因子C施加加权?2,1-范数正则化，以强制差异图像中不同能量区间间的群体稀疏性。同时，自适应学习正交基L1和L2以捕获全局低秩结构，共同增强低秩和群体稀疏特征的提取。
(3) 最后，我们将全局光谱低秩性、空间NSS和群体稀疏性整合到一个统一的优化框架中，构建NDLRTF模型。开发了近端交替最小化（Attouch等人，2013年）（PAM）算法来高效求解NDLRTF模型，并证明了其全局收敛到临界点的理论性。在模拟和临床数据上的比较实验表明，所提出的方法在混合噪声抑制、空间结构恢复和光谱特征保留方面显著优于几种最先进算法。

本文的其余部分组织如下：第2节介绍NDLRTF模型的符号、预备知识和优化过程，并提供理论收敛性保证。第3节介绍算法收敛性分析、参数选择、复杂性分析和实验结果。最后，第4节对本文进行总结。

**符号和预备知识**
为简洁起见，我们在表1中解释了主要符号和定义，其中一些来自Lu等人（2020年）。基于t-积，我们引入了张量奇异值分解（T-SVD），它是加权张量核范数（WTNN）的基础。

**定理1 T-SVD（Lu等人，2020年）**
设X∈Rn1×n2×n3，则X的T-SVD表示为：X=U?S?VT，
其中U∈Rn1×n1×n3和V∈Rn2×n2×n3是正交的，S∈Rn1×n2×n3是一个F-对角张量。T-SVD的详细步骤如下：首先进行FFT。

**实验配置**
在实验中，我们使用模拟数据和临床数据评估所提出的NDLRTF方法：小鼠胸部模体、胸部/腹部患者数据和水模体。选择SART（Andersen和Kak，1984年）、WTTV（Chu等人，2013年）、ASSIST（Xia等人，2019年）和NWTCTV-WTTV（Liu等人，2024年）作为基准方法与所提出的方法进行比较。这些比较方法的参数配置基于它们各自提供的建议。

**结论和讨论**
在本文中，我们提出了一种基于DLRTF框架的NDLRTF算法，用于去除高斯和泊松混合噪声。该方法使用DLRTF框架来表征全局光谱低秩特性，同时在空间因素中引入了非局部和群稀疏性约束。这种设计结合了光谱CT图像的噪声敏感性和低秩特性，以及光谱差异图像的共享稀疏性，从而实现了更强的噪声消除和更好的细节恢复。

**收敛性分析**
在这一部分，我们给出了所提出的NDLRTF算法的全局收敛定理。定理2的证明受到了Attouch等人（2013年）中定理6.2的启发。我们通过证明L({M,B,C,X,L1,L2})满足以下三个条件，证明了序列{Mt,Bt,Ct,Xt,L1t,L2t}收敛到L({M,B,C,X,L1,L2})的临界点：
1. L({Mt,Bt,Ct,Xt,L1t,L2t})在每个点上都满足K?性质；
2. 序列{Mt,Bt,Ct,Xt,L1t,L2t}满足充分下降条件；
3. ...

**参数分析**
为了选择胸腔图像的最佳参数值，我们绘制了在不同参数设置下65个投影视角和1×10^6光子数下的PSNR和SSIM变化曲线，如图18所示。我们清楚地观察到，当α=1×10^-5且β=1时，PSNR和SSIM获得了最满意的结果。图18(e)–(n)展示了平衡因子μ、ρ、ρ1、λ和η对PSNR和SSIM的影响。从图中可以看出，当μ=5×10^-4、ρ=5×10^-3、ρ1=5×10^3、λ=5×10^-3和η=1×10^-6时，PSNR和SSIM表现最佳。

**正交性和初始化分析**
基于方程(26)和(28)中L1和L2的解，正交性约束（L1?L1=I和L2?L2=I）在每次迭代中通过正交Procrustes解被理论性地强制执行。为了验证这种强制执行的数值稳定性，我们使用1×10^6光子和65个投影视角下的胸腔图像对NDLRTF方法的优化过程进行了正交性分析。每个因素的正交性误差...

**在1×10^4光子下的临床数据研究**
在本小节中，我们进一步评估了所提出的NDLRTF方法在低光子条件下的性能，使用胸腔和腹部CT图像进行对比实验。首先，我们使用NDLRTF方法和其他竞争算法重建光谱CT图像，并选择了一个细节丰富的小区域进行评估。然后，通过差异图分析了不同方法在噪声抑制和纹理细节恢复方面的性能。

**挑战性案例研究**
为了评估NDLRTF方法在高噪声水平（5×10^3、1×10^4、1×10^5光子）和65个投影视角下的鲁棒性，我们使用腹部图像进行了一个挑战性案例研究。腹部图像包含高对比度的复杂组织，如软组织、骨骼和脂肪，这使得重建变得非常困难。有限的投影数量和高噪声进一步加剧了这一挑战。如图27和表14所示，第8能级的定性和定量结果...

**高分辨率重建分析**
为了进一步评估所提出方法在视觉质量方面的性能，我们在1×10^6光子和65个投影视角下进行了额外的重建实验，分辨率为512×512。重建图像的大小为512×512×8，图像阵列覆盖了40×40 mm2的区域。X射线源到旋转中心的距离设置为50 mm，而X射线源到探测器的距离为100 mm。探测器由1024个元素组成，总长度为90...

**CRediT作者贡献声明**
刘春燕：撰写——原始草稿、可视化、方法论、调查、形式分析、数据管理。
吴同乐：验证、方法论、形式分析。
王宏：撰写——审阅与编辑、验证、形式分析。
胡殿林：撰写——审阅与编辑、验证。
张斌：撰写——审阅与编辑、监督。
王建军：撰写——审阅与编辑、监督、项目管理、资金获取。

**利益冲突声明**
作者声明他们没有已知的竞争性财务利益或个人关系可能影响本文报告的工作。

**致谢**
本工作得到了重庆市自然科学基金（CSTB2023NSCQ-LZX0044, cstc2021ycjh-bgzxm0015）、国家重点研发计划（2021YFB3101500）和中国自然科学基金（12071380）的支持。