MATLAB 曲面擬合
阿新 • • 發佈:2020-09-10
這裡用到的還是最小二乘方法,和上一次這篇文章原理差不多。
就是首先構造最小二乘函式,然後對每一個係數計算偏導,構造矩陣乘法形式,最後解方程組。
比如有一個二次曲面:z=ax^2+by^2+cxy+dx+ey+f
首先構造最小二乘函式,然後計算係數偏導(我直接手寫了):
解方程組(下圖中A矩陣後面求和符號我就沒寫了啊),然後計算C:
程式碼如下:
1 clear all;
2 close all;
3 clc;
4
5 a=2;b=2;c=-3;d=1;e=2;f=30; %係數
6 n=1:0.2:20;
7 x=repmat(n,96,1 );
8 y=repmat(n',1,96);
9 z=a*x.^2+b*y.^2+c*x.*y+d*x+e*y +f; %原始模型
10 surf(x,y,z)
11
12 N=100;
13 ind=int8(rand(N,2)*95+1);
14
15 X=x(sub2ind(size(x),ind(:,1),ind(:,2)));
16 Y=y(sub2ind(size(y),ind(:,1),ind(:,2)));
17 Z=z(sub2ind(size(z),ind(:,1),ind(:,2)))+rand(N,1)*20; %生成待擬合點,加個噪聲
18
19 hold on;
20 plot3(X,Y,Z,'o');
21
22 A=[N sum(Y) sum(X) sum(X.*Y) sum(Y.^2) sum(X.^2);
23 sum(Y) sum(Y.^2) sum(X.*Y) sum(X.*Y.^2) sum(Y.^3) sum(X.^2.*Y);
24 sum(X) sum(X.*Y) sum(X.^2) sum(X.^2.*Y) sum(X.*Y.^2) sum(X.^3);
25 sum(X.*Y) sum(X.*Y.^2) sum(X.^2.*Y) sum(X.^2.*Y.^2) sum(X.*Y.^3 ) sum(X.^3.*Y);
26 sum(Y.^2) sum(Y.^3) sum(X.*Y.^2) sum(X.*Y.^3) sum(Y.^4) sum(X.^2.*Y.^2);
27 sum(X.^2) sum(X.^2.*Y) sum(X.^3) sum(X.^3.*Y) sum(X.^2.*Y.^2) sum(X.^4)];
28
29 B=[sum(Z) sum(Z.*Y) sum(Z.*X) sum(Z.*X.*Y) sum(Z.*Y.^2) sum(Z.*X.^2)]';
30
31 C=inv(A)*B;
32
33 z=C(6)*x.^2+C(5)*y.^2+C(4)*x.*y+C(3)*x+C(2)*y +C(1); %擬合結果
34
35 mesh(x,y,z)
結果如下,深色曲面是原模型,淺色曲面是用噪聲資料擬合的模型: