Keywords

Strap-Down; INS; Accelerometer; Gyro.

Acceleration Measurement

Different researches and text books had discussed the inertial navigation system, or some of its components such as [1-4]. The inertial navigation system accelerometer consists of a proof mass m suspended from a case by a pair of springs, as illustrated in Figure. 1. Referring to Figure. 1, the input axis is indicated by the arrow. If the body to which the system is attached is accelerating along this axis by an acceleration, a, the proof mass will be displaced from its equilibrium position by an amount, x. From Newton’s law, the force affecting on the proof mass is equal to the mass times the acceleration. And from Hook’s law, the force generated by the springs are equal to the equivalent spring constant, K, times the displacement x. Thus, from Newton’s law we find that

---

---

F = m a = K x ------(1)

Therefore the measured acceleration is

a = (K/m) x ------(2)

Consequently, acceleration in the plus direction causes the proof mass to move downward indicating positive acceleration. The displacement x is sensed using a pick-off and scaled to provide the value of acceleration. The proof mass equilibrium position is calibrated for zero acceleration. If the accelerometer is set on a bench in the earth’s environment, the gravity force will shift the proof mass downward w.r.t the case indicating positive acceleration. However, the gravitational acceleration is downward. For this reason, equation( 1 ) becomes

F = m (a-g) = K x ------(3)

And the measured acceleration is

a-g = (K/m) x = f ------(4)

Where, f is the specific force.

Commonly used Reference Frames

The earth centered earth fixed (ECEF) coordinates are defined such that the XE axis points towards the intersection between the Greenwich meridian and the earth’s equatorial plane, this point has a longitude and latitude (0,0) deg respectively [5]. The ZE axis of the ECEF axis is the earth’s rotational axis about itself. The YE axis is chosen to complete the right hand rule system.

The local-level (navigation) frame is defined such that X axis (E) points toward the east, Y axis (N) points to the north, and the Z axis (UP) points in the direction from the earth center to the INS origin. The ECEF and the navigation frame are illustrated in Fig. 2. The Geocentric-Inertial system has its origin at the earth’s center. The fundamental plane is the equatorial plane. The positive Z axis points in the direction of the North Pole, the positive X axis points in the direction of the vernal equinox on the first day of spring, and the positive Y axis completes the right handed set of coordinates. It is important to keep in mind that the inertial coordinate system is nonrotating w.r.t the stars (except for precession of the equinoxes).

The angle between the vernal equinox unit vector and Greenwich meridian is called αg –the “Greenwich sidereal time.” What we need is a convenient way to calculate the angle αg for any date and time of day. If we knew what αg was on a particular day and time we could calculate αg for any future time since we know that in one day the earth turns through 1.0027379093 complete rotations on its axis.

Suppose we take the value of αg at 0h UT on 1 January of a particular year and call it αgo. If, in addition, we express time in decimal fractions of a day, we can convert a particular day and time into a single number of days which have been elapsed since our “time zero.” If we call this number D, then

αg = αgo + 1.0027379093 x 360ox D [In Degrees]------(5)

or

αg = αgo + 1.0027379093 x 2π x D [In Radians]------(6)

Where

αgo = 1.74933340 rad at 1/1/1970 0h:0m:0s.------(7)

---

---

The transformation matrix from inertial to ECEF coordinates is now computed from, The body system of axes is defined in Figure. 4.

---

---

---

---

Theoretical measurements of a moving gyro triad

The actual measurements of the accelerometers are the sum of the nominal measurements and the errors. If the accelerometers are stationary and miss-leveled by r, p, and θ errors about the x, y, and z body axes, they will measure a part of the gravitational force according to the relation.



Where, Arpθ is the miss-alignment matrix. The nominal terrestrial usage of the gyros involves sensing the earth rotational components. As illustrated in Figure. 5, it is evident that the relation between the east velocity component and the longitude could be expressed as

---

---

Referring to Fig. 6, we could also write

