# 华北重磁网格的多源特征融合：空间分块交叉验证下，融合没有增益，随机划分会造出一个假增益

**问题**：基于公开重磁网格数据运行多模态特征融合基线模型，给出交叉验证结果、运行日志与基线性能对比。
**落成的口径**：华北 105–120°E / 35–50°N，2 弧分网格（约 3.7 km），188,294 个有标注格；
重力 = GGMplus 地面重力扰动；磁 = NOAA EMAG2v3 上延 4 km 磁异常；标注 = USGS 概化地质图 geo3al 落成 4 类；
模型 = 梯度提升树（默认超参数，不调参），仅重力 / 仅磁 / 早期融合 / 后期融合，外加多数类平凡基线。
**这一轮真跑了**：3 种划分 × 3 个种子 × 4 个模型的 5 折交叉验证；两种零模型（标注置换 9 次、环形平移 27 次）；
正负两把尺子对照各 3 个种子；按 1° 空间块整块重抽 1000 次的区间与配对 p；纯标准库独立复算。纯 CPU。

---

## 0. 结论

**三句话：**

1. **融合没有带来实用增益。** 1° 空间分块下，早期融合（重力 + 磁 8 个特征拼接）平衡准确率（BA）为 **0.4546**，
   仅重力为 **0.4533**，差 **+0.0012**。三个种子的差分别是 +0.0005 / +0.0000 / +0.0032，
   按空间块重抽的 95% 区间都跨 0，配对 p = 0.944 / 0.995 / 0.604。后期融合（两个单源模型概率取平均）差 −0.0048。
   换成更保守的 2.5° 分块，早期融合反而比仅重力**低 0.0179**（三个种子 p = 0.019 / 0.009 / 0.028）。
2. **如果用像元随机划分来评估，会看到一个其实不存在的融合增益。** 同样的模型和数据，随机划分下早期融合 BA 为 **0.5517**，
   比 1° 分块高 **0.0971**；此时它比仅重力高 **+0.0743**（三个种子 p 均为 0.001，区间下界都在 0.062 以上）。
   随机划分只把仅重力抬高了 0.0240，却把融合抬高了 0.0971 —— 特征越多，模型越容易靠「记住测试格子的邻居」得分。
   **在地学网格上比较融合方案时，随机划分下看到的增益，在本轮数据里几乎全部来自空间泄漏。**
3. **能测到的信号来自重力，磁数据在这把尺子上几乎测不出东西。** 仅重力比保留空间结构的环形平移零模型高 0.1076 BA；
   仅磁只高 0.0057（z = 0.48），而且低于零模型 9 次复本里的最大值 0.2916。
   已知答案的尺子证明这套评估测得出融合增益：标注同时依赖重力与磁时，融合比单源高 +0.36（p = 0.001）；
   标注只依赖重力时，融合与单源的差为 0.0000。所以「融合没用」不是评估方法测不出来，
   而是在这个尺度上，这份磁数据没有给出重力之外的新信息。

**这一轮不能支持的推断：**

- **不能**推出「重磁融合在地学分类中普遍无用」。本轮只有一种磁数据（EMAG2v3 在这一带是偏平滑的长波长编译数据）、
  一个区域、一种模型、一份 1:数百万比例尺的概化地质图。换成分辨率更高的航磁数据，结果可能不同。
- **不能**推出「重力能识别岩性」。GGMplus 是未做布格改正的重力扰动，里面混着地形；新生代覆盖层多在盆地低处，
  重力的一部分本事可能是在看地形。本轮没有高程数据把这两者拆开。
- **不能**外推到矿区尺度、成矿预测或三维深部结构。
- 约 0.45 的 BA 不是可用的地质填图精度，它只用于同一把尺子上的相对比较。

**跑之前写下的预测与判定：**

