伍華杰 甘佺 呂永強(qiáng) 汪啟宏 伍一凡 宋婧 鳳麟核團(tuán)隊(duì)

1（中國(guó)科學(xué)院合肥物質(zhì)科學(xué)研究院 合肥 230031）

2（中國(guó)科學(xué)技術(shù)大學(xué) 合肥 230026）

3（中科超睿技術(shù)有限公司 青島 266041）

4（中子科學(xué)國(guó)際研究院 青島 266199）

5（山東省中子科學(xué)技術(shù)重點(diǎn)實(shí)驗(yàn)室 青島 266199）

中子照相技術(shù)作為一種對(duì)物體內(nèi)部結(jié)構(gòu)進(jìn)行無(wú)損檢測(cè)的手段，因其具有穿透能力強(qiáng)［1］、輕元素靈敏、成分識(shí)別精準(zhǔn)［2］及抗干擾強(qiáng)等獨(dú)特的優(yōu)勢(shì)，在航天航空［3-4］、鋰電池［5-6］、混凝土建筑和非開(kāi)箱安檢［7］等方面檢測(cè)中都具有不可替代的重要價(jià)值。隨著中子照相技術(shù)的發(fā)展，包括基于反應(yīng)堆、大型加速器、小型加速器和微型中子源等的中子照相系統(tǒng)陸續(xù)投入研發(fā)和市場(chǎng)，因此中子照相無(wú)損檢測(cè)具有良好的發(fā)展前景。雖然中子照相系統(tǒng)逐漸向微小型化發(fā)展，但由于基本的成像結(jié)構(gòu)未變，中子照相圖像依然存在幾何不銳度，亟須開(kāi)展關(guān)于幾何不銳度修正的方法研究。

中子照相圖像幾何不銳度的修正，屬于中子照相圖像的復(fù)原問(wèn)題，且一直是一個(gè)難點(diǎn)。金煒等［8］提出一種基于小波降噪的Richardson-Lucy（wavelet-RL）算法，假設(shè)中子照相圖像中的少量噪聲為高斯白噪聲，采用BivaShrink雙變量去噪模型，刻畫(huà)小波變換后圖像相同位置處不同尺度子帶系數(shù)間的相關(guān)性，顯著提升圖像的去噪性能，從而提高圖像幾何不銳度的修正效果。Souza 等［9-10］提出搖擺曲線（Rocking Curve，RC）的概念，RC 半峰全寬（Full Width at Half Maximum，F(xiàn)WHM）表示幾何不銳度，并采用Richardson-Lucy（RL）算法進(jìn)行復(fù)原，圖像的前景邊緣部分變得更加銳利，但是未考慮γ白斑噪聲（Gamma Spots Noise，GSN）的去除。喬雙等［11］將調(diào)制核引入到正則化RL方法中，提出一種用于中子照相圖像去模糊的方法SK-TV-RL，能夠抑制中子照相圖像弱噪聲的放大，并能復(fù)原圖像中的細(xì)節(jié)信息。Yazid 等［12］對(duì)RTP（PUSPATI TRIGA Mark II）反應(yīng)堆得到的中子照相圖像，采用中值濾波與RL結(jié)合的方法進(jìn)行復(fù)原，然而直接采用中值濾波去除GSN會(huì)損壞圖像前景的全局細(xì)節(jié)與邊緣信息，降低圖像的整體復(fù)原效果。

目前，有關(guān)幾何不銳度修正的研究大多集中于含有弱噪聲的中子照相圖像，對(duì)于含有強(qiáng)噪聲的中子照相圖像的研究較少，而在采用RL方法修正幾何不銳度的過(guò)程中，中子照相圖像中的GSN難以去除并會(huì)不斷放大，嚴(yán)重影響圖像的質(zhì)量。因此，本文提出一種基于改進(jìn)RL 的中子照相圖像幾何不銳度修正算法，并與以往有關(guān)中子照相圖像的幾何不銳度修正算法作對(duì)比，驗(yàn)證采用改進(jìn)的RL算法進(jìn)行幾何不銳度修正的效果。

本文采取的方法主要分為兩個(gè)過(guò)程：第一步是根據(jù)中子照相圖像的幾何不銳度，對(duì)用于表示幾何不銳度的點(diǎn)擴(kuò)散函數(shù)（Point Spread Function，PSF）建立數(shù)學(xué)模型；第二步是結(jié)合第一步得出的PSF，采用改進(jìn)的RL 方法對(duì)圖像中的前景邊緣細(xì)節(jié)進(jìn)行復(fù)原，以達(dá)到更好的幾何不銳度修正效果，提升圖像的分辨率。

