盲解卷积和详细程序

实验3:基于最佳维纳滤波器的盲解卷积算法算法原理：概论：反褶积是通过压缩地震记录中的基本地震子波，压制交混回响和短周期多次波，从而提高时间分辨率，再现地下地层的反射系数。反褶积通常应用于叠前资料，也可广泛用于叠后资料。理想的反褶积应该压缩子波并消除多次波，在地震地道内只留下地层反射系数。子波压缩可以通过将反滤波器作为反褶积算子来实现，它与地震子波做褶积时，反滤波器可以将地震子波转变成尖脉冲。当应用于地震合成记录时，反滤波输出应为地层脉冲响应，精确的反滤波器设计可用最小平方模型来实现。反褶积处理的基本假设是震源子波为最小相位。褶积模型：假设1：地层是由具有常速的水平层组成；假设2：震源产生一个平面压缩波(P波)，法向入射到层边界上，在这种情况下，不产生剪切波(S波)；假设3：震源波形在地下传播过程中不变，即它是稳定的；数学上，褶积模型由下式给出：x(t)=w(t)*e(t)+n(t)(3-1)式中：x(t)代表地震记录，w(t)为基本地震子波，e(t)为震源信号，n(t)为随机噪声，*表示褶积。反褶积试图从地震记录中恢复反射系数序列(严格的说是脉冲响应)。假设4：噪音成分为零，于是式(3-1)变为x(t)=w(t)*e(t)(3-2)假设5：震源波形是已知的；假设6：反射系数序列是一个随机过程。这意味着地震记录具有地震子波的特征，即它们的自相关和振幅谱是相似的；假设7：地震子波是最小相位的，因此，它有一个最小相位的逆。最佳维纳滤波器：维纳滤波器是以最小平方误差为准则的，即要使下式最小：设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳一霍夫(Wiener-Hopf)方程。滤波器的维纳一霍夫方程如下：莒。升...rn_2口i幻—J_如一式中，r,a和g(i=0,1,2,...n-1)分别为输入子波的自相关、维纳滤波系iii数和期望输出与输入子波的互相关。下图-1为维纳滤波器的设计和应用流程图:图-1为维纳滤波器的设计和应用流程图确定维纳滤波器的系数需要求解维纳一霍夫方程，由方程可以看到自相关矩阵是对称的。这个特殊矩阵称作Toeplitz矩阵，可用莱文逊递归法求解。最佳维纳滤波器(a,a,a，…,a)是最佳的，是指它的实际输出与期望输012n-1出之间的最小平方误差最小。当期望输出是零延迟尖脉冲(1,0,0，...，0)时，维纳滤波器与最小平方滤波器相同，即后者是前者的特例。维纳滤波器可以考虑任一种期望输出而不仅限于零延迟尖脉冲。期望输出可以有5类选择：类型1：零延迟尖脉冲；类型2：任一延迟尖脉冲；类型3：时间提前了的输入序列；类型4：零相位子波；类型5：任意期望波形。常用子波：地震资料处理中常用的子波有以下几种：Ricker子波：时域表达式：w(t)=(1-2兀2f212)e-兀2fM12(3-4)M

（3-5）频域表达式：W（f）=土fe（3-5）其中，八为子波的主频。下图为主频fM=50Hz的Ricker的时域和频域波形。FourierSpectrumDfFourierSpectrumDfRicker图-2Ricker子波时域/频域波形（2）Berlage子波：(3-6)时域表达式为：w(t)=AH(t)tne-bcos(2兀ft+^)下图为主频为30Hz的Berlage(3-6)Berlage干波觥闾油形Berime子波时嫩波影10.2T|n>&/m&Fr&quenc^H2图-3Berlage子波时域/频域波形（3）一种常用的模拟子波：Berlage干波觥闾油形Berime子波时嫩波影10.2T|n>&/m&Fr&quenc^H2图-3Berlage子波时域/频域波形时域表达式：w（t）=Asin（旦）e-0.12t-t0（3-7）6.4该子波对应的最小相位、最大相位、零相位和混合相位子波如下4图所示。

(b)最大相位子波(a)最小相位子波(b)最大相位子波(a)最小相位子波(d)混合相位子波(c)零相位子波图-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('输入子波');gridon;y=zeros(1,N);y(1)=1;%期望输出subplot(2,2,2);plot(y);title('期望输出');gridon;%求维纳滤波器的系数[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('维纳滤波器系数');gridon;yy=conv(x,h);%实际输出yy=yy(1:N);subplot(2,2,4);plot(yy);title('实际输出');gridon;

error=norm(yy-y)八error=norm(yy-y)八2%误差的累积能量结果分析：（一）五种不同的期望输出：用以上的程序验证算法的正确性，以下五幅图为五种不同类型的期望输出下得到的维纳滤波器系数和实际的输出。（输入子波为式（3-7）中的最小相位子波，如图4（a）所示）下图5的期望输出为零延迟尖脉冲，实际输出与期望输出的误差的平方和（能量）L=2.77x10-5.oooO224■-cpnwdE<-60-0.52004006008001000oooO224■-cpnwdE<-60-0.52004006008001000Time/ms维纳波波器系数200400600BOOWOOTime/ms实际输出输入子波（最小相位）期里输出（零延退单位脉冲）5oQ省m_-dm<02004006008001000Time/ms8642Jo.0.O0>prl三dE<2004006008001000下图6的期望输出为延迟100ms的尖脉冲，实际输出与期望输出的误差的平方和（能量）L=4.93x10-5.

输入子披（最小相惊）200400600800输入子披（最小相惊）20040060080010001—I■■08060.402°200400600BOO1000TimMrs实际馆出0.806040.2QJ.——2004006008001000TimMrs图6-延迟100ms尖脉冲下图7的期望输出为时间提前了的输入子波，实际输出与期望输出的误差的输入了波（最小捌位一）靳咀输出〔时间提值了的输入序刊）200JOO6008001000实味输出平方和（能量）L=1.3x10-3.输入了波（最小捌位一）靳咀输出〔时间提值了的输入序刊）200JOO6008001000实味输出图7-时间提前了的输入子波下图8的期望输出为零相位输入子波，实际输出与期望输出的误差的平方和（能量）L=3.85x10-4.

锥熟源赧羿系数推入子波（最小相仲）0.S0.60.4020-02-042D0400锥熟源赧羿系数推入子波（最小相仲）0.S0.60.4020-02-042D04006008001000Time/ms200JOO6008001000O1..234tlpnwdE<下图9的期望输出为主频为10Hz的Ricker子波，实际输出与期望输出的误差的平方和（能量）L=2.3x10-4.输入子波（最小捌位）5OosspnMELE<-0.5—输入子波（最小捌位）5OosspnMELE<20D4口06003001000Timem期营输出(RickerT-^)2004006008001000Time/ms实际输出;5Otupm=dulvI-0.05200JOO600800100066期营输出(RickerT-^)2004006008001000Time/ms实际输出;5Otupm=dulvI-0.05200JOO6008001000664200O.O.O.命pm_o.E<-02■0.42004006008001000Time/ms从以上五种期望输出的结果可以看出，期望输出与实际输出的误差很小，从而说明最佳维纳滤波器的正确性。（二）子波的相位对结果的影响:1.期望输出为Ricker子波：期营输出(Ri如子波〉200400600800Time/ms实际输出-SotuTFm云mv输人子波20Q4口06003001000Timeyms推纳源波普系数5oosQDpnMeLE<10001000400600800Time/ms2008.6.4.2Q0.O.0.O.命pnwE<-02■0.4-005200JOO6008001000(a)最小相位误差的平方和乙=2.3x10-4期堕鹿出(HidwrTfe)200400前0800Time/ms实际输出5oaupm=dmv1000Time/ms场入子波20D4口06003001000Time/mso5-O.如pn-i—«--eLE<推纳凋波羿系敖200JOO6008001000输人子波SOS□.o.一spnl-a.E<200400前0800Time/ms实环箱出期箜输出(RickerT-Wl)2004Q0600800Timeynns

赧纳胡波掘系数31000100021O-1-2

由pn*一一dm<0>pnwE<.5o.WOO400600800误差的平方和£=0.0669物氏子波Q.51000200JU400600800Time/ms1000(c)混合相位220040060030。1000Timeyrns罪.纳旅波牌系疝o5o.-aipnl--dE<2004006008001000Time/nrs实廊箱出5050-1-ll>pm_D.E<误差的平方和£=26.12842004006003001000Time/ms(d)混合相位3-0.5物人子波50o.QDpnl--ELE<4006003001000Timem始汕源波普系徂lol0.0.-心pm=dE<njo1o知00:5l.00J002期营输击(Ricker7^)200Time/rns6008001000实廊输出O2-名2=dE<(e)最大相位误差的平方和£=29.91972.期望输出为零延迟尖脉冲：输氏子波期堡输出CRicke「芋波)-0.5—11120Q4QO6003001000Timem始油源波牌系徂2D04006008D0-0.5—11120Q4QO6003001000Timem始油源波牌系徂2D04006008D01000Time/ms实际输出.8.6.4,20.0.。注apm=dmv20ooO24-■瑁3_q.e<-6020040060080010002004006008001000Time/ms.864.20.0O注mpaKEV误差的平方和£=2.77x10-5期堕蛔CRiCfc&l子波》0.520Q4QQ60D3001000Timem

推汕源波普系数2D04006008D01000Time/ms实际输出50g3Q_E<■10D02004006008001000Time/ms0.520Q4QQ60D3001000Timem

推汕源波普系数2D04006008D01000Time/ms实际输出50g3Q_E<■10D02004006008001000Time/ms8.B4.2O0..O.0.0>p3=dE<200JOO6008001000误差的平方和£=0.0442输氏子波0.5oSO.■sipnl_-dE<(b)混合相位1期堑输他LRieke「子波)10.60.40.220040060。30。1000Timeyrns始汕滤波牌系E:2D04006003DQ1000Time/ms实际输出20oo-2a>pm_Q.E<0.302010-01-02200JOO600MO10002004006008001000Time/ms(c)混合相位2误差的平方和£=0.6166

