矿体的平面投影轮廓线是矿山工程布置与采矿设计的基础，常用于矿体结构分析、模型重建、矿山开拓设计和开采计划制定等业务场景。在模型重建中，矿体轮廓线的精度直接影响重建后矿体模型的精度，进而影响矿体资源储量的精细化估算、管理以及后续基于矿体模型的采准和爆破设计。因此，研究矿体投影轮廓线对于矿山资源储量管理以及生产设计均具有重要意义。

目前，常用的国内外数字采矿软件多采用三角网格模型来表达矿体的空间位置和大小（谭正华等，2017；荆永滨等，2016；李丹等，2013）。三角网格模型是用一系列的空间三角形面片组成的网格来逼近三维实体（侯志刚，2020），随着模型的不断细化，三角形面片的数量逐渐增多，对算法速度有了更高的要求。近年来，研究人员在轮廓线提取相关算法上开展了一些研究。例如，在三维空间中对多边形进行联合求并（张华鑫等，2011；黄志等，2014）、求交（贾晓彦等，2011；姚晓等，2018）、求差（姜晓琴等，2015；王慧青等，2016）和快速融合（杨洋等，2016；阮孟贵等，2010）等处理，对三角网格曲面求交（李宁等，2013；陈振等，2016；徐敬华等，2018；史永丰等，2019）以及将三角网格投影至二维平面，通过布尔运算对三角形面片进行递归求并（姜旭东等，2016；吴付坤等，2019；Zhong et al.，2019），得到最终的轮廓线，边界线的精度能够与模型完全匹配。网格划分是该类算法常用的数据处理方法，该方法通过将空间分解成大小相同的网格，对落入相应网格中的数据添加索引（赖广陵等，2018），从而建立数据之间的空间关联，提高运行效率。上述研究在矿体轮廓线提取算法上取得了较好的成果，提取的投影轮廓线精度较高，但也存在一定的局限，如计算量大、运行效率低等，这些不足在三角形面片个数较多时尤为明显。

在矿山的实际工作中，目前主流的数字采矿软件在求取形态简单的矿体投影轮廓线时，速度尚可接受，但对于一些地质条件复杂、三角形面片多达几十万甚至上百万个的情况下，这些软件提取投影轮廓线的速度极慢。在日常设计工作中，经常需要求取投影轮廓线，因此投影轮廓线提取速度慢在一定程度上会影响工作效率。基于此，本文提出了一种基于网格划分的多边形布尔运算方法，将三维模型用投影的方式降低至二维，然后用网格分割三角形面片集合，将网格单元划分为外部单元、内部单元和边界单元3种类型，通过处理减少需要计算的三角形数量后，以网格单元为单位对三角形面片进行递归布尔运算，得到最终的投影轮廓线。该方法极大地减少了计算工作量，提高了算法的运行效率，同时保证了投影轮廓线的精度，已在DIMINE数字采矿软件中得到应用。

三角形面片的数量是决定三角集合布尔运算耗费时间的关键。对于三角网格模型，形成最终轮廓线的是位于边缘的三角形面片的边，内部的三角形面片在计算过程中只起到向外填充的作用。因此，整合内部的三角形面片可极大地减少计算量，提高运行效率。

本算法的基本原理：将网格划分成三角网格后，找出位于矿体内部的网格单元，使用大小等同于网格单元的矩形代替单元中的三角形面片进行计算。算法基本流程如下：

（1）将构成矿体空间模型的三角形集合降维投影至选定的二维平面形成二维三角集合，过程中删去因投影产生的奇异三角形。

（2）使用自适用规则网格划分平面中的三角集合。

（3）根据每个网格单元的位置及其邻接单元中的数据及外部单元、内部单元和边界单元3种类型对每个单元进行划分。

（4）使用大小等同于网格单元的矩形代替内部单元中的三角形面片参与计算。

（5）清理边界单元中的冗余三角形。

（6）以网格单元为单位递归求并，得到结果。

算法流程如图1所示。

与目前常用的投影轮廓线算法相比，本算法采用邻接单元搜索方法将网格单元进行分类，并用单元的边界矩形代替内部单元中的三角形面片参与运算，该方法有效地提高了算法速度。由于跨越了多个网格单元的三角形面片通常存在于复杂矿体模型中，因此将其划分至与之相交的三角形面片所在的边界单元中，有利于减少冗余计算，同时可为得到精度较高的投影轮廓线提供保障。但是，该算法在对含有细小孔洞的矿体模型求取投影轮廓线时，效果不太理想。

