# Mathematical Model and Matlab Simulation of Strapdown Inertial Navigation System

Authors: Wen Zhang; Mounir Ghogho; Baolun Yuan
Journal: Modelling and Simulation in Engineering
 (2012)

---

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


---

## Body

## 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 Figure1. 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 are drawn. Mathematical details are given in Appendices A–D.Figure 1
The schema of the proposed mathematical model and Matlab simulation of SINS.


## 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 Figure2. 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 (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 Pxb axis, roll Pyb axis, and yaw Pzb 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 Oxi, Oyi, and Ozi, with Ozi 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 Oxe, Oye, and Oze. The axis Oze is the Earth polar axis. The axis Oxe 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 Ozi.(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.Figure 2
Theb-frame illustration and the definition of axis rotations.Figure 3
The Reference Frames.In this paper, we choose thep-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 PyN axis pointing to the true north. In order to achieve this condition worldwide, the n-frame must rotate at much greater rates about its PzU 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 PzU axis rapidly in order to maintain the PyN 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 ωepp is always set to zero, that is, ωepzp=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 PyN axis and Pyp 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 thep-frame is described.Along the same lines as in [3], a navigation equation for a wander azimuth system can be constructed as follows:(1)v̇ep=Cbpfb-(2Cepωiee+ωepp)vep+gp,
where Cbp 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:(2)Ċbp=CbpΩpbb,
where Ωpbb is the skew symmetric form of ωpbb, 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, vep. This is then used to generate the turn rate of the wander frame with respect to the Earth frame, ωepp. The direction cosine matrix which relates the wander frame to the Earth frame, Cep, may be updated using the following equation:(3)Ċpe=CpeΩepp,(4)(Ċpe)T=(Ωepp)T(Cpe)T=-Ωepp(Cpe)T,
where the superscript T means matrix transposition.Since(Ċpe)T=Ċep, (Cpe)T=Cep and skew symmetric matrix is (Ωepp)T=-Ωepp (see Appendix A), (4) can be rewritten as(5)Ċep=-ΩeppCep,
where Ωepp is a skew symmetric matrix formed from the elements of the angular rate vector ωepp; we could have ωepp=-ωpee when the rotation angles are reciprocal. Because the z-component of ωepp is set to zero, ωepzp=0, the matrix expression of ωepp is ωepp=[ωepxpωepyp0]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 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 DCMCbp

The definition of the rotation sequence fromp-frame to b-frame is (see Figure 4)(6)xpypzp→zp,ψGxe′ye′ze′→yp′′,θxe′′ye′′ze′′→yp′′,γxbybzb,
where ψG 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:(7)Cpb=C3C2C1=[cos⁡γ0-sin⁡γ010sin⁡γ0cos⁡γ][1000cos⁡θsin⁡θ0-sin⁡θcos⁡θ]×[cos⁡ψGsin⁡ψG0-sin⁡ψGcos⁡ψG0001].Figure 4
The Relation Betweenb-Frame and p-Frame.The vehicle attitude DCMTbp is then obtained as(8)Cbp=(Cpb)T=[cos⁡γcos⁡ψG-sin⁡γsin⁡θsin⁡ψG-cos⁡θsin⁡ψGsin⁡γcos⁡ψG+cos⁡γsin⁡θsin⁡ψGcos⁡γsin⁡ψG+sin⁡γsin⁡θcos⁡ψGcos⁡θcos⁡ψGsin⁡γsin⁡ψG-cos⁡γsin⁡θcos⁡ψG-sin⁡γcos⁡θsin⁡θcos⁡γcos⁡θ].For thep-frame system, the angle between the grid north yp and the true north yN is the wander azimuth angle α. So the angle between the horizontal projection along yp′ axis of the vehicle’s vertical axis zb and the real north yN is the heading angle ψ. We have that(9)ψ=ψG+α.
So the direction cosine matrix Cbn from b-frame to n-frame is(10)Cbn=[cos⁡γcos⁡ψ-sin⁡γsin⁡θsin⁡ψ-cos⁡θsin⁡ψsin⁡γcos⁡ψ+cos⁡γsin⁡θsin⁡ψcos⁡γsin⁡ψ+sin⁡γsin⁡θcos⁡ψcos⁡θcos⁡ψsin⁡γsin⁡ψ-cos⁡γsin⁡θcos⁡ψ-sin⁡γcos⁡θsin⁡θcos⁡γcos⁡θ].The gimbal anglesψ, θ, and γ and the gimbal rates ψ̇, θ̇, and γ̇ are related to the body rate ωnbb, which is the turn rate of the b-frame with respect to n-frame and measured in b-frame as follows:(11)[ωnbxbωnbybωnbzb]=[0γ̇0]+C3[θ̇00]+C3C2[00ψ̇]=[cos⁡γθ̇-sin⁡γcos⁡θψ̇γ̇+sin⁡θψ̇sin⁡γθ̇+cos⁡γcos⁡θψ̇].


