人体四肢成像是骨科疾病检查中的重要环节^[1-4]。目前临床上常用的四肢成像手段有常规X线摄影和全身CT扫描^[5]。但是，常规X线摄影难以对隐匿性骨折或者细微骨折进行精确诊断^[6]，而全身CT的空间分辨率有限^[7]，同时还面临着四肢负重成像困难、辐射剂量高等问题^[8]。

针对上述成像手段中存在的问题，四肢锥束CT应运而生。作为一种用于肢体骨科检查的低剂量3D成像技术^[9]，四肢锥束CT具有高空间分辨率、低辐射剂量和易安装等优点^[10]，并且可以提供优质的骨组织图像。因此，基于上述优点，四肢锥束CT在国外受到了极大的关注。Koskinen等^[11]人提出使用CBCT进行腕关节造影成像；Tuominen EK等人提出一种专用于下肢端骨科成像的CBCT扫描仪^[12]；Demehri S提出适于幼儿足部或脚踝异常等疾病检查的CBCT成像设备^[13]，它们存在成像空间分辨率低、构造复杂或不能实现四肢全面检查等缺点^[14]。但是在国内，四肢锥束CT的相关研究却尚无文献报道，为填补国内相关研究空白，作者所在实验室针对四肢骨骼成像高分辨特殊需求，对四肢成像专用的锥束CT（ECBCT）进行了深入研究，采用高分辨率X-线平板探测技术和GPU加速快速成像重建算法，搭建了用于四肢骨骼成像的ECBCT整机系统，该系统主要包括X射线曝光控制模块、机械控制模块、数据采集模块、数据校正模块、几何校正模块与图像重建模块。在上述模块中，几何校正模块与图像重建模块是实现优质图像重建和剂量优化的关键。在实验中，使用该系统的原型设备扫描多种体模（笔筒体模和IQ体模^[15]）以及动物肢体并进行重建，以评估该原型机系统的稳定性及成像性能。

1 系统和技术

本文所研制的ECBCT整机系统分为硬件部分和软件部分。其中硬件部分包括：X射线曝光控制模块、机械控制模块和数据采集模块；软件部分包括：数据校正模块、几何校正模块、图像重建模块。图 1展示了该系统的硬件结构和整机外观。图 2展示了系统的数据处理流程，在此流程中，由硬件部分进行数据采集，对采集的投影数据进行暗场校正和增益校正，并采用几何校正所获得的参数进行图像重建（包括解析和迭代重建）。

1.1 X射线曝光控制模块

X射线曝光控制模块包括集成的X射线管、双输出高压电源和有控制电路的灯丝电源。X射线球管采用固定阳极，其阳极倾角为5°，焦点尺寸为0.5 mm（标称值），高压电源的标称电压为50~90 KV，管电流范围为1~10 mA。系统最大功率为900 W，脉冲频率为40 HZ，扫描最大时间周期为45 s。

1.2 机械控制模块

机械控制模块主要硬件包括机架和滑块，其主要作用为控制整机系统运动。机架主要作用为保持系统采样的平稳性，其转盘放置在固定底座上（图 1），转盘的最大旋转速率为18°/s，可在20 s内实现360°扫描（本文中转盘的速率设置9°/s，实现40 s内360°扫描）。机架升降范围0~30 cm，能够实现人体下肢端脚部到膝盖不同部位的成像。

传统螺旋CT均采用滑环结构的扫描架^[16]，可实现360°扫描，但结构复杂、成本高、难以维护。鉴于此，本文提出滑块限位技术，该滑块安装在机架转盘起始处，可控制旋转一周的机架停止运动，并进行自动复位，不仅可实现360°扫描，而且具有成本低、结构简单、性能稳定等优势，滑块位于旋转盘之下，其结构如图 3所示。

1.3 数据采集模块

数据采集所使用的探测器型号为Varian PaxScan 2530DX（Varian, Palo Alto CA, USA）^[17]。该探测器采用非晶硅半导体材料，其尺寸为209.09×302.46 mm²，极限分辨率为3.60 lp/mm，可满足正常人体四肢的全覆盖、高分辨率成像。为满足临床不同需求及成像速度，探测器单元可采取组合方式为1×1，2×2，3×3和4×4，对应的像素尺寸为0.139 × 0.139 mm²，0.278 × 0.278 mm²，0.417×0.417 mm²和0.556×0.556 mm²。

