CN104062681A - 一种基于分数阶导数的地震层位追踪预处理方法 - Google Patents

Abstract

本发明提供了一种基于分数阶导数的地震层位追踪预处理方法，属于石油天然气地球物理勘探领域。所述方法包括以下步骤：(1)选取地震图像的中心像素点，然后确定分数阶微分的阶数；(2)选取各向同性扩散滤波方向；(3)对中心像素点进行同相轴连续性增强；(4)取下一个中心像素点，重复步骤(1)至步骤(3)，直到处理完所有像素点为止；(5)进行各向同性扩散滤波处理；(6)对地震图像的纵向分辨率进行调整，至此就完成了地震层位追踪预处理。利用本发明方法能同时控制同相轴横向上连续性以及纵向上的分辨率，为地震层位自动追踪提供了一种简单、计算效率高、效果明显的预处理方法。

Description

一种基于分数阶导数的地震层位追踪预处理方法

技术领域

本发明属于石油天然气地球物理勘探领域，具体涉及一种基于分数阶导数的地震层位追踪预处理方法。

背景技术

层位追踪在地震资料解释中有着举足轻重的作用，其准确性对后续地震资料的处理和解释起着重要作用。地震层位追踪可以简单地分为人工追踪，自动追踪，以及人工半自动追踪，近年来商业地球物理软件以及一些地球物理工作者研发了地震层位追踪方法，最简单的方法是按照地震波形的波峰或者波谷追踪，其他层位追踪的方法如《基于高阶统计量的层位自动追踪技术》(彭文，熊晓军等，新疆石油地质，2006年12月)，甚至一些商业地球物理软件中出现了三维的层位自动追踪技术。所有这些自动追踪方法目的是为了减轻繁重的人工层位追踪，当地震同相轴出现严重不连续的时候，自动追踪层位的结果会出现错误，目前折中的做法是用人工干预的半自动层位追踪方法跳过不连续点，但是实际生产中发现若同相轴连续性太差半自动追踪只能改为人工进行，极大耗费解释人员的精力，因此欲提高层位自动追踪的精度必须在一定程度上增强同相轴的连续性。

为了提高层位自动追踪的精度，必须增强地震同相轴的连续性，常规的同相轴连续性增强方法是使用高斯平滑滤波器，这是一种低通滤波器，削弱地震的高频成分，让低频成分相对增多来增强同相轴的连续性，这种方法由于削弱了地震信号的高频成分使得地震图像在纵向上分辨率降低，会降低层位自动追踪的精度，甚至出现“串层”现象。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种基于分数阶导数的地震层位追踪预处理方法，能同时控制同相轴横向上连续性以及纵向上的分辨率，为地震层位自动追踪提供一种方法简单、计算效率高、效果明显的预处理方法。

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

一种基于分数阶导数的地震层位追踪预处理方法，包括以下步骤：

(1)选取地震图像的中心像素点，然后确定分数阶微分的阶数；

(2)选取各向同性扩散滤波方向；

(3)对中心像素点进行同相轴连续性增强；

(4)取下一个中心像素点，重复步骤(1)至步骤(3)，直到处理完所有像素点为止；

(5)进行各向同性扩散滤波处理；

(6)对地震图像的纵向分辨率进行调整，至此就完成了地震层位追踪预处理。

所述步骤(1)中选取地震图像的中心像素点是这样实现的：

地震图像的像素点是按照行和列排列，使用时窗依次在地震图像上从左到右滑动，滑动的步长为两个像素点之间的间距，当时窗滑到某一行的最右边时，时窗移到地震图像的下一行像素点，每个时窗内的中心点像素就是要选取的中心像素点。

所述时窗的大小是3*3或5*5或7*7。

所述阶数是小于1的。

所述步骤(2)中选取各向同性扩散滤波方向是这样实现的：

取关于中心像素点对称的8个方向，即 $\mathrm{\θ}=0,\frac{\mathrm{\π}}{4},\frac{\mathrm{\π}}{2},\frac{3\mathrm{\π}}{4},\mathrm{\π},\frac{5\mathrm{\π}}{4},\frac{3\mathrm{\π}}{2},\frac{7\mathrm{\π}}{4}.$

所述步骤(3)是这样实现的：

对中心像素点按照(8)式进行同相轴连续性增强：

u(x，y)＝ηu₀(x，y)+Δu₀(x，y) (8)

其中，u(x，y)为增强后的中心像素点，u₀(x，y)为图像增强前的中心像素点值，▽u₀(x，y)为图像增量，η为调解因子，是由用户给出的，其范围为[0 1]。

