简介一份面向地震勘探与地球物理研究人员的 Matlab SEGY 数据处理工具包。SEGY 是地震数据存储与交换的标准格式SegyMAT 通过读取与解析函数帮助用户在 Matlab 中完成文件头解析、道数据提取与显示大幅降低地震数据读取门槛。压缩包共 75 个文件以 m 源码为主63 个另有 fig 界面文件、dat 与 segy 示例数据以及 license、revision 说明整体仅 516KB轻量易用。内部包含完整源码、图形交互界面、演示脚本及测试数据使用者可直接运行演示程序体验读取流程也可基于底层工具函数进行二次开发实现批量转换、剖面绘制与地震成像等进阶分析。工具包还梳理了 SEGY 头部信息与道头字段的解析逻辑适合用作学习地震数据格式的参考实现。已有 645 人学习下载适合具备一定 Matlab 基础、需要处理 SEGY 数据的地球物理专业学生与工程师。 做地震勘探或者物探数据处理的同学应该都体会过这种尴尬对方发来一个.sgy文件用专业软件打开毫无压力但你想在 MATLAB 里做点个性化的滤波、道集显示或者频谱分析第一关读文件就卡住了。SEGY 这种格式跟普通文本或者 xls 完全不一样它是一个带元数据的二进制块状结构如果对它没有一个基本概念工具包拿到手里也容易踩坑、调半天不知道问题出在哪。这篇内容我把自己用 MATLAB 读 SEGY 的几个常用工具包和完整实操流程整理了一遍同时把最容易翻车的字节序、IBM 浮点问题也单独拎出来讲了看完照着做至少能解决 90% 的常规读取需求。1. SEGY格式老底子搞不懂这个结构工具包用起来也是盲人摸象1.1 一个SEGY文件由三块组成SEGYSEG-Y是地震数据交换标准格式由勘探地球物理学家学会提出。一个标准 SEGY 文件在逻辑上分成三大块文件卷头、扩展文本头可选、道数据块。其中卷头固定占用 3600 字节前 3200 字节通常是 EBCDIC 编码的文本信息记录工区名、观测日期等作业描述后 400 字节是二进制卷头采样点数、采样间隔、数据格式码这些关键参数都固定写在这 400 字节的特定位置上。从第 3601 字节开始就是一条一条的地震道每条道由 240 字节道头加若干个采样点组成采样点可以是 4 字节浮点、整型等多种方式。很多新手拿工具包读出来一堆乱码问题往往出在对这三块的了解不够。比如老标准 SEG-Y rev0 只有一个文本卷头而 rev1 之后可以在标准卷头后面扩展多个 3200 字节的文本卷头生产数据库导出的文件里经常带着一堆扩展文本头。如果你只知道3600 字节之后就是道数据这个固定公式再碰上扩展头就会读错位置。把它想象成一本砌好了目录的书先看目录卷头再逐页翻道每一页上还贴着一张便签道头工具包就是帮你把目录和便签解析出来的助手。1.2 卷头里最常用的字段二进制卷头里字段很多但实际处理中我几乎只关心三个采样点数、采样间隔、数据格式码。这几个字段在 SEG-Y rev1 标准中的字节位置如下字节位置长度含义3217-32182 字节采样间隔单位微秒3221-32222 字节每道采样点数3223-32242 字节数据格式码数据格式码是判断后续数据体怎么读的关键。常见格式码包括1 表示 IBM 浮点2 表示 32 位整型3 表示 16 位整型5 表示 IEEE 浮点8 表示 8 位整型。这个字段决定了你该用float32读还是int32读也决定了老数据那套 IBM 浮点转换要不要上。拿到一个陌生文件第一步永远是把这个格式码看清楚后面 80% 的读取坑都能提前避开。1.3 为什么道头比数据体更容易被忽略道数据的 240 字节道头是另一个重点。道头里保存着每道的位置信息、炮点号、CDP 号、偏移距、采样点数等后续做道集分析、属性提取全靠它。很多人用工具包读回一个数据矩阵后就忘掉了道头结果画图时想标 CDP 坐标又得重新读文件。按 SEG-Y rev1 标准道头内第 21-24 字节是 CDP 号第 37-40 字节是偏移距第 115-116 字节是本道采样点数第 117-118 字节是本道采样间隔。不同处理系统可能会对非标准字段做自定义但这几个标准字段基本是通用的。如果你用的工具包能把这些信息解析成结构体一定要保留下来而不是只拿 Data 矩阵。2. 工具包盘点SEGYMAT、SegMAT和自写方案的取舍2.1 SEGYMAT老牌且轻量SEGYMAT 是 MATLAB 社区里流传最广的 SEGY 读写工具包之一函数接口简单对新手非常友好。核心调用就一行[Data, Seis, TextHeader] ReadSEGY(line.sgy);读取结果里 Data 是一个矩阵每一列是一道Seis 里保存道头信息TextHeader 是文本卷头。它的主要特点是一把梭方便但内存消耗大适合中小规模文件。如果你的文件只有几百 MB 或者稍微上 GB用 SEGYMAT 快速读出来做验证完全没有问题。2.2 SegMAT面向大文件和现代标准SegMAT 是另一个经常被提到的工具包比 SEGYMAT 更注重对 SEG-Y rev1/rev2 的支持而且支持分开读取道头和采样数据遇到几十 GB 的叠前数据不至于一次性把内存打满。它适合在生产数据上做先读道头选道、再按需取数据这种流式处理。不过它的接口命名在不同版本里有一些调整用之前先打开包里的 README 或者 demo 脚本看一眼确认你下载的那个版本具体函数名是什么。2.3 SeisLab读取只是整个地震处理库的一角如果你读完之后还要做带通滤波、反褶积、道均衡、相干体分析这些常规地震处理可以考虑 SeisLab。它由加拿大卡尔加里大学的 CREWES 项目组维护把 SEG-Y 读写纳入了完整的地震处理流程工具箱里有大量现成的地震显示和处理函数。缺点是包体大、依赖多如果单纯只是为了读个文件用它有点杀鸡用牛刀。2.4 自写脚本的适用场景有些时候引包反而不划算。比如文件是某个采集系统自己改造过的道头结构或者你只需要从一千个文件里快速抽取卷头某个字段做质量检查这时候自己用fopenfread写一段小脚本更可控。别觉得写底层读取很难搞清楚字节序和偏移量之后读 SEGY 其实就是一道指针定位 按格式读的算术题。工具包维护活跃度主要特点适合场景SEGYMAT较老函数简洁直接读全文件中小规模数据、快速上手SegMAT较活跃支持分道读取、rev1/rev2大文件、生产数据处理SeisLab较活跃完整地震处理库读写只是其中一部分需要配套滤波、属性分析等整套流程自写脚本—完全可控依赖少特殊结构文件、关键字段检查3. SEGYMAT实操把SEGY读进MATLAB工作区3.1 安装与路径设置在 MATLAB File Exchange 里搜 SEGYMAT 就能找到源码包GitHub 上也有镜像仓库。下载后解压把整个文件夹加入 MATLAB 搜索路径即可。如果只读不写几十个.m文件足够用不需要编译任何 C 代码这也是它受欢迎的原因之一。安装完成后可以先用一个小文件做冒烟测试。我这里以某个陆地工区的单炮记录为例文件不大几万道还没到正好适合 SEGYMAT 一次性读取。3.2 核心调用与结果解析[Data, Seis, TextHeader] ReadSEGY(shot1.sgy); size(Data)从size(Data)能直接看到数据矩阵的维度通常返回的是ns × ntrns 是每道采样点数ntr 是总道数。Seis是道头结构体不同版本字段命名差异比较大比如采样间隔可能是dt也可能写成si或者SampleInt。不要凭记忆去猜读完先用fieldnames(Seis)看一眼有哪些字段再对照你要用的信息去取。有一个新手经常忽略的细节TextHeader是一个字符矩阵或者 cell 数组里面装的是文本卷头内容。有些老文件是 EBCDIC 编码SEGYMAT 在读取时会自动帮你转换成可读的 ASCII 字符串省去了自己处理 EBCDIC 的麻烦。3.3 不引包也能读卷头的备选方案不管用不用工具包下面这段脚本都建议存一份。它能在几秒钟内确认文件的基本参数排查问题的时候尤其管用fid fopen(shot1.sgy, r, ieee-be); fseek(fid, 3217-1, bof); dt fread(fid, 1, uint16); fseek(fid, 3221-1, bof); ns fread(fid, 1, uint16); fseek(fid, 3223-1, bof); fmt fread(fid, 1, uint16); fprintf(dt%d us, ns%d, fmt%d\n, dt, ns, fmt); fclose(fid);fseek的偏移量从 0 开始所以标准里写的第 3217 字节实际偏移是 3216。读取时指定了ieee-be这是大端字节序的关键稍后会详细讲。有了这组参数你再去决定后面是继续引包读全文件还是自己写逐道读取的脚本。4. 从道头到地震剖面别急着用imagesc4.1 先构好时间轴拿到 Data 矩阵后如果卷头里的采样间隔 dt 单位是微秒时间轴可以这样构造t (0:size(Data,1)-1) * dt * 1e-6;很多初学者直接imagesc(Data)结果横轴是道数、纵轴是样点序号看着像地震剖面但拿不出手。正确的做法是纵轴用时间或深度横轴最好换成道号或 CDP 号figure; imagesc(1:size(Data,2), t, Data); axis ij; xlabel(Trace No.); ylabel(Time (s)); colormap(gray); colorbar;这一步看起来简单但对后面所有分析都影响很大。时间轴搞错了频谱、速度分析全都会跟着错。4.2 数据显示前的归一化与裁剪地震数据动态范围非常大直接imagesc通常只能看到几条强反射弱信号全被淹没了。常规做法是设定一个裁剪阈值 clip把超出正负 clip 的值截断。clip 可以取最大绝对值的一定比例也可以用百分位数来设定clip prctile(abs(Data(:)), 90); DataClip max(min(Data, clip), -clip); imagesc(1:size(Data,2), t, DataClip);这只是显示层面的处理拿 clip 之后的数据去做频谱分析和反演是错误思路会直接改变振幅关系。另外还有一种是按道做归一化把每道除以该道最大值适合看弱信号轮廓但道间相对幅值关系会被破坏做 AVO 分析时一定要慎用。4.3 wiggle图与灰度图的简单做法imagesc的灰度图适合快速浏览整个剖面而 wiggle 图变波形图是地震解释里更常用的显示方式。MATLAB 没有现成的 wiggle 函数自己实现也不复杂每一道画一条曲线在横向上按道间距做平移正负极填充颜色区分。几百道的范围内用循环逐道plot性能还可以几千道就会明显变慢。如果只是做质量检查灰度图已经足够如果要做正式汇报或者解释再用 wiggle。我的经验是先灰度图快速看全局再抽关键道段画 wiggle 图瞧细节两个配合效率最高。5. 最常见的两处翻车字节序和IBM浮点5.1 字节序为什么读出来全是天文数字SEGY 标准默认是大端字节序也就是高位字节存在前面。而 PC 上的 MATLAB 默认按小端解释数据如果直接用fopen(fid, r)再fread一个int32读出来的值相当于把每个字节倒着念了一遍原本正常的数值会变成几亿几十亿的乱码。解决办法是打开文件时显式指定大端字节序fid fopen(shot1.sgy, r, ieee-be);使用工具包时一般不用自己操心这个问题因为工具包内部已经处理好了。但如果你在某个国产处理系统导出的文件上发现读出来数值异常先不要怀疑工具包自己用上面的方法读一下卷头和前几道确认字节序到底是不是标准大端。有的系统会自作聪明地输出小端 SEGY这种情况只能手动绕开工具包按小端去读。5.2 IBM浮点老数据的隐藏老坑数据格式码为 1 的 IBM 浮点是最容易坑人的地方。IBM 浮点和 IEEE 浮点虽然都占 4 字节但指数基数是 16 而不是 2指数偏置是 64 而不是 127直接用ieee-be读 IBM 浮点会得到完全错误的数值。老工区、老磁带导出的 SEGY 里这种格式很常见。转换思路是先按 4 字节无符号整数读出符号位看第 1 字节最高位指数看第 1 字节低 7 位尾数是后 3 字节组成的二进制分数。转换代码如下function v ibm2float(ibmBytes) % ibmBytes: N x 4 uint8按大端排列的IBM浮点 b double(ibmBytes); sign -2 * double(bitget(ibmBytes(:,1), 8)) 1; expo bitshift(bitand(b(:,1), 127), 0); frac b(:,2)/256 b(:,3)/65536 b(:,4)/16777216; v sign .* frac .* 16 .^ (expo - 64); end用的时候先读uint8再调用这个函数不要直接用fread读float32。SEGYMAT 内部有的版本能自动处理格式码为 1 的情况有的版本处理得不对所以拿到老数据时最好随机抽几道做数值合理性验证。5.3 拿到文件先做的五个检查我在处理一个陌生 SEGY 文件时会按这个顺序做快速体检先读二进制卷头确认 dt、ns、fmt再打印前两道道头检查 CDP 号、偏移距是否合理然后读第一道数据做统计看最大值、最小值、均值的量级如果是老数据随机抽一道做 IBM 浮点转换对比最后再看整个文件大小与ntr (filesize - 3600) / (240 ns*4)估算的道数是否吻合。这套流程走下来绝大多数异常情况都能提前暴露避免后面盲目处理造成大量返工。6. 大型SEGY文件与后续扩展6.1 问题场景单文件几十GB怎么办三维地震勘探的单文件动辄几十 GBReadSEGY一把梭把全部数据读进内存基本会直接卡死或者报内存不足。这时候要么用 SegMAT 的分道读取功能要么自己写循环逐道处理。我的习惯是先用一个小脚本把道头读完确认总道数和每道采样点数再按道区间分块读取把数据块交给下游处理函数用完一块丢一块。总道数的估算可以用文件大小来计算d dir(big.sgy); ns 2000; % 示例每道2000个采样点 ntrTotal floor((d.bytes - 3600) / (240 ns*4));如果文件里有多个扩展文本头这个估算会偏大但对道数级确认够用。6.2 逐道读取示例下面这段代码演示了如果因为某些原因必须自己逐道读取应该怎么写。这里假设数据格式码为 5也就是 IEEE 浮点fid fopen(big.sgy, r, ieee-be); blockLen 240 ns * 4; for itr 1:ntrTotal fseek(fid, 3600 (itr-1)*blockLen, bof); fread(fid, 240, uint8); % 跳过道头 trace fread(fid, ns, float32); % 按需处理当前道 end fclose(fid);如果数据格式码是 1把fread改成读uint8之后转成 N×4 矩阵再调用上面写的ibm2float函数。逐道循环在 MATLAB 里确实不算快但对于超大数据来说这是最稳的办法而且可以随时用fseek跳到任意一道灵活性远高于一次性全读。6.3 后续扩展方向读进 MATLAB 之后能做的事情就很多了叠前道集抽取、振幅补偿、频谱分析、FK 滤波、道均衡、时深转换都可以在此基础上按需实现。如果你的流程是深度学习方向的还可以把 SEGY 导成.npy或者 HDF5 格式再喂给 Python 框架避免两边来回折腾格式。工具包只是入口后面能走多远取决于你对自己的数据有多了解。最后说个自己的体会读 SEGY 这件事工具包一抓一大把但真正让项目不翻车的往往不是工具本身而是肯不肯花半小时搞清楚文件头和字节序。拿到一个读起来很怪的文件先别急着骂工具包回头检查格式码是什么、文件是不是被采集系统改过道头。搞清楚这些再复杂的 SEGY 文件也就是一道循环的事。本文还有配套的精品资源点击获取