基于Python 编程的航空伽马全谱分析方法的应用研究

2022-09-01冯春圆王会波周子阳蔡文军杨春雷朱林

世界核地质科学 2022年2期

冯春圆，王会波，周子阳，蔡文军，杨春雷，朱林

（1.核工业航测遥感中心河北省航空探测与遥感技术重点实验室中核铀资源地球物理勘查技术中心（重点实验室），河北石家庄 050002；2.中核地矿科技集团有限公司，北京100013；3.华业（北京）科技有限公司，北京100013；4.中国冶金地质总局一局，河北三河 065201）

传统的航空伽马能谱数据处理方法主要为国际原子能机构（1991）［1］推荐的“标准三窗法”，简称“三窗法”（TWM），即所谓的窗数据处理法，这种方法在地表勘查和浅部铀矿床勘探过程中曾发挥了重要作用，而且仍在地质矿产勘查等许多领域得到较为广泛的应用。该方法也是众多数据处理软件常用的数据处理方法，如oasis montaj、AGRS GeoProbe 等软件，均提供标准三窗法数据处理模块。随着航空伽马能谱测量向放射性本底调查、环境测量、核事故应急航空伽马能谱测量、地下水及温泉调查以及大气氡评价和修正等方面的应用［2-9］。数据处理方法也从单一的窗数据处理走向谱数据处理发展。

航空伽马全能谱分析方法，简称“全谱法”（FSA），是基于全能谱数据进行的一种数据处理方法。就是利用航空伽马能谱测量中所有能谱数据，而不是常规的窗数据，直接进行飞机和仪器本底、宇宙辐射本底、航测高度等的各项修正，而将大气氡视为变量，与地面辐射源的钾、铀、钍元素含量以及其他核污染源（人工核素）的浓度或活度一起进行计算的分析方法［10］。该方法能有效解决窗数据处理中大气氡修正引入较大不确定度问题，前人已在原理和应用方面作了大量的研究［11-14］。核工业航测遥感中心于2003年在Windows XP 操作系统上，基于Microsoft.NET Framework 模块通过C#调用MATLAB软件进行了针对航空伽马全能谱分析方法的编程。受限于MATLAB 与Microsoft.NET Framework的通信、MATLAB与C#通信过程中的数据转换等影响，以致软件整体的运行速度较慢、可移植性和兼容性较低。

Python 是一种面向对象的解释性的跨平台高级计算机语言，具有语法简单、可移植性强的特点。同时，它是一种通用性语言并具有极为丰富的类库，可应用于数据库、多媒体、科学计算、网络、游戏等诸多领域。

作者从航空伽马全能谱分析方法的原理出发，以计算天然放射性核素为例，利用Python 的科学计算库，实现了航空伽马全能谱分析方法计算天然放射性核素。通过对云南昆明某试验区航空伽马能谱测量数据的处理，发现利用Python 编程处理数据，全过程无需调用MATLAB、SCILAB 等大型数值计算软件，即可实现了航空伽马全能谱分析方法的全过程，数据处理方法可行、结果可靠，处理效率大幅提升，为航空伽马全能谱分析方法计算以及其他全谱数据处理提供了一种可实现的方式。

1 航空伽马全能谱分析方法原理

对于以天然放射性核素为目的的航空伽马能谱测量，测量源项包括：宇宙射线、飞机本底、大气氡及陆地天然放射性核素40K、铀系和钍系等。测量能谱的计数率可由下式表示：

式中：N—一个（1×256）观测能谱计数率数据矩阵；B—一个（1×256）飞机和仪器辐射本底矩阵；C—一个（1×256）宇宙射线本底矩阵；Q—一个（1×4）放射性核素（天然放射性K、U、Th 及大气氡Rn）的含量（或质量活度）矩阵；M—一个（256×4）单位能谱矩阵（又称灵敏度矩阵），MT为M的转置矩阵，M-T为MT的逆矩阵。其中，计数率数据矩阵N为实际测量值，属于已知项，飞机和仪器辐射本底矩阵B、宇宙射线本底矩阵C、单位能谱矩阵U为修正系数，需在实际飞行前通过高高度测量、动态带飞行测量、航模模型刻度等获取，各测量要求详见航空伽马能谱测量规范［15］，矩阵Q为天然放射性K、U、Th 及大气氡（Rn）含量或质量活度的待求矩阵。

航空伽马全能谱分析方法是基于全谱来分析计算，在谱分析之前须检查实测谱数据和雷达、温度等数据，以确保用于全能谱分析方法计算的准确性。主要的检验内容有数据检验、峰漂修正等。数据检验包括雷达、温度、气压等数据是否完整、正常，谱数据是否有峰漂现象等；峰漂修正是针对有峰位漂移的谱进行的一种修正，使得全谱信息更符合实际测量情况。

