电气设备高速旋转物体的CT图像重建方法

王朝华 陈德峰 刘荣海

（1．国网河南省电力公司电力科学研究院，河南 郑州 450052；2．检测成像北京高校工程研究中心，北京 100048；3．云南电网有限责任公司电力科学研究院，云南 昆明 650000）

0 引言

CT系统工作时，X射线被前准直器准直成一个狭窄的射束，通过被检测的工件后被探测器组接收，得到多个角度下射线的强度信号，这些电信号经过光电倍增管或光二极管放大后，再经过模数转换器转换成数字信号，通过计算机运用一定的数学方法处理后，在显示器上得到一个重建的、完整的二维层析图像。由于在同一断面上各点材料的密度不同，因此各点的线性衰减系数也不相同，这样，吸收X射线的量也不同，这些差异可在图像上显示出来，为研究、判断分析提供了可靠的依据。

国内外现有工业X射线CT还无法实现对高速旋转物体（例如高速运转中的航空发动机、高速离心机以及燃气轮机等）的CT成像[1-2]。有关高速旋转物体CT成像的研究尚处在探索阶段。2003年提出了基于随机正交投影识别和重排的高速旋转物体CT成像方法，并对该方法进行了仿真数据验证，但是该方法要求所采集的数据仍是准静止且不发生数据混叠的。2014年相关学者分析了研制航空发动机在线成像CT系统面临的难点问题，给出了系统设计、数据采集方案以及图像重建方法，提出了降低由数据混叠引起的图像模糊方法，并构建高能CT系统进行模拟成像测试。高速旋转物体CT成像的主要难点是受X射线源流强和探测器数据采集速度的限制，无法在物体高速旋转的条件下采集准静止的CT数据，因此经典的CT成像模型不再适用，如果仍采用经典CT成像模型将导致重建图像在沿旋转方向出现模糊。

1 数学模型

CT常用算法大体上可分为3类，第一类为系列重建法或称系列（级数）展开法，例如联合迭代重建法；第二类为变换重建法，例如傅里叶变换重建法等；第三类为其他重建法。目前，变换重建法中的卷积反投影法应用较多，常用的有平行束卷积反投影法和扇形束卷积反投影法。

该文在若干假设下，针对电网设备高速旋转物体CT成像给出了数据采集时间的选取准则，建立了射线混叠的CT成像数学模型，并提出了相应的图像迭代重建算法。当数据存在噪声时，射线混叠的CT成像数学模型具有高病态性。经模拟CT数据和实际CT数据验证，该文所提出的算法在数值上具有较好的鲁棒性，可以正确重建高速旋转物体的CT图像，且不会在旋转方向产生图像模糊。

每个CT扫描数据Ij可以看作是入射流强I由被测物体x经过一定调制过程Fj得到的数据，即Ij=Fj（I，x）；而CT图像重建过程则可以看作是由一系列关系Ij=Fj（I，x）解调出x（j∈J，J是扫描数据指标集）。在此意义上，该文的思想与利用编码孔径和利用准直器提高分辨率的思想有相似之处，而不同之处在于相应的调制过程和解调方法，类似的思想还有基于虚拟焦点的超分辨重建方法。该文对高速旋转物体CT成像做了一些较为理想的假设，而实际系统可能会有旋转不匀速、旋转轴晃动的情况，但相应的CT成像问题本质上仍是调制解调问题，只是在调制关系Ij=Fj（I，x，Θ）中增加了参数族Θ，因此需要同时研究Θ的测量或估计方法和图像重建（解调）算法。

假设：被检测物体是匀速转动的；旋转轴不发生晃动；探测器的采样时间是连续可选的；X射线是单能的；忽略多次散射光子以及电子对淹没产生的光子。在该假设下，高速旋转物体断层CT成像如公式（1）所示。

2 高速转动数据快速处理方法

高速转动图像重建具有数据量大的特点，因此在工程应用中考虑重建算法的计算效率也是CT技术的一个重要方面。目前对反投影部分的CUDA加速策略主要如下：1）在线程分配时考虑反投参数的Z轴相关性。每个线程各自完成对应Z轴上一系列点的重建任务，这种线程分配方式可以减少GPU的重复计算量。2）纹理存储器的使用。纹理存储器能够被缓存，因此使用纹理存储器可以提高访问速度，并且纹理存储器带有优化过的寻址系统，支持方便、高效的二维插值操作，可用来完成反投影中的双线性插值运算。3）对全局存储器釆用合并访问优化。在满足合并访问条件的情况下，只需要一次传输便可以处理来自一个half-warp的16个线程对全局内存的访问请求。4）对三角函数计算的优化。在GPU中完成一次三角函数运算需要32个时钟周期，而一次乘加只需要4个时钟周期，因此在计算三角函数时釆用递归方法，将不同角度下的三角函数运算转换为一次乘加运算，便可在一定程度上缩短GPU的运算时间。由于GPU架构本身在高性能计算方面的优势，在应用上述优化策略后， 便可获得CPU上百倍的加速比。在上述并行策略的基础上，做了更细致的优化，使优化后的加速效果获得大幅提升。

