52matlab技术网站，matlab教程，matlab安装教程，matlab下载

标题: 捷联惯导程序 [打印本页]

作者: matlab的旋律 时间: 2020-5-15 05:24
标题: 捷联惯导程序
捷联惯导程序，依据加表和陀螺仪的输出数据来求解飞行器的姿态

clc;
clear;
format long;  %设置数据精度为15位小数

Data=importdata('temp.txt'); % 导入实验所采集的数据，以矩阵形式赋给Data变量，temp.txt必须与该M文件在同一个文件夹中
Px=Data(:,3);  % Px,Py,Pz为陀螺仪的输出值
Py=Data(:,4);
Pz=Data(:,5);

Nx=Data(:,6);  % Nx,Ny，Nz为加速度计的输出值
Ny=Data(:,7);
Nz=Data(:,8);
% 陀螺仪模型参数标定如下：
Sx  = -4.085903e-006 ; Sy  = -4.085647e-006 ; Sz  = -4.085170e-006 ;
Mxy = 5.059527e-003  ; Mxz = -1.031103e-003 ; Myx = -3.355451e-003 ;
Myz = 3.508468e-003  ; Mzx = -1.266671e-003 ; Mzy = -2.318244e-004 ;
Dx  = -2.009710e-006 ; Dy  = 8.156346e-007  ; Dz  = -5.749059e-007 ;
GyroCali_A = [ 1 -Mxy -Mxz ; -Myx 1 -Myz ; -Mzx -Mzy 1 ];
% 加速度计模型参数标定如下：
Kx  = 9.272930e-004  ; Ky  = 9.065544e-004  ; Kz  = 9.443748e-004  ;
Ixy = 6.533872e-003  ; Ixz = 9.565992e-004  ; Iyx = -6.319376e-003 ;
Iyz = -6.902339e-004 ; Izx = -1.144549e-003 ; Izy = -3.857963e-004 ;
Bx  = -3.400847e-002 ; By  = -8.916341e-003 ; Bz  = -9.947414e-003 ;
AccCali_A = [1 -Ixy -Ixz ; -Iyx 1 -Iyz ; -Izx -Izy 1 ];

Delta_t = 0.05; %采样时间为0.05秒

Delta_Theta_x = 0;
Delta_Theta_y = 0;
Delta_Theta_z = 0;  %定义陀螺仪输出的角度增量

Delta_Vx = 0;
Delta_Vy = 0;
Delta_Vz = 0;    %定义加速度计输出的速度增量

L = zeros(1,12001);
L(1)= 45.7328*pi/180 ; %纬度用L表示，纬度的初始值划为弧度形式，因为后面计算位置矩阵更新
L(2)= 45.7328*pi/180 ; %时需要用到前两次的L值来计算当前L值，所以在此定义2个初始L值
Lamda = 126.6287*pi/180 ; %经度用Lamda表示，经度的初始值划为弧度形式
h = 136 ;                %高度用h表示
V = [ 0 ; 0 ; 0 ];  %导航坐标系中的东北天初始速度都为0
Vx = 0;          %方便后面的速度计算与速度更新
Vy = 0;
Vz = 0;

Theta = 0;
Gama  = 0;
Fai = 0;       %初始姿态角（俯仰角/倾斜角/航向角）都为0,此处均为弧度

Re = 6378254 ;  Rp = 6356803 ;%定义地球的半长轴与半短轴
e  = (Re - Rp)/Re ;          %定义旋转椭球扁率（椭球度）
Wie = 15.04107/180*pi ;    %定义地球自转角速度，地球坐标系相对于惯性坐标系的角速度

Theta_Matrix = zeros(1,12000); %定义姿态角矩阵，供画图用
Gama_Matrix  = zeros(1,12000);
Fai_Matrix = zeros(1,12000);
L_Matrix    = zeros(1,12001); %定义经纬度矩阵，供画图用，L的特殊性决定了其数据个数为12001
L_Matrix(1)  = 45.7328;
Lamda_Matrix = zeros(1,12000);
Ve_Matrix = zeros(1,12000); %定义速度矩阵，供画图用
Vn_Matrix = zeros(1,12000);
Vu_Matrix = zeros(1,12000);

