凸度測量過程中點(diǎn)擴(kuò)展函數(shù)的蒙特卡羅模擬與解析表達(dá)研究

胡克敏，吳志芳，苗積臣，談春明

（清華大學(xué) 核能與新能源技術(shù)研究院 核檢測技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100084）

清華大學(xué)核能與新能源技術(shù)研究院研制成功的凸度測量系統(tǒng)是利用X 射線法高精度測量板帶材（鋼板）橫斷面厚度分布，進(jìn)而計(jì)算出板帶材凸度、楔度及中心線厚度等參數(shù)的大型多功能軋鋼自動(dòng)化裝備。為保證±0.1%的測量精度，對測量過程中康普頓散射的影響分析及校正是一必需的關(guān)鍵環(huán)節(jié)。凸度測量系統(tǒng)中散射的影響主要分為兩部分：1）被測鋼板自身康普頓散射造成的影響，也是散射影響的主要原因；2）凸度測量系統(tǒng)各結(jié)構(gòu)部件造成的散射影響，此部分的影響較小，且錯(cuò)綜復(fù)雜，較難定量分析。本文主要針對第一部分即鋼板自身的散射影響作深入研究。

目前國內(nèi)外散射校正的方法［1］包括3類：1）硬件方法，基本原理為采用物理方法消除或減少進(jìn)入探測單元的散射射線來實(shí)現(xiàn)；2）軟件較正方法，包括卷積法、反卷積法、蒙特卡羅（MC）方法、理想模型估計(jì)等；3）混合校正方法，通過散射校正板、頻率調(diào)制［2］等硬件的設(shè)計(jì)取得散射影響分布，再融合軟件方法進(jìn)一步得到精確信息的方法。本文采用蒙特卡羅方法模擬計(jì)算凸度測量系統(tǒng)中點(diǎn)擴(kuò)展函數(shù)（PSF）與板材厚度、射線能量間的關(guān)系曲線，采用非線性擬合方法確定PSF理論模型參數(shù)，得到單能射線下凸度測量系統(tǒng)的PSF理論公式，為進(jìn)一步實(shí)現(xiàn)凸度測量系統(tǒng)的散射校正奠定基礎(chǔ)。

1 點(diǎn)擴(kuò)展函數(shù)的蒙特卡羅模擬計(jì)算［3］

1.1 計(jì)算模型

以凸度測量系統(tǒng)為參照，模擬計(jì)算能量為E的δ脈沖射線垂直照射厚度為D 的鐵板，探測面至被測鋼板下表面的距離H 為60cm（圖1）。

圖1 PSF MC模擬計(jì)算模型Fig.1 MC simulation model of PSF

本文采用MCNP 軟件對上述模型在不同射線能量、不同鐵板厚度情況下的PSF 進(jìn)行相關(guān)計(jì)算分析。另外，本文所述的PSF僅指到達(dá)探測器的散射光子數(shù)占射線源光子總數(shù)的概率（或稱散射光子的貢獻(xiàn)），不包括透射光子。

1.2 射線能量、鋼板厚度對PSF的影響

參照凸度測量系統(tǒng)，射線能量范圍設(shè)置為0～0.3MeV，測量厚度范圍為0～16mm，探測面到鋼板下表面的距離為60cm。由射線與物質(zhì)的相互作用宏觀經(jīng)驗(yàn)公式可得到透射光子概率為：

其中：μ 為線性衰減系數(shù)；D 為鋼板厚度；μm 為質(zhì)量衰減系數(shù)；Dm為質(zhì)量厚度。

1）當(dāng)入射射線為單一能量時(shí)，PSF 隨射線能量的變化

鋼板厚度設(shè)置為8mm，分別模擬計(jì)算出不同能量的射線作用下，探測器接收的散射、透射光子概率分布情況，計(jì)算結(jié)果如圖2～4所示。

圖2 不同射線能量下PSF沿探測面的分布曲線Fig.2 Distribution curve of PSF with detection surface under different ray energy

圖3 透射光子概率隨射線能量的變化Fig.3 Transmission probability vs ray energy

