这里用到的还是最小二乘方法,和上一次这篇文章原理差不多。

就是首先构造最小二乘函数,然后对每一个系数计算偏导,构造矩阵乘法形式,最后解方程组。

比如有一个二次曲面:z=ax^2+by^2+cxy+dx+ey+f

首先构造最小二乘函数,然后计算系数偏导(我直接手写了):

解方程组(下图中A矩阵后面求和符号我就没写了啊),然后计算C:

代码如下:

 1 clear all; 2 close all; 3 clc; 4  5 a=2;b=2;c=-3;d=1;e=2;f=30;   %系数          6 n=1:0.2:20; 7 x=repmat(n,96,1); 8 y=repmat(n',1,96); 9 z=a*x.^2+b*y.^2+c*x.*y+d*x+e*y +f;      %原始模型     10 surf(x,y,z)11 12 N=100;13 ind=int8(rand(N,2)*95+1);14 15 X=x(sub2ind(size(x),ind(:,1),ind(:,2)));16 Y=y(sub2ind(size(y),ind(:,1),ind(:,2)));17 Z=z(sub2ind(size(z),ind(:,1),ind(:,2)))+rand(N,1)*20;       %生成待拟合点,加个噪声18 19 hold on;20 plot3(X,Y,Z,'o');21 22 A=[N sum(Y) sum(X) sum(X.*Y) sum(Y.^2) sum(X.^2);23    sum(Y) sum(Y.^2) sum(X.*Y) sum(X.*Y.^2) sum(Y.^3) sum(X.^2.*Y);24    sum(X) sum(X.*Y) sum(X.^2) sum(X.^2.*Y) sum(X.*Y.^2) sum(X.^3);25    sum(X.*Y) sum(X.*Y.^2) sum(X.^2.*Y) sum(X.^2.*Y.^2) sum(X.*Y.^3) sum(X.^3.*Y);26    sum(Y.^2) sum(Y.^3) sum(X.*Y.^2) sum(X.*Y.^3) sum(Y.^4) sum(X.^2.*Y.^2);27    sum(X.^2) sum(X.^2.*Y) sum(X.^3) sum(X.^3.*Y) sum(X.^2.*Y.^2) sum(X.^4)];28 29 B=[sum(Z) sum(Z.*Y) sum(Z.*X) sum(Z.*X.*Y) sum(Z.*Y.^2) sum(Z.*X.^2)]';30 31 C=inv(A)*B;32 33 z=C(6)*x.^2+C(5)*y.^2+C(4)*x.*y+C(3)*x+C(2)*y +C(1);           %拟合结果34 35 mesh(x,y,z)

结果如下,深色曲面是原模型,浅色曲面是用噪声数据拟合的模型: