CN109875593B - X射线多谱分离成像方法、存储介质和装置 - Google Patents

Abstract

本发明公开一种X射线多谱分离成像方法、存储介质和装置，该方法包括：步骤11：设定被测物的成像电压范围[V1,V2]；步骤12：在[V1,V2]的范围内取N个不同的扫描电压，N大于被测物所包含的材质总数K，且任意两个扫描电压之间的差值大于第一预设值，在每个扫描电压下对被测物进行全角度扫描成像，获得符合成像灰度要求的N个投影图像序列；步骤13：基于N个投影图像序列，求解满足泊松极大似然原理的变能量多谱衰减表达式中的被测物的所有材质对应的投影厚度矩阵D；步骤14：基于D，计算X射线中每个窄能谱段的分离投影p_r，r＝1,2…R；步骤15：重建p_r得到每个窄能谱段的重建图像a_r。本发明的方法可以实现低成本、能谱分辨率高的多能谱CT成像。

Description

X射线多谱分离成像方法、存储介质和装置

技术领域

本发明涉及仪器领域，特别涉及一种X射线多谱分离成像方法、存储介质和装置。

背景技术

在X射线成像过程中，由于X射线的多谱特性以及物质的衰减系数通常随光子能量的增加而降低，射线在穿过物体时其光谱中心随着穿过厚度的增加逐渐向高能方向移动，即产生硬化效应。硬化效应破坏了衰减系数与投影数据间的线性关系，使得常规重建方法所得CT图像产生硬化伪影，出现“同影异物”和“同物异影”现象，不利于定量化分析，如新型材料的组分含量、矿石孔隙率等微观结构的定量化表征。

针对上述问题，出现了多谱CT成像技术的相关研究。典型的双能CT通过获取两个不同能谱的投影数据，能够实现对投影数据的非线性重建，它主要包括有投影双能分解和图像重建两个步骤，目前常见的双能CT***比较常用的重建算法是基于基物质分解模型的双物质分解算法，由于引入了能谱的信息，双能CT较传统CT有更好的物质区分能力。但双能CT只有两个能谱，区分度不高，于是人们提出了多谱CT成像技术。其中，应用和研究较多的是基于光子计数型探测器和基于同步辐射源的多能谱CT成像技术，光子计数型探测器可以实现多个窄能谱成像，缺点是其时间/空间分辨力较低、技术成本较高；同步辐射源具有良好的单色性，但其设备庞大，造价昂贵，难以在实际工程中应用。

如何实现低成本、能谱分辨率高的多能谱(多谱)CT成像技术是目前亟待解决的技术问题。

发明内容

有鉴于此，本发明提供一种X射线多谱分离成像方法、存储介质和装置，以解决如何实现低成本、能谱分辨率高的多能谱CT成像技术的问题。

本发明提供一种X射线多谱分离成像方法，该方法包括：

步骤11：设定被测物的成像电压范围[V1,V2]，当V∈[V1,V2]时，被测物的穿透率符合第一预设条件、且被测物的任意两种材质在V时的线性衰减系数的相对偏差满足预设要求对应的能量范围符合第二预设条件；

步骤12：在[V1,V2]的范围内取N个不同的扫描电压，N大于被测物所包含的材质总数K，且任意两个扫描电压之间的差值大于第一预设值，在每个扫描电压下对被测物进行全角度扫描成像，获得符合成像灰度要求的N个投影图像序列，全角度扫描包括T个不同的成像角度；

步骤13：基于N个投影图像序列，求解满足泊松极大似然原理的变能量多谱衰减表达式中的被测物的所有材质对应的投影厚度矩阵D；

步骤14：基于D，计算X射线中每个窄能谱段的分离投影p_r，r＝1,2…R，R为X射线的窄能谱段总数；

步骤15：重建p_r得到每个窄能谱段的重建图像a_r。

本发明还提供一种非瞬时计算机可读存储介质，非瞬时计算机可读存储介质存储指令，指令在由处理器执行时使得处理器执行上述的X射线多谱分离成像方法中的步骤。

本发明还提供一种X射线多谱分离成像装置，包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，处理器执行计算机程序时实现上述的X射线多谱分离成像方法中的步骤。

