肖灑 何偉波 帥茂兵 蘭明聰 韋孟伏

?

基于天然μ子散射特征評(píng)估鈾板內(nèi)水平狹縫寬度的方法

肖灑 何偉波 帥茂兵 蘭明聰 韋孟伏

（中國(guó)工程物理研究院材料研究所 江油 621907）

針對(duì)經(jīng)典μ子成像方法用于特殊對(duì)象結(jié)構(gòu)探測(cè)的技術(shù)難點(diǎn)，提出了一種通過(guò)直接比較對(duì)象發(fā)生結(jié)構(gòu)變化前后μ子散射特征的差異來(lái)進(jìn)行形變類型和尺度判定的新思路，稱為正向參比法。相比于經(jīng)典的逆向反演法，該方法在原理上對(duì)于初始結(jié)構(gòu)已知的特殊核對(duì)象的結(jié)構(gòu)變化探測(cè)更有優(yōu)勢(shì)。從初步模擬結(jié)果來(lái)看，在合理的天然μ子通量下，該方法可實(shí)現(xiàn)鈾板內(nèi)亞毫米尺度水平狹縫的存在及其尺度的準(zhǔn)確判定。

μ子，重核材料，多重庫(kù)侖散射，非參數(shù)檢驗(yàn)，接受者操作特征分析

μ子斷層成像技術(shù)憑借其無(wú)需額外放射源、安全系數(shù)高、穿透能力強(qiáng)、對(duì)重核材料靈敏等眾多優(yōu)勢(shì)成為了輻射成像領(lǐng)域一個(gè)新興的研究熱點(diǎn)[1?3]。在利用輻射成像手段進(jìn)行目標(biāo)對(duì)象結(jié)構(gòu)診斷相關(guān)領(lǐng)域中，常需要研究重核材料（例如鈾）及其元件經(jīng)過(guò)一段時(shí)間放置或特定工藝處理后的內(nèi)部結(jié)構(gòu)變化問(wèn)題，并且目標(biāo)對(duì)象的初始結(jié)構(gòu)組成通常是可以準(zhǔn)確預(yù)知的。對(duì)于此類問(wèn)題，傳統(tǒng)的μ子散射成像反演法的基本思路是：分別對(duì)變化前、后兩個(gè)目標(biāo)對(duì)象進(jìn)行μ子散射特征測(cè)量并利用特定算法反演出其內(nèi)部三維結(jié)構(gòu)圖像，通過(guò)對(duì)比二者圖像的差異即可獲知目標(biāo)對(duì)象之間的結(jié)構(gòu)差異。這是一種直接、自然的方法，對(duì)于某些類型的目標(biāo)對(duì)象也是行之有效的，但該方法往往涉及不確定性大、復(fù)雜性高的圖像反演過(guò)程。如果目標(biāo)對(duì)象組成結(jié)構(gòu)較為復(fù)雜、尤其是具有大量孔洞、狹縫或低密度內(nèi)填充物等各種類空心結(jié)構(gòu)時(shí)，傳統(tǒng)的成像反演方法極有可能受到嚴(yán)重的偽影干擾而難以給出清晰、準(zhǔn)確的圖像，因而使其應(yīng)用范圍大為受限。基于上述原因，亟需發(fā)展一種適用范圍更廣、實(shí)現(xiàn)難度更低的基于μ子散射特征來(lái)識(shí)別目標(biāo)對(duì)象之間結(jié)構(gòu)差異的新方法。

在這種背景下，我們提出了一種基于目標(biāo)對(duì)象發(fā)生結(jié)構(gòu)變化前后對(duì)天然宇宙射線μ子的散射特征變化來(lái)識(shí)別目標(biāo)對(duì)象所發(fā)生的結(jié)構(gòu)變化類型以及評(píng)估該特定結(jié)構(gòu)變化的尺度的新方法。由于該方法不涉及逆向的圖像反演（重建）過(guò)程，而是采用了與預(yù)先建立的參比數(shù)據(jù)庫(kù)進(jìn)行比對(duì)的方法。并且，由于該方法規(guī)避了繁瑣且易被干擾的圖像重建過(guò)程，與目標(biāo)對(duì)象的初始結(jié)構(gòu)復(fù)雜程度無(wú)關(guān)，因而對(duì)于重核材料對(duì)象結(jié)構(gòu)變化診斷相關(guān)領(lǐng)域而言，其適用性要比傳統(tǒng)成像反演方法更為廣泛。

