基于數值分析法的多點源核料位計放射源強度分布計算及測量的非線性修正

程曉龍，趙修良，周劍良，徐繼圓，呂 洋

(南華大學 核科學技術學院，湖南 衡陽 421001)

以γ射線技術為基礎的核子料位計已廣泛應用于工程實際。其基本原理為物料會吸收穿過的射線，射線強度隨著料位的移動而改變，通過核儀表將變化的射線強度轉換成電信號或其他形式信號來顯示和控制料位。根據不同的測量對象，放射源和探測器的配置方式可采用點源-點探測器、點源棒探測器、棒源點探測器和棒源棒探測器。對于測量范圍大、壁厚和介質對射線吸收強的對象，由于源強和測量精度的限制一般采用線狀源。線狀源最大的優點是能滿足測量要求同時實現核料位計線性刻度，在實際應用時只需一點或兩點標定即可。德國伯托公司生產的LB系列線狀源核子料位計在工程實際中應用較多，國內很早就有關于滿足線性刻度的線狀源源強分布的理論研究[1-15]，最新專利CN102735312A公開了一種60Co棒狀源的生產工藝[16]，但目前得到線源分布函數極為復雜，很難指導實際生產。

本研究根據線源的理論，提出了一種計算近似線性輸出的多點源強度分布的方法，再結合Excel和Matlab軟件進行殘差分析，進行線性輸出的修正。

1 多點源核子料位計放射源的設計

由于線源制作工藝難度大，多點源核子料位計是目前取代線源實現大尺寸、復雜工況下精度測量的最佳方法。多點源強度分布的計算可以借助線源分布函數。下面根據工程實際說明多點組合源的設計過程。

1.1 理論計算模型的建立及相關參數設置

理論計算模型相關參數及坐標系建立示意圖示于圖1。

容器尺寸為100 cm×100 cm×100 cm，壁厚5 cm,碳鋼材質，容器內物料密度為7.83 g/cm3；容器內溫度為20 ℃～120 ℃；容器內壓力為-0.05 MPa～0.1 MPa。測量技術指標：測量范圍0～90 cm；測量精度小于5 mm；料位測量上盲區小于5 cm，下盲區小于5 cm；響應時間10～60 s。

圖1 核料位計測量系統示意圖Fig.1 Nuclear gauge system schematic diagram

設定μ1為容器材料對所用γ射線的質量衰減系數，經查表[17]60Co源在鐵介質下的質量衰減系數為0.053 1 cm2/g，設定μ2為被測介質的質量衰減系數，與鐵介質近似取0.053 cm2/g。 圖1中,L為源與探測器間的水平距離，L=125 cm；D為探測器到容器左內壁的水平距離，D=105 cm；H為容器內高，H=90 cm。

1.2 多點源強度分布計算

將頂部點源和料位上限置于同一水平線，將測量范圍(5～85 cm)以8 cm均分為10等分，在左側分布11個點源。利用線源的理論，設定線源的分布函數為S(x)，則坐標為x點的點源強度為S(x)。當料面位于h高度時，探測器接收到的射線來自料面以下和以上的點源。設定探測器有效截面為S，探測效率為ε、放射源單次衰變γ光子產額η，由此得到計數率與料位的表達式：

(1)

其中Q=3.7e7εηS，B1、B2分別為寬束γ射線穿過容器壁和被測物介質的射線照射量累積因子，可查的表達式為[17]：

B=A1e-α1μR+(1-A1)e-α2μR

(2)

t1(x)、t2(x)分別是射線穿過容器壁和測量介質的實際厚度，A1與源、介質的種類和射線強度有關,可查表得到。用60Co源計算如下：

B1=23.123e(0.057 21×0.053 1×7.8×t1(x))-

22.123e(0.019 79×0.053×7.8×t1(x))

B2=23.123e(0.057 21×0.053×7.83×t2(x))-

22.123e(0.019 79×0.053×7.83×t2(x))

雖然式(1)無法得到嚴格的線性輸出，但在實際測量中一般μ2ρD比較大，針對本測量對象μ2ρD>>5，因此料面下介質對射線吸收嚴重，這樣式(1)中的第一項遠大于第二項，而且可以選擇適當的S(x)函數使式(1)第一項實現嚴格線性，從而得到計數率與料位間的近似線性函數。從式(1)可見，當?。?/p>

