0.   引言

对地质、水利水电工程、矿山、水文、地下工程、地下水和城市等众多应用领域都亟需对地下的地质结构进行精确的几何表达^[1]，由于地下钻孔采样费用昂贵，通常需要根据稀疏的钻孔数据来快速准确重构的地下三维模型，因而一种准确反映地质构造情况又支持快速更新的建模方法，成为解决这些问题的关键^[2]

围绕三维地质建模这个核心问题，大量而有成效的研究集中在利用四面体、八叉树、栅格和广义三棱柱等体元相邻无缝组合在一起来表达地质体的几何形状^[3]。广义三棱柱体元专门为地质构造设计的，直接采用钻孔数据构建层状地质体三维地质模型，保留了3个钻孔之间的内在联系，可以区分地层的顺层和穿层方向，这让广义三棱柱在三维地质建模方法中有着独特的优势，因而适合于需要快速简单，局部频繁更新和大区域的地质建模^[4]。

吴立新等研究了广义三棱柱体元模型的拓扑关系类型和表达方法^[5]；齐安文等^[6]、陈学习等^[7]、Li等^[8]、李青元等^[9]、车德福等^[10]、Wang等^[11]研究了包含复杂地质构造（断层、褶皱、透镜体，缺失层）的广义三棱柱体元地质模型的构建方法，并提出了多种解决方案；林冰仙等研究了广义三棱柱模型中虚拟钻孔的添加方法，并用于三维地质建模中去。但是目前相关理论研究还比较分散，大多关注于某一个应用的背景问题。在解决利用广义三棱柱体元在工程应用时，钻孔数据稀疏，模型精度低的问题，还有待深入研究。这种方式更具备普适性，目前直接基于广义三棱柱体元构建地质模型的精度取决于钻孔采样的分布密度且无法直接应用地统计插值方法^[12]，当钻孔采样间距大、高程变化大时，广义三棱柱模型地层界面变化不平滑，不能满足实际应用的一些需求，实用价值大大降低。

本课题研究广义三棱柱三维地质模型的插值与平滑的理论和方法，基于钻孔构建的广义三棱柱体元模型，建立一种自适应插值方法得到精细化广义三棱柱三维地质建模的理论体系框架和数学方法，实现对现有层面粗糙、生硬变化、钻孔相连构建的广义三棱柱地质模型进行插值和平滑处理，在不增加新的数据和不改变三棱柱体元数据模型的情况下，得到一个更精确和平滑变化的广义三棱柱三维地质模型。

1.   三棱柱平滑度定义

1.1   定义

针对广义三棱柱的体元几何特征，本文拟定研究将三棱柱模型进行自动精细化为主题，本文提出了一种自适应三棱柱插值的方法，这个关键是评估模型的哪些区域变化粗糙需要插值，这在传统的三棱柱构建方法是制约方法更多应用的瓶颈，如图 1所示。

图  1  广义三棱柱地质模型

Figure  1.  GTP model built by boreholes

下载: 全尺寸图片 幻灯片

三棱柱的平滑度就定义为三棱柱地层相对于周围的高程起伏变化的大小，因而评估地层面的平滑度是实现自适应插值的关键，本文针对广义三棱柱构建的三维地质模型提出基于兼顾顶底面平滑度的三棱柱插值方法。广义三棱柱的平滑度尝试将计算机几何造型中的曲面几何曲率引入到广义三棱柱的曲率计算中，研究计划引入高斯曲率作为主要突破点，高斯曲率的几何意义，即球面上的面积/曲面局部面积的极限，高斯曲率确实反映了曲面局部的弯曲程度。主曲率定义为最大曲率K_max，垂直于极大曲率面的曲率为极小值，最小曲率K_min，K表达为曲面的第二基本型和第一基本型之比是依赖曲面z=r(u, v)的点（u, v）和在该点的切方向du/dv的函数，称为曲面在该点沿切方向du/dv的法曲率K表达为

高斯曲率反映当前点的弯曲程度，定义为

