MATLAB读取SAC地震数据：rdsac.m脚本实现与实战

简介：一个用于处理 SAC 格式地震数据的 MATLAB 脚本包，面向地震学、地球物理学领域的科研人员与技术人员，解决在 MATLAB 环境中直接读取和分析 SAC 文件的需求。压缩包内共 1 个文件，为 rdsac.m 脚本，体积仅 1KB。该脚本围绕 SAC 数据格式封装了二进制读取、时间序列转换、文件头元数据解析与基础错误检查等功能，输出 MATLAB 友好数据结构，便于后续滤波、频谱分析、事件定位等二次开发。已有 247 人学习或下载这份资源，适合具备一定 MATLAB 基础、希望将 SAC 数据导入统一分析流程的初学者快速上手。借助该脚本，用户可省去手动解析二进制格式的繁琐步骤，将更多精力集中在信号处理与可视化上，从而提升地震数据处理效率。

1. 一个 rdsac.m 脚本为什么值得留着

处理 SAC 格式地震数据的人,几乎都遇到过同一个尴尬:仪器厂家给的文件是 SAC,手里最顺手的工具却是 MATLAB。直接load肯定不行,textread更不可能,因为 SAC 是二进制格式,数据结构由固定的头部区和波形数据区组成。rdsac.m 这类脚本的价值就在于把这段二进制解析封装成一次函数调用,让 MATLAB 用户不用关心字节序、头部偏移和数据类型,直接拿到时间序列数组和采样率。这篇博文我会从 SAC 的二进制布局讲起,给出 rdsac.m 的完整实现思路,再落到批量处理台站数据、排查读取失败和做简单事件检测的具体用法。适合刚接触 SAC 格式的 MATLAB 用户,也适合写过读取代码但没见过完整边界情况的工程师。

2. SAC 文件头的二进制布局与 rdsac.m 的解析策略

2.1 为什么不能跳过文件头直接读波形

SAC 文件由两部分构成:固定长度的头部区(Header)和紧随其后的波形数据区。头部区的长度在 SAC 规范中是固定的 632 字节,这 632 字节又被细分为三类字段:数值型字段(浮点与整型)、逻辑型字段(字符)和文本型字段。对 rdsac.m 来说,真正影响数据读取的参数都集中在数值型字段里:delta(采样间隔,秒)、npts(采样点数)、b(起始时刻,秒)、e(结束时刻)、stla/stlo(台站经纬度)、kcmpnm(分量名称)等。

如果跳过文件头直接读数据,连最基本的npts都不知道,自然无法确定要读多少个采样点。更麻烦的是,不同仪器产出的 SAC 文件可能使用不同的字节序,少数非规范文件甚至会截断头部。这些信息只能从头部字段中判断,所以 rdsac.m 的第一步动作永远是解析头部,而不是急着取波形,这是所有 SAC 读取工具的共同逻辑。

2.2 头部字段顺序与字节偏移

SAC 头部的排列顺序是有规范的,首先是 70 个浮点型字段,然后是 15 个整型字段,接着是 20 个字符型字段,最后是 176 字节的文本型字段。由于 MATLAB 的fread按顺序读取二进制流,解析头部时只要按照这个顺序逐个读出来即可,不需要手动计算每个字段的绝对偏移。为了让你对映射关系有个直观印象,我把它压缩成一张常用字段表:

注意表格中的“相对位置”是从 0 开始计数的。rdsac.m 在读取时必须严格保持这个顺序,否则解析出来的delta可能落到别人的stla上,这种错误极难排查。

2.3 rdsac.m 的解析策略:一次性读完还是分步读

常见的做法是在 MATLAB 中分三次fread:第一次读 70 个 float 和 15 个 int,第二次读 20 个字符字段,第三次读剩余的文本区。这样做的理由是分步读取时字段定位不容易出错,代码也更接近 SAC 规范本身的描述。我倾向于把 float 和 int 一次读完,因为它们加起来才 85 个数值,开两个临时数组分别存,可读性并不差。顺序上先判断文件是否包含正确的npts,再决定是否继续读取数据区。

这里还有一个细节:MATLAB 默认按 IEEE 小端方式解析二进制,而部分 UNIX 工作站产出的 SAC 文件是大端。rdsac.m 内部一般会提供字节序参数,默认值可以设成小端,遇到异常数据时再切换。字节序的判断依据通常不是文件扩展名,而是头部npts读出来是否合理,比如是否为正数、是否小于文件大小除以数据字长。

3. rdsac.m 的读取实现:从 fopen 到时间序列数组

3.1 核心读取代码与数据结构设计

