第一次拿到这套中国地震震相表的时候我确实原地愣了好一会儿。文件名看起来挺直白压缩包解出来以后里面躺着几十万行带时间、带台站代码、带震相名称的表格。但真要把这些字段变成能用于研究、能写出论文、能支撑工程分析的数据中间隔着的距离比想象中大得多。震相表这玩意儿表面上是哪个台站在哪个时刻记到了哪种波实际上背后是一个区域的地壳速度结构、台站布局质量、地震事件定位精度甚至还有当年记录仪器和人工分析员的水平痕迹。这篇文章我想从一个实际使用者的角度把拿到这套数据之后应该怎么做、应该想什么、容易在哪翻车尽量完整地梳理一遍。1. 打开震相表之前先搞清楚这几个基础概念1.1 什么是震相什么人会需要这张表地震发生时震源会向四面八方辐射弹性波。这些波穿过不同介质、经过不同路径到达地表台站在波形记录上表现为一个个起跳点。每个起跳点对应一种特定类型的波就叫一个震相seismic phase。专业分析员每天的工作很大一部分就是在地震图上把这一个个起跳点标出来记录下精确到十分之一秒的到时再把它们汇编成表。这套中国地震震相表就是把散落在各省台网、各个历史时期报告里的震相到时记录做成了统一格式的表格数据。它适合谁用三类人最需要。第一类是做区域地震学研究的需要用真实震相数据研究地壳结构、断层活动第二类是搞地震监测和预警算法开发的需要真实样本验证定位精度第三类是数据科学领域的人想拿一套有物理含义的真实时空数据做可视化、做特征工程。无论哪一类前提都是先把数据本身读明白。1.2 常见震相名称与物理含义速查震相名称不是随便起的每个字母都有讲究。第一个字母代表波的类型P是纵波S是横波。下标或后缀代表传播路径。整理一个我在实际处理中经常遇到的速查表震相代码含义典型特征Pg地壳内传播的直达P波近震中最清晰、最先到达Pn沿莫霍面滑行的折射P波震中距较大时先于Pg到达Sg地壳内传播的直达S波P波之后约数秒至数十秒到达Sn沿莫霍面滑行的折射S波与Pn成对出现P穿过地幔的体波P波远震记录中的主震相S穿过地幔的体波S波远震中P之后的大幅度震相L / R面波Love波/Rayleigh波速度最慢、振幅最大、周期长为什么搞清这些基础概念很重要因为直接拿P和S去匹配数据遇到Pn、Pg就会全乱了。我见过有同事把Pn到时当P波拿去算震中距结果定位偏离真实震中几十公里。类似这种问题不是数据本身错了而是使用者对震相体系的物理意义理解不够。2. 字段逐项拆解与第一次数据体检2.1 每个字段是什么、怎么用把震相表打开之后一般会看到这样几类核心字段事件标识、发震时刻、震中坐标、震级、台站代码、震相名称、到时、震中距、方位角、振幅、周期等。我挑几个最容易用错或者最容易忽略的详细说。事件标识这一列看着人畜无害多套数据合并时却最容易出错。不同来源的事件编号规则不一样有的按时间编码有的是数据库自动增长的ID还有一个事件在不同台网里可能有多个编号。我拿到数据后做的第一件事永远是检查同一个事件标识下是否有重复行、是否有两个事件标识指向同一条记录。这类脏数据在早期手写汇编转录入的年代特别常见。震中距和方位角这两列很多人在分析时会直接用它但一定要弄清楚单位。震中距有的给角度度有的给公里有的按地球平均半径做了不同椭球修正。方位角则是从台站看震中的方向角跟从震中看台站是反的用错方向会直接影响后续反演。振幅和周期字段则要格外小心。同一个事件不同台站记录的振幅能否直接拿来对比取决于仪器响应是否做过校正。如果原始记录里混用了不同时期、不同增益的仪器数据这些字段只能谨慎使用甚至更适合直接丢弃。2.2 时间格式与坐标口径——两个最隐蔽的坑处理震相表最阴险的问题不是数据内容本身而是时间格式和坐标口径不统一。时间这一列有些记录写的是北京时间有些写的是UTC甚至同一张表里不同时期的数据混用。哪怕只有八小时的偏差配合震相走时分析就是致命的P波从震中到几百公里外台站也就几十秒八小时的系统性偏差会让所有走时曲线彻底失去意义。坐标也是这样。震中经纬度有的用地理坐标有的用地心坐标台站坐标有的用WGS84有的用北京54或者西安80等历史坐标系统。在地壳尺度的研究中几百米的坐标系统差异可能不算大但在地震定位和震相到时残差分析里一两公里的差异就能造成明显的系统偏差。检查方法也很朴素把表里所有事件的经纬度画到地图上如果看到某些事件点整体偏移到海里去了或者与其他事件之间存在奇怪的系统性错位十有八九是坐标口径出了问题。我自己的习惯是拿到表格先做一轮极简体检统计行数、检查缺失值比例、画出经纬度散点图、查看时间范围是否合理、检查每个事件平均关联台站数量。这套流程大概半小时能筛掉一半以上的低级问题。3. 震相表为什么值钱从记录到地球内部结构3.1 走时曲线的原理一张图看出地下结构震相表最有价值的地方是它可以用于构建走时曲线。把同一个地震事件在所有台站的P波到时减去发震时刻得到P波的走时再把每个台站的震中距作为横轴、走时作为纵轴把这些点连起来就得到一条走时曲线。地下介质的速度结构会直接反映在走时曲线的斜率上。如果地下是均匀介质走时曲线就是一条过原点的直线斜率是速度的倒数。真实地球不是均匀的P波在地壳浅层慢、在深部快穿过莫霍面进入地幔之后更快所以走时曲线会出现明显的拐点和分段特征。Pn震相在图中表现为一条斜率较缓的直线段因为它以较高的速度沿莫霍面滑行。只看一张走时曲线图你就能大致判断这个区域的地壳平均速度和莫霍面深度。这也让走时曲线天然成为数据质量检查的利器。把理论走时曲线叠在实际数据上偏离过大的点基本都有问题要么震相标错要么定位误差太大要么台站钟有偏差。3.2 用震相数据做速度研究的三种典型思路有了可靠的震相表可以做的事情远不止画图围观。第一种思路是走时反演。把已知的地震事件、已知的台站位置和已知的震相到时组合起来反演出地壳分层速度结构。这是最经典的做法Pn走时反演可以得到莫霍面之上的平均P波速度S-P到时差结合P波走时可以推出Vp/Vs比值进而估算泊松比、判断岩性。第二种思路是震源参数精化。通过多台站的震相到时做联合定位可以把单个地震事件的位置和发震时刻整体优化。尤其是那些记录年代较早、定位精度较低的事件用现代方法重新处理一遍往往能得到明显更好的结果。第三种思路是做地震活动性统计。震相表里的事件分布、震级频度关系b值、时间序列特征反映着一个地区的应力积累状态和活动水平。这类研究对单个事件的绝对精度要求不那么高但对数据完整性要求很高处理时要特别注意台网检测能力随时间的变化。从震相记录到地下结构认识本质上是一个反演问题。正演是从模型算记录反演是从记录推模型。震相表就是整个链条中最基础、最原始的一环它值钱的原因在于原始波形数据量太大、格式不统一而震相表已经完成了从波形到特征数字的关键一步压缩。4. 清洗与质检把表格变成可用数据集的完整流程4.1 第一步格式统一与缺失值处理实际操作中我拿到手的数据先要过三道格式关卡。第一关是时间格式。不管原始表哪列是文本、哪列是数值全部统一成ISO格式的UTC时间字符串再转成时间戳存储。这一步看似简单但有一种情况特别折磨人记录里只写了日期和时间没有标注时区而当时的历史报告里默认写的北京时间。处理这种数据只能结合发震时刻和台站位置做间接验证没有绝对可靠的办法所以我建议在清洗文档里标注时间基准为推测可能存在系统性偏差。第二关是震相名称的归一化。不同历史时期对同一震相的叫法不完全一样有的写成Pn有的写成P n有的写成小写pn。还有更麻烦的情况把初至震相统一标成P但实际可能是Pg这种信息丢失是无法恢复的只能按现有标记使用。第三关是缺失值处理。震相表最常缺的是振幅、周期这类辅助参数其次是某个台站对某个事件没有记录。我一般的处理原则是缺失辅助字段的行可以保留后续分析按列忽略缺失核心字段震相到时、台站代码的行直接删除极少数行连震相名称都是空的这类行大概率是录入错误没有补全的意义。4.2 第二步用理论走时残差做自动化质检表清洗干净之后下一步就是把坏数据挑出来。最有效的方法不是人眼看而是做理论走时残差检查。原理是这样的每一个地震事件都有位置和发震时刻每一个台站都有坐标。用一个已知的地球速度模型比如通用的全球一维模型iasp91或者某个区域的粗略分层模型可以正演出理论上某个震相应该在哪一秒到达某个台站。把理论到时减去表中实际到时得到的就是走时残差。正常情况下Pn震相的区域台网走时残差应该在一两秒以内Pg震相也在几秒量级。如果你看到一个标记为P的震相残差达到几十秒甚至上百秒原因通常只有一个这个震相被标错了。最常见的错法是S波被标成P波或者某个后续震相被当成了初至。我处理过一次数据发现一个台网早期记录的P波残差系统性偏大10秒以上后来比对了相邻台站才发现那批记录实际标的是P之后第一个强振幅震相而不是初至P。具体的做法可以写成简单的脚本逻辑按事件遍历所有记录对每条记录用该事件的位置、发震时刻和台站坐标计算理论到时然后与实际到时作差。对残差绝对值超过阈值的记录标记为可疑。阈值的选取要看区域地壳薄、台站密集的地区阈值可以收紧到2秒高原和山区路径上速度结构复杂阈值放宽到5秒更合理。4.3 第三步可视化复核——走时图一眼看出问题自动化检查能筛掉大量明显的错误但总有些看起来合理的错误混在合格区间里。这时我会把所有干净数据绘成走时图用眼睛再扫一遍。画法也简单横轴是震中距纵轴是P波和S波的走时用不同颜色区分震相。一张合格的地壳尺度走时图上你能看到清晰的两簇点P簇从原点出发以较陡斜率向外延伸S簇斜率较缓地跟在后面。如果某个色点孤零零地落在两条主带之间的空白区域它大概率是个错标震相。如果某一片区域的点整体偏移那往往是定位系统误差或者速度模型不适用造成的。我第一次画图的时候图上一堆S波的点散落在P波带里密密麻麻根本分不出两簇。查了很久才发现是数据里的震相字母大小写混用了S被写成了小写s在分组统计时跟大写S没分到一起加上部分台站把Pg标成了p。统一格式之后再画两簇就清清楚楚了。吃一堑长一智可视化不只是给文章配图用的它本身就是最好的数据质量检查工具。5. 三种最常见的应用落地方式5.1 走时曲线绘制与模型对比先说最简单也最常用的一种用清洗后的数据绘制走时曲线并跟理论模型对比。实现思路很直接。给每个事件算出每个震相的走的时和震中距然后散点图把P、S、Pn、Pg分别画出来。画图用的库不用纠结Python的matplotlib就够。再调用某个走时计算库比如用pyrocko或者obspy里的taup接口拿到iasp91的理论走时曲线叠在图上对比。我常用的一段流程大致是这样import pandas as pd import matplotlib.pyplot as plt # 读取已清洗好的震相表 df pd.read_csv(clean_phase_table.csv, parse_dates[phase_time]) # 走时计算phase_time - origin_time df[travel_time] (df[phase_time] - df[origin_time]).dt.total_seconds() # 按震相分组画散点 fig, ax plt.subplots(figsize(8, 6)) for phase_name, color in [(P, #4444cc), (S, #cc4444)]: sub df[df[phase_name] phase_name] ax.scatter(sub[epi_distance_km], sub[travel_time], s8, labelphase_name, colorcolor, alpha0.6) ax.set_xlabel(震中距 (km)) ax.set_ylabel(走时 (s)) ax.legend() ax.grid(alpha0.3) plt.savefig(travel_time_check.png, dpi150)画完图之后对比区域性的走时特征跟全球模型的差异本身就能说明很多问题。比如某一区域的实测Pn斜率和iasp91差异偏大说明那一片上地幔顶部的速度跟全球平均不一样。差距越大说明这个区域的地壳地幔结构越特殊这往往是深入研究选题的起点。5.2 单事件定位从P-S到时差到台网交汇另一个很有教学价值的应用是用震相表实现手动版地震定位。虽然专业定位软件一秒钟就能算完但亲手做一遍能极大加深对定位原理的理解。基本原理不复杂。S波和P波在同一路径上的速度不同S波后到二者到时差S-P与震中距正相关。已知区域平均Vp和Vs时震中距可以近似写成[ \Delta (t_S - t_P) \cdot \frac{Vp \cdot Vs}{Vp - Vs} ]一个台站只能告诉你震源在距离它\Delta公里的某个圆上。至少三个台站的圆交汇才能把震源位置约束到一个区域内。实际操作中我把事件关联的所有台站都找出来按上述公式给每个台站画一个圆然后找出最密集的交点区域。这个方法的误差来源很直观Vp/Vs的估计偏离真实值圆半径就偏台站分布在震源一侧时圆与圆的交角很小定位精度急剧下降。亲手体会过这些误差之后你就会明白现代定位算法里为什么要引入残差加权和阻尼最小二乘——无非是想办法让不同台站的信息合理组合并压制坏数据的影响。如果想把震中距算得更准确就不该用简单的平均速度而是要结合走时曲线。在走时图上反查P波、S波走时对应的震中距相当于用实测走时曲线做了更精细的模板匹配。当年用人工方式做定位本质上就是在做这件事。5.3 进阶Pn速度反演与双差定位的探索方向如果你对反演有兴趣震相表还能更进一步。Pn速度反演是相对成熟的做法。Pn震相沿莫霍面滑行其走时可以写成震中距除以Pn速度再加上一个与莫霍面深度、上地幔速度相关的截距项。把大量Pn走时数据放在一起通过拟合走时-距离开的直线的斜率和截距就能估计Pn速度和莫霍面性质。做这类分析时震相表的台站覆盖、方位角分布、震相识别的稳定性都直接影响结果。双差定位double-difference relocation则是把同一地区邻近事件与同一台站的走时差进行比较消掉路径上的公共误差从而获得高精度相对位置。它对震相到时质量要求很高但对绝对理论模型的依赖较弱。用震相表里同一簇小地震的数据做双差定位能显著提升震群空间分布的清晰度用来识别断层精细结构非常合适。这些进阶分析不一定每个人都需要做但了解它们能帮你判断手里的震相表还有多大的潜力没挖掘。6. 实操中真正值得记住的几条经验6.1 台网与仪器响应对震相数据的影响震相表数据看起来是客观记录但实际上它带着当时台网水平和分析员处理习惯的烙印。早期台站密度低记录仪灵敏度有限大量小震根本没被记录到震相表里小震事件存在明显的空间不均匀性。台站少的地方一个事件可能只有一个台站有数据这种事件定位精度很差用来做统计分析时要格外小心。不同台站的仪器响应不同也会影响震相辨认。短周期仪器对高频P波敏感长周期仪器更容易记录到面波宽频带仪器记录S波更清楚。如果你发现某个台站记录的S波数量系统性偏少不一定是台站附近地壳结构异常可能单纯是仪器频带覆盖不够。所以在分析中我建议记录每个事件的有效台站数分析时用这一列做加权或者做筛选效果会好很多。一个事件只有一两个台站记录时它只能告诉你这里发生过地震任何更精细的解释都不可靠。6.2 版本管理与元数据记录震相表这类数据产品会不定期更新。旧版本里缺失的记录新版本可能补上了早期目录中的定位结果新版可能用更优算法重新处理过。如果你的论文或项目用到了这套数据强烈建议记录下载日期和数据版本号并把清洗过程中的中间文件都留一份方便追溯。我自己在本地处理时的习惯是所有原始文件放进一个readonly目录清洗脚本输出到一个processing目录最终用于分析的干净数据放在analysis目录。每一步产生的文件都带时间戳。这样即使三个月后重新整理分析也能一眼看出当时的处理链路。还有一个小细节如果你在论文里使用了这批数据引用时除了注明数据版本还要详细说明震相挑选准则和清洗步骤。原因很简单别人要复现你的研究最不确定的环节往往不是算法本身而是数据处理的那部分主观决定。6.3 不要把质量等级当免检金牌有些数据产品会在每条记录上给出质量等级比如A、B、C或者好、中、差。很多人处理数据时习惯性地只保留高等级记录把低等级记录全部过滤掉。这个做法在样本量足够大时是合理的但对于小样本研究就要慎重。质量等级往往是根据震相清晰度和信噪比判定的它不完全等于物理可靠性。B级记录可能只是噪声稍大但到时标定并不偏差C级记录中可能包含某些特殊路径下的清晰震相只是波形形态不那么标准而被打低分。我把过一遍全区走时数据后最大的体会是任何等级字段都只是先验信息真正可靠的验证还是要靠走时残差一致性检查和可视化复核来把关。6.4 对完整保持怀疑震相表从汇编完成那一刻起就已经比真实世界丢失了一部分信息。没有哪个台网能记录到所有地震没有哪个分析员能识别所有震相。拿到一张几十万行的震相表第一反应不应该是数据真多而应该是这里头缺了什么。台网的检测能力随时间和空间变化——某年之前仪器老旧某年之后增加了新台站——这不叫数据缺陷而是数据的观测系统属性。做区域地震活动性对比时要先把检测能力曲线估算出来否则被观测能力迷惑会把台网升级误判为地震活动增强。这套中国地震震相表我前后用了将近两年最大的收获反而不在某个具体结论上而是养成了处理任何地理时空数据时先问观测系统长什么样的习惯。再把话说回来数据质量永远是相对的把使用目标和数据特点对齐比找到完美数据更现实。如果你正准备拿这套数据做研究不妨从画第一张走时图开始——一张图不够就画十张数据里的问题会自己浮到表面来的。