Customizable Stochastic High-Fidelity Model of the Sensors and Camera Onboard a Fixed Wing Autonomous Aircraft

The navigation systems of autonomous aircraft rely on the readings provided by a suite of onboard sensors to estimate the aircraft state. In the case of fixed wing vehicles, the sensor suite is usually composed by triads of accelerometers, gyroscopes, and magnetometers, a Global Navigation Satellite System (GNSS) receiver, and an air data system (Pitot tube, air vanes, thermometer, and barometer), and it is often complemented by one or more digital cameras. An accurate representation of the behavior and error sources of each of these sensors, together with the images generated by the cameras, is indispensable for the design, development, and testing of inertial, visual, or visual–inertial navigation algorithms. This article presents realistic and customizable models for each of these sensors; a ready-to-use C++ implementation is released as open-source code so non-experts in the field can easily generate realistic results. The pseudo-random models provide a time-stamped series of the errors generated by each sensor based on performance values and operating frequencies obtainable from the sensor’s data sheets. If in addition, the simulated true pose (position plus attitude) of the aircraft is provided, the camera model generates realistic images of the Earth’s surface that resemble those taken with a real camera from the same pose.


Introduction and Outline
The sensors onboard an autonomous aircraft measure various aspects of the aircraft real or actual trajectory x = x TRUTH and provide these measurements to the aircraft guidance, navigation and control (GNC) system.The outputs of these sensors, collectively known as the sensed trajectory x = x SENSED , represent the only link between the real but unknown actual trajectory and the GNC system in charge of achieving an actual trajectory that deviates as little as possible from the guidance targets.

SENSORS
x (t t ) = x TRUTH (t t ) x (t s ) = x SENSED (t s ) Note that although the actual aircraft state varies continuously in the real world, in simulation it is usually the outcome of a high frequency discrete integration process [2] that results in x (t t ) = x TRUTH (t t ), where t t = t • ∆t TRUTH .
The sensed trajectory x (t s ) = x SENSED (t s ), where t s = s • ∆t SENSED , is however intrinsically discrete, although the working frequency of the different sensors may vary.This article considers that all sensors operate at the same rate of ∆t SENSED , with the exception of the GNSS receiver and the onboard camera, which work at ∆t GNSS and ∆t IMG , respectively.It is also assumed that all sensors are fixed to the aircraft structure in a strapdown configuration and that their measurement processes are instantaneous and time synchronized with each other at their respective frequencies.Table 1 provides the default frequencies employed in this article as well as in the corresponding software implementation [1], where all values can be adjusted by the user.

Discrete
Table 1: Default frequencies of the states in the actual and sensed trajectories The sensed trajectory can be defined as a time stamped series of state vectors x = x SENSED that groups the measurements provided by the different onboard sensors, comprising the only view of the actual trajectory at the disposal of the GNC system: Note that when present, the super index represents the frame or reference system in which a certain variable is viewed; if two sub indexes are present, it implies that the vector goes from the first frame to the second.Table 2 makes use of the body frame F B , rigidly attached to the aircraft structure with origin in its center of mass [3], the NED frame F N also centered on the aircraft center of mass with axes in the North-East-Down directions [3], and the inertial frame F I , usually considered as centered in the Sun with axes fixed with respect to other stars [4].Additionally, the specific force f IB is defined as the non gravitational acceleration experienced by the aircraft body with respect to an inertial frame [5].The main objective of this article is to properly model the sources of error incurred by the different sensors onboard the aircraft, with special emphasis on the inertial sensors (accelerometers and gyroscopes) that provide the essential information on which the navigation system of low SWaP vehicles relies.The model parameters 1 OSP stands for Outside Static Pressure. 2OAT stands for Outside Static Temperature. 3TAS stands for True Air Speed. 4 AOA stands for Angle Of Attack. 5AOS stands for Angle of Sideslip. 6CAM stands for Camera.

Components
2 need to be customizable so the user can employ the values that better resemble the performances of the specific equipment being modeled.They also need to be stochastic to properly represent the nature of the different random processes involved, while ensuring that the time variation of the errors generated by each sensor can be repeated if so desired.Finally, the primary use of the sensor models developed in this article is to be part of a high fidelity flight simulation such as that described in [2], so the models need to be realistic and include as few simplifications as possible.
Section 2 describes the objectives of the simulation with special emphasis on its stochastic properties and how it can be customized for different sensors.It is followed by detailed explanations of the models representing the errors present in the measurements of the different sensors.Section 3 focuses on the inertial sensors and explains their characteristics and error sources, resulting in a model that captures their stochastic properties.A similar but less detailed approach is taken in section 4 with the magnetometers, air data system, and GNSS receiver.
Section 5 describes the characteristics of the digital camera employed for visual navigation, and the tool employed to generate realistic images that resemble what a real camera would view if located at the same position and attitude; although the camera differs from all other sensors in that it does not provide a measurement or reading but a digital image, in this article it is indeed considered a sensor as it provides the GNC system with information about its surroundings that can be employed for navigation.Last, section 6 describes sensor calibration activities that are indispensable for the determination of various of the parameters present in the sensor models.

Customizable and Stochastic Simulation
The different sensor models described in this article have been implemented in a C++ simulation (the open-source code is available at [1]) with the objective of properly simulating their stochastic nature while ensuring that the simulated variation with time of the errors can be repeated at a later time if so desired.The implementation hence follows these steps: 1. Initialize a discrete uniform distribution with any seed (any value is valid, so 1 was employed by the author), which produces random integers where each possible value has an equal likelihood of being produced.Call this distribution a number of times equal or higher than twice the maximum number of runs to be executed (each run provides the variation of time of all sensors for an unlimited amount of time, and corresponds to a single aircraft flight), divide them into two groups of the same size, and store the results for later use.
These values, called Υ i,A and Υ j,F , are respectively the aircraft seeds and the flight seeds, where "i" is the aircraft number representing given fixed error realizations 7 , and "j" is the trajectory number representing run-to-run and in-run error realizations 8 .The stored aircraft and flight seeds become the initialization seeds for each of the simulation executions or runs, so this step does not need to be repeated.
2. Every time the simulator needs to obtain the errors generated by the different sensors (which usually correspond to a given flight), it is initialized with a given aircraft seed Υ i,A together with a flight seed Υ j,F .
As these seeds are the only inputs required for all the stochastic processes within the sensors, the results of a given run can always be repeated by employing the same seeds.

Type
Error Sources Seeds Aircraft i fixed υ i,A,ACC , υ i,A,GYR , υ i,A,MAG , υ i,A,PLAT , υ i,A,CAM Flight j run-to-run & in-run υ j,F,ACC , υ j,F,GYR , υ j,F,MAG , υ j,F,OSP , υ j,F,OAT υ j,F,TAS , υ j,F,AOA , υ j,F,AOS , υ j,F,GNSS The selected aircraft and flight seeds are then employed to initialize two different discrete uniform distributions.One is executed five times to provide the fixed sensor seeds (υ i,A,ACC , υ i,A,GYR , υ i,A,MAG , υ i,A,PLAT , and υ i,A,CAM ), while the other is realized nine times to obtain the run sensor seeds (υ j,F,ACC , υ j,F,GYR , υ j,F,MAG , 7 Fixed error contributions vary from aircraft to aircraft but are constant for all flights of that aircraft. 8Run-to-run and in-run error contributions vary from one flight of a given aircraft to the next. υ j,F,OSP , υ j,F,OAT , υ j,F,TAS , υ j,F,AOA , υ j,F,AOS , and υ j,F,GNSS ).These seeds hence become the initialization seeds for each of the different sensors described throughout this article.
3. Each sensor relies on either one or two standard normal distributions N (0, 1), depending on whether its error model is based exclusively on run-to-run and in-run error contributions or it also contains fixed error sources.The normal distributions of every sensor are initialized with the corresponding seeds (υ i,F and υ j,R ) for that sensor.
4. Upon initialization, the fixed normal distribution of every sensor is employed to generate all the values corresponding to scale factors, cross couplings, hard iron magnetism, and mounting errors.The run normal distribution in turn is employed to generate the required bias offsets.