| 编号 | 跑之前写下的预测 | 判定 | 连续量 |
|---|---|---|---|
| P1 | 1° 分块下，早期融合相对最强单源的 ΔBA 均值落在 [−0.02, +0.02] | **成立** | ΔBA 均值 +0.0012；三个种子的区间都跨 0；p = 0.944 / 0.995 / 0.604 |
| P2 | 随机划分把早期融合的 BA 抬高至少 0.10 | **部分成立**（落在 0.05–0.10 档） | 抬高 +0.0971 |
| P3 | 仅重力与早期融合都高于全部零模型复本，且高出零模型均值 ≥ 0.05 | **成立** | 六组里最小的高出量为 0.0958（仅重力，种子 0，零模型最大值 0.3913） |
| P4 | 仅磁比它自己的环形平移零模型均值高不到 0.05 | **成立** | +0.0057，z = 0.48 |
| 尺子 | 正对照测得出（ΔBA ≥ 0.10 且区间下界 > 0），负对照测不出 | **通过**，P1 的结论有效 | 正对照 +0.3627 / +0.3617 / +0.3610，p = 0.001；负对照 0.0000 / −0.0001 / −0.0001，p = 0.504 / 0.109 / 0.195 |

**一处披露**：为估算机时，写预测之前用 1 个种子、1° 分块跑过一次计时试验，看到的是仅重力 BA 0.451、
仅磁 0.280、拼接融合 0.449。P1 和 P4 是看过这组数字之后写的，**不是盲预测**；严格意义上的事前预测是 P2、P3 与两把尺子。
这一条在预注册里原样写着。

## 1. 数据、标注与网格

| 角色 | 数据 | 分辨率 | 来源 |
|---|---|---|---|
| 重力 | GGMplus 地面重力扰动（Hirt et al. 2013） | 原生 0.002°（约 220 m），按格内全部像元求均值到 2 弧分 | https://ddfe.curtin.edu.au/gravitymodels/GGMplus/data/dg/ |
| 磁 | NOAA EMAG2v3（2017）上延 4 km 版 | 原生 2 弧分 | https://www.ngdc.noaa.gov/geomag/data/EMAG2/EMAG2_V3_20170530/ |
| 标注 | USGS OFR 97-470F《远东概化地质图》geo3al | 矢量多边形，1:数百万比例尺 | https://www.sciencebase.gov/catalog/item/60abc7f9d34ea221ce51e5ee |

**为什么不是 WGM2012**：WGM2012 的网格提取页带人机验证，脚本无法匿名直取。提取页地址：
https://bgi.obs-mip.fr/data-products/outils/wgm2012-visualizationextraction/
改用的 GGMplus 与 WGM2012 的长波部分同样来自 EGM2008，并补了地形短波。
区别要说清楚：WGM2012 提供布格异常，GGMplus 只提供自由空气类的重力扰动，**地形信号没有扣除**。

**标注是什么**：geo3al 的每个多边形带岩性类型（TYPE）与地质时代（GEN_GLG）。在任何模型运行之前，按下面的规则落成 4 类：

| 类别 | 规则 | 可用格数 |
|---|---|---:|
| 新生代覆盖层 | 时代为 Q / N / Pg / T / Cz，以及风成沉积 | 79,713 |
| 中生界–古生界地层 | 时代为 K / J / Tr / P / C / D / S / O / Cm 及其组合 | 54,586 |
| 前寒武系 | 时代为 A / Pt / pCm | 13,221 |
| 岩浆岩 | 岩性类型为侵入 / 火山（任何时代） | 40,774 |

区域内 202,500 个格子里，4,704 个不在任何多边形内，7,449 个标注无法归类（水体、时代跨界如 PzpCm、未分等），
其余再去掉重磁任一缺值的格子，得到 188,294 个可用格。完整映射表在 `work/data/grid_summary.json`。

**特征**：两种数据源各 4 个完全对称的特征 —— 原值、梯度模、3×3 格局部标准差、减去 15×15 格（约 55 km）区域均值后的剩余场。

