基于马尔科夫随机场的岩性识别方法

田玉昆，周 辉，袁三一

中国石油大学油气资源与探测国家重点实验室，CNPC物探重点实验室，北京 102249

1 引 言

近年来，中国陆上含油气盆地勘探开发已经日趋成熟，在勘探思路上已经由构造油气藏向岩性油气藏转变.在这种条件下，准确识别岩性、流体，对油气勘探开发具有重要意义.

自1984年Ostrander［1］提出了从地震资料中估算岩石弹性特征的AVO技术之后，通过地震资料进行岩性识别的方法得到了飞速发展.目前，通过地震数据反演得到弹性参数，再通过岩石物理模型的转换得到含流体饱和度、孔隙度等储层物性参数是进行岩性识别、储层评价的常规手段.直接通过弹性参数进行岩性识别的方法主要有两种：一是建立弹性参数量板，即统计各种岩性的弹性参数范围，然后对反演得到的弹性参数进行判断，落于某种岩性量板范围内岩性的弹性参数归于该种岩性.由于不同岩性的弹性参数常常存在一定范围的重叠，所以这种方法的识别精度不高，而且，由于地区间的差异，每次识别都要建立不同的量板，工作量大［2］；二是通过交会图的方式，采用几种对岩性敏感的弹性参数进行交会，选取在交会图上能够使岩性区分开的弹性参数或者岩石物理参数进行岩性的识别，这种方法通常是选取测井资料中已有的弹性参数进行交会，然后用交会的结果识别地震反演剖面上的岩性，同样的，这种方法也会受限于弹性参数重叠的情况，如果通过交会图，找不到合适的参数用于岩性识别，那么对反演剖面进行岩性识别就会陷入困境.此外，在交会图进行区域划分方面，一般都采取粗略地描述，或者手动勾绘，这种方法存在很大的不确定性及人为因素［3］.但是相对于直接通过地震剖面解释的结果判别岩性，这种方法有一定的数据分析作为基础，无疑要比仅定性分析的地震剖面识别方法更具有可信度.

基于概率统计理论进行岩性识别的方法在国外已趋成熟，这种方法通过统计岩石物理模型，从叠前地震数据反演得到储层物性参数，进而进行岩性、流体等的预测.Mukerji等［4］给出了统计岩石物理方法应用流程，考虑实际数据误差，并给出模型精确度描述.Mukerji等［4］和 Eidsvik等［5］研究了利用统计岩石物理模型预测储层物性参数的方法，并给出了不确定性分析.Bachrach［6］基于 Gassmann方程建立统计岩石物理模型，计算孔隙度和含水饱和度用于岩性、流体预测.Larsen［7］、Gunning和 Glinsky［8］、Buland［9］、Eidsvik［10］以及 Marit等［11－14］使用叠前反演得到的弹性参数，经由岩石物理模型，通过马尔科夫链、马尔科夫随机场的条件概率形式作为岩性流体的垂向先验分布，成功利用该方法进行岩性流体预测.邓继新和王尚旭［15］实现了基于统计岩石物理模型计算孔隙度和含水饱和度；胡华锋［16］以地震资料和测井资料为反演基础资料，统计岩石物理模型，反演得到孔隙度、含水饱和度和泥质含量信息，利用这些信息来进行储层的预测.

本文试图直接利用叠前反演得到的弹性参数，通过对测井资料进行统计，运用概率统计的方法对岩性进行分析.这种方法在早期是通过岩性量板或者交会图，建立起岩性同弹性参数之间的关系来进行岩性识别.但是由于弹性参数精度不高，并且在参数范围上存在大量重叠，导致识别准确度不理想，同时，不论岩性量板法还是交会图分析法都需要进行大量的准备工作，识别效率较低.在这种情况下，本文通过统计测井资料，建立起描述不同岩性的似然模型，然后用马尔科夫随机场建立先验模型.在贝叶斯的框架下进行岩性识别，有利于综合各种不同的先验信息.该方法不仅仅限于本文提到的三种弹性参数，可以根据实际情况，采用对岩性敏感的弹性参数或者岩石物理参数约束岩性识别过程.本方法的优势在于，能够通过马尔科夫随机场先验模型，建立起点与点之间的相互关系作为空间约束，通过这种邻域内的相互关系，增强识别的准确性和空间上的连续性.该方法的意义在于，能够克服参数范围存在重叠情况下导致的识别不准问题，这就意味着，即使弹性参数反演的结果存在一定误差，通过邻域系统地调整，也能得到较准的识别结果.简而言之，该方法不仅仅得到一个点的岩性，而且能够描述相邻位置的相互关系和相互作用，得到更符合实际情况的岩性剖面；同时，该方法加入了测井的指导，只要测井中出现的岩性，均能在剖面上识别.所以，通过对多口井的统计，能够全面得到各种岩性的弹性参数分布范围，从而准确得到岩性识别的结果.