5.
Once the simulation has been initialized, it is able to estimate the errors generated by each sensor working at the required sensor rate.As the simulation time advances, every time a sensor is called to provide a measurement, its already initialized and used run normal distribution is called to generate the corresponding random walk increments and white noises.
6. Additionally, when provided with a time series of the aircraft position and attitude, the simulation generates realistic images of the Earth surface that resemble those taken with an onboard camera.
All the error models described in the following sections rely on a series of parameters whose value varies depending on the specific sensors represented in the simulation.Even though the user is expected to specify the value of these parameters, this article suggests the values used by the author to represent the sensors that could nowadays be installed onboard a small size fixed wing low SWaP autonomous aircraft, and which from this point on are referred as the "Baseline" configuration.

Inertial Sensors
Inertial sensors comprise accelerometers and gyroscopes, which measure specific force and inertial angular velocity about a single axis, respectively [6].An inertial measurement unit (IMU) encompasses multiple accelerometers and gyroscopes, usually three of each, obtaining three dimensional measurements of the specific force and angular rate [3] viewed in the platform frame F P (section 3.5).The individual accelerometers and gyroscopes however are not aligned with the F P axes, but with those of the non orthogonal accelerometers F A and gyroscopes F Y frames, also defined in section 3.5.The output of the inertial sensors must hence first be transformed from the F A and F Y frames to F P , as described in sections 3.6 and 3.7, and then from the F P frame to the body frame F B as explained in section 3.8, where they can be employed by the navigation system.The accelerometers and gyroscopes are assumed to be infinitesimally small and located at the IMU reference point (section 3.8), which coincides with the origin of these three frames The IMU is physically attached to the aircraft structure in a strapdown configuration, so both the displacement T B BP and the Euler angles φ BP = [ψ P , θ P , ξ P ] T that describe the relative position and rotation between the body F B and platform F P frames are constant.Accelerometers can be divided by their underlying technology into pendulous and vibrating beams, while gyroscopes are classified into spinning mass, optical (ring laser or fiber optic), and vibratory [5].Current inertial sensor development is mostly focused on micro machined electromechanical system (MEMS) sensors 9 , which makes direct use of the chemical etching and batch processing techniques used by the electronics integrated circuit industry to obtain sensors with small size, low weight, rugged construction, low power consumption, low price, high reliability, and low maintenance [7].On the negative side, the accuracy of MEMS sensors is still low, although tremendous progress has been achieved in the last two decades and more is expected in the future.
There is no universal classification of inertial sensors according to their performance, although they can be broadly assigned into five different categories or grades: marine (submarines and spacecraft), aviation (commercial and military), intermediate (small aircraft and helicopters), tactical (unmanned air vehicles and guided weapons), and automotive (consumer) [5].The full range of grades covers approximately six orders of magnitude of gyroscope performance and only three for the accelerometers, but higher performance is always associated with bigger size, weight, and cost.Tactical grade IMUs cover a wide range of performance values, but can only provide a standalone navigation solution for a few minutes, while automotive grade IMUs are unsuitable for navigation.
The different errors that appear in the measurements provided by accelerometers and gyroscopes are described in section 3.1.Section 3.2 presents a model for the measurements of a single inertial sensor, while sections 3.3 and 3.4 focus on how to obtain the specific values for white noise and bias on which the model relies from the documentation.Section 3.5 describes the reference systems required to represent the IMU measurements.
Additional errors appear when three accelerometers or gyroscopes are employed together, and these are modeled in section 3.6 for accelerometers and section 3.7 for gyroscopes.The analysis of the inertial sensors concludes with sections 3.9 and 3.10, which provide a comprehensive error model for the IMU measurements, as well as the specifications of the baseline configuration.The final model also depends on the relative position of the IMU with respect to the body frame, which is described in section 3.8.

Inertial Sensor Error Sources
In addition to the accelerometers and gyroscopes, an IMU also contains a processor, storage for the calibration parameters, one or more temperature sensors, and a power supply.As described below, each sensor has several error sources, but each of them has four components: fixed contribution, temperature dependent variation, runto-run variation, and in-run variation [5,8].The first two can be measured at the laboratory (at different temperatures) and the calibration results stored in the IMU so the processor can later compensate the sensor outputs based on the reading provided by the temperature sensor.Calibration however increases manufacturing costs so it may be absent in the case of inexpensive sensors.The run-to-run variation results in a contribution to a given error source that varies every time the sensor is employed but remains constant within a given run.It can not be compensated by the IMU processor but can be calibrated by the navigation system every time it is turned on with a process known as fine-alignment [3,5,9].The in-run contribution to the error sources slowly varies during execution and can not be calibrated in the laboratory nor by the navigation system.Let's now discuss the different sources of error that influence an inertial sensor [5,10,11,12]:

IMU Grade Accelerometer Bias
• The bias is an error exhibited by all accelerometers and gyroscopes that is independent of the underlying specific force or angular rate being measured, and comprises the dominant contribution to the overall sensor error.It can be defined as any nonzero output when the sensor input is zero [8], and can be divided into its static and dynamic components.The static part, also known as fixed bias, bias offset, turn-on bias, or bias repeatability, comprises the run-to-run variation, while the dynamic component, known as in-run bias variation, bias drift, bias instability (or stability), is typically about ten percent of the static part and slowly varies over periods of order one minute.As the bias is the main contributor to the overall sensor error, its value can be understood as a sensor quality measure.Table 4 provides approximate values for the inertial sensor biases according to the IMU grade [5]: While the bias offset can be greatly reduced through fine-alignment [3,5,9], the bias drift can not be determined and needs to be modeled as a stochastic process.It is mostly a warm up effect that should be almost non existent after a few minutes of operation, and corresponds to the minimum point of the sensor Allan curve [11,13].It is generally modeled as as a random walk process obtained by the integration of a white noise signal, coupled with limits that represent the conclusion of the warm up process.
• The scale factor error is the departure of the input output gradient of the instrument from unity following unit conversion at the IMU processor.It represents a varying relationship between sensor input and output caused by aging and manufacturing tolerances.As it is a combination of fixed contribution plus temperature dependent variation, most of it can be eliminated through calibration (section 6.1).
• The cross coupling error or non orthogonality error is a fixed contribution that arises from the misalignment of the sensitive axes of the inertial sensors with respect to the orthogonal axes of the platform frame due to manufacturing limitations, and can also be highly reduced through calibration.The scale factor and cross coupling errors errors are in the order of 10 −4 and 10 −3 [−] for most inertial sensors, although they can be higher for some low grade gyroscopes.The cross coupling error is equal to the sine of the misalignment, which is listed by some manufacturers.
• System noise or random noise is inherent to all inertial sensors and can combine electrical, mechanical, resonance, and quantization sources.It can originate at the sensor itself or at any other electronic equipment that interferes with it.System noise is a stochastic process usually modeled as white noise because its noise spectrum is approximately white, and can not be calibrated as there is no correlation between past and future values.A white noise process is characterized by its spectral power density PSD, which is constant as it does not depend on the signal frequency.It corresponds to the 1 [sec] crossing of the sensor Allan curve [11,13].
System noise is sometimes referred to as random walk, which can generate confusion with the bias.The reason is that the inertial sensor outputs are always integrated to obtain ground velocity in the case of accelerometers and aircraft attitude in the case of gyroscopes.As the integration of a white noise process is indeed a random walk, the later term is commonly employed to refer to system noise.and vibratory gyroscopes to specific force), scale factor nonlinearity, and higher order errors (spinning mass gyroscopes and pendulous accelerometers).

