引用本文
喻思羽, 李少华, 段太忠, 王鸣川. 基于局部各向异性的非平稳多点地质统计学算法[J]. 物探与化探, 2017,41(2): 262-269.
YU Si-Yu, LI Shao-Hua, DUAN Tai-Zhong, WANG Ming-Chuan. Non-stationary multiple-point geostatistics algorithm base on local anisotropy[J]. Geophysical and Geochemical Exploration, 2017,41(2): 262-269.  
Doi:10.11720/wtyht.2017.2.11
Permissions
Copyright©2017, 《物探与化探》编辑部
《物探与化探》编辑部 所有
基于局部各向异性的非平稳多点地质统计学算法
喻思羽1, 李少华1, 段太忠2, 王鸣川2
1.长江大学 地球科学学院,武汉 430100
2.中石化石油勘探开发研究院,北京 100083
李少华(1972-),男,博士生导师,教授,主要从事地质统计学、地质建模方面的研究和教学工作。Email:534354156@qq.com

作者简介: 喻思羽(1987-),男,博士研究生,现从事地质统计学研究及相关软件研发工作。

收稿日期: 2016-07-25
基金: 国家自然科学基金(41572121); 湖北省自然科学基金创新群体“储层精细表征与建模”(2016CFA024); 国家重大科技专项(2016ZX05033-003-007)
摘要

基于平稳性假设的传统多点地质统计学在非平稳建模中存在不足。提出一种基于局部各向异性分区的多点地质统计建模方法,利用不同方向变差函数变程组合的玫瑰花图定量表征训练图像的局部各向异性,结合经典多维尺度分析和K均值聚类分析对非平稳训练图像进行自动分区,各子区域相互独立,分区内运用传统平稳性多点建模算法模拟,最终实现非平稳性建模。以裂缝网络建模为例,新方法较好地模拟并再现训练图像的非平稳性,对多点地质统计建模算法的非平稳建模具有很好借鉴意义。

关键词: 各向异性; 玫瑰花图; 多维尺度分析; 多点地质统计学; 非平稳性
中图分类号:P613.4 文献标志码:A 文章编号:1000-8918(2017)02-0262-08
Non-stationary multiple-point geostatistics algorithm base on local anisotropy
YU Si-Yu1, LI Shao-Hua1, DUAN Tai-Zhong2, WANG Ming-Chuan2
1.College of Geosciences,Yangtze University,Wuhan 430100,China
2.Petroleum Exploration and Production Research Institute,Sinopec,Beijing 100083,China
Abstract

Classical multiple-point geostatistics based on stationary assumption has shortcomings in non-stationary modeling application.This paper proposes a new non-stationary simulation based on segmentation by using anisotropy.The new method uses rose diagram consisting of different directions' variogram range to quantify characterization of train image's local anisotropy.In combination with multidimensional scaling and K-means,the proposed method can automatically segment non-stationary train image into several sub-regions which are independent of each other.In the sub-region,classical MPS algorithms are used to simulate sub-regions' geological model which eventually form a whole non-stationary realization.Using fracture network modeling as an example,the proposed method simulated and reproduced non-stationary train image perfectly.The results obtained by the authors provide reference for non-stationary modeling of multiple-point geostatistics.

Keyword: anisotropy; rose diagram; multidimensional scaling; MPS; non-stationary
0 引言

