基于细胞自动机的断控岩溶过程数值模拟方法

刘彦锋, 张文彪, 段太忠, 赵华伟

(中国石化石油勘探开发研究院， 北京 100083)

随着塔河油田勘探开发的不断深入，发现在塔河油田托甫台以及顺北地区缝洞储集体的发育与大型走滑断裂之间具有较好的匹配关系，且溶蚀发育具有向“纵深”发展的趋势，呈现出“断控岩溶”特征。据此，鲁新便等[1-2]提出“断溶体”圈闭的概念，并通过多年的矿场实践总结出断溶体油藏高产井主要沿大型断裂带分布的特征。近年来，中国一些专家学者针对断溶体的描述、分类、地震预测及三维表征等做了较多的工作。周文等[3]根据断裂与溶洞发育关系，提出“控洞断裂”概念，认为塔河地区大-中型断裂对岩溶洞穴的控制作用明显，奠定了断控岩溶作用的基础。商晓飞等[4]根据断控岩溶的发育特点，按照溶蚀强度的变化划分为不同的溶蚀相带，指出溶蚀过程可以根据主控因素进行追踪和描述。鲁新便等[5]、胡向阳等[6-7]提出分级相控岩溶建模的思路，采用多尺度分类建模的思路，将不同发育规模的缝洞储集体进行融合。已有关于断溶体的描述及表征主要集中在对断溶体圈闭形态的刻画、连通性分析、成因模式总结等，对于断溶体的地质成因过程方面主要停留在定性描述阶段，缺乏一种定量化的模拟方法来推动断控岩溶过程的定量分析。已有岩溶缝洞型储层地质建模技术，目前仍主要以传统的地质统计学和地震雕刻法为主，尚缺乏一种基于岩溶形成过程的模拟技术来刻画储层的分布形态。

关于岩溶过程模拟，已有研究多数集中在小尺度孔隙的溶蚀过程物理模拟和数值计算，缺乏对大尺度特别是油藏尺度溶蚀孔洞发育规律的模拟研究。何治亮等[8]、彭军等[9]根据深层高温高压的特点，通过物理模拟实验的方式探索了白云岩与灰岩溶蚀速率的变化，指出随着温度压力增加，溶蚀速率均明显增强的特征。范维等[10]通过室内试验研究差异溶蚀作用，表明除了温度、压力和水溶液介质等因素的，碳酸盐岩因为岩石成分和结构的差异也会被表现出不同的溶蚀速率。Kang等[11]提出基于格子玻尔兹曼数值模拟方法进行孔隙级别的溶蚀-沉淀过程模拟，在溶蚀过程数值模拟方面迈出了重要一步，但该方法还主要停留在微观孔隙的溶蚀模拟阶段，对于大型溶蚀孔洞级别的模拟未曾涉及。

细胞自动机，作为一种自然界中复杂系统的完全离散化数学模型，非常适合描述在离散时间维上演化的动力学过程。Burgess[12]以碳酸盐岩为例，分别考虑其平面演化及纵向生长的总体特征，并加入沉积物搬运及海平面变化等因素，采用细胞自动机的基本规则，再现了碳酸盐岩沉积演化过程；Liang等[13-14]提出了基于细胞自动机的地貌动力学模型模拟河流-三角洲沉积过程，演化规则考虑了水体、沉积物、搬运、剥蚀等因素。说明细胞自动机方法在沉积物演化方面具有较好的适用性，表明溶蚀改造过程也可借鉴该方法进行数值模拟研究。

基于以上，现将细胞自动机方法引入到碳酸盐岩溶蚀过程模拟研究，在断控岩溶过程主控因素分析基础上，设计溶蚀过程基本规则，形成基于细胞自动机的断控岩溶数值模拟方法。

1 断控岩溶过程地质特征

断溶体的形成过程是大气水沿近垂直的断裂系统下渗至碳酸盐岩地层，岩溶水优先沿宽度大的主干断裂下渗，同时以较低的速度向侧向渗流，表现为以深大断裂带为核心发生溶蚀扩大作用[15-16]。基于大量野外露头、钻井、岩心、测录井及地震资料，分析断溶体发育位置、形态、规模，与大型断裂、裂缝、地层岩性之间的控制作用关系，确定断控岩溶过程划分为4个阶段(图1)。

