CN103926623A - 一种压制逆时偏移低频噪音的方法 - Google Patents

Abstract

本方法公开一种压制逆时偏移低频噪音的方法，该方法对基于声波方程和震源归一化互相关成像条件得到的Walkaway VSP资料逆时偏移成果数据，应用拉普拉斯去噪后的剩余的噪音能量进行压制，噪音压制采用的方法是改进的非局部均值滤波法，其相似系数的计算主要依据所输入的用于逆时偏移的深度域层速度模型。该方法应用于经拉普拉斯滤波后的成像体处理中，可以有效的压制低频噪音，从而使地下的地质体成像更加清晰、准确。

Description

一种压制逆时偏移低频噪音的方法

技术领域

本发明属于地球物理勘探技术领域，涉及一种压制逆时偏移低频噪音的方法，其具体是Walkaway VSP地震资料逆时偏移处理中的去噪方法，尤其适用于对Walkaway VSP逆时偏移方法得到的最后成像数据体进行去噪处理，从而使成像结果更加清晰。

背景技术

逆时偏移技术是目前地震勘探领域内成像精度最高的地震偏移方法。在勘探目标日趋复杂多样的今天，逆时偏移技术在复杂构造成像方面的优势将有助于地震勘探向深层隐蔽性油气藏探索。逆时偏移技术最早出现在1983年，当时Whitmore(1983)、Baysal(1983)、Loewenthal(1983)、Levin(1984)、Hildebrand(1987)等地球物理学家分别对逆时偏移的概念、基本原理、应用范围等展开研究，为逆时偏移的发展奠定了基础。在逆时偏移过程中，成像条件的选择对最后的成像效果有着至关重要的影响，Hu和McMechan(1986)与Chang和MeMchan(1986)将激发时刻成像条件应用于叠前逆时偏移处理中；Hu(1986)基于直达波振幅能量最强利用有限差分法进行震源波场沿拓计算激发时刻成像条件。Whitmore和Lines(1986)利用震源波场和检波点波场进行互相关。Kaelin和Guitton(2006)在互相关的基础上对震源归一化，显著提高了成像数据体反射信息的精度。因互相关成像条件会给成像体带来严重的低频噪音，Yu Zhang等(2009)提出拉普拉斯算子滤波，相当于对成像波场进行低角度域衰减，对低频噪音的去除取得了明显的效果。

非局部均值滤波(NLM)最早由Buades(2005)提出，算法利用成像点窗口与临近创建间的相似性来加强构造信息，从而有效地衰减随机噪音。由于NLM算法对每一个成像点降噪时需要计算全部成像点的相似系数，因此具有较大的计算量。许多学者为减少计算时间做了不同的研究，Sheng(2009)将算法在GPU上运行，取得了较好的效果。Mahmoudi和Sapiro(2005)，Brox等(2008)为减少每个成像点的计算时间，将求取相似系数限制在成像点为中心的一个范围内。这些算法较原本的算法节省了十几倍、甚至上百倍的计算时间，使得NLM算法成功的应用于医学数据(Coupe et al.,2008)、雷达数据(Deledalle et al.,2011)、音频数据(Zoican,2010)、显微成像(Wei and Yin,2010)等领域。

现有技术中，虽说逆时偏移技术已经应用于地球物理勘探技术领域，但是在地震噪声的处理上一直没有取得进一步的突破，尤其在针对Walkaway VSP逆时偏移中噪音压制处理环节。

发明内容

针对上述技术难题，本发明人经过多年研究，针对Walkaway VSP逆时偏移中噪音压制处理环节，提供了一种改进的滤波方法对成像体进行去噪处理，也就是提供一种压制逆时偏移低频噪音的方法。

依据本发明的技术方案，压制逆时偏移低频噪音的方法包括以下步骤：

1)在井中垂向等间隔布置检波器，地表设置过井口的炮线，炮点等间隔分布，人工激发地震波，将井中接收到的Walkaway VSP地震信号记录到磁带上；

2)从磁带读取地震记录，做常规的噪音压制、反褶积、波场分离和速度分析等预处理，得到预处理后的炮记录和速度模型；

3)读取Walkaway VSP炮记录，进行傅立叶变换，得到炮记录的振幅谱，分析得到炮记录的主频，记为f_p，单位为Hz，按照如下公式计算地震子波：

$f\left(t\right)=[1-2{\left(\mathrm{\π}{f}_{p}t\right)}^2]e^{-{\left(\mathrm{\π}{f}_{p}t\right)}^2}---\left(1\right)$

