魏 星 王 鋒 王 英 肖無云 艾憲蕓 張 斌 徐紅鵑 張 磊 馬新華

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

γ射線成像技術(shù)引入核輻射監(jiān)測領(lǐng)域后，在發(fā)現(xiàn)熱點、確定沾染區(qū)域、定位放射源、輻射防護(hù)評估等方面得到廣泛應(yīng)用，具有巨大的技術(shù)優(yōu)勢。以非接觸方式在一定距離外獲取 γ射線源的精細(xì)圖像，這是其他核輻射監(jiān)測技術(shù)無法實現(xiàn)的。為滿足當(dāng)前核輻射監(jiān)測任務(wù)對高性能γ射線成像系統(tǒng)的需求，提出了研制一種探測效率高、能量響應(yīng)范圍寬、能區(qū)分不同放射性核素圖像的γ射線成像系統(tǒng)。

γ射線無法聚焦，因此人們發(fā)明了各種非聚焦技術(shù)來實現(xiàn)對γ射線源的成像。其中基于空間編碼的編碼孔(Coded Aperture, CA)成像技術(shù)[1]因其高靈敏度、高分辨率的優(yōu)點，成為目前核輻射監(jiān)測中應(yīng)用最多的成像技術(shù)。為滿足 Nyquist采樣定律，對CA成像系統(tǒng)探測器位置分辨率的基本要求是不能超過小孔尺寸的一半，因此一個CA系統(tǒng)通常需要大量讀出通道以達(dá)到足夠的位置計算精度。在像距一定的情況下，CA系統(tǒng)的角度分辨率由探測器位置分辨率決定。因此，保持角度分辨率的同時擴(kuò)展探測器面積，必然導(dǎo)致讀出通道數(shù)急劇增加。這意味著整個系統(tǒng)的復(fù)雜度和成本大大提高，對所有通道的一致性校正等工作也變得非常困難。CA系統(tǒng)的另一個缺點是探測器厚度受限制，無法提高對高能射線的探測效率。這是因為高能射線在探測器中的平均自由程更大，如果探測器較厚，斜入射光子就容易跨過像素邊界發(fā)生作用，造成像素間串?dāng)_，引起圖像模糊。為克服CA系統(tǒng)存在的這些問題，提出研制基于旋轉(zhuǎn)調(diào)制器(Rotating Modulator, RM)成像技術(shù)[2]的γ射線成像系統(tǒng)。這是因為RM成像技術(shù)是一種時間調(diào)制成像技術(shù)，在時間域超采樣，可以重建超過系統(tǒng)幾何角度分辨率的圖像。因此，RM系統(tǒng)可選尺寸更大的探測器，達(dá)到和CA系統(tǒng)相當(dāng)?shù)慕嵌确直媛省Ｓ钟捎赗M系統(tǒng)采用獨立的非位置靈敏探測器，不存在像素間串?dāng)_，因此可以增加探測器厚度以提高對高能射線的探測效率。

1 RM技術(shù)簡介

RM 成像系統(tǒng)構(gòu)成主要包括一個柵條狀旋轉(zhuǎn)準(zhǔn)直器和一個非位置靈敏探測器陣列，其角度分辨率由柵條寬度和像距決定[3]。準(zhǔn)直器旋轉(zhuǎn)過程中，入射γ射線會在探測器平面形成柵條陰影。因為柵條寬度等于探測器直徑，因此陰影對探測器入射面的遮擋比例在0%?100%變化，又因為準(zhǔn)直器具有180°對稱性，因此變化周期為p。這使得探測器輸出呈現(xiàn)以p為周期的起伏變化，這一過程即時間調(diào)制過程。假設(shè)目標(biāo)空間劃分為N個像素點，S(n)表示目標(biāo)空間第n個像素的強(qiáng)度；系統(tǒng)包含D個探測器，d表示探測器序號；采樣時間周期為T，單位s，整個測量過程持續(xù)M個時間周期，Od(m)表示第m個時間周期內(nèi)探測器d的計數(shù)。假設(shè)系統(tǒng)傳遞函數(shù)為P，則目標(biāo)空間第n個像素在第m個時間周期向探測器d的傳遞系數(shù)為Pd(n,m)。為簡化公式書寫，約定式(2)中d的取值范圍是1?D，n和n￠的取值范圍是1?N，m的取值范圍是1?M。由此得到調(diào)制過程的離散數(shù)學(xué)表達(dá)式寫為：

