Stochastic High Fidelity Simulation and Scenarios for Testing of Fixed Wing Autonomous GNSS-Denied Navigation Algorithms
Eduardo Gallo[The author holds a MSc in Aerospace Engineering by the Polytechnic University of Madrid and has twenty-two years of experience working in aircraft trajectory prediction, modeling, and flight simulation. He is currently a Senior Trajectory Prediction and Aircraft Performance Engineer at Boeing Research & Technology Europe (BR&TE), although he is publishing this article in his individual capacity and time as part of his PhD thesis titled “Autonomous Unmanned Air Vehicle GNSS-Denied Navigation”, advised by Dr. Antonio Barrientos within the Centre for Automation and Robotics of the Polytechnic University of Madrid.] [Contact<EMAIL_ADDRESS><https://orcid.org/0000-0002-7397-0425>]

§ ABSTRACT

Autonomous unmanned aerial vehicle (UAV) inertial navigation exhibits an extreme dependency on the availability of global navigation satellite systems (GNSS) signals, without which it incurs in a slow but unavoidable position drift that may ultimately lead to the loss of the platform if the GNSS signals are not restored or the aircraft does not reach a location from which it can be recovered by remote control. This article describes an stochastic high fidelity simulation of the flight of a fixed wing low SWaP (size, weight, and power) autonomous UAV in turbulent and varying weather intended to test and validate the GNSS-Denied performance of different navigation algorithms. Its open-source  implementation is available in [1].

Onboard sensors include accelerometers, gyroscopes, magnetometers, a Pitot tube, an air data system, a GNSS receiver, and a digital camera, so the simulation is valid for inertial, visual, and visual inertial navigation systems. Two scenarios involving the loss of GNSS signals are considered: the first represents the challenges involved in aborting the mission and heading towards a remote recovery location while experiencing varying weather, and the second models the continuation of the mission based on a series of closely spaced bearing changes.

All simulation modules have been modeled with as few simplifications as possible to increase the realism of the results. While the implementation of the aircraft performances and its control system is deterministic, that of all other modules, including the mission, sensors, weather, wind, turbulence, and initial estimations, is fully stochastic. This enables a robust evaluation of each proposed navigation system by means of Monte-Carlo simulations that rely on a high number of executions of both scenarios.

§ OUTLINE

This article begins with a short introduction to GNSS-Denied navigation, followed by an explanation of the approach taken to introduce randomness into the simulation by means of multiple stochastic parameters. The next section focuses on the high-fidelity modeling of the different simulation modules (with the exception of the navigation system), which is followed by a detailed description of both scenarios and the definition of the metrics employed for evaluation. Future articles will propose and evaluate different navigation systems based on the stochastic high fidelity simulation described in this article.

§ GNSS-DENIED NAVIGATION

UAVs can be classified into remotely controlled aircraft, in which the ground human operator continuously provides control commands to the platform, and autonomous aircraft, which rely on onboard computers to execute the previously uploaded mission objectives. In the case of remotely piloted platforms, the UAV also employs the communications channel to provide the ground with its current status (position, velocity, attitude, etc.) and in some cases photos or video, on which the ground operator relies to continuously generate the control instructions. Autonomous aircraft can also employ the communications channel to provide information to the ground, which may decide to update the mission objectives based on it, but as its own name implies, they usually operate without any kind of communication with its operator. They just continue executing its mission until the flight concludes or communications are restored.

The number, variety, and applications of UAVs have grown exponentially in the last few years, and the trend is expected to continue in the future. This is particularly true in the case of low SWaP vehicles because their reduced cost makes them suitable for a wide range of applications, both civil and military. A significant percentage of these vehicles are capable of operating autonomously. With small variations, these platforms rely on a suite of sensors that continuously provides noisy data about the airframe state, a navigation algorithm to estimate the aircraft pose (position plus attitude), and a control system that, based on the navigation output, adjusts the aircraft control mechanisms to successfully execute the preloaded mission.

This article focuses on fixed wing autonomous platforms, which are generally equipped with a GNSS receiver, accelerometers, gyroscopes, magnetometers, an air data system, a Pitot tube, and one or multiple cameras[A single camera is employed in the simulation.]. The combination of accelerometers and gyroscopes is known as the Inertial Measurement Unit (IMU). The errors introduced by all sensors grow significantly as their SWaP decreases, in particular in the case of the IMU. The recent introduction of solid state accelerometers and gyroscopes has dramatically improved the performance of low SWaP IMUs, with new models showing significant improvements when compared to those fabricated only a few years ago. The problem of noisy measurement readings is compounded in the case of low SWaP vehicles by the low mass of the platforms, which results in less inertia and hence more high frequency accelerations and rotations caused by the atmospheric turbulence.

Aircraft navigation has traditionally relied on the measurements provided by accelerometers, gyroscopes, and magnetometers, incurring in an slow but unbounded position drift that could only be stopped by triangulation with the use of external navigation (radio) aids. More recently, the introduction of satellite navigation has completely removed the position drift and enabled autonomous inertial navigation in low SWaP platforms [2, 3, 4]. On the negative side, low SWaP inertial navigation exhibits an extreme dependency on the availability of GNSS signals. If not present, inertial systems rely on dead reckoning, which results in velocity and position drift, with the aircraft slowly but steadily deviating from its intended route.

The availability of GNSS signals cannot be guaranteed by any means. In addition to the (unlikely) event of one or various GNSS satellites ceasing to broadcast (voluntarily or not), the GNSS signals can be accidentally blocked by nearby objects (mountains, buildings, trees), corrupted by their own reflections in those same objects, or maliciously interfered with by broadcasting a different signal in the same frequency. Any of the above results in what is known as GNSS-Denied navigation. In that event, the vehicle is unable to fly its intended route or even return to a safe recovery location, which leads to the uncontrolled loss of the airframe if the GNSS signal is not recovered before the aircraft runs out of fuel (or battery in case of electric vehicles). [5] provides a summary of the challenges of GNSS-Denied navigation and the research efforts intended to improve its performance.

The extreme dependency on GNSS availability is not only one of the main impediments for the introduction of small UAVs in civil airspace, where it is not acceptable to have uncontrolled vehicles causing personal or material damage, but it also presents a significant drawback for military applications, as a single hull loss may compromise the onboard technology. At this time there are no comprehensive solutions to the operation of low SWaP autonomous UAVs in GNSS-Denied scenarios, although the use of onboard cameras to provide an additional relative pose measurement seems to be one of the most promising routes. Bigger and more expensive UAVs, this is, with less stringent SWaP requirements, can rely to some degree on much more accurate IMUs (at the expense of SWaP), and additional communications equipment to overcome this problem, but for most autonomous UAVs, the permanent loss of the GNSS signal is equivalent to losing the airframe in an uncontrolled way.

Visual navigation relies on the images taken by one or more onboard cameras to estimate the variation with time of the aircraft position and attitude [6, 7]. Based  on visual Odometry (VO), it also incurs in a slow but unbounded position drift that can only be eliminated if the vehicle revisits a given zone for a second time, a technique known as simultaneous localization and mapping (SLAM). If this is not the case, the drift can also be removed through image registration, in which the images are compared with those stored in a database. Inertial and visual navigation techniques can be integrated to improve their combined GNSS-Denied navigation capabilities, resulting in a visual inertial navigation system or visual inertial odometry (VIO).

§ STOCHASTIC SIMULATION

The trajectory flown by an aircraft is the result of a combination of physical and electronic processes, some deterministic and others stochastic[The aircraft performances and its control system are both deterministic, as given inputs will always results in the same outputs. On the other hand, the atmospheric turbulence and the errors introduced by the onboard sensors are stochastic and continuously varying.]. This duality also occurs in the high fidelity simulation, in which the deterministic or stochastic nature of each process is present in the simulation module that models it, hence increasing the realism of the results. As such, a single simulated flight is not sufficient to evaluate the performance of a given GNSS-Denied navigation system, as the randomness implicit to all stochastic processes may result in conditions that are favorable or unfavorable for the algorithm being evaluated, leading to results that are too optimistic or pessimistic.

The most rigorous way of testing any proposed navigation algorithm is to evaluate it by means of a Monte Carlo simulation, in which a sufficiently high number of executions or runs () of a certain scenario are performed, and the results are aggregated to obtain meaningful and realistic metrics that describe the algorithm behavior. It is also necessary to ensure that the stochastic parameters, in addition to randomly changing between different Monte Carlo executions, are also repeatable if a given run needs to be executed a second time.

Given the high number of stochastic parameters influencing each flight, it is necessary to employ a high number  of executions or runs in each Monte Carlo simulation, with 100 considered as the default. The process begins with the initialization of a discrete uniform distribution with a given seed (any value is valid, so 1 was employed by the author), which generates pseudo random integers in which each possible value has an equal likelihood of being produced. This distribution is realized  times to generate the trajectory seeds ,   j = 1 …. The stored trajectory seeds  become the initialization seeds for each of the flight simulation executions, so this step does not need to be repeated.

Every time the Monte Carlo simulation needs to generate a trajectory, it is initialized with the corresponding trajectory seed . As this seed is the only input required for all the stochastic processes within the simulation, the results of a given run can always be repeated by employing the same seed. At the beginning of each execution, the trajectory seed  is employed to initialize a second discrete uniform distribution, which is realized seventeen times (always in the same order to maintain repeatability) to generate the seeds employed by the different simulation modules. These seeds are listed in table <ref>, and their purpose is explained when applicable later in this article.

Seed     Purpose

, ,   		    Accelerometers, gyroscopes, and magnetometers

, ,  		    Barometer, thermometer, and GNSS receiver

, ,   		    Pitot tube and air vanes

IMU position and attitude

Digital camera

, ,      Wind, weather, and turbulence

Mission guidance objectives

Gravity and magnetic field realism

Fine alignment results

tableSeeds required to initialize the different simulation stochastic processes

§ HIGH FIDELITY SIMULATION

This article presents a stochastic high fidelity simulation of the flight in turbulent conditions and varying weather of an autonomous fixed wing low SWaP vehicle, with the purpose of testing the capabilities of different inertial, visual, and visual inertial navigation algorithms in GNSS-Denied conditions. The simulation open source  code is available in [1]. It is intended to serve as a test bench for the development and evaluation of different GNSS-Denied navigation techniques. For that purpose, the various components of the simulation have been modeled with as few simplifications as possible to increase the realism of the results under two different scenarios, one representing the challenges involved in aborting the mission and heading towards a remote recovery location while experiencing varying weather, and the second modeling the continuation of the mission based on a series of closely spaced bearing changes.

[auto, node distance=2cm,>=latex']
[coordinate](x0input) ;
[coordinate, below of=x0input, node distance=3.0cm] (deltaCNTRinput);
[blockyellow, right of=x0input, minimum width=1.8cm, node distance=3.0cm, align=center, minimum height=1.25cm] (FLIGHT) ;
[coordinate, below of=FLIGHT, node distance=1.75cm] (midblocks);
[blockyellow, right of=midblocks, minimum width=1.8cm, node distance=1.25 cm, align=center, minimum height=1.25cm] (EARTH) ;
[blockyellow, left of=midblocks, minimum width=1.8cm, node distance=1.25 cm, align=center, minimum height=1.25cm] (AIRCRAFT) ;
[coordinate, right of=FLIGHT, node distance=2.8cm] (crosspoint1);
[coordinate, below of=crosspoint1, node distance=2.7cm] (crosspoint2);
[coordinate, left of=crosspoint2, node distance=5.5cm] (crosspoint3);
[blockgreen, right of=FLIGHT, minimum width=2.6cm, node distance=5.0cm, align=center, minimum height=1.25cm] (SENSORS) ;
[blockgreen, right of=SENSORS, minimum width=2.6cm, node distance=5.0cm, align=center, minimum height=1.25cm] (NAVIGATION) ;
[coordinate, right of=NAVIGATION, node distance=2.6cm] (x0estinput);
[coordinate, below of=NAVIGATION, node distance=1.5cm](tmp1) ;
[blockgreen, below of=SENSORS, minimum width=2.6cm, node distance=3.0cm, align=center, minimum height=1.25cm] (CONTROL) ;
[blockgreen, below of=NAVIGATION, minimum width=2.6cm, node distance=3.0cm, align=center, minimum height=1.25cm] (GUIDANCE) ;
[coordinate, right of=GUIDANCE, node distance=2.6cm](xrefinput) ;

[->] ($(EARTH.north)-(0.50cm,0cm)$) – ($(FLIGHT.south)+(0.75cm,0cm)$);
[->] (EARTH.west) – (AIRCRAFT.east);
[->] ($(AIRCRAFT.north)+(0.50cm,0cm)$) – ($(FLIGHT.south)-(0.75cm,0cm)$);
[->] (x0input) – node[pos=0.4]  (FLIGHT.west);
[->] (x0estinput) – node[pos=0.4]  (NAVIGATION.east);
[->] (FLIGHT.east) – node[pos=0.5] =  (SENSORS.west);
[->] (SENSORS.east) – node[pos=0.5] =  (NAVIGATION.west);
[black] (tmp1) circle [radius=1pt];
[->] (NAVIGATION.south) – node[pos=0.95] =  (tmp1) -| (CONTROL.north);
[->] (tmp1) – (GUIDANCE.north);
[->] (xrefinput) – node[pos=0.4]  (GUIDANCE.east);
[->] (GUIDANCE.west) – node[pos=0.5]  (CONTROL.east);
[->] (CONTROL.west) -| node[pos=0.25]  (deltaCNTRinput) |- ($(AIRCRAFT.west)+(0cm,0.4cm)$);
[black] (crosspoint1) circle [radius=1pt];
[->] (crosspoint1) |- (EARTH.east);
[->] (crosspoint1) – (crosspoint2) – (crosspoint3) |- ($(AIRCRAFT.west)-(0cm,0.4cm)$);

Components of the high fidelity simulation

It is necessary to remark that the focus is on low SWaP autonomous UAVs, ruling out the use of high quality sensors, which in general are bigger and have more weight, heavier platforms with more inertia against atmospheric turbulence, as well as the assistance of any kind of communications between the platform and the ground. Platforms that generate their lift by means of rotating blades (helicopters and multirotors) are also excluded from the simulation. Vertical Take Off and Landing (VTOL) as well as Short Take Off and Landing (STOL) vehicles, which use their rotors to generate lift when the airspeed is low, but otherwise are capable of flying like conventional fixed wing aircraft, can be included but only when behaving like fixed wing platforms. Although the modeled platform employs a piston engine to power the propeller, there would be no significant difference if it were replaced by an electric motor.

The objective is to develop and implement a sufficiently realistic simulation so as to enable the implementation and testing of new navigation algorithms to verify if they result in improvements when operating in GNSS-Denied conditions, and if they do, to accurately evaluate the benefits. The simulation models the influence on the resulting aircraft trajectory of many different factors, such as the atmospheric conditions, the wind field, the air turbulence, the aircraft aerodynamic and propulsive performances, its onboard sensors and their error sources, the guidance objectives that make up the mission, the control system that moves the throttle and the aerodynamic controls so the trajectory conforms to the guidance objectives, and the navigation system that processes the data obtained by the sensors and feeds the control system. Figure <ref> shows the different components of the simulation.

The simulation consists on two distinct processes. The first, represented by the yellow blocks on the left of figure <ref>, focuses on the physics of flight and the interaction between the aircraft and its surroundings that results in the actual or real aircraft trajectory = ; the second, represented by the green blocks on the right, contains the aircraft systems in charge of ensuring that the resulting trajectory adheres as much as possible to the mission objectives. It includes the different sensors whose output comprise the sensed trajectory = , the navigation system in charge of filtering it to obtain the estimated trajectory = , the guidance system that converts the reference objectives  into the control targets , and the control system that adjusts the position of the throttle and aerodynamic control surfaces  so the estimated trajectory  is as close as possible to . As shown in the figure, the two parts of the simulation are not independent. The total error or difference between  and  is the combination of the navigation system error (difference between  and ) and the flight technical error (difference between  and ). The different simulation components are described below:

§.§.§ Earth

The Earth model conforms the underlying frame in which the aircraft motion and its interaction with the environment takes place:

* The World Geodetic System 1984 (WGS84) [8] ellipsoid is employed for the Earth surface.

* The Somigliana ellipsoidal formula provided by [8] is adopted for the gravity acceleration  as an approximation to the combination of the Earth Gravitational Model 1996 (EGM96) [9] gravitation and the Earth centrifugal acceleration.

* An spherical Earth of radius  [10] is employed to obtain the ratio between geodetic and geopotential altitudes (h and H) [11].

* The wind is modeled as the combination of its low and high frequency components (wind field and turbulence, respectively). The wind field provides the horizontal wind (intensity and direction) as a function of time and geodetic coordinates  (longitude λ, latitude φ, and altitude h) and considers zero vertical wind. The Dryden model [12] is employed for the turbulence or high frequency wind.

* A generalization of the International Civil Aviation Organization (ICAO) Standard Atmosphere (ISA) [11] is used to model the atmosphere (relationships between pressure p, temperature T, pressure altitude , and geopotential altitude H). Called ICAO Non Standard Atmosphere or INSA [13], it relies on the same hypotheses as ISA but extends it with two parameters, the temperature and pressure offsets ( and ), which are zero for ISA. In this way, the temperature in standard mean sea level conditions (= 0) is +  instead of  for ISA, while the pressure in mean sea level conditions (H = 0) is +  instead of  for ISA.

* The weather model provides the temperature and pressure offsets as functions of time and horizontal position (longitude and latitude).

* The World Magnetic Model (WMM) is employed for the Earth magnetic field [14].

* In order to represent the EGM96 and WMM model errors as well as local effects that result in the onboard models never being fully equal to the real magnetism and gravity experienced by the aircraft, the magnetic and gravity models employed when simulating the aircraft motion are slightly different than those employed by the navigation filter, with the differences stochastically taken from [14, 15].

The ellipsoid, gravity field, INSA atmospheric relationships, and Earth magnetism are all deterministic and do not vary. The atmospheric turbulence encountered by the aircraft during its flight and the local variations between the real gravity and magnetic fields experienced by the aircraft with the models employed by the navigation system are stochastic, being randomly generated for each execution  based on their respective seeds  and . Finally, although the weather and wind fields consist of user provided inputs, the two scenarios described below define them according to stochastic parameters randomly generated based on the  and  seeds.

§.§.§ Aircraft

The aircraft model provides the aerodynamic and propulsive forces and moments generated by the platform, together with the fuel consumption, the center of gravity position, and the inertia matrix. The modeled platform is a high wing aircraft (maximum mass 19.715 kg, length 2.32 m, wing span 2.68 m) equipped with a ∧ tail, a two stroke atmospheric gasoline engine with a maximum power of 4.18 kW, and a two blade fixed pitch pusher propeller with a diameter of 0.51 m.

Aircraft overall view

* Mass and inertia analysis is employed to provide the quasi stationary relationship between the inertia matrix and the center of gravity position with the fuel load.

* Athena Vortex Lattice [16] is used to obtain the dependencies of the platform aerodynamic force and moment coefficients with the angles of attack α and sideslip β, the position of the control surfaces  (elevator, ailerons, and rudder), and the platform angular velocity components [17]. The results have been verified for longitudinal and lateral stability [18]. Incompressible air flow is employed to compute the aerodynamic forces and moments.

* The propulsive force (thrust) and moment (propeller torque) are obtained by combining a model that provides the engine power and fuel consumption as functions of the throttle position and atmospheric pressure and temperature [19], with a propeller model that provides its thrust and power coefficients based on the advance ratio [20].

The aircraft performances are deterministic and do not vary for the different executions of the Monte Carlo simulation.

§.§.§ Flight

Given the initial conditions  and the variation with time of the throttle and control surfaces t, the actual or real trajectory t is the outcome of the interaction between the aircraft and the environment that surrounds it. Considering as inertial a reference frame centered at the Sun with axes fixed with respect to the stars, it is possible to establish a system of first order differential equations composed by the kinematic, force, rotation, moment, and mass equations, which is based on a state vector comprised by the aircraft geodetic coordinates , absolute velocity viewed in the body frame , attitude from NED to body represented by its unit quaternion , inertial angular velocity viewed in body , and mass m [18]. Note that the Earth and aircraft performance algebraic equations are indirectly included in the system of differential equations, as shown in figure <ref>.

\begin{eqnarray}
\nm{\xvec\lrp{t}} & = & \nm{\lrsb{\xEgdt\lrp{t}, \, \vB\lrp{t}, \, \qNB\lrp{t}, \, \wIBB\lrp{t}, \, m\lrp{t}}^T}\label{eq:flight_state_vector} \\
\nm{\xvecdot\lrp{t}} & = & \nm{\vec f \big(\xvec\lrp{t}, \, t, \, \deltaCNTR\lrp{t}\big)}\label{eq:flight_differentials} \\
\nm{\xvec\lrp{0}} & = & \nm{\xveczero}\label{eq:flihgt_initial}
\end{eqnarray}

A fourth order Runge-Kutta scheme at a frequency of 500 hz is employed to integrate it and obtain the actual trajectory t_t = t_t = t ·. A standard Runge-Kutta integration scheme implies considering the unit quaternions as 4-vectors that belong to ℝ^4 instead of 𝕊𝕆(3) and operating with the addition and scalar multiplication of the ℝ^4 vector space; the propagated 4-vectors hence do not belong to 𝕊𝕆3, this is, are not unit quaternions, and must be normalized after each step to eliminate the deviation from the three dimensional manifold 𝕊𝕆3 incurred by the ℝ^4 propagation.

Although at 500 hz the errors are relatively small, it is more rigorous and precise to integrate the unit quaternion respecting the constraints of its Lie group, this is, ensuring that  propagates without deviating from the 𝕊𝕆(3) manifold. To do so, the fourth order Runge-Kutta scheme can be modified so its relies on the 𝕊𝕆(3) ⊕ operator and the relationship between the  unit quaternion and the angular velocity  belonging to its local tangent space 𝔰𝔬(3) [21]. Note that  can easily be obtained from  and the rest of the state vector components.

There are no stochastic parameters influencing the integration of the equations of motion, which is a fully deterministic process.

§.§.§ Sensors

The onboard sensors measure various aspects of the actual trajectory and provide these measurements to the navigation system. It is assumed that all sensors are attached to the aircraft structure in a strapdown configuration, provide instantaneous measurements, and are time synchronized with each other. With the exception of the GNSS receiver that operates at 1 hz and the digital camera that works at 10 hz, all other sensors obtain measurements at a frequency of 100 hz. The sensed trajectory =  is composed by the specific force  measured by the accelerometers, the inertial angular velocity  provided by the gyroscopes, the magnetic field  measured by the magnetometers, the geodetic coordinates _GDT and ground velocity  provided by the GNSS receiver (when operative), the atmospheric pressure p and temperature T provided by the barometer and thermometer, the true airspeed  measured by the Pitot tube, the angles of attack α and sideslip β provided by the air vanes, and the digital image of the ground I⃗ obtained by the camera:
t_s = s · = ,   ,   ,   _GDT,   ,   p,   T,   ,   α,   β,   I⃗^Teq:sensors

This section provides a summary of the modeling of the error sources of the various sensors, which are described in detail in [22]. Inertial sensors (accelerometers and gyroscopes) are modeled taking into consideration fixed error sources (scale factor, cross coupling, lever arm, IMU body attitude), run-to-run bias offset, and in-run bias drift and system noise [2, 3, 4]. The bias and system noise models are based on [23, 24]. Inertial sensors are assumed to be lab calibrated, which reduces their scale factor and cross coupling errors.

The simulation considers that the IMU lever arm or position with respect to the aircraft structure is fixed or deterministic. However, its orientation with respect to the body frame is stochastic and hence varies randomly for each execution  based on . The IMU knowledge about its own position and attitude with respect to body, which are necessary inputs to obtain  and , is not perfect, and the errors are also stochastically controlled by . Standard normal distributions applied to the inertial sensors specifications are initialized with the  and  seeds and employed to obtain the value of all fixed, run-to-run, and in-run error contributions.

The magnetometers model is simpler than that of the inertial sensors as it does not include bias drift, although in exchange it needs to consider the negative effects of hard and soft iron magnetism. Additionally, their position and attitude with respect to the aircraft structure have no influence in the readings as long as they are mounted away from permanent magnets and electrical equipment. The scale factor, cross coupling, and hard iron magnetism errors can be reduced if the magnetometers are calibrated through swinging once installed on the aircraft. All error sources are stochastically generated by applying standard normal distributions initialized with  to the magnetometer specifications.

The barometer, thermometer, Pitot tube, and air vanes (measuring the angles of attack and sideslip) are modeled by a combination of system noise and bias offset exclusively. All parameters vary stochastically in each execution of the Monte Carlo simulation, and are controlled by their respective seeds (,   ,   ,   ,   ).

Before the GNSS signals are lost, the navigation system can rely on the position and velocity readings provided by a GNSS receiver. The position error is modeled as the sum of a white noise process plus slow varying ionospheric effects [25], which in turn are modeled as the sum of the bias offset plus a random walk. In the case of the ground velocity error, only system noise is considered. As in previous cases, the simulation errors are obtained by applying standard normal distributions initialized with the  seed to the receiver specifications.

The simulation also includes a single onboard camera working at a frequency of 10 [hz] and rigidly attached to the aircraft structure. The simulated image generation process is instantaneous and relies on the shutter speed being sufficiently high that all images are equally sharp. In addition, the camera  setting remains constant during the flight, all generated images are noise free, the visible spectrum radiation reaching all patches of the Earth surface remains constant, and in addition the terrain is considered Lambertian, so its appearance at any given time does not vary with the viewing direction. These assumptions imply that a given terrain object is represented with the same luminosity in all images, even as its relative pose (position and attitude) with respect to the camera varies [26]. Geometrically, the simulation adopts a perspective projection or pinhole camera model, which in addition is perfectly calibrated and hence shows no distortion [26].

As in the case of the inertial sensors, the simulation considers that the location of the camera onboard the aircraft is deterministic but its orientation is stochastic, as is the navigation system estimation of both the camera position and attitude. All these parameters vary for each execution  based on the  seed.

The digital camera differs from all other sensors in that it does not return a sensed variable  consisting of its real value  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 trajectory = . When provided with the camera pose 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 position and attitude. To do so, it relies on three software libraries:

*  [27] is an open source high performance 3D graphics toolkit written in  and , 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.

*  [28] is a dynamic and scalable 3D Earth surface rendering toolkit that relies on , 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,  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[This is analogous to stating that  images are always aligned with the horizon.].

*  is a modification to  implemented by the author so it is also capable of accepting the bank angle of the camera with respect to the  axes.  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 <ref> shows two different views of a volcano in which the dome of the mountain, having very steep slopes, is properly rendered.

Example of  images

§.§.§ Guidance and Control

The guidance and control systems work together so the estimated trajectory =  resembles the reference objectives  as much as possible. In the simulation they both operate at 50 hz. The guidance system supplies the control system with the proper targets , and this in turn adjusts the position of the throttle and aerodynamic control surfaces. As the aircraft has four controls (throttle, elevator, ailerons, and rudder), the guidance system output is at all times composed by four targets and a trigger . The trigger is a boolean expression that is continuously evaluated based on ; when its sign switches, the current control targets are replaced by the next targets in the pipeline.

\begin{eqnarray}
\nm{\xREF} & = & \nm{\lrb{\vec{\delta}_{\sss{TARGET,1}}, \, \vec{\delta}_{\sss{TARGET,2}}, \, \ldots , \vec{\delta}_{\sss{TARGET,n}}}}\label{eq:guidance1} \\
\nm{\deltaTARGET\lrp{t_c}} & = & \nm{\deltaTARGET\lrp{c \cdot \DeltatCNTR} = \lrsb{\deltaTARGETT, \, \deltaTARGETE, \,  \deltaTARGETA, \, \deltaTARGETR, \, \deltaTRG}^T}\label{eq:guidance2} \\
\nm{\deltaCNTR\lrp{t_c}} & = & \nm{\deltaCNTR\lrp{c \cdot \DeltatCNTR} = \lrsb{\deltaT, \ \deltaE, \ \deltaA, \ \deltaR}^T}\label{eq:control}
\end{eqnarray}

The control system is implemented as four primary proportional integral derivative (PID) loops, each associated to a certain variable and control mechanism (true airspeed  for the throttle , body pitch angle θ for the elevator , body bank angle ξ for the ailerons , and sideslip angle β for the rudder ). In addition, some secondary PID loops are also employed (pressure altitude  and aerodynamic path angle  for the elevator, bearing χ for the ailerons), in which the output of the secondary loop acts as the input for the primary one. All loops contain a low pass filter to smooth the derivative inputs [29] and use set point ramping to reduce overshooting.

The guidance and control algorithms have been adjusted to the aircraft performances and are hence deterministic. The reference objectives  are however user provided inputs that in the two scenarios described below are randomly generated based on the  seed.

§.§.§ Navigation

The mission of the navigation system is to obtain the estimated trajectory =  based on the sensor outputs =  and an initial estimation , as shown in figure <ref>. As the main objective of this simulation is to act as a bench test to analyze the behavior of different inertial and visual navigation systems in GNSS-Denied conditions, no navigation system is included by default as part of the simulation. Follow-on articles will include navigation algorithms proposed by the author as well as the results obtained when applied to Monte Carlo simulations of the two scenarios defined below.

Both inertial and visual navigation systems require an initial estimation of their respective state vectors . While in some occasions these initial estimations can be directly taken from the sensor measurements, other observed trajectory variables are not measured by any sensors, and their initial estimation relies on the pre-flight procedures, and in particular in the outputs of the fine alignment process [2, 3, 4].

When the initial estimation can be directly taken from a given sensor, as is the case with the initial position (provided by the GNSS receiver) required by both inertial and visual systems, no additional randomness is required by the simulation. However, the  seed is employed together with a standard normal distribution to estimate the results of the fine alignment process, which  results in initial estimations for the aircraft attitude (required by both inertial and visual systems) as well as the initial errors of the accelerometers, gyroscopes, and magnetometers.

§ SCENARIO #1

Scenario #1 has been selected with the objective of adequately representing the challenges faced by an autonomous fixed-wind UAV that suddenly cannot rely on the GNSS signals and hence changes course to reach a predefined recovery location situated at approximately one hour of flight time. In the process, in addition to executing an altitude and airspeed adjustment, the autonomous aircraft faces significant weather and wind field changes that make its GNSS-Denied navigation even more challenging. The following sections describe the different parameters that together define the scenario, with special emphasis on the stochastic ones whose values vary in the different executions or runs of a given Monte-Carlo simulation.

§.§.§ Mission

The mission, guidance targets, or reference objectives of scenario #1 consists in the six flight segments graphically shown in figure <ref>. The scenario is intended to showcase the ability of an autonomous aircraft to reach a predefined recovery location once it is flying in GNSS-Denied conditions:

* The first segment starts with the aircraft flying at constant airspeed , constant pressure altitude , constant bearing , and no sideslip β= 0. As the navigation filter is employing the GNSS signals, the solution is stable and hence the difference between the observed trajectory =  and the actual trajectory =  does not increase with time. For this reason, there is no lack in generality in the assumption that the GNSS signal gets lost at = 100   sec[The results would be similar if the aircraft flew with GNSS for a much longer period of time.]. It is also realistic to assume that the aircraft will continue flying in the same way for a short period of time until confirming that it can not recover the GNSS signal, and then change its mission at  by performing a bearing change in an attempt to reach a predefined recovery point.

* The turn is executed at the same airspeed  and pressure altitude  with a body bank angle of = ±10   deg, with the sign given by the shortest turn necessary to reach the final bearing , upon which the turn concludes.

(center) at (+0.0,+0.0);
(pointC) at ($(center) + (210:0.5cm)$);
(pointD) at ($(center) + (110:0.5cm)$);
(pointB) at ($(pointC) + (-60:1.0cm)$);
(pointA) at ($(pointC) + (-60:3.0cm)$);
(pointE) at ($(pointD) + (20:3.0cm)$);
(pointF) at ($(pointD) + (20:5.0cm)$);
(pointG) at ($(pointD) + (20:6.0cm)$);
(pointH) at ($(pointD) + (20:10.0cm)$);

[blue, name path = segmentAB, ultra thick] [->] (pointA) – (pointB);
[blue, name path = segmentBC, ultra thick] [->] (pointB) – (pointC);
(pointC) [blue, name path = segmentCD, ultra thick] arc[start angle=210, end angle=110, radius=0.5cm];
[blue, name path = segmentDE, ultra thick] [->] (pointD) – (pointE);
[blue, name path = segmentEF, ultra thick] [->] (pointE) – (pointF);
[blue, name path = segmentFG, ultra thick] [->] (pointF) – (pointG);
[blue, name path = segmentGH, ultra thick] [->] (pointG) – (pointH);

(pointAref) at ($(pointA) + (90:1.2cm)$);
[name path = Aref, thin] [-] (pointA) – (pointAref);
($(pointA) + (90:1.0cm)$) [name path = pointAangle, thin] arc[start angle=90, delta angle=30, radius=1.0cm] [->] node [pos=0.5, above=6pt] - ;

(pointDref) at ($(pointD) + (90:1.2cm)$);
[name path = Dref, thin] [-] (pointD) – (pointDref);
($(pointD) + (90:1.0cm)$) [name path = pointDangle, thin] arc[start angle=90, delta angle=-70, radius=1.0cm] [->] node [pos=0.7, above=10pt] ;

(pointAA) at ($(pointA) + (+20:1.0cm)$);
[black!50, thin] [->] (pointA) – (pointAA) node [black, right=1pt, text width=4.5cm] Start: targets are , , and ;

(pointBB) at ($(pointB) + (+260:1.0cm)$);
[black!50, thin] [->] (pointB) – (pointBB) node [black, left=1pt] GNSS signal lost at ;

(pointCC) at ($(pointC) + (-15:2.0cm)$);
[black!50, thin] [->] (pointC) – (pointCC) node [black, right=1pt, text width=5cm] Initiates turn at  with bank angle ;

(pointDD) at ($(pointD) + (+0:2.0cm)$);
[black!50, thin] [->] (pointD) – (pointDD) node [black, right=1pt, text width=4.5cm] Concludes turn when reaching new bearing target ;

(pointEE) at ($(pointE) + (+130:1.5cm)$);
[black!50, thin] [->] (pointE) – (pointEE) node [black, left=0pt, text width=3.8cm] After , changes airspeed to ;

(pointFF) at ($(pointF) + (+330:1.0cm)$);
[black!50, thin] [->] (pointF) – (pointFF) node [black, right=1pt, text width=4.5cm] After , initiates change of altitude with path angle ;

(pointGG) at ($(pointG) + (+145:2.0cm)$);
[black!50, thin] [->] (pointG) – (pointGG) node [black, left=0pt, text width=3.8cm] Maintains altitude after reaching ;

(pointHH) at ($(pointH) + (+180:1.0cm)$);
[black!50, thin] [->] (pointH) – (pointHH) node [black, left=1pt] Scenario concludes at ;

Schematic view of the scenario #1 mission

Note that the real objective of the aircraft is to reach a predefined recovery location, not to fly a given bearing . However, one leads to the other, as the aircraft needs to be able to follow a bearing in order to reach a location; and if it is capable of reaching a location, it means that it can also follow a bearing. From the point of view of setting up the scenario and evaluating the results, it is more convenient for each seed to imply the use of a given bearing  for the same length of time (see below) than to convert each seed to a different longitude and latitude which the aircraft would take different times to reach due to wind variations.

* Once it captures the bearing that should take it to its intended destination, the aircraft follows it at the same airspeed  and pressure altitude  for a given period of time .

* At this time the aircraft changes its airspeed to  without modifying its altitude  or course . Note that a change of airspeed is not necessary to reach the intended destination but is included to show that the proposed systems are capable of accelerating or decelerating the aircraft without the use of GNSS signals. These conditions are maintained for a given period of time .

* The aircraft then modifies its pressure altitude without changing its airspeed  or bearing  to avoid the terrain or to get an improved exposure to radio control signals with which it may be recovered from its destination. As in the previous segment, note that this change in pressure altitude is in theory not necessary but is included to prove that it can be achieved in GNSS-Denied conditions. The maneuver is executed with an aerodynamic path angle of = ±2  deg, with the sign given by whether the final pressure altitude  is higher or lower than the existing one.

* The last segment involves the aircraft flying at the existing airspeed, pressure altitude, and bearing until time = 3800   sec. The use of such a long flight time (in the order of one hour) in GNSS-Denied conditions ensures that no instabilities or other negative effects in the navigation system are hidden by not allowing the aircraft to fly long enough in GNSS-Denied conditions.

The above description of the guidance targets or reference trajectory can also be expressed in a more formal way according to the previously introduced notation:

\begin{eqnarray}
\nm{\xREF}           & = & \nm{\lrb{\deltaTARGETone, \, \deltaTARGETtwo, \, \deltaTARGETthr, \, \deltaTARGETfou, \, \deltaTARGETfiv, \, \deltaTARGETsix}} \label{eq:Sim_Scenario_Mission_RT} \\
\nm{\deltaTARGETone} & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \chiINI, \, \beta = 0, \, \tTURN}^T}\label{eq:Sim_Scenario_Mission_target1} \\
\nm{\deltaTARGETtwo} & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \xiTURN, \, \beta = 0, \, \chiEND}^T}\label{eq:Sim_Scenario_Mission_target2} \\
\nm{\deltaTARGETthr} & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \chiEND, \, \beta = 0, \, \DeltatOPTAS}^T}\label{eq:Sim_Scenario_Mission_target3} \\
\nm{\deltaTARGETfou} & = & \nm{\lrsb{\vtasEND, \, \HpINI, \,  \chiEND, \, \beta = 0, \, \DeltatOPHp}^T}\label{eq:Sim_Scenario_Mission_target4} \\
\nm{\deltaTARGETfiv} & = & \nm{\lrsb{\vtasEND, \, \gammaTASCLIMB, \,  \chiEND, \, \beta = 0, \, \HpEND}^T}\label{eq:Sim_Scenario_Mission_target5} \\
\nm{\deltaTARGETsix} & = & \nm{\lrsb{\vtasEND, \, \HpEND, \,  \chiEND, \, \beta = 0, \, \tEND}^T}\label{eq:Sim_Scenario_Mission_target6}
\end{eqnarray}

