曲面擬合原理與實(shí)例

1、問題:給定一組坐標(biāo)(x , y , z ), g =i,2,n 多項(xiàng)式函數(shù)對(duì)所給的坐標(biāo)進(jìn)行擬合：表示有n個(gè)點(diǎ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則函數(shù)又可表示為f (x, y) = xTAy，擬合的目標(biāo)就是求出系數(shù)矩陣A。最小二乘法：構(gòu)造關(guān)于系數(shù)a的多元函數(shù):s(a ,L , a )= E f (x , y ) - z 2 = (E a xi-i yj-i - z )pq2gg gggijg=ig=ii=ij=i點(diǎn)(a11a )是多元函數(shù)s(a ,Lii , a )的極小點(diǎn)其中為權(quán)函數(shù)，默認(rèn)為i，所以點(diǎn)(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上式實(shí)際共有px q個(gè)等式，可將這px q個(gè)等式寫成矩陣的形式有: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 階矩陣，實(shí)現(xiàn)函數(shù)為 function A=leftmatrix（x,p,y,q）； V 為長(zhǎng) pq 的列 向量，實(shí)現(xiàn)函數(shù)為function B=rightmatrix（x,p,y,q,z）。這樣就可以算出列矩陣a, 然后轉(zhuǎn)化成A。例子：某地區(qū)有一煤礦，為估計(jì)其儲(chǔ)量以便于開采，先在該地區(qū)進(jìn)行勘探。假設(shè)該地區(qū)是一長(zhǎng)方形區(qū)域，長(zhǎng)為4公里，寬為5公里。經(jīng)勘探得到如下數(shù)據(jù):煤礦勘探數(shù)據(jù)表編號(hào)12345678910橫坐標(biāo)（公里）1111122222縱坐標(biāo)（公里）1234512345煤層厚度（米）13.7225.808.4725.2722.3215.4721.331

6、4.4924.8326.19編號(hào)11121314151617181920橫坐標(biāo)（公里）3333344444縱坐標(biāo)（公里）1234512345煤層厚度（米）23.2826.4829.1412.0414.5819.9523.7315.3518.0116.29請(qǐng)你估計(jì)出此地區(qū)內(nèi)（2 x 4,1 y 以的和% x,y向量長(zhǎng)度相同% p,q分別為x,y的冪的次數(shù)m=length(x);if (narginm 求和end下面一段程序先進(jìn)行擬合，然后驗(yàn)證擬合的效果，具體操作：先輸入 x=.y=.z=.p=.q=(注意 x,y,z 是向量)；擬合得到系數(shù)a，也就是得到了擬合的函數(shù)；根據(jù)擬合函數(shù)計(jì)算給定點(diǎn)(xx

7、, yy)的函數(shù)值z(mì)z=f(xx, yy)并進(jìn)行畫圖檢驗(yàn)。程序保存于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;%把長(zhǎng)為p*q的列向量a_n轉(zhuǎ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代入所擬合的函數(shù)求出的函數(shù)值for i=1 : p% 函數(shù)為aa(i,j)*x人i*y人j, (i=1.p,j=1.q)for j=1 : q% aa pxq的系數(shù)的矩陣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