基于运动参数解耦与运动约束驱动的锥形束计算机断层扫描刚性运动伪影校正方法

Motion parameter decoupling and motion constraint-driven optimization for correcting rigid motion artifacts in cone-beam computed tomography

摘要

目的 为了解决锥形束计算机断层扫描（CBCT）长时间扫描过程中因患者运动和机架抖动导致的刚性运动伪影问题，我们提出一种基于运动参数解耦与运动约束驱动的CBCT刚性运动校正方法。 方法 采用基于3D-2D刚性配准的运动估计框架，将刚性运动参数划分为层间运动和层内运动，并设计分步优化顺序解耦参数间的相互干扰。针对迭代过程中伪影表征由多轮廓重叠向边缘模糊的动态演变，构建渐进式代价函数，实现优化目标从投影数据一致性约束向结构细节恢复的适应性过渡。为解决3D-2D配准的多解性问题及迭代过程中可能引入的图像空间位置偏移，引入运动估计约束机制消除全局偏差，从而增强算法收敛性。 结果 所提算法能够准确估计刚性运动轨迹并运动补偿修复CBCT图像。在头部仿真数据中，本文方法在3种运动程度下均取得最优量化指标。相较于次优方法，峰值信噪比（PSNR）平均提升约2.1%，结构相似性（SSIM）平均提升约2.9%，均方根误差（RMSE）平均降低约6.5%。在膝关节仿真与真实猪肢干数据上的实验结果进一步验证了算法的有效性与泛化能力。 结论 本文提出的刚性运动伪影校正算法在估计运动轨迹与抑制图像伪影方面具有良好性能，为临床CBCT成像中刚性运动伪影的抑制提供了一种可行且稳定的解决方案。

Abstract

Objective To solve the problem of rigid motion artifacts caused by patient movement and gantry vibration during long-duration cone-beam computed tomography (CBCT) using geometric parameter decoupling and motion-constrained optimization. Methods Using a motion estimation framework based on 3D-2D rigid registration, the rigid motion parameters were categorized into out-of-plane and in-plane motions, and a stepwise optimization sequence was designed to decouple the mutual interference among the parameters. In response to dynamic evolution of artifact characteristics from multi-contour overlap to edge blurring during the iterative process, a progressive cost function was formulated to facilitate an adaptive transition of the optimization objective from projection data consistency constraints to structural detail recovery. To address the issue of multi-solution in 3D-2D registration and the potential spatial misalignment introduced during iteration, a motion estimation constraint mechanism was incorporated to eliminate global bias, thereby enhancing the algorithm's convergence. Results The proposed algorithm accurately estimated rigid motion trajectories and restored the images via motion compensation. On head simulation data, this model achieved optimal quantitative metrics across 3 motion levels, improved PSNR by 2.1% and SSIM by 2.9%, and reduced RMSE by 6.5% compared to the suboptimal methods. Additional validation on knee simulation and real porcine data further demonstrated its efficacy and generalization capability. Conclusion The proposed rigid motion artifact correction algorithm demonstrates good performance in estimating motion trajectories and compensating for image artifacts, thus providing a viable and robust solution for suppressing rigid motion artifacts in clinical CBCT imaging.

Graphical abstract

关键词

刚性运动伪影校正 / CBCT / 3D-2D刚性配准 / 参数解耦 / 运动约束

Key words

rigid motion artifact correction / cone-beam computed tomography / 3D-2D rigid registration / parameter decoupling / motion constraint

锥形束计算机断层扫描（CBCT）在牙科、颌面外科、耳鼻喉科、放射治疗以及图像引导介入等诸多医学专业中具有广泛应用^［1］。它能够提供扫描部位的高分辨率三维图像^［2］，为疾病诊断和治疗规划提供重要的影像学依据。然而，CBCT的机架旋转速度与平板探测器采集帧率受硬件性能制约，导致扫描时间延长，因此不可避免地增加了扫描过程中患者运动和机架抖动的风险^［3］。经典的CT重建算法基于静态扫描目标和理想扫描轨迹的假设。患者运动和机架振动破坏了预设的成像几何关系，导致测量投影信号在CBCT图像重建过程中沿着错误轨迹反向投影，从而在重建图像上产生刚性运动伪影。这些伪影在图像上表现为条纹伪影、多重轮廓伪影、结构失真等，严重干扰临床诊断和治疗^［4］。此外，为获取诊断质量的图像，常需重复扫描，这将使患者暴露于额外辐射剂量下，增加辐射累积效应及潜在长期健康风险。因此，需要一种可靠的方法来防止运动伪影。

为了降低运动发生的可能性，常使用头带、牙托、固定咬合块等约束装置辅助患者进行CBCT成像^［5，6］。但这些装置无法完全限制患者运动，残余的轻微运动仍可能在重建图像中产生伪影。在紧急头部外伤等情况下，使用约束装置可能不切实际。对于某些患者，如儿科患者，可采用全身麻醉或镇静，但这可能增加扫描期间及之后的不良反应风险^［7］。另一方面，CBCT制造商探究更快的成像方法，如提高机架旋转速度、相控阵列CBCT等^{［8， 9］}，然而，这些方法受限于当前硬件和技术能力。

为了克服CBCT成像中的刚性运动伪影，发展出多种刚性伪影校正方法，主要分为两类。其一是通过增强投影数据一致性或优化图像质量度量来消除伪影，逼近干净图像。例如，Lee团队通过对齐圆形轨迹的镜像投影来修复运动引起的投影错位^［10］；Berger等^［11］基于弦图域的傅里叶一致性约束抑制伪影；Zhang等^［12］则将运动退化图像建模为干净图像与模糊核的卷积，借助盲去卷积算法联合估计二者。随着深度学习技术的广泛应用，众多数据驱动与模型融合的运动修复方法也被提出。例如，Lee团队^［13］基于U-WGAN对抗网络学习运动伪影特征，以修复运动受损图像。Ali等^［14］利用深度学习提取投影数据中无运动部分并结合加权重建方法实现图像复原。然而，运动修复方法虽然能够在一定程度上减轻图像中的刚性运动伪影，但未能从根本上解决投影数据与成像几何不匹配的问题，导致其在恢复结构细节方面存在不足，甚至可能引入二次失真。

