編碼成像MLEM互補算法的參數(shù)優(yōu)化

侯英偉， 宋玉收， 孫世杰， 柳若琦， 黃麗萍， 胡力元， 劉輝蘭

(1.哈爾濱工程大學(xué) 核科學(xué)與技術(shù)學(xué)院，黑龍江 哈爾濱 150001； 2.中國核動力研究設(shè)計院，四川 成都 610213)

基于編碼孔成像的伽馬相機可以提供環(huán)境中伽馬射線源的位置信息與能譜，快速對伽馬源項進行定位與分析，因此編碼成像技術(shù)已成功地應(yīng)用于核材料監(jiān)管和防擴散領(lǐng)域。MURA[1]編碼具有開孔率大(約50%)、在有限大小內(nèi)的自相關(guān)特性好、可實現(xiàn)中心反對稱的方形排列等優(yōu)點，被應(yīng)用在各種成像技術(shù)中[2]，同時也是現(xiàn)在伽馬編碼孔成像最常用的編碼方式。

當(dāng)使用編碼孔進行近場成像時，編碼孔成像的數(shù)學(xué)模型證明[3]只有在遠場條件下，實驗設(shè)置具有理想相關(guān)特性成像關(guān)系才完美符合點擴散函數(shù)。在近場條件下，點源所放出的伽馬射線不能看作是相互平行，因此源項在探測器上的投影具有不均勻性，經(jīng)典反卷積解碼后的重構(gòu)圖像會產(chǎn)生偽影[4]。最大似然期望最大化(maximum likelihood expectation maximization，MLEM)算法[5]可以進一步提高編碼孔成像質(zhì)量，但是當(dāng)投影圖像所含噪聲較大時，迭代后噪聲的影響也相應(yīng)放大；另一方面為消除近場成像所產(chǎn)生的偽像與噪聲，Accorsi等[6]提出了編碼孔成像的雙模測量法，采用硬件實現(xiàn)，不需要耗費大量計算時間，操作方式簡便不會增加成像難度。為進一步消除偽像與噪聲提高信噪比，李漢平等[7]提出基于互補碼的改進MLEM算法，該算法在傳統(tǒng)MLEM算法公式中加入編碼板修正因子[8]，在迭代算法中進行互補編碼板相減計算以達到消除噪聲的目的?；パa編碼MLEM雖然能夠有效去除偽像提高信噪比，但是互補碼方法的修正因子與迭代次數(shù)是重要的成像參數(shù)，如果修正因子與迭代次數(shù)設(shè)置不正確，反而會引入噪聲[9]，干擾源項位置的判斷。并且如果源項在視野范圍(field of view，F(xiàn)OV)內(nèi)的不同位置，最佳的參數(shù)設(shè)置有所不同。

本文基于MLEM互補算法，利用Geant4對編碼孔成像過程進行模擬，通過模擬不同源項位置的成像過程，研究互補編碼MLEM成像法的重構(gòu)圖像結(jié)果質(zhì)量與修正參數(shù)的關(guān)系，選取合適的修正因子大小，對不同源項位置所需迭代次數(shù)給出擬合經(jīng)驗公式并加以驗證。

1 基于互補編碼的MLEM成像算法

為減小由于投影不均勻性產(chǎn)生的重構(gòu)偽影，MLEM被應(yīng)用在編碼孔成像中。MLEM算法是基于泊松統(tǒng)計Poisson-Likehood概率函數(shù)對數(shù)最大值的圖像重建方法，MLEM算法為：

(1)

式中：fk(x，y)表示第k次迭代后放射源圖像的估計值；y(x，y)表示探測器得到的投影數(shù)據(jù)；h(x，y)表示編碼函數(shù)。迭代過程中的初始迭代估計值的取值有2種方式：1)f0(x,y)≡1；2)將相關(guān)解碼得到的重建圖像作為初始迭代估計值。

另一方面，成像時由于偽影的存在，將正編碼模式和反編碼模式對同一放射源進行編碼孔成像，并采用相同的相關(guān)算法去解碼重建圖像[10]。將2幅重建圖像的圖像相加，就能消除編碼孔近場成像時由于二階項導(dǎo)致的偽影，從而增強重建圖像的質(zhì)量。

在使用MLEM迭代算法進行圖像重建過程中，如果原始的編碼板投影結(jié)果包含較強的噪聲，隨著迭代次數(shù)的增加噪聲對圖像的影響也會相應(yīng)增大。為了抑制噪聲對重建圖像的影響，理論證明互補編碼去噪方法對MLEM方法同樣起作用。對迭代式(1)中的編碼函數(shù)進行了修正，引入了編碼板修正因子β。即：