## 2. 怎么读这个结果

- **平衡准确率（BA）** = 4 类召回率的平均。随机猜或永远猜多数类都是 0.25。用它而不是普通准确率，是因为前寒武系只占约 7%。
- **空间自相关怎么处理的**：相邻格子的重磁特征和地质标注都很像。如果随机把格子分到训练集和测试集，
  测试格子的邻居几乎一定在训练集里，模型靠「抄邻居」就能得高分。主分析因此用 **1° × 1° 空间分块**（30 × 30 格，约 110 km）：
  整块随机分到 5 折里，同一块的格子永远不跨训练 / 测试。另报更保守的 2.5° 分块，以及明知有泄漏的像元随机划分，
  后者只用来量化泄漏有多大。
- **区间与 p 值**：ΔBA 在同一种子、同一折划分上配对计算。95% 区间按 1° 空间块（224 块）整块重抽 1000 次得到；
  p 值按块随机互换两个模型的预测 1000 次得到，最小可报告值约 0.001。
- **两种零模型**：① 标注在格子之间随机置换，切断一切对应关系，也抹掉空间结构 —— 这是「随便猜」的水平；
  ② 整张标注图相对特征图做环形平移（两个方向各至少 3°），**保留**两张图各自的空间结构、只切断对应关系 ——
  这是「只靠空间结构能蒙对多少」的水平，比 0.25 高，是更难越过的门槛。
  每个模型的环形平移零模型只有 9 次复本，经验 p 的下限就是 0.1，所以这里看的是 z 值与「是否高于全部复本的最大值」。
- **两把尺子**：把真实的重磁特征按已知规则造一份合成标注，看这套流程能不能测出已知的答案。正对照的标注同时依赖重力与磁
  （应当测出融合增益），负对照的标注只依赖重力（应当测不出）。任一把尺子失灵，P1 的结论按预注册作废。

## 3. 各模型在三种划分下的成绩

平衡准确率，括号内为种子 0 / 1 / 2：

| 模型 | 1° 空间分块（主分析） | 2.5° 空间分块 | 像元随机划分（有泄漏） |
|---|---|---|---|
| 多数类平凡基线 | 0.25 | 0.25 | 0.25 |
| 仅磁 | 0.2831（0.2844 / 0.2827 / 0.2821） | 0.2649 | 0.3477 |
| 仅重力 | 0.4533（0.4499 / 0.4544 / 0.4558） | 0.4302 | 0.4774 |
| 早期融合 | 0.4546（0.4504 / 0.4544 / 0.4590） | 0.4124 | 0.5517 |
| 后期融合 | 0.4486 | 0.4203 | 0.5043 |

1° 分块下的 macro-F1：仅磁 0.2664，仅重力 0.4009，早期融合 0.4183，后期融合 0.4041。
每折成绩在 `work/results/experiment.json`，每次拟合的日志在 `work/logs/01_experiment_main.log`。

**融合 − 仅重力（ΔBA，按种子）：**

| 划分 | 种子 0 | 种子 1 | 种子 2 | 均值 |
|---|---|---|---|---|
| 1° 分块 · 早期融合 | +0.0005 [−0.0098, +0.0107] p=0.944 | +0.0000 [−0.0110, +0.0097] p=0.995 | +0.0032 [−0.0068, +0.0124] p=0.604 | **+0.0012** |
| 1° 分块 · 后期融合 | −0.0036 [−0.0103, +0.0030] p=0.331 | −0.0052 [−0.0126, +0.0022] p=0.232 | −0.0055 [−0.0113, −0.0000] p=0.102 | −0.0048 |
| 2.5° 分块 · 早期融合 | −0.0165 [−0.0277, −0.0058] p=0.019 | −0.0190 [−0.0291, −0.0093] p=0.009 | −0.0182 [−0.0315, −0.0055] p=0.028 | −0.0179 |
| 2.5° 分块 · 后期融合 | −0.0084 [−0.0154, −0.0012] p=0.043 | −0.0107 [−0.0163, −0.0052] p=0.003 | −0.0106 [−0.0173, −0.0036] p=0.010 | −0.0099 |
| 随机划分 · 早期融合（事后补充） | +0.0750 [+0.0630, +0.0859] p=0.001 | +0.0741 [+0.0635, +0.0853] p=0.001 | +0.0739 [+0.0625, +0.0848] p=0.001 | +0.0743 |