由圖2可看出，隨射線能量從60keV 增加到300keV，PSF的幅值逐漸增大，但增長幅度漸緩，且曲線形狀保持軸對稱不變。從物理分析的角度，雖然隨著能量的增加，康普頓散射的截面Σs減小，發(fā)生散射的概率減小，但由于低能段光子更易發(fā)生光電效應(yīng)和多次散射，被衰減吸收的較多，導(dǎo)致最終出射的散射光子減少，因此PSF 的幅值隨能量的增大稍有增大。由圖3可看出，透射光子的概率隨能量的變化呈多次曲線形狀，且MCNP 模擬計(jì)算值與理論值完全符合。由圖4可看出，探測位置離入射中心點(diǎn)越近，散射光子概率越大，即所接收發(fā)生散射的光子越多；當(dāng)射線能量在0.06～0.18 MeV之間時(shí)，E 越大，散射量越大，也呈多次曲線形狀。

圖4 不同探測位置處散射量隨射線能量的變化曲線Fig.4 Curve of scattering probability vs ray energy at different detection positions

2）單能射線與連續(xù)能譜射線作用下PSF的計(jì)算

實(shí)際上，凸度測量系統(tǒng)所用的X 射線源具有連續(xù)能譜，因此，研究在連續(xù)能譜射線作用下PSF的計(jì)算很有必要。但直接對連續(xù)能譜射線進(jìn)行分析計(jì)算較困難，且管電壓和管電流不穩(wěn)定，也會(huì)造成能譜的較大差值，因此，本文設(shè)想采用多個(gè)單能計(jì)算結(jié)果的線性組合來代替對整個(gè)連續(xù)能譜的分析計(jì)算，即首先選擇若干典型的單能射線，分別代表整個(gè)連續(xù)能譜的不同能段；然后，分析計(jì)算出每個(gè)單能射線下的MCNP模擬結(jié)果（如散射光子分布）；再將這些單能射線下的計(jì)算結(jié)果進(jìn)行加權(quán)求和，作為整個(gè)連續(xù)能譜作用下的MCNP計(jì)算結(jié)果。

針對最大能量為180keV 的X 光機(jī)連續(xù)能譜，本文選取一簡單模型（僅選4個(gè)單能值）進(jìn)行分析（表1），鋼板厚度取8mm。分別計(jì)算單能情況下加權(quán)求和的結(jié)果和直接用整個(gè)連續(xù)能譜計(jì)算的結(jié)果（圖5、6），可看出，兩者基本吻合，但存在約10%的差距，對8mm 厚度鋼板，10%的差距對散射的影響在萬分之幾范圍，小于凸度測量系統(tǒng)千分之一的精度要求，且通過選用更多的單能射線及優(yōu)化權(quán)重因子，可獲得更好的模擬結(jié)果，因此，在實(shí)際應(yīng)用中，可考慮用單能的線性組合代替實(shí)際能譜計(jì)算。

表1 能譜分段數(shù)據(jù)Table 1 Division data of energy spectrum

圖5 散射光子概率分布結(jié)果對比Fig.5 Contrast of scattering photon probability distribution

圖6 總探測光子概率分布結(jié)果對比Fig.6 Contrast of total detection photon probability distribution

3）PSF隨鋼板厚度變化分析

選用射線能量為0.18 MeV，分別計(jì)算鋼板厚度D 為1～16mm 的PSF，計(jì)算結(jié)果如圖7～9所示。由圖7 可看出，隨著鐵板厚度從1mm增大到10mm，散射的PSF（即散射光子的概率分布）逐漸增大，且隨著厚度的增加，PSF增大的幅度卻逐漸減小，主要原因是隨著厚度的增大，光子發(fā)生多次散射而被吸收的比例增多，導(dǎo)致散射后最終出射的光子減少。由圖8可看出，MCNP模擬計(jì)算的透射概率和理論值較吻合，符合射線衰減規(guī)律。由圖9可看出，在0.18 MeV 射線下，PSF 隨鋼板厚度D的增大先增大后減小，距入射中心點(diǎn)越近，幅值越大，散射光子越多。

圖7 PSF沿探測面的分布曲線Fig.7 Distribution curve of PSF with detection surface

圖8 透射光子概率隨鋼板厚度的變化Fig.8 Transmission probability vs steel thickness

2 PSF的解析表達(dá)