本发明的X射线多谱分离成像方法可以实现X射线投影图像盲分离，根据物体的材质、尺寸，选取合理的成像参数获取X射线投影图像序列；然后依据投影图像序列，求解满足泊松极大似然原理的变能量多谱衰减表达式中的被测物的所有材质对应的投影厚度矩阵D；再根据矩阵D实现多谱投影图像的投影分离。相比于基于硬件和其他多能谱重建的X射线CT成像技术，该方法从投影域实现了多谱X射线图像的能谱分离，无需对现有的高动态CT成像设备进行改进，具有成本低、能谱分辨率高等优点，对于推广定量化CT成像技术的应用有重要意义。

附图说明

图1为本发明X射线多谱分离成像方法的流程图；

图2为一个扫描电压获得的投影图像序列示意图；

图3为图2的子单元细分示意图；

图4为图1中步骤13的求解D的第一实施例；

图5为图1中步骤13的求解D的第二实施例；

图6为本发明X射线多谱分离成像装置的结构图。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面结合附图和具体实施例对本发明进行详细描述。

如图1所示，本发明提供一种X射线多谱分离成像方法，包括：

步骤11(S11)：设定被测物的成像电压范围[V1,V2]，当V∈[V1,V2]时，被测物的穿透率符合第一预设条件、且被测物的任意两种材质在V时的线性衰减系数的相对偏差满足预设要求对应的能量范围符合第二预设条件；

例如，第一预设条件：依据被测物的材质、尺寸、探测器成像动态范围，并避免欠曝光和过曝光现象造成信息的缺失，选取合适的电压区间[V₁,V₂]使物体的穿透率达到30％～70％或其他范围。

第二预设条件为：从被测物组分分析的角度出发，依据材质的线性衰减系数曲线，满足任意电压V∈[V₁,V₂]时任意第k₁、k₂种材质满足式(1)的能量范围Ω与V所产生的能量范围Ω₁的交集占Ω₁的60％以上，以保证各组分在选取能量范围上具有可区分性。

式中μ_k1(E)和μ_k2(E)分别表示第k₁、k₂种材质在能量E时的线性衰减系数。根据预设要求的严格与否，式(1)中的0.05也可以换成其他值。

步骤12(S12)：在[V1,V2]的范围内取N个不同的扫描电压，N大于被测物所包含的材质总数K，且任意两个扫描电压之间的差值大于第一预设值，在每个扫描电压下对被测物进行全角度扫描成像，获得符合成像灰度要求的N个投影图像序列，全角度扫描包括T个不同的成像角度；

N个不同的扫描电压可以是等间隔电压序列，也可以为非等间隔电压序列，无论是否等间隔，都需满足，N大于被测物所包含的材质总数K，且任意两个扫描电压之间的差值大于第一预设值，第一预设值根据穿过被测物的同一射线路径上射线强度与入射强度比值的最小区分要求设定，例如最小区分要求为比值的相对偏差大于2％。

同样地，T个成像角度也可以为等间隔或非等间隔，例如T个不同的成像角度包括：初始成像角度为θ₀、角度间隔为△θ的T个成像角度。T的取值由步骤13决定，通常T≥180均可满足步骤13求解D的要求。

每个扫描电压的每个角度成像时，通过调节X射线管的电流和探测器的积分时间使采集到的X射线图像灰度值在探测器的有效成像范围内，即成像灰度符合要求。

步骤13(S13)：基于N个投影图像序列，求解满足泊松极大似然原理的变能量多谱衰减表达式中的被测物的所有材质对应的投影厚度矩阵D；

步骤14(S14)：基于D，计算X射线中每个窄能谱段的分离投影p_r，r＝1,2…R，R为X射线的窄能谱段总数；

步骤13得到的矩阵D＝(d_k,m)_K×M，其中M为每个扫描电压的投影图像序列的子单元总数，通过式(2)可得分离投影p_r：

步骤15(S15)：重建p_r得到每个窄能谱段的重建图像a_r。

例如通过对p_r进行一般线性重建(闫镔,李磊.CT图像重建算法[M].科学出版社,2014)得到每个窄能谱段的重建图像a_r。

以下给出步骤13的原理说明。

