SEGY課程設(shè)計

1、成 績評閱人成都理工大學(xué)工程技術(shù)學(xué)院地震資料分析課程設(shè)計結(jié)業(yè)報告任務(wù)號任務(wù)名稱提交資料1SEGY格式文件讀寫1文本報告、程序（電子版）2瞬時屬性分析技術(shù)文本報告、程序（電子版）3SEGY格式數(shù)據(jù)操作2程序（電子版）4計算地震道振幅譜程序（電子版）班 級：2012級勘查技術(shù)與工程2班姓 名：李南吉學(xué) 號：201220911222提交日期：2015年5月11日 JOB 1 SEGY格式文件讀寫1一、 實驗內(nèi)容1. 了解常見的磁帶記錄格式2. 編程實現(xiàn)SEG-Y格式讀寫，并提取inline為（110+學(xué)號）的橫剖面和cross-line（400-學(xué)號）的縱剖面。二、 實驗相關(guān)軟件和數(shù)據(jù)說明1. Se

2、gy文件分析工具2. Fimage地震資料顯示程序3. stack.sgy 三維地震資料(in-line范圍105-164,cross-line范圍：251-450，time范圍500-1500ms)三、 實驗內(nèi)容和步驟1. SEGY格式地震數(shù)據(jù)讀寫地震勘探主要包括三個階段，分別為野外地震數(shù)據(jù)采集、室內(nèi)地震資料處理和地震資料解釋。且現(xiàn)代資料處理和解釋都是在計算機上進(jìn)行的，因此數(shù)據(jù)必須通過相應(yīng)的介質(zhì)交換和存儲，常見的存儲介質(zhì)包括磁帶和盤類（硬盤、光盤、U盤等）。由于野外工作環(huán)境和早期的工藝限制，地震隊野外施工采集到的數(shù)據(jù)都是存儲在磁帶上的，近些年由于工藝進(jìn)步，室內(nèi)處理解釋中數(shù)據(jù)存儲和交換介質(zhì)基本

3、都是通過硬盤。不管是磁帶介質(zhì)還是盤類介質(zhì)，其數(shù)據(jù)都是以一定的二進(jìn)制格式存儲的。所謂格式，指的是地震數(shù)據(jù)以及各種信息在文件內(nèi)部的存放方式及順序。常見的數(shù)據(jù)格式有SEGY、SEGD、SEG2等，這些格式都是由勘探地球物理學(xué)家協(xié)會（Society of Exploration Geophysicists，SEG）制定的（全世界統(tǒng)一的地震數(shù)據(jù)格式）。SEGD和EEG2是野外（時序）記錄格式，其中SEGD格式是野外使用最多的格式，SEG2格式僅適用于PC工作環(huán)境，目前只有工程地震勘探中在使用。SEGY是室內(nèi)(道序)記錄格式（解編后的炮集記錄，CMP道集，疊加和偏移剖面均為此格式），也是室內(nèi)處理解釋作業(yè)最

4、為常用的標(biāo)準(zhǔn)格式，目前各種商業(yè)處理軟件都提供SEGY格式數(shù)據(jù)加載功能，本講義也是也是以SEGY格式為標(biāo)準(zhǔn)進(jìn)行講解的。圖1 SEGY文件格式1) SEGY數(shù)據(jù)格式SEGY數(shù)據(jù)格式由三部分組成（如圖1所示）。第一部分是3200字節(jié)的EBCDIC卡片頭（編碼方式不同于ASCII）（圖2），稱其為說明信息（text header），包括40個對磁帶進(jìn)行描述的文本數(shù)據(jù)（例如，每行包括08個字符的40行文本）。這部分一般包括工區(qū)信息、數(shù)據(jù)采集者、處理者和解釋者等信息。第二部分是400字節(jié)的二進(jìn)制卷頭（Reel Header 或 binary header）（表1），描述了該磁帶卷上內(nèi)容的有關(guān)信息。SEGY