Single Axis Inertial Sensor Error Model
Since all sensors modeled in this article are required to provide measurements at equispaced discrete times t s = s • ∆t SENSED = s • ∆t, this section focuses on obtaining a discrete model for the bias and white noise errors of a single axis inertial sensor.The results obtained here will be employed in the following sections to generate a comprehensive IMU model.
Let's consider a sensor in which the difference between its measurement at any given time x (t) and the real value of the physical magnitude being measured at that same time x (t) can be represented by a zero mean white noise Gaussian process η v (t) with spectral density σ 2 v : Dividing (2) by ∆t SENSED = ∆t and integrating results in: Assuming that the measurement and real value are both constant over the integration interval10 [14] yields This expression results in the white noise sensor error w (t), which is the difference between the sensor measurement x (t) and the true value x (t).Its mean and variance can be readily computed: E [w (t 0 + ∆t)] = 0 (6) Based on these results, the white noise error can be modeled by a discrete random variable identically distributed to the above continuous white noise error, this is, one that results in the same mean and variance, where ) is a standard normal random variable: Let's now consider a second model in which the measurement error or bias is given by a first order random walk process or integration of a zero mean white noise Gaussian process η u (t) with spectral density σ 2 u : Its mean and variance can be quickly computed: These results indicate that the bias can be modeled by a discrete random variable identically distributed (this is, resulting in the same expected value and variance) to the continuous random walk above: where N u ∼ N (0, 1) is a standard normal random variable.Operating with the above expression results in the final expression for the discrete bias, as well as its mean and variance: A comprehensive single axis sensor error model without scale factor can hence be constructed by adding together the influence of the system noise provided by (8) and the bias given by ( 13) [15], while assuming that the standard normal random variables N u and N v are uncorrelated 11 : E [e BW (s ∆t)] = 0 Var The discrete sensor error or difference between the measurement provided by the sensor x (s ∆t SENSED ) = x (s ∆t) at any given discrete time s ∆t and the real value of the physical variable being measured at that same discrete time x (s ∆t) is the combination of a bias or first order random walk and a white noise process, and depends on three parameters: the bias offset B 0 , the bias instability σ u , and the white noise σ v .The contributions of these three different sources to the sensor error as well as to its first and second integrals 12 are very different and inherent to many of the challenges encountered when employing accelerometers and gyroscopes for inertial navigation, as explained below.with time of the sensor error mean (left side) and standard deviation (right side) given by ( 17) and ( 18) together with the average of fifty different runs.In addition, the left figure also includes ten of those runs to showcase the variability in results implicit to the random variables 13 , while the right side shows the theoretical contribution to the standard deviation of each of the three components.In addition to the near equivalence between the theory and the average of several runs, the figure shows that the bias instability is the commanding long term factor in the deviation between the sensor measurement and its zero mean (the standard deviation of the bias instability grows with the square root of time while the other two components are constant).As discussed in section 3.1, the bias drift or bias instability is indeed the most important quality parameter of an inertial sensor.This is also the case when the sensor output is integrated, as discussed below.
Let's integrate the sensor error over a timespan s ∆t to evaluate the growth with time of both its expected value combination of the two continuous error models given by ( 5) and ( 9).This is the case even if considering that the two zero mean white noise Gaussian processes ηu and ηv are independent and hence uncorrelated.It is however possible to obtain a discrete model whose discrete bias and white noise components are not only identically distributed that those of their continuous counterparts [14], even adding the equivalence of covariance between the bias and the sensor error, but this results in a significantly more complex model that behaves similarly to the one above at all but the shortest time samples after sensor initialization.The author has decided not to do so in the model described in this article, reducing complexity with little or no loss of realism.
12 Gyroscopes measure angular velocity and their output needs to be integrated once to obtain attitude, while accelerometers measure specific force and are integrated once to obtain velocity and twice to obtain position. 13Although the data is generated at 100 [hz], the figure only employs 1 out of every 1000 points, so it appears far less noisy than the real data.
and its variance 14 :  has been employed in the experiment), but the growth in the standard deviation (over a non zero mean) is much quicker than before.The bias instability continues to be the dominating factor but now increases with a power of t 3/2 , while the bias offset and white noise contributions also increase with time, although with powers of t and t 1/2 respectively.Let's continue the process and integrate the error a second time: ) Figure 4 shows the same type of figures but applied to the second integral of the error (g BW (0) = 1.5 has been employed in the experiment).In this case the degradation of the results with time is even more intense to the point were the measurements are useless after a very short period of time.Unless corrected by the navigation system, 14 As the interest lies primarily in s 1, a simple integration method such as the rectangular rule is employed.
this is equivalent to the error in position obtained by double integrating the output of the accelerometers.Let's summarize the main points of the single axis inertial sensor discrete error model developed in this section, that includes the influence of the bias and the system error but not that of the scale factor and cross coupling errors included in the three dimensional error model of section 3.9.The error e BW (s ∆t), which applies to specific force for accelerometers and inertial angular velocity in the case of gyroscopes, depends on three factors: bias offset B 0 , bias drift σ u , and white noise σ v .Its mean is always zero, but the error standard deviation grows with time (∝ t 1/2 ) due to the bias drift with constant contributions from the bias offset and the white noise.
When integrating the error to obtain f BW (s ∆t), equivalent to ground velocity for accelerometers and attitude for gyroscopes, the initial speed error or initial attitude error f BW (0) becomes the fourth contributor, and an important one indeed, as it becomes the mean of the first integral error at any time.The standard deviation, which measures the spread over the non zero mean, increases very quickly with time because of the bias instability (∝ t 3/2 ), with contributions from the offset (∝ t) and the white noise (∝ t 1/2 ).If integrating a second time to obtain g BW (s ∆t), or position in case of the accelerometer, the initial position error g BW (0) turns into the fifth contributor.The expected value of the position error grows linearly with time due to the initial velocity error with an additional constant contribution from the initial position error, while the position standard deviation (measuring spread over a growing average value) grows extremely quick due mostly to the bias instability (∝ t 5/2 ), but also because of the bias offset (∝ t 2 ) and the white noise (∝ t 3/2 ).Table 6 shows the standard units of the different sources of error for both accelerometers and gyroscopes.

Obtainment of System Noise Values
This section focuses on the significance of system or white noise error σ v and how to obtain it from sensor specifications, which often refer to the integral of the output instead of the output itself.As the integral of white noise is a random walk process, the angle random walk of a gyroscope is equivalent to white noise in the angular rate output, while velocity random walk refers to specific force white noise in accelerometers [8].The discussion that follows applies to gyroscopes but is fully applicable to accelerometers if replacing angular rate by specific force and attitude or angle by ground velocity.
Angle random walk, measured in [rad/sec 1/2 ], [deg/hr 1/2 ], or equivalent, describes the average deviation or error that occurs when the sensor output signal is integrated due to system noise exclusively, without considering other error sources such as bias or scale factor [16].If integrating multiple times and obtaining a distribution of end points at a given final time s ∆t, the standard deviation of this distribution, containing the final angles at the final time, scales linearly with the white noise level σ v , the square root of the integration step size ∆t, and the square root of the number of steps s, as noted by the last term of (21).This means that an angle random walk of 1 [deg/sec

Obtainment of Bias Drift Values
This section describes the meaning of bias instability σ u (also known as bias stability or bias drift) and how to obtain it from sensor specifications.As in the previous section, the discussion is centered on gyroscopes but it is fully applicable to accelerometers as well.Bias instability can be defined as the potential of the sensor error to stay within a certain range for a certain time [17].A small number of manufacturers directly provide sensor output changes over time, which directly relates with the bias instability (also known as in-run bias variation, bias drift, or rate random walk) per the second term of (18).If provided with an angular rate change of x [deg/sec] (1σ) in t [sec], then σ u can be obtained as follows [14,18]: As the bias drift is responsible for the growth of sensor error with time (figure 2), manufacturers more commonly provide bias stability measurements that describe how the bias of a device may change over a specified period of time [15] 25) can be used to obtain σ u .As the bias behaves as a random walk over time whose standard deviation grows proportionally to the square root of time, the bias stability is sometimes referred as a bias random walk.
In reality bias fluctuations do not really behave as a random walk.If they did, the uncertainty in the bias of a device would grow without bound as the timespan increased, which is not the case.In practice the bias is constrained to be within some range, and therefore the random walk model is only a good approximation to the true process for short periods of time [15].

