一种基于相位CT的材料分解与识别方法

摘要

本发明公开了一种基于相位CT的材料分解与识别方法，所述方法包括：步骤1，令m＝0,初始化第一基材图像f(m)和第二基材图像g(m)；步骤2，根据式(6)计算离散化吸收投影数据的第m轮迭代中的估计值和对应的残差以及离散化微分相位投影数据在第m轮迭代中的估计值和对应的残差步骤3，根据式(5)迭代求解第一基材图像f(m)在第m+1轮迭代中的值f(m+1)和第二基材图像g(m)在第m+1轮迭代中的值g(m+1)；步骤4，对f(m+1)，g(m+1)进行正则化操作；步骤5，若停止条件不满足，令m＝m+1，并返回步骤2。本发明方法联合吸收投影信息及微分相位投影信息进行基材分解，得到的基材图像噪声更小，后续计算的等效原子序数图像的噪声也比常规的更小。

权利要求

1.一种基于相位CT的材料分解与识别方法，其特征在于，包括：(m) (m)

步骤1，令m＝0,初始化第一基材图像为f 和第二基材图像为g ；

步骤2，根据式(6)计算离散化吸收投影数据 的第m轮迭代中的估计值 和对应的残差 以及离散化微分相位投影数据 在第m轮迭代中的估计值 和对应的残差(m+1)步骤3，根据式(5)迭代求解第一基材图像f在第m+1轮迭代中的值f 和第二基材图像(m+1)g在第m+1轮迭代中的值g ；

(m+1) (m+1)

步骤4，对f ，g 进行正则化操作；

步骤5，若停止条件不满足，令m＝m+1，并返回步骤2；

其中：

表示第一基材的fj和第二基材的gj在投影角度 下对第u条射线路径的贡献；

表示第一基材的fj在第m+1轮迭代中的值；

表示第一基材的fj在第m轮迭代中的值；

表示第二基材的gj在第m+1轮迭代中的值；

表示第二基材的gj在第m轮迭代中的值；

λ表示松弛因子；

表示离散化吸收投影数据 的第m轮迭代中的估计值；

表示离散化微分相位投影数据 在第m轮迭代中的估计值；

表示第m轮迭代中的吸收投影残差 的第u个分量；

表示第m轮迭代中的微分相位投影残差 的第u个分量。

2.如权利要求1所述的基于相位CT的材料分解与识别方法，其特征在于，采用式(7)表示的第二优化模型求解 的最优值式中， 表示 和当前相位投影残差 的二范数距离最小时的相位投影残差

3.如权利要求1或2所述的基于相位CT的材料分解与识别方法，其特征在于，离散化吸收投影数据 和离散化微分相位投影数据 表示式(3)：其中：

表示投影角度 下不同探测器接收到的离散化吸收投影数据，其中的元素 表示第u个探测器像素采集的吸收投影值，U表示探测器的总像素个数，τ表示向量的转置；

表示投影角度 下的投影矩阵，其中的J表示每幅图像的总像素个数；

τ

f＝(f1,f2...,fJ) 表示第一基材的分布函数图像f(x,y)的离散化图像，其中的元素fj为f(x,y)在被测样品上第j个像素上的采样值；

τ

g＝(g1,g2...,gJ) 表示第二基材的分布函数图像g(x,y)的离散化图像，其中的元素gj为g(x,y)在被测样品上第j个像素上的采样值；

表示投影角度 下不同探测器的离散化微分相位投影数据，其中的元素 表示第u个探测器像素采集的微分相位投影值。

4.如权利要求3所述的基于相位CT的材料分解与识别方法，其特征在于，通过求解如下式(4)提供的第一优化模型，得到最后的分解图像：其中：

*

f表示第一优化模型取得最小时第一基材的分布函数图像f；

*

g表示第一优化模型取得最小时第二基材的分布函数图像g；

R(·)是正则化算子。

5.如权利要求3所述的基于相位CT的材料分解与识别方法，其特征在于，还包括：(m) (m)

步骤6，根据基材图像f 和g ，利用公式(1)得到吸收图像和相位图像；

其中：

μ(x,y)表示被测样品的点(x,y)处的吸收图像；

δ(x,y)表示被测样品的点(x,y)处的相位图像；

μ1为第一基材的线性吸收图像；

μ2为第二基材的线性吸收图像；

δ1为第一基材的相位图像；

δ2为第二基材的相位图像；

f(x,y)为第一基材的分布函数图像；

g(x,y)为第二基材的分布函数图像；

步骤7，利用式(9)求解等效原子序数：其中：

Cp、CE、CZ是通过扫描与被测物体相近的模体的数据进行标定求解的参数常量；

