看懂火山 EnKF｜中文互动教材

Chapter 00

先别管公式：它是在做“地下侦探”

看不见地下岩浆房，就先提出很多种可能，再用地表观测一步步淘汰不合理的猜测。EnKF 的核心其实就是这么朴素。

整份代码只有一条主线：猜测地下状态 → 预测地表变形 → 和真实观测比较 → 修正猜测 → 第二天继续。

01地下状态

位置、深度、半径、压力，共 5 个参数。

02Mogi 模型

把地下参数翻译成地表会移动多少。

03GPS / InSAR

仪器告诉我们地表实际上移动了多少。

04EnKF 修正

预测和观测差得越多，猜测就需要改得越多。

05连续 100 天

每天吸收新证据，形成一条参数演化轨迹。

为什么不是只猜一个答案？

学生

直接算出“岩浆源在这里、压力这么大”不就行了吗？

老师

观测有噪声，模型也不完美，而且不同参数可能产生相似的地表变形。一个答案会假装我们“非常确定”；100 个答案可以同时表达最可能值和不确定性。

学生

所以每个集合成员就是一个地下世界？

老师

对。每个成员都是一套完整的地下源参数。把 100 个成员都送进 Mogi 模型，就得到 100 套地表变形预测。

先记住四个词

A

状态 state

我们真正想知道、但不能直接测到的地下参数。本教程的状态有 5 维。

例： \(\boldsymbol\theta=[x_0,y_0,z_0,a,\Delta P]^\mathsf T\)

B

观测 observation

仪器真正给出的数据。GPS 是三方向位移，InSAR 是雷达视线方向的位移。

C

前向模型 forward model

输入地下参数，输出地表预测。本教程使用 Mogi 模型。

D

同化 assimilation

把模型预测和观测合并，用新证据修正状态估计。

Chapter 01

只补真正会用到的数学

不需要先学完整的线性代数和概率论。下面四件事够你读懂大部分公式：向量、均值与方差、协方差、正态分布。

01 / VECTOR

向量：把一组数打包

一个状态向量不是神秘对象，只是按固定顺序放在一起的 5 个数字。顺序一旦改变，含义就会改变。

\(\boldsymbol\theta=[x_0,y_0,z_0,a,\Delta P]^\mathsf T\)

02 / MATRIX

矩阵：把很多向量并排

如果每一列是一套地下参数，100 列就是 100 个可能世界。这里矩阵大小是 \(5\times100\)。

行 = 参数；列 = 集合成员。

03 / VARIANCE

方差：大家分得有多散

100 个成员越分散，说明我们越不确定；越集中，说明估计越明确。标准差是方差开平方，更接近原始单位。

04 / COVARIANCE

协方差：两个量会不会一起变

如果“压力高”的成员通常也预测“地表抬升大”，压力与抬升就有正协方差。EnKF 正是利用这种共同变化来修正参数。

协方差为什么是 EnKF 的发动机？

假设观测告诉我们“真实抬升比预测更大”。EnKF 会查看 100 个成员：哪些成员能产生更大的抬升？如果这些成员往往有更高的压力，它就会提高压力；如果它们往往更浅，它也可能减小深度。

协方差不是告诉模型“正确答案是什么”，而是告诉模型“应该沿哪个方向修改参数”。

正态分布在这里做什么？

记号 \(X\sim\mathcal N(\mu,\sigma^2)\) 表示：数值大多落在均值 \(\mu\) 附近，离均值越远越少见，分散程度由标准差 \(\sigma\) 决定。

直觉

\[x_0\sim\mathcal N(0,1^2)\]

源的东西位置以 0 km 为中心，标准差 1 km

不是说一定在 0。而是说大多数初始成员在 0 附近，少部分会更远。

Chapter 02

仪器到底测到了什么？

GPS 和 InSAR 都在测地表变形，但它们“看世界”的方式完全不同。一个连续、稀疏、三维；另一个间隔较长、密集、只有一个投影方向。

比较	GPS	InSAR
本教程频率	每天一次	每 10 天一次
空间分布	9 个离散站点	133 个像元，覆盖范围更密
每个位置给多少数	东、北、上 3 个分量	沿雷达视线 LOS 的 1 个分量
优势	时间连续、方向完整	空间覆盖广，容易看出变形形状
指定噪声	水平 2.5 mm；垂向 5 mm	LOS 8 mm