Component     Variable     Definition     Unit     Restrictions

Initial Airspeed                         N29, 1.5^2           [m/sec]     24 < < 34

Initial pressure altitude                  N2700, 200^2         [m]

Initial bearing                           U-179, 180           [deg]

Sideslip angle                β              0                          [deg]

GNSS-Denied time		      		        100		   		          [sec]

Turn time				      		        + N(30, 50^2)       [sec]       > + 15

Turn bank angle			      	        ±10			          [deg]

Final bearing                             U-179, 180           [deg]       |- |> 10

First time interval                  N500, 100^2          [sec]       > 150

2*Final airspeed     2*     2*N, 1.5^2     2*[m/sec]     24 < < 34

|- |> 0.5

2*Second time interval     2*     2*N500, 100^2     2*[sec]     > 150

+ < 2500

Path angle				           ±2                      [deg]

Final pressure altitude                    N, 300^2       [m]         |- |> 100

Final time                                  3800					      [sec]

tableDeterministic and stochastic components of scenario #1 mission

Table <ref> shows the definition of the deterministic and stochastic parameters that result in the scenario #1 mission. Note that it has been developed with the objective of including as much variability as possible, as the only deterministic variables are those that do not result in loss of generality. For each Monte Carlo simulation run , the  seed is employed to initialize a standard normal distribution N and a discrete uniform distribution U, which are then realized several times each to obtain the different parameters listed in table <ref>. This ensures that the results are repeatable if the same seed is employed again.

