CN103091714B - 一种自适应面波衰减方法 - Google Patents

Abstract

本发明提供了一种自适应面波衰减方法，属于石油地球物理勘探领域。所述方法首先对地震记录进行低通滤波，获得面波能量占优的记录；然后利用Hilbert变换将记录转换成解析信号、利用包络振幅差异原理寻找每个样点处最优的视速度；再在该样点处沿最优视速度选取一时空窗，利用奇异值分解选取前1～2个最大奇异值重构时窗信号，只保留当前样点处的振幅值，逐点完成后得到面波模型；然后将面波模型与原始记录进行匹配滤波，进一步修正面波模型的波形差异；最终将面波模型从原始记录中减掉得到最终的处理结果。本发明方法实现了每一点面波方向的精确定位、在一定程度上克服了地震记录中面波纵横向振幅变化，具有自适应性。

Description

一种自适应面波衰减方法

技术领域

本发明属于石油地球物理勘探领域，具体涉及一种自适应面波衰减方法。

背景技术

在地震勘探中，面波通常是很强的干扰波，它存在于地震记录中，会大大降低资料的信噪比，特别是复杂区地震记录，面波的表现形式比较复杂，对于资料的后续处理非常不利。因此，对复杂区地震记录面波剔除与压制方法的研究是非常必要的。

在野外采集时，常用组合法压制面波，它利用面波的方向特性，将几个检波器按一定方式放置在一起，使面波有一定时差，而反射波时差较小，按一定权重相加后，可将面波削弱，而有效波增强。但是组合法有一个明显的缺点，即压制高频，降低分辨率，并且只有在采集时使用，不利于信号数字处理手段的应用。

在室内资料处理时，目前常用的方法有：F-K滤波，τ-p滤波，中值滤波，K-L变换，分频带处理等。它们主要是依靠面波在地震记录中的低频、低视速度和线性特征。其中F-K滤波、τ-p滤波是将信号变换到F-K域、τ-p域，根据面波与有效反射波在传播方向上的差别滤波；中值滤波、K-L变换利用面波在记录中的线性特征以及反射波的非线性特征滤波；分频带处理利用面波与反射波在频率上的差别，与其它方法相结合剔除面波。

但是，这些处理方法存在以下不足：(1)对于时间上和面波不重合的反射波部分会有影响，例如，τ-p滤波、中值滤波、K-L变换仅考虑面波的方向特性，当地表、地下构造复杂时，面波和部分有效反射波的视速度相差不大，因此，剔除与压制面波的同时，势必压制了这部分有效反射波，而分频带与其它方法相结合，仅考虑频带特性，对于时间上与面波不重合的有效信号的低频带范围，也会有影响；(2)不适合复杂区的地震记录面波剔除和压制。对于均匀地表介质，面波在地震记录中的表现形式与地表及观测系统有关。当地表水平或单斜时，二维地震勘探记录中面波呈线性分布，三维地震勘探的二维剖面中面波呈双曲线分布，地表起伏及弯线勘探时，面波在记录中的分布与地表几何形态有关，其特征更加复杂。因此，上述方法中根据其线性特征剔除和压制面波会产生误差。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种自适应面波衰减方法，对地震记录中的面波线性化后利用面波的极化特征，提取更为精确的面波分量，然后将其从原始资料中减掉。

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

一种自适应面波衰减方法，所述方法根据面波在复数域中的极性特征，利用奇异值分解将面波模型从原始地震记录中提取出来，然后将所述面波模型从原始地震记录中减去就得到了去噪后的地震记录。

所述方法包括以下步骤：

(1)自动方向追踪；

(2)局部时空窗数据选择；

(3)奇异值分解；

(4)自适应衰减。

所述步骤(1)具体如下：首先根据面波的优势频带对原始地震记录进行低通滤波，获得面波能量占优的地震记录；然后利用Hilbert变换将所述面波能量占优的地震记录转换成复数域中的解析信号；最后利用包络振幅差异原理寻找每个计算样点处的最优视速度；

所述利用Hilbert变换将所述面波能量占优的地震记录转换成复数域中的解析信号是这样实现的：

在某一地震道的某个计算样点处取n个视速度，并取该地震道前后各m道，经Hilbert变换成解析信号，具体如下；

设该计算样点处的径向分量、切向分量和垂直分量分别为u_r(t)、v_r(t)和w_r(t)，利用Hilbert变换将各分量转换成解析信号：

u(t)＝u_r(t)+iH(u_r(t))

v(t)＝v_r(t)+iH(v_r(t))

w(t)＝w_r(t)+iH(w_r(t))

式中H表示Hilbert变换。

所述利用包络振幅差异原理寻找每个计算样点处的最优视速度是这样实现的：