式中：$L = {r_{{\text{uu}}}}n$；$N={r}_{\text{vv}}n$；$M={r}_{\text{uv}}n$；$E={r}_{\text{u}}{n}_{\text{u}}$；$F={r}_{\text{u}}{n}_{\text{v}}$；$F={r}_{\text{u}}{n}_{\text{v}}$；$F={r}_{\text{v}}{n}_{\text{v}}$；r是曲面的偏微分。E，F，G是曲面的第一不变量，L，M，N为曲面的第二不变量。

因此在三维CAD软件中都把高斯曲率分析作为分析曲面造型，当曲面的高斯曲率变化比较大比较快的时候表明曲面内部变化比较大也就意味这曲面的光滑程度越低。

高斯曲率可以反映曲面的局部起伏程度，尝试利用Taubin方法给出面积加权和的计算方法计算${N_{{v_i}}}$。如图 2所示，那么可以定义广义三棱柱的顶面和底面的平滑度为

图  2  广义三棱柱平滑度

Figure  2.  Smoothness of GTP

下载: 全尺寸图片 幻灯片

一个广义三棱柱是否平滑和精细，决定了其是否需要继续细分和插值，这取决于顶面曲率N₁、底面曲率N₂和三维Delaunay法则的遵循度来确定的，用G₃来计算，计算方法为，设定分别对应的权值为w₁，w₂，w₃，那么定义单一三棱柱的整体平滑度可以计算为

式中：n为当前钻孔的层数；h为当前三棱柱的垂直高度；i为迭代三棱柱所有层数的当前层（0 < i < n）。

1.2   计算方法

整体的广义三棱柱模型的插值，单一广义三棱柱的插值计算主要是用于评估一个三棱柱和相邻三棱柱的之间的平滑连续程度，确定是否需要进行广义三棱柱插值，通过计算广义三棱柱的平滑度，设定一个插值阈值，满足阈值的情况下，就自动进行曲率计算，然后，可以实现针对三维地质模型的自适应插值计算。

（1）迭代顶面和高斯曲率计算

从模型的边缘开始，按照拓扑关系迭代所有的三棱柱，针对每个三棱柱进行迭代计算。首先采用高斯曲率的方式迭代所有三棱柱计算顶底面的曲率，然后针对每个三棱柱计算每一层的子三棱柱的曲率。

（2）广义三棱柱平滑度计算

根据广义三棱柱平滑度的定义迭代计算所有三棱柱的平滑度。根据地质建模精度的要求设定每一个三棱柱的插值分类阈值，针对不同的精细度，设定不同的插值参数，控制一个三棱柱的细分程度和细分策略方法，满足实际需要。

（3）特殊构造处理

在地质构造的边界上的不平滑和不连续，不能参与插值和计算，因而需要提前进行标记，主要包括：断层两侧、缺失层的内部边缘和尖灭等内部边界。

2.   虚拟钻孔设计

2.1   设计方法

在上节研究了地质模型的精细度的计算方法，按照阈值控制和模型的精细度要求，可以通过迭代计算得到需要进行插值的区域，本文需要研究一种基于精细度控制的模型精细化方法。①设定计算精细度w，确定是否需要插值细分。②根据细分个数和计算方法，确定虚拟钻孔的位置。③虚拟钻孔的计算，主要是针对钻孔和地层相交的层位点的高程计算，计算方法主要利用已有的均分法。④特殊地质构造边界处理，针对断层、褶皱等特殊地质构造，按照边界点插入虚拟钻孔完成特殊构造处理。

2.2   三棱柱细分方法

根据上一节计算后的平滑度值，要实现针对三棱柱模型的平滑化和精细化，采取的细分方法和对应的光滑度与插值法制紧密相关，假定N_GTP的阈值设定为N，对应的细分的三棱柱个数为n，得到的精度为c，每一种精度选择不同的最优方法。

（1）迭代内部分解，兼顾已有的三棱柱三角形形状，选择在三角形的重心位置插入点后裂解为3个子三棱柱，可以在顶面和底面保证曲面更平滑。通过在内部子三棱柱再插入虚拟钻孔的方式进行多次细分。

（2）边缘分解，按照均分插入钻孔的方式，参考三角形的均匀细分三角形的构建方法，二次细分和三次细分的位置计算采用Delaunay规则细分法。

（3）Delaunay规则细分法，在插入新钻孔的几种方案中，针对每个面四边形的Delaunay规则符合度进行比较，选择最优作为最终方案。

