September 2015

NASA/TM–2015-218803

A Short Tutorial on Inertial Navigation System 

and Global Positioning System Integration

Kyle M. Smalling 

Northrop Grumman, Hampton, Virginia

Kenneth W. Eure

Langley Research Center, Hampton, Virginia
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the auspices 

of the Agency Chief Information Officer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NTRS Registered and its public interface, the 

NASA Technical Reports Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in both 

non-NASA channels and by NASA in the NASA STI 

Report Series, which includes the following report 

types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase of 

research that present the results of NASA 

Programs and include extensive data or theoretical 

analysis. Includes compilations of significant 

scientific and technical data and information 

deemed to be of continuing reference value. 

NASA counter-part of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations. 

• TECHNICAL MEMORANDUM. 

Scientific and technical findings that are 

preliminary or of specialized interest, 

e.g., quick release reports, working 

papers, and bibliographies that contain minimal 

annotation. Does not contain extensive analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. 

Collected papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or 

co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from NASA 

programs, projects, and missions, often 

concerned with subjects having substantial 

public interest. 

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI program, 

see the following: 

• Access the NASA STI program home page at 

http://www.sti.nasa.gov 

• E-mail your question to help@sti.nasa.gov 

• Phone the NASA STI Information Desk at 

757-864-9658 

• Write to: 

NASA STI Information Desk 

Mail Stop 148 

NASA Langley Research Center 

Hampton, VA 23681-2199 
September 2015

NASA/TM–2015-218803

A Short Tutorial on Inertial Navigation System and 

Global Positioning System Integration

Kyle M. Smalling

Northrop Grumman, Hampton, Virginia

Kenneth W. Eure

Langley Research Center, Hampton, Virginia

National Aeronautics and 

Space Administration

Langley Research Center 

Hampton, Virginia 23681-2199
Available from: 

NASA STI Program / Mail Stop 148 

NASA Langley Research Center 

Hampton, VA 23681-2199 

Fax: 757-864-6500 

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 

official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics and 

Space Administration. 
A Short Tutorial on Inertial Navigation System and Global Positioning 

System Integration 

Kyle Smalling, Northrop Grumman

Kenneth W. Eure, NASA Langley Research Center

Abstract

The purpose of this document is to describe a simple method of integrating Inertial Navigation

System (INS) information with Global Positioning System (GPS) information for an improved estimate of

vehicle attitude and position. A simple two dimensional (2D) case is considered. The attitude estimates

are derived from sensor data and used in the estimation of vehicle position and velocity through dead

reckoning within the INS. The INS estimates are updated with GPS estimates using a Kalman filter. This

tutorial isintended for the novice user with a focus on bringing the reader from raw sensor measurements

to an integrated position and attitude estimate. An application is given using a remotely controlled ground

vehicle operating in assumed 2D environment. The theory is developed first followed by an illustrative

example. 

1. Introduction

With advancements in Micro-Electro-Mechanical Systems (MEMS) the estimate of vehicle

position and attitude has become affordable and increasingly accurate as sensors have become smaller

and cheaper. The integration of sensor outputs to achieve reasonable estimates remains an ongoing area

of research. Integrating INS with GPS is most commonly done with tight coupling, which integrates sensor

data at the measurement level, or loose coupling, which uses outputs of INS and GPS without knowledge

of internal variables. The approach presented here is loose coupling. There is an abundance of literature

pertaining to strapdown inertial navigation systems and integration with GPS [1]. However, much of this

literature is difficult to follow for the uninitiated. In order to greatly simplify complexity, the work

presented here will be in two dimensions. Additional simplifications include a flat, nonrotating earth and

short vehicle trajectories. These assumptions allow one to assume constant gravity, elevation, and

magnetic field, all on a two-dimensional flat earth. 

The remainder of this paper is contained in the following 4 sections. Section 2 covers the

mathematical development for attitude estimation around a single axis and INS-GPS fusion for derivation

of a platform’s 3 degrees-of-freedom (3DOF): (xnorth, xeast, ψ). Next, Section 3 describes the

experimental setup used to collect data for algorithm evaluation. Lastly, Sections 4 and 5 respectively

present results and discuss conclusions. 

2. Mathematical development

Reference frames play a large role in navigation. Before vehicle position and attitude can be