2 方法原理

Marit和 Henning［12－13］介绍了利用马尔科夫随机场的条件概率形式作为先验模型进行岩性/流体反演，受此启发，结合图像分割中的一些方法，本文利用反演的弹性参数作为特征场，以马尔科夫模型作为先验的岩性标号场对岩性识别进行了研究.考虑到地下储层边缘几何特征，弹性参数在空间上的分布等因素，用马尔科夫随机场模型描述待识别的地下岩性的先验分布符合沉积规律和地下构造实际情况.马尔科夫随机场方法是建立在马尔科夫随机场模型和Bayesian理论的基础上，马尔科夫随机场模型提供了不确定性描述与先验知识联系的纽带，并利用叠前反演得到的弹性参数，根据统计决策和估计理论中的最优准则确定岩性分类问题的目标函数，求解满足这些条件的最大可能分布，从而将分类问题转化为最优化问题.

马尔科夫随机场模型在储层属性随机模拟中应用已久，它利用变量状态的转移概率矩阵对随机过程的未来状态进行分析预测，是一种半定量的研究方法［17－20］.简而言之，某点的性质，只与其邻域系统有关.基于马尔科夫随机场的地质统计模型可以很好地反映复杂空间的连续性，显然这种结果对岩性识别是有利的.

2.1 马尔科夫随机场

马尔科夫随机场是对一维马尔科夫过程的扩展.对于二维平面上的岩性剖面，可以将其视为一个马尔科夫随机场.

图1 不同阶邻域系统及其势团Fig.1 Different order neighborhood system and its clique

设关于岩性π的随机场是关于邻域系统δ的马尔科夫随机场，ω为关于岩性π的分类标号.马尔科夫随机场的性质可以描述为

式（2）称为马尔科夫随机场的条件概率形式，事实上，它很难求取、Marit和Henning利用从测井资料统计得到的计数矩阵对其进行了描述，反映了岩性剖面的局部性质，但是这种方法无法对整个剖面进行定义.Hammersley－Clifford定理揭示了马尔科夫随机场与Gibbs分布之间的等价性，把马尔科夫随机场与Gibbs分布对应起来，通过单个位置点及其邻域的简单特性得到整个剖面全局特性.

Hammersley－Clifford定理：设S是一个邻域系统，当且仅当p（π）是一个关于S的Gibbs分布时，π是关于S的马尔科夫随机场.Gibbs的分布形式为［21］

其中，Z是归一化常数，T表示温度，U（ω）表示能量函数.

2.2 MRF—MAP框架下的岩性识别

本文岩性识别的方法是在Bayesian框架下，利用马尔科夫随机场对岩性的先验信息进行描述，反演目标为岩性类别，已知参数为叠前反演得到的弹性参数.这里 MAP是指最大后验概率，它是Bayesian框架下反演最常用的最优化准则.

如之前定义，S为定义在岩性剖面上的空间位置集合，假设反演得到的弹性参数场为R，那么R的任意一次实现可以记为r＝ （r1，r2，…，rM×N），岩性反演问题需要求解的是满足最大后验概率准则条件下对于每一个位置点岩性的分类，这个分类过程通过对岩性进行标号完成，将标号结果称为标号场，记为ωk，表示共有k个类别.根据贝叶斯公式［22］，有

其中p（r）为弹性参数场的先验分布，通常情况下是一个常数，不参与计算，一般不予考虑.p（ω）是岩性标号场的先验Gibbs分布，满足马尔科夫性.在本研究中，采用马尔科夫随机场的potts模型［23］，其能量函数具体表达式为

其中，δ（ωi，ωj）为 Kronecker delta，U（ω）为式（3）中的能量函数，通过其能量来描述某点属于某种岩性的概率，最大概率的状态对应于邻域点为同一岩性的状态，通过这种描述，使得在一定范围内，相近的弹性参数被划归为同一种岩性，保证空间上的连续性.p（r｜ω）是似然函数，一般情况下，假定弹性参数之间彼此独立，即

