CN104714251B - 用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法 - Google Patents

Abstract

本发明用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，包括如下步骤，1.将叠前CSP数据用小波变换得到信号的解析部分；2.根据步骤1得到的解析部分得到信号的瞬时振幅；3.在瞬时振幅的剖面上进行局部线性Radon变换，拾取瞬时振幅的峰值振幅，得到倾斜叠加峰值振幅剖面；4.基于倾斜叠加峰值振幅剖面得到边缘检测剖面，用于在边缘检测剖面上拾取同相轴及走时。在瞬时振幅剖面上进行倾斜叠加变换，因为瞬时振幅剖面相较原始地震记录更加圆滑，包含更多的低频分量，且可滤除部分高频毛刺噪声，而二者的同相轴位置相同。有效压制了噪声，克服了现有方法难以在大偏移距内拾取连续同相轴的困难，且操作方便，实现了同相轴自动拾取。

Description

用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法

技术领域

本发明涉及地震勘探技术领域，具体为一种用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法。

背景技术

走时是衰减和速度层析成像的重要参数之一，我们需要拾取每个同相轴不同偏移距处反射波的走时，包括大偏移距处的同相轴。可以看出，走时拾取的前提是同相轴拾取，同相轴拾取也是地震解释的重要步骤之一。但由于噪声的影响，大偏移距处道集的信噪比低，拾取连续同相轴比较困难。为了拾取正确的反射波走时，我们需要在大偏移距范围内拾取连续同相轴，以便从记录上读取对应射线的走时。同相轴在地震数据解释中起重要作用。现有的同相轴拾取方法有人工拾取方法及自动拾取方法等，人工拾取方法计算速度慢、精度低，且受到主观因素的影响；自动拾取方法有互相关法、模式识别法、神经网络法等，这些方法很难拾取连续同相轴，且抗噪性差。

发明内容

针对现有技术中存在的问题，本发明提供一种信号的信噪比高，边缘清晰准确，完整度及连续性好的，用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法。

本发明是通过以下技术方案来实现：

本发明用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，包括如下步骤，

步骤1，将叠前CSP数据用小波变换得到信号的解析部分；

步骤2，根据步骤1得到的解析部分得到信号的瞬时振幅；

步骤3，在瞬时振幅的剖面上进行局部线性Radon变换，拾取瞬时振幅的峰值振幅，得到倾斜叠加峰值振幅剖面；

步骤4，基于倾斜叠加峰值振幅剖面得到边缘检测剖面，用于在边缘检测剖面上拾取同相轴及走时。

优选的，步骤1中，将叠前CSP数据s(t)用小波变换得到信号的解析部分H[s(t)]如下，

$H\left[s\left(t\right)\right]=\mathrm{Im}[{\∫}_{-\∞}^{\∞}S(b,a)a^{-1}\mathrm{da}/{\∫}_0^{\∞}{\hat{g}}_R\left(\mathrm{\ω}\right){\mathrm{\ω}}^{-1}\mathrm{d\ω}];$

其中，S(b,a)是地震道集数据的小波变换，为小波函数g(t)的Fourier变换的实部，a为尺度因子，b为平移因子。

进一步，s(t)关于g(t)的小波变换定义为：

$S(b,a)=\frac{1}{a}{\∫}_{-\∞}^{+\∞}s\left(t\right)\stackrel{\‾}{g}\left(\frac{t-b}{a}\right)\mathrm{dt};$

式中，t,b∈R，a＞0；g(t)∈L¹(R,dt)∩L²(R,dt)，是g(t)的复共轭。

优选的，步骤2中，所述的瞬时振幅如下，

A(t)＝|s(t)+i·H[s(t)]|；

其中，H[s(t)]为用小波变换计算的信号的解析部分，A(t)为瞬时振幅的模，s(t)为叠前CSP数据。

优选的，步骤3中所述的局部线性Radon变换方法具体步骤如下，

步骤3.1，在瞬时振幅剖面上选定参考道；

步骤3.2，对于参考道上的某个时间截距τ_j，将参考道附近的几道沿n_p个具有不同斜率p_j(j＝1,2,…,n_p)的直线进行叠加，斜率以Δp为间隔采样；