LOS：把三维位移“压扁”到一条视线

雷达不是分别告诉我们东、北、上移动多少，而是看位移在雷达视线方向上投影了多少。就像手电筒照一个三维物体，墙上只留下一个方向的影子。

N24–N25

\[\widehat{\boldsymbol\ell}=[\ell_E,\ell_N,\ell_U]^\mathsf T,\qquad u_{\mathrm{LOS}}=\ell_Eu_E+\ell_Nu_N+\ell_Uu_U\]

三维位移与 LOS 单位向量做点积

点积的通俗解释：分别计算东、北、上三个方向对雷达视线贡献了多少，再把贡献相加。

为什么两种数据要一起用？

GPS 像每天写一次日记，告诉你变化的时间过程；InSAR 像每十天拍一张广角照片，告诉你变形的空间形状。前者补时间，后者补空间。

同化的价值不只是“数据更多”，而是把不同数据各自擅长的部分放在同一个物理模型中。

Chapter 03

从地下球体到地表位移

Mogi 模型把岩浆房简化成弹性半空间里的小球形压力源。它不是真实火山的全部，却能用很少参数抓住“增压导致地表向外、向上鼓起”的主要形状。

调参数，观察地表怎么变

源深度 depth3.5 km

源半径 radius0.8 km

压力变化 pressure60 MPa

源强度 \(a^3\Delta P\)30.72

中心最大抬升—

剪切模量 \(\mu\)20 GPa

拖动画面可以旋转。箭头为了看清被大幅放大，数值读数仍按真实 Mogi 公式计算。

Mogi 公式按因果顺序拆开

N15 → N21

N15

\[\mu=\frac{E}{2(1+\nu)}\]

先把杨氏模量换成剪切模量

为什么需要它：Mogi 位移公式用 \(\mu\) 描述岩石抵抗剪切变形的能力。\(E\) 越大，岩石越硬，同样压力造成的位移越小。

N16

\[\rho_i=\sqrt{(x_i-x_0)^2+(y_i-y_0)^2}\]

计算观测点离源正上方有多远

这是平面上的勾股定理。先算东西差，再算南北差。\(\rho_i=0\) 表示观测点正好位于源中心上方。

N17

\[B=\frac{(1-\nu)a^3\Delta P}{\mu}\]

把半径、压力和岩石弹性合成“变形强度”

最关键的是 \(a^3\Delta P\)。半径扩大一点会被三次方放大，所以半径与压力很容易互相补偿。

N18

\[u_{E,i}=\frac{B(x_i-x_0)}{(z_0^2+\rho_i^2)^{3/2}},\quad u_{N,i}=\frac{B(y_i-y_0)}{(z_0^2+\rho_i^2)^{3/2}},\quad u_{U,i}=\frac{Bz_0}{(z_0^2+\rho_i^2)^{3/2}}\]

把强度分配到东、北、上三个方向

分子决定方向，分母决定衰减。离源越远，三维距离 \(\sqrt{z_0^2+\rho_i^2}\) 越大，位移会很快减小。

分三步理解分母为什么是三次方

源到观测点的三维距离是 \(R=\sqrt{z_0^2+\rho_i^2}\)。
点源引起的场强随距离大致按 \(1/R^2\) 衰减。
再乘方向余弦，例如垂向分量为 \(z_0/R\)，得到 \(z_0/R^3\)。

N19

\[V_0=\frac{4}{3}\pi a^3\]

球形源的参考体积

这是普通球体体积公式。它再次说明为什么半径以三次方进入问题。

N20

\[\Delta V=\frac{V_0}{K}\Delta P=\frac{4\pi}{3K}a^3\Delta P\]

压力变化对应多少源体积变化

体积变化 \(\Delta V\) 仍然依赖 \(a^3\Delta P\)，所以它往往比“半径”和“压力”两个单独参数更稳定。

N21

\[K=\frac{E}{3(1-2\nu)}\]

体积模量：岩石抵抗整体压缩的能力

\(K\) 越大，岩石越难被压缩；相同压力对应的体积变化越小。

Mogi 模型的假设与边界