h(x,y)=ha(x,y)+βhb(x,y)

(2)

式中：ha(x,y)是正編碼函數(shù)；hb(x,y)是正編碼函數(shù)旋轉(zhuǎn)90°后得到的反編碼函數(shù)，β為修正因子，取值范圍(-0.8, 0)。式(2)代入式 (1)中可得：

(3)

由式 (2)可知，互補編碼的MLEM成像算法是利用反編碼修正后的編碼函數(shù)，對迭代結(jié)果fk(x，y)進行投影后與原始投影相比較，并進行反投影運算。在“投影”與“反投影”過程中均利用反編碼函數(shù)對噪聲進行抑制，達到有效減少噪聲的目的。

2 編碼孔成像過程模擬

在利用互補編碼的MLEM成像算法進行迭代成像過程中，需要選取合適的重構(gòu)參數(shù)——迭代次數(shù)N與修正因子β，使得利用算法進行重構(gòu)后的圖像具有更好的圖像質(zhì)量。為了尋找合適的迭代參數(shù)，使用4.10.4版本Geant4進行蒙卡模擬，Geant4的PhysicsList采用QGSP_BERT，其包含所有的典型物理過程。

模擬采用137Cs伽馬點狀放射源，放出的伽馬能量為662 keV。MURA編碼板設(shè)置為19×19維擴展編碼板，單質(zhì)鎢作為編碼孔材料，鎢的原子序數(shù)是74，密度為19.30 g/cm3，其對能量為662 keV的伽馬射線的線性衰減系數(shù)約為1.7 cm-1，編碼板幾何厚度設(shè)置為1.5 cm。如圖1所示為本次模擬的編碼板具體設(shè)計形狀，單位小孔大小為2 mm×2 mm。

圖1 嵌套MURA (19×19)編碼板開孔情況Fig.1 MURA (19×19) coding board

CsI晶體作為位置靈敏探測器，單位編碼孔的采樣率為4尺寸為47.5 mm×47.5 mm×50 mm，即探測器在XOY平面被平均分為76×76個像素，通過Geant4的靈敏探測器功能獲取伽馬射線入射至CsI晶體的位置信息，并劃歸至對應(yīng)的CsI像素前端面中心位置，以獲得離散化的投影圖像。

將伽馬源所在位置設(shè)置為坐標(biāo)原點，伽馬點源、編碼板、探測器沿Z軸依次放置，源項所在位置即探測系統(tǒng)的FOV為坐標(biāo)系的XOY平面，原點為探測器與編碼板的中心投影位置。伽馬源距編碼板中心80 cm，編碼板中心與探測器前端面的距離為20 cm。投影成像的放大倍數(shù)為1.25倍，此編碼板探測系統(tǒng)的FOV為19 cm×19 cm。圖2為編碼孔成像模擬的幾何設(shè)置示意圖。

圖2 編碼孔成像幾何設(shè)置示意Fig.2 Diagrammatic of coding aperture imaging

3 成像結(jié)果分析與優(yōu)化方法

以極坐標(biāo)觀察伽馬源所在XOY平面，改變伽馬源在成像平面的位置，將源項分別設(shè)置于極角θ為0°、15°、30°、45°處，距原點距離R即源項所在極徑分別設(shè)置為50、80、160、200、240和290 mm處分別進行成像模擬。使用基于互補編碼板的MLEM方法進行源項位置重構(gòu)，測量不同源項位置合適迭代次數(shù)N與修正因子β的關(guān)系。

3.1 成像結(jié)果質(zhì)量評價參數(shù)

為更好地判斷重構(gòu)參數(shù)的選取是否正確，引入峰噪比(peak contrast to noise ratio，PCNR)與源項重構(gòu)強度的半高全寬(full width half maximum，F(xiàn)WHM)作為成像質(zhì)量的評價參數(shù)。峰噪比PCNR為：

(4)

式中：S為重構(gòu)圖像中計數(shù)強度的最大值；B為圖像背景均值；σ為圖像背景標(biāo)準(zhǔn)差。

PCNR越大證明圖像的峰值相對于背景噪聲的對比度越大，在成像結(jié)果源項定位位置正確的前提下PCNR間接表示了圖像的成像質(zhì)量。在進行互補編碼板的MLEM迭代計算過程中，為選取合適的迭代次數(shù)N與修正因子β，選取峰噪比PCNR與FWHM作為判斷標(biāo)準(zhǔn)。如圖3為源項位置位于(141.4, 141.4)時，不同修正因子β設(shè)置條件下FWHM與PCNR隨迭代次數(shù)增加的變化情況，不同形狀點代表不同的β取值，每個圖像的右側(cè)端點為首次迭代開始點。

