一种谐波噪声的压制处理方法及装置

摘要

本发明提供了一种谐波噪声的压制处理方法及装置，其中，该方法包括：首先，利用波形过零点次数的谐波噪声识别法识别待处理的地震道数据的多个样点中含谐波噪声的样点，将识别出的含谐波噪声的样点作为谐波样点，并对谐波样点进行标记；再利用分频的半波的平均振幅识别法确定各个频带中已标记的谐波样点是否需要压制，并对确定出的各个频带中需要压制的谐波样点进行谐波压制处理；最后，对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据。在本发明提供的实施例中，采用两种谐波噪声的识别方式对地震高效采集带来的谐波噪声进行有效识别，再对识别出的谐波噪声进行了有效压制，从而提高了地震资料的信噪比和品质。

权利要求

1.一种谐波噪声的压制处理方法，其特征在于，所述方法包括：

在待处理的地震道数据的多个样点中选取一个样点作为切除点，在选取的所述切除点处创建时窗长度为nwindows的时窗；根据所述谐波噪声的最大周期和所述时窗长度确定所述时窗内谐波噪声的过零点个数nmin；根据地震道数据中相邻的两个样点对应的振幅是否同号确定所述时窗内地震道数据的过零点个数ntrace；按照预设的时窗移动长度nmove将所述时窗向所述待处理的地震道数据的道结束点处方向移动；判断所述时窗内地震道的过零点个数ntrace是否大于所述谐波噪声的过零点个数nmin，若是，则将移出的样点确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识；直到所述时窗的下边缘点移动至所述地震道数据的道结束点；判断所述时窗内地震道的过零点个数ntrace是否大于所述谐波噪声的过零点个数nmin，若是，则将所述时窗内的所有样点均确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识；

对所述待处理的地震道数据进行频带划分处理，以划分得到多个不同频率段的频带；

从划分得到的多个所述频带中选取一个频带作为标准频带，并将所述标准频带之后的多个频带作为待压制频带；从多个所述待压制频带中，选取所述标准频带之后的两个频带作为当前待压制频带；在选取的两个所述当前待压制频带的已标记的谐波样点中，分别逐一判断两个所述当前待压制频带的各个所述谐波样点的所述当前待压制频带的半波的平均振幅是否大于所述标准频带的半波的平均振幅与插值函数的乘积；若是，则对所述谐波样点的地震道数据进行谐波噪声的压制处理；将压制后的两个所述当前待压制频带中与所述标准频带相邻的频带作为下一个标准频带；直到当前选取的下一个标准频带为最后一个频带；

对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据。

2.根据权利要求1所述的方法，其特征在于，所述在待处理的地震道数据的多个样点中选取一个样点作为切除点之前，还包括：设置谐波噪声的最小频率 并根据所述谐波噪声的最小频

率确定与所述最小频率对应的所述谐波噪声的最大周期 其中，fmin表示谐波噪声的最小频率，fbegin表示地震道正常子波振幅谱的开始频率，fend表示地震道正常子波振幅谱的结束频率，fharm表示开始识别的谐波噪声频率，nzero表示谐波噪声的最大周期。

3.根据权利要求1所述的方法，其特征在于，所述在选取的所述当前待压制频带的已标记的谐波样点中，分别逐一判断各个所述谐波样点的所述当前待压制频带的半波的平均振幅是否大于所述标准频带的半波的平均振幅与插值函数的乘积，包括：当所述谐波样点位于地震道数据的浅层，则判断各个所述谐波样点的所述当前待压制频带的半波的平均振幅是否大于所述标准频带的半波的平均振幅与浅层插值函数p1(t)的乘积；

当所述谐波样点位于地震道数据的中深层，则判断各个所述谐波样点的所述当前待压制频带的半波的平均振幅是否大于所述标准频带的半波的平均振幅与中深层插值函数p2(t)的乘积。

4.根据权利要求3所述的方法，其特征在于，所述浅层插值函数为

所述中深层插值函数为；

其中，dtime表示从切除时刻开始压制的最小门槛值，etime表示距道结束时刻压制的最小门槛，ldelaytm表示地震道数据的浅层与中深层间的临界点，lmute表示地震道数据中的选取的切除点，llength表示地震道数据中道结束点。

5.根据权利要求4所述的方法，其特征在于，对所述谐波样点的地震道数据进行谐波噪声的压制处理，包括：当所述谐波样点位于地震道数据的浅层时，则根据公式 对与所述标

准频带k-1相邻的所述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式 对所述标准频带k-1之后的第二个所述待压制频带k+1中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

当所述谐波样点位于地震道数据的中深层时，则根据公式 对与所

述标准频带k-1相邻的所述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式 对所述标准频带k-1之后的第二个所述待压制频带k+1中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

其中，t的取值为lmute至llength，yk-1(t)表示标准频带k-1中样点t的半波的平均振幅，yk(t)表示与标准频带k-1相邻的频带k中样点t的半波的平均振幅，xk(t)表示与标准频带k-1相邻的频带k中样点t的原始地震道数据，yk+1(t)表示标准频带k-1之后的第二个的频带k+1中样点t的半波的平均振幅，xk+1(t)表示标准频带k-1之后的第二个的频带k+1中样点t的原始地震道数据，γ表示振幅谱上振幅衰减系数。

6.根据权利要求1-5任一项所述的方法，其特征在于，对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据，包括：根据公式 确定压制后的地震道数据，其中，t表示样点号变量， 表示

频带k中样点t对应的压制后地震道数据，nfnum表示划分得到的频带总数。

7.一种谐波噪声的压制处理装置，其特征在于，所述装置包括：

第一识别模块，用于在待处理的地震道数据的多个样点中选取一个样点作为切除点，在选取的所述切除点处创建时窗长度为nwindows的时窗；根据所述谐波噪声的最大周期和所述时窗长度确定所述时窗内谐波噪声的过零点个数nmin；根据地震道数据中相邻的两个样点对应的振幅是否同号确定所述时窗内地震道数据的过零点个数ntrace；按照预设的时窗移动长度nmove将所述时窗向所述待处理的地震道数据的道结束点处方向移动；判断所述时窗内地震道的过零点个数ntrace是否大于所述谐波噪声的过零点个数nmin，若是，则将移出的样点确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识；直到所述时窗的下边缘点移动至所述地震道数据的道结束点；判断所述时窗内地震道的过零点个数ntrace是否大于所述谐波噪声的过零点个数nmin，若是，则将所述时窗内的所有样点均确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识；

第二识别模块，用于对所述待处理的地震道数据进行频带划分处理，以划分得到多个不同频率段的频带；从划分得到的多个所述频带中选取一个频带作为标准频带，并将所述标准频带之后的多个频带作为待压制频带；从多个所述待压制频带中，选取所述标准频带之后的两个频带作为当前待压制频带；在选取的两个所述当前待压制频带的已标记的谐波样点中，分别逐一判断两个所述当前待压制频带的各个所述谐波样点的所述当前待压制频带的半波的平均振幅是否大于所述标准频带的半波的平均振幅与插值函数的乘积；若是，则对所述谐波样点的地震道数据进行谐波噪声的压制处理；将压制后的两个所述当前待压制频带中与所述标准频带相邻的频带作为下一个标准频带；直到当前选取的下一个标准频带为最后一个频带；