式中，t为时间，单位是s；

4)读取速度模型，将步骤3)得到的子波作为震源子波，应用双程波方程计算正传波场；

5)读取Walkaway VSP炮记录和速度模型，将炮记录作为输入，应用双程波方程计算反传波场；

6)将步骤4)和步骤5)从第一炮开始循环，逐炮做数据处理，对正传和反传波场采用震源归一化互相关成像条件进行成像；

7)对成像数据体进行拉普拉斯滤波，得到滤波后的地震数据；

8)读取速度模型，求取非局部均值滤波计算过程需要的相似系数；

9)对滤波后的地震数据进行非局部均值滤波处理，得到最终的输出成果数据，成果数据将指示地质层位划分和确定微小断裂体系的解释方案，同时用于储层预测和含石油天然气有利区的识别。

优选地，步骤3)所采用的Walkaway VSP逆时偏移需要的震源子波设计为雷克子波。

优选地，步骤4)和步骤5)的波场正、反传所采用的声波方程可写为：

$\frac{{\∂}^{2}w}{\∂t^2}=v^2(\frac{{\∂}^{2}w}{\∂x^2}+\frac{{\∂}^{2}w}{\∂z^2})$

式中：w为波场值，是空间坐标(x，y，z)和时间t的函数，v是相应空间位置的速度，单位是m/s。

步骤4)读取速度模型v后，令x和z方向的空间网格步长Δx＝Δz＝Δd，用下式确定时间采样步长。

$\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;d}}<1/\sqrt{2\underset{m=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{b}_{2m-1}}$

式中，N为差分阶数，N₁为小于N的最大奇数；Δt为时间步长，单位是s，b为差分系数。

步骤6)中对每一时刻的震源波场、检波点波场互相关，累加所有时刻的相关结果，并与所有时刻震源波场自相关的累加值做比。震源归一化互相关成像条件的计算公式：

$M_{{\mathrm{RS}}_{i,j}}=\frac{{{\&Integral;W}_S}_{i,j}^t*{W_R}_{i,j}^t\mathrm{dt}}{{{\&Integral;W}_S}_{i,j}^t*{W_S}_{i,j}^t\mathrm{dt}}$

式中，为震源波场，为检波点波场，i，j分别为x，z方向的坐标，单位是m，t为时间，单位是s。是先对每一时刻的震源波场与检波点波场互相关，再对震源波场自相关，最后将累加所有时刻的互相关结果与累加所有时刻的自相关结果做比。

步骤7)的拉普拉斯滤波对步骤6)的成像数据体进行相当于角度域的滤波处理。

步骤8)的非局部均值滤波先对速度模型求取相似系数，对经拉普拉斯滤波后的数据体进行非局部均值滤波，以压制地震偏移数据中的剩余噪音。

相对于现有技术，本发明的压制逆时偏移低频噪音的方法具有下述优点：

1)本方法逆时偏移后的成像数据体几乎不含任何假设，具有较强的适应性；

2)本发明实施简单，可操作性强，可根据实际情况选择不同的参数；

3)本方法可明显改善拉普拉斯滤波后成像数据体的成像效果。

附图说明

附图1为成像点间相似系数的计算原理图。

附图2为含低频噪音的成像地震剖面。

附图3为经拉普拉斯滤波后的地震剖面。

附图4为速度模型。

附图5为非局部均值滤波后的成像地震剖面。

附图6为图3的局部放大图。

附图7为图5的局部放大图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。另外地，不应当将本发明的保护范围仅仅限制至下述具体步骤或具体参数。

本发明的压制逆时偏移低频噪音的方法基于以下研究：

(1)震源归一化互相关成像条件是在原本互相关成像条件下进行震源归一化，互相关成像条件可表示为如下形式：

$M_{{\mathrm{RS}}_{i,j}}={\&Integral;W_S}_{i,j}^t*{W_R}_{i,j}^t\mathrm{dt}$

式中，为震源波场，为检波点波场。其中，i，j分别为x，z方向的坐标，单位是m，t为波场的不同时刻，单位是s，成像值为所有时刻相关结果的累加。

(2)互相关成像条件本身的保幅效果较差，震源归一化成像条件与之相比更具有振幅保真性，其形式如下：

$M_{{\mathrm{RS}}_{i,j}}=\frac{{{\&Integral;W}_S}_{i,j}^t*{W_R}_{i,j}^t\mathrm{dt}}{{{\&Integral;W}_S}_{i,j}^t*{W_S}_{i,j}^t\mathrm{dt}}$