从探测器接收到信号首先经过对数放大器进行对数压缩，以使后面的电路只需工作在一个窄的范围内。压缩后的信号再经过积分器计算光子总和，积分器输出信号经多路混合器变成一路，使用共同的模数变换器将连续模拟时域信号转变为离散的数字序列，最终送往计算机。

1.4 数据校正模块

在实际成像中，探测器在没有激活时也具有输出，该输出被称为探测器的暗电流^[18]；此外，探测器各个单元的响应并不完全一致^[19]。因此，需要对原始的投影数据进行暗场校正和增益校正^[20]。

（1）暗场校正：在不进行曝光的情况下采集多次投影数据并求其均值，获得暗场数据，该数据即为探测器不接收信号时的暗电流。进行暗场校正时将原始投影数据与暗场数据相减，其公式表达如下：

$ {P}'(x, y)={{P}_{0}}(x, y)-D $

(1)

其中，(x, y)为二维投影数据的索引值，P₀、P'分别为暗场校正前后投影数据$D=\frac{1}{n}\text{ }\sum\limits_{i=\text{ }1}^{n}{{}}{{D}_{i}}(x, y) $为暗场数据，n为暗场数据采集次数。

（2）增益校正：不放置物体，使用该系统在同等扫描参数下进行多次曝光，获得空扫数据，并求其均值，获得增益场数据，使用上述暗场校正方法校正该增益场数据，最后进行投影的增益校正，其公式表达为：