数据重构模块，用于对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据。

8.根据权利要求7所述的装置，其特征在于，所述装置还包括：

参数设置模块，用于设置谐波噪声的最小频率 并根据所

述谐波噪声的最小频率确定与所述最小频率对应的所述谐波噪声的最大周期其中，fmin表示谐波噪声的最小频率，fbegin表示地震道正常子波振幅谱的开始频率，fend表示地震道正常子波振幅谱的结束频率，fharm表示开始识别的谐波噪声频率，nzero表示谐波噪声的最大周期。

9.根据权利要求7所述的装置，其特征在于，所述谐波样点识别与压制单元包括：

浅层谐波样点识别子单元，用于当所述谐波样点位于地震道数据的浅层，则判断各个所述谐波样点的所述当前待压制频带的半波的平均振幅是否大于所述标准频带的半波的平均振幅与浅层插值函数p1(t)的乘积；

中深层谐波样点识别子单元，用于当所述谐波样点位于地震道数据的中深层，则判断各个所述谐波样点的所述当前待压制频带的半波的平均振幅是否大于所述标准频带的半波的平均振幅与中深层插值函数p2(t)的乘积。

10.根据权利要求9所述的装置，其特征在于，所述浅层插值函数为

所述中深层插值函数为；

其中，dtime表示从切除时刻开始压制的最小门槛值，etime表示距道结束时刻压制的最小门槛，ldelaytm表示地震道数据的浅层与中深层间的临界点，lmute表示地震道数据中的选取的切除点，llength表示地震道数据中道结束点。

11.根据权利要求10所述的装置，其特征在于，所述谐波样点识别与压制单元包括：

浅层谐波样点压制子单元，用于当所述谐波样点位于地震道数据的浅层时，则根据公式 对与所述标准频带k-1相邻的所述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式 对所

述标准频带k-1之后的第二个所述待压制频带k+1中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

中深层谐波样点压制子单元，用于当所述谐波样点位于地震道数据的中深层时，则根据公式 对与所述标准频带k-1相邻的所述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式

对所述标准频带k-1之后的第二个所述待压制频带k+1中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

其中，t的取值为lmute至llength，yk-1(t)表示标准频带k-1中样点t的半波的平均振幅，yk(t)表示与标准频带k-1相邻的频带k中样点t的半波的平均振幅，xk(t)表示与标准频带k-1相邻的频带k中样点t的原始地震道数据，yk+1(t)表示标准频带k-1之后的第二个的频带k+1中样点t的半波的平均振幅，xk+1(t)表示标准频带k-1之后的第二个的频带k+1中样点t的原始地震道数据，γ表示振幅谱上振幅衰减系数。

12.根据权利要求7-11任一项所述的装置，其特征在于，所述数据重构模块包括：

数据重构单元，用于根据公式 确定压制后的地震道数据，其中，t表示样点号变量， 表示频带k中样点t对应的压制后地震道数据，nfnum表示划分得到的频带总数。

说明书

一种谐波噪声的压制处理方法及装置

技术领域

[0001] 本发明涉及地震勘探技术领域，具体而言，涉及一种谐波噪声的压制处理方法及装置。

背景技术

[0002] 现阶段的岩性地层和复杂构造勘探已经对地震勘探提出了薄储层小断层识别、精细构造圈闭目标刻画的高要求，因此，常规的地震勘探技术已经难以满足需求。“两宽一高”的地震勘探技术代表着地球物理技术的发展方向，它采用可控震源高效采集技术，具有较好的经济可行性，并且，通过大带道能力、功能强大的地震仪器及宽频可控震源来实现。

[0003] 目前，高密度三维地震采集技术普遍采用可控震源滑动扫描方式，虽然大幅度提高了采集效率，但是，在可控震源地震勘探的过程中，由于可控震源机械系统和液压伺服系统的非线性震动、以及震源与大地耦合等因素的影响，产生了谐波畸变。这种畸变以扫描信号频率范围的倍数出现，形成了谐波干扰。对于相关后的地震记录数据而言，谐波干扰将对本炮其自身产生高次谐波干扰(在初至前)或低次谐波干扰(在初至后的数据中)。尤其是当可控震源地震勘探采用滑动扫描方式时，后一炮不等前一炮扫描结束就已经开始激发，导致前后两炮甚至几炮数据部分叠加在一起，后一炮的高次谐波干扰将会产生其自身的谐波干扰(在初至前)和前一炮(甚至是前几炮)数据的邻炮干扰，前一炮的低次谐波干扰可能会产生其自身的谐波干扰(在初至后的数据中)和后一炮(甚至是后几炮)数据的邻炮干扰，从而形成了自身谐波干扰和邻炮谐波干扰。以下将自身谐波干扰和邻炮谐波干扰统称为谐波干扰，即谐波噪声。

[0004] 在可控震源地震勘探的过程中，这种谐波噪声严重降低了地震资料的信噪比，影响了地震资料的品质。

发明内容

[0005] 有鉴于此，本发明实施例的目的在于提供一种谐波噪声的压制处理方法及装置，通过有效识别地震道数据中的谐波噪声并对识别出的谐波噪声进行有效压制，可以提高地震资料的信噪比和品质。

[0006] 第一方面，本发明实施例提供了一种谐波噪声的压制处理方法，该方法包括：

[0007] 利用波形过零点次数的谐波噪声识别法识别待处理的地震道数据的多个样点中含谐波噪声的样点，将识别出的上述含谐波噪声的样点作为谐波样点，并对上述谐波样点进行标记；

[0008] 利用分频的半波的平均振幅识别法确定各个频带中已标记的谐波样点是否需要压制，并对确定出的各个频带中需要压制的谐波样点进行谐波压制处理；

[0009] 对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据。

[0010] 结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，上述利用波形过零点次数的谐波噪声识别法识别待处理的地震道数据的多个样点中含谐波噪声的样点之前，还包括：

[0011] 设置谐波噪声的最小频率 并根据上述谐波噪声的最小频率确定与上述最小频率对应的上述谐波噪声的最大周期 其中，fmin表示
谐波噪声的最小频率，fbegin表示地震道正常子波振幅谱的开始频率，fend表示地震道正常子波振幅谱的结束频率，fharm表示开始识别的谐波噪声频率，nzero表示谐波噪声的最大周期，单位为毫秒。

[0012] 结合第一方面的第一种可能的实施方式，本发明实施例提供了第一方面的第二种可能的实施方式，其中，上述利用波形过零点次数的谐波噪声识别法识别待处理的地震道数据的多个样点中含谐波噪声的样点，将识别出的上述含谐波噪声的样点作为谐波样点，并对上述谐波样点进行标记，包括：

[0013] 步骤a：在待处理的地震道数据中,以切除点作为开始样点，并在该切除点处创建时窗长度为nwindows的时窗；

[0014] 步骤b：根据上述谐波噪声的最大周期和上述时窗长度确定上述时窗内谐波噪声的过零点个数nmin；

[0015] 步骤c：根据地震道数据中相邻的两个样点对应的振幅是否同号确定上述时窗内地震道数据的过零点个数ntrace；

