我国地震勘探工作越来越多地涉及到复杂地表问题。一般而言, 复杂地表地区是指包含沙漠、戈壁、山地、黄土塬、水网等地貌, 以及风化层速度、厚度、地表高程等剧烈变化的地区^[1]。在复杂地表地区, 既可能存在高速地层、岩层直接出露, 又可能存在低速、疏松的冲洪积物、风积物堆积。对于地震勘探而言复杂地表问题的影响主要表现在近地表地震波速度剧烈变化, 强散射波、面波等干扰波发育, 引起地震波场、反射波同相轴畸变等^[2-3]。

在我国西部某复杂地表工区, 受崎岖地表和近地表多尺度非均匀体等因素的影响, 地震勘探采集的地震数据散射干扰波发育, 信噪比较低。地震散射波是当散射障碍体的尺度与入射波波长相比很小或相近时, 由于干涉效应而产生的与散射障碍体几何形状有关的相干或不相干的波列, 且往往以众多散射波互相叠加的复杂形态出现^[4]。散射干扰波在时间分布、频率特征等方面均与有效反射波组接近, 采用傅氏变换、小波变换^[5-6]等传统的时频域滤波方法难以达到理想的去噪效果。

目前公认的好的去噪方法是分析噪声与有效地震信号在不同"域"中的表现形式和特征差异, 进而采取最有针对性的去噪方法^[7]。本文在深入分析复杂地表散射干扰波和有效反射波信号特征差异的基础上, 针对散射干扰波与有效反射波在传播方向、出射角度方面不同的特点, 引进并改良了一种能够利用噪声与反射波信号的出射角度偏差, 从角度域入手进行滤波处理的矢量分解去噪方法^[8-9]。本文在简述方法原理和具体实现过程的基础上, 针对在算法实现过程中应采用能够尽可能提高窗口内地震数据相关度的"曲边时窗"模式等问题进行了深入论述, 并对比分析了包含随机干扰波、近地表散射干扰波的模拟地震数据去噪实验结果。最后通过实际单炮地震数据和地震叠加剖面成像效果验证了本文方法的有效性。

角度域矢量分解去噪方法的理论基础包括两个方面: ①合振动可看作若干个不同方向的振动共同作用的结果; ②相邻道地震信号具有相关性。利用相邻多道地震信号(图 1)统计反射地震信号出射时的方向信息, 保留与该方向一致的成分、压制偏离该方向的成分, 即可实现压制干扰噪声, 提高地震数据信噪比的目的^[10]。

图 1(Fig. 1)

图 1 某实际地震记录中的相邻地震道 a相邻3道; b相邻5道

矢量分解去噪法首先采用"地震信号矢量"概念来描述地震信号的方向特征^[11-12]。可以想象, 若将地震记录中相邻N道上某一时刻t的振幅值作为空间坐标系中不同坐标轴方向上的坐标值, 便可构建一个从原点出发的地震信号矢量, 其数学描述为:

其中, A_t为地震信号矢量, i₁, i₂, …, i_N为不同坐标轴方向上的单位矢量, x_t1, x_t2, …, x_tN为相邻N道地震信号在相应时刻的振幅值。从定义上可以发现, 地震信号矢量携带着地震波经地下传播后返回地表、被检波器接收时的出射角度信息。

若划分窗口进行统计计算, 则地震信号矢量的方向信息可进一步利用"单位相关矢量"来描述。计算某时窗单位相关矢量的过程如下: ①首先构建一个初始单位相关矢量d₀(如d₀=(i₁+i₂+…+i_N)/ $\sqrt N $), 根据地震信号矢量与d₀的内积是否小于0可将地震信号矢量划分为两个相互对立的集合; ②将某一集合内的所有地震信号矢量关于坐标原点进行翻转; ③将此时所有的地震信号矢量相加并归一化, 得到新的单位相关矢量d₁; ④用d₁代替d₀重复前述步骤进行迭代计算, 当d₁与d₀方向差别不大时终止迭代(通常在2~3次迭代之后可实现)。为了表述形象化, 现利用图 1a和图 1b中的地震数据将三维和五维情况下, 根据地震信号矢量计算单位相关矢量的过程通过三维直角坐标系和五维平行坐标系进行描述, 如图 2和图 3所示。更高维度的情况与它们类似。

图 2(Fig. 2)

图 2 三维地震信号单位相关矢量求取过程 a三维信号矢量划分; b三维信号矢量翻转; c矢量求和并归一化

