SEGY课程设计

1、成 绩评阅人成都理工大学工程技术学院地震资料分析课程设计结业报告任务号任务名称提交资料1SEGY格式文件读写1文本报告、程序（电子版）2瞬时属性分析技术文本报告、程序（电子版）3SEGY格式数据操作2程序（电子版）4计算地震道振幅谱程序（电子版）班 级：2012级勘查技术与工程2班姓 名：李南吉学 号：201220911222提交日期：2015年5月11日 JOB 1 SEGY格式文件读写1一、 实验内容1. 了解常见的磁带记录格式2. 编程实现SEG-Y格式读写，并提取inline为（110+学号）的横剖面和cross-line（400-学号）的纵剖面。二、 实验相关软件和数据说明1. Se

2、gy文件分析工具2. Fimage地震资料显示程序3. stack.sgy 三维地震资料(in-line范围105-164,cross-line范围：251-450，time范围500-1500ms)三、 实验内容和步骤1. SEGY格式地震数据读写地震勘探主要包括三个阶段，分别为野外地震数据采集、室内地震资料处理和地震资料解释。且现代资料处理和解释都是在计算机上进行的，因此数据必须通过相应的介质交换和存储，常见的存储介质包括磁带和盘类（硬盘、光盘、U盘等）。由于野外工作环境和早期的工艺限制，地震队野外施工采集到的数据都是存储在磁带上的，近些年由于工艺进步，室内处理解释中数据存储和交换介质基本

3、都是通过硬盘。不管是磁带介质还是盘类介质，其数据都是以一定的二进制格式存储的。所谓格式，指的是地震数据以及各种信息在文件内部的存放方式及顺序。常见的数据格式有SEGY、SEGD、SEG2等，这些格式都是由勘探地球物理学家协会（Society of Exploration Geophysicists，SEG）制定的（全世界统一的地震数据格式）。SEGD和EEG2是野外（时序）记录格式，其中SEGD格式是野外使用最多的格式，SEG2格式仅适用于PC工作环境，目前只有工程地震勘探中在使用。SEGY是室内(道序)记录格式（解编后的炮集记录，CMP道集，叠加和偏移剖面均为此格式），也是室内处理解释作业最

4、为常用的标准格式，目前各种商业处理软件都提供SEGY格式数据加载功能，本讲义也是也是以SEGY格式为标准进行讲解的。图1 SEGY文件格式1) SEGY数据格式SEGY数据格式由三部分组成（如图1所示）。第一部分是3200字节的EBCDIC卡片头（编码方式不同于ASCII）（图2），称其为说明信息（text header），包括40个对磁带进行描述的文本数据（例如，每行包括08个字符的40行文本）。这部分一般包括工区信息、数据采集者、处理者和解释者等信息。第二部分是400字节的二进制卷头（Reel Header 或 binary header）（表1），描述了该磁带卷上内容的有关信息。SEGY

5、格式的第三部分是实际地震道数据。每道含有240字节的道头（Trace Header）（表2），数据记录在道头后面，可能使用4种32位IBM浮点型格式中的一种。（注意，IBM格式与现代普通IBM PC机上所通用的IEEE格式不同，在PC机进行数据处理前需要先进行数据转换）。2）SEGY格式编程思路3200字节说明部分，主要记录采集处理信息，一般处理地震数据不用处理这部分的内容。这部分内容只需要读出来再写入新的文件中去即可，编程时我们只需要定义一个长度为3200的char类型数组即可。卷头中包含了采样间隔和采样点数的重要信息，编程时定义为short int类型数组，就可以得到这些参数了。数据总道数

6、判断，用fseek函数分别指向文件起始和末尾，偏移量就为文件大小（字节）LEN，通过如下公式可以求出总道数总道数=（LEN-3200-400）/（每道采样点数*4+240）三维数据体比二维数据体多了一维，多的这一维代表线号，和二维数据体相比较，如果二维数据体代表的是一条线的数据，那三维数据体就是代表一个工作面的数据。由此我们也知道通常情况下，三维数据体的数据通常情况下会比二维数据体大很多。三维数据体里面的数据格式如下（1） 卷头部分，3200字节+400字节（2） 第一线(inline 1)数据：第一道（cross-line 1）道头部分，第一道数据部分 第二道（cross-line 2）道头

7、部分，第二道数据部分 ， （3） 第二线(inline 2)数据：第一道道头部分，第一道数据部分 第二道道头部分，第二道数据部分 ， （4） ： ， （5） 最后一线数据：第一道道头部分，第一道数据部分 第二道道头部分，第二道数据部分 ， 四、程序实现#include<stdio.h>void main()char SM3200; /定义说明short int reelheader200; /定义卷头 float traceheader60; /定义道头int i,LEN,N,n; / N为总道数 LEN为文件大小（字节）FILE *fp,*fp1,*fp2;fp=fopen(&q