1 方法原理

為了解決傳統(tǒng)的μ子散射成像反演法用于特殊目標(biāo)對(duì)象結(jié)構(gòu)差異診斷上的難題，我們嘗試了一種新的方法，如圖1所示。其基本思想是：首先，建立目標(biāo)對(duì)象標(biāo)準(zhǔn)狀態(tài)下的理論模型(P1)，然后充分考慮目標(biāo)對(duì)象所有可能發(fā)生的一系列結(jié)構(gòu)變化（包括形變類型、尺度、位置等多個(gè)參量），建立一組或多組與之分別對(duì)應(yīng)的微擾模型(P2)。在這一步驟中，微擾模型體系必須具有很好的完備性。其次，通過(guò)蒙特卡羅模擬獲得上述各個(gè)微擾模型對(duì)應(yīng)的μ子散射特征量(P3)，然后將這些散射特征量與相應(yīng)的結(jié)構(gòu)變化形成一一對(duì)應(yīng)的映射關(guān)系，從而組成一套完備的參比數(shù)據(jù)庫(kù)待用(P4)。P1?P4這幾個(gè)步驟（圖1中虛線框內(nèi)）構(gòu)成了該方法的先決條件。而在日常應(yīng)用（步驟1?3）中，對(duì)于實(shí)際待測(cè)目標(biāo)對(duì)象，采用μ子成像平臺(tái)實(shí)驗(yàn)測(cè)量方法獲得其μ子散射特征量，然后將該散射特征與數(shù)據(jù)庫(kù)中的數(shù)據(jù)進(jìn)行比對(duì)，尋找出相似度最高的條目，即可在給定的顯著性水平下，判定待測(cè)目標(biāo)對(duì)象發(fā)生了相應(yīng)的結(jié)構(gòu)變化。由于該方法摒棄了較為繁瑣、易引入誤差的圖像反演過(guò)程，試圖通過(guò)直接比較兩種狀態(tài)下的目標(biāo)對(duì)象對(duì)μ子的散射特征的差異來(lái)獲取目標(biāo)對(duì)象的結(jié)構(gòu)變化信息，因而僅涉及正向的散射測(cè)量過(guò)程，故稱之為“正向參比法”。

圖1 正向參比法的實(shí)現(xiàn)及其所涉及的關(guān)鍵技術(shù)

2 正向參比法的實(shí)現(xiàn)

以水平放置的鈾板的GEANT4模擬結(jié)果為例，對(duì)正向參比法用于鈾板內(nèi)部水平狹縫寬度識(shí)別的實(shí)現(xiàn)以及其中的關(guān)鍵技術(shù)作詳細(xì)說(shuō)明。

2.1 μ子散射特征量參比數(shù)據(jù)庫(kù)的建立

2.1.1 微擾模型的建立

本文采用的鈾板狹縫模型如圖2所示，左、右兩邊分別為理想鈾板模型和基于該模型擴(kuò)展得到的一組分別帶有寬度為d（1≤≤，為指定的有限正整數(shù)）的狹縫的鈾板模型。為簡(jiǎn)單起見(jiàn)，假定鈾板只可能產(chǎn)生水平方向的狹縫，并且所有可能的狹縫寬度的集合={d}為有限集，至少是可以近似看作有限集。由于微擾模型必須具有完備性，因此需要針對(duì)集合中的每一個(gè)元素d分別設(shè)計(jì)相應(yīng)的狹縫模型，然后進(jìn)行后續(xù)的μ子散射蒙特卡羅模擬。

