摘要
多旋翼无人机航磁是近年来发展起来的一种低飞航磁系统，具有仿地飞行能力，可在特殊景观区等人力不可企及的区域开展工作，但其受限于续航时间，通常以设置多个起飞点的方式分区块施工，在实际测量中往往需频繁更换起飞点。低飞航磁易受风的影响，不可能严格按照规划航线飞行，在飞行测量过程中遇强电磁干扰可能出现不正常的磁场值，甚至导致磁力仪信号失锁而采集不到数据，因此有必要在完成每架次航磁测量工作后及时评价数据质量，对不合格的测线进行问题分析并及时补测。以往航磁质量评价过程较繁琐，需截取各测线数据后对飞行质量和航磁质量分别进行评价。本文通过在GMT软件平台下编写脚本，开发出适用于多旋翼无人机航磁的质量评价系统，实现了测线数据提取和飞行质量、航磁质量的统计计算，结合卫星影像和DEM数据输出为矢量化图件，便于野外人员快速分析和判断航磁工作质量，对提升航磁野外工作效率、保障施工连续性方面起到了促进作用。
Abstract
The Aeromagnetic system of multirotor is a recently developed low-flying aeromagnetic system, which possesses ground-like flight capabilities and can operate in inaccessible terrains. However, due to its limited endurance, it often requires multiple takeoff points for measurement purposes, resulting in frequent replacements of these points. Additionally, the low-flying nature of this system makes it susceptible to wind disturbances and deviations from planned routes. Strong electromagnetic interference may lead to abnormal magnetic field values during flight measurements, potentially causing loss of magnetometer signal and data collection failure. Therefore, timely evaluation of data quality is essential after each aeromagnetic measurement flight to identify any unqualified survey lines and conduct supplementary measurements accordingly. Previously, the process of aeromagnetic quality evaluation was cumbersome as both flight quality and aeromagnetic quality had to be assessed separately after intercepting data from each survey line. In this study, we have developed a script-based quality evaluation system for multirotor aeromagnetic measurements using the GMT software platform. This system enables efficient extraction of survey line data along with statistical calculations for assessing both flight quality and aeromagnetic quality.
Keywords
0 引言
航磁测量是发展最早、应用最广泛的航空物探方法，广泛应用于金属矿产勘查、油气资源勘探、基础地质研究、地球深部探测等领域（熊盛青，2020）。近年来，随着国内外无人机技术的快速发展和磁力仪小型化的实现，无人机航磁技术得到了不断的开发和完善。无人机平台以固定翼、直升机和多旋翼为主（崔志强等，2014）。中、大型固定翼无人机（李文杰等，2014；陈晓云等，2019；闫东和周乃恩， 2020；郭华等，2022；周普志等，2022）和直升机无人机（黄岩等，2019；西永在等，2019）具有航时长、载荷大的特点。多旋翼无人机可以在距地面更近的飞行高度工作，采样密度和数据分辨率更高（Parvar et al.，2018；Walter et al.，2020），具备仿地飞行能力，可在特殊景观区、中高山地区等人力不可企及区域开展工作（张文杰等，2021）。但其受载荷、电池容量等因素制约（Mu et al.，2020），续航时间较短（刘双等，2023），通常将工区分割成若干区块，利用多个起飞点完成整个工区的航磁测量工作，野外工作中需频繁转场，航磁测量过程中遇强电磁干扰可能出现不正常的磁场值，甚至导致磁力仪信号失锁而采集不到数据，因此现场评价航磁资料质量是十分必要的。
航磁资料的质量评价是基于测线数据进行的，然而工作过程中通常在地面将磁力仪调试正常后再进行飞行测量，故航磁原始数据既包含航线规划的测线数据，也包含测线外的数据，因此提取测线数据是进行质量评价的前提。
航磁资料的质量评价分为飞行质量和航磁数据质量两方面内容。飞行质量包括导航（偏航距） 和飞行高度两部分，多旋翼无人机按规划的航线自主飞行（崔志强等，2016；严月浩等，2022），飞行姿态易受风的干扰（Wang et al.，2019），航线的高程也并非严格按照设定的高度（常宜峰等，2016），而飞行高度对分辨弱磁异常的影响很大（梁盛军等， 2017），因此飞行质量关系到航磁资料整体质量。航磁数据质量通过每条测线的动态噪声进行评价。目前航磁工作质量评价多基于人工或多软件平台结合的后处理，不利于野外工作人员快速分析数据质量。由于多旋翼无人机航磁测量机动施工的特点，开发可用于野外的质量评价系统具有现实意义。
本文以GMT软件为平台编写脚本，航磁原始数据经预处理后提取测线数据，统计每条测线偏航距和磁测数据的动态噪声水平，以DEM数据为参考计算飞行高度。加入卫星影像图和地形图，以平面图、剖平图和统计图表的形式呈现测线数据质量情况，野外人员据此可快速评判航磁工作质量。
1 系统设计思路
航磁工作质量评价分为飞行质量和航磁数据质量 2 部分，质量评价工作建立在提取规划航线范围内的有效数据基础上，因此要对航磁原始数据进行整理并拆分保存为测线数据，逐个评价每条测线数据质量。飞行质量评价包括导航精度和飞行高度质量 2 方面内容。首先利用 GMT 软件的数据处理模块以规划航线提取规划航线范围数据并拆分为单独测线，计算各测点偏航距，统计每条测线偏航距范围和平均值，达到评价各测线导航质量的目的；再利用 GMT 软件由 DEM 网格文件取得各测点地表高程，进而计算出飞行高度，统计每条测线飞行高度范围和平均值，对各线飞行高度质量作出评价。航磁数据质量通过动态噪声水平评价，按四阶差分方法计算每条测线的动态噪声水平。
以上对航磁原始数据的统计计算，得到了各线的航磁动态噪声水平、导航精度和飞行高度质量统计结果。但统计值对数据的反映不够直观，野外人员仅能了解到数据是否正常，而不清楚问题出在何处。因此需要把偏航距、飞行高度和磁总场数据绘制于平面图上，借助卫星图片、地形等高线和已知干扰因素等外部参考数据，利用GMT软件的绘图模块生成图表样式的质量评价统计结果，对飞行质量和航磁数据质量进行直观评判。及时发现不合格的数据所在的测线和平面位置，分析不合格磁异常资料是否由踏勘明确的干扰源引起。并将每条测线的统计结果以统计表和直方图的形式依次展示在平面图之下，从而达到了对每架次、每条测线航磁数据进行综合表述。
在利用 GMT 脚本编程实现对航磁原始数据的预处理和质量统计的同时，也统计出各测线的起止坐标、起止时间、采集时长和无人机平均飞行速度等信息，有助于野外人员快速评价航磁原始数据质量，及时发现问题，分析影响因素，对不合格的测线及时返工，避免重复的进出起飞点，完善多旋翼无人机航磁施工流程，提升整体工作效率和工作质量。
2 系统设计与实现
2.1 开发环境和语言
GMT（Generic Mapping Tools），中文名“通用制图工具”，是地球科学使用最广泛的制图软件之一，具有强大的绘图功能和数据处理功能（Wessel et.al.，2019），本文以 GMT6.4 为开发平台，在 Windows 下通过 Git for Windows 运行 bash 脚本，也可以在 Linux或macOS下执行bash脚本。
GMT 软件提供了 UnixTools 脚本工具，包括 awk、paste、ls、grep等Linux下的常用命令。awk是适于提取数据和标准化输出的文本处理工具（Aho et al.，1979），在脚本编程中应用最多，本文使用 GNU awk 5.2 版本进行数据提取和输出（Arnold，2022）； paste 命令用于按列合并文件；ls和 grep 用于从目录文件夹中提取指定的格式的数据。
2.2 系统功能
系统功能主要分为预处理和图形显示 2 部分 （图1）。
预处理：主要完成航磁数据的测线数据提取、飞行质量和航磁数据质量统计等功能。
图形显示：根据航磁资料质量评价统计结果，结合卫星图片、DEM 数据、已知干扰因素等绘制图表，输出保存为矢量化图件。
图1无人机航磁野外资料质量评价系统功能
表1数据文件格式说明
图2无人机航磁野外资料质量评价系统工作流程图
XTK—偏航距；AGL—飞行高度；TMI—磁总场值
2.3 程序的实现
利用GMT软件对航磁原始数据进行预处理，提取测线数据，并对其进行飞行质量和航磁数据质量的统计和评价，以图表的样式输出。无人机航磁野外资料质量评价系统工作流程见图2，数据处理过程中的文件格式见表1。
2.3.1 预处理
预处理包括测线数据提取、飞行质量统计和航磁质量统计，测线数据提取后生成 txt 文件，各项质量统计结果赋值给不同的变量，GMT绘图时调用这些变量即可将参数和数值绘制于各类图表上。
（1）测线数据提取
由于航磁系统组成不同，可能包括磁补偿器和飞行高度测量设备采集的数据，这部分数据对于质量评价没有意义，因此需要对航磁数据进行简化，提取与工作质量有关的数据列，以便于拆分测线数据和质量评价。但需注意原始数据中测线的单一性，以免错误读取数据。
①原始数据 awk 命令包括经纬度坐标、海拔高程、磁场值和时间戳，保存为 ini.tmp 文件（表1），该文件包含所有航迹，可追溯无人机起降和工作过程。
②根据规划航线文件（Plan.txt），使用 select 模块-L 选项筛选测线范围数据，通过管道传递，使用split 模块以给定方位角（-A 选项）或方位角变化量 （-C 选项）进一步筛选并拆分为测线数据（mag. tmp），代码如下：
gmt select ini. tmp-fg-i0: 3-Lplan. txt+d10e+p | gmt split-fg-C10-D2k >mag.tmp
由于select模块只能输出5列数据，时间戳需另外输出为 time.tmp文件，再通过 paste命令合并 mag. tmp、time.tmp 文件，经过管道传递，使用 awk 命令将所需数据列保存为sor.tmp文件。
③利用mapproject模块计算并增加偏航距列，通过管道传递到 grdtrack 模块，利用 DEM 网格文件提取坐标点对应的地表高程数据，最后用 awk 命令计算飞行高度、编辑线号并生成 profile.txt文件。该文件汇总了所有的测线数据，并包含有偏航距和飞行高度数据，用于各质量参数剖平图绘制，代码如下：
# 设定数据范围
R=`gmt info sor.tmp-I0.01/0.005+e`
gmt mapproject sor. tmp-fg-Lplan. txt+ue | gmt grdtrack-GDEM.grd 
>profile.txt
>profile.txt
利用 convert 模块将 profile.txt 文件拆分为各条线的nc格式数据（profile_0.nc、profile_1.nc等），用于每条测线数据统计。以0起始的拆分后测线可明确表示航磁测量的先后顺序，所以没有按照规划航线线号编排，代码如下：
gmt convert profile.txt-Dprofile_%01d.nc-fo0x, 1y
（2）飞行质量统计
飞行质量评价分为导航质量（偏航距）和飞行高度两部分，利用 ls和 grep 命令读取拆分后的测线数据（nc 格式），生成列表，遍历该列表内的所有测线文件，并依次对各条测线进行统计和评价。nc格式测线数据中第1行为线号，第2行起由8列数据组成（表1）。
数据统计过程中需大量定义和使用变量，主要通过 awk 命令和 math 模块进行变量的定义，以 awk 命令取测线开始纬度和结束纬度为例，代码如下：
#起点纬度：某线文件第 2 列第 2 行（第 1 行为线号）
lat_start= 
#终点纬度：某线文件第2列最后1行
lat_end= 
通过使用 math 模块取各线文件的偏航距 （XTK）和飞行高度（AGL）的极值和平均值。如偏航距的最小值、最大值和平均值可由以下代码完成：
#最小偏航距：某线文件第7列极小值
xtk_min=`gmt math $nc-i6-S LOWER =`
#最大偏航距：某线文件第7列极大值
xtk_max=`gmt math $nc-i6-S UPPER =`
#平均偏航距：某线文件第7列算术平均值
xtk_mean=`gmt math $nc-i6-S MEAN =`
（3）航磁质量统计
通过航磁总场值变化范围可以快速判断测线是否出现明显不正常的数据，同样利用math模块进行极值的求取。
航磁数据动态噪声水平按四阶差分方法评价，主要使用GMT的math模块和awk、paste命令。由于航磁所用磁力仪采样率通常远高于 2 Hz，故先利用 awk对数据进行抽稀为 0.5 s间隔处理，利用 math模块下的DIFF求取相邻点的距离差值和磁总场差值，再用awk命令按Ti-2、Ti-1、Ti、Ti+1、Ti+2分别取出磁总场值，使用paste合并以上5列数据，得到同一行下5个时间的磁总场数据，利用 awk 命令进行四阶差分计算，剔除磁场变化大于600 nT/km数据后，使用math 模块求取总测点数和平均磁场值，最终统计出测线动态噪声水平。
2.3.2 图形显示
通过对航磁原始的预处理，得到了测线范围内的航磁数据和每条测线上飞行质量、航磁质量的统计结果，再通过GMT软件将实测数据和质量统计结果输出为矢量化图件。
图形显示的成果表达，采用“3平面图+N统计图表”的形式布局。“3 平面图”位于图像顶端，由飞行航迹图、飞行高度剖平图和航磁总场剖平图组成，用于直观反映飞行质量和航磁数据质量。“N 统计图表”位于“3平面图”之下，为 N条测线的质量评价统计表和直方图，用于明确各条测线的工作质量。
“3 平面图”为范围相同的平面图，利用 subplot 模块设置为3张子图并依次绘制。由于航磁原始数据中的测线条数是不确定的，而GMT软件绘图原点位于左下角，不改变绘图原点会造成“3 平面图”下的“N 统计图表”显示不全。因此需在 subplot begin 中设置-Y 选项，沿 Y 方向增加一个比较大的高度 （如-Y100c即绘图原点向上提升100 cm），可防止测线质量评价图表被截断的情况发生。
由于测线的条数是可变的，因此单条测线的质量评价图表需用循环语句和subplot模块完成绘制。首先 ls 和 grep 命令新建包含所有测线（nc 格式）的列表，使用 whlie 循环语句遍历该列表文件，再用 subplot 模块逐个绘制每条测线的质量评价统计结果。
（1）飞行航迹图
飞行航迹图为以卫星影像图为底图的平面图，叠加规划航线和实测航迹。利用grdcut模块截取平面图范围内卫星影像图片，并用 grdimage 模块显示该图像，再用 plot 模块依次添加规划航线（Plan.txt） 和实测航迹（ini.tmp），最后利用 convert 和 text 模块将profile.txt文件中的线号标于图上。
（2）飞行高度剖平图
飞行高度剖平图以地形等值线平面图为底图，叠加飞行高度剖面异常。利用 grdimage 和 grdcon‐ tour模块将 DEM 网格转化为地形等值线平面图，以设计飞行高度为参考，利用 wiggle 模块将 profile.txt 文件中的飞行高度列叠加于图上，实测航迹和线号亦可标于图上。
DEM 的分辨率会直接影响到质量评价的精度 （翟延亮等，2022）。目前开源 DEM 数据众多，Li et al.（2022）评价了NASADEM、AW3D DEM和Coperni‐cus DEM 等 DEM 数据的质量，以上三者的垂向精度均在 6.59～6.73 m，区别不大。DZ/T 0405-2022《无人机航空磁测数据采集技术要求》中对DEM垂向精度的要求是优于设计飞行高度的 20%，因此一般条件下以上几种 DEM 数据均可用于评价飞行高度。但需注意的是绝大部分地表高程数据均为 DSM 数据，包含了地貌因素，而航磁工作评价的飞行高度实际指的是离地高度，因此选择基于 Copernicus DEM 的去地貌数据——FABDEM 可能是更好的选择，研究发现，与 Copernicus DEM 数据集相比， FABDEM 的高程 RMSE 降低了 24%，偏差降低了 135%（Christopher et al.，2023）。因此应尽量选用 DEM数据进行飞行高度的计算。
（3）航磁总场剖平图
航磁总场剖平图主要分析航磁总场值是否出现突变。以 math 模块统计的所有测线磁总场平均值为参考，利用 wiggle 模块将 profile.txt 文件中的航磁总场值列以剖平图形式绘制，再添加航迹、标注线号。此外，可通过将踏勘整理的已知干扰因素保存为 Info.txt 文件，使用 plot、text 模块将干扰位置以点、线和注释的形式标于图上，便于分析不正常的数据是否由已知干扰引起。
（4）测线质量评价图表
由测线质量统计表和直方图（偏航距、飞行高度、磁总场值等）组成。质量统计表包括当前测线的起止坐标、测线长度、起止时间、采集时长、飞行平均速度、动态噪声等统计结果。直方图反映各参数的分布情况，横轴为参数变化范围，纵轴为百分比，内注参数值范围和平均值。利用subplot模块设置子图，再通过 text 模块将统计结果得到的变量标注为文字，使用histogram模块绘制直方图，并用leg‐ end模块在直方图中插入统计值。
2.4 系统特点及优势
通过在GMT平台下编写脚本，开发出多旋翼无人机航磁质量评价系统，对快速判断航磁工作质量具有现实意义，该质量评价系统具有以下优势：
（1）脚本运行速度可满足快速评价要求
通常脚本在1～2 min完成运行，生成图像，达到了快速评价的目的。除计算机本身性能外，影响脚本运行速度的主要因素为卫星图片的切割和图片的导入，可通过适当降低图片的分辨率或提前准备当前工作日内数据的方式，减少数据的读取和写入时间。
（2）卫星影像和DEM数据有助于直观判断
卫星影像和地形等高线的加入，对分析飞行质量有指导意义，通过不同地形、地貌上的偏航距和飞行高度信息，可以了解航磁系统的工作状态，判断因地形突变导致的风力作用对航磁系统的影响，同时可对航线规划中出现的不合理现象分析，进而在后续工作中加以纠正。
（3）航磁质量评价图表有利于判断工作质量
每条测线的飞行质量和航磁质量统计结果以图表的形式输出，便于使用者判断出现问题的具体测线，是平面图的有力补充。
（4）航磁测线数据可直接用于后处理
拆分后的原始数据是后续数据处理的基础，标记线号的测线数据包含了坐标、时间、磁总场等数据，可直接用于后处理工作。
3 应用实例
“大兴安岭及胶辽重要找矿区矿产调查”课题所属中国地质调查局发展研究中心2020年度“重要锡、锰等矿集区矿产地质调查”二级项目，课题组在山东莱芜市开展了 1∶50000 低飞航磁工作，本文以其中 2架次航磁野外数据为例，利用编写的 GMT脚本对航磁原始数据进行预处理和质量评价。
3.1 总体质量情况
该份航磁原始数据由 4 条测线（2 次往返飞行） 组成，起飞点位于测线南部，图3中展示了飞行航迹图、飞行高度剖平图、航磁总场剖平图及4条测线的质量统计图表。
从飞行航迹图看，实际飞行航迹与规划航线偏差很小，没有出现大的偏航。无人机从起飞点起飞后沿直线飞行到 0 线南端工作，向北飞行完成设计航线任务后转向西侧 500 m 的 1线北端点向南完成测量后手控飞行至起飞点降落，更换电池后沿直线飞至 2 线南端向北飞行，至北端点后转向东 500 m 的 3 线进行测量，最后手控降落于起飞点。从飞行高度剖平图上，以总体设计飞行高度100 m为基准，测线南部地形较平缓，无人机飞行高度与设计飞行高度符合较好，北部山区飞行高度变化较大，总体表现为山脊飞行高度低，山谷飞行高度高。其中 3 线因地形复杂，飞行高度较总体设计飞行高度偏高（以紫色为主）。从航磁总场剖面图看，0 线北部出现明显的数据错误，其余各线总场值变化平稳。
图3航磁野外工作质量评价在山东莱芜项目中的应用
3.2 测线质量评价
从各线的飞行质量和航磁数据质量统计结果看，4条测线的飞行高度和偏航距均符合规程要求，其中西侧 2、3 线地形复杂区飞行高度略高（平均飞高 103.3 m、110.9 m），平均偏航距均小于 3 m。除 0 线外各线总场值变化范围不大，动态噪声最高 0.07 nT，满足一级质量标准。0 线总场值数据范围明显不正常，这一点在航磁总场剖平图上已有显示，其动态噪声达104 nT，据此判定0线需进行补充测量。
4 结论
本文以GMT软件为平台，开发了针对多旋翼无人机航磁测量的质量评价系统，实现了对飞行高度、偏航距和动态噪声的质量统计。通过加入卫星影像、DEM 数据等外部参考资料，丰富了图件成果的表达形式，可读性强，便于使用者做出正确的判断，及时对不合格的测线进行返工，从而避免重复进出场造成时间的浪费。同时，经过预处理完成了航磁原始数据到测线数据的转换，可直接用于航磁资料后处理。基于 GMT 脚本的多旋翼无人机航磁质量评价系统针对导出后的航磁数据进行计算和统计，可为航磁资料质量做出基础评价，为完善航磁施工流程，提高工作效率提供了保障。随着无线数据传输技术的发展，无人机飞行数据和航磁数据实时发送到地面处理终端已逐步实现，飞行质量和航磁数据质量的实时评价将成为未来发展的趋势。