李漢平 王鋒 艾憲蕓

?

編碼板成像系統(tǒng)MLEM算法優(yōu)化

李漢平 王鋒 艾憲蕓

（國民核生化災(zāi)害防護(hù)國家重點實驗室防化研究院 北京 102205）

為增強(qiáng)γ輻射編碼成像系統(tǒng)中最大似然期望最大化(Maximum Likelihood Expectation Maximization, MLEM)算法對噪聲的抑制，提高算法重建質(zhì)量，基于互補(bǔ)編碼板成像相減去除噪聲的思想，對原有MLEM算法進(jìn)行改進(jìn)，提出了MLEM互補(bǔ)算法，并引入修正因子對MLEM互補(bǔ)算法的收斂程度進(jìn)行控制。通過模擬及實驗方式驗證了MLEM互補(bǔ)算法的有效性，給出了不同情況下MLEM互補(bǔ)算法達(dá)到最優(yōu)值時修正因子的擬合曲線。結(jié)果表明，MLEM互補(bǔ)算法可有效抑制噪聲，提高重建圖像質(zhì)量。修正因子的擬合曲線可以有效確定MLEM算法的最佳值。

編碼板，圖像重建，噪聲抑制，最大似然期望最大化法，蒙特卡羅模擬

在γ射線編碼板成像系統(tǒng)中，重建算法的好壞直接影響重建圖像的質(zhì)量。對于修正均勻冗余陣列(Modified Uniformly Redundant Arrays, MURA)編碼模式下的編碼板成像系統(tǒng)，研究人員對一些重建算法進(jìn)行過對比研究[1]，而最大似然期望最大化(Maximum Likelihood Expectation Maximization, MLEM)算法由于能對噪聲進(jìn)行很好的抑制，且在投影數(shù)據(jù)不完全時可以對圖像進(jìn)行很好的重建而成為主流算法。MLEM算法的缺點是：在編碼板投影圖像所含噪聲較大時，隨著迭代次數(shù)的增加，噪聲對重建圖像的影響會相應(yīng)放大。目前對于MLEM算法的改進(jìn)主要集中在核醫(yī)學(xué)領(lǐng)域，如有序子集期望最大化(Ordered Subset Expectation Maximization, OSEM)算法、子集序列期望最大化(Subset Sequence Expectation Maximization, SSEM)算法以及計算調(diào)控有序子集期望最大化(Count Regulated Ordered Subset Expectation Maximization, CROSEM)算法等[2?3]。由于編碼板成像與醫(yī)學(xué)成像在數(shù)據(jù)采集方式上的不同，醫(yī)學(xué)成像中的相關(guān)算法不能直接應(yīng)用于編碼板成像。在編碼板成像中，有關(guān)MLEM算法改進(jìn)的論述較少。

對于編碼板成像系統(tǒng)，噪聲大小直接影響其圖像質(zhì)量。為了抑制噪聲對重建圖像的影響，通過理論證明[4]，可以用兩個互補(bǔ)的編碼孔經(jīng)重建所得圖像相減來減小噪聲對重建圖像的影響。本文基于該思想對編碼板成像中的MLEM算法進(jìn)行改進(jìn)，提出MLEM互補(bǔ)算法。該算法在傳統(tǒng)MLEM算法公式中加入編碼板修正因子，在迭代算法中進(jìn)行互補(bǔ)編碼板相減計算。MLEM互補(bǔ)算法中最佳修正因子值的選取與編碼板線性衰減系數(shù)與厚度的乘積有關(guān)。通過蒙特卡羅模擬給出算法最佳修正因子的t擬合公式，并進(jìn)行實驗驗證。實驗結(jié)果證明，與傳統(tǒng)MLEM算法相比，MLEM互補(bǔ)算法在成像時可以更有效地抑制噪聲。采用模擬所得-擬合公式，可以給出不同情況下MLEM互補(bǔ)算法最佳值。

1 算法簡介

1.1 MLEM算法介紹

MLEM算法是基于泊松統(tǒng)計模型的最大貝葉斯后驗概率的圖像復(fù)原方法，該算法由Shepp等[5]于1982年提出。MLEM算法主要分為兩步：1) 計算完全投影數(shù)據(jù)的似然函數(shù)在測定數(shù)據(jù)及當(dāng)前參數(shù)估計值下的期望；2) 求出使期望最大化的參數(shù)估計值。