圖2 鈾板狹縫模型

2.1.2 μ子散射的蒙特卡羅模擬

為了獲得足夠多的μ子散射徑跡，我們編制了一套基于GEANT4的μ子成像系統(tǒng)來(lái)模擬μ子在目標(biāo)對(duì)象中的多重庫(kù)侖散射(Multiple Coulomb Scattering, MCS)行為。天然宇宙射線μ子源項(xiàng)采用CRY軟件包[4]來(lái)生成。模擬計(jì)算中考慮了μ子與材料的庫(kù)侖散射、電離以及軔致輻射等主要物理作用。

2.1.3 數(shù)據(jù)庫(kù)的構(gòu)建

根據(jù)Moliere經(jīng)驗(yàn)公式，μ子散射極角概率分布函數(shù)()對(duì)不同目標(biāo)對(duì)象具有較好的特異性，故本文選取()作為μ子散射特征量。然而，由于天然μ子能譜范圍極廣、入射方向各異，()往往難以采用精確的解析函數(shù)進(jìn)行描述。并且由于μ子與物質(zhì)多重庫(kù)侖散射的隨機(jī)性，即使對(duì)于同一目標(biāo)，兩次實(shí)驗(yàn)得到的散射角分布也可能存在差異?？紤]這兩方面的因素，實(shí)際應(yīng)用中，需要進(jìn)行同一條件下的多次重復(fù)實(shí)驗(yàn)來(lái)構(gòu)建總體()的多個(gè)獨(dú)立樣本。

2.2 目標(biāo)對(duì)象結(jié)構(gòu)變化的判定——基于假設(shè)檢驗(yàn)的二元分類器

根據(jù)上述分析，利用正向參比法評(píng)估目標(biāo)物體結(jié)構(gòu)變化的核心問(wèn)題可歸結(jié)為待測(cè)目標(biāo)對(duì)象μ子散射特征量與預(yù)建數(shù)據(jù)庫(kù)中μ子散射特征量之間的比對(duì)問(wèn)題?？紤]到實(shí)際采用的是多次重復(fù)實(shí)驗(yàn)獲得的多個(gè)相互獨(dú)立的散射角分布樣本，因此問(wèn)題的實(shí)質(zhì)是：在總體分布形式()未知的條件下，檢驗(yàn)兩個(gè)獨(dú)立樣本是否取自同一總體的非參數(shù)檢驗(yàn)問(wèn)題。這里采用比較累積頻數(shù)分布函數(shù)的Kolmogorov- Smirnov (KS)檢驗(yàn)方法[5]，其具體過(guò)程如下：

令為單個(gè)μ子的散射角度量：

= (Δθ)2+ (Δθ)2(1)

式中：Δθ、Δθ分別為和平面上的累積偏轉(zhuǎn)角。假定兩次獨(dú)立實(shí)驗(yàn)的目標(biāo)物體對(duì)應(yīng)的天然μ子散射角總體分別服從1()和2()分布，則兩次實(shí)驗(yàn)測(cè)得的散射角可以看作分別來(lái)自這兩個(gè)總體的樣本。考慮假設(shè)檢驗(yàn)問(wèn)題：

在給定的顯著性水平下，計(jì)算KS檢驗(yàn)的值（零假設(shè)0可被拒絕的最小顯著性水平）。如果≤，則拒絕0，認(rèn)為兩次實(shí)驗(yàn)?zāi)繕?biāo)對(duì)象存在差異；反之，如果>，則接受0，認(rèn)為兩次實(shí)驗(yàn)?zāi)繕?biāo)對(duì)象完全相同。

