0 引言

在诸多薄层识别方法中，较早应用的是时频分析技术和谱分解技术。时频分析技术主要分为线性和非线性时频分析两类。1982年Morlet等^[1-2]提出小波变换并将其应用于地震数据分析，为后续的时频分析建立了基础。谱分解技术是由Partyka等^[3]在1999年率先提出并进行了实际的应用测试，通过求取地震数据中不同频率和相位的切片，准确刻画了不连续地质体，证明谱分解技术可用于薄层识别。

基于测不准原理，地球物理学家研究了小波变换、S变换、Wigner-Ville分布等方法的时频分析技术，并且在薄层识别中取得了较好的效果^[4-5]。Sun等^[6]将时频分析应用于烃类检测。Castagna等^[7]指出低频阴影与油气聚集的关系。Sinha等^[5]和张繁昌等^[8]讨论了地震数据谱分解方法。Liu^[9]与Wei等^[10]研究了谱分解在层序地层解释和烃类检测中的应用。时频分析与谱分解技术的主要思想均是最大限度地提取局部信号进而实现地层岩性解释与流体识别。但是，受调谐作用的影响，解释精度受制于地震资料的分辨率极限。

随着反演理论的不断完善和反演算法的不断发展，众多学者将反演算法引入到薄层识别中，并取得了较好的效果^[11-14]。近年来，高分辨率随机反演因其可识别调谐尺度的薄互层而广受关注^[15]。以地质统计学为基础，以测井资料为约束条件的地震随机反演方法的分辨率高于确定性反演^[16-19]，但因为计算效率低、反演结果随机性高而制约了该类方法的普及。

针对地震随机反演计算效率低的问题，王保丽等^[20]通过构建基于傅里叶滑动平均谱模拟随机反演方法，采用逐步变形算法对目标函数进行快速迭代构造，简化求解过程，大幅度提高了计算效率。张广智等^[21]通过构建一种自适应的建议分布，显著加快了马尔科夫链的收敛速度，提高了整体计算效率。李坤等^[22]提出一种建立在差分进化蒙特卡洛—马尔科夫链(MCMC)随机模型假设上的相约束叠前地震概率化反演方法，通过改变待反演参数实现对后验概率密度分布的有效模拟。

针对地震随机反演结果随机性高的问题，陈彦虎^[23]提出波形指示反演方法，该方法具有分辨率高和适用性强的特点，得到广泛应用^[24-26]。熊晓军等^[27]将波形指示反演方法应用于川中二叠系栖霞组薄储层预测。胡玮等^[28]利用该方法进行了超深层致密砂岩的薄储层预测。黄勇等^[29]通过使用多点地质统计学建模方法，优化了波形指示反演方法的建模过程，取得了更为精细的储层描述效果。鲁秀芹等^[30]使用叠前AVO数据并利用波形指示反演进行了原生煤与构造煤的空间分布预测，预测结果符合生产实际。周爽爽等^[31]提出一种基于地震波波形约束的MCMC随机反演算法，该算法充分利用地球物理映射关系，以波形相关系数为约束，通过伪普通克里金插值模拟建立一个具有地震波波形指示的初始模型，从而降低随机性。但地球物理反演固有的非线性特征仍使反演结果存在多解性。

地震波波形是地震勘探中地层结构的最直观的表现形式，蕴含着处理、解释过程所需的丰富的地下信息。对于地震波波形的研究主要包括两部分：①地震波波形属性提取及应用；②地震波波形分类技术与应用。

地震波波形属性的应用可以追溯至1963年，Churlin等^[32]通过对振幅信息的研究提出了“亮点”、“暗点”和“平点”技术。徐伯勋等^[33]在时频域提取到18种波形特征参数，初步形成了一套地震波波形特征参数提取、分析和储层预测的方法。随着计算机性能和人工智能技术的发展，地震波波形分类技术也得到了迅猛发展。Sibille等^[34]首先提出了基于波形分类进行属性分析的基本原理。Zahraa等^[35]在沉积相研究中使用了地震波波形分类技术。Priezzhev等^[36]运用Kohonen 3D神经网络和RGB技术实现了地震波波形分类的可视化应用。郑和忠^[37]开发了支持向量机的地震波波形分类应用平台。蔡涵鹏等^[38]提出了一种半监督地震波波形分类方法，利用深度学习技术分析地震相。这些研究表明地震波波形中蕴含了大量的地层岩性信息。