(3)

(4)

1.3 多點源總源強確定

料位測量有一定的誤差，設定為E，儀表采樣時間設定為t，就放射源本身而言，它所帶來的誤差只有放射性統計漲落引起的誤差，取其極限誤差：

(5)

(6)

(7)

為了保證有足夠的源強以達到要求的測量精度，把各個點源作為一個點源來考慮，并安裝在A1處，將射線束看做窄束則有：

(8)

(8)式中θ1是A1對探測器的仰角；t為二次儀表采樣時間，計算出總源強為23.421 GBq。

1.4 計算結果

利用式(3)解出各點點源的強度，與探測器組合成的計數率液位計算式為：

(9)

S(i)是對應點源的強度，x為點源的坐標，用60Co做源NaI(TI)閃爍晶體探測器，S=3 cm×3 cm探測效率為ε=30%，60Co源的η=2，通過對10段料位進行計算得到相關數據示于表1和圖2、圖3。

表1 60Co 多點源相關參數

圖2 60Co計數率料位曲線Fig.2 60Co source counting rate-material level curve

圖3 多點源棒源中各點源強度分布Fig.3 The point source activity distribution of the multipoint source

1.5 輸出特性的理論分析

由非線性分布的多點源組成的料位測量系統，其輸出并沒有達到嚴格的線性刻度。由于儀器的標定多是采用兩點標定，將料位的下限和上限所對應的探測器計數率標定料位儀表的零點和滿量程，中間按線性刻度。這樣造成儀表指示與實際料位的偏差，因此不能直接應用到對測量精度要求高的工程實際中，需要經過非線性修正。

2 非線性修正

實現計數率液位線性響應流程圖示于圖4。本研究提出三種非線性修正方法：1)對點源的源強進行補償；2)利用殘差進行軟件補償；3)用線性方程理論計算源強分布。非線性修正路線圖示于圖4。

2.1 改善源強進行線性補償

60Co計數率料位曲線示于圖2，多點源棒源中各點源強度分布示于圖3。由圖2可以看出，隨著料位h的增加，計數率的增加率降低。有資料表明[8-15]，底部點源對探測器的計數率貢獻是隨料位降低呈指數上升的，因此通過對底部點源強度補償就可得到的更好的線性輸出。具體做法是：利用Excel在原來的數據基礎上，分別自下而上給點源乘以一系列遞減的倍數，然后根據繪制出的圖形再次進行局部調整，經過調整得出較好的線性輸出，計算結果如表2和圖5、圖6所示。

圖4 實現計數率液位線性響應流程圖Fig.4 The flow chart used for realizing the counting rate-level

點源序號源位置坐標點/cm源強度/×37 MBq不同液位坐標/cm計數率/s-11-11.1942.4855 222.282-1.6730.41134 727.6337.8623.16214 234.26417.3817.84293 741.36526.9013.93373 247.78636.4311.04452 754.40745.958.90532 263.58855.487.31611 779.25965.006.11691 307.041074.525.2277854.34 1184.055.0585429.94

圖5 修正后的60Co計數率料位曲線Fig.5 60Co source counting rate-material level curve after correction

圖6 多點源棒源中各點源強度分布Fig.6 The point source activity distribution of the multipoint source

由圖5可以看出，料位計的輸出已經滿足較高的線性度，R2=0.999 6。為了達到料位計在5 cm～85 cm量程內實現5 mm的測量精度要求，根據擬合的一次函數y=-60.392x+5 496計算部分料位處所對應的計數率，結果列于表3。

表3 根據公式計算相應料位的計數率

理論上根據計數變化分辨液位的變化，該計數變化值應該達到放射性統計漲落引起的絕對偏差δ的3倍左右，否則不能準確測量液位。針對5 mm測量精度的要求，由表3可以看出,當料位從5 cm變化到5.5 cm時，引起的計數率差約為31，而此時絕對誤差約為73。為了達到精度要求可以通過增加放射源活度和測量時間來實現，當源活度擴大3.7倍為633 MBq×37 MBq，測量時間定為14 s，得到的料位變化引起的計數率差約為1 550，絕對誤差約為516,達到了可準確區分的程度。