5、格式的第三部分是實際地震道數(shù)據(jù)。每道含有240字節(jié)的道頭（Trace Header）（表2），數(shù)據(jù)記錄在道頭后面，可能使用4種32位IBM浮點型格式中的一種。（注意，IBM格式與現(xiàn)代普通IBM PC機上所通用的IEEE格式不同，在PC機進(jìn)行數(shù)據(jù)處理前需要先進(jìn)行數(shù)據(jù)轉(zhuǎn)換）。2）SEGY格式編程思路3200字節(jié)說明部分，主要記錄采集處理信息，一般處理地震數(shù)據(jù)不用處理這部分的內(nèi)容。這部分內(nèi)容只需要讀出來再寫入新的文件中去即可，編程時我們只需要定義一個長度為3200的char類型數(shù)組即可。卷頭中包含了采樣間隔和采樣點數(shù)的重要信息，編程時定義為short int類型數(shù)組，就可以得到這些參數(shù)了。數(shù)據(jù)總道數(shù)

6、判斷，用fseek函數(shù)分別指向文件起始和末尾，偏移量就為文件大?。ㄗ止?jié)）LEN，通過如下公式可以求出總道數(shù)總道數(shù)=（LEN-3200-400）/（每道采樣點數(shù)*4+240）三維數(shù)據(jù)體比二維數(shù)據(jù)體多了一維，多的這一維代表線號，和二維數(shù)據(jù)體相比較，如果二維數(shù)據(jù)體代表的是一條線的數(shù)據(jù)，那三維數(shù)據(jù)體就是代表一個工作面的數(shù)據(jù)。由此我們也知道通常情況下，三維數(shù)據(jù)體的數(shù)據(jù)通常情況下會比二維數(shù)據(jù)體大很多。三維數(shù)據(jù)體里面的數(shù)據(jù)格式如下（1） 卷頭部分，3200字節(jié)+400字節(jié)（2） 第一線(inline 1)數(shù)據(jù)：第一道（cross-line 1）道頭部分，第一道數(shù)據(jù)部分 第二道（cross-line 2）道頭

7、部分，第二道數(shù)據(jù)部分 ， （3） 第二線(inline 2)數(shù)據(jù)：第一道道頭部分，第一道數(shù)據(jù)部分 第二道道頭部分，第二道數(shù)據(jù)部分 ， （4） ： ， （5） 最后一線數(shù)據(jù)：第一道道頭部分，第一道數(shù)據(jù)部分 第二道道頭部分，第二道數(shù)據(jù)部分 ， 四、程序?qū)崿F(xiàn)#include<stdio.h>void main()char SM3200; /定義說明short int reelheader200; /定義卷頭 float traceheader60; /定義道頭int i,LEN,N,n; / N為總道數(shù) LEN為文件大小（字節(jié)）FILE *fp,*fp1,*fp2;fp=fopen(&q

8、uot;112","rb"); /打開原文件后轉(zhuǎn)換成pc格式fp1=fopen("橫剖面","wb"); /橫剖面fp2=fopen("縱剖面","wb"); /縱剖面for(i=0;i<3200;i+)fscanf(fp,"%c",&SMi); /輸入說明fread(reelheader,400,1,fp); /輸入卷頭for(i=0;i<200;i+)n=reelheader10;fseek(fp,0L,2);LEN=ftell(fp);N=

9、(LEN-3200-400)/(n*4+240); printf("%dn",N); /總道數(shù) fseek(fp,16609200L,0);for(i=0;i<200;i+)fread(traceheader,240,1,fp);fwrite(traceheader,240,1,fp1);fread(traceheader,2004,1,fp);fwrite(traceheader,2004,1,fp1); /橫剖面fseek(fp,387324L,0);for(i=0;i<60;i+)fread(traceheader,240,1,fp);fwrite(tra

10、ceheader,240,1,fp2);fread(traceheader,2004,1,fp);fwrite(traceheader,2004,1,fp2);fseek(fp,446556L,1); /縱剖面 JOB 2 SEGY數(shù)據(jù)格式操作2一、實驗內(nèi)容：1. 地震數(shù)據(jù)再采樣，將原始地震數(shù)據(jù)進(jìn)行兩倍抽稀2. 讀取IBM浮點數(shù)格式數(shù)據(jù)2. 嘗試實現(xiàn)IEEE和IBM浮點數(shù)格式轉(zhuǎn)換二、數(shù)據(jù)說明stackieee.sgy IEEE格式數(shù)據(jù)（微機格式）stackieeeibm.sgy IBM格式數(shù)據(jù)（工作站格式）三、實驗內(nèi)容和步驟1. 地震數(shù)據(jù)再采樣地震資料處理中，經(jīng)常需要對地震數(shù)據(jù)進(jìn)行再采樣，一般