下面这份代码是我按 rdsac.m 的思路重写的可运行版本,兼容 MATLAB R2016b 之后的绝大多数版本。它把文件路径作为入参,输出一个包含头部字段和时间序列的结构体,这也是地震数据读取工具最通用的约定。

function sac = rdsac(filename, byteorder) % rdsac 读取 SAC 格式地震数据文件 % 输入: % filename - SAC 文件路径 % byteorder - 字节序,'ieee-le' 或 'ieee-be',默认 'ieee-le' % 输出: % sac - 结构体,包含头部字段和时间序列 if nargin < 2 byteorder = 'ieee-le'; end fid = fopen(filename, 'rb', byteorder); if fid == -1 error('rdsac:fileOpenFailed', '无法打开文件: %s', filename); end % 读取数值型头部:70 个 float + 15 个 int float_head = fread(fid, 70, 'float32'); int_head = fread(fid, 15, 'int32'); % 读取字符型头部:20 个 8 字节字符串 char_head = fread(fid, 20 * 8, 'uint8')'; char_head = char(char_head); % 读取文本型头部,直到头部区结束 text_body = fread(fid, 176, 'uint8')'; text_body = char(text_body); % 跳过头部位剩余部分(如有),再读取波形数据 data = fread(fid, inf, 'single'); fclose(fid);

代码先把 70 个 float 字段读入float_head,再把 15 个 int 字段读入int_head,最后用fread(fid, inf, 'single')读取全部剩余数据。这样头部解析和波形提取共用同一个文件句柄,避免了多次打开文件带来的句柄管理问题。这里的关键在于fread的类型参数必须与 SAC 规范严格一致:float32对应 SAC 的浮点字段,int32对应整型字段,single对应波形数据的 4 字节浮点。如果 MATLAB 版本较老,不支持float32写法,可以用'float32'的等价形式'single'代替。

3.2 字段归档与时间轴重建

头部解析出来是一堆裸数组,直接丢给用户并不友好。继续处理:把 float 字段按 SAC 规范的顺序归档,并重建时间轴。

% 从数值头部中提取关键字段 sac.delta = float_head(1); sac.b = float_head(6); sac.e = float_head(7); sac.stla = float_head(32); sac.stlo = float_head(33); sac.stel = float_head(34); sac.npts = int_head(10); % 从字符头部中提取台站名与分量名 sac.kstnm = strtrim(char_head(1:8)); sac.kcmpnm = strtrim(char_head(7*8-7:7*8)); % 数据长度检查 expected = round((sac.e - sac.b) / sac.delta) + 1; if abs(expected - sac.npts) > 1 warning('rdsac:nptsMismatch', '头部 b/e/delta 推算点数与 npts 不一致'); end % 只保留有效数据长度 data = data(1:min(sac.npts, numel(data))); % 重建时间轴 sac.t = sac.b + (0:sac.npts-1) * sac.delta; sac.data = data;

这里有几个容易忽略的细节:字符字段在 SAC 规范中是 8 字节定长字符串,不足部分用空格填充,所以读取后必须strtrim去掉首尾空格,否则台站名里会带一串空字符,后续做字符串比较时永远匹配不上。npts的检查也很有必要,因为它反映了一个常见情况:头部声明的采样点数与实际数据区长度不一致,这通常由文件截断或写入错误引起。兼容的做法是取两者中的较小值,而不是直接报错,这样至少能把能用的那部分数据读出来。

3.3 数据长度与采样率的边界场景

SAC 文件中的npts是按“实际写入的采样点”计的,但有些 pre-event 数据会在b为负值的情况下记录,起始时刻偏移并不影响采样点数。rdsac.m 中重建时间轴时直接以b + (0:npts-1) * delta计算,即便b为负数也不会有问题。另一个常见边界是delta为 0 或负数,这种文件多半是损坏的,遇到时应该在读取阶段就抛错,否则后续绘图时横轴会出现无穷大。我在实际处理台阵数据时见过delta被写成 0.0 的文件,症状是所有波形挤成一条竖线,排查了很久才发现是头部位错误,所以 rdsac.m 里值得加一个判断:如果delta <= 0,直接报错提示文件损坏。

4. 用 rdsac.m 批量处理台站数据的实战路径与参数调优

4.1 批量读取目录下所有 SAC 文件

单文件读取只是基础,实际项目里最常见的诉求是批量处理一个台站目录下连续多天的记录。以下脚本使用dir枚举目录中的全部 SAC 文件,再逐条调用 rdsac.m,把返回的结构体存入 cell 数组。这样后续无论是画波形还是提取振幅,都不用再关心文件打开与关闭的逻辑。