式中，Bd(m)是探測器d在第m個時間周期的本底計數(shù)。系統(tǒng)參數(shù)已知的情況下，可以預(yù)先計算Pd(n,m)[4]。因為準(zhǔn)直器旋轉(zhuǎn)的周期性，因此將時間軸映射到準(zhǔn)直器轉(zhuǎn)角坐標(biāo)來計算其歸一化值。

圖像重建的過程就是根據(jù)式(1)求解 S(n)的過程。求解的第一步是通過互相關(guān)將時間域的測量數(shù)據(jù)Od(m)轉(zhuǎn)換到圖像空間：

式(3)是典型的投影方程，A可看做大小為N×N的傳遞矩陣，A的計算由式(3)給出。

互相關(guān)圖像C中，在真實源圖像周圍出現(xiàn)波紋狀起伏邊緣，因此需要進(jìn)一步求解式(2)來消除這些邊緣。迭代法適合于求解式(2)。受統(tǒng)計噪聲影響，代數(shù)迭代算法容易出現(xiàn)噪聲放大，在重建圖像中產(chǎn)生偽像，因此提出用基于泊松統(tǒng)計模型的極大似然-期望最大化(Maximum Likelihood Expectation Maximization, MLEM)迭代算法[5]重建圖像。MLEM算法迭代公式為：

根據(jù)探測器計數(shù)，可以反推出目標(biāo)空間一點的源強(qiáng)。假設(shè)測量過程中本底計數(shù)率恒為b，min?1；沒有準(zhǔn)直器情況下，位于像素n的源在探測器d處引起的計數(shù)率恒為s，min?1。設(shè)：

根據(jù)式(6)求出計數(shù)率 s，假設(shè)被測目標(biāo)距離較遠(yuǎn)，則源強(qiáng)S(n)可由式(7)計算：

式中，l1是物距；ε是探測效率；θ是源所處位置的極角；r是探測器入射面半徑；φ是對應(yīng)能量光子的產(chǎn)額。

2 樣機(jī)設(shè)計

根據(jù)實際條件，RM樣機(jī)研制選用了7個濱松的CH292 NaI(Tl)探測器。該型號探測器靈敏區(qū)尺寸為Φ40 mm×40 mm。7個探測器排列方式可參考文獻(xiàn)[3]。1號探測器位于旋轉(zhuǎn)軸與像面交點，即探測器陣列的中心；2?7號探測器間隔60°均勻排列于1號探測器四周，其與 1號探測器的圓心距離均為7.1cm。這種排列可以最大化利用整個視場范圍，并使系統(tǒng)對目標(biāo)空間的靈敏度更加均勻。準(zhǔn)直器由8條寬4 cm、厚1 cm的鉛條組成。鉛條平行排列，間隔4 cm。每個鉛條由2 mm鋁包裹。準(zhǔn)直器旋轉(zhuǎn)中心位于第4條和第5條柵條之間，計算出準(zhǔn)直器旋轉(zhuǎn)半徑為30 cm。由此確定鉛條長度為60 cm。所有鉛條固定于兩個鋁制圓環(huán)（內(nèi)環(huán)）之間，兩個鋁制圓環(huán)通過鋼珠方式與另兩個鋁制圓環(huán)（外環(huán)）耦合，外環(huán)再與鋁合金框架固定連接。這樣，內(nèi)環(huán)就可帶動鉛條轉(zhuǎn)動，內(nèi)環(huán)的轉(zhuǎn)動由一個步進(jìn)電機(jī)驅(qū)動。準(zhǔn)直器平面和探測器平面間距離80 cm，樣機(jī)的幾何角度分辨率約2.9°。