基于上述KS檢驗(yàn)對(duì)0拒絕/接受的結(jié)論，可以構(gòu)建一個(gè)判定目標(biāo)對(duì)象是否發(fā)生變化的二元分類器。令顯著性水平在一定范圍內(nèi)變化，即可在接受者操作特征(Receiver Operating Characteristic, ROC)空間內(nèi)繪制出相應(yīng)的ROC曲線。ROC曲線可以對(duì)該分類器的性能進(jìn)行可視化，給出分類器的真陽(yáng)性率(True Positive Rate, TPR)（收益）與假陽(yáng)性率(False Positive Rate, FPR)（開(kāi)銷）之間的相對(duì)平衡關(guān)系[6?7]。下文中，我們會(huì)通過(guò)ROC曲線來(lái)評(píng)估該基于KS檢驗(yàn)的二元分類器性能。

3 結(jié)果與討論

3.1 不同μ子通量下可識(shí)別的最小狹縫寬度

由于所用μ子通量不同會(huì)導(dǎo)致散射樣本容量不同，進(jìn)而可能影響假設(shè)檢驗(yàn)的判定結(jié)果，因此有必要對(duì)上述方法在各種不同μ子通量條件下對(duì)目標(biāo)對(duì)象是否存在差異的識(shí)別能力進(jìn)行研究。

這里以厚度1 cm的鈾質(zhì)均勻平板作為研究對(duì)象，在平板中心水平面上分別設(shè)置寬度為1= 0.1mm、2=0.2 mm、3=0.5 mm和4=1.0 mm的4種水平狹縫。以40組理想鈾板（狹縫寬度0=0）的μ子散射模擬結(jié)果獲取的40個(gè)散射角分布樣本，以及40組帶有寬度d（為[1, 4]范圍內(nèi)某一指定正整數(shù)）狹縫的鈾板的μ子散射模擬結(jié)果獲取的40個(gè)散射角分布樣本共同組成參比數(shù)據(jù)庫(kù)，以帶有寬度d狹縫的鈾板作為待測(cè)目標(biāo)對(duì)象，依次選取μ子通量為50 cm?2、125 cm?2、250 cm?2和1000 cm?2，分別采用上述正向參比法考察在一定的顯著性水平下，各種μ子通量條件下該方法對(duì)狹縫的識(shí)別能力。對(duì)于給定的，一對(duì)(TPR, FPR)就對(duì)應(yīng)了ROC空間的一個(gè)點(diǎn)。在(0, 1)范圍內(nèi)改變，就可以得到一系列ROC空間點(diǎn)組成的一條ROC曲線，如圖3?5所示。根據(jù)ROC曲線的含義可知，曲線越靠近左上角，即曲線下面積(Area Under Curve, AUC)越接近1，相應(yīng)的分類器性能也越好。由圖3可知，50cm?2的μ子通量下該分類器已可實(shí)現(xiàn)對(duì)寬度低至1.0 mm的水平狹縫的完美分辨（AUC約為1），相應(yīng)的閾值≈ 0。

圖3 μ子通量為50 cm?2時(shí)，分類器對(duì)于帶有寬度di的水平狹縫的鈾板的識(shí)別性能ROC曲線 (a) d1=0.1 mm，AUC 0.60，(b) d2=0.2 mm，AUC 0.75， (c) d3=0.5 mm，AUC 0.96，(d) d4=1.0 mm，AUC 1.00

隨著進(jìn)一步增大入射μ子通量，每個(gè)散射樣本的容量也會(huì)增大，從而對(duì)判定結(jié)果也可能產(chǎn)生影響。圖4中結(jié)果顯示，μ子通量增大到125cm?2時(shí)，選擇適當(dāng)?shù)拈撝?，該分類器可近似完美識(shí)別最小寬度為0.5mm的狹縫（AUC約為1）。

圖4 μ子通量為125 cm?2時(shí)，分類器對(duì)于帶有寬度di的水平狹縫的鈾板的識(shí)別性能ROC曲線 (a) d1=0.1 mm，AUC 0.55，(b) d2=0.2 mm，AUC 0.86， (c) d3=0.5 mm，AUC 1.00

圖5表明，當(dāng)μ子通量增大到1000 cm?2時(shí)，選擇適當(dāng)?shù)拈撝?，該分類器可較完美識(shí)別最小寬度為0.2 mm的狹縫（AUC約為0.99）。