首先需将步骤12获得的N个投影图像序列转化为F＝(f_n,m)_N×M；M为每个扫描电压的投影图像序列的子单元总数，由于每个子单元对应一个射线路径，M也可以理解每个扫描电压下X射线的射线路径总数。

图2为某个扫描电压下获得的投影图像序列，将其中每幅图像都如图3所示以图像像素作为子单元，并为每个子单元建立索引，则M＝T×n1×m1。

将每个扫描电压下获得的投影图像序列都以图像像素作为子单元，并为每个子单元建立索引，就可以得到F矩阵如下。

F是一个N×M的矩阵。其中f_n,m为第n个扫描电压投影图像序列中第m条路径上的射线强度(对应图像中该单元的灰度值)与入射强度比值。入射强度，一般指未穿过物体区域的射线路径上的灰度值。

在X射线成像过程中，X射线的多谱特性使得探测器所获得的投影图像数据中包含着多谱段的衰减，这为利用多谱投影图像进行投影分解提供了依据。

设X射线路径L上的强度为I，初始强度为I₀，归一化光谱为S(E)，最大能量为E_max，物体在空间位置x处衰减系数为μ(x,E)，根据多谱衰减表达式有

令E_max＝max{E_1,max,…,E_n,max,…,E_N,max}，E_n,max表示第n个扫描电压的最大光子能量。依据实际情况将能量区间[0,E_max]划分为R个窄能谱段子区间[E'_r-1,E'_r]，r＝1,2…R，其中0＝E'₀<E'₁<…<E'_R＝E_max，通常等间隔取

区间内能量E_r＝t·E'_r+(1-t)·E'_r-1，一般取t＝0.5。

若成像物质由K种材质组成，假设能量区间[E'_r-1,E'_r]充分窄使得第k种材质在该区间上的衰减系数恒定为μ_k(E_r)，E_r取定后，各材质的线性衰减系数μ_k(E_r)参照美国NIST网站可得。

令在路径L上被测物的第k种材质的厚度为d_k，于是式(3)可转化为

其中

称为光谱分解系数，满足

它通常是未知的，而

即为待求的分解投影。

而d_k与射线路径有关，即不同射线路径上每种材质的厚度不同，要求得单一路径上多种材质的厚度仅靠单扫描电压下的多谱投影是严重欠定问题，所以需要采集多个扫描电压下不同能量的投影图像。同时由于光谱分解系数未知，故一般需要采集扫描电压总数N大于物质所包含材质的数目K。

每个扫描电压的投影图像序列中共包含M条射线路径，f_n,m为第n个扫描电压的投影图像序列中第m条路径上的射线强度与入射强度比值，s_n,r为第n个扫描产生X射线谱上第r个光谱分解系数，d_k,m为第m条路径上第k种材质的厚度。假设每条路径上的射线强度与入射强度比值是独立分布的泊松随机变量，即

令S＝(s_n,r)_N×R为光谱分解系数矩阵，R为X射线的窄能谱段总数；D＝(d_k,m)_K×M为投影厚度矩阵。

记

为理想的射线强度与入射强度比值。

则步骤13为根据已知的F＝(f_n,m)_N×M，求解满足泊松极大似然原理的变能量多谱衰减表达式

中的被测物的所有材质对应的投影厚度矩阵D。

以下给出步骤13的第一实施例。

对于N×M个泊松多谱投影测量值f_n,m，以广义KL散度形式的负对数似然函数来衡量实际数据f_n,m与理想值

之间的差异，其似然函数形式如下：

其中，S、D为自变量，F的每一行是以投影图像序列像素分布依次排列并包含不同投影角度的投影数据。于是得到如下基于泊松先验的变能量X射线投影图像第一分解模型：

其中：

S＝(s_n,r)_N×R为光谱分解系数矩阵，D＝(d_k,m)_K×M为投影厚度矩阵。

μ_k(E_r)为第k种材质在E_r时的线性衰减系数，E_r为第r个窄能谱段的光子能量。

即求解满足泊松极大似然原理的变能量多谱衰减表达式

中D和S，也可以转换为求解满足第一分解模型的D和S。

根据光谱分解系数和射线路径上材质的厚度均为非负值，得到求解变量S和D满足非负约束，记S≥0,D≥0表示矩阵内的每一个元素s_n,r≥0，d_k,m≥0，以S·1_R＝1_N表示每个电压下的光谱衰减系数和为1，1_R表示R维全1列向量。