[0016] 步骤d：按照预设的时窗移动长度nmove将上述时窗向上述待处理的地震道数据的道结束点处方向移动；

[0017] 步骤e：判断上述时窗内地震道的过零点个数ntrace是否大于上述谐波噪声的过零点个数nmin，若是，则将移出的样点确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识；

[0018] 依次循环步骤b至步骤e，直到上述时窗的下边缘点移动至上述地震道数据的道结束点；判断上述时窗内地震道的过零点个数ntrace是否大于上述谐波噪声的过零点个数nmin，若是，则将上述时窗内的所有样点均确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识。

[0019] 结合第一方面，本发明实施例提供了第一方面的第三种可能的实施方式，其中，上述利用分频的半波的平均振幅识别法确定各个频带中已标记的谐波样点是否需要压制，并对确定出的各个频带中需要压制的谐波样点进行谐波压制处理，包括：

[0020] 步骤1：对上述待处理的地震道数据进行频带划分处理，以划分得到多个不同频率段的频带；

[0021] 步骤2：从划分得到的多个上述频带中选取一个频带作为标准频带，并将上述标准频带之后的多个频带作为待压制频带；

[0022] 步骤3：从多个上述待压制频带中，选取上述标准频带之后的两个频带作为当前待压制频带；

[0023] 步骤4：在选取的两个上述当前待压制频带的已标记的谐波样点中，分别逐一判断两个上述当前待压制频带的各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与插值函数的乘积；若是，则对上述谐波样点的地震道数据进行谐波噪声的压制处理；将压制后的两个上述当前待压制频带中与上述标准频带相邻的频带作为下一个标准频带；

[0024] 依次循环步骤3至步骤4，直到当前选取的下一个标准频带为最后一个频带。

[0025] 结合第一方面的第三种可能的实施方式，本发明实施例提供了第一方面的第四种可能的实施方式，其中，上述在选取的上述当前待压制频带的已标记的谐波样点中，分别逐一判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与插值函数的乘积，包括：

[0026] 当上述谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)，则判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与浅层插值函数p1(t)的乘积；

[0027] 当上述谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)，则判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与中深层插值函数p2(t)的乘积。

[0028] 结合第一方面的第四种可能的实施方式，本发明实施例提供了第一方面的第五种可能的实施方式，其中，上述浅层插值函数为

[0029] 上述中深层插值函数为；

[0030] 其中，dtime表示从切除时刻开始压制的最小门槛值，etime表示距道结束时刻压制的最小门槛，ldelaytm表示地震道数据的浅层与中深层间的临界点，lmute表示地震道数据中的切除点，llength表示地震道数据中道结束点。

[0031] 结合第一方面的第三种可能的实施方式，本发明实施例提供了第一方面的第六种可能的实施方式，其中，对上述谐波样点的地震道数据进行谐波噪声的压制处理，包括：

[0032] 当上述谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)时，则根据公式对与上述标准频带k-1相邻的上述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式 对上述
标准频带k-1之后的第二个上述待压制频带k+1中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

[0033] 当上述谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)时，则根据公式对与上述标准频带k-1相邻的上述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式 对上述
标准频带k-1之后的第二个上述待压制频带k+1中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

[0034] 其中，t的取值为lmute至llength，yk-1(t)表示标准频带k-1中样点t的半波的平均振幅，yk(t)表示与标准频带k-1相邻的频带k中样点t的半波的平均振幅，xk(t)表示与标准频带k-1相邻的频带k中样点t的原始地震道数据，yk+1(t)表示标准频带k-1之后的第二个频带k+1中样点t的半波的平均振幅，xk+1(t)表示标准频带k-1之后的第二个频带k+1中样点t的原始地震道数据，γ表示振幅谱上振幅衰减系数。

[0035] 结合第一方面至第一方面的第六种可能的实施方式中的任一种可能的实施方式，本发明实施例提供了第一方面的第七种可能的实施方式，其中，对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据，包括：

[0036] 根据公式 确定压制后的地震道数据，其中，t表示样点号变量，表示频带k中样点t对应的压制后地震道数据，nfnum表示划分得到的频带总数。

[0037] 第二方面，本发明实施例还提供了一种谐波噪声的压制处理装置，该装置包括：

[0038] 第一识别模块，用于利用波形过零点次数的谐波噪声识别法识别待处理的地震道数据的多个样点中含谐波噪声的样点，将识别出的上述含谐波噪声的样点作为谐波样点，并对上述谐波样点进行标记；

[0039] 第二识别模块，用于利用分频的半波的平均振幅识别法确定各个频带中已标记的谐波样点是否需要压制，并对确定出的各个频带中需要压制的谐波样点进行谐波压制处理；

[0040] 数据重构模块，用于对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据。

[0041] 结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，上述装置还包括：

[0042] 参数设置模块，用于设置谐波噪声的最小频率 并根据上述谐波噪声的最小频率确定与上述最小频率对应的上述谐波噪声的最大周期其中，fmin表示谐波噪声的最小频率，fbegin表示地震道正常子波振幅谱的开
始频率，fend表示地震道正常子波振幅谱的结束频率，fharm表示开始识别的谐波噪声频率，nzero表示谐波噪声的最大周期，单位为毫秒。

[0043] 结合第二方面的第一种可能的实施方式，本发明实施例提供了第二方面的第二种可能的实施方式，其中，上述第一识别模块包括：

[0044] 时窗创建单元，用于在待处理的地震道数据中,以切除点作为开始样点，并在该切除点处创建时窗长度为nwindows的时窗；

[0045] 第一确定单元，用于根据上述谐波噪声的最大周期和上述时窗长度确定上述时窗内谐波噪声的过零点个数nmin；

[0046] 第二确定单元，用于根据地震道数据中相邻的两个样点对应的振幅是否同号确定上述时窗内地震道数据的过零点个数ntrace；

[0047] 时窗移动单元，用于按照预设的时窗移动长度nmove将上述时窗向上述待处理的地震道数据的道结束点处方向移动；

[0048] 谐波噪声样点标记单元，用于判断上述时窗内地震道的过零点个数ntrace是否大于上述谐波噪声的过零点个数nmin，若是，则将移出的样点确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识；直到上述时窗的下边缘点移动至上述地震道数据的道结束点；判断上述时窗内地震道的过零点个数ntrace是否大于上述谐波噪声的过零点个数nmin，若是，则将上述时窗内的所有样点均确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识。

[0049] 结合第二方面，本发明实施例提供了第二方面的第三种可能的实施方式，其中，上述第二识别模块包括：

[0050] 频带划分单元，用于对上述待处理的地震道数据进行频带划分处理，以划分得到多个不同频率段的频带；

[0051] 标准频带选取单元，用于从划分得到的多个上述频带中选取一个频带作为标准频带，并将上述标准频带之后的多个频带作为待压制频带；

[0052] 当前待压制频带确定单元，用于从多个上述待压制频带中，选取上述标准频带之后的两个频带作为当前待压制频带；

[0053] 谐波样点识别与压制单元，用于在选取的两个上述当前待压制频带的已标记的谐波样点中，分别逐一判断两个上述当前待压制频带的各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与插值函数的乘积；若是，则对上述谐波样点的地震道数据进行谐波噪声的压制处理；将压制后的两个上述当前待压制频带中与上述标准频带相邻的频带作为下一个标准频带；直到当前选取的下一个标准频带为最后一个频带。

