基于等径外切圆的工业CT系统空间分辨率测试

齐子诚，倪培君，姜 伟，郭智敏，高红俐

（1.浙江工业大学机械工程学院，浙江杭州310014；2.中国兵器科学研究院宁波分院，浙江宁波315103）

1 引 言

工业计算机断层成像技术（Industrial Com⁃puterized Tomography，ICT）具有限制条件少、成像直观、分辨率高的优点，广泛应用于航空航天、兵器及汽车等行业的产品内部质量检测［1-3］。目前，国内工业CT系统已成系列，包括微焦点（纳米）工业CT、15MeV大型工业CT系统等，系统的空间分辨率覆盖40 lp/mm（微焦点CT）到0.5 lp/mm（大型加速器CT），密度分辨率最小可达0.1%。空间分辨率作为表征工业CT系统性能的一个重要参数，与缺陷检测能力密切相关。因此，实现工业CT系统空间分辨率的有效测试对CT系统设计、研制、验收和应用等方面有着重要意义。

目前，工业CT系统空间分辨率测试方法按原理［4］可分为周期性结构测试法和调制传递函数（Modulation Transfer Function，MTF）测试法。周期性结构测试法［5-6］是在CT图像中直接观察线对、圆（方）孔阵列等周期性结构，利用可识别的最小结构尺寸来确定系统的空间分辨率，具有结果直观、检测效率高等优点，但是测试能力和精度完全依赖于模体中周期性结构的空间频率范围、间隔及加工精度。随着系统空间分辨率的提高，测试模体的加工要求、难度显著提升，影响了该方法的实用性。MTF法具有模体结构简单、加工要求低、检测分辨率高等优点。根据模体形状不同，该方法可分为圆盘法［7-9］和刃边法［10-15］，这两种方法均利用模体在CT图像上的边缘灰度信息，利用计算机提取边缘图像灰度分布进行微分运算形成边缘扩散函数（Edge Spread Function，ESF），再经过傅里叶变换计算获得MTF曲线，用于表征工业CT系统的空间分辨率。在实际应用中边缘图像灰度提取、平均降噪、插值拟合等过程较为复杂，影响了方法的实用性。基于上述研究现状，本文使用一对直径相同的圆柱体形成连续变间距测试模型，建立了一种可等价于周期性结构（线对卡模体）的工业CT系统空间分辨率测试方法。

2 空间分辨率测试模型

2.1 周期性结构测试模型

工业CT扫描过程如图1所示。将周期性结构（线对卡模体）固定于转台上，射线源产生X射线穿透模体并被探测器接收，模体随转台旋转一周（360°），与此同时探测器实时采集射线穿透模体后的衰减信号，并传输至上位机实现CT图像重建。

图1 条形模体在工业CT扫描示意图Fig.1 Scanning of strip phantom（line-to-card）in indus⁃trial CT

图2 CT成像过程Fig.2 Process of CT imaging

工业CT扫描形成的检测图像是空间域里物体经过系统调制后得到的结果。因此，可以利用调制度反映物体细节在背景中的对比度，进而表征CT系统的空间分辨率。线对卡模体中不同空间频率结构的CT成像前后对比如图2所示。可见经过CT成像后结构的空间频率不变，对比度降低了，形成了衰减的近似正弦波。在各个频段的周期性结构中，分别将图像最亮处（灰度最大值）记为Imax，图像最暗处（灰度最小值）记为Imin，则调制度MT可表示为［16-17］：

为探求线对卡模体的CT成像过程，先建立周期性条形结构的数学模型。在理想情况下，条形结构（线对卡）可表示为矩形脉冲函数。为研究方便，分别对波峰、波谷进行讨论，如图3所示。条形结构中的波峰p(x)、波谷q(x)分别表示为［18］：

图3 条形模体表现形式Fig.3 CT manifestation of strip phantom

其中：n∈Z，2ai为波峰（波谷）矩形的脉冲宽度，ε(x)为 单 位 阶 跃 信 号 ，定 义 为由于工业CT成像系统是一种线性非时变系统［12］，因此其成像导致的对比度下降可近似为高斯退化函数模型g(x)。成像过程可表示为输入函数与高斯退化函数的卷积，则条形结构中波峰的工业CT成像过程表示为：

同理，波谷成像过程表示为：

