CN106249288A - 基于Shearlet域的极化滤波面波压制方法 - Google Patents

Abstract

本发明涉及一种基于Shearlet域的极化滤波面波压制方法。本发明方法在shearlet域内进行极化分析，可以使得不同视速度的面波与反射波分布在不同的方向下，并且在各方向下进行极化分析，可以避免假频面波的干扰。本发明的方法对不同方向的shearlet系数采用不同的椭球率阈值，从而在反射波的方向上，椭球率阈值偏大，尽可能的保留反射波，而在面波的方向上椭球率阈值偏小，尽可能的滤除面波。

Description

基于Shearlet域的极化滤波面波压制方法

技术领域

本发明涉及地震勘探领域，是地震数据处理的一个环节,具体为一种基于Shearlet域的极化滤波面波压制方法。

背景技术

在常规地震勘探中，面波是一种常见的规则干扰波，面波能量通常较强，它的存在严重降低了地震剖面的信噪比。低频、低视速度以及椭圆极化是面波的主要特性。据此，目前主要的面波压制的技术有：频率域滤波、视速度滤波、极化滤波。

一、视速度滤波

视速度滤波就是基于面波与有效反射波之间的视速度差异，将地震数据换到另外一个域内，然后在变换域内能够较方便地得到视速度信息，从而可以利用二者的差异进行滤波。FK变换是一种较常用的视速度滤波，将地震数据变换到FK域中，FK正反变换如式(1)和式(2)所示：

$D(f,k)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}d(t,x)e^{-ift-ikx}dtdx---\left(1\right);$

$d(t,x)=\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}\underset{-\mathrm{\∞}}{\overset{+\mathrm{\∞}}{\∫}}D(f,k)e^{ift+ikx}dfdk---\left(2\right);$

视速度可以表达为：

$v=\frac{f}{k}---\left(3\right);$

则可以合理的对数据进行滤波；

式(1)—(3)中，f为频率，k为波数。

例如，如图1所示，一个二维地震数据，横向为偏移距(可以理解为空间域)，纵向为时间。可见地震数据内较为平缓的弧线为需要保留的有效信息，它的速度较大，因此可以在较短的时间传播到其它空间点上；而面波集中在中间部分，视觉上比较陡，这是因为它的速度比较小，因此在一定的时间只能传播到附近的空间点上。若对该地震数据进行FK变换，如图2所示，根据视速度的表达式(3)，底部弯折的线条视速度值较低，为需要滤除的面波部分，而中间的为需要保留的反射波部分。可见此时有效信息与面波的重叠相对于图1减少了不少，有利于二者的分离，但是稍微仔细的观察可以发现折叠的面波(底部弯折部分，该现象在地震勘探中称为假频)与有效波依然存在一些重叠，因此不能完全的将二者分离。对面波进行切除后的FK谱如图3所示，可见圈内仍有少部分面波残留，由于该部分与反射波重叠，因此无法完全去除。图4为滤波后的地震数据，大部分面波还是被滤除了，而少部分假频面波能量仍然残留，如图4中圆形区域内所示。

上述例子可见，当面波存在假频(底部FK能量折叠)的时候，假频面波与反射波会存在重叠，因此在FK域无法将二者分离。滤波后会有如图4所示的假频面波的残留。

二、极化滤波

极化滤波目前已经发展到了时频域极化滤波，该方法不受假频面波的影响。

1、极化滤波基本概念

同一类型的地震波至引起质点的运动，其轨迹类似于椭球体。如图5，根据极化的概念，这个椭球的长轴v₁即是粒子运动的极化方向。不同的波在空间的质点运动轨迹是不同的，利用多分量记录，在以t时刻为中心，时窗长度为T₀的滑动时窗中构建协方差矩阵，计算出描述波的粒子运动轨迹的特征参数。利用这些参数可以设计出合适的极化滤波函数，保留那些满足特定要求的质点运动(葛勇，韩立国，韩文明，等.1996.极化分析研究及其在波场分离中的应用[J].长春地质学院学报，26(1)：83-88.黄中玉，高林，徐亦鸣，等.1996.三分量数据的偏振分析及其应用[J].石油物探，35(2)：10-16.)。

极化滤波最早由Flinn提出的时间域固定时窗极化滤波(Flinn E A,1965,Signalanalysis using rectilinearity and direction of particle motion[J].Proceedingsof the IEEE.53(12):1874-1876)，发展到Diallo等人提出的时间域自适应极化滤波，再到后来Pinnegar、Diallo等人提出的时频域极化滤波。

