基于VBA技术自动生成铀矿矿体水平投影图应用研究

2021-12-29徐志强王龙辉

世界核地质科学 2021年4期

徐志强,王龙辉

（核工业二〇八大队，内蒙古包头 014010）

我国砂岩型铀矿储量计算主要使用矿体水平投影图，现有矿体水平投影图制作，首先是利用地质技术人员提供的矿段揭穿点深度，结合测斜成果，计算得到揭穿点的坐标，把揭穿点坐标通过测绘软件CASS 7.0展绘到图形里，按照矿化信息，把工业孔、矿化孔、异常孔的符号手工添加进图形；然后，将平米铀含量、厚度、品位等参数按一定格式，通过计算机输入对应钻孔。依据基点基线位置，绘出勘探网，结合钻孔的坐标值范围生成图框；最后圈定储量范围，传统方法是根据钻孔坐标及矿化信息、勘探线位置，通过AutoCAD中画圆及offset命令画辅助线，按照矿体圈定的范围，通过辅助线的交点来确定矿体边界的拐点，并通过Pline命令手动连接矿体拐点，同时给每一个矿体加上参数标示框，完成矿体水平投影图的绘制。

在实际工作中，需及时绘制矿体水平投影图，实时更新估算资源储量，在制作期间应注意矿段划分原则，块段外推、块段圈定及图件绘制等各种事项［1］。矿体投影图不仅是为了更好的服务于工作［2］，也是进行下一步勘查、矿床评价的重要依据［3］，为储量的计算起到重要的作用［4］。在实际工作中，生产一线工作者使用CAD进行矿体水平投影图的绘制，电子的图件有助于提高效率、减少人为误差、方便后期保存及处理，但通过CAD制图依然需要手动鼠标绘制工作效率仍然不高。矿产勘查、储量核实中块段多、数据量大，容易使部分块段选漏，操作过程单一，出错机率大，容易造成眼误或手误，质量自检性差，检查过程费时费力，难以从不同角度分析矿体并绘制矿体水平投影图。因此，开发一种快速、自动生成砂岩型铀矿矿体水平投影图的程序，能够进一步提高铀矿地质人员的工作效率，对铀矿勘查工作具有重要意义。

1 应用开发环境介绍

1.1 应用开发环境的确定

程序以Microsoft office的Excel为数据平台，用Excel VBA对数据进行处理，以AutoCAD为图形平台，通过AutoCAD VBA进行图形绘制。现有铀矿生产数据的电子数据以E xcel数据格式为主，如坐标数据、测斜数据、矿体数据、铀矿矿化数据等，故在Excel平台上用Excel VBA进行数据处理即是最实用也是最直接快捷，Auto CAD图形平台图形绘制方便、使用广泛、绘制图形精度高，特别是图形面积计算精确，所以矿体水平投影图的图形平台选用Auto CAD。

1.2 Excel VBA简介

Excel软件是0ffice办公系列软件之一，具有强大的数据录入和数据处理能力［5］。虽然Excel可以满足通常大部分数据的处理需求，但要对一些数据做专业性较强的处理，就要在Excel中使用VBA。VBA诞生于VB之后，继承了VB语言的很多对象、属性和方法，语法结构与VB基本相同。通过VBA能够编写程序实现可以数据处理的自动化，极大提升数据处理效率［6］。

Excel应用程序的开发需要在Visual Basic编辑器（即VBE）中进行。VBE能够与Excel无缝结合，但VBE环境不能独立打开，必须在启动Excel或其它宿主程序后才能打开运行，使用Visual Basic编辑器可以完成VBA过程创建、VBA用户窗体创建、对象属性查看或修改以及VBA程序调试等任务。

1.3 Auto CAD介绍

AutoCAD是目前工程设计中应用最广泛的软件之一，被广泛应用于机械、建筑、管线施工、家装、地质等行业［7］。AutoCAD具有强大的二维绘图功能，如绘制平面图、剖面图、进行尺寸标注等，并可以进行二次开发。

1.3.1 Auto CAD Active X Automation技术

AutoCAD Active X Automation将AutoCAD对象显示到“外部世界”，使得许多编程语言或应用程序（例如VB、C++、C#、Java、Word VBA、Excel VBA等）可以访问和编辑CAD对象。VBA通过Auto CAD Active X Automation接口将消息发送到AutoCAD，并通过接口对AutoCAD进行控制。Auto CAD ActiveX Automation技术简化了程序设计过程，便于对AutoCAD进行二次开发，同时为其它Windows应用程序（如Excel和Word）之间的通信提供了公共平台［8］。

