基于优化的CT系统标定与滤波反投影图像重建

李江珊，邹义华，姚 竹，王志勇

(1.电子科技大学 信息与通信工程学院，四川 成都 611731；2.电子科技大学 数学科学学院，四川 成都 611731)

本文是对2017年全国大学生数学建模竞赛题目的求解，题目和数据由组委会提供[1]。现有某类型的CT系统，平行的X射线入射到样本上，经过衰减后在另一端被探测器接收，探测器有512个等距的探测单元，可测量接收到的X射线的能量。将整个发射探测系统绕一个固定的旋转中心旋转180个方向，经增益处理后可得到180组数据[2]。为减小CT系统的误差，论文建立了基于几何分析和优化的标定模型。在此基础上，通过滤波反投影法对图像进行了重建。最后，利用计算机模拟分析了标定模板的精度和稳定性，并进行了新模板设计。

1 模型的建立与求解

1.1 基于直线方程和曲线方程的CT系统标定模型

CT系统的标定随着CT系统的发展不断发展，其大致可分为两类，一类是自标定方法，另一类则是基于标定模板的标定方式[3]。自标定方式是指对CT系统参数进行标定时仅根据被检测物体本身的投影。而不依靠其他数据。该方法增加了计算量，无法获得实际分辨率并且具有局限性，所以使用较少。而基于标定模板的标定方式是通过在特定的运动轨迹下采集投影，然后与CT系统参数建立联系，从而得到CT系统的几何参数。本文就是建立在第一代CT系统的基础上，利用基于标定模板的标定方式对系统进行标定。标定的具体方法有基于统计特征的提取法，最小二乘法，自适应遗传算法等[4]，其中自适应遗传算法涉及的参数较多，计算量较大[5]，并且遗传算法不能保证得到的解为最优解；最小二乘法则是需要建立多个目标规划模型[6]，较为复杂。本文将问题分割为确定增益系数和探测器单元间距、确定CT系统X射线的180个方向、确定旋转中心位置等多个模块，相较于自适应遗传算法减少了参数，降低了计算难度，并且在计算CT系统方向部分直接利用几何代数进行求解，相较于最小二乘法，也更为简单。

为了更加直观地观察题目给出接收信息，画出接收信息随组数、探测器单元的变化曲线，如图1所示。显然，中间类似花瓶形状的凸起是由模板中的椭圆导致的，弯曲曲线是由模板中的圆导致的。基于以上认识，下面给出基于直线方程和曲线方程CT系统标定模型。

1.1.1 确定增益系数和探测器单元间距

根据文献[7]得到比赛数据附件2[1]中的接收信息η与相应X射线通过介质的长度l成正比：

图1 接收信息随组数、探测器单元的变化曲线

η=Al

(1)

取出位于圆心同一侧的经过了圆介质的3条相邻X射线的接收信息，假设圆心到其中一条射线的距离为d0，d表示相邻探测器单元之间的距离，R表示圆的半径，可得：

(2)

(3)

(4)

1.1.2 确定CT系统X射线的180个方向

分析得到射线的方向和探测器的初始位置大致如图2所示，以正方形托盘的中心为原点建立如图2所示的坐标系。

推知第k个方向上编号为n的射线方程式为：

ykn=(tanθk)x+bkn

(5)

图2 CT系统示意图

为了求出各组数据对应的平行射线与x轴正方向的夹角，取比赛数据附件2中两条经过椭圆的射线对应的接收信息值ηk4、ηk5，根据式(1)和椭圆弦长公式可以得到以下两组方程：

(6)

(7)

其中，

ik、jk分别表示第k组数据中经过椭圆的两条射线的编号。

1.1.3 确定旋转中心位置

根据直线绕点旋转的性质可知，旋转中心到旋转过后直线的距离与旋转中心到未旋转直线的距离相等。

(8)

式中，(x0,y0)表示旋转中心坐标，ki、bi分别表示第i条射线方程的斜率和截距。

设定旋转中心位于正方形托盘内，旋转中心到直线的理论距离值必不大于正方形托盘的内最长线段(对角线)的长度，即得到优化模型的约束条件：

(9)

式中，(x0,y0)表示旋转中心坐标，r表示旋转中心到直线的理论距离值。求解优化模型即可得到旋转中心的位置。

1.1.4 模型结果与分析

将增益系数值和探测器单元间距代入模型求解，求解出在180个方向上对应的180组平行射线与水平方向夹角结果，将其绘制成柱状图形式如图3所示。平行射线与图2所示坐标系的x轴正方向的夹角范围为(-60.385 3°,118.695 9°)，相邻的两个方向的角度差集中在1°附近。

利用遗传算法得到旋转中心坐标(在图2所示的坐标系中)的10组数据，将这10组数据的均值作为最终的旋转中心坐标，则最后确定的旋转中心坐标为(-9.268 9,6.274 3)，位于第二象限，较为靠近原点。

图3 CT系统X射线的180个方向示意图

1.2 基于滤波反投影法的CT图形重建模型

根据文献[2]采用滤波反投影模型重建附件3[1]的CT图形。

论文的滤波反投影法分为4个部分：滤波器选择、线性内插，反投影重建和确定吸收率。

1.2.1 滤波器选择

论文分别采用R-L滤波器[8]、S-L滤波器[9]、R-L和S-L混合滤波器、R-L和NEW混合滤波器进行滤波处理，比较滤波效果。

