常用SEISMIC BINARY数据的读取

1.使用MATLAB读取外部SEISMIC BINARY DATA

由于SEISMIC BINARY DATA在工作中经常碰到。例如SEGY, SEG2格式的地震数据在石油勘探、天然地震以及工程物探中经常被存储在磁盘文件中。科技工作者们需要将这些磁盘中的SEISMIC BINARY DATA导入到计算机内存中。因此，这个工作可以利用FILE OPEN 和 READ, WRITE语句以及CLOSE语句实现。编程过程比较死板，需要一步一步来，遇到有重复步骤的可以稍微调用循环语句实现。为了对编程过程数据格式进行保密或者不触犯相关保密条例。以下SEISMIC BINARY DATA一律采用某数据格式来表达。

某二进制SEISMIC DATA的读取

path(path,'D:');
id=fopen('c4','r','ieee-le');%ieee-le为little endian;ieee-be为big endian
fseek(id,0,'bof');
fseek(id,4,'bof');
ntr=fread(id,1,'int16');
fseek(id,0,'bof');
fseek(id,10,'bof');
nt=fread(id,1,'int16');
data=zeros(nt,ntr);
data1=zeros(nt,ntr);
fseek(id,0,'bof');
fseek(id,14,'bof');
dt=fread(id,1,'int32');
for i=1:ntr
fseek(id,0,'bof');
fseek(id,512+nt*4*(i-1),'bof');
for j=1:nt
data(j,i)=fread(id,1,'int32');
end

for j=1:nt
data(j,i)=bitshift(int32(data(j,i)),0);%32位不做任何偏移，0偏移
end
end
%d=quantizer([32 16]);
%data1=bin2num(d,data);
data1=data.*1e10;
%data1=hex2dec(data);

image(data1(:,1:ntr));
fclose(id);
以上代码可以顺利将外部磁盘中的SEISMIC BINARY DATA读取到内存中。特别是MATLAB的MAT二进制格式是处理所有数据的基础。

在matlab version 5中，MAT文件由一个128字节的文件头和若干个数据单元组成。每个数据单元有一个8个字节的tag，用于说明数据单元的占用的字节数（不包括tag的8个字节）和数据类型。

文件头header里有124字节的文本描述区域和4个字节的flag。flag中的前2个字节说明version，后两个字节是endian indicator。文本描述区域主要说明MAT文件的版本，创建于哪个平台，创建时间。flag中的version说明的是创建这个MAT文件的matlab的版本。edian indicator包括两个字符M和I。

     char mat_data_fhead1[51] = 

           {"MATLAB 5.0 MAT-file, Platform: PCWIN, Created on: "};   

     char mat_data_fhead2[51] =   {"                                                  "};    

     char mat_data_fhead3[4] = {0, 0x01, 0x49, 0x4d};   

     char* datetime = NULL;   

     time_t ltime;   

     tm* today;   

     time(ltime);   

     today = localtime(ltime);   

     datetime = asctime(today);   

     fwrite(mat_data_fhead1, 1, 50, fp);   

     fwrite(datetime, 1, 24, fp);   

     fwrite(mat_data_fhead2, 1, 50, fp);   

     fwrite(mat_data_fhead3, 1, 4, fp);

关于edian：endian: The ordering of bytes in a multi-byte number.
定义：在计算机系统体系结构中用来描述在多字节数中各个字节的存储顺序。相关概念还有MSB（Most Significant Bit/Byte）和LSB（Least Significant Bit/Byte）。在所有的介绍字节序的文章中都会提到字节序分为两类：Big-Endian和Little-Endian。引用标准的Big-Endian和Little-Endian的定义如下：
a) Little-Endian就是低位字节排放在内存的低地址端，高位字节排放在内存的高地址端。
b) Big-Endian就是高位字节排放在内存的低地址端，低位字节排放在内存的高地址端。
c) 网络字节序：TCP/IP各层协议将字节序定义为Big-Endian，因此TCP/IP协议中使用的字节序通常称之为网络字节序。
PS：有些文章中称低位字节为最低有效位，高位字节为最高有效位。

如果edian indicator中的值为MI，则读取MAT数据时应该用IM的顺序。若对于16bit的数据，则要进行两个字节数据的交换。

数据单元的格式

每个数据单元开头都有8个字节的tag用于说明数据单元存储的数据类型和字节数（不包括tag的8个字节）。version5支持多种数据类型。data type中的值为1到14。除了用数值表示某种类型外，还用标识单词联系一种类型。例如data type中存储的是数值1时，代表8bit singed，它的标识单词就为miINT8，方便了用户记忆。值14的标志单词是miMATRIX，代表一种矩阵数据。

数据单元tag中的字节数是每个数据单元不包括8个字节tag的数据字节个数。

接下来的就是存储的数据。数据需要64bit对齐，不够时要补齐到64bit。数据类型是miMATRIX时，数据单元tag中字节数包括矩阵中每个padding的数据个数。其他数据类型时，字节数不包括padding的个数。

当存储的数据不超过4个字节时，还可以采用压缩的数据单元格式。用4个字节存储tag，另外4个字节存储数据。在编程时，tag的前两个字节不为零时，则说明采用的是压缩的数据单元格式。在把数据写入MAT文件中时，压缩的数据单元格式是优先选择的。

datatype值为14的数据类型是：array data，包括了各种类型的array，如数值矩阵，字符矩阵，稀疏矩阵。是一种复合类型结构。字节数包括所有subelement字节数之和。每个subelement都有自己的tag。主要有array flags， dimensions array subelement， array name subelement， real part（pr）subelement， image part subelement。

下图为MAT文件格式的FORTRAN,C等语言的读写程序：

2.python中的obspy及sgyio库的调用

下面以obspy库的调用为例，打开一个sgy文件的命令：

import obspy

t=obspy.read('Model+GathEP-Z.sgy')

tr=t[0]

 print(tr.stats)
         network:
         station:
        location:
         channel:
       starttime: 1970-01-01T00:00:00.000000Z
         endtime: 1970-01-01T00:00:00.485200Z
   sampling_rate: 2500.0
           delta: 0.0004
            npts: 1214
           calib: 1.0
         _format: SEGY
            segy: AttribDict({'trace_header': LazyTraceHeaderAttribDict({'unpacked_header': b'\x00\x00\x00\x01\x00\x00\x00\x01\x00\x00\x00\x01\x00\x00\x00\x01\x00\x00\x00\x01\x00\x00\x00\x01\x00\x00\x00\x01\x00\x01\x00\x00\x00\x00\x00\x01\xff\xff\xfe\x0c\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\xff\x9c\xff\x9c\x00\x00\xc3P\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x04\xbe\x01\x90\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00', 'endian': '>', 'sample_interval_in_ms_for_this_trace': 400, 'year_data_recorded': 0})})
>>> tr.plot()

                                                                          单道地震数据显示

对于地震多道剖面数据显示，需要计算偏移距

>>> i=0
>>> for tr in t:
...   i=i+1
...   tr.stats.distance=i*10
...
>>> t.plot(type='section')

>>pip install ipykernel

>>python -m ipykernel install --user --name python33 --display-name python33

 

地震剖面数据显示

读取地震事件的函数为obspy.readevents('*.ndx')