序贯指示模拟法在三维地质模型中的应用——以西辽河平原为例

杨明武 龙文华

1.中煤科工集团北京华宇工程有限公司 北京 100120

2.广东省地质环境监测总站 广东 广州 510510

地质体三维建模这一概念由加拿大的学者Simon W.Houdling在1993年提出后，近年来随着芯片技术的飞速发展，计算机软件运算能力大幅提高，因而，计算机图形学理论的日趋渐完善，地学领域的地质专家对三维地质结构模型的研究和计算机领域的研发专家对三维可视化的研究也取得了极大的突破。目前,国内外专家开发出了诸多的地学方面商业化应用软件,如:灵图的VRMap系统（中国）、加拿大金康公司开发的GemCom软件、美国考格尼塞斯公司的GeoSec3D软件、T-surf公司的GOCAD软件（法国）、以色列帕拉代姆地球物理公司的GeoDepth系列软件、美国CTech公司开发的环境地质可视化系统(EVS&MVS)及地下水模拟系统GMS(美国)等。从这些软件的使用功能来看，它们大多都能利用区域地震资料、场区钻井资料和地表地形资料生成三维地质模型解释成果,且各有各的特色。但是以上软件多针对油气藏、物探及采矿行业，仅EVS&MVS和GMS在水文地质、环境地质三维建模方面较为成熟，但应用上多集中于局部范围内的研究，对于区域松散层的建模研究较少。

GOCAD软件作为新一代地质建模软件的代表,具有通用地质模型所设想的功能——全三维、全拓扑、使用空间插值算法、拥有完善的地质统计分析功能等特征，达到了半智能化地质建模的世界最高水平,在国外受到众多用户青睐[1]。

由于在大区域且存在复杂相变的松散沉积盆地的三维模型构建时，用非常有限的离散点数据进行数学插值的方法描述空间岩性分布格局时会产生较大的模型失真。本文试图应用GOCAD中的序贯指示模拟功能对西辽河平原地层结构建立三维地质模型，探索一种能够真实反应在大区域盆地内构建三维地质模型有效的方法，同时对西辽河平原的地质、水文地质结构提供可视化的三维认识，为各含水层组的分布、变化特征的深入研究提供了依据，为水文地质概念模型的建立奠定了基础。

1 序贯指示模拟原理

序贯模拟是地质统计学的一个重要组成部分，是将序贯思路与克立格插值相结合的条件模拟方法，将沿着随机路径序贯求出各地质网格结点的条件累积分布函数值作为主要开发思路，并从中取得数值模拟值。常用的有序贯高斯模拟和序贯指示模拟。序贯高斯模拟主要适用于满足高斯分布的连续型数据场，而序贯指示模拟则主要适用于离散型数据场，也可用于离散化的连续型数据场，其基本原理为[2,3]：

利用线性无偏估计克里格法,来估计条件概率的公式为：

基于上，由于K阈值已经给定，那么对于任何一个待估位置，各个阈值都有一个方程组和其对应，在变量Z(x)的变化范围内，用K阈值对该估计范围离散化，通过求解每个阈值对应的克里格方程组，可以得到每个阈值下对应的累积分布函数F{Zk,x|(n)}，然后，采用Monto Carlo法来求得随机函数Z(x)在该位置的一个具体实现值。当随机变量Z(x)为分类变量(如岩相、沉积微相等)时,只需将指示指标变量的定义变为:

式中，tk为分类变量类型k，对应于Z(x)的k种分类变量类型由阈值的个数K表示,而估计的条件概率公式变为:

2 序贯指示数值模拟的步骤[4]

首先将模拟变量划分为K个阈值(设有n个初始值,L个网格节点),同时，定义一个通过所有网格节点的随机路径，然后按如下步骤进行序贯指示模拟的实现。

(1)定义查找范围内的条件数据，条件数据包括邻域内的原始指示数据和模拟得到的指示值；

(2)对任意k=1,2,…,首先要由K所对应的阈值指示协方差数值模型建立普通克里格模拟方程组,然后求解,最后由指示条件数据的线性组合计算该阈值的分布范围函数；

(3)由Monto Carlo法得到每个节点上的模拟值；

(4)按上述已定义的K个阈值,将模拟值转换为K个指示估算值；

(5)将序贯模拟结果归入条件指示数据集中；

(6)重复上述(1)-(5)项，沿着随机路径对下一个节点进行序贯模拟，直到所有的结点都被模拟。

