一种自动计算地震数据光缆耦合噪声强弱程度的振幅谱方法与流程

本发明属于地震勘探数据处理领域，特别涉及一种光缆耦合噪声强弱程度的自动计算方法。

背景技术：

分布式声学传感(distributedacousticsensing,das)是在一种出现后得到快速发展的技术，它的工作原理是基于散射效应，并且集合了传感和传输两种功能，能够实现远距离的检测和测量。das系统的光纤仪器抗电磁干扰、耐高温高压，可以不挪动井下仪器而重复试验，从而增加产量。das系统记录得到的vsp(verticalseismicprofile)数据，在成像、地层检测等方面有很广泛的应用，信噪比能够满足数据成像的要求。

但是，由于光纤缆绳在采集与传输过程中没有贴靠到井壁上并受到振动事件的影响，导致产生强烈的振荡，而且振荡波的接触点在整个观测过程中相当稳定，因此采集的vsp记录受到光缆耦合噪声干扰，降低了地震数据的信噪比，影响数据的后续分析处理。而计算光缆耦合噪声强弱程度，可以更好地压制光缆耦合噪声，增强数据处理的效果。

现有技术：

经验法。该方法主要是从业者根据地震数据处理经验，通过观察和比对，来判断确定地震数据光缆耦合噪声的强弱程度。

现有技术的缺点：

1、该方法判断确定的指标，属于一个主观估计值，不够精确，在精确性方面有所欠缺。

2、该方法针对多道数据确定一个综合强弱程度的指标，无法针对每一道数据确定强弱程度。

技术实现要素：

本发明的目的在于提供一种自动计算地震数据光缆耦合噪声强弱程度的振幅谱方法。本发明通过将单道时域数据转换到频率域，计算频率域离散采样值的方差或峰度，并作为判断衡量光缆耦合噪声强弱程度的指标。这种非人工的自动判断方法提高了判断衡量单道数据光缆耦合噪声强弱程度的精确性，是自适应压制光缆耦合噪声的基础。

为了实现上述目的，本发明采用如下技术方案：

一种自动计算地震数据光缆耦合噪声强弱程度的振幅谱方法，包括以下步骤：

步骤01：读取das(distributedacousticsensing)采集的vsp(verticalseismicprofile)记录的单道数据；

步骤02：对步骤01读取的单道数据做傅里叶变换，从时域转换到频率域；

步骤03：将频率域离散采样值x[k]归一化处理；

步骤04：计算归一化振幅谱的方差或峰度。

进一步的，重复步骤01-04直到vsp记录的所有道数据处理完成。

进一步的，步骤02中对地震资料单道数据做傅里叶变换，从时域转换到频率域，具体包括：

离散的单道数据采样点值记为x[n]，x[k]为x[n]的离散傅里叶变换，将单道地震数据从时域变换到频率域，使用快速傅里叶变来求取离散的频率值：

x[k]＝dft(x[n]),

假设单道地震数据x[n]的采样点个数为n，则通过fft算法得到的频率域离散采样点个数也为n。由于fft算法得到的频谱以nyquist频率对称，因此考虑前n/2个频率域采样值，也就是0-nyquist频率范围内的频谱。由于实际地震信号的频带有限，因此本发明只考虑0-一半的nyquist频率之间的频谱。假设ωk为频谱的第k点的离散频率，则有下式：

其中，dt为采样间隔，则为采样频率，采样频率记为fn，采样频率一半记为fn/2。

进一步的，步骤03中将频率域离散采样值x[k]归一化处理，包括：

归一化后的频率域离散采样值记为y[k]，则有：

y[k]＝x[k]/m,

式中，m为振幅谱离散采样值绝对值的最大值，即：

