先说说我为什么会对这份数据上心。做地震学研究的人都知道震相表是绕不开的基础数据之一。大到地震定位、走时层析成像小到一次课程设计里的震相到时拾取都要和“某个台站在某时某刻记录到了某个震相”这种记录打交道。但现实是想拿到一份覆盖中国区域、格式统一、字段齐全的震相表并不容易台网目录、论文附录、机构内部资料各有各的格式连震相命名都可能不统一。CnOpenData 的中国地震震相表正好把这些麻烦事压缩到了一个数据集里。我第一次拿到这份表的时候第一反应是“干净”第二反应是“终于不用自己从十几个网页里扒数据了”。这篇文章不搞虚的直接把我用这份数据的完整思路、字段理解、实操代码和踩过的坑都摆出来想研究地震、做数据分析、或者写毕业论文的人都可以按着这份经验走一遍。这份数据适合谁简单说两类人最受益一是做地球物理、地震学研究的学生和科研人员需要用到真实震相到时做定位或者走时模拟二是搞数据科学和机器学习的开发者想拿真实地震数据训练震相自动拾取模型但又不想花时间处理原始事件波形数据。对这两类人来说CnOpenData 中国地震震相表提供的价值不是“一张表”那么简单它把观测端最核心的震相信息结构化、标准化了省掉了大量清理和校验工作。后面我会一步步拆解这张表里到底有什么坑、怎么用才能发挥它的价值。1. 内容整体设计与思路拆解1.1 震相表到底是一张什么表先建立共识。地震震相表通俗讲就是一张“地震到时的体检记录单”。每次地震发生后布设在各地的地震台站会记录到地震波到达的时间最常见的两大震相是 P 波纵波跑得快和 S 波横波跑得慢。震相表就是把这些记录逐条列出来哪个地震哪个台站在什么时刻读到了什么震相以及这个台站相对震中大概多远。有了这张表你就能做很多事。最经典的用法是地震定位近震用 P、S 到时差算震中距远震用多个台站的 P 波走时做交切再结合多个台站就能反推出地震发生的位置和深度。没有震相表这些工作就得自己从波形文件里手动量取到时效率极低而且不同人量取的标准还会有偏差。CnOpenData 整理的中国地震震相表本质上就是把这一步前置了它把散落的观测记录打包成标准结构化表格所以我一直把这类数据称作“研究的原材料”。1.2 CnOpenData 这份数据为什么值得用我自己用过的震相数据来源很多包括一些国际数据库和国内台网的重做目录但多数时候需要自己拼接、转换时间格式、统一震相命名。CnOpenData 的中国地震震相表有几个很实际的优势。首先覆盖范围聚焦中国区域。这一点对做区域构造、地壳速度结构研究的用户很关键数据集中在国内台网记录空间关联更容易建立。其次字段设计明显是奔着“直接可分析”去的除了震相到时还包括事件信息、台站信息、震级、距离、方位角等不用再到处关联外部表。最后从文件格式和编码处理来看它给的是规整表格直接读进 pandas 或者 ArcGIS、GMT 都能接得住。对我这种习惯用 Python 处理数据的人这已经是很省心的结构了。这里必须多提一句不同批次的表字段名可能漂移比如时间列可能叫time也可能叫arrival_time。我下面写的都是基于我拿到的某个快照版本大家使用前先df.info()看一眼再动手不要死记列名。1.3 怎么判断这份数据靠不靠谱任何数据集拿回来第一件事不是跑模型而是验数据质量。震相表的核心是“到时”和“震相名称”我通常用两个土办法验证。第一个办法抽取一个已知的大地震事件比如某个 6 级以上地震把震相到时减去发震时间得到走时再根据台站距离做一条理论走时曲线看规律。P 波走时应随距离稳定增加S 波走时同样随距离增加二者之间的差也随距离增加。如果画出来的点乱七八糟说明表里可能混入了错检的震相。第二个办法检查震相类型的可用性。近震台站主要记 P、S远震台站可能记 Pn、 Sn、 Pdiff 等。如果一个“地震震相表”里全是同一种震相且距离跨度很小那说明它可能只提取了某个单一目录的子集。CnOpenData 这份表我看下来震相种类还算丰富能覆盖到从近震到远震的常见情况。当然数据采集过程中的拾取误差一定存在所以后续分析必须加入质量控制步骤。这也是下面我要重点展开的实操内容。2. 核心细节解析与实操要点2.1 主要字段结构与解读拿到表之后别急着画图先把列名一个个拆开看。以我手里的版本为例典型的列结构可以分为四类事件信息、台站信息、震相信息、计算派生信息。事件信息包括事件编号、发震时刻、震中纬度、震中经度、震源深度、震级、震级类型。这些字段告诉你是哪一次地震。发震时刻通常精确到秒甚至毫秒震中经纬度一般用十进制度。震源深度单位一般给的是公里但有些表会写成米需要认真看元数据。台站信息包括台网代码、台站代码、台站纬度、台站经度、台站高程。台站代码是全局唯一的例如同一个台站名在不同台网下可能是两个台站因此台网代码和台站代码需要一起用。震相信息是这张表的核心包括震相名称和震相到时。震相名称一般是标准 IASPEI 命名比如 P、S、Pn、Sn、PmP、SKS 等。震相到时是台站实际记录到的波到达时刻单位通常是 UTC 时间格式可能是2023-01-01 12:34:56.78这种字符串。派生字段则包括震中距可能有度和公里两种单位、方位角、慢度等。这类字段通常是基于某个参考地球速度模型计算出来的比如用走时表算出理论到时偏差。看到这些列时先确认它的计算单位是弧度、度还是公里不然画图时很容易错得离谱。下表是我整理过的一个字段对照实际使用时可以直接当速查手册。字段类别典型字段名单位/格式说明事件信息event_id字符串地震事件唯一编号事件信息origin_timeUTC 时间发震时刻事件信息ev_lat, ev_lon十进制度震中经纬度事件信息depthkm震源深度事件信息mag无震级数值台站信息net字符串台网代码台站信息sta字符串台站代码台站信息st_lat, st_lon十进制度台站经纬度震相信息phase字符串震相名称震相信息phase_timeUTC 时间震相到时派生信息dist_deg度震中距派生信息azim度台站相对震中的方位角派生信息travel_time秒走时phase_time - origin_time2.2 震相命名规则先读懂名字再分析很多初学者拿到震相表看到一列phase第一反应是“应该就是 P 和 S 吧”结果一数里面居然有十几种名字一下就懵了。其实震相命名有一套逻辑掌握了就是肌肉记忆。P 和 S 是基础体波。P 波是纵波S 波是横波。在近震范围内直达波最先到的是 Pg 和 Sg这里的“g”表示地壳内传播的直达波。Pn 和 Sn 表示莫霍面折射波可以理解为地壳底部的首波。PmP 是地壳底部反射波在近震震相多样性分析里经常见到。远震情况下P 波穿过地幔还会分出 PKP、SKS、Pdiff 等核幔边界相关震相。我处理表里数据的一般原则是近震研究只保留 Pg、Sg、Pn、Sn如果表里没细分到 Pg/Sg直接把 P 和 S 作为近似远震研究则可以放宽到 P、S、Pdiff、PKP 等。换句话说不要看到非 P/S 就觉得是脏数据先确认自己的研究尺度再决定保留哪些震相。注意如果某个事件的震相表里同时出现 P 和 Pn它们不是重复记录。P 表示直达 P 波Pn 表示浅层折射首波两者走时特征不同混用会导致定位偏差。2.3 数据精度和单位的坑震相表最容易被忽视的是精度问题。发震时刻和震相到时的精度直接影响走时计算如果精度只到秒定位误差可能会到几公里甚至十几公里。理论上现代地震目录能到 0.01 秒甚至毫秒级别但人工修编和自动拾取混合的数据集里依然能看到不少只精确到秒的记录。做精细定位前最好先看看时间字段的长度如果都是秒级建议在方法上接受相应误差不要强求真秒级定位。单位问题同样要命。震源深度大部分表给 km但也有部分表给 m如果直接拿去做走时计算整个深度会差一千倍。震中距也一样有的表给的是度有的表给公里。处理时如果你需要统一可以按地球半径把公里转成度dist_deg dist_km / 111.19纬度方向上一度大约 111.19 公里但严格说需要按球面距离公式计算。另一个高频坑是震级类型同一次地震mb、Ms和Mw数值差异可能很大筛选大事件时不要只看mag列还要看mag_type。2.4 数据下载与文件格式说明CnOpenData 的数据下载后一般是压缩包解压后常见的是 CSV 或 Excel 格式。CSV 文件需要注意编码中文标注如果有乱码用 UTF-8-SIG 读取基本能解决import pandas as pd df pd.read_csv(CnOpenData_China_Seismic_Phase.csv, encodingutf-8-sig) print(df.shape) print(df.columns.tolist())如果文件较大建议在读取时就指定需要的列和 dtype避免把大字段全读进内存。比如只需要事件、台站和震相字段可以用usecols参数df pd.read_csv( CnOpenData_China_Seismic_Phase.csv, usecols[event_id, origin_time, phase, phase_time, sta, sta_lat, sta_lon, mag, depth], parse_dates[origin_time, phase_time], encodingutf-8-sig )这也是我踩过的一个小坑一开始直接全字段读取一个接近千万行的表把内存吃掉了大半后来精简字段才顺利跑动。数据量大的时候能用多少列就只读多少列这是处理表格数据的基本素养。3. 实操过程与核心环节实现3.1 数据清洗与基础筛选拿到原始表以后我先做一轮系统性清洗。第一步是检查空值和重复值特别是核心字段event_id、phase、phase_time不能有空。第二步是时间字段统一origin_time和phase_time转成 pandas 的datetime64[ns]类型并去除明显不合逻辑的记录比如到时早于发震时刻的记录。下面是我常用的一段清洗流程import pandas as pd import numpy as np df pd.read_csv(CnOpenData_China_Seismic_Phase.csv, encodingutf-8-sig) df[origin_time] pd.to_datetime(df[origin_time], errorscoerce) df[phase_time] pd.to_datetime(df[phase_time], errorscoerce) # 去掉核心列缺失的记录 df df.dropna(subset[event_id, phase, phase_time]) # 去掉到时晚于发震时刻太多或早于发震时刻的异常记录 df df[(df[phase_time] df[origin_time]) (df[phase_time] - df[origin_time] pd.Timedelta(hours3))] # 删除同一事件、同一台站、同一震相的重复记录 df df.drop_duplicates(subset[event_id, sta, phase, phase_time]) # 计算走时秒 df[travel_time] (df[phase_time] - df[origin_time]).dt.total_seconds()这段代码里最关键的是“到时不能早于发震”和“重复记录去重”。我在实际数据里见过因为格式错位导致到时比发震时刻早了几个月的记录这通常是不同系统的时间基准没对齐造成的不洗干净后面每一步都会受影响。筛选地震尺度时我一般按震级过滤# 只看 5 级以上地震方便走时规律观察 msel df[df[mag] 5.0] print(msel[event_id].nunique())3.2 用走时曲线快速评估整张表的质量走时曲线是验证震相表质量最直观的手段。原理很简单P 波速度比 S 波快所以同一台站 S 波到时会晚于 P 波并且两种震相的走时都随台站距离增大而递增。如果表里大量点在走时曲线图上“满天飞”要么是距离算错了要么是震相名称标错了。我通常按事件随机抽几个出来画图import matplotlib.pyplot as plt import pandas as pd sample_events df[event_id].drop_duplicates().sample(5, random_state42) fig, axes plt.subplots(2, 3, figsize(14, 8), sharexTrue, shareyTrue) axes axes.flatten() for i, ev in enumerate(sample_events): ev_data df[df[event_id] ev] ax axes[i] for ph, color in zip([P, Pn, S, Sn], [blue, cyan, red, orange]): sub ev_data[ev_data[phase] ph] if not sub.empty: ax.scatter(sub[dist_deg], sub[travel_time], s12, labelph, colorcolor) ax.set_title(fEvent {ev}) ax.legend(fontsize8) ax.grid(alpha0.3) for ax in axes.flat: ax.set_xlabel(Distance (deg)) ax.set_ylabel(Travel time (s)) plt.tight_layout() plt.show()正常数据画出来会呈现清晰的两条“带状曲线”上面是 S 波下面是 P 波。如果某个事件只有两三个台站记录在图上可能只有稀疏的两三个点这不代表数据有问题只是台站密度不够。如果点分布完全不合理就要考虑该事件的震相拾取是否混入了噪声。3.3 简单地震定位用到时差估算震中距震相表最核心的应用就是地震定位。这里我不准备讲复杂的非线性反演先演示一个新手也能上手的单事件双台定位思路。近震条件下假设地壳平均 P 波速度vp 6.0 km/sS 波速度vs 3.5 km/s即波速比vp/vs 1.71。同一台站记录的 S-P 到时差dt ts - tp震中距d可以近似为d dt / (1/vs - 1/vp)换算成常用写法就是vp 6.0 vs 3.5 dt ts - tp # 单位为秒 dist_km dt / (1/vs - 1/vp)在代码里我先选一个事件找到同时记录到 P 和 S 的台站计算每个台站的 S-P 到时差再换算距离# 选取一个事件 ev_id sample_events.iloc[0] ev df[(df[event_id] ev_id) (df[phase].isin([P, S]))] # 透视出 P 和 S 到时 pivot ev.pivot_table(indexsta, columnsphase, valuesphase_time, aggfuncfirst) pivot[dt] (pivot[S] - pivot[P]).dt.total_seconds() # 用平均波速估算震中距 vp, vs 6.0, 3.5 pivot[dist_km] pivot[dt] / (1/vs - 1/vp) pivot[dist_deg] pivot[dist_km] / 111.19 print(pivot[[dt, dist_km, dist_deg]].head())注意这种简易方法假设波速均匀适合岩石圈尺度比较小的区域。实际定位中要结合更多台站和更复杂的速度模型但作为理解震相表物理含义的入门练习已经够用。3.4 结合台站坐标做几何可视化有了台站经纬度和估算的震中距你其实可以在地图上画圆以每个台站为圆心以估算距离为半径多个圆的交点就是大致震中位置。用 Python 的cartopy画一个示例import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt fig plt.figure(figsize(8, 8)) ax fig.add_subplot(1, 1, 1, projectionccrs.PlateCarree()) ax.set_extent([100, 110, 28, 38], crsccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle:) # 台站 震中距圆的示意 for _, row in pivot.iterrows(): ax.scatter(row[sta_lon], row[sta_lat], marker^, colorblack, s40) circle plt.Circle( (row[sta_lon], row[sta_lat]), row[dist_deg], colorblue, alpha0.15, transformccrs.PlateCarree() ) ax.add_patch(circle) plt.show()从图上你会直观看到震中距圆是否收敛在一定范围。如果圆交不到一起可能是速度模型给得不准也可能是某条震相被读成了错误类别。这就是震相表里的常见问题之一下一节我会集中聊这些坑。4. 常见问题与排查技巧实录4.1 时间格式和时区不统一震相表这种多来源数据最容易出问题的就是时间。我遇到过四种情况时间列是字符串没有转 datetime时间列带时区标记但部分行没有发震时刻用的 UTC震相到时用的本地时间还有个别行phase_time比发震时刻早明显是异常值。处理时我的建议是先统一成 UTC不做任何本地化转换。用 pandas 解析时加errorscoerce无效时间会转为NaT方便后面定位df[phase_time] pd.to_datetime(df[phase_time], utcTrue, errorscoerce)解析完再查一遍isna()的数量如果异常比例偏高说明原文件可能存在其他时间编码格式需要回到源头确认。4.2 震相列出现非标准编码自动拾取和人工修正混合的数据震相列偶尔会出现PA,SB,XP之类的非标准编码或者同一个实际震相被写成两种叫法。我的处理思路是先看不同取值的频数再决定是保留还是清洗。print(df[phase].value_counts())如果频数高且明确可辨识就做映射phase_map { P: P, p: P, Pg: P, Pb: P, S: S, s: S, Sg: S, Sb: S, Pn: Pn, Sn: Sn, PmP: PmP, SmS: SmS } df[phase_clean] df[phase].map(phase_map).fillna(OTHER)把分析范围以外的震相统一标记成OTHER这样既保住了数据的完整性也避免它们在定位里捣乱。4.3 台站坐标缺失或错位震相表通常会直接给出台站经纬度但不同批次可能存在台站名一样、坐标不一致的问题。比如同一个代码在不同时期经过迁移位置从 A 点搬到了 B 点。这种坑很隐蔽表面看不出来一旦你画震中距圆就会露馅。我会用两个办法校验一是和公开台站台账做关联二是检查同一台站经纬度的标准差。如果某个台站的经纬度离散度明显偏大就重点检查是不是用了不同时期的位置。另一类问题是经纬度列名混有sta_lat与sta_lat大小写之类的情况读取后注意统一列名。4.4 走时差为负或异常偏大这是震相表里最影响结果的一类数据质量陷阱。走时差为负说明到时早于发震时刻逻辑上不可能除非时钟不同步或震相拾取错误。走时异常偏大比如距离不到 100 公里却用了 200 秒的走时通常说明该记录是远震误标成了近震。遇到这类记录我建议根据研究尺度做阈值过滤。比如只研究区域震时可以直接排除走时大于 500 秒的记录如果研究远震则反过来排除走时过短的记录。不要试图修复这类异常值删除或标记为异常即可。4.5 常见问题速查表症状可能原因处理方案时间列读出来是 object 字符串未加 parse_dates 参数用pd.to_datetime(..., errorscoerce)部分到时早于发震时间时区或时钟偏差直接过滤保留逻辑上成立的行同一台站同一震相出现多条数据源合并重复按 event_id sta phase phase_time 去重phase 列有很多奇怪编码自动拾取程序输出非标准名映射到标准震相无法映射则归为 OTHERdist_deg 与 dist_km 数值对不上单位理解错误以元数据为准统一换算成目标单位台站坐标画图时出现飞点台站搬迁或用错经纬度对照公开台站台账校正4.6 实操心得永远保留一份原始数据副本这是我在处理各种震相表时最深刻的经验。无论做了多么精细的清洗原始数据一定要留一个只读副本。因为清洗逻辑一旦出错你可能需要回到最初状态重新处理。我用一个很土但有效的流程下载后先复制一份raw_origin.csv所有清洗操作都在新的 DataFrame 上做并且每一次转换都写成独立的代码段。这样分析结果出问题的时候能倒查是清洗逻辑的问题还是原始表本身的问题。强烈建议任何做数据方向的人养好这个习惯。5. 延伸应用与影响范围分析5.1 成为震相自动拾取模型的训练集现在地震波形数据越来越多靠人工一个台站一个台站读震相已经不现实很多组把自动拾取模型作为研究方向。CnOpenData 中国地震震相表最大的潜在价值之一就是可以作为训练集的标注来源。具体做法是把震相表中的事件、台站、到时作为标签再对应到台站的波形数据切成 P 窗口和 S 窗口喂给深度学习模型。表里字段越标准构造训练数据集越省力。不过我提醒一句模型训练前必须做严格质量控制因为自动拾取器的训练数据如果本身就包含错标震相模型学到的就是错的规律。可以先按波速比一致性做一轮筛选比如vp/vs在 1.65~1.80 范围外的记录要重点复核。5.2 走时层析成像的观测方程构建走时层析成像是地震学里的常规武器用大量 P、S 波走时反演地下速度结构。使用震相表时你需要把每个“震相到时-发震时刻”转成走时再把台站位置和震源位置写成观测方程。表里如果已经给了走时或慢度列这一步会更快但反演前依然要做射线路径检查避免因为表里的无效记录污染反演结果。5.3 与其它数据产品联合挖掘震相表单独用能做定位但和地震目录、台站元数据、波形数据放在一起能做的事就完全不同了。比如提取天然地震的背景噪声、约束衰减参数、分析断层带非均质性等等。这类数据交叉分析的关键在于主键设计和时间对齐。我一般会用event_id关联目录用net.sta关联台站元数据尽量满足三方数据无缝衔接。5.4 扩展思路从区域震相表走向自主分析管道拿到一张震相表只是第一步。我后来做过的一个项目是在这份数据基础上搭了一条半自动管道输入一个区域、一个震级范围自动筛选事件自动计算走时残差再输出异常事件清单。整套逻辑其实不复杂核心就是把上文提到的清洗、校验、画图、定位这几步串起来写成脚本。这样做的好处是什么第一每次更新数据后不用人工重复检查第二可以形成统一的 QC 报告哪些事件台站覆盖好、哪些事件明显有问题一眼就能看到。我已经把这套经验用在了几个数据集上效果稳定。如果读者们有兴趣后续我也可以详细写写这条管道里的模块设计细节。回看我自己用 CnOpenData 中国地震震相表的整个流程最大的体会就是“结构化数据 严格质控”带来的效率提升远远超过预期。很多人拿到表后第一反应是赶紧跑代码出图但实际上先花半小时验验数据、清清洗洗、理理字段后面会省出几个小时的排查时间。另外一个小技巧把它和强震目录、连续波形数据放在一起管理时尽量用同一个事件编号体系哪怕需要自己写映射表也绝对值得。数据不是拿来“跑一下”的而是要在稳定的数据基础上反复打磨才会让后续每个分析都更踏实。