旷 敏，牟廉明，郑云升，叶 薇，邱利琼

（1.重庆理工大学 理学院，重庆 400054；2.内江师范学院 数学与信息科学院，四川 内江 641100；3.西华大学 理学院，成都 610039；4.顾县中学，四川 广安 638000）

1 背景

CT技术可在不破坏物体的前提下，通过断层扫描清晰直观地观测到物体的结构，其扫描部分由X射线、介质样品和一定数量的探测器组成，是一种利用X光对物体进行扫描来检测目标断层的投影技术［1－2］。Radon提出投影重建原理后，Hounsfield于1972年在英国放射学会上首次公布第1台CT系统扫描装置，从此CT技术得到迅猛发展［3－4］。而CT系统在安装时往往存在误差，导致成像质量差，因此提高CT系统几何参数的准确性和精确度的研究显得尤为重要。

目前，对CT系统几何参数标定方法的研究比较活跃。1999年，Gullberg等［5］利用质点在探测器上运动形成的投影标定出椭圆形参数，从而得到CT系统的几何参数。2000年，Noo等［6］利用最小二乘法拟合出标定模体的椭圆形投影轨迹，并与系统的几何参数组合成关系式，从而确定CT系统的几何参数。2006年，Liu等［7］比较了心图法、对角插值法、曲线拟合和几何法4种方法来测定CT系统旋转中心（COR）的精确性。同年，孙怡等［8］选择基于4个点状目标的金属板作为定标模型，利用标准模型的某一固定角度下的投影数据标定几何参数。2009年，李保磊等［9］提出利用正弦图数据的冗余性采用总和均分的方式来确定COR的位置。2013年，孟凡勇等［10］利用中心投影具有π分割数据一致性的特点，提出利用原始数据精确确定COR的方法。2018年，龙建辉等［11］基于图像轮廓曲线的提取，通过建立数学模型来标定CT系统参数。2019年，夏彬珂等［12］和李江珊等［13］采用最小二乘法拟合优化函数，并利用遗传算法求几何参数的值。以上方法侧重提高COR的精确性，十分依靠算法，并在特定条件下才可以求出精确解，条件改变后，方法就不一定适用。据此，为提高X射线平行入射的二维CT系统几何参数的精确性，并适用于不同模板，从多目标规划角度，提出一种多目标优化模型。当条件增加或减少时，模型精确度高、适用性强且稳定。

2 不同模板的几何参数标定优化问题及建模

对于典型的CT系统，平行入射的X射线垂直于探测器平面，将每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射－接收系统绕某固定的旋转中心逆时针旋转180次。对每个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。

2.1 第1类模板：正方形托盘上放置均匀介质圆

在100 mm×100 mm正方形托盘上放置1个均匀固体介质（半径为25 mm圆）组成的模板，模板的几何信息如图1所示。求CT系统的几何参数，即X射线入射方向（IDX）、CT系统旋转中心（COR）和探测器单元间距（DUS）。

图1 第1类模板示意图

1）确定任意1条X射线方程

将模板的正方形托盘中心与圆中心重合，建立平面直角坐标系。为进一步确定COR与512条X射线之间的几何关系，任取1个点作为旋转中心点o′（x0，y0），任取X射线的1个初始方向为θ。由于探测器单元是等距排列的，当平行光入射时，旋转中心点到每条射线的距离都可以表示，如图2。假设X射线从不同方向入射的DUS不同，用Dj来表示，则任意1条X射线到旋转中心o′的距离…，180。

图2 X射线任意1个入射角度示意图

根据图3，作经过点o′（x0，y0）的射线l0，设射线l0与y轴交点的纵坐标为a0，则射线l0的斜截式为y0＝tanθjx0＋a0，即a0＝y0－tanθjx0。任意1条射线在线段l0的下方或上方，可以得到图3所示的几何关系，则构成一个直角三角形△A0A′iAi，可以求得ai。设点Ai到A0的距离为a′i，则a′i＝因此任意1条X射线与y轴的相交点Ai的纵坐标为ai＝a0－a′i。综上，可以得到任意1条X射线的直线方程是y＝tanθjx＋ai。即：