假设	通俗含义	可能带来的偏差
球形、小源	岩浆房用一个小球代替	真实岩床、岩脉或不规则源可能不适合
均匀、各向同性	四周岩石性质一样	分层、断层、不同岩性会改变变形形状
线性弹性	压力加倍，位移大致加倍	塑性、破裂和黏弹性行为无法表示
平坦半空间	地表当作无限平面	陡峭火山地形会产生误差

Chapter 04

代码怎样“造出”100 天观测？

合成实验的好处是我们知道真正答案。先用真实参数算无噪声位移，再加入已知噪声，最后看 EnKF 能不能把答案找回来。

从真实压力到观测向量

N22 → N25

N22

\[\Delta P_k^{\mathrm{true}}=\dot P_{\mathrm{true}}t_k,\qquad \dot P_{\mathrm{true}}=1\ \mathrm{MPa/day}\]

真实压力每天增加 1 MPa

第 1 天是 1 MPa，第 100 天是 100 MPa。几何参数保持不变。

N23

\[\mathbf y_k^{\mathrm{GPS}}=[u_{E,1},u_{N,1},u_{U,1},u_{E,2},u_{N,2},u_{U,2},\ldots]^\mathsf T\]

把 9 个 GPS 站的三分量排成 27 维向量

顺序不能错。每个站连续放东、北、上三个数，然后才轮到下一个站。

N24

\[\widehat{\boldsymbol\ell}=[\ell_E,\ell_N,\ell_U]^\mathsf T,\qquad \ell_E^2+\ell_N^2+\ell_U^2=1\]

定义雷达看向地面的单位方向

本教程使用约 \([-0.38,0.08,0.92]\)，即 LOS 对垂向最敏感。

N25

\[u_{\mathrm{LOS}}=\ell_Eu_E+\ell_Nu_N+\ell_Uu_U\]

把三维位移投影为一个 LOS 数

每个 InSAR 像元只贡献一个观测值，但像元数量多，能提供空间形状。

一天的数据是怎样拼起来的？

10GPS + InSAR

20GPS + InSAR

30GPS + InSAR

40GPS + InSAR

50GPS + InSAR

60GPS + InSAR

70GPS + InSAR

80GPS + InSAR

90GPS + InSAR

100GPS + InSAR

普通日：\(m_k=27\)。InSAR 日：\(m_k=27+133=160\)。EnKF 的状态维度始终是 5，但观测维度会随日期变化。

为什么一定要加噪声？

如果把完美的无噪声数据直接交给滤波器，问题会比真实世界容易得多。加入已知标准差的高斯噪声，是在模拟 GPS 和 InSAR 仪器的不确定性，也让我们能检验 nRMSE 是否接近 1。

Chapter 05

用 100 个成员表达“不确定”

集合不是重复计算 100 次而已。它用成员之间的分散程度表达不确定性，用成员之间的共同变化估计协方差。

初始集合与物理边界

N26 → N28

N26

\[\theta_{p,j}^{\mathrm{initial}}\sim\mathcal N(\mu_{p,0},\sigma_{p,0}^2)\]

每个参数从指定分布抽样

\(p\) 表示参数，\(j\) 表示成员。每个成员都独立抽一套参数。

N27

\[\mathbf A_{\mathrm{initial}}=[\boldsymbol\theta_1,\boldsymbol\theta_2,\ldots,\boldsymbol\theta_N]\in\mathbb R^{5\times N}\]

把 100 个状态并排组成矩阵

矩阵第 \(j\) 列就是第 \(j\) 个可能世界。

N28

\[-8\le x_0,y_0\le8,\quad1\le z_0\le8,\quad0.2\le a\le1.3,\quad0\le\Delta P\le160\]

给参数设置物理允许范围

越界成员不是简单截断，而是像小球碰墙一样“反射”回来，避免成员全堆在边界。

初始分布长什么样？

参数	代码中的初始分布	真实值（第 1 天）	直觉
\(x_0\)	\(\mathcal N(0,1^2)\) km	0.25 km	横向位置先猜在火山中心附近
\(y_0\)	\(\mathcal N(0,1^2)\) km	−0.50 km	允许约公里级偏移
\(z_0\)	\(\mathcal N(4,0.65^2)\) km	3.50 km	故意与真实值略错开
\(a\)	\(\mathcal N(1,0.12^2)\) km	0.80 km	用较有信息的先验限制半径
\(\Delta P\)	\(\mathcal N(5,2^2)\) MPa	1 MPa	代码实际使用正态分布