图 1 与图 2（`figs/`）把这两张表画在一起。

## 4. 随机划分把成绩抬高了多少

| 模型 | 1° 分块 BA | 随机划分 BA | 抬高量 |
|---|---:|---:|---:|
| 仅重力 | 0.4533 | 0.4774 | +0.0240 |
| 仅磁 | 0.2831 | 0.3477 | +0.0646 |
| 早期融合 | 0.4546 | 0.5517 | **+0.0971** |
| 后期融合 | 0.4486 | 0.5043 | +0.0557 |

读法：

- 泄漏**不是**对所有模型一视同仁。单源 4 个特征的仅重力只被抬高 0.0240，8 个特征的早期融合被抬高 0.0971，
  于是随机划分下两者的差从 +0.0012 被放大到 +0.0743。**评估方式本身就能凭空造出一个「融合有效」的结论。**
- 预注册的 P2 写的是「至少 0.10」，实测 0.0971 落在 0.05–0.10 的「部分成立」档，按规则如实记为部分成立。
- 反过来，2.5° 分块比 1° 分块更严，早期融合从 0.4546 掉到 0.4124，掉得比仅重力（0.4533 → 0.4302）多。
  外推距离越远，多出来的磁特征越像是在拖后腿。

**分类别召回率（1° 分块，三种子平均；跑完之后加的描述性补充，不参与判定）：**

| 模型 | 新生代覆盖层 | 中生界–古生界地层 | 前寒武系 | 岩浆岩 |
|---|---:|---:|---:|---:|
| 仅重力 | 0.497 | 0.309 | 0.595 | 0.412 |
| 仅磁 | 0.351 | 0.208 | 0.232 | 0.342 |
| 早期融合 | 0.546 | 0.339 | 0.504 | 0.429 |

融合在覆盖层、地层、岩浆岩三类上略好，在前寒武系上差了约 0.09，四类一平均几乎抵消。这只是描述，没有做显著性检验。

## 5. 零模型

| 模型 | 1° 分块实测 BA | 环形平移零模型：均值 ± 标准差（9 次） | 零模型最大值 | 高出均值 | z |
|---|---:|---:|---:|---:|---:|
| 仅重力 | 0.4533 | 0.3457 ± 0.0258 | 0.3913 | +0.1076 | 4.17 |
| 早期融合 | 0.4546 | 0.3509 ± 0.0175 | 0.3789 | +0.1037 | 5.92 |
| 仅磁 | 0.2831 | 0.2774 ± 0.0118 | 0.2916 | +0.0057 | 0.48 |

标注置换零模型（早期融合，9 次）：BA 在 0.2485 到 0.2526 之间，与 0.25 无差别。

读法：保留空间结构的环形平移零模型能「蒙」到 0.35 左右，远高于 0.25 —— 这正是为什么只和随机猜比会高估信号。
仅重力与早期融合在每个种子上都高于全部零模型复本（P3 成立）；仅磁与它自己的零模型分不开（P4 成立）。

## 6. 这把尺子测得准吗（正负对照）

