在地震资料处理中,高度保持地震波的真振幅特征,尽量提高地震记录的信噪比和分辨率,即称为“三高”处理,这一直是地震资料处理人员追求的目标。因为“三高”处理的质量直接影响到岩性参数提取以及地震勘探的精度和效果。
10.3.1 真振幅恢复
保持地震波的真振幅特征(简称保幅处理),从广义讲应包含两大方面内容:即真振幅恢复(或称振幅补偿)和其他各项处理中的振幅保持问题。本节主要讨论真振幅恢复的方法,而对其他各项处理中凡要影响到振幅特征的处理方法,则要采取相应的措施,尽可能的使振幅的相对关系保持不变。
地震记录经增益恢复处理后,其振幅特征已与地表检波器所接收到的地震波振幅特征一致。这种振幅仍不称为真振幅,我们所谓的真振幅是指由地层波阻差而产生的反射波振幅,即能反映地层岩性变化的振幅。在地表所接收到的振幅除有地层波阻抗的变化因素外,还有球面扩散因素以及非弹性衰减的因素,因此需要消除球面扩散和非弹性衰减的影响,恢复地震波的真振幅特征。
球面扩散是当波离开震源传播时由于波前扩展造成的振幅衰减。这样的振幅衰减(A)与传播距离r成反比
勘查技术工程学
其中v是界面上覆介质的平均速度;t是反射的记录时间。对球面扩散作校正需要用时变函数vt乘以数据。
非弹性衰减是弹性波能量在岩石中传播时,由于内摩擦而耗散为热被地层吸收的结果。原理部分已说明这种衰减是频率和传播距离的指数形式的函数
勘查技术工程学
其中α为非弹性衰减系数(吸收系数)
勘查技术工程学
所以,用函数eαvt乘以数据就可校正非弹性衰减。至此,真振幅恢复处理完成。
系数α可从增益恢复及球面扩散校正后的振幅-时间函数来测定。为了得到α的较好统计估计,要用一组地震道测定能量来求得衰减曲线。
还有另一种真振幅恢复的方法,这时不需要速度信息。在增益恢复之后,假设振幅衰减是指数函数。因此,按照最小平方法,用指数函数拟合增益校正后的记录,就得到真振幅校正函数(即包括球面扩散和非弹性衰减校正两者)。
前述已知,波前发散因子K为
勘查技术工程学
式中r和t分别为波的传播距离和传播时间,C和a为与地层速度有关的常数。
吸收衰减因子是
勘查技术工程学
式中α是吸收系数;b是待定的常数。波前发散和吸收衰减总的影响是
勘查技术工程学
求取a和b的方法如下。
从地震记录上读取反射波的振幅极值(波峰或波谷),以(10.3-5)为回归方程,得
勘查技术工程学
式中:ut=lnAi-lnti;Ai,ti为振幅极值及其对应的时间;N为振幅极值的点数。校正函数是a-1tebt。
为了获得有代表性的真振幅恢复参数,所选的地震道应是没有多次波及有较高的信噪比。对地质条件稳定地区,一组参数就可代表全区。在工区内地质条件有较大变化时,这些参数要重新计算。
10.3.2 提高信噪比的数字滤波处理
在地震勘探中,用于解决地质任务的地震波称为有效波,而其他波统称为干扰波。压制干扰,提高信噪比是一项贯穿地震勘探全过程的任务。除在野外数据采集中采取相应措施压制干扰外,在地震资料数字处理中数字滤波也是一项非常重要的提高信噪比的措施。
数字滤波方法是利用有效波和干扰之间频率和视速度方面的差异来压制干扰的,分别称为频率滤波和视速度滤波。又因频率滤波只需对单道数据进行运算,故称为一维频率滤波。实现视速度滤波需同时处理多道数据,故称为二维视速度滤波。本节主要介绍这两种滤波方法。
10.3.2.1 一维频率滤波
所谓一维数字滤波是指用计算机实现对单变量信号的滤波,该单变量可以是时间或频率,也可以是空间或波数。以时间或频率为例讨论一维数字滤波,其他原理相同。
(1)一维数字滤波原理
设地震记录x(t)是由有效波s(t)和干扰波n(t)组成,即
勘查技术工程学
其频谱为
勘查技术工程学
式中:X(f)为 x(t)的频谱;S(f)、N(f)分别为 s(t)、n(t)的频谱。如果 X(f)的振幅谱|X(f)|可用图10-6表示。说明有效波的振幅谱|S(f)|处在低频段,而干扰波的振幅谱处于高频段。
图10-6 有效波和干扰波频谱分布示意图
若设计一频率域函数 H(f)的振幅谱为|H(f)|,
勘查技术工程学
其图形为图10-7(a)所示。
令
勘查技术工程学
及
勘查技术工程学
在时间域有(利用傅里叶变换的褶积定理)
勘查技术工程学
称 H(f)为一维滤波器频率响应,(10.3-9)式为频率域滤波方程,h(t)为 H(f)的时间域函数,称为一维滤波器滤波因子(图10-7(b))。(10.3-11)为时间域滤波方程,y(t)和 Y(f)分别为滤波后仅存在有效波的地震记录及频谱,φx(f)、φy(f)、φh(f)分别为滤波前、滤波后地震记录及滤波器的相位谱,以上滤波主要是利用了有效波和干扰波的频率差异消除干扰波,故也称为频率滤波。
图10-7 滤波频率响应及滤波因子
以上所述的滤波器称为理想低通滤波,根据有效波和干扰波的频段分布不同,还可将滤波器分为理想带通滤波器、理想高通滤波器等。所谓理想是指滤波器的频率响应是一个矩形门,门内的有效波无畸变地通过,称为通频带,而门外的干扰波全部消除。在数字滤波中这一点实际是做不到的。因为数字滤波时所能处理的滤波因子只能是有限长,而由间断函数组成的理想滤波器的滤波因子是无限长的。实际应用中只能截断为有限长,截断后就会出现截断效应,即截断后的滤波因子所对应的频率响应不再是一个理想的矩形门,而是一条接近矩形门,但有振幅波动的曲线,这种现象称为吉普斯现象。
由于频率响应曲线在通频带内是波动的曲线,滤波后有效波必定会发生畸变。另外,在通频带外也是波动的曲线,必定不能有效地压制干扰。为了避免吉普斯现象,可采用若干方法,其中之一是镶边法。
10.3.2.2 二维视速度滤波
(1)二维视速度滤波的提出
在地震勘探中,有时有效波和干扰波的频谱成分十分接近甚至重合,这时无法利用频率滤波压制干扰,需要利用有效波和干扰波在其他方面的差异来进行滤波。如果有效波和干扰波在视速度分布方面有差异,则可进行视速度滤波。这种滤波要同时对若干道进行计算才能得到输出,因此是一种二维滤波。
地表接收的地震波动实际上是时间和空间的二维函数g(t,x),即是振动图和波剖面的组合,二者之间通过
勘查技术工程学
发生内在联系。式中k为空间波数,表示单位长度上波长的个数,f为频率,描述单位时间内振动次数,v为波速。
实际地震勘探总是沿地面测线进行观测,上述波数和速度应以波数分量kx和视速度v*代入。则有
勘查技术工程学
既然地震波动是空间变量x和时间变量t的二维函数,且空间和时间存在着密切关系,无论单独进行哪一维滤波都会引起另一维特性的变化(例如单独进行频率滤波会改变波剖面形状,单独进行波数滤波会影响振动图形,产生频率畸变),产生不良效果。那么只有根据二者的内在联系组成时间空间域(或频率波数域)滤波,才能达到压制干扰,突出有效波的目的。因此,应该进行二维滤波。
(2)二维视速度滤波的原理
二维滤波原理是建立在二维傅里叶变换基础上的。沿地面直测线观测到的地震波动g(t,x)是一个随时间和空间变化的波,通过二维正、反傅里叶变换得到其频率波数谱G(ω,kx)和时空函数。
勘查技术工程学
上式说明,g(t,x)是由无数个角频率为ω=2πf、波数为kx的平面简谐波所组成,它们沿测线以视速度v*传播。
如果有效波和干扰波的平面简谐波成分有差异,有效波的平面谐波成分以与干扰波的平面谐波成分不同的视速度传播如图10-8,则可用二维视速度滤波将它们分开,达到压制干扰,提高信噪比的目的。
(3)二维滤波的计算
图10-8 有效波和干扰波以不同成分平面简谐波的传播
二维线性滤波器的性质由其空间-时间特性h(t,x)或频率-波数特性H(ω,kx)所确定。同一维滤波一样,在时-空域中,二维滤波由输入信号g(t,x)与滤波
算子h(t,x)的二维褶积运算实现,在频率-波数域中,由输入信号的谱G(ω,kx)与滤波器的频率波数特性H(ω,kx)相乘来完成。
勘查技术工程学
由于地震观测的离散性和排列长度的有限性,必须用有限个(N个)记录道的求和来代替对空间坐标的积分。
勘查技术工程学
式中,n为原始道号;m为结果道号。
由(10.3-15)式可见,二维褶积可归结为对一维褶积的结果再求和。故测线上任一点处二维滤波的结果可由N个地震道的一维滤波结果相加得到。这时每一道用各自的滤波器处理,其时间特性hm-n(t)取决于该道与输出道之间的距离。沿测线依次计算,可以得到全测线上的二维滤波结果(图10-9)。
与理想一维滤波一样,理想二维滤波也要求在通放带内频率-波数响应的振幅谱为1,在通放带外为0,相位谱亦为0,即零相位滤波。因此,二维理想滤波器的频率-波数响应是正实对称函数(二维对称,即对两个参量均对称),空间时间因子必为实对称函数。二维滤波同样存在伪门现象和吉普斯现象,也可采用镶边法和乘因子法解决。因是二维函数,情况复杂得多,通常只采用减小采样间隔(包括时间采样间隔Δt 和频率采样间隔Δf)和增大计算点数(包括时、空二方向上的点数 M 和N)的方法。
图10-9 二维滤波计算示意图(N=5)
(4)扇形滤波
最常用的二维滤波是扇形滤波。它能滤去低视速度和高频的干扰。其频率波数响应为
勘查技术工程学
图10-10 扇形滤波器的频率波数响应
通放带在f-kx平面上构成由坐标原点出发,以f轴和kx轴为对称的扇形区域(图10-10)。因此这种滤波器称为扇形滤波器。
利用傅里叶反变换可求出其因子为
勘查技术工程学
当在计算机上实现运算时,需要离散化。对时间采样:t=nΔ,n=0,±1,±2,……,Δ为时间采样间隔,Δ=1/2fc。空间采样间隔即输入道的道间距Δx。
由标准扇形滤波器可以组构出既压制高视速度干扰,又压制低视速度干扰的切饼式滤波器,进而还可组构出同时压制高、低频干扰的带通扇形滤波器和带通切饼式滤波器。
在叠加前应用扇形滤波,压制的目标可以是面波、散射波、折射波或电缆振动产生的波。至于在叠加后的应用,则可压制从倾斜界面上产生的多次反射或侧面波。
10.3.3 提高纵向分辨率的反滤波处理
由地震波的传播理论可知,在介质中地震波是以地震子波的形式在地下传播。地面接收到的反射波地震记录是地层反射系数与地震子波的褶积。因此,地层相当一个滤波器,使反射系数序列变成了由子波组成的地震记录,降低了地震勘探的纵向分辨率。反滤波的目的就是要设计一个反滤波器,再对地震记录滤波,消除地层滤波的作用,提高地震记录的纵向分辨率。
由前所述,地震记录是地层反射系数序列r(t)与地震子波b(t)的褶积,即
勘查技术工程学
由于子波的问题,使高分辨率的反射系数脉冲序列变成了低分辨的地震记录,b(t)就相当地层滤波因子。为提高分辨率,可设计一个反滤波器,设反滤波因子为a(t),并要求a(t)与b(t)满足以下关系
勘查技术工程学
用a(t)对地震记录x(t)反滤波
勘查技术工程学
其结果为反射系数序列。以上即为反滤波的基本原理。
反滤波在具体实现时,核心是确定反滤波因子a(t)。由于地震子波的不确定性以及地震记录中噪音干扰的存在,实际中要确定精确的a(t)是非常困难的,甚至是不可能的。为此在不同的近似假设条件下,相继研究了很多种确定反滤波因子a(t)的方法,这些方法基本可以分为两大类:一类是先求取地震子波b(t),再根据b(t)求a(t);另一大类是直接从地震记录中求a(t)。每一类中又有很多不同的方法(就仅反滤波方法之多,说明了反滤波处理的难度)。下面就反滤波方法中具有代表性的几种反滤波进行讨论。
10.3.3.1 地层反滤波
地层反滤波属于先求子波b(t),再求a(t)的方法。该方法要求有测井资料以及较好的井旁地震记录道。首先由声波测井资料转换与井旁地震记录道x(t)相匹配的地层反射系数序列r(t),对r(t)及x(t)求其频谱可得频率域方程为
勘查技术工程学
即有
勘查技术工程学
式中B(ω)为子波b(t)的谱,再由子波与反滤波因子的关系有
勘查技术工程学
经反傅里叶变换得 a(t)。式中 A(ω)为反滤波因子的频谱。写成 z 变换,为 A(z)=,可见A(z)是一个有理分式,要使A(z)具有稳定性,分母多项式B(z)的根必须在单位圆外,即要求子波b(t)为最小相位。
利用测井和井旁地震道求取子波及反滤波因子,即可用该反滤波因子对测线的其他道进行反滤波。
10.3.3.2 最小平方反滤波
最小平方反滤波是最小平方滤波(或称维纳滤波、最佳滤波)在反滤波领域中的应用。
最小平方反滤波的基本思想在于设计一个滤波算子,用它把已知的输入信号转换为与给定的期望输出信号在最小平方误差的意义下是最佳接近的输出。
设输入信号为x(t),它与待求的滤波因子h(t)相褶积得到实际输出y(t),即y(t)=x(t)*h(t)。由于种种原因,实际输出y(t)不可能与预先给定的期望出(t)完全一样,只能要求二者最佳地接近。判断是否最佳接近的标准很多,最小平方误差准则是其中之一,即当二者的误差平方和为最小时,则意味着二者为最佳地接近。在这个意义下求出滤波因子h(t)所进行的滤波即为最小平方滤波。
若待求的滤波因子是反滤波因子a(t),对输入子波b(t)反滤波后的期望输出为d(t),实际输出为y(t),按最小平方原理,使二者的误差平方和为最小时求得的反滤波因子称为最小平方反滤波因子。用它对地震记录x(t)进行的反滤波为最小平方反滤波。
设输入离散信号为地震子波b(n)={b(0),b(1),…,b(m)},待求的反滤波因子a(n)={a(m0),a(m0+1),a(m0+2),……,a(m0+m)},m0为a(t)的起始时间,(m+1)为a(t)的延续长度,b(n)与a(n)的褶积为实际输出y(n),即
勘查技术工程学
为地震子波与期望输出的互相关函数。
根据最小平方原理,经推导即可得到最小平方反滤波的基本方程:
勘查技术工程学
式中,
勘查技术工程学
为地震子波的自相关函数,
勘查技术工程学
为地震子波与期望输出的互相关函数。
(10.3-24)式是一个线性方程组,写成矩阵形式为
勘查技术工程学
式中利用了自相关函数的对称性。该方程中,系数矩阵为一种特殊的正定矩阵,称为一般的托布里兹矩阵,该矩阵方程可用莱文森递推算法快速求解。
式(10.3-27)适应子波b(n)为最小相位、最大相位和混合相位。式中反滤波因子a(n)的起始时间m0与子波的相位有关,其取值规则由子波及反滤波因子的z变换确定。
10.3.3.3 预测反滤波
预测问题是对某一物理量的未来值进行估计,利用已知的该物理量的过去值和现在值得到它在未来某一时刻的估计值(预测值)的问题。它是科学技术中十分重要的问题。天气预报、地震预报、反导弹的自动跟踪等都属于这类问题。预测实质上也是一种滤波,称为预测滤波。
(1)预测反滤波原理
根据预测理论,若将地震记录x(t)看成一个平稳的时间序列,地震子波b(t)为物理可实现的最小相位信号,反射系数r(t)为互不相关的白噪声,由地震记录的褶积模型,在(t+α)时的地震记录x(t+α)为
勘查技术工程学
分析(10.3-28)式的第一项
勘查技术工程学
可见这一项是由反射系数r(t)的将来值决定的。若令第二项为
勘查技术工程学
^x(t+α)是 t 和t 以前时刻的r(t)值决定的,也就是说(t+α)可由现在和过去的资料预测,称(t+α)为预测值。求 x(t+α)与(t+α)的差值为
勘查技术工程学
ε(t+α)称为预测误差,或称为新记录。比较(10.3-28)及(10.3-29)两式,当预测值已知时,从原记录x(t+α)中减去预测值(t+α)后形成的新记录ε(t+α)中比原记录中涉及的反射系数少,与子波褶积后波形的干涉程度轻,波形易分辨,即分辨率提高了。
在上式中α称为预测距或预测步长。当α=1时,
勘查技术工程学
即有
勘查技术工程学
这时(t+1)时刻的预测误差与反射系数之间仅差一个常数b(0)。
因此,选预测距α=1,预测误差为反射系数,达到了反滤波的目的,此时称为预测反滤波。
当α>1时,预测误差为预测滤波结果,预测滤波主要用于消除多次波,尤其是消除海上鸣震。
(2)计算预测值(t+α)的方法
在预测滤波及预测反滤波中,关键是计算预测值(t+α),其方法如下。
由反滤波方程
勘查技术工程学
代入预测值(t+α)的表达式
勘查技术工程学
式中令τ=s-j,c(s)=b(j+α)a(s-j)称为预测因子。a(t)为反滤波因子。预测值(t+α)为预测因子 c(s)与地震记录的褶积。
现在需设计一个最佳预测因子c(s),使求取的预测值(t+α)与x(t+α)最接近,即使预测误差的平方和(误差能量)
勘查技术工程学
为最小。根据最小平方原理,可得线性方程组
勘查技术工程学
式中Rxx(τ)为地震记录的自相关函数
勘查技术工程学
T为相关时窗长度,m+1是预测因子长度。将(10.3-34)写成矩阵形式为
勘查技术工程学
解此方程组即可求得预测滤波因子c(t),用它对地震记录x(t)褶积可以求出未来时刻(t+α)时的最佳预测值(t+α)。