[0054] 结合第二方面的第三种可能的实施方式，本发明实施例提供了第二方面的第四种可能的实施方式，其中，上述谐波样点识别与压制单元包括：

[0055] 浅层谐波样点识别子单元，用于当上述谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)，则判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与浅层插值函数p1(t)的乘积；

[0056] 中深层谐波样点识别子单元，用于当上述谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)，则判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与中深层插值函数p2(t)的乘积。

[0057] 结合第二方面的第四种可能的实施方式，本发明实施例提供了第二方面的第五种可能的实施方式，其中，上述浅层插值函数为 上述中深层插值函数为；

[0058] 其中，dtime表示从切除时刻开始压制的最小门槛值，etime表示距道结束时刻压制的最小门槛，ldelaytm表示地震道数据的浅层与中深层间的临界点，lmute表示地震道数据中的切除点，llength表示地震道数据中道结束点。

[0059] 结合第二方面的第三种可能的实施方式，本发明实施例提供了第二方面的第六种可能的实施方式，其中，上述谐波样点识别与压制单元包括：

[0060] 浅层谐波样点压制子单元，用于当上述谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)时，则根据公式 对与上述标准频带k-1相邻的上述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式对上述标准频带k-1之后的第二个上述待压制频带k+1中需要压制的
谐波样点t对应的原始地震道数据进行谐波压制处理；

[0061] 中深层谐波样点压制子单元，用于当上述谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)时，则根据公式 对与上述标准频带k-1相邻的上述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式 对上述标准频带k-1之后的第二个上述待压制频带k+1中需要
压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

[0062] 其中，t的取值为lmute至llength，yk-1(t)表示标准频带k-1中样点t的半波的平均振幅，yk(t)表示与标准频带k-1相邻的频带k中样点t的半波的平均振幅，xk(t)表示与标准频带k-1相邻的频带k中样点t的原始地震道数据，yk+1(t)表示标准频带k-1之后的第二个频带k+1中样点t的半波的平均振幅，xk+1(t)表示标准频带k-1之后的第二个频带k+1中样点t的原始地震道数据，γ表示振幅谱上振幅衰减系数。

[0063] 结合第二方面至第二方面的第六种可能的实施方式中的任一种可能的实施方式，本发明实施例提供了第二方面的第七种可能的实施方式，其中，上述数据重构模块包括：

[0064] 数据重构单元，用于根据公式 确定压制后的地震道数据，其中，t表示样点号变量， 表示频带k中样点t对应的压制后地震道数据，nfnum表示划分得到的频带总数。

[0065] 在本发明提供的实施例中，首先，利用波形过零点次数的谐波噪声识别法识别待处理的地震道数据的多个样点中含谐波噪声的样点，将识别出的含谐波噪声的样点作为谐波样点，并对谐波样点进行标记；再利用分频的半波的平均振幅识别法确定各个频带中已标记的谐波样点是否需要压制，并对确定出的各个频带中需要压制的谐波样点进行谐波压制处理；最后，对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据。在本发明提供的实施例中采用两种谐波噪声的识别方式对地震高效采集带来的谐波噪声进行有效识别，再对识别出的谐波噪声进行了有效压制，从而提高了地震资料的信噪比和品质。

[0066] 为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

[0067] 为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

[0068] 图1示出了本发明实施例所提供的谐波噪声的压制处理方法的流程示意图；

[0069] 图2a示出了本发明实施例所提供的利用波形过零点次数的谐波噪声识别法识别过程中时窗移动的示意图；

[0070] 图2b示出了本发明实施例所提供的利用分频的半波的平均振幅识别法识别过程中分频后各个频带的地震道数据的示意图；

[0071] 图2c示出了本发明实施例所提供的利用分频的半波的平均振幅识别法识别过程中分频后各个频带的平均振幅的示意图；

[0072] 图2d示出了本发明实施例所提供的利用分频的半波的平均振幅识别法识别过程中各个频带中相邻频带重叠的示意图；

[0073] 图2e示出了本发明实施例所提供的利用分频的半波的平均振幅识别法识别过程中浅层与中深层的差值函数曲线的示意图；

[0074] 图3示出了本发明实施例所提供的通过实验得到的原始单炮数据、按照相关技术的方法进行谐波噪声压制后的单炮数据、以及按照本发明实施例的方法进行谐波噪声压制后的单炮数据的对比示意图；

[0075] 图4中示出了本发明实施例所提供的与图3中AC区域的三个单炮数据对应的频谱分析图；

[0076] 图5中示出了本发明实施例所提供的与图3中BCD区域的三个单炮数据对应的频谱分析图；

[0077] 图6中示出了本发明实施例所提供的谐波压制前的单炮数据的叠加剖面、按照本发明实施例的方法进行谐波噪声压制后的单炮数据的叠加剖面、以及叠加后的谐波噪声剖面示意图；

[0078] 图7示出了本发明实施例所提供的谐波噪声的压制处理装置的结构示意图。

具体实施方式

[0079] 下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0080] 考虑到在可控震源地震勘探的过程中，由于可控震源机械系统和液压伺服系统的非线性震动、以及震源与大地耦合等因素的影响，会产生谐波畸变，从而形成谐波噪声，这种谐波噪声会严重降低地震资料的信噪比，影响地震资料的品质。基于此，本发明实施例提供了一种谐波噪声的压制处理方法，以及一种可以应用该方法的谐波噪声的压制处理装置。

[0081] 谐波噪声除了具有频率范围成倍增长的特性之外，还具有相关特性、时间特性及振幅特性。需要对谐波噪声的特性进行分析，搞清谐波噪声的分布特性、频率特性及振幅特性，才能找出识别谐波噪声并进行压制的方法，从而提高地震资料信噪比。

[0082] 由于大地的吸收衰减作用，地震记录的主要能量通常集中在12～28Hz范围内，而采用滑动扫描方式得到的地震道数据的谐波噪声的频率通常大于30Hz；此外，滑动扫描每组可控震源的震动要间隔一个滑动时间，而震动间隔时间愈长，谐波噪声和有效信号的能量差异就愈大，这就为利用谐波频率和谐波的半波的平均振幅识别谐波噪声并对谐波噪声进行压制提供了条件。基于此，如图1所示，本发明实施例提供的一种谐波噪声的压制处理方法，该方法包括步骤S102-S106，具体如下：

[0083] 首先，设置谐波噪声的最小频率 并根据上述谐波噪声的最小频率确定与上述最小频率对应的上述谐波噪声的最大周期 其中，
fmin表示谐波噪声的最小频率，fbegin表示地震道正常子波振幅谱的开始频率，fend表示地震道正常子波振幅谱的结束频率，fharm表示开始识别的谐波噪声频率，nzero表示谐波噪声的最大周期，nzero的单位是毫秒；

[0084] 然后，对待处理的地震道数据中的多个样点进行谐波噪声的识别与压制处理，具体包括：