1 中子照相圖像復(fù)原方法

實(shí)際上，中子照相系統(tǒng)只能獲取到降質(zhì)后的中子照相圖像，可以大致將中子照相圖像的降質(zhì)過(guò)程描述為一個(gè)離散線性系統(tǒng)［13］：

式中：Idegraded表示降質(zhì)中子照相圖像；I表示原始清晰中子照相圖像；k表示PSF；n表示加性噪聲；*為卷積。中子照相圖像的復(fù)原需要同時(shí)考慮PSF和加性噪聲，既要準(zhǔn)確估計(jì)PSF，也應(yīng)盡可能地去除加性噪聲的干擾。

本文提出的方法中，首先根據(jù)中子源出口到像探測(cè)器的距離、中子源出口直徑和被檢品中心到像探測(cè)器的距離，精準(zhǔn)地估計(jì)出中子照相圖像的PSF，然后根據(jù)中子照相圖像中GSN的特點(diǎn)，提出一種改進(jìn)的RL（improved-RL）算法。算法總體流程如圖1所示。

圖1 算法總體流程圖Fig.1 Block diagram of the improved-RL algorithm

1.1 中子照相圖像的PSF估計(jì)

一般地，中子照相圖像的PSF 估計(jì)分為三個(gè)過(guò)程：第一步是計(jì)算中子照相圖像的幾何不銳度；第二步是為PSF 選擇合適的分布模型；第三步是根據(jù)幾何不銳度和分布模型估計(jì)出PSF。

通常情況下，在中子照相系統(tǒng)中，通過(guò)中子源出口的中子束穿過(guò)被檢品某點(diǎn)后，在像探測(cè)器處會(huì)產(chǎn)生一圓形光斑，如圖2 所示。其中，E表示被檢品邊緣上的點(diǎn)，C表示被檢品的中心點(diǎn)。該圓形光斑被稱為幾何不銳度，其直徑大小一般為［14］：

圖2 中子照相系統(tǒng)結(jié)構(gòu)示意圖Fig.2 Structural diagram of the neutron imaging system

式中：L表示中子源出口到像探測(cè)器距離；D表示中子源出口直徑；d為被檢品中心到像探測(cè)器的距離。

與X 射線衍射場(chǎng)相似，中子源各點(diǎn)處發(fā)射的中子束強(qiáng)度分布形狀為RC［15］。同理，所有穿過(guò)被檢品邊緣某一點(diǎn)后的中子束強(qiáng)度分布形狀也為RC。將RC 的FWHM 表示為中子照相圖像的幾何不銳度G［8］。如圖3所示，其中，W表示RC的FWHM。用于表示PSF 的函數(shù)模型包括高斯分布［16］、柯西分布［12，17］等。本文采用二維柯西分布模型來(lái)表示PSF，表達(dá)式為：

圖3 搖擺曲線示意圖Fig.3 Sketch of the rocking curve

式中：Cauchy(x，y)表示二維柯西分布函數(shù)；ξ表示為分散度參數(shù)。幾何不銳度G與參數(shù)ξ的關(guān)系為：

G確定后，根據(jù)式（4）可得參數(shù)ξ的值，代入式（3）中可得表示中子照相圖像的PSF。

1.2 改進(jìn)的RL算法

RL 算法是Richardson 和Lucy 提出的一種非盲反卷積圖像復(fù)原方法，算法表達(dá)式如式（5）所示：

式中：It表示第t次迭代的結(jié)果；k表示PSF；kT表示k在水平方向的翻轉(zhuǎn)。

由于照相得到的中子照相圖像中存在灰度值較高的GSN，若直接采用RL算法進(jìn)行復(fù)原，得到的復(fù)原圖像會(huì)出現(xiàn)嚴(yán)重的振鈴效應(yīng)。因此，建立如下數(shù)學(xué)模型：

式中：γ(·)為被GSN 污染的區(qū)域，γ?(·)為未被GSN 污染的區(qū)域，在該數(shù)學(xué)模型下需要檢測(cè)出所有的γ(·)。由于GSN具有孤立性，可采用拉普拉斯算子檢測(cè)出γ(·)，并配合中值濾波去除γ(·)內(nèi)的GSN，具體方法如式（7）所示：

