CN101158653A - 一种锥束ct系统的散射测定和校正方法 - Google Patents

Abstract

本发明公开了一种锥束CT系统的散射测定和校正方法，设定参数后采集空气投影图像和射束衰减网格投影图像，对被检测物体进行圆周扫描，采集带射束衰减网格和被检测物体的投影图像集I和被检测物体的投影图像集II，计算射束衰减网格中每个金属小球中心的投影位置，采用射束衰减网格校正方法计算与投影图像集I中的投影图像一一对应的散射场分布图像，将投影图像集I减去对应的散射场分布图像得到散射校正后的投影图像集III；由此通过滤波反投影重建算法重建出散射校正后的序列切片图像。本发明能实用于从低能到中、高能的锥束CT系统散射校正，是锥束CT系统的一种简单、有效的散射校正方法。

Description

一种锥束CT系统的散射测定和校正方法

技术领域

本发明属于无损检测领域，涉及锥束CT系统散射测定和校正方法。

背景技术

高分辨率锥束CT(Cone Beam Computed Tomography，CBCT)由X射线源产生锥形束射线，利用面阵探测器采集被检测物体的投影图像，是近年来国际上解决无损检测问题最有发展前途的一种高新技术。与传统的二维CT相比，CBCT具有很高的扫描速度，在一个扫描周期内，可完成大量断层图像的扫描，获得XYZ三向同性的高精度空间分辨率图像，并且重建出的切片图像具有切片序列连续、切片内和切片间的空间分辨率相同、精度高等特点。

目前，制约锥束CT系统广泛应用的技术瓶颈是重建切片的图像质量，散射(Scatter)是影响图像质量的一个重要因素。一般而言，被检测零件的密度越大，X射线穿透物体所需的扫描电压越高，由散射引起的伪影越严重。对于散射问题，许多学者做了大量的研究。影响散射的因素主要有X射线的能量、零件的材料及几何形状、零件的厚度等。目前较为实用的散射校正方法有以下几种：

(1)基于三能量窗(TEW)的散射校正方法：通过主能量窗获取总光子数量，并采用主能量窗两边的两个子能量窗估计散射光子数量，由总光子数量减去散射光子数量得到透射光子数量。基于三能量窗的散射校正大量工作放在了主能量窗和两个子能量窗的位置和宽度上。两个子能量窗尽可能窄和近的靠近主能量窗，以便精确估计当前主能量窗的散射部分。然而，子窗口宽度越小，探测器探测到的光子数量越小，噪声的估计准确性越差。

(2)基于多分辨分析的散射校正方法：该方法首先获取原始扫描图像f；然后对投影图像f采用小波多分辨分析，每次只对低频分量逐级进行小波分解，得到第n级逼近A_nf；只对n级进行重构得到的图像f′，作为散射的估计；然后对投影图像做散射校正f＝f-f′；最后对校正后的投影图像进行滤波反投影重建。这种方法的典型处理见于公开文献：李永利，刘贵忠，潘德恒，X射线成像中图像的综合校正方法研究，兵工学报，2002，23(2)：196-200。但是基于多分辨分析的散射校正方法对小波类型和小波分解级数要求较高，选择不当难以实现，而且方法本身具有较大的误差。

(3)基于Beam Stop的散射校正方法：这种方法见于公开文献Ruola Ning，Xiangyang Tang，D.L.Conover，X-Ray scatter suppression algorithm for cone beam volume CT.Proc.SPIE vol.4682，1605-7422，2002。该方法采用反散射网格(Beam Stop)铅球吸收X射线，阻止X射线透过物体被探测器接受，这样在探测器铅球中心投影位置上得到该位置的散射值，利用这些散射值拟合散射场，进行散射校正。反散射网格置于X射线源和检测物体之间，由许多小铅球组成的阵列，该装置要求合理选择铅球的尺寸，在能保证X射线能被铅球完全吸收的同时使得其投影图像的阴影越小越好。实验中需要采集一组带反散射网格的投影图像I和另一组不带反散射网格的投影图像II。图像I中探测器探测到的铅球阴影位置上的光子数假定为是X射线的散射光子，然后采用三次样条曲面插值估计散射场分布，将图像II减去对应的散射图像得到散射校正后的投影图像，最后进行滤波反投影重建得到切片图像。该方法要求反散射网格完全挡住X射线光子透过物体，对医学上较低扫描电压时有较好的效果，但在工业CT检测中，通常需较高的扫描电压，该方法难以保证实验中铅球尺寸的需求。