[0085] 步骤S102：利用波形过零点次数的谐波噪声识别法识别待处理的地震道数据的多个样点中含谐波噪声的样点，将识别出的上述含谐波噪声的样点作为谐波样点，并对上述谐波样点进行标记，其中，待处理的地震道数据为WA(t)，样点t对应的谐波噪声指示标志bflag(t)，将各个样点的初始bflag(t)设置为0；

[0086] 步骤S104：利用分频的半波的平均振幅识别法确定各个频带中已标记的谐波样点是否需要压制，并对确定出的各个频带中需要压制的谐波样点进行谐波压制处理；

[0087] 步骤S106：对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据。

[0088] 在本发明提供的实施例中采用两种谐波噪声的识别方式对地震高效采集带来的谐波噪声进行有效识别，再对识别出的谐波噪声进行了有效压制，从而提高了地震资料的信噪比和品质。

[0089] 具体的，上述利用波形过零点次数的谐波噪声识别法识别待处理的地震道数据的多个样点中含谐波噪声的样点，将识别出的上述含谐波噪声的样点作为谐波样点，并对上述谐波样点进行标记，包括：

[0090] 步骤a：在待处理的地震道数据中,以切除点作为开始样点，并在该切除点处创建时窗长度为nwindows的时窗；

[0091] 步骤b：根据上述谐波噪声的最大周期和上述时窗长度确定上述时窗内谐波噪声的过零点个数nmin，其中， nmin是一个固定值；

[0092] 步骤c：根据地震道数据中是否同号确定上述时窗内地震道数据的过零点个数ntrace，具体的，将地震道数据中相邻的两个样点所对应的振幅依次两两相乘，即WA(t)*WA(t+1)，若相乘结果小于零(相邻的两个样点对应的振幅异号，即一正值一负值)，表示该相邻两个样点之间过零点，则ntrace计数一次(数值加一)；若相乘结果大于零(相邻的两个样点对应的振幅同号，即同为正值或者同为负值)，表示该相邻两个样点之间未过零点，则ntrace数值不变；直到t＝该时窗中的最后一个样点，此时ntrace的数值标识该时窗内地震道数据的过零点个数；

[0093] 步骤d：按照预设的时窗移动长度nmove将上述时窗向上述待处理的地震道数据的道结束点处方向移动；

[0094] 步骤e：判断上述时窗内地震道的过零点个数ntrace是否大于上述谐波噪声的过零点个数nmin，若是(ntrace大于nmin)，则将移出的样点确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识(将移除的nmove个样点t至t+nmove对应的谐波噪声指示标志bflag(t)至bflag(t+nmove)设置为99)；若ntrace小于nmin，则移除的nmove个样点t至t+nmove对应的谐波噪声指示标志bflag(t)至bflag(t+nmove)的数值不变，仍为0；

[0095] 依次循环步骤c至步骤e，直到上述时窗的下边缘点移动至上述地震道数据的道结束点；判断上述时窗内地震道的过零点个数ntrace是否大于上述谐波噪声的过零点个数nmin，若是，则将上述时窗内的所有样点(样点llength-nwindows至llength)均确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识，其中，由于考虑到深层谐波干扰严重，因而，对于最后一个时窗进行单独处理，如最后一个时窗内ntrace大于nmin，则将该时窗内的所有样点均标记为含有谐波噪声的样点，将样点llength-nwindows至llength对应的bflag(llength-nwindows)至bflag(llength)设置为99。

[0096] 需要说明的是，考虑到对于每个时窗而言，谐波噪声的最大周期nzero、时窗长度nwindows均相同，因而，每个时窗内谐波噪声的过零点个数nmin是一个固定值，即每个时窗中nmin大小是不变的，因而，依次循环步骤c至步骤e即可。

[0097] 例如，如图2a所示，待处理的地震道数据各个样点t的取值为lmute至llength，每个时窗长度为nwindows，时窗移动长度为nmove，上述切除点为lmute(假如未预先设置切除点时，则认为从t＝1的样点作为开始识别的样点，并从t＝1的样点处开始创建时窗长度为nwindows的时窗)，举例说明，若待处理的地震道数据的时间序列是7000毫秒，采样间隔是2毫秒，因而地震道数据中包含有7000/2＝3500个样点，即地震道数据中道结束点llength＝3500；将样点1至样点lmute-1切除，其主要目的为切除初至波等干扰波，从而提高信噪比，可以将第500个样点作为切除点，即地震道数据中的切除点lmute＝500，将作为lmute＝500创建时窗的起始点，基于单位时间内过零点次数的谐波频率识别法从第500至3500的样点中确定出哪些样点含有谐波噪声，将含有谐波噪声的样点对应的谐波噪声指示标志bflag(t)的值设置为99。

[0098] 在本发明提供的实施例中，若确定出任一样点t含有谐波噪声，可以将该样点t对应的谐波噪声指示标志bflag(t)的值设置为某一预设值，例如，设置为99，用以指示该样点t含有谐波噪声，bflag(t)未设置预设值之前可以是0。通过不断的移动时窗，比较时窗内地震道的过零点个数ntrace与谐波噪声的过零点个数nmin的大小关系，进而确定待处理的地震道数据中各个样点中是否含有谐波噪声，完成基于波形过零点次数的谐波噪声识别法的谐波噪声识别，即识别完成后，含由谐波噪声的样点t对应的bflag(t)＝99。

[0099] 具体的，上述利用分频的半波的平均振幅识别法确定各个频带中已标记的谐波样点是否需要压制，并对确定出的各个频带中需要压制的谐波样点进行谐波压制处理，包括：

[0100] 步骤1：对上述待处理的地震道数据进行频带划分处理，以划分得到多个不同频率段的频带；

[0101] 其中，根据公式 确定划分得到的频带总数，其中，nfnum表示划分的频带个数，nreserve表示保留频带个数(nreserve＝2为常数参量)，fcutoff表示谐波噪声的截至频率，fharm表示开始识别的谐波噪声的频率，fbegin表示地震道数据正常子波振幅谱的开始频率；计算得到nfnum之后，根据公式 对原始
地震道数据进行分频处理，划分得到nfnum个带宽相等的频带，其中，hk(t)表示滤波因子，x(t)表示原始地震道数据，xk(t)表示滤波后数据，即分频后的数据， 为卷积云算符，t表示样点号变量，k表示频带变量；

[0102] 其中，仍以若待处理的地震道数据的时间序列是7000毫秒，采样间隔是2毫秒为例，因而地震道数据中包含有7000/2＝3500个样点，则分成nfnum个频带之后，每一个频带也包含有3500个样点。

[0103] 步骤2：从划分得到的多个上述频带中选取一个频带作为标准频带，并将上述标准频带之后的多个频带作为待压制频带；

[0104] 步骤3：从多个上述待压制频带中，选取上述标准频带之后的两个频带作为当前待压制频带；

[0105] 步骤4：在选取的两个上述当前待压制频带的已标记的谐波样点中，分别逐一判断两个上述当前待压制频带的各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与插值函数的乘积；若是，则对上述谐波样点的地震道数据进行谐波噪声的压制处理；将压制后的两个上述当前待压制频带中与上述标准频带相邻的频带作为下一个标准频带；其中，正常子波振幅谱是将标准能量频带内以半波为窗口的平均绝对值能量作为该半波中每一点的标准能量值，根据公式确定频带k中样点t的半波的平均振幅，其中，m1和m2分别表示频带k的起始样点号和结束样点号，xk(j)表示频带k中的样点j对应的原始数据；

