所谓MCMC岩性反演，主要是利用贝叶斯理论作为反演方法的理论框架，协同地震、测井、岩性等多尺度地质信息定量描述反演问题的不确定性(Kaipio and Somersalo，2004；Stuart，2010；Hammer and Tjelmel and ，2011，田玉昆等，2013)，通过MCMC算法从参数后验分布中获取样本，综合岩石和储层参数信息基于马尔科夫过程进行地震振幅反演，该方法的优点是能够使反演结果更忠实地震数据，且相比其它波阻抗反演方法能够实现更高精度的岩性反演结果.

在20世纪九十年代，利用地震波属性参数进行岩性反演开始引起国际学者的重视，(Mavko et al.，1998)，由此产生利用地质统计学方法实现地震波振幅变化的岩性参数反演，这种随机反演方法的不确定性主要来自两方面：地震数据与储层参数相关性的不确定性和储层参数与岩性参数相关性的不确定性(张爱印等，2004；Gunning and Glinsky，2004).近年来，国内外学者趋于论述在贝叶斯框架下利用MCMC方法估算储层参数(如波阻抗、岩石密度等连续型参数)(

杨培杰和印兴耀，2008；Andrieu et al.，2010； 张广智等，2011；Beskos et al.，2013；张繁昌等，2014)，这是一种基于地震波振幅信息的参数反演方法，其优点是能够使反演结果更加忠实于作为似然约束条件的地震数据体，再利用MCMC算法从参数后验分布中获取样本(Liu，2001)，为增快计算速度，也可以利用改进的算法从近似后验分布中取样(Larsen et al.，2006；Hammer and Tjelmel and ，2011)，但这种近似算法易将一些后验分布稀疏化进而忽略掉其中的重要信息，这推动了从非近似后验分布中模拟取样的方法研究(Bul and et al.，2008).

利用MCMC岩性反演方法在计算储层参数数据体的同时能够随机生成相应的岩性数据体，在实际运用中，人们往往更关心这些更为直观的岩性数据体的计算原理以及其实现过程的不确定性(张广智等，2014).在岩性反演算法中通常假设岩性参数具有固定的空间相关性，用以反映储层中岩性的分布特征(张宏兵等，2005；杨锴等，2012).岩性反演惯用的手段是采用贝叶斯方法整合不同类型的多尺度信息(如地震、测井、岩性等)，来求取地震弹性参数和岩性参数的后验分布(Scales and Tenorio，2001；陈建江和印兴耀，2007)，而对于一些“反问题”，考虑到计算空间的有限容量，需结合MCMC法(或其它迭代抽样算法)从后验分布空间进行随机采样(魏超等，2008；万玲等，2013；王保丽等，2015).目前岩性反演的主要问题依然是各参数之间相互关联的不确定性(Bosch et al.，2007，2009).

本文分‘贝叶斯框架的构建原理’和‘反演算法原理’两部分系统阐述在贝叶斯框架下基于MCMC算法预测储层砂体的地质统计学原理，将反演过程运用的统计学符号赋予地质含义，进而分析MCMC法岩性反演的主要实现过程，具体包括建立基于马尔科夫过程的岩性参数先验函数(Krumbein and Dacey，1969)、通过岩石物理分析建立岩性和储层参数的联系(Avseth et al.，2005)、建立似然函数反映储层参数与地震数据的联系(Bul and and Omre，2003)、构建能反映空间相关性的马尔科夫链，实现基于蒙特卡罗迭代算法在复杂后验分布空间的公平采样四个方面.同时应用该算法对朝阳沟油田扶余油层未开发区朝65井所在岩性剖面进行反演预测，基于岩性参数的空间结构与储层参数的空间联系，以地震波振幅信息作为似然条件进行岩性参数反演，利用多次反演实现结果验证MCMC法的可靠性和不确定性.