2、极化滤波计算式

设s(t)为地震信号，选定分析时窗长度T₀，则时窗内的协方差矩阵可由式(4)—(6)表达：

$M\left(t\right)=\left[\begin{array}{ccc}{I}_{xx}\left(t\right)& I_{xy}\left(t\right)& I_{xz}\left(t\right)\\ I_{yx}\left(t\right)& I_{yy}\left(t\right)& I_{yz}\left(t\right)\\ I_{zx}\left(t\right)& I_{zy}\left(t\right)& I_{zz}\left(t\right)\end{array}\right]---\left(4\right);$

$I_{km}\left(t\right)=\frac{1}{T_0}{\∫}_{-T_0/2}^{T_0/2}\[s_k(t+\mathrm{\τ})-u_k\left(t\right)\]\[s_m(t+\mathrm{\τ})-u_m\left(t\right)\]d\mathrm{\τ},(k,m=x,y,z)---\left(5\right);$

$u_k\left(t\right)={\∫}_{-T_0/2}^{T_0/2}{s}_k(t+\mathrm{\τ})d\mathrm{\τ},(k=x,y,z)---\left(6\right);$

其中，M(t)即为所构建的协方差矩阵，u(t)(即u_k(t)和u_m(t))为时窗内信号的均值，I_km(t)为各分量之间的协方差。求解得到式(4)中M(t)的特征向量v₁(t)，v₂(t)，v₃(t)以及相对应特征值λ₁(t)＞λ₂(t)＞λ₃(t)，这分别代表着时窗内信号振动的三个主要方向和这三个方向的轴的大小，称v₁为主极化方向。当地震波是体波即纵波或横波的时候，会有λ₁(t)＞＞λ₂(t)＞λ₃(t)，即能量会主要集中在一个方向上。而若地震波为椭圆极化如面波，则会有λ₁(t)≈λ₂(t)＞λ₃(t)，即地震波的能量主要集中在两个方向上，且两个方向能量大小相近。据此，根据式(7)-(8)定义极化率或者椭球率来对极化轨迹形状作描述，从而可以对面波与反射波进行分离。

$e\left(t\right)=\frac{{\mathrm{\λ}}_2\left(t\right)}{{\mathrm{\λ}}_1\left(t\right)+{\mathrm{\ϵ}}^2}---\left(7\right);$

$T\left(t\right)=\frac{{\[1-{\mathrm{\λ}}_2\left(t\right)/{\mathrm{\λ}}_1\left(t\right))}^2+{(1-{\mathrm{\λ}}_3(t)/{\mathrm{\λ}}_1(t\left)\right)}^2+{\left({\mathrm{\λ}}_2\right(t)/{\mathrm{\λ}}_1(t)-{\mathrm{\λ}}_3(t)/{\mathrm{\λ}}_1(t)\]}^2}{2{\[1-{\mathrm{\λ}}_2\left(t\right)/{\mathrm{\λ}}_1\left(t\right)+{\mathrm{\λ}}_3\left(t\right)/{\mathrm{\λ}}_1\left(t\right)\]}^2}---\left(8\right);$

式(7)和(8)中的e(t)为椭球率，ε为防止分母为零而引入的小数，其取值为0.0001，T(t)为极化率。

椭球率代表的是次轴与主轴之间的模长比值，地震信号线性极化程度越高，这个比值越趋近于零，线性极化程度越低，这个比值越接近于1。极化率则刚好相反，时窗内信号线性极化程度越高，极化率越趋近于1，线性极化程度越低，极化率越趋近于零。

自适应极化滤波同样基于协方差矩阵，利用Hilbert变换得到信号在t时刻的瞬时振幅、瞬时频率及瞬时相位信息，依据瞬时频率来给定时窗大小，并用这三个瞬时信息构造谐波来近似信号(Diallo,et al.,2006)。

其中，时间域自适应极化分析方法具体如(9)—(13)所示：

$T_0\left(t\right)=\frac{2\mathrm{\π}}{{\mathrm{\Ω}}_{xyz}\left(t\right)}=\frac{6\mathrm{\π}}{{\mathrm{\Ω}}_x\left(t\right)+{\mathrm{\Ω}}_y\left(t\right)+{\mathrm{\Ω}}_z\left(t\right)}---\left(9\right);$

$c_i\left(t\right)=s_i\left(t\right)+{\mathrm{js}}_i^h\left(t\right)---\left(10\right);$