首先，所述解析信号共2m+1道，沿其中一个视速度方向，在彼此相邻的道上选取包络值，再将两者的包络值相减得到相邻道之间的振幅包络差δ_amp，然后将2m个δ_amp求和得到Amp：

${\mathrm{\δ}}_{\mathrm{amp}}^i=h_i-h_{i-1}$

$\mathrm{Amp}=\underset{i=1}{\overset{2m}{\mathrm{\Σ}}}{\mathrm{\δ}}_{\mathrm{amp}}^i$

上式中，h_i、h_i-1分别表示表示相邻道i和i-1的包络值；

然后同理求得所有n个视速度对应的Amp值，最小的Amp值对应的视速度即为该计算样点处的最优视速度，这样就完成了自动方向追踪。

所述步骤(2)具体如下：以所述计算样点为中心点选取局部时空窗，其中时窗为p个样点，道窗取2m+1道，在该时空窗内，用沿步骤(1)获得的最优视速度的方向截取所述解析信号的数据，将其放到矩阵中，得到的数据矩阵为(2m+1)×p；然后每一道均与中间道(即第m+1道)做互相关，求取剩余相对时差δt，利用所述剩余相对时差δt对各道进行静态时移，进一步校正相邻道之间的时差。

所述步骤(3)具体如下：对步骤(2)中得到的数据矩阵进行奇异值分解，得到该数据矩阵的所有奇异值，将所有奇异值按从大到小排序，然后取排在前面的1到2个奇异值重构该数据矩阵中相关程度最强的水平同相轴的能量，即重构了面波能量，在重构后的数据中取中心点处的振幅值即得到了所述计算样点处的噪音模型。

所述步骤(4)具体如下：逐点计算得到所有计算样点处的噪音模型，即对所有计算样点重复步骤(1)-(3)，得到所有计算样点处的噪音模型，所有计算样点处的噪音模型构成了面波模型model(t)；然后对得到的面波模型model(t)进行分时窗匹配滤波，最后将经过分时窗匹配滤波后的面波模型从原始地震记录中减掉，得到去噪后的有效地震记录。

所述分时窗匹配滤波是这样实现的：

$e=\underset{t}{\mathrm{\Σ}}\underset{i}{\mathrm{\Σ}}{({\mathrm{win}}_i\left(t\right)\·d\left(t\right)-{\mathrm{\ω}}_i\left(t\right)*{\mathrm{win}}_i\left(t\right)\·\mathrm{mode}l\left(t\right))}^2$

式中t为时间，i为时窗序号，win_i(t)为时窗加权函数，ω_i(t)为第i个时窗的匹配系数，*表示褶积，“.”表示乘积。

与现有技术相比，本发明的有益效果是：

1)解析域完成过程，充分利用了面波在解析域中强能量的极化特征；

2)逐点搜索最优视速度方向及确定面波模型，实现了每一点面波方向的精确定位；

3)在逐点构建的小时空窗内提取面波模型，在一定程度上克服了地震记录中面波纵横向振幅变化，具有自适应性。

附图说明

图1是Hilbert变换示意图。

图2是本发明自适应面波衰减方法的步骤框图。

图3是本发明自适应面波衰减方法的步骤(2)中局部时空窗数据选择的示意图。

图4是本发明实施例中的处理效果图，其中图4-1是原始单炮(即地震记录)，图4-2是利用本发明方法开发的模块处理的结果，图4-3是得到的噪音模型。

具体实施方式

下面结合附图对本发明作进一步详细描述：

本发明方法的原理如下：

(1)面波的极性特征

在多波多分量地震资料中，有多种极性化分析方法，如协方差矩阵和奇异值分解法。Vidale(1986)提出了利用复数地震道方法研究多分量资料的极性化特征。其极性化分析是在Hilbert复数域中进行的，如图1所示。

设径向分量、切向分量和垂直分量分别为u_r(t)、v_r(t)和w_r(t)，将各分量转换成解析信号：

u(t)＝u_r(t)+iH(u_r(t))

v(t)＝v_r(t)+iH(v_r(t))

w(t)＝w_r(t)+iH(w_r(t))

式中H表示Hilbert变换。

对于Rayleigh波而言，变换到Hilbert域中其极性化轨迹为椭圆。实际地震记录接收到地表噪音不仅包括Rayleigb波，还包括Love波和折射波。多种波的存在降低了极性分析的可靠度，极化轨迹不再是椭圆。但是在Hilbert域中面波依然表现为强能量噪音，这就提供了利用奇异值分解(SVD)提取噪音模型的物理基础。

(2)奇异值分解

令A∈R^M×N，则存在正交矩阵U∈R^M×M、V∈R^N×N和对角矩阵∑使得