### 3.2. Vehicle Position DCMCep

Position matrixCep is the DCM from e-frame to p-frame. It has the following rotating sequence (see Figure 5):(12)xeyeze→ze,λxe′ye′ze′→yp′′,90∘-φxe′′ye′′ze′′→zp′′,90∘xEyNzU→zU,αxpypzp,
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:(13)Cep=[-sin⁡αsin⁡φcos⁡λ-cos⁡αsin⁡λ-sin⁡αsin⁡φsin⁡λ+cos⁡αcos⁡λsin⁡αcos⁡φ-cos⁡αsin⁡φcos⁡λ+sin⁡αsin⁡λ-cos⁡αsin⁡φsin⁡λ-sin⁡αcos⁡λcos⁡αcos⁡ψcos⁡φcos⁡λcos⁡φsin⁡λsin⁡φ].Figure 5
The Relation Betweene-Frame and p-Frame.In Section4, 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.


## 3.1. Vehicle Attitude DCMCbp

The definition of the rotation sequence fromp-frame to b-frame is (see Figure 4)(6)xpypzp→zp,ψGxe′ye′ze′→yp′′,θxe′′ye′′ze′′→yp′′,γxbybzb,
where ψG 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:(7)Cpb=C3C2C1=[cos⁡γ0-sin⁡γ010sin⁡γ0cos⁡γ][1000cos⁡θsin⁡θ0-sin⁡θcos⁡θ]×[cos⁡ψGsin⁡ψG0-sin⁡ψGcos⁡ψG0001].Figure 4
The Relation Betweenb-Frame and p-Frame.The vehicle attitude DCMTbp is then obtained as(8)Cbp=(Cpb)T=[cos⁡γcos⁡ψG-sin⁡γsin⁡θsin⁡ψG-cos⁡θsin⁡ψGsin⁡γcos⁡ψG+cos⁡γsin⁡θsin⁡ψGcos⁡γsin⁡ψG+sin⁡γsin⁡θcos⁡ψGcos⁡θcos⁡ψGsin⁡γsin⁡ψG-cos⁡γsin⁡θcos⁡ψG-sin⁡γcos⁡θsin⁡θcos⁡γcos⁡θ].For thep-frame system, the angle between the grid north yp and the true north yN is the wander azimuth angle α. So the angle between the horizontal projection along yp′ axis of the vehicle’s vertical axis zb and the real north yN is the heading angle ψ. We have that(9)ψ=ψG+α.
So the direction cosine matrix Cbn from b-frame to n-frame is(10)Cbn=[cos⁡γcos⁡ψ-sin⁡γsin⁡θsin⁡ψ-cos⁡θsin⁡ψsin⁡γcos⁡ψ+cos⁡γsin⁡θsin⁡ψcos⁡γsin⁡ψ+sin⁡γsin⁡θcos⁡ψcos⁡θcos⁡ψsin⁡γsin⁡ψ-cos⁡γsin⁡θcos⁡ψ-sin⁡γcos⁡θsin⁡θcos⁡γcos⁡θ].The gimbal anglesψ, θ, and γ and the gimbal rates ψ̇, θ̇, and γ̇ are related to the body rate ωnbb, which is the turn rate of the b-frame with respect to n-frame and measured in b-frame as follows:(11)[ωnbxbωnbybωnbzb]=[0γ̇0]+C3[θ̇00]+C3C2[00ψ̇]=[cos⁡γθ̇-sin⁡γcos⁡θψ̇γ̇+sin⁡θψ̇sin⁡γθ̇+cos⁡γcos⁡θψ̇].


## 3.2. Vehicle Position DCMCep

Position matrixCep is the DCM from e-frame to p-frame. It has the following rotating sequence (see Figure 5):(12)xeyeze→ze,λxe′ye′ze′→yp′′,90∘-φxe′′ye′′ze′′→zp′′,90∘xEyNzU→zU,αxpypzp,
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:(13)Cep=[-sin⁡αsin⁡φcos⁡λ-cos⁡αsin⁡λ-sin⁡αsin⁡φsin⁡λ+cos⁡αcos⁡λsin⁡αcos⁡φ-cos⁡αsin⁡φcos⁡λ+sin⁡αsin⁡λ-cos⁡αsin⁡φsin⁡λ-sin⁡αcos⁡λcos⁡αcos⁡ψcos⁡φcos⁡λcos⁡φsin⁡λsin⁡φ].Figure 5
The Relation Betweene-Frame and p-Frame.In Section4, 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.


## 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 Section2 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.Figure 6
SINSENU-Frame Mechanization.


### 4.1. The Initial Parameters

