基于平均熵值的多点地质统计建模参数优选方法

尹舒祚， 喻思羽， 李少华， 杜 川， 金 浩

(长江大学地球科学学院， 武汉 430100)

1962年，法国数学家Matheron[1]提出了区域化变量的概念，创立了地质统计学。传统的地质统计建模方法[2]，如克里金法使用变差函数对空间场进行估值，但是克里金法作为确定性建模方法，对具有结构性和随机性的复杂地质特征表达存在一定的局限性。斯坦福大学Journel将随机性特征引入地质统计学储层建模，提出了随机模拟方法。1992年多点地质统计学(multi-point geostatistical simulation, MPS)被提出，并且以2002年提出的Snesim算法为关键点，迅速成为随机建模的研究热点[3]。使用“训练图像”的多点地质统计建模方法可很好地表达复杂地质变量的空间结构性，综合了基于象元和基于目标方法的优点[4]，克服了传统基于变差函数地质统计学在刻画储层形态结构特征方面的不足，也可以很好忠于硬数据，展现了良好的应用前景。目前多点地质统计学的建模方法有迭代法和非迭代法两类[5]，迭代法存在难以满足收敛条件、计算速度过于缓慢等问题，导致难以在实际建模过程中得到应用，故逐渐被非迭代法取代。非迭代法包含了Strebelle[6]提出的基于概率的Snesim算法，Arpat[7]提出的基于数据样式的Simpat算法，以及之后发展的Filtersim[8]、DS-MPS[9]等多种多点地质统计算法。

多点地质统计建模方法已经在众多领域得到应用，研究学者[10-14]应用不同的多点算法对三角洲前缘微相、河流相储层分布、储层沉积相及数字岩心等领域进行地质模型的建立和探讨，证实了多点地质统计建模方法的实用性和可行性。使用MPS的一个重要问题是建模质量不仅取决于所使用的训练图像，还取决于一系列算法参数[15]。建模参数作用于空间多点数据事件的提取与重建完整流程，因此选用建模参数是否合适不仅影响建模效果，而且存在建模效率及内存开销[16-17]问题。为此，中外相关学者做了大量的研究，包括以弯曲河道作为训练图像，利用Snesim算法对重要的建模参数(条件数据点最大数量、搜索邻域、多重网格等)进行了敏感性分析[18-20]，并对如何合理设置相关参数提出了建议。Straubhaar等[21]在Impala中用列表结构及索引树取代搜索树，实现了算法计算并行化，提高了模拟效率；喻思羽等[22]通过深入剖析Simpat算法原理及特点，针对其在平衡计算效率和内存空间问题上存在的不足，提出了基于样式降维聚类的多点地质统计建模算法；Yang[23]和Honarkhah等[24]对最佳模板大小，利用熵提出了不同的模板优选方法。Meerschman等[25]对并查集(disjoint-set，DS)算法的综合敏感性分析，提供了DS参数设置的指南；Baninajar等[26]提出了MPS自动参数优化器，这是一种基于随机优化的通用方法，可以快速逼近任何MPS方法和不同类型设置的最优参数。Bai等[27]提出了一种估算Snesim最小多重网格数量的定量方法，使用空间关联度来分析寻找训练图像具有统计意义的最大规模的结构信息，以此来估计多重网格的合适取值。对建模者在训练图像优选上存在的主观标准，近几年有学者[28-30]提出了客观的训练图像评估和优选方法。

每个MPS算法都有一组特定的参数，这些参数设置对多点地质建模中模拟效果的好坏至关重要，需在建模之前进行选择。虽然MPS算法在最佳参数选择上已有一些筛选方法，但包含了人为主观因素，缺乏一定的准确客观评价原则。到目前为止，寻找最佳参数的常见做法是进行灵敏度分析，但是该方法过于烦琐。为了解决这一问题，现提出基于平均熵值的多点地质统计建模参数优选方法。在地质统计学中，熵是一种随机性的统计度量，可用来表征训练图像的纹理，通过计算不同尺寸模式分解下每个模型的平均熵序列，并将训练图像和模拟图像的平均熵曲线差异用Hsim相似度进行量化表示，得到基于平均熵值的空间相关性评价指标与建模参数的关系曲线，通过对量化数据进行分析，找到Hsim相似度曲线趋于平稳的拐点值，将该值作为最佳模拟参数。

