📚 课程目录第4章使用segpy读取与写入SEG-Y

4. 使用segpy读取与写入SEG-Y

segpy 是一个纯 Python 实现的 SEG-Y 读写库，它的核心价值在于把 SEG-Y 文件拆成三个可编程对象：文件级文本头（3200 字节 EBCDIC）、二进制文件头（400 字节）和逐道道头（240 字节）加道数据。相比直接用 struct 手工解析字节流，segpy 提供了字段名映射、字节序自动处理和按道迭代的接口，适合做批量读取、头字段筛选和简单文件的生成。本节围绕“读—改—写”这条主线，把 segpy 的接口用法落到可执行的代码上。

一、segpy 的对象模型与读取入口

segpy 的读取入口是 segy_reader，它返回一个 SeismicData 对象，该对象同时持有文本头、二进制头和道集合。道集合支持按索引访问，也支持迭代，每道是一个 Trace，其 header 属性是道头字段的字典式容器，data 属性是 numpy 数组。理解这个层级关系，是后续自定义头字段和写文件的基础。

from segpy.reader import create_reader
from segpy.writer import write_segy
from segpy.encoding import EBCDIC

with open('line12.sgy', 'rb') as f:
    reader = create_reader(f, encoding=EBCDIC)
    print(reader.num_traces())          # 道数
    print(reader.sample_interval())     # 采样间隔，单位微秒
    print(reader.num_samples())         # 每道样点数
    print(reader.data_format())         # 数据格式码，如 1=IBM浮点, 5=IEEE浮点

    for i, trace in enumerate(reader.traces()):
        h = trace.header
        print(h.trace_sequence_number_within_line,
              h.source_coordinate_x,
              h.cdp_x,
              trace.data[:5])
        if i >= 2:
            break

这里有两个容易踩坑的点。第一，create_reader 需要以二进制模式打开文件，且文件对象必须在读取期间保持打开；如果提前关闭，迭代道时会抛异常。第二，编码参数决定文本头如何解码，国内老资料多为 EBCDIC，新导出文件常见 ASCII，读错编码不会报错但文本头会变成乱码，应先用 reader.text_header() 打印前几行确认。

二、道头字段的命名与自定义

segpy 内置了 SEG-Y rev1 标准道头字段名，例如 trace_sequence_number_within_line、source_coordinate_x、cdp_x、offset、elevation_scalar 等，字段名与字节位置一一对应。读取时直接按名取值，避免了记字节偏移。但实际工区常把自定义信息塞进标准字段未使用的字节区间（如 181–184、233–240 字节），这时需要自定义字段。

segpy 通过 HeaderField 和 Header 子类机制支持自定义。做法是定义一个继承 Header 的类，用 HeaderField 声明字段名、字节偏移、字节数和数据类型，然后在读写时把该类传给 header_class 参数。下面示例在 181–184 字节处定义一个 4 字节整型的“采集班组号”字段：

from segpy.header import Header, HeaderField, FieldType

class MyTraceHeader(Header):
    crew_id = HeaderField(offset=180, size=4,
                          field_type=FieldType.INT32,
                          description='采集班组号')

# 读取时指定自定义头类
with open('line12.sgy', 'rb') as f:
    reader = create_reader(f, encoding=EBCDIC,
                           trace_header_class=MyTraceHeader)
    t0 = next(iter(reader.traces()))
    print(t0.header.crew_id)

偏移量以 0 为起点、按道头起始字节计算，size 必须与字段实际占用字节一致，否则相邻字段会错位。自定义字段只影响 Python 侧的解析视图，不会改变文件本身的字节布局；写入时只要用同一个头类，字段就会落到对应字节位置。

三、写入一个简单 SEG-Y 文件

写入使用 write_segy，需要提供输出文件对象、文本头、二进制头和道序列。文本头是 40 行、每行 80 字符的字符串列表；二进制头用 BinaryReelHeader 构造，至少要设置采样间隔、样点数和数据格式；道序列由 Trace 对象组成，每道包含头和数据。下面生成一个 3 道、每道 100 个样点、IEEE 浮点格式的最小文件：

import numpy as np
from segpy.binary_reel_header import BinaryReelHeader
from segpy.trace_header import TraceHeaderRev1
from segpy.trace import Trace
from segpy.writer import write_segy
from segpy.encoding import EBCDIC

text_header = ['C01 示例工区 line12 测试文件'.ljust(80)] + \
              ['C%02d' % i + ' ' * 77 for i in range(2, 41)]

bin_header = BinaryReelHeader(
    sample_interval=2000,      # 2 ms
    num_samples=100,
    data_format=5,             # 5 = IEEE 浮点
)

traces = []
for i in range(3):
    th = TraceHeaderRev1(
        trace_sequence_number_within_line=i + 1,
        cdp_x=1000 + i * 25,
        source_coordinate_x=1000 + i * 25,
        crew_id=7,             # 自定义字段
    )
    data = np.sin(np.linspace(0, 4 * np.pi, 100)).astype(np.float32)
    traces.append(Trace(header=th, samples=data))

with open('out.sgy', 'wb') as f:
    write_segy(f, bin_header, text_header, traces,
               encoding=EBCDIC, trace_header_class=MyTraceHeader)

写入完成后必须回读校验，这是判断“写对了”的唯一可靠方式。校验至少覆盖三项：道数是否等于 3、sample_interval 是否为 2000、第 2 道的 cdp_x 是否为 1025、crew_id 是否为 7，以及样点数组与写入前是否逐点一致。若回读时 crew_id 为 0 或异常大值，通常是偏移量或字节数写错，或写入时未传同一个 trace_header_class。

自定义字段的偏移和长度建议先在纸上对照 SEG-Y rev1 道头表核对一遍，再写进代码。标准字段占用的字节区间不要复用，否则回读时标准字段会被覆盖成错误值，且这种错误在数据体上完全看不出来。

四、批量处理中的实用约定

把上述接口用于批量任务时，建议固定几条约定：读取阶段只做筛选和统计，不修改原文件；需要改头字段时，读出全部道到内存或流式写出到新文件，避免原地覆盖；每个文件处理完立即回读关键字段并记录到日志，形成“文件—道数—采样率—自定义字段值”的核对表。对于道数上万的大文件，用 reader.traces() 迭代而不是一次性 list()，可以显著降低内存占用。掌握这套读写接口后，后续的批量筛选、头字段修正和格式转换都可以在此基础上组合实现。