---

---

Where h is the spacecraft height. Finally, the local-level frame angular velocity components are

In this treatment we have assumed that the body frame and the inertial frame are aligned together. The relationship between the navigation frame and the INS body frame is usually established through either a stationary alignment process, or by continuous external velocity information from GPS receiver. Thus, the attitude angles are used in generating the rotation between the b-frame and the navigation frame. The rotation rates measured by the gyros are used to constantly update the transformation matrix between the body frame and the navigation frame. This transformation is used to transform the measured acceleration from the body frame to the navigation frame (in the case of strap down INS). Once transformed, integrating the acceleration twice will provide the INS position difference w.r.t. an initial point. However, as mentioned before, accelerometer measurements include both the platform and gravity, and hence knowledge of the gravitational field is mandatory to determine the platform acceleration. It is obvious that the INS is fundamentally dependent on an accurate knowledge of the initial position, velocity, and attitude of the moving platform prior to the start of INS operation.

Modeling the INS

Mathematical notations

The spacecraft center of mass position is given by re = [xe ye ze]T. Transformation matrices are used to transform any navigation vector from one computational frame to another. As an example, the position vector in the e-frame (earth fixed frame), re, is transformed to the i-frame (inertial frame) to give ri through the transformation matrix Rie, and hence

ri = Rie, re------(15)

Transformation from the i-frame to the e-frame is established throughout the inverse transformation, Therefore,

re=(Rie)-1 ri= Rie ri------(16)

And the angular velocity of the e frame w.r.t. the i frame is denoted as ωie . If this angular velocity is expressed in the e-frame, it will be denoted ωiee . The rotation between two coordinate frames can be expressed as a sum of rotations between different coordinate frames as

ωiee=ωile+ωlee------(17)

Rotation can be represented in two ways, an angular velocity vector, and a skew symmetric matrix comprising the same vector components. For example, the angular velocity vector ωiee may be represented either in a vector form ωiee = [ ωx ωy ωz]T, or in the skew matrix form



Similarly, the velocity vector could be expressed as



It is also important to note that coordinate transformations apply also to angular velocity vectors. As an example, the angular velocity vector, , could be transformed to the inertial frame using the relation

ωiee=Rieωiee-----(20)

The equivalent transformation between two skew-symmetric matrices has the special form

Ωiie = Rie Ωiie Rei------(21)

If we have two vectors a and b, their cross product a b can be expressed in terms of their skew symmetric matrices ( A or B) as

a x b = A b = -Ba -----(22)

Time derivative of the rotation matrix

For a fixed point in a certain coordinate frame (body-frame) the position vector (r) can be transformed from the body frame, b, to the inertial frame, i, by the transformation

ri = Rib rb------(23)

Differentiation of both sides of equation ( 23 ) assuming that rb is constant (i.e. the point in spacecraft body is assumed to have a pure rotational motion without any translation) will give



The velocity of the body could be computed in the inertial frame from



Using equation ( 22 )we could write equation ( 25 ) in the form



Substitution using equation ( 21 ), into equation ( 26 ) gives



From equation ( 23 ) we could write equation ( 27 ) in the form



Since the matrix Rbi is the inverse of the Rib matrix, then RbiRib =I and equation ( 28 ) is reduced to



Comparing equation ( 24 ) with equation ( 29 ) we could easily deduce that

Calculation of position and velocity derivatives

The vector ra is transformed from the a coordinate frame to the b coordinate frame according to the relation

rb = Rb + Rba ra------(31)

Where: Rb is the vector directed from the b-frame origin to the a-frame origin, and Rba is the matrix that transforms from the a-frame to the b-frame. If the a-frame rotates w.r.t. the b-frame, then differentiation of equation ( 31 ) gives



Substitution using equation ( 30 ), by replacing the subscript b with a and the super script i with b into equation ( 32 ) gives



