一种动态匹配特征子波拾取地震同相轴的技术

彭仁艳，徐振旺，刘文峰，余 锋，董旭光

1.中国石化石油工程地球物理有限公司科技研发中心，江苏 南京 210000

2.中国石油辽河油田分公司勘探开发研究院，辽宁 盘锦 124010

3.中国石化石油工程地球物理有限公司华东分公司，江苏 南京 210000

4.中国石化石油工程地球物理有限公司国际业务发展中心，江苏 南京 210000

引言

地震同相轴自动识别与拾取在速度分析、资料解释等方面具有重要的应用价值，但对于大规模三维地震数据，人工拾取工作量巨大，效率低下，无法满足工业生产的需要[1-8]。在地震剖面上自动识别与拾取同相轴方面，前人已进行了大量研究，基本思想为基于数据的空间分布特征，利用局部相似性寻找线性相关的同相轴[9-14]。根据计算局部相似性算法的不同对拾取方法进行分类，目前，工业界中使用最广泛的是基于倾角信息的拾取算法，即根据地震数据的空间局部相似性，估计整个数据空间的倾角场，然后，根据参考道或种子点，沿着估计的倾角方向依次追踪相邻道，最后，将拾取的点进行模式分类，构造出相应的反射层位。其中，描述局部相关性，计算倾角的方法是整个自动拾取的关键步骤，其计算方法也有很多种[15]，如利用瞬时波数方向估计反射层的倾角[16]、利用平面波分解或预测计算局部空间的倾角[17]、利用结构张量的方法估计倾角[18]等。平面波分解假设同相轴是局部线性的，计算精度不高；结构张量对数据比较敏感，易受噪音的干扰。另外一种比较常用的自动拾取方法是基于多项式拟合的拾取算法[19]，此类方法大多用在成像道集上的剩余时间（深度）差拾取，通过构造多项式系数拟合剩余曲率曲线，来完成拾取的工作，这种方法基于多项式假设，不能适用于任意地震剖面拾取情况。

动态波形匹配是广泛用于声波信号匹配、字典匹配等领域的自动匹配算法，它提供了良好的求取两列信号最佳匹配的方法。动态波形匹配没有考虑地震信号的特点，本文结合地震信号的子波特征，对地震信号进行矢量化特征表达，在此基础上结合动态波形匹配算法，提出一种同相轴自动识别与拾取方法。介绍了动态波形匹配的基本原理，在此基础上，将子波特征矢量化的数据利用动态波形匹配算法进行层位自动追踪，实现基于子波特征的动态波形匹配同相轴自动拾取。给出了理论数据测试实例和实际资料测试实例，证明了方法的正确性和实用性。

1 动态波形匹配法基本原理

动态波形匹配是互相关运算的推广，是波形特征提取的重要算法，早期在声波匹配、字典对比等方面应用较为普遍。近年来，地球物理学家将该算法推广到地震数据处理中，并提出动态时间规整方法，在去除地震数据的静校正量、校平地震道集及匹配纵横波成像深度差等方面均得到成功应用[20-25]。动态波形匹配通过计算波形匹配最短路径，把波形匹配最短路径的大小作为相邻地震道波形之间差异大小程度的度量。

动态波形匹配的目的是找到一个时移量的序列，使得两个序列f和g 相匹配后的误差和达到最小。令f为参考道，g 与f相匹配，定义eili表示对时间序列g 的第i个采样点进行时移li后与时间序列f第i个采样点之间的距离，其表达式为

动态时间规整DTW 算法通过动态规划法进行求解，分为累积和回溯两步计算。

累积过程中，利用校正误差eili递归地计算距离累积误差dili，假设时间应变的约束条件为|ui−ui−1|≤1，则累积过程为

在累积过程的最后一步，即i=N−1 时，计算出所有时差l对应的距离函数后，找到最小的距离，就可以得到最终的最小总距离函数

通过往前递归地回溯，最后可以找到每一个采样点对应的时差，即最短路径。

2 基于子波特征表达的动态波形匹配地震剖面自动拾取