1.3.2 Auto CAD对象模型

Application对象位于AutoCAD对象模型最顶层，其它层状结构的对象是它的子对象。Preferences是系统参数对象，可用于AutoCAD的参数设置。Documents、Document分别是图形文件对象集合、图形文件对象，通过Document对象可实现对AutoCAD文件的新建、打开、输入、输出等操作。Layer、Layers是图层对象、图层集合对象，用于控制图层文件对象，如添加、删除图层，设置图层的线型、尺寸和颜色等。

2 基本思路与方法

2.1 基本思路

本方法流程如图1所示，利用钻孔孔口坐标表、勘探线基线起止点坐标、方向、基础孔距、基础线距通过VBA图框及勘探线计算程序得到图框、坐标及勘探网数据表，并通过VBA自动绘制到图形空间［9］；利用钻孔孔口坐标表、钻孔测斜数据表、矿体揭穿深度表计算揭穿点坐标，再结合揭穿点坐标、勘探类型、钻孔矿化信息表，将钻孔位置及矿化信息展上CAD图，并由上述数据计算矿体拐点坐标；将矿体拐点坐标上图，由程序进行初步连接，得到矿体水平投影图。

图1 基于VBA技术自动生成铀矿矿体水平投影图基本流程图Fig.1 The basic flow chart of automatically generating horizontal projection map of uranium ore body based on VBA technology

2.2 程序实现方法

2.2.1 EXCEL VBA数据预处理模块实现

2.2.1.1 矿体揭穿点计算

由测斜数据、揭穿矿体中心深度及孔口坐标

计算矿体揭穿点平面坐标。先从测斜成果中取得揭穿中点深度的上一测斜点和下一测斜点成果数据，用距离加权法计算出X和Y值，再根据孔口坐标计算出矿体揭穿中点的坐标。距离加权的插值算法（图2）为：

图2 距离加权示意图Fig.2 Diagram of distance weighting

式中：X、Y—揭穿点的坐标，X1、Y1、X2、Y2分别为揭穿点上下两个测点的坐标；D—揭穿点的深度，D1、D2分别为上下两测点的深度。

2.2.1.2 矿体水平投影图的图框和坐标格网及勘探线的计算

1）图框计算

利用Excel公式计算出所有点X、Y坐标的极大和极小值，分别在X和Y坐标的极大值、极小值坐标值再加、减一个固定值，得到图框四个角点的坐标。

2）坐标格网的计算

根据作图比例尺，计算出坐标网格的坐标间距，对极值坐标按坐标间距取整，计算出每条网格线的坐标。

3）勘探网的计算

以基点坐标为中心，根据勘探线方位角及基线方位角、孔距和线距，计算出每条勘探线的起点和终点坐标。

2.2.1.3 矿体拐点的计算

1）判断每一个工业孔在勘探线方向、基线方向共四个方向距离最近钻孔的矿化类型。

以工业孔为中心计算每个钻孔与该工业孔的距离及方位角，判断其与勘探线两个方向和基线两个方向距离最近的钻孔孔位，每个方向的方向角范围设为30°，距离按地质块段法要求，根据矿化类型判断，每个方向有五种可能：无钻孔、无矿孔、异常孔、矿化孔及工业孔。

推算算法如下：

两点的方位角及距离计算方法如下：

K＝（Y2－Y1）/（X2－X1）

如果X2－X1＜0，方位角α＝arctan（K）＋180如果X2－X1＞0，且Y2－Y1＞0，α＝arctan（K）＋180

如果X2－X1＞0，且Y2－Y1＜0，α＝arctan（K）＋360

算法的VBA代码如下：

dis=（（x1-x2）*（x1-x2）+（y1-y2）*（y1-y2））^（0.5）’两点间距离

k=（y2-y1）/（x2-x1）’求斜率

Ifx2-x1＜0 Then angle=Application.WorksheetFunction.Degrees（Atn（k））+180

Ifx2-x1＞0 Andy2-y1＞0 Then

angle = Application.Worksheet Function.Degrees（Atn（k））’求方位角

Ifx2-x1＞0 Andy2-y1＜0 Then

angle=Application.Worksheet Function.Degrees（Atn（k））+360

If angle＞K D-anglefw And angle＜K D+anglefw Then

If dis＜disKN And dis＜disfw Then……

2）根据四个方向关系计算矿体拐点坐标

