EMR方法实战：从台网目录到监测能力曲线的最小算例

简介这份资源面向地震学研究者、地震台网运维人员及相关专业学生聚焦利用EMR经验震级关系方法估算地震台网的最小完整性震级Mc为台网监测能力评估与布局优化提供可复用的计算工具。压缩包共13个文件全部为m脚本整体约17KB涵盖核心Mc计算、KS检验与bootstrapping误差评估、震级-频度分布拟合、正态分布与泊松概率辅助函数以及示例运行脚本构成一套从数据统计到模型验证的完整流程。已有196人学习下载说明其在同类工具中具有一定参考价值。读者可借助这些脚本输入自有地震目录完成震级与观测参数经验关系的建立、Mc阈值推算及不确定性评估并据此判断台网对低震级事件的探测能力适合作为地震监测能力分析与台站优化研究的基础代码框架。1. 从一张台网震级底图说起EMR 方法到底在算什么如果你手头有一份区域地震台网的观测报告想回答“这套台网在本地到底能监测到多小的地震”最直接的办法不是翻规范而是把 EMR 方法跑一遍。EMR 是“震级—距离”经验关系Empirical Magnitude–Range的缩写核心思路很朴素用台网历史记录中每个台站实际能检测到的最小震级按震中距分档统计再拟合出一条随距离衰减的监测能力曲线。它不依赖理论噪声模型也不要求你知道每个台基的场地响应只吃“已经记录到的目录”这一份数据。对做地震台网运维、台址勘选、监测能力评估的人来说这套方法最大的价值是结论直接对应“这套台网现在能报出几级地震”而不是停留在信噪比公式里。标题里的“EMR-example”就是一份可复现的最小算例把台网目录、震级、震中距三样东西串起来输出一张监测能力底图。适合谁适合手里有台网目录、想快速评估监测下限、又不想从零推导检测概率的从业者。2. EMR 方法的输入准备目录、震级与震中距怎么对齐2.1 为什么 EMR 只认“已记录到”的震相EMR 的统计对象是台网目录里每一条“某台站对某次地震有震相拾取”的记录。它不关心这次地震有没有被定位只关心这个台站有没有在这个震中距上留下可用的震级。常见做法是从台网观测报告中导出“台站—地震”配对表字段至少包含台站代码、震中距、震级、震相类型。这里有个容易翻车的点很多人直接拿地震目录的震级去配对但目录震级是台网平均震级不是单台震级。EMR 要的是单台震级否则距离分档后震级会被平均掉曲线整体偏低。我一般会从震相报告里取单台震级如果只有目录震级就退而求其次用“台站参与定位”的记录但要在文档里注明这个近似。2.2 震中距分档别用等间距用等对数震中距从几公里到几百公里跨度很大等间距分档会让近距档样本爆炸、远距档样本稀疏。常见做法是按对数等间距分档比如 0–10 km、10–20 km、20–50 km、50–100 km、100–200 km、200–500 km。每档内统计该档所有记录的最小震级或者取第 5 百分位震级作为“可检测下限”。为什么不用最小值因为单条异常记录会把下限拉低第 5 百分位更稳。下面这段 Python 就是做分档和下限统计的最小骨架。import numpy as np import pandas as pd # 读取台站-地震配对表字段station, dist_km, mag, phase df pd.read_csv(emr_pairs.csv) # 对数等间距分档边界 bins [0, 10, 20, 50, 100, 200, 500] labels [0-10, 10-20, 20-50, 50-100, 100-200, 200-500] df[dist_bin] pd.cut(df[dist_km], binsbins, labelslabels, rightFalse) # 每档取第5百分位震级作为检测下限 result df.groupby(dist_bin)[mag].quantile(0.05).reset_index() result.columns [dist_bin, mag_lower] print(result)逻辑说明pd.cut按给定边界分档rightFalse表示左闭右开避免边界值重复计入。quantile(0.05)取第 5 百分位比最小值抗异常。参数说明bins要根据你台网的实际震中距分布调整如果台网最远记录只有 300 km最后一档写到 500 没意义改成 300。phase字段可以用来过滤只保留 Pg/Sg 这类定位震相避免面波震级混进来。2.3 单台震级缺失时怎么补如果震相报告里只有振幅和周期没有直接给单台震级就得自己算。常见做法是用地方震震级公式M log(A) 1.11*log(R) 0.00189*R - 2.09其中 A 是最大振幅微米R 是震中距公里。这个公式是华南地区的经验式换区域要换系数。补完之后一定要做一致性检查把算出来的单台震级和目录震级做散点如果系统偏差超过 0.3 级说明公式系数不匹配得换。这一步没有捷径血泪经验是宁可少用几条记录也不要用错公式把整条曲线带偏。3. 用 EMR-example 跑通监测能力曲线拟合、绘图与参数3.1 最小可复现流程从配对表到能力曲线拿到分档下限后下一步是拟合一条“震级随距离变化”的曲线。EMR 常用线性拟合M(R) a * log10(R) b其中 R 是震中距a 和 b 是拟合系数。为什么用 log10(R)因为震级本身是对数量纲距离取对数后线性关系更稳。下面这段代码完成拟合并输出系数。import numpy as np from scipy.optimize import curve_fit # 取分档中心距离 bin_centers np.array([5, 15, 35, 75, 150, 350]) mag_lower result[mag_lower].values def emr_func(r, a, b): return a * np.log10(r) b popt, pcov curve_fit(emr_func, bin_centers, mag_lower) a, b popt print(f拟合系数: a{a:.3f}, b{b:.3f}) # 预测 50 km 处的监测下限 print(f50 km 监测下限: {emr_func(50, a, b):.2f})逻辑说明curve_fit用最小二乘拟合popt是系数pcov是协方差。参数说明bin_centers要和你实际分档一致如果分档是 0–10、10–20中心取 5 和 15 没问题但如果分档是对数等间距中心应该取几何平均而不是算术平均。拟合完看pcov对角线如果某个系数方差特别大说明该档样本太少考虑合并相邻档。3.2 绘图时最容易忽略的三个参数画监测能力图时横轴是震中距对数纵轴是震级。三个参数决定图能不能用第一横轴范围要覆盖台网实际监测范围别默认 0–500 km如果台网最远 200 km画到 500 就是空白第二纵轴下限要低于最小震级 0.5 级否则曲线贴底看不清第三分档点要画出来不能只画拟合线否则审图的人不知道你用了多少数据。下面这段是绘图骨架。import matplotlib.pyplot as plt plt.figure(figsize(8, 5)) plt.scatter(bin_centers, mag_lower, colorred, label分档下限) r_smooth np.logspace(np.log10(5), np.log10(350), 100) plt.plot(r_smooth, emr_func(r_smooth, a, b), b-, labelEMR拟合) plt.xscale(log) plt.xlabel(震中距 (km)) plt.ylabel(震级) plt.ylim(mag_lower.min() - 0.5, mag_lower.max() 0.5) plt.legend() plt.grid(True, whichboth, ls--) plt.savefig(emr_capability.png, dpi300)逻辑说明logspace生成对数等间距的平滑横轴避免折线感。whichboth让主次网格都显示方便读对数坐标。参数说明dpi300是投稿或报告用的分辨率如果只是内部看150 就够。ylim下限减 0.5 是为了让最低点不贴边。3.3 拟合系数的物理含义与边界a 值反映监测能力随距离衰减的快慢a 越大远距离能监测到的震级越高说明台网在远距越吃力。b 值是 1 km 处的截距但实际没有 1 km 的记录所以 b 只作为拟合参数不要单独解释。常见范围区域台网 a 在 1.0–2.0 之间b 在 -1.0–0.5 之间。如果 a 超过 2.5检查是不是远距档样本太少导致拟合被拉偏如果 a 小于 0.5检查是不是近距档震级普遍偏高把曲线压平了。这些边界值不是规范是我跑了几十个台网总结出来的经验区间超出就要回头查数据。4. 避坑与排查EMR 算例里最容易翻车的五件事4.1 现象曲线整体比预期低 0.5 级 → 原因用了目录震级而非单台震级 → 解决回震相报告取单台震级这是最常见的翻车。目录震级是多个台站平均后的结果天然比单台震级平滑且偏低。如果你直接拿目录震级配对分档下限会被拉低曲线整体下移。解决方法是回到震相报告取每个台站自己的震级。如果震相报告里没有单台震级用振幅周期自己算但要用对区域公式。4.2 现象远距档拟合线翘尾 → 原因远距档样本太少被个别大震级记录主导 → 解决合并档位或设最小样本数远距档往往只有几条记录如果这几条恰好都是较大地震下限就偏高拟合线在远端上翘。解决方法是设一个最小样本数比如每档至少 10 条记录不够就和相邻档合并。合并后档中心要重新算别直接用原来的。4.3 现象近距档下限异常低 → 原因混入了非定位震相或爆破记录 → 解决按震相类型和事件类型过滤近距档如果混入爆破、塌陷或非定位震相震级会异常低把下限拉下去。过滤方法是只保留 Pg/Sg 震相事件类型只保留天然地震爆破和塌陷单独处理。如果目录里没有事件类型字段用震源深度和波形特征辅助判断但这一步比较费时建议在数据准备阶段就做好。4.4 现象拟合不收敛 → 原因分档中心距离取了算术平均而非几何平均 → 解决对数分档用几何中心对数等间距分档时档内距离不是均匀分布算术平均会偏向大值。比如 10–20 km 档算术平均是 15但几何平均是 14.1。差别不大但如果所有档都偏拟合就会不收敛。解决方法是统一用几何平均sqrt(下限*上限)。4.5 现象换区域后曲线完全不可用 → 原因震级公式系数没换 → 解决按区域换地方震震级公式地方震震级公式是区域性的华南的系数拿到西北用系统偏差可能超过 0.5 级。换区域时先找该区域已发表的地方震震级公式如果没有用该区域台网的历史记录重新拟合系数。这一步没有后悔药只能老老实实做。5. 进阶技巧用 EMR 曲线反推台站增补优先级跑完 EMR 曲线只是第一步真正有用的是拿它做台站布局决策。具体做法把现有台网的 EMR 曲线画出来再模拟“如果某个方位增补一个台站该方位的监测下限能降多少”。常见做法是按方位角分扇区每个扇区单独跑 EMR找出监测能力最弱的扇区优先增补。下面这个表格是我常用的扇区评估模板填完就能看出哪个方位最需要补台。扇区记录数50 km 下限100 km 下限是否达标N1201.22.1是NE451.82.9否E881.42.3是SE322.03.2否填表时注意记录数少于 50 的扇区下限置信度低先补数据再下结论。达标标准按你的台网任务定如果任务是监测 1.5 级以上50 km 下限低于 1.5 才算达标。增补优先级按“下限差值 × 扇区面积”排序差值越大、面积越大越优先。最后说一个我自己的习惯每次跑完 EMR我都会把分档下限和拟合系数存成 CSV下次换区域时直接对比看是数据问题还是方法问题。这个习惯帮我省了很多重复排查的时间。希望帮到你。本文还有配套的精品资源点击获取