图1 断控岩溶演化过程示意图Fig.1 Schematic diagram of fault controlled karst evolution process

(1)断溶初期：受走滑构造运动影响，形成从地层深部到浅部的立体断裂结构，走滑断裂纵向延伸长度较大，在主干断裂附近会伴随一系列小规模的断裂，整体呈“V”字形，是断控岩溶过程最重要的格架基础。

(2)断溶早期：大气水沿近似垂直的断裂系统下渗，优先沿着多组裂缝交汇处、裂缝宽度大、密度大的区域下渗，裂缝的扩溶作用是断控岩溶的最重要表现形式，断溶体发育形成雏形。

(3)断溶中期：随着岩溶作用增强，孔洞规模逐渐加大，裂缝不断溶蚀而拓宽，在裂缝密集的地方，孔洞开始连片、合并，形成纵向连通的溶蚀孔洞和小规模溶洞；除了断裂系统，背景岩性也对溶蚀过程有一定影响。

(4)断溶晚期：溶洞继续扩大，周围溶蚀孔洞也继续增多增大，最后演化至现今大型洞穴，是岩溶最发育的区域，周围发育大规模溶蚀孔和溶蚀扩大缝，向外逐渐过渡到裂缝储集区，最后伴随上覆地层沉积覆盖，地表水引起的溶蚀作用逐渐停止。

2 细胞自动机岩溶模拟方法

2.1 细胞自动机基本原理

细胞自动机是一种离散迭代演化模型，可用于模拟复杂动态系统随时间的变化，具有运算速度快、程序编码简单的优势[17-19]。该演化模型由无限(或有限)个相同的网格组成，每个网格处于一种有限状态，每个网格在t时刻的状态由t-1时的有限集合的状态确定，且每次演化时，每个网格遵从同一规则一齐演进，大量细胞通过领域细胞的简单相互作用而构成动态系统的演化，它可以通过简单的演化规则模拟宏观的复杂系统。目前细胞自动机已成为一个应用广泛的研究领域，在图像处理、加密、压缩、计算机仿真模拟等方面都有广泛的应用[20-22]。

针对特定的物理过程需要在细胞自动机的基本假设下设计相应的演化模型。细胞自动机通过刻画最小单元的局部作用及其多个最小单元的相互作用模拟实际物理过程。一个完整的细胞自动机模型由细胞、细胞空间、邻域和细胞规则4部分组成。细胞是离散的组成细胞自动机的最小单元；全部的细胞组成了细胞空间，代表所描述的整个物理系统；邻域表示相互作用的局部细胞团体；细胞演化规则刻画邻域内细胞的作用机理，是细胞自动机的灵魂。对特定物理过程的细胞自动机模拟的好坏关键在于演化规则是否反映了客观事物内在的本质特征。

2.2 岩溶过程细胞自动机建模

细胞自动机模拟断溶体储层的形成以走滑断裂形成的裂缝网络为初始条件，模拟地层在有一定酸浓度的大气淡水在地层岩性差异、裂缝网络和重力等因素的共同作用下的演化过程。细胞自动机采用自下而上的方法建模，对溶蚀过程建模的过程主要包括确定溶蚀过程的性质、对断控岩溶系统进行网格状分割、确定细胞的初始状态、设计溶蚀过程的演化规则。

(1)确定溶蚀过程的性质。细胞自动机最早提出时用于对自组织系统的模拟，不受外界因素影响，只考虑局部因素，目前细胞自动机已经可以用来描述含有外部因素的系统[23]。断溶体形成过程中不断有岩溶水从地表向下供给，因此，溶蚀过程是具有外部因素影响的三维动态系统。

(2)断控岩溶系统网格状分割。采用地质模型网格的方法对系统进行三维结构化网格划分，整个网格系统组成细胞空间，每个网格代表一个待模拟的细胞单元，后续描述中的网格即代表细胞，演化规则考虑三维邻域状态。

