› Kriging 代理模型 模拟器 返回
代理模型 / 机器学习

Kriging 代理模型 模拟器 — 高斯过程回归

用 Kriging（高斯过程回归）将高成本的 CFD 或 FEM 评估替换为少量训练点的代理模型。通过改变 SE / Matérn 核函数和相关长度 ℓ、块金方差 σ_n² 等，直观学习预测均值和95%置信区间、对数边际似然、RMSE 的行为。

参数设置
核函数
真实函数光滑性的先验假设
训练点数 N
高成本评估点数（CFD/FEM/实验）
相关长度 ℓ
核函数宽度。越大越假设光滑的函数
信号方差 σ_f²
函数的振幅尺度
块金方差 σ_n²
观测噪声的方差。确定性 CAE 设为 1e-6~1e-3
真实函数振幅 A
f(x) = A·sin(2πx) 的振幅
计算结果
—
训练点数 N
—
平均采样间隔
—
预测标准差 σ_pred
—
95% 置信区间宽度
—
RMSE 预测误差
—
对数边际似然
—
Kriging 插值 — 真函数、训练点、预测±95%CI

白线：真函数 A·sin(2πx)，红圆：训练点（含观测噪声），蓝线：Kriging 预测均值，蓝带：95% 置信区间。观测噪声动态波动显示。

Kriging 预测 + 置信区间（Chart.js）
对数边际似然 logML vs 相关长度 ℓ
理论和主要公式

$$\hat\mu(x_*) = k_*^T K^{-1} y,\qquad \hat\sigma^2(x_*) = k(x_*,x_*) - k_*^T K^{-1} k_*$$

K = 核矩阵（N×N，K_{ij}=k(x_i,x_j)+σ_n²δ_{ij}），k_* = 预测点与训练点的协方差向量。同时提供后验均值 μ̂ 和方差 σ̂²。

$$k_{SE}(r) = \sigma_f^2 \exp\!\left(-\frac{r^2}{2\ell^2}\right)$$

平方指数（RBF）核函数。r=|x−x'|，ℓ=相关长度，σ_f²=信号方差。假设无限次可微的「最光滑的」函数族。

$$k_{M5/2}(r) = \sigma_f^2\left(1+\frac{\sqrt{5}\,r}{\ell}+\frac{5r^2}{3\ell^2}\right)\exp\!\left(-\frac{\sqrt{5}\,r}{\ell}\right)$$

Matérn 5/2 核函数（贝叶斯优化的默认选择）。2阶可微，比 SE 能表现更急剧的变化。

$$\log p(y\mid X,\theta) = -\tfrac12 y^T K^{-1} y - \tfrac12 \log|K| - \tfrac{N}{2}\log(2\pi)$$

对数边际似然。第1项=数据拟合，第2项=模型复杂度惩罚。通过对 θ=(ℓ,σ_f²,σ_n²) 最大化来估计超参数。

Kriging 代理模型（高斯过程回归）