Equation ( 33 ) has the form of the coriolis theorem which states that, If two coinciding frames of reference experience relative angular rotation Ωaba(t) with rb(t)= Rbara(t), and Rba = RbaΩaba(t), then the time rate of change of the two coordinate systems are related by



Taking the second derivative of equation ( 33 ) gives



Substitution using equation ( 30 ) into equation ( 35 ) results in

place formula

From which we could write

Mechanization of inertial frame

The inertial body velocity can be expressed using coriolis equation using the following relation



Differentiating equation ( 37 ) and letting d/dt re= ve, we get



Application of equation ( 37 ) to the second term in the right hand side of equation ( 38 ) gives



But the earth turn rate is constant, therefore [d/dt]=0, and equation ( 39 ) becomes



Equation ( 40 ) could now be substituted instead of the second term of equation ( 38 ) and this substitution results in



From which it is evident that



Equation ( 4 ) says that a = f + g, where a = d2/dt2 ri. Consequently, equation (42 ) becomes



Looking carefully into equation (43 ) we could easily see that:

1-f : Is the specific force acceleration to which the INS is subjected.
2-ωie x ve : Is the coriolis acceleration.
3-ωie x [ωie x r]: Is the centripetal acceleration caused by the earth rotation with respect to an inertial frame.
4-ωie x [ωie x r] + g : is the centrifugal acceleration representing the local gravity vector, gl. So equation (43) could be written as



This equation could be expressed in the inertial frame of reference as

Mechanization of earth centered fixed frame

Using Coriolis equation, the ground velocity expressed in ECEF coordinate system, Vee , is differentiated with respect to time to give



Substitution from equation ( 45 ) results in



Expressing equation ( 48 ) in the ECEF axes gives

Mechanization of navigation frame

Mechanization of the navigation (local level) frame proceeds as the earth centered fixed frame mechanization. The derivatives of the ground velocity in the (local level) frame and the inertial frame of references are related by



Substitution from equation ( 45 ) into equation ( 50 ) gives



Equation ( 51 ) is expressed in the local level frame by

Mechanization equations for land based users

Referring to Figure.6, The earth rotation angular velocity vector is in the same direction of λ° , Thus for earth based users equations ( 12 ), ( 13 ), and ( 14 ) include the components of the earth’s angular velocity vector. As a result, let’s define



So equations ( 12 ), ( 13 ), and ( 14 ) could be put in the form

ωill=ωiel+ωell-----(54)

The first important note is related to the attitude mechanization equation. Applying equation ( 30 ) in the local level navigation frame results in:

Rl0b= Rlb Ωibb-----------------------(55)

But we know that

Ωlbb=Ωlib+Ωibb=-Ωilb+Ωibb=Ωibb-Ωilb------(56)

And

Ωilb=Ωieb+Ωelb-----(57)

Substitution from equation ( 57 ) into equation ( 56 ) gives

Ωlbb=Ωibb-(Ωebb+-Ωelb)------(58)

Consequently the attitude mechanization equations become

Rl0b= (Ωibb-Ωebb- Ωelb)------(59)

Where the Ωibb matrix is obtained from gyro measurements, Ωelb is the skew symmetric matrix corresponding to the vector Ωelb which could be obtained from



And similarly we could compute the skew symmetric matrix Ωieb from the corresponding angular velocity Ωieb vector which could be computed from



And Rib is the transformation matrix from the local level to the body frame defined by



With P, the pitch angle (inclination), r is the roll angle, A is the azimuth angle.

Finally, the transformation matrix from the local level to the earth fixed frame is given by

Quaternion attitude representation

Quaternion parameters are introduced to describe the rotation of the b-frame w.r.t. the local level frame. The transformation matrix that transforms a vector from the reference frame to the body frame is a proper real orthogonal one. This means that it preserves the lengths of vectors and the angles between them, and also the inverse of the matrix is simply its transpose. So if A is the transformation matrix that transforms a vector from the reference to the body coordinate system, its transpose AT will transform that vector back to the reference system.

