地震波形指示反演方法及其应用

陈彦虎，毕建军，邱小斌，陈有兵，杨辉， 曹佳佳，邸永香，赵海山，李志向

（1. 北京中恒利华石油技术研究所，北京 100102；2. 中国矿业大学（北京），北京 100083）

0 引言

地震反演技术在过去40 余年间得到了长足发展，赵政璋[1]、印兴耀[2]、撒利明[3]、甘利灯[4]等专家学者系统总结了地震反演技术的发展现状。目前石油工业界应用较为广泛的地震反演技术主要包括3 大类：第1类是基于褶积模型的地球物理反演[5]，其发展经历了从叠后“亮点”技术、波阻抗反演技术[6]到叠前AVO 技术[7]、弹性波阻抗反演技术[8]，近年来随着宽方位地震技术的发展，地震反演进一步拓展到了OVT（Offset Vector Tile，炮检距向量片）域[9]；第2 类是非线性反演技术，包括神经网络、支持向量机、贝叶斯、模式识别和遗传算法等[3]；第3 类是地质统计学反演，通过随机模拟实现提高纵向分辨率[10]。

上述反演技术在储集层预测领域取得了良好的效果，但由于各自不同的方法原理，均存在一定局限性：基于褶积模型的地球物理反演无法突破地震分辨率的限制[11]；非线性反演得到的解往往非全局最优，反演结果地质规律性较差[12]；地质统计学反演高频成分来自随机模拟，反演结果随机性强，横向分辨率低[13]。

中国陆相沉积储集层厚度小、横向变化快[14]，目前对储集层识别精度的需求已达1 m 左右[15]，现有的地震反演技术无法满足这一精度要求。针对该难题，本文在详细分析地震波形与测井高频信息内在联系的基础上，系统阐述了地震波形指示反演的方法原理，利用地震波形横向相似性驱动测井高频信息，实现了高分辨率反演。通过模型验证以及大庆长垣薄互层识别实例分析，认为波形指示反演能够较好地解决陆相薄互层砂体预测的难题。

1 地震波形指示反演方法的理论基础

1.1 地震纵向分辨率和横向分辨率探讨

地震分辨率分为纵向分辨率和横向分辨率两个方面。纵向分辨率的概念最早由Rayleigh[16]给出，即相邻两个反射界面的分辨率极限为1/4 波长，厚度小于1/4 波长的地质体即可以定义为薄层，利用地震数据无法分辨。为解决薄储集层预测问题，许多学者开展了大量研究，Widess[17]提出在理想情况下，根据振幅横向变化能够识别任意厚度的薄层；曾洪流[18]提出了“横向检测率”的概念，利用地层切片可以检测厚度小于1/4 波长的地质体横向变化。

通过一个正演模型来讨论薄层的纵向分辨率和横向分辨率。地质模型为：在泥岩背景下，连续发育一套厚度为3 m 的砂岩，砂岩之上发育一套厚度为0～3 m 的透镜状薄砂岩（见图1a），泥岩夹层厚度为3 m，砂泥岩地震波速度分别为3 500 m/s 和2 800 m/s，密度分别为2.65 g/cm3和2.26 g/cm3。利用35 Hz 零相位雷克子波进行褶积。从得到的正演剖面（见图1b）可以看出，纵向上无法直接分辨出上覆薄砂岩。但由于薄砂岩的发育，地震波形横向上发生了变化，提取图1b 中最大波峰处的振幅和频率（见图1c），可以看出地震振幅和频率横向上发生了很大变化，即地震波形包含了薄层信息。因此可以得出结论：纵向上地震振幅无法“分辨”薄层，但可以利用横向波形变化“识别”薄层。

1.2 地震波形与测井高频信息的内在联系

地震波形包含了地震运动学和动力学多种信息，是地质沉积作用、岩性岩相组合，储集层物性以及流体等多种地质信息的综合响应。曾洪流[19]提出了“地震沉积学”的概念，充分利用地震资料的横向高分辨率，提高了薄层的预测精度，在薄储集层预测中得到了广泛应用。

本文在地震沉积学技术的基础上，引入测井资料，将地震资料的横向高分辨率和测井曲线的纵向高分辨率有机结合起来，首先要解决的问题是如何建立地震 波形和测井高频信息之间的内在联系。

图1 地质模型与正演剖面

在实际资料中，相似的沉积特征往往具有相似的岩性组合，相似的岩性组合往往具有相似的地震波形特征。选取一个陆相薄储集层实际资料进行分析，图2a 为A、B 两口井的地震剖面，该地震资料信噪比高，但频率较低，主频约为20 Hz，频宽约为10～35 Hz。提取过两口井同一层段井旁道地震波形，由于两口井存在深度差异，因此首先对两口井的地震波形经过顶面对齐和厚度一致性校正，然后进行叠合对比（见图2b），从图中可以看出，两者相关系数达94%。因此，两口井的井旁地震道具有相似的地震波形特征。

图2 A 井和B 井地震波形特征

再进一步分析测井曲线的相似性特征：首先分析与地震反射直接相关的纵波波阻抗曲线，两口井原始纵波波阻抗曲线相关系数为75%（见图3a），远低于地震波形的相似性，分析其原因，测井曲线具有非常高的频率，高频成分的差异降低了其相似性。依次对测井曲线进行滤波，分别保留0～500 Hz、0～400 Hz、0～300 Hz、0～200 Hz 以及0～100 Hz（见图3b—图3f），相关系数逐步提高到88%、92%、93%、94%和95%，达到了非常高的相似性。由此可见，由于测井曲线高频信息的差异性导致原始测井曲线相似性较低，对测井曲线逐步降低频率进行滤波，测井曲线相似性逐步提高，当测井曲线频率达200～300 Hz 时，相关系数就达到了地震波形的相似性，因此可以建立起地震波形与测井高频信息之间的联系，为高分辨率地震反演提供了依据。开展地震波形指示反演时，可以通过分析不同频率下测井曲线相似性和目标储集层的厚度共同确定反演的最佳截止频率。

当储集层特征复杂、纵波波阻抗无法区分岩性时，需要借助自然伽马等曲线区分岩性，但是基于褶积模型的地震反演无法预测自然伽马曲线。借助分析纵波波阻抗曲线相似性的思路来分析自然伽马曲线，如图4a，两口井原始自然伽马曲线相关系数只有71%，依次对自然伽马曲线进行滤波，保留0～500 Hz、0～400 Hz、0～300 Hz、0～200 Hz 以及0～100 Hz 频率成分（见图4b—图4f），相关系数逐步提高到88%、90%、92%、93%和95%，达到了和纵波波阻抗基本一致的相关性。

图3 A、B 井纵波波阻抗曲线不同频带范围滤波

从沉积学的角度分析，在等时地层格架内，地震波形的变化反映了岩性组合的变化，岩性组合是沉积 相或者地震相的表现形式。从测井曲线的角度分析，波阻抗和自然伽马等其他曲线均可以反映岩性组合或沉积相的变化。因此，低频地震信息与高频测井信息的对应关系不仅适用于波阻抗曲线，同样适用于自然伽马等非波阻抗曲线，这就奠定了利用低频地震资料和高频测井资料信息进行高分辨率反演（或模拟）的基础。需要说明的是，对于波阻抗曲线为地震波形指示反演，而对于自然伽马等非波阻抗曲线为地震波形指示模拟。

2 地震波形指示反演方法

上节实际数据地震波形与测井曲线的分析表明，相似波形对应的测井曲线在较宽频带内呈现较高的相似性，因此利用地震波形横向相似性驱动高频测井信息实现高分辨率反演[20]，建立了地震波形指示反演方法（Seismic Meme Inversion，简称SMI）。其中，“Meme”的涵义在《牛津英语词典》中定义为通过模仿传递信息，即通过相似地震波形的模仿，传递地震波形代表的高频测井曲线信息，实现高分辨率反演。

地震波形指示反演实现过程中，首先通过奇异值分解实现井旁地震道波形动态聚类分析，建立地震波形结构与测井曲线结构的映射关系，生成不同类型波形结构（代表不同类型的地震相）的测井曲线样本集；然后通过分析不同类型波形结构对应的样本集分布，分别建立不同地震相类型的贝叶斯反演框架；然后在不同贝叶斯框架下，分别优选样本集的共性部分作为初始模型进行迭代反演；最后在反演迭代过程中，以样本集的最佳截止频率为约束条件，得到高分辨率的反演结果。地震波形指示反演方法的基本原理主要包括以下3 个方面。

2.1 利用奇异值分解实现地震波形聚类

地震波形参数和井点属性的对应关系可以定义为一个n m× 阶矩阵A，则有：

式中，U 和V 为地震波形数据和井点属性数据，分别为n n× 阶和m m× 阶正交矩阵，VT为V 的共轭转置，Λ为n m× 阶非负实数对角矩阵，表示地震波形数据和井点数据的相关性。

奇异值分解即对A 进行正交分解，当矩阵的秩为r 时，则矩阵A 可分解为r 个本征向量的代数和，则矩阵A 的总能量可表示为：

