一种针对气枪记录的三阶段初至震相拾取方案

朱逸馨 张云鹏

1.北京大学地球与空间科学学院地球物理学系, 北京 100871; 2.中国地震局地球物理研究所地震观测与地球物理成像重点实验室, 北京 100081; † E-mail： zhuyixin@pku.edu.cn

随着主动源相关技术的发展, 21世纪以来, 人们开始使用主动源来探查和解决地震领域的问题。虽然要解决全球尺度下的问题存在很大的困难, 但在局部地区的研究中卓有成效[1-2]。作为主动源震源之一, 气枪早在1963年就已问世, 但是直至专利失效后, 才在地震勘探中广泛使用[3]。此后, 因低频信息丰富、重复性高、激发能量大、与水耦合性好、衰减慢以及环境友好、性价比高等优点, 气枪震源成为主动源中最重要的方法之一[4]。

2015年10月, 以探测长江地下深部结构为目的,中国地震局组织在长江安徽段进行激发的气枪震源探测实验, 即“安徽实验”[5]。此次气枪震源探测实验台站密度较高, 覆盖范围较大, 若能合理处理该次实验数据, 既可以对气枪震源在长江区域使用的效果进行评估[6], 又有利于陆上气枪震源探测地下深部结构等研究[7]。

地震资料分析中, 初至拾取是至关重要的一项基础工作。初至拾取的准确性直接决定静校正、反演以及层析成像等后续处理的质量[8]。最早的初至拾取方法是20世纪70年代由Peraldi等[9]提出的通过对相邻道进行互相关来确定相对初至时差的方法。Stevenson[10]1976年提出STA/LTA方法, 利用有效信号与噪声的振幅差异, 通过两个长度不同的滑动时窗来判断初至到达时间, 该方法目前广泛应用于以微地震为代表的初至拾取中[11]。Coppens[12]1985年根据有效信号与噪声的能量差异, 提出通过不同长度时窗内的能量比来确定初至的方法, 至今仍被广泛使用并不断改进[13]。总之, 初至拾取的结果在很大程度上取决于数据的信噪比。

本文选用的测线为“安徽实验”中沿长江布设的短周期流动测线, 位于安徽省南部, 靠近江苏、浙江等人口密度较大的省份, 测线附近大小不同的城市、村庄密集, 工业发达。此次数据由埋深较浅的流动台站记录, 背景噪声对数据的影响较大, 信噪比普遍较低[5], 初至拾取十分困难。针对这个问题,本文设计一种针对气枪记录的三阶段初至震相拾取方法： 首先根据测线排列和采集方式与地震勘探测线类似的特点, 引入地震勘探的信号处理方法, 对数据进行预处理; 然后, 分别以平均能量比、熵值和分形维度作为特征值, 计算单道数据的特征曲线;最后, 针对难以设定阈值的问题, 结合保边平滑方法[14], 对特征曲线进行平滑, 根据导数最值位置获得单道数据的初至时间。本研究还基于Sabbione等[15]提出的选择修正步骤, 根据数据的特点, 设计可自动拾取和判别有效初至的流程, 对“安徽实验”的主动源数据进行自动拾取。与传统方法相比, 本文采用的自动拾取流程对数据信噪比要求更低, 并可为后续处理提供更多可用资料。

1 气枪源数据分析

“安徽实验”的激发点和流动测线分布见图 1, 20个固定激发点(五角星)由原用于海上的气枪震源船激发[7]。本文使用的流动测线(三角形)是沿长江布设的L0测线记录到的叠后数据。L0测线的直线长度约为253 km, 包含97个流动台站。地震仪的采样率为200 Hz, 最小偏移距为0.576 km, 最大偏移距为252.208 km。

