基于基效应分解的多谱投影序列盲分离算法研究

引用本文

赵耀霞, 韩焱, 陈平. 基于基效应分解的多谱投影序列盲分离算法研究[J]. 光谱学与光谱分析, 2019,39(9): 2763-2768.
ZHAO Yao-xia, HAN Yan, CHEN Ping. Blind Separation of Multi-Voltage Projection Sequence Based on Fundamental Effect Decomposition[J]. Spectroscopy and Spectral Analysis, 2019,39(9): 2763-2768.
Doi:10.3964/j.issn.1000-0593(2019)09-2763-06 复制到剪切板

Permissions

《光谱学与光谱分析》期刊社所有

基于基效应分解的多谱投影序列盲分离算法研究

赵耀霞, 韩焱, 陈平^*

中北大学信息与通信工程学院, 信息探测与处理山西省重点实验室, 山西太原 030051

*通讯联系人 e-mail: pc0912@163.com

作者简介: 赵耀霞, 1979年生, 中北大学信息与通信工程学院博士研究生 e-mail: zhao@nuc.edu.cn

收稿日期: 2018-08-02 接受日期: 2018-12-10

基金: 国家自然科学基金项目(61571404, 61801437, 61601412), 山西省高等学校优秀青年学术带头人课题资助

摘要

在常规CT成像系统中, 发出的X射线是连续能谱的, 导致重建图像出现硬化伪影, 影响了材料组分区分, 无法进行定量表征。解决这一问题的关键在于实现多能谱CT成像, 利用多个窄能谱段或单能量CT图像, 提高组分与图像灰度的对应性。相对于传统CT, 能谱CT具有更强的组分区分能力, 有利于实现物体组分的定量分析。现有的基于光子计数探测器多谱CT在成像时间分辨率和空间分辨率存在局限, 基于能谱滤波的多谱CT能谱区分度受限。而基于变电压多谱投影序列盲分离的多谱CT, 通过分解连续能谱投影, 获取窄能谱投影, 进而实现能谱CT成像, 确保物质组分与重建图像灰度值的对应性, 实用性较强。但是由于X射线能谱和物体组分的未知, 在盲分离过程中, 衰减系数未知, 并且能谱划分是不确定的, 导致窄能谱CT重建图像的能量指向性不强, 对应能量值与参考能量偏差偏大, 影响组分定量分析。因此, 针对盲分离中能谱划分不确定性和重建图像能量指向性问题的开展研究。利用衰减系数的光电效应和康普顿效应分解, 构建能量约束, 消除能量的不确定性, 降低分解所得投影重建图像的能量与参考能量值的偏差。在基于以残差的局部方差和最小为优化目标的分解模型中, 将分解模型中的衰减系数按光电效应和康普顿效应分解为能量项和材料项, 利用能量项的可预知性, 依据预先划分的窄能谱段设置其值, 固定各分解投影对应的窄能谱段, 作为对能量的约束条件。求解所得各分解投影为能量已知投影, 对其重建可得到能量确定的各窄能谱段的图像。选择衰减系数相近的硅铝材质构成外硅内铝圆柱体进行实验验证, 在有能量约束的求解结果中, 硅铝衰减系数与参考值偏差小于无能量约束, 所得重建图像中硅铝变化率与理论值趋势较一致, 能量指向性强, 与参考能量偏差降低。结果表明, 所提方法解决了基于变电压序列盲分离多谱CT成像的能量指向问题, 能谱分辨率更高, 组分表征更准确。

关键词: 多谱CT; 多电压投影序列; 盲分离; 窄能谱投影; 能量指向

中图分类号:O434.1 文献标志码:A

Blind Separation of Multi-Voltage Projection Sequence Based on Fundamental Effect Decomposition

ZHAO Yao-xia, HAN Yan, CHEN Ping^*

School of Information and Communication Engineering, Shanxi Key Laboratory of Signal Capturing & Processing, North University of China, Taiyuan 030051, China

*Corresponding author

Abstract

