曲面拟合原理与实例

问题：给定一组坐标(,,)g g g x y z ，1,2,g =…,n ，表示有n 个点。

要求用以下二元多项式函数对所给的坐标进行拟合：,11111,111(,)p qp qi j i j ijij ij i j f x y a xya x y ----=====∑∑∑即211112131212122232111211123111211123(,)q q q q i i i i q i i i iq p p p p q p p p pq a a y a y a y f x y a x a xy a xy a xy a x a x ya x y a x y a x a x y a x y a x y ------------++++=+++++++++++++++L L M L M L设1112121222221211,,q q p p pq p q a a a x y a a a x y a a a x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦x y A L LM M O M M M L则函数又可表示为(,)T f =x y x Ay ，拟合的目标就是求出系数矩阵A 。

最小二乘法：构造关于系数ij a 的多元函数：2112111111(,,)[(,)]()pqnni j pq g g g g g ij g g g i j s a a f x y z a x y z ωω--=====-=-∑∑∑∑L点（11a ，…，pq a ）是多元函数11(,,)pq s a a L 的极小点，其中g ω为权函数，默认为1，所以点（11a ，…，pq a ）必须满足方程组0ijsa ∂=∂ 在1g ω=的情况下，有{}21111111111[(,)]2[(,)][(,)]2[(,)]2(,)nggg g ij ijng g g g g g ij ni j g g g g g g n i j i j g gg g g g g g s f x yz a a f x y z f x y a f x y z x y x y f x y x y z ==--=----=∂∂=-∂∂⎧⎫∂⎪⎪=-⎨⎬∂⎪⎪⎩⎭=-⎡⎤=-⎣⎦∑∑∑∑因此可得111111(,)nni j i j ggg g g g g g g xyf x y x y z ----===∑∑1111111111p qnni j i j ggggg g g g g xya x y x y z αβαβαβ------=====∑∑∑∑ ,11111111,11p qnni j i j ggg gg g g g g xya x y x y z αβαβαβ------====∑∑∑ ,1111111,111()p qn n i j i j g g g g g g g g g a x y x y x y z αβαβαβ------===⎡⎤=⎢⎥⎣⎦∑∑∑ 令11111(,)()ni j g g g gg u i j x y x yαβαβ----==∑，111(,)ni j gg g g v i j x y z --==∑ 则,1,1(,)(,)p qa u i j v i j αβαβαβ==∑(,)(11),(,)i j p q =，…, 上式实际共有p q ⨯个等式，可将这p q ⨯个等式写成矩阵的形式有：111111(1,1)(1,1)(1,1)(,)(,)(,)pq pq pq u u a v u p q u p q a v p q ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦L M O M M M L 也就是U*a=V 的形式，其中1111(1,1)(1,1)(,)(,)pq pq u u u p q u p q ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦U L M OM L ，11pq a a ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦a M ，(1,1)(,)v v p q ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦V MU 为pq pq ⨯阶矩阵，实现函数为function A=leftmatrix(x,p,y,q)；V 为长pq 的列向量，实现函数为function B=rightmatrix(x,p,y,q,z)。

这样就可以算出列矩阵a ,然后转化成A 。

例子：某地区有一煤矿，为估计其储量以便于开采，先在该地区进行勘探。

假设该地区是一长方形区域，长为4公里，宽为5公里。

经勘探得到如下数据：请你估计出此地区内（51,42≤≤≤≤y x ）煤的储量，单位用立方米表示，并用电脑画出该煤矿的三维图象。

如果直接画出三维曲面图形：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 4 Y =1 1 1 12 2 2 23 3 3 34 4 4 45 5 5 5 Z =13.7200 15.4700 23.2800 19.950025.8000 21.3300 26.4800 23.73008.4700 14.4900 29.1400 15.350025.2700 24.8300 12.0400 18.010022.3200 26.1900 14.5800 16.2900粗略计算体积：底面积乘以平均高度。

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;for j=21:61for i=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.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!');endU_length=p*q; % U 为p*q阶方阵U=zeros(U_length,U_length); % 赋值0，目的是分配内存for i=1 : p*qfor j= 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.mfunction V=rightmatrix(x,p,y,q,z)% U*a=V% V 为一个列向量长为p*q% x y z 为点的坐标%p q 分别为x y幂的最高次数if nargin~=5error('error check check! rightmatrix')endV=zeros(p*q,1);for i=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.mfunction sh=quotient(x,y)% sh 为 x/y 的商sh=(x-mod(x,y))/y;qiuhe.mfunction he=qiuhe(x,p,y,q,z)% he x^p*y^q 从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);endhe=0;for i=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.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=U\V;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)*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.2678 0.2132 -0.0624 0.0058-0.7287 1.3972 -0.9275 0.2412 -0.02100.4416 -0.8415 0.5487 -0.1407 0.0122-0.0680 0.1295 -0.0839 0.0214 -0.0018注意：权函数在拟合的函数非常重要，不过她只能按你遇到的具体问题来取，我这里为1；当p，q越大时，拟合的函数与原数据的方差越小，但是有可能函数本身抖动非常厉害，可以画图看出来。