图3 任意1条射线与l0两种位置关系示意图

2）建立超定方程组

根据定义，圆的标准方程为x2＋y2＝252，单位为mm。联立射线与圆的方程组：

3）建立第1类模板的CT系统参数标定优化模型

需要求解该模板下CT系统的几何参数，即IDX、COR和DUS。于是确定决策变量为IDX范围为θj、COR的坐标位置为o′（x0，y0）和DUS为Dj。

512条X射线穿过模板后，会在探测器接收器上显示接收信号，因此DUS需要满足一定的范围，从而使得X射线从不同方向入射时DUS最小，因此确定目标函数为180组方向下DUS的平均方差和最小，即。为减少旋转中心的偏差，当旋转中心越接近坐标原点时，512根X射线在180组方向下更容易检测。于是选择目标函数为旋转中心点o′到坐标原点o的距离最小，即

有3类条件对其约束。

约束1：X射线与圆不相交。当X射线与物体相交时，X射线强度会衰减，因此探测器上接收的信号弱。满足Δ≤0条件，即：

约束2：当不同入射方向的X射线发生强度衰减时，确定探测器单元间距范围。已知射入的X射线有512组，X射线入射角度有180组，因此不同的入射角度会影响DUS。为精确确定DUS范围，根据几何关系，用图4中3种角度入射计算得DUS范围为：

图4 X射线3种入射角度示意图

约束3：自身范围限制。在平面上，平行束IDX范围为一周，即0≤θj≤2π。由于圆在正方形托盘中，则COR一定在模板上，即

综上，得到第1类模板的CT系统多目标参数优化模型，i＝1，2，…，512，j＝1，2，…，180。

2.2 第2类模板：正方形托盘中放置非均匀介质椭圆和圆

在100 mm×100 mm正方形托盘上放置由2个非均匀固体介质（长半轴为15 mm、短半轴为40 mm的椭圆和半径为4 mm的圆）组成的模板，模板的几何信息如图5所示。其他条件相同，求CT系统的几何参数。

图5 第2类模板示意图

1）确定任意1条X射线方程

在本模板中固体介质为非均匀的，且从1个上升为2个，即增加了椭圆。为确定模板下的多目标优化模型，利用相同的建模思路。X射线不论射入椭圆还是圆都会发生衰减，根据图7的几何关系确定任意1条X射线方程过程相同，则可以确定任意1条X射线方程为y＝tanθjx＋bi，其中2，…，512，j＝1，2，…，180。

2）建立超定方程组

图6 X射线任意1个入射角度示意图

图7 射线和l0两种位置关系示意图

同样，利用弦长公式组成超定方程组，则椭圆和圆的弦长分别为：

3）建立第2类模板的CT系统参数标定优化模型

其中：k＝2，3，i＝1，2，…，512，j＝1，2，…，180。

3 适用于不同模板下平行束CT系统参数标定的多目标优化模型

为提高CT系统几何参数的准确性和精确度，并保证在样品数量改变的情况下精确度相对较高，根据式（4）和式（8）两个模型，在构成模板样品数量不同、形状不同、介质不同时，提出平行束CT系统参数标定的多目标优化模型。则利用射线与模板上样品之间的几何关系，以COR的坐标（x0，y0）、探测器单元间距（DUS）和X射线的入射方向（IDX）为决策变量，以X射线从180个方向入射时DUS方差和的平均值最小和旋转中心点o′到坐标原点的距离最小为目标函数，以X射线与介质样品不相交、3种方向入射时（即水平、竖直、东北45°方向）确定IDX范围和自身限制为约束条件，建立适用于不同模板的平行束CT系统参数标定的多目标优化模型。

