當我們手中握有大量的資料時，對於二維的資料，我們會對他們進行直線擬合、對數擬合，圓曲線的擬合等等。這些擬合的方法都是運用的了非常古老而又非常有效的方法，即最小二乘法。 
今天給大家介紹一種三維球面資料的擬合方法，該方法也是運用的最小二乘的方法。旨在使擬合的半徑在均方意義下誤差達到最小。

公式推導

設擬合後的球面的球心為(x_0,y_0,z_0)及半徑r。 
對於每一點擬合後估計的值與實際值的差值為: 
 
則誤差的平方和為: 
 
注意E是x_0,y_0,z_0,r的函式。因此令E分別對x_0,y_0,z_0,r的偏導數等於0，即可求出x_0,y_0,z_0,r，有: 
 
令： 
 
則有： 
 
由(1)-(4)得 
 
由(2)-(4)得 
 
由(3)-(4)得 

Matlab模擬

首先生成若干個球面資料，然後加入一定能量的噪聲。最後利用上述的公式計算擬合後的球心座標和球面半徑，下面給出的是matlab模擬程式碼:

%最小二乘的方法進行擬合
clear all;
close all
clc;
R = 2;         %球面半徑
x0 = 100;      %球心x座標
y0 = 1;        %球心y座標
z0 = 76;       %球心z座標
%********************************生成隨機球面資料********************
 
****************
alfa = 0:pi/50:pi;
sita = 0:pi/50:2*pi;
num_alfa = length(alfa);
num_sita = length(sita);
x = zeros(num_alfa,num_sita);
y = zeros(num_alfa,num_sita);
z = zeros(num_alfa,num_sita);
for i = 1:num_alfa
    for j = 1:num_sita
        x(i,j) = x0+R*sin(alfa(i))*cos(sita(j));
        y(i,j) = y0+R*sin
 
(alfa(i))*sin(sita(j));
        z(i,j) = z0+R*cos(alfa(i));
    end
end

x = reshape(x,num_alfa*num_sita,1);
y = reshape(y,num_alfa*num_sita,1);
z = reshape(z,num_alfa*num_sita,1);
figure;
plot3(x,y,z,'*');
title('生成的沒有噪聲的球面資料');
%加入均值為0的高斯分佈噪聲 
amp = 0.1;
x = x + amp*rand(num_alfa*num_sita,1);
y = y + amp*rand(num_alfa*num_sita,1);
z = z + amp*rand(num_alfa*num_sita,1);
figure;
plot3(x,y,z,'*');
title('加入噪聲的球面資料');
%*******************************************************************************************
%球面擬合算法
num_points = length(x);
x_avr = sum(x)/num_points;
y_avr = sum(y)/num_points;
z_avr = sum(z)/num_points;

xx_avr = sum(x.*x)/num_points;
yy_avr = sum(y.*y)/num_points;
zz_avr = sum(z.*z)/num_points;
xy_avr = sum(x.*y)/num_points;
xz_avr = sum(x.*z)/num_points;
yz_avr = sum(y.*z)/num_points;

xxx_avr = sum(x.*x.*x)/num_points;
xxy_avr = sum(x.*x.*y)/num_points;
xxz_avr = sum(x.*x.*z)/num_points;
xyy_avr = sum(x.*y.*y)/num_points;
xzz_avr = sum(x.*z.*z)/num_points;
yyy_avr = sum(y.*y.*y)/num_points;
yyz_avr = sum(y.*y.*z)/num_points;
yzz_avr = sum(y.*z.*z)/num_points;
zzz_avr = sum(z.*z.*z)/num_points;
%計算求解線性方程的係數矩陣
A = [xx_avr - x_avr*x_avr,xy_avr - x_avr*y_avr,xz_avr - x_avr*z_avr;
     xy_avr - x_avr*y_avr,yy_avr - y_avr*y_avr,yz_avr - y_avr*z_avr;
     xz_avr - x_avr*z_avr,yz_avr - y_avr*z_avr,zz_avr - z_avr*z_avr];
b = [xxx_avr - x_avr*xx_avr + xyy_avr - x_avr*yy_avr + xzz_avr - x_avr*zz_avr;
     xxy_avr - y_avr*xx_avr + yyy_avr - y_avr*yy_avr + yzz_avr - y_avr*zz_avr;
     xxz_avr - z_avr*xx_avr + yyz_avr - z_avr*yy_avr + zzz_avr - z_avr*zz_avr];
b = b/2;

resoult = inv(A)*b;

x00 = resoult(1);     %擬合出的x座標
y00 = resoult(2);     %擬合出的y座標
z00 = resoult(3);     %擬合出的z座標
r = sqrt(xx_avr-2*x00*x_avr+x00*x00 + yy_avr-2*y00*y_avr+y00*y00 + zz_avr-2*z00*z_avr+z00*z00);   %擬合出的球半徑r

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16
• 17
• 18
• 19
• 20
• 21
• 22
• 23
• 24
• 25
• 26
• 27
• 28
• 29
• 30
• 31
• 32
• 33
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42
• 43
• 44
• 45
• 46
• 47
• 48
• 49
• 50
• 51
• 52
• 53
• 54
• 55
• 56
• 57
• 58
• 59
• 60
• 61
• 62
• 63
• 64
• 65
• 66
• 67
• 68
• 69
• 70
• 71
• 72
• 73
• 74
• 75
• 76
• 77

執行的結果如下： 
 
 
擬合後的x00=100.0502, y00=1.0491, z00=76.0512, r=2.0009與真實的x0=100, y0=1, z0=76, R=2非常的接近。

當然如果得到的球面資料不是在整個球面均勻分佈的也可以得到很不錯的擬合結果，當得到的資料如下圖所示： 
 
 

擬合後的x00=100.0510, y00=1.0537, z00=76.0540, r= 1.9952與真實值依然很接近。

原文地址

http://blog.csdn.net/hj199404182515/article/details/53462512

但是個人感覺他的求導推導有問題

歡迎大家一起來討論一下：

紅色標記部分

歡迎高手給以解釋，非常感謝

掃碼關注我們：跟著數理化走天下

獲得更多的資訊哦，一起交流，一起成長哦：微訊號：跟著數理化走天下，純屬個人的交流，無盈利目的