基于多点统计学的地质体三维重构及其在地铁工程中的应用*

刘成军，刘恒光，侯卫生，3，4，王典，陈晓丹，张华，陈勇华

1. 广州地铁设计研究院股份有限公司，广东 广州 510010

2. 中山大学地球科学与工程学院，广东 珠海 519082

3. 南方海洋科学与工程广东省实验室（珠海），广东 珠海 519082

4. 广东省地球动力学作用与地质灾害重点实验室，广东 珠海 519082

经济发展促进了城市尤其是大城市、城市群的快速扩张。大力发展地下轨道交通是解决城市高速发展中所面临的交通拥挤、实现土地集约利用等问题的一个良好途径。然而，漫长地质过程不仅形成了形态复杂多变的地质体，也改变了岩石的结构和力学性质。建立合理的、高精度、大比例尺的三维地质结构模型，明确地下轨道交通沿线的地质条件，是地铁工程的设计和施工中一个重要的数据框架和基础。

三维地质模型是基于现代空间信息理论和技术建立的一种具有地质意义的数学表达，给出了易于理解的地质体空间结构三维透视图，为基础地质调查、工程设计、施工和决策等提供基础数据和决策支持［1-6］。假设数据分布服从二阶平稳条件，基于变差图的空间分布函数，两点地质统计学方法可充分利用软数据和硬数据约束，实现复杂地质空间结构的三维重构［7-9］。然而，基于变差图的空间关系函数缺乏对数据整体空间相关性的度量，使得两点地质统计学方法在模拟非线性地质结构上存在显著偏差，难以准确揭示地质结构的空间变化［7，10-13］。

与两点统计学相对应的是，多点统计学（MPS，multiple-point statistics）方法从训练图像（TI，training image）中提取多点间不同尺度模式结构，通过随机模拟过程重组这些空间特征。MPS 着重关注空间中多点之间相关性的表达，可在模拟过程中引入地质知识，使模拟结果不仅忠于建模数据，还能有效描述地质体的复杂空间几何形态，更好地服务于各种科学研究和应用领域，如油气藏勘探［14-17］、煤田地质［18］、水文地质［19-21］、遥感制图［22，23］、矿产资源评价［24］、微观结构特征［25-28］等。MPS 模拟方法有两个关键步骤：一是当前数据事件搜索，即从TI 中搜索出符合当前数据事件条件的空间模式；二是数据事件更新，即对所搜索出的模式进行更新和粘贴。从粘贴像素的数目来看，已有MPS 方法大致可以分成基于点模式和基于模板的方法两大类。前者在每次模拟过程中仅粘贴一个像素点，常用方法主要有SNESIM［29］、DS［7］、HOSIM［30-31］和IMPALA［32］等。基于模板的方法在每次模拟中粘贴一个模板大小的数据，如GOSIM［33］、SIMPAT［34］、FILTERSIM［35］、DISPAT［36］和MS-CCSIM［37］等。除GOSIM 外，其余方法均以序贯模拟为基础。这种基于序贯模拟的方法在数据事件更新过程中所产生的误差难以得到修正，并在模拟过程中不断累积以致模拟结果偏差较大。引入具有更低敏感度的栅格路径方法可以减小误差，但误差累积仍蕴含其内。结合图像重构中的多尺度迭代计算方法，GOSIM 克服了传统迭代MPS 方法难以收敛的问题，并在模拟过程中不断修正误差，获得了更符合实际的模拟结果［33］。

本质上，三维地质建模的核心在于寻找并明确已知地质数据之间的空间关系，推测未采样位置的地质属性值。与区域地质调查不同的是，受城市建设的影响，地铁工程地质调查的横向宽度相对较窄，钻孔数据却相对密集，形成了具有数据量少、但局部密度高的数据特点。常规的三维地质建模方法难以从已知数据中提取出地质属性的空间分布结构。因此，数据量少、局部密度高数据条件下的地质结构三维重建成为三维地质模型在地铁工程中深入应用的一个瓶颈。本文以Extended GOSIM 算法［38］为核心，结合地铁工程数据特点，以二维地质剖面作为建模的基础数据集，从分级融合的角度探讨面向地铁工程应用的地质结构三维重建方法。

1 基于MPS 的地质结构三维重建方法

1.1 MPS的基本概念

MPS 使用数据模板扫描TI 获取地质数据的空间模式以及地质变量的统计量，并利用这些统计量随机模拟出与TI 相似的特征。TI 是对现实世界的一个随机过程的抽象描述，是描述地下介质结构的二维或三维图像。可见，TI 的来源也可以多种多样，如露头照片的解释，地质剖面图，或者是基于对象或基于过程的模拟实现。本研究采用了二维钻孔地质剖面作为TI进行三维地质模拟。

