CN106125132A - 含单频干扰地震道的迭代识别和压制方法 - Google Patents

Abstract

本发明涉及一种含单频干扰地震道的迭代识别和压制方法。主要解决了无法精确估算单频干扰的频率、无法有效精确识别待处理地震道是否含有单频干扰等问题。该方法能够逐道自动识别待处理地震道是否含有单频干扰，不需要人工参与识别，不受地震道处理域的限制，而且能够精确估算单频干扰的振幅、相位及其频率信息，从而实现单频干扰的有效压制。

Description

含单频干扰地震道的迭代识别和压制方法

技术领域

本发明属于油气勘探中的地震资料数字处理领域，主要针对油田勘探成熟区原始采集地震资料中的噪声压制处理，具体一种含单频干扰地震道的一种迭代识别和自适应压制方法。

背景技术

在油气勘探程度较为成熟的地区，密集分布的电网对布置在附近的地震检波器形成电磁场感应，从而在接收到的地震资料中叠加一种或多种频率成分的单频干扰，即工业电干扰，且从浅层到深层其振幅、相位和频率基本保持稳定不变，同时这些频率成分的单频信号与地震资料的优势频带范围相重叠，因此，高保真、有效地压制这种干扰信号是当前叠前地震资料预处理的重要环节。

早期采用频域陷波方法来压制单频干扰，其结果将损失限制频率附近的有效信号频率成分，同时选取参数不当也会引入吉布斯现象和边界效应等问题。为此，地球物理工作者们提出了多种时间域或频率域的单频干扰压制方法，如基于维纳滤波的工业电干扰去噪方法、线性调频谱法消除工业电干扰、余弦逼近工业电干扰的压制方法等，调研可知，上述文献方法只能对已知含有单频干扰的地震道进行单频干扰的识别和压制，并要求单频干扰的频率基本固定，而实际单频干扰的频率成分是在一定范围内变化的。

另外，业界普遍采用共检波点域均方根振幅能量法来识别待处理地震道是否含有单频干扰，但这种方法需人工参与识别、计算流程复杂，且受单一共检波点道集内不同地震道能量的累加效应影响，对识别和压制含较弱能量的单频干扰无能为力，经能量补偿后这些地震道中的单频干扰噪声将被放大，从而影响后续地震预处理道集的品质。

发明内容

针对上述背景技术无法精确估算单频干扰的频率、无法有效精确识别待处理地震道是否含有单频干扰等问题，本发明公开了一种时间域寻优的含单频干扰地震道迭代识别和压制方法，该方法能够逐道自动识别待处理地震道是否含有单频干扰，不需要人工参与识别，不受地震道处理域的限制，而且能够精确估算单频干扰的振幅、相位及其频率信息，从而实现单频干扰的有效压制。

技术方案：一种含单频干扰地震道的迭代识别和压制方法，其步骤如下：

Ⅰ、结合地表电缆分布情况调查和原始地震数据频谱分析确定处理工区单频干扰初始频率f₀。

Ⅱ、建立单频干扰的离散数学模型。以余弦函数为例:

y_i＝C cos[2πfiΔt-α]，其中i＝1,2,…,M (1)

式(1)中，C和α分别为单频干扰的振幅和相位，Δt为时间采样间隔，i代表采样点序号，M为地震道的总样点个数；

为避免同时求解式(1)中与y_i存在非线性关系的频率f和相位α，将式(1)进一步分解为正余弦函数的加权组合形式：

y_i＝A cos(2πfiΔt)+B sin(2πfiΔt) (2)

式(2)中，A和B分别为余弦和正弦函数的系数，且均与y_i成线性关系。

Ⅲ、建立目标函数。

根据待处理地震道和单频干扰模型的残差能量建立目标函数Q：

$Q=\underset{i=1}{\overset{M}{\Σ}}{\[S_i-y_i\]}^2---\left(3\right)$

Ⅳ、根据最小二乘法建立数值离散公式，使得待处理地震道和单频干扰模型的残差能量实现最小化。

经换算得到如下离散计算公式：

$A=\frac{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\sin }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)-\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\cos }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\sin }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)-{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}^2}---\left(5\right)$

$B=\frac{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)-\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\cos }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}{{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}^2-\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\cos }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\sin }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}---\left(6\right)$

根据式(2)可进一步推得式(1)中单频干扰的振幅C和相位α如下：

和α＝arctan(B/A)＝2πfτΔt (7)

