CN113189636A - 基于多能量响应的航空γ能谱分析方法 - Google Patents

Abstract

本申请涉及一种基于多能量响应的航空γ能谱分析方法。其中，首先采集不同能量(也即多能量)γ射线的能谱数据，之后依次对能谱数据中的各计数率进行本底修正、能窗剥离修正和高度衰减修正，再采用灵敏度换算法将高度衰减修正后的各计数率统一到相同单位，从而得到分析结果。如此设置，可以提高无人机航空γ能谱仪对多能量γ射线核素的探测能力，解决地面和空中测量时，原始数据修正系数的确定问题，给出地面和空中测量的结果；还可进一步推广应用于地面和地下钻井工作中的γ能谱测量。与现有方法相比，本方法可同时分析多个γ射线能窗，测量精度更高；与全谱分析法相比，本方法可以避免一些干扰谱段对测量结果的影响。

Description

基于多能量响应的航空γ能谱分析方法

技术领域

本申请涉及能谱分析技术领域，尤其涉及一种基于多能量响应的航空γ能谱分析方法。

背景技术

航空γ能谱测量是将大型γ能谱仪安装在飞机、汽艇等移动平台上，采用边飞行、边测量的方式获取移动平台航迹路线上的γ射线能谱及其对应的辅助数据的测量方法。该方法早期主要用于铀矿勘查，之后逐渐发展到基础地质调查、辐射环境调查，核应急监测等领域，对航空γ能谱测量的探测精度要求越来越高，测量源项也越来越复杂。

近几年来，新出现了以无人机为载体的无人机航空γ能谱仪，从固定翼飞机、旋翼直升机等大型有人驾驶平台向无人驾驶的智能平台(即：无人机)方向发展，对航空γ能谱仪的小型化、轻型化和智能化提出了新的挑战。相对于原来的有人机载航空γ能谱仪，无人机航空γ能谱仪的探测灵敏度显著降低。虽然测量飞行的高度能够从原来的60～180m下降到0～120m，但还是导致了单能量γ能谱分析法计算出的核素含量或活度浓度值的准确性和稳定性的明显下降。因此，亟需一种具有更高准确性和稳定性的γ能谱分析方法。

发明内容

本申请提供一种基于多能量响应的航空γ能谱分析方法，以解决目前的能谱分析方法准确性和稳定性不高的问题。

本申请的上述目的是通过以下技术方案实现的：

本申请实施例提供一种基于多能量响应的航空γ能谱分析方法，其包括：

使用多种放射性核素对无人机航空γ能谱仪进行测试，以采集不同能量γ射线的能谱数据；其中，所述能谱数据包括多个不同能窗的计数率；

对各计数率进行本底修正，以扣除仪器本身、无人机平台、宇宙射线、大气氡在能谱仪采集时产生的计数率；

对本底修正后的各计数率进行能窗剥离修正；

采用对数拟合法对能窗剥离修正后的各计数率进行高度衰减修正；

采用灵敏度换算法将高度衰减修正后的各计数率统一到含量或活度浓度量值单位，得到分析结果。

可选的，所述使用多种放射性核素对无人机航空γ能谱仪进行测试，以采集不同能量γ射线的能谱数据，包括：

对能谱仪进行设置，具体包括：

根据测试结果中能谱仪对不同能量γ射线的特征峰峰位和能量分辨率，确定能谱仪的能量刻度系数和分辨率刻度系数；

由所述能谱仪刻度系数计算能谱仪的窗宽，计算得到多能量响应的能窗位置和宽度；

基于所述能窗位置和宽度对能谱仪进行设置。

可选的，所述多种放射性核素包括K系核素、U系核素和Th系核素；K系核素对应包括K能窗，U系核素对应包括U1、U2和U3能窗，Th系核素对应包括Th1和Th2能窗；

进行能窗剥离修正时，在传统能窗剥离法的基础上，将U系和Th系核素的低能特征γ射线能窗计数率纳入计算，得到无人机飞行高度上的能窗净计数率；公式如下：

式中，n_k.k、n_u.u、n_t.t为修正后的K、U、Th能窗净计数率；n_k、n_u1、n_u2、n_u3、n_t1、n_t2为本底修正后的能窗计数率；c_k.u1、c_k.u2、c_k.u3、c_k.t1、c_k.t2为K窗对U1、U2、U3、Th1、Th2能窗的剥离系数，△c_k.u1(h_stp)、△c_k.u2(h_stp)、△c_k.u3(h_stp)、△c_k.t1(h_stp)、△c_k.t2(h_stp)为K窗剥离系数随高度的变化函数；c_u1.k、c_u1.u2、c_u1.u3、c_u1.t1、c_u1.t2为U1窗对K、U2、U3、Th1、Th2能窗的剥离系数，△c_u1.k(h_stp)、△c_u1.u2(h_stp)、△c_u1.u3(h_stp)、△c_u1.t1(h_stp)、△c_u1.t2(h_stp)为U1窗剥离系数随高度的变化函数；c_t1.k、c_t1.u1、c_t1.u2、c_t1.u3、c_t1.t2为Th1窗对K、U1、U2、U3、Th2能窗的剥离系数，△c_t1.k(h_stp)、△c_t1.u1(h_stp)、△c_t1.u2(h_stp)、△c_t1.u3(h_stp)、△c_t1.t2(h_stp)为Th1窗剥离系数随高度的变化函数；其中，c_t1.u1、c_t1.k、c_u1.k与传统剥离法中的α、β、γ一致，c_k.u1、c_k.t1、c_u1.t1与传统剥离法中的a、b、g一致。

可选的，所述剥离系数的计算过程包括：

将无人机航空γ能谱仪置于不同的航空放射性测量模型标准装置上，分别对本底模型AP-B、钾模型AP-K、铀模型AP-U、钍模型AP-Th、混合模型AP-M进行测试；

根据测试结果，按照下式计算无人机航空γ能谱仪在地面的剥离系数和灵敏度：