在不考虑噪声的情况下，卷积退化后形成的波峰（波谷）极值必然在x=0处，则波峰Imax（波谷Imin）表示为：

分析各个条形结构之间的影响关系，简化波峰Imax（波谷Imin）的表达式。假设工业CT成像过程的高斯退化函数模型为：

根据标准正态分布规律可知，在3倍标准偏差的区间内，g(x)的积分函数面积约占整体面积的99.73%。当波峰（波谷）宽度2ai等于标准差σ时，在3倍标准偏差的区间内包含了3组线对，此时的MT值约为1.1%，可见对波峰（波谷）幅值影响较为显著的是相邻的两组线对。因此，当1.1%≤MT≤100%时，波峰Imax可近似为：

将式（6）、式（8）代入式（1），即可获得线对卡CT成像后调制度MT的表达式：

2.2 等径外切圆结构模型分析

双球法［19］是Carmignato等人提出的一种三维CT系统空间分辨率评估方法，主要是利用一个简单结构的参考试块（两个具有相同公称直径的点接触校准球体），对该参考试块进行CT重建，观察接触点的失真情况。在指定误差范围内，利用可测量的最小球体直径来确定CT系统的空间分辨率。该方法的优点在于参考试块结构简单，但是对于具有微米分辨率的工业CT系统则需要直径在微米级的球体，球体试块制造、校准及使用难度显著提升。Zanini等［20］提出了一种基于双球的改进空间分辨率测试方法，采用较大直径的球体，利用局部自适应算法提取球体边缘，计算接触点扭曲区域的高度来表征CT系统的空间分辨率，但是在水平方向上由于射束硬化、噪声等原因，接触点附近灰度分布特别复杂，影响扭曲区域高度的测量精度和稳定性，测量结果发生振荡。本文构造了一种垂吊式圆柱体线接触模体，用于改善CT检测效果，测量结果等价于线对卡法，更适用于线阵工业CT系统的空间分辨率测试。

如图4所示，两个具有相同公称直径(D)的校准柱体，使它们在一条线上相互接触。由于CT重建过程造成接触区域图像畸变，该畸变区域的面积随着CT系统空间分辨率的提高而减小。因此，两个相互接触的校准柱体可作为一种狭缝宽度连续变化的矩形波测试卡［21］。

图4 双圆模体与空间分辨率的关系Fig.4 Relationship between double circular phantom and spatial resolution

图5 双圆模体平面Fig.5 Double circular phantom plane

如图5所示，两个相邻的圆柱形结构平行于圆心线段的间距可表示为：

与圆心线段距离为h的灰度分布函数可表示为：

其中：ah为与圆心线段距离为h的双圆间距半宽。需要注意，f(x，ah)表示归一化的灰度值，且0≤f(x，ah)≤1。

将x=0代入式（11）可得不同间距下中轴线上的CT图像灰度分布，有：

由于正态分布函数g(x)的积分无解析解，引入g(x)积分近似计算公式［22］：

其中W≥1。

2.3 等径外切圆与周期性结构等效分析

将线对卡CT成像后的调制度MT引入g(x)积分近似计算公式G(x)，即将式（13）代入式（9），则MT可表示为：

如图6所示，当线对卡的线对波峰（波谷）宽度2ai与双圆模体的间距2ah相等时，将式（14）代入式（15）可得：

则线对波峰（波谷）宽度为2ah（2ai）对应的空间频率f（单位：线对数/毫米）可表示为：

图6 双圆模体与线对卡宽度相等Fig.6 Width of double-circle phantom equals to the one of line-to-card

其中k为修正系数，约为1.1～1.2。

式（17）转换可得到：

将式（18）与式（16）联立即可获得利用双圆模体测得的等效线对卡MTF曲线。

在工程应用中，通常采用10%调制度下的线对宽度来比较工业CT系统的空间分辨率，因此将MT=10%代入式（16）并进行简化，得到：

对式（19）进行多项式求根，取大于1的实数解，可得W≈2.810 49。

将W≈2.810 49代入式（14），求得：

因此，10%的调制度下CT系统的空间分辨率为：

2.4 等径外切圆MTF自动测量算法流程

