Python实现SEG-Y地震数据读取与频谱分析全流程实战

简介本资源是一套面向地震勘探领域初学者与科研人员的MATLAB SEGY数据处理工具集聚焦于地震数据读取、可视化与频谱分析等核心环节。压缩包共含6个.m文件总大小仅4KB轻量实用涵盖SEGY标准格式解析含IBM浮点数转换、地震道wiggle图绘制、GUI交互响应、数据预处理及频率谱计算等功能模块完整覆盖从原始数据加载到特征提取的基础流程。已有555人学习下载适用于高校地球物理专业课程实践、野外数据快速质检或科研项目前期分析。用户可直接调用segy_Read.m解析二进制SEGY头信息与道数据结合Wiggle.m生成专业波形图并通过getFrequencySpectrum.m开展频域特征识别所有脚本结构清晰、注释充分具备良好可读性与二次开发基础。1. 项目整体设计与思路拆解1.1 SEG-Y 格式到底是个什么“盒子”在地震勘探里野外采集回来的数据绝大多数都以 SEG-Y 格式存放它由 SEGSociety of Exploration Geophysicists勘探地球物理学家学会标准化用来记录地震道数据和相关的道头信息。做地震数据处理的人几乎每天都要跟这种二进制文件打交道。很多人第一次拿到 SEG-Y 文件试图用文本编辑器打开看到的是一片乱码和零星的 ASCII 字符然后就开始头大。这个文件结构其实并不复杂但你得先把它想成一个“有固定规则的盒子”最前面是 3600 字节的文本卷头EBCDIC/ASCII 文本头之后是 3200 字节的二进制卷头再往后就是一条条记录道每道由 240 字节的道头和若干字节的数据体组成。SEG-Y 标准实际上有多个版本最常见的是 SEG-Y Rev0、Rev1以及 2002 年之后逐步推广的 Rev2。Rev0 是最早期版本数据结构相对简单Rev1 增加了扩展文本头和更多道头属性的规范Rev2 则扩展了头段、坐标系统以及数据格式码范围。很多实际生产的文件虽然标称符合 SEG-Y 标准但厂商会在道头里放入自定义的扩展信息甚至使用非标准的坐标存储位置。这就是为什么新手按标准偏移量读数据有时会遇到“读数对不上”的问题——不是你的代码错了而是文件本身使用了厂商变体必须先摸清头段信息。标题里的“segy_Read”和“segyread”本质上就是围绕这套格式去做读取、显示和后续分析。很多编程环境里都有名为 segyread 的现成入口比如 Seismic Unix 里的 segyread 命令、MATLAB 里的 segyread 函数Python 生态里也有 segyio、ObsPy 这类封装好的库。自己动手写读取代码不是为了重复造轮子而是为了真正理解这个盒子的内部结构这样遇到格式变体时才不会慌。1.2 为什么选择“先读格式、再做显示、最后频谱分析”这条路线一个完整的地震数据处理流程通常从格式解析开始到波形显示、频谱分析、滤波、增益恢复再到偏移成像等。这个项目只聚焦“读取显示 频谱分析”这一段是因为它处于整个流程的最前端也是数据质量评估的基础。你不会想在一堆还没弄明白采样率、道数、极性是否正确的数据上直接做高端的处理那样只会放大错误。我个人的做法是第一步先把文本卷头打印出来看采集参数、观测系统、处理流程描述第二步读二进制卷头确认采样点数、采样间隔、数据格式码第三步扫描所有道的道头统计道号、坐标、采样点数的分布快速发现坏道和缺道第四步才是把数据体抽出来画波形图最后再对感兴趣的道段做频谱分析。这个顺序可以最大化减少数据误读也为后续处理提供干净的输入。为什么不用现成的软件一把梭因为很多商业软件把读取过程包装得太“黑盒”了。你导入一个 SEG-Y 文件看到漂亮的剖面但如果数据有异常比如某个版本的 IBM 浮点转换错误、某几道采样点数不一致这些软件往往只会用默认参数强行解释容易掩盖问题。写一遍读取脚本相当于给数据做了一个“入厂体检”每一步都落在明处后面出了问题也知道往哪里排查。2. 核心细节解析与实操要点2.1 卷头和道头里必须“抠”出来的关键参数文本卷头占 3200 字节通常是 EBCDIC 编码有的文件也写成 ASCII里面是人工可读的施工说明、处理流程等信息。从第 3201 字节开始是 3200 字节的二进制卷头里面的信息全部是二进制数端序为大端Big-Endian。你需要重点读取的偏移量包括字节 3217–3218采样间隔单位是微秒μs。比如值为 2000表示采样间隔 2000 μs也就是 2 ms对应的采样率是 500 Hz。字节 3221–3222每道采样点数。这个值如果为 0 或非法值说明这个文件可能没有遵守规范可能需要在道头里重新获取每道采样点数。字节 3225–3226数据格式码这个是最关键的参数。1 表示 4 字节 IBM 浮点数2 表示 4 字节定长整型3 表示 2 字节定长整型5 表示 4 字节 IEEE 浮点数8 表示 1 字节整型9 表示 8 字节 IEEE 浮点数等。读错格式码数据基本就是一堆垃圾值。240 字节的道头里也有几个关键位置字节 1–4 是道序号字节 29–30 是道内采样点数有些文件用这个值覆盖卷头的采样点数字节 109–110 是道内采样间隔单位同样是微秒字节 1157–1158 在某些自定义实现中也可能是采样点数。需要提醒的是SEG-Y 道头偏移量的标准在不同版本间有细微差异生产文件又经常塞入厂商扩展属性所以最稳妥的验证方法是先按标准偏移读一批数据出来对前几道的波形做可视化检查振幅是否连续、道间是否有明显跳变再回头调整读取规则。另外千万不要忽略坐标信息。道头里通常会存放炮点坐标和检波点坐标常见的偏移位置位于字节 71–78炮点 X/Y和字节 81–88检波点 X/Y附近但不同厂商定义不同有的放在 181–188、191–198 等扩展位置。坐标影响后续的观测系统定义和静校正如果一开始就记错位置后面的空间位置关系就全乱了。读取脚本里最好把可能用到的坐标偏移都打印出来或者写一个交互式检查函数快速查看文件里说“X 坐标”的那个字段到底是什么量级、什么单位。2.2 从“读出来”到“显示出来”归一化和显示方式的选择地震数据的振幅动态范围极大浅层强反射和深部弱反射可能相差 100 倍以上。如果直接按原始振幅绘制剖面上往往只剩下一两条亮轴深层信息全部被“压扁”。所以在显示之前通常要做振幅归一化或增益处理。一种简单有效的方法是做全局百分位归一化统计所有样点振幅的绝对值取 98% 分位数的绝对值作为归一化因子然后把数据除以这个因子将振幅压缩到 [0, 1] 区间附近。这个策略比简单的最大值归一化更鲁棒因为最大值可能是某个尖峰噪声用它归一化会把其他信号压得太暗。想要更进一步的显示增强可以试试自动增益控制AGC在一个时间窗口内计算 RMS均方根振幅用当前采样点的局部 RMS 去归一化原始振幅。窗口长度一般选取 200–500 ms窗口太短会破坏波形相对振幅窗口太长则增强效果不明显。显示方式上最常用的是变密度图灰度/彩色填充和波形图wiggle trace。波形图适合展示单道细节能直观看到极性、相位和同相轴连续性变密度图适合展示整个剖面的宏观结构通过颜色映射关注振幅强弱。实际项目里我喜欢把两种方式叠加在变密度图的底色上用半透明波形叠加既能看全局又能看细节。颜色表建议避开彩虹色使用“蓝白红”或者“黑红白”这类地震行业常见的规范化色标因为彩虹色会造成不同振幅之间的“等距错觉”不便于解释同相轴的连续性。这里也涉及一个“显示”与“分析”的边界问题显示归一化只是为了让人眼看得清楚用于频谱分析和反演的仍应该是原始振幅数据。如果你把归一化后的数据拿去做频谱分析得到的频谱将无法反映真实的能量衰减趋势。项目里最好把“显示数据”和“计算数据”两条路径分开读取时保留原始数据显示时复制一份再归一化这样两头都不耽误。2.3 频谱分析之前要理解地震信号的几个“脾气”要做频谱分析先得知道地震信号长什么样。地震子波的主频通常在 10–80 Hz 之间陆上可控震源或炸药震源激发的地震波浅层传播时高频成分较多随着传播距离增加大地对高频的吸收衰减明显深层信号的主频会明显下降。所以一条地震道的频谱不会是平直的而是呈现“带通”特征低频端有少量能量中间有主频峰高频端逐渐衰减。采样率决定了奈奎斯特频率超过奈奎斯特频率的部分会发生假频混叠。比如采样间隔是 2 ms采样率就是 500 Hz奈奎斯特频率是 250 Hz。地震有效信号能量通常集中在 100 Hz 以内250 Hz 的奈奎斯特频率足够覆盖但如果采样间隔是 4 ms采样率降到 250 Hz奈奎斯特频率只有 125 Hz此时如果数据里有高频干扰就很容易被“折叠”到低频段干扰你对有效信号的判断。做频谱分析时还需要注意“窗”的作用。对整道长信号直接做 FFT频谱会因为端点截断产生频谱泄漏导致本来尖而窄的主频峰变得“胖乎乎”。更科学的做法是先对信号做去均值和去趋势处理再乘一个汉宁窗或布莱克曼窗这样信号两端平滑过渡到零频谱上的旁瓣会大幅减少。窗口长度也很关键。对 2 ms 采样、共 1000 个采样点的信号直接做 1024 点 FFT频率分辨率约为 0.49 Hz但如果信号长度只有几百毫秒频率分辨率可能就有好几赫兹这会糊掉精细谱结构。建议在做频谱分析前先打印当前道的时间和频率分辨率做到心里有数。3. 实操过程与核心环节实现3.1 环境搭建与数据准备这个项目我用的是 Python 生态原因是数据分析上手快可视化生态完整。核心依赖如下numpy数组运算和 FFT 基础segyio高性能 SEG-Y 读取库支持大文件内存映射matplotlib绘制波形图、频谱图scipy用于信号处理welch 谱估计、窗函数、去趋势。如果你更倾向 ObsPy 的地震学工具链它内部也封装了 SEG-Y 读写接口但 ObsPy 主要面向天然地震数据处理对勘探 SEG-Y 的道头标准覆盖并不全。我的建议是勘探数据优先用 segyio天然地震数据再考虑 ObsPy。安装 segyio 很简单pip install segyio没有真实 SEG-Y 数据时可以自己构造一个程序生成的“模拟 SEG-Y”。使用 segyio 创建文件需要先定义采样点数、道数、采样间隔和数据格式码然后用segyio.create生成文件并向里面写入数据。这样不仅能验证读取流程还能内置已知频谱特征比如 30 Hz 主频的雷克子波用来检验频谱分析代码是否准确。3.2 用 Python 实现 SEG-Y 读取、头段解析与显示下面这段代码是 segy_Read 项目里的核心读取函数它接受一个 SEG-Y 文件名返回卷头信息、道头数组和数据矩阵import segyio import numpy as np def read_segy(file_path): with segyio.open(file_path, r, strictFalse) as f: # 读取二进制卷头核心参数 samples int(f.bin[segyio.BinField.Samples]) # 每道采样点数 interval int(f.bin[segyio.BinField.Interval]) # 采样间隔微秒 # 读取文本卷头 text_header f.text[0] # 读取所有道数据形状为 (道数, 每道采样点数) traces f.trace.raw[:].copy() # 读取道头中常用的属性 trace_num f.attributes(segyio.TraceField.TRACE_SEQUENCE_LINE)[:] sx f.attributes(segyio.TraceField.SourceGroupScalar)[:] # 采样时间轴 time_axis np.arange(samples) * interval / 1000.0 # 转毫秒 return { samples: samples, interval_us: interval, text_header: text_header, traces: traces, trace_num: trace_num, time_axis: time_axis, }关键点有两个一是strictFalse。有些 SEG-Y 文件存在非标准道头严格模式会直接抛出异常非严格模式则能容忍并继续读取这对于处理“野路子”生产数据很有用二是f.trace.raw[:]它按原始存储值返回数据不会帮你做格式码转换之外的处理。如果你打算做归一化显示建议用f.trace[:]的副本如果你要保留原始振幅做频谱分析用raw更合适。显示部分我写了两个函数一个画单道波形另一个画多道变密度剖面。单道波形直接以时间为横轴、振幅为纵轴把曲线画出来即可。变密度图可以用matplotlib的imshow纵轴是时间横轴是道号颜色用地震行业常用的蓝白红色标绘制时记得把每一道做百分位归一化避免强能量道压暗整个剖面的显示。3.3 频谱分析的实现与结果解读读进来一炮数据之后就可以对目标道做频谱分析了。这里的“目标道”可以是炮记录中的某一道也可以是叠加后的某一道。下面用一个完整的函数实现对输入信号做去均值、去趋势、加窗、补零然后计算单边振幅谱并返回频率轴和振幅。import numpy as np from scipy import signal def compute_amplitude_spectrum(data, dt_ms2.0): data: 一维地震道数据 dt_ms: 采样间隔毫秒 fs 1000.0 / dt_ms # 采样率Hz n len(data) # 去均值和线性趋势 detrended signal.detrend(data, typelinear) # 汉宁窗 window np.hanning(n) win_data detrended * window # 补零到下一个2的幂次提升显示分辨率 nfft int(2 ** np.ceil(np.log2(n))) spectrum np.fft.rfft(win_data, nnfft) amp np.abs(spectrum) / np.sum(window) * 2.0 freq np.fft.rfftfreq(nfft, ddt_ms / 1000.0) return freq, amp # 使用示例 freq, amp compute_amplitude_spectrum(traces[100], dt_ms2.0)这段代码里有几个经验值得展开说。去线性趋势非常重要因为地震数据如果存在低频漂移直接加窗做 FFT 会在低频段形成一大团虚假能量淹没真实低频信号。窗函数方面汉宁窗的旁瓣衰减很好适合对地震子波这种带宽有限的信号如果你更关心幅度估计的保真度也可以选择矩形窗但旁瓣会大一些。补零操作并不会提高真实频率分辨率它只是让频谱曲线看起来更光滑便于读数。频谱结果出来后怎么解读第一步看主频峰在哪个频率。勘探数据的主频一般落在 20–50 Hz 区间如果主频跑到 100 Hz 以上而你处理的是深层反射数据那大概率是噪声或格式误读。第二步看高频衰减趋势如果高频端突然有一个尖锐的凸起可能是 50 Hz 工频干扰或机械谐振噪声。第三步看低频端如果 0–5 Hz 附近能量异常高可能就是数据漂移或去趋势不彻底。另外如果某道频谱和其他道的频谱明显不同比如整体噪声电平抬高说明这道可能是坏道或受到强烈干扰。对于多道数据只分析单道样本不够全面我通常还会做平均振幅谱。对所有道逐一计算振幅谱然后取中位数不是均值避免坏道拉高得到的平均谱反映整炮数据的总体频带特征。这能帮助判断全炮数据质量是否一致也能作为后续时频分析或滤波参数设计的依据。3.4 关键参数的计算过程采样率、奈奎斯特频率与频率分辨率很多初学者容易在采样率和奈奎斯特频率的换算上栽跟头。项目里常见的二进制卷头采样间隔单位是微秒比如常见的 2000 表示 2000 μs也就是 2 ms。换算成采样率采样率 fs 1 / 0.002 s 500 Hz奈奎斯特频率f_nyquist fs / 2 250 Hz这意味着频谱图横轴最高只能画到 250 Hz画到 500 Hz 反而会误导因为 Nyquist 以上全是镜像伪频。频率分辨率则取决于参与 FFT 的信号实际长度df fs / NN 是参与 FFT 的实际样本数补零前比如 N1000、fs500 Hz则 df 0.5 Hz。如果你需要分辨相隔 2 Hz 的两个相邻谱峰就必须保证 df 小于 2 Hz否则峰会被合并成一个大包。还有一个容易忽视的参数是“时窗长度与频谱精度”的关系。对 2 ms 采样的数据如果想得到 1 Hz 的频率分辨率需要的实际信号长度至少是 1 秒也就是 500 个采样点。如果只截取 200 ms 的子段频率分辨率只有 5 Hz。所以在做频谱分析时先确认目标频段内需要分辨的最小间隔再决定截取多长的数据窗这个思路比盲目把整道拿去做 FFT 更严谨。4. 常见问题与排查技巧实录4.1 读取结果显示成“花屏”或杂乱无章的毛刺这个问题我在实际项目里遇到过很多次常见原因有两个。第一是字节序搞错。SEG-Y 头段是 Big-EndianIBM 浮点数据也是 Big-Endian如果误用小端序解析读取出来的数值会非常离谱波形会呈现剧烈的锯齿状或全乱码。排查方法很简单打印前几个原始值对比 IBM 浮点在同一段内存里的预期值比如某个已知的强振幅道前几个样点应该在 -1 到 1 之间如果出现 1e20 之类的天文数字就要检查端序。第二是 IBM 浮点转 IEEE 浮点。IBM 浮点是一种老式计算体系使用的浮点格式它由 1 位符号位、7 位基 16 的指数部分和 24 位尾数组成。很多现代语言没有内置 IBM 浮点转换函数需要自己实现def ibm2ieee(raw_int): if raw_int 0: return 0.0 sign (raw_int 31) 0x01 exponent (raw_int 24) 0x7F mantissa raw_int 0x00FFFFFF if exponent 0 and mantissa 0: return 0.0 value mantissa / float(1 24) exponent_val exponent - 64 value value * (16.0 ** exponent_val) return -value if sign else value好在 segyio 会根据格式码自动做格式码转换所以项目里我直接使用segyio而不手写底层转换但知道原理依然重要——因为总有一些非主流格式码不会被自动处理。4.2 读出来的数据全零或者只有个别道有信号全零数据通常是文件头信息与实际存储不一致。最常见的情况是采样点数读错如果实际道里每个采样点占用字节数超过了你配置的“每道长度”读取会错位导致大部分数据都落在错误的位置看起来就像是全零或大片无效值。遇到这种情况我会先用一个十六进制查看器翻一下文件中某一道数据的起始位置看看前几个字节是不是可解析的浮点数如果看起来不是多半是格式码或道头偏移有问题。还有一种情况是 SEG-Y 文件不是标准的一段数据紧跟一个道头而是先把所有道头批量存储再存所有道体这是 SEG-Y Rev2 和非标准扩展中可能出现的存储布局。用默认顺序读取时会把道头当成数据体解释出来的自然全是“噪声”。排查时用segyio.tracecount和segyio.samples对比文件大小快速估算文件总字节数与“道数 × 每道字节数”是否吻合能在几秒钟内定位这种布局问题。4.3 频谱图“毛刺太多”主峰不突出如果频谱画出来后像一把锯齿梳子而不是平滑的峰型首先检查是否做了去均值和加窗。没去均值时频谱零频附近会出现很大的直流分量没加窗或窗口太窄时谱泄漏会让谱线出现周期性的波动。把汉宁窗加上、把窗长度选到信号主频周期的 3–5 倍以上毛刺会明显减少。此外如果信号本身信噪比低比如含有大量随机噪声频谱基线抬高是很正常的不要误判为“有两个主频峰”最好同时画出多道平均谱来做对比。还有一个容易被忽略的细节多道平均谱要用振幅谱的平方功率谱平均而不是直接对振幅谱做算术平均。因为振幅谱的平均会受相位差异影响功率谱平均能更稳定地体现能量分布。我通常用scipy.signal.welch对多段数据做 Welch 平均功率谱估计然后用它的平方根转换为振幅谱。这样得到的谱线更平滑比单次 FFT 的结果稳定得多。4.4 大数据体读取太慢或内存直接爆掉一个三维观测的 SEG-Y 文件可能高达几十 GB直接用f.trace.raw[:]一次性读取内存开销会非常夸张。正确的做法是用内存映射或逐道读取with segyio.open(file_path, r, strictFalse) as f: for i in range(f.tracecount): trace_i f.trace.raw[i] # 对单道做分析、显示或保存结果segyio在内部使用内存映射memory-mapped filef.trace[i]只加载需要的那一道速度很快。对于频谱分析这种需要全道统计的任务可以分块读取比如每次读 1000 道处理完就释放避免一次性把整个数据体塞进内存。另外如果只需要做显示可以只抽取部分道做抽稀比如每 5 道取 1 道对宏观认识整条剖面足够用了。项目里我还会把解析出来的二进制卷头、道头属性先缓存成 JSON 或 numpy 数组避免每次处理都重新解析整个文件。这样做之后二次处理时只读数据体部分速度提升非常明显。5. 进一步扩展时频分析与数据质控频谱分析解决的是“哪些频率有能量”以及“能量强弱如何”的问题但如果想知道不同时间段的频率变化就需要用时频分析。地震波传播过程中高频成分衰减更快导致浅层剖面的主频较高、深层主频较低这种“随时间变化”的频谱特征在平均频谱上看不出来只有通过短时傅里叶变换STFT或小波变换才能展现。我建议在这个项目里追加一个简单的 STFT 模块用scipy.signal.stft就可以实现f, t, Zxx signal.stft(data, fs500.0, nperseg256, noverlap200)然后把时频谱用pcolormesh画出来就能看到不同时间窗内的频率能量分布。在地震数据质控中这尤其有用比如浅层强反射在时频谱上应该呈现集中在 30–60 Hz 的能量条带如果某个时间段出现了宽频高频能量说明可能有面波噪声或异常干扰如果 50 Hz 处始终存在一条亮线基本可以确定是工频干扰后续滤波设计就有了明确目标。另外把频谱分析结果做成“批量质控报告”非常实用。项目里我会写一个循环对每一炮数据提取平均主频、有效频带宽度、信噪比估计然后把这些指标汇总成一张表格用于快速扫描整个数据集的品质分布。这个思路比一炮一炮人工看波形要高效得多尤其适合资料处理前的大规模数据筛选。对于新手先把单道频谱分析理解透再逐步扩展到这种批量质控就不会被海量数据淹没。6. 写在最后的几点实操体会我做了多年地震数据处理回头来看segy_Read 这类“底层工具”其实比花哨的成像算法更考验细节处理能力。数据格式没吃透、道头参数摸不清后面的步骤做得再漂亮也是空中楼阁。所以我一直跟身边的新人强调拿到一个 SEG-Y 文件先别急着灌进商业软件自己写段脚本打印一下文本卷头、二进制卷头和各道道的道头统计把这些“元信息”通读一遍你会对数据有完全不一样的理解。另外处理过程中要养成“随时保存中间结果”的习惯。比如解析出来的头段信息、归一化参数、频谱指标在每一步都留下可追溯的记录。这样一旦后面出了问题可以快速定位是哪个环节引入的误差。我在项目中还经常会加一个“数据自检”函数检查采样点数是否一致、道间振幅是否存在系统性突变、频谱主频是否在预期范围内任何一项不通过都会给出警告。这个小习惯帮我省下了大量排查无效数据的精力。如果你刚开始接触 SEG-Y我建议先用自己程序生成的模拟数据练手因为你知道它的真实主频和道头内容可以验证自己的代码是否正确。等脚本稳定了再用真实野外数据测试遇到格式变体的时候你会有底气去排查问题。最后再分享一个小技巧绘制频谱图时横轴不要直接线性坐标到底可以同时画线性坐标和对数坐标两幅图。线性坐标便于读出主频峰对数坐标便于观察高频衰减趋势两者搭配会让你的数据质量判断快一个量级。本文还有配套的精品资源点击获取