a b L (a 1) 2 (b ) 2 ( ) 2 2 2
（7）
目的是寻找系数(a,b)使L最小， 这要求变量L随系数(a,b)而变并使之为 零对上式进行简化，取L对a和b的偏导数，并使其为零，得到：
5 a b 2 2 a 5 b 0 2
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'; z=conv(x,h'); subplot(212); plot(z) %绘制满足最小均方误差的实际输出结果 %h为维纳滤波器系数 %计算实际输出信号 %输入和期望输出的互相关 %计算子波的自相关
期望输出很近似， 在输入子波幅度很小的区域得到的实际输出结果误 差很大。 3.5 得到均方误差最小的实际输出 运行结果： j =13; min1=0.0015.
4．讨论地震子波的相位对处理结果的影响 地震数据处理的目地是将野外采集的地震记录用处理模块进行 处理后得到成像好，分辨率高的地震剖面，地震记录可描述为地震子 波与地层脉冲响应或地下反射系数的褶积。就某种意义是那个讲，地 震数据处理实际上就是一个对地震子波不断改造的过程。 地震子波经 过傅立叶变换之后可以得到振幅谱和相位谱， 因此在地震记录中可通 过拓宽地震子波的振幅谱来提高地震剖面的分辨率， 也可以通过改变
(0,1,0))。首先，计算互相关。我们知道输入子波的自相关。将互相关 和子波自相关的结果代入前述矩阵方程，得到：
(13) 求解滤波器系数，结果为 (a, b) (
16 2 , ) 。这个滤波器用于输入子 21 21
波，如下表所示。和我们期望的一样，输出结果与最小平方滤波的输 出结果一样。
输入想得到的提前了的时间:10
可见，提前的时间越长，得到的结果误差越来越大。
3.3 得到0相位子波 运行结果：
从对比图也可以看出在输入序列能量集中的地方得到的实际输出和 期望输出误差比较小，其它区域误差较大。 3.4 得到任意期望输出：
从对比图也可以看出， 在输入子波能量集中的地方得到的实际结果与
盲信号实验报告
盲解卷积算法
姓 系 专 学
名： 别： 业： 号：
丁宪成 电信学院 电磁场与微波 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)； x(t) = w(t)* f (t)* x(t)； δ(t) = w(t)* f (t)； （1） （2） （3） （4）
可以看出，对于延迟尖脉冲，实际输出与期望输出的最小平方误差 为0.190；而对期望输出的零延迟尖脉冲，最小平方误差为0.762。这 表明， 将子波(－0.5,1)转换为延迟尖脉冲(0,1,0)比零延迟尖脉冲(1,0,0) 误差小。 2. 实验内容 地震记录信号是地震子波与地层脉冲响应相卷积的结果， 如果我 们要想得到脉冲信号必须进行盲解卷积，地震子波一般可以近似测 得，我们要做的就是求出地震子波对应的逆信号，理想情况它们相卷 积的结果为单位脉冲。地震子波有零相位子波，最小相位子波，混合
（8）
有两个方程和两个未知数即滤波器系数(a,b)，可变成下列普通矩阵形 式：
5 / 2 1 a 2 1 5 / 2 b 0
求解滤波系数，我们得到(a,b)为(0.95,0.38)。 1.5 最佳维纳滤波器
（9）
回顾研究反滤波和最小平方滤波时讨论的期望输出， 即零延迟尖脉冲 (1,0,0)。重写最小平方反滤波的求解方程如下：
（5）
F ( z)
1 1 1 1 z z 2 ... 1 2 4 ) 的系数 (1, , ,...) 代表逆滤波算子f(t)有关的时间序列。可以看 2 4
出它有无限多个系数， 然而它们递减的很快。 如同任何滤波过程一样，
实际应用的算子都是被截断的。 1.4 最小平方反滤波 当输入子波良好，其z变换的逆可以用一收敛序列表示，则上面所描 述的反滤波将得到一个很好的近似于尖脉冲的输出.将下面问题列出 方程式：给定输入子波(1－0.5)，寻找一个二项滤波器(a,b)，使实际 输出与期望输出(1,0,0)之间的误差在最小二乘的意义上最小。将滤波 器(a,b)与输入子波(1,-0.5)褶积以计算实际输出。误差的累积能量L定 义为实际的与期望的输出的系数之差的平方和：
相位子波和最大相位子波，如下图所示：
程序说明：在下面的程序中，我们任意给定一个模拟地震子波，通过 一系列的操作得到想要得到的结果。算法如下： 第一步：产生模拟地震子波信号，给出预白百分比； 第二步：给定期望输出，计算输入的自相关序列以及输入和期望输出 的自相关序列； 第三步：根据输入信号的自相关序列和预白百分比产生Toeplitz矩阵，
(12) 方程已用作了归一化。 在前面讨论过的最小平方反滤波器有与矩阵方 程相同的形式。因此，脉冲反褶积在数学上是与最小平方反滤波相同 的。而在实际上两类滤波器有一些区别。在脉冲反褶积(统计反褶积) 的情况， 方程在左边的自相关矩阵由输入地震记录(假设6)求得， 而在 最小平方反滤波(确定性反褶积)的情况则直接由已知的震源子波求 得。 结论：如果输入子波不是最小相位，则脉冲反褶积不能将它 转化为完全的零延迟尖脉冲。虽然振幅谱实质上是平的，但 输出的相位谱不是最小相位的。 而且脉冲反褶积算子是输入子波的最 小相位对应的子波的逆，这个子波可以是也可以不是最小相位的。 1.7 整形子波滤波器处理 脉冲反褶积将子波(-0.5,1)压缩为零延迟尖脉冲(1,0,0)时有些困难。从 能量分布来说，这个输入子波更接近于延迟尖脉冲，如(0,1,0)，而不 是零延迟尖脉冲(1,0,0)。或者，将子波(-0.5,1)转换为延迟尖脉冲的滤 波器会比将它整形为零延迟尖脉冲的滤波器产生更少的误差。 根据最 佳维纳滤波器的流程图重新设计和应用一个滤波器(期望输出为
根据Toeplitz矩阵和互相关序列求出维纳滤波器的滤波系数； 第四步：通过将输入序列和滤波系数卷积求的实际输出。 3. 实验结果： 3.1 得到任意延迟脉冲信号： 输入想得到的延迟脉冲的延迟时间:0（零延迟脉冲）
输入想得到的延迟脉冲的延迟时间:5
任意输出延迟时间，可以得到任意延迟的脉冲。
3.2 得到任意提前的输入序列 运行结果： 输入想得到的提前了的时间:2
f (t ) (t ) *
1 w(t )
用流程图表示为：
1.3 震源反子波 计算震源反子波在数学上是利用 z 变换来实现的。例如，假设基本子 波为两点时间序列（1，－0.5)
1 1 1 1 w( z ) 1 z F ( z ) 1 z z 2 ... ； 1 2 2 4 1 z 2
地震子波的相位谱来达到提高分辨率的目的。 分辨率高的地震子波在 频率域上表现为振幅谱尽可能的拓宽。 而具有同样振幅谱的地震淄博 中，最小相位的地震子波分辨率最高。实际地震数据处理中对地震子 波相位特性有要求和影响的处理模块主要包括可控震源地震子波最 小化，预滤波，叠前反褶积和叠后地震子波零相位化。我们知道地震 数据处理的母的是提高地震记录的信噪比和分辨率，得到高信噪比， 高分辨率的地震剖面。 在频率域中表现为使得振幅频带尽可能宽且剩 余子波是零相位的， 在处理中有些模块对地震数据的子波相位要求有 假设前提条件。一般情况下，叠前反褶积输入的地震数据要求子波是 最小相位的，若是零相位的需要进行最小相位化，根据相同的振幅谱 零相位的子波分辨率最高的原则， 需要对叠后数据的子波进行零相位 化处理。 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; %预白百分比
（11） 式中 ri 、ai 和 gi (i＝0,l,2,3,…,n-1)分别为输入子波自相关、 维纳滤波系 数和期望输出与输入子波的互相关。 其计算过程的流程图可以用下面 图形表示为：
与最小平方滤波的联系：最佳维纳滤波器 (a0 , a1 , a2 ,...an1 ) 是最佳的， 是指它的实际输出与期望输出之间的最小平方误差最小。 当期望输出 是零延迟尖脉冲(1,0,0,…,0)时，维纳滤波与最小平方滤波相同，即后 者是前者的特例。 1.6 脉冲反褶积 第一类期望输出(零延迟尖脉冲)的处理称为脉冲反褶积。期望脉冲 (1,0,0,…,0)与输入子波 ( x0 , x1 ,..., xn1 ) 的互相关得到序列(1,0,0,…,0)。 正则方程的一般形式改写成特殊形式：
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 得到零相位子波 零相位子波：所谓零相位子波，就是子波序列进行傅立叶变换得到的 %期望输出