发明内容

为了克服现有技术不能有效应用于工业无损检测中的不足，本发明提供了一套锥束CT系统的散射测定方案及其校正算法，以解决中高能锥束CT系统在工业应用中的散射校正问题，从而获得高质量高精度的重建切片图像。

本发明解决其技术问题所采用的技术方案包括以下步骤：

1.根据被检测物体确定对其进行锥束CT扫描的投影放大比、平板探测器的数据采集方式、X射线源的电压与电流等扫描参数，这些参数在以下涉及投影图像数据采集的步骤中保持不变；

2.不放置任何物体，通过平板探测器采集若干幅空气投影图像，取其平均值作为待用的空气投影图像；

3.将射束衰减网格(Beam Attenuation Grid，BAG)置于X射线源和旋转工作台之间，通过平板探测器采集若干幅射束衰减网格的投影图像，取其平均值作为待用的射束衰减网格投影图像；

4.将被检测物体置于旋转工作台上，进行圆周扫描，扫描过程中旋转工作台可采用步进旋转或连续旋转，通过平板探测器采集带射束衰减网格和被检测物体的投影图像集I；

5.移除射束衰减网格，仅对被检测物体进行与步骤4完全相同的圆周扫描，其中需要特别注意的是被检测物体的起始扫描投影位置与步骤4的起始扫描投影位置相同，通过平板探测器采集被检测物体的投影图像集II，投影图像集II与投影图像集I的投影幅数相同，且零件各个投影位置对应相同；

6.通过步骤3得到的射束衰减网格投影图像计算射束衰减网格中每个金属小球中心的投影位置；

7.通过金属小球中心的投影位置和步骤2、4、5得到的投影图像，采用射束衰减网格校正方法，计算与投影图像集I中的投影图像一一对应的散射场分布图像；

8.将投影图像集I中的各幅投影图像象素灰度减去与其一一对应的散射场分布图像象素灰度，得到散射校正后的投影图像集III；

9.由投影图像集III通过滤波反投影重建算法重建出散射校正后的序列切片图像。

上述各步骤采集的投影图像在获取过程中已经过必要的暗场校正、坏像素校正和增益校正，校正手段可采用平板探测器厂方配套程序进行，也可以自行根据公知技术开发相应程序进行。

上述方法第3步中采用的射束衰减网格需添加在X射线源和被检测物体之间，如图3所示。射束衰减网格的构成方法是在低吸收系数的薄板上嵌入均匀分布的金属小球阵列。对射束衰减网格的一般要求是：在满足X射线能穿透金属小球和被检测物体的前提下，金属小球直径越小越好，并使金属小球之间具有适当的间距，一般为金属小球直径的1～2倍。显然，该网格比Beam Stop网格易于实现，且能满足在中高能电压下扫描的需求。实际中比较容易制作的射束衰减网格的构成方法是在低吸收系数的有机玻璃板上嵌入均匀分布的钢珠阵列。

射束衰减网格的放置，有如下两种可供选择：一种是放在被检测物体的前方(位置1)；另一种是放在被检测物体的后方(位置2)。对于位置1，X射线先透过射束衰减网格的金属小球，被小球吸收后的X射线再透过被检测物体。对于位置2，X射线先透过被检测物体，被检测物体吸收后的X射线再透过射束衰减网格金属小球。无论哪种位背，射束衰减网格本身也产生散射，本发明提出的散射校正方法假设：在相同的扫描条件下，射束衰减网格加上被检测物体的扫描和仅有被检测物体的扫描到达探测器的散射光子数不变。对于位置1射束衰减网格离探测器较远，且产生的散射部分被物体吸收，对本发明提出的校正方法影响较小，因此，本发明选择射束衰减网格置于被检测物体的前方。

