地震子波反褶積程序

1、精選優(yōu)質(zhì)文檔-傾情為你奉上 盲信號(hào)實(shí)驗(yàn)報(bào)告 盲解卷積算法姓 名： 丁憲成 系 別： 電信學(xué)院 專 業(yè)： 電磁場(chǎng)與微波 學(xué) 號(hào)： 指導(dǎo)教師： 陳文超 2011年07月13日盲解卷積算法 1. 原理：幾個(gè)重要概念：1.1 褶積模型假設(shè)：假設(shè)1：地層是由具有常速的水平層組成；假設(shè)2：震源產(chǎn)生一個(gè)平面壓縮波(P波)，法向入射到層邊界上，在這種情況下，不產(chǎn)生剪切波(S波)；假設(shè)3：震源波形在地下傳播過程中不變，即它是穩(wěn)定的；假設(shè)4：噪音成分是零；假設(shè)5：震源波形是已知的；假設(shè)6：反射系數(shù)序列是一個(gè)隨機(jī)過程。這意味著地震記錄具有地；震子波的特征，即它們的自相關(guān)和振幅譜是相似的；假設(shè)7：地震子波是最小相位的

2、，因此，它有一個(gè)最小相位的逆。1.2 反濾波 如果定義濾波算子為f (t),則f (t) 與已知地震記錄x(t)的褶積得到一個(gè)對(duì)地層脈沖響應(yīng)e(t)的估計(jì)e(t) = f (t) x(t)； （1）x(t) = w(t)* f (t)* x(t)； （2）(t) = w(t)* f (t)； （3） （4）用流程圖表示為：1.3 震源反子波計(jì)算震源反子波在數(shù)學(xué)上是利用z變換來實(shí)現(xiàn)的。例如，假設(shè)基本子波為兩點(diǎn)時(shí)間序列（1，0.5)； （5） （6）的系數(shù)代表逆濾波算子f(t)有關(guān)的時(shí)間序列?？梢钥闯鏊袩o限多個(gè)系數(shù)，然而它們遞減的很快。如同任何濾波過程一樣，實(shí)際應(yīng)用的算子都是被截?cái)嗟摹?.4 最

3、小平方反濾波當(dāng)輸入子波良好，其z變換的逆可以用一收斂序列表示，則上面所描述的反濾波將得到一個(gè)很好的近似于尖脈沖的輸出.將下面問題列出方程式：給定輸入子波(10.5)，尋找一個(gè)二項(xiàng)濾波器(a,b)，使實(shí)際輸出與期望輸出(1,0,0)之間的誤差在最小二乘的意義上最小。將濾波器(a,b)與輸入子波(1,-0.5)褶積以計(jì)算實(shí)際輸出。誤差的累積能量L定義為實(shí)際的與期望的輸出的系數(shù)之差的平方和： （7）目的是尋找系數(shù)(a,b)使L最小，這要求變量L隨系數(shù)(a,b)而變并使之為零對(duì)上式進(jìn)行簡(jiǎn)化，取L對(duì)a和b的偏導(dǎo)數(shù)，并使其為零，得到： （8）有兩個(gè)方程和兩個(gè)未知數(shù)即濾波器系數(shù)(a,b)，可變成下列普通矩陣

4、形式： （9）求解濾波系數(shù)，我們得到(a,b)為(0.95,0.38)。1.5 最佳維納濾波器回顧研究反濾波和最小平方濾波時(shí)討論的期望輸出，即零延遲尖脈沖(1,0,0)。重寫最小平方反濾波的求解方程如下： （10）兩邊同時(shí)除以2得到：上述結(jié)論被維納普遍化以推導(dǎo)出將輸入轉(zhuǎn)換為任意期望輸出的濾波器(Robinson和Treltel，1980)。一個(gè)濾波器長(zhǎng)度為n的矩陣方程的普遍形式是：（11）式中、和(i0,l,2,3,n-1)分別為輸入子波自相關(guān)、維納濾波系數(shù)和期望輸出與輸入子波的互相關(guān)。其計(jì)算過程的流程圖可以用下面圖形表示為：與最小平方濾波的聯(lián)系：最佳維納濾波器是最佳的，是指它的實(shí)際輸出與期望