为了更好地分辨噪声与信号, 以便有针对性地压制噪声, 首先要观察并掌握有效信号和主要噪声的特征。本文采用图 2 所示的作图方法, 同时对数据的时域波形、频谱和时频谱进行观察和分析, 给出时频分布图。另外, 本文用“D07.L0.133”的数据形式来表示当D07点激发时, L0测线上133号台站记录的数据。

根据气枪震源的实际频率特征[16]和震源台站记录特征, 可以确定气枪震源激发的信号主要频率范围为4～6 Hz。图2(a)展示“安徽实验”中的一道原始数据, 其时域图(右下图)中 0 s 附近出现的波形为所要拾取的初至波, 本文视其为有效信号, 在时频谱(右上图)中表现为较为集中的能量团, 对比频谱(左图)可以发现有效信号的主频为5 Hz 左右。根据有效信号的频率特点, 观察其他数据, 我们将数据中存在的噪声分为单频噪声、低频噪声和高频噪声, 图2(b)～(d)为原始数据中噪声特征明显的示例。

图2(b)显示一道包含单频噪声的数据。单频噪声的特点是频带极窄, 能量很强, 一般是其他信号的数倍甚至数十倍, 对信号影响极大, 但是可以得到其在频率轴上的具体位置。

图2(c)显示一道包含低频噪声的数据。由于要拾取的初至波主频为5 Hz左右, 本文中低频噪声指频率小于3 Hz的信号。

图2(d)显示一道包含高频噪声的数据。本文中高频噪声指频率大于6 Hz的信号。由于低频和高频噪声与有效信号的差别主要是频带范围不同, 因此本文直接对低频和高频噪声进行处理。

对比不同道的数据可以发现, 共接收点的数据频谱有一定的相似性。以189号台站记录的部分数据为例, 虽然激发点不同, 但都在4 Hz和6 Hz附近存在明显的单频噪声, 且低频与高频噪声的频率分布及强度几乎一致(图3)。

2 数据预处理方法

本文使用数据的信噪比普遍不高, 若直接进行初至拾取, 则难度较大, 结果的可靠性较差, 不利于后续工作。因此, 本文在进行初至拾取前, 先对数据进行噪声压制预处理, 在一定程度上提高数据的信噪比, 减小噪声对初至的影响。具体地, 由于部分单频噪声的强度是其他信号的数倍甚至数十倍, 严重影响数据的信噪比, 本文采用的流程是先压制单频噪声, 后压制低频(高频)噪声。

图1 气枪震源实验中激发点和流动测线分布Fig.1 Source and receiver’s distribution in air-gun source experiments

2.1 压制单频噪声

单频噪声是一种连续干扰, 主要特点是频带极窄, 能量很强, 对信号的影响极大, 因此在初至拾取前必须进行单频噪声的压制。我们根据单频噪声频带极窄和可以得到具体频率的特点, 利用中值滤波可以滤除尖脉冲[17]的性能, 在频率域中对单频噪声的频带进行中值滤波, 从而压制单频噪声。

中值滤波的基本原理是, 通过将一个点的数值替换为邻域中各个点的中值来消除孤立的噪声点。以图 4 中数据为例, 图4(a)显示原始数据的波形, 其频谱如图4(d)所示。在图4(a)和(d)中几乎无法分辨出有效信号的波形和频率, 可以发现在3～4 Hz的频率范围内存在3个明显的单频噪声。首先, 根据单频噪声的频率的振幅值是其他信号的数倍甚至数十倍的特点, 通过在频率域查找最大值的方法, 确定单频噪声的频带。然后, 依次以频带中每个点为中心, 在频率轴上取一个有21个采样点的频率窗, 用窗内频率振幅值的中值替换原频率的振幅值, 即对单频噪声的频带进行中值滤波。对频带外的频率不进行任何处理, 得到图4(e)所示的频谱, 单频噪声的频率的振幅值经过中值滤波得到平滑, 但还存在很强的低频信号。将处理后的频谱反变换回时域后, 得到的波形如图4(b)所示, 信噪比仍然较低, 根据图4(e)所示的频谱可知, 信号中还存在较强的低频信号, 所以在此基础上还要对低频和高频噪声进行压制。