所述图像增量▽u₀(x，y)是利用(7)式计算得到的：

$\mathrm{\Δu}(x,y)=\underset{\mathrm{\θ}=0}{\overset{2\mathrm{\π}}{\mathrm{\Σ}}}{u}_{\mathrm{\θ}}^{\mathrm{\γ}}(x,y)...\left(7\right)$

其中 $\mathrm{\θ}=0,\frac{\mathrm{\π}}{4},\frac{\mathrm{\π}}{2},\frac{3\mathrm{\π}}{4},\mathrm{\π},\frac{5\mathrm{\π}}{4},\frac{3\mathrm{\π}}{2},\frac{7\mathrm{\π}}{4},$ 为沿着θ方向的γ阶偏微分。

所述步骤(4)中是依次将所有地震像素点作为中心像素点进行处理的。

所述步骤(5)是这样实现的：

由用户选择扩散时间，然后按照(5)式进行各向同性扩散滤波处理，如果地震同相轴的连续性没有达到用户的要求，则扩大扩散时间，重复步骤(5)，如果地震同相轴的连续性达到用户的要求，则停止各向同性扩散滤波处理；

$\left\{\begin{array}{c}\frac{\∂u}{\∂t}=\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2}\\ u(x,y,0)=u\left(t_0\right)\end{array}\right....\left(5\right).$

所述步骤(6)是这样实现的：

按照(9)式对地震图像的纵向分辨率进行调整：

$u^{\′}(x,y)=\underset{\mathrm{\θ}\≠0,\mathrm{\π}}{\mathrm{\Σ}}F\left[u(x,y)\right]...\left(9\right)$

其中F为算子，F＝2f(n)-αf(n-1)-βf(n+1)，f(n)为地震信号，α，β为调解因子，是由用户给出的，其范围在[0，1]之间。

与现有技术相比，本发明的有益效果是：利用本发明方法能同时控制同相轴横向上连续性以及纵向上的分辨率，为地震层位自动追踪提供了一种简单、计算效率高、效果明显的预处理方法。

附图说明

图1是不同阶数偏微分方程的图像增强方法对应的幅频特征。

图2是偏微分掩模示意图。

图3是本发明实施例中的地震剖面及层位追踪结果。

图4图3地震剖面层位追踪结果的单独显示。

图5采用本发明方法对图3中的地震剖面处理效果及层位追踪结果。

图6是图5地震剖面层位追踪结果的单独显示。

具体实施方式

下面结合附图对本发明作进一步详细描述：本发明是一种半自动层位追踪中的地震预处理技术，预处理的目的是为了增强地震同相轴的横向连续性，提高层位追踪的精度。为了达到此目的，本发明引入分数阶微分的图像处理方法，非线性地增加同相轴的低频成分，同时借鉴各向同性扩散滤波的思想使得基于分数阶微分的滤波方法能够连续地进行直到满足要求为止，扩散滤波完成后地震资料的横向分辨率下降这对层位追踪不利，为此本发明用一个改进后的局部高频滤波器控制地震剖面的纵向分辨率。实际地震资料表明利用本发明方法，大大改善了地震剖面的横向连续性，提高了层位追踪的精度，能够大幅度减轻人工追踪层位的劳动强度。

具体如下：

(1)基于分数阶微分的增强地震同相轴连续性方法

近年来基于偏微分方程的图像处理技术越来越引起人们的重视，地震数据处理技术也引进了这种技术，比如基于二阶微分的拉普拉斯图像增强技术，再如J.weickert的各向异性扩散滤波技术广泛应用于地震图像的保边去噪，上述这些方法都基于整数阶微分，在数学分支—碎形学中将整数阶微分推广到分数阶微分，使用分数阶微增强图像会产生一些不同于整数阶的结果，以下是分别选取0.3阶，0.7阶，1阶，2阶偏微分说明其在地震图像增强处理中对地震频率成分的不同影响。仅以时间方向上的信号为例。设地震图像数据为u(t)，对应的分数阶微分为：

$\frac{d^{\mathrm{\γ}}u\left(t\right)}{d{t}^{\mathrm{\γ}}}\≈u\left(t\right)+(-\mathrm{\γ})u(t-1)+\frac{(-\mathrm{\γ})(-\mathrm{\γ}+1)}{2}u(t-2)+....\frac{\mathrm{\Γ}(n-1-\mathrm{\γ})}{(n-1)!\mathrm{\Γ}(-\mathrm{\γ})}...\left(1\right)$

