0 前言

经过数十年的发展，三维地质建模已经成为矿产勘探^[1-2]，石油、天然气储层模拟^[3-4]以及土木工程^[5]等诸多领域的重要工具。本质上，三维地质建模从稀疏分布的地质数据中提取空间结构特征，并通过数学方法重建这些特征。人机交互建模方法能很好地重构地质结构间空间关系^[6-7]，但是这类方法费时较长且难以灵活更新模型。空间插值是一种有效获取地质数据空间结构模式重要且被广泛应用的手段^[8-9]。结合空间插值方法，利用隐式或显式函数来表征未知点和已知点之间的空间分布关系，例如距离函数、变差函数和体积函数^[10-13]，这类方法因一套参数只能生成一个确定的地质模型而被称为确定性方法。改变建模参数方能生成不同的空间模式和模型。

风险决策需要多个形态特征和空间分布相对逼真且具差异性的模型支撑，对快速构建和更新复杂地质模型提出了新的要求^[14-15]。通过扰动输入数据来生成多个地质模型^[16-17]的方法低估了数据稀疏区域的空间变异性。Røe等^[18]使用基于变差函数的随机模拟方法来形成多个地质构造界面，较好地反映了空间变异性。相比于基于变差函数的方法，即两点统计学(TPS)方法，多点统计学(MPS)方法能更好地捕捉复杂空间特征，并易于结合先验地质信息。因此，在过去二十年中，多点统计学得到了长足发展^[19-22]，并出现了多个面向三维地质结构模拟的多点统计学算法^[23-26]。

多点统计学三维模拟的一个重要核心在于如何获取三维训练图像。在一些简单地质现象如古河道的模拟中^[27-28]，通常使用基于过程或基于对象的方法建立三维训练图像。针对特殊情况需要其他方法来获取^[29-30]，例如在重建三维孔隙介质时，将利用X射线计算机断层扫描技术得到的三维图像作为训练图像^[31-32]。在矿产领域中，大多采用已经建立好的三维地质模型作为训练图像^{[23, 25, 33-34]}。在整个建模区域几乎没有可重复模式的情况下，Boucher等^[35]提出只模拟在用户定义的不确定区范围内的地质接触带，但是这样的不确定区需要额外的存储空间。可见，建立一个理想的初始三维地质模型是非常困难的。二维地质剖面图经常作为野外地质调查成果来表达地质结构。因此，获取训练图像的一个有效途径就是直接使用二维地质剖面图作为训练图像，而不是将一个其他方法生成的初始三维模型作为训练图像。Comunian等^[29]提出了一种二维序贯模拟算法，比较发现，其比概率集成和列表集成的方法效果要好，但是随着网格中已模拟部分的增加，模拟质量会下降。Gueting等^[36]提出了一种结合序贯二维模拟算法^[29]和直接采样法^[37]的算法，获得了比单独使用这两种算法更好的效果；由于该算法是对序贯二维模拟的结果再使用直接采样法，所以两种方法之间需要选择一个合适的切换点。

此外，在多点统计学模拟中如何集成从多个训练图像中提取的结构信息是另一个难题。大多数多点统计学算法采用单个训练图像进行模拟，不能有效提取和利用空间特征分布的不均一性和变异性。同时使用多个训练图像能够让多点统计学模拟结果更加合理^[24]，集成多个训练图像中的模式信息也就成为多点统计学模拟中的重要一环。Caers^[38]使用tau模型^[39]集成二维训练图像中的统计信息。Comunian等^[29]对比分析了用概率集成和合并列表的方法集成多个二维训练图像的信息。Silva等^[24]提出用条件概率和多点熵方法合并两个训练图像，但本质上这个方法仍然使用一个合并而成的训练图像进行模拟。因此，在模拟过程中仍需要关注如何整合多个训练图像的不同模式。

大多数序贯模拟多点统计学算法难以确定合适的模拟顺序。由于整个模拟网格只访问一次，所以模拟过程中产生的错误无法得到矫正；但是迭代方法可以有效矫正模拟错误，提高模拟质量^[40-41]。GOSIM(global optimization simulation)是一种基于全局优化的多点统计学迭代算法，有较好的模式复制能力^[41]。本文改进了GOSIM算法，提出了一种以二维钻孔地质剖面为训练图像的三维地质结构多点统计学模拟方法。

1 扩展GOSIM的MPS模拟方法