式中，n_k.k、n_k.u1、n_k.u2、n_k.t1、n_k.t2为探测器在AP-K模型上测到的能窗计数率扣除AP-B模型上计数率后的K、U1、U2、U3、Th1、Th2窗计数率结果，n_u.k、n_u.u1、n_u.u2、n_u.u3、n_u.t1、n_u.t2、为探测器在AP-U模型上测到的能窗计数率扣除AP-B模型上计数率后的窗计数率结果，n_t.k、n_t.u1、n_t.u2、n_t.u3、n_t.t1、n_t.t2为探测器在AP-T模型上测到的能窗计数率扣除AP-B模型上计数率后的窗计数率结果；

q_k.k、q_k.u、q_k.t为使用标准配套能谱仪定值测量得到的AP-K模型中的K、U、Th含量扣除AP-B模型中K、U、Th含量的量值，q_u.k、q_u.u、q_u.t为AP-U模型中的K、U、Th含量扣除AP-B模型中K、U、Th含量的量值，q_t.k、q_t.u、q_t.t为AP-Th模型中的K、U、Th含量扣除AP-B模型中K、U、Th含量的量值；

s_k为探测器的K窗对K元素的探测灵敏度，s_u为探测器的U1窗对U元素的探测灵敏度，s_t为探测器的Th1窗对Th元素的探测灵敏度；

c_k.u1、c_k.u2、c_k.u3、c_k.t1、c_k.t2为U1、U2、U3、Th1、Th2窗对K元素的灵敏度与K窗对K元素灵敏度的比值，c_u1.k、c_u1.u2、c_u1.u3、c_u1.t1、c_u1.t2为K、U2、U3、Th1、Th2窗对U元素的灵敏度与U1窗对U元素灵敏度的比值，c_t1.k、c_t1.u1、c_t1.u2、c_t1.u3、c_t1.t2为K、U1、U2、U3、Th2窗对Th元素的灵敏度与Th1窗对Th元素灵敏度的比值，这些参数亦称为剥离系数。

可选的，所述采用对数拟合法对能窗剥离修正后的各计数率进行高度衰减修正，包括：

采用如下公式计算高度衰减修正系数：

式中，a、b、c、d为多项式拟合系数，令：

G(h)＝ah³+bh²+ch

则有

ln(100×n_0.i)＝ln(100×n_h.i)+G(h₀)-G(h_stp)。

可选的，所述采用对数拟合法对能窗剥离修正后的各计数率进行高度衰减修正，之后还包括：

根据蒙卡模拟计算确定高度衰减修正系数的拟合误差，具体包括：

根据无人机航空γ能谱仪的几何结构，以及各结构部件制作材料的元素成分、密度，建立蒙特卡罗模拟数学模型；

根据航空放射性测量模型标准装置的物质成分、密度，建立一个半径大于2.5km，厚度在50cm以上的饱和体源数学模型；

将饱和体源数学模型设定为核素分布均匀的模拟源，其中，模拟源发射的γ射线能量、强度及子体核素分支比等参数参照国际原子能机构核数据库提供的U系、Th系子体核素，以及⁴⁰K核素；

采用方差减小技术对无人机航空γ能谱仪和饱和体源组成的数学模型进行优化，在大型科学计算平台上模拟计算无人机航空γ能谱仪在饱和体源模型中心点正上方0.135、1、5、10、20、30、60、90、120、150、180、210、240、270、300m高度时的模拟能谱数据；其中，0.135m为航空放射性测量模型标准装置静态测试时探测器的离地高度；

根据无人机航空γ能谱仪的实际刻度和能窗设置参数对模拟能谱数据进行高斯展宽，进而计算K、U1、U2、U3、Th1、Th2窗计数率，并根据在航空放射性测量模型标准装置上测得的灵敏度，对0～300m高度上的模拟数据进行归一化；

按下式计算能谱仪在饱和体源上h高度时的剥离系数：

式中，s_h.k.j、s_h.u.j、s_h.t.j为第j窗对K、U、Th元素的灵敏度；s_h.k.k为K窗对K元素的灵敏度，s_h.u.u1为U1窗对U元素的灵敏度，s_h.t.t1为Th1窗对Th元素的灵敏度；c_h.k.j、c_h.u.j、c_h.t.j为第j窗对K、U、Th元素的剥离系数；

按下式计算无人机航空γ能谱仪灵敏度随高度变化的修正参数：

式中，s_i.j(h_stp)为能谱仪在饱和体源模型h高度上，第j窗对第i元素的灵敏度；a_i.j、b_i.j、c_i.j、d_i.j为第j窗对第i元素的灵敏度修正参数；

参照IAEA-1363技术报告，根据不同高度的模拟计算结果，采用线性回归拟合得到剥离系数随高度变化的修正参数。

可选的，所述方法还包括：

采用示值误差检验法，对无人机航空γ能谱仪在地面的灵敏度和剥离系数进行检验；

以及，

采用示值误差检验法，对无人机航空γ能谱仪灵敏度随高度变化的修正参数和剥离系数随高度变化的修正系数进行检验。

可选的，所述采用示值误差检验法，对无人机航空γ能谱仪在地面的灵敏度和剥离系数进行检验，具体包括：

根据无人机航空γ能谱仪在混合模型AP-M上的测量数据，采用如下公式计算模型含量或活度浓度值；

式中，n_m.k、n_m.u1、n_m.u2、n_m.t1、n_m.t2为探测器在AP-M模型上测到的能窗计数率扣除AP-B模型上计数率后的K、U1、U2、U3、Th1、Th2窗计数率；Q_i为第i种核素的含量或活度浓度量值，对于天然放射性而言，K含量的单位一般为％、U为10^-6、Th为10^-6；对于活度浓度量值而言，K、U、Th的单位一般为Bq/kg；s_i为第i种核素的灵敏度，对于天然放射性核素而言，K含量灵敏度的单位一般为cps/％、U为cps/10^-6、Th为cps/10^-6；K、U、Th活度浓度灵敏度的单位为cps/(Bq/kg)；

根据下式计算能谱仪地面测量结果与AP-M模型参考量值的示值误差：

