可控震源因其具有安全、环保、低耗、高效、方式灵活等优点而备受青睐, 在全球陆上石油地震勘探采集中占有很大的比重。可控震源采集与处理技术得到了广泛的研究与应用^[1-5], 但与之相关的谐波干扰制约了可控震源采集技术的应用。谐波干扰产生的原因主要有两个:一是由于可控震源机械装置和震动装置以及液压伺服系统的非线性导致的谐波干扰; 二是震板与大地的耦合效应产生的谐波干扰。谐波干扰损害了可控震源采集资料的品质, 降低了资料的分辨率^[6-8]。根据谐波的特征差异, 可分为力信号内含谐波和表层响应谐波。李振春等^[9]通过分析可控震源正演模拟数据和实际数据, 揭示了力信号内含谐波和表层响应畸变谐波的产生机理和特征。

国内外学者针对谐波干扰的压制进行了大量的研究。采用设计合理的采集方式可以在一定程度上压制谐波干扰, 例如WOMACK等^[10-11]论述了利用相位编码消除可控震源同步激发产生的谐波干扰。SILVERMAN等^[12]提出变相位分离单炮数据并进行谐波压制的方法。WARD等^[13]分析了SILVERMAN变相位方法的缺点, 认为这种方法理论上不能消除偶次谐波干扰, 并提出通过增加不同相位的扫描次数来消除偶次谐波干扰。曹务祥^[14]提出了利用组合扫描压制谐波畸变的技术, 通过分段设计各子扫描频宽构建扫描组达到压制谐波干扰的目的。蓝益军等^[15]提出了现场压制滑动扫描谐波干扰的方法, 利用谐波干扰出现的时间, 合理设计采集参数, 避免滑动扫描谐波干扰的产生。涂伟等^[16]详细分析了采集中分段扫描的谐波干扰压制效果。不同于常规噪声, 谐波干扰噪声压制难度更大, 因此发展了一系列针对谐波干扰的室内压制方法。主要包括相移法、频率域滤波法和预测相减法。相移法是针对基波信号与谐波信号的相位差异进行压制的。LI等^[17]基于线性频率扫描方式, 提出了纯相移滤波(PPSM)方法来压制谐波干扰, 该方法简单高效、直观稳定, 特别适用于VSP资料, 但不适用于滑动扫描记录。因此, 黄建平等^[18]对其进行改进, 提出了一种可以适用于滑动扫描技术的相移滤波法。频率域相减法则针对有效信号与谐波噪声在频率域中的时间差异进行谐波压制。RAS等^[19]对滑动扫描纪录中的谐波干扰进行正演模拟, 并在频率域通过选择合适的滑动扫描时间分离谐波干扰。随后, FLEURE^[20]改善了频率域分离谐波并进行压制的方法。预测滤波法采用的方式是对谐波噪声进行预测并进行相减, 得到压制谐波噪声的炮记录, 该方法理论基础强, 压制谐波噪声准确, 受到了越来越多学者的关注。MEUNIER等^[21]采用的方法是从最后一炮开始向前进行处理, 通过互相关算法预测各次谐波干扰, 然后减去谐波干扰。MOERIG^[22]对道记录进行相关, 利用负时间域的信号来作为对前一炮谐波干扰的估计, 然后预测并消除谐波。MARTIN等^[23]提出了一种谐波预测分离技术。田新琦等^[24]利用基波与谐波相关前后的特征设计谐波算子, 并将其作用于相关前原始地震记录以压制谐波干扰。曹务祥等^[25]提出利用可控震源地面力信号压制谐波干扰的方法。曲英铭等针对力信号内含谐波干扰采用自适应的谐波干扰压制方法, 在地面力记录不准确的情况下对谐波干扰进行准确压制, 但该方法仅适用于力信号中包含的谐波压制, 当地表耦合性较差时, 表层响应畸变产生的谐波能量很强, 无法采用该方法进行压制^[26]。倪宇东^[27]提出了地面力信号滤波法, 并详细论述了该方法的基本原理和实际应用效果。韩文功等^[28]提出了一种基于反射系数和地面力信号的谐波干扰消除法。除此之外, 不同的滤波方法也被应用于压制谐波干扰。李凤磊等^[29]提出了基于褶积滤波的滑动扫描谐波干扰压制方法。李合群^[30]提出一种滤波法谐波干扰压制方法。尚永生^[31]提出一种基于稀疏反演的谐波干扰压制方法。林娟等^[32]应用独立同步扫描邻炮干扰压制方法滤除滑动扫描资料上的谐波干扰。曲英铭等^[33]对远距离同步滑动扫描(DS4)采集方式可控震源产生的谐波干扰和同步交涉干扰同时进行了压制。