在γ輻射編碼成像中MLEM算法[6?8]表示為：

1.2 基于互補(bǔ)碼板的MLEM算法

基于互補(bǔ)編碼板相減去除噪聲的思想，在進(jìn)行MLEM迭代時，對迭代公式(1)中的編碼函數(shù)(,) 進(jìn)行修正。令(,)=′(,)+h″(,)，其中：′(,)為正常編碼函數(shù)；″(,)為′(,)逆時針旋轉(zhuǎn)90°所得編碼函數(shù)；為算法修正因子，取值范圍為(?0.8,0)。將(,)帶入迭代公式，可得：

在計算第次放射源強(qiáng)度分布的迭代估計值和編碼函數(shù)所計算得出的投影估計值時使用互補(bǔ)編碼板相減去噪的方式對噪聲進(jìn)行抑制。由式(2)可知，由于在計算時對進(jìn)行了噪聲抑制，使得實際測定值與修正理論投影值的比值對重建圖像的修正效果更好。將該值與修正后的(,)做相關(guān)運(yùn)算，使得迭代算法再次對噪聲進(jìn)行抑制，進(jìn)一步減少噪聲對重建圖像的影響。

2 模擬及實驗結(jié)果分析

2.1 模擬條件及模擬結(jié)果分析

使用蒙特卡羅模擬軟件MCNP5 (Monte Carlo N Particle Transport Code 5)對γ輻射編碼板成像過程進(jìn)行模擬，得到實驗所需編碼板投影數(shù)據(jù)。建立的γ輻射編碼板成像模型示意圖如圖1所示。編碼板以11×11擴(kuò)展模式的MURA編碼板為基礎(chǔ)，幾何尺寸為50mm×50mm×10mm，采用方孔模式，孔徑2mm，材料選用鉛；點源位于視野中心軸線，距編碼板500mm處，能量為662keV；探測器材料選用鍺酸鉍(Bi2O3-GeO2, BGO)閃爍晶體，幾何尺寸為50mm×50mm×15mm，距編碼板20mm，計數(shù)時，將BGO晶體模型按100×100×1進(jìn)行網(wǎng)格劃分，并采用*F8計數(shù)卡對每個網(wǎng)格內(nèi)沉積能量進(jìn)行記錄。通過蒙特卡羅模擬所得理想點源圖像及其投影圖像如圖2所示。

圖1 編碼板模型示意圖

圖2 理想點源(a)和點源投影(b)圖像

采用蒙特卡羅模擬所得編碼板投影數(shù)據(jù)，分別通過MLEM算法和MLEM互補(bǔ)算法對圖像進(jìn)行重建。圖3為MLEM算法與MLEM互補(bǔ)算法(=?0.6)在迭代次數(shù)=5時的重建圖像。由圖3可知，兩種算法都可對模擬編碼圖像進(jìn)行有效重建。采用MLEM互補(bǔ)算法對模擬圖像進(jìn)行重建，其對噪聲的抑制效果更好，重建圖像更接近模擬的點源圖像。通過直接觀察法可知，MLEM互補(bǔ)算法的成像效果要優(yōu)于傳統(tǒng)的MLEM算法。

對重建圖像進(jìn)行評估時，主觀評價往往很難判斷算法的優(yōu)劣。為進(jìn)一步驗證算法的有效性，在對重建圖像的質(zhì)量進(jìn)行評價時，采用歸一化均方誤差(Normalized Root Mean Square Error, NRMSE)以及相關(guān)系數(shù)(Correlation Coefficient, CC)這兩種定量方法對重建圖像的質(zhì)量進(jìn)行評價。用歸一化均方誤差衡量算法的接近程度，相關(guān)系數(shù)衡量重建圖像與原圖像相似程度[9]。NRMSE與CC評價標(biāo)準(zhǔn)的值和定義分別為：

表1 MLEM算法和MLEM互補(bǔ)算法NRMSE與CC值

2.2 實驗條件及實驗結(jié)果分析