Platform, Accelerometers, and Gyroscopes Frames
The following sections make use of three different reference frames to describe the readings of accelerometers and gyroscopes: • The platform frame F P is a Cartesian reference system with its origin located at the IMU reference point (section 3.8) and its three axes {i P 1 , i P 2 , i P 3 } forming a right hand system that is loosely aligned with the aircraft body axes, so they point in the general directions of the aircraft fuselage (forward), aircraft wings (rightwards), and downward, respectively [3,9,19].
A proper definition of the platform frame is indispensable for navigation, as the calibrated outputs of the accelerometers and gyroscopes are based on it (sections 3.6 and 3.7).The F P frame can be obtained from the body frame F B by a rotation best described by the Euler angles φ BP = [ψ P , θ P , ξ P ] T15 followed by a translation T B BP (section 3.8) from the aircraft center of mass to the IMU reference point.• The accelerometers frame F A is a non orthogonal reference system also centered at the IMU reference point [3,9,19].The basis vectors {i A 1 , i A 2 , i A 3 } are aligned with each of the three accelerometers sensing axes16 (section 3.6), but they are not orthogonal among them due to manufacturing inaccuracies.This implies that the angles between the F A and F P axes are very small.
It is always possible, with no loss of generality, to impose that i P 1 coincides with i A 1 and that i P 2 is located in the plane defined by i A 1 and i A 2 .If this is the case, i A 1 ⊥ i P 2 , i A 1 ⊥ i P 3 , and i A 2 ⊥ i P 3 , and the relative attitude between the F P and F A frames can be defined by three independent small rotations.
-The i A 2 axis can be obtained from i P 2 by means of a small rotation α ACC,3 about i P 3 .
-The i A 3 axis can be obtained from i P 3 by two small rotations: α ACC,1 about i P 1 and α ACC,2 about i P 2 .Although the exact relationships can be obtained [9], and given that the angles are very small, it is possible to consider cos α ACC,i = 1, sin α ACC,i = α ACC,i , and α ACC,i • α ACC,j = 0 ∀ i, j ∈ {1, 2, 3} , i = j, resulting in the following transformations between free vectors viewed in the platform (v P ) and accelerometer (v A ) frames respectively 17 : • The gyroscopes frame F Y is similar to the accelerometers frame F A defined above, but aligned with the gyroscopes sensing axes instead of those of the accelerometers [3,9,19].It is also a non orthogonal reference system centered at the IMU reference point, but no simplifications can be made about the relative orientation of its axes {i Y 1 , i Y 2 , i Y 3 } with respect to those of F P , so their relative attitude is defined by six small rotations α GYR,ij ∀ i, j ∈ {1, 2, 3} , i = j, where α GYR,ij is the rotation of i Y i about i P j .
An approach similar to that employed with accelerometers leads to the following transformations between free vectors viewed in the platform (v P ) and gyroscope (v Y ) frames:

Accelerometer Triad Sensor Error Model
An IMU is equipped with an accelerometer triad composed by three individual accelerometers, each of which measures the projection of the specific force over its sensing axis as described in section 3.2 while incurring in an error e BW,ACC that can be modeled as a combination of bias offset, bias drift, and white noise (16).The three accelerometers can be considered infinitesimally small and located at the IMU reference point, which is defined as the intersection between the sensing axes of the three accelerometers.As the accelerometer frame F A is centered at the IMU reference point and its three non orthogonal axes coincide with the accelerometers sensing axes, (30) joins together the measurements of the three individual accelerometers: where f A IA is the specific force viewed in the accelerometer frame F A , f A IA represents its measurement also viewed in F A , e A BW,ACC is the error introduced by each accelerometer (16), and S ACC is a square diagonal matrix containing the scale factor errors {s ACC,1 , s ACC,2 , s ACC,3 } for each accelerometer (section 3.1).It is however preferred to obtain an expression in which the specific forces are viewed in the orthogonal platform frame F P instead of the accelerometer frame F A .As both share the same origin, where R PA and R AP , defined by ( 26) and ( 27), contain the cross coupling errors {α ACC,1 , α ACC,2 , α ACC,3 } generated by the misalignment of the accelerometer sensing axes.The scale factor and cross coupling errors contain fixed and temperature dependent error contributions (refer to section 3.1) that can be modeled as normal random variables: where s ACC and α ACC can be obtained from the sensor specifications.Equation ( 31) can be transformed to make it more useful by defining the accelerometer scale and cross coupling error matrix M ACC : Considering that the scale and cross coupling errors are uncorrelated and very small, and taking into account the expressions for the mean and variance of the sum and product of two random variables [20], the different components m ACC,ij of M ACC can be obtained as follows ∀ i, j ∈ {1, 2, 3}: Let's also define the accelerometer error transformation matrix N ACC as A process similar to that employed above leads to: Taking into account the expressions for the mean and variance of the sum and product of two random variables [20], and knowing that the cross coupling errors are very small 1 + α 2 ACC ≈ 1 , it can be proven that the bias and white noise errors viewed in the platform frame F P respond to a expression similar to ( 16): e P BW,ACC = e P BW,ACC (s ∆t SENSED ) = e P BW,ACC (s ∆t) where each N u,ACC and N v,ACC is a random vector composed by three independent standard normal random variables N (0, 1).Note that as the bias drift is mostly a warm up process that stabilizes itself after a few minutes of operation, the random walk within (42) is not allowed to vary freely but is restricted to within a band of width ± 100 σ uACC ∆t 1/2 .The final model for the accelerometer measurements viewed in F P results in where M ACC is described in ( 34) through ( 37) and e P BW,ACC is provided by (42).This model relies on inputs for the bias offset B 0ACC , bias drift σ uACC , white noise σ vACC , scale factor error s ACC , and cross coupling error m ACC , which are obtained from the accelerometer specifications shown in table 10 within section 3.10.

Gyroscopes Triad Sensor Error Model
The IMU is also equipped with a triad of gyroscopes, each of which measures the projection of the inertial angular velocity over its sensing axis as described in section 3.2.The obtainment of the gyroscope triad model is fully analogous to that of the accelerometers in the previous section, with the added difficulty that the transformation between the gyroscope frame F Y and platform frame F P relies on six small angles instead of three.The starting point hence is: where ω P IP is the inertial angular velocity viewed in the platform frame F P , ω P IP represents its measurement also viewed in F P , e Y BW,GYR is the error introduced by each gyroscope (16), S GYR is a square diagonal matrix containing the scale factor errors {s GYR,1 , s GYR,2 , s GYR,2 }, and R PY and R YP , defined by ( 28) and ( 29), contain the cross coupling errors α GYR,12 , α GYR,21 , α GYR,13 , α GYR,31 , α GYR,23 , α GYR,32 generated by the misalignment of the gyroscope sensing axes.
Operating in the same way as in section 3.6 leads to: e P BW,GYR = e P BW,GYR (s ∆t SENSED ) = e P BW,GYR (s ∆t) where each N ui,GYR and N v,GYR is a random vector composed by three independent standard normal random variables N (0, 1).As in the case of the accelerometers, the random walk within (45) representing the bias drift is not allowed to vary freely but is restricted to within a band of width ± 100 σ uGYR ∆t 1/2 .This model relies on inputs for the bias offset B 0GYR , bias drift σ uGYR , white noise σ vGYR , scale factor error s GYR , and cross coupling error m GYR , which are obtained from the gyroscope specifications shown in table 9 within section 3.10.The gyroscope scale and cross coupling error matrix M GYR responds to:

Mounting of Inertial Sensors
Equations ( 43) and ( 46) contain the relationships between the specific force f P IP and inertial angular velocity ω P IP and their measurements f P IP , ω P IP , when evaluated and viewed in the platform frame F P .However, from the point of view of the navigation system, both magnitudes need to be evaluated and viewed in the body frame F B instead of F P .These equations thus need to be modified so they relate f B IB with f B IB as well as ω B IB with ω B IB , respectively, as described in section 3.9 below.To do that, it is necessary to define the relative pose (position plus attitude) between the F P and F B frames.Note that the IMU, represented by the platform frame F P , should be mounted in the aircraft as close as possible to the center of gravity (this reduces errors, as described in section 3.9), and loosely aligned with the aircraft body axes.
To provide realism to the stochastic simulation, this article assumes that the real displacement T B BP between the two frames is deterministic, while the relative rotation φ BP = [ψ P , θ P , ξ P ] T is stochastic.In this way, each simulated aircraft represented by its specific aircraft seed Υ i,A exhibits a slightly different IMU platform attitude with respect to the aircraft body: • As the IMU reference point is fixed with respect to the structure but the aircraft center of mass slowly moves as the fuel load diminishes, it is possible to establish a linear model that provides the displacement between the origins of both frames according to the aircraft mass 18 : • The platform Euler angles respond to the stochastic model provided by (51), in which each specific Euler angle is obtained as the product of the standard deviations (σ ψ P , σ θ P , σ ξ P ) by a single realization of a standard normal random variable N (0, 1) (N ψ P , N θ P , and N ξ P ).
φ BP = σ ψ P N ψ P , σ θ P N θ P , σ are six realizations of a standard normal random variable N (0, 1).Note that table 7 lists the default values employed in the simulation, which can be adjusted by the user.
T BP can be considered quasi stationary as it slowly varies based on the aircraft mass, and the relative position of their axes φ BP remains constant because the IMU is rigidly attached to the aircraft structure.Although Euler angles have been employed in this section, from this point on it is more practical to employ the rotation matrix R BP to represent the rotation between two different frames [21].The time derivatives of T BP and R BP are hence zero:

Comprehensive Inertial Sensor Error Model
Two considerations are required to establish the measurement equations for the inertial sensors viewed in the body frame F B .First, let's apply the composition rules of section A.6 considering F I as F 0 , F B as F 1 , and F P as F 2 , which results in: Second, it is also necessary to consider that as RBP is a rotation matrix in which all rows and columns are unitary vectors, the projection of the F P frame bias and white noise errors e P BW,ACC and e P BW,GYR onto the F B frame does not change their stochastic properties: As the inertial angular velocity does not change when evaluated in the F B and F P frames (55), its measurement in the body frame can be derived from (46) by first projecting it from F B to F P based on the real rotation matrix R BP and then projecting back the measurement into F B based on the estimated rotation matrix RBP .The bias and white noise error is also projected according to (60): The expression for the specific force measurement is significantly more complex because the back and forth transformations of the specific force between the F B and F P frames need to consider the influence of the lever arm T BP , as indicated in (58).The additional terms introduce errors in the measurements, so as indicated in section 3.8 it is desirable to locate the IMU as close as possible to the aircraft center of mass.
Note that this expression can not be directly evaluated as the estimated values for the inertial angular velocity and acceleration ( ωB IB , αB IB ) are unknown by the IMU until obtained by the navigation filter.The IMU can however rely on the gyroscope readings, directly replacing ωB IB with ω B IB and computing α B IB based on the difference between the present and previous ω B IB readings, resulting in:  T BP displacement between the body and platform frames, which is deterministic), as expressions (61, 63) rely on the errors provided by (59, 60), the scale and cross coupling matrices given by (34,47), and the transformations given by (50, 51, 52, 53).
In the case of the accelerometer triad, the stochastic nature of the fixed and run-to-run error contributions to the model relies on three realizations of normal distributions for the bias offset, three for the scale factor errors, three for the cross coupling errors, and nine for the mounting errors, while the in-run error contributions require three realizations each for the bias drift and system noise at every discrete sensor measurement.The gyroscope triad is similar but requires six realizations to model the cross coupling errors instead of three, while using the same six realizations that the accelerometer triad to model the true and estimated rotation between the F B and Expressions (61, 63) can be rewritten to show the measurements as functions of the full errors (E ACC , E GYR ), which represent all the errors introduced by the inertial sensors with the exception of white noise.

Baseline Inertial Sensors
The value of all parameters present in the previous models can be specified by the user when employing the simulation [1,2].This article however provides default values that represent the inertial sensors that could nowadays be installed onboard a small size fixed wing low SWaP autonomous aircraft, and which are referred to as the baseline configuration.In the case of the gyroscopes, they correspond to the MEMS gyroscopes installed inside the Analog Devices ADIS16488A IMU [22].Table 9 shows its performances, which have been taken from the data sheet when possible, and corrected when suspicious.A calibration process as that described in section 6.1 is assumed to eliminate 95% of the scale factor and cross coupling errors.The previous table show the performances in three columns of data.The left most column ("Spec") corresponds to data taken directly from the sensors specifications, which get converted in the middle column ("Value") to the parameters shown in sections 3.2 through 3.9 20 .The right column ("Calibration") contains the specifications employed in the simulation after the calibration process, which reduces the scale factor and cross coupling errors by 95%.The baseline accelerometers also correspond to the MEMS accelerometers installed inside the Analog Devices ADIS16488A IMU [22].All the specifications shown in table 10 have been taken from the data sheet.As in the case of the gyroscopes, a calibration process as that described in section 6.1 is assumed to eliminate 95% of the scale factor and cross coupling errors 21 .

Non Inertial Sensors
This section describes the different non-inertial sensors usually installed onboard a fixed wind autonomous aircraft, such as a triad of magnetometers to measure the Earth magnetic field, a GNSS receiver that provides absolute position and velocity measurements, and the air data system, which in addition to the pressure altitude and temperature also provides a measurement of the airspeed and the airflow angles.

