基于最小二乘法的航空γ能譜解析

吳和喜 葛良全 羅耀耀 劉義保 余 飛

1（東華理工大學 核技術應用教育部工程研究中心 南昌 330013）2（東華理工大學 核科學與工程學院 南昌 330013）3（成都理工大學 成都 610059）

基于最小二乘法的航空γ能譜解析

吳和喜1,2,3葛良全3羅耀耀3劉義保1,2余 飛2

1（東華理工大學 核技術應用教育部工程研究中心 南昌 330013）2（東華理工大學 核科學與工程學院 南昌 330013）3（成都理工大學 成都 610059）

在天然輻射環(huán)境下，無法找到僅含單一天然放射性系列的近似無限大地層用于獲取航空γ譜儀的響應譜，這是導致無法利用最小二乘法解析航空γ能譜的關鍵問題所在。本文基于介質互換原理、地層模塊化與分能區(qū)多次模擬的思想，通過構建Monte Carlo組合模型，從而獲得了含不同放射性物質的任意形狀地層上空航空γ譜儀的響應譜。依據(jù)航空γ能譜測點上γ輻射來源與石家莊動態(tài)帶分析，建立了一套可靠的最小二乘法航空γ能譜解析模型。通過對野外10 km測線分析發(fā)現(xiàn)，當?shù)貙雍看笥?×10-6時，該方法與標準方法分析結果相當。同時該方法為人工放射性核素響應譜的剝離提供參考。

航空γ譜儀，Monte Carlo模擬，最小二乘法

1976年，國際原子能機構(International Atomic Energy Agency, IAEA)將逆矩陣法作為航空γ能譜分析地表層天然放射性核素的標準方法[1]，該方法的核心思想[2]為：在扣除其它因素（宇宙射線、飛行器與攜帶的儀器設備和大氣氡）的影響后，天然γ輻射環(huán)境下航空γ能譜某一道區(qū)內的計數(shù)為地表層平衡235U&238U（含量比按天然組成）系、平衡釷系和40K對該道區(qū)的貢獻之和。在能量線性且不存在譜線漂移的情況下，可將逆矩陣法中某一能區(qū)計數(shù)響應規(guī)律引申到各道，這就是逐道最小二乘法[2]，核心思想為采用最小二乘法擬合解析航空γ能譜。該方法已被廣泛應用于能量色散X熒光分析儀中，如Campbell[3]、Gardner[4]和Gysel[5]等建立起能量色散X熒光分析儀對單元素響應譜的擬合函數(shù)，通過最小二乘法與基體效應校正方案獲得待測樣品中各元素的含量[5]。通過分析，上述方法能否有效運用于解析任意測點上的航空γ能譜，其關鍵在于準確獲得僅含上述三種天然放射性之一且近似無限大地層情況下的不同高度的航空γ譜儀響應譜，但實際中不存在此類地層。

蒙特卡羅(Monte Carlo, MC)數(shù)值模擬方法是國際上公認的解決粒子輸運問題的有效技術，并被廣泛應用于輻射監(jiān)測、核反應堆設計和核技術應用等領域[6-8]。作者在MCNP5軟件平臺上建立球殼模型初步解決了含任意放射性核素的近似無限大地層上不同高度航空γ譜儀的響應譜模擬問題[9]。深入研究發(fā)現(xiàn)，該模型未考慮到遠距離地層所發(fā)射的γ射線對航空γ能譜低能段的影響。為此，本文基于介質互換原理、模塊化地層與能量分區(qū)多次模擬的思想建立組合模型來解決上述問題。依據(jù)航空γ探測點上γ輻射組成，建立起天然輻射環(huán)境下基于MC模擬響應譜與最小二乘法擬合的航空γ能譜解析計算方案，并以野外實驗驗證分析該模型的優(yōu)劣。

1 組合模型的建立