colormap name = bluered,
cycle list=[of colormap],
x dir=reverse,
axis lines*=none,
axis line style = draw=transparent,line width=0.0001pt,
xticklabel style=anchor=--90,
axis line style=-stealth,
legend entries=  deg,   deg,
legend style=font=,
legend cell align=right,
figs/seeds_mission_bearing.txt+ [scatter, only marks, mark=*,       scatter src=explicit] table [header=false, x index=1,y index=2, meta index=0] ;
+ [scatter, only marks, mark=square*, scatter src=explicit] table [header=false, x index=3,y index=4, meta index=0] ;

colormap name = bluered,
cycle list=[of colormap],
x dir=reverse,
axis lines*=none,
axis line style = draw=transparent,line width=0.0001pt,
xticklabel style=anchor=--90,
axis line style=-stealth,
legend entries=-   deg,
legend style=font=,
legend cell align=right,
figs/seeds_mission_bearing.txt+ [scatter, only marks, mark=pentagon*,  scatter src=explicit] table [header=false, x index=5,y index=4, meta index=0] ;
[black,dashed] coordinates (+10,+0.0) (+10,+2.2);
[black,dashed] coordinates (-10,+0.0) (-10,+2.2);

Randomness in scenario #1 mission initial and final bearings

Figures <ref> and <ref> are intended to show the wide range of conditions covered by the scenario #1 mission. The left hand plot of figure <ref> shows the initial bearing  (circles) and final bearing  (squares) for each of the different  executions or seeds  (one color per seed), while the right hand plot shows the bearing change (pentagons) at the turn, with dashed lines representing the restrictions listed in table <ref>. The simulation hence covers a wide array of cases, from those in which the aircraft barely changes its bearing after loosing the GNSS signals to those in which it turns around and returns with a bearing close to opposite the initial one. Note that the radius has no meaning in this figure, and it is only intended to improve the visualization of the different bearings.