It is also known that a proper real orthogonal 3x3 matrix has at least one eigenvector with eigenvalue unity [7] and hence there exists a unit vector, ê that isn’t changed by A so that

Aê=ê------(64)

Which means that the vector ê has the same components along both body and reference axes, so ê is the axis of rotation. This demonstrates Euler’s theorem: “the most general displacement of a rigid body with one point fixed is a rotation about some axis”.

Referring to the transformation matrices described earlier, for a single rotation about single axis, we find that: for a rotation with single angle Φ the trace (sum of all diagonal elements) of the corresponding transformation matrix is

tr (A(ɸ))=1 + 2 cos ɸ------(65)

Generalizing this rule for any sequence of rotations, we may transform any successive rotations to a single rotation by an angle Φ about a unit vector (axis of rotation) such that

cosΦ= ½ [tr(A) -1]------(66)

Examining this equation, we find that it has two solutions for the Euler axis / angle Φ, which differ in sign only. The two solutions have the rotation axis vectors ê in opposite directions.

This means that a rotation about the vector ê by an angle Φ is equivalent to a rotation about the vector -ê by -Φ.

The Euler symmetric parameters are defined in terms of Euler axis / angle and the components of the axis of rotation ê as follows

q1 = e1 sin (Φ/2 )------(67-a)
q2 = e2 sin (Φ/2 )------(67-b)
q3 = e3 sin (Φ/2 )------(67-c)
q4 = cos (Φ/2 ) ------(67-d)

The four Euler symmetric parameters satisfy also the constraint

q21 + q22 + q23 + q24 = 1------(68)

The four parameters could be regarded as the components of a quaternion vector,



The direction cosine matrix could be expressed in terms of quaternion as in equation (70).

And also we could inversely express Euler symmetric parameters in terms of the direction cosine matrix as follows

q1 = ±0.5 (1+ A11- A22 - A33)½------(71-a)
q2 = ±0.5 (1- A11+ A22 - A33)½------(71-b)
q3 = ±0.5 (1- A11- A22 + A33)½------(71-c)
q4 = ±0.5 (1+ A11+ A22 + A33)½------(71-d)

Attitude mechanization equations in terms of quaternion

Based on the pre-discussed principals we could write

ωibb=ωilb+ωlbb------(72)

Where ωibb is sensed by the gyroscopes. Thus,

ωlbb=ωibb+ωilb------(73)

And similarly,

ωilb=ωieb+ωelb------(74)

we know that, ωieb=Rblωiel and ωelb=Rblωell . Thus we could write equation (74) as



Substitution from equation (75) into equation (73) gives



In order to combine the parameters for two individual rotations we find that

[A(q'') ] = [A(q')][A(q) ]------(77)

The resulting quaternion q" could be found from the final attitude matrix, or it could be in another different way [7]. If q" represents the orientation of the body frame with respect to the reference frame at time t + Δt ,q represents it at time t then represents the orientation of the body frame relative to the position that it occupied at time t, thus we could write

q'1=eu sin (Δø/2)------(78-a)
q'2=ev sin (Δø/2)------(78-b)
q'3=ew sin (Δø/2)------(78-c)
q'v=cos(Δø/2)------(78-d)

Where eu, ev and ew are the components of the rotation axis unit vector (i.e. they are in the direction of p, q and r), ΔΦ = ωΔt =√p2+q2+r2 Δt (this is valid only for an infinitesimal time interval Δt), also we may write

q"=qq'

Or in a matrix form



Replacing q' by the change in Euler axis/angle



Thus we may possibly write



Where I is the fourth order identity matrix. In addition, all the right hand side terms are evaluated at time t. The former equation is mostly useful if the axis of rotation doesn’t change over the time interval Δt. For the case of attitude dynamics and prediction, equation (81) may be converted to a differential equation, though it may be used as it is to reduce the integration effort. Also small angles approximation will be used with the assumption of infinitesimal time period Δt to get

cos (ΔΦ/2) ≈1 and sin (ΔΦ/2) ≈ (ΔΦ/2) ≈ (ωΔt/2)

Noting that ω=ω ê=[ωx ωy ωz]T is the instantaneous angular velocity vector Finally we could write

q(t+Δt) ≈ [ID + Ω (Δt/2)]q(t)------(82)

And Ω is the skew-symmetric matrix



Where ωx, ωy, and ωz are the components of the angular velocity vector ωblb.

To go over the main points, the INS mechanization equations for land based users are summarized as

Position mechanization equation

The position of a moving platform described in the l-frame is expressed in terms of the curvilinear coordinates (Φ,λ,h), where Φ is the latitude, λ is the longitude, h is the altitude. Thus, we could write



And from equations (10) and (11) we could write

Velocity mechanization equation

Using equation (52), we could write

Thus,



Equation (86) implies the existence of a gravity model. To simplify the problem, we present the simplest possible gravity model defined by [6]

g(h) = g(0) / (1+h/R)2------(87)

Where, g(0) = 9.780318 m/s2, and g(h) is the magnitude of earth gravitational vector defined in the local level frame by g = [0 0 –g(h)]T. The local gravity vector is thus

gll = g-ωie x[ωie x r]------(88)

Attitude mechanization equations

The mechanization equations of attitude are given in equation (59) as R°lb = ( Ωbib-Ωbie-Ωbel ) or alternatively in equation (82) as q(t + Δt)≈[ID +Ω(Δt/2)]q(t), and Ωblb is given by its corresponding form of ωblb given by equation (76).

Modification of the mechanization equations for space based users

Looking carefully into equations (52), (53), and (54) we could deduce some important notes. First, the term ωlie, represents a component of the earth’s rotation sensed by a stationary INS (i.e. VE = VN = VUP = 0 m/sec) on the earth’s surface. For space based users, this term won’t be sensed by the INS at all. So, ωlie =0, and ωlil = ωlel. The second important note is related to the attitude mechanization equation. Applying equation (30) in the local level navigation frame results in:

Rl0b= Rlb Ωibb------(89)

But we know that

Ωblb=Ωbli+Ωbib=-Ωbil+Ωbib=Ωbib-Ωbil------(90)

And

Ωbil=Ωbie+Ωbel------(91)

Substitution from equation (91) into equation (90) gives

Ωblb=Ωbib-(Ωbie+Ωbel)------(92)

But Ωbie is the skew-symmetric matrix corresponding to ωbie which is equal to Rbl, and we also know that ωlie =0. Thus, equation (92) reduces to

Ωblb=Ωbib-Ωbel------(93)

Consequently the attitude mechanization equations become

Rlob = Rlb(Ωbib-Ωbel)------(94)

Where the Ωbib matrix is obtained from gyro measurements, Ωbel is the skew symmetric matrix corresponding to the vector ωbel which could be obtained from

Velocity mechanization equation

Using equation (52) in addition to that, ωlie=0 , we could write


Thus,



The gravity model given by equation (88) should now be modified to reflect the absence of earth rotation sensing by the accelerometers ( recalling that ωlie =0). So the gravity model becomes

gll = g −ωie × [ωie × r ]= g------(97)

Attitude mechanization equations

The mechanization equations of attitude are given in equation (94) as Rlb= Rlb(Ωbib-Ωbel) or alternatively in equation (1.99) as q(t + Δt)≈[ID + Ωblb(Δt / 2) ]q(t), and Ωblb is now given by equation (93) instead of its corresponding form of Ωblb given by equation (76).

Conclusion

The most important conclusion coming out from this research is that the application through which the SDINS is used affects the associated mechanization equations and the gravity model. For terrestrial applications, the SDINS model is developed clearly with the appropriate systems of axes. If the SDINS is utilized by a user who isn’t resting on the earth the SDINS model must be modified according to section (7).