3 应用实例

3.1 研究区概况

研究区位于内蒙古自治区的中东部，包括吉林省和辽宁省的部分地区，属西辽河平原。区内的总体地形走势为南北高、中间低；最高点高程为670.0m，最低点高程119.0m。

区内地层主要由第四系粘性土、砂和砾石层等构成；其地层厚度变化规律为：中心厚，边缘薄，做为沉积中心的开鲁一带，地层厚度可达200m。地层岩性总体变化规律为：水平方向，自山前向平原，从上游向下游颗粒由粗变细，变化过程为砂砾石—中粗砂—细砂、粉砂，且粘土夹层增多，厚度增大。

上部地层由全新统、上更新统顾乡屯组、排头营子组的中细砂、粉砂、细砂和泥质细砂组成；中部地层由中更新统大青沟组冲洪积、冲湖积的粗砂、中砂、细砂和淤泥质粉砂组成；下部地层则由中更新统白土山组冰水、冰碛的砾石、砂砾石、泥质砂砾石和粗砂、中细砂等组成。

3.2 三维地质模型的建立

3.2.1 岩性概化

本次工作收集到研究区共641个钻孔，通过分析发现，钻孔年代不一，岩性描述复杂且缺乏统一性，考虑到序贯数值模拟要与水文地质模型相结合的原则，根据第四系地层岩性的颗粒粒径、透水性等水文地质特性，将模拟区内的第四系地层岩性概化为粘土、粉土、细砂、粗砂和砾石四类。

3.2.2 建立三维地质体模型

根据地面高程的DEM数据，确定三维地质体模型的顶板，同时利用钻孔资料及地质剖面图确定其底板，由此建立研究区的三维地质体模型（见图1），该模型为进行岩性随机模拟确定出空间范围，这样可以从很大程度上避免顶底板交叉的问题。

图1 研究区三维地质体模型

3.2.3 序贯指示模拟

3.2.3.1 频率及阈值计算

将第四系岩性作为类型变量进行指示变异函数计算。根据岩性概化，分别用模型变量1、2、3、4来代表粘土类、粉土类、细砂类、粗砂和砾石类。在模拟进行指示变异函数计算时，分别用模拟类型变量1、2、3、4作为阈值，见表1。

表1 各类地层岩性频率及其阈值

3.2.3.2 变异函数参数计算

选择平面上的45°方向，135°方向，垂直方向的滞后距0.01，水平方向滞后距0.025、容差角20°、带宽0.2，分别计算四个类型变量的指示变异函数参数，计算结果见表2。

表2 各类地层岩性变异函数参数计算结果表

3.2.3.3 序贯指示模拟

将计算所得的数值模拟参数输入到GOCAD的序贯指示模拟模块中进行序贯指示模拟，即可得到西辽河平原区松散层三维地质模型。图2可视化呈现了模拟区第四系地层岩性分布，图3以栅栏图的方式呈现了研究区第四系三维地质结构。从图2、图3可见研究区内地层岩性颗粒粗细分布与变化，区内地层岩性以细砂类、粗砂类为主；向东到中部岩性粒径由粗变细，粘粒物质的含量逐渐增加，砂层与粘土层交错分布。

图2 研究区三维地质模型

图3 研究区立体栅栏图

3.2.4 模型验证

为了验证序贯指示模拟地质模型的效果与精度，本文选择了一条实际勘探地质剖面图与本次数值模拟模型结果进行比较（见图4）。从模拟地质剖面中可以看出，研究区西部的第四系地层岩性以细砂类、粗砂类为主，向东粘粒物质的含量逐渐增大，逐渐变成砂土层与粘土层交错沉积，而实际地质剖面的西南部以细砂、中细砂及中砂为主，向北东粒度逐渐变细，岩性变为细砂、中细砂夹多层粉土、粉质粘土，这说明模型模拟的结果与实际情况基本相符,所建立的三维地质模型能够比较真实地反映研究区的第四系地质结构。

图4 序贯指示模拟剖面与实测地质剖面对比图

4 结论

利用GOCAD的序贯指示模拟模块建立了西辽河平原第四系三维地质模型，真实、直观的反映了区内第四系的地质结构，解决了以往利用二维剖面和等值线图不能真实、直观的反应西辽河平原第四系地层结构的问题，同时为该区分析水文地质条件、建立水文地质概念模型和地下水流动态评价模型提供了地质依据。