矩阵A 经过（2）式进行奇异值分解之后，可以用前r 个非零奇异值对应的奇异向量表示矩阵A 的主要特征，即实现地震波形动态聚类分析，得到不同储集层类型地震波形与测井曲线特征的对应关系，建立初始的样本集。

2.2 测井曲线小波变换

小波变换将数字信号从时间域变换到频率域，可以表征信号在时间域和频率域的特征，能够同时定量预测信号的低频稳定部分和高频突变部分。

对前文建立的样本集中的测井曲线分别用不同的截止频率l 利用（3）式开展离散小波变换，将测井曲线分解为低—中频宏观特征信息、高频细节信息和超高频噪声信息3 部分，需要说明的是，由于测井曲线具有非常高的频率，以声波时差曲线为例，其频率高达20 kHz[21]，因此，这里定义的低—中频大致为100～200 Hz，甚至300～400 Hz 的频率范围，即本文中相对于地震频率范围定义的“测井高频信息”，提取的低—中频宏观特征信息即该样本集中所有测井曲线的共性结构，该共性结构可以作为波形指示反演的初始模型。

2.3 贝叶斯框架下的波形指示反演

地震反演的基础是褶积模型为：

假设噪声n 满足高斯分布：

将（5）式代入（4）式，建立地震数据似然函数：

在贝叶斯反演中，假设弹性参数模型m 也符合高斯分布，可以得到模型的先验分布为：

将数据条件概率分布与模型先验概率分布的乘积作为模型的后验概率分布函数，其具体表达式为：

这样对于给定的地震波形d，可按照后验概率应用Gibbs 抽样法计算得到模型m 的期望值。

（8）式中，概率最大时得到的解，即为反演的最终解。对（8）式两端求对数，得到目标函数：

为了使后验概率最大，对式（8）关于模型参数Δm求导得到：

使用迭代模型扰动量的方法逼近样本数据，得到最终的高分辨率反演结果。

地震波形指示反演方法具有以下特点：①地震波形指示反演利用地震波形横向变化驱动高频测井信息特征实现高分辨率反演。反演结果纵向上与测井高频信息吻合，具有纵向高分辨率（当测井曲线高频共性结构为200～300 Hz 时，薄层预测的分辨率为2～3 m）。同时反演结果横向上遵循地震波形的变化，也具有横向高分辨率，因此，地震波形指示反演可以同时提高反演结果的纵横向分辨率，是一种高精度的反演方法。②井间唯一存在的数据为地震数据，地震波形的横向变化体现了岩相组合的变化。地震波形指示反演通过地震波形横向变化代替变差函数空间域插值模拟，实现了地震相自动控制下的反演，克服了传统相控反演需要人为事先给定沉积相而导致的主观性，是真正意义上的相控反演。③地震波形指示反演方法不仅可以实现高分辨率波阻抗反演，还可以实现自然伽马、电阻率和孔隙度等非波阻抗参数的相控高分辨率模拟，突破了地震反演只能得到波阻抗结果的局限性，对于利用地震信息进行储集层参数模拟是一次巨大的进步。

3 方法验证

为了验证波形指示反演方法的合理性和对薄互层的识别能力，建立了薄互层正演模型进行反演实验。

地质模型为（见图5）：泥岩背景中发育4 组薄互层砂体，第1 组为3 个厚度为3 m 的砂岩，泥岩夹层厚度为3 m，后面3 组分别去掉其中1 个砂体。其中，砂泥岩速度分别为3 500 m/s 和2 800 m/s，密度分别为2.65 g/cm3和2.26 g/cm3。为了便于开展波形指示反演，在该模型上建立W1—W9 共9 口虚拟井，代表不同的储集层特征。

图5 薄互层地质模型

利用主频30 Hz 到120 Hz 的零相位雷克子波对地质模型进行褶积，选取主频为30，60，120 Hz 分别代表低频、中频和高频的地震正演剖面开展不同反演方法实验。图6a 至图6c 分别为主频30，60，120 Hz 的地震正演剖面，黄色曲线为纵波波阻抗曲线。从图中可以看出，30 Hz 地震资料由于分辨率较低，完全无法识别薄互层砂体，但由于砂体组合不同，地震波形差异也较大；60 Hz 地震资料只能识别出第3 组砂体；直到120 Hz 地震资料才可以识别4 组砂体的每一个薄 砂体。

