三维地质建模是地质数据可视化、地质空间分析和数字矿山建设的关键技术之一（张申等，2007；吴立新等，2012）。隐式建模是一种新兴的建模方法，其利用空间插值方法对地勘数据进行插值，并结合计算机三维模型可视化技术，通过多边形网格化方法得到实体表面模型（Zhong et al.，2019b；毕林等，2018）。经过20多年的发展，学者们提出了多种高效的插值方法，包括离散光滑插值法（Frank et al.，2007）、距离幂次反比法（Lu et al.，2008）、局部多项式法（Baker et al.，1975）、最邻近点法（Olivier et al.，2012）、移动最小二乘法（Fleishman et al.，2005）和径向基插值法（Cuomo et al.，2017；Skala，2017）。然而这些方法主要应用于稠密的采样数据，难以较好地适应稀疏不均匀采样且需要满足多种地质规则约束的地勘数据。

为了寻找可以融合多种地质约束的隐式建模方法，国内外学者进行了大量的研究。例如：Calcagno et al.（2008）发展了一些地质规则用于插值包含多个地质序列的几何域。Guo et al.（2018）利用Hermite变换将剖面线转换为坐标和法向量，通过剖面线约束控制地层界面的几何形态。Zhong et al.（2019a）根据钻孔数据的实际地质条件和矿化域的走向，构造几种约束规则对隐式面形状进行约束。综上所述，通过添加一些约束可获得一个具有更高可靠性的地质模型。然而，对于基于钻孔数据的地质建模，在建模之前需要进行大量人工处理，即建模所需约束数据多为地质人员手动提取（Thornton et al.，2018），地质规则约束获取难度大且需要花费地质人员大量的时间和精力，基于手动提取进行地质建模的操作为半自动化隐式建模，效率不高。因此，为了保证全自动化的隐式建模，同时充分考虑地质特征条件，减少前期地质人员的工作，需要研究自动化提取隐式建模所需地质规则约束的方法。

为了实现直接基于钻孔数据的地质体隐式建模约束规则自动构造，本文提出了一种自动化提取与量化地层特征参数的方法。首先基于原始地质钻孔数据拟合地层平面并计算平面产状，提出顾及地质规则约束的地层边界点影响半径，利用三维椭球体进行局部搜索并构建邻域钻孔群，计算地质界面特征参数来自动构造插值约束；进一步分析了地质体建模规则约束，之后运用基于径向基和克里金的2种插值方法对地层曲面进行自动插值求解，实现了三维地质模型的自动化构建。

地质体隐式建模的流程主要分为3个步骤（Wang et al.，2018；郭甲腾等，2016；Lajaunie et al.，1997）：构造地质规则约束、隐式函数求解和隐式曲面重构。隐式曲面的插值求解是建立在地质采样点集所表示的插值约束的基础上。而在隐式曲面的构建问题中，原始的采样点集最终都需要转化为相应的插值约束和建模规则。为了便于表述，本文将隐式建模所需数据集划分为两类，即插值约束和建模规则，统称为地质规则约束。

插值约束为构造隐式函数插值方程的约束条件，比如用于限制地质界面和形状的点约束、梯度约束和切向约束等。从地质建模的角度来看，插值约束由地层的产状特征来表述。

隐式建模插值所需采样点集（x_i，G_j），i,j=1，2，…，N，一般为2种类型，即点数据x_i和梯度数据G_j。点数据x_i为不同地质界面之间的分界点，梯度数据G_j为垂直于每个地质界面的极化单位矢量（McInerney et al.，2005），等效为3个偏导数∂Z（x_j）/∂x，∂Z（x_j）/∂y，∂Z（x_j）/∂z，表达形式（Aug，2004）如下：

式中：Gjx、Gjy和Gjz分别为地质空间任一点x_j处在x、y和z方向上的梯度分量。

梯度数据用来进一步约束地质模型的三维空间形态，梯度分量可由地质空间中一点处的方位信息（方位角、倾角和极性）表示，如图1所示，二者转换关系（Aug，2004）如下：

式中：D为倾角，即水平线与最大倾线之间的角度；A为方位角，即最大倾线相对于磁北方向的夹角；Polarite为极性（+1或-1）信息，用来确保方向的一致性，当沉积层极性为正，代表从老到新的方向；当侵入体极性为正，代表从外部到内部的侵入。

方位信息由地质空间任一点处的方位角、倾角和极性组成，因此按照传统地层产状的计算方法可得到空间任一点处的方位角和倾角，之后根据地层属性（沉积层或侵入体）获得相关极性，进而可转化为梯度数据并构建相应的插值约束。