A＝U∑V^T

式中

$\mathrm{\Σ}=\left[\begin{array}{cc}{\mathrm{\Σ}}_1& 0\\ 0& 0\end{array}\right]$

且∑₁＝diag(δ₁，δ₂，......，δ_r)，其对角元素按照顺序δ₁≥δ₂≥......≥δ_r＞0，r＝rank(A)排列。这就是矩阵的奇异值分解定理，称δ_i(i＝1，2，…，r-1，r)为矩阵A的奇异值，而A＝U∑V^T称为矩阵A的奇异值分解式。

当N道地震数据为线性无关时，它的秩r＝N。此时所有的δ_i均不为零，因此要完整地重构X就需要把所有的特征图像进行加权求和，即

$x=\underset{I=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\δ}}_i{u}_i{v}_i^T$

如果X的各道地震数据相似，即N道记录全部线性相关时，X的秩r＝1，即只用一个特征图像加权，即就可以完全重构原数据体X。一般情况存在一个p＜r，按公式对特征图像进行加权求和来重构地震记录X。

如果仅用前p个特征图像来重构X，其重构误差为

$\mathrm{\ϵ}=\underset{i=p+1}{\overset{r}{\mathrm{\Σ}}}{\mathrm{\δ}}_i^2$

因此，重构地震记录X所需要的特征图像的个数依赖于X的道与道之间的线性相关性。相关程度越高，所需要的特征图像的个数就越少。且最大奇异值对应的子矩阵横向相关性最大，较小奇异值对应的子矩阵横向相关性较小。

如图2所示，本发明方法包括以下步骤：

(1)自动方向追踪

首先对原始地震记录进行低通滤波(例如0-12Hz，这是根据资料的特点人为选取的，也可以为其它频率范围。)，得到面波能量占优的地震记录(即图2中的面波优势频带)。在某一道的某个计算样点处取n个视速度(所有的视速度都小于面波最大速度)，并取该道前后各m道，经Hilbert变换成解析信号(利用原理(1)中的三个公式)，共2m+1道，沿一个视速度方向，在彼此相邻的道上选取包络值，然后将两者的包络值相减得到相邻道之间的振幅包络差δ_amp，然后将2m个δ_amp求和得到Amp，

${\mathrm{\δ}}_{\mathrm{amp}}^i=h_i-h_{i-1}$

$\mathrm{Amp}=\underset{i=1}{\overset{2m}{\mathrm{\Σ}}}{\mathrm{\δ}}_{\mathrm{amp}}^i$

式中h_i、h_i-1分别表示表示相邻道i和i-1的包络值。

同理求得所有n个视速度对应的Amp值，最小的Amp值对应的视速度即为该计算样点处最优视速度方向，这样就完成了自动求取最优速度，也就是图2中的视速度扫描。

(2)局部时空窗数据选择

如图3所示，以计算样点为中心点选取局部时空窗，其中时窗为p个样点，道窗取2m+1道，即2m+1道，每道p个样点，在以计算样点为中心的邻区内，用沿步骤(1)得到的最优视速度方向截取解析信号数据，将其放到矩阵中，得到数据矩阵为(2m+1)×p。由于截取的数据有可能存在不平的现象，所有每一道均与中间道(第m+1道)做互相关，求取剩余相对时差δt，利用所述剩余相对时差δt对各道进行静态时移，进一步校正相邻道之间的时差，利用互相关将其彻底拉平。

(3)奇异值分解

将步骤(2)得到的的数据矩阵进行奇异值分解，然后取前1到2个奇异值重构该数据中相关程度最强(奇异值分解中头两个最大的奇异值对应的信号，相关程度最强，这是奇异值分解的原理)的水平同相轴的能量，即重构面波能量，在重构的数据中仅保留计算样点(即中心点)处的振幅值，作为计算样点处的噪音模型。具体来说，经过Hilbert变换后，信号变成了复数信号，它的实部就是实际数据，而虚部数据是Hilbert变换后得到的数据，SVD重构后的信号还是复数的，因此该信号的实部就是振幅值，也就得到了该计算样点处的噪音模型。

(4)自适应衰减

逐点计算得到所有计算样点处的噪音模型，即对所有计算样点重复步骤(1)-(3)，得到所有样点处的噪音模型，即面波模型model(t)，逐点计算可以避免信号的纵横向变化。由于得到面波模型与原始记录中的面波存在一定的差异，需要进行分时窗匹配滤波：

$e=\underset{t}{\mathrm{\Σ}}\underset{i}{\mathrm{\Σ}}{({\mathrm{win}}_i\left(t\right)\·d\left(t\right)-{\mathrm{\ω}}_i\left(t\right)*{\mathrm{win}}_i\left(t\right)\·\mathrm{mode}l\left(t\right))}^2$

