地震波形约束的蒙特卡洛—马尔科夫链随机反演方法

周爽爽 印兴耀* 裴 松 杨亚明

(①中国石油大学(华东)地球科学与技术学院，山东青岛 266580；②海洋国家实验室海洋矿产资源评价与探测技术功能实验室，山东青岛 266071)

0 引言

地震反演是获取地下介质弹性参数，揭示地下储层展布规律、物性及含油气性的有效途径[1-4]。随着油田勘探、开发的不断推进，对储层预测的精度要求也越来越高，因此高精度储层反演技术在薄互层识别中的作用越来越重要[5-6]。地震反演分为确定性反演与随机反演两大类，其中确定性反演以正则化方法为基础[7-10]，通过弱化反演过程的不适定性，求解一个逼近真实解的最优解或近似解，但是该方法受地震分辨率及带限子波的限制，在识别调谐尺度内薄层时具有一定的局限性。

与传统的确定性反演方法相比，随机反演识别薄层具有一定优势[11]。以地质统计学为基础、以测井资料为条件数据的地震随机反演方法的分辨率高于常规确定性反演[12-15]，因此迅速得到广泛应用，但是提高计算效率以及消除随机性一直是难点。基于序贯高斯、模拟退火等地质统计学算法的随机反演利用垂向密集的井点数据求取空间变差函数进行随机建模[16]，然而在求取各个方向的变差函数变程时人为因素影响较大，导致反演结果随机性较强[17-19]。王保丽等[20]构建了基于FFT-MA(fast Fourier transform-moving average)谱模拟的快速随机反演方法，采用逐步变形算法快速迭代求解重新构造的目标函数，将复杂的求解问题简单化，有效提高了运算效率。针对随机采样，以蒙特卡洛—马尔科夫链(Monte Carlo Markov chain，MCMC)算法为例，在获取模型参数后验概率密度分布随机样本的基础上，采用统计分析手段获取模型参数的统计特征(如均值、协方差等)[21-23]。Hong等[24]采用多尺度标准遗传算法的MCMC方法不断优化已知模型以获得最优参数，提高了MCMC算法的计算效率以及估值精度；基于Metropolis算法[25]，Hastings[26]提出了Metropolis-Hastings(M-H)算法，介绍了蒙特卡洛估计误差的相关理论、应用技术以及误差评估方法；Smith等[27]利用MCMC方法获取模型参数的后验分布，尝试利用MCMC算法逼近贝叶斯统计分析中的积分等运算，可提供大量与未知参数有关的全局信息；张广智等[28-32]提出构建一种自适应的建议分布，能够自适应并不断更新构成的马尔科夫链，有效地加快了马尔科夫链的收敛速度；李坤等[33]在待反演参数服从混合概率先验模型的前提下，提出了基于差分进化MCMC随机模型的相约束叠前地震概率化反演方法，充分模拟了待反演模型的后验概率密度分布(如均值、方差及置信区间)。

针对上述方法存在的随机性强、运算效率低等问题，本文提出了基于地震波形约束的MCMC随机反演方法。充分利用观测地震数据和待反演参数之间的地球物理映射关系，应用相关系数，根据已知的地震波形之间的相似性特征指导井数据进行伪普通克里金插值模拟，建立具有地震波形指示的初始模型；在此基础上，进一步在贝叶斯框架下构建观测地震数据和测井数据协同约束的后验概率密度分布，结合M-H采样算法多次随机模拟具有地震波形指示的初始模型参数，利用后验均值作为模型参数的最优解。该方法有效地提高了反演稳定性和横向连续性，降低了随机性，有效地弱化了地震噪声对反演结果的影响，并且极大地加快了马尔科夫链的收敛速度，有效地提高了该方法的运算效率和估算精度。

1 理论与方法

1.1 初始模型

地震振幅的变化间接地反映了地下储层参数的空间变化，由于观测地震数据和储层参数之间存在一定的地球物理映射关系[33]，因此可以根据已知观测地震数据指导已知井信息完成地下未知储层参数的插值运算。基于地震波形约束的MCMC随机反演方法借鉴传统的地质统计学建模思想，建立具有地震波形指示的初始模型的基本思想是：在工区所有井数据的基础上，统计所有井的井旁地震道以及井上与井旁地震道对应位置的预测参数曲线，依据统计得到的待预测点处的地震波形与井旁地震道所有样本点处的地震波形之间的相关系数，对所有井进行相关度排序，优选与待预测点关联度高的井样本建立初始模型；利用普通伪克里金插值方法对含有高频成分的井数据进行无偏最优估计，最终反演的地震波形与原始地震特征一致，其中包含确定性较强的宽频(高频成分)弹性参数。