r0是经典原子半径，h是普朗克常数，c是光速常数；

CKN(E)为Klein‑Nishina截面；

μ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的吸收图像；

δ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的相位图像；

ρe(x,y)表示被测样品的点(x,y)处的电子密度；

Z(x,y)表示被测样品的点(x,y)处的等效原子序数。

说明书

一种基于相位CT的材料分解与识别方法

技术领域

[0001] 本发明涉及材料识别技术领域，特别是关于一种基于相位CT的材料分解与识别方法。

背景技术

[0002] X射线成像技术被广泛应用在很多领域，例如生物医学、工业和安全检查等等。传统X射线吸收CT是基于不同物质或者组织之间吸收辐射的差异性进行成像，对于轻元素构成的材料(如肌肉、血管和乳腺等)成像不够理想。而X射线相位衬度CT是基于物质对穿过的X射线相位的改变而成像，通常其对轻物质的成像比吸收CT有较大的优势。X射线相衬成像方法里常用的是光栅相衬方法和衍射增强法，这两种方法都能同时提供吸收和相位信息。利用吸收信息和相位信息进行的物质分解和识别是相位CT的一个重要应用，在低能范围内，其效果比传统的利用双能CT进行物质分解和识别的方法更好。

[0003] 目前，常用的基于相位CT的材料分解和识别技术主要是：1)基于投影域的分解技术；2)基于重建图像域的分解技术。由于采集到的数据的相位信息是相位投影的微分，需要经过积分运算才能得到相位投影。这个积分过程容易引起整体的低频噪声，使得计算的相位投影不准确，进而导致在投影域的分解误差扩大，所以在相位CT的材料分解重建中，通常以图像域的分解为主。另外，这两种方法将重建和分解的步骤分离，得到的分解图像的质量并不是太好，会影响到后续的定量分析。因此，为提高分解图像的质量，发展同时利用吸收和相位信息的重建分解算法是有意义的。

[0004] 基于重建图像域的分解技术方法是基于吸收图像和相位图像可通过两种基材图像线性组合表示，在分别重建出吸收图像和相位图像后，选择所需要的两个基材，查找当前等效能量下两基材的吸收和相位信息作为线性组合的系数，然后求解得到两种基材的分布函数图像。

[0005] 但是，在图像域中进行一次性的分解，存在噪声放大问题。分解的结果会受到重建图像的影响，而重建图像的质量主要由初始的投影数据的噪声水平和重建算法的选取决定，没有联合吸收及相位的信息进行重建和分解，得到的结果欠佳。

发明内容

[0006] 本发明的目的在于提供一种基于相位CT的材料分解与识别方法来克服或至少减轻现有技术的上述缺陷中的至少一个。

[0007] 为实现上述目的，本发明提供一种基于相位CT的材料分解与识别方法，所述方法包括：

[0008] 步骤1，令m＝0,初始化第一基材图像f(m)和第二基材图像g(m)；

[0009] 步骤2，根据式(6)计算离散化吸收投影数据 的第m轮迭代中的估计值 和对应的残差 以及离散化微分相位投影数据 在第m轮迭代中的估计值 和对应的残差

[0010] 步骤3，根据式(5)迭代求解第一基材图像f(m)在第m+1轮迭代中的值f(m+1)和第二基(m) (m+1)材图像g 在第m+1轮迭代中的值g ；

[0011] 步骤4，对f(m+1)，g(m+1)进行正则化操作；

[0012] 步骤5，若停止条件不满足，令m＝m+1，并返回步骤2；

[0013]

[0014]

[0015] 其中：

[0016] 表示第一基材的fj和第二基材的gj在投影角度 下对第u条射线路径的贡献；

[0017] 表示第一基材的fj在第m+1轮迭代中的值；

[0018] 表示第一基材的fj在第m轮迭代中的值；

[0019] 表示第二基材的gj在第m+1轮迭代中的值；

[0020] 表示第二基材的gj在第m轮迭代中的值；

[0021] λ表示松弛因子；

[0022] 表示离散化吸收投影数据 的第m轮迭代中的估计值；

[0023] 表示离散化微分相位投影数据 在第m轮迭代中的估计值；

[0024] 表示第m轮迭代中的吸收投影残差 的第u个分量；

[0025] 表示第m轮迭代中的微分相位投影残差 的第u个分量。

[0026] 进一步地，采用式(7)表示的第二优化模型求解 的最优值

[0027]

[0028] 式中， 表示 和当前相位投影残差 的二范数距离最小时的相位投影残差

[0029] 进一步地，离散化吸收投影数据 和离散化微分相位投影数据 表示式(3)：

