简介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 取在 30200 之间看输入维度和平滑程度太少了欠拟合太多了单步更新又开始变慢。提示诱导点的位置不是随便撒的。均匀铺在训练区间上是最省事的做法想要更好效果就用 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 Pm200 时就是 4 万个数还要做矩阵求逆Python 里跑几千步就明显吃力。EnKF 换了个思路不维护协方差维护一个样本集。用 M 个成员 {w⁽¹⁾, …, w⁽ᴹ⁾} 去近似后验分布均值就是样本均值协方差就是样本协方差P ≈ (1/(M-1)) Σᵢ (w⁽ⁱ⁾ − w̄)(w⁽ⁱ⁾ − w̄)ᵀ每步只需要 O(M·m) 的矩阵乘法而且天然并行m 再大也不怕。代价是采样误差——M 太小协方差估不准这就需要后面讲的「观测扰动」和「协方差膨胀」来补救。2.3 EnKF 与粒子滤波、递归最小二乘的差别同样是序贯更新这三者的取舍值得说清楚选错了会白写一堆代码。方法后验表示单步代价非线性观测典型问题递归最小二乘 (RLS)单一均值 协方差矩阵O(m²)不支持只给点估计无不确定度粒子滤波 (PF)带权粒子集O(M·m)支持高维退化严重需重采样EnKF等权集成O(M·m)近似支持需要扰动观测防塌缩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, ls1.0, sf1.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, ls1.2, sf1.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, sizeW.shape) * Q return W noise注意rng.normal(..., sizeW.shape) * Q这里是逐元素乘Q是长度为 m 的向量会广播到每一行——这允许你对不同诱导点给不同的漂移强度比如边缘区域给大一点。3.3 用样本协方差算卡尔曼增益这是整个项目的心脏。函数接收新观测(x_new, y_new)返回更新后的集成。def enkf_update(W, Z, x_new, y_new, sigma_n, ls, sf, infl1.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.01.1 之间。为什么要膨胀因为有限集成估出来的协方差天然偏小滤波器会过度自信、拒绝吸收新观测加一点膨胀能显著提升跟踪能力。这个系数是我踩坑踩出来的后面第 5 章细说。Cwy是权重和观测的样本协方差m×1Cyy是观测的样本方差标量。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), sizelen(w), pw) # 多项式重采样 TH TH[idx] rng.normal(0.0, 0.005, sizeTH.shape) # 抖动防退化 W W[idx] return W, THmember_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 目录怎么摆原始论文代码结构往往按实验组织工程化落地时我一般会拆成下面这种形状职责分明方便替换其中任意一块文件职责改动频率kernel.pyRBF / Matern 核超参数用字典传入低gp_enkf.pypredict / update / resample 三个核心方法中data_stream.py从 CSV / 传感器读数据产出 (x, y) 迭代器高train_online.py主循环日志与指标中config.yamlM、m、Q、sigma_n、膨胀因子极高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(ddof1) # 集成均值与标准差 history.append({t: t, mu: mu, sd: sd, rmse: np.sqrt((mu - y_t) ** 2)})y_hat.mean()就是 GP 后验均值y_hat.std(ddof1)就是后验标准差——注意用ddof1样本协方差要用无偏估计否则 M 小的时候会系统性低估不确定度画出来的置信带偏窄验证时容易误判成“模型挺自信的”。评估别只看 RMSE。在线 GP 的价值在于它同时给出不确定度所以两个指标一起看RMSE点预测准不准反映均值质量NLPD负对数预测密度0.5 * log(2πσ²) (y − μ)² / (2σ²)同时惩罚均值偏差和方差误标定NLPD 是我最看重的指标。只看 RMSE 的话一个把不确定度写成恒定值 1.0 的模型也能拿好看分数但它在数据稀疏区的置信带完全不可信做异常检测会直接翻车。第一次跑建议先用合成数据y sin(x) 0.1 * noisex 从 −5 到 5 依次喂进去。前 20 个点会看到置信带很宽先验占主导到 100 个点之后带子收窄并贴合曲线如果带子一直不收窄十有八九是sigma_n设得太大或者膨胀因子过大。5. 避坑与排查集成塌缩、噪声标定和诱导点退化这套代码跑通不难跑准很难。下面五条是我在几轮实验里真踩过的按「现象 → 原因 → 解决」记录。5.1 集成侧的三个高频故障现象一跑几百步之后所有成员的预测值几乎一模一样置信带宽度趋近于零新观测进来几乎不改变输出。原因是集成塌缩。EnKF 用样本协方差代替理论协方差如果更新时所有成员都被同一个干净观测拉向同一点成员之间的距平就会指数衰减几轮下来 M 个成员退化成 1 个之后的增益全是噪声。解决就是第 3.3 节里那个y_pert每个成员必须配一个独立扰动的观测扰动量级取sigma_n这是理论上的正确做法不是权宜之计。同时把膨胀因子infl从 1.0 提到 1.051.1给距平一点额外补偿。现象二一开始收敛挺好跑到中途误差突然飙升然后稳定在一个明显偏高的水平。多半是过程噪声 Q 太小。Q 小意味着模型极度信任历史一旦真实函数发生漂移设备老化、工况切换滤波器跟不上只能靠微小的增益慢慢爬表现出来的就是「卡住」。诊断方法很简单打印每步的增益范数np.linalg.norm(K)如果它持续小于 1e-3就是 Q 给得太保守。往上调一个数量级通常就活了。现象三预测方差比实际残差小一大截也就是过度自信验证集里真实值频繁落在 95% 置信带外面。两个可能。一是ddof用了 0样本方差系统性偏小改成ddof1。二是膨胀因子加过头了反而把协方差推得过大导致滤波器过度反应、残差分布变厚。这两种情况现象相似但方向相反排查方法是把当前步的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 核在sf1e3这种量级下算指数容易溢出。现象五诱导点全挤在一小块区域另一大片输入空间完全没有覆盖那片区域的预测直接退化成先验均值。这是诱导点退化。原因有两种一是初始Z就用linspace铺得太窄二是后面如果做了什么自适应调整比如按梯度移动 Z点会被吸引到数据密集区稀疏区就荒了。诊断方式是拿Z和训练集 X 画在同一张图上一眼就能看出来。解决就是把 Z 初始化成覆盖整个输入范围或者用 k-means 在数据上聚一次同时保留几个「哨兵点」固定在大范围两端不动。提示不要指望 EnKF 的随机游走去自动修正 Z 的位置。过程噪声是加在权重上的诱导点位置是固定的它不会自己跑。位置选错了后面所有步骤都在补偿这个错误。6. 二次开发换核、批量观测与多步预测的改造点工程里原样照搬的情况其实很少多数时候你要改三四个地方。这一章把最常见的改造点列出来都附上落点和验证方式。换核函数是最容易的。所有核只要满足「输入 (n,d) 和 (m,d) 返回 (n,m)」这个签名就能直接替换rbf。做设备退化这种带趋势的信号我一般换 Matern 3/2核适用场景注意点RBF平滑函数、温度曲线ls过小会退化成插值Matern 3/2有拐点、轨迹跟踪实现要处理√3 r / ls的数值稳定性周期核日周期、工况循环周期超参数要单独在线学换完之后别急着看 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_{t1}先做enkf_predict得到先验集成但不要做 update没有观测直接算member_predict再拿这个集成做下一步预测如此递推。关键是递推时 Q 要累积三阶预测的方差应该明显大于一步如果两者差不多说明enkf_predict里噪声加错了位置——我见过有人在循环外只加一次看起来代码跑通了置信带却一条直线排查了半天才发现。验证改造是否成功我固定跑三个检查合成数据上的覆盖率达到名义水平95% 区间实际覆盖 93%97%、标准化残差标准差接近 1、以及单步更新耗时随 t 不增长这是稀疏 GP 的核心承诺如果耗时往上飘多半是哪里的矩阵在悄悄变长。说个自己的教训。早期我做在线实验图省事把sigma_n、infl、Q三个参数一起扫48 组配置跑了一下午结果最好那组在另一段数据上直接崩了。后来我改成一次只动一个参数而且每个配置都在三个不同随机种子上跑才看清infl和Q之间是强耦合的——它们都在控制「模型多愿意接受新观测」同时调等于在两个旋钮上追同一个目标。从那以后我每次改这类在线滤波器都强制走一遍「固定随机种子、单参数扫描、标准化残差诊断」这三步再急也不跳。希望帮到你。本文还有配套的精品资源点击获取
)