多点地质统计建模方法(multiple-point geostatistics, MPS)是储层地质建模研究领域的前沿[1, 2]。基于空间多点相关性, 多点地质统计建模方法预先用数据样板扫描提取训练图像的数据样式, 即重复观察训练图像所有位置的样本点, 得到训练图像内在的几何形态及非均质特征, 完成未知区域的预测[3]。平稳性假设是经典地质统计学的理论基础[4], 即区域化变量均值、方差必须平稳且与位置无关。因此训练图像必须具备统计平稳性[5]。然而自然界没有绝对的平稳现象, 平稳是相对的, 非平稳是绝对的。构造运动、水进(退)、物源等因素不断变化, 控制着沉积体系朝着非平稳方向演变[6]。典型的非平稳地质现象有冲积扇、辫状河三角洲等。针对非平稳建模, 目前有3种基于非平稳训练图像的建模方法。第1种方法将平稳的训练图像进行旋转、缩放变换, 得到不同角度和尺度的训练图像, 并作为空间独立不相关的多点统计信息库, Zhang[7]和Arpat[8]把这种思路分别运用于SNESIM和SIMPAT算法上, 取得一定效果。但建立若干旋转、缩放变换的训练图像的多点统计信息库, 不仅增加内存开销同时效率较低, 此外旋转、缩放参数的选取较为困难。第2种方法是基于分区方法, 2011年, Mehrdad Honarkhah通过形态学特征, 例如Gabor纹理过滤器[9], 把非平稳的训练图像分割为多个子区域, 然后提取各个子区域内部的空间多点相关性, 再用于相应子区域的模拟。由于纹理过滤器Gabor是二维图像分析工具, 其难以用于三维建模。第3种方法是基于距离函数[10]建立储层模式与位置的关联性, 再现具有空间变化性的非平稳性沉积储层空间分布。由于模式与空间位置的关联函数对距离值的大小非常敏感, 距离过小使随机性较弱, 距离过大难以保证非平稳性的再现, 距离值难以给定。

针对现有方法不足, 笔者提出基于局部各向异性自动分区的非平稳多点地质统计学算法ASNSIM(anisotropy segmentation-based non-stationary simulation)。通过引入玫瑰花图定量表征局部各向异性, 对训练图像中所有玫瑰花图进行降维和聚类分析, 将非平稳训练图像自动分割成若干平稳的子区域。子区域内部沿用经典多点建模算法完成模拟计算。以扇三角洲和离散裂缝网络模型为例, 新方法较好再现非平稳训练图像的空间变化特征。

1 ASNSIM的基本原理

经典多点地质统计建模方法难以再现非平稳训练图像的形态展布特征, 如图1所示, 图1a为扇三角洲的训练图像[11], 训练图像网格维度为200× 200, 对比A区(扇根)和B、C和D区(扇缘的不同位置), 扇三角洲分支河道的宽度和走向发生明显变化。经典多点地质统计建模方法(如Simpat)难以再现训练图像的非平稳特征, 如图1b所示。

真实地质现象通常符合非平稳是绝对的、平稳是相对的客观规律, 局部的平稳随着空间距离增加逐渐变为全局的非平稳。各向异性是反映(非)平稳程度的指标, 理论上不同位置的各向异性相似性越强, 则该区域平稳性越强, 反之非平稳性越强。ASNSIM使用玫瑰花图定量表征训练图像的局部各向异性, 进而对非平稳训练图像进行识别分区。

图1 非平稳的训练图像[11](a)和Simpat模拟实现(b)

1.1 基于玫瑰花图表征各向异性

1965年Matheron提出基于矩估计的变差函数[12], 定义为在方向α 、相距|h|的区域化变量Z(x)与Z(x+h)增量的方差之半[13], 定义公式:

r(x, h)=12Var[z(x)-z(x+h)]。（1）

变差函数能反映空间结构性和随机性[14], 可用于评价区域化变量在某个方向上某一距离范围内的变化程度, 变程是变差函数的主要参数。在变程范围内, 区域化变量具有相关性, 范围之外不再具有相关性。

不同方向的变程能够反映不同方向的空间变异性, 多个不同方向的变程组合称为玫瑰花图[15], 八方向玫瑰花图定义为:

Rose={ai|i=0, 45, 90, 135, 180, 225, 270, 315}, （2）

ai是方向为i的变差函数拟合曲线的变程。变差函数具有方向对称性, 即方位角0° 与180° 的变程值是相同的, 以此类推, 实际中只需计算0° 、45° 、90° 与135° 这4个方位角的变差函数即可。若区域化变量具有各向异性, 玫瑰花图呈现椭圆形, 其长轴方向表示最大连续方向, 短轴为最小连续方向[16]。

1.2 局部各向异性与整体非平稳性