其中，n为参与计算的弹性参数的种类数目.假设p（r｜ω）服从高斯分布，即每个类都可以由弹性参数的期望μi和标准差σi来唯一确定其高斯分布：

给定岩性类别的似然函数与先验函数以后，通过先验函数能够保证空间上的连续性，同时，通过似然函数的调整，能够使每一点的岩性在多种弹性参数共同作用的情况下尽可能准确地分类.

综上所述，在MAP准则下求解分类结果，可以归结为求使式（4）最大时候ω的估计，记为

整个建模过程具体如下：

（1）根据问题的性质，定义S上的邻域系统，选择合适的阶数得到相应的势团；

（2）定义势函数的具体形式并写出Gibbs分布；

（3）定义似然函数，选择合适的分布形式；

（4）写出后验概率最简化的表示形式；

（5）估计模型参数并优化寻找MAP准则下的解.

2.3 算法流程

针对本文研究的岩性反演目标函数，采用条件迭代（Iterated Conditional Mode简写为ICM）的算法求解该反问题.假定通过叠前反演得到的弹性参数场为R，岩性标号场为ω，具体算法流程如下：

（1）通过统计测井资料，得到不同岩性下参数与岩性反演的各弹性参数的均值μi与标准差σi；

（2）根据似然函数最大准则求取最大似然分类，以此作为初始分类0，即对剖面上的每一点取其最大似然概率下的标号，遍历整个剖面得到初始分类；

（4）比较第k次与第k－1次的分类结果，如果满足迭代终止条件则将第k次结果作为最终分类结果，否则取k＝k＋1，转到第3步继续计算，直到满足迭代终止条件为止.

将上述过程表示成流程图的形式如图2所示.

图2 ICM算法流程图Fig.2 ICM algorithm flow chart

3 模型试算

为了测试上述方法，本文设计了两个模型进行检验：（1）水平和楔形地层模型，用以检验本文方法在弹性参数范围有重叠情况下识别岩性的能力和描述薄层的能力；（2）Marmousi II模型，用以检验本文方法描述地下复杂构造空间形态的能力.

3.1 水平和楔形地层模型

图3为一个垂向上为500个样点，横向上为100道的岩性剖面，整个剖面上一共分为4种岩性.岩性为块状分布，在水平层位置，没有横向上弹性参数的变化.在垂直位置101～200的范围内，有一楔形模型，用来验证对薄层形状的刻画能力.

假设在第50道有一口井，并有纵波速度、横波速度和密度曲线.通过自上而下统计各岩性的弹性参数范围，求取高斯分布的均值和标准差.各岩性的弹性参数范围如图4所示：

通过图4可以看到，不同岩性弹性参数范围之间存在一定的重叠，这样，仅仅通过观察弹性参数剖面很难区分出不同的岩性.在实际资料中，也存在这样的情况，所以，通过弹性参数直接识别岩性存在一定的困难.

图5是用本文方法得到的识别结果，从图中可以看到，楔形模型得到了完美的识别，表明本文方法能很好地描述地下的构造形态，并尽可能地保留反演结果的分辨率.这对于刻画砂体展布、河道形态和储层范围具有重要意义.

但是在图5中也可以看到，在垂向坐标值51～55中间有一个薄层未被识别，这是因为在测井资料统计的过程中，是按照点数统计，该岩性的这个薄层所占点数很少，所以对于最终的概率分布函数的形态影响不大，另外，该部分与其上覆的岩性在弹性参数数值上有较大的重叠，所以在识别结果中归为上覆的地层.在实际中，这种情况并不多见，对于一种岩性，通过统计其弹性参数得到的概率分布函数，能够全面均匀地包含各点信息，因此，能够准确地描述各点属于该种岩性的概率.

图5中，除去垂向坐标值51～55中的一层外，每种岩性都得到了准确分类，表明本文提出的岩性识别方法合理有效.当地层弹性参数在某范围内大致均匀分布时，即使与其它岩性的弹性参数有较大的重叠，在多参数联合参与识别以及先验函数共同调整下，这些岩性也能得到很好地识别.

3.2 Marmousi II模型

为了进一步检验该方法对于复杂地下构造形态的识别能力，以下采用Marmousi II模型进行检验.

