曲面拟合原理与实例

1、问题:给定一组坐标（Xg,yg,Zg），g=1,2,n ,表示有n个点。要求用以下二元多项式函数对所给的坐标进行拟合:p,qf (X, y)八 aijXij占1i J j A.i丄 j丄aj xyf (x, y) =a.a12ya21Xa22xya.3 y22a23 xyq 4a.qy,qSqXyi _±ai.xi _1ai2X yi _L 2ai3X yi A q 4aiqX yapixp-ap2"yap3X丄p_Lq Aa pq X y1q2qxpqf (x, y)=x T A y则函数又可表示为，拟合的目标就是求出系数矩阵 A。00最小二乘法：构造关于系数aj的多元函

2、数：nnp q2i i j i2s（aii, ,apq）二'gf （Xg,yg） Zgg（二二 aij X y - Zg）g =Lgi=t j=t点（a.,，a pq ）是多兀函数s（a.i,，apq）的极小点，其中,为权函数，默 认为i,所以点（aii，apq ）必须满足方程组:s0n2f (Xgg) -Zg.'aij ；'aij g ±n .：二' 2 f(Xg,yg)Zg f(Xg,yg)n八辽f (Xgyg) Zglx；丄g企ni丄 j丄i ± j丄-=2、|Xg yg f (Xg, yg)Xg ygZgg =1 因此可得nnl i

3、丄 j 丄 i 4 j 4'Xgygf(Xg,yg) J Xg g Zgg ±g -1np q-n、 X；丄yg 一 - a.Xg7ygJ=嘉 X/yZg g 11 |，：1g =1np ,qT- i丄j 4亍Ct丄 B丄Xg ygIa_：Xg jgg 土: : ±i 'n寸 i丄 j 4 =Xg YgZgg 土P,qn：nT- . 丁 /Of_LP_Li_Lj_Lp i-Lj_La母送(Xgyf Xg yg )=瓦 Xg ygZg母壬1 g£g壬nU (i, j)八(g=L、誥1Xg y|.L i 1 j 1 xg Xg Yg )ni -1 j

4、-1V(i, J) =Xg yg Zgg =Lp,qa：u：(i, j)二 v(i, j) J：；： 土i(i, j) =(1,1),，(p,q)上式实际共有p q个等式，可将这p q个等式写成矩阵的形式有:、1(1,1) u pq (1,1)an 1_v(1,1) |a+alam:I :=:®1(p,q) Upq (p,q) apq 一(p,q)也就是U*a=V的形式，其中"un (1,1)11 (p, q)Upq (1,1)U pq ( p , q)anaapq 一v(1,1)aV ='.v( p,q) jU为pq pq阶矩阵，实现函数为function A=l

5、eftmatrix（x,p,y,q） ； V为长pq的列 向量，实现函数为function B=rightmatrix（x,p,y,q,z）。这样就可以算出列矩阵 a, 然后转化成A。例子：某地区有一煤矿，为估计其储量以便于开采，先在该地区进行勘探。假设 该地区是一长方形区域，长为 4公里，宽为5公里。经勘探得到如下数据：煤矿勘探数据表编号121 31 41 5678910横坐标（公里）1111122222纵坐标（公里）1234512345煤层厚度（米）13.7225.808.47 25.27 22.3215.4721.3314.4924.8326.19编号111213141516171819

6、20横坐标（公里）3333344444纵坐标（公里）1234512345煤层厚度（米）:23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29请你估计出此地区内（2_x_4,1_y_5 ）煤的储量，单位用立方米表示,并用电脑画出该煤矿的三维图象。如果直接画出三维曲面图形：clear;x=1:4;y=1:5;X,Y=meshgrid(x,y)Z=13.7225.808.4725.2722.32;15.4721.3314.4924.8326.19;23.2826.4829.1412.0414.58;19.9523.7315.3518.0116.29&#

7、39;surf(X ,Y,Z);X =12341234123412341234Y =11112222333344445555Z =13.7200 15.4700 23.2800 19.950025.800021.330026.480023.73008.470014.490029.140015.350025.270024.830012.040018.010022.320026.190014.580016.290030 *粗略计算体积：底面积乘以平均高度。p=sum(Z);q=P(：,2,3,4);h=sum(q')/15v=2000*4000*hh =20.0773v =1.6062e+

8、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);进行三次多项式插值：41);XI,YI=meshgrid(xi,yi);xi=li nspace(1,4,31 );yi=li nspace(1,5,Zl=interp2(X,Y,Z,XI,YI,'cubic ');surf(XI,YI,ZI);进行插值后计算体积：底面积乘以平均高度xi=li nspace(1,4,6

9、1);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.8222V =1.6658e+008上面是插值的方法解题，下面用拟合的方法解题。 为此编写了几个M函数：leftmatrix.mfun cti onU=leftmatrix(x,p,y,q)% U*a=V a为系

10、数列矩阵，长度为p*q% U为左边p*q乘p*q矩阵% x,y为长度一致的列矩阵，给定点的坐标% p,q为拟合的函数中x,y的幕的最高次数m=le ngth(x);if (n argi n=4) & (m=le ngth(y)error( 'error check check!');endU_le ngth=p*q;U=zeros(U_le ngth,U_le ngth);for i=1 : 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)=

11、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 check! rightmatrix'endV=zeros(p*q,1);for i=1 : p*qx_z=quotie nt(i-1,q);y_z=mod(i-1,q)

12、;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的幕的次数m=le ngth(x);%输入量至少为四，x,y行向量长度必需一样)；%没有z ,默认为元素全部为1的向量% 1->m 求和if (nargin<4 )&(m=length(y)error(

13、9;error check check!'endif n arg in=4z=on es(m,1);endhe=0;for i=1:mhe=he+x(i)Ap * y(i)Aq*z(i);end下面一段程序先进行拟合，然后验证拟合的效果，具体操作：先输入 x=.、= 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.47

14、21.3314.4924.8326.19;23.2826.4829.1412.0414.58;19.9523.7315.3518.0116.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求

15、商jj=mod(i-1,q)+1;aa(ii,jj)=a_n(i,1);endaa m=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);是 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