SG2地震数据转DAT：segdat.m实现二进制到文本的格式转换

简介面向地震数据处理与分析人员这份资源包聚焦地震记录由二进制采集格式向文本数据格式转换的需求提供基于MATLAB的转换脚本。压缩包内仅包含1个m格式的MATLAB脚本文件压缩包整体大小仅约1KB结构单一轻量便于直接查看和调试。脚本围绕格式转换的核心流程设计依次实现SG2二进制数据解析、采样率与通道信息提取、按列整理地震参数、输出可读DAT文本文件等关键步骤同时保留了必要的参数调整位置可适配不同来源的地震记录。目前已有607人学习/下载适合地震监测、工程地震及相关科研人员使用能够快速打通两类数据之间的通道减少重复编写底层解析代码的时间帮助使用者更清晰地理解不同存储方式对后续分析的影响。1. 从sg2到datsegdat.m让地震数据转换不再卡壳流动台站收回的原始记录是SG2二进制道文件而震相分析脚本只认DAT文本这种情况在临时观测和应急数据处理中几乎每周都会碰到。SG2保留了最完整的地震波形适合采集端和存储端使用但一旦进入跨软件流程交接解析不了二进制格式后面的分析工作就只能停下来。segdat.zip里的segdat.m就是为这个场景准备的它读出SG2道头里的采样率、通道号、时间戳和采样点数再把每个采样点按行写成DAT文本同时附带可读的元数据注释。需要调整设备差异时只改函数入口的几个字段即可不必为每个台站重写脚本。它适合台站运维人员、数据分析工程师以及刚进入地震数据处理方向的学生。2. SG2与DAT格式的地震数据存储差异2.1 SG2的二进制封装道头、数据区与采样点排列SG2本质上是按道组织的二进制记录。一个文件可以包含一道或多道波形每道前面是固定长度的道头后面紧接数据区。道头的长度在不同采集系统间并不统一常见的规格有512字节、1024字节数据区才是真正按采样率排列的波形序列。道头里记录事件时间、采样点数、采样间隔、通道序号以及增益等参数。因为字段是排列在连续字节中的解析时不能靠名字只能靠偏移和变量类型。fread读回来的字节数组不会告诉你哪几个字节是采样点数必须根据格式文档或者实际观测去验证。即便同一个厂家的设备固件版本不同也可能调整字段位置所以在segdat.m里我习惯把关键偏移放在文件开头集中管理。数据区通常按定点数或浮点数排列。定点数常见的是int16或int32浮点数则可能是single或double。选择哪种类型取决于采集设备的ADC位数。int16的优点是省空间动态范围大约为96dBint32能给出更大的动态范围但文件体积翻倍。为了后续分析口径统一segdat.m会先按原始类型读入再统一转成double振幅。2.2 DAT文本格式的行组织与常用字段布局地震处理中的DAT文件没有全球统一标准但多数情况下是纯文本每行一个采样点也可以用多列表示时间、幅值和多个分量。读取端只需要用readmatrix或load就能加载这让DAT成为数据交换时的首选文本载体。segdat.m输出的DAT采用注释头加数据区的结构。文件开头以井号开头记录事件时间和采样率数据区每行一个双精度浮点数。如果要输出多个通道最好为每个通道单独生成一个dat文件而不是把几百个通道塞进同一个文本。单列布局的好处是后续做互相关或谱分析时按文件读取不需要处理复杂的分隔符。DAT格式的主要代价在体积和性能。实测中int32的SG2文件写成文本后会膨胀4到6倍读取时间也从毫秒级变成秒级。它更适合用于分析结果的归档和平台间交换不适合做长期原始数据仓库。2.3 转换的四项对齐关系与验证方法把SG2转成DAT前先要明确四项对齐信息采样点数对齐、采样间隔对齐、通道号对齐、时间戳对齐。下表按segdat.m默认实现列了一组最常用的对应关系。地震参数SG2中的典型位置DAT输出位置转换注意点采样点数道头内固定字偏移数据区行数读取时按字节类型转换采样间隔道头内固定字偏移头部注释# sample_rate注意单位是毫秒还是微秒通道号道头内固定字偏移头部注释# channel多道数据建议拆分为多个文件事件时间道头ASCII区头部注释# event_time保留原始时区信息数据体道头后的数据区每行一个采样值int16/int32转浮点时要标定拿到一个陌生SG2文件时可以先用MATLAB做一次快速检查fid fopen(event_20240708.sg2, rb, ieee-le); hdr fread(fid, 64, *uint16); disp(hdr(1:20)); fclose(fid);这里先用小端字节序读取前64个字观察前20个数值是否符合设备已知参数。如果数值全都跳动得没有规律或者采样点数是一个天文数字多半是字节序选错应将ieee-le换成ieee-be再试。fread中的*uint16表示按两个字节一个字读取道头里很多字段是16位或32位先按16位观察结构比较容易定位问题。3. segdat.m 的实现与转换流程3.1 入口参数设计把字节序、道头长度和标定交给调用者写转换脚本最容易犯的错误是把所有参数写死在代码里。道头长度、字节序、采样间隔单位这三个值在不同采集器之间几乎没有规律可循。segdat.m把这三个值作为可选项暴露在函数签名里调用者可以根据实测数据传参不用每次编辑源码。function datfile segdat(sg2file, datfile, opts) arguments sg2file (1,1) string datfile (1,1) string output.dat opts.byteOrder (1,1) string ieee-le opts.headBytes (1,1) double 1024 opts.scale (1,1) double 1.0 end fid fopen(sg2file, rb, opts.byteOrder); if fid 0 error(segdat:OpenFailed, 无法打开文件: %s, sg2file); end % 后续道头解析、数据读取与写出逻辑 fclose(fid); end代码中的arguments是MATLAB R2019b之后支持的输入校验语法。调用segdat(abc.sg2)会走默认参数调用segdat(abc.sg2,out.dat,byteOrder,ieee-be)可以快速切换大小端。scale参数用于传感器标定读到的整数计数乘以它才能变成物理单位。为什么headBytes要独立设置因为数据区起点不一定总是等于道头长度。有的设备在道头后面还会附加一段量纲说明这时需要把headBytes设置成道头长度加额外偏移。提供一个可调参数比把取数据区位置的逻辑写死在代码里要干净得多。3.2 道头字段读取typecast、fseek 与物理量换算道头解析依赖fseek和typecast。先定位道头开始位置读取固定字节再用typecast把字节组合成数值。MATLAB里不能直接把uint8数组当uint32用必须显式做类型转换。hdr fread(fid, 2048, *uint8); % 以下偏移按常见SG2结构示例真实文件需对照格式手册核对 npts double(typecast(hdr(17:20), uint32)); dtMicro double(typecast(hdr(21:24), uint32)); channelId double(typecast(hdr(25:26), uint16)); % 定位到数据区并读取全部采样点 fseek(fid, opts.headBytes, bof); raw fread(fid, npts, *int32); phys double(raw) * opts.scale;这段代码假设采样点数记录在第17到20字节采样间隔记录在第21到24字节通道号记录在第25到26字节。这是很多SG2类文件通用的一套布局但并非硬标准。实际使用时要拿真实文件验证偏移值。typecast(hdr(17:20),uint32)会把第17字节作为低位第20字节作为高位在小端模式下正好得到正确数值如果打开时指定了大端字节顺序就要反过来理解。phys的计算包含了scale。仪器响应标定值一般是一个极小的浮点数例如1.6e-9 m/s/count乘以整型计数后得到的是m/s单位的速度序列。如果只关心相对波形形态把scale设为1即可。数据区如果是single浮点而不是int32应把fread的数据类型改成*single随后再做一次double(raw)。判断依据是道头里的数据类型标识和每条采样点的字节宽度。3.3 DAT写出与分块写盘机制写出步骤建议分块。用一个较大的缓冲区把数据格式化后一次性fprintf避免在循环里逐点调用造成写盘过慢。下面是最常用的单通道输出写法。out fopen(datfile, w); fprintf(out, # event_time: 2024-07-08T12:34:56\n); fprintf(out, # sample_rate: %.6f Hz\n, 1e6 / dtMicro); fprintf(out, # channel: %d\n, channelId); fprintf(out, # data_start\n); block 524288; for s 1:block:npts e min(s block - 1, npts); fprintf(out, %.8e\n, phys(s:e)); end fclose(out);前三行fprintf负责写注释头后续读取端可以用readmatrix(out.dat,CommentStyle,#)直接跳过井号行。采样率从dtMicro微秒换算成Hz公式是1 / (dtMicro * 1e-6)。如果dtMicro实际单位是毫秒而不是微秒换算因子需要改成1e-3。分块写的作用是让长时间连续记录转成DAT时内存不会暴涨。每个block块取52万行左右即可输出数据按行排列每行一个科学计数法浮点数。%.8e保留8位有效数字对地震观测数据来说已经足够同时能让文件保持紧凑。3.4 输出参数与文件体积预期输出DAT的体积是这类文本转换绕不开的问题。下面按每秒采样1000点、float64文本格式估算每行约15字节。一天连续记录会产生约130MB文本虽然远小于同等时长的原始SG2但在批量处理时仍要考虑磁盘空间。输出参数建议值说明precision%.8e8位有效数字浮点输出blockSize524288每批最多写52万行separateChannelstrue多道时按通道分别写文件includeHeadertrue保留井号注释行转换完成后用文本编辑器打开文件头部核对注释行和数据行数是否匹配。如果行数比采样点数少很多优先检查3.2里的headBytes是否偏大把数据区起点定到了波形的中段。4. 参数配置、边界处理与常见排错4.1 字节序、采样间隔与道头偏移的排查用segdat.m处理一批外部数据时先跑一个样本文件看输出的头部注释是否符合常识。采样率若在0.0001附近大概率是采样间隔单位被误当成微秒而原始数据里的单位其实是秒采样率若变成几千万Hz则多半是字节序颠倒。症状可能原因排查动作采样点数是一个巨大整数字节序或偏移错误打开十六进制查看前4字节采样率接近0单位换算错确认dtMicro代表ms还是us输出波形只有一条平线数据区偏移错误从headBytes附近逐段试探输出文件体积少一半数据区是int16但按int32读检查道头数据类型标识排查时不要直接相信读出的头字段。把道头前256字节导出成十六进制文本手动计算前几个字段的十进制值和设备面板上显示的采样率做对比。若差2倍或4倍基本是字节序或字宽选错若差一个固定倍数则多半是采样间隔单位换算错。4.2 坏道、丢帧与饱和数据的统计方式地震记录中因为采集器死机、雷击、传感器升压不到位会出现坏道。转换时统计非有限值和满幅值比例并写进注释后续处理就能知道这段数据的真实质量。% 统计NaN和Inf数量 bad ~isfinite(phys); % 统计接近满幅的饱和样点 sat abs(phys) abs(opts.scale) * 0.999; fprintf(out, # bad_count: %d\n, nnz(bad)); fprintf(out, # sat_count: %d\n, nnz(sat));bad数组标记NaN和Inf点位。遇到大量NaN时不建议直接删除因为不知道这些位置原始采样值是否还有分析价值。segdat.m只统计并原样保留后续拆帧或滤波时再决定补零还是做数据修补。饱和点同理保留原值并标记比插值更可靠因为插值会引入虚假低频成分。4.3 地震dat与其他同名格式的辨别在网上一搜dat格式会出现微信dat文件查看器、佳能dat文件恢复这样的内容。微信dat是聊天图片缓存佳能dat文件恢复针对存储卡残留文件它们与地震时间序列毫无关系。判断dat属于哪一类不能看扩展名要看文件内部结构。segdat.m生成的DAT应当是首行带井号注释、接下来每行一个浮点数的纯文本而微信dat加密后是乱码佳能dat文件恢复依赖文件系统上下文。拿地震脚本去读这些文件会在fopen之后的fread阶段直接报错。理解这个区别对运维场景有价值存放数据时采用规范命名把地震dat单独放在events目录并让批量脚本自动生成日志。否则半年后看到一批.dat文件很难定位哪一个对应哪条测线。4.4 转换结果与原始波形的快速对比转换结束后用MATLAB把DAT读回截取一段1000个采样点与SG2原始波形叠加对比。如果两条曲线重合说明偏移和标定都正确。d readmatrix(output.dat, CommentStyle, #); % 从原始SG2中读取同一位置的数据段作对比 fid fopen(event_20240708.sg2, rb, ieee-le); fseek(fid, opts.headBytes, bof); ref fread(fid, 1000, *int32); fclose(fid); plot(1:1000, double(ref) * opts.scale, x); hold on; plot(1:1000, d(1:1000, 2), -);这里使用readmatrix读取DAT时默认会把第一列当作行号因此数据值落在第二列d(:,2)中。若输出是单列数据则需要用load(output.dat)或者readmatrix(output.dat,NumHeaderLines,4)来处理注释头。5. 批量处理与结果验证的进阶技巧5.1 一键批量转换台站目录台站数据通常按天或按小时切成很多小文件一次要跑几十个。用dir列出所有SG2文件再循环调用segdat即可。files dir(raw/*.sg2); for k 1:numel(files) in fullfile(files(k).folder, files(k).name); out fullfile(converted, strrep(files(k).name, .sg2, .dat)); segdat(in, out, headBytes, 1024); end这里strrep只替换文件名后缀不会误改目录部分。如果某份数据的道头偏移不同循环里要加一条日志记录失败文件再单独处理。5.2 把segdat扩展成dat转sg2回写接口很多分析软件只接受DAT但数据交换时对方又要求sg2。可以把segdat.m逆过来写先从DAT注释头中读取采样率等元数据再把采样值按整型编码写回。关键在于scale因子写回前把振幅乘以合适系数转成int32确保动态范围不浪费。d readmatrix(output.dat, CommentStyle, #); fid fopen(restored.sg2, wb); fwrite(fid, int32(d(:, 2) * opts.scale), int32); fclose(fid);这里没有重建道头只适合实验验证。正式回写时还要把事件时间和通道号写进道头不能用占位符糊弄。5.3 把segdat.m接进数据质量日志实际项目中我会把segdat.m封装为批处理脚本让它顺手输出每个文件的质量摘要包括采样点数、坏道占比、饱和占比和输出文件大小。这样从采集端拿到原始SG2经过一次转换就能同时得到DAT文件和可供后续分析的质量清单。质量日志采用制表符分隔每行一个文件字段顺序固定file_name _ npts _ sample_rate _ bad_ratio _ sat_ratio _ dat_kbfprintf(logFid, %s\t%d\t%.1f\t%.4f\t%.4f\t%.1f\n, ... name, npts, sampleRateHz, badRatio, satRatio, fileSizeKb);后续做频谱分析时我会先加载这份质量日志把坏道比例大于5%的文件过滤掉再进入计算流程。本文还有配套的精品资源点击获取