8、uot;112","rb"); /打开原文件后转换成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); /总道数 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数据格式操作2一、实验内容：1. 地震数据再采样，将原始地震数据进行两倍抽稀2. 读取IBM浮点数格式数据2. 尝试实现IEEE和IBM浮点数格式转换二、数据说明stackieee.sgy IEEE格式数据（微机格式）stackieeeibm.sgy IBM格式数据（工作站格式）三、实验内容和步骤1. 地震数据再采样地震资料处理中，经常需要对地震数据进行再采样，一般

11、进行地震数据的抽稀操作，这样可以降低数据大小，加快处理速度。抽稀过程中需要注意抽稀后信号采样率应该满足采样定理。实现方法：1） 改变卷头中采样间隔（17 -18字节，单位：微妙）和每道采样点数（21-22）。2） 将数据道进行抽稀。2. 读取IBM浮点数格式数据参考文献：1）32位IEEE和IBM浮点数结构及其转换方法2）三维地震数据体任意剖面抽取算法及C语言程序实现3. IEEE和IBM浮点数格式转换特定计算机硬件平台，只能处理与自己硬件平台相同格式的浮点数。地震资料处理和解释通常使用IBM、SUN等大型计算机系统。不同体系架构的计算机系统之间的浮点数交换与共享，对于地震数据处理具有重要意义

12、。32位IEEE和IBM浮点数是SEGY地震数据中两种常用的数据格式，因此有必要掌握这两种格式之间的转换方法。参考文献：1）32位IEEE和IBM浮点数结构及其转换方法2）三维地震数据体任意剖面抽取算法及C语言程序实现四、主要的程序实现#include<stdio.h>void main()char SM3200; /定义说明short int reelheader200; /定义卷头short int traceheader120; /定义道头int i,j=0,k;float a501,yy251; / a为原采样点数,yy为更改后采样点数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字节-采样间隔reelheader10=251

14、; /21-22字节-每道采样点数 fwrite(reelheader,400,1,fp1); / 写入卷头 fseek(fp,3600L,0); for(i=0;i<12000;i+) fread(traceheader,240,1,fp); /读出道头 traceheader57=251; / 115-116字节-每道采样点数 traceheader58=4000; / 116-117字节-采样间隔 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); /输出抽稀后数据 fclose(fp); fclose(fp1); JOB 3 瞬时属性分析技术要求：计算地震数据体三瞬属性，并分别提取inline为（120+学号）的横剖面和cross-line（375-学号）的纵剖面的瞬时振幅、瞬时相位、瞬时频率剖面。快速傅里叶变换（FFT）程序：MCMPFFT.C说明：基于Cooley-Tukey radix-2 DIT算法的FFT20世纪70年代后期， Taner等人（1979）在全面研究复数地震道三瞬参数分析原理的基础上，开始利用复数地震道分析地震信号的瞬时振幅（Instantaneous Ampli

16、tude）、瞬时相位（Instantaneous Phase）和瞬时频率（Instantaneous Frequency） ，这就是著名的三瞬参数分析，用于描述地震道中谱信息随时间的变化属性（Taner和Sheriff等，1977）。对实际地震资料的复数道分析做出了开创性贡献。迄今为止，利用复数道分析（解析信号）得到的三瞬时参数又衍生出了许多属性，因此，基于复数道的瞬时属性分析仍然是地震资料解释中最重要的属性之一。1. Hilbert（希尔伯特）变换信号x(t)的希尔伯特变换定义为记为或。显然是x(t)通过LTI系统的输出，该系统的频率响应为如图1(b)所示，它的幅频特性全为1。图1(a)至(

17、c)给出了其频域作用。图1 希尔伯特变换希尔伯特变换的实质是对信号的正、负频率部分分别实施相移和。由实信号x(t)与它的希尔伯特变换构成的复信号为称为x(t)的解析信号(analytic signal)或信号预包络(pre-envelope)，有时也将解析信号称为信号的希尔伯特变换。可以看出而式中是频域的阶跃函数。因此，如果信号x(t)是确定信号，则解析信号也是确定的，且其频谱为实部和虚部可看成是解析信号（复数道）在两个互相正交的平面上的投影（如图2）。图2 解析信号的直观理解由傅里叶变换的特性知，实信号频谱的负频率部分完全是冗余的，它可以从频谱的正频率部分获得。希尔伯特变换是在通信理论中引入

