把地震目录灌进SQLite做近场筛选,我踩过的三个坑

手头有一份约 12 万条的区域地震目录(CSV,字段:发震时刻、经纬度、深度、震级、震源位置精度),跑在 4 核 8G 的机器上,Python 3.11 + SQLite 3.40。需求很简单:给定一个台站坐标,筛出 50 公里内的所有事件,按震级降序,再按 10 公里深度分箱统计。第一版写得很快,结果错得也很快。

第一步是建表。别用默认的浮点比较去筛距离,先把经纬度存成 REAL,再建一个覆盖索引:

CREATE TABLE eq(
  id INTEGER PRIMARY KEY,
  t TEXT NOT NULL,
  lon REAL NOT NULL,
  lat REAL NOT NULL,
  depth REAL,
  mag REAL,
  locfix TEXT
);
CREATE INDEX idx_lonlat ON eq(lon, lat);

候选集用经纬度矩形框粗筛,再在 Python 里用 haversine 精算,这一步比在 SQL 里算球面距离快了一个量级:

import math
R = 6371.0
def hav(lat1, lon1, lat2, lon2):
    p1, p2 = math.radians(lat1), math.radians(lat2)
    dp = p2 - p1
    dl = math.radians(lon2 - lon1)
    a = math.sin(dp/2)**2 + math.cos(p1)*math.cos(p2)*math.sin(dl/2)**2
    return 2*R*math.asin(math.sqrt(a))

粗筛框的经度跨度要除以 cos(纬度)。我第一次直接用了固定 ±0.5 度,在高纬台站上漏掉了东侧事件——这不是精度问题,是框根本没覆盖到。

第二个坑是时间字段。CSV 里是 2023-05-12 03:41:07.2 这种带小数的字符串,我原样入库,然后按字符串排序,结果 .10 排在了 .2 前面。后来统一在入库时转成 epoch 秒(float)另存一列,字符串只留显示用。

第三个坑是深度。目录里空深度和 10 公里固定深度混在一起,分箱时我把 NULL 静默丢掉了,统计出来的浅源比例偏高。现在的做法是先查 SELECT depth IS NULL, count(*) FROM eq GROUP BY 1,确认缺失比例再决定是否单独列一组。

代价说清楚:矩形粗筛在台站靠近 180 度经线或极点时会失效,得退化成全表扫描;12 万条全扫大约 0.3 秒,还能接受,但上百万条就得换 R 树或直接上空间库。另外 haversine 假设地球是球,50 公里尺度误差在 0.3% 量级,够用,但别拿它去算震中距做走时反演。

最后一个没解决的问题:目录里同一事件有多个机构给出的位置,locfix 字段只标了精度等级,没有统一的事件 ID,去重时我按 30 秒时间窗加 100 公里距离聚类,边界上还是会出现一拆二。你们处理多源目录时是按什么规则合并的?

0
赞同度 100% · 0 条评论 · 5 阅读

0 条评论

支持 Markdown · 回复会实时推送给对方

还没有评论,写下第一条。