${\stackrel{\‾}{s}}_i(t+\mathrm{\τ})\≈\left|c_i\left(t\right)\right|cos\[{\mathrm{\Ω}}_i\left(t\right)\mathrm{\τ}+\arg c_i\left(t\right)\]---\left(11\right);$

Ω_i(t)为s_i(t)的瞬时频率值，Ω_xyz(t)为三个分量的瞬时频率的平均值，T₀(t)为求得的最优滤波时窗大小；c_i(t)表示的是s_i(t)的解析信号，为地震数据s_i(t)的Hilbert变换；τ代表的是t附近的时间增量，arg表示相位角算子。将时间t附近[t-T₀/2,t+T₀/2]的信号用近似，则自适应协方差矩阵中各元素可以根据式(9)、(10)直接计算，从而避免了大量的求和运算。其中，表示复数的实部。

$\begin{array}{c}{I}_{ij}\left(t\right)=\frac{1}{T_0}\underset{-T_0/2}{\overset{T_0/2}{\∫}}\left({\stackrel{\‾}{s}}_i\right(t+\mathrm{\τ})-{\stackrel{\‾}{u}}_i(t\left)\right)\left({\stackrel{\‾}{s}}_j\right(t+\mathrm{\τ})-{\stackrel{\‾}{u}}_j(t\left)\right)d\mathrm{\τ}=\left|c_i\left(t\right)\right|\left|c_j\left(t\right)\right|\\ \×\{\sin c\[\frac{{\mathrm{\Ω}}_i\left(t\right)-{\mathrm{\Ω}}_j\left(t\right)}{2}{T}_0\left(t\right)\]\cos \[\arg c_i\left(t\right)-\arg c_j\left(t\right)\]+\sin c\[\frac{{\mathrm{\Ω}}_i\left(t\right)+{\mathrm{\Ω}}_j\left(t\right)}{2}{T}_0\left(t\right)\]\cos \[\arg c_i\left(t\right)+\arg c_j\left(t\right)\]\}-{\stackrel{\‾}{u}}_i\left(t\right){\stackrel{\‾}{u}}_j\left(t\right)\end{array}---\left(12\right);$

值得注意的是，此时的与M(t)已经有本质区别。M(t)需要统计一段时间[t-T₀/2,t+T₀/2]内的粒子极化情况，而对于的计算，时窗内的信号已经由t时刻的瞬时信息近似得到，因此只包含了t时刻的瞬时信息，得到的是每一个时间点的瞬时极化信息。并且自适应极化滤波的时间窗口大小T₀(t)随着各分量的平均瞬时频率Ω_xyz(t)的变化而变化，使得每个点的时间窗口总能调整到信号的一个子波周期大小，从而使得窗口最优化。

Diallo等提出在将自适应极化滤波延伸到小波域中，在小波域内计算统计地震数据中质点振动的极化分布。该方法同样基于协方差矩阵，并用一个近似方程来计算时窗内的协方差矩阵，时窗大小由时频点上的瞬时频率决定，具体如下：

首先，将一维信号s_k(t)与经过伸缩、平移之后的母小波g(t)做内积，得到小波系数WS_k(b,a)，如式(14)所示：

${\mathrm{WS}}_k(b,a)=<g,s_k\left(t\right)>={\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}\frac{1}{a}{g}^{*}\left(\frac{t-b}{a}\right)s_k\left(t\right)dt---\left(14\right);$

再由(15)式得到信号在时频域中的瞬时频率：

${\mathrm{\&Omega;}}_k(t,a)=\frac{\&part;}{\&part;t}\arg {\mathrm{WS}}_k(t,a)---\left(15\right);$

尺度a下的小波变换WS_k(t,a)在时间t附近的局部近似值为

这个值可以用小波系数的模与瞬时相位来表示：

${\hat{s}}_k(t+\mathrm{\&tau;},a)\&ap;\left|{\mathrm{WS}}_k(t,a)\right|cos\&lsqb;{\mathrm{\&Omega;}}_k(t,a)\mathrm{\&tau;}+\arg {\mathrm{WS}}_k(t,a)\&rsqb;---\left(16\right);$

则时频域中的协方差矩阵可改写为(17)-(18)式：

$M(t,a)=\left[\begin{array}{ccc}{I}_{xx}(t,a)& I_{xy}(t,a)& I_{xz}(t,a)\\ I_{yx}(t,a)& I_{yy}(t,a)& I_{yz}(t,a)\\ I_{zx}(t,a)& I_{zy}(t,a)& I_{zz}(t,a)\end{array}\right]---\left(17\right);$