2.3   虚拟钻孔计算

虚拟钻孔的计算方法主要包括虚拟钻孔的位置计算和钻孔的层位计算，在上一节中得到了虚拟钻孔插入的个数和位置。虚拟钻孔的层位计算方法采用伯恩斯坦-贝塞尔Bernstein-Bézier二次三角形的计算方法得到。这种方式在局部插值中类似于样条曲面的三角化扩展，计算速度快，并且计算精度稳定。N次三角形有(n+1)(n+2)/2个控制点${v_{1, 2 \cdots n}}$，将三角形的三边分别进行n等分，连接平行于各边的等分线段，得到一个三角形形成的三角网格，这个网格就是控制网格或者叫Bézier网格，那么根据Bernstein-Bézier多项式理论可知，n次Bernstein多项式是线性无关的, 它们组成的n次多项式，因此也称为P_n的Bernstein基函数。并且满足单位分解性，

所有n次Bernstein多项式组成的行向量记为Bⁿ，重新换算到笛卡尔直角坐标系中，对于任意n元二次多项式有，

3.   自适应精细化建模方法

具体流程如图 3。

图  3  精细化三维地质模型构建方法

Figure  3.  Building method for smooth GTP model

下载: 全尺寸图片 幻灯片

（1）广义三棱柱地质模型构建，根据原始钻孔构建广义三棱柱几何模型，对于复杂地质模型通常有很多不整合情况，包括透镜体、缺失层、断层等情况，需要保留缺失地层的广义三棱柱。

（2）平滑度计算，根据提出的平滑度研究方案，计算一个广义三棱柱的平滑度，然后按照顺序从上到下迭代三个相同钻孔构建的其他层对应的广义三棱柱，分别计算平滑度。

（3）自适应添加虚拟钻孔，首先，确定预先设置虚拟钻孔的插值停止阀值，这个数值根据模型的精度要求、数据量、计算时间事先设定一个合理值，需要根据多次试验得到。虚拟钻孔的位置分三棱柱的中心或者重心和三角形的Bernstein-Bézier细分方法。

（4）迭代插值过程，重新计算经过插值的初始广义三棱柱的平滑度，如果仍然大于阀值，则重复（2）和（3）两步直至满足要求为止，就完成了整个地质模型的插值过程。即可得到一个精细化的三维地质模型。

4.   应用案例

为了验证本文提出的广义三棱柱自适应插值方法的可行性，本文的研究区域是辽宁省的某个地质勘探数据，为了简化问题，只选择了部分区域数据，共有52个钻孔，如图 4，探测到了6个地层，在未测定位置和虚拟钻孔位置则选择用空间统计插值方法确定。

图  4  基于52个钻孔构建的地质曲面和模型，增加138个虚拟钻孔后的地质曲面和本文方法构建的精细化模型

Figure  4.  The GTP surface and model based on 52 borehole data, the surface based on the 138 virtual borehole added by the adaptive GTP interpolation method and the model

下载: 全尺寸图片 幻灯片

根据钻孔勘探数据，利用剖面数据和区域知识理解，在断层边界处添加5个虚拟钻孔，不整合地层边界处添加6个虚拟钻孔，利用三棱柱几何模型构建方法，构建一个的广义三棱柱地层模型，共524个独立的三棱柱，如图 4，尽可能详细的描述了地层几何限定，包括断层和不整合层等。设置计算的平滑度计算，还有25个三棱柱需要多次插值，其中利用了中心点分解法插入了138个虚拟钻孔。通过多次插值后，平滑度都满足要求，相邻三棱柱和钻孔之间过度平滑，同时三棱柱的形状都更优，证实满足Delaunay三角形法则。

5.   结语

根据本文提出的一种三维地质模型构建方法，针对三棱柱模型中精细化方法缺失的问题，提出一种基于模型平滑度定义和计算方法的三棱柱地质模型精细化构建方法，首先建立了模型精细度的计算标准，其次研究了基于设定精细度的模型自适应细分构建方法，这既考虑了传统建模方法的优势，又考虑了精细化三棱柱模型的构建需要。未来，在细分过程中的虚拟钻孔计算方法可以深入研究，从而为用户提供更少的设置参数，以便更自动化地完成建模过程。