面波是反射波地震勘探中常见的噪声干扰, 在炮点附近呈扇形分布, 视速度低、能量强, 且具有频散的特点, 严重影响了中深层的有效反射, 降低了地震资料的信噪比。因此, 有效的面波压制对后续地震资料的处理和解释具有重要作用。目前有多种面波压制方法, 如一维频率域滤波法、f-k滤波法、Hilbert-Huang变换法、小波变换法、曲波变换法等。

一维频率域滤波法将每一道地震记录进行傅里叶变换后通过设计合适的高通滤波器来滤除面波成分^[1], 但高通滤波同时会丢失低频部分的反射波信号, 产生吉布斯现象; f-k滤波法考虑了面波频率低、视速度低的特点^[2], 但当面波能量远强于反射波时, 会造成有效波严重失真, 并在滤波器的边界产生新的相干噪声; 孔庆丰^[3]利用Hilbert-Huang变换较高的时频分辨率来压制面波, 但对于有效波与面波重叠较多的区域不能很好地识别区分; 张华等^[4]、岳龙等^[5]通过小波变换的多尺度特性识别出面波区域, 减小了去除面波时对有效信号的影响, 但小波变换只对处理点奇异有效, 对于地震记录的曲线奇异特征不能有效地表示; 曲波变换对线和超平面的奇异性表示优于小波变换^[6], ZHENG等^[7]、董烈乾等^[8]对曲波域中主要含面波的尺度和方向的系数进行阈值处理, 也能够有效去除面波干扰。

实际地震数据通常由多种信号成分组成, 只采用单一变换不能有效地表现数据的内部结构特征。近年来稀疏表示理论迅速发展, 并被不断引入到地震数据处理领域^[9-11]。STARCK等^[12-13]提出了以信号的稀疏性和形态多样性为基础的形态成分分析(Morphological Component Analysis, MCA)方法, WANG等^[14]将其引入地震记录面波压制处理, 提出了基于单道处理的压制面波的改进MCA法。陈文超等^[15]和李海山等^[16]根据反射波与面波形态结构的差异, 将一维局部离散余弦变换与一维小波变换分别作为反射波和面波的稀疏表示字典, 通过求解稀疏优化问题实现了两者的分离。以上方法都是以单道记录为处理对象, 没有充分考虑到道与道之间的相关性, 增加了数据的计算复杂度。XU等^[17]对整个二维地震记录进行处理, 联合二维非抽样离散小波变换(2D-Undecimated Discrete Wavelet Transform, 2D-UDWT)和二维局部离散余弦变换(2D-Local Discrete Cosine Transform, 2D-LDCT)构成超完备字典, 将其应用到MCA方法分离反射波和面波分量的问题中, 也得到了较好的压制效果。

鉴于小波变换自身的局限性, 本文在XU等^[17]研究工作的基础上, 将表示性能更好的离散曲波变换作为面波的稀疏表示字典, 将2D-LDCT作为反射波的稀疏表示字典, 并对2D-LDCT字典的内部结构进行优化, 更好地分离出了面波分量和反射波分量, 实现了对面波的有效压制。最后用理论模型数据和实际地震资料进行了应用测试。

稀疏表示能够在变换域内用尽可能少的原子的线性组合来逼近原信号, 对于信号x=[x₁, x₂, …, x_N]^T, 若只有m个元素为非0值, 绝大部分元素为零, 则称该信号是稀疏的, 其稀疏度为m^[18]。如果信号在时域不满足稀疏性特征, 可以通过适当的变换使其稀疏化。一般地, 一个信号x可由K个原子d_k的线性叠加表示:

(1)

其中, 字典D是N×K矩阵, 若D的列数大于行数, 即K>N, 称此字典为超完备字典; α_k(k=1, 2, …, K)为x在字典D=[d₁, d₂, …, d_K]下的稀疏表示系数, α的稀疏度可由α的0范数‖α‖₀表示, 即α中非0分量的个数。若α是稀疏的, 则信号x可在字典D下得到稀疏表示, 即该信号可表示为一小部分原子的线性叠加。

