尧图网页设计 
简介面向地震勘探数据处理需求这套基于MATLAB的工具包聚焦SEGY格式地震数据的读取、显示与频谱分析适合地震数据处理人员、科研人员及高校相关专业学生使用。压缩包内共6个文件全部为.m脚本整体体积仅4KB包含了从SEGY文件头与道数据解析、波形绘制、IBM浮点格式转换到数据预处理和频谱计算等多项功能构成了一套轻量级且可扩展的基础处理流程。目前已吸引557人学习说明该资源在实践教学与项目研究中具备一定参考价值。借助这些脚本用户可将二进制SEGY数据高效解析为矩阵形式通过Wiggle图直观查看地震道振幅变化利用频谱分析识别地下介质的频率响应进而为地质构造解释、岩性识别以及后续去噪、标定等环节提供有效支撑有助于提升地震资料解释的效率和可靠性。1. segy 读取不只是打开文件从 segyread 到频谱分析30 分钟跑通最小流程拿到一块工区的 SEGY 炮集文件很多人的第一反应是赶紧画波形图。真正做过资料品质控制的人知道最花时间的反而是第一关把 segy 文件读进来。segyread这个动作看起来只是「读取」实际上它决定了后面所有显示和频谱分析结果是否可信。标题把「segy 读取显示」和「频谱分析」放在一起正好点出了地震勘探数据处理链条上最常见的诉求把磁带或磁盘上的二进制道数据变成能看的炮集图、能算的振幅谱然后判断这套资料能不能用、频带够不够。本文按「格式布局 → 读取 → 显示 → 频谱 → 排查」的顺序展开写给需要自己动手吃透 SEGY 的从业者。新手能照着跑熟手能拿去核对参数边界。2. SEGY 二进制布局卷头、道头和道数据segyread 第一步都在这里错位SEGY 文件不像 CSV 那样一行一个逗号。它是一串连续二进制先文件头再逐道「道头 道数据」。文件头一共 3600 字节其中 3200 字节是文本卷头很多是 EBCDIC 编码打印出来像乱码但本来就不是设计给人直接读的接着 400 字节二进制卷头里面写着采样间隔、采样点数、数据格式码。再往后的每个道固定 240 字节道头后面紧跟道数据。segyread 要干的本质工作就是按这个布局一步步走。很多人把 SEGY 当普通二进制直接np.fromfile出来的 shape 天然是错的——因为道头和数据是交错的不是纯数值矩阵。2.1 SEGY 的三级结构3200 字节文本卷头、400 字节二进制卷头与 240 字节道头先建立坐标系。SEG-YRev1 标准的文件结构可以用一张表说清楚这张表也是后面所有代码的偏移依据位置字段名字节数说明文件头偏移 0-3199文本卷头3200EBCDIC/ASCII记录工区名、观测系统等文件头偏移 3200-3599二进制卷头400采样间隔、采样点数、格式码等二进制卷头偏移 17-18采样间隔2单位微秒常见 20002 ms、40004 ms二进制卷头偏移 21-22采样点数2每道采样点数如 1500、3000二进制卷头偏移 25-26数据格式码21IBM 浮点232 位定点316 位定点5IEEE 浮点道头偏移 1-4道序列号4全文件内道顺序第一道应为 1道头偏移 21-24CDP 号4注意部分老资料此字段为 0道头偏移 37-40偏移距4带符号整数单位米道头偏移 115-116采样点数2该道实际采样点数可与卷头不一致这里最容易犯的错是把「3200 字节文本卷头」当成整个文件头。读文件时只跳 3200 字节结果 400 字节的二进制卷头被当成第一道道头后面的字段全部错位。另一个高频坑是格式码老磁带资料大多是 1IBM 浮点这种格式的内部表达是「符号 尾数 基 16 指数」和现代 IEEE 浮点完全不同直接按 4 字节浮点读出来就是天文数字。如果你在代码里看到振幅 1e30 之类先查格式码别改算法。2.2 segyread 的三种读法ObsPy、Seismic Unix 与自写解析不同场景选不同读取路线我一般按下面三选一路线适用场景主要限制ObsPyread(file.sgy)快速验证、绘图、单文件分析整文件进内存超大文件吃力Seismic Unixsegyread命令批量处理、管道流水线、老资料需要编译环境头字段固定按标准表自写struct解析特殊格式码、深度学习数据入口偏移错一个后面全错位ObsPy 的好处是把道头字段解析成 dict直接按名字取Seismic Unix 则适合把segyread、sugain、suxwigb、sufft串成管道炮集批量处理非常顺手。自写解析最可控尤其当你需要把数据喂给 PyTorch 或 TensorFlow 时没必要引入头部字段到 numpy 的中间转换。但自写解析要求你对上一小节的偏移表烂熟于心没把握时先用 ObsPy 打印一份道头字典对照。2.3 最小读取验证先打印道头再谈频谱不管用哪条路线读进来的第一件事不是画图而是验证道头。下面的脚本用 ObsPy 读取并打印第一个道的关键信息from obspy import read st read(shot1.sgy) tr st.traces[0] print(tr.stats) # delta, npts, sampling rate h tr.stats.segy.trace_header print(tracl , h[tracl]) # 道序列号第一道应为 1 print(fldr , h[fldr]) # 原始记录号即炮号 print(cdp , h[cdp]) # CDP 号 data tr.data.astype(float) # FFT 之前先转浮点 print(采样点数:, tr.stats.npts, 采样间隔(s):, tr.stats.delta) print(道数据前 5 个样点:, data[:5])逻辑说明tr.stats.delta已经换算成秒是 ObsPy 做好的换算h[tracl]是文件内道序号第一道应为 1如果打出来是 0 或乱码说明卷头跳读有问题。fldr是炮号用于确认文件里的道是按炮组织还是按道组织。参数说明真正跑批量前建议连续打印前 3 道的tracl和fldr看是否按预期递增。如果tracl1,2,3但fldr乱跳说明炮集文件里道排列不是按采集顺序后面显示前必须先排序。再给一个不依赖 ObsPy 的自写解析版本帮助你理解偏移的物理意义import struct with open(shot1.sgy, rb) as f: f.seek(3200) # 跳过 3200 字节文本卷头 bin_header f.read(400) # 400 字节二进制卷头 dt_us struct.unpack(h, bin_header[17:19])[0] # 采样间隔微秒 ns struct.unpack(H, bin_header[21:23])[0] # 采样点数 fmt struct.unpack(h, bin_header[25:27])[0] # 数据格式码 print(dt(us) , dt_us, ns , ns, format , fmt) f.seek(3600) # 文本卷头 二进制卷头 tr_header f.read(240) # 第一道道头 tracl struct.unpack(i, tr_header[0:4])[0] print(first trace tracl , tracl)逻辑说明SEGY 普遍是大端字节序所以struct.unpack用前缀。bin_header[17:19]对应二进制卷头第 17-18 字节换算到文件绝对位置就是 3217-3218。第二段跳转到 3600正好跳过完整文件头然后读 240 字节道头取前 4 字节的道序列号。参数说明h是 2 字节有符号短整数H是无符号短整数采样点数用无符号更稳i是 4 字节有符号整数。如果你把i写成h读到高位字节时会得到一个完全错误但看起来「合理」的数字这是自写解析最阴险的错误。提示道头第 115-116 字节有该道自己的采样点数标准允许它覆盖卷头值。读多道时每道的实际长样点数以道头为准不要全信文件头这是老资料最常见的隐藏差异。3. segy 读取显示把炮集变成能看的波形与变面积图读对了之后显示是另一道坎。勘探软件里的「显示」不只是把数据画出来而是把同一炮几十上百道按道排开让人肉眼能看出同相轴是否连续、初至是否正常、有没有坏道。决定成败的两个参数是增益和道顺序。增益不对深层的弱反射信号全被浅层强能量压没道顺序不对同相轴呈锯齿状再好的数据看起来都是坏的。3.1 波形与变面积显示增益没调对之前图全是假的波形图wiggle适合看单道振幅细节变面积图variable area把正振幅涂黑适合看同相轴横向连续性。工业软件默认两种显示同时开自己在本地跑用 Seismic Unix 一行就能出图segyread tapeshot1.sgy shot1.su suxwigb shot1.su perc97 titleShot 1逻辑说明segyread把 SEGY 读成 SU 内部二进制格式suxwigb读取并弹出波形窗口。perc97表示按 97% 分位的振幅做归一化而不是按最大绝对值这样个别异常道不会把整个图的显示范围拽坏。参数说明perc是这类显示工具里最重要的参数取值建议 95-99 之间。97 适用于大多数正常采集资料如果工区里有明显的强噪声道调到 95 能压住干扰但弱反射信息也会同步变淡。无图形界面的服务器上跑suxwigb会报 X11 错误这种环境建议直接用 Python 输出 PNG。3.2 道头里的排列陷阱CDP、道号与偏移距的顺序问题炮集文件里的道不一定是按偏移距排好的。有的按接收道号有的按采集顺序有的中间混了废道。直接画会出现横向道序来回跳同相轴呈锯齿状。我一般先读道头偏移距字段按偏移距升序重排这是带项目时的血泪经验import struct import numpy as np from obspy import read st read(shot1.sgy) # 直接从原始文件同步读道头偏移距避免不同库字段名差异 with open(shot1.sgy, rb) as f: f.seek(3200 400) # 跳到第一个道头 offsets [] for tr in st.traces: tr_header f.read(240) offsets.append(struct.unpack(i, tr_header[36:40])[0]) ns tr.stats.npts f.seek(ns * 4, 1) # 跳过当前道数据4 字节/样点 order np.argsort(offsets) st.traces [st.traces[i] for i in order] print(排序后前 10 个偏移距(米):, [offsets[i] for i in order[:10]])逻辑说明这段代码保持 ObsPy 对象不变只调整了st.traces的顺序。文件指针按「240 字节道头 ns×4 字节数据」步进每道的偏移距从道头第 37-40 字节读出。排序后打印前 10 个偏移距应该呈现单调递增或递减。参数说明ns * 4是假定每个样点 4 字节。如果格式码是 316 位定点要改成ns * 2如果格式码是 1 或 54 字节没问题。走通后可以把offsets和cdp一起打印核对观测系统定义是否与道头一致。3.3 单炮显示最小脚本加载、抽取、归一化与输出在 Python 里做变面积风格的炮集显示最快的方式是用imshow把道集矩阵画成灰度/彩色图import matplotlib.pyplot as plt import numpy as np from obspy import read st read(shot1.sgy) n_traces len(st.traces) dt st.traces[0].stats.delta # 秒 npts st.traces[0].stats.npts # 堆成矩阵行 时间列 道 matrix np.stack([tr.data.astype(float) for tr in st.traces], axis1) matrix matrix / np.max(np.abs(matrix)) # 归一化到 [-1, 1] plt.figure(figsize(12, 6)) plt.imshow(matrix, aspectauto, cmapseismic, extent[0, n_traces, npts * dt, 0]) plt.xlabel(Trace number) plt.ylabel(Time (s)) plt.title(Shot gather - variable area style) plt.colorbar(labelNormalized amplitude) plt.tight_layout() plt.savefig(shot1_gather.png, dpi150)逻辑说明np.stack(..., axis1)把n_traces条长度npts的一维地震道组成二维矩阵行是时间列是道。extent把横轴映射为道号、纵轴映射为时间秒。注意imshow纵轴默认从上往下所以extent里 y 范围写[npts*dt, 0]而不是[0, npts*dt]这样 t0 在顶部符合地震显示习惯。参数说明cmapseismic红蓝代表正负振幅类似变面积图的极性想更接近工业软件的纯黑白风格可以换成cmapgray并自己控制填充阈值。dpi150适合屏幕查看和报告插图正式出版级图建议 300。振幅差异大的炮集先做归一化还不够需要用滑动窗口 AGC否则深层信号依然看不见。4. 频谱分析把地震道从时间域变到频率域主频和频带一眼看穿频谱分析在勘探里干两件事看有效频带、看主频。采集资料的激发能量、检波器耦合、环境噪声、吸收衰减都会在振幅谱上留下痕迹。FFT 本身不是难点难在输入数据是否干净——直接从 segyread 拿到的原始道几乎不能直接做 FFT。4.1 FFT 之前的三件事去均值、去趋势与选时窗直接对整道做 FFT会看到三个典型假象0 Hz 附近的直流尖峰、端点不连续导致的频谱泄漏、低频段能量被噪声污染。原因分别是均值不为零、时域首尾振幅不连续、低速噪声能量集中在最低频。处理办法固定三步去均值、去趋势、加时窗。import numpy as np def trace_spectrum(data, dt, taper0.05): 单道振幅谱。 data: 1D float 数组 dt: 采样间隔单位秒 taper: 两端 cosine 渐变比例 data data - np.mean(data) # 去均值 n len(data) # 两端 cosine taper保留主体能量 if taper 0: n_taper int(n * taper) x np.linspace(0, 1, n_taper) edge_win 0.5 * (1 - np.cos(np.pi * x)) win np.ones(n) win[:n_taper] edge_win win[-n_taper:] edge_win[::-1] data data * win spec np.fft.rfft(data) freqs np.fft.rfftfreq(n, ddt) return freqs, np.abs(spec)逻辑说明去均值用data - np.mean(data)把直流分量清掉taper 用余弦渐变而不是整道乘汉宁窗是为了只软化首尾、保留主体信号能量。rfft只输出 0 到奈奎斯特频率的谱线对实信号足够计算量减半。参数说明dt的单位必须是秒。SEGY 头部里存的是微秒ObsPy 已经转好但自写解析时如果直接把dt_us传进来频率轴会缩小 100 万倍这是频谱分析里最常见的单位翻车点。taper建议 0.03-0.1短道或信号占满整道的记录用 0.03。4.2 振幅谱与相位谱从 segyread 读出的数据里能看到什么单道振幅谱抖动很大看整体频带必须做多道平均from obspy import read import numpy as np st read(shot1.sgy) dt st.traces[0].stats.delta npts st.traces[0].stats.npts freqs np.fft.rfftfreq(npts, ddt) amp_sum np.zeros_like(freqs) for tr in st.traces: d tr.data.astype(float) d d - d.mean() # 去均值 spec np.fft.rfft(d * np.hanning(npts)) amp_sum np.abs(spec) avg_amp amp_sum / len(st.traces) peak_idx 1 np.argmax(avg_amp[1:]) # 跳过 0 Hz print(平均主频: %.2f Hz % freqs[peak_idx])逻辑说明每道先做均值去除再乘汉宁窗然后累加振幅谱最后除以道数得到平均振幅谱。主频的取法是跳过 0 Hz 后找最大值避免直流分量抢走峰值。参数说明np.hanning(npts)是整道加窗与前面的 taper 思路不同——整道加窗适合多道平均时压制边瓣代价是两端信号权重低如果你关注初至附近的波形频谱应该用 4.1 的自定义 taper。单道谱主频没有统计意义至少 20 道平均才稳定。振幅谱主要看三点主频位置、-6 dB 带宽、陷波点。比如 50 Hz 附近出现明显下凹基本是工业电干扰主频明显低于工区设计值说明激发或近地表吸收有问题。相位谱在常规 QC 里用得少主要在做最小相位化或仪器响应校正前看起跳点而且 FFT 裸算出来的相位是折叠的不能直接统计必须先做 unwrap。4.3 频谱图与主频统计一套可复现的品质参数单炮的频谱要往下游用不能只说「看起来还行」。我习惯把每一炮的主频和有效带宽落成数字这样横向对比才有依据lo, hi 3.0, 100.0 # 工区有效频带上下限 mask (freqs lo) (freqs hi) band_amp avg_amp[mask] band_freqs freqs[mask] pk band_freqs[np.argmax(band_amp)] half band_amp.max() / 2 idx_half np.where(band_amp half)[0] if idx_half.size 2: bw band_freqs[idx_half[-1]] - band_freqs[idx_half[0]] else: bw 0.0 print(带内主频 %.1f Hz, -6dB 带宽 %.1f Hz % (pk, bw))逻辑说明mask先限定有效频带避免地滚波和直流附近的能量把峰值引到最低频。带内主频是有效频带里的最大振幅对应频率-6dB 带宽取振幅谱大于主峰一半的频率范围宽度。参数说明lo3, hi100是典型的陆上可控震源资料范围。低频端 3 Hz 以下的地滚波能量会严重干扰主频统计高频端到奈奎斯特的 70% 左右比较保险4 ms 采样奈奎斯特 125 Hz取 100 Hz 刚好避开高频噪声尾巴。不同工区按资料频率特征调整这两个值即可。5. 读取与频谱分析避坑从道头错位到频率轴翻倍的 5 个典型翻车现场前面三章把主流程走通真正到资料上才是考验。下面 5 个坑是我在项目和工区带教里反复见到的按现象、原因、解决三层写方便排查时对号入座。5.1 道头和格式解析的坑格式码误判与卷头跳读坑一IBM 浮点被当 IEEE 读。现象是振幅变成 1e30 级别的天文数字波形显示全是毛刺频谱完全不可信。原因是老磁带资料的格式码经常是 1IBM 浮点它的内部表达是「符号 尾数 基 16 指数」和现代 IEEE 浮点完全两套规则。解决方法是先读二进制卷头第 25-26 字节拿到格式码遇到 1 先做 IBM 到 IEEE 的转换。ObsPy 对常用格式码会处理但自写解析时这是必踩点转换算法有标准 C 实现可以直接移植别手工猜公式。坑二卷头只跳 3200 字节。现象是第一道道头字段错乱tracl不是 1 而是负数或大数后续所有道全部错位。原因是文件头实际是 3200 字节文本卷头加 400 字节二进制卷头共 3600 字节代码里只跳了文本卷头。解决方法是把跳读位置改成 3600然后打印第一道tracl验证是否为 1。另一个隐藏版本是有些老资料带扩展文本卷头文件头比 3600 更大此时先打印文件前几 KB 判断卷头实际长度再决定跳读字节数。5.2 频谱与显示异常的坑直流塔、频率轴翻倍与道序断裂坑三0 Hz 附近一座直流塔把主频压没。现象是频谱图左侧顶点冲天带内主频完全看不清。原因是数据均值不为 0常见来源是检波器零漂、去噪工具残留的常数项、或时窗起止点振幅不连续。解决方法是 FFT 前先去均值再看时域波形尾部有没有台阶或斜线有就先把尾巴切掉或加 taper。这个坑在短记录和高低频能量悬殊的资料上尤其严重。坑四频率轴整体翻倍。现象是工区资料主频应该在 30 Hz 附近算出来 60 Hz或者 4 ms 采样的资料横轴显示到 250 Hz而奈奎斯特只有 125 Hz。原因是采样间隔单位没转对最常见的是把道头里的采样点数ns误当成采样率或者把 4000 微秒的采样间隔在代码里写死成 2000 微秒。解决方法是打印dt_us确认工区设计采样率2 ms 就是 20004 ms 就是 4000传给rfftfreq之前先除以 1e6 转成秒。我在老资料上吃过一次这个亏主频标到整一倍排查两天最后发现是单位问题不是算法问题。坑五显示出来的炮集同相轴呈锯齿状。现象是波形连续但横向道序来回跳初至和折射波完全对不上。原因是 SEGY 道在文件里的排列不是按观测系统的炮检距顺序有的按接收道号有的混了废道。解决方法是读取道头第 37-40 字节的偏移距字段排序后再显示。排序后如果偏移距仍然不是单调变化就要怀疑观测系统定义和道头字段对应关系有问题这事比显示更严重会影响后续静校正和速度分析。6. 把 segyread 接进采集质量检查批量频谱 QC 与结果验证单文件跑通之后下一步是批量。一条测线几百炮手工挨个画图不现实。我把 segyread、频谱分析和 CSV 输出串成一个 QC 脚本一炮一分钟能跑完。6.1 批量 QC 最小实现读 SEGY、算主频、落 CSVimport csv import numpy as np from obspy import read def qc_one_shot(path, lo3.0, hi100.0): st read(path) dt st.traces[0].stats.delta npts st.traces[0].stats.npts freqs np.fft.rfftfreq(npts, ddt) amp np.zeros_like(freqs) for tr in st.traces: d tr.data.astype(float) d d - d.mean() d d * np.hanning(npts) amp np.abs(np.fft.rfft(d)) amp / len(st.traces) mask (freqs lo) (freqs hi) pk freqs[mask][np.argmax(amp[mask])] energy float(np.sum(amp[mask] ** 2)) return pk, energy results [] for shot_file in [shot_001.sgy, shot_002.sgy, shot_003.sgy]: pk, energy qc_one_shot(shot_file) results.append([shot_file, round(pk, 2), round(energy, 2)]) with open(qc_freq.csv, w, newline) as f: w csv.writer(f) w.writerow([shot, peak_freq_hz, band_energy]) w.writerows(results)逻辑说明qc_one_shot对每一炮做多道平均振幅谱输出带内主频和带内能量。主频相邻炮之间跳跃超过 3-5 Hz说明震源激发一致性有问题能量异常高或异常低说明存在强噪声道、坏道或废炮需要回看单炮图。参数说明lo, hi沿用 4.3 的有效频带设置同一工区全测线用同一组值数值才能横向比较。汉宁窗会压低首尾振幅但 QC 关心的是炮间相对比较同一窗函数下结论是稳定的。验证方法要在真实资料之前做生成一个已知主频 30 Hz 的 Ricker 子波采样间隔 2 ms喂给qc_one_shot看返回是否接近 30 Hz。如果对不上先查dt单位再查时窗处理。这一步骤能救下整条测线的 QC 结论——频谱分析代码一旦有 bug所有炮的主频都是错的事后很难回溯。最后一件事每次 segyread 完第一件事就是把dt_us和ns打进日志。我早年在一批 4 ms 老资料上把采样间隔头读错频谱主频标到整一倍追了两天最后定位在单位换算。从那以后我坚持让任何一次分析都能回溯到最原始的头部字段。这个习惯比任何调试技巧都管用希望帮到你。本文还有配套的精品资源点击获取 SEGY读取显示与频谱分析：从道头解析到批量QC最小流程
简介面向地震勘探数据处理需求这套基于MATLAB的工具包聚焦SEGY格式地震数据的读取、显示与频谱分析适合地震数据处理人员、科研人员及高校相关专业学生使用。压缩包内共6个文件全部为.m脚本整体体积仅4KB包含了从SEGY文件头与道数据解析、波形绘制、IBM浮点格式转换到数据预处理和频谱计算等多项功能构成了一套轻量级且可扩展的基础处理流程。目前已吸引557人学习说明该资源在实践教学与项目研究中具备一定参考价值。借助这些脚本用户可将二进制SEGY数据高效解析为矩阵形式通过Wiggle图直观查看地震道振幅变化利用频谱分析识别地下介质的频率响应进而为地质构造解释、岩性识别以及后续去噪、标定等环节提供有效支撑有助于提升地震资料解释的效率和可靠性。1. segy 读取不只是打开文件从 segyread 到频谱分析30 分钟跑通最小流程拿到一块工区的 SEGY 炮集文件很多人的第一反应是赶紧画波形图。真正做过资料品质控制的人知道最花时间的反而是第一关把 segy 文件读进来。segyread这个动作看起来只是「读取」实际上它决定了后面所有显示和频谱分析结果是否可信。标题把「segy 读取显示」和「频谱分析」放在一起正好点出了地震勘探数据处理链条上最常见的诉求把磁带或磁盘上的二进制道数据变成能看的炮集图、能算的振幅谱然后判断这套资料能不能用、频带够不够。本文按「格式布局 → 读取 → 显示 → 频谱 → 排查」的顺序展开写给需要自己动手吃透 SEGY 的从业者。新手能照着跑熟手能拿去核对参数边界。2. SEGY 二进制布局卷头、道头和道数据segyread 第一步都在这里错位SEGY 文件不像 CSV 那样一行一个逗号。它是一串连续二进制先文件头再逐道「道头 道数据」。文件头一共 3600 字节其中 3200 字节是文本卷头很多是 EBCDIC 编码打印出来像乱码但本来就不是设计给人直接读的接着 400 字节二进制卷头里面写着采样间隔、采样点数、数据格式码。再往后的每个道固定 240 字节道头后面紧跟道数据。segyread 要干的本质工作就是按这个布局一步步走。很多人把 SEGY 当普通二进制直接np.fromfile出来的 shape 天然是错的——因为道头和数据是交错的不是纯数值矩阵。2.1 SEGY 的三级结构3200 字节文本卷头、400 字节二进制卷头与 240 字节道头先建立坐标系。SEG-YRev1 标准的文件结构可以用一张表说清楚这张表也是后面所有代码的偏移依据位置字段名字节数说明文件头偏移 0-3199文本卷头3200EBCDIC/ASCII记录工区名、观测系统等文件头偏移 3200-3599二进制卷头400采样间隔、采样点数、格式码等二进制卷头偏移 17-18采样间隔2单位微秒常见 20002 ms、40004 ms二进制卷头偏移 21-22采样点数2每道采样点数如 1500、3000二进制卷头偏移 25-26数据格式码21IBM 浮点232 位定点316 位定点5IEEE 浮点道头偏移 1-4道序列号4全文件内道顺序第一道应为 1道头偏移 21-24CDP 号4注意部分老资料此字段为 0道头偏移 37-40偏移距4带符号整数单位米道头偏移 115-116采样点数2该道实际采样点数可与卷头不一致这里最容易犯的错是把「3200 字节文本卷头」当成整个文件头。读文件时只跳 3200 字节结果 400 字节的二进制卷头被当成第一道道头后面的字段全部错位。另一个高频坑是格式码老磁带资料大多是 1IBM 浮点这种格式的内部表达是「符号 尾数 基 16 指数」和现代 IEEE 浮点完全不同直接按 4 字节浮点读出来就是天文数字。如果你在代码里看到振幅 1e30 之类先查格式码别改算法。2.2 segyread 的三种读法ObsPy、Seismic Unix 与自写解析不同场景选不同读取路线我一般按下面三选一路线适用场景主要限制ObsPyread(file.sgy)快速验证、绘图、单文件分析整文件进内存超大文件吃力Seismic Unixsegyread命令批量处理、管道流水线、老资料需要编译环境头字段固定按标准表自写struct解析特殊格式码、深度学习数据入口偏移错一个后面全错位ObsPy 的好处是把道头字段解析成 dict直接按名字取Seismic Unix 则适合把segyread、sugain、suxwigb、sufft串成管道炮集批量处理非常顺手。自写解析最可控尤其当你需要把数据喂给 PyTorch 或 TensorFlow 时没必要引入头部字段到 numpy 的中间转换。但自写解析要求你对上一小节的偏移表烂熟于心没把握时先用 ObsPy 打印一份道头字典对照。2.3 最小读取验证先打印道头再谈频谱不管用哪条路线读进来的第一件事不是画图而是验证道头。下面的脚本用 ObsPy 读取并打印第一个道的关键信息from obspy import read st read(shot1.sgy) tr st.traces[0] print(tr.stats) # delta, npts, sampling rate h tr.stats.segy.trace_header print(tracl , h[tracl]) # 道序列号第一道应为 1 print(fldr , h[fldr]) # 原始记录号即炮号 print(cdp , h[cdp]) # CDP 号 data tr.data.astype(float) # FFT 之前先转浮点 print(采样点数:, tr.stats.npts, 采样间隔(s):, tr.stats.delta) print(道数据前 5 个样点:, data[:5])逻辑说明tr.stats.delta已经换算成秒是 ObsPy 做好的换算h[tracl]是文件内道序号第一道应为 1如果打出来是 0 或乱码说明卷头跳读有问题。fldr是炮号用于确认文件里的道是按炮组织还是按道组织。参数说明真正跑批量前建议连续打印前 3 道的tracl和fldr看是否按预期递增。如果tracl1,2,3但fldr乱跳说明炮集文件里道排列不是按采集顺序后面显示前必须先排序。再给一个不依赖 ObsPy 的自写解析版本帮助你理解偏移的物理意义import struct with open(shot1.sgy, rb) as f: f.seek(3200) # 跳过 3200 字节文本卷头 bin_header f.read(400) # 400 字节二进制卷头 dt_us struct.unpack(h, bin_header[17:19])[0] # 采样间隔微秒 ns struct.unpack(H, bin_header[21:23])[0] # 采样点数 fmt struct.unpack(h, bin_header[25:27])[0] # 数据格式码 print(dt(us) , dt_us, ns , ns, format , fmt) f.seek(3600) # 文本卷头 二进制卷头 tr_header f.read(240) # 第一道道头 tracl struct.unpack(i, tr_header[0:4])[0] print(first trace tracl , tracl)逻辑说明SEGY 普遍是大端字节序所以struct.unpack用前缀。bin_header[17:19]对应二进制卷头第 17-18 字节换算到文件绝对位置就是 3217-3218。第二段跳转到 3600正好跳过完整文件头然后读 240 字节道头取前 4 字节的道序列号。参数说明h是 2 字节有符号短整数H是无符号短整数采样点数用无符号更稳i是 4 字节有符号整数。如果你把i写成h读到高位字节时会得到一个完全错误但看起来「合理」的数字这是自写解析最阴险的错误。提示道头第 115-116 字节有该道自己的采样点数标准允许它覆盖卷头值。读多道时每道的实际长样点数以道头为准不要全信文件头这是老资料最常见的隐藏差异。3. segy 读取显示把炮集变成能看的波形与变面积图读对了之后显示是另一道坎。勘探软件里的「显示」不只是把数据画出来而是把同一炮几十上百道按道排开让人肉眼能看出同相轴是否连续、初至是否正常、有没有坏道。决定成败的两个参数是增益和道顺序。增益不对深层的弱反射信号全被浅层强能量压没道顺序不对同相轴呈锯齿状再好的数据看起来都是坏的。3.1 波形与变面积显示增益没调对之前图全是假的波形图wiggle适合看单道振幅细节变面积图variable area把正振幅涂黑适合看同相轴横向连续性。工业软件默认两种显示同时开自己在本地跑用 Seismic Unix 一行就能出图segyread tapeshot1.sgy shot1.su suxwigb shot1.su perc97 titleShot 1逻辑说明segyread把 SEGY 读成 SU 内部二进制格式suxwigb读取并弹出波形窗口。perc97表示按 97% 分位的振幅做归一化而不是按最大绝对值这样个别异常道不会把整个图的显示范围拽坏。参数说明perc是这类显示工具里最重要的参数取值建议 95-99 之间。97 适用于大多数正常采集资料如果工区里有明显的强噪声道调到 95 能压住干扰但弱反射信息也会同步变淡。无图形界面的服务器上跑suxwigb会报 X11 错误这种环境建议直接用 Python 输出 PNG。3.2 道头里的排列陷阱CDP、道号与偏移距的顺序问题炮集文件里的道不一定是按偏移距排好的。有的按接收道号有的按采集顺序有的中间混了废道。直接画会出现横向道序来回跳同相轴呈锯齿状。我一般先读道头偏移距字段按偏移距升序重排这是带项目时的血泪经验import struct import numpy as np from obspy import read st read(shot1.sgy) # 直接从原始文件同步读道头偏移距避免不同库字段名差异 with open(shot1.sgy, rb) as f: f.seek(3200 400) # 跳到第一个道头 offsets [] for tr in st.traces: tr_header f.read(240) offsets.append(struct.unpack(i, tr_header[36:40])[0]) ns tr.stats.npts f.seek(ns * 4, 1) # 跳过当前道数据4 字节/样点 order np.argsort(offsets) st.traces [st.traces[i] for i in order] print(排序后前 10 个偏移距(米):, [offsets[i] for i in order[:10]])逻辑说明这段代码保持 ObsPy 对象不变只调整了st.traces的顺序。文件指针按「240 字节道头 ns×4 字节数据」步进每道的偏移距从道头第 37-40 字节读出。排序后打印前 10 个偏移距应该呈现单调递增或递减。参数说明ns * 4是假定每个样点 4 字节。如果格式码是 316 位定点要改成ns * 2如果格式码是 1 或 54 字节没问题。走通后可以把offsets和cdp一起打印核对观测系统定义是否与道头一致。3.3 单炮显示最小脚本加载、抽取、归一化与输出在 Python 里做变面积风格的炮集显示最快的方式是用imshow把道集矩阵画成灰度/彩色图import matplotlib.pyplot as plt import numpy as np from obspy import read st read(shot1.sgy) n_traces len(st.traces) dt st.traces[0].stats.delta # 秒 npts st.traces[0].stats.npts # 堆成矩阵行 时间列 道 matrix np.stack([tr.data.astype(float) for tr in st.traces], axis1) matrix matrix / np.max(np.abs(matrix)) # 归一化到 [-1, 1] plt.figure(figsize(12, 6)) plt.imshow(matrix, aspectauto, cmapseismic, extent[0, n_traces, npts * dt, 0]) plt.xlabel(Trace number) plt.ylabel(Time (s)) plt.title(Shot gather - variable area style) plt.colorbar(labelNormalized amplitude) plt.tight_layout() plt.savefig(shot1_gather.png, dpi150)逻辑说明np.stack(..., axis1)把n_traces条长度npts的一维地震道组成二维矩阵行是时间列是道。extent把横轴映射为道号、纵轴映射为时间秒。注意imshow纵轴默认从上往下所以extent里 y 范围写[npts*dt, 0]而不是[0, npts*dt]这样 t0 在顶部符合地震显示习惯。参数说明cmapseismic红蓝代表正负振幅类似变面积图的极性想更接近工业软件的纯黑白风格可以换成cmapgray并自己控制填充阈值。dpi150适合屏幕查看和报告插图正式出版级图建议 300。振幅差异大的炮集先做归一化还不够需要用滑动窗口 AGC否则深层信号依然看不见。4. 频谱分析把地震道从时间域变到频率域主频和频带一眼看穿频谱分析在勘探里干两件事看有效频带、看主频。采集资料的激发能量、检波器耦合、环境噪声、吸收衰减都会在振幅谱上留下痕迹。FFT 本身不是难点难在输入数据是否干净——直接从 segyread 拿到的原始道几乎不能直接做 FFT。4.1 FFT 之前的三件事去均值、去趋势与选时窗直接对整道做 FFT会看到三个典型假象0 Hz 附近的直流尖峰、端点不连续导致的频谱泄漏、低频段能量被噪声污染。原因分别是均值不为零、时域首尾振幅不连续、低速噪声能量集中在最低频。处理办法固定三步去均值、去趋势、加时窗。import numpy as np def trace_spectrum(data, dt, taper0.05): 单道振幅谱。 data: 1D float 数组 dt: 采样间隔单位秒 taper: 两端 cosine 渐变比例 data data - np.mean(data) # 去均值 n len(data) # 两端 cosine taper保留主体能量 if taper 0: n_taper int(n * taper) x np.linspace(0, 1, n_taper) edge_win 0.5 * (1 - np.cos(np.pi * x)) win np.ones(n) win[:n_taper] edge_win win[-n_taper:] edge_win[::-1] data data * win spec np.fft.rfft(data) freqs np.fft.rfftfreq(n, ddt) return freqs, np.abs(spec)逻辑说明去均值用data - np.mean(data)把直流分量清掉taper 用余弦渐变而不是整道乘汉宁窗是为了只软化首尾、保留主体信号能量。rfft只输出 0 到奈奎斯特频率的谱线对实信号足够计算量减半。参数说明dt的单位必须是秒。SEGY 头部里存的是微秒ObsPy 已经转好但自写解析时如果直接把dt_us传进来频率轴会缩小 100 万倍这是频谱分析里最常见的单位翻车点。taper建议 0.03-0.1短道或信号占满整道的记录用 0.03。4.2 振幅谱与相位谱从 segyread 读出的数据里能看到什么单道振幅谱抖动很大看整体频带必须做多道平均from obspy import read import numpy as np st read(shot1.sgy) dt st.traces[0].stats.delta npts st.traces[0].stats.npts freqs np.fft.rfftfreq(npts, ddt) amp_sum np.zeros_like(freqs) for tr in st.traces: d tr.data.astype(float) d d - d.mean() # 去均值 spec np.fft.rfft(d * np.hanning(npts)) amp_sum np.abs(spec) avg_amp amp_sum / len(st.traces) peak_idx 1 np.argmax(avg_amp[1:]) # 跳过 0 Hz print(平均主频: %.2f Hz % freqs[peak_idx])逻辑说明每道先做均值去除再乘汉宁窗然后累加振幅谱最后除以道数得到平均振幅谱。主频的取法是跳过 0 Hz 后找最大值避免直流分量抢走峰值。参数说明np.hanning(npts)是整道加窗与前面的 taper 思路不同——整道加窗适合多道平均时压制边瓣代价是两端信号权重低如果你关注初至附近的波形频谱应该用 4.1 的自定义 taper。单道谱主频没有统计意义至少 20 道平均才稳定。振幅谱主要看三点主频位置、-6 dB 带宽、陷波点。比如 50 Hz 附近出现明显下凹基本是工业电干扰主频明显低于工区设计值说明激发或近地表吸收有问题。相位谱在常规 QC 里用得少主要在做最小相位化或仪器响应校正前看起跳点而且 FFT 裸算出来的相位是折叠的不能直接统计必须先做 unwrap。4.3 频谱图与主频统计一套可复现的品质参数单炮的频谱要往下游用不能只说「看起来还行」。我习惯把每一炮的主频和有效带宽落成数字这样横向对比才有依据lo, hi 3.0, 100.0 # 工区有效频带上下限 mask (freqs lo) (freqs hi) band_amp avg_amp[mask] band_freqs freqs[mask] pk band_freqs[np.argmax(band_amp)] half band_amp.max() / 2 idx_half np.where(band_amp half)[0] if idx_half.size 2: bw band_freqs[idx_half[-1]] - band_freqs[idx_half[0]] else: bw 0.0 print(带内主频 %.1f Hz, -6dB 带宽 %.1f Hz % (pk, bw))逻辑说明mask先限定有效频带避免地滚波和直流附近的能量把峰值引到最低频。带内主频是有效频带里的最大振幅对应频率-6dB 带宽取振幅谱大于主峰一半的频率范围宽度。参数说明lo3, hi100是典型的陆上可控震源资料范围。低频端 3 Hz 以下的地滚波能量会严重干扰主频统计高频端到奈奎斯特的 70% 左右比较保险4 ms 采样奈奎斯特 125 Hz取 100 Hz 刚好避开高频噪声尾巴。不同工区按资料频率特征调整这两个值即可。5. 读取与频谱分析避坑从道头错位到频率轴翻倍的 5 个典型翻车现场前面三章把主流程走通真正到资料上才是考验。下面 5 个坑是我在项目和工区带教里反复见到的按现象、原因、解决三层写方便排查时对号入座。5.1 道头和格式解析的坑格式码误判与卷头跳读坑一IBM 浮点被当 IEEE 读。现象是振幅变成 1e30 级别的天文数字波形显示全是毛刺频谱完全不可信。原因是老磁带资料的格式码经常是 1IBM 浮点它的内部表达是「符号 尾数 基 16 指数」和现代 IEEE 浮点完全两套规则。解决方法是先读二进制卷头第 25-26 字节拿到格式码遇到 1 先做 IBM 到 IEEE 的转换。ObsPy 对常用格式码会处理但自写解析时这是必踩点转换算法有标准 C 实现可以直接移植别手工猜公式。坑二卷头只跳 3200 字节。现象是第一道道头字段错乱tracl不是 1 而是负数或大数后续所有道全部错位。原因是文件头实际是 3200 字节文本卷头加 400 字节二进制卷头共 3600 字节代码里只跳了文本卷头。解决方法是把跳读位置改成 3600然后打印第一道tracl验证是否为 1。另一个隐藏版本是有些老资料带扩展文本卷头文件头比 3600 更大此时先打印文件前几 KB 判断卷头实际长度再决定跳读字节数。5.2 频谱与显示异常的坑直流塔、频率轴翻倍与道序断裂坑三0 Hz 附近一座直流塔把主频压没。现象是频谱图左侧顶点冲天带内主频完全看不清。原因是数据均值不为 0常见来源是检波器零漂、去噪工具残留的常数项、或时窗起止点振幅不连续。解决方法是 FFT 前先去均值再看时域波形尾部有没有台阶或斜线有就先把尾巴切掉或加 taper。这个坑在短记录和高低频能量悬殊的资料上尤其严重。坑四频率轴整体翻倍。现象是工区资料主频应该在 30 Hz 附近算出来 60 Hz或者 4 ms 采样的资料横轴显示到 250 Hz而奈奎斯特只有 125 Hz。原因是采样间隔单位没转对最常见的是把道头里的采样点数ns误当成采样率或者把 4000 微秒的采样间隔在代码里写死成 2000 微秒。解决方法是打印dt_us确认工区设计采样率2 ms 就是 20004 ms 就是 4000传给rfftfreq之前先除以 1e6 转成秒。我在老资料上吃过一次这个亏主频标到整一倍排查两天最后发现是单位问题不是算法问题。坑五显示出来的炮集同相轴呈锯齿状。现象是波形连续但横向道序来回跳初至和折射波完全对不上。原因是 SEGY 道在文件里的排列不是按观测系统的炮检距顺序有的按接收道号有的混了废道。解决方法是读取道头第 37-40 字节的偏移距字段排序后再显示。排序后如果偏移距仍然不是单调变化就要怀疑观测系统定义和道头字段对应关系有问题这事比显示更严重会影响后续静校正和速度分析。6. 把 segyread 接进采集质量检查批量频谱 QC 与结果验证单文件跑通之后下一步是批量。一条测线几百炮手工挨个画图不现实。我把 segyread、频谱分析和 CSV 输出串成一个 QC 脚本一炮一分钟能跑完。6.1 批量 QC 最小实现读 SEGY、算主频、落 CSVimport csv import numpy as np from obspy import read def qc_one_shot(path, lo3.0, hi100.0): st read(path) dt st.traces[0].stats.delta npts st.traces[0].stats.npts freqs np.fft.rfftfreq(npts, ddt) amp np.zeros_like(freqs) for tr in st.traces: d tr.data.astype(float) d d - d.mean() d d * np.hanning(npts) amp np.abs(np.fft.rfft(d)) amp / len(st.traces) mask (freqs lo) (freqs hi) pk freqs[mask][np.argmax(amp[mask])] energy float(np.sum(amp[mask] ** 2)) return pk, energy results [] for shot_file in [shot_001.sgy, shot_002.sgy, shot_003.sgy]: pk, energy qc_one_shot(shot_file) results.append([shot_file, round(pk, 2), round(energy, 2)]) with open(qc_freq.csv, w, newline) as f: w csv.writer(f) w.writerow([shot, peak_freq_hz, band_energy]) w.writerows(results)逻辑说明qc_one_shot对每一炮做多道平均振幅谱输出带内主频和带内能量。主频相邻炮之间跳跃超过 3-5 Hz说明震源激发一致性有问题能量异常高或异常低说明存在强噪声道、坏道或废炮需要回看单炮图。参数说明lo, hi沿用 4.3 的有效频带设置同一工区全测线用同一组值数值才能横向比较。汉宁窗会压低首尾振幅但 QC 关心的是炮间相对比较同一窗函数下结论是稳定的。验证方法要在真实资料之前做生成一个已知主频 30 Hz 的 Ricker 子波采样间隔 2 ms喂给qc_one_shot看返回是否接近 30 Hz。如果对不上先查dt单位再查时窗处理。这一步骤能救下整条测线的 QC 结论——频谱分析代码一旦有 bug所有炮的主频都是错的事后很难回溯。最后一件事每次 segyread 完第一件事就是把dt_us和ns打进日志。我早年在一批 4 ms 老资料上把采样间隔头读错频谱主频标到整一倍追了两天最后定位在单位换算。从那以后我坚持让任何一次分析都能回溯到最原始的头部字段。这个习惯比任何调试技巧都管用希望帮到你。本文还有配套的精品资源点击获取