问题：
给定一组坐标(xg,yg,zg)，g?1,2,…,n，表示有n个点。要求用以下二元多项式函数对所给的坐标进行拟合：
p,qpqijf(x,y)??ij?1,1aijxi?1yj?1???ai?1j?1xi?1yj?1
即
f(x,y)?a11?a21x??ai1x??ap1xp?1i?1?a12y?a22xy?ai2x?ap2xi?1?a13y22????2?a1qyq?1q?1?a23xyy?ai3x?ap3xi?1?a2qxy?aiqx?apqxi?1y??yq?1
p?1yp?1y2??p?1yq?1设
?1??1??a11????xy?????a2122?????x?x,y?y,A???????????????ap1pq?????xy????a12a22?ap2????a1q??a2q? ???apq??则函数又可表示为f(x,y)?xTAy，拟合的目标就是求出系数矩阵A。
最小二乘法：
构造关于系数aij的多元函数：
nngpgqijs(a11,?,apq)???g?1[f(xg,yg)?zg]?2??(??ag?1i?1j?1xi?1yj?1?zg)
2点（a11，…，apq）是多元函数s(a11,?,apq)的极小点，其中?g为权函数，默认为1，所以点（a11，…，apq）必须满足方程组
?s?aij?0
在?g?1的情况下，有
?s?aij???aijnn?[f(xg?1g,yg)?zg]2????????2[f(xg,yg)?zg][f(xg,yg)]??aijg?1???????2[f(xg?1ng?1n
g,yg)?zg]xgygi?1j?1?i?1j?1i?1j?1?2??xyf(x,y)?xyzg?gggg?gg?因此可得
nni?1g?xg?1nyj?1gf(xg,yg)?pq?xg?1i?1gygzg
nj?1?xg?1i?1gyj?1g??a??x??1??1p,q??1gyg??1??xg?1i?1gygzg
j?1nn?xg?1p,qi?1gyj?1g????1,1na??xg??1yg??1??xg?1i?1gygzg
j?1???a?????1,1??(xg?1??1gygx??1i?1gyj?1g?)???n?xg?1i?1gygzg
j?1令
nnu??(i,j)??(xg?1??1gyg??1xi?1gyj?1g)，v(i,j)??xg?1i?1gygzg
j?1则
p,q????1,1a??u??(i,j)?v(i,j) (i,j)?(1，1),…,(p,q)
上式实际共有p?q个等式，可将这p?q个等式写成矩阵的形式有：
?u11(1,1)????u(p,q)?11???upq(1,1)????upq(p,q)???a11??v(1,1)???????? ?????a??v(p,q)???pq??也就是U*a=V的形式，其中
?u11(1,1)?U????u(p,q)?11???upq(1,1)????upq(p,q)???a11???，a?????a??pq??v(1,1)???，V????
??v(p,q)??U为pq?pq阶矩阵，实现函数为function A=leftmatrix(x,p,y,q)；V为长pq的列向量，实现函数为function B=rightmatrix(x,p,y,q,z)。这样就可以算出列矩阵a,然后转化成A。
例子：
某地区有一煤矿，为估计其储量以便于开采，先在该地区进行勘探。假设该地区是一长方形区域，长为4公里，宽为5公里。经勘探得到如下数据：
煤矿勘探数据表 编号 横坐标（公里） 纵坐标（公里） 煤层厚度（米） 1 1 1 13.72 11 3 1 23.28 2 1 2 25.80 12 3 2 26.48 3 1 3 8.47 13 3 3 29.14 4 1 4 25.27 14 3 4 12.04 5 1 5 22.32 15 3 5 14.58 6 2 1 15.47 16 4 1 19.95 7 2 2 21.33 17 4 2 23.73 8 2 3 14.49 18 4 3 15.35 9 2 4 24.83 19 4 4 18.01 10 2 5 26.19 20 4 5 16.29 编号 横坐标（公里） 纵坐标（公里） 煤层厚度（米） 请你估计出此地区内（2?x?4,并用电脑画出该煤矿的三维图象。
如果直接画出三维曲面图形：
1?y?5）煤的储量，单位用立方米表示，
clear;x=1:4;y=1:5;[X,Y]=meshgrid(x,y) Z=[13.72 25.80 8.47 25.27 22.32; 15.47 21.33 14.49 24.83 26.19; 23.28 26.48 29.14 12.04 14.58; 19.95 23.73 15.35 18.01 16.29]' surf(X,Y,Z); X =
1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 Y =
1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 Z =
13.7200 15.4700 23.2800 19.9500
25.8000 21.3300 26.4800 23.7300 8.4700 14.4900 29.1400 15.3500 25.2700 24.8300 12.0400 18.0100 22.3200 26.1900 14.5800 16.2900
30252015105543211234
粗略计算体积：底面积乘以平均高度。
p=sum(Z); q=p(:,[2,3,4]); h=sum(q')/15 v=2000*4000*h h =
20.0773 v =
1.6062e+008
进行线性插值：
xi=linspace(1,4,31);yi=linspace(1,5,41);[XI,YI]=meshgrid(xi,yi); ZI=interp2(X,Y,Z,XI,YI,'linear'); surf(XI,YI,ZI);
进行三次多项式插值：
xi=linspace(1,4,31);yi=linspace(1,5,41);[XI,YI]=meshgrid(xi,yi); ZI=interp2(X,Y,Z,XI,YI,'cubic'); surf(XI,YI,ZI);
进行插值后计算体积：底面积乘以平均高度。
曲面拟合原理与实例