如图4所示，在步骤13的第一实施例中，计算满足第一分解模型的被测物所有材质对应的投影厚度矩阵D包括：

步骤21：设定S和D的初始值：

E_n,max为第n个扫描电压所产生的最大光子能量。

步骤22：固定D，使用非负矩阵分解方法对满足第一分解模型的S进行更新；

步骤23：固定S，使用加速近端梯度算法对满足第一分解模型的D进行更新；

步骤24：计算D值更新前后的差值，判断差值是否满足第三预设条件；如果是，输出更新后的D，如果否，则返回步骤22。

例如，第三预设条件为

或

或||Dⁱ-D^i-1||_F＜ε；ε值根据实际情况选定，如ε＝10^-3。

其中，步骤22中使用非负矩阵分解方法对满足第一分解模型的S进行更新包括：

其中：

i表示迭代次数，初始值为1。

其中，步骤23中使用加速近端梯度算法对满足第一分解模型的D进行更新包括：

步骤41：Cⁱ＝0；

步骤42：更新

其中：

tⁱ为搜索步长，

第四预设条件为第(i-η)～第(i-1)次的搜索步长不变；

步骤43：判断

是否成立，如果是，执行步骤24，如果否，则令Cⁱ＝Cⁱ+1，返回步骤42。

以下给出步骤13的第二实施例。

记f_n,m和f_n,m+1表示相邻射线路径上的射线强度与入射强度比值，由于物体组分呈连续性分布，故式(7)中相邻射线路径上每种材质的投影厚度值d_k,m和d_k,m+1差异性较小，所以可以对投影厚度矩阵添加光滑先验，建立如下惩罚函数

式中D_k＝(d_k,m)_M×1，Q表示二维或一维离散梯度算子，

为Tikhonov(L2)正则化函数，β_k为正则化参数，使用该参数可以为每种材质添加光滑性程度不同的先验。

在式(7)中增加式(9)可得第二分解模型如下，式中α为正则化参数。

基于此，步骤13还可以包括：

将N个投影图像序列转化为F＝(f_n,m)_N×M输入第二分解模型，M为每个扫描电压的投影图像序列的子单元总数，计算满足第二分解模型的被测物所有材质对应的投影厚度矩阵D；

第二分解模型为

其中：

S＝(s_n,r)_N×R为光谱分解系数矩阵；D＝(d_k,m)_K×M为投影厚度矩阵；

μ_k(E_r)为第k种材质在E_r时的线性衰减系数，E_r为第r个窄能谱段的光子能量；

α为正则化参数，可设为10^-7≤α≤10^-3；

D_k＝(d_k,m)_M×1，Q为二维或一维离散梯度算子；

为Tikhonov(L2)正则化函数，β_k为正则化参数。

如图5所示，在步骤13的第二实施例中，计算满足第二分解模型的被测物所有材质对应的投影厚度矩阵D包括：

步骤31：设定S和D的初始值：

步骤32：固定D，使用非负矩阵分解方法对满足第二分解模型的S进行更新；

步骤33：固定S，使用加速近端梯度算法对满足第二分解模型的D进行更新；

步骤34：计算D值更新前后的差值，判断差值是否满足第三预设条件；如果是，输出更新后的D，如果否，则返回步骤32。

例如，第三预设条件为

或

或||Dⁱ-D^i-1||_F＜ε；ε值根据实际情况选定，如ε＝10^-3。

其中，步骤32中使用非负矩阵分解方法对满足第二分解模型的S进行更新包括：

其中：

i表示迭代次数，初始值为1。

其中，步骤33中使用加速近端梯度算法对满足第一分解模型的D进行更新包括：

步骤51：Cⁱ＝0；

步骤52：更新

其中：

tⁱ为搜索步长，

第四预设条件为第(i-η)～第(i-1)次的搜索步长不变；

为梯度算子Q的转置矩阵的第m行向量。

步骤53：判断

是否成立，如果是，执行步骤34，如果否，则令Cⁱ＝Cⁱ+1，返回步骤52。