$\begin{array}{c}{I}_{km}(t,a)=\frac{1}{T_{km}(t,a)}{\&Integral;}_{-T_{km}(t,a)/2}^{T_{km}(t,a)/2}\&lsqb;{\hat{s}}_k(t+\mathrm{\&tau;},a)-u_{km}\&rsqb;\&lsqb;{\hat{s}}_m(t+\mathrm{\&tau;},a)-u_{mk}\&rsqb;d\mathrm{\&tau;}\\ =\left|{\mathrm{WS}}_k(b,a)\right|\left|{\mathrm{WS}}_m(b,a)\right|\\ \&times;\{\sin c\&lsqb;\frac{{\mathrm{\&Omega;}}_k\left(t\right)-{\mathrm{\&Omega;}}_m\left(t\right)}{2}{T}_{km}\left(t\right)\&rsqb;\cos \cos \&lsqb;\arg {\mathrm{WS}}_k(b,a)-\arg {\mathrm{WS}}_m(b,a)\&rsqb;\\ +\sin c\&lsqb;\frac{{\mathrm{\&Omega;}}_k\left(t\right)+{\mathrm{\&Omega;}}_m\left(t\right)}{2}{T}_{km}\left(t\right)\&rsqb;\cos \cos \&lsqb;\arg {\mathrm{WS}}_k(b,a)+\arg {\mathrm{WS}}_m(b,a)\&rsqb;\}\\ -u_{km}{u}_{mk}\end{array}---\left(18\right);$

其中，T_km(t,a)为时频点[t,a]的时窗大小(式(19))，它的取值取决于三个分量在各个时频点上的瞬时频率大小。u_km(b,a)为均值，定义为(20)式：

$T_{km}(t,a)=\frac{6\mathrm{\&pi;}}{{\mathrm{\&Omega;}}_x(t,a)+{\mathrm{\&Omega;}}_y(t,a)+{\mathrm{\&Omega;}}_z(t,a)}---\left(19\right);$

同样地，求取协方差矩阵M(t,a)的特征值λ₁(t,a)＞λ₂(t,a)＞λ₃(t,a)与特征向量v₁(t,a)，v₂(t,a)，v₃(t,a)，则小波域中的椭球率可以表示为：

$e(t,a)=\frac{{\mathrm{\&lambda;}}_2(t,a)}{{\mathrm{\&lambda;}}_1(t,a)+{\mathrm{\&epsiv;}}^2}---\left(21\right).$

3、极化滤波应用于面波压制

时频域极化滤波将传统的极化滤波变换到时频域(如小波域、S域)进行极化分析，其优势在于能够对时间域重合而频率域不重合的面波与有效波进行分离，同时相比于视速度滤波，该方法对分道进行处理，能够避免假频面波的干扰。

同样以图1中的地震数据为例，时频域极化滤波对各道单独处理，此处以第85道(地震数据的第85列)为例作简要说明。

首先，将地震数据变换到小波域中，变换后的小波谱如图6中的(a)和(b)所示。从图6的(a)和(b)可以看出面波频率随着时间线性增加，与反射信号频率越来越接近，重叠区域不断增加，尤其是后两个反射信号，重叠更加严重。

然后，在小波域计算信号的椭球率，本例椭球率如图7所示，由于面波能量较强，所以在面波与反射波重叠的区域椭球率也呈现出高值的特点。根据椭球率设置阈值，压制掉椭球率e＞0.5的部分，滤波后的小波谱如图8中的(a)和(b)所示。可见滤波后小波谱的反射波区域已经残缺了部分信息，这是由于面波与反射波在小波域中重叠，导致重叠部分的椭球率偏高，因此重叠部分的反射波信息被滤除。

该技术针对单道进行处理，能够避免视速度引起的假频面波干扰。但是由于时频域仍然不能完全将面波与反射信号完全分开，因此无法通过时频域极化滤波方法将面波滤除。如上述例子所示，若对面波进行压制容易滤除反射波信息。

现有方法无论是FK滤波或者是时频域极化滤波，均无法对面波与反射波进行有效分离，究其原因在于这些方法使用的变换无法将面波与反射波完全分开，因此滤波无法彻底。