11、進(jìn)行地震數(shù)據(jù)的抽稀操作，這樣可以降低數(shù)據(jù)大小，加快處理速度。抽稀過程中需要注意抽稀后信號采樣率應(yīng)該滿足采樣定理。實現(xiàn)方法：1） 改變卷頭中采樣間隔（17 -18字節(jié)，單位：微妙）和每道采樣點數(shù)（21-22）。2） 將數(shù)據(jù)道進(jìn)行抽稀。2. 讀取IBM浮點數(shù)格式數(shù)據(jù)參考文獻(xiàn)：1）32位IEEE和IBM浮點數(shù)結(jié)構(gòu)及其轉(zhuǎn)換方法2）三維地震數(shù)據(jù)體任意剖面抽取算法及C語言程序?qū)崿F(xiàn)3. IEEE和IBM浮點數(shù)格式轉(zhuǎn)換特定計算機硬件平臺，只能處理與自己硬件平臺相同格式的浮點數(shù)。地震資料處理和解釋通常使用IBM、SUN等大型計算機系統(tǒng)。不同體系架構(gòu)的計算機系統(tǒng)之間的浮點數(shù)交換與共享，對于地震數(shù)據(jù)處理具有重要意義

12、。32位IEEE和IBM浮點數(shù)是SEGY地震數(shù)據(jù)中兩種常用的數(shù)據(jù)格式，因此有必要掌握這兩種格式之間的轉(zhuǎn)換方法。參考文獻(xiàn)：1）32位IEEE和IBM浮點數(shù)結(jié)構(gòu)及其轉(zhuǎn)換方法2）三維地震數(shù)據(jù)體任意剖面抽取算法及C語言程序?qū)崿F(xiàn)四、主要的程序?qū)崿F(xiàn)#include<stdio.h>void main()char SM3200; /定義說明short int reelheader200; /定義卷頭short int traceheader120; /定義道頭int i,j=0,k;float a501,yy251; / a為原采樣點數(shù),yy為更改后采樣點數(shù)FILE *fp,*fp1;fp=fo

13、pen("stackieee.sgy","rb");fp1=fopen("抽稀.sgy","wb"); fseek(fp,0L,0);for(i=0;i<3200;i+)fscanf(fp,"%c",&SMi); fprintf(fp1,"%c",SMi); /輸入說明fseek(fp,3200L,0);fread(reelheader,400,1,fp); /輸入卷頭reelheader8=4000; /17-18字節(jié)-采樣間隔reelheader10=251

14、; /21-22字節(jié)-每道采樣點數(shù) fwrite(reelheader,400,1,fp1); / 寫入卷頭 fseek(fp,3600L,0); for(i=0;i<12000;i+) fread(traceheader,240,1,fp); /讀出道頭 traceheader57=251; / 115-116字節(jié)-每道采樣點數(shù) traceheader58=4000; / 116-117字節(jié)-采樣間隔 fwrite(traceheader,240,1,fp1); fread(a,2004,1,fp); for(k=0,j=0;k<501;k+) if(k%2=0) yyj+=ak

15、; fwrite(yy,1004,1,fp1); /輸出抽稀后數(shù)據(jù) fclose(fp); fclose(fp1); JOB 3 瞬時屬性分析技術(shù)要求：計算地震數(shù)據(jù)體三瞬屬性，并分別提取inline為（120+學(xué)號）的橫剖面和cross-line（375-學(xué)號）的縱剖面的瞬時振幅、瞬時相位、瞬時頻率剖面。快速傅里葉變換（FFT）程序：MCMPFFT.C說明：基于Cooley-Tukey radix-2 DIT算法的FFT20世紀(jì)70年代后期， Taner等人（1979）在全面研究復(fù)數(shù)地震道三瞬參數(shù)分析原理的基礎(chǔ)上，開始利用復(fù)數(shù)地震道分析地震信號的瞬時振幅（Instantaneous Ampli