本文结合地震波波形匹配和波形库思想，提出基于地震波波形相似的薄互层识别方法。首先，将具有纵向高分辨率特性的测井数据引入地震波波形库的构建中，通过测井数据中的局部地层结构与井旁地震道中复合反射波形的对应关系建立两者之间的联系，即用不同的地震波波形代表不同的地层结构，形成井震对应波形库。进而基于所构建的波形库，采用波形匹配方法查找与地震道局部相似的波形，并采用相似性约束下的拟合方法重构地震道，得到匹配波形的权重值。最后根据波形库中的井震对应关系由权重值构建阻抗曲线。以地震波波形特征相似性为桥梁，实现了数据驱动下的高分辨率反演，在确保结果可靠的同时，大大提高了薄互层识别能力。

1 方法原理 1.1 建立井旁道—测井敏感曲线波形库

完备的波形库对波形分类等研究非常重要^[39]。通常有两种构建波形库的方法：一是通过解释数据约束，如井数据、层位、断层等，截取特定的地层结构模型及其对应的地震波形构成波形库；二是以井旁地震道上的复合波形为基准，截取相应的测井曲线代表特定的地层结构，然后通过时窗、步长等参数约束，自适应地截取井旁地震道的波形构建一个完善的波形库，以满足后续地震道波形匹配需求。

本文提出基于截取步长和重叠长度的双参数波形库构建方法，地震道匹配效果较好，提高了反演的准确性，同时能兼顾算法的可实现性。图 1为双参数波形截取示意图。

建立波形库的目的是为后续波形匹配提供足够可靠的样本库，因此步长参数的选择十分重要。截取步长应选择与待反演实际数据的子波长度相近；重叠长度控制了波形库中的波形数量，选取时应同时考虑反演精度及计算效率。

本文采用一种自动选取截取步长的方法，具体实现步骤如下：

(1) 扫描某一井旁道，计算该道第一组极值点的间隔L_max，则截取步长取值范围为[Trunc(L_max/4), 2L_max]；

(2) 设置重叠长度为0，生成多组波形库数据；

(3) 抽取单道地震数据，并按截取步长取值范围分别进行波形匹配；

(4) 对波形匹配结果进行相关性分析，取相关性最大的截取步长作为该工区的截取步长。

一般认为截取步长与重叠长度相差为1时，构建的是完备波形库；重叠长度小于截取步长的一半时，构建的是欠完备波形库。

地震波波形是特定地层结构下地层厚度、岩性、孔隙特征等多因素共同作用的综合响应，因此地震波波形蕴含了丰富的地下介质信息。在一个工区内部，特别是目标储层段，相似的波形代表相似的地层结构，可以用波形的相似关系推断地层岩性参数，例如波阻抗。

当目标储层对波阻抗参数不敏感时，以阻抗曲线作为测井特征曲线进行反演得到的数据体无法有效识别薄储层。因此，应选择对薄互层敏感的测井曲线，或者重构能够有效识别薄互层的阻抗曲线。图 2为井旁道—测井敏感(阻抗)曲线波形库构建流程。

为便于后续波形匹配的计算，构建波形库之前需要将井旁地震道与经过敏感性优选和重采样后的测井曲线以坐标编码的形式建立对应关系，组成井旁道—测井敏感曲线数据库，然后采用本文提出的双参数波形库构建方法建立所需波形库。

1.2 波形匹配

基于波形库的波形匹配过程是在波形库中优选与待匹配波形相近的样本波形，其实质为求取不同波形间的相似度。常见的计算相似度的方法有线性相关系数法、Manhattan距离法、欧氏距离法、切比雪夫距离法、闵可夫斯基距离法、标准化欧氏距离法、马氏距离法等^[40]。其中，线性相关系数法对地震波形状敏感性较高，对振幅值敏感性较低，Manhattan距离法对振幅差异更加敏感。波形相似性既要考虑形状的相似性，也要考虑振幅相似性。为此，联合应用Manhattan距离分析法与线性相关系数法，提高了波形相似计算的合理性。

给定两个长度为l的波形x和y，线性相关系数X的计算公式为

$ X=\frac{\sum\limits_{k=1}^{l}({x}_{k}-\overline{x})({y}_{k}-\overline{y})}{\sqrt{\sum\limits_{k=1}^{l}({x}_{k}{-\overline{x})}^{2}\cdot \sum\limits_{k=1}^{l}({y}_{k}{-\overline{y})}^{2}}} $