根據(jù)介質互換原理[10]，在保證如圖1(a)所示的航空γ能譜探測模型中替換空氣的介質質量吸收系數(shù)與探測高度H的乘積與空氣的相應量乘積相一致時，航空γ譜儀的響應譜將保持不變。在實踐中，常采用木板等材料替代空氣進行不同高度的航空γ譜儀刻度[10-11]，但這僅能保證對單一能量特征γ射線的全能峰計數(shù)響應規(guī)律一致。而MC模擬時可通過增加空氣的密度來降低H值，此時能保證介質互換原理對整個測量能區(qū)內所有特征γ射線均成立。

圖1 航空γ能譜儀測量MC模擬模型演化圖(a) 航空γ能譜探測幾何模型，(b) 地層縱向模塊化，(c) 地層橫向模塊化Fig.1Evolutionary process of MC model of airborne gamma-ray spectrometry surveying. (a) Detection geometry model of airborne gamma-ray spectrometry, (b) Vertical modularization of ground surface, (c) Horizontal modularization of ground surface

可以看出，原來H=120 m、地層厚度L=1.5 m、地層半徑R≈10×(H+L)=1 215 m的近似無限大航空γ探測地層，用密度提高100倍的空氣介質互換變?yōu)镠=1.2 m、L=1.5 m、R≈10×(H+L)=27 m的地層。利用MCNP5軟件（本文均采用F8卡記錄，能譜刻度參數(shù)詳見文獻[9]）對介質互換前后的模型分別模擬并比對發(fā)現(xiàn)，因地層發(fā)射的γ射線在空氣中的衰減一致，模擬結果的統(tǒng)計誤差改善微弱。但該方案的優(yōu)勢在于模擬近似無限大地層時，所需的地層半徑R很小；劣勢在于介質互換原理在點探測器下成立，倘若介質互換后探測器大小不變，需通過探測效率值來校正模擬譜。

根據(jù)中心極限定理可知，模擬結果和真值間的偏差與模擬源粒子數(shù)目的根號值成反比，說明計算粒子越多結果越精確。根據(jù)圖1(b)將模擬地層以ΔL為間隔劃分成多層，將每一層按圖1(c)劃分成多個邊長為ΔR的正六棱柱模塊（容易構成近似圓柱體地層，且采用MCNP5內虛圓柱體源抽樣時源粒子的損失最少）。從圖1(c)可以看出，按ΔL=10 cm、ΔR=2 m將上述H=120 m下的變換地層劃分成1 905個正六棱柱模塊。按照基于點探測器的γ場理論，僅需模擬105個正六棱柱模塊。但實際使用的航空γ能譜儀探測器為長方體NaI(Tl)晶體，考慮到不同位置地層帶來的角響應差異和長方體探測器的軸對稱性，需模擬一個象限內510個正六棱柱模塊。

同時依據(jù)文獻[12]，對平衡235U&238U系而言，取衰變幾率大于1%的28條γ射線作為輸入源，模擬每個正六棱柱模塊時將28條射線分成7簇（保證衰變幾率大于5%的特征γ射線分開模擬，每一簇的總衰變幾率相近）單獨模擬，則得到僅含平衡235U&238U系的上述地層上120 m航空γ譜儀響應譜需模擬計算3570次，模擬結束后將所有結果線性疊加后乘4即可。該方法極大提高了模擬地層的源抽樣粒子數(shù)，但計算量很大。基于實驗室32臺計算機（4核處理器Inter i5-2400，內存2.99 GB，晶振頻率3.10 GHz）搭建4個并行MCNP5軟件計算平臺[13]（一個計算平臺32個CPU核，單核速度效率約為80%），計算上述3 570個響應譜（源抽樣粒子數(shù)為3.57×1012，總計數(shù)不確定度小于1.98%）合計21天19.25小時。對平衡釷系（源抽樣粒子數(shù)為3.06×1012，總計數(shù)不確定度小于1.67%）和40K（源抽樣粒子數(shù)為5.1×1011，總計數(shù)不確定度小于0.78%）同樣采用上述計算方案。

