多电压X射线图像分解的高对比度CT成像

魏交统，陈 平，韩 焱



多电压X射线图像分解的高对比度CT成像

魏交统，陈 平，韩 焱

( 中北大学 信息探测与处理山西省重点实验室，太原 030051 )

针对常规CT系统中X射线是连续宽能谱而导致重建图像中衰减特性相近组分对比度较低的问题，论文提出一种分解多电压X射线图像序列获取新投影的方法，用新投影重建的图像有效地提高了对衰减特性相近组分的区分能力。该方法依据X射线成像规律，将各电压下图像看作多个窄能谱图像的叠加，借鉴非负矩阵分解，以分解的误差平方和最小为目标对多电压X射线图像分解。用硅铝材质的圆柱体实验，新投影重建的图像相比直接重建图像，硅与铝对比度明显提高，表明了方法的有效性。

信号处理；CT成像；X射线图像；多电压；非负矩阵分解；对比度

0 引 言

常规CT成像系统的X射线是连续宽能谱，导致CT图像出现硬化伪影，对比度降低。重建图像灰度不仅与材质、能量有关，还与像素位置有关，使衰减特性相近的组分无法有效区分[1]。解决这类问题，理想情况是通过单能CT实现组分衰减系数与能量的唯一对应，如同步辐射CT[2]，但目前实际工程应用的CT系统尚难以实现单能X射线源[3]，一般通过双能或多能谱CT提高图像不同组分的对比度[4]。在医疗领域，双能CT能有效提高软组织与其他物质的对比度[5]，如2009年GE公司推出的宝石能谱CT。为进一步提高图像质量，人们引入了更多的能量信息，提出了多谱成像技术。近年来出现了具有单光子能量分辨能力的光子计数型探测器。光子计数型探测器能将接收到的光子按能量分为不相交的能带区间，进行窄能谱CT成像，实现组分的定量分析[6]。然而该探测器时间分辨力受限，工程应用存在局限[7]。借鉴这种多窄能谱成像思想，牛素鋆等提出了基于能谱滤波分离的多谱CT成像方法，添加不同滤波片获取多个窄能谱CT图像，但实际中添加滤波片使能谱变窄程度有限，不同能谱区分度不够，无法实现组分的有效分离[8]。

借鉴光子计数型探测器分能谱段成像的思想，论文对多个电压下X射线图像分解，以期获取窄能谱图像。该方法不需要硬件条件的改变，将单电压下X射线图像看作多个窄能谱图像的加权和，不同电压下图像仅是对应的窄能谱图像加权系数的不同。以各电压下X射线图像分解的总误差平方和最小为目标函数建立分解模型。求解模型获取新投影，重建图像，提高衰减系数相近组分的对比度。

1 变电压X射线成像模型

X射线管发出的是连续能谱X射线，成像模型可写成[4,9-10]：

其中：0表示射线初始强度，是射线衰减后强度，分别以实际成像中的背景灰度和图像灰度等效，表示X射线穿过的路径，(，)表示路径上处，能量为的射线的衰减系数，()表示归一化的等效能谱，满足：

与X射线初射能谱、探测器的闪烁体等有关[9]。式(1)中积分为线性运算，离散变为求和，等式两边同除以0：

其中：E为第个窄能谱段能量，(E)为第个窄能谱图像的加权系数，为物体中材质种类，u为第种材质在E的衰减系数，d为X射线穿过路径上第种材质的长度。其中，(E)满足：

当物体与成像系统几何位置确定后，d固定，能量划分固定后，u固定，只(E)随电压改变，因此，可将窄能谱图像看作基图像，不同电压X射线图像视作基图像的不同线性叠加。

单个电压下X射线图像是多个窄能谱图像以未知方式混合的结果。由于混合方式的未知，从一个电压下图像中无法获取窄能图像信息。因此，借鉴盲源分离思想，采集多个电压的X射线图像，通过不同的混合方式表现各窄能谱图像更多的信息，构造合适的目标函数，对其优化实现对各窄能谱投影的估计[11]。