🙋
「代理模型」就是用来替代真实 CFD 的，对吧？但多项式回归不也可以吗？为什么要用 Kriging？
🎓
很好的问题。Kriging 与其他回归的本质区别是「不仅输出预测值，还同时输出预测的不确定性（方差）」。多项式或 RBF 插值只返回「某个值」，但 Kriging 在每个点 x_* 都返回「均值 μ̂(x_*)」和「标准差 σ̂(x_*)」。这样就能量化「这个采样点附近自信」还是「这个地方很稀疏所以不确定」。这正是贝叶斯优化中期望改进（Expected Improvement）的基础——它用这个不确定性来决定下一个评估点。
🙋
明白了。那左边的「相关长度 ℓ」和「块金方差 σ_n²」怎么确定呢？我随便改一下 RMSE 就变了很多，不知道什么是对的…
🎓
答案就是「最大化对数边际似然（logML）的 ℓ」。右下角的图表就是把 logML 画成 ℓ 的函数，找到峰值的位置就行。实际工作中用 BFGS 或 L-BFGS 数值优化，但观测点太少时容易陷入局部最优，所以标准做法是「多起点优化」——从 ℓ=0.05, 0.5, 5 这样的多个初始值开始。GPyOpt、scikit-learn 的 GaussianProcessRegressor、SMT、Dakota 这些库都在内部帮你自动做了这个优化。
🙋
SE 和 Matérn 核的预测曲线差别很大，SE 看起来最光滑…什么时候要用 Matérn 呢？
🎓
SE 是「无限次可微」的假设，对很光滑的现象特别强。比如机翼的升力系数随攻角变化，SE 就足够了。但座屈、接触、相变这种「某处行为急剧变化」的物理，SE 的无限光滑假设反而让它在不连续点周围振荡。这时候 Matérn 3/2（只1阶可微）或 Matérn 5/2（2阶可微）更「诚实」。贝叶斯优化圈子里 Matérn 5/2 几乎是标配，不确定时就从 5/2 开始，再用交叉验证和残差图和 SE 比较。
🙋
N=8 时 RMSE 有 1.07，判定是 NG。增加到 N=20 就会好吗？要多少个点才够？
🎓
没错，加到 N≈20，sin(2πx) 这样的函数 RMSE 就会掉到 0.2 以下。经验法则是「1维需要10~20点、2维30~50点、5维100~200点」。所需点数约按 10^d 这样随维数增长，这叫「维数灾难」。实际工作中维数超过10时，就该考虑从纯 Kriging 切到 PCE、深度核学习这些方向了。但 CAE 设计优化的典型规模（5~10维、100~500点）用 Kriging 正好。
🙋
最后一个问题。CFD 结果这种「确定性、没有噪声」的输出，也要加块金方差 σ_n² 吗？
🎓
从数值稳定性看，一定要加。理论上确定性输出 σ_n²=0 就够，但那样 K 矩阵几乎奇异，K^{-1} 会爆炸。所以必须加「jitter」——个小数字（1e-6~1e-3），这就是块金的作用，是正则化项。而实验数据的话，就直接用测量误差的方差。CAE 和实验混合的多保真 Kriging 中，每类数据都配一个独立的块金，这是标准做法。

常见问题