2 航空伽马全能谱分析的程序设计

根据航空伽马全能谱分析方法的原理，程序实现过程分为数据读取与整理、预处理、拟合计算、结果输出4 个部分（图1）。

图1 航空伽马全谱分析法Python 编程设计简图Fig.1 The design diagram of Python programming for airborne gamma-ray full spectrum analysis

数据读取与整理包括本底谱数据、单位谱（K、U、Th、Rn）数据以及实测谱数据的读取，并对实测谱数据进行初步数据检查。检查内容主要包括雷达、温度、气压数据等是否齐全和可靠，如有跳点、假值等情况时需进行数据剔除和插值，尤其是雷达数据的连续性较差时，可采用5 点平滑滤波等方法进行处理。预处理主要包括死时间修正、能谱滤波、峰漂修正、计算STP 等效高度、本底扣除等，最后将飞机仪器本底谱、宇宙射线本底谱和单位谱数据分别整理成矩阵B、C、M（K、U、Th、Rn 单位谱及修正系数）。其中死时间修正、能谱滤波、峰漂修正、计算STP 等效高度均是对实测谱数据做的处理。拟合计算将净计数率谱矩阵N 按矩阵表达式（3）进行拟合即可得到K、U、Th、Rn 结果。结果输出根据数据读取的内容，将坐标信息（x，y）、雷达高度、K、U、Th、Rn 等重新整理匹配后输出。

3 航空伽马全能谱分析的编程

本次Python 编程主要引用了Csv（csv 数据格式）、NumPy（数组、矩阵计算）、Pandas（数据整理）、SciPy（工程科学计算库）等第三方类库。

3.1 数据读取及整理

航空伽马全能谱分析方法计算所用数据分为两类，一类为实测谱数据，一般包含了测线、基点号、测量日期、时间、测量活时间、温度、雷达高度、气压高度、K、U、Th、cos 窗数据以及0～255 道下测谱数据等；其次为参数文件，包括本底谱数据和K、U、Th、Rn 单位谱数据等，本次以“csv”数据格式为例。宇宙射线系数与飞机仪器本底参数合并为一个文件，包含道址（Char）、宇宙射线系数（cos_xs）、飞机仪器本底（air_bd）；单位谱数据包括道址（Char）、K 单位谱（K_COUNT）和衰减系数（Kb）、U 单位谱（U_COUNT）和衰减系数（Ub）、Th 单位谱（TH_COUNT）和衰减系数（THb）；大气氡单位谱文件包括道址（Char）、大气氡单位谱（Rn_unit）和衰减系数（Rnb）。

3.2 数据预处理

预处理包括死时间修正、能谱滤波、峰漂修正、计算STP 等效高度、宇宙射线本底计算和本底扣除以及STP 等效高度下单位谱计算等内容。

1）死时间修正：利用每道计数率与活时间对每道进行修正。

2）能谱滤波：由航空伽马全能谱分析方法的原理可知，该方法是基于全谱进行数据处理，那么谱线的噪声水平如何，对于最终的计算结果有很大的影响。国内外学者对能谱滤波已有较多的研究，如多项式拟合法、最小二乘法、小波滤波、Kalman（卡尔曼滤波）及NASVD（奇异值分解）方法等。其中小波滤波、Kalman（卡尔曼滤波方法）及NASVD 方法是目前较为流行的能谱降噪技术［16-22］，本次谱线滤波选择了NASVD 方法，基于前人的研究结果对原始谱线性滤波处理，通过程序编程完成了滤波处理，核心部分编程如下：

其中，np.linalg.svd 为奇异值分解函数，np.sqrt（）为平方根计算函数，np.dot（）为矩阵的乘积函数，diag（）为创建对角矩阵函数，A_NA 为利用实测谱数据计算所得的单位方差矩阵，U_p 为奇异值分解产生的m×m 阶酉矩阵，T_p 为奇异值分解产生的m×n 阶对角矩阵，V_pT 是奇异值分解产生的n×n阶酉矩阵；St_unit为实测谱数据归一化行向量，V_p 为V_pT 的转置矩阵。

3）峰漂修正：能谱峰漂一直是航空伽马谱测量中最为关键的一项指标，本次将能谱峰漂大于±0.5 道的谱均进行修正。此次修正方法基于前人关于峰漂修正的理论编程完成［23-24］。

4）STP 等效高度计算：当读入一个实测谱数据后，利用温度、雷达和气压数据根据公式（4），计算该测点对应的STP 等效高度。

式中：Hstp—某一测点对应的STP 等效高度，m；Hr—实测雷达高度，m；T是温度，℃；Hb是气压高度，m。

5）计算宇宙射线本底：根据实测宇宙射线窗计数与宇宙射线本底修正系数计算宇宙射线本底谱。计算公式如下：