数据模板T是由空间中n个向量组成的一种几何结构。假定模板中心位置为u，那么模板就可以表示为：T=｛u，u1，u2，…，un-1｝=｛u+h0，u+h1，u+h2，…，u+hm｝（如图1所示），其中h0，h1，h2，…，hm为滞后距离。

图1 一个二维的5×5的数据模板示意图Fig.1 Schematic diagram of two dimensional data template with size of 5×5

MPS 在获取TI 的空间模式时，从TI 的某个位置开始，以数据模板大小为窗口，获取窗口内的数据属性分布状态；然后，滑动窗口，再获取新窗口范围内的属性分布，直至TI内所有的结构均有效获取，其过程如图2（a）所示。扫描结束后，MPS即可获得如图2（b）所示的空间模式集。通常而言，对TI 进行搜索的时间相对较长。为了加速搜索模式，本研究在获取空间模式集后，构建了模式数据库，避免了模拟过程中反复对TI进行全局搜索。

图2 MPS的空间模式获取过程示意图Fig.2 Schematic diagram of extracting spatial patterns from TI in the MPS

考虑到一个MPS 的实现（模拟结果）R与训练图像集合UTI在空间结构上具有一定的相似性，那么R中就应该有尽可能多的模式来自UTI。令PRi表示R中的一个模式，PTI为UTI中的模式，R和UTI之间的差异就可以定义为［39］

其中D（·）为广义距离函数。对于任一个PiR∈R，找到与它最接近（相似）的PTI∈UTI，R与UTI相对应的模式距离之和就代表了R和UTI之间的整体不相似程度。可见，d（R，UTI）是MPS 模拟中需要最小化的目标函数。式（1）给出了基于MPS的地质结构三维重建过程中的全局核函数，相应的目标函数就是实现式（1）的最小化，即

其中Rfinal为最终的实现结果。令P1TI，...，PnTI表示UTI中的所有模式。对于PR∈R和PTI∈UTI，P_TI是PR的“最近邻”，如果它们满足

式（3）表明P_TI是与PR最相似的模式。如果P_TI与PR的距离大于PR与它的最近邻的距离值，则称P_TI是PR的“近似最近邻”。这样，我们就可以把问题从实现的整体相似性转换到空间模式相似性的求解上。

1.2 基于MPS的地质结构三维重建方法

Rfinal值的求取，需要同时知道R中的模式以及UTI中与它们最相似的模式。然而，式（2）中的R和Rfinal都是未知的，简单方法难以实现该问题的求解。前人研究表明，期望最大化（EM，expectationmaximization） 算法非常适合解决这类优化问题［40-41］。EM 算法采用迭代计算的思路来实现最小化目标函数过程，包含了两个主要步骤：E步骤和M 步骤。E 步骤为每个PR∈R找到最相似的PTI∈UTI，即模拟网格中每个数据事件在UTI中的近似最相似模式；M 步骤则根据找到的包含该点的近似最近邻的对应点计算加权平均值，更新该点的属性值。EM 算法实现基本过程如图3 所示。由此，EM 迭代可以实现最小化目标函数。

图3 EM迭代算法实现流程Fig.3 The implementation process of the iterative EM algorithm

地铁工程大多采用钻孔勘察方式来直接获取地质体的空间分布。这些钻孔一般沿地铁线路呈条带状分布，在线路的横截面上分布较少。仅通过钻孔数据难以获得足够的地质结构的空间延展信息。依托已有地质调查成果，根据钻孔岩芯数据所绘制的二维地质剖面不仅再现了单个钻孔内的地质结构信息，还融合了地质专家对地质对象空间展布特征的认知。因此，本文以二维钻孔地质剖面作为建模数据集进行地质结构三维重建。

然而，二维剖面并不能直接提供地质对象的三维空间展布模式。Extended GOSIM算法可以有效实现以二维钻孔地质剖面为UTI的三维地质结构重构［38］。考虑到破碎带中存在不同类型的岩性，以及地层可能存在的沉积旋回，本研究以Extended GOSIM算法为基础，提出了一个融合分级建模思想的多尺度MPS地质结构三维重建算法（见图4）。即，首先将具有分级特征的地质对象作为一个整体进行特征结构的重建，获得该地质对象分布大致的范围，然后，以其分布范围为约束，再进行次级结构的细化模拟。

图4 融合分级建模思想的多尺度MPS地质结构三维重建算法实现过程Fig.4 The multiple-scale MPS implementation process of 3D geological structure reconstruction with hierarchical modeling idea