樣機(jī)的電子學(xué)和軟件部分采用自主研發(fā)方式。設(shè)計了7通道讀出電路。各通道可循環(huán)讀出設(shè)定時間采樣周期內(nèi)的計數(shù)值。時間采樣周期在1?108ms可調(diào)。為驗證該成像系統(tǒng)具有能量分辨能力，能夠獲取不同放射性核素的圖像，采用了自主研發(fā)的4 k道數(shù)字化多道分析器(Digital Multi-Channel Analyzer, DMCA)。此外，電路部分還集成了步進(jìn)電機(jī)控制卡和驅(qū)動器。讀出電路通過以太網(wǎng)口與計算機(jī)通信，電機(jī)控制器則通過PCI口與計算機(jī)通信，相應(yīng)地開發(fā)了專用的計算機(jī)數(shù)據(jù)采集處理軟件和電機(jī)控制軟件。為避免長期轉(zhuǎn)動過程中，可能出現(xiàn)的皮帶輪打滑以及電機(jī)轉(zhuǎn)速不均勻，數(shù)據(jù)采集電路設(shè)計了霍爾開關(guān)接口。在準(zhǔn)直器外環(huán)上安裝了一個固定的霍爾開關(guān)，在內(nèi)環(huán)上安裝了一個磁鐵。這樣，內(nèi)環(huán)每旋轉(zhuǎn)一周，磁鐵會觸發(fā)霍爾開關(guān)一次，據(jù)此對內(nèi)環(huán)旋轉(zhuǎn)位置和時間的關(guān)系進(jìn)行校正。最后，在樣機(jī)前方安裝了可見光攝像頭，用于捕捉可見光圖像。樣機(jī)實物如圖1所示。

圖1 樣機(jī)實物圖Fig.1 Photograph of the RM prototype.

3 實驗

3.1 能量分辨率

首先對探測器能量分辨率進(jìn)行測量。在沒有準(zhǔn)直器情況下，測量了活度為6.29×104Bq的137Cs源和活度為 5.8×104Bq的60Co的混合能譜，其中 1號探測器測得的能譜如圖 2，兩種核素的特征峰非常明顯。計算得到探測器能量分辨率為7%@662keV。

圖2 實測能譜Fig.2 Measured spectrum.

3.2 探測效率

為了保證7路輸出數(shù)據(jù)的一致性，需對7個探測器的效率進(jìn)行一致性校正。這里定義絕對探測效率為：

式中，I是記錄的粒子數(shù)，Iin是入射粒子數(shù)。效率測量采用一枚活度為1.85×108Bq的137Cs密封源。其結(jié)構(gòu)符合 GB/T 13366-2009規(guī)定，活性區(qū)尺寸Φ4mm× 4mm，外殼尺寸Φ8 mm×10 mm。

首先測量了7路本底輸出，再在沒有準(zhǔn)直器的情況下對置于6 m外的源進(jìn)行測量。其中本底測量時間2 h，源測量時間30 min。測得的計數(shù)率和探測效率如表 1所示，7個探測器平均探測效率為22.2%。

表1 探測器計數(shù)及效率Table 1 Detectors count rate and efficiency.

3.3 靈敏度

RM系統(tǒng)的歸一化靈敏度定義為：

式中，η是超分辨因子，對RM系統(tǒng)，η=1。式中的取平均值表示在時間域上的平均。根據(jù)式(9)計算可知：在約±14.6°范圍內(nèi)，平均靈敏度約 0.31，超過這個范圍后，受調(diào)制的探測器數(shù)量減少，靈敏度下降，直至沒有探測器受調(diào)制，靈敏度下降為 0。因此該樣機(jī)有效視場范圍約±14.6°。

3.4 成像實驗

成像實驗和效率測量使用同一放射源。該樣機(jī)對置于6 m外的放射源進(jìn)行成像測量，放射源位于極坐標(biāo) θ=3°、方位角 ψ=60°。準(zhǔn)直器轉(zhuǎn)速0.86r·min?1，時間采樣周期 200 ms，共測量 300 s，圖3(a)是2號探測器隨時間變化的輸出曲線。從探測器輸出數(shù)據(jù)中截取了整圈的部分，映射到準(zhǔn)直器轉(zhuǎn)角坐標(biāo)后，有效測量時間為256 s，輸出曲線如圖3(b)。圖4是用MLEM算法重建的圖像。

為驗證輻射圖像中源的位置是否準(zhǔn)確，將輻射圖像與可見光圖像用加權(quán)平均法[6]進(jìn)行融合。圖 5是將源置于不同位置用同一配準(zhǔn)矩陣得到的兩幅融合圖像，該配準(zhǔn)矩陣通過基于外部控制點的方法[7]得到。圖5中輻射圖像能夠很好地與實際位置匹配，證明了重建源位置的準(zhǔn)確性。

