椭球拟合简介
MESE IMＵ中，x,y,z轴的度量单位并不相同，假设各轴之间相互直。
加计静止状态（也就是只受重力的状态下），各个姿态只受重力的，x,y,z轴值（假设x,y,z轴相互垂直并且度量单位都一致，如三轴的度量单位都是2000（1g），量程固定的情况下），在三维空间中，重力点都在一个球面上，但各轴之间的度量单位都会有偏差，所以各姿态重力点都落在一个椭球面上，椭球的中心，就是加速度的偏移量，也就是校准值。
图1为标定原理公式
第一个椭球公式A、B、C，X0、Y0、Z0对应：Scale参数（比例因子）、Bias参数。
图1
令E对于a,b,c,d,e,f的各一阶偏导数等于0，通过求解偏导数方程组，得到a,b,c,d,e,f参数。进而得到A、B、C（Scale参数，比例因子），X0、Y0、Z0（Bias参数）。A、B、C，X0、Y0、Z0对应第一个椭球公式。
详细的展开数学原理以及偏导公式请参考：空间二次曲面数据拟合算法推导及仿真分析。
实测分析
采集19组加计静止数据。椭球拟合-标定加计（bias、scale）结果如图2所示：
可见加计静止数据恰好拟合成椭球面。
拟合数据分析：
1、采用加计输出模值应为重力加速度为评价标准。
rawData = 1.001998g
fitData = 1.00000g
2、补偿误差为 0.001998 * 9.8035 = 0.019587m/s2
19组加计静止数据如下：
ax ay az
-17.538 25.758 2006.695
-1997.923 -27.4 -25.143
-18.408 26.256 2006.927
-555.95 16.058 1925.793
-46.025 384.178 1969.542
37.619 -507.349 -1938.779
11.991 450.054 -1956.484
583.968 59.653 -1916.821
-313.098 -22.123 -1982.363
-16.711 26.097 2006.502
-13.122 -409.013 1963.334
581.886 -28.452 1921.501
-12.521 433.975 1959.581
30.725 -1985.632 211.399
-30.957 1991.537 198.426
-1998.108 -25.868 -3.205
-1970.602 -311.407 99.16
-1973.875 312.339 36.82
802.364 -92.147 1837.896
椭球拟合matlab代码如下：
clear all;
close all
clc;mat_path = './Fit.csv';fprintf('opening the mat file.\n')
data_imu = csvread(mat_path);imu_ax = data_imu(:,1);
imu_ay = data_imu(:,2);
imu_az = data_imu(:,3);
##deta_T = data_imu(:,9);x = imu_ax(1:1:end);
y = imu_ay(1:1:end);
z = imu_az(1:1:end);%转换成一维的数组便于后续的处理
##x = [87 301 274 312 -3805 4389 261 327 -1963 3024];
##y = [-52 -45 4088 -4109 -24 6 2106 -2047 -13 18];
##z = [-4454 3859 -303 -305 -390 -228 -3848 -3880 -3797 -3449];figure;
plot3(x,y,z,'r*');
hold on;
##string = ['加入噪声的椭球面数据，SNR=',num2str(SNR),'dB'];
title('Ellipsoid Specific Fitting');%*******************************************************************************************
%空间二次曲面拟合算法
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\n',x00);
fprintf('y0 = %f\n',y00);
fprintf('z0 = %f\n',z00);
fprintf('A = %f\n',AA);
fprintf('B = %f\n',BB);
fprintf('C = %f\n',CC);[x1, y1, z1] = ellipsoid(x00,y00,z00,AA,BB,CC);
##figure
mesh(x1, y1, z1)
legend('rawData', 'fitElli')
axis equalrawNoise = sqrt(x.^2 + y.^2 + z.^2) / 2000;
rawNoise = sum(rawNoise) / length(rawNoise);
fprintf('rawNoise = %f\n',rawNoise);
fitNoise = sqrt(((x-x00)/AA).^2 + ((y-y00)/BB).^2 + ((z-z00)/CC).^2);
fitNoise = sum(fitNoise) / length(fitNoise);
fprintf('fitNoise = %f\n',fitNoise);
