Abstract
Basic principles of the strapdown inertial navigation system (SINS) using the outputs of strapdown gyros and accelerometers are explained, and the main equations are described. A mathematical model of SINS is established, and its Matlab implementation is developed. The theory is illustrated by six examples which are static status, straight line movement, circle movement, s-shape movement, and two sets of real static data.
1. Introduction
Many navigation books and papers on inertial navigation system (INS) provide readers with the basic principle of INS. Some also superficially describe simulation methods and rarely provide the free code which can be used by new INS users to help them understand the theory and develop INS applications. Commercial simulation software is available but is not free. The objective of this paper is to develop an easy-to-understand step-by-step development method for simulating INS. Here we consider the most popular INS which is the strapdown inertial navigation system (SINS). The mathematical operations required in our work are mostly matrix manipulations and more generally basic linear algebra [1]. In this paper, Matlab [2] is chosen as the simulation environment. It is a popular computing environment to perform complex matrix calculations and to produce sophisticated graphics in a relatively easy manner. A large collection of Matlab scripts are now available for a wide variety of applications and are often used for university courses. Matlab is also becoming more and more popular in industrial research centers in the design and simulation stages.
The main purposes of this paper are to establish a mathematical model and to develop a comprehensive Matlab implementation for SINS. The structure of the proposed mathematical model and Matlab simulation of SINS is shown in Figure 1. In Section 2, the INS-related orthogonal coordinates (the body frame, the inertial frame, the Earth frame, the navigation frame, the ENU-frame, and the wander azimuth navigation frame) are described and figures to illustrate the relationship between the frames are provided. The basic principle of SINS is described in the wander azimuth navigation frame (-frame). In Section 3, two important direction cosine matrices (DCMs), the vehicle attitude DCM and the position DCM, and the related important attitude and position angles are defined. In Section 4, the simulation for data generation of gyros and accelerometers is described in ENU-frame. Instead of -frame, ENU-frame is chosen because the outputs of gyros and accelerometers are easier to obtain in this frame. The Matlab implementation is given and described step by step. Four kinds of scenarios (static, straight, circle, and s-shape) are set as examples of different kinds of vehicle trajectories. In Section 5, the mathematical model of SINS is set up and the calculation steps in -frame are provided. In Section 6, the required initial parameters and other initial data calculation for the SINS model are given for the different simulation scenarios. In Section 7, Matlab implementation code functions are listed and described. Further, simulation results for the four above-mentioned scenarios are presented; two examples from real SINS experiment data are also provided to verify the validity of the developed codes. Finally, conclusions are drawn. Mathematical details are given in Appendices A–D.
2. Principles
A fundamental aspect of inertial navigation is the precise definition of a number of Cartesian coordinate reference frames. Each frame is an orthogonal, right-handed, coordinate frame or axis set. For all the coordinate frames used in this paper, a positive rotation about each axis is taken to be in a clockwise direction looking along the axis from the origin, as indicated in Figure 2. A negative rotation corresponds to an anti-clockwise direction. This convention is used throughout this paper. It is also worth pointing out that a change in attitude of a body, which is subjected to a series of rotations about different axes, is not only a function of the rotation angles, but also on the order in which the rotations occur. In this paper, the following coordinate frames are used [3].(1)The body frame (-frame): the -frame, depicted in Figure 2, is an orthogonal axis set which has its origin at the center of the vehicle, point , and is aligned with the pitch axis, roll axis, and yaw axis of the vehicle in which the navigation system is installed.(2)The inertial frame (-frame): the -frame, depicted in Figure 3, has its origin at the center of the Earth and its axes nonrotating with respect to fixed stars; these axes are denoted by , , and , with being coincident with the Earth polar axis.(3)The Earth frame (-frame): the -frame, depicted in Figure 3, has its origin at the center of the Earth and axes nonrotating with respect to the Earth; these axes are denoted by , , and . The axis is the Earth polar axis. The axis is along the intersection of the plane of the Greenwich meridian and the Earth equatorial plane. The Earth frame rotates with respect to the inertial frame at a rate about the axis .(4)The navigation frame (-frame): the -frame, depicted in Figure 3, is a local geographic navigation frame which has its origin at the location of the navigation system, point (the navigation system is fixed inside the vehicle and we assume that the navigation system is located exactly at the center of the vehicle), and axes aligned with the directions of east , north and the local vertical up . When the -frame is defined in this way, it is called the “ENU-frame.” The turn rate of the navigation frame with respect to the Earth-fixed frame, , is governed by the motion of the point with respect to the Earth. This is often referred to as the transport rate.(5)The wander azimuth navigation frame (-frame): the -frame, depicted in Figure 3, may be used to avoid the singularities in the computation which occur at the poles of the -frame. Like the -frame, it is of a local level but rotates through the wander angle about the local vertical. Here we do not call this frame -frame ( for wander) for notation clarity since and may look similar when printed. Letter in -frame stands for platform; indeed the wander azimuth navigation frame is of a local level and thus forms a horizontal platform.
In this paper, we choose the -frame as the navigation frame for vehicle trajectory calculation, for the following reason. In the local geographic navigation frame mechanization, the -frame is required to rotate continuously as the system moves over the surface of the Earth in order to keep its axis pointing to the true north. In order to achieve this condition worldwide, the -frame must rotate at much greater rates about its axis as the navigation system moves over the surface of the Earth in the polar regions, compared to the rates required at lower latitudes. It is clear that near the polar areas the local geographic navigation frame must rotate about its axis rapidly in order to maintain the axis pointing to the pole. The heading direction will abruptly change by when moving past the pole. In the most extreme case, the turn rate becomes infinite when passing over the pole. One way of avoiding the singularity, and also providing a navigation system with worldwide capability, is to adopt a wander azimuth mechanization in which the -component of is always set to zero, that is, . A wander axis system is a local level frame which moves over the Earth surface with the moving vehicle, as depicted in Figure 3. However, as the name implies, the azimuth angle between axis and axis varies with the vehicle position on Earth. This variation is chosen in order to avoid discontinuities in the orientation of the wander frame with respect to Earth as the vehicle passes over either the north or south pole.
In the remainder of this section, the main principle of SINS in the -frame is described.
Along the same lines as in [3], a navigation equation for a wander azimuth system can be constructed as follows: where is the direction cosine matrix used to transform the measured specific force vector in -frame into -frame. This matrix propagates in accordance with the following equation: where is the skew symmetric form of , the -frame angular rate with respect to the -frame.
Equation (1) is integrated to generate estimates of the vehicle speed in the wander azimuth frame, . This is then used to generate the turn rate of the wander frame with respect to the Earth frame, . The direction cosine matrix which relates the wander frame to the Earth frame, , may be updated using the following equation: where the superscript means matrix transposition.
Since , and skew symmetric matrix is (see Appendix A), (4) can be rewritten as where is a skew symmetric matrix formed from the elements of the angular rate vector ; we could have when the rotation angles are reciprocal. Because the -component of is set to zero, , the matrix expression of is . This process is implemented iteratively and enables any singularities to be avoided.
In the next section, the two important DCMs, the vehicle attitude DCM and vehicle position DCM, are defined, as well as the vehicle-attitude-related attitude angles and vehicle-position-related position angles.
3. Direction Cosine Matrices (DCMs)
In this section, the vehicle attitude DCM with the corresponding attitude angles and the vehicle position DCM with the corresponding position angles are described separately.
3.1. Vehicle Attitude DCM
The definition of the rotation sequence from -frame to -frame is (see Figure 4) where is the gird azimuth angle (0–360°), is the pitch angle (−90°–90°), and is the roll angle (−180°–180°). The above rotation can be written in the following matrix form:
The vehicle attitude DCM is then obtained as
For the -frame system, the angle between the grid north and the true north is the wander azimuth angle . So the angle between the horizontal projection along axis of the vehicle’s vertical axis and the real north is the heading angle . We have that So the direction cosine matrix from -frame to -frame is
The gimbal angles , , and and the gimbal rates , , and are related to the body rate , which is the turn rate of the -frame with respect to -frame and measured in -frame as follows:
3.2. Vehicle Position DCM
Position matrix is the DCM from -frame to -frame. It has the following rotating sequence (see Figure 5): where is the longitude angle (−180°–180°), is the latitude angle (−90°–90°), and is the wander azimuth angle (0–360°). The above rotation can be written in the following matrix form:
In Section 4, a trajectory simulation method in the ENU-frame is described step by step to generate sensor data. In Section 5, a trajectory and attitude simulator method in the -frame is described step by step to derive the desired trajectory and attitude from the simulated sensor data or real sensor data; Section 6 provides the initial parameters and initial data calculation.
4. Sensor Data Generator
The Purpose of Trajectory Simulation is to Generate Data of the 3 Orthogonal Gyros and the 3 Orthogonal Accelerometers According to the Designed Trajectory. It is Mentioned in Section 2 That -Frame is Set up to Avoid the Singularities When the Vehicle Passes Over Either the North or South Pole. But in Most Applications, The SINS Systems Are Seldom Operated under This Extreme Environment. The ENU-Frame Can be Implemented Easier Than -Frame, so it is Chosen as the Navigation Frame. Figure 6 Shows the Whole Process of the SINS Principal in the ENU-Frame Mechanization. First, The Vehicle Trajectory in the ENU-Frame is Set. Then, The Sensor Ideal Output is Derived using the Inverse Principle of INS. The Sensor Simulation Data Can be Obtained by Adding Noise to the Ideal Data. Then, we Use the Simulated Sensor Data to Derive the Noise-Corrupted Simulated Trajectory. Besides, The Difference Between the Ideal and Simulated State Vectors Can be Set as the Input for the Observed Measurements in the Kalman Filter.
4.1. The Initial Parameters
For the designed trajectory, the initial parameters are(1)initial position, latitude , longitude , height ;(2)initial velocity ;(3)the designed variation of acceleration , which varies with time according to the designed trajectory;(4)the designed variations of the attitude angles, pitch , roll , and heading , and attitude angle rates, , , and , which vary with time according to the designed trajectory.
4.2. The Update of Velocity
The velocity is updated as
where is the time step.
4.3. The Update of Position
The position is updated as
4.4. The Update of Attitude
The attitude angles are updated as
The attitude rates are updated as
The expressions for , , , , , and depend on the designed trajectory.
The direction cosine matrix can be calculated using matrix expression (10). We have that .
4.5. Gyro Data Generator
The output of the gyros is
where is the simulated actual output, is the unit matrix, is the diagonal matrix whose diagonal elements correspond to the 3 gyros’ scale factor errors, and is the gyro’s drift and can be simulated as the sum of a constant noise and a random white noise:
In a static base, is equal to zero, whereas, in a moving base it is obtained as
is related to velocity and can be expressed as
4.6. Accelerometer Data Generator
The measurement of the accelerometer is the specific force:
where is the simulated actual output, is the unit matrix. is the diagonal matrix whose diagonal elements correspond to the 3 accelerometers’ scale factor errors, is the bias considered as the sum of a constant noise and a random white noise . , and , where is the 9th element of and is the vehicle altitude.
4.7. Examples
For four examples of static, straight line, circle, and s-shape situations, details will be given next under the conditions that the vehicle is moving on the surface of the Earth with no attitude change except for the heading angle, which means that the pitch angle, roll angle, and altitude are constants during the simulation process:
The calculation method for the other parameters for the four situations is described as follows.(1)Static:(2)Straight line:(3)Circle:(4)S-shape:
5. Mathematical Model and Trajectory Calculation Steps
After the Gyro and Accelerometer Data Are Simulated using the Method Described in the Previous Section under the Designed Scenario, The Next Step we Have to do is to Figure Out the Mathematical Model of SINS and the Calculation Steps to Process the Sensor Data to Get the Calculated Trajectories. Based on the Basic Principles of Strapdown Inertial Navigation System [4], we Draw the Mathematical Model in the -Frame Mechanization in Figure 7. The Calculation Steps Are Described Below. Although the Situation That the Vehicle Passes Over Either the North or South Pole Seldom Happens, The Universal -Frame is Still Chosen Instead of the Simpler ENU-Frame to Give a Navigation Illustration in a Different Frame.
5.1. Quaternion Update and Optimal Normalization
There are three kinds of strapdown attitude representations: DCM, Euler angle, and quaternion. In this paper, we choose quaternion. The reason why quaternion is chosen is explained in [3].
The quaternion formed by a rotating body frame around the platform frame is The update for the quaternion can be obtained by solving the following quaternion differential equation: Based on the Euclide norm minimized indicator [4], the optimal normalization for the quaternion is
5.2. Calculation
is vehicle attitude DCM which transforms the measured angle in the -frame to the -frame, with its 9 components . (Here we use to distinguish it from the components of which is used below.)
After obtaining , , , and using (29), can be calculated as
5.3. Specific Force Transformation from in -Frame to in -Frame
The specific force in the -frame can be transformed to in the -frame by multiplication with DCM :
5.4. Velocity Calculation
The velocity update can be obtained by solving the following differential equation: The ground speed is the vehicle velocity projection on the horizontal plane:
5.5. Position Matrix Update
The update for the position matrix can be obtained by solving the following differential equation, noticing that :
5.6. Position Angular Velocity Update
In the chosen wander azimuth navigation frame, we have , and where where the elements of position matrix can be obtained using (35).
5.7. Earth Angular Velocity and Attitude Angular Velocity Calculation
We Have That
5.8. Attitude Angle Calculation
The relation between attitude matrix and the three attitude angles, grid azimuth angle , pitch angle , and roll angle , is
Thus, the principal values of , , and are Considering the defined range of the angles, the expressions of the real values of , , , and are
5.9. Position Angle Calculation
The relation between position matrix and the 3 position angles, longitude , latitude , and wander azimuth angle , is
Thus, the principal values of , , and are Considering the defined range of the angles, the expressions of the real values of , , and are
5.10. Heading Angle Calculation
The heading angle is calculated as To make sure that will not be out of range, we should determine it according to
5.11. Velocity in -Frame Calculation
We Have That
5.12. Altitude Calculation
For the calculation of the altitude, damped methods should be used because it diverges with time. To simplify problems, in our simulations, we set the altitude to zero, that is, surface of the Earth.
5.13. Local Gravity g Calculation
The local gravity is calculated as [5] where , is the latitude and is the altitude above sea level.
Before we carry out the implementation of the above described mathematical model of SINS, we have to know the initial parameters of the system, which will be described in the following Section.
6. Initial Parameters and Initial Data Calculation
For the calculations in Section 5, we first need to know the given initial parameters and the corresponding initial data.
6.1. Initial Parameters
(1)Initial position, latitude , longitude , height . The values of these parameters should be the same as the corresponding ones in Section 4.1.(2)Initial wander azimuth angle . We could choose at the very beginning. The value should be the same as the corresponding ones in Section 4.1.(3)Initial velocity , , .(4)If barometric altimeter applied, initial external reference height can be supplied.6.2. Initial Alignment Data
(1)Initial attitude matrix is determined by initial alignment process . when .(2)Initial position matrix is determined by initial alignment process . when .6.3. Initial Data Calculation
(1)Initial attitude angles , , and determination: The initial attitude angles , , and can be calculated using (41) and (42). Because , heading angle .(2)Initial quaternion calculation: From the diagonal elements in (31) and the quaternion constraint equation, we have that The solution to (50) Assuming to be positive, according to (31), we have that (3)Initial position matrix : Substituting initial position, latitude , longitude and initial wander azimuth into (43), we can obtain the initial position matrix .(4)Initial Earth angular velocity and initial attitude angular velocity calculations: use (38) and (39).(5)Initial position angular velocity calculation: use (36) and (37).(6)Initial gravity calculation: use (49) and element in .(7)Initial ground velocity calculation: use (34).At this point, the whole SINS model, including sensor data generator and initial parameters, is fully described. The following Section will provide a Matlab implementation of the SINS theory.
7. Matlab Implementation and Simulation Examples
First, the Matlab program structure and the main codes are given. The Matlab implementation is illustrated using six examples: static, straight, circle, s-shape, and the other two from real SINS experimental data.
7.1. Matlab Implementation and Codes
The program structure is given in Figure 8. The program starts from “Begin” and ends at “Stop.” The gyro and accelerometer data are obtained either from a sensor data generator described in Section 4 or from the real SINS experiment logged files. Processing the sensor data with the initial parameters, using the method described in Section 5, we get the attitude, velocity and position values of the system at specific times. After all data are processed, the program will stop and the results will be provided.
The main Matlab codes are presented next.(1)Initial settings:(a)initSettings.m contains initial parameters and constants used in the simulation project.(2)Trajectory part:(a)initialCalculation_static.m gives the initial calculation for the static situation;(b)trajectorySimulater_static.m simulates gyro and accelerometer data for the static situation;(c)initialCalculation_straight.m gives the initial calculation for the straight line situation;(d)trajectorySimulater_straight.m simulates gyro and accelerometer data for the straight line situation;(e)initialCalculation_cirlce.m gives the initial calculation for the circle situation;(f)trajectorySimulater_circle.m simulates gyro and accelerometer data for the circle situation;(g)initialCalculation_Sshape.m gives the initial calculation for the s-shape situation;(h)trajectorySimulater_Sshape.m simulates gyro and accelerometer data for the s-shape situation;(3)Simulation part:(a)INSmain.m is the main program; the simulation starts from here;(b)AltitudeParamete.m calculates the four damping parameters to damp the altitude error according to the input parameters and , to be used with the external reference altitude;(c)InitializePosition.m gives the initial position initLong, initLat, initAlt, the external reference altitude extAlt, and the wander azimuth angle wanderAzimuth; it calculates the initial position matrix and then orthogonalizes the matrix;(d)InitializeAttitude.m gives the initial alignment error and calculates the attitude matrix (strapdown matrix);(e)InitializeQuaternion.m calculates the quaternion according to the input attitude matrix;(f)ComputeAngularVelocity.m calculates the position angular velocity, earth angular velocity, and position angle increment in the -frame and resets the gyroscopes and accelerometers;(g)ComputeQuaternionRungeKutta.m computes the quaternion using Runge-Kutta method [6]; see Appendices B and C;(h)ComputeAttitudeMatrix.m computes the attitude matrix and transfers the raw data of the accelerometers to the -frame;(i)ComputeVelocity.m computes the velocity, in the wander azimuth frame (-frame) and ENU-frame, the ground velocity and altitude;(j)ComputePositionMatrix.m computes the position matrix.(k)ComputePosition.m computes latitude, longitude and wanderAzimuth;(l)ComputeAttitudeAngle.m computes the attitude angle of pitchAngle, tiltAngle, gridAzimuth and courseAngle;(m)OrthogonalizeMatrix.m computes matrix orthogonalization; see Appendix D;(n)QuaCofMatrix.m is called by ComputeQuaternionRungeKutta.m;(o)PlotResult.m plots the results of the simulation project.
7.2. Simulation Examples
In this subsection, there are 6 SINS simulation examples. Example 1 is the static situation simulation, where the vehicle trajectory in the -frame is a fixed point. Example 2 is the straight line situation simulation, where the vehicle trajectory in the -frame is a straight line. Example 3 is the circle situation simulation, where the vehicle trajectory in the -frame is a circle. Example 4 is the s-shape situation simulation, where the vehicle trajectory in the -frame is an s-shape line. Here, high-accuracy SINS simulation is applied to the four situations. The initial latitude and longitude errors are set to be 1 minute. The simulation time is set to 3600 seconds. The initial positions are dependent on the designed trajectories.
In order to verity the validity of the Matlab codes further, two sets of real static data are used, and we refer to these as Examples 5 and 6. The two sets of real data, set A and set B, are collected from the same SINS in the same place but at different times. The 2 data sets are 24 hours long.
All the errors (the angle error, the velocity error, and the position error) will contribute to the distance error in the INS trajectory calculation. Thus, the distance error is a key index of an INS system. The distance error will increase with time, so it is always associated with a time stamp.
Example 1 (Static situation simulation). The static situation is the most basic and simple situation where the output of the gyro is the Earth rotating angular velocity and the output of the accelerometer is the gravity. Figure 9 shows the designed true trajectory. Figure 10 shows the difference between the calculated angle and the true angle. Figure 11 shows the differences between the calculated PV (position and velocity) and the true PV. The maximum value of the distance error in 1 hour is 3.5 nm (nautical mile).
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(d)
Example 2 (Straight line situation simulation). The straight line situation corresponds to a vehicle moving along the northwest direction. Figure 12 shows the designed true trajectory. Figure 13 shows the difference between the calculated angle and the true angle. Figure 14 shows the differences between the calculated PV and the true PV. The maximum value of the distance error in 1 hour is 3.7 nm.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(d)
Example 3 (Circle situation simulation). The circle situation corresponds to a vehicle moving along a circle. Figure 15 shows the designed true trajectory. Figure 16 shows the difference between the calculated angle and the true angle. Figure 17 shows the difference between the calculated PV and the true PV. The maximum value of the distance error in 1 hour is 3.0 nm.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(d)
Example 4 (S-shape situation simulation). The s-shape situation corresponds to a vehicle moving along an s-shaped line. Figure 18 shows the designed true trajectory of s-shape situation simulation. Figure 19 shows the difference between the calculated angle and the true angle. Figure 20 shows the differences between the calculated PV and the true PV. The maximum value of the distance error in 1 hour is 3.3 nm.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(d)
Example 5 (Real static data set A simulation). First, we process data set A [7]. Figure 21 shows the trajectory for the real data set A; from the figure we can conclude that the system is static. In Figure 22, the red line corresponds to the three attitude angle errors of the real system, while the blue line corresponds to the three attitude angle errors processed by the Matlab code. We can also show that the difference between the red and blue lines is negligible. In Figure 23, the red line corresponds to the position and velocity errors of the real system, while the blue line corresponds to the position and velocity errors processed by the Matlab code. We can also see that the difference between the red and blue lines is negligible and this validates the correctness of the Matlab code. The error described by the red lines (output from the real system) is slightly smaller than that described by the blue lines (simulation). This is due to the fact that the real system is processed in a much higher rate and thus its input is more accurate than the simulated system.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(d)
Example 6 (Real static data set B simulation). Figure 24 shows the trajectory of the real data set B; from the figure we can conclude that the system is static too. In Figure 25, the red line corresponds to the three attitude angle errors of the real system, while the blue line corresponds to the three attitude angle errors obtained by the Matlab code when applied to the real raw sensor data set B. We can also see that the difference between the red and blue lines is negligible. In Figure 26, the red line corresponds to the position and velocity errors of the real system, while the blue line corresponds to the position and velocity errors obtained by the Matlab code when applied to the real raw sensor data set B. We can also see that the difference between the red and blue lines is negligible, and this further validates the correctness of the Matlab code.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(d)
8. Conclusions
In this paper, a mathematical model for the strapdown inertial navigation system (SINS) is built and its Matlab implementation is developed. First, a number of Cartesian coordinate reference frames that relate to SINS are introduced, the basic principle of SINS in the wander azimuth navigation frame (-frame) is explained, and the main equations are described. Second, the important attitude direction cosine matrix and position direction cosine matrix in the -frame are defined in detail. Third, the mathematical model for SINS simulation is described in detail. Fourth, a trajectory simulator model is set up to generate data from three orthogonal gyros and three orthogonal accelerometers. The initial parameters and initial data calculations for the mathematical model are also carried out. Finally, a Matlab implementation of SINS is developed. The proposed simulation method is illustrated with four examples, static, straight line, circle, and s-shape trajectories; details are given under the condition that the pitch angle, roll angle, and altitude are constant during the simulation process. Further, two sets of real experimental data are processed to verify the validity of the Matlab code.
Appendices
A. Symmetric Matrix Basic Operation
For a vector , its skew symmetric matrix is
We can easily show that
B. Fourth-Order Runge-Kutta Method
For numerical analysis, the fourth-order Runge-Kutta method is an important iterative method for the approximation of solutions of ordinary differential equations. Here, in this paper, the fourth-order Runge-Kutta method is adopted to update the quaternion.
The steps for the fourth-order Runge-Kutta method are the following.(1)Calculate slope , the slope at the beginning of the interval, to determine the value of at the point using the Euler method: where is the time step between time and time , .(2)Calculate slope , the slope at the midpoint of the interval, to determine the value of at the point using Euler’s method:(3)Calculate slope , again the slope at the midpoint, to determine the value: (4)Calculate slope , the slope at the end of the interval, with its value determined using : (5)Average the four slopes; greater weights are given to the slopes at the midpoint: (6)Finally, using the average slope , the value of is
C. Angular Velocity Extraction
From Appendix B and (29), we need to provide the attitude angular velocity in a period of to update the quaternion. By inspecting the expression of in (39), we know that the variations of and are slow, while changes quickly. So, only needs to be given in a period of . We know that (we next use to simplify notation) is the output of gyro which gives data in the form of angle increment during the time interval . For first-order angular velocity extraction, we have that In order to provide and , we need to do second-order angular velocity extraction: where is the angle increment from time to and is the angle increment from time to time .
D. Matrix Orthogonalization Method
For the direction cosine matrix , the optimal orthogonalization method is to get which makes the following Euclidian function have the minimum value [8]: The expression for is thus where the superscript means the transpose operator. It is difficult to solve the above equation directly. Instead, we use an iterative method. Assume to be initial matrix, and to be the matrix obtained after iterations. The iteration process is as follows: If at the step, the following function: satisfies (e.g., ), then the iteration procedure can be stopped and is taken to be the final result.
Abbreviations and Symbols
| SINS: | Strapdown inertial navigation system |
| DCM: | Direction cosine matrix |
| : | Center of the Earth |
| : | Center of the vehicle |
| : | 3 orthogonal axes or the 3 components of a Cartesian coordinate |
| : | Body frame |
| : | Inertial frame |
| : | Earth frame |
| : | Navigation frame |
| ENU: | East-North-UP navigation frame, which is identical to the -frame in this paper |
| : | Wander azimuth frame |
| : | Transpose of matrix |
| : | Velocity vector measured in -frame with respect to -frame |
| : | Vehicle attitude DCM used to transform the measured angle in-frame to -frame, with its 9 components |
| : | Transpose of is used to transform the measured vector in -frame to -frame |
| : | Vehicle position DCM used to transform the measured vector in-frame to -frame, with its 9 components |
| : | Transpose of is used to transform the measured vector in -frame to -frame |
| : | Vehicle attitude DCM used to transform the measured angle in-frame to -frame |
| : | Specific force vector measured in -frame |
| : | Specific force vector measured in -frame |
| : | Specific force vector measured in -frame; the output of the 3 accelerometers |
| : | Constant value of the turn rate of the Earth, rad/s |
| : | Turn rate of the Earth measured in -frame |
| : | Turn rate of the -frame with respect to -frame, which is measured in -frame; the output of the 3 gyros |
| : | Transport rate of the -frame with respect to -frame, which is measured in -frame |
| : | Turn rate of the -frame with respect to -frame, which is measured in -frame |
| : | Turn rate of the -frame with respect to -frame, which is measured in -frame |
| : | Turn rate of the -frame with respect to -frame, which is measured in -frame |
| : | Turn rate of the -frame with respect to -frame, which is measured in -frame |
| : | Turn rate of the -frame with respect to -frame, which is measured in -frame |
| : | Skew matrix form of |
| : | Skew matrix form of |
| : | Gravity vector measured in -frame |
| : | Local gravity scalar |
| : | Local gravity scalar at sea level |
| : | Grid azimuth angle of the vehicle in -frame with respect to -frame |
| : | Wander azimuth angle of -frame with respect to -frame |
| : | Heading angle of the vehicle in -frame with respect to -frame |
| : | Grid pitch angle of the vehicle in -frame with respect to -frame or-frame |
| : | Grid roll angle of the vehicle in -frame with respect to -frame or-frame |
| : | Increase of the heading angle |
| : | Increase of the grid pitch angle |
| : | Increase of the grid roll angle |
| : | Longitude of the vehicle |
| : | Latitude of the vehicle |
| : | Altitude of the vehicle above the sea level of the Earth |
| , , : | Initial vehicle position (latitude, longitude, height) |
| : | Time step |
| : | Vehicle acceleration |
| : | Initial vehicle velocity (east, north, up) |
| Vehicle velocity (east, north, up) | |
| : | Vehicle ground velocity |
| : | Vehicle position measured in -frame with respect to -frame |
| : | Major eccentricity of the ellipsoid of the Earth |
| : | Length of the semi-major axis of the Earth |
| : | Meridian radius of curvature of the Earth |
| : | Transverse radius of curvature of the Earth |
| : | Free curvature radiuses |
| : | Turn torsion |
| : | Quaternion |
| : | Four components of the quaternion |
| : | Period of the circle trajectory in simulation |
| : | Period of the s-shape trajectory in simulation |
| : | Amplitude of the s-shape trajectory in simulation |
| PV: | Position and velocity. |