李晓晴，刘锡明，苗积臣，吴志芳

(清华大学 核能与新能源技术研究院 核检测技术北京市重点实验室，北京 100084)

计算机断层成像技术(CT)是一种先进的非接触无损检测技术，能清晰、准确地对物体成像，可反映物体的内部结构、密度、缺陷位置和形状等特性[1]。经过了近50年的发展，CT技术已广泛应用于医学诊断、工业探伤、港口安检等诸多领域。随着平板探测器技术的不断成熟，基于平板探测器的锥束CT成为了重要的工业无损检测手段[2]。

影响锥束CT成像质量的因素较多，其中CT系统的几何参数精确获取程度是影响CT图像重建质量的关键因素。按照是否需要特定模体，现有的参数校正方法可分为无模体法[3-6]与模体法[7-10]两大类。用于获取几何参数的标定模体的结构材质各异，典型的标定模体有西北工业大学使用的细琴弦和带豁口的钢尺[11]、Ford与Cho等[4,12]使用的将多个小钢球等角排列而成的双钢珠平行圆环。相较于通过反复迭代重建图像以获取几何参数的无模体法，模体法具有参数获取速度较快、结果不受初值影响等优点，使用更为广泛，但模体法存在模体不易制作或集成度差，使用时可能需要更换模体、操作繁琐等问题。

为解决上述问题，本文设计一种易于制作、操作简便、集成化的模体，提出一种参数快速获取方法：平行双丝法。本文首先介绍该方法的原理，利用几何投影关系推导锥束CT几何参数的计算公式，然后通过仿真实验评估方法的准确性与鲁棒性，最后使用真实模体，在高精度锥束CT系统平台上进行实验验证。

1 参数定义

对于锥束CT，本文定义射线源到探测器的距离为DSD、射线源到旋转中心的距离为DSO，旋转中心在探测器上的投影位置坐标为(u0,v0)，探测器绕轴u=u0的旋转角度记为σ，探测器绕轴v=v0的旋转角度记为ζ，探测器在u-v平面内绕点(u0,v0)旋转的角度记为η。图1为锥束CT几何参数的定义。这些参数对重建图像的质量影响程度不同。根据文献[13-14]，影响最大的两个参数分别为u0和η，本文将重点关注这两个参数。默认探测器的旋转角度σ和ζ为0°。

2 平行双丝法原理

本文提出的方法使用的模体是两根平行金属丝，故将方法命名为平行双丝法。两根平行丝记为l1、l2，设计时l1、l2的距离d已知。

a——CT系统DSD、DSO参数定义；b——探测器绕轴u=u0的旋转角度σ；c——探测器绕轴v=v0的旋转角度ζ；d——探测器绕点(u0，v0)的旋转角度η图1 锥束CT几何参数定义Fig.1 Geometric parameter of cone-beam CT

平行双丝法的基本思想为：将模体放在工件转台上，保持双丝与转台平面垂直，转台旋转1周，获取模体360°投影图像，对投影图像进行数据处理，计算出DSD、DSO、u0、η4个参数。关键步骤推导如下。

2.1 u0、DSD计算方法

操作时将模体置于转台上，两丝垂直于转台表面，使丝与旋转轴平行。由几何知识可知，在旋转1周的过程中，有且仅有两个位置使得两平行丝的投影重合，且这两个位置的投影射线与中心线夹角相等，记为α。图2为平行丝投影重合位置示意图，图2中选取过射线中心线且垂直于旋转轴的平面作为截面。A1、A2分别为第1根细丝l1的两个重合位置，B1、B2分别为第2根细丝l2的两个重合位置，两个虚线圆分别为l1、l2旋转1周的轨迹，O为旋转中心，D为射线源与旋转中心连线在探测器上的投影位置，S为射线源的位置。用A0代表细丝l1的初始位置，转台从初始位置旋转到双丝重合的第1个位置旋转角度为θ1、旋转到第2个重合位置的角度为θ2。