其中：用Dj表示DUS，n表示整个发射－接收系统绕旋转中心的次数，j＝1，2，…，n。探测器的个数为m，i＝1，2，…，m。利用判别式Δrij表示X射线是否射入样品中，X射线被截断的弦长用lrij表示，r为X射线射入不同样品的数量，r＝1，2，3，…。探测器间距与正方形托盘和样品大小有关，因此最大间距和最小间距分别用c和d表示。COR在正方形托盘内，因此横纵坐标确定。用θj表示IDX，考虑发射－接收系统按逆时针旋转是在1周内，即θj∈［0，2π］。

4 定量分析

参考文献［21］的部分数据，拟合并多次迭代求解第2类模板的几何参数，得到结果：CT系统的旋转中心（COR）为（－9.213 4，6.035 1），探测器单元间距（DUS）为0.276 2 mm，X射线入射初始方向（IDX）为120.210 0°，每次逆时针旋转约0.997 1°，并得到180组方向。

为验证模型的精度，以第2类模板为例，将几何参数结果和9个样本［12－20］的初始IDX、COR和DUS与样本1的数据［22］进行误差分析。分别得到10个样本初始IDX、COR和DUS的误差数据，具体见表1～3。

表1 IDX误差分析

表2 COR误差分析

表3 DUS误差分析

由表1～3可知：第2类模板的初始IDX为0.563 7°＜0.617 7°；COR的误差距离0.248 5 mm＜0.640 1 mm；DUS的误差为0.000 6 mm＜0.000 8mm。3部分都低于平均值，且COR的误差距离仅次于样本3、4，而本文中DUS的误差仅为0.000 6 mm，误差非常小。从整体上看，利用提出的优化模型求解出的CT系统的参数比现有大部分结果更精确，因此所提出的优化模型可行度高，精确度也较高。

5 模型应用

为考察模型的实际使用情况，检验模型的稳定性，选择第2类模板进行应用验证。根据计算出的几何参数和文献数据［21］重建图像，并与精度较高的样本3的图像进行对比，检查成像的质量，进而说明模型的稳定性。

由于模板上的样品对X射线有吸收和散射作用，导致X射线发生衰减。将模板上介质样品的各单位近似均匀处理，得到X射线衰减系数公式根据文献［23］，基于Radon变换和反变换原理，利用Matlab内置图像重构函数iradon，采用滤波反投影算法重建图像。此时，图像需要校正，采用简单的平移、旋转，平移前后坐标用（x，y）与（x′，y′）表示，在水平和竖直方向的偏移量分别用d x，d y表示，旋转角度用α表示。

由图8、9可知，校正后的图像不仅消除了伪影，而且图像的边缘更光滑，轮廓更清晰。将本文的校正图像与样本3的校正图像对比，虽然图像重建的算法不同、模型不同，但两者差别甚微。由此可见，所提出的多目标优化模型重建的图像能达到其他方法重建的效果。换言之，说明本文中模型在保证成像质量的同时精度高，表明模型的稳定性强。

图8 第1个图校正前、本文方法校正后和 样本3校正后的图像

图9 第2个图校正前、本文方法校正后和 样本3校正后的图像

6 结论

为提高平行束CT系统几何参数的精确性，从多目标规划角度，提出一种平行束CT系统参数标定的多目标优化模型。该模型适用于不同数量、不同形状、不同介质构成的模板。利用超定方程组方程多、求解值少且便于迭代优化的特点和几何关系确定相应模板的多目标优化模型。通过分析发现，在不同数量、不同形状、不同介质构成模板的条件下，两类模板的模型结构类似，因此总结提出了一种平行束CT系统参数标定的多目标优化模型。通过对第2类模板的计算，并将9组样本分别与样本1数据进行定量分析发现，本文结果优于大部分样本，精度较高。另外，将第2类模板的几何参数结合相关数据重构图像，发现相较于样本3所重构的图像差别不大。因此，本文中所提出的模型不但适用于该类问题，还能保证精度较高、稳定性较强。

由于未对图像重建算法做深入研究，在多目标规划领域，第2类模板的计算结果是全局最优解，而实际问题中可能不存在这样的解。因此，需要进一步从多目标规划问题出发求有效解，确保最优解在实际中一定存在。