简介pickmt是一套基于MATLAB实现的地震矩张量相位拾取与反演工具面向从事震源机制分析的地球物理研究者及高年级本科生、研究生核心目标是利用三分量地震记录拾取P波和S波到时并通过反演获得震源矩张量的六个独立分量以刻画地震破裂的力学性质。压缩包共含31个文件其中26个.m源代码组织了相位拾取、合成地震图生成、滤波绘图、矩张量反演等主要环节2个.mat文件提供示例地震波形数据另有2个txt说明文档和1个README帮助用户快速上手整包大小约30.57MB。相位拾取模块负责从波形中自动或半自动标记波至反演模块则采用数值优化方法拟合观测数据代码流程清晰便于追踪每一步计算。目前已有728人学习下载适合具备一定MATLAB基础和地震学背景的读者参考使用。使用者拿到的是完整源码与示例数据可直接运行复现处理流程也可以根据实际观测资料修改输入格式或扩展反演目标函数具备较好的二次开发空间。 做微震监测或者诱发地震研究的朋友应该都经历过这种场景晚上拖着高信噪比波形看半天只为手动拾几个P波和S波的到时好不容易拾完了还得导到专门的矩张量反演程序里重新调参数一折腾就是半宿。pickmt这套Matlab代码把面向矩张量反演的相位拾取和反演合成一条流水线从原始三份量波形到最终的震源机制解也就是大家常说的沙滩球图一个流程跑完。它重点解决的正是我日常处理数据时最头疼的两件事震相拾取的主观性以及反演流程的割裂感。如果你正在做微震监测、诱发地震、矿区冲击地压或者近震震源机制研究或者你是刚接触矩张量反演方向的研究生这篇文章都值得花十分钟看完。我会把相位拾取的原理、反演中的关键参数、以及实际跑数据时容易踩的坑一次讲清楚。1. 项目定位为什么要把拾取和反演绑在一起1.1 传统流程的痛点常规做法是两步走先拾取震相再反演矩张量。第一步用人工拾取或者通用拾取器第二步用ISOLA、TDMT_INV、FOCI之类的工具。看似清晰实际上衔接很痛苦。首先是格式转换各个工具要求的输入格式不一样其次是拾取结果到反演环节的传递完全脱节拾取出来的到时和极性没有办法直观地反馈到反演结果上。更深层的问题是矩张量反演对震相拾取的精度要求比定位高得多。定位允许0.2~0.3秒的到时误差反演振幅时哪怕0.1秒的偏差截取的波形段就不同振幅拟合残差立马变大。手动拾取的主观性在这里体现得特别明显同一个人隔天捡同一组波形P波到时差个一两帧太正常了。所以把拾取和反演放进同一个框架里不光是省事更重要的是让拾取结果可以直接作为反演的质控输入逻辑上形成闭环。1.2 为什么用Matlab实现这个工具选择Matlab我认为是基于实际的。波形处理本质上就是矩阵运算Matlab的向量化操作写滤波、滑动窗、特征值分解都非常顺手。第二是生态信号处理工具箱、优化工具箱、绘图功能齐全震相拾取里用到的带通滤波、AIC计算、极化分析所需的协方差矩阵特征分解都是几行代码的事。第三是交互Matlab的图形窗口非常适合做拾取结果的人工校对画波形、画理论到时、标注极性一个figure就能搞定。当然Matlab的缺点是循环慢但事件级的数据量其实很小单事件几十个台站、每个台站几千个采样点向量化之后毫秒级就算完了完全不是瓶颈。2. 相位拾取的三级结构STA/LTA、AIC与极化分析2.1 STA/LTA先粗筛STA/LTA是最经典的自动触发算法原理就是用短时窗平均值和长时窗平均值的比值找突变。短窗反映瞬时能量变化长窗反映背景噪声水平当地震波到达时能量突然增大这个比值会冲高超过阈值就认为有震相到达。实际实现时特征函数可以选择原始振幅绝对值或者能量。我推荐用能量对高频震相的响应更灵敏。参数上STA窗长0.5~1秒LTA窗长5~10秒触发阈值2.5~4。阈值设置就看你的信噪比条件高信噪比的天然地震可以把阈值打到4低信噪比的微震得压到2.2左右但阈值低了误触发率也会上升所以STA/LTA只能做粗筛不能单独用它定到时。2.2 AIC精确定时STA/LTA告诉你“这里有信号”AIC告诉你“信号精确从哪个采样点开始”。AIC的原理是把一段波形在某个分界点k处分成前后两段计算两段数据各自拟合自回归模型的AIC值之和当k正好在震相到时时前后两段的统计特征差异最大AIC取得最小值。具体计算时不需要在整个数据窗上搜索全局AIC最小值那样容易把到时定在远处的尾波大振幅上。正确做法是先用STA/LTA得到一个粗略到时然后在这个到时前后各取一段范围比如到时前0.5秒到后1秒在这个局部窗口内计算AIC最小值位置就是精确到时。这样既快又稳。2.3 极化分析分离P波和S波P波和S波的重要区别在于质点运动方向P波质点沿射线路径振动近似径向且偏垂直S波质点垂直射线路径振动偏水平。利用三分量记录做极化分析可以自动给震相分类。实现方法是取事件窗内Z、N、E三分量数据组成N×3的矩阵计算3×3协方差矩阵再做特征分解。最大特征值对应的特征向量就是质点的极化主轴。如果主轴方向接近台站到震源的射线方向判定为P波如果主轴方向接近垂直于射线方向判定为S波。同时计算极化度也就是最大特征值占总能量的比例极化度大于0.6说明质点运动线性度好这个震相的质量可靠可以用于反演。这里有个技巧极化分析的时间窗不要开太大取0.1~0.2秒就够太长了容易把后续的转换波、反射波混进来反而把主轴的指向搞乱。3. 矩张量反演的关键参数与稳定性控制3.1 从观测方程到最小二乘矩张量反演的物理基础是远场位移可以表示为矩张量分量与格林函数空间导数的线性组合。反演时我们把每个台站观测到的P波和S波位移振幅写成向量d把由速度模型、射线路径、辐射花样决定的系数写成矩阵G矩张量的6个独立分量写成向量m于是问题简化为线性方程组d G·m。理论上的解就是最小二乘解。实际处理中还有个常见约束纯剪切破裂的矩张量迹为零也就是ISO分量为0。如果震源有体积变化比如流体注入导致的膨胀就得放开这个约束同时解ISO、CLVD和DC三个部分。这里我给个建议先把ISO约束为零跑一遍如果拟合残差明显大再放开约束。反过来一上来就放飞六分量容易反演出物理解释不了的震源机制。3.2 台站几何与阻尼做矩张量反演最怕的不是噪声是台站分布不好。台站如果全部集中在一个方位角范围反演矩阵就是病态的解对噪声极其敏感。我自己实测方位角覆盖小于180度时条件数轻松破千反演出来的沙滩球几乎每天都不一样换个滤波频段就变样。所以实操层面有两条硬性要求一是参与反演的台站最好不少于6个并且方位角尽量均匀展开二是反演前检查G矩阵的条件数条件数超过1000就考虑加阻尼正则化把解变成m(GᵀGλI)⁻¹Gᵀd。λ的取值可以通过L曲线法确定日常处理我给个大致范围0.01到1之间信噪比越低λ适当取大。还有一个容易被忽略的点近台站振幅大、信噪比高但离震源越近格林函数对速度模型误差越敏感。我一般会给近台适当降权防止一两个台主导整个反演结果。3.3 结果解读与质量检验反演得到6个矩张量分量后需要分解成物理意义明确的三个部分双力偶DC、补偿线性矢量偶极CLVD和各向同性ISO。DC分量代表剪切破裂是地壳地震的主要机制ISO代表体积变化常见于火山活动或流体注入CLVD则常与复杂裂隙或非双力偶源有关。怎么判断反演结果可信我习惯看三个指标一是波形或振幅拟合残差拟合相关系数至少0.8以上二是走时残差拾取的到时和理论到时差控制在0.1秒内三是分解后DC占比是否落在合理区间。如果反演出一个ISO高达60%但同时又带大量DC的结果先别急着写论文大概率是速度模型不对或者某个台站的极性标错了。4. 实操流程从波形到沙滩球4.1 数据准备与预处理建议每个事件单独建一个目录waveforms、metadata、results三个子目录分开。波形先统一转成SAC或者miniSEED格式Matlab下用现成的文件读取函数导入。第一步做去均值、去线性趋势、去仪器响应。去仪器响应这一步很多新手会跳过去实际上它对振幅反演是致命的——不同台站的仪器响应不一致反演出来的振幅比就是错的。接下来是带通滤波。微震数据我常用2~8Hz的带通如果事件很小、高频丰富可以提到5~15Hz。滤波之后检查所有台站的采样率是否统一不统一先重采样。数据是速度记录的话反演前要积分转成位移Matlab里用cumtrapz即可注意频域积分的话要除以2πf。4.2 拾取参数配置拾取参数建议写在一个配置脚本里方便批量跑事件。我的常用配置是这样fs 200; % 采样率按实际修改 stalta_len_short 0.5*fs; % STA窗长 0.5s stalta_len_long 10*fs; % LTA窗长 10s thresh_on 3.0; % 触发阈值 aic_range [0.5 1.5]*fs; % AIC搜索范围到时前0.5s到后1.5s pol_window 0.1*fs; % 极化分析窗长如果波形噪声大把触发阈值降到2.2同时把AIC搜索范围适当拉大不然容易把到时的候选段弄丢。跑完自动拾取之后我的习惯是用Matlab的图形窗口把所有台站的波形画出来把自动拾取的到时标注在图上快速扫一眼。这一步大概花两三分钟但能避免后面反演出结果后一脸茫然。4.3 反演参数配置反演需要速度模型、震源位置、参与反演的台站列表以及时窗长度。速度模型至少要有三层沉积层、结晶地壳、上地幔格式类似深度(km) Vp(km/s) Vs(km/s) 密度(g/cm³) 0.0 5.8 3.4 2.6 5.0 6.2 3.6 2.7 20.0 6.8 3.9 2.9时窗长度上P波段取到时前0.05秒到后0.45秒S波段取到时前0.1秒到后0.9秒。窗口开太长会把后续震相接进来开太短又截不全有效振幅。反演带宽要和拾取滤波保持一致否则振幅关系对不上。4.4 一次完整运行的流程记录完整跑一次事件大概是这样的顺序读取config文件载入台站坐标、速度模型、事件初步定位结果。对每个台站做预处理、带通滤波、仪器响应校正。STA/LTA粗拾取得到可能的震相区间。AIC精确定时得到P波或S波的精确到时。极化分析对震相分类并同时输出质点极化方向和极化度。汇总所有台站的到时、振幅、极性写入phase文件。计算理论格林函数系数组装G矩阵。最小二乘反演做DC/CLVD/ISO分解。输出沙滩球图、拟合残差、走时残差。整个过程如果自动跑单事件几分钟内搞定加上人工校对拾取结果一般也就五分钟到十分钟一个事件比手工流程快了一个量级。5. 常见问题与排查技巧5.1 拾取问题和反演问题速查问题现象可能原因排查与解决P波到时被拾到S波上信噪比低STA/LTA误触发提高触发阈值或增加极化分析二次判别AIC到时总偏向早窗口内混入前一个事件的尾波缩小AIC搜索范围检查事件间隔反演的ISO分量异常高速度模型过于粗略或台站方位角覆盖差更新速度模型增加台站或固定ISO0重跑拟合残差集中在个别台站时窗包含了转换震相或反射震相缩小反演时窗避开已知转换波到达时间反演结果对滤波频段特别敏感台站几何差导致解不稳定检查G矩阵条件数加阻尼正则化沙滩球和已知构造应力场明显矛盾某个台站通道极性标反了用已知事件的P波初动方向标定极性表5.2 我给新手的三个独家小技巧第一个技巧S波不要在全波形上找。很多工具默认在全波形上搜索S波误触率非常高。我的做法是先用P波到时截一段压制P波能量后在P波到时的后续窗口里做短时能量比值扫描这样找到的S波到时可靠得多。第二个技巧反演前花五分钟做极性自检。把所有台站的垂直分量初始运动方向画出来和已知的震源机制、射线路径对比一下如果整体反了一个方向说明台站极性定义有问题。这个问题在微震监测台网中尤其常见更换地震计、接线松动都可能造成极性翻转不做自检往往会得出一个完美但错误的沙滩球。第三个技巧别盲目追求低残差。反演拟合残差低不意味着结果就是对的过度拟合噪声的低残差结果反而危险。我更倾向于在台站覆盖良好、相位质量明确的数据子集上做反演哪怕少用几个台站也比把所有台站一股脑丢进去要稳。最后再分享一点体会用这套流程大半年下来我最深的感受是拾取质量直接决定反演上限。很多人花大量时间调反演参数却忽略了相位拾取的精度实际上0.1秒的到时误差就能让振幅反演结果明显偏掉。pickmt把拾取和反演放进一个框架最大价值不是给你省掉手动操作而是让拾取结果可控可查出了异常能一眼定位到是拾取的问题还是反演的问题。后续可以扩展的方向也不少比如接入深度学习拾取器做初筛、在台阵数据上做相对矩张量反演、根据反演残差迭代更新速度模型。如果手头正好有微震或者近震数据建议先用高信噪比事件的完整流程跑通再慢慢处理低信噪比的难题这个工具会越用越顺手。本文还有配套的精品资源点击获取