For the designed trajectory, the initial parameters are(1)
initial position, latitudeφ0, longitude λ0, height h0;(2)
initial velocityv=[vE0,vN0,vU0];(3)
the designed variation of accelerationa, 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(14)v⟵v+aΔt,whereΔt is the time step.


### 4.3. The Update of Position

The position is updated as(15)latitude:L⟵L+vNΔtRN,longitude:λ⟵λ+vEΔtsecLREaltitude:h⟵h+vUΔt.,


### 4.4. The Update of Attitude

The attitude angles are updated as(16)pitch:θ⟵θ+Δθ,roll:γ⟵γ+Δγ,heading:ψ⟵ψ+Δψ.The attitude rates are updated as(17)pitch:θ̇⟵θ̇+Δθ̇,roll:γ̇⟵γ̇+Δγ̇,heading:ψ̇⟵ψ̇+Δψ̇.The expressions forΔθ, Δγ, Δψ, Δθ̇, Δγ̇, and Δψ̇ depend on the designed trajectory.The direction cosine matrixCbn can be calculated using matrix expression (10). We have that Cnb=(Cbn)T.


### 4.5. Gyro Data Generator

The output of the gyros is(18)ωibb=(I-Sg)(Cnb(ωien+ωenn)+ωnbb)+ɛb,whereωibb is the simulated actual output, I is the 3×3 unit matrix, Sg 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: ɛb=ɛconstb+ɛrandomb:(19)ωien=[0ωiecos⁡Lωiesin⁡L].In a static base,ωnbb is equal to zero, whereas, in a moving base it is obtained as(20)ωnbb=[cos⁡γθ̇-sin⁡γcos⁡θψ̇γ̇+sin⁡θψ̇sin⁡γθ̇+cos⁡γcos⁡θψ̇].ωenn is related to velocity v=[vE,vN,vU]T and can be expressed as(21)ωenn=[-vNRNvEREvEtanLRN].


### 4.6. Accelerometer Data Generator

The measurement of the accelerometer is the specific force:(22)fb=(I-Sa)Cnbfn+ηb,fn=a+(2ωien+ωenn)×v-g,wherefb is the simulated actual output, I is the 3×3 unit matrix. Sa 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=ηconstb+ηrandomb. g=[00g]T, and g=9.7803+0.051799C332-0.94114×10-6h(m/s2), where C33 is the 9th element of Cep and h 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:(23)Δθ=0,Δγ=0,Δθ̇=0,Δγ̇=0,The calculation method for the other parameters for the four situations is described as follows.(1)
Static:(24)L=constant,λ=constant,vE=constant,vN=constant,Δψ=0,Δψ̇=0.(2)
Straight line:(25)aE=constant,aN=constant,vE=vE+aEΔt,vN=vN+aNΔt,ψ=ψ0+arctan(vEvN),ψ̇=aNvE-aEvNvE2+vN2.(3)
Circle:(26)vg=constant,Δψ=mod⁡(2πΔtTcircle,2π),Δψ̇=2πTcircle,aE=-2πvgcos⁡ψTcircle,aN=-2πvgsin⁡ψTcircle.(4)
S-shape:(27)vg=constant,aE=-vgcos⁡(ψ0+Asshapesin⁡(2πt/Tsshape))Tsshape,⋅2πAsshapecos⁡(2πt/Tsshape)Tsshape,aN=-vgsin⁡(ψ0+Asshapesin⁡(2πt/Tsshape))Tsshape⋅2πAsshapecos⁡(2πt/Tsshape)Tsshape,ψ=ψ0+Asshapesin⁡(2πtTsshape).ψ̇=2πAsshapecos⁡(2πt/Tsshape)Tsshape.


## 4.1. The Initial Parameters

For the designed trajectory, the initial parameters are(1)
initial position, latitudeφ0, longitude λ0, height h0;(2)
initial velocityv=[vE0,vN0,vU0];(3)
the designed variation of accelerationa, 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(14)v⟵v+aΔt,whereΔt is the time step.


## 4.3. The Update of Position

The position is updated as(15)latitude:L⟵L+vNΔtRN,longitude:λ⟵λ+vEΔtsecLREaltitude:h⟵h+vUΔt.,


## 4.4. The Update of Attitude

The attitude angles are updated as(16)pitch:θ⟵θ+Δθ,roll:γ⟵γ+Δγ,heading:ψ⟵ψ+Δψ.The attitude rates are updated as(17)pitch:θ̇⟵θ̇+Δθ̇,roll:γ̇⟵γ̇+Δγ̇,heading:ψ̇⟵ψ̇+Δψ̇.The expressions forΔθ, Δγ, Δψ, Δθ̇, Δγ̇, and Δψ̇ depend on the designed trajectory.The direction cosine matrixCbn can be calculated using matrix expression (10). We have that Cnb=(Cbn)T.


