張 斌，王 英，艾憲蕓，魏 星

（防化研究院 國(guó)民核生化災(zāi)害防護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102205）

編碼孔γ相機(jī)作為一種新型的核輻射監(jiān)測(cè)設(shè)備，能通過生成可視化圖像使操作人員快速識(shí)別并定位視野中放射源目標(biāo)，可極大縮短作業(yè)時(shí)間，減少人員受到的照射劑量，在環(huán)境輻射監(jiān)測(cè)、核廢料處理、通關(guān)口岸貨物放射性檢測(cè)等領(lǐng)域具有廣泛應(yīng)用。與普通光學(xué)相機(jī)不同的是，編碼孔γ相機(jī)并不能直接對(duì)放射性物質(zhì)進(jìn)行成像，而要在成像過程中進(jìn)行所謂的“編碼”和“解碼”：放射源發(fā)射出的γ射線經(jīng)過編碼孔準(zhǔn)直器調(diào)制后在探測(cè)器上形成投影圖像，經(jīng)過離散采樣后傳入計(jì)算機(jī)中，再通過特定圖像重建算法恢復(fù)出可視化的放射性強(qiáng)度二維分布圖像。線性算法和迭代算法是編碼孔γ相機(jī)圖像重建中的兩類主要算法。線性算法如交叉相關(guān)法［1］等具有算法簡(jiǎn)單、計(jì)算速度快的優(yōu)點(diǎn)，但在投影圖像受噪聲影響較大的情況下重建圖像質(zhì)量較差。以MLEM 算法為代表的迭代算法能很好地抑制噪聲，對(duì)不完全數(shù)據(jù)亦能重建，然而MLEM 算法的缺點(diǎn)是收斂速度慢，且在投影圖像中包含噪聲較大時(shí)，迭代次數(shù)過高會(huì)使重建圖像中噪聲隨迭代次數(shù)的增加而加劇。針對(duì)這種情況，本文將直接解調(diào)法（direct demodulation method）［2］思想應(yīng)用于MLEM 算法中，在迭代過程中加入先驗(yàn)物理約束條件，并以交叉相關(guān)法計(jì)算結(jié)果作為迭代初始值，加快MLEM迭代的收斂速度，提高重建圖像的質(zhì)量。

1 算法介紹

1.1 二維MLEM 算法

極大似然估計(jì)法是根據(jù)極大似然準(zhǔn)則來求解似然函數(shù)中未知參數(shù)的一種參數(shù)估計(jì)法，解決這類問題的基本方法是求偏導(dǎo)解方程，但在圖像重建中，未知參數(shù)眾多，計(jì)算如此繁復(fù)的方程組幾乎是不可能的。用期望最大化法實(shí)現(xiàn)極大似然圖像重建最早由Shepp等［3］于1982年提出，最終形成MLEM 算法。期望最大化法以迭代法為基礎(chǔ)，通過迭代更新方程組中參數(shù)的估計(jì)值，使似然函數(shù)不斷增大，直至逼近其最大值，最終獲得極大似然參數(shù)的估計(jì)值，這樣就把原本極其復(fù)雜的極大似然估計(jì)問題轉(zhuǎn)化為簡(jiǎn)單的求解期望值最大化過程。MLEM 算法是一種基于泊松統(tǒng)計(jì)模型的貝葉斯后驗(yàn)概率最大的圖像復(fù)原方法，該算法以符合實(shí)際情況的模型為基礎(chǔ)，重建圖像的分辨率和本底噪聲抑制能力較傳統(tǒng)線性算法提升很多，因此在核醫(yī)學(xué)等領(lǐng)域得到了廣泛應(yīng)用。編碼孔γ相機(jī)圖像重建中的二 維MLEM 算 法［4］可 表 示 為：

式中：“＊”表示卷積運(yùn)算；“?”表示相關(guān)運(yùn)算；fk（x，y）為第k次迭代后的圖像估計(jì)值；p（x，y）為γ射線透過編碼板后在探測(cè)器上的實(shí)際投影圖像；h（x，y）為編碼函數(shù)。

1.2 基于約束的MLEM 算法（DD-MLEM）

在理想無噪聲的情況下，MLEM 迭代的解是收斂的，重建出的圖像與目標(biāo)圖像最終趨于一致。但實(shí)際成像中所測(cè)得的投影圖像數(shù)據(jù)中必然包含有各種噪聲，這使得迭代不能完全收斂：在迭代次數(shù)較低時(shí)，圖像質(zhì)量逐漸提高；當(dāng)?shù)螖?shù)超過某一值后，圖像質(zhì)量反而變差。這主要是因?yàn)镸LEM 迭代的修正準(zhǔn)則是要求投影估計(jì)值與實(shí)測(cè)投影越接近越好，但對(duì)于重建出的圖像卻無任何約束。為改善MLEM 迭代在實(shí)際使用中的收斂性，本文借鑒了直接解調(diào)法思想來優(yōu)化MLEM算法，提高重建圖像質(zhì)量。直接解調(diào)法最初是一種應(yīng)用于高能射線天文望遠(yuǎn)鏡中的成像方法，該方法與其他算法的顯著區(qū)別是利用先驗(yàn)的物理知識(shí)作為約束條件，并添加到迭代求解的過程中。在迭代過程中加入的物理約束是對(duì)于被測(cè)目標(biāo)的已知信息，所以使用這種方法不僅能有效抑制投影數(shù)據(jù)中本底噪聲的影響，使迭代收斂性得到明顯改善，而且最終重建出的圖像較單純使用MLEM 算法能更加趨近于真實(shí)情況。