对于油田的未开发区块，已知信息具有不确定性和较大局限性，需考虑利用多尺度信息综合预测储层砂体，以降低钻井风险.为协同已知地质信息，通常以贝叶斯法作为理论框架(Scales and Tenorio，2001; Duijndam，2003；苑书金，2007；吴媚等，2008；杨锴等，2012).贝叶斯法主要根据最小风险代价或最大似然比进行判决，是依据贝叶斯准则进行判别分析的一种多元统计分析法(陈建江和印兴耀，2012；杨培杰，2014).

针对地下空间的某一样本点的储层砂体预测，贝叶斯判别公式可以衍生为在已知地质信息的约束下，计算该样本点的储层参数和岩性参数的概率密度函数.图 1表示MCMC反演方法的贝叶斯框架，其中X表示储层或岩性参数，E(Evidence)表示验证性地质信息(真实地质数据，包括地震、测井等真实信息)，H(Hypothesis)表示假设性地质信息(经地质参数分析后得到预测性信息)，贝叶斯判别公式表示为

图 1

Fig. 1

图 1 贝叶斯框架Fig. 1 The Bayesian framework

其中P(H|E)表示符合验证性地质信息的同时，能够符合假设性地质信息的概率分布，可近似看做是全概率事件(P(H|E)~1)；P(X|H)表示经地质参数分析(如已知样本点概率密度函数分析和变差函数分析等)后得到的预测性储层参数的先验函数；P(E|X)表示在先验函数取值空间内，符合实际地震资料、测井资料的似然概率密度函数；P(X|H，E)表示在先验和似然函数的约束下，得到的储层和岩性参数的后验概率密度函数.

建立式中先验函数需要区分岩性对储层参数进行岩石物理分析，通常认为不同岩性内的储层参数的概率分布具有相互独立性，而同一岩体内的储层参数会有一定相关性(Kj?nsberg et al.，2010)，这是由于同一地层内相同岩性的成岩条件具有一定相似性，建立先验函数的难点是将岩性参数同储层参数联系起来，两种参数相匹配的不确定性是MCMC反演结果具有不确定性的主要原因(Bosch et al.，2007，2009)；似然函数通常是针对地震信息的残差分析，但如果要使反演实现完全符合测井数据，可将测井解释结果或者波阻抗曲线等地质信息作为似然函数进行约束，通常测井信息更多用于构建先验函数或作为反演成果的验证条件，而非一种硬性约束信息(杨立强，2003).

在实践中，由于未知参数多为高维、复杂的后验概率分布，且高维积分计算十分困难，使得贝叶斯推断方法在实践中的应用受到很大的限制.MCMC方法突破了这一计算问题，它能够建立反映岩性参数及储层参数空间相关性的马尔科夫链，并通过转移核规则引导马尔科夫链扰动过程，在多约束条件下得到与地震数据匹配的全局最优解.

本文以一条垂向线上的样本点为例，通过一维视角阐述该反演算法.在实现过程中，如图 2所示，将地质模型视为多样本点空间，其中垂向线上样本点被赋予岩性参数和储层参数(岩性代码和纵波阻抗值)，样本点编码为i∈(1，…，n)，垂向线上各个样本点MCMC反演的核心是通过构建一条马尔科夫链拟合岩性参数和储层参数之间的垂向相关性，下面着重介绍马尔科夫链的构建原理.

图 2

Fig. 2

图 2 基于马尔科夫过程构建先验函数Fig. 2 The prior construction based on the Markov process

将岩性概括为砂、泥岩两类，每个样本点均被赋值岩性参数代码，样本点的厚度(采样间隔)决定岩性参数的分辨率.垂向线上的岩性参数特征表示为各岩性出现的相对频率、垂向上出现相同岩性的平均间隔厚度及各岩性排列易出现的次序.

构建马尔科夫链，若要定义能够保留岩性分类特征的参数，需要假设垂线上各点岩性仅取决于临近点的岩性，如图 2所示，垂线上岩性参数被定义为

f=(f₁，…，f_n)，