Shearlet变换(剪切波变换)，继承了Curvelet变换和Contourlet变换的优点，是一种新型多尺度几何分析工具。它通过对基函数进行缩放、平移和剪切等仿射变换来生成具有不同特征的Shearlet函数，对于包含C2奇异曲线或曲面的高维信号具有最优逼近特性。对于二维信号，它不仅能够检测到所有的奇异点，而且能够自适应跟踪奇异曲线的方向，并且随着尺度参数变化，能准确的描述函数的奇异性特征。与其它多尺度分析工具相比，Shearlet变换具有更简单的数学结构，对尺度方向表征更加细腻。并且对于一个N×N图像，经过Shearlet变换后的各尺度各方向上的系数仍然是N×N，从而有利于极化分析的进行。

发明内容

针对现有技术中存在的缺陷，本发明的目的在于提供一种基于Shearlet域极化滤波的面波压制方法，该方法能够在更高维的域(Shearlet域)中分离面波与反射波，并用极化分析来检测面波部分，对面波部分进行剔除，从而更加彻底的去除面波干扰。

为达到以上目的，本发明采取的技术方案是：

一种基于Shearlet域的极化滤波面波压制方法，其特征在于，包括如下步骤：

1)Shearlet变换：采用Shearlet变换将地震数据分解为各尺度各方向上的Shearlet系数，(Shearlet系数除了最粗尺度只有1个方向外，其余尺度i均有2^[j/2]+2个方向，因此一共可以得到个不同尺度不同方向的Shearlet系数；n为分解层数或尺度数，根据地震数据大小选定，n越大结果越准确，但计算速度越慢；j取2、3、4、……或n；i为具体的尺度，i取1、2、3、……或n)

2)极化分析：通过极化分析方法计算各尺度不同方向上的Shearlet系数的椭球率e，

3)面波压制：对各尺度不同方向上的Shearlet系数采用不同的椭球率阈值e₀进行滤波，以最大限度地保留反射波并最大限度地压制面波。

在上述基于Shearlet域的极化滤波面波压制方法中，还包括步骤4)：对面波压制后的Shearlet系数进行反变换，获得滤波结果。

在上述基于Shearlet域的极化滤波面波压制方法中，步骤3)中所述不同的椭球率阈值e₀为：

在反射波出现的方向上设定Shearlet系数的椭球率阈值e₀＝e₁，以最大限度地保留反射波；

在只含有少量或者不含有反射波的方向上设定Shearlet系数的椭球率阈值e₀＝e₂，以最大限度地压制面波；

且e₁>e₂,以在面波与反射波重叠的部分尽可能减少对反射波的损害。

在上述基于Shearlet域的极化滤波面波压制方法中，

所述反射波出现的方向为小于面波的临界角度的方向(即小角度方向)；

所述只含有少量或者不含有反射波的方向为大于面波的临界角度的方向(即大角度方向)。

在所述小角度方向即视速度较大的方向上，大多是反射波也有一些假频面波，为了最大限度保留反射波，对这部分的Shearlet系数设定的椭球率阈值较高；

在所述大角度方向即视速度较小的方向上，只含有少量或者不含有反射波，大多是面波，为了最大限度压制面波，对这部分的Shearlet系数设定的椭球率阈值较低。

在上述基于Shearlet域的极化滤波面波压制方法中，所述e₂为大角度方向上所设的椭球率阈值，该值的选取参照实际大角度方向上椭球率谱适当选取。

所述e₁为为小角度方向上所设阈值，该值的选取参照实际小角度方向上椭球率谱适当选取。

在上述基于Shearlet域的极化滤波面波压制方法中，所述滤波为滤除e＞e₀的部分，得到各尺度不同方向上滤波后的Shearlet系数。

在上述基于Shearlet域的极化滤波面波压制方法中，步骤2)所述极化分析方法为自适应极化分析方法。

在上述基于Shearlet域的极化滤波面波压制方法中，所述自适应极化分析方法为时间域自适应极化分析方法。

在上述基于Shearlet域的极化滤波面波压制方法中，所述时间域自适应极化分析方法具体可使用式(9)-(13)进行：

$T_0\left(t\right)=\frac{2\mathrm{\&pi;}}{{\mathrm{\&Omega;}}_{xyz}\left(t\right)}=\frac{6\mathrm{\&pi;}}{{\mathrm{\&Omega;}}_x\left(t\right)+{\mathrm{\&Omega;}}_y\left(t\right)+{\mathrm{\&Omega;}}_z\left(t\right)}---\left(9\right);$

