曲面拟合原理与实例

1、问题:给定一组坐标(xgygZg) , g 1,2,n ,表示有n个点。要求用以下二元0多项式函数对所给的坐标进行拟合:f(x,y)i 1 j3j x yij 1,11aii 1 j 1i 1 j 1 iX yf (x, y)a11a12y2a13yLq 1a21xa22xy2a23xyLq 1a2qxyM1i 1i 12Li 1 q 1ai1xy33x yaiqXyqMp 1p 1p 12Lp 1 qap1Xap2X yap3XyapqX yp,qpq即1x2x xMxp,yy2yMyq,Aa12La1qa22La2qMOMap2Lapqaiia21Map1则函数又可表示为f (x, y)x

2、TAy ,拟合的目标就是求出系数矩阵 A。最小二乘法:构造关于系数aj的多元函数:p q/i 1 j 1g(ajX yi 1 j 1nzg)22s( ai1，L，apq)gf(xg，yg)Zg g 1点（311，apq）是多元函数s（a11,L ,apq）的极小点，其中 g为权函数，默认为1,所以点（811，apq ）必须满足方程组s3ij在g 1的情况下，有aijaij2f (Xg,yg)Zgg i2f(Xg,yg)叩石If)2f (Xg,yg)i 1 j 1ZgXg ygg inxg1yg 1f(xg,yg) xg+gz2g i因此可得u (i, j)n(Xgg 11yg1 i 1Xgyg

3、1), v(i,j)ni 1 j 1Xg yg Zgg 1nni 1 j 1i 1 j 1Xg ygf(Xg,yg)Xg yg Zgg 1g 1np qniXg1yg1 1a Xg yg1i 1 jXg Yg1Zgg 11 1g 1np,qniXg1yg1 1 1aXgygi 1 j 1Xg ygzgg 11,1g 1p,qnnai 1 j(xg yg Xg yg1、i 1)Xgyg 1zg1,1g 1g 1p,qa u (i, j) v(i, j) (i, j)(1,1),(p,q)1,1上式实际共有p q个等式，可将这比1(1,1) LUpq(1,1)anMO M MUn(p,q) L U

4、pq(p,q) apq也就是U*a=V的形式，其中Un(1,1)LUpq (1,1)UMOM，Un(P,q)LUpq(p,q)p q个等式写成矩阵的形式有:v(1,1)Mv( p, q)anv(1,1)a M ，V Mapqv(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公里。经勘探

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请你估计出此地区内（2x4,1 y 5 ）煤的储量，单位用立方米表示,并用电脑画出该煤矿的三维图象。如果直接画出三维曲面图形:clear;

6、x=1:4;y=1:5;X,Y=meshgrid(x,y)Z=13.7225.80 8.4725.27 22.32;15.47 21.3314.49 24.83 26.19;23.28 26.48 29.14 12.0414.58;19.95 23.73 15.35 18.0116.29surf(X ,Y,Z);X =12341234123412341234Y =11112222333344445555Z =13.7200 15.4700 23.2800 19.950025.800021.330026.480023.73008.470014.490029.140015.350025.27002

7、4.830012.040018.010022.320026.190014.580016.2900粗略计算体积：底面积乘以平均高度。p=sum(Z);q=P(：,2,3,4);h=sum(q)/15v=2000*4000*hh =20.0773v =1.6062e+008进行线性插值：xi=li nspace(1,4,31);yi=li nspace(1,5,41);XI,YI=meshgrid(xi,yi);ZI=i nterp2(X,Y,Z,XI,YI,li near);surf(XI,YI,ZI);1 1进行三次多项式插值：41);XI,YI=meshgrid(xi,yi);xi=li n

8、space(1,4,31);yi=li nspace(1,5,Zl=interp2(X,Y,Z,XI,YI,cubic );surf(XI,YI,ZI);30、进行插值后计算体积：底面积乘以平均高度xi=li nspace(1,4,61);yi=li nspace(1,5,81);XI,YI=meshgrid(xi,yi);Zl=i nterp2(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

9、.8222V =1.6658e+00830.上面是插值的方法解题，下面用拟合的方法解题。 为此编写了几个M函数：leftmatrix.mfun cti onU=leftmatrix(x,p,y,q)% U*a=V a为系数列矩阵，长度为p*q% U为左边p*q乘p*q矩阵% x,y为长度一致的列矩阵，给定点的坐标% p,q为拟合的函数中x,y的幕的最高次数m=le ngth(x);if (nargin=4) & (m=le ngth(y)error(enderror check check!);U_le ngth=p*q;U=zeros(U_le ngth,U_le ngth);for i=1

10、 : p*qfor j= 1 : p*qx_z=quotie nt(j-1,q)+quotie nt(i-1,q); y_z=mod(j_1,q)+mod(i_1,q);U(i,j)=qiuhe(x,x_z,y,y_z);endend% U为p*q阶方阵%赋值0，目的是分配内存% x 的幕的次数，quotie nt%y的幕的次数为求商rightmatrix.mfun cti onV=rightmatrix(x,p,y,q,z)% U*a=V% V为一个列向量长为p*q% x y z 为点的坐标%p q 分别为xy幕的最高次数if n arg in=5error( error check che

11、ck! rightmatrixendV=zeros(p*q,1);for i=1 : p*qx_z=quotie nt(i-1,q);y_z=mod(i-1,q);V(i,1)=qiuhe(x,x_z,y,y_z,z);endquuotie nt.mfun cti on sh=quotie nt(x,y)% sh 为x/y 的商sh=(x-mod(x,y)/y;qiuhe.mfun cti on he=qiuhe(x,p,y,q,z)% he xAp*yAq从 1 m 的和% x,y 向量长度相同% p,q分别为x,y的幕的次数%输入量至少为四，x,y行向量长度必需一样默认为元素全部为1的向量

12、m=le ngth(x);if (narginm 求和end下面一段程序先进行拟合，然后验证拟合的效果，具体操作：先输入x=._、= 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.7225.808.4725.2722.32;15.4721.3314.4924.8326.19;23.2826.4829.1412.0414.58;19.9523.7315.3518.011

13、6.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=i nv(U)*V;a_n=UV;for i=1 : length(a_n)%把长为p*q 的列向量a_n转换成p*q的矩阵aaii=quotie nt(i-1,q)+1; % quotie nt求商jj=mod(i-1,q)+1;aa(ii,jj)=a_n(i,1);endaa m=31;n=41; %m=4;n=5; XI,YI=mes

14、hgrid(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);是 xx,yy 代入所 拟 合的函数求出的函数 值函数为艺 aa(i,j)*xAi*yAj, (i=1.p,j=1.q)为 pxq 的系数的矩 阵%zz=0; % zz for i=1 : p %for j=1 : q % aaxt=xx.A(i-1);yt=yy.A(j-1);xy=xt.*yt; zz=zz+aa(i,j).*xy;endendZI=reshape(zz,n,m); surf(XI,YI,ZI); %axis(1 4 1 5 0 30) aa =1.0e