Magnetometers
Magnetometers measure magnetic field intensity along a given direction and are very useful for estimating the aircraft heading.Although other types exist, magnetoinductive and magnetoresistive sensors are generally employed for navigation due to their accuracy and small size [5,7].As with the inertial sensors, three orthogonal magnetometers are usually employed in a strapdown configuration to measure the magnetic field with respect to the body frame F B .
Unfortunately magnetometers do not only measure the Earth magnetic field B, but also that generated by the aircraft permanent magnets and electrical equipment (known as hard iron magnetism), as well as the magnetic field disturbances generated by the aircraft ferrous materials (soft iron magnetism).For that reason, the magnetometers should be placed in a location inside the aircraft that minimizes these errors.On the positive side, magnetometers do not exhibit the bias instability present in inertial sensors, and the error of an individual sensor can be properly modeled by the combination of bias offset and white noise.A triad of magnetometers capable of measuring the magnetic field in three directions adds the same scale factor (nonlinearity) and cross coupling (misalignment) errors as those present in the inertial sensors, together with the transformation between the magnetic axes and the body ones.
Modeling the behavior of a triad of magnetometers is simpler but less precise than that of inertial sensors, as the effect of the fixed hard iron magnetism is indistinguishable from that of the run-to-run bias offset, while the fixed effect of soft iron magnetism is indistinguishable from that of the scale factor and cross coupling error matrix.This has several consequences.First of all is that magnetometers can not be calibrated at the laboratory before being mounted in the aircraft as in the case of inertial sensors (section 6.1), but are instead calibrated once attached to the aircraft by a process known as swinging (section 6.3), which is less precise as the aircraft attitude during swinging can not be determined with so much accuracy as it would be in a laboratory setting.Second is that defining a magnetic platform frame to then transform the results into body axes serves no purpose, as the magnetometer readings are only valid, this is, contain the effects of hard and soft iron magnetism, once they are attached to the aircraft, and then they can be directly measured in body axes.And third is that percentage wise the errors induced by the magnetometers are bigger than those of the inertial sensors.The simulation relies on the following model: where B B is the measurement viewed in F B , B 0,MAG is the run-to-run bias offset, B HI,MAG is the fixed hard iron magnetism, M MAG is a fixed matrix combining the effects of soft iron magnetism with the scale factor and cross coupling errors, and B N is the real magnetic field including local anomalies.The baseline magnetometer features employed in the simulation are shown in table 11 in the same format as those of the inertial sensors, where the white noise has been taken from [5] and the rest of the parameters correspond to the magnetometers present in the Analog Devices ADIS16488A IMU [22].Although the value of hard and soft iron magnetism in aircraft is rather small, the author has not been able to obtain trusted values for them.

MAG Baseline
To avoid eliminating sources of error, the author has decided to increase by 50% the values for bias offset, scale factor, and cross coupling errors found in the literature, as shown in the column named "Compensation".As both result in a similar effect, the author expects that the realism of the results will not be adversely affected.
In the case of the bias, the author has assigned most of the error to the fixed hard iron error B HI,MAG and the remaining to the run-to-run bias offset B 0,MAG .
A swinging process as that described in section 6.3 is assumed to eliminate 90% of the fixed error contributions, this is, the hard iron magnetism, the scale factor, and the cross coupling error 22 .This number is inferior to the 95% reduction achieved by the calibration of inertial sensors because although the lack of magnetometer bias drift facilitates calibration, the determination of the body attitude during calibration is inherently less accurate working with the complete aircraft than just with the IMU in a laboratory setting.The final results can be found in the rightmost column above.
The final model for a triaxial magnetometer thus includes contributions from the hard iron magnetism, bias offset, system noise, soft iron magnetism, scale factor, and cross coupling errors: where N HI,MAG , N u0,MAG and N vs,MAG are uncorrelated normal vectors of size three each composed of three uncorrelated standard normal random variables N (0, 1).The soft iron, scale factor and cross coupling matrix M MAG does not vary with time and is computed as follows: In this expression N m,MAG contains nine outputs of a standard normal random variable N (0, 1), and the symbol • represents the Hadamart or element-wise matrix product.
Table 12 lists the error sources contained in the magnetometer model represented by (67), noting that soft iron magnetism is also included in both the scale factor and cross coupling errors.The first two columns list the different error sources, while the third column specifies their origin according to the criterion established in the first paragraph of section 3.1.The last two columns list the seeds employed to ensure the results variability for different aircraft (Υ i,A ) as well as different flights (Υ j,F ): and that all errors are modeled as stochastic variables or processes.The stochastic nature of the fixed and run-to-run error contributions to the magnetometer model relies on three realizations of normal distributions for the hard iron magnetism, three for the bias offset, three for the scale factor errors, and six for the cross coupling errors, while the in-run error contributions require three realizations for system noise at every discrete sensor measurement.
Expression (67) can be rewritten to show the measurements as functions of the magnetometer full error E MAG , which represents all the errors introduced by the magnetometers with the exception of white noise: It is important to note that because of the heavy processing required to fix a position based on the satellite signals, GNSS receivers are not capable of working at the high frequencies characteristic of inertial and air data sensors, with 1 [hz] being the default value in this article (∆t GNSS = 1 [sec]).The position error of a GNSS receiver can be modeled as the sum of a zero mean white noise process plus slow varying ionospheric effects [23] modeled as the sum of the bias offset plus a random walk.This random walk is modeled with a frequency of 1/60 [hz]

Global Navigation Satellite System Receiver
(∆t ION = 60 [sec]), and linearly interpolated in between.The ground velocity error is modeled exclusively with a white noise process.
e GNSS,POS (g ∆t GNSS ) = x GDT − x GDT = σ GNSS,POS N g,GNSS,POS + e GNSS,ION (g ∆t GNSS ) (70) e GNSS,ION (g ∆t GNSS ) = e GNSS,ION (i ∆t ION ) + r f ION e GNSS,ION (i + 1) ∆t ION − e GNSS,ION (i ∆t ION ) (72) where σ GNSS,POS , σ GNSS,ION , B 0,GNSS,ION , and σ GNSS,VEL are taken from the table below, and N g,GNSS,POS , N g,GNSS,VEL , N u0,GNSS,ION , and N j,GNSS,ION and are uncorrelated normal vectors of size three each composed of three uncorrelated standard normal random variables N (0, 1).Also note that as both g and f ION are integers, the quotient remainder theorem guarantees that there exist unique integers i and r that comply with (73) [24].The horizontal position accuracy in the table above corresponds to the U-blox NEO-M8 receiver data sheet [25],

GNCC
where CEP stands for circular error probability.As CEP is equivalent to 1.18 standard deviations [26], it enables the obtainment of σ GNSS,POS,HOR .The author has not found any reference for GNSS vertical accuracy, but has determined through conversations with several knowledgeable individuals that it is at least 50% higher than the horizontal one.A conservative value for σ GNSS,POS,VER of twice that of σ GNSS,POS,HOR has been selected.The ionospheric effects have also been obtained from these conversations.The velocity accuracy also originates at the U-blox NEO-M8 receiver data sheet [25].Assuming that it corresponds to a per axis error of ±0.05 [m/sec] instead of CEP, and knowing that the 50% mark of a normal distribution lies at 0.67448 standard deviations, it is possible to obtain σ GNSS,VEL .Note that all errors are modeled as stochastic variables or processes.Three realizations of a normal distribution are required for the run-to-run error contributions, while the in-run error contributions require three realizations each for position and velocity at every discrete sensor measurement, plus an extra three when corresponding for the ionospheric error.