局部区块的玫瑰花图相似度越大, 则它们属于相同子区域可能性越大。如图1a所示, 局部区块A~D的几何形态区别较大。为定量表征这种差异性, 计算局部区块A~D在0° 、45° 、90° 和135° 这4个方位角的实验变差函数[17](图2), 并用最小二乘法拟合其球状模型[18]。局部区块A的玫瑰花图(图2b)在方位角90° 的变程大于方位角0° 、45° 和135° 的变程, 并且与区块B~D相比, 区块A的各向同性较强; 区块B(图2d)玫瑰花图的长轴方向为90° , 表现出较强的各向异性; 区块C(图2f)的玫瑰花图主方位角为0° , 河道走向与B垂直。综上可知局部区块A~D的玫瑰花图特征(图2)与训练图像(图1)表现的非平稳性是一致的。

图2 训练图像A~D区的拟合变差函数曲线及玫瑰花

1.3 基于多维尺度分析降低玫瑰花图维度

玫瑰花图由8个方向变程组成, 直接从8个维度上分析玫瑰花图的相似(异)度较困难[19]。笔者采用经典多维尺度分析把玫瑰花图从高维降为低维, ASNSIM算法将玫瑰花图从8维降为2维, 2维点间距反映了玫瑰花图的相似(异)度。如图3所示, 以局部区块扫描训练图像提取局部空间数据, 局部区块网格大小为51× 51。计算局部采样区8个方向的实验变差函数并拟合出球状模型, 构建玫瑰花图, 训练图像的全部玫瑰花图集合称为玫瑰花图库, 其公式为:

RoseLibrary={Rose(loc1), Rose(loc2), …, Rose(locn)}, (3)

式中:Rose(locn)是训练图像中位置为locn的玫瑰花图, RoseLibrary是训练图像所有局部玫瑰花图组合的玫瑰花图库。

进行多维尺度分析前先要构建玫瑰花图库的距离矩阵(Distance matrix), 矩阵元素是两个玫瑰花图之间的曼哈顿距离, 其公式:

dij=∑k=18|roseik-rosejk|, (4)

式中:dij是第i个和第j个玫瑰花图之间的曼哈顿距离, ros eik是玫瑰花图库中第i个玫瑰花图的第k个变程。表1是根据图3训练图像构建的玫瑰花图库距离矩阵(部分)。

图3 局部区块扫描训练图像示意图

表1 玫瑰花图的距离矩阵(部分)

经典多维尺度分析方法只需要做特征值分解就可以实现 2022, 首先定义一个距离矩阵A, 矩阵A的行索引为i、列索引为j的元素是

aij=-12dij2。（5）

建立伪中心化内积阵

B=HAH, （6）

其中H是中心化矩阵,

H=I-1n11T∈Rn×n, （7）

其中I是n× n的单位矩阵, 1是n× 1的列向量[1 1 … 1]T∈ Rn× 1。然后对矩阵B进行特征值分解

B=VBΛBVTB, （8）

其中, VB是矩阵B的特征向量, Λ B是矩阵B的特征值。对特征值的大小进行排序, 然后选择前d个的特征向量VB, d左乘特征值Λ B, d的平方根即可以获得玫瑰花图库在d维度上最佳构图, d通常为2, 计算如下

Xd=VB, dΛB, d12∈Rd。（9）

其中:Xd是d维欧式空间中的点, 也称之为古典MDS的解。

基于经典多维尺度分析方法对玫瑰花图库距离矩阵降维得到二维欧式空间的感知图(Perceptual map), 如图4所示。感知图的点距反映了玫瑰花图的相似(异)性, 也就反映了局部各向异性的相似(异)性。

1.4 基于K均值聚类的非平稳分区

按照ASNSIM原理, 如图4所示, 感知图的点距越近, 其对应玫瑰花图的相似性越高, 代表其局部各向异性越相似, 相应局部区块属于相同子区域的概率越大。因此对感知图的离散点进行K均值聚类即可自动将非平稳训练图像分割为若干子区域。

K均值是一种基于距离划分的聚类算法[23], 点距越近, 两者相似性越大。该算法认为簇是由距离相近的对象组成, 把最终紧凑且独立的簇作为最终目标。具体步骤:① 从n个点随机选取k个点作为质心; ② 对剩余的每个点测量其到每个质心的距离, 并把它归到最近的质心的类; ③ 重新计算已经得到的各个类的质心; ④ 迭代2~3步直至新的质心与原来质心相等或小于指定阈值, 聚类结束。

图4 基于K-means的二维感知图聚类结果

