曲面拟合原理与实例

问题：给定一组坐标(,,)gggxyz，1,2,g…,n，表示有n个点。要求用以下二元多项式函数对所给的坐标进行拟合：,11111,111(,)pqpqijijijijijijfxyaxyaxy即211112131212122232111211123111211123(,)qqqqiiiiqiiiiqppppqppppqaayayayfxyaxaxyaxyaxyaxaxyaxyaxyaxaxyaxyaxy设1112121222221211,,qqpppqpqaaaxyaaaxyaaaxyxyA则函数又可表示为(,)TfxyxAy，拟合的目标就是求出系数矩阵A。最小二乘法：构造关于系数ija的多元函数：2112111111(,,)[(,)]()pqnnijpqgggggijgggijsaafxyzaxyz点（11a，…，pqa）是多元函数11(,,)pqsaa的极小点，其中g为权函数，默认为1，所以点（11a，…，pqa）必须满足方程组0ijsa在1g的情况下，有21111111111[(,)]2[(,)][(,)]2[(,)]2(,)nggggijijnggggggijnijggggggnijijggggggggsfxyzaafxyzfxyafxyzxyxyfxyxyz因此可得111111(,)nnijijgggggggggxyfxyxyz1111111111pqnnijijgggggggggxyaxyxyz,11111111,11pqnnijijgggggggggxyaxyxyz,1111111,111()pqnnijijgggggggggaxyxyxyz令11111(,)()nijggggguijxyxy，111(,)nijggggvijxyz则,1,1(,)(,)pqauijvij(,)(11),(,)ijpq，…,上式实际共有pq个等式，可将这pq个等式写成矩阵的形式有：111111(1,1)(1,1)(1,1)(,)(,)(,)pqpqpquuavupqupqavpq也就是U*a=V的形式，其中1111(1,1)(1,1)(,)(,)pqpquuupqupqU，11pqaaa，(1,1)(,)vvpqVU为pqpq阶矩阵，实现函数为functionA=leftmatrix(x,p,y,q)；V为长pq的列向量，实现函数为functionB=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请你估计出此地区内（51,42yx）煤的储量，单位用立方米表示，并用电脑画出该煤矿的三维图象。如果直接画出三维曲面图形：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]'surf(X,Y,Z);X=12341234123412341234Y=11112222333344445555Z=13.720015.470023.280019.950025.800021.330026.480023.73008.470014.490029.140015.350025.270024.830012.040018.010022.320026.190014.580016.290012341234551015202530粗略计算体积：底面积乘以平均高度。p=sum(Z);q=p(:,[2,3,4]);h=sum(q')/15v=2000*4000*hh=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;forj=21:61fori=1:81H=H+ZI(i,j);n=n+1;endendnH=H/nS=2000*4000;V=S*Hn=3321H=20.8222V=1.6658e+008上面是插值的方法解题，下面用拟合的方法解题。为此编写了几个M函数：leftmatrix.mfunctionU=leftmatrix(x,p,y,q)%U*a=Va为系数列矩阵，长度为p*q%U为左边p*q乘p*q矩阵%x,y为长度一致的列矩阵，给定点的坐标%p,q为拟合的函数中x,y的幂的最高次数m=length(x);if(nargin~=4)&(m~=length(y))error('errorcheckcheck!');endU_length=p*q;%U为p*q阶方阵U=zeros(U_length,U_length);%赋值0，目的是分配内存fori=1:p*qforj=1:p*qx_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.mfunctionV=rightmatrix(x,p,y,q,z)%U*a=V%V为一个列向量长为p*q%xyz为点的坐标%pq分别为xy幂的最高次数ifnargin~=5error('errorcheckcheck!rightmatrix')endV=zeros(p*q,1);fori=1:p*qx_z=quotient(i-1,q);y_z=mod(i-1,q);V(i,1)=qiuhe(x,x_z,y,y_z,z);endquuotient.mfunctionsh=quotient(x,y)%sh为x/y的商sh=(x-mod(x,y))/y;qiuhe.mfunctionhe=qiuhe(x,p,y,q,z)%hex^p*y^q从1－m的和%x,y向量长度相同%p,q分别为x,y的幂的次数m=length(x);if(nargin4)&(m~=length(y))%输入量至少为四，x,y行向量长度必需一样error('errorcheckcheck!');endifnargin==4%没有z,默认为元素全部为1的向量z=ones(m,1);endhe=0;fori=1:mhe=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.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);%a_n=inv(U)*V;a_n=U\V;fori=1:length(a_n)%把长为p*q的列向量a_n转换成p*q的矩阵aaii=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代入所拟合的函数求出的函数值fori=1:p%函数为Σaa(i,j)*x^i*y^j，(i=1...p,j=1...q)forj=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([1415030])aa=1.0e+003*0.1465-0.26780.2132-0.06240.0058-0.72871.3972-0.92750.2412-0.02100.4416-0.84150.5487-0.14070.0122-0.06800.1295-0.08390.0214-0.0018注意：权函数在拟合的函数非常重要，不过她只能按你遇到的具体问题来取，我这里为1；当p，q越大时，拟合的函数与原数据的方差越小，但是有可能函数本身抖动非常厉害，可以画图看出来。