基于MapGis平台可视化筛选无人值守航空物探数据

杨金政， 邱崇涛， 田 宇， 高国林， 祁 程

(1.核工业航测遥感中心，石家庄 050002;2.中核集团公司 铀资源地球物理勘查技术中心(重点实验室)，石家庄 050002)

0 前言

航空地球物理勘探简称航空物探，具有效率高、成本低、便于大面积工作、探测深度较大等优点，是基础性和公益性地质调查、战略性矿产勘查的重要手段，是地质勘查现代化的标志之一[1]。为了降低人员成本与风险，2011年，中国国土资源部试行了航空物探测量无人值守工作方式，2013年由中国核工业集团公司和中国地质研究院又成功进行了无人机航空测量试验。无人值守的工作方式与传统方式的一个显著差异是飞机从起飞到落地，仪器使用同一个测线编号持续记录数据，首先需要对航迹数据进行测线识别、裁剪，目前，许多国内、外著名的地球物理数据处理软件(如Montaj Oasis、geoProbe等)，在此环节上，处理起来并不十分“舒畅”。笔者利用国内地质行业最常用的MapGis地理信息系统，对航测数据实现可视化的数据编辑与筛选，不仅适用于无人值守飞行数据，而且对于常规飞行数据，同样具有针对性强、工作效率高的特点。

MapGis地理信息系统是武汉中地信息工程有限公司开发的GIS基础平台软件系统，集地图输入、数据库管理及空间分析的综合系统；

不仅图形编辑功能强大，而且具有较完善的二次开发功能[2]；是国内地质、物探领域中应用最广泛的软件之一。属性管理子系统是MapGis的重要组成部分，可实现对属性数据管理[3]。本程序是利用MapGis图元的属性结构(包括点、线和区属性结构)，通过录入、修改属性数据，借助MapGis强大的二次开发和图形编辑功能，交互式地进行数据的识别、筛选，最终获得符合质量要求的测量数据。

1 设计背景

1.1 航空物探数据简介

我中心承担的航测任务中，大多以航磁、航放同步测量为主。航空数据采集是以飞行架次为单位。数据记录为一个二进制文件，先将二进制文件转换为3个ASCII文件(航磁、航放窗数据、航放全谱)；每个测量数据均包括文件头和测量数据两个部分；文件头主要记录仪器型号、参数设置、测量时间、地点、架次等参数。测量数据记录了该架次所有实测数据，包括早晚测试(能谱)、早、晚基线(磁、放)、测线(磁、放)等测量数据。图1为无人机航放窗数据ASCII文件的数据格式。从图1中可看出，测量数据包括测线号、基点号、经纬度、坐标值、高度值、温度值、放射性元素测量值、日期和时间等内容。

1.2 数据质量控制

航测数据的采集受诸多因素的影响，在确保天气、仪器性能、地面磁日变站等正常前提下，通过测区范围、飞行速度、离地飞行高度、偏航距等指标衡量数据质量；在物探项目中，项目设计测线是野外数据采集最基本的依据。离地飞行高度(离地高度)，在航空伽玛能谱测量中是一项十分重要的参数。一般情况下，飞机应沿地形起伏飞行，高度保持在80 m～ 120 m。由于放射性计数率随高度呈指数衰减，离地高度过高，导致计数率偏低，将会使地面弱异常被遗漏，同时增大统计涨落误差；离地高度过低，无法将测线间区域全覆盖[4]，同时飞行安全指数大大降低。偏航距指实际飞行航迹与设计测线的偏移距离，依据测量比例尺确定最大偏航距[5]。在实际飞行测量中，一条测线结束后，飞机转入下一条测线时，测量数据虽处于项目设计之外，但在实际数据处理时，应尽量多地保留测区外符合质量要求的数据。在测区中，受高大建筑物、地形等因素影响，往往也会出现“超高”、“偏离”的现象，此时必须及时了解具体情况，做出下一步飞行方案。

2 程序实现

为了尽可能地保留测区外的有效数据，同时及时准确掌握各架次的飞行质量，可视化的数据筛选显得十分必要。

2.1 设计流程

总体上，利用此程序进行数据编辑遵循“数据—图形—数据”变换过程，分为5个步骤，除第4步为手工操作外，其余均由程序自动完成(图2)。

1)设计测线绘制。根据设计坐标，生成具有属性结构的点、线、区三个文件，将各自设计测线信息录入其中，作为自动识别依据。

图2 程序设计流程图Fig.2 The flow chart of program

图3 无人机某一架次物探作业航迹图Fig.3 Map of a sortie of flight track of airborne geophysical survey using drone

2)航迹线的绘制。将实测数据绘制成图，线文件为航迹线，包含“测线编号”的属性。点文件包括：①子图点文件，表明“测点位置”和“高度等级”；②注释点文件，包括“测线编号”、“飞行高度”信息(图3)。

3)依据设计测线，将航迹线自动分离、识别并补充录入“偏航距”信息。

4)在MapGis中，依据测点信息(点文件)，手工编辑航迹线，剔除不合格数据。