| 尺子 | 标注构造 | 仅重力 BA | 仅磁 BA | 早期融合 BA | ΔBA（融合 − 最强单源），按种子 |
|---|---|---|---|---|---|
| 正对照 | 2×[重力剩余场 > 中位数] + [磁剩余场 > 中位数]，25% 随机翻转 | 0.4467 / 0.4497 / 0.4489 | 0.4439 / 0.4477 / 0.4430 | 0.8093 / 0.8114 / 0.8099 | +0.3627 [+0.3571, +0.3684]；+0.3617 [+0.3558, +0.3673]；+0.3610 [+0.3556, +0.3669]；p 均为 0.001 |
| 负对照 | 重力剩余场四分位，25% 随机翻转 | 0.8096 / 0.8103 / 0.8098 | 0.2656 / 0.2641 / 0.2620 | 0.8096 / 0.8102 / 0.8097 | 0.0000、−0.0001、−0.0001；p = 0.504 / 0.109 / 0.195 |

两把尺子都按预注册的门槛通过：真有融合增益时，同一套模型、同一套 1° 分块能测出 +0.36；真没有时，测出来就是 0。
主结论 P1 因此有效。合成标注只用于这两个对照，没有混进任何真实结果。

## 7. 这一轮**不能**支持什么，以及已声明的混杂

1. **地形混杂（最要紧的一条）。** GGMplus 是地面重力扰动，未做布格改正。重力「看得准」的一部分可能是在看地形起伏：
   盆地低处多是新生代覆盖层，山区多出露前寒武系与岩浆岩。本轮没有引入高程数据，因此**不能**把 0.4533 这个数读成「重力识别岩性的能力」。
2. **分辨率不对等。** 重力原生约 220 m、求均值到 3.7 km；磁数据在这一带本来就是偏平滑的长波长编译。
   融合增益为零，可能是这份磁数据的分辨率问题，而不是融合方法本身的问题。
3. **标注的天花板。** 1:数百万比例尺的概化地质图，边界粗、时代分组粗，4 类划分也是我们定的。更细的标注可能改变结论。
4. **只有一种模型、一种区域、默认超参数。** 深度学习、图神经网络等更复杂的融合架构与超参数搜索都不在本轮范围内。
5. **不能**外推到矿区尺度、成矿预测或三维深部结构 —— 公开网格分辨率和公开标注都撑不到那一层。
6. **P1、P4 不是盲预测**（见第 0 节披露）。

如果要往下做，最能改变结论的三件事是：引入高程数据把地形从重力里拆出去；拿到 WGM2012 布格异常替换当前重力；
换一份分辨率更高的磁数据，看融合增益是否出现。

## 8. 怎么复现

```
python3 work/code/00_build_grid.py <原始数据目录> work/data/grid.npz work/data/grid_summary.json
python3 work/code/01_experiment.py work/data/grid.npz main   work/results
python3 work/code/01_experiment.py work/data/grid.npz nulls  work/results
python3 work/code/01_experiment.py work/data/grid.npz rulers work/results
python3 work/code/02_stats.py work/results
python3 work/code/02b_supplement.py work/results
python3 work/code/03_recheck.py work/results
python3 work/code/04_figs.py work/results work/data/grid.npz figs
python3 audit_delivery.py
```

- 环境：Python 3.11，numpy 2.5，scikit-learn 1.9.1，matplotlib 3.11；统计与复算步骤也在 Python 3.14 / numpy 2.4 下跑过，结果一致。
- `01_experiment.py` 的三个部分互不依赖、并行运行，每次拟合都打一行日志（`work/logs/`）。
- `02_stats.py` 在导出逐格预测表那一段做过一次纯提速的改动（原写法对每个格子重新解压一次整个数组），
  统计部分未动；改动前后生成的 `experiment.json` 逐字节一致。
- `02b_supplement.py`（随机划分下的区间与 p、分类别召回率）与 `04_figs.py` 是跑完之后加的，不参与 P1–P4 判定。
- `03_recheck.py` 只用 Python 标准库，从逐格预测重新数混淆矩阵；各臂 BA、ΔBA 与抬高量和主分析的差异不超过 1e-15。
- `audit_delivery.py` 同样只用标准库：检查预注册指纹、重算头条数字、复核 P1–P4 与尺子判定。