简介一套基于MATLAB的地震数据处理工具箱专注从互相关环境地震噪声中估计格林函数并测量地震数据中的时移面向地球物理与地震学方向的研究者和学生可用于背景噪声监测与地下介质变化分析。压缩包共含三百五十个文件涵盖地震波形数据、脚本函数、仪器响应文件、辅助脚本及说明文档等类型整体大小约三百二十三兆字节目录结构清晰。目前已有六十六人学习。工具箱实现了从数据预处理、噪声源分离、信号增强到格林函数估计与时移测量的完整流程算法成熟支持长序列噪声分析。借助数值计算与可视化能力使用者可高效提取地下介质响应并分析温度、压力及流体分布等因素引起的波速变化。配套文档与示例文件有助于快速上手与二次开发。1. 先把思路说清楚噪声互相关与时移测量做地震学研究的人基本都绕不开环境噪声互相关。以前处理台站数据时最烦的就是那些连续不断的背景振动风、海浪、车流、工厂机器它们把地震事件波形搅得一团糟。后来被动源成像逐渐普及我们才发现这些“没用的噪声”里其实藏着一条稳定可信的传播路径信息——把两个台站的长时间记录做互相关得到的波形竟然接近二者之间传播的经验格林函数。基于这个思路不少课题组会开发自己的Matlab工具箱用来从互相关环境地震噪声中估计格林函数并在不同时间段之间测量地震数据中的时移从而推断地下介质波速的微弱变化。这套方法最大的价值在于“被动”。不需要人工震源也不用等天然地震只要台站持续记录就能通过噪声互相关获得台站之间的传播信息。加上地震仪长期连续观测的成本相对可控所以非常适合做重复性强、周期长的监测工作。对于做背景噪声成像、地壳介质物性监测、火山活动跟踪或者地下储层变化研究的人来说一个顺手的Matlab工具箱能省掉大量重复劳动。新手用它快速上手上手老手则可以通过改参数扩展自己的科研流程。下面我按自己的使用习惯把这类工具箱从原理到实操逐层拆开讲。1.1 为什么互相关能“变废为宝”先解决一个最基本的问题两个台站记录的噪声凭什么能互相“勾兑”出格林函数这里面的理论假设是扩散场模型。当噪声源在空间上足够均匀并且长时间随机激励出复杂的散射波场时任意两个接收点之间的互相关函数会趋近于两点之间格林函数的时间导数形态。换句话说噪声场里已经包含了传播路径的响应信息只是被纷乱的随机信号掩盖了互相关运算就是把这些随机相位抹掉把稳定的路径响应刷出来。我用一个生活场景做类比。你在一个大会场两端各放一支录音笔现场人来人往、说话声混乱单独听某一支录音什么都没有但如果把两支录音笔的长时间录音做互相关就能大致提取出声波从会场一头传到另一头的脉冲响应。哪怕每个人说话都是随机的足够长的公共声场信息也能把这条路“刷”出来。实际地震台站记录里的噪声源主要是海洋波浪、风和人类活动虽然不满足严格的均匀分布但通过长时间叠加、谱白化、时域归一化等预处理可以在很大程度上逼近扩散场条件。拿到经验格林函数之后需要注意一个问题它和理论格林函数的相位信息基本可靠但振幅并不是绝对尺度。因为噪声源的强度、方向和分布会直接影响互相关波形振幅所以做时移测量比做绝对衰减研究更可靠。这也是为什么这类工具箱通常把重心放在“测时间延迟”而不是“测绝对振幅”上。1.2 时移测量到底在测什么“时移”这个名字听着抽象翻译成人话就是对比两段不同时期的互相关波形看看同一组台站对之间的传播波形是不是整体提前或滞后了。假如地下介质波速下降波走完同一段距离需要更长时间互相关波形就会整体向后延迟反过来介质硬化或波速升高波形就会提前。我们把这个相对时间变化写下来通常用 dt/t 表示它和波速相对变化 dv/v 近似满足dt/t ≈ -dv/v也就是说测出千分之一的走时变化基本等于测出负的千分之一波速变化。别小看这个量级地壳介质在地震前后、火山活动期、降雨入渗或地下水开采过程中波速变化通常只有百分之零点几甚至更低普通走时分析根本分辨不出来。噪声互相关的方法优点在于重复性特别好参考波形相对固定通过长时间窗口的统计平均可以把噪声压制得很低从而稳定分辨极小的时移。实际应用里这个技术已经被用到地震后断层愈合监测、火山喷发前兆观测、地下水储量变化跟踪等场景。比如一次强震后断层带附近破碎导致波速明显下降之后随着应力恢复和裂缝愈合波速又慢慢回升时移曲线会看到一条“下降—回升”的弧线。这类信号用传统走时方法很难稳定捕捉但借助噪声互相关和基于Matlab的时移测量工具箱可以变成一套比较成熟的日常分析流程。1.3 一个工具箱该具备的模块骨架我见过很多版本的噪声互相关Matlab工具箱模块划分五花八门但核心骨架基本一致。我整理了一个比较通用的模块表方便你对照自己手里的代码或准备入手的工具箱模块主要任务关键输入/输出数据读取模块读取SAC、miniSEED等连续波形处理通道和元数据原始数据文件整理后的波形矩阵预处理模块去均值、去趋势、带通滤波、时域归一化、谱白化连续波形干净可用的波形片段互相关计算模块分段互相关、单日叠加、多日叠加、信噪比计算两段波形互相关函数序列时移测量模块移动窗口互相关或拉伸法估计相对时移并评估误差参考波形和目标波形dt/t、dv/v曲线可视化与输出模块绘制互相关热图、时移曲线输出参数表中间结果最终图形和文件模块化的好处是便于替换和调试。比如数据读取格式变了只需要改第一层又想试另一种波速变化估计方法不需要动前面的互相关计算直接换掉时移模块就行。这也是我推荐大家在接受别人的工具箱时先按这五个模块梳理一遍代码结构的原因否则后面改参数容易越改越乱。2. 工具箱核心功能参数背后的门道工具拿到手最容易踩坑的是“什么参数都敢改改完不知道发生什么”。实际上这套流程的每个核心参数都有明确物理含义弄懂了原理参数选择就不再是拍脑袋。2.1 预处理参数先去掉“明显的错误”预处理阶段的目标不是让波形变好看而是把非平稳、非噪声源的瞬态干扰压下去尽量让剩余记录逼近随机平稳噪声场。最基础的是去均值和去趋势目的是消除仪器直流偏置和长周期衰减。紧接着是带通滤波这一步决定了你关注的是哪个频段的传播波。陆地环境噪声的长周期部分以海洋波浪的次生微震为主通常在0.05~0.3 Hz短周期部分则受风和人类活动影响大比如1~10 Hz。跑长时间地壳波速监测时我比较常用0.02~0.2 Hz或0.1~1 Hz两个频段交叉验证避免单一频段被异常噪声源污染。时域归一化和谱白化是预处理里最有“手艺”的部分。时域归一化最狠的一招是“one-bit”也就是把每个采样点直接变成1或-1只保留符号不保留振幅。这一招能把地震事件、仪器脉冲等大幅瞬态信号直接压平但代价是丢失振幅信息。谱白化则是在频率域做归一化让不同频带振幅相对均衡避免某个频带特别强的噪声源主导最后结果。一般建议顺序是去均值、去趋势、带通滤波、谱白化、再做一次时域归一化。具体组合要看原始数据质量没有通吃方案。实际操作中有个容易忽视的细节滤波器不要一味追求高陡度。高阶级联滤波器虽然频带边界更干净但容易产生振铃效应把本不存在的周期信号“造”出来。我更推荐用二阶到四阶巴特沃斯滤波器配合零相位滤波也就是Matlab里的 filtfilt 思路避免相位偏移影响后面的时移测量准确度。2.2 互相关与叠加参数信噪比是逼出来的互相关阶段最关键的参数是分段长度、叠加方式和信噪比控制。分段长度一般选一天原因很现实一天的连续记录既能保证足够长的噪声平均又能提供较好的时间分辨率方便逐天追踪波速变化。如果时间分辨率要求更高可以压到几小时一段但相应地单段信噪比会下降如果想提高信噪比可以把多天甚至整月数据叠成一条参考波形牺牲时间分辨率换稳定性。叠加方式上线性叠加最直接也最常用直接把每天的互相关函数求平均。问题在于如果某一天出现强干扰或大幅异常这一天的互相关会被拉偏从而污染整个叠加结果。所以我更倾向于用相位加权叠加PWS把瞬时相位一致性好的样本赋予更高权重在一定程度上压制异常样本。需要提醒的是PWS并不是越用越好如果数据信噪比本身很高线性叠加和PWS结果差别不大没必要为了“高级”而强行加权。互相关函数的信噪比一般这么算先在正负延时方向找到信号窗口比如台站间距除以面波速度所在的延时区间再取远离信号窗口的“纯噪声”区域用两个区域能量比作为评估指标。低于某个阈值时这一天结果就别纳入叠加了硬加进去只会拉低整体质量。另外实际噪声源分布往往不是均匀的互相关函数正负两个分支经常不对称这是正常现象不要强行把两支平均成完美对称否则可能掩盖真实的传播信息。2.3 时移估计参数MWCS方法细节时移测量我自己最常用的是移动窗口互相关Moving Window Cross-SpectralMWCS方法。思路很朴素把参考波形和目标波形切成多个短窗口每个窗口分别做互相关得到该窗口对应的局部时移再把所有窗口的时移结果沿窗口中心时间做线性拟合斜率和截距就对应整体的相对时移和频率相关的色散信息。这里的窗口长度和步长非常考验经验。窗口太短单个互相关结果被噪声主导局部时移跳动剧烈窗口太长又等于把不同时间段的信号平均在一起丢失了空间分辨率。通常做法是以中心周期的5~10倍为窗口长度下限比如工作在0.1 Hz附近主周期约10秒窗口长度可以选100~200秒步长则选窗口长度的50%左右。最大允许时移一般设置为窗户长度的一半超过这个范围的结果直接归为异常。设置完这些参数后算法会输出每个窗口的时移量、相关系数和误差最后按相关系数或误差阈值剔除不靠谱的点再做一次稳健回归得到最终的 dt/t 估计。还有一类方法叫拉伸法Stretching原理是对目标波形施加时间轴的压缩或拉伸寻找与参考波形相关系数最大的伸缩因子然后由伸缩因子直接换算 dv/v。两种方法在不同信噪比下表现不同MWCS对高频成分更敏感拉伸法对窄带信号更稳。成熟工具箱一般同时保留两条路线方便结果交叉验证。我只强调一点无论是哪种方法参考时段的选择非常关键。尽量选介质状态相对稳定、波形信噪比最高的时段如果把大震后波速快速变化的那段当作参考后续时移很可能出现整体偏移误判为匀速漂移。3. 实操流程从原始数据到时移曲线理论再清楚最后还是得落在“能跑出图”上。这一节我用一套典型的Matlab流程把从原始数据到时移曲线的关键步骤串起来。3.1 环境准备与数据整理环境准备的第一步不是写代码而是把数据整理干净。Matlab版本建议R2016b以上因为很多工具箱代码会用隐式扩展和新式语法。除了Matlab本身需要确认信号处理工具箱、并行计算工具箱已经安装。顺带提醒一句网上搜索“Matlab工具箱”时很容易混进各种硬件检测、系统优化相关的名字比如图吧工具箱之类跟地震数据处理完全不是一回事别下错东西。我们这里说的工具箱本质是一批把读写、滤波、互相关、时移测量封装成函数的脚本集合用addpath添加目录后就能调用。数据格式上国内很多台阵数据是SAC格式或miniSEED格式。Matlab读取SAC比较常见可以直接用rdmseed或SAC工具箱读成结构体如果要处理长时间连续记录最好先用文件清单把每天的SAC文件对应到连续时间轴上。元数据至少要有采样率、台站经纬度、通道方向和起止时间后续计算台间距、选择延时窗口都需要这些信息。我习惯在工作目录下建三个子文件夹raw_data放原始记录process存放每日互相关中间结果figures输出图形。这样批量调试时不会把不同阶段的文件搞混。3.2 一个可落地的Matlab示例流程下面这段代码是示意框架实际运行时需要根据自己的数据接口替换读取函数但整体流程可以直接照搬。第一步是设置参数并读取两台站垂直分量波形% 参数设置 fs 20; % 采样率单位Hz按实际数据修改 freqRange [0.02 0.2]; % 关注频带 maxLag 300; % 最大互相关延时单位秒 winLen 86400 * fs; % 分段长度一天 % 读取两台站连续波形 [d1, meta1] readMySAC(TA.001..HHZ.sac); [d2, meta2] readMySAC(TA.002..HHZ.sac); % 基础预处理 d1 detrend(d1 - mean(d1)); d2 detrend(d2 - mean(d2)); d1 butterworth_bandpass(d1, fs, freqRange); d2 butterworth_bandpass(d2, fs, freqRange);这里用到了自定义函数 readMySAC 和 butterworth_bandpass前者负责读文件并返回波形和元数据后者用零相位滤波器实现带通。接下来分段做互相关并把所有天的结果存进一个矩阵% 分段互相关并叠加 numSeg floor(min(length(d1), length(d2)) / winLen); ccDaily zeros(numSeg, 2 * maxLag * fs 1); for i 1:numSeg idx (i-1)*winLen 1 : i*winLen; seg1 time_normalize(d1(idx)); % 时域归一化 seg2 time_normalize(d2(idx)); [cc, lag] xcorr(seg1, seg2, maxLag*fs, coeff); ccDaily(i, :) cc; end ccStack mean(ccDaily, 1);时域归一化函数 time_normalize 可以用滑动绝对均值法也可以用one-bit符号化我建议先试one-bit看波形稳定后再换更温和的方法。得到每日互相关矩阵之后选择前一段时间作为参考时段然后对后续每天做移动窗口互相关% 移动窗口互相关估计时移 refStack mean(ccDaily(1:30, :), 1); for k 31:numSeg [dtAtWin, dtErr, coef] mwcs(refStack, ccDaily(k, :), ... Fs, fs, ... FreqRange, freqRange, ... WinLen, 120, % 窗口长度单位秒 Step, 30); % 步长单位秒 dvvCurve(k-30) mean(dtAtWin, omitnan); % 简单平均得到近似dvv end这里 mwcs 是封装好的函数返回每个窗口的局部时移、误差和相关系数。日常分析时不要把 dvvCurve 直接当成最终结果至少先看一眼 coef 序列把低相关结果筛掉。这个流程跑通后中间结果建议先保存为 .mat 文件后续画图就不要再重新读原始波形了。3.3 结果怎么看从波形图到时移曲线拿到互相关矩阵和时移曲线后的第一件事不是急着解释物理意义而是先看波形形态是否符合预期。我习惯画一张“互相关热图”横轴是延时纵轴是日期颜色表示互相关振幅。信噪比高的台站对能看到在正负延时方向各有一条清晰的同相轴位置大致对应台站间距除以瑞利波群速度。如果整张图都是散粒噪声没有任何稳定同相轴多半是频带选错、数据缺失或者台站间距过大。时移曲线一般画成“横轴时间、纵轴dt/t或dv/v”的带误差条散点图。质量好的结果应呈现均匀波动误差条小且不同频率段的结果趋势一致。如果曲线像楼梯一样一段一段跳很可能是数据分段之间拼接问题或钟漂如果误差条异常大说明当前频带的信噪比不足需要考虑缩小频带或延长叠加天数。时间分辨率与信噪比之间存在取舍比如想监测降雨引发的地表波速变化可能需要用几天滑动窗口叠加来保证信噪比代价是变化事件的起止时间被模糊。跑数据时一定要在脚本头部写清楚“用了几天的叠加窗口、参考时段是哪几天”否则几个月后再回看参数全忘干净了。4. 常见问题与排查技巧不管是新手还是老手跑这类流程都免不了遇到“结果很丑”的时刻。我把自己遇到过的典型问题整理成一张速查表再补充几个排查思路希望能帮你少走弯路。4.1 互相关结果一团糟先排掉这四类坑现象常见原因排查与解决互相关图看不到清晰同相轴频带选择不匹配尝试更低频段如0.005~0.1Hz或更高频段0.1~1Hz每日互相关序列不连续数据缺失、触发事件没剔除检查SAC文件时间轴剔除低信噪比时段只在正延时或负延时分支有信号噪声源方向性太强属于正常现象不做对称平均只用稳定分支波形形态每天剧烈变化未做谱白化或归一化方法不合适增加谱白化先用one-bit测试稳定性这里最容易被忽略的是数据缺失。SAC文件本身可能存在时间轴不连续、采样点跳变或零值填充如果不检查直接进入互相关会在特定延时位置引入假信号。我自己习惯在每个预处理步骤后都把波形画出来扫一眼哪怕只是看一秒钟的局部波形也能发现很多“看起来能跑但实际有毒”的问题。4.2 时移曲线“假漂移”怎么识别时移曲线出现大范围线性漂移时不要立刻联想到地下介质变化先排查设备问题。时钟漂移是最常见的一种假漂移如果所有台站对的时移曲线都呈现同步线性增长而且增长速率大致相同大概率是记录仪内部时钟与GPS时间失步导致每天波形整体被“压扁”或“拉伸”。这种时候可以拿某个已知走时的地震事件做交叉验证直接用天然地震的P波或S波到时检查时钟偏差。另一种假漂移源是参考时段选得不合适。如果参考波形本身包含地震破裂后的异常阶段后续曲线会整体出现恒定偏差看起来像匀速漂移。此外温度变化也可能造成仪器电子元件的相位响应变化这种情况下的时移曲线往往带有明显季节性和当地气温曲线高度相关。排查时可以把不同频段的时移结果叠在一起看如果不同频段趋势不一致更大概率是设备和参考问题而不是地下介质各向异性变化。4.3 大数据量下的效率经验噪声互相关看起来只是两个波形相乘求和但如果台站多、时间长计算量会迅速膨胀。一百个台站就能组合出近五千个台对每个台对每天处理几条数据累积起来非常可观。我的经验是先做一次数据检查把所有台站的互相关计算任务划分成独立模块优先用 parfor 并行循环。注意 parfor 里尽量减少共享内存变量把每天结果单独写到一个临时单元循环结束后再统一合并否则内存反复拷贝反而拖慢速度。中间结果及时落盘也是一条重要经验。计算单日互相关可以很快但后期叠加和时移测量都需要反复读取这些中间结果与其每次重新读SAC重新滤波不如在第一次跑完后就把 ccDaily 矩阵按台站对存成 .mat 文件。为了节省磁盘空间和内存可以在保存前转换为 single 类型精度损失对时移分析影响很小。还有一点先拿一个台站对跑通全流程确认参数和代码逻辑没问题再批量启动所有台站对。批量跑了半天之后才发现滤波参数写错这种错误很常见也非常浪费时间。4.4 一个值得保留的习惯最后分享一个我自己的固定习惯拿到任何新的噪声互相关工具箱我不会直接跑全流程而是先挑一个台站对、一个窄频带、一天的数据把“读取—预处理—互相关—时移”整条链路跑一遍确认每一级输出的图形、尺寸、量纲符合预期。这个过程通常只需要十分钟但能避免后面大量无效计算。我见过有人花一天时间把上百个台站对全跑完最后才发现参考波形选错了所有结果都要重算。数据处理的耐心往往体现在这些看起来“多余”的检查上。参数记录也非常重要我会在每次运行的脚本头部用注释写清楚频带、归一化方式、参考时段和叠加窗口方便几周后回看时还能还原当时的处理逻辑。本文还有配套的精品资源点击获取