In X-ray CT imaging system, the polychromatic X-ray results in beam hardening artifacts, which influence the material composition distinction. It can’t realize quantitative characterization. The multi-spectrum CT can reach the higher correspondence between the composition and image gray by narrow-energy-width or monoenergetic CT. Compared with traditional CT, spectral CT can distinguish the different component. The implement method based on the photon counting detector is limited in temporal resolution and spatial resolution. The implement method based on filter doesn’t have enough discrimination. And the implement method based on blind separation of multi-voltage projection sequences has better practicability. It guarantees the correspondence between the materials and the gray of reconstructed image. However, the attenuation coefficients are unknown and energy spectrum division is uncertain, the analytic energy value of reconstructed images is ambiguous. The error between the analytic energy and referenced energy is higher. Then it will influence the precision of the multi-component quantitative analysis. For this problem, an improved method is proposed. It utilizes the attenuation coefficient decomposition of compton effect and photoelectric effect as an energy constraint to eliminate the uncertainty energy partition. The error between the energy value of reconstructed images with decomposed projection and the reference energy value is reduced. In the decomposition model, the optimum object function is local variance sum of residual error minimum. The attenuation coefficient is decomposed as energy dependency term and material dependency term according to compton effect and photoelectric effect. The energy dependency term can be known in advance. It can be as energy constraint and used to fix the energy value of narrow-energy-width. Then the energy of decomposed projection is known and the energy of corresponding reconstructed images also is known. A cylinder composed of aluminum and silicon is used in the verification experiment since aluminium and silicon have approximate attenuation coefficients. The error between the attenuation coefficients of reconstructed images with the energy constraint is less than the result of reconstructed images without the energy constraint. The contrast tendency of silicon and aluminium with energy is close to the theoretical value. Also the difference with reference energy is reduced. The result shows that the proposed method solves the energy directivity problem of multi-spectrum CT based on the blind separation of multi-voltage projection sequences. The energy spectrum resolution ratio is higher. The composition representation is more accurate.

Keyword: Multi-spectrum CT; Multi-voltage projection sequences; Blind separation; Narrow-spectrum projection; Energy order

文章图片

引言

多谱CT相对于常规CT, 考虑了X射线的连续性, 构建了基于多谱X射线的投影观测模型, 通过重建, 抑制了硬化伪影, 既能实现缺陷检测的定性分析, 也能实现组分区分的定量分析^{[1, 2]}。双能CT是一种典型多谱CT, 是在获取高低能投影基础上, 采用基效应或基材料分解投影, 重建分解系数, 得到重建图像^[3]。但仍受射线多谱特性影响, 一般为“ 伪双能CT” , 能谱区分度不够, 组分区分度有限^{[4, 5]}。光子计数型探测器的出现, 可以在接收端实现能谱分离, 得到窄能谱段投影, 实现窄能谱CT重建以及组分定量分析^[6]。但是, 由于计数率不足、成像效率低、分辨率低等问题, 实际应用受限^[7]。借鉴光子计数型探测器的分能谱段成像思想, 提出基于能谱滤波分离的多谱CT成像方法, 通过在射线出束口添加不同厚度、材质滤波片方法减小能谱宽度, 实现窄能谱CT成像^[8]。另一种是基于盲分离思想, 建立优化模型, 对多电压X射线图像分解, 获取窄能谱投影。目前主要有以加权误差平方和最小和以偏差的局部方差和最小的优化模型, 重建图像硬化伪影能够得到抑制, 符合窄能谱图像特性, 但是两模型中都缺乏能量约束使得所得重建图像序列顺序混乱, 对应能量不明确, 与参考值偏差较大, 不利于后续的定量分析应用^{[9, 10]}。

针对多电压X射线图像分解所得投影的重建图像能量不明确, 对应能量值与参考能量偏差偏大问题, 拟改进局部方差和的分离模型, 并借助衰减系数的光电效应和康普顿效应分解, 修正重建图像序列的能量顺序, 降低重建图像对应能量值的偏差。

1 多电压投影序列分解模型

1.1 多电压X射线成像模型

实际X射线成像原理的离散形式可以表述为

$I = I_{0} \overset{R}{\sum_{r = 1}} S (E_{r}) \exp (- \overset{K}{\sum_{k = 1}} u_{k} (E_{r}) d_{k})$ (1)