式中：M(·)表示中值濾波；U(x，y)表示點(diǎn)(x，y)的鄰域；?2Idegraded(x，y)表示Idegraded二階梯度之和；T為指定閾值。

上述過(guò)程只會(huì)對(duì)含有GSN 的區(qū)域進(jìn)行去噪處理，即使在多次使用的情況下，仍然不會(huì)損壞圖像中的高頻信息。因此，將以上步驟引入到RL 算法中，以對(duì)圖像復(fù)原過(guò)程中的振鈴效應(yīng)進(jìn)行抑制。改進(jìn)的RL算法的表達(dá)式為：

改進(jìn)的RL 算法在中子照相圖像反卷積前對(duì)GSN進(jìn)行去除，抑制因GSN產(chǎn)生的振鈴效應(yīng)。

2 實(shí)驗(yàn)結(jié)果和分析

本文中子照相系統(tǒng)的中子源類型為熱中子（能量區(qū)間為0.005～0.5 eV），成像端采用LiF 轉(zhuǎn)換屏和電荷耦合元件（Charge Coupled Device，CCD）結(jié)合的數(shù)字成像方式。其中，該CCD 相機(jī)為Atik Infinity Camera，分辨率為1 392×1 040。為降低CCD相機(jī)的暗電流，將CCD相機(jī)放置于密封的真空箱內(nèi)，溫度降至-10 ℃。另外，在真空箱外包裹一層鉛屏蔽層以防止中子束或次級(jí)粒子的輻射。本文選取線對(duì)樣品進(jìn)行實(shí)驗(yàn)［18-19］，線對(duì)樣品實(shí)物如圖4 所示。線對(duì)樣品擁有多組粗細(xì)不同的線狀區(qū)域，有利于檢測(cè)中子照相系統(tǒng)的分辨率。

圖4 線對(duì)樣品實(shí)物圖Fig.4 Photograph of the line-pair sample

首先，測(cè)量中子照相系統(tǒng)的參數(shù)值，得到中子源出口到像探測(cè)器距離L與中子源出口直徑D的比值，即L/D=80，被檢品中心到像探測(cè)器的距離d=2 cm，轉(zhuǎn)換屏大小為20 cm×20 cm。拍照方式設(shè)置為單次曝光，曝光300 s后記錄并保存圖像。由于保存后的圖像尺寸較大，且線對(duì)樣品的特征信息主要集中在圖像的中心區(qū)域，需要進(jìn)行裁剪，裁剪后的圖像尺寸為512 pixel×512 pixel；然后根據(jù)式（2）計(jì)算出中子照相圖像的幾何不銳度，再根據(jù)式（3）、（4）計(jì)算出對(duì)應(yīng)的PSF；最后對(duì)線對(duì)中子照相圖像分別采用RL 算法、小波降噪RL算法、調(diào)制核正則化RL算法、文獻(xiàn)［12］方法和本文改進(jìn)的RL算法進(jìn)行圖像復(fù)原，并繪制復(fù)原后圖像的行灰度曲線，實(shí)驗(yàn)結(jié)果如圖5 所示。從實(shí)驗(yàn)結(jié)果可看出，相對(duì)于RL 算法、小波降噪RL 算法和調(diào)制核正則化RL 算法，經(jīng)本文改進(jìn)的RL算法處理后，復(fù)原后圖像噪聲的波動(dòng)幅度明顯減?。幌鄬?duì)于文獻(xiàn)［12］方法，復(fù)原后圖像行灰度曲線的峰谷值更大、波峰更窄，幾何不銳度修正程度有較大的提升。

圖5 線對(duì)樣品的中子照相圖像、復(fù)原結(jié)果與對(duì)應(yīng)的行灰度曲線圖(a) 原始中子照相圖像，(b) RL算法，(c) 小波降噪RL，(d) 調(diào)制核正則化RL，(e) 文獻(xiàn)[12]方法，(f) 本文方法，(g) 原始中子照相圖像對(duì)應(yīng)行灰度曲線圖，(h) RL算法對(duì)應(yīng)行灰度曲線圖，(i) 小波降噪RL對(duì)應(yīng)行灰度曲線圖，(j) 調(diào)制核正則化RL對(duì)應(yīng)行灰度曲線圖，(k) 文獻(xiàn)[12]方法對(duì)應(yīng)行灰度曲線圖，(l) 本文方法對(duì)應(yīng)行灰度曲線圖Fig.5 Neutron photographic image of line-pair sample, restored results and corresponding line grey profile(a) Original neutron photographic image, (b) RL algorithm, (c) Wavelet-RL, (d) SK-TV-RL, (e) Method in Ref.[12], (f) Our method,(g) Corresponding line grey profile of original neutron photographic image, (h) Corresponding line grey profile of RL algorithm, (i)Corresponding line grey profile of wavelet-RL, (j) Corresponding line grey profile of SK-TV-RL, (k) Corresponding line grey profile of method in Ref.[12], (l) Corresponding line grey profile of our method