为了快速获取双圆模体在中轴线上的灰度分布，设计了一种图像处理及测量算法实现中轴线灰度分布的提取及MTF曲线计算。基于等径外切圆的工业CT系统MTF自动测试方法主要包含几个重要的步骤：CT图像二值化；模体质心计算；双圆模体圆心计算；圆心中轴线上图像灰度提取；计算CT系统的MTF曲线。数据处理流程图如图7所示。

图7 CT图像数据处理流程Fig.7 Data processing of CT images

假设双圆模体所在CT图像f(x，y)的尺寸为m×n，其中x为图像数组横坐标，y为图像数组纵坐标。如图8（a）和8（b）所示，利用自动阈值分割方法（clustering）转换为二值化图像b(x，y)，则b(x，y)的质心C位置坐标(xc，yc)可表示为：

其中num为满足条件的i或j的数量。

为了能够准确识别双圆心位置，需要对两圆进行有效分离。如图8（c）所示，对b(x，y)进行二值形态学处理，获得双圆边缘图像e(x，y)，表示为：

其中：Θ表示腐蚀，即b(x，y)被结构B所腐蚀采用K-means聚类方法计算双圆圆心坐标，其主要步骤如下：

（1）首先，搜索e(x，y)中值为1且与质心C距离最远的点，设为N。初始化类A={C}，并将质心C的位置作为类A的初始中心位置；初始化类B={N}，并将点B的位置作为类B的初始中心位置；

（2）在e(x，y)中寻找值为1的点，即双圆的边缘位置，计算该点与类A和类B的中心点的距离。当该点与类A的中心距离小于类B的中心距离时，将该点计入类A，重新计算类A中所有点的中心位置；反之亦然。

图8 双圆模体圆心自动提取结果Fig.8 Automatic extraction of center of double circle phantom

（3）迭代计算步骤（2）直至e(x，y)中所有值为1的点处理完毕，即可获得圆心坐标分别为L=(xL，yL)和R=(xR，yR)。

如图8（d）所示，利用圆心坐标(xL，yL)和(xR，yR)计算中轴线l上的点坐标(xl，yl)，可表示为：

计算中轴线l上每个点到圆心连线段LR的距离hl，则有：

其 中：A=xR−xL，B=(yL−yR)/(xL−xR)，

将hl代入式（17）即可获得相应的线对空间频率。

在双圆模体CT图像f(x，y)中，采用数值微分法提取点坐标(xl，yl)对应位置的像素灰度值f(xl，yl)。将灰度值f(xl，yl)归一化并进行分段线性拟合处理，代入式（14）、式（16）即可获得调制度MT。图9分别表示中轴线的灰度分布和MTF曲线的计算结果。

图9 中轴线灰度分布及MTF曲线Fig.9 Gray scale distribution of central axis and MTF curve

3 实 验

3.1 试块的制作

选用牌号为GCr15的轴承钢作为双圆模体材料，该材质具有淬透性好、硬度高、耐磨性强和尺寸稳定等优点。首先在原料条上进行断切，得到两个外形尺寸均为Ф9 mm×16 mm的坯料；进行研磨加工及表面抛光，得到两个Ф8 mm×15 mm的圆柱模体，接着用SGM15/14型高温炉进行淬火及低温回火，使其硬度达到45～50；对其表面进行防锈处理，最后采用三坐标仪对模体尺寸进行计量，验证其尺寸公差是否满足工业CT空间分辨率测试的要求。采用垂吊方式保证模体圆柱侧面紧密贴合，如图10所示。

图10 双圆模体示意图Fig.10 Double circular phantom

实验前，将双圆模体置于20℃条件下恒温2 h；实验时，将模体置于CT设备转台上，选取模体的中部区域进行工业CT扫描成像。实验所用的设备为高能加速器工业CT成像系统，由6 MeV加速器射线源、608通道线阵探测器、双立柱运动平台和Ф1 m转台组成。设备探测器单元水平方向尺寸为0.3 mm（固定），相邻通道水平间隔为1.3 mm，探测器单元垂直方向尺寸（切片厚度）可利用准直器进行调节，范围为0.25～5 mm，最大剂量率约为800 cGy/min×1 m，焦点尺寸为2 mm，采用三代CT扫描方式。

3.2 有效性分析

为验证本文方法的有效性，将它与线对卡法［16-17］进行比较，研究本文方法的测试结果受空间分辨率、噪声水平、稳定性和双圆模体直径变化的影响，以及与线对卡法测试结果的差异。如图11所示，线对卡宽度分别为0.1/0.15/0.2/0.25/0.3/0.4/0.5 mm。