GOSIM算法^[41]首先从训练图像中随机抽取模式来初始化模拟网格，即获得初始模型，然后使用EM(expectation-maximization)迭代方法来缩小训练图像和模拟网格之间的差异。迭代分为E和M两个步骤，在E步骤(随机搜索)中，找到模拟网格中每个数据事件在训练图像中的近似最近邻；在M步骤(网格更新)中，对模拟实现中每个点的所有近似最近邻的对应点计算加权平均值，以此更新每个点的属性值。然而，GOSIM方法不能采用二维训练图像进行三维模拟。

本文提出的算法结合了基于模板的序贯模拟和基于GOISM的多尺度迭代算法。由于GOSIM不能直接使用二维剖面来构建三维地质模型，本算法先采用序贯模拟方法构建初始模型，然后将序贯模拟结果作为初始模型进行多尺度迭代。本文算法有3个主要步骤(图 1)。第一步是设置模拟参数，包括多尺度迭代的尺度数目、每个尺度的迭代次数、模板大小以及模拟网格的尺寸等。第二步是建立三维训练图像、三维模式库以及生成初始模型。在模拟网格降采样到最粗尺度时，扩展二维剖面以生成三维训练图像；然后，从三维训练图像中提取三维模式，根据重叠区分类之后建立三维模式库；从三维模式库中提取模式来进行序贯模拟，生成一个有粗糙地层结构的初始模型。第三步使用GOSIM算法进行地质结构模拟，通过多次迭代实现地质结构的三维重构。与GOSIM相同，本文算法每次迭代可以分为E步骤和M步骤。E步骤中，对于模拟网格R中的每一个模式P_R∈R，在训练图像I中搜索与之近似最相近的模式P_I∈I；在M步骤中，根据P_I对模拟网格R进行更新。整个过程结合了多尺度方法。最终输出最精细尺度的R作为最终结果。

1.1 三维训练图像构建

正如前人^{[25, 42]}指出的，合适的训练图像是多点统计学模拟获得高质量结果的关键因素之一。本文提出一种用二维地质剖面构建三维训练图像的方法。将4个闭合的二维剖面作为模拟区域的边界(图 2a)。直接从二维训练图像提取的二维模式不能较好地表征三维空间特征，不利于重构三维地质结构。为了获得地质数据的三维空间模式，需要将二维剖面图的维度提升到三维。本文将二维剖面朝模拟网格的中心方向扩展，扩展的范围约等于(大于或等于)一个模板的宽度(图 2b)，使得后续模拟能够从三维训练图像中提取完整的三维模式。

地质剖面从二维扩展到三维的过程是逐层进行的，图 2c、d是图 2a、b中训练图像局部扩展的放大示意图。当对内层网格进行扩展时：如果当前扩展点是顶点，例如图 2c中的P，点P的属性值就会从外层中两个与P相邻的点的属性值中随机挑选, 如黑色箭头指示，图中P选中了黄色；如果当前扩展点不是顶点，如图 2d中点Q，其属性从外层3个相邻点的属性值中随机挑选。如此依次进行，在填充完一层后再填充下一层，按照设定的层数，从外层训练图像向内逐层扩展，最终得到的三维训练图像是具有一个模板厚度的地质剖面，如图 2b图所示。

1.2 模拟网格初始化

在建立地层模型时，没有合适初始模型的迭代随机搜索方法会使模拟结果不能很好地再现地层分布的方向延展性特征。GOSIM在进行EM迭代之前通过随机函数生成一个初始模型。对于二值图像如河道分布图的模拟，通过EM迭代能获得不错的效果，但是对于具有多种属性的地层模型，随机初始化无法重建有空间延展性的地层。因此，在应用EM迭代方法之前，需要获得合理的初始模型。模拟网格是通过序贯模拟实现初始化的。

本方法中使用了多尺度的方法，模拟从最粗糙尺度开始，逐步进行到最精细尺度，每一个模拟尺度的边长是前一个尺度的两倍。先将带有训练图像的模拟网格降采样到最粗尺度，用1.1节中的方法建立三维训练图像。初始化是一个从已知区域向未知区域推进的序贯模拟过程，即从三维模式库中搜索并提取合适的模式，用这些模式依照路径顺序填充模拟网格中的属性空白区域(图 3)。在每次搜索中，根据模拟网格中当前数据事件与已知区域的重叠区，在模式库中找到与重叠区对应区域相似度符合要求的模式，从模式库中取出然后粘贴到当前数据事件的位置。这种路径有利于网格内部的属性被外围训练图像约束，保证地层层序、分布和形态的一致性。