QDpH--eLE<输氏于波20Q4口06003001000Timeynns■纳止■波牌系数期枣输出(Ri如子波》2D04006008001000Time/ms.6.4,20,0..O.apm=dmv200如0600QDpH--eLE<输氏于波20Q4口06003001000Timeynns■纳止■波牌系数期枣输出(Ri如子波》2D04006008001000Time/ms.6.4,20,0..O.apm=dmv200如06008001000oo0>pa_Q.E<200400600知。1000Time/ms误差的平方和£=0.8019睢纳匪波羿系故sipn-i—»一-言<期堕蛔(Rickert)1。80.60.40.20-IIIL—,一[M'+'睢纳匪波羿系故sipn-i—»一-言<期堕蛔(Rickert)1。80.60.40.20-IIIL—,一[M'+'200JOO6008001000050gm_o.E<2D04006003DQ1QOOTime/msYinJ实际精由20200400600300WOOTime/ms名m=dul<误差的平方和£=99573.期望输出为延迟尖脉冲:崎人了做20Q4口06003001000Timem■纳止■波牌系数505o.o.■QDpnl-a.E<期堕鹿出■:Rieke「了-波)2D0400600崎人了做20Q4口06003001000Timem■纳止■波牌系数505o.o.■QDpnl-a.E<期堕鹿出■:Rieke「了-波)2D0400600Time/ms实际■出3001000ooDo224--H瑁b_d.e<o6-o知000040020>p3=-dE<o.e060.402□10002004006008001000Time/ms误差的平方和£=0.0201物入子波o5o.■sipn-i—,--dE<误差的平方和£=0.0201物入子波o5o.■sipn-i—,--dE<期期出(Rick改子波〉1{J.30.60..40.20200400600800Timeyrns蜡纳渡波掘系菽5010002D04006003DQ1000Time/ms实际输出200如。600凯0oo-5mpnwEV.8.64.2DO.O.QO.0>p