基于改進Richardson-Lucy算法的中子照相圖像幾何不銳度修正方法研究

伍華杰 甘佺 呂永強 汪啟宏 伍一凡 宋婧 鳳麟核團隊

1（中國科學院合肥物質科學研究院 合肥 230031）

2（中國科學技術大學 合肥 230026）

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

4（中子科學國際研究院 青島 266199）

5（山東省中子科學技術重點實驗室 青島 266199）

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

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

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

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

1 中子照相圖像復原方法

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

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

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

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

1.1 中子照相圖像的PSF估計

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

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

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

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

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

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

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

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

1.2 改進的RL算法

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

式中：It表示第t次迭代的結果；k表示PSF；kT表示k在水平方向的翻轉。

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

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

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

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

改進的RL 算法在中子照相圖像反卷積前對GSN進行去除，抑制因GSN產生的振鈴效應。

2 實驗結果和分析

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

圖4 線對樣品實物圖Fig.4 Photograph of the line-pair sample

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

圖5 線對樣品的中子照相圖像、復原結果與對應的行灰度曲線圖(a) 原始中子照相圖像，(b) RL算法，(c) 小波降噪RL，(d) 調制核正則化RL，(e) 文獻[12]方法，(f) 本文方法，(g) 原始中子照相圖像對應行灰度曲線圖，(h) RL算法對應行灰度曲線圖，(i) 小波降噪RL對應行灰度曲線圖，(j) 調制核正則化RL對應行灰度曲線圖，(k) 文獻[12]方法對應行灰度曲線圖，(l) 本文方法對應行灰度曲線圖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

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

表1 文獻[12]方法和改進的RL算法復原圖像后的平均梯度和空間頻率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的表達式如下：

一般來說，AG和SF值越大，表明圖像的前景紋理層次越多、信息結構越復雜。圖5中可看出，原始中子照相圖像、經RL算法、小波降噪RL算法和調制核正則化RL 算法復原后的中子照相圖像仍然存在大量的GSN，而當圖像中存在較多的GSN 時，指標會異常增大，影響AG 和SF 對圖像前景信息復原效果的判定，因此，只計算文獻［12］方法和本文方法復原后圖像的AG和SF指標，如表1所示。

從表1可以看出，與文獻［12］方法相比，經過本文提出的算法復原后，線對中子照相圖像AG 和SF值分別提高了60.23%和29.90%，圖像銳度有明顯的提升。

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

3 結語

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

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