上述方法第6步中射束衰减网格金属小球中心投影位置的计算方法为：对射束衰减网格投影图像经过二值化、轮廓提取、轮廓跟踪得到每个小球的投影轮廓，对这些轮廓采用圆拟合可计算出小球投影的近似圆形区域的圆心位置和直径。由于采用锥束CT，小球投影轮廓不是准确的圆形区域，因此，再以该圆心为中心、直径为边长的正方形区域内寻找投影图像灰度最小的位置作为小球中心的投影位置。

上述方法第7步中散射场分布图像的计算方法为：首先定义第4步中的扫描为扫描1(图4(a))，定义第5步中的扫描为扫描2(图4(b))，并假设射线源发出的X射线经预过滤后近似为单能射线，然后由扫描1和扫描2采用以下方法计算射束衰减网格中一个金属小球中心投影位置处(通常是平板探测器中一个接收象元，设该象元为G)的散射值。

对扫描1，设初始入射X射线强度为I₀，I₀透过小球后的透射强度为I₁，I₀透过小球并透过物体被象元G接收的透射强度为I₂，象元G接收到的所有散射射线强度为S，象元G接收到的总射线强度为C₂，显然C₂＝I₂+S。根据Beer定理，以小球为研究对象，有

$I_1=I_0{e}^{-u_{1}d}---\left(1\right)$

以被检测物体为研究对象，有

$I_2=C_2-S=I_1{e}^{-u_{2}l}---\left(2\right)$

对扫描2，设初始入射X射线强度仍为I₀，象元G接收到的所有散射射线强度仍为S，I₀透过被检测物体后的透射强度为I₃，象元G接收到的总射线强度为C₁，显然C₁＝I₃+S。以被检测物体为研究对象，有

$I_3=C_1-S=I_0{e}^{-u_{2}l}---\left(3\right)$

以上各式中，d为小球直径，l为射线在象元G位置穿越被检测物体的厚度，μ₁为小球的线性衰减系数，μ₂为被检测物体的线性衰减系数。

式(2)除以式(3)，得

$\frac{I_1}{I_0}=\frac{C_2-S}{C_1-S}---\left(4\right)$

将式(1)代入式(4)，有

$I_1=I_0{e}^{-u_{1}d}=I_0\frac{C_2-S}{C_1-S}---\left(5\right)$

则 $S=\frac{C_2-C_1{e}^{-u_{1}d}}{1-e^{-u_{1}d}}---\left(6\right)$

一般可认为平板探测器象元的输出灰度值与其所接收到的X射线强度成正比，因此可用扫描1时象元G的灰度代替C₂，用扫描2时象元G的灰度代替C₁。根据上式可计算出一个小球中心投影位置上象元G处的散射值S，对其余小球中心投影位置象元进行完全相似的计算，可以得到覆盖整个成像区域的很多象元处的散射值，然后采用双三次样条对其进行插值拟合得到整幅投影图像的散射场分布。

本发明的有益效果是：提出了一种新颖的锥束CT散射校正方法，该方法克服了基于BeamStop网格的散射校正方法在工业CT应用中难以满足较高扫描电压下铅球尺寸需求的缺点，实用于从低能到中、高能的锥束CT系统散射校正。对经散射校正后的投影图像进行锥束CT重建，所得切片中的杯状伪影相对于原切片明显减少，并可提高重建切片图像的对比度，改善图像质量，是锥束CT系统的一种简单、有效的散射校正方法。

下面结合附图和实施例对本发明进一步说明。

附图说明

图1是锥束CT扫描成像示意图；

图2是基于射束衰减网格的散射校正方法流程图；

图3是射束吸收网格校正方法实验装置图；

图4是透射散射X射线照射路径示意图；

图5是圆柱检测零件探测器输出原始投影、校正投影和散射强度之间的线性灰度关系图；

图6是圆柱检测零件未散射校正重建切片和散射校正后重建切片相同位置的线性灰度比较图。

具体实施方式

对一个直径为40mm，材料为铝质的圆柱检测零件进行锥束CT扫描，X射线源采用德国YXLON Y.TU 450-D02，平板探测器采用美国Varian PaxScan2520，应用本发明方法对其进行散射校正，如图2所示，执行以下步骤：