$c_i\left(t\right)=s_i\left(t\right)+{\mathrm{js}}_i^h\left(t\right)---\left(10\right);$

${\stackrel{\&OverBar;}{s}}_i(t+\mathrm{\&tau;})\&ap;\left|c_i\left(t\right)\right|cos\&lsqb;{\mathrm{\&Omega;}}_i\left(t\right)\mathrm{\&tau;}+\arg c_i\left(t\right)\&rsqb;---\left(11\right);$

$\begin{array}{c}{I}_{ij}\left(t\right)=\frac{1}{T_0}\underset{-T_0/2}{\overset{T_0/2}{\&Integral;}}\left({\stackrel{\&OverBar;}{s}}_i\right(t+\mathrm{\&tau;})-{\stackrel{\&OverBar;}{u}}_i(t\left)\right)\left({\stackrel{\&OverBar;}{s}}_j\right(t+\mathrm{\&tau;})-{\stackrel{\&OverBar;}{u}}_j(t\left)\right)d\mathrm{\&tau;}=\left|c_i\left(t\right)\right|\left|c_j\left(t\right)\right|\\ \&times;\{\sin c\&lsqb;\frac{{\mathrm{\&Omega;}}_i\left(t\right)-{\mathrm{\&Omega;}}_j\left(t\right)}{2}{T}_0\left(t\right)\&rsqb;\cos \&lsqb;\arg c_i\left(t\right)-\arg c_j\left(t\right)\&rsqb;+\sin c\&lsqb;\frac{{\mathrm{\&Omega;}}_i\left(t\right)+{\mathrm{\&Omega;}}_j\left(t\right)}{2}{T}_0\left(t\right)\&rsqb;\cos \&lsqb;\arg c_i\left(t\right)+\arg c_j\left(t\right)\&rsqb;\}-{\stackrel{\&OverBar;}{u}}_i\left(t\right){\stackrel{\&OverBar;}{u}}_j\left(t\right)\end{array}---\left(12\right);$

式(9)—(13)中，

Ω_i(t)为s_i(t)的瞬时频率值，

Ω_xyz(t)为三个分量的瞬时频率的平均值，

T₀(t)为求得的最优滤波时窗大小；

c_i(t)表示的是s_i(t)的解析信号，

为地震数据s_i(t)的Hilbert变换；

τ代表的是t附近的时间增量，

arg表示相位角算子，

表示复数的实部，

将时间t附近[t-T₀/2,t+T₀/2]的信号用近似，则自适应协方差矩阵中各元素可以根据式(9)、(10)直接计算，从而避免了大量的求和运算。

在上述基于Shearlet域自适应极化滤波的面波压制方法中，面波的临界角度θ₁的计算方法如式(22)所示：

${\mathrm{\&theta;}}_1=acr\cot \left(\frac{v_{0}dt}{dx}\right)---\left(22\right);$

其中，v₀为面波的最大视速度，dt为采样时间间隔，dx为空间采样间隔。

在上述基于Shearlet域自适应极化滤波的面波压制方法中，步骤1)使用如下式(23)进行：

SH_ψf(j,k,t,x)＝〈f,ψ_j,k,t,x〉(23)，

其中，SH(j,k,t,x)代表不同尺度不同方向上的Shearlet系数，j代表尺度，k代表方向，t代表时间，x代表偏移距，f代表频率。

本发明的有益效果如下：

本发明所提供的一种基于Shearlet域的极化滤波面波压制方法，在shearlet域内进行极化分析，可以使得不同视速度的面波与反射波分布在不同的方向下，并且在各方向下进行极化分析，可以避免假频面波的干扰。本发明的方法可以对不同方向的shearlet系数采用不同的椭球率阈值，从而在反射波的方向上，椭球率阈值偏大，尽可能的保留反射波，而在面波的方向上椭球率阈值偏小，尽可能的滤除面波。

附图说明

本发明有如下附图：

图1为原始地震数据；横坐标为偏移距(m)，纵坐标为时间(s)；

图2为图1地震数据的FK谱；横坐标为波数(m-¹)，纵坐标为频率(Hz)；图中色柱的颜色代表FK谱值；

图3为切除面波后的FK谱；横坐标为波数(m-¹)，纵坐标为频率(Hz)；图中色柱的颜色代表FK谱值；

图4为滤波后的地震记录；横坐标为偏移距(m)，纵坐标为时间(s)；

图5为三分量极化示意图；

