简介：GP-EnKF资源围绕在线高斯过程回归与集合卡尔曼滤波的融合实现，面向需要处理动态数据流、实时状态估计与不确定性量化的机器学习学习者及工程人员，重点解决传统GPR在在线数据更新中计算成本高、难以实时适应新观测的问题。内容结合Fusion 2018论文背景，涉及高斯过程先验、集合近似后验、预测—更新迭代及Python实现思路，可用于环境监测、控制工程、信号处理等场景的代码复现与实验参考。压缩包约22KB，文件总数与类型明细暂未提供（显示为0/暂无数据），使用前需自行解压核对目录与依赖。已有307人学习，适合具备Python与概率基础、希望理解GP-EnKF算法流程并借鉴实现结构的读者。

1. GP-EnKF 是什么：用一组粒子权重替高斯过程跑在线回归

产线上一条温度曲线每秒吐一个点，你想让回归模型跟着数据走，而不是每来一个观测就把整个高斯过程重训一遍——精确 GP 每步 O(n³)，跑半天就卡死。GP-EnKF 正是为这个场景准备的：它用一组集成成员去表示诱导点上的权重，观测到来时只做一次集合卡尔曼滤波的递推更新，全程只有矩阵乘法，Python 加上 NumPy 就能跑起来。这套代码来自一篇在线高斯过程学习的论文实现，适合做流式回归、传感器标定、设备退化趋势跟踪这类「数据一条条来、模型要一直在线」的活。你要是只会批量fit/predict，读它能补上在线贝叶斯这一课。

2. 模型骨架：诱导点、权重状态空间与 EnKF 的两步递推

先把整体结构说清楚，不然后面看代码会一头雾水。这套实现的核心思路只有一句话：别去更新函数值本身，去更新一组诱导点上的权重。函数值有无穷多个，权重只有 m 个（m 是诱导点数量，通常几十到几百），状态维度一下就固定住了。

2.1 为什么不用精确 GP：诱导点把 O(n³) 压成 O(m³)

精确 GP 回归的后验是

mu(x*) = k*ᵀ (K + σ²I)⁻¹ y

每来一个新观测，K 就长大一行一列，求逆是 O(n³)。在线场景里 n 一直涨，这条路走不通。

稀疏 GP 的做法是引入 m 个伪输入（诱导点）Z = [z₁, …, z_m]，假设函数值可以由这些点上的权重线性组合近似：

f(x) ≈ φ(x)ᵀ w, φ(x) = k(x, Z)

φ(x) 就是 x 和所有诱导点的核相似度向量，长度固定为 m。这样一来，未知量从「无穷维函数」变成「m 维权重向量 w」，复杂度从 O(n³) 降到 O(m³)，而且和已观测样本数彻底解耦。我一般把 m 取在 30～200 之间，看输入维度和平滑程度，太少了欠拟合，太多了单步更新又开始变慢。

提示：诱导点的位置不是随便撒的。均匀铺在训练区间上是最省事的做法，想要更好效果就用 k-means 在已有数据上聚一聚再当初始 Z，后面靠过程噪声去吸收位置偏差。

2.2 状态空间写法：权重是状态，观测是函数值

有了 φ(x)，模型就变成了一个标准的线性状态空间：

• 状态方程：w_t = w_{t-1} + η，η ~ N(0, Q)，随机游走，允许模型随时间缓慢漂移
• 观测方程：y_t = φ(x_t)ᵀ w_t + ε，ε ~ N(0, σ_n²)

如果直接用卡尔曼滤波，需要维护一个 m×m 的协方差矩阵 P，每步做 P ← P − K H P，m=200 时就是 4 万个数，还要做矩阵求逆，Python 里跑几千步就明显吃力。

EnKF 换了个思路：不维护协方差，维护一个样本集。用 M 个成员 {w⁽¹⁾, …, w⁽ᴹ⁾} 去近似后验分布，均值就是样本均值，协方差就是样本协方差：

P ≈ (1/(M-1)) Σᵢ (w⁽ⁱ⁾ − w̄)(w⁽ⁱ⁾ − w̄)ᵀ

每步只需要 O(M·m) 的矩阵乘法，而且天然并行，m 再大也不怕。代价是采样误差——M 太小协方差估不准，这就需要后面讲的「观测扰动」和「协方差膨胀」来补救。

2.3 EnKF 与粒子滤波、递归最小二乘的差别

同样是序贯更新，这三者的取舍值得说清楚，选错了会白写一堆代码。

