北航惯性导航综合实验二实验报告.doc

惯性测量单元安装误差系数标定实验

二零一三年六月十日
2.1 惯性测量单元安装误差系数标定试验
一、实验目的
1、掌握惯性测量单元（inertial measurement unit，IMU）的标度系数、安装误差、零偏的标定方法；
2、利用现有实验条件实现实验过程的设计。
二、实验内容
利用单轴速率转台，进行IMU的安装误差系数标定，并通过公式计算该安装误差系数。
三、实验系统组成
单轴速率位置转台、MEMS惯性测量单元、稳压电源、数据采集系统。
四、实验原理
IMU安装误差系数的计算方法
通常，惯导系统至少需要三个陀螺和三个加速度计，用以感知载体的三轴角速度和加速度变化。将这些陀螺和加计按照敏感轴两两正交的方式集成在一起，安装在一个结构框架上，便构成了一个能感知完整惯性测量信息的小型系统，称之为惯性测量单元。对惯性测量单元进行标定时，除了要对其中的陀螺、加速度计进行常规标定外，还要考虑由于安装时不能严格保证敏感轴两两正交所带来的交叉耦合误差，即，要对IMU的安装误差进行标定，测量出不正交角。因此，在考虑IMU的安装误差、标度因数误差、零偏误差的情况下，建立东北天坐标系下IMU的角速度通道误差方程。
(1)
式中为惯性系统i轴向陀螺输出角速度，为i轴向的输入角速度；为i轴向陀螺零偏；为i轴向陀螺标度因数；为角速度通道的安装误差系数，i和j为坐标轴X，Y，Z的统称。
设输入矩阵为，输出矩阵为，则标度因数、安装误差系数与陀螺漂移组成的矩阵可按最小二乘法估计为：

类似，可计算加速度计的标度因数、安装误差系数与加计零偏。
设输入矩阵为，输出矩阵为，则标度因数、安装误差系数与陀螺漂移组成的矩阵可按最小二乘法估计为：

五、实验内容
1、陀螺安装误差测试实验
1) 速率转台处于“停止”状态，接通电源，预热至IMU工作稳定；
2) 分别以10°/s，20°/s，40°/s，60°/s，80°/s的速率转动转台，打开监控计算机中的数据采集软件。在每一个旋转速率下，转台正转，旋转稳定后，采集转台旋转360°的过程中IMU的输出数据，停转，存储数据；转台反转，如上再次采集IMU输出数据，停转，存储数据；
3) 翻转工装，依次使得陀螺敏感轴X、Y轴依次平行于转台旋转轴，在每个位置上重复上述步骤，稳定后记录转动相应敏感轴的角速度当量均值并保存数据；
六、实验结果

陀螺数据求平均处理

10°/s
20°/s
40°/s
60°/s
80°/s

X轴反向

-9.9123

-19.9246

-39.8253

-59.6601

-79.1389

Y轴反向

-10.3317

-19.9998

-40.4727

-59.7430

-79.2054
Z轴反向
-9.9690
-19.8442
-39.8682
-59.4884
-79.0262
X轴正向
9.9243
19.9380
39.8195
59.6146
79.3660
Y轴正向
10.0445
19.9458
39.9983
59.9388
79.2885
Z轴正向
10.0398
19.9238
39.9615
59.7230
79.3576

按公式（2）计算陀螺标度因数、安装误差和零偏结果如下。

2.2加速度计安装误差测试
一、实验步骤
1) 接通电源，预热至IMU工作稳定，启动数据采集软件；
2) 摇动转台手柄使IMU安装台面垂直，顺时针旋转垂直方向的转台，每隔20°作为一个实验测试位置，直到转过360°回到原位置，再依次逆时针旋转垂直方向的转台，分别记录18组Z向加速度计输出数据；
3) 将转台台面调至水平，安装IMU使加速度计的Z向垂直于水平面内。摇动转台手柄使IMU安装台面垂直，调整转台Z向加速度计处于水平位置，此位置记为初始位置。
4) 顺时针旋转垂直方向的转台，每隔20°作为一个实验测试位置，直到转过360°回到原位置，再依次逆时针旋转垂直方向的转台，分别记录18组X、Y向加速度计输出数据。
二、实验结果
5) 按公式（3）计算加速度计标度因数、安装误差和零位误差。
标度因数及安装误差阵

加速度计的安装误差系数和零偏