% 批量读取示例 sacDir = '/data/seismic/2025/091'; files = dir(fullfile(sacDir, '*.SAC')); nFiles = length(files); sacList = cell(nFiles, 1); for i = 1:nFiles filePath = fullfile(sacDir, files(i).name); try sacList{i} = rdsac(filePath); fprintf('[OK] %s npts=%d delta=%.3f\n', ... files(i).name, sacList{i}.npts, sacList{i}.delta); catch ME fprintf('[FAIL] %s %s\n', files(i).name, ME.message); sacList{i} = []; end end

对批量处理来说,加一层try-catch非常值得:某个文件损坏时,脚本可以跳过它继续处理其余文件,同时打印出失败文件名,方便事后单独排查。这里sacList{i} = []的意义在于保持 cell 数组长度不变,后续循环读取时不会因为缺少元素而报 index 越界错误。delta打印出来可以直接检查采样率是否异常,比如短周期仪器写出的 delta 通常是 0.01 秒,长周期或强震仪可能是 0.005 或 0.1 秒。

4.2 分析参数与输出结构的选择

批量处理完成后,紧接着的操作往往是筛选有效台站或去除坏道。rdsac.m 输出结构体里的npts、delta、stla、stlo这组字段,正好可以作为筛选条件。常见的筛选维度有两个:时间跨度是否满足事件分析窗口,以及台站坐标是否落在目标区域范围内。我习惯把筛选逻辑抽成一个独立的函数,输入sac结构体,输出布尔值,而不是在读取循环里写一堆if判断,这样后续调整阈值只改一个地方。

% 筛选条件示例:采样率不低于 50Hz,时长不低于 30 秒 valid = false(nFiles, 1); for i = 1:nFiles if isempty(sacList{i}), continue; end fs = 1 / sacList{i}.delta; dur = sacList{i}.t(end) - sacList{i}.t(1); valid(i) = (fs >= 50) && (dur >= 30); end selectedFilenames = {files(valid).name};

这段代码里fs是采样率,由delta取倒数得到;dur由时间轴首尾相减得到,和npts / fs等价,但用时间轴算可以规避npts头部异常的情况。筛选逻辑单独写成valid数组后,与文件列表一一对应,无论是画在平面图上还是写入文本清单都方便。

4.3 绘图与质量检查结合

批量读取之后,可视化是最快的质量检查手段。下面代码把所有文件的波形画在同一张图上,纵轴按文件名排列,一眼就能看出哪些通道是平的、哪些有尖峰。

figure('Color', 'w'); hold on; offset = 0; for i = 1:length(sacList) if isempty(sacList{i}), continue; end s = sacList{i}; plot(s.t, s.data + offset, 'LineWidth', 0.6); offset = offset + max(abs(s.data)) * 2; end xlabel('Time (s)'); ylabel('Amplitude + offset'); title('Batch SAC waveform viewer');

这里的offset每次累加当前道最大振幅的两倍,目的是让各条波形从上到下依次排开,避免互相重叠。如果不加偏移直接叠合绘图,振幅小的近震波形会被远震大振幅完全盖住。需要注意plot(s.t, s.data + offset)中s.t和s.data必须等长,而 rdsac.m 在数据长度不匹配时已经做了截断处理,所以这里不会报维度错误。

5. 常见读取失败场景与 rdsac.m 的兼容性边界

5.1 字节序错误是最隐蔽的问题

字节序导致读取失败的案例在地震数据处理中属于高频问题。SAC 规范允许按小端或大端存储,而 Linux 和 Windows 平台默认都是小端,但部分 Solaris 或 AIX 时代的老工作站会把 SAC 写成大端。rdsac.m 默认按小端解析,遇到大端文件时npts会变成一个巨大的正数或负数,进而导致数据区读取异常。判断方法很简单:用默认字节序读取,如果npts大于 100000000 或为负数,更换字节序再读一次,结果一般就正常了。

% 自动判断字节序 try s = rdsac(filename, 'ieee-le'); if s.npts <= 0 || s.npts > 1e8 error('rdsac:badNpts', 'npts invalid, retry with big-endian'); end catch s = rdsac(filename, 'ieee-be'); end

这种“先小端失败再大端重试”的策略在批量读取老数据时很实用,代价只是第一次读取可能白读了几十毫秒,但避免了人工逐个文件试错。需要说明的是,SAC 头部里没有明确的字节序标志位,所以任何自动判断都是基于字段合理性推断,纯启发式逻辑。如果你想在 rdsac.m 里直接实现这个逻辑,可以在解析头部后立即校验npts的数值范围。

5.2 头部字段缺失与数据截断的处理

