魏交統(tǒng)，陳 平，韓 焱

?

多電壓X射線圖像分解的高對(duì)比度CT成像

魏交統(tǒng)，陳 平，韓 焱

( 中北大學(xué) 信息探測(cè)與處理山西省重點(diǎn)實(shí)驗(yàn)室，太原 030051 )

針對(duì)常規(guī)CT系統(tǒng)中X射線是連續(xù)寬能譜而導(dǎo)致重建圖像中衰減特性相近組分對(duì)比度較低的問題，論文提出一種分解多電壓X射線圖像序列獲取新投影的方法，用新投影重建的圖像有效地提高了對(duì)衰減特性相近組分的區(qū)分能力。該方法依據(jù)X射線成像規(guī)律，將各電壓下圖像看作多個(gè)窄能譜圖像的疊加，借鑒非負(fù)矩陣分解，以分解的誤差平方和最小為目標(biāo)對(duì)多電壓X射線圖像分解。用硅鋁材質(zhì)的圓柱體實(shí)驗(yàn)，新投影重建的圖像相比直接重建圖像，硅與鋁對(duì)比度明顯提高，表明了方法的有效性。

信號(hào)處理；CT成像；X射線圖像；多電壓；非負(fù)矩陣分解；對(duì)比度

0 引 言

常規(guī)CT成像系統(tǒng)的X射線是連續(xù)寬能譜，導(dǎo)致CT圖像出現(xiàn)硬化偽影，對(duì)比度降低。重建圖像灰度不僅與材質(zhì)、能量有關(guān)，還與像素位置有關(guān)，使衰減特性相近的組分無法有效區(qū)分[1]。解決這類問題，理想情況是通過單能CT實(shí)現(xiàn)組分衰減系數(shù)與能量的唯一對(duì)應(yīng)，如同步輻射CT[2]，但目前實(shí)際工程應(yīng)用的CT系統(tǒng)尚難以實(shí)現(xiàn)單能X射線源[3]，一般通過雙能或多能譜CT提高圖像不同組分的對(duì)比度[4]。在醫(yī)療領(lǐng)域，雙能CT能有效提高軟組織與其他物質(zhì)的對(duì)比度[5]，如2009年GE公司推出的寶石能譜CT。為進(jìn)一步提高圖像質(zhì)量，人們引入了更多的能量信息，提出了多譜成像技術(shù)。近年來出現(xiàn)了具有單光子能量分辨能力的光子計(jì)數(shù)型探測(cè)器。光子計(jì)數(shù)型探測(cè)器能將接收到的光子按能量分為不相交的能帶區(qū)間，進(jìn)行窄能譜CT成像，實(shí)現(xiàn)組分的定量分析[6]。然而該探測(cè)器時(shí)間分辨力受限，工程應(yīng)用存在局限[7]。借鑒這種多窄能譜成像思想，牛素鋆等提出了基于能譜濾波分離的多譜CT成像方法，添加不同濾波片獲取多個(gè)窄能譜CT圖像，但實(shí)際中添加濾波片使能譜變窄程度有限，不同能譜區(qū)分度不夠，無法實(shí)現(xiàn)組分的有效分離[8]。

借鑒光子計(jì)數(shù)型探測(cè)器分能譜段成像的思想，論文對(duì)多個(gè)電壓下X射線圖像分解，以期獲取窄能譜圖像。該方法不需要硬件條件的改變，將單電壓下X射線圖像看作多個(gè)窄能譜圖像的加權(quán)和，不同電壓下圖像僅是對(duì)應(yīng)的窄能譜圖像加權(quán)系數(shù)的不同。以各電壓下X射線圖像分解的總誤差平方和最小為目標(biāo)函數(shù)建立分解模型。求解模型獲取新投影，重建圖像，提高衰減系數(shù)相近組分的對(duì)比度。

1 變電壓X射線成像模型

X射線管發(fā)出的是連續(xù)能譜X射線，成像模型可寫成[4,9-10]：