可以得到：

根据式(5)～式(7)可知，在给定单频干扰的频率f情况下，可以精确获得单频干扰的振幅C和相位α；

Ⅴ、根据牛顿迭代法寻优逼近单频干扰频率的精确解。

根据式(8)可知无法得到频率f的解析表达式，因此，引入牛顿迭代法来解决频率精确解的逼近问题。

式(9)中，N为牛顿迭代法的迭代次数。根据前后两次迭代的频率误差小于或等于给定误差的范围时，即满足|f_n+1-f_n|≤ξ，终止迭代计算。式(9)经推导，得到如下公式：

$f_{n+1}=f_n-\frac{\underset{i=1}{\overset{M}{\Σ}}\[2{\mathrm{iS}}_{i}sin(2{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-{\mathrm{\α}}_n)\]-\underset{i=1}{\overset{M}{\Σ}}\[{\mathrm{iC}}_{n}sin(4{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-2{\mathrm{\α}}_n)\]}{\underset{i=1}{\overset{M}{\Σ}}\[4{\mathrm{\πi}}^2{S}_i\mathrm{\Δ}t\cos (2{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-{\mathrm{\α}}_n)\]-\underset{i=1}{\overset{M}{\Σ}}\[4{\mathrm{\πi}}^2{C}_n\mathrm{\Δ}tcos(4{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-2{\mathrm{\α}}_n)\]}---\left(10\right)$

根据式(10)可知，单频干扰的频率可以通过迭代逐渐收敛到局部最小值。

Ⅵ、建立含单频干扰地震道的识别准则。

根据待处理地震道和预测得到的单频干扰模型进行归一化互相关计算，根据归一化互相关系数ψ及给定的单频干扰识别门槛值Ψ，实现含单频干扰地震道的逐道自动识别。

归一化互相关系数：

式(11)中ψ的值域范围为[0,1]。当归一化互相关系数ψ小于给定的门槛值Ψ时，自动认为该地震道不含单频干扰，反之亦然，这里的Ψ一般取值为0.2。

通过迭代寻优，实现含单频干扰地震道的自动识别与自适应压制具体方法

的处理步骤如下：

Ⅰ、根据单频干扰的初始频率f₀、式(5)和式(6)，再结合式(7)，计算得到单频干扰初始模型的振幅C₀和相位α₀；

Ⅱ、根据式(11)计算得到归一化互相关系数ψ₀，并与单频干扰识别门槛值Ψ进行比较，若ψ₀>Ψ，则认为该地震道含有单频干扰，进入下一步骤，否则，认为该地震道不含单频干扰，因此，不做任何处理，可结束该地震道的识别和压制；

Ⅲ、根据式(8)进行最优化迭代计算，得到迭代更新后的单频干扰频率f₁；

Ⅳ、根据条件判断|f_n+1-f_n|≤ξ，如满足该判断条件，则根据步骤Ⅰ计算得到单频干扰最终模型的振幅C和相位α，并将预测得到的单频干扰模型道从待处理地震道中减去，即x_i＝S_i-y_i，得到经自适应压制单频干扰后的高保真地震道；如不满足该判断条件，则根据计算得到单频干扰模型的振幅C₁和相位α₁，再根据步骤Ⅱ和步骤Ⅲ完成迭代计算，直到满足条件为止；

Ⅴ、对于多种频率成分单频干扰的迭代识别和压制问题，重新设置单频干扰的初始频率，进行迭代处理即可，也可以结合多核CPU并行处理能够有效提高处理效率。

有益效果是：通过采用本发明的技术方案，该方法能够逐道自动识别待处理地震道是否含有单频干扰，不需要人工参与识别，不受地震道处理域的限制，而且能够精确估算单频干扰的振幅、相位及其频率信息，从而实现单频干扰的有效压制，因此，本发明方法可为当前地震资料高保真预处理提供重要的理论和方法参考。

附图说明

图1a为本发明由不同振幅、频率和相位合成的单频干扰模型数据。

图1b为本发明单频干扰模型数据对应的单频干扰振幅、频率和相位变化示意图。

图1c本发明方法自动识别和压制掉的单频干扰数据。

图1d本发明方法压制单频干扰后的剩余数据。

图2a本发明实际含单频干扰的共炮点记录及其频谱。

图2b本发明方法压制单频干扰后的共炮点记录及其频谱。

图2c本发明方法压制掉的单频干扰及其频谱。

图2d本发明方法用于识别含单频干扰地震道的归一化自相关函数。

图2e本发明方法识别出含单频干扰地震道的单频干扰频率。

图2f本发明方法识别出含单频干扰地震道的单频干扰振幅。

图2g本发明方法识别出含单频干扰地震道的单频干扰相位。

图3a实际含单频干扰的共检波点记录及其频谱。

图3b本发明方法压制单频干扰后的共检波点记录及其频谱。

图3c本发明方法压制掉的单频干扰及其频谱。

具体实施方式

下面结合附图对本发明做进一步说明。

1、一种含单频干扰地震道的迭代识别和压制方法，其步骤如下：

Ⅰ、结合地表电缆分布情况调查和原始地震数据频谱分析确定处理工区单频干扰初始频率f₀。

Ⅱ、建立单频干扰的离散数学模型。以余弦函数为例:

y_i＝C cos[2πfiΔt-α]，其中i＝1,2,…,M (1)

式(1)中，C和α分别为单频干扰的振幅和相位，Δt为时间采样间隔，i代表采样点序号，M为地震道的总样点个数；

为避免同时求解式(1)中与y_i存在非线性关系的频率f和相位α，将式(1)进一步分解为正余弦函数的加权组合形式：

y_i＝A cos(2πfiΔt)+B sin(2πfiΔt) (2)

式(2)中，A和B分别为余弦和正弦函数的系数，且均与y_i成线性关系。

Ⅲ、建立目标函数。

根据待处理地震道和单频干扰模型的残差能量建立目标函数Q：

$Q=\underset{i=1}{\overset{M}{\Σ}}{\[S_i-y_i\]}^2---\left(3\right)$

Ⅳ、根据最小二乘法建立数值离散公式，使得待处理地震道和单频干扰模型的残差能量实现最小化。

经换算得到如下离散计算公式：

$A=\frac{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\sin }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)-\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\cos }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\sin }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)-{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}^2}---\left(5\right)$

