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

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

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

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

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

1 基于互補(bǔ)編碼的MLEM成像算法

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

(1)

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

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

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

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)可知，互補(bǔ)編碼的MLEM成像算法是利用反編碼修正后的編碼函數(shù)，對(duì)迭代結(jié)果fk(x，y)進(jìn)行投影后與原始投影相比較，并進(jìn)行反投影運(yùn)算。在“投影”與“反投影”過程中均利用反編碼函數(shù)對(duì)噪聲進(jìn)行抑制，達(dá)到有效減少噪聲的目的。

2 編碼孔成像過程模擬

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

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

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

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

將伽馬源所在位置設(shè)置為坐標(biāo)原點(diǎn)，伽馬點(diǎn)源、編碼板、探測(cè)器沿Z軸依次放置，源項(xiàng)所在位置即探測(cè)系統(tǒng)的FOV為坐標(biāo)系的XOY平面，原點(diǎn)為探測(cè)器與編碼板的中心投影位置。伽馬源距編碼板中心80 cm，編碼板中心與探測(cè)器前端面的距離為20 cm。投影成像的放大倍數(shù)為1.25倍，此編碼板探測(cè)系統(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平面，改變伽馬源在成像平面的位置，將源項(xiàng)分別設(shè)置于極角θ為0°、15°、30°、45°處，距原點(diǎn)距離R即源項(xiàng)所在極徑分別設(shè)置為50、80、160、200、240和290 mm處分別進(jìn)行成像模擬。使用基于互補(bǔ)編碼板的MLEM方法進(jìn)行源項(xiàng)位置重構(gòu)，測(cè)量不同源項(xiàng)位置合適迭代次數(shù)N與修正因子β的關(guān)系。

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

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

(4)

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

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

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

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

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

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

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

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

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

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

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

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

3.3 優(yōu)化方法驗(yàn)證與對(duì)比

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

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

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

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

4 結(jié)論

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

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