动态波形匹配通过定义点与点之间距离，求解点与点之间的最短路径。该方法没有考虑地震信号的特征。根据褶积模型，地震数据可以视为反射系数序列与子波的褶积加随机噪音组成，即

地震剖面拾取就是识别主要的反射同相轴，即提取主要反射系数序列。地震子波及噪音的存在对自动拾取带来较大的干扰，增加了自动拾取的难度。地震子波数学上就是一小段波动序列，由于薄层调谐效应、吸收衰减以及噪音的影响，在实际地震剖面中子波发生较大变化。为了减少子波对自动拾取的影响，对地震子波进行特征提取，将地震子波抽象为由波峰位置、波宽、振幅及特征斜率（主瓣与副瓣幅值比）等属性组成的序列。抽象过程表示为矢量形式为

式中：p--子波波峰/波谷极值位置；q--子波波宽（子波持续时间/深度范围）；a--子波波峰/波谷极值大小；v--子波特征斜率；iii，jjj，kkk，lll--单位向量。

式（9）为子波特征矢量。对地震剖面s(t) 进行子波提取，并用子波特征矢量表示

式中：

s(t)--长度为t的地震剖面；ccctc长度为tc的特征矢量剖面数据，即选取适当振幅（波峰/波谷）阈值，通过子波特征矢量表达将长度为t的地震剖面转化为长度为tc的矢量数据。

式（10）将地震数据转化为子波特征矢量数据，将一维地震道转化为四维特征序列数据。需要指出的是，式（10）中tc通常远小于t，虽然将一维数据扩展为四维，但是转换后的数据规模得到大大压缩，能够显著提高数据的后续处理效率。子波特征矢量表达对地震剖面进行了特征的提取，降低了剖面中子波以及噪音的影响，对地震同相轴的识别和提取具有重要作用。

为了进行地震同相轴的自动拾取，需要对特征表达后的地震数据进行自动追踪。定义emn,klnk为特征表达矢量数据序列cmn和cm+k,n+lnk之间的距离，根据矢量距离公式，定义为

式中：

ωp，ωq，ωa，ωvp,q,a,v 相应的加权系数；m--道号；j对应道的样点数；k--道号移动量；lnk采样点移动量，通常，|k|=1 表示相邻两道进行波形匹配；|l|＜∆--相邻子波的距离，根据不同的数据特征进行选择，一般大于同一地震道相邻子波的最大间距。

在实际处理中，会逐地震道依次迭代处理。不失一般性，后面的公式中忽略m，并令k=1，emn,klnk简化为en,ln。

用动态规划法求解满足en,ln的最短路径。首先，递归地计算距离累积误差dn,ln，误差积累过程，表达式如下

式中：ε--误差阈值。

当两个子波之间的误差小于阈值ε 时，认为该两个子波处在同一个同相轴上，令其累积误差最小，就可以把这两个子波先匹配起来。当输入数据信噪比低，子波特征变化不明显时，可适当减小ε 值。对于每一个参考位置子波，都存在一个最短路径找到另一道上与其波形特征最为逼近的子波，因此，这里把一定动态范围内的匹配距离存储起来，以便于后续的误差回代。

下一步为误差回代，找到一个时移量的序列，使得两个序列相匹配后的误差和达到最小，找到最短路径。

式中：

L子波特征化后地震道样点数。

至此，两个序列中的子波动态匹配工作已经完成。对整个地震剖面，依次进行序列的两两匹配，把相互匹配的子波连接起来，就构成了同相轴，从而实现地震剖面的自动拾取。

3 数值实验

3.1 性能测试

首先，采用理论数据测试本文方法的正确性以及抗噪性；图1a 为理论地震剖面，构造一高斯子波，使其在横向上沿正弦函数分布，纵向上子波间的间隔为50 个网格，整个剖面共有11 个平行同相轴。图1b 是利用本文方法在图1a 所示数据上自动拾取的结果（图中红色圆点）。可以看出，图1b 中的层位均是由同一个子波构造而成，信噪比很高，非常易于拾取，当前大部分自动拾取方法都能做到。

图1 合成地震剖面及自动拾取结果覆盖于地震剖面上Fig.1 The synthetic seismic section and the automatic picking result covering on the seismic section