1.2.2 线性内插

根据文献[10-11]，本文利用每个单元格集合中心的坐标来表征整个单元格，记为(i,j)，但(i,j)可能在n0射线和n0+1射线之间，因此需要进行线性内插预测经过点(i,j)的射线的接收信息，即：

p1(xij,θm)=

(10)

其中

(11)

Δ=xij-(n0+N0)d

(12)

式中，n0表示512条射线中，位于点(i,j)左侧且距离点(i,j)最近的射线编号；Δ表示(i,j)到n0射线的距离。由此即可求出在任意方向上虚拟的探测器单元接收到的经过点(i,j)的射线信息。

1.2.3 反投影重建

(13)

式中，p1(xij,θm)=p(xij,θm)f(xij)，p(xij,θm)表示接收信息，f(xij)表示滤波冲激响应函数。根据上述分析即可计算得到规模为256×256的单元的具体衰减系数。

1.2.4 确定吸收率

论文认为衰减系数与吸收率成正比[17]，通过建立优化模型来求得最优的系数值。

利用滤波反投影法对附件2[1]数据进行重建，得到重建后的吸收率衰减系数矩阵与附件1中的实际数据进行比较。采用实际吸收率与理论吸收率的之差的绝对值之和为目标函数，则优化模型的目标是：

(11)

式中，k0(i,j)表示单元格(i,j)的实际吸收率。该优化模型是无约束优化模型。

结果如表1所示，R-L和NEW混合滤波器的滤波效果最好，因此采用R-L和NEW混合滤波器来重构附件3[1]和附件5[1]的图像。

表1 不同滤波器的滤波效果

利用滤波反投影模型得到附件3[1]和附件5[1]数据重建后样本的几何形状和在正方形托盘中的位置分别如图4所示，其中白色部分表示有介质。而附件4[1]中10个位置处的吸收率具体值如表2所示。

图4 附件3和附件5样本形状与位置

表2 10个位置的吸收率

1.3 标定模型精度、稳定性分析与模板改进

1.3.1 标定模型的精度分析

模型的精度是指计算值与真实值的接近程度。论文通过模拟已知的实际情况，设定发射-接收系统的旋转中心，相邻探测器单元间距和180个方向角为θk，利用计算机模拟出接收信息。再根据模型2.1利用模拟得到的接收信息进行CT系统标定，将计算结果与实际结果进行对比，得到探测器单元之间的间距的相对误差约为3×10-6，X射线与水平方向夹角与实际值的误差平方和值为0.020 9。根据上述结果，可以看出参数标定模型的精度很高。

1.3.2 标定模型的稳定性分析

因为接收信息的获得过程是未知的，所以接收信息具有不确定性，需要检验外界扰动对标定结果的影响，即分析标定模型的稳定性。为了描述接收信息的不确定性，给每个接收信息加上一个在区间[0,0.01]之间的随机扰动。根据模型2.1利用新的接收信息来计算增益系数、间距等参数。通过未扰动数据计算的结果对比，旋转中心的横纵坐标的绝对误差分别为-0.128 8、0.057 5，相对误差为1.41%、0.92%。探测器单元之间的间距的相对误差为-0.83%，增益系数A′的相对误差为0.24%。计算得到X射线与水平方向夹角与实际值的误差平方和值为0.006 3。根据上述结果，可以看出参数标定模型的稳定性较好。

1.3.3 标定模板与标定模型的改进

模板的尺寸和位置如图5所示。

(a)长方-圆模板

(b)长方-椭圆-圆模板

利用模型2.1确定增益系数、探测器单元间距、CT系统X射线的180个方向和旋转中心位置的方法，重新建立相关方程，并对新标定模板重新分析精度和稳定性。

对于长方形-圆形模板，其计算增益系数和探测器单元间距的方法与原模板相同，则与原模板相比，增益系数和探测器单元间距的精度和稳定性不变。而其算出的X射线与水平方向夹角与实际值的误差平方和值为0.229 9，比原模板的误差平方和大，这是因为采用新模板计算得到在90°附近的夹角值与实际值相差很大，但在其他位置计算值与实际值相差十分小，有些甚至没有误差。因此，可以认为新模板的精度有所提升。给模拟的接收信息增加一个在区间[0,0.01]之间的随机扰动，利用扰动后的数据对CT系统重新进行标定，得到其与未扰动结果之间的误差平方和为0.002 8，比原模板小，即模型的稳定性更高。

对于长方形-椭圆-圆形模板，其计算增益系数和探测器单元间距的方法与原模板相同，则与原模板相比，增益系数和探测器单元间距的精度和稳定性不变。其计算出的X射线与水平方向夹角与实际值的误差平方和值为0.006 8，比原模板和长方形-圆模板的误差平方和小，说明新模板的精度比原模板和长方形-圆模板的高。加上一个在区间[0,0.01]之间的随机扰动，利用扰动后的数据对CT系统重新进行标定，得到其与未扰动结果之间的误差平方和为0.038 4，比原模板大，即模型的稳定性更低。

2 结束语

模型2.1先通过圆确定增益和相邻探测单元的间隔，减少优化模型的决策变量，简化模型，便于计算。模型2.2利用多种滤波器对附件二求解吸收系数，通过对比选择了较优的滤波器处理附件三和附件五。

模型2.2中预测未知点的接收信息时，可以考虑多项式插值、样条插值等方法，提高图形精度。