首先利用正演地震数据开展常规稀疏脉冲反演，图6d—图6f 分别为主频30，60，120 Hz 的稀疏脉冲反演剖面，从图中可以看出，稀疏脉冲反演结果表现出了和地震资料近似的分辨率：30 Hz 反演结果完全无法识别薄互层砂体；60 Hz 反演结果可以分辨出第3 组砂体的两个砂体；直到120 Hz 反演结果才可以分辨出 4 组砂体的每一个单砂体。

图6 不同主频正演地震数据及不同方法反演剖面

然后再利用正演地震数据开展地震波形指示反演，图6g—图6i 分别为主频30、60 和120 Hz 的波形指示反演结果，从图中可以看出，3 个结果分辨率基本相同，都能清晰地识别4 组砂体的每个单砂体。

为进一步验证地震波形指示反演的抗噪性，对图6a 中主频为30 Hz 的地震道分别加入5%、10%和20%的随机噪声（见图7a—图7c），然后分别进行地震波形指示反演。从得到的纵波波阻抗结果（见图7d—图f）可以看出，在3 种不同的噪声水平下，反演结果均可以准确地识别出3 m 厚的薄砂岩储集层，只是在泥岩背景下，出现了不同程度的噪声，但不影响目标砂岩的准确识别，证实了地震波形指示反演方法具有较强的抗噪性。

模型验证结果表明，稀疏脉冲反演无法识别厚度小于1/4 波长的薄互层砂体。地震波形指示反演能够突破1/4 波长的限制，只要地震波形横向有差异，即使在地震资料分辨率较低的情况下，也能识别不同的薄互层砂体。同时加随机噪声反演测试结果表明，地震波形指示反演具有较强的抗噪性。

4 应用实例

4.1 储集层特征

研究区位于松辽盆地大庆长垣北部白垩系姚家组葡萄花油层。前人研究认为，该地层沉积时期大庆长垣北部地区受北部水系控制形成大型河流—三角洲复合体，广泛发育近南北向水下分流河道沉积[22]。砂体埋藏深度为900～1 200 m，单砂体厚度为0.2～15.0 m，47%的砂体厚度为1～3 m（见图8a），呈现出典型的薄互层特征。油田进入开发后期阶段，挖潜的主要对象从厚层砂体转变为厚度小、横向变化快的砂体，因此需要精细刻画薄互层砂体。

图7 加随机噪声的30 Hz 正演地震数据与地震波形指示反演剖面

图8 砂体厚度特征和岩性解释图版

利用测、录井资料开展岩性测井解释，建立砂泥岩解释图版（见图8b）。砂岩纵波波阻抗小于7 600 m·g/(s·cm3)，因此，纵波波阻抗可以识别砂泥岩，可以开展地震波形指示反演得到纵波波阻抗结果预测砂岩分布。

4.2 地震波形指示反演方法的应用

本次研究地震三维工区面积6.25 km2，工区内有钻井39 口，且分布较为均匀。为验证波形指示反演算法的合理性，开展不同数量井参与反演实验：第1 组实验39 口井全部参与反演（见图9a、图9b）；第2 组实验模拟勘探阶段少井的情况，从39 口井中随机选取10 口井参与反演，其余29 口井为验证井（见图9c、图9d）；第3 组实验模拟评价阶段不均匀井网的情况，北边25 口井参与反演，南部14 口井为验证井（见图9e、图9f）。从上述3 组实验的反演剖面（见图9a、图9c、图9e）上可以看出，无论是参与井还是验证井，反演结果与钻井吻合程度都比较高，并且3 组实验反演结果高度相似。按照岩性直方图确定的砂岩门槛值，提取了上油层组砂岩（反演剖面中黑色虚线框，见图9a、图9c、图9e）厚度图（见图9b、图9d、图9f），从图中可见，3 个厚度图形态基本一致，都表现出了近南北向发育的两个河道砂体分支，主河道边界清晰，河道展布形态符合地质规律。

为进一步定量描述反演预测的精度，提取反演预测的上油组砂体厚度，与测井解释砂岩厚度进行对比分析。图10a 为39 口井全部参与反演对比结果，反演预测砂岩厚度与测井解释厚度相关系数达90.8%；图10b 为模拟勘探阶段随机选取10 口井参与反演的结果，红色点为参与反演井，相关系数达91.2%，蓝色点为验证井结果，相关系数达80%。图10c 模拟评价阶段北边25 口井参与反演结果，红色点为参与反演井，相关系数达90.2%，蓝色点为验证井结果，相关系数达80.8%。

通过上述不同方法反演结果对比以及不同数量井参与反演结果表明，波形指示反演结果具有两方面的优势：①分辨率高，能识别厚度为2～3 m 的薄互层，同时反演精度高，参与井吻合率达90%，验证井吻合率达80%；②地震波形驱动测井曲线实现高分辨率反演，具有相控思想，不受井数量和井位分布影响，反演结果符合地质规律。