colormap name = bluered,
cycle list=[of colormap],
xmin=24, xmax=34, xtick=24,26,...,34,
xlabel=  m/ sec,
ymin=+1650, ymax=+3650, ytick=+1650,+2050,...,+3650,
ylabel=  m,
axis lines=left,
axis line style=-stealth,
legend entries=INI, END,
legend style=font=,
legend cell align=left,
figs/seeds_mission_alt_speed.txt+ [scatter, only marks, mark=*,       scatter src=explicit] table [header=false, x index=4,y index=1, meta index=0] ;
+ [scatter, only marks, mark=square*, scatter src=explicit] table [header=false, x index=5,y index=2, meta index=0] ;

colormap name = bluered,
cycle list=[of colormap],
xmin=-4, xmax=4, xtick=-4,-3,...,4,
xlabel=-   m/ sec,
ymin=-750, ymax=+750, ytick=-750,-500,...,750,
ylabel=-   m,
axis lines=left,
axis line style=-stealth,
figs/seeds_mission_alt_speed.txt+ [scatter, only marks, mark=pentagon*, scatter src=explicit] table [header=false, x index=6,y index=3, meta index=0] ;
[black,dashed] coordinates (-8,+100) (-0.5,+100) (-0.5,+750);
[black,dashed] coordinates (-8,-100) (-0.5,-100) (-0.5,-750);
[black,dashed] coordinates (+8,+100) (+0.5,+100) (+0.5,+750);
[black,dashed] coordinates (+8,-100) (+0.5,-100) (+0.5,-750);

Randomness in scenario #1 mission initial and final airspeed and pressure altitude

The left hand plot of figure <ref> shows the initial (circle) and final (square) pressure altitude and airspeed of the scenario #1 mission for each of the different  executions, while the right hand plot shows  the changes (pentagons) included in the maneuvers described above, with dashed lines representing the restrictions listed in table <ref>. As in the previous case, the simulation includes cases with small changes in altitude and airspeed, and others with very significant variations. Note that the default airspeed value of 29   m/sec shown in table <ref>, as well as the low and high speed limits of 24 and 34   m/sec, have been selected according to the aircraft performances, while the default pressure altitude of 2700   m is a compromise between the thrust force diminution with altitude experienced by piston engine aircraft and the need to fly at a sufficiently elevated height over the terrain to ensure that the images generated by the  tool are sharp.

cycle list=red,no markers,very thick,
blue,no markers,very thick,
orange!50!yellow,no markers,very thick,
violet,no markers,very thick,
green,no markers,very thick,
magenta,no markers,very thick,
olive,no markers,very thick,
x dir=reverse,
axis lines*=none,
axis line style = draw=transparent,line width=0.0001pt,
xticklabel style=anchor=--90,

axis line style=-stealth,

figs/error_filter_pos_hor_m.txttable [header=false, x index= 1,y index= 0] ;
table [header=false, x index= 5,y index= 4] ;
table [header=false, x index= 9,y index= 8] ;
table [header=false, x index=13,y index=12] ;
table [header=false, x index=17,y index=16] ;

table [header=false, x index=21,y index=20] ;
table [header=false, x index=25,y index=24] ;
table [header=false, x index=29,y index=28] ;
table [header=false, x index=33,y index=32] ;
table [header=false, x index=37,y index=36] ;

table [header=false, x index=41,y index=40] ;
table [header=false, x index=45,y index=44] ;
table [header=false, x index=49,y index=48] ;
table [header=false, x index=53,y index=52] ;
table [header=false, x index=57,y index=56] ;

table [header=false, x index=61,y index=60] ;
table [header=false, x index=65,y index=64] ;
table [header=false, x index=69,y index=68] ;
table [header=false, x index=73,y index=72] ;
table [header=false, x index=77,y index=76] ;

table [header=false, x index=81,y index=80] ;
table [header=false, x index=85,y index=84] ;
table [header=false, x index=89,y index=88] ;
table [header=false, x index=93,y index=92] ;
table [header=false, x index=97,y index=96] ;

table [header=false, x index=101,y index=100] ;
table [header=false, x index=105,y index=104] ;
table [header=false, x index=109,y index=108] ;
table [header=false, x index=113,y index=112] ;
table [header=false, x index=117,y index=116] ;

table [header=false, x index=121,y index=120] ;
table [header=false, x index=125,y index=124] ;
table [header=false, x index=129,y index=128] ;
table [header=false, x index=133,y index=132] ;
table [header=false, x index=137,y index=136] ;

table [header=false, x index=141,y index=140] ;
table [header=false, x index=145,y index=144] ;
table [header=false, x index=149,y index=148] ;
table [header=false, x index=153,y index=152] ;
table [header=false, x index=157,y index=156] ;

table [header=false, x index=161,y index=160] ;
table [header=false, x index=165,y index=164] ;
table [header=false, x index=169,y index=168] ;
table [header=false, x index=173,y index=172] ;
table [header=false, x index=177,y index=176] ;

table [header=false, x index=181,y index=180] ;
table [header=false, x index=185,y index=184] ;
table [header=false, x index=189,y index=188] ;
table [header=false, x index=193,y index=192] ;
table [header=false, x index=197,y index=196] ;

table [header=false, x index=201,y index=200] ;
table [header=false, x index=205,y index=204] ;
table [header=false, x index=209,y index=208] ;
table [header=false, x index=213,y index=212] ;
table [header=false, x index=217,y index=216] ;

table [header=false, x index=221,y index=220] ;
table [header=false, x index=225,y index=224] ;
table [header=false, x index=229,y index=228] ;
table [header=false, x index=233,y index=232] ;
table [header=false, x index=237,y index=236] ;

table [header=false, x index=241,y index=240] ;
table [header=false, x index=245,y index=244] ;
table [header=false, x index=249,y index=248] ;
table [header=false, x index=253,y index=252] ;
table [header=false, x index=257,y index=256] ;