(3)确定细胞初始状态。初始状态是走滑作用下形成的断裂系统，该状态由裂缝网络和基质网格组成，裂缝发育的网格f=1表示已溶蚀，基质处溶蚀相f为0，表示未溶蚀。溶蚀状态的取值只能为0或1，中间状态通过下面叙述的其他参数表示。

(4)设计溶蚀演化规则。细胞自动机演化规则的确定要考虑细胞邻域、响应、属性及时间和空间上的相互影响。根据断溶体地质研究认识，断控岩溶以断裂破碎带为核心发生溶蚀扩大，其溶蚀程度的差异受断裂分布、裂缝密度、地层岩性、大气淡水或底部热液接触时间等综合因素影响。由于这些因素涉及地质历史时期的地质过程，参数多，实际测量难度大，且量纲不同，因此，采用基于概率的这种无量纲方法设计溶蚀演化规则。

溶蚀演化规则的基本假设：以自身网格为基础，三维空间相邻的网格为邻域，每个网格的属性有酸浓度c、溶蚀状态f、连通概率pd、不同溶蚀状态接触时连通概率pf、裂缝发育程度pp、背景岩石相的溶蚀容易程度pr，所有参数均为0～1的概率值。

演化规则：如果待模拟网格周围的酸浓度比自己大，则该网格的浓度变大，周围网格的浓度变小，同时满足质量守恒要求。酸的扩散考虑不同方向的连通概率、网格的溶蚀状态、裂缝发育概率、背景岩石相溶蚀容易程度。通过不同类型、不同方向的系数体现不同类型地质因素对溶蚀过程的影响，基于这些基本规则可以设计具体的溶蚀过程算法。

2.3 溶蚀过程算法设计

假设模拟的总时间步为N，前N1个时间步地层处于暴露状态，通过大气降水向地层表面提供含有一定浓度酸的岩溶水，从N1到N时间步为非暴露，不再有地表岩溶水供给，地层内部的岩溶水继续溶蚀和演化，一直达到平衡状态，不再溶蚀。

为了描述具体的细胞自动机演化算法，需要进一步定义一些参数。假设某个网格酸浓度为c0(无量纲方式，最大值为1，最小值为0)，溶蚀状态为f0，累积酸溶蚀为acu，其邻域为前、后、左、右、上、下6个方向网格，上部网格酸浓度为c1，下部网格浓度为c2，左侧网格酸浓度为c3，右侧网格酸浓度为c4，前方网格酸浓度为c5，后方网格酸浓度为c6，上述几个方向的连通概率分别为pd1、pd2、pd3、pd4、pd5、pd6，上述几个方向不同溶蚀状态接触时的连通概率分别为pf1、pf2、pf3、pf4、pf5、pf6(图2)。不同方向的连通系数和不通过溶蚀状态接触的连通概率需要预先给定，由于重力的影响，通常设置为pd1>pd3=pd4=pd5=pd6>pd2，背景岩石相影响的溶蚀概率为pr，已溶蚀与未溶蚀细胞接触时连通概率远小于已溶蚀与已溶蚀细胞接触时的连通概率。每个时间步，暴露在地表的网格酸浓度重置为1，表示有充足的酸供给量，然后其他网格按演化规则同步更新。任意网格在时间维度上累积酸浓度表示为该网格的溶蚀程度，是模拟的主要目标。

图2 不同方向酸浓度、连通概率及相接触连通性示意图Fig.2 Schematic diagram of acid concentration, connectivity coefficient and contact connectivity in different directions

演化的核心规则是在每个时间步每个网格计算酸浓度的增加量dc并更新累计酸溶蚀和溶蚀相f。酸浓度增加量的初始值为dc=0，依次考虑6个方向的浓度影响：

如果c1>c0，则dc+=pd1pp1pf1pr1c1，c1=(1-pd1pp1pf1pr1)c1；

如果c2>c0，则dc+=pd2pp2pf2pr2c2，c2=(1-pd2pp2pf2pr2)c2；

如果c3>c0，则dc+=pd3pp3pf3pr3c3，c3=(1-pd3pp3pf3pr3)c3；

