0 引言

自William Smith绘制了第一张现代意义下的地质图以来^[1]，人们一直在寻求更合适和更精细的方法来表达地质现象。随着科学计算可视化的发展，从三维甚至更高维的角度来理解和重构地质对象已经成为工程地质评价、矿产资源评价、油气藏模拟等领域中定量分析的手段和方法。为了获得高精度的三维地质结构模型，需综合利用多元数据如钻孔、剖面地质图和平面地质图等^[2-5]。其中，由于能提供地质体的纵向分布、横向对比以及相互关联关系，交叉的三维剖面成为建模中的重要基础数据^[6-9]。然而，测量中的偏差、采样的有限性、认知的局限性等因素使地质剖面数据不可避免地存在误差和不确定性^[10]。已有研究表明，建模数据的不确定性会影响三维模型的几何形态和属性分布，并最终影响到各类决策^[11-12]。因此，需深入探讨和分析建模所需地质剖面数据的不确定性。

从制图的角度来看，地质剖面数据通常被抽象成地质点和地质界线数据。其中，地质界线为地质体边界线或构造线，地质点为地质界线上的节点，部分节点同时也是钻孔中的分层位置。由此可见，剖面上地质界线的不确定性决定着整个地质剖面数据的不确定性空间分布。因此，本文着重探讨地质剖面数据中地质界线位置不确定性的分析方法。在GIS(geographic information system)领域，众多学者对线实体空间不确定性做过深入探讨和研究^[13-16]。依托点类型数据不确定性分析方法，现有方法的基本思想是：假设构成线实体点的误差服从正态分布且误差相互独立，先构建点误差(椭)圆及置信体，引入协方差矩阵，再从线段两个端点导出整个线段内任意点的分布，以不同参数表征线元误差带模型，如ε-带、g-带及H带等^[13-16]，从而获得整个线实体的三维置信体及不确定性可视化度量。

然而，由于地质体特征复杂、多样，以及数据源的多元性，不同数据的误差分布有所差异；即使同一类型数据，如钻孔分层数据，其误差也有可能服从不同的分布^{[10, 17]}。可见，单一的正态误差分布假设掩盖了其他误差分布特征对模型构建的影响。建模数据的不确定性导致了地质体结构和属性模拟过程、结果的随机性。因此，用随机模拟方法分析地质数据的不确定性是一个很好的思路。本文使用Monte Carlo方法模拟剖面数据地质界线上不同概率分布(如正态分布，均匀分布)情况下的抽样采集，结合地层分界空间位置与地质属性的耦合关系，利用地质属性概率分布简化地质界线空间不确定性的可视化，探讨多种概率分布条件下剖面数据中地质界线的空间不确定性。

自Monte Carlo方法的概念提出以来^[18]，它已经发展成为一种能处理复杂问题的随机模拟方法。Monte Carlo方法是一种依据输入变量预定义的概率分布函数进行随机数值模拟的方法，主要用于模拟输入变量存在不确定性的情形以及准确获取由输入数据不确定性而带来的结果不确定性分布^[19]。主要包括两个步骤：构造概率模型和实现从已知概率分布抽样。

剖面上地质界线由多个相互连接的线段构成，其空间不确定性分布由构成地质界线的无数个离散点联合概率分布函数确定，线段间任一点的概率分布函数则由线段两端点的概率分布函数获得^{[11, 20-21]}。在计算地质界线不确定性分布之前，需要根据先验知识指定地质界线上各节点误差的分布密度函数，例如，正态分布、均匀分布和离散分布3种分布函数分别表示了连续变化的地层分界、缺失了一段的地层分界和间隔出现的地层分界的不确定性^[10]。那么，组成地质界线的节点空间不确定性可以用这些一维或二维误差分布函数表示，两个节点间的地质界线则由两个节点的连线组成。可见，剖面上地质界线的不确定性主要是由节点的误差分布密度函数决定的。

利用Monte Carlo模拟方法模拟地质界线不确定性时，需要利用计算机对节点的误差分布函数进行随机抽样生成随机变量样本值，即生成误差分布函数区间(0, 1) 上的伪随机数。在保证拓扑关系不变的前提下，每次生成的节点重新连接就生成了新的分界线。这些新的地质界线就构成了相应的不确定性。