1 基于平均熵值的多点地质统计建模参数优选方法

1.1 熵值计算原理

熵的概念由德国物理学家Clausius于1865年提出，最初是用来描述“能量退化”的物质状态参数之一。随着统计物理、信息论等一系列学科的发展和完善，熵在不同领域得到了重要应用。在地质统计学中，熵作为一种随机性的统计度量，可用来表征训练图像的纹理。根据香农第三定理——信源编码定理，香农的熵是衡量消息中所包含的信息，将该概念应用于训练图像上，可以确定所需的最小信息，以便可靠地表示由二进制数字编码的模式。所以熵在模板大小优选方面，通过静态增加训练图像中的模板大小，可以观察到熵值会先急剧增加(编码或压缩图像的平均信息数在增加)，之后随着模板大小增加到表示训练图像平稳特征的最佳大小，熵值会以很慢的速度趋于平稳(编码具有一些平稳特征的图像所需信息近似训练图像本身的平稳性特征)。

熵在模板大小选择中，假设某一特定对象(如泥岩相)在训练图像上大量重复，最佳模板大小将是捕获泥岩对象的模板大小。当模板大小逐渐增加时，会有如下情况：①模板大小低于最佳值，提取模式中包含的信息量将增加；②模板大小达到最佳值，作为泥岩对象的训练图像特征被很好地编码在图像中；③模板大小超过最佳模板，提取模式中的信息量不再增加。

根据上述概念，随着模板大小增加，以熵值为特征的信息量逐渐增加，当信息量增加到与训练图像相似时，熵值趋于稳定。尺寸大小为nx×ny×nz的图像熵值计算公式为

(1)

式(1)中：K为随机变量的可能结果数；pi为概率质量函数(即直方图)。熵值越高，意味着更多的随机性。

1.2 实现步骤

对于任意随机模拟方法而言，建模参数值发生递增或者递减变化时，生成随机模型的表征变化属于渐变，而不是突变。提出的建模参数优选方法建立在此前提条件下，以多点地质建模参数——数据样板尺寸为例，随着数据样板尺寸逐渐增大，生成随机模型和训练图像的空间相关性和几何形态越来越相似。

采用有序递增排列的建模参数集模拟一组随机模型集，计算在不同尺寸模式分解方式下训练图像和所有随机模型的平均熵，并对相同建模参数值实现的100个随机模型模式均熵求均值，利用基于相似性度量的高维数据聚类Hsim算法描述随机模型与训练图像(training image, TI)的平均熵曲线差异，选取拟合曲线上开始趋于平稳的拐点所对应参数值作为优选参数。以建模参数“数据模板尺寸”为例，图1为基于平均熵值的建模参数优选方法具体流程，具体步骤如下。

步骤一输入训练图像TI，给定MPS关键参数——数据模板尺寸的集合C={1,2,…,K},集合C为连续单调递增数值，K表示数据模板大小。

步骤二采用MPS，使用TI和C生成随机模型集M，其中第k个参数Ck对应随机模型集定义为Mk，每个Mk中包含m个随机模型，模型集M中共计有K×m个随机模型。

步骤三选取参数Ck的模型集Mk和TI，计算不同尺寸模式分解下的平均熵序列，计算公式为

(2)

步骤四计算M模型集K组模型的模式平均熵与训练图像的平均熵曲线差异，依次比对模型和训练图像平均熵值曲线和Hsim相似性，空间相似性计算公式为

(3)

高维数据相似性度量函数Hsim()计算公式为

(4)

式(4)中：S(X,Y)为对象X和Y在空间上的相似程度，其值越大，表明数据之间的相似性程度越高；d为空间维度。

步骤五k每次增加1，若k≤K-1，则进入步骤三、步骤四，否则进入步骤六。

图1 优选方法流程图Fig.1 Flow chart of the preferred method

步骤六建立参数k与相似度的拟合曲线，分析并选取曲线变化趋于平稳时的拐点参数值k作为最优参数。随着建模参数值的递增，基于建模参数值的随机模型与TI相似度逐渐增加，并达到难以人为区分的程度，即当建模参数大于最佳值后，随机模型的空间结构特征变化趋于稳定。