GP-EnKF 的位置很清楚：状态维度 m 可能有上百维，粒子滤波在这个维度上直接失效（权重集中到一个粒子上），而 RLS 给不出预测方差、也没法在线学超参数。EnKF 用「高斯假设 + 样本协方差」换来了高维下的稳定性，代价是后验只能是高斯的——对 GP 回归来说这刚好不是什么损失。

顺带一提，EnKF 的更新式和卡尔曼滤波长得一模一样，只是把理论协方差换成了样本协方差：

K = C_wy (C_yy + σ_n²)⁻¹ w⁽ⁱ⁾ ← w⁽ⁱ⁾ + K (y + ε⁽ⁱ⁾ − φ(x)ᵀ w⁽ⁱ⁾)

看懂这个式子，整个项目就通了七成。

3. 核心代码拆解：一步 predict 加一步 update 的完整链路

这一章把关键函数一段段拆开。代码是 NumPy 风格，和论文里的符号一一对应，你可以直接抄进自己的工程里改。

3.1 核矩阵与诱导点特征

所有东西都建立在核函数上。RBF 是最常用的，两个超参数：长度尺度ls控制平滑程度，信号方差sf控制幅值。

import numpy as np def rbf(A, B, ls=1.0, sf=1.0): """A: (n, d), B: (m, d) -> (n, m) 的 RBF 核矩阵""" d2 = (A[:, None, :] - B[None, :, :]) ** 2 # 广播出 (n, m, d) 的平方差 return sf ** 2 * np.exp(-0.5 * d2.sum(-1) / ls ** 2)

这里刻意用了显式广播而不是cdist，因为后面要按成员各自的ls/sf分别算，写成显式公式改起来方便。d2.sum(-1)把输入维度 d 消掉，留下样本对之间的平方距离矩阵，再套指数。

诱导点特征矩阵就是一次核计算：

m = 30 Z = np.linspace(-5, 5, m)[:, None] # 诱导点位置，单维输入，形状 (m, 1) Phi = rbf(X_train, Z, ls=1.2, sf=1.0) # 训练集的 (n, m) 特征矩阵

Phi的每一行就是 φ(xᵢ)，长 m。如果后面要做批量回归，这个矩阵算一次就能复用，别在循环里反复算——我见过有人把它放进 for 循环，2000 个点跑了三分半。

3.2 集成初始化与过程噪声

集成成员从先验里采样。GP 先验下权重服从 N(0, K_zz)，工程上直接按标准正态撒一遍也够用，反正几十步之后就被数据纠正了。

rng = np.random.default_rng(0) M = 100 # 集成成员数，50~200 之间比较稳 W = rng.normal(0.0, 1.0, size=(M, m)) # 每个成员一组诱导点权重 (M, m) # 过程噪声：控制模型跟随数据的快慢 Q = np.full(m, 0.01) # 每个诱导点权重的随机游走标准差

M是最需要调的参数。太小（比如 20）时样本协方差抖动厉害，增益算出来会乱跳；太大（比如 500）单步就变慢，而且精度提升不明显。100 左右是个甜蜜点。

Q决定「模型有多不信任上一刻的自己」。Q 大，模型追得快但噪声也学进去；Q 小，平滑但跟不上突变。做设备退化这种缓变量，我给 1e-3 量级；做快速变化的信号，给 1e-2 甚至更大。这个参数没有万能的取值，只能拿验证集扫。

预测步就是加噪声，简单到不像话，但它承担了整个非平稳建模：

def enkf_predict(W, Q, rng): """随机游走预测步：每个成员独立加过程噪声""" noise = rng.normal(0.0, 1.0, size=W.shape) * Q return W + noise

注意rng.normal(..., size=W.shape) * Q这里是逐元素乘，Q是长度为 m 的向量，会广播到每一行——这允许你对不同诱导点给不同的漂移强度，比如边缘区域给大一点。

3.3 用样本协方差算卡尔曼增益

这是整个项目的心脏。函数接收新观测(x_new, y_new)，返回更新后的集成。

def enkf_update(W, Z, x_new, y_new, sigma_n, ls, sf, infl=1.05): M = W.shape[0] phi = rbf(x_new[None, :], Z, ls, sf) # (1, m) 观测点的特征 y_hat = W @ phi.T # (M, 1) 每个成员的预测值 # 协方差膨胀：把距平放大一点，防止集成过度收敛 A = (W - W.mean(0)) * infl # (M, m) Y = (y_hat - y_hat.mean(0)) * infl # (M, 1) Cwy = A.T @ Y / (M - 1) # (m, 1) 权重-观测协方差 Cyy = Y.T @ Y / (M - 1) + sigma_n ** 2 # (1, 1) 观测方差 + 噪声 K = Cwy @ np.linalg.inv(Cyy) # (m, 1) 卡尔曼增益 # 扰动观测：给每个成员配一个带噪的 y，防止集成塌缩 y_pert = y_new + rng.normal(0.0, sigma_n, size=(M, 1)) W_new = W + (K @ (y_pert - y_hat).T).T return W_new