步骤3.3，计算该时间截距处瞬时振幅沿不同方向叠加的和，将该叠加值记录在τ-p坐标轴相应的位置(τ_j,p_j)上，当选取的叠加斜率与同相轴的斜率接近或相等时，t-x域中的记录沿该直线的叠加值最大；

步骤3.4，将最大叠加值的平均值放置在t-x域中的对应位置(τ_j,x_m)上，能够构造一个超道集以增加信噪比，称该剖面为倾斜叠加峰值振幅剖面，对应于最大叠加值的斜率同样放置在t-x域的对应位置上，构成射线梯度剖面；其中，t为时间，x为距离域。

优选的，步骤3中，Radon变换频率域的离散形式为：

$M(f,p)=\underset{m=1}{\overset{\mathrm{nx}}{\mathrm{\Σ}}}\hat{A}(f,x_m)e^{j2\mathrm{\πfp}{x}_m}$

式中，M(f,p)＝∫m(τ,p)e^-j2πfτdτ，为瞬时振幅A(t,x_m)的傅里叶变换，x_m为参考道，m(τ,p)是时间域的Radon变换。

优选的，步骤4中得到的边缘检测剖面由仿真软件Matlab中的edge函数利用canny微分算子得到。与现有技术相比，本发明具有以下有益的技术效果：

本发明通过小波变换、Radon变换以及构造超道集途径提高了信号的信噪比，实现了同相轴的自动拾取。倾斜叠加变换即Radon变换，其作用类似滤波器，可以压制噪声。我们在瞬时振幅剖面上进行倾斜叠加变换，因为瞬时振幅剖面相较原始地震记录更加圆滑，包含更多的低频分量，且可滤除部分高频毛刺噪声，而二者的同相轴位置相同。有效压制了噪声，克服了现有方法难以在大偏移距内拾取连续同相轴的困难，且操作方便，实现了同相轴自动拾取。

进一步的，该边缘检测剖面由Matlab中的edge函数利用canny微分算子得到，用该算子提取的边缘清晰且精确，完整度及连接性较好，从边缘检测剖面上可拾取局部连续同相轴及对应的走时。

附图说明

图1本发明实例中所述的Radon变换示意图。

图2a本发明实例中所述的实际CSP道集。

图2b本发明实例中所述实际CSP道集的瞬时振幅剖面。

图3a本发明实例中所述的倾斜叠加峰值振幅剖面。

图3b本发明实例中所述倾斜叠加峰值振幅剖面的边缘检测剖面。

图4本发明实例中所述的同相轴拾取流程图。

具体实施方式

下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。

本发明用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，如图4所示，其包括以下步骤：

(1)、用小波变换计算得到信号的解析部分，其中叠前CSP数据作为原始信号；

(2)、根据解析部分计算得到信号的瞬时振幅：

A(t)＝|s(t)+i·H[s(t)]| (1)

其中， $H\left[s\left(t\right)\right]=\mathrm{Im}[{\∫}_{-\∞}^{\∞}S(b,a)a^{-1}\mathrm{da}/{\∫}_0^{\∞}{\hat{g}}_R\left(\mathrm{\ω}\right){\mathrm{\ω}}^{-1}\mathrm{d\ω}]$ 为用小波变换计算的信号的解析部分，A(t)为瞬时振幅的模，s(t)为叠前CSP数据。

(3)、在瞬时振幅剖面上计算局部线性Radon变换，拾取瞬时振幅的峰值振幅，得到倾斜叠加峰值振幅剖面。

其中，局部Radon变换即倾斜叠加，又称为τ-p变换，其频率域的离散形式为：

$M(f,p)=\underset{m=1}{\overset{\mathrm{nx}}{\mathrm{\Σ}}}\hat{A}(f,x_m)e^{j2\mathrm{\πfp}{x}_m}---\left(2\right)$

式中，M(f,p)＝∫m(τ,p)e^-j2πfτdτ，为瞬时振幅A(t,x_m)的傅里叶变换，x_m为参考道，m(τ,p)是时间域的Radon变换。

