地震子波反褶积程序

盲信号实验报告盲解卷积算法姓名：丁宪成系别：电信学院专业：电磁场与微波学号：3110035012指导教师：陈文超2011年07月13日盲解卷积算法1.原理：几个重要概念：1.1褶积模型假设：假设1：地层是由具有常速的水平层组成；假设2：震源产生一个平面压缩波(P波)，法向入射到层边界上，在这种情况下，不产生剪切波(S波)；假设3：震源波形在地下传播过程中不变，即它是稳定的；假设4：噪音成分是零；假设5：震源波形是已知的；假设6：反射系数序列是一个随机过程。这意味着地震记录具有地；震子波的特征，即它们的自相关和振幅谱是相似的；假设7：地震子波是最小相位的，因此，它有一个最小相位的逆。1.2反滤波如果定义滤波算子为f(t),则f(t)与已知地震记录x(t)的褶积得到一个对地层脉冲响应e(t)的估计e(t)=f(t)∗x(t)；（1）x(t)=w(t)*f(t)*x(t)；（2）δ(t)=w(t)*f(t)；（3）1()()*()fttwt（4）用流程图表示为：1.3震源反子波计算震源反子波在数学上是利用z变换来实现的。例如，假设基本子波为两点时间序列（1，－0.5)1()12wzz2111()1...12412Fzzzz；（5）2111()1...12412Fzzzz（6）()Fz的系数11(1,,,...)24代表逆滤波算子f(t)有关的时间序列。可以看出它有无限多个系数，然而它们递减的很快。如同任何滤波过程一样，实际应用的算子都是被截断的。1.4最小平方反滤波当输入子波良好，其z变换的逆可以用一收敛序列表示，则上面所描述的反滤波将得到一个很好的近似于尖脉冲的输出.将下面问题列出方程式：给定输入子波(1－0.5)，寻找一个二项滤波器(a,b)，使实际输出与期望输出(1,0,0)之间的误差在最小二乘的意义上最小。将滤波器(a,b)与输入子波(1,-0.5)褶积以计算实际输出。误差的累积能量L定义为实际的与期望的输出的系数之差的平方和：222(1)()()22abLab（7）目的是寻找系数(a,b)使L最小，这要求变量L随系数(a,b)而变并使之为零对上式进行简化，取L对a和b的偏导数，并使其为零，得到：522502abab（8）有两个方程和两个未知数即滤波器系数(a,b)，可变成下列普通矩阵形式：5/21215/20ab（9）求解滤波系数，我们得到(a,b)为(0.95,0.38)。1.5最佳维纳滤波器回顾研究反滤波和最小平方滤波时讨论的期望输出，即零延迟尖脉冲(1,0,0)。重写最小平方反滤波的求解方程如下：5/41/2221/25/40ab（10）两边同时除以2得到：上述结论被维纳普遍化以推导出将输入转换为任意期望输出的滤波器(Robinson和Treltel，1980)。一个滤波器长度为n的矩阵方程的普遍形式是：（11）式中ir、ia和ig(i＝0,l,2,3,…,n-1)分别为输入子波自相关、维纳滤波系数和期望输出与输入子波的互相关。其计算过程的流程图可以用下面图形表示为：与最小平方滤波的联系：最佳维纳滤波器0121(,,,...)naaaa是最佳的，是指它的实际输出与期望输出之间的最小平方误差最小。当期望输出是零延迟尖脉冲(1,0,0,…,0)时，维纳滤波与最小平方滤波相同，即后者是前者的特例。1.6脉冲反褶积第一类期望输出(零延迟尖脉冲)的处理称为脉冲反褶积。期望脉冲(1,0,0,…,0)与输入子波011(,,...,)nxxx的互相关得到序列(1,0,0,…,0)。正则方程的一般形式改写成特殊形式：(12)方程已用作了归一化。在前面讨论过的最小平方反滤波器有与矩阵方程相同的形式。因此，脉冲反褶积在数学上是与最小平方反滤波相同的。而在实际上两类滤波器有一些区别。在脉冲反褶积(统计反褶积)的情况，方程在左边的自相关矩阵由输入地震记录(假设6)求得，而在最小平方反滤波(确定性反褶积)的情况则直接由已知的震源子波求得。结论：如果输入子波不是最小相位，则脉冲反褶积不能将它转化为完全的零延迟尖脉冲。虽然振幅谱实质上是平的，但输出的相位谱不是最小相位的。而且脉冲反褶积算子是输入子波的最小相位对应的子波的逆，这个子波可以是也可以不是最小相位的。1.7整形子波滤波器处理脉冲反褶积将子波(-0.5,1)压缩为零延迟尖脉冲(1,0,0)时有些困难。从能量分布来说，这个输入子波更接近于延迟尖脉冲，如(0,1,0)，而不是零延迟尖脉冲(1,0,0)。或者，将子波(-0.5,1)转换为延迟尖脉冲的滤波器会比将它整形为零延迟尖脉冲的滤波器产生更少的误差。根据最佳维纳滤波器的流程图重新设计和应用一个滤波器(期望输出为(0,1,0))。首先，计算互相关。我们知道输入子波的自相关。将互相关和子波自相关的结果代入前述矩阵方程，得到：(13)求解滤波器系数，结果为162(,)(,)2121ab。这个滤波器用于输入子波，如下表所示。和我们期望的一样，输出结果与最小平方滤波的输出结果一样。可以看出，对于延迟尖脉冲，实际输出与期望输出的最小平方误差为0.190；而对期望输出的零延迟尖脉冲，最小平方误差为0.762。这表明，将子波(－0.5,1)转换为延迟尖脉冲(0,1,0)比零延迟尖脉冲(1,0,0)误差小。2.实验内容地震记录信号是地震子波与地层脉冲响应相卷积的结果，如果我们要想得到脉冲信号必须进行盲解卷积，地震子波一般可以近似测得，我们要做的就是求出地震子波对应的逆信号，理想情况它们相卷积的结果为单位脉冲。地震子波有零相位子波，最小相位子波，混合相位子波和最大相位子波，如下图所示：程序说明：在下面的程序中，我们任意给定一个模拟地震子波，通过一系列的操作得到想要得到的结果。算法如下：第一步：产生模拟地震子波信号，给出预白百分比；第二步：给定期望输出，计算输入的自相关序列以及输入和期望输出的自相关序列；第三步：根据输入信号的自相关序列和预白百分比产生Toeplitz矩阵，根据Toeplitz矩阵和互相关序列求出维纳滤波器的滤波系数；第四步：通过将输入序列和滤波系数卷积求的实际输出。3.实验结果：3.1得到任意延迟脉冲信号：输入想得到的延迟脉冲的延迟时间:0（零延迟脉冲）输入想得到的延迟脉冲的延迟时间:5任意输出延迟时间，可以得到任意延迟的脉冲。3.2得到任意提前的输入序列运行结果：输入想得到的提前了的时间:2输入想得到的提前了的时间:10可见，提前的时间越长，得到的结果误差越来越大。3.3得到0相位子波运行结果：从对比图也可以看出在输入序列能量集中的地方得到的实际输出和期望输出误差比较小，其它区域误差较大。3.4得到任意期望输出：从对比图也可以看出，在输入子波能量集中的地方得到的实际结果与期望输出很近似，在输入子波幅度很小的区域得到的实际输出结果误差很大。3.5得到均方误差最小的实际输出运行结果：j=13;min1=0.0015.4．讨论地震子波的相位对处理结果的影响地震数据处理的目地是将野外采集的地震记录用处理模块进行处理后得到成像好，分辨率高的地震剖面，地震记录可描述为地震子波与地层脉冲响应或地下反射系数的褶积。就某种意义是那个讲，地震数据处理实际上就是一个对地震子波不断改造的过程。地震子波经过傅立叶变换之后可以得到振幅谱和相位谱，因此在地震记录中可通过拓宽地震子波的振幅谱来提高地震剖面的分辨率，也可以通过改变地震子波的相位谱来达到提高分辨率的目的。分辨率高的地震子波在频率域上表现为振幅谱尽可能的拓宽。而具有同样振幅谱的地震淄博中，最小相位的地震子波分辨率最高。实际地震数据处理中对地震子波相位特性有要求和影响的处理模块主要包括可控震源地震子波最小化，预滤波，叠前反褶积和叠后地震子波零相位化。我们知道地震数据处理的母的是提高地震记录的信噪比和分辨率，得到高信噪比，高分辨率的地震剖面。在频率域中表现为使得振幅频带尽可能宽且剩余子波是零相位的，在处理中有些模块对地震数据的子波相位要求有假设前提条件。一般情况下，叠前反褶积输入的地震数据要求子波是最小相位的，若是零相位的需要进行最小相位化，根据相同的振幅谱零相位的子波分辨率最高的原则，需要对叠后数据的子波进行零相位化处理。5附录：5.1得到想要的任意延迟脉冲t=0:1:30;x=sin(pi*(t-2)/6.4).*exp(-0.12*abs(t-2));%产生模拟子波subplot(211);plot(t,x);m=length(x);p=0.005;%预白百分比i=input('输入想得到的延迟脉冲的延迟时间:');y=[zeros(1,i),1,zeros(1,2*m-2-i)];%期望输出n=length(y);r_xx=xcorr(x);%计算子波的自相关A=fliplr(r_xx(1:m));R=toeplitz([(1+p/100)*A(1),A(2:length(A))]);%产生Toeplitz矩阵r_xy=xcorr(x,y);%输入和期望输出的互相关G1=fliplr(r_xy(1:n));G2=G1(1:m);h=inv(R)*G2';%h为维纳滤波器系数z=conv(x,h');%计算实际输出信号subplot(212);plot(z)%绘制满足最小均方误差的实际输出结果5.2得到时间提前了的输入序列说明：假设输入是（1，2，3，4，5，6），提前时间为2个单位，那么得到的期望输出序列为（3，4，5，6，0，0，…），所以上面的程序只要修改一下期望输出即可，将下面两行程序：i=input('输入想得到的延迟脉冲的延迟时间:');y=[zeros(1,i),1,zeros(1,2*m-2-i)];%期望输出修改为：i=input('输入想得到的提前了的时间:');y=[x(i+1:m),zeros(1,m+i-1)];%期望输出5.3得到零相位子波零相位子波：所谓零相位子波，就是子波序列进行傅立叶变换得到的结果全是正实数，因为正实数的相位为零，我们可以给定一个实偶序列，其傅立叶变换是正偶实数，则它肯定是零相位的，例如余弦序列，其傅立叶变换的结果就是正实数，所以是零相位的，那么程序就要修改为：t=0:1:30;x=sin(pi*(t-1)/6.4).*exp(-0.12*abs(t-1));%产生模拟子波subplot(311);plot(t,x);m=length(x);p=0.005;%预白百分比t1=-30:1:30;y=cos(pi*t1/30);%期望输出subplot(312)plot(y);n=length(y);r_xx=xcorr(x);%计算子波的自相关A=fliplr(r_xx(1:m));R=toeplitz([(1+p/100)*A(1),A(2:length(A))]);%产生Toeplitz矩阵r_xy=xcorr(x,y);%输入和期望输出的互相关G1=fliplr(r_xy(1:n));G2=G1(1:m);h=inv(R)*G2';%h为维纳滤波器系数z=conv(x,h');%计算实际输出信号subplot(313);plot(z)%绘制满足最小均方误差的实际输出5.4得到任意期望的输出期望输出可以用随机数产生，程序改为：t=0:1:30;x=sin(pi*(t-2)/6.4).*exp(-0.12*abs(t-2));%产生模拟子波subplot(311);plot(t,x);m=length(x);p=0.005;%预白百分比y=rand(1,2*m-1);%期望输出subplot(312);plot(y);n=length(y);r_xx=xcorr(x);%计算子波的自相关A=fliplr(r_xx(1:m));R=toeplitz([(1+p/100)*A(1),A(2:length(A))]);%产生Toeplitz矩阵r_xy=xcorr(x,y);%输入和期望输出的互