## 4.5. Gyro Data Generator

The output of the gyros is(18)ωibb=(I-Sg)(Cnb(ωien+ωenn)+ωnbb)+ɛb,whereωibb is the simulated actual output, I is the 3×3 unit matrix, Sg 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: ɛb=ɛconstb+ɛrandomb:(19)ωien=[0ωiecos⁡Lωiesin⁡L].In a static base,ωnbb is equal to zero, whereas, in a moving base it is obtained as(20)ωnbb=[cos⁡γθ̇-sin⁡γcos⁡θψ̇γ̇+sin⁡θψ̇sin⁡γθ̇+cos⁡γcos⁡θψ̇].ωenn is related to velocity v=[vE,vN,vU]T and can be expressed as(21)ωenn=[-vNRNvEREvEtanLRN].


## 4.6. Accelerometer Data Generator

The measurement of the accelerometer is the specific force:(22)fb=(I-Sa)Cnbfn+ηb,fn=a+(2ωien+ωenn)×v-g,wherefb is the simulated actual output, I is the 3×3 unit matrix. Sa 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=ηconstb+ηrandomb. g=[00g]T, and g=9.7803+0.051799C332-0.94114×10-6h(m/s2), where C33 is the 9th element of Cep and h 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:(23)Δθ=0,Δγ=0,Δθ̇=0,Δγ̇=0,The calculation method for the other parameters for the four situations is described as follows.(1)
Static:(24)L=constant,λ=constant,vE=constant,vN=constant,Δψ=0,Δψ̇=0.(2)
Straight line:(25)aE=constant,aN=constant,vE=vE+aEΔt,vN=vN+aNΔt,ψ=ψ0+arctan(vEvN),ψ̇=aNvE-aEvNvE2+vN2.(3)
Circle:(26)vg=constant,Δψ=mod⁡(2πΔtTcircle,2π),Δψ̇=2πTcircle,aE=-2πvgcos⁡ψTcircle,aN=-2πvgsin⁡ψTcircle.(4)
S-shape:(27)vg=constant,aE=-vgcos⁡(ψ0+Asshapesin⁡(2πt/Tsshape))Tsshape,⋅2πAsshapecos⁡(2πt/Tsshape)Tsshape,aN=-vgsin⁡(ψ0+Asshapesin⁡(2πt/Tsshape))Tsshape⋅2πAsshapecos⁡(2πt/Tsshape)Tsshape,ψ=ψ0+Asshapesin⁡(2πtTsshape).ψ̇=2πAsshapecos⁡(2πt/Tsshape)Tsshape.


## 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 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.Figure 7
SINSp-Frame Mechanization.


### 5.1. QuaternionQ 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(28)Q=q0+q1ib+q2jb+q3kb.
The update for the quaternion can be obtained by solving the following quaternion differential equation:(29)[q̇0q̇1q̇2q̇3]=12[0-ωpbxb-ωpbyb-ωpbzbωpbxb0ωpbzb-ωpbybωpbyb-ωpbzb0ωpbxbωpbzbωpbyb-ωpbxb0][q0q1q2q3].
Based on the Euclide norm minimized indicator [4], the optimal normalization for the quaternion is(30)Q⟵Qq02+q12+q22+q32.


### 5.2.Cbp Calculation

Cbp is vehicle attitude DCM which transforms the measured angle in the b-frame to the p-frame, with its 9 components Tij,i,j=1,2,3. (Here we use Tij to distinguish it from the components Cij,i,j=1,2,3 of Cep which is used below.)After obtainingq0, q1, q2, and q3 using (29), Cbp can be calculated as(31)Cbp=[T11T12T13T21T22T23T31T32T33]=[q02+q12-q22-q322(q1q2-q0q3)2(q1q3+q0q2)2(q1q2+q0q3)q02-q12+q22-q322(q2q3-q0q1)2(q1q3-q0q2)2(q2q3+q0q1)q02-q12-q22+q32].


### 5.3. Specific Force Transformation fromfb in b-Frame to fp in p-Frame

The specific forcefb in the b-frame can be transformed to fp in the p-frame by multiplication with DCM Cbp:(32)fp=Cbpfb,[fxpfypfzp]=Cbp[fxbfybfzb],


### 5.4. Velocityvep Calculation

The velocityvep update can be obtained by solving the following differential equation:(33)v̇ep=fp-(2ωiep+ωepp)vep+gp,[v̇xv̇yv̇z]=[fxpfypfzp]-[00g]+[02ωiezp-(2ωieyp+ωepyp)-2ωiezp02ωiexp+ωepxp2ωieyp+ωepyp-(2ωiexp+ωepxp)0]×[vxvyvz].
The ground speed is the vehicle velocity projection on the horizontal plane:(34)vg=vx2+vy2.