1.根据被检测物体确定对其进行锥束CT扫描的投影放大比为1.18，平板探测器全分辨率采集速度3幅/秒，扫描电压280KV，扫描电流0.25mA，这些参数在以下涉及投影图像数据采集的步骤中保持不变；

2.不放置任何物体，通过平板探测器采集6幅空气投影图像，取其平均值作为待用的空气投影图像；

3.将射束衰减网格置于X射线源和旋转工作台之间，通过平板探测器采集6幅射束衰减网格的投影图像，取其平均值作为待用的射束衰减网格投影图像，射束衰减网格由低吸收系数的有机玻璃板填上均匀分布的钢珠阵列组成，钢珠直径为3mm，相邻的两个小球中心间隔为5mm；

4.将被检测物体置于旋转工作台上，进行圆周扫描，扫描过程中旋转工作台连续旋转，通过平板探测器采集带射束衰减网格和被检测物体的投影图像集I，共360幅；

5.移除射束衰减网格，仅对被检测物体进行与步骤4完全相同的圆周扫描，其中需要特别注意的是被检测物体的起始扫描投影位置与步骤4的起始扫描投影位置相同，通过平板探测器采集被检测物体的投影图像集II，投影图像集II与投影图像集I的投影幅数相同，且零件各个投影位置对应相同；

6.通过步骤3得到的射束衰减网格投影图像计算射束衰减网格中每个金属小球中心的投影位置，以其中一个小球为例，计算得到的中心位置为(918，1023)；

7.通过金属小球中心的投影位置和步骤2、4、5得到的投影图像，采用射束衰减网格校正方法，计算与投影图像集I中的投影图像一一对应的散射场分布图像。以步骤(6)中计算得到的(918，1023)位置为例，其上空气投影图像值为2361.2417，射束衰减网格投影图像上的图像值为1397.1904，由式(1)计算得到小球的透射系数 $e^{-u_{1}d}=$ $0.5917.$ 在起始投影方位，即进行第一幅投影时，投影图像集II中第一幅投影(918，1023)位置的灰度值为C1＝642.9714，投影图像集I中第一幅投影(918，1023)位置的灰度值为C2＝509.4459，则该位置上的散射值为 $S=\frac{C_2-C_1{e}^{-u_{1}d}}{1-e^{-u_{1}d}}=315.9368.$ 然后采用双三次样条对所有金属小球中心的投影位置上的散射值进行插值拟合得到该投影方位上整幅投影图像的散射场分布，如点(895，1049)上的散射值为335.2801；

8.将投影图像集I中的各幅投影图像象素灰度减去与其一一对应的散射场分布图像象素灰度，得到散射校正后的投影图像集III，如(918，1023)位置上散射校正后的图像值为642.9714-315.9368＝327.0346；

9.由投影图像集III通过滤波反投影重建算法重建出散射校正后的序列切片图像。

对被检测零件进行圆轨迹锥束CT扫描时，从Varian公司的PaxScan2520平板探测器获取的投影图像数据已采用平板探测器厂方配套程序进行必要的暗场校正、坏像素校正和增益校正。

图5显示了圆柱检测零件探测器输出原始投影图像(P+S)、散射校正后投影图像(P)和散射场图像(S)在相同位置上的线性灰度分布关系，从图中可明显看出，经散射校正后的投影图像灰度值得到一定程度的降低。

图6显示了圆柱检测零件未校正的重建切片图像和散射校正后的重建切片图像在位置上的线性灰度分布情况比较。很明显，通过本发明的散射校正方法，圆柱检测零件的杯状伪影得到了较好的改进，同时图像的对比度也得到了一定提高。

Claims (5)

1.一种锥束CT系统的散射测定和校正方法，其特征在于包括下述步骤：

(a)根据被检测物体确定对其进行锥束CT扫描的投影放大比、平板探测器的数据采集方式、X射线源的电压与电流等扫描参数，这些参数在以下涉及投影图像数据采集的步骤中保持不变；

(b)不放置任何物体，通过平板探测器采集若干幅空气投影图像，取其平均值作为待用的空气投影图像；

(c)将射束衰减网格置于X射线源和旋转工作台之间，通过平板探测器采集若干幅射束衰减网格的投影图像，取其平均值作为待用的射束衰减网格投影图像；

