球面数据拟合算法+椭球数据拟合算法

当我们手中握有大量的数据时,对于二维的数据,我们会对他们进行直线拟合、对数拟合,圆曲线的拟合等等。这些拟合的方法都是运用的了非常古老而又非常有效的方法,即最小二乘法。
今天给大家介绍一种三维球面数据的拟合方法,该方法也是运用的最小二乘的方法。旨在使拟合的半径在均方意义下误差达到最小。

公式推导

设拟合后的球面的球心为(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)得
在这里插入图片描述
写成矩阵的形式:
在这里插入图片描述
求解该矩阵得到x_0,y_0,z_0值,然后带入(4)式中得到r的值。

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

运行的结果如下:
在这里插入图片描述
在这里插入图片描述
拟合后的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与真实值依然很接近。

空间二次曲面数据拟合算法

将上面球形情况进行更一般的推广,即取消了球面数据这一限制,数据可以是椭球面形式的,也就是说是任意的空间二次曲面形式的,可见球面的数据是它的一个特例。下面给出的是空间二次曲面的标准表达式:
在这里插入图片描述
一共有6个未知的参数,x0, y0, z0, A, B, C。写成一般式如下所示:
在这里插入图片描述
其中
在这里插入图片描述
对于有N个三维椭球面样本对其进行椭球面拟合,我们只需要对参数a,b,c,d,e,f进行估计,从而就可以得到x0, y0, z0, A, B, C。那么怎么利用样本去估计这些参数呢?这实际上就是模型参数估计的内容,模型参数估计有很多种方法,其中最基本的方法就是最小二乘法(Least Squares Method),由于它的原理直观,算法简单,收敛性能好,且不要求先验的统计知识,因而被广泛应用。最小二乘法是在1795年由大数学家高斯(C.F.Gauss)研究天体运动轨道问题提出的,它的基本原理是实际观测值与模型计算值的误差的平方和最小原理,由此而得名“最小二乘”法。应该注意的是最小二乘法是一种思想,由它衍生出来的公式可以有很多种。本次的空间二次数据拟合算法的推导利用的也是最小二乘法,不同的是本次是先估计参数a,b,c,d,e,f,然后间接的得到参数x0, y0, z0, A, B, C,这样做给公式推导带来了很大的方便,而且结果与直接推导的完全一样。

基于最小二乘的拟合方法公式推导

下面来利用最小二乘的思想进行推导,首先对于每个样本点(x_i,y_i,z_i),其误差可表示为:
在这里插入图片描述
则误差平方和为:
在这里插入图片描述
我们令E对于a,b,c,d,e,f的各一阶偏导数等于0,则有:
在这里插入图片描述
令:
在这里插入图片描述
带入化简得:
在这里插入图片描述
写成矩阵的形式有:
在这里插入图片描述
求解该六元一次方程组即可得到a,b,c,d,e,f的估计值。进而得到x0, y0, z0, A, B, C。

Matlab代码分析与仿真

代码主要有两个部分,第一部分是生成带噪声的椭球面数据,第二部分是利用上述推导的算法来估计相关的参数,下面给出相关的代码

clear all;
close all
clc;
A = 300;     % x方向上的轴半径
B = 400;     % y方向上的轴半径
C = 500;     % z方向上的轴半径
x0 = -120;   %椭球球心x坐标
y0 = -67;    %椭球球心y坐标
z0 = 406;    %椭球球心z坐标
SNR = 30;    %信噪比

%********************************生成随机椭球球面数据************************************
%利用在球坐标系下的参数方程来生成各个椭球的样本点
num_alfa = 100;
num_sita = 50;
alfa = (0:num_alfa-1)*1*pi/num_alfa;
sita = (0:num_sita-1)*2*pi/num_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 + A*sin(alfa(i))*cos(sita(j));
        y(i,j) = y0 + B*sin(alfa(i))*sin(sita(j));
        z(i,j) = z0 + C*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 = 10^(-SNR/20)*A;
x = x + amp*rand(num_alfa*num_sita,1);
y = y + amp*B/A*rand(num_alfa*num_sita,1);
z = z + amp*C/A*rand(num_alfa*num_sita,1);
figure;
plot3(x,y,z,'*');
string = ['加入噪声的椭球面数据,SNR=',num2str(SNR),'dB'];
title(string);
%*******************************************************************************************
%空间二次曲面拟合算法
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;
%四次项统计平均
yyyy_avr = sum(y.*y.*y.*y)/num_points;
zzzz_avr = sum(z.*z.*z.*z)/num_points;
xxyy_avr = sum(x.*x.*y.*y)/num_points;
xxzz_avr = sum(x.*x.*z.*z)/num_points;
yyzz_avr = sum(y.*y.*z.*z)/num_points;


%计算求解线性方程的系数矩阵
A0 = [yyyy_avr yyzz_avr xyy_avr yyy_avr yyz_avr yy_avr;
     yyzz_avr zzzz_avr xzz_avr yzz_avr zzz_avr zz_avr;
     xyy_avr  xzz_avr  xx_avr  xy_avr  xz_avr  x_avr;
     yyy_avr  yzz_avr  xy_avr  yy_avr  yz_avr  y_avr;
     yyz_avr  zzz_avr  xz_avr  yz_avr  zz_avr  z_avr;
     yy_avr   zz_avr   x_avr   y_avr   z_avr   1;];
%计算非齐次项
b = [-xxyy_avr;
     -xxzz_avr;
     -xxx_avr;
     -xxy_avr;
     -xxz_avr;
     -xx_avr];

resoult = inv(A0)*b;
%resoult = solution_equations_n_yuan(A0,b);

x00 = -resoult(3)/2;                  %拟合出的x坐标
y00 = -resoult(4)/(2*resoult(1));     %拟合出的y坐标
z00 = -resoult(5)/(2*resoult(2));     %拟合出的z坐标

AA = sqrt(x00*x00 + resoult(1)*y00*y00 + resoult(2)*z00*z00 - resoult(6));   % 拟合出的x方向上的轴半径
BB = AA/sqrt(resoult(1));                                                    % 拟合出的y方向上的轴半径
CC = AA/sqrt(resoult(2));                                                    % 拟合出的z方向上的轴半径

fprintf('拟合结果\n');
fprintf('a = %f, b = %f, c = %f, d = %f, e = %f, f = %f\n',resoult);
fprintf('x0 = %f, 相对误差%f\n',x00,abs((x00-x0)/x0));
fprintf('y0 = %f, 相对误差%f\n',y00,abs((y00-y0)/y0));
fprintf('z0 = %f, 相对误差%f\n',z00,abs((z00-z0)/z0));
fprintf('A = %f,  相对误差%f\n',AA,abs((A-AA)/A));
fprintf('B = %f,  相对误差%f\n',BB,abs((B-BB)/B));
fprintf('C = %f,  相对误差%f\n',CC,abs((C-CC)/C));

运行该程序得到如下结果:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
可以看到拟合的误差很小,同时我们也可以看到拟合后的参数A,B,C的相对误差比x0,y0,z0要小。

转自:https://blog.csdn.net/HJ199404182515/article/details/59480954
https://blog.csdn.net/HJ199404182515/article/details/53462512

版权声明:本文为Do_Not_Ask_Me原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/Do_Not_Ask_Me/article/details/103372105