根据块段法推算拐点坐标。每一个方向对应五种不同的可能，它们所对应的距离d值不同。若对应方向无钻孔，d为基本工程间距的四分之一；若对应方向最近钻孔为无矿孔，d取值为两孔间距的四分之一或基本工程间距的四分之一；若对应方向最近钻孔为矿化孔，d取值为两孔间距的二分之一或基本工程间距的二分之一；对应方向最近钻孔为工业孔，取d值为0；若距离最近钻孔为异常孔，按无矿孔处理（图3-图7，红、蓝、黄点分别表示工业孔、矿化孔、异常孔）。

图3 对应钻孔d值Fig.3 The d value assignment diagram for counterpart borehole

图4 对应无矿孔d值Fig.4 The d value assignment diagram for barren bore hole

图5 对应矿化孔d值Fig.5 The d value assignment diagram for mineralized hole

图6 对应异常孔d值Fig.6 The d value assignment diagram for abnormal borehole

图7 对应工业孔d值为0Fig.7 The d value assignment diagram for the industrial bore

然后编程由已知d1与d2计算拐点坐标（图8，红、蓝点分别表示工业孔、矿化孔），算法如下：

图8 由两个d值计算拐点坐标Fig.8 Calculating the inflection point coordinates from two d values

式中：X、Y—计算拐点坐标；Y0、X0—工业孔的坐标；ΔY、ΔX为向Y方向和X方向推算距离。

推算出相互关联的工业孔组，根据每一个工业孔四个方向的钻孔类型，把相互连接的工业孔挑选出来，并且生成一张矿体分类表，并且运用EXCEL的公式把计算出来的相应工业孔的矿体拐点坐标复制到表格相应的工业孔的行里。

2.2.2 Auto CAD VBA绘图模块实现

2.2.2.1 VBA实现钻孔和平米铀含量展图

1）图符号的制作及钻孔揭穿点的展图

按照现行铀矿规范制作铀矿化信息的图块，将其放在图形文件目录下，用矿化信息类别命名，如“工业铀矿孔.dwg”，在自动绘制时直接作为块展绘到图形空间。循环读取揭穿点表的坐标和孔号，按照表内的矿化信息读取对应的图块，将图块和孔号插入到图形空间中。

关键代码如下：

2）平米铀含量展图

在钻孔号后，需要按照特定格式标明钻孔的所属块段号码、水平厚度、平均品位、米百分含量，按照已知比例尺算得每个字的长、宽进，按照规范要求自动绘制。

关键代码如下：

2.2.2.2 VBA图框和坐标格网及勘探线的自动生成

1）图框和坐标格网的自动生成

利用EXCEL计算出的区域内X的最大及最小值和Y的最大最小值，再加上图框的增加值得到图框的四个角点，并对四个角坐标取整公里网格，绘制出图框的内外框，并根据Y的极值和X最大值，绘制图名并放到中间；根据图形比例尺绘制公里网格并标注。

关键代码如下：

2）勘探网的自动生成

根据勘探网的基点、基线方向、勘探线方向、线号距、孔号距及勘探网号绘制勘探网。首先，根据基点和基线方向绘制0号基线，再根据基点、勘探线方向和线号距计算每条勘探线的起点和终点坐标（算法见下）以及勘探线号，然后绘制勘探线，并在勘探线起点和终点标明勘探线号。

勘探线起点和终点计算方法（图9）为：

图9 勘探线端点计算示意Fig.9 Calculation of end point of exploration line

式中：X1、Y1—计算勘探线起点坐标；Y0、X0—基点的坐标，ΔY、ΔX—向Y方向和X方向推算距离；β—勘探线方位角。

勘探线自动生成的主要代码如下：

2.2.2.3 VBA实现矿体边界的初步圈定

矿体边界由拐点坐标连接构成，并且需要确定拐点连接方向及顺序。使用GIS点在线段左边判断原理算法结合两点间距离的限差进行初步圈定。GIS左手算法，已知两点A（X1，Y1）和B（X2，Y2），线外一点P（X0，Y0），（X2-X1）（Y0-Y1）-（Y2-Y1）（X0-X1）＞0则点P在线段AB的左侧（图10）。

图10 点在线段左边判断Fig.10 Judgment diagram of point on the left side of line segment

3 矿体水平投影图自动生成实例应用及与手工图对比分析

3.1 实例应用

本研究方法已经用于实际生产中，现以某铀矿详查项目为例，对本程序进行实例验证。

3.1.1 数据预处理模块的验证

3.1.1.1 揭穿中点计算