estimated, one must know what the estimates are relative to. Here, we use two reference frames: the

1
vehicle body frame (x,y,z), and the navigation frame which is fixed to the earth (North, East, Down), (NED). 

In the body frame, we have the origin as the rotational center of the vehicle, x-axis out the front, and yaxis out the right side. The frame of the vehicle is assumed to always be coplanar with the north and east 

vectors; the z axis is normal to this plan pointing down. Since the distances traveled by the vehicle are 

very small in relation to the earth, we will assume a flat earth with perpendicular lines of latitude and 

longitude. Also, we will assume there are no variations in the local gravitational and magnetic fields. These 

assumptions together with loose coupling of the GPS-INS greatly simplify the mathematical development. 

Also, the development presented here is for a vehicle that is assumed to be traveling in 2D. This is a 

reasonable assumption since the experimental vehicle trajectory was over a localized flat surface. 

 Estimation of vehicle position and attitude will be done in the navigation frame. It is assumed here 

that the senor frame and body frame are co-located and aligned. In order to be used for state estimation, 

the acceleration measurements onboard the vehicle must be converted from the measured body frame 

to the navigation frame. This may be accomplished by the simple matrix operations shown below in 

Algorithm 1 [2]. 

Algorithm 1

Figure 1 Diagram of two-dimensional vehicle body frame and navigation frame. Body frame is front of vehicle and right side (x, 

y) and navigation frame is North and East (N, E). 

2
|m| = √mx

(k)

2 + my

(k)

2 |1

mx

(k) =

mx

(k)

|m|

 my

(k) =

my

(k)

|m|

ψ(k) = atan2 (

−my(k)

mx(k) )

Cb

n

(k) = [

cosψ(k) − sinψ(k)

sinψ(k) cosψ(k) ]

f

b

(k) = [ fx

b(k)

 fy

b(k)

]

f

n(k) = [

fnorth

n

(k)

feast

n (k)

] = Cb

n

(k) ∙ f

b

(k)

Figure 1 shows the vehicle body frame x and y and the navigation frame of north and east. The 

angular displacement is the vehicle yaw, ψ, and the angular rate is ω, following the right hand rule with z

pointing down. In this algorithm the outputs from the two axis magnetometer, mx

(k) and my

(k), are 

normalized by division of the magnetic vector magnitude, |m|. The yaw angle is found using the MATLAB 

atan2 function which preserves the correct quadrant. Care should be taken that a negative sign is placed 

on the magnetic y measure to match yaw as defined in the figure. The Directional Cosine Matrix (DCM), 

Cb

n

, may then be formed using ψ. The DCM may then be used to compute the accelerations in the 

navigation frame. Throughout this algorithm, k denotes the integer time step, indicating that this 

algorithm is carried out every time step. The computed accelerations in the navigation frame, f

n

, may 

then be used to update the state vector in the Kalman filter, later described in Algorithm 3. Also, the yaw 

angle is combined with the gyro measures in the complementary filter to form an improved estimate of 

the vehicle attitude in Algorithm 2. 

 There are numerous schemes to estimate vehicle attitude based on sensor data. One common 

technique is the implementation of a Bayesian estimator such as an Extended Kalman Filter (EKF), an 

Unscented Kalman Filter (UKF), or a particle filter (PF). Because of its simplicity, stability, and accuracy, a 

special form of the complementary filter developed by Mahony [3] will be used. This algorithm is 

computationally efficient and well suited for real time processing using limited computational resources; 

it is shown below in Algorithm 2. 

Algorithm 2 

ωb(k) = [

0 −ω(k)

ω(k) 0

]

πa

(k) =

1

2

∙ (Ĉ

b

n

T

(k) ∙ Cb

n

(k) − Cb

nT

(k) ∙ Ĉ

b

n

(k))

Ω(k) = ωb(k) + Kp ∙ πa

(k)

 

1

 Normalization of the magnetic measurements is not necessary; however it may improve noise immunity. 

3
Ĉ

b

b

(k, k + 1) = I −

Ω(k) ∙ sin(|Ω(k)| ∙ Ts

)

|Ω(k)| +

Ω(k)

2

∙ (1 − cos(|Ω(k)| ∙ Ts

))

|Ω(k)|2

Ĉ

b

n

(k + 1) = Ĉ

b

n

(k) ∙ Ĉ