對重建圖像實際能達(dá)到的角度分辨率進(jìn)行了測試。測試采用兩枚活度為1.85×108Bq的137Cs密封源，源的結(jié)構(gòu)尺寸與效率測量所用的源相同。在6 m距離處，將兩個源逐漸靠近，得到相應(yīng)的重建圖像。以峰的半高寬處為分辨標(biāo)準(zhǔn)，當(dāng)兩個源間隔 2°時（距離20.9 cm），兩個峰的半高寬處相互分開，當(dāng)兩個源間隔1.8°時（間隔18.9 cm），兩個峰的半高寬處相互重疊。圖6是重建得到的兩個點源的峰。

圖3 2號探測器隨時間輸出曲線(a)和隨轉(zhuǎn)角變化曲線(b)Fig.3 Plot of detector 2 output (counts) against time (a) and rotating angle (b).

圖6 兩個點源的重建峰間隔：(a) 2°，(b) 1.8°Fig.6 Reconstructed peaks for two point sources.Separation in angle: (a) 2°, (b) 1.8°

3.5 源強(qiáng)計算

將放射源置于不同距離，測量 5次，反推得到的源強(qiáng)見表2。

圖4 MLEM算法的點源探測重建圖像Fig.4 Reconstructed image of a point source by MLEM.

表2 源強(qiáng)計算Table 2 Calculated source intensity.

圖5 融合圖像Fig.5 Superimposed images.

4 結(jié)語

RM成像技術(shù)是一種時間調(diào)制的γ射線成像技術(shù)。與采用CA技術(shù)的成像系統(tǒng)比較，RM系統(tǒng)可以采用更厚的探測器而不存在像素間串?dāng)_問題，因此探測效率更高，有更寬的能量響應(yīng)范圍。此外，RM 系統(tǒng)能夠同時完成能譜測量和成像測量，具有分辨不同放射性核素強(qiáng)度分布圖像的能力。通過超分辨重建圖像，RM系統(tǒng)能夠達(dá)到和CA系統(tǒng)相當(dāng)?shù)慕嵌确直媛省?/p>

研制的RM成像樣機(jī)完成了對RM成像技術(shù)的驗證。樣機(jī)采用鉛制旋轉(zhuǎn)準(zhǔn)直器和NaI(Tl)探測器陣列，可同時測量能譜和放射性核素強(qiáng)度分布圖像，能量分辨率達(dá)到7%@662 keV。用MLEM算法重建圖像能分辨間隔2°的兩個點源。作為驗證樣機(jī)，其性能仍有很大的提高空間，例如用LaBr3(Ce)探測器以提高探測效率和能量分辨率。用鎢合金材料以減少準(zhǔn)直器厚度，減少穿透射線成份和圖像重建過程中引入物理約束條件等。

1 Jean in 't Zand. Coded aperture camera imaging concept[OL]. http://astrophysics.gsfc.nasa.gov/cai/coded_intr.html, 1996

2 Durouchoux P, Hudson H, Hurford G, et al. Gamma-ray imaging with a rotating modulator[J]. Astronomy and Astrophysics, 1983, 120(1): 150?155

3 魏星, 王鋒, 肖無云, 等. 一種旋轉(zhuǎn)調(diào)制器成像技術(shù)研究[J]. 核技術(shù), 2014, 37(2): 020401 WEI Xing, WANG Feng, XIAO Wuyu, et al. Study on rotating modulator imaging technique[J]. Nuclear Techniques, 2014, 37(2): 020401

4 Dadurkevicius V, Ralys D A. Modulation patterns of astronomical imaging system based on rotating grids[J].Astrophysics and Space Science, 1985, 113: 233?247

5 Shepp L A, Vardi Y. Maximum likelihood reconstruction for emission tomography[J]. IEEE Transaction on Medical Imaging, 1982, MI-1: 133?122

6 郭佳. 紅外與可見光圖像融合技術(shù)的研究[D]. 西安:西安工業(yè)大學(xué), 2009 GUO Jia. Research on fusion technology of infrared and visible image[D]. Xi'an: Xi'an Industrial University, 2009

7 陳顯毅. 圖像配準(zhǔn)技術(shù)及其MATLAB編程實現(xiàn)[M]. 北京: 電子工業(yè)出版社, 2009 CHEN Xianyi. Image registration technology and its MATLAB programming[M]. Beijing: Publishing House of Electronics Industry, 2009