将第个电压下X射线图像的第个像素对应的/0记为f，=1,2,…,，=1,2,…,。由式(3)，个电压的X射线成像写成矩阵形式为

其中：=(f)，=(s)，=(u)，=(d)。的行是构成各电压X射线图像(除以背景后)的窄能谱图像加权系数。的一列对应一材质在各窄谱段的衰减系数，的相应行是同种材料在各像素上的衰减长度。物体中由多成分均匀混合的材料视为一种材料。的每一行为一个窄能谱投影，为最终所求投影。

2 X射线图像分解模型

在式(5)中，为已知量，其每一个元素对应一个方程，共个方程，，，为未知量，共有++个未知量，要求式(5)对应的方程个数大于未知量的个数，即

，，都是非负矩阵，式(5)是一个特殊的非负矩阵分解模型。非负矩阵分解中通常是将问题转化为优化模型，可利用Karush-Kuhn-Tucker (KKT)条件求解模型[12-13]。式(5)中>，不能直接将其转化成两个非负矩阵分解问题。参照非负矩阵求解方法，考虑到实际成像中存在误差[9]，将式(5)的求解转化成以误差平方和最小为目标函数的优化模型求解，如下：

同样借助于KKT条件可得式(7)的迭代求解公式为

其中：“°”表示矩阵的Hadamard乘积(对应元素相乘)，“¾”表示矩阵的Hadamard除法(对应元素相除)。

对于式(7)的解和，给一任意´可逆矩阵，则¢=，¢=-1满足¢¢=，即¢和¢也是式(7)的解。由于最终需要的为，而上述多解性不影响最终量，因此，可以不考虑这种多解性，但其会对迭代过程造成影响，为减少多解性的影响，在迭代初始化中，使的每一列值都是从大到小的，并且每一次迭代中对的列进行归一化[13]。依据式(4)，行和为1，迭代中对的行归一化。

，，是求解前确定的量，而值需要在求解时设定。值越大，成像中能谱划分越窄，但极窄的 能谱是不必要的，并且会增加求解的计算量，因此，值不需要特别大。得到和，计算得新投影。

综上，式(5)求解算法流程图如图1，迭代停止准则为相邻两次迭代的目标函数值相对差异小于。

3 实验结果及分析

选取衰减系数接近的铝与硅构成的圆柱体进行验证实验如图2，铝圆柱直径30 mm，硅圆环内径30 mm，外径40 mm。依据美国国家标准技术研究院(NIST)数据库，铝与硅的线性衰减系数曲线如图3所示，两者在40 keV到140 keV范围内衰减系数差异小于10%，在能谱较宽的单电压CT和双能CT中无法有效区分。

实验的CT系统由450 kV的GE射线源(ISOVOLT-450)，选用焦斑尺寸为3.6 mm，16位的PerkinElmer平板探测器(XRD1621 AN14 ES)组成。平板探测器中心距射线源140 cm，物体旋转中心距探测器中心20 cm。依据圆柱体直径，电压选择120 kV、130 kV、140 kV，电流3 mA。投影间隔2º，共180个角度。为与论文方法所得投影图比较，给出了120 kV和140 kV对数变换后的投影图及其一行的灰度变化，如图4(a)和(b)。

期望以10 keV为一个窄能谱段，依据最高140 kV电压，划分为14个能谱段，即=14，取=0.01。

与基本非负矩阵分解类似，模型求解可能会落入局部最优点，因此，，初始值可依据先验信息设置，使其靠近真解，若随机选取初始值，需重复多次，选择最优解。以[0,1]随机值为初始值，依据先验，有中120 kV、130 kV下缺失的窄能谱图像的加权系数为0，对的每列值按从大到小重排，进行10 000次重复，最小目标函数值为2.101 8，对应迭代次数为73。选择分解后噪声最小，中等，最大的三个投影，如图4(c)、(d)、(e)，同时给出了其中间一行的灰度变化。