(3)与互相关成像条件相比，震源归一化互相关成像条件可以在一定程度内衰减浅层能量，提高深层能量。

震源归一化互相关成像条件在成像过程中会不可避免的引入低频噪音，这给成像结果带来了一定程度的干扰。拉普拉斯滤波是普遍采用的用来去除低频噪音的滤波方法，它直接对成像体进行处理，对低频噪音的去除效果明显。

拉普拉斯算子可以表示为：

${\&dtri;}^2={\&PartialD;}_x^2+{\&PartialD;}_z^2$

式中，分别表示对x，z的二阶偏导数。

其在波数域可表示为：

$\mathrm{FFT}\left({\&dtri;}^2\right)=-(k_x^2+k_z^2)=-{\left|\stackrel{\&RightArrow;}{k_I}\right|}^2$

式中，是波数矢量，有FFT代表傅立叶变换。

由余弦定理可知：

${\left|\stackrel{\&RightArrow;}{k_I}\right|}^2={\left|\stackrel{\&RightArrow;}{k_R}\right|}^2+{\left|\stackrel{\&RightArrow;}{k_S}\right|}^2-2\left|\stackrel{\&RightArrow;}{k_R}\right|\left|\stackrel{\&RightArrow;}{k_S}\right|\cos (\mathrm{\&pi;}-2\mathrm{\&theta;})=\frac{{4\mathrm{\&omega;}}^2}{v^2}{\cos }^2\mathrm{\&theta;}$

式中，θ为入射角，分别为检波点波场与震源波场的波数矢量。

成像数据体经拉普拉斯滤波后，仍存在一定的噪音残余，非局部均值滤波可以很好的滤除掉噪音。非局部均值滤波的原理如下：

一个含有噪音的成像结果I可以表示为如下形式：

I＝U+N

式中，U为无噪音的成像结果，N为噪音。

对于成像结果中的任意位置m，去噪后的成像结果可以表示为：

$\hat{I}\left(m\right)=\underset{k}{\mathrm{\&Sigma;}}w(m,k)I\left(k\right)$

式中，w(m，k)为m和k之间成像值的相似性，其满足0≤w(m，k)≤1且Σ_kw(m，k)＝1。其中，每个m都存在与k相关的相似系数。为了定量表示相似性，定义N_m为以m为中心的窗口。为计算N_m和N_k之间的相似性，采用高斯加权欧几里得距离表示：

$\begin{array}{c}{D}^2(m,k)={\left|\right|I\left(N_m\right)-I\left(N_k\right)\left|\right|}_{2,a}^2\\ =\underset{l}{\overset{\mathrm{nl}}{\mathrm{\&Sigma;}}}{\left[G_u\left(l\right)(I\left(N_m\left(l\right)\right)-I\left(N_k\left(l\right)\right))\right]}^2\end{array}$

式中，为高斯加权欧几里得距离的平方，G_u为高斯函数，其形式为：

$G_u(x,y)=\exp (-\frac{{(x-x_0)}^2+{(y-y_0)}^2}{2u})$

式中，x₀，y₀为高斯函数中心，其中x，y与上式中的l项对应，u为偏离参数。

N_m和N_k之间的相似系数可以表示为：

$w(m,k)=\frac{1}{Z\left(m\right)}\exp \left(\frac{-D^2(m,k)}{h^2}\right)---(*)$

其中，为规则化系数，以确保Σ_kw(m，k)＝1。

下面对发明进一步阐述，本发明不局限于Walkaway VSP逆时偏移，对地面地震逆时偏移同样适用，并且可对实际资料进行处理，具有广泛的适应性。

本发明的压制逆时偏移低频噪音的方法包括以下步骤：

1)在井中垂向等间隔布置检波器，地表设置过井口的炮线，炮点等间隔分布，人工激发地震波，将井中接收到的Walkaway VSP地震信号记录到磁带上；

2)从磁带读取地震记录，做常规的噪音压制、反褶积、波场分离和速度分析等预处理，得到预处理后的炮记录和速度模型；

3)读取Walkaway VSP炮记录，进行傅立叶变换，得到炮记录的振幅谱，分析得到炮记录的主频，记为f_p，单位为Hz，按照如下公式计算地震子波：