Air Data System
The mission of the air data system is to measure the aircraft pressure altitude H P [27,28] by means of the atmospheric pressure p, the outside air temperature T, the airspeed v TAS , and the angles of attack and sideslip that provide the orientation of the aircraft structure with respect to the airflow.
A barometer or static pressure sensor, generally part of the Pitot tube as explained below [7], measures atmospheric pressure, which can be directly translated into pressure altitude [27,28].The σ OSP value shown in table 15 below, where OSP stands for outside static pressure, comes from the ± 10 [m] altitude error listed in the specifications of the Aeroprobe air data system [29], which translates into ± 100 [pa] at a pressure altitude of 1500 [m].Although not present in the documentation, the author has decided to also include a bias offset B 0OSP for added realism, so the barometer error is modeled as a combination of bias offset and random noise: where N 0,OSP and N s,OSP are uncorrelated standard normal random variables N (0, 1).
Air data systems are also equipped with a thermometer to measure the external air temperature T. As in the case of the barometer above, the specifications of the Analog Devices ADT7420 temperature sensor [30]  Temperature Error (3σ) ± 0.15 Table 15: Performance of "Baseline" atmospheric sensors A Pitot probe is a tube with no outlet pointing directly into the undisturbed air stream, where the values of the air variables (temperature, pressure, and density) at its dead end resemble the total or stagnation variables of the atmosphere prior to its deceleration inside the Pitot [31].Total atmospheric variables (temperature T t , pressure p t , and density ρ t ) are those obtained if a moving fluid with static atmospheric variables (T, p, ρ) and speed v decelerates until it has no speed through a process that is stationary, has no viscosity nor gravity (gravitation plus inertial) accelerations, is adiabatic, and presents fixed boundaries for the analyzed control volume [32,33].
Such a process maintains the fluid total enthalpy, as well as its entropy, and hence complies with the Bernoulli equation [32]: where κ = 1.4 is the air adiabatic index.In addition to the static pressure and temperature sensors discussed above, a Pitot tube is also equipped with a dynamic pressure sensor located at the dead end to measure the air flow total pressure p t .The air data system then estimates the aircraft true airspeed based on the following expression, which results from applying (78) at the Pitot dead end as well as at the static ports: The errors induced when measuring the true airspeed this way can also be modeled by a combination of bias offset and system noise, where σ TAS is the error standard deviation taken from table 16, B 0TAS has been added by the author for increased realism, and N 0,TAS and N s,TAS are two uncorrelated standard normal random variables The Aeroprobe air data system specifications [29] list a maximum airspeed error of 1 [m/sec], which can be interpreted as 3 σ, and hence results in the σ TAS value shown below: Geometrically, the simulation adopts a perspective projection or pinhole camera model [34], which in addition is perfectly calibrated and hence shows no distortion.Table 18 contains the default camera parameters employed in the simulation, which can be modified by the user if so desired.

Mounting of Camera
The digital camera can be located anywhere on the aircraft structure as long as its view of the terrain is unobstructed by other platform elements.It is desirable that the lever arm or distance between the camera optical center and the aircraft center of mass is as small as possible to reduce the negative effects of any camera alignment error.With respect to its orientation, the camera should be facing down to show a balanced view of the ground during level flight, but minor deviations are not problematic.
As in the case of the IMU platform, the simulation considers that the camera location is deterministic but its orientation stochastic.The expressions below are hence analogous to those employed in section 3.
The problem however is not the real translation and rotation between the F B and F C frames given by the previous equations, but the accuracy with which they are known to the navigation system.are six realizations of a standard normal random variable N (0, 1).The translation T B BC between the origins of the F B and F C frames can be considered quasi stationary as it slowly varies based on the aircraft mass (83), and the relative position of their axes φ BC remains constant because the camera is rigidly attached to the aircraft structure (84).

Earth Viewer
The camera model differs from all other sensor models described in this chapter in that it does not return a sensed variable x consisting of its real value x plus a sensor error E, but instead generates a digital image simulating what a real camera would record based on the aircraft position and attitude as given by the actual or real trajectory x = x TRUTH .When provided with the camera pose with respect to the Earth at equally time spaced intervals, the simulation is capable of generating images that resemble the view of the Earth surface that the camera would record if located at that particular pose.To do so, it relies on three software libraries: • OpenSceneGraph [35] is an open source high performance 3D graphics toolkit written in C++ and OpenGL, used by application developers in fields such as visual simulation, games, virtual reality, scientific visualization and modeling.The library enables the representation of objects in a scene by means of a graph data structure, which allows grouping objects that share some properties to automatically manage rendering properties such as the level of detail necessary to faithfully draw the scene, but without considering the unnecessary detail that slows down the graphics hardware drawing the scene.
• osgEarth [36] is a dynamic and scalable 3D Earth surface rendering toolkit that relies on OpenSceneGraph, and is based on publicly available orthoimages of the area flown by the aircraft.Orthoimages consist of aerial or satellite imagery geometrically corrected such that the scale is uniform; they can be used to measure true distances as they are accurate representations of the Earth surface, having been adjusted for topographic relief, lens distortion, and camera tilt.When coupled with a terrain elevation model, osgEarth is capable of generating realistic images based on the camera position as well as its yaw and pitch, but does not accept the camera roll (in other words, the osgEarth images are always aligned with the horizon).
• Earth Viewer is a modification to osgEarth implemented by the author so it is also capable of accepting the bank angle of the camera with respect to the NED axes.Earth Viewer is capable of generating realistic Earth images as long as the camera height over the terrain is significantly higher than the vertical relief present in the image.As an example, figure 5 shows two different views of a volcano in which the dome of the mountain, having very steep slopes, is properly rendered.This section describes various calibration processes required for the determination of the fixed and run-to-run error contributions to the accelerometers, gyroscopes, magnetometers, and onboard camera.These procedures only need to be executed once and do not need to be repeated unless the sensors are replaced or their position inside the aircraft modified 23 .
The calibration procedures include the laboratory calibration of the accelerometers and gyroscopes described in section 6.1, the determination of the pose between the platform and body frames explained in section 6.2, the magnetometer calibration or swinging described in section 6.3, and the determination of the pose between the camera and body frames explained in section 6.4.Their main objective is the determination of the fixed contributions to the sensor error models 24 , this is, the scale factor and cross coupling errors of both inertial sensors and magnetometers ( MACC , MGYR , MMAG ) 25 , the magnetometers hard iron magnetism BHI,MAG , the body to platform transformation ( TB BP , φBP ), and the body to camera transformation ( TB BC , φBC ).These procedures also provide estimations for the run-to-run error contributions ( B0ACC , B0GYR , B0,MAG ), but these need to be discarded as they change every time the aircraft systems are switched on.

Inertial Sensors Calibration
Calibration is the process of comparing instrument outputs with known references to determine coefficients that force the outputs to agree with the references over a range of output values [9].The IMU inertial sensors need to be calibrated to eliminate the fixed errors originated from manufacturing, and also to determine their temperature sensitivity [19].The calibration process requires significant material and time resources, but it greatly reduces the measurement errors.While high grade IMUs are always factory calibrated, low cost ones generally are not, so it is necessary to calibrate the IMU at the laboratory before mounting it on the aircraft [5].
The calibration process is executed at a location where the position x GDT and the gravity26 vector g c have been previously determined with great precision [9].It relies on a three axis table, which enables rotating the IMU with known angular velocities into a set of predetermined precisely controlled orientations [19,37].Accelerometer and gyroscope measurements are then compared to reference values (gravity for the accelerometers, torquing rate plus Earth angular velocity for the gyroscopes), and the differences employed to generate corrections [9].
During calibration, the amount of time that the IMU is maintained stationary at each attitude, as well as the time required to rotate it between two positions, are trade offs based on two opposing influences.On one side, longer periods of time are preferred as the negative influence of the system noise in the measurements tends to even out over time, while on the other, shorter times imply smaller variations of the bias drift over the measurement interval.
It is worth noting that as the calibration is performed before the IMU is installed on the aircraft, it relies on the platform frame F P and the models contained in sections 3.6 and 3.7.Although it is possible to use a calibration strategy based on selecting platform orientations that isolate sensor input onto a single axis (for example, gravity will only be sensed by the accelerometer that is placed vertically with respect to the Earth surface) to then apply least squares techniques, in real life it is better to employ state estimation techniques27 to obtain estimates of the inertial sensors scale factors, cross coupling errors, and bias offsets [19,5].The process is repeated at different temperatures so the IMU processor can later apply the correction based on the IMU sensor temperature [5].

