文档介绍：实验四最小平方反褶积实验题目:已知两个地震模型,反射系数已给出。Re1为厚层模型,Re2为薄层模型,时间域采样间隔dt=,然后根据最小平方反褶积的原理编写脉冲反褶积程序,并将反褶积后的结果与已知模型比较。可以对模型添加不同水平的随机噪声,检验其对反褶积算法的影响。实验内容:把延续几十至l00ms的地震子波压缩成原来的震源脉冲形式,地震记录变为反映反射系数序列的窄脉冲组合,这就是反滤波所要完成的工作。反褶积的目的就是为了把地震子波压缩成尖脉冲,使实际的地震记录变成反射系数序列。假设地震记录为(1—1)其中为有效信号,为干扰波。首先假设不存在干扰波,即:(1—2)对两边求傅氏变换,则得到频率域的地震记录表示式:(1—3)式中,、和分别为地震频谱、子波频谱和反射系数的频谱。显然:(1—4)如果令:(1—5)则有:(1—6)再对(1—6)式做反傅氏变换至时间域,就可得到:(1—7)式中,为的时间函数。根据(1—7)式知:(1—8)因为为地震子波,而和之间又存在着频谱互为倒数的关系(即),由此可知,如已知地震子波,利用数学方法求出,再利用(1—7)式让反子波与地震记录做褶积,就可以求出反射系数序列,即(1—9)经过这样的处理,就可以达到把地震子波压缩成尖脉冲,从而达到提高地震记录纵向分辨能力的目的。脉冲反褶积的基本思想在于设计一个滤波算子,用它把已知的输入信号转换为与给定的期望输出信号在最小平方误差的意义下是最佳接近的输出。若将地震子波作为反滤波的输入,期望输出则为尖脉冲。若设计另一滤波器输入信号是某滤波器的输出,而期望输出是该滤波器的输入,则按此思想求得的滤波因子即称为脉冲反滤波因子,用它进行的滤波就是脉冲反滤波,即脉冲反褶积。先假设期望输出为窄脉冲,在子波已知的情况下,设待求的反滤波因子起始时刻为,延续长度为。即当已知输入——地震子波时,实际输出为实际输出与期望输出的误差平方和为(1—10)要使Q为最小,数学上就是求Q的极值问题,即求满足(1—11)的滤波因子。为地震子波的自相关函数,而为地震子波与期望输出的互相关函数,故(1—11)式可写为(1—12)此方程系数矩阵即为拖布利兹矩阵。若期望输出是脉冲,则互相关为(1—13)基本方程(1—12)变为(1—14)一般情况下,地震子波为未知的,为在未知子波的情况下求出反滤波因子,必须对地震子波及反射系数序列加上一定的假设条件,他们包括:,即其自相关为(1—15)。根据假设A,地震子波的自相关可以用地震记录的自相关代替。根据假设B,可知地震子波的Z变换的零点全部在单位圆外,也即反滤波因子的Z变换的分母多项式的零点全在单位圆外,故是稳定的、物理可实现的。因此,,自由项变为。又因必为物理可实现的,故,,,。令,则基本方程变为(1—16)这就是脉冲反褶积的基本方程,其系数矩阵中各元素可直接由地震记录求得。求出的反滤波因子仅与相差常数倍,不影响压缩子波、提高分辨率的反滤波作用。当求取了反褶积因子后,令其与地震记录进行褶积运算,即,则即为经过脉冲反褶积之后输出地新的地震记录。源程序:fm=50;r=3;n=200;dt=;t=[0:1:n-1]*dt;w=exp(-(2*pi*fm/r*t).^2).*sin(2*pi*fm.*t);figure(1),plot(t,w,'k')title('wavelet')R1=load('D:\study\re1-');R2=load('D:\study\re2-');A1=conv(w,R1);A2=conv(w,R2);figure(2),plot(A1,'k'),holdonfigure(3),plot(A2,'r'),holdoffy1=xcorr(w);y1=y1(floor((length(y1)+1)/2):floor((length(y1)+1)/2)+length(w)-1);y2=xcorr(w);y2=y2(floor((length(y2)+1)/2):floor((length(y2)+1)/2)+length(w)-1);y1=syme(y1);y2=syme(y2);d1=[1,zeros(1,floor(length(w))-1)];d2=[1,zeros(1,floor(length(w))-1)];a1=LDLt(y1,d1);a2=LDLt(y2,d2);record1=conv(a1,A1);record2=conv(a2,A2);figure(4),plot(record1,'k'),holdonfigure(5),plot(record2,'r'),holdoff运