5)依据编辑后的航迹线端点坐标，筛选实测数据。

6)换名存储航迹线，作为下一架次数据编辑的参照。

2.2 程序实现

此程序基于MapGis软件提供的COM组件[6]，使用C#语言编写。有关MapGis的二次开发环境、组件注册与引用以及点、线、区的写入等技术[7-8]。笔者以2015年在新疆喀什地区开展的无人机航空物探综合站测量技术研发与应用示范项目[9-10]中的航测数据为依据，得到简述程序实现过程。

2. 2.1 设计测线成图

1)准备工作。在程序运行前，建立带有属性结构的点、线、区模板文件。由于该程序设计仅针对于飞行质量，并不对测量值进行修正、处理，所以线、区文件只添加“测线编号”属性；点文件添加“基点号”和“测线编号”两项属性(图4)，其他数值未写入属性结构中。上述三个模板文件均为空文件，随后生成的MapGis文件均基于此模板完成[3]。设置方法可参见参考文献[3]。

为提高运算效率，坐标位置也未写入属性结构中，而是通过图中绘制的测线、测点的坐标值换算获得。

2)设计测线的绘制。根据基线、设计测线的端点坐标，在MapGis中绘制设计测线图(图5)，包括测线编号(点文件)、设计测线(线文件)、测线识别范围(区文件)。测线识别范围是指单条设计测线两侧外扩接近1/2线距，两端外扩1 km～ 2 km的区域。由于线、区文件均是基于上述模板文件生成，所以它们中的图元均具有“测线编号”的属性结构。线文件，属性录入代码为：

myLinArea.Load(modelWL);//加载线模板文件

Record att2 = new Record(); //属性定义并实例化

myLinArea.att.Get(LineCounts, out att2); //获取原有属性

att2.Value[2] = LineNo; //获取属性(线号)

myLinArea.att.Write(LineCounts, att2);//写入属性

……

//另存线文件

myLinArea.Save(WorkArea23 + @"设计测线.wl");

//区属性定义与线属性类似，(略)

另外，基于测区范围外1 km～2 km，小于飞机转弯区绘制一个范围框，称之为“航迹分离框”，用于实测测线的分离(图4)。

图4 点、线、区文件属性结构参数Fig.4 Parameters of attribute structure of point, line and region files(a)点文件属性结构;(b)线文件属性结构;(c)区文件属性结

图5 在MapGis中航迹编辑界面图Fig.5 The edit interface of flight track on MapGis

2.2.2 实测数据成图

本次测量中，航磁测量采样率为10 Hz，航放测量为1 Hz[9]。为了提高运算效率，优先选用数据量较小的能谱数据，利用线号、基点号、坐标位置(x，y)、离地高度等数据列，绘制具有属性结构的测点(子图)、注释点、航迹(线)(图3、图5)。依据离地高度等级，将测点和注释点赋予不同的颜色，直观显示数据质量。测点及注释点的绘制代码如下：

private void writingCedian(PntArea PntAreaCedian, PntArea PntAreaCDnote, DotInfo5 myDot2, Pnt_Info PntInfoNor, Pnt_Info PntInfoNote) //测点、注释点子函数

{

double fontSize = Convert.ToDouble(textBox11.Text);

Record att3 = new Record();//属性定义并实例化

D_Dot pntXY3cd = new D_Dot();

pntXY3cd.x = myDot2.x; //测点X坐标

pntXY3cd.y = myDot2.y; //测点Y坐标

D_Dot pntXY3cdAlt = new D_Dot();

pntXY3cdAlt.x = myDot2.x+1; //注释点X坐标

pntXY3cdAlt.y = myDot2.y - fontSize*0.75; //注释点X坐标

PntAreaCedian.Append(pntXY3cd, "a", PntInfoNor);

int DotCounts = PntAreaCedian.count - 1;

PntAreaCedian.att.Get(DotCounts, out att3); //取点属性

att3.Value[1] = myDot2.fn;//基点号赋值

att3.Value[2] = myDot2.lineNo;//测线号赋值

att3.Value[3] = myDot2.ralt;//雷达高度号赋值

PntAreaCedian.att.Write(LineCounts, att3);//写入属性值

//写入测点

PntAreaCDnote.Append(pntXY3cd,myDot2.fn,PntInfoNote);

//写入注释点,高度值,偏航距为零,测线识别后,再写入

PntAreaCDnote.Append(pntXY3cdAlt,myDot2.ralt.ToString() + "-" + myDot2.dis, PntInfoNoteNor);

}

程序中myDot2、PntInfoNor和PntInfoNote分别为测点、注释点的参数值，以不同的点参数值标示飞行质量等级。

2.2.3 测线识别

首先，利用航迹分离框将本架次航迹线分离为独立的航迹线段(图5)。代码如下：

double m_Radiu = 0.0001;//模糊半径

MAPGISBASCOM1Lib.IAnalysis pClip = new MAPGISBASCOM1Lib.Analysis();

//myFrame-分离框rawLinArea-原始线文件desLinArea-结果线文件

