CnOpenData中国地震震相表解析：从数据清洗到地震定位

先说说我为什么会对这份数据上心。做地震学研究的人都知道，震相表是绕不开的基础数据之一。大到地震定位、走时层析成像，小到一次课程设计里的震相到时拾取，都要和“某个台站在某时某刻记录到了某个震相”这种记录打交道。但现实是，想拿到一份覆盖中国区域、格式统一、字段齐全的震相表并不容易，台网目录、论文附录、机构内部资料各有各的格式，连震相命名都可能不统一。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这种字符串。

派生字段则包括震中距（可能有度和公里两种单位）、方位角、慢度等。这类字段通常是基于某个参考地球速度模型计算出来的，比如用走时表算出理论到时偏差。看到这些列时，先确认它的计算单位是弧度、度还是公里，不然画图时很容易错得离谱。

下表是我整理过的一个字段对照，实际使用时可以直接当速查手册。

2.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', encoding='utf-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'], encoding='utf-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', encoding='utf-8-sig') df['origin_time'] = pd.to_datetime(df['origin_time'], errors='coerce') df['phase_time'] = pd.to_datetime(df['phase_time'], errors='coerce') # 去掉核心列缺失的记录 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(hours=3))] # 删除同一事件、同一台站、同一震相的重复记录 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_state=42) fig, axes = plt.subplots(2, 3, figsize=(14, 8), sharex=True, sharey=True) 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'], s=12, label=ph, color=color) ax.set_title(f"Event {ev}") ax.legend(fontsize=8) ax.grid(alpha=0.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/s，S 波速度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(index='sta', columns='phase', values='phase_time', aggfunc='first') 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, projection=ccrs.PlateCarree()) ax.set_extent([100, 110, 28, 38], crs=ccrs.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='^', color='black', s=40) circle = plt.Circle( (row['sta_lon'], row['sta_lat']), row['dist_deg'], color='blue', alpha=0.15, transform=ccrs.PlateCarree() ) ax.add_patch(circle) plt.show()

从图上你会直观看到震中距圆是否收敛在一定范围。如果圆交不到一起，可能是速度模型给得不准，也可能是某条震相被读成了错误类别。这就是震相表里的常见问题之一，下一节我会集中聊这些坑。

4. 常见问题与排查技巧实录

4.1 时间格式和时区不统一

震相表这种多来源数据，最容易出问题的就是时间。我遇到过四种情况：时间列是字符串没有转 datetime；时间列带时区标记，但部分行没有；发震时刻用的 UTC，震相到时用的本地时间；还有个别行phase_time比发震时刻早，明显是异常值。

处理时我的建议是，先统一成 UTC，不做任何本地化转换。用 pandas 解析时，加errors='coerce'，无效时间会转为NaT，方便后面定位：

df['phase_time'] = pd.to_datetime(df['phase_time'], utc=True, errors='coerce')

解析完再查一遍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 常见问题速查表

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 中国地震震相表的整个流程，最大的体会就是“结构化数据 + 严格质控”带来的效率提升远远超过预期。很多人拿到表后第一反应是赶紧跑代码出图，但实际上，先花半小时验验数据、清清洗洗、理理字段，后面会省出几个小时的排查时间。另外一个小技巧：把它和强震目录、连续波形数据放在一起管理时，尽量用同一个事件编号体系，哪怕需要自己写映射表，也绝对值得。数据不是拿来“跑一下”的，而是要在稳定的数据基础上反复打磨，才会让后续每个分析都更踏实。