0 引言

离散点插值计算在模型构建中常发挥重要作用。在解决特定科学问题时，所获取到的离散点数据往往是有限的，利用插值算法可以很好地解决未知点的赋值问题，如对海洋风场、地下水流场、地质实体等的模型构建。

随着可视化建模软件的广泛应用，三维地质建模技术在采矿行业、建筑地基勘查、城市规划、排水管网布设众多行业发挥了重要作用。在插值算法方面：Chen Guoxing等^[1]通过对剪切波速度随深度变化的关系研究，结合钻孔信息，利用克里金插值算法对苏州城区的地质情况进行了三维模型构建；He Zhenwen等^[2]对简单地质模型中的网格简化方法进行了探讨；Hou Huikun等^[3]基于矿区之间的平面图和剖面图数据建立了类比右三角棱镜网络模型，并将该模型应用于常村煤矿；Ítalo Gomes Gonçalves等^[4]对用于地质结构隐式建模的势场方法的机器学习算法进行了探讨；Mohamed Arfaoui等^[5]对克里格插值算子在El Kef-Ouargha地区重力面中的预测优势进行了研究。在地质建模方面：Xiong Zuqiang等^[6]通过对地质建模中数据形式和建模软件实现原理的分析，论述了三维地质建模在工程中的应用进展；Hu Jinhu等^[7]对基于钻孔数据的地质建模在工程中的应用进行了研究；Zhou Zhao等^[8]通过数据间的拓扑关系对三维地质建模进行了尝试；s.I.Bilibin等^[9]利用DV—Ge对碳酸盐储层的地质情况进行了刻画；Qu Honggang等^[10]通过对不同城市功能区进行地质分区，探讨了针对这些不同功能区的地质建模方法。在成矿机理、地质灾害防治和矿产潜力分析方面：Liu Jibo等^[11]通过VC++和OpenGL搭建了一个用于地质建模的平台，将三维地质建模应用到了矿区塌陷的预测当中；Li C L等^[12]基于主要层位控制的思路，对大庆油田进行了地质建模研究，为油田的二次开发提供了理论支撑；Wu Qiang等^[13]通过对矿区进行三维地质建模，实现了矿区的数字化。在对特定地质构造刻画方面：Saffarzadeh Sadegh等^[14]通过对地震测量地质构造的主控参数进行分析，并利用基于这些主控参数下的地震数据对研究区的地层、断裂进行了刻画；Li Xutao等^[15]对在缺乏大量数据支持情况下的碳酸盐储层三维地质实现方法进行了研究。

上述研究表明，现阶段的地质建模主要利用建模软件对地质实体进行三维表示，鲜有文献报道不同插值算法之间分析和对比的研究成果。插值算法与插值效果密切相关，选择不当会导致插值结果和实际情况存在较大偏差，影响模型精度。因此，笔者利用GMS(groundwater modeling system)软件提供的两种经典插值算法(反距离权重法和自然邻域法)对研究区进行地质建模，对比了两种插值算法的建模效果，归纳了两种插值算法的适用条件，以期为地质建模的插值算法选择提供参考。

1 插值算法基本原理

根据数据基本特点的不同，插值算法选取也会随之不同，有基于数学统计型的插值算法，也有基于几何拓扑关系的插值算法。随着对地质构造认识的深入，地质建模的插值算法也呈现出不断向精细化发展的趋势。由初期的拉格朗日插值算法、高斯插值算法、样条曲线插值到克里金插值法、反距离权重插值法、自然邻域插值法、谢别德插值法，插值算法呈现出多元化发展的趋势。现阶段出现了一些改进型的插值算法，但其基本的理论基础并没有太大变化。本文基于区域断裂空间展布特征，选取GMS软件中的反距离权重插值算法和自然邻域插值算法为三维实现的插值算法，对三维地质建模效果进行对比研究。

1.1 反距离权重插值法

反距离权重^[16-18](inverse distance weighted，IDW)是一种确定性方法，用于已知散点集的多变量插值。未知点分配值用已知点坐标的加权平均值计算。反距离权重法权重的分配与插值点到已知点距离的倒数相关。

假设预期结果是未知函数的离散赋值U，U是x和y的函数，而x，y都是在研究区D内取任意实数值的变量，其数学表达如下：

如果x和y在研究区D内分别有n个值，则有U(x₁, y₁)，U(x₂, y₂)，…，U(x_n, y_n)。为了便于表达，令U₁ =U(x₁，y₁), U₂=U(x₂，y₂), …, U_n=U(x_n，y_n)，则该插值算法的基本形式如下：

式中：(x，y)为插值点；(x_i, y_i )为已知点；d_i是从已知点到插值点的给定距离(度量算子)；λ_i为插值点权重；n为用于插值的已知点总数。权重的赋值与插值点到已知点的距离负相关，即距离较大处的已知点对插值点贡献较小，反之则对插值点的贡献较大。