建模规则为用于构建完整三维实体模型各个子势场间的组合规则。从地质建模的角度来看，一般采用地层序列和地层关系来表示完整地质体各个地质界面间的组合规则。建模规则是限制多个地质界面间延展趋势的主要约束条件。

Calcagno et al.（2008）定义了3种基本的势场组合规则来表示地层序列。①具有局部平行特征的地质界面可以分组为一个地质序列；②在地层序列中根据地质年代定义地质序列所对应模型的时间顺序；③在地层序列中采用上覆（Onlap）和侵入（Erode）2种基本关系来表示地质序列间的地质事件关系，其中，被指定为上覆关系的地质序列不能改变旧的地质序列的几何形状，而被指定为侵入关系的地质序列则可以切割、截断旧的地质序列的几何形状。

地质界面特征参数是由边界点坐标和局部产状参数组成。由于边界点坐标为不同地层之间的分界点，基于原始钻孔数据可直接获取，因此，对于隐式建模所需地质规则约束的自动构建，关键是局部产状参数的提取和量化，具体方法流程如图2所示。

（1）根据原始地质钻孔数据，提取可反映全局产状的沉积层中的点坐标，利用最小二乘法拟合地层所在空间平面，任取平面上三点来计算产状。

（2）依据全局产状信息，分别沿倾向、走向和厚度3个方向计算空间变异函数，得到试验变异函数并拟合理论变异函数曲线，获得变异函数在3个方向的变程值，即地层边界点的影响半径。

（3）基于地质钻孔数据分别提取各岩性段的上下分界点，即边界点。

（4） 利用三维椭球体进行局部搜索，构建邻域钻孔群并按照不同的策略提取边界点集。

（5）基于不同类型的边界点集分别计算地质界面上的局部产状参数，最后生成隐式建模所需的地质规则约束。

基于原始钻孔数据进行局部产状参数的提取与量化，需利用局部椭球体搜索顾及地质规则的地层边界点。该过程存在2个问题：①地理空间数据存在各向异性，椭球体搜索方向的设置在地质空间中各个方向上对空间岩性的影响不同，随意选取研究方向可能导致所选方向的岩性不具备代表全局各向异性的分布规律，无法反映地层整体的赋存特征。②椭球体半径设置过大，会导致选取的钻孔距离较大，提取的钻孔边界点之间不具备相关性，无法反映地层局部变化情况，进而影响所建模型的准确度。

因此，为解决上述问题，需通过确定地层全局产状信息，获得地层整体分布情况；然后沿倾向、走向和厚度3个方向计算空间变异函数获取变程值，从而得到局部椭球体搜索方向及半径。

地层全局产状的确定一般有2种方法：一是通过原始钻孔勘探信息人为指定地层走向，即地层走向方向与勘探线方向垂直；二是通过拟合研究区地层平面，计算平面产状来确定整个地层的产状。第2种方法的具体步骤如下：

①地层平面拟合。地层平面拟合是依据原始的地质钻孔中可代表全局产状的沉积层起止点坐标，提取沉积层中点坐标（x_i，y_i，z_i），i=0，1，2，…，n-1，对于要拟合的地层平面ax+by+cz+d=0，移项后可得z=a₀x+a₁y+a₂，对一系列的n个点（x_i，y_i，z_i），只要使S=∑i=0n-1(a0x+a1y+a2-z)2最小，即可求出a₀，a₁，a₂，最后获得拟合的地层平面。

②地层产状计算。在计算地层产状之前，需要任取平面内三点坐标，计算方法为通过平面内三点坐标求该面产状（董兆岗，2000）。

为了获取地层边界点影响半径，需计算空间变异函数来获取变程值。由于精确到3个研究方向上的岩性数据点较少，需通过离散化钻孔数据，获得足够多的数据点来计算试验变异函数值。具体步骤如下：①钻孔各岩性段空间离散。在各岩性段起止点之间生成更多的点，以每个岩性段的起点坐标为第一个计算点，每间隔0.5 m生成一个新的数据点，最后不足0.5 m按照实际距离赋值，计算新增点的三维坐标，构建三维钻孔离散点集。②试验变异函数的拟合。基于钻孔各岩性段离散点，对步骤①中的离散点集进行拟合，通过变异函数计算得到3个方向的变程值即为地层边界点的影响半径。

根据第2.1小节计算的地层全局产状及边界点影响半径，通过局部椭球体搜索满足每个钻孔邻近范围内的所有钻孔构建邻域钻孔群，之后确定邻域钻孔群内所有钻孔之间的拓扑关系，根据不同的策略获得具有相关性的边界点，进而构建邻域边界点集。该方法具体流程如下：