图2 平行丝投影重合位置示意图Fig.2 Coincidence position of projection for parallel wires

在理想位置下，SD⊥DP1、SD⊥DP2、DP1=DP2。因而射线源中心线在探测器上的投影横坐标u0等于细丝l1在探测器的投影横坐标uP1与细丝l2在探测器的投影横坐标uP2的平均值：

(1)

射线源到探测器的距离为：

(2)

2.2 DSO计算方法

找到双丝投影重合位置后，寻找平行丝所在平面与中心线垂直时的投影位置。由几何知识可知，旋转1周的过程中，这样的位置有且仅有两个，如图3所示。A3、A4分别为双丝所在平面与中心线垂直时细丝l1的两个位置，B3、B4分别为垂直时细丝l2的两个位置。两虚线圆分别为l1、l2旋转1周的轨迹。用A0代表细丝l1的初始位置，记转台从初始位置旋转到平行丝所在平面与中心线垂直的第1个位置时角度为θ3，旋转到第2个位置时角度为θ4。

图3 平行丝连线垂直于中心射线时投影示意图Fig.3 Projection diagram of perpendicular to central ray position for parallel wires

在旋转1周的过程中，两平行丝距离是不变的，A1B1=A2B2=A3B3=A4B4=d，旋转中心到两丝连线的距离r不变。由等比定理可得：

(3)

2.3 不同投影位置的角度关系

通过对图像中两根丝投影重合度的判断可得到θ1、θ2，需要通过θ1、θ2计算出α、θ3、θ4，下面给出角度关系的推导公式。为了使推导最为简易，选取旋转中心到双丝连线的垂足位置进行计算。图4为角度关系示意图，其中圆为旋转中心到双丝连线垂足的旋转轨迹，C0为初始时垂足的位置，C1、C2为双丝投影重合时垂足的位置，对应的转台旋转角度为θ1、θ2，C3、C4为双丝平面与中心线垂直时垂足的位置，对应的转台旋转角度为θ3、θ4；圆的两条切线表示双丝投影重合时穿过双丝的射线，与中心线的夹角为α。

图4 转台旋转角度关系Fig.4 Rotation angle relation of turntable diagram

依据初始位置的不同，角度关系有4种情况：

(4)

2.4 η的计算方法

探测器绕点(u0,v0)旋转角度为η时，投影情况如图5所示，图5中虚线表示理想探测器位置，实线表示旋转后探测器位置。可知，旋转后两丝的投影平行，与探测器v′轴的夹角等于探测器绕点(u0,v0)的旋转角度η。

图5 探测器绕点(u0,v0)旋转前后投影对比示意图Fig.5 Projection before and after detector rotate about point (u0,v0)

实验获得双丝投影后，计算投影图中两丝与探测器v轴的夹角，即得到探测器绕点(u0,v0)的旋转角度η。按照式(5)进行坐标变换，变换后的几何参数和理想情况下是一致的。

(5)

考虑实际情况下丝的放置位置与转台旋转轴可能不完全平行，为了消除这种影响，计算双丝旋转1周获得的所有投影图像的η值，并计算这些η值的平均值，平均值即为系统的η值。

在实际测量校正中，由于旋转中心投影横坐标v0对重建图像质量影响很小，所以不对其进行精确计算，简单采用探测器坐标v轴中值作为v0值。

3 基于仿真方法的误差分析

为验证在探测器绕轴u=u0旋转，或绕轴v=v0旋转情况时方法的鲁棒性，利用Matlab进行仿真模拟计算。模拟计算的参数为：DSD=1 150 mm，DSO=900 mm，d=150 mm，旋转中心到两丝连线的距离为100 mm。

3.1 σ误差分析

图6a为探测器绕轴u=u0旋转时细丝的投影，虚线表示转前探测器的位置，实线表示转后探测器的位置；图6b为Matlab模拟出的探测器无旋转和绕轴u=u0旋转角度σ为5°时的双丝投影。可看出，两丝投影依然平行，但投影在探测器上的长度、两投影间距离均有改变。

