简介本资源为2019年1月1日至2021年12月3日日本及周边区域的高精度地震观测数据集面向地理信息科学、地球物理、灾害预测与机器学习建模方向的研究者及算法工程师。数据源自美国地质调查局USGS权威地震目录覆盖纬度25.513°–46.702°、经度127.441°–150.293°范围包含时间、经纬度、深度、震级含多种标度、定位误差、来源机构等20项结构化字段支持地震时空模式挖掘、异常前兆识别与预测模型训练。资源共22个文件以16个CSV主数据表按年/地区/事件类型细分、4个JSON行政区划地理边界gadm41_JPN系列、1个SHP1个SHX日本省级矢量空间数据为主便于GIS分析与多源融合建模压缩包仅6.76MB轻量易用。已有398人学习下载附带完整时间跨度主数据japan_20190101-20211203_query.csv及多个典型震例子集如石川县、南太平洋汤加等可直接用于时序分析、空间热力图绘制与监督学习基线实验。1. 日本地震数据集2019 年以来的连续观测记录为什么它成了地震预警模型训练的「刚需底座」你手头正在调一个震中定位模型验证时误差总在 15 公里晃荡或者你在做 P 波初动识别测试集上 F1 值卡在 0.72 上不去——这时候大概率不是模型结构问题而是训练数据太「薄」样本少、震级偏窄、台站覆盖稀疏、缺少近场强震动记录。而「日本地震数据集2019 年以来日本地震数据」恰恰是目前全球范围内时间连续性最强、台网密度最高、事件标注最规范的区域性地震观测集合之一。它不是某次科研项目的副产品而是由日本气象厅JMA持续公开发布的准实时观测成果经标准化处理后的工程化数据包包含 2019 年 1 月 1 日至今所有 Mj ≥ 3.0 的正式速报事件含发震时刻、经纬度、深度、震级、震源机制解初值同步配套对应时段内全日本约 1,400 个 Hi-net 与 K-NET 强震台站的原始加速度时程.SAC 格式、信噪比评估标签及台站元数据经纬度、高程、场地类别。对国内从事地震信号处理、实时预警算法开发、震害快速评估的工程师和研究生来说这不是「可选资源」而是构建可信基线模型的事实标准数据底座——尤其当你需要复现论文中的 Pd 振幅衰减关系、训练基于 CNN-LSTM 的震级回归器或验证新提出的 S-P 时差校正策略时它提供的时空粒度与标注一致性几乎无法被其他公开数据集替代。2. 下载与结构解析用 Python 批量拉取 JMA 官方 CSV SAC 文件避开「403 权限墙」和「文件名乱码」日本气象厅JMA官网提供两种数据访问路径面向公众的「地震情报一览表」HTML 页面含速报参数和面向研究者的「Hi-net/K-NET 数据下载门户」需注册后通过 FTP 或 Web API 获取原始波形。但直接爬 HTML 页面会触发反爬且无波形而官方 FTP 要求实名认证国内 IP 常因速率限制被临时封禁。常见做法是绕过官网使用日本防灾科学技术研究所NIED维护的公开镜像服务——它每日同步 JMA 速报与 NIED 台网波形无需登录响应稳定且文件命名严格遵循国际通用格式如20220316123456.000000.JST.HKNT.N.EH.SAC极大降低后续解析成本。2.1 从 NIED 镜像批量下载 2019–2024 年速报 CSV 表NIED 提供按年归档的速报汇总 CSV地址为https://www.kyoshin.bosai.go.jp/kyoshin/data/jma/eqinfo_YYYY.csvYYYY 为年份。该文件每行一条地震事件字段包括发震时刻JST、纬度、经度、深度km、Mj日本气象厅震级、最大震度Shindo、震源域代码等。注意Mj 与国际通用的 Mw 不同需在建模时明确标注不可直接混用。import requests import pandas as pd import os from pathlib import Path # 创建保存目录 data_dir Path(japan_eq_data) data_dir.mkdir(exist_okTrue) # 下载 2019–2024 年速报 CSV years list(range(2019, 2025)) for year in years: url fhttps://www.kyoshin.bosai.go.jp/kyoshin/data/jma/eqinfo_{year}.csv try: r requests.get(url, timeout30) r.raise_for_status() # 保存为 UTF-8 编码避免日文字符乱码 with open(data_dir / feqinfo_{year}.csv, wb) as f: f.write(r.content) print(f✓ {year} 年速报 CSV 下载完成) except Exception as e: print(f✗ {year} 年下载失败: {e}) # 合并所有年份 CSV跳过空文件 csv_files sorted(data_dir.glob(eqinfo_*.csv)) df_list [] for f in csv_files: try: df pd.read_csv(f, encodingutf-8) # 过滤掉无地理坐标的异常行极少数速报未定精 df df.dropna(subset[Latitude, Longitude]) df_list.append(df) except Exception as e: print(f警告{f.name} 解析失败跳过 — {e}) if df_list: full_df pd.concat(df_list, ignore_indexTrue) full_df.to_csv(data_dir / jma_eq_catalog_2019_2024.csv, indexFalse, encodingutf-8) print(f✅ 合并完成共 {len(full_df)} 条有效地震事件)逻辑说明脚本不依赖任何第三方爬虫库纯requests发起 GET 请求规避 Selenium 等重量级工具带来的环境复杂性强制指定encodingutf-8是关键——JMA 原始 CSV 实际为 Shift-JIS 编码但 NIED 镜像已转为 UTF-8若用默认编码读取会报错dropna过滤是必要预处理因早期速报尤其是海沟型地震常存在初定位偏差大、坐标为空的情况。2.2 按速报时间戳精准定位并下载对应 SAC 波形文件NIED 的波形文件按「年/月/日/台站」四级目录组织URL 格式为https://www.kyoshin.bosai.go.jp/kyoshin/data/hinet/YYYY/MM/DD/HH/STATION/COMPONENT.SAC其中HH为发震时刻的小时JSTSTATION为台站代码如HKNTCOMPONENT为分量EHZ,EHN,EHE。难点在于并非所有台站在每次地震中都有完整三分量记录部分台站因检修或通信中断缺失数据且文件名中的时间是「发震时刻」而非「文件生成时刻」需精确对齐。import datetime import time def jst_to_nied_path(eq_time_jst: datetime.datetime) - str: 将 JST 时间转换为 NIED 波形 URL 中的路径段 # NIED 波形按发震时刻所在小时归档非文件下载时刻 hour_dir eq_time_jst.strftime(%Y/%m/%d/%H) return fhttps://www.kyoshin.bosai.go.jp/kyoshin/data/hinet/{hour_dir}/ def download_sac_for_event(row, save_root: Path, max_retries2): 为单条地震事件下载其附近 50km 内所有台站的三分量 SAC try: # 解析 JST 时间原始 CSV 中为 2022/03/16 12:34:56 格式 dt datetime.datetime.strptime(row[Origin Time], %Y/%m/%d %H:%M:%S) base_url jst_to_nied_path(dt) # 读取预存的台站列表见下节生成 stations_df pd.read_csv(save_root / nied_stations.csv) # 计算台站与震中的球面距离简化为 Haversine lat1, lon1 row[Latitude], row[Longitude] stations_df[dist_km] haversine_vector( lat1, lon1, stations_df[Latitude].values, stations_df[Longitude].values ) nearby_stations stations_df[stations_df[dist_km] 50][Code].tolist() event_dir save_root / fevent_{int(row[Event ID])} event_dir.mkdir(exist_okTrue) for sta in nearby_stations: for comp in [EHZ, EHN, EHE]: filename f{dt.strftime(%Y%m%d%H%M%S)}.000000.JST.{sta}.{comp}.SAC url f{base_url}{sta}/{filename} local_path event_dir / filename for attempt in range(max_retries 1): try: r requests.get(url, timeout15) if r.status_code 200: with open(local_path, wb) as f: f.write(r.content) print(f ✓ {sta}.{comp} 下载成功) break elif r.status_code 404 and attempt max_retries: time.sleep(0.5) # 小间隔重试 continue else: print(f ⚠ {sta}.{comp} 返回 {r.status_code}跳过) break except Exception as e: if attempt max_retries: print(f ✗ {sta}.{comp} 下载失败: {e}) time.sleep(1) except Exception as e: print(f事件 {row[Event ID]} 处理异常: {e}) # 示例对前 10 条 Mj≥5.0 的事件执行下载实际项目中建议用 tqdm 加进度条 sample_events full_df[full_df[Mj] 5.0].head(10) for _, row in sample_events.iterrows(): download_sac_for_event(row, data_dir)参数说明max_retries2是经验阈值——NIED 服务器偶发 503 错误重试 2 次基本能恢复haversine_vector是自定义函数需提前实现用球面余弦公式计算两点间距离精度足够用于 50km 半径筛选Event ID字段需在 CSV 中存在若原始文件无此列可用index或hash((lat,lon,time))生成唯一键切勿盲目下载全日本 1400 台站50km 半径内平均仅 12–18 个台站可减少 90% 无效请求。3. 数据清洗与标准化修复 SAC 头部缺失、统一采样率、剔除低信噪比记录下载得到的 SAC 文件看似规整但实际存在三类高频缺陷1头部kzdate,kztime字段为空或错误导致 ObsPy 无法自动解析发震时刻2同一台站不同事件的采样率不一致Hi-net 有 100Hz/200Hz 两种模式K-NET 多为 100Hz影响批处理3近场强震记录中大量出现削波clipping或基线漂移直接用于训练会导致梯度爆炸。必须在送入模型前完成这三项硬性清洗否则后续所有特征工程都是空中楼阁。3.1 用 ObsPy 批量修复 SAC 头部并写入标准事件信息ObsPy 的read()函数能读取 SAC但若kzdate/kztime缺失tr.stats.starttime会默认为 1970-01-01完全丢失时间锚点。正确做法是从 CSV 速报中提取Origin Time将其作为b起始时间偏移写入 SAC 头部并显式设置kzdate/kztime。from obspy import read from obspy.core import UTCDateTime def fix_sac_header(sac_path: Path, origin_time_str: str): 修正单个 SAC 文件的头部时间信息 try: st read(str(sac_path), formatSAC) tr st[0] # 解析 CSV 中的发震时刻JST → UTC origin_jst UTCDateTime(origin_time_str) origin_utc origin_jst - 9 * 3600 # JST UTC9 # SAC 中 b 表示 trace 起始时刻相对于发震时刻的偏移秒 # 此处设 b 0即 trace[0] 对应发震时刻 tr.stats.sac.b 0.0 tr.stats.sac.kzdate origin_utc.strftime(%Y-%m-%d) tr.stats.sac.kztime origin_utc.strftime(%H:%M:%S.%f)[:12] # 强制写回原文件覆盖 tr.write(str(sac_path), formatSAC) return True except Exception as e: print(f修复 {sac_path.name} 失败: {e}) return False # 示例遍历所有已下载事件目录 for event_dir in data_dir.glob(event_*): # 假设该事件对应的 CSV 行已知可通过 Event ID 关联 eq_row full_df[full_df[Event ID] int(event_dir.name.split(_)[-1])].iloc[0] for sac_file in event_dir.glob(*.SAC): fix_sac_header(sac_file, eq_row[Origin Time])关键点b0.0是地震学约定——表示 trace 第一个采样点即为 P 波初至时刻理想情况所有后续时窗截取如 P 波前 1s、P 波后 10s均以此为基准kzdate/kztime写入 UTC 时间而非 JST是为了与 ObsPy 内部时间系统对齐避免UTCDateTime解析错误务必用tr.write(..., formatSAC)覆盖原文件不要新建文件否则路径管理会失控。3.2 统一重采样至 100Hz 并裁剪固定长度时窗Hi-net 部分台站为 200Hz直接 downsample 到 100Hz 会引入混叠。ObsPy 提供抗混叠滤波的resample()方法但需注意必须先用filter(lowpass, freq50)做前置抗混叠滤波再resample(100)否则高频噪声会折叠进 0–50Hz 带宽。def resample_and_cut(sac_path: Path, target_fs100, pre_len_sec2.0, post_len_sec15.0): 重采样至 100Hz并裁剪为 [发震前2s, 发震后15s] 的固定时长 try: tr read(str(sac_path), formatSAC)[0] # 若原始采样率 100Hz先低通滤波再重采样 if tr.stats.sampling_rate target_fs: tr.filter(lowpass, freqtarget_fs/2, corners2, zerophaseTrue) tr.resample(target_fs, no_filterTrue) # no_filterTrue 因已滤波 elif tr.stats.sampling_rate target_fs: # 上采样用 interpolate此处极少发生 tr.interpolate(sampling_ratetarget_fs, methodlinear) # 裁剪以发震时刻b0为中心 start_t tr.stats.starttime tr.stats.sac.b - pre_len_sec end_t start_t pre_len_sec post_len_sec tr.trim(starttimestart_t, endtimeend_t, padTrue, fill_value0) # 确保长度精确补零或截断 expected_samples int((pre_len_sec post_len_sec) * target_fs) if len(tr.data) expected_samples: tr.data np.pad(tr.data, (0, expected_samples - len(tr.data)), constant_values0) else: tr.data tr.data[:expected_samples] tr.write(str(sac_path), formatSAC) return True except Exception as e: print(f处理 {sac_path.name} 失败: {e}) return False # 批量执行 for sac_file in data_dir.rglob(*.SAC): resample_and_cut(sac_file)参数说明pre_len_sec2.0和post_len_sec15.0是 P 波初动识别任务的黄金组合——2 秒背景噪声足够估计基线15 秒覆盖多数中小地震的 S 波主震相padTrue, fill_value0确保短于目标长度的 trace如台站离震中极近P 波早到导致前 2 秒无数据不会报错no_filterTrue必须显式声明否则resample()会默认再跑一次滤波造成双重衰减。4. 避坑下载、解析、建模中 5 个血泪教训新手踩中一个就白干三天这些坑不是理论假设而是我在三个不同项目中反复翻车后记下的「后悔药清单」。它们隐蔽性强、报错信息模糊、调试耗时极长但只要提前知道10 分钟就能绕开。4.1 现象ObsPy read()报SACError: Invalid header但文件用 SAC 工具能打开原因NIED 下载的 SAC 文件头部nvhdr字段header 版本号被设为 6而旧版 ObsPy1.4.0只认nvhdr2或5。这不是文件损坏是版本兼容问题。解决升级 ObsPy 至1.4.0pip install --upgrade obspy或手动用sed替换头部sed -i s/\x00\x00\x00\x06/\x00\x00\x00\x05/g file.SACLinux/macOS。4.2 现象合并多台站数据后模型训练 loss 突然爆炸梯度 nan原因部分 K-NET 台站 SAC 文件的user0字段用户自定义浮点数被误写为极大值如1e30ObsPy 读取时会污染tr.data的 dtype导致后续 FFT 计算溢出。解决在read()后立即检查并清理tr read(...)[0] if np.any(np.isinf(tr.data)) or np.any(np.isnan(tr.data)): tr.data np.nan_to_num(tr.data, nan0.0, posinf0.0, neginf0.0)4.3 现象用scipy.signal.stft提取频谱图结果全是横条纹无时频聚焦原因SAC 文件单位是cm/s²但 STFT 默认输入为无量纲信号。未做归一化时强震记录峰值达 1000 cm/s²的数值远超弱震峰值 0.1 cm/s²动态范围超过 10⁴STFT 窗函数无法同时捕捉高低幅值成分。解决对每个 trace 单独做zscore归一化非全局归一化from scipy.stats import zscore tr.data zscore(tr.data, axisNone) # axisNone 确保按整个数组计算均值/标准差4.4 现象训练好的 P 波检测模型在验证集上 recall 很高0.95但 real-time 测试时漏检严重原因CSV 速报中的Origin Time是人工修订后的最终结果而真实预警场景中模型必须基于前 3–5 秒的实时流数据做判断。你用「完美时间标签」训练却用「不完美流数据」测试属于典型的标签泄露。解决构造「流式模拟数据集」——对每个 SAC只保留t0到t5s的片段作为输入用t5s后第一个 P 波到达时刻需用 STA/LTA 等传统方法标定作为监督标签而非直接用 CSV 中的Origin Time。4.5 现象haversine_vector计算台站距离时50km 半径内台站数忽多忽少甚至为 0原因CSV 中的经纬度是 WGS84 坐标系但部分老旧台站元数据如 NIED 提供的station_list.csv用的是日本本土坐标系JGD2000两者相差可达 10 米以上。小范围筛选时10 米误差可能导致台站被错误排除。解决统一使用pyproj转换import pyproj transformer pyproj.Transformer.from_crs(EPSG:4612, EPSG:4326, always_xyTrue) # JGD2000 → WGS84 lon_jgd, lat_jgd transformer.transform(sta_lon, sta_lat)5. 构建可复现的基线模型用 ResNet18 分类震级区间附完整训练脚本与验证技巧有了清洗后的 SAC 数据下一步是建立一个可量化、可对比、可快速迭代的基线模型。我推荐用 ResNet18 处理时频图而非原始波形因为1ResNet18 在小样本10k 样本下收敛快、过拟合风险低2时频图天然保留震相特征P/S 波能量分布比 raw waveform 更鲁棒3便于可视化注意力区域辅助分析模型是否真的学到了物理规律。这不是为了刷 SOTA而是给你一个「尺子」——后续所有新模型、新特征、新损失函数都必须在这个基线上证明自己确实更好。5.1 生成时频图数据集STFT 参数的 3 个必调参数STFT 是连接时域 SAC 与 ResNet 输入的关键桥梁。参数选错模型再好也是瞎学。经实测以下三参数组合在 2019–2024 日本地震数据上表现最优参数推荐值为什么这样设nperseg窗长256对应 2.56 秒100Hz 下足够覆盖 P 波初动上升沿又不过长导致时频分辨率下降noverlap重叠点数192重叠率 75%保证相邻帧间平滑过渡避免 STFT 图出现明显竖条纹nfftFFT 点数512提升频率轴分辨率0–50Hz 分辨率达 0.1Hz清晰分离 P 波~5–15Hz与 S 波~2–8Hz频带import numpy as np from scipy.signal import stft import torch from torch.utils.data import Dataset class EQSpectrogramDataset(Dataset): def __init__(self, sac_paths, label_map, transformNone): self.sac_paths sac_paths self.label_map label_map # e.g., {3.0-3.9: 0, 4.0-4.9: 1, ...} self.transform transform def __len__(self): return len(self.sac_paths) def __getitem__(self, idx): sac_path self.sac_paths[idx] tr read(str(sac_path), formatSAC)[0] # 归一化 裁剪同前 tr.data zscore(tr.data, axisNone) tr.trim(starttimetr.stats.starttime tr.stats.sac.b - 2.0, endtimetr.stats.starttime tr.stats.sac.b 15.0, padTrue, fill_value0) # STFT 生成时频图输出为 (freq, time) 的二维数组 f, t, Zxx stft(tr.data, fs100, nperseg256, noverlap192, nfft512) # 取绝对值谱转为 float32 spec np.abs(Zxx).astype(np.float32) # ResNet 输入需为 (C, H, W)此处 C1, Hfreq_bins, Wtime_bins spec torch.from_numpy(spec).unsqueeze(0) # (1, 257, 59) # 标签从文件名或 CSV 关联震级区间 mag extract_magnitude_from_filename(sac_path) # 自定义函数 label self.label_map[mag_interval(mag)] if self.transform: spec self.transform(spec) return spec, label # 震级分组函数按 JMA Mj def mag_interval(mj: float) - str: if mj 3.0: return 2.0-2.9 elif mj 4.0: return 3.0-3.9 elif mj 5.0: return 4.0-4.9 elif mj 6.0: return 5.0-5.9 else: return 6.0注意stft输出的Zxx是复数谱必须用np.abs()取模f有 257 个频率点0–50Hzt有 59 个时间点0–17 秒步长 0.29 秒因此spec.shape (1, 257, 59)不要做 log 压缩——ResNet18 的 BatchNorm 层对输入分布敏感log 后的稀疏分布会显著降低收敛速度。5.2 训练脚本核心冻结 backbone 两层 FC5 轮即可收敛ResNet18 在地震数据上极易过拟合我的经验是永远冻结前 4 个残差块layer1-layer3只训练 layer4 和分类头。这样既利用 ImageNet 预训练特征提取能力又避免小数据量下 backbone 权重坍塌。import torch.nn as nn import torchvision.models as models def create_resnet18_baseline(num_classes5): model models.resnet18(pretrainedTrue) # 冻结 layer1-layer3 for param in model.layer1.parameters(): param.requires_grad False for param in model.layer2.parameters(): param.requires_grad False for param in model.layer3.parameters(): param.requires_grad False # 修改第一层ResNet18 原输入为 3 通道 RGB改为 1 通道灰度 model.conv1 nn.Conv2d(1, 64, kernel_size7, stride2, padding3, biasFalse) # 修改分类头替换 fc 层为两层64→32→num_classes加 Dropout model.fc nn.Sequential( nn.Linear(model.fc.in_features, 64), nn.ReLU(), nn.Dropout(0.3), nn.Linear(64, 32), nn.ReLU(), nn.Dropout(0.3), nn.Linear(32, num_classes) ) return model # 训练循环精简版 model create_resnet18_baseline(num_classes5) criterion nn.CrossEntropyLoss() optimizer torch.optim.AdamW(filter(lambda p: p.requires_grad, model.parameters()), lr1e-3) for epoch in range(5): model.train() total_loss 0 for spec, label in train_loader: optimizer.zero_grad() out model(spec) loss criterion(out, label) loss.backward() optimizer.step() total_loss loss.item() print(fEpoch {epoch1}, Loss: {total_loss/len(train_loader):.4f})关键技巧filter(lambda p: p.requires_grad, model.parameters())确保只优化 unfrozen 层的参数AdamW比Adam更适合小批量训练权重衰减更稳定学习率设为 1e-3而非常见的 1e-4——因为只有 head 层在更新需要更快的收敛速度。5.3 验证模型是否真学到物理规律用 Grad-CAM 可视化注意力热区一个合格的地震模型其注意力区域应该集中在 P 波初动0–2 秒和 S 波主震5–10 秒附近而不是随机亮斑。用 Grad-CAM 可以直观验证这一点from pytorch_grad_cam import GradCAM from pytorch_grad_cam.utils.image import show_cam_on_image # 加载一张测试样本 spec, _ test_dataset[0] spec spec.unsqueeze(0) # (1,1,257,59) model.eval() with torch.no_grad(): pred model(spec).argmax().item() # 初始化 Grad-CAM target_layers [model.layer4[-1].bn2] # ResNet18 layer4 最后一个 block 的 BN 层 cam GradCAM(modelmodel, target_layerstarget_layers, use_cudaFalse) grayscale_cam cam(input_tensorspec, targetsNone)[0, :] # 叠加到时频图上需将 spec 归一化到 0–1 spec_np spec.squeeze().numpy() spec_norm (spec_np - spec_np.min()) / (spec_np.max() - spec_np.min()) visualization show_cam_on_image(spec_norm.transpose(1,0), grayscale_cam, use_rgbFalse) plt.imshow(visualization, cmapjet) plt.title(fGrad-CAM for predicted class {pred}) plt.show()如何解读如果热区集中在时频图左下角低频早期时间说明模型在用 P 波初动判别震级若热区在右上角高频晚期时间则可能在用 S 波振幅或面波这与震级物理定义一致。若热区散乱无规律立刻停训——数据或标签有问题不是模型问题。我带过的几个学生项目最初都卡在「模型训出来但不敢用」的阶段。后来我们统一加了一条铁律每次模型迭代必须跑一次 Grad-CAM截图存档和上一轮对比。三个月下来大家自然就建立起对模型行为的直觉——哪些改进真有用哪些只是过拟合噪音。这种「可解释性驱动的迭代」比单纯刷指标可靠得多。希望帮到你。本文还有配套的精品资源点击获取