简介这是一套基于Kriging克里金方法的代理模型MATLAB实现代码包面向需要开展空间插值预测、代理模型构建与自适应建模的科研人员和工程师。资源以KrigingModelCode为核心完整覆盖协方差函数选择、模型参数估计、Kriging线性系统求解、预测误差分析以及自适应迭代更新流程帮助用户实现从静态空间插值到动态优化代理模型的完整链路。压缩包共18个文件包含16个.m程序脚本、1份dace.pdf文档和1个changelog更新说明整体约838KB其中.m脚本覆盖高斯、指数、球面等常用协方差函数以及采样、回归、拟合、预测等关键流程模块PDF提供理论推导与DACE工具箱使用指南便于边读边练。该资源已有1640人学习下载适合具备一定MATLAB基础、希望在项目中快速集成Kriging代理模型或深入理解空间统计插值原理的读者。通过阅读代码与文档可掌握自适应Kriging的建模思路并可直接复用采样、回归、参数估计与预测模块降低二次开发成本。1. 别把昂贵仿真当黑盒Kriging代理模型与DACE源码包的定位一次CFD仿真跑一个算例要三四个小时而优化算法动不动就要评估上千个设计点。传统响应面在强非线性问题上精度崩塌神经网络又需要海量样本才能收敛。Kriging代理模型也叫高斯过程回归、DACE用二三十个样本就能同时给出预测值和预测方差是工程优化里性价比最高的插值方案。这份kriging.zip打包了DACE工具箱的完整MATLAB源码dacefit.m负责模型拟合predictor.m负责预测lhsamp.m和gridsamp.m负责实验设计regpoly0/1/2.m和corrgauss.m、correxp.m、corrspherical.m等相关系数函数全部开放。适合做结构优化、可靠性分析、参数标定的工程师也适合想把自适应采样逻辑接进自己优化框架的算法研究者。2. Kriging插值的数学原理与DACE核心设计相关函数与回归多项式的取舍2.1 从回归残差到随机过程Kriging模型的基本形式Kriging把目标函数拆成两部分一个显式的回归趋势项和一个隐含的随机偏差项写成y(x) f(x)^T β z(x)。回归项f(x)^T β描述全局趋势z(x)是均值为零的平稳高斯过程用来捕捉局部偏差。任何两个点x_i、x_j的偏差之间都有一个由相关函数定义的协方差Cov(z(x_i), z(x_j)) σ² R(x_i, x_j; θ)这个R只与两点之间的距离和超参数θ有关。DACE的拟合过程就是利用观测数据估计回归系数β、过程方差σ²和超参数θ。θ的作用非常直接它决定了每个输入维度上距离多远才算不相关。某个维度上的θ越大模型在这个方向上的相关长度越短拟合出来的曲面局部起伏越剧烈。dacefit.m内部通过最大似然估计逐维搜索θ这也是整个工具包里计算开销最大的部分。2.2 相关函数家族corrgauss、correxp、corrlin、corrcubic、corrspherical怎么选相关函数是Kriging模型的灵魂它决定插值曲面的光滑程度。DACE在corr*.m这组文件里实现了7种相关函数最常用的是高斯型、指数型和样条型。函数文件数学形式光滑度典型适用场景corrgauss.mexp(-Σθ_k·d_k²)无穷阶光滑大部分连续仿真问题首选correxp.mexp(-Σθ_k·d_k)连续但不可导噪声大、响应剧烈变化的问题correxpg.mexp(-Σθ_k·d_k^p)0p≤2由p控制需要额外控制光滑度时corrlin.mmax(0, 1-Σθ_k·d_k)分段线性快速试算、低精度场景corrcubic.m三次多项式核一阶光滑工程中兼顾精度与稳定corrspherical.m球状核超过变程归零分段光滑地质统计背景的插值corrspline.mB样条核分段光滑样本点较多、局部波动明显实际项目里我一般固定用corrgauss。它的无穷阶光滑性质对大多数CAE仿真响应都成立而且correxpg在p2时就是高斯核没必要多引入一个待估参数。只有看到预测曲面出现不合理的过度振荡时才换corrcubic或corrspherical。2.3 回归多项式regpoly0/1/2扮演的角色regpoly0.m是常数回归regpoly1.m是线性回归regpoly2.m是带交叉项和平方项的二次回归。很多人误以为回归阶数越高越好实际上在样本量少的时候高阶回归会和随机偏差项争夺解释权。DACE官方文档也建议初始建模用regpoly0或regpoly1只有确认响应存在明显趋势或者建模点数超过100时再考虑regpoly2。回归项与相关函数是配合关系。趋势项已经把平滑部分吃掉了偏差项只需要处理剩余残差。如果把回归设成regpoly2相关函数的θ往往会倾向于很小的值因为趋势项已经解释掉大部分变化偏差项只剩白噪声。2.4 从文件列表看DACE接口dacefit与predictor的参数契约整个工具箱的调用关系非常精简核心只有两个函数。dacefit.m做训练predictor.m做预测。% 训练模型 [dmodel, perf] dacefit(S, Y, regr, corr, theta0, lob, upb); % 预测新点 [y, or] predictor(x, dmodel);dacefit的前四个参数分别是样本点矩阵Sn行d列、响应值Yn行1列、回归函数句柄和相关函数句柄。后三个参数都是关于θ的theta0是搜索初值lob和upb分别是下界和上界向量维度与输入维度一致。参数含义设置建议Sn×d 样本矩阵推荐用lhsamp先生成各行不能重复Yn×1 响应向量最好做归一化数值范围过大容易让似然函数溢出regr回归函数句柄默认regpoly0趋势明显用regpoly1corr相关函数句柄默认corrgauss换其他核需重新调θ边界theta0θ初值向量取lob和upb的几何中心比如0.1*ones(1,d)lob/upbθ上下界下界不能取0上界取5到10足够dacefit返回的dmodel是一个结构体里面保存了回归系数beta、过程方差sigma2、相关参数theta以及原始样本数据。predictor做预测时会把x和dmodel里的样本重新算一遍相关矩阵所以样本量过大时预测速度也会明显变慢。3. MATLAB实战用lhsamp采样、dacefit训练、predictor预测一个二维Kriging模型3.1 实验设计lhsamp生成拉丁超立方样本代理模型的样本点质量直接决定拟合上限。均匀设计保证空间覆盖但每个维度上均匀分布产生的组合数量会随维度爆炸。拉丁超立方采样LHS把每个维度分成等概率区间在每个区间里恰好取一个点既保证覆盖又控制样本量。DACE的lhsamp.m实现的就是这个方法。% 定义测试函数 f (x) sin(6*x(:,1)) cos(4*x(:,2)) 0.5*x(:,1).*x(:,2); % 生成30个二维样本每个维度取值区间[0,1] S lhsamp(30, 2); Y f(S); % 查看样本分布是否覆盖均匀 scatter(S(:,1), S(:,2), 36, Y, filled); colorbar;lhsamp(n, k)的第一个参数是样本数第二个是维度数输出是n×k矩阵每一列独立地在[0,1]区间内做分层抽样。层数等于样本数n所以n太小比如小于10时分层效果不明显。工程上二维问题一般取20到50个初始样本再配合后面讲的自适应加点。3.2 用dacefit拟合超参数θ拿到样本和响应后需要给dacefit指定θ的搜索范围。θ决定了相关函数在各自维度上的衰减速率给太宽会让搜索陷入平坦区域给太窄又会错过最优值。% 设置theta初值和边界 theta0 0.2 * ones(1, 2); % 初值取几何中心附近 lob 1e-3 * ones(1, 2); % 下界略大于0 upb 5 * ones(1, 2); % 上界经验值 % 训练Kriging模型 [dmodel, perf] dacefit(S, Y, regpoly1, corrgauss, theta0, lob, upb); % 查看拟合得到的超参数 disp(dmodel.theta);这段代码里theta0给的是0.2这意味着初始认为两个维度在距离0.2左右相关性开始明显衰减。lob给1e-3而不是0是因为相关函数在θ趋近0时几乎变成常数函数似然函数会退化导致搜索不稳定。upb给5对应相关长度约0.45已经足够覆盖[0,1]区间的一半。dmodel.theta输出的就是最大似然估计的结果如果两个维度数值差异大说明响应在这两个方向上的变化剧烈程度不同。3.3 用predictor做预测并计算MSEKriging相比其他插值方法的核心优势就是能输出预测方差。predictor的第二个返回值or里带了一个mse字段表示预测值的均方误差。这个MSE不是误差估计而是模型对自身不确定性的度量样本点附近MSE小远离样本点MSE迅速变大。% 用gridsamp生成网格化预测点 bounds [0 1; 0 1]; Xg gridsamp(bounds, 40); % 预测响应值和均方误差 [Yg, or] predictor(Xg, dmodel); MSEg or.mse; % 可视化预测曲面和不确定性 subplot(1,2,1); surf(reshape(Xg(:,1), 41, 41), reshape(Xg(:,2), 41, 41), reshape(Yg, 41, 41)); title(预测曲面); subplot(1,2,2); surf(reshape(Xg(:,1), 41, 41), reshape(Xg(:,2), 41, 41), reshape(MSEg, 41, 41)); title(预测方差);gridsamp(bounds, q)按每个维度划分成q段返回(q1)^d个网格点第一维变化最快。这里的reshape(..., 41, 41)就是把展平的点阵还原成网格方便surf绘图。预测方差图上样本点所在位置会形成明显的谷地远离样本的位置方差抬升这个信息就是后面自适应加点策略的核心依据。3.4 模型验证用dsmerge和真实误差评估模型建完不能直接用先验证。最直接的方式是用留一交叉验证每次拿掉一个样本用剩余样本训练预测被拿掉的点统计误差。DACE的dsmerge.m虽然主要是数据合并工具但配合交叉验证可以让数据管理更干净。% 计算真实误差在网格点上对比预测值和真实函数 Ytrue f(Xg); rmse sqrt(mean((Yg - Ytrue).^2)); maxerr max(abs(Yg - Ytrue)); fprintf(RMSE: %.4f, Max Error: %.4f\n, rmse, maxerr); % 手动做一个简单的留一验证 n size(S, 1); cv_err zeros(n, 1); for i 1:n St S; Yt Y; St(i,:) []; Yt(i) []; dm_tmp dacefit(St, Yt, regpoly1, corrgauss, theta0, lob, upb); [yp, ~] predictor(S(i,:), dm_tmp); cv_err(i) abs(yp - Y(i)); end fprintf(LOO平均误差: %.4f\n, mean(cv_err));RMSE反映整体精度留一验证能看出模型有没有过度依赖某几个样本点。如果某个点的留一误差远大于其他点说明这个点附近的响应变化剧烈当前样本密度不够后续加点应该优先补这一带。交叉验证的代价是每轮都要重新训练样本100个以内时可以接受超过200就比较慢了。4. 自适应Kriging加点策略U函数与EI期望改进的完整迭代流程4.1 初始模型为什么不够从预测方差看改进空间固定样本量的Kriging模型像一张绷在样本点上的膜样本点处完全贴合离开样本点就松弛下来。膜松弛的程度就是前面看到的MSE。自适应Kriging的思路是不断寻找预测方差大、或者对目标函数影响最敏感的位置把这些位置补进样本集重新训练让膜逐渐绷紧。代理模型这几年在AI辅助工程优化里重新热起来但很多人直接把神经网络当代理样本一少就过拟合。Kriging这类统计代理模型反而更稳因为它自带不确定性估计优化算法可以明确知道哪里还没学好。这个性质让基于Kriging的EGO、AK-MCS等方法成为可靠性分析和昂贵仿真的标准做法。4.2 U函数与EI准则两种最经典的加点策略加点策略解决下一批样本点取在哪里的问题。U函数用于可靠性分析当极限状态函数g(x)0定义失效面时定义U(x) |ŷ(x)| / σ(x)σ(x)是预测标准差。U值越小这个点越有可能被误分类到失效面错误的一侧。U2对应的误分类概率约为2.28%U3时降到0.13%。期望改进准则EI用于全局优化EI(x) (ymin - ŷ(x))·Φ((ymin - ŷ(x))/σ(x)) σ(x)·φ((ymin - ŷ(x))/σ(x))其中ymin是当前样本中的最优响应Φ是标准正态CDFφ是PDF。EI同时考虑预测值相对当前最优的改进量和预测方差在开发和探索之间自动平衡。4.3 AK-MCS完整迭代流程从候选池里选点加入把这两种策略落到DACE上完整流程是先生成大规模蒙特卡洛候选池反复用当前模型预测候选池选出U最小或EI最大的点仿真后用dsmerge合并数据并重新训练直到满足停止条件。% 初始实验设计 S0 lhsamp(12, 2); Y0 f(S0); [dmodel, ~] dacefit(S0, Y0, regpoly1, corrgauss, 0.2*ones(1,2), 1e-3*ones(1,2), 5*ones(1,2)); % 生成候选池 Sc lhsamp(100000, 2); % AK-MCS 迭代 S S0; Y Y0; for iter 1:100 % 预测候选池 [Yc, or] predictor(Sc, dmodel); sigma_c sqrt(or.mse); % 计算U值找最小U U abs(Yc) ./ sigma_c; [Umin, idx] min(U); % 停止条件最小U大于2认为误分类风险足够低 if Umin 2 fprintf(迭代收敛共加点%d次\n, iter - 1); break; end % 把最不确定的点加入训练集 S [S; Sc(idx, :)]; Y [Y; f(Sc(idx, :))]; % 重新训练 [dmodel, ~] dacefit(S, Y, regpoly1, corrgauss, dmodel.theta, 1e-3*ones(1,2), 5*ones(1,2)); end这段代码里有一个关键细节重新训练时theta0不再用固定值而是用上一轮dmodel.theta。这样做是因为样本微调后最优θ不会突变用上一轮的θ做初值能显著加快收敛。预测候选池那一步是计算瓶颈10万个候选点乘几百个训练样本每轮都要算一次相关矩阵训练集超过200个点后循环会明显变慢。4.4 停止准则与计算量大时的替代方案准则公式适用场景说明U函数停止min(U(x)) ≥ 2可靠性分析误分类概率约2.28%工程上常用U函数严格停止min(U(x)) ≥ 3高可靠度需求误分类概率降至0.13%EI最大改进max(EI(x)) ε全局优化改进量低于阈值就停预算上限累计仿真次数Nmax所有场景兜底准则必须加计算量大时还有一个替代办法不要每加一个点就重新训练。DACE的dacefit是全局重新拟合和真正的增量更新不一样。如果候选池有10万点可以每轮选2到3个U值最小的点一起加入减少重训练次数代价是每轮选中的点可能集中在同一区域降低样本多样性。5. 工程落地的坑与技巧相关矩阵奇异、参数边界与数据归一化5.1 相关矩阵奇异或接近奇异dacefit最常报的错来自相关矩阵奇异。原因通常是样本点重复或过于接近两个几乎重合的点会让相关矩阵两行向量近乎线性相关矩阵求逆直接爆掉。这个坑在自适应加点后期特别容易出现U值最小的点往往落在已有样本附近加进去后新老样本距离小于1e-6矩阵立刻奇异。排查方法分三步先用dsmerge对样本去重再检查最小成对距离最后把lob下限调高。我一般把lob设在1e-3到1e-2之间太小会让相关矩阵在距离0附近过于平坦数值上更容易出问题。dsmerge.m的默认容差是1e-9如果两次仿真的输入完全一致但响应有微小差异需要手动打乱容差让其中一个点被过滤掉。5.2 theta初值与边界的设置技巧theta0、lob、upb这三个参数是最容易被忽略的。lob和upb定义了θ的搜索盒dacefit在盒子里逐维做fminbnd搜索。上界给太大比如100搜索空间里大部分区域的似然函数值几乎为零优化器很容易在平坦区停滞上界给太小真实相关长度超出搜索范围拟合出来的曲面明显欠光滑。实用的设置方式是先跑一遍默认参数theta00.1, lob1e-3, upb10看dmodel.theta落在哪个量级再把搜索范围收窄到该值周围半个数量级。比如第一次拟合得到theta[2.5, 0.4]第二次就把lob设为[0.5, 0.05]、upb设为[5, 0.5]。重拟合后对比perf.cv和perf里的交叉验证误差如果下降不明显说明第一次的结果已经接近全局最优。5.3 数据归一化与增量数据管理DACE的lhsamp生成的样本天然落在[0,1]区间但换成网格采样或外部数据源后各个维度可能数量级差异很大。比如一个是温度300到800一个是压力1e5到1e6如果不归一化相关函数的θ会在不同维度上互相牵制两个维度之间的θ数值失去可比性模型对量纲更大的维度自动赋予更高权重这是不对的。% 对样本做min-max归一化 S_min min(S0); S_max max(S0); S_norm (S0 - S_min) ./ (S_max - S_min); % 训练 [dmodel, ~] dacefit(S_norm, Y, regpoly1, corrgauss, 0.2*ones(1,d), 1e-3*ones(1,d), 5*ones(1,d)); % 预测时新样本也要用同一组min/max变换 x_new_norm (x_new - S_min) ./ (S_max - S_min); [y, ~] predictor(x_new_norm, dmodel);归一化时S_min和S_max是训练集上算出来的预测阶段必须沿用同一组值不能对预测点重新算否则分布不一致预测结果直接偏掉。dsmerge.m在自适应迭代里也要配合使用它能在合入新样本时顺带做去重和排序避免样本集越来越大后重复点堆积。最后留一个排查路径先拿工具包里的demo数据跑通dacefit和predictor把theta0从0.1换到1观察perf.cv变化再去掉一个样本看MSE是否在删掉的位置合理放大最后才接自己的采样函数和真实仿真。这样每一步出错都能定位到是采样、训练还是预测环节的问题。本文还有配套的精品资源点击获取