利用組合模型及球殼模型獲得三種天然放射性系列地層上120 m航空γ譜儀響應譜如圖2所示（為方便顯示與比對，將平衡235U&238U系按1764.494keV峰區(qū)最大計數(shù)校正為10；平衡釷系按2614.533 keV峰區(qū)最大計數(shù)校正為1000；40K的1460.83 keV峰區(qū)最大計數(shù)校正為1），可以看出兩種模擬模型獲得的譜線受統(tǒng)計漲落影響很小，30-256道計數(shù)吻合，但在1-29道計數(shù)差異明顯。將圖2中兩種模型模擬結果對應道計數(shù)相減后歸一化，結果如圖3所示，可以看出三種地層的1-29道差異相同，高能段的些許差異主要源于統(tǒng)計漲落的影響。

由于“譜平衡”成分與入射能量無關，采用上述模擬方案得到圖1中第17-25圈正六棱柱地層（僅含1 460.83 keV的特征γ射線）上120 m航空γ譜儀響應譜，將其歸一化后各道計數(shù)乘1.34（通過最大值相等的原則計算獲得）如圖3所示，其與上述兩種模擬模型結果差異譜歸一化基本吻合，證實上述兩種模型差異的主要來源為遠距離地層產生的“譜平衡”成分。

圖2 兩種模型航空γ譜儀響應譜模擬結果比對Fig.2Comparison of the response spectra of airborne gamma-ray spectrometry by the two MC models.

圖3 兩種模型航空γ譜儀響應譜模擬結果差異來源Fig.3The source of difference among the response spectra of airborne gamma-ray spectrometry by the two MC models.

2 最小二乘法航空γ能譜解析

根據(jù)天然輻射環(huán)境下航空γ譜儀探測點上的γ輻射來源[14]，可知無限大地層上空H高度的航空γ能譜中第x道計數(shù)yH,x表達式如下：式中：Sx為飛行器與攜帶的儀器設備對航空γ能譜第x道計數(shù)的貢獻；DH,x為宇宙射線對H探測高度上航空γ能譜第x道計數(shù)的貢獻；RH,x為大氣氡對H探測高度上航空γ能譜第x道計數(shù)的貢獻；UH,x、TH,x和KH,x分別表示MC模擬得到僅含平衡235U&238U系、平衡釷系和40K的地層上H探測高度上航空γ譜儀響應譜（特征峰區(qū)最大計數(shù)歸一化）第x道的計數(shù)；cU、cTh和cK為待擬合參數(shù)；εx為航空γ能譜第x道實測計數(shù)與上述各組成總計數(shù)率間的差值。

式(1)中Sx+DH,x+RH,x可直接采用測區(qū)內或周邊大面積水域上空相同高度航空γ能譜第x道計數(shù)來表征[14]。根據(jù)最小二乘法原理，cU、cTh和cK參數(shù)可通過感興趣連續(xù)譜段x∈[i, j]內2εx/yH,x之和最小獲得，擬合程序中選用粒子群與模擬退火優(yōu)化混合算法[15]來求解。

通過對石家莊動態(tài)帶8個探測高度上航空γ能譜數(shù)據(jù)（300 s累積平均譜）擬合，典型高度的擬合效果如圖4所示，可見各道數(shù)據(jù)高度吻合；得到參數(shù)cU、cTh和cK隨H的變化規(guī)律如圖5所示，利用負指數(shù)函數(shù)擬合（擬合優(yōu)度大于97.36%）得到三者的高度校正系數(shù)μU、μTh和μK分別為0.004 671 m-1、0.005948 m-1和0.007 325 m-1，與運用文獻[14]所述方法分析得到的結果(0.005163 m-1、0.005706 m-1和0.007107 m-1[16])間的相對偏差分別為-9.53%、4.24%和3.06%，間接證實本文所述解析方法具備相當高的準確度。

圖4 石家莊動態(tài)帶90 m高空航空γ能譜擬合效果Fig.4Fitting effect of airborne gamma-ray spectrometry at the dynamic scale model in Shijiazhuang City when H=90 m.

圖5 石家莊動態(tài)帶上參數(shù)cU、cTh和cK的擬合結果Fig.5Fitting values of cU, cThand cKat the dynamic scale model in Shijiazhuang City.