為驗證模擬結(jié)果的正確性，通過實驗對算法的有效性進(jìn)行驗證。實驗中編碼板以11×11擴(kuò)展的MURA編碼模式為基礎(chǔ)，編碼板幾何尺寸分別為50mm×50mm×10mm、50mm×50mm×20mm、50mm×50mm×30mm，孔型為方孔，孔徑2mm，材料選用鉛；探測器由50mm×50mm×10mm的BGO陣列晶體，耦合濱松H8500位置靈敏光電倍增管組成，探測器距編碼板20mm。點源采用強(qiáng)度3.7×104Bq的137Cs源，位于視野中心軸線，距編碼板500mm處。

采用實驗所得編碼板投影數(shù)據(jù)，分別通過MLEM算法和MLEM互補(bǔ)算法對圖像進(jìn)行重建。圖4為MLEM算法和MLEM互補(bǔ)算法在迭代次數(shù)=5時對三種厚度編碼板實測投影重建圖像。由圖4(a)可知，隨著編碼板厚度的增加，MLEM算法重建圖像的背景噪聲越來越小，圖像質(zhì)量逐步提高。由圖4(b)可知，采用MLEM互補(bǔ)算法進(jìn)行圖像重建時，在編碼板厚度相同的情況下，MLEM互補(bǔ)算法的重建結(jié)果比MLEM算法的重建結(jié)果對噪聲的抑制更明顯。

圖4 MLEM算法(a)與MLEM互補(bǔ)算法(b)實驗數(shù)據(jù)重建圖像

表2為MLEM算法和MLEM互補(bǔ)算法在迭代次數(shù)=5時對三種厚度編碼板實測投影重建圖像的NRMSE和CC值。由表2可知，在編碼板厚度相同的情況下，MLEM互補(bǔ)算法的NRMSE和CC值皆優(yōu)于傳統(tǒng)的MLEM算法，實驗結(jié)果與模擬結(jié)果相一致。與MLEM算法相比，采用MLEM互補(bǔ)算法對編碼圖像進(jìn)行重建，可以更有效地抑制噪聲，提高重建圖像的質(zhì)量。

表2 MLEM算法和MLEM互補(bǔ)算法的NRMSE、CC值

2.3值對MLEM互補(bǔ)算法的影響

通過MLEM互補(bǔ)算法對上述實驗及模擬投影數(shù)據(jù)進(jìn)行重建時，=?0.6。為研究不同值對MLEM互補(bǔ)算法重建效果的影響，選取不同值對實驗所得厚度=10mm、=20mm、=30mm三種厚度編碼板的實測投影值進(jìn)行互補(bǔ)MLEM算法重建。圖5為三種厚度下互補(bǔ)MLEM算法的NRMSE和CC值隨修正因子值的變化。

圖5 三種厚度鉛制編碼板NRMSE (a)和CC (b)值隨a變化

由圖5可知，隨著值在(?0.8,0)區(qū)間內(nèi)不斷增大，重建圖像的質(zhì)量先隨之提高，當(dāng)達(dá)到某一特定值后，重建圖像的質(zhì)量開始下降。對于MLEM算法存在最佳值，使得改進(jìn)后MLEM算法效果最優(yōu)。編碼板厚度為1 cm、2 cm、3 cm時，最佳值分別為?0.7、?0.61、?0.55?？梢婋S著編碼板準(zhǔn)直器厚度的增加，MLEM互補(bǔ)算法修正因子的最佳值逐漸減小，最佳值的選取與物質(zhì)阻止本領(lǐng)有關(guān)。由γ射線在物質(zhì)中的衰減規(guī)律，采用上述編碼板模型，模擬不同t值下，算法最佳值的變化。圖6為最佳值隨t的變化。由圖6可知，最佳值隨t值一同增加，當(dāng)t值達(dá)到某一特定值時，值的變化趨于平緩。計算所得t擬合公式為：

=?0.6+0.12421?0.01252()2(6)

圖6 最佳a(bǔ)值隨mt值變化擬合曲線

表3為三種厚度下鉛制編碼板實測投影數(shù)據(jù)所得最佳值與模擬數(shù)據(jù)擬合公式計算所得最佳值。

表3 三種厚度鉛制編碼板實驗及模擬數(shù)據(jù)最佳a(bǔ)值

由表3可知，由于噪聲影響，實驗數(shù)據(jù)所得值比模擬數(shù)據(jù)所得的要大。無法通過模擬數(shù)據(jù)所得t擬合公式直接給出實測情況下MLEM互補(bǔ)算法最佳值。在使用t擬合公式給出實際成像系統(tǒng)所需最佳值時，需對模擬所得最佳值乘以相應(yīng)比例系數(shù)。