(4)、基于倾斜叠加峰值振幅剖面计算边缘检测剖面，用于在边缘检测剖面上拾取同相轴及走时。该边缘检测剖面由Matlab中的edge函数利用canny微分算子得到，可以看出，用该算子提取的边缘清晰且精确，完整度及连接性较好，从边缘检测剖面上可拾取局部连续同相轴及对应的走时。该同相轴拾取方法有效压制了噪声，克服了现有方法难以在大偏移距内拾取连续同相轴的困难，且操作方便，实现了同相轴自动拾取。

具体步骤如图4所示。

瞬时振幅的计算公式如下：

A(t)＝|s(t)+i·H[s(t)]| (1)

其中， $H\left[s\left(t\right)\right]=\mathrm{Im}[{\∫}_{-\∞}^{\∞}S(b,a)a^{-1}\mathrm{da}/{\∫}_0^{\∞}{\hat{g}}_R\left(\mathrm{\ω}\right){\mathrm{\ω}}^{-1}\mathrm{d\ω}]$ 为用小波变换计算的信号的解析部分，A(t)为瞬时振幅的模，s(t)为叠前CSP数据，S(b,a)是地震道集数据的小波变换，也就是CSP数据的小波变换，为小波函数g(t)的Fourier变换的实部，a为尺度因子，b为平移因子，s(t)关于g(t)的小波变换定义为：

$S(b,a)=\frac{1}{a}{\∫}_{-\∞}^{+\∞}s\left(t\right)\stackrel{\‾}{g}\left(\frac{t-b}{a}\right)\mathrm{dt}---\left(2\right)$

式中，t,b∈R，a＞0；g(t)∈L¹(R,dt)∩L²(R,dt)，是g(t)的复共轭。以实际数据为例，图2a为某油田的共炮点道集(CSP)，图2b为对应于共炮点道集的瞬时振幅剖面。该共炮点道集共595道，最小偏移距为90m，相邻检波器之间的距离是10m。共炮点道集是同一炮激发、不同检波器接收的地震道的集合。从瞬时振幅剖面上可以看出，其信噪比优于原始剖面，且同相轴位置对应于原始剖面的同相轴位置。

倾斜叠加即局部Radon变换。本发明利用局部线性Radon变换求取射线参数。线性Radon变换的积分路径是线性的，又称为τ-p变换，其频率域的离散形式为：

$M(f,p)=\underset{m=1}{\overset{\mathrm{nx}}{\mathrm{\Σ}}}\hat{A}(f,x_m)e^{j2\mathrm{\πfp}{x}_m}---\left(3\right)$

式中，M(f,p)＝∫m(τ,p)e^-j2πfτdτ，为瞬时振幅A(t,x_m)的傅里叶变换，x_m为参考道，m(τ,p)是时间域的Radon变换。局部线性Radon变换的计算过程如图1所示：在瞬时振幅剖面上选定参考道，对于参考道上的某个时间截距τ_j，将参考道附近的几道沿n_p个具有不同斜率p_j(j＝1,2,…,n_p)的直线(斜率以Δp为间隔采样)进行叠加，计算该时间截距处瞬时振幅沿不同方向叠加的和，将该叠加值记录在τ-p坐标轴相应的位置(τ_j,p_j)上，当选取的叠加斜率与同相轴的斜率接近或相等时，t-x域中的记录沿该直线的叠加值最大。将最大叠加值的平均值放置在t-x域(时间-距离域)中的对应位置(τ_j,x_m)上，可以构造一个超道集以增加信噪比，称该剖面为倾斜叠加峰值振幅剖面，对应于最大叠加值的斜率同样放置在t-x域的对应位置上，构成射线梯度剖面。线性Radon变换的计算方法有很多种，本发明采用高精度频率-空间域矩阵相乘法计算。进行Radon变换前，将时间域的瞬时振幅数据补零以增加频率分辨率，本发明中补零长度为数据长度的三倍。

由局部Radon变换得到的倾斜叠加峰值振幅剖面如图3a所示，由该剖面可以看出，与瞬时振幅剖面相比，二者的同相轴位置相吻合，倾斜叠加峰值振幅剖面的同相轴更清晰，干扰被有效压制，这是由于当沿梯度方向叠加参考道邻域内的几道数据时，噪声被进一步过滤的缘故。

