《盲解卷积和详细程序》由会员分享，可在线阅读，更多相关《盲解卷积和详细程序（16页珍藏版）》请在金锄头文库上搜索。
1、实验3:基于最佳维纳滤波器的盲解卷积算法一. 算法原理：1. 概论：反褶积是通过压缩地震记录中的基本地震子波，压制交混回响和短周期多次 波，从而提高时间分辨率，再现地下地层的反射系数。反褶积通常应用于叠前资 料，也可广泛用于叠后资料。理想的反褶积应该压缩子波并消除多次波，在地震地道内只留下地层反射系 数。子波压缩可以通过将反滤波器作为反褶积算子来实现，它与地震子波做褶积 时，反滤波器可以将地震子波转变成尖脉冲。当应用于地震合成记录时，反滤波 输出应为地层脉冲响应，精确的反滤波器设计可用最小平方模型来实现。反褶积处理的基本假设是震源子波为最小相位。2. 褶积模型：假设1：地层是由具有常速的水平层
2、组成；假设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：反射系数序列是一个随机过程。这意味着地震记录具有地震子波的 特征，即它
3、们的自相关和振幅谱是相似的；假设7：地震子波是最小相位的，因此，它有一个最小相位的逆。3. 最佳维纳滤波器：维纳滤波器是以最小平方误差为准则的，即要使下式最小：设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或 传递函数的表达式，其实质就是解维纳一霍夫(Wiener-Hopf)方程。滤波器的维纳一霍夫方程如下：莒。升. rn_2口 i幻J_如一式中，r , a和g ( i = 0,1,2,.n-1 )分别为输入子波的自相关、维纳滤波系 i i i数和期望输出与输入子波的互相关。下图-1为维纳滤波器的设计和应用流程图:图-1为维纳滤波器的设计和应用流程图确定维纳滤波器的系数需要求
4、解维纳一霍夫方程，由方程可以看到自相关矩 阵是对称的。这个特殊矩阵称作Toeplitz矩阵，可用莱文逊递归法求解。最佳维纳滤波器(a ,a ,a，,a )是最佳的，是指它的实际输出与期望输 0 1 2n-1出之间的最小平方误差最小。当期望输出是零延迟尖脉冲(1,0,0，.，0)时， 维纳滤波器与最小平方滤波器相同，即后者是前者的特例。维纳滤波器可以考虑任一种期望输出而不仅限于零延迟尖脉冲。期望输出可以有5类选择：类型1：零延迟尖脉冲；类型2：任一延迟尖脉冲；类型3：时间提前了的输入序列；类型4：零相位子波；类型5：任意期望波形。二. 常用子波：地震资料处理中常用的子波有以下几种：(1) Ric
5、ker 子波：时域表达式：w(t) = (1-2兀 2f 212)e-兀2fM12(3-4)M（3-5）频域表达式：W（f） = 土 f e-f其中，八为子波的主频。下图为主频fM =50Hz的Ricker的时域和频域波形。Fourier Spectrum Df Ricker图-2 Ricker子波时域/频域波形（2）Berlage子波：(3-6)时域表达式为：w(t) = AH(t)tne-b cos(2兀 f t + )下图为主频为30Hz的Berlage子波的时域和频域波形:Berlage干波觥闾油形Berime子波时嫩波影 10.2T|n&/m&Fr&quencH2图-3 Berlag
6、e子波时域/频域波形（3）一种常用的模拟子波：时域表达式：w（t） = A sin（旦）e-0.12t-t0（ 3-7）6.4该子波对应的最小相位、最大相位、零相位和混合相位子波如下4图所示。(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
7、(输入子波);grid on;y=zeros(1,N);y(1)=1;% 期望输出subplot(2,2,2);plot(y);title(期望输出);grid on;%求维纳滤波器的系数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(维纳滤波器系数);grid on;yy=conv(
8、x,h); % 实际输出yy=yy(1: N);subplot(2,2,4);plot(yy);title(实际输出);grid on;error=norm(yy-y)八2%误差的累积能量四. 结果分析：（一）五种不同的期望输出：用以上的程序验证算法的正确性，以下五幅图为五种不同类型的期望输出下 得到的维纳滤波器系数和实际的输出。（输入子波为式（3-7）中的最小相位子波， 如图4（a）所示）下图5的期望输出为零延迟尖脉冲，实际输出与期望输出的误差的平方和（能量）L = 2.77 x 10-5.o o o O2 2 4 - cpnwdE-60-0.5200 400 600 800 1000Tim
9、e/ms维纳波波器系数200 400 600 BOO WOO Time/ms 实际输出输入子波（最小相位）期里输出（零延退单位脉冲）5 oQ省 m_-dmprl三dE200 400 600 800 1000图5-零延迟尖脉冲下图6的期望输出为延迟100ms的尖脉冲，实际输出与期望输出的误差的平方和（能量）L = 4.93 x10-5.输入子披（最小相惊）200 400 600 800 1000期望输出（延退卞00单位聊1 I 0 80 60.40 2200 400 600 BOO 1000Tim Mrs实际馆出0.80 60 40.2Q J. 200 400 600 800 1000Tim M
10、rs图6-延迟100ms尖脉冲下图7的期望输出为时间提前了的输入子波，实际输出与期望输出的误差的输入了波（最小捌位一）靳咀输出时间提值了的输入序刊）200 JOO 6008001000实味输出平方和（能量）L = 1.3 x 10-3.图7-时间提前了的输入子波下图8的期望输出为零相位输入子波，实际输出与期望输出的误差的平方和（能量）L = 3.85x 10-4.锥熟源赧羿系数推入子波（最小相仲）0.S0.60.40 20-02-042D04006008001000Time/ms200 JOO 6008001000O 1.23 4tlpnwdE图8-零相位子波下图9的期望输出为主频为10Hz的
11、Ricker子波，实际输出与期望输出的误差的平方和（能量）L= 2.3x10-4.输入子波（最小捌位）5 O osspnMELE-0.5 w1:20D4口06003001000Timem期营输出(RickerT-)2004006008001000Time/ms实际输出;5O tupm=dulvI-0.05200 JOO 60080010006 6 4 2 0 0O.O.O.命 pm_o.E-020.4 2004006008001000Time/ms服汕凋波普系徂图9-Ricker子波从以上五种期望输出的结果可以看出，期望输出与实际输出的误差很小，从 而说明最佳维纳滤波器的正确性。（二）子波的相
12、位对结果的影响:1 .期望输出为Ricker子波：期营输出(Ri如子波200400600800Time/ms 实际输出-So tuTFm云 mv输人子波20Q4口06003001000Timeyms 推纳源波普系数5 o osQDpnMeLE10001000400600800Time/ms2008.6.4.2Q 0.O.0.O.命 pnwE-0 20.4-0 05200 JOO 6008001000(a)最小相位误差的平方和乙=2.3 x 10-4期堕鹿出(HidwrTfe)200400 前 0800Time/ms 实际输出5 o aupm=dmv1000Time/ms场入子波20D4口060
13、03001000Time/mso 5 -O. 如 pn-i-eLE推纳凋波羿系敖200 JOO 6008001000输人子波SOS .o.一 spnl-a.E200400 前 0800Time/ms 实环箱出期箜输出(RickerT-Wl)2004Q0600800Timeynns赧纳胡波掘系数3100010002 1 O-1-2由 pn* 一一 dmpnwE.5o.WOO 400600800误差的平方和 = 0.0669物氏子波Q.51000200JU400600800Time/ms1000(c)混合相位220040060030。1000Timeyrns罪.纳旅波牌系疝o 5 o.-aipnl-dEpm_D.E误差的平方和 = 26.12842004006003001000Time/ms(d)混合相位3-0.5物人子波5 0 o.QDpnl-ELE4006