2.2 压制低频和高频噪声

2.2.1 小波变换基本理论

针对低频和高频噪声, 本文进行基于小波变换的带通滤波。小波变换是20世纪80年代发展起来的一种具有较好时频局部化特性的方法[18], 由于可以同时在时间域和频率域进行分析而成为去除信号噪声的有力工具之一[19-20]。对离散信号d(n)进行小波分解的过程可以理解为用小波函数Ψ(n)和相应的尺度函数φ(n), 通过线性组合的方式来展开表达[18]：

图2 原始数据的时频分布Fig.2 Time-frequency distribution of data

图3 189号台站接收的不同激发点数据的频谱Fig.3 Frequency spectrum of data recorded by station 189

其中,k为平移因子;j为伸缩因子, 对应第j层的小波分解。式(1)中的尺度系数a和小波系数b可以分别通过信号与小波函数和尺度函数的内积获得：

若期望将离散信号d(n)按式(1)的方式分解为j层,可以通过式(2)和(3)获得各层的小波系数以及尺度系数。

可以这样理解, 小波变换在时间域将信号分解为高频和低频两个部分： 高频部分表示为小波系数;低频部分表示为尺度系数, 并作为下一次分解的输入信号。最终, 可以将数据表示为多个小波系数和一个尺度系数。将这些小波系数和尺度系数重新组合后, 可重构被分解的数据。因此, 小波变换本质上是一个多通道的带通滤波器, 可以压制与有效信号频带分离的噪声, 但无法压制与有效信号频带重叠的噪声(如白噪声、单频噪声等)。本文利用小波变换的分解和重构过程, 通过对数据进行带通滤波来压制与有效信号频带分离的低频和高频噪声。

2.2.2 基于小波变换的带通滤波

本文选取db15小波作为小波基。dbN小波表示Daubechies小波, 是由Daubechies等[20]以2的整次幂为尺度构造的一种小波函数。dbN小波存在离散小波变换, 且具有正交性、双正交性以及紧支撑性等特征, 其正则性随阶数N的增加而增加。

仍以图4中数据为例, 对压制单频噪声后的数据进行基于小波变换的带通滤波。滤波前的数据波形如图4(b)所示, 本文基于小波变换对数据进行带通滤波的主要步骤如下。

图4 预处理结果Fig.4 Result of pre-treatment

1)小波分解。利用式(2)和(3), 通过小波变换将数据分解为5层, 获得的小波系数和尺度系数如图5所示。其中, 第j层小波系数对应对数据的第j次小波分解, 且所对应的中心频率大致是上一层的1/2, 而尺度系数相当于分解后剩余的低频部分。

2)对系数进行处理。选择与有效信号中心频率相近的小波系数, 其余系数置 0, 本文选择的是第3～5层的小波系数。

3)小波重构。将处理后的系数重新组合构建成信号, 获得基于小波变换进行带通滤波后的波形,如图4(c)所示, 背景噪声中的低频信号已被有效地滤除。

3 初至拾取方法原理

3.1 平均能量比法

平均能量比法的基本思想是Coppens[12]提出的能量法, 主要是利用有效信号与背景噪声的能量差异[11,13]。本文设置两个长度不同的时窗, 利用式(4)和(5)分别计算两个时窗内信号的平均能量, 然后通过式(6)获得平均能量比。

其中,Es(n)表示长度为ns的时窗内信号的平均能量,Et(n)表示长度为n的时窗内信号的平均能量,e(n)表示Es(n)与Et(n)的平均能量比。

如图6所示, 本文将长时窗取为第一个数据到短时窗中最后一个数据的长度, 即其长度随短时窗的平移而增大, 因此需要设置的参数有且仅有短时窗的长度ns。

