企业官网建设流程全解析

简介：面向地震勘探数据处理需求，这套基于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-Y（Rev1 标准）的文件结构可以用一张表说清楚，这张表也是后面所有代码的偏移依据：

这里最容易犯的错是把「3200 字节文本卷头」当成整个文件头。读文件时只跳 3200 字节，结果 400 字节的二进制卷头被当成第一道道头，后面的字段全部错位。另一个高频坑是格式码：老磁带资料大多是 1（IBM 浮点），这种格式的内部表达是「符号 + 尾数 + 基 16 指数」，和现代 IEEE 浮点完全不同，直接按 4 字节浮点读出来就是天文数字。如果你在代码里看到振幅 1e30 之类，先查格式码，别改算法。

2.2 segyread 的三种读法：ObsPy、Seismic Unix 与自写解析

不同场景选不同读取路线，我一般按下面三选一：

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，看是否按预期递增。如果tracl=1,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 tape=shot1.sgy > shot1.su suxwigb < shot1.su perc=97 title="Shot 1"

逻辑说明：segyread把 SEGY 读成 SU 内部二进制格式，suxwigb读取并弹出波形窗口。perc=97表示按 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 字节。如果格式码是 3（16 位定点），要改成ns * 2；如果格式码是 1 或 5，4 字节没问题。走通后可以把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], axis=1) matrix = matrix / np.max(np.abs(matrix)) # 归一化到 [-1, 1] plt.figure(figsize=(12, 6)) plt.imshow(matrix, aspect='auto', cmap='seismic', extent=[0, n_traces, npts * dt, 0]) plt.xlabel('Trace number') plt.ylabel('Time (s)') plt.title('Shot gather - variable area style') plt.colorbar(label='Normalized amplitude') plt.tight_layout() plt.savefig('shot1_gather.png', dpi=150)

逻辑说明：np.stack(..., axis=1)把n_traces条长度npts的一维地震道组成二维矩阵，行是时间，列是道。extent把横轴映射为道号、纵轴映射为时间秒。注意imshow纵轴默认从上往下，所以extent里 y 范围写[npts*dt, 0]而不是[0, npts*dt]，这样 t=0 在顶部，符合地震显示习惯。

参数说明：cmap='seismic'红蓝代表正负振幅，类似变面积图的极性；想更接近工业软件的纯黑白风格，可以换成cmap='gray'并自己控制填充阈值。dpi=150适合屏幕查看和报告插图，正式出版级图建议 300。振幅差异大的炮集，先做归一化还不够，需要用滑动窗口 AGC，否则深层信号依然看不见。

4. 频谱分析：把地震道从时间域变到频率域，主频和频带一眼看穿

频谱分析在勘探里干两件事：看有效频带、看主频。采集资料的激发能量、检波器耦合、环境噪声、吸收衰减都会在振幅谱上留下痕迹。FFT 本身不是难点，难在输入数据是否干净——直接从 segyread 拿到的原始道，几乎不能直接做 FFT。

4.1 FFT 之前的三件事：去均值、去趋势与选时窗

直接对整道做 FFT，会看到三个典型假象：0 Hz 附近的直流尖峰、端点不连续导致的频谱泄漏、低频段能量被噪声污染。原因分别是均值不为零、时域首尾振幅不连续、低速噪声能量集中在最低频。处理办法固定三步：去均值、去趋势、加时窗。

import numpy as np def trace_spectrum(data, dt, taper=0.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, d=dt) 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, d=dt) 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 带宽取振幅谱大于主峰一半的频率范围宽度。

参数说明：lo=3, hi=100是典型的陆上可控震源资料范围。低频端 3 Hz 以下的地滚波能量会严重干扰主频统计；高频端到奈奎斯特的 70% 左右比较保险，4 ms 采样奈奎斯特 125 Hz，取 100 Hz 刚好避开高频噪声尾巴。不同工区按资料频率特征调整这两个值即可。

5. 读取与频谱分析避坑：从道头错位到频率轴翻倍的 5 个典型翻车现场

前面三章把主流程走通，真正到资料上才是考验。下面 5 个坑是我在项目和工区带教里反复见到的，按现象、原因、解决三层写，方便排查时对号入座。

5.1 道头和格式解析的坑：格式码误判与卷头跳读

坑一：IBM 浮点被当 IEEE 读。现象是振幅变成 1e30 级别的天文数字，波形显示全是毛刺，频谱完全不可信。原因是老磁带资料的格式码经常是 1（IBM 浮点），它的内部表达是「符号 + 尾数 + 基 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 就是 2000，4 ms 就是 4000，传给rfftfreq之前先除以 1e6 转成秒。我在老资料上吃过一次这个亏，主频标到整一倍，排查两天最后发现是单位问题，不是算法问题。

坑五：显示出来的炮集同相轴呈锯齿状。现象是波形连续但横向道序来回跳，初至和折射波完全对不上。原因是 SEGY 道在文件里的排列不是按观测系统的炮检距顺序，有的按接收道号，有的混了废道。解决方法是读取道头第 37-40 字节的偏移距字段，排序后再显示。排序后如果偏移距仍然不是单调变化，就要怀疑观测系统定义和道头字段对应关系有问题，这事比显示更严重，会影响后续静校正和速度分析。

6. 把 segyread 接进采集质量检查：批量频谱 QC 与结果验证

单文件跑通之后，下一步是批量。一条测线几百炮，手工挨个画图不现实。我把 segyread、频谱分析和 CSV 输出串成一个 QC 脚本，一炮一分钟能跑完。

6.1 批量 QC 最小实现：读 SEGY、算主频、落 CSV

import csv import numpy as np from obspy import read def qc_one_shot(path, lo=3.0, hi=100.0): st = read(path) dt = st.traces[0].stats.delta npts = st.traces[0].stats.npts freqs = np.fft.rfftfreq(npts, d=dt) 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 老资料上把采样间隔头读错，频谱主频标到整一倍，追了两天最后定位在单位换算。从那以后，我坚持让任何一次分析都能回溯到最原始的头部字段。这个习惯比任何调试技巧都管用，希望帮到你。

本文还有配套的精品资源，点击获取