曲面擬合原理與實例

1、問題:給定一組坐標（xg,yg,zg） , g1,2,n，表示有n個點。要求用以下二元f(x, y)p,qi 1 j ajX yij 1,1p q1aiji 1 j 1i 1 j 1x y即f (x, y)a11ai2y2a13ya21xa22xy2a23xy多項式函數(shù)對所給的坐標進行擬合:q 1aiqyq 1a2qxy0i 1i 1i 12i 1 q 1aiixa2x yai3x y aiqx yqp 1p 12ap2Xy ap3X yIIIp 1 q 1apqX y設11XyX2x ,y2 y1ap1xp 1pqxya11ai2illa1qa21a22a2q,A11111 f+n*ap1a

2、p2IIIapq則函數(shù)又可表示為f （x, y）xTAy ,擬合的目標就是求出系數(shù)矩陣 A最小二乘法:構(gòu)造關(guān)于系數(shù)aj的多元函數(shù):p q/i 1 j 1g(ajX yi 1 j 1nzg)2s( ai1 | , a pq )g f(xg, yg)zg g 1點（an，apq）是多兀函數(shù)s（an,川,apq）的極小點，其中g(shù)為權(quán)函數(shù)，默認為1,所以點（an，apq ）必須滿足方程組saij在g 1的情況下，有aij2f (Xg,yg)Zgaj g 1n2f(Xg,yg)g 1n2f (Xg,yg)g 1n2xg1yg 1f(Xg,yg)兒匕比g 1因此可得叩石f(Xg，yg)i 1 j 1ZgX

3、g ygnxg1yg 1f(Xg, yg)g 1ni 1Xgg 1yg 1zgp qg1 a1 1ni 1 j 1X ygg 1p,qp,qa1,1u (i, j)1Xg ygni 1 j 1X yg Zgg 11 1a Xgyg1,1i 1 j 1X yg Zg1nz 11 i 1 j 1、(Xgyg Xg yg )g 1ni 1 j 1Xg ygg 111 i 1 j 1、yg xg yg)，ni 1 jXg yg Zgg 1p,qa u1,1(i,j)v(i, j)(i,j)上式實際共有p q個等式，可將這p q個等式寫成矩陣的形式有:Uii(1,1)II Upq(1,1)牛v(1,1)

4、4v(p,q)Un(p,q)III Upq(p,q) apq也就是U*a=V的形式，其中Un")| Upq(1,1)Uu”(P,q)|丨 Upq(p,q)311a -apqv(1,1)，v :v(p,q)U為pq pq階矩陣，實現(xiàn)函數(shù)為function A=leftmatrix（x,p,y,q） ； V為長pq的列 向量，實現(xiàn)函數(shù)為function B=rightmatrix（x,p,y,q,z）。這樣就可以算出列矩陣 a, 然后轉(zhuǎn)化成A。例子：某地區(qū)有一煤礦，為估計其儲量以便于開采，先在該地區(qū)進行勘探。假設該地區(qū)是一長方形區(qū)域，長為 4公里，寬為5公里。經(jīng)勘探得到如下數(shù)據(jù): 煤礦勘

5、探數(shù)據(jù)表編號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請你估計出此地區(qū)內(nèi)（2x4,1 y 5 ）煤的儲量，單位用立方米表示,并用電腦畫出該煤礦的三維圖象。如果直接畫出三維曲面圖形:clear;x=1:4;y=1:5

6、;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.29'surf(X ,Y,Z);X =12341234123412341234Y =11112222333344445555Z =13.7200 15.4700 23.2800 19.950025.800021.330026.480023.73008.470014.490029.140015.350025.270024.8300

7、12.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

8、);xi=li nspace(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

9、=2000*4000;V=S*H n =3321H =20.8222V =1.6658e+00830.上面是插值的方法解題，下面用擬合的方法解題。 為此編寫了幾個M函數(shù)：leftmatrix.mfun cti onU=leftmatrix(x,p,y,q)% U*a=V a為系數(shù)列矩陣，長度為p*q% U為左邊p*q乘p*q矩陣% x,y為長度一致的列矩陣，給定點的坐標% p,q為擬合的函數(shù)中x,y的幕的最高次數(shù)m=le ngth(x);if (nargin=4) & (m=le ngth(y)error(end'error check check!');U_le ng

10、th=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)=qiuhe(x,x_z,y,y_z);endend% U為p*q階方陣%賦值0，目的是分配內(nèi)存% x 的幕的次數(shù)，quotie nt%y的幕的次數(shù)為求商rightmatrix.mfun cti onV=rightmatrix(x,p,y,q,z)% U*a=V% V為一個列向量 長為p*q% x y z 為點的坐標%p q 分

11、別為xy幕的最高次數(shù)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);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 的和%

12、x,y 向量長度相同% p,q分別為x,y的幕的次數(shù)%輸入量至少為四，x,y行向量長度必需一樣)；%沒有z ,默認為元素全部為1的向量% 1->m 求和m=le ngth(x);if (nargin<4 )&(m=length(y)error('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=._、= z=. p=q=(注意x,y,z是向量); 擬合得到系

13、數(shù)a,也就是得到了擬合的函數(shù)； 根據(jù)擬合函數(shù)計算給定點(xx, yy)的函數(shù)值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);%

14、 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轉(zhuǎ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=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 代入所 擬 合的函數(shù) 求出的函數(shù) 值函數(shù)為藝 aa(i,j)*xAi*yAj, (i=1.p,j=1.q)為 pxq 的系數(shù)的矩 陣%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 3