第二类方法通过估计患者运动轨迹参数实现运动补偿，分为有基准标记方法和与无基准标记方法。有基准标记方法通常借助外部光学系统或基准标记物来追踪患者运动^［15-18］。然而，基于跟踪器的运动测量引入了额外硬件，而且在机架旋转过程中易受视线遮挡干扰，其几何精度也难以满足亚毫米级校正需求^［19］。鉴于外部跟踪器的局限性，研究者们转而通过跟踪投影内部不透明的基准标记物来估计运动。虽然这种方法更具优势，但放置基准标记物可能干扰临床工作流程，且标记物附着结构可能存在弹性变形，无法准确反映患者运动^［20］。无基准标记方法主要包括投影数据一致性方法、3D-2D配准方法及自校准方法。投影数据一致性方法调整成像几何参数以最大化投影数据间的一致性度量^［21-24］。3D-2D配准方法模拟CT前向模型生成虚拟投影并与实际采集的 2D 投影进行细致比对，寻找相似度最大化所对应的成像几何^［25-27］。自校准方法以量化形式评价图像优劣，优化成像几何参数，促使重建图像质量提升^［28-31］。然而，上述方法多基于随机优化策略，存在以下固有缺陷：首先，其收敛性严重依赖参数设置，缺乏严格的收敛保证；其次，在高维运动参数空间中计算效率显著降低，难以满足临床实时性要求；此外，不同扫描案例需重新调整优化参数，严重影响了方法的鲁棒性和可重复性。这些局限性使得基于随机优化的运动估计方法在精度、效率和稳定性方面难以兼顾，制约了其在临床实践中的推广应用。

鉴于此，本文提出了一种基于3D-2D配准的运动估计和补偿方法。为避免陷入局部最小值，提出基于参数解耦的分层确定性搜索策略进行参数优化，而非随机优化。该策略将刚性运动参数划分为层间运动（引起相邻层间信息混叠）和层内运动（导致同层图像内容混叠），并精心设计分层优化顺序降低参数间的耦合干扰；同时，针对不同阶段刚性运动伪影表征定制了渐进式代价函数，以实现优化目标从投影数据一致性约束向结构细节恢复的适应性过渡。此外，鉴于3D-2D刚性配准的多解性以及运动导致的图像空间位置偏差，我们引入一种运动估计约束机制，以增强算法的收敛性。基于3D-2D配准的运动估计和补偿通常需要频繁执行3D重建与前投影，计算成本高昂。针对此问题，利用刚性运动特性，结合GPU并行优化，提出了一个快速高效的解决方案。

1 资料和方法

1.1 CBCT刚性运动补偿

CBCT三维图像重建的本质是将多角度采集的X射线投影数据重建成表征组织衰减系数的空间分布。图像重建的准确性取决于对X射线源、探测器与患者之间空间几何关系的精确估计。理想情况下，患者保持静止且CBCT扫描系统沿预设轨迹稳定运行，预估计的成像几何参数能够准确描述三维空间点与二维探测器坐标之间的映射关系，并有多种解析和代数重建算法稳定求解重建问题。然而，在实际临床扫描中，患者不自主运动和CBCT机架抖动造成预设成像几何和实际投影几何之间出现差异。为了解决这一问题，运动轨迹参数估计方法通常采用两阶段处理流程：首先基于投影与图像数据估计刚性运动轨迹，随后基于估计结果修正成像几何并完成运动补偿重建。处理流程（图1）。

CBCT扫描过程中患者的刚性运动及机架抖动均可等价地建模为对扫描系统施加刚性变换。因此，我们寻求等效刚性变换来校正成像几何，随后利用校正后的成像几何结合测量投影数据进行运动补偿重建。成像几何由X射线源位置、平板探测器中心位置及探测器行列方向单位向量共同描述，并在CBCT扫描前已知。获得等效刚性变换后，成像几何通过以下公式更新：

\begin{array}{c}{\mathit{S}}_j^{\text{'}}=\#{\mathit{T}}_j\bullet {\mathit{S}}_j\end{array}（1）

\begin{array}{c}{\mathit{D}}_j^{\text{'}}=\#{\mathit{T}}_j\bullet {\mathit{D}}_j\end{array}（2）

\begin{array}{c}{\widehat{\mathit{u}}}_j^{\text{'}}=\#{\mathit{R}}_j\bullet {\widehat{\mathit{u}}}_j\end{array}（3）

\begin{array}{c}{\widehat{\mathit{v}}}_j^{\text{'}}=\#{\mathit{R}}_j\bullet {\widehat{\mathit{v}}}_j\end{array}（4）

其中，{\mathit{S}}_j\mathrm{、}{\mathit{D}}_j分别表示第j个投影角度下X射线源与探测器中心的原始位置，{\widehat{\mathit{u}}}_j\mathrm{、}{\widehat{\mathit{v}}}_j为探测器行列方向的单位向量。{\mathit{T}}_j为第j个投影角度下成像系统的等效刚性变换矩阵，{\mathit{R}}_j为其旋转分量。{\mathit{S}}_j^{\mathrm{\text{'}}}\mathrm{、}{\mathit{D}}_j^{\mathrm{\text{'}}}\mathrm{、}{\widehat{\mathit{u}}}_j^{\mathrm{\text{'}}}\mathrm{、}{\widehat{\mathit{v}}}_j^{\mathrm{\text{'}}}分别对应校正后的X射线源位置、探测器中心位置以及探测器行列方向单位向量。

1.2 CBCT刚性运动估计

1.2.1 运动建模

实验采用的世界坐标系(x,y,z)如图2所示。CBCT绕z轴旋转，且初始时刻X射线源位于x轴正方向。为便于分析，本文引入了旋转坐标系，其原点与世界坐标系原点重合，且每一投影角度下：其u轴平行于探测器像素(\mathrm{0},\mathrm{0})到(\mathrm{0},\mathrm{1})的矢量方向，v轴平行于像素(\mathrm{0},\mathrm{0})到(\mathrm{1},\mathrm{0})的矢量方向，而w轴与探测器法向量方向一致。

我们采用运动参数向量\mathit{p}=[{\mathit{p}}_{\mathrm{1}},\dots ,{\mathit{p}}_j,\dots ,{\mathit{p}}_N{]}^T描述CBCT扫描过程中各投影角度对应的等效刚性变换，其中N为总投影数。第j个投影角度刚性运动参数定义为{\mathit{p}}_j=[t_{jw},t_{ju},t_{jv},{\phi }_{jw},{\phi }_{ju},{\phi }_{jv}{]}^T。此处t_{jw},t_{ju},t_{jv}分别表示沿旋转坐标系w,u,v轴的平移量，{\phi }_{jw},{\phi }_{ju},{\phi }_{jv}为绕对应轴的旋转角度。刚性运动参数{\mathit{p}}_j通过以下公式转化为等效刚性变换矩阵{\mathit{T}}_j：

\begin{array}{c}{\mathit{T}}_j=\left[\begin{array}{cc}{\mathit{R}}_j& {\mathit{t}}_j\\ \mathrm{0}& \mathrm{1}\end{array}\right]=\left[\begin{array}{cc}{\mathit{R}}_{\theta }& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right]\left[\begin{array}{cc}{\mathit{R}}_{{\mathit{p}}_j}& {\mathit{t}}_{{\mathit{p}}_j}\\ \mathrm{0}& \mathrm{1}\end{array}\right]{\left[\begin{array}{cc}{\mathit{R}}_{\theta }& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\end{array}\right]}^T\end{array}（5）

{\mathit{t}}_{{\mathit{p}}_j}={\left[t_{jw},t_{ju},t_{jv}\right]}^T（6）

\begin{array}{c}{\mathit{R}}_{{\mathit{p}}_j}=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{jv}& -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{jv}& \mathrm{0}\\ \mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{jv}& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{jv}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{ju}& \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{ju}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{ju}& \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{ju}\end{array}\right]\left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{jw}& -\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{jw}\\ \mathrm{0}& \mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{jw}& \mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{jw}\end{array}\right]\end{array}（8）

其中，{\mathit{R}}_j和{\mathit{t}}_j分别{\mathit{T}}_j的旋转分量和平移分量；{\mathit{R}}_{\theta }是CBCT机架固有的旋转，用于将旋转坐标系中运动变换至世界坐标系；{\mathit{R}}_{{\mathit{p}}_j}为由3个旋转角度{\phi }_w,{\phi }_u,{\phi }_v构成的旋转矩阵，旋转顺序依次为绕w,u,v轴；{\mathit{t}}_{{\mathit{p}}_j}为平移向量。

1.2.2 渐进式代价函数设计

我们采用基于3D-2D配准的运动估计方法，通过模拟CBCT前向投影过程，以最小化重投影与测量投影的差异为目标估计运动参数，进而通过迭代优化实现成像几何校正与图像恢复。在迭代优化过程中，成像几何误差与CBCT刚性运动伪影表征呈现显著相关性。当成像几何误差较大时，重建图像出现明显的多轮廓重叠与解剖结构扭曲；随着误差逐渐减小，伪影表征逐步转变为局部结构失真，主要表现为边缘模糊与拖影等现象。针对上述伪影表征的演变规律，我们设计了渐进式代价函数。该复合函数采用动态加权策略：在迭代初期（{\lambda }_i=\mathrm{1}）侧重基于投影数据一致性的全局配准；在迭代后期（{\lambda }_i=\mathrm{0}）则强化边缘结构特征的局部恢复。构建的优化目标函数如下：

\underset{\mathit{p},\mathit{x}}{\mathrm{A}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left({\lambda }_i{L}_{\mathrm{1}}+\left(\mathrm{1}-{\lambda }_i\right)L_{\mathrm{2}}\right)\left(\mathrm{8}\right)

L_{\mathrm{1}}={‖\mathit{y}-A_{\mathit{p}}\mathit{x}‖}_{\mathrm{1}}\left(\mathrm{9}\right)

L_{\mathrm{2}}={‖\mathrm{L}\mathrm{O}\mathrm{G}(\mathit{y})-\mathrm{L}\mathrm{O}\mathrm{G}(A_{\mathit{p}}\mathit{x})‖}_{\mathrm{1}}\left(\mathrm{10}\right)

其中A_{\mathit{p}}表示经运动参数\mathit{p}更新的成像几何模型，\mathit{x}为重建图像，\mathit{y}为测量投影数据，{\lambda }_i为随迭代次数i变化的超参数，\mathrm{L}\mathrm{O}\mathrm{G}为高斯-拉普拉斯算子。超参数{\lambda }_i初始设为1，并在优化过程中根据运动参数估计的收敛状态自适应调整。具体而言，定义第i次迭代与第i-\mathrm{1}次迭代之间各刚性运动参数的变化量为（∆t_{jw},∆t_{ju},∆t_{jv},∆{\phi }_{jw},∆{\phi }_{ju},∆{\phi }_{jv}）。考虑到t_{jw}（沿中心射线方向的平移）仅引起投影的微弱缩放效应，难以准确估计，故其变化未纳入切换依据。当所有投影角度下其余参数变化量的绝对值之和满足

（11）

时，判定运动估计趋于稳定，此时残余误差较小，刚性运动伪影表征为边缘模糊，因此将{\lambda }_i由1切换为0，转而强化边缘结构细节的恢复。

1.2.3 参数求解

采用全参数空间网格搜索策略同时优化6个运动参数需要极大的计算资源与内存开销，在实际应用中难以实现。逐步优化虽然可以有效减少资源消耗，但面临运动参数间耦合干扰的固有问题。基于刚性运动的物理特性，可将参数划分为层间运动({\phi }_w,{\phi }_u,t_v)和层内运动({\phi }_v,t_w,t_u)，并据此实施分步校正。首先估计层间运动参数，即先估计并校正v轴方向平移t_v，随后同步优化w轴与u轴的旋转参数{\phi }_w与{\phi }_u。对于剩余的层内运动参数，其在等中心平面附近主要表现为层内伪影。基于这一特性，可通过3D-2D配准该区域内有限探测器行的投影数据来实现层内运动参数的快速同步校正。各校正步骤均采用网格搜索优化相应的运动参数（依次为t_v； {\phi }_w，{\phi }_u； {\phi }_v， t_w， t_u），其搜索范围与步长依据大数据统计先验设定，从而在保证全局搜索能力前提下有效约束解空间，显著降低计算负担。方法整体流程如图3所示。

为了实现高效可靠的运动校正，需克服网格搜索策略带来的巨大计算需求和内存压力。基于刚体运动的整体性特征，我们将3D-2D配准范围约束于等中心层临近区域，以降低内存占用，同时避免锥角伪影干扰。此外，我们进一步结合稀疏角度配准与GPU并行计算。首先按稀疏倍数对投影角度均匀采样，构建稀疏子集，并在该子集上估计各角度的刚性运动参数。由于各角度的参数估计相互独立，可基于GPU并行框架执行多角度的刚性运动参数搜索。随后，采用一维线性插值结合窗口平滑滤波，对稀疏子集估计结果进行拟合，从而重构出全投影角度下的连续运动轨迹。

1.3 运动估计约束

在基于3D-2D配准的CBCT运动估计中，由于恢复的是患者与成像系统之间的相对几何，其解空间表现出刚性等价特性，即对估计运动轨迹施加任意刚性变换并相应调整成像区域，利用相同投影数据可以重建出完全一致的CBCT图像。这一特性导致迭代优化过程中可能引入患者与成像系统的同步空间偏移。随着迭代次数增加，图像内容可能逐渐偏移出预设的成像区域，造成运动参数估计误差累积与图像质量持续退化，最终形成难以收敛的恶性循环。为解决该问题，需对运动轨迹施加有效约束。鉴于对运动轨迹施加约束在抑制空间偏移的同时也会引入额外的估计误差，我们仅对影响最显著的v轴方向平移进行约束。在每次运动估计之后，所有投影角度的t_v参数减去首帧投影对应的t_v参数值，以消除空间位置偏差。

t_{jv}^{\text{'}}=t_{jv}-t_{\mathrm{0}v}（12）

1.4 实验设计

1.4.1 数据准备及实验设置

仿真研究基于商用牙科CBCT扫描仪（Matrix5000，Fussen）采集的无伪影的头部CBCT图像进行刚性运动伪影仿真。采用标准圆周扫描轨迹，利用8控制节点的Akima样条模型随机生成刚性扰动(t_w,t_u,t_v,{\phi }_w,{\phi }_u,{\phi }_v)，据此更新X射线源位置、探测器中心位置以及探测器行列方向单位向量并对原始CBCT图像进行前投影。获得的投影数据添加泊松噪声（每条射线入射光子数为1×10⁵）后作为测量投影并采用标准圆周成像几何重建出CBCT刚性运动伪影图像。CBCT成像参数设置为：源-等中心距离425 mm，源-探测器距离635 mm，平板探测器尺寸400×500像素，像素间距0.6 mm，且系统理想对齐（平板探测器无偏置和倾斜）。为了评估对于不同CBCT刚性运动伪影表征的校正性能，我们基于Akima样条模型，在3种幅度范围内随机采样控制点并仿真相应的刚性运动伪影：［-0.5，0.5］、［-1.0，1.0］和［-3.0，3.0］（\mathrm{m}\mathrm{m}/°），分别对应轻度（Minor artifact）、中度（Moderate artifact）和重度（Severe artifact）。此外，仿真研究额外采用1例膝关节CBCT扫描数据，在保持相同运动模型与成像参数的条件下，验证所提算法对于不同部位刚性运动伪影校正的泛化性能。

真实数据使用西安交通大学CBCT实验平台采集，扫描对象为猪肢干。实验平台及样本设置如图4所示，其中A为西安交通大学CBCT实验平台，B为放置在CBCT旋转台上带控制云台的猪肢干样本。CBCT设备配置1536×1536像素的平板探测器，像素尺寸为0.1 mm，几何标定参数为：源-等中心距379.146 mm，源-探测器距528.377 mm，探测器左偏0.802mm、上偏48.379 mm，xyz三轴倾角分别为0.217°、-0.281°、-0.630°。在80 kVp/8 mA条件下，通过云台施加刚性运动，采集576帧投影数据。并通过对比标准静止扫描的重建结果，验证算法实效。

本文所提算法在配备Intel Xeon E5-2686 v4处理器与NVIDIA RTX 4090显卡的计算平台上运行。代价函数切换阈值\tau设置为5。运动参数搜索范围设置为［-5，5］（mm）和［-5，5］（°），搜索步长为0.2。该搜索设置依据CBCT扫描过程中患者不自主运动的统计数据^［1］，即绝大多数头部刚性运动的平移和旋转幅度均在±5 mm/±5°范围内，同时0.2的步长在保证配准精度的前提下兼顾计算效率。稀疏角度倍数为5。经过10次迭代优化后，选取代价函数最小值对应的结果作为最终输出。

1.4.2 对比方法及评价指标

本文选取5种方法作为对照：包括Sun等 ^［25］与Chen等 ^［27］提出的运动估计方法，以及Restormer^［32］、LIT-Former^［33］和TT-UNet^［34］3种基于深度学习的图像恢复方法。其中，对于深度学习方法，我们基于原始头部CBCT图像，通过刚性运动伪影仿真构建包含36 864组配对图像的数据集用于模型训练。本文所提方法与其他对比方法均在18例独立的运动仿真患者数据上进行性能测试。本文采用3项标准化指标进行定量评价：结构相似度（SSIM）衡量图像结构特征保持能力，峰值信噪比（PSNR）表征信号重建质量，均方根误差（RMSE）量化像素级重建精度。各指标数学定义如下：

SSIM\left(I_{correct},I_{gt}\right)=\frac{\left(\mathrm{2}{\mu }_{correct}{\mu }_{gt}+C_{\mathrm{1}}\right)\left(\mathrm{2}{\sigma }_{correct\_gt}+C_{\mathrm{2}}\right)}{\left({\mu }_{correct}^{\mathrm{2}}+{\mu }_{gt}^{\mathrm{2}}+C_{\mathrm{1}}\right)\left({\sigma }_{correct}^{\mathrm{2}}+{\sigma }_{gt}^{\mathrm{2}}+C_{\mathrm{2}}\right)}（13）

其中，I_{correct},I_{gt} 分别为校正图像和真值图像，{\mu }_{correct},{\mu }_{gt} 为对应均值，{\sigma }_{correct},{\sigma }_{gt} 为对应标准差，{\sigma }_{correct\_gt} 为校正图像I_{correct}和真值图像I_{gt}间的协方差，常数C_{\mathrm{1}}=(K_{\mathrm{1}}{L)}^{\mathrm{2}}、C_{\mathrm{2}}=(K_{\mathrm{2}}{L)}^{\mathrm{2}}取 K_{\mathrm{1}}=\mathrm{0.01}、K_{\mathrm{2}}=\mathrm{0.03}、L=\mathrm{255}。

\begin{array}{l}\begin{array}{c}\begin{array}{l}PSNR\left(I_{correct},I_{gt}\right)=\\ \mathrm{10}lo{g}_{\mathrm{10}}\frac{MA{X}_{I_{correct}}^{\mathrm{2}}}{\frac{\mathrm{1}}{MNT}{\displaystyle \sum _{i=\mathrm{1}}^M}{\displaystyle \sum _{j=\mathrm{1}}^N}{\sum }_{k=\mathrm{1}}^T{\left[I_{correct}\left(i,j,k\right)-I_{gt}\left(i,j,k\right)\right]}^{\mathrm{2}}}\end{array}\end{array}\\ \end{array}（14）

\begin{array}{c}\begin{array}{l}RMSE\left(I_{correct},I_{gt}\right)=\\ \sqrt[]{\frac{\mathrm{1}}{MNT}{\displaystyle \sum _{i=\mathrm{1}}^M}{\displaystyle \sum _{j=\mathrm{1}}^N}{\displaystyle \sum _k^T}{\left[I_{correct}\left(i,j,k\right)-I_{gt}\left(i,j,k\right)\right]}^{\mathrm{2}}}\end{array}\end{array}（15）

其中，M,N,T为CBCT图像的3个维度大小，（i,j,k）为体素索引，MA{X}_{I_{correct}} 是校正图像的最大体素值。考虑到刚性运动导致图像空间位置偏差，所有校正图像在计算定量指标之前均需与真值图像进行刚性配准，以确保指标计算的准确性。

2 结果

2.1 仿真数据实验结果

图5~7分别从横截面、矢状面与冠状面3个视角，对比展示不同方法在头部仿真数据上的复原效果。各图首行为真值图像及3种运动程度下的受损图像，其余各行为对比方法及本文方法的复原结果；图中关键细节由黄色虚线框标注出，并在右下角附有对应的局部放大图。图像显示窗均为［-1000，2500］HU。可以观察到，随着运动程度加剧，刚性运动伪影依次呈现为边缘模糊、结构变形乃至多重轮廓重叠，图像失真状况显著加剧。在严重刚性运动伪影情况下，Chen和Sun的方法校正性能下降，复原图像中仍存在边缘模糊或结构失真。网络恢复方法借助大数据先验在一定程度上恢复轻微伪影下的牙釉质边界，但对结构失真类伪影修复效果有限。而本文方法对3种运动程度下的受损图像均有良好的复原效果，解剖结构细节被成功修复，如高亮牙釉质层。图8展示了Chen、Sun及本文方法对图5~7中严重刚性运动伪影例子的等效刚性运动轨迹估计结果。t_w偏移量仅引起投影放大比例的微小变化，其误差可以忽略不计。而对于其它参数的估计，本文方法的估计结果与真值最为接近。表1展示了不同运动伪影校正算法在头部仿真数据上的量化指标结果。本文方法在PSNR、RMSE和SSIM三项指标上均取得最优。具体而言，在严重伪影情况下，本文方法的PSNR（30.57）显著高于对比方法，RMSE（7.55）则最低，同时SSIM（0.8018）也保持领先，这表明本文方法在提升图像信噪比、降低误差以及保持结构相似性方面均优于其他对比算法。

图9展示了运动估计方法在膝关节仿真数据上的复原效果。图像按3行排列，分别对应横断面、矢状面及冠状面视图；首列为真值图像，第2列为运动受损图像，第3~5列为运动估计方法的复原结果。图像显示窗均为［-1000，1200］HU。在膝关节运动受损图像中可见明显的多重轮廓叠加现象，具体表现为骨边界重复显现。Chen、Sun及本文方法均能有效校正由刚性运动导致的轮廓混叠与模糊伪影。其中本文方法在整体复原质量上表现最优，在肌肉、韧带等软组织的纹理还原方面展现出更高的保真度，重建结果与真值解剖结构的一致性最高。不同运动估计方法在膝关节仿真数据上的定量比较结果（表2）。本文方法在PSNR（35.7197）、RMSE（4.1740）和SSIM（0.9531）上均优于Chen和Sun的方法，进一步验证了所提算法在不同解剖部位数据上的有效性和泛化能力。

2.2 真实数据实验结果

图10展示了不同运动估计方法在真实猪肢干数据上的复原结果。图像按2行排列，分别显示横截面和矢状面视图。首列为无运动干扰下的真值图像；第2列为为通过云台施加刚性运动所产生的受损图像；第3~5列为各运动估计方法对应的复原结果。图像显示窗均为［-1200，1200］HU。CBCT扫描过程中引入的刚性运动导致图像出现明显的解剖结构混叠：边缘模糊、结构失真、条形伪影。Chen和Sun的方法虽能恢复样本的大致轮廓，但在细节重建方面存在不足。相较之下，本文方法实现了解剖结构的完整重建，尤其在骨缝等精细结构处恢复准确，复原结果与无运动真值图像高度保真。在真实猪肢干数据上的定量评估结果如表3所示。本文方法取得了最高的PSNR（35.5427）和SSIM（0.8949），以及最低的RMSE（4.2599），说明其在真实运动场景下仍能有效恢复图像质量，复原结果与无运动真值图像最为接近，展现了良好的临床应用鲁棒性。

2.3 运动估计约束的必要性及有效性验证

图11对比了不同运动估计约束策略下的图像复原效果。从左至右依次为：真值图像、运动受损图像、全参数约束复原结果、无约束复原结果，以及仅对t_v参数施加约束的复原结果。各图像右下角放大框清晰呈现了黄色虚线框标注区域的局部细节。结果显示，全参数约束及无约束方案虽能恢复组织结构的基本轮廓，但在牙齿凹槽形态及牙根管结构等解剖特征处仍存在明显模糊。仅对t_v参数施加约束的方案在保持牙齿细微结构完整性方面表现最优。图12对比了不同约束策略下的刚性运动轨迹估计结果。无约束方案与仅对t_v参数约束的方案的估计轨迹均接近真值，但后者吻合度更高，这表明对t_v参数施加约束具有积极作用。

2.4 渐进式代价函数的必要性及有效性验证

图13对比了不同代价函数对图像复原效果的影响。从左至右依次展示了真值图像、运动受损图像、基于L1代价函数的复原结果，以及采用渐进式代价函数的复原结果。结果可见，复原图像均消除大部分运动伪影，但在细节复原质量上存在差异。L1代价函数未能完整恢复高亮牙釉质层的边缘结构及牙根管的形态特征，而渐进式代价函数则在保留这些解剖结构方面表现出更优的性能，牙釉质边缘清晰，牙根管形态也得到了准确的恢复。该结果表明渐进式代价函数设计能够有效适应不同迭代阶段中伪影特征的变化，从而提升图像复原的精度与结构保真度。

3 讨论

本文针对CBCT扫描过程中因患者运动和机架抖动导致的刚性运动伪影问题，提出了一种基于参数解耦的分层确定性搜索策略优化的CBCT刚性运动估计与补偿算法。该算法通过精心设计的优化顺序解耦运动参数间的相互干扰，并针对不同阶段刚性运动伪影表征定制了渐进式代价函数，有效解决了传统随机优化方法易陷入局部最小值的问题。

实验结果表明，基于深度学习的运动修复方法在轻微运动伪影场景下能够借助数据先验恢复部分结构边界，但其校正能力存在一定局限。这类方法主要在图像域进行操作，通过学习到的图像先验来修复被伪影破坏的结构信息。然而，由于未能解决投影数据与成像几何失配这一根本问题，当面对中、重度运动导致的结构失真与多重轮廓重叠等复杂伪影时，其校正效果会受到限制。投影数据沿错误的反投影路径重建所引入的混叠伪影具有病态特性，仅从图像外观上进行修复难以完全恢复真实的解剖结构，甚至可能引入非真实的细节或因过度平滑损失关键纹理，导致复原图像保真度不足。

在仿真数据实验中，本文方法在头部和膝关节数据上均表现出优异的性能。随着运动程度加剧，刚性运动伪影从边缘模糊逐步演变为结构变形乃至多重轮廓重叠，本文方法在三种运动程度下均能有效恢复解剖结构细节，特别是在牙釉质层等高亮结构的恢复方面表现突出。定量分析结果表明，本文方法在PSNR、RMSE和SSIM三项指标上均优于对比方法。运动轨迹估计结果进一步验证了本文方法在参数估计精度方面的优势。

在真实猪肢干数据实验中，本文方法展现出了良好的临床应用潜力。相较于Chen和Sun的方法，本文方法在骨缝等精细结构的恢复方面更为准确，重建结果与无运动真值图像高度一致。

本文方法的主要优势在于：首先，提出基于参数解耦的分层确定性搜索策略避免了传统随机优化方法对参数设置的敏感性，提供了更为稳定的收敛保证；其次，通过分阶段优化和参数解耦设计，结合GPU并行计算与稀疏角度配准策略，显著提升了运算效率。在实际测试中，单例数据的完整运动估计与补偿过程仅需数分钟即可完成，具备良好的临床适用性；再者，渐进式代价函数的设计契合了刚性运动伪影表征的演变规律，增强了算法的适应性与结构保真度。

需要指出的是，运动估计约束策略的有效性高度依赖于具体的算法框架。 本研究在验证约束机制时发现，全参数约束方案在部分投影角度下的估计结果与真值存在较大偏差，这可能源于其引入的额外误差超过了约束本身的校正能力。因此，在本算法中，仅对影响最显著的t_v平移参数施加约束展现了最佳性能。这一现象揭示了一个更深层次的原理，即约束本质上是参数校正精度与引入误差风险之间的权衡。不同的运动估计算法因其模型和优化策略的差异，可能导致误差传播特性截然不同。因此，约束策略的最优形式不能一概而论，而需根据具体算法的复原结果与误差特性进行定制。

然而，本研究仍存在一些局限性，有待在后续工作中深入探究。首先，当前运动参数搜索采用固定步长，虽在典型临床场景下兼顾了计算效率与估计精度，但对亚毫米级微小运动（如生理性震颤）的分辨能力有限。未来将引入分层自适应搜索机制，即先在稀疏角度子集上以粗步长快速定位运动大致范围，再在完整角度集上围绕初步估计构建局部搜索窗口，并动态细化步长，以提升微小运动的估计精度。其次，本文方法基于全局刚性运动假设，难以适用于存在局部非刚性形变的场景（如下颌张闭口、软组织滑移等）。在此类情况下，单一刚性变换无法准确刻画不同解剖区域的运动差异，可能导致校正不充分。为此，未来将探索“全局–局部”混合运动建模框架，在全局刚性补偿后，进一步对残差投影数据建模局部运动，以拓展方法对非刚性运动的适应能力。此外，在临床实际场景中，金属种植体、低剂量扫描及稀疏角度成像等因素均会显著加剧投影数据的噪声与不完备性。此类条件下，重建图像严重退化，导致3D-2D配准中的一致性度量失真，进而影响运动参数估计的准确性。后续工作将融合图像域恢复网络，在3D-2D配准前对配准图像进行结构增强与伪影抑制，以提供更可靠的图像先验，从而增强算法在复杂临床场景下的鲁棒性。

综上所述，本文提出的基于参数解耦的分层确定性搜索策略优化的CBCT刚性运动伪影校正算法能够准确估计刚性运动轨迹并通过运动补偿有效恢复CBCT图像质量，同时具备较高的计算效率，为解决临床扫描中的运动伪影问题提供了一种可靠且实用的技术方案。

参考文献

[1] Spin-Neto R, Wenzel A. Patient movement and motion artefacts in cone beam computed tomography of the dentomaxillofacial region: a systematic literature review[J]. Oral Surg Oral Med Oral Pathol Oral Radiol, 2016, 121(4): 425-33. doi：10.1016/j.oooo.2015.11.019

[2] Loubele M, Bogaerts R, Van Dijck E, et al. Comparison between effective radiation dose of CBCT and MSCT scanners for dento-maxillofacial applications[J]. Eur J Radiol, 2009, 71(3): 461-8. doi：10.1016/j.ejrad.2008.06.002

[3] Nemtoi A, Czink C, Haba D, et al. Cone beam CT: a current overview of devices[J]. Dentomaxillofac Radiol, 2013, 42(8): 20120443. doi：10.1259/dmfr.20120443

[4] Spin-Neto R, Costa C, Salgado DM, et al. Patient movement characteristics and the impact on CBCT image quality and interpretability[J]. Dentomaxillofac Radiol, 2018, 47(1): 20170216. doi：10.1259/dmfr.20170216

[5] Nardi C, Taliani GG, Castellani A, et al. Repetition of examination due to motion artifacts in horizontal cone beam CT: comparison among three different kinds of head support[J]. J Int Soc Prev Community Dent, 2017, 7(4): 208-13. doi：10.4103/jispcd.jispcd_17_17

[6] Hanzelka T, Dusek J, Ocasek F, et al. Movement of the patient and the cone beam computed tomography scanner: objectives and possible solutions[J]. Oral Surg Oral Med Oral Pathol Oral Radiol, 2013, 116(6): 769-73. doi：10.1016/j.oooo.2013.08.010

[7] Malviya S, Voepel-Lewis T, Eldevik OP, et al. Sedation and general anaesthesia in children undergoing MRI and CT: adverse events and outcomes[J]. Br J Anaesth, 2000, 84(6): 743-8. doi：10.1093/oxfordjournals.bja.a013586

[8] Beeres M, Wichmann JL, Paul J, et al. CT chest and gantry rotation time: does the rotation time influence image quality?[J]. Acta Radiol, 2015, 56(8): 950-4. doi：10.1177/0284185114544242

[9] Berland LL, Smith JK. Multidetector-array CT: once again, technology creates new opportunities[J]. Radiology, 1998, 209(2): 327-9. doi：10.1148/radiology.209.2.9807555

[10] Hernandez D, Eldib ME, Hegazy MAA, et al. A head motion estimation algorithm for motion artifact correction in dental CT imaging[J]. Phys Med Biol, 2018, 63(6): 065014. doi：10.1088/1361-6560/aab17e

[11] Berger M, Maier A, Xia Y, et al. Motion compensated fan-beam CT by enforcing fourier properties of the sinogram[C]//2014 Third International Conference on Image Formation in X-ray Computed Tomography. Salt Lake City, Utah, USA. 2014: 329-32.

[12] Zhang Y, Zhang LY. A rigid motion artifact reduction method for CT based on blind deconvolution[J]. Algorithms, 2019, 12(8): 155. doi：10.3390/a12080155

[13] Hu ZL, Jiang CH, Zhang QY, et al. Wasserstein generative adversarial networks for motion artifact removal in dental CT imaging[C]//Medical Imaging 2019: Physics of Medical Imaging. February 16-21, 2019. San Diego, USA. SPIE, 2019: 113. doi：10.1117/12.2511818

[14] Ali ASRA, Fusiello A, Landi C, et al. Motion artifacts detection in short-scan dental CBCT reconstructions[EB/OL]. 2023: arXiv: 2304.10154. doi：10.1109/nssmicrtsd49126.2023.10337961

[15] Weisenberger AG, Gleason SS, Goddard J, et al. A restraint-free small animal SPECT imaging system with motion tracking[J]. IEEE Trans Nucl Sci, 2005, 52(3): 638-44. doi：10.1109/tns.2005.851399

[16] Herzog H, Tellmann L, Fulton R, et al. Motion artifact reduction on parametric PET images of neuroreceptor binding[J]. J Nucl Med, 2005, 46(6): 1059-65.

[17] Kim JH, Nuyts J, Kyme A, et al. A rigid motion correction method for helical computed tomography (CT)[J]. Phys Med Biol, 2015, 60(5): 2047-73. doi：10.1088/0031-9155/60/5/2047

[18] Kyme AZ, Se S, Meikle SR, et al. Markerless motion estimation for motion-compensated clinical brain imaging[J]. Phys Med Biol, 2018, 63(10): 105018. doi：10.1088/1361-6560/aabd48

[19] Reaungamornrat S, Otake Y, Uneri A, et al. Tracker-on-C for cone-beam CT-guided surgery: evaluation of geometric accuracy and clinical applications[J]. Med Imag 2012 Image Guid Proced Rob Interv Model, 2012, 8316: 831609. doi：10.1117/12.911454

[20] Ouadah S, Jacobson M, Stayman JW, et al. Correction of patient motion in cone-beam CT using 3D-2D registration[J]. Phys Med Biol, 2017, 62(23): 8813-31. doi：10.1088/1361-6560/aa9254

[21] Yu HY, Wang G. Data consistency based rigid motion artifact reduction in fan-beam CT[J]. IEEE Trans Med Imaging, 2007, 26(2): 249-60. doi：10.1109/tmi.2006.889717

[22] Berger M, Xia Y, Aichinger W, et al. Motion compensation for cone-beam CT using Fourier consistency conditions[J]. Phys Med Biol, 2017, 62(17): 7181-215. doi：10.1088/1361-6560/aa8129

[23] Aichert A, Berger M, Wang J, et al. Epipolar consistency in transmission imaging[J]. IEEE Trans Med Imaging, 2015, 34(11): 2205-19. doi：10.1109/tmi.2015.2426417

[24] Preuhs A, Maier A, Manhart M, et al. Double your views–exploiting symmetry in transmission imaging[M]//Medical Image Computing and Computer Assisted Intervention – MICCAI 2018. Cham: Springer International Publishing, 2018: 356-64. doi：10.1007/978-3-030-00928-1_41

[25] Sun T, Kim JH, Fulton R, et al. An iterative projection-based motion estimation and compensation scheme for head X-ray CT[J]. Med Phys, 2016, 43(10): 5705. doi：10.1118/1.4963218

[26] Berger M, Müller K, Aichert A, et al. Marker-free motion correction in weight-bearing cone-beam CT of the knee joint[J]. Med Phys, 2016, 43(3): 1235-48. doi：10.1118/1.4941012

[27] Chen MY, He P, Feng P, et al. General rigid motion correction for computed tomography imaging based on locally linear embedding[J]. Opt Eng, 2018, 57(2): 1. doi：10.1117/1.oe.57.2.023102

[28] Huang H, Siewerdsen JH, Zbijewski W, et al. Reference-free learning-based similarity metric for motion compensation in cone-beam CT[J]. Phys Med Biol, 2022, 67(12): 125020. doi：10.1088/1361-6560/ac749a

[29] Sisniega A, Stayman JW, Yorkston J, et al. Motion compensation in extremity cone-beam CT using a penalized image sharpness criterion[J]. Phys Med Biol, 2017, 62(9): 3712-34. doi：10.1088/1361-6560/aa6869

[30] Kingston A, Sakellariou A, Varslot T, et al. Reliable automatic alignment of tomographic projection data by passive auto-focus[J]. Med Phys, 2011, 38(9): 4934-45. doi：10.1118/1.3609096

[31] Maur S, Stsepankou D, Hesser J. Auto-calibration by locally consistent contours for dental CBCT[J]. Phys Med Biol, 2018, 63(21): 215018. doi：10.1088/1361-6560/aae66d

[32] Zamir SW, Arora A, Khan S, et al. Restormer: efficient transformer for high-resolution image restoration[C]//2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). June 18-24, 2022, New Orleans, LA, USA. IEEE, 2022: 5718-29. doi：10.1109/cvpr52688.2022.00564

[33] Chen ZH, Niu C, Gao Q, et al. LIT-former: linking in-plane and through-plane transformers for simultaneous CT image denoising and deblurring[J]. IEEE Trans Med Imaging, 2024, 43(5): 1880-94. doi：10.1109/tmi.2024.3351723

[34] Deng ZH, Zhang WK, Chen KL, et al. TT U-Net: temporal transformer U-Net for motion artifact reduction using PAD (pseudo all-phase clinical-dataset) in cardiac CT[J]. IEEE Trans Med Imaging, 2023, 42(12): 3805-16. doi：10.1109/tmi.2023.3310933

基金资助

国家重点研发计划(2024YFC2417801)

国家自然科学基金(62571422)

广州市科技计划项目(202206010148)

RIGHTS & PERMISSIONS

版权所有©《南方医科大学学报》编辑部2021