$ P(x, y)=\frac{{P}'(x, y)}{R-D}\times {{(R-D)}_{median}} $

(2)

其中，$R=\frac{1}{m}\sum\limits_{j=\text{ }1}^{m}{{}}{{R}_{j}}(x, y) $为增益场数据，m为空扫次数；(∙)_median为一个算子，其作用为获取括号内数据的中值；P为增益校正后的投影数据。采用上述算法，将暗场校正、增益校正集成于投影数据采集过程。数据采集的软件界面如图 4所示，其中Offest Calibration为暗场校正，Gain Calibration为增益校正。

1.5 几何校正模块

ECBCT系统组装完成之后，需要对几何参数进行标定，即系统的几何校正^[21-22]。本文采用基于钢珠点的方法进行几何校正^[23]。该方法对钢珠点进行扫描，获取其椭圆轨迹投影，并建立世界坐标系及投影坐标系，进而建立钢珠点空间坐标与其投影坐标之间的映射公式，通过求解坐标映射公式获得实际几何参数。

图 5中（A）和（B）分别展示了理想以及实际情况下对钢珠点进行扫描时的坐标映射关系。在进行图像重建时，所需要用到的几何参数有：（1）焦点到探测器的距离R；（2）焦点到旋转中心的距离R_F；（3）焦点到探测器的垂足点坐标(u₀, v₀)；（4）探测器绕着x，y，z三个轴的偏转角度分别为θ，η，φ。由此，钢珠点投影坐标与其空间坐标之间映射关系可表达为：

$ {{u}_{\alpha }}\hat{u}\text{ }+{{v}_{\alpha }}\hat{v}\text{ }+{R}'=\mu ({{r}_{\alpha }}-s) $

(3)

其中($ {{u}_{\alpha }}, {{v}_{\alpha }}$)为真实未校正探测器上投影轨迹坐标，$\tilde{x}, \tilde{z} $为x，z方向的单位向量，$ {{r}_{\alpha }}=({{x}_{\alpha }}, {{y}_{\alpha }}, z)$, α为当前投影视角，${R}'\equiv o-s=({{d}_{x}}, R+{{d}_{y}}, {{d}_{z}}) $定义了（真实）探测器相对于在焦点s处o点的位置。$ \tilde{u}\text{ }=O\tilde{x}\text{ }, \tilde{v}\text{ }=O\tilde{z}$, O为3 × 3的正交矩阵，代表探测器坐标的旋转，其公式表达为：

$ O=\left( \begin{matrix} (\text{cos}\eta \text{cos}\varphi -\text{ sin}\eta \text{sin}\theta \text{sin}\varphi )&-\text{cos}\theta \text{sin}\varphi &(-\text{cos}\varphi \text{sin}\theta -\text{ cos}\eta \text{sin}\theta \text{sin}\varphi ) \\ (\text{cos}\varphi \text{sin}\eta \text{sin}\theta +\text{ cos}\eta \text{sin}\varphi )&\text{cos}\theta \text{cos}\varphi &(\text{cos}\eta \text{cos}\varphi \text{sin}\theta -\text{ sin}\eta \text{sin}\varphi ) \\ \text{cos}\theta \text{sin}\eta &-\text{sin}\theta &\text{cos}\eta \text{cos}\theta \\ \end{matrix} \right) $

(4)

式（4）的求解采用基于傅里叶变换的解析求解方式，求解过程可参阅^[23]。

1.6 图像重建模块

本文所介绍的ECBCT系统分别采用解析算法（FDK）^[24]和迭代算法（SIRT）^[25]进行图像重建（图 6）。

（1）解析重建：FDK算法公式表述如下：

$ f(r, \varphi , z)=\frac{1}{2}\int\limits_{0}^{2\text{ }\!\!\pi\!\!\text{ }}{{}}\frac{D}{D-s}\int\limits_{-\infty }^{\infty }{{}}\frac{D}{\sqrt{{{D}^{2}}+{{l}^{2}}+\hat{z}}}g(l, \hat{z}, \beta )h({l}'-l)\text{d}l\text{d}\beta $

(5)

式中，h(l)是一维斜坡滤波器的卷积核，D是焦点到旋转中心的距离，$ g(l, \hat{z}, \beta )$是锥束投影数据，l是探测器上的线性坐标，β为投影角度，z为重建点到xy平面的距离，$ {\hat{z}}$为穿过重建点击中探测器上的像素点到l轴的距离，r为重建点到旋转中心的距离，φ为r方向与x轴方向的夹角，s为重建点到虚拟探测器的距离即$ s=r\text{sin}(\varphi -\beta )$, 及${l}'=\frac{Dr\text{cos}(\varphi -\beta )}{D-r\text{sin}(\varphi -\beta )} $。$ \frac{D}{\sqrt{{{D}^{2}}+{{l}^{2}}+\hat{z}}}$是锥形束射线入射角的余弦函数。图 5为FDK算法所用的坐标系。

（2）迭代重建：SIRT算法公式表示如下：

$ {{\mu }^{k+\text{ }1}}={{\mu }^{k}}+\alpha \frac{{{A}^{T}}(A{{\mu }^{k}}-p)}{{{A}^{T}}A} $

(6)

其中k为迭代次数，μ为重建图像，A为向前投影算子，A^T将投影图形反投到重建区域。p为投影图像。SIRT算法以FDK重建结果为初值，迭代200步。

（3）CT值校准：使用上述算法进行重建后的图像为衰减值，因此，需要对其进行CT值校准^[26]，其校准公式为：

$ CT值=\frac{({{\mu }_{物}}-{{\mu }_{水}})}{{{\mu }_{水}}}\times \text{ }1000 $

(7)

其中μ_物，μ_水分别为物体和水的衰减系数，水的衰减系数可通过对水模进行扫描重建而获得。本文所采用的水模如图 7所示，其中所填充的水为纯净水。

2 实验方法 2.1 扫描对象

采用物理体模（笔筒体模、IQ体模）和动物肢体行ECBCT成像系统性能评估实验。其中图 8为IQ体模，用于测试图像均匀性及系统的噪声功率谱；笔筒体模用于测试系统几何校正参数的准确性，动物肢体用于评估测试系统骨组织成像性能。

2.2 图像评估

通过两个方面进行定量评估。第一个方面为图像均匀性测试^[27]：根据IQ体模的重建图像，选择感兴趣区域（regions of interest, ROI），测量平均CT值及方差，根据均值、方差间的差异评估重建图像均匀性。

另一方面为噪声功率谱（noise power spectrum, NPS）^[28]，NPS以傅里叶变换为基础，将信号从时域转换至频域进行测量分析，其计算公式为：

$ NPS(u,v)=\frac{({{x}_{0}}{{y}_{0}})}{{{N}_{x}}{{N}_{y}}}|\sum\limits_{i,j}{{}}\Delta \alpha [i,j]\cdot {{e}^{-2\text{ }\!\!\pi\!\!\text{ }i\left( ui{{x}_{0}}+vj{{y}_{0}} \right)}}{{|}^{2}} $

(8)

其中，图像上的任一点(i, j)的信号强度为α[i, j], Δα[i, j]表示两幅图像做差后的噪声图像，N_x，N_y为ROI内沿着坐标轴的像素个数，x₀，y₀为采样间距，u，v为空间频率。

3 结果 3.1 系统几何校正

采用钢珠点几何校正体模，得到校正之后的几何参数:焦点到探测器的距离为743.13 mm；焦点到旋转中心得距离为538.24 mm；探测器绕着x，y，z三个轴的偏转角度分别为0.62，1.94，0.18 rad；探测器沿u，v方向分别平移6.09，39.31 mm。对笔筒进行扫描，图 9所示为笔筒扫描重建图像，其中，（A）为采用理想几何参数重建的结果，（B）为校正几何参数重建的结果。从图中可以看出，使用理想几何参数重建的图像中含有严重伪影，笔芯显示模糊；使用校正后的几何参数进行重建，图像伪影得到消除，图像细节部分能够清晰显示。

3.2 图像均匀性与噪声功率谱

通过扫描IQ体模的基底部分测量图像的均匀度和噪声功率谱，并在4个方向选取4个等轴的ROI，分别标记为ROI 1~4，其直径为FOV直径的20%，如图 10所示。ROI1，ROI2，ROI3，ROI4的平均CT值分别为1.01，1.00，0.92，0.93；方差分别为13.32，12.93，12.75，13.02，其数值波动较小，体现出图像具有较好的均匀性。

在相同扫描条件（80 KV, 9 mA）、IQ体模位置不变的情况下，根据连续两次扫描所得到重建图像，得到图像的NPS曲线，如图 11所示。图中两次扫描计算的NPS曲线具有良好的一致性，体现出系统良好的噪声稳定性。

3.3 动物肢体扫描重建

对动物肢体（猪脚）进行扫描，其重建图像的横断面、冠状面、矢状面如图 12所示，图中左、右侧图像分别为FDK及SIRT算法重建结果。为了显示系统的骨成像效果，使用阈值分割方法^[29]对各猪脚图像中的骨进行分割，并对其伪彩显示。两者的骨纹理均能够清晰显示，具有较高分辨率。从上述各图中可以看出，FDK的重建结果具有更好的图像分辨率，SIRT具有更低的噪声水平。

4 讨论

对于骨性关节炎、骨质疏松症及其并发症骨折等骨科疾病的检查，常规X线摄影，全身CT都存在不足，如常规X线摄影难以对隐匿性骨折或者细微骨折进行精确成像，全身CT骨结构成像的空间分辨率仍有限且四肢负重成像困难、辐射剂量高。鉴于此，作者所在实验室研发了针对于四肢骨成像的ECBCT系统。

本文对上述ECBCT系统进行了详细介绍。该系统主要包括X射线曝光控制模块、机械控制模块、数据采集模块、数据校正模块、图像重建模块与几何校正模块等6大模块。其中，X射线曝光控制模块、机械控制模块和数据采集模块为系统的硬件部分，其主要功能是进行目标物体的扫描及投影数据的采集。硬件部分的模块化与集成化可带来系统体积小、灵活性高、可升级性强等诸多优点。此外，机架旋转部分采用了滑块限位技术，该技术可实现全角度扫描及自动复位，同时也具有成本低廉、结构简单、实用性强等优点。

数据校正模块、几何校正模块与图像重建模块为系统的软件部分，其主要功能为进行投影数据的处理、获取准确的系统几何参数及图像重建。其中几何校正采用钢珠点坐标映射的方法，并使用傅里叶级数展开的方法求解映射方程，在保证其鲁棒性的同时，极大减轻了参数的计算量。笔筒体模重建结果表明，使用本文介绍的几何校正方法，可获得准确的几何参数。该系统提供了两种重建方法，分别为解析重建及迭代重建。动物肢体数据结果表明，两种方法均能获得高质量骨组织图像。

除却诸多优点，本文所介绍的ECBCT系统存在一些其不完善的地方，如：机架尚未添加翻转功能，在扫描人体上肢端时存在不便；未对投影数据做射束硬化效应校正及散射校正。在未来的工作，将增加整机移动及机架翻转功能，以方便上肢部位的扫描成像，此外，开发射束硬化效应校正及散射线校正算法，提高图像质量。

综上所述，本文所介绍的ECBCT系统稳定性好、分辨率高、辐射剂量低，能够对骨组织进行优质成像，为骨关节炎及骨质疏松症的诊断提供重要的临床信息。该系统的成功研制推动了国产高端骨科CT的自主化研发进程。