圖5 μ子通量為1000 cm?2時(shí)，分類器對(duì)于帶有寬度d的水平狹縫的鈾板的識(shí)別性能ROC曲線 (a) d1=0.1 mm，AUC 0.83，(b) d2=0.2 mm，AUC 0.99

3.2 多種狹縫寬度可能性情況下的識(shí)別性能評(píng)估

根據(jù)§3.1的結(jié)果，在μ子通量1000cm?2的條件下，采用正向參比法可完美分辨出1cm鈾板內(nèi)寬度最小可達(dá)0.2mm的狹縫。然而，在實(shí)際應(yīng)用中，除了需要對(duì)目標(biāo)對(duì)象是否產(chǎn)生狹縫進(jìn)行簡(jiǎn)單定性的Y/N判斷以外，目標(biāo)對(duì)象的結(jié)構(gòu)變化可能性比較多，即使對(duì)于同一類型和位置的變化（如中心層的水平狹縫），往往同時(shí)存在多種尺度（如狹縫寬度）的可能性。為了評(píng)估這種情況下正向參比法對(duì)目標(biāo)對(duì)象結(jié)構(gòu)變化的識(shí)別性能，構(gòu)建如下場(chǎng)景：設(shè)計(jì)狹縫寬度分別為234的三個(gè)鈾板作為待測(cè)目標(biāo)對(duì)象，以狹縫寬度分別為234的鈾板各40組（共計(jì)120組）的模擬結(jié)果獲取的120個(gè)散射角分布樣本組成參比數(shù)據(jù)庫(kù)，選取μ子通量為1000 cm?2，考察該情形下正向參比法對(duì)于不同狹縫寬度的識(shí)別能力。其結(jié)果如圖6所示。

圖6 μ子通量1000 cm?2條件下，鈾板目標(biāo)對(duì)象狹縫寬度同時(shí)存在圖中三種可能時(shí)，正向參比法對(duì)狹縫寬度di的目標(biāo)對(duì)象的識(shí)別性能ROC曲線 (a) d2=0.2 mm，AUC 0.99，(b) d3= 0.5 mm，AUC 1.00， (c) d4=1.0 mm，AUC 1.00

由圖6可知，在μ子通量達(dá)到1000 cm?2的條件下，對(duì)于任意一種狹縫寬度為d(=2, 3, 4)的鈾板對(duì)象來(lái)說(shuō)，其ROC曲線的AUC值（分別為0.99、1和1）都接近1的理想值。由此可見(jiàn)，在該鈾板對(duì)象結(jié)構(gòu)變化的可能性被限定為上述三種寬度的水平狹縫的前提下，正向參比法可以較為理想地實(shí)現(xiàn)其中任意一種寬度的狹縫與另外兩種狹縫的區(qū)分與識(shí)別。這表明，在適當(dāng)?shù)南薅l件下，正向參比法可以實(shí)現(xiàn)目標(biāo)對(duì)象特定結(jié)構(gòu)變化尺度的準(zhǔn)確評(píng)估。

4 結(jié)語(yǔ)

相對(duì)于經(jīng)典μ子成像方法用于特殊對(duì)象結(jié)構(gòu)探測(cè)的技術(shù)難點(diǎn)，本文提出的正向參比法更適用于初始結(jié)構(gòu)已知的特殊對(duì)象內(nèi)部結(jié)構(gòu)變化探測(cè)，有望實(shí)現(xiàn)更好的空間分辨率和時(shí)間效率。從帶有0.1?1.0mm范圍內(nèi)不同寬度水平狹縫的1 cm厚度鈾板模型的GEANT4模擬結(jié)果來(lái)看，在合理的天然μ子通量條件下，該方法對(duì)于目標(biāo)對(duì)象的微小結(jié)構(gòu)變化尺度進(jìn)行準(zhǔn)確判定和評(píng)估，未來(lái)有望在復(fù)雜大型核對(duì)象的無(wú)損探測(cè)等相關(guān)領(lǐng)域取得顯著成效。