由式(4)可知：当插值点与已知点重合时，已知点处的函数值为插值后的值；当插值点与已知点不重合时，则插值点处的函数值为各已知点权值与其函数值之积的求和值。

1.2 自然邻域插值法

自然邻域法^[19]可找到距查询点最近的输入样本子集，并基于区域大小按比例对这些样本应用权重来进行插值^[20]，该插值也称为Sibson或“区域占用(area-stealing)”插值。自然邻域简单来说是具有共用边的区域，而在共用边区域内的离散点为自然邻点，所有点的自然邻域都与邻近泰森(Voronoi)多边形相关(图 1a)。首先，泰森多边形由所有指定点构造而成，并由可识别的多边形表示。然后会在插值点周围创建新泰森多边形。显然A点的自然邻点有6个，分别为1，2，3，4，5，6。这些邻点可作为插值网格的节点，将这些节点按一定的规则顺次连接起来就构成了Delaunay三角网(图 1b)。新的多边形与原始多边形之间的重叠比例将被用作权重。如定义邻点坐标N_i(x)的权重为第二顺序泰森多边形面积与第一顺序泰森多边形面积之比，在图 1c中即为边afghe所围成的面积与边abcde所围成的面积的比值。从权重的角度来看，自然邻域插值法实际上是一种将权重的计算限定在最近范围内的一种插值法。当求某一未知点的预测值时，有如下数学表达式：

式中：N_i(x)为x点关于其自然邻域节点下的权重；f_i为已知点的坐标值。

由上述分析可知，该插值算法的基本属性是其具有局部性，仅使用查询点周围的样本子集，且保证插值高度不会超过样本子集内部的最大高度，对局部数据特性的继承性能较好。但是，在同一地层，自然邻域插值法必须要有4个控制点，在地层缺失严重的情况下难以实现插值计算。

2 实例研究 2.1 插值实现软件

为实现对自然邻域插值法和反距离权重插值法插值效果的对比，本文选用美国杨伯翰大学环境模拟实验室和美国陆军水路试验站合作开发的GMS三维实现软件^[22-25]。该软件包含了用于三维实现的TINs、Borehole、Solid等模块，可对获取到的钻孔地层信息进行插值，实现对三维实体模型任意角度的切割，便于对建模效果进行查看和分析。并且，GMS软件平台可以利用自身模块进行地下水运移模拟，也可以通过与其他软件的数据共享实现热量传递、污染物扩散等过程的模拟。

2.2 研究区地质背景概况

本次研究区为华北地台潘庄凸起构造单元，研究区中部为凸起区，以F₁、F₂、F₃、F₄边界断裂为界，四周为凹陷区。

受中生代燕山运动及新生代喜山运动的影响，研究区在中生代至新近纪处于隆起剥蚀环境，中生界及下古近系大多缺失。新近纪—第四纪以来，接受新生界沉积。其中：凸起区主要地层由中上元古界和古生界组成，上覆新生界，厚1 400~1 900 m，缺失古近系；周边凹陷区相对完整，仅缺失中生界部分地层。

从研究区搜集到的300余口钻孔资料看，从老到新揭示的主要地层有：蓟县系雾迷山组(Jxw)，青白口系(Qb)，寒武系(∈)，奥陶系(O)，石炭系—二叠系(C—P)，白垩系(K)，古近系沙河街组(Es)、东营组(Ed)，新近系馆陶组(Ng)、明化镇组(Nm)，第四系(Q)。

受强烈新构造运动影响，区内基岩断裂较为发育，除边界断裂外，还有发育大、小断裂20余条(图 2)。由于断裂主要形成于隆升环境，导致断裂大多为正断层。

2.3 建模实现

将GMS软件整个建模过程归纳为以下几步：

1) 钻孔编录数据的整理，按照岩性、层位等指标整理为GMS软件的输入数据。

2) 将研究区的边界信息导入模型，并利用map模块对边界进行三角网格建立。

3) 利用Borehole模块对导入的钻孔信息进行识别和计算，选取反距离权重/自然邻域插值法进行插值计算，生成地质实体。

研究区插值点的数目为TINs模块下生成的三角网的节点数，总共为82 960个，两种插值算法下的插值密度保持一致。最终生成研究区地质实体图(图 3)。从图 3中可以看出，反距离权重插值算法对地层起伏状态的描述更明显，在有断裂分布处的插值更能反映地层的起伏情况。

3 结果对比分析