图 3(Fig. 3)

图 3 五维单位相关矢量的求取过程 a五维信号矢量划分; b五维信号矢量翻转; c矢量求和并归一化

矢量分解去噪法中所涉及的方向和角度域概念与当前地震波逆时偏移成像研究中的角度域概念在本质上一致, 又略有不同。逆时偏移成像中的角度通常指重建后的震源波场和反向传播的检波点波场在地下空间某处相遇时地震波传播方向的夹角, 可以在地震波数值模拟过程中进行波场分离计算^[13-14]。而矢量分解去噪法中的角度概念则更倾向于表示地震波经过地下传播后出射地表时的角度信息, 可以直接根据地震记录计算, 通过地震信号矢量、单位相关矢量等概念进行描述。

执行前述步骤并逐点滑动时窗可计算出各时窗中心点处(时窗中间道的中间采样点位置)的单位相关矢量d_t=a_t1i₁+a_t2i₂+…+a_tNi_N, |d_t|=1, 然后便可通过矢量分解计算实现对地震数据的降噪处理, 具体步骤如下。

根据矢量投影和矢量加减法计算法则将A_t分解为沿d_t方向的相关成分A_st和垂直于d_t方向的非相关成分A_nt:

式中: |A_t|表示矢量A_t的模; 〈A_t, d_t〉表示矢量A_t与d_t之间的夹角; D=|A_t||d_t|cos〈A_t, d_t〉= $\mathop \sum \limits_{i = 1}^N {x_{{t_i}}}{a_{{t_i}}}$, 为矢量A_t与d_t的标量积(点积)。可以认为, 在相关成分A_st中有效地震信号占主要部分; 在非相关成分A_nt中干扰噪声占主要部分。

将A_st作为输出结果, 其各个分量对应时窗内各道地震数据t时刻的新的振幅值, 根据时窗逐点滑动计算的特点, 每次只更新时窗中心点处的计算结果。

通过前文对于矢量分解去噪法基本方法原理和实现过程的描述, 可知其在计算过程中并不破坏地震数据在横向和纵向上的相对振幅关系, 因此是一种保幅的去噪方法。分析可知, 提高求取单位相关矢量的精度是改善该方法去噪效果、防止损害有效地震信号成分的关键。

在复杂情况下地震信号的矢量方向受炮检距、地层倾角、干扰噪声等多种因素的影响, 在计算单位相关矢量之前最好选择最具相关性的地震数据。为了达到该目的, 目前较好的做法是结合倾角扫描技术^[15-16]。但这种做法需要在小的时窗范围内影响因素具有单一的方向特征, 这一前提通常在地震记录的很多位置上难以得到保证。

如果引入数学上曲边图形^[17]的概念构建一种"曲边时窗", 使时窗内各道地震数据均与时窗中间道达到最大相关, 无疑将提高求出的单位相关矢量的精度, 更重要的是这种做法可以淡化矢量分解法与倾角扫描技术相结合时较为严苛的假设前提条件。

为了阐述直观, 如果仅从时窗的角度进行分析, 经典矢量分解去噪法可以认为采用了"矩形时窗"模式; 与倾角扫描技术相结合的矢量分解去噪法采用了一种能够自动调整倾斜角度、确保时窗内各道数据达到最大互相关的"倾斜时窗"模式; 而本次研究采用的"曲边时窗"模式则能够进一步针对各道地震数据自动调整时移量, 确保窗口内各道均与时窗中间道达到最大相关。3种不同时窗模式的示意图如图 4所示。对添加了频带范围与有效地震信号主频接近的随机噪声的模拟地震数据, 分别采用矩形时窗、倾斜时窗、曲边时窗3种模式所得到的去噪结果如图 5所示, 滤出的噪声如图 6所示。由图 5和图 6可以看出, 对于同一套地震数据采用不同时窗模式的矢量分解去噪法进行处理时, 曲边时窗方法通常可以在高效去噪的同时使有效反射信号得到更好的保护。因此本次研究主要采用曲边时窗矢量分解去噪法进行复杂地表散射干扰波的去噪处理。

图 4(Fig. 4)

图 4 矢量分解方法不同时窗示意 a矩形时窗; b倾斜时窗; c曲边时窗

图 5(Fig. 5)