三，实验分析
1， Matlab在求逆的时候容易出现非奇异，导致不能求逆，应当适当调整输入阵的元素顺序，避免奇异。
2.3旋转调制原理验证实验
一、实验目的
1、通过认识旋转调制技术，实现理论课学习范畴的拓展；
2、验证旋转调制技术的效果，加强学生对旋转调制技术的理解。
二、实验内容
观摩单轴旋转调制系统工作过程，学习旋转调制原理，验证旋转调制技术对陀螺、加速度计性能的调制效果。
三、实验系统组成
旋转调制式捷联惯导系统、稳压电源、数据采集系统。
四、实验原理
旋转调制是陀螺漂移的自补偿技术，设X向陀螺的漂移为，加计零偏为，Y向陀螺的漂移为，加计零偏为,平台绕Z轴以的角速度旋转，如图1所示，则地理坐标系下的等效东向和北向陀螺漂移和加速度计零偏的表达式有：

图1 旋转调制捷联惯导的组成框图
（1）
（2）
旋转调制技术可以将惯性器件引起的误差大大降低。
五、实验内容
1．教师讲解
l 旋转调制技术简介
l 旋转调制系统及显控系统简介
2．操作演示
l 演示旋转调制系统工作过程
l 控制电机分别处于旋转和锁定状态，输出并存储陀螺和加速度计数据。
3．理论探究
对两种状态下的陀螺漂移角速率分别求姿态误差角，探究旋转调制对陀螺性能改善的作用。
旋转调制后 ：

未经调制的误差角：

六、实验分析
由图可以看出，经过旋转调制的姿态误差角从原理上讲应该比未经过旋转调制的姿态误差角小，但是由于其他干扰，只能看出经过旋转调制后相对于未经过旋转调制的情况幅值确定而不是发散。
四，实验源程序
1，陀螺仪安装误差标定源程序
%%%%%%%%%%%X轴反向
x10n=load('E:\惯性器件综合实验\2\陀螺标定数据\x10n.txt');
x20n=load('E:\惯性器件综合实验\2\陀螺标定数据\x20n.txt');
x40n=load('E:\惯性器件综合实验\2\陀螺标定数据\x40n.txt');
x60n=load('E:\惯性器件综合实验\2\陀螺标定数据\x60n.txt');
x80n=load('E:\惯性器件综合实验\2\陀螺标定数据\x80n.txt');
x101=x10n(2:7201,1);
x201=x20n(2:3601,1);
x401=x40n(2:1801,1);
x601=x60n(2:1201,1);
x801=x80n(28:929,1);
x101=mean(x101);
x201=mean(x201);
x401=mean(x401);
x601=mean(x601);
x801=mean(x801);
x1=[x101 x201 x401 x601 x801];
%%%%%%%%%%%Y轴反向
y10n=load('E:\惯性器件综合实验\2\陀螺标定数据\y10n.txt');
y20n=load('E:\惯性器件综合实验\2\陀螺标定数据\y20n.txt');
y40n=load('E:\惯性器件综合实验\2\陀螺标定数据\y40n.txt');
y60n=load('E:\惯性器件综合实验\2\陀螺标定数据\y60n.txt');
y80n=load('E:\惯性器件综合实验\2\陀螺标定数据\y80n.txt');
y101=y10n(2:7201,2);
y201=y20n(2:3601,2);
y401=y40n(2:1801,2);
y601=y60n(2:1201,2);
y801=y80n(28:929,2);
y101=mean(y101);
y201=mean(y201);
y401=mean(y401);
y601=mean(y601);
y801=mean(y801);
y1=[y101 y201 y401 y601 y801];
%%%%%%%%%%%Z轴反向
z10n=load('E:\惯性器件综合实验\2\陀螺标定数据\z10n.txt');
z20n=load('E:\惯性器件综合实验\2\陀螺标定数据\z20n.txt');
z40n=load('E:\惯性器件综合实验\2\陀螺标定数据\z40n.txt');
z60n=load('E:\惯性器件综合实验\2\陀螺标定数据\z60n.txt');
z80n=load('E:\惯性器件综合实验\2\陀螺标定数据\z80n.txt');
z101=z10n(2:7201,3);
z201=z20n(2:3601,3);
z401=z40n(2:1801,3);
z601=z60n(2:1201,3);
z801=z80n(28:929,3);
z101=mean(z101);
z201=mean(z201);
z401=mean(z401);
z601=mean(z601);
z801=mean(z801);
z1=[z101 z201 z401 z601 z801];
%%%%%%%%%%%X轴正向
x10p=load('E:\惯性器件综合实验\2\陀螺标定数据\x10p.txt');
x20p=load('E:\惯性器件综合实验\2\陀螺标定数据\x20p.txt');
x40p=load('E:\惯性器件综合实验\2\陀螺标定数据\x40p.txt');
x60p=load('E:\惯性器件综合实验\2\陀螺标定数据\x60p.txt');
x80p=load('E:\惯性器件综合实验\2\陀螺标定数据\x80p.txt');
x101=x10p(2:7201,1);
x201=x20p(2:3601,1);
x401=x40p(2:1801,1);
x601=x60p(2:1201,1);
x801=x80p(28:929,1);
x101=mean(x101);
x201=mean(x201);
x401=mean(x401);
x601=mean(x601);
x801=mean(x801);
x0=[x101 x201 x401 x601 x801];
%%%%%%%%%%%%%%%%Y轴正向%%%%%%%%%%%%%%%%%%%%%
y10p=load('E:\惯性器件综合实验\2\陀螺标定数据\y10p.txt');
y20p=load('E:\惯性器件综合实验\2\陀螺标定数据\y20p.txt');
y40p=load('E:\惯性器件综合实验\2\陀螺标定数据\y40p.txt');
y60p=load('E:\惯性器件综合实验\2\陀螺标定数据\y60p.txt');
y80p=load('E:\惯性器件综合实验\2\陀螺标定数据\y80p.txt');
y101=y10p(2:7201,2);
y201=y20p(2:3601,2);
y401=y40p(2:1801,2);
y601=y60p(2:1201,2);
y801=y80p(28:929,2);
y101=mean(y101);
y201=mean(y201);
y401=mean(y401);
y601=mean(y601);
y801=mean(y801);
y0=[y101 y201 y401 y601 y801];
%%%%%%%%%%%%%%%%Z轴正向%%%%%%%%%%%%%%%%%%%%%
z10p=load('E:\惯性器件综合实验\2\陀螺标定数据\z10p.txt');
z20p=load('E:\惯性器件综合实验\2\陀螺标定数据\z20p.txt');
z40p=load('E:\惯性器件综合实验\2\陀螺标定数据\z40p.txt');
z60p=load('E:\惯性器件综合实验\2\陀螺标定数据\z60p.txt');
z80p=load('E:\惯性器件综合实验\2\陀螺标定数据\z80p.txt');
z101=z10p(2:7201,3);
z201=z20p(2:3601,3);
z401=z40p(2:1801,3);
z601=z60p(2:1201,3);
z801=z80p(28:929,3);
z101=mean(z101);
z201=mean(z201);
z401=mean(z401);
z601=mean(z601);
z801=mean(z801);
z0=[z101 z201 z401 z601 z801];
k=[1 1 1 1 1 ]
x=[x1,x0];
y=[y1,y0];
z=[z1,z0];
W1=[10 20 40 60 80 -10 -20 -40 -60 -80;
80 60 40 20 10 -10 -20 -40 -60 -80;
60 40 20 10 80 -10 -20 -40 -60 -80;
1 1 1 1 1 1 1 1 1 1 ];
W0=[ 9.9243 19.9380 39.8195 59.6146 79.3660 -9.9123 -19.9246 -39.8253 -59.6601 -79.1389;
79.2885 59.9388 39.9983 19.9458 10.0445 -10.3317 -19.9998 -40.4727 -59.7430 -79.2054;
59.7230 39.9615 19.9238 10.0398 79.3576 -9.9690 -19.8442 -39.8682 -59.4884 -79.0262]
Sg= W0*(W1)'*inv(W1*(W1)');

