曲面拟合原理与实例.doc
问题：给定一组坐标，表示有n个点。要求用以下二元多项式函数对所给的坐标进行拟合：即设则函数又可表示为，拟合的目标就是求出系数矩阵A。最小二乘法：构造关于系数的多元函数：点（，）是多元函数的极小点，其中为权函数，默认为1，所以点（，）必须满足方程组 在的情况下，有因此可得令，则 上式实际共有个等式，可将这个等式写成矩阵的形式有：也就是U*a=V的形式，其中，U为阶矩阵，实现函数为function A=leftmatrix(x,p,y,q)；V为长的列向量，实现函数为function B=rightmatrix(x,p,y,q,z)。这样就可以算出列矩阵a,然后转化成A。例子：某地区有一煤矿，为估计其储量以便于开采，先在该地区进行勘探。假设该地区是一长方形区域，长为4公里，宽为5公里。经勘探得到如下数据：煤矿勘探数据表编号12345678910横坐标（公里）1111122222纵坐标（公里）1234512345煤层厚度（米）13.7225.808.4725.2722.3215.4721.3314.4924.8326.19编号11121314151617181920横坐标（公里）3333344444纵坐标（公里）1234512345煤层厚度（米）23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29请你估计出此地区内（）煤的储量，单位用立方米表示，并用电脑画出该煤矿的三维图象。如果直接画出三维曲面图形：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 4Y = 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5Z = 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 粗略计算体积：底面积乘以平均高度。p=sum(Z);q=p(:,2,3,4);h=sum(q')/15v=2000*4000*h h = 20.0773v = 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); 进行插值后计算体积：底面积乘以平均高度。xi=linspace(1,4,61);yi=linspace(1,5,81);XI,YI=meshgrid(xi,yi);ZI=interp2(X,Y,Z,XI,YI,'cubic'); surf(XI,YI,ZI); H=0;n=0;for j=21:61for i=1:81H=H+ZI(i,j);n=n+1;endendnH=H/nS=2000*4000;V=S*H n = 3321H = 20.8222V = 1.6658e+008 上面是插值的方法解题，下面用拟合的方法解题。为此编写了几个M函数：leftmatrix.mfunction U=leftmatrix(x,p,y,q)% U*a=V a为系数列矩阵，长度为p*q% U为左边p*q乘p*q矩阵% x,y 为长度一致的列矩阵，给定点的坐标% p,q 为拟合的函数中x,y的幂的最高次数 m=length(x);if (nargin=4) & (m=length(y) error('error check check!');end U_length=p*q; % U 为p*q阶方阵U=zeros(U_length,U_length); % 赋值0，目的是分配内存for i=1 : p*q for j= 1 : p*q x_z=quotient(j-1,q)+quotient(i-1,q); % x 的幂的次数，quotient为求商 y_z=mod(j-1,q)+mod(i-1,q); % y 的幂的次数 U(i,j)=qiuhe(x,x_z,y,y_z); endendrightmatrix.mfunction V=rightmatrix(x,p,y,q,z)% U*a=V % V 为一个列向量 长为p*q% x y z 为点的坐标 %p q 分别为x y幂的最高次数 if nargin=5 error('error check check! rightmatrix')end V=zeros(p*q,1);for i=1 : p*q x_z=quotient(i-1,q); y_z=mod(i-1,q); V(i,1)=qiuhe(x,x_z,y,y_z,z);endquuotient.mfunction sh=quotient(x,y)% sh 为 x/y 的商 sh=(x-mod(x,y)/y;qiuhe.mfunction he=qiuhe(x,p,y,q,z)% he xp*yq 从1>m的和% x,y 向量 长度相同% p,q分别为x,y的幂的次数m=length(x);if (nargin<4 )&(m=length(y) %输入量至少为四，x,y行向量长度必需一样 error('error check check!');endif nargin=4 %没有 z , 默认为元素全部为1的向量 z=ones(m,1);end he=0;for i=1:m he=he+x(i)p * y(i)q*z(i); % 1->m 求和end下面一段程序先进行拟合，然后验证拟合的效果，具体操作：先输入x=y=z=p=q= (注意x,y,z是向量)；拟合得到系数a，也就是得到了拟合的函数；根据拟合函数计算给定点(xx, yy)的函数值zz=f(xx, yy)并进行画图检验。程序保存于M文件fit.m。fit.mclear;X,Y=meshgrid(1:4,1:5);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'x=reshape(X,20,1);y=reshape(Y,20,1);z=reshape(Z,20,1);p=4;q=5;U=leftmatrix(x,p,y,q); % U*a_n=VV=rightmatrix(x,p,y,q,z); %a_n=inv(U)*V;a_n=UV; for i=1 : length(a_n) % 把长为p*q 的列向量a_n转换成p*q的矩阵aa ii=quotient(i-1,q)+1; % quotient求商 jj=mod(i-1,q)+1; aa(ii,jj)=a_n(i,1);endaam=31;n=41;%m=4;n=5;XI,YI=meshgrid(linspace(1,4,m),linspace(1,5,n);xx=reshape(XI,m*n,1);yy=reshape(YI,m*n,1); zz=zeros(m*n,1); xy=zeros(m*n,1); xt=zeros(m*n,1); yt=zeros(m*n,1);%zz=0; % zz 是 xx,yy 代入所拟合的函数 求出的函数值for i=1 : p % 函数为aa(i,j)*xi*yj，(i=1.p,j=1.q)for j=1 : q % aa 为pxq的系数的矩阵 xt=xx.(i-1); yt=yy.(j-1); xy=xt.*yt; zz=zz+aa(i,j).*xy; endend ZI=reshape(zz,n,m);surf(XI,YI,ZI); %axis(1 4 1 5 0 30) aa = 1.0e+003 * 0.1465 -0.2678 0.2132 -0.0624 0.0058 -0.7287 1.3972 -0.9275 0.2412 -0.0210 0.4416 -0.8415 0.5487 -0.1407 0.0122 -0.0680 0.1295 -0.0839 0.0214 -0.0018 注意：权函数在拟合的函数非常重要，不过她只能按你遇到的具体问题来取，我这里为1；当p，q越大时，拟合的函数与原数据的方差越小，但是有可能函数本身抖动非常厉害，可以画图看出来。