图 5 不同时窗模式去噪结果对比 a含噪模拟地震数据; b矩形时窗模式去噪结果; c倾斜时窗模式去噪结果; d曲边时窗模式去噪结果

图 6(Fig. 6)

图 6 不同时窗模式滤出的噪声对比 a加入的噪声数据; b矩形时窗模式滤出的噪声; c倾斜时窗模式滤出的噪声; d曲边时窗模式滤出的噪声

为了检测矢量分解去噪方法滤除近地表散射干扰波的效果, 参考前人研究成果^[18]建立了一个混杂大量非均匀地质体的水平层状介质模型, 如图 7所示。在该模型基础上采用交错网格有限差分法求解粘声波方程, 正演模拟得到了包含大量散射干扰噪声和少量全频带随机干扰噪声的模拟地震记录, 如图 8a所示。对该地震记录利用曲边时窗矢量分解去噪方法进行滤波处理, 得到的滤波结果和滤除的干扰噪声分别如图 8b和图 8c所示。可以看到, 图 8b中有效反射波的清晰程度和连续性均得到了改善, 原本存在于图 8a中的"斜"散射干扰噪声得到了较好的分离; 图 8c中散射干扰波特征明显, 不包含有效反射地震信号。据此可以认为, 本文方法能够有效滤除数值模拟数据近地表的散射干扰波。

图 7(Fig. 7)

图 7 混杂非均匀地质体的水平层状介质模型

图 8(Fig. 8)

图 8 近地表散射干扰矢量分解法滤波实验(曲边时窗) a含散射干扰波数据; b滤波后的数据; c滤除的噪声

在我国西部某地表形态复杂、近地表存在多尺度非均匀体的工区进行地下结构探测工作时^[19], 采集的地震数据中散射干扰波发育, 经过试验难以采用常规滤波方法达到理想的去噪效果。

工区关键采集参数为: 道距4 m、炮距24 m、接收道300、记录长度3 s、采样间隔1 ms, 滚动施工。区内典型的散射干扰波发育的野外地震记录如图 9a所示, 图 9b和图 9c分别为对地震数据进行角度域矢量分解去噪后得到的处理结果以及滤除的噪声。对比图 9a和图 9b可以看出, 原始地震记录上有效反射波周围的散射干扰波和随机噪声均得到了较好的滤除, 被保留下来的有效反射波同相轴在清晰度和连续性方面均得到了明显改善。在图 9c上, 被滤除的噪声散射特征明显, 经过分析认为其主要成分为与有效反射信号传播方向不同的散射干扰波以及传播方向无规律的随机干扰噪声。在图 9c上难以观察到有效反射波同相轴, 结合前文对矢量分解去噪算法原理的分析, 可以认为在去噪过程中有效反射波信号得到了合理的保护。

图 9(Fig. 9)

图 9 滤除散射干扰波前、后的单炮地震记录以及滤除的噪声 a滤除散射干扰波前的地震记录; b滤除散射干扰波后的地震记录; c滤除的噪声

在其它处理步骤和参数一致的情况下对L2线共704炮滤除散射干扰波前、后的地震数据进行叠加处理, 得到的反射波叠加成像剖面如图 10所示。对比图 10a和图 10b可以看出, 滤除近地表散射干扰波之后的地震剖面, 质量有很大提升, 尤其是在图中红色箭头处原本受散射干扰波严重影响的近地表成像区域, 地震波同相轴的连续性、可追踪性、分辨率均得到了明显改善。这种改善有利于连续追踪地震目标岩体的边界和顶界。

图 10(Fig. 10)

图 10 L2线滤除散射干扰波前(a)、后(b)的叠加成像剖面

本次研究提出了将角度域矢量分解去噪方法应用于复杂地表地震勘探中近地表散射干扰波压制的方法, 通过理论分析以及模拟地震数据、实际叠前地震数据和叠加成像剖面对比, 证实了角度域矢量分解去噪方法能够针对近地表散射干扰波实现良好的信噪分离效果。

作为一种能够利用噪声与反射地震信号出射地表时的角度偏差, 从角度域入手进行去噪的滤波方法, 当噪声信号与有效信号在分布和时频特征等方面接近时, 可以作为时频域、小波域等去噪方法的有力补充。这种方法还具有保持振幅以及可用于叠前地震数据处理等优点, 分析研究过程和研究成果可以发现, 该方法能够较好地从实际地震数据中分离出散射地震波, 据此推测该方法可用于地震散射波成像等方面的研究工作。