逐段说明：

phi是新观测点在诱导点上的特征，长度 m，和权重的维度对齐。y_hat = W @ phi.T这一步把所有成员各自的预测都算出来了，形状 (M, 1)，这是后面所有统计量的基础。

A和Y是距平（去掉样本均值），infl是协方差膨胀因子，取值 1.0～1.1 之间。为什么要膨胀？因为有限集成估出来的协方差天然偏小，滤波器会过度自信、拒绝吸收新观测，加一点膨胀能显著提升跟踪能力。这个系数是我踩坑踩出来的，后面第 5 章细说。

Cwy是权重和观测的样本协方差（m×1），Cyy是观测的样本方差（标量）。np.linalg.inv(Cyy)在 1×1 的情况下有点杀鸡用牛刀，直接1.0 / Cyy更快，写成 inv 是为了和论文公式对照，也方便以后扩展成多维观测。

y_pert这一步最容易被忽略。如果所有成员都用同一个干净的y_new，更新后它们的距平会收缩到接近零，M 个成员变成一个，之后再也没法估协方差——这就是集成塌缩（ensemble collapse）。加上与观测噪声同量级的扰动，距平才能维持住。

最后的更新W + K @ (y_pert - y_hat).T是标准的 EnKF 增量式，形状对齐后再转置回来。sigma_n是观测噪声标准差，它直接决定了滤波器多信任新数据：sigma_n大，增益小，更新保守；sigma_n小，增益接近 1，几乎完全采纳观测。

3.4 超参数在线学习：状态增广加重采样

到上一节为止，ls、sf、sigma_n都是写死的。论文里一个重要贡献是把这三个超参数也学起来——做法是状态增广：把log(ls)、log(sf)、log(sigma_n)拼到权重后面，一起当成状态。

# 每个成员除了 m 个权重，再多带 3 个超参数（取对数保证正数） TH = np.tile(np.log([1.2, 1.0, 0.1]), (M, 1)) \ + rng.normal(0.0, 0.02, size=(M, 3)) # (M, 3)

EnKF 本身只处理「观测 = H·状态 + 噪声」这种形式，超参数没有直接观测，不能直接用增益公式更新。常见做法是两步走：先用每个成员自己的超参数算预测似然，按似然做加权重采样，让好成员留下、差成员被替换；再加一点微小扰动维持多样性。

def member_predict(W, TH, Z, x_new): """每个成员用自己的 ls/sf 算预测值，返回 (M,)""" out = np.empty(W.shape[0]) for i in range(W.shape[0]): ls, sf = np.exp(TH[i, 0]), np.exp(TH[i, 1]) phi = rbf(x_new[None, :], Z, ls, sf)[0] out[i] = W[i] @ phi return out def resample_theta(W, TH, y_hat, y_new, rng): """按预测似然对集成做系统重采样""" sig = np.exp(TH[:, 2]) loglik = -0.5 * ((y_new - y_hat) / sig) ** 2 - np.log(sig) # (M,) w = np.exp(loglik - loglik.max()) w /= w.sum() idx = rng.choice(len(w), size=len(w), p=w) # 多项式重采样 TH = TH[idx] + rng.normal(0.0, 0.005, size=TH.shape) # 抖动防退化 W = W[idx] return W, TH

member_predict里的 for 循环看起来低效，但 M 只有 100 左右，每次几十微秒，比向量化后的一堆索引操作还清楚。loglik是各成员的高斯对数似然，-log(sig)那项是归一化常数，少了它重采样会系统性偏向小噪声成员。

rng.choice用似然做概率重采样，好成员被复制、差成员被淘汰。后面那行+ normal(0, 0.005)千万别省——这是超参数版的「扰动观测」，不加的话几轮之后所有成员的超参数就完全一致了，学习过程提前结束。

注意：重采样会让权重和超参数同时被复制，等于给整个状态加了方差。如果你的 Casella 稳定性要求高（比如要做长期误差带），建议每 20 步重采样一次，而不是每步都做。

4. 跑起来：环境、目录结构与在线流式回归的闭环

聊完原理，该动手跑了。这一章从零搭一个能看结果的最小工程。

4.1 依赖与环境