原文字与代码有一处不一致

说明文字写“初始压力从 0–2 MPa 的均匀分布抽样”，但真正执行的代码使用 \(\mathcal N(5,2^2)\)，均匀分布那一行被注释掉了。理解程序时，最终应以实际执行的代码为准。

Chapter 06

EnKF 的 14 步数学地图

不要把 N1–N14 当成 14 个孤立公式。它们分成四件事：组织猜测、组织观测、比较预测与观测、根据差异更新状态。

A. 描述 100 个预测状态

N1–N3

N1

\[\mathbf A_k^f=[\boldsymbol\theta_{1,k}^f,\ldots,\boldsymbol\theta_{N,k}^f]\in\mathbb R^{n\times N}\]

预测集合矩阵

一句话：把第 \(k\) 天的 100 套预测参数并排放好。这里 \(n=5,N=100\)。

N2

\[\overline{\mathbf A}_k^f=\mathbf A_k^f\mathbf1_N,\qquad\mathbf A_k^{\prime f}=\mathbf A_k^f-\overline{\mathbf A}_k^f\]

集合平均与成员偏差

先求大家的平均答案，再记录每个成员偏离平均多少。撇号 \('\) 就表示这种偏差。

N3

\[\mathbf P_{e,k}^f=\frac{\mathbf A_k^{\prime f}(\mathbf A_k^{\prime f})^\mathsf T}{N-1}\]

从成员偏差估计协方差

矩阵对角线是每个参数的方差；非对角线是参数之间的共同变化。除以 \(N-1\) 是样本协方差的标准做法。

B. 描述观测与观测噪声

N4–N7

N4

\[\mathbf d_k=[d_{1,k},d_{2,k},\ldots,d_{m_k,k}]^\mathsf T\]

当天真正测到的观测向量

普通日有 27 个数，InSAR 日有 160 个数。

N5

\[\mathbf R_k=\operatorname{diag}(\sigma_{1,k}^2,\ldots,\sigma_{m_k,k}^2)\]

观测误差协方差

这里假设各观测误差相互独立，所以只有对角线有值。\(\sigma\) 越小，观测越可信。

N6

\[\mathbf d_{j,k}=\mathbf d_k+\boldsymbol\epsilon_{j,k},\qquad\boldsymbol\epsilon_{j,k}\sim\mathcal N(\mathbf0,\mathbf R_k)\]

每个成员配一份略有不同的观测

随机 EnKF 不让所有成员看到完全一样的数据，而是在合理误差范围内扰动观测。

N7

\[\boldsymbol\Upsilon_k=[\boldsymbol\epsilon_{1,k},\ldots,\boldsymbol\epsilon_{N,k}],\qquad\mathbf D_k=\mathbf d_k\boldsymbol1_N^\mathsf T+\boldsymbol\Upsilon_k\]

把 100 份扰动观测组成矩阵

\(\boldsymbol\Upsilon_k\) 的每一行还会减去自身平均值，确保扰动整体均值为零，不把观测系统性推高或推低。

C. 把状态变成预测观测

N8–N9

N8

\[\mathbf Y_k^f=[h_k(\boldsymbol\theta_{1,k}^f),\ldots,h_k(\boldsymbol\theta_{N,k}^f)]\]

对每个成员运行非线性 Mogi 前向模型

这里写成 \(\mathbf H_k\mathbf A_k^f\) 只是传统记号，代码并没有拿一个固定 \(\mathbf H\) 直接乘，因为 Mogi 是非线性的。

N9

\[\mathbf Y_k'=\mathbf Y_k^f-\overline{\mathbf Y_k^f},\qquad\mathbf D_k'=\mathbf D_k-\mathbf Y_k^f\]

预测观测偏差与创新

创新 innovation 就是“观测 − 预测”。它告诉滤波器每个成员错在哪里、错多少。

D. 求更新量、决定是否保留

N10–N14

N10

\[[(N-1)\mathbf I_N+(\mathbf Y_k')^\mathsf T\mathbf R_k^{-1}\mathbf Y_k']\mathbf W_k=(\mathbf Y_k')^\mathsf T\mathbf R_k^{-1}\mathbf D_k'\]

在集合空间里求分析权重

左边：原集合稳定性 + 预测观测相对于噪声的变化。
右边：预测变化与观测创新的对应关系。
解出的 \(\mathbf W_k\) 告诉我们该怎样组合状态偏差。

为什么转到集合空间求解？

直接在观测空间求解，需要处理 \(m_k\times m_k\) 矩阵；InSAR 日这里是 160×160。
集合空间只处理 \(N\times N\) 矩阵，这里是 100×100。
当观测数量大于集合数时，这样通常更省计算。

N11

\[\mathbf A_k^a=\mathbf A_k^f+\mathbf A_k^{\prime f}\mathbf W_k\]

预测集合 + 观测给出的修正 = 分析集合

这是整套 EnKF 最关键的更新式。\(\mathbf A'\mathbf W\) 把观测信息映射回位置、深度、半径和压力。

N12

\[\boldsymbol\theta_{j,k}^{f}=\boldsymbol\theta_{j,k-1}^{\mathrm{selected}}\]

把前一天保留的集合原样带到下一天

这叫持续性预测 persistence forecast。本代码没有在预测步骤主动增加压力，也没有加过程噪声。

N13

\[\operatorname{nRMSE}_k=\sqrt{\frac1{m_k}\sum_{i=1}^{m_k}\left(\frac{\overline y_{i,k}-d_{i,k}}{\sigma_{i,k}}\right)^2}\]

把残差除以各自噪声后再算均方根

nRMSE 约为 1，表示模型误差大致和观测噪声一样大；远大于 1，表示拟合还不够好。

N14

\[\mathbf A_k^{\mathrm{selected}}=\begin{cases}\mathbf A_k^a,&\mathrm{nRMSE}_k^a\le1\ \text{或}\ \mathrm{nRMSE}_k^a<(1-\delta)\mathrm{nRMSE}_k^f\\\mathbf A_k^f,&\text{否则}\end{cases}\]

RMSE 门控：分析真的更好才保留

标准 EnKF 总是保留分析集合；门控版先检查拟合质量。原代码还有一条：如果预测本来已经 \(\le1\)，就直接保留预测。

Chapter 07

把公式一行行翻译成 Python

下面保留最关键的原代码结构。你可以切换模块、复制代码，再对照每段下方的中文解释。

def mogi_displacement(x, y, x0, y0, z0, a, dP, E, nu):
    # 1. 岩石弹性：公式 N15
    mu = E / (2.0 * (1.0 + nu))

    # 2. 观测点相对源中心的位置：公式 N16
    dx = x - x0
    dy = y - y0
    rho2 = dx**2 + dy**2

    # 3. 源的综合变形强度：公式 N17
    B = (1.0 - nu) * a**3 * dP / mu

    # 4. 距离衰减，并返回东、北、上位移：公式 N18
    denominator = (z0**2 + rho2)**1.5
    return B*dx/denominator, B*dy/denominator, B*z0/denominator

输入是什么？

观测点坐标、源位置/深度/半径/压力，以及岩石弹性参数。

输出是什么？

每个观测点的东、北、向上位移，单位为米。

def draw_observation_perturbations(sigma, N, rng):
    sigma = np.asarray(sigma, float)

    # 每个观测、每个成员各抽一次随机误差
    Upsilon = rng.normal(
        0.0,
        sigma[:, None],
        size=(sigma.size, N)
    )

    # 每一行减去集合均值，让扰动整体均值为零
    return Upsilon - Upsilon.mean(axis=1, keepdims=True)

为什么要中心化？

避免随机扰动恰好整体偏正或偏负，给分析带来额外系统偏差。

矩阵形状

\(m_k\times N\)：每一列对应一个集合成员的观测扰动。

def stochastic_enkf_analysis(A_f, d, sigma, forward_ensemble, rng):
    n, N = A_f.shape

    # N8：让每个预测成员通过非线性前向模型
    Y_f = forward_ensemble(A_f)

    # N2、N9：状态偏差和预测观测偏差
    A_prime = A_f - A_f.mean(axis=1, keepdims=True)
    Y_prime = Y_f - Y_f.mean(axis=1, keepdims=True)

    # N6–N9：构造扰动观测与创新
    Upsilon = draw_observation_perturbations(sigma, N, rng)
    D_prime = d[:, None] + Upsilon - Y_f

    # N10：在集合空间求权重 W
    Rinv_Y = Y_prime / (sigma[:, None]**2)
    Rinv_D = D_prime / (sigma[:, None]**2)
    system = (N-1)*np.eye(N) + Y_prime.T @ Rinv_Y
    rhs = Y_prime.T @ Rinv_D
    W = np.linalg.solve(system, rhs)

    # N11：把观测修正映射回物理状态
    A_a = A_f + A_prime @ W
    return A_a

@ 是什么？

Python 的矩阵乘法。它不是逐元素相乘。

solve 而不是求逆

直接解线性方程通常比显式计算矩阵逆更稳定、更高效。

if not gated:
    # 标准 EnKF：总是接受候选分析
    accept = True

elif forecast_nrmse <= target_nrmse:
    # 预测已经达到噪声水平，不再改它
    accept = False

else:
    # 候选达到目标，或至少改善 min_improve
    accept = (
        candidate_nrmse <= target_nrmse
        or candidate_nrmse
           < forecast_nrmse * (1.0 - min_improve)
    )

A = A_candidate if accept else A_f

门控改了什么？

只改“保留分析还是预测”的选择，不改前面的卡尔曼分析公式。

本教程阈值

target_nrmse=1.0，门控运行时 min_improve=0.002，即至少改善 0.2%。

for day in DAYS:
    obs = obs_by_day[int(day)]
    forward = daily_forward_operator(obs["include_insar"])

    # 1. 预测：把昨天保留的集合原样带过来
    A_f, _ = forecast_ensemble(A)

    # 2. 计算预测拟合
    forecast_metrics = observation_metrics(
        A_f, obs["d"], obs["sigma"], forward
    )

    # 3. 用当天观测做一次候选分析
    A_raw, _ = stochastic_enkf_analysis(
        A_f, obs["d"], obs["sigma"], forward, rng
    )
    A_candidate, repairs = safeguard_physical_ensemble(A_raw)

    # 4. 标准版总接受；门控版按 N14 决定
    A = A_candidate if accept else A_f

    # 5. 保存均值、标准差、预测和诊断量
    out["mean"].append(A.mean(axis=1))
    out["std"].append(A.std(axis=1, ddof=1))

每天只分析一次

同一天的数据不会反复塞进滤波器，避免重复使用同一份证据。

InSAR 日的区别

只是当天观测向量更长；分析函数本身不需要换一套。

Chapter 08

标准 EnKF 与 RMSE 门控

两个工作流使用相同的观测、随机种子、Mogi 模型和卡尔曼分析。差别只有一个：候选分析是否一定被保留。

步骤	标准随机 EnKF	RMSE 门控 EnKF
预测	前一天保留集合原样带到今天
候选分析	用同一套 N8–N11 计算一次
是否保留	总是保留候选分析	达到噪声目标或有足够改善才保留
优点	简单、连续吸收所有信息	避免一次不好的分析把状态推向更差位置
风险	坏分析也会被带到下一天	门控可能拒绝某些有用但改善不大的更新

互动：你来判断这一天要不要接受

预测 nRMSE

候选分析 nRMSE

最低改善比例 \(\delta\)

修改数值，右侧会立即按原代码逻辑重新判断。

接受分析

候选分析更好

候选 nRMSE 相比预测有足够改善，因此保留分析集合。

最低需要小于1.397

互动：看 100 个成员怎样逐日收敛

下图画的是“半径—压力”平面。绿色曲线是保持真实源强度 \(S=a^3\Delta P\) 不变的组合；很多不同半径和压力都能落在相近曲线上，这就是不可辨识性。

集合成员演化 conceptual animation

同化日 Day 0

圆点是集合成员，星号是真实值。动画用于解释收敛与参数权衡，不是重新运行完整 Python 模型。

重要：实际预测函数没有“每天加 1 MPa”

真实合成源的压力按 N22 每天增加，但滤波器的 forecast_ensemble(A) 只是 return A.copy(), 0。也就是说，滤波器不知道真实压力增长规律，只能靠每天的新观测把压力重新推高。

Chapter 09

结果不是“一个正确答案”

要同时看集合均值、集合分散、观测拟合、参数边界和物理可辨识性。某个数字接近真实值，并不代表所有参数都独立可靠。

本机从头运行得到的第 100 天结果

参数	真实值	标准 EnKF 均值	RMSE 门控均值
\(x_0\) (km)	0.250	0.294	0.294
\(y_0\) (km)	−0.500	−0.463	−0.427
深度 \(z_0\) (km)	3.500	3.365	3.338
半径 \(a\) (km)	0.800	1.294	1.296
压力 \(\Delta P\) (MPa)	100.0	19.96	19.01

位置和深度相对接近真实值，但半径偏大、压力偏小。这不是简单的“模型失败”，而是 Mogi 模型主要看见 \(a^3\Delta P\) 的典型表现。

为什么半径和压力会互相替代？

核心

\[S=a^3\Delta P\quad\Longrightarrow\quad\Delta P=\frac{S}{a^3}\]

保持同样的 \(S\)，大半径可以配小压力

如果半径从 0.8 km 增加到 1.2 km，为保持相同源强度，压力只需要原来的 \((0.8/1.2)^3\approx0.30\)。

实际运行还暴露了两处值得核对的地方

标准 EnKF 的分析边界保护累计触发 1470 次，门控版触发 1678 次；因此“边界保护保持不活动”的文字结论与实际运行不一致。
半径上界是 1.3 km，而最终均值约 1.295 km，非常靠近上界。解释压力和半径时必须明确说明边界正在影响结果。

一份合格的科学解读应该怎么写？

DO

可以说

数据较好约束源的水平位置、深度和综合源强度；两种工作流得到相近位置。

DON'T

不要直接说

岩浆房半径就是 1.295 km、压力就是 19.96 MPa，而且两者都同样可靠。

CHECK

必须检查

nRMSE、残差空间分布、集合范围、边界触发次数，以及半径—压力散点图。

加入更真实的预测模型、过程噪声、非球形源、地形/分层弹性或独立地质先验。

Chapter 10

确认自己真的理解了

先自己回答，再点开解释。如果能用自己的话讲清这些问题，你已经能读懂原 Notebook 的主干。

因为成员群体同时表达最可能状态和不确定性，而且成员之间的共同变化提供了状态—观测协方差，告诉滤波器该怎样修正参数。

地表变形振幅主要依赖 \(a^3\Delta P\)。许多“较大半径 + 较小压力”和“较小半径 + 较大压力”的组合会产生近似变形。

GPS 每天提供稀疏但完整的三分量时间序列；InSAR 每 10 天提供密集的 LOS 空间图。一个补时间连续性，一个补空间形状。

模型平均残差大致与规定的观测噪声同量级。它不等于模型“完全正确”，只表示误差没有明显超过噪声预期。

没有。N8–N11 的候选分析完全相同；门控只在分析之后决定保留候选分析还是原预测。

不会。真实合成数据按每天 1 MPa 增长，但 forecast_ensemble 只复制集合。滤波器必须从新观测中推断压力增长。

中英术语速查

中文	英文	在本教程中的意思
集合	ensemble	100 套可能的地下源状态
预测	forecast	吸收当天观测前的状态/观测预测
分析	analysis	吸收当天观测后的更新状态
扰动 / 异常量	perturbation / anomaly	成员偏离集合平均的部分
创新	innovation	扰动观测减去对应预测
协方差	covariance	两个量在集合中共同变化的程度
前向算子	forward operator	把源参数变成 GPS/InSAR 预测的函数
视线向	line of sight, LOS	雷达观测敏感的方向
均方根误差	root mean square error, RMSE	综合衡量预测与观测差异
可辨识性	identifiability	数据能否把不同参数独立区分出来

读原 Notebook 的推荐顺序

01先跑一遍

先看图和输出，不急着读每个矩阵。

02读 Mogi

确认 5 个参数怎样产生位移。

03读 N8–N11

抓住“预测观测 → 创新 → 权重 → 更新”。

04读主循环

看每天的预测、分析、选择和保存。

05查结果边界

不要只看最终均值，要看不确定性和保护次数。

看懂火山EnKF