数字滤波与反滤波

广义而言，滤波可以看作为对某一信号的改造作用。改造之前的信号称为滤波输入，改造之后的信号称为滤波输出。输入散局、输出、滤波器就构成了滤波三要素（图5-4-6）。反射波地震资料数字处理中需改造的信号是含有干扰的地震信号，输出是不含干扰或干扰减少的地震信号。改造信号的滤波既可以用物理过程实现，也可以用数学运算实现。前者是所谓的电滤波器，后者即为数字滤波。

图5-4-6 滤波三要素

滤波器的设计一般在频率域中进行。首先利用傅里叶变换分析有效波、干扰波的频率成分。据此确定有效波、干扰波的频谱范围。例如，如图5-4-7所示，有效波的频率成分在中间频率范围内，干扰波分布在高、低频范围内。据此确定一个带通滤波器的频率响应H（ω），然后经傅里叶反变换得到h（t）在时间域中实现滤波，或者直接在频率域中进行滤波后再反变换得到滤波输出。如图5-4-7，所设计的滤波器称为门式滤波，边缘十分陡，滤波参数只有二个：高截止频率和低截止频率。实际使用的滤波器边缘较为平缓，有一定坡度，故有4个滤波参数：f1、f2、f3和f4。f1和f2之间为低频过渡带；f3和f4之间为高频过渡带；小于f1或大于f4时H（ω）＝0，f2和f3之间H（ω）＝1。设计一维频率滤波器就是根据干扰波和有效波的频谱分布确定这四个滤波参数。

图5-4-7 滤波器的设计H（ω）

所谓反滤波，仍然是一个滤波过程，只不过是一种特殊的滤波过程。这个滤波过程是针对另外某一个滤波过程而设计的，其作用恰好与另外那个滤波过程的作用相反隐绝，将该滤波过程的作用抵消。

设x（t）是滤波器h（t）的冲携让输入信号，y（t）为其输出。若设计一个滤波器a（t），使得当输入信号为y（t）时的输出正好是x（t）（图5-4-8），滤波器a（t）的作用与滤波器h（t）的作用正好相反，a（t）将h（t）的作用抵消［即若将h（t）和a（t）看作为一个滤波器的话，则输入为x（t），输出仍为x（t）］。因此，a（t）就称为h（t）的反滤波器。当然，h（t）也是a（t）的反滤波器。

图5-4-8 反滤波的概念

由此可见，反滤波实际上是一种特殊的滤波，有极强的针对性，必须有“反”的对象。如果失去了“反”的对象，则反滤波就失去了意义。反射波地震资料数字处理中正是根据反射地震勘探的实际情况，针对不同的对象设计不同的反滤波以达到不同的目的。

%用phantom函数可以获得仿体图像；

%用randon可获得不同角度的一维投影；

clear all

P = phantom('Modified Shepp-Logan',256)

R=radon(P)

figureimshow(R,[])

figure

imshow(P,[])title('仿体图')

%直接反投影法

l = pow2(nextpow2(size(R,1))-1)%重构图像的大小

P_1 = zeros(l,l)%用于存放重构后的图像

for i=1 : size(R,2)

    tmp = imrotate( repmat(R(:,i),1,size(R,1)),i-1,'bilinear' )

    tmp = tmp(floor(size(tmp,1)/2-l/2)+1:floor(size(tmp,1)/2+l/2),floor(size(tmp,2)/2-l/2)+1:floor(size(tmp,2)/2+l/2))

    P_1=P_1+tmp

end

P_1=P_1/size(R,2)

P_1=rot90(P_1)

figureimshow(P_1,[])title('直接反投影法')

%滤波反投影法

N=180

%滤波

H=size(R,1)

h=zeros((H*2-1),1)

for i=0:H-1

    if i==0

        h(H-i)=1/4

    elseif rem(i,2)==0

        h(H-i)=0

  姿御      h(H+i)=0

    else

        h(H-i)=-1/(i*pi)^2

        h(H+i)=-1/(i*pi)^2

  迹昌岩  end

end

x=zeros(H,N)

for i=1:N

    s=R(:,i)

    xx=conv(s',h')

    x(:,i)=xx(H:2*H-1)

end

%反投影

P_3=zeros(l,l)

for i=1:l

    for j=1:l

        for k=1:180

            theta=k/180*pi

            t=(j-l/2-0.5)*cos(theta)+(l/2+0.5-i)*sin(theta)+(H+1)/2

            t1=floor(t)

            t2=floor(t+1) 

            P_3(i,j)=P_3(i,j)+(t2-t)*x(t1,k)+(t-t1)*x(t2,k)

        end

    end

end

P_3=pi/N*P_3

figureimshow(P_3,[])title('滤波反投迅扮影法')

(一)一般反滤波

所谓反滤波仍然是一个滤波过程，这种滤波过程的作用恰好与某个其他滤波过程的作用相反。

设x(t)是时间函数为h(t)的滤波器的输入，y(t)为输出，则有

地震勘探

现在设计一滤波器a(t)，使得当y(t)作为其输入时，得到的输出一定是x(t)，即

地震勘探

则a(t)就是h(t)的反滤波。此过程用图4-26表示。

图4-26 反滤波的概念

将(4-48)式代入(4-47)式得

y(t)=y(t)*a(t)*h(t)

根据δ函数理论，有

y(t)=y(t)*δ(t)

因此a(t)*h(t)=δ(t)

此式给出了滤波因子h(t)与反滤波禅拦因子a(t)之间的时间域关系。在频率域中

A(ω)H(ω)=1

地震勘探

写成z变换的形式，有

A(z)H(z)=1

和 这是一个有理分式形式的z变换，H(z)为一多项式。由简单的代数知识可知，A(z)=1/H(z)一般是一个无穷级数。那么，A(t)一般应是一个无穷序列。

如上一节的分析，当H(z)的所有零点均不在单位圆上时，A(z)是稳定戚敬的，但只有当H(z)的所有零点均在单位圆外时，它才是物理可实现的。H(z)的所有零点均在单位圆外意味着h(t)是最小相位的。

(二)地震勘探反滤波

地震勘探反滤波的主要任务是抵消大地滤波作用，其中包括地震记录道中各种装备(如检波器、记录仪，它们都可以看成是一种滤波装置)对地震子波的滤波作用，从而提高纵向分辨率，这是本节主要要讨论的内容。

当然，某些规则干扰波(如海上交混回响)的形成过程也可看作滤波过程。针对这种滤波过程设计的反滤波目的在于压制它们，提高信噪比。这类特殊问题将在本节最后作些简介。

研究反滤波就是研究如何设计一个滤波器去抵消另一个滤波器的作用。通常有两种方法用来设计反滤波器，即确定性方法和统计方法。式(4-50)可以作为确定地震勘探反滤波因子的一个十分简单的代数方法，是一种确定性方法。但是，利用此式确定反滤波因子存在两个问题:①求A(z)必须事先已知B(z)，即大地滤波因子b(t)，在地震勘探中这一点往往难以做到②从理论上说，a(t)一般是一个无穷序列，实际工作中只能取有限项，所以反滤波因子只能近似地确定。如何近似可以有不同的办法。利用(4-50)式取近似只能人为地截断，其结果不一定最好。因此，在地震勘探实际中往往利用统计方法求取最佳的滤波因子。下面介绍几种高袭慎常用的统计反滤波方法。

---