4.3 不同反演方法对比

地震波形指示反演对大庆长垣薄互层的识别精度达到了2 m，通过地震波形指示反演与主流的确定性反演方法——稀疏脉冲反演结果的对比分析，进一步证明 地震波形指示反演具有更高的精度。

图9 不同数量井参与反演结果对比

图10 反演预测砂体厚度与测井解释砂岩厚度交会图

目的层三维地震资料主频为40 Hz，纵波速度为 3 000 m/s 左右，按照地震分辨率1/4 波长理论，地震资料可识别的最大厚度约为18.75 m，显然不满足薄互层砂体识别的要求。从过井的地震剖面上（见图11a）可以看到，地震资料无法识别葡萄花油层的薄互层砂体。首先开展常规稀疏脉冲反演，从得到的纵波波阻抗剖面上（见图11b）可以清楚地看到，稀疏脉冲反演分辨率与地震资料大致相当，只能大致识别大套的砂层组，无法预测薄互层砂体。然后开展波形指示反演，从得到的纵波波阻抗反演剖面上（见图11c），可以直观地看到，波形指示反演结果的纵向分辨率明显高于 稀疏脉冲反演，能够识别薄互层砂岩，同时，砂岩横向变化特征符合地震特征，能够较好地刻画砂体横向边界，即波形指示反演也具有较高的横向分辨率。

图11 地震剖面与不同方法反演效果对比

为进一步定量论证两种反演方法的精度，分别从两个反演数据体上提取井点处曲线（以图11 中W3 井为例），如图12 中第1、2 道红色曲线为实测纵波波阻抗曲线（为便于和反演结果对比，实际纵波波阻抗曲线以2 m 采样率进行了滤波），可以看出，波形指示反演结果与实测曲线吻合度远高于稀疏脉冲反演。分别对两种反演结果按照砂岩门槛值进行岩性解释，与实际测井岩性解释结果对比（见图12）可以看出，波形指示反演解释岩性与测井解释岩性吻合度高，可识别厚度为1.8～3.0 m 的砂岩；而稀疏脉冲反演结果与测井解释岩性吻合度低，可识别的砂体厚度为10 m。由 此可见，波形指示反演无论是吻合率还是反演精度都远高于稀疏脉冲反演。

图12 稀疏脉冲反演与波形指示反演结果对比

5 结论

地震波形指示反演通过地震波形动态聚类分析，建立地震波形结构与高频测井曲线结构的映射关系，提高了反演结果的纵向分辨率和横向分辨率；通过构建不同地震相类型的贝叶斯反演框架，实现了真正意义上的相控反演。

地震波形指示反演的模型验证结果和大庆长垣的实例应用表明，地震波形指示反演对薄互层的识别精度达到了2 m，反演结果与钻井的吻合率超过80%。地震波形指示反演能够提高储集层预测精度，为陆 相盆地横向相变快的薄互层砂岩识别提供了借鉴和 参考。

同时，大庆长垣实例应用还表明，地震波形指示反演适用于勘探、评价以及开发等各个阶段。但是本方法的基本原理要求已钻井尽可能覆盖研究区主要沉积相类型，因此，实际应用中需要一定数量的已钻井，不适用于勘探早期阶段井少的情况。地震波形指示反演同样适用于碳酸盐岩和火山岩储集层，当碳酸盐岩和火山岩储集层为层状介质时，反演的要求和注意事项与碎屑岩储集层一样；当碳酸盐岩或火山岩储集层为丘状介质时，需要通过精细的层位解释构建出丘状介质的格架模型，才能获得理想的反演结果。 符号注释：

A——地震波形数据和井点属性数据矩阵；d——地震数据矩阵；G——地震子波矩阵；GR——自然伽马曲线，API；I——先验信息矩阵；l——测井曲线共性结构相关截止频率，Hz；m——求解的弹性参数模型矩阵；n——噪声矩阵；N——样本集的总数；r——矩阵的秩；R——反演预测砂体厚度与测井解释砂岩厚度相关系数，%；U——地震波形数据矩阵；V——井点属性数据矩阵；W——样本井集测井曲线；W ——样本井集测井曲线平均值；Λ——地震波形数据和井点属性数据的相关性矩阵；δi——AAT特征值的非负平方根；Δm——模型参数的扰动量； ( , )tψ ω ——小波函数；σ——地震数据的协方差；σm——模型的方差；σΔm——模型参数扰动量的方差。