[0106] 依次循环步骤3至步骤4，直到当前选取的下一个标准频带为最后一个频带。

[0107] 例如，如图2b所示，将待处理的地震道数据划分为6个带宽相等的频带，即nfnum＝6，k的取值为1至6，k＝1的低频频带不进行基于分频的谐波平均振幅识别法的谐波噪声识别与压制，如图2c所示，将k＝2的频带作为首次选取的标准频带(正常子波振幅谱)，对k＝3至k＝6的频带进行基于分频的谐波平均振幅识别法的谐波噪声识别与压制；正常子波振幅谱是以半波为窗口的平均振幅作为该半波中每一点的标准值，为通过比较频带间的半波的平均振幅来进行谐波噪声的识别与压制提供了基础，避免了相关技术中基于单点能量识别谐波干扰的缺陷。

[0108] 需要说明的是，待压制频带可以是划分得到的所有频带，也可以是划分得到的频率较高的频带。优选的，考虑到低频不被破坏，因此，频率最低的M-1个频带不作为待压制的频带，将k＝M-1(如k＝2)的频带作为标准频带，从k＝M(即k＝3)的频带开始进行压制，即将k＝M(k＝3)的频带至k＝nfnum(k＝6)的频带作为待压制的频带，其中，M为预设正整数。例如，当M＝3时，第1个频带不进行处理，第2个频带作为正常子波振幅谱标准频带，而第3个频带作为开始压制的谐波噪声频带。其中，划分得到的nfnum个频带，按照频率由低到高的顺序进行编号，编号越小的频带的频率越低，即第1个频带的频率最低，第nfnum个频带的频率最高。

[0109] 具体的，如图2d所示仍以将待处理的地震道数据划分为6个频带为例，将k＝2的频带作为标准频带，将k＝3至k＝6的频带作为待压制频带，即k＝2的频带的后1频带(k＝3)、后2频带(k＝4)、后3频带(k＝5)、后4频带(k＝6)；

[0110] 1)第一次识别与压制：将k＝2的频带之后的两个频带作为当前待压制频带即后1频带(k＝3)、后2频带(k＝4)，在后1频带(k＝3)、后2频带(k＝4)的已标记的谐波样点中，分别逐一对k＝3的频带和k＝4的频带中的已标记的谐波样点进行基于分频的谐波平均振幅识别法的识别与压制，即将各个谐波样点的半波的平均振幅与标准频带的半波的平均振幅与插值函数的乘积进行比较，根据比较结果判断谐波样点知否需要压制；然后，将压制后的k＝3的频带作为下一个标准频带；

[0111] 2)第二次识别与压制：将k＝3的频带之后的两个频带作为当前待压制频带即后2频带(k＝4)、后3频带(k＝5)，在后2频带(k＝4)、后3频带(k＝5)的已标记的谐波样点中，分别逐一对k＝4的频带和k＝5的频带中的已标记的谐波样点进行基于分频的谐波平均振幅识别法的识别与压制，即将各个谐波样点的半波的平均振幅与标准频带的半波的平均振幅与插值函数的乘积进行比较，根据比较结果判断谐波样点知否需要压制；然后，将压制后的k＝4的频带作为下一个标准频带(k＝4的频带经过两次识别与压制)；

[0112] 3)第三次识别与压制：将k＝4的频带之后的两个频带作为当前待压制频带即后3频带(k＝5)、后4频带(k＝6)，在后3频带(k＝5)、后4频带(k＝6)的已标记的谐波样点中，分别逐一对k＝5的频带和k＝6的频带中的已标记的谐波样点进行基于分频的谐波平均振幅识别法的识别与压制，即将各个谐波样点的半波的平均振幅与标准频带的半波的平均振幅与插值函数的乘积进行比较，根据比较结果判断谐波样点知否需要压制；然后，将压制后的k＝5的频带作为下一个标准频带(k＝5的频带经过两次识别与压制)；

[0113] 4)第四次识别与压制：将k＝5的频带之后的两个频带作为当前待压制频带即后4频带(k＝6)，由于k＝5的频带之后只有一个待处理的频带，在后4频带(k＝6)的已标记的谐波样点中，分别逐一对k＝6的频带中的已标记的谐波样点进行基于分频的谐波平均振幅识别法的识别与压制，即将各个谐波样点的半波的平均振幅与标准频带的半波的平均振幅与插值函数的乘积进行比较，根据比较结果判断谐波样点知否需要压制；然后，将压制后的k＝6的频带作为下一个标准频带(k＝6的频带经过两次识别与压制)，所有频带均识别与压制结束。

[0114] 综上可知，当将k-1频带作为标准频带，k频带(k-1频带的后1频带)和k+1频带(k-1频带的后2频带)分别为当前待压制频带：

[0115] 1)第一次识别与压制：k-1频带作为标准频带，对k频带和k+1频带进行基于分频的谐波平均振幅识别法的识别与压制；

[0116] 2)第二次识别与压制：k频带作为标准频带，对k+1频带和k+2频带进行基于分频的谐波平均振幅识别法的识别与压制；

[0117] 3)第三次识别与压制：k+1频带作为标准频带，对k+2频带和k+3频带进行基于分频的谐波平均振幅识别法的识别与压制；

[0118] 4)第四次识别与压制：若k+3为最后一个频带，k+2频带作为标准频带，只对k+3频带进行基于分频的谐波平均振幅识别法的识别与压制；其中，k+1频带、k+2频带、k+3频带均进行了两次进行基于分频的谐波平均振幅识别法的识别与压制。

[0119] 在本发明实施例中，采用两种谐波噪声的识别方式对地震高效采集带来的谐波噪声进行有效识别，即先采用基于单位时间内过零点次数的谐波频率识别法，再采用基于分频的谐波平均振幅识别法，其中，在采用基于分频的谐波平均振幅识别法进行识别时，实际上对于选取的标准频带之后的第二个频带起到最后一个频带(以划分6个频带，且选第2频带为标准频带为例，第4频带至第6频带)而言，基于分频的谐波平均振幅识别法进行了两次识别与压制，即高频频带需要进行两次识别和压制，可以对两个频带中重叠部分中的样点进行有效的压制，从而取得更好的谐波噪声压制效果。

[0120] 其中，考虑到谐波样点位于地震道数据的浅层与中深层(包括中层和深层两部分)时的谐波噪声干扰能力存在一定的差异，因而，如图2e所示，在比较半波的平均振幅时谐波样点位于地震道数据的浅层或者中深层引入不同的差值函数，提高了对谐波噪声进行识别与压制的准确度，基于此，上述在选取的上述当前待压制频带的已标记的谐波样点中，分别逐一判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与插值函数的乘积，包括：

[0121] 当上述谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)，则判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与浅层插值函数p1(t)的乘积；

[0122] 当上述谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)，则判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与中深层插值函数p2(t)的乘积。