在利用上述方案對其它區(qū)域測點的航空γ能譜分析得到cU、cTh和cK的擬合值后，根據(jù)相對分析方法可計算測點下地層中eU、eTh和40K的含量MU、MTh和MK（單位為1），具體如下：

式中：常量2.03×10-6、10.44×10-6和0.017 9分別為石家莊動態(tài)帶上eU、eTh和40K的含量[14]；2.604、5.5934和26.672 1分別為圖5中cU、cTh和cK的值擬合線在縱軸上的截距。

3 野外應用效果分析

2006年8月中國國土資源航空物探遙感中心利用GR-820航空γ譜儀對內蒙古白云鄂博地區(qū)進行了高品質的航測。由于天然放射性核素的半衰期很長且在土壤中的遷移速度很慢，9年前后天然放射性核素的含量及分布可認為不變。選擇航測區(qū)域內穿越牧場、淺表層無人類挖掘痕跡的兩條共計約10km的測線進行比對分析。利用文獻[14]所述方法和本文擬合方法對航空γ能譜進行分析，結果分別如圖6中虛線和實線所示。同時以25 m為間隔在測線上進行地面γ能譜儀測量，結果如圖6中點所示?？梢钥闯觯瑑煞N航空γ能譜分析方法得出的天然放射性核素含量與地面γ能譜實測結果變化趨勢基本一致，兩條測線均在石英巖上存在明顯的低能異常。

因航空γ能譜單譜測量時間僅為1 s，當?shù)貙臃派湫院枯^低時特征峰區(qū)計數(shù)很少，導致利用全譜數(shù)據(jù)擬合時結果偏離真實，如小于4×10-6的eU含量擬合結果；但在地層含量大于4×10-6時，擬合結果基本能與文獻[14]所述的標準方法媲美，如eTh含量的擬合結果；而地層中放射性含量越高，擬合結果的跟隨能力越強，如40K含量的擬合結果。

4 結語

本文基于介質互換原理、地層模塊化及分能區(qū)多次模擬的思想建立任意形狀地層上空不同高度航空γ譜儀響應譜的MCNP5模擬組合模型，該模型充分解決了球殼模型未考慮遠距離地層發(fā)射的γ射線對低能譜段的影響。依據(jù)航空γ能譜測點上γ輻射來源，建立了航空γ能譜最小二乘法解析模型。對石家莊動態(tài)帶上8個探測高度航空γ能譜分析發(fā)現(xiàn)，擬合譜與實測譜間高度吻合且高度衰減規(guī)律與標準方法分析結果在±10%內吻合，證實組合模型及最小二乘擬合模型準確可靠。

依據(jù)動態(tài)帶含量與擬合參數(shù)的高度衰減函數(shù)，結合相對分析方法構建出地層內eU、eTh和40K的含量計算模型。通過對內蒙古白云鄂博航測區(qū)內兩條5 km的測線分析發(fā)現(xiàn)：地層放射性含量大于4×10-6時，擬合結果基本能與航空γ能譜標準分析方法媲美；隨著含量的升高，擬合結果對地層含量的跟隨能力越強；而在地層含量低于4×10-6時，解析結果偏離真實，后續(xù)研究時可僅對高能譜段(1.36-2.77 MeV)進行擬合，即采用最小二乘-逆矩陣法提高解譜精度。

圖6 內蒙古實驗區(qū)兩條測線三種方法的分析結果比對(a) A測線eU含量，(b) A測線eTh含量，(c) A測線40K含量，(d) B測線eU含量，(e) B測線eTh含量，(f) B測線40K含量Fig.6 Comparison of the values using three analyzing methods at an experimental zone in Inner Mongolia. (a) Content of eU in line A, (b) Content of eTh in line A, (c) Content of40K in line A, (d) Content of eU in line B, (e) Content of eTh in line B, (f) Content of40K in line B

因該方法能獲得測點上空航空γ能譜擬合結果，后續(xù)可用其提取人工放射性核素的航空γ譜儀響應譜，為準確圈定人工放射性輻射污染區(qū)域提供技術支撐。

致謝感謝中國國土資源航空物探遙感中心的萬建華高工提供的航空γ能譜測量數(shù)據(jù)，感謝曾為該數(shù)據(jù)付出辛勞的全體成員。