5、輸出之間的最小平方誤差最小。當(dāng)期望輸出是零延遲尖脈沖(1,0,0,0)時(shí)，維納濾波與最小平方濾波相同，即后者是前者的特例。1.6 脈沖反褶積第一類期望輸出(零延遲尖脈沖)的處理稱為脈沖反褶積。期望脈沖(1,0,0,0)與輸入子波的互相關(guān)得到序列(1,0,0,0)。正則方程的一般形式改寫成特殊形式：(12)方程已用作了歸一化。在前面討論過的最小平方反濾波器有與矩陣方程相同的形式。因此，脈沖反褶積在數(shù)學(xué)上是與最小平方反濾波相同的。而在實(shí)際上兩類濾波器有一些區(qū)別。在脈沖反褶積(統(tǒng)計(jì)反褶積)的情況，方程在左邊的自相關(guān)矩陣由輸入地震記錄(假設(shè)6)求得，而在最小平方反濾波(確定性反褶積)的情況則直接由已知

6、的震源子波求得。結(jié)論：如果輸入子波不是最小相位，則脈沖反褶積不能將它轉(zhuǎn)化為完全的零延遲尖脈沖。雖然振幅譜實(shí)質(zhì)上是平的，但輸出的相位譜不是最小相位的。而且脈沖反褶積算子是輸入子波的最小相位對(duì)應(yīng)的子波的逆，這個(gè)子波可以是也可以不是最小相位的。1.7 整形子波濾波器處理脈沖反褶積將子波(-0.5,1)壓縮為零延遲尖脈沖(1,0,0)時(shí)有些困難。從能量分布來說，這個(gè)輸入子波更接近于延遲尖脈沖，如(0,1,0)，而不是零延遲尖脈沖(1,0,0)。或者，將子波(-0.5,1)轉(zhuǎn)換為延遲尖脈沖的濾波器會(huì)比將它整形為零延遲尖脈沖的濾波器產(chǎn)生更少的誤差。根據(jù)最佳維納濾波器的流程圖重新設(shè)計(jì)和應(yīng)用一個(gè)濾波器(期望輸

7、出為(0,1,0)。首先，計(jì)算互相關(guān)。我們知道輸入子波的自相關(guān)。將互相關(guān)和子波自相關(guān)的結(jié)果代入前述矩陣方程，得到：(13)求解濾波器系數(shù)，結(jié)果為。這個(gè)濾波器用于輸入子波，如下表所示。和我們期望的一樣，輸出結(jié)果與最小平方濾波的輸出結(jié)果一樣。可以看出，對(duì)于延遲尖脈沖，實(shí)際輸出與期望輸出的最小平方誤差為0.190；而對(duì)期望輸出的零延遲尖脈沖，最小平方誤差為0.762。這表明，將子波(0.5,1)轉(zhuǎn)換為延遲尖脈沖(0,1,0)比零延遲尖脈沖(1,0,0)誤差小。2. 實(shí)驗(yàn)內(nèi)容地震記錄信號(hào)是地震子波與地層脈沖響應(yīng)相卷積的結(jié)果，如果我們要想得到脈沖信號(hào)必須進(jìn)行盲解卷積，地震子波一般可以近似測(cè)得，我們要做的

8、就是求出地震子波對(duì)應(yīng)的逆信號(hào)，理想情況它們相卷積的結(jié)果為單位脈沖。地震子波有零相位子波，最小相位子波，混合相位子波和最大相位子波，如下圖所示：程序說明：在下面的程序中，我們?nèi)我饨o定一個(gè)模擬地震子波，通過一系列的操作得到想要得到的結(jié)果。算法如下：第一步：產(chǎn)生模擬地震子波信號(hào)，給出預(yù)白百分比；第二步：給定期望輸出，計(jì)算輸入的自相關(guān)序列以及輸入和期望輸出的自相關(guān)序列；第三步：根據(jù)輸入信號(hào)的自相關(guān)序列和預(yù)白百分比產(chǎn)生Toeplitz矩陣，根據(jù)Toeplitz矩陣和互相關(guān)序列求出維納濾波器的濾波系數(shù)；第四步：通過將輸入序列和濾波系數(shù)卷積求的實(shí)際輸出。3. 實(shí)驗(yàn)結(jié)果：3.1 得到任意延遲脈沖信號(hào)：輸入想得