本质上是同一种方法。Kriging 是20世纪50年代由南非矿业工程师 Krige 提出的地质统计学术语，用于补间地下矿床品质，而机器学习社区中同样的方程被称为高斯过程回归（GPR）。两者都是「通过观测点的加权和生成预测值，同时输出预测方差」的线性贝叶斯估计，通过核函数 k(x,x') 定义光滑性。实现上的区别只是 Kriging 用变异函数 γ(h) 表示，而 GPR 用协方差函数 k(r) 表示。在 CAE 和优化领域通常称为 Kriging，在机器学习和不确定性量化领域通常称为 GPR。
SE（平方指数、RBF）核函数假设「最光滑的」无限次可微函数。当真实响应光滑时（如空气动力学分析或热问题），首先选择 SE。Matérn 3/2 表示只有一阶可微的函数（具有角或不连续导数），适用于结构失稳或接触问题等响应急剧变化的情况。Matérn 5/2 介于两者之间，是贝叶斯优化的事实标准。建议从 Matérn 5/2 开始，通过残差图和交叉验证与 SE 比较。
通常通过最大化对数边际似然（log marginal likelihood）来估计超参数。本工具的第二个图表显示了对数边际似然作为 ℓ 的函数，其峰值位置就是最大似然估计值。在实际应用中使用 BFGS 或 L-BFGS 进行数值优化，但当观测点较少时容易陷入局部最优，因此标准做法是从多个初始值（ℓ=0.05, 0.5, 5 等）开始的多起点优化。块金方差 σ_n² 代表观测噪声的方差：对于确定性 CAE 输出，使用较小的值（1e-6~1e-3）；对于实验数据，直接使用测量误差的方差。
标准 GPR 的计算复杂度为 O(N³)，内存复杂度为 O(N²)，所以当训练点数 N 超过约 2000 时，普通 PC 会遇到困难。通过 ARD（自动相关度确定）学习每个维度的相关长度，可以处理约 20 维的问题；超过此维度会出现「维数灾难」，性能下降。对于高维大规模数据，应切换到稀疏 GP（诱导点法）、SVGP、深度核学习、PCE（多项式混沌展开）或神经网络等方法。对于典型的 CAE 设计优化情形（5~10 维、100~500 个点），标准 Kriging 是最佳选择。

实际应用

贝叶斯优化（Bayesian Optimization）：在 CFD 或 FEM 的每次评估需要数小时的设计优化中，将 Kriging 代理模型与期望改进（Expected Improvement）获得函数结合，提出「下一步评估哪里能获得最大改进」。常见于飞机翼型形状优化、火箭喷嘴设计、化工厂运行条件优化等，通常用数十至数百次高成本评估达到全局最优。Google Vizier、Facebook Ax、scikit-optimize 等开源工具广泛使用。

CAE 结果的代理模型与稳健设计：汽车碰撞分析、电磁场分析、流体热耦合分析等每个案例需要数十分钟至数天的模拟，用 Kriging 建立代理模型，再用蒙特卡洛方法在厚度、材料常数、荷载等的不确定性下流运 1 万个案例的 UQ（不确定性量化）。直接运行 CAE 1 万次不现实，但代理模型能达到 0.1 秒/案例，一晚上就完成。Dakota、SMT、UQLab 等是标准工具。

地质统计学与资源勘探：Kriging 的发源地。从数十个钻井取得的矿石品质或地下水位数据出发，三维补间矿床品质分布图。广泛应用于石油天然气、稀土矿探测、地下水污染扩散预测等。通常与变异函数分析结合，确定空间相关结构，这是行业标准做法。

机器学习的不确定性量化：自动驾驶障碍物检测、医疗图像诊断等「预测可能错了要发警告」的场景，用 GPR 作辅助模型。深度核学习（DKL）用神经网络特征加 GP 后层的方式，即使对图像这样的高维输入也能给出预测方差。NASA 太空船导航、丰田驾驶辅助、基因泰克药物筛选等安全关键领域在推进应用。

常见误区与注意事项

最大的陷阱是「把 Kriging 的置信区间当成绝对真理」。GP 的事后方差 σ̂² 是「假设了正确的核和观测数据」前提下的方差，核与真实函数族不匹配就会系统偏离。比如把本来有不连续的阶梯函数用 SE（无限可微）核硬套，不连续处预测大错特错，但置信区间却说「很有把握」。实务必须通过交叉验证检查「实测值落在95%置信区间内的频率」是否接近公称的95%（calibration test）。

其次是「固定超参数一用到底」。在优化初期估的 ℓ 到了设计空间其他地方可能完全不适用。每加一个新评估点都要重新优化超参数才是正确做法。反过来，如果每轮都大幅波动，说明观测噪声设太小，或初始采样不好（应该用拉丁超立方采样 LHS 或 Sobol 序列）。

还有「Kriging 擅长外推」的错误认知。实际上完全相反——训练数据包围盒外的区域，预测均值快速回到先验均值（通常是0），方差发散到 σ_f²，也就是「我啥都不知道」。这本身是诚实的，但如果优化算法提议了边界外的点，就白浪费一次评估。贝叶斯优化要硬性约束设计空间，或在靠近边界处衰减期望改进。另一种方案是混合 PCE（多项式混沌展开）这类「全局多项式基」的方法，形成 PC-Kriging 混合体。

使用指南

  1. 在 1~50 范围内设置训练数据点数（numSamples），对应高成本分析（有限元法、流体分析等）的评估点数
  2. 选择核函数（平方指数/Matérn3-2/Matérn5-2），调整相关长度（lengthScale）在 0.1~5.0 范围内以控制函数光滑性
  3. 输入信号方差（signalVariance）和块金方差（nuggetVariance）来量化模型的不确定性，查看预测值、95%置信区间宽度、RMSE、对数边际似然

具体计算示例

飞机机翼弯曲刚度优化中，从 12 个有限元分析评估点（翼梁厚度范围 1.2~4.5 mm）构建 Kriging 模型。参数设为相关长度=0.8、信号方差=850、块金方差=2.5，在厚度 2.8 mm 处得到：预测弯曲刚度=245 N·m²、95%置信区间宽度=18.5 N·m²、RMSE=3.2 N·m²、对数边际似然=-8.6。由此不仅能预测未评估点的响应，还能同时获取不确定性，为下一个评估点的最优配置（拉丁超立方计划等）奠定基础。

实际工作中的注意事项

  1. 块金方差反映计算误差或制造偏差；在 0.5~5.0 范围内设置，太小容易过拟合，需谨慎
  2. Matérn5-2 核函数适合自动车身应力分析等光滑响应曲面；Matérn3-2 适合复合材料板屈曲荷载等中等光滑度应用
  3. 训练点较少（N<8）时置信区间变宽，全设计空间信度偏低，需增加采样点
  4. 对数边际似然低于 -15 时说明核参数过度调优，应检查输入数据的归一化（均值0、标准差1）