為客觀衡量幾何不銳度修正效果，通過(guò)無(wú)參考圖像質(zhì)量評(píng)價(jià)指標(biāo)平均梯度（Average Gradient，AG）和空間頻率（Spatial Frequency，SF）對(duì)中子照相圖像復(fù)原質(zhì)量進(jìn)行分析，如表1 所示。其中，AG 可以反映圖像中前景信息的紋理變化特征，SF能夠反映圖像的整體活躍程度，均可衡量圖像的清晰程度。

表1 文獻(xiàn)[12]方法和改進(jìn)的RL算法復(fù)原圖像后的平均梯度和空間頻率Table 1 AGs and SFs of restored images based on the algorithm presented in Ref.[12] and improved-RL algorithm

AG的定義如下：

式中：M和N分別表示圖像的寬度和高度；D為梯度算子。

SF的表達(dá)式如下：

一般來(lái)說(shuō)，AG和SF值越大，表明圖像的前景紋理層次越多、信息結(jié)構(gòu)越復(fù)雜。圖5中可看出，原始中子照相圖像、經(jīng)RL算法、小波降噪RL算法和調(diào)制核正則化RL 算法復(fù)原后的中子照相圖像仍然存在大量的GSN，而當(dāng)圖像中存在較多的GSN 時(shí)，指標(biāo)會(huì)異常增大，影響AG 和SF 對(duì)圖像前景信息復(fù)原效果的判定，因此，只計(jì)算文獻(xiàn)［12］方法和本文方法復(fù)原后圖像的AG和SF指標(biāo)，如表1所示。

從表1可以看出，與文獻(xiàn)［12］方法相比，經(jīng)過(guò)本文提出的算法復(fù)原后，線對(duì)中子照相圖像AG 和SF值分別提高了60.23%和29.90%，圖像銳度有明顯的提升。

因此，通過(guò)結(jié)合圖5 和表1 的結(jié)果可看出，本文中改進(jìn)的RL 算法相比于RL 算法、小波降噪RL 算法和調(diào)制核正則化RL 算法，能抑制GSN 在中子照相圖像復(fù)原過(guò)程中產(chǎn)生的振鈴效應(yīng)，更有利于圖像前景細(xì)節(jié)信息的復(fù)原；相比于文獻(xiàn)［12］方法，有更好的幾何不銳度修正效果。

3 結(jié)語(yǔ)

本文提出了一種基于改進(jìn)RL 的中子照相圖像幾何不銳度修正算法對(duì)幾何不銳度進(jìn)行修正。通過(guò)準(zhǔn)確計(jì)算表示幾何不銳度的PSF，結(jié)合改進(jìn)的RL算法，修正了因幾何不銳度引起的圖像降質(zhì)，同時(shí)避免了GSN對(duì)圖像復(fù)原過(guò)程的影響。實(shí)驗(yàn)結(jié)果表明，本文提出的方法可以獲得高質(zhì)量的復(fù)原圖像，能夠進(jìn)行高質(zhì)量的中子照相無(wú)損檢測(cè)，對(duì)提升中子照相系統(tǒng)的檢測(cè)效率與精度意義重大。

作者貢獻(xiàn)聲明伍華杰負(fù)責(zé)圖像數(shù)據(jù)的獲取、算法的設(shè)計(jì)、程序的編寫、文章的編寫與修訂；甘佺負(fù)責(zé)實(shí)驗(yàn)裝置的搭建、文章的編寫與修訂；呂永強(qiáng)提供部分算法思路、文章的修訂；汪啟宏負(fù)責(zé)提供部分算法思路；伍一凡負(fù)責(zé)驗(yàn)證算法的正確性；宋婧負(fù)責(zé)提供實(shí)驗(yàn)裝置的指導(dǎo)。