如果c4>c0，则dc+=pd4pp4pf4pr4c4，c4=(1-pd4pp4pf4pr4)c4；

如果c5>c0，则dc+=pd5pp5pf5pr5c5，c5=(1-pd5pp5pf5pr5)c5；

如果c6>c0，则dc+=pd6pp6pf6pr6c6，c6=(1-pd6pp6pf6pr6)c6。

更新c0、acu、f：

c0+=dc，acu+=c0；

如果c0>1，则f+=1。

基于细胞自动机的断溶体溶蚀模拟方法的步骤总结如下(图3)。

图3 基于细胞自动机的断溶体模拟主要步骤Fig.3 Main steps of fault-controlled karst simulation based on cellular automata

步骤1设置初始溶蚀状态、各个方向连通系数、各个方向不同溶蚀相接触时的连通性、每个网格的裂缝发育密度，模拟总的时间步数N，以及暴露的时间步数N1，其中N1

步骤2在时间步t，设置表层网格的酸浓度，如果时间步t

步骤3生成随机路径。

步骤4沿随机路径访问所有网格，如果是表层网格跳过，则否执行下列步骤。

步骤5根据上述演化规则计算该网格的浓度c0，并更新累计酸溶蚀acu及溶蚀状态f。

步骤6判断是否访问完所有节点，如果是则执行步骤7，否则再次执行步骤4～步骤6。

步骤7判断时间步是否等于N，如果是则结束模拟，否则继续执行步骤2～步骤7。

3 应用实例

3.1 地质背景

研究聚焦在塔河油田10区西南部与托甫台区东北部交汇处的一个典型断溶体单元，深大断裂与局部构造共同控制，储集体整体发育，规模大，油气富集程度高。该井组经过多年的开采，表现出累产高、水体能力强、连通关系复杂的特征。基于前期刻画的储集体形态，虽对井间关联性有了一定的认识，但仍存在许多目前地质认识无法解释的动态现象。为进一步挖掘剩余油，优化开发生产对策，开展溶蚀过程模拟有助于明确缝洞储集体空间配置关系和连通样式，为方案调整提供地质依据。

3.2 溶蚀过程模拟

地质模型维数为170×251×410，网格大小为50 m×50 m×1 m。通过地震蚂蚁体、断裂系统分析等手段获取离散裂缝网络，并把离散裂缝粗化到网格系统作为初始状态，设置裂缝所穿越网格为溶蚀相，其他位置为非溶蚀相，然后通过基于细胞自动机的方法模拟断溶体的溶蚀过程演化。

设置总的模拟时间步为N=300，前期暴露时间步为N1=250，后期50个时间步处于埋藏状态。设置不同方向的岩溶水溶蚀性连通概率，pd1=0.9、pd2=0.001、pd3=0.1、pd4=0.1、pd5=0.1、pd6=0.1，设置溶蚀相与周围网格的连通概率pf1=0.9、pf2=0.01、pf3=0.5、pf4=0.5、pf5=0.5，pf6=0.5，设置非溶蚀相与周围网格的连通性均为0.005。模拟结果基本可以再现岩溶水在断裂系统等因素控制下对碳酸盐岩地层的溶蚀过程，溶蚀程度演化过程复合地质规律，可以清楚反映不同溶蚀程度储层的发育位置，对认识优势储层分布、储集体连通关系具有重要意义(图4，图5)。

图4 不同时间步剖面演化特征Fig.4 Profile characteristic of different evolution step

图5 不同模拟时间步平面演化特征Fig.5 Plane characteristic of different evolution step

4 结论

(1)根据矿场实践特征总结了断控岩溶演化规则，并首次将细胞自动机的方法引入断控岩溶过程模拟，设计规则均基于实际地层演化特征，能够客观反映形成过程本质，具有较强的实用性。

(2)传统地质统计学方法相比，该方法考虑了成因过程，对于井资料少、地震资料品质低的断溶体发育区，具有预测准确度高的优势。

(3)初始条件及溶蚀参数对模拟结果有重要影响，为了使模拟结果与实际数据吻合，需要借助反演或优化的方法，进行下一步深入研究。