由圖3可知無論選取多大的修正因子，隨著迭代次數(shù)的增加，PCNR逐漸變大與FWHM逐漸變小也就是重構(gòu)圖像的質(zhì)量逐漸提高。

圖3 不同修正因子β設(shè)置條件下FWHM與PCNR隨迭代次數(shù)增加的變化情況Fig.3 With different correction factors β, the change of FWHM and PCNR with the increase of iteration number

隨著迭代次數(shù)的增加，圖像質(zhì)量的提升幅度尤其是FWHM變小的幅度越來越小，過多的迭代次數(shù)會極大降低計算效率消耗計算性能，因此計算過程中不推薦使用過多的迭代次數(shù)。在某些修正因子下當(dāng)重構(gòu)圖像質(zhì)量提升至一定大小后，如圖3中β>0.4情況，圖4所示繼續(xù)進行迭代計算將導(dǎo)致圖像出現(xiàn)紊亂，使得重構(gòu)的源項位置出現(xiàn)偏差，并且圖像質(zhì)量降低，甚至圖像變形。圖4為源項位置位于(141.4, 141.4)，β取值為-0.5時不同迭代次數(shù)的成像結(jié)果。

圖4 不同迭代次數(shù)的成像結(jié)果Fig.4 Imaging results of different iterations

為避免成像結(jié)果的質(zhì)量出現(xiàn)震蕩及提高計算效率，在模擬伽馬源成像的過程中，設(shè)置當(dāng)成像結(jié)果第一次達到PCNR>30時，便達到了合適的迭代次數(shù)。由于本文算研究內(nèi)容的均是單點源，PCNR是對整個圖像的評價參數(shù)，而FWHM是對重構(gòu)的源項位置分辨的評價參數(shù)，所以沒有將FWHM參數(shù)作為迭代次數(shù)選取的限制條件。

3.2 迭代修正參數(shù)的計算方法

為研究修正參數(shù)與源項所在極徑的關(guān)系，將同一極角的數(shù)據(jù)進行對比分析。圖5為源項所在位置的極徑不同時，最佳迭代次數(shù)N隨著修正因子β取值的變化情況。圖5(a)與(b)分別為源項所在位置極角0°與30°，不同形狀點為不同源項所在位置的極徑大小。由圖5可知，當(dāng)采用的修正因子絕對值太小時，如圖中當(dāng)β>-0.3時，噪聲的抑制效果不明顯，因此需要更多的迭代次數(shù)才能獲得較好的成像結(jié)果。需要注意的是，過多的迭代次數(shù)不僅有可能放大噪聲對成像的影響，而且增加了成像過程的不穩(wěn)定性。當(dāng)源項所在位置極徑較大位于FOV邊緣時，源項在探測器上的投影更不均勻，相應(yīng)重構(gòu)后的偽影也較強。因此當(dāng)修正參數(shù)絕對值較低時較弱的去噪效果不能完全消除偽影，無論經(jīng)過多少次迭代計算也無法獲得合適的重構(gòu)圖像結(jié)果。如圖5中當(dāng)R=290 mm，β>-0.3時，無法達到PCNR>30質(zhì)量要求，沒有有效的數(shù)據(jù)點。

圖5 源項所在位置的極徑不同時，最佳迭代次數(shù)N隨著修正因子β的變化情況Fig.5 When the polar diameter of the source is different, the optimal number of iterations N changes with the correction factor β

當(dāng)修正因子的絕對值很大時，偽影消除效果明顯，所需的迭代次數(shù)也相應(yīng)減少。但如圖3所示過強的偽影消除效果，導(dǎo)致迭代次數(shù)增加至一定數(shù)量反而使得成像質(zhì)量出現(xiàn)震蕩，甚至無法成像定位。且每次迭代結(jié)果相差較大不易控制迭代次數(shù)，不易尋找到普適性的迭代次數(shù)。通過對比成像結(jié)果可知，由于偽像主要由源項在探測器上投影不均勻造成的，因此相比于極徑，角度對參數(shù)選取的影響較小。

總的來說，相比于源項所在極角，極徑對于迭代參數(shù)的選取影響更大；不應(yīng)選取絕對值過小的修正系數(shù)，會降低迭代計算效率以及無法重構(gòu)出較大的極徑位置的源項；也不應(yīng)選取絕對值過大的修正系數(shù)，使得所需迭代次數(shù)過小而不易控制成像質(zhì)量并且導(dǎo)致圖像質(zhì)量變差。

