简介针对SAC格式地震数据的读取需求这份压缩包提供了一段MATLAB脚本rdsac.m面向地震学、地球物理研究领域的科研人员、学生及波形数据分析工程师。SAC是地震学中广泛采用的标准数据格式可记录时间序列、频谱和事件元数据该脚本的核心任务是从二进制SAC文件中提取波形数据与头部信息将其转换为MATLAB可直接运算的数组同时解析采样率、起始时间、台站位置等参数为后续分析提供规范的数据结构。压缩包仅包含1个m脚本大小约1KB轻量简洁适合在MATLAB环境下直接调用借助它使用者无需从零编写底层读取代码即可快速完成数据导入并进一步利用MATLAB丰富的信号处理工具箱开展滤波去噪、频谱分析、震相识别和波形可视化等工作批量处理大量SAC文件时也能显著提升效率无论是科研分析还是日常教学演示都非常方便。目前已有247人学习使用适合作为SAC数据处理与MATLAB交互的入门工具也便于根据实际需求自行修改和扩展。1. 地震波形数据与SACrdsac究竟把什么交到你手上地震台站记录到的连续波形数据绝大多数以SAC格式保存。SAC是Seismic Analysis Code定义的一种二进制格式头段固定占512字节存放元数据后面紧贴单精度浮点数据。直接拿Matlab读这个文件并不复杂但头段里几十个字段的偏移、字节序、字符段填充规则一旦记错读出来的波形时间轴漂移、幅度全乱。rdsac.m做的事情就是把SAC头段解析成一个Matlab结构体并取出数据让读文件变成一行调用。这里从rdsac.m.tar.gz这个包里最常见的rdsac函数说起按“格式背景 → 安装配置 → 单文件读取 → 批量处理 → 进阶写回”的顺序把用Matlab处理地震SAC数据时最实际、最容易出问题的路径走一遍。适合用Matlab做台网数据后处理、需要快速预览波形、或者要把SAC数据接进自己信号处理流程的工程师。2. 认识SAC头段rdsac读取前的格式背景与配置检查2.1 SAC文件的512字节头段里存了什么SAC文件的物理结构其实非常简单前面是固定长度的头段后面是波形数据。头段512字节由三块组成浮点区、整数区、字符区分别存放连续型参数、整型参数和台站名、分量名这类字符串。采样间隔delta、起始时间b、采样点数npts、台站名kstnm、分量名kcmpnm这些字段在做地震数据处理时几乎每次都会用到。rdsac这类Matlab脚本的本质就是打开文件后按SAC头段布局逐字段fread出来再根据npts字段读取紧随其后的数据体。也就是说它不依赖任何外部库或mex文件纯Matlab实现这也是为什么它能在R2019到R2026b这么宽的版本范围内稳定运行不要求额外工具箱。SAC头段里最常用的字段和rdsac返回结构体字段的对应关系如下SAC字段含义rdsac返回字段常见写法用途采样间隔秒DELTA / delta时间轴换算的步长起始时间秒B / b相对于参考时刻的起点偏移采样点数NPTS / npts数据长度读取数据体的依据台站名KSTNM / kstnm事件归组、台站筛选分量名KCMPNM / kcmpnm区分BHZ/BHN/BHE三分量台网名KNETWK / knetwk按台网批量筛选参考年/日/时/分/秒NZYEAR, NZJDAY等还原波形绝对UTC时间震源经纬度EVLA, EVLO与台站坐标算震中距注意rdsac不同流传版本返回的字段名大小写不完全一致后面第4章会专门说这个坑。拿到任何一份rdsac.m先不要猜字段名动手打印一次最稳妥。2.2 解压rdsac.m.tar.gz并让Matlab找到函数标题里那个rdsac.m.tar.gz本质就是把rdsac.m及若干辅助脚本打包后的分发文件。在Linux或macOS下解压非常直接tar -xzf rdsac.m.tar.gz ls -l rdsac*Windows下没有直接支持tar.gz右键解压的用7-Zip或WinRAR解压两次先解出.tar再解出文件夹。解压后不要急着双击打开rdsac.m先把文件目录告诉Matlab。% 把解压后的目录加入MATLAB搜索路径 addpath(/home/seismo/tools/rdsac); % 保存路径设置保证下次启动MATLAB仍然生效 savepath;addpath把rdsac所在目录临时加入工作路径savepath把当前路径设置写入pathdef保存避免每次启动重新添加。如果解压目录会移动不建议写死在这个位置实际工程里我一般会在自己的启动脚本startup.m里用mfilename(fullpath)动态拼出上级目录再addpath那个目录这样整个工具库随身携带不会断。2.3 用which和help确认rdsac已就绪路径配置完成之后先验证rdsac确实可见再跑正式代码。which rdsac help rdsac % 若which返回空说明路径没生效 % 若help提示找不到函数检查文件名是否确实为rdsac.mwhich返回的是一串完整路径看到路径就代表Matlab能定位到这个函数。help则输出rdsac的函数签名和调用格式不同版本的rdsac这时的输出会略有差异但一般都会写明[data, header] rdsac(filename)或类似形式。这一步花不了十秒钟却能省下后面“明明写了addpath却还是Undefined function”的排查时间。3. 用rdsac把SAC波形读进Matlab并画出地震记录3.1 单文件读取的最小命令与头段字段核对假设目录下有一个2010年某事件的垂直向波形文件文件名按IRIS常见命名方式带上了开始时间但最终准确的时间信息在头段里。读取它只需要一行[sig, hdr] rdsac(2010.001.12.00.00.0000.BHZ.SAC);sig是长度为hdr.NPTS的列向量保存实际波形振幅hdr是结构体保存SAC头段里的全部元数据。读进来之后第一件事不是画图是检查头段字段名。fieldnames(hdr) % 查看关键参数 dt hdr.DELTA; % 采样间隔单位秒 t0 hdr.B; % 起始时间偏移单位秒 npts hdr.NPTS; % 采样点数 fprintf(dt%.4f s, t0%.2f s, npts%d\n, dt, t0, npts);代码里用大写字段名是常见版本的习惯如果你的rdsac返回小写字段把对应位置改成hdr.delta、hdr.b即可。dt、t0、npts分别从结构体里取出来再赋给更短的局部变量后续处理时不用反复敲hdr前缀也方便在调试时直接看这些关键值。fprintf把采样间隔、起点、点数一次打印出来和SAC头段里手动查到的数值做核对。3.2 把采样点换算成时间轴delta与b的组合画波形图不能直接用“第1个点、第2个点”做横轴需要把采样点索引换算成物理时间。rdsac读出来的b字段是起始时间的偏移量单位秒把这个偏移和delta、点数组合在一起就能得到相对时间轴t t0 (0:npts-1) * dt; figure(Color, w); plot(t, sig, k-, LineWidth, 0.5); xlabel(Time (s) from event reference); ylabel(Amplitude (counts)); title(sprintf(%s %s, hdr.KSTNM, hdr.KCMPNM)); grid on;横轴直接使用秒数需要绝对UTC时刻时把SAC头段里的参考日期和b字段合并换算startUTC datetime(hdr.NZYEAR, 1, 1) ... days(hdr.NZJDAY - 1) ... hours(hdr.NZHOUR) minutes(hdr.NZMIN) seconds(hdr.NZSEC); tt startUTC seconds(t); plot(tt, sig); xlabel(UTC Time); xtickformat(HH:mm:ss);这段代码的关键在于NZJDAY是“年内第几天”用datetime(year,1,1)先构造当年1月1日零点再加天数偏移天然处理了闰年问题比手动按每月天数累加可靠得多。绘图横轴变成datetime数组后Matlab会自动按时间刻度显示网格xtickformat控制横轴格式为时分秒。3.3 读出来的波形全是“大数”时的量纲检查做台网数据的人经常会遇到这种情况rdsac读出来的振幅不是-1到1之间的标准值而是几十万甚至上百万的量级。这通常是正常的。SAC文件里存的是数据采集器的原始计数单位是counts量程、增益、仪器响应都记录在其他头段里rdsac不会帮你做去仪器响应。处理时注意两点一是做数值分析前先确认自己需要的是counts还是物理量要转加速度或速度得结合台站灵敏度参数做反卷积二是别在counts量级下贸然做幅值比较不同台站、不同仪器的背景噪声水平可能差几个数量级。如果只是想快速看一眼事件是否清晰先把均值去掉再画图会更直观这个操作放在第4章的批处理代码里一并处理。4. 批量处理SAC数据的三分量记录与三个常见坑4.1 批量读目录下所有SAC文件真实地震事件往往不是单文件一个事件对应多个台站、每个台站又有BHZ、BHN、BHE三个分量。批量读取的第一步是用dir拿到文件列表再用rdsac逐个读取folder events/2010/001; filelist dir(fullfile(folder, *.SAC)); % 数据长度不一时使用cell数组保存避免矩阵拼接报错 D cell(numel(filelist), 1); H struct([]); for i 1:numel(filelist) filepath fullfile(filelist(i).folder, filelist(i).name); [D{i}, H(i)] rdsac(filepath); % 顺手检查填充值 if any(D{i} -12345) warning(%s contains undefined sample points, filelist(i).name); end enddir返回的filelist包含文件夹和文件两部分fullfile把目录和文件名拼成完整路径。数据用cell数组保存是为了容纳不同长度地震目录里三个分量的记录时长可能不齐直接塞进统一矩阵会触发维度错误。H(i)直接赋值给结构体数组每个元素保存一个文件的完整头段后面按台站名、分量名筛选时随时可以取用。在这个循环里顺便检查-12345填充值等于把第4.2节的坑提前暴露出来。坑现象对策填充值-12345波形出现孤立的巨大尖峰读取后查找并置NaN字段名大小写变化hdr.DELTA报错fieldnames核对实际字段名采样率不一致拼接矩阵后波形错位按台站统一网格重采样4.2 坑1SAC填充值-12345会伪造成真实振幅SAC格式规定头段字段和数据体中的无效采样点必须用-12345填充。数据体里出现-12345意味着该点无有效记录但直接读进Matlab后它就是一个普通数画图时表现为冲向图幅边缘的尖峰。更危险的是做功率谱时一个-12345会被当成真实能量污染整个频段。正确处理方式是读取后把它替换成NaN或者根据前后样本插值% 读取后立即处理填充值 data D{i}; data(data -12345) NaN; % 若后续滤波函数不支持NaN用interp1线性插值恢复 validIdx ~isnan(data); if any(~validIdx) sum(validIdx) 1 data interp1(find(validIdx), data(validIdx), 1:numel(data), linear, 0); enddata -12345生成逻辑索引把无效点全部置为NaN这样在绘图和统计时能直观看到空白区间。interp1那一段是针对后处理阶段的限制很多滤波函数遇到NaN直接报错所以在线性插值前先用find找出有效点位置用interp1按邻近有效值填补。注意插值只适用于零星坏点如果整段连续都是-12345这段记录本身就该废弃。4.3 坑2头段字段大小写随rdsac版本变化rdsac在多年流传过程中出现过不止一个维护分支不同分支对头段字段名的大小写处理不同。有的版本返回hdr.DELTA有的版本返回hdr.delta还有的版本输出是hdr.HDR.DELTA这样的嵌套结构。写脚本时如果按某个版本的习惯硬编码字段名换一台机器跑可能直接报错。避免这个问题的习惯是每次拿到新环境先执行一次fieldnames(hdr)并把输出贴到日志里或者用一个小工具函数做字段名适配function val getHdrField(hdr, name) % 同时兼容大写和小写字段名 if isfield(hdr, upper(name)) val hdr.(upper(name)); elseif isfield(hdr, lower(name)) val hdr.(lower(name)); else error(Header field %s not found, name); end end这个函数的逻辑很直白要取字段delta时先尝试hdr.DELTA再尝试hdr.delta两层都不存在就报错。封装成独立函数后主脚本里只需要写getHdrField(hdr, delta)版本差异被隔离在函数内部。团队协作或长期维护数据流水线时这种防呆写法能省掉大量换环境后的报错排查。4.4 坑3采样率不一致导致矩阵拼接错位同一个台站的BHZ、BHN、BHE三分量理论上采样率一致、时间对齐但实际数据集里经常出现某一个分量因通道故障被重采样过或者文件在归档前做了降采样而另两个没做。直接把三个向量纵向拼成矩阵强制对齐的采样点对应的物理时刻错开后续做旋转或极化分析时结果全错。统一采样率用resample函数处理% 对某一分量做重采样统一到台站公共采样率 p 1; % 重采样后采样率的分母 q 1; % 重采样后采样率的分子 % 假设目标采样率为40 Hz当前为20 Hz % 需要把数据长度扩展2倍 [p, q] rat(40 / 20); hn resample(hn, p, q); % 记下新采样率 hn.DELTA 1 / 40;resample的参数p和q是整数表示先把信号按p倍上采样再按q倍抽取。rat函数把目标采样率与当前采样率的比值近似成最简整数比后传给resample。重采样之后的头段里delta不会自动更新必须手动把hdr.DELTA改成新的采样间隔否则后面做时间轴时又会错位。这一处是批处理里最容易被忽略的隐性错误。5. 进阶rdsac读出ZNE三分量后的数据拼装与写回5.1 把BHZ/BHN/BHE拼成ZNE矩阵并做时间对齐拿到同一个台站的三分量常见需求是拼成南-北-东分量矩阵做旋转、粒子运动图或极化分析。三个分量各自读取后先做时间对齐再做矩阵拼接[z, hz] rdsac(TA.1234.BHZ.SAC); [n, hn] rdsac(TA.1234.BHN.SAC); [e, he] rdsac(TA.1234.BHE.SAC); % 对齐到垂直向的时间网格 n resample(n, hz.NPTS, hn.NPTS); e resample(e, hz.NPTS, he.NPTS); ZNE [z(:), n(:), e(:)];resample里hz.NPTS和hn.NPTS相除得到整数倍重采样比例这里隐含假设两个分量时长相同、只是采样率不同。如果时长本身就不同先做公共时间网格对齐再用interp1插值。ZNE矩阵的行是时间索引列分别对应垂直、北向、东向后续做偏振分析时直接取ZNE任意两列作为输入。5.2 保留原头段字节写回SAC文件对波形做过滤波或校正后要交回SAC生态继续处理最稳妥的写回方式是保留原始头段字节只替换数据体% 读取原始文件头段注意SAC文件是大端字节序 fid fopen(TA.1234.BHZ.SAC, rb, ieee-be); headerBytes fread(fid, 512, uint8); fclose(fid); % 写入新文件原头段 新数据 fid fopen(TA.1234.BHZ.filtered.SAC, wb, ieee-be); fwrite(fid, headerBytes, uint8); fwrite(fid, single(newData), float32); fclose(fid);fread读入的是头段原始字节不经过任何解析因此浮点区、整数区、字符区的字段偏移绝对不会出错这就绕开了手工构造头段时最容易犯的字段顺序错误。fwrite按float32写入单精度数据保持SAC数据体的标准格式。需要注意如果滤波改变了时长必须同步修改头段里的NPTS和相关字段此时建议用能正确更新头段的写SAC工具替代这种字节搬运法。写回后用rdsac把原文件和回写文件分别读进来求差信号最大值[orig, ~] rdsac(TA.1234.BHZ.SAC); [new, ~] rdsac(TA.1234.BHZ.filtered.SAC); disp(max(abs(orig - new)));如果差值接近零说明写回成功。处理后的波形与原始波形振幅谱存在刻度差异时优先检查数据体写入时是否用了float32、字节序是否与原始文件一致这两处是写回后数据异常的根源。本文还有配套的精品资源点击获取
