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

吴和喜 葛良全 罗耀耀 刘义保 余 飞

1（东华理工大学 核技术应用教育部工程研究中心 南昌 330013）2（东华理工大学 核科学与工程学院 南昌 330013）3（成都理工大学 成都 610059）

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

吴和喜1,2,3葛良全3罗耀耀3刘义保1,2余 飞2

1（东华理工大学 核技术应用教育部工程研究中心 南昌 330013）2（东华理工大学 核科学与工程学院 南昌 330013）3（成都理工大学 成都 610059）

在天然辐射环境下，无法找到仅含单一天然放射性系列的近似无限大地层用于获取航空γ谱仪的响应谱，这是导致无法利用最小二乘法解析航空γ能谱的关键问题所在。本文基于介质互换原理、地层模块化与分能区多次模拟的思想，通过构建Monte Carlo组合模型，从而获得了含不同放射性物质的任意形状地层上空航空γ谱仪的响应谱。依据航空γ能谱测点上γ辐射来源与石家庄动态带分析，建立了一套可靠的最小二乘法航空γ能谱解析模型。通过对野外10 km测线分析发现，当地层含量大于4×10-6时，该方法与标准方法分析结果相当。同时该方法为人工放射性核素响应谱的剥离提供参考。

航空γ谱仪，Monte Carlo模拟，最小二乘法

1976年，国际原子能机构(International Atomic Energy Agency, IAEA)将逆矩阵法作为航空γ能谱分析地表层天然放射性核素的标准方法[1]，该方法的核心思想[2]为：在扣除其它因素（宇宙射线、飞行器与携带的仪器设备和大气氡）的影响后，天然γ辐射环境下航空γ能谱某一道区内的计数为地表层平衡235U&238U（含量比按天然组成）系、平衡钍系和40K对该道区的贡献之和。在能量线性且不存在谱线漂移的情况下，可将逆矩阵法中某一能区计数响应规律引申到各道，这就是逐道最小二乘法[2]，核心思想为采用最小二乘法拟合解析航空γ能谱。该方法已被广泛应用于能量色散X荧光分析仪中，如Campbell[3]、Gardner[4]和Gysel[5]等建立起能量色散X荧光分析仪对单元素响应谱的拟合函数，通过最小二乘法与基体效应校正方案获得待测样品中各元素的含量[5]。通过分析，上述方法能否有效运用于解析任意测点上的航空γ能谱，其关键在于准确获得仅含上述三种天然放射性之一且近似无限大地层情况下的不同高度的航空γ谱仪响应谱，但实际中不存在此类地层。

蒙特卡罗(Monte Carlo, MC)数值模拟方法是国际上公认的解决粒子输运问题的有效技术，并被广泛应用于辐射监测、核反应堆设计和核技术应用等领域[6-8]。作者在MCNP5软件平台上建立球壳模型初步解决了含任意放射性核素的近似无限大地层上不同高度航空γ谱仪的响应谱模拟问题[9]。深入研究发现，该模型未考虑到远距离地层所发射的γ射线对航空γ能谱低能段的影响。为此，本文基于介质互换原理、模块化地层与能量分区多次模拟的思想建立组合模型来解决上述问题。依据航空γ探测点上γ辐射组成，建立起天然辐射环境下基于MC模拟响应谱与最小二乘法拟合的航空γ能谱解析计算方案，并以野外实验验证分析该模型的优劣。

1 组合模型的建立

根据介质互换原理[10]，在保证如图1(a)所示的航空γ能谱探测模型中替换空气的介质质量吸收系数与探测高度H的乘积与空气的相应量乘积相一致时，航空γ谱仪的响应谱将保持不变。在实践中，常采用木板等材料替代空气进行不同高度的航空γ谱仪刻度[10-11]，但这仅能保证对单一能量特征γ射线的全能峰计数响应规律一致。而MC模拟时可通过增加空气的密度来降低H值，此时能保证介质互换原理对整个测量能区内所有特征γ射线均成立。

图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倍的空气介质互换变为H=1.2 m、L=1.5 m、R≈10×(H+L)=27 m的地层。利用MCNP5软件（本文均采用F8卡记录，能谱刻度参数详见文献[9]）对介质互换前后的模型分别模拟并比对发现，因地层发射的γ射线在空气中的衰减一致，模拟结果的统计误差改善微弱。但该方案的优势在于模拟近似无限大地层时，所需的地层半径R很小；劣势在于介质互换原理在点探测器下成立，倘若介质互换后探测器大小不变，需通过探测效率值来校正模拟谱。

根据中心极限定理可知，模拟结果和真值间的偏差与模拟源粒子数目的根号值成反比，说明计算粒子越多结果越精确。根据图1(b)将模拟地层以ΔL为间隔划分成多层，将每一层按图1(c)划分成多个边长为ΔR的正六棱柱模块（容易构成近似圆柱体地层，且采用MCNP5内虚圆柱体源抽样时源粒子的损失最少）。从图1(c)可以看出，按ΔL=10 cm、ΔR=2 m将上述H=120 m下的变换地层划分成1 905个正六棱柱模块。按照基于点探测器的γ场理论，仅需模拟105个正六棱柱模块。但实际使用的航空γ能谱仪探测器为长方体NaI(Tl)晶体，考虑到不同位置地层带来的角响应差异和长方体探测器的轴对称性，需模拟一个象限内510个正六棱柱模块。

同时依据文献[12]，对平衡235U&238U系而言，取衰变几率大于1%的28条γ射线作为输入源，模拟每个正六棱柱模块时将28条射线分成7簇（保证衰变几率大于5%的特征γ射线分开模拟，每一簇的总衰变几率相近）单独模拟，则得到仅含平衡235U&238U系的上述地层上120 m航空γ谱仪响应谱需模拟计算3570次，模拟结束后将所有结果线性叠加后乘4即可。该方法极大提高了模拟地层的源抽样粒子数，但计算量很大。基于实验室32台计算机（4核处理器Inter i5-2400，内存2.99 GB，晶振频率3.10 GHz）搭建4个并行MCNP5软件计算平台[13]（一个计算平台32个CPU核，单核速度效率约为80%），计算上述3 570个响应谱（源抽样粒子数为3.57×1012，总计数不确定度小于1.98%）合计21天19.25小时。对平衡钍系（源抽样粒子数为3.06×1012，总计数不确定度小于1.67%）和40K（源抽样粒子数为5.1×1011，总计数不确定度小于0.78%）同样采用上述计算方案。

利用组合模型及球壳模型获得三种天然放射性系列地层上120 m航空γ谱仪响应谱如图2所示（为方便显示与比对，将平衡235U&238U系按1764.494keV峰区最大计数校正为10；平衡钍系按2614.533 keV峰区最大计数校正为1000；40K的1460.83 keV峰区最大计数校正为1），可以看出两种模拟模型获得的谱线受统计涨落影响很小，30-256道计数吻合，但在1-29道计数差异明显。将图2中两种模型模拟结果对应道计数相减后归一化，结果如图3所示，可以看出三种地层的1-29道差异相同，高能段的些许差异主要源于统计涨落的影响。

由于“谱平衡”成分与入射能量无关，采用上述模拟方案得到图1中第17-25圈正六棱柱地层（仅含1 460.83 keV的特征γ射线）上120 m航空γ谱仪响应谱，将其归一化后各道计数乘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 最小二乘法航空γ能谱解析

根据天然辐射环境下航空γ谱仪探测点上的γ辐射来源[14]，可知无限大地层上空H高度的航空γ能谱中第x道计数yH,x表达式如下：式中：Sx为飞行器与携带的仪器设备对航空γ能谱第x道计数的贡献；DH,x为宇宙射线对H探测高度上航空γ能谱第x道计数的贡献；RH,x为大气氡对H探测高度上航空γ能谱第x道计数的贡献；UH,x、TH,x和KH,x分别表示MC模拟得到仅含平衡235U&238U系、平衡钍系和40K的地层上H探测高度上航空γ谱仪响应谱（特征峰区最大计数归一化）第x道的计数；cU、cTh和cK为待拟合参数；εx为航空γ能谱第x道实测计数与上述各组成总计数率间的差值。

式(1)中Sx+DH,x+RH,x可直接采用测区内或周边大面积水域上空相同高度航空γ能谱第x道计数来表征[14]。根据最小二乘法原理，cU、cTh和cK参数可通过感兴趣连续谱段x∈[i, j]内2εx/yH,x之和最小获得，拟合程序中选用粒子群与模拟退火优化混合算法[15]来求解。

通过对石家庄动态带8个探测高度上航空γ能谱数据（300 s累积平均谱）拟合，典型高度的拟合效果如图4所示，可见各道数据高度吻合；得到参数cU、cTh和cK随H的变化规律如图5所示，利用负指数函数拟合（拟合优度大于97.36%）得到三者的高度校正系数μ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 石家庄动态带90 m高空航空γ能谱拟合效果Fig.4Fitting effect of airborne gamma-ray spectrometry at the dynamic scale model in Shijiazhuang City when H=90 m.

图5 石家庄动态带上参数cU、cTh和cK的拟合结果Fig.5Fitting values of cU, cThand cKat the dynamic scale model in Shijiazhuang City.

在利用上述方案对其它区域测点的航空γ能谱分析得到cU、cTh和cK的拟合值后，根据相对分析方法可计算测点下地层中eU、eTh和40K的含量MU、MTh和MK（单位为1），具体如下：

式中：常量2.03×10-6、10.44×10-6和0.017 9分别为石家庄动态带上eU、eTh和40K的含量[14]；2.604、5.5934和26.672 1分别为图5中cU、cTh和cK的值拟合线在纵轴上的截距。

3 野外应用效果分析

2006年8月中国国土资源航空物探遥感中心利用GR-820航空γ谱仪对内蒙古白云鄂博地区进行了高品质的航测。由于天然放射性核素的半衰期很长且在土壤中的迁移速度很慢，9年前后天然放射性核素的含量及分布可认为不变。选择航测区域内穿越牧场、浅表层无人类挖掘痕迹的两条共计约10km的测线进行比对分析。利用文献[14]所述方法和本文拟合方法对航空γ能谱进行分析，结果分别如图6中虚线和实线所示。同时以25 m为间隔在测线上进行地面γ能谱仪测量，结果如图6中点所示。可以看出，两种航空γ能谱分析方法得出的天然放射性核素含量与地面γ能谱实测结果变化趋势基本一致，两条测线均在石英岩上存在明显的低能异常。

因航空γ能谱单谱测量时间仅为1 s，当地层放射性含量较低时特征峰区计数很少，导致利用全谱数据拟合时结果偏离真实，如小于4×10-6的eU含量拟合结果；但在地层含量大于4×10-6时，拟合结果基本能与文献[14]所述的标准方法媲美，如eTh含量的拟合结果；而地层中放射性含量越高，拟合结果的跟随能力越强，如40K含量的拟合结果。

4 结语

本文基于介质互换原理、地层模块化及分能区多次模拟的思想建立任意形状地层上空不同高度航空γ谱仪响应谱的MCNP5模拟组合模型，该模型充分解决了球壳模型未考虑远距离地层发射的γ射线对低能谱段的影响。依据航空γ能谱测点上γ辐射来源，建立了航空γ能谱最小二乘法解析模型。对石家庄动态带上8个探测高度航空γ能谱分析发现，拟合谱与实测谱间高度吻合且高度衰减规律与标准方法分析结果在±10%内吻合，证实组合模型及最小二乘拟合模型准确可靠。

依据动态带含量与拟合参数的高度衰减函数，结合相对分析方法构建出地层内eU、eTh和40K的含量计算模型。通过对内蒙古白云鄂博航测区内两条5 km的测线分析发现：地层放射性含量大于4×10-6时，拟合结果基本能与航空γ能谱标准分析方法媲美；随着含量的升高，拟合结果对地层含量的跟随能力越强；而在地层含量低于4×10-6时，解析结果偏离真实，后续研究时可仅对高能谱段(1.36-2.77 MeV)进行拟合，即采用最小二乘-逆矩阵法提高解谱精度。

图6 内蒙古实验区两条测线三种方法的分析结果比对(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

因该方法能获得测点上空航空γ能谱拟合结果，后续可用其提取人工放射性核素的航空γ谱仪响应谱，为准确圈定人工放射性辐射污染区域提供技术支撑。

致谢感谢中国国土资源航空物探遥感中心的万建华高工提供的航空γ能谱测量数据，感谢曾为该数据付出辛劳的全体成员。

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

2 庞巨丰. γ能谱数据分析[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: 航空伽玛能谱测量规范[S]. 中华人民共和国核行业标准, 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]. 计算机技术与发展, 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 万建华. 复杂地形区航空γ能谱数据处理方法研究[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年于东华理工大学获硕士学位，现为博士研究生，讲师，主要从事核探测及其数据处理研究

葛良全，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