實(shí)際成像的圖像重建過程中主要使用本底約束條件：若fk（x，y）＜fb（x，y），則令fk（x，y）＝fb（x，y），其中fb（x，y）為本底約束值。如果本底噪聲水平較低，可選擇fb（x，y）＝0；如果本底噪聲較強(qiáng)，為了提高成像靈敏度，需先得到背景輻射強(qiáng)度分布db，db可取估計(jì)值或由實(shí)際觀測(cè)數(shù)據(jù)中剔除放射源影響后獲取，然后再求解調(diào)制方程Hfb＝db（H 為系統(tǒng)的傳遞函數(shù)），最終得到迭代的本底約束fb。

此外，在進(jìn)行圖像重建迭代時(shí)，一般選擇均勻的灰度圖像作為迭代的初始值，本文仿真中采用了交叉相關(guān)法的重建結(jié)果作為迭代的初始值，實(shí)驗(yàn)證明這種方法可在一定程度上加快迭代收斂速度。交叉相關(guān)法的計(jì)算過程可由下式表示：

其中：O′（x，y）為對(duì)放射源目標(biāo)O（x，y）的估計(jì)圖像；P（x，y）為投影圖像；G（x，y）為解碼函數(shù)（由成像系統(tǒng)的編碼函數(shù)決定）；N（x，y）為噪聲項(xiàng)。

由此，本文中基于約束的MLEM 算法計(jì)算過程可分為以下幾個(gè)步驟：

1）令k＝0，初始化fk（x，y）（一般將初始化圖像f0（x，y）設(shè)置為各像素相等的均勻灰度圖像，本文中將其設(shè)置為交叉相關(guān)法的重建圖像）；

2）應(yīng)用式（1）計(jì)算fk＋1（x，y）；

3）對(duì)fk＋1（x，y）施加一定的物理約束條件；

4）根據(jù)一定的準(zhǔn)則判斷迭代重建結(jié)果是否達(dá)到收斂要求，若已達(dá)到要求，則停止迭代；若未達(dá)到要求，則令k＝k＋1，將更新后的fk（x，y）再次轉(zhuǎn)入執(zhí)行步驟2。

2 仿真實(shí)驗(yàn)與分析

2.1 仿真條件

仿真實(shí)驗(yàn)以19×19的改良均勻冗余陣列（modified uniformly redundant array，MURA）編碼孔準(zhǔn)直器為基礎(chǔ)，成像系統(tǒng)參數(shù)選擇為：物距500mm，像距80mm，編碼板小孔單元邊長(zhǎng)2mm，編碼板厚度5mm，材料為鎢。探測(cè)目標(biāo)為位于視野中心軸線上的1 個(gè)點(diǎn)源，活度為1 MBq，能量為662keV。由上述條件模擬的理想點(diǎn)源圖像和投影圖像如圖1所示。迭代時(shí)在模擬投影圖像中加入隨機(jī)高斯噪聲，得到含噪聲的投影數(shù)據(jù)進(jìn)行圖像重建。

2.2 評(píng)價(jià)準(zhǔn)則

在對(duì)比兩種算法重建圖像的效果時(shí)，有時(shí)通過主觀觀察很難判斷，為了驗(yàn)證所提出方法的有效性，采用歸一化均方誤差（normalized root mean square error，NRMSE）衡量算法的收斂速度，用相關(guān)系數(shù)（correlation coefficient，CC）評(píng)價(jià)重建圖像的質(zhì)量，它們的定義［5］分別為：

式中：xij和分別為真實(shí)圖像和重建圖像的像素灰度值和分別為真實(shí)圖像和重建圖像的平均像素灰度值。歸一化均方誤差表示重建圖像與真實(shí)圖像的接近程度，其值越小代表重建效果越好，而相關(guān)系數(shù)評(píng)價(jià)的是重建圖像與真實(shí)圖像的相似度，其值越大表示重建圖像與真實(shí)圖像越接近。

圖1 模擬的點(diǎn)源圖像和投影圖像Fig.1 Simulated point source image and projection image

2.3 仿真結(jié)果分析

在實(shí)際成像中，當(dāng)成像系統(tǒng)計(jì)數(shù)率足夠高時(shí)，投影圖像中疊加的噪聲趨近于高斯分布。仿真實(shí)驗(yàn)采用了這一噪聲模型，在DD-MLEM迭代中構(gòu)造基于高斯分布的本底并進(jìn)行圖像重建，將重建后的背景圖像作為迭代約束條件，采用交叉相關(guān)法計(jì)算結(jié)果作為迭代初值。圖2和圖3分別為MLEM 算法和DD-MLEM 算法迭代次數(shù)n＝30、90、270時(shí)的重建圖像。通過直接觀察可看出，在選取合適迭代次數(shù)的情況下，兩種算法均能得到相對(duì)較好的重建效果，但在迭代次數(shù)過高時(shí)，MLEM 算法重建的圖像局部區(qū)域噪聲被放大，原本光滑的圖像出現(xiàn)離散的斑點(diǎn)，在視覺上難以接受；而DD-MLEM 迭代則不然，在過高迭代次數(shù)下依然能很好地控制圖像噪點(diǎn)，保持圖像的可辨識(shí)性。