图6为第85道小波谱；横坐标为时间(s)，纵坐标为频率(Hz)；(a)X分量；(b)Z分量；每张图中上部的三个箭头尾处所指为反射波；下部的一个箭头尾处所指为面波；图中色柱的颜色代表小波频谱值；

图7为第85道椭球率；横坐标为时间(s)，纵坐标为频率(Hz)；图中色柱的颜色代表椭球率的值；

图8为第85道滤波后小波谱；横坐标为时间(s)，纵坐标为频率(Hz)；(a)X分量；(b)Z分量；图中色柱的颜色代表小波频谱值；

图9为尺度4大角度方向上Shearlet极化滤波结果，(a)为滤波前的子剖面，(b)为滤波前的椭球率，(c)为滤波后的子剖面；横坐标为偏移距(m)，纵坐标为时间(s)；图中色柱的颜色代表椭球率的值；

图10为尺度4小角度方向上Shearlet极化滤波结果，(a)为滤波前的子剖面，(b)滤波后的子剖面，(c)为滤波前的椭球率；横坐标为偏移距(m)，纵坐标为时间(s)；图中色柱的颜色代表椭球率的值；

图11为Shearlet极化滤波剖面，(a)(b)分别为滤波前的X/Z分量剖面，(c)(d)分别为Shearlet极化滤波后的剖面；横坐标为偏移距(m)，纵坐标为时间(s)；

图12为Shearlet滤波剖面FK谱，(a)X分量，(b)Z分量；横坐标为波数(m^-1)，纵坐标为频率(Hz)；图中色柱的颜色代表Shearlet系数FK谱值；

图13为Shearlet极化滤波第85道时频谱对比，(a)(b)分别为原始不含面波模拟第85道X、Z分量时频谱；(c)(d)分别为Shearlet极化滤波后第85道X、Z分量时频谱；横坐标为时间(s)，纵坐标为频率(Hz)；图中色柱的颜色代表椭球率的值；

图14为第85道数据分别采用小波域极化滤波与Shearlet域自适应极化滤波前后数据对比，横坐标为时间(s)，纵坐标为振幅。

具体实施方式

以下结合附图对本发明作进一步详细说明。

实施例1

本发明所述的一种基于Shearlet域的极化滤波面波压制方法，具体操作步骤如下：

1)Shearlet变换：将图1所示的地震数据分解为5个尺度(即n＝5)，采用Shearlet变换(式(23))共得到49个各尺度各方向上的Shearlet系数，即49个副子图像(m＝49)。

2)极化分析：通过时间域自适应极化分析方法(式(9)至(13))，计算各尺度不同方向上的Shearlet系数的椭球率e。

3)面波压制：对各尺度不同方向的Shearlet系数采用不同的椭球率阈值e₀进行滤波，滤除e＞e₀的部分，得到各尺度不同方向上滤波后的Shearlet系数，以最大限度地保留反射波并最大限度地压制面波。

所述不同的椭球率阈值e₀为：

在小于面波的临界角度的方向(即小角度方向)设定Shearlet系数的椭球率阈值e₀＝e₁＝0.6，以最大限度地保留反射波；

在大于面波的临界角度的方向(即大角度方向)设定Shearlet系数的椭球率阈值e₀＝e₂＝0.4，以最大限度地压制面波；

且e₁>e₂。

所述e₂比面波与反射波重叠的部分中反射波的椭球率最高值稍大，如在图10(c)中的面波与反射波重叠的部分中反射波的椭球率最高值为0.35，因此选取e₂为0.4。

面波的临界角度θ₁的计算方法如式(22)所示，具体如下：

面波的最大视速度v₀为1500m/s，已知采样时间间隔dt为0.001s，空间采样间隔dx为3m，代入式(22)，得出面波的临界角度θ₁为63.4°

据此判断反射波出现在小角度方向上，即分布于第29-48系数上，使用一个较高的椭球率阈值即0.6，滤除e＞0.6的部分。

在面波与反射波重叠较少的方向上即第1-28以及第49个系数上使用一个较低的椭球率阈值即0.4，滤除e＞0.4的部分。

4)对面波压制后的Shearlet系数进行反变换，获得滤波结果。

以第4尺度斜向为例，如图9中的(a)所示，面波与反射波出现在同一个方向中，但是二者几乎不重叠。计算得到该方向椭球率如图9中的(b)所示，在面波的区域，椭球率呈现高值，而反射波区域为低值。该方向经本发明方法滤波后的结果如图9中的(c)所示，低速的面波被压制的非常彻底且反射波信号基本没有损失。