MCA方法将多种具有不同原子特征的字典联合起来构成超完备冗余字典, 从而获得更加稀疏的信号表示方式, 实现对不同形态信号分量的有效分离^[12-13]。假定输入的待处理信号s∈R^N×1由N个形态特征不同的分量s_n线性组成, 即:

(2)

式中:s_n(n=1, 2, …, N)表示不同形态成分的信号分量; 超完备字典T_n(n=1, 2, …, N)为s_n对应的稀疏表示字典; α_n(n=1, 2, …, N)为s_n在相应字典下的稀疏表示系数, 同时MCA理论假定对于混合信号中的每一种信号分量s_n, 其在相应的字典T_n下都能得到非常稀疏的表示, 而在字典T_k≠n下不能得到有效的稀疏表示。混合信号s的稀疏表示可以描述为求解下列问题:

(3)

按照基追踪(Basis Pursuit, BP)算法和匹配追踪(Matching Pursuit, MP)算法^[19]的思想, 对(3) 式中的约束条件适当松弛, 可转化为求解如下问题:

(4)

式中:‖α_n‖₁表示第n个分量稀疏表示系数的L₁范数稀疏惩罚项; 表示L₂范数保真约束项; λ为拉格朗日乘子, 控制着字典合成信号与原信号s相比的失真度。(4) 式是一个典型的凸优化问题, ELAD等^[20]将块协调松弛(Block Coordinate Relaxation, BCR)算法用于MCA问题的求解, 成功地分离出了不同形态的信号分量。

假设含面波地震记录的每道信号是有效反射波和面波的线性叠加, 则单炮地震记录可表示为:

(5)

式中:s_r, s_g, n分别表示有效反射波分量、面波分量和随机噪声。在该模型中加入随机噪声分量为0均值的高斯白噪声, 可以提高模型的普适性。

对于反射波分量s_r, 用超完备字典T_r∈M^N×L_r进行稀疏表示时, 有:

(6)

得到的表示系数α_r非常稀疏, 字典T_r能够有效地稀疏表示反射波分量s_r。

对于面波分量s_g, 用超完备字典T_g∈M^N×L_g进行稀疏表示时, 有:

(7)

得到的表示系数α_g也非常稀疏, 则字典T_g能够有效地稀疏表示面波分量s_g，如果字典T_r不能很好地稀疏表示面波分量, 字典T_g不能很好地稀疏表示反射波分量，则单炮地震记录s可利用由T_r和T_g联合组成的超完备字典进行稀疏表示, 并通过求解如下的最优化问题实现面波与反射波的分离:

(8)

MCA算法的关键环节是字典组合的选择。根据MCA理论, 2种字典不仅需要能够分别稀疏地表示各信号分量, 它们的原子波形也需要具有明显的差异性。由于反射波和面波的形态结构差别很大, 因此只要为这2种信号分量找到相应的稀疏表示字典, 就能够在相应的稀疏域中对其进行有效地稀疏表示。

地震记录中的面波分量频率低、视速度低, 具有频散特性, 且面波与反射波同相轴斜率差异较大, 而离散曲波变换具有非常好的局部性、各向异性和很强的方向性, 能够有效地表示波前和具有线状特征的同相轴^[21], 因而我们选择离散曲波变换作为面波的稀疏表示字典。

离散曲波变换^{[6, 22]}是一种多尺度变换方法, 可以有效地表示二维信号。在二维笛卡尔坐标系中, 离散曲波变换采用同中心的楔形窗将频率域划分成不同的尺度和角度, 代表频率域内不同尺度和方向的曲波，该楔形窗定义为:

(9)

式中:分别表示径向窗函数和角度窗函数。引入一系列相同间隔的斜率tanθ_l=l·2^－[j/2], l=－2^[j/2], …, 2^[j/2]－1, 则楔形窗可进一步表示为:

(10)

式中:S_θl是剪切矩阵, S_θl=10－tanθ_l1。则离散曲波系数可表示为:

(11)

本文采用基于快速wrapping算法的离散曲波变换作为面波的稀疏表示字典, 对于反射波频带和面波频带交叠的区域, 离散曲波变换能够充分利用良好的方向识别能力, 有效地分离出面波分量。

