一种考虑散射光子影响的X射线能谱估计方法

一种考虑散射光子影响的x射线能谱估计方法
技术领域
1.本发明涉及x射线ct成像技术领域，特别是一种考虑散射光子影响的x射线能谱估计方法。


背景技术：

2.x射线计算机断层(x ray computed tomography，简称x射线ct)成像技术能够无损地展示被测物体的内部结构信息，在医学诊疗、工业缺陷检测、制造误差分析以及逆向工程等方面得到越来越广泛的应用。ct成像的准确性依赖于x射线多色正投影建模的准确性，而x射线能谱是影响x射线多色正投影建模准确性的重要因素。x射线能谱在能谱ct图像重建、ct图像硬化校正等应用中扮演着重要角色。ct系统中x射线源产生的x射线流强很大，x射线能谱难以直接利用x射线光谱仪等设备进行测量；散射大小与被测物体的属性以及ct系统所处的环境相关，难以准确估计，且散射光子对能谱估计的准确性有一定影响。因此，想要对x射线多色正投影过程进行精确建模，需要对x射线能谱、散射进行准确估计。


技术实现要素：

3.本发明实施例的目的在于提供一种考虑散射光子影响的x射线能谱估计方法，来克服或至少减轻现有技术的上述缺陷中的至少一个。
4.为实现上述目的，本发明实施例提供一种考虑散射光子影响的x射线能谱估计方法，包括：
5.步骤1，设置散射光子信号为常数，构建关于x射线能谱和散射常数的非线性方程组：
[0006][0007]
其中，p
i
是第i条射线路径下获得的x射线多色投影值；i表示射线路径的编号，取值范围从1开始，最大值为探测器的单元个数与物体扫描角度数的积；为x射线路径的指标集；s(e)＝(s1,s2,
…
，s
m
)是未知的x射线能谱的离散形式，其中，m表示x射线能谱的最大能量；δ是能谱离散间隔；φ(e)＝(φ1,φ2,
…
,φ
m
)是被测物体的质量衰减系数的离散形式；r
i
＝(r
i1
,r
i2
,
…
,r
ij
)是第i条射线路径的投影向量，r
ij
表示第j个像素对第i条射线路径的投影贡献；f＝(f1,f2,
…
,f
j
)
t
表示离散化的图像，f
j
为图像f在第j个像素上的采样值；
[0008]
步骤2，设置所述非线性方程组中的x射线能谱s(e)和所述散射常数sc的初始值s(e)
(0)
和sc
(0)
；其中，所述s(e)
(0)
通过x射线能谱软件仿真预设电压下的x射线能谱得到，所述sc
(0)
为范围0.001
‑
0.1内的预设值；
[0009]
步骤3，在每一次迭代点(s(e)
(n)
,sc
(n)
)处将所述非线性方程组(1)进行一阶泰勒展开，得到线性方程组为：
[0010][0011]
其中：
[0012][0013][0014][0015]
其中，s(e)
(n)
和sc
(n)
分别表示第n次迭代结束后s(e)和sc的值；
[0016]
步骤4，将线性方程组(2)式转化为形式ax＝b，利用em(expectation
‑
maximization，期望最大化)算法求解所述线性方程组ax＝b，得到所述x射线能谱和所述散射常数的估计值；
[0017]
其中，
[0018][0019]
x＝(x1,x2,
…
,x
n
)
t
＝(s(e),sc)
t
ꢀꢀꢀ
(7)
[0020][0021]
采用如下公式利用所述em算法求解所述线性方程组ax＝b：
[0022][0023]
其中，x表示待求的未知数，x
j
表示第j个未知数，x
j(k)
表示第k次迭代的x
j
的结果，m表示所述线性方程组中方程的个数，n表示所述线性方程组中未知数的个数，a
ij
表示第i个方程第j个未知数对应的系数，b
i
表示第i个方程的右端项的数值；
[0024]
步骤5，判断所述步骤4得到的所述x射线能谱和所述散射常数的估计值是否收敛，如果判断结果为否，执行下一次迭代，重复所述步骤3和4；当所述x射线能谱和所述散射常数的估计值收敛时，执行步骤6；
[0025]
步骤6，根据所述x射线能谱和所述散射常数的估计值得到最终的x射线能谱。
[0026]
可选的，所述散射常数sc的初始值为0.01。
[0027]
可选的，s(e)
(0)
＞0，sc
(0)
＞0，且：
[0028]
可选的，判断所述步骤4得到的所述x射线能谱和所述散射常数的估计值是否收敛包括：
[0029]
利用所述x射线能谱和所述散射常数的估计值按照所述公式(3)计算得到估计的多色投影值p
i(n)
；若p
i(n)
与测量得到的多色投影值p
i
之差的二范数的平方小于给定阈值ε1，则判断收敛，即：
[0030][0031]
其中ε1＞0。
[0032]
可选的，判断所述步骤4得到的所述x射线能谱和所述散射常数的估计值是否收敛包括：
[0033]
相邻两次迭代估计出的x射线能谱之差的二范数的平方小于给定阈值ε2，则判断
收敛，即：
[0034][0035]
其中ε2＞0。
[0036]
与现有技术相比，本发明至少具有下述优点：
[0037]
由于本发明方法在估计x射线能谱时考虑了低频的散射信息，提高了x射线能谱在高能段的准确性；在估计x射线能谱和散射值时，将非线性方程组线性化之后利用em算法进行逐步迭代更新，降低了em算法对初值的依赖程度；本发明方法对标定模体精度要求不高，采用单材质的圆柱体即可满足需要，简单易用。
附图说明
[0038]
图1示出本发明实施例提供的考虑散射光子影响的x射线能谱估计方法的流程图；
[0039]
图2a为本发明一实施例所提供的待估计能谱和所设置的初始值示意图；
[0040]
图2b为本发明一实施例所提供的用于测试的圆形铝块模体示意图；
[0041]
图3为本发明一实施例中采用考虑散射光子影响的x射线能谱估计方法得到的x射线能谱结果示意图；
[0042]
图4a为本发明另一实施例所提供的用于估计x射线能谱的直径为4cm的圆柱形铝块模体的照片；
[0043]
图4b为本发明另一实施例对图4a中铝块进行扫描后利用sart算法直接重建得到的重建图像；
[0044]
图4c为本发明另一实施例对图4b中重建图像做阈值分割和二值化操作后的图像；
[0045]
图4d为本发明另一实施例对图4c做ct正投影操作后得到的铝块模体交线长与多色投影对应关系的散点图；
[0046]
图4e为本发明另一实施例对图4d中同一交线长下对应的多色投影值做加权平均操作后得到的铝块模体交线长与多色投影对应关系曲线；
[0047]
图5a为本发明另一实施例所提供的待估计x射线能谱的初始值以及采用本发明方法所估计的x射线能谱结果示意图；
[0048]
图5b为本发明另一实施例利用图5a中估计出的x射线能谱计算的多色投影曲线与实际得到的多色投影曲线拟合程度示意图；
[0049]
图6a为本发明一示例所提供的用于硬化伪影去除实验的直径为8cm的圆柱形铝块模体的照片；
[0050]
图6b为本发明一示例对图6a中铝块进行扫描后利用sart算法直接重建得到的重建图像，也是硬化伪影校正前的图像；
[0051]
图6c为本发明一示例利用图5a中估计出的x射线能谱去硬化后得到的结果；
[0052]
图6d为本发明一示例对图6c图像取中心行得到的线密度曲线图像。
具体实施方式
[0053]
下面结合附图和实施例对本发明进行详细的描述。
[0054]
图1示出本发明实施例提供的考虑散射光子影响的x射线能谱估计方法的流程图。
如图1所示，本发明实施例提供的考虑散射光子影响的x射线能谱估计方法包括以下步骤：
[0055]
步骤110，设置散射光子信号为常数，结合x射线多色投影正过程的非线性模型，构建关于x射线能谱和散射常数的非线性方程组。
[0056]
其中，假定散射光子信号为低频信号，设置为常数sc，结合x射线多色投影正过程的非线性模型，构建的关于x射线能谱和散射常数的非线性方程组为：
[0057][0058]
其中，p
i
是第i条射线路径下获得的x射线多色投影值，其值通过测量得到。i表示射线路径的编号，其取值范围从1开始，最大值为探测器的单元个数乘以物体扫描角度数。为x射线路径的指标集，其值由ct系统确定，在扫描参数确定的情况下，的取值是固定的。s(e)＝(s1,s2,
…
,s
m
)是未知的x射线能谱的离散形式，其中，m表示x射线能谱的最大能量，其值由ct系统的电压决定，假设当前ct系统电压为140kv，那么m取值为140。δ是能谱离散间隔，取值固定，在本实施例中取值为1。φ(e)＝(φ1,φ2,
…
,φ
m
)是被测物体的质量衰减系数的离散形式，当被测物体确定时，该φ(e)的形式确定。r
i
＝(r
i1
,r
i2
,
…
,r
ij
)是第i条射线路径的投影向量，r
ij
表示第j个像素对第i条射线路径的投影贡献，r
i
和r
ij
的值由ct系统确定，在扫描参数确定的情况下，r
i
和r
ij
值是固定的。f＝(f1,f2,
…
,f
j
)
t
表示离散化的图像，f
j
为图像f在第j个像素上的采样值，在扫描参数确定的情况下，f＝(f1,f2,
…
,f
j
)
t
的值是固定的。
[0059]
通过上式(1)，建立x射线能谱和散射常数的非线性关系。
[0060]
步骤120，为非线性方程组中的x射线能谱和散射赋初值。
[0061]
在一个实施例中，令x射线能谱初值为s(e)
(0)
，散射初值为sc
(0)
；其中s(e)
(0)
＞0，sc
(0)
＞0，且
[0062]
通常可以利用x射线能谱开源软件spectrum gui生成相应电压下的仿真能谱作为s(e)
(0)
的初始值。s(e)
(0)
的初始值可以由经验值选取，大概取值范围在0.001
‑
0.1之间。
[0063]
步骤130，将非线性方程组进行一阶泰勒展开，得到相应的线性方程组。
[0064]
在一个实施例中，在每一次迭代点(s(e)
(n)
,sc
(n)
)处将非线性方程组(1)进行一阶泰勒展开，其中，s(e)
(n)
和sc
(n)
分别表示第n次迭代结束后s(e)和sc的值。经过泰勒展开后得到的线性方程组为：
[0065][0066]
其中：
[0067][0068][0069][0070]
步骤140，利用em(expectation
‑
maximization)算法求解上述线性方程组，更新x射线能谱和散射的估计值。
[0071]
在一个实施例中，将线性方程组(2)式转化为形式ax＝b，
[0072][0073]
x＝(x1,x2,
…
,x
n
)
t
＝(s(e),sc)
t
ꢀꢀꢀ
(7)
[0074][0075]
利用em算法求解线性方程组ax＝b，利用的公式为：
[0076][0077]
其中，x表示待求的未知数，x
j
表示第j个未知数，x
j(k)
表示第k次迭代的x
j
的结果，m表示线性方程组中方程的个数，n表示线性方程组中未知数的个数，a
ij
表示第i个方程第j个未知数对应的系数，b
i
表示第i个方程的右端项的数值。
[0078]
在一个实施例中，本步骤中按照(9)式求解线性方程组ax＝b得到s(e)
(n)
和sc
(n)
的值，利用求解得到的值更新x射线能谱和散射的估计值，然后将估计的x射线能谱与散射值进行归一化处理，即：
[0079][0080]
以约束散射值sc的大小。
[0081]
步骤150，判断估计得到的x射线能谱与散射值是否收敛，如果否，则重复步骤130至140，如果收敛，则执行步骤160。
[0082]
在一个实施例中，判断估计的x射线能谱与散射值是否收敛有两个条件，二者满足其一即可停止迭代。
[0083]
条件1，利用估计的x射线能谱与散射值按照(3)式计算得到估计的多色投影值p
i(n)
；若p
i(n)
与测量得到的多色投影值p
i
之差的二范数的平方小于给定阈值ε1，则认为已经收敛，即：
[0084][0085]
其中ε1＞0，本实施例中取值为0.001。
[0086]
条件2，相邻两次迭代估计出的x射线能谱之差的二范数的平方小于给定阈值ε2，则认为已经收敛，即：
[0087][0088]
其中ε2＞0，本实施例中取值为0.001。
[0089]
容易理解，还可以选择其他收敛条件，包括但不限于最大迭代次数，以及利用估计的x射线能谱生成估计的多色投影曲线与测量得到的多色投影曲线的视觉吻合程度等判断是否满足收敛条件。
[0090]
步骤160，根据收敛的x射线能谱与散射值得到最终的x射线能谱。
[0091]
下面通过一个具体实施例来说明本发明提供的考虑散射光子影响的x射线能谱估计方法的具体实现过程。
[0092]
图2a为本实施例所提供的待估计能谱和所设置的初始值示意图。其中，待估计能
谱利用x射线能谱开源软件spectrum gui生成，仿真了hiray 7射线源在140kv电压下的能谱，并对其进行了归一化处理，如图2a中虚线所示；能谱初始值同样利用spectrum gui生成，仿真了oxford series6000射线源在140kv电压下的能谱，并对其进行归一化处理，如图2a中实线所示。图2b为本实施例所提供的用于测试的圆形铝块模体示意图，铝的质量衰减系数可以从例如美国国家标准技术研究院(nist)网站获得。x射线能谱和物质的质量衰减系数的采样间隔均为1kv。
[0093]
利用待估计能谱对图2b的圆形铝块进行ct扫描，按照公式(1)计算得到720个角度的在1536个探测器单元下的投影数据p，即i＝720*1536，公式(1)中散射常数设置为sc＝0.01。扫描参数设置如下：射线源到转台中心的距离(sod)为500mm，射线源到探测器的距离(sdd)为1000mm。线探测器单元数目为1024个，每个探测器单元的尺寸为0.2mm。对获得的投影数据p添加初始入射强度为i0＝105的泊松噪声，噪声添加公式如下：
[0094][0095]
其中，p
noisy
为添加噪声后的投影数据，poissrnd为生成泊松分布随机数的函数。
[0096]
本示例中提供的考虑散射光子影响的x射线能谱估计方法的具体实现过程包括：
[0097]
1)采用上述初始值给上述步骤110公式(1)中的x射线能谱赋初值为s(e)
(0)
，散射初值设置为sc
(0)
＝0.001。
[0098]
2)假设经过了n次迭代，此时x射线能谱和散射的估计值为(s(e)
(n)
,sc
(n)
)。将公式(1)关于(s(e),sc)在(s(e)
(n)
,sc
(n)
)处进行一阶泰勒展开，得到线性方程组如上述步骤130中的(2)式；其中p
i(n)
、q
i(n)
和φ的计算公式参照上述公式(3)
‑
(5)。
[0099]
3)根据公式(6)
‑
(8)将线性方程组(2)对应为ax＝b的形式，利用迭代格式即公式(9)求解并更新x射线能谱与散射值，并对更新后的x射线能谱与散射值按照公式(10)进行归一化处理。
[0100]
4)若步骤3)中得到的x射线能谱与散射值收敛，则停止迭代，根据收敛的x射线能谱与散射值得到最终的x射线能谱；否则执行下一次迭代(即n+1次迭代)，重复步骤2)和3)中的操作，直到x射线能谱与散射值收敛。
[0101]
图3为本实施例中采用考虑散射光子影响的x射线能谱估计方法得到的x射线能谱结果图。如图3所示，本文算法可以从含有噪声且含有散射的投影数据中有效恢复出待估计的x射线能谱，具有一定的抗噪性。
[0102]
下面通过一个具体实施例来说明本发明在实际扫描过程中的运用及其在去除硬化伪影中的应用。
[0103]
图4a为本实施例所提供的用于估计x射线能谱的直径为4cm的圆柱形铝块模体的照片，铝的质量衰减系数从美国国家标准技术研究院(nist)网站获得。利用x射线能谱开源软件spectrum gui生成，仿真了oxford series6000射线源在140kv电压下的能谱，并对其进行归一化处理，作为初值s(e)
(0)
。x射线能谱和物质的质量衰减系数的采样间隔均为1kv。利用工业ct设备对图4a中的铝块进行扫描，扫描参数配置如下：射线源到转台中心的距离(sod)为315.6mm，射线源到探测器的距离(sdd)为754.0mm；射线源电压设置为140kv，电流设置为0.14ma，射线源前添加1.5mm的铝滤波片；线探测器单元数目为1920个，每个探测器单元的尺寸为0.2mm；共采集1440个角度的数据，每个角度曝光时间0.2s，最终得到大小为
1920
×
1440的投影数据p
[m]
。
[0104]
在运用本发明所提供的方法前，需要对扫描得到的投影数据p
[m]
进行处理。首先，利用sart算法直接对投影数据p
[m]
进行重建得到重建图像，如图4b所示；然后对重建图像用阈值分割处理后再进行二值化处理，得到的二值化图像如图4c所示；按照相同的扫描配置，对二值化的图像进行ct正投影操作，得到图4a模体的交线长数据；将交线长数据与多色投影数据进行配对，可以得到图4d所示的模体交线长与多色投影对应关系的散点图；对同一交线长下的多色投影数据值进行加权平均，每个交线长可以得到唯一的多色投影值p，最终得到模体交线长与多色投影对应关系曲线，如图4e所示。
[0105]
然后，采用本发明提供的考虑散射光子影响的x射线能谱估计方法进行能谱估计，具体实施步骤如下：
[0106]
1)采用上述初始值给公式(1)中x射线能谱赋初值为s(e)
(0)
，散射初值设置为sc
(0)
＝0.001。
[0107]
2)假设经过了n次迭代，此时x射线能谱和散射的估计值为(s(e)
(n)
,sc
(n)
)。将公式(1)关于(s(e),sc)在(s(e)
(n)
,sc
(n)
)处进行一阶泰勒展开，对应处理后的多色投影值p，得到线性方程组如步骤130中的公式(2)；其中p
i(n)
、q
i(n)
和φ的计算公式可参照公式(3)
‑
(5)。
[0108]
3)根据公式(6)
‑
(8)将线性方程组(2)对应为ax＝b的形式，利用迭代格式即公式(9)求解并更新x射线能谱与散射值，并对更新后的x射线能谱与散射值按照(10)式进行归一化处理。
[0109]
4)若步骤3)中得到的x射线能谱与散射值收敛，则停止迭代，根据收敛的x射线能谱与散射值得到最终的x射线能谱；否则执行下一次迭代(即n+1次迭代)，重复步骤2)和3)中的操作，直到x射线能谱与散射值收敛。
[0110]
图5a为本实施例中待估计x射线能谱的初始值以及采用本发明方法所得到的x射线能谱结果示意图。实采数据真实的x射线能谱未知，可以利用估计出的x射线能谱结果按照公式(1)生成相应的多色投影值，将其与实采得到的多色投影p
[m]
进行对比。图5b为利用图5a中估计出的x射线能谱计算的多色投影曲线与实际得到的多色投影曲线拟合程度示意图，可以看到，利用本发明提供的考虑散射光子影响的x射线能谱估计方法估计出的x射线能谱生成的多色投影曲线，可以很好的拟合实采数据，从而证明了本发明方法的可行性。
[0111]
下面再通过一示例进一步说明采用本发明提供的考虑散射光子影响的x射线能谱估计方法估计出的x射线能谱在去除硬化伪影中的应用。
[0112]
图6a为本示例所提供的用于硬化伪影去除实验的直径为8cm的圆柱形铝块模体的照片。采用扫描图5a所述模体时同样的扫描设置，对其进行扫描，利用sart算法直接重建得到的重建图像如图6b所示，图6b也是硬化伪影校正前的图像。利用图5a中估计出的x射线能谱去硬化后得到的结果如图6c所示，对图6c图像取中心行得到的线密度曲线如图6d所示。可以看到，图6c中圆心部分的硬化伪影有了明显改善，图6d原本因射线硬化而导致的曲线凹陷被修正。图6c与图6d说明了利用本发明提供的方法估计的x射线能谱一定程度上可以去除硬化伪影。
[0113]
本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。
[0114]
最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本
领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。