%以下计算捷联矩阵的初始值，捷联矩阵的初始值仅仅由Theta，Gama，Fai的初始值决定
T = [ cos(Gama)*cos(Fai)-sin(Gama)*sin(Theta)*sin(Fai)  -cos(Theta)*sin(Fai)  sin(Gama)*cos(Fai)+cos(Gama)*sin(Theta)*sin(Fai) ;
      cos(Gama)*sin(Fai)+sin(Gama)*sin(Theta)*cos(Fai) cos(Theta)*cos(Fai)  sin(Gama)*sin(Fai)-cos(Gama)*sin(Theta)*cos(Fai) ;
            -sin(Gama)*cos(Theta)                            sin(Theta)                cos(Gama)*cos(Theta)                ];
%由捷联矩阵的初始值计算初始四元数值，为捷联矩阵的实时更新做准备
if(T(3,2)-T(2,3)>0)
      Q1 = 0.5*sqrt(1+T(1,1)-T(2,2)-T(3,3));
else if (T(3,2)-T(2,3)==0)
      Q1 = 0;
   else Q1 = -0.5*sqrt(1+T(1,1)-T(2,2)-T(3,3));
      end          %求解Q1
end

if(T(1,3)-T(3,1)>0)
      Q2 = 0.5*sqrt(1-T(1,1)+T(2,2)-T(3,3));
else if (T(1,3)-T(3,1)==0)
   Q2 = 0;
   else Q2 = -0.5*sqrt(1-T(1,1)+T(2,2)-T(3,3));
      end          %求解Q2
end

if(T(2,1)-T(1,2)>0)
      Q3 = 0.5*sqrt(1-T(1,1)-T(2,2)+T(3,3));
   else if (T(2,1)-T(1,2)==0)
   Q3 = 0;
      else Q3 = -0.5*sqrt(1-T(1,1)-T(2,2)+T(3,3));
      end          %求解Q3
end

Q0 = 0.5*sqrt(1-Q1*Q1-Q2*Q2-Q3*Q3); %求解Q0
Q = [Q0 ; Q1 ; Q2 ; Q3];  %四元数初始值
Q = Q / sqrt(Q0*Q0+Q1*Q1+Q2*Q2+Q3*Q3); %四元数的初始归一化，为得到最小漂移误差

%以下求位置矩阵的初始值，通过位置矩阵更新后，反过来算运载体所在的经纬度
%位置矩阵仅仅与经纬度有关系，Ce2n表示把地球坐标系转换为导航坐标系的转换矩阵
Ce2n = [    -sin(Lamda)          cos(Lamda)             0 ;
         -sin( L(1) )*cos(Lamda) -sin(L(1))*sin(Lamda)    cos( L(1) );
            cos( L(1) )*cos(Lamda) cos( L(1) )*sin(Lamda) sin( L(1) )    ];

%大循环，共执行12000次，实时更新捷联矩阵，速度矩阵，位置矩阵，保存作图所需数据
for k = 1:12000;

GyroCali_B = [Sx*Px(k)-Dx*Delta_t ; Sy*Py(k)-Dy*Delta_t ; Sz*Pz(k)-Dz*Delta_t ];
Delta_Theta = GyroCali_A * GyroCali_B ; %计算陀螺仪输出的角度增量
Delta_Theta_x = Delta_Theta(1);
Delta_Theta_y = Delta_Theta(2);
Delta_Theta_z = Delta_Theta(3);
Delta_Theta_Module = sqrt( Delta_Theta_x * Delta_Theta_x + Delta_Theta_y * Delta_Theta_y + Delta_Theta_z * Delta_Theta_z );

AccCali_B = [Kx*Nx(k)-Bx*Delta_t ; Ky*Ny(k)-By*Delta_t ; Kz*Nz(k)-Bz*Delta_t ];
Delta_V = AccCali_A * AccCali_B ;    %计算加速度计输出的速度增量
Delta_Vx = Delta_V(1);
Delta_Vy = Delta_V(2);
Delta_Vz = Delta_V(3);
Delta_V_Module = sqrt( Delta_Vx * Delta_Vx + Delta_Vy * Delta_Vy + Delta_Vz * Delta_Vz );