式中，Er_i为无人机航空γ能谱仪的地面示值误差；q_i、Q_i为第i种元素含量的测量值和参考量值；Q_b.i为AP-B模型中第i种元素的参考量值。

可选的，所述采用示值误差检验法，对无人机航空γ能谱仪灵敏度随高度变化的修正参数和剥离系数随高度变化的修正系数进行检验，具体包括：

使用无人机航空γ能谱仪在地面半径200m，放射性分布均匀，无障碍物遮挡的陆地上空进行不同高度的悬停测量，以得到无人机航空γ能谱测量原始数据；其中，每个高度测量飞行时间不少于10分钟；

获得无人机航空γ能谱测量原始数据后，采用前述步骤的能窗剥离修正方法和高度衰减修正方法进行数据处理，得到不同飞行高度的测量结果；

采用下式计算无人机航空γ能谱仪测量结果的示值误差：

式中，Er_h.i为无人机航空γ能谱仪的示值误差；q_h.i为无人机航空γ能谱仪在h高度上对第i种元素含量的测量值；Q_i为标准配套能谱仪对地表第i种元素测定的参考量值。

可选的，在所述对本底修正后的各计数率进行能窗剥离修正的步骤中，以及在所述采用对数拟合法对能窗剥离修正后的各计数率进行高度衰减修正的步骤中，还对各计数率进行STP高度修正，以消除不同环境温度和大气压力的影响。

本申请的实施例提供的技术方案可以包括以下有益效果：

本申请的实施例提供的技术方案中，首先采集不同能量(也即多能量)γ射线的能谱数据，之后依次对能谱数据中的各计数率进行本底修正、能窗剥离修正和高度衰减修正，再采用灵敏度换算法将高度衰减修正后的各计数率统一到相同单位，从而得到分析结果。如此设置，可以提高无人机航空γ能谱仪对多能量γ射线核素的探测能力，解决地面和空中测量时，原始数据修正系数的确定问题，给出地面和空中测量的结果；还可进一步推广应用于地面和地下钻井工作中的γ能谱测量。与现有方法相比，本方法可同时分析多个γ射线能窗，测量精度更高；与全谱分析法相比，本方法可以避免一些干扰谱段对测量结果的影响。

应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本申请。

附图说明

此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本申请的实施例，并与说明书一起用于解释本申请的原理。

图1为本申请实施例提供的一种基于多能量响应的航空γ能谱分析方法的整体流程示意图；

图2为本申请实施例提供的基于多能量响应的航空γ能谱分析方法的具体实现流程示意图；

图3为基于多能量响应的航空γ能谱分析方法中的能谱剥离系数和探测灵敏度的确定流程示意图；

图4为多能量响应航空γ能谱仪对K元素的灵敏度随高度变化曲线和拟合误差曲线示意图；

图5为多能量响应航空γ能谱仪对U元素的灵敏度随高度变化曲线和拟合误差曲线示意图；

图6为多能量响应航空γ能谱仪对Th元素的灵敏度随高度变化曲线和拟合误差曲线示意图。

具体实施方式

这里将详细地对示例性实施例进行说明，其示例表示在附图中。下面的描述涉及附图时，除非另有表示，不同附图中的相同数字表示相同或相似的要素。以下示例性实施例中所描述的实施方式并不代表与本申请相一致的所有实施方式。相反，它们仅是与如所附权利要求书中所详述的、本申请的一些方面相一致的装置和方法的例子。

航空γ能谱测量是一种覆盖范围大、测量速度快、经济可靠的放射性测量技术，在我国铀矿资源勘查、辐射环境调查、核应急、核反恐、放射源搜寻中发挥着重要作用。随着无人机技术的快速发展，智能化、无人员损伤风险的无人机航空γ能谱测量技术越来越受到人们的重视。

本发明的目的在于针对小型化、轻型化的航空γ能谱测量技术，提供一种新的多能量响应的γ射线能谱分析方法，以提高探测器对发射多能量γ射线的核素的探测能力。以下通过实施例进行详细说明。

实施例

参照图1和图2，图1为本申请实施例提供的一种基于多能量响应的航空γ能谱分析方法的整体流程示意图，图2为本申请实施例提供的基于多能量响应的航空γ能谱分析方法的具体实现流程示意图。

如图1和图2所示，该方法至少包括以下步骤：

S101：使用多种放射性核素对无人机航空γ能谱仪进行测试，以采集不同能量γ射线的能谱数据；其中，所述能谱数据包括多个不同能窗的计数率；

本实施例中，放射性核素使用K系核素、U系核素和Th系核素，比如¹³⁷Cs、⁴⁰K、²⁰⁸Tl、²²⁸Ac、²¹⁴Bi等点状放射源，并且，以核工业航测遥感中心研制的UGRS-5型无人机航空γ能谱仪测量数据为例进行详细说明。该航空γ能谱仪配置1条CsI(Na)晶体探测器，晶体尺寸4″×4″×16″，能量分辨率在10％左右。

具体的，使用多种放射性核素对无人机航空γ能谱仪进行测试，以采集不同能量γ射线的能谱数据时，需要对能谱仪进行设置，具体设置过程包括：

根据测试结果中能谱仪对不同能量(例如，本实施例中采用661.7keV、1460.8keV、2614.5keV、911.2keV、583.2keV、1764.5keV、609.3keV、2204.1keV等)γ射线的特征峰峰位和能量分辨率(FWHM)，确定能谱仪的能量刻度系数和分辨率刻度系数；由能谱仪刻度系数计算能谱仪的窗宽，计算得到多能量响应的能窗位置和宽度；再基于能窗位置和宽度对能谱仪进行设置。

其中，1)针对能谱仪刻度：

①采用高斯函数对测量的能谱进行最小二乘拟合，得到每种γ射线的特征峰峰位和标准偏差：

式中，x为能谱道值；μ为特征γ射线全能峰中心峰位；σ为峰位的标准偏差；k为全能峰面积；a为线性函数斜率；b为线性函数截距。

则能量分辨率FWHM为：

②能谱仪能量刻度系数采用线性回归方程计算：

E(Ch)＝aCh+b (3)

式中，E(Ch)为γ射线特征能量值；Ch为γ射线特征峰峰位，即高斯函数计算出的μ值；a、b为能谱仪能量刻度系数。

③能谱仪分辨率刻度系数采用二次回归方程计算：

FWHM(Ch)＝aCh²+bCh+c (4)

式中，FWHM(Ch)为能谱仪对γ射线的能量分辨率；Ch为γ射线特征峰峰位，即高斯函数中的μ值；a、b、c为能量分辨率刻度系数。

2)针对能窗设置：

根据国际原子能机构(IAEA)核数据库提供的天然放射性核素特征γ射线能量，采用公式(3)的逆函数计算γ射线的特征峰峰位，即：

式中，Ch_i为第i种γ射线的特征峰峰位；E_i为第i种γ射线的特征能量值。

能谱仪各能窗的窗宽按下式计算：

LLCh_i＝Ch_i-FWHM_i (6)

RLCh_i＝Ch_i+FWHM_i (7)

式中，LLCh_i、RLCh_i分别表示第i能窗的左边界和有边界能谱道值；Ch_i为第i能窗的γ射线特征峰峰位；FWHM_i为第i能窗峰位的能量分辨率。

本实施例中，根据测试数据，航空γ能谱仪的能量刻度系数为：

E(Ch)＝11.928Ch-12.256；

分辨率刻度系数为：

FWHM(Ch)＝2.73e-5Ch²+6.38e-2Ch+1.952。

且由于采用CsI(Na)晶体的航空γ能谱仪测量天然放射性核素²¹⁴Bi、²⁰⁸Tl时，需要考虑天然放射性核素⁴⁰K的影响，因此采用如下表的参数对能谱仪进行设置，其中，K系核素对应包括K能窗，U系核素对应包括U1、U2和U3能窗，Th系核素对应包括Th1和Th2能窗。

通过上述步骤得到多个不同能窗的计数率等能谱数据后，则按照后续步骤进行分析。

S102：对各计数率进行本底修正，以扣除仪器本身、无人机平台、宇宙射线、大气氡在能谱仪采集时产生的计数率；

也即，将仪器本身、无人机平台、宇宙射线、大气氡在能谱仪探测器中产生的计数率作为一个整体，采用综合本底法对其进行扣除：

n_i＝N_i-N_b (8)

式中，N_i为无人机航空γ能谱仪对地测量飞行过程中，记录的第i能窗计数率；N_b为第i能窗的综合本底；N_i为第i能窗的本底修正结果。

进一步的，进行本底修正时，通过如下方法确定本底修正系数(也即各能窗的综合本底)：

在选定的水面上方进行不同高度的悬停测量；其中，要求水上动态测试点距离水岸1km以上，水深大于20m，每个高度测量飞行时间不少于10分钟；此外，本实施例选定的水面为黄壁庄水库，测试的不同高度包括30m、50m、70m、90m、120m等等；

根据测量的原始数据计算每个飞行高度上的能窗计数率平均值，进而计算不同飞行高度的平均值作为航空γ能谱仪本底修正系数。

根据本实施例的实际测试数据，航空γ能谱仪的本底修正系数见下表：

S103：对本底修正后的各计数率进行能窗剥离修正；

为了充分利用U系、Th系核素产生的特征γ射线计数率，更好地提取能谱数据中的有用信息，减小统计涨落和仪器噪声对测量结果的影响。本实施例在进行能窗剥离修正时，在传统能窗剥离法(也即IAEA-1363报告)的基础上，将U系和Th系核素的低能特征γ射线能窗计数率纳入计算，得到无人机飞行高度上的能窗净计数率(也即修正结果)；公式如下：

式中，n_k.k、n_u.u、n_t.t为修正后的K、U、Th能窗净计数率；n_k、n_u1、n_u2、n_u3、n_t1、n_t2为本底修正后的能窗计数率；c_k.u1、c_k.u2、c_k.u3、c_k.t1、c_k.t2为K窗对U1、U2、U3、Th1、Th2能窗的剥离系数，△c_k.u1(h_stp)、△c_k.u2(h_stp)、△c_k.u3(h_stp)、△c_k.t1(h_stp)、△c_k.t2(h_stp)为K窗剥离系数随高度的变化函数；c_u1.k、c_u1.u2、c_u1.u3、c_u1.t1、c_u1.t2为U1窗对K、U2、U3、Th1、Th2能窗的剥离系数，△c_u1.k(h_stp)、△c_u1.u2(h_stp)、△c_u1.u3(h_stp)、△c_u1.t1(h_stp)、△c_u1.t2(h_stp)为U1窗剥离系数随高度的变化函数；c_t1.k、c_t1.u1、c_t1.u2、c_t1.u3、c_t1.t2为Th1窗对K、U1、U2、U3、Th2能窗的剥离系数，△c_t1.k(h_stp)、△c_t1.u1(h_stp)、△c_t1.u2(h_stp)、△c_t1.u3(h_stp)、△c_t1.t2(h_stp)为Th1窗剥离系数随高度的变化函数；其中，c_t1.u1、c_t1.k、c_u1.k与传统剥离法中的α、β、γ一致，c_k.u1、c_k.t1、c_u1.t1与传统剥离法中的a、b、g一致。

其中，如图3所示，剥离系数(也即c_k.u1等)和灵敏度的计算过程包括：

将无人机航空γ能谱仪置于不同的航空放射性测量模型标准装置(简称：航放模型)上，分别对本底模型AP-B、钾模型AP-K、铀模型AP-U、钍模型AP-Th、混合模型AP-M进行测试；其中，每个模型测量2次，每次测量10分钟；同时，参照JJG(军工)26-2012航空γ能谱仪检定规程，使用标准配套能谱仪对航空放射性测量模型标准装置进行参考量值测定；

根据测试结果，按照下式计算无人机航空γ能谱仪在地面的剥离系数和灵敏度：

式中，n_k.k、n_k.u1、n_k.u2、n_k.t1、n_k.t2为探测器在AP-K模型上测到的能窗计数率扣除AP-B模型上计数率后的K、U1、U2、U3、Th1、Th2窗计数率结果，n_u.k、n_u.u1、n_u.u2、n_u.u3、n_u.t1、n_u.t2、为探测器在AP-U模型上测到的能窗计数率扣除AP-B模型上计数率后的窗计数率结果，n_t.k、n_t.u1、n_t.u2、n_t.u3、n_t.t1、n_t.t2为探测器在AP-T模型上测到的能窗计数率扣除AP-B模型上计数率后的窗计数率结果；

q_k.k、q_k.u、q_k.t为使用标准配套能谱仪定值测量得到的AP-K模型中的K、U、Th含量扣除AP-B模型中K、U、Th含量的量值，q_u.k、q_u.u、q_u.t为AP-U模型中的K、U、Th含量扣除AP-B模型中K、U、Th含量的量值，q_t.k、q_t.u、q_t.t为AP-Th模型中的K、U、Th含量扣除AP-B模型中K、U、Th含量的量值；

s_k为探测器的K窗对K元素的探测灵敏度，s_u为探测器的U1窗对U元素的探测灵敏度，s_t为探测器的Th1窗对Th元素的探测灵敏度；

c_k.u1、c_k.u2、c_k.u3、c_k.t1、c_k.t2为U1、U2、U3、Th1、Th2窗对K元素的灵敏度与K窗对K元素灵敏度的比值，c_u1.k、c_u1.u2、c_u1.u3、c_u1.t1、c_u1.t2为K、U2、U3、Th1、Th2窗对U元素的灵敏度与U1窗对U元素灵敏度的比值，c_t1.k、c_t1.u1、c_t1.u2、c_t1.u3、c_t1.t2为K、U1、U2、U3、Th2窗对Th元素的灵敏度与Th1窗对Th元素灵敏度的比值，这些参数亦称为剥离系数。

根据本实施例中航空放射性测量模型标准装置静态测试数据，无人机航空γ能谱仪在地面的剥离系数和灵敏度计算结果见下表：

S104：采用对数拟合法对能窗剥离修正后的各计数率进行高度衰减修正；

具体的，前述步骤测得的是无人机航空γ能谱仪在地面的剥离系数和灵敏度结果，但实际上，剥离系数和灵敏度是随高度衰减的，因此，为了保证最终结果的精确性，还需要对前述步骤的能窗剥离修正结果进行高度衰减修正。

更具体的，传统的高度衰减修正采用能窗计数率与飞行高度呈指数衰减的规律，一般采用下式计算归一化高度(h₀)上的能窗净计数率：

式中，h_stp为测量飞行离地高度换算到STP状态时的高度值；n_0.i为第i能窗在STP高度为h₀时的净计数率。对于有人机航空γ能谱测量而言，一般取h₀＝120m；n_h.i为第i能窗在h_stp高度上的净计数率；μ_i为第i能窗的高度衰减系数。

参照EJ/T 1032-2005航空伽玛能谱测量规范，上述方法一般采用60～300m高度的实测数据进行拟合。根据蒙卡模拟计算结果，若飞行超出该高度范围，拟合误差将超过±10％。而对于本实施例采用的新型的无人机平台，航空γ能谱测量的飞行高度一般可以控制50m左右，上述方法的修正误差较大。因此，本实施例采用对数拟合法进行高度衰减修正，即：

式中，a、b、c、d为多项式拟合系数，令：

G(h)＝ah³+bh²+ch (13)

则有

ln(100×n_0.i)＝ln(100×n_h.i)+G(h₀)-G(h_stp) (14)

根据上述方法对能窗剥离修正后的n_k.k、n_u.u、n_t.t进行高度衰减修正，可得到归一化高度(h₀)上的计数率。之后还可以根据蒙卡模拟计算确定高度衰减修正系数的拟合误差，验证拟合效果。其中，对高度衰减修正后的结果进行蒙卡模拟计算的具体过程包括：

①根据无人机航空γ能谱仪的几何结构，以及各结构部件制作材料的元素成分、密度，建立蒙特卡罗模拟数学模型；其中，各部件制作材料的元素成分可采用全元素取样分析法对材料进行取样分析，几何结构可采用成像扫描和数字化成图的方式提高数学建模的准确性；

②根据航空放射性测量模型标准装置的物质成分、密度，建立一个半径大于2.5km，厚度在50cm以上的饱和体源数学模型；

③将饱和体源数学模型设定为核素分布均匀的模拟源，其中，模拟源发射的γ射线能量、强度及子体核素分支比等参数参照国际原子能机构(IAEA)核数据库提供的U系、Th系子体核素，以及⁴⁰K核素；

④采用方差减小技术对无人机航空γ能谱仪和饱和体源组成的数学模型进行优化，在大型科学计算平台上模拟计算无人机航空γ能谱仪在饱和体源模型中心点正上方0.135、1、5、10、20、30、60、90、120、150、180、210、240、270、300m高度时的模拟能谱数据；其中，0.135m为航空放射性测量模型标准装置静态测试时探测器的离地高度；

⑤根据无人机航空γ能谱仪的实际刻度和能窗设置参数对模拟能谱数据进行高斯展宽，进而计算K、U1、U2、U3、Th1、Th2窗计数率，并根据在航空放射性测量模型标准装置上测得的灵敏度，对0～300m高度上的模拟数据进行归一化；也即：

式中，s_h.i.j为模拟高度为h时，模拟第i元素在j窗的灵敏度；s_0.i.j为航空放射性测量模型标准装置静态测试时，实测第i元素在j窗的灵敏度；n_0.i.j为能谱仪在饱和体源上0.135m高度时，模拟第i元素在j窗的计数率；n_h.i.j为能谱仪在饱和体源上h高度时，模拟第i元素在j窗的计数率。

⑥按下式计算能谱仪在饱和体源上h高度时的剥离系数：

式中，s_h.k.j、s_h.u.j、s_h.t.j为第j窗对K、U、Th元素的灵敏度；s_h.k.k为K窗对K元素的灵敏度，s_h.u.u1为U1窗对U元素的灵敏度，s_h.t.t1为Th1窗对Th元素的灵敏度；c_h.k.j、c_h.u.j、c_h.t.j为第j窗对K、U、Th元素的剥离系数；