在实际的野外生产工作中，钻探施工期间，钻孔往往不是直孔，均会在一定程度上朝一定的方位进行偏移，所以孔口坐标并不能代表钻孔矿体的所在位置，一般根据孔口坐标及钻孔的偏斜方位及偏斜距计算钻孔矿体的实际坐标。本文由EXCEL表格的矿体顶板和底板数据及揭穿中点的上一点和下一点的测斜坐标数据计算出揭穿中点的坐标（图11，数据单位：m）。

图11 揭穿点计算示意图Fig.11 Schematic diagram of intersection point calculation

3.1.1.2 图框和坐标格网及勘探线的计算

在制作矿体水平投影图时，必须先计算好图框和坐标网格及勘探线，才能开展下一步工作，可以由EXCEL公式统计出X坐标、Y坐标的极大和极小值，并输入矿床的勘探网的基点坐标、勘探线方位角、孔距、线距、勘探网号等参数，点击按钮直接生成图框和勘探线（图12，数据单位：m）。

图12 图框及坐标网格计算示意图Fig.12 Schematic diagram of Frame and coordinate grid calculation

3.1.1.3 矿体拐点的计算

图框和坐标格网及勘探线生成后，便可开展下一步工作，对矿体的拐点进行计算，才能对矿体的范围进行圈定。由揭穿点坐标、勘探线方位角、基线方位角、勘探线方向基本工程间距、基线方向基本工程间距，利用VBA编写的函数推算出工业孔四个方向所对应情况（图13，数据单位：m）。再由推算的情况计算每一个工业孔的矿体拐点（图14，数据单位：m）。

图13 四个方向点判断示意图Fig.13 Schematic diagram of four directional points

图14 拐点计算示意Fig.14 Indication of inflection point calculation

3.1.2 Auto CAD绘图模块的验证

EXCEL模块对需要绘制图形数据处理好了后，通过程序的执行按扭自动生成矿体水平投影图。

3.1.2.1 钻孔和平米铀含量展图

通过程序进行计算后，对数据进行预处理，生成每个钻孔的矿体揭穿点（图15，数据单位：m），在利用揭穿点数据在Auto CAD中进行揭穿点的展布，形成揭穿点展图效果图，可以确定揭穿点的位置及钻孔的矿化类型，与以往模式一致，但速度较快，较为便捷（图16，数据单位：m）。

图15 揭穿点表意图Fig.15 List of intersection the point of intent

图16 揭穿点展图效果图Fig.16 Intersection point rendering

3.1.2.2 图框和坐标格网及勘探线的自动生成

在前期数据比较齐全的情况下，点击按钮直接生成图框和勘探线，图框和坐标网格会根据图面的大小进行适当的匹配，不需要再进行人工计算，勘探线也会根据钻孔的位置进行生成，如图17所示，图框及坐标格网基本覆盖钻孔范围，勘探线及基线也进行了自动绘制，勘探线基本为北西向南东，勘探线数量为13条，走向为135°，基线基本为北东向南西，基线为1条，走向为45°，为下一步矿体的圈定做好了准备工作。

图17 图框和坐标网格及勘探线效果图Fig.17 Overview of frame and coordinate grid of exploration line

3.1.2.3 矿体边界的初步圈定

程序可根据区内每个钻孔的实际情况，在计算出矿体拐点坐标的基础上，圈定出矿体的边界范围，对无法连接的矿体进行分割，圈定出相应的矿体块段，如图18所示，矿体基本走向为北东南向西向，且被划分为多个块段，初步圈定了矿体的边界，与人工制作图件对比基本无较大的差异，圈定的边界准确可靠。

图18 矿体初步圈定效果图Fig.18 Overview of preliminary delineated orebody

3.2 程序自动成图与手工成图对比分析

3.2.1 效率对比

矿体水平投影图制作，程序绘制效率要远比传统方法绘制的效率高，程序成图所需时间基本为15 min左右，而人工用时基本为665 min，为程序成图时间的44倍，这也体现了程序成图的优势，所需时间较短，提高了工作效率，能够对图件进行及时的更新（表1）。

表1 程序自动成图与手工成图效率对比Table 1 Comparison of efficiency between automatic mapping and manual mapping

3.2.2 质量对比

因程序是以Excel输入的数据为基础的，所以只需要检查输入数据有没有错误，程序自动生成错误率很低。而手工绘制矿体水平投影图，绘制过程避免了人为误差，错误率较高。但从矿体水平投影图的内容来看，矿体的圈定上，手工圈定的相比程序圈定的精准，程序是初步圈定，可以先给地质研究区圈定大致的矿体范围，进而精确圈定，随着程序在日后与生产实践结合逐步完善，最终一定能实现矿体精确圈连。