其中样本点i的岩性被定义为f_i；0…L分别对应不同岩性，对于任意样本点，均有其相应的岩性相对概率.利用岩性参数的概率分布空间，垂向上按由上至下(或者由下至上)的顺序定义马尔科夫链，对于任意样本点f_i，在f_i-1的岩性概率已知的条件下，该样本点的岩性概率只依赖f_i-1的岩性概率分布，这是马尔科夫链的基本性质.基于以上描述，构建马尔科夫链的过程(图 2)可被表述为

公式中P(f₁)表示顶部样本点岩性概率的分布，P(f_i|f_i-1)是已知前一样本点f_i-1岩性概率分布的前提下，当前样本点f_i的岩性分布，也就是马尔科夫链的转移矩阵(亦称转移核).该转移核受分辨率影响，且分辨率取决于样本点厚度，若样本点厚度改变，转移核将随之变化，然而转移核虽有变化，却能大体保留主要岩性参数特征.

通过以上方法可以构建出一个简单的转移核(见图 2)，该转换核只是用来定义垂线上的岩性参数的先验分布.如果先验分布受地震波的振幅信息的约束，马尔科夫链就不再是简单的一步转换，我们将储层参数(本文特指纵波阻抗值)定义为m=(m₁，…，m_n)，样本点的岩性参数和储层参数的关系通过岩石物理分析得到：在模型中每一个样本点的储层参数都会被定义为：m_i，m_j，m_k(i、j、k代表三维空间的三个方向)，由于本文分析的是一条垂向线，可以只考虑垂向方向.将P(m_i|f_i)表示样本点i、岩性为f时对应的储层参数概率分布，它可以通过模拟该样本点的高斯分布得到，高斯分布由期望E{m_i|f_i} = μ_fi和方差σ_fi定义为

由此可见，任意岩性参数都对应不同的储层参数概率分布，这导致储层参数的分布具有多模态性质.

储层参数的空间相关性并非由上述参数分析得到，如图 2右侧所示：当我们考虑两个样本点i和j，如果它们来自同一块砂体，那么这两个样本点的储层参数会有一定的相关性.这种特性可以通过空间相关性来模拟，先假设储层参数在不同岩性区间内的储层参数是不相关的，而在相同岩性区间内是相关的，如此假设的结果是同一块砂体或同一块泥岩区域内不同样本点的储层参数概率密度分布函数会具有一定相关性，如果在i和j之间所有样本点都具相同岩性，我们可以得出其协方差：