与上述优化策略相比，该文所提出的方法与现有方法的不同主要体现在以下3个方面：1）线程分配方式。2）通过使用常数存储器存储计算反投影参数时的中间量，减少GPU的重复计算量。3）通过优化一个kernel中反投影的张数来减少访问全局存储器的次数。

3 数值试验

首先，采用仿真数据验证该文提出的数据采集和重建算法的有效性。其次，利用实际工业CT系统采集的数据验证该文算法对多能CT数据的有效性。

3.1 仿真数据验证

根据以下条件模拟高速旋转物体的扇形束CT数据：射线源为能量8 MeV的单能射线源；探测器为线阵列探测器，由512个探测器单元组成，各单元的尺寸为0.60 mm。被测物体模型为如图1（a）所示的铁材质工件，环形的外径为208.30 mm，内径为158.70 mm，其最长弦为134.91 mm，铁对8 MeV的光子的线性衰减系数为0.023 63 mm-1。射线源到旋转中心的距离为1 000.00 mm，射线源到探测器中心的距离为1 200.00 mm。射线混叠数据由旋转17圈（6 120 °）采集的720 个投影数据得到，即每旋转8.5 °采集1幅投影数据。扫描数据中加入泊松噪声的方法如下：由已知光子流强I0和被测物体模型，根据公式（6）数值计算；再以为泊松分布的均值产生一个随机值，并用该随机值替换Ik，j。 由此可知，数据的噪声水平与I0和物体吸收率是相关的。对给定物体来说，I0越小，数据噪声水平越高。

图1 理想仿真数据以及重建图像

分别考察重建算法对无噪声数据和有噪声数据的重建效果。图1分别为用3种方法对无噪声的混叠扫描数据重建的CT图像，重建图像矩阵为512×512。使用FBP算法直接重建的图像如图1（a）所示，第 4.1 节所述的直接解调算法重建的图像如图1（b）所示，第3.2节提出的迭代算法重建的图像如图1（c）所示。图1（a）在旋转方向存在明显的图像模糊，远离旋转中心的像素模糊程度相应加重。在无噪声的情况下，由直接解调方法和该文提出的迭代算法均能重建满意的图像。

3.2 实际数据验证

该文在射线混叠CT数学模型中假设射线是单能的，然而工业CT的射线源所发出的射线是多能的。为验证该文方法对多能射线的有效性，用工业CT设备扫描的数据进行试验验证。采用延长探测器采样时间的方法来模拟扫描数据混叠。试验所用的工业CT系统探测器单元个数为2 604个，探测器单元宽度为0.25 mm，探测器的采样时间为60 ms。射线源焦点到转台中心的距离为600.00 mm，射线源焦点到探测器的距离为1 139 mm。将叶轮放置于转台旋转中心处进行CT扫描。转台共旋转7圈（即2 520°），探测器的采样角度数为1 800 个，即每旋转1.4 °采集1幅投影数据。叶轮实物照片和扫描数据如图2所示。在叶轮中加入4根铅笔芯作为标记物，以考察重建图像沿旋转方向的模糊程度。

图2 叶轮实物照片和混叠扫描数据

图3（a）、 图3（c）和 图3（e）分别为直接FBP算法、直接解调算法和该文迭代重建算法重建的图像，图3（b）、图3（d）和图3（f）分别为它们的局部放大图像。由直接FBP算法重建的图像发生沿旋转方向的模糊，直接解调算法已难以重建出叶轮的信息，而该文迭代重建算法仍可以重建出质量较好的图像。

图3 叶轮实际扫描数据重建图像

4 结语

该文在若干假设下，针对高速旋转物体的CT成像问题建立了射线混叠的CT数学模型，提出了CT采样时间的选取准则及其相应的迭代重建算法，并通过仿真数据和实际数据验证了该方法的可行性。尽管该方法做了射线单能的假设，但是实际工业CT数据验证表明，将该方法用于多能数据也具有很好的近似性。该文是对高速旋转物体CT的初步研究，在物体转速有扰动以及旋转轴有晃动等情况下，还需要进一步深入研究如何重建高速旋转物体的CT图像。