其中, r为第r个离散的窄能谱段, k为构成物体的第k种材质, d_k为射线穿过路径上第k种材质的长度, u_k(E_r)为第k种材质在第r个能量段上的衰减系数, S(E_r)为第r个窄能谱段对应的加权系数, 其与能谱等因素有关, 满足

$\overset{R}{\sum_{r = 1}} S (E_{r}) = 1$ (2)

对式(1)两边同时除以I₀, 可得

$\frac{I}{I_{0}} = \overset{R}{\sum_{r = 1}} S (E_{r}) \exp (- \overset{K}{\sum_{k = 1}} u_{k} (E_{r}) d_{k})$ (3)

对多个电压下的X射线成像, 将第n个电压下第m个像素对应的I/I₀记为f_nm, 则多个电压X射线成像可写成矩阵形式

$F = S \exp (- UD) + σ$ (4)

其中 $F = (f_{nm})_{NM}, S = (s_{nr})_{NR}, U = (u_{rk})_{RK}, D = (d_{km})_{KM}, σ$ 表示散射等造成的偏差。矩阵S的行是构成各电压X射线图像的窄能谱段的加权系数。 U的一列对应一种材质在各窄谱段的衰减系数, D的相应行是同种材料在各像素上对应的衰减长度。

造成X射线衰减的因素有瑞利散射、光电效应、康普顿散射和电子对效应, 在十几 keV到1 MeV之间时, 衰减主要由光电效应和康普顿效应造成, 衰减系数可以按照这两种基效应分解为两者的线性组合

$u_{k} (E_{r}) = a_{c} (k) f_{KN} (E_{r}) + a_{p} (k) f_{p} (E_{r})$ (5)

式中, “ c” 和“ p” 分别表示康普顿效应和光电效应, f_KN(E)和f_p(E)只与能量有关, 其值可以预先确定, a_c(k)和a_p(k)则只与材料有关。令

$U' = [\begin{array}{l} f_{KN} (E_{1}) & f_{p} (E_{1}) \\ ︙ & ︙ \\ f_{KN} (E_{R}) & f_{p} (E_{R}) \end{array}]$ (6)

$D' = AD = [\begin{array}{l} a_{c 1} & \dots & a_{cK} \\ a_{p 1} & \dots & a_{pK} \end{array}] D$ (7)

U'是R× 2矩阵, D'是2× M矩阵, 则模型(5)可以改写为

$F = S \exp (- U'D') + σ$ (8)

若先不考虑σ , 从解方程角度考虑模型(8), 要求方程个数不少于未知量个数, 即

$NM \geq NR + 2 M$ (9)

1.2 投影分解模型

显然, 由于σ 未知, 式(8)只能通过优化模型求解。在X射线成像中, 散射表现为低频信号。低频, 说明变化缓慢, 信号相邻采样点的差异小, 可以采用方差作为对散射图像低频特性的描述。考虑到散射图像与物体相关, 对X射线图像每一点, 取其一个局部图像, 以其局部图像的方差作为其变化快慢程度的衡量, 对所有点的局部方差求和作为图像整体的低频描述。变电压X射线图像的分解模型改为如式(10)形式

$\begin{array}{l} \min \overset{N}{\sum_{n = 1}} \overset{M - T}{\sum_{m = T + 1}} [Var {[F - S \exp {(- U'D')]}_{nm 局部}}] \\ s.t. S \geq 0, U' \geq 0, D' \geq 0 \end{array}$ (10)

其中, T是与局部图像有关的参数, 需要在求解前设置。以二维重建为例, X射线图像的局部图像大小取为1× (2T+1), 当前点在局部图像的中心。利用KKT条件, 可推导出模型的迭代求解公式^[11]。此外, 由式(2)可知, 在迭代求解中需要对S的行进行归一化。

由于U'中f_KN(E)和f_p(E)的值只由能量决定, 因此, 只要确定了窄能谱段划分方式, 将能谱段的中间能量代入f_KN(E)和f_p(E), 计算可得U'值, 分解所得的每个投影都有明确的能量对应, 即通过f_KN(E)和f_p(E)约束了分解所得投影对应的能量, 避免了能量的混乱。

2 结果与讨论

