📚 课程目录第3章使用ObsPy读取SEG-Y文件

3. 使用ObsPy读取SEG-Y文件

在上一节中我们已经了解了SEG-Y文件的文本头、二进制头与道头结构。手工解析这些字节固然可行，但实际项目中更常见、更稳妥的做法是借助成熟的地震学库。ObsPy 是 Python 生态中最常用的地震数据处理库之一，它把 SEG-Y 的字节解析、字节序判断、道头映射和道数据读取都封装成了简洁的接口，让我们可以把精力放在数据本身而不是二进制细节上。

本节的核心对象是 ObsPy 中的 obspy.io.segy 模块，最常用的入口是 read_segy() 函数。它返回一个 Stream 对象，其中每个 Trace 对应 SEG-Y 中的一道。理解“文件 → Stream → Trace → 头信息 + 数据数组”这条对象链，是掌握 ObsPy 读取 SEG-Y 的关键。

基本读取流程

最简读取只需一行代码。假设我们有一个示例文件 demo_line1.sgy，其中包含 120 道、每道 1000 个采样点、采样间隔 2 ms：

from obspy.io.segy.segy import read_segy

stream = read_segy("demo_line1.sgy")
print(len(stream))          # 120，道数
tr = stream[0]
print(tr.data.shape)        # (1000,)
print(tr.stats.sampling_rate)  # 500.0 Hz（2 ms 采样间隔）
print(tr.stats.segy.trace_header)

这里 tr.data 是一个 NumPy 一维数组，保存该道的振幅序列；tr.stats 保存采样率、道号、起始时间等元信息；而 SEG-Y 特有的道头字段则挂在 tr.stats.segy.trace_header 下，以字典形式提供，例如 trace_header.trace_sequence_number_within_line、trace_header.source_coordinate_x 等。文本头可通过 stream[0].stats.segy.textual_header 获取，二进制头则在 stream[0].stats.segy.binary_file_header 中。

ObsPy 的 SEG-Y 读取默认会尝试自动判断字节序。如果数据来自不同厂商的采集系统，建议显式指定 endian 参数，避免“看起来能读、数值却不对”的隐性错误。

读取头信息与道数据的实用写法

实际处理中，我们往往需要批量提取道头字段做质量控制。下面的示例把每道的道号、CDP 号和 X 坐标整理成列表，便于后续检查：

import numpy as np

records = []
for tr in stream:
    th = tr.stats.segy.trace_header
    records.append({
        "trace_no": th.trace_sequence_number_within_line,
        "cdp": th.ensemble_number,
        "x": th.source_coordinate_x,
        "ns": tr.stats.npts,
    })

print(records[:3])
# 示例输出（合理假设）：
# [{'trace_no': 1, 'cdp': 101, 'x': 500000, 'ns': 1000}, ...]

若要一次性取出所有道数据组成二维矩阵（道 × 采样点），可用 np.stack([tr.data for tr in stream])，得到形状为 (120, 1000) 的数组，这正是后续滤波、增益、成像的常用输入形式。

处理常见读取错误

字节序错误是最常见的问题。SEG-Y 标准早期多为大端（big-endian），而部分现代系统输出小端。若读出的振幅出现 1e+30 量级的异常值或全为 0，应怀疑字节序。可显式指定：

stream = read_segy("demo_line1.sgy", endian="<")  # 小端
stream = read_segy("demo_line1.sgy", endian=">")  # 大端

文件截断或道数不符时，ObsPy 可能抛出 unexpected end of file 类错误。此时应先核对文件实际大小与二进制头中记录的道数是否一致，必要时用 os.path.getsize() 与“3600 + 道数 × (240 + 采样点数 × 字节数)”做粗算比对。

非标准或损坏的文本头会导致文本头出现乱码，但通常不影响道数据读取。若连道数据也无法解析，说明文件可能并非标准 SEG-Y（例如某些采集系统的私有格式），此时可尝试 segyio 等替代库，或先用十六进制工具确认 3200 字节处是否为文本头起始。

不要用 try/except 简单吞掉读取异常。捕获后应打印文件路径、文件大小和异常类型，否则批量处理时很难定位是哪一口井或哪一条测线出了问题。

批量读取的稳健封装

面向整条测线或整个目录时，建议把读取逻辑封装成函数，统一处理字节序回退与异常记录。下面是一个可直接参考的示例结构：

import os
from obspy.io.segy.segy import read_segy

def safe_read(path):
    for endian in (">", "<"):
        try:
            st = read_segy(path, endian=endian)
            if len(st) > 0 and abs(st[0].data).max() < 1e20:
                return st, endian
        except Exception as e:
            last_err = e
    return None, str(last_err)

results = {}
for f in os.listdir("segy_dir"):
    if f.lower().endswith((".sgy", ".segy")):
        st, info = safe_read(os.path.join("segy_dir", f))
        results[f] = info

这段代码先用大端尝试，失败或数值异常时回退到小端，并把每个文件的读取结果记录下来。判断“数值是否合理”用 abs(data).max() < 1e20 是一个简单而有效的启发式规则，能过滤掉大部分字节序错误导致的荒谬振幅。