简介这是一份基于MATLAB实现的曲面拟合程序源码面向刚接触数值计算的新手以及有一定经验的开发人员可帮助读者快速掌握曲面拟合的基本流程与核心算法并迁移到数据可视化、参数估计等实际任务中。资源包共6个文件以5个m脚本文件为主代码按功能拆分分别负责主流程、矩阵构建、求和计算等环节且带有注释便于阅读和二次开发另有1个doc说明文档可配合源码对照理解降低上手门槛。整个压缩包仅22KB轻量易用。目前已有841人学习下载适合需要借助MATLAB完成曲面插值或拟合分析任务的读者参考。通过这份源码不仅能了解曲面拟合模型的构造思路、系数矩阵求解步骤与结果验证方法还能借鉴其模块化组织方式提高编写数值计算程序的规范性和效率减少从零编码的试错成本。1. MATLAB曲面拟合程序到底在拟合什么MATLAB 的曲面拟合解决的是这样一类问题只有一堆离散的 (x, y, z) 观测点却要在任意坐标处估计 z 的取值。地形高程、光学镜面形貌、温湿度场分布、传感器标定误差面本质上是同一个数学模型图像处理里用像素坐标和深度值重建深度图也属于这个范畴。初次拿到「曲面拟合程序源码」的人最容易把插值当成拟合。插值让曲面精确穿过观测点拟合则让曲面在整体误差最小的意义上逼近数据。一份合格的程序源码交付标准就两条系数能脱离原始数据导出任意新点都能求值模型复杂度与数据噪声之间有明确取舍而不是无脑提高多项式次数。一套能落地的源码通常由四块拼成数据读取与预处理、设计矩阵构造与系数求解、网格化求值与可视化、指标与残差评估。下面按曲面拟合的三种实现路径、手写最小二乘的完整流程、质量评估与调参、稳健扩展的顺序把整条链路走通。2. 曲面拟合的三种实现路径fit函数、手写最小二乘与自定义方程2.1 拟合与插值的边界griddata 什么时候不能用griddata做的是插值基于 Delaunay 三角剖分让曲面在每个数据点上精确穿过观测值。它的优势是不需要任何模型假设点足够密时能紧紧跟住复杂几何形状。代价是没有解析表达式、没有可导出的系数观测噪声被原样吞进曲面向外推到数据范围之外时结果完全不可控。判断场景的标准很直接数据来自确定性几何或仿真、需要逼真还原时用插值数据带测量噪声、需要参数或对未知点做泛化预测时走拟合。后者在 MATLAB 里一行就能起步F scatteredInterpolant(x, y, z)。而拟合要回答函数族、次数、求解方式和评估指标四个问题下文默认讨论的就是最小二乘意义上的回归拟合。2.2 路径一用 fit 函数做内置多项式曲面装了 Curve Fitting Toolbox 时fit是最省事的入口它把最小二乘求解、系数命名和归一化选项全部封装好。poly22表示含 6 个系数的完整二次曲面% 演示数据真实模型为二次曲面加高斯噪声 rng(7); x 2 * rand(150, 1) - 1; y 2 * rand(150, 1) - 1; z 0.3 0.8*x - 0.4*y 1.2*x.^2 - 0.6*x.*y 0.9*y.^2 0.15*randn(150, 1); % 拟合二次多项式曲面 f fit([x, y], z, poly22);fit的第一个参数把 x、y 并成一个矩阵z 必须是同长度的列向量。poly22的完整形式是z p00 p10*x p01*y p20*x^2 p11*x*y p02*y^2求解在最小二乘意义下完成返回的f是一个 cfit 对象。在任意新点上求值直接写z_new f(0.3, -0.2)配合 meshgrid 时f(Xq, Yq)一次算出整张网格前提是两个矩阵尺寸一致。系数和模型文本这样取coef coeffvalues(f); % 数值向量顺序与 coeffnames 对应 names coeffnames(f); % {p00; p10; p01; p20; p11; p02} formula(f) % 打印模型方程的文本形式poly33、poly44 依此类推指数对按「总次数递增、同一总次数内先 x 后 y」排列。工具箱额外提供 Normalize、Robust、Weights 等选项后面章节逐个展开。2.3 路径二手写设计矩阵的最小二乘多项式没有工具箱或者想完全掌控求解过程时常见做法是把多项式拟合拆成三步枚举单项式、排成设计矩阵、用反斜杠求解。对二次曲面设计矩阵就是六列n numel(z); % 基函数顺序1, x, y, x^2, x*y, y^2 A [ones(n, 1), x, y, x.^2, x.*y, y.^2]; c A \ z; % 最小二乘意义下的系数向量这个思路是社区里流传的 polyfitn 一类工具的共同内核把曲面拟合降维成线性回归。A 的每一列是一个基函数在全部样本点上的取值等式A * c ≈ z是超定方程组\自动按最小二乘求解。系数含义明确c(1) 是常数项c(2)、c(3) 是一次项c(4)、c(5)、c(6) 是三个二次项。次数提升后的通用枚举逻辑放在第三章这里先记住「为什么能拆成线性回归」这个关键认知。2.4 路径三lsqcurvefit 拟合自定义非线性方程多项式不够用、模型里必须出现指数或三角函数时走 Optimization Toolbox 的lsqcurvefit。函数句柄的第一个输入是待估参数向量 c第二个输入是自变量矩阵 XY% 自定义曲面z a * exp(-b*x^2) * cos(c*y) d fun (c, XY) c(1) * exp(-c(2) * XY(:, 1).^2) .* cos(c(3) * XY(:, 2)) c(4); c0 [1, 1, 1, 0]; % 初值对结果影响很大 c lsqcurvefit(fun, c0, [x, y], z); % 不加边界约束 % 新点求值zq fun(c, [xq, yq])句柄里 XY(:, 1) 是 x、XY(:, 2) 是 y整个表达式必须整体向量化不能写成只吃标量的形式。c0是最大变量同一个模型换初值可能收敛到完全不同的局部解稳妥做法是围绕物理量级撒 510 组初值取目标函数最小的一次。Curve Fitting Toolbox 的fittype也能定义自定义方程同样对初值敏感原理一致。三条路径的取舍按下表判断比较维度fit 内置多项式手写最小二乘lsqcurvefit方程形式poly11 到 poly55 的固定族任意整数次多项式任意自定义非线性函数工具箱依赖Curve Fitting Toolbox仅基础 MATLABOptimization Toolbox系数可解释性按 coeffnames 直接解读自定指数表对应参数有物理意义典型场景快速对比次数、GUI 探索教学、嵌入式环境二次开发指数/三角混合的物理模型主要坑高次在数据范围外震荡归一化和病态需自己处理初值敏感易陷局部最优提示三个方案返回的是同一类数学对象——让残差平方和最小的参数。差别只在函数族范围、工具箱依赖和可控程度没有绝对优劣。3. 最小二乘曲面拟合的完整实现从散点到拟合面3.1 读入 CSV 并构造任意次数的设计矩阵散点数据通常以 CSV 存放三列依次是 x、y、z。读入用readmatrix比旧的csvread对缺失值和混合列更宽容data readmatrix(surface_points.csv); % 默认从 A1 开始读全部数值 x data(:, 1); y data(:, 2); z data(:, 3);列顺序必须与文件一致文件带表头时给readmatrix加NumHeaderLines, 1。接下来的核心是通用基函数枚举。给定最高总次数 deg所有单项式满足i j degdeg 3; n numel(z); k (deg 1) * (deg 2) / 2; % 系数总数 exps zeros(k, 2); % 每行存放一个指数对 [i, j] idx 0; for d 0:deg for i 0:d idx idx 1; exps(idx, :) [i, d - i]; end end A ones(n, k); for t 1:k A(:, t) x.^exps(t, 1) .* y.^exps(t, 2); end c A \ z; % 得到长度 k 的系数向量指数枚举顺序是 (0,0)、(1,0)、(0,1)、(2,0)、(1,1)、(0,2)……即「总次数从低到高、同次数先 x 后 y」。这个顺序决定系数向量 c 的哪个位置对应哪个单项式后续预测、写报告都要以 exps 表为准。系数个数随次数增长很快总次数 deg123456系数个数 k3610152128n 远大于 k 时拟合才有统计意义。300 个点拟合到 deg415 个系数是常见配置只有几十个点却上 deg5就是在拟合噪声。3.2 求解方式与条件数为什么永远是 A\zA\z是 MATLAB 对超定最小二乘系统的推荐入口内部走带列主元的 QR 分解必要时切换 SVD。不要用教科书里的法方程形式c (A*A) \ (A*z)法方程把条件数平方本来 cond1e6 的问题会变成 1e12double 精度下有效位数所剩无几。更不要用inv(A*A) * A * zinv 从来不该出现在求解路径里。求解前扫一眼矩阵状态fprintf(cond(A) %g\n, cond(A));条件数超过 1e8 时系数的小数位基本不可信即使 R² 看着还挺高。低次数下条件数通常可控次数一高问题立刻暴露。3.3 网格求值与 surf 可视化系数求出来后在规则网格上重建拟合面并叠加原始散点[Xq, Yq] meshgrid(linspace(-1, 1, 80)); Zq zeros(size(Xq)); for t 1:k Zq Zq c(t) * Xq.^exps(t, 1) .* Yq.^exps(t, 2); end figure; surf(Xq, Yq, Zq, FaceAlpha, 0.65, EdgeColor, none); hold on; plot3(x, y, z, .r, MarkerSize, 8); xlabel(x); ylabel(y); zlabel(z); legend(拟合面, 原始散点, Location, best);meshgrid 生成的 Xq、Yq 是同尺寸矩阵循环累加时用.^逐元素幂不能写成Xq^i。FaceAlpha半透明是为了看清散点与曲面的贴合关系点很密时把散点换成scatter3(x, y, z, 6, z, filled)用颜色再叠加一层 z 值信息。3.4 高次多项式必然遇到的病态先归一化再拟合基函数 x^i y^j 在 x、y 量级不统一时列与列之间数值范围差异巨大x 在 1e3 量级时x^6 就是 1e18设计矩阵条件数指数级上涨。实测里次数从 2 涨到 5cond(A) 差出四五个数量级是常态。正确做法是先中心化再缩放把两个变量都压到 [-1, 1]mx mean(x); sx std(x); my mean(y); sy std(y); xs (x - mx) / sx; ys (y - my) / sy; % 用 xs、ys 替代 x、y 走 3.1 的构造流程 % 求值时对查询点同样做归一化 Xqs (Xq - mx) / sx; Yqs (Yq - my) / sy;归一化之后的系数对应缩放后的变量写论文报告时按x xs * sx mx反变换回物理坐标。工具箱的fit里对应选项是Normalize, on效果相同。一条经验deg 4 且数据范围很宽时不做归一化的拟合结果可以直接判定为不可信。4. 曲面拟合质量评估与调优R²、RMSE、残差与次数选择4.1 四个指标一口气算完拟合完第一件事是算指标不是看图。基于第 3 章的设计矩阵 A 和系数 cpred A * c; % 训练点回代拟合值 residual z - pred; n numel(z); k size(A, 2); % k 是参数个数 SSE residual * residual; SST (z - mean(z)) * (z - mean(z)); R2 1 - SSE / SST; RMSE sqrt(SSE / n); adjR2 1 - (1 - R2) * (n - 1) / (n - k); AIC n * log(SSE / n) 2 * k; BIC n * log(SSE / n) k * log(n);各指标的口径和判读经验指标作用判读经验R²模型解释的方差比例物理测量数据 0.95 以上算可用低于 0.9 优先怀疑缺项而非噪声调整 R²惩罚参数个数与 R² 差距明显时说明模型存在冗余项RMSE残差标准差与 z 同量纲与 std(z) 对比小于 10% 才算抓住主要信息AIC / BIC模型间比较只可比同一份数据AIC 偏预测导向BIC 偏爱简洁模型提示R² 对过拟合几乎没有分辨力次数往上涨 R² 必然单调不降。选次数要看调整 R² 和交叉验证误差别盯着 R² 那一列。4.2 残差图的三类判读模式残差等于观测值减拟合值。把残差画在三维散点上模式和成因一一对应figure; scatter3(x, y, residual, 24, residual, filled); colorbar; xlabel(x); ylabel(y); zlabel(残差);三种典型形态要能一眼分清残差沿某个方向呈弯曲带状说明缺了对应的高次项或交互项需要加次数残差随位置呈喇叭口扩大说明方差非齐性需要权重拟合大部分点残差在零附近随机分布、极少数点明显跳出去则是离群点用稳健拟合处理。配合histogram(residual, 30)看分布形态正态性检验可加lillietest(residual)p 值小于 0.05 说明残差偏离正态通常意味着模型结构有问题。4.3 用交叉验证选次数比 R² 靠谱得多把第 3 章的构造逻辑收成函数这是源码包里最常见的核心文件之一function [c, exps] fit_poly_surface(x, y, z, deg) n numel(z); k (deg 1) * (deg 2) / 2; exps zeros(k, 2); idx 0; for d 0:deg for i 0:d idx idx 1; exps(idx, :) [i, d - i]; end end A ones(n, k); for t 1:k A(:, t) x.^exps(t, 1) .* y.^exps(t, 2); end c A \ z; end留出 20% 数据做测试集其余做训练对次数从 1 到 6 遍历rng(11); idx randperm(n); ntr floor(n * 0.8); tr idx(1:ntr); te idx(ntr 1:end); for deg 1:6 [c, exps] fit_poly_surface(x(tr), y(tr), z(tr), deg); pred_te zeros(numel(te), 1); for t 1:size(exps, 1) pred_te pred_te c(t) * x(te).^exps(t, 1) .* y(te).^exps(t, 2); end rmse_val(deg) sqrt(mean((z(te) - pred_te).^2)); end [best, bestDeg] min(rmse_val); fprintf(最优次数 %d, 测试集 RMSE %.4f\n, bestDeg, best);训练集误差会随次数单调下降测试集误差先降后升最低点对应的次数就是当前数据量下的合理选择。数据集不大时把单次 8:2 切分换成重复 K 折避免切分运气左右结论。4.4 次数上界的经验值多项式曲面拟合不是次数越高越好Runge 现象在二维同样存在高次多项式在数据范围边缘剧烈震荡采样点之间出现毫无物理意义的波浪。经验值供参考光滑物理量用 deg3 或 4 收尾带明显噪声的数据封顶 deg3deg5 以上只在点密度极高、边界外推需求几乎为零时考虑。上万点或形状复杂的数据建议换径向基插值或深度学习回归——多项式曲面本身的表达能力有限硬加次数只会放大边缘误差。5. 曲面拟合进阶稳健拟合、权重与模型导出5.1 用 Bisquare 稳健拟合压制离群点数据里混入少数坏点丢帧、跳变、粗大误差时普通最小二乘会被这些点的平方残差牵着走。fit的 Robust 选项专门处理这类场景f0 fit([x, y], z, poly33); % 普通拟合 f1 fit([x, y], z, poly33, Robust, Bisquare); % 迭代加权稳健拟合Bisquare 按残差大小给样本重新分配权重残差大的点权重趋近于零LAR 以最小化绝对残差为目标对离群点更狠。对比两个模型的 RMSE 和系数变化量如果差异明显说明数据里不止一两个坏点。5.2 已知测量误差时的权重拟合各批数据测量精度不同时权重取误差方差的倒数误差大的点权重小。fit直接支持sigma 0.05 * ones(size(z)); sigma(1:40) 0.4; % 前 40 个点来自低精度设备 w 1 ./ sigma.^2; fw fit([x, y], z, poly22, Weights, w);手写实现里对应一行变换对设计矩阵 A 和观测 z 分别乘以sqrt(w)再照常反斜杠求解。权重写成相对值也成立不必是真实方差的倒数。5.3 模型落盘与预测区间让源码可以复用拟合完把模型对象存成 .mat下次直接加载不用重新读散点save(surface_fit_poly33.mat, f1, x, y, z); % 新会话里 S load(surface_fit_poly33.mat); z_new S.f1(x_new, y_new); % x_new、y_new 需同为列向量或同尺寸矩阵取系数和公式用于论文或导出给其他语言names coeffnames(f1); coef coeffvalues(f1); fprintf(%s %.6g\n, names{k}, coef(k)); % 逐项打印 formula(f1)对 cfit 对象调用predint(f1, [xq, yq])可拿到 95% 置信区间generateCode(f1)生成一段不依赖工具箱的求值代码适合交给只装基础 MATLAB 的同事。最后补一步把预测区间画出来观察区间宽度是否随查询点偏离数据重心而扩张——扩张速度夸张的模型外推预测基本没有参考价值。本文还有配套的精品资源点击获取 
长河君 · 2026/9/11 17:39:20 · 阅读 1.2k MATLAB曲面拟合实战：从最小二乘到模型评估与调参
简介这是一份基于MATLAB实现的曲面拟合程序源码面向刚接触数值计算的新手以及有一定经验的开发人员可帮助读者快速掌握曲面拟合的基本流程与核心算法并迁移到数据可视化、参数估计等实际任务中。资源包共6个文件以5个m脚本文件为主代码按功能拆分分别负责主流程、矩阵构建、求和计算等环节且带有注释便于阅读和二次开发另有1个doc说明文档可配合源码对照理解降低上手门槛。整个压缩包仅22KB轻量易用。目前已有841人学习下载适合需要借助MATLAB完成曲面插值或拟合分析任务的读者参考。通过这份源码不仅能了解曲面拟合模型的构造思路、系数矩阵求解步骤与结果验证方法还能借鉴其模块化组织方式提高编写数值计算程序的规范性和效率减少从零编码的试错成本。1. MATLAB曲面拟合程序到底在拟合什么MATLAB 的曲面拟合解决的是这样一类问题只有一堆离散的 (x, y, z) 观测点却要在任意坐标处估计 z 的取值。地形高程、光学镜面形貌、温湿度场分布、传感器标定误差面本质上是同一个数学模型图像处理里用像素坐标和深度值重建深度图也属于这个范畴。初次拿到「曲面拟合程序源码」的人最容易把插值当成拟合。插值让曲面精确穿过观测点拟合则让曲面在整体误差最小的意义上逼近数据。一份合格的程序源码交付标准就两条系数能脱离原始数据导出任意新点都能求值模型复杂度与数据噪声之间有明确取舍而不是无脑提高多项式次数。一套能落地的源码通常由四块拼成数据读取与预处理、设计矩阵构造与系数求解、网格化求值与可视化、指标与残差评估。下面按曲面拟合的三种实现路径、手写最小二乘的完整流程、质量评估与调参、稳健扩展的顺序把整条链路走通。2. 曲面拟合的三种实现路径fit函数、手写最小二乘与自定义方程2.1 拟合与插值的边界griddata 什么时候不能用griddata做的是插值基于 Delaunay 三角剖分让曲面在每个数据点上精确穿过观测值。它的优势是不需要任何模型假设点足够密时能紧紧跟住复杂几何形状。代价是没有解析表达式、没有可导出的系数观测噪声被原样吞进曲面向外推到数据范围之外时结果完全不可控。判断场景的标准很直接数据来自确定性几何或仿真、需要逼真还原时用插值数据带测量噪声、需要参数或对未知点做泛化预测时走拟合。后者在 MATLAB 里一行就能起步F scatteredInterpolant(x, y, z)。而拟合要回答函数族、次数、求解方式和评估指标四个问题下文默认讨论的就是最小二乘意义上的回归拟合。2.2 路径一用 fit 函数做内置多项式曲面装了 Curve Fitting Toolbox 时fit是最省事的入口它把最小二乘求解、系数命名和归一化选项全部封装好。poly22表示含 6 个系数的完整二次曲面% 演示数据真实模型为二次曲面加高斯噪声 rng(7); x 2 * rand(150, 1) - 1; y 2 * rand(150, 1) - 1; z 0.3 0.8*x - 0.4*y 1.2*x.^2 - 0.6*x.*y 0.9*y.^2 0.15*randn(150, 1); % 拟合二次多项式曲面 f fit([x, y], z, poly22);fit的第一个参数把 x、y 并成一个矩阵z 必须是同长度的列向量。poly22的完整形式是z p00 p10*x p01*y p20*x^2 p11*x*y p02*y^2求解在最小二乘意义下完成返回的f是一个 cfit 对象。在任意新点上求值直接写z_new f(0.3, -0.2)配合 meshgrid 时f(Xq, Yq)一次算出整张网格前提是两个矩阵尺寸一致。系数和模型文本这样取coef coeffvalues(f); % 数值向量顺序与 coeffnames 对应 names coeffnames(f); % {p00; p10; p01; p20; p11; p02} formula(f) % 打印模型方程的文本形式poly33、poly44 依此类推指数对按「总次数递增、同一总次数内先 x 后 y」排列。工具箱额外提供 Normalize、Robust、Weights 等选项后面章节逐个展开。2.3 路径二手写设计矩阵的最小二乘多项式没有工具箱或者想完全掌控求解过程时常见做法是把多项式拟合拆成三步枚举单项式、排成设计矩阵、用反斜杠求解。对二次曲面设计矩阵就是六列n numel(z); % 基函数顺序1, x, y, x^2, x*y, y^2 A [ones(n, 1), x, y, x.^2, x.*y, y.^2]; c A \ z; % 最小二乘意义下的系数向量这个思路是社区里流传的 polyfitn 一类工具的共同内核把曲面拟合降维成线性回归。A 的每一列是一个基函数在全部样本点上的取值等式A * c ≈ z是超定方程组\自动按最小二乘求解。系数含义明确c(1) 是常数项c(2)、c(3) 是一次项c(4)、c(5)、c(6) 是三个二次项。次数提升后的通用枚举逻辑放在第三章这里先记住「为什么能拆成线性回归」这个关键认知。2.4 路径三lsqcurvefit 拟合自定义非线性方程多项式不够用、模型里必须出现指数或三角函数时走 Optimization Toolbox 的lsqcurvefit。函数句柄的第一个输入是待估参数向量 c第二个输入是自变量矩阵 XY% 自定义曲面z a * exp(-b*x^2) * cos(c*y) d fun (c, XY) c(1) * exp(-c(2) * XY(:, 1).^2) .* cos(c(3) * XY(:, 2)) c(4); c0 [1, 1, 1, 0]; % 初值对结果影响很大 c lsqcurvefit(fun, c0, [x, y], z); % 不加边界约束 % 新点求值zq fun(c, [xq, yq])句柄里 XY(:, 1) 是 x、XY(:, 2) 是 y整个表达式必须整体向量化不能写成只吃标量的形式。c0是最大变量同一个模型换初值可能收敛到完全不同的局部解稳妥做法是围绕物理量级撒 510 组初值取目标函数最小的一次。Curve Fitting Toolbox 的fittype也能定义自定义方程同样对初值敏感原理一致。三条路径的取舍按下表判断比较维度fit 内置多项式手写最小二乘lsqcurvefit方程形式poly11 到 poly55 的固定族任意整数次多项式任意自定义非线性函数工具箱依赖Curve Fitting Toolbox仅基础 MATLABOptimization Toolbox系数可解释性按 coeffnames 直接解读自定指数表对应参数有物理意义典型场景快速对比次数、GUI 探索教学、嵌入式环境二次开发指数/三角混合的物理模型主要坑高次在数据范围外震荡归一化和病态需自己处理初值敏感易陷局部最优提示三个方案返回的是同一类数学对象——让残差平方和最小的参数。差别只在函数族范围、工具箱依赖和可控程度没有绝对优劣。3. 最小二乘曲面拟合的完整实现从散点到拟合面3.1 读入 CSV 并构造任意次数的设计矩阵散点数据通常以 CSV 存放三列依次是 x、y、z。读入用readmatrix比旧的csvread对缺失值和混合列更宽容data readmatrix(surface_points.csv); % 默认从 A1 开始读全部数值 x data(:, 1); y data(:, 2); z data(:, 3);列顺序必须与文件一致文件带表头时给readmatrix加NumHeaderLines, 1。接下来的核心是通用基函数枚举。给定最高总次数 deg所有单项式满足i j degdeg 3; n numel(z); k (deg 1) * (deg 2) / 2; % 系数总数 exps zeros(k, 2); % 每行存放一个指数对 [i, j] idx 0; for d 0:deg for i 0:d idx idx 1; exps(idx, :) [i, d - i]; end end A ones(n, k); for t 1:k A(:, t) x.^exps(t, 1) .* y.^exps(t, 2); end c A \ z; % 得到长度 k 的系数向量指数枚举顺序是 (0,0)、(1,0)、(0,1)、(2,0)、(1,1)、(0,2)……即「总次数从低到高、同次数先 x 后 y」。这个顺序决定系数向量 c 的哪个位置对应哪个单项式后续预测、写报告都要以 exps 表为准。系数个数随次数增长很快总次数 deg123456系数个数 k3610152128n 远大于 k 时拟合才有统计意义。300 个点拟合到 deg415 个系数是常见配置只有几十个点却上 deg5就是在拟合噪声。3.2 求解方式与条件数为什么永远是 A\zA\z是 MATLAB 对超定最小二乘系统的推荐入口内部走带列主元的 QR 分解必要时切换 SVD。不要用教科书里的法方程形式c (A*A) \ (A*z)法方程把条件数平方本来 cond1e6 的问题会变成 1e12double 精度下有效位数所剩无几。更不要用inv(A*A) * A * zinv 从来不该出现在求解路径里。求解前扫一眼矩阵状态fprintf(cond(A) %g\n, cond(A));条件数超过 1e8 时系数的小数位基本不可信即使 R² 看着还挺高。低次数下条件数通常可控次数一高问题立刻暴露。3.3 网格求值与 surf 可视化系数求出来后在规则网格上重建拟合面并叠加原始散点[Xq, Yq] meshgrid(linspace(-1, 1, 80)); Zq zeros(size(Xq)); for t 1:k Zq Zq c(t) * Xq.^exps(t, 1) .* Yq.^exps(t, 2); end figure; surf(Xq, Yq, Zq, FaceAlpha, 0.65, EdgeColor, none); hold on; plot3(x, y, z, .r, MarkerSize, 8); xlabel(x); ylabel(y); zlabel(z); legend(拟合面, 原始散点, Location, best);meshgrid 生成的 Xq、Yq 是同尺寸矩阵循环累加时用.^逐元素幂不能写成Xq^i。FaceAlpha半透明是为了看清散点与曲面的贴合关系点很密时把散点换成scatter3(x, y, z, 6, z, filled)用颜色再叠加一层 z 值信息。3.4 高次多项式必然遇到的病态先归一化再拟合基函数 x^i y^j 在 x、y 量级不统一时列与列之间数值范围差异巨大x 在 1e3 量级时x^6 就是 1e18设计矩阵条件数指数级上涨。实测里次数从 2 涨到 5cond(A) 差出四五个数量级是常态。正确做法是先中心化再缩放把两个变量都压到 [-1, 1]mx mean(x); sx std(x); my mean(y); sy std(y); xs (x - mx) / sx; ys (y - my) / sy; % 用 xs、ys 替代 x、y 走 3.1 的构造流程 % 求值时对查询点同样做归一化 Xqs (Xq - mx) / sx; Yqs (Yq - my) / sy;归一化之后的系数对应缩放后的变量写论文报告时按x xs * sx mx反变换回物理坐标。工具箱的fit里对应选项是Normalize, on效果相同。一条经验deg 4 且数据范围很宽时不做归一化的拟合结果可以直接判定为不可信。4. 曲面拟合质量评估与调优R²、RMSE、残差与次数选择4.1 四个指标一口气算完拟合完第一件事是算指标不是看图。基于第 3 章的设计矩阵 A 和系数 cpred A * c; % 训练点回代拟合值 residual z - pred; n numel(z); k size(A, 2); % k 是参数个数 SSE residual * residual; SST (z - mean(z)) * (z - mean(z)); R2 1 - SSE / SST; RMSE sqrt(SSE / n); adjR2 1 - (1 - R2) * (n - 1) / (n - k); AIC n * log(SSE / n) 2 * k; BIC n * log(SSE / n) k * log(n);各指标的口径和判读经验指标作用判读经验R²模型解释的方差比例物理测量数据 0.95 以上算可用低于 0.9 优先怀疑缺项而非噪声调整 R²惩罚参数个数与 R² 差距明显时说明模型存在冗余项RMSE残差标准差与 z 同量纲与 std(z) 对比小于 10% 才算抓住主要信息AIC / BIC模型间比较只可比同一份数据AIC 偏预测导向BIC 偏爱简洁模型提示R² 对过拟合几乎没有分辨力次数往上涨 R² 必然单调不降。选次数要看调整 R² 和交叉验证误差别盯着 R² 那一列。4.2 残差图的三类判读模式残差等于观测值减拟合值。把残差画在三维散点上模式和成因一一对应figure; scatter3(x, y, residual, 24, residual, filled); colorbar; xlabel(x); ylabel(y); zlabel(残差);三种典型形态要能一眼分清残差沿某个方向呈弯曲带状说明缺了对应的高次项或交互项需要加次数残差随位置呈喇叭口扩大说明方差非齐性需要权重拟合大部分点残差在零附近随机分布、极少数点明显跳出去则是离群点用稳健拟合处理。配合histogram(residual, 30)看分布形态正态性检验可加lillietest(residual)p 值小于 0.05 说明残差偏离正态通常意味着模型结构有问题。4.3 用交叉验证选次数比 R² 靠谱得多把第 3 章的构造逻辑收成函数这是源码包里最常见的核心文件之一function [c, exps] fit_poly_surface(x, y, z, deg) n numel(z); k (deg 1) * (deg 2) / 2; exps zeros(k, 2); idx 0; for d 0:deg for i 0:d idx idx 1; exps(idx, :) [i, d - i]; end end A ones(n, k); for t 1:k A(:, t) x.^exps(t, 1) .* y.^exps(t, 2); end c A \ z; end留出 20% 数据做测试集其余做训练对次数从 1 到 6 遍历rng(11); idx randperm(n); ntr floor(n * 0.8); tr idx(1:ntr); te idx(ntr 1:end); for deg 1:6 [c, exps] fit_poly_surface(x(tr), y(tr), z(tr), deg); pred_te zeros(numel(te), 1); for t 1:size(exps, 1) pred_te pred_te c(t) * x(te).^exps(t, 1) .* y(te).^exps(t, 2); end rmse_val(deg) sqrt(mean((z(te) - pred_te).^2)); end [best, bestDeg] min(rmse_val); fprintf(最优次数 %d, 测试集 RMSE %.4f\n, bestDeg, best);训练集误差会随次数单调下降测试集误差先降后升最低点对应的次数就是当前数据量下的合理选择。数据集不大时把单次 8:2 切分换成重复 K 折避免切分运气左右结论。4.4 次数上界的经验值多项式曲面拟合不是次数越高越好Runge 现象在二维同样存在高次多项式在数据范围边缘剧烈震荡采样点之间出现毫无物理意义的波浪。经验值供参考光滑物理量用 deg3 或 4 收尾带明显噪声的数据封顶 deg3deg5 以上只在点密度极高、边界外推需求几乎为零时考虑。上万点或形状复杂的数据建议换径向基插值或深度学习回归——多项式曲面本身的表达能力有限硬加次数只会放大边缘误差。5. 曲面拟合进阶稳健拟合、权重与模型导出5.1 用 Bisquare 稳健拟合压制离群点数据里混入少数坏点丢帧、跳变、粗大误差时普通最小二乘会被这些点的平方残差牵着走。fit的 Robust 选项专门处理这类场景f0 fit([x, y], z, poly33); % 普通拟合 f1 fit([x, y], z, poly33, Robust, Bisquare); % 迭代加权稳健拟合Bisquare 按残差大小给样本重新分配权重残差大的点权重趋近于零LAR 以最小化绝对残差为目标对离群点更狠。对比两个模型的 RMSE 和系数变化量如果差异明显说明数据里不止一两个坏点。5.2 已知测量误差时的权重拟合各批数据测量精度不同时权重取误差方差的倒数误差大的点权重小。fit直接支持sigma 0.05 * ones(size(z)); sigma(1:40) 0.4; % 前 40 个点来自低精度设备 w 1 ./ sigma.^2; fw fit([x, y], z, poly22, Weights, w);手写实现里对应一行变换对设计矩阵 A 和观测 z 分别乘以sqrt(w)再照常反斜杠求解。权重写成相对值也成立不必是真实方差的倒数。5.3 模型落盘与预测区间让源码可以复用拟合完把模型对象存成 .mat下次直接加载不用重新读散点save(surface_fit_poly33.mat, f1, x, y, z); % 新会话里 S load(surface_fit_poly33.mat); z_new S.f1(x_new, y_new); % x_new、y_new 需同为列向量或同尺寸矩阵取系数和公式用于论文或导出给其他语言names coeffnames(f1); coef coeffvalues(f1); fprintf(%s %.6g\n, names{k}, coef(k)); % 逐项打印 formula(f1)对 cfit 对象调用predint(f1, [xq, yq])可拿到 95% 置信区间generateCode(f1)生成一段不依赖工具箱的求值代码适合交给只装基础 MATLAB 的同事。最后补一步把预测区间画出来观察区间宽度是否随查询点偏离数据重心而扩张——扩张速度夸张的模型外推预测基本没有参考价值。本文还有配套的精品资源点击获取