利用相关系数根据地震波形相似性优选井样本(图1)。根据观测地震数据和待反演参数之间的地球物理映射关系，充分利用地震波形相似性由优选的井样本模拟地下未知模型参数。用普通伪克里金插值公式表示待反演参数的建模过程

(1)

式中：Z(x0)为待模拟插值参数x0处的值；Z(xi)为由井旁地震数据筛选的已预测井样本点xi处的对应值；λi为Z(xi)在待插值处的权重；n为参与井的个数。本文根据普通克里金的定义，推导了一种新的利用地震波形相似性指导井数据进行无偏最优估计的插值方法，并建立了具有波形指示的初始模型。

在普通克里金定义中，假设对于某区域内任意点p(x，y)的期望c与方差σ2都是相同的，即对于任意一点p(x，y)，都有

图1 地震波形和井参数优选示意图[5]①为优选样本，②为统计样本E[p(x，y)]=E(p)=c

(2)

Δ[p(x，y)]=σ2

(3)

无偏估计条件为

(4)

普通克里金插值公式为

(5)

(6)

化简可得

(7)

相应地，估计方差为

(8)

令Cij=Cov(zi，zj)，则有

(9)

相关系数为

(10)

式中rij为地震数据点zi地震波形与地震数据点zj地震波形的相关系数。将式(10)代入式(9)，则有

(11)

我们的目标是寻找使J最小的一组λi，且J是λi的函数。因此将J对λi求偏导数并令其为0，即

(12)

(13)

(14)

由于rkj=rjk，则式(14)变为

(15)

由式(15)得到了求解权重系数λj的方程组

(16)

写成矩阵形式

对矩阵进行逆运算即可求出已知点处地震数据在待预测点处所占的权重，根据观测地震数据和待反演参数之间的地球物理映射关系，以此矩阵求得的权重λj代替式(1)中的λi，从而完成初始模型的建立(图2)。

图2 初始模型建立流程

从初始模型建立流程(图2)可知，应用地震波形相似性指导已知井参数信息建立初始模型的过程与常规序贯采样方式不同，在每一个待预测点都经过严格的波形优选进行伪克里金插值，已插值点并不参与下一个未知点的插值过程，每个插值点的模拟都从实际优选样本出发进行无偏最优估计，有效降低了地震噪声的影响和随机误差的累积。

建立具有波形指示的初始模型过程中需要注意地震波形样本选取的时窗和样本井的最高截止频率这两个重要参数的选取。地震波形样本选取的时窗用于对比待预测点地震波形和井旁地震道波形组成点的个数，下文将该参数简称为有效样本数。参与井样本的最高截止频率代表参与井的目标层井数据达到一定相似性时的频率；当参与井目标层的井旁地震道地震波形达到一定的相似度时，通过不断对参与井的目标层纵波阻抗数据进行高频滤波，发现参与井的目标层井数据之间的相关系数随着最高截止频率的降低而增大，当其与井旁地震道的相关系数相同时，频率范围远远大于地震数据的有效频宽[34]。最后，利用普通伪克里金插值方法在地震波形相似性的指导下对含有高频成分的优选井样本数据进行无偏最优估计，并将优选井的纵波阻抗数据作为先验信息，从而建立具有地震波形指示的初始模型，该模型包含了确定性较强(大于地震有效频宽)的宽频纵波阻抗参数。

1.2 M-H算法

MCMC算法的主要目的是构成一条平稳分布的马尔科夫链，同时该链是不可约的[21]。M-H算法是构造马尔科夫链的一种常用方法。本文采用M-H算法根据某个建议分布不断迭代后验概率密度函数，使其达到平稳条件。

首先，设q(θ，v)是与目标后验概率密度分布π(θ|D)较接近的某已知分布D，一般称其为建议分布，也称为备选生成密度，q(θ，v)表示参数θ的潜在转移θ→v，满足

(17)

(1)任意选择一个初值θ(t)，并令t=0 ；

(2)从q[θ(t)，v]中生成一个备选值θ*，并从均匀分布U(0，1)中抽取随机数μ；

(3)若μ≤α[θ(t)，θ*]，则θ(t+1)=θ*，否则令θ(t+1)=θ(t)，此时接受概率[26](转移核概率)为

(18)

(4)令t=t+1 ，并返回到步骤(2)，反复迭代直至平稳状态。

很容易证明，以此种方式构建的接受概率α(θ，v)满足细致平衡条件。

1.3 基于地震波形约束的MCMC随机反演方法