图5 小波分解后的尺度系数与小波系数Fig.5 Scale coefficient and wavelet coefficients after wavelet decomposition

图6 平均能量比法原理示意图Fig.6 Schematic diagram of average energy ratio method

3.2 熵值法

熵的概念最早出现在热力学中, Denis等[21]在一维信号的时间序列中引入这个概念, 此时, 熵可以理解为信号d(n)波动的强度, 计算公式如下：

其中,Ld(n)表示长度为ne的时窗中波形曲线的长度,则式(7)可以写成

3.3 分形法

分形的概念最早由Mandelbrot等[22]提出, 用来描述物体的不规则程度。分形维度是用来描述这种不规则程度的重要参数, 其中平面曲线的分形维度在[1, 2]范围内。Boschetti[23]将这个概念引入初至拾取方法中, 可以理解为信号的波形越杂乱无章,分形维度就越大, 否则反之。值得注意的是, 在有效信号到达时, 信号的分形维度会突然降低, 这与前面介绍的平均能量比法和熵值法不同。

根据定义, 地震道集是一种自仿射分形[24], 且满足函数：

变形后可以写成

可以将H视为logF(k)/logk的斜率。F(k)表示信号d(n)经过间距为k的采样后, 其平方差的期望值[15,25], 即

同时, 自仿射分形的分形维数D满足

首先根据式(11), 分别计算当采样间隔k=1, 2,3, 4时, 期望值F(k)的取值; 然后利用直线拟合得方法获取logF(k)/logk的斜率H; 根据H值和式(12),很容易得到分形维度D。

需要说明的是, 高斯白噪声具有随机性, 复杂度较高, 理论上分形维度为 2, 有效信号是由激发子波与地下介质通过褶积得到的, 与高斯白噪声相比维度较低。实际计算分形维度时, 背景噪声并非全部是随机噪声, 背景噪声和有效信号的维度差异并不大, 目前常用的做法是先在数据中添加与背景噪声能量相近的低振幅高斯白噪声[15,25], 一方面可以破坏相干噪声的相干性, 另一方面可以使分形维度的变化更明显。因此, 本文在计算分形维度前, 先利用式(13)估计数据的信噪比SNR：

其中,d为信号,dn为背景噪声。然后利用式(14)得到高斯白噪声d′：

其中, random为随机数。最后, 通过式(15)将低振幅的高斯白噪声d′添加到原数据中：

4 保边平滑方法

设置阈值是目前判定初至时间的常用方法[26],即先确定一个阈值, 将每道数据的特征值到达并超过该阈值的瞬间作为初至时间。值得注意的是, 理论上, 有效信号的分形维度小于背景噪声, 因此当特征值为分形维度时, 我们将特征值到达并小于该阈值的瞬间作为初至时间。以图7(a)中的单道数据为例, 以平均能量比、熵值和分形维度作为特征值,得到的特征曲线分别如图7(b)～(d)中黑色曲线所示, 通过阈值法拾取的初至时间如图7(a)中倒三角形所示。

观察图7可以发现, 特征曲线的变化不是突变的, 因此很难确定最佳阈值。另外, 若有效信号到达前存在特征值小于有效信号, 但大于阈值的噪声(特征值为分形维度时, 则情况相反), 则会严重影响初至时间的拾取。

保边平滑是一种在保护边缘信息的基础上, 对图像进行平滑的方法, 常用来提高地震剖面的信噪比[14]。本文使用简单的一维保边平滑方法, 通过设置长度为np个采样点的滑动时窗, 对特征曲线进行平滑, 以期获取特征曲线的主要变化趋势。以np=3为例, 简单来说, 对数据d(i)进行保边平滑就是先获取包含数据d(i)的所有时窗。

然后, 在所有时窗中找到标准差最小的时窗, 用该时窗数据的平均值替换原数据的d(i)。