(1)

式中$ \overline{x} $、$ \overline{y} $分别是波形x、$ y $的平均值。

取两个等长的拥有m个采样点的波形a和波形b，累加相应振幅差的绝对值即可得到Manhattan距离，即

$ \mathrm{M}\_\mathrm{D}=\sum\limits_{k=1}^{m}\left|{a}_{k}-{b}_{k}\right| $

(2)

线性相关系数的取值范围为$ [-\mathrm{1, 1}] $，Manhattan距离的取值范围为$ [0, +\mathrm{\infty }) $。Manhattan距离需要归一化后才能与线性相关系数联合使用。归一化Manhattan距离计算公式为

$ \mathrm{M}\_\mathrm{D}\mathrm{\text{'}}=\frac{\mathrm{M}\_\mathrm{D}}{\sum\limits_{k=1}^{m}(\left|{a}_{k}\right|+\left|{b}_{k}\right|)} $

(3)

并得到联合相关系数

$ \mathrm{M}\_\mathrm{X}=X-\mathrm{M}\_\mathrm{D}\mathrm{\text{'}} $

(4)

通常无法用单个波形与目标波形匹配，为了减小匹配误差，采用上述算法在波形库中查找多个相似度较高的波形，之后本文采用普通克里金插值法进行拟合，有

$ \widehat{W}\left(k\right)=\sum\limits_{j=1}^{n}{\lambda }_{j}{w}_{j}\left(k\right) $

(5)

式中：$ \widehat{W} $为拟合波形第$ k $点处的值；$ w $为波形库中相似度最高的前$ n $个波形；$ \lambda $为权重，且满足

$ \sum\limits_{j=1}^{n}{\lambda }_{j}=1 $

(6)

普通克里金插值需要满足的假设条件为$ w $空间属性具均一性。对于空间任意一点$ (\alpha , \beta ) $，都有同样的期望c和方差$ {\sigma }^{2} $。

为了提高拟合的准确度，本文引入拉格朗日乘子，导出无偏最优估计权重求取方法。无偏估计期望值为

$ E(\widehat{W}-W)=0 $

(7)

式中W为待拟合波形。相应地，估计方差为

$ \begin{array}{l}J=\mathrm{\Delta }(\widehat{W}-W)\\ \ \ \ =\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}{\lambda }_{i}{\lambda }_{j}\mathrm{C}\mathrm{o}\mathrm{v}({w}_{i}, {w}_{j})-\\ \ \ \ \ \ \ \ 2\sum\limits_{i=1}^{n}{\lambda }_{i}\mathrm{C}\mathrm{o}\mathrm{v}({w}_{i}, W)+{\sigma }^{2}\end{array} $

(8)

式中$ \mathrm{\Delta } $代表取协方差。记$ {w}_{0}=W $，$ {C}_{ij}=\mathrm{C}\mathrm{o}\mathrm{v}({w}_{i}, {w}_{j}) $，则式(8)可改写为

$ J=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}{\lambda }_{i}{\lambda }_{j}{C}_{ij}-2\sum\limits_{i=1}^{n}{\lambda }_{i}{C}_{i0}+{\sigma }^{2} $

(9)

利用协方差公式可以将线性相关系数表示为

$ {R}_{ij}=\frac{{C}_{ij}}{\sqrt{\mathrm{\Delta }\left({w}_{i}\right)}\sqrt{\mathrm{\Delta }\left({w}_{j}\right)}} $

(10)

式中$ {R}_{ij} $为波形$ {w}_{i} $与$ {w}_{j} $的相关系数。则式(9)可写为

$ \begin{array}{l}J=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}{\lambda }_{i}{\lambda }_{j}{R}_{ij}\sqrt{\mathrm{\Delta }\left({w}_{i}\right)}\sqrt{\mathrm{\Delta }\left({w}_{j}\right)}-\\ \ \ \ \ \ \ \ 2\sum\limits_{i=1}^{n}{\lambda }_{i}{R}_{i0}\sqrt{\mathrm{\Delta }\left({w}_{i}\right)}\sqrt{\mathrm{\Delta }\left({w}_{0}\right)}+{\sigma }^{2}\\ \ \ \ =\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}{\lambda }_{i}{\lambda }_{j}{R}_{ij}{\sigma }^{2}-2\sum\limits_{i=1}^{n}{\lambda }_{i}{R}_{i0}{\sigma }^{2}+{\sigma }^{2}\end{array} $