$B=\frac{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)-\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\cos }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}{{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}^2-\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\cos }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\sin }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}---\left(6\right)$

根据式(2)可进一步推得式(1)中单频干扰的振幅C和相位α如下：

和α＝arctan(B/A)＝2πfτΔt (7)

可以得到：

根据式(5)～式(7)可知，在给定单频干扰的频率f情况下，可以精确获得单频干扰的振幅C和相位α；

Ⅴ、根据牛顿迭代法寻优逼近单频干扰频率的精确解。

根据式(8)可知无法得到频率f的解析表达式，因此，引入牛顿迭代法来解决频率精确解的逼近问题。

式(9)中，N为牛顿迭代法的迭代次数。根据前后两次迭代的频率误差小于或等于给定误差的范围时，即满足|f_n+1-f_n|≤ξ，终止迭代计算。式(9)经推导，得到如下公式：

$f_{n+1}=f_n-\frac{\underset{i=1}{\overset{M}{\Σ}}\[2{\mathrm{iS}}_{i}sin(2{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-{\mathrm{\α}}_n)\]-\underset{i=1}{\overset{M}{\Σ}}\[{\mathrm{iC}}_{n}sin(4{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-2{\mathrm{\α}}_n)\]}{\underset{i=1}{\overset{M}{\Σ}}\[4{\mathrm{\πi}}^2{S}_i\mathrm{\Δ}t\cos (2{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-{\mathrm{\α}}_n)\]-\underset{i=1}{\overset{M}{\Σ}}\[4{\mathrm{\πi}}^2{C}_n\mathrm{\Δ}tcos(4{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-2{\mathrm{\α}}_n)\]}---\left(10\right)$

根据式(10)可知，单频干扰的频率可以通过迭代逐渐收敛到局部最小值；

Ⅵ、建立含单频干扰地震道的识别准则；

根据待处理地震道和预测得到的单频干扰模型进行归一化互相关计算，根据归一化互相关系数ψ及给定的单频干扰识别门槛值Ψ，实现含单频干扰地震道的逐道自动识别；

归一化互相关系数：

式(11)中ψ的值域范围为[0,1]。当归一化互相关系数ψ小于给定的门槛值Ψ时，自动认为该地震道不含单频干扰，反之亦然，这里的Ψ一般取值为0.2；

通过迭代寻优，实现含单频干扰地震道的自动识别与自适应压制具体方法

的处理步骤如下：

Ⅰ、根据单频干扰的初始频率f₀、式(5)和式(6)，再结合式(7)，计算得到单频干扰初始模型的振幅C₀和相位α₀；

Ⅱ、根据式(11)计算得到归一化互相关系数ψ₀，并与单频干扰识别门槛值Ψ进行比较，若ψ₀>Ψ，则认为该地震道含有单频干扰，进入下一步骤，否则，认为该地震道不含单频干扰，因此，不做任何处理，可结束该地震道的识别和压制；

Ⅲ、根据式(8)进行最优化迭代计算，得到迭代更新后的单频干扰频率f₁；

Ⅳ、根据条件判断|f_n+1-f_n|≤ξ，如满足该判断条件，则根据步骤Ⅰ计算得到单频干扰最终模型的振幅C和相位α，并将预测得到的单频干扰模型道从待处理地震道中减去，即x_i＝S_i-y_i，得到经自适应压制单频干扰后的高保真地震道；如不满足该判断条件，则根据计算得到单频干扰模型的振幅C₁和相位α₁，再根据步骤Ⅱ和步骤Ⅲ完成迭代计算，直到满足条件为止；

Ⅴ、对于多种频率成分单频干扰的迭代识别和压制问题，重新设置单频干扰的初始频率，进行迭代处理即可，也可以结合多核CPU并行处理能够有效提高处理效率。

由图1a‐图1d可知，本发明方法适用于含不同频率、振幅和相位情况下的合成模型数据的单频干扰压制处理，且能够精确识别单频干扰的频率、振幅和相位信息，因此本发明方法处理结果具有较高的信噪比，不存在边界效应等问题。

根据图2a‐图2g可知，本发明方法可以处理含一种单频干扰频率成分情况下的实际共炮点道集地震数据如图2a。通过地震数据频谱分析可知，单频干扰主要集中在50Hz附近如图2b，于是选取50Hz作为单频干扰初始频率，由此计算得到图2d所示的归一化互相关系数，并与门槛值进行对比，本例中选取0.2；对识别出含单频干扰的地震道作进一步的单频干扰振幅、频率和相位的寻优逼近迭代计算，最终得到的单频干扰的振幅、频率和相位如图2e、图2f、图2g，并实现自适应压制。从图2e可知，识别出的单频干扰频率在50Hz附近一定范围内变化，从图2b和图2c可知，本发明方法能够有效压制单频干扰，识别出的单频干扰数据中不存在任何的有效信号，同时对含异常值的地震道不做任何处理，有效信号频谱特征得到保持。

根据图3a‐图3c可知，本发明方法可以处理含多种单频干扰频率成分情况下的实际共检波点道集地震数据，由于本发明方法是逐道处理，因此不受处理域的限制，如图3a可知通过地震数据频谱分析可知，单频干扰主要集中在50Hz、150Hz、250Hz、350Hz附近，于是分别选取50Hz、150Hz、250Hz、350Hz作为单频干扰的初始频率，由此计算相应的归一化互相关系数，本例中选取0.2，并与门槛值进行对比，对识别出含单频干扰的地震道作进一步的单频干扰振幅、频率和相位的寻优逼近迭代计算，并实现自适应压制。从图3b和图3c可知，本发明方法能够有效压制多种频率成分的单频干扰，道集信噪比明显得到改善，有效信号频谱特征得到保持。从图3c还可知，本发明方法识别出的单频干扰道集数据中不存在任何有效信号，从频谱上分析，单频干扰的频率主要集中在50Hz、150Hz、250Hz、350Hz位置附近。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (2)

1.一种含单频干扰地震道的迭代识别和压制方法，其步骤如下：

Ⅰ、结合地表电缆分布情况调查和原始地震数据频谱分析确定处理工区单频干扰初始频率f₀。

Ⅱ、建立单频干扰的离散数学模型。以余弦函数为例:

y_i＝C cos[2πfiΔt-α]，其中i＝1,2,…,M (1)

式(1)中，C和α分别为单频干扰的振幅和相位，Δt为时间采样间隔，i代表采样点序号，M为地震道的总样点个数；

为避免同时求解式(1)中与y_i存在非线性关系的频率f和相位α，将式(1)进一步分解为正余弦函数的加权组合形式：

y_i＝A cos(2πfiΔt)+B sin(2πfiΔt) (2)

式(2)中，A和B分别为余弦和正弦函数的系数，且均与y_i成线性关系。

Ⅲ、建立目标函数。

根据待处理地震道和单频干扰模型的残差能量建立目标函数Q：

$Q=\underset{i=1}{\overset{M}{\Σ}}{\[S_i-y_i\]}^2---\left(3\right)$

Ⅳ、根据最小二乘法建立数值离散公式，使得待处理地震道和单频干扰模型的残差能量实现最小化。

其中

(4)

经换算得到如下离散计算公式：

$A=\frac{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\sin }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)-\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\cos }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\sin }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)-{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}^2}---\left(5\right)$