m＝max(abs(x[k])。

进一步的，步骤04中计算归一化振幅谱的方差或振幅谱的峰度，包括：

首先计算归一化振幅谱方差，记为v，即计算0-一半的nyquist频率之间的离散采样值方差：

其中，为y[k]的均值：

计算归一化振幅谱峰度，记为p，即计算0-一半的nyquist频率之间的离散采样值的峰度：

进一步的，振幅谱方差越小，光缆耦合噪声越强；振幅谱峰度越大，光缆耦合噪声越强。

本发明称根据振幅谱方差或者峰度来计算das系统采集地震记录中光缆耦合噪声强弱程度的方法为“振幅谱方法”。

相对于现有技术，本发明具有以下有益效果：利用本发明的自动计算光缆耦合噪声强弱程度的振幅谱方差或振幅谱峰度方法，得到的方差或峰度指标都能够很好地衡量光缆耦合噪声的强弱程度。同时在振幅谱方差或振幅谱峰度的基础上，可以确定稀疏优化压制光缆耦合噪声的算法参数，取得较好的压噪结果。

附图说明

图1为das系统采集vsp记录数据；

图2a为图1所示数据的第350道数据振幅谱；图2b为图1所示数据的第350道数据时频谱；

图3a为图1所示数据的第400道数据振幅谱；图3b为图1所示数据的第400道数据时频谱；

图4为实际地震数据；

图5a为图4所示数据抽取第750道数据的振幅谱；图5b为图4所示数据抽取第1110道数据的振幅谱；图5c为图4所示数据抽取第400道数据的振幅谱；

图6为地层模型结构；

图7a为有效信号模型；图7b为光缆耦合噪声模型；图7c为有效信号与五倍噪声合成模型数据；图7d为有效信号与十倍噪声合成模型数据；

图8为图7c与图7d所示合成地震模型的振幅谱方差；

图9为图7c与图7d所示合成地震模型的振幅谱峰度；

图10为实际地震数据；

图11为图10实际地震数据振幅谱方差；

图12为图10实际地震数据振幅谱峰度；

图13a为未受光缆耦合噪声干扰模拟数据；图13b为受光缆耦合噪声干扰模拟数据；图13c为图13b所示数据分离的有效信号；图13d为图13b所示数据分离的光缆耦合噪声；

图14a为实际地震数据；图14b为图14a所示实际地震数据分离的有效信号；图14c为图14a所示实际地震数据分离的光缆耦合噪声；

图15为本发明的流程图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步详细的说明。

本发明为一种自动计算地震数据光缆耦合噪声强弱程度的振幅谱方法，通过将单道地震数据从时间域转换到频率域，并计算归一化振幅谱的方差或峰度，并将此值作为衡量光缆耦合噪声强弱程度的指标。

请参阅图15所示，本发明为一种自动计算地震数据光缆耦合噪声强弱程度的振幅谱方法，对地震记录信号实施步骤01-步骤04，具体包括：

步骤01：读取das采集的vsp记录的单道数据；

步骤02：对步骤01读取的单道数据做傅里叶变换，从时域转换到频率域；

步骤03：将频率域离散采样值x[k]归一化处理；

步骤04：计算归一化振幅谱的方差或峰度；

重复步骤01-04直到vsp记录的所有道数据处理完成。

如图1所示为das系统采集记录的实际地震数据，共2000道，采样点数为2000，采样间隔为1ms。可以看到，该实际数据受到了强烈的光缆耦合噪声干扰。光缆耦合噪声降低了地震数据的信噪比，严重时掩盖了有效信号的同相轴，影响后续的波场分离、数据成像、属性分析等应用，因此能否有效压制该噪声是das技术在地震勘探领域下一步发展的关键因素。而能否有效计算光缆耦合噪声强弱程度，是自适应压制光缆耦合噪声的一个重要内容。

抽取图1所示实际地震记录中光缆耦合噪声比较明显的第350道与几乎不含光缆耦合噪声的第400道数据，展示其振幅谱与时频谱，结果如图2a、图2b与图3a、图3b所示。可以看到，有效信号的时频谱表现为一定的脉冲特征，而光缆耦合噪声的振幅谱上表现为若干个近似脉冲特征，时频谱上表现为若干条水平直线的特征，也就是说，光缆耦合噪声主要表现为若干个单一频率信号的叠加。

进一步的，步骤02中对地震资料单道数据做傅里叶变换，从时域转换到频率域，具体包括：

离散的单道数据采样点值记为x[n]，x[n]的离散傅里叶变换记为x[k]，使用fft来求取离散的频率值：

x[k]＝dft(x[n]).

单道地震数据x[n]的采样点个数为n，通过fft算法得到的频率域离散采样点x[k]个数也为n。只需要考虑前n/2个频率域采样值，也就是0-一半的nyquist频率范围内的频谱。假设ωk为频谱的第k点的离散频率，则有下式：

其中，dt为采样间隔，则为采样频率，采样频率记为fn，采样频率一半记为fn/2。

进一步的，步骤03中将频率域离散采样值x[k]归一化处理，具体包括：

归一化后的频率域离散值记为y[k]，则有：

y[k]＝x[k]/m,

式中，m为步骤02得到的振幅谱离散采样值绝对值的最大值，即：

m＝max(abs(x[k]).

如图4所示，为das实际地震资料，从该地震资料中选取三道数据，即第750道、1110道、400道，很明显这三道数据分别含有较强、中等、较弱光缆耦合噪声。如图5a、图5b、图5c所示为这三道数据的振幅谱图。可以发现，光缆耦合噪声越强，振幅谱越表现为若干个单峰特征，光缆耦合噪声越弱，振幅谱的这种单峰特征越不明显。

进一步的，步骤04中计算归一化振幅谱的方差或峰度，包括：

首先计算归一化振幅谱方差，记为v,即计算0-一半的nyquist频率之间的频率域离散采样点的方差，则有：

其中，为y[k]的均值：

接下来计算归一化振幅谱峰度，也称作四阶累积量，即计算0-一半的nyquist频率之间的频率域离散采样点的峰度，记为p，则有：

需要说明的是，计算振幅谱峰度可以不对振幅谱做归一化，与归一化后的振幅谱峰度结果是相同的；但是，计算振幅谱方差之前必须对振幅谱归一化，否则将不能有效衡量光缆耦合噪声强弱程度。振幅谱方差越小，光缆耦合噪声越强；振幅谱峰度越大，光缆耦合噪声越强。

使用图6所示的地层模型可以得到图7a所示的合成零偏移距vsp记录，包括上行波场与下行波场。针对das光缆耦合噪声的地震记录特征，给出了光缆耦合噪声的合成记录，如图7b所示。

将如图7a所示的有效信号与如图7b所示的光缆耦合噪声模型进行叠加，包括将一倍的有效信号加上五倍的噪声得到如图7c所示的含噪数据，将一倍的有效信号加上十倍的噪声得到如图7d所示的含噪数据，很明显，图7d所示数据的光缆耦合噪声要强于图7c。

分别计算图7c与图7d所示数据的振幅谱方差，结果如图8所示。可以看到：第一，无论对于含有强噪声还是弱噪声的数据，在21-60道的区域振幅谱方差都比其它道区域小，这是因为模型数据只在21-60道区域受到了光缆耦合噪声干扰；第二，在21-60道受到噪声干扰的区域，强噪声数据各道数据的振幅谱方差比弱噪声的振幅谱方差要小。因此，可以验证，光缆耦合噪声越强，振幅谱方差越小，振幅谱方差是自动计算光缆耦合噪声强弱程度的有效方法。

分别计算图7c与图7d所示数据的振幅谱峰度，结果如图9所示。可以看到：第一，无论对于含强噪数据还是含弱噪数据，在21-60道的区域峰度都比其它区域大，同样是因为模型数据只在21-60道区域受到光缆耦合噪声干扰；第二，在21-60道受到噪声干扰的区域，含强噪声模型的振幅谱峰度比含弱噪声模型的振幅谱峰度要大。因此，模型数据验证表明，光缆耦合噪声越强，振幅谱峰度越大，振幅谱峰度可以自动计算衡量光缆耦合噪声强弱程度。

图10是das系统采集的实际vsp记录，计算该实际地震数据的振幅谱方差，结果如图11所示。从图11可以看到振幅谱方差可以很精确地反映实际地震资料各道数据中光缆耦合噪声的强弱程度，即噪声越强，振幅谱方差越小。接下来，抽取图10所示实际地震资料的第120道、220道、280道数据，这三道数据受到光缆耦合噪声干扰越来越强，计算其振幅谱方差，结果如表1所示，可以发现，随着光缆耦合噪声强度的增强，振幅谱方差逐渐减小。因此，振幅谱方差是一种可以自动计算光缆耦合噪声强弱程度的有效方法。

表1抽取单道数据的振幅谱方差

计算图10所示实际地震数据所有道的振幅谱峰度，结果如图12所示。从图12可以发现，振幅谱峰度可以很精确地反映实际地震资料各道数据中光缆耦合噪声的强弱程度，即噪声越强，振幅谱峰度越大。接下来，我们同样抽取图10所示实际地震资料的第120道、220道、280道数据，计算振幅谱峰度，计算结果如表2所示，可以发现，随着光缆耦合噪声强度的增强，振幅谱峰度逐渐增大。因此，振幅谱峰度是一种可以自动计算光缆耦合噪声强弱程度的有效方法。

表2抽取单道数据的振幅谱峰度

本发明具有如下有益效果：

1)本发明方法对单道数据计算振幅谱方差或者峰度，作为衡量光缆耦合噪声强弱程度的指标，具有很强的自适应性，可以作为压制光缆耦合噪声的有效参考；

2)本发明方法使用快速傅里叶变换将地震信号从时域转换到频率域，计算效率很高；