[0030]

[0031] 其中：

[0032] 表示投影角度 下不同探测器接收到的离散化吸收投影数据，其中的元素 表示第u个探测器像素采集的吸收投影值，U表示探测器的总像素个数，τ表示向量的转置；

[0033] 表示投影角度 下的投影矩阵，其中的J表示每幅图像的总像素个数；

[0034] f＝(f1,f2...,fJ)τ表示第一基材的分布函数图像f(x,y)的离散化图像，其中的元素fj为f(x,y)在被测样品上第j个像素上的采样值；

[0035] g＝(g1,g2...,gJ)τ表示第二基材的分布函数图像g(x,y)的离散化图像，其中的元素gj为g(x,y)在被测样品上第j个像素上的采样值；

[0036] 表示投影角度 下不同探测器的离散化微分相位投影数据，其中的元素 表示第u个探测器像素采集的微分相位投影值。

[0037] 进一步地，通过求解如下式(4)提供的第一优化模型，得到最后的分解图像：

[0038]

[0039] 其中：

[0040] f*表示第一优化模型取得最小时第一基材的分布函数图像f；

[0041] g*表示第一优化模型取得最小时第二基材的分布函数图像g；

[0042] R(·)是正则化算子。

[0043] 进一步地，所述方法还包括：

[0044] 步骤6，根据基材图像f(m)和g(m)，利用公式(1)得到吸收图像和相位图像；

[0045]

[0046] 其中：

[0047] μ(x,y)表示被测样品的点(x,y)处的吸收图像；

[0048] δ(x,y)表示被测样品的点(x,y)处的相位图像；

[0049] μ1为第一基材的线性吸收图像；

[0050] μ2为第二基材的线性吸收图像；

[0051] δ1为第一基材的相位图像；

[0052] δ2为第二基材的相位图像；

[0053] f(x,y)为第一基材的分布函数图像；

[0054] g(x,y)为第二基材的分布函数图像；

[0055] 步骤7，利用式(9)求解等效原子序数：

[0056]

[0057] 其中：

[0058] Cp、CE、CZ是通过扫描与被测物体相近的模体的数据进行标定求解的参数常量；

[0059] r0是经典原子半径，h是普朗克常数，c是光速常数；

[0060] CKN(E)为Klein‑Nishina截面；

[0061] μ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的吸收图像；

[0062] δ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的相位图像；

[0063] ρe(x,y)表示被测样品的点(x,y)处的电子密度；

[0064] Z(x,y)表示被测样品的点(x,y)处的等效原子序数。

[0065] 由于本发明根据分解基材图像的需要，建立了直接由微分相位和吸收投影数据重建基材图像的、能抑制噪声的优化成像模型，结合图像迭代重建方法，提出了一种利用微分相位和吸收的物质分解方法，该方法建立在优化模型的基础上，采用的迭代算法可直接从吸收和微分相位的投影中重建基材成分，进一步地计算可得到等效原子序数，具有抑制噪声和分解准确的优点。

附图说明

[0066] 图1为Image‑based方法和本发明方法分解的各种基材图像结果对比示意图。

[0067] 图2为使用双能采集和本发明方法计算的等效原子序数图像对比示意图。

具体实施方式

[0068] 下面结合附图和实施例对本发明进行详细的描述。

[0069] 根据材质分解方法，在单能的情况下，被测样品的线性吸收图像可以用两种基材的线性吸收图像近似线性表示。同理，被测样品的折射率实部衰减率也可以用两种基材的折射率实部衰减率近似线性表示。于是，可以列出下列方程组(1)：

[0070]

[0071] 其中：

[0072] μ(x,y)表示被测样品的点(x,y)处的吸收图像。

[0073] δ(x,y)表示被测样品的点(x,y)处的折射率实部衰减率系数，即文中的相位图像。

[0074] μ1为第一基材的线性吸收图像。

[0075] μ2为第二基材的线性吸收图像。

[0076] δ1为第一基材的相位图像。

[0077] δ2为第二基材的相位图像。

[0078] f(x,y)为第一基材的分布函数图像。

[0079] g(x,y)为第二基材的分布函数图像。

[0080] (x,y)表示样品坐标系中任一点的坐标，样品坐标系通常选取样品的几何中心为原点，X轴和Y轴相互垂直，其方向可根据实际需求进行选取。

[0081] 基于光栅干涉仪的X射线微分相位成像通过光栅步进的采集方式，可以分离出每个投影角度 下样品对X射线的吸收和微分相移信息表示为式(2)：

[0082]

[0083] 其中：

[0084] 为 坐标下接收到的吸收投影数据。