在本研究中，有效获取地质对象的三维模式是实现地质结构合理重建的基础和前提。利用Extended GOSIM 算法所提出的二维剖面扩展步骤［38］，我们以随机的模式将二维剖面扩展到一个模板大小厚度的三维剖面，然后再从扩展剖面中提取出地质结构的三维模式，并建立三维空间模式数据库。地质剖面从二维扩展到三维的过程是逐层进行的。首先将二维剖面按照实际位置放置于三维模拟网格中。此时，除剖面所占据的网格外，其余网格节点的属性均未知。当对二维剖面进行扩展时，待扩展点的属性值就从外层中与当前待扩展点相邻的已赋值或剖面点的属性值中随机挑选一个。如此依次进行，在填充完一层后再填充下一层，按照设定的层数，从外层训练图像向内逐层扩展，最终得到的三维训练图像是具有一个模板大小厚度的地质剖面，如图5 所示。可见，在扩展的三维地质剖面时，并不是简单地将二维剖面属性复制到模拟网格节点中，而是以二维剖面为基础的随机赋属性值过程。图6给出了剖面⑤在西、中和东侧3个位置上的切片，其中黑色椭圆的虚线范围标识了3个切片的属性值差异位置。在剖面的不同位置上，属性分布差异主要出现在岩性变化的交界处。

图5 扩展后的三维地质剖面图Fig.5 Three-dimensional(3D)geological cross-sections expended from 2D cross-sections

图6 剖面⑤经过拓展后在不同位置的切片图（虚线椭圆标注了不同切片的差异）Fig.6 Slices in different positions that cut from different position of the extended geological cross-sections(Elliptic curves in dash line show differences between slices)

为了加快从模式库中搜索模式，在完成二维剖面扩展后，利用移动立方体扫描三维地质剖面，建立三维模式数据库。所采用的移动立方体大小为建模过程中所采用的模板一致，其过程与1.1中所阐述的基本一致。所不同的是，移动立方体所得到的是三维空间结构。在建立模式数据库时，先利用汉明距离［42］来计算三维模式之间的相似度，通过聚类的方法，将具有较小汉明距离的模式划分至相同的子模式类中。当所有模式都被聚类并建立索引后，就可获得建模所需的三维模式数据库。

2 地铁站点地质数据分析

本研究以广州市轨道交通十一号线江泰路站址为研究实例（见图7）。研究区位于晓港-赤岗碎屑岩台地与冲-海积平原地貌交界地带［43］。钻孔揭示的地层主要为白垩系和第四系。第四系主要有杂填土、第四纪残积土（Q4el）、全新统（Q4）和上更新统（Q3）河流三角洲相和海陆交互相淤泥质土、砂层、粉质黏土，其中全新统一般几米到十几米，第四系总厚度可达20 余米。研究区揭露的白垩系地层主要为上白垩统三水组康乐段（K2s1）和下白垩统白鹤洞组广钢段（K1b2）。康乐段（K2s1）属内陆湖泊相为主的粗砂-细砂碎屑岩建造，岩性主要为砂岩、含砾粗砂岩、砾岩，以及钙质粉砂岩和钙质泥岩。在颗粒组成上表现为下粗上细，泥质胶结为主，粉细粒结构，中厚层状构造。受岩体胶结物及颗粒组成不同影响，风化岩中常出现不同风化程度的夹层，导致部分地段不同风化程度的风化岩交替出现。广钢段（K1b2）为内陆湖泊相为主的粗砂-细砂碎屑岩夹不纯碳酸盐岩，岩性主要为粉砂质泥岩、泥质粉砂岩、粉砂岩、细砂岩、含砾砂岩、粉砂质泥灰岩和泥灰岩互层。中厚层状构造，泥质、钙质胶结，并含微-薄层及星点状-网状石膏。该段顶部被三水组康乐段平行不整合覆盖。

图7 研究区基岩地质图（改自文献［44］）Fig.7 Bedrock geological map of the study area(modified from ref. [44])

钻孔岩芯还揭露出英安斑岩、安山岩等火成岩。火成岩主要从上白垩统三水组康乐段与下白垩统白鹤洞组广钢段之间侵入，呈柳叶刀形，后受广三断裂切割，呈两段分离的扇叶型。其中，英安斑岩（ηγ53-3）岩石呈斑状结构，基质为霏细结构，斑晶25%～40%，其中石英为15%，其余为钾长石、斜长石、角闪石、黑云母等。安山岩青灰色夹暗紫红色，岩体破碎，岩芯呈碎块状夹块状，斑状结构，块状构造，见气泡，岩质较坚硬。

前人研究表明，广三断裂在本站与线路相交，走向基本上为近东西向或略有偏转，倾向南或略有偏转，倾角50º～85º［45］。在该站址的多处钻孔岩芯中发现有断层泥、断层角砾和碎裂岩等断层现象。构造岩以碎裂岩为主，该断裂南东端明显错开白垩系地层，呈逆时针扭动方向，属斜冲正断裂。已有研究表明，广三断裂在第四纪早期曾有过多次活动，但晚更新世以来未见地表或近地表的活动迹象。