table [header=false, x index=261,y index=260] ;
table [header=false, x index=265,y index=264] ;
table [header=false, x index=269,y index=268] ;
table [header=false, x index=273,y index=272] ;
table [header=false, x index=277,y index=276] ;

table [header=false, x index=281,y index=280] ;
table [header=false, x index=285,y index=284] ;
table [header=false, x index=289,y index=288] ;
table [header=false, x index=293,y index=292] ;
table [header=false, x index=297,y index=296] ;

table [header=false, x index=301,y index=300] ;
table [header=false, x index=305,y index=304] ;
table [header=false, x index=309,y index=308] ;
table [header=false, x index=313,y index=312] ;
table [header=false, x index=317,y index=316] ;

table [header=false, x index=321,y index=320] ;
table [header=false, x index=325,y index=324] ;
table [header=false, x index=329,y index=328] ;
table [header=false, x index=333,y index=332] ;
table [header=false, x index=337,y index=336] ;

table [header=false, x index=341,y index=340] ;
table [header=false, x index=345,y index=344] ;
table [header=false, x index=349,y index=348] ;
table [header=false, x index=353,y index=352] ;
table [header=false, x index=357,y index=356] ;

table [header=false, x index=361,y index=360] ;
table [header=false, x index=365,y index=364] ;
table [header=false, x index=369,y index=368] ;
table [header=false, x index=373,y index=372] ;
table [header=false, x index=377,y index=376] ;

table [header=false, x index=381,y index=380] ;
table [header=false, x index=385,y index=384] ;
table [header=false, x index=389,y index=388] ;
table [header=false, x index=393,y index=392] ;
table [header=false, x index=397,y index=396] ;

GNSS-Based horizontal tracks for scenario #1 seeds Υ_1 through Υ_100

Figure <ref> shows the horizontal tracks obtained with the aid of GNSS when executing the = 100 runs of scenario #1. Note that the variability in the flown distance is caused by both the airspeed variations shown in table <ref> as well as those of the wind field described in table <ref> below.

§.§.§ Weather and Wind

The simulation adopts the INSA model [13] to define the atmospheric conditions encountered by the aircraft during flight, in which the temperature and pressure offsets (,   ) are provided as functions of horizontal position and time. Similarly, the low frequency wind field is also defined as two user provided functions[Both north and east wind, as well as magnitude and direction can be employed.] based on horizontal position, altitude, and time. In order to properly represent the influence that weather and wind field changes may have in the resulting GNSS-Denied trajectory, scenario #1 is defined as containing constant weather and wind field both at the beginning (, , ) and at the end (, , ) of the trajectory, with linear variations in between, as shown in figure <ref>.

(origin)     at (+0.0,+0.0);
(tGNSS)      at (+0.5,+0.0);
(tINIDeltaT) at (+2.0,+0.0);
(tENDDeltaT) at (+5.9,+0.0);
(tINIDeltap) at (+2.6,+0.0);
(tENDDeltap) at (+4.7,+0.0);
(tINIWIND)   at (+3.3,+0.0);
(tENDWIND)   at (+6.8,+0.0);
(tEND)       at (+8.0,+0.0);

(DeltaTINI)  at (+0.0,+0.2);
(DeltaTEND)  at (+0.0,+2.7);
(DeltapINI)  at (+0.0,+2.4);
(DeltapEND)  at (+0.0,+0.9);
(vwindINI)   at (+0.0,+0.6);
(vwindEND)   at (+0.0,+2.2);
(chiWINDINI) at (+0.0,+1.6);
(chiWINDEND) at (+0.0,+1.3);
(yEND)       at (+0.0,+2.9);

(corner)       at ($(tEND)       + (yEND)$);
(XXDeltaTINI)  at ($(tINIDeltaT) + (DeltaTINI)$);
(XXDeltaTEND)  at ($(tENDDeltaT) + (DeltaTEND)$);
(XXDeltaT111)  at ($(tEND)       + (DeltaTEND)$);
(YYDeltapINI)  at ($(tINIDeltap) + (DeltapINI)$);
(YYDeltapEND)  at ($(tENDDeltap) + (DeltapEND)$);
(YYDeltap111)  at ($(tEND)       + (DeltapEND)$);
(ZZvwindINI)   at ($(tINIWIND)   + (vwindINI)$);
(ZZvwindEND)   at ($(tENDWIND)   + (vwindEND)$);
(ZZvwind111)   at ($(tEND)       + (vwindEND)$);
(TTchiWINDINI) at ($(tINIWIND)   + (chiWINDINI)$);
(TTchiWINDEND) at ($(tENDWIND)   + (chiWINDEND)$);
(TTchiWIND111) at ($(tEND)       + (chiWINDEND)$);

[name] [->] (origin) – (tEND);
[name] [->] (origin) – (yEND);
[name] [-]  (tEND) – (corner);
[name] [-]  (yEND) – (corner);
node at($(tGNSS)     +(+0.0,-0.2)$) ;
node at($(tINIDeltaT)+(+0.0,-0.2)$) ;
node at($(tENDDeltaT)+(+0.0,-0.2)$) ;
node at($(tINIDeltap)+(+0.0,-0.2)$) ;
node at($(tENDDeltap)+(+0.0,-0.2)$) ;
node at($(tINIWIND)  +(+0.0,-0.2)$) ;
node at($(tENDWIND)  +(+0.0,-0.2)$) ;
node at($(tEND)      +(+0.0,-0.2)$) ;
node at($(DeltaTINI) +(-0.3,-0.0)$) ;
node at($(DeltaTEND) +(-0.35,-0.0)$) ;
node at($(DeltapINI) +(-0.3,-0.0)$) ;
node at($(DeltapEND) +(-0.35,-0.0)$) ;
node at($(vwindINI)  +(-0.4,-0.0)$) ;
node at($(vwindEND)  +(-0.45,-0.0)$) ;
node at($(chiWINDINI)+(-0.45,-0.0)$) ;
node at($(chiWINDEND)+(-0.5,-0.0)$) ;

[black,  ultra thin, dashed] [-] (tGNSS)      – ($(tGNSS) + (yEND)$);
[black,  ultra thin, dashed] [-] (tINIDeltaT) – (XXDeltaTINI);
[black,  ultra thin, dashed] [-] (tENDDeltaT) – (XXDeltaTEND);
[black,  ultra thin, dashed] [-] (tINIDeltap) – (YYDeltapINI);
[black,  ultra thin, dashed] [-] (tENDDeltap) – (YYDeltapEND);
[black,  ultra thin, dashed] [-] (tINIWIND)   – (ZZvwindINI);
[black,  ultra thin, dashed] [-] (tENDWIND)   – (ZZvwindEND);
[black,  ultra thin, dashed] [-] (tINIWIND)   – (TTchiWINDINI);
[black,  ultra thin, dashed] [-] (tENDWIND)   – (TTchiWINDEND);
[black,  ultra thin, dashed] [-] (DeltaTEND)  – (XXDeltaTEND);
[black,  ultra thin, dashed] [-] (DeltapEND)  – (YYDeltapEND);
[black,  ultra thin, dashed] [-] (vwindEND)   – (ZZvwindEND);
[black,  ultra thin, dashed] [-] (chiWINDEND) – (TTchiWINDEND);

[blue,    name path = DeltaT, thick] [-] (DeltaTINI) – (XXDeltaTINI) – (XXDeltaTEND) – (XXDeltaT111);
[red,     name path = Deltap, thick] [-] (DeltapINI) – (YYDeltapINI) – (YYDeltapEND) – (YYDeltap111);
[green,   name path = vwind,  thick] [-] (vwindINI)  – (ZZvwindINI)  – (ZZvwindEND)  – (ZZvwind111);
[magenta, name path = chiWIND,thick] [-] (chiWINDINI)– (TTchiWINDINI)– (TTchiWINDEND)– (TTchiWIND111);

Schematic view of scenario #1 weather and wind field

The values determining the initial and final weather and wind field, together with the times at which the linear variations start and conclude, are defined stochastically as shown in table <ref>.

Component     Variable     Definition     Unit     Restrictions

Temp ramp initial time                     N400, 600^2                    [sec]                     > 50

Temp ramp final time                       N+ 1200, 600^2     [sec]                     > + 600

Initial temperature offset                  N0, 10^2                       ^∘K

Final temperature offset                    N, 3^2               ^∘K

Pressure ramp initial time                 N400, 600^2                    [sec]                     > 50

Pressure ramp final time                   N+ 1200, 600^2     [sec]                     > + 600

Initial pressure offset                     N0, 1500^2                     [pa]

Final pressure offset                       N, 300^2             [pa]

Wind ramp initial time                       N400, 600^2                    [sec]                     > 50

Wind ramp final time                         N+ 1200, 600^2     [sec]                     > + 300

Initial wind speed			                 N0, 7^2		   		         [m/sec]

Final wind speed			                 N, 3^2			     [m/sec]

Initial wind bearing                       U-179, 180				     [deg]

Final wind bearing                         N, 15^2		     [deg]

tableDeterministic and stochastic components of scenario #1 weather and wind field

The values defining the weather changes have been selected by the author so they are significantly more intense than what a normal flight would experience, with the intention of highlighting any negative effects that weather and wind field changes may have in the accuracy of the GNSS-Denied navigation algorithms. On average they imply a temperature change of 3   ^∘K[Note that this temperature change does not include altitude changes and is exclusively restricted to atmosphere warming or cooling. This is intended to represent positive changes during the morning and negative during the evenings.], a pressure change of 300   pa[As in the temperature case, this pressure change does not include altitude changes, and represents the change produced by the aircraft flying into a high or low pressure weather system. To provide more clarity, note that 108600   pa is the highest sea level pressure ever recorded, 103000   pa corresponds to a strong high pressure system, 101325   pa is the standard sea level pressure [11], 100000   pa corresponds to a typical low pressure system, 95000   pa to a category III hurricane, and 87000   pa is the lowest ever recorded.], and a wind change of 3   m/sec and 15   deg in bearing, all over 20   min.

colormap name = bluered,
cycle list=[of colormap],
xmin=-5000, xmax=5000, xtick=-5000,-3000,...,5000,
ymin=-25, ymax=+26, ytick=-25,-15,...,25,
ylabel=  ^∘K,
axis lines=left,
axis line style=-stealth,
legend entries=,   , ,   ,
legend style=font=,
legend cell align=center,
legend pos=north west,
figs/seeds_weather_offsets.txt+ [scatter, only marks, mark=*,       scatter src=explicit] table [header=false, x index=7,y index=2, meta index=0] ;
+ [scatter, only marks, mark=square*, scatter src=explicit] table [header=false, x index=9,y index=4, meta index=0] ;

colormap name = bluered,
cycle list=[of colormap],
xmin=-800, xmax=1000, xtick=-800,-600,...,1000,
xlabel=- pa,
ymin=-6, ymax=+6, ytick=-6,-4,...,6,
ylabel=-   ^∘K,
axis lines=left,
axis line style=-stealth,
figs/seeds_weather_offsets.txt+ [scatter, only marks, mark=pentagon*,  scatter src=explicit] table [header=false, x index=10,y index=5, meta index=0] ;

Randomness in scenario #1 temperature and pressure offsets

The seeds  and  corresponding to each trajectory seed  are the only inputs required to obtain all the values listed in table <ref>, hence ensuring that the results are repeatable if the same seed is employed again. Figure <ref> showcases the wide variety of weather conditions encountered during the simulation. As in previous figures, the left hand plot shows the initial and final temperature and pressure offsets (, ) for each of the  different executions, while the right hand plot shows the changes in offsets between the start and the end of the trajectory. Note that although most values can be considered normal, some executions involve weather conditions that deviate substantially from standard ones, which is the objective behind the values selected for table <ref>.