在实际的 SAC 文件中,头部字段缺失的情况比想象中常见。例如,部分弱震仪器写入的 SAC 文件把b、e置为 0,只在数据区里存波形,此时时间轴的起点就无法忠实还原。rdsac.m 对这种文件不应该直接报错,而应以b=0为默认起点继续读取,把判断交给上层调用者。还有一种更普遍的问题:数据区长度小于npts声明长度,这通常由网络传输中断或磁盘写满导致。我在代码里用data(1:min(sac.npts, numel(data)))截断数据,并把实际长度输出供用户检查,比直接抛异常更友好。如果整个目录里出现多个文件同时截断,基本可以判断是采数软件批量写入时出了问题,而不是单个文件损坏。

5.3 读取结果的正确性验证

读完了不代表读对了。我建议用两种方式交叉验证 rdsac.m 的结果:一是对比文件大小与npts的关系,二是将读取结果与 SAC 官方工具saclhdr的输出对比头部字段。文件大小校验的逻辑是:

% 验证文件大小与 npts 是否匹配 fileInfo = dir(filename); headerLen = 632; dataBytes = fileInfo.bytes - headerLen; expectSamples = dataBytes / 4; % 4 字节浮点 fprintf('npts=%d, 按文件大小推算=%d\n', sac.npts, expectSamples);

如果npts与按文件大小推算出来的采样点数不一致,说明头部或数据区有异常。用saclhdr对比时,重点关注delta、npts、stla、stlo四个字段,只要这四个字段一致,读取结果就可以放心用于后续分析。需要注意,SAC 文件可能含有多个通道数据,此时头部字段npts表示单个通道的采样点数,文件大小与npts的关系还取决于通道数,不能用上面的简单公式直接套。

6. 结合 rdsac.m 做事件检测的快速验证技巧

6.1 用读取结果直接计算 STA/LTA 特征函数

SAC 数据读取成功后,最值得先做的事是做一个简单的 STA/LTA 事件检测。STA/LTA(Short-Term Average / Long-Term Average)是地震学里最经典的事件触发算法,原理是短时窗振幅均值与长时窗振幅均值的比值超过阈值时判定为事件。rdsac.m 返回的t和data正好构成 STA/LTA 需要的完整时间序列输入。

% STA/LTA 事件检测示例 s = rdsac('/data/seismic/tail.sac'); dt = s.delta; x = s.data; % 去除均值,避免直流偏置影响 x = x - mean(x); % 短窗与长窗长度 nsta = round(0.5 / dt); nlta = round(10.0 / dt); % 计算特征函数(绝对振幅) cf = abs(x); % 滑动均值方式计算 STA 与 LTA sta = movmean(cf, nsta); lta = movmean(cf, nlta); % 添加极小值防止除零 ratio = sta ./ (lta + eps); % 用阈值检测并记录触发时刻 thr = 3.0; triggerIdx = find(ratio > thr);

这段代码里的movmean是 MATLAB 内置滑动均值函数,比手写循环快很多,而且自带边界填充处理。nsta和nlta是根据采样率动态计算的:短窗取 0.5 秒,长窗取 10 秒,这个配比适合短周期近震事件。如果数据是远震,短窗可以放到 1 秒以上,长窗 20 秒左右,STA/LTA 的响应会平滑很多。

6.2 触发时刻与波形对齐验证

检测到triggerIdx后,需要把索引换算成绝对时间,再与原始波形对比验证。这里直接用s.t(triggerIdx(1))就能得到第一个触发时刻对应的秒数。实际处理中,一次事件的触发索引可能连续有好几百个,常见做法是取第一个索引作为事件起点,然后向后搜索一个不小于短窗长度的空白段作为事件终点。验证时叠加绘制波形和 STA/LTA 曲线,确认高比值区间与实际波形高振幅段吻合。

% 找到连续触发段的终点 gap = round(1.0 / dt); trigStart = triggerIdx(1); endIdx = trigStart; for i = 2:length(triggerIdx) if triggerIdx(i) - triggerIdx(i-1) > gap break; end endIdx = triggerIdx(i); end % 输出事件时间窗口 fprintf('事件窗口: %.2fs - %.2fs\n', s.t(trigStart), s.t(endIdx));

这种事件检测思路不追求精确到事件 P 波或 S 波到时,而是作为数据质量预筛查,在批量处理中快速挑出包含有效震相的 SAC 文件。rdsac.m 读取出的数据结构与movmean的组合,让整个过程只用了十几行 MATLAB 代码,却覆盖了从文件解析到事件识别的完整链路。如果你在真实台站数据上触发率偏低,优先检查delta是否正确、数据里是否有大量尖峰噪声,这两个因素最容易干扰 STA/LTA 比值。

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