基于FILTERSIM算法的油藏沉积微相模拟研究

孟欣然，梁堰波，孟宪海，杨 钦

(1.北京航空航天大学软件开发环境国家重点实验室，北京100191;2.北京航空航天大学计算机学院，北京100191)

0 引言

地质统计学建模方法从传统上分为两个流派:基于象元的两点地质统计学方法和基于目标的地质统计学方法［10］。第一种方法是以变差函数为工具，利用两点变差函数来反映储层空间结构，能够很好地忠实于实际数据的约束，但很难表征复杂的几何形态 (如弯曲河道)和空间结构，因此更适用于海相而非陆地相的地质结构。而第二种方法则更符合地质工作者的经验，但是较难与实际数据相结合。多点地质统计学 (Multiple-point Statistics，以下简称MPS)［1］模拟方法结合了基于象元的模拟方法和基于对象的模拟方法两者的优点。首先，多点地质统计学是一种改进的基于象元的方法，继承了基于象元的模拟方法的灵活性，每一步只对一个像素点进行建模，因此很容易满足控制点的要求;同时，由于引入了训练图像，使得建模结果的曲线连续性有了很大提高，尤其擅长模拟复杂的河道地质结构，更符合实际的需求。

MPS第一个比较重要的算法是SNESIM算法［2］，此处不再赘述。SNESIM算法只能用于油藏沉积微相的建模，而不能用于数值连续的油藏属性建模，同时，如果沉积相类型的数目过多，SNESIM算法的效率不甚理想。当训练图像较复杂或者分类较多时，SNESIM算法对内存的需求也是巨大的。因此，2006年，一种新的 MPS算法——FILTERSIM［3］诞生了，该算法能够在合理的内存需求下完成沉积微相建模和属性建模。FILTERSIM算法引入了一系列不同方向的线性滤波器 (linear filter)把训练图像进行分类，类似的训练图像模式被归为一类，并以这些模式的均值作为这一类的标志，称为这个类的“原型” (prototype)，在模拟过程中，原型与当前约束数据最相近的那个类被选为当前网格点的实现。SNESIM算法把所有备选的训练图像都存储在一个搜索树里，因此对内存要求较高，而FILTERSIM算法只存储每个训练图像的中心点，因此，极大的降低了内存需求。

同时，油藏最重要的属性参数，如孔隙率、渗透率、含油饱和度等都是连续数据，SNESIM算法无法实现模拟连续数据的需求，而FILTERSIM算法可以分为FILTERSIMCATE和FILTERSIM_CONT，后者可以对连续数据进行精确模拟，弥补了SNESIM算法的不足。

1 FILTERSIM算法概要

FILTERSIM算法主要分为3步:滤波分 (filter score)计算、模式分类和模式模拟。此外，对于已经进行过地震数据处理的油藏，FILTERSIM算法可以综合地震处理数据进行上述3个过程。

1.1 计算滤波分

FILTERSIM算法一个非常重要的概念是搜索模板(search template)τn，用来搜索训练图像。搜索模板的大小用3个整数表示(nx，ny，nz)来表示，并且三个参数必须都是奇数。搜索模板的每一个节点都记录一个其与中心节点的偏移量，ua=u+ha(a=1，2，…，n)。如图1所示，ha代表了 (a)中除去中心点u之外的其他80个点较中心点的偏移向量、(b)中除去中心点u之外26个点较中心点u的偏移向量。将搜索模板应用于训练图像就获得了一个图像模式集。图2显示了将搜索模板应用于训练图像一个点获得的一个图像模式的过程。

FILTERSIM算法里的滤波器 (filter)定义为一系列数据的集合TJ，集合的元数为 J，TJ={u0;hj，j=1，…，J}，集合中的每个节点uj都用与中心节点u0的偏移向量hi=(x，y，z)i来表示，同时赋予该节点一个特定的加权量(weight)fj。偏移向量的元素x、y、z都是整数。为了简化存储，FILTERSIM提出了一个滤波分 (filter score)ST(u)的概念，公式是

其中，pat(u+hi)是节点的属性值，J=nx×ny×nz。一套滤波器是无法捕捉训练图像的全部特征的，因此，一个滤波器集合F共同用于发掘训练图像的特征。上面的公式可以改写为

