简介这是一份面向数据分析、仿真优化与不确定性量化研究者的Kriging模型与拉丁超立方抽样LHS学习资料包。资源以MATLAB源码为主体包含拉丁超立方抽样、Kriging拟合、相关函数、回归多项式、预测及网格生成等17个m文件另有1个示例数据mat文件与1篇DACE工具箱说明PDF共19个文件压缩包仅1.48MB适合熟悉MATLAB基础、希望掌握代理模型与空间插值的工程技术人员。目前已有434人学习下载。通过阅读代码与说明文档可系统理解LHS如何高效探索高维设计空间掌握Kriging建模流程、相关函数选择与预测精度评估并将其应用于实验设计、参数估计和复杂系统响应面构建等实际场景。1. LHS-Kriging 把仿真费用锐减的最直接思路在有限元仿真、CFD 和优化算法里代价最大的往往不是模型实现本身而是获得足够多的高可信样本。地质统计学里走出来的 Kriging 模型恰好扮演了“用插值替代一次重仿真”的角色而拉丁超立方抽样LHS则是在样本数量受控时保住空间覆盖度的一种采样策略。把这两者放到同一个压缩包 LHS-Kriging.zip 里意味着只需要 lhsu.m 生成一组均匀分布的训练点再用 dacefit.m 把 Kriging 模型的回归项、相关函数和超参数一起拟合出来后续点位无需重新跑昂贵仿真。一个常被忽略的事实是普通随机抽样在维数升高后边际覆盖并不均匀而 LHS 通过每维等概率分层让几十个样本点也能稳定支撑 Kriging 的参数估计。这套组合适合做代理模型、不确定性量化以及多目标优化前的响应面构建。2. DACE 工具箱里 Kriging 模型的拟合链路regpoly、corr 与 dacefitKriging 的工程表述是一个回归项加一个随机过程项y(x)∑f_j(x)β_jz(x)。z(x) 的协方差由相关函数控制DACE 工具箱把这条链路拆成四个可替换的部分回归函数 regpoly0/1/2、相关函数 corr*.m、全局优化拟合 dacefit.m以及后续插值用的 predictor.m。文件清单里的 corrcubic.m、corrgauss.m、corrspherical.m、corrlin.m、corrspline.m、correxpg.m、correxp.m 都是可插拔的相关函数。初学的人容易把这堆文件当成不同模型实际上它们只是同一个 Kriging 骨架上的零件。2.1 相关函数与回归基函数如何决定 Kriging 的“手感”回归基函数的数学形式很简单regpoly0 是常数regpoly1 是线性regpoly2 是二次多项式。回归项用于描述全局趋势相关函数用于描述局部偏差。实际工程里响应面变化剧烈时我会优先试 corrgauss因为它无限平滑能抓住连续迅速变化的趋势但高斯相关函数对 theta 初始值敏感局部优化容易陷入一个过大的尺度参数。correxp 是指数相关会让预测面更偏向分段变化适合不连续响应corrspherical 来自地统计学适合在有限距离外相关性降为零的场景。在 dacefit 内部参数估计通常通过极大似然或交叉验证完成也就是选择一组 theta 让模型对观测数据的重现概率最大。负对数似然函数里包含两项第一项来自拟合残差第二项来自相关矩阵的行列式。行列式项对 theta 的变化非常敏感这就是 lob/upb 设得过宽会导致优化器在非凸曲面上来回震荡的原因。相关函数表达式特征典型场景corrgaussexp(-theta*d^2)光滑响应面、CFD 代理模型correxpexp(-theta*d)带尖角、不连续响应correxpg广义指数指数族扩展可手动控制指数corrlin线性衰减简单插值、低维验证corrspherical球状变差函数地质与空间数据corrcubic三次多项式衰减中等平滑、鲁棒性较高corrspline样条类相关样本分布不规则我在实际项目中选相关函数不只看数学性质还会看 dacefit 收敛日志里的 theta 是否落在搜索边界上。如果某个 theta 一直贴着 upb说明该维度响应较弱或样本点不足以识别相关性此时把 regpoly 降为常数模型或换 corrspline 一般能改善条件数。需要注意的是DACE 内部的相关矩阵求逆对重复样本点非常敏感同一坐标出现两次会让 R 矩阵奇异dacefit 直接报错或返回毫无意义的 theta。2.2 dacefit 的典型调用与 theta 参数边界dacefit 的标准入口是把样本点 S、响应 Y、回归函数句柄、相关函数句柄以及 theta 的初始值和上下界一次性传入。一个最小调用长这样% 训练数据来自某个 2 维黑箱仿真 S [0.1 0.2; 0.3 0.4; 0.5 0.6]; Y [12.3; 15.6; 17.2]; theta0 [0.5 0.5]; % 每维一个初始值 lob [0.01 0.01]; % theta 下界 upb [20 20]; % theta 上界 [dmodel, perf] dacefit(S, Y, regpoly0, corrgauss, theta0, lob, upb);这段代码把 theta 初始值设为 0.5上下界按 0.01 到 20 给定。在 DACE 内部theta 是相关函数里的尺度参数theta 越大相关长度越小模型越倾向于只在样本点邻近生效。lob 不能设成 0否则相关矩阵求逆时数值不稳定upb 设得过大则会放大局部最优点与全局最优点的梯度差导致优化器在毫无特征的平坦区域上浪费迭代。perf 里一般保存了优化迭代信息或误差指标用来确认拟合过程是否收敛。如果仿真输出本身带噪声拟合出的模型会把噪声也“插值”进去此时不要迷信训练集的插值误差而是要看留出集上的预测误差。2.3 dsmerge.m 与多组样本合并的常见陷阱压缩包里还包含 dsmerge.m用于把两组样本和响应合并成同一数据集。我在做增量采样时经常用到它但踩过一个坑dsmerge 并不会重新对样本去重只是把 S、Y 等按行拼接。如果两组 LHS 样本恰好生成同一行坐标合并后仍然会出现重复点直接进入 dacefit 会触发相关矩阵奇异。因此合并前要先用 uniquetol 或去重逻辑对坐标做处理。% 用 lhsu 生成两组样本合并前做精确去重 X1 lhsu([0 0], [1 1], 10); X2 lhsu([0 0], [1 1], 10); S_all [X1; X2]; [S_dedup, ia] uniquetol(S_all, ByRows, true); Y_all Y_all(ia); % 同步保留响应这里 uniquetol 的 ByRows 选项把精度控制权交给默认容差如果输入变量的量级差异极大需要先对每列做归一化再判断重复否则容差会被大数值维吞掉。合并后的样本点数量越多dacefit 的矩阵分解越慢但相对地Kriging 的预测不确定性也会下降。3. LHS 初始化实战lhsu.m 与 latin_hs.m 的高维采样拉丁超立方抽样的核心思想是把每个输入维的区间等分成 n 个互不重叠的子区间在每一维中随机抽取一个代表值然后随机配对成 n 个样本点。由于每一维上每个子区间恰好出现一次样本在各维上的边际投影是均匀的。对比用 rand 生成的完全随机样本LHS 的方差更小尤其在 n 不大、维度高于 3 时优势明显。很多资料把 LHS 说成“随机拉丁超立方”这容易让人误以为它和蒙特卡洛只有工程实现上的差别实际上两者在样本结构上的差异是本质性的。3.1 LHS vs 蒙特卡洛随机抽样在低样本量下的差距随机抽样靠大数定律逼近分布但工程中能负担的样本往往只有 30-100 个。此时随机抽样容易出现跑偏某些区域样本密集另一些区域完全空白。Kriging 模型对这种空白区域几乎没有约束力因为它本质上是基于距离加权的插值器。LHS 则不会漏掉任何一段等概率区间所以在这个场景下经常被选作默认采样器。在机器学习的超参数调优里这也叫“lhs 初始化”常见的贝叶斯优化库多数会用 LHS 作为第一轮评估点的生成方式。一个三维问题如果直接用全因子设计需要 3^327 个点LHS 同样用 27 个点能覆盖更多水平组合如果上升到 6 维全因子设计需要 729 个点而 30 个 LHS 点就能把每维 30 个水平都覆盖这对低成本代理模型来说非常关键。采样方式边际均匀性相关性控制适用样本量典型用途rand 随机差无法控制1000MC 统计模拟lhsdesign好通过迭代优化30-500代理模型lhsu.m好默认分层排列20-200DACE 模型latin_hs.m好默认分层排列同左兼容旧项目lhsdesign 是 MATLAB 自带的优化后 LHS而 lhsu.m 更轻量不依赖并行迭代最小化。二者生成的矩阵都可以直接喂给 dacefit区别主要在于额外引入的相关性优化对 Kriging 的 theta 参数估计是否有帮助。样本量很小时过度优化 LHS 的相关性可能把采样过程中偶然出现的正相关消除掉导致后续 Kriging 评估的响应面过于平坦我一般用 lhsu 做初版模型等有第二轮补点需求时再考虑 lhsdesign。3.2 从 lhsu 到 dacefit 的一体化初始化流程只谈理论不落代码没有意义。下面是一个完整的 LHS-Kriging 初始化流程把压缩目录里两个核心函数串起来n 25; % 训练点数量 dim 3; % 输入维度 lb [-5 -5 -5]; ub [5 5 5]; % 第一步拉丁超立方采样生成均匀覆盖的设计矩阵 X lhsu(lb, ub, n); % 第二步把每一行送入黑箱仿真获得响应 Y Y zeros(n, 1); for i 1:n Y(i) expensiveSimulation(X(i, :)); end % 第三步用 dacefit 初始化 Kriging 模型 theta0 ones(1, dim) * 0.3; dmodel dacefit(X, Y, regpoly0, corrgauss, theta0, ... ones(1, dim) * 1e-3, ones(1, dim) * 50);这里的关键参数是 theta0。很多人会直接填 1但在输入尺度从 -5 到 5 时theta1 过小模型几乎对所有点都高度相关拟合出的响应面经常是一条平滑的直线。我把初始 theta0 设为 0.3配合 1e-3 到 50 的搜索边界让 dacefit 在这段区间内自己找合适的相关长度。如果响应的数值尺度差异极大最好先对 Y 做标准化避免不同量纲的响应放大数值误差。标准化的做法是对 Y 减去均值再除以标准差预测完成后按同一变换逆回去。3.3 latin_hs.m 的兼容性检查latin_hs.m 的调用方式和 lhsu 几乎一样都是生成 n×dim 矩阵区别仅在于内部随机排列的实现细节。旧项目从 latin_hs 迁移到 lhsu 时只需替换函数名后续 dacefit 的使用不需要变化。我建议新代码直接用 lhsu并给随机数生成器固定种子这样同事复现你的实验时不会因为采样路径不同得到完全不同的 Kriging 超参数。如果发现 latin_hs 在某个 MATLAB 版本下提示缺少 Statistics Toolbox不要惊慌检查它是否调用了 randperm 或 unidrnd多数情况下把外部依赖替换成 randperm 就能继续用。4. predictor.m 与 gridsamp.mKriging 预测输出与误差置信区间拟合不是终点模型要被拿去预测和寻优。DACE 将预测封装在 predictor.m 中输入是一组候选点 x 和 dmodel输出是预测值和均方根误差。gridsamp.m 则用来生成规则网格把预测结果画成云图或等值线快速看响应面的形貌。这两者在 DACE 里的关系可以理解为gridsamp 负责生成“在哪问”predictor 负责回答“值是多少、有多不确定”。4.1 predictor 的返回值与均方根误差含义predictor 的常见调用是[yhat, mse] predictor(x, dmodel)。yhat 是条件期望mse 是该点预测方差把它开方后得到的是预测标准差反映模型在该位置的置信程度。当候选点远离训练样本时mse 会显著变大这正是 Kriging 比普通多项式回归更实用的地方它能明确告诉你哪些区域还不该被信任。优化迭代里我通常不看单点的 yhat而是直接用 mse 做加点准则在 mse 最大的位置补一个样本下一次拟合后误差云图会逐渐平稳。当 theta 很大时任何两个点的相关性都趋近零mse 在非训练点位置快速抬升这就是为什么观测到 mse 云图满是“尖刺”时几乎可以断定 theta 被优化到了 upb 边界附近。4.2 gridsamp 网格密度与误差热区的匹配gridsamp 的语法是x gridsamp(range, q)range 是一个 dim×2 的矩阵每行表示该维的上下界q 是每维的网格数。它生成的坐标矩阵可以直接送给 predictor。网格密度越大响应面越细腻但计算成本也线性上升。下面是生成二维网格预测并绘制误差热区的例子% 二维输入x1 为 [0,1]x2 为 [0,1]每维 40 个网格点 range [0 1; 0 1]; xgrid gridsamp(range, 40); [yhat, mse] predictor(xgrid, dmodel); x1 reshape(xgrid(:,1), 40, 40); x2 reshape(xgrid(:,2), 40, 40); mseMap reshape(mse, 40, 40); surf(x1, x2, mseMap); % 误差高值区就是需要补充采样的位置这段代码里 40 是每维的网格点数网格太密会拖慢 surf 渲染太疏又看不出局部误差突起。网格密度选择的依据是相邻网格间距不大于 theta 相关长度的三分之一否则误差热区的边界会失真。对于 40×40 的网格predictor 单次调用只需执行几百次相关矩阵运算交互式检查很流畅。用途gridsamp 参数推荐做法快速看趋势q20只查看 yhat 等值线误差热区定位q50结合 mse 判断加样位置高维寻优不用 gridsamp用随机 LHS 候选点避免网格爆炸这里有一个常见误区直接用 gridsamp 对三维以上空间做全网格预测点数会爆炸。正确的做法是先用 lhsu 生成一批候选点再调用 predictor 找出 mse 最大或 yhat 最优的位置而不是生成完整网格。比如三维空间每维 100 点就是一百万行predictor 构建相关向量时内存占用会突然拉到上百兆Kriging 的并行优势被网格遍历抵消。4.3 用 predictor 做一次快速全局寻优既然 predictor 的单次评估足够便宜就可以在它外面套一个数值优化器搜出当前替代模型上的最优点再把这个点放到真实仿真上去验证。常见的做法是用 fmincon 从多个起点出发避免陷入单峰。lb [0 0]; ub [1 1]; f (x) predictor(x(:), dmodel); x0 lhsu(lb, ub, 10); % 多起点也来自拉丁超立方 ybest inf; for i 1:size(x0, 1) [xopt, yopt] fmincon(f, x0(i,:), [], [], [], [], lb, ub); if yopt ybest xbest xopt; ybest yopt; end end这里x(:)的作用是把列向量转成行向量因为 predictor 要求每一行是一个样本点。多起点的数量取 10 到 20 足够再多也是在同一片平滑曲面上反复移动。优化器返回的 xbest 只是当前 Kriging 代理上的最优解必须经过真实仿真确认后才能进入下一轮样本集。5. 用 data1.mat 验证 LHS-Kriging 的边界data1.mat 是压缩包里自带的数据文件通常用于快速演示 dacefit 与 predictor 的效果。加载后用 who 查看其中的变量一般会看到样本矩阵和响应向量这里按 S 和 Y 来处理实际变量名不同时只需把代码里的 S/Y 替换掉。把它当作验证集而不是训练集的常见做法是取其中一部分训练、一部分留出然后比较 yhat 与 Y。这样比只看插值误差更能说明泛化能力因为 Kriging 的插值误差在训练点上通常接近零无法反映真实预测能力。5.1 留出验证脚本load data1.mat; % 读入 S, Y rng(42); % 固定随机种子保证结果可复现 idx randperm(length(Y)); train idx(1:floor(length(Y)*0.8)); test idx(floor(length(Y)*0.8)1:end); dmodel dacefit(S(train,:), Y(train,:), regpoly0, corrgauss, ... ones(1,size(S,2))*0.2, ... ones(1,size(S,2))*1e-3, ... ones(1,size(S,2))*100); [yhat, mse] predictor(S(test,:), dmodel); rmse sqrt(mean((Y(test) - yhat).^2));这段代码里的 rng(42) 是随机种子让每次抽出的训练/测试划分保持确定。训练集占 80%测试集占 20%theta0 取遍每维 0.2 比直接取 1 更温和搜索上界 100 给大尺度解留了空间。如果 rmse 相对 Y 的量程超过 5%优先怀疑的不是相关函数选错而是样本点数量不足或 LHS 没有覆盖响应突变区域。5.2 判断模型边界的两个信号第一个信号是 theta 收敛到上界。theta 越大相关长度越短模型几乎只在样本点周围有预测能力这意味着样本间距过大Kriging 在空白区域会快速退化为回归项。第二个信号是 mse 在测试点上的均值远大于训练点。Kriging 是插值模型训练点上的 mse 通常接近零测试点上的 mse 才会反映未探索区域的不确定性。因此把 predictor 返回的 mse 当成下一轮采样的热力地图而不是简单丢弃是这套方法能持续迭代的关键。LHS-Kriging 的边界在于它假设响应是平滑且确定性的如果仿真本身含随机扰动需要先多次重复取均值再进入 dacefit。5.3 批量加点时如何更新 LHS-Kriging当 mse 热区比较分散时单点加点容易顾此失彼我会在热区附近生成一组小型 LHS 样本而不是只补一个最优位置。% 从训练数据范围缩小到热区局部范围 hotZone [0.3 0.7; 0.2 0.8]; Xadd lhsu(hotZone(:,1), hotZone(:,2), 8); Yadd expensiveSimulation(Xadd); S [S; Xadd]; Y [Y; Yadd]; dmodel dacefit(S, Y, regpoly0, corrgauss, dmodel.theta, lob, upb);这里把热区范围写死成常数便于演示实际使用中应该根据第 4 章画出的 mseMap 动态生成上下界。dmodel.theta 作为下一轮拟合的初始值能明显减少 dacefit 的迭代次数。验证完 data1.mat 之后可以把同样的流程替换成自己的黑箱仿真入口只要保证 S 每行与 Y 一一对应代码链路不需要改动。本文还有配套的精品资源点击获取