綜上所述，過高的的或過低的修正因子選取，都會相應(yīng)的產(chǎn)生圖像重構(gòu)問題，因此綜合考慮，在當(dāng)前模擬設(shè)置下，選取-0.5作為固定的修正因子，此數(shù)值與文獻[8]中提供的修正因子計算公式的計算結(jié)果接近。圖6為修正因子為-0.5時，迭代次數(shù)隨源項所在極徑的變化情況，不同形狀的數(shù)據(jù)點為不同的源項所在極角大小。趨勢再次證明，在確定的修正因子條件下，迭代次數(shù)與源項極角關(guān)系很小，且隨著極徑的增加所需迭代次數(shù)也相應(yīng)增加。

利用3次多項式對源項所在極徑R與最佳迭代次數(shù)N進行擬合，擬合結(jié)果如圖6所示，關(guān)系式為：

圖6 最佳迭代次數(shù)N與源項所在極徑的關(guān)系Fig.6 The relationship between the optimal number of iterations N and the polar diameter of the source

N=[-1.40+0.15R-1.1×10-3R2+

2.77×10-6R3]

(5)

式中R為源項位置的極徑大小，需要注意的是由于迭代次數(shù)不能是小數(shù)，因此需要將迭代次數(shù)的計算結(jié)果四舍五入取整。并且由于編碼板以及探測器是都是正方形，并且編碼板形狀是以正方形對角線為對稱軸的軸對稱圖形，因此通過模擬0～45°的源項位置可以代表大多數(shù)的源項位置情況。利用此關(guān)系式，即可完成針對此設(shè)置下的不同成像位置的迭代參數(shù)選取。

3.3 優(yōu)化方法驗證與對比

為驗證這一方法的準(zhǔn)確性，在不改變其他實驗?zāi)M設(shè)置的條件下，模擬(40 mm, 40 mm)與(200 mm, 200 mm)位置的伽馬點源的成像結(jié)果，其位置所在極徑分別為56.57、282.84 mm，若β取值為-0.5，通過式 (5)計算得到的理想迭代次數(shù)為4次與16次，計算如圖7所示，圖7(a)與(b)為源項位置(40 mm, 40 mm)時，未進行迭代計算與迭代4次之后的成像結(jié)果，圖7(c)為對應(yīng)利用傳統(tǒng)δ解碼[1]成像結(jié)果，圖7(d)與(e)為源項位置(200 mm, 200 mm)時，未進行迭代計算與迭代16次之后的成像結(jié)果，圖7(f)為對應(yīng)利用傳統(tǒng)δ解碼[11-12]成像結(jié)果。評價參數(shù)如表1所示。

表1 成像結(jié)果評價參數(shù)PCNR對比表Table 1 Comparison table of evaluation parameter PCNR of imaging results

圖7 利用計算得到的迭代次數(shù)所重構(gòu)的圖像結(jié)果Fig.7 The result of image reconstruction using the number of iterations calculated

由圖7可知，經(jīng)過計算公式所得次數(shù)進行迭代后的位置分辨率達到1.2 mm，并且如表 1所示，雖然未經(jīng)迭代的成像結(jié)果PCNR參數(shù)均小于傳統(tǒng)δ解碼方法，但經(jīng)過特定次數(shù)的迭代之后，成像結(jié)果的PCNR參數(shù)明顯提高，成像質(zhì)量相應(yīng)有效增加。結(jié)果說明利用此方法計算得到的迭代次數(shù)與修正因子的大小，針對本次模擬條件，進行基于MLEM互補算法成像可以有效提高成像質(zhì)量。

4 結(jié)論

1)在進行編碼近場成像時，可以利用MLEM互補算法消除偽像。本文通過對大量不同修正參數(shù)的成像結(jié)果的圖像性能參數(shù)進行比較，針對MLEM互補算法的修正參數(shù)提出選取方法，為FOV內(nèi)任意位置點源的圖像重構(gòu)提供最優(yōu)的參數(shù)選擇。與未經(jīng)迭代的圖像結(jié)果以及傳統(tǒng)δ解碼的重構(gòu)結(jié)果進行對比，可有效提高成像質(zhì)量。

2)結(jié)果顯示MLEM互補算法的最優(yōu)成像參數(shù)與源項所在位置的極徑關(guān)系密切，而與所在角度關(guān)聯(lián)性不大。造成這一現(xiàn)象的原因是：偽像主要是由源項在探測器上投影的不均勻造成的。而源項所在極徑不同造成的偽像差異性大于角度不同形成的差異性，因此相比于極徑，角度對參數(shù)選取的影響較小。