$f\left(t\right)=[1-2{\left(\mathrm{\&pi;}{f}_{p}t\right)}^2]e^{-{\left(\mathrm{\&pi;}{f}_{p}t\right)}^2}---\left(1\right)$

4)读取速度模型，将步骤3)得到的子波作为震源子波，应用双程波方程计算正传波场；

5)读取Walkaway VSP炮记录和速度模型，将炮记录作为输入，应用双程波方程计算反传波场；

6)将步骤4)和步骤5)从第一炮开始循环，逐炮做数据处理，对正传和反传波场采用震源归一化互相关成像条件进行成像；

7)对成像数据体进行拉普拉斯滤波，得到滤波后的地震数据；

8)读取速度模型，求取非局部均值滤波计算过程需要的相似系数；

9)对滤波后的地震数据进行非局部均值滤波处理，得到最终的输出成果数据，成果数据将指示地质层位划分和确定微小断裂体系的解释方案，同时用于储层预测和含石油天然气有利区的识别。

优选地，步骤3)所采用的Walkaway VSP逆时偏移需要的震源子波设计为雷克子波。

优选地，步骤4)和步骤5)的波场正、反传所采用的声波方程可写为：

$\frac{{\&PartialD;}^{2}w}{\&PartialD;t^2}=v^2(\frac{{\&PartialD;}^{2}w}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}w}{\&PartialD;z^2})$

式中：w为波场值，是空间坐标(x，y，z)和时间t的函数，v是相应空间位置的速度，单位是m/s。

步骤4)读取速度模型v后，令x和z方向的空间网格步长Δx＝Δz＝Δd，用下式确定时间采样步长。

$\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;d}}<1/\sqrt{2\underset{m=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{b}_{2m-1}}$

式中，N为差分阶数，N₁为小于N的最大奇数；Δt为时间步长，b为差分系数。

步骤6)中对每一时刻的震源波场、检波点波场互相关，累加所有时刻的相关结果，并与所有时刻震源波场自相关的累加值做比。震源归一化互相关成像条件的计算公式：

$M_{{\mathrm{RS}}_{i,j}}=\frac{{{\&Integral;W}_S}_{i,j}^t*{W_R}_{i,j}^t\mathrm{dt}}{{{\&Integral;W}_S}_{i,j}^t*{W_S}_{i,j}^t\mathrm{dt}}$

式中，为震源波场，为检波点波场，i，j分别为x，z方向的坐标，单位是m，t为时间，单位是s。是先对每一时刻的震源波场与检波点波场互相关，再对震源波场自相关，最后将累加所有时刻的互相关结果与累加所有时刻的自相关结果做比。

步骤7)的拉普拉斯滤波对步骤6)的成像数据体进行相当于角度域的滤波处理。

步骤8)的非局部均值滤波先对速度模型求取相似系数，对经拉普拉斯滤波后的数据体进行非局部均值滤波，以压制地震偏移数据中的剩余噪音。

下面将结合附图及实例对发明进一步示意阐述。

参考附图1，图中显示的是以1框中圆点为中心的窗口与以周围点为中心的窗口之间计算相似系数，如计算1框中心点与2框中心点之间的相似系数，首先确定窗口大小，即以点为中心的所包括的数据大小，然后对1框与2框之间的数据计算高斯欧几里得距离，最后得到两框中心点之间的相似系数。非局部均值滤波原本是对全局计算相似系数，如计算1框中心点与成像体内其余所有点之间的相似系数，但这会产生较大的计算量，因此仅在以1框中心点为中心的虚线边框内的所有点来计算相似系数。

参考附图2为Walkaway VSP逆时偏移互相关成像后的结果，可以看出具有明显的低频噪音，影响最后的成像效果。

参考附图3为经拉普拉斯滤波后的成像体，可以看出低频噪音得到较好的压制，说明拉普拉斯滤波具有较为明显的去噪效果。

参考附图4是成像体的速度模型，为实际验证非局部均值滤波的效果，将实际速度模型进行平滑，作为非局部均值滤波求取相似系数的输入图形。

参考附图5是经非局部均值滤波后的成像体，可以看出噪音明显减少，地震成像更清晰。

参考附图6和图7分别是图3和图5的局部放大对比图，对比白色箭头所示位置可以看出，经非局部均值滤波后，成像数据的噪音得到进一步的压制，局部构造更清晰。

显然上述实施例仅为清楚的说明本发明所做的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上，还可以做出其他不同形式的变化或变动，这里无需也无法对所有实施方式予以穷举。由此所引申的显而易见的变化或变动仍处于本发明创造的保护范围之中。