Marmousi II模型是一个垂向、横向上弹性参数变化比较大的模型.在如图6所示横向上第100道、第700道的位置处，假设有两口井，对100道位置的井根据其弹性参数范围，人为的分为5种岩性，并使用同样的方法流程统计弹性参数范围和计算未知参数.700道位置处的井不参与统计，用于后面同识别结果进行对比.最终得到基于两种不同邻域系统的分类结果（如图6）.图6a为基于二阶邻域系统的分类结果，图6b为基于一阶邻域系统的分类结果.

从图6a中可以发现，整个剖面上分类层次清楚，5种岩性都得到了准确地识别，结果的形态符合Marmousi II模型.每种岩性层内平滑，不存在分类错误的噪声点.垂向上的薄层得到了较好辨别，横向上的结果比较连续，并且图中背斜、断层的形态都得到了准确地刻画，这对于油气的勘探是很有利的.表明本文方法在地下弹性参数连续且有较大重叠的情况下，也能比较准确地获得理想的分类结果.

图6b是一阶邻域系统对应的分类结果，与二阶邻域系统相比较差别不大，主要区别是在空间上的连续性不如二阶.

图7为横向第100道、700道位置处的真实岩性与识别岩性结果的对比.在两个位置，两条曲线都基本重合，在岩性发生剧烈变化的位置处，识别结果有一些偏差，但是这种偏差很小，不影响最终整个剖面上的分类结果.这表明本文岩性识别的方法精度还是很高的.

图8描述该方法的收敛速度，取误差的对数显示.从图中可以看到，本文方法收敛很快，经过10次迭代以后，基本得到收敛，因此在计算效率上是可以接受的.

为了讨论反演处理不确定性对该方法的影响，在Marmousi II模型中，分别加入了1%和5%的随机误差（5%的误差分布如图9所示），以检验处理过程中产生的不确定性对岩性识别的影响.

图8 收敛速度Fig.8 Convergence rate

图10为加入1%和5%的误差后岩性识别的结果.从图中可以看出，随着误差的增大，各种岩性大致上还是比较清楚地区分开，同时，岩层的空间结构也得到较好地刻画.这表明该方法在反演结果不太准确的条件下，也有比较准确的识别结果.但同时，对比图10a和图10b可以发现，随着误差的增大，岩性识别的结果会受到一定影响，一些岩性层变得不连续，产生模糊的边界.这表明反演的不确定性会对岩性识别产生一定程度的影响，并且随着误差比例的增加，对岩性识别结果的影响也会增强.

4 结 论

本文提出了一种利用反演的弹性参数识别岩性的方法.该方法主要的优势在于，对于有一定范围重叠或者误差的弹性参数，本方法能够在最大程度上避免识别错误，且高效快捷.该方法的岩性识别主要通过马尔科夫随机场邻域系统的调整来实现.事实上马尔科夫随机场本身并不是一种岩性识别方法，但它可以作为一个先验模型嵌入到某种特定的算法中，以达到影响岩性识别结果的目的.具体而言，其通过调整相邻点的相互关系和相互作用，从全局角度描述岩性，将各点的相互关系加以传播，从而得到空间上有一定连续性的识别结果.同时，该方法的识别过程不损害反演结果已有的分辨率，能够最有效地保护反演分辨率.

从本文中的例子可以看到，马尔科夫随机场能够很好的描述剖面上弹性参数表示出的轮廓和结构，空间上有一定的连续性，能够细致刻画地下构造空间上的形态，得到合适的分类结果.马尔科夫随机场模型在反映岩性剖面随机性的同时，又能反映剖面的一些潜在固有的空间结构，可以更有效地描述地下构造.本方法采用纵横波速度和密度参数进行岩性识别，事实上，在具体使用中，并不限于这三种弹性参数，采用对岩性更为敏感的参数，能够使本方法的识别结果更精确.此外，该方法也存在一定的局限性，马尔科夫随机场先验模型产生的是一个块的识别结果，如果识别的地层需要强调一种连续的变化，该方法会有一定的局限.

（References）

［1］ Ostrander W J.Plane－wave reflection coefficients for gas sands at nonnormal angles of incidence.Geophysics，1984，49（10）：1637－1648.

［2］ 何又雄，姚姚.基于参考道的岩性识别与岩性剖面非线性反演.石油勘探与开发，2005，32（3）：61－63.He Y X，Yao Y.Lithological profile non－linear inversion by reference trace－based lithological identification.Petroleum Exploration and Development （in Chinese），2005，32（3）：61－63.

