USGS地震数据抓取与可视化：从GeoJSON到Plotly的完整实战

简介USGS地震数据获取与预处理工具包面向地理信息或数据科学初学者解决从USGS平台批量下载地震JSON并将其转换为结构化CSV进行后续分析的问题。压缩包共4个文件包括Jupyter Notebook转换脚本HW2.ipynb、精简后的地震明细CSV、原始JSON数据以及README说明文档整体体积仅1.44MB便于快速部署本地环境。目前已有650人学习下载适合需要掌握公开地震数据接口调用、JSON解析与表格化整理的入门者。通过该包可完整了解使用wget命令行抓取USGS数据、借助JSONView检查返回字段以及利用Notebook执行数据清洗与降维存储的流程配套的CSV文件还可用作后续可视化或统计建模的直接输入节省从零检索API文档和调试格式的时间。1. 全球地震目录五分钟一更新USGS_earthquake 到底是干什么用的如果只盯着标题里的 USGS 三个字母很容易以为这是一个研究板块运动的学术项目实际上USGS_earthquake 更像是一个长期可用的地震数据总线。美国地质调查局的公开地震接口每 5 分钟同步一次全球台网捕捉到的震动事件任何人在不申请密钥的情况下就能拿到近实时目录、历史事件、震级修正记录和基础波形元数据。这个项目真正解决的问题不是帮你预测地震而是把散落在多个专用接口里的地震事件变成你本地数据管道里一行行结构化的记录。很多人第一次拉数据时会反直觉地发现两件事第一USGS 目录里大量事件是负震级或 0 级以下的微震直接画图会得到一片密密麻麻的噪点第二每条震例的时间字段实际上是一个 UTC 毫秒时间戳附带的时区偏移量常常为空新手用本地时间直接解析整张时间轴都会错位。USGS_earthquake 的价值就在这些细节里选对接口、理解字段语义、做对过滤和可视化。这篇文章按我自己的落地路径来写先讲数据源选型再给可复现的抓取和绘图脚本最后列出反复踩过的坑以及一个用古登堡-里克特定律校准数据质量的方法。适合的读者不是地质学专家而是要做监测大屏、灾害数据归档、地理分析或者单纯想用真实公开数据练手的数据工程师。2. 先分清楚数据源GeoJSON Feed、FDSN 事件服务与 CSV 出口怎么选USGS 对外提供了不止一个地震数据入口新手最容易犯的错误是拿到一个 URL 就开始写代码跑通之后就把它固化到生产环境里从不质疑这个接口到底覆盖什么范围、稳不稳定。我这边的做法是先做一次数据源选型把接口的实时性、回溯范围、字段完整度列成表再按业务场景决定用哪个。常见的选择有三个实时 GeoJSON Feed、FDSN 事件服务、以及带 formatcsv 的快捷导出口。三者都能返回震例列表但覆盖范围和使用限制差异非常大。2.1 三个常用端点的实时性与覆盖范围GeoJSON Feed 是大多数人最先接触的入口地址在 earthquake.usgs.gov/earthquakes/feed/v1.0/summary 下按时间窗口预生成了 all_hour、all_day、all_week、all_month 几个汇总文件格式固定为 GeoJSON。它的定位是给网页和轻量应用做实时展示数据来自自动定位系统通常在震动发生后几分钟内就会出现在 Feed 里更新频率大概每 5 分钟一次。好处是无需参数、无需分页、一个 URL 直接拉全量坏处是它只保留最近 30 天左右的数据并且所有事件都没有经过人工复核震级、位置随时可能被后续更新覆盖。FDSN 事件服务则是另一套更正式的查询接口路径在 /fdsnws/event/1/query遵循国际数字地震台网联合会的通用规范支持通过 starttime、minmagnitude、经纬度范围等参数做针对性查询也能回溯到几十年前的历史目录。研究级场景里USGS 最终发布的事件目录以它为准。这个接口返回的震例里每条事件还可能附带后续修订出来的多版本震级记录适合做震级的时序分析。它的代价是语法更复杂且对单次查询返回的记录数有限制不能像 Feed 一样一把梭拉全量。第三个 CSV 快捷出口本质上就是 FDSN 事件服务的一种输出格式URL 上指定 formatcsv 就能直接拿到逗号分隔文本。它的优势是省去 JSON 解析用 Excel 或 pandas 读取都很顺手缺点是字段会少掉一部分比如缺少详细的事件详情链接和多种震级信息做深度分析时会被迫回头再查单个事件。我平时只把它当作快速探查工具比如先确认某个时间窗内到底有多少条记录真正入库时还是会切回 GeoJSON。三个入口的适合场景可以这样概括实时监测页面用 Feed归档和回溯用 FDSN 事件服务临时看数据用 CSV 出口。选型时还有一个容易忽略的细节Feed 里的 all_month 虽然能拿到 30 天数据但它的数据源和 FDSN query 并不完全同步具体表现是同一个事件在同一时刻查出的震级可能差 0.1 到 0.3这个问题后面避坑章节会专门讲。2.2 过滤器参数时间、震级、经纬度框的边界条件如果说数据源选型决定你能不能拿到数据参数过滤就决定你拿到的数据是不是你要的那批。FDSN 事件服务最常用的参数我在下表里列一下方便写请求时对照。需要注意这些参数全部是 HTTP query string用 requests 时放进 params 字典即可别自己拼 URL。参数名作用常见取值示例starttime / endtime限定事件发生时间范围ISO 8601 格式2024-01-01T00:00:00minmagnitude / maxmagnitude按目录首选震级过滤minmagnitude5.5minlatitude / maxlatitude纬度范围南纬为负-45 到 45minlongitude / maxlongitude经度范围西经为负-180 到 180mindepth / maxdepth震源深度范围单位千米0 到 100eventtype事件类型可选 earthquake 或 otherearthquakeorderby结果排序time、time-asc、magnitudelimit / offset分页控制单页上限约 20000limit10000format输出格式geojson、csv、kml这里有三个边界条件需要特别说明。第一是时间范围FDSN 接口能回溯的数据可以到上世纪七十年代但超过三十年以上的早期目录震级类型以体波震级和面波震级为主跟现在的矩震级不可直接混用我在做长周期趋势分析时会刻意分时间段建模。第二是震级过滤的语义minmagnitude 过滤的是目录首选震级USGS 每次只选择一个震级作为展示代表可能是 Mww、Ml、Mb 中的任何一种如果你想专门筛矩震级需要在事件详情里再过滤一次。第三是经纬度框不支持跨日期变更线环绕地球比如把西经 175 度和东经 175 度放在同一个查询里接口会认为经度范围超过 180 度而报错或返回空结果正确做法是拆成两个查询再合并结果。深度参数同样是常见误区mindepth 和 maxdepth 的单位是千米不是米。全球事件里震源深度可以到几百千米浅源地震集中在 0 到 30 千米如果你想把浅源地震单独排查出来直接 minlatitude 和 maxlatitude 双框往往更稳定。最后提醒一句limit 参数默认值对 FDSN 来说通常是 10000 左右单次超过 20000 条记录容易触发服务端错误拉取历史数据时我习惯按周切片而不是按月一次拉完。2.3 返回字段的解析约定毫秒时间戳、首选震级和坐标顺序确认了接口和过滤条件下一步是理解响应里的字段语义。GeoJSON 格式的响应整体是一个符合 GeoJSON 规范的对象核心数据在 features 数组下的每个元素里每个元素代表一条地震事件。properties 对象里存放震级、地点描述和时间信息geometry 对象存放经纬度坐标和深度这样的嵌套结构让不少只处理过二维表的数据工程师第一次上手时很懵。最容易看走眼的字段是 time它不是一个字符串而是从 1970-01-01T00:00:00Z 开始计数的毫秒数所有时间戳都是 UTC 时刻本身不带时区。properties.tz 提供的是该事件当地时区相对 UTC 的分钟偏移量但这个字段在某些震例里是 null所以解析时不能依赖它做时间轴换算。我的处理习惯是统一先把毫秒转成 pandas 的 UTC datetime只在最终展示时按项目需要的时区做转换中间分析流程一概使用 UTC这样能避免一半以上的时间错乱问题。另一个关键字段是 properties.mag注意它是首选震级。USGS 会根据数据质量自动选择一个它认为最可靠的震级类型放到这里同一事件在不同时刻访问这个字段值可能因为人工修订而改变。如果你要研究震级稳定性或者做余震序列分析更可靠的做法是取 properties.detail 指向的单事件详情接口那里面列出了所有台网计算出的震级类型和对应数值。此外 geometry.coordinates 的顺序是经度在前、纬度在后、深度在最后跟常见地图 API 的 lat/lon 顺序相反解析时务必写成 lon, lat, depth否则所有震中会变成海里的异常点。字段语义这块我吃过一次亏在一张时间序列图上有 30% 的事件时间点对不上排查到大半夜才发现是 time 字段被当成了本地时间直接取小时而不是先转换到 UTC。建议拿到数据后的第一件事不是画图而是先打印前三条记录的原始字段做人工确认这一步能省下后面无数排错时间。3. 本地跑通最小抓取脚本从实时 Feed 到存量回填的两种写法选型做完接下来就是把数据真正抓到本地。整个抓取链路可以分为三层第一层用 requests 访问公开接口拿到 JSON第二层把嵌套的 GeoJSON 解析成扁平化的 DataFrame第三层处理历史数据回填时的分页和时间切片。很多教程只演示最小请求导致读者在真实场景下遇到数据量变大、请求超时、内存暴涨时完全不知道怎么扩展。这一章按我实际使用的顺序来写从最小可用脚本开始一步一步加到能反复执行的批量抓取任务。3.1 最小请求脚本请求头、超时与返回状态判断最小请求不需要任何第三方库直接用标准库 urllib 也能做但我普遍推荐 requests因为它处理 params 编码和 JSON 解析更省心。先拉最近 24 小时全球目录代码如下。这个脚本可以作为所有后续代码的骨架我一开始做某次全球事件可视化 Demo 时就是从这个版本出发的。import requests # 最近24小时全球地震目录GeoJSON格式 url https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_day.geojson headers { User-Agent: Seismo-Data-Lab/0.1 (contact: dataexample.com), Accept: application/json, } resp requests.get(url, headersheaders, timeout15) print(status_code:, resp.status_code) print(content_type:, resp.headers.get(Content-Type)) data resp.json() features data.get(features, []) print(event count:, len(features)) if features: props features[0][properties] print(sample properties:, props)请求头里的 User-Agent 在这里是必须养成的习惯。USGS 对没有 User-Agent 的匿名请求一视同仁但你带上一个能标明用途的标识出问题排查时更容易在服务日志里定位自己的请求。timeout 我设置成 15 秒因为全球目录文件不算大但偶尔会遇到网络波动15 秒已经能覆盖正常响应。status_code 只检查等不等于 200 还不够我通常还会看一眼 Content-Type 是不是 application/json有时候服务端会返回一个 HTML 错误页但状态码仍然是 200直接解析 JSON 就会抛异常。第一次跑这个脚本时如果看到 event count 是空或者只有几十条不用怀疑这是正常的。所有震例的 properties.mag 分布从负数到 7 以上都有全部画出来就是一大片噪点后续过滤要靠 minmagnitude 参数或者本地做条件筛选。另外这里的 features 数组顺序并不保证按时间排列需要按 properties.time 重新排序不要假设索引就是时间顺序。3.2 把嵌套 GeoJSON 解析成扁平 DataFrame拿到原始 JSON 之后下一步是把它整理成能进 pandas 的表格结构。GeoJSON 的嵌套层级深直接用 pandas.json_normalize 可以省不少代码但我在生产环境里更倾向于手写一个简单的提取函数原因是字段重命名和空值处理更可控而且代码可读性更好。以下是解析函数。import pandas as pd def parse_usgs_geojson(data): rows [] for f in data.get(features, []): props f.get(properties, {}) g f.get(geometry, {}) # GeoJSON坐标顺序为经度、纬度、深度 coords g.get(coordinates, [None, None, None]) coords (coords [None] * 3)[:3] # 防止坐标字段缺失 rows.append({ event_id: f.get(id), time_epoch_ms: props.get(time), tz_offset_min: props.get(tz), mag: props.get(mag), place: props.get(place), detail_url: props.get(detail), lon: coords[0], lat: coords[1], depth_km: coords[2], }) df pd.DataFrame(rows) # 毫秒时间戳统一转成UTC datetime if not df.empty: df[time_utc] pd.to_datetime(df[time_epoch_ms], unitms, utcTrue) return df # 接前面接口返回的 data df parse_usgs_geojson(data) print(df[[event_id, time_utc, mag, place, lat, lon, depth_km]].head(5))这段代码里处理了两个边界问题。第一个是 geometry 字段可能缺失某些特殊事件类型在目录中并不附带坐标我通过坐标列表填充 None 来保证列数对齐。第二个是 time_epoch_ms 转换为 time_utc 时显式指定 unitmspandas 默认是纳秒如果把毫秒时间戳直接传给 pd.to_datetime 而不带 unit结果会变成 1970 年附近的巨大时间值这是抓数据后最常见的低级错误之一。解析完成后一定要做一次数据质量检查我一般会先看 df.info() 和 df.describe()。一个正常的 24 小时全球目录应该有几千到上万条记录lat 范围在 -90 到 90lon 范围在 -180 到 180。如果发现 lat 出现 100 多这种越界值说明坐标解析时顺序写反了立即回头改代码不要带着脏数据往下走。tz_offset_min字段经常是整列的 None这符合预期不要因为它空就删掉整列在图谱上地震所在地理位置的时区分布仍可能用得到。3.3 翻历史目录FDSN 事件服务的分页与时间切片回填实时 Feed 只能拿最近 30 天做历史分析时必须切换到 FDSN 事件服务。这个接口支持 starttime、endtime 参数可以精确到秒也支持 limit 和 offset 做结果分页。限流方面虽然没有明确写死每分钟请求次数但官方对单次请求的返回量有上限体感上单次超过 2 万条记录容易超时或报错所以我更倾向于按时间切片让每次请求落在几千到一万条左右。import requests import time base https://earthquake.usgs.gov/fdsnws/event/1/query headers {User-Agent: Seismo-Data-Lab/0.1} # 复用前面的定义 def fetch_by_date_range(starttime, endtime, minmagnitude4.5): all_features [] offset 1 page_size 10000 while True: params { format: geojson, starttime: starttime, endtime: endtime, minmagnitude: minmagnitude, orderby: time, limit: page_size, offset: offset, } resp requests.get(base, paramsparams, headersheaders, timeout30) resp.raise_for_status() payload resp.json() feats payload.get(features, []) all_features.extend(feats) # 返回记录数小于单页上限说明没有下一页了 if len(feats) page_size: break offset 1 time.sleep(0.5) # 简单限速避免连续请求打满服务 return all_features # 示例取2024年1月全球4.5级以上事件 features fetch_by_date_range(2024-01-01T00:00:00, 2024-01-31T23:59:59) print(total events:, len(features))这个函数有几个值得细说的点。其一orderby 建议固定为 time这样分页的 offset 是基于时间顺序的游标不至于在翻页时因排序不稳定导致重复或遗漏。其二limit 直接设成 10000是为了让单页数据量适中既不会频繁翻页也不会因为太大而触发服务端限制。其三while True 加 len(feats) page_size 的退出条件当最后一页不足上限时自然停止逻辑直观且不会出现把空页当结果的情况。时间切片还有个优化技巧把一个月拆成四个周区间分别请求能减少单次请求的数据量同时避免因为临时网络问题丢掉整个月的任务。每个切片请求完成后把 features 合并再统一传给前面的 parse_usgs_geojson 函数生成 DataFrame。这个模式也适合做成定时任务每天跑一次前一天的时间窗把新增震例追加到自己的本地数据库里形成稳定的增量归档管道。4. 把震中从表格变成地图Plotly 散点地图的两个关键映射抓下来的数据最终要让人能一眼看懂。地震数据最直观的可视化形式就是震中分布图用气泡大小表示震级、颜色表示深度既能看到全球地震带轮廓也能突出目标事件。可视化实现我一般直接用 Plotly Express因为它交互性好、能导出 HTML对新手最友好而且散点地图 API 的设计和 pandas 接得很顺。但如果闭眼就画很容易遇到负震级导致气泡大小报错、跨日期变更线时地图中心错乱之类的坑。这一章只讲两个关键映射怎么把震级正确映射成气泡大小怎么让地图自动居中到数据范围。4.1 用震级映射气泡大小、用深度映射颜色的实现先做数据清洗再画图。直接从 DataFrame 画 scatter_geo 时如果 mag 列里有负值Plotly 会报错因为气泡 size 属性不允许负数。所有负震级事件在地图上都应该被视为最小可见点不能直接删掉否则全球小震分布信息就丢了。我的处理方式是做一次归一化把最小震级映射为 1这样所有点都成为合法正数。import plotly.express as px # 去除无法定位的行这是画图前最基础的清洗 plot_df df.dropna(subset[lat, lon, mag]).copy() # 把震级归一化到正数区间避免Plotly气泡尺寸报错 min_mag plot_df[mag].min() plot_df[bubble_size] (plot_df[mag] - min_mag) 1.5 plot_df[bubble_size] plot_df[bubble_size] * 2 # 放大整体气泡 fig px.scatter_geo( plot_df, latlat, lonlon, sizebubble_size, colordepth_km, hover_nameplace, color_continuous_scaleViridis, projectionnatural earth, title全球震中分布颜色深度大小震级, ) fig.show()这段代码里bubble_size 是显式构造的列用来绕过 Plotly 对 marker.size 必须为正数的限制。乘数 2 只是视觉上的放大系数实际项目中我会根据数据点的密集程度调整比如只画 5 级以上时就调成 3 或 4让气泡更有辨识度。color 用 depth_km 而不是震级是因为震级已经由气泡大小表达了颜色再放震级会造成信息冗余深度则能体现不同地震带的发震特征比如环太平洋深源震和浅源震在图上的颜色差异非常明显。需要注意px.scatter_geo 对数据量有一定的性能瓶颈两三千个点交互还流畅但如果拉了一整个月的全球小震数据几万点一起渲染浏览器会很卡。我的实践是把可视化和数据筛选分离画图前先在 DataFrame 上做一次条件过滤比如只保留 mag 3 的事件或者按自己关注的区域做经纬度切片。这样图上的信息密度反而更高比全量渲染更加可读。4.2 跨日期变更线时让地图自动居中的技巧地图自动居中是个容易被忽略的问题。Plotly 的默认地图中心是 0 度经线附近当数据集中在西半球或跨越大西洋时默认视图会显得很偏。手动指定 center 也不是不行但如果数据分布动态变化每次手动改经纬度就是在找死路。Plotly 提供了一个 fitbounds 参数可以让地图自动匹配当前数据点范围效果类似点图自适应。fig.update_geos( fitboundslocations, projection_typenatural earth, showcountriesTrue, showcoastlinesTrue, ) # 如果数据量小可以关闭地图背景只留数据点 fig.update_layout(margin{r: 0, t: 40, l: 0, b: 0}) fig.show()fitbounds 的取值是 locations它会根据 scattergeo trace 里的经纬度范围自动计算中心点这比手动指定 center 要省心得多。但是有一个坑当数据跨越日期变更线时比如同时包含 175°E 和 -175°W 的事件Plotly 会认为经度差接近 350 度缩放到整个世界地图整个图就变得很小。处理方式是在抓取时就把跨线区域的事件经度统一转换到 [0, 360] 区间或者直接按两半分别画两张子图。我的习惯是优先拆成子图因为转换经度会让图上位置和原始数据对不上反而引起歧义。出图之前还有一道工序把生成的可视化导出成 HTML方便浏览器里查看和分享。Plotly 导出 HTML 是fig.write_html(global_map.html)这个文件是自包含的不需要再装环境适合把结果发给同事看。如果要嵌入到自己的 Web 服务里也可以换成fig.to_json()传给前端渲染但那是另一个话题了这里不展开。5. USGS_earthquake 实战避坑五个翻车现场与排查清单我在初期处理 USGS 数据时踩过不少坑有些坑当时花了一整天才排查清楚。这些坑几乎没有例外全集中在对字段语义的理解、对不同接口特性的把握和对数据分布的错误预期上。这一章把它们单独拿出来按现象、原因、解决的顺序写清楚。你可以把这个章节当成一份排查清单每遇到一次异常就往这里对号入座。5.1 时间轴整体偏移几个小时时区字段还忽有忽无现象画出的震中分布图没问题但按时间序列统计每小时事件数时峰值时刻和官方页面完全不对齐整个曲线平移了几个小时。我曾在一次数据看板调试里发现所有事件时间都慢了八小时整整一个下午的时间都在质疑自己的统计代码。原因USGS 返回的 properties.time 是 UTC 毫秒时间戳本身是绝对时刻。我在解析后直接调用了本地时间的小时属性相当于把 UTC 时刻当成了东八区去读所以整条时间轴全都偏移。properties.tz 字段只在部分事件中有值不能指望用它来统一校准。解决所有内部处理统一用 UTC转换成df[time_utc] pd.to_datetime(df[time_epoch_ms], unitms, utcTrue)只有在图表展示的最终环节才把 UTC 转成目标时区转换时用 pandas 的dt.tz_convert不要用字符串替换之类的土办法。这样无论项目部署在哪台机器上结果都不会因为系统时区不同而再次漂移。5.2 GeoJSON Feed 和 FDSN 查同一事件震级对不上现象实时监控页面上显示某事件震级 5.3第二天用 FDSN 事件服务查同一个 event id震级变成了 5.6位置也有几百米的偏移团队以为数据管道出了问题。原因这是两类接口的定位差异。Feed 来自自动处理系统速度快但精度有限FDSN 目录里包含的是经过人工复核和重新定位的正式结果震级修正很常见。自动定位的震级波动 0.3 到 0.5 并不奇怪尤其对远震和大震。解决明确接口分工实时展示层允许用 Feed但所有进入历史库的数据必须以 FDSN 查询为准。入库时如果按 event id 做增量更新要设计成更新覆盖语义而不是不存在才插入。这样当同一事件出现修订后的震级时库里保存的是最新可靠值而不是第一次抓到的临时值。不要试图在 Feed 和 FDSN 之间做数据对齐两边本来就是两种口径。5.3 minmagnitude0 拉了一个月数据量大到本地直接卡死现象请求参数只写了 starttime 和 endtime没加震级过滤返回的 features 数量超过十几万解析成 DataFrame 后内存占用飙升后续画图完全卡住浏览器直接白屏。原因USGS 目录里包含大量零级以下的小震和微震全球每天可以达到几千条。不加过滤时一个月的数据量轻松达到几十万条远超地图可视化的合理承载范围。这不是接口的 bug而是对数据规模缺乏预估。解决先跑一次带 formatcsv 的快速查询确认总量再决定时间切片粒度。比如先查 2024-01-01 到 2024-01-07 的记录数如果一周就超过两万条就把窗口切成按天拉取。可视化前至少先做 minmagnitude 过滤我通常是从 2.5 或 3.0 起步再根据分析目的逐步下调。高速下载时空闲时间可以加一句time.sleep(0.3)既不惹恼服务器也给自己留调控余地。5.4 震中全部落在海里还在赤道两边对称现象地图上散点几乎全部分布在太平洋和赤道附近大量点落在海中呈一条直线陆地上的明显地震带却基本没有点。第一反应是 USGS 数据有误但网页上明明显示得很正常。原因这是 GeoJSON 坐标系顺序被写成 lat, lon 的结果。GeoJSON 规范明确规定 coordinates 是 [经度, 纬度, 深度]但很多用习惯了地图 API 的人在解析时习惯性先填 latitude 再填 longitude一交换所有点都被投影到了错误的位置。赤道两侧对称是因为纬度被填进了经度槽经度值超过 90 的部分被强制包裹后出现了镜像效果。解决在解析函数里固定坐标映射顺序并加一道断言防线。解析完成后检查df[lat].between(-90, 90)和df[lon].between(-180, 180)的通过率如果发现越界值直接抛异常终止流程。这样即使未来换了人维护代码、改了字段名也能在最早阶段发现坐标错位而不是等到地图上出现诡异图形才回头查。5.5 URL 参数没编码带空格的请求返回 400现象手写 URL 拼 request 时用了 starttime2024-01-01 00:00:00中间带空格接口直接返回 HTTP 400日志里只写了 Bad Request没有任何更具体的提示。查了一圈文档也没找到原因。原因HTTP GET 请求的 query string 里空格必须编码为%20或。很多教程里用 requests.get(url, paramsparams) 传入字典requests 会自动处理编码于是没人注意到这个问题。自己改成字符串拼接后站点的合法性检查就会卡在未编码的空格和冒号上。解决统一用 requests 的 params 字典参数不手拼 URL。如果确实需要生成带签名的 URL用 urllib.parse.urlencode 对参数做编码处理再把编码后的 query string 拼到基础地址后面。这个坑不大但会让人在排错时白白浪费不少时间放在这里提醒一句。6. 进阶用古登堡-里克特定律拟合 b 值反向检查数据完整性抓下来的数据质量到底怎么样不能光靠肉眼扫地图也值得用一个简单的定量手段来验证。地震学里的古登堡-里克特定律描述震级和频次的关系log10N a - bM其中 N 是大于等于震级 M 的事件数b 值正常范围大约在 0.8 到 1.2。b 值异常偏大通常意味着小震级事件在目录里大量缺失数据截止不完整。利用这个关系做一次线性回归就能找到自己数据集的完整性下限 Mc也能验证接口是否把该收的事件都收进来了。import numpy as np mags plot_df[mag].dropna().values # 震级分桶统计每档事件数量再计算累积频次 bin_edges np.arange(np.floor(mags.min()), np.ceil(mags.max()), 0.1) counts, _ np.histogram(mags, binsbin_edges) counts counts[::-1] # 从大到小累积 cum_counts np.cumsum(counts) # 取中值作为该震级档的代表 bin_centers (bin_edges[:-1] bin_edges[1:]) / 2 logN np.log10(cum_counts) # 只拟合大于完整性下限Mc的数据段这里以3.0为例 mask bin_centers 3.0 coeff np.polyfit(bin_centers[mask], logN[mask], 1) b_value -coeff[0] print(fb值: {b_value:.2f})这段代码的核心是最小二乘拟合np.polyfit 返回的斜率取负号就是 b 值。拟合前要把所有低震级段剔除否则低震级缺失带来的曲线弯曲会拉高斜率b 值失真到 1.5 甚至更高。完整性下限 Mc 的判断方法很简单把拟合直线和实际累积曲线画在同一张图上曲线开始明显偏离直线的那个震级就是 Mc。低于 Mc 的数据在统计上不可靠做趋势分析时应丢弃。我自己的习惯是每个区域数据集都算一次 b 值如果某个月的 b 值突然跳到 1.4 以上就会先检查抓取任务是不是断了再怀疑数据源本身的问题。这一步做下来USGS_earthquake 从拉数到验证的完整链路就通了。项目日后真正起作用的不是那条地图而是你建立的这套数据处理习惯先确认接口语义再写防御性解析最后用科学的统计方法做自查。我从一开始直接拿全量震级拟合跑出过 1.8 的荒唐 b 值后来才意识到低震级段的漏报全都被带进了回归里。这个错犯过一次之后我再也不会跳过 Mc 判断直接套公式。希望这篇笔记里的思路和脚本能帮你少走这段弯路也希望你在自己的数据管道里早日跑出稳定的 b 值。本文还有配套的精品资源点击获取