技术领域
[0001] 本
发明属于
地震资料处理领域，具体涉及一种可控震源次谐波噪声的压制方法及系统。
背景技术
[0002] 可控震源勘探过程中，由于震源机械和震动装置的非线性，以及震板与大地的耦合问题，使得可控震源输出的
信号产生畸变。由于畸变信号除了产生震源期望的基波
频率外，还会产生谐波或次谐波频率。这些谐波在基波和谐波频率相关过程中，由于基波频率和谐波频率的相关性，会在相关后产生强周期性噪声。次谐波噪声通常作为混响存在于可控震源震地震数据中，这些噪声表现为与主要事件平行且在主要事件之前或之后产生的横向相干带。次谐波噪声与有效的基波信号频率重叠部分，采用线性升频扫描互相关后会将次谐波噪声压缩到较晚时间区域内，影响深层反射波
信噪比。
[0003] 国内外学者提出的压制方法主要针对高次谐波干扰。Li等(1995)针对线性扫描方式提出了纯
相移滤波方法，可有效地压制近场和远场的谐波畸变，它是针对传统的单炮点采集提出的，不适用于可控震源滑动扫描采集方式。针对滑动扫描数据，Fleure(2002)提出了一种在频率域分离谐波并对其压制的技术，Meunier和Bianchi等(2002)提出了一种通过预测谐波，再在每道记录中依次减去谐波的方法。Moerig(2007)提出通过对相关后数据中负时间域的谐波噪声进行统计，来估计地震记录中谐波
能量，进而消除谐波畸变的方法。Abd El-Aal(2010)和Claudio Bagaini(2010)分别在时间域和频率域推导谐波预测算子来预测相关后地震记录中的谐波噪声，但它需要用到相关前数据，无法直接处理相关后数据。
Charles Sicking(2009)利用地面
力信号和扫描起始时间设计
滤波器，将滤波器与基波能量褶积来估计该炮对上一炮记录产生的谐波噪声，该方法可直接用来处理相关后记录。
[0004] 中国
专利公开文献CN201510278108.1公开了一种基于自适应匹配滤波算子的可控震源谐波压制方法；中国专利公开文献CN201510278789.1公开了一种基于预测滤波法和纯相移法的可控震源谐波压制方法；中国专利公开文献CN201210563556.2公开了一种压制可控震源谐波干扰的方法；中国公开文献“Hilbert-Huang变换去除可控震源谐波畸变”(清华大学学报(自然科版)，2011年06期)；中国公开文献“利用可控震源的力信号压制谐波干扰”(石油物探，2011年01期)；中国公开文献“自适应匹配预测滤波压制可控震源谐波”(石油地球物理勘探，2016年06期)。目前公开文献中的方法主要是用于压制可控震源资料中的高阶谐波噪声，对于低次谐波噪声压制方法涉及较少。
发明内容
[0005] 本发明的目的在于解决上述
现有技术中存在的难题，提供一种可控震源次谐波噪声的压制方法及系统，将次谐波噪声变换到其它域进行压制，在有效保留原始
地震道有效信号的情况下，消除震源信号中的次谐波干扰噪声。
[0006] 本发明是通过以下技术方案实现的：
[0007] 一种可控震源次谐波噪声的压制方法，所述方法将含次谐波噪声的原始地震记录从空间-时间域变换到次谐波域，使得次谐波噪声与震源基波信号分离开，对次谐波噪声进行
切除，再从次谐波域反变换到空间-时间域，得到压制次谐波噪声的相关前记录。
[0008] 所述方法包括：
[0009] 步骤S1：输入含次谐波噪声的原始地震记录sc(t)，并根据扫描信号参数生成次谐波信号qc(t)；
[0010] 步骤S2：将次谐波信号qc(t)与含次谐波噪声的原始地震记录sc(t)进行互相关处理，得到变换到次谐波域的相关记录；
[0011] 步骤S3：利用面波时差对步骤S2得到的变换到次谐波域的相关记录进行线性时差校正，经过线性时差校正后的次谐波噪声集中在一个
水平时窗内；
[0012] 步骤S4：对步骤S3得到的集中在水平时窗内的次谐波噪声开取时窗，对时窗内的次谐波噪声进行置零，得到压制噪声后的地震记录；
[0013] 步骤S5：利用面波时差将步骤S4得到的压制噪声后的地震记录进行反线性时差校正得到反线性时差校正后的地震记录，然后利用次谐波信号qc(t)对所述反线性时差校正后的地震记录进行反相关运算，将其从次谐波域变换到空间-时间域，得到压制次谐波噪声的相关前记录。
[0014] 所述步骤S1中输入的含次谐波噪声的原始地震记录sc(t)是实际采集到的地震记录，对其假设如下：
[0015] sc(t)＝r(t)*(q(t)+qc(t))，
[0016] 其中，qc(t)为震源基波信号，qc(t)为次谐波信号，r(t)为反射系数；
[0017] 所述步骤S1中的根据扫描信号参数生成次谐波信号qc(t)的操作包括：
[0018] 所述扫描信号参数包括：起始频率、终止频率和扫描长度；
[0019] 将所述起始频率的一半作为次谐波信号的起始频率，将所述终止频率的一半作为次谐波信号的终止频率；
[0020] 将所述扫描长度作为次谐波信号的扫描长度；
[0021] 根据所述次谐波信号的起始频率、次谐波信号的终止频率、扫描长度生成次谐波信号。
[0022] 所述步骤S2的操作包括：
[0023] 利用下式进行互相关处理：
[0024] sc(t)*qc(-t)＝(r(t)*q(t)+r(t)*qc(t))*qc(-t)
[0025] ＝r(t)*q(t)*qc(-t)+r(t)*qc(t)*qc(-t)
[0026] ＝r(t)*ccf{q(t)，qc(t)}+r(t)*acf{qc(t)}
[0027] 其中，qc(-t)为低次谐波，ccf{q(t)，qc(t)}为震源基波信号和次谐波信号的互相关，acf{qc(t)}为次谐波信号的自相关。
[0028] 所述步骤S3中的利用面波时差对步骤S2得到的变换到次谐波域的相关记录进行线性时差校正的操作包括：
[0029] 根据不同接收道的偏移距和面波视速度，计算各道面波的到达时间，利用各道面波的到达时间计算得到面波时差，利用面波时差对步骤S2得到的变换到次谐波域的相关记录做线性时差校正。
[0030] 所述步骤S4中的时窗的宽度为次谐波噪声的一个
波长的宽度。
[0031] 一种可控震源次谐振噪声的压制系统，包括：
[0032] 输入单元：用于输入含次谐波噪声的原始地震记录sc(t)，并根据扫描信号参数生成次谐波信号qc(t)；
[0033] 互相关单元：与输入单元连接，用于将次谐波信号qc(t)与含次谐波噪声的原始地震记录sc(t)进行互相关处理，得到变换到次谐波域的相关记录，并将其发送给线性时差校正单元；
[0034] 线性时差校正单元：与互相关单元连接，利用面波时差对互相关单元发送来的变换到次谐波域的相关记录进行线性时差校正，经过线性时差校正后的次谐波噪声集中在一个水平时窗内；
[0035] 切除噪声单元：与线性时差校正单元连接，对线性时差校正单元得到的集中在水平时窗内的次谐波噪声开取时窗，对时窗内的次谐波噪声进行置零，得到压制噪声后的地震记录，并将其发送给输出单元；
[0036] 输出单元：与切除噪声单元连接，利用面波时差将切除噪声单元发送来的压制噪声后的地震记录进行反线性时差校正得到反线性时差校正后的地震记录，然后利用次谐波信号qc(t)对所述反线性时差校正后的地震记录进行反相关运算，将其从次谐波域变换到空间-时间域，得到压制次谐波噪声的相关前记录，并将其输出。
[0037] 所述输入单元对输入的含次谐波噪声的原始地震记录sc(t)假设如下：
[0038] sc(t)＝r(t)*(q(t)+qc(t))，
[0039] 其中，q(t)为震源基波信号，qc(t)为次谐波信号，r(t)为反射系数；
[0040] 所述扫描信号参数包括：起始频率、终止频率和扫描长度；所述输入单元将所述起始频率的一半作为次谐波信号的起始频率，将所述终止频率的一半作为次谐波信号的终止频率；将所述扫描长度作为次谐波信号的扫描长度；根据所述次谐波信号的起始频率、次谐波信号的终止频率、扫描长度生成次谐波信号qc(-t)。
[0041] 所述互相关单元利用下式进行互相关处理：
[0042] sc(t)*qc(-t)＝(r(t)*q(t)+r(t)*qc(t))*qc(-t)
[0043] ＝r(t)*q(t)*qc(-t)+r(t)*qc(t)*qc(-t)
[0044] ＝r(t)*ccf{q(t)，qc(t)}+r(t)*acf{qc(t)}
[0045] 其中，qc(-t)为低次谐波，ccf{q(t)，qc(t)}为震源基波信号和次谐波信号的互相关，acf{qc(t)}为次谐波信号的自相关。
[0046] 所述切除噪声单元开取的时窗的宽度为次谐波噪声的一个波长的宽度。
[0047] 本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机可执行的至少一个程序，所述至少一个程序被所述计算机执行时使所述计算机执行可控震源次谐振噪声的压制方法中的步骤。
[0048] 与现有技术相比，本发明的有益效果是：本发明提出的可控震源次谐波噪声压制方法，将噪
声波列变换到“次谐波域”，使得谐波噪声成分集中起来，并与基波信号分离开，使得谐波信号更容易被检测到并得到压制，在有效保留原始地震道有效信号的情况下，消除了震源信号中的次谐波干扰噪声。试验结果表明：该方法对模拟资料中低次谐波噪声进行消除，效果较好，可进一步应用到实际资料处理中。
附图说明
[0049] 图1含次谐波干扰的原始相关前记录；
[0050] 图2次谐波信号与原始记录的互相关记录；
[0051] 图3线性时差后记录；
[0052] 图4切除干扰波能量后的记录；
[0053] 图5反线性时差后记录；
[0054] 图6压制次谐波噪声后的相关前记录；
[0055] 图7a次谐波噪声压制前的相关记录；
[0056] 图7b次谐波噪声压制后的相关记录；
具体实施方式
[0058] 下面结合附图对本发明作进一步详细描述：
[0059] 本发明提供了一种次谐波噪声压制方法，通过设计次谐波算子，将含噪声记录从空间-时间域变换到“次谐波域”，利用面波视速度进行线性时差，次谐波同相轴能量聚焦成像，基波同相轴能量分散，开取时窗对识别出的次谐波干扰切除，再反变换到空间-时间域，得到压制次谐波噪声后的记录。
[0060] 本发明方法的基本思路是利用扫描信号参数(扫描起始/终止频率、扫描长度)设计次谐波信号，把次谐波信号作为滤波器算子，与含次谐波噪声信号作互相关处理，将含噪声记录从空间-时间域变换到“次谐波域”，此时次谐波干扰同相轴能量聚焦成像，有效信号同相轴能量分散，由此能够选择性识别出次谐波干扰并进行切除，再反变换到空间-时间域，得到压制次谐波噪声后的记录。
[0061] 假设：原始未相关记录sc(t)＝r(t)*(q(t)+qc(t))，q(t)为震源基波信号，qc(t)为次谐波信号，r(t)为反射系数，sc(t)为含次谐波噪声的地震记录，是实际采集到的地震信号，其包括了震源基本信号和次谐波噪声。
[0062] 如图8所示，本发明可控震源次谐波噪声的压制方法的具体实现步骤为：
[0063] 步骤S1：根据扫描信号参数(起始/终止频率，扫描长度)生成次谐波信号qc(t)(次谐波信号的起始/终止频率为扫描信号的起始/终止频率的1/2,扫描长度、
采样率与扫描信号的相同)。输入含次谐波噪声的地震记录sc(t)和次谐波信号qc(t)。步骤S1得到的次谐波信号包括了高次谐波和低次谐波，本发明方法只对低次谐波进行处理。
[0064] 步骤S2：将次谐波信号qc(t)与原始记录sc(t)进行互相关处理，得到变换到次谐波域的相关记录。
[0065] sc(t)*qc(-t)＝(r(t)*q(t)+r(t)*qc(t))*qc(-t)
[0066] ＝r(t)*q(t)*qc(-t)+r(t)*qc(t)*qc(-t)
[0067] ＝r(t)*ccf{q(t)，qc(t)}+r(t)*acf{qc(t)}
[0068] 其中，qc(-t)为低次谐波，ccf{q(t)，qc(t)}为震源基波信号和次谐波信号的互相关，acf{qc(t)}为次谐波信号的自相关。
[0069] 步骤S3：利用面波视速度对变换到次谐波域的相关记录进行线性时差校正。具体如下：根据不同接收道的偏移距值和面波视速度，计算各道面波达到时间，然后计算得到面波时差，利用面波时差对次谐波域道集做线性时差校正，线性时差校正后的次谐波噪声能量会集中在一个较窄水平时窗内。
[0070] 步骤S4：对步骤S3中集中在较窄时窗内的次谐波噪声能量开取时窗，时窗大小可选为次谐波噪声一个波长宽度，对时窗内噪声能量进行置零，达到压制噪声目的。
[0071] 步骤S5：再利用所述面波时差，将压制噪声后的数据进行反线性时差(利用现有的反线性时差校正方法实现的)，并利用次谐波信号qc(t)对压噪数据进行反相关运算，将次谐波域记录变换到x-t域(空间-时间域)，得到压制次谐波噪声的相关前记录。
[0072] 一种可控震源次谐振噪声的压制系统，包括：
[0073] 输入单元：用于输入含次谐波噪声的原始地震记录sc(t)，并根据扫描信号参数生成次谐波信号qc(t)；
[0074] 互相关单元：与输入单元连接，用于将次谐波信号qc(t)与含次谐波噪声的原始地震记录sc(t)进行互相关处理，得到变换到次谐波域的相关记录，并将其发送给线性时差校正单元；
[0075] 线性时差校正单元：与互相关单元连接，利用面波时差对互相关单元发送来的变换到次谐波域的相关记录进行线性时差校正，经过线性时差校正后的次谐波噪声集中在一个水平时窗内；
[0076] 切除噪声单元：与线性时差校正单元连接，对线性时差校正单元得到的集中在水平时窗内的次谐波噪声开取时窗，对时窗内的次谐波噪声进行置零，得到压制噪声后的地震记录，并将其发送给输出单元；
[0077] 输出单元：与切除噪声单元连接，利用面波时差将切除噪声单元发送来的压制噪声后的地震记录进行反线性时差校正得到反线性时差校正后的地震记录，然后利用次谐波信号qc(t)对所述反线性时差校正后的地震记录进行反相关运算，将其从次谐波域变换到空间-时间域，得到压制次谐波噪声的相关前记录，并将其输出。
[0078] 所述输入单元对输入的含次谐波噪声的原始地震记录sc(t)假设如下：
[0079] sc(t)＝r(t)*(q(t)+qc(t))，
[0080] 其中，q(t)为震源基波信号，qc(t)为次谐波信号，r(t)为反射系数；
[0081] 所述扫描信号参数包括：起始频率、终止频率和扫描长度；所述输入单元将所述起始频率的一半作为次谐波信号的起始频率，将所述终止频率的一半作为次谐波信号的终止频率；将所述扫描长度作为次谐波信号的扫描长度；根据所述次谐波信号的起始频率、次谐波信号的终止频率、扫描长度生成次谐波信号qc(t)。
[0082] 所述互相关单元利用下式进行互相关处理：
[0083] sc(t)*qc(-t)＝(r(t)*q(t)+r(t)*qc(t))*qc(-t)
[0084] ＝r(t)*q(t)*qc(-t)+r(t)*qc(t)*qc(-t)
[0085] ＝r(t)*ccf(q(t)，qc(t)}+r(t)*acf{qc(t)}
[0086] 其中，qc(-t)为低次谐波，ccf{q(t)，qc(t)}为震源基波信号和次谐波信号的互相关，acf{qc(t)}为次谐波信号的自相关。
[0087] 所述切除噪声单元开取的时窗的宽度为次谐波噪声的一个波长的宽度。
[0088] 为了验证这种压制次谐波干扰方法可行性，对模拟资料进行测试。设计震源
频率范围为5～70hz，采样率为2ms，扫描时间10s的线性升频基波扫描信号。同时设计频率范围为2.5～35hz，采样率为2ms，扫描时间10s的次谐波信号。基波扫描信号和次谐波信号同时作为震源信号，与反射系数褶积，得到相关前合成记录如图1所示。现要对合成记录中的谐波干扰采用本发明的次谐波压制方法进行处理。
[0089] 首先，将次谐波信号与含次谐波干扰的相关前记录进行相关处理，将记录变换到次谐波域，得到相关后记录，如图2所示。
[0090] 再利用面波视速度对相关后记录进行线性时差。在校正后记录中，次谐波干扰能量集中在一个较窄时窗内(如图3中的细长方形框所示)，开取时窗，对这部分干扰波能量进行切除，如图4所示。
[0091] 最后利用面波视速度对切除干扰波能量后的记录进行线性反校正(如图5所示)，并用次谐波信号将记录从“次谐波域”反变换到常规空间-时间域，如图6所示。
[0092] 为了对比次谐波压制前后效果，用基波扫描信号分别与次谐波压制前、后的记录进行互相关处理，如图7a和图7b所示。可以看出次谐波噪声得到了很好压制。
[0093] 本发明方法将噪声波列变换到“次谐波域”，使得谐波噪声成分集中起来，并与基波信号分离开，使得谐波信号更容易被检测到并得到压制。
[0094] 上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的
基础上，很容易做出各种类型的改进或
变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。