[0123] 在本发明提供的实施例中，考虑到浅层和中深层谐波噪声的能量有差异，浅层弱，中深层强，因此，为了有效压制谐波噪声，而不伤及浅层和中深层的有效信号，本发明实施例的方法中采用不同的差值函数来识别浅层、中深层，相较于相关技术的方法采用浅层和中深层的识别参数不变的方法，本发明实施例的方法取得的效果更优。

[0124] 具体的，上述浅层插值函数为

[0125] 上述中深层插值函数为；

[0126] 其中，dtime表示从切除时刻开始压制的最小门槛值，etime表示距道结束时刻压制的最小门槛，ldelaytm表示地震道数据的浅层与中深层间的临界点，lmute表示地震道数据中的切除点，llength表示地震道数据中道结束点。

[0127] 其中，对上述谐波样点的地震道数据进行谐波噪声的压制处理，包括：

[0128] 当上述谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)时，则根据公式对与上述标准频带k-1相邻的上述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式 对上述
标准频带k-1之后的第二个上述待压制频带k+1中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

[0129] 当上述谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)时，则根据公式对与上述标准频带k-1相邻的上述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式 对上述
标准频带k-1之后的第二个上述待压制频带k+1中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

[0130] 其中，t的取值为lmute至llength，yk-1(t)表示标准频带k-1中样点t的半波的平均振幅，yk(t)表示与标准频带k-1相邻的频带k中样点t的半波的平均振幅，xk(t)表示与标准频带k-1相邻的频带k中样点t的原始地震道数据，yk+1(t)表示标准频带k-1之后的第二个的频带k+1中样点t的半波的平均振幅，xk+1(t)表示标准频带k-1之后的第二个的频带k+1中样点t的原始地震道数据，γ表示振幅谱上振幅衰减系数。

[0131] 具体的，将k-1频带作为首次选取的标准频带，当频带k中的谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)时，判断yk(t)是否大于yk-1(t)×p1(t)，若yk(t)＞yk-1(t)×p1(t)，则根据公式 对与标准频带k-1后1频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；以及，当频带k+1中的谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)时，判断yk+1(t)是否大于yk-1(t)×p1(t)，若yk+1(t)＞yk-1(t)×p1(t)，则根据公式 对与标准频带k-1后2频带k+1中需要压制的谐波样点
t对应的原始地震道数据进行谐波压制处理；

[0132] 当频带k中的谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)时，判断yk(t)是否大于yk-1(t)×p2(t)，若yk(t)＞yk-1(t)×p2(t)，则根据公式 对与标准频带k-1后1频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；以及，当频带k+1中的谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)时，判断yk+1(t)是否大于yk-1(t)×p2(t)，若yk+1(t)＞yk-1(t)×p2(t)，则根据公式
对与标准频带k-1后2频带k+1中需要压制的谐波样点t对应的原始地
震道数据进行谐波压制处理。

[0133] 具体的，对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据，包括：

[0134] 根据公式 确定压制后的地震道数据，其中，t表示样点号变量，表示频带k中样点t对应的压制后地震道数据，nfnum表示划分得到的频带总数。

[0135] 在本发明提供的实施例中采用两种谐波噪声的识别方式对地震高效采集带来的谐波噪声进行有效识别，再对识别出的谐波噪声进行了有效压制，从而提高了地震资料的信噪比和品质；进一步的，在采用基于分频的谐波平均振幅识别法进行谐波噪声识别时，对高频频带中的谐波样点进行了两次基于分频的谐波平均振幅识别法的识别和压制，因而可以对两个频带的重叠部分中的样点进行有效的压制，从而取得更好的谐波噪声压制效果；更进一步的，考虑到谐波样点位于地震道数据的浅层与中深层时的谐波噪声干扰能力存在一定的差异，因而，在比较半波的平均振幅时谐波样点位于地震道数据的浅层或者中深层引入不同的差值函数，提高了对谐波噪声进行识别与压制的准确度。

[0136] 通过实验得到的实验结果表明，相较于相关技术的方法，本发明实施例的方法能够降低地震资料的信噪比，取得更好的地震资料品质。图3中示出了原始单炮数据、按照相关技术的方法进行谐波噪声压制后的单炮数据、以及按照本发明实施例的方法进行谐波噪声压制后的单炮数据的对比示意图，图4中示出了图3中AC区域对应的原始单炮数据的频谱分析、按照相关技术的方法进行谐波噪声压制后的单炮数据的频谱分析、以及按照本发明实施例的方法进行谐波噪声压制后的单炮数据的频谱分析的对比示意图。图5中示出了图3中BCD区域对应的原始单炮数据的频谱分析、按照相关技术的方法进行谐波噪声压制后的单炮数据的频谱分析、以及按照本发明实施例的方法进行谐波噪声压制后的单炮数据的频谱分析的对比示意图。从图3、图4和图5中可以明显看出，与相关技术的方法相比，按照本发明实施例的方法进行谐波噪声压制后的数据更加干净。图4所示的频谱分析图表明，由于谐波噪声能量较强，30Hz-90Hz频率之间的谐波高频噪声比较突出；在低频1Hz-5Hz中，按照本发明实施例的方法进行谐波噪声压制后的单炮数据的频谱的有效成分要高于原始单炮数据的频谱以及按照相关技术的方法进行谐波噪声压制后的单炮数据的频谱；在5Hz-30Hz中，按照相关技术的方法进行谐波噪声压制后的单炮数据的频谱与按照本发明实施例的方法进行谐波噪声压制后的单炮数据的频谱的有效成分较接近；在30Hz-90Hz中主要是谐波噪声，按照本发明实施例的方法进行谐波噪声压制后的单炮数据中的谐波噪声得到了有效的压制。图5所示的频谱分析图表明，由于谐波能量较强，同样30Hz-90Hz频率之间谐波高频噪声比较突出。按照相关技术的方法进行谐波噪声压制后的单炮数据的频谱与按照本发明实施例的方法进行谐波噪声压制后的单炮数据的频谱的低频有效成分得到加强，同时10Hz-30Hz有效成分也得到加强。但在1Hz-30Hz内，本发明实施例的方法的保幅程度高于相关技术的方法；在30Hz-90Hz谐波噪声按照本发明实施例的方法得到了有效压制，从图4和图5中可以看出，与相关技术的方法相比，按照本发明实施例的方法进行谐波噪声压制后的单炮数据中的谐波噪声得到了有效的压制，并且按照本发明实施例的方法进行谐波噪声压制后的单炮数据的频谱的光滑型更高，同时，单频干扰压制要更好，整个频谱包络线也更加理想。图6中示出了谐波压制前的单炮数据的叠加剖面、按照本发明实施例的方法进行谐波噪声压制后的单炮数据的叠加剖面、以及叠加后的谐波噪声剖面示意图，从图中可以看出，压制前与压制后相比，按照本发明实施例的方法进行谐波噪声压制后的单炮数据的叠加剖面，其谐波噪声得到了有效的压制，信噪比提高。从主要勘探目的层成像效果改善明显，取得了令人满意的效果。

[0137] 本发明实施例还提供了一种谐波噪声的压制处理装置，如图7所示，该装置包括：