其中Γ为伽马函数，表达式为：

$\mathrm{\Γ}\left(\mathrm{\α}\right)=\underset{0}{\overset{+\∞}{\∫}}{e}^{-x}{x}^{\mathrm{\α}-1}\mathrm{dx}=(\mathrm{\α}-1)!...\left(2\right)$

时间域偏微分算子与频率域的对应关系为偏微分算子的图像信号增强频率域形式为：

Δu(w)＝u(w)w^r (3)

图1为不同阶数偏微分方程的图像增强方法对应的幅频特征，从图1中可以看出当偏微分的阶数为分数时，基于偏微分的处理方法表明图像的低频和高频都会非线性增加，但是低频增强的快，高频增加的慢。而当阶数大于1的时候低频和高频成分仍然非线性地增加，但是高频增加的快，低频增加的慢。

进行地震层位自动追踪的时候，在不考虑断层的情况下同相轴的连续性越好，追踪的结果越理想，当地震资料横向连续性变差时追踪的结果变差，因此对于层位的自动追踪同相轴的连续性很重要，生产实际中发现当地震信号中的低频成分相对增加的时候同相轴的连续性会变好，有上述分析可知分数阶微分对于低频增加有利，因此基于分数阶微分的图像增强方法一种增强同相轴连续的有效方法。

(2)基于分阶微分的各向同性扩散滤波方法

扩散滤波的思想使得滤波过程可以连续进行直到达到满意为止，Weickert提出的各向异性扩散滤波方程为：