3)本发明方法在实现时，各道数据之间没有联系，可以并行处理，适合工业数据处理。

下面将基于本发明的稀疏优化压噪方法应用到das采集vsp记录中有效信号和光缆耦合噪声的分离。应用结果表明，基于光缆耦合噪声强弱程度的判断，光缆耦合噪声压制效果比较理想。

比较图13a所示的vsp模型数据与图13b所示的受到光缆耦合噪声干扰的vsp数据，可以发现，有效信号受到很强的光缆耦合噪声干扰，部分有效信号信息甚至被掩盖。使用振幅谱峰度来计算光缆耦合噪声强弱程度，并根据噪声强弱程度选取稀疏优化算法的参数，得到如图13c所示的有效信号，和如图13d所示的光缆耦合噪声。对比图13a与图13c可知，有效地压制了光缆耦合噪声。

利用实际das地震数据(图14a所示)进行实验验证。该炮集数据总共2000道，采样时间间隔为2ms，记录长度为6.0s。从剖面上可以看出，该记录中的信号受到很强的光缆耦合噪声干扰，导致有效信号被噪声覆盖，信噪比降低，对地震资料的分析以及解释造成严重影响。使用振幅谱峰度来计算光缆耦合噪声强弱程度，并根据噪声强弱程度选取稀疏优化算法的参数，获得有效信号和光缆耦合噪声，分别如图14b和图14c所示。比较图14a、图14b发现，可以有效地压制地震记录中的光缆耦合噪声，并且被光缆耦合噪声掩盖的有效信号同相轴在有效信号剖面中清晰地显示。

以上的算例实验中，利用本发明的振幅谱方法可以自动计算das采集vsp记录中的光缆耦合噪声强弱程度，该指标为后续资料的分析奠定基础，同时本发明方法具有很强的自适应性。

最后需要说明的是，以上模型和实际资料算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围，在本发明的精神和原则之内，所做的任何修改，改进或等同替换等，均应在本发明的保护范围内。