colormap name = bluered,
cycle list=[of colormap],
x dir=reverse,
axis lines*=none,
axis line style = draw=transparent,line width=0.0001pt,
xticklabel style=anchor=--90,
ymin=0, ymax=20.0, ytick=0,5,...,15,
legend entries=  m/sec,deg,   m/sec,deg,
legend style=font=,
legend cell align=right,
figs/seeds_weather_wind.txt+ [scatter, only marks, mark=*,       scatter src=explicit] table [header=false, x index=5,y index=3, meta index=0] ;
+ [scatter, only marks, mark=square*, scatter src=explicit] table [header=false, x index=6,y index=4, meta index=0] ;

colormap name = bluered,
cycle list=[of colormap],
x dir=reverse,
axis lines*=none,
axis line style = draw=transparent,line width=0.0001pt,
xticklabel style=anchor=--90,
ymin=0, ymax=8.0, ytick=0,2,...,6,
axis line style=-stealth,
legend entries=|- |,
legend style=font=,
legend cell align=right,
figs/seeds_weather_wind.txt+ [scatter, only marks, mark=pentagon*,  scatter src=explicit] table [header=false, x index=8,y index=7, meta index=0] ;

Randomness in scenario #1 wind speeds

Figure <ref> shows the wind field for each of the different  executions. As in previous figures, the left hand plot represents the initial and final wind speed, while the right hand plots shows the wind change experienced by the aircraft during its flight. In contrast with figure <ref>, these plots show both the wind field intensity and bearing. Note the relatively high wind speeds and wind changes experienced by some executions when compared with the aircraft airspeed values discussed in table <ref>.

The stochastic nature of the scenario #1 weather and wind fields does not only reside in their initial and final values represented in figures <ref> and <ref>, but in the velocity at which the weather and wind field transition from the initial values to their final ones. Although not plotted, the velocity of these changes is specified in table <ref>.

§.§.§ Terrain Type

The type of terrain overflown by the aircraft has a significant influence in the performance of the visual navigation algorithms. In particular, the terrain texture (or lack of) and its elevation relief are the two most important characteristics in this regard. For this reason, each of the scenario #1  Monte Carlo simulation runs is repeated four times with the only change of the initial geodetic coordinates, so the resulting flights occur in four different zones or types of terrain. The zones described below are intended to represent a wide array of terrain types[Alpine terrain with elevated vertical relief is not included as piston engines aircraft can not fly high enough for the previously described  tool to generate realistic images in this type of terrain.]. Note that changing the area of the world were the flight takes place also modifies the gravitational and magnetic fields, but as all considered zones are far from the Earth magnetic poles, this does not imply any significant change in the performance of the inertial algorithms. Images representative of each zone as viewed by the onboard camera are also included.

* The “desert” (DS) zone is located in the Sonoran desert of southern Arizona (USA) and northern Mexico, with an initial position given by λ_DS = 248.001185 deg, φ_DS = 32.157903 deg, and h_GROUND,DS = 661.0 m. It is characterized by a combination of bajadas (broad slopes of debris) and isolated very steep mountain ranges. There is virtually no human infrastructure or flat terrain, as the bajadas have sustained slopes of up to 7 [deg]. The altitude of the bajadas ranges from 300 to 800 [m] above , and the mountains reach up to 800 [m] above the surrounding terrain. Texture is abundant because of the cacti and the vegetation along the dry creeks.

Typical “desert” (DS) terrain view

* The “farm” (FM) zone is located in the fertile farmland of southeastern Illinois and southwestern Indiana (USA), with an initial position given by λ_FM = 272.122371 deg, φ_FM = 38.865625 deg, and h_GROUND,FM = 144.0 m. A significant percentage of the terrain is made of regular plots of farmland, but there also exists some woodland, farm houses, rivers, lots of little towns, and roads. It is mostly flat with an altitude above  between 100 and 200 [m], and altitude changes are mostly restricted to the few forested areas. Texture is nonexistent in the farmlands, where extracting features is often impossible.

Typical “farm” (FM) terrain view

* The “forest” (FR) zone is located in the deciduous forestlands of Vermont and New Hamspshire (USA), with an initial position given by λ_FR = 287.490805 deg, φ_FR = 43.354486 deg, and h_GROUND,FR = 200.0 m. The terrain is made up of forests and woodland, with some clearcuts, small towns, and roads. There are virtually no flat areas, as the land is made up by hills and small to medium size mountains that are never very steep. The valleys range from 100 to 300 [m] above , while the tops of the mountains reach 500 to 900 [m]. Features are plentiful in the woodlands.

Typical “forest” (FR) terrain view

* The “mix” (MX) zone is located in northern Mississippi and extreme southwestern Tennessee (USA), with an initial position given by λ_MX = 270.984538 deg, φ_MX = 34.720636 deg, and h_GROUND,MX = 133.0 m. Approximately half of the land consists of woodland in the hills, and the other half is made up by farmland in the valleys, with a few small towns and roads. Altitude changes are always presents and the terrain is never flat, but they are smaller than in the DS and FR zones, with the altitude oscillating between 100 and 200 [m] above .

Typical “mix” (MX) terrain view

§ SCENARIO #2

Scenario #1 is intended to represent the challenges involved in reaching a predefined recovery location without the usage of GNSS signals and while experiencing dynamic weather and wind conditions. Although it is quite generic, there are two reasons why it can not constitute the only scenario employed to evaluate the behavior of different GNSS-Denied navigation algorithms. On one side, although the scenario #1 mission includes three maneuvers (a change of bearing, a change of airspeed, and a change of pressure altitude), for the most part it consists of a long straight level flight in which the aircraft systems have lots of time to recover from estimation errors induced by the maneuvers, which are executed quite apart from each other. Scenario #2 includes several maneuvers and avoids long straight level segments to ensure that the GNSS-Denied navigation algorithms are valid for different types of missions.

Additionally, the GNSS-Denied position errors incurred by some inertial navigation systems may exhibit a significant dependency on the weather and wind changes that occur since the aircraft can no longer employ the GNSS signals for navigation. The prevalence of these errors coupled with the significant scenario #1 weather and wind changes prevents the detection of smaller error sources when analyzing the Monte Carlo simulation results.

These two reasons indicate the need for a scenario #2, which represents the challenges involved in continuing with the original mission upon the loss of the GNSS signals, executing a series of continuous maneuvers over a relatively short period of time with no weather or wind variations.

§.§.§ Mission

The mission, guidance targets, or reference trajectory of scenario #2 consists in the seventeen flight segments defined in (<ref>) through (<ref>) below. It is intended to showcase the ability of an autonomous aircraft to continue executing its mission while in GNSS-Denied conditions, specially in a congested, urban, or mountainous environment in which continuous maneuvers are required to avoid other aircraft, infrastructure, or terrain.

* The whole trajectory is executed at constant airspeed , constant pressure altitude , and no sideslip β= 0.

* The first segment coincides with that of scenario #1, as the aircraft flies at constant bearing , first with the assistance of GNSS signals until = 100 sec and then in GNSS-Denied conditions until t_TURN,1.

* The rest of the mission consists in a series of eight turning maneuvers connecting straight segments flown at constant bearings χ_i. Turns are executed with body bank angles of = ±10   deg, with the sign given by the shortest turn necessary to reach the next bearing χ_i, upon which the turn concludes, and the duration of the straight segments is given by . The short time between consecutive maneuvers is intended to avoid providing the inertial navigation system with sufficient time to stabilize itself in case the maneuvers have a destabilizing effect on the estimation errors.

* The scenario concludes at = 500   sec, which for some seeds may not be sufficient to conclude the eight turns.

\begin{eqnarray}
\nm{\xREF}                 & = & \nm{\lrb{\deltaTARGETone, \, \deltaTARGETtwo, \dots, \deltaTARGETseventeen}} \label{eq:Sim_Scenario_Alter_Mission_RT} \\
\nm{\deltaTARGETone}       & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \chiINI, \, \beta = 0, \, t_{{\sss TURN},1}}^T}\label{eq:Sim_Scenario_Alter_Mission_target1} \\
\nm{\deltaTARGETtwo}       & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \xi_1,   \, \beta = 0, \, \chi_1}^T}\label{eq:Sim_Scenario_Alter_Mission_target2} \\
\nm{\deltaTARGETthr}       & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \chi_1,  \, \beta = 0, \, \Delta t_{{\sss TURN},2}}^T}\label{eq:Sim_Scenario_Alter_Mission_target3} \\
\nm{\cdots}                &   & \nm{\cdots}\nonumber \\
\nm{\deltaTARGETfourteen}  & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \xi_7,   \, \beta = 0, \, \chi_7}^T}\label{eq:Sim_Scenario_Alter_Mission_target14} \\
\nm{\deltaTARGETfifteen}   & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \chi_7,  \, \beta = 0, \, \Delta t_{{\sss TURN},8}}^T}\label{eq:Sim_Scenario_Alter_Mission_target15} \\
\nm{\deltaTARGETsixteen}   & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \xi_8,   \, \beta = 0, \, \chi_8}^T}\label{eq:Sim_Scenario_Alter_Mission_target16} \\
\nm{\deltaTARGETseventeen} & = & \nm{\lrsb{\vtasINI, \, \HpINI, \,  \chi_8,  \, \beta = 0, \, \tEND}^T}\label{eq:Sim_Scenario_Alter_Mission_target17}
\end{eqnarray}

Table <ref> shows the definition of the deterministic and stochastic variables on which the scenario #2 mission relies. Note that it has been developed with the objective of including as much variability as possible, as the only deterministic variables are those that do not result in loss of generality. The initial bearing, altitude, and airspeed coincide with those of scenario #1 shown in the left hand side of figures <ref> and <ref>, represented by circles.

Component     Variable     Index     Definition     Unit     Restrictions

Airspeed                                                         N29, 1.5^2         [m/sec]     24 < < 34

Pressure altitude	                                               N2700, 200^2       [m]

Initial bearing           = χ_0                               U-179, 180         [deg]

Sideslip angle            β                                          0                        [deg]

GNSS-Denied time	      		                                    100		   		         [sec]

Turn bank angles	      	           i = 1, …, 8     ±10			         [deg]

Intermediate bearings     χ_i                i = 1, …, 7     U-179, 180         [deg]       |χ_i - χ_i-1 |> 10

Final bearing             = χ_8                               U-179, 180         [deg]       |χ_8 - χ_7 |> 10

Initial turn time	      t_TURN,1                              + N(30, 50^2)     [sec]       > + 15

Turn time intervals                 i = 2, …, 8     U10, 50            [sec]

Final time                                                          500					     [sec]

tableDeterministic and stochastic components of the scenario #2 mission

cycle list=red,no markers,very thick,
blue,no markers,very thick,
orange!50!yellow,no markers,very thick,
violet,no markers,very thick,
green,no markers,very thick,
magenta,no markers,very thick,
olive,no markers,very thick,
x dir=reverse,
axis lines*=none,
axis line style = draw=transparent,line width=0.0001pt,
xticklabel style=anchor=--90,

axis line style=-stealth,

figs/error_filter_pos_alter_hor_m.txttable [header=false, x index= 1,y index= 0] ;
table [header=false, x index= 5,y index= 4] ;
table [header=false, x index= 9,y index= 8] ;
table [header=false, x index=13,y index=12] ;
table [header=false, x index=17,y index=16] ;

table [header=false, x index=21,y index=20] ;
table [header=false, x index=25,y index=24] ;
table [header=false, x index=29,y index=28] ;
table [header=false, x index=33,y index=32] ;
table [header=false, x index=37,y index=36] ;

table [header=false, x index=41,y index=40] ;
table [header=false, x index=45,y index=44] ;
table [header=false, x index=49,y index=48] ;
table [header=false, x index=53,y index=52] ;
table [header=false, x index=57,y index=56] ;

table [header=false, x index=61,y index=60] ;
table [header=false, x index=65,y index=64] ;
table [header=false, x index=69,y index=68] ;
table [header=false, x index=73,y index=72] ;
table [header=false, x index=77,y index=76] ;

table [header=false, x index=81,y index=80] ;
table [header=false, x index=85,y index=84] ;
table [header=false, x index=89,y index=88] ;
table [header=false, x index=93,y index=92] ;
table [header=false, x index=97,y index=96] ;