首先选择投影平面，获取平面法向方向及平面坐标原点后，使用转换矩阵将当前坐标系转换至平面坐标系；将所有三角形面片投影至平面，在平面上构建新的二维三角形面片集；计算新集合的X、Y坐标极值，得到三角形面片集合的外包矩形，并以该外包矩形为基准构建二维N*N规则网格划分区域，设三角形面片总数为S，考虑到冗余度以及内存空间占用，N的取值设置为S/10（取整），即不考虑冗余的情况下，平均每个网格单元中存储10个数据。遍历三角形面片集，将每一个三角形面片根据各自空间位置划进相应的网格单元中，具体划分方式如下：先求出三角形面片的最小外包矩形，根据最小外包矩形的坐标确定可能包含该三角形面片的网格单元，再对这些网格单元进行判断，得到划分结果。以图2为例，矿体模型以垂直于纸面的方向投影，得到二维三角集合后，以网格划分三角集合。

将网格单元划分为外部单元、内部单元和边界单元。内部单元：在轮廓线生成后，位于轮廓线内部的网格单元，即被三角形面片完全填充的单元，如图2中的A。外部单元：在轮廓线生成后，位于轮廓线外部或孔洞的网格单元，即不含三角形面片的单元，如图2中空白阴影部分。边界单元：位于实体区域和外部区域之间的区域，即含有三角形面片但不被完全填充的单元，如图2中的B、C。网格单元的分类流程如下：

（1）创建搜索因子M，M的取值将决定后续步骤中网格单元的搜索范围与数量。

（2）遍历所有网格单元，如果单元中包含三角形面片，便将其暂时划分为内部单元，如果区域中不含任何三角形面片，则划分为外部单元。

（3）上述步骤完成后，遍历内部单元，使用一种搜索邻接单元的方法对每个内部单元进行判断。如图3所示，对一个内部单元周围的邻接单元进行搜索，设要对单元A进行搜索，当M=1时，搜索周围一层单元；当M=2时，搜索周围两层单元，M值继续增大时，依此类推。

如果被搜索单元中存在外部单元，则将该内部单元重新划分为边界单元；如果被搜索单元不含外部单元，则不做改变。一些单元由于处于网格的外围，无法进行完整搜索，因此在进行邻接单元的搜索之前，应判断该内部单元是否位于网格外围。设单元的行、列索引值为i、j，网格行、列的最大最小索引值分别为I₁、I₂、J₁和J₂，外围网格的判断语句为：

if（（i+M>I₁||i-M<I₂）&&（j+M>J₁||j-M<J₂））

｛

Return true；

｝

else

｛

Return false；

｝

以图2中的网格为例，设M=1，按照上述方式对A、B、C 3个内部单元进行判断。A的被搜索单元中不含外部单元，所以A为内部单元；B右下角的一个被搜索单元为外部单元，于是B被划分为边界单元。按照公式，C位于网格外围，无法进行完整的搜索，因此C被划分为边界单元。其余内部单元按同样处理方式。

（4）在遍历完成后，在每个内部单元中添加一个大小、坐标均等同于对应单元的矩形面片，并删除其中所有的三角形面片索引。

网格单元分类的流程图如图4所示。

网格的划分一定会有三角形被分配到多个网格单元，造成重复运算。对网格单元进行分类后，可对边界单元中的三角形面片进行去重操作。具体方式为：对于一个位于边界单元，且跨越了2个以上分区的三角形面片，在其所属边界单元中寻找与其相交或相邻的另外一个三角形面片，选取二者任意一个交点的坐标位置，交点的坐标如果位于任意边界单元，就将这2个三角形面片划分至该边界单元。

以图5为例：U1、U2和U3为边界单元，Z1为内部单元。2、3、4为分属于边界单元中的3个三角形面片。1分别与2、3、4有A、B、C 3个交点，遍历周围的三角形面片，假如先找到三角形面片2，判断相交于点A，A位于U1单元中，于是将三角形面片1划分至单元U1。

算法处理主要包括降维投影、网格划分、单元的分类、冗余三角形面片清理、多边形判断相交以及布尔运算求并。设矿体模型三角形面片总数为n，则降维投影、网格划分的时间复杂度为O（n），设生成的网格单元总数为m，则单元分类的时间复杂度为O（m）；设单元分类后剩余的三角形面片数为z，矩形的数量为v，边界单元的数量为s，内部单元的数量为l，冗余三角形面片清理的效率为O（s），设布尔运算时的多边形平均边数为p，判断相交效率为O［p*（z+v）log（z+v）］，多边形求并的效率为O［（z+v）*p*logp］。因此算法的时间复杂度为O［p*（z+v）*（logp+log（z+v）］。

