Mathematical Model and Matlab Simulation of Strapdown Inertial Navigation System

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.


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 (p-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 p-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 p-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

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 anticlockwise 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 (b-frame): the b-frame, depicted in Figure 2, is an orthogonal axis set which has its origin at the center of the vehicle, point P, and is aligned with the pitch Px b axis, roll P y b axis, and yaw Pz b axis of the vehicle in which the navigation system is installed.
(2) The inertial frame (i-frame): the i-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 Ox i , Oy i , and Oz i , with Oz i being coincident with the Earth polar axis.
(3) The Earth frame (e-frame): the e-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 Ox e , Oy e , and Oz e .The axis Oz e is the Earth polar axis.The axis Ox e 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 ω ie about the axis Oz i .
(4) The navigation frame (n-frame): the n-frame, depicted in Figure 3, is a local geographic navigation frame which has its origin at the location of the navigation system, point P (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 PE, north PN and the local vertical up PU.When the n-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, ω en , is governed by the motion of the point P with respect to the Earth.This is often referred to as the transport rate.(5) The wander azimuth navigation frame (p-frame): the p-frame, depicted in Figure 3, may be used to avoid the singularities in the computation which occur at the poles of the n-frame.Like the n-frame, it is of a local level but rotates through the wander angle about the local vertical.Here we do not call this frame w-frame (w for wander) for notation clarity since w and ω may look similar when printed.Letter p in p-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 p-frame as the navigation frame for vehicle trajectory calculation, for the following reason.In the local geographic navigation frame mechanization, the n-frame is required to rotate continuously as the system moves over the surface of the Earth in order to keep its P y N axis pointing to the true north.In order to achieve this condition worldwide, the n-frame must rotate at much greater rates about its Pz U 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 Pz U axis rapidly in order to maintain the P y N axis pointing to the pole.The heading direction will abruptly change by 180 • 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 z-component of ω p ep is always set to zero, that is, ω p epz = 0.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 P y N axis and P y p 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 p-frame is described.
Along the same lines as in [3], a navigation equation for a wander azimuth system can be constructed as follows: where C p b is the direction cosine matrix used to transform the measured specific force vector in b-frame into p-frame.This matrix propagates in accordance with the following equation: where Ω b pb is the skew symmetric form of ω b pb , the b-frame angular rate with respect to the p-frame.
Equation ( 1) is integrated to generate estimates of the vehicle speed in the wander azimuth frame, v p e .This is then used to generate the turn rate of the wander frame with respect to the Earth frame, ω p ep .The direction cosine matrix which relates the wander frame to the Earth frame, C p e , may be updated using the following equation: T .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 vehicleposition-related position angles.

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.
The vehicle attitude DCM T p b is then obtained as For the p-frame system, the angle between the grid north y p and the true north y N is the wander azimuth angle α.So the angle between the horizontal projection along y p axis of the vehicle's vertical axis z b and the real north y N is the heading angle ψ.We have that So the direction cosine matrix C n b from b-frame to n-frame is The gimbal angles ψ, θ, and γ and the gimbal rates ψ, θ, and γ are related to the body rate ω b nb , which is the turn rate of the b-frame with respect to n-frame and measured in b-frame as follows: 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 p-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.

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 p-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 p-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.

The Initial Parameters.
For the designed trajectory, the initial parameters are (1) initial position, latitude ϕ 0 , longitude λ 0 , height h 0 ; (2 (3) the designed variation of acceleration a, 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.

The Update of Velocity.
The velocity is updated as where Δt is the time step.

The Update of Position.
The position is updated as latitude: 4.4.The Update of Attitude.The attitude angles are updated as pitch: The attitude rates are updated as The expressions for Δθ, Δγ, Δψ, Δ θ, Δ γ, and Δ ψ depend on the designed trajectory.
The direction cosine matrix C n b can be calculated using matrix expression (10).We have that

Gyro Data Generator. The output of the gyros is
where ω b ib is the simulated actual output, I is the 3 × 3 unit matrix, S g is the 3 × 3 diagonal matrix whose diagonal elements correspond to the 3 gyros' scale factor errors, and ε b 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, ω b nb is equal to zero, whereas, in a moving base it is obtained as T and can be expressed as 4.6.Accelerometer Data Generator.The measurement of the accelerometer is the specific force: where f b is the simulated actual output, I is the 3 × 3 unit matrix.S a is the 3 × 3 diagonal matrix whose diagonal elements correspond to the 3 accelerometers' scale factor errors, η b is the bias considered as the sum of a constant noise and a random white noise η b = η b const + η b random .g = [0 0 g] T , and g = 9.7803 + 0.051799C 2 33 − 0.94114 × 10 −6 h (m/s 2 ), where C 33 is the 9th element of C p e and h is the vehicle altitude.

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.

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 p-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 p-frame is still chosen instead of the simpler ENU-frame to give a navigation illustration in a different frame.

Quaternion Q 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 is vehicle attitude DCM which transforms the measured angle in the b-frame to the pframe, with its 9 components T i j , i, j = 1, 2, 3. (Here we use T i j to distinguish it from the components C i j , i, j = 1, 2, 3 of C p e which is used below.) After obtaining q 0 , q 1 , q 2 , and q 3 using (29), C p b can be calculated as 2 q 1 q 2 − q 0 q 3 2 q 1 q 3 + q 0 q 2 2 q 1 q 2 + q 0 q 3 q 2 0 − q 2 1 + q 2 2 − q 2 3 2 q 2 q 3 − q 0 q 1 2 q 1 q 3 − q 0 q 2 2 q 2 q 3 + q 0 q 1 q 2 0 − q 2 1 − q 2 2 + q 2 3 ⎤ ⎥ ⎦. (31)

Specific Force
The ground speed is the vehicle velocity projection on the horizontal plane:   where where the elements of position matrix C p e can be obtained using (35).Thus, the principal values of ψ G , θ, and γ are

Earth Angular Velocity ω
Considering the defined range of the angles, the expressions of the real values of ψ G , θ, γ, and are (42)

Position Angle Calculation. The relation between position matrix C
p e 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 ϕ ←− ϕ principal , (45)

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 (47) 5.11.Velocity v n e in n-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.

Local Gravity g Calculation.
The local gravity g is calculated as [5] g = 9.7803 + 0.051799C 2 33 − 0.94114 × 10 −6 h m/s 2 , (49 where C 33 = sin ϕ, ϕ is the latitude and h 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.

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.
The values of these parameters should be the same as the corresponding ones in Section 4.1.
(2) Initial wander azimuth angle α 0 .We could choose α 0 = 0 at the very beginning.The value should be the same as the corresponding ones in Section 4.1.
(4) If barometric altimeter applied, initial external reference height h ref0 can be supplied.

Simulation Examples.
In this subsection, there are 6 SINS simulation examples.Example 1 is the static situation simulation, where the vehicle trajectory in the n-frame is a fixed point.Example 2 is the straight line situation simulation, where the vehicle trajectory in the n-frame is a straight line.Example 3 is the circle situation simulation, where the vehicle trajectory in the n-frame is a circle.Example 4 is the s-shape situation simulation, where the vehicle trajectory in the n-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).
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.
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.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.
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.
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.

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 (p-frame) is explained, and the main equations are described.Second, the important attitude direction cosine matrix and position direction cosine matrix in the p-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.

A. Symmetric Matrix Basic Operation
For a vector ω = [ω x ω y ω z ] T , 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 k 1 , the slope at the beginning of the interval, to determine the value of y i+1/2 at the point t i+1/2 using the Euler method: where τ is the time step between time t i and time t i+1 , τ = t i+1 − t i .
(2) Calculate slope k 2 , the slope at the midpoint of the interval, to determine the value of y i+1/2 at the point t i+1/2 using Euler's method: (3) Calculate slope k 3 , again the slope at the midpoint, to determine the y i+1 value: ) Calculate slope k 4 , the slope at the end of the interval, with its y i+1 value determined using k 3 : (5) Average the four slopes; greater weights are given to the slopes at the midpoint:

C. Angular Velocity Extraction
From Appendix B and (29), we need to provide the attitude angular velocity ω b pb in a period of τ/2 to update the quaternion.By inspecting the expression of ω b pb in (39), we know that the variations of ω In order to provide ω(t i )ω(t i+1/2 ) and ω(t i+1 ), we need to do second-order angular velocity extraction: where Δθ i1 is the angle increment from time t i to t i+1/2 and Δθ i2 is the angle increment from time t i+1/2 to time t i+1 .

D. Matrix Orthogonalization Method
For the direction cosine matrix C, the optimal orthogonalization method is to get C which makes the following Euclidian function have the minimum value [8]: The expression for C is thus where the superscript T means the transpose operator.It is difficult to solve the above equation directly.Instead, we use an iterative method.Assume C 0 to be initial matrix, and C n to be the matrix obtained after n iterations.The iteration process is as follows: satisfies f n+1 − f n ≤ (e.g., = 10 −10 ), then the iteration procedure can be stopped and C n+1 is taken to be the final result.

Abbreviations and Symbols
SINS: Strapdown inertial navigation system DCM: Direction cosine matrix O: Center of the Earth P: Center of the vehicle x, y, z: 3 orthogonal axes or the 3 components of a Cartesian coordinate b: Bodyframe i: Inertial frame e: Earth frame n: Navigation frame ENU: East-North-UP navigation frame, which is identical to the n-frame in this paper p: Wander azimuth frame

Figure 1 :Figure 2 :
Figure 1: The schema of the proposed mathematical model and Matlab simulation of SINS.
skew symmetric matrix formed from the elements of the angular rate vector ω p ep ; we could have ω p ep = −ω e pe when the rotation angles are reciprocal.Because the z-component of ω p ep is set to zero,

3. 1 .
Vehicle Attitude DCM C p b .The definition of the rotation sequence from p-frame to b-frame is (see Figure 4) x p y p z p zp, ψG −−−→ x e y e z e y p , θ − −− → x e y e z e y p , γ − −− → x b y b z b , (6) where ψ G is the gird azimuth angle (0-360 • ), θ is the pitch angle (−90 • -90 • ), and γ is the roll angle (−180 • -180 • ).The Modelling and Simulation in Engineering above rotation can be written in the following matrix form:

) 3 . 2 .
Vehicle Position DCM C p e .Position matrix C p e is the DCM from e-frame to p-frame.It has the following rotating sequence (see Figure 5): x e y e z e ze, λ −−→ x e y e z e y p ,90 • −ϕ −−−−−→ x e y e z e z p , 90 •

Figure 5 :
Figure 5: The relation between e-frame and p-frame.

Figure 9 :
Figure 9: The designed trajectory of static simulation.

Figure 10 :
Figure 10: Angle error of static simulation.

Figure 11 :
Figure 11: Position and velocity error of static simulation.

Figure 12 :
Figure 12: The designed trajectory of straight line simulation.

) 5 . 4 .
Velocity v p e Calculation.The velocity v p e update can be obtained by solving the following differential equation: vp e = f p − 2ω

) 5 . 5 .C
Position Matrix C p e Update.The update for the position matrix C p e can be obtained by solving the following differential equation, noticing that ω 11 C 12 C 13

) 5 . 6 .
Position Angular Velocity ω p ep Update.In the chosen wander azimuth navigation frame, we have ω p epz = 0, and

C
11 C 12 C 13

Figure 19 :
Figure 19: Angle error of s-shape simulation.

Figure 21 :Figure 22 :Figure 23 :
Figure 21: The trajectory of real data set A.

Figure 24 :Figure 25 :Figure 26 :
Figure 24: The trajectory of real data set B.
Finally, using the average slope k, the value of y i+1 isy i+1 = y i + τk.(B.6) p ep and ω p ie are slow, while ω b ib changes quickly.So, only ω bib needs to be given in a period of τ/2.We know that ω b ib (we next use ω to simplify notation) is the output of gyro which gives data in the form of angle increment Δθ i during the time interval τ.For first-order angular velocity extraction, we have that

(D. 3 )
If at the n + 1 step, the following function:

A
T : Transpose of matrix A v p e : Velocity vector measured in p-frame with respect to e-frame C p b : Vehicle attitude DCM used to transform the measured angle in b-frame to p-frame, with its 9 components T i j , i, j = 1, 2, 3 C b p : Transpose of C p b is used to transform the measured vector in p-frame to b-frame C p e : Vehicle position DCM used to transform the measured vector in e-frame to p-frame, with its 9 components C i j , i, j = 1, 2, 3 C e p : Transpose of C p e is used to transform the measured vector in p-frame to e-frame C n b : Vehicle attitude DCM used to transform the measured angle in b-frame to n-frame f p : Specific force vector measured in p-frame f n : Specific force vector measured in n-frame f b : Specific force vector measured in b-frame; the output of the 3 accelerometers ω ie : Constant value of the turn rate of the Earth, ω ie = 7.2921151467 × 10 −5 rad/s ω n ie : Turn rate of the Earth measured in n-frame ω b ib : Turn rate of the b-frame with respect to i-frame, which is measured in b-frame; the output of the 3 gyros ω n en : Transport rate of the n-frame with respect to e-frame, which is measured in n-frame ω e ie : Turn rate of the e-frame with respect to i-frame, which is measured in e-frame ω p ep : Turn rate of the p-frame with respect to e-frame, which is measured in p-frame ω e pe : Turn rate of the e-frame with respect to p-frame, which is measured in e-frame ω b pb : Turn rate of the b-frame with respect to p-frame, which is measured in b-frame ω b nb : Turn rate of the b-frame with respect to n-frame, which is measured in b-frame Ω b pb : Skew matrix form of ω b pb Ω p ep : Skew matrix form of ω p ep g p : Gravity vector measured in p-frame g: Local gravity scalar g 0 : Local gravity scalar at sea level