其中：0表示射線初始強(qiáng)度，是射線衰減后強(qiáng)度，分別以實(shí)際成像中的背景灰度和圖像灰度等效，表示X射線穿過的路徑，(，)表示路徑上處，能量為的射線的衰減系數(shù)，()表示歸一化的等效能譜，滿足：

與X射線初射能譜、探測(cè)器的閃爍體等有關(guān)[9]。式(1)中積分為線性運(yùn)算，離散變?yōu)榍蠛停仁絻蛇呁?：

其中：E為第個(gè)窄能譜段能量，(E)為第個(gè)窄能譜圖像的加權(quán)系數(shù)，為物體中材質(zhì)種類，u為第種材質(zhì)在E的衰減系數(shù)，d為X射線穿過路徑上第種材質(zhì)的長(zhǎng)度。其中，(E)滿足：

當(dāng)物體與成像系統(tǒng)幾何位置確定后，d固定，能量劃分固定后，u固定，只(E)隨電壓改變，因此，可將窄能譜圖像看作基圖像，不同電壓X射線圖像視作基圖像的不同線性疊加。

單個(gè)電壓下X射線圖像是多個(gè)窄能譜圖像以未知方式混合的結(jié)果。由于混合方式的未知，從一個(gè)電壓下圖像中無法獲取窄能圖像信息。因此，借鑒盲源分離思想，采集多個(gè)電壓的X射線圖像，通過不同的混合方式表現(xiàn)各窄能譜圖像更多的信息，構(gòu)造合適的目標(biāo)函數(shù)，對(duì)其優(yōu)化實(shí)現(xiàn)對(duì)各窄能譜投影的估計(jì)[11]。

將第個(gè)電壓下X射線圖像的第個(gè)像素對(duì)應(yīng)的/0記為f，=1,2,…,，=1,2,…,。由式(3)，個(gè)電壓的X射線成像寫成矩陣形式為

其中：=(f)，=(s)，=(u)，=(d)。的行是構(gòu)成各電壓X射線圖像(除以背景后)的窄能譜圖像加權(quán)系數(shù)。的一列對(duì)應(yīng)一材質(zhì)在各窄譜段的衰減系數(shù)，的相應(yīng)行是同種材料在各像素上的衰減長(zhǎng)度。物體中由多成分均勻混合的材料視為一種材料。的每一行為一個(gè)窄能譜投影，為最終所求投影。

2 X射線圖像分解模型

在式(5)中，為已知量，其每一個(gè)元素對(duì)應(yīng)一個(gè)方程，共個(gè)方程，，，為未知量，共有++個(gè)未知量，要求式(5)對(duì)應(yīng)的方程個(gè)數(shù)大于未知量的個(gè)數(shù)，即

，，都是非負(fù)矩陣，式(5)是一個(gè)特殊的非負(fù)矩陣分解模型。非負(fù)矩陣分解中通常是將問題轉(zhuǎn)化為優(yōu)化模型，可利用Karush-Kuhn-Tucker (KKT)條件求解模型[12-13]。式(5)中>，不能直接將其轉(zhuǎn)化成兩個(gè)非負(fù)矩陣分解問題。參照非負(fù)矩陣求解方法，考慮到實(shí)際成像中存在誤差[9]，將式(5)的求解轉(zhuǎn)化成以誤差平方和最小為目標(biāo)函數(shù)的優(yōu)化模型求解，如下：

同樣借助于KKT條件可得式(7)的迭代求解公式為

其中：“°”表示矩陣的Hadamard乘積(對(duì)應(yīng)元素相乘)，“?”表示矩陣的Hadamard除法(對(duì)應(yīng)元素相除)。