对图7(a)中单道数据的特征曲线进行保边平滑后, 曲线的变化更加“干脆”(图7(b)～(d)中灰色曲线)。对平滑后的曲线进行求导, 结果如图7(f)～(h)所示。当特征值为平均能量比或熵值时, 我们选取图7(f)～(h)中导数最大值的时间作为初至时间。同理, 当特征值为分形维度时, 我们选取图7(h)中导数最小值的时间作为初至时间, 图7(e)中倒三角标注的是基于3种不同特征曲线获得的初至时间, 与图7(a)中利用阈值法获得的初至时间完全相同。

图7 单道数据的特征曲线及初至拾取结果Fig.7 Characteristic curve and result of pick

利用保边平滑方法进行初至拾取, 可以避免对阈值大小的讨论, 有利于多道数据的自动拾取, 且只要背景噪声的特征值小于有效信号(分形维度大于有效信号), 便不会发生误拾情况, 可以有效地避免因背景噪声特征值较大(或分形维度较小)而发生的误拾, 从而降低拾取结果对数据信噪比的依赖。

为了验证利用保边平滑方法的有效性和准确性, 我们选择同样为气枪震源的炮集进行初至拾取,该数据来自Yilmaz等[27]提供的40个炮集中的记录17。拾取结果如图8所示, 其中蓝色实线为阈值法拾取的初至时间, 红色虚线为利用保边平滑法拾取的初至时间。该数据为勘探数据, 信噪比高且初至能量强, 阈值法和保边平滑法拾取的初至结果十分接近。

5 数据实例

我们通过以下3个阶段, 对“安徽实验”中L0测线记录的叠后数据进行初至拾取。

阶段1： 对原始数据进行预处理, 包括利用中值滤波压制单频噪声以及基于小波变换压制低频和高频噪声。

阶段2： 分别以平均能量比、熵值和分形维度为特征值, 计算数据预处理后的特征曲线。

阶段3： 结合保边平滑方法, 对特征曲线进行平滑, 使其变化更“干脆”; 根据保边平滑结果的导数最值位置来确定初至时间, 当特征值为平均能量比或熵值时, 初至时间取导数最大值时间, 当特征值为分形维度时, 初至时间取导数最小值时间。

在实际初至拾取过程中, 由于初至波能量较小,或者有效信号完全被噪声淹没等原因, 会导致自动拾取结果出现错误。Sabbione等[15]在处理勘探数据时, 对炮集数据进行了综合考虑, 最终通过对拾取结果进行直线拟合, 将所得的直线作为参考初至,对拾取结果进行一定的约束和判断。由于错误的初至时间会严重影响后续处理, 基于“宁缺毋滥”的原则, 我们利用直线拟合来去除误拾的初至。通过三阶段拾取初至后, 对近偏移距的离散度较小的初至进行直线拟合, 并利用拟合直线对初至进行判断,将与直线偏离太大的近偏移距初至及离散度较大的远偏移距初至都视为无效初至。

图8 阈值法(蓝)和保边平滑法(红)拾取的结果Fig.8 Result of threshold method (blue)and the edge-preserving smoothing method (red)

然后, 分别使用传统的设置阈值方法(简称传统方法)与本文方法, 对同一数据进行初至拾取, 比较两种方法的拾取结果。以D16激发点数据为例,图9中蓝色圆点显示在没有对数据进行预处理的情况下, 使用传统方法拾取的初至时间。

图9中红色圆点显示利用本文三阶段方法进行初至拾取, 并利用曲线拟合进行选择修正后, 在拟合曲线附近的初至, 红色“×”号显示不在拟合曲线附近的初至。为了更明显地显示传统方法与本文方法所得初至的差异, 将21, 23和25道集的拾取结果放大展示于图的右侧。经过综合考虑和多次尝试后, 本文中不同时窗的长度取值如表 1 所示。其中,nT表示有效信号第一个周期的采样点数, 根据多个数据的观察, 有效信号第一个周期大致为0.2 s, 为了便于计算和对比,nT取 40 个采样点数。

观察图9可以发现, 使用传统方法拾取的结果中存在很多错误, 例如图9(a)中第47道、图9(b)中第45～53道以及图9(c)中第45和47道等近偏移距数据, 这些错误结果是由初至到达前存在特征值小于初至但大于阈值的噪声造成的。同时, 图9(a)中第31和35道数据的拾取结果比实际初至延后, 原因是该道初至波能量较弱, 但后至波能量较强。这些出现错误结果的道, 使用本文方法都可以拾取到有效初至。但是, 对于远偏移距的道集, 本文方法同样无法获得有效初至。

如表 2 所示, 对 L0 测线的所有数据进行处理后, 使用本文方法拾取的有效初至数量(近偏移距的、初至离散度较小的范围内)远多于传统方法。比较 3 种特征值的结果可以发现, 本文方法对于分形法的改善效果远不如另外两种特征值方法。分析其原因主要是, 信号的分形维度仅在[1, 2]范围内变化, 传统的阈值法便可以较好地拾取大部分数据的初至时间, 而平均能量比和熵值的变化范围较大,并且, 不同道集间的变化范围差异也较大, 因而在自动拾取中, 利用保边平滑法可以更好地适应多道数据。

图9 传统方法(蓝)与本文方法(有效初至(红点)、无效初至(红叉))的拾取结果Fig.9 Result of traditional method (blue)and suggested method (red)

表1 不同时窗的长度取值Table 1 The width of different time windows

表2 本文方法及与传统方法拾取有效初至数量Table 2 Results of the traditional and proposed method

对比 3 种不同特征值的拾取方法, 由于分形维度反映的是波形复杂度, 基于分形维度的拾取方法在初至波能量较弱时表现最好, 拾取的有效初至在3种方法中最准确, 但背景噪声中存在相干噪声时对基于分形维度的拾取结果影响最大。基于平均能量比的特征曲线受相干噪声的影响较小, 拾取的有效初至数量最多, 但在初至波能量较弱的情况下,准确性较差, 容易发生延时情况。基于熵值的特征曲线变化最平缓, 因此根据熵值确定的初至时间最不稳定, 经常比实际初至提前或延后。

6 结论

本文提出利用以下三阶段方法对“安徽实验”中L0测线记录的叠后数据进行初至拾取。阶段1： 对原始数据进行预处理, 包括利用中值滤波压制单频噪声以及基于小波变换压制低频和高频噪声; 阶段2： 计算数据预处理后的特征曲线; 阶段3： 结合保边平滑方法, 对特征曲线进行平滑, 根据保边平滑结果的导数最值位置来确定初至时间。在对实际数据进行初至拾取时, 我们还利用曲线拟合对拾取结果进行选择修正。本文方法在近偏移距的初至拾取中获得较好的结果, 但无法在远偏移距的数据中获得有效初至。

本文还对比了3种常用特征值, 即平均能量比、熵值和分形维度对拾取结果的影响。其中, 基于分形维度的方法在初至波能量较弱时的拾取结果最准确, 但容易受到相干噪声的影响; 基于平均能量比的特征曲线受到噪声的影响较小, 拾取的有效初至数量最多, 但准确性较差, 容易在初至波能量较弱的情况下发生延迟; 基于熵值的初至时间最不稳定,经常比实际初至提前或延后。

目前, 不少自动拾取方法都比较成熟。但是,针对用于勘探数据的讨论较多, 对研究地球深部问题数据的讨论较少。本文利用“安徽实验”数据, 对以探测长江地下深部构造为目的的气枪震源数据进行初至拾取的初步尝试, 可为今后的相关研究提供可借鉴的经验。

致谢 北京大学宁杰远教授对本研究给予很大的帮助, 并提出建设性的意见; 中国地震局地球物理研究所王宝善研究员提供有关气枪震源实验的关键数据。在此一并致以诚挚的谢意。