%使用毕卡法求解四元数更新矩阵，即捷联矩阵
Bika = zeros(4);
Bika(1,1) = cos(0.5 * Delta_Theta_Module);
Bika(1,2) = -Delta_Theta_x / Delta_Theta_Module * sin(0.5 * Delta_Theta_Module);
Bika(1,3) = -Delta_Theta_y / Delta_Theta_Module * sin(0.5 * Delta_Theta_Module);
Bika(1,4) = -Delta_Theta_z / Delta_Theta_Module * sin(0.5 * Delta_Theta_Module);
Bika(2,1) = -Bika(1,2);
Bika(2,2) =  Bika(1,1);
Bika(2,3) = -Bika(1,4);
Bika(2,4) =  Bika(1,3);
Bika(3,1) = -Bika(1,3);
Bika(3,2) = -Bika(2,3);
Bika(3,3) =  Bika(1,1);
Bika(3,4) = -Bika(1,2);
Bika(4,1) = -Bika(1,4);
Bika(4,2) = -Bika(2,4);
Bika(4,3) = -Bika(3,4);
Bika(4,4) =  Bika(1,1);
Q = Bika * Q; % 每循环一次，更新一次四元素Q值，为求捷联矩阵
Q = Q / sqrt(Q0*Q0+Q1*Q1+Q2*Q2+Q3*Q3);  %四元数的归一化，为得到最小漂移误差
Q0 = Q(1);
Q1 = Q(2);
Q2 = Q(3);
Q3 = Q(4);

%捷联矩阵的四元数表达式
T = [ Q0*Q0+Q1*Q1-Q2*Q2-Q3*Q3    2*(Q1*Q2-Q0*Q3)    2*(Q1*Q3+Q0*Q2)
         2*(Q1*Q2+Q0*Q3)    Q0*Q0-Q1*Q1+Q2*Q2-Q3*Q3 2*(Q2*Q3-Q0*Q1)
         2*(Q1*Q3-Q0*Q2)    2*(Q2*Q3+Q0*Q1)    Q0*Q0-Q1*Q1-Q2*Q2+Q3*Q3 ];
%*********************************************************************
%********************求三个姿态角Theta，Gama和 Fai ********************
%*********************************************************************
Theta_Main = asin( T(3,2) );
Gama_Main = atan( -T(3,1) / T(3,3));
Fai_Main  = atan( -T(1,2) / T(2,2));
Theta = Theta_Main;

if (T(3,3)>0)
      Gama = Gama_Main ;
else if (T(3,3)<0 && Gama_Main > 0)
         Gama = Gama_Main + pi;
      else  Gama = Gama_Main - pi;  %此处用else实为不妥，不过为了程序的完善性，只能这样了
      end
end

if ( T(2,2)<0 )
      Fai = Fai_Main + pi ;
else if (T(2,2)==0)
         Fai = pi/2;
      else if ( Fai_Main>0)
         Fai = Fai_Main;
         else Fai = Fai_Main + 2*pi ;
         end
      end
end

%以下存储姿态角到三个矩阵里面，为画图做准备
Theta_Matrix(k) = Theta*180/pi; %作图用矩阵,以角度表示
Gama_Matrix(k) = Gama*180/pi;    %作图用矩阵,以角度表示
if (Fai<2*pi)
      Fai_Matrix(k) = Fai*180/pi;
