说明书用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法
技术领域
本发明涉及地震勘探技术领域，具体为一种用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法。
背景技术
走时是衰减和速度层析成像的重要参数之一，我们需要拾取每个同相轴不同偏移距处反射波的走时，包括大偏移距处的同相轴。可以看出，走时拾取的前提是同相轴拾取，同相轴拾取也是地震解释的重要步骤之一。但由于噪声的影响，大偏移距处道集的信噪比低，拾取连续同相轴比较困难。为了拾取正确的反射波走时，我们需要在大偏移距范围内拾取连续同相轴，以便从记录上读取对应射线的走时。同相轴在地震数据解释中起重要作用。现有的同相轴拾取方法有人工拾取方法及自动拾取方法等，人工拾取方法计算速度慢、精度低，且受到主观因素的影响；自动拾取方法有互相关法、模式识别法、神经网络法等，这些方法很难拾取连续同相轴，且抗噪性差。
发明内容
针对现有技术中存在的问题，本发明提供一种信号的信噪比高，边缘清晰准确，完整度及连续性好的，用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法。
本发明是通过以下技术方案来实现：
本发明用于同相轴自动拾取的倾斜叠加峰值振幅处边缘检测法，包括如下步骤，
步骤1，将叠前CSP数据用小波变换得到信号的解析部分；
步骤2，根据步骤1得到的解析部分得到信号的瞬时振幅；
步骤3，在瞬时振幅的剖面上进行局部线性Radon变换，拾取瞬时振幅的峰值振幅，得到倾斜叠加峰值振幅剖面；
步骤4，基于倾斜叠加峰值振幅剖面得到边缘检测剖面，用于在边缘检测剖面上拾取同相轴及走时。
优选的，步骤1中，将叠前CSP数据s(t)用小波变换得到信号的解析部分H[s(t)]如下，
H[s(t)]=Im[∫-∞∞S(b,a)a-1da/∫0∞g^R(ω)ω-1dω];]]>
其中，S(b,a)是地震道集数据的小波变换，为小波函数g(t)的Fourier变换的实部，a为尺度因子，b为平移因子。
进一步，s(t)关于g(t)的小波变换定义为：
S(b,a)=1a∫-∞+∞s(t)g‾(t-ba)dt;]]>
式中，t,b∈R，a＞0；g(t)∈L1(R,dt)∩L2(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，将参考道附近的几道沿np个具有不同斜率pj(j＝1,2,…,np)的直线进行叠加，斜率以Δp为间隔采样；
步骤3.3，计算该时间截距处瞬时振幅沿不同方向叠加的和，将该叠加值记录在τ-p坐标轴相应的位置(τj,pj)上，当选取的叠加斜率与同相轴的斜率接近或相等时，t-x域中的记录沿该直线的叠加值最大；
步骤3.4，将最大叠加值的平均值放置在t-x域中的对应位置(τj,xm)上，能够构造一个超道集以增加信噪比，称该剖面为倾斜叠加峰值振幅剖面，对 应于最大叠加值的斜率同样放置在t-x域的对应位置上，构成射线梯度剖面；其中，t为时间，x为距离域。
优选的，步骤3中，Radon变换频率域的离散形式为：
M(f,p)=Σm=1nxA^(f,xm)ej2πfpxm]]>
式中，M(f,p)＝∫m(τ,p)e-j2πfτdτ，为瞬时振幅A(t,xm)的傅里叶变换，xm为参考道，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[s(t)]=Im[∫-∞∞S(b,a)a-1da/∫0∞g^R(ω)ω-1dω]]]>为用小波变换计算的信号的解析部分，A(t)为瞬时振幅的模，s(t)为叠前CSP数据。
(3)、在瞬时振幅剖面上计算局部线性Radon变换，拾取瞬时振幅的峰值振幅，得到倾斜叠加峰值振幅剖面。
其中，局部Radon变换即倾斜叠加，又称为τ-p变换，其频率域的离散形式为：
M(f,p)=Σm=1nxA^(f,xm)ej2πfpxm---(2)]]>
式中，M(f,p)＝∫m(τ,p)e-j2πfτdτ，为瞬时振幅A(t,xm)的傅里叶变换，xm为参考道，m(τ,p)是时间域的Radon变换。
(4)、基于倾斜叠加峰值振幅剖面计算边缘检测剖面，用于在边缘检测剖面上拾取同相轴及走时。该边缘检测剖面由Matlab中的edge函数利用canny微分算子得到，可以看出，用该算子提取的边缘清晰且精确，完整度 及连接性较好，从边缘检测剖面上可拾取局部连续同相轴及对应的走时。该同相轴拾取方法有效压制了噪声，克服了现有方法难以在大偏移距内拾取连续同相轴的困难，且操作方便，实现了同相轴自动拾取。
具体步骤如图4所示。
瞬时振幅的计算公式如下：
A(t)＝|s(t)+i·H[s(t)]| (1)
其中，H[s(t)]=Im[∫-∞∞S(b,a)a-1da/∫0∞g^R(ω)ω-1dω]]]>为用小波变换计算的信号的解析部分，A(t)为瞬时振幅的模，s(t)为叠前CSP数据，S(b,a)是地震道集数据的小波变换，也就是CSP数据的小波变换，为小波函数g(t)的Fourier变换的实部，a为尺度因子，b为平移因子，s(t)关于g(t)的小波变换定义为：
S(b,a)=1a∫-∞+∞s(t)g‾(t-ba)dt---(2)]]>
式中，t,b∈R，a＞0；g(t)∈L1(R,dt)∩L2(R,dt)，是g(t)的复共轭。以实际数据为例，图2a为某油田的共炮点道集(CSP)，图2b为对应于共炮点道集的瞬时振幅剖面。该共炮点道集共595道，最小偏移距为90m，相邻检波器之间的距离是10m。共炮点道集是同一炮激发、不同检波器接收的地震道的集合。从瞬时振幅剖面上可以看出，其信噪比优于原始剖面，且同相轴位置对应于原始剖面的同相轴位置。
倾斜叠加即局部Radon变换。本发明利用局部线性Radon变换求取射线参数。线性Radon变换的积分路径是线性的，又称为τ-p变换，其频率域的离散形式为：
M(f,p)=Σm=1nxA^(f,xm)ej2πfpxm---(3)]]>
式中，M(f,p)＝∫m(τ,p)e-j2πfτdτ，为瞬时振幅A(t,xm)的傅里叶变换，xm为参考道，m(τ,p)是时间域的Radon变换。局部线性Radon变换的计算过程如图1所示：在瞬时振幅剖面上选定参考道，对于参考道上的某个时间截距τj，将参考道附近的几道沿np个具有不同斜率pj(j＝1,2,…,np)的直线(斜率以Δp为间隔采样)进行叠加，计算该时间截距处瞬时振幅沿不同方向叠加的和，将该叠加值记录在τ-p坐标轴相应的位置(τj,pj)上，当选取的叠加斜率与同相轴的斜率接近或相等时，t-x域中的记录沿该直线的叠加值最大。将最大叠加值的平均值放置在t-x域(时间-距离域)中的对应位置(τj,xm)上，可以构造一个超道集以增加信噪比，称该剖面为倾斜叠加峰值振幅剖面，对应于最大叠加值的斜率同样放置在t-x域的对应位置上，构成射线梯度剖面。线性Radon变换的计算方法有很多种，本发明采用高精度频率-空间域矩阵相乘法计算。进行Radon变换前，将时间域的瞬时振幅数据补零以增加频率分辨率，本发明中补零长度为数据长度的三倍。
由局部Radon变换得到的倾斜叠加峰值振幅剖面如图3a所示，由该剖面可以看出，与瞬时振幅剖面相比，二者的同相轴位置相吻合，倾斜叠加峰值振幅剖面的同相轴更清晰，干扰被有效压制，这是由于当沿梯度方向叠加参考道邻域内的几道数据时，噪声被进一步过滤的缘故。
边缘检测用来识别重要属性和结构变化，它能够剔除不相关的数据，减少数据信息量，是特征提取中重要的研究方向。图3b是基于倾斜叠加峰值振幅剖面得到的边缘检测剖面。该边缘检测剖面由Matlab中的edge函数利用canny微分算子得到，可以看出，用该算子提取的边缘清晰且精确，完整度及连接性较好，从边缘检测剖面上可拾取局部连续同相轴及对应的走时。 该同相轴拾取方法有效压制了噪声，克服了现有方法难以在大偏移距内拾取连续同相轴的困难，且操作方便，实现了同相轴自动拾取。