最小平方反褶积实验报告

实验四最小平方反褶积实验题目：已知两个地震模型，反射系数已给出。Re1为厚层模型，Re2为薄层模型，时间域采样间隔dt=1ms.要求用50Hz雷克子波先分别生成两个模型的合成地震道，然后根据最小平方反褶积的原理编写脉冲反褶积程序，并将反褶积后的结果与已知模型比较。可以对模型添加不同水平的随机噪声，检验其对反褶积算法的影响。实验内容：把延续几十至l00ms的地震子波压缩成原来的震源脉冲形式，地震记录变为反映反射系数序列的窄脉冲组合，这就是反滤波所要完成的工作。反褶积的目的就是为了把地震子波压缩成尖脉冲，使实际的地震记录变成反射系数序列。假设地震记录为0()()()()()()xtStntbtnt（1—1）其中()St为有效信号，()nt为干扰波。首先假设不存在干扰波()nt，即：()()()()xtStbtt（1—2）对两边求傅氏变换，则得到频率域的地震记录表示式：()()()XB（1—3）式中，()X、()B和()分别为地震频谱、子波频谱和反射系数的频谱。显然：1()()()XB（1—4）如果令：1()()AB（1—5）则有：()()()AX（1—6）再对（1—6）式做反傅氏变换至时间域，就可得到：()()()()()()tatxtatbtt（1—7）式中，()at为()A的时间函数。根据（1—7）式知：()()()atbtt（1—8）因为()bt为地震子波，而()at和()bt之间又存在着频谱互为倒数的关系（即1()()AB），由此可知，如已知地震子波，利用数学方法求出()at，再利用（1—7）式让反子波()at与地震记录()xt做褶积，就可以求出反射系数序列()t，即()()()taxt（1—9）经过这样的处理，就可以达到把地震子波压缩成尖脉冲，从而达到提高地震记录纵向分辨能力的目的。脉冲反褶积的基本思想在于设计一个滤波算子，用它把已知的输入信号转换为与给定的期望输出信号在最小平方误差的意义下是最佳接近的输出。若将地震子波作为反滤波的输入，期望输出则为尖脉冲。若设计另一滤波器输入信号()gt是某滤波器的输出，而期望输出()t是该滤波器的输入，则按此思想求得的滤波因子()at即称为脉冲反滤波因子，用它进行的滤波就是脉冲反滤波，即脉冲反褶积。先假设期望输出为窄脉冲()dt，在子波已知的情况下，设待求的反滤波因子()at起始时刻为0m，延续长度为(1)m。即0000()[(),(1),(2),,()]atamamamamm当已知输入——地震子波()[(0),(1),,()]btbbbn时，实际输出为00()()()()()mmmctatbtabt实际输出与期望输出的误差平方和为00002[()()()]mmnmmtmmQabtdt（1—10）要使Q为最小，数学上就是求Q的极值问题，即求满足0()Qal000(,1,,)lmmmm（1—11）的滤波因子()at。00()()()mmnbbtmbtbtlrl为地震子波的自相关函数，而00()()()mmnbdtmdtbtlrl为地震子波与期望输出的互相关函数，故（1—11）式可写为mmrmrmrmmamamarmrmrmrrrmrrrbdbdbdbbbbbbbbbbbbbbbbbb00000010110110（1—12）此方程系数矩阵即为拖布利兹矩阵。若期望输出是脉冲，则互相关为00()()()()mmnbdtmrltbtlbl（1—13）基本方程（1—12）变为mmbmbmbmmamamarmrmrmrrrmrrrbbbbbbbbbbbbbbbbbb000000110110110（1—14）一般情况下，地震子波为未知的，为在未知子波的情况下求出反滤波因子，必须对地震子波及反射系数序列加上一定的假设条件，他们包括：A.假设反射系数序列()Rt是随机的白噪序列，即其自相关为10()()0,RRr，其他（1—15）B.假设地震子波是最小相位的。根据假设A，地震子波的自相关()bbr可以用地震记录()xt的自相关()xxr代替。根据假设B，可知地震子波的Z变换()Bz的零点全部在单位圆外，也即反滤波因子()at的Z变换1()()AzBz的分母多项式的零点全在单位圆外，故()at是稳定的、物理可实现的。因此，00m，自由项变为[(0),(1),,()]Tbbbm。又因()bt必为物理可实现的，故(1)0b，(2)0b，，()0bm。令'()()(0)atatb，则基本方程变为001100110110'''maaarmrmrmrrrmrrrxxxxxxxxxxxxxxxxxx（1—16）这就是脉冲反褶积的基本方程，其系数矩阵中各元素可直接由地震记录求得。求出的反滤波因子'()at仅与()at相差常数倍，不影响压缩子波、提高分辨率的反滤波作用。当求取了反褶积因子()at后，令其与地震记录()xt进行褶积运算，即'()()()Statxt，则()St即为经过脉冲反褶积之后输出地新的地震记录。源程序：fm=50;r=3;n=200;dt=0.001;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-2.txt');R2=load('D:\study\re2-2.txt');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运行程序得图：地震子波合成地震记录（厚层模型、薄层模型）反褶积得到的反射系数与模型比较（厚层、薄层）利用合成地震记录自相关构造Toeplitz矩阵进行反褶积运算,得到的反射系数与模型比较：（厚层、薄层）在模型上添加白噪，使更符合地震记录自相关前提条件，进行反褶积运算自上而下：添加噪声(0.001*randn(size(R)))、(0.01*randn(size(R)))、(0.1*randn(size(R)))噪声较小（0.001）时，地震记录自相关得到的反射系数含有较多的毛刺，后面的部分甚至掩盖了模型中的反射系数；中间的图（0.1），后面的几个主要反射系数在计算结果上的体现，幅度要比毛刺大，辨别程度比上图高；而当白噪添加较多时，实际上全是白噪声（掩盖了R1）计算得到的反射系数与模型相差不大。地震记录添加噪声对反褶积结果的影响：噪声程度较小时，就已经影响到了反射系数的计算，但大体还能够反应部分反射系数的信息；而噪声程度较大时，计算得到的反射系数已经完全不能反映模型的情况。结论分析：反射系数与白噪相差较大时，用地震记录自相关代替子波自相关进行反褶积运算，得到的反射系数误差较大。厚层模型更明显，t增大，反褶积得到的反射系数与模型相差更大。与厚层模型相比，薄层模型由于较接近白噪（波形跳动较小，接近连续变化），地震记录自相关得到的反射系数与模型偏差较小。记录中的噪声对地震记录自相关反褶积计算的结果影响很大，信噪比较低时，尤其对于深层信息差别较大。