1）提取每个钻孔沉积层的下边界点和侵入体的上下边界点坐标构成钻孔点数据集｛A₁，A₂，…，A_n｝，根据各个钻孔的开口位置，将其投影到XY二维平面上，如图3所示。

2）设提取的原始钻孔点坐标为A（x，y，z），ω为方位角，对原始钻孔点数据集｛A₁，A₂，…，A_n｝按照顺时针的方向进行θ度的旋转，其中，θ=180-ω，将数据旋转到新的坐标系下，得到新的钻孔点数据集P。

3）遍历新的钻孔点数据集P中每条钻孔上的开口点，提取其坐标y作为特征值。比较特征值y的大小并按照从大到小的顺序对所有钻孔进行排序编号，若y值大小相等则从大到小比较x值，之后构建钻孔数据集｛B₁，B₂，…，B_n｝，如图4所示。对每个钻孔按照编号顺序进行局部椭球体搜索，椭球体搜索方向和半径由地层产状及变程值确定。定义与当前钻孔B_i通过局部椭球体搜索得到的所有钻孔｛B_i，B_i+1，…，B_n｝，i=1，2，…，n，即为当前钻孔的邻域钻孔群，据此按照构网关系建立邻域钻孔搜索树。

4）遍历数据集｛B₁，B₂，…，B_n｝中的每条钻孔，对于当前钻孔B_i，邻域钻孔集为｛B_m｝，判断邻域钻孔集B_m中的钻孔数m的大小，根据m的大小分别采取单孔（m=1）、双孔（m=2）和三孔（m≥3）的策略进行邻域边界点集（O₁、O₂和O₃）的提取。算法步骤如下：

①若m=1，则满足局部椭球体搜索的钻孔只有当前钻孔一个，只进行当前钻孔B_i的边界点集提取，即提取整个钻孔的边界点S_j构成由单个边界点组成的单边界点集｛O₁｝。

②若m=2，则先进行当前钻孔B_i的单边界点集提取，再进行2个钻孔B_i和B_m之间的边界点集提取。如图5（a）所示，遍历B_i中每个边界点S_j，若该点为沉积层的下边界点，分别提取2个钻孔中同为沉积层且岩性相同的下边界点。若该点为侵入体的上边界点，首先判断该侵入体所属岩性，其次判断该岩性的上边界点在当前钻孔B_i上的数量L，若L>1，则按照z坐标从大到小在该点处依次对另一个钻孔B_m上具有相同岩性的上边界点进行沿厚度方向的搜索，提取满足椭球体沿厚度方向搜索半径h内的边界点。当该点为侵入体的下边界点时，同理得之。由此可得2个相同岩性为一组的边界点构成的双边界点集｛O₂｝。

③若m≥3，分别进行当前钻孔以及当前钻孔和其余各个钻孔之间的边界点集提取，之后进行3个钻孔边界点集提取。如图5（b）所示，进行3个钻孔边界点集提取之前，首先确定邻域钻孔群内各钻孔之间的拓扑关系，以当前钻孔B₁为中心，根据钻孔数m的大小，若m=3，则提取△B₁B₂B₃中的3个钻孔；若m>3，则按照钻孔编号排列组合的关系分别选取△B₁B₂B₃，△B₁B₂B₄，…，△B₁B_kB_m（2<k<m）中的3个钻孔；若3个钻孔在同一直线上，则不进行当前钻孔的提取。确定完3个钻孔的组合后，遍历B₁中每个边界点S_j，按照步骤②中双钻孔之间的搜索方法，分别进行当前钻孔B₁对另外2个钻孔的局部椭球体沿厚度方向的搜索，再进行另外2个钻孔之间的椭球体搜索，当3个钻孔中相同的岩性互相满足椭球体搜索半径时，则提取出同时满足3个钻孔的岩性边界点，由此获得由3个边界点为一组数据的三边界点集｛O₃｝。

5）遍历钻孔数据集｛B₁，B₂，…，B_n｝，分别获得所有钻孔的邻域边界点集，由于提取的单个、双个和3个边界点组成的点集｛O₁｝、｛O₂｝、｛O₃｝为旋转后的边界点坐标，需按逆时针旋转θ度获得复原后的点集｛O₁｝*、｛O₂｝*和｛O₃｝*。

根据第1.1小节描述的地质规则约束，分别按照插值约束和建模规则约束来构建插值整个地质体模型所需的输入数据集。