Claims (7)

1.一种压制逆时偏移低频噪音的方法，其包括以下步骤：

1)在井中垂向等间隔布置检波器，地表设置过井口的炮线，炮点等间隔分布，人工激发地震波，将井中接收到的Walkaway VSP地震信号记录到磁带上；

2)从磁带读取地震信号记录，做常规的噪音压制、反褶积、波场分离和速度分析的预处理，得到预处理后的Walkaway VSP炮记录和速度模型；

3)读取Walkaway VSP炮记录，进行傅立叶变换，得到炮记录的振幅谱，分析得到炮记录的主频，记为f_p，单位为Hz，按照如下公式(1)计算地震子波：

$f\left(t\right)=[1-2{\left(\mathrm{\&pi;}{f}_{p}t\right)}^2]e^{-{\left(\mathrm{\&pi;}{f}_{p}t\right)}^2}---\left(1\right)$

式中，t为时间，单位是s；

4)读取速度模型，将步骤3)得到的地震子波作为震源子波，应用双程波方程计算正传波场；

5)读取Walkaway VSP炮记录和速度模型，将炮记录作为输入，应用双程波方程计算反传波场；

6)将步骤4)和步骤5)从第一炮开始循环，逐炮做数据处理，对正传波场和反传波场采用震源归一化互相关成像条件进行成像，得到成像数据体；

7)对成像数据体进行拉普拉斯滤波，得到滤波后的地震数据；

8)读取速度模型，求取非局部均值滤波计算过程需要的相似系数；

9)对滤波后的地震数据进行非局部均值滤波处理，得到最终的输出成果数据；

10)将得到成果数据用于指示地质层位划分和确定微小断裂体系的解释方案，同时用于储层预测和含石油天然气有利区的识别。

2.根据权利要求1所述的压制逆时偏移低频噪音的方法，其特征在于，所述步骤3)所采用的Walkaway VSP逆时偏移需要的震源子波设计为雷克子波。

3.根据权利要求1所述的压制逆时偏移低频噪音的方法，其特征在于，步骤4)中的正传波场和步骤5)中的反传波场所采用的声波方程如下：

$\frac{{\&PartialD;}^{2}w}{\&PartialD;t^2}=v^2(\frac{{\&PartialD;}^{2}w}{\&PartialD;x^2}+\frac{{\&PartialD;}^{2}w}{\&PartialD;z^2})$

式中：w为波场值，是空间坐标(x，y，z)和时间t的函数，v是相应空间位置的速度，单位是m/s。

4.根据权利要求1所述的压制逆时偏移低频噪音的方法，其特征在于，所述步骤4)从炮集数据中读取震源位置坐标、检波点位置坐标，读取速度模型v，选取正方形空间计算网格，令x和z方向的空间网格步长Δx＝Δz＝Δd，用下式(2)确定时间采样步长：

$\frac{\mathrm{v\&Delta;t}}{\mathrm{\&Delta;d}}<1/\sqrt{2\underset{m=1}{\overset{N_1}{\mathrm{\&Sigma;}}}{b}_{2m-1}}---\left(2\right)$

上式(2)中，N为差分阶数，N₁为小于N的最大奇数；Δt为时间步长，b为差分系数。

5.根据权利要求1所述的压制逆时偏移低频噪音的方法，其特征在于，所述步骤6)中对每一时刻的震源波场、检波点波场互相关，累加所有时刻的相关结果，并与所有时刻震源波场自相关的累加值做比；震源归一化互相关成像条件的计算公式如下：

$M_{{\mathrm{RS}}_{i,j}}=\frac{{{\&Integral;W}_S}_{i,j}^t*{W_R}_{i,j}^t\mathrm{dt}}{{{\&Integral;W}_S}_{i,j}^t*{W_S}_{i,j}^t\mathrm{dt}}$

式中，为震源波场，为检波点波场，i，j分别为x，z方向的坐标，单位是m(米)，t为时间，单位是s(秒)。

6.根据权利要求1所述的压制逆时偏移低频噪音的方法，其特征在于，所述步骤7)的拉普拉斯滤波对步骤6)的成像数据体进行相当于角度域的滤波处理。

7.根据权利要求1所述的压制逆时偏移低频噪音的方法，其特征在于，所述步骤8)的非局部均值滤波先对速度模型求取相似系数，再对步骤7)经拉普拉斯滤波处理后的数据进行滤波处理。