b

b

(k, k + 1)

T

 ψ(k) = atan2(Ĉ

b

n

2,1(k + 1), Ĉ

b

n

1,1(k + 1))

f

n(k) = [

fnorth

n

(k)

feast

n (k)

] = Ĉ

b

n

(k + 1) ∙ f

b

(k)

 In Algorithm 2, the matrix ωb contains the output of the z-axis gyro sensor in the body frame. The 

matrix Cb

n

 is that formed in Algorithm 1. The matrix Ĉ

b

n

 is an estimate of the true Directional Cosine Matrix 

formed by combining Cb

n

 of Algorithm 1 with the gyro reading at time step k. Ts is the sample rate. From 

Ĉ

b

n

 the yaw is computed every time step. The subscripts on Ĉ

b

n

2,1 and Ĉ

b

n

1,1 represent the second row, first 

column and the first row, first column values respectively. The equation for Ĉ

b

b

(k, k + 1) forms the time 

step update. If a simpler form of integration is desired, one can use Euler integration with little loss of 

accuracy. A more accurate estimation of the acceleration in the navigation frame is obtained by using 

Algorithm 2 compared to Algorithm 1 since both magnetometer and gyro inputs are considered. The gain, 

Kp, is set to one or greater and determines the relative contributions of the gyro and DCM to the attitude 

solution. In summary the inputs to Algorithm 2 are the gyro readings and DCM from Algorithm 1. The 

outputs are the attitude and an improved estimate of the DCM. The matrix Ĉ

b

n

 is used to transform the 

accelerations, which are measured in the body frame, to the navigation frame. Once in the navigation 

frame, the accelerations may be used as inputs into Algorithm 3 to help estimate vehicle displacement 

and velocity. 

 The acceleration expressed in the navigation frame may be combined with GPS measurements of 

position and velocity to perform state estimation of the vehicle in the navigation frame. This is done with 

a Kalman filter. In order to construct the Kalman filter, the vehicle dynamic equations are written in state 

space form. The state vector consist of the displacement, velocity, and acceleration in each component 

direction; that is north and east. 

x(k) = x(k − 1) + ẋ(k − 1) ∙ Ts + 

1

2

∙ ẍ(k − 1) ∙ Ts

2

ẋ(k) = ẋ(k − 1) + ẍ(k − 1) ∙ Ts

ẍ(k) = ẍ(k − 1) |2

In matrix form we have 

[

x(k)

ẋ(k)

ẍ(k)] = [

1 Ts

1

2

⁄ ∙ Ts

2

0 1 Ts

0 0 1 ] ∙ [

x(k − 1)

ẋ(k − 1)

ẍ(k − 1)]

The above equation may be written in compact notation. 

x(k) = Φ ∙ x(k − 1)

 

2

 The acceleration is assumed not to change over one time step. 

4
Here x is the state vector of displacement, velocity, and acceleration. 

In state space form, the Kalman filter may be expressed in orthogonal north and east directions as in 

Algorithm 3. Here, two channels are used to independently evaluate the state estimates, north and east. 

Algorithm 3

Prediction: 

xnorth

− (k) = Φ ∙ xnorth(k − 1)

Pnorth

− = Φ ∙ Pnorth(k − 1) ∙ ΦT + Q

Update: 

Δεnorth(k) = [

LatGps(k)

VelGpsnorth(k)

fnorth

n

(k)

] − xnorth

−

(k)

Knorth = Pnorth

−

∙ (Pnorth

− + R)−1

xnorth(k) = xnorth

− (k) + Knorth ∙ Δεnorth(k)

Pnorth(k) = (I − Knorth) ∙ Pnorth

−

xeast

− (k) = Φ ∙ xeast(k − 1)

Peast

− = Φ ∙ Peast(k − 1) ∙ ΦT + Q

Δεeast(k) = [

LongGps(k)

VelGpseast(k)

feast

n (k)

] − xeast

− (k)

Keast = Peast

− ∙ (Peast

− + R)−1

xeast(k) = xeast

− (k) + Keast ∙ Δεeast(k)

Peast(k) = (I − Keast) ∙ Peast

−

In Algorithm 3, the state to be estimated is x, and the measurements LatGps and LongGps are the latitude 

and longitude in meters displaced from an arbitrary starting point. These values are derived from the GPS 