⑦由于灵敏度与能窗计数率成正比(参照下文公式(19))，则灵敏度随高度的变化规律与航空γ能谱仪能窗计数率相同，因此按下式计算无人机航空γ能谱仪灵敏度随高度变化的修正参数：

式中，s_i.j(h_stp)为能谱仪在饱和体源模型h高度上，第j窗对第i元素的灵敏度；a_i.j、b_i.j、c_i.j、d_i.j为第j窗对第i元素的灵敏度修正参数；

本实施例中，无人机航空γ能谱仪灵敏度随高度变化的修正参数见下表：

⑧参照IAEA-1363技术报告，根据不同高度的模拟计算结果，采用线性回归拟合得到剥离系数随高度变化的修正参数。

本实施例中，无人机航空γ能谱仪剥离系数随高度变化的修正参数见下表：

根据蒙卡模拟计算结果，在0～300m高度上高度衰减修正系数的拟合误差在±5％以内，明显好于指数衰减拟合函数的效果(参见图4～6)。

S105：采用灵敏度换算法将高度衰减修正后的各计数率统一到含量或活度浓度量值单位，得到分析结果。

也即，将不同飞行高度上的能窗计数率归一到h₀高度上以后，采用灵敏度换算法，将计数率统一到含量或活度浓度量值单位：

式中，Q_i为第i种核素(K、U、Th)的含量或活度浓度量值，对于天然放射性而言，K含量的单位一般为％、U为10^-6、Th为10^-6；对于活度浓度量值而言，K、U、Th的单位一般为Bq/kg；s_i为第i种核素的灵敏度，对于天然放射性核素而言，K含量灵敏度的单位一般为cps/％、U为cps/10^-6、Th为cps/10^-6；K、U、Th活度浓度灵敏度的单位为cps/(Bq/kg)。

本申请的实施例提供的技术方案中，首先采集不同能量(也即多能量)γ射线的能谱数据，之后依次对能谱数据中的各计数率进行本底修正、能窗剥离修正和高度衰减修正，再采用灵敏度换算法将高度衰减修正后的各计数率统一到相同单位，从而得到分析结果。如此设置，可以提高无人机航空γ能谱仪对多能量γ射线核素的探测能力，解决地面和空中测量时，原始数据修正系数的确定问题，给出地面和空中测量的结果；还可进一步推广应用于地面和地下钻井工作中的γ能谱测量。与现有方法相比，本方法可同时分析多个γ射线能窗，测量精度更高；与全谱分析法相比，本方法可以避免一些干扰谱段对测量结果的影响。

此外，由于无人机航空γ能谱测量过程中，大气层的密度会受环境温度和大气压力影响，造成探测器对地面核素的响应发生变化。因此步骤S104和S105中，应对飞行离地高度进行标准温度和大气压高度(STP)计算，从而得到更加稳定的STP离地高度。也即，对各计数率进行STP高度修正，以消除不同环境温度和大气压力的影响，参照如下公式：

式中，h_stp为环境温度0℃(即273.15K)、大气压力101.325kPa状态下，探测器的离地高度。h_obs为实际测量环境中，探测器的离地高度值；T为实际测量环境中的大气温度值；P为测量环境中的大气压力值。

其中，当采用海拔高度表示大气压力环境时，需要进行大气压力换算，即：

式中，H为实际测量环境的海拔高度；R为摩尔气体常量，取8314.472(m³ kPa)/(mol K)；T为热力学温度；g为重力加速度，通常为9.8m/s²；M为气体的摩尔质量，干燥空气的摩尔质量为28.9634g/mol；P₀为标准大气压，取为101.325kPa；P为所在海拔高度上的气压值。

此外，为了保证分析结果的准确性，还可以采用示值误差检验法，对无人机航空γ能谱仪在地面的灵敏度和剥离系数进行检验；以及，采用示值误差检验法，对无人机航空γ能谱仪灵敏度随高度变化的修正参数和剥离系数随高度变化的修正系数进行检验。

其中，对无人机航空γ能谱仪在地面的灵敏度和剥离系数进行检验的具体过程包括：

根据无人机航空γ能谱仪在混合模型AP-M上的测量数据，采用如下公式计算模型含量或活度浓度值；

根据无人机航空γ能谱仪在混合模型AP-M上的测量数据，采用如下公式以及公式(16)计算模型含量或活度浓度值；

式中，n_m.k、n_m.u1、n_m.u2、n_m.t1、n_m.t2为探测器在AP-M模型上测到的能窗计数率扣除AP-B模型上计数率后的K、U1、U2、U3、Th1、Th2窗计数率；

根据下式计算能谱仪地面测量结果与AP-M模型参考量值的示值误差：

式中，Er_i为无人机航空γ能谱仪的地面示值误差；q_i、Q_i为第i种元素含量的测量值和参考量值；Q_b.i为AP-B模型中第i种元素的参考量值。

参照JJG(军工)26-2012航空γ能谱仪检定规程的要求，地面示值误差应控制在±5％以内

此外，对无人机航空γ能谱仪灵敏度随高度变化的修正参数和剥离系数随高度变化的修正系数进行检验的具体过程包括：

使用无人机航空γ能谱仪在地面半径200m，放射性分布均匀，无障碍物遮挡的陆地上空进行不同高度的悬停测量，以得到无人机航空γ能谱测量原始数据；其中，每个高度测量飞行时间不少于10分钟；

获得无人机航空γ能谱测量原始数据后，采用前述步骤的能窗剥离修正方法和高度衰减修正方法进行数据处理，得到不同飞行高度的测量结果；

采用下式计算无人机航空γ能谱仪测量结果的示值误差：

式中，Er_h.i为无人机航空γ能谱仪的示值误差；q_h.i为无人机航空γ能谱仪在h高度上对第i种元素含量的测量值；Q_i为标准配套能谱仪对地表第i种元素测定的参考量值。

参照JJG(军工)26-2012航空γ能谱仪检定规程的要求，无人机航空γ能谱测量在不同高度上测量结果的示值误差应控制在±30％以内。