9、到的延遲脈沖的延遲時(shí)間:0（零延遲脈沖）輸入想得到的延遲脈沖的延遲時(shí)間:5任意輸出延遲時(shí)間，可以得到任意延遲的脈沖。3.2 得到任意提前的輸入序列運(yùn)行結(jié)果：輸入想得到的提前了的時(shí)間:2輸入想得到的提前了的時(shí)間:10可見，提前的時(shí)間越長(zhǎng)，得到的結(jié)果誤差越來越大。3.3 得到0相位子波運(yùn)行結(jié)果：從對(duì)比圖也可以看出在輸入序列能量集中的地方得到的實(shí)際輸出和期望輸出誤差比較小，其它區(qū)域誤差較大。3.4 得到任意期望輸出：從對(duì)比圖也可以看出，在輸入子波能量集中的地方得到的實(shí)際結(jié)果與期望輸出很近似，在輸入子波幅度很小的區(qū)域得到的實(shí)際輸出結(jié)果誤差很大。3.5 得到均方誤差最小的實(shí)際輸出運(yùn)行結(jié)果：j =13;m

10、in1=0.0015.4討論地震子波的相位對(duì)處理結(jié)果的影響地震數(shù)據(jù)處理的目地是將野外采集的地震記錄用處理模塊進(jìn)行處理后得到成像好，分辨率高的地震剖面，地震記錄可描述為地震子波與地層脈沖響應(yīng)或地下反射系數(shù)的褶積。就某種意義是那個(gè)講，地震數(shù)據(jù)處理實(shí)際上就是一個(gè)對(duì)地震子波不斷改造的過程。地震子波經(jīng)過傅立葉變換之后可以得到振幅譜和相位譜，因此在地震記錄中可通過拓寬地震子波的振幅譜來提高地震剖面的分辨率，也可以通過改變地震子波的相位譜來達(dá)到提高分辨率的目的。分辨率高的地震子波在頻率域上表現(xiàn)為振幅譜盡可能的拓寬。而具有同樣振幅譜的地震淄博中，最小相位的地震子波分辨率最高。實(shí)際地震數(shù)據(jù)處理中對(duì)地震子波相位特