在图1a 的基础上，对其添加高斯白噪声，所使用的信噪比公式为

式中：

r(t)--原始有效信号；n(t)--高斯白噪声；T地震数据纵向采样点数。

按式（14）添加信噪比，如图2a 所示，加入随机噪音后的信噪比为5 dB。由于噪音的影响，剖面变得更加模糊，层位也不如图1a 明显，但是依然有明显地震同相轴的模式。采用本文自动拾取算法，选取振幅阈值为0.8（整个地震剖面最大振幅的0.8倍）、纵向动态匹配范围10 个采样点，作用于图2a的数据上，拾取结果见图2b，红色圆点表示拾取到的地震同相轴。拾取的同相轴总体上正弦函数变化的趋势不变，局部在纵向上随机抖动。这是由于噪音的影响使得子波的形状也变得上下抖动，波峰极值的位置发生偏移，拾取的结果才会出现图2b 中所示的情况。进一步降低信噪比，如图3a 所示为信噪比1 dB 的合成地震剖面，由于信噪比更低，剖面变得更加模糊，相同参数下自动拾取的地震同相轴的连续性自然不如图2b，但是依然能够看到明显的同相轴趋势。对比图1b、图2b 和图3b 中的拾取结果，可见本文自动拾取算法具有一定的抗噪能力。

对图1a 中高斯子波的位置（时间）和相位做随机扰动，最大扰动范围为8 个网格，扰动结果见图4a，相邻道的层位随机错动，相应的子波的相位也发生变化，这都给自动拾取带来了一定难度。在图4a 上采用本文方法进行自动拾取，选取振幅阈值和动态匹配范围与图2a 一致，拾取结果如图4b 中红色圆点所示，所有的层位都被拾取出来了，因为层位的错动拾取结果不如图1b 中那样光滑，层位线段表现的曲折不平，局部出现较大的抖动。在此基础上，加大最大扰动范围为16 个网格，对于复杂的实际数据来说，同一层位的相邻点错动距离也很少达到16 个网格，一般情况下很少出现这种剧烈的陡构造。扰动结果见图5a，相邻道的层位随机错动，子波的相位也会发生更大变化，使得自动拾取难度更大。同样参数设置下，使用本文算法拾取图5a中的地震同相轴，拾取结果如图5b 中红色圆点所示，可以看出，所有的层位都被拾取出来了，因为层位的错动更加严重，拾取结果不如图4b 中那样光滑，层位线段严重的曲折不平，局部位置出现较大的抖动，甚至出现了尖刺和沟谷。通过这两组实验说明，本文的方法对于子波位置和相位扰动具有一定的适应性，而且在实际数据中相邻点出现16 个网格的错动是比较罕见的，因此，即使在复杂情形下，本文拾取算法的精度仍然能够得到满足。

图2 5 dB 含噪地震剖面及自动拾取结果覆盖于地震剖面上Fig.2 The noisy seismic section(5 dB)and the automatic picking result covering on the seismic section

图3 1 dB 含噪地震剖面及自动拾取结果覆盖于地震剖面上Fig.3 The noisy seismic section(1 dB)and the automatic picking result covering on the seismic section

图4 8 个网格时移地震剖面及自动拾取结果覆盖于地震剖面上Fig.4 The time-lapse seismic profile(8 grids)and automatic picking results covering on the seismic section

图5 16 个网格时移地震剖面及自动拾取结果覆盖于地震剖面上Fig.5 The time-lapse seismic profile(16 grids)and automatic picking results covering on the seismic section

前文只考虑了子波位置和相位随机扰动和随机噪音的单一影响，将二者同时考虑进行拾取方法的稳定性测试也是很有必要的。在地震数据处理中，位置扰动对应着子波道间时差的变化，噪音大小则决定了数据的质量（信噪比和分辨率），二者是影响拾取结果好坏的关键因素。