2.2 利用殘差進行軟件非線性補償

利用Matlab對圖2進行分析，得出殘差分布曲線示于圖7。分析殘差曲線，進行三次多項式擬合結果示于圖8。

圖7 殘差分布曲線Fig.7 Residual distribution curve

圖8 殘差擬合曲線Fig.8 The residual fitting curve graph

設定通過多點源進行理論計算的計數率液位曲線函數為y1=f1(x),經最小二乘法擬合的線性函數為y2=f2(x)；殘差曲線函數為y3=f3(x)。由于殘差曲線是由理論計算曲線和擬合直線的相減得到的[18]，則有y3=y1-y2，進而線性函數y2=y1-y3，則非線性補償具體步驟如下：

(1)根據理論計算模型和設定的參數得出理論計算曲線函數y1；

(2)用最小二乘法擬合y1得到線性函數y2；

(3)通過Matlab進行殘差分析得出y3；

(4)進行兩點標定：任取兩個不同的料位(一般取測量下限和上限)x1和x2，理論計算得出f1(x1)和f2(x2)，再代入擬合的殘差方程y3，得出f3(x1)和f3(x2)，進行修正y1修=f1(x1)-f3(x1)，y2修=f1(x2)-f3(x2)；

通過在主機內設定相應的軟件程序，即可達到在不改變源活度分布的情況下，實現測量的線性輸出。根據標定曲線進行性能驗證，結果列于表4。

由表4可以看出，當料位從5 cm變化到5.5 cm時，引起的計數率差約為20，而此時絕對誤差約為62。通過增加放射源活度和測量時間來提高精度，當源活度擴大5.8倍為633 MBq×37 MBq，測量時間定為15 s，得到的料位變化引起的計數差約為1 740，絕對誤差約為576,達到了可準確區分的程度。

表4 由標定曲線計算相應料位計數率

2.3 利用線性方程理論計算點源強度分布

計算模型示于圖1，由圖1所示的坐標關系，料位上升的方向與x軸方向相同，當物料上升到x=hL時，單位時間探測器接收到的粒子數nL為[19]：

(10)

但是多數情況下，射線是寬束，因此理論計算nL=nL·B1·B2，u1、u2是容器壁和物料對射線的線衰減系數。根據線性測量的要求，探測器接收的粒子數n滿足方程ni=k·hi+n0，式中ni是料位處于hi時探測器接收到的粒子數。由于各個點源等距離分布，則hL-1-hL=hL-hL+1；

nL-1-nL=nL-nL+1=常數

(11)

聯合式(10)、(11)得到：

