汶川地震泥石流数据集复现

Kanwuqing2026-09-262026-09-27

复现了 Wang 等人 2022 年发在 Scientific Data 上的一个数据论文, 数据在 Zenodo 上, 内容是 2008 年汶川地震之后龙门山一带的泥石流事件和降雨记录

论文: Two multi-temporal datasets to track debris flow after the 2008 Wenchuan earthquake
论文 doi: 10.1038/s41597-022-01658-y
数据 doi: 10.5281/zenodo.6891244

代码在 kanwuqing/debris-flow-reproduction

这篇论文本身只发数据, 不给模型也不给阈值, 所以下面的模型和阈值都是我自己写的, 不是抄论文的

数据

压缩包解压放到 data/raw/, 里面是这样

data/raw/
├─ dataset1/Dataset1 Debris flow and Triggering rain/
│  ├─ Dataset1_Debris flow.csv              186 条事件
│  ├─ Debris_flow_catchment.{shp,dbf,prj}   79 个流域多边形
│  └─ ID_RG_DF_DATA/ID_RG_DF_<DF_ID>/*.txt  465 个逐小时降雨序列
├─ dataset2/Dataset2 Dam/Dataset2_Dam.csv   拦砂坝
├─ Read me.pdf                              官方字段说明
└─ References of data sources.pdf

跑通了之后统计出来是 186 条事件, 79 个流域, 465 个降雨文件, 分布在 147 个文件夹里, 涉及 39 个雨量站, 拦砂坝 55 座分布在 23 个流域, 流域面积加起来 851.0

逐年事件数差别很大

年份	2008	2009	2010	2011	2012	2013	2016	2019	2020
事件数	29	6	98	22	1	16	1	12	1

2010 年占了 98 次, 一半还多, 对应 8·13 那次群发. 2014 2015 2017 2018 四年一次都没有

环境

Windows 11, Python 3.12, 全是传统机器学习, 没有深度学习

1	pandas, numpy, scikit-learn, matplotlib, seaborn, xgboost(可选)

随机种子固定成 20220901, 所有划分和所有模型的 random_state 都从它派生

读数据遇到的坑

坐标列反了

事件表 Dataset1_Debris flow.csv 里, 叫 Latitude 的那一列装的是 103.x, 叫 Longitude 的装的是 31.x

汶川在四川, 纬度 31 经度 103, 所以这两列是名字和内容对调了

要是不管它直接拿去做空间归属, 所有点都会跑到地球另一边去. 读进来的时候把两列名字换回 lon / lat, 原始的留着备查

一半事件没有时刻

186 条里有 95 条的 Time(24 hour) 是空的

用当天 00:00 顶上, 另外标一个 time_missing = True. 因为时刻不准, 判断"这场雨有没有触发泥石流"的时候就用整个日历日当窗口, 不精确到小时

事件名和流域名对不上

事件表里一个事件是一个沟名, 比如 Mayangdian

流域表里一个流域是好几个沟用 / 连起来, 比如 Douyaping/Mayangdian/Maliuwan/Qingling/Guanshan 对应 W16

所以先按名字劈开建索引做最长匹配, 匹配不上的再用点在多边形里判断. 最后 186 条里 167 条找到了流域, 占 89.8%

切雨场

每个事件文件夹里是周围雨量站从事件前 7 天到事件后 1 天的逐小时降雨, 这段时间里不可能只有一场雨

这其实是好事, 因为切成很多场之后, 同一批数据里天然就有负样本了, 不用自己去造

切分规则都写成常量放在 features.py 顶部

常量	值	意思
`WET_MM`	0.1	这个小时雨量超过 0.1mm 才算湿小时
`DRY_GAP_H`	6	连着 6 小时不下雨就切成两场
`MIN_TOTAL_MM`	1.0	总雨量不到 1mm 的场次丢掉
`API_K`	0.9	前期降水指数的小时衰减系数

每场雨算 16 个特征, 总雨量, 历时, 湿小时数, 最大 1/3/6/12/24/48/72 小时雨量, 平均雨强, 峰值雨强, 前期 24/72/168 小时雨量, 还有 API

和泥石流当天重叠的场次标成正样本, 同一窗口里其它场次是负样本, 最后 1597 行, 正 444 负 1153

流域特征

流域是多边形, 自己写了个 SHP 读取器解析, 没装 geopandas

一开始是想装 geopandas 的, 但是担心别人复现的时候卡在装 GDAL 这种二进制依赖上, 就自己写了, 顺便把 DBF 也写了

从多边形能算出来面积, 周长, 紧凑度, 包围盒长宽, 质心经纬度, 支沟个数, 加上数据集 2 的坝数量和建成年份

数据集里没有 DEM, 所以没做坡度高程这些特征, 缺就是缺, 不编

任务 A: 降雨阈值

行和行不独立

一开始我直接分层随机切 7:3, 跑出来随机森林 0.9625, 感觉挺好

后来发现不对劲. 一次泥石流事件会贡献好几个雨量站的文件, 同一场雨在好几个站上都被算了一遍, 随机切的话这一场雨可能一半在训练集一半在测试集里, 这就是泄漏

所以改成同时报三种切法

切法	怎么切	作用
`random`	分层随机 7:3	有泄漏, 只当乐观上界看
`group_event`	按 DF_ID 分组 GroupKFold	主协议, 同一次事件不跨折
`temporal`	2013 年前训练, 2016 年后测试	看能不能跨时间

结果

主协议下

模型	accuracy	precision	recall	F1	ROC-AUC
logistic_regression	0.9036	0.8059	0.8604	0.8322	0.9225
random_forest	0.9562	0.9329	0.9077	0.9201	0.9832
svm_rbf	0.9180	0.8323	0.8829	0.8568	0.9340
gradient_boosting	0.9555	0.9388	0.8986	0.9183	0.9735
xgboost	0.9580	0.9415	0.9054	0.9231	0.9802

随机森林混淆矩阵是 [[1124, 29], [41, 403]], 漏报 41 个

随机切法随机森林 0.9625, 只比主协议高 0.63 个百分点, 说明泄漏确实有但是没到让结论崩掉的程度

特征重要性

排前面的

特征	重要性
`n_wet_hours`	0.1739
`duration_h`	0.1715
`max_24h`	0.0899
`max_72h`	0.0723
`total_mm`	0.0681
`api`	0.0643

单小时峰值 max_1h 只有 0.0217, 排到 14 名去了

看分布也是这个意思, 触发场总雨量中位数 49.10mm, 非触发场 9.00mm, 差 5.5 倍. 但是平均雨强 2.62 和 2.15, 只差 1.2 倍

所以判断会不会触发, 主要是看下了多久下了多少, 不是看下得多急

temporal 切法不能信

2016 年之后只有 57 行测试数据

跑出来逻辑回归 0.9825 随机森林 0.8947, 差得离谱. 57 行里 23 个正例, 一个样本就值 1.8 个百分点, 指标乱跳很正常

这个协议我不用, 也没法拿它论证"阈值随震后恢复"这种事

I-D 阈值

历时怎么定义

泥石流阈值一般是这个形式

我一开始拿整场雨的时长当 , 雨强就是总雨量除以时长. 跑出来不行

因为龙门山这边季风雨经常连着下几十个小时, 干间隔 6 小时切不开, 一场雨历时中位数 20 小时最长 82 小时, 点云被少数超长样本拽着走, log-log 相关只有 0.274, 斜率完全是那几个极值点定的

这里还有个事, threshold.py 注释里原来写的是"几乎零相关 r = -0.08", 图 1 标题里写的是 r = 0.01, 这俩数我都复算不出来, 实测是 +0.274. 应该是早期草稿写上去忘了改, 都改掉了

改成固定窗长: 对取 1 3 6 12 24 48 72 小时, 每个窗长找任意小时窗口里的最大雨量,

这样一个雨场贡献 7 个点, 一共 11179 个对

拟合结果

一共拟合了 6 条, 和机器学习用同一批行同一批标签打分

名称	拟合	accuracy	precision	recall	F1
`ols_all_points`	0.3503	0.6456	0.4361	0.9369	0.5951
`ols_trigger_only`	0.4855	0.6807	0.4446	0.5968	0.5096
`ols_trigger_minus_1sd`	0.4855	0.6606	0.4473	0.9369	0.6055
`quantile_p50_trigger`	0.9909	0.7120	0.4881	0.7365	0.5871
`quantile_p90_trigger`	0.9649	0.7733	0.9271	0.2005	0.3296
`quantile_p90_all`	0.9911	0.7558	0.6154	0.3243	0.4248

这里有个事挺意外. 把 11179 个点直接丢进最小二乘, 只有 0.3503, 只用触发点也才 0.4855

但是先按窗长取分位数再拟合, 就跳到 0.9649 到 0.9911

所以幂律只在分位数包络上成立, 直接对原始点云做 OLS 得到的那条线形式对参数不可靠, 同一个的说法在两种做法下差别很大

最高的 quantile_p90_trigger accuracy 0.7733 看着不错, precision 0.9271 也挺好, 但是 recall 只有 0.2005, 444 个触发场漏了 355 个, 当预警阈值没法用

把线压低倒是能把 recall 拉到 0.9369, 但是 precision 掉到 0.4361, 报十次警五次半是空的

经验曲线 F1 最高 0.6055, 随机森林 0.9201, 差了 0.31. 毕竟经验曲线只用和两个量, 随机森林用了 16 个特征

一个失败的尝试

试了下从随机森林的 0.5 概率等高线反推一条 I-D 曲线, 做法是对每个窗长造合成雨场然后扫峰值深度, 找预测概率第一次超过 0.5 的地方

结果返回空字典

因为合成样本只有一个特征被赋值, 其它 15 个填 0, 模型在这种输入上概率不单调, 找不到稳定的穿越点

1	ml_equivalent_threshold -> {}

没去调参把它凑出来, 失败就写成失败

任务 B: 易发性

79 个流域一行一个, 标签是事件数大于中位数的算易发. 中位数是 1.0, 所以正类 35 个负类 44 个

跑之前要把 n_events first_event last_event total_volume_m3 event_span_years 这几列删掉, 它们都是标签算出来的, 留着就是自己骗自己

RepeatedStratifiedKFold 5 折乘 10 次

模型	accuracy	precision	recall	F1	ROC-AUC
logistic_regression	0.7598 ± 0.0970	0.7543	0.7429	0.7307 ± 0.1079	0.8082
random_forest	0.7277 ± 0.0847	0.6920	0.7457	0.7001 ± 0.1125	0.7951
svm_rbf	0.7115 ± 0.1129	0.6833	0.7057	0.6738 ± 0.1504	0.7822
gradient_boosting	0.7140 ± 0.1091	0.7043	0.6800	0.6737 ± 0.1262	0.7992
xgboost	0.7326 ± 0.0803	0.7101	0.7029	0.6927 ± 0.1070	0.8186

标准差都有 0.08 到 0.11, 模型之间 accuracy 差 0.7115 到 0.7598, 还不到一个标准差

所以不能说逻辑回归比随机森林好, 只能说这些特征大概有 0.78 到 0.82 的排序能力, 而且不稳

特征重要性里 centroid_lat 最高 0.1543, n_gullies 第二 0.1272. 逻辑回归里 n_gullies 系数最大是 +1.95, 支沟越多越易发, 这个符合直觉

经纬度重要性高大概率是因为它俩在代替"离震中多远", 不是纬度本身有什么作用, 不往因果上解读

坝有没有用

坝数	流域数	事件数	平均事件数
0	56	105	1.88
1	9	25	2.78
2	3	9	3.00
>=3	11	28	2.55

有坝的平均事件数 2.55 到 3.00, 无坝的 1.88, 有坝反而更高

这个不能读成"坝导致了泥石流", 是选址偏差, 坝本来就被建在灾害多的地方

另外两个数能佐证, 面积和有没有坝相关 0.5133, 大流域更可能有坝; 面积和事件数相关 -0.0175, 几乎无关. 所以有坝流域事件多这件事不能拿面积解释掉

1 2	prevented = False 62 个流域 134 次事件平均 2.16 prevented = True 17 个流域 33 次事件平均 1.94

prevented 这两行差了 0.22 次, 79 个流域的样本量下这个差不算证据, 不解读

想评估坝的效果得做建坝前后的对比, 不是横截面比. 数据里有 built_year, 但是建成年集中在 2011 到 2014, 建完之后观测年份太少, 功效不够, 做不出结论

数字是怎么对上的

中间吃过亏, 报告草稿里有几个数字是我自己记岔了写上去的, 和结果文件对不上

后来定了个规矩, 报告里引用的每个数字都必须能从 CSV 或者 JSON 里查到, 不凭记忆写

1 2	scripts/explore/report_numbers.py -> docs/_numbers.md scripts/explore/audit_artefacts.py -> docs/_audit.txt

改过分析方法之后就重跑一遍再核一遍. 这轮核出来改了 6 处

位置	原来写的	实测
报告 OLS	0.0200 / 0.0197	0.3503 / 0.4855
报告分位数	0.9204	0.9909 / 0.9649 / 0.9911
`threshold.py` 注释	r = -0.08	+0.274
图 1 标题	r = 0.01	+0.274
正样本数	439	444
无时刻事件数	96	95

make_figures.py 里图 1 的标题原来是写死的字符串, 现在改成运行时算了

一键复现

1
2
3

python -m venv .venv
.\.venv\Scripts\python.exe -m pip install -r requirements.txt
.\.venv\Scripts\python.exe run_all.py

run_all.py 顺序跑 5 个阶段, 最后自动做 SHA-256 校验

跑了两遍, 16 个产物全部一样

1 2	manifest artefacts: 16 \| present: 16 OK: every artefact reproduces bit-for-bit.

整个流程 49 到 62 秒

踩的坑

后台跑 python 不能接管道

1	ResourceUnavailable: Program 'python.exe' failed to run: 拒绝访问

后台跑的时候只要加了 | Tee-Object 或者 *> 就会这样, 不加就没问题. 换成 Start-Process 加重定向也可以

无头 Edge 导 PDF 被拦

1	Program 'msedge.exe' failed to run: 拒绝访问

加上 --user-data-dir 指定到仓库里的临时目录之后就能跑了

结果文件

results/metrics_summary.csv              全部模型的指标和混淆矩阵
results/task_a_rainfall_threshold.json   任务 A 全部结果, 含三种切法和 6 条曲线
results/task_b_susceptibility.json       任务 B 全部结果, 含 79 个流域明细
results/figures/fig01..07*.png           7 张图
results/repro_manifest.json              16 个产物的 SHA-256
docs/report.md   docs/report.pdf         实验报告
docs/repro_log.md                        这篇

报告里所有图上的数字都能在 CSV 和 JSON 里查到, 不用去读图

引用

Wang L., Chang M., Le J., Xiang L., Ni Z. (2022). Two multi-temporal datasets to track debris flow after the 2008 Wenchuan earthquake. Scientific Data 9:525

数据 doi: 10.5281/zenodo.6891244, 许可 CC BY 4.0

这个形式是 Caine 1980 的经典写法

代码都是自己写的, 没抄别人的实现. 固定窗长最大雨量构造对, 按事件分组防泄漏, 用中位数切易发性标签, 这几样是这行的通行做法

报告里所有阈值系数都是这次数据自己拟合的. 之前翻到过一些汶川震区的 I-D 系数, 但是只拿到标题和链接, 没找到全文, 就没拿来当基线