因此，训练图像的存储规模由nx×ny×nz缩减为F。比如说，一个11×11×3像素的3D训练图像可以用FILTERSIM算法生成的9个F值来表示。如图3所示。

图3 利用滤波器简化过程

对于连续的训练图像 (训练图像像素的值代表属性值)来说，F个滤波分直接用来代表训练图像元素本身。对于分类型训练图像 (像素值只有几种类型，代表不同的沉积微相)，假设沉积相种类为K，训练图像首先要转化为一系列二进制指示向量 Ik(u)，k=0，…，K －1，u∈Ti

图4给出了二维图形常用的六个滤波器，公式如下:

(a)南北方向均值，中心行权值较重

(b)东西方向均值，中心列权值较重

(c)南北方向梯度，边缘行权值较重

(d)东西方向梯度，边缘列权值较重

(e)南北方向曲率

(f)东西方向曲率

图4 适用于二维训练图像的滤波器组合

1.2 图像模式分类

将F个滤波器应用于K个沉积相的训练图像可以获得F×K个滤波分图，这个维度远远小于图像模式的像素数，大大减小了维度。相似的图像模式会得到差不多的滤波分，因此，可以将滤波分进行分类，相似的滤波分为一类，以均值作为这一类的原型 (prototype)，记为prot，公式如下

分类策略最常用的有两类，一类是交叉分类 (cross partition)，一类是K-均值聚类 (K-mean cluster)，此处不详细阐述。前者的优点是效率较快，但是分类结果较粗糙。后者是最简单的非监督学习模型，能够获得较好的分类效果。

1.3 图像模式模拟

图像模拟过程就是不断寻找与当前点最相近的模式的过程，具体算法如下:

算法1

(1)利用一定的滤波器获得滤波分

(2)在滤波分空间 (F×K维)进行模式分类

(3)在模拟网格G中定位控制数据

(4)定义一个访问G的随机路径

(5)for路径中的每一个节点u，do

(6)提取以u为中心点的约束数据集dev

(7)找到与dev最相近的父原型protp

(8)if protp存在子原型集，then

(9)找到与dev最接近的子原型protc

(10)从protc中随机选择一个模式pat

(11)else

(12)从protp中随机选择一个模式pat

(13)end if

(14)将pat黏贴到已经模拟的点集中，并冻结dev中的所有点

(15)end for

2 利用FILTERSIM算法进行油藏沉积微相建模

由于SNESIM算法在解决油藏沉积微相方面的能力很强 (前提是沉积相种类较少)，为了评估FILTERSIM算法在沉积微相建模方面的性能，我们对Stanford V数据分别使用SNESIM算法和FILTERSIM进行了沉积微相模拟。如图5所示，Stanford V数据分为4个部分，其中grid_2d是一个用于进行属性建模的100*130*1笛卡尔网格，并且提供了三个地震属性数据，P(crevasse|seis)2、P(mud|seis)2和P(sand|seis)2代表了3种沉积相的分布概率。TI_2D_large提供了算法需要的训练图像，well提供约束的井数据。图6是Standford V数据中的训练图像和地震数据。

图6中，(a)是二维的三相训练图像，其中浅色代表砂质相 (sand channel)，深色代表泥相 (mud face)，小点代表断层 (crevasse);(b)是目标的泥相概率分布图像。

在本实例中，地震处理提供的3种沉积相分布被用于距离计算，前述的算法可以被改写成如下过程 (其中tau模型是一个合成概率场的工具):

算法2

(1)对于随机路径上的每一个节点u，同时从以实现的模拟数据和地震数据中获得约束数据集dev和sdev(后者来源于地震数据)

(2)if dev是空集then

(3)用sdev代替dev，并利用新的dev获取最近的模式

(4)else

(5)用dev获取最近的模式prot

(6)用tau模型将模式prot和地震数据sdev合成一个新控制数据集dev*，作为本地概率分布模型

(7)寻找与dev*最接近的模式，继续模拟过程

(8)end if