1 Borozdin K, Hogan G, Morris C,. Surveillance: radiographic imaging with cosmic-ray muons[J]. Nature, 2003, 422(6929): 277. DOI: 10.1038/422277a.

2 Schultz L, Borozdin K, Gomez J,. Image reconstruction and materialdiscriminationcosmic ray muon radiography[J]. Nuclear Instruments and Methods in Physics ResearchA, 2004, 519(3): 687?694. DOI: 10.1016/j.nima.2003.11.035.

3 Morris C, Alexander C, Bacon J,. Tomographic imaging with cosmic ray muons[J]. Science & Global Security, 2008, 16(1?2): 37?53. DOI:10.1080/ 08929880802335758.

4 Hagmann C, Lange D, Wright D. Monte Carlo simulation of proton-induced cosmic-ray cascades in the atmosphere (UCRL-TM-229452)[R]. USA: Lawrence Livermore National Laboratory, 2007.

5 以恒冠. 面向特殊對(duì)象的宇宙線繆子散射成像方法研究[D]. 北京: 清華大學(xué), 2015. YI Hengguan. Research on cosmic ray muon scattering tomography methods for special objects inspection[D]. Beijing: Tsinghua University, 2015.

6 Fawcett T. An introduction to ROC analysis[J]. Pattern Recognition Letters, 2006, 27: 861?874. DOI:10.1016/j. patrec.2005.10.010.

7 鄒洪俠, 秦鋒, 程澤凱, 等. 二類分類器的ROC曲線生成算法[J]. 計(jì)算機(jī)技術(shù)與發(fā)展, 2009, 19(6): 109?112. ZOU Hongxia, QIN Feng, CHENG Zekai,. Algorithm for generating ROC curve of two-classifier[J]. Computer Technology & Development, 2009, 19(6): 109?112.

Width evaluation of horizontal slit in uranium slabs using a forward comparison method based on cosmic-ray muon scattering

XIAO Sa HE Weibo SHUAI Maobing LAN Mingcong WEI Mengfu

(Institute of Materials, China Academy of Engineering Physics, Jiangyou 621907, China)

Background: Muon tomography based on multiple Coulomb scattering (MCS) has been widely used in structural detection of nuclear objects. However, for some particular objects, the tomographic methodwould be thwarted by failing to imaging them. Purpose:In view of the technical difficulties in structural detection of those objects, a novel method named “Forward Comparison Method (FCM)” was proposed. Methods:Both the deformation and its scale were identified effectively by directly evaluating the evolution of muon scattering properties with respect to the object. This method, instead of the traditional inverse one, showed a methodological superiority in detection of nuclear objects with prior knowledge of the original structure.Results:Preliminary results based on a computational simulation demonstrated that FCM could be used to identify the horizontal slit in a uranium slab and evaluate its width on the magnitude of sub-millimeters, on condition that the muon flux was reasonable.Conclusion:The FCM method can be effectively applied in structural detection of particular nuclear objects of hollow structure on a satisfactory timescale.

Muon, High-Z material, MCS, Non-parametric test, Receiver operating characteristic analysis

TP312，TL81

10.11889/j.0253-3219.2017.hjs.40.050202

肖灑，男，1983年出生，2011年畢業(yè)于中國(guó)工程物理研究院，現(xiàn)為博士研究生，主要從事射線成像技術(shù)研究

帥茂兵，E-mail: shuaimaobing@caep.cn

2016-11-15，

2017-03-17

Supported by Science and Technology Development Fund of China Academy of Engineering Physics (No.2015B0103014)

中國(guó)工程物理研究院科學(xué)技術(shù)發(fā)展基金(No.2015B0103014)資助

XIAO Sa, male, born in 1983, graduated from China Academy of Engineering Physics in 2011, doctoral student, mainly engaged in the research of radiography

SHUAI Maobing, E-mail: shuaimaobing@caep.cn

2016-11-15, accepted date: 2017-03-17