最小二乘法曲面拟合原理
数学基础
对于离散数据点 \((x_i, y_i, z_i)\)，我们要找到曲面函数 \(z = f(x,y)\) 的最佳拟合。
多项式曲面模型：
\[z = \sum_{i=0}^{m}\sum_{j=0}^{n} a_{ij}x^iy^j
\]
其中 \(m\) 和 \(n\) 是多项式次数，\(a_{ij}\) 是待求系数。
最小二乘原理
最小二乘法通过最小化残差平方和来确定系数：
\[\min \sum_{k=1}^{N} [z_k - f(x_k,y_k)]^2
\]
实现
1. 基础曲面拟合函数
function [coefficients, fitted_surface, rmse] = surface_fit_least_squares(x, y, z, m, n)% 基于最小二乘法的曲面拟合% 输入：% x, y, z: 数据点坐标% m, n: x和y方向的拟合多项式次数% 输出：% coefficients: 拟合系数矩阵% fitted_surface: 拟合曲面函数句柄% rmse: 均方根误差% 数据预处理x = x(:); y = y(:); z = z(:);% 构建设计矩阵AA = [];for i = 0:mfor j = 0:nif i + j <= max(m, n) % 限制总阶数A = [A, (x.^i) .* (y.^j)];endendend% 求解最小二乘问题coefficients = (A' * A) \ (A' * z);% 计算拟合值z_fit = A * coefficients;% 计算误差rmse = sqrt(mean((z - z_fit).^2));% 创建拟合曲面函数fitted_surface = @(xq, yq) evaluate_surface(xq, yq, coefficients, m, n);fprintf('曲面拟合完成！\n');fprintf('多项式次数: %d×%d\n', m, n);fprintf('均方根误差(RMSE): %.6f\n', rmse);fprintf('确定系数(R²): %.6f\n', 1 - sum((z - z_fit).^2) / sum((z - mean(z)).^2));
endfunction zq = evaluate_surface(xq, yq, coeffs, m, n)% 评估拟合曲面在给定点的值xq = xq(:); yq = yq(:);Aq = [];idx = 1;for i = 0:mfor j = 0:nif i + j <= max(m, n)Aq = [Aq, (xq.^i) .* (yq.^j)];idx = idx + 1;endendendzq = Aq * coeffs;% 保持输出形状与输入一致if ~isvector(xq) || ~isvector(yq)zq = reshape(zq, size(xq));end
end
2. 不同基函数的曲面拟合
function [coefficients, fitted_surface] = surface_fit_general(x, y, z, basis_functions)% 通用基函数的最小二乘曲面拟合% basis_functions: 基函数细胞数组，例如 {@(x,y)1, @(x,y)x, @(x,y)y, @(x,y)x.*y}x = x(:); y = y(:); z = z(:);% 构建设计矩阵A = zeros(length(x), length(basis_functions));for k = 1:length(basis_functions)A(:, k) = basis_functions{k}(x, y);end% 求解系数coefficients = (A' * A) \ (A' * z);% 创建拟合函数fitted_surface = @(xq, yq) evaluate_general_surface(xq, yq, coefficients, basis_functions);
endfunction zq = evaluate_general_surface(xq, yq, coeffs, basis_funcs)xq = xq(:); yq = yq(:);zq = zeros(size(xq));for k = 1:length(coeffs)zq = zq + coeffs(k) * basis_funcs{k}(xq, yq);endif ~isvector(xq) || ~isvector(yq)zq = reshape(zq, size(xq));end
end
3. 带正则化的稳健拟合
function [coefficients, fitted_surface] = surface_fit_regularized(x, y, z, m, n, lambda)% 带正则化的曲面拟合（解决过拟合问题）% lambda: 正则化参数x = x(:); y = y(:); z = z(:);% 构建设计矩阵A = [];for i = 0:mfor j = 0:nif i + j <= max(m, n)A = [A, (x.^i) .* (y.^j)];endendend% 带正则化的最小二乘[U, S, V] = svd(A, 'econ');s = diag(S);% Tikhonov 正则化s_reg = s ./ (s.^2 + lambda^2);coefficients = V * diag(s_reg) * U' * z;% 创建拟合函数fitted_surface = @(xq, yq) evaluate_surface(xq, yq, coefficients, m, n);
end
4. 完整的演示示例
function demo_surface_fitting()% 曲面拟合演示函数% 生成示例数据[x, y, z] = generate_sample_data();% 可视化原始数据figure('Position', [100, 100, 1200, 800]);% 原始数据subplot(2, 3, 1);scatter3(x, y, z, 40, z, 'filled');title('原始离散数据');xlabel('X'); ylabel('Y'); zlabel('Z');colorbar; grid on;% 不同次数的拟合比较degrees = [1, 2, 3, 4];errors = zeros(size(degrees));for i = 1:length(degrees)m = degrees(i);n = degrees(i);subplot(2, 3, i+1);% 进行曲面拟合[coeffs, surface_func, rmse] = surface_fit_least_squares(x, y, z, m, n);errors(i) = rmse;% 生成网格用于绘制曲面[Xq, Yq] = meshgrid(linspace(min(x), max(x), 50), ...linspace(min(y), max(y), 50));Zq = surface_func(Xq, Yq);% 绘制拟合曲面surf(Xq, Yq, Zq, 'FaceAlpha', 0.7);hold on;scatter3(x, y, z, 40, 'r', 'filled');title(sprintf('%d阶多项式拟合 (RMSE=%.4f)', m, rmse));xlabel('X'); ylabel('Y'); zlabel('Z');grid on;colorbar;end% 误差比较subplot(2, 3, 6);plot(degrees, errors, 'o-', 'LineWidth', 2, 'MarkerSize', 8);xlabel('多项式次数');ylabel('RMSE');title('拟合误差 vs 多项式次数');grid on;% 显示最佳拟合参数[min_error, best_idx] = min(errors);fprintf('\n最佳拟合次数: %d阶\n', degrees(best_idx));fprintf('最小RMSE: %.6f\n', min_error);
endfunction [x, y, z] = generate_sample_data()% 生成示例数据 - 添加噪声的二次曲面rng(42); % 设置随机种子以便重现% 生成网格点[X, Y] = meshgrid(-2:0.3:2, -2:0.3:2);x = X(:); y = Y(:);% 真实的二次曲面z_true = 2 + 0.5*X - 0.8*Y + 1.2*X.^2 - 0.6*X.*Y + 0.9*Y.^2;% 添加噪声noise_level = 0.5;z = z_true(:) + noise_level * randn(size(z_true(:)));% 添加一些异常值outlier_indices = randperm(length(z), round(0.05*length(z)));z(outlier_indices) = z(outlier_indices) + 3 * noise_level * randn(size(outlier_indices));
end
5. 交互式GUI曲面拟合工具
function surface_fit_gui()% 创建GUI界面进行曲面拟合fig = figure('Name', '曲面拟合工具', ...'NumberTitle', 'off', ...'Position', [100, 100, 1400, 800]);% 控制面板uipanel('Parent', fig, ...'Title', '控制面板', ...'Position', [0.02, 0.02, 0.2, 0.96]);% 结果显示区域axes('Parent', fig, ...'Position', [0.25, 0.55, 0.35, 0.4]);title('原始数据与拟合曲面');axes('Parent', fig, ...'Position', [0.65, 0.55, 0.35, 0.4]);title('拟合残差');axes('Parent', fig, ...'Position', [0.25, 0.05, 0.35, 0.4]);title('误差分析');axes('Parent', fig, ...'Position', [0.65, 0.05, 0.35, 0.4]);title('交叉验证结果');% 创建控件create_controls(fig);
endfunction create_controls(fig)% 创建GUI控件% 多项式次数选择uicontrol('Parent', fig, ...'Style', 'text', ...'String', 'X方向次数:', ...'Position', [30, 700, 100, 20]);uicontrol('Parent', fig, ...'Style', 'popupmenu', ...'String', {'1','2','3','4','5','6'}, ...'Position', [140, 700, 60, 20], ...'Tag', 'x_degree');uicontrol('Parent', fig, ...'Style', 'text', ...'String', 'Y方向次数:', ...'Position', [30, 650, 100, 20]);uicontrol('Parent', fig, ...'Style', 'popupmenu', ...'String', {'1','2','3','4','5','6'}, ...'Position', [140, 650, 60, 20], ...'Tag', 'y_degree');% 数据生成按钮uicontrol('Parent', fig, ...'Style', 'pushbutton', ...'String', '生成示例数据', ...'Position', [30, 600, 170, 30], ...'Callback', @generate_data_callback);% 拟合按钮uicontrol('Parent', fig, ...'Style', 'pushbutton', ...'String', '执行曲面拟合', ...'Position', [30, 550, 170, 30], ...'Callback', @fit_surface_callback);% 结果显示区域uicontrol('Parent', fig, ...'Style', 'text', ...'String', '拟合结果:', ...'Position', [30, 450, 100, 20], ...'FontWeight', 'bold');uicontrol('Parent', fig, ...'Style', 'edit', ...'String', '', ...'Max', 3, ...'Position', [30, 300, 170, 150], ...'Tag', 'result_text', ...'HorizontalAlignment', 'left', ...'Style', 'listbox');
endfunction generate_data_callback(~, ~)% 生成数据回调函数[x, y, z] = generate_sample_data();% 存储数据setappdata(gcf, 'x_data', x);setappdata(gcf, 'y_data', y);setappdata(gcf, 'z_data', z);% 绘制原始数据axes_handle = findobj(gcf, 'Type', 'axes', 'Position', [0.25, 0.55, 0.35, 0.4]);axes(axes_handle);scatter3(x, y, z, 40, z, 'filled');title('原始离散数据');xlabel('X'); ylabel('Y'); zlabel('Z');colorbar; grid on;
endfunction fit_surface_callback(~, ~)% 曲面拟合回调函数% 获取数据x = getappdata(gcf, 'x_data');y = getappdata(gcf, 'y_data');z = getappdata(gcf, 'z_data');if isempty(x)errordlg('请先生成数据或加载数据！', '错误');return;end% 获取拟合参数x_degree_popup = findobj(gcf, 'Tag', 'x_degree');y_degree_popup = findobj(gcf, 'Tag', 'y_degree');m = get(x_degree_popup, 'Value');n = get(y_degree_popup, 'Value');% 执行拟合[coeffs, surface_func, rmse, r_squared] = surface_fit_least_squares(x, y, z, m, n);% 更新结果显示result_text = findobj(gcf, 'Tag', 'result_text');result_str = sprintf('拟合完成！\n多项式次数: %d×%d\nRMSE: %.6f\nR²: %.6f\n系数数量: %d', ...m, n, rmse, r_squared, length(coeffs));set(result_text, 'String', result_str);% 可视化结果visualize_fit_results(x, y, z, surface_func, coeffs, m, n);
endfunction visualize_fit_results(x, y, z, surface_func, coeffs, m, n)% 可视化拟合结果% 生成拟合曲面[Xq, Yq] = meshgrid(linspace(min(x), max(x), 50), ...linspace(min(y), max(y), 50));Zq = surface_func(Xq, Yq);% 计算拟合值z_fit = surface_func(x, y);residuals = z - z_fit;% 绘制拟合曲面axes_handle1 = findobj(gcf, 'Position', [0.25, 0.55, 0.35, 0.4]);axes(axes_handle1);cla;surf(Xq, Yq, Zq, 'FaceAlpha', 0.7, 'EdgeColor', 'none');hold on;scatter3(x, y, z, 40, 'r', 'filled');title(sprintf('%d×%d 多项式曲面拟合', m, n));xlabel('X'); ylabel('Y'); zlabel('Z');legend('拟合曲面', '原始数据', 'Location', 'best');colorbar; grid on;% 绘制残差axes_handle2 = findobj(gcf, 'Position', [0.65, 0.55, 0.35, 0.4]);axes(axes_handle2);scatter3(x, y, residuals, 40, abs(residuals), 'filled');title('拟合残差');xlabel('X'); ylabel('Y'); zlabel('残差');colorbar; grid on;% 绘制残差直方图axes_handle3 = findobj(gcf, 'Position', [0.25, 0.05, 0.35, 0.4]);axes(axes_handle3);histogram(residuals, 20);title('残差分布');xlabel('残差值'); ylabel('频数');grid on;% QQ图检验正态性axes_handle4 = findobj(gcf, 'Position', [0.65, 0.05, 0.35, 0.4]);axes(axes_handle4);qqplot(residuals);title('残差QQ图');grid on;
end
使用示例
基本使用方法
% 运行演示
demo_surface_fitting();% 或者直接使用拟合函数
[x, y, z] = generate_sample_data();
[coeffs, surface_func, rmse] = surface_fit_least_squares(x, y, z, 2, 2);% 在新点评估拟合曲面
x_new = 0.5; y_new = -0.3;
z_pred = surface_func(x_new, y_new);
fprintf('在(%.1f, %.1f)处的预测值: %.4f\n', x_new, y_new, z_pred);
启动GUI工具
% 启动交互式曲面拟合工具
surface_fit_gui();
参考代码 基于最小二乘法的离散数据的曲面拟合 www.youwenfan.com/contentcni/63812.html
关键特性
- 多种基函数支持：多项式、自定义基函数
- 正则化选项：防止过拟合
- 完整误差分析：RMSE、R²、残差分析
- 交互式可视化：3D曲面、残差图、QQ图
- 用户友好界面：GUI工具便于操作
应用建议
- 数据预处理：确保数据质量，处理异常值
- 模型选择：从低阶开始，避免过拟合
- 交叉验证：评估模型泛化能力
- 正则化：高维数据时使用正则化防止过拟合
