最小平方反褶積實(shí)驗(yàn)報(bào)告

1、精選優(yōu)質(zhì)文檔-傾情為你奉上實(shí)驗(yàn)四 最小平方反褶積實(shí)驗(yàn)題目：已知兩個(gè)地震模型，反射系數(shù)已給出。Re1為厚層模型，Re2為薄層模型，時(shí)間域采樣間隔dt=1ms.要求用50Hz雷克子波先分別生成兩個(gè)模型的合成地震道，然后根據(jù)最小平方反褶積的原理編寫(xiě)脈沖反褶積程序，并將反褶積后的結(jié)果與已知模型比較?？梢詫?duì)模型添加不同水平的隨機(jī)噪聲，檢驗(yàn)其對(duì)反褶積算法的影響。實(shí)驗(yàn)內(nèi)容：把延續(xù)幾十至l00ms的地震子波壓縮成原來(lái)的震源脈沖形式，地震記錄變?yōu)榉从撤瓷湎禂?shù)序列的窄脈沖組合，這就是反濾波所要完成的工作。反褶積的目的就是為了把地震子波壓縮成尖脈沖，使實(shí)際的地震記錄變成反射系數(shù)序列。假設(shè)地震記錄為 （11）其中為有

2、效信號(hào)，為干擾波。首先假設(shè)不存在干擾波，即： （12）對(duì)兩邊求傅氏變換，則得到頻率域的地震記錄表示式： （13）式中，、和分別為地震頻譜、子波頻譜和反射系數(shù)的頻譜。顯然： （14）如果令： （15）則有： （16）再對(duì)（16）式做反傅氏變換至?xí)r間域，就可得到： （17）式中，為的時(shí)間函數(shù)。根據(jù)（17）式知： （18）因?yàn)闉榈卣鹱硬ǎ椭g又存在著頻譜互為倒數(shù)的關(guān)系（即），由此可知，如已知地震子波，利用數(shù)學(xué)方法求出，再利用（17）式讓反子波與地震記錄做褶積，就可以求出反射系數(shù)序列，即 （19）經(jīng)過(guò)這樣的處理，就可以達(dá)到把地震子波壓縮成尖脈沖，從而達(dá)到提高地震記錄縱向分辨能力的目的。脈沖反褶積的

3、基本思想在于設(shè)計(jì)一個(gè)濾波算子，用它把已知的輸入信號(hào)轉(zhuǎn)換為與給定的期望輸出信號(hào)在最小平方誤差的意義下是最佳接近的輸出。若將地震子波作為反濾波的輸入，期望輸出則為尖脈沖。若設(shè)計(jì)另一濾波器輸入信號(hào)是某濾波器的輸出，而期望輸出是該濾波器的輸入，則按此思想求得的濾波因子即稱(chēng)為脈沖反濾波因子，用它進(jìn)行的濾波就是脈沖反濾波，即脈沖反褶積。先假設(shè)期望輸出為窄脈沖，在子波已知的情況下，設(shè)待求的反濾波因子起始時(shí)刻為，延續(xù)長(zhǎng)度為。即當(dāng)已知輸入地震子波時(shí)，實(shí)際輸出為實(shí)際輸出與期望輸出的誤差平方和為 （110）要使Q為最小，數(shù)學(xué)上就是求Q的極值問(wèn)題，即求滿(mǎn)足 （111）的濾波因子。為地震子波的自相關(guān)函數(shù)，而為地震子波與

4、期望輸出的互相關(guān)函數(shù)，故（111）式可寫(xiě)為 （112）此方程系數(shù)矩陣即為拖布利茲矩陣。若期望輸出是脈沖，則互相關(guān)為 （113）基本方程（112）變?yōu)?（114）一般情況下，地震子波為未知的，為在未知子波的情況下求出反濾波因子，必須對(duì)地震子波及反射系數(shù)序列加上一定的假設(shè)條件，他們包括：A.假設(shè)反射系數(shù)序列是隨機(jī)的白噪序列，即其自相關(guān)為 （115）B.假設(shè)地震子波是最小相位的。根據(jù)假設(shè)A，地震子波的自相關(guān)可以用地震記錄的自相關(guān)代替。根據(jù)假設(shè)B，可知地震子波的Z變換的零點(diǎn)全部在單位圓外，也即反濾波因子的Z變換的分母多項(xiàng)式的零點(diǎn)全在單位圓外，故是穩(wěn)定的、物理可實(shí)現(xiàn)的。因此，自由項(xiàng)變?yōu)?。又因必為物理可?shí)

5、現(xiàn)的，故，。令，則基本方程變?yōu)?（116） 這就是脈沖反褶積的基本方程，其系數(shù)矩陣中各元素可直接由地震記錄求得。求出的反濾波因子僅與相差常數(shù)倍，不影響壓縮子波、提高分辨率的反濾波作用。當(dāng)求取了反褶積因子后，令其與地震記錄進(jìn)行褶積運(yùn)算，即，則即為經(jīng)過(guò)脈沖反褶積之后輸出地新的地震記錄。源程序：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:studyr

6、e1-2.txt');R2=load('D:studyre2-2.txt');A1=conv(w,R1);A2=conv(w,R2);figure(2),plot(A1,'k'),hold onfigure(3)，plot(A2,'r'),hold offy1=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)+len

7、gth(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'),hold onfigure(5),plot(record2,'r'),hold off運(yùn)行程序得圖：地震子波合成地震記錄（厚層模型、薄層模型）反褶積得到的反射系數(shù)與模型比較（厚

8、層、薄層）利用合成地震記錄自相關(guān)構(gòu)造Toeplitz矩陣進(jìn)行反褶積運(yùn)算,得到的反射系數(shù)與模型比較：（厚層、薄層）在模型上添加白噪，使更符合地震記錄自相關(guān)前提條件，進(jìn)行反褶積運(yùn)算自上而下：添加噪聲(0.001*randn(size(R)、(0.01*randn(size(R)、(0.1*randn(size(R)噪聲較?。?.001）時(shí)，地震記錄自相關(guān)得到的反射系數(shù)含有較多的毛刺，后面的部分甚至掩蓋了模型中的反射系數(shù)；中間的圖（0.1），后面的幾個(gè)主要反射系數(shù)在計(jì)算結(jié)果上的體現(xiàn)，幅度要比毛刺大，辨別程度比上圖高；而當(dāng)白噪添加較多時(shí)，實(shí)際上全是白噪聲（掩蓋了R1）計(jì)算得到的反射系數(shù)與模型相差不大。地震記錄添加噪聲對(duì)反褶積結(jié)果的影響：噪聲程度較小時(shí)，就已經(jīng)影響到了反射系數(shù)的計(jì)算，但大體還能夠反應(yīng)部分反射系數(shù)的信息；而噪聲程度較大時(shí)，計(jì)算得到的反射系數(shù)已經(jīng)完全不能反映模型的情況。結(jié)論分析：反射系數(shù)與白噪相差較大時(shí)，用地震記錄自相關(guān)代替子波自相關(guān)進(jìn)行反褶積運(yùn)算，得到的反射系數(shù)誤差較大。厚層模型更明顯，