当前数据事件被更新之后，已模拟区域得到了扩充。在模拟下一个位置时，根据数据事件与新已知区域的重叠区搜索模式。在网格中不同位置模拟时，当前数据事件与已模拟区域的重叠区形状是不同的(图 4)。例如，在模拟的起点位置，即模拟网格内部空白区域边缘的一个角上，三维数据事件只与三维训练图像构成重叠区，并且这个重叠区为两个具有一定厚度的相交面，重叠区厚度为预先设定的参数。在本算法中，最初的已知区域是模拟网格边缘的三维训练图像，于是模拟就从网格的边缘向内部逐层进行，从模拟网格中空白区域与三维训练图像相交的位置开始。

为了加快序贯模拟搜索模式的速度，在生成的三维训练图像中提取出所有的三维模式，分类之后建立三维模式库。分类时，首先提取重叠区形状类型，然后利用重叠区相似度对每类重叠区的模式进行分类。相似度计算使用Hamming距离^[43]，即统计所有模式同一重叠区内属性值相同点所占的比例。在搜索时，数据事件先与每一类模式的代表模式比较，如果相似度符合设定阈值，则在这一类中随机挑选一个模式作为结果；如果遍历所有分类都没有找到相似度满足阈值的模式，则将阈值降低，再重新与每一类的代表模式比较。考虑到模拟网格和模式的尺寸，按照预设的重叠区宽度决定的模拟间隔，序贯模拟并不一定能刚好填满整个网格，所以在1.1节建立三维训练图像的过程中，可以计算模拟网格边长和模式大小的关系，适当增加训练图像扩展时的层数，确保随后的序贯模拟能够填满模拟网格的所有未知区域。

1.3 随机搜索和网格更新

通过简单模式拼贴式的序贯模拟得到的模型不是合理的最终模型，因为初始化模型存在大量拼贴痕迹(图 3c)。这是由于依据重叠区相似度在模式库中搜索模式时，难以找到与数据事件相似度满足阈值的模式，所以每次循环结束时都要降低阈值重新搜索，使得拼贴的两个相邻模式之间的过渡并不自然。GOSIM算法中的核心步骤是EM迭代，同时结合多尺度方法，虽然不能有效重建有方向延展性的地层特征，但是可以有效平滑模型中的拼贴痕迹，优化初始模型，让模型特征更接近训练图像特征。

以最粗糙尺度模拟网格中完成的初始化结果作为迭代的初始模型，每一个尺度迭代完成后，结果作为下一精细尺度的初始模型。迭代步骤如图 5所示。在每一个尺度的模拟中，都会执行若干次EM迭代，每一次迭代都包括E、M两个步骤。E步骤中，先为模拟网格中的每个像素赋予一个从训练图像中随机提取的模式，并且计算这些模式与初始模型中数据事件的距离。然后搜索分为两步进行：第一步是对模拟网格中所有像素进行遍历，在访问每一个像素时，将赋予该像素的模式和赋予临近像素的模式进行对比，在这些候选模式中找到一个与当前模拟网格位置数据事件最相近的模式，然后赋予当前像素，替换原有的模式；第二步是随机窗口搜索，通过在三维训练图像中设置范围由大到小的窗口，在每一级窗口随机提取模式与当前模式比较，如果更相似再替换当前模式。E步骤完成之后，每一个像素点都有一个和当前位置数据事件近似最相似的模式，当前位置数据事件是以当前访问位置像素点为中心的数据事件。M步骤根据E步骤搜索的结果更新模拟网格中的属性值，由于每个像素都包含一个模式，所以将每个像素包含的模式粘贴到以该像素位置为中心的数据事件处时，每个像素上会有多个相邻模式重叠的属性值，对这些重叠的属性值进行加权后，粘贴到该像素的位置。一次EM迭代中，通常E步骤执行多次，因为需要找到近似最相似模式，而M步骤只需执行一次。

2 模拟实例