1）对于控制地层几何形态的插值约束，可通过提取的邻域边界点集｛O₁｝*、｛O₂｝*和｛O₃｝*计算各个地层的约束数据。如图6所示，假设原始钻孔为垂直钻孔，因此单边界点集｛O₁｝*在点x₁处的方位角和倾角都取0；双边界点集｛O₂｝*按照两点确定一条倾伏线，首先求出该倾伏线的产状，再以全局产状作为另一条倾伏线，由2条倾伏线获得双边界点中点x₂处的方位角和倾角；三边界点集｛O₃｝*则通过三点成面，通过三点坐标可计算平面中点处x₃的方位角和倾角；极性则按照岩性所属地层属性取值，沉积层取+1，侵入体上边界点取-1，下边界点取+1；通过方位信息和梯度数据之间的转换公式，可获得梯度约束数据，点约束数据则为不同地层之间的分界点，遍历所有的邻域边界点集进行计算，最终获得由2种约束数据构成的采样点集（x_i，G_j）。

2）对于构建整体模型的建模规则约束，主要为使用地质规则来组合多个势场：①根据地层属性，设定采样点集（x_i，G_j）中沉积层序列关系为Onlap，侵入体的序列关系为Erode。②按照地质年代顺序，钻孔中的岩性由下至上为从老到新，设定采样点集（x_i，G_j）中各个岩性所属地层的时间序列从下至上按1，2，…，N排列，即代表生成模型时从下向上按时间顺序依次构建。③对于约束数据所属地层，沉积层按从下向上赋予较前的时间序列值，侵入体在沉积层之后赋予，若侵入体某个地层中间夹有其他岩性的地层，则最后赋予夹杂地层。

在上述算法流程的基础上，分别对新疆某地下铀矿进行三维地质建模试验，在此区域的三维地质模型中共包含40条钻孔，如图7（a）所示，钻孔中包含5种不同的岩性，使用不同颜色来表示。针对原始地质钻孔数据，按照上述方法分别进行全局产状参数和局部影响半径的计算，在进行局部椭球体搜索时，沿走向、倾向和厚度方向的半径分别设定为40，40，10 m，依照上述方法获得插值约束所需邻域边界点集，如图7（b）所示，之后分别提取不同地层的采样点集，共包含546个点数据和2 041个梯度数据。

最后通过势场法和Hermite径向基插值进行单个地层曲面的求解，借助DIMINE平台进行可视化。图8（a）为泥岩层势场法和HRBF法叠加显示效果图；图8（b）为粗砂岩层势场法和HRBF法叠加显示效果图。

以上为单一地层的构建，对于复杂形态地质模型的构建，需要进行多势场组合约束。对于新疆某地下铀矿，主要含有5种岩性的地质体，分别为煤、泥岩、细砂岩、中砂岩和粗砂岩。为了对整体模型进行建模，需要按照地质年代顺序以及各地层覆盖关系对各个地层进行排列组合，之后通过开源代码Gempy（Miguel et al.，2019）进行复杂形态模型的求解。具体建模方法如下：如图9所示，首先在Gempy代码中输入提取的点数据和梯度数据，然后设置所建模型的边界大小以及地质体的精度（分辨率），按照地质年代顺序设置模型构建层序，指定各地层序列关系（表1），最后通过势场法进行插值求解。通过开源代码Gempy进行三维地质模型的求解，可获得表征地层表面的结构数据（Delaunay三角网的顶点坐标和序列号），之后把地层表面数据转换为DIMINE可识别的.off文件，在DIMINE平台进行可视化操作，建模结果如图10所示。

从实际建模效果来看，模型比较符合地质人员的认识。通过勘探线剖面，可以观察该三维地质模型各岩层的厚度和位置，模型地层尖灭和透镜体情况如图11所示。

（1）提出了一种基于钻孔数据自动化提取表征地质界面特征参数的方法。该方法兼顾地质规则及约束，以地统计学的空间变异函数计算变程值作为局部椭球体搜索半径，以保证邻域钻孔群构建的准确性，根据不同情况对领域边界点集进行提取并计算产状参数，之后构建地质约束规则，最后实现自动化的三维地质层面求解。

（2）本文方法提取的数据适用于复杂形态地质模型，在新疆某地下铀矿的三维地质建模中得到了应用，证明了该方法在实际工程中的适用性，尤其是对于尖灭和夹层数据的识别和处理，具有较好的效果。

（3）在HRBF插值求解单个地质模型以及Gempy算法构建复杂地质体模型的过程中，在输入过多提取数据的情况下，建模过程需要花费大量时间，导致建模效率过低，因此需要结合地勘数据的特性，充分考虑地质体建模的实际情况，进一步优化提取地质规则约束数据的方法。

（4）由于地下空间的复杂性，本文提出的地质体隐式建模约束规则自动构造方法具有一定的局限性；对于脉状、囊状和不规状等更为复杂的地层，局部产状参数提取并量化的难度很大，如何针对这些地层进行自动化建模，还有待进一步深入研究。