图4 投影及标记处的灰度曲线

对比图4各投影图，分解后的投影噪声比原始投影增大，但灰度变化趋势与原始投影一致，即结构信息未改变，保证了两者重建图像结构的一致。同时，分解所得投影灰度变化的陡峭程度较原始投影发生了改变，并且，这种变化并不是原始投影的等比例扩大或缩小，这意味着重建后图像的对比度要发生改变。

采用解析法重建图像。为减小噪声对结果比较的影响，对重建图像做中值滤波，滤波窗口5´5。图4(a)、(b)、(c)、(d)各投影的重建图像及其中间一行的灰度变化如图5。图4(e)噪声太大，不再对其重建。

由于射线硬化影响，重建图像出现硬化伪影，同种材质灰度值表现为中心低，边缘高，图5(a)和(b)在铝和硅界面上呈现出边缘增强现象，图5(c)和(d)由于硬化伪影减弱，边缘增强现象消失。对比图5各图像，(c)、(d)中硅铝差异比(a)、(b)中明显，(a)、(b)受硬化伪影影响，同种材质灰度差异较大，铝与硅的灰度重叠，无法依据灰度区分硅与铝，硅和铝对比度较小，而图(c)、(d)中铝明显高于硅的灰度，两者灰度无重叠，差异明显。按照：

图5 重建图像及标记处的灰度曲线

计算图5各图像铝与硅的对比度，其中Al和Si分别为铝材质区域和硅材质区域的平均灰度。各图像对比度如表1，相比于直接重建图像的对比度0.026 9，利用新投影重建的图像对比度有明显提高，达到0.2122。这说明利用对多电压X射线图像分解获取的新投影重建的图像，提高了对衰减特性相近材质的区分能力。

Table 1 Contrast of aluminium and silicon of the reconstructed images in Fig.5

4 结 论

针对常规CT系统X射线是连续宽能谱造成重建图像中衰减特性相近组分对比度较低的问题，论文提出了一种分解多电压X射线图像获取新投影的方法，用新投影重建的图像，有效地提高了对衰减系数相近材质的区分能力。该方法无需改变系统硬件，依据X射线成像规律，以图像分解为窄能谱图像的总误差平方和最小为目标建立优化模型，借鉴非负矩阵分解算法求解。在外硅内铝圆柱CT成像实验中，用新投影重建的图像可有效区分硅和铝，相比直接重建图像，对比度明显提高。进一步需研究新投影与真正窄能谱投影的差异，改进模型，获取更准确的窄能谱投影及其对应能量范围，在常规CT系统实现定量CT测量。

[1] JOSHUA D E，BRUCE R W，DAVID G P，. Experimental implementation of a polyenergetic statistical reconstruction algorithm for a commercial fan-beam CT scanner [J]. Physica Medica(S1120-1797)，2013，29(5)：500-512.

[2] TSUCHIY A，UESUGI K，NAKANO T，. Quantitative evaluation of attenuation contrast of X-ray computed tomography images using monochromatized beams [J]. American Mineralogist(S0003-004X)，2005，90(1)：132-142.

[3] MASETTI S，FIASCHETTI M，TURCO A，. Development of a Multi-Energy CT for Small Animals：Characterization of the Quasi-Monochromatic X-Ray Source [J]. IEEE Transactions on Nuclear Science(S0018-9499)，2009，56(1)：29-35.

[4] SEMERCI O，HAO N，KILMER M E，. Tensor-based formulation and nuclear norm regularization for multi-energy computed tomography [J]. IEEE Transactions on Image Processing(S1057-7149)，2014，23(4)：1678-1693.

[5] 郝佳，张丽，邢宇翔，等. 基于同步辐射光源的双能CT成像方法 [J]. 清华大学学报，2011，51(4)：457-461.
HAO Jia，ZHANG Li，XING Yuxiang，. Dual-energy CT imaging method using synchrotron radiation [J]. Journal of Tsinghua University，2011，51(4)：457-461.