本研究建模实例位于广东省广州市某地，区内地形平坦，地貌类型单一。地表为第四系所覆盖，钻孔揭露主要地层有下石炭统石磴子组(C₁s)，上白垩统大塱山组第一段(K₂dl¹)、上白垩统大塱山组二段(K₂dl²)，全新统桂洲组(Qhg)和人工填土(Q^ml)。石磴子组主要岩性为灰岩，大塱山组岩性主要为灰质砾岩。石磴子组与大塱山组呈角度不整合接触，由于溶洞发育普遍，接触边界犬牙交错。桂洲组为海陆交互相冲淤积层，不整合于大塱山组和石磴子组之上。可见，研究区晚石炭世后地壳抬升，地层接受风化剥蚀；晚白垩世时整个地区下沉，接受沉积；此后，地层再次抬升、剥蚀，直至全新世下沉接受沉积。

本研究依据钻孔绘制的4条闭合地质剖面(图 6)作为建模用训练图像，模拟区域划分为566×220×215的网格。模拟采用了3个尺度，每个尺度边长是后一个尺度的1/2，分别是141×55×53，283×110×107和566×220×215。

2.1 无条件模拟

最精细尺度的无条件模拟结果(图 7)表明：地层接触界面基本光滑，层序正常，相互间比例合理，无地层缺失，地层连续性及起伏形态与二维训练图像相似，地层在模拟结果中的分布位置受到它们在训练图像中分布位置的约束。例如人工填土层、全新统桂洲组和下石炭统石磴子组在二维训练图像中水平方向上较为连续，在模拟结果中，这些地层在三维模拟区域内的水平方向上分布相对连续。当某一个地层在训练图像中分布较少时，则在模拟结果中该地层一般只会出现在边界上靠近包含该层的训练图像的位置，例如上白垩统大塱山组第一段出现在训练图像的西北角和南剖面的两个局部位置，模拟结果中该地层基本被限制在了模拟网格的西北角或模拟网格内靠近南剖面那两个位置的地方，很少会出现在模拟网格的其他地方。

图 8是参数相同情况下4次无条件模拟的结果，将人工填土层与桂洲组作透明处理。可以看出在参数不变的情况下，多次模拟会产生不同结果，这些结果中地层的形态和分布位置不尽相同，但是整体的形态特征和训练图像保持一致，并且与边缘二维训练图像能够自然衔接。

2.2 条件模拟

在模拟区域内选择了3个钻孔J702，J703和ZK504-3作为约束数据。以此3个钻孔作为约束条件的模拟结果(图 9)表明，条件数据对于模拟网格中该位置的地层属性有较强的约束作用，钻孔周围地层属性与钻孔属性一致。同时也能发现单一条件数据的影响范围有限，如果建模区域内有更多钻孔数据，整体地层的构建就会更接近真实情况。

为了比较模拟结果与训练图像之间的相似程度，抽取了两个无条件模拟和两个条件模拟结果，分别统计了训练图像和模拟结果中各地层占模拟单元数的比例(表 1)。三维训练图像是从二维地质剖面扩展而来，为保证训练图像的有效性，二者在地层占比上应具有相同的趋势。从统计结果来看，尽管三维训练图像和二维地质剖面中的各地层之间比率有所差异，但地层分布的总体趋势相似，石磴子组占最大比例，人工填土占比最小。

模拟结果中各地层比例基本相似，地层属性分布具有较高的相似性，但又有所区别，因此模拟结果合理且具有一定的空间不确定性。但是，需要注意的是，大塱山组第一段在模拟结果中有大幅减小的趋势，而石磴子组占比均比训练图像中的高。从二维地质剖面上地层分布的情况来看，石磴子组在建模区域大规模分布。因此，本文提出的算法对大规模分布地层有可能存在过度拟合的情形。

3 结论

1) 针对三维地质建模中三维训练图像获取困难的问题，本文提出了一种用二维剖面作为训练图像的多点统计学三维模拟方法。本方法结合了序贯模拟和迭代的方法，将二维剖面扩展为三维训练图像，并从中提取三维模式建立模式库，利用序贯模拟方法生成初始模型，采用GOSIM的迭代方法对初始模型进行迭代优化，使模拟结果更为合理。

2) 建模实例结果表明，本论文所提出的方法能够准确模拟研究区的地层层序，很好地再现二维地质剖面所反映的地层结构关系。

3) 本文所提出的方法确保了训练图像对内部模拟网格的约束，保证了模型中地层层序的完整有序，并且消除了网格中有拼贴痕迹、不连续的错误，但是对于大规模分布的地层可能存在过度拟合的情形。