这套东西对环境的胃口极小，NumPy 是唯一的硬依赖，画图用 Matplotlib，数据加载用 Pandas 就够。Python 3.8 以上都行，我实测 3.9 和 3.11 都没问题。

python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate pip install numpy matplotlib pandas

不需要 GPU，不需要 PyTorch/TensorFlow。整个算法就是几十行 NumPy 矩阵运算，这是它相比深度 GP 方案最大的工程优势——部署到边缘设备、工控机、甚至树莓派上都毫无压力。

有个小细节值得提醒：NumPy 的随机数生成器别再用np.random.seed那套老 API，统一用np.random.default_rng(seed)。新接口的生成器可以独立实例化，多线程或者多次实验时不会互相污染，我因为这个问题对一个不收敛的实验排查了两个小时。

4.2 目录怎么摆

原始论文代码结构往往按实验组织，工程化落地时我一般会拆成下面这种形状，职责分明，方便替换其中任意一块：

config.yaml里的参数全抽出来，是因为在线实验里你至少会扫十几次参数组合。硬编码在代码里的下场就是每改一次都要 git diff 一遍，很容易顺手把别的东西改坏。

4.3 主循环与评估指标

数据流式进来，每一步只做「预测 → 更新 → 评估」三件事：

import numpy as np W = rng.normal(0.0, 1.0, size=(M, m)) TH = np.tile(np.log([1.2, 1.0, 0.1]), (M, 1)) history = [] for t, (x_t, y_t) in enumerate(stream): W = enkf_predict(W, Q, rng) # 1) 预测：加过程噪声 W = enkf_update(W, Z, x_t, y_t, sigma_n, ls, sf) # 2) 更新：吸收观测 if t % 50 == 0: # 3) 每 50 步评估一次 y_hat = member_predict(W, TH, Z, x_t) # 各成员预测 (M,) mu, sd = y_hat.mean(), y_hat.std(ddof=1) # 集成均值与标准差 history.append({'t': t, 'mu': mu, 'sd': sd, 'rmse': np.sqrt((mu - y_t) ** 2)})

y_hat.mean()就是 GP 后验均值，y_hat.std(ddof=1)就是后验标准差——注意用ddof=1，样本协方差要用无偏估计，否则 M 小的时候会系统性低估不确定度，画出来的置信带偏窄，验证时容易误判成“模型挺自信的”。

评估别只看 RMSE。在线 GP 的价值在于它同时给出不确定度，所以两个指标一起看：

• RMSE：点预测准不准，反映均值质量
• NLPD（负对数预测密度）：0.5 * log(2πσ²) + (y − μ)² / (2σ²)，同时惩罚均值偏差和方差误标定

NLPD 是我最看重的指标。只看 RMSE 的话，一个把不确定度写成恒定值 1.0 的模型也能拿好看分数，但它在数据稀疏区的置信带完全不可信，做异常检测会直接翻车。

第一次跑建议先用合成数据：y = sin(x) + 0.1 * noise，x 从 −5 到 5 依次喂进去。前 20 个点会看到置信带很宽（先验占主导），到 100 个点之后带子收窄并贴合曲线，如果带子一直不收窄，十有八九是sigma_n设得太大或者膨胀因子过大。

5. 避坑与排查：集成塌缩、噪声标定和诱导点退化

这套代码跑通不难，跑准很难。下面五条是我在几轮实验里真踩过的，按「现象 → 原因 → 解决」记录。

5.1 集成侧的三个高频故障

现象一：跑几百步之后，所有成员的预测值几乎一模一样，置信带宽度趋近于零，新观测进来几乎不改变输出。

原因是集成塌缩。EnKF 用样本协方差代替理论协方差，如果更新时所有成员都被同一个干净观测拉向同一点，成员之间的距平就会指数衰减，几轮下来 M 个成员退化成 1 个，之后的增益全是噪声。解决就是第 3.3 节里那个y_pert，每个成员必须配一个独立扰动的观测，扰动量级取sigma_n，这是理论上的正确做法，不是权宜之计。同时把膨胀因子infl从 1.0 提到 1.05～1.1，给距平一点额外补偿。

现象二：一开始收敛挺好，跑到中途误差突然飙升，然后稳定在一个明显偏高的水平。

多半是过程噪声 Q 太小。Q 小意味着模型极度信任历史，一旦真实函数发生漂移（设备老化、工况切换），滤波器跟不上，只能靠微小的增益慢慢爬，表现出来的就是「卡住」。诊断方法很简单：打印每步的增益范数np.linalg.norm(K)，如果它持续小于 1e-3，就是 Q 给得太保守。往上调一个数量级通常就活了。