16、tude）、瞬時相位（Instantaneous Phase）和瞬時頻率（Instantaneous Frequency） ，這就是著名的三瞬參數(shù)分析，用于描述地震道中譜信息隨時間的變化屬性（Taner和Sheriff等，1977）。對實際地震資料的復(fù)數(shù)道分析做出了開創(chuàng)性貢獻(xiàn)。迄今為止，利用復(fù)數(shù)道分析（解析信號）得到的三瞬時參數(shù)又衍生出了許多屬性，因此，基于復(fù)數(shù)道的瞬時屬性分析仍然是地震資料解釋中最重要的屬性之一。1. Hilbert（希爾伯特）變換信號x(t)的希爾伯特變換定義為記為或。顯然是x(t)通過LTI系統(tǒng)的輸出，該系統(tǒng)的頻率響應(yīng)為如圖1(b)所示，它的幅頻特性全為1。圖1(a)至(

17、c)給出了其頻域作用。圖1 希爾伯特變換希爾伯特變換的實質(zhì)是對信號的正、負(fù)頻率部分分別實施相移和。由實信號x(t)與它的希爾伯特變換構(gòu)成的復(fù)信號為稱為x(t)的解析信號(analytic signal)或信號預(yù)包絡(luò)(pre-envelope)，有時也將解析信號稱為信號的希爾伯特變換?？梢钥闯龆街惺穷l域的階躍函數(shù)。因此，如果信號x(t)是確定信號，則解析信號也是確定的，且其頻譜為實部和虛部可看成是解析信號（復(fù)數(shù)道）在兩個互相正交的平面上的投影（如圖2）。圖2 解析信號的直觀理解由傅里葉變換的特性知，實信號頻譜的負(fù)頻率部分完全是冗余的，它可以從頻譜的正頻率部分獲得。希爾伯特變換是在通信理論中引入

18、的（Gabor，1946），在無線通信中，常將實信號頻譜的負(fù)頻率部分去掉，只保留其正頻率部分，這樣頻譜就不再具有共軛對稱性，這種頻譜對應(yīng)的時域信號就是解析信號。已成為通信理論研究的一個重要工具。從上面分析可以看出，解析信號的算法有兩種，分別為時域算法和頻域算法。本課程設(shè)計中我們講解頻域算法。信號x(t)的解析信號z(t)的頻域算法如下1. 計算信號x(t)的快速傅里葉變換（FFT）X()。2. 創(chuàng)建數(shù)組h與x()大小相同，其值為3. 計算ht和X的點積（Z =ht.* X）4. 計算Z 的傅里葉反變換即為x(t)的希爾伯特變換。2. 三瞬參數(shù)計算瞬時屬性定義如下：1） 瞬時振幅：2） 瞬時相位

19、：3） 瞬時頻率：為避免對瞬時相位求導(dǎo)帶來的誤差（因為相位與振幅無關(guān)，微弱的噪聲都會對相位產(chǎn)生嚴(yán)重影響，會帶來偏差，如果用標(biāo)準(zhǔn)公式中的相位對時間求導(dǎo)計算瞬時頻率，會使誤差發(fā)生傳遞和放大）,故利用解析信號本身求瞬時頻率：式中:Imag為復(fù)數(shù)取虛部算子，T為采樣間隔。三瞬屬性剖面示例如下圖所示。圖3 三瞬屬性示例2. 三瞬屬性的地質(zhì)意義通過各種數(shù)學(xué)變換計算出的地震屬性多達(dá)上百種，并不是所有的屬性都有明確的地質(zhì)意義。對于地震屬性的分類，目前并沒有統(tǒng)一的標(biāo)準(zhǔn)，按照不同的劃分方法，同一屬性的歸屬類別也不一樣，此處不詳究某一屬性的歸屬問題，只是簡單介紹一下常用的地震時頻屬性及其它們的物理和地質(zhì)意義。1）

20、瞬時振幅地層巖性的變化或儲層中流體性質(zhì)、孔隙度的改變能夠引起層速度的變化，而這些變化在地震剖面上會有一定的反映，出現(xiàn)“平點”、“亮點”、“暗點”等現(xiàn)象，這些參量間的相互關(guān)系構(gòu)成了振幅屬性在儲層預(yù)測應(yīng)用中的理論基礎(chǔ)。瞬時振幅能夠指示地層中巖性或油氣聚集引起的變化，在亮點和暗點型油氣藏和儲層的分析預(yù)測中有比較好的應(yīng)用效果。2） 瞬時頻率頻率屬性與振幅屬性有相通的地方，它的變化實際上也是由層速度的改變引起的，而引起層速度變化的原因不外乎也是儲層特性發(fā)生改變造成的。眾所周知，地震波穿過含油氣的儲層時高頻成分會迅速衰減，使接收到的反射波具有明顯的低頻特征，所以頻率屬性在油氣預(yù)測中具有不可低估的作用。瞬時

