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

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

?

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

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

( 中北大學 信息探測與處理山西省重點實驗室，太原 030051 )

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

信號處理；CT成像；X射線圖像；多電壓；非負矩陣分解；對比度

0 引 言

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

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

1 變電壓X射線成像模型

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

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

與X射線初射能譜、探測器的閃爍體等有關[9]。式(1)中積分為線性運算，離散變?yōu)榍蠛?，等式兩邊同除?：

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

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

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

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

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

2 X射線圖像分解模型

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

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

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

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

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

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

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

3 實驗結果及分析

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

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

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

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

圖4 投影及標記處的灰度曲線

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

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

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

圖5 重建圖像及標記處的灰度曲線

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

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

4 結 論

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

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

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