此外，图1中步骤13的求解也可以参照非负矩阵模型基本求解方法的推导，该方法对初始值的依赖性较强，而且在极小值附近收敛缓慢。

以上是对本发明X射线多谱分离成像方法的说明。

本发明的X射线多谱分离成像方法可以实现X射线投影图像盲分离，根据物体的材质、尺寸，选取合理的成像参数获取X射线投影图像序列；然后依据投影图像序列，求解满足泊松极大似然原理的变能量多谱衰减表达式中的被测物的所有材质对应的投影厚度矩阵D；再根据矩阵D实现多谱投影图像的投影分离。相比于基于硬件和其他多能谱重建的X射线CT成像技术，该方法从投影域实现了多谱X射线图像的能谱分离，无需对现有的高动态CT成像设备进行改进，具有成本低、能谱分辨率高等优点，对于推广定量化CT成像技术的应用有重要意义。

本发明还提供一种非瞬时计算机可读存储介质，非瞬时计算机可读存储介质存储指令，指令在由处理器执行时使得处理器执行上述的X射线多谱分离成像方法中的步骤。

本发明还提供一种X射线多谱分离成像装置，包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，处理器执行计算机程序时实现上述的X射线多谱分离成像方法中的步骤。

如图6所示，该装置包括：

扫描电压范围设定模块：设定被测物的成像电压范围[V1,V2]，当V∈[V1,V2]时，被测物的穿透率符合第一预设条件、且被测物的任意两种材质在V时的线性衰减系数的相对偏差满足预设要求对应的能量范围符合第二预设条件；

变能量全扫描模块：在[V1,V2]的范围内取N个不同的扫描电压，N大于被测物所包含的材质总数K，且任意两个扫描电压之间的差值大于第一预设值，在每个扫描电压下对被测物进行全角度扫描成像，获得符合成像灰度要求的N个投影图像序列，全角度扫描包括T个不同的成像角度；

D解算模块：基于N个投影图像序列，求解满足泊松极大似然原理的变能量多谱衰减表达式中的被测物的所有材质对应的投影厚度矩阵D；

分离投影计算模块：基于D，计算X射线中每个窄能谱段的分离投影p_r，r＝1,2…R，R为X射线的窄能谱段总数；

分离投影图像重建模块：重建p_r得到每个窄能谱段的重建图像a_r。

其中，D解算模块还可以包括：

将N个投影图像序列转化为F＝(f_n,m)_N×M输入第一分解模型，M为每个扫描电压的投影图像序列的子单元总数，计算满足第一分解模型的被测物所有材质对应的投影厚度矩阵D；

第一分解模型为

其中：

S＝(s_n,r)_N×R为光谱分解系数矩阵；

D＝(d_k,m)_K×M为投影厚度矩阵；

μ_k(E_r)为第k种材质在E_r时的线性衰减系数，E_r为第r个窄能谱段的光子能量。

进而，计算满足第一分解模型的被测物所有材质对应的投影厚度矩阵D包括：

初值设定模块1：设定S和D的初始值：

S更新模块1：固定D，使用非负矩阵分解方法对满足第一分解模型的S进行更新；

D更新模块1：固定S，使用加速近端梯度算法对满足第一分解模型的D进行更新；

判断模块1：计算D值更新前后的差值，判断差值是否满足第三预设条件；如果是，输出更新后的D，如果否，则返回S更新模块1。

S更新模块1中：使用非负矩阵分解方法对满足第一分解模型的S进行更新包括：

其中：

i表示迭代次数，初始值为1。

D更新模块1中：使用加速近端梯度算法对满足第一分解模型的D进行更新包括：

初值设定模块2：Cⁱ＝0；

D更新模块2：更新

其中：

i表示迭代次数，初始值为1；

tⁱ为搜索步长，

第四预设条件为第(i-η)～第(i-1)次的搜索步长不变；

判断模块2：判断

是否成立，如果是，执行判断模块1，如果否，则令Cⁱ＝Cⁱ+1，返回D更新模块2。

其中，D解算模块还可以包括：

将N个投影图像序列转化为F＝(f_n,m)_N×M输入第二分解模型，M为每个扫描电压的投影图像序列的子单元总数，计算满足第二分解模型的被测物所有材质对应的投影厚度矩阵D；

