融合钻孔与地质剖面的三维地质混合插值方法

李 健，王心宇，刘沛溶，马玉荣，王广印

(1.郑州大学 地球科学与技术学院 河南 郑州 450001；2.郑州大学 水利科学与工程学院 河南 郑州 450001；3.郑州大学 图书馆 河南 郑州 450001)

0 引言

三维地质模型能够直观地展示地层的地质结构、岩性、空间拓扑等信息，不仅为专业人员分析地质构造、断层分布等工作提供准确的信息，也能为地下空间的开发利用提供可靠依据[1-5]。而三维地质空间建模的主要难点在于如何利用有限、离散的地质数据来准确表达复杂的地质体内部结构与信息[6]。在传统地质建模方法中，钻孔是最主要的数据来源[7-9]，但存在成本较高、数据稀少，无法提供大范围精确的地质构造信息等缺陷。随着技术的发展，包括物探、化探、磁法、电法在内的新兴勘测方法被广泛使用，为地质勘测数据提供了有利的补充。

为了丰富数据和提高建模精度，一些学者开展了融合多源数据的地质三维建模的研究，如罗智勇等[10]在融合钻孔与地质剖面数据的基础上，手工修正地层分界线，克服了以往单纯依靠钻孔数据建模，结果不精确且难以修正的问题，但其插值方法单一，缺少对地质数据插值方法选择的考虑。唐丙寅等[11]提出了一种基于钻孔并结合剖面数据的点→线→面→体快速递进三维地质建模方法，成功降低了地层间拓扑构建的复杂度，但在插值算法选择上仍存在一定局限性。周良辰等[12]在构建地质剖面图的环节中引入虚拟钻孔，有效地简化了算法复杂度，但该方法在大范围地质建模中效率不高且准确度难以保证。Jrgensen等[13]和Hyer等[14]则是结合了钻孔数据与航空电磁(airborne electromagnetic,AEM)数据，利用AEM获得了密集的空间信息，同时通过测量的电阻率获得了岩性信息。以上这些方法都取得了较好的效果，但数据融合过程大多基于已有的成熟软件来实现[15-17]，缺乏对图像类原始数据与钻孔数据统一格式相互融合的方法研究。另外，在对采样数据进行插值的过程中，对所有地层只使用单一插值算法，忽视了地质原理及相关制约条件对插值方法的要求和制约，无法获得最优的插值结果。

针对以上问题，本文以钻孔数据与地质剖面图数据为基础，提出了一种基于融合数据的混合插值方法。该方法可以将两种源数据融合为具有空间拓扑信息的离散点数据，提高了地质勘探资料利用率，融合过程高效透明，基于融合后的数据，通过选取各地层拟合效果最佳的算法进行插值加密，进而实现高精度三维地质模型的构建，生成的模型能够反映实际地质情况。

1 数据源融合方法

钻孔数据目前是三维地质建模的主要数据来源[18-20]，单纯依靠其得到的模型构建结果往往精度低且难以修正。而浅层地质剖面数据作为三维地质建模的另一常用数据源，获取成本低且数据相对丰富连续，可以很好地控制区域地层在水平方向的延伸，弥补钻孔数据的缺陷。但目前地质剖面图的有效利用率很低，究其主要原因，是因为图片格式无法很好地支持三维模型构建，且二维数据难以准确表达三维空间信息。如图1为浅层地质深度剖面图，横坐标表示距离，纵坐标表示深度，图中的线为地层分界线。

图1 浅层地质深度剖面图Figure 1 Shallow geologic depth profile

针对当前钻孔与地质剖面数据格式不统一、难以无缝融合的问题，本文将二者都处理成具有空间坐标和分层属性的矢量点来实现协同建模，融合过程主要分为四个步骤：提取钻孔分层节点；地质剖面图离散转化；图像配准；统一组织数据。

1.1 提取钻孔分层节点

城市地下地质构造是多期次沉积、剥蚀交替作用的结果，空间分布上有很强的不均匀性，仅基于岩性分层，其结果往往不准确，因此需要对钻孔按照年代、层序进行分层，并综合时间与岩性等属性，最终提取出钻孔分层节点[21]。

在层序分层中，首先利用控制性钻孔提取分层节点，然后对钻孔的成果进行测年工作，根据每个孔的测年成果，考虑地层岩性、颜色、形成因素，参照“河南省平原第四纪地层统一划分标准”，将每个钻孔划分为全新统(约1.2万年)、上更新统(约15万年)、中更新统(约73万年)和下更新统(约248万年)四个时期。

1.2 地质剖面图离散转化

为了解决图片格式的数据在建模中无法直接应用的问题，本文使用了一种新的剖面图处理方式，将剖面图中的地层分界信息处理为离散点，便于与钻孔数据进行融合。剖面图信息提取流程如图2所示。