$\left\{\begin{array}{c}\frac{\∂u}{\∂t}=\mathrm{div}(\mathrm{fD}\▿u)\\ u(x,y,0)=u\left(t_0\right)\end{array}\right....\left(4\right)$

这里把断层因子f以及扩散张量D去掉使其退化为增强同相轴连续性的各向同性扩散滤波方法，相应的方程为：

$\left\{\begin{array}{c}\frac{\∂u}{\∂t}=\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2}\\ u(x,y,0)=u\left(t_0\right)\end{array}\right....\left(5\right)$

将偏微分的二阶导数换为分数阶微分得：

$\left\{\begin{array}{c}\frac{\∂u}{\∂t}=\frac{{\∂}^{\mathrm{\γ}}u}{\∂x^{\mathrm{\γ}}}+\frac{{\∂}^{\mathrm{\γ}}u}{\∂y^{\mathrm{\γ}}}\\ u(x,y,0)=u\left(t_0\right)\end{array}\right....\left(6\right)$

实际数值计算中选择8方向的微分掩模设图像中心点的灰度值为u(x，y)图像灰度值的增量

$\mathrm{\Δu}(x,y)=\underset{\mathrm{\θ}=0}{\overset{2\mathrm{\π}}{\mathrm{\Σ}}}{u}_{\mathrm{\θ}}^{\mathrm{\γ}}(x,y)...\left(7\right)$

其中 $\mathrm{\θ}=0,\frac{\mathrm{\π}}{4},\frac{\mathrm{\π}}{2},\frac{3\mathrm{\π}}{4},\mathrm{\π},\frac{5\mathrm{\π}}{4},\frac{3\mathrm{\π}}{2},\frac{7\mathrm{\π}}{4}$ 为沿着θ方向的γ阶偏微分。图2为偏微分掩模示意图。最终地震同相轴连续性增强的方程为：

u(x，y)＝ηu₀(x，y)+Δu₀(x，y) (8)

u(x，y)为增强后的像素点，u₀(x，y)为图像增强前的像素点值，▽u₀(x，y)为图像增量，η为调解因子范围为[0 1]，用此因子来调解图像时的高频和低频成分，以兼顾地震同相轴的连续性和时间分辨率。

(3)纵向分辨率处理

当基于分数阶的各向同性扩散滤波完成后，由于在8个方向上进行的增加低频成分，时间方向上分辨率必然下降，这对层位追踪的精度非常不利，此时需要对所有道地震信号时间方向上“压缩”，不失横向上的连续性，实验证明算子F＝2f(n)-αf(n-1)-βf(n+1)f(n)为地震信号，α，β为调解因子范围在[0，1]之间，设中心像素点为u(x，y)，将滤波算子作用在除了水平方向上的其他方向上即算子方向为：

$\mathrm{\θ}=\frac{\mathrm{\π}}{4},\frac{\mathrm{\π}}{2},\frac{3\mathrm{\π}}{4},\frac{5\mathrm{\π}}{4},\frac{3\mathrm{\π}}{2},\frac{7\mathrm{\π}}{4},$

作用后的像素点 $u^{\′}(x,y)=\underset{\mathrm{\θ}\≠0,\mathrm{\π}}{\mathrm{\Σ}}F\left[u(x,y)\right]...\left(9\right)$

(4)本发明方法的步骤

本发明方法包括以下步骤：

(1)选取地震图像的中心像素点，然后确定分数阶微分的阶数，本发明取的阶数要小于1，比如取0.3，0.7；

具体来说，选取地震图像的中心像素点是这样实现的：

地震图像的像素点是按照行和列排列，本实施例中使用5×5(但并不拘泥于5*5，一般建议3*3，5*5，7*7)的时窗依次在地震图像上从左到右滑动，滑动的步长为两个像素点之间的间距，当时窗滑到某一行的最右边时，时窗移到地震图像的下一行像素点，每个时窗内的中心点像素就是要选取的中心像素点。

(2)选取各向同性扩散滤波方向，各个方向是关于中心像素点对称的，本发明中取8个方向即 $\mathrm{\θ}=0,\frac{\mathrm{\π}}{4},\frac{\mathrm{\π}}{2},\frac{3\mathrm{\π}}{4},\mathrm{\π},\frac{5\mathrm{\π}}{4},\frac{3\mathrm{\π}}{2},\frac{7\mathrm{\π}}{4};$

(3)对中心像素点按照(8)式进行同相轴连续性增强；

u(x，y)＝ηu₀(x，y)+Δu₀(x，y) (8)

其中，u(x，y)为增强后的中心像素点，u₀(x，y)为图像增强前的中心像素点值，▽u₀(x，y)为图像增量(该增量是利用(7)式计算得来的)，η为调解因子(这个参数可以选取0到1之间的任意数值，这个因子值大的时候有利于保持同相轴的时间分辨率，反之有利于增强同相轴的连续性，数值由用户自己给出)，其范围为[01]；

(4)取下一个中心像素点，重复步骤(1)-(3)，直到处理完所有像素点为止(所有地震像素点都依次被作为中心像素点)；

(5)当地震图像所有像素点处理完成后，选择扩散时间(扩散时间只能由用户根据经验来选取，一般起始时间为0时间步长为0.5，扩散时间1.5，即扩散三次)按照(5)式进行各向同性扩散滤波处理，直到地震同相轴的连续性达到要求为止(用户观察处理过的地震图像，用户感觉达到自己要求了就可以停止滤波了，一般扩散三次)；

$\left\{\begin{array}{c}\frac{\∂u}{\∂t}=\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2}\\ u(x,y,0)=u\left(t_0\right)\end{array}\right....\left(5\right)$

(5)式中，u为任意一个地震图像像素点的值，t为扩散时间，t₀为初始扩散时间，x，y分别为图像的横坐标纵坐标。

扩散时间选得越长连续性增强的效果越好但同时同相轴的时间分辨率会下降，因此要兼顾同相轴连续性与时间分辨率，因此只要逐步加大扩散时间并观察处理的地震图像效果，如果满足用户要求就停止对图像的处理，不满足就继续加大扩散时间再处理。

(6)完成扩散滤波后按照(9)式选择适当的参数调整(这个参数也是用户自己经验选择的，默认值α，β参数为0.9，0.9)地震图像的纵向分辨率：

作用后的像素点 $u^{\′}(x,y)=\underset{\mathrm{\θ}\≠0,\mathrm{\π}}{\mathrm{\Σ}}F\left[u(x,y)\right]...\left(9\right)$

其中F为算子，F＝2f(n)-αf(n-1)-βf(n+1)，f(n)为地震信号，α，β为调解因子，其范围在[0，1]之间，设中心像素点为u(x，y)(即增强后的中心像素点)。

选择α，β参数后将(9)式作用于扩散滤波后的地震图像就达到了调整分辨率的目的，至此就完成了地震层位追踪预处理。

图3为选择验证本发明效果的地震剖面，此地震剖面横向上连续性较差。在此剖面上选择一个种子点A，采用一种自动追踪算法，追踪的结果已经投影在了图3的地震剖面(如图3中的深黑线所示)，追踪的结果已经严重偏离了真实层位，在实际生产中这一部分的层位追踪必须依靠人工追踪了。图4为追踪出来的层位，曲线出现突变点，不连续。

图5为运用本发明方法处理后的地震剖面及其追踪出的层位，选取与图3中相同的种子点A，自动追踪算法也相同。对比真实的地震层位，大部分自动追踪出的层位是准确的，不准确的部分可以选择新的种子点追踪。图6展示了追踪出的层位，相比预处理前追踪的层位，准确、曲线连续。实际资料表明本发明的地震同相轴横向连续性增强方法，大大改善了地震资料的品质，使得半自动的层位追踪准确可靠，大大减轻了人工追踪层位的劳动强度。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (10)

1.一种基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述方法包括以下步骤：

(1)选取地震图像的中心像素点，然后确定分数阶微分的阶数；

(2)选取各向同性扩散滤波方向；

(3)对中心像素点进行同相轴连续性增强；

(4)取下一个中心像素点，重复步骤(1)至步骤(3)，直到处理完所有像素点为止；

(5)进行各向同性扩散滤波处理；

(6)对地震图像的纵向分辨率进行调整，至此就完成了地震层位追踪预处理。

2.根据权利要求1所述的基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述步骤(1)中选取地震图像的中心像素点是这样实现的：

地震图像的像素点是按照行和列排列，使用时窗依次在地震图像上从左到右滑动，滑动的步长为两个像素点之间的间距，当时窗滑到某一行的最右边时，时窗移到地震图像的下一行像素点，每个时窗内的中心点像素就是要选取的中心像素点。

3.根据权利要求2所述的基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述时窗的大小是3*3或5*5或7*7。

4.根据权利要求3所述的基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述阶数是小于1的。

5.根据权利要求4所述的基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述步骤(2)中选取各向同性扩散滤波方向是这样实现的：

取关于中心像素点对称的8个方向，即 $\mathrm{\θ}=0,\frac{\mathrm{\π}}{4},\frac{\mathrm{\π}}{2},\frac{3\mathrm{\π}}{4},\mathrm{\π},\frac{5\mathrm{\π}}{4},\frac{3\mathrm{\π}}{2},\frac{7\mathrm{\π}}{4}.$

6.根据权利要求5所述的基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述步骤(3)是这样实现的：

对中心像素点按照(8)式进行同相轴连续性增强：

u(x，y)＝ηu₀(x，y)+Δu₀(x，y) (8)

其中，u(x，y)为增强后的中心像素点，u₀(x，y)为图像增强前的中心像素点值，▽u₀(x，y)为图像增量，η为调解因子，是由用户给出的，其范围为[0 1]。

7.根据权利要求6所述的基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述图像增量▽u₀(x，y)是利用(7)式计算得到的：

$\mathrm{\Δu}(x,y)=\underset{\mathrm{\θ}=0}{\overset{2\mathrm{\π}}{\mathrm{\Σ}}}{u}_{\mathrm{\θ}}^{\mathrm{\γ}}(x,y)...\left(7\right)$

其中 $\mathrm{\θ}=0,\frac{\mathrm{\π}}{4},\frac{\mathrm{\π}}{2},\frac{3\mathrm{\π}}{4},\mathrm{\π},\frac{5\mathrm{\π}}{4},\frac{3\mathrm{\π}}{2},\frac{7\mathrm{\π}}{4},$ 为沿着θ方向的γ阶偏微分。

8.根据权利要求7所述的基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述步骤(4)中是依次将所有地震像素点作为中心像素点进行处理的。

9.根据权利要求8所述的基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述步骤(5)是这样实现的：

由用户选择扩散时间，然后按照(5)式进行各向同性扩散滤波处理，如果地震同相轴的连续性没有达到用户的要求，则扩大扩散时间，重复步骤(5)，如果地震同相轴的连续性达到用户的要求，则停止各向同性扩散滤波处理；

$\left\{\begin{array}{c}\frac{\∂u}{\∂t}=\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2}\\ u(x,y,0)=u\left(t_0\right)\end{array}\right....\left(5\right).$

10.根据权利要求9所述的基于分数阶导数的地震层位追踪预处理方法，其特征在于：所述步骤(6)是这样实现的：

按照(9)式对地震图像的纵向分辨率进行调整：

$u^{\′}(x,y)=\underset{\mathrm{\θ}\≠0,\mathrm{\π}}{\mathrm{\Σ}}F\left[u(x,y)\right]...\left(9\right)$

其中F为算子，F＝2f(n)-αf(n-1)-βf(n+1)，f(n)为地震信号，α，β为调解因子，是由用户给出的，其范围在[0，1]之间。