用C语言读写SGY格式的地震数据文件

Author:Yangwqcumt用用CC语语言言读读写写SSGGYY格格式式的的地地震震数数据据Yangwqcumt百度专用，转载需授权地震勘探野外采集的数据，以及经过资料处理获得的三维数据体，只要是放在计算机里，都是以二进制文件的形式存放的。这些文件的处理显示等工作，一般都可以用商业化的软件进行。但是作为一个从事地震勘探研究的技术人员，有时会有些想法，有某种灵感，但是原有的软件又不允许你去做某种试验以验证你的想法。这时候，自编个小程序显然有必要，而且弄好后你的成就感会很强烈。1.SEGY格式地震数据文件地震数据，是以各种格式存放的。所谓格式，指的是地震数据以及各种信息在文件内部的存放方式及顺序。常见的地震数据格式，有segy格式、seg2格式、segd格式等。同样的格式，还有微机版、工作站版及其它版本。本文仅是入门级材料，我们仅就微机版segy格式进行分析。Segy格式的地震数据文件，属于典型的流式文件，它的信息和数据都是按字节顺序一个个地存放的，每个字节都有其特定的含义。这种格式的文件，由文件头部的3600字节以及地震道组成。文件头前部的3200字节共分为40行，每行80个字符，但这些字符不是ascii码，是一种称为ebcdic的编码。一般这部分都不去读，或者只能显示出来查看其中的内容。接下来是400字节的二进制部分。这里面有长整型数和短整型数，其具体含义参见附录一。Author:Yangwqcumt每个地震道由道头240字节（参见附录一）和地震数据组成。地震数据的个数和类型（指它是浮点数整数还是什么）文件头中有定义。此处我们假定所有的数据都是微机的四字节浮点数。2.文件头部3200字节特殊编码部分的读取。该部分共分为40行，每行80个字符，但这些字符不是ascii码，是一种称为ebcdic的编码。转换成ascii码可以采用查表的方法进行。一般处理地震数据不用读这部分的内容。#includestdio.h#includestdlib.hvoidmain(){unsignedcharE2A[256]={0,1,2,3,156,9,134,127,151,141,142,11,12,13,14,15,16,17,18,19,157,133,8,135,24,25,146,143,28,29,30,31,128,129,130,131,132,10,23,27,136,137,138,139,140,5,6,7,144,145,22,147,148,149,150,4,152,153,154,155,20,21,158,26,32,160,161,162,163,164,165,166,167,168,91,46,60,40,43,33,38,169,170,171,172,173,174,175,176,177,93,36,42,41,59,94,45,47,178,179,180,181,182,183,184,185,124,44,37,95,62,63,186,187,188,189,190,191,192,193,194,96,58,35,64,39,61,34,195,97,98,99,100,101,102,103,104,105,196,197,198,199,200,201,202,106,107,108,109,110,111,112,113,114,203,204,205,206,207,208,209,126,115,116,117,118,119,120,121,122,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,123,65,66,67,68,69,70,71,72,73,232,233,234,235,236,237,125,74,75,76,77,78,79,80,81,82,238,239,240,241,242,243,92,159,83,84,85,86,87,88,89,90,244,245,246,247,248,249,48,49,50,51,52,53,54,55,56,57,250,251,252,253,254,255};inti,j;unsignedcharf3200[3200];FILE*f1;chars1[200];printf(InputSgyFileName:\n);scanf(%s,s1);//输入文件名f1=fopen(s1,rb);fread(f3200,3200,1,f1);fclose(f1);for(i=0;i40;i++){for(j=0;j80;j++)printf(%c,E2A[f3200[i*80+j]]);printf(\n);}}Author:Yangwqcumt3.文件头400字节二进制部分的读取可以把400字节作为若干个长整型及短整型数据读入：#includestdio.h#includestdlib.hvoidmain(){FILE*f1;inti,l;inttraces;unsignedcharf3200[3200];charFileName[200];longints1[100];shortints2[200];printf(输入地震文件名[*.sgy]:);scanf(%s,FileName);f1=fopen(FileName,rb);if(f1==NULL)//文件打开不成就显示错误信息，然后退出。{printf(FileOpenerror!\n);exit(0);}fread(f3200,3200,1,f1);//读入前面的字节fread(s1,400,1,f1);//读入接着的字节，作为长整型数fseek(f1,3200,0);//退回去，fread(s2,400,1,f1);//再以短整型数读入那些个字节//////////////////////////////////for(i=0;i3;i++)printf([%d]:%d\n,i+1,s1[i]);//显示文件头格式说明中的二进制部分的前三项for(i=3;i15;i++)//显示4-1,4-2，。。。，15-1,15-2{printf([%d-1]:%d\n,i+1,s2[2*i]);printf([%d-2]:%d\n,i+1,s2[2*i+1]);}///////////////////////////////////*在文件头中，最重要的是下面几项：5-1:这个文件的地震道的时间采样间隔，单位是微秒；6-1：每个地震道的样点数；7-1:数据的格式。现在一般都是四字节的浮点数，格式取.*/printf(时间采样间隔[微秒]：%d\n,s2[8]);printf(地震道的样点数：%d\n,s2[10]);printf(数据格式代码：%d\n,s2[12]);/*有了这些信息，这个文件含有的地震道数是多少呢？Author:Yangwqcumt若文件长度:l则道数：traces=(l-3600)/(240+s2[10]*4)*/fseek(f1,0,2);//l=ftell(f1);//l此时就是文件的字节数；printf(文件的长度=%d\n,l);traces=(l-3600)/(240+4*s2[10]);printf(地震道数是：%d\n,traces);fclose(f1);}也可以把那个400字节作为一个结构体，该结构体定义见附录二。//这个例子使用了结构体来读取文件头信息#includestdio.h#includestdlib.h#includesegyHeader.hvoidmain(){FILE*f1;inti,l;inttraces,trace_length;charFileName[]=100.sgy;structSegyReelHdrStructFileHeader;f1=fopen(FileName,rb);if(f1==NULL)//文件打开不成就显示错误信息，然后退出。{printf(FileOpenerror!\n);exit(0);}fread(&FileHeader,sizeof(FileHeader),1,f1);trace_length=FileHeader.hns;printf(%s:%d\n,时间采样间隔[微秒],FileHeader.hdt);printf(%s:%d\n,每道的样点数,FileHeader.hns);printf(%s:%d\n,数据的格式代号,FileHeader.format);fseek(f1,0,2);//l=ftell(f1);//l此时就是文件的字节数；printf(文件的长度[字节数]=%d\n,l);traces=(l-3600)/(240+4*trace_length);fclose(f1);printf(%s:%d\n,地震道数,traces);}4.某道数据的读取Author:Yangwqcumtincludestdio.h//包含头文件#includestring.h//#includestdlib.h//voidmain(){FILE*f1,*f2;//定义两个文件指针变量inti,l;charFileName[200];//字符数组，存放文件名intTraces,Trace_length,Trace2read;//顾名思义，这几个变量用来存放:地震道数、地震道的长度、要读的那一个地震道的序号float*TraceData;//定义一个浮点型的指针，一会儿开辟内存后放一个地震道的数据/////////////////////////获得地震那个文件名///////////printf(输入地震文件名[*.sgy]:);scanf(%s,FileName);//当然上面两行，可用一个这样的语句代替：strcpy(FileName,100.sgy);f1=fopen(FileName,rb);//打开文件，打开形式为：二进制读rbif(f1==NULL)//判断打开了没有，不成功就返回吧。{printf(Cannotopeninputfile!\n);//显示信息exit(0);//退出}Trace2read=430;//要读取哪一道，可以键盘输入，为方便在这就给定了，设为第430道Trace_length=800;//一个地震道里面数据的个数,可以从文件头获得，这里先拿来用了Traces=631;//该文件地震道的个数，也可以通过文件头里面的信息设法获得，也是先拿来用l=3600L+(240+Trace_length*4L)*(Trace2read)+240;//要读取的那一个地震道的，数据部分，在这个文件中的开始的位置，这个很重要，下面做点说明/*一般来说，地震勘探的数据文件是有格式的。所谓格式，就是地震数据、及相关的信息在地震文件中的存放顺序。现在研究的这个叫做sgy格式。这种格式比较普遍。它由文件头和地震道组成，文件头字节，每个地震道由字节的道头和数据组成。所以一个sgy格式的文件的组成为：3600字节+240字节+第一道数据+240字节+第二道数据+240字节+第三道数据+。。。。。。。文件头部的字节以及每道的字节道头有信息，这个先不研究。对于我们给的这个数据，道数据的长度800，是指个浮点数，所以其字节数要乘以,得到，加上字节的道头，共计字节，所以第道的开始位置，就是：+430*3440，然后由于我们不读这个地震道的道头，数据部分的位置还要加，l=3600+3440*430+240;请时刻注意，c的数组以及其他什么东西的编号多是从0开始的。*/fseek(f1,l,0);//定位到那个地震道的数据的开始位置TraceData=newfloat[Trace_length];Author:Yangwqcumt//先开一块内存，要不数据没地方放，你看事先定好可以吗，比如：floatTraceData[800];frea