$B=\frac{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)-\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\cos }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{S}_i\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}{{\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\cos \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\sin \left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}^2-\left(\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\cos }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\×\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\sin }^2\left(2\mathrm{\π}fi\mathrm{\Δ}t\right)\right)}---\left(6\right)$

根据式(2)可进一步推得式(1)中单频干扰的振幅C和相位α如下：

和α＝arctan(B/A)＝2πfτΔt (7)

可以得到：

根据式(5)～式(7)可知，在给定单频干扰的频率f情况下，可以精确获得单频干扰的振幅C和相位α；

V、根据牛顿迭代法寻优逼近单频干扰频率的精确解。

根据式(8)可知无法得到频率f的解析表达式，因此，引入牛顿迭代法来解决频率精确解的逼近问题。

式(9)中，N为牛顿迭代法的迭代次数。根据前后两次迭代的频率误差小于或等于给定误差的范围时，即满足|f_n+1-f_n|≤ξ，终止迭代计算。式(9)经推导，得到如下公式：

$f_{n+1}=f_n-\frac{\underset{i=1}{\overset{M}{\Σ}}\[2{\mathrm{iS}}_{i}sin(2{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-{\mathrm{\α}}_n)\]-\underset{i=1}{\overset{M}{\Σ}}\[{\mathrm{iC}}_{n}sin(4{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-2{\mathrm{\α}}_n)\]}{\underset{i=1}{\overset{M}{\Σ}}\[4{\mathrm{\πi}}^2{S}_i\mathrm{\Δ}t\cos (2{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-{\mathrm{\α}}_n)\]-\underset{i=1}{\overset{M}{\Σ}}\[4{\mathrm{\πi}}^2{C}_n\mathrm{\Δ}tcos(4{\mathrm{\πf}}_{n}i\mathrm{\Δ}t-2{\mathrm{\α}}_n)\]}---\left(10\right)$

根据式(10)可知，单频干扰的频率可以通过迭代逐渐收敛到局部最小值。

VI、建立含单频干扰地震道的识别准则。

根据待处理地震道和预测得到的单频干扰模型进行归一化互相关计算，根据归一化互相关系数ψ及给定的单频干扰识别门槛值Ψ，实现含单频干扰地震道的逐道自动识别。

归一化互相关系数：

式(11)中ψ的值域范围为[0，1]。当归一化互相关系数ψ小于给定的门槛值Ψ时，自动认为该地震道不含单频干扰，反之亦然，这里的Ψ一般取值为0.2。

2.根据权利要求1所述的一种含单频干扰地震道的迭代识别和压制方法，其特征在于：通过迭代寻优，实现含单频干扰地震道的自动识别与自适应压制具体方法的处理步骤如下：

I、根据单频干扰的初始频率f₀、式(5)和式(6)，再结合式(7)，计算得到单频干扰初始模型的振幅C₀和相位α₀；

II、根据式(11)计算得到归一化互相关系数ψ₀，并与单频干扰识别门槛值Ψ进行比较，若ψ₀＞Ψ，则认为该地震道含有单频干扰，进入下一步骤，否则，认为该地震道不含单频干扰，因此，不做任何处理，可结束该地震道的识别和压制；

III、根据式(8)进行最优化迭代计算，得到迭代更新后的单频干扰频率f₁；

IV、根据条件判断|f_n+1-f_n|≤ξ，如满足该判断条件，则根据步骤1计算得到单频干扰最终模型的振幅C和相位α，并将预测得到的单频干扰模型道从待处理地震道中减去，即x_i＝S_i-y_i，得到经自适应压制单频干扰后的高保真地震道；如不满足该判断条件，则根据计算得到单频干扰模型的振幅C₁和相位α₁，再根据步骤2和步骤3完成迭代计算，直到满足条件为止；

V、对于多种频率成分单频干扰的迭代识别和压制问题，重新设置单频干扰的初始频率，进行迭代处理即可，也可以结合多核CPU并行处理能够有效提高处理效率。