基于地震波形约束的MCMC随机反演方法，充分利用地震数据之间的相关系数代替测井数据之间的相关系数，构建具有地震波形指示的初始模型。以此为基础，在贝叶斯框架下构建观测地震数据和测井数据协同约束的后验概率密度分布，并将其作为目标函数，结合M-H算法不断扰动模型参数，利用后验均值作为模型参数的最优解，并且取多次有效实现的均值作为期望值输出，有效提高了反演稳定性和横向连续性。

反演问题[35]可以表述为

p(m|d)=f(m)+e

(19)

式中：p(m|d)为待反演参数的后验概率密度函数，m为地震波形指示初始模型约束，d为观测数据；f为正演算子；e为观测噪声。在贝叶斯框架下联合似然函数p(d|m)和先验信息p(m)的p(m|d)可写为[36]

p(m|d)=p(m)·p(d|m)

(20)

(21)

(22)

式中N为采样点的个数。式(20)变为

(23)

式中：ε为可调因子；λ为常数。

为了得到所估参数的后验概率密度分布p(m|d)，利用M-H算法生成马尔科夫链，其建议分布函数取为高斯分布

(24)

其中q(mk，m*)表示采样点为k时的初始模型的值mk的潜在转移(mk→m*)，转移概率[22]为

(25)

为了使马尔科夫链最终收敛于未知参数m的后验概率分布，后验概率密度分布即为平稳分布p(m)·p(d|m)。令

(26)

则转移概率[28]为

α(mk，m*)=exp[α′(mk，m*)]

(27)

其中

α′(mk，m*)=min{0，J(m*)+lg[q(m*，mk)]

-J(mk)-lg[q(mk，m*)]}

(28)

最后，在M-H采样算法下不断扰动模型参数，多次随机模拟具有地震波形指示的初始模型参数，利用后验均值作为模型参数的最优解。

2 二维模型试算

为了进一步验证基于地震波形约束的MCMC随机反演效果，抽取部分Marmousi2模型反演试算叠后纵波阻抗IP。模型时窗为1000～1240ms，从部分Marmousi2模型(图3a)中随机抽取11道纵波阻抗数据作为已知伪井数据(图3b)，利用伪井数据与30Hz的雷克子波合成的地震数据(图3c)建立初始模型(图4)。可见：有效样本数越小，初始模型的横向连续性越差，纵向分辨率越高，反之亦然；有效样本数为10时初始模型的纵向分辨率和横向分辨率相对较高(图4b)，间接反映了原始阻抗模型(图3a)的大致特征，因此取该模型作为模型测试的最终初始模型。

为了更好地测试基于地震波形约束的MCMC随机反演的抗噪性，分别对含噪和无噪地震数据进行二维反演测试。图5为初始模型第42道数据的基于地震波形约束的MCMC随机反演结果。由图可见，地震噪声对反演结果影响不大。图6为反演结果及其合成记录。由图可见：无噪数据、含噪数据地震波形约束的MCMC随机反演结果(图6c、图6e)与原始模型(图6a)相近，表明地震波形约束的MCMC随机反演方法具有较好的抗噪性和可行性；无噪数据、含噪数据的地震波形约束的MCMC随机反演结果合成记录(图6d、图6f)均与原始地震记录(图3c)基本一致。

图3 纵波阻抗模型以及原始地震记录(a)纵波阻抗； (b)纵波阻抗伪井数据； (c)原始地震记录

图4 有效样本数分别为5(a)、10(b)、15(c)和20(d)时建立的初始模型

图5 初始模型第42道数据的基于地震波形约束的MCMC随机反演结果(a)无噪； (b)信噪比为10

图6 反演结果及其合成地震记录(a)原始模型； (b)地震记录(信噪比为10)； (c)地震波形约束的MCMC随机反演结果(无噪)； (d)图c的 合成记录； (e)地震波形约束的MCMC随机反演结果(信噪比为10)； (f)图e的合成记录

3 实际资料应用

将基于地震波形约束的MCMC随机反演方法应用于中国东部F区实际数据。该区面积约为20km2，区内有11口井，反演时窗为2000～2600ms，采样率为2ms。将经过标准化和精细井震标定的11口井作为样本井，建立具有波形指示的初始模型。任意选取一个过两口井的剖面为研究资料。根据样本井纵波阻抗曲线与对应的井旁地震道分析结果，为了尽可能使井曲线与地震波形具有较好的映射关系，以保证目标层样本井曲线间的相关系数大于目标层井旁地震道间的相关指数，对样本井作低通滤波，选取最高截止频率为130Hz，从而建立地震波形指示的初始模型。基于地震波形约束的MCMC随机反演的目的是提高随机反演方法的确定性和反演结果的分辨率，同时降低反演结果的随机性。