图6 探测器绕轴u=u0旋转前后平行丝投影对比Fig.6 Projection before and after detector rotate around axis u=u0

本文计算了探测器旋转角度σ在0.5°～10°范围内变化时，对几何参数u0、DSD、DSO的影响，表1列出几何参数的计算结果与精确值的相对误差。

3.2 ζ误差分析

图7a为探测器绕v=v0轴旋转时细丝的投影，虚线表示旋转前探测器位置，实线表示旋转后探测器位置；图7b为计算机模拟出的探测器无旋转和绕轴v=v0旋转角度ζ为5°时双丝的投影对比。可看出，两丝的投影不再平行，会因探测器的倾斜而产生一定的夹角。

表1 不同σ下几何参数相对误差Table 1 Relative error of geometric parameter in different σ

图7 探测器绕v=v0轴旋转前后平行丝投影对比Fig.7 Projection before and after detector rotate around axis v=v0

本文计算了探测器旋转角度ζ在0.5°～10°范围内变化时，对几何参数u0、DSD、DSO的影响，表2列出了几何参数计算结果与精确值的相对误差。

表2 不同ζ下几何参数的相对误差Table 2 Relative error of geometric parameter in different ζ

由上述模拟结果可知，本方法受探测器绕u=u0轴旋转角度的影响略高于受绕v=v0轴旋转角度的影响，但影响均非常小，在肉眼可见的角度偏差内对参数计算的影响可忽略不计，证明该方法鲁棒性较好。

4 实验验证

4.1 实验方案

为检验方法的实际效果，实际制作了模体并在高精度工业CT系统上进行实验验证。选取0.3 mm钽丝作为模体细丝，外壳采用有机玻璃，两丝距离d为100 mm。系统参数及实验参数为：X射线源能量，320 keV；探测器像素尺寸，0.2 mm；探测器像素数，2 048×2 048；扫描投影数，2 000。

4.2 操作步骤

实验操作步骤为：1) 将模体垂直放置于转台平面上，转台旋转360°进行投影，得到各角度的投影值；2) 根据投影图中双丝与v轴的夹角，计算探测器绕点(u0,v0)的旋转角度η；3) 找到双丝投影重合时的投影图，记录转台旋转到两个重合位置的旋转角度θ1、θ2，记录投影图中双丝投影的横坐标uP1、uP2；4) 根据式(4)计算两平行丝所在平面垂直于中心线时的转台旋转角度α、θ3、θ4；5) 找出转台旋转角度为θ3的投影图，记录投影图中两丝的距离P3P4；6) 找出转台旋转角度为θ4的投影图，记录投影图中两丝的距离P4P5；7) 根据式(1)～(3)，计算出几何参数DSD、DSO和旋转中心投影横坐标u0。

4.3 实验结果

由于探测器v轴像素为2 048，故将v0设为1 024。表3列出几何参数设计值与实验值的对比，可看出，实验值与设计值吻合较好。

表3 几何参数设计值与实验值对比Table 3 Comparison between designed and experimental geometric parameters

图8为实验值重建图形，可看出，图像质量清晰度较好，说明平行双丝法获取的参数满足实际需求。

图8 实验值重建图形Fig.8 Reconstruction image of experiment result

5 结论

针对现有锥束CT参数获取模体法存在的模体不集成、制作要求较高、参数获取过程操作较为繁琐等问题，本文设计了一种操作简便的集成化参数获取方法：平行双丝法。该方法将模体旋转1周，获取360°的投影图像，即可准确获取DSD、DSO、u0、η。研究结果表明，该方法所使用的双丝模体结构简单、易于制作、操作简单、实用性强，参数获取结果精确度高、鲁棒性好，能较好地满足高精度锥束工业CT实际应用需求。