Estimation # Coefficients
This article assumes that the bias offset is exclusively a run-to-run source of error that varies every time the IMU is switched on, so the B0ACC and B0GYR coefficients obtained by calibration are discarded as they have no relation to the offsets that occur during flight.Modeling the calibration process results implies reducing the scale factor and cross couplings errors found on the inertial sensors specifications by an arbitrary amount of 95%.This reduction is already included in the baseline tables contained in section 3.10.To summarize, instead of applying (87, 88) to the measurements obtained by (61, 63), the simulation directly employs (61, 63) with reduced M GYR and M ACC taken from the tables of sections 3.10.

Determination of the Platform Frame Pose
The true relative pose between the body and platform frames (F B , F P ), given by T B BP and φ BP , as well as their estimated values TB BP and φBP , play a key role in the readings generated by the inertial sensors, as explained in section 3.9.Let's see how they can be estimated.
Considering that the position of the aircraft center of mass is known (in both full and empty tank configurations), the true displacement T B BP can be determined with near exactitude based on the IMU attachment point to the aircraft, resulting in the very small error σ TB BP assigned in table 7 for the estimation of TB BP in (52).With regards to the attitude φ BP , after mounting the IMU platform so two of its axes are approximately aligned with the forward and down directions of an approximate aircraft plane of symmetry (with no particular need for accuracy), measure the angular deviation φ BP by means of the self-alignment [5], resulting in the small error σ φBP assigned in table 7 for the estimation of φBP in (53).

Swinging or Magnetometer Calibration
Magnetometer calibration is inherently more complex than that of the inertial sensors as it must be performed with the sensors already mounted on the aircraft, as otherwise it would not capture the fixed contributions of the hard iron and soft iron magnetisms (section 4.1).The calibration process, known as swinging, relies on obtaining magnetometer readings while the aircraft is positioned at different attitudes that encompass a wide array of heading, pitch, and roll values [5], and is executed at a location where the magnetic field is precisely known.
The accuracy of the results is very dependent of the precision with which the different aircraft attitudes can be determined during swinging.This can be done with self-alignment procedures [5] or with the use of expensive static instruments.In any case, attitude accuracy is always going to be inferior to that obtained with a three axis table during inertial sensor calibration.Once the magnetic field readings are obtained, they are compared to the real magnetic field values, and expression (67) employed with least squares techniques to obtain estimations of the bias (sum of hard iron magnetism B HI,MAG and offset B 0,MAG ), and the scale factor and cross coupling matrix M MAG , which also includes the soft iron magnetism.The process can be repeated several times to isolate the influence of hard iron magnetism (a fixed effect that does not change) from the offset, which is a run-to-run error source that changes every time the magnetometers are turned on.

Figure 2 :Figure 2
Figure 2: Propagation with time of sensor error mean and standard deviation

Figure 3 :Figure 3
Figure 3: Propagation with time of first integral of sensor error mean and standard deviation Figure 3 follows the same pattern as figure 2 but applied to the error integral instead of to the error itself.It would represent the attitude error resulting from integrating the gyroscope output, or the velocity error expected when integrating the specific force measured by an accelerometer.The conclusions are the same as before but significantly more accentuated.Not only is the expected value of the error constant instead of zero (f BW (0) = 3

Figure 4 :
Figure 4: Propagation with time of second integral of sensor error mean and standard deviation

A
GNSS receiver enables the determination of the aircraft position and absolute velocity based on signals obtained from various constellations of satellites, such as GPS, GLONASS, and Galileo.The position is obtained by triangulation based on the accurate satellite position and time contained within each signal.Instead of derivating the position with respect to time, which introduces noise, GNSS receivers obtain the vehicle absolute velocity by measuring the Doppler shift between the constant satellite frequencies and those measured by the receiver.

Figure 5 :
Figure 5: Example of Earth Viewer images

Table 2 :
Components of sensed trajectory

Table 4 :
Typical inertial sensor biases according to IMU grade

Table 6 :
Units for single axis inertial sensor error sources 1/2 ] translates into an standard deviation for the error of 1 [deg] after 1 [sec], 10 [deg] after 100 [sec], and 1000 1/2 ≈ 31.6 [deg] after 1000 [sec].Manufacturers often provide this information as the power spectral density PSD of the white noise process in [deg 2 /hr 2 /hz] or equivalent, where it is necessary to take its square root to obtain σ v , or as the root PSD in [deg/hr/hz1/2] that is equivalent to σ v .Sometimes it is even provided as the PSD of the random walk process, not the white noise, in units [deg/hr] or equivalent.It is then necessary to multiply this number by the square root of the sampling interval ∆t or divide it by the square root of the sampling frequency to obtain the desired σ v value.

Table 7 :
ξ P N ξ P IMU mounting accuracy values Once the real pose between the F P and F B frames is established, it is necessary to specify its estimation employed T (52) φBP = φ BP + σ φBP N ψP , σ φBP N θP , σ φBP N ξP T (53) where the default standard deviations σ TB BP and σ φBP are shown in table 7, and N ψP , N θP , N ξP , N TP

Table 8
lists the error sources contained in the comprehensive inertial sensor error model represented by (61, 63).The first two columns list the different error sources, while the third column specifies their origin according to the criterion established in the first paragraph of section 3.1.The section where each error is described appears on the fourth column, followed by the seeds employed to ensure the results variability for different aircraft (Υ i,A ) as well as different flights (Υ j,F ):

Table 8 :
Inertial sensor error sourcesNote that B 0ACC , B 0GYR , σ uACC , σ uGYR , σ vACC , σ vGYR , s ACC , s GYR , m ACC , and m GYR are taken from the right hand column of tables 9 and 10 (section 3.10), while the remaining error sources come from section 3.8.It is worth pointing out that all errors are modeled as stochastic variables or processes (with the exception of the

Table 12 :
Magnetometer error sourcesNote that B HI,MAG , B 0,MAG , σ v,MAG , s MAG , and m MAG are taken from the right column of the specs listed above,

Table 13 :
Performances of GNSS receiver

Table 14
lists the error sources contained in the GNSS receiver model represented by (70, 71).The first two columns list the different error sources, while the third column specifies their origin according to the criterion established in the first paragraph of section 3.1.The last two columns list the seeds employed to ensure the results variability for different flights (Υ j,F ):

Table 14 :
GNSS receiver error sources 15ly include system noise, but the author has decided to also include a bias offset for added realism.The model contained in (77) hence relies on the bias offset B 0OAT and system noise σ OAT shown in table15, where OAT stands for outside air temperature, plus two uncorrelated standard normal random variables N (0, 1) (N 0,OAT and N s,OAT ): e OAT (s ∆t SENSED ) = e OAT (s ∆t) = T (s ∆t) − T (s ∆t) = B 0OAT N 0,OAT + σ OAT N s,OAT

Table 16 :
Performance of "Baseline" Pitot tube and air vanes

Table 18 :
Camera parameters 8, where each specific camera Euler angle is obtained as the product of the standard deviations (σ ψ C , σ θ C , σ ξ C ) contained in table 19 by a single realization of a standard normal random variable N (0, 1) (N ψ C , N θ C , and N ξ C ): BC + σ φBC N ψC , σ φBC N θC , σ φBC N ξC 19and σ φBC are shown in table19, and N ψC , N θC , N ξC , N TB

Table 19
lists the default values employed in the simulation, which can be adjusted by the user:

Table 19 :
Camera mounting accuracy values

Table 20 :
MACC 6 {ŝ ACC,i , mACC,ij } MGYR 9 {ŝ GYR,i , mGYR,ij } Results of calibration processThe twenty-one coefficients estimated in the calibration process are listed in table 20.Once the coefficients have been estimated, they can be introduced into the IMU processor so it automatically performs the corrections

Table 21 :
Results of swinging processThe fifteen coefficients estimated in the swinging process are listed in table 21.Once the coefficients have been estimated, they can be introduced into the processor so it automatically performs the corrections shown in (89):