2 Simpat建模参数优选实例

多点地质统计储层随机建模方法依赖于训练图像，训练图像作为地质先验模型，其质量的好坏对模拟结果至关重要。使用基于样式的多点地质统计算法Simpat以channel模型作为训练图像进行非条件模拟，对Simpat的重要建模参数——数据模板尺寸开展研究。

图2 channel训练图像Fig.2 The channel training image

图2中，channel为河流相数字模型，包含砂岩和泥岩两种相类型，训练图像大小为101×101。使用Simpat对TI进行非条件模拟，针对模拟模板尺寸参数的优选问题，设置数据样板尺寸大小从1递增到15，在每个数据样板参数下，模拟生成100个随机模型。如图3所示为不同模板大小生成的随机模型，随着模板尺寸逐渐增加，生成随机模型的河道连续性快速改善，形态表征与TI越来越相似，当模板大小增加到3或者4后，河道随机模型的纹理表征改善越来越低，人为很难区分模拟结果的纹理表征，若通过人为主观视觉来确定改善效果的临界值，给定一个最优参数，会增加建模过程中的不确定性。

根据图1的参数优选方法流程，采用式(2)计算每个随机模型和训练图像在不同尺寸分解模式下的平均熵序列，图4为模板大小为1～15的随机模型与TI的平均熵序列拟合曲线，随着模板大小增加到3后，随机模型平均熵序列和训练图像拟合曲线逐渐接近，根据地质统计学中熵值的概念，表示随机生成模型的空间形态结构达到稳定，即模拟结果和TI越来越相似。

再求相同模板参数模拟的100个随机模型平均熵序列均值，由式(4)计算训练图像channel和随机模型平均熵的Hsim相似度，来量化表示两者之间的差异度。如图5所示，随着模拟样板尺寸从1开始增加，随机模型和TI的相似度逐渐增加并趋于稳定，即100个随机模拟结果的平均熵和TI的平均熵曲线差异越来越小并趋于稳定。当模拟样板尺寸达到3时，空间相关性评价指标逐渐平稳，模拟结果质量不再有明显提升，从而确定出本例最优模板参数大小。

图3 基于Simpat不同模板尺寸的channel随机模型Fig.3 The channel random model based on different template sizes of Simpat

图4 随机模型和TI的平均熵序列拟合曲线Fig.4 Random model and TI mean entropy sequence fitting curves

图5 基于平均熵值的Hsim相似度变化曲线(channel)Fig.5 Hsim similarity change curves based on average entropy(channel)

使用多点地质统计算法Simpat对另一训练图像circle进行非条件模拟，同时对Simpat的重要建模参数——数据模板尺寸开展研究。图6为circle数字模型，包含砂岩和泥岩两种相类型，训练图像大小为101×101。

图6 circle训练图像Fig.6 The circle training image

使用Simpat对TI进行模拟，模式实现流程和参数设置同河道训练图像，图7为模板大小1-15生成的circle随机模型。根据图1参数优选方法流程，随机模型、熵值计算、相似度计算及量化拟合熵值差异曲线同河道训练图像，得到Hsim相似度曲线如图8所示，同channel，随着模板大小增加，随机模型和训练模型空间相似度越来越相似，由空间相似度评价指标可得模板大小3为最优参数。

图8 基于平均熵值的Hsim相似度变化曲线(circle)Fig.8 Hsim similarity change curve based on average entropy (circle)

图7 基于Simpat不同模板尺寸的circle随机模型Fig.7 The circle random model based on different template sizes of Simpat

3 结论

(1)选取最佳参数适用于离散的训练图像，模式平均熵值计算采用单点计算代替两点熵值计算，节约了计算时间,提高了参数优选效率。

(2)结合多点地质统计随机模型空间特征和建模参数的相关性认识，建模参数值在递增或递减时，生成随机模型特征属于渐变而非突变，提出一种基于平均熵值的多点地质统计建模参数优选方法。以建模参数——模板大小为例，随着模板尺寸增加，随机模型与TI的空间结构、形态特征越来越相似，利用模式均熵来量化其空间特征，结合Hsim相似度度量方法对随机模型和TI模式均熵差异进行定量化分析，给出客观的优选参数值，对多点建模参数选择和优化建模效果具有实际意义。