21、頻率能夠反映地層巖性、流體性質(zhì)等物性的改變，也能用于分析地層接觸關(guān)系等，相對來說瞬時頻率的約束條件較少，但應(yīng)注重低頻成分的獲取和高頻成分的保持。3） 瞬時相位在對地層橫向不連續(xù)性和巖層巖性尖滅的研究中，瞬時相位應(yīng)用效果比較顯著，它是刻畫同相軸連續(xù)性的標(biāo)尺，且不依賴于能量的強弱，地震剖面上的有效波無論振幅大小，在相位剖面上都能很好地顯示出來。當(dāng)?shù)卣鸩ㄔ诰鶆蚪橘|(zhì)中傳播時，它的相位是連續(xù)變化的；但當(dāng)?shù)卣鸩ń?jīng)過有異常存在的介質(zhì)時，它的相位值會在異常對應(yīng)的位置發(fā)生顯著變化，失去連續(xù)性。因此，利用瞬時相位能夠較好的對地下分層和地下異常進(jìn)行辨別，當(dāng)相位圖中出現(xiàn)明顯的幅值變化，相位曲線不再連續(xù)時，就可以判斷該

22、處存在分層或異常。5、 主要的程序?qū)崿F(xiàn)3-1橫剖面 #include <stdio.h>6、 #include <stdlib.h>7、 #include <math.h>8、 #include <msp.h>9、 void mcmpfft(complex x,int n,int isign)10、 11、 complex t,z,ce;12、 float pisign;13、 int mr,m,l,j,i,nn;14、 15、 for(i=1;i<=16;i+)16、 17、 nn=pow(2,i);18、 if(n=nn) break;

23、19、 20、 if(i>16)21、 22、 printf(" N is not a power of 2 ! n");23、 return;24、 25、 z.real=0.0;26、 pisign=4*isign*atan(1.);27、 mr=0;28、 for(m=1;m<n;m+)29、 l=n;30、 while(mr+l>=n)l=l/2;31、 mr=mr%l+l;32、 if(mr<=m)33、 continue;34、 t.real=xm.real;t.imag=xm.imag;35、 xm.real=xmr.real;xm.i

24、mag=xmr.imag;36、 xmr.real=t.real;xmr.imag=t.imag;37、 38、 /*-*/39、 l=1;40、 while(1)41、 42、 if(l>=n)43、 44、 if(isign=-1)45、 return;46、 for(j=0;j<n;j+)47、 48、 xj.real=xj.real/n;49、 xj.imag=xj.imag/n;50、 51、 return;52、 53、 for(m=0;m<l;m+)54、 55、 for(i=m;i<n;i=i+2*l)56、 57、 z.imag=m*pisign/l;

25、58、 ce=cexp(z);59、 t.real=xi+l.real*ce.real-xi+l.imag*ce.imag;60、 t.imag=xi+l.real*ce.imag+xi+l.imag*ce.real;61、 xi+l.real=xi.real-t.real;62、 xi+l.imag=xi.imag-t.imag;63、 xi.real=xi.real+t.real;64、 xi.imag=xi.imag+t.imag;65、 66、 67、 l=2*l;68、 69、 /傅里葉變換70、 void main()71、 72、 73、 74、 short int reelhe