else  Fai_Matrix(k) = Fai*180/pi-360;  %作图用矩阵,以角度表示
end
%到此为止，姿态角的求解完毕，以下先求速度
%*********************************************************************
%********************求飞行器相对于东北天的速度*************************
%*********************************************************************
Rm = Re*( 1-2*e+3*e*sin( L(k+1) )^2 );
Rn = Re*( 1+e*sin(L(k+1))^2 );
LL = 3/2*L(k+1) - 1/2*L(k) ;
F = [ 0  -1/(Rm + h)  0 ; 1/(Rn + h) 0  0 ; tan( LL )/(Rn + h)  0  0];
g = 9.7803+0.051799*sin(L(k+1))^2-0.94114e-006*h;
G = [0;0;-g];
Wen2n = F*V;
Wie2n = [0 ; Wie*cos(L(k+1)); Wie*sin(L(k+1))];
W = 2*Wie2n + Wen2n;
W_X = [ 0 -W(3) W(2) ; W(3)  0  -W(1) ; -W(2)  W(1) 0 ]; %此式中W_X为 2*Wie+Wen2n 的反对称矩阵
V = V + T*( Delta_V+ [ 0 -Delta_Theta_z Delta_Theta_y ; Delta_Theta_z  0  -Delta_Theta_x ; -Delta_Theta_y  Delta_Theta_x 0 ]*Delta_V ) + Delta_t*(G-W_X*V);

Vx = V(1);
Vy = V(2);
Vz = V(3);
Ve_Matrix(k) = V(1);
Vn_Matrix(k) = V(2);
Vu_Matrix(k) = V(3);

%*********************************************************************
%********************求飞行器所在的经纬度******************************
%*********************************************************************

Epsilon = F*V*Delta_t;
Ce2n =( eye(3) - [ 0 -Epsilon(3) Epsilon(2) ; Epsilon(3)  0  -Epsilon(1) ; -Epsilon(2)  Epsilon(1) 0 ] )*Ce2n ; %位置矩阵实时更新
%下面通过位置矩阵来实时更新经纬度
L(k+2) = asin( Ce2n(3,3)); %由于L本来就是以矩阵形式定义的，下面定义一个把L用角度表示的矩阵
L_Matrix(k+1) = L(k+2)*180/pi;
Lamda_Main = atan( Ce2n(3,2)/Ce2n(3,1) ); %计算出来的L和Lamda都是弧度制的
if (Ce2n(3,1)>0)
      Lamda = Lamda_Main;
else if (Lamda_Main<0)
         Lamda = Lamda_Main + pi;
      else Lamda = Lamda_Main - pi;
      end
end
Lamda_Matrix(k) = Lamda*180/pi;  %作图用矩阵，以角度表示

end
%**************************************************************************
%********************以下是画图程序*****************************************
%**************************************************************************
k=1:1:12000;             %绘制三轴姿态变化图线-绿色
figure(1);
plot(k/20,Theta_Matrix(k),'g');
xlabel('Time(second)');
ylabel('Angle(degree)');
title('Theta(俯仰角)');
grid on;
figure(2);
plot(k/20,Gama_Matrix(k),'g');
xlabel('Time(second)');
ylabel('Angle(degree)');
title('Gama（滚转角）');
grid on;
figure(3);
plot(k/20,Fai_Matrix(k),'m');
xlabel('Time(second)');
ylabel('Angle(degree)');
title('Fai（偏航角）');
grid on;
                        %绘制东北天个方向的速度变化曲线-红色
figure(4);
plot(k/20,Ve_Matrix(k),'r');
xlabel('Time(second)');
ylabel('Velocity(m/s)');
title('Ve（东向速度）');
grid on;
figure(5);
plot(k/20,Vn_Matrix(k),'r');
xlabel('Time(second)');
ylabel('Velocity(m/s)');
title('Vn（北向速度）');
grid on;
figure(6);
plot(k/20,Vu_Matrix(k),'r');
xlabel('Time(second)');
ylabel('Velocity(m/s)');
title('Vu（天向速度）');
grid on;
                        %绘制飞行器所在经纬度曲线-蓝色
figure(7);
plot(k/20,L_Matrix(k),'b');
xlabel('Time(second)');
ylabel('Degree');
title('Latitude L(纬度)');
grid on;

figure(8);
plot(k/20,Lamda_Matrix(k),'b');
xlabel('Time(second)');
ylabel('Degree');
title('Longitude Lamda(经度)');
grid on;

作者: matlab的旋律 时间: 2020-5-15 05:25
捷联惯导程序，依据加表和陀螺仪的输出数据来求解飞行器的姿态

作者: 豹纹096 时间: 2020-6-7 07:23

欢迎光临 52matlab技术网站，matlab教程，matlab安装教程，matlab下载 (http://www.52matlab.com/)