图2 剖面图信息提取流程Figure 2 Section information extraction process

考虑到图片格式在计算机中的兼容性问题，需先将cdr格式转换为jpg格式，然后移除其中的标题和刻度等文字内容，以降低后续图像处理的复杂程度，导出时设置适宜每英寸点数(dots perinch,DPI)，保证剖面图精度。图像在预处理之后被抽象为二维矩阵，该矩阵与图像分辨率具有相同的行与列，对应像素的灰度值在矩阵中表示为数值。之后设定阈值对图像进行二值化，将其处理为灰度值只有0和255的二值化图像，但其地层界限的宽度经常是大于一个像素的，为了准确基于线来提取分层信息，需要对图像进行进一步细化，将线条细化为单位像素宽度。

1.3 图像配准

在完成剖面图的离散转化后，基于测线的起始点和高程对图像进行配准，统计剖面图横向和纵向上单位刻度所占像素数量，计算单个像素所代表的实际距离和深度，然后从每条地层界线上等距离提取点，提取距离可以根据建模精度需要自动调节。最后基于行列号和测线起始端点的坐标，计算提取点的坐标和高程。

1.4 统一组织数据

钻孔数据层序分层与剖面图离散转化获取的地层节点，具有相同的地质分层属性，可以直接作为三维地质建模的基础数据。经过融合处理之后的点数据，具有空间位置信息和地质分层的属性信息。为了统一管理基础数据以便后期建模，设计通用数据库模型用于存储两类点数据。

在点存储结构的设计中，不仅记录有每个点的坐标值，还使用layerId字段记录点所在的地层界面编号。为了区分两类不同的点，设置枚举字段pointType记录点的来源(钻孔/剖面)，同时通过fileId字段绑定来源文件。为了方便对文件资料的管理，本文同时设计了钻孔和剖面的存储结构。剖面存储结构记录了起始点坐标、长度和日期，钻孔存储结构记录了坐标、高程、深度、孔径和日期，两类原始文件都通过filePath字段记录其在计算机中存储的物理路径。

2 三维地质空间混合插值

空间插值的实质是通过对原始采样点的空间分布进行分析，拟合出一种函数关系使其尽可能地逼近采样点，并且通过该函数关系可以确定区域内任意位置的属性值。在三维地质建模中，通常利用插值算法对离散的采样点进行插值加密，确定地层之间的分界面，并在此基础上进行实体模型的构建。因此，插值加密的结果不仅影响地层界面的准确性，同时决定了最终实体模型的可靠性[22]。众多学者在空间插值方面进行了积极的探索，目前常用的插值方法有克里金、反距离加权和径向基函数插值[23]。其中：克里金插值方法与空间变量间的相关性有密切联系，只有当处于一定空间的数据点存在某种相关性才能使用克里金插值方法，其精确度高，可以有效弥补数据集中存在的聚类影响，但计算步骤繁琐，插值速度慢；反距离加权插值方法简便易行，但对权重函数中幂次的选择十分敏感，易受采样集群影响，造成“牛眼”效应；径向基函数插值方法具有形式简单和向高维扩展的优点，但其难以处理数据量较大的点集[24]。因此在地层模型构建过程中，需要对采样数据进行仔细研究和分析，在各地层中找到一种误差较小、适用性较高的插值方法。

在实际使用中，进行插值前需要根据工程地质条件、数据特征和数据量来灵活、合理地选择空间插值方法，从而建立精确逼真的三维地质模型。为了正确合理地对分层数据进行插值，首先需要对有限的地层分界点进行探索性空间数据分析，判断数据间的空间自相关性，然后选择合适的插值方法进行插值实现。如图3所示，本文按照以下流程对地层界面逐层进行数据插值。

2.1 探索性空间数据分析(EDA)

探索性空间数据分析(exploratory data analysis,EDA)是指对已有数据在尽量少的先验假设下通过作图、制表、方程拟合、计算特征量等手段探索数据的结构和规律的一种数据分析方法。在本文的研究中，旨在通过探索性分析发现每层采样点的数据关系，即是否具有空间自相关性，再选择适合的插值方法进行插值实现。在地学计算领域中，变异函数通常被作为描述区域化变量的主要工具[25]。当空间点x在区域上变化时，任意一方向a，相距|h|，区域化变量Z(x)在点x和(x+h)处的值Z(x)与Z(x+h)差的方差一般定义为区域化变量Z(x)在空间a方向上的变异函数，记为γ(x,h),其计算公式为

图3 插值流程图Figure 3 Interpolation flow chart

(1)

通常利用采样点及变异函数的计算公式得出样本点的实验变异函数，拟合后的曲线为经验变异函数曲线，如图4所示。

图4 变异函数曲线Figure 4 Variogram curve

