Python实现txt格式文件转segy格式文件

---

1 Segy 格式文件介绍

地震数据处理常用的文件格式为segy格式。

标准segy文件一般包括三部分，第一部分是EBCDIC文件头，长度为3200字节，包括40条记录，每条记录80字节。用来保存一些对地震数据体进行描述的信息；

第二部分是二进制文件头，长度为400字节，用来存储描述segy文件的一些关键信息，包括segy文件的数据格式、采样点数、采样间隔、测量单位等一些信息，这些信息一般存储在二进制文件头的固定位置上；

第三部分是实际的地震道，每条地震道都包含240字节的道头信息和地震道数据。道头数据中一般保存该地震道对应的线号、道号、采样点数、大地坐标等信息，但一些关键的参数位置（如线号、道号在道头中的位置）并不固定。

2 Python 编程实现思路

我们将txt数据转为segy数据时，需要自己添加文件头信息，本文添加文件头信息的思路是：先准备一份完整的包含文件头信息的segy数据，然后在此基础上修改相关信息。

第一部分3200字节的EBCDIC文件头信息不需要修改，可以直接写入。

第二部分400字节的二进制文件头，需要修改的信息包括：13-14字节，每个记录的数据道数（每炮道数或总道数）；17-18字节，采样间隔（us）；21-22字节，样点数/每道（道长）。

第三部分240字节的道头信息，需要修改的信息包括：1-4字节，一条测线中的道顺序号。如果一条测线有若干卷带，顺序号连续递增；9-12字节，原始的野外记录号（炮号）；13-16字节，在原始野外记录中的道号；115-116字节，本道的采样点数；117-118字节，本道的采样间隔，以us表示。

3 Python 代码

import numpy as np
import struct
import os

if __name__ == '__main__':
    # 读取地震数据
    nSample = 500    # 每道采样点数
    nTrace = 600    # 地震道数
    S = np.zeros((nSample, nTrace))    # 地震数据
    f = open('Marmousi2/seismic_marmousi-ii.txt')
    data = f.readlines()
    row_z = 0
    for i in range(nSample):
        S[row_z, :] = data[i].split(' ')[0:nTrace]
        row_z += 1
    f.close()

    # 写入地震数据，并添加卷头、道头等信息
    # nTrace = 2301
    # nSample = 751
    fSegy = open("Marmousi2/example_for_segy.segy", "rb")    # 打开一份包含完整文件头信息的segy格式数据
    # fSegy_size = os.path.getsize("Marmousi2/example_for_segy.segy")    # 计算文件所占字节数
    # print((fSegy_size - 3600) / (500 * 4 + 240))    # 计算地震道数
    fileHandle = open('Marmousi2/seismic_marmousi-ii.segy', 'wb')
    cardHeader = fSegy.read(3200)    # 读取3200字EBCDIC文件头
    fileHandle.write(cardHeader)     # 写入EBCDIC文件头
    reelHeader = fSegy.read(400)     # 读取400字节二进制文件头

    # 修改相关参数
    reelHeader0 = bytearray(reelHeader)
    reelHeader0[12:14] = struct.pack('>i', nTrace)[2:4]    # 修改地震道数
    reelHeader0[16:18] = struct.pack('>i', 2000)[2:4]      # 修改采样时间间隔 单位为微秒(us)
    reelHeader0[20:22] = struct.pack('>i', nSample)[2:4]   # 修改地震道道采样点数

    # 测试修改是否正确
    tempValue = reelHeader0[20:22]
    hexValue = tempValue.hex()    # 转为十六进制
    decValue = int(hexValue, 16)
    print("每道采样点数：")
    print(decValue)

    fileHandle.write(bytes(reelHeader0))    # 写入二进制文件头

    # 写入数据体
    for i in range(nTrace):
        # fSegy.seek(240 * i + 3600 + nSample * 4 * i, 0)    # 每一道都跳过240字节道头
        fSegy.seek(3600, 0)    # 跳过3600字节文件头
        traceHeader = fSegy.read(240)
        traceHeader0 = bytearray(traceHeader)
        traceHeader0[0:4] = struct.pack('>i', i + 1)
        traceHeader0[8:12] = struct.pack('>i', i + 1)
        traceHeader0[12:16] = struct.pack('>i', i + 1)
        traceHeader0[114:116] = struct.pack('>i', nSample)[2:4]  # 数据道采样点数
        traceHeader0[116:118] = struct.pack('>i', 2000)[2:4]     # 微秒（us）形式的采样间隔
        fileHandle.write(bytes(traceHeader0))
        for j in range(nSample):
            fileHandle.write(struct.pack('>f', S[j, i]))
    fSegy.close()
    fileHandle.close()

    f = open("Marmousi2/seismic_marmousi-ii.segy", "rb")
    f_size = os.path.getsize("Marmousi2/seismic_marmousi-ii.segy")    # 计算文件所占字节数
    print("地震道数：")
    print((f_size - 3600) / (nSample * 4 + 240))    # 计算地震道数
    f.close()

运行结果：