[6] Cynthia H McCollough，LENG Shuai，YU Lifeng，. Dual- and Multi-Energy CT：Principles，Technical Approaches，and Clinical Applications [J]. Radiology(S0033-8419)，2015，276(3)：637–653.

[7] RAKVONGTHAI Y，WORSTELL W，FAKHRI G，. Spectral CT Using Multiple Balanced K-Edge Filters [J]. IEEE Transactions on Medical Imaging(S0278-0062)，2015，34(3)：740-747.

[8] 牛素鋆，潘晋孝，陈平. 基于能谱滤波分离的多谱计算机层析成像方法 [J]. 光学学报，2014，34(10)：1034001-1-1034001-7.
NIU Suyun，PAN Jinxiao，CHEN Ping. Multi-Spectrum Computed Tomography Imaging Method Based on Energy Spectrum Filtering Separation [J]. Acta Optica Sinica，2014，34(10)：1034001-1-1034001-7.

[9] 张慧滔，张朋. 基于计算层析成像扫描数据的X射线能谱估计方法 [J]. 光学学报，2013，33(11)：1134001-1-1134001-9.
ZHANG Huitao，ZHANG Peng. X-Ray Spectrum Estimation Method from Scanning Data of Computed Tomography Phantoms [J]. Acta Optica Sinica，2013，33(11)：1134001-1-1134001-9.

[10] GAO Hao，YU Hengyong，Osher Stanley，. Multi-energy CT based on a prior rank，intensity and sparsity model (PRISM) [J]. Inverse Problems(S0266-5611)，2011，27(11)：115012-115033.

[11] KEZIOU A，FENNIRI H，GHAZDALI A，. New blind source separation method of independent/ dependent sources [J]. Signal Processing(S0165-1684)，2014，104(6)：319-324.

[12] GILLIS N，GLINEUR F. Using underapproximations for sparse nonnegative matrix factorization [J]. Pattern Recognition (S0031-3203)，2009，43(4)：1676-1687.

[13] XU Wei，LIU Xin，GONG Yihong. Document clustering based on non-negative matrix factorization [C]// Proceedings of the 26th annual international ACM SIGIR conference on Research and development in information retrieval，New York，July 28–August 1，2003：267-273.

Applying Decomposition of Multi-voltage X-ray Images to Accomplish
High Contrast Computed Tomography Imaging

WEI Jiaotong，CHEN Ping，HAN Yan

( Shanxi Key Laboratory of Signal Capturing & Processing, North University of China, Taiyuan 030051, China )

The X-ray tube of usual computed tomography imaging system produces X-ray beam with a polychromatic spectrum, which results in low contrast of the materials with approximate attenuation coefficients in the reconstructed image. A method was proposed to obtain a new projection by decomposition of multi-voltage X-ray images. The distinction of the materials was enlarged in the reconstructed image by the new projection. In this method, the X-ray image of each voltage was considered as a sum of many narrow-energy-width X-ray images according to the law of X-ray imaging. Referring to non-negative matrix factorization, the new projection was got by minimizing the sum of squared error of multi-voltage X-ray images’ decomposition. A cylinder composed of silicon and aluminum was used in the experiment. The contrast of silicon and aluminum is evidently improved in the new reconstructed image compared to the direct reconstructed image. It shows the effectiveness of the proposed method.

information processing; computed tomography imaging; X-ray images; multi-voltage; non-negative matrix factorization; contrast

1003-501X(2016)08-0059-05

TP391

A

10.3969/j.issn.1003-501X.2016.08.010

2015-12-03；

2016-03-11

国家自然科学基金(61227003, 61301259, 61471325, 61571404)；山西省自然科学基金(2015021099)；山西省研究生优秀创新项目(20143081)

魏交统(1990-)，男(汉族)，山东临沂人。博士研究生，主要研究CT图像重建。E-mail:wei92837465@126.com。