［3］ 李国福.多参数储层预测及流体识别方法研究［硕士论文］.成都：成都理工大学，2011.Li G F. Multi－parameter reservoir prediction and fluid identification method research（in Chinese）［Master′s thesis］.Chengdu：Chengdu University of Technology，2011.

［4］ Mukerji T，Avseth T，Mavko G，et al.Statistical rock physics：Combining rock physics，information theory，and geostatistics to reduce uncertainty in seismic reservoir characterization.The Leading Edge，2001，20（3）：313－319.

［5］ Eidsvik J，Avseth P，Omre H，et al.Stochastic reservoir characterization using prestack seismic data.Geophysics，2004，69（4）：978－993.

［6］ Bachrach R.Joint estimation of porosity and saturation using stochastic rock－physics modeling.Geophysics，2006，71（5）：O53－O63.

［7］ Larsen A L，Ulvmoen M，Omre H，et al.Bayesian lithology/fluid prediction and simulation on the basis of a Markov－chain prior model.Geophysics，2006，71（5）：R69－R78.

［8］ Gunning J，Glinsky M E.Detection of reservoir quality using Bayesian seismic inversion.Geophysics，2007，72（3）：R37－R49.

［9］ Buland A，Kolbjørnsen O， Hauge R，et al.Bayesian lithology and fluid prediction from seismic prestack data.Geophysics，2008，73（3）：C13－C21.

［10］ Eidsvik J，Omre H，Mukerji T，et al.Seismic reservoir prediction using Bayesian integration of rock physics and markov random fields：A North Sea example.The Leading Edge，2002，21（3）：290－294.

［11］ Spikes K，Mukerji T，Dvorkin J，et al.Probabilistic seismic inversion based on rock－physics models.Geophysics，2007，72（5）：R87－R97.

［12］ Ulvmoen M，More H.Improved resolution in Bayesian lithology/fluid inversion from prestack seismic data and well observations：Part 1－Methodology.Geophysics，2010，75（2）：R21－R35.

［13］ Ulvmoen M，More H，Buland A.Improved resolution in Bayesian lithology/fluid inversion from prestack seismic data and well observations：Part 2－Real case study.Geophysics，2010，75（2）：B73－B82.

［14］ Ulvmoen M，Hammer H.Bayesian lithology/fluid inversioncomparison of two algorithms.Computational Geosciences，2010，14（2）：357－367.

［15］ 邓继新，王尚旭.基于统计岩石物理的含气储层饱和度与孔隙度联合反演.石油天然气学报，2009，31（1）：48－53.Deng J X，Wang S X.Joint inversion of saturation and porosity in gas reservoirs based on statistical rock physics.Journal of Oil and Gas Technology （in Chinese），2009，31（1）：48－53.

［16］ 胡华锋.基于叠前道集的储层参数反演方法研究［硕士论文］.青岛：中国石油大学（华东），2009.Hu H F.The Study of Petrophysical－Properties Inversion Base on Pre－stack Seismic Data （in Chinese）［Master′s thesis］.Qingdao：China University of Petroleum，2009.

［17］ Elfeki A，Dekking M.A Markov chain model for subsurface characterization：Theory and applications. Mathematical Geology，2001，33（5）：569－589.

［18］ Carle S F，Fogg G E.Modeling spatial variability with one and multidimensional continuous－Lag Markov chains.Mathematical Geology，1997，29（7）：891－918.

［19］ Weissmann G S， Fogg G E. Multi－scale alluvial fan heterogeneity modeled with transition probability geostatistics in a sequence stratigraphic framework.Journal of Hydrology，1999，226（1－2）：48－56.

［20］ Norberg T，Rosén L，Baran A，et al.On modeling discrete geological structures as Markov random fields.Mathematical Geology，2002，34（1）：63－77.

［21］ 李旭超，朱善安.图像分割中的马尔可夫随机场方法综述.中国图像图形学报，2007，12（5）：789－798.Li X C，Zhu S A.A survey of the Markov random field method for image segmentation.Journal of Image and Graphics （in Chinese），2007，12（5）：789－798.

［22］ Yuan S Y，Wang S X，Li G F.Random noise reduction using Bayesian inversion.Journal of Geophysics and Engineering，2012，9（1）：60－68.

［23］ Pérez P.Markov random fields and images.CWI Quarterly，1998，11（4）：413－437.