table [header=false, x index=101,y index=100] ;
table [header=false, x index=105,y index=104] ;
table [header=false, x index=109,y index=108] ;
table [header=false, x index=113,y index=112] ;
table [header=false, x index=117,y index=116] ;

table [header=false, x index=121,y index=120] ;
table [header=false, x index=125,y index=124] ;
table [header=false, x index=129,y index=128] ;
table [header=false, x index=133,y index=132] ;
table [header=false, x index=137,y index=136] ;

table [header=false, x index=141,y index=140] ;
table [header=false, x index=145,y index=144] ;
table [header=false, x index=149,y index=148] ;
table [header=false, x index=153,y index=152] ;
table [header=false, x index=157,y index=156] ;

table [header=false, x index=161,y index=160] ;
table [header=false, x index=165,y index=164] ;
table [header=false, x index=169,y index=168] ;
table [header=false, x index=173,y index=172] ;
table [header=false, x index=177,y index=176] ;

table [header=false, x index=181,y index=180] ;
table [header=false, x index=185,y index=184] ;
table [header=false, x index=189,y index=188] ;
table [header=false, x index=193,y index=192] ;
table [header=false, x index=197,y index=196] ;

table [header=false, x index=201,y index=200] ;
table [header=false, x index=205,y index=204] ;
table [header=false, x index=209,y index=208] ;
table [header=false, x index=213,y index=212] ;
table [header=false, x index=217,y index=216] ;

table [header=false, x index=221,y index=220] ;
table [header=false, x index=225,y index=224] ;
table [header=false, x index=229,y index=228] ;
table [header=false, x index=233,y index=232] ;
table [header=false, x index=237,y index=236] ;

table [header=false, x index=241,y index=240] ;
table [header=false, x index=245,y index=244] ;
table [header=false, x index=249,y index=248] ;
table [header=false, x index=253,y index=252] ;
table [header=false, x index=257,y index=256] ;

table [header=false, x index=261,y index=260] ;
table [header=false, x index=265,y index=264] ;
table [header=false, x index=269,y index=268] ;
table [header=false, x index=273,y index=272] ;
table [header=false, x index=277,y index=276] ;

table [header=false, x index=281,y index=280] ;
table [header=false, x index=285,y index=284] ;
table [header=false, x index=289,y index=288] ;
table [header=false, x index=293,y index=292] ;
table [header=false, x index=297,y index=296] ;

table [header=false, x index=301,y index=300] ;
table [header=false, x index=305,y index=304] ;
table [header=false, x index=309,y index=308] ;
table [header=false, x index=313,y index=312] ;
table [header=false, x index=317,y index=316] ;

table [header=false, x index=321,y index=320] ;
table [header=false, x index=325,y index=324] ;
table [header=false, x index=329,y index=328] ;
table [header=false, x index=333,y index=332] ;
table [header=false, x index=337,y index=336] ;

table [header=false, x index=341,y index=340] ;
table [header=false, x index=345,y index=344] ;
table [header=false, x index=349,y index=348] ;
table [header=false, x index=353,y index=352] ;
table [header=false, x index=357,y index=356] ;

table [header=false, x index=361,y index=360] ;
table [header=false, x index=365,y index=364] ;
table [header=false, x index=369,y index=368] ;
table [header=false, x index=373,y index=372] ;
table [header=false, x index=377,y index=376] ;

table [header=false, x index=381,y index=380] ;
table [header=false, x index=385,y index=384] ;
table [header=false, x index=389,y index=388] ;
table [header=false, x index=393,y index=392] ;
table [header=false, x index=397,y index=396] ;

GNSS-Based horizontal tracks for scenario #2 seeds Υ_1 through Υ_100

Figure <ref> shows the horizontal tracks obtained with the aid of GNSS when executing the = 100 runs of scenario #2.

§.§.§ Weather and Wind

The influence of the weather and wind changes since the time the GNSS signals are lost is properly captured by scenario #1. In scenario #2, the weather and wind remain constant throughout the trajectory, which enables the identification of other navigation error sources that are obscured in scenario #1 by the prevalence of the weather and wind changes. In contrast with scenario #1, which was designed to represent with as much realism as possible the challenges encountered by an autonomous aircraft trying to reach a predefined recovery location, scenario #2 is unrealistic as it completely eliminates the most important error sources for the inertial estimation of the aircraft position. A minimum amount of weather and wind changes is always present in real flights, so employing scenario #2 will always results in predictions that are optimistic.

Component     Variable     Definition     Unit

Temperature offset                  N0, 10^2		    ^∘K

Pressure offset                     N0, 1500^2	    [pa]

Wind speed			                 N0, 7^2		    [m/sec]

Wind bearing                       U-179, 180	    [deg]

tableDeterministic and stochastic components of scenario #2 weather and wind field

Scenario #2 hence is defined as containing constant weather and wind field (, , ). These values are determined stochastically as shown in table <ref>, and coincide with the initial conditions of scenario #1 (table <ref>). These are shown in the left hand side of figures <ref> and <ref>, represented by circles.

§.§.§ Terrain Type

The short duration and continuous maneuvering of scenario #2 enables the use of two additional zones or terrain types in its Monte Carlo simulations. These two zones are not employed in scenario #1 because the author could not locate wide enough areas with a prevalence of this type of terrain (note that scenario #1 trajectories can conclude up to 125 [km] in any direction from its initial coordinates, but only 12 [km] for scenario #2, as shown in figures <ref> and <ref>).

* The “prairie” (PR) zone is located in the Everglades floodlands of southern Florida (USA), with an initial position given by λ_PR = 279.088834 deg, φ_PR = 25.855172 deg, and h_GROUND,PR = 10.0 m. It consists of flat grasslands, swamps, and tree islands located a few meters above , with the only human infrastructure being a few dirt roads and landing strips, but no settlements. Features may be difficult to obtain in some areas due to the lack of texture.

Typical “prairie” (PR) terrain view

* The “urban” (UR) zone is located in the Los Angeles metropolitan area, with an initial position given by λ_UR = 241.799731 deg, φ_UR = 33.924426 deg, and h_GROUND,UR = 26.0 m. It is composed by a combination of single family houses and commercial buildings separated by freeways and streets. There is some vegetation but no natural landscapes, and the terrain is flat and close to .

Typical “urban” (UR) terrain view

Although scenario #2 also makes use of the four terrain types listed for scenario #1 (DS, FM, FR, and MX), it is worth noting that the variability of the terrain is significantly higher for scenario #1 because of the bigger land extension covered. The altitude relief, abundance or scarcity of features, land use diversity, and presence of rivers and mountains is hence more varied when executing a given run of scenario #1 on a certain type of terrain, than when executing the same run for scenario #2. From the point of view of the influence of the terrain on the visual navigation algorithms, scenario #1 should theoretically be more challenging than #2.

§ METRICS

This section defines the metrics employed to evaluate the performance of a given navigation algorithm in GNSS-Denied conditions based on the results obtained in Monte Carlo simulations applied to the two scenarios defined above. To be useful, the metrics employed must be as simple as possible so as to facilitate the comparison process between two different algorithms, but also be representative enough of the performances that the navigation system should exhibit.

Let's imagine an error function f_X ,   ,   t = f_X ,   ,   t that evaluates the difference between a certain variable of the observed trajectory X_EST and the same variable in the actual trajectory X_TRUTH, both at the same time t:
f_X ,   ,   t = f_X t = X_ESTt - X_TRUTHt eq:Sim_Metrics_function_f

This error function is evaluated simultaneously with the navigation system outputs (t = t_n = n ·). A good navigation system is that which diminishes the navigation system error, so the output of the error function f_X is as small as possible. Note that for each scenario #1 execution or seed , there exist = / + 1 = 380001 evaluations of the error function f_X, which is renamed f_Xj to specify that it applies to execution . In the case of scenario #2, the number of evaluations is 50001.
f_Xj _j,   _j,   t_n = f_Xj t_n = X_ESTjt_n - X_TRUTHjt_n eq:Sim_Metrics_function2_f

§.§.§ Trajectory Metrics

The error function f_Xj can be considered a random or stochastic process, and hence its mean X̂ and standard deviation X̂ can be readily computed. In addition, it is also interesting to obtain X̂, this is, the maximum absolute value[The maximum absolute value metric of a group of values returns the signed value with the highest absolute value.] of the function f_Xj:
X̂  =  ∑_n=1^ f_Xj t_n       √(∑_n=1^ (f_Xj t_n - X̂)^2)       max_n=1^  |f_Xj t_n| ^T eq:Sim_Metrics_mean_std_max_j

The above mean, standard deviation, and absolute value maximum provide a good understanding of how good the navigation system is at continuously estimating X_EST when compared with its real value X_TRUTH for an specific trajectory execution or seed , and for that reason they are known as trajectory metrics or seed metrics.

§.§.§ Final State Value

In some ocassions it may be necessary to ignore the f_Xj error function evolution with time to only focus on its value at the end of the scenario (t_n = ). This value is provided by the final state value:
X̂ = f_Xj eq:Sim_metrics_final_traj

§.§.§ Aggregated Metrics

The behavior of the navigation system can only be properly evaluated when it is applied to a sufficiently big number of executions (), and it is for this purpose that the aggregated metrics (means, standard deviations, and absolute value maximums) are computed[Note that the aggregated mean of maximums and the aggregated standard deviation of maximums only make sense if applied to the absolute value of the maximums instead of the maximum themselves. They are defined accordingly in (<ref>, <ref>).]. To obtain the aggregated metrics, the trajectory mean X̂, standard deviation X̂, and absolute value maximum X̂ for each seed are considered realizations of a random variable:

\begin{eqnarray}
\nm{\begin{bmatrix} \nm{\mumu{\hat{X}}} \\ \nm{\musigma{\hat{X}}} \\ \nm{\mumax{\hat{X}}} \end{bmatrix}} & = & \nm{\begin{bmatrix} \nm{\cfrac{\displaystyle\sum_{j=1}^{\nEX} \muj{\hat{X}}}{\nEX}} & \nm{\cfrac{\displaystyle\sum_{j=1}^{\nEX} \sigmaj{\hat{X}}}{\nEX}} & \nm{\cfrac{\displaystyle\sum_{j=1}^{\nEX}{\lvert \maxj{\hat{X}} \rvert}}{\nEX}} \end{bmatrix}^T} \label{eq:Sim_Metrics_mean_mean_std_max} \\
\nm{\begin{bmatrix} \nm{\sigmamu{\hat{X}}} \\ \nm{\sigmasigma{\hat{X}}} \\ \nm{\sigmamax{\hat{X}}} \end{bmatrix}} & = & \nm{\begin{bmatrix} \nm{\sqrt{\cfrac{\displaystyle\sum_{j=1}^{\nEX} \Big(\muj{\hat{X}} - \mumu{\hat{X}}\Big)^2}{\nEX}}} & \nm{\sqrt{\cfrac{\displaystyle\sum_{j=1}^{\nEX} \Big(\sigmaj{\hat{X}} - \musigma{\hat{X}}\Big)^2}{\nEX}}} & \nm{\sqrt{\cfrac{\displaystyle\sum_{j=1}^{\nEX} \Big({\lvert \maxj{\hat{X}} \rvert} - \mumax{\hat{X}}\Big)^2}{\nEX}}} \end{bmatrix}^T} \label{eq:Sim_Metrics_std_mean_std_max} \\
\nm{\begin{bmatrix} \nm{\maxmu{\hat{X}}} \\ \nm{\maxsigma{\hat{X}}} \\ \nm{\maxmax{\hat{X}}} \end{bmatrix}} & = & \nm{\begin{bmatrix} \nm{\displaystyle\max_{j=1}^{\nEX} \ \lvert \muj{\hat{X}} \rvert} & \nm{\displaystyle\max_{j=1}^{\nEX} \ \lvert \sigmaj{\hat{X}} \rvert} & \nm{\displaystyle\max_{j=1}^{\nEX} \ \lvert \maxj{\hat{X}} \rvert} \end{bmatrix}^T} \label{eq:Sim_Metrics_max_mean_std_max}
\end{eqnarray}