(11)

最终目标是寻找一组$ {\lambda }_{i} $，使得$ J $最小。因为$ J $是$ {\lambda }_{i} $的函数，故该问题可化简为求解$ J $对$ {\lambda }_{i} $求偏导数为零的方程，即

$ \frac{\partial J}{\partial {\lambda }_{i}}=0 $

(12)

对此，可采用等式约束—拉格朗日乘数法求解。引入拉格朗日乘子$ \vartheta $，则目标函数为

$ F=J+2\vartheta \left(\sum\limits_{i=1}^{n}{\lambda }_{i}-1\right) $

(13)

对该目标函数求解一组$ \vartheta $和$ \lambda $值，使$ J $最小，即

$ \frac{\partial \left[J+2\vartheta \left(\sum\limits_{i=1}^{n}{\lambda }_{i}-1\right)\right]}{\partial {\lambda }_{j}}=0\ \ \ \ j=\mathrm{1, 2}, \cdots , n $

(14)

$ \frac{\partial \left[J+2\vartheta \left(\sum\limits_{i=1}^{n}{\lambda }_{i}-1\right)\right]}{\partial \vartheta }=0 $

(15)

化简可得

$ \left\{\begin{array}{l}{R}_{11}{\lambda }_{1}+{R}_{12}{\lambda }_{2}+\cdots +{R}_{1n}{\lambda }_{n}+\frac{\vartheta }{{\sigma }^{2}}={R}_{10}\\ {R}_{21}{\lambda }_{1}+{R}_{22}{\lambda }_{2}+\cdots +{R}_{2n}{\lambda }_{n}+\frac{\vartheta }{{\sigma }^{2}}={R}_{20}\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ⋮\\ {R}_{n1}{\lambda }_{1}+{R}_{n2}{\lambda }_{2}+\cdots +{R}_{nn}{\lambda }_{n}+\frac{\vartheta }{{\sigma }^{2}}={R}_{n0}\end{array}\right. $

(16)

用矩阵形式表示为

$ \left[\begin{array}{ccccc} R_{11} & R_{12} & \cdots & R_{1 n} & 1 \\ R_{21} & R_{22} & \cdots & R_{2 n} & 1 \\ & & \vdots & & \\ R_{n 1} & R_{n 2} & \cdots & R_{n n} & 1 \\ 1 & 1 & \cdots & 1 & 0 \end{array}\right]\left[\begin{array}{c} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_n \\ \frac{\vartheta}{\sigma^2} \end{array}\right]=\left[\begin{array}{c} R_{n 0} \\ R_{n 0} \\ \vdots \\ R_{n 0} \\ 1 \end{array}\right] $

(17)

求解该矩阵即可得到权重值$ \lambda $，进而完成曲线拟合计算。图 3为欠完备波形库下使用不同波形相似度计算方法所得曲线拟合结果，对比可知，联合相关系数法曲线拟合结果（拟合曲线2）的效果优于线性相关系数法（拟合曲线1）。

1.3 阻抗反演与技术流程

基于井旁道波形匹配拟合获得权重系数，由波形库中对应的阻抗曲线计算波阻抗，即可实现阻抗反演，计算公式为

$ \stackrel{}{\widehat{\mathrm{I}\mathrm{m}\mathrm{p}}}\left(k\right)=\sum\limits_{j=1}^{n}{\lambda }_{j}\mathrm{I}\mathrm{m}\mathrm{p}_{j}\left(k\right) $

(18)

式中：$ \widehat{\mathrm{I}\mathrm{m}\mathrm{p}}\left(k\right) $为待反演地震道局部波形在第k个点处的阻抗值；$ {\lambda }_{j} $为通过优选得到的第j个井旁道波形对应的拟合权重；$ \mathrm{I}\mathrm{m}{\mathrm{p}}_{j}\left(k\right) $为第j个井旁道波形对应的阻抗波形上第k个点处的阻抗值；$ k=\mathrm{1, 2}, \cdots , \mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p} $，step为截取步长。

基于地震波波形相似的薄互层识别方法技术流程如图 4所示。

（1）数据预处理。利用测井解释成果及测井曲线进行敏感曲线分析，根据需要进行敏感曲线重构，并进行精确井—震标定。

（2）构建井旁道—测井敏感曲线波形库。自适应求取截取步长，根据实际需要设置重叠长度，采用双参数对井旁地震道和对应测井敏感曲线进行波形裁剪及坐标编码，形成井旁道—测井敏感曲线波形库。