readings3

. Likewise VelGps represents the respective velocities as measured from GPS. The matrices Q

and R are the process and measurement noise covariance, respectively. In practice Q and R are tuned for 

desired performance. Algorithm 3 is an iterative process; it is executed every time a new GPS reading is 

acquired. It consist of two steps: a prediction step (indicated by a minus superscript), which computes an 

estimate of the state x before the latest GPS measurements are used, and an update step which uses the 

latest GPS and acceleration measurements. An error term, Δε, is formed between the measured and the 

predicted state. The Kalman gain, K, informs the algorithm how much weight to place on the error update. 

This weighted error term is used to update the state estimation. Since this process is repeated every time 

a new GPS reading arrives, the algorithm tracks the evolution of the state through time. 

3. Experimental Setup 

Algorithms 1 to 3 were used to post-process data obtained from a small radio controlled (R/C) car 

shown in Figure 2. The car was driven in loops around a local area and data from the INS and GPS were 

recorded. The purpose of this experiment was to demonstrate the use of the algorithms and explore 

implementation issues for practical vehicle attitude and position estimation. 

 

3 Because of the local flat earth approximation, the conversion from degrees of Latitude and Longitude is scalar. 

For our location in Hampton, VA we have 110,980 meters per degree Latitude and 88,892 meters per degree 

Longitude. 

5
Figure 2. Radio controlled car used to gather data for navigation algorithms. Red circuit card is the VN200. 

