Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Accuracy Improvement by Boundary Conditions for Inertial Navigation

Accuracy Improvement by Boundary Conditions for Inertial Navigation Hindawi Publishing Corporation International Journal of Navigation and Observation Volume 2010, Article ID 869127, 10 pages doi:10.1155/2010/869127 Research Article Accuracy Improvement by Boundary Conditions for Inertial Navigation Tuukka Nieminen, Jari Kangas, Saku Suuriniemi, and Lauri Kettunen Tampere University of Technology, Department of Electronics, Unit of Electromagnetics, P.O. Box 692, 33101 Tampere, Finland Correspondence should be addressed to Tuukka Nieminen, tuukka.nieminen@tut.fi Received 16 September 2009; Revised 8 January 2010; Accepted 12 May 2010 Academic Editor: Paul Cross Copyright © 2010 Tuukka Nieminen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The term inertial navigation is often automatically associated with the term initial value problem. However, there are many applications where it is possible to end up with a boundary value problem (BVP) as well. We show that in case of a BVP, the finite element method that incorporates boundary conditions can be efficiently used to compute position and velocity estimates not prone to accumulation of errors. For further accuracy enhancements, a method of combining inertial measurements with additional constraints is proposed. This way, we can model sensor errors, known to limit the accuracy of the system. The capabilities of the proposed methods are demonstrated with real-life examples. 1. Introduction denoted as v(t). Notice that the form (1) follows from the more general navigation situation by neglecting the rotation Typically, inertial measurements are made to have estimates and “curvature” of the Earth. of current position and velocity in real time. The set Now, from the theory of ordinary differential equations of equations used to compute the position and velocity (ODEs), we know that we need exactly two independent estimates out of the actual measurements depends greatly on constraints to unambiguously solve a 2nd order ODE (1). the application in hand. Equations for a general navigation Particularly, if these two independent constraints are both application are presented, for example, in [1]and [2–5]. In given at time t = T ,wehavean initial value problem (IVP) this study, however, we shall concentrate on the growing [6]: given a(t), v(T ),and r(T ), find r(t) such that 0 0 market of consumer applications employing inertial sensors within a suitable price (and quality) range. Thus, we can r ¨(t) = a(t), simplify the equations needed to compute the position and v(T ) = v , 0 0 (2) velocity estimates. This is so, because in this case the errors are more likely to be determined by the limited accuracy of r(T ) = r 0 0 the sensors rather than the accuracy of the used equations. In terms of the simplified equations, we basically need hold. Problems that are not IVPs based on the above to solve the following problem: given a(t): R → R ,find definition are called boundary value problems (BVPs). A r(t): R → R such that r satisfies special case of a BVP, convenient to our purposes, is stated as follows: given a(t), one of the constraints r ¨(t) = a(t). (1) In (1), “¨” represents the second time derivative, r(t)is v(T ) = v or r(T ) = r , 0 0 0 0 (3) the position of the object in a suitable coordinate frame, a(t) represents the acceleration of the object in the same and one of the constraints coordinate frame, and t is time. The velocity of the object (r ˙ ), given in the same coordinate frame as a(t)and r(t), is v(T ) = v or r(T ) = r , (4) 1 1 1 1 2 International Journal of Navigation and Observation find r(t)for T ≤ t ≤ T such that it satisfies how to exploit 1D finite element method (FEM) to solve the 0 1 underlying BVP with various possible choices of boundary r ¨(t) = a(t), conditions. At this stage, we assume that no additional constraints are given and that the measurements are exact. In v(T ) = v ,or r(T ) = r 0 0 0 0,(5) Section 5, we face the reality with faulty measurements and exploit linear additional constraints to enhance the accuracy v(T ) = v or r(T ) = r . 1 1 1 1 of the results. The underlying BVP is treated as exact, but We know for a fact that a problem of the form (5)has a the measurements are corrected using a linear sensor error model. As a result, we get two systems of linear equations: unique solution whenever the position is fixed at least at one boundary: This knowledge comes from the fact that (7)isa an exact one for the BVP and possibly an overdetermined one-dimensional case of a certain class of partial differential one for the parameters of the sensor error model. We will equations called Poisson problems, properties of which are solve these together to yield the corrected results. Finally, two well known [7]. real-life examples are discussed in Sections 6 and 7 before the Notice that there is a natural reason why problems (2) conclusions in Section 8. and (5) are distinguished. Namely, the solution methods for IVPs and BVPs differ substantially from each other 2. Preliminaries [6, 8]. This is also where this paper differs from numerous research articles considering inertial navigation: previously, Let us first motivate the chosen approach by means of the problems have been given in the form (2), whereas we examples: an application suitable to our approach is ski consider problems of the form (5). jumping, which has been the prime motivation of this study The key assumption of our approach is the following: [9]. As an inertial navigation problem, it includes the use inertial measurements related to a certain time period, includ- of consumer grade sensors with knowledge of the boundary ing a set of boundary values and possibly some additional values (in this case, position at both ends of the event). See constraints, are all available when processed. The length of the Section 7 for further details. time period (i.e., domain)[T , T ]in(5) can be anything 0 1 It turns out that especially sports applications tend to have from fractions of a second to some minutes. There are properties that are well suited to the considered approach: different ways we may end up with a BVP of the form the included actions are often periodic, in such a way that the (5): firstly, the underlying problem may naturally be a BVP. certain short-term action (e.g., a single step) repeats several Secondly, we may have an IVP, which we can pose as a BVP. times or some longer-term action a few times (e.g., a single This, of course, requires that we also know something about lap or a single jump as discussed earlier) during a certain the result at the end of the interesting event and that we event. In these kinds of applications, it is natural to encounter can afford to wait for the results until the end of the event. problems of the form (5) rather than (2): for example, in long In the first case, we generally do not have any additional jump (considering only the jump part), at time T (“take constraints we could use, but in the second case we end up off ”), the velocity of the shoe is known but position is not with a BVP and at least one additional constraint. We will and at time T (“landing”), the position of the shoe when it discuss examples of both cases in the next section. hits the surface of the sand can be accurately measured but The case where we only know the velocity at both ends the velocity is unknown. is, however, special and deserves some attention. According For a more general view, even a GPS-assisted inertial to the discussion above, we do not have a unique solution navigation system can be considered as a series of separate for this kind of a BVP although it fits the definition of our navigation periods with given boundary values rather than a model BVP (5). It is neither a valid IVP according to (2). As single event with additional constraints given at certain time it turns out, it is possible to solve also this without loss of instances. This is an example of an IVP, which can be posed generality. This is based on the fact that we can always fix the as a BVP. position at one boundary without changing the “shape” of the solution. Then, we will come up with a well-defined BVP with an additional constraint. 2.1. Some Remarks. There are two main classes of numerical The overall scope of this paper is to show how to treat methods for BVP’s. One class includes so-called shooting inertial navigation problems that are naturally (or know- methods and the other class methods of weighted residuals ingly) posed as a BVP of the form (5). Attitude computations such as FEM [8, 10]. We concentrate on the latter because are not considered, and where necessary, attitude is assumed of its property to minimize an error norm over the whole to be available. The form of problem (5)allowsustoconsider integration interval rather than minimizing only the local it as a set of three one-dimensional (1D) problems, rather error [7]. Another tempting property is that it concerns the than one three-dimensional (3D) problem. Notice that two position directly, giving us a possibility to more easily handle boundary values per a dimension are required to obtain a various types of additional constraints we will encounter unique solution. On the other hand, there is no need for the later on. boundary values for different dimensions to be of the same In many applications, inertial measurements are not the type. only source of information. In practice, however, the number This paper is organized as follows: in Section 2,pre- or the type of these additional constraints—combined with liminaries are discussed. In Sections 3 and 4, we will show the different kinds of a solution method—does not suggest International Journal of Navigation and Observation 3 the use of the traditional filtering (see Section 5 for details). [2, 15], whereas the proposed approach is a single- For these situations, we will introduce a computationally phase method with only few computations per an cheap and easily exploitable method to enhance the accuracy unknown. of the position and velocity estimates. The proposed method is based on sensor error modeling and is characterized by the Due to the significant differences in these two approach- following assumptions. es, we will in this context concentrate on the proposed method. Obviously, there are situations where the two (i) Sensor errors are modeled as constant errors. While methods could both be used. Comparison of the methods the behavior of a certain sensor error is in real life a in such a situation is interesting, although not addressed in stochastic process, one is usually able to fairly model this paper. it at least momentarily as a constant error. A suitable Finally, recall that problem (5) can be considered as a mathematical tool to characterize this is the Allan set of three 1D problems. Thus, let us focus on the 1D variance [11]. case for a while. Details of the more prevalent 3D case are considered later on. Notice that because of this, the symbols (ii) Considering consumer grade sensors, causes of the will be changed a bit: r , v ∈ R whereas r , v ∈ R and 0 0 0 0 most significant errors are usually known (e.g., bias so on. In the following treatment, it is assumed that the and scale factor error, both changing from turn-on to accelerometer samples represent an instantaneous value of turn-on [2]). the specific force. In other words, it is assumed that the sensor (iii) In particular, the sensor noise is not modeled. This output is not processed in any way before the “navigation is a conscious modeling decision to prevent unneces- computer.” sary “smoothing” of data. (iv) Additional constraints are treated as “exact”. That 3. Solution of a BVP is, the overall error in the additional constraints is assumed to be smaller than the error caused by the The main goal of this section is to form a linear system of simplification of the navigation equations. equations of the form Ax = b for the position estimates x. These are now expressed as a vector containing the position When additional constraints are available, problems at each discrete time instant (referred to as x ), where the resembling the ones considered here have previously been vector b is a function of a(t). In the following treatment, the resolved using the means of fixed interval smoothing (or total number of the samples is N , t ∈ [T , T ] is the value of i 0 1 “Kalman smoother”, if the type of the filter is fixed) [12– time instant i ∈ Z (1 ≤ i ≤ N ), and h = t − t . i i+1 i 15]. Considering inertial navigation, the most used methods Now, let us rephrase problem (5) as a variational of solving the problems are the two-filter smoother [16] equation and Rauch-Tung-Striebel smoother [17], used for example in [15, 18–20]. In both cases, the problem itself is posed as an T T 1 1 IVP and the basic idea is to run the filter in the forward ru ¨ dt = au dt ∀u ∈ U , (6) T T 0 0 direction as a “predictor” in phase one and then to run the filter in the backward direction while combining these two where U is a Hilbert space [7]. Second derivative of r can be results to yield the corrected result in phase two. Between the eliminated by integrating the left-hand side of (6) by parts, proposed method and the fixed interval smoothing method, which yields the fundamental difference is that in the proposed method the dynamics model is based on the BVP formulation and T 1 T T 1 1 1 in the previous approaches, on the IVP formulation with ru ¨ dt = ru ˙ − r˙u ˙ dt = au dt ∀u ∈ U , additional smoothing. T T T 0 T 0 0 Comparing the fixed interval smoothing technique to the (7) proposed method, in addition to the previously mentioned points, the most significant differences are as follows. where the notation “ ” stands for substitution. By evaluating (i) As fixed interval smoothing is run in both directions, the substitution term and rearranging (7), we get the filter needs two process models, which can be problematic [21]. As the proposed method is based T T 1 1 on the BVP formulation, it does not make a distinc- r˙u ˙ dt =− au dt + r˙(T )u(T ) 1 1 tion between forward and backward directions. T T 0 0 (8) (ii) In the filtering approach, each additional constraint − r˙(T )u(T ) ∀u ∈ U. 0 0 is assumed to be attached to a single time instance [2], while the proposed method does not make such While (8) is otherwise in a convenient form for our a restriction (see Section 5 for details). purposes, it is not well suited for practical computations, (iii) Fixed interval smoothing is a two-phase method because U is an infinite-dimensional space. Thus, let us requiring numerous computations per time step approximate U with a finite-dimensional space spanned by 4 International Journal of Navigation and Observation 1/h i−1 1/h i+2 φ φ φ φ 1 i i+1 N i−1 −1/h −1/h ··· ··· t i+1 h h h i−1 i i+1 t t t t t t t t 1 2 i−1 i i+1 i+2 N−1 N Figure 2: Time derivative of the lowest order basis functions φ Figure 1: Lowest order basis functions φ . (dashed line) and φ (dotted line). i+1 piecewise affine basis functions (“affine function = linear Now that we have discretized the problem, we are getting function + a constant”) closer to the equation Ax = b stated as our goal. In fact, 0, t< t , i−1 (12) is a system of linear equations. Thus, let us start by ⎪ assembling the matrix A. Knowing that the index i in (12) (t − t ) i−1 ⎨ , t ≤ t< t , i−1 i refers to a certain column and j to a certain row, the elements i−1 φ = (9) of A are given as ⎪ (t − t ) 1 − , t ≤ t< t , ⎪ i i+1 ⎪ h i T1 ˙ ˙ (13) A i, j = φ φ dt ∀i, j = 1, 2,... , N , 0, t ≥ t , i j i+1 often referred to as the “hat” functions. Functions (9)for few which are easy to compute by substitution of (11) into (13). values of i are shown in Figure 1. Also, let us use the same In practice, A is going to have only few nonzero elements basis functions φ to discretize u and r . This choice is often all of them at the diagonal, subdiagonal, and superdiagonal referred to as the Galerkin method [7]. Note that it is also (assuming the “obvious” indexing). For equally spaced nodes possible to choose basis functions φ different from the lowest with step size h, for example, we have elements 2/h at the order approximation used here, when considered necessary. diagonal and −1/h at the sub-and superdiagonal. Given the basis functions φ , notice that the position r Our next task is to compute the vector b.Asseenfrom can be approximated as a piecewise affine function (12), the task is to integrate term aφ over the interval [T , T ]. To do this, we fit a piecewise affine function to a(t) 0 1 r ≈ x φ . (10) i i between every node as seen in Figure 3 with dashed line, i=1 which yields From (8), we see that the derivatives of the basis functions j +1 j−1 (9) are also needed. It holds that aφ dt = a t +2a t j j−1 j ⎧ t j−1 (14) 0, t< t , ⎪ i−1 ⎪ 1 + 2a t + a t ⎪ j j +1 , t <t <t , i−1 i 6 ˙ h i−1 φ = (11) ⎪ for every j ∈ [2, N − 1]. For nodes j = 1and j = N,weget ⎪ − , t <t <t , i i+1 0, t> t . h i+1 1 aφ dt = [2a(t ) + a(t )], j 0 1 From Figure 2, one can verify that function φ is a piecewise i (15) constant function, discontinuous at points t , t ,and t . h i−1 i i+1 N−1 aφ dt = [a(t ) +2a(t )], j N−1 N From (8) we see that we need to integrate a similar term with N−1 discontinuities at the nodes over the domain. In other words, respectively. these discontinuities do not matter, which is a well-known Finally, let us consider how to apply the different types fact from integral theory. of boundary values into (12). At first, notice that terms In total, we are now in the position to discretize equation φ (T )and φ (T )seenin(12) will be nonzero (evaluating (8), which yields j 1 j 0 to value one) only for the values j = 1and j = N , T T 1 1 respectively. If r˙(T ) = v or r˙(T ) = v of (5)isgiven, 0 0 1 1 ˙ ˙ x φ φ dt =− aφ dt + r˙(T )φ (T ) i i j j 1 j 1 the respective boundary condition is called a Neumann T T 0 0 i=1 (12) boundary condition [7]. These can be applied directly by ˙( ) ( ) adding the given values into b. − r T φ T ∀j = 1, 2,... , N. 0 j 0 International Journal of Navigation and Observation 5 a(t ) i−1 An element-wise representation of C(t)is ⎡ ⎤ c (t) c (t) c (t) 11 12 13 ⎢ ⎥ ⎢ ⎥ ( ) ( ) ( ) C(t) = ⎢ c t c t c t ⎥ . (18) a(t ) 21 22 23 i+1 ⎣ ⎦ c (t) c (t) c (t) 32 32 33 a(t ) In this paper, we will assume that an estimate of C(t)for each time instant is available. Note that C(t) depends only on the attitude, which is independent of the position and velocity, when the assumptions considered in the first section ··· ··· t n are valid [2]. Furthermore, let vector g be the acceleration h h i−1 i due to the gravity represented in the navigation frame. With t t t t t t t 1 i−2 i−1 i i+1 i+2 N these notions, the specific force measurements f made by the Figure 3: Basis function φ along with the corresponding accelera- accelerometers can be “converted” into accelerations as tion values. n b b b n a = c f + c f + c f + g , 1 11 1 12 2 13 3 1 n b b b n a = c f + c f + c f + g , (19) 2 21 1 22 2 23 3 2 The other possibility for the boundary values is to have n b b b n r (T ) = r or r (T ) = r fixed thus having a Dirichlet 0 0 1 1 a = c f + c f + c f + g , 31 32 33 3 1 2 3 3 boundary condition [7]. As the value at the corresponding boundary is already known, it does not need to be solved. where the time dependencies have not been explicitly stated. Now, following from (19) and the right-hand side of Thus, we can move all terms depending on it to b reducing the number of unknown terms by one. For the reduced (12), we have an equation system, the corresponding term of the last two terms of (12) will be zero, as mentioned above. Recall that it is also n b j =− a φ dt i j possible to have a Dirichlet condition on one boundary and T (20) a Neumann condition on the other boundary. b b b n We have now means to assembly an equation of the form =− c f + c f + c f + g φ dt i1 i2 i3 j 1 2 3 i Ax = b (16) for the j th component of vector b ∀ i = [1,2,3]. Thus, the linear equation for the 3D position is ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ for the (1D) position of the object. Depending on the type of A 00 x b 1 1 1 the applied boundary conditions, the number of unknowns ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ n×n ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ n is equal to N − 1or N − 2. Matrix A ∈ R is known to be 0A 0 x b ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ . (21) 2 2 2 ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ symmetric and positive definite [22]. With the chosen basis 00 A x b 3 3 3 functions, A is also a tridiagonal matrix. In practice, these properties guarantee that (16) can be solved with linear time As discussed in the previous section, matrix A is a tridi- complexity [8]. In other words, when doubling the amount agonal matrix. Then, also the system matrix of (21)isnot of unknowns N , the time needed to solve the system is also only tridiagonal, but a block diagonal matrix. For instance, doubled (approximately), which is certainly not true for a if boundary conditions are fixed in the same way for all general system of linear equations. dimensions, it follows that A = A = A holds. 1 2 3 From the position data x,itisnow astraightforward With the equations derived in this section, it is possible task to compute the velocity using a suitable numerical to solve the 3D BVP using the presented FEM method. The differentiation formula. Since the position data is relatively solution of this more general problem can also be found with smooth due to the “double integration” process, we have not linear time complexity, as in the 1D case. experienced any problems in computing the derivative with an adequate accuracy. 5. Using a Number of Additional Constraints to Model Sensor Errors 4. Generalization to 3D So far we have constructed a method to compute position In this section, we will generalize the method described in and velocity data with certain boundary conditions. We will the previous section to the 3D case. For this, we introduce a now propose a way to exploit a number of additional con- coordinate transformation matrix C(t) used to transform the straints concerning the velocity or the position information b n specific force measurements f (t) into the accelerations a (t) to estimate sensor errors. For this, let us now precisely define represented in a suitable navigation frameasfollows: the term “additional constraint”: an additional constraint is a linear equation that bounds the position or the velocity of the n b a (t) = C(t)a (t). (17) object at an arbitrary number of time instances t ∈ [T , T ]. i 0 1 6 International Journal of Navigation and Observation At simplest, the previous definition could simply mean that they should be treated as unknowns. Other typical errors are v (t ) is fixed. On the other hand, a less intuitive equation caused by misalignment of the axes, changing temperature 1 i and drifting bias. [r (t )+ r (t )] = z is also a valid constraint. This way, 1 i 2 i 1 i=1 As a practical example used in the example 1 later in we can for example exploit constraints like x = x for any i i j Section 6, let us model the sensor errors as follows: and j without saying anything about the absolute position at these points. It could be hard to exploit these kind of f (t) = sf (t) + b, (25) constraints properly in traditional filtering problems. In this section, we will use constraints systematically to enhance accuracy. To make the concepts presented in this where f (t) is the corrected specific force measurement, s is section clear, let us first present the equations for the 1D case. some constant scaling factor, and b is some constant bias term correcting the erroneous measurement. By replacing the measured f (t) treated in the previous section with the 5.1. Treating Additional Constraints. Let us first assume N×N that A ∈ R holds. That is, given a Dirichlet bound- corrected acceleration f (t), it is easy to form an equation of ary condition, we add an auxiliary equation to the system the form instead of removing the corresponding boundary value from b = Fl (26) the system. This is a necessary procedure when exploiting additional constraints. for vector b seen in (16). Matrix F is in this case an N × In general, all linear constraints can be represented in the 2 matrix, whose elements are formed by computing the form T T 1 1 integrals − aφ and − φ from (12). Vector l is simply j j T T 0 0 Dx = e, (22) sb [ ] . Let us now replace the right-hand side of (16)with p×N p where D ∈ R and e ∈ R hold, where p is the number of the constraint equations. For a practical example closely b + b = b + Fl, (27) 0 0 related to the example 1, consider that (5)issolvedwith Dirichlet boundary values. If the Neumann boundary values where b depends only on the given boundary conditions are also known and the samples are equally spaced, one can (which are treated as exact) and corresponding boundary construct (22)as values. This distinction is necessary, since one should ⎡ ⎤ not modify the given boundary values by applying error 3 2 1 − − 0 ··· 00 0 0 modeling on them. ⎢ ⎥ 2h h 2h ⎢ ⎥ D = (23) In general, a problem of linear constraints used to model ⎣ ⎦ 1 2 3 000 0 ··· 0 − linear sensor errors can be stated as 2h h 2h Ax = b + Fl, ⎡ ⎤ (28) ⎣ ⎦ Dx = e, e = . (24) N×q q where F ∈ R , l ∈ R hold, and q is the number of Matrix D is formed using a three-point differentiation modeled sensor error terms. method based on polynomial interpolation [10]todeter- Now, since A is known to be invertible, one gets an mine Neumann boundary values from the displacement equation data. −1 −1 These constraints could be exploited just by computing DA Fl = e − DA b (29) the linear least squares problem constructed by adding these −1 additional equations to (16). Usually, however, it is favorable for l. Note that the matrix DA F has dimension p × q, that to use the additional constraints to model sensor errors,at is, it is typically a very small matrix compared to A.With l least when one has any knowledge of the types of the errors known, it is easy to solve for x. included in the measurements. For more information on In the case where p> q (more constraints than the assumptions related to the sensor error modeling, recall model parameters), one can also take reliability of different Section 2.1. In this text, we will present a linear sensor error measurements into account by solving a weighted least model. squares problem (WLS) [22]. In general, (29) has a unique solution l only when p = q holds and the row rank of D is full. Typically, making sure that p ≥ q and that the row 5.2. Sensor Error Model. In typical situations, bias offset rank of D is at least q,(29) has either a unique or a least (i.e., the sensor shows nonzero output when no forces are squares solution for l. In each case, the upper equation of acting upon it) and scaling factor are the two terms which (28) is treated as an exact equation. have to beknown very accurately in order to have any realistic position or velocity estimate as a solution. Because of the run-to-run variations of these errors, they cannot be 5.3. Additional Constraints in 3D. Let us now consider assumed to be constant between two separate events. Thus, the use of additional constraints in a 3D case. Using the International Journal of Navigation and Observation 7 notions form the previous section and (22) as an example, constraints can be represented as ⎡ ⎤ a true ⎢ ⎥ ⎢ ⎥ D⎢ x ⎥ = e. (30) meas ⎣ ⎦ ω(t) In addition to the sensor errors discussed in the 1D example, it is now also possible to model errors like attitude errors, unknown value of the acceleration due to the gravity, and cross-correlation of different sensors. A particularly useful method is to treat the acceleration due to the gravity Figure 4: Measurement setup of example 1. as an unknown three-dimensional vector, which is then subtracted from the specific force measurements given in the global coordinates. This reduces the systems sensitivity to where only the relevant indices are shown. Thus, we get size initial attitude errors. ⎡ ⎤ As an example, let us now derive the equations for sensor p p p p p p p0 0 1111 111 1211 121 1311 131 ⎢ ⎥ error of (25) for each axis in addition to the unknown accel- ⎢ ⎥ F = ⎢ p p p p p p 0p0⎥ , (34) 2111 211 2211 221 2311 231 eration due to the gravity discussed above. Thus, we have a ⎣ ⎦ total of 9 unknown model parameters. As one could expect, p p p p p p 00 p 3111 311 3211 321 3311 331 we must take the effects of the coordinate transformation matrix into account when deriving the necessary equations. n n n (35) l = s b s b s b g g g . 1 1 2 2 3 3 1 2 3 By plugging (25) into (21), we have For the nine unknown model parameters we need at least b b nine constraints. This can be covered, for example, with the b j =− c s f + b + c s f + b i i1 1 1 i2 2 2 1 2 knowledge of the velocity of the object at three different points. In total, we have equations identical to (28), with b n − c s f + b + g φ dt i3 3 3 j 3 i dimension N replaced by 3N . T T 1 1 =−s c f φ dt − b c φ dt 1 i1 j 1 i1 j 6. Example 1 T T 0 0 (31) T T To demonstrate the use of the proposed method, an example 1 1 − s c f φ dt − b c φ dt 2 i2 j 2 i2 j 2 using readily available consumer grade accelerometers [23] T T 0 0 was created. An accelerometer was mounted on a horizon- T T 1 1 tally rotating rate table, whose angular velocity (rotation − s c f φ dt − b c φ dt 3 i3 j 3 i3 j rate) can be controlled. The accelerometer was mounted in T T 0 0 such a way that its measurement direction was approximately the same as the direction of the tangential acceleration caused − g φ dt. by angular acceleration of the rate table, as seen in Figure 4. The aim of this example was not to get position and The several integrals in (31) are scalars for each j ∈ velocity as accurately as possible, but to compare different [1, 2,...N ], easily evaluated with (14)and (15), since no solution methods in a situation where one needs to get unknown terms appear inside the integrals. Let us now reasonable position and velocity estimates regardless of the gather the results into vectors fact that the measurement contains significant errors. The angular velocity of the rate table was set to follow function m b illustrated in the Figure 5 a certain number of times. With p j =− c f φ dt (32) klmn j kl l this kind of a setup, the true acceleration, velocity, and displacement of the mounted sensor were known up to the for all j ∈ [1, 2,... , N ]. With these notions, we get, for accuracy of the motor rotating the rate table. The errors example caused by the rate table were observed using the angular rate output of the motor and found to be negligible compared to p j =− c f φ dt, other error sources. kl11 kl j In first test, the rate table was set to repeat exactly the T1 function representedinFigure 5 ten times, leading to seven (33) p j =− c φ dt, kl1 kl j full revolutions (or 10.4meters) in T = 20 seconds. In the other test, the function of the same form was scaled in such T1 a way that 50 repeats lead to total of 62 full revolutions (or p j =− φ dt, 92.3meters) in T = 100 seconds. The sampling frequency of 0 8 International Journal of Navigation and Observation Table 1: Comparison of different methods computing object’s velocity and position. “I“ stands for time-stepping, “II“ for FEM, and “III“ for FEM with sensor error modeling. Test 1 T est 2 Vel. error [m/s] t t t t t t 1 2 3 1 2 3 I0.07 0.10 0.12 33.24 65.47 97.17 II 0.01 0.04 0.06 −0.67 0.50 1.15 III 0.04 0.05 0.04 0.33 0.54 0.44 Pos. error [m] t t t t t t 1 2 3 1 2 3 I0.10 0.49 1.00 510.96 1992.99 4434.70 II −0.31 −0.29 −0.15 −42.34 −43.42 −16.59 III −0.10 −0.07 −0.05 −6.78 −4.72 −0.61 9 10 −5 −10 11.52 2.53 00.51 1.52 Time (s) Time (s) Figure 6: Ideal (thick line) and measured acceleration (thin line). Figure 5: Angular velocity of the rate table. the sensor was set to 1000 Hz and an accurate (16 bit) analog three separate points by computing the difference between to digital converter was used. Accelerations were measured the computed and the real value. Point t was located at T/3, using a consumer grade 12 g accelerometer [23]. t at T/2, and t at 2T/3. 2 3 The raw accelerometer data was mapped to accelerations As seen in Table 1, method I seems to increase the error with the scale factor given by the manufacturer. From this with increasing time, as expected. Method II on the contrary, acceleration data, position and velocity were computed by a thanks to the basic property of the variational technique, number of different methods: does not increase the error but distributes it over the whole time period. The difference between these two methods is (i) “Traditional” IVP (double integration with trapezoid clearly seen in the velocity and position errors of the longer rule), test (test 2), where method II gives much better estimates (ii) FEM with Dirichlet boundary conditions (BCs), than method I. In each test, method III clearly outperforms (iii) FEM with Dirichlet BCs combined with additional methods I and II, which is expected due to the provided two Neumann BCs to supply the error model (25). additional constraints. Figures 7 and 8 demonstrate the differences between In methods I and II, the bias error of the sensor was velocity and position estimates given by methods I and II estimated by averaging the output while the sensor was at during test 1. Estimates given by method III coincide with rest. This was done in order to make method I (and to the reference plots. Figure 7 shows only the last spin of the some extent, method II) comparable to the method III by test 1 for better view of the differences. reduction of the large bias error. This is rarely possible in general and only method III can be used to reliably detect any remaining bias (example 2 in Section 7 is an example of 7. Example 2 this case). As a reference, the ideal (ideal driving motor) and the measured accelerations (accelerometer bias removed) of This example considers the computation of the velocity and the first spin are plotted in Figure 6. position of a ski jumper during a single jump. As compared Table 1 shows the main results of the tests. In each test, to the previous example, this is a more realistic and general the computed results were compared to the known value in inertial navigation problem with six degrees of freedom. ω (rad/s) Acceleration (m/s ) International Journal of Navigation and Observation 9 2.5 40 1.5 −20 0.5 −40 0 −60 −0.5 −80 19 19.520 20.5 −50 0 50 100 Time (s) x (m) Figure 7: Velocity estimates given by methods I (thin black line) Figure 9: Computed two-dimensional trajectories of two indepen- and II (thin gray line) compared to the ideal velocity (thick line). dent events (thin gray and thin black lines) along with the known profile of the hill (thick black line). take off and eventually land at some point of the landing hill. Drawn circles at both ends of the trajectories demonstrate the start and the end of the navigation period and the given Dirichlet boundary conditions. Drawn squares represent the landing points, in this case quite different for the two events, agreeing well with the recorded jump lengths. Notice that the trajectories are computed using two distinct measurement systems, each containing unique error sources. Figure 10 shows that the vertical component of the velocity of the same two jumps is plotted as a function of horizontal displacement. Notice how the absolute values −2 of the vertical velocity substantially differ while the small 0 5 10 15 20 25 changes in the velocity are practically identical. One might Time (s) first consider the small changes as typical stochastic errors caused by the inertial navigation system, especially when Figure 8: Position estimates given by methods I (thin black line) and II (thin gray line) compared to the ideal position (thick solid dealing with low performance sensors. line). Given that two independent measurement systems show the same variations in the velocity at the same locations, it is evident that there is actually only a negligible amount of stochastic errors present. Instead, the small variations are Computation of the attitude of the object is based on the data caused by deterministic sources, namely, in this particular given by three standard consumer grade gyroscopes [24]. application the uneven inrun hill. The position and velocity were then measured with three Unfortunately, the estimates cannot be compared with standard consumer grade accelerometers [23]and computed a reference trajectory, because such data are not available. with the proposed method with a sensor error model similar Thus,wecannotgivethe exactamountoferror presentinthe to the one presented in Section 5.3. The needed additional position and velocity estimates. We do however claim that the constraints contained information about achieved accuracy is something one does not typically expect (i) the location of the jumper at five points evenly spread from consumer grade inertial sensors. on the inrun hill (in Figure 9, the part of thick black line with negative x-and positive y-coordinates) 8. Conclusion (ii) the trajectory of the jumper after the landing, which The work was motivated by applications, where it is natural should coincide with the linearization of the landing to encounter BVPs instead of IVPs. In many cases, it is also hill (in Figure 9, the part of thick black line with x- possible to formulate an IVP as a BVP, given that the results coordinates greater than 100 m). are not required in real time. In Figure 9, two-dimensional trajectories of two jumpers Finite element method is utilized to solve inertial navi- are plotted along with the known profile of the hill. At first, gation problems formulated as BVPs. As a result, we get a the trajectories follow the profile of the hill until the jumpers linear system of equations for the position estimates, whose Velocity (m/s) Distance (m) y (m) 10 International Journal of Navigation and Observation −4 [4] M.S.Grewal, L. R. Weill, andA.P.Andrews, Global Positioning Systems, Inertial Navigation, and Integration, John Wiley & Sons, New York, NY, USA, 2001. [5] J.A.Farrell andM.Barth, The Global Positioning System & −6 Inertial Navigation, McGraw–Hill, New York, NY, USA, 1999. [6] L. F. Shampine, I. Gladwell, and S. Thompson, Solving ODEs with Matlab, Cambridge University Press, Cambridge, UK, −8 [7] S. Larsson and V. Thome, Partial Differential Equations with Numerical Methods, Springer, New York, NY, USA, 2003. [8] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetter- −10 ling, Numerical Recipes in C (The Art of Scientific Computing), Cambridge University Press, Cambridge, UK, 2nd edition, −12 [9] M. Virmavirta, T. Nieminen, T. Tarhasaari, and L. Kettunen, −70 −60 −50 −40 −30 −20 −10 0 “High precision inertial measurement for tracking the trajec- x (m) tories of ski jumpers “Smart Boot Project”,” in Proceedings of the 13th Congress European College of Sport Science, p. 572, Figure 10: Vertical velocity (v ) as a function of horizontal dis- Estoril, Portugal, July 2008. placement (x) of the two independent measurements. [10] D. Kincaid and W. Cheney, Numerical Analysis: Mathematics of Scientific Computing, BROOKS/COLE, Florence, Ky, USA, 3rd edition, 2002. [11] N. El-Sheimy, H. Haiying, and N. Xiaoji, “Analysis and mod- solution can be found with linear time complexity. It is eling of inertial sensors using allan variance,” EEE Transactions demonstrated that solving a BVP rather than an “equivalent” on Instrumentation and Measurement, vol. 57, no. 1, pp. 140– 149, 2008. IVP gives more accurate results. [12] A. Gelb, Applied Optimal Estimation, MIT Press, Cambridge, For further accuracy enhancements, an efficient way of Mass, USA, 1974. combining inertial measurements with possible additional [13] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with constraints is created. This gives us a possibility to model Applications to Tracking and Navigation, Chapman & Hall / constant sensor errors, known to limit the achievable CRC, London, UK, 2004. accuracy of the system. While the error model significantly [14] J. L. Crassidis and J. L. Junkins, Optimal Estimation of Dynamic enhances the accuracy of the system, it is kept computation- Systems, Chapman & Hall / CRC, London, UK, 2004. ally simple and easily adoptable. [15] P. D. Groves, Principles of GNSS, Inertial, and Multisensor In practice, the accuracy improvements allow us to Integrated Navigation Systems, Artech House, Norwood, Mass, exploit inertial sensors of certain performance level in more USA, 2008. challenging applications. For this, it is necessary to see that [16] D. Fraser and J. Potter, “The optimum linear smoother as the concept of inertial navigation does not invariably imply combination of two optimum linear filters,” IEEE Transactions on Automatic Control, vol. 14, no. 4, pp. 387–390, 1969. an IVP, but a BVP as well. Then, the use of FEM will provide [17] H. E. Rauch, F. Tung, and C. T. Striebel, “Maximum likelihood an efficient way to compute position and velocity estimates estimates of linear dynamic systems,” AIAA Journal, vol. 3, no. not prone to the accumulation of errors. 8, pp. 1445–1450, 1965. In larger scale, the current paper serves as an introduc- [18] K. Gade, “Navlab, a generic simulation and postprocessing tion to the idea of formulating inertial navigation problems tool for navigation,” European Journal of Navigation, vol. 2, no. as BVPs. As a consequence, further studies are needed to 4, pp. 21–59, 2004. address problems to which the presented tools do not provide [19] Y. Yang, Z. Jin, W. Tian, and F. Qian, “Application of fixed an obvious solution. These include, for example, stochastic interval smoothing to gps/dr integrated navigation system,” in errors, reliability of the possible additional constraints (as Proceedings of the IEEE Intelligent Transportation Systems, vol. compared to the accuracy of the IMU), and coupling of 2, pp. 1027–1031, October 2003. position and attitude errors. [20] A. B. Willumsen and ∅. Hegrenæs, “The joys of smoothing,” in Proceedings of the IEEE Bremen: Balancing Technology with Future Needs (OCEANS ’09 ), pp. 1–7, Bremen, Germany, May References [21] M. Klaas, M. Briers, N. de Freitas, A. Doucet, S. Maskell, and D. Lang, “Fast particle smoothing: if i had a million particles,” [1] I. Y. Bar-Itzhack, “Navigation computation in terrestrial strapdown inertial navigation systems,” IEEE Transactions on in Proceedings of the 23rd International Conference on Machine Learning (ICML ’06), vol. 148, pp. 481–488, Pittsburgh, Pa, Aerospace and Electronic Systems, vol. 13, no. 6, pp. 679–689, 1977. USA, 2006. [22] J. W. Demmel, Applied Numerical Linear Algebra, SIAM, [2] D.H.Titterton andJ.L.Weston, Strapdown Inertial Navigation Technology, Institution of Engineering and Technology, Lon- Philadelphia, Pa, USA, 1st edition, 1997. [23] VTI. VTI SCA620-CHCV1A Datasheet, 2006. Revision 2/2, don, UK, 2nd edition, 2004. [3] A. B. Chatfield, Fundamentals of High Accuracy Inertial Nav- Checked on 22.12.2009. [24] Analog Devices. AD ADXRS300ABG Datasheet, 2004. Revision igation, American Institute of Aeronautics and Astronautics, Reston, Va, USA, 1997. B, Checked on 22.12.2009. v (m/s) y International Journal of Rotating Machinery International Journal of Journal of The Scientific Journal of Distributed Engineering World Journal Sensors Sensor Networks Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Volume 2014 Journal of Control Science and Engineering Advances in Civil Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Submit your manuscripts at http://www.hindawi.com Journal of Journal of Electrical and Computer Robotics Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 VLSI Design Advances in OptoElectronics International Journal of Modelling & Aerospace International Journal of Simulation Navigation and in Engineering Engineering Observation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2010 Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com http://www.hindawi.com Volume 2014 International Journal of Active and Passive International Journal of Antennas and Advances in Chemical Engineering Propagation Electronic Components Shock and Vibration Acoustics and Vibration Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png International Journal of Navigation and Observation Hindawi Publishing Corporation