記實驗值與模擬值的比值為比例系數(shù)，根據(jù)t擬合公式變化趨勢，對隨t值變化進(jìn)行擬合。圖7為比例系數(shù)隨t值的變化。

圖7 比例系數(shù)N隨最佳mt值變化擬合曲線

計算所得t擬合公式為：

=1.0156+0.4462?0.0651()2(7)

在求MLEM互補(bǔ)算法的最佳值時，可以通過模擬所得t擬合公式計算結(jié)果乘以得到。

為驗證上述結(jié)論，采用1cm鎢板進(jìn)行實驗，圖8為鎢板的NRMSE和CC值隨值的變化。由圖8可知，此時MLEM互補(bǔ)算法最佳值為0.65。在入射γ射線能量為662keV時，鎢的線性衰減系數(shù)=1.84cm?1[10]；將=1cm帶入擬合公式(6)，得出模擬情況下算法最佳值為0.41，由擬合公式(7)，此時為1.6。最終用于1cm鎢制編碼板的最佳=1.6×0.41=0.65。與實驗測得數(shù)據(jù)所得最佳值相一致。在實際應(yīng)用中，可以用上述方法確定MLEM互補(bǔ)算法最佳值。

圖8 1cm厚鎢制編碼板NRMSE (a)和CC (b)值隨a變化

3 結(jié)語

實驗證明，將互補(bǔ)編碼板相減去噪思想應(yīng)用于MLEM算法中是可行的。通過對蒙特卡羅模擬所得投影數(shù)據(jù)及實測所得投影數(shù)據(jù)的重建結(jié)果可知，相比于MLEM算法，在編碼板成像中，MLEM互補(bǔ)算法可以有效抑制噪聲的影響，提高圖像重建質(zhì)量；MLEM互補(bǔ)算法中的最佳值的選取與t有關(guān)；通過模擬得出t擬合公式所得模擬最佳值與t擬合公式所得的乘積，可以得到實測數(shù)據(jù)的最佳值，模擬結(jié)果與實驗結(jié)果相一致。采用該方法可以確定MLEM互補(bǔ)算法最佳值；在編碼板成像系統(tǒng)中，采用蒙特卡羅模擬方法對其進(jìn)行分析具有一定的可行性。

1 趙翠蘭, 陳立宏, 李勇平. MURA編碼輻射成像系統(tǒng)的解碼方法[J]. 核技術(shù), 2014, 37(8): 080401. DOI: 10.11889/j.0253-3219.2014.hjs.37.080401.

ZHAO Cuilan, CHEN Lihong, LI Yongping. Decoding process of a radiation imaging system using MURA coded aperture collimator[J]. Nuclear Techniques, 2014, 37(8): 080401. DOI: 10.11889/j.0253-3219.2014.hjs.37.080401.

2 楊娟, 王明泉, 石浪, 等. OSEM重建算法及其改進(jìn)算法的研究和比較[J]. 計算機(jī)工程與設(shè)計, 2015, 36(9): 2524?2526. DOI: 10.16208/j.issn1000-7024.2015.09.040. YANG Juan, WANG Mingquan, SHI Lang,. Research and comparison on OSEM and its improved reconstruction algorithms[J]. Computer Engineering and Design, 2015, 36(9): 2524?2526. DOI: 10.16208/j. issn1000-7024.2015.09.040.

3 金永杰, 馬天予. 核醫(yī)學(xué)儀器與方法[M]. 哈爾濱: 哈爾濱工程大學(xué)出版社, 2010.

JIN Yongjie, MA Tianyu. Nuclear medical instrument and method[M]. Harbin: Harbin Engineering University Press, 2010.

4 Barrett H H, Swindell W. 放射成像[M]. 莊天戈, 周頌凱, 譯. 北京: 科學(xué)出版社, 1988.

Barrett H H, Swindell W. Radiation imaging[M]. ZHUANG Tiange, ZHOU Songkai, trans. Beijing: Science Press, 1988.

5 SheppL A,Vardi Y. Maximum likelihood reconstruction for emission tomography[J]. IEEE Transactions on Medical Imaging, 1982, 1(2):113?122. DOI: 10.1109/TMI.1982.4307558.