[0085] 为 坐标下微分相位投影数据。

[0086] 表示测量坐标系中任一点的坐标。

[0087] 将式(2)可以离散化为式(3)：

[0088]

[0089] 其中：

[0090] 表示投影角度 下不同探测器接收到的离散化吸收投影数据，其中的元素 表示第u个探测器像素采集的吸收投影值，U表示探测器的总像素个数，τ表示向量的转置。

[0091] 表示投影角度 下的投影矩阵，其中的元素 表示fj和gj在投影角度 下对第u条射线路径的贡献，J表示每幅图像的总像素个数。

[0092] f＝(f1,f2...,fJ)τ表示f(x,y)的离散化图像，其中的元素fj为f(x,y)在被测样品上第j个像素上的采样值。

[0093] g＝(g1,g2...,gJ)τ表示g(x,y)的离散化图像，其中的元素gj为g(x,y)在被测样品上第j个像素上的采样值。

[0094] 表示投影角度 下不同探测器的离散化微分相位投影数据，其中的元素 表示第u个探测器像素采集的微分相位投影值；D是微分算子 离散后的矩阵形式。

[0095] 本发明实施例解决的是根据离散化吸收投影数据 和离散化微分相位投影数据 重建两种基材的分布函数图像f和g。

[0096] 本发明实施例是通过求解如下式(4)提供的第一优化模型，得到最后的分解图像：

[0097]

[0098] 其中：

[0099] f*表示第一优化模型取得最小时第一基材的分布函数图像；

[0100] g*表示第一优化模型取得最小时第二基材的分布函数图像；

[0101] R(·)是正则化算子，如TV、L0或图像滤波操作；

[0102] ||·||2表示二范数的平方。

[0103] 也可以采用其它的模型得到最后分解的图像，比如：采用其它范数代替上述第一优化模型的二范数的平方，并通过在吸收项的范数前面加一个参数，控制吸收项和相位项的比重，使得由f和g计算的吸收项和微分相位项尽可能接近探测器探测到的吸收值和微分相位值(保真)。

[0104] 在式(4)的保真项中，应用如下式(5)表示的迭代步骤：

[0105]

[0106] 其中：

[0107]

[0108] 表示第一基材的fj在第m+1轮迭代中的值。

[0109] 表示第一基材的fj在第m轮迭代中的值。

[0110] 表示第二基材的gj在第m+1轮迭代中的值。

[0111] 表示第二基材的gj在第m轮迭代中的值。

[0112] λ表示松弛因子，其取值在(0，1)之间，可抑制噪声的影响。

[0113] 表示J个 的总和。

[0114] 表示U个 的总和。

[0115] 表示离散化吸收投影数据 的第m轮迭代中的估计值。

[0116] 表示离散化微分相位投影数据 在第m轮迭代中的估计值。

[0117] 表示第m轮迭代中的吸收投影残差 的第u个分量。

[0118] 表示第m轮迭代中的微分相位投影残差 的第u个分量。

[0119] 计算相位投影残差 需要利用微分算子的逆，即D‑1，而该算子对噪声有放大作用。为了克服噪声对重建图像的影响，可以采用式(7)表示的第二优化模型求解：

[0120]

[0121] 式中， 表示 和当前相位投影残差 的二范数距离最小时的相位投影残差

[0122] 在另一个实施例中，可以采用式(8)表示的第三优化模型求解：

[0123]

[0124] 式中， 表示D(e)和当前相位投影残差 距离的二范数与e的正则化最小时的相位投影残差e。

[0125] 在又一个实施例中，还可以采用距离范数和对变量e进行正则化约束，来确定相位投影残差

[0126] 实际迭代中采用的是固定步数的梯度下降法来进行求解。当然，也可以直接求解或其它求解方式。

[0127] 对于第一模型的正则化项，本实施例使用的策略是在对所有投影数据使用式(5)(m+1) (m+1)更新完f ，g 后，对其进行一次正则化操作。正则化算子主要选取TV算子，从而可以较好地抑制噪声的影响。

[0128] 综上，本发明实施例的迭代步骤如下：

[0129] 步骤1，初始化：f(0)＝0，g(0)＝0，m＝0。

[0130] 步骤2，根据式(6)计算吸收投影的估计值 和对应的残差

[0131] 步骤3，根据式(6)计算相位投影的估计值 由式(7)计算其残差

[0132] 步骤4，根据式(5)迭代求解f(m+1)，g(m+1)。

[0133] 步骤5，对f(m+1)，g(m+1)进行正则化操作。

