随着全球陆上油气勘探不断向复杂构造区域、深层和超深层以及非常规储层延伸，地震信号采集效率与数据质量之间的矛盾日益突出。可控震源(vibroseis)技术凭借其卓越的可控性、良好的环境适应性和强大的重复激发能力，成为替代炸药震源的工业化标准方案，并得到广泛认可和应用[1-3]。近年来，可控震源高效采集技术取得显著进展，形成交替扫描、滑动扫描、多点同步扫描以及高保真采集等一系列先进作业模式，使日采集效率从原来的几百炮跃升至上万炮，有力地推动陆上可控震源勘探进入高效率和高密度采集的新阶段[4-6]。然而，高效率的采集模式不可避免地引入多种背景噪声，对数据质量造成严重的干扰。这些背景噪声包括如下几个方面: 1) 由震源液压伺服系统的非线性振动引起的谐波干扰，其频率为基波扫描信号的整数倍，并在时频域内与有效信号产生混叠[7-8]; 2) 滑动扫描模式下，相邻炮信号在时间和空间上的混叠，经过相关处理后形成伪同相轴，导致构造解释出现多解性[9-10]; 3) 可控震源机械系统引起的地面振动，产生沿近地表传播的低频机械噪声，这种噪声频带较宽，初始能量强，对资料信噪比构成显著的负面影响[11]。上述背景噪声成分均呈现显著的非平稳特性(如时变振幅、频散相位以及非线性调频特征)，对基于平稳性假设的传统噪声压制方法构成根本性的挑战[12-13]。因此，如何有效地消除可控震源在高效采集过程中伴生的背景噪声，是获取高质量地震数据的关键环节。
针对可控震源背景噪声压制问题，国内外学者提出多种方法体系，主要分为室外采集技术优化和室内数据处理两大方向。在室外采集技术优化方面，研究者通过改进震源系统和优化扫描信号设计，从源头减少谐波干扰。可控震源技术的进步显著地推动了非线性扫描信号的应用，特别是在提升地震数据质量方面取得突破性进展。Bagaini[14]和 Ziolkowski[15]通过改进扫描信号设计，使可控震源能够有效地激发低频信号并降低畸变。蓝加达[16]通过调整非线性扫描速率，延长高频段扫描时间来补偿衰减，显著地提升哈萨克斯坦 S区的勘探分辨率。陶知非[17]通过改进液压伺服阀响应速度和振动器机械结构，提高可控震源高频输出上限，在一定程度上改善了可控震源高频输出信号的品质。魏福吉等[18]基于目标层频谱特征设计扫描信号，并利用频谱曲线反演时频函数，虽然提高了华北地区 1500 m深层资料的信噪比，但该方法对有效频宽构成一定的限制，不利于后期处理。徐雷良等[19]提出阻尼雷克子波模型，在保持窄旁瓣特性的同时，扩展了低频成分。
背景噪声抑制是处理数据的核心内容，现有技术主要分为两类。一类为传统滤波方法。Li等[20]通过深入地分析扫频信号，构建纯相移滤波算子，实现对单道数据中谐波噪声的有效压制。黄建平等[21]在此基础上进行拓展，通过统计地面力信号，设计适用于多炮数据的相移滤波器，但该方法仅适用于相关前数据中的高次谐波压制。Sicking等[22]提出基于地面力信号分解的频域滤波方法，通过精确地计算基波与谐波的频谱比值来设计滤波器，但处理过程中需从最后一炮数据开始，无法实现单炮操作。另一类为稀疏优化与自适应改进方法。毛海波等[23]提出基于自适应字典学习的优化方法，通过数据驱动，动态地调整字典，有效地克服固定字典在处理复杂地震资料时的局限性，显著地提升数据保真性和噪声压制效果。
针对上述传统方法存在的局限性，本文旨在解决可控震源高效采集过程中非平稳背景噪声(谐波干扰和机械振动等)难以被有效地压制且有效信号易失真的核心问题，技术路线如下: 首先，构建基于扫描信号与未相关可控震源数据中有效信号时频先验关系的非平稳信号分解模型，利用扫描信号与有效信号的时频关联性建立约束条件; 然后，采用正交匹配追踪算法，自适应地分离未相关可控震源数据中的有效信号与复杂非平稳背景噪声，从而突破传统方法压制效果不佳的问题; 最后，通过合成数据测试验证方法在复杂噪声干扰下的有效性，结合塔里木盆地戈壁低信噪比实际数据，与 F-K域预测滤波法进行对比，验证本文方法的优越性。
1 方法原理
1.1 频率域自回归模型
假设已获取一个大小为 N×M的未相关的可控震源地震记录 d(x，t)，其中 x ∈{1，2，…，N}和 t ∈{1，2，…，M}分别为空间和时间采样点。首先，对时间域数据 d(x，t)进行傅里叶变换，得到频率域地震数据 D(x，f ) ( f ∈(f1，f2，…，fI))，I为频率点个数。然后，对频率域地震数据进行单一频率切片处理，对于给定频率值fk，k ∈( 1，2，…，K)，该频率值对应的频率切片 D(x，fi)表示为
对单个频率切片 D(x，fi)，分别实施前向和后向 AR自回归模型预测。
前向预测: 利用前 L个数据点预测当前点右侧信号，计算公式为
其中，L是自回归滤波器阶数(空间预测的步数);al(fk)为前向滤波器系数，通过最小二乘法求解:
式中，E为单位矩阵，μ 为预白化系数，
为共轭转置，[·] T表示转置运算。Gforward代表前向设计矩阵，为 Toeplitz结构，如下式所示:
后向预测: 利用后 L个数据点预测当前点左侧信号，计算公式为
其中，L是自回归滤波器长度(空间预测的步数);bl(fk)为前后滤波器系数，求解方式与向后 AR自回归模型预测类似:
式中，
为共轭转置。Gbackward代表后向设计矩阵，为 Hankel结构，如下式所示:
最后，将向前和向后预测结果进行叠加，并对中间区域取平均值，以便减少边界误差，同时进行傅里叶逆变换，获得时空域数据 Dout( x， t):
1.2 非平稳信号分解技术压制背景噪声
假设已有一个大小为 N×M，未与扫描信号相关前的可控震源地震记录 d(x，t)，其中 x ∈{1，2，…，N}和 t ∈{ 1，2，…，M}分别为空间方向与时间方向采样点，非平稳信号分解技术对时空域地震数据中每一列信号进行处理。对其中第 i列数据进行如下处理:
可控震源的采集通常将预设好的正弦线性升频扫描信号 s(t)作为震源激发信号，其表达式为
其中，扫描信号 s(t)的振幅 A(t)和频率 f(t)均为时间函数，f(t)的表达式为
式中，fs代表扫描信号的初始频率，fd代表扫描信号的截至频率，T代表扫描信号的时间长度。
在理想情况下，检波器接收到的地震信号是多个地下反射界面与扫描信号褶积后的总和，而单个地下反射界面与扫描信号褶积后的信号与扫描信号之间只有振幅以及起始时间不同。因此，未与扫描信号相关前的可控震源地震单列记录可以表示为
式中，K为地下介质反射系数的数量; ak(t)是第 k个地下反射界面与扫描信号褶积后信号的未知变化的振幅; f (t)是第 k个地下反射界面与扫描信号褶积后信号的已知的变化频率; φk是第 k个地下反射界面与扫描信号褶积后信号的初始相位，也就是与扫描信号的差异时间延迟; n(t)为背景噪声。
因此，未与扫描信号相关前的可控震源地震单列记录可以表示为如下矩阵形式:
其中，di为地震数据中第 i列信号，矩阵 S中的元素Sk = sin[f ( t ) +φk]，矩阵 a =[a1 ( t )，a2 (t )，…，a K(t)]。
由于矩阵 Sk已知，问题转化为从 di中求解振幅系数 ak(t)。因此，可以构造如下目标函数:
式中，λ 为大于 0的正则化参数，用于调解模型的泛化能力与噪声抑制效果。可以利用正交匹配追踪算法解决此问题。当残差的二范数小于人为设置的噪声阈值ε 时，表明信号分解较为完全，可表示为
本文方法处理可控震源地震数据的流程如图1所示。
图1 非平稳信号分解技术流程
Fig.1 Flowchart of technology based on nonstationary signal decomposition
1.3 信噪比计算
为了更好地验证本文方法的有效性，将信噪比(SNR)作为定量评价指标，计算公式为
其中，Ps代表信号的平均功率，Pn代表噪声的平均功率。
2 模型测试
非平稳信号分解技术专门用于处理未相关的可控震源数据，该技术能够从夹杂复杂背景噪声的未相关可控震源数据中分离并提取出有效信号，有助于提高经与扫描信号互相关处理后的地震数据质量。本研究通过设计合成数据来验证非平稳信号分解技术提取复杂背景噪声下未相关可控震源数据中有效信号的性能。
首先，设计一个包含 5个不同深度波阻抗界面的合成数据模型，界面相对阻抗值分别为 0.5，0.4，0.1，-0.2和-0.4，对应的地层反射系数序列如图2所示。震源子波采用线性升频扫描信号，其参数设置如下: 时长为 20 s，采样间隔为 1 ms，频率范围为1.5～84 Hz，其波形如图3所示。为避免吉布斯效应，在扫描信号两端施加余弦斜坡函数进行加窗处理(图3中红色虚线椭圆所示)。
图2 合成波阻抗界面示意图
Fig.2 Schematic diagram of synthetic acoustic impedance interfaces
图3 扫描信号波形(上图)与局部放大(下图)
Fig.3 Sweep signal waveform diagram (up) and part zoom (down)
在理想条件(即忽略地层介质对能量的吸收等因素)下，通过将扫描信号与合成的波阻抗界面进行褶积运算，得到未相关的可控震源数据。图4(a)为在未考虑谐波等背景噪声干扰情况下，扫描信号与合成波阻抗界面褶积后的地震数据波形图(即经过未相关处理的可控震源数据)，由于扫描信号的长时特性，难以清晰地区分 5个波阻抗界面。然而，图4(b)所示未相关可控震源数据的时频域图中清晰地揭示，5个有效信号的频率斜率与扫描信号的频率斜率保持一致，相互独立，但能量存在明显的差异。进一步对比可知，时频域 5个倾斜信号的起始点与合成波阻抗界面 5个系数的时间位置完全吻合，同时能量的差异也与这些系数的绝对值相对应。
图4 扫描信号谐波干扰前波形频率
Fig.4 Frequency of waveform before harmonic interference of sweep signal
为真实地模拟实际可控震源数据采集过程中复杂背景噪声的干扰情况，首先在扫描信号的基础上叠加一、二、三阶谐波分量(图5(a)～(c))，这些谐波的频率分别是扫描信号频率的 1，2和 3倍(图5(d)～(f))，其能量强度依次递减为扫描信号的 80%，50%和 40%。从图5可以清晰地观察到，随着谐波频率斜率逐渐增加，其能量强度逐渐减弱。
图5 一、二、三阶谐波的波形和时频
Fig.5 First，second and third harmonic waveforms and time-frequency
在扫描信号中添加谐波干扰后，时域波形出现明显的畸变(图6(b)中红色虚线椭圆所示)。为了清晰地展示谐波的影响，取 0～2 s的波形进行展示(图6(a)和(b))，在时频域可观察到高频干扰成分(图6(c)和(d))。
图6 扫描信号谐波干扰后波形频率
Fig.6 Frequency of waveform after harmonic interference of sweep signal
在原有谐波干扰的基础上，进一步给扫描信号添加随机噪声干扰。在不考虑地层介质吸收等复杂情况的前提下，将带背景噪声的扫描信号与合成波阻抗界面进行褶积运算，从而获得带背景噪声的未相关可控震源数据。图7(a)展示带背景噪声的扫描信号与合成波阻抗界面褶积后的总响应(即带背景噪声未相关可控震源数据)。从图4(a)波形图中难以清晰地分辨 5个波阻抗界面，同时相较于图7(a)，可以明显地观察到背景噪声的干扰。这表明背景噪声使得信号复杂度增加，有效信息被掩盖。图7(b)展示带背景噪声未相关可控震源数据的时频域，除有效信号外，还可见两类噪声，一类为斜率更大的谐波噪声(图7(b)中红色箭头所示)，另一类为覆盖全局的随机噪声(图7(b)中红色虚线方框所示)。这些噪声严重地影响有效信号的识别。
图7 带噪扫描信号与地层系数褶积
Fig.7 Noisy sweep signal and wave impedance interface after convolution
本文对合成数据测试中正则化参数 λ 与噪声阈值 ε 的选取遵循“在保护有效信号的情况下尽量多压制噪声”的原则，结合合成数据中谐波噪声(能量为基波的 80%，50%和 40%)与全局随机噪声的整体强度，经过多次测试，确定 λ 的取值为 0.001。该取值可有效地避免过拟合，同时充分压制背景噪声。考虑到合成数据的噪声功率大小，确定噪声阈值 ε的取值为 1.2×10-3。对带背景噪声的未相关可控震源数据进行非平稳信号分解处理，结果如图8所示。观察波形图和时频图可知，经非平稳信号分解处理后，可控震源中有效信号得以完整地保留下来，处理结果清晰地显示出与图4所示结果的一致性，证明了非平稳信号分解技术在受到复杂背景噪声干扰的未相关可控震源数据中提取有效信号卓越能力，也就是说，该技术在压制未相关可控震源背景噪声方面表现出显著的效果。
图8 带噪未相关可控震源数据非平稳信号分解结果
Fig.8 Results of nonstationary signal decomposition of uncorrelated noisy vibroseis data
同时，通过式(18)计算处理前后的信噪比，进一步证明非平稳信号分解技术的有效性。经过计算，原始数据的信噪比为-6.25 dB，经过非平稳信号分解后的信噪比为 7.67 dB。
3 实际资料处理
为了更深入地验证非平稳信号分解技术在复杂背景噪声下的未相关可控震源数据中提取有效信号及压制背景噪声的优越性，相较于 F-K域预测反褶积技术，选取塔里木盆地塔西南戈壁工区的可控震源高效采集实际数据作为测试样本。该工区地表主要由砂砾石构成，中浅层地质构造复杂，导致原始单炮数据的信噪比低，背景噪声种类繁多且非常活跃。实际数据采集过程中采用扫描采集技术。
将本文选取的未相关可控震源数据作为测试数据，每炮包含 230道，每道记录时长为 30 s，采样间隔为 1 ms。从图9可以清晰地看出，由于激发震源是长扫描信号，具有长时性，无法从中直接提取类似炸药震源产生的地震数据信息，因此需要与扫描信号进行互相关，之后才能得到类似炸药震源的具有可解释性的地震数据。
图9 未相关可控震源数据
Fig.9 Uncorrelated vibroseis data
对选取的未相关可控震源数据，分别采用 F-K域预测反褶积技术以及非平稳信号分解技术进行处理。
F-K域预测反褶积技术的参数设置如下: 预测步长为 2 (戈壁工区地震道间距小(高密度采集)，道间相关性强，小步长能在压制噪声的同时，避免过度平滑导致同相轴失真)，频率截断为 85 Hz (完全匹配扫描信号的最高频率 84 Hz，精准地剔除 85 Hz以上的高频谐波噪声(可控震源的谐波是基波整数倍，85 Hz以上无有效信号))，预白化因子为 0.001 (极小值保证反褶积过程不引入额外畸变，适配低信噪比数据的稳定处理)，预测算子长度为 0.5 (折中“噪声压制效果”和“信号保真性”，避免过长导致过拟合，过短导致噪声压制不彻底)。
对于非平稳信号分解技术，设定频率斜率为扫描信号频率变化斜率即可，而非平稳信号分解技术中选择正则化参数 λ 与噪声阈值 ε 时遵循“在保护有效信号的情况下尽量地多压制噪声”原则。结合工区背景噪声(高频谐波和机械振动噪声)的分布特征以及原始数据信噪比水平，经过测试，当 λ 取值为0.002时，能够平衡好实际数据背景噪声压制程度与避免过度压制有效信号。对于噪声阈值，选取 ε =2×10-3，能够保证不进行无意义的过多迭代，达到终止效果。
从图10(a)可以看出，仅仅经过互相关处理后的野外数据受到各种复杂因素干扰，包含大量背景噪声，使得野外数据的可解释性较差。如图10(b)所示，F-K域预测反褶积技术的背景噪声抑制效果不尽如人意，在近中偏移距浅深层地震数据中仍然存在背景噪声干扰，同相轴不够连续。如图10(c)所示，经过非平稳信号分解处理后，全局地震数据同相轴清晰且连续，证明对背景噪声的压制较为彻底。
图10 可控震源数据处理结果对比
Fig.10 Comparison of the vibroseis data processing results
(a)为互相关处理后结果，(b)为经过 F-K 域预测反褶积技术处理后结果，(c)为经过非平稳信号分解处理后结果
为进一步展示本文方法处理效果的细节，将结果中 3～5 s时段进行局部放大。放大视图(图11(a))中红色方框区域显示，图11(a)中背景噪声显著，尤以道号 140附近为甚，同相轴几乎无法识别; 图11(b)和(c)中可以看出非平稳信号分解技术得到的地震数据同相轴更加连贯和清晰，而 F-K域预测反褶积技术得到的地震数据中仍然残留部分背景噪声的干扰。总体而言，非平稳信号分解技术处理后的地震数据在信噪比方面的表现更优异。
图11 可控震源数据局部放大(3～5 s)结果对比
Fig.11 Comparison of local amplification (3-5 s) results of the vibroseis data
(a)为互相关处理后结果，(b)为经过 F-K 域预测反褶积技术处理后结果，(c)为经过非平稳信号分解处理后结果，(d)为 F-K 域预测反褶积技术处理后结果与互相关处理后结果的残差，(e)为非平稳信号分解技术处理后结果与互相关处理后结果的残差
4 结论
本文针对可控震源高效采集中的背景噪声问题，提出一种基于非平稳信号分解的噪声压制方法。通过构建扫描信号与未相关可控震源数据中有效信号的时频先验约束模型，并借助正交匹配追踪算法，实现对谐波干扰和机械振动等非平稳背景噪声与有效信号的自适应分离，突破传统滤波技术在噪声平稳性假设方面的限制，为复杂噪声的压制提供了新的思路。合成数据验证表明，在包含多阶谐波干扰(能量达到基波的 80%)以及全局随机噪声干扰的复杂背景噪声条件下，该方法仍然能够准确地提取未相关可控震源数据中的有效信号。塔里木盆地实际数据的应用结果表明，与 F-K域预测滤波法相比，经过非平稳信号分解技术处理后的地震数据同相轴更加连续和清晰，证明其能够有效地压制背景噪声并保留有效信号。
[1]佟训乾，林君，姜弢，等.陆地可控震源发展综述.地球物理学进展，2012，27(5): 1912-1921
[2]倪宇东，王井富，马涛，等.可控震源地震采集技术的进展.石油地球物理勘探，2011，46(3): 349-356
[3]周大同，周恒，张慕刚，等.可控震源施工效率估算方法.石油地球物理勘探，2008，43(增刊 2): 50-54
[4]周如义，魏铁，张新峰，等.可控震源高效交替扫描作业技术及应用.石油地球物理勘探，2008，43(增刊 2): 15-20
[5]宁宏晓，唐东磊，皮红梅，等.国内陆上“两宽一高”地震勘探技术及发展.石油物探，2019，58(5): 645-653
[6]汪恩华，赵邦六，王喜双，等.中国石油可控震源高效地震采集技术应用与展望.中国石油勘探，2013，18(5): 24-34
[7]曹务祥，张慕刚.滑动扫描谐波分析.石油地球物理勘探，2005，40(5): 499-503
[8]李振春，曲英铭，韩文功，等.可控震源两种谐波产生机理与特征研究.石油物探，2016，55(2): 159-172
[9]曲英铭.可控震源特征噪音压制及混叠数据成像研究[D].青岛: 中国石油大学(华东)，2018
[10]张慕刚，祝杨，董烈乾，等.可控震源超高效混叠采集技术及应用.地球物理学进展，2021，36(3):1176-1182
[11]王彦峰，周旭，段孟川，等.一种可控震源机械振动噪声的压制方法、装置和设备: CN202310683 328.7[P].2024-12-10: 3-4
[12]王立歆，陈国金，崔树果，等.基于地震干涉法的可控震源黑三角噪声压制技术及应用.石油物探，2021，60(6): 925-936
[13]魏周宏，马红伟，潘远洲.可控震源低频谐波畸变压制技术.物探装备，2015，25(5): 293-296
[14]Bagaini C.Low-frequency vibroseis data with maximum displacement sweeps.The Leading Edge，2008，27(5): 582-591
[15]Ziolkowski A.Review of vibroseis data acquisition and processing for better amplitudes: adjusting the sweep and deconvolving for the time-derivative of the true groundforce.Geophysical Prospecting，2010，58(1):41-53
[16]蓝加达.可控震源非线性扫描在高分辨率地震采集中的应用.石油物探，2008，47(2): 208-211
[17]陶知非.改善可控震源高频信号输出品质的探讨.物探装备，2008，18(2): 71-77
[18]魏福吉，赵国勇，张剑中，等.基于目的层频谱的可控震源非线性扫描信号的设计方法.CN20151009 6211.4[P].2015-09-30: 7-9
[19]徐雷良，张剑，赵国勇.基于阻尼雷克子波的可控震源非线性扫描信号设计方法.石油地球物理勘探，2022，57(3): 540-549
[20]Li X P，Soellner W，Hubral P.Elimination of harmonic distortion in vibroseis data.Geophysics，1995，60(2):503-516
[21]黄建平，周学锋，郭军，等.滑动扫描记录中压制谐波干扰方法.中国石油大学学报(自然科学版)，2012，36(2): 81-85
[22]Sicking C，Fleure T，Nelan S，et al.Slip sweep harmonic noise rejection on correlated shot data // SEG Technical Program Expanded Abstracts 2009.Tulsa: Society of Exploration Geophysicists，2009: 36-40
[23]毛海波，马俊彦，王晓涛，等.基于自适应字典学习的可控震源数据谐波噪声压制方法.石油物探，2020，59(5): 725-735