Accuracy Improvement by Boundary Conditions for Inertial Navigation

Loading next page...
 
/lp/hindawi-publishing-corporation/accuracy-improvement-by-boundary-conditions-for-inertial-navigation-5MrsIJ4vIo

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Hindawi Publishing Corporation
Copyright
Copyright © 2010 Tuukka Nieminen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
ISSN
1687-5990
DOI
10.1155/2010/869127
Publisher site
See Article on Publisher Site

Abstract

Hindawi Publishing Corporation International Journal of Navigation and Observation Volume 2010, Article ID 869127, 10 pages doi:10.1155/2010/869127 Research Article Accuracy Improvement by Boundary Conditions for Inertial Navigation Tuukka Nieminen, Jari Kangas, Saku Suuriniemi, and Lauri Kettunen Tampere University of Technology, Department of Electronics, Unit of Electromagnetics, P.O. Box 692, 33101 Tampere, Finland Correspondence should be addressed to Tuukka Nieminen, tuukka.nieminen@tut.fi Received 16 September 2009; Revised 8 January 2010; Accepted 12 May 2010 Academic Editor: Paul Cross Copyright © 2010 Tuukka Nieminen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The term inertial navigation is often automatically associated with the term initial value problem. However, there are many applications where it is possible to end up with a boundary value problem (BVP) as well. We show that in case of a BVP, the finite element method that incorporates boundary conditions can be efficiently used to compute position and velocity estimates not prone to accumulation of errors. For further accuracy enhancements, a method of combining inertial measurements with additional constraints is proposed. This way, we can model sensor errors, known to limit the accuracy of the system. The capabilities of the proposed methods are demonstrated with real-life examples. 1. Introduction denoted as v(t). Notice that the form (1) follows from the more general navigation situation by neglecting the rotation Typically, inertial measurements are made to have estimates and “curvature” of the Earth. of current position and velocity in real time. The set Now, from the theory of ordinary differential equations of equations used to compute the position and velocity (ODEs), we know that we need exactly two independent estimates out of the actual measurements depends greatly on constraints to unambiguously solve a 2nd order ODE (1). the application in hand. Equations for a general navigation Particularly, if these two independent constraints are both application are presented, for example, in [1]and [2–5]. In given at time t = T ,wehavean initial value problem (IVP) this study, however, we shall concentrate on the growing [6]: given a(t), v(T ),and r(T ), find r(t) such that 0 0 market of consumer applications employing inertial sensors within a suitable price (and quality) range. Thus, we can r ¨(t) = a(t), simplify the equations needed to compute the position and v(T ) = v , 0 0 (2) velocity estimates. This is so, because in this case the errors are more likely to be determined by the limited accuracy of r(T ) = r 0 0 the sensors rather than the accuracy of the used equations. In terms of the simplified equations, we basically need hold. Problems that are not IVPs based on the above to solve the following problem: given a(t): R → R ,find definition are called boundary value problems (BVPs). A r(t): R → R such that r satisfies special case of a BVP, convenient to our purposes, is stated as follows: given a(t), one of the constraints r ¨(t) = a(t). (1) In (1), “¨” represents the second time derivative, r(t)is v(T ) = v or r(T ) = r , 0 0 0 0 (3) the position of the object in a suitable coordinate frame, a(t) represents the acceleration of the object in the same and one of the constraints coordinate frame, and t is time. The velocity of the object (r ˙ ), given in the same coordinate frame as a(t)and r(t), is v(T ) = v or r(T ) = r , (4) 1 1 1 1 2 International Journal of Navigation and Observation find r(t)for T ≤ t ≤ T such that it satisfies how to exploit 1D finite element method (FEM) to solve the 0 1 underlying BVP with various possible choices of boundary r ¨(t) = a(t), conditions. At this stage, we assume that no additional constraints are given and that the measurements are exact. In v(T ) = v ,or r(T ) = r 0 0 0 0,(5) Section 5, we face the reality with faulty measurements and exploit linear additional constraints to enhance the accuracy v(T ) = v or r(T ) = r . 1 1 1 1 of the results. The underlying BVP is treated as exact, but We know for a fact that a problem of the form (5)has a the measurements are corrected using a linear sensor error model. As a result, we get two systems of linear equations: unique solution whenever the position is fixed at least at one boundary: This knowledge comes from the fact that (7)isa an exact one for the BVP and possibly an overdetermined one-dimensional case of a certain class of partial differential one for the parameters of the sensor error model. We will equations called Poisson problems, properties of which are solve these together to yield the corrected results. Finally, two well known [7]. real-life examples are discussed in Sections 6 and 7 before the Notice that there is a natural reason why problems (2) conclusions in Section 8. and (5) are distinguished. Namely, the solution methods for IVPs and BVPs differ substantially from each other 2. Preliminaries [6, 8]. This is also where this paper differs from numerous research articles considering inertial navigation: previously, Let us first motivate the chosen approach by means of the problems have been given in the form (2), whereas we examples: an application suitable to our approach is ski consider problems of the form (5). jumping, which has been the prime motivation of this study The key assumption of our approach is the following: [9]. As an inertial navigation problem, it includes the use inertial measurements related to a certain time period, includ- of consumer grade sensors with knowledge of the boundary ing a set of boundary values and possibly some additional values (in this case, position at both ends of the event). See constraints, are all available when processed. The length of the Section 7 for further details. time period (i.e., domain)[T , T ]in(5) can be anything 0 1 It turns out that especially sports applications tend to have from fractions of a second to some minutes. There are properties that are well suited to the considered approach: different ways we may end up with a BVP of the form the included actions are often periodic, in such a way that the (5): firstly, the underlying problem may naturally be a BVP. certain short-term action (e.g., a single step) repeats several Secondly, we may have an IVP, which we can pose as a BVP. times or some longer-term action a few times (e.g., a single This, of course, requires that we also know something about lap or a single jump as discussed earlier) during a certain the result at the end of the interesting event and that we event. In these kinds of applications, it is natural to encounter can afford to wait for the results until the end of the event. problems of the form (5) rather than (2): for example, in long In the first case, we generally do not have any additional jump (considering only the jump part), at time T (“take constraints we could use, but in the second case we end up off ”), the velocity of the shoe is known but position is not with a BVP and at least one additional constraint. We will and at time T (“landing”), the position of the shoe when it discuss examples of both cases in the next section. hits the surface of the sand can be accurately measured but The case where we only know the velocity at both ends the velocity is unknown. is, however, special and deserves some attention. According For a more general view, even a GPS-assisted inertial to the discussion above, we do not have a unique solution navigation system can be considered as a series of separate for this kind of a BVP although it fits the definition of our navigation periods with given boundary values rather than a model BVP (5). It is neither a valid IVP according to (2). As single event with additional constraints given at certain time it turns out, it is possible to solve also this without loss of instances. This is an example of an IVP, which can be posed generality. This is based on the fact that we can always fix the as a BVP. position at one boundary without changing the “shape” of the solution. Then, we will come up with a well-defined BVP with an additional constraint. 2.1. Some Remarks. There are two main classes of numerical The overall scope of this paper is to show how to treat methods for BVP’s. One class includes so-called shooting inertial navigation problems that are naturally (or know- methods and the other class methods of weighted residuals ingly) posed as a BVP of the form (5). Attitude computations such as FEM [8, 10]. We concentrate on the latter because are not considered, and where necessary, attitude is assumed of its property to minimize an error norm over the whole to be available. The form of problem (5)allowsustoconsider integration interval rather than minimizing only the local it as a set of three one-dimensional (1D) problems, rather error [7]. Another tempting property is that it concerns the than one three-dimensional (3D) problem. Notice that two position directly, giving us a possibility to more easily handle boundary values per a dimension are required to obtain a various types of additional constraints we will encounter unique solution. On the other hand, there is no need for the later on. boundary values for different dimensions to be of the same In many applications, inertial measurements are not the type. only source of information. In practice, however, the number This paper is organized as follows: in Section 2,pre- or the type of these additional constraints—combined with liminaries are discussed. In Sections 3 and 4, we will show the different kinds of a solution method—does not suggest International Journal of Navigation and Observation 3 the use of the traditional filtering (see Section 5 for details). [2, 15], whereas the proposed approach is a single- For these situations, we will introduce a computationally phase method with only few computations per an cheap and easily exploitable method to enhance the accuracy unknown. of the position and velocity estimates. The proposed method is based on sensor error modeling and is characterized by the Due to the significant differences in these two approach- following assumptions. es, we will in this context concentrate on the proposed method. Obviously, there are situations where the two (i) Sensor errors are modeled as constant errors. While methods could both be used. Comparison of the methods the behavior of a certain sensor error is in real life a in such a situation is interesting, although not addressed in stochastic process, one is usually able to fairly model this paper. it at least momentarily as a constant error. A suitable Finally, recall that problem (5) can be considered as a mathematical tool to characterize this is the Allan set of three 1D problems. Thus, let us focus on the 1D variance [11]. case for a while. Details of the more prevalent 3D case are considered later on. Notice that because of this, the symbols (ii) Considering consumer grade sensors, causes of the will be changed a bit: r , v ∈ R whereas r , v ∈ R and 0 0 0 0 most significant errors are usually known (e.g., bias so on. In the following treatment, it is assumed that the and scale factor error, both changing from turn-on to accelerometer samples represent an instantaneous value of turn-on [2]). the specific force. In other words, it is assumed that the sensor (iii) In particular, the sensor noise is not modeled. This output is not processed in any way before the “navigation is a conscious modeling decision to prevent unneces- computer.” sary “smoothing” of data. (iv) Additional constraints are treated as “exact”. That 3. Solution of a BVP is, the overall error in the additional constraints is assumed to be smaller than the error caused by the The main goal of this section is to form a linear system of simplification of the navigation equations. equations of the form Ax = b for the position estimates x. These are now expressed as a vector containing the position When additional constraints are available, problems at each discrete time instant (referred to as x ), where the resembling the ones considered here have previously been vector b is a function of a(t). In the following treatment, the resolved using the means of fixed interval smoothing (or total number of the samples is N , t ∈ [T , T ] is the value of i 0 1 “Kalman smoother”, if the type of the filter is fixed) [12– time instant i ∈ Z (1 ≤ i ≤ N ), and h = t − t . i i+1 i 15]. Considering inertial navigation, the most used methods Now, let us rephrase problem (5) as a variational of solving the problems are the two-filter smoother [16] equation and Rauch-Tung-Striebel smoother [17], used for example in [15, 18–20]. In both cases, the problem itself is posed as an T T 1 1 IVP and the basic idea is to run the filter in the forward ru ¨ dt = au dt ∀u ∈ U , (6) T T 0 0 direction as a “predictor” in phase one and then to run the filter in the backward direction while combining these two where U is a Hilbert space [7]. Second derivative of r can be results to yield the corrected result in phase two. Between the eliminated by integrating the left-hand side of (6) by parts, proposed method and the fixed interval smoothing method, which yields the fundamental difference is that in the proposed method the dynamics model is based on the BVP formulation and T 1 T T 1 1 1 in the previous approaches, on the IVP formulation with ru ¨ dt = ru ˙ − r˙u ˙ dt = au dt ∀u ∈ U , additional smoothing. T T T 0 T 0 0 Comparing the fixed interval smoothing technique to the (7) proposed method, in addition to the previously mentioned points, the most significant differences are as follows. where the notation “ ” stands for substitution. By evaluating (i) As fixed interval smoothing is run in both directions, the substitution term and rearranging (7), we get the filter needs two process models, which can be problematic [21]. As the proposed method is based T T 1 1 on the BVP formulation, it does not make a distinc- r˙u ˙ dt =− au dt + r˙(T )u(T ) 1 1 tion between forward and backward directions. T T 0 0 (8) (ii) In the filtering approach, each additional constraint − r˙(T )u(T ) ∀u ∈ U. 0 0 is assumed to be attached to a single time instance [2], while the proposed method does not make such While (8) is otherwise in a convenient form for our a restriction (see Section 5 for details). purposes, it is not well suited for practical computations, (iii) Fixed interval smoothing is a two-phase method because U is an infinite-dimensional space. Thus, let us requiring numerous computations per time step approximate U with a finite-dimensional space spanned by 4 International Journal of Navigation and Observation 1/h i−1 1/h i+2 φ φ φ φ 1 i i+1 N i−1 −1/h −1/h ··· ··· t i+1 h h h i−1 i i+1 t t t t t t t t 1 2 i−1 i i+1 i+2 N−1 N Figure 2: Time derivative of the lowest order basis functions φ Figure 1: Lowest order basis functions φ . (dashed line) and φ (dotted line). i+1 piecewise affine basis functions (“affine function = linear Now that we have discretized the problem, we are getting function + a constant”) closer to the equation Ax = b stated as our goal. In fact, 0, t< t , i−1 (12) is a system of linear equations. Thus, let us start by ⎪ assembling the matrix A. Knowing that the index i in (12) (t − t ) i−1 ⎨ , t ≤ t< t , i−1 i refers to a certain column and j to a certain row, the elements i−1 φ = (9) of A are given as ⎪ (t − t ) 1 − , t ≤ t< t , ⎪ i i+1 ⎪ h i T1 ˙ ˙ (13) A i, j = φ φ dt ∀i, j = 1, 2,... , N , 0, t ≥ t , i j i+1 often referred to as the “hat” functions. Functions (9)for few which are easy to compute by substitution of (11) into (13). values of i are shown in Figure 1. Also, let us use the same In practice, A is going to have only few nonzero elements basis functions φ to discretize u and r . This choice is often all of them at the diagonal, subdiagonal, and superdiagonal referred to as the Galerkin method [7]. Note that it is also (assuming the “obvious” indexing). For equally spaced nodes possible to choose basis functions φ different from the lowest with step size h, for example, we have elements 2/h at the order approximation used here, when considered necessary. diagonal and −1/h at the sub-and superdiagonal. Given the basis functions φ , notice that the position r Our next task is to compute the vector b.Asseenfrom can be approximated as a piecewise affine function (12), the task is to integrate term aφ over the interval [T , T ]. To do this, we fit a piecewise affine function to a(t) 0 1 r ≈ x φ . (10) i i between every node as seen in Figure 3 with dashed line, i=1 which yields From (8), we see that the derivatives of the basis functions j +1 j−1 (9) are also needed. It holds that aφ dt = a t +2a t j j−1 j ⎧ t j−1 (14) 0, t< t , ⎪ i−1 ⎪ 1 + 2a t + a t ⎪ j j +1 , t <t <t , i−1 i 6 ˙ h i−1 φ = (11) ⎪ for every j ∈ [2, N − 1]. For nodes j = 1and j = N,weget ⎪ − , t <t <t , i i+1 0, t> t . h i+1 1 aφ dt = [2a(t ) + a(t )], j 0 1 From Figure 2, one can verify that function φ is a piecewise i (15) constant function, discontinuous at points t , t ,and t . h i−1 i i+1 N−1 aφ dt = [a(t ) +2a(t )], j N−1 N From (8) we see that we need to integrate a similar term with N−1 discontinuities at the nodes over the domain. In other words, respectively. these discontinuities do not matter, which is a well-known Finally, let us consider how to apply the different types fact from integral theory. of boundary values into (12). At first, notice that terms In total, we are now in the position to discretize equation φ (T )and φ (T )seenin(12) will be nonzero (evaluating (8), which yields j 1 j 0 to value one) only for the values j = 1and j = N , T T 1 1 respectively. If r˙(T ) = v or r˙(T ) = v of (5)isgiven, 0 0 1 1 ˙ ˙ x φ φ dt =− aφ dt + r˙(T )φ (T ) i i j j 1 j 1 the respective boundary condition is called a Neumann T T 0 0 i=1 (12) boundary condition [7]. These can be applied directly by ˙( ) ( ) adding the given values into b. − r T φ T ∀j = 1, 2,... , N. 0 j 0 International Journal of Navigation and Observation 5 a(t ) i−1 An element-wise representation of C(t)is ⎡ ⎤ c (t) c (t) c (t) 11 12 13 ⎢ ⎥ ⎢ ⎥ ( ) ( ) ( ) C(t) = ⎢ c t c t c t ⎥ . (18) a(t ) 21 22 23 i+1 ⎣ ⎦ c (t) c (t) c (t) 32 32 33 a(t ) In this paper, we will assume that an estimate of C(t)for each time instant is available. Note that C(t) depends only on the attitude, which is independent of the position and velocity, when the assumptions considered in the first section ··· ··· t n are valid [2]. Furthermore, let vector g be the acceleration h h i−1 i due to the gravity represented in the navigation frame. With t t t t t t t 1 i−2 i−1 i i+1 i+2 N these notions, the specific force measurements f made by the Figure 3: Basis function φ along with the corresponding accelera- accelerometers can be “converted” into accelerations as tion values. n b b b n a = c f + c f + c f + g , 1 11 1 12 2 13 3 1 n b b b n a = c f + c f + c f + g , (19) 2 21 1 22 2 23 3 2 The other possibility for the boundary values is to have n b b b n r (T ) = r or r (T ) = r fixed thus having a Dirichlet 0 0 1 1 a = c f + c f + c f + g , 31 32 33 3 1 2 3 3 boundary condition [7]. As the value at the corresponding boundary is already known, it does not need to be solved. where the time dependencies have not been explicitly stated. Now, following from (19) and the right-hand side of Thus, we can move all terms depending on it to b reducing the number of unknown terms by one. For the reduced (12), we have an equation system, the corresponding term of the last two terms of (12) will be zero, as mentioned above. Recall that it is also n b j =− a φ dt i j possible to have a Dirichlet condition on one boundary and T (20) a Neumann condition on the other boundary. b b b n We have now means to assembly an equation of the form =− c f + c f + c f + g φ dt i1 i2 i3 j 1 2 3 i Ax = b (16) for the j th component of vector b ∀ i = [1,2,3]. Thus, the linear equation for the 3D position is ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ for the (1D) position of the object. Depending on the type of A 00 x b 1 1 1 the applied boundary conditions, the number of unknowns ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ n×n ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ n is equal to N − 1or N − 2. Matrix A ∈ R is known to be 0A 0 x b ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ . (21) 2 2 2 ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ symmetric and positive definite [22]. With the chosen basis 00 A x b 3 3 3 functions, A is also a tridiagonal matrix. In practice, these properties guarantee that (16) can be solved with linear time As discussed in the previous section, matrix A is a tridi- complexity [8]. In other words, when doubling the amount agonal matrix. Then, also the system matrix of (21)isnot of unknowns N , the time needed to solve the system is also only tridiagonal, but a block diagonal matrix. For instance, doubled (approximately), which is certainly not true for a if boundary conditions are fixed in the same way for all general system of linear equations. dimensions, it follows that A = A = A holds. 1 2 3 From the position data x,itisnow astraightforward With the equations derived in this section, it is possible task to compute the velocity using a suitable numerical to solve the 3D BVP using the presented FEM method. The differentiation formula. Since the position data is relatively solution of this more general problem can also be found with smooth due to the “double integration” process, we have not linear time complexity, as in the 1D case. experienced any problems in computing the derivative with an adequate accuracy. 5. Using a Number of Additional Constraints to Model Sensor Errors 4. Generalization to 3D So far we have constructed a method to compute position In this section, we will generalize the method described in and velocity data with certain boundary conditions. We will the previous section to the 3D case. For this, we introduce a now propose a way to exploit a number of additional con- coordinate transformation matrix C(t) used to transform the straints concerning the velocity or the position information b n specific force measurements f (t) into the accelerations a (t) to estimate sensor errors. For this, let us now precisely define represented in a suitable navigation frameasfollows: the term “additional constraint”: an additional constraint is a linear equation that bounds the position or the velocity of the n b a (t) = C(t)a (t). (17) object at an arbitrary number of time instances t ∈ [T , T ]. i 0 1 6 International Journal of Navigation and Observation At simplest, the previous definition could simply mean that they should be treated as unknowns. Other typical errors are v (t ) is fixed. On the other hand, a less intuitive equation caused by misalignment of the axes, changing temperature 1 i and drifting bias. [r (t )+ r (t )] = z is also a valid constraint. This way, 1 i 2 i 1 i=1 As a practical example used in the example 1 later in we can for example exploit constraints like x = x for any i i j Section 6, let us model the sensor errors as follows: and j without saying anything about the absolute position at these points. It could be hard to exploit these kind of f (t) = sf (t) + b, (25) constraints properly in traditional filtering problems. In this section, we will use constraints systematically to enhance accuracy. To make the concepts presented in this where f (t) is the corrected specific force measurement, s is section clear, let us first present the equations for the 1D case. some constant scaling factor, and b is some constant bias term correcting the erroneous measurement. By replacing the measured f (t) treated in the previous section with the 5.1. Treating Additional Constraints. Let us first assume N×N that A ∈ R holds. That is, given a Dirichlet bound- corrected acceleration f (t), it is easy to form an equation of ary condition, we add an auxiliary equation to the system the form instead of removing the corresponding boundary value from b = Fl (26) the system. This is a necessary procedure when exploiting additional constraints. for vector b seen in (16). Matrix F is in this case an N × In general, all linear constraints can be represented in the 2 matrix, whose elements are formed by computing the form T T 1 1 integrals − aφ and − φ from (12). Vector l is simply j j T T 0 0 Dx = e, (22) sb [ ] . Let us now replace the right-hand side of (16)with p×N p where D ∈ R and e ∈ R hold, where p is the number of the constraint equations. For a practical example closely b + b = b + Fl, (27) 0 0 related to the example 1, consider that (5)issolvedwith Dirichlet boundary values. If the Neumann boundary values where b depends only on the given boundary conditions are also known and the samples are equally spaced, one can (which are treated as exact) and corresponding boundary construct (22)as values. This distinction is necessary, since one should ⎡ ⎤ not modify the given boundary values by applying error 3 2 1 − − 0 ··· 00 0 0 modeling on them. ⎢ ⎥ 2h h 2h ⎢ ⎥ D = (23) In general, a problem of linear constraints used to model ⎣ ⎦ 1 2 3 000 0 ··· 0 − linear sensor errors can be stated as 2h h 2h Ax = b + Fl, ⎡ ⎤ (28) ⎣ ⎦ Dx = e, e = . (24) N×q q where F ∈ R , l ∈ R hold, and q is the number of Matrix D is formed using a three-point differentiation modeled sensor error terms. method based on polynomial interpolation [10]todeter- Now, since A is known to be invertible, one gets an mine Neumann boundary values from the displacement equation data. −1 −1 These constraints could be exploited just by computing DA Fl = e − DA b (29) the linear least squares problem constructed by adding these −1 additional equations to (16). Usually, however, it is favorable for l. Note that the matrix DA F has dimension p × q, that to use the additional constraints to model sensor errors,at is, it is typically a very small matrix compared to A.With l least when one has any knowledge of the types of the errors known, it is easy to solve for x. included in the measurements. For more information on In the case where p> q (more constraints than the assumptions related to the sensor error modeling, recall model parameters), one can also take reliability of different Section 2.1. In this text, we will present a linear sensor error measurements into account by solving a weighted least model. squares problem (WLS) [22]. In general, (29) has a unique solution l only when p = q holds and the row rank of D is full. Typically, making sure that p ≥ q and that the row 5.2. Sensor Error Model. In typical situations, bias offset rank of D is at least q,(29) has either a unique or a least (i.e., the sensor shows nonzero output when no forces are squares solution for l. In each case, the upper equation of acting upon it) and scaling factor are the two terms which (28) is treated as an exact equation. have to beknown very accurately in order to have any realistic position or velocity estimate as a solution. Because of the run-to-run variations of these errors, they cannot be 5.3. Additional Constraints in 3D. Let us now consider assumed to be constant between two separate events. Thus, the use of additional constraints in a 3D case. Using the International Journal of Navigation and Observation 7 notions form the previous section and (22) as an example, constraints can be represented as ⎡ ⎤ a true ⎢ ⎥ ⎢ ⎥ D⎢ x ⎥ = e. (30) meas ⎣ ⎦ ω(t) In addition to the sensor errors discussed in the 1D example, it is now also possible to model errors like attitude errors, unknown value of the acceleration due to the gravity, and cross-correlation of different sensors. A particularly useful method is to treat the acceleration due to the gravity Figure 4: Measurement setup of example 1. as an unknown three-dimensional vector, which is then subtracted from the specific force measurements given in the global coordinates. This reduces the systems sensitivity to where only the relevant indices are shown. Thus, we get size initial attitude errors. ⎡ ⎤ As an example, let us now derive the equations for sensor p p p p p p p0 0 1111 111 1211 121 1311 131 ⎢ ⎥ error of (25) for each axis in addition to the unknown accel- ⎢ ⎥ F = ⎢ p p p p p p 0p0⎥ , (34) 2111 211 2211 221 2311 231 eration due to the gravity discussed above. Thus, we have a ⎣ ⎦ total of 9 unknown model parameters. As one could expect, p p p p p p 00 p 3111 311 3211 321 3311 331 we must take the effects of the coordinate transformation matrix into account when deriving the necessary equations. n n n (35) l = s b s b s b g g g . 1 1 2 2 3 3 1 2 3 By plugging (25) into (21), we have For the nine unknown model parameters we need at least b b nine constraints. This can be covered, for example, with the b j =− c s f + b + c s f + b i i1 1 1 i2 2 2 1 2 knowledge of the velocity of the object at three different points. In total, we have equations identical to (28), with b n − c s f + b + g φ dt i3 3 3 j 3 i dimension N replaced by 3N . T T 1 1 =−s c f φ dt − b c φ dt 1 i1 j 1 i1 j 6. Example 1 T T 0 0 (31) T T To demonstrate the use of the proposed method, an example 1 1 − s c f φ dt − b c φ dt 2 i2 j 2 i2 j 2 using readily available consumer grade accelerometers [23] T T 0 0 was created. An accelerometer was mounted on a horizon- T T 1 1 tally rotating rate table, whose angular velocity (rotation − s c f φ dt − b c φ dt 3 i3 j 3 i3 j rate) can be controlled. The accelerometer was mounted in T T 0 0 such a way that its measurement direction was approximately the same as the direction of the tangential acceleration caused − g φ dt. by angular acceleration of the rate table, as seen in Figure 4. The aim of this example was not to get position and The several integrals in (31) are scalars for each j ∈ velocity as accurately as possible, but to compare different [1, 2,...N ], easily evaluated with (14)and (15), since no solution methods in a situation where one needs to get unknown terms appear inside the integrals. Let us now reasonable position and velocity estimates regardless of the gather the results into vectors fact that the measurement contains significant errors. The angular velocity of the rate table was set to follow function m b illustrated in the Figure 5 a certain number of times. With p j =− c f φ dt (32) klmn j kl l this kind of a setup, the true acceleration, velocity, and displacement of the mounted sensor were known up to the for all j ∈ [1, 2,... , N ]. With these notions, we get, for accuracy of the motor rotating the rate table. The errors example caused by the rate table were observed using the angular rate output of the motor and found to be negligible compared to p j =− c f φ dt, other error sources. kl11 kl j In first test, the rate table was set to repeat exactly the T1 function representedinFigure 5 ten times, leading to seven (33) p j =− c φ dt, kl1 kl j full revolutions (or 10.4meters) in T = 20 seconds. In the other test, the function of the same form was scaled in such T1 a way that 50 repeats lead to total of 62 full revolutions (or p j =− φ dt, 92.3meters) in T = 100 seconds. The sampling frequency of 0 8 International Journal of Navigation and Observation Table 1: Comparison of different methods computing object’s velocity and position. “I“ stands for time-stepping, “II“ for FEM, and “III“ for FEM with sensor error modeling. Test 1 T est 2 Vel. error [m/s] t t t t t t 1 2 3 1 2 3 I0.07 0.10 0.12 33.24 65.47 97.17 II 0.01 0.04 0.06 −0.67 0.50 1.15 III 0.04 0.05 0.04 0.33 0.54 0.44 Pos. error [m] t t t t t t 1 2 3 1 2 3 I0.10 0.49 1.00 510.96 1992.99 4434.70 II −0.31 −0.29 −0.15 −42.34 −43.42 −16.59 III −0.10 −0.07 −0.05 −6.78 −4.72 −0.61 9 10 −5 −10 11.52 2.53 00.51 1.52 Time (s) Time (s) Figure 6: Ideal (thick line) and measured acceleration (thin line). Figure 5: Angular velocity of the rate table. the sensor was set to 1000 Hz and an accurate (16 bit) analog three separate points by computing the difference between to digital converter was used. Accelerations were measured the computed and the real value. Point t was located at T/3, using a consumer grade 12 g accelerometer [23]. t at T/2, and t at 2T/3. 2 3 The raw accelerometer data was mapped to accelerations As seen in Table 1, method I seems to increase the error with the scale factor given by the manufacturer. From this with increasing time, as expected. Method II on the contrary, acceleration data, position and velocity were computed by a thanks to the basic property of the variational technique, number of different methods: does not increase the error but distributes it over the whole time period. The difference between these two methods is (i) “Traditional” IVP (double integration with trapezoid clearly seen in the velocity and position errors of the longer rule), test (test 2), where method II gives much better estimates (ii) FEM with Dirichlet boundary conditions (BCs), than method I. In each test, method III clearly outperforms (iii) FEM with Dirichlet BCs combined with additional methods I and II, which is expected due to the provided two Neumann BCs to supply the error model (25). additional constraints. Figures 7 and 8 demonstrate the differences between In methods I and II, the bias error of the sensor was velocity and position estimates given by methods I and II estimated by averaging the output while the sensor was at during test 1. Estimates given by method III coincide with rest. This was done in order to make method I (and to the reference plots. Figure 7 shows only the last spin of the some extent, method II) comparable to the method III by test 1 for better view of the differences. reduction of the large bias error. This is rarely possible in general and only method III can be used to reliably detect any remaining bias (example 2 in Section 7 is an example of 7. Example 2 this case). As a reference, the ideal (ideal driving motor) and the measured accelerations (accelerometer bias removed) of This example considers the computation of the velocity and the first spin are plotted in Figure 6. position of a ski jumper during a single jump. As compared Table 1 shows the main results of the tests. In each test, to the previous example, this is a more realistic and general the computed results were compared to the known value in inertial navigation problem with six degrees of freedom. ω (rad/s) Acceleration (m/s ) International Journal of Navigation and Observation 9 2.5 40 1.5 −20 0.5 −40 0 −60 −0.5 −80 19 19.520 20.5 −50 0 50 100 Time (s) x (m) Figure 7: Velocity estimates given by methods I (thin black line) Figure 9: Computed two-dimensional trajectories of two indepen- and II (thin gray line) compared to the ideal velocity (thick line). dent events (thin gray and thin black lines) along with the known profile of the hill (thick black line). take off and eventually land at some point of the landing hill. Drawn circles at both ends of the trajectories demonstrate the start and the end of the navigation period and the given Dirichlet boundary conditions. Drawn squares represent the landing points, in this case quite different for the two events, agreeing well with the recorded jump lengths. Notice that the trajectories are computed using two distinct measurement systems, each containing unique error sources. Figure 10 shows that the vertical component of the velocity of the same two jumps is plotted as a function of horizontal displacement. Notice how the absolute values −2 of the vertical velocity substantially differ while the small 0 5 10 15 20 25 changes in the velocity are practically identical. One might Time (s) first consider the small changes as typical stochastic errors caused by the inertial navigation system, especially when Figure 8: Position estimates given by methods I (thin black line) and II (thin gray line) compared to the ideal position (thick solid dealing with low performance sensors. line). Given that two independent measurement systems show the same variations in the velocity at the same locations, it is evident that there is actually only a negligible amount of stochastic errors present. Instead, the small variations are Computation of the attitude of the object is based on the data caused by deterministic sources, namely, in this particular given by three standard consumer grade gyroscopes [24]. application the uneven inrun hill. The position and velocity were then measured with three Unfortunately, the estimates cannot be compared with standard consumer grade accelerometers [23]and computed a reference trajectory, because such data are not available. with the proposed method with a sensor error model similar Thus,wecannotgivethe exactamountoferror presentinthe to the one presented in Section 5.3. The needed additional position and velocity estimates. We do however claim that the constraints contained information about achieved accuracy is something one does not typically expect (i) the location of the jumper at five points evenly spread from consumer grade inertial sensors. on the inrun hill (in Figure 9, the part of thick black line with negative x-and positive y-coordinates) 8. Conclusion (ii) the trajectory of the jumper after the landing, which The work was motivated by applications, where it is natural should coincide with the linearization of the landing to encounter BVPs instead of IVPs. In many cases, it is also hill (in Figure 9, the part of thick black line with x- possible to formulate an IVP as a BVP, given that the results coordinates greater than 100 m). are not required in real time. In Figure 9, two-dimensional trajectories of two jumpers Finite element method is utilized to solve inertial navi- are plotted along with the known profile of the hill. At first, gation problems formulated as BVPs. As a result, we get a the trajectories follow the profile of the hill until the jumpers linear system of equations for the position estimates, whose Velocity (m/s) Distance (m) y (m) 10 International Journal of Navigation and Observation −4 [4] M.S.Grewal, L. R. Weill, andA.P.Andrews, Global Positioning Systems, Inertial Navigation, and Integration, John Wiley & Sons, New York, NY, USA, 2001. [5] J.A.Farrell andM.Barth, The Global Positioning System & −6 Inertial Navigation, McGraw–Hill, New York, NY, USA, 1999. [6] L. F. Shampine, I. Gladwell, and S. Thompson, Solving ODEs with Matlab, Cambridge University Press, Cambridge, UK, −8 [7] S. Larsson and V. Thome, Partial Differential Equations with Numerical Methods, Springer, New York, NY, USA, 2003. [8] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetter- −10 ling, Numerical Recipes in C (The Art of Scientific Computing), Cambridge University Press, Cambridge, UK, 2nd edition, −12 [9] M. Virmavirta, T. Nieminen, T. Tarhasaari, and L. Kettunen, −70 −60 −50 −40 −30 −20 −10 0 “High precision inertial measurement for tracking the trajec- x (m) tories of ski jumpers “Smart Boot Project”,” in Proceedings of the 13th Congress European College of Sport Science, p. 572, Figure 10: Vertical velocity (v ) as a function of horizontal dis- Estoril, Portugal, July 2008. placement (x) of the two independent measurements. [10] D. Kincaid and W. Cheney, Numerical Analysis: Mathematics of Scientific Computing, BROOKS/COLE, Florence, Ky, USA, 3rd edition, 2002. [11] N. El-Sheimy, H. Haiying, and N. Xiaoji, “Analysis and mod- solution can be found with linear time complexity. It is eling of inertial sensors using allan variance,” EEE Transactions demonstrated that solving a BVP rather than an “equivalent” on Instrumentation and Measurement, vol. 57, no. 1, pp. 140– 149, 2008. IVP gives more accurate results. [12] A. Gelb, Applied Optimal Estimation, MIT Press, Cambridge, For further accuracy enhancements, an efficient way of Mass, USA, 1974. combining inertial measurements with possible additional [13] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with constraints is created. This gives us a possibility to model Applications to Tracking and Navigation, Chapman & Hall / constant sensor errors, known to limit the achievable CRC, London, UK, 2004. accuracy of the system. While the error model significantly [14] J. L. Crassidis and J. L. Junkins, Optimal Estimation of Dynamic enhances the accuracy of the system, it is kept computation- Systems, Chapman & Hall / CRC, London, UK, 2004. ally simple and easily adoptable. [15] P. D. Groves, Principles of GNSS, Inertial, and Multisensor In practice, the accuracy improvements allow us to Integrated Navigation Systems, Artech House, Norwood, Mass, exploit inertial sensors of certain performance level in more USA, 2008. challenging applications. For this, it is necessary to see that [16] D. Fraser and J. Potter, “The optimum linear smoother as the concept of inertial navigation does not invariably imply combination of two optimum linear filters,” IEEE Transactions on Automatic Control, vol. 14, no. 4, pp. 387–390, 1969. an IVP, but a BVP as well. Then, the use of FEM will provide [17] H. E. Rauch, F. Tung, and C. T. Striebel, “Maximum likelihood an efficient way to compute position and velocity estimates estimates of linear dynamic systems,” AIAA Journal, vol. 3, no. not prone to the accumulation of errors. 8, pp. 1445–1450, 1965. In larger scale, the current paper serves as an introduc- [18] K. Gade, “Navlab, a generic simulation and postprocessing tion to the idea of formulating inertial navigation problems tool for navigation,” European Journal of Navigation, vol. 2, no. as BVPs. As a consequence, further studies are needed to 4, pp. 21–59, 2004. address problems to which the presented tools do not provide [19] Y. Yang, Z. Jin, W. Tian, and F. Qian, “Application of fixed an obvious solution. These include, for example, stochastic interval smoothing to gps/dr integrated navigation system,” in errors, reliability of the possible additional constraints (as Proceedings of the IEEE Intelligent Transportation Systems, vol. compared to the accuracy of the IMU), and coupling of 2, pp. 1027–1031, October 2003. position and attitude errors. [20] A. B. Willumsen and ∅. Hegrenæs, “The joys of smoothing,” in Proceedings of the IEEE Bremen: Balancing Technology with Future Needs (OCEANS ’09 ), pp. 1–7, Bremen, Germany, May References [21] M. Klaas, M. Briers, N. de Freitas, A. Doucet, S. Maskell, and D. Lang, “Fast particle smoothing: if i had a million particles,” [1] I. Y. Bar-Itzhack, “Navigation computation in terrestrial strapdown inertial navigation systems,” IEEE Transactions on in Proceedings of the 23rd International Conference on Machine Learning (ICML ’06), vol. 148, pp. 481–488, Pittsburgh, Pa, Aerospace and Electronic Systems, vol. 13, no. 6, pp. 679–689, 1977. USA, 2006. [22] J. W. Demmel, Applied Numerical Linear Algebra, SIAM, [2] D.H.Titterton andJ.L.Weston, Strapdown Inertial Navigation Technology, Institution of Engineering and Technology, Lon- Philadelphia, Pa, USA, 1st edition, 1997. [23] VTI. VTI SCA620-CHCV1A Datasheet, 2006. Revision 2/2, don, UK, 2nd edition, 2004. [3] A. B. Chatfield, Fundamentals of High Accuracy Inertial Nav- Checked on 22.12.2009. [24] Analog Devices. AD ADXRS300ABG Datasheet, 2004. Revision igation, American Institute of Aeronautics and Astronautics, Reston, Va, USA, 1997. B, Checked on 22.12.2009. v (m/s) y International Journal of Rotating Machinery International Journal of Journal of The Scientific Journal of Distributed Engineering World Journal Sensors Sensor Networks Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Volume 2014 Journal of Control Science and Engineering Advances in Civil Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Submit your manuscripts at http://www.hindawi.com Journal of Journal of Electrical and Computer Robotics Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 VLSI Design Advances in OptoElectronics International Journal of Modelling & Aerospace International Journal of Simulation Navigation and in Engineering Engineering Observation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2010 Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com http://www.hindawi.com Volume 2014 International Journal of Active and Passive International Journal of Antennas and Advances in Chemical Engineering Propagation Electronic Components Shock and Vibration Acoustics and Vibration Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

Journal

International Journal of Navigation and ObservationHindawi Publishing Corporation

Published: Jun 30, 2010

References