搜索模板的大小为11×11×1(搜索模板大小可以自行设置，取决于对精度的需求，一般为训练图像的1/10-1/5)。图7是分别利用SNESIM算法和FILTERSIM算法获得100组结果选取的三个结果。

图7 利用FILTERSIM算法和SNESIM算法实现二维三相沉积微相分类结果

图7 中，(a)、(c)、(e)是由FILTERSIM算法实现的3个沉积微相分类结果;(b)、(d)、(f)是由SNESIM算法实现的3个沉积微相分类结果。

图8是利用FILTERSIM算法和SNESIM算法实现的沉积微相模拟结果生成的mud概率分布图像。可以看出FILTERSIM算法获得的概率模型边缘比较模糊，而SNESIM算法获得的较清晰，但是两者都较好地反映了地震数据揭示的地址概率模型。

图8 利用FILTERSIM算法和SNESIM算法实现的沉积微相模拟结果生成的泥相概率分布图像

图8 中，(a)是由100个FILTERSIM算法模拟结果生成的泥相概率分布;(b)是由100个SNESIM算法模拟结果生成的泥相概率分布。

3 利用FILTERSIM算法模拟连续属性数据

为了验证FILTERSIM算法在连续属性数据模拟方面的性能，我们在Stanford VI［4］数据集上进行了实验。本实验的训练图像是一个80×80×40像素的三维图像，像素值的范围是0-1，代表油藏的孔隙度，如图9所示。孔隙度均值为0.6522，方差为0.078。抽取其中的0.5%数据作为约束数据，如图10所示。利用连续FILTERSIM算法获得的模拟结果之一如图11所示。表1是随机选取的三个结果中孔隙度均值和方差数据。可以看出，无论是图形模式、均值还是方差，FILTERSIM的模拟结果都与训练图像十分接近。

图11 连续FILTERSIM算法获得的模拟结果之一

表1 FILTERSIM算法随机选取的3个结果中均值和方差数据

4 结束语

FILTERSIM算法能够在较少的CPU和内存消耗下获得与SNESIM算法近似的结果，并且能处理沉积相建模和属性建模，这是SNESIM算法无法实现的。然而，在沉积相建模质量方面来看，FILTERSIM算法略逊于SNESIM算法，尤其是在边缘轮廓方面比较模糊。这在井间距离较大或者井数据较少的地质条件下尤其有用。此外，FILTERSIM还能用于连续属性数据模拟，并且模拟结果能较好地反映训练图像的模式 (均值和方差均与训练图像十分接近)。

［1］ZHANG T F.Incorporating geological conceptual models and interpretations into reservoir modeling using multiple-point geo-statistics［J］.Earth Science Frontiers，2008，15(1):26-35.

［2］Boucher A.Considering complex training images with search tree partitioning［J］.Computers ＆ Geosciences，2009，35(6):1151-1158.

［3］ZHANG T F.Filter-based training pattern classification for spatial pattern simulation［D］.Ph D dissertation，USA:Stanford University，2006:11-26.

［4］Remy N，Boucher A，WU JB.Applied geostatistics with SGeMS:A Users'guide［M］.Cambridge University Press， 2009:108-134.

［5］WU J B.4D seismic and multiple-point pattern data integration using geostatistics［D］.Ph D dissertation， USA:Stanford University，2007:163-173.

［6］Chugunova T L，HU L Y.Multiple-point simulations constrained by continuous auxiliary data［J］.Mathematical Geosciences 2008，40(2):133-146.

［7］Arpat G B，Caers J.Conditional simulation with patterns［J］.Mathematical Geology，2007，39(2):177-203.

［8］Sanchez-Brea L M，Bernabeu E.Uncertainty estimation by convolution using spatial statistics［J］.IEEE Transactions on Image Processing，2006，15(10):3131-3137.

［9］LU D T，ZHANGT J，YANGQ，et al.A reconstruction method of porous media integrating soft data with hard data［J］.Chinese Science Bulletin，2009，54(11):1876-1885.

［10］LUO Yang.Using MPSin the river reservoir simulation［J］.Geological Science and Technology Information，2008，27(3):68-72(in Chinese).［骆杨.多点地质统计学在河流相储层建模中的应用［J］.地质科技情报，2008，27(3):68-72.］