而对于小角度方向，反射波与假频面波存在重叠，我们采用的方法是设置一个较高的椭球率阈值0.6，压制e＞0.6的部分。图10中的(a)为第4尺度小角度方向的Shearlet系数，可见面波与反射波存在重叠区域也有互相分离区域。本发明通过极化分析得到该方向椭球率如图10中的(c)所示，在面波域反射波未重叠的部分，面波椭球率呈现高值，反射波椭球率呈现出低值，而在二者重叠的部分如图中红色圆形区域，由于反射波能量仍然占主导，因此椭球率仍然比纯面波区域低，从而可以通过设置高椭球率阈值使得该圆形区域的反射波得到保留。对此方向的系数采用本发明方法的滤波结果如图10中的(b)所示，反射信号基本没有损失，而纯面波区域的面波也全部被滤除。

对处理完之后的各方向Shearlet系数进行Shearlet反变换，得到的最终的滤波结果如图11所示，可见Shearlet域极化滤波剖面不含有面波的假频部分，面波被压制得较为彻底，且反射同相轴缺失很少。滤波后地震剖面的FK谱如图12所示，经Shearlet极化滤波后的FK谱没有出现如图3所示的假频信息残留，与反射波重叠的假频面波被压制得较为彻底。同样抽取Shearlet极化滤波后的第85道信号进行小波变换，并与原始不含面波的反射记录的小波谱对比如图13所示，可以发现二者只有非常细微的差别，反射信号基本上都被保留下来。同样对Shearlet域极化滤波前后的第85道数据进行对比如图14所示，可见滤波后的保留的信号与原始反射信号非常吻合，只有少量的精度损失。这说明了本发明方法能够在彻底压制面波的同时不会对反射波造成损害，且能够抗假频干扰。

本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (8)

1.一种基于Shearlet域的极化滤波面波压制方法，其特征在于，包括如下步骤：

1)Shearlet变换：采用Shearlet变换将地震数据分解为各尺度各方向上的Shearlet系数，

2)极化分析：通过极化分析方法计算各尺度不同方向上的Shearlet系数的椭球率e，

3)面波压制：对各尺度不同方向上的Shearlet系数采用不同的椭球率阈值e₀进行滤波，以最大限度地保留反射波并最大限度地压制面波。

2.如权利要求1所述的基于Shearlet域的极化滤波面波压制方法，其特征在于：还包括步骤4)：

对面波压制后的Shearlet系数进行反变换，获得滤波结果。

3.如权利要求1所述的基于Shearlet域的极化滤波面波压制方法，其特征在于：

步骤3)中所述不同的椭球率阈值e₀为：

在反射波出现的方向上设定Shearlet系数的椭球率阈值e₀＝e₁，以最大限度地保留反射波；

在只含有少量或者不含有反射波的方向上设定Shearlet系数的椭球率阈值e₀＝e₂，以最大限度地压制面波；

且e₁>e₂,以在面波与反射波重叠的部分尽可能减少对反射波的损害。

4.如权利要求3所述的基于Shearlet域的极化滤波面波压制方法，其特征在于：

所述反射波出现的方向为小于面波的临界角度的方向；

所述只含有少量或者不含有反射波的方向为大于面波的临界角度的方向。

5.如权利要求1所述的基于Shearlet域的极化滤波面波压制方法，其特征在于：

所述滤波为滤除e＞e₀的部分，得到各尺度不同方向上滤波后的Shearlet系数。

6.如权利要求1所述的基于Shearlet域的极化滤波面波压制方法，其特征在于：步骤2)所述极化分析方法为时间域自适应极化分析方法。

7.如权利要求4所述的基于Shearlet域极化滤波面波压制方法，其特征在于：

面波的临界角度θ₁的计算方法如式(22)所示：

${\mathrm{\&theta;}}_1=acr\cot \left(\frac{v_{0}dt}{dx}\right)---\left(22\right);$

其中，v₀为面波的最大视速度，dt为采样时间间隔，dx为空间采样间隔。

8.如权利要求1所述的基于Shearlet域的极化滤波面波压制方法，其特征在于：步骤1)使用如下式(23)进行：

SH_ψf(j,k,t,x)＝<f,ψ_j,k,t,x> (23)，

其中，SH(j,k,t,x)代表不同尺度不同方向上的Shearlet系数，j代表尺度，k代表方向，t代表时间，x代表偏移距，f代表频率。