现象三：预测方差比实际残差小一大截，也就是过度自信，验证集里真实值频繁落在 95% 置信带外面。

两个可能。一是ddof用了 0，样本方差系统性偏小，改成ddof=1。二是膨胀因子加过头了，反而把协方差推得过大，导致滤波器过度反应、残差分布变厚。这两种情况现象相似但方向相反，排查方法是把当前步的sd和(y − mu)都记下来，算一下标准化残差z = (y − mu) / sd，它的标准差应该接近 1。明显小于 1 说明方差偏大，明显大于 1 说明偏小，这个诊断量比看 RMSE 直观得多。

5.2 数值与超参数上的坑

现象四：跑几十步后出现 NaN，或者增益矩阵里冒出 1e15 量级的数。

Cyy + sigma_n²接近零时求逆就炸了。什么时候会接近零？sigma_n设成 0，同时集成塌缩导致Y全零。两个因素一叠加，分母就是 0。修法是加个保底：Cyy = max(Cyy + sigma_n**2, 1e-8)，同时确保sigma_n不为 0（实在想表示无噪声，给 1e-6 而不是 0）。另外sf别设太大，RBF 核在sf=1e3这种量级下算指数容易溢出。

现象五：诱导点全挤在一小块区域，另一大片输入空间完全没有覆盖，那片区域的预测直接退化成先验均值。

这是诱导点退化。原因有两种：一是初始Z就用linspace铺得太窄，二是后面如果做了什么自适应调整（比如按梯度移动 Z），点会被吸引到数据密集区，稀疏区就荒了。诊断方式是拿Z和训练集 X 画在同一张图上，一眼就能看出来。解决就是把 Z 初始化成覆盖整个输入范围，或者用 k-means 在数据上聚一次，同时保留几个「哨兵点」固定在大范围两端不动。

提示：不要指望 EnKF 的随机游走去自动修正 Z 的位置。过程噪声是加在权重上的，诱导点位置是固定的，它不会自己跑。位置选错了，后面所有步骤都在补偿这个错误。

6. 二次开发：换核、批量观测与多步预测的改造点

工程里原样照搬的情况其实很少，多数时候你要改三四个地方。这一章把最常见的改造点列出来，都附上落点和验证方式。

换核函数是最容易的。所有核只要满足「输入 (n,d) 和 (m,d) 返回 (n,m)」这个签名，就能直接替换rbf。做设备退化这种带趋势的信号，我一般换 Matern 3/2：

换完之后别急着看 RMSE，先看 NLPD。核不匹配时 RMSE 可能还行，但方差会明显误标定，NLPD 立刻暴露问题。

批量观测是第二个高频需求。传感器常常一次回传一个窗口而不是单点，这时候把phi从 (1, m) 换成 (k, m)，y_hat变成 (M, k)，Cyy从标量变成 k×k 矩阵，np.linalg.inv就派上用场了。代码本身不用大改，因为第 3.3 节里已经是矩阵形式写的。唯一要注意的是观测扰动的维度要跟着改，size=(M, k)而不是(M, 1)，而且协方差要对角——不同传感器之间的噪声通常不相关，别用一个标量硬套。

多步预测需要在时间轴上迭代。给定 x_{t+1}，先做enkf_predict得到先验集成，但不要做 update（没有观测），直接算member_predict；再拿这个集成做下一步预测，如此递推。关键是递推时 Q 要累积，三阶预测的方差应该明显大于一步，如果两者差不多，说明enkf_predict里噪声加错了位置——我见过有人在循环外只加一次，看起来代码跑通了，置信带却一条直线，排查了半天才发现。

验证改造是否成功，我固定跑三个检查：合成数据上的覆盖率达到名义水平（95% 区间实际覆盖 93%～97%）、标准化残差标准差接近 1、以及单步更新耗时随 t 不增长（这是稀疏 GP 的核心承诺，如果耗时往上飘，多半是哪里的矩阵在悄悄变长）。

说个自己的教训。早期我做在线实验图省事，把sigma_n、infl、Q三个参数一起扫，48 组配置跑了一下午，结果最好那组在另一段数据上直接崩了。后来我改成一次只动一个参数，而且每个配置都在三个不同随机种子上跑，才看清infl和Q之间是强耦合的——它们都在控制「模型多愿意接受新观测」，同时调等于在两个旋钮上追同一个目标。从那以后我每次改这类在线滤波器，都强制走一遍「固定随机种子、单参数扫描、标准化残差诊断」这三步，再急也不跳。希望帮到你。

本文还有配套的精品资源，点击获取