1 International Atomic Energy Agency. Airborne gamma ray spectrometer surveying[R]. Technical Reports Series No.323, Vienna, 1991

2 龐巨豐. γ能譜數(shù)據(jù)分析[M]. 西安: 山西科學技術出版社, 1990: 424-485

PANG Jufeng. Gamma-ray energy spectrum data processing[M]. Xi’an: Shanxi Scientific and Technical Publishers, 1990: 424-485

3 Campbell J L, Mcdonald L, Hopman T, et al. Simulations of Si(Li) X-ray detector response[J]. X-ray Spectrometry, 2001,30(4): 230-241. DOI: 10.1002/xrs.492

4 Gardner R P, Xu L B. Status of the Monte Carlo library least-squares (MCLLS) approach for non-linear radiation analyzer problems[J]. Radiation Physics & Chemistry, 2009,78(10): 843-851. DOI: 10.1016/j.radphyschem. 2009.04.023

5 Gysel M V, Lemberge P, Espen P V. Implementation of a spectrum fitting procedure using a robust peak model[J]. X-ray Spectrometry, 2003,32: 434-441. DOI: 10.1002/ xrs.666

6 Liu L P, Liu Y B, Wang J, et al. The influence of reactor core parameters on effective breading coefficient keff[J]. Chinese Physics B, 2008,17(3): 996-1000. DOI: 10.1088/1674-1056/17/3/026

7 魏強林, 楊波, 王毅, 等. 含硼鋼對慢中子衰減性能的蒙特卡羅模擬[J]. 核技術, 2010,33(5): 367-369

WEI Qianglin, YANG Bo, WANG Yi, et al. M-C simulation of slow neutron attenuation in boroncontaining stainless steel[J]. Nuclear Techniques, 2010,33(5): 367-369

8 劉合凡, 葛良全, 謝希成, 等. 基于蒙特卡羅方法的XRF探測器立體角分析[J]. 核技術, 2015,38(6): 060502. DOI: 10.11889/j.0253-3219.2015.hjs.38.060502

LIU Hefan, GE Liangquan, XIE Xicheng, et al. Calculation of desired X-ray collection angle on XRF analyzer designed by Monte Carlo method[J]. Nuclear Techniques, 2015,38(6): 060502. DOI: 10.11889/ j.0253-3219.2015.hjs.38.060502

9 吳和喜, 楊秀英, 葛良全, 等. 求解航空γ能譜模擬中的深穿透問題[J]. 核技術, 2016,39(7): 070201. DOI: 10.11889/j.0253-3219.2016.hjs.39.070201

WU Hexi, YANG Xiuying, GE Liangquan, et al. Solving the deep penetration problem in Monte Carlo simulation of airborne gamma-ray spectrum[J]. Nuclear Techniques, 2016,39(7): 070201. DOI: 10.11889/j.0253-3219.2016. hjs.39.070201

10 盧存恒, 劉慶成, 韓長青. 空間γ場的彈性變化及應用[M]. 北京: 原子能出版社, 2006: 1-25

LU Cunheng, LIU Qingcheng, HAN Zhangqing. Elasticity change of spatial gamma field and application[M]. Beijing: Atomic Energy Press, 2006: 1-25

11 胡明考, 張積運, 王新興, 等. 大型人工核素平面源校準裝置的研建及應用[R]. 中國科學技術進展報告(第一卷), 2009

HU Mingkao, ZHANG Jiyun, WANG Xinxing, et al. Establishment and application of a large calibration device of artificial radionuclide plane source[R]. Progress Report on China Nuclear Science & Technology (Vol.1), 2009

12 Heath R L. Development of new expanded editions of the Gamma-ray Spectrum Catalogues in electronic format[J]. Journal of Radioanalytical & Nuclear Chemistry, 1998,233(1): 81-89. DOI: 10.1007/BF02389651

13 楊磊, 劉玉娟, 吳和喜, 等. 用于蒙特卡羅并行計算的小型Windows PC集群建設方案[J]. 核技術, 2012,35(1): 55-59