本研究在江泰路站址13 288 m2范围内收集了38 个钻孔的岩芯数据，并对钻孔数据进行分层解释，结合浅层地震和井间电阻率CT 解释剖面，绘制了6 条交叉剖面（见图8）。其中，剖面⑤为钻孔CT 测量位置，剖面⑥为联络剖面以约束断裂形态和走向。

图8 研究区地质数据分布示意图Fig.8 Geological data collected in the study area

3 模拟结果分析

以图8中的钻孔地质剖面作为建模数据，本研究利用前述方法重建了建模区域内的三维地质结构。模拟过程采用了2个尺度，每个尺度的网格大小是后一个尺度的8倍。在最精细尺度下，模拟网格数为288×352×604，各方向上网格精度为0.25 m。

除火成岩和破碎带外，其余地质体呈近水平分布，地层之间没有出现倒转和穿插现象（图9）。其中，杂填土、粉土和粉质黏土呈层状水平连续分布于研究区域中，其平均厚度分别为4.09、3.94和3.23 m（见图9：d～f）。中、粗砂以及淤泥质土少量、零星分布在模型中（见图9f）。破碎带主要分布在建模区的东北部，主体呈东西方向延展，在中部向下尖灭于模型底部。破碎带顶部起伏较大，高程变化范围为-27.86～-0.58 m。破碎带中依岩性特征，可以进一步细分为断层角砾、断层泥和碎裂岩。其中，断层角砾是破碎带主要岩性，占总体体积的81.23%；断层泥和碎裂岩并没有大面积、连续性分布，仅零星分布于破碎带中（见图10：b～d）。可见，以已建立的破碎带整体三维结构为约束，再利用随机模拟方法可以很好地重构破碎带次级岩性分布。

图9 研究区三维地质模型（单位：m）Fig.9 Three-dimensional（3D)geological model of the study area(unit:m)

为验证模型的有效性，本研究给出了地质剖面和模拟结果中各地质体的占比情况（见表1）。除了破碎带与或火成岩外，各地质体在模型中整体占比与在剖面中的占比基本接近。破碎带所包括的断层角砾、断层泥和断层岩整体分布在建模区的北部地区，分布并不均匀，模型占比均小于其在剖面整体中的占比。此外，在剖面中只出现在破碎带周围的火成岩，在建模结果中也表现出相似特性，其占比也远小于剖面中整体中占比。

表1 研究区三维模型与地质剖面中各地质对象占比Table 1 The proportions of different geological bodies in cross-sections and the final 3D model %

从钻孔所揭露的岩性看，工程力学性质相对较差的主要为破碎带和风化火成岩。站址与破碎带模型的空间位置关系（图10a）表明：站址范围内涵盖了部分破碎带，破碎带体积约占站址范围的3.2%。站址的中东部底板与破碎带相交。为更好分析站址底部岩性分布情况，图11 给出了站址底部以下20 m 范围内的岩性分布情况。在100 m 以东区域，站址底板分布有大量断层角砾、断层泥和火成岩。在70 m 以东区域，站址下方具有厚度不等的破碎带分布，以西区域以砂岩为主。可见，在进行站址底板设计过程中，应结合站址底部岩性变化特征，分段设计底板构筑物结构。

图10 破碎带模型及地铁站址关系（单位：m）Fig.10 The spatial relationship of the metro station and fault zone model(unit:m)

图11 站址底部以下20 m内岩性分布情况（单位：m）Fig.11 Lithology distribution within 20 m beneath the bottom of the metro station(unit:m)

4 结 论

针对地铁工程应用中数据集中、但难以获取地质对象空间分布模式的难题，本论文以二维钻孔地质剖面为建模数据集，依托Extended GOSIM算法，提出了一种融合分级建模思想的多尺度MPS 地质结构三维重建方法。广州某地铁站址三维地质模型的构建实例表明：分级建模的思想可以很好地实现断裂中不同类型岩性的空间分布，也可以将其推广到沉积旋回地层的三维重建中；所提出的三维建模方法能充分利用有限数据内所蕴含的地质对象间空间关系，有效识别、提取和重建地质对象的三维空间模式。可见，基于MPS的隐函数建模方法能构建出高精度的三维地质结构，可以为地铁工程设计、施工提供精细数据基础。以随机方式将二维地质剖面扩展成三维剖面在一定程度上可以解决三维地质结构难以获取的问题，由于扩展范围只有一个模板大小，简单随机外推方法所获得的三维地质剖面与二维地质剖面间的差异过小，也难以获得地质对象的全局结构特征。