地震记录反射波分量局部相关性强, 具有明显的波动特征, 而局部离散余弦变换(LDCT)是类似于傅里叶变换的一种实值变换, 它对具有高度相关性的信号有非常好的能量聚集性, 能够有效地捕捉信号的局部结构特征, 准确地描述不同时段的信号特性, 自适应地跟踪信号的变化^[23-24], 因而我们选择2D-LDCT作为反射波分量的稀疏表示字典。

LDCT是一种基于重叠正交基或双正交基的变换, 本文采用的是Ⅳ型局部余弦基, 其定义如下:

(12)

其中, h_i=a_i+1－a_i表示信号的第i段, g_i(t)是定义在区间I_i=[a_i－η_i, a_i+1+η_i+1]上相互重叠的光滑窗函数, 它将信号分成分段信号, 且满足a_i+η_i≤a_i+1－η_i+1, η_i和η_i+1分别为左右两侧的重叠半径。LDCT一般通过先对信号做折叠分段处理, 再对其进行快速Ⅳ型LDCT得到。原始的2D-LDCT字典只能处理方阵数据, 本文对其进行了优化, 无需在数据变换时将其扩充为方阵, 使其能自适应地处理任意大小的地震数据, 从而实现对反射波有效的描述和稀疏表示。

我们利用两组实验来对比说明不同字典对反射波和面波信号的稀疏表示性能。实验选取经过滤波后不含面波的实际地震记录作为实验反射波信号, 该记录共240道, 每道采样点数为1024, 如图 1所示。选取受到强面波干扰的实际地震记录切片作为实验面波信号, 该记录共64道, 每道采样点数为1024, 如图 2所示。选取的2D-LDCT字典为离散Ⅳ型余弦变换, 窗口宽度为10;选取的离散曲波变换字典为基于wrapping算法的快速离散曲波变换。为了更好说明离散曲波变换字典和2D-LDCT字典对两种信号分量稀疏表示的优越性, 并与2D-UDWT字典稀疏表示的性能进行对比, 其小波母函数采用消失矩为10的Symmlet小波。对2组实验信号分别用这3种字典进行稀疏表示, 将得到的系数根据幅值进行排序, 并按相同的比例系数将其中较小幅值的系数置0, 再分别通过3种变换的逆变换求出各自的重构信号, 计算重构信号与原始信号的重构相对误差, 以此来评价各字典对反射波和面波的稀疏表示能力, 重构误差越小, 则稀疏表示能力越强。

图 1(Fig. 1)

图 1 实验反射波信号

图 2(Fig. 2)

图 2 实验面波信号

反射波信号在各变换字典下的重构误差随置0比例系数的变化曲线如图 3所示, 可以看出, 在置0比例系数相同的情况下, 2D-LDCT字典重构误差较小, 即能用更少的系数稀疏表示反射波分量, 2D-UDWT字典的稀疏表示能力适中, 而离散曲波变换字典重构误差最大, 即不能很好的稀疏表示反射波分量; 面波信号在各变换字典下的重构误差随置0比例系数的变化曲线如图 4所示, 可以看出, 在置0比例系数相同的情况下, 离散曲波变换字典重构误差较小, 即能用更少的系数稀疏表示面波分量, 而2D-LDCT字典的重构误差最大, 即不能很好的稀疏表示面波分量, 2D-UDWT字典的重构误差介于两者之间。离散曲波变换字典和2D-LDCT字典体现出的对2种信号成分稀疏表示能力的差异, 使得基于这2种字典组合的MCA方法可用于地震记录中面波信号和反射波信号的分离。

图 3(Fig. 3)

图 3 反射波信号在3种不同字典下稀疏性比较

图 4(Fig. 4)

图 4 含强能量面波地震信号在3种不同字典下稀疏性比较

在实际处理中, 求解(8) 式中的最优化问题得到的系数维度较高, 需占用大量内存, 因此, 将(8) 式进一步转化:

(13)

式中:r_r, r_g表示字典T_r, T_g零空间中的任一向量; T_r⁺=T_r^T(T_rT_r^T)^－1表示T_r的伪逆; T_g⁺表示T_g的伪逆。(13) 式将求解表示系数α_r, α_g转化为求解两类信号分量s_r和s_g, 同时引入残余向量r_r和r_g, 根据边界函数、BCR算法^[25]和最大后验概率(Maximum A Posteriori, MAP)估计, 可以设置r_r=0, r_g=0^[12], 则最优化问题变为:

(14)

上述问题的求解以BCR算法为基础, 分别固定一种信号成分不变, 更新另一种信号分量, 通过迭代阈值法最小化该分量的惩罚项和约束项得到更新系数, 再与其相应的表示字典相乘即可得到更新后的信号成分。在对面波成分进行更新时做相应的修改, 利用离散曲波变换字典对其进行5层分解后, 保留曲波域中主要含面波的尺度和方向上的系数不变(主要分布在Fine尺度和Detail尺度的低频高波数方向), 对其它系数进行阈值处理, 使得反射波分量与面波分量的分离效果更加明显。

图 5a为利用雷克子波合成的有效反射波记录, 该记录共240道, 每道采样点数为1000, 道间距为30m, 包含3个反射同相轴, 主频为40Hz, 最大振幅为3, 信号采样间隔为5ms, 记录长度为5s;图 5b为加入的面波, 面波利用扫描信号合成, 频带范围为5~15Hz, 最大振幅为12;图 5c为合成的受面波干扰的地震记录, 其信噪比为-22.46dB, 在面波与反射波重叠的区域, 反射波几乎被全部掩盖, 受损严重。应用本文方法对图 5c含面波合成地震记录进行处理, 设置迭代次数为50, 分离出的反射波分量和面波分量如图 6a和图 6b所示。为了进一步说明本文方法的有效性, 分别采用基于2D-UDWT与2D-LDCT字典组合的MCA法和曲波变换法对该地震记录进行面波压制处理。其中, 2D-UDWT字典和2D-LDCT字典分别稀疏表示面波分量和反射波分量, 曲波变换使用快速离散wrapping算法。图 6c和图 6d为基于2D-UDWT与2D-LDCT字典组合的MCA法得到的反射波分量和面波分量, 图 6e和图 6f分别为曲波变换法处理后得到的反射波分量和面波分量。为了更清晰地显示3种方法的压制效果, 分别选取合成反射波记录和各方法处理后得到的反射波记录的第100道数据进行波形和振幅谱分析, 得到的结果如图 7所示。

图 5(Fig. 5)

图 5 合成地震记录 a合成反射波记录; b合成面波记录; c受面波干扰的合成地震记录

图 6(Fig. 6)

图 6 合成地震记录面波压制效果对比 a本文方法得到的反射波; b本文方法得到的面波; c基于2D-UDWT与2D-LDCT字典组合的MCA法得到的反射波; d基于2D-UDWT与2D-LDCT字典组合的MCA法得到的面波; e曲波变换法得到的反射波; f曲波变换法得到的面波

图 7(Fig. 7)

图 7 合成地震记录单道波形及振幅谱 a 图 5a第100道记录的波形; b 图 7a的振幅谱; c 图 6a第100道记录的波形; d 图 7c的振幅谱; e 图 6c第100道记录的波形; f 图 7e的振幅谱; g 图 6e第100道记录的波形; h 图 7g的振幅谱

由图 6a和图 6b可以看出, 合成地震数据中的面波和反射波得到了有效的分离, 反射波分量中只有极少部分的面波残留; 分离出的面波分量中也几乎不含反射波, 图 7c中分离出的反射波第100道记录的波形和合成的反射波单道记录的波形非常接近。由图 6c和图 6d的分离结果可以看出, 该方法能较好地去除面波分量, 但同时也损失了少量反射波信息, 图 7e中反射波单道记录的波形中还存在部分干扰。由图 6e和图 6f中的处理结果可以看出, 该方法能较好地去除面波分量, 同时能很好地保留反射波分量的信息, 面波分量中几乎不含反射波, 但图 7g中反射波单道记录的波形中还含有较多干扰。由图 7中各种方法分离出的反射波单道记录的振幅谱可以看出, 相比于本文方法, 基于2D-UDWT与2D-LDCT字典组合的MCA法和曲波变换法都损失了部分低频成分的有效信号。以上各种方法的面波压制效果可以用处理后的信噪比和均方根误差(RMSE)来客观评价, 得到的结果如表 1所示。从表 1中可以看出, 采用本文基于离散曲波变换和2D-LDCT字典组合的MCA方法压制面波后, 地震数据的信噪比提高最多, 均方根误差最小, 这表明本文方法能够对地震记录中的面波进行更有效地去除。