[0138] 第一识别模块702，用于利用波形过零点次数的谐波噪声识别法识别待处理的地震道数据的多个样点中含谐波噪声的样点，将识别出的上述含谐波噪声的样点作为谐波样点，并对上述谐波样点进行标记；

[0139] 第二识别模块704，用于利用分频的半波的平均振幅识别法确定各个频带中已标记的谐波样点是否需要压制，并对确定出的各个频带中需要压制的谐波样点进行谐波压制处理；

[0140] 数据重构模块706，用于对各个频带对应的压制后地震道数据进行重构，得到压制后的地震道数据。

[0141] 在本发明提供的实施例中采用两种谐波噪声的识别方式对地震高效采集带来的谐波噪声进行有效识别，再对识别出的谐波噪声进行了有效压制，从而提高了地震资料的信噪比和品质。

[0142] 进一步的，上述装置还包括：

[0143] 参数设置模块，用于设置谐波噪声的最小频率 并根据上述谐波噪声的最小频率确定与上述最小频率对应的上述谐波噪声的最大周期其中，fmin表示谐波噪声的最小频率，fbegin表示地震道正常子波振幅谱的开
始频率，fend表示地震道正常子波振幅谱的结束频率，fharm表示开始识别的谐波噪声频率，nzero表示谐波噪声的最大周期，单位为毫秒。

[0144] 进一步的，上述第一识别模块702包括：

[0145] 时窗创建单元，用于在待处理的地震道数据中,以切除点作为开始样点，并在该切除点处创建时窗长度为nwindows的时窗；

[0146] 第一确定单元，用于根据上述谐波噪声的最大周期和上述时窗长度确定上述时窗内谐波噪声的过零点个数nmin；

[0147] 第二确定单元，用于根据地震道数据中相邻的两个样点对应的振幅是否同号确定上述时窗内地震道数据的过零点个数ntrace；

[0148] 时窗移动单元，用于按照预设的时窗移动长度nmove将上述时窗向上述待处理的地震道数据的道结束点处方向移动；

[0149] 谐波噪声样点标记单元，用于判断上述时窗内地震道的过零点个数ntrace是否大于上述谐波噪声的过零点个数nmin，若是，则将移出的样点确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识；直到上述时窗的下边缘点移动至上述地震道数据的道结束点；判断上述时窗内地震道的过零点个数ntrace是否大于上述谐波噪声的过零点个数nmin，若是，则将上述时窗内的所有样点均确定为含有谐波噪声的样点，并对确定出的各个含有谐波噪声的样点进行标识。

[0150] 进一步的，上述第二识别模块704包括：

[0151] 频带划分单元，用于对上述待处理的地震道数据进行频带划分处理，以划分得到多个不同频率段的频带；

[0152] 标准频带选取单元，用于从划分得到的多个上述频带中选取一个频带作为标准频带，并将上述标准频带之后的多个频带作为待压制频带；

[0153] 当前待压制频带确定单元，用于从多个上述待压制频带中，选取上述标准频带之后的两个频带作为当前待压制频带；

[0154] 谐波样点识别与压制单元，用于在选取的两个上述当前待压制频带的已标记的谐波样点中，分别逐一判断两个上述当前待压制频带的各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与插值函数的乘积；若是，则对上述谐波样点的地震道数据进行谐波噪声的压制处理；将压制后的两个上述当前待压制频带中与上述标准频带相邻的频带作为下一个标准频带；直到当前选取的下一个标准频带为最后一个频带。

[0155] 进一步的，上述谐波样点识别与压制单元包括：

[0156] 浅层谐波样点识别子单元，用于当上述谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)，则判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与浅层插值函数p1(t)的乘积；

[0157] 中深层谐波样点识别子单元，用于当上述谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)，则判断各个上述谐波样点的上述当前待压制频带的半波的平均振幅是否大于上述标准频带的半波的平均振幅与中深层插值函数p2(t)的乘积。

[0158] 进一步的，上述浅层插值函数为 上述中深层插值函数为；

[0159] 其中，dtime表示从切除时刻开始压制的最小门槛值，etime表示距道结束时刻压制的最小门槛，ldelaytm表示地震道数据的浅层与中深层间的临界点，lmute表示地震道数据中的切除点，llength表示地震道数据中道结束点。

[0160] 进一步的，上述谐波样点识别与压制单元包括：

[0161] 浅层谐波样点压制子单元，用于当上述谐波样点位于地震道数据的浅层(lmute＜t＜ldelaytm)时，则根据公式 对与上述标准频带k-1相邻的上述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式对上述标准频带k-1之后的第二个上述待压制频带k+1中需要压制的
谐波样点t对应的原始地震道数据进行谐波压制处理；

[0162] 中深层谐波样点压制子单元，用于当上述谐波样点位于地震道数据的中深层(ldelaytm≤t≤llength)时，则根据公式 对与上述标准频带k-1相邻的上述当前待压制频带k中需要压制的谐波样点t对应的原始地震道数据进行谐波压制处理；并根据公式 对上述标准频带k-1之后的第二个上述待压制频带k+1中需要
压制的谐波样点t对应的原始地震道数据进行谐波压制处理；

[0163] 其中，t的取值为lmute至llength，yk-1(t)表示标准频带k-1中样点t的半波的平均振幅，yk(t)表示与标准频带k-1相邻的频带k中样点t的半波的平均振幅，xk(t)表示与标准频带k-1相邻的频带k中样点t的原始地震道数据，yk+1(t)表示标准频带k-1之后的第二个的频带k+1中样点t的半波的平均振幅，xk+1(t)表示标准频带k-1之后的第二个的频带k+1中样点t的原始地震道数据，γ表示振幅谱上振幅衰减系数。

[0164] 进一步的，上述数据重构模块706包括：

[0165] 数据重构单元，用于根据公式 确定压制后的地震道数据，其中，t表示样点号变量， 表示频带k中样点t对应的压制后地震道数据，nfnum表示划分得到的频带总数。

[0166] 在本发明提供的实施例中采用两种谐波噪声的识别方式对地震高效采集带来的谐波噪声进行有效识别，再对识别出的谐波噪声进行了有效压制，从而提高了地震资料的信噪比和品质；进一步的，在采用基于分频的谐波平均振幅识别法进行谐波噪声识别时，对高频频带中的谐波样点进行了两次基于分频的谐波平均振幅识别法的识别和压制，因而可以对两个频带的重叠部分中的样点进行有效的压制，从而取得更好的谐波噪声压制效果；更进一步的，考虑到谐波样点位于地震道数据的浅层与中深层时的谐波噪声干扰能力存在一定的差异，因而，在比较半波的平均振幅时谐波样点位于地震道数据的浅层或者中深层引入不同的差值函数，提高了对谐波噪声进行识别与压制的准确度。

[0167] 本发明的谐波噪声的压制处理装置中包括的各个模块可以为设备上的特定硬件或者安装于设备上的软件或固件等。所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，前述描述的模块和单元的具体工作过程，均可以参考上述方法实施例中的对应过程。在本申请所提供的几个实施例中，应该理解到，所揭露的谐波噪声的压制处理方法及装置，可以通过其它的方式实现。以上所描述的谐波噪声的压制处理装置的实施例仅仅是示意性的，例如，上述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

[0168] 所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。

[0169] 所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

[0170] 以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。