式中：Cos—某一测点对应宇宙射线本底谱；cos—某一测点实测宇宙射线窗计数率，cos_xs—经过高高度飞行得到的宇宙射线修正系数谱。

6）本底扣除：将实测谱数据减掉宙射线本底C和飞机仪器本底B，最后得到净计数率谱。

7）STP 等效高度下单位谱计算：根据IAEA 推荐的公式（6），将地面1 m 处的单位谱（K、U、Th）根据其随高度的变化系数换算至任一STP 等效高度上，具体公式如下：

式中：Nℎ—经换算修正至任一STP 等效高度后的单位谱；N0—地面上1 m 处的单位计数率谱；μ—K、U、Th、Rn 的单位谱随高度变化的衰减系数谱；h—某一测点的STP 等效高度，m。

3.3 拟合计算

将扣除本底后的净计数率谱数据与STP 高度下的单位谱进行矩阵拟合计算，即可得到与K、U、Th、Rn 元素含量的结果矩阵，具体拟合函数定义如下：

其中，spe 为实测净计数率的谱数据（只有下测0-255 道）；M 为某测点STP 等效高度上的单位谱矩阵；np.dot（）为求取两个矩阵的乘积函数；.T 为求取某矩阵的转置矩阵；np.linalg.inv（）为求取某矩阵的逆矩阵；rs为结果矩阵；包含K、U、Th、Rn 4个元素的含量。

最后从读取的数据中整理提取x坐标、y坐标、雷达高度、窗数据等，从结果矩阵中提取K、U、Th、Rn 元素含量形成新数组，写入到结果文件中。

4 应用效果分析

根据程序设计进行了软件编程，其中核心代码不到100 行，即完成了数据预处理和拟合计算两大部分。本次选择云南昆明-昆阳磷矿试验区进行了数据处理。

4.1 程序运行效果分析

试验区长、宽约30 km，线距为1 km，共计2.3万余组能谱数据（256 道）。在i5CPU、8GBRAM 配置的Win 8系统中，该程序耗时约1 min 30 s左右；在相同配置的Win7、Win10 系统中，数据处理耗时误差均不超过15 s，并在各系统上运行稳定，无退出、死机等情况。相比以往通过调用MATLAB、SCILAB 等数值计算软件来拟合计算的编程，该程序运行效率有大幅提升，且兼容性较好。

该程序在不同版本的windonws 系统上安装简单、便捷，无需单独安装不同版本的Microsoft.NET Framework 模块以及MATLAB 或SCILAB 等大型软件，仅安装开源的Anaconda 编程开发工具即可运行程序，具有较好的移植性。

4.2 数据处理效果分析

将本程序计算的全谱法钾、铀、钍含量与传统三窗法计算的钾、铀、钍含量对比分析，结果如见表1 和图2。总体上，两种方法的钾、铀、钍分布特征及变化规律基本一致，钾元素分布特征最为接近，钍次之，铀元素变化略大。如在县街镇西部以及海口镇以西一带，三窗法航放铀含量偏高晕、高值晕规模较大、连续性较好，全谱法航放铀含量连续性较差；其次在安宁市南部7 km 处，三窗法航放铀含量呈南北向带状分布的偏高晕，而全谱法铀含量则为偏低值特征。结合飞行高度数据分析，三窗法的航放铀含量偏高是由航测飞行超高引起。

表1 试验区两种数据处理方法结果统计Table 1 Statistics of processed result of test area by two processing methods

图2 试验区两种数据处理结果对比Fig.2 Result comparison of two processing methods in the test area

为进一步验证本程序数据处理结果的可靠性以及与传统三窗法的差别，选取了L2033 线中南段约14 km 测线和L1891 测线数据进行了对比分析。两测线全谱法与三窗法计算的钾、铀、钍含量结果如表2所示，图3、图4分别为L2033测线和L1891 测线钾、铀、钍含量结果对比图。

L2033 线穿过昆阳磷矿区（164），钾和钍含量变化趋势一致（图3）；全谱法铀含量整体略大于三窗法铀含量，并在最小值上体现较为明显，但均值相差较小，相对误差仅为6.1%（表2）。

图3 L2033 线三窗法（TWM）与全谱法（FSA）计算结果对比Fig.3 Results comparison of L2033 Calculated by TWM and FSA

L1891 线穿过了观音山磷矿区（114），两种方法计算的结果基本一致（图4），仅在171-203测点（灰线区间）内，钾、铀含量结果变化较大，钾含量最大值相差可达44.4%，铀含量最大值相差为25.9%，且全谱法计算结果均低于三窗法（表2）。经分析该区段为陡峭沟谷，飞行高度均大于140 m，局部高达425.8 m，下探测器接收到的计数极少，三窗法剥离及修正处理中产生了局部“假值”，而全谱法计算结果低于三窗法计算结果，综合分析更符合该区地质情况。