摘要
当前航磁测量朝着更高精度和空间分辨率方向发展，多旋翼无人机具备低空仿地飞行的优势，成为大比例尺地质调查工作中的重要手段。多旋翼无人机按照预装在飞控系统中的规划路径执行飞行任务，进而完成航磁测量工作，因此航线规划成为实现低飞仿地航磁数据采集的关键因素。常见的路径规划软件普遍以覆盖测量为主，具有系统封闭、坐标系统单一和对DEM兼容性差等不足，开发适用于航磁测量的航线规划软件是十分必要的。本文以GMT开源软件为平台，依据其模块功能，提出平面测网布置和空中路径规划的工作流程，编写bash脚本程序。利用工区范围坐标、测线间距、测线方位和测线外延长度等已知信息生成平面测网，以平面测网为基础在测线上插值坐标点，通过DEM数据得到规划空中路径点高程，给定公差值对空中航路点进行简化处理，最终生成适用于无人机遥控器的KML格式文件，为多旋翼无人机执行仿地航磁测量工作提供了具有现实意义的解决方案。
Abstract
The current aeromagnetic survey is developing towards higher accuracy and spatial resolution. Multi rotor UAV have the advantage of low-altitude terrain following flight, and have become an important tool in large-scale geological survey work. Multi rotor UAV perform flight tasks according to the planned path preinstalled in the flight control system, thereby completing the aeromagnetic survey work. Therefore, route planning becomes a key factor in achieving low-flying terrain following aeromagnetic data acquisition. Common path planning software generally focuses on coverage measurement, with limitations such as closed systems, single coordinate system, and poor compatibility with DEM. It is necessary to develop a route planning software suitable for aeromagnetic survey. This article uses the GMT open source software as a platform, based on its module functions, proposes a workflow for plane survey network layout and air path planning, and writes bash script programs. Using known information such as work area coordinates, line spacing, line orientation, and line extension length to generate a plane survey network, interpolate coordinate points on the survey line based on the plane survey network, obtain the elevation of planned air path points through DEM data, simplify the air route points with given tolerance values, and finally generate KML format files suitable for UAV remote control. This provides a practical solution for Multi rotor UAV to perform terrain following aeromagnetic survey work.
Keywords
0 引言
航空磁测是一种快速、高效的航空物探方法（王庆乙等，2016），近年来，以多旋翼无人机为平台的航磁测量系统在地质调查和矿产勘查中得到了广泛的应用。多旋翼无人机可以在距地面更近的飞行高度工作，适用于复杂地形条件下的大比例尺航磁数据采集（Douglas and Peacker，1973），进而对航磁异常分辨率提出了更高的要求（陈斌等，2010；李月臣等， 2020）。飞行高度是分辨弱磁异常的重要因素（梁盛军等，2017），低飞测量可以获得更高的空间分辨率航磁数据（Parshin et al.，2018；Walter et al.，2020），保持恒定的离地高度有利于减少地形引起的磁响应 （Fairhead，2015），因此仿地飞行是低飞航磁测量发展的必经之路，也是当前航磁工作的难点问题。
路径规划是无人机应用研究中的重要内容 （Gugan and Haque，2023）。无人机路径规划与应用场景相关，可分为静态路径规划和动态实时路径规划两类。静态路径规划通过对已知环境信息分析得到优化的飞行路径，将其装载在无人机上使其沿预定航线飞行（高晖等，2001）；动态路径规划以最优路径选择、避障、目标检测与识别为研究内容 （Vashisth et al.，2021）。航磁的仿地飞行属于静态路径规范范畴，无人机航磁系统沿一定高度地形趋势面进行测量（王猛等，2022）。通过高分辨率数字高程模型（DEM）和测线端点可以取得测线上的路径点的高程数据（武雷超等，2023），把路径点高程提升到设计飞行高度，并对路径点简化即可得到规划航线，将规划航线文件传送到多旋翼无人机遥控器系统，机载航磁系统在规划路径上执行低飞仿地飞行数据采集，从而达成高分辨率航磁测量的目的 （张文杰等，2021）。
当前存在的无人机空中路径规划软件主要为无人机厂商及商业公司开发，通常侧重于覆盖工区的航线设计，以经纬度坐标规划航线，不利于基于投影坐标的航磁工作需求。部分软件可以根据 DEM数据实现仿地飞行功能，但对数据源和数据量均存在限制，随着航磁测量朝着低空—超低空方向发展，对DEM数据的精度和分辨率提出了更高的要求，不利于执行低飞测量。与 PIX4D、DroneDeploy 等商业软件和无人机厂商开发的航线规划软件相比，基于开源平台编写脚本具有更大的开放性与灵活性。本文以GMT软件为平台，提出了航线规划脚本的开发流程，利用工区范围、测线距、设计飞行高度等已知要素，结合 DEM 数据编写脚本，生成适用于无人机遥控器的KML格式文件，适用于各型多旋翼无人机执行低空仿地航磁测量工作任务。
1 系统设计思路
航磁测量通常为面积性工作，工区范围为多个拐点构成的多边形呈现，航磁工作开展前，需根据测线方位、测线间距和测线外延长度等参数设计覆盖工区的航磁测网（Coombes et al.，2018），航磁测网由多条航线组成，无人机采用往复作业模式进行数据采集（Li et al.，2023），航磁测线不仅要考虑测线上的路径，也要规划测线间的路径以及起飞点与测线起点、降落点与测线终点间的路径。根据工区范围坐标，通过GMT软件求取工区中心点位置和最大测线长度，利用测线方位经过几何计算得到覆盖工区的测网，再通过使用工区范围多边形文件截断、外延测线长度等处理，可得到适合工区范围多边形的航磁测网拐点平面坐标文件，利用DEM数据对以上得到的测网文件进行适当的加密、赋高程和提升高度处理后，对测线长度和飞行高度组成的多段线进行简化处理，转换为经纬度坐标并输出为KML格式文件，就完成了低飞仿地测网的生成工作。
2 系统设计与实现
2.1 开发语言
GMT（Generic Mapping Tools）是一款开源地图绘制软件，具有强大的绘图功能和数据处理功能 （Wessel et al.，2019），本文以GMT6.4为开发平台编写 bash 脚本，可在 Windows、Linux 和 macOS 等系统下运行。
2.2 系统功能
系统主要实现生成平面测网坐标和空中规划航线两部分内容（图1）。平面测网坐标包括各测线端点坐标，用于绘制实际材料图及后期针对航磁原始数据的预处理工作。空中规划航线为简化后的航路点构成的KML格式文件，将其导入无人机遥控器即可以实现仿地飞行的工作目的。
图1工作流程图
2.3 程序的实现
航磁测量的航线规划涉及到平面测网布置和空中路径规划两部分内容。
2.3.1 平面测网布置
平面测网的生成包含 4 个步骤：①由工区范围多边形文件计算中心点坐标和最大测线长度，以最大测线长度计算主测线点坐标；②计算通过主测线测点的测网端点坐标；③利用工区范围多边形文件截断测网文件得到各线性在多边形上的坐标；④测线按设定的长度进行外延，取得相应的坐标。以下对计算过程进行简述：
（1）中心点坐标计算
商业软件一般基于经纬度坐标组成的工区多边形数据进行测网布设，在测线方位为 0°、90°等条件下，测线难以通过投影坐标公里网格布设，因此需以工区多边形投影坐标为坐标计算的基础，对求取的中心点进行适当的偏移。如以南北向（0°）布置测线时，通过 spatial 模块计算得到投影多边形数据 （polygon.txt）中心点坐标后，传递给info模块，调整X 坐标（
x_mc）为测线间距（line_dis）的倍数，这样 X 坐标就可以落在以线距为倍数的公里网格上，对于 Y 坐标（y_mc）则没有限制，此时取得的不是中心点，可称为似中心点（图2a），南北向测线中心点坐标求取代码如下：
x_mc）为测线间距（line_dis）的倍数，这样 X 坐标就可以落在以线距为倍数的公里网格上，对于 Y 坐标（y_mc）则没有限制，此时取得的不是中心点，可称为似中心点（图2a），南北向测线中心点坐标求取代码如下：
x_mc=`gmt spatial polygon.txt-Q+p | gmt info-I line_dis-C-o0`
line_dis-C-o0`
y_mc= `gmt spatial polygon. txt-Q+p | gmt info-I--C-o2`
反之，对于正东西向测线，需对Y坐标进行倍数调整，X坐标不做限制。对于测线方位不以正南北、正东西出现的测线，中心点坐标都无需取测线间距倍数。实际使用中，可以增加 if语句对测线方位进行判断，进而让脚本自行选择相应的中心点求取公式。
（2）计算覆盖半径
为了使生成的测网能完整覆盖工区范围，需用 mapproject 模块中心点到工区多边形文件的距离，再使用 math 模块统计出最大距离 lmax（图2a），并调整为满足测线间距的倍数的最短长度（ex_mc），代码如下：
ex_mc），代码如下：
ex_mc= `gmt mapproject polygon. txt-G x_mc/ y_mc+uc | gmt math STDIN-C2-S-o2 UPPER line_dis DIV CEIL line_dis MUL =`
x_mc/ y_mc+uc | gmt math STDIN-C2-S-o2 UPPER line_dis DIV CEIL line_dis MUL =`
（3）生成主测线点坐标
使用 project 模块以中心点（x_mc/y_mc）为中心，按测线间距（line_dis）沿垂直测线方向（azi‐ muth）按最大长度（-ex_mc/ex_mc）生成主测线上点坐标（图2b），代码如下：
x_mc/y_mc）为中心，按测线间距（line_dis）沿垂直测线方向（azi‐ muth）按最大长度（-ex_mc/ex_mc）生成主测线上点坐标（图2b），代码如下：
gmt project-C x_mc/y_mc-A azimuth-G line_dis-N-L-ex_mc/ex_mc-o0, 1 >point.gmt
x_mc/y_mc-A azimuth-G line_dis-N-L-ex_mc/ex_mc-o0, 1 >point.gmt
（4）生成覆盖多边形工区的测网
GMT通常采用 project模块生成测网，也可以使用 grdtrack 模块生成，该模块更适合生成航磁测量所需的往复线，根据主测线点坐标，在主测线全长 ex_mc_2（ex_mc的 2倍）上以测线间距（line_dis） 仅生成主测线半长（ex_mc）的测线坐标，该坐标包含测线两个端点和中心点（各主测线点），得到覆盖工区的正方形测网，代码如下：
ex_mc_2（ex_mc的 2倍）上以测线间距（line_dis） 仅生成主测线半长（ex_mc）的测线坐标，该坐标包含测线两个端点和中心点（各主测线点），得到覆盖工区的正方形测网，代码如下：
gmt grdtrack point. gmt-G dem-C ex_mc_2/ex_mc/line_dis+a >line.gmt
dem-C ex_mc_2/ex_mc/line_dis+a >line.gmt
（5）以工区多边形截断测网并延长测线
使用spatial模块以工区范围多边形文件截断测网文件（图2c），再将每条测线两端按设计的长度 （line_ex）延长，可利用simplify模块排除中心点后，得到工区测线端点文件（图2d），代码如下：
line_ex）延长，可利用simplify模块排除中心点后，得到工区测线端点文件（图2d），代码如下：
gmt spatial line. gmt-Tpolygon. txt | gmt spatial-Wline_ex | gmt simplify-T1 >line.txt
line_ex | gmt simplify-T1 >line.txt
图2根据工区范围生成测线过程示意图
a—生成主测线点；b—生成扩边测线；c—截断测线；d—测线等距外扩
工区范围主要分为凸多边形和凹多边形，测线距 500 m，外延 100 m，联络线距 2 km，外延 1 km 在凸多边形生成测线效果见图3，在凹多边形生成测线效果见图4。
2.3.2 空中航线规划
规划航线是在平面测网基础上，通过为测线加密坐标，利用 DEM 文件取得各点地表高程，由设计飞行高度可得航路点海拔高度，简化航路点平面坐标投影长度和航路点高程组成的多段线，转换为经纬度坐标后输出为 KML 格式文件。具体流程详述如下：
图3凸多边形工区测线生成结果
a—测线方位：0°；b—测线方位：40°
①与平面测网中的测线不同，规划航线还包括了测线间的线段，因此必须将多段往复测线中增加测线间线段的坐标，可通过 awk 命令整体复制一次端点坐标，再插入多段线段头标识符的方式生成； ②根据生成的测线文件，通过 sample1d模块按一定的间距插值坐标，插值间距主要与DEM数据分辨率有关，一般可选择 DEM 网度作为插值间距；③用 grdtrack 模块根据 DEM 数据取得各点的地表高程值；④传递给 split 模块增加累积距离列；⑤用 sim‐ plify 模块以累积距离和飞行高度（高程与设计飞行高度height之和）为横、纵坐标，以给定公差（toler‐ ance）对多段线做简化，得到简化后的规划航线；⑥ 用 mapproject模块将简化后的规划航线转换为经纬度坐标；⑦使用 ogr2ogr模块转换为 KML 格式文件，该文件即为工区规划航线文件，可供无人机遥控器使用，对于不同的机型还需根据具体需求修改KML 文件内容。
height之和）为横、纵坐标，以给定公差（toler‐ ance）对多段线做简化，得到简化后的规划航线；⑥ 用 mapproject模块将简化后的规划航线转换为经纬度坐标；⑦使用 ogr2ogr模块转换为 KML 格式文件，该文件即为工区规划航线文件，可供无人机遥控器使用，对于不同的机型还需根据具体需求修改KML 文件内容。
图4凹多边形工区测线生成结果
a—测线方位：0°；b—测线方位：40°
在实际使用中，仅需使用部分航线，并将航线与起飞点和降落点组合成航线，对于起降点与测线端点间为简单开阔地形时，可以直接截取使用所需的航线规划文件；对于复杂地形，可将起降点与所需测线端点建立新的测线，使起降点到测线端点间同样执行仿地飞行，当然这部分不处于测线范围内的路径也可以单独生成飞行高度较高的航路点并拼接到规划航线上，但需注意确保各段航线正确衔接。
2.4 关于地形简化的算法
地形简化的关键在于使用少量特征点，根据简化程度表达多尺度地形（张娜等，2022）。
GMT的 simplify 模块基于 Douglas-Peucker 算法对复杂线段进行简化，Douglas-Peucker 算法使用点-边距离作为衡量误差的标准（Douglas and Peacker，1973），是线状要素抽稀的经典算法（张慧莹等，2019），前人将该方法应用于无人机飞行航线规划（梁文勇等，2020；王炳乾等，2020）。图5给出了该算法的计算过程，首先连接原始线段的首尾点形成直线，计算所有点到直线的距离，选出距离该直线最远的点，如其到直线的距离大于给定的公差，则添加到简化的结果中，形成新的多段线，重复以上的距离计算和公差比较的过程，并不断更新简化结果（图5b、c），最终得到简化后的线段（图5d）。
图5Douglas-Peucker算法对线段简化过程
a—首尾连接取最大偏差点；b—分段后分别取最大偏差点；c—进一步分段；d—完成最终的简化线段
2.5 脚本运行中的问题和解决办法
本文以 GMT 为平台开发适用于航磁仿地低飞航线规划脚本，其工作流程是可行的，但 GMT 下的模块的部分功能还存在不足和瑕疵，如 spatial 模块的一些选项在当前（6.4）或更早的版本存在错误。因此尽管脚本设计是可行的，为了使程序顺利运行，还需对代码做一些合理的替代性修改。
2.5.1 关于多边形裁剪
最早出现于 GMT5中对点、线和多边形进行地理空间处理的 spatial 模块，具备强大的分析与计算能力，但其在处理多边形截断测网文件时容易出现问题，当测网中的测线恰好通过多边形文件的某一条边时，该条测线与多边形的交点计算会出错，得不到正确的截断后的测线端点坐标。
为了提取多边形和正方形测网的完整交点坐标，可以对多边形范围文件进行微扩边，多边形文件数值小数位为 n，将其扩大至 10-（n+1），如多边形坐标为整数时扩大 0.1 m，1 位小数时扩大 0.01 m。以微扩边后的多边形截断测网，再通过控制输出结果精度，即可得到正确的交点坐标，代码如下：
gmt spatial line. gmt-Tpolygon_exp. txt--FORMAT_FLOAT_OUT=%.0f
其中，line.gmt 为正方形测网文件，ploygon_exp. txt 为微扩边后的文件，临时设置参数 FORMAT_FLOAT_OUT 用于控制输出数据小数位，需根据原始多边形坐标小数位修改。
2.5.2 关于缓冲区计算
建立在微扩边多边形文件的基础上的对测网文件截断并取得坐标是可行的，微扩边文件同样使用 spatial 模块计算，然而遗憾的是该模块用于计算缓冲区的-Sb 选项同样存在瑕疵，计算复杂多边形的缓冲区会出现飞点（图6a），虽然可以通过手工去除，但不利于脚本连贯运行。GMT 引入了 GDAL/ OGR 开源库，以 ogr2ogr 模块的形式使用，如二者的 implify 模块同样基于 Douglas-Peucker 算法，处理结果相当，可以利用 ogr2ogr 模块中的 st_buffer 函数完成缓冲区的计算，代码如下：
ogr2ogr-f "GMT" polygon_exp. gmt polygon. gmt-dialect sqlite-sql "select st_buffer (geometry, exp) from polygon"
exp) from polygon"
其中，polygon.gmt为原始工区范围多边形文件， polygon_exp.gmt 为微扩边后的多边形文件，ogr2ogr 利用 sql 数据库语言中的 st_buffer 函数进行缓冲区的计算，geometry为几何计算，exp为缓冲区宽度（微扩边长度）。ogr2ogr得到的缓冲区文件包括外扩多边形和内缩多边形（图6b），再利用 select 模块选择在原始多边形外的闭合多边形就得到了微扩边多边形文件。
2.6 脚本特点及优势
以 GMT 软件为平台开发的仿地飞行航线规划脚本，助推实现高分辨率航磁数据采集，其具有以下几点优势：
图6spatial和ogr2ogr多边形等距外扩对比图
a—spatial扩边结果；b—ogr2ogr扩边结果
（1）脚本执行效率高
基于GMT平台的脚本存在交互上的不利因素，但其利用的模块相对较少，使用难度不大，输入已知变量即可在数秒内生成任意形状工区的航线文件。
（2）以投影坐标布设航线
常见航线规划软件普遍基于经纬度坐标生成飞行路径，对航磁采集常见的正南北（东西）向测线支持不足，本脚本通过条件语句对预设测线方位进行判断，对于正南北（东西）向测线，利用预设线距使航线按投影坐标公里网格布置，有利于规范制图和成果数据阅读。
（3）支持多格式、数据源DEM/DSM模型
目前已有无人机可通过遥控器下载 DEM 数据实现仿地飞行，但一般都指定了 DEM 数据源，且有数据量限制。随着飞行高度的降低，低精度的DEM 已不适应低空仿地飞行的需要，飞行高度的降低对 DEM数据的精度和分辨率提出了更高的要求，对于树木、建筑较多的工区，还需要使用高现势性的 DSM 数据，甚至还需预先航测才能确保飞行安全。 GMT可以直接读取网格格式的DEM数据，对于更为广泛的数据格式，可通过GDAL程序进行转换，从而满足GMT脚本对坐标点高程的提取需求，因此GMT 脚本程序在DEM数据的选择上更具开放性。
（4）航路点没有数量限制
随着航磁勘探比例尺的加大，势必通过降低飞行高度来取得高分辨率资料，在复杂地形条件下，为满足仿地低飞的要求需加大航路点密度，对于工区航线来说航路点总数是极大的。脚本对生成的航路点没有总数限制，无需分块生成航线。
3 应用实例
本文以“大兴安岭及胶辽重要找矿区矿产调查”课题为例，生成航磁工作所需航线规划文件。为便于显示，将线距抽稀为 1 km，以 Douglas-Peucker算法按不同公差计算得到不同的航路点（图7b～d），对比工区 DEM 数据（图7a）可见，工区南北两侧为山区，航路点密集，中部地势平缓区域航路点稀疏。同时，随着公差的增大，航路点被逐步抽稀。
规划飞行高度100 m的航线在三维地形图和影像图上的显示见图8，山区以 10 m 公差简化后的航点仍然较多，说明地形高差变化较大，因此更换测线时加密航路点是十分必要的。此外，对于复杂多边形，仅用 1 条航线可能无法完成形成正规测网覆盖工区（图7东北部测线为独立航线），尽管目前有很多覆盖规划路径的算法，但其基于拐点坐标完成的方式势必截断航线，对于航磁工作是不利的，通过选择线号范围分别存储航线文件的方式仍然是最有效的解决办法。
图7不同公差的航路点分布图
a—某工区地形图；b—公差5 m航路点分布；c—公差10 m航路点分布；d—公差20 m航路点分布
图8某工区航线设计
a—规划航线在三维地形上的显示；b—规划航线在三维卫星图上的显示
4 结论
（1）本文以 GMT 软件为开发平台，提出了适用于无人机航磁系统的航线规划系统开发流程，以脚本程序生成KML格式文件，可为各类无人机遥控器提供空中航线引导，为实现高精度航磁资料提供了可靠的路径规划支持。
（2）脚本程序应用了成熟的线要素简化算法，在保证无人机实现仿地飞行的前提下，对飞行路径进行了优化，有助于提高无人机飞控系统的稳定运行。
（3）程序对DEM数据源和数据格式具有良好的兼容性，能够处理复杂地区包含海量航路点的航线生成任务，提高了航线规划的灵活性和适用性。
（4）本文指出了GMT软件在地理空间处理方面的一些不足，并提出了简单有效的解决方案和代码，为后续研究和开发提供了重要参考。
综上所述，本文开发的基于GMT软件的无人机航磁航线规划脚本程序在航线规划精度、兼容性和解决方案创新等方面具有显著优势。该系统为提高无人机航磁工作测量精度和空间分辨率提供了重要支持和技术手段，对于推动无人机航磁领域的进一步发展具有重要意义。