目前，国外一些开源算法库中也有投影轮廓线的类似算法。CGAL算法库（Computational Geometry Algorithms Library）是国外开源的C++几何算法库，提供了一系列高效、可靠的几何算法，主要包含与几何相关的数据结构和算法，例如三角网格化、Voronoi图、多边形和多面体的布尔运算等，广泛应用于几何计算相关的领域。CGAL：join是CGAL算法库中一种典型的用于求解二维布尔运算求并的投影算法，常用于投影轮廓线处理，该算法能得到高精度的投影轮廓线，但求解速度较慢，一些成熟的商业软件在求取投影轮廓线时都在使用该算法。考虑到CGAL：join算法应用较广泛，且提取的投影轮廓线精度较高，具有一定的代表性，因此，选取该算法与本文算法进行精度和求解速度的比较。在Windows平台下，用Visual Studio 2013开发工具集实现了上述2种算法。选取了一系列三角形面片数量不同的矿体模型，三角形面片个数分别为10 000、50 000、100 000、200 000、500 000以及1 000 000。在搜索因子M取不同的值时，分别研究其对投影轮廓线精度和求解速度的影响。

在上述三角形面片数量不同的矿体模型中，任意选取一个模型进行研究，这里选取三角形面片个数为50 000的模型，该模型形态简单，M分别取1、2、3、4、5，同时与CGAL：join算法提取的投影轮廓线进行比较，如图6所示。

分析图6可知，当M=1、M=2时，得到的投影轮廓线与模型的边界在拐角处不吻合，主要存在于凹角处，但随着M值的增大，不吻合的情况有所减少；当M=3、M=4、M=5时，得到的投影轮廓线基本一致，与模型边界吻合，且与CGAL：join算法得到的投影轮廓线重合。

当矿体形态简单时，本文算法和CGAL：join算法在精度上都能满足需求。为此，本文选取了形态较复杂的矿体模型进行对比试验，本文算法M取3，结果如图7所示。

分析图7可知，2种方法的结果总体上大致相同，但CGAL：join算法在局部地方会出现与模型三角形面片不吻合的情况，而在相同部位，本文方法的结果与模型三角形面片完全吻合。因此，本文算法求取复杂矿体投影轮廓线的结果较CGAL：join算法更为精确，当M=3时，投影轮廓线已具有较高的精度，能够满足矿山绝大部分的日常需求。

考虑到矿体模型中间存在夹石的情况，即模型中间会出现孔洞，而网格的大小取决于三角形面片的总数，选取同一个孔洞对应不同三角形面片个数的模型进行对比试验，此时M取相同的值，投影轮廓线结果如图8所示。

分析图8可知，当孔洞与网格的尺寸之比大于1时（理论上越大越好），CGAL：join与本文算法的结果都与模型三角形面片吻合；反之，2种方法的结果都有问题。对于本文算法，同一个孔洞，三角形面片个数越多，网格越大，当孔洞较小时，会出现网格比孔洞大的情况，即孔洞与网格的尺寸之比小于1，理论上，孔洞内至少要包含一个网格单元。因此，在该情况下，投影轮廓线的精度得不到保证。

选取上述提到的不同三角形面片个数的6个矿体模型，分别设置M=1、M=2、M=3，同时运行CGAL：join算法，提取的轮廓线和原始模型（M=3时）如图9所示。其中，N表示三角形面片的个数，图中黑色的边线为程序提取后的投影轮廓线，黄色的面片为模型，轮廓线和模型三角形面片基本完全吻合。

程序在运行时，分别记录了筛除三角形面片的数量和求解速度，如表1和表2所示。分析表1和表2可知，CGAL：join算法不筛除内部三角形面片的数量，而本算法则对内部三角形面片进行了处理，且搜索因子M越大，筛除的三角形面片越少，图9（a）中的三角形面片个数为10 000，当M为1、2、3时，筛除的三角形面片个数分别为8 233、6 519和5 093。CGAL：join算法由于不筛除内部三角形面片，提取投影轮廓线时，耗时普遍比本算法要长，随着面片数量的增加，2种算法的耗时差距更加明显。本算法中，搜索因子M越大，耗时有所增加，但差距不大，各算法效率对比如图10所示。

（1）本文算法基于网格划分方法，减少了参与运算的三角形，剩余面片的总数远小于三角网格的三角形面片，相较于普通的求并方式，该方法提高了筛除效率和求解速度，收敛稳定，速度较快，效率提升十分显著，已在DIMINE数字采矿软件中得到应用。

（2）本文算法若想获得更精确的轮廓线，需要增大M的取值。影响算法精度的部位主要是凹角较大的凹边，一般来说，当M=3时，精度已经可以满足矿山绝大部分日常使用需求。根据自适应网格的特性，当三角形面片增多时，网格数量会增多，三角形面片越多，筛除的效率越高，M值的增大对算法时间的影响也相对减少，因此算法对三角网格数量较多的模型效果会更好。

（3）算法对孔洞的处理能力取决于孔洞与网格的大小比例，关键在于孔洞部位能否完全包含至少一个网格单元。这样在单元划分后，孔洞部位才能存在外部单元，并在搜索后得到周围的边界单元。由于矿体建模的特性，大部分矿体不会含有太多过小的孔洞。但对于少部分含有细小孔洞的矿体模型，需开展进一步的研究。