边缘检测用来识别重要属性和结构变化，它能够剔除不相关的数据，减少数据信息量，是特征提取中重要的研究方向。图3b是基于倾斜叠加峰值振幅剖面得到的边缘检测剖面。该边缘检测剖面由Matlab中的edge函数利用canny微分算子得到，可以看出，用该算子提取的边缘清晰且精确，完整度及连接性较好，从边缘检测剖面上可拾取局部连续同相轴及对应的走时。该同相轴拾取方法有效压制了噪声，克服了现有方法难以在大偏移距内拾取连续同相轴的困难，且操作方便，实现了同相轴自动拾取。

Claims (6)

1.用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，其特征在于，包括如下步骤，

步骤1，将叠前CSP数据用小波变换得到信号的解析部分；

步骤2，根据步骤1得到的解析部分得到信号的瞬时振幅；

步骤3，在瞬时振幅的剖面上进行局部线性Radon变换，拾取瞬时振幅的峰值振幅，得到倾斜叠加峰值振幅剖面；

步骤4，基于倾斜叠加峰值振幅剖面得到边缘检测剖面，用于在边缘检测剖面上拾取同相轴及走时；

步骤3中所述的局部线性Radon变换方法具体步骤如下，

步骤3.1，在瞬时振幅剖面上选定参考道；

步骤3.2，对于参考道上的某个时间截距τ_j，将参考道附近的几道沿n_p个具有不同斜率p_j(j＝1,2,…,n_p)的直线进行叠加，斜率以Δp为间隔采样；

步骤3.3，计算该时间截距处瞬时振幅沿不同方向叠加的和，将叠加值记录在τ-p坐标轴相应的位置(τ_j,p_j)上，当选取的叠加斜率与同相轴的斜率接近或相等时，t-x域中的记录沿该直线的叠加值最大；

步骤3.4，将最大叠加值的平均值放置在t-x域中的对应位置(τ_j,x_m)上，能够构造一个超道集以增加信噪比，称该剖面为倾斜叠加峰值振幅剖面，对应于最大叠加值的斜率同样放置在t-x域的对应位置上，构成射线梯度剖面；其中，t为时间，x为距离域。

2.根据权利要求1所述的用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，其特征在于，步骤1中，将叠前CSP数据s(t)用小波变换得到信号的解析部分H[s(t)]如下，

$H\[s\left(t\right)\]=\mathrm{Im}\[{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}S(b,a)a^{-1}da/{\∫}_0^{\mathrm{\∞}}{\hat{g}}_R\left(\mathrm{\ω}\right){\mathrm{\ω}}^{-1}d\mathrm{\ω}\];$

其中，S(b,a)是地震道集数据的小波变换，为小波函数g(t)的Fourier变换的实部，a为尺度因子，b为平移因子。

3.根据权利要求2所述的用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，其特征在于，s(t)关于g(t)的小波变换定义为：

$S(b,a)=\frac{1}{a}{\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}s\left(t\right)\stackrel{\‾}{g}\left(\frac{t-b}{a}\right)dt;$

式中，t,b∈R，a＞0；g(t)∈L¹(R,dt)∩L²(R,dt)，是g(t)的复共轭。

4.根据权利要求1所述的用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，其特征在于，步骤2中，所述的瞬时振幅如下，

A(t)＝|s(t)+i·H[s(t)]|；

其中，H[s(t)]为用小波变换计算的信号的解析部分，A(t)为瞬时振幅的模，s(t)为叠前CSP数据。

5.根据权利要求1所述的用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，其特征在于，步骤3中，Radon变换频率域的离散形式为：

$M(f,p)=\underset{m=1}{\overset{nx}{\Σ}}\hat{A}(f,x_m)e^{j2{\mathrm{\πfpx}}_m}$

式中，M(f,p)＝∫m(τ,p)e^-j2πfτdτ，为瞬时振幅A(t,x_m)的傅里叶变换，x_m为参考道，m(τ,p)是时间域的Radon变换。

6.根据权利要求1所述的用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，其特征在于，步骤4中得到的边缘检测剖面由仿真软件Matlab中的edge函数利用canny微分算子得到。