These metrics provide an overall view of the navigation system behavior under different circumstances, as they take into consideration different combinations of the stochastic parameters that define the two scenarios.

§.§.§ Aggregated Final State Metrics

The aggregated metrics are adequate to evaluate the navigation system performance as long as each evaluation of f_Xj t_n can be considered a realization of a stochastic process, which is the case for those variables that can be tracked (with more or less accuracy) without drift or accumulated errors by the navigation system. If a certain variable presents a drift, this is, the estimation error has a tendency to increase with time, then both the trajectory as well as the aggregated metrics are meaningless as the results depend on the execution length . In this case, it makes sense to compute the aggregated final state metrics, which rely on the final state value of every execution :
X̂  =  ∑_j=1^ X̂     √(∑_j=1^ (X̂ - X̂ )^2)     max_j=1^  |X̂ | ^Teq:Sim_Metrics_mean_std_max_end

§.§.§ Time Aggregated Metrics

The time aggregated metrics consider that all evaluations of the error function f_Xj for different runs at a given time instant t_n together make up a random vector f_Xn of size , and then computes its mean X̂ and standard deviation X̂. They are primarily intended for graphical representation:

\begin{eqnarray}
\nm{\mun{\hat{X}}}    & = & \nm{\cfrac{\displaystyle\sum_{j=1}^{\nEX} f_{Xj} \lrp{t_n}}{\nEX}}\label{eq:Sim_Metrics_mean_s} \\
\nm{\sigman{\hat{X}}} & = & \nm{\sqrt{\cfrac{\displaystyle\sum_{j=1}^{\nEX} \Big(f_{Xj} \lrp{t_n} - \mun{\hat{X}}\Big)^2}{\nEX}}}\label{eq:Sim_Metrics_std_s}
\end{eqnarray}

§.§.§ Drift, Bounded, and Biased Estimations

When evaluating the capability of a navigation system to estimate or track a given variable, it is necessary to determime whether the estimation experiences a drift, this is, the estimation error grows with time, or not, in which case the estimation error never exceeds a certain limit and is hence bounded. If drift is not present, the individual executions or runs can be ordered in accordance with their trajectory metrics, and the system estimation capabilities evaluated by means of the aggregated metrics. If drift can not be avoided, then the final state value constitutes the best indicator for each simulation run, while the system performance needs to be evaluated by means of the aggregated final state metrics.

The bounded estimation of a given variable is unbiased or zero mean if the aggregated metrics for the means (X̂, X̂, X̂) are significantly smaller than both the standard deviations (X̂, X̂, X̂) and the maximums (X̂, X̂, X̂), and their ratios are more pronounced as the number of executions  grows. If this is not the case, the estimation is known as biased. When drift is present, an unbiased or zero mean estimation is that in which the aggregated final state mean X̂ is significantly smaller than both the standard deviation X̂ and the maximum X̂.

every axis legend/.append style=

cycle list=blue,no markers,very thick,
red,no markers,very thick,
blue,dashed,no markers,thin,blue,dashed,no markers,thin,
red,dashed,no markers,thin,red,dashed,no markers,thin,
xmin=0, xmax=1000, xtick=0,100,...,1000, xticklabels=,,,
ymin=-6, ymax=6, ytick=-6,0,6, yticklabels=,0,,

axis lines=left,
axis line style=-stealth,
legend entries=X̂ ±X̂  unbiased (zero mean), bounded,
X̂ ±X̂  biased, bounded,
legend columns=2,
legend style=font=,
legend cell align=left,
[] (0.0,0.0) – (1000.0,0.0);
figs/metric_type.txttable [header=false, x index=0,y index=1]  ;
table [header=false, x index=0,y index=4]  ;
table [header=false, x index=0,y index=2]  ;
table [header=false, x index=0,y index=3]  ;
table [header=false, x index=0,y index=5]  ;
table [header=false, x index=0,y index=6]  ;

every axis legend/.append style=

cycle list=violet,no markers,very thick,
orange,no markers,very thick,
violet,dashed,no markers,thin,violet,dashed,no markers,thin,
orange,dashed,no markers,thin,orange,dashed,no markers,thin,
xmin=0, xmax=1000, xtick=0,100,...,1000, xticklabels=,,,
ymin=-6, ymax=6, ytick=-6,0,6, yticklabels=,0,,

axis lines=left,
axis line style=-stealth,
legend entries=X̂ ±X̂  unbiased (zero mean), drift,
X̂ ±X̂  biased, drift,
legend columns=2,
legend style=font=,
legend cell align=left,
[] (0.0,0.0) – (1000.0,0.0);
figs/metric_type.txttable [header=false, x index=0,y index=7]  ;
table [header=false, x index=0,y index=10] ;
table [header=false, x index=0,y index=8]  ;
table [header=false, x index=0,y index=9]  ;
table [header=false, x index=0,y index=11] ;
table [header=false, x index=0,y index=12] ;

every axis legend/.append style=

cycle list=magenta,no markers,very thick,
green,no markers,very thick,
magenta,dashed,no markers,thin,magenta,dashed,no markers,thin,
green,dashed,no markers,thin,green,dashed,no markers,thin,
xmin=0, xmax=1000, xtick=0,100,...,1000, xticklabels=,,,
xlabel=t s,
ymin=0, ymax=6, ytick=0,3,6, yticklabels=0,,,

axis lines=left,
axis line style=-stealth,
legend entries=|X̂ | ±|X̂ |  bounded,
|X̂ | ±|X̂ |  drift,
legend columns=2,
legend style=font=,
legend cell align=left,
figs/metric_type.txttable [header=false, x index=0,y index=13]  ;
table [header=false, x index=0,y index=16]  ;
table [header=false, x index=0,y index=14]  ;
table [header=false, x index=0,y index=15]  ;
table [header=false, x index=0,y index=17]  ;
table [header=false, x index=0,y index=18]  ;

Examples of time aggregated metrics results

If a variable corresponds to a scalar or vector component, this is, it can hold positive and negative values[Each of the body attitude Euler angles, and the geometric altitude, for example.], its estimation can be either biased or unbiased, and can also drift or be bounded. A bias is generally indicative of a problem or lack of balance in either the estimation algorithm or the input data itself, while drift may caused by defects in the algorithm or lack of observability. Examples of the four possibilities are shown in the first two plots of figure <ref>.

On the other hand, if the variable corresponds with an absolute value < |·|> or the norm of a vector <  · >, this is, it can only hold positive values[The attitude rotation vector and the horizontal position, for example.], its estimation is always biased (unless it is perfect and the error is zero) and the main consideration is whether drift is present or the estimation is bounded. The third plot within figure <ref> shows these two possibilities.

The same classification can be applied to the individual estimations for a given run instead of the time aggregated metrics. It is common and often unavoidable that the individual estimations are biased or exhibit drift while these dissapear when aggregated. Although not ideal, this generally indicates the the source of the bias or drift is random and is external to the algorithm itself, originating within the observations or input data.

§ CONCLUSIONS

This article describes an stochastic high fidelity simulation of the flight of a fixed wing low SWaP autonomous UAV in turbulent and varying weather, for which the open-source  implementation is available in [1]. The simulation is intended to test the behavior of different inertial, visual, and visual-inertial navigation systems in GNSS-Denied conditions, in which current systems experience an unbounded position drift or position error growth with time. The article describes the different simulation modules, defines two scenarios representative of GNSS-Denied conditions, and presents the metrics employed for the evaluation of the navigation performance. Future articles will make use of this simulation infrastructure to test novel GNSS-Denied navigation algorithms proposed by the author.

[1]
E. Gallo, “High Fidelity Flight Simulation for an Autonomous Low SWaP Fixed
Wing UAV in GNSS-Denied Conditions.”
<https://github.com/edugallogithub/gnssdenied_flight_simulation>, 2020.
C++ open source code.

[2]
J. A. Farrell, “Aided Navigation, GPS with High Rate Sensors,” McGraw-Hill, Electronic Engineering Series, 2008.

[3]
P. D. Groves, “Principles of GNSS, Inertial, and Multisensor Integrated
Navigation Systems,” Artech House, GNSS Technology and Application
Series, 2008.

[4]
A. B. Chatfield, “Fundamentals of High Accuracy Inertial Navigation,” American Institute of Aeronautics and Astronautics, Progress in Astronautics
and Aeronautics, vol 174, 1997.

[5]
C. Tippitt, A. Schultz, and W. Procino, “Vehicle Navigation: Autonomy Through
GPS-Enabled and GPS-Denied Environments,” state of the art report
DSIAC-2020-1328, Defense Systems Information Analysis Center, December 2020.

[6]
D. Scaramuzza and F. Fraundorfer, “Visual Odometry Part 1: The First 30 Years
and Fundamentals,” IEEE Robotics & Automation Magazine, 2011.

[7]
F. Fraundorfer and D. Scaramuzza, “Visual Odometry Part 2: Matching,
Robustness, Optimization, and Applications,” IEEE Robotics &
Automation Magazine, 2012.

[8]
“Department of Defense World Geodetic System 1984,” tech. rep., National
Imagery and Mapping Agency, January 2000.
NIMA TR8350.2 Amendment 1, edition.

[9]
“The Development of the Joint NASA GSFC and NIMA Geopotential Model EGM96,”
tech. rep., National Aeronautics and Space Administration, July 1998.

[10]
“Smithsonian Meteorological Tables,” tech. rep., 1963.
R.J. List, edition.

[11]
“Manual of the ICAO International Standard Atmosphere,” tech. rep.,
International Civil Aviation Organization, 2000.
ICAO DOC-7488/3, edition.

[12]
“Flying qualities of piloted aircraft,” tech. rep., Department of Defense,
November 1980.
US Military Handbook MIL-F-8785C.

[13]
E. Gallo, “Quasi Static Atmospheric Model for Aircraft Trajectory Prediction
and Flight Simulation,” 2021.

[14]
“The US/UK World Magnetic Model for 2015-2020,” tech. rep., National
Geophysical Data Center, NOAA, 2015.

[15]
M. M. Madden, “Gravity Modeling for Variable Fidelity Environments,” AIAA Modeling and Simulation Technologies Conference and Exhibit, 2006.

[16]
M. Drela and H. Youngren, “Athena vortex lattice.”

[17]
M. E. Eshelby, Aircraft Performance, Theory and Practice.
Arnold, 2000.

[18]
B. Etkin, Dynamics of Atmospheric Flight.
John Wiley & Sons, 1972.

[19]
J. B. Heywood, Internal Combustion Engine Fundamentals.
Mechanical Engineering Series, McGraw-Hill, 1988.

[20]
P. Hill and C. Perterson, Mechanics and Thermodynamics of Propulsion.
Addison Wesley Longman,  ed., 1992.

[21]
Joan Sola, Quaternion Kinematics for the Error-State Kalman Filter,
October 2017.

[22]
E. Gallo, “Customizable Stochastic High Fidelity Model of the Sensors and
Camera onboard a Low SWaP Fixed Wing Autonomous Aircraft,” 2021.

[23]
J. L. Crassidis, “Sigma-Point Kalman Filtering for Integrated GPS and
Inertial Navigation,” IEEE Transactions on Aerospace and Electronic
Systems, vol. 42, no. 2, pp. 750–756, 2006.

[24]
O. J. Woodman, “An Introduction to Inertial Navigation,” Tech. Rep. 696,
University of Cambridge Computer Laboratory, August 2007.

[25]
M. Kayton and W. R. Fried, Avionics Navigation Systems.
John Wiley & Sons,  ed., 1997.

[26]
Y. Ma, S. Soatto, J. Kosecka, and S. S. Sastry, An Invitation to 3-D
Vision, From Images to Geometric Models.
Imaging, Vision, and Graphics, Springer, November 2001.

[27]
Open Scene Graph.

[28]

[29]
T. Hagglund, “Signal Filtering in PID Control,” IFAC Proceedings
Volumes, vol. 45, no. 3, pp. 1–10, 2012.