凸度測量過程中點擴展函數的蒙特卡羅模擬與解析表達研究

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

（清華大學 核能與新能源技術研究院 核檢測技術北京市重點實驗室，北京 100084）

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

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

1 點擴展函數的蒙特卡羅模擬計算［3］

1.1 計算模型

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

圖1 PSF MC模擬計算模型Fig.1 MC simulation model of PSF

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

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

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

其中：μ 為線性衰減系數；D 為鋼板厚度；μm 為質量衰減系數；Dm為質量厚度。

1）當入射射線為單一能量時，PSF 隨射線能量的變化

鋼板厚度設置為8mm，分別模擬計算出不同能量的射線作用下，探測器接收的散射、透射光子概率分布情況，計算結果如圖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減小，發生散射的概率減小，但由于低能段光子更易發生光電效應和多次散射，被衰減吸收的較多，導致最終出射的散射光子減少，因此PSF 的幅值隨能量的增大稍有增大。由圖3可看出，透射光子的概率隨能量的變化呈多次曲線形狀，且MCNP 模擬計算值與理論值完全符合。由圖4可看出，探測位置離入射中心點越近，散射光子概率越大，即所接收發生散射的光子越多；當射線能量在0.06～0.18 MeV之間時，E 越大，散射量越大，也呈多次曲線形狀。

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

2）單能射線與連續能譜射線作用下PSF的計算

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

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

表1 能譜分段數據Table 1 Division data of energy spectrum

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

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

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

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

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

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

2 PSF的解析表達

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

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

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

其中，C1、C2、S1、S2均為參數。

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

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

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

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

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

3 結論

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

［1］ 張峰，閆鑌，李建新，等.工業X-CT 散射校正技術綜述［J］.CT 理論與應用研究，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］ 裴鹿成，張孝澤.蒙特卡羅方法及其在粒子輸運問題中的應用［M］.北京：科學出版社，1980.

［4］ 魏彪，馬曉昕，金煒，等.定量中子數字成像散射校正的蒙特卡羅模擬［J］.強激光與粒子束，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.