[0134] 步骤6，若停止条件不满足，令m＝m+1，并返回步骤2。其中，停止条件可以保真项的二范数是否小于某个值，或者设置固定的迭代步数。

[0135] 步骤7，返回f(m)和g(m)。

[0136] 相位CT及本发明方法的一个应用是对物质进行定量分析，利用上述方法得到基材(m) (m)图像f 和g ，代入式(1)计算吸收图像μ(x,y)和相位图像δ(x,y)，即式(9)中的μ(E,x,y)和δ(E,x,y)，再通过比如查表的方式获得等效能量E。接着，利用以下公式(9)可求解得到等效原子序数：

[0137]

[0138] 其中：

[0139] Cp、CE、CZ是通过扫描与被测物体相近的模体的数据进行标定求解的参数常量，比如：根据被测物体，选取基材相近的已知模体，那么，Z(x,y)、μ(E,x,y)和δ(E,x,y)均为已知量，再查表得到等效能量E，最后可以拟合出Cp、CE、CZ动具体数值。

[0140] r0是经典原子半径，h是普朗克常数，c是光速常数，各具体数值可通过查表获得。

[0141] CKN(E)为Klein‑Nishina截面，通常使用公式(10)计算：

[0142]

[0143] 其中a＝E/511keV是相对电子质量，r0是经典原子半径。

[0144] μ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的吸收图像。

[0145] δ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的相位图像。

[0146] ρe(x,y)表示被测样品的点(x,y)处的电子密度。

[0147] Z(x,y)表示被测样品的点(x,y)处的等效原子序数。

[0148] 此流程得到的吸收图像和相位图像是同时利用了吸收信息和相位信息，而并非原来单纯的由吸收信息得到吸收图像、相位信息得到相位图像。联合两种信息可减少原来单一信息直接计算吸收图像和相位图像的噪声影响，计算的等效原子序数图像的噪声也会更小。

[0149] 下面是一组实验数据：在上海光源BL13W实验站上，使用基于光栅干涉仪的X射线微分相位成像装置上采集。光栅干涉仪由周期为2.396um的π/2相位光栅和周期为2.4um的吸收光栅构成，两块光栅之间的距离为46.38mm。实验所用的光子能量为20keV，采用有效像素阵列为2048×600的sCMOS探测器，像素尺寸为6.5μm。采集的数据是使用8步步进方法，即对相位光栅在一个周期内移动8步分别采集数据，完成一个角度的采集，而在180度内等间隔采集540个角度的数据。最后根据步进方法提取吸收和微分相位公式计算吸收投影数据和微分相位投影数据。

[0150] 如图1所示，a是实物照片，模体包含四种成分：低密度聚乙烯(LDPE)、聚甲基丙烯酸甲酯(PMMA)、聚四氟乙烯(PTFE)和水。这几种成分被放在一个试管中的外直径10.7mm的聚乙烯塑料试管，PTFE，LDPE和PMMA的直径分别是2.0mm，4.0mm和5.6mm。

[0151] 在图1中的b的吸收图像断层重建图中，圆环代表聚乙烯塑料容器，其中放置三种成分不同的小球，即圆环中间的三个灰度不同的圆面积代表三种成分的小球，左边最暗的是LDPE小球,下边最亮的是PTFE小球,而右上方的是PMMA小球，其余部分是水。

[0152] 在图1中的c的折射率实部衰减率断层重建图中，LDPE和水无法区分开。这表明，利用相位信息重建的断层图像并不一定比利用吸收信息重建的断层图像的衬度高。

[0153] 图1中的第二第三行是两种方法分解的PTFE基和水基结果，图1中的e和f分别是ImageBase和PAMD‑SART方法分解的PTFE基图像，图1中的h和i分别是ImageBase和PAMD‑SART方法分解的水基图像，图1中的d和g分别是两种方法下PTFE基和水基结果的剖线对比图。

[0154] 定量分析的结果如图2所示，与相位CT对比的是使用双能采集计算的结果，在图上可看到相位CT计算的等效原子序数图像是要优于双能的。

[0155] 图2，分解结果对比图，其中的a20keV下计算的原子序数断层图像，其中的b是12keV的虚拟吸收图像断层图像，其中的c是20keV的虚拟吸收图像断层图像，其中的d是双能计算的原子序数断层图像，其中的e是12keV的真实吸收图像断层图像，其中的f是20keV的真实吸收图像断层图像。

[0156] 基于投影域的分解技术，基于重建图像域的分解技术和DEXA都可以实现材料分解和识别。但前两种方法属于两步法，存在噪声放大的问题。DEXA当被测样品由弱吸收的低原子序数物质构成时,其物质区分能力往往难以满足要求。

[0157] 最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。