式中t为时间，i为时窗序号，win_i(t)为时窗加权函数，ω_i(t)为第i个时窗的匹配系数，*表示褶积，“.”表示乘积。

面波模型model(t)经过匹配滤波以后再从原始记录中减掉，得到去噪后的有效记录。

根据本发明方法的步骤，在Cygwin下利用fortran语言进行了实现(就是图2中的自主开发模块)，并得到了如图4所示的处理结果，其中图4-1是原始单炮(即地震记录)，图4-2是利用本发明方法开发的模块处理的结果，图4-3是得到的噪音模型，从图4可以看出，利用本方法的面波模型克服了面波频散的特点，构建的面波模型较其他视速度滤波技术更加接近实际资料中的面波形态。

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

Claims (6)

1.一种自适应面波衰减方法，其特征在于：所述方法根据面波在复数域中的极性特征，利用奇异值分解将面波模型从原始地震记录中提取出来，然后将所述面波模型从原始地震记录中减去就得到了去噪后的地震记录；

所述方法包括以下步骤：

(1)自动方向追踪；

(2)局部时空窗数据选择；

(3)奇异值分解；

(4)自适应衰减；

所述步骤(1)具体如下：首先根据面波的优势频带对原始地震记录进行低通滤波，获得面波能量占优的地震记录；然后利用Hilbert变换将所述面波能量占优的地震记录转换成复数域中的解析信号；最后利用包络振幅差异原理寻找每个计算样点处的最优视速度；

所述利用包络振幅差异原理寻找每个计算样点处的最优视速度是这样实现的：

首先，所述解析信号共2m+1道，沿其中一个视速度方向，在彼此相邻的道上选取包络值，再将两者的包络值相减得到相邻道之间的振幅包络差δ_amp，然后将2m个δ_amp求和得到Amp：

上式中，h_i、h_i-1分别表示相邻道i和i-1的包络值；

然后同理求得所有n个视速度对应的Amp值，最小的Amp值对应的视速度即为该计算样点处的最优视速度，这样就完成了自动方向追踪。

2.根据权利要求1所述的自适应面波衰减方法，其特征在于：所述利用Hilbert变换将所述面波能量占优的地震记录转换成复数域中的解析信号是这样实现的：

在某一地震道的某个计算样点处取n个视速度，并取该地震道前后各m道，经Hilbert变换成解析信号，具体如下；

设该计算样点处的径向分量、切向分量和垂直分量分别为u_r(t)、v_r(t)和w_r(t)，利用Hilbert变换将各分量转换成解析信号：

u(t)＝u_r(t)+iH(u_r(t))

v(t)＝v_r(t)+iH(v_r(t))

w(t)＝w_r(t)+iH(w_r(t))

式中H表示Hilbert变换。

3.根据权利要求2所述的自适应面波衰减方法，其特征在于：所述步骤(2)具体如下：以所述计算样点为中心点选取局部时空窗，其中时窗为p个样点，道窗取2m+1道，在该时空窗内，用沿步骤(1)获得的最优视速度的方向截取所述解析信号的数据，将其放到矩阵中，得到的数据矩阵为(2m+1)×p；然后每一道均与中间道做互相关，求取剩余相对时差δt，利用所述剩余相对时差δt对各道进行静态时移，进一步校正相邻道之间的时差。

4.根据权利要求3所述的自适应面波衰减方法，其特征在于：所述步骤(3)具体如下：对步骤(2)中得到的数据矩阵进行奇异值分解,得到该数据矩阵的所有奇异值，将所有奇异值按从大到小排序，然后取排在前面的1到2个奇异值重构该数据矩阵中相关程度最强的水平同相轴的能量,即重构了面波能量,在重构后的数据中取中心点处的振幅值即得到了所述计算样点处的噪音模型。

5.根据权利要求4所述的自适应面波衰减方法，其特征在于：所述步骤(4)具体如下：逐点计算得到所有计算样点处的噪音模型，即对所有计算样点重复步骤(1)－(3)，得到所有计算样点处的噪音模型，所有计算样点处的噪音模型构成了面波模型model(t)；然后对得到的面波模型model(t)进行分时窗匹配滤波，最后将经过分时窗匹配滤波后的面波模型从原始地震记录中减掉，得到去噪后的有效地震记录。

6.根据权利要求5所述的自适应面波衰减方法，其特征在于：所述分时窗匹配滤波是这样实现的：

式中t为时间，i为时窗序号，win_i(t)为时窗加权函数，ω_i(t)为第i个时窗的匹配系数，*表示褶积，“.”表示乘积。