在图1a 的基础上，对高斯子波的位置（时间）沿纵向做随机的扰动，如图6a 所示，最大的扰动范围为10 个网格，同时加入随机噪音，信噪比为5 dB。将本文的自动拾取算法应用在图6a 的合成地震剖面上，参数（振幅阈值和动态匹配范围）与前面保持一致，拾取结果如图6b 所示，红色圆点表示拾取到的地震同相轴。拾取的同相轴总体上保持正弦函数变化的趋势，局部存在纵向随机抖动，同相轴拾取结果有一定的道间时差，可见地震剖面在随机扰动的同时加上随机噪音的影响，使得层位曲线非常粗糙，波峰极值的位置发生偏移，拾取的结果才会出现图6b 中所示的情况。在此随机扰动的基础上，如图7a 所示，降低地震剖面的信噪比到1 dB。拾取结果展示在图7b 中，红色圆点表示拾取到的同相轴，相比图6b 结果，发现拾取结果同相轴连续降低，是因为随机噪音和位置扰动使得子波难以维持其原有的形状和变化趋势，子波能量分布不均衡，基于子波波形特征的自动拾取方法对子波变化比较敏感，如遇到子波变形过于剧烈的同相轴，则其拾取精度会下降。

图6 5 dB+10 个网格含噪+时移地震剖面及自动拾取结果覆盖于地震剖面上Fig.6 Noisy&time-lapse seismic section(5 dB+10 grids)and automatic picking results covering on the seismic section

图7 1 dB+10 个网格含噪+时移地震剖面及自动拾取结果覆盖于地震剖面上Fig.7 Noisy&time-lapse seismic section(1 dB+10 grids)and automatic picking results covering on the seismic section

3.2 实际资料测试

利用东部某探区实际地震成像剖面测试本文自动拾取方法的实用性。针对该探区地震剖面，选取振幅阈值为0.6、纵向动态匹配范围10 个采样点，图8a、图8b 分别为该探区Inline 测线方向和Crossline 测线方向的成像剖面以及本文算法自动拾取效果的展示（图中红色圆点所示）。从图中可以看出，通过给定合适的参数，文中拾取算法能够较为准确地识别出地震同相轴位置，为后续基于成像道集的层析成像、层位解释提供了一定帮助。

图8 东部某探区实际地震剖面同相轴自动拾取结果Fig.8 Automatic picking up result of in-phase axis of actual seismic section in an exploration area in the east

4 讨论

通过数值实验分析，针对低信噪比和道间时差地震数据，本文算法仍然可以实现地震同相轴的自动拾取工作。根据子波特征矢量的定义和算法原理，整个自动拾取过程中，振幅阈值和动态匹配范围两个参数需要注意。振幅阈值越大，子波特征表达后，保留能量越大的特征子波，因此，自动拾取结果越稀疏；而动态匹配范围，则是控制匹配追踪时纵向扫描的范围影响自动拾取的精度与效率，特别是如果引入种子点后，拾取特定地震同相轴时，动态匹配范围过大，可能会导致错误的拾取结果。

本文主要是给出了一种地震同相轴自动拾取的方法，并用于地震剖面同相轴的自动拾取中，其目的是为成像道集层析反演提供层位信息。而本文算法在地震处理解释中的具体应用理论上应该还有一定的扩展空间。比如，引入种子点概念后，可以用于特定反射界面的追踪、地震初至信息的获取、角道集RMO 拾取等等。但是，难以满足对特定断层的解释与识别工作。

5 结论

（1）基于褶积模型地震数据可以表达为反射系数序列和子波的褶积，地震同相轴拾取其实质是拾取主要反射系数序列的过程。

（2）子波以及噪音存在对地震同相轴的自动拾取带来很大的干扰，影响自动拾取的精度。基于子波形态特征表达，将地震数据转化为子波特征矢量数据，减弱了子波以及噪音对自动拾取的影响。通过定义矢量距离，利用动态规划法计算道与道之间的最小距离，从而实现的地震同相轴的自动拾取。同时通过特征化子波，可以对地震数据进行压缩，有利于提升海量地震数据自动拾取的计算效率。

（3）理论和实际数据的测试结果表明，本文算法可以适应低信噪比和存在道间时差的地震数据。但是，本文未具体讨论考虑种子点后算法的应用情况，这也是今后进一步研究需要探讨的问题。