曲面拟合原理与实例

1、问题:给定一组坐标(x , y , z ), g =i,2,n 多项式函数对所给的坐标进行拟合：表示有n个点。要求用以下二元f (x, y) = axi-iyj-i =ijj=i,iaxi-i yj-iiji=1 j=1f (x, y) =aii+a xM21+ai2 j+a xy+a y 2+a xy 2+L+L+a y q-iiq+a xyq-i2q+a xi-iiiM+a xi-i ji 2+a xi-i y 2i3+a xi-iyq-iiq+a xp-iP1+a xp-i yp2+a xp-i y2P 3+a xp-i yq-ipqTiLraaaxyii12Liq,A =aaax 2,

2、y =y 221M22MO2 qMMMLaaa_ xP _ yq _L plp2pq则函数又可表示为f (x, y) = xTAy，拟合的目标就是求出系数矩阵A。最小二乘法：构造关于系数a的多元函数:s(a ,L , a )= E f (x , y ) - z 2 = (E a xi-i yj-i - z )pq2gg gggijg=ig=ii=ij=i点(a11a )是多元函数s(a ,Lii , a )的极小点其中为权函数，默认为i，所以点(a厂。叫)必须满足方程组ds=0daij在s = i的情况下，有 g6s _ Q Ei互=/乙f(1七)一Jij ij g =1=E2f(x ,y )

3、-z f(x ,y ) g g g oag gijg = 1=E zf(x , y ) - z xi-1 yj-J g g g g g g = 1=2E xi y j-1f (x , y ) - x- y j-1 z = 1因此可得xi-i yj-i f (x , y ) = xi-i yj-1 zg=1g=1Exi-1 yj-1 E Ea弗 xa-1 y P-1=Exi-1 yj-1 z g=1a=1 p =1g=1.EE.xi-1 yj-1 a xa-1 yP-1 =xi-1 yj-1 z g =1ap=1,1g=1aP=1,1令a E (xa-1 y P-1 xi-1 yj-1)g=1=

4、Exi-1 yj-1 zg=1u p (i, j) = E (xa-1 yP-1xi-1yj-1)，v(i, j) = En x yj-1 zg = 1g = 1归弗弗(i，j)= v(i，j)0j)=(),(p,q) aP=1,1上式实际共有px q个等式，可将这px q个等式写成矩阵的形式有:u11(1,1) LUpq (1)Ma M= v(1,1)-MMO(p, q)LU (p, q)a1- pqv( p, q) _也就是U*a=V的形式，其中u (1,1)Lu (1,1)a11 v(1,1)-U =11MOpqM，a =M，V=MU(p, q )Lupq ( p, q )aL pqJv

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

6、4.4924.8326.19编号11121314151617181920横坐标（公里）3333344444纵坐标（公里）1234512345煤层厚度（米）23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29请你估计出此地区内（2 x 4,1 y 以的和% x,y向量长度相同% p,q分别为x,y的幂的次数m=length(x);if (narginm 求和end下面一段程序先进行拟合，然后验证拟合的效果，具体操作：先输入 x=.y=.z=.p=.q=(注意 x,y,z 是向量)；拟合得到系数a，也就是得到了拟合的函数；根据拟合函数计算给定点(xx

7、, 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.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);%

8、a_n=inv(U)*V;%把长为p*q的列向量a_n转换成？*4的矩阵aa% quotient 求商a_n=UV;for i=1 : length(a_n) ii=quotient(i-1,q)+1;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=

9、zeros(m*n,1);%zz=0;% zz是xx,yy代入所拟合的函数求出的函数值for i=1 : p% 函数为aa(i,j)*x人i*y人j, (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;endendZI=reshape(zz,n,m);surf(XI,YI,ZI); %axis(1 4 1 5 0 30)aa =1.0e+003 火0.1465-0.26780.2132-0.06240.0058-0.72871.3972-0.92750.2412-0.02100.4416-0.84