AL=[(eμ2·Δh/sinθ1+e-μ2·Δh/sinθ1-2)·

(12)

(12)式中，μ1、μ2、A1、θ1、Δh是已知量，L=2時可以求出A2和A1的比值，以此類推求出：

A1∶A2∶A3……A11=a∶b∶c……k

(13)

結合計算模型計算，先算出總源強為633×37 MBq再求出各項參數，結果如表5和圖9、圖10所示。

表5 多點源相關參數

圖9 計數率料位曲線Fig.9 The counting rate-material level curve

圖10 多點組合源中各點源強度分布Fig.10 The point source activity distribution of the multipoint source

由式(9)、(10)、(11)計算出的點源線性度能達到R2=0.998 2，不需要進行非線性修正。為了滿足5 cm～85 cm量程內5 mm的測量精度要求，根據擬合的線性函數y=-95.823x+8 843.7進行了部分計算，結果列于表6。

由表6可以看出,當料位從5 cm變化到5.5 cm時，引起的計數率差約為48，而此時絕對誤差約為91。通過增加放射源強度和測量時間達到精度要求，當源強度不變時，測量時間定為33 s，得到的料位變化引起的計數差約為1 584，絕對誤差約為523,達到了可準確區分的程度。

表6 根據公式計算相應液位的計數率

3 輻射防護計算

放射源的輻射防護是核子料位計關鍵，γ射線和物質作用時有三種主要形式，即光電效應、康普頓效應和電子對效應。γ射線在這三種作用的過程中逐漸被吸收，或變成另一種能量較小的光子。理論和實踐證明，屏蔽γ射線效果比較好的材料是“高原子序數的重物質”。常使用的有鉛、鐵、混凝土、巖石等，所需的厚度也各不相同。工程實際中，對核子料位計的輻射防護，鉛是用得較多的防護材料之一。

要使防護達到經濟實用的效果，屏蔽防護材料必須要有合適的厚度，本工作用一種確定厚度的方法來進行計算，在選擇好屏蔽防護位置以后，首先要確定該防護位置的最大照射量率(C/kg)，理論上可以忽略設備的遮擋，通過(14)式求得近似值。

X0=

(14)

(14)式中，X0為屏蔽防護位置離放射源最短距離R厘米處的照射量率；A為放射源的放射性活度；Γ為照射量率常數，由資料查得[16]137Cs是3.28(R·cm2/h·mCi)、60Co是13.2(R·cm2/h·mCi)。

其次確定加屏蔽防護以后要達到的照射量率標準。根據國家標準[20]，生產裝置一般執行在“距源容器的1 m區域內很少有人停留的場所”的標準，屏蔽防護后達到的照射量率小于25 μSv/h即可(如果執行其他標準，計算方法不變)。假定該值為X，那么，可以求減弱倍數K:

(15)

通過K值和放射源的γ射線的能量E，在資料[17]中可查防護層厚度，以鉛為例的結果列于表7。

表7 各向同性點源γ射線減弱K倍所需的鉛屏厚度Table 7 The lead thickness when the isotropic point source γ-ray weaken k times/cm

棒狀多點源示于圖11，由于是棒狀多點源，所以對某一點進行防護計算時要考慮所有點源的照射量率貢獻，計算圖11中各點的照射量率，并確定出最佳屏蔽層厚度，分別計算了多點組合源四周和兩端，結果列于表8。

表8 各計算點總減弱倍數K

由表8可知，側面K最大值分別為315、309和323，由表7進行插值得到對應鉛板的厚度分別為11.4 cm、11.3 cm和11.4 cm，頂部的鉛板厚度分別為10.1 cm、10.3 cm和9.8 cm，底部的鉛板厚度分別為11 cm、11 cm和11.2 cm。

4 方法對比

以上分別用三種方法，對多點源料位計的線性輸出進行了理論計算，通過調整都滿足了性能指標，相關結果列于表9。

表9 三種計算多點源滿足線性刻度的參數對比

圖11 防護計算點分布Fig.11 The distribution of the protection calculation points

由表9可以看出，三種方法計算出的源強分布都可以滿足5 cm～85 cm量程內5 mm測量精度的要求。但是在歸一化的條件下，方法1、2與3相比，可以用更短的測量時間達到更高的計數率；而方法2較方法1可以在相近的時間內達到更高的計數率。

從計算過程來看，方法3計算過程復雜但不需要修正，直接實現近似線性輸出；方法1需要通過Excel對多點組合源強度分布進行修正，雖然沒有明確的修正方法，但是通過少數次的嘗試就可以達到理想效果；方法2通過軟件來進行修正，方法明確、計算相對簡便，適用于實際工程中。

5 分析與總結

射線強度隨著所穿過的介質層厚度變化按指數規律變化，為了滿足輸出信號與料位變化成線性的要求，根據積分方程理論，多點組合源的源強分布服從如下規律：1)只有滿足線性刻度方程的源強分布，才能實現線性測量；2)需要有足夠強的底源，底源可以是一個或多個強度相同的點源；3)滿足線性刻度要求的源強分布是近似符合指數分布的積分方程。

本研究利用數值分析方法，從計數率-液位的表達式及線性計算式入手，通過等距取點的方式分別計算出11個坐標點處點源的強度值，并且根據設定的計算模型用三種方法進行了計算，通過調整相關參數，實現了計數率-液位較好的線性刻度，通過擴大源總強度和增加測量時間達到了性能指標。計算結果表明，這三種計算方法可以在理論上設計出能夠滿足線性輸出的多點組合源。

參考文獻：

[1]Wm Scott Graeme. Method of Determining Fluid Le-vel and Interface in a Closed Vessel：US，6104033[P]. 2000-08-15.

[2]Jung-Hoom Chun， Richard C.Lanza. Solid/Liquid In-terface Detection in Casting Processes by Gamma-Ray Attenuation：US，5673746[P].1997-10-07.