11、性有要求和影響的處理模塊主要包括可控震源地震子波最小化，預(yù)濾波，疊前反褶積和疊后地震子波零相位化。我們知道地震數(shù)據(jù)處理的母的是提高地震記錄的信噪比和分辨率，得到高信噪比，高分辨率的地震剖面。在頻率域中表現(xiàn)為使得振幅頻帶盡可能寬且剩余子波是零相位的，在處理中有些模塊對(duì)地震數(shù)據(jù)的子波相位要求有假設(shè)前提條件。一般情況下，疊前反褶積輸入的地震數(shù)據(jù)要求子波是最小相位的，若是零相位的需要進(jìn)行最小相位化，根據(jù)相同的振幅譜零相位的子波分辨率最高的原則，需要對(duì)疊后數(shù)據(jù)的子波進(jìn)行零相位化處理。5 附錄：5.1 得到想要的任意延遲脈沖t=0:1:30;x=sin(pi*(t-2)/6.4).*exp(-0.12*a

12、bs(t-2); %產(chǎn)生模擬子波subplot(211);plot(t,x);m=length(x);p=0.005; %預(yù)白百分比 i=input(輸入想得到的延遲脈沖的延遲時(shí)間:);y=zeros(1,i),1,zeros(1,2*m-2-i); %期望輸出n=length(y);r_xx=xcorr(x); %計(jì)算子波的自相關(guān)A=fliplr(r_xx(1:m);R=toeplitz(1+p/100)*A(1),A(2:length(A); %產(chǎn)生Toeplitz矩陣r_xy=xcorr(x,y); %輸入和期望輸出的互相關(guān)G1=fliplr(r_xy(1:n);G2=G1(1:m);h

13、=inv(R)*G2; %h為維納濾波器系數(shù)z=conv(x,h); %計(jì)算實(shí)際輸出信號(hào)subplot(212);plot(z) %繪制滿足最小均方誤差的實(shí)際輸出結(jié)果5.2 得到時(shí)間提前了的輸入序列說明：假設(shè)輸入是（1，2，3，4，5，6），提前時(shí)間為2個(gè)單位，那么得到的期望輸出序列為（3，4，5，6，0，0，），所以上面的程序只要修改一下期望輸出即可，將下面兩行程序：i=input(輸入想得到的延遲脈沖的延遲時(shí)間:);y=zeros(1,i),1,zeros(1,2*m-2-i); %期望輸出修改為：i=input(輸入想得到的提前了的時(shí)間:);y=x(i+1:m),zeros(1,m+i-

14、1); %期望輸出5.3 得到零相位子波零相位子波：所謂零相位子波，就是子波序列進(jìn)行傅立葉變換得到的結(jié)果全是正實(shí)數(shù)，因?yàn)檎龑?shí)數(shù)的相位為零，我們可以給定一個(gè)實(shí)偶序列，其傅立葉變換是正偶實(shí)數(shù)，則它肯定是零相位的，例如余弦序列，其傅立葉變換的結(jié)果就是正實(shí)數(shù)，所以是零相位的，那么程序就要修改為：t=0:1:30;x=sin(pi*(t-1)/6.4).*exp(-0.12*abs(t-1); %產(chǎn)生模擬子波subplot(311);plot(t,x);m=length(x);p=0.005; %預(yù)白百分比 t1=-30:1:30;y=cos(pi*t1/30); %期望輸出subplot(312)pl

15、ot(y);n=length(y);r_xx=xcorr(x); %計(jì)算子波的自相關(guān)A=fliplr(r_xx(1:m);R=toeplitz(1+p/100)*A(1),A(2:length(A); %產(chǎn)生Toeplitz矩陣r_xy=xcorr(x,y); %輸入和期望輸出的互相關(guān)G1=fliplr(r_xy(1:n);G2=G1(1:m);h=inv(R)*G2; %h為維納濾波器系數(shù)z=conv(x,h); %計(jì)算實(shí)際輸出信號(hào)subplot(313);plot(z) %繪制滿足最小均方誤差的實(shí)際輸出5.4 得到任意期望的輸出期望輸出可以用隨機(jī)數(shù)產(chǎn)生，程序改為：t=0:1:30;x=si

16、n(pi*(t-2)/6.4).*exp(-0.12*abs(t-2); %產(chǎn)生模擬子波subplot(311);plot(t,x);m=length(x);p=0.005; %預(yù)白百分比 y=rand(1,2*m-1); %期望輸出subplot(312);plot(y);n=length(y);r_xx=xcorr(x); %計(jì)算子波的自相關(guān)A=fliplr(r_xx(1:m);R=toeplitz(1+p/100)*A(1),A(2:length(A); %產(chǎn)生Toeplitz矩陣r_xy=xcorr(x,y); %輸入和期望輸出的互相關(guān)G1=fliplr(r_xy(1:n);G2=G1

17、(1:m);h=inv(R)*G2; %h為維納濾波器系數(shù)z=conv(x,h); %計(jì)算實(shí)際輸出信號(hào)subplot(313);plot(z) %繪制滿足最小均方誤差的實(shí)際輸出結(jié)果5.5 得到均方誤差最小的實(shí)際輸出上面幾個(gè)程序都是在期望輸出給定的情況下求的的，但實(shí)際上如果輸入子波為最大相位，期望輸出又是無延遲的脈沖，得到的實(shí)際結(jié)果往往和期望輸出誤差相差很大，所以如果要得到較好的實(shí)際輸出，我們可以通過對(duì)期望輸出做一定的延遲來得到最優(yōu)化的結(jié)果，這里的最優(yōu)化就是最小均方誤差。可以通過改動(dòng)一下上面的程序，增加判別均方誤差最小的步驟得到：t=0:1:30;x=sin(pi*(t-1)/6.4).*exp(-0.12*abs(t-1); %產(chǎn)生模擬子波subplot(211);plot(t,x);m=length(x);min1=10000;p=0.005; %預(yù)白百分比for i=0:m-1 y=zeros(1,i),1,zeros(1,2*m-2-i); %期望輸出n=length(y);r_xx=xcorr(x); %計(jì)算子波序列自相關(guān)A=fliplr(r_xx(1:m);R=toeplitz