圖2 MLEM 迭代重建圖像Fig.2 Reconstructed image by MLEM iteration

圖3 DD-MLEM 迭代重建圖像Fig.3 Reconstructed image by DD-MLEM iteration

為進(jìn)一步評(píng)價(jià)重建圖像的質(zhì)量，利用式（3）和式（4）給出的評(píng)價(jià)準(zhǔn)則對(duì)兩種算法的重建圖像效果進(jìn)行定量分析。圖4示出兩種算法所重建圖像的歸一化均方誤差隨迭代次數(shù)的變化情況。由圖4a可看出，在迭代次數(shù)相對(duì)較低時(shí)，兩種算法的歸一化均方誤差均隨迭代次數(shù)的增加而減小，重建圖像質(zhì)量不斷提高，但DD-MLEM算法的收斂速度顯然要快得多；當(dāng)?shù)螖?shù)繼續(xù)增加時(shí)，由圖4b可知，MLEM 算法的歸一化均方誤差在達(dá)到一最小值后迅速增大，重建圖像的質(zhì)量也隨之變差，而DD-MLEM 算法在迭代次數(shù)過高的情況下歸一化均方誤差只是略微變大，重建圖像的質(zhì)量仍能保持得比較好。圖5示出兩種算法的相關(guān)系數(shù)隨迭代次數(shù)的變化情況，可更加明顯地看出DD-MLEM 算法收斂速度不但更快，而且提高了算法最佳迭代次數(shù)，重建圖像的質(zhì)量也優(yōu)于原算法。在迭代次數(shù)過高時(shí)，MLEM迭代的重建圖像質(zhì)量迅速下降，而DD-MLEM 迭代則能將圖像質(zhì)量一直保持在較高水平上。

圖4 歸一化均方誤差隨迭代次數(shù)的變化Fig.4 NRMSE of two algorithms versus iteration time

圖5 相關(guān)系數(shù)隨迭代次數(shù)的變化Fig.5 CC of two algorithms versus iteration time

3 結(jié)論

在編碼孔γ相機(jī)圖像重建中，MLEM 迭代算法能較為有效地抑制統(tǒng)計(jì)噪聲，較傳統(tǒng)交叉相關(guān)法恢復(fù)出的圖像質(zhì)量要高得多，但缺點(diǎn)是收斂速度較慢且迭代次數(shù)過高時(shí)圖像質(zhì)量變差。本文提出一種基于約束的MLEM 算法，在迭代過程中加入物理約束條件，并用交叉相關(guān)法計(jì)算值作為迭代初始值。實(shí)驗(yàn)結(jié)果證明，DD-MLEM 算法不僅收斂速度加快，且重建圖像質(zhì)量和迭代的收斂性與收斂速度均得到顯著提高。使用這種方法可更加充分地利用被觀測(cè)目標(biāo)的已知信息，從而有效抑制由于噪聲等因素所引起的迭代解的病態(tài)特征，相比僅用迭代法或相關(guān)法能更好地重建出反映真實(shí)情況的圖像。值得注意的是，此方法重建圖像的質(zhì)量對(duì)約束條件有較強(qiáng)的依賴，對(duì)于復(fù)雜的測(cè)量環(huán)境，應(yīng)根據(jù)現(xiàn)場(chǎng)實(shí)際情況選取更為合理的物理約束條件以保證成像質(zhì)量。

［1］ 何會(huì)林，董永偉，吳伯冰，等.空間硬X 射線編碼孔成像望遠(yuǎn)鏡樣機(jī)的研制及初步實(shí)驗(yàn)結(jié)果［J］.高能物理與核物理，2007，31（5）：437-441.HE Huilin，DONG Yongwei，WU Bobing，et al.A prototype of spatial hard X-ray coded aperture imaging telescope and the primary laboratory experimental results［J］.High Energy Physics and Nuclear Physics，2007，31（5）：437-441（in Chinese）.

［2］ 李惕碚，吳枚.高能天文中成像和解譜的直接方法［J］.天體物理學(xué)報(bào)，1993，13（3）：215-224.LI Tibei，WU Mei.The direct method for spectral and image restoration in high energy astronomy［J］.Acta Astrophysica Sinica，1993，13（3）：215-224（in Chinese）.

［3］ SHEPP L A，VARDI Y.Maximum likelihood restoration for emission tomography［J］.IEEE Transactions on Medical Imaging，1982，1（2）：113-122.

［4］ 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.

［5］ 何佳偉，劉東升，桂志國(guó).可變有序子集PML 算法在PET 中的應(yīng)用［J］.中北大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，31（6）：646-650.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（in Chinese）.