本文插值效果对比研究的整体思路是，将原始数据和插值结果进行对比，阐明反距离插值法和自然邻域插值法的插值效果，并分析出各自的插值精度，总结这两种方法各自的适用条件，为地质建模插值算法的选择提供参考。在原始数据的预处理阶段，首先，对原始数据各地层的离散程度进行分析，评价钻孔分布的不均匀性；其次，对钻孔揭示的地层进行均方差分析，确定真实地层在插值前的整体起伏状态。在插值结果分析处理阶段，利用导出的地层插值数据提取每个地层的底板散点数据，并计算其均方差，进行纵向分析和横向对比。纵向分析是对各地层的插值均方差进行分析，得出地层在垂向上的起伏情况，并与断裂的情况进行对比；横向对比是对某一地层原始数据的离散程度和插值后的离散程度作对比，分析两种方法对整个研究区地层的继承情况。以对比分析结果为依据，总结得出两种插值算法的优劣及各自的适用条件。

3.1 均方差对比

因为研究区没有钻孔完全揭穿雾迷山组，对其对比分析不具有代表性，所以本研究着重考虑第四系到青白口系的地层。从图 3的建模结果可以看出：研究区地层分布较为复杂，第四系、奥陶系和青白口系在整个研究区内都有分布，其他地层在研究区都存在不同程度的缺失，东营组分布范围较小。从表 1和图 3可看出：随着地层深度的增加，可以发现每个地层的起伏趋势在增长，以石炭系—二叠系以下地层的起伏最为显著。从两种插值算法得到的各地层的均方差来看，二者差别不大，均表现出了对原始数据的继承效果(图 4)。但从馆陶组到白垩系，由于揭露的钻孔资料较少，自然邻域插值得到的地层均方差比反距离权重插值法得到的地层均方差大，表现出了较大的地层起伏。而且，其插值点也较反距离权重插值法得到的点多，表明自然邻域模型中这几个地层的分布较广。

3.2 相对误差对比

为了更精细地对比两种插值算法的插值精度，挑选研究区比较有代表性的15个钻孔(图 2)，将钻孔位置处的插值数据和钻孔原始数据进行对比。选取了4个比较具有代表性的地层进行了相对误差分析(图 5)。图 5中的钻孔统一用数字代号1—15表示，数字1—15对应的钻孔编号依次为SR35D、DL-76B、SR37D、SR34、SR22D、HD-25、HDR30、HDR-34D、HDR-32D、DL-53、DL-55B、SR13、SR-33D、DL-79B、SR-21；曲线没有点控制的地方表明有地层缺失。从图 5a 、b 可以看出两种插值算法对第四系和明化镇组的刻画都比较接近，均方差表现一致。在对第四系的刻画中，除了自然邻域插值法在DL-79B钻孔处产生了较大误差外，其他位置处的误差都控制在了±0.05之间(图 5a)；在对明化镇组的刻画中，除了反距离权重插值法在SR35D产生了-0.20的误差、自然邻域在SR-33D和SR13处产生了-0.20和-0.13的误差外，其他钻孔位置处的插值误差都控制在了±0.05之间(图 5b)。馆陶组(图 5c)缺失严重，插值结果的相对误差有较大的波动，尤其是自然邻域插值法的SR-33D处出现了较大误差，其误差高于70%。同时，奥陶系也存在相类似情况(图 5d)。从4个地层的相对插值误差来看，反距离权重插值法的整体误差较自然邻域插值法小，有较高的插值精度。

综上，自然邻域法在有地层缺失的层位插值结果的离散程度较大，表明对地层的插值计算结果有较大误差；反距离插值法的计算结果只在断裂附近出现地层起伏(图 3)，有地层缺失的层位的插值结果误差在可控范围内，更能反映整个地层的整体趋势。

4 结论与讨论

1) 反距离权重插值法对断裂处起伏地层的刻画更明显。

2) 反距离权重插值法具有很好的普适性，在地层完整和有缺失的情况下都能实现插值，且插值误差比自然邻域插值法更小。

3) 自然邻域法在地层有缺失的地方，插值得到的点更多，一旦研究区地层缺失太严重，势必造成较大的插值失真。

总之，反距离权重插值法几乎适用于所有情况，地层有缺失和钻孔分布极不均匀的情况都适用；而自然邻域插值法适用于数据分布均匀，某一地层至少有4个控制点的情况。在实际应用中，需要根据具体的地质条件和钻孔数据情况综合分析选择最优的插值算法。在地层较为完整，没有出现大的断裂和严重地层缺失的情况下，反距离权重和自然邻域插值法都可以取得很好的建模效果，甚至在某些方面自然邻域可以表现出优于反距离权重法的光滑性。

本文仅对两种插值算法的插值效果和对地质数据的适用性进行了比较研究。而地层数据在空间上的离散程度，以及在平面上的聚集状态等因素没做深入研究，而这些因素很可能会对插值算法中的权重产生影响，从而影响插值算法的插值精度。至于其影响程度如何，以及如何改进方法来减少数据聚集对权重的影响将有待进一步的研究。