6 張斌, 王英, 艾憲蕓, 等. 基于約束的MLEM圖像重建算法[J]. 原子能科學(xué)技術(shù), 2014, 48(增1): 668?672. DOI: 10.7538/yzk.2014.48.S0.0668.

ZHANG Bin, WANG Ying, AI Xianyun,. MLEM image reconstruction algorithm based on physical constraints[J]. Atomic Energy Science and Technology, 2014, 48(Suppl 1): 668?672. DOI: 10.7538/yzk.2014.48. S0.0668.

7 洪俊杰. γ相機(jī)中編碼孔經(jīng)準(zhǔn)直器設(shè)計及數(shù)字圖像重建[D]. 湖北: 華中科技大學(xué), 2006. HONG Junjie. Design of coded aperture collimator and digital image reconstruction in gamma-ray camera[D]. Hubei: Huazhong University of Science and Technology, 2006.

8 Mu Z P, Liu Y H. Aperture collimation correction and maximum-likelihood image reconstruction for near-field coded aperture imaging of single photon emission computerized tomography[J]. IEEE Transactions on Medical Imaging, 2006, 25(6): 701?711. DOI: 10.1109/ TMI.2006.873298.

9 何佳偉, 劉東升, 桂志國. 可變有序子集PML算法在PET中的應(yīng)用[J]. 中北大學(xué)學(xué)報(自然科學(xué)版), 2010, 31(6): 646?650. DOI: 10.3969/j.issn.1673-3193.2010.06. 021.

HE Jiawei, LIU Dongsheng, GUI Zhiguo. Application of modified subsets PML algorithm to PET image reconstruction[J]. Journal of North University of China (Natural Science Edition), 2010, 31(6): 646?650. DOI: 10.3969/j.issn.1673-3193.2010.06.021.

10 方杰. 輻射防護(hù)導(dǎo)論[M]. 北京: 原子能出版社, 1991.FANG Jie. Introduction to radiation protection[M]. Beijing: Atomic Press, 1991.

Algorithm optimization of MLEM in coded aperture imaging system

LI Hanping WANG Feng AI Xianyun

(State Key Laboratory of NBC Protection for Civilian, Research Institute of Chemical Defense, Beijing 102205, China)

Background: Ingamma-ray imager, reconstruction algorithm directly affects the quality of the reconstructed image. The maximum likelihood expectation maximization (MLEM) iteration algorithm is widely used inmodified uniformly redundant arrays(MURA) coded aperture for the reason of satisfactory performance of suppressing noise, improving signal noise ratio (SNR) and reducing distortion. But MLEM iteration algorithm also has shortcomings, like amplifies the noise with the increase of the iterative times. Purpose: This study aims to improve the ability of suppressing noise for MLEM iteration algorithm, and increase precision and resolution of coded aperture imaging system.Methods:Based on the de-noising method of complementary coded aperture, a correction factorfor MLEM algorithm modification is proposed to control the convergence rate of the complementary MLEM iteration algorithm.Both Monte Carlo simulation and experimental date are used to verify the effectiveness of the complementary MLEM iteration algorithm. Results: The results of Monte Carlo simulation and experiment prove the complementary MLEM iteration algorithm is efficient and practicable in coded aperture image. Its efficiency is relatedto the modifying factor. The fitting formula of-inthe complementary MLEM iteration algorithm is=?0.6+0.12421?0.01252()2. Conclusions: The coded aperture gamma camera is capable of imaging gamma-rays instantly. It can form radioactive 2D-map of different radionuclides. The reconstruction algorithm is important for coded aperture gamma camera. With the development of coded aperture, more and more study of algorithm modification will focus on it.

Coded aperture, Reconstruction of image, Suppress noise, MLEM, Monte Carlo simulation

LI Hanping, male, born in 1992, graduated from Shanghai Jiao Tong University in 2014, master student, major in radiation protection and environmental protection

2016-11-20, accepted date: 2016-12-14

TL812

10.11889/j.0253-3219.2017.hjs.40.020404

李漢平，男，1992年出生，2014年畢業(yè)于上海交通大學(xué)，現(xiàn)為碩士研究生，輻射防護(hù)與環(huán)境保護(hù)專業(yè)

2016-11-20，

2016-12-14