简介这套 Matlab SEGY 读取工具包面向物探与地震数据处理人员基于 SegyMAT 实现 SEGY 格式地震数据的读取、解析、写入与转换。包内共有 75 个文件以 63 个 m 脚本和 7 个 fig 界面文件为主体同时提供 1 个 segy 示例数据、2 个 dat 数据及许可证与版本说明文件压缩后仅 516KB小巧完整。核心函数 ReadSegy、WriteSegy 支持批量导入导出配套多个 GUI 交互界面便于无命令行基础者操作同时包括 Su2Segy、Segy2Su 等格式转换工具以及 wiggle、绘图模块可直接用于地震剖面显示与简单分析。目前已有 646 人学习下载适合需要快速读取标准地震道数据、开展格式转换或二次开发的 Matlab 用户。借助该工具包可省去从零编写底层解析逻辑的繁琐过程直接调用现成接口完成道头信息提取、数据矩阵获取和基础可视化。1. 从SEGY到地震剖面这个MATLAB工具包到底解决什么问题拿到野外采集回来的SEGY文件直接用文本编辑器打开满屏都是乱码用商业解释软件看剖面又得装一套笨重的license。这个MATLAB工具包干的事很简单把一个SEGY文件完整拆成文本头、二进制头、道头和数据体再组装成MATLAB矩阵让你在脚本里自由处理地震数据。它适合两类人一类是做勘探数据处理的学生想快速把SEGY读进MATLAB做滤波、频谱分析另一类是搞算法验证的工程师不想被重型软件绑架只想命令行里一条命令把数据拽出来。SEGY这套格式在地下地球物理领域用了三十多年看似原始但看懂它的结构很多数据处理痛点就迎刃而解。2. SEGY格式与工具包设计先把二进制黑匣子拆开2.1 文本卷头与二进制卷头3200 400 字节的元信息SEGY文件开头是一块固定长度的元信息区。文件最前面的3200字节是文本卷头Textual Header按标准规定每行80字节、共40行存放文件描述、采集参数、处理记录等人工可读信息。早期的SEGY文件这段文本用EBCDIC编码后来不少文件直接用ASCII所以读取时得做两种编码的兼容。接着是400字节的二进制卷头Binary Header所有关键数值参数都藏在这里。它里面有几个字段做读取时绕不开采样间隔单位微秒、每道采样点数、数据格式码。采样点数决定了一个道有多长格式码则告诉你怎么解释后续的数据体字节。格式码是SEGY里最容易翻车的字段之一——1代表IBM浮点、2代表32位定点整数、3代表16位定点整数、4才是IEEE浮点。很多早期的采集数据用IBM浮点不转换直接用fread读出来就是一堆莫名其妙的数。2.2 道头与道数据每条地震道的独立信息块二进制卷头之后整个文件剩余部分被划分为一条条的道trace。每条道内部又分成两部分240字节的道头Trace Header和若干字节的采样数据。道头里记录的是一次激发对应的野外参数——道号、炮点坐标、检波点坐标、偏移距、采样点数等。这些信息在后期做道编辑、切除、动校正时非常有用。数据体的大小怎么算每道数据字节数 道头240字节 采样点数 × 单采样点字节数。单采样点字节数由格式码决定IBM浮点和32位整数都是4字节16位整数是2字节IEEE浮点可能是4字节或8字节。工具包读取时先解析二进制头拿到采样点数和格式码算出单道大小然后就能按这个固定步长在整个文件里循环跳读。注意标准SEGY还有一种扩展文本卷头Extended Textural Header插在二进制卷头之后、第一道之前占用的也是3200字节遇到这种文件总偏移量就要多加3200。2.3 工具包目录结构与模块划分这套工具包按功能拆成几个独立的m文件目录结构大致长这样segy_tools/ ├── read_segy.m # 主入口函数返回结构体 ├── parse_text_header.m # 文本卷头解析支持EBCDIC/ASCII ├── parse_binary_header.m # 二进制卷头关键字段提取 ├── parse_trace_header.m # 单道道头解析 ├── read_trace_data.m # 按格式码读取一组道数据 ├── ibm2ieee.m # IBM浮点转IEEE浮点 ├── plot_seismic.m # 地震剖面绘图 └── example_demo.m # 完整调用示例主入口read_segy的返回值是一个结构体S里面包含textHeader、binHeader、traceHeaders、data四个字段。data是二维矩阵行对应采样点、列对应道号这也是后续绝大部分处理脚本最舒服的输入格式。工具包把格式解析、数据读取、可视化分开你想替换某一段实现时不用动其它模块。接口上保持一个函数干一件事方便你在此基础上扩展。3. 核心模块拆解一个SEGY文件从字节流到数据矩阵的完整链路3.1 文本头解析EBCDIC与ASCII的自动识别文本卷头在标准SEGY里是EBCDIC编码但实际流传的数据很多是ASCII所以解析函数得先探测编码类型。常见的做法是统计字节分布EBCDIC文本里大量字节落在0x81到0xC9、0xE5到0xF9区间ASCII文本则集中在0x20到0x7E之间。function textStr parse_text_header(fid) % fid 为已打开的文件句柄当前游标位于文件起始处 % 读取3200字节文本卷头原始字节 rawBytes fread(fid, 3200, uint8); % 统计可打印ASCII字符比例超过一半则判定为ASCII编码 asciiRatio sum(rawBytes 32 rawBytes 126) / 3200; if asciiRatio 0.5 textStr native2unicode(rawBytes, US-ASCII); else textStr ebcdic2ascii(rawBytes); % 转义EBCDIC end % 按80字符一行切分共40行 textStr reshape(textStr, 80, 40); textStr cellstr(textStr); end这段代码先读原始字节通过统计可打印ASCII字符占比来决定走哪条解码路径。asciiRatio大于0.5的判断依据是EBCDIC编码中可读字符的码点分布与ASCII差异极大误判概率很低。native2unicode按US-ASCII解码后reshape成80列40行的矩阵每行正好是一条文本卷头记录。返回的cellstr数组方便你逐行查看文件来源和采集参数。3.2 二进制头参数提取与数据格式分发二进制卷头里真正影响后续读取的只有那么四五个字段但它们的字节偏移是固定的直接用fseek跳过文本卷头再按偏移读取最稳妥。function binInfo parse_binary_header(fid) % fid: 文件句柄已越过3200字节文本卷头区域 % 采样间隔偏移3217-3218单位微秒 fseek(fid, 3216, bof); % 3217-13216MATLAB从0开始计偏移 sampleInterval fread(fid, 1, int16); % 每道采样点数偏移3221-3222 fseek(fid, 3220, bof); samplesPerTrace fread(fid, 1, int16); % 数据格式码偏移3225-3226 fseek(fid, 3224, bof); formatCode fread(fid, 1, int16); binInfo struct(sampleInterval, sampleInterval, ... samplesPerTrace, samplesPerTrace, ... formatCode, formatCode); end注意fseek的第三个参数bof表示相对文件开头计算偏移。偏移3216对应规范里的字节3217因为MATLAB计数从0开始而规范从1开始。采样间隔用的是微秒换算成秒要除以1e6。formatCode是后面读取数据体的分派依据常见的1、3、4三种值分别对应IBM浮点、16位整数、IEEE浮点。如果你的数据是别的格式码比如8字节IEEE浮点这个工具包需要自行扩展数据格式分支。3.3 道数据循环读取与IBM浮点转换读道数据是整个工具包里最重的一段逻辑下面这个函数完成单道数据的读取和格式转换function [traceData, traceHeaders] read_trace_data(fid, binInfo, traceCount) % 读取全部道数据返回矩阵和道头结构 samples binInfo.samplesPerTrace; formatCode binInfo.formatCode; % 根据格式码确定单采样点的字节数 switch formatCode case 1 % IBM浮点4字节 sampleBytes 4; readFunc ibm2ieee; case 2 % 32位定点整数 sampleBytes 4; readFunc (x) double(typecast(x, int32)); case 3 % 16位定点整数 sampleBytes 2; readFunc (x) double(typecast(x, int16)); case 4 % IEEE浮点4字节 sampleBytes 4; readFunc (x) double(typecast(x, single)); otherwise error(不支持的数据格式码: %d, formatCode); end traceData zeros(samples, traceCount); % 预分配矩阵 traceHeaders cell(traceCount, 1); for k 1:traceCount traceHeaderBytes fread(fid, 240, uint8); traceHeaders{k} parse_trace_header(traceHeaderBytes); % 读取采样点数个原始字节 rawData fread(fid, samples * sampleBytes, uint8); traceData(:, k) readFunc(rawData); end end这段代码把数据读取拆成两步先读240字节道头再读定长的采样数据。预分配traceData矩阵是关键一步如果边读边拼接矩阵数据量大时性能下降非常明显。readFunc按格式码分派到不同的转换函数IBM浮点走专门的ibm2ieee其它格式直接用typecast重解释字节。typecast不会修改底层数据只是换一种数据类型去看同一段内存所以比用fread直接指定精度更可控。3.4 IBM浮点转IEEE浮点的bit级运算IBM浮点格式是SEGY数据的老朋友结构上跟IEEE完全不一样1位符号、7位指数基数为16、24位尾数。转换的核心逻辑就是对这32个bit做重新解释。function value ibm2ieee(bytesVector) % bytesVector: 4字节的uint8数组大端序 % 将4字节拼成一个32位无符号整数 ibmInt uint32(0); for i 1:4 ibmInt bitor(bitshift(ibmInt, 8), uint32(bytesVector(i))); end % 提取符号位、指数、尾数 signVal bitget(ibmInt, 32); exponent bitshift(bitand(ibmInt, uint32(hex2dec(7F000000))), -24); fraction bitand(ibmInt, uint32(hex2dec(00FFFFFF))); % 零值特殊处理 if exponent 0 fraction 0 value 0; return; end % IBM浮点: (-1)^sign * (fraction/2^24) * 16^(exponent-64) mantissa double(fraction) / 2^24; expVal 16^(double(exponent) - 64); value (1 - 2 * double(signVal)) * mantissa * expVal; endIBM浮点的指数是以16为底的这是最容易踩坑的点——如果你用2为底去还原指数转换出来的振幅值会差出好几个数量级。代码里bitor和bitshift配合把4个字节按大端序拼成uint32这是因为SEGY文件规定大端字节序。如果你读到的数据在PC上表现异常优先怀疑字节序问题而不是IBM转换本身。4. 跑通一个真实SEGY文件从命令行到地震剖面4.1 环境准备与快速调用工具包对MATLAB版本要求不高R2016b以上就行核心逻辑只用到了fopen、fread、fseek、typecast这些基础函数不需要额外的工具箱。拿到工具包目录后先把segy_tools文件夹加进MATLAB搜索路径addpath(你的路径/segy_tools); savepath; % 保存路径设置下次启动不用重复添加然后找一份SEGY数据文件调用read_segy主函数。注意文件路径尽量用绝对路径相对路径在脚本跨目录执行时容易找不到文件。调用方式如下S read_segy(C:/data/line_a001.segy);返回的S是一个结构体S.data是二维数组行数等于每道采样点数列数等于总道数。S.binInfo里有采样间隔和格式码S.traceHeaders里每一行存储了对应道的野外采集参数。想快速检查文件有没有读对直接看数据矩阵的尺寸和值域size(S.data) % 期待是 [采样点数, 道数] min(S.data(:)) % 看看振幅范围是否合理 max(S.data(:))如果发现数据范围大得离谱先看格式码对不对如果范围小到几乎全是0看看是不是IBM浮点转换那条路径出了问题。4.2 数据清洗与道编辑读取只是第一步实际处理前通常要做数据体检。SEGY文件里偶尔混进坏道——采样点数不对、振幅全部为0、或者某个道的数据是重复采集。下面的代码利用道头信息做一道简单的筛选% 提取每道的道号供后续参考 traceNums cellfun((h) h.traceNumber, S.traceHeaders); % 找出振幅全为0的死道 deadTraces find(max(abs(S.data)) 0); % 去掉死道后重建数据矩阵 if ~isempty(deadTraces) S.data(:, deadTraces) []; S.traceHeaders(deadTraces) []; fprintf(移除了 %d 条死道\n, length(deadTraces)); end这段逻辑里max(abs(S.data))按列取最大值所有值都为零的列就是死道。cellfun从traceHeaders结构体数组里提取traceNumber字段方便你后续按道号排序或者查道。直接从矩阵里删列是最粗暴但有效的坏道处理手段处理完记得同步删道头否则道号和数据对不上。4.3 地震剖面可视化工具箱工具包里带了一个简易绘图函数plot_seismic用法很简单传入数据矩阵和采样间隔就能出图plot_seismic(S.data, S.binInfo.sampleInterval);内部实现用的是imagesc加自定义色彩映射不是海底底图那种复杂的渲染function plot_seismic(data, sampleIntervalUs) [nSamples, nTraces] size(data); % 构建时间轴单位毫秒 timeAxis (0:nSamples-1) * sampleIntervalUs / 1000; % 用imagesc画图并设置坐标轴 imagesc(1:nTraces, timeAxis, data); set(gca, YDir, normal); % 时间轴从上往下增大 colormap(flipud(gray)); % 地震剖面惯用黑-白-黑颜色 xlabel(道号); ylabel(时间 (ms)); end这里用imagesc而不是直接用plot的wiggle方式好处是画超多道数据的剖面时性能好很多坏处是密集剖面看起来没有波形图那么直观。你真要做出版级图形建议在这基础上改用传统的variable-density绘图。set(gca, YDir, normal)很关键默认imagesc会翻转Y轴不设置的话时间方向会被倒过来这是新手的常见翻车点。5. SEGY踩坑指南字节序、编码与扩展道头5.1 在别人机器上正常的文件我读出来全是乱数现象同一个SEGY文件同事在Linux工作站上读出来剖面正常你在Windows笔记本上读出来数据值巨大或者全是噪声。原因SEGY标准规定大端字节序而x86平台默认小端字节序。很多商业软件读取时会自动做字节序转换但自己写的fread如果直接指定int16或者singleMATLAB会按本机小端序去解释高低字节颠倒。这时候采样间隔可能变成几千、采样点数变成负数数据体就更不用说了。解决读取所有多字节字段时在fread里显式指定机器格式或者读完用typecast配合swapbytes。例如读采样间隔时改成fseek(fid, 3216, bof); rawBytes fread(fid, 2, uint8); sampleInterval double(typecast(swapbytes(uint16(typecast(rawBytes, uint16))), int16));更省事的做法是fread直接加l结尾的精度字符串比如int16改成int16加b参数big-endianfread(fid, 1, int16, 0, b)这是MATLAB里绕过字节序问题的正道。5.2 文本卷头全是乱码一行可读内容都没有现象文本卷头打印出来是空行、特殊符号或者看不太懂的字符序列。原因早期采集设备生成的SEGY文件用EBCDIC编码而非ASCII。EBCDIC里字母和数字的码位和ASCII完全不同按ASCII去解码自然是一团乱码。解决用解析函数里的ebcdic2ascii做转换或者先检测编码再决定路径。我在工具包里直接实现了简化版EBCDIC映射表覆盖空格、数字、大写字母、小写字母这些常用字符。需要说明的是EBCDIC也有多种方言变体如果遇到转换后个别字符还是不对多半是变体差异补映射表就行。5.3 读出来的采样点数跟道头里写的不一致现象数据矩阵尺寸跟预期的采样点数对不上或者读完一道后文件游标位置错乱后面所有道的道头都解析错误。原因二进制卷头里的每道采样点数字段是全局设定但具体某一道的道头里也有采样点数字段。有些处理系统会修改局部道头里的采样点数导致文件内不同道的实际字节长度不一致不能按固定步长跳读。解决读取每道时不要信任二进制卷头的全局采样点数而是以当前道的道头字段为准做动态计算。工具包在read_trace_data里先解析道头、取出该道实际的采样点数再按这个数读取数据体。代价是每道都要多读一次240字节道头但换来的是对非标准文件的鲁棒性值。5.4 文件里有扩展文本卷头读到的数据全是偏移的现象读完文本头和二进制头之后fread读道头时第一个道的240字节里全是乱码或者首次读就报错文件过早结束。原因SEG-Y rev1以后允许在二进制卷头后面插入一个额外的3200字节扩展文本卷头而不少采集系统默认就写这个扩展区。你按标准的3600字节起始偏移去找第一道自然就偏了3200字节。解决读完二进制卷头后再往后检查有没有扩展文本卷头。最直接的办法是读取二进制头第3517-3518字节的扩展文本卷头计数或直接看文件偏移3600字节处的头部是否像道头。稳妥的做法是工具包里加一个选项读二进制头后先尝试解析偏移3600处的240字节看道号字段是否合理不合理就把起始游标往后挪3200字节再解析一次。这类玄学问题说白了就是SEGY标准衍生出多个版本处理时多留条后路总没错。6. 二次开发把工具包改造成自己的地震数据处理管线的三个方向工具包到手之后最常见的改造需求是把读取结果喂给后续的数据处理脚本。我建议重点做三件事批量读取、格式转换、并行加速。批量读取的场景很典型——一个工区几百个SEGY文件每个文件对应一条测线想一次性全部载入。改造思路不复杂写一个外层脚本遍历文件列表逐个调用read_segy把每个文件的S.data存到cell数组里。fileList dir(C:/data/raw/*.segy); allData cell(length(fileList), 1); for i 1:length(fileList) S read_segy(fullfile(fileList(i).folder, fileList(i).name)); allData{i} S; end如果文件太多导致内存吃紧可以只读道头信息不读数据体或者按道索引分段读取。工具包目前的read_segy是整文件读入做内存映射改造时你可以把fread改成按需跳读思路跟read_trace_data一致。格式转换是另一个常用改造点。SEGY是勘探行业通用格式但深度学习、信号处理社区更常用npy或h5。在read_segy返回的S.data基础上直接调用MATLAB的save或者写一个转h5的脚本就能对接下游Sw S.data; writematrix(Sw, seismic_data.csv); % 道行为主方便导入Python并行加速方面read_trace_data里逐道读取的循环是天然的并行候选。用parfor替换for注意每个迭代里fread是独立文件游标操作MATLAB的fread并不会自动互斥所以需要先预分配游标或者改用每道单独fopen。我踩过的坑是直接parfor跑同一个fid会随机出错最后拆成每道开一次文件句柄才稳定。从这里往后我每次改完并行逻辑都会强制跑一遍全工区数据比对结果确认振幅和道头信息没被并行搞乱这个习惯救了我好几次。希望帮到你。本文还有配套的精品资源点击获取