對(duì)于式(7)的解和，給一任意′可逆矩陣，則￠=，￠=-1滿足￠￠=，即￠和￠也是式(7)的解。由于最終需要的為，而上述多解性不影響最終量，因此，可以不考慮這種多解性，但其會(huì)對(duì)迭代過程造成影響，為減少多解性的影響，在迭代初始化中，使的每一列值都是從大到小的，并且每一次迭代中對(duì)的列進(jìn)行歸一化[13]。依據(jù)式(4)，行和為1，迭代中對(duì)的行歸一化。

，，是求解前確定的量，而值需要在求解時(shí)設(shè)定。值越大，成像中能譜劃分越窄，但極窄的 能譜是不必要的，并且會(huì)增加求解的計(jì)算量，因此，值不需要特別大。得到和，計(jì)算得新投影。

綜上，式(5)求解算法流程圖如圖1，迭代停止準(zhǔn)則為相鄰兩次迭代的目標(biāo)函數(shù)值相對(duì)差異小于。

3 實(shí)驗(yàn)結(jié)果及分析

選取衰減系數(shù)接近的鋁與硅構(gòu)成的圓柱體進(jìn)行驗(yàn)證實(shí)驗(yàn)如圖2，鋁圓柱直徑30 mm，硅圓環(huán)內(nèi)徑30 mm，外徑40 mm。依據(jù)美國(guó)國(guó)家標(biāo)準(zhǔn)技術(shù)研究院(NIST)數(shù)據(jù)庫，鋁與硅的線性衰減系數(shù)曲線如圖3所示，兩者在40 keV到140 keV范圍內(nèi)衰減系數(shù)差異小于10%，在能譜較寬的單電壓CT和雙能CT中無法有效區(qū)分。

實(shí)驗(yàn)的CT系統(tǒng)由450 kV的GE射線源(ISOVOLT-450)，選用焦斑尺寸為3.6 mm，16位的PerkinElmer平板探測(cè)器(XRD1621 AN14 ES)組成。平板探測(cè)器中心距射線源140 cm，物體旋轉(zhuǎn)中心距探測(cè)器中心20 cm。依據(jù)圓柱體直徑，電壓選擇120 kV、130 kV、140 kV，電流3 mA。投影間隔2o，共180個(gè)角度。為與論文方法所得投影圖比較，給出了120 kV和140 kV對(duì)數(shù)變換后的投影圖及其一行的灰度變化，如圖4(a)和(b)。

期望以10 keV為一個(gè)窄能譜段，依據(jù)最高140 kV電壓，劃分為14個(gè)能譜段，即=14，取=0.01。

與基本非負(fù)矩陣分解類似，模型求解可能會(huì)落入局部最優(yōu)點(diǎn)，因此，，初始值可依據(jù)先驗(yàn)信息設(shè)置，使其靠近真解，若隨機(jī)選取初始值，需重復(fù)多次，選擇最優(yōu)解。以[0,1]隨機(jī)值為初始值，依據(jù)先驗(yàn)，有中120 kV、130 kV下缺失的窄能譜圖像的加權(quán)系數(shù)為0，對(duì)的每列值按從大到小重排，進(jìn)行10 000次重復(fù)，最小目標(biāo)函數(shù)值為2.101 8，對(duì)應(yīng)迭代次數(shù)為73。選擇分解后噪聲最小，中等，最大的三個(gè)投影，如圖4(c)、(d)、(e)，同時(shí)給出了其中間一行的灰度變化。

圖4 投影及標(biāo)記處的灰度曲線

對(duì)比圖4各投影圖，分解后的投影噪聲比原始投影增大，但灰度變化趨勢(shì)與原始投影一致，即結(jié)構(gòu)信息未改變，保證了兩者重建圖像結(jié)構(gòu)的一致。同時(shí)，分解所得投影灰度變化的陡峭程度較原始投影發(fā)生了改變，并且，這種變化并不是原始投影的等比例擴(kuò)大或縮小，這意味著重建后圖像的對(duì)比度要發(fā)生改變。

采用解析法重建圖像。為減小噪聲對(duì)結(jié)果比較的影響，對(duì)重建圖像做中值濾波，濾波窗口5′5。圖4(a)、(b)、(c)、(d)各投影的重建圖像及其中間一行的灰度變化如圖5。圖4(e)噪聲太大，不再對(duì)其重建。

由于射線硬化影響，重建圖像出現(xiàn)硬化偽影，同種材質(zhì)灰度值表現(xiàn)為中心低，邊緣高，圖5(a)和(b)在鋁和硅界面上呈現(xiàn)出邊緣增強(qiáng)現(xiàn)象，圖5(c)和(d)由于硬化偽影減弱，邊緣增強(qiáng)現(xiàn)象消失。對(duì)比圖5各圖像，(c)、(d)中硅鋁差異比(a)、(b)中明顯，(a)、(b)受硬化偽影影響，同種材質(zhì)灰度差異較大，鋁與硅的灰度重疊，無法依據(jù)灰度區(qū)分硅與鋁，硅和鋁對(duì)比度較小，而圖(c)、(d)中鋁明顯高于硅的灰度，兩者灰度無重疊，差異明顯。按照：

圖5 重建圖像及標(biāo)記處的灰度曲線

計(jì)算圖5各圖像鋁與硅的對(duì)比度，其中Al和Si分別為鋁材質(zhì)區(qū)域和硅材質(zhì)區(qū)域的平均灰度。各圖像對(duì)比度如表1，相比于直接重建圖像的對(duì)比度0.026 9，利用新投影重建的圖像對(duì)比度有明顯提高，達(dá)到0.2122。這說明利用對(duì)多電壓X射線圖像分解獲取的新投影重建的圖像，提高了對(duì)衰減特性相近材質(zhì)的區(qū)分能力。

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

4 結(jié) 論

針對(duì)常規(guī)CT系統(tǒng)X射線是連續(xù)寬能譜造成重建圖像中衰減特性相近組分對(duì)比度較低的問題，論文提出了一種分解多電壓X射線圖像獲取新投影的方法，用新投影重建的圖像，有效地提高了對(duì)衰減系數(shù)相近材質(zhì)的區(qū)分能力。該方法無需改變系統(tǒng)硬件，依據(jù)X射線成像規(guī)律，以圖像分解為窄能譜圖像的總誤差平方和最小為目標(biāo)建立優(yōu)化模型，借鑒非負(fù)矩陣分解算法求解。在外硅內(nèi)鋁圓柱CT成像實(shí)驗(yàn)中，用新投影重建的圖像可有效區(qū)分硅和鋁，相比直接重建圖像，對(duì)比度明顯提高。進(jìn)一步需研究新投影與真正窄能譜投影的差異，改進(jìn)模型，獲取更準(zhǔn)確的窄能譜投影及其對(duì)應(yīng)能量范圍，在常規(guī)CT系統(tǒng)實(shí)現(xiàn)定量CT測(cè)量。

[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，F(xiàn)IASCHETTI 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]. 清華大學(xué)學(xué)報(bào)，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，F(xiàn)AKHRI G，. Spectral CT Using Multiple Balanced K-Edge Filters [J]. IEEE Transactions on Medical Imaging(S0278-0062)，2015，34(3)：740-747.

[8] 牛素鋆，潘晉孝，陳平. 基于能譜濾波分離的多譜計(jì)算機(jī)層析成像方法 [J]. 光學(xué)學(xué)報(bào)，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] 張慧滔，張朋. 基于計(jì)算層析成像掃描數(shù)據(jù)的X射線能譜估計(jì)方法 [J]. 光學(xué)學(xué)報(bào)，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，F(xiàn)ENNIRI 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

國(guó)家自然科學(xué)基金(61227003, 61301259, 61471325, 61571404)；山西省自然科學(xué)基金(2015021099)；山西省研究生優(yōu)秀創(chuàng)新項(xiàng)目(20143081)

魏交統(tǒng)(1990-)，男(漢族)，山東臨沂人。博士研究生，主要研究CT圖像重建。E-mail:wei92837465@126.com。