### 5.5. Position MatrixCep Update

The update for the position matrixCep can be obtained by solving the following differential equation, noticing that ωepzp=0:(35)Ċep=-ΩeppCep,Cep=[C11C12C13C21C22C23C31C32C33],[Ċ11Ċ12Ċ13Ċ21Ċ22Ċ23Ċ31Ċ32Ċ33]=[00-ωepyp00ωepxpωepyp-ωepxp0][C11C12C13C21C22C23C31C32C33].


### 5.6. Position Angular Velocityωepp Update

In the chosen wander azimuth navigation frame, we haveωepzp=0, and(36)[ωepxpωepyp]=[-1τa-1Ryp1Rxp1τa][vexpveyp],
where(37)1Ryp=1Re(1-eC332+2eC232),1Rxp=1Re(1-eC332+2eC132),1τa=2eReC13C23,
where the elements of position matrix Cep can be obtained using (35).


### 5.7. Earth Angular Velocityωiep and Attitude Angular Velocity ωpbb Calculation

We Have That(38)ωiep=Cepωiee=[C11C12C13C21C22C23C31C32C33][00ωie]=[ωieC13ωieC23ωieC33],(39)ωpbb=ωibb-ωipb=ωibb-(Cbp)-1(ωiep+ωepp),


### 5.8. Attitude Angle Calculation