为验证图像窄能谱特性, 选择衰减系数相近的硅铝材质构成的内铝外硅圆柱体为实验对象, 横截面直径为40 mm, 内部铝圆柱直径30 mm。根据NIST数据库, 在10~140 keV范围内, 硅与铝的衰减系数随能量变化, 逐渐变小, 约60 keV衰减系数反转变化[如图1(a)所示], 对应算出硅铝衰减系数相对差异如图1(b)所示。

实验系统是由450kV的GE射线机(ISOVOLT-450)和平板探测器(PerkinElmer XRD1621)组成的CT系统。探测器与射线源距离为140 cm, 物体中心距探测器20 cm。选择成像电压为120, 130和140 kV, 电流为3 mA。

实验中各电压下投影及重建图像如图2所示。

	Figure Option View Download New Window
	图2 各电压下投影及重建图像第一行为120 kV, 第二行为130 kV, 第三行为140 kVFig.2 Projection and reconstruction of images at each voltage First row 120 kV, second row 130 kV, third row 140 kV

从图2中, 可以看出由于X射线的多谱效应, 图像硬化伪影较严重, 影响组分分析。根据本研究提出的方法, 针对图2的投影, 以10 keV为能量间隔划分能谱, 选择每个能量段中间能量值作为参考值。取经验参数T=1, 当两次迭代的目标函数值相对差异小于0.01时, 停止迭代, 最大迭代次数设置为500。

以偏差的局部方差和最小为优化目标的无能量约束模型中, 分解120, 130和140 kV下投影, 所得到部分重建图像如图3所示。

	Figure Option View Download New Window
	图3 无能量约束模型方法从上到下依次为第2, 7, 12投影及CT图像Fig.3 Model without energy constraint From the top down: the second, seventh and 12th projections and CT images

从分解重建的图像(图3)可以看出, 通过投影序列盲分离可以有效实现近似窄谱成像, 硬化伪影得到有效抑制。利用所提出的分解方法分解结果及重建后图像如图4。

	Figure Option View Download New Window
	图4 有能量约束方法从上到下依次为第2, 7, 12投影及CT图像Fig.4 Model with energy constraint From the top down: the second, seventh and 12th projections and CT images

从分解的重建结果来看, 能够达到与无能量约束模型一致的结果。进一步针对能谱指向性进行分析。无能量约束模型解得的衰减系数与参考值的差异图如图5所示, 有能量约束模型解得的衰减系数与参考值的差异图如图6所示。

	Figure Option View Download New Window
	图5 无能量约束方法的衰减系数 (a): 铝; (b): 硅Fig.5 Attenuation coefficient without energy constraints method (a): Aluminum; (b): Silicon

	Figure Option View Download New Window
	图6 有能量约束方法的衰减系数 (a): 铝; (b): 硅Fig.6 Attenuation coefficient with energy constraints method (a): Aluminum; (b): Silicon

从图5(a, b)和图6(a, b)看, 通过本文方法, 由重建图像计算出的铝衰减系数与参考衰减系数的相对偏差较小。进一步分析由重建序列得到的硅铝对比度如图7(a, b)所示。

	Figure Option View Download New Window
	图7 重建图像硅铝对比度 (a): 无能量约束; (b): 有能量约束Fig.7 The Contrast of silicon and aluminum in the reconstructed images (a): Without energy constraint; (b): With energy constraint

从图7(a)可看到, 无能量约束的解析方法对比度趋势杂乱, 与图1(b)的趋势不符, 而本文提出的分解方法, 其对比度变化趋势[图7(b)]与图1(b)的趋势较吻合, 差异更小。

	Figure Option View Download New Window
	图1 硅铝衰减系数及其相对差异 (a): 硅铝衰减系数; (b): 硅铝衰减系数相对差异Fig.1 Attenuation coefficient and their relative difference of aluminium and silicon (a): Attenuation coefficient of aluminium and silicon; (b): Attenuation coefficient relative difference of aluminium and silicon

因此, 从上述实验结果来看, 相对于无能量约束模型, 本研究提出的分解方法所得到的窄谱图像经CT重建, 既能有效抑制多谱硬化伪影, 同时也保持了图像序列的能量顺序性, 硅铝衰减系数对应能量值与参考值偏差更小。

3 结论