2，计算加速度计安装误差源程序
xd=load('E:\惯性器件综合实验\3\accData\xd.txt');
xu=load('E:\惯性器件综合实验\3\accData\xu.txt');
yd=load('E:\惯性器件综合实验\3\accData\yd.txt');
yu=load('E:\惯性器件综合实验\3\accData\yu.txt');
zd=load('E:\惯性器件综合实验\3\accData\zd.txt');
zu=load('E:\惯性器件综合实验\3\accData\zu.txt');
xd1=mean(xd)/1000;
xu1=mean(xu)/1000;
yd1=mean(yd)/1000;
yu1=mean(yu)/1000;
zd1=mean(zd)/1000;
zu1=mean(zu)/1000;
xd11=xd1(4:6)';
xu11=xu1(4:6)';
yd11=yd1(4:6)';
yu11=yu1(4:6)';
zd11=zd1(4:6)';
zu11=zu1(4:6)';
Ao=[xd11 yd11 zd11 xu11 yu11 zu11]
AI=[1 0 0 -1 0 0;
0 1 0 0 -1 0;
0 0 1 0 0 -1;
1 1 1 1 1 1];
C=inv(AI*AI');
Sacce=Ao*AI'*C;

旋转调制与未经过旋转调制的源程序
Q=load('E:\惯性器件综合实验\惯性导航试验数据\2\旋转调制实验数据\data.xls');
%%%%%%%求转换矩阵%%%%%%%%
Wx1=Q(160002:170002,1);
Wy1=Q(160002:170002,2);%%%提取转换矩阵的数据
Wz1=Q(160002:170002,3);
Wx1=mean(Wx1)*0.56534653/3600*1000*pi/180 ; %%%%转换为度每秒
Wy1=mean(Wy1)*0.57063519/3600*1000*pi/180 ;
Wz1=mean(Wz1)*(-0.57213617)/3600*1000*pi/180;
W=[Wx1 Wy1 Wz1];
Ax=Q(160002:170002,4);
Ay=Q(160002:170002,5);
Az=Q(160002:170002,6);
Ax=mean(Ax)*(3.1507301e-4)*1000;
Ay=mean(Ay)*(3.1618922e-4)*1000;
Az=mean(Az)*(3.1207785e-4)*1000;
A=[Ax Ay Az];
g=9.8;
wie=7.292115147e-5;
L = 40.0211142228246*pi/180;
C3=A/g;

c3=C3';
C12=(Wx1-C3(1)*wie*sin(L))/(wie*cos(L));
C22=(Wy1-C3(2)*wie*sin(L))/(wie*cos(L));
C32=(Wz1-C3(3)*wie*sin(L))/(wie*cos(L));
C2=[C12 C22 C32];
c2=C2';
C1=cross(C2,C3);
c=C1';
Ctb=[c,c2,c3];
%%%%%%%%%%%%旋转调制转换%%%%%%%%
%%%%%%%%%%%%%%%%提取数据%%%%%%%%%%%%%%%%%%
Wx=Q(32001:160001,1);
Wy=Q(32001:160001,2);%%%提取旋转调制的数据
Wz=Q(32001:160001,3);
JD=Q(32001:160001,7);
Wx=Wx*0.56534653/3600*1000*pi/180 ; %%%%转换为hu度每秒
Wy=Wy*0.57063519/3600*1000*pi/180 ;
Wz=Wz*(-0.57213617)/3600*1000*pi/180;
D=JD*0.00022865853658536584*pi/180; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5转化为弧度
wie=7.292115147e-5;
L = 40.0211142228246*pi/180;
wibbx=Wx.*cos(D)+Wy.*sin(D);
wibby=Wx.*sin(D)+Wy.*cos(D);
wibbz=Wz;
Wibb=[wibbx wibby wibbz];
wittx=0;
witty=wie*cos(L);
wittz=wie*sin(L);
Witt=[wittx witty wittz]';
Witb=Ctb* Witt;
wtbbx= wibbx-Witb(1);
wtbby= wibby-Witb(2);
wtbbz=wibbz-Witb(3);
nx=zeros(160001-32001,1);
ny=zeros(160001-3200,1);
nz=zeros(160001-32001,1);
nx(1)=wtbbx(1,1)*1/1000;
ny(1)=wtbby(1,1)*1/1000;
nz(1)=wtbbz(1,1)*1/1000;
for i=2:(160001-32001)
nx(i)=nx(i-1)+wtbbx(i,1)*1/1000;
ny(i)=ny(i-1)+wtbby(i,1)*1/1000;
nz(i)=nz(i-1)+wtbbz(i,1)*1/1000;
u(i)=i;
end
plot(u,nx);
plot(u,ny);
plot(u,nz);
%%%%%%%%%%%%%%%%%%%%%%%%%%%未经过旋转调制的果%%%%%%%%%%%
g=9.8;
wie=7.292115147e-5;
L = 40.0211142228246*pi/180;
Wx2=Q(160001:225458,1);
Wy2=Q(160001:225458,2);
Wz2=Q(160001:225458,3);
JD1=Q(160001:225458,7);
%%%%%%%%%%%%5转换单位%%%%%%%%
Wxx=Wx2*0.56534653/3600*1000*pi/180 ; %%%%转换为度每秒
Wyy=Wy2*0.57063519/3600*1000*pi/180 ;
Wzz=Wz2*(-0.57213617)/3600*1000*pi/180;
D=JD1*0.00022865853658536584*pi/180; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5转化为弧度
wibbx1=Wxx;
wibby1=Wyy;
wibbz1=Wzz;
Wibb=[wibbx1 wibby1 wibbz1];
%%%%%%%%%%%%%求转换矩阵
wittx=0;
witty=wie*cos(L);
wittz=wie*sin(L);
Witt=[wittx witty wittz]';
Witb=Ctb* Witt;
wtbbx1= wibbx1-Witb(1);
wtbby1= wibby1-Witb(2);
wtbbz1=wibbz1-Witb(3);
nx1(1)=wtbbx1(1,1)*1/1000;
ny1(1)=wtbby1(1,1)*1/1000;
nz1(1)=wtbbz1(1,1)*1/1000;
for i=2:(225458-160001)
nx1(i)=nx1(i-1)+wtbbx1(i,1)*1/1000;
ny1(i)=ny1(i-1)+wtbby1(i,1)*1/1000;
nz1(i)=nz1(i-1)+wtbbz1(i,1)*1/1000;
u(i)=i;
end
plot(u,ny1);