pClip.ClipLin(myFrame, rawLinArea, desLinArea, m_Radiu, Enum_Clip_Type.gisOVLY_INCLIP);

然后，遍历各航迹线段，利用测线识别范围(区文件)，逐一赋予“测线编号”属性。传递属性关键代码如下：

for (int intK = 1; intK < myLinArea.count; intK++)

{ //遍历所有实测测线

//获取实测测线上的点集合

int i = myLinArea.Get(intK, out dataSet, out inf, out dima);

…….(在测线寻找3个合适的点坐标，代码较长，略)

myLinArea.att.Get(intK, out att2); //获取实测测线属性

for (int intN = 1; intN < myRegArea.count; intN++)

{ //遍历所有设计测线范围

//计算实测测点与设计测线范围的最大距离

double minDis1 = myRegArea.MaxDistOfPntToReg(dot1,intN); double minDis2 =myRegArea.MaxDistOfPntToReg(dot2,intN);

if (minDis1 == 0 && minDis2 == 0)

{ //最大距离为0时，表明测线在设计测线范围内

//获取设计测线范围的属性，即线号

myRegArea.RegAtt.Get(intN, out myRecord);

string LineNo = myRecord.Value[3]; //赋值

att2.Value[2] = LineNo; // 传递给线属性(线号)

myLinArea.att.Write(intK, att2); // 录入线属性

break;

}

}

}

最后，根据设计测线位置，计算每个航迹点的偏航距并绘制更新偏航距的标示(图5)。

2.2.4 编辑航迹

借助MapGis软件平台，使用不同的颜色、符号标示出每个测点数据质量等级，将实测数据以图形展现，可以很清晰了解本架次的数据质量。参照以往架次的航迹线和测点标识，利用MapGis图形编辑功能，逐一修剪本架次航迹线，确定需保留数据范围。

2.2.5 数据筛选与输出

数据编辑基于MapGis中确定的航迹线端点，从源数据文件中提取质量合格的数据。

读取每条航迹线的端点坐标和航迹点(点文件)属性结构，从航迹点中获取“测线编号”和“基点号”。以“基点号”为依据，比对、筛选航磁、航放源数据，更换“测线编号”，获得最终数据。从航迹线中提取信息的代码如下：

for (int intL = 1; intL < myLinArea3.count; intL++)

{

//获取图元位置

myLinArea3.Get(intL, out lineSet, out linA, out dim);

myLinArea3.att.Get(intL, out attLine); //获取图元属性

for (int intI = 0; intI < lineSet.count; intI++)

{ //遍历线图元节点

tempXY.x = Convert.ToDouble (lineSet[intI].x.ToString ("0.00"));//提取X坐标

tempXY.y = Convert.ToDouble (lineSet[intI].y.ToString ("0.00"));//提取Y坐标

tempXY.lineNo =attLine.Value[2]; //提取线号

lineLinInfo.Add(tempXY); //写入泛型集合

}

}

从航迹点中提取信息的代码如下：

for (int intK = 1; intK < myPntArea3.count; intK++)

{ //遍历点图元

myPntArea3.att.Get(intK, out myAtt);//获取图元属性

myPntArea3.GetPos(intK, out myXY); //获取图元位置

//提取X坐标

myDot1.x = Convert.ToDouble(myXY.x.ToString("0.00"));

//提取Y坐标

myDot1.y = Convert.ToDouble(myXY.y.ToString("0.00"));

myDot1.lineNo =myAtt.Value[2]; //提取线号

myDot1.Fn =myAtt.Value[1]; //提取基点号

//提取离地高度值

myDot1.ralt = Convert.ToDouble(myAtt.Value[3]);

myDot2.Add(myDot1); //写入泛型集合

}

图6为数据筛选后的结果数据。从图6中可以看出，数据格式未变，测线编号已将原来的“005”变为了实际测线号，并删除了不符合质量要求的数据。

另外，该程序还进行了本架次的数据质量、工作量等信息计算、统计。

将本架次的航迹线更换不同颜色，以飞行日期和架次命名保存此次航迹线，供下一架次的数据裁截做参照对比。

2.3 有人值守航测数据的处理

该程序同样适用于有人值守航测数据筛选。在正常情况下，由于每个架次的测线编号已事先录入，可直接手工编辑航迹线，再根据线号、基点号比对，筛选实测数据。如若在仪器发生故障、操作员记录错误等情况下，容易圈定错误范围，利用航迹线属性修改、图形编辑筛选数据，可明显提高工作效率。

3 结语

MapGis软件是国内地质行业使用范围最广的软件之一，基于其属性结构进行二次开发，具有开发效率高，操作简单快捷等优点。利用程序提供的图形界面，交互式地编辑图形，获取合格的实测数据，不仅提高了数据裁切的针对性，而且可以熟知飞行数据质量和工作进度。同时，在数据截取后，对飞行数据质量进行统计，避免了重复工作，提高了工效。此程序已应用于新疆喀什地区无人机试验项目(2015年)、黑龙江省完达山—太平岭地区航测项目(2016年)和黑龙江省龙江—嫩江地区航测项目(2016年)中，经不断完善，实用效果较好。