图中块金效应指距离(h)为0时的变异函数值，基台值指变异函数随着距离的增加趋于稳定的最大值，变程描述了变异函数的有效范围，超出范围的变量不再具有相关关系。通过观察变异函数的曲线特征，匹配地统计学提供的理论模型的线性组合进行拟合。

本文的数据分析结果主要通过绘制变异函数云图的方式来展示，将数据集各个方向和位置对的变异函数值显示在一张图上，图中每个点代表一个数据对，横坐标表示两点之间的距离，纵坐标表示变异函数值。如果区域内的采样点具有空间自相关性质，云图会表现为明显的逐渐升高的趋势。

2.2 留一法交叉验证

为了使插值拟合结果更加科学可观，减小基于视觉效果选择方法的主观判断影响，采用留一法对插值结果分层进行交叉验证，逐地层从有N个样本的数据中选第一号样点，将其余N-1个样点作为建模样本，选用不同的插值方法拟合1号点，计算拟合值与实测值的误差，然后将1号样本放回到整个样本集，选取另一个不参与建模的样本。如此循环，直到全部样本都被用于检验一次，交叉验证的结果是每个样本点都有一个实测值和拟合值。

2.3 量化混合插值精度

在绘制折线图的基础上，选取平均相对误差和均方根误差作为最后衡量插值精度的标准。拟合值与实测值之间的误差有正有负，通常使用计算误差绝对值的平均值来消除正负的影响，以更好地反映误差的趋势性，平均相对误差数值越小，拟合效果越好。均方根误差也称标准误差，是实际值与不同拟合值误差的平方和验证点个数的比值的平方根，均方根误差越小，整体拟合精度越高。

3 实验与分析

基于上述多源数据融合混合插值方法，本文设计了一套数据处理→融合分析→插值实现的工作流程，以郑州市某区域的地质勘测资料为基础进行了试验性研究。钻孔数据取自《郑州活断层探测工作报告》，共收集33个工程地质钻孔，较浅的钻孔深度为60多米，共16个，其余钻孔深度均在100米以上。钻孔区位置示意图如图5所示，钻孔分布如图6所示。

图5 研究区位置示意图Figure 5 Location map of the study area

图6 钻孔分布示意图Figure 6 Schematic diagram of borehole distribution

8条浅层地质勘探测线，有6条测段落在研究范围内，记录有测段名称，起止点经度、纬度和测线长度等信息，原始资料中记录了对应测段的地质深度剖面。

之后提取钻孔分层节点，在收集的33个钻孔中，8号、12号、25号和34号钻孔为控制性钻孔，通过这4点匹配其他钻孔，逐一对其余29个钻孔进行地层划分。最终浅层钻孔上确定了3个地层分界点，深层钻孔上确定了5个地层分界点。钻孔分层节点提取结果如图7所示。

图7 钻孔分层结果Figure 7 Drilling layering result

随后对地质剖面图进行离散转化，将去除无关信息的图像进行二值化处理，导出时设置DPI为300，生成灰度只有0和255的二值化图像，并进一步细化。图8展示了原始矩阵→二值化→细化过程的矩阵表达，红色框表示目标轮廓。

图8 图像细化过程的矩阵表达Figure 8 Matrix expression of image thinning process

然后基于测线起始点坐标对剖面图进行配准,从每条地层上等距离提取点。图9为离散点提取结果。

图9 剖面图信息提取示意图Figure 9 Schematic diagram of profile information extraction

经过以上流程，每一幅剖面图都提取出有限个离散点，将点数据导入数据库模型中统一管理，并绘制变异函数云图来分析数据间的自相关性。图10展示了每个地层采样点的变异函数云图。

图10 各地层变异函数云图Figure 10 Cloud map of variogram

由图10可知：第1层、第2层变异函数值随着样本点距离的增加而逐渐增大，点云从低逐渐升高，符合空间自相关的云图特点。另外4个地层的云图比较散乱，没有明显的上升趋势，故这4个地层样本点之间不具有空间自相关性。

为了方便描述，6个地层分别用Z0、Z1、Z2、Z3、Z4、Z5表示，由于Z0、Z1地层具有空间自相关性，故可以使用克里金插值法。对采样点分别使用克里金插值、反距离加权插值、薄板样条和规则样条4种方法进行插值，插值结果对比如下。

Z0、Z1地层插值结果如表1所示。从插值结果可以看出，Z0、Z1两个地层总体呈现左(西)高右(东)低的地势，区域内高低起伏变化不剧烈。其中，克里金插值得到的曲面平缓，高点和低点之间能够平滑过渡；反距离加权插值得到的曲面，在高点和低点附近有明显的凸包和凹包，符合反距离插值的特点；薄板样条插值的结果整体偏高，在左上部分有轻微过渡不平滑的现象，且出现了褶皱；规则样条则对薄板样条出现褶皱这一缺点做了矫正，但是有轻微“牛眼”现象。