（3）计算拟合权重。计算相似度，优选前n个相似度最高的波形，或者通过设置阈值优选波形，利用无偏最优估计求取权重值。

（4）输出高分辨率反演体。根据井—震坐标关系对所需测井敏感参数值进行解码，利用第（3）步所得的权重值计算目标点处的反演参数值，最终得到高分辨率反演数据体。

2 二维模型试算

截取局部Marmousi模型作为理论模型进行实验分析。模型纵向为370个样点，采样率为1 ms，时间长度为370 ms，横向为1000个CDP。

图 5为Marmousi模型速度场及利用主频30 Hz的零相位雷克子波与模型纵波阻抗进行褶积得到的原始合成地震记录。利用本文方法对该合成记录进行测试，具体试算过程及参数设置如下。

（1）抽取第100、第200、…、第800道纵波阻抗数据作为伪井数据。

（2）设置截取步长为55个样点、重叠长度为40个样点，建立井旁道—阻抗波形库。

（3）输入截取步长长度的待反演地震道波形，选用联合相关系数法计算相似度；若相似度存在大于-0.9的值则将其对应的地震波波形取出作为待拟合波形，否则取3个相似度最大值对应的地震波波形作为待拟合波形。

（4）利用式(17)计算拟合权重。

（5）解码地震波波形对应的阻抗波形，然后按照步骤（4）所求权重进行阻抗反演。

（6）重复步骤（3）~步骤（5），直至该待反演地震道所有采样点对应的阻抗值求取完毕，然后进行下一地震道计算，直至所有地震道求取完毕，即可得到阻抗剖面。

图 6为第550道阻抗反演结果，与原始阻抗值对比可见，利用本文算法较好地还原了阻抗曲线。图 7为原始纵波阻抗剖面（图 7a）、反演阻抗剖面（图 7b）及其分别对应的合成记录（图 7c、图 7d）对比，虽然图 7b相对图 7a多出了一些薄层，但这些薄层的阻抗差非常小，可以忽略，反演阻抗合成地震记录与原始地震记录基本一致，表明基于地震波波形相似的薄互层识别方法具有较高可行性。

3 实际资料应用

实际工区内共有39口有效井，采样率为1 ms，采样点为901个，每条测线有110道，反演时窗为550~1450 ms。将经过精细井—震标定的39口井的波阻抗曲线及井旁地震道作为样本，选用截取步长为151、重叠长度为90建立波形库。选取过5口井的Inline13剖面为研究目标对比分析处理效果。从图 8a原始剖面可以看到，测井曲线长度并未覆盖整个反演时窗。本方法利用波形库思想的优势，可以在少井或测井深度未覆盖的区域实现薄互层的识别。因此与原始地震剖面对比可以看出，反演阻抗剖面（图 8b）的分辨率得到明显提升。

图 9为第50道反演阻抗曲线求取反射系数后与40 Hz雷克子波褶积得到的合成地震记录与原始地震记录对比。通过对比可以清楚看到两者趋势基本一致，验证了本文方法的可靠性。

通过层位约束截取目标层段，应用本文方法对目标层反演，前、后剖面对比如图 10所示。将测井解释得到的砂、泥岩岩性柱投影在反演剖面上（图 10b）。与测井解释结论对比可见：反演剖面上的层位信息与砂、泥岩钻井分层吻合度较高，相对于原始地震数据（图 10a）有效提高了薄互层识别能力，但是反演剖面的分辨率距离测井的纵向分辨率还有较大的差距。

4 结论

(1) 本文提出基于地震波波形相似的薄互层识别方法，采用波形库思想和波形匹配技术，充分发挥了测井、地震数据分别具有的纵、横向高分辨率优势，以波形相似性为驱动，避免求取地震子波，降低了反演结果的多解性，大大提高了处理结果的可靠性和稳定性。

(2) 模型试算和实际资料测试结果表明，基于地震波波形相似的薄互层识别方法具有较好的可行性和适用性，有效提高了纵向分辨率，对识别调谐尺度内的薄互层具有现实意义。

(3) 本文方法的核心是建立地震波波形与测井阻抗曲线(薄互层地层结构)的对应关系，这种关系在稳定的沉积相内部或者目的层段是可靠的，如果垂向上时间跨度太大，或者横向上沉积环境发生大的变化，则反演效果可能会下降甚至无效。因此测井信息的时空分布是本方法成功应用的关键。