A Study of the Effects of Stochastic Inertial Sensor Errors in Dead-Reckoning Navigation Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classifled information. John H. Wall Certiflcate of Approval: George T. Flowers Professor Mechanical Engineering David M. Bevly, Chair Assistant Professor Mechanical Engineering John Y. Hung Professor Electrical and Computer Engineering Joe F. Pittman Interim Dean Graduate School A Study of the Effects of Stochastic Inertial Sensor Errors in Dead-Reckoning Navigation John H. Wall A Thesis Submitted to the Graduate Faculty of Auburn University in Partial Fulflllment of the Requirements for the Degree of Master of Science Auburn, Alabama August 4, 2007 A Study of the Effects of Stochastic Inertial Sensor Errors in Dead-Reckoning Navigation John H. Wall Permission is granted to Auburn University to make copies of this thesis at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. Signature of Author Date of Graduation iii Thesis Abstract A Study of the Effects of Stochastic Inertial Sensor Errors in Dead-Reckoning Navigation John H. Wall Master of Science, August 4, 2007 (B.S.M.E., Christian Brothers University, 2005) 144 Typed Pages Directed by David M. Bevly The research presented in this thesis seeks to quantify the error growth of navigation frame attitude, velocity, and position as solely derived from acceleration and rotation- rate measurements from a strapdown Inertial Measurement Unit (IMU). The wide-spread availability of the Global Positioning System (GPS) and increased technological advances in Inertial Navigation Systems (INS) technology has made possible the use of increasingly afiordable and compact GPS/INS navigation systems. While the fusion of GPS and inertial sensing technology ofiers exceptional performance under nominal conditions, the accuracy of the provided solution degrades rapidly when traveling under bridges, dense foliage, or in urban canyons due to loss of communication with GPS satellites. The degradation of the navigation solution in this inertial dead-reckoning mode is a direct result of the numerical integration of stochastic errors exhibited by the inertial sensors themselves. As the accuracy of the GPS/INS combined system depends heavily on the standalone performance of the INS, flrm quantiflcation of the performance of inertial dead-reckoning is imperative for system selection and design. iv To provide quantiflcation of the accuracy of inertial dead-reckoning, stochastic mod- els are selected which approximate the noise and bias drift present on a wide variety of both accelerometers and rate-gyroscopes. The stochastic identiflcation techniques of Al- lan variance and experimental autocorrelation are presented to illustrate the extraction of process parameters from experimental data using the assumed model forms. The selected models are then used to develop analytical expressions for the variance of subse- quent integrations of the stochastic error processes. The resulting analytical expressions are validated using Monte Carlo simulations. The analytical analysis is extended to a simple navigation scenario in which a vehicle is constrained to travel on a planar surface with no lateral velocity. Monte Carlo simulation techniques are employed to exemplify and compare the expected results of inertial navigation in higher dynamic scenarios. v Acknowledgments I would like to thank my thesis advisor, Dr. David Bevly, for providing me with the invaluable opportunity to work as a graduate research assistant in the GPS and Vehicle Dynamics Lab (GAVLab) at Auburn University. His ever enthusiastic vigor for new research and his comfortable attitude made my experience as a master?s student most rewarding. I also give thanks to Dr. George Flowers for his warm welcome to Auburn and support throughout my time in the program. I owe much of my success as a student and researcher to all the students in the GAVLab. My sincere gratitude is due to my fellow labmates for their helping hands in the theoretical and experimental and the willingness with which they shared their expertise. I thank my brother, Michael, for his multi-dimensional support as a roommate, friend, and technical advisor. I thank my parents, who have given their unconditional love and support. I thank Julia, for her earnest ear and ardent companionship that has kept my spirit high throughout this entire experience. vi Style manual or journal used Journal of Approximation Theory (together with the style known as \aums"). Bibliography follows the IEEE Transactions format. Computer software used The document preparation package TEX (speciflcally LATEX) together with the departmental style-flle aums.sty. vii Table of Contents List of Figures x 1 Introduction 1 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Prior Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 Inertial Sensor Modeling 8 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Simple Sensor Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Stochastic Sensor Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.1 Gaussian White Noise . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.2 The Gauss-Markov Process . . . . . . . . . . . . . . . . . . . . . . 12 2.3.3 A Simple Stochastic Model . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Stochastic Identiflcation Techniques . . . . . . . . . . . . . . . . . . . . . 16 2.4.1 The Allan Variance . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.2 Experimental Autocorrelation . . . . . . . . . . . . . . . . . . . . . 21 2.4.3 Implementation Issues . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.5 Experimental Quantiflcation and Identiflcation . . . . . . . . . . . . . . . 24 2.5.1 Manufacturer Sensor Speciflcations . . . . . . . . . . . . . . . . . . 24 2.5.2 Experimental Approach . . . . . . . . . . . . . . . . . . . . . . . . 25 2.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3 Covariance Propagation of Stochastic Errors 28 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Simplifled Navigation Scenario: Single Axis . . . . . . . . . . . . . . . . . 29 3.3 General Characterization of Raw Sensor Measurement . . . . . . . . . . . 30 3.4 Characterization of Integrated Sensor Measurement . . . . . . . . . . . . . 31 3.5 Single Axis Stochastic Error Contributions . . . . . . . . . . . . . . . . . . 32 3.5.1 Variance of Integrated Wide Band noise . . . . . . . . . . . . . . . 33 3.5.2 Variance of Double Integrated Wide Band noise . . . . . . . . . . . 35 3.5.3 Variance of 1st order Gauss-Markov process . . . . . . . . . . . . . 36 3.5.4 Variance of Integrated 1st order Gauss-Markov process . . . . . . . 38 3.5.5 Variance of Double Integrated 1st order Gauss-Markov process . . 41 3.5.6 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.6 Validation of the Error Propagation . . . . . . . . . . . . . . . . . . . . . 46 3.6.1 Propagation of Wide-Band Noise Process . . . . . . . . . . . . . . 47 viii 3.6.2 Propagation of Gauss-Markov Process . . . . . . . . . . . . . . . . 49 3.7 Application Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.8 Illustration of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.8.1 Relative Magnitudes . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.8.2 Efiect of Time Constant . . . . . . . . . . . . . . . . . . . . . . . . 58 3.9 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4 Two Dimensional Error Propagation 61 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 Velocity Error in Navigation Frame Under No-Slip Planar Motion . . . . . 65 4.2.1 Mean and Variance of East Velocity . . . . . . . . . . . . . . . . . 67 4.2.2 Mean and Variance of North Velocity . . . . . . . . . . . . . . . . 70 4.2.3 Cross Covariance of North and East Velocity . . . . . . . . . . . . 73 4.2.4 Probabilistic Characterization of Velocity Errors . . . . . . . . . . 75 4.2.5 Validation of the Velocity Error Characterization . . . . . . . . . . 76 4.3 Position Error in No-Slip Planar Motion . . . . . . . . . . . . . . . . . . . 80 4.4 Propagation of Error in Planar Motion with Slip . . . . . . . . . . . . . . 85 4.4.1 Acceleration Error in Navigation Frame for Slip-Case . . . . . . . 85 4.5 Comparison of Slip and No-Slip Mechanizations . . . . . . . . . . . . . . . 89 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5 Six DOF Analysis 96 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.2 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.2.1 Orientation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.2.2 Translation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.3 Comparison to Planar Mechanization . . . . . . . . . . . . . . . . . . . . . 100 6 Conclusions 108 6.1 Overall Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.2 Di?culties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 Bibliography 114 Appendices 117 A Stochastic Parameter Identification with an Automotive-Grade IMU 118 B Demonstration of 6-DOF Mechanization of Experimental Inertial Measurements Taken at Talledega Superspeedway 125 ix List of Figures 2.1 Sample Plot of Wide Band Noise, 2 = 1 . . . . . . . . . . . . . . . . . . 11 2.2 Sample Plot of Gauss-Markov process, 2b = 1 , ? = 100 . . . . . . . . . 14 2.3 Allan Variance of Simulated Data 5 Hz 2rw = 1:2, 2b = 4, ? = 300 . . . 20 2.4 Sample Autocorrelation: 1.7x105 time units, fs = 5 Hz, 2b = 4, ? = 200sec . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.1 Standard Deviation of Wide-Band Noise Process: 10Hz, ! = 1, 2000 Monte Carlo iterations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Standard Deviation of Integrated Wide-Band Noise Process: 10Hz, ! = 1, 2000 Monte Carlo iterations . . . . . . . . . . . . . . . . . . . . . . . . 48 3.3 Standard Deviation of Double Integrated Wide-Band Noise Process: 10Hz, ! = 1, 2000 Monte Carlo iterations . . . . . . . . . . . . . . . . . 49 3.4 Standard Deviation of Gauss-Markov Process: 10Hz, b = 2, ? = 30, 2000 Monte Carlo iterations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5 Standard Deviation of Integrated Gauss-Markov Process: 10Hz, b = 2, ? = 200, 2000 Monte Carlo iterations . . . . . . . . . . . . . . . . . . . . 50 3.6 Standard Deviation of Double Integrated Gauss-Markov Process: 10Hz, b = 2, ? = 200, 2000 Monte Carlo iterations . . . . . . . . . . . . . . . . 51 3.7 Body Constrained to Travel in One Direction . . . . . . . . . . . . . . . . 52 3.8 Simulated Acceleration Proflle: 10Hz, ! = 0.5, b = 0.25, ? = 200 . . . . 54 3.9 Velocity with Bounds from Accel Proflle . . . . . . . . . . . . . . . . . . 55 3.10 Position with Bounds from Accel Proflle . . . . . . . . . . . . . . . . . . 55 3.11 Integrated Sensor 1- Bounds for Ratio from 0.1 to 1 . . . . . . . . . . . . 57 3.12 Integrated Sensor 1- Bounds for Ratio from 0.01 to 0.1 . . . . . . . . . . 57 x 3.13 Integrated Sensor 1- Bounds for Ratio from 0.001 to 0.1 . . . . . . . . . 58 3.14 Integrated Sensor 1- Bounds for Fixed Ratio of 0.1 and ? Between 200 to 1000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.15 Value of Integrated 1- Bounds at time index of 120 for Fixed Ratio of 0.1 and ? Between 200 to 1000 . . . . . . . . . . . . . . . . . . . . . . . . 59 4.1 Simplifled Coordinate Frame . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.2 Simplifled Coordinate Frame 2D . . . . . . . . . . . . . . . . . . . . . . . 63 4.3 Navigation Relationships for Side-Slip Vehicle . . . . . . . . . . . . . . . . 63 4.4 Navigation Relationships for Vehicle with No-Slip . . . . . . . . . . . . . . 64 4.5 Simulated Position Trajectory . . . . . . . . . . . . . . . . . . . . . . . . 77 4.6 Standard Deviation of East Velocity, 2000 Monte Carlo iterations . . . . 79 4.7 Standard Deviation of North Velocity, 2000 Monte Carlo iterations . . . . 79 4.8 Deflned Position Trajectory . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.9 Velocity Trajectory With No Side Slip . . . . . . . . . . . . . . . . . . . . 91 4.10 Yaw Angle Trajectory With No Side Slip . . . . . . . . . . . . . . . . . . 91 4.11 3- Bounds on Simulated Yaw Angle . . . . . . . . . . . . . . . . . . . . . 92 4.12 RMS 3- Bounds on Velocity . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.13 RMS 3- Bounds on Position . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.1 Simplifled Coordinate Frame . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.2 Mechanization of IMU Measurements . . . . . . . . . . . . . . . . . . . . . 100 5.3 Deflned Position Trajectory . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.4 Velocity Trajectory for Planar Motion . . . . . . . . . . . . . . . . . . . . 102 5.5 Yaw Angle Trajectory for Planar Motion . . . . . . . . . . . . . . . . . . . 102 5.6 3- Bounds on Yaw Angle Mechanization Comparison . . . . . . . . . . . 104 xi 5.7 RMS 3- Bounds on Velocity Mechanization Comparison . . . . . . . . . 104 5.8 RMS 3- Bounds on Position Mechanization Comparison . . . . . . . . . 105 5.9 3- Bounds on Pitch Angle DOF Comparison . . . . . . . . . . . . . . . . 106 5.10 3- Bounds on Roll Angle DOF Comparison . . . . . . . . . . . . . . . . . 106 A.1 Filtered Accelerometer Outputs . . . . . . . . . . . . . . . . . . . . . . . . 119 A.2 Filtered Gyro Outputs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 A.3 Raw and Filtered Data Within Selected Section . . . . . . . . . . . . . . . 120 A.4 Filtered Data Within Selected Section . . . . . . . . . . . . . . . . . . . . 121 A.5 Autocorrelation of Filtered Gyro Data . . . . . . . . . . . . . . . . . . . . 122 A.6 Allan Variance of Gyro Data . . . . . . . . . . . . . . . . . . . . . . . . . 123 B.1 Position of Track from Starflre GPS . . . . . . . . . . . . . . . . . . . . . 126 B.2 Velocity from GPS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 B.3 Yaw Angle (Heading) from IMU . . . . . . . . . . . . . . . . . . . . . . . 128 B.4 Roll and Pitch Angles from IMU . . . . . . . . . . . . . . . . . . . . . . . 129 B.5 North and East Component Velocities . . . . . . . . . . . . . . . . . . . . 130 B.6 North and East Component Positions . . . . . . . . . . . . . . . . . . . . 132 xii Chapter 1 Introduction 1.1 Overview With the wide-spread availability of the Global Positioning System (GPS) and con- tinued technological advances in Inertial Navigation Systems (INS), the fleld of navi- gation continues to exploit these ever expanding resources. One of the most popular methods of current-day navigation involves the synergy of these two systems and is commonly known as GPS/INS integration. The basic idea of the GPS/INS navigation approach is to utilize the complementary strengths and weaknesses of each system in order to navigate with accuracy superior to that of either component system?s stand alone performance. GPS provides precise global position, velocity, and heading infor- mation, but updates at a relatively slow rate, is subject to interference, and requires a clear view of the sky. The INS most commonly consists of an inertial measurement unit (IMU) which is typically attached rigidly to the frame of the navigating body. The IMU is a device in which three accelerometers and three rate-gyroscopes (referred hence- forth simply as rate-gyros) are orthogonally mounted in a sealed unit. Its purpose is to measure the body accelerations and rotation rates in the corresponding component directions and about the corresponding axes as deflned by its alignment on the vehicle. The inertial sensors which make up an IMU are available in a wide-variety of grades determined by the basic principle of operation, quality of materials, and integrity of methods/design. The accelerometers and rate-gyros, as with many electronic sensors, are corrupted by stochastic-type disturbances (noise) which manifest themselves on the 1 sampled output of the device as random variables. The characteristics of these random variables depend upon the principle of operation and integrity of materials/methods used to produce the device. To achieve navigation information from such an INS, the outputs must be numerically integrated and transformed to attain the desired attitude, velocity, and position information. Due to the stochastic errors on the sensors, the integral values exhibit an ever increasing variance. Therefore, in contrast to GPS, Inertial Navigation Systems provide high update information, but will digress with time without corrections. In summary, the GPS/INS combined system uses the high update inertial measurements of the INS to boost the slow rate of the position information from GPS. The most common GPS/INS approach employs an optimal estimation technique known as the Kalman Filter [1]. This Navigation Kalman Filter numerically blends the high-rate, short term accuracy of the INS with the long term accuracy of the GPS in an optimal way. Many variations and approaches of the GPS/INS Kalman Filtering method have been researched with the intent of providing the best possible navigation solution. The loosely-coupled GPS/INS uses the common measurements from a standard receiver and requires at least four satellites for normal operation. More advanced approaches such as tightly-coupled and deeply-integrated GPS/INS probe deeper within the GPS receiver and blend the INS with more of the available GPS information. Both of these methods intend to improve some of the weaknesses of the common GPS/INS system for more robust navigation performance [2]. In addition to the GPS/INS only approach, re- searchers have also developed and implemented navigation systems which include other aiding measurements such as vision [3, 4], odometry [5], and laser-scanners [6]. These ad- ditional measurements supplement the existing GPS/INS system to improve the overall navigation performance, especially during GPS outages. 2 For the case when GPS is unavailable, the INS and other sensors become the sole means of navigation, and the accuracy of the attitude, velocity and position informa- tion as derived from body-frame-only measurements degrades with time. This GPS-less navigation mode is known as dead-reckoning whereby the body in motion navigates only with information on board itself. The characteristic error growth in the dead-reckoning mode depends upon the amount of time since the GPS signal outage, the equations used in providing the navigation solution (kinematics of vehicle motion), and the integrity of the sensors used to dead-reckon. The particular focus of the research presented in this text is to perform an analysis of inertial-only navigation in which both GPS and other sensors are not used. This analysis provides a quantiflcation of the dead-reckoning accuracy of the GPS/INS system when GPS is unavailable. As there are many difierent types and grades of inertial measure- ment units, the navigation system designer is charged with the tasks of selecting an IMU appropriate for the environment in which it will be employed, the budget of the project, and speciflc objectives of the flnal system. The typical manufacturer speciflcations of the raw inertial measurement errors do not fully depict its performance as a navigation means. That is, they fail to provide the information necessary to flrmly quantify the nav- igation performance of an IMU within the common scenarios encountered. The common speciflcations from the sensor manufacturer provide only a rough quantiflcation of the raw output of the sensors; no information of the propagated error is available. There- fore, to accurately compare and distinguish the dead-reckoning performance of the many varieties of IMUs, an analysis of the integrated inertial measurements is imperative. 3 1.2 Prior Art GPS/INS integration is a very widely implemented and studied scheme as it provides high rate (100 Hz) and high accuracy positioning (1-10 meters) information to a navigat- ing vehicle at a reasonable cost [2]. As GPS is susceptible to signal loss, the performance of the GPS/INS system is often quantifled by its INS-only (i.e. dead-reckoning) perfor- mance. Approaches to the quantiflcation of the INS-only navigation performance have been based heavily on experimental results. Techniques such as parametric modeling of IMU-derived positioning data [7] and comparisons based on the use of vehicle dynamics in mechanization algorithms [8] all intend to quantify the performance of the GPS/INS system during GPS outages. These post-processing analyses, while giving precise infor- mation about the INS-only system performance for particular experimental cases, fail to provide information in support of the general performance of inertial dead-reckoning. In support of the more general approach to the quantiflcation of performance of iner- tial dead-reckoning, research has been conducted in which the analysis commences with initial investigation of the inertial sensor errors themselves. Such approaches have inves- tigated the characteristics of the stochastic behavior of the raw inertial sensor outputs by using experimental identiflcation techniques of Allan variance [9, 10, 11, 12, 13], and autocorrelation [11, 14]. Based on stochastic error models selected using the identifl- cation techniques, the research has been extended to provide the exact error growth of the integral values of the stochastic process [11]. Other researchers following the same path, have investigated the analytical propagation of error into the navigation solution when mechanized in simple scenarios [14, 15]. These propagation analyses, however, have restricted their analytical study to the the efiect of white (wide-band) noise on 4 the navigation solution thereby neglecting signiflcant efiects due to sensor bias-drift [14]. The research presented in this thesis has expanded the previous analytical work to in- clude analytical results of the in uence of both wide-band noise and sensor drift on the integral sensor outputs. Additional studies of the sensor errors in more complex vehicle motion and mechanization methods are presented as well. 1.3 Contributions To accomplish the goals of this research, this thesis presents a study of the stochastic inertial sensor error characteristics to provide a practical quantiflcation of the uncertainty growth when dead-reckoning with inertial measurements. This quantiflcation provides the a priori information on the performance of GPS/INS/AuxSensor navigation in sup- port of the system designer?s goal of achieving the best navigation performance within the criteria of given objectives. The research presented provides quantiflcation of inertial navigation through: ? Proposing a candidate sensor model su?cient to capture the stochastic behavior of inertial sensor outputs and detailing the means by which the model parameters may be identifled. ? Development of expressions for the variance of the integrated values of the inertial sensor error sources. ? An analysis on the relative in uence of sensor error model parameters on the variance of integrated outputs. ? A study of the propagation of inertial sensor errors into the navigation equations for a vehicle kinematically constrained in a two-dimensional scenario. 5 ? A comparison of the in uence of the navigation equations on the performance of inertial navigation under various vehicle motion assumptions. 1.4 Thesis Organization Chapter 2 introduces a simple inertial error model consisting of the sum of two stationary Gaussian random processes. The processes model the stochastic behavior of white noise and sensor bias drift observed on both accelerometers and rate-gyros. The means of experimental identiflcation of the speciflc error process parameters with the use of Allan variance and Autocorrelation techniques are presented. Chapter 3, the heart of this research, utilizes the stochastic models presented in Chapter 2 to derive analytical expressions for the variance of numerically integrated and double-integrated error processes. The derived expressions are validated using Monte Carlo computer simulations with flxed parameters. A simple example is used to illustrate the applicability of the integrated error variance expressions for a simple single-axis navigation scenario. This chapter concludes its contribution with a study of the efiects of the stochastic error process parameters on the integrated sensor error growth. Chapter 4 applies the information from Chapter 3 to study a simple navigation scenario in which a vehicle travels on a at plane with no lateral velocity (i.e no side- slip). Expressions for the variance of 2-D velocity error are derived in closed-form and the framework for 2-D position error is shown. The analysis is then extended to the more general 2-D case where a vehicle experiences side-slip and exhibits lateral velocity. A Monte Carlo computer simulation is used to show the additional error when no-slip trajectory is processed with the assumed-slip equations. 6 Chapter 5 introduces the six degrees-of-freedom (6-DOF) navigation equations as required for the inertial navigation of a body experiencing motion in the most general sense. The characteristics of the error for the 6-DOF inertial navigation scenario are exemplifled with the use a computer simulation. In the simulation, the inertial navigation error from the 6-DOF method is compared to the planar navigation method for a planar vehicle trajectory. The purpose of this chapter is to show the efiect of additional sensors and transformations when they are not required. Chapter 6 concludes the work of this thesis with a summary of the results, discus- sionofthevariousassumptionsandapproximationsemployedinthetext, andsuggestions for future research. 7 Chapter 2 Inertial Sensor Modeling 2.1 Introduction A preliminary step in describing the behavior of inertial navigation error is to deter- mine appropriate error models of the inertial sensors themselves. This chapter presents a simple inertial sensor model su?cient to describe the approximate stochastic nature of the outputs of both accelerometers and rate-gyros across the range of available devices. The identiflcation of the models of two random processes presented in this chapter, have been well researched by many authors including [9, 11, 16] and the applicable techniques and limitations of these methods are presented proceeding the models in the following sections. 2.2 Simple Sensor Models Simple stochastic sensor models are selected to provide approximate representations to the behavior of that observed on a wide range of inertial measurement sensors. These simple models apply to both accelerometers and rate-gyros and lend themselves well for deriving analytical results which quantify the propagation of sensor errors into their processed states. Although more advanced stochastic models could be developed, the simpler models provide approximations in support of the goals of this thesis: closed-form expressions for the error growth when dead-reckoning with inertial measurements. The general inertial measurement model used in this study has the form of Equation (2.1) and is suitable for either an accelerometer or rate-gyro. 8 ymeas = (SF)y +?+c (2.1) y is the true sensor value ? are stochastic ?static? terms c is the constant bias SF is the scale factor The scale factor, though shown in [17] to exhibit some stochastic behavior through- out the range of sensor motion and temperature, is assumed to have flxed relationship between the sensor?s input and output. The bias, c is also a flxed quantity, exhibits a de- terministic behavior with respect to the motion the sensor is sensing. Simple calibration techniques can easily ascertain the scale-factor and constant bias and are thus removable prior to implementation. The stochastic terms however, are not predictable and remain even on a calibrated sensor output. Therefore, a fully calibrated static sensor output as measured can be described by Equation (2.2). ymeas = ? (2.2) The stochastic terms, ?, are, at this point in the heuristic approach, unknown. The following section introduces simple stochastic processes which intend to capture the unpredictable stochastic behavior of the calibrated sensor output. 9 2.3 Stochastic Sensor Models 2.3.1 Gaussian White Noise Inspection of a short-duration static sensor output reveals that the value of the signal jumps from one value to the next in an unpredictable, or random, manner. Successive samples in time of such a signal are described as uncorrelated from one time-step to the next. A random process whose successive values in time are uncorrelated is known as a white noise process. While the conceptual abstraction of white noise implies that there is inflnite frequency content in the process, the sensor output is better described by wide- band noise. As as the frequency content of a signal is limited by sampling equipment, materials, and other phenomena, wide-band noise is, like white-noise, uncorrelated, but for all measurable or applicable frequencies. Closer observations of the short-duration static output indicate that the values of the static output at any given time are scattered in higher density about some mean value and becoming scarcer with increasing value from the mean. A good statistical model for this scattering of the data is the Gaussian or Normal Distribution, whose probability density function is given by Equation (2.3). fX(x) = 1p2? x exp " ?12 x?m x x ?2# (2.3) mx is the mean of random variable x x is the standard deviation of random variable x 10 Complete probabilistic characterization of a random variable, x, with a Gaussian distribution requires knowledge only of the mean value, mx, and variance 2x. Therefore, for the uncorrelated wide-band Gaussian noise process, the sensor output value at any timeisdescribedbyaGaussiandistribution withthesamemeanandvariance. Therefore, each successive time sample of the process is independent (i.e. no relationship exists between the values from one step to the next). An example of a zero mean, unit-variance Gaussian wide-band noise process is shown in Figure 2.1. 0 200 400 600 800 1000?5 0 5 Value Time Index output ?? +? Figure 2.1: Sample Plot of Wide Band Noise, 2 = 1 For high grade sensors, the static output for short time durations (order of minutes) can solely be described by a wide-band noise process. However as the cost of sensor decreases, the integrity of the its output also decreases and other stochastic phenomena can be observed on the measurement. The output of lower-grade sensors exhibit both uncorrelated (wide-band) noise and additional stochastic noises which exhibit varying levels of time correlation, discontinuity, and other irregularities. The discontinuities and irregularities pose a problem in stochastic characterization as their presence dominates in 11 less-expensive sensor outputs. As a result, a signiflcant amount of study has been dedi- cated to the development of sophisticated models that accurately capture such low-grade sensor behaviors. However, for the purpose of this research, a simple stochastic model is selected to provide a conservative approximation to the bias drift in accordance with the practical utility of the modeling approach of this thesis. The dominant phenomena on the low grade sensors often materialize as a slowly wandering error in the output and are often referred to as a sensor?s drifting bias, or bias drift. The following section presents an approximate model for this stochastic drifting behavior. 2.3.2 The Gauss-Markov Process The 1st-order Gauss-Markov process has been used extensively in the navigation and estimation community to model the various stochastic drift characteristics present on many types of navigation system outputs [11, 18]. For the purpose of the research goals of this thesis, this process provides a conservative approximation to the observed bias drift on many inertial measurements. A random process is said to be Markovian if its probability density function (PDF) at any point in future time can be completely specifled with the knowledge of the pro- cess PDF at the current time. The Markovian property is analogous to the concept in linear state-space systems whereby the future states can be ascertained by current states and inputs. A Gauss-Markov process is a stochastic process whose underlying random phenomenon that drives the process is as a Gaussian sequence of random variables [19]. A commonly used 1st-order Gauss-Markov process model is simply the output of a low-pass fllter with a zero-mean white noise input. The governing stochastic linear difierential equation for such a process is expressed as 12 _b = ?b ? +!b (2.4) ? is the time constant !b is a zero mean Gaussian random variable with variance, 2!b The process given by Equation (2.4), in continuous time, is somewhat of an abstrac- tion of the more applicable discrete version of the Gauss-Markov process (or sequence). This stochastic sequence is given by the linear stochastic difierence Equation (2.5). bk = 'kbk?1 +wbk 'k = e??tk? (2.5) ' is the state-transition matrix for the process wb is a zero mean Gaussian random variable with variance, 2wb As seen by Equation (2.5) the process output is the sum of the Gaussian driving noise, !, and past values of itself. The output of the process therefore can also be de- scribed by a Gaussian distribution and can be completely probabilistically characterized by its mean and variance functions. In contrast to wide-band noise, however, the Markov process exhibits a non-zero time correlation for any given realization of the process be- cause of its dependence on its past values. This correlation characteristic is what causes 13 the process to appear as a slowly drifting bias, which is the desired model. One such example of the Gauss-Markov process is shown in Figure 2.2. 0 200 400 600 800 1000?2 ?1 0 1 2 3 Value Time Index output ?? +? Figure 2.2: Sample Plot of Gauss-Markov process, 2b = 1 , ? = 100 Since the driving process is zero mean and Gaussian, the output of the process is zero mean and Gaussian with a transformed variance. Since the time constant is a value greater than the sampling frequency (generally much larger), the characteristic root of the 1st-order Markov dynamics is inside the unit circle, and thus the process is stable. The stability of the process indicates that the variance of its output will reach a steady state value after some initial settling time at which point the process is considered to be stationary. Once steady state has reached, the process autocorrelation function takes the form of Equation (2.6). Rbb(T) = b2e?T=? (2.6) 14 ? is the time constant or correlation time b is the variance of the process 2!k The simple autocorrelation expression lends itself well to experimental identiflcation as it is straightforward to extract the magnitude and time constant of the process from a given plot of (2.6). Due to the form of Equation (2.6), the Gauss-Markov process is also referred to as exponentially correlated noise [16]. 2.3.3 A Simple Stochastic Model For short time intervals an inertial sensor?s output appears as uncorrelated noise (white noise); however for longer time intervals the sensor exhibits a correlated noise (drifting bias). An approximate model of this behavior is simply the sum of the two random processes introduced in the previous sections. The error due to the stochastic behavior of an inertial sensor at time step, k, can thus be described by Equation (2.7). ?k = !k +bk (2.7) !k is uncorrelated wide-band noise with zero mean and variance 2! bk is a 1st-order Gauss-Markov process with time constant, ? and variance 2b The model above assumes that both processes are zero-mean and no correlation exists between the white noise process and the Gauss-Markov process; they are indepen- dent. The model consists of three parameters which allow for three degrees of freedom 15 in describing the sensor behavior of any given accelerometer or rate-gyro. The following section introduces methodologies used in determining the best set of parameters from experimental data of a given inertial sensor. 2.4 Stochastic Identiflcation Techniques Several techniques exist to process experimental sensor data and determine the various types of error sources present. Given that the raw sensor output is modeled by Equation (2.7), such identiflcation techniques allow the extraction of the stochastic model parameters !, 2b, and ?. Of the array of stochastic modeling techniques, the Allan variance technique has become widely popular in the inertial navigation community due to its intuitive implementation, straightforward interpretation, and in most cases unique indication of various stochastic disturbances [20]. This section shows its use as a comprehensive means to determine approximate parameters to describe a given inertial sensor. Following the Allan variance technique, the experimental autocorrelation method is presented. Supplemental to the Allan variance, the autocorrelation method is used for speciflc identiflcation of the parameters associated with the Gauss-Markov process, 2b, and ?. 2.4.1 The Allan Variance The Allan variance technique was introduced by David Allan in the 1960s to char- acterize the frequency stability of high-precision atomic clocks [21]. The Allan variance technique is sometimes analogized to the time domain version of the Fourier transform as it provides a measure of the dominance of a stochastic process over various ranges of time. The technique computes the variance of arrangements of successive averages of a 16 time set of data which indicates the contribution or dominance of various error sources as a function of averaging time. Inspection of plots of the Allan variance versus averag- ing time indicate the presence of speciflc stochastic processes present on a measurement. The square root of the Allan variance (root Allan variance) for many common stochastic processes appear as straight lines on a log-log scale. As a result, simple visual inspec- tion of the root Allan variance provides immediate information about a sensor?s overall stochastic behavior. The precise deflnition of the Allan variance is repeated below as adapted from [9]. Given a set of N inertial measurements, ?, sampled at a rate of fs Hz, deflne a vector of averaging times, T, ranging from T0 seconds to up to half the total time length of the data set ( N2fs) as shown in Equation (2.8). T = ? T0 T0 +fs T0 +2fs ::: N2fs ? (2.8) For each averaging time, T, deflne K = N=M clusters where M is the number of samples per cluster (M = Tfs). Compute cluster averages with Equation (2.9) where k is the time index of the raw data. ??(T) = 1 M MX i=1 ?(k?1)M+i; k = 1;:::;K (2.9) The Allan variance is then computed using Equation (2.10), with an approximation to the true ensemble average [22]. 2AV (T) = 12(K ?1) K?1X k=1 ??? k+1(M)? ??k(M) ?2 (2.10) 17 The Allan variance can be expressed as a simple sum of the Allan variance contri- butions of each of the dominant processes. For example, the Allan variance as a function the contribution of process 1, 2, 3 and so on is expressed as 2AV = 2p1 + 2p2 + 2p3 +::: (2.11) Well-known analytical expressions of the power spectral density (PSD) of common stochastic processes have been related to analytical expressions for the Allan variance [9, 22]. The result of this relationship has yielded Allan variance equations for each noise process in terms of its processes parameters and the averaging time, T. The component Allan variances can be summed using Equation (2.11) to yield the total Allan variance curve as a function of averaging time. A full description of the method used to compute the Allan variance from experimental data can be found in [20]. The inertial sensor model introduced in Section 1.3 is the sum of the two stochas- tic processes: wide-band noise and exponentially correlated noise (Gauss-Markov). The wide-bandnoiseisquantifledbythe\randomwalk"parameter, rw, whichisthestandard- deviation of the process normalized to the square-root of the sampling frequency, fs in Hertz. The Allan variance expression for wide-band noise is thus given by Equation (2.12). 2rw(T) = wp fs ?2 (2.12) 18 The Allan variance expression for a flrst-order Markov process is given in terms of its time constant, ?, and driving noise variance, 2!b, as shown in Equation (2.13). 2bias(T) = ( !b?) 2 T h 1? ?2T ? 3?4exp?T? +exp?2T? ?i (2.13) The resulting analytical expression for the Allan variance of a static sensor sensor output with a stochastic model as in Equation (2.7) is given as the sum of equations (2.12) and (2.13) as shown in equation (2.14). 2AV (T) = 2rwT + ( !b?) 2 T h 1? ?2T ? 3?4exp?T? +exp?2T? ?i (2.14) The root Allan variance as is most commonly plotted is simply the square root of the summed quantities in Equation (2.11). Figure 2.3 shows a sample root Allan variance plot as computed from a simulated sensor output with speciflcations as labeled. 19 100 101 102 103 104 10?1 100 Averaging time, T Allan Root Variance, ? A V(T) Data Equation Bound Figure 2.3: Allan Variance of Simulated Data 5 Hz 2rw = 1:2, 2b = 4, ? = 300 The analytical expressions shown above have been used by others to curve-flt ex- perimental Allan variance data for the purpose of extracting the underlying stochastic process parameters and thus characterizing the process. [9, 22, 13]. For example, the random walk parameter is straightforward to extract as it is simply the value of the Allan variance at the averaging time equal to one second. However, since the expression for the Markov process in Equation (2.13) is a non-linear function of the process parameters, estimation becomes di?cult and employment of manual methods may be necessary to provide an estimate of the process magnitude and time constant. One such problem in practical utilization of the Allan variance for sensor character- ization is that its accuracy is limited by the time-length of the experimental data set. The range of averaging time is computed for half the range of experimental data. Ad- ditionally, the number of segments of averaged data for which the variance is computed 20 become smaller as the averaging time increases and thus yields statistically less signifl- cant values. As a result, the bounds on the accuracy of the Allan variance increase with averaging time at rates dependent upon the total length of the experimental data set. This behavior can been seen in the sample Allan variance in Figure 2.3. Such bounds are discussed further in [9]. Due to this requirement, accurate Allan variance identiflcation of an inertial sensor exhibiting slow drift characteristics requires a sample data set of length many times greater the time constant of the particular drift of interest. For many inertial measurements, the times required for reasonable accuracy can be on the order of several days. Additionally, the time constant for the sensor is simply unknown since it is to be identifled. Some rules of thumb are suggested for the length and required sample frequency of the sensor data set to ensure accurate Allan variance identiflcation for time-varying processes in [20]. 2.4.2 Experimental Autocorrelation The Allan variance technique, while su?cient to extract the parameter associated with the white noise process, remains a di?cult method for the bias drift characteriza- tion. As the Gauss-Markov process in steady-state has an autocorrelation in the form of Equation (2.6), the bias time constant, ?, and magnitude of drift, b, can be extracted from the experimental autocorrelation. However, since the sensor is modeled as the sum of the two processes, successful identiflcation of the Markov process requires its isolation. Approximate isolation can be performed by low-pass flltering the raw sensor output. The flltering removes the higher-frequencies for which the output is uncorrelated while leav- ing the correlated low-frequency data of interest. However, despite the attenuation of the 21 higher frequencies of the wide-band noise process, its lower-frequencies still exist as ad- ditive noise on the drifting bias, increasing the di?culty of identiflcation. Furthermore, appropriate selection of the isolating fllter cut-ofi frequency requires some knowledge of the approximate time constant of the drift process. As it is the time constant that is to be identifled, the autocorrelation technique is an approximate and iterative process. Equation (2.6) gives the characteristic autocorrelation function of the Gauss-Markov process, which is a simple exponential decay with an initial magnitude of 2b and time constant, ?. Figure 2.4 shows a sample autocorrelation plot of a simulated Markov process and its analytical curve. As is shown by the dotted lines on the plot, the time constant is found by reading the time-shift for which the autocorrelation value decays to 1 e of its initial magnitude. The initial magnitude of the Gauss-Markov process is simply the y-intercept, or variance of the flltered process.0 200 400 600 800 1000 1200 ?1 0 1 2 3 4 Lags, t R x x(t) Data Equation ?2b e?1 3?? Bound Figure 2.4: Sample Autocorrelation: 1.7x105 time units, fs = 5 Hz, 2b = 4, ? = 200sec 22 Unfortunately, the experimental autocorrelation shares the same practical drawback with the Allan variance: the accuracy of the autocorrelation is dependent upon the length of the dataset. In general, the data set must be a length su?ciently longer than the time constant of interest. Experiments performed by the author in an attempt to extract reasonable parameters from raw gyro and accelerometer data proved very di?cult as long data sets often carried un-modeled disturbances (such as temperature efiects, long initial settling time). For a stochastic model with an autocorrelation in the form of Equation 2.6, the upper bound of the variance of the empirically-derived autocorrelation, Rbbexp at any lag-value can be expressed by Equation (2.15) [23]. VAR[Rbbexp] <= 2 4b? Td (2.15) Where, 2b is the variance, ? is the time-constant, and Td is the time-length of the experimental data set. The bounds plotted in Figure 2.4 show the large uncertainty in the Gauss-Markov process experimental autocorrelation for a data set of length much longer than its time constant. For the speciflc sample of data generated in the flgure, the time length of data set was more than 150 times the time constant. 2.4.3 Implementation Issues As discussed for both methods, accurate identiflcation of the process for the slowly- varying stochastic behavior of inertial sensors requires a data set many times longer than the time constant of the process. For many inertial sensors the bias drift is a slow process 23 (large time constant) yet still contributes signiflcantly to the output. Efiective identifl- cation for such a process requires a very long data set, sometimes outside the range of logging capabilities or a feasible window of time. Additionally, since the sensor parame- ters are unknown, the required length of the data set is therefore unknown and several iterations may be necessary. Additionally, the approximate Gaussian models presented may not provide a su?cient characterization of their observed behavior, especially for lower grade sensors. The identifler of these irregular sensors must then resort to highly conservative parameter estimates which simply give rough values comparable to that supplied by the sensor?s manufacturer. Discussion on the supplied sensor speciflcations and the need for accurate identiflcation follows in the next section. 2.5 Experimental Quantiflcation and Identiflcation 2.5.1 Manufacturer Sensor Speciflcations Sensor speciflcation sheets give a simple overview of a sensor?s operating charac- teristics including measurement range, input power consumption, data format, as well as expected accuracy. The accuracy of the sensor, as it is limited by the magnitude of the stochastic processes on the output, is quantifled in varying detail across the broad spectrum of sensors and their listed speciflcations. The speciflcations often indicate pa- rameters which coarsely bound the expected accuracy based on two assumed stochastic characteristics: noise and bias. The noise is quantifled by the random walk parameter as introduced in the Allan variance identiflcation section. The random walk parameter, rw is deflned by Equation (2.12). Its units are usually given in two forms shown in Equation 24 (2.16). rw = [output units]pHz rw = R[output units]dt pHr (2.16) The random walk parameter is easily extracted from experimental Allan variance data as shown in the previous section. A wider array of speciflcations are listed to quantify the bias drift. As a minimum, sensor manufacturers publish a conservative maximum value or maximum standard devi- ation within which the properly calibrated sensor is expected to output. Others give bias drift quantiflcation in terms of a value ascertained from the Allan variance chart. Ref- erence [24] indicates that the bias variation or bias instability parameter listed in some sheets is the lowest point on the Allan variance chart. In addition to the magnitude of the bias drift, some manufacturers of higher grade devices include some indication of the speed at which this bias drifts from the mean value by a correlation time. In any case, inertial sensor manufacturers provide very little information in support of full stochastic characterization of the outputs. 2.5.2 Experimental Approach Upon review of the available literature, manufacturers provide only conservative bounds on a sensor?s stochastic characteristics. This information gives only a rough starting point for the accurate characterization of this research. This section presents the basic methodology by which su?cient identiflcation can be performed. A general methodology can be roughly outlined by the following steps: 25 1. Examine available speciflcations by manufacture 2. Determine required sampling frequency and duration of data set based on specs 3. Acquire completely static data set (on level surface) 4. Remove constant bias (subtract the mean) 5. Run the Allan variance of the data set, extract the wide-band noise magnitude 6. Filter(zero-phase fllter) the data to reveal the underlying moving bias 7. Process the flltered data in the autocorrelation and extract time constant The success of experimental identiflcation is often di?cult in practice due to many of the reasons discussed in the preceding sections. The main di?culties are that the sensor drift model is only an approximation, the approximate model parameters are unknown, and the bias can not be fully isolated. As a result, the general methodology remains a long and highly iterative process. The general process of stochastic sensor error parameter identiflcation is demonstrated in Appendix A following the chapters of this thesis. The appendix presents the use of the techniques of this chapter to identify the stochastic model parameters with experimental data from an automotive-grade IMU. 2.6 Conclusion In this chapter, a simple inertial sensor error model that provides an approximation to the stochastic behavior observed on a a wide-range of inertial sensors grades and types has been presented. The model consists of the sum of two stationary, Gaussian random processes which describe the short-term and long-term stochastic behavior of static iner- tial sensor outputs. The techniques of Allan variance and autocorrelation were presented 26 as means to identify the three parameters required to describe the assumed model form from experimental sensor data. The experimental procedure was then outlined and dif- flculties of the method presented. In the following chapters, each of the error processes as detailed in this chapter are used for the purpose of quantifying the error in position, velocity, and attitude when dead-reckoning with an IMU in various kinematic scenarios. 27 Chapter 3 Covariance Propagation of Stochastic Errors 3.1 Introduction Bounded accuracy in inertial navigation depends upon regular position, velocity, and attitude measurements to compensate for the error growth of the integrated IMU signals. Many navigation methods employ regular measurements from GPS sensors, vision, odometry, and other sources of velocity, position, and attitude data. These measurements are often fused together in a navigation Kalman fllter resulting in an optimal estimate of the vehicle?s state. As GPS requires an unobstructed line-of-sight to at least four satellites, it fails to provide accurate data when traveling under bridges, heavily wooded areas, and in downtown city streets where tall buildings bound the path of the receiver. Under such conditions when the GPS data becomes unavailable, the Kalman fllter reduces to a simple algorithm in which the navigation states are derived solely from the integrated outputs of the inertial sensors initialized to the last \best estimate". As all inertial sensors are inherently corrupted with stochastic type errors (as introduced in Chapter 2), the integration of these signals cause the uncertainty in the resulting navigation states to increase with each step in time. As a result, the error in the estimated position, velocity, and attitude states grow with time. It is the goal of this chapter to quantify the the error growth due to the integrated stochastic errors present on the inertial measurements. By using the stochastic models from Chapter 2, this quantiflcation is achieved by deriving expressions for the variance of the integrated sensor errors using a technique 28 modeled after [11]. The variance expressions are then validated using a Monte Carlo simulation. For a single-axis sensor, the variance describes the error expected when dead-reckoning in one degree of freedom motion. The chapter concludes with a sensitivity analysis illustrating the in uence of the stochastic model parameters on the variance of the integrated sensor. 3.2 Simplifled Navigation Scenario: Single Axis The motion of a navigating body in the inertial or navigation frame can be derived from body-flxed inertial measurements of acceleration, a, and rotation rate, g. To attain the vehicle states of velocity, orientation, and position in the navigation frame in the general sense, the body-frame measurements are transformed using nonlinear difieren- tial equation relationships (shown later in Chapter 5). As a building block in providing an analysis of the general navigation scenario, preliminary attention is flrst turned to a simple one degree-of-freedom (1-DOF) motion scenario in which a single axis gyro or single axis accelerometer is integrated to provide the navigation frame states in its com- ponent direction. Equations (3.1-3.3) show the single-axis vehicle states of orientation, velocity, and position (?, V, P) as derived from the integrations of the corresponding inertial measurements. ? = Z gdt (3.1) V = Z adt (3.2) P = Z V dt (3.3) 29 It is the task of the following sections to quantify the error resulting from the use of Equations (3.1-3.3) with inertial sensors modeled as shown in Chapter 2. 3.3 General Characterization of Raw Sensor Measurement Recall the simple sensor model from Chapter 2 as described by Equation (2.1) ymeas = (SF)y +?+b As stochastic terms are assumed zero-mean, the mean value of the sensor output is simply the deterministic terms. E[ymeas(t)] = E[(SF)y(t)]+E[?(t)]+E[b] = (SF)y(t)+b (3.4) The variance of the sensor output is equal to the variance of the stochastic error. VAR[ymeas] = E[y2meas(t)]?E[ymeas(t)]2 = E[((SF)y +?(t)+b)((SF)y(t)+?(t)+b)]?((SF)y(t)+?(t)+b)2 ymeas2(t) = ?2(t) (3.5) Since the error sources are Gaussian and uncorrelated, the variance in the sensor output can be expressed as the sum of the variances of the two contributing error sources: wide band noise, !, and Gauss-Markov process, b, as introduced in Chapter 2. ?2(t) = !2(t)+ b2(t) (3.6) 30 3.4 Characterization of Integrated Sensor Measurement Given a sensor modeled with Equation (2.1) and whose stochastic errors are charac- terized as the sum of the independent noise sources introduced in Chapter 2, the purely integrated sensor is characterized as follows. The mean of the integrated sensor is simply the value of the integrated deterministic terms. E[ Z ymeas(t)dt] = E ?Z [(SF)y(t)+b]dt ? +E ?Z ?(t)dt ? = Z [(SF)y(t)+b]dt+0 = SF Z y(t)dt+bt (3.7) The variance of the integrated sensor output is the variance of the integrated independent stochastic error sources. E[ Z y2meas(t)dt] = E[ Z (SF)y +?(t)+bdt ? Z (SF)y(t)+?(t)+bdt ? ] ? Z (SF)y(t)+?(t)+bdt ?2 = Z E[?2(t)]dt (3.8) The variance of the integrated sensor output is equal to the sum of the variances of the independent integrated error sources. R ymeas2(k) = R !2(k)+ R b2(k) (3.9) 31 Applying the same methods to the double integrated case yields. E[ ZZ ymeas(t)dt] = E ?ZZ [(SF)y(t)+b]dt2 ? +E ?ZZ ?(t)dt2 ? = ZZ [(SF)y(t)+b]dt2 +0 = SF ZZ y(t)dt2 +bt2 (3.10) Performing analogous operations yields the resulting variance for the double-integrated sensor output. RR ymeas2(k) = RR !2(k)+ RR b2(k) (3.11) As shown in Equations (3.9) and (3.11), the variance in the integrated sensor output is equal to the sum of the variances of the integrated stochastic error sources. In the following sections the expressions of the individual variance functions of the random error processes and their integrals are derived. 3.5 Single Axis Stochastic Error Contributions Assuming that the numerical integration of the sensor is performed using an Euler approximation, the resulting integral values are simply scaled sums of the inertial values. As the scaled sum is a linear operation and the stochastic processes are Gaussian, the integrated stochastic processes are also Gaussian with transformed mean and variance functions. The following sections in this chapter derive the variance functions of the integrated and double-integrated wide-band noise and Gauss-Markov processes for nu- merical integration using the Euler method. The straightforward time-domain technique 32 modeled after [11] is employed here for derivations with the error models presented in Chapter 2. 3.5.1 Variance of Integrated Wide Band noise A derivation of the variance of integrated wide-band noise using the technique as follows has been shown in [11]. It is repeated here as an instructive example of the methodology by which the proceeding expressions are derived. Let ?y represent a wide-band noise process with variance 2!. ?y = ! (3.12) Integrating the ?y yields its integral value. _y = Z !dt (3.13) The above integration can be approximated by Euler?s method (left hand sum) with the initial condition _y0 = 0. _yk = _yk?1 +?t!k?1 = _y0 +?t k?1X i=0 !i (3.14) 33 Square both sides to obtain _yk _yk = ? ?t k?1X i=0 !i !? ?t k?1X i=0 !i ! (3.15) _yk _yk = ?t2 ?k?1X i=0 !i !?k?1X i=0 !i ! Take the expected value of the squared expression. E[_yk _yk] = E " ?t2 ?k?1X i=0 !i !?k?1X i=0 !i !# (3.16) The expectation of all of cross-terms of !i are equal to zero, as successive ! values in time are completely uncorrelated. The expression therefore reduces to E[_yk _yk] = ?t2E "k?1X i=0 !2i # = ?t2 k?1X i=0 E[!i!i] (3.17) The flnal result is an expression for the variance of integrated wide-band noise as a function of the variance of the wide-band noise, time index, and sampling interval. 2_y = 2!?t2k (3.18) 34 3.5.2 Variance of Double Integrated Wide Band noise Double integrating the wide-band noise process, ?y, to yield its double-integral value is shown below y = Z _ydt = ZZ !dt2 (3.19) Approximating the double-integration by Euler?s method (left hand sum), the following substitution is made yk = yk?1 +?t_yk?1 = y0 +?t2 k?1X j=0 ?j?1X i=0 !i ! (3.20) Simpliflcation of the double summation yields a single summation with an indexed coef- flcient. yk = ?t2 k?1X j=0 (k?j ?1)!j (3.21) Squaring both sides gives ykyk = ?t4 0 @ k?1X j=0 (k?j ?1)!j 1 A 0 @ k?1X j=0 (k?j ?1)!j 1 A (3.22) Taking the expected value of both sides with knowledge that successive values of !j are uncorrelated results in E[ykyk] = ?t4 k?1X j=0 (k?j ?1)2E[!j!j] (3.23) 35 Expansion of the summation leads to E[ykyk] = ?t4 0 @(k?1)2 k?1X j=0 (1)?2(k?1) k?1X j=0 j + k?1X j=0 j2 1 AE[!j!j] (3.24) Using the analytic solutions for power series summations the expression reduces to E[ykyk] = ?t4 k(k?1)2 ?2(k?1)12k(k +1)+ 16k(k +1)(2k +1) ? E[!j!j] (3.25) Further simpliflcation yields an expression for the variance of double integrated wide- band noise as a function of its variance, time index, and sampling interval. y2 = ?t4 !2 1 3k 3 + 1 2k 2 + 1 6k ? (3.26) 3.5.3 Variance of 1st order Gauss-Markov process The difierential equation for the 1st order Gauss-Markov process as given in Equation (2.4) can be realized using an Euler approximation bk = bk?1 +?t_bk?1 = bk?1 +?t?bk?1? +?t!bk?1 = 1? ?t? ? bk?1 +?t!bk?1 (3.27) For clarity in derivation, deflne A = 1? ?t? ? to get bk = Abk?1 +?t!k?1 (3.28) 36 The expression can be written as the following summation where the initial condition of the process is assumed to be zero, b0 = 0. bk = Ak?1b0 +?t k?1X i=0 Ak?i?1!bi bk = ?t k?1X i=0 Ak?i?1!bi (3.29) While this zero initial condition assumption simplifles the analysis below, the result- ing process may take time to settle into steady state depending on the size of the time constant. As this stochastic process is only an approximation to the observed sensor phenomenon, the signiflcance of the initial condition is uncertain. Squaring both sides obtains bkbk = ? ?t k?1X i=0 Ak?i?1!bi !? ?t k?1X i=0 Ak?i?1!bi ! (3.30) Applying the expectation operator to both sides with the knowledge that successive !bi values in time are completely uncorrelated and exhibit an identical variance results in E[bkbk] = ?t2 k?1X i=0 A2(k?i?1)E[!bi!bi] = ?t2A2k?2 k?1X i=0 A?2iE[!bi!bi] = ?t2A2k?2 k?1X i=0 A?2i 2!b (3.31) 37 Using the solution to the geometric series yields the following analytical expression E[bkbk] = ?t2A2k?2 1?A?2k 1?A?2 ? 2!b (3.32) Further simpliflcation results in the following expression for the variance of a 1st order Gauss-Markov process as a function of the variance of the driving noise, 2!b, time index, k, and sampling interval, ?t. 2b = ?t2 2!b A2k ?1 A2 ?1 ? (3.33) Note that for positive values of ?, A is less than one and therefore 2b will reach a steady-state value. 3.5.4 Variance of Integrated 1st order Gauss-Markov process Let ?x represent the bias drift as modeled by the Gauss-Markov process. Assume the process is realized by an Euler integration with zero initial condition as in Equation (3.27). ?xk = bk = Abk?1 +?t!bk?1 (3.34) 38 Next, approximate the integral using Euler?s method with initial condition _x0 = 0 _x = Z ?xdt = Z bdt = _xk?1 +?t?xk?1 = _x0 +?t?A?xk?2 +?t!bk?2? = ?t2 k?1X j=0 ?j?1X i=0 Aj?i?1!bi ! (3.35) The summation can be rewritten as _xk = ?t2 k?2X j=0 ? jX i=0 AjA?i!bi ! (3.36) Expand the summation for k = 5 and collect the !bi terms _x5 = ?t2 4X j=0 ?j?1X i=0 Aj?i?1!bi ! (3.37) = ?t2?!b0 +A!b0 +!b1 +A2!b0 +A!b1 +!b2 +A3!b0 +A2!b1 +A!b2 +!b3? = ?t2?!b0?A0 +A1 +A2 +A3?+!b1?A0 +A1 +A2?+!b2?A0 +A1?+!b3?A0?? Investigation of the expansion and rearrangement of Equation (3.36) results in simplifl- cation to the double summation _xk = ?t2 k?2X i=0 !bi 0 @ k?1?iX j=0 Aj 1 A (3.38) 39 Using the solution of the geometric series, the expression is further reduced to a single summation _xk = ?t2 k?2X i=0 1?Ak?1?i 1?A ? !bi = ?t 2 1?A k?2X i=0 ? 1?Ak?1?i ? !bi (3.39) Then both sides are squared to obtain _xk = ?t 4 (1?A)2 ?k?2X i=0 ? 1?Ak?1?i ? !bi !?k?2X i=0 ? 1?Ak?1?i ? !bi ! (3.40) Taking the expected value of both sides with knowledge that successive !bi values in time are completely uncorrelated results in E[_xk _xk] = ?t 4 (1?A)2 k?2X i=0 ? 1?Ak?1?i ?2 E[!bi!bi] = ?t 4 (1?A)2 k?2X i=0 1? 2A k AAi + A2k A2A2i ? 2!b = ?t 4 (1?A)2 ?k?2X i=0 (1)? 2A k A k?2X i=0 A?i + A 2k A2 k?2X i=0 A?2i ! 2!b (3.41) Using the solutions to the geometric series to simplify the summations gives E[_xk _xk] = ?t 4 (1?A)2 (k?1)? 2A k A 1?A1?k 1?A?1 + A2k A2 1?A2?2k 1?A?2 ? !2 (3.42) 40 This results in an expression for the variance of single integrated 1st order Gauss-Markov process in terms of the variance of the driving noise, 2!b, time index, k, and sampling interval, ?t. _x2 = ?t4 2!b 1+2A?2Ak ?2A1+k +A2k ?k +kA2 ?1+2A?2A3 +A4 ? (3.43) The above equation can be expressed in a condensed form _x2 = ?t4 2!b ? ?a1 +a2Ak ?a3A2k +a4k ? (3.44) Where the constants of Equation (3.44) are a1 = 1+2A?1+2A?2A3 +A4 a2 = ?2?2A?1+2A?2A3 +A4 a3 = 1?1+2A?2A3 +A4 a4 = A 2 ?1 ?1+2A?2A3 +A4 (3.45) 3.5.5 Variance of Double Integrated 1st order Gauss-Markov process Let x represent the double integration of the bias drift, b x = Z _xdt = ZZ bdt2 As before, Euler integration is used to realize the Markov process, _x, using the initial condition x0 = 0. The double integration the process is represented below as a series of 41 nested summations xk = ?t3 k?1X m=0 m?1X j=0 j?1X i=0 Aj?i?1!bi (3.46) Expanding the summations for k = 5 and collect the ! terms results in x5 = ?t3 4X m=0 m?1X j=0 j?1X i=0 Aj?i?1!bi = ?t3?A2!b0 +(!b1 +2!b0)A+!b2 +2!b1 +3!b0? = ?t3?!b0(A2 +2A1 +3A0)+!b1(A1 +2A0)+!b2(A0)? (3.47) Investigation of the preceding expansion leads to the simpliflcation of the triple summa- tion expression into two nested summations in which !i is removed from the innermost summation xk = ?t3 k?3X i=0 !bi k?3?iX j=0 (j +1)Ak?3?i?j = ?t3 k?3X i=0 !biAk?3?i k?3?iX j=0 (jA?j +A?j) (3.48) By substituting the following derivatives into Equation (3.48), the following expression is obtained d dA ?A?j? = ?1 A ?jA?j? jA?j = ?A ddA ?A?j? 42 where, the iA?i term can be decoupled as shown below xk = ?t3 k?3X i=0 !biAk?3?i 0 @ k?3?iX j=0 A?j ?A ddA k?3?iX j=0 A?j 1 A (3.49) Using the solution to the geometric series of A?i, the expression becomes xk = ?t3 k?3X i=0 !biAk?3?i 1?A?k+2+i 1?A?1 ? ?A ddA 1?A?k+2+i 1?A?1 ?? (3.50) Evaluation of the derivative allows further simpliflcation, yielding xk = ?t3 k?3X i=0 !biAk?3?i ? 1?A?k+2+i 1?A?1 ? ?A (?k +2+i)A?k+2+i A(1?A?1) ? 1?A?k+2+i A2(1?A?1)2 ?? = ?t3 k?3X i=0 !bi (A k?1?i +A?2?kA+k +iA?i) (A?1)2 (3.51) With the knowledge that successive!bi values in time are completely uncorrelated, squar- ing and applying the expectation operator gives E[xkxk] = ?t6 k?3X i=0 E[!bi!bi](A k?1?i +A?2?kA+k +iA?i)2 (A?1)4 = ?t3E[!bi!bi](A?1)4 k?3X i=0 A2k?2A?i +(A?2?kA+k)2 +(A?1)i2)+ 2Ak?1(A?2?kA+k)A?i +2Ak?1(A?1)Aii+ 2(A?2?kA+k)(A?1)i ? (3.52) 43 A distribution of the above summation yields E[xkxk] = ?t6E[!bi!bi](A?1)4 ? A2k?2 k?3X i=0 A?i +(A?2?kA+k)2 +(A?1)2 k?3X i=0 i2 + 2Ak?1(A?2?kA+k) k?3X i=0 A?i +2Ak?1(A?1) k?3X i=0 iAi + 2(A?2?kA+k)(A?1) k?3X i=0 i ! (3.53) Using the solutions to the geometric and power series allows the summations to be simplifled further E[xkxk] = ?t3E[!bi!bi](A?1)4 ? A2k?21?A ?2k+4 1?A?2 +(A?2?kA+k) 2(k?2)+ (A?1)216(k?3)(k?2)(2k?5)+ 2Ak?1(A?2?kA+k)1?A ?k+2 1?A?1 + 2Ak?1(A?1)(?A) ddA1?A ?k+2 1?A?1 + (A?2?kA+k)(A?1)(k?3)(k?2) ! (3.54) The flnal result is an expression for the variance of a double integrated Gauss-Markov process in terms of the in terms of the variance of the driving noise, 2!b, time index, k, 44 and sampling interval, ?t. x2 = ?t6 !b2 ? (?2?4A+2A4 +4A3 +2)k3 + (9?12A?6A2 +12A3 ?3A4)k2 + (?13+8A?8A3 +A4 +12Ak ?12A2+k)k + (6?12A2 ?12Ak +12A2+k +6A2k) ? (3.55) The above equation can also be expressed as x2 = ?t6 !b2 ? c1k3 +c2k2 + (3.56) (c3 +12Ak ?12A2+k)k + (c4 ?12Ak +12A2+k +6A2k) ? where, c1 = ?2?4A+2A4 +4A3 +2 c2 = 9?12A?6A2 +12A3 ?3A4 c3 = ?13+8A?8A3 +A4 c4 = 6?12A2 3.5.6 Summary of Results The preceding sections derived the variance expressions for the raw, integrated, and double integrated stochastic error processes of wide-band noise and exponentially- correlated noise. Tables 3.1 and 3.2 show the resulting expressions in summary. The 45 left-most column represents the level of integration and the right column indicates the corresponding variance functions. Recall that k is the time index and ?t is the sample interval. Table 3.1: Variance Contributions of Wide-Band Noise Integrals State Variance, 2(k) ! 2!R ! dt 2!?t2kRR ! dt2 ?t4 !2?13k3 + 12k2 + 16k? Table 3.2: Variance Contributions of 1st-Order Gauss-Markov Process Integrals State Variance, 2(k) b ?t2 !b2 ? A2k?1 A2?1 ? R bdt ?t4 2 !b ??a 1 +a2Ak ?a3A2k +a4k ? RR bdt2 ?t6 !b2 ? c1k3 +c2k2 +(c3 +12Ak ?12A2+k)k +(c4 ?12Ak +12A2+k +6A2k) ? 3.6 Validation of the Error Propagation In order to validate the derived expressions for the variance functions of the inte- grated stochastic processes, a Monte Carlo simulation was employed. The basic idea of the simulation is to generate a large number of independent stochastic processes using flxed parameters for a given window of time and then integrate (and double integrate) each simulated process over the duration. The variance over all the simulated runs is 46 computed for each time step. The computed variance functions represent the variance versus time of the integral processes. The empirically deduced variances are then com- pared to the derived expressions in Tables 3.1 and 3.2. The following plots show the analytical variance functions compared against the Monte Carlo results for each of the stochastic processes. In each example flgure, the simulated variance matches well to the analytic expression thus validating the expressions derived in the preceding sections. 3.6.1 Propagation of Wide-Band Noise Process Figures 3.1, 3.2, and 3.3 show the validation of standard deviation functions for a wide-band noise process, its integral, and its double integral, respectively. The plots show that the equations as derived and listed in the Tables match the variance generated in the Monte Carlo simulation. 47 0 10 20 30 40 50 60?3 ?2 ?1 0 1 2 3 Time ?(t) Sample Run Equation MonteCarlo Figure 3.1: Standard Deviation of Wide-Band Noise Process: 10Hz, ! = 1, 2000 Monte Carlo iterations 0 10 20 30 40 50 60 ?2 ?1.5 ?1 ?0.5 0 0.5 1 1.5 2 Time ?(t) Sample Run Equation MonteCarlo Figure 3.2: Standard Deviation of Integrated Wide-Band Noise Process: 10Hz, ! = 1, 2000 Monte Carlo iterations 48 0 10 20 30 40 50 60 ?80 ?60 ?40 ?20 0 20 40 60 80 Time ?(t) Sample Run Equation MonteCarlo Figure 3.3: Standard Deviation of Double Integrated Wide-Band Noise Process: 10Hz, ! = 1, 2000 Monte Carlo iterations 3.6.2 Propagation of Gauss-Markov Process Figures 3.4, 3.5, and 3.6 show the validation of standard deviation functions for a Gauss-Markov process, its integral, and its double integral. The plots show that the equations as derived and listed in the Tables match the variance achieved through the Monte Carlo simulation. Note that for the non-integrated Gauss-Markov process shown in Figure 3.4, the process variance reaches steady state as determined by Equation (3.33). 49 0 10 20 30 40 50 60?3 ?2 ?1 0 1 2 3 Time ?(t) Sample Run Equation MonteCarlo Figure 3.4: Standard Deviation of Gauss-Markov Process: 10Hz, b = 2, ? = 30, 2000 Monte Carlo iterations 0 10 20 30 40 50 60 ?40 ?30 ?20 ?10 0 10 20 30 40 Time ?(t) Sample Run Equation MonteCarlo Figure 3.5: Standard Deviation of Integrated Gauss-Markov Process: 10Hz, b = 2, ? = 200, 2000 Monte Carlo iterations 50 0 10 20 30 40 50 60 ?1000 ?500 0 500 1000 Time ?(t) Sample Run Equation MonteCarlo Figure 3.6: Standard Deviation of Double Integrated Gauss-Markov Process: 10Hz, b = 2, ? = 200, 2000 Monte Carlo iterations 3.7 Application Example While the results obtained in this chapter are in direct support of the more general navigation scenarios presented in later chapters, the expressions for the propagation of the errors in this chapter can be directly applied to a single-axis navigation scenario. It is the purpose of this section to illustrate the use of such expressions with such an example. Suppose a body is constrained to move in a straight line trajectory as depicted by Figure 3.7. Suppose additionally that the sensitive axis of an accelerometer is coincident with the traveling direction of the body. The accelerometer speciflcations of random walk, bias drift variance, time constant, and sample frequency are known and the sensor has been fully calibrated to remove any constant bias or efiects due to temperature. The 51 output of the calibrated measurement is integrated to obtain the velocity and position for a given acceleration proflle along a straight line position trajectory. For each time step, in addition to the values of velocity and position, the results of this chapter can provide the expected accuracy of the acceleration, velocity, and position. Figure 3.7: Body Constrained to Travel in One Direction Assume the sensor speciflcations for the accelerometer are given in Table 3.3. These speciflcations represent a low grade accelerometer with an exaggerated bias magnitude. Table 3.3: Sample Accelerometer Speciflcations Speciflcation Value fs 10 Hz 2! 0.5 ms2 2b 0.25 ms2 ? 200 s For any acceleration proflle, a, in a single direction, the mean value of the veloc- ity, V and position P in the same direction are simply the integral values of the true 52 acceleration. The mean values represent the Euler integration of true sensor outputs. E[V(k)] = ?t kX i=0 ai meters per second (3.57) E[P(k)] = ?t2 kX j=0 jX i=0 ai meters (3.58) The variance of the the velocity, V, is the sum of the variances of the error con- tributions from integrated wide-band noise and integrated Markov process as listed in Tables 3.1 and 3.2. 2V (k) = 2R ?(k) (3.59) = 2R !(k)+ 2R b(k) (3.60) The variance of the the position, P, is the sum of the variances of the error contri- butions from double integrated wide-band noise and double integrated Markov process as listed in Tables 3.1 and 3.2. 2P(k) = 2RR ?(k) (3.61) = 2RR !(k)+ 2RR b(k) (3.62) To demonstrate these results, consider the sinusoidal acceleration proflle as shown in Figure 3.8. This noisy accelerometer when integrated gives the velocity is shown in Figure 3.9. Another step of integration gives the position shown in Figure 3.9. As is evident by the plots, the 1- bounds show that the error growth becomes larger in time and with each level of integration. The bounds shown in the plots can be thought of 53 as a time-dependent corridor in which the integral values of velocity and position are expected to reside. As all error processes are Gaussian, the 1- bounds plotted in Figures 3.9 and 3.10 specify the region where approximately 66.7 percent of all trajectories are expected to travel. 0 10 20 30 40 50 60 ?2 ?1.5 ?1 ?0.5 0 0.5 1 1.5 2 Time, [s] Acceleration, [m/s 2 ] Sensor True Figure 3.8: Simulated Acceleration Proflle: 10Hz, ! = 0.5, b = 0.25, ? = 200 54 0 10 20 30 40 50 60 0 5 10 15 20 25 30 Time [s] Velocity, [m/s] Sample Run 1?? Bounds MonteCarlo Figure 3.9: Velocity with Bounds from Accel Proflle0 10 20 30 40 50 600 200 400 600 800 1000 Time, [s] Position, [m] Sample Run 1?? Bounds MonteCarlo Figure 3.10: Position with Bounds from Accel Proflle 55 3.8 Illustration of Results Inertial sensors exhibit varying magnitudes of each of their component error pro- cesses as well as the characteristic time constant of the drifting component. The follow- ing sections illustrate the efiect of the component stochastic process parameters on the growth of the variance of the total integrated sensor values. 3.8.1 Relative Magnitudes The two stochastic error processes, when integrated, each uniquely contribute to variance function of the integrated sensor output. A sensor?s wide-band noise component will dominate the error for short integration intervals while the drifting bias dominates for longer durations. This relative efiect of each can be investigated by observing the efiect of adjusting the ratio of the process magnitudes, b w, on the variance function of the integrated output for a flxed time constant. Setting the wide-band noise standard deviation to 1, Figures 3.11, 3.12, and 3.13 and show the 1- bounds of the integrated sensor output for various ranges of the bias drift standard deviation. Since the rate of increase of the integrated bias drift variance is higher than that of the integrated wide-band noise for any relative ratio, even small relative bias magnitudes will cause the bias to eventually dominate the variance growth. The general efiect of the relative ratio is that an increase Gauss-Markov process magnitude, b, causes it to dominate sooner, resulting in a larger rate of error growth in the integrated sensor output. 56 0 20 40 60 80 100 1200 2 4 6 8 10 12 14 16 Time ?(t) 1?? Wide?Band 1?? Sensor Ratio = 0.1 Ratio = 1 Figure 3.11: Integrated Sensor 1- Bounds for Ratio from 0.1 to 10 20 40 60 80 100 1200 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Time ?(t) 1?? Wide?Band 1?? Sensor Ratio = 0.01 Ratio = 0.1 Figure 3.12: Integrated Sensor 1- Bounds for Ratio from 0.01 to 0.1 57 0 20 40 60 80 100 1200 0.5 1 1.5 2 2.5 3 Time ?(t) 1?? Wide?Band 1?? Sensor Ratio = 0.001 Ratio = 0.01 Figure 3.13: Integrated Sensor 1- Bounds for Ratio from 0.001 to 0.1 3.8.2 Efiect of Time Constant For a flxed relative magnitude of 0.1, Figure 3.14 illustrates the shape of the bounds for various values of the Markov model time constant, ?. For the range of time constants shown, larger values of ? cause a slower increase in the rate of error propagation, while lower values indicate a faster increase. For a flxed relative magnitude of 0.1, Figure 3.15 illustrates the efiect of a larger range of ? on the value of the integrated sensor error bounds at a particular time of 120 units. For very small values of ?, the initial conditions of the Markov model dominate over the input noise. As the time constant increases, the maximum error value peaks and then levels ofi in a nonlinear fashion. For most inertial sensors, the time constant is usually much longer than values corresponding to the peak. The result of Figure 3.14 best describes the efiect of time constant within its expected range. 58 0 20 40 60 80 100 1200 1 2 3 4 5 6 7 Time ?(t) 1?? Wide?Band 1?? Sensor tau = 200 tau = 1000 Figure 3.14: Integrated Sensor 1- Bounds for Fixed Ratio of 0.1 and ? Between 200 to 1000 0 200 400 600 800 1000 1200 1400 1600 1800 3 4 5 6 7 8 9 ? ?(120) Figure 3.15: Value of Integrated 1- Bounds at time index of 120 for Fixed Ratio of 0.1 and ? Between 200 to 1000 59 3.9 Conclusion In this chapter, analytical variance expressions have been derived which quantify the error growth of subsequent integrations of inertial sensors exhibiting the assumed stochastic model forms of Chapter 2. A Monte Carlo simulation of the stochastic pro- cesses was used to validate the analytical results and further simulations illustrated the use of the expressions in the quantiflcation of accuracy for the single-axis case. The flnal sections of this chapter showed the relative and total efiects of the three stochastic model parameters on the resulting variance functions of the integrated sensor. The results of this chapter are used in direct support of derived expressions and analysis for the planar navigation scenario studied in Chapter 4. 60 Chapter 4 Two Dimensional Error Propagation 4.1 Introduction In order to describe the most general motion of a navigating vehicle, six degrees of freedom are required. An inertial measurement unit (IMU) attached to the navigating body takes measurements of these six degrees of freedom which include three orthogo- nal accelerations (ax,ay,az) and three orthogonal rotation rates (gx,gy,gz). In order to navigate within a suitable frame of reference such as on the surface of the earth, these measurements must be transformed and integrated to values of orientation, velocity, and position in that frame. For vehicles traveling within short ranges on the earth, a suit- able frame of reference is a simple cartesian coordinate system in which North, East, Down (NED) axes are aligned according to the right-hand rule. The orientation of the vehicle in this navigation frame can be described by roll (`), pitch ( ), and yaw (?) angles as deflned about the North, East, and Down axes, respectively. Figure 4.1 is a three-dimensional diagram depicting the body frame and navigation frame coordinate systems. It the task of the inertial navigator to perform the necessary operations on the body frame measurements to achieve the desired navigation frame values. This process of operating on the inertial measurements is referred to as inertial mechanization. For the 6- DOF scenario the mechanization equations are non-linear and require several calculations involving multiple measurements (see Chapter 5). If, however, the vehicle operates under some kinematic constraints, the resulting governing equations may be greatly simplifled. 61 Figure 4.1: Simplifled Coordinate Frame This chapter focuses its attention on the planar navigation scenario in which a vehicle is constrained to move on a North-East plane. Figure (4.2) depicts the simplifled motion of the constrained body. The body can translate in the north and east directions and rotate only about the direction orthogonal to the plane. The rotation rate sensed by a gyro about the z-axis on the constrained body is the same as the rotation rate about the down axis in the navigation frame. The yaw angle, ?, as measured positive from north about the down axis in the navigation frame requires no transformation to the navigation frame and is directly related to the sensed yaw rate, gz, by simple integration. The vehicle in this planar navigation scenario, depending upon its capability and trajectory, can operate under additional kinematic constraints. For general motion, the vehicle experiences side-slip, in which the body experiences velocity in both the direction it is pointing, body frame x known as heading (longitudinal direction), and the body frame y direction (lateral direction). This is the case on many vehicles in which high dynamic maneuvers force the body to point in a direction incoincident with its path of motion. Figure 4.3 shows the dynamic equations for use in navigating a 62 Figure 4.2: Simplifled Coordinate Frame 2D side-slipping vehicle with inertial measurements. As is evident by the diagram, the gyro is flrst integrated and then used in transforming the acceleration measurements to the navigation frame. The transformed accelerations are then integrated once for velocity, and twice for position. Figure 4.3: Navigation Relationships for Side-Slip Vehicle For many four-wheeled vehicles performing more moderate maneuvers, the velocity of the body can be assumed to be strictly coincident with its heading (body frame x). In this scenario, the body moves only in the direction it is pointing and is considered to be operating under the no-slip condition. For the no-slip case described above, the dynamic 63 relationships between the inertial measurements and navigation frame vehicle states of velocity and position simplify even further. These simplifled relationships for the no-slip case are shown in Figure 4.4. Similar to the case with side-slip, the navigation frame values of velocity and position are ascertained from the subsequent integrations. How- ever, in this no-slip case, the velocity is solely derived from integrating the longitudinal accelerometer, ax. Consequently, the acceleration, requires no transformation with the integrated gyro and the lateral accelerometer, ay, is not necessary for navigation. As a result, inertial navigation for the no-slip case requires one less integration and one less measurement. This observation is shown later in this chapter to ofier some potential dead-reckoning improvement when such kinematic assumptions are valid. Figure 4.4: Navigation Relationships for Vehicle with No-Slip Using results from and a similar approach to Chapter 3, this chapter presents a derivation of the propagation of navigation-frame velocity errors for the planar no-slip case as shown in Figure 4.4. For bodies that operate under such kinematic constraints, the resulting variance expressions provide the accuracy expected when dead-reckoning with the applicable dynamic equations and choice of inertial sensors. The resulting velocity expressions are validated with the same type of Monte Carlo simulations as used 64 in Chapter 3. Following the no-slip case is a discussion of the additional requirements for the case where a vehicleis expected to slip. Variance expressions for the transformation of accelerations for the no-slip case are presented. This chapter concludes by demonstrating the additional error induced when employing the slip equations for a non-slip trajectory. 4.2 Velocity Error in Navigation Frame Under No-Slip Planar Motion This section derives expressions for the variance of the 2-D velocity as derived from the minimum amount of inertial measurements necessary. The basic method used below deflnes the simplest governing inertial navigation relationships, substitutes the inertial measurement errors, and applies the expectation operator to the squared expressions. Using small angle approximations, the expressions can be simplifled to show a linear propagation of sensor errors from the IMU to the vehicle states. For the navigating body in the planar scenario, the yaw angle of the vehicle is derived by direct integration of the rotation rate sensed about the axis aligned orthogonal to the plane, gz. ? = Z gz dt (4.1) The velocity, V, of the body is always tangential to its path under the no-slip condition, and therefore can be derived from the acceleration sensed along the its x axis (ax accelerometer) by the integral relationship V = Z ax dt (4.2) 65 Assuming Euler integration of the above quantities, the preceding integrals are re- duced to summations Vk = ?t k?1X i=0 ax +V0 (4.3) ?k = ?t k?1X i=0 gz +?0 (4.4) The rate gyro and accelerometer inertial measurements are corrupted by stochastic error sources, ?g, and ?a as introduced and quantifled in earlier chapters. The total velocity and yaw angle can therefore be expressed as Vk = ?t k?1X i=0 (a+?a)+V0 = ?t k?1X i=0 a+?t k?1X i=0 ?a +V0 (4.5) ?k = ?t k?1X i=0 (g +?g)+?0 = ?t k?1X i=0 g +?t k?1X i=0 ?g +?0 (4.6) Fortheclarityinthederivationsfollowingitishelpfultoredeflnetheaboveequations with simpler notation. Note that despite the removal of subscript, k, the values are still 66 functions of time. A = ?t k?1X i=0 a+V0 (4.7) B = ?t k?1X i=0 ?a (4.8) fi = ?t k?1X i=0 g +?0 (4.9) fl = ?t k?1X i=0 ?g (4.10) In summary, A and fi represent the velocity and yaw angle as derived from Euler integration of the mean sensor outputs and B and fl represent the error due to the integrated accelerometer and rate gyro, respectively. 4.2.1 Mean and Variance of East Velocity The velocity in the east direction under this scenario is computed by taking the resultant velocity from the accelerometer and transforming it into the east component direction using the resultant yaw angle from the rate gyro. VEASTk = V sin(?) = ? ?t k?1X i=0 a+V0 +?t k?1X i=0 ?a ! sin(?t k?1X i=0 g +?0 +?t k?1X i=0 ?g) = (A+B)sin(fi+fl) (4.11) Using a trigonometric identity, the sine factor is expanded as shown below VEASTk = (A+B)(sin(fi)cos(fl)+cos(fi)sin(fl)) (4.12) 67 Provided that the integrated gyro error, fl, is su?ciently small, small-angle approx- imations of Equations (4.13) and (4.14) can be made. cosfl ? 1 (4.13) sinfl ? fl (4.14) In general, this small-angle approximation is valid within the typical range of inter- est. As the integrated sensor error grows outside the region for which the assumption is valid (jflj < 10 degrees), the resulting velocity and position errors exceed the range of accuracy in which this research seeks to quantify. The small angle approximation results in the simplifled expression for the east ve- locity VEASTk = (A+B)(sin(fi)+cos(fi)fl) (4.15) The mean function of velocity in the east direction is found by taking the expected value. E[VEASTk] = E[(A+B)(sin(fi)+cos(fi)fl)] = Asinfi+sinfiE[B]+AcosfiE[fl]+cosfiE[Bfl] (4.16) Since the stochastic errors are zero mean, the mean of the east velocity is simply the true value ascertained from the true part of the inertial measurements. E[VEASTk] = Acosfi (4.17) 68 To derive the variance of the east velocity, it is flrst squared as shown below (VEASTk)2 = ?Asinfi+sinfiB +Acosfifl +cosfiBfl?2 = A2 sin2(fi)+A2 cos2 fifl2 +sin2 fiB2 +cosfi2B2fl2 +2A2 sinficosfifl +2ABsin2 fi+4ABsinficosfifl +2ABcos2 fifl2 +2B2 sinficosfifl (4.18) Next, the expected value of the squared expression is taken as follows E?(VEASTk)2? = A2 sin2(fi)+A2 cos2 fiE?fl2?+sin2 fiE?B2?+cosfi2E?B2fl2? +2A2 sinficosfiE?fl?+2Asin2 fiE?B?+4AsinficosfiE?Bfl? +2Acos2 fiE?Bfl2?+2sinficosfiE?B2fl? (4.19) Assuming the two stochastic sources are independent, the expression reduces to E?(VEASTk)2? = A2 sin2(fi)+A2 cos2 fiE?fl2? +sin2 fiE?B2?+cos2 fiE?B2?E?fl2? (4.20) Computing the variance of the expression gives the following result VAR?(VEASTk)2? = E?(VEASTk)2??E?(VEASTk)?2 = A2 cos2 fiE?fl2?+sin2 fiE?B2?+cos2 fiE?B2?E?fl2? (4.21) 69 Recall from Equations (4.8) and (4.10), that B is the integrated accelerometer error and B is the integrated gyro error. The variances of the integrated inertial sensors were derived in Chapter 3 and shown again here E?B2? = 2R ?a (4.22) E?fl2? = 2R ?g (4.23) Where, 2R ?a is the variance of the integrated accelerometer errors and 2R ?g is the variance of the integrated gyro errors. The expressions are substituted below to yield the variance function of east velocity 2VEAST(k) = V 2 cos2 ? 2R ?g +sin2 ? 2R ?a +cos2 ? 2R ?a 2R ?g (4.24) Equation (4.24) shows that the variance is a three-termed expression consisting of the trajectory (velocity and heading) and the variances of the integrate accelerometer and integrated gyro. The last term includes the product of the two latter variances indicating that the variance of east velocity is not Gaussian distributed. However, for su?ciently large values of velocity, V , the flrst term dominates and the east velocity is approximately Gaussian. 4.2.2 Mean and Variance of North Velocity The velocity in the North direction under this scenario is computed by taking the resultant velocity from the accelerometer and transforming it into the North component 70 velocity using the resultant yaw angle from the rate gyro. VNORTHk = V cos(?) = ? ?t k?1X i=0 a+V0 +?t k?1X i=0 ?a ! cos(?t k?1X i=0 g +?0 +?t k?1X i=0 ?g) = (A+B)cos(fi+fl) (4.25) Using a trigonometric identity, the cosine factor is expanded as shown below VNORTHk = (A+B)(cos(fi)cos(fl)?sin(fi)sin(fl)) (4.26) Using the same small angle approximations shown before in Equations (4.13) and (4.14), the expression reduces to VNORTHk = (A+B)(cos(fi)+sin(fi)fl) (4.27) The mean function of velocity in the north direction is found by taking the expected value E[VNORTHk] = E[(A+B)(cos(fi)+sin(fi)fl)] = Acosfi+cosfiE[B]+AsinfiE[fl]+sinfiE[Bfl] (4.28) Since the stochastic errors are zero mean, the mean of the north velocity is simply the true value ascertained from the true part of the inertial measurements E[VNORTHk] = Acosfi (4.29) 71 To derive the variance of the north velocity, it is flrst squared as shown below (VNORTHk)2 = ?Acosfi+cosfiB +Asinfifl +sinfiBfl?2 = A2 cos2(fi)+A2 sin2 fifl2 +cos2 fiB2 +sinfi2B2fl2 +2A2 cosfisinfifl +2ABcos2 fi+4ABcosfisinfifl +2ABsin2 fifl2 +2B2 cosfisinfifl (4.30) Next, taking the expected value of the squared expression gives E?(VNORTHk)2? = A2 cos2(fi)+A2 sin2 fiE?fl2?+cos2 fiE?B2?+sinfi2E?B2fl2? +2A2 cosfisinfiE?fl?+2Acos2 fiE?B?+4AcosfisinfiE?Bfl? +2Asin2 fiE?Bfl2?+2cosfisinfiE?B2fl? (4.31) Assuming the two stochastic sources are independent and zero mean, the expression reduces to E?(VNORTHk)2? = A2 cos2(fi)+A2 sin2 fiE?fl2? +cos2 fiE?B2?+sin2 fiE?B2?E?fl2? (4.32) Computing the variance of the expression gives the terms VAR?(VNORTHk)2? = E?(VNORTHk)2??E?(VNORTHk)?2 = A2 sin2 fiE?fl2?+cos2 fiE?B2?+sin2 fiE?B2?E?fl2?(4.33) 72 Where, fi is the integrated accelerometer error and B is the integrated gyro error. The variances of these integrated inertial sensors are derived in Chapter 3 and shown again here E?B2? = 2R ?a (4.34) E?fl2? = 2R ?g (4.35) Where, 2R ?a is the variance of the integrated accelerometer errors and 2R ?g is the variance of the integrated gyro errors. Substituting the known integrated error variances gives the variance function for the north velocity. 2VNORTH(k) = V 2 sin2 ? 2R ?g +cos2 ? 2R ?a +sin2 ? 2R ?a 2R ?g (4.36) As shown above for the east velocity in Equation (4.24) it is evident that for rela- tively large velocities, the north velocity is approximately Gaussian due to the dominance of the flrst term of Equation (4.36) over the third term. 4.2.3 Cross Covariance of North and East Velocity To supplement the characterization of the velocity error for this planar no-slip sce- nario, the cross-covariance between the north and east components is derived as follows. The cross covariance is deflned as COV?VNORTHkVEASTk? = E?(VNORTHk ?E[VNORTH])(VEASTk ?E[VEAST])? (4.37) 73 Using the linear approximation from Equations (4.13) and (4.14) results in COV?VNORTHkVEASTk? = E?(Bcosfi?Asinfifl ?Bsinfifl)(Bsinfi+Acosfifl +Bcosfifl)? (4.38) Next, the expectation operator is expanded and shown below COV?VNORTHkVEASTk? = E h B2 cosfisinfi+Acos2 fifl +cos2 fiB2fl ?Asin2 fiBfl ?A2 sinficosfifl2 ?2AsinficosfiBfl2 ?sin2 fiB2fl ? sinficosfiB2fl2 i = cosfisinfiE?B2?+Acos2 fiE?fl?+cos2 fiE?B2fl? ?Asin2 fiE?Bfl??A2 sinficosfiE?fl2? (4.39) ?2AsinficosfiE?Bfl2??sin2 fiE?B2fl?? sinficosfiE?B2fl2? (4.40) Since the integrated accelerometer error, B, and integrated gyro error, fl, are zero mean and uncorrelated, the expression reduces to COV?VNORTHkVEASTk? = cosfisinfiE?B2??A2 sinficosfiE?fl2?? sinficosfiE?B2fl2? (4.41) 74 The flnal result is an expression for the cross-covariance of the velocity in the north and east component directions as shown below COV?VNORTHkVEASTk? = sinficosfi ? 2R ?a ?A2 2R ?g ? 2R ?a 2R ?g ? (4.42) 4.2.4 Probabilistic Characterization of Velocity Errors As discussed above, if the velocity is relative large, the north and east velocities ascertained from the no-slip mechanization are approximately Gaussian. Under this condition, the mean, variance, and cross covariance functions derived above completely characterize the propagation of the stochastic sensor errors into the velocity errors for this planar no-slip scenario. The resulting probabilistic characterization is captured for each step in time by a two-dimensional Gaussian surface as described by Equation (4.43). fV (VN;VE) = 12? VN VE p1??2 exp h 1 2(1??2) ? VN 2VN + VE 2VE ? 2?VNVE VN VE !i (4.43) Where, ?, the correlation coe?cient is computed by ? = COV ?V NORTHkVEASTk ? VN VE (4.44) and 2VN and 2VE are the variances of velocity in the north and east component directions as derived in the previous sections. Note that the Gaussian surface of Equation (4.43) is centered about zero and thus represents only the deviation from the mean values, which are known from the analysis of the preceding sections in Equations (4.17) and (4.29). 75 The expressions derived for the variance of the velocity error in the planar no-slip scenario indicate that the shape of the Gaussian surface described by Equation (4.43) is completely characterized at any given time by the instantaneous velocity and yaw angle. The width of its spread is dependent upon both the parameters of the stochastic processes present on the measurements and the time since integration. The stochastic characterization of the two-dimensional no-slip velocity error is therefore, at any given time, independent of the trajectory history. In other words, the variance in the velocity error is a function only of (in terms of vehicle trajectory) the instantaneous velocity and yaw angle of the navigating body. This implies that the velocity error for a vehicle that is inertially navigating under the no-slip assumption with the same initial and flnal conditions can be fully characterized by the variance expressions of Chapter 3. 4.2.5 Validation of the Velocity Error Characterization The same Monte Carlo validation methodology as for the single axis case in Chapter 3 was used to validate the velocity variance expressions of the preceding sections. The basic idea of the simulation is to flrst deflne a North/East Position trajectory in time and numerically difierentiate the trajectory to obtain the North/East Velocities. From the velocity, derive a yaw angle assuming the body points tangential to the path and then difierentiate to obtain the body frame longitudinal acceleration and yaw rate (under a no-slip assumption). The stochastic errors are then added to the inertial values to simulate the body-frame sensor measurements. Next, the simulated measurements are transformed and integrated to attain the now-corrupted attitude and 2-D velocity. After collecting a su?cient number of simulated integrations, the mean, covariance, and cross- covariance values of the outputs are computed for each step in time. These resulting 76 variance functions are compared to the analytical expressions using the known simulation parameters. Since the mean of the resulting velocity components is simply the processed value of the true measurements, only the errors about the mean need be considered. To both illustrate and validate the expressions derived for the velocity accuracy, a candidate position trajectory is chosen as shown in Figure 4.5.0 10 20 30 40 50 60 ?10 0 10 20 30 40 50 60 70 East Position [m] North Position, [m] Start Finish Figure 4.5: Simulated Position Trajectory Sensor speciflcations for the accelerometer and gyro chosen relative to trajectory and are listed in Table 4.1. These speciflcations are taken from a rough identiflcation of a 6-DOF Crossbow IMU-400CD, a medium grade ($3K-$4K) inertial measurement unit. Appendix A demonstrates sensor parameter identiflcation on this particular device. 77 Table 4.1: Simulated Sensor Speciflcations Gyro Spec Value Accel Spec Value fs 10 Hz fs 10 Hz 2rw 0.14948 deg/s/sqrt(Hz) 2rw 0.000412 g/sqrt(Hz) 2b 0.0061183 deg/s 2b 0.000100 g ? 1300 seconds ? 500 seconds The deflned position trajectory is processed as described above, and the standard deviation of the velocity errors are computed. Figure 4.6 and 4.7 show comparisons of the computed velocity standard deviation to derived expression for the deflned trajectory for the east and north component velocities, respectively. The y-axis is the 3? value of the component velocity representing the corridor of near perfect certainty (greater than 99% probability) within where the velocities are expected to reside. As can be seen in the plots, the derived variance functions match well with the simulated data. The velocity variance expressions can be observed in Figures 4.6 and 4.7 to exhibit an oscillating type behavior. As can be seen in the derived Equations (4.24) and (4.36), the variance of the component velocities depends upon the instantaneous value of the Velocity, V, as well as the heading, ?. As ? is oscillating according to the sinusoidal position trajectory, the variance of the component velocities thus exhibits the efiect above. 78 0 10 20 30 40 50 60 ?0.3 ?0.2 ?0.1 0 0.1 0.2 0.3 Time, [s] 3? ? V EAST , [m/s] Sample Run Equation MonteCarlo Figure 4.6: Standard Deviation of East Velocity, 2000 Monte Carlo iterations0 10 20 30 40 50 60 ?0.1 ?0.05 0 0.05 0.1 Time, [s] 3? ? V NORTH , [m/s] Sample Run Equation MonteCarlo Figure 4.7: Standard Deviation of North Velocity, 2000 Monte Carlo iterations 79 4.3 Position Error in No-Slip Planar Motion Due to the added analytical complexity in computing the variance of the planar position for the no-slip scenario, this analysis has not been not completed in closed- form. Instead, this section will show by example of the derivation of east position, the limitations of the current analysis approach and lack of necessary information. The end results of the example derivation make clear the need for more complete statistical characterization of the integrated sensor errors, and is an avenue of future work as described in Chapter 6. Equation (4.45) below gives the velocity in the east direction using the small angle approximation from the earlier Equation (4.16). VEASTk = (A+B)(sin(fi)+cos(fi)fl) = Asinfi+sinfiB +Acosfifl +cosfiBfl (4.45) The east position can be derived using the Euler method of numerical integration. The result is expressed below in summation form with the initial condition PEAST0 = 0. PEASTk = ?t k?1X j=0 (Asinfi+sinfiB +Acosfifl +cosfiBfl) (4.46) The mean value of the east position is straightforward since the stochastic processes, B and fl are zero-mean and independent. The mean value is simply the numerical 80 integration of the deterministic terms in the summation. E[PEASTk] = ?t k?1X j=0 (Asinfi+sinfiE[B]+AcosfiE[fl]+cosfiE[Bfl]) = ?t k?1X j=0 Asinfi (4.47) Using the same approach from the previous sections, in order to derive the east position variance the expression is flrst squared. PEASTk2 = 0 @?t k?1X j=0 (Asinfi+sinfiB +Acosfifl +cosfiBfl) 1 A 2 = ?t2 ?k?1X j=0 Acosfi k?1X j=0 Acosfi+ k?1X j=0 sinfiB k?1X j=0 sinfiB + k?1X j=0 Acosfifl k?1X j=0 Acosfifl + k?1X j=0 cosfiBfl k?1X j=0 cosfiBfl + 2 k?1X j=0 Asinfi k?1X j=0 sinfiB +2 k?1X j=0 Asinfi k?1X j=0 Acosfifl + 2 k?1X j=0 Asinfi k?1X j=0 cosfiBfl +2 k?1X j=0 sinfiB k?1X j=0 Acosfifl + 2 k?1X j=0 sinfiB k?1X j=0 cosfiBfl + 2 k?1X j=0 AcosfiB k?1X j=0 cosfiBfl ! (4.48) 81 Since the processes are independent and zero mean taking the expected value of the squared expression cancels out all the \cross-terms" leaving the following expression E?PEASTk2? = ?t2 ? E 2 4 k?1X j=0 Acosfi k?1X j=0 Acosfi 3 5+E 2 4 k?1X j=0 sinfiB k?1X j=0 sinfiB 3 5+ E 2 4 k?1X j=0 Acosfifl k?1X j=0 Acosfifl 3 5+ E 2 4 k?1X j=0 cosfiBfl k?1X j=0 cosfiBfl 3 5 ! (4.49) Taking the variance of the above expression reduces to the following three terms VAR?PEASTk2? = ?t2 ? E 2 4 k?1X j=0 sinfiB k?1X j=0 sinfiB 3 5+ E 2 4 k?1X j=0 Acosfifl k?1X j=0 Acosfifl 3 5+ E 2 4 k?1X j=0 cosfiBfl k?1X j=0 cosfiBfl 3 5 ! (4.50) 82 Observing that each of the three terms in Equation (4.50) have the same basic form, a single expansion is illustrated with an expansion of the flrst term for k = 5. ?t2 ? E 2 4 k?1X j=0 sinfifl k?1X j=0 sinfifl 3 5 ! = ?t2 ? sin2 fi1B12 +sin2 fi2B22 +sin2 fi3B32 +sin2 fi4B42 + 2sinfi1 sinfi2B1B2 +2sinfi1 sinfi2B1B3 +2sinfi1 sinfi3B1B4 + 2sinfi1 sinfi2B2B3 +2sinfi1 sinfi2B2B4 + 2sinfi1 sinfi2B3B4 ? (4.51) Recognizing the pattern of expansion, this flrst term can be generalized to the following expression ?t2 ? E 2 4 k?1X j=0 sinfifl k?1X j=0 sinfifl 3 5 ! = ?t2E 2 4 k?1X l=0 sin2 filBl2 + k?1X j=0 2sinfij k?1X i=j+1 sinfiiBjBi 3 5 (4.52) Distributing the expectation operator yields the following expression for the expan- sion of the flrst term of Equation (4.50) ?t2 ? E 2 4 k?1X j=0 sinfifl k?1X j=0 sinfifl 3 5 ! = ?t2 ?k?1X l=0 sin2 filE?Bl2?+2 k?1X j=0 sinfij k?1X i=j+1 sinfiiE[BjBi] ! (4.53) 83 The variables within the argument of the expectation operators are the integrated accelerometer error, B. While Chapter 3 provides the integrated sensor error variance, E?Bl2?intheflrstsummation, noanalysishasbeendonetoascertaintheautocorrelation, E[BjBi] expression in the second double-summation term. While further simpliflcation may be possible by substituting the original integrated stochastic model form of B into the second term summation, the indexed coe?cients still remain. Provided that the expectations of Equation (4.53) can be ascertained, the simpliflcation of the expression still remains a di?cult task due to the fact that that summations of the coe?cients sinfi and sin2 fi may not be able to be simplifled or represented by a closed-form expression. As the goal of this thesis is to provide closed-form expressions which bound the propagation of error in terms of ascertainable vehicle states (whether computed or measured), the above analysis remains to be explored in future work. By example, the sample expansion of the flrst term of the variance expression of Equation (4.50) has shown di?culty in achieving a solid closed-form expression for the planar position for the no-slip case. In future work, other approaches of analysis may provide a more complete and successful characterization of the errors as propagated in the dynamic equations. However, with the use of Monte Carlo simulations, much insight can be gained regarding the behavior of inertial navigation in various scenarios. The remainder of this thesis will use simulation to exemplify claims resulting from the analysis of inertial navigation in more complex kinematic scenarios. The next section extends the current planar motion error analysis to a vehicle experiencing side-slip. 84 4.4 Propagation of Error in Planar Motion with Slip The above sections have derived variance functions that quantify the propagation of inertial sensor errors when the no-slip assumption is valid and corresponding dynamic relationships employed. Under the simplifled motion, the dynamic relationships between the body frame and navigation frame vehicle states are such that only one accelerometer and one gyro are required to describe the motion of the vehicle. However, a vehicle expe- riencing side-slip requires an additional measurement to describe all states of its motion. For the planar case with side-slip, two accelerometers mounted in the body frame x and y directions and a single gyro mounted orthogonal to the plane in the z direction can be used to derive the values of 2-D velocity and position (see Figure 4.2). Due to the kinematics of the slipping motion, the navigation frame states require a coordinate transformation of the two accelerometer measurements (see Figure 4.3). As this coor- dinate transformation is achieved with the use of the integrated gyro measurement, the resulting navigation frame acceleration components exhibit an unbounded error growth in time. As the states of velocity and position require additional integrations of the transformed accelerations, their accuracy will grow at even faster rates than that of the transformed accelerations. As a result, the slip-case navigation-frame acceleration, ve- locity, and position will exhibit much higher error growth rates than under the no-slip assumption. 4.4.1 Acceleration Error in Navigation Frame for Slip-Case In order to illustrate the error growth when the no-slip assumption can not be made, the following characterization of the transformed accelerations is shown. The 85 acceleration in the planar North and East directions when the vehicle experiences side- slip can be derived from body frame inertial measurements ax0 and ay0. The navigation frame accelerations are related to the measurements by the following relationships. aNORTH = ax0cos? ?ay0sin? (4.54) aEAST = ax0sin? +ay0cos? (4.55) Where, as in the no-slip case, the the yaw angle of the vehicle can be derived by Euler integration of the gz measurement. ?k = ?t k?1X i=0 gz +?0 (4.56) As the inertial measurements include a true acceleration, a, and a stochastic error component, ?a, the expressions are expanded as aNORTH = (ax +?ax)cos? ??ay +?ay?sin? (4.57) aEAST = (ax +?ax)sin? +?ay +?ay?cos? (4.58) Using the deflnition for the Euler-integrated yaw rate (Equation (4.9)) and Euler- integrated gyro error (Equation (4.10)) the expression becomes aNORTH = (ax +?ax)cos(fi+fl)??ay +?ay?sin(fi+fl) (4.59) aEAST = (ax +?ax)sin(fi+fl)+?ay +?ay?cos(fi+fl) (4.60) 86 As in the no-slip velocity derivation, the trigonometric factors can be expanded by using the same identities and small-angle approximation. These operations give aNORTH = a0x cos(fi)?a0x sin(fi)fl ?a0y sin(fi)?a0y cos(fi)fl (4.61) aEAST = a0x sin(fi)?a0x cos(fi)fl +a0y cos(fi)?a0y sin(fi)fl (4.62) In order to calculate the variance of the north velocity, the simplifled expression is flrst squared aNORTH2 = ? a0x cos(fi)?a0x sin(fi)fl ?a0y sin(fi)?a0y cos(fi)fl ?2 = a0x2 cos2 (fi)+a0x2 sin2 (fi)fl2 +a0y2 sin2 (fi)+a0y2 cos2 (fi)fl2 ? 2a0x2 cos(fi)sin(fi)fl +2a0xa0y cos(fi)sin(fi)+ 2a0xa0y cos2 (fi)fl +2a0xa0y sin2 (fi)fl + 2a0xa0y cos(fi)sin(fi)fl2 +2a0y2 cos(fi)sin(fi)fl (4.63) Anticipating the expectation operator, the cross terms are removed from the ex- pression as the error sources are independent and zero-mean. The accelerometer errors, ?a are substituted and the expectation is applied to the resulting expression. E?aNORTH2? = ax2 cos2 (fi)+ax2 sin2 (fi)E?fl2?+ay2 sin2 (fi)+ay2 cos2 (fi)E?fl2?+ E??ax2?cos2 (fi)+E??ax2?sin2 (fi)E?fl2?+ E??ay2?sin2 (fi)+E??ay2?cos2 (fi)E?fl2?+ (4.64) 87 Subtracting the mean squared value yields the variance of the acceleration in the North direction. E?aNORTH2? = ax2 sin2 (fi)E?fl2?+ay2 cos2 (fi)E?fl2?+ cos2 (fi)E??ax2?+sin2 (fi)E??ax2?E?fl2?+ sin2 (fi)E??ay2?+cos2 (fi)E??ay2?E?fl2?+ (4.65) The results from Chapter 3 can be substituted to yield the flnal variance of the North acceleration as resulting from the transformation of the acceleration measurements. VAR?aNORTH2? = ? ax2 sin2 (^?)+ay2 cos2 (^?) ? R ?gz2 + cos2 (^?) ?ax2 +sin2 (^?) ?ax2 R ?gz2 + sin2 (^?) ?ay2 +cos2 (^?) ?ay2 R ?gz2 (4.66) Using the same procedure as above, the variance of the east acceleration reduces to the following expression VAR?aEAST2? = ? ax2 cos2 (^?)+ay2 sin2 (^?) ? R ?gz2 + sin2 (^?) ?ax2 +cos2 (^?) ?ax2 R ?gz2 + cos2 (^?) ?ay2 +sin2 (^?) ?ay2 R ?gz2 (4.67) Where ^? is the integrated gyro measurement. In the resulting variance expressions for the slip case accelerations, ?ax2 and ?ay2 are simply the output variances of the accelerometers. While these values are bounded 88 in time, the integrated gyro variance R ?gz2 is not. This unbounded acceleration er- ror causes an even faster error growth in the integrated velocity and double-integrated position values. For the no-slip case, no such acceleration transformation is necessary. Consequently, the no-slip velocity and position error will grow at a slower rate than that for the slip case. This suggests that when such vehicle constraints are valid, better dead-reckoning performance is achieved when using the fewest possible measurements with the fewest possible integrations. The following sections seeks to support this claim with a simulation example. 4.5 Comparison of Slip and No-Slip Mechanizations Since for the no-slip case, only a single accelerometer is required and no transfor- mation is necessary, the velocity and position values derived with the simpler dynamic equations and therefore exhibit a slower error growth. The conclusion then is that if and only if the no-slip condition exists, better dead-reckoning performance can be achieved by using as few measurements as possible and few sensor integrations as possible. By using computer simulation tools, these claims are demonstrated as follows. A comparison of the two methods can be realized with the use of a Monte Carlo simulation as used for the validation of the variance expressions earlier in this chapter. Within this particular simulation, a position trajectory is deflned and the resultant vehi- cle orientation, velocities, and accelerations are computed under the no slip assumption. A Monte Carlo simulation is performed (3000 iterations) in which simulated inertial errors with assumed sensor error parameters are added to the inertial values of accelera- tions and rotation rates. The simulated inertial measurements are sent to two calculation routines: one in which the no-slip condition is assumed, the other in which the side-slip 89 assumption is assumed. For 3000 iterations, the variance of the velocity, position, and attitude for each time step is calculated and stored for comparison. Figure 4.8 shows a sample North vs. East position trajectory deflned in the navi- gation frame. Under the no-slip assumption, the 2-D velocity is computed and shown in Figure 4.9; the corresponding no-slip yaw angle is shown in Figure 4.10. Using the same sensor speciflcations as listed in Table 3.1 (with a faster sample rate of 100Hz), the inertial measurements are simulated and processed for both the no-slip case in Figure 4.4 and the slip case in Figure 4.3 for 2000 iterations. The variance functions of the yaw angle, 2-D velocity, and 2-D position are then computed for each case over the 3000 iterations. 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 8 9 North Position [m] East Position [m] Figure 4.8: Deflned Position Trajectory 90 0 10 20 30 40 50 600 0.2 0.4 0.6 0.8 1 1.2 1.4 Velocity [m/s] Time [s] Figure 4.9: Velocity Trajectory With No Side Slip0 10 20 30 40 50 60?50 0 50 100 150 200 Yaw Angle [deg] Time [s] Figure 4.10: Yaw Angle Trajectory With No Side Slip 91 Figure 4.11 shows the 3- bounds for the yaw angle under each kinematic assump- tion. Since the processing is identical, the attitude of the body under each assumption exhibits identical variance growth curves. However, Figures 4.12 (Velocity) and 4.13 (Po- sition) show a difierent result. Here, the growth in the integral states derived using the no-slip mechanization is slower than that using the slip mechanization thus exemplifying the efiect claimed. 0 10 20 30 40 50 60 0 5 10 15 20 Yaw Angle 3? ? [deg] Time [s] No Slip Mechanization Side Slip Mechanization Figure 4.11: 3- Bounds on Simulated Yaw Angle 92 0 10 20 30 40 50 600 0.5 1 1.5 Velocity RMS 3? ? [m/s] Time [s] No Slip Mechanization Side Slip Mechanization Figure 4.12: RMS 3- Bounds on Velocity0 10 20 30 40 50 600 10 20 30 40 50 Position RMS 3? ? [m] Time [s] No Slip Mechanization Side Slip Mechanization Figure 4.13: RMS 3- Bounds on Position 93 This example shows that when the no-slip assumption is valid, use of the no-slip dynamic equations results in more accurate dead-reckoning performance. In contrast, use of the side-slip equations result in a higher rate of error growth due to the unnecessary acceleration measurement and additional step of integration in the governing equations. It is emphasized again that the no-slip equations only provide a valid result when such a no-slip assumption can be made. 4.6 Conclusion This chapter has presented an analysis of the propagation of the stochastic inertial sensor errors into the position, velocity, and attitude vehicle states for a body restricted to planar motion. For this planar case, it is shown that the vehicle?s z axis is always aligned to the navigation frame ? axis and therefore its errors are simply the integrated gyro errors. These integrated sensor errors can be quantifled by direct application of the variance expressions derived in Chapter 3. Within this planar navigation scenario, two kinematiccasesarestudied: side-slipandnoside-slip. Forthenoside-slipcase, avehicle?s velocity vector is coincident with its heading, and velocity is ascertained from a single accelerometer integration. However, in the side-slip case the vehicle points away from its direction of travel and the velocity must be derived from an additional measurement and additional integration step. This chapter shows through analytical results and a simulation example that when a vehicle has no side-slip, better dead-reckoning accuracy is achieved by employing the simpler equations with fewer measurements for the no-slip mechanization, as compared to the added measurements of the side-slip mechanization. This chapter has also shown through the position derivation for the no-slip case, the di?culty in the current analysis approach. It is suggested that further analytical analysis 94 be performed using techniques which provide more comprehensive statistical information on the propagated inertial errors to analytically quantify the position errors in planar mechanization. 95 Chapter 5 Six DOF Analysis 5.1 Introduction As mentioned in Chapter 4, six degrees of freedom (6-DOF) are required to kine- matically describe the most general motion of a body in space. In modern-day inertial navigation these six degrees of freedom are typically measured with an inertial measure- ment unit (IMU) rigidly attached to the body. The IMU measures three orthogonal accelerations (ax,ay,az) and three orthogonal rotation rates (gx,gy,gz) for a total of 6 degrees of freedom. For the purpose of navigation, such measurements need to be trans- formed into a coordinate system suitable for navigating. Figure 5.1 shows a common navigation-frame cartestian coordinate system as commonly employed on many inertial navigation systems, especially ground vehicles [25]. In one such cartesian frame, the navigating vehicle?s state is described by its velocity, and position in terms of coordi- nates North, East, and Down (NED) and by its roll (`), pitch ( ), and yaw (?) angles about the NED axes, respectively. The diagram from Chapter 4 is shown again in here in Figure 5.1 to illustrate the body and navigation frame axes. This chapter presents the equations necessary to use body frame measurements from a 6-DOF IMU to describe a vehicle?s state in the navigation frame. This 6-DOF navi- gation scheme is applicable to all of the scenarios introduced in the previous chapters as it captures the most general vehicle motion. As Chapters 3 and 4 have illustrated, the accuracy of any inertial navigation system degrades with time, sensor integrity, trajec- tory, and dynamic relationships employed without aid from external measurements. In 96 Figure 5.1: Simplifled Coordinate Frame speciflc regards to the efiect of dynamic relationships, Chapter 4 showed a comparative example showing the relative decline in inertial navigation performance when unneces- sary measurements and integrations are used in deriving the navigation-frame states. This chapter will extend this point with a similar example in which inertial navigation of a planar trajectory is compared using the planar equations of Chapter 4 and the general 6-DOF method presented in this chapter. 5.2 Equations of Motion 5.2.1 Orientation Two common conventions are used to describe a body?s orientation in the NED navigation frame: Euler angles and Quaternions. The Euler angle representation, while intuitive and straightforward to implement, exhibits a matrix singularity for pitch angles at 90 degrees. As this particular orientation is rarely encountered on many vehicles which use the NED frame, Euler angles remain as popular choice for ground vehicles. The 97 Quaternion approach involves a set of four linear difierential equations which describe orientation in three dimensional space. The method has the advantage of being immune to any particular singularities and therefore is numerically stable for all orientations. However, as the quaternion values themselves lack strong physical meaning, they are commonly transformed from and to Euler angles using a nonlinear relationship. The advantage of the quaternion approach for ground vehicle applications therefore is mainly in numerical computation and the relative accuracy compared to the Eulerian angle method is negligible. A more complete comparison and discussion of the two methods can be found in [26]. For the simple study in this thesis, the Euler angle representation is presented as follows. The Eulerian angular velocities are described in terms of the body frame rotation rates by the following set of flrst order difierential equations 2 66 66 64 _` _ _? 3 77 77 75 = 1 cos 2 66 66 64 cos sin`sin cos`sin 0 cos`cos ?sin`cos 0 sin` cos` 3 77 77 75 2 66 66 64 gx gy gz 3 77 77 75 (5.1) where gx, gy, and gz are the rotation rates as aligned to orthogonal axes on the body and (`), ( ), and (?) are the rotation rates about the navigation frame axes. The nonlinear relationships of orientation given by Equations (5.1) are numerically integrated to obtain the resulting Euler angles describing the attitude of the navigating body in space. It is instructive to note that when two angles are zero, the angular rate correspond- ing to the remaining angle holds a one-to-one relationship from the body frame to the navigation frame. Additionally, if the body frame rotation rates gx and gy are zero then 98 the relationship between the navigation-frame yaw rate, ?, and body-frame rate gz ex- hibit a one-to-one relationship. The latter situation is precisely the planar case studied in Chapter 3 in which the orientation of the body was constrained to rotate only about its z-axis. 5.2.2 Translation The accelerations as measured in the body-frame, must be transformed into the navigation frame using the Euler angles obtained from the orientation calculations. The Euler angles are used to construct the direction-cosine rotation matrix [26], which simply re-orients the three accelerations as measured in the body frame, to the navigation frame North, East, Down directions. The relationships between body frame and navigation frame accelerations are shown by Equation (5.2). 2 66 66 64 aN aE aD 3 77 77 75 = 2 66 66 64 cos cos? ?cos`sin? +sin`sin cos? sin`sin? +cos`sin cos? cos sin? cos`cos? +sin`sin sin? sin`cos? +cos`sin sin? ?sin sin`cos cos`cos 3 77 77 75 2 66 66 64 ax ay az 3 77 77 75(5.2) Once the accelerations are transformed, the gravity component is subtracted from the down acceleration to yield the kinematic acceleration of the body in the navigation frame. The resulting velocity and positions as described by the North, East, Down coordinate system are then derived by direct integration of the transformed accelerations. This process in which the IMU outputs are transformed and integrated into usable navigational quantities is known as mechanization. The mechanized IMU rigidly at- tached to the navigating body is considered the Inertial Navigation System (INS). The 6-DOF IMU mechanization algorithm as introduced above is summarized by Figure 5.2. 99 The simple mechanization shown here neglects the coriolis, centripetal accelerations, and other efiects experienced by an IMU moving on the rotating earth. For the short range (short time) for which many ground vehicles travel, these efiects are small and the mechanization discussed here is su?cient to support the typical requirements. Figure 5.2: Mechanization of IMU Measurements See Appendix B for a demonstration of the mechanization equations as presented. This appendix presents the mechanization of a medium grade IMU 6-DOF and compares its performance to position and velocity from a high-accuracy difierential GPS receiver. 5.3 Comparison to Planar Mechanization In Chapter 4 it was shown that for a planar vehicle trajectory where the body experiences no side-slip, the position and velocity error growth using the assumed-slip planar mechanization exhibited a faster variance growth than results from the no-slip planar method. The faster error growth observed with the slip equations was due to two contributing factors. The slip case required both an acceleration measurement and an additional step of numerical integration in the computation of its velocity and position 100 values. The 6-DOF navigation equations presented above are a natural extension of the 2-D slip case to the 6-DOF system. When the planar assumption can be made, using the 6-DOF equations will add unnecessary error into the system. As a result, the error growth for the 6-DOF case will be much worse. In the following, an example is shown using the same trajectory of Chapter 4 to show the amount of additional error induced when the 6-DOF method is employed. Figures 5.3 5.4 and 5.5 show a sample position, velocity, and yaw angle trajectory. Like the identical trajectory of the chapter 4 slip/no-slip comparison, the velocity, and yaw angle are derived from the deflned position trajectory.0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 8 9 North Position [m] East Position [m] Figure 5.3: Deflned Position Trajectory 101 ?0.1 0 0.1 0.2 0.3 0.4 ?1 ?0.5 0 0.5 1 North Velocity [m/s] East Velocity [m/s] Figure 5.4: Velocity Trajectory for Planar Motion0 10 20 30 40 50 60?50 0 50 100 150 200 Yaw Angle [deg] Time Figure 5.5: Yaw Angle Trajectory for Planar Motion 102 Using the sensor speciflcations as shown in Table 5.1, a Monte Carlo simulation was performed by computing 6-DOF position, velocity, and attitude with simulated acceleration measurements, ax, ay, az and rotation rates gx, gy, gz. Another Monte Carlo simulation was then performed with the planar side-slip equations of Chapter 4 using simulated ax, ay and gz with the same sensor speciflcations. Both simulations in this example were performed with 1200 iterations. The resulting standard deviations of the two dimensional rms velocity, rms position and the attitude angle are then plotted for comparison. Table 5.1: Simulated Sensor Speciflcations (Comperable to Crossbow IMU-400C) Gyro Spec Value Accel Spec Value fs 100 Hz fs 100 Hz 2rw 0.14948 deg/s/sqrt(Hz) 2rw 0.000412 g/sqrt(Hz) 2b 0.61183 deg/s 2b 0.000100 g ? 1300 seconds ? 500 seconds Figure 5.6 compares the standard deviation results for the yaw angle computed with the planar equations and then with the 6-DOF equations. It is evident that the error in the yaw angle for this trajectory is approximately the same. As the 6-DOF equations \take out" the efiect of the roll and pitch angles on the yaw angle, the yaw orientation is the same. Figures 5.7 and 5.8 compare the standard deviation results for the 2-D rms velocity and 2-D rms position, respectively. It is clear from these two plots that the additional measurements in the 6-DOF scheme cause a much more severe error growth as contrasted to the 3-DOF case. 103 0 10 20 30 40 50 600 5 10 15 20 Yaw Angle 3? ? [deg] Time [s] Planar Mechanization 6?DOF Mechanization Figure 5.6: 3- Bounds on Yaw Angle Mechanization Comparison0 10 20 30 40 50 600 10 20 30 40 50 60 70 80 Velocity RMS 3? ? [m/s] Time [s] Planar Mechanization 6?DOF Mechanization Figure 5.7: RMS 3- Bounds on Velocity Mechanization Comparison 104 0 10 20 30 40 50 600 200 400 600 800 1000 1200 1400 Position RMS 3? ? [m] Time [s] Planar Mechanization 6?DOF Mechanization Figure 5.8: RMS 3- Bounds on Position Mechanization Comparison Figure 5.9 and 5.10 show the standard deviation results for pitch and roll angles, respectively. These plots show precisely why the position and velocity growth is so large: the integrated stochastic errors on the additional roll and pitch measurements cause the accelerations to be rotated to an incorrect orientation. As a result, the mis- oriented accelerations are integrated along their incorrect directions to give a much more inaccurate North/East velocity and position. 105 0 10 20 30 40 50 600 2 4 6 8 10 12 Pitch Angle 3? ? [deg] Time [s] Planar Mechanization 6?DOF Mechanization Figure 5.9: 3- Bounds on Pitch Angle DOF Comparison0 10 20 30 40 50 600 2 4 6 8 10 12 Roll Angle 3? ? [deg] Time [s] Planar Mechanization 6?DOF Mechanization Figure 5.10: 3- Bounds on Roll Angle DOF Comparison 106 In the example above, it is clear that the body constrained to travel on the at plane not only lacked the beneflt of the extra gyro measurements - they greatly increased the overall error in the desired states. While this sample simulation has demonstrated the efiect of unnecessary inertial measurements in dead-reckoning navigation, it should be re-iterated that if the body was indeed traveling in motion that required a 6-DOF characterization, no less than six inertial measurements can be employed to correctly compute the vehicle?s trajectory. In other words, better navigation performance will most likely not be achieved by employing kinematic relationships which are simpler than the motion of the navigating body. In conclusion, when a vehicle in motion is under kinematic constraints, few measurements and integrations as required to completely describe the vehicle will yield the smallest amount of error. 107 Chapter 6 Conclusions 6.1 Overall Contributions This thesis has presented an analysis of the expected accuracy of inertial dead- reckoning in a variety of navigation scenarios. First, stochastic models were proposed as approximations to the observed behavior of many common types and grades of ac- celerometers and rate gyros. The selected inertial sensor error model included the sum of two independent Gaussian processes characterized by three parameters. Parameter iden- tiflcation methods of Allan variance and experimental autocorrelation were presented as the means of extracting the model parameters from experimental data. Derivations of the variance of subsequent integrations of each sensor error component process were per- formed and validated using Monte Carlo simulations. An application of the integrated error source variance expressions were demonstrated for the single-axis navigation sce- nario in which a single accelerometer or single gyro is used to ascertain integral navigation states in its flxed direction. The single-axis results were then expanded to derive expres- sions for the propagation of the inertial sensor error into the mechanization equations used for planar navigation in which a traveling body experiences no side-slip. The re- sults of the no-slip position were discussed and shown to have limitations due to the analysis techniques used. The planar no-slip inertial navigation mechanization was then compared to the planar mechanization with slip for the purpose of showing the errors induced by additional integrations and measurements required of the slip mechanization. The flnal chapter concluded the quantiflcation of inertial navigation by presenting the 108 6-DOF kinematic relationships as applied in the navigation of a body unconstrained in three dimensional space. The concluding point of Chapter 4 was reiterated again with a simulation example illustrating the additional dead-reckoning error induced by employ- ing unnecessary additional integrations and measurements in the 6-DOF equations for a simple planar trajectory. In summary, this research has provided 1. Approximate stochastic models which capture the necessary behavior of accelerom- eter and rate-gyroscope outputs as measured from various grade inertial measure- ment units. 2. Derivationsofthevarianceofthenumericallyintegratedvaluesoftheinertialsensor error sources of wide-band noise and exponentially correlated noise (Gauss-Markov process) from the error models. 3. Derivations of the variance of the 2-D (North/East) velocity error for a planar nav- igation scenario in which a vehicle experiences no side-slip. The position derivation is performed to show the need for methods of analysis superior to that of this thesis. 4. Comparison of the no-slip and slip and general 6-DOF inertial navigation equations for the planar case. 5. The six degree of freedom equations used for the most general navigation scenario. In addition to the derivations, this thesis has demonstrated the claim that the fewest number of integrations and measurements required to provide valid vehicle states in an inertial navigation system yields the minimal amount of error. It was shown through a simulation comparison of the planar no-slip, planar slip, and 6-DOF general 109 motion methods, that the inertial navigation with the simplest system had the best dead- reckoning performance (when the real vehicle trajectory was described by the simplest set of equations). 6.2 Di?culties The quantity and nature of approximations presented in this thesis warrant a dis- cussion of the limitations of the contributing results. It should flrst be noted that the variance expressions for the integrated sensor errors in Chapter 3 are derived based on the assumption that inertial sensors can be modeled with the approximations in Chapter 2. The experimental identiflcation and validation of of simple models with real experi- mental data is a di?cult process due to the extremely large amount of data required and the inadequacy of the assumed model forms. Based on these di?culties, the accuracy of the variance expressions based on the assumed model form should only re ect the accu- racy of the experimental identiflcation techniques by employing conservative parameter estimates. The techniques used to ascertain the variance of the integrated inertial sensor errors in Chapter 3 proved useful in its straightforward and intuitive approach. However, as shown for the planar no-slip position variance derivation in Chapter 4, the autocorrela- tion information was needed to simplify the expression to closed-form. The straightfor- ward method of derivation in Chapter 3, while yielding the desired variance expressions, did not allow for simple identiflcation of higher order measures of statistics. Additional information such as the cross-correlation and autocorrelation are di?cult to attain with the techniques employed in this thesis. In conclusion, other techniques are suggested for 110 use in quantifying the propagation of stochastic errors through various transformations and integration. In Chapter 3, the variance of the integral values of the Gauss-Markov process as- sumes that the process is realized with a zero-initial condition at the onset of integration. This zero initial condition for some values of the Markov model parameters may not ac- curately re ect the desired nature of the process as it may exhibit some initial transients in reaching its stationary status. The developed expressions in Chapter 4 for the variance of the planar case errors were based on small angle linear approximations for trigonometric operations on the integrated gyro error. For errors outside the range for which the approximations are valid, the error will propagate non-linearly and the resulting value will no longer be a Gaussian variable. The variance expression will then be inadequate to fully characterize the error of the nonlinearly processed inertial errors. However, the size of integrated gyro errors necessary to break the linear approximation is well outside the range for which a system designer is interested. Therefore the approximation is rarely an issue for the use of the variance expressions as practical analysis tools. 6.3 Future Work In order to more fully and successfully characterize the stochastic outputs of ac- celeration and rotation rate outputs of many grade IMUs, it is suggested that much experimental data be taken and studied to ensure the feasibility of the approximations introduced in this thesis. As the stochastic models introduced in Chapter 2 do not in- clude terms for temperature or range of motion, a study of such efiects should be done to ensure that their in uence is negligible. Given that the stochastic models from Chapter 2 111 are su?cient to match the behavior of the sensor errors, their corresponding parameters should be identifled with higher confldence using appropriately the techniques of Allan variance and autocorrelation. Once the models are validated and parameters identifled, the resulting variance expressions in Chapters 3 and 4 should be tested with experimental data to verify their use in real inertial navigation scenarios. As discussed above, the methodology used to propagate the sensor errors has room for improvement. Other more comprehensive statistical analysis techniques may provide a better and more complete characterization of the inertial system errors and perfor- mance in dead-reckoning. Since the inertial navigation performance using the 6-DOF mechanization scheme depends highly on non-linear relationships, the Gaussian sensor errors do not propagate linearly into the velocity, position, and orientation states and thus exhibit non-Normal distributions. More complete probabilistic characterizations of the propagation of sensor errors into the applicable nonlinear equations may provide broader results to the navigation systems mechanized for general motion. Using validated stochastic models and improved methods of error propagation anal- ysis, future research may expand its scope to include analyses of dead-reckoning with inertial measurements in a larger range of scenarios. As other sensors such as odome- try, laser scanners, and vision, are continually being fused into the navigation system, knowledge of the dead-reckoning performance of these systems becomes increasingly de- sirable. In addition, more sophisticated GPS/INS algorithms such as Tightly-Coupled GPS/INS allow for dead-reckoning aid from any available satellites when the number visible are less than that required for a position flx. Firm quantiflcation of the increased performance from any type of vehicle constraint, GPS, or auxiliary sensors in the large 112 range of applicable mechanizations provides much valuable information to the navigation community at large. 113 Bibliography [1] R. E. Kalman, \A new approach to linear flltering and prediction problems," ASME: Journal of Basic Engineering, vol. 82, pp. 34{45, 1960. [2] M. Newlin, \Design and development of a gps intermediate frequency and imu data acquisition system for advanced integrated architectures," Master?s thesis, Auburn University, December 2006. [3] J. M. Clanton, D. M. Bevly, and A. S. Hodel, \Highway lane tracking using gps in conjunction with onboard imu and vision-based lane tracking measurements." Fort Worth, TX: Institute of Navigation: ION-GNSS Conference, September 2006. [4] M. M. Veth and J. Raquet, \Fusion of low-cost imaging and inertial sensors for navigation." Fort Worth, TX: Institute of Navigation: ION-GNSS Conference, September 2006. [5] C. R. Carlson, \Estimation with applications for automobile dead reckoning and control," Ph.D. dissertation, Stanford University, April 2004. [6] W. Travis, A. T. Simmons, and D. M. Bevly, \Corridor navigation with a lidar/ins kalman fllter solution." Las Vegas, Nevada: IEEE Intelligent Vehicle Conference, 2005. [7] S. Nassar, \Improving the inertial navigation system (ins) error model for ins and ins/dgps applications," Ph.D. dissertation, University of Calgary, November 2003. [8] S. Godha and M. E. Cannon, \Integration of dgps with a low cost mems-based inertial measurement unit (imu) for land vehicle navigation application." Long Beach, CA: Institute of Navigation GNSS, September 2005. [9] M. M. Tehrani, \Ring laser gyro data analysis with cluster sampling technique," vol. 412. SPIE: International Society for Optical Engineering, 1983, pp. 207{220. [10] H. Hou and N. El-Sheimy, \Inertial sensors errors modeling using allan variance." Portland, OR: Institute of Navigation GNSS, September 2003. [11] D. Gebre-Egziabher, \Design and performance analysis of a low-cost aided dead reckoning navigator," Ph.D. dissertation, Stanford, 2004. [12] H. Kim, J. G. Lee, and C. G. Park, \Performance improvement of gps/ins integrated system using allan variance analysis." Sydney, Australia: Internation Symposium on GNSS/GPS, December 2004. 114 [13] B. E. Grantham and M. A. Bailey, \A least-squares normalized error regression algorithm with application to the allan variance noise analysis method." Coronado, CA: Institute of Navigation: Position Location and Navigation Symposium, April 2006. [14] W.Flenniken, \Modelinginertialmeasurementunitsandanalyzingtheefiectoftheir errors in navigation applications," Master?s thesis, Auburn University, December 2005. [15] D. Bevly, \High speed, dead reckoning, and towed implement control for auto- matically steered farm tractors using gps," Ph.D. dissertation, Stanford University, August 2001. [16] A. Gelb, Applied Optimal Estimation, A. Gelb, Ed. M.I.T. Press, 1974. [17] S. J. Pethel, \Test and evaluation of high performance micro electro-mechanical system based inertial measurement units." Coronado, CA: Institute of Navigation: Position Location and Navigation Symposium, April 2006. [18] J. Rankin, \Gps and difierential gps: An error model for sensor simulation." IEEE, 1994. [19] K. J. Astrom and B. Wittenmark, Computer Controller Systems, 2nd ed. Prentice Hall, 1990. [20] K. Breitfelder, \Ieee standard speciflcation format guide and test procedure for single-axislasergyros," IEEEAerospaceandElectronicSystemsSociety, Tech.Rep., September 1995. [21] D. W. Allan, \Statistics of atomic frequency standards," vol. 54, no. 2. IEEE, February 1966. [22] L. C. Ng, \On the application of allan variance method for ring laser gyro per- formance characterization," Lawrence Livermore National Laboratory," Informal Report, October 1993. [23] R. G. Brown and P. Y. C. Hwang, Introduction to Random Signals and Applied Kalman Filtering, 3rd ed. Wiley, 1997. [24] W. Stockwell, \Bias stability measurement: Allan variance," Crossbow Technology, Inc., Tech. Rep. [25] T. N. Gillespie, Fundamentals of Vehicle Dynamics. Warrendale, PA: SAE, 1992. [26] P. H. Zipfel, Modeling and Simulation of Aerospace Vehicle Dynamics, J. Przemie- niecki, Ed. Reston, VA: AIAA, 2000. 115 [27] W. Travis, \Methods for minimizing navigation errors induced by ground vehicle dynamics," Master?s thesis, Auburn University, May 2006. 116 Appendices 117 Appendix A Stochastic Parameter Identification with an Automotive-Grade IMU This appendix serves to demonstrate the parameter identiflcation techniques of Al- lan variance and autocorrelation on automotive-grade inertial measurements. The tech- niques are used to identify parameters of the assumed stochastic model from accelerom- eters and rate-gyros from a Crossbow IMU-400C logged using a PC with Windows XP via an RS-232 serial connection. The accelerometers and gyros were logged at a sample frequency of fs = 5 Hz for approximately 48 hours in a climate controlled o?ce while resting on a level desk. The meanisremovedfromeachsensorlogandtheoutputsarefllteredtorevealanyunderlying drifting bias. Figure A.1 shows the three Crossbow accelerometers zero-phase flltered (Matlab flltfllt() function) with a second-order low-pass Butterworth fllter with cutofi frequency of 0.001 Hz. Figure A.2 shows the zero-phase flltered rate-gyros (second-order Butterworth with 0.001 Hz cutofi frequency). 118 0 10 20 30 40 50?2 ?1.5 ?1 ?0.5 0 0.5 1 1.5 2x 10 ?3 time [hr] Filtered Output [g] ax ay az Figure A.1: Filtered Accelerometer Outputs0 10 20 30 40 50?0.08 ?0.06 ?0.04 ?0.02 0 0.02 0.04 time [hr] Filtered Output [deg/s] gx gy gz Figure A.2: Filtered Gyro Outputs 119 As can be seen in both sets of flltered outputs, the sensors exhibit a very slow drifting behavior over the 48 hour logging period. While this slow drift may efiect navigation system performance for long spans of time, it is of primary interest to understand the drifting behavior within short time intervals (the intervals for which the sensors will be used to dead-reckon). In order to employ the identiflcation techniques using the assumed Gauss-Markov model, focus is turned to the most steady of the outputs: the z-axis gyro, gz, within its most constant interval from hour 20 to 45. Figure A.3 shows the raw and flltered data for the gz gyro within the interval selected. The data within this interval is flltered to remove the high frequency content of the wide-band noise, leaving the driflng bias for identiflcation with the autocorrelation method. The flltered data shown in the plot was processed with the same low-pass fllter as before with a cutofi frequency of 0.0005 Hz, a value determined after multiple iterations to yield the best identiflcation results with the data used.0 5 10 15 20 25?1.5 ?1 ?0.5 0 0.5 1 1.5 time [hr] Gyro Data [deg/s] RawFiltered Figure A.3: Raw and Filtered Data Within Selected Section 120 Figure A.3 shows the flltered data by itself. This plot indicates the relatively long term stability of the selected data, while revealing the characteristic bias drift of interest.0 5 10 15 20 25?0.02 ?0.015 ?0.01 ?0.005 0 0.005 0.01 0.015 Filtered Gyro Data [deg/s] Time [hr] Figure A.4: Filtered Data Within Selected Section 121 The sectioned and flltered data is then processed with an autocorrelation function. The result is shown in Figure A.5. The run reveals an exponential function such as that exhibited by the Gauss-Markov model. From this data, the time constant is extracted by selecting the intersection on the 1e horizontal line and the variance from the y-intercept. Note that these parameters are extracted only with the confldence re ected by the bound as shown which is referenced around the best-flt line.0 1000 2000 3000 4000 5000 ?1 0 1 2 3 4 5x 10 ?5 Time [s] Autocorrelation, R xx (t) Data Best Fit ?2b e?1 1?? Bound Figure A.5: Autocorrelation of Filtered Gyro Data Now that rough Gauss-Markov parameters are identifled, an Allan variance calcula- tion is performed on the raw experimental data within the section selected. Figure A.6 shows the results of the Allan variance. 122 100 101 102 103 104 105 10?4 10?3 10?2 10?1 Averaging Time, T [s] Root Allan Variance ? AV Exp Fit 3?? Bound Figure A.6: Allan Variance of Gyro Data Simple inspection of the Allan variance allows extraction of the random walk pa- rameter (y-intercept of the graph); this parameter quantifles the wide-band noise. Using the identifled parameters from the autocorrelation and the extracted random walk pa- rameter, a best flt line is constructed and plotted to represent the Allan variance of the experimental gyro data with the assumed model form. The 3- bounds, based on the averaging time and total length of the data set (25 hours @ 5Hz) are plotted along with the flt and experimental data. The flt, while showing a fair agreement with the experimental data, is only a conservative representation as re ected by the large and spreading bounds. 123 Using the techniques above with the gyro data set, the parameters used to represent the rate-gyro of the Crossbow IMU-400C are shown in Table A.1. Table A.1: Results of Crossbow Gyro Identiflcation Gyro Spec Value fs 5 Hz 2rw 0.14963 deg/s/sqrt(Hz) 2b 0.0061 deg/s ? 1100 seconds As shown on the plots of Allan variance and autocorrelation, the bounds on the accuracy of each are very large. These bounds need to be smaller to re ect high con- fldence in the parameter estimation. Higher levels of confldence require longer sets of data. However, as shown by the long-term behavior of the flltered IMU outputs, it is di?cult to ascertain a long set of static data stable enough to employ the techniques with the assumed model forms. 124 Appendix B Demonstration of 6-DOF Mechanization of Experimental Inertial Measurements Taken at Talledega Superspeedway This appendix serves to demonstrate the 6-DOF mechanization of an automotive grade 6-DOF IMU in a medium-sized vehicle traveling around a track at high speeds. The experimental setup is as follows. An Inflniti G-35 sedan was equipped with a single antenna Novatel StarflreTM DGPS ( < 10cm accuracy) receiver mounted on a roof rack roughly in the planar center of gravity (CG) of the car. An automotive-grade Crossbow IMU-400C 6-DOF IMU was attached rigidly to the console of the vehicle, which is roughly located at the vehicle?s CG. GPS position, velocity, and course was logged at 5 Hz from the GPS, and inertial measurements were logged at 133Hz from the IMU via RS-232 serial connections to a PC in the trunk running Windows XP. For a more complete description of the vehicle test-bed and data acquisition system see [27]. The data used in this appendix was logged with the experimental vehicle and sensors while driving 3 laps around the 2.66 mile track at Talledega superspeedway (Talledega, AL) for approximately 6.5 minutes. The track has bank angles of about 33 degrees in the turns, 16.5 degrees in the tri-oval, and 3 degrees in the straights. The inertial data logged during the laps was post-processed using the 6-DOF mechanization equations presented in Chapter 5 to provide position, velocity, and attitude in the North, East, Down (NED) coordinate system. The inertially-derived values are compared to the high-accuracy GPS information to show its performance in dead-reckoning. 125 Figure B.1 shows the North and East position from the 6-DOF mechanization of the IMU as compared to the GPS position. It is immediately evident that the dead- reckoning position from the IMU quickly drifts away from the true position of the track (represented by the sub-10cm GPS position). As time elapses for each lap around the track, the IMU-derived position, while coarsely resembling the ovular shape of the true trajectory, drifts with increasing variance from the actual position. ?1000 ?500 0 500 1000 ?500 0 500 1000 1500 Position from IMU North Position [m] East Position [m] StarFireCrossbow StartFinish Figure B.1: Position of Track from Starflre GPS 126 Figure B.2 shows the velocity of the runs around the track as reported by the GPS receiver. As is shown, the large track allowed for velocities higher than are generally allowed on public highways.0 1 2 3 4 5 6 7 0 20 40 60 80 100 120 GPS Velocity, [mph] Time [min] Figure B.2: Velocity from GPS 127 Figure B.3 shows a comparison of the yaw angles from the IMU and GPS as the vehicle travels around the track. For the high speed and moderate banked turns around the track, the true heading and course are expected to be the same (no vehicle side-slip). As time progresses however, the heading derived IMU drifts from the GPS-reported orientation due to the integration of the stochastic errors present on the three gyros used to compute the yaw angle.0 100 200 300 400 500 0 50 100 150 200 250 300 350 400 Angle [deg] Time [s] Course from GPS Heading from IMU Figure B.3: Yaw Angle (Heading) from IMU 128 Figure B.4 shows the roll and pitch angles as purely derived from the IMU. As only a single antenna GPS was logged during the test, no un-biased estimate of these angles was available and so only the IMU-derived angles are shown. With knowledge of bank angles of the track, however, it is evident that the values of roll match fairly well to the slope of the track in the turns and straights.0 50 100 150 200 250 300 350 400 ?5 0 5 10 15 20 25 30 35 40 Angle [deg] Time [s] roll angle , ? pitch angle, ? Figure B.4: Roll and Pitch Angles from IMU 129 Figure B.5 shows of comparison of the North and East component velocities from the IMU and GPS. As shown, the IMU component velocities resemble the shape of the trajectory but drift with increasing distance from the GPS values. 0 50 100 150 200 250 300 350 400?40 ?20 0 20 40 60 North Velocity [m/s] Time [s] Inertial GPS 0 50 100 150 200 250 300 350 400?50 0 50 East Velocity [m/s] Time [s] Inertial GPS Figure B.5: North and East Component Velocities 130 Figure B.6 shows of comparison of the North and East component positions from the IMU and GPS. As shown, the IMU component positions resemble the shape of the trajectory but drift with increasing distance from the GPS values. It is evident that when compared to the velocity subplots of Figure B.5, the position error grows at a faster rate. Due to the additional level of integration (once more than velocity) of the stochastic sensor errors, this accelerated rate in position error growth is expected. 131 0 50 100 150 200 250 300 350 400?1000 ?500 0 500 1000 1500 2000 North Position [m] Time [s] Inertial GPS 0 50 100 150 200 250 300 350 400?1500 ?1000 ?500 0 500 1000 1500 East Position [m] Time [s] Inertial GPS Figure B.6: North and East Component Positions 132