第二分解模型为

其中：

S＝(s_n,r)_N×R为光谱分解系数矩阵；

D＝(d_k,m)_K×M为投影厚度矩阵；

μ_k(E_r)为第k种材质在E_r时的线性衰减系数，E_r为第r个窄能谱段的光子能量；

α为正则化参数；

D_k＝(d_k,m)_M×1，Q为二维或一维离散梯度算子；

为Tikhonov(L2)正则化函数，β_k为正则化参数。

进而，计算满足第二分解模型的被测物所有材质对应的投影厚度矩阵D包括：

初值设定模块3：设定S和D的初始值：

S更新模块3：固定D，使用非负矩阵分解方法对满足第二分解模型的S进行更新；

D更新模块3：固定S，使用加速近端梯度算法对满足第二分解模型的D进行更新；

判断模块3：计算D值更新前后的差值，判断差值是否满足第三预设条件；如果是，输出更新后的D，如果否，则返回S更新模块3。

S更新模块3中：使用非负矩阵分解方法对满足第二分解模型的S进行更新包括：

其中

i表示迭代次数，初始值为1。

D更新模块3中：使用加速近端梯度算法对满足第二分解模型的D进行更新包括：

初值设定模块4：Cⁱ＝0；

D更新模块4：更新

其中：

i表示迭代次数，初始值为1；

tⁱ为搜索步长，

第四预设条件为第(i-η)～第(i-1)次的搜索步长不变；

为梯度算子Q的转置矩阵的第m行向量。

判断模块4：判断

是否成立，如果是，执行判断模块3，如果否，则令Cⁱ＝Cⁱ+1，返回D更新模块4。

可选地，T大于等于180。

可选地，T个不同的成像角度包括：初始成像角度为θ₀、角度间隔为△θ的T个成像角度。

需要说明的是，本发明的X射线多谱分离成像装置的实施例，与X射线多谱分离成像方法的实施例原理相同，相关之处可以互相参照。

以上所述仅为本发明的较佳实施例而已，并不用以限定本发明的包含范围，凡在本发明技术方案的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (13)

1.一种X射线多谱分离成像方法，其特征在于，所述方法包括：

步骤11：设定被测物的成像电压范围[V1,V2]，当V∈[V1,V2]时，所述被测物的穿透率符合第一预设条件、且所述被测物的任意两种材质在所述V时的线性衰减系数的相对偏差满足预设要求对应的能量范围符合第二预设条件；其中，所述第一预设条件：依据被测物的材质、尺寸、探测器成像动态范围，并避免欠曝光和过曝光现象造成信息的缺失，选取的电压区间[V₁,V₂]使物体的穿透率达到30％～70％；所述第二预设条件为：从被测物组分分析的角度出发，依据材质的线性衰减系数曲线，满足任意电压V∈[V₁,V₂]时任意第k₁、k₂种材质满足能量范围Ω与V所产生的能量范围Ω₁的交集占Ω₁的60％以上，以保证各组分在选取能量范围上具有可区分性；

步骤12：在[V1,V2]的范围内取N个不同的扫描电压，所述N大于所述被测物所包含的材质总数K，且任意两个扫描电压之间的差值大于第一预设值，在每个所述扫描电压下对所述被测物进行全角度扫描成像，获得符合成像灰度要求的N个投影图像序列，所述全角度扫描包括T个不同的成像角度；

步骤13：基于所述N个投影图像序列，求解满足泊松极大似然原理的变能量多谱衰减表达式中的所述被测物的所有材质对应的投影厚度矩阵D；

步骤14：基于所述D，计算所述X射线中每个窄能谱段的分离投影p_r，r＝1,2…R，R为所述X射线的窄能谱段总数；

步骤15：重建所述p_r得到所述每个窄能谱段的重建图像a_r。

2.根据权利要求1所述的方法，其特征在于，所述步骤13包括：

将所述N个投影图像序列转化为F＝(f_n，m)_N×M输入第一分解模型，M为每个扫描电压的投影图像序列的子单元总数，计算满足所述第一分解模型的所述被测物所有材质对应的投影厚度矩阵D；

所述第一分解模型为

其中：

S＝(s_n，r)_N×R为光谱分解系数矩阵；