基于K-means聚类方法对图4的二维感知图进行聚类分析, 感知图被分成k个类簇k-Clusters, 本次研究中k等于4, 其中每个点的分类标记解释为训练图像的子区域标记。把4-Clusters的全部点的分类信息作为子区域标记赋给局部区块, 得到非平稳训练图像分区结果, 如图5所示。基于ASNSIM方法以局部空间各向异性为基础, 综合运用经典多维尺度和K-means聚类分析, 计算效率高, 分区结果真实客观。

图5 基于ASN策略的非平稳训练图像分区结果

2 ASNSIM的实现步骤

ASNSIM增强了多点地质统计算法进行非平稳建模的能力。使用ASNSIM实现非平稳训练图像的自动分区, 子区域内部可采用传统MPS算法完成建模, 如SNESIM、SIMPAT等。图6是基于SIMPAT改进的ASNSIM算法流程。

SIMPAT的核心思想是基于相似度比较训练图像的数据样式与数据事件的相似程度, 用数据样板扫描训练图像获取全部数据样式加入样式数据库。计算数据样式与数据事件之间的相似度, 匹配样式数据库中与数据事件最相似的数据样式并覆盖、冻结数据事件直到模拟完成。

在SIMPAT基础上, ASNSIM对训练图像进行分区处理。模拟子区域的时候, 用于匹配数据事件的不是训练图像全部数据样式, 而是子区域内部的数据样式。基于分区网格, 以数据样板扫描训练图像建立各子区域的样式数据库, 其中子区域样式数据库的数据样式之间具备很强平稳性。

按照图6编制ASNSIM算法的计算机程序。输入训练图像(图1a)模拟2个随机实现(图7), 模拟结果相比图1b 的SIMPAT模拟实现较好地再现了扇三角洲的河道形态及分布变化特征。

图6 ASNSIM算法处理流程

图7 基于ASNSIM算法模拟结果
a— 基于ASNSIM模拟实现1; b— 基于ASNSIM模拟实现2

3 实例研究

以离散裂缝建模为例检验ASNSIM算法对非平稳模型的再现能力。裂缝网络的训练图像(图8a)发育裂缝情况较复杂[24], 裂缝走向具有明显的空间变化性, 左部裂缝走向为北偏西、右部为北偏东, 左部、中部和右部的裂缝产状有明显差异。从多点地质统计学角度上看, 训练图像具较强非平稳性。基于ASNSIM对裂缝网络训练图像做自动分区(分区数目设定为3)。图8a研究区被分成左、中、右三块, 分区范围与图8a训练图像的空间非平稳性较符合。

图8 裂缝网络训练图像[24](a)及基于ASNSIM的分区(b)

对比SIMPAT和ASNSIM算法的裂缝网络模型, 发现SIMPAT模拟实现(图9a)与训练图像差异明显, 难以刻画裂缝发育方向和分布变化规律; 图9b~d是基于ASNSIM的3个随机实现, 模拟结果具有随机性同时较好再现了裂缝网络的空间分布变化规律。

图9 裂缝网络的非条件模拟实现
a— Simpat模拟实现; b~d— ASMSIM模拟实现

4 结论

1) 计算非平稳训练图像中局部区块在8个方向的变差函数变程, 构建定量表征局部各向异性的玫瑰花图。以扇三角洲为例, 扇三角洲分流河道的走向及宽度与玫瑰花图形态有较强相关性。进一步提出基于局部各向异性定量表征训练图像非平稳性的策略, 使用局部区块扫描训练图像提取局部空间数据。计算玫瑰花图, 采用多维尺度分析将玫瑰花图进行降维, 将玫瑰花图投影到二维平面, 最后进行k均值聚类, 同一类簇的玫瑰花图代表的局部区域属于同一子区域, 子区域之间相互独立, 实现非平稳训练图像的自动分区。

2) 以离散裂缝网络建模为例, 用ASNSIM将非平稳裂缝训练图像分割为3个平稳的子区域, 然后提取子区域的样式数据库用以完成子区域内部的模拟。结果表明新方法能够有效将非平稳区域分割为若干相对平稳的子区域, 并且在经典多点地质统计建模方法基础上实现非平稳建模。

The authors have declared that no competing interests exist.

