简介pickmt 是一套基于 MATLAB 开发的地震矩张量分析工具面向具备一定 MATLAB 编程基础与地震学背景的研究者用于解决震源机制分析中相位拾取与参数反演的问题。其核心功能分为两部分相位拾取负责从地震波形中识别并提取精确的到时信息反演则借助数学优化方法求解震源矩张量等震源参数可应用于地震事件定位与震源特性研究。资源包共 31 个文件以 26 个 .m 脚本为主体另含 .mat 数据、.txt 说明与 README.md 文档压缩包约 30.57MB目录涵盖合成数据生成、滤波绘图、相位拾取与反演等模块结构清晰便于按流程查阅。目前已有 729 人学习下载。借助 MATLAB 的数值计算与可视化能力读者可获得一套可运行、可修改的震源分析脚本并依据示例数据与说明文档快速上手按自身研究需求扩展功能。1. 从一段微震波形说起pickmt 到底在拾什么、反什么如果你手上有三分量微震记录想从一段看不出 P 波起跳的波形里把震源机制抠出来pickmt 这类 Matlab 代码就是干这个的。名字拆开看pick 是相位拾取mt 是 moment tensor矩张量合起来就是「先拾相位、再反演矩张量」的一条流水线。它解决的不是「有没有事件」而是「这个事件的破裂面朝哪、是剪切还是张裂、震级多大」。适合做矿山微震、水力压裂监测、室内声发射、岩石力学实验的人——只要你拿得到至少 6 个通道的三分量波形并且大致知道台站坐标和速度模型就能跑。反直觉的地方在于很多人以为矩张量反演靠的是振幅其实相位拾取的精度才是那条最短的木板拾错 2 个采样点机制解可能整个翻面。2. 相位拾取从 STA/LTA 到 AICpickmt 里到底用哪套2.1 为什么不能只靠 STA/LTA 交差短时平均比长时平均STA/LTA是最省事的拾取器但它给的是「能量突变点」不是「初动到达点」。在低信噪比微震里STA/LTA 经常提前半个周期触发因为 P 波前面的噪声被窗口平滑掉了。pickmt 这类代码通常把 STA/LTA 当粗定位再用 AIC赤池信息准则在粗定位附近做精拾取。AIC 的逻辑是把信号分成两段一段是噪声、一段是信号找让两段模型拟合误差之和最小的那个点。它对初动敏感但对窗口长度也敏感——窗口开太大AIC 最小值会漂到波形中段。我一般会这么配STA 窗口取 0.51 个主频周期LTA 窗口取 510 倍 STA触发阈值 35 倍均方差。粗触发后再用 AIC 在触发点前后各 2050 个采样点内搜最小值。下面是一段可直接抄的 Matlab 骨架% 输入data 为单通道列向量fs 为采样率 % 输出tp 为拾取到的采样点索引 function tp pick_aic(data, fs) % --- STA/LTA 粗触发 --- sta_len round(0.01 * fs); % STA 窗口 10 ms lta_len round(0.10 * fs); % LTA 窗口 100 ms thr 4.0; % 触发阈值 n length(data); ratio zeros(n,1); for i lta_lensta_len : n sta mean(data(i-sta_len1:i).^2); lta mean(data(i-lta_len-sta_len1:i-sta_len).^2); ratio(i) sta / (lta eps); end k find(ratio thr, 1, first); if isempty(k), tp []; return; end % --- AIC 精拾取在粗触发点前后各 30 点内搜索 --- w 30; seg data(max(1,k-w) : min(n,kw)); m length(seg); aic zeros(m,1); for j 2:m-1 v1 var(seg(1:j)); v2 var(seg(j1:end)); aic(j) j * log(v1eps) (m-j-1) * log(v2eps); end [~, jmin] min(aic(2:end-1)); tp max(1,k-w) jmin; % 换算回全局索引 end逻辑说明STA/LTA 负责「大概在哪」AIC 负责「精确到点」。参数上sta_len和lta_len是核心主频越高窗口越短thr调低会多触发、调高会漏触发微震建议从 4 开始试。w是搜索半径太大容易把 AIC 最小值引到 S 波或尾波上太小则可能框不住真实初动。跑完一定要画图核对别信数字。2.2 三分量联合拾取别让单通道骗了你单通道拾取在 Z 分量上做遇到 P 波初动弱、S 波强的记录很容易把 S 波当成 P 波。pickmt 的常见做法是三分量联合先算三个分量的能量包络再算偏振参数比如协方差矩阵的特征值比用「线性度」和「入射角」辅助判断。P 波到达时质点运动接近线性、方向沿射线S 波到达时线性度下降、横向分量能量上升。实操里我会把 Z 分量拾取结果作为候选再用水平两分量的偏振线性度做二次确认线性度低于 0.6 的候选直接丢掉。% 三分量偏振线性度检查 % E,N,Z 为等长列向量fs 采样率tp 为候选拾取点 function ok check_polarization(E, N, Z, fs, tp) w round(0.02 * fs); % 20 ms 窗口 idx max(1,tp-w) : min(length(Z), tpw); X [E(idx), N(idx), Z(idx)]; X X - mean(X); % 去均值 C cov(X); % 3x3 协方差矩阵 lambda sort(eig(C), descend); linearity 1 - lambda(2)/(lambda(1)eps); % 线性度 ok linearity 0.6; % 阈值可调 end参数说明窗口w取 12 个主频周期太短协方差估计不稳太长会把 S 波混进来。linearity阈值 0.6 是经验值噪声大时降到 0.5信号干净时提到 0.7。这一步不产生新拾取只做「否决」能挡掉相当一部分误拾。2.3 拾取结果怎么存、怎么传给反演拾取完不要直接塞进反演先存成结构体或表格字段至少包含台站编号、通道、P 拾取采样点、S 拾取采样点如果有、拾取质量标记。质量标记建议分三档A 表示人工核对过B 表示自动拾取且偏振检查通过C 表示自动拾取但偏振存疑。反演时只喂 A 和 BC 档当备选。这个习惯能省掉后面大量返工——我见过太多人把 C 档数据混进去反演出来的机制解自己都不敢信。3. 矩张量反演从观测到六个分量的最小二乘3.1 矩张量到底是什么为什么是六个分量矩张量 M 是一个 3×3 对称矩阵九个元素里六个独立。它描述的是震源处等效体力偶物理上对应「破裂面滑动 体积变化」。反演的目标就是解出这六个分量再分解成 ISO各向同性、DC双力偶、CLVD补偿线性矢量偶极三部分。对微震来说DC 占比高说明是剪切破裂ISO 占比高可能是张裂或塌陷。pickmt 的反演核心是用拾取到的 P/S 到时和振幅构建观测向量 d再用格林函数矩阵 G 建立 d G·m最后最小二乘解 m。3.2 格林函数怎么来速度模型和射线追踪格林函数是反演的「黑匣子」——它把震源到台站的传播效应算出来。常见做法是一维层状速度模型 射线追踪算出每个台站每个相位的走时和振幅衰减。速度模型不用很精细但分层不能太离谱层数太少走时误差大层数太多反演不稳定。我一般用 35 层P 波速度从 3.0 到 6.0 km/s 递增S 波速度按 Vp/Vs 1.7 换算。如果研究区有测井或以往定位结果优先用实测值别硬套通用模型。% 一维层状模型射线追踪简化版仅算直达 P 走时 % vp 为各层 P 波速度向量h 为各层厚度向量src 和 rec 为 [x y z] function [tt, amp] trace_1d(vp, h, src, rec) dz abs(src(3) - rec(3)); dh sqrt((src(1)-rec(1))^2 (src(2)-rec(2))^2); dist sqrt(dz^2 dh^2); % 简化用平均速度算走时 v_avg sum(vp .* h) / sum(h); tt dist / v_avg; amp 1 / (dist eps); % 几何扩散 end逻辑说明这是最简版本只算直达波走时和几何扩散。真实反演里还要考虑层间折射、自由表面放大、仪器响应。参数上vp和h必须和实际模型一致src和rec坐标单位统一都用米或都用千米。如果走时残差普遍大于 5 ms先回头查速度模型别急着调反演。3.3 最小二乘求解与约束观测方程 d G·m 里G 是台站数×6 的矩阵d 是观测到时/振幅向量。台站少于 6 个时方程欠定必须加约束。常见约束有两种一是矩张量迹为零假设无体积变化二是固定 ISO 分量比例。pickmt 一般用阻尼最小二乘m (GᵀG λI)⁻¹Gᵀdλ 是阻尼因子控制解的稳定性。λ 太大解被压扁太小解震荡。实操里从 0.1 开始试看残差和范数的 trade-off 曲线拐点附近取值。% 阻尼最小二乘反演矩张量 % G: n x 6 格林函数矩阵, d: n x 1 观测向量, lambda: 阻尼因子 function m invert_mt(G, d, lambda) [n, ~] size(G); if n 6 warning(台站数少于 6解不唯一建议加约束); end I eye(6); m (G * G lambda * I) \ (G * d); % 输出残差供检查 res d - G * m; fprintf(残差范数: %.4f\n, norm(res)); end参数说明lambda是唯一需要反复调的参数建议扫 0.01、0.1、1、10 四个量级画残差-范数曲线选拐点。G的列顺序要和m的分量顺序一致通常 Mxx, Myy, Mzz, Mxy, Mxz, Myz别搞混。反演完一定要做分解看 DC/ISO/CLVD 比例比例离谱说明数据或模型有问题。4. 避坑与排查拾取和反演里最容易翻车的五件事4.1 现象拾取点整体偏移反演残差普遍偏大原因仪器响应没去除或者采样率换算错了。很多微震记录仪有内置滤波直接拿原始计数做拾取初动会被延迟。解决先做去仪器响应和去均值确认时间轴单位是秒不是采样点。检查方法很简单——拿一个已知到时的事件看拾取误差是否在 12 个采样点内。4.2 现象反演出的机制解每次都不一样换个台站组合就翻面原因台站方位角覆盖太差所有台站挤在一个象限里。矩张量反演对方位角覆盖极敏感覆盖差时解不唯一。解决至少保证台站分布在三个以上象限如果做不到加约束比如固定 DC 比例并明确标注解的不确定性。别硬报一个「精确解」。4.3 现象AIC 拾取在低信噪比下跳到 S 波上原因搜索窗口w开太大或者 STA/LTA 粗触发点本身就偏后。解决缩小w到 1020 个采样点同时把 STA/LTA 阈值降到 3 左右让粗触发提前。如果还是跳改用三分量偏振线性度做二次否决线性度低的候选直接弃用。4.4 现象反演残差很小但机制解物理上不合理原因过拟合。阻尼因子lambda太小解在拟合噪声。解决增大lambda看残差是否只是轻微上升而解变得稳定。另一个可能是格林函数矩阵条件数太大检查台站坐标是否有重复或过近的。4.5 现象代码跑通但结果和以往定位对不上原因坐标系不一致。有的用 ENU东-北-天有的用 NED北-东-下混用会让矩张量分量错位。解决在代码开头统一声明坐标系所有输入输出都按同一套来。我一般用 ENU因为 Matlab 里画图直观。转换时注意 Z 分量符号别把向下当向上。5. 进阶技巧用合成测试验证整条流水线反演代码写完别急着上真实数据。先做合成测试自己造一个已知矩张量用正演算出各台站的理论波形加噪声再走一遍拾取和反演看能不能把已知解找回来。这一步是「后悔药」——真实数据里你永远不知道真值合成数据里你知道。具体做法取一个纯 DC 源Mxy 1其余为 0用 3.2 的射线追踪算走时用 Ricker 子波合成波形加 10% 高斯噪声然后跑拾取和反演。% 合成测试已知 DC 源验证反演能否恢复 M_true [0 1 0; 1 0 0; 0 0 0]; % 纯 DCMxy1 % 正演对每个台站算理论到时和振幅合成 Ricker 子波 % ...正演代码略核心是 G * m_true 得到 d_syn d_syn G * [M_true(1,1); M_true(2,2); M_true(3,3); ... M_true(1,2); M_true(1,3); M_true(2,3)]; d_noisy d_syn 0.1 * randn(size(d_syn)); % 10% 噪声 m_est invert_mt(G, d_noisy, 0.1); % 对比 m_est 和真值看 DC 占比是否接近 1参数说明噪声水平从 5% 开始加逐步到 20%看反演在哪个噪声水平下开始失稳。这个失稳阈值就是你真实数据的可信边界。如果 10% 噪声下解就翻面说明台站布局或速度模型有问题先解决再上真实数据。合成测试还有一个好处能标定拾取误差对反演的影响。故意把拾取点偏移 1、2、3 个采样点看机制解怎么变心里就有数了。我自己的习惯是每换一个研究区先跑合成测试再跑一个已知事件最后才批量处理。这套流程救过我很多次——有一次合成测试里 15% 噪声下解就崩了回头查发现是台站坐标里有一个点录错了改过来之后 25% 噪声下都稳。别嫌麻烦拾取和反演这条链上任何一个环节的疏忽都会在最终解里放大。希望帮到你。本文还有配套的精品资源点击获取