18、的（Gabor，1946），在无线通信中，常将实信号频谱的负频率部分去掉，只保留其正频率部分，这样频谱就不再具有共轭对称性，这种频谱对应的时域信号就是解析信号。已成为通信理论研究的一个重要工具。从上面分析可以看出，解析信号的算法有两种，分别为时域算法和频域算法。本课程设计中我们讲解频域算法。信号x(t)的解析信号z(t)的频域算法如下1. 计算信号x(t)的快速傅里叶变换（FFT）X()。2. 创建数组h与x()大小相同，其值为3. 计算ht和X的点积（Z =ht.* X）4. 计算Z 的傅里叶反变换即为x(t)的希尔伯特变换。2. 三瞬参数计算瞬时属性定义如下：1） 瞬时振幅：2） 瞬时相位

19、：3） 瞬时频率：为避免对瞬时相位求导带来的误差（因为相位与振幅无关，微弱的噪声都会对相位产生严重影响，会带来偏差，如果用标准公式中的相位对时间求导计算瞬时频率，会使误差发生传递和放大）,故利用解析信号本身求瞬时频率：式中:Imag为复数取虚部算子，T为采样间隔。三瞬属性剖面示例如下图所示。图3 三瞬属性示例2. 三瞬属性的地质意义通过各种数学变换计算出的地震属性多达上百种，并不是所有的属性都有明确的地质意义。对于地震属性的分类，目前并没有统一的标准，按照不同的划分方法，同一属性的归属类别也不一样，此处不详究某一属性的归属问题，只是简单介绍一下常用的地震时频属性及其它们的物理和地质意义。1）

20、瞬时振幅地层岩性的变化或储层中流体性质、孔隙度的改变能够引起层速度的变化，而这些变化在地震剖面上会有一定的反映，出现“平点”、“亮点”、“暗点”等现象，这些参量间的相互关系构成了振幅属性在储层预测应用中的理论基础。瞬时振幅能够指示地层中岩性或油气聚集引起的变化，在亮点和暗点型油气藏和储层的分析预测中有比较好的应用效果。2） 瞬时频率频率属性与振幅属性有相通的地方，它的变化实际上也是由层速度的改变引起的，而引起层速度变化的原因不外乎也是储层特性发生改变造成的。众所周知，地震波穿过含油气的储层时高频成分会迅速衰减，使接收到的反射波具有明显的低频特征，所以频率属性在油气预测中具有不可低估的作用。瞬时

21、频率能够反映地层岩性、流体性质等物性的改变，也能用于分析地层接触关系等，相对来说瞬时频率的约束条件较少，但应注重低频成分的获取和高频成分的保持。3） 瞬时相位在对地层横向不连续性和岩层岩性尖灭的研究中，瞬时相位应用效果比较显著，它是刻画同相轴连续性的标尺，且不依赖于能量的强弱，地震剖面上的有效波无论振幅大小，在相位剖面上都能很好地显示出来。当地震波在均匀介质中传播时，它的相位是连续变化的；但当地震波经过有异常存在的介质时，它的相位值会在异常对应的位置发生显著变化，失去连续性。因此，利用瞬时相位能够较好的对地下分层和地下异常进行辨别，当相位图中出现明显的幅值变化，相位曲线不再连续时，就可以判断该

22、处存在分层或异常。5、 主要的程序实现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; /学号提出的横剖面第一道数据77、78、79、 float shu256;80、 float traceheader120; /道头信息81、 82、 struct s83、 84、 float real; 85、 float imag;86、 z1256; /定义结构体，储存复数的实部和虚部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字节-卷头信息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字节-采样点数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; /学号提出的横剖面第一道数据189、 float shuju1256; /进行傅里叶变换的数据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字节-卷头信息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字节-采样点数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+学号、cross-line为397-学号的地震道，通过傅里叶变换求取该地震道的振

37、幅谱，并用Excel作图。方法：通过连续傅里叶变换可以求取信号的振幅谱，而实际计算机处理中我们是通过快速傅里叶变换来计算的信号的连续傅里叶变换，进而获得信号的振幅谱。下面通过国立台湾大学丁建均副教授的教程介绍如何通过快速傅里叶变换计算信号的连续傅里叶变换。参考文献：Implementing the Continuous Fourier Transform by the Discrete Fourier Transform二、主要的程序实现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; /学号提出的横剖面第一道数据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字节-卷头信息 /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字节-采样点数 / 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 fprintf(fp1,"%ft",bi);4-2纵剖面振幅#include <stdio.h>#inclu