YANG Lei, LIU Yujuan, WU Hexi, et al. A solution of small-scaled PC cluster for Monte Carlo parallel computation based on Windows OS[J]. Nuclear Techniques, 2012,35(1): 55-59

14 EJ/T 1032-2005: 航空伽瑪能譜測量規(guī)范[S]. 中華人民共和國核行業(yè)標準, 2005: 1-25

EJ/T 1032-2005: Specification for airborne gamma-ray spectrometry[S]. The People’s Republic of China Nuclear Trade Standards, 2005: 1-25

15 鄭申海, 胡小兵, 鄭滿滿, 等. 改進粒子群和模擬退火混合算法及其應用[J]. 計算機技術與發(fā)展, 2013,23(7): 26-30. DOI: 10.3969/j.issn.1673-629X.2013.07.007

ZHENG Shenhai, HU Xiaobing, ZHENG Manman, et al. An improved hybrid algorithm based on particle swarm optimization and simulated annealing and its application[J]. Computer Technology and Development, 2013,23(7): 26-30. DOI: 10.3969/j.issn.1673-629X. 2013.07.007

16 萬建華. 復雜地形區(qū)航空γ能譜數(shù)據(jù)處理方法研究[D].北京: 中國地質大學(北京), 2012

WAN Jianhua. Study on airborne gamma-ray spectrometry data processing method in complex terrain area[D]. Beijing: China University of Geosciences (Beijing), 2012

Analyzing airborne gamma-ray spectrum by the least square method

WU Hexi1,2,3GE Liangquan3LUO Yaoyao3LIU Yibao1,2YU Fei2
1(Engineering Research Center of Nuclear Technology Application,Ministry of Education,East China University of Technology,Nanchang 330013,China) 2(School of Nuclear Science and Engineering,East China University of Technology,Nanchang 330013,China) 3(Chengdu University of Technology,Chengdu 610059,China)

Background：It is difficult to establish the least square method for parsing airborne gamma-ray spectrum because the approximately infinite stratum does not exist when it contains only235U and238U series,232Th series or40K.Purpose：This study aims to acquire the response spectrum of airborne gamma-ray spectrometry (AGS) by Monte Carlo simulation and establish the above-described method.Methods：Based on the principle of media interchange, the basic thought of modular stratum and energy partitioning for characteristic gamma-ray, the combination model is used to obtain the AGS’ response spectrum of any shape of stratum by Monte Carlo simulation. And the least square model for solving airborne gamma-ray spectrum is established by the above response spectrum.Results：The corrected model of eU, eTh and40K content in stratum is acquired by analyzing the eight airborne gamma-ray spectra from the dynamic scale model in Shijiazhuang City. The calculation values of 10 km survey lines coincide with the specification for AGS when the content of natural radioactivity nuclide is larger than 4×10-6.Conclusion：The results verify that the calculation model is a valid analytical technique of AGS for naturalradioactivity nuclide in stratum. And it provides a valid method to acquire the AGS’ response spectrum of artificial radioactivity nuclide.

Airborne gamma-ray spectrometry, Monte Carlo simulation, Least square method

WU Hexi, male, born in 1985, graduated from East China University of Technology with a master’s degree in 2010, doctoral student, lecturer,

GE Liangquan, E-mail: glq@cdut.edu.cn

TL817+.2

10.11889/j.0253-3219.2016.hjs.39.110201

國家自然科學基金(No.11505027、No.11665001、No.41604116)、核技術應用教育部工程研究中心基金(No.HJSJYB2014-7&8)資助

吳和喜，男，1985年出生，2010年于東華理工大學獲碩士學位，現(xiàn)為博士研究生，講師，主要從事核探測及其數(shù)據(jù)處理研究

葛良全，E-mail: glq@cdut.edu.cn

Supported by National Natural Science Foundation of China (No.11505027, No.11665001, No.41604116), Open-ended Foundation from Engineering

Research Center of Nuclear Technology Application, Ministry of Education (No.HJSJYB2014-7&8)

mainly engaged in the research of nuclear detection and data processing

2016-09-05，

2016-09-26