曲面擬合原理與實例(共10頁)

1、精選優(yōu)質文檔-傾情為你奉上問題：給定一組坐標，表示有n個點。要求用以下二元多項式函數對所給的坐標進行擬合：即設則函數又可表示為，擬合的目標就是求出系數矩陣A。最小二乘法：構造關于系數的多元函數：點（，）是多元函數的極小點，其中為權函數，默認為1，所以點（，）必須滿足方程組 在的情況下，有因此可得令，則 上式實際共有個等式，可將這個等式寫成矩陣的形式有：也就是U*a=V的形式，其中，U為階矩陣，實現函數為function A=leftmatrix(x,p,y,q)；V為長的列向量，實現函數為function B=rightmatrix(x,p,y,q,z)。這樣就可以算出列矩陣a,然后轉化成A。

2、例子：某地區(qū)有一煤礦，為估計其儲量以便于開采，先在該地區(qū)進行勘探。假設該地區(qū)是一長方形區(qū)域，長為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請你估計出此地區(qū)內

3、（）煤的儲量，單位用立方米表示，并用電腦畫出該煤礦的三維圖象。如果直接畫出三維曲面圖形：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 4Y = 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5

4、 5 5 5Z = 13.7200 15.4700 23.2800 19.9500 25.8000 21.3300 26.4800 23.7300 8.4700 14.4900 29.1400 15.3500 25.2700 24.8300 12.0400 18.0100 22.3200 26.1900 14.5800 16.2900 粗略計算體積：底面積乘以平均高度。p=sum(Z);q=p(:,2,3,4);h=sum(q')/15v=2000*4000*h h = 20.0773v = 1.6062e+008 進行線性插值：xi=linspace(1,4,31);yi=linsp

5、ace(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

6、);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*H n = 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 為長度一致的列矩陣，給

7、定點的坐標% p,q 為擬合的函數中x,y的冪的最高次數 m=length(x);if (nargin=4) & (m=length(y) error('error check check!');end U_length=p*q; % U 為p*q階方陣U=zeros(U_length,U_length); % 賦值0，目的是分配內存for i=1 : p*q for j= 1 : p*q x_z=quotient(j-1,q)+quotient(i-1,q); % x 的冪的次數，quotient為求商 y_z=mod(j-1,q)+mod(i-1,q); % y 的

8、冪的次數 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=5 error('error check check! rightmatrix')end V=zeros(p*q,1);for i=1 : p*q x_z=quotient(i-1,q); y_z=mod(i-1,q); V(i,1)=qiuhe(x,x_z,y,y_z,z);end

9、quuotient.mfunction sh=quotient(x,y)% sh 為 x/y 的商 sh=(x-mod(x,y)/y;qiuhe.mfunction he=qiuhe(x,p,y,q,z)% he xp*yq 從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

10、,1);end he=0;for i=1:m he=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.

11、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=UV; 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;

12、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)*xi*yj，(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; endend ZI=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