地质剖面是三维地质建模中重要的数据源之一。它可以是地球物理的地质解释剖面，也可以是由钻孔分层信息构建的地质剖面。前者的不确定性由地球物理数据决定，后者即钻孔地质剖面既反映了原始采样数据的信息，也包含了地质专家对研究区域地质过程的认知。本论文所探讨的地质剖面即为钻孔地质剖面。

钻孔地质剖面通常由两类约束数据构成：钻孔分层数据和钻孔间地质界线。钻孔分层数据是利用钻孔取芯分析得来的信息，为硬约束数据。钻孔间地质界线是地质专家根据经验和知识绘制出来的，用于表征钻孔间的地层分界或断裂位置，为软约束数据。但是，由于地质现象的复杂性如地层尖灭等，部分钻孔间地质界线并非简单直接连接。例如，图 1a为本文研究所采用的一个三维剖面，从空间数据的组织角度来看，钻孔间的地质界线是由若干线段组成的。

a为钻孔地质剖面，圆柱形为钻孔，黑色曲线为基岩分界线，红色箭头所指为地层尖灭位置；b、c、d分别为按条件①、条件②、条件③ 模拟10 000次的结果，紫色线为给定端点分布函数后基岩线的模拟结果。 图 1 三维剖面及其基岩线不确定性Monte Carlo方法模拟结果 Figure 1 Cross-section in 3D and uncertainty simulation results of its basement rock boundaries with Monte Carlo

图选项

本论文选取了广东某地区的一条三维钻孔地质剖面作为实例数据检测方法的可行性和实用性。该剖面由不在同一平面的23个深度不一的钻孔分层信息勾勒而成，共由人工填土、灯笼沙段、万顷沙段、横栏段、杏坛段、残积层、宝月组和布心组等8个地层组成，其中宝月组和布心组为古近纪地层，其余均为第四纪地层(图 1a)。图 1a中黑色曲线为古近纪地层与第四纪地层的分界位置，在模拟过程中被看成是基岩分界线的期望位置。

算法测试平台为CPU酷睿i5-2300，内存16 GB的便携式计算机。模拟数据为一条含有86个节点的基岩初始分界线，以及23个钻孔轨迹线。在模拟基岩线空间分布不确定性时，根据钻孔轨迹线将图 1a所示的基岩分界线按长度投影到二维平面上，对构成基岩分界线的节点按照误差分布函数进行重采样，在二维空间上利用Monte Carlo方法进行模拟，再将模拟结果投影到三维空间上。经过反复对比分析，当模拟次数达到10 000次左右时，模拟结果的精确程度能够良好地反映地质剖面的真实状况。在此基础上，减少模拟次数会降低模拟结果的精确程度，而增加模拟次数不但不能显著地提高模拟结果的精确程度，反而会增加时间成本，降低效率。因此，本文对比分析了基岩分界线节点在服从不同误差分布函数参数时模拟10 000次的结果，如图 1所示，图中的紫色部分代表模拟得到的10 000条基岩分界线。不同误差分布函数的参数如表 1所示。剖面两个端点的误差分布函数可视作一维正态分布，其函数参数为参数4(表 1)；其余节点分别满足如下条件：

条件①，基岩分界线各节点均为正态分布(参数1)；

条件②，钻孔中基岩分界位置的空间误差服从正态分布(参数1)，其余基岩分界位置的空间误差服从圆域均匀分布(参数2)；

条件③，钻孔中基岩分界位置的空间误差服从正态分布(参数1)，其余基岩分界位置的空间误差服从矩形均匀分布(参数3)。

当基岩分界线节点空间误差服从不同的分布函数时，基岩线分布模拟结果特征不尽相同。由于钻孔上基岩分界位置的空间误差分布函数为正态分布，这类节点模拟结果呈现由中心向四周密度变低的现象。在非钻孔位置，当节点的空间位置误差服从均匀分布时，这类节点处的模拟结果中的点密度是均匀一致的^[10]。从基岩线的整体模拟结果来看，当分界线上各节点空间位置误差服从相同分布函数如正态分布时(图 1b)，基岩分界线分布相对平滑、均匀；当分界线上各节点空间位置误差服从不同分布函数时(图 1c、d))，模拟所得基岩分界线在节点位置的特征变化明显，节点间线段分布范围小于节点所在位置。然而，Monte Carlo方法模拟出来的结果虽然可以给出若干可能的地层分界线，但是不能清晰地表达不确定性的空间分布，也不利于不确定性的定量分析。

在3种条件下完成10 000次模拟的平均用时分别为171 ms、187 ms、156 ms。