D＝(d_k，m)_N×M为投影厚度矩阵；

μ_k(E_r)为第k种材质在E_r时的线性衰减系数，E_r为第r个窄能谱段的光子能量。

3.根据权利要求2所述的方法，其特征在于，所述计算满足所述第一分解模型的所述被测物所有材质对应的投影厚度矩阵D包括：

步骤21：设定所述S和D的初始值：

步骤22：固定所述D，使用非负矩阵分解方法对满足所述第一分解模型的S进行更新；

步骤23：固定所述S，使用加速近端梯度算法对满足所述第一分解模型的D进行更新；

步骤24：计算D值更新前后的差值，判断所述差值是否满足第三预设条件；如果是，输出更新后的D，如果否，则返回步骤22。

4.根据权利要求3所述的方法，其特征在于，所述使用非负矩阵分解方法对满足所述第一分解模型的S进行更新，包括：

其中：

i表示迭代次数，初始值为1。

5.根据权利要求1所述的方法，其特征在于，所述步骤13包括：

将所述N个投影图像序列转化为F＝(f_n，m)_N×M输入第二分解模型，M为每个扫描电压的投影图像序列的子单元总数，计算满足所述第二分解模型的所述被测物所有材质对应的投影厚度矩阵D；

所述第二分解模型为

其中：

S＝(s_n，r)_N×R为光谱分解系数矩阵；

D＝(d_k，m)_K×M为投影厚度矩阵；

μ_k(E_r)为第k种材质在E_r时的线性衰减系数，E_r为第r个窄能谱段的光子能量；

α为正则化参数；

D_k＝(d_k，m)_M×1，Q为二维或一维离散梯度算子；

为Tikhonov(L2)正则化函数，β_k为正则化参数。

6.根据权利要求5所述的方法，其特征在于，所述计算满足所述第二分解模型的所述被测物所有材质对应的投影厚度矩阵D包括：

步骤31：设定所述S和D的初始值：

步骤32：固定所述D，使用非负矩阵分解方法对满足所述第二分解模型的S进行更新；

步骤33：固定所述S，使用加速近端梯度算法对满足所述第二分解模型的D进行更新；

步骤34：计算D值更新前后的差值，判断所述差值是否满足第三预设条件；如果是，输出更新后的D，如果否，则返回步骤32。

7.根据权利要求6所述的方法，其特征在于，所述使用非负矩阵分解方法对满足所述第二分解模型的S进行更新，包括：

其中：

i表示迭代次数，初始值为1。

8.根据权利要求3所述的方法，其特征在于，所述使用加速近端梯度算法对满足所述第一分解模型的D进行更新包括：

步骤41：Cⁱ＝0；

步骤42：更新

其中：

i表示迭代次数，初始值为1；

tⁱ为搜索步长，

ξ∈(0，1)，

所述第四预设条件为第(i-η)～第(i-1)次的搜索步长不变；

步骤43：判断

是否成立，如果是，执行步骤24，如果否，则令Cⁱ＝Cⁱ+1，返回步骤42。

9.根据权利要求6所述的方法，其特征在于，所述使用加速近端梯度算法对满足所述第二分解模型的D进行更新包括：

步骤51：Cⁱ＝0；

步骤52：更新

其中：

i表示迭代次数，初始值为1；

tⁱ为搜索步长，

ξ∈(0，1)，

所述第四预设条件为第(i-η)～第(i-1)次的搜索步长不变；

为梯度算子Q的转置矩阵的第m行向量；

步骤53：判断

是否成立，如果是，执行步骤34，如果否，则令Cⁱ＝Cⁱ+1，返回步骤52。

10.根据权利要求1所述的方法，其特征在于，所述T大于等于180。

11.根据权利要求10所述的方法，其特征在于，所述T个不同的成像角度包括：初始成像角度为θ₀、角度间隔为△θ的T个成像角度。

12.一种非瞬时计算机可读存储介质，所述非瞬时计算机可读存储介质存储指令，其特征在于，所述指令在由处理器执行时使得所述处理器执行如权利要求1至11中任一所述的X射线多谱分离成像方法中的步骤。

13.一种X射线多谱分离成像装置，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至11中任一所述的X射线多谱分离成像方法中的步骤。

Images