26、ader200; /卷頭信息75、 char SM3200; /說明信息76、 float shuju501; /學(xué)號提出的橫剖面第一道數(shù)據(jù)77、78、79、 float shu256;80、 float traceheader120; /道頭信息81、 82、 struct s83、 84、 float real; 85、 float imag;86、 z1256; /定義結(jié)構(gòu)體，儲存復(fù)數(shù)的實部和虛部87、88、 int i,j,n=256,isign=1,k;89、 FILE *fp,*fp1;90、91、 fp=fopen("stackieee.sgy","

27、rb"); /打開原文件92、 fp1=fopen("橫剖面振幅.xls","wb");93、 fseek(fp,0L,0);94、95、 fread(SM,3200,1,fp);96、 fwrite(SM,3200,1,fp1); /寫入說明信息97、 98、 fread(reelheader,400,1,fp);99、 reelheader10=256; /卷頭第21-22字節(jié)-卷頭信息100、 fwrite(reelheader,400,1,fp1); /寫入卷頭信息101、102、 fseek(fp,16609200L,0); /取第2

28、7個橫刨面103、 for(i=0;i<200;i+)104、 105、 fread(traceheader,240,1,fp);106、 traceheader57=256; /115-116字節(jié)-采樣點數(shù)107、 fwrite(traceheader,240,1,fp1); /寫入道頭信息108、 fread(shuju,2004,1,fp);109、 for(j=0;j<256;j+)110、 111、 z1j.real=shujuj;112、 z1j.imag=0.0;113、 114、 mcmpfft(z1, 256, -1); /傅里葉變換 115、 for(k=0;k

29、<256;k+)116、 shuk=z1k.real;117、 fwrite(shu,1024,1,fp1);118、 119、 3-2縱剖面120、 #include <stdio.h>121、 #include <stdlib.h>122、 #include <math.h>123、 #include "msp.h"124、 void mcmpfft(complex x,int n,int isign)125、 126、 complex t,z,ce;127、 float pisign;128、 int mr,m,l,j,i,n

30、n;129、 130、 for(i=1;i<=16;i+)131、 132、 nn=pow(2,i);133、 if(n=nn) break;134、 135、 if(i>16)136、 137、 printf(" N is not a power of 2 ! n");138、 return;139、 140、 z.real=0.0;141、 pisign=4*isign*atan(1.);142、 mr=0;143、 for(m=1;m<n;m+)144、 l=n;145、 while(mr+l>=n)l=l/2;146、 mr=mr%l+l;1

31、47、 if(mr<=m)148、 continue;149、 t.real=xm.real;t.imag=xm.imag;150、 xm.real=xmr.real;xm.imag=xmr.imag;151、 xmr.real=t.real;xmr.imag=t.imag;152、 153、 /*-*/154、 l=1;155、 while(1)156、 157、 if(l>=n)158、 159、 if(isign=-1)160、 return;161、 for(j=0;j<n;j+)162、 163、 xj.real=xj.real/n;164、 xj.imag=xj.

32、imag/n;165、 166、 return;167、 168、 for(m=0;m<l;m+)169、 170、 for(i=m;i<n;i=i+2*l)171、 172、 z.imag=m*pisign/l;173、 ce=cexp(z);174、 t.real=xi+l.real*ce.real-xi+l.imag*ce.imag;175、 t.imag=xi+l.real*ce.imag+xi+l.imag*ce.real;176、 xi+l.real=xi.real-t.real;177、 xi+l.imag=xi.imag-t.imag;178、 xi.real=xi

33、.real+t.real;179、 xi.imag=xi.imag+t.imag;180、 181、 182、 l=2*l;183、 184、 /傅里葉變換185、 void main()186、 short int reelheader200; /卷頭信息187、 char SM3200; /說明信息188、 float shuju501; /學(xué)號提出的橫剖面第一道數(shù)據(jù)189、 float shuju1256; /進(jìn)行傅里葉變換的數(shù)據(jù)190、 float shu256;191、 float traceheader120; /道頭信息192、 int i,j,k,n=256,isign=1;1

34、93、 struct s194、 195、 float real; 196、 float imag;197、 z1256;198、 FILE *fp,*fp1;199、 fp=fopen("stackieee.sgy","rb"); /打開原文件200、 fp1=fopen("縱剖面.sgy","wb");201、 fseek(fp,0L,0);202、 fread(SM,3200,1,fp);203、 fwrite(SM,3200,1,fp1); /寫入說明信息204、 fread(reelheader,400,

35、1,fp);205、 reelheader10=256; /卷頭第21-22字節(jié)-卷頭信息206、 fwrite(reelheader,400,1,fp1); /寫入卷頭信息207、 fseek(fp,387324L,0); /取第109個縱刨面208、 for(i=0;i<60;i+)209、 210、 fread(traceheader,240,1,fp);211、 traceheader57=256; /115-116字節(jié)-采樣點數(shù)212、 fwrite(traceheader,240,1,fp1); /寫入道頭信息213、 fread(shuju,2004,1,fp);214、

36、for(j=0;j<256;j+)215、 216、 z1j.real=shujuj;217、 z1j.imag=0.0;218、 219、 mcmpfft(z1,256, -1); /傅里葉變換 220、 for(k=0;k<256;k+)221、 shuk=z1k.real;222、 fwrite(shu,1024,1,fp1);223、 fseek(fp,387324L,1); /跳到下一條線該道224、 225、226、227、 JOB 4 計算地震道振幅譜一、要求：提取in-line為110+學(xué)號、cross-line為397-學(xué)號的地震道，通過傅里葉變換求取該地震道的振

37、幅譜，并用Excel作圖。方法：通過連續(xù)傅里葉變換可以求取信號的振幅譜，而實際計算機處理中我們是通過快速傅里葉變換來計算的信號的連續(xù)傅里葉變換，進(jìn)而獲得信號的振幅譜。下面通過國立臺灣大學(xué)丁建均副教授的教程介紹如何通過快速傅里葉變換計算信號的連續(xù)傅里葉變換。參考文獻(xiàn)：Implementing the Continuous Fourier Transform by the Discrete Fourier Transform二、主要的程序?qū)崿F(xiàn)4-1橫剖面振幅#include <stdio.h>#include <stdlib.h>#include <math.h>