The relation between attitude matrixCbp and the three attitude angles, grid azimuth angle ψG, pitch angle θ, and roll angle γ, is(40)Cbp=[cos⁡γcos⁡ψG-sin⁡γsin⁡θsin⁡ψG-cos⁡θsin⁡ψGsin⁡γcos⁡ψG+cos⁡γsin⁡θsin⁡ψGcos⁡γsin⁡ψG+sin⁡γsin⁡θcos⁡ψGcos⁡θcos⁡ψGsin⁡γsin⁡ψG-cos⁡γsin⁡θcos⁡ψG-sin⁡γcos⁡θsin⁡θcos⁡γcos⁡θ].Thus, the principal values ofψG, θ, and γ are(41)θprincipal=sin⁡-1T32,γprincipal=tan-1-T31T33,φGprincipal=tan-1-T12T22.
Considering the defined range of the angles, the expressions of the real values of ψG, θ, γ, and are(42)θ⟵θprincipal,γ⟵{γprincipal,ifT33>0,γprincipal+180∘,ifT33<0,γprincipal<0,γprincipal-180∘,ifT33<0,γprincipal>0,ψG⟵{ψGprincipal,ifT22>0,ψGprincipal>0,ψGprincipal+360∘,ifT22>0,ψGprincipal<0,ψGprincipal+180∘,ifT22<0.


### 5.9. Position Angle Calculation

The relation between position matrixCep and the 3 position angles, longitude λ, latitude φ, and wander azimuth angle α, is(43)Cep=[-sin⁡αsin⁡φcos⁡λ-cos⁡αsin⁡λ-sin⁡αsin⁡φsin⁡λ+cos⁡αcos⁡λsin⁡αcos⁡φ-cos⁡αsin⁡φcos⁡λ+sin⁡αsin⁡λ-cos⁡αsin⁡φsin⁡λ-sin⁡αcos⁡λcos⁡αcos⁡ψcos⁡φcos⁡λcos⁡φsin⁡λsin⁡φ].Thus, the principal values ofφ, λ, and α are(44)φprincipal=sin⁡-1C33,λprincipal=tan-1C32C31,αprincipal=tan-1C13C23.
Considering the defined range of the angles, the expressions of the real values of φ, λ, and α are(45)φ⟵φprincipal,λ⟵{λprincipal,ifC31>0,λprincipal+180∘,ifC31<0,λprincipal<0,λprincipal-180∘,ifC31<0,λprincipal>0,α⟵{αprincipal,ifC23>0,αprincipal>0,αprincipal+360∘,ifC23>0,αprincipal<0,αprincipal+180∘,ifC23<0.


### 5.10. Heading Angle Calculation

The heading angleψ is calculated as(46)ψ=ψG+α.
To make sure that ψ will not be out of range, we should determine it according to(47)ψ⟵{ψ,ifψ<360∘,ψ-360∘,ifψ≥360∘.


### 5.11. Velocityven in n-Frame Calculation

We Have That(48)ven=[vEvNvU]=[veypcos⁡α-vexpsin⁡αveypsin⁡α+vexpcos⁡αvezp].


### 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 gravityg is calculated as [5](49)g=9.7803+0.051799C332-0.94114×10-6h(m/s2),
where C33=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.


## 5.1. QuaternionQ 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(28)Q=q0+q1ib+q2jb+q3kb.
The update for the quaternion can be obtained by solving the following quaternion differential equation:(29)[q̇0q̇1q̇2q̇3]=12[0-ωpbxb-ωpbyb-ωpbzbωpbxb0ωpbzb-ωpbybωpbyb-ωpbzb0ωpbxbωpbzbωpbyb-ωpbxb0][q0q1q2q3].
Based on the Euclide norm minimized indicator [4], the optimal normalization for the quaternion is(30)Q⟵Qq02+q12+q22+q32.


## 5.2.Cbp Calculation

Cbp is vehicle attitude DCM which transforms the measured angle in the b-frame to the p-frame, with its 9 components Tij,i,j=1,2,3. (Here we use Tij to distinguish it from the components Cij,i,j=1,2,3 of Cep which is used below.)After obtainingq0, q1, q2, and q3 using (29), Cbp can be calculated as(31)Cbp=[T11T12T13T21T22T23T31T32T33]=[q02+q12-q22-q322(q1q2-q0q3)2(q1q3+q0q2)2(q1q2+q0q3)q02-q12+q22-q322(q2q3-q0q1)2(q1q3-q0q2)2(q2q3+q0q1)q02-q12-q22+q32].


## 5.3. Specific Force Transformation fromfb in b-Frame to fp in p-Frame

The specific forcefb in the b-frame can be transformed to fp in the p-frame by multiplication with DCM Cbp:(32)fp=Cbpfb,[fxpfypfzp]=Cbp[fxbfybfzb],


## 5.4. Velocityvep Calculation

The velocityvep update can be obtained by solving the following differential equation:(33)v̇ep=fp-(2ωiep+ωepp)vep+gp,[v̇xv̇yv̇z]=[fxpfypfzp]-[00g]+[02ωiezp-(2ωieyp+ωepyp)-2ωiezp02ωiexp+ωepxp2ωieyp+ωepyp-(2ωiexp+ωepxp)0]×[vxvyvz].
The ground speed is the vehicle velocity projection on the horizontal plane:(34)vg=vx2+vy2.


## 5.5. Position MatrixCep Update

The update for the position matrixCep can be obtained by solving the following differential equation, noticing that ωepzp=0:(35)Ċep=-ΩeppCep,Cep=[C11C12C13C21C22C23C31C32C33],[Ċ11Ċ12Ċ13Ċ21Ċ22Ċ23Ċ31Ċ32Ċ33]=[00-ωepyp00ωepxpωepyp-ωepxp0][C11C12C13C21C22C23C31C32C33].


## 5.6. Position Angular Velocityωepp Update

In the chosen wander azimuth navigation frame, we haveωepzp=0, and(36)[ωepxpωepyp]=[-1τa-1Ryp1Rxp1τa][vexpveyp],
where(37)1Ryp=1Re(1-eC332+2eC232),1Rxp=1Re(1-eC332+2eC132),1τa=2eReC13C23,
where the elements of position matrix Cep can be obtained using (35).


## 5.7. Earth Angular Velocityωiep and Attitude Angular Velocity ωpbb Calculation

We Have That(38)ωiep=Cepωiee=[C11C12C13C21C22C23C31C32C33][00ωie]=[ωieC13ωieC23ωieC33],(39)ωpbb=ωibb-ωipb=ωibb-(Cbp)-1(ωiep+ωepp),


## 5.8. Attitude Angle Calculation

The relation between attitude matrixCbp and the three attitude angles, grid azimuth angle ψG, pitch angle θ, and roll angle γ, is(40)Cbp=[cos⁡γcos⁡ψG-sin⁡γsin⁡θsin⁡ψG-cos⁡θsin⁡ψGsin⁡γcos⁡ψG+cos⁡γsin⁡θsin⁡ψGcos⁡γsin⁡ψG+sin⁡γsin⁡θcos⁡ψGcos⁡θcos⁡ψGsin⁡γsin⁡ψG-cos⁡γsin⁡θcos⁡ψG-sin⁡γcos⁡θsin⁡θcos⁡γcos⁡θ].Thus, the principal values ofψG, θ, and γ are(41)θprincipal=sin⁡-1T32,γprincipal=tan-1-T31T33,φGprincipal=tan-1-T12T22.
Considering the defined range of the angles, the expressions of the real values of ψG, θ, γ, and are(42)θ⟵θprincipal,γ⟵{γprincipal,ifT33>0,γprincipal+180∘,ifT33<0,γprincipal<0,γprincipal-180∘,ifT33<0,γprincipal>0,ψG⟵{ψGprincipal,ifT22>0,ψGprincipal>0,ψGprincipal+360∘,ifT22>0,ψGprincipal<0,ψGprincipal+180∘,ifT22<0.


## 5.9. Position Angle Calculation

The relation between position matrixCep and the 3 position angles, longitude λ, latitude φ, and wander azimuth angle α, is(43)Cep=[-sin⁡αsin⁡φcos⁡λ-cos⁡αsin⁡λ-sin⁡αsin⁡φsin⁡λ+cos⁡αcos⁡λsin⁡αcos⁡φ-cos⁡αsin⁡φcos⁡λ+sin⁡αsin⁡λ-cos⁡αsin⁡φsin⁡λ-sin⁡αcos⁡λcos⁡αcos⁡ψcos⁡φcos⁡λcos⁡φsin⁡λsin⁡φ].Thus, the principal values ofφ, λ, and α are(44)φprincipal=sin⁡-1C33,λprincipal=tan-1C32C31,αprincipal=tan-1C13C23.
Considering the defined range of the angles, the expressions of the real values of φ, λ, and α are(45)φ⟵φprincipal,λ⟵{λprincipal,ifC31>0,λprincipal+180∘,ifC31<0,λprincipal<0,λprincipal-180∘,ifC31<0,λprincipal>0,α⟵{αprincipal,ifC23>0,αprincipal>0,αprincipal+360∘,ifC23>0,αprincipal<0,αprincipal+180∘,ifC23<0.


## 5.10. Heading Angle Calculation

The heading angleψ is calculated as(46)ψ=ψG+α.
To make sure that ψ will not be out of range, we should determine it according to(47)ψ⟵{ψ,ifψ<360∘,ψ-360∘,ifψ≥360∘.


## 5.11. Velocityven in n-Frame Calculation

We Have That(48)ven=[vEvNvU]=[veypcos⁡α-vexpsin⁡αveypsin⁡α+vexpcos⁡αvezp].


## 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 gravityg is calculated as [5](49)g=9.7803+0.051799C332-0.94114×10-6h(m/s2),
where C33=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.


## 6. Initial Parameters and Initial Data Calculation

For the calculations in Section5, we first need to know the given initial parameters and the corresponding initial data.


### 6.1. Initial Parameters

(1)
Initial position, latitudeφ0, longitude λ0, height h0. 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.(3)
Initial velocityvE0, vN0, vU0.(4)
If barometric altimeter applied, initial external reference heighthref0 can be supplied.


### 6.2. Initial Alignment Data

(1)
Initial attitude matrix is determined by initial alignment processCb0p. Cb0p=Cb0n when α0=0.(2)
Initial position matrix is determined by initial alignment processCe0p. Ce0p=Ce0n when α0=0.


### 6.3. Initial Data Calculation

(1)
Initial attitude anglesφ0, λ0, and α0 determination: The initial attitude angles ψG0, θ0, and γ0 can be calculated using (41) and (42). Because α0=0, heading angle ψ0=ψG0.(2)
Initial quaternion calculation: From the diagonal elements in (31) and the quaternion constraint equation, we have that
(50)q02+q12-q22-q32=T11,q02-q12+q22-q32=T22,q02-q12-q22+q32=T33,q02+q12+q22+q32=1,
The solution to (50)
(51)|q1|=121+T11-T22-T33,|q2|=121-T11+T22-T33,|q3|=121-T11-T22+T33,|q0|=1-q12-q22-q32.
Assuming q0 to be positive, according to (31), we have that
(52)sign⁡(q0)=sign⁡(1),sign⁡(q1)=sign⁡(T32-T23),sign⁡(q2)=sign⁡(T13-T31),sign⁡(q3)=sign⁡(T21-T12).(3)
Initial position matrixCe0p: Substituting initial position, latitude φ0, longitude λ0 and initial wander azimuth α0=0 into (43), we can obtain the initial position matrix Ce0p.(4)
Initial Earth angular velocityωie0p and initial attitude angular velocity ωpb0b calculations: use (38) and (39).(5)
Initial position angular velocityωep0p calculation: use (36) and (37).(6)
Initial gravityg0 calculation: use (49) and element C33 in Ce0p.(7)
Initial ground velocityvg0 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.


## 6.1. Initial Parameters

(1)
Initial position, latitudeφ0, longitude λ0, height h0. 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.(3)
Initial velocityvE0, vN0, vU0.(4)
If barometric altimeter applied, initial external reference heighthref0 can be supplied.


## 6.2. Initial Alignment Data

(1)
Initial attitude matrix is determined by initial alignment processCb0p. Cb0p=Cb0n when α0=0.(2)
Initial position matrix is determined by initial alignment processCe0p. Ce0p=Ce0n when α0=0.


## 6.3. Initial Data Calculation

(1)
Initial attitude anglesφ0, λ0, and α0 determination: The initial attitude angles ψG0, θ0, and γ0 can be calculated using (41) and (42). Because α0=0, heading angle ψ0=ψG0.(2)
Initial quaternion calculation: From the diagonal elements in (31) and the quaternion constraint equation, we have that
(50)q02+q12-q22-q32=T11,q02-q12+q22-q32=T22,q02-q12-q22+q32=T33,q02+q12+q22+q32=1,
The solution to (50)
(51)|q1|=121+T11-T22-T33,|q2|=121-T11+T22-T33,|q3|=121-T11-T22+T33,|q0|=1-q12-q22-q32.
Assuming q0 to be positive, according to (31), we have that
(52)sign⁡(q0)=sign⁡(1),sign⁡(q1)=sign⁡(T32-T23),sign⁡(q2)=sign⁡(T13-T31),sign⁡(q3)=sign⁡(T21-T12).(3)
Initial position matrixCe0p: Substituting initial position, latitude φ0, longitude λ0 and initial wander azimuth α0=0 into (43), we can obtain the initial position matrix Ce0p.(4)
Initial Earth angular velocityωie0p and initial attitude angular velocity ωpb0b calculations: use (38) and (39).(5)
Initial position angular velocityωep0p calculation: use (36) and (37).(6)
Initial gravityg0 calculation: use (49) and element C33 in Ce0p.(7)
Initial ground velocityvg0 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 Figure8. 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.Figure 8
SINS Program Structure.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 k4 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 p-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 p-frame;(i)
ComputeVelocity.m computes the velocity, in the wander azimuth frame (p-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. Example1 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 Examples5 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. Figure9 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).The designed trajectory of static simulation.
(a)(b)(c)Angle error of static simulation.
(a)(b)(c)Position and velocity error of static simulation.
(a)(b)(c)(d)Example 2 (Straight line situation simulation).
The straight line situation corresponds to a vehicle moving along the northwest direction. Figure12 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.The designed trajectory of straight line simulation.
(a)(b)(c)Angle error of straight line simulation.
(a)(b)(c)Position and velocity error of straight line simulation.
(a)(b)(c)(d)Example 3 (Circle situation simulation).
The circle situation corresponds to a vehicle moving along a circle. Figure15 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.The designed trajectory of circle simulation.
(a)(b)(c)Angle error of circle simulation.
(a)(b)(c)Position and velocity error of circle simulation.
(a)(b)(c)(d)Example 4 (S-shape situation simulation).
The s-shape situation corresponds to a vehicle moving along an s-shaped line. Figure18 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.The designed trajectory of s-shape simulation.
(a)(b)(c)Angle error of s-shape simulation.
(a)(b)(c)Position and velocity error of s-shape simulation.
(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.The trajectory of real data set A.
(a)(b)(c)Angle error of real data set A.
(a)(b)(c)Position and velocity error of real data set A.
(a)(b)(c)(d)Example 6 (Real static data set B simulation).
Figure24 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.The trajectory of real data set B.
(a)(b)(c)Angle error of real data set B.
(a)(b)(c)Position and Velocity Error of Real Data Set B.
(a)(b)(c)(d)


## 7.1. Matlab Implementation and Codes

The program structure is given in Figure8. 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.Figure 8
SINS Program Structure.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 k4 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 p-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 p-frame;(i)
ComputeVelocity.m computes the velocity, in the wander azimuth frame (p-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. Example1 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 Examples5 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. Figure9 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).The designed trajectory of static simulation.
(a)(b)(c)Angle error of static simulation.
(a)(b)(c)Position and velocity error of static simulation.
(a)(b)(c)(d)Example 2 (Straight line situation simulation).
The straight line situation corresponds to a vehicle moving along the northwest direction. Figure12 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.The designed trajectory of straight line simulation.
(a)(b)(c)Angle error of straight line simulation.
(a)(b)(c)Position and velocity error of straight line simulation.
(a)(b)(c)(d)Example 3 (Circle situation simulation).
The circle situation corresponds to a vehicle moving along a circle. Figure15 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.The designed trajectory of circle simulation.
(a)(b)(c)Angle error of circle simulation.
(a)(b)(c)Position and velocity error of circle simulation.
(a)(b)(c)(d)Example 4 (S-shape situation simulation).
The s-shape situation corresponds to a vehicle moving along an s-shaped line. Figure18 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.The designed trajectory of s-shape simulation.
(a)(b)(c)Angle error of s-shape simulation.
(a)(b)(c)Position and velocity error of s-shape simulation.
(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.The trajectory of real data set A.
(a)(b)(c)Angle error of real data set A.
(a)(b)(c)Position and velocity error of real data set A.
(a)(b)(c)(d)Example 6 (Real static data set B simulation).
Figure24 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.The trajectory of real data set B.
(a)(b)(c)Angle error of real data set B.
(a)(b)(c)Position and Velocity Error of Real Data Set B.
(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 (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.



---

*Source: 264537-2012-01-05.xml*