The R/C car was an HPI Trophy Truggy Flux (part #:107018). A modification made to the vehicle 

was the replacement of the factory plastic body and wing with pieces of plywood in order to have an easy 

surface to mount the electronics. Suspension modifications were made to accommodate the extra weight 

from additional batteries and electronics. The factory 2-channel 2.4GHz radio was replaced with a 7-

channel Spektrum DX-7s; this allowed additional channels to be used in conjunction with the autopilot on 

board, which was not used for this experiment. The factory off-road tires and rims were replaced with 

slightly smaller but much smoother and more rigid street tires in order to increase traction on the paved 

surfaces. With all batteries and modifications, the R/C car weighed 12.45 lbs., which is approximately 2 

lbs. heavier than stock. 

For this experiment a BeagleBoard Xm and a VN-200 GPS/INS were used to collect data. The VN200 was configured to output a binary message stream which was saved to a USB flash drive for postprocessing in MATLAB. To reduce this noise in the accelerometer data, a layer of Dubro 1/4 inch latex 

foam rubber was mounted in between the VN-200 and BeagleBoard. Additional noise reduction was 

realized by using a noise reduction filter within the VN-200 architecture. 

6
4. Experimental Results

 The data recorded from the radio controlled car was post-processed to obtain estimates of vehicle 

attitude and position. Estimates were compared with those given by the VN200. Figure 3 shows the 

attitude of the car during a segment of the test. Since the test car was driven on a flat track, only yaw is 

shown. Pitch and roll are assumed to be negligible. 

Figure 3. Yaw angle in degrees. Blue is from Directional Cosine Matrix, red is from VN200, and green is from complementary 

filter. 

0 50 100 150 200 250 300 350

-200

-150

-100

-50

0

50

100

150

200

Time in Seconds

Yaw

DCM

VN200

Comp

7
Figure 4. Yaw angle in degrees using unwrap 

 As can be seen in Figure 3 and Figure 4, the estimates obtained from the DCM using the 

magnetometer only and those of the complementary filter closely follow the yaw computed by the VN200. 

Part of the difference may be attributed to the fact that the VN200 uses true north as the yaw reference 

while the estimates of the DCM and complementary filter use magnetic north. This would account for a 

difference of approximately 10°42′ for the test track location. A zoomed in section is shown in Figure 5. 

Here, one can see by comparing the DCM solution with the complementary filter that the complementary 

filter is combining the magnetometer and gyro measures to produce a smooth attitude estimate. 

0 50 100 150 200 250 300 350

-1200

-1000

-800

-600

-400

-200

0

200

Time in Seconds

Yaw

DCM

VN200

Comp

8
Figure 5. Zoomed in section of yaw estimation showing smoothing action of complementary filter (green line). 

 The attitude information is used by the inertial navigation system as described in Algorithm 3 to 

compute estimates of vehicle position. Vehicle position was computed in two independent channels using 

a Kalman filter to integrate the dynamic model with GPS readings to obtain a smooth position estimate. 

The plots shown in Figure 6 and Figure 7 are for the eastward velocity and position and those of Figure 8

and Figure 9 are northward. 

110 115 120 125 130 135 140 145 150 155

70

75

80

85

90

95

100

105

110

115

Time in Seconds

Yaw

DCM

VN200

Comp

9
Figure 6. Eastward Velocity 

Figure 7. Eastward Position 

0 50 100 150 200 250 300 350

-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

Time in seconds

Velocity Easting, m/s

Eastward Velocity

From GPS

Filter Estimate

0 50 100 150 200 250 300 350

-25

-20

-15

-10

-5

0

5

10

15

20

25

Time in seconds

Meters Easting

Eastward Position, Longitude From GPS

Filter Estimate

10
Figure 8. Northward Velocity 

Figure 9. Northward Position. 

0 50 100 150 200 250 300 350

-1.5

-1

-0.5

0

0.5

1

1.5

2

Time in seconds

Velocity Northing, m/s

Northward Velocity

From GPS

Filter Estimate

0 50 100 150 200 250 300 350

10

15

20

25

30

35

40

45

50

Time in seconds

Merers Northing

Northward Position, Lattitude

From GPS

Filter Estimate

11
 Figures 6-9 show the velocity and position tracking of the Kalman filter over a large portion of the 

test path. As can be seen, close tracking is maintained throughout. In order to see an advantage gained 

by incorporating INS with GPS, a zoomed in section of the vehicle location is shown in Figure 10 and Figure 

11. 

Figure 10. Figure illustrating improvement of Eastward position estimates, red, by Kalman filter. 

Figure 11. Figure illustrating improvement of Northward position estimates, red, by Kalman filter. 

594 596 598 600 602 604 606 608 610 612

9

10

11

12

13

14

15

16

17

18

19

Time in seconds Meters Easting

Eastward Position, Longitude

From GPS

Filter Estimate

418 420 422 424 426 428 430 432 434 436 438

35

36

37

38

39

40

41

Time in seconds Merers Northing

Northward Position, Lattitude

From GPS

Filter Estimate

12
Figure 10 and Figure 11 illustrate the ability of the Kalman filter to use the dynamic model to fill 

in position measurements between GPS updates. 

5. Conclusions 

 This tutorial describes a simple method for the integration of INS and GPS under a restrictive set 

of assumptions. The assumptions of a flat nonrotating earth with constant fields brought much 

simplification to the navigation problem. The purpose of this document is to offer a simple introduction 

into the complex field of strapdown inertial navigation systems and GPS/INS integration. The basic 

principles covered here may be extended to more complex navigation and filtering problems. 

Appendix A 

 In this appendix the Kalman filter used in Algorithm 3 will be derived from Bayes’ theorem [4]. 

Consider the state space equation show below. 

xk = Axk−1 + qk−1

yk = Hxk + rk

Here xk is the n-dimension state at time step k, A is the state transition matrix, yk is the vector of sensor 

measurements and H describes the connection between the states and measurements. The vectors 

qk−1 and rk are independent zero mean Gaussian process and measurement noise vectors respectively. 

Since the system described above is linear and excited by Gaussian noise, the state and measurement 

vectors must also be Gaussian. The probability density function, pdf, for a multivariate Gaussian 

distribution is shown below where we use the notation N( random variable ∣ mean, covarience ). 

N( x ∣ m, P ) =

1

(2π)

n

2

⁄

|P|

1

2

⁄

exp (−

1

2

(x − m)

TP−1(x − m))

In the above |P| is the determinant of the covariance matrix P and m is the mean state vector. The Kalman 

filter has two steps: prediction and update. In the prediction an estimate of the state at time step k is 

computed based on past measurements y1:k−1 before measurement yk is taken. In the update the 

measurement yk is used to update the state estimation. The prediction and update steps occur at every 

time step making the Kalman filter a recursive algorithm. The recursive formulas are derived below. 

 For the prediction at time step k, consider the joint pdf shown below conditioned on past 

measurements y1:k−1; here the equality follows from conditional probability. 

p( xk, xk−1 ∣ y1:k−1 ) = p( xk ∣ xk−1, y1:k−1 )p( xk−1 ∣ y1:k−1 )

Because the first term on the right hand side has no dependence on the measurements y, we may write 

the above equation as shown below. 

p( xk, xk−1 ∣ y1:k−1 ) = p( xk ∣ xk−1 )p( xk−1 ∣ y1:k−1 )

Both sides may be integrated over xk−1 to obtain the Chapman-Kolmogorov equation. 

13
∫ p( xk, xk−1 ∣ y1:k−1 )dxk−1

∞

−∞

= ∫ p( xk ∣ xk−1 )p( xk−1 ∣ y1:k−1 )

∞

−∞

dxk−1

The joint distribution now becomes the marginal distribution conditioned on past measurements y1:k−1. 

One can solve this equation to obtain the desired pdf for the predicted state mean mk

− and covariance Pk−. 

p( xk ∣ y1:k−1 ) = ∫ p( xk ∣ xk−1 )p( xk−1 ∣ y1:k−1 )

∞

−∞

dxk−1

The pdfs in the above equation are Gaussian and may be expressed as shown below. 

N( xk ∣ mk

−, Pk− ) = ∫ N( xk ∣ Axk−1, Q )N( xk−1 ∣ mk−1, Pk−1 )dxk−1

∞

−∞

Here Q is the process noise convenience matrix. Given the past values of mk−1and Pk−1from the previous 

time step the above equation may be solved for the current predicted mean and covariance values 

mk

− and Pk−. These values are estimates based on measurements up to and including k-1, but do not 

include the measurement yk . The solution is shown below4

. 

mk

− = Amk−1

Pk

− = APk−1A

T + Q

The above equations are predictions of the mean of the state xk and covariance at time step k, but without 

using the measurement yk. This is to say all measurements before k were used to give a best estimate of 

the state and covariance based on past measurements, state estimations and knowledge of the state and 

measurement distributions. The next step is to update this estimate using the information obtained in 

measurement yk. 

Bayes’ rule is used in the update step. Recall Bayes’ as shown below for random variables A and 

B. 

p( A ∣ B ) =

p( B ∣ A )p(A)

p(B)

In general Bayes’ rule can also be expressed with an additional conditional variable C as shown below. 

p( A ∣ B, C ) =

p(B ∣ A, C )p(A ∣ C)

p(B ∣ C)

For our derivation, we can use Bayes’ rule as shown below. Here we have an estimation of xk that 

includes yk for the update step. 

p( xk ∣ yk, y1:k−1 ) =

p( yk ∣ xk, y1:k−1 )p(xk ∣ y1:k−1)

p(yk ∣ y1:k−1)

 

4

 For the Gaussian linear system case the Chapman-Kolmogorov equation may be explicitly solved. Alternatively 

note that the expected value E{xk} = E{Axk−1 + qk−1} = Amk−1. Also the covariance cov{Axk−1 + qk−1} = 

E{(Axk−1 + qk−1 − Amk−1)(Axk−1 + qk−1 − Amk−1)

T

} = AE{(xk−1xk−1

T − mk−1mk−1

T

)}AT + Q = APk−1A

T + Q

where E is the expectation operator and the time dependence of Q has been suppressed. 

14
Since our derivation involves Gaussian distributions, we are only interested in computing the mean and 

covariance of the state update. To this end we will only examine the exponent of the distribution; which 

contains the mean and covariance, xk and Pk respectively. We can express this as the proportion shown 

below after recognizing that p( yk ∣ xk, y1:k−1 ) = p( yk ∣ xk

) since there is no dependency on y1:k−1 in 

our state space model. Also we have written yk, y1:k−1 in compact notation as y1:k. 

p( xk ∣ y1:k ) ∝ p( yk ∣ xk

)p(xk ∣ y1:k−1)

The pdfs in the right hand side of this equation are known. The first term is simply the measurement 

equation pdf and the second is that determined in the previous prediction step. We may now express the 

update step in terms of the respective Gaussian distributions. 

N( xk ∣ mk, Pk

) ∝ N( yk ∣ Hxk, R )N(xk ∣ mk

−, Pk−)

Here R is the measurement noise covariance matrix and mk

−, Pk− are the estimates computed in the 

prediction step. The goal is to find mk, Pk which are the updated mean and covariance based on the 

current measurement yk. One can do this by completing the square in the product of the two Gaussian 

pdfs given above. 

N( xk ∣ mk, Pk

) ∝

1

(2π)

n

2

⁄

|R|

1

2

⁄

exp (−1

2

(yk − Hxk)

TR−1(yk − Hxk

))

1

(2π)

n

2

⁄

|Pk

−|

1

2

⁄

exp (−1

2

(xk − mk

−)

T

(Pk

−

)

−1(xk − mk

−))

Consider the exponent only. The expansion is shown below. In this expansion, care must be taken to maintain the 

order since this is a matrix equation and multiplication does not commute. 

−

1

2

(yk

TR

−1yk − yk

TR

−1Hxk − xk

TH

TR−1yk + xk

TH

TR−1Hxk + xk

T

(Pk

−)

−1xk − xk

T

(Pk

−)

−1mk

−

− (mk

−)

T

(Pk

−)

−1xk + (mk

−)

T

(Pk

−)

−1mk

−

)

We need to shape the above expression into the form (xk − mk

)

TPk

−1(xk − mk

) in order to identify the 

updated mean mk and covariance Pk. In order to do this we gather like terms and drop any term that 

does not contain xk. The dropped terms can be ignored since they do not contribute to the expression for 

the Gaussian mean and covariance. Factoring out xk and recalling that the transpose of a covariance 

matrix is just the matrix itself we obtain the expression below. 

−

1

2

{xk

T

(H

TR−1H + (Pk

−

)

−1)xk − xk

T

(H

TR−1yk + (Pk

−

)

−1mk

−

) − (yk

TR

−1H + (mk

−

)

T

(Pk

−

)

−1)xk}

In order to make the math more tractable while completing the square we define some variables. 

V = HTR

−1H + (Pk

−

)

−1

b = HTR

−1yk + (Pk

−

)

−1mk

−

Using the above, we can write our expression as shown. 

xk

TVxk − xkTb − bTxk

In order to complete the square we can add a term to this expression as long as it does not contain xk. 

15
xk

TVxk − xkTb − bTxk + bTV

−1b

Also, recall that VV−1 = I can be inserted into any term to help complete the square. 

xk

TVxk − xkTVV−1b − bTVV−1xk + bTV

−1VV−1b = (xk − V−1b)

TV(xk − V−1b)

The fact that the above equality holds can be verified by expanding the right side to get the left. We now 

have the desired form and can identify the updated mean and convenience matrices. 

mk = V−1b = (H

TR−1H + (Pk

−

)

−1)

−1(H

TR−1yk + (Pk

−

)

−1mk

−

)

Pk = V−1 = (H

TR−1H + (Pk

−

)

−1)

−1

To make the computations more efficient, we use the matrix inverse lemma shown below where A,U,C

and V are simply placeholders to demonstrate the lemma. 

(A + UCV)

−1 = A−1 − A−1U(C−1 + VA−1U)

−1VA−1

Pk = ((Pk

−

)

−1 + HTR

−1H)

−1 = Pk

− − Pk−H

T

(R + HPk

−H

T

)

−1HPk

−

mk = {Pk

− − Pk−H

T

(R + HPk

−H

T

)

−1HPk

−}(H

TR−1yk + (Pk

−

)

−1mk

−

)

mk = Pk

−H

TR−1yk + mk

− − Pk

−H

T

(R + HPk

−H

T

)

−1HPk

−H

TR−1yk − Pk

−H

T

(R + HPk

−H

T

)

−1Hmk

−

We need to simplify the yk term. 

Pk

−H

TR−1 − Pk

−H

T

(R + HPk

−H

T

)

−1HPk

−H

TR−1 =

Pk

−H

TR−1 − Pk

−H

T

(R + HPk

−H

T

)

−1(R + HPk

−H

T

)R

−1 + Pk

−H

T

(R + HPk

−H

T

)

−1 =

Pk

−H

T

(R + HPk

−H

T

)

−1

Replacing the yk coefficient with the above simplified expression one may now write mk as shown below. 

The previously derived prediction equations and covariance update are also shown for completeness. 

Prediction step: 

mk

− = Amk−1

Pk

− = APk−1A

T + Q

Update step: 

mk = mk

− + Pk

−H

T

(R + HPk

−H

T

)

−1(yk − Hmk−)

Pk = Pk

− − Pk

−H

T

(R + HPk

−H

T

)

−1HPk

−

The equations are iterated every time step as in Algorithm 3. 

Appendix B 

 The purpose of this section is to provide MATLAB code to implement Algorithms 2 and 3. The code 

snippets are partial and are not intended to be functioning code, but are presented only as an aid to the 

16
reader. Both blocks of code run in a loop that executes every time step. The code for Algorithm 2 is shown 

below. 

 mnorm = norm([mx(k) my(k)]); % mx(k) my(k) from magnetometer

 mx(k) = mx(k)/mnorm; my(k) = my(k)/mnorm;

 psi = atan2(-my(k),mx(k)); 

 DCM = [cos(psi) -sin(psi);

 sin(psi) cos(psi)];

 YAW(k) = 180/pi*psi; %YAW from magnetometer

 gyro = [ 0 -wz(k); % wz(k) from z axis gyro

 wz(k) 0 ]; 

 omaga = 0.5*(Rh'*DCM - DCM'*Rh);

 OM = gyro + kp*omaga; % kp >= 1, gain constant

 nOM = norm(OM); %For exact solution

 A = eye(2) - OM.*sin(nOM*Ts)./nOM + (OM^2).*(1-cos(nOM*Ts))/nOM^2;

 Rh = Rh*A';

 yaw(k) = 180/pi*atan2(Rh(2,1),Rh(1,1)); %yaw, Complementary Filter

The above code returns Rh which is an improved estimate of the DCM. From Rh the yaw angle is 

calculated. 

The code for Algorithm 3 is shown below. This block is for the longitude channel. The inputs are the time, 

GPS longitude, GPS velocity, and eastward acceleration in the navigation frame. 

 Ts = Time(k)-Time(k-1); %Time step

 A = [1 Ts 0.5*Ts^2; 0 1 Ts; 0 0 1];

 x(:,k) = A*x(:,k-1); %Predict

 P = A*P*A' + Q;

 err(:,k) = [LongGpsE(k); VelNedE(k); accel_nav(2,k)] - x(:,k);

 S = P + R;

 K = P*(S\I);

 x(:,k) = x(:,k) + K*err(:,k); %Update

 P = (I - K)*P;

The estimated state x(:,k) combines both the INS and GPS measurements for improved estimates of 

position, velocity, and acceleration. The latitude channel follows in like manner. 

Bibliography 

[1] D. Titterton and J. Weston, Strapdown Inertial Navigation Techology, London, United Kingdom: The 

Institution of Engineering and Technology, 2004. 

[2] [Online]. Available: http://aerostudents.com/files/avionics/InertialNavigationSystems.pdf.

[3] R. Mahony, M. Euston, J. Kim, P. Coote and T. Hamel, "A Non-Linear Observer for Attitude 

Estimation of a Fixed-Wing Unmanned aerial Vehicle Without GPS Measurements," Transactions of 

the Institute of Measurement and Control, pp. 699-717, 2010. 

[4] S. Sa ̈rkka ̈, Bayesian Filtering and Smoothing, Cambridge, United Kingdom: Cambridge University 

Press, 2013. 

17
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

A Short Tutorial on Inertial Navigation System and Global Positioning 

System Integration 

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Smalling, Kyle M.; Eure, Kenneth W.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20602

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 04

Availability: NASA STI Program (757) 864-9658

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

The purpose of this document is to describe a simple method of integrating Inertial Navigation System (INS) information with 

Global Positioning System (GPS) information for an improved estimate of vehicle attitude and position. A simple two 

dimensional (2D) case is considered. The attitude estimates are derived from sensor data and used in the estimation of vehicle 

position and velocity through dead reckoning within the INS. The INS estimates are updated with GPS estimates using a 

Kalman filter. This tutorial is intended for the novice user with a focus on bringing the reader from raw sensor measurements 

to an integrated position and attitude estimate. An application is given using a remotely controlled ground vehicle operating in 

assumed 2D environment. The theory is developed first followed by an illustrative example. 

15. SUBJECT TERMS

GPS; Kalman Filtering; Navigation; Tutorial

18. NUMBER

 OF 

 PAGES

22

19b. TELEPHONE NUMBER (Include area code)

(757) 864-9658

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

 999182.02.60.07.03 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA-TM-2015-218803

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 09 - 2015