图11 线对卡模体及CT图像Fig.11 Wire pair card phantom and CT image

微动次数是指探测器单元通过小范围机械运动形成过采集插值实现超分辨率重建，是影响高能工业CT系统空间分辨率的主要因素。微动装置采用丝杆螺母副传动、交叉滚子滑动模块导向实现高精度小范围平移，将探测器相邻单元间距按微动次数进行等分，利用光栅尺测量探测器单元的实际位置用于软件算法的插值计算。因此，通过改变工业CT系统检测时的微动次数可形成不同空间分辨率的CT图像。图12（a）～12（c）分别给出了2次、5次和10次微动获得的线对卡和双圆模体CT图像，实验采用的切片厚度为1 mm、图像分辨率为4 096×4 096（像元×像元）、成像范围为200 mm×200 mm。采用两种方法测得10%调制度下的空间分辨率如图13所示，可见两种方法获得的10%MTF的空间分辨率基本吻合，最大误差在10%以内。

研究不同CT图像分辨率对MTF测试结果的影响，实验采用2次微动、切片厚度为1 mm、成像范围为200 mm×200 mm，如图14所示。增加图像分辨率有助于提高检测系统的空间分辨率，但是当空间分辨率接近系统极限时，增加图像分辨率对系统空间分辨率的提高效果不明显且会大量消耗硬件内存。从测量结果可见，两种方法均体现出图像分辨率对MTF测试结果的影响，在高图像分辨率下（8 192×8 192），本文方法测量值（10%调制度）略高于线对卡测量值，原因在于线对卡法不同空间分辨率线对的间隔较大，10%调制度下的实际空间分辨率需要通过插值拟合计算获得，而双圆法的等效线宽间隔在10%调制度下远小于线对卡法。因此，双圆法的测量精度高于线对卡法，同时在日常检测时应选择合适的图像分辨率以获得最佳的空间分辨率及最优的系统开销。

图12 不同微动条件下CT图像Fig.12 CT images with different fretting frequencies

图13 不同微动频率下测试结果比较Fig.13 Comparison of test results with different fretting frequencies

图14 不同图像分辨率下测试结果比较Fig.14 Comparison of test results with different image resolutions

图15 不同切片厚度下测试结果比较Fig.15 Comparison of test results with different slice thicknesses

研究不同噪声水平对MTF测试结果的影响，实验采用2次微动、图像分辨率为4 096×4 096（像元×像元）、成像范围为200 mm×200 mm，如图15所示。切片厚度是影响CT图像噪声水平的主要因素，前期实验表明切片厚度越大，图像噪声水平越低。可见，噪声对本文方法具有一定的影响，噪声越大与线对卡法的结果差异越大，原因在于本文方法通过对中轴线上灰度分布进行分段线性拟合来降低图像噪声的干扰，降噪效果有限。通过分析，两者误差控制在10%以内，基本满足检测要求。

图16 不同直径模体的测试结果比较Fig.16 Comparison of test results of phantom with differ⁃ent diameters

研究不同直径模体对MTF测试结果的影响，实验采用5次微动、图像分辨率为4 096×4 096（像元×像元）、成像范围为200 mm×200 mm、切片厚度为1 mm，双圆模体直径分别为3，6，8，20，26，32 mm，如图16所示。可见，圆柱体直径的变化不会显著影响空间分辨率测量结果，噪声、伪影是引起不同直径模体测量结果波动的主要原因。

4 结 论

本文针对传统线对卡法测量工业CT系统空间分辨率的局限性，推导了线对卡法测量模型及等径外切圆变间距模型，建立了两者的对应等价关系，构建了基于等径外切圆的MTF测量方法，分析了测量结果受空间分辨率、图像分辨率、噪声水平及模体直径的影响规律，并采用GCr15轴承钢模体实验验证了本方法的有效性。实验结果表明，本方法的测量结果接近线对卡法，误差控制在10%以内，具有模体加工要求低、测量过程简单的优点，可推广用于CT系统的设计、研制、验收和应用等方面。下一步将针对等径外切圆法的降噪方法及伪影对测量结果的影响开展研究，使本方法的适用性更广。