公式中σ_f是样本点所在位置对应的储层参数的方差，(t是样本点厚度(采样间隔)，R是通过已知参数的样本点进行变差分析得到的垂向变程，垂向变程决定了样本点参数在垂向相关性的消减速度.R值偏小会导致同一岩性区间内相关性降低，R值偏大会导致储层参数空间相关性增大.

据此可以推断：在样本点岩性确定的情况下，P(m|f)同样具有马尔科夫性质，即在岩性参数作为已知条件时，当前样本点m_i的储层参数的分布与前一样本点m_i-1的储层参数分布具有一定相关性，由此构建其马尔科夫链(图 2)：

将具有空间相关性的P(f)(公式2)和P(m|f)(公式3)马尔科夫链运用到贝叶斯框架中，即可成为先验函数P(X|H)的组成要素(公式6)，此函数可作为假设性先验函数运用到样本点的蒙特卡洛迭代中.公式(6)为

MCMC反演方法中似然函数的主要作用是将储层参数与地震振幅信息联系起来，它通过将储层参数(本文主要针对纵波阻抗)转换为反射系数，再与子波褶积形成能够与地震波比较的合成记录.通常叠后反演研究的是水平叠加和偏移剖面，在界面倾角不大时，可以近似认为是垂直入射.将地震波的振幅信息定义为d=(d₁，…，d_n)，代表垂线上各样本点的地震信息、m代表储层参数、e代表观察噪音、g(.)代表正演算子，该问题可以表述如下：d=g(m)+e.在实际情况中，观测噪音通常为白噪声，似然函数一般采用正态分布.假设噪音是均值为0，方差为(_e的高斯随机噪音，n表示观测数据的个数，可以得到观测地震记录的似然函数为

上述似然函数表示的是观测数据与模型参数的拟合程度，似然函数的值的大小与二者拟合程度具有正相关(即似然函数越大，观测数据与模型参数拟合得越好，反之亦然)，由此可知地震误差(即地震资料的不确定性)主要来源于储层参数与地震资料相匹配的不确定性.

由公式(1)可推断后验函数(公式(8))，其中P(f)P(m|f)为先验函数，P((d|m))为似然函数，该公式主要是推断岩性参数和储层参数的后验分布，它是一个较为复杂的多维函数，甚至难于用具体公式来表达，MCMC法所解决的就是如何从复杂的多维分布空间取样，针对样本点建立迭代条件，抽取能够收敛于后验分布函数的样本.公式(8)为

通过MCMC算法可以从后验函数生成大量f，m(岩性和储层参数)的实现，这些实现代表后验函数的概率分布，通过从先验分布P(f)P(m|f)中随机抽取岩性参数和储层参数的初始值进行首次迭代，进而每一次迭代过程的建议分布均由前一次迭代的实现得出.其中关键是通过Metropolis-Hastings算法构造转移核，建立一个马尔科夫链使其具有与目标分布(后验函数)一致的平稳分布，Metropolis-Hastings算法是一种迭代算法，利用它能够在马尔科夫过程中产生新的岩性和储层参数(m^*，f^*)的建议分布，并使之与前一步的岩性和储层参数(m，f)的概率分布相比较，如果建议被接受，(m^*，f^*)就成为当前分布，否则回归之前的(m，f)再重新进行迭代.建议分布是否被接受由一个随机的接受概率来决定，公式为

其中，P(f^*，m^*|d)和P(f，m|d)是在公式(8)中分别由当前建议分布和前一步实现的后验函数的计算值；q(m，f|m^*，f^*)是由前一步的实现产生新的建议分布的概率判定准则，可称为正向转移核； q(m^*，f^* | m，f)是反向转移核，是由当前建议分布产生前一步实现的概率判定准则，通过构建接受概率 P_接受，使马尔科夫链的极限概率分布收敛于后验概率分布P(f，m|d)(Green，1995).

通常采用混合高斯分布作为建议机制，根据后验函数迭代产生新的实现，该过程主要包含三个步骤：

(1)在前一步岩性参数的概率分布中随机抽取建议岩性参数为(Geyer and M?ller，1994)

(2)在岩性参数和地震数据的约束下随机抽取建议储层参数为

(3)通过接受概率判定是否接受建议分布，其公式为

其中P(f^*，m^*|d)和P(f，m|d)是来自公式(9)的后验分布；q(f^* | f)和q(f | f^*)代表岩性参数的正向和逆向转移核；q(m|f，d)和q(m^*|f^*，d)代表储层参数在地震数据和当前岩性条件约束下的正向和逆向转移核，如果P_接受判定接受，则将当前状态设为(f^*，m^*)，否则维持前一状态(f，m).为简化计算，通常假设马尔科夫链是可逆的，即正、反向转移核一致，将公式(3)、(7)、(8)分别代入公式(12)可得到公式(13)，经过多次迭代，马尔科夫链的极限收敛于后验分布，能够得到样本点的岩性参数和储层参数的平稳实现.公式(13)为

以朝阳沟油田朝84-6区块未开发区朝65井为例，图 3上部显示朝65井所在地震剖面，反演时窗选择在1.1~1.24 s之间，地震资料的重采样间隔为1 ms，图 3下部为朝65井的井震标定成果，从左至右分别为：时间道-地震道-合成记录道-岩性参数道-储层参数道(波阻抗曲线)-深度道，经过朝84-6区块的岩石组构分析，岩性主要以浅水枝状三角洲分流河道相的长石质岩屑砂岩和分流河道间相的泥岩为主，合成记录是通过褶积公式计算(该公式与MCMC反演中所用到的褶积公式一致)，比较原始地震数据与合成地震数据，可以看出朝65井井震匹配较好.图中地震解释反射层T₀、T₁、T₂、T₃分别对应地质分层：FI₁₁(扶I砂组顶部)、FII₁(扶II砂组顶部)、FIII₁₁(扶III砂组顶部)、FD(扶III砂组底部).

图 3

Fig. 3

图 3 朝65井的井震标定示意图Fig. 3 The well-seismic calibration of the Chao65 well

由地震资料预测岩性分布的不确定主要体现在两方面：一是储层参数与地震数据不匹配，二是岩性参数与储层参数不匹配.图 3中可见朝65井的井震匹配较好，即储层参数与地震数据匹配较好.在此基础上，对储层砂体反演结果进行分析，通过分析不同模型参数(采样间隔、参数概率分布等)对反演结果的影响，来验证MCMC反演预测储层砂体的可靠性.

朝65井的井信息可以给出贝叶斯框架的先验分布，图 4显示针对不同采样间隔，在先验函数中岩性参数的转移矩阵(即fi的转移核：P(fi│fi-1)(见公式2))和岩性参数的先验分布(图 4中的固定概率向量)，以及在各自先验条件下得到的朝65井所在剖面的砂岩概率实现(转移矩阵和固定概率向量均可由软件统计计算得到，反演所用地震资料分别为进行0.5 ms和2 ms重采样，以保证地震和测井的采样间隔一致).图中不同实现结果对比可知，选取不同采样间隔，岩性概率转移核会有所不同，但其参数分析结果及岩性概率反演实现结果依然能大体保留主要岩性分类特征.由于研究区储层致密，总体呈砂泥薄互层的特征，选取0.5 ms作为采样间隔，使先验函数能更精确描述先验信息的概率分布.

图 4

Fig. 4

图 4 不同采样间隔反演实现结果对比图Fig. 4 The comparison of inversion realizations with different measuring interval

图 5显示以岩性参数和储层参数的概率分布作为先验条件，分扶I砂组、扶II砂组、扶III砂组分别计算的3个MCMC反演实现结果，图下部是朝65井所在反演剖面，井柱显示为实钻砂体(黄色为砂岩，蓝色为泥岩)，将其叠置于反演剖面用于检验反演成果.观测3个实现结果能够发现：对于储层参数，观察各个反演剖面左侧的储层参数曲线(蓝色为测井数据，红色为反演结果)可知其符合率较高，这是由于井震匹配较好(见图 3)，储层参数与地震数据的不确定性相对较低.

图 5

Fig. 5

图 5 MCMC法岩性反演结果Fig. 5 The post-stack MCMC lithology inversion realizations

对于岩性参数，扶I砂组、扶II砂组反演出的砂体与实钻砂体较为匹配，部分砂岩厚度与实钻砂厚略有差异，这一方面说明地震资料与储层参数匹配较好，能够在地震数据作为似然条件进行约束时，生成与实钻砂体较为匹配的实现结果；另一方面说明井上储层参数的先验分布被砂岩和泥岩区分较好(见图 5左上部，这里储层参数指由声波时差和地层密度曲线合成的纵波阻抗曲线，储层参数先验分布的横坐标是经过对数转换后的相对数值，数值越大对应波阻抗越大)；观察扶III砂组的储层参数先验分布，明显砂岩和泥岩的高斯概形叠置较多，岩性参数与储层参数的不匹配，导致反演结果的不确定性增大(图 5中红框内反演砂体明显与实钻砂体不一致，并且三个实现也各有差异).对于两者不匹配的原因，可能是由于扶III砂组接近钻井底部深度，测井数据易出现不稳定，加之其砂地比较低，岩石物理分析不易出现规律性，使不同岩性的储层参数概率分布有部分叠置情况.

观察反演结果的岩性参数分布(可看做是MCMC反演岩性参数的后验分布，对应图 5的右上部)，对比各砂组岩性参数后验分布与先验分布无明显规律性(扶II砂组先后较为一致，扶I砂组和扶III砂组变化较大)，这主要是由于MCMC反演方法具有全局寻优的特点，由于先验函数只是根据已知数据推测出的假设性地质信息(见图 1)，而在地震数据作为似然约束条件时，经过多次迭代后，反演结果会跳出先验函数的限制，更忠实于构建验证性似然函数的地震数据.

观察三个实现结果各砂组对应的地震剖面(图 5中地震剖面的绿色、红色和蓝色线是对照反演实现结果在各砂组进行的单砂体刻画)，可以发现扶I砂组的砂岩大多对用地震反射波偏下的部位，地震波自下而上振幅增强，地震波横向连续性较好，导致反演预测砂体具有薄层大面积连续分布的特点，反演结果符合前人对该地区的沉积研究成果，该沉积时期朝84-6井区所在区域平稳沉降，地形起伏较小，湖水快速水进，由于河湖共同作用导致砂体呈薄层席状分布，垂向上砂泥间互，侧向上砂厚变化幅度小，连续性好；扶II砂组的地震剖面的地震波侧向连续性较差，振幅强度中等，反演砂体主要对应正极性地震波的中下部，砂体上倾尖灭和砂体叠加现象明显，这是由于该沉积时期属于三角洲平原沉积，以枝状河流亚相沉积为主，由于河流摆动、河道迁移等作用导致河道砂体频繁切叠、废弃而出现三角洲平原分流河道的沉积特征；扶III砂组的测井数据和实钻砂体数据吻合并不好，在与地震数据同时参与MCMC岩性反演时，反演结果更符合地震数据(图 5中反演剖面的红圈处对应扶III砂组蓝色虚线刻画的单砂体)，虽然岩性数据并未在红框处有砂体显示，但反演结果依旧根据地震波属性特征反演出砂体，只是砂体的厚度和连续性具有更大的不确定性，这是由于储层参数与岩性参数不匹配以及地震数据与岩性参数不匹配的双重不确定性所决定.由此可见，MCMC反演是更忠实于地震数据的岩性反演，对于扶III砂组这种不确定性较高的反演结果，通常需要进一步检查测井、地质和地震资料的合理性和准确性，以减小反演结果的不确定性.

针对反演结果的不确定性，通常将多次实现的结果进行计算，公式为

其中P(f_i=l|d)是地震约束下样本点i=1，…，n及岩性l=1，…，L的岩性概率，N为实现个数，f_i，k是样本点i处的岩性实现k，I(f_i，k=1)是指示方程：如果f_i，k=1则I(f_i，k=1)=1，否则I(f_i，k=1)=0.

将由公式(14)计算得到的岩性概率作为最终的反演结果.图 5中反演实现结果右侧蓝框内从左至右显示地震道和岩性概率曲线.结合地震剖面和反演实现结果，可进一步验证地震数据作为约束条件去构建似然函数的重要性，同时观察反演结果的符合率，岩性概率曲线和实钻岩性参数道对应较好，这说明MCMC方法能够实现忠实于地震数据的岩性反演.

本文从贝叶斯理论出发，利用假设性地质信息作为先验条件、验证性地质信息作为似然条件，构建能够反映空间相关性的马尔科夫链，实现在多条件约束下的储层和岩性参数反演预测，并保证反演结果能够忠实于地震振幅信息.朝65井的实例分析展示了不同采样间隔对模型参数和反演结果的影响，同时通过对比分析朝65井的多个反演实现结果，分析了MCMC反演方法的不确定性因素，证实了MCMC岩性反演方法能够协同地质、地震、测井等多尺度信息，并使反演结果更忠实于地震数据.

致 谢    感谢各位审稿专家给出的宝贵修改意见，感谢东北石油大学地球科学学院赵荣教授对本文的指导，感谢编辑部主编的热心帮助.