表1 Z0、Z1地层插值结果Table 1 Interpolation results of Z0，Z1

Z2、Z3、Z4、Z5地层插值结果如表2所示。

由插值结果可以看出，Z2、Z3、Z4、Z5 4个地层总体呈现左(西)高右(东)低的地势，区域内高低起伏变化剧烈。反距离加权插值的曲面在Z2和Z3地层上表现出了明显的“牛眼”，Z4和Z5地层则不明显，这是因为Z2和Z3地层局部出现了频繁的高低起伏；薄板样条的插值结果最为平滑，但在边角区域出现了褶皱和翘边现象；规则样条插值的曲面介于薄板样条和反距离插值之间，曲面比较平滑，但是也有局部“牛眼”现象。

表2 Z2、Z3、Z4、Z5地层插值结果Table 2 Interpolation results of Z2，Z3，Z4，Z5

从插值效率上来看，克里金插值首先要对样本数据进行探索性空间数据分析，优化变异函数的参数，同时还要顾及无偏和最优两个限定条件，计算较为复杂，插值效率低；反距离加权插值仅考虑距离的变化对估值的影响，所以效率最快；薄板样条和规则样条作为一种函数插值，需要求解待定系数，插值效率中等。

在插值结果对比之后，使用留一法对其进行交叉验证，其结果如图11所示。数据点与红色实线的距离代表误差。

图11 留一法验证结果对比图Figure 11 Comparison chart of leave-one-out verification results

从Z0和Z1地层的折线图中可以看出，克里金插值结果与实测值最为接近，而且趋势与图中实测值一致，表明了当数据具有空间自相关特点时候，基于地学统计的克里金插值要优于其他插值方法。从Z2、Z3、Z4、Z5地层的折线图可以看出，由于地形起伏变化较为剧烈，同属某一地层的采样点对应的高程相差较大，所有的插值方法与实测线都有较大的误差。总体上看，薄板样条插值的误差最为明显，规则样条插值和反距离加权插值二者的误差呈现出大致相同的趋势。

最后，选取平均相对误差和均方根误差来衡量插值精度，如表3所示，黑体数字表示均方根最小，该项插值结果最优。Z0和Z1地层选择克里金插值取得了最好的插值结果；Z3和Z5两个地层采用反距离加权插值取得了较好的结果；Z4地层的验证结果比较特殊，反距离插值结果具有最小的相对误差平均值，规则样条具有最小的均方根误差。在6个地层的插值结果中，薄板样条虽然能够得到较为平滑的插值表面，但由于过度拟合现象，使曲面偏离了实际地层，从而使精确度不如其他插值方法。

根据交叉验证的结果，分别选用误差最小的插值方法对每一个地层进行插值，Z0、Z1选择克里金插值法；Z2、Z4选择规则样条插值法；Z3、Z5选择反距离加权插值法。在插值过程中，通过规则格网的控制实现等距离插值，如图12所示构建出地层点模型。在点模型的基础上生成三角网，实现地层面模型的构建，最终结果如图13所示。

表3 交叉验证结果Table 3 Cross-validation results

图12 地层点模型Figure 12 Stratum point model

图13 地层面模型Figure 13 Stratigraphic interface model

本节针对三维地质建模最适合插值方法的选择进行了对比研究。首先利用探索性空间分析对融合后的插值离散点的相关性进行检验，发现第1、第2层具有空间自相关性的特点，可选择克里金插值法，随后结合图表与误差分析，对比了4种不同插值算法在各地层的插值效果，最终逐层选择误差最小的插值方法对离散点进行插值，有效地实现了高精度地质模型的构建。

4 结论

本文提出了一种融合钻孔与地质剖面的三维地质空间混合插值方法，该方法的主要优势有以下几个方面：① 实现了一种新的数据融合建模方式，与传统三维地质建模方法相比较，该方法充分地利用了已有的地质勘探数据，发挥了有限地质数据的价值，提高了建模精度；为了解决剖面图无法直接作为建模源数据的问题，创新性地提出了图像→矢量化→离散点的处理方式，最终将剖面图转换为离散点，完成与钻孔数据的无缝融合。② 在充分研究经典插值算法原理的基础上，结合探索性空间分析对离散点在不同地层中的分布差异进行研究，并考虑不同插值的优缺点，分别选用最适合的插值方法对每一个地层进行插值，最终生成的模型满足城市地下三维地质的高精度建模要求。

目前，本文提出的融合钻孔与地质剖面的三维地质空间混合插值方法可以在城市平原地区三维地质建模中得到应用，但还缺乏对具有复杂地质构造如断层、尖灭等情况的研究。在今后还需进一步深入该方法在这一方面的研究，使其能够适应褶皱、断层等复杂地质情况，以扩大该方法的适用范围。