上述PSF 的分析研究是基于數(shù)值表達(dá)的，若能把PSF用解析形式表達(dá)出來，將更有利于實(shí)際應(yīng)用。通過對PSF的曲線形狀分析可知，其近似于高斯分布，若用1 個(gè)高斯函數(shù)表述，其結(jié)果與實(shí)際偏差較大，因此，本文選擇采用兩個(gè)高斯函數(shù)的線性組合來擬合PSF，這一點(diǎn)在SPECT、PET 散射校正中已得到驗(yàn)證，在中子成像中也有較好的效果［4］。用兩個(gè)高斯函數(shù)的線性組合模型來近似PSF，其二維形式［5］為：

圖9 不同探測位置處散射量隨鋼板厚度的變化曲線Fig.9 Curve of scattering probability with steel thickness under different detection positions

對PSF沿y 軸方向積分，得到PSF的一維積分形式為：

其中，C1、C2、S1、S2均為參數(shù)。

對于理想單能準(zhǔn)直射線束，參數(shù)C1、C2、S1、S2與射線能量、被測鋼板性質(zhì)、被測鋼板厚度、被測鋼板到探測面的距離（空氣隙）相關(guān)。

針對凸度測量系統(tǒng)，探測器到鋼板下表面的距離為設(shè)定值60cm，假定射線能量為180keV，則參數(shù)C1、C2、S1、S2僅與鋼板厚度相關(guān)，即：

通過蒙特卡羅方法模擬，分別計(jì)算出厚度為1～16mm 的鋼板PSF的結(jié)果，然后采用非線性最小二乘原理和方法，用式（5）進(jìn)行曲線擬合，求出分別對應(yīng)于不同厚度鋼板的參數(shù)C1、C2、S1、S2的最佳擬合估計(jì)值（表2），均方差R 越趨近于1，相干系數(shù)RMSE 越趨近于0，所得擬合值就越好。各參數(shù)與鋼板厚度的曲線關(guān)系如圖10所示。對各曲線采用多項(xiàng)式擬合的方法，可得到凸度儀系統(tǒng)在射線能量為180keV 時(shí)PSF各參數(shù)與鋼板厚度的關(guān)系式，據(jù)此，就可計(jì)算出任意鋼板厚度下的PSF分布，從而進(jìn)行散射校正，提高測量精度。

表2 參數(shù)C1、C2、S1、S2 擬合估計(jì)值及其均方差和相干系數(shù)Table 2 Fitting parameters，mean square error and correlation coefficient

圖10 各參數(shù)與鋼板厚度的擬合曲線Fig.10 Fitting curve of parameters and steel thickness

3 結(jié)論

本文通過蒙特卡羅模擬方法，研究了鋼板凸度測量中散射影響的PSF及其影響因素，在此基礎(chǔ)上，采用雙高斯函數(shù)模型，對PSF 的解析表述方法進(jìn)行了分析探究，給出180keV 射線情況下點(diǎn)擴(kuò)展函數(shù)的解析表達(dá)形式，這對鋼板凸度測量中的散射校正問題提供了技術(shù)支持。

［1］ 張峰，閆鑌，李建新，等.工業(yè)X-CT 散射校正技術(shù)綜述［J］.CT 理論與應(yīng)用研究，2009，18（4）：34-43.ZHANG Feng，YAN Bin，LI Jianxin，et al.Review of scatter correction on X-ray industrial computed tomography［J］.CT Theory and Applications，2009，18（4）：34-43（in Chinese）.

［2］ ZHU L，BENNETT N R.Scatter correction method for X-ray CT using primary modulation：Theory and preliminary results［J］.IEEE Transactions on Medical Imaging，2006，25（12）：1 573-1 587.

［3］ 裴鹿成，張孝澤.蒙特卡羅方法及其在粒子輸運(yùn)問題中的應(yīng)用［M］.北京：科學(xué)出版社，1980.

［4］ 魏彪，馬曉昕，金煒，等.定量中子數(shù)字成像散射校正的蒙特卡羅模擬［J］.強(qiáng)激光與粒子束，2008，20（9）：1 545-1 550.WEI Biao，MA Xiaoxin，JIN Wei，et al.Monte-Carlo simulation of scattering correction forquantitative neutron imaging system［J］.High Power Laser and Particle Beams，2008，20（9）：1 545-1 550（in Chinese）.

［5］ NUYTS J，BOSMANS H，SUETENS P.An analytical model for compton scatter in a homogeneously attenuating medium［J］.IEEE Transactions on Medical Imaging，1993，12（3）：421-429.