表 1(Tab. 1)

表 1 不同方法压制面波效果对比

表 1 不同方法压制面波效果对比

为验证本文方法对实际地震记录处理的有效性, 将此方法用于压制图 8所示的地震记录中的面波, 该地震记录共240道, 每道采样点数为1024, 采样间隔为4ms。图 9a和图 9b展示了其第115道记录的波形和振幅谱, 分析可知, 面波主要分布在低频部分0~20Hz范围内。设置迭代次数为30次, LDCT窗口宽度为16, 分离后的反射波和面波如图 10a和图 10b所示。为了进一步说明本文方法的有效性, 分别应用基于2D-UDWT与2D-LDCT字典组合的MCA法和曲波变换法对该地震记录进行处理。其中, 2D-UDWT字典和2D-LDCT字典分别稀疏表示面波分量和反射波分量, 曲波变换使用快速离散wrapping算法。图 10c和图 10d为基于2D-UDWT与2D-LDCT字典组合的MCA法处理后得到的反射波分量和面波分量, 图 10e和图 10f为曲波变换法处理后得到的反射波分量和面波分量。为了更清晰地显示3种方法的压制效果, 分别选取利用各方法处理后得到的反射波记录的第115道记录进行波形和振幅谱分析, 如图 11所示。

图 8(Fig. 8)

图 8 含面波干扰的实际地震记录

图 9(Fig. 9)

图 9 实际地震记录第115道记录的波形(a)及其振幅谱(b)

图 10(Fig. 10)

图 10 实际地震记录面波压制效果对比 a本文方法得到的反射波; b本文方法得到的面波; c基于2D-UDWT与2D-LDCT字典组合的MCA法得到的反射波; d基于2D-UDWT与2D-LDCT字典组合的MCA法得到的面波; e曲波变换法得到的反射波; f曲波变换法得到的面波

图 11(Fig. 11)

图 11 实际单道地震记录的波形及振幅谱对比 a 图 10a第115道记录的波形; b 图 11a的振幅谱; c 图 10c第115道记录的波形; d 图 11c的振幅谱; e 图 10e第115道记录的波形; f 图 11e的振幅谱

从图 10中可以看出, 应用本文方法分离出来的反射波分量同相轴较为清晰, 并且较好地保留了原地震记录的有效信息, 分离出的面波只含有微弱的反射波; 利用基于2D-UDWT与2D-LDCT字典组合的MCA法和曲波变换法得到的反射波中都几乎不含面波信号, 但分离出来的面波中都含有少量明显的反射波信号。从图 11中单道记录的波形及振幅谱可以看出, 相比于本文方法得到的结果, 基于2D-UDWT与2D-LDCT字典组合的MCA法和曲波变换法得到的反射波单道记录的波形都存在不同程度的幅度失真, 并且两种方法得到的反射波都损失了部分低频成分的信息。这表明, 本文方法压制面波的效果较好, 能够很好地将反射波和面波分量分离, 同时对有效信号的损害较小。

本文利用MCA方法压制地震数据中的面波, 根据面波和有效反射波的形态结构差异, 分别选取表示性能更好的离散曲波变换和优化后的2D-LDCT字典作为两种信号分量的稀疏表示字典, 建立地震记录在两种字典下的面波分离模型, 在采用块协调松弛算法来求解该模型时, 将曲波域内主要含面波的尺度和方向上的系数保持不变, 对其他系数进行阈值处理, 获得了明显的压制效果。和基于2D-UDWT与2D-LDCT字典组合的MCA法、曲波变换法相比, 本文方法在有效压制面波干扰的同时能很好地保护反射波信号的波形特征, 是一种更精确有效的面波压制方法。然而实际地震数据中往往含有多种干扰, 使用变换字典组合的方法进行面波压制, 会存在将原有其它干扰转变为新的干扰噪声的问题, 如何有效规避此问题还需要进一步研究。