前人研究的谐波干扰压制方法主要针对的是力信号中包含的谐波, 而没有考虑表层响应产生的谐波, 实际上, 这两种类型的谐波产生机理与特征不同^[9], 谐波压制方法也存在较大差异。因此本文针对两种谐波干扰产生机理和特征采用不同的压制方法进行联合压制。最后, 对存在2种类型谐波干扰的实际资料进行了谐波压制实验。

1 两种谐波干扰联合压制方法

首先, 根据李振春等^[9]分析的两种谐波产生机理与特征, 谐波干扰可以分为两种:①可控震源非线性畸变使得震板输出的力信号与输入的扫描信号不同产生的谐波干扰, 这种谐波干扰无法避免, 但可以通过震板输出的力信号记录下来, 称之为力信号内含谐波; ②可控震源力信号向下传播时由于震板与地表存在不耦合及可能发生的谐振现象, 导致输出的力信号经过地表后发生畸变, 与真实向下传播的信号不同, 称之为表层响应谐波, 此类谐波不能被震板记录下来。这两种谐波不仅产生机理不同, 而且特征也有较大的不同。力信号内含谐波干扰的特征为:①频带范围为有效波频带范围的整数倍; ②当采用线性升频扫描信号时, 出现在负时间轴上, 主要对前一炮有效能量产生干扰; ③与有效波在相位上存在明显的差异; ④力信号内含谐波遍布整个剖面, 但近偏移距影响更严重。表层响应畸变产生的谐波干扰特征为:①频带范围比较窄; ②当采用线性升频扫描信号时, 出现在正时间轴上, 主要污染当前炮有效能量; ③能量非常强, 分布于近偏移距处; ④在炮集上表现出明显的随机性, 有的单炮受表层响应谐波污染严重, 有的单炮则几乎不存在表层响应谐波。

由两种谐波干扰产生的机理及特征可以看出, 力信号中包含的谐波与表层响应畸变产生的谐波在频率特征、出现位置、采集方式、出现时间、相位特征上都存在明显的差异, 因此其压制方法也必然有较大的不同。通过合理的野外采集参数优化, 可以减少两种谐波的产生, 本文研究的重点是在室内进行谐波压制。因为可控震源数据相关前所占内存巨大, 相关后所占内存大大减少, 所以通常可控震源采集的海量数据只存储相关后数据。本文讨论了对相关后可控震源数据中的力信号内含谐波和表层响应谐波的联合压制方法。

1.1 力信号内含谐波的压制方法

力信号内含谐波信号可以通过震板输出的力信号记录下来, 因此我们可以通过力信号预测出力信号内含谐波噪声。

基于SERIFF等^[34]对谐波的理论研究, 以线性升频扫描信号为例, 定义参考信号如下:

$ {{\mathit{s}}_{\mathit{k}}}\left( \mathit{t} \right)\text{=}{{\mathit{a}}_{\mathit{k}}}\left( \mathit{t} \right)\text{sin }\!\![\!\!\text{ 2}\mathit{\pi kf}\left( \mathit{t} \right)\mathit{t}\text{ }\!\!]\!\!\text{ } $

(1)

式中:t为扫描时间; k为正整数, 表示谐波阶次; a_k(t)为各次谐波的振幅值; f(t)为瞬时频率。a_k(t)由下式给出:

$ {{\mathit{a}}_{\mathit{k}}}\left( \mathit{t} \right)\text{=}\left\{ \begin{align} &{{\mathit{c}}_{\mathit{k}}}\left( \mathit{t} \right)\ \ \ \ \ \ \ \ \ \text{0}\le \mathit{t}\le \mathit{T} \\ &\ \ \ \ \text{0}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \text{其它} \\ \end{align} \right. $

(2)

对于线性降频扫描信号, 公式(1)中的f(t)可表示为:

$ \mathit{f}\left( \mathit{t} \right)\text{=}{{\mathit{f}}_{\text{H}}}\text{-}\frac{{{\mathit{f}}_{\text{H}}}\text{-}{{\mathit{f}}_{\text{L}}}}{\text{2}\mathit{T}}\mathit{t} $

(3)

式中:T为扫描长度; f_L为扫描起始频率; f_H为扫描终止频率。通过测试对比, 公式(2)中c_k(t)取如下窗函数:

$ {{\mathit{c}}_{\mathit{k}}}\left( \mathit{t} \right)\text{=}\left\{ \begin{matrix} \frac{{{\mathit{A}}_{\mathit{k}}}}{2}\left[\text{1+cos }\!\!\pi\!\!\text{ }\left( \frac{\mathit{t}}{{{\mathit{T}}_{\text{1}}}}\text{+1} \right) \right]&\text{0}\le \mathit{t<}{{\mathit{T}}_{\text{1}}} \\ {{\mathit{A}}_{\mathit{k}}}&{{\mathit{T}}_{\text{1}}}\le \mathit{t}\le \mathit{T}\text{-}{{\mathit{T}}_{\text{1}}} \\ \frac{{{\mathit{A}}_{\mathit{k}}}}{2}\left[\text{1-cos }\!\!\pi\!\!\text{ }\left( \frac{\mathit{T}\text{-}\mathit{t}}{{{\mathit{T}}_{\text{1}}}} \right) \right]&\mathit{T}\text{-}{{\mathit{T}}_{\text{1}}}\mathit{<t}\le {{\mathit{T}}_{\text{1}}} \\ \end{matrix} \right. $

(4)

式中:T₁表示陡坡长度; A_k表示振幅。

由力信号内含谐波信号的特征可以得出:力信号内含谐波信号的频率是基波信号频率的整数倍。谐波畸变信号是基波信号和各次谐波分量的线性叠加, 则谐波畸变信号可表达为:

$ {{\mathit{s}}_{\text{D}}}\left( \mathit{t} \right)\text{=}\sum\limits_{\mathit{k}\text{=1}}^{\mathit{M}}{{{\mathit{s}}_{\mathit{k}}}}\left( \mathit{t} \right)\text{=}\sum\limits_{\mathit{k}\text{=1}}^{\mathit{M}}{{{\mathit{a}}_{\mathit{k}}}}\left( \mathit{t} \right)\text{sin(2}\mathit{\pi kft}\text{)} $

(5)

式中:M表示谐波的最高阶次。

基于褶积模型可知, 第i炮某道检波器所接收到的连续相关前记录是该炮地面力信号与该炮炮点和检波点组成的共中心点处反射系数的褶积:

$ {{\mathit{x}}_{\mathit{i}}}\left( \mathit{t} \right)\text{=}{{\mathit{g}}_{\mathit{i}}}\left( \mathit{t} \right)\text{*}{{\mathit{r}}_{\text{im}}}\left( \mathit{t} \right) $

(6)

式中:x_i(t), g_i(t)和r_im(t)分别表示第i炮相关前记录、地面力信号和共中心点位置的反射系数; “*”表示褶积运算。

而地面力信号是基波信号和各阶谐波分量的线性叠加, 因此方程(6)中的g_i(t)可表示为:

$ {{\mathit{g}}_{\mathit{i}}}\left( \mathit{t} \right)\text{=}{{\mathit{b}}_{\mathit{i}}}\left( \mathit{t} \right)\text{+}\sum\limits_{\mathit{k}\text{=2}}^{\mathit{M}}{{{\mathit{h}}_{\mathit{ik}}}}\left( \mathit{t} \right) $

(7)

式中:b_i(t)表示第i炮地面力信号分离出来的基波信号; h_ik(t)表示第i炮地面力信号分离出来的k次谐波信号。分离方程(7)中基波信号b_i(t)和谐波信号h_ik(t)的方法有很多^[35-37], 我们根据基波信号和谐波信号的差异采用相移法进行分离^[26]。

这样最终相关后记录等于相关前记录与震源扫描信号的互相关:

$ {{{\mathit{{x}'}}}_{\mathit{i}}}\left( \mathit{t} \right)\text{=}{{\mathit{x}}_{\mathit{i}}}\left( \mathit{t} \right)\otimes \mathit{s}\left( \mathit{t} \right) $

(8)

式中:x′_i(t)为第i炮某道地震相关记录; ⊗表示互相关运算。

由公式(6)、公式(7)和公式(8)得到相关记录:

$ \begin{align} &{{{\mathit{{x}'}}}_{\mathit{i}}}\left( \mathit{t} \right)\text{=}{{\mathit{b}}_{\mathit{i}}}\left( \mathit{t} \right)\text{*}{{\mathit{r}}_{\text{im}}}\left( \mathit{t} \right)\otimes \mathit{s}\left( \mathit{t} \right)\text{+} \\ &\ \ \ \ \ \ \ \ \ \ \ \sum\limits_{\mathit{k}\text{=2}}^{\mathit{M}}{{{\mathit{h}}_{\mathit{ik}}}}\left( \mathit{t} \right)\text{*}{{\mathit{r}}_{\text{im}}}\left( \mathit{t} \right)\otimes \mathit{s}\left( \mathit{t} \right) \\ \end{align} $

(9)

式中, 等号右侧第一项表示有效波项, 第二项表示谐波项。将公式(9)中的谐波项变换到频率域, 得:

$ {{\mathit{X}}_{\mathit{ih}}}\left( \mathit{f} \right)\text{=}\sum\limits_{\mathit{k}\text{=2}}^{\mathit{M}}{{{\mathit{H}}_{\mathit{ik}}}}\left( \mathit{f} \right){{\mathit{R}}_{\text{im}}}\left( \mathit{f} \right)\mathit{S}\left( \text{-}\mathit{f} \right)\text{ }\!\!~\!\!\text{ } $

(10)

式中, X_ih(f)表示第i炮某道所产生的全部谐波干扰, 简记为:

$ {{\mathit{X}}_{\mathit{ih}}}\left( \mathit{f} \right)\text{=}\mathit{H}\left( \mathit{f} \right)\mathit{R}\left( \mathit{f} \right)\mathit{S}\left( \text{-}\mathit{f} \right) $

(11)

公式(11)就是基于褶积模型的频率域谐波干扰计算公式。可以看出:谐波干扰是谐波分量与反射系数的褶积再与扫描信号的互相关。

由于反射系数在实际资料处理过程中难以准确求取, 本文根据公式(6)、公式(8)和公式(11)可以得到如下谐波干扰计算公式, 从而避免求取反射系数。

$ {{\mathit{X}}_{\mathit{ih}}}\left( \mathit{f} \right)\text{=}\frac{{{{\mathit{{X}'}}}_{\mathit{i}}}\left( \mathit{f} \right)\mathit{H}\left( \mathit{f} \right)\mathit{S}\left( \text{-}\mathit{f} \right)}{\mathit{G}\left( \mathit{f} \right)\mathit{S}\left( \text{-}\mathit{f} \right)}\text{=}\mathit{C}\left( \mathit{f} \right){{{\mathit{{X}'}}}_{\mathit{i}}}\left( \mathit{f} \right)\text{ }\!\!~\!\!\text{ } $

(12)

由公式(12)可知:利用地面力信号G(f)、谐波分量H(f)、扫描信号S(f)设计滤波器C(f), 然后作用于相关记录X′_i(f), 可求出该炮产生的谐波干扰, 最后从含谐波噪声的炮记录中减去, 即可得到压制谐波噪声后的炮记录。

由上述推导过程可以看出, 该方法基于地下介质为完全弹性的假设, 但在实际情况下, 此条件很难得到满足。特别是在沙漠地区, 可控震源在地表以上激发, 表层的衰减非常严重, 而地层的粘滞性对有效波和谐波噪声的衰减不同, 因此用预测的结果不能完全压制力信号内含谐波。可以采用两种方法对本算法进行改进:①对炮记录进行反Q滤波, 以改善地层粘滞性对力信号内含谐波压制的影响; ②引入最小二乘反馈机制, 在最小二乘理论意义下通过计算基于L₂模的滤波因子最小化来预测谐波噪声和实际谐波噪声的差^[26], 压制流程如图 1所示。

1.2 表层响应谐波的压制方法

表层响应谐波干扰的能量通常比有效波能量强得多。从空间分布上来看并不是所有炮都存在表层响应谐波干扰, 这与地表情况和震源状态等有关, 总的来看表层响应谐波干扰在炮集上存在一定的随机性。表层响应谐波频率通常出现在某个狭窄的频率段上。能量强、窄频带和随机性构成了表层响应谐波的3个特征, 本文根据这3个特征研究了基于分频滤波的表层响应谐波压制技术。

将采集数据转换到共中心点(CMP)域之后, 谐波干扰仍然分布在近偏移距附近，且能量很强, 但谐波干扰并没有连成片, 中间夹杂着没有受谐波干扰的正常数据道, 如图 2所示。应用基于分频处理的谐波压制技术能够很好地对谐波进行压制, 并且不损失有效信号, 保幅性较好。

基于分频处理的谐波压制技术原理如图 3所示, 实现步骤如下。

1) 在数据上划分纵向时窗, 在时窗内进行傅里叶变换, 将数据变换到频率域, 然后进行分频带分析。图 3a所示为4 Hz一个频带范围进行的分析, 其中红色曲线表示异常振幅道。

2) 在每个频带内计算各道的平均振幅。即:

$ \mathit{E}\text{=}\frac{\text{1}}{\mathit{n}}\sum\limits_{\mathit{i}\text{=1}}^{\mathit{n}}{{{\mathit{A}}_{\mathit{i}}}} $

(13)

式中:E为各道的平均振幅; A_i为瞬时振幅; n为道数。

3) 在每个频带内按照振幅属性应用空间中值滤波器来识别强能量谐波干扰。

4) 将识别出的谐波干扰道清零, 然后利用周围道插值获得新的道, 将数据合并后转换到时间域, 完成一个时窗内的处理过程。

按照上述步骤完成所有时窗的处理后完成一个道集的处理。

需要指出的是, 在不同地表条件下, 表层响应谐波都满足窄频带和随机性特征, 但不一定满足强能量特征。因此实际应用该方法的条件为谐波干扰能量要明显强于有效波能量, 如果不满足, 直接通过叠加就可以实现对炮集上具有随机性的表层响应谐波的压制。同时, 在应用该方法前, 需要对炮记录做地表一致性校正。

1.3 谐波联合压制流程

本文提出的谐波联合压制方法的实现流程如图 4所示。在输入相关后的炮记录时, 首先要判断该炮记录受何种谐波干扰的污染。当只存在一种谐波干扰时, 采用对应的谐波干扰压制方法即可, 但当两种谐波干扰都存在进行联合压制时, 首先要进行力信号内含谐波干扰压制, 然后再进行表层响应畸变的谐波干扰压制。

此外, 进行力信号内含谐波压制时, 因为力信号内含谐波干扰产生的谐波是在负时间轴上, 对前一炮的有效能量造成干扰, 因此, 需要用后一炮的力信号预测前一炮产生的谐波干扰。不过在实际采集中, 两炮之间扫描信号与力信号内含谐波信号不会有太大的差异, 因此可以用当前炮的力信号设计谐波压制滤波器。

2 实际资料试算 2.1 力信号内含谐波压制

首先, 我们对采用滑动扫描采集方式得到的二维试验数据进行试算。采集道数为250道, 记录长度为8 s, 激发方式为一台可控震源一次激发, 滑动时间为8 s, 扫描频率为6~72 Hz, 扫描长度为16 s, 驱动幅度为70%, 斜坡长度为500 ms, 扫描方式为线性升频。

对于滑动扫描采集方式, 当采用线性升频扫描信号时, 力信号内含谐波干扰出现在负时间轴上, 该谐波包含两部分:一部分是当前炮深层反射产生的谐波干扰, 这部分谐波出现在浅层, 能量较弱; 另一部分是后一炮产生的谐波干扰, 对前一炮的有效信号造成干扰^[9], 这部分谐波能量强, 因此图 5只展示前一炮, 即受力信号内含谐波干扰影响严重的炮。从图 5a所示谐波压制前炮记录可以看出, 谐波主要出现在红色椭圆区域, 并不在本炮的近偏移距处, 原因是该谐波干扰不是当前炮产生的, 而是后一炮产生在负时间轴作用于当前炮的结果, 因此该谐波出现在后一炮的近偏移距处, 远离当前炮的偏移距。从图 5a还可以看出, 该炮记录只存在力信号内含谐波, 因此采用本文提出的基于力信号的预测滤波方法进行压制, 结果如图 5b所示。可以看出, 经过谐波压制后, 炮记录上谐波干扰得到了准确的压制, 而深部有效能量没有受到损害。图 5c为压制的谐波干扰, 可以看出, 本方法准确地压制了谐波干扰。图 5d为滤波算子c(t)在时间域的显示。对受谐波干扰影响的第一组炮记录进行水平叠加成像, 图 6为叠加剖面部分区域的放大图, 图 6a为谐波压制前的剖面, 可以看出因为谐波干扰的存在, 剖面杂乱无章, 几乎没有同相轴, 信噪比很低; 图 6b为谐波压制后的叠加剖面, 可以看出, 虽然二维试验数据覆盖次数较低, 但经过谐波压制后出现了明显的同相轴, 谐波干扰得到了明显的压制, 有效能量得到凸显。

2.2 表层响应谐波压制

对我国西部某探区实际资料采用分频滤波方法压制了表层响应谐波干扰。该资料采用交替扫描采集方式, 横向采样炮数为400, 记录长度为8 s, 激发方式为三台可控震源一次激发, 滑动时间为24 s, 扫描频率为5~84 Hz, 扫描长度为16 s, 驱动幅度为70%, 斜坡长度为500 ms, 扫描方式为线性升频。

图 7为采用分频滤波法压制表层响应谐波干扰前后的炮记录。从交替扫描的定义可知, 后一炮产生在负时间轴上的力信号内含谐波不会对当前炮造成影响, 而且当前炮产生的力信号内含谐波应出现在浅层区域, 能量较弱, 因此可以判断出图 7红框区域所示的谐波干扰不是力信号内含谐波, 而是表层响应谐波。对比两图中红框区域可以看出, 谐波压制前炮记录存在明显的谐波干扰, 谐波能量比有效波能量强得多。在深部区域, 谐波干扰与有效波混叠在一起, 覆盖了炮记录的有效信息。谐波压制后, 谐波干扰得到了很好的消除, 而且基本上没有造成有效能量的损失。图 8展示了谐波压制前后的叠加剖面, 从图 8a可以看出, 在红框所示的中深层区域存在明显的谐波干扰噪声, 造成剖面中部区域同相轴不连续, 降低了剖面的信噪比。而采用本文提出的分频滤波方法处理后(图 8b), 表层响应谐波干扰得到了明显的压制, 同相轴更加连续, 信噪比得到了很大的提升。通过该资料试算可以看出, 基于分频处理的谐波干扰压制技术能够很好地压制和去除表层响应谐波干扰, 单炮记录和叠加剖面的信噪比都得到了较大程度的提高。

2.3 两种谐波联合压制

图 9为同时存在力信号内含谐波和表层响应畸变谐波的相关后炮记录。该记录采用滑动扫描采集方式得到, 因此后一炮产生的力信号内含谐波干扰会出现在前一炮的有效能量区域, 对前一炮的有效能量造成干扰。与此同时, 该地区由于地表耦合条件差, 产生了明显的表层响应谐波。采集方式为滑动扫描, 滑动时间为8 s, 其它采集参数与2.1节中的相同。图 9中, 红框区域所示的谐波为力信号内含谐波(从左到右分别记作窗口1—窗口5), 蓝框区域所示的谐波为表层响应谐波(从左到右分别记作窗口1—窗口5), 为了进一步证明这一点, 图 10给出该区域的频谱图, 不同颜色的曲线表示不同窗口内的频谱。其中图 10a为力信号内含谐波, 从频谱图中可以看出, 该类型的谐波没有明显地集中到某一频率段内, 而图 10b所示的表层响应谐波能量集中于10~20 Hz区域内。

首先应用本文提出的力信号内含谐波压制方法压制力信号内含谐波, 结果如图 11所示。从图 11可以看出, 红框所示的力信号内含谐波得到了很好的压制, 而蓝框所示的表层响应谐波没有得到压制, 此结果也进一步证明了表层响应谐波不能通过地面力信号设计滤波器对其进行压制。力信号内含谐波压制之后, 利用本文提出的表层响应谐波压制方法对表层响应谐波进行压制, 结果如图 12所示, 可以看出, 经过进一步表层响应谐波压制后, 蓝框所示区域的谐波噪声也得到了很好压制。经过联合压制后, 图 9炮记录中的谐波噪声都得到了很好的压制, 也进一步证明了本文联合压制方法能够准确压制力信号内含谐波和表层响应谐波。

3 结论与认识

在可控震源野外采集的地震数据中存在两种谐波:力信号内含谐波和表层响应谐波, 这两种谐波产生机理及特征都存在很大的不同, 需要采用不同的谐波压制方法进行压制。本文针对这两类谐波提出了两种不同的谐波压制方法:基于地面力信号设计的谐波滤波器和基于共中心点域设计的分频滤波器。前者主要针对力信号内含谐波, 通过记录的力信号采用最小二乘的思想不断匹配预测真实的力信号内含谐波干扰, 达到准确消除力信号内含谐波干扰的目的, 此方法中可控震源力信号是野外可控震源采集及室内谐波压制的关键; 后者主要针对表层响应谐波, 这种谐波在某个频带范围内存在较强的谐波能量, 通过分频段在共中心点道集完成压制。两种谐波干扰压制方法组合, 能够准确压制谐波噪声。

需要注意的是, 两种谐波并不一定同时存在。力信号内含谐波产生在负时间轴, 该谐波包含两部分:一部分是当前炮深层反射产生的谐波干扰, 这部分谐波出现在浅层, 能量较弱; 另一部分是后一炮产生的谐波干扰, 对前一炮的有效信号造成干扰, 这部分谐波能量强, 且主要出现在滑动扫描及更高效的采集方式中。而表层响应谐波主要产生在正时间轴, 谐波干扰能量强弱主要取决于地表与震板耦合性以及谐波的阶次。在压制谐波之前需要根据两种谐波的特征对谐波类型进行判断并采用对应的谐波压制方法进行谐波压制, 一个最直接有效的方式就是对谐波的频谱进行分析。如果两种谐波都存在, 需要采用两种方法组合的方式, 但需要注意的是要先压制力信号内含谐波, 后压制表层响应谐波。

此外, 传统的谐波干扰压制方法是在相关前进行, 本文联合压制方法是针对相关后数据进行压制。可控震源相关前数据量巨大, 存储每一炮所需要的内存空间远远超过只存储相关后数据所需的内存空间, 对相关前数据进行谐波压制需要更大的计算量, 因此本文提出的相关后谐波干扰压制方法可以极大地提高计算效率、减少存储需求。

可控震源采集技术在未来地震勘探中具有广阔的应用前景, 其发展趋势向高效采集方向发展, 多炮之间的滑动时间越来越短, 而谐波干扰成为制约可控震源高效采集技术发展的瓶颈之一。因此, 谐波干扰的准确快速压制是推动可控震源采集技术发展的关键技术手段之一。