38、;#include "msp.h"void mcmpfft(complex x,int n,int isign)complex t,z,ce; float pisign; int mr,m,l,j,i,nn; for(i=1;i<=16;i+) nn=pow(2,i); if(n=nn) break; if(i>16)printf(" N is not a power of 2 ! n");return;z.real=0.0; pisign=4*isign*atan(1.); mr=0; for(m=1;m<n;m+) l=n; whi

39、le(mr+l>=n)l=l/2; mr=mr%l+l; if(mr<=m) continue; t.real=xm.real;t.imag=xm.imag; xm.real=xmr.real;xm.imag=xmr.imag; xmr.real=t.real;xmr.imag=t.imag; /*-*/l=1;while(1) if(l>=n) if(isign=-1) return; for(j=0;j<n;j+) xj.real=xj.real/n; xj.imag=xj.imag/n; return; for(m=0;m<l;m+) for(i=m;i&l

40、t;n;i=i+2*l) z.imag=m*pisign/l; ce=cexp(z); t.real=xi+l.real*ce.real-xi+l.imag*ce.imag; t.imag=xi+l.real*ce.imag+xi+l.imag*ce.real; xi+l.real=xi.real-t.real; xi+l.imag=xi.imag-t.imag; xi.real=xi.real+t.real; xi.imag=xi.imag+t.imag; l=2*l; /傅里葉變換void main() short int reelheader200; /卷頭 char SM3200; /

41、說明float shuju501; /學(xué)號提出的橫剖面第一道數(shù)據(jù)float shu256; float traceheader120; /道頭 float b128; struct sfloat real; float imag;z1256; int i,j,n=256,isign=1,k; FILE *fp,*fp1; fp=fopen("stackieee.sgy","rb"); /打開原文件 fp1=fopen("橫剖面振幅.xls","w+"); fseek(fp,0L,0); fread(SM,3200,

42、1,fp); /fwrite(SM,3200,1,fp1);/寫入說明信息 fread(reelheader,400,1,fp); reelheader10=256; /卷頭第21-22字節(jié)-卷頭信息 /fwrite(reelheader,400,1,fp1);/寫入卷頭信息 fseek(fp,9428400L,0); /取第27個橫刨面 for(i=0;i<200;i+)fread(traceheader,240,1,fp); traceheader57=256; /115-116字節(jié)-采樣點數(shù) / fwrite(traceheader,240,1,fp1); /寫入道頭信息 frea

43、d(shuju,2004,1,fp);for(j=0;j<256;j+)z1j.real=shujuj;z1j.imag=0.0;mcmpfft(z1, 256, -1); /傅里葉變換 for(k=0;k<128;k+) shuk=sqrt(z1k.imag*z1k.imag+z1k.real*z1k.real);/橫剖面振幅 fprintf(fp1,"%ft",shuk); fprintf(fp1,"n"); for(i=0;i<128;i+) bi=500.0/128.0*(i+1); /將EXCEL中橫坐標(biāo)改為bi fprintf(fp1,"%ft",bi);4-2縱剖面振幅#include <stdio.h>#inclu