基于最小二乘法的离散数据曲面拟合MATLAB实现方法

最小二乘法曲面拟合原理

数学基础

对于离散数据点 \((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(:);
    
    % 构建设计矩阵A
    A = [];
    for i = 0:m
        for j = 0:n
            if i + j <= max(m, n)  % 限制总阶数
                A = [A, (x.^i) .* (y.^j)];
            end
        end
    end
    
    % 求解最小二乘问题
    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));
end

function zq = evaluate_surface(xq, yq, coeffs, m, n)
    % 评估拟合曲面在给定点的值
    xq = xq(:); yq = yq(:);
    Aq = [];
    idx = 1;
    
    for i = 0:m
        for j = 0:n
            if i + j <= max(m, n)
                Aq = [Aq, (xq.^i) .* (yq.^j)];
                idx = idx + 1;
            end
        end
    end
    
    zq = 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);
end

function 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);
    end
    
    if ~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:m
        for j = 0:n
            if i + j <= max(m, n)
                A = [A, (x.^i) .* (y.^j)];
            end
        end
    end
    
    % 带正则化的最小二乘
    [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);
end

function [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);
end

function 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');
end

function 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;
end

function 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);
end

function 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

关键特性

1. 多种基函数支持：多项式、自定义基函数
2. 正则化选项：防止过拟合
3. 完整误差分析：RMSE、R²、残差分析
4. 交互式可视化：3D曲面、残差图、QQ图
5. 用户友好界面：GUI工具便于操作

应用建议

• 数据预处理：确保数据质量，处理异常值
• 模型选择：从低阶开始，避免过拟合
• 交叉验证：评估模型泛化能力
• 正则化：高维数据时使用正则化防止过拟合