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

2、，因此，它有一個最小相位的逆。1.2 反濾波 如果定義濾波算子為f (t),則f (t) 與已知地震記錄x(t)的褶積得到一個對地層脈沖響應e(t)的估計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 震源反子波計算震源反子波在數學上是利用z變換來實現的。例如，假設基本子波為兩點時間序列（1，0.5)； （5） （6）的系數代表逆濾波算子f(t)有關的時間序列。可以看出它有無限多個系數，然而它們遞減的很快。如同任何濾波過程一樣，實際應用的算子都是被截斷的。1.4 最

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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