针对利用多电压X射线投影序列盲分离, 获取窄能谱投影重建图像能量指向不明确, 对应能量值偏差偏大问题, 提出了一种以偏差的局部方差和最小为优化目标的基于基效应分析的投影序列盲分离模型, 利用衰减系数的光电效应和康普顿效应分解作为约束条件, 修正分解能量顺序, 使得所求解窄能谱重建图像对比度较高, 硬化伪影小, 而且有了明确的能量指向, 降低了重建图像对应能量值与参考值的偏差。并以硅铝圆柱体为模体的验证实验, 验证了方法的有效性。

参考文献

文献列表

[1]	Fang L, Yin X, Wu L, et al. International Journal of Pharmaceutics, 2017, 531(2): 658. [本文引用:1]
[2]	Ketcham R A, Hanna R D. Computers & Geosciences, 2014, 67(4): 49. [本文引用:1]
[3]	HAO Jia, ZHANG Li, XING Yu-xiang, et al(郝佳, 张丽, 邢宇翔, 等). Journal of Tsinghua University·Science and Technology(清华大学学报·自然科学版), 2011, 51(4): 457. [本文引用:1]
[4]	LI Bao-lei, ZHANG Ping-yu, LI Bin, et al(李保磊, 张萍宇, 李斌, 等). Acta Optica Sinica(光学学报), 2017, 37(10): 1034001. [本文引用:1]
[5]	Goodsitt M M, Christodoulou E G, Larson S C. Medical Physics, 2011, 38(4): 2222. [本文引用:1]
[6]	CHEN Zhi-qiang, ZHANG Li, JIN Xin(陈志强, 张丽, 金鑫). Chines Science Bulletin(科学通报), 2017, 62(13): 1350. [本文引用:1]
[7]	Wang S, Gao H, Zhang L, et al. Physics in Medicine & Biology, 2017, 62(6): 2208. [本文引用:1]
[8]	NIU Su-su, PAN Jin-xiao, CHEN Ping(牛素鋆, 潘晋孝, 陈平). Acta Optica Sinica(光学学报), 2014, 34(10): 1. [本文引用:1]
[9]	Wei J, Han Y, Chen P. Measurement Science & Technology, 2016, 27(2): 025402. [本文引用:1]
[10]	Wei J, Han Y, Chen P. International Journal of Biomedical Imaging, 2017, 2017, (2): 1. [本文引用:1]
[11]	Gillis N, Glineur F. CORE Discussion Papers/Center for Operations Research and Econometrics, 2008, 1: 1. [本文引用:1]

2017

0.0

... 引言多谱CT相对于常规CT, 考虑了X射线的连续性, 构建了基于多谱X射线的投影观测模型, 通过重建, 抑制了硬化伪影, 既能实现缺陷检测的定性分析, 也能实现组分区分的定量分析^[1,2] ...

2014

0.0

2011

0.0

... 双能CT是一种典型多谱CT, 是在获取高低能投影基础上, 采用基效应或基材料分解投影, 重建分解系数, 得到重建图像^[3] ...

2017

0.0

... , 能谱区分度不够, 组分区分度有限^[4,5] ...

2011

0.0

... , 能谱区分度不够, 组分区分度有限^[4,5] ...

2017

0.0

... 光子计数型探测器的出现, 可以在接收端实现能谱分离, 得到窄能谱段投影, 实现窄能谱CT重建以及组分定量分析^[6] ...

2017

0.0

... 但是, 由于计数率不足、成像效率低、分辨率低等问题, 实际应用受限^[7] ...

2014

0.0

... 借鉴光子计数型探测器的分能谱段成像思想, 提出基于能谱滤波分离的多谱CT成像方法, 通过在射线出束口添加不同厚度、材质滤波片方法减小能谱宽度, 实现窄能谱CT成像^[8] ...

2016

0.0

... 目前主要有以加权误差平方和最小和以偏差的局部方差和最小的优化模型, 重建图像硬化伪影能够得到抑制, 符合窄能谱图像特性, 但是两模型中都缺乏能量约束使得所得重建图像序列顺序混乱, 对应能量不明确, 与参考值偏差较大, 不利于后续的定量分析应用^[9,10] ...

2017

0.0

2008

0.0

... 利用KKT条件, 可推导出模型的迭代求解公式^[11] ...