实验3:基于最佳维纳滤波器的盲解卷积算法
一．算法原理：
1.概论：
反褶积是通过压缩地震记录中的基本地震子波，压制交混回响和短周期多次波，从而提高时间分辨率，再现地下地层的反射系数。
反褶积通常应用于叠前资料，也可广泛用于叠后资料。
理想的反褶积应该压缩子波并消除多次波，在地震地道内只留下地层反射系数。
子波压缩可以通过将反滤波器作为反褶积算子来实现，它与地震子波做褶积时，反滤波器可以将地震子波转变成尖脉冲。
当应用于地震合成记录时，反滤波输出应为地层脉冲响应，精确的反滤波器设计可用最小平方模型来实现。
反褶积处理的基本假设是震源子波为最小相位。
2.褶积模型：
假设1：地层是由具有常速的水平层组成；
假设2：震源产生一个平面压缩波(P波)，法向入射到层边界上，在这种情况下，不产生剪切波(S波)；
假设3：震源波形在地下传播过程中不变，即它是稳定的；
数学上，褶积模型由下式给出：
x t w t e t n t
=+（3-1）
()()*()()
式中：()
n t为随机x t代表地震记录，()
e t为震源信号，()
w t为基本地震子波，()
噪声，*表示褶积。
反褶积试图从地震记录中恢复反射系数序列（严格的说是脉冲响应）。
假设4：噪音成分为零，于是式（3-1）变为
=（3-2）
x t w t e t
()()*()
假设5：震源波形是已知的；
假设6：反射系数序列是一个随机过程。
这意味着地震记录具有地震子波的特征，即它们的自相关和振幅谱是相似的；
假设7：地震子波是最小相位的，因此，它有一个最小相位的逆。
3.最佳维纳滤波器：
维纳滤波器是以最小平方误差为准则的，即要使下式最小：
设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳－霍夫（Wiener－Hopf）方程。
滤波器的维纳－霍夫方程如下：
（3-3）
式中，i r ，i a 和i g （0,1,2,...1i n =-）分别为输入子波的自相关、维纳滤波系数和期望输出与输入子波的互相关。
下图-1为维纳滤波器的设计和应用流程图：
图-1为维纳滤波器的设计和应用流程图
确定维纳滤波器的系数需要求解维纳－霍夫方程，由方程可以看到自相关矩阵是对称的。
这个特殊矩阵称作Toeplitz 矩阵，可用莱文逊递归法求解。
最佳维纳滤波器（0121,,,...,n a a a a -）是最佳的，是指它的实际输出与期望输出之间的最小平方误差最小。
当期望输出是零延迟尖脉冲（1,0,0，...，0）时，维纳滤波器与最小平方滤波器相同，即后者是前者的特例。
维纳滤波器可以考虑任一种期望输出而不仅限于零延迟尖脉冲。
期望输出可以有5类选择：
类型1：零延迟尖脉冲；
类型2：任一延迟尖脉冲；
类型3：时间提前了的输入序列；
类型4：零相位子波；
类型5：任意期望波形。
二． 常用子波：
地震资料处理中常用的子波有以下几种：
（1）Ricker 子波：
时域表达式：222222()(12)M f M w t f t e t ππ-=- （3-4）
频域表达式：22232()M f f M f W f e f π-= （3-5）
其中，M f 为子波的主频。
下图为主频M f =50Hz 的Ricker 的时域和频域波形。
图-2 Ricker 子波时域/频域波形
（2）Berlage 子波：
时域表达式为： 00()()cos(2)n t w t AH t t e f t απφ-=+ （3-6） 下图为主频为30Hz 的Berlage 子波的时域和频域波形：
图-3 Berlage 子波时域/频域波形
（3）一种常用的模拟子波：
时域表达式：00.12()sin()6.4t t t
w t A e π--= （3-7）
该子波对应的最小相位、最大相位、零相位和混合相位子波如下4图所示。
（a）最小相位子波（b）最大相位子波
（c）零相位子波（d）混合相位子波
图-4 模拟子波四种不同相位的时域波形
三．Matlab源程序及说明：
fs=10;ts=1/fs;%采样频率
N=1000;t=ts*(0:N-1);
t0 = 4; % 最小相位子波¨
x=sin(pi*t/6.4).*exp(-0.12*abs(t-t0));%输入子波
figure(1);subplot(2,2,1);plot(x);title('输入子波');grid on;
y=zeros(1,N);y(1)=1;%期望输出
subplot(2,2,2);plot(y);title('期望输出');grid on;
%--------求维纳滤波器的系数------
[Rx,lags]=xcorr(x);%输入信号的子相关函数
Rxx=toeplitz(Rx(N:2*N-1));%对称化自相关函数矩阵使之成为toeplitz矩阵Rxy=xcorr(x,y);%输入信号与期望信号的互相关函数
Rxy=Rxy(N:2*N-1);
h=(inv(Rxx)*Rxy')';%维纳滤波器系数
subplot(2,2,3);plot(h);title('维纳滤波器系数');grid on;
yy=conv(x,h); %实际输出
yy=yy(1: N);
subplot(2,2,4);plot(yy);title('实际输出');grid on;
error=norm(yy-y)^2 %误差的累积能量
四． 结果分析：
（一）五种不同的期望输出：
用以上的程序验证算法的正确性，以下五幅图为五种不同类型的期望输出下得到的维纳滤波器系数和实际的输出。
（输入子波为式（3-7）中的最小相位子波，如图4（a ）所示）
下图5的期望输出为零延迟尖脉冲，实际输出与期望输出的误差的平方和（能量）52.7710L -=⨯.
图5-零延迟尖脉冲
下图6的期望输出为延迟100ms 的尖脉冲，实际输出与期望输出的误差的平方和（能量）54.9310L -=⨯.
图6-延迟100ms 尖脉冲
下图7的期望输出为时间提前了的输入子波，实际输出与期望输出的误差的平方和（能量）31.310L -=⨯.
图7-时间提前了的输入子波
下图8的期望输出为零相位输入子波，实际输出与期望输出的误差的平方和（能量）43.8510L -=⨯.
图8-零相位子波
下图9的期望输出为主频为10Hz的Ricker子波，实际输出与期望输出的误差的平方和（能量）4
L-
=⨯.
2.310
图9-Ricker子波
从以上五种期望输出的结果可以看出，期望输出与实际输出的误差很小，从而说明最佳维纳滤波器的正确性。
（二）子波的相位对结果的影响：
1．期望输出为Ricker子波：
（a ）最小相位
误差的平方和4
2.310L -=⨯
（b ）混合相位1
误差的平方和0.0086L =
（c）混合相位2
L=
误差的平方和0.0669
（d）混合相位3
L=
误差的平方和26.1284
（e ）最大相位
误差的平方和29.9197L =
2．期望输出为零延迟尖脉冲：
（a ）最小相位
误差的平方和52.7710L -=⨯
（b）混合相位1
L=
误差的平方和0.0442
（c）混合相位2
L=
误差的平方和0.6166
（d）混合相位3
L=
误差的平方和0.8019
（e）最大相位
L=
误差的平方和9957
3．期望输出为延迟尖脉冲：
（a ）最小相位
误差的平方和0.0201L =
（b ）混合相位1
误差的平方和4
7.9010L -=⨯
（c ）混合相位2
误差的平方和4
5.2510L -=⨯
（d ）混合相位3
误差的平方和0.9671L =
（e）最大相位
L
误差的平方和0.9995
从以上三个例子可以看出，当输入子波的相位从最小相位逐渐变为最大相位时，处理结果会有变化：当输入子波的相位在最小相位附近时，实际输出和期望输出误差很小，即维纳滤波器效果较好；当输入子波的相位逐渐接近最大相位时，实际输出和期望输出误差较大，且越靠近最大相位，误差越大，即维纳滤波器效果较差。
因此在实际的地震资料处理时，总是假设震源子波为最小相位。