图7为第57道某点的马尔科夫链，可见经过10000次迭代后纵波阻抗的马尔科夫链基本收敛。鉴于本文利用后验均值作为模型参数的最优解，为了确保反演结果的精度，尽可能降低反演解的随机性，在二维实际资料反演时设置迭代次数为12000，并取后2000次的反演结果均值作为1次随机反演结果输出。

图8为第57道不同初始模型的反演结果。由图可见：对比不同方法的1次反演结果(图8a、图8c)表明，基于地震波形约束的MCMC随机反演结果的稳定性高于确定性反演结果，且后者使用不同的初始模型对反演结果影响不大；对比不同方法的20次反演结果(图8b、图8d)表明，不同模型的基于地震波形约束的MCMC随机反演效果好于确定性反演。为了对比基于地震波形约束的MCMC随机反演和确定性反演的效果，在下文前者采用具有地震波形指示的初始模型，后者采用普通低频模型。

图9为第57道不同样本数的反演结果。由图可见，基于地震波形约束的MCMC随机反演结果的高、低频信息丰富(图9b、图9d)，有效缓解了确定性反演精度对初始模型的依赖，且反演结果的纵向分辨率较高。

图10为基于地震波形约束的MCMC随机反演结果及初始模型。由图可见：①有效样本数为10的初始模型(图10e)效果好于有效样本数为20的初始模型(图10f)，且能反映20次反演结果(图10b)的基本特征。②有效样本数为10(图10b)、20(图10d)的20次反演结果趋势和井曲线趋势大致吻合，且前者的反演精度更高。对比原始地震记录(图11a)与有效样本数为10的反演结果合成记录(图11b)可见，两者几乎一致，也证明基于地震波形约束的MCMC随机反演结果的可靠性。

图7 第57道某点的马尔科夫链

图8 第57道不同初始模型的反演结果(a)常规初始模型1次反演结果； (b)常规初始模型20次反演结果； (c)有效样本数为10的地震 波形指示初始模型1次反演结果； (d)有效样本数为10的地震波形指示初始模型20次反演结果

图9 第57道不同样本数的反演结果(a)有效样本数为10的合成记录与原始地震记录； (b)有效样本数为10的20次反演结果； (c)有效样本数为20的合成记录与原始地震记录； (d)有效样本数为20的20次反演结果

图10 基于地震波形约束的MCMC随机反演结果及初始模型(a)有效样本数为10的1次反演结果； (b)有效样本数为10的20次反演结果； (c)有效样本数为20的1次反演结果； (d)有效样本数为20的20次反演结果； (e)有效样本数为10的初始模型； (f)有效样本数为20的初始模型 将两口验证井的纵波阻抗曲线投射在反演剖面中，为了更好地了解井曲线变化趋势，对井曲线进行了130Hz的低通滤波

图12为不同反演结果。由图中的白色虚线标记区域可见：基于地震波形约束的随机反演方法(图12a)和常规随机反演方法(图12c)均有效地提高了反演精度，对识别调谐尺度内的薄储层具有一定优势；前者在提高纵向分辨率的同时也提高了横向分辨率，且在引入波形指示初始模型的约束下，反演结果具有丰富的高、低频信息，有效地缓解了常规稀疏脉冲反演(图12b)精度对初始模型(图12d、图12e)的依赖。

图11 原始地震记录(a)与有效样本数为10的反演结果合成记录(b)

图12 不同反演结果(a)有效样本数为10的基于地震波形约束的MCMC 20次随机反演结果； (b)常规稀疏脉冲反演 结果； (c)常规随机反演结果； (d)基于波形指示初始模型的稀疏脉冲反演结果； (e)常规初始模型

4 结论

(1)本文采用基于地震波形约束的MCMC随机反演方法，利用地震波形的相似性指导测井数据建立具有波形指示的初始模型；此外，在贝叶斯框架下构建了观测地震数据和测井数据协同约束的后验概率密度分布，引入M-H采样算法多次随机模拟初始模型参数，利用后验均值作为模型参数的最优解，有效地提高了反演的稳定性和横向连续性。

(2)考虑到MCMC随机反演方法的固有性质，其反演精度和马尔科夫链迭代效率受不同建议分布以及不同模型约束条件的影响较大。为了尽可能降低这些不确定因素的影响，提高反演的稳定性，本文在MCMC随机反演方法的基础上加入波形指示的初始模型约束，可以使马尔科夫链快速收敛，间接提高了反演效率，同时也改善了随机反演精度和横向连续性。

(3)模型试算和实际资料反演效果表明，基于地震波形约束的MCMC随机反演方法具有较好的抗噪性，有效提高了反演精度，对识别调谐尺度内薄储层具有一定优势，在提高纵向分辨率的同时也提高了横向分辨率。