[3]Pogmore，Christopher Edward.Metal Level Detection Method and Apparatus：EP，0859223A1[P]. 1998-10-02.

[4]Alain Bonnemay.Apparatus for Measuring the Le-vel of a Liquid in an Enclosure：US，4591719[P]. 1986-04-27.

[5]David MW.Apparatus for and Methods of Remot-ely Monitoring Outputs of Geiger-Mueller Tubes：US, 4251810[P].1981-01-17.

[6]David J-Shoemaker.Linearization of Radiation Gaugi-ng for Measuring the Quantity of Fluidic Materials in Containers：US，3594579[P].1971-07-20.

[7]Jack GC. Tank Level Measurement System： US，3098154[P].1963-07-16.

[8]王其俊. 同位素儀表[M]. 北京： 原子能出版社， 1984:38-63.

[9]劉元襄， 侯建華， 王樹順. 輻射料位計用長棒形多點源的設計與應用[J].同位素，1998，11(1)：1-7.

Liu Yuanxiang, Hou Jianhua, Wang Shushun. The design and application of the multi-point source for the radiative level[J]. Journal of Isotopes, 1998, 11(1):1-7(in Chinese).

[10]楊介甫. LB3120型γ-液位計的線性刻度與源強分布的定量分析[J].原子能科學技術，1980，5： 588-593.

Yang Jiefu. The linear scale and the quantitative analysis of source activity distribution of the LB3120 type γ-liquid level gauge[J]. Atomic Energy Science and Technology, 1980, 5:588-593(in Chinese).

[11]楊介甫.線性刻度型γ-料位計的源活度分布檢驗[J].原子能科學技術，1985，2： 167-172.

Yang Jiefu. The source activity distribution test of the linear scale type γ-material level gauge[J]. Atomic Energy Science and Technology, 1985, 2: 167-172(in Chinese).

[12]王其俊.獲得和輻射液位計線性刻度的方法[J].原子能科學技術，1964(1):93-106.

Wang Qijun. The method of realizing radioactive liquid level gauge linear scale[J]. Atomic Energy Science and Technology, 1964(1): 93-106(in Chinese).

[13]肖德濤.結晶液面計棒狀137Cs源的設計[J].原子能科學技術，1994，28(4):378-381.

Xiao Detao. Design of the rod-like 137Cs source used in the crystallizer liquidometer[J]. Atomic Energy Science and Technology, 1994, 28(4): 378-381(in Chinese).

[14]吳禮才.料位測量用線狀放射源的線強度分布及線性測量誤差的估算[J].原子能科學技術，1982，6： 703-708.

Wu Licai. Material level measuring linear source line intensity distribution and linear measurement error estimate[J]. Atomic Energy Science and Technology, 1982, 6: 703-708(in Chinese).

[15]趙靜安.獲得直接測量輻射強度的液位計的線性刻度的兩種方法[J].原子能科學技術，1963，5：388-396.

Zhao Jingan. Two methods for achieving nuclear level gauge linear scale directly[J]. Atomic Energy Science and Technology, 1963, 5: 388-396(in Chinese).

[16]曾松柏，黃勇，雷嗣煩，等.低活度鈷-60棒狀料位源及其制備方法：CN，102735312[P].2012-10-17.

[17]李紅星.輻射防護基礎[M].北京：原子能出版社，1982.

[18]朱新巖，史忠科.基于殘差特性分析的野值檢測與剔除方法[J].飛行力學，2008，26(6)：79-83.

Zhu Xinyan, Shi Zhongke. Outlier detection method based on characteristic analyzing of residue[J]. Flight Dynamics, 2008, 26(6): 79-83(in Chinese).

[19]王德全.多點源γ-液位計源強的線性刻度分析與設計[J].同位素，1992，5(2).

Wang Dequan.The design and analysis of the linear graduation of multipoint-source intensity in γ-ray liquid level meter[J]. Journal of Isotopes, 1992, 5(2)：84-87 (in Chinese).

[20]中華人民共和國衛生部.GBZ114-2006 密封放射源及密封γ放射源容器的放射衛生防護標準[S].北京：中華人民共和國衛生部，2006：1-3.