(d)将被检测物体置于旋转工作台上，进行圆周扫描，扫描过程中旋转工作台采用步进旋转或连续旋转，通过平板探测器采集带射束衰减网格和被检测物体的投影图像集I；

(e)移除射束衰减网格，仅对被检测物体进行与步骤d完全相同的圆周扫描，被检测物体的起始扫描投影位置与步骤d的起始扫描投影位置相同，通过平板探测器采集被检测物体的投影图像集II，投影图像集II与投影图像集I的投影幅数相同，且零件各个投影位置对应相同；

(f)通过步骤c得到的射束衰减网格投影图像计算射束衰减网格中每个金属小球中心的投影位置；

(g)通过金属小球中心的投影位置和步骤b、d、e得到的投影图像，采用射束衰减网格校正方法，计算与投影图像集I中的投影图像一一对应的散射场分布图像；

(h)将投影图像集I中的各幅投影图像象素灰度减去与其一一对应的散射场分布图像象素灰度，得到散射校正后的投影图像集III；

(i)由投影图像集III通过滤波反投影重建算法重建出散射校正后的序列切片图像。

2.根据权利要求1的一种锥束CT系统的散射测定和校正方法，其特征在于：上述各步骤采集的投影图像在获取过程中已经过暗场校正、坏像素校正和增益校正。

3.根据权利要求1的一种锥束CT系统的散射测定和校正方法，其特征在于：上述步骤c中采用的射束衰减网格需添加在X射线源和被检测物体之间，射束衰减网格的构成方法是在低吸收系数的薄板中填上均匀分布的金属小球阵列。

4.根据权利要求1的一种锥束CT系统的散射测定和校正方法，其特征在于：上述步骤f中对射束衰减网格投影图像经过二值化、轮廓提取、轮廓跟踪得到每个小球的投影轮廓，对这些轮廓采用圆拟合可计算出小球投影的近似圆形区域的圆心位置和直径；再以该圆心为中心、直径为边长的正方形区域内寻找投影图像灰度最小的位置作为小球中心的投影位置。

5.根据权利要求1的一种锥束CT系统的散射测定和校正方法，其特征在于上述步骤g中散射场分布图像的计算方法为：

定义步骤d中的扫描为扫描1，定义步骤e中的扫描为扫描2，假设射线源发出的X射线经预过滤后近似为单能射线；

对扫描1，设初始入射X射线强度为I₀，I₀透过小球后的透射强度为I₁，I₀透过小球并透过物体被象元G接收的透射强度为I₂，象元G接收到的所有散射射线强度为S，象元G接收到的总射线强度为C₂，C₂＝I₂+S，根据Beer定理，以小球为研究对象，有

$I_1=I_0{e}^{-u_{1}d}---\left(1\right)$

以被检测物体为研究对象，有

$I_2=C_2-S=I_1{e}^{-u_{2}l}---\left(2\right)$

对扫描2，设初始入射X射线强度仍为I₀，象元G接收到的所有散射射线强度仍为S，I₀透过被检测物体后的透射强度为I₃，象元G接收到的总射线强度为C₁，C₁＝I₃+S，以被检测物体为研究对象，有

$I_3=C_1-S=I_0{e}^{-u_{2}l}---\left(3\right)$

以上各式中，d为小球直径，l为射线在象元G位置穿越被检测物体的厚度，μ₁为小球的线性衰减系数，μ₂为被检测物体的线性衰减系数；

式(2)除以式(3)，得

$\frac{I_1}{I_0}=\frac{C_2-S}{C_1-S}---\left(4\right)$

将式(1)代入式(4)，有

$I_1=I_0{e}^{-u_{1}d}=I_0\frac{C_2-S}{C_1-S}---\left(5\right)$

则 $S=\frac{C_2-C_1{e}^{-u_{1}d}}{1-e^{-u_{1}d}}---\left(6\right)$

根据上式计算出一个小球中心投影位置上象元G处的散射值S，对其余小球中心投影位置象元进行完全相似的计算，可以得到覆盖整个成像区域的很多象元处的散射值，然后采用双三次样条对其进行插值拟合得到整幅投影图像的散射场分布。

Images