参考文献
[1] 石书缘, 尹艳树, 冯文杰. 多点地质统计学建模的发展趋势[J]. 物探与化探, 2012, 36(4): 655-660. [本文引用:1]
[2] 尹艳树, 张昌民, 李少华. 多点地质统计学原理、方法及应用[M]. 北京: 地质出版社, 2013: 2-4. [本文引用:1]
[3] Honarkhah M. Stochastic simulation of patterns using distance-based pattern modeling[D]. Paloma Alto: Stanford University, 2011. [本文引用:1]
[4] 张仁铎. 空间变异理论及应用[M]. 北京: 科学出版社, 2005. [本文引用:1]
[5] Journel A G. The deterministic side of geostatistics[J]. Journal of the International Association for Mathematical Geology, 1985, 17(1): 1-15. [本文引用:1]
[6] Strebelle, Zhang T. Non-stationary multiple-point geostatistical models in: Leuangthong O, Deutsch C V. Quantitative Geology and Geostatistics[M]. Springer Netherland s, 2005: 235-244. [本文引用:1]
[7] Zhang T. Multiple-point simulation of multiple reservoir facies[D]. Stanford: Stanford University. [本文引用:1]
[8] Arpat G B. Sequential simulation with patterns[D]. Stanford: Stanford University, 2005. [本文引用:1]
[9] Honarkhah M, Caers J. Direct pattern-based simulation of non-stationary geostatistical models[J]. Mathematical Geosciences, 2012, 44(6): 651-672. [本文引用:1]
[10] 尹艳树, 张昌民, 李少华, 等. 一种基于沉积模式的多点地质统计学建模方法[J]. 地质论评, 2014, 1: 216-221. [本文引用:1]
[11] Arpat G B. Sequential simulation with patterns[D]. Stanford: Stanford University, 2005. [本文引用:1]
[12] Matheron G F. The Theory of regionalized variables and its application[J]. école Nationale Supérieure des Mines de Paris, 1971: 5. [本文引用:1]
[13] 王仁铎, 胡光道. 线性地质统计学[M]. 北京: 地质出版社, 1988. [本文引用:1]
[14] 侯景儒, 尹镇南, 李维明, 等. 实用地质统计学[M]. 北京. 地质出版社, 1998: 36-45. [本文引用:1]
[15] 靳松, 朱筱敏, 钟大康. 变差函数在沉积微相自动识别中的应用[J]. 石油学报, 2006, 3: 57-60. [本文引用:1]
[16] 施小清, 姜蓓蕾, 卞锦宇, 等. 以地质统计方法推估上海第三承压含水层渗透系数的分布[J]. 工程勘察, 2009, 1: 36-41. [本文引用:1]
[17] Deutsch C V, Journel A G. GSLIB: Geostatistical software libarary and user's guide[M]. New York: Oxford University Press, 1992: 39-53. [本文引用:1]
[18] 刘焕荣. 实验变异函数及其理论模型的计算与拟合研究[D]. 昆明: 昆明理工大学, 2014. [本文引用:1]
[19] Ding C, He X, Zha H, et al. Adaptive dimension reduction for clustering high dimensional data[C]//IEEE International Conference on Data Mining IEEE Computer Society, 2002: 147. [本文引用:1]
[20] Bécavin C, Tchitchek N, Mintsa-Eya C, et al. Improving the efficiency of multidimensional scaling in the analysis of high-dimensional data using singular value decomposition[J]. Advance Access publication, 2011, 27(10): 1413-1421. [本文引用:1]
[21] Kruskal J B. Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis[J]. Psychometrika, 1964, 29(1): 1-27. [本文引用:1]
[22] Mardia K V. Some properties of clasical multi-dimesional scaling[J]. Communications in Statistics-Theory and Methods, 1978, 7(13): 1233-1241. [本文引用:1]
[23] Kanungo T, Mount D M, Netanyahu N S, et al. An efficient k-means clustering algorithm: Analysis and implementation[J]. Pattern Analysis and Machine Intelligence, IEEE Transactions, 2002, 24(7): 881-892. [本文引用:1]
[24] Boucher A. Considering complex training images with search tree partitioning[J]. Computers and Geosciences, 2009, 35(6): 1151-1158. [本文引用:1]