通过Monte Carlo模拟方法可以得到剖面上若干可能地质界线的位置，然而同时显示这些可能分界线位置的建模数据空间不确定性表征方法难以进行定量分析和可视化(图 1)。为此，本文引入了“地质属性概率”^[11]的概念来分析和可视化建模数据的空间不确定性。

地质属性概率指的是空间上一点属于某一地质属性(如地层)的概率值，它与地层分界位置的空间不确定性具有密切联系。图 2示意性地说明了两个地层交界线附近空间点所属地层的概率计算方法。图 2区域由地层F和地层M组成，AB和CD分别为两次模拟得到的可能地层分界位置。若地层分界位置为AB，则点Q属于地层F；若地层分界位置为CD，则点Q属于地层M。假设边界线的模拟次数为N，点Q属于地层F的次数为n，则点Q属于地层F和M的概率分别为：

(1)

图 2 地质属性概率计算示意图 Figure 2 Schematic diagram of geological attribute probability calculation

图选项

将剖面离散化，根据式(1) 可以计算出各离散点所属地层的属性概率值，得到地质属性概率的分布，即剖面上各离散点属于基岩的可能性大小(图 3)。与用若干个可能的边界线来表示模拟结果的空间展布不同(图 1)，地质属性概率等值区的方式可定量可视化空间中每一个点属于某一地层的可能性，也可以直观地看到多种误差分布函数对地层界线不确定性的影响(图 3)。

a、b、c分别为以概率值方式显示的条件①，条件②、条件③ 的模拟结果。蓝色线为提取概率值位置。 图 3 Monte Carlo模拟结果地质属性概率三维分布 Figure 3 3D visualization of Monte Carlo simulation results with probability distribution

图选项

为了更细致地对比不同条件下地层分界位置的变化特征，将三维地质属性概率分布根据钻孔轨迹线投影到二维平面上，在y=-21 m处沿水平方向和x=4 002.85 m处沿垂直方向(三维空间对应位置见图 3)分别取不同条件下的平面地质属性概率剖面，可以得到不同条件下横向和纵向的地质属性概率值曲线(图 4)。通过概率分布曲线可以定量分析不同模拟次数、不同分布函数参数条件下概率分布的区别。从图 3和图 4来看，在组成分界线节点处的误差概率分布函数存在差异时，节点附近的属性概率值变化较大，但总体趋势尤其是在非节点上的概率值分布趋势大致一致，然而局部区域仍存在很大差异(图 4a)。例如：在图 4a中x=800 m附近，条件③ 计算得到y = -21 m处属于基岩的概率达到0.5，而其他两个条件计算得到的结果接近0.7；在x=1 500~1 700 m以及3 800~4 000 m附近，3个条件计算得到的结果趋势大致一致，但局部具体值还有细微差别。在分界线节点处为已知的误差概率分布，而在非节点处，其不确定性的纵向分布与两端节点的误差分布函数相关，不同分布函数对其有一定的影响，但是总体趋势一致(图 4b)。

a.y=-21 m；b.x=4 002.8 5 m。 图 4 不同误差概率分布函数组合模拟结果对比 Figure 4 Simulation results comparison of different error probability distribution functions

图选项

可见，当分界线节点误差服从不同分布时，分界线附近的地质属性概率分布是不同的，即分界线的不确定性分布是不同的。

假设同一地质要素不同位置上的误差服从不同分布函数，将建模数据的不确定性问题转化为地质属性概率分布问题，本文提出了一种基于Monte Carlo的三维剖面地质界线空间分布不确定性分析方法。Monte Carlo根据误差分布函数分别对不同位置的空间节点模拟采集过程，可以综合同一分界线上不同空间的误差特征，突破解析法难以集成多种误差分布函数的瓶颈，适用于分析多元地质数据条件下的空间不确定性。地质属性概率的引入，可以定量表征地质界线的空间不确定性分布，简化了Monte Carlo不确定性模拟结果的分析过程。

本研究中假定各节点误差分布相互独立，难以准确表征地质界线整体的不确定性。在使用本文所提出的方法分析地质数据的空间不确定性时，应当根据实际样本的数据量合理地选取随机模拟的次数，在满足模拟结果精确度的条件下尽可能地节省时间成本，以提高效率。此外，由于先验的误差分布函数决定了地质界线的不确定性分布，地质专家根据数据特征、实际地质背景以及自身经验等因素会直接影响到不确定性分析的效果。