曲面擬合原理與實例

1、問題：給定一組坐標，表示有n個點。要求用以下二元多項式函數(shù)對所給的坐標進行擬合：即設則函數(shù)又可表示為，擬合的目標就是求出系數(shù)矩陣a。最小二乘法：構造關于系數(shù)的多元函數(shù)：點（，）是多元函數(shù)的極小點，其中為權函數(shù)，默認為1，所以點（，）必須滿足方程組 在的情況下，有因此可得令，則 上式實際共有個等式，可將這個等式寫成矩陣的形式有：也就是u*a=v的形式，其中，u為階矩陣，實現(xiàn)函數(shù)為function a=leftmatrix(x,p,y,q)；v為長的列向量，實現(xiàn)函數(shù)為function b=rightmatrix(x,p,y,q,z)。這樣就可以算出列矩陣a,然后轉化成a。例子：某地區(qū)有一煤礦，為估

2、計其儲量以便于開采，先在該地區(qū)進行勘探。假設該地區(qū)是一長方形區(qū)域，長為4公里，寬為5公里。經(jīng)勘探得到如下數(shù)據(jù)：煤礦勘探數(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ū)內（）煤的儲量，單位用立方米

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.29surf(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 5 5 5z = 13.7200

4、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=linspace(1,5,41);xi,yi=meshg

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

6、,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函數(shù)：leftmatrix.mfunction u=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=length(x);if (nargin=4

7、) & (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 的冪的次數(shù)，quotient為求商 y_z=mod(j-1,q)+mod(i-1,q); % y 的冪的次數(shù) u(i,j)=qiuhe(x,x_z,y,y_z); endendrightmatrix.mfunction v=right

8、matrix(x,p,y,q,z)% u*a=v % v 為一個列向量 長為p*q% x y z 為點的坐標 %p q 分別為x y冪的最高次數(shù) 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);endquuotient.mfunction sh=quotient(x,y)% sh 為 x/y 的商 sh=(x-mod(x,y)/y;qiuhe.mfun

9、ction he=qiuhe(x,p,y,q,z)% he xp*yq 從1m的和% x,y 向量 長度相同% p,q分別為x,y的冪的次數(shù)m=length(x);if (narginm 求和end下面一段程序先進行擬合，然后驗證擬合的效果，具體操作：先輸入x=y=z=p=q= (注意x,y,z是向量)；擬合得到系數(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.72 25.80 8.47 25.27 22.32; 15.47

10、 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=uv; for i=1 : length(a_n) % 把長為p*q 的列向量a_n轉換成p*q的矩陣aa ii=quotient(i-1,q)+1;

11、 % 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 代入所擬合的函數(shù) 求出的函數(shù)值for i=1 : p % 函數(shù)為aa(i,j)*xi*yj，(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; 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.3972