可以理解的是，上述各实施例中相同或相似部分可以相互参考，在一些实施例中未详细说明的内容可以参见其他实施例中相同或相似的内容。

需要说明的是，在本申请的描述中，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性。此外，在本申请的描述中，除非另有说明，“多个”的含义是指至少两个。

流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本申请的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本申请的实施例所属技术领域的技术人员所理解。

应当理解，本申请的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。

本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。

此外，在本申请各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。

上述提到的存储介质可以是只读存储器，磁盘或光盘等。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本申请的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管上面已经示出和描述了本申请的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本申请的限制，本领域的普通技术人员在本申请的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (10)

1.一种基于多能量响应的航空γ能谱分析方法，其特征在于，包括：

使用多种放射性核素对无人机航空γ能谱仪进行测试，以采集不同能量γ射线的能谱数据；其中，所述能谱数据包括多个不同能窗的计数率；

对各计数率进行本底修正，以扣除仪器本身、无人机平台、宇宙射线、大气氡在能谱仪采集时产生的计数率；

对本底修正后的各计数率进行能窗剥离修正；

采用对数拟合法对能窗剥离修正后的各计数率进行高度衰减修正；

采用灵敏度换算法将高度衰减修正后的各计数率统一到含量或活度浓度量值单位，得到分析结果。

2.根据权利要求1所述的方法，其特征在于，所述使用多种放射性核素对无人机航空γ能谱仪进行测试，以采集不同能量γ射线的能谱数据，包括：

对能谱仪进行设置，具体包括：

根据测试结果中能谱仪对不同能量γ射线的特征峰峰位和能量分辨率，确定能谱仪的能量刻度系数和分辨率刻度系数；

由所述能谱仪刻度系数计算能谱仪的窗宽，计算得到多能量响应的能窗位置和宽度；

基于所述能窗位置和宽度对能谱仪进行设置。

3.根据权利要求1所述的方法，其特征在于，所述多种放射性核素包括K系核素、U系核素和Th系核素；K系核素对应包括K能窗，U系核素对应包括U1、U2和U3能窗，Th系核素对应包括Th1和Th2能窗；

进行能窗剥离修正时，在传统能窗剥离法的基础上，将U系和Th系核素的低能特征γ射线能窗计数率纳入计算，得到无人机飞行高度上的能窗净计数率；公式如下：

式中，n_k.k、n_u.u、n_t.t为修正后的K、U、Th能窗净计数率；n_k、n_u1、n_u2、n_u3、n_t1、n_t2为本底修正后的能窗计数率；c_k.u1、c_k.u2、c_k.u3、c_k.t1、c_k.t2为K窗对U1、U2、U3、Th1、Th2能窗的剥离系数，△c_k.u1(h_stp)、△c_k.u2(h_stp)、△c_k.u3(h_stp)、△c_k.t1(h_stp)、△c_k.t2(h_stp)为K窗剥离系数随高度的变化函数；c_u1.k、c_u1.u2、c_u1.u3、c _u1.t1、c_u1.t2为U1窗对K、U2、U3、Th1、Th2能窗的剥离系数，△c_u1.k(h_stp)、△c_u1.u2(h_stp)、△c_u1.u3(h_stp)、△c _u1.t1(h_stp)、△c _u1.t2(h_stp)为U1窗剥离系数随高度的变化函数；c_t1.k、c_t1.u1、c_t1.u2、c_t1.u3、c_t1.t2为Th1窗对K、U1、U2、U3、Th2能窗的剥离系数，△c_t1.k(h_stp)、△c_t1.u1(h_stp)、△c_t1.u2(h_stp)、△c_t1.u3(h_stp)、△c_t1.t2(h_stp)为Th1窗剥离系数随高度的变化函数；其中，c_t1.u1、c_t1.k、c_u1.k与传统剥离法中的α、β、γ一致，c_k.u1、c_k.t1、c_u1.t1与传统剥离法中的a、b、g一致。

4.根据权利要求3所述的方法，其特征在于，所述剥离系数的计算过程包括：

将无人机航空γ能谱仪置于不同的航空放射性测量模型标准装置上，分别对本底模型AP-B、钾模型AP-K、铀模型AP-U、钍模型AP-Th、混合模型AP-M进行测试；

根据测试结果，按照下式计算无人机航空γ能谱仪在地面的剥离系数和灵敏度：

式中，n_k.k、n_k.u1、n_k.u2、n_k.t1、n_k.t2为探测器在AP-K模型上测到的能窗计数率扣除AP-B模型上计数率后的K、U1、U2、U3、Th1、Th2窗计数率结果，n_u.k、n_u.u1、n_u.u2、n_u.u3、n_u.t1、n_u.t2、为探测器在AP-U模型上测到的能窗计数率扣除AP-B模型上计数率后的窗计数率结果，n_t.k、n_t.u1、n_t.u2、n_t.u3、n_t.t1、n_t.t2为探测器在AP-T模型上测到的能窗计数率扣除AP-B模型上计数率后的窗计数率结果；

q_k.k、q_k.u、q_k.t为使用标准配套能谱仪定值测量得到的AP-K模型中的K、U、Th含量扣除AP-B模型中K、U、Th含量的量值，q_u.k、q_u.u、q_u.t为AP-U模型中的K、U、Th含量扣除AP-B模型中K、U、Th含量的量值，q_t.k、q_t.u、q_t.t为AP-Th模型中的K、U、Th含量扣除AP-B模型中K、U、Th含量的量值；

s_k为探测器的K窗对K元素的探测灵敏度，s_u为探测器的U1窗对U元素的探测灵敏度，s_t为探测器的Th1窗对Th元素的探测灵敏度；

c_k.u1、c_k.u2、c_k.u3、c_k.t1、c_k.t2为U1、U2、U3、Th1、Th2窗对K元素的灵敏度与K窗对K元素灵敏度的比值，c_u1.k、c_u1.u2、c_u1.u3、c_u1.t1、c_u1.t2为K、U2、U3、Th1、Th2窗对U元素的灵敏度与U1窗对U元素灵敏度的比值，c_t1.k、c_t1.u1、c_t1.u2、c_t1.u3、c_t1.t2为K、U1、U2、U3、Th2窗对Th元素的灵敏度与Th1窗对Th元素灵敏度的比值，这些参数亦称为剥离系数。

5.根据权利要求4所述的方法，其特征在于，所述采用对数拟合法对能窗剥离修正后的各计数率进行高度衰减修正，包括：

采用如下公式计算高度衰减修正系数：

式中，a、b、c、d为多项式拟合系数，令：

G(h)＝ah³+bh²+ch

则有

ln(100×n_0.i)＝ln(100×n_h.i)+G(h₀)-G(h_stp)。

6.根据权利要求5所述的方法，其特征在于，所述采用对数拟合法对能窗剥离修正后的各计数率进行高度衰减修正，之后还包括：

根据蒙卡模拟计算确定高度衰减修正系数的拟合误差，具体包括：

根据无人机航空γ能谱仪的几何结构，以及各结构部件制作材料的元素成分、密度，建立蒙特卡罗模拟数学模型；

根据航空放射性测量模型标准装置的物质成分、密度，建立一个半径大于2.5km，厚度在50cm以上的饱和体源数学模型；

将饱和体源数学模型设定为核素分布均匀的模拟源，其中，模拟源发射的γ射线能量、强度及子体核素分支比等参数参照国际原子能机构核数据库提供的U系、Th系子体核素，以及⁴⁰K核素；

采用方差减小技术对无人机航空γ能谱仪和饱和体源组成的数学模型进行优化，在大型科学计算平台上模拟计算无人机航空γ能谱仪在饱和体源模型中心点正上方0.135、1、5、10、20、30、60、90、120、150、180、210、240、270、300m高度时的模拟能谱数据；其中，0.135m为航空放射性测量模型标准装置静态测试时探测器的离地高度；

根据无人机航空γ能谱仪的实际刻度和能窗设置参数对模拟能谱数据进行高斯展宽，进而计算K、U1、U2、U3、Th1、Th2窗计数率，并根据在航空放射性测量模型标准装置上测得的灵敏度，对0～300m高度上的模拟数据进行归一化；

按下式计算能谱仪在饱和体源上h高度时的剥离系数：

式中，s_h.k.j、s_h.u.j、s_h.t.j为第j窗对K、U、Th元素的灵敏度；s_h.k.k为K窗对K元素的灵敏度，s_h.u.u1为U1窗对U元素的灵敏度，s_h.t.t1为Th1窗对Th元素的灵敏度；c_h.k.j、c_h.u.j、c_h.t.j为第j窗对K、U、Th元素的剥离系数；

按下式计算无人机航空γ能谱仪灵敏度随高度变化的修正参数：

式中，s_i.j(h_stp)为能谱仪在饱和体源模型h高度上，第j窗对第i元素的灵敏度；a_i.j、b_i.j、c_i.j、d_i.j为第j窗对第i元素的灵敏度修正参数；

参照IAEA-1363技术报告，根据不同高度的模拟计算结果，采用线性回归拟合得到剥离系数随高度变化的修正参数。

7.根据权利要求6所述的方法，其特征在于，还包括：

采用示值误差检验法，对无人机航空γ能谱仪在地面的灵敏度和剥离系数进行检验；

以及，

采用示值误差检验法，对无人机航空γ能谱仪灵敏度随高度变化的修正参数和剥离系数随高度变化的修正系数进行检验。

8.根据权利要求7所述的方法，其特征在于，所述采用示值误差检验法，对无人机航空γ能谱仪在地面的灵敏度和剥离系数进行检验，具体包括：

根据无人机航空γ能谱仪在混合模型AP-M上的测量数据，采用如下公式计算模型含量或活度浓度值；

式中，n_m.k、n_m.u1、n_m.u2、n_m.t1、n_m.t2为探测器在AP-M模型上测到的能窗计数率扣除AP-B模型上计数率后的K、U1、U2、U3、Th1、Th2窗计数率；Q_i为第i种核素的含量或活度浓度量值，对于天然放射性而言，K含量的单位一般为％、U为10^-6、Th为10^-6；对于活度浓度量值而言，K、U、Th的单位一般为Bq/kg；s_i为第i种核素的灵敏度，对于天然放射性核素而言，K含量灵敏度的单位一般为cps/％、U为cps/10^-6、Th为cps/10^-6；K、U、Th活度浓度灵敏度的单位为cps/(Bq/kg)；

根据下式计算能谱仪地面测量结果与AP-M模型参考量值的示值误差：

式中，Er_i为无人机航空γ能谱仪的地面示值误差；q_i、Q_i为第i种元素含量的测量值和参考量值；Q_b.i为AP-B模型中第i种元素的参考量值。

9.根据权利要求7所述的方法，其特征在于，所述采用示值误差检验法，对无人机航空γ能谱仪灵敏度随高度变化的修正参数和剥离系数随高度变化的修正系数进行检验，具体包括：

使用无人机航空γ能谱仪在地面半径200m，放射性分布均匀，无障碍物遮挡的陆地上空进行不同高度的悬停测量，以得到无人机航空γ能谱测量原始数据；其中，每个高度测量飞行时间不少于10分钟；

获得无人机航空γ能谱测量原始数据后，采用前述步骤的能窗剥离修正方法和高度衰减修正方法进行数据处理，得到不同飞行高度的测量结果；

采用下式计算无人机航空γ能谱仪测量结果的示值误差：

式中，Er_h.i为无人机航空γ能谱仪的示值误差；q_h.i为无人机航空γ能谱仪在h高度上对第i种元素含量的测量值；Q_i为标准配套能谱仪对地表第i种元素测定的参考量值。

10.根据权利要求1-9任一项所述的方法，其特征在于，在所述对本底修正后的各计数率进行能窗剥离修正的步骤中，以及在所述采用对数拟合法对能窗剥离修正后的各计数率进行高度衰减修正的步骤中，还对各计数率进行STP高度修正，以消除不同环境温度和大气压力的影响。

Images