Abstract
Safety certification of GNSS-aided inertial navigation systems (INS) in civil aircraft requires thorough testing to ensure proper operation, even in worst-case conditions. One error that must be considered is that of gravity compensation on accelerometer measurements. Prior to the work described in this paper, no stochastic models existed with the Gaussian bounding of the tails required to ensure integrity performance. This paper describes a method to determine efficient stochastic models of the error of current high-order gravity models such as EGM2008. The stochastic and high-order models are combined to achieve a high-fidelity model suitable for use in testing systems designed for low-approach operations such as RNP-AR. This paper also describes a method to determine efficient stochastic models for low-order gravity models such as the WGS-84 ellipsoidal model. Such models may be used in testing systems designed for operations with less stringent lateral requirements.
1 INTRODUCTION
Safety certification of GNSS-aided inertial navigation systems (INS) in civil aircraft requires thorough testing to ensure proper operation, even in worst-case conditions. The heart of an INS is the orthogonal triad of specific-force sensors called accelerometers that measure the vector sum of Newtonian acceleration and its reaction to gravity. The INS computes estimates of gravity in order to isolate the Newtonian acceleration from the specific-force measurements, thereby enabling the subsequent determination of velocity and position. A wide range of gravity models and databases exist that vary in fidelity. Even high-order spherical harmonics-based gravity models contain errors that result in acceleration, velocity, and positioning errors. Previous research by the authors has shown that the mis-modeling of gravity in a GNSS-aided navigation-grade INS can result in horizontal positioning errors of nearly 200 meters after coasting (i.e., continued operation after loss of GNSS) for only 10 minutes (Needham & Braasch, 2017).
Standards-making bodies, such as RTCA in the United States, develop and publish, in cooperation with certification authorities, requirements and testing guidance for manufacturers of aviation equipment. Current guidance for so-called gravity compensation testing in RTCA documents DO-229 (RTCA, 2020a) and DO-316 (RTCA, 2009) is limited in applicability to enroute phases of flight in which high-fidelity gravity modeling is not required. Thus, the current testing methods are not appropriate for low-approach operations such as Required Navigation Performance–Authorization Required (RNP-AR). Furthermore, although general guidance has been given in the aforementioned documents, those desiring to simulate gravity for less stringent flight operations are still faced with the task of deriving an acceptable model. Such models must mimic the behavior of the actual gravity field even in low-probability cases. Although the Earth’s gravity field is obviously deterministic, we hypothesize that the variations of the true gravity vector from the gravity models can be viewed/modeled as stochastic processes.
As will be shown in the remainder of this paper, the aforementioned hypothesis is reasonable and, as a result, computationally efficient stochastic models can be derived for use in certification simulations. The first contribution is a method to derive and apply efficient stochastic models of the error of existing high-order gravity models such as EGM2008 (Pavlis et al., 2012). As mentioned earlier, even high-order gravity models contain errors and, for the purposes of safety certification testing, these errors cannot be ignored. In Monte Carlo simulations performed during certification testing, the derived models can be used along with the existing high-order models to simulate nominal and worst-case gravity environments.
The second contribution is a method to derive stochastic models of the error of commonly used low-order gravity models such as the WGS-84 ellipsoidal gravity model (National Geospatial-Intelligence Agency, 2014). For equipment certification involving phases of flight with requirements less stringent than RNP-AR, this lower fidelity simulation method has the advantage of being significantly less computationally intensive than the first method. The two methodologies presented herein may thus be utilized to test various systems designed to operate in various phases of flight.
The authors have served on RTCA Special Committee 159 (SC-159) Working Group 2C (WG2C) whose work was recently published as DO-384: Minimum Operational Performance Standards (MOPS) for GNSS Aided Inertial Systems (RTCA, 2020b). The authors of the present article were the primary developers of the gravity modeling appendix in that document. Much of the work presented in this paper was performed in support of that effort.
The remainder of the paper is organized as follows: Section 2 describes gravity compensation and the impact of mis-modeling. Section 3 provides a graphical explanation of the little-understood phenomenon whereby the observed spatial spectra of gravity deflections differ significantly in along-track versus cross-track directions. Section 4 presents a high-fidelity gravity simulation architecture suitable to support certification testing in phases of flight such as RNP-AR. Section 5 describes the technical approach used to derive the stochastic model portions of the high-fidelity simulation and an example model derived with this approach is given in Section 6. Section 7 presents a low-fidelity gravity simulation architecture suitable to support certification in less-stringent phases of flight such as RNP X where X >= 0.3. Also described is the technical approach and an example. Finally, conclusions and recommendations are presented in Section 8.
2 GRAVITY COMPENSATION AND MIS-MODELING
The well-known accelerometer-sensing equation illustrates the relationship between sensed/measured specific-force Newtonian acceleration and its reaction to gravity:
1
in which all three terms are vectors. Note that the right-hand side of the equation is the sum of Newtonian acceleration and the reaction to gravity. It is not acceleration minus gravity.
Nevertheless, gravity compensation is a simple rearrangement of this equation to isolate the desired Newtonian acceleration:
2
Thus, gravity compensation requires the determination of a gravity vector at the point in space where the INS is operating. Gravity is determined either through a software algorithm/model or a database. No gravity model or database is perfect, however, and mis-modeling results in an error in the determined Newtonian acceleration. The gravity modeling error is typically described as the difference between the true and modeled gravity vectors. The difference between the vectors is characterized by the anomaly (difference in vector magnitude) and the deflection of the vertical or DOV (difference in vector angle).
DOV is commonly decomposed into angular components in the local north-south and east-west planes. It is well known that the gravity anomaly affects vertical positioning performance but has negligible impact on horizontal positioning. Conversely, as will be described shortly, the impact of the DOV is most pronounced in the horizontal. Since barometric, not inertial, altitude governs the key points in instrument approach operations (e.g., decision height; ICAO, 2009), the main concern in gravity mis-modeling for civil aviation is DOV rather than the anomaly.
As shown in Figure 1, DOV results in a horizontal component of the true gravity vector. An error in the modeling/compensation of DOV results in a net horizontal acceleration error of approximately 5 micro-g per arc-second of uncompensated deflection. Unsurprisingly, large deflections occur in the vicinity of large mountains, ocean trenches, and continental shelves. Two examples of these were described by the authors previously (Needham & Braasch, 2020) and are summarized again here.
Illustration showing the DOV as the angle of the effective gravity vector relative to the ellipsoidal normal
Figure 2 depicts a simulation of the required navigation performance (RNP) 0.3 approach into Runway 02 at the Kahului Airport on the island of Maui. Figure 3 illustrates the north-south component of the DOV along this approach. Note that over the span of approximately 60 nautical miles, the DOV changes from a maximum value of 55 arc-seconds to a minimum value of nearly –30 arc-seconds. The deflection, thus, changes by nearly 85 arc-seconds, or equivalently, more than 400 micro-gs.
Trajectory of an RNP 0.3 approach into Runway 02 of the Kahului Airport
North-south DOV for an approach into Runway 02 of Kahului Airport
Figure 4 depicts a simulation of the RNP approach into Runway 33 at Ted Stevens Anchorage International Airport in Alaska. The circling maneuver near the end of the approach is an actual part of the published procedure and is required in order to lose altitude after passing over the mountains.
Simulated approach onto Runway 33 at Ted Stevens Anchorage International Airport
Figure 5 illustrates the north-south component of the DOV along this approach. The large excursion that has a negative peak approximately 190 nautical miles from the airport is due to the effects of the continental shelf, whereas the high spatial frequencies observed closer to the airport are due to the Kenai Mountains to the South and the Chugach mountains to the East. As will be described in more detail later, high-order gravity models such as EGM2008 capture low spatial frequencies, but do not fully capture the high spatial frequencies. Low-order models such as WGS-84 ellipsoidal gravity do not capture either high or low spatial frequencies.
North-south DOV for an approach onto Runway 33 of Ted Stevens Anchorage International Airport
3 EFFECTS OF A STATIONARY MASS ON A MOVING OBJECT
As described earlier, DOV is typically decomposed into two orthogonal components which typically are chosen in the north-south and east-west directions. However, it has been known for over 50 years that the observed spatial spectra of DOV components are dependent upon direction of travel (Harriman & Harrison, 1986; Shaw et al., 1969).
Specifically, the observed spatial power spectrum of the DOV along-track component is significantly different than that of the cross-track component. This phenomenon must be taken into account when stochastic models of DOV are being developed for moving vehicles such as aircraft. Since this phenomenon is not well understood within the civil aviation community, a brief graphical explanation is presented here.
First, recall the vector form of the inverse-square law of universal gravitation:
3
where G is the gravitational constant, m1 and m2 are the masses of the two objects, r21 is the distance vector between them, and 
is the unit vector from Object 1 to Object 2.
Now, consider an object (Object 1) moving in a straight line, passing a stationary object (Object 2) offset from the moving object’s line of motion as shown in Figure 6. In this diagram, the positive y-axis represents the positive along-track direction, and the positive x-axis represents the cross-track direction. At times t−∞ and t∞, the moving object is sufficiently far from Object 2 so as to consider the line-of-sight between them to be parallel with Object 1’s direction of motion. At times t1 and t3, the gravitational force on Object 1 has components in both the along-track and cross-track directions. At time t2, Object 2 is perpendicular to Object 1’s line of motion and, thus, no along-track gravitational force exists.
Visualization of gravitational forces acting on Object 1 as it moves on a path offset from Object 2
However, at this same time, the cross-track component reaches a maximum. As the object continues in the positive y direction, the along-track component grows to a maximum (as the line-of-sight projects more and more onto the along-track direction) but eventually decays due to the 1/r2 factor. The along-track and cross-track components of the force, plotted as a function of the y position coordinate, are shown in Figure 7.
Normalized along-track and cross-track components of the gravitational force experienced as Object 1 passes Object 2 that is offset from the trajectory
Figure 7 clearly demonstrates the fact that the along-track and cross-track components have different properties. First, the along-track component is zero when the cross-track component has its maximum value. Furthermore, the cross-track component has a non-zero mean (as opposed to the along-track component) and, in fact, is non-negative (given the chosen coordinate axes). It is, thus, not surprising that the observed spatial spectra of these components are significantly different. This is demonstrated in the power spectra shown in Figure 8. As expected with a zero-mean signal, the along-track component had a null at zero frequency, whereas the cross-track component had a maximum at zero frequency. The results are consistent with those of a similar analysis seen in Harriman and Harrison (1986).
Normalized spatial power spectra demonstrating differences between the along-track and cross-track components of the gravitational forces
It is important to emphasize that the differences in the spectra are not inherent to the gravity field, itself, but rather how the physical forces project onto the along-track and cross-track directions of a moving object. The impact of this phenomena in the present context is that the gravity mis-modeling error components depend on the direction of travel of the moving vehicle. Thus, as will be described in a later section, separate error simulation models must be developed for each component.
4 HIGH-FIDELITY GRAVITY SIMULATION ARCHITECTURE
The development of efficient, yet robust, error models is needed to enable the large-scale Monte Carlo simulations typically used in the certification testing of safety-critical navigation systems. The gravity error models must generate simulated data that is statistically representative of the errors (also referred to as residuals) that remain after gravity compensation.
Figure 9 illustrates a high-fidelity gravity model error simulation architecture. The input to the gravity model error simulator is the vehicle trajectory being simulated and the output is the gravity DOV residuals (comprised of two orthogonal components). The residuals are comprised of deterministic and stochastic components.
Block diagram demonstrating the generation of the total simulated gravity model error
A high-accuracy reference model (e.g., EGM2008) was used to form the deterministic portion by differencing it with the estimate formed by the system-under-test (navigator gravity model). The stochastic component simulates the error of the reference model with autoregressive moving average (ARMA) models that are described in the next section.
This architecture has distinct advantages over simpler alternatives. Previous research has shown that high-fidelity stochastic modeling of the total DOV cannot be done efficiently due to the non-stationary nature of the data (Needham & Braasch, 2020). As described in a later section, low-fidelity stochastic models can be derived, but they are not appropriate for low-approach operations with the most stringent lateral requirements. Alternately, generating/simulating the DOV solely with the high-accuracy reference model is not acceptable, either, due to known errors in these reference models. The approach presented here solves both of these problems by utilizing the high-accuracy reference model to account for the known deterministic portion of the DOV, with the addition of a stochastic model to account for the errors in the reference model.
5 HIGH-FIDELITY STOCHASTIC MODEL DEVELOPMENT APPROACH
The first step in developing the stochastic model in question was to determine the errors in the high-accuracy reference model. Reference models such as EGM2008 are valid over the entire surface of the Earth and above the surface, up to altitudes far above the operations of any civilian aircraft. Currently, there is no source of truth value that can be used to evaluate the accuracy of these reference models over their entire volume of applicability.
However, higher accuracy regional databases do exist that are valid only at the surface of the Earth. For example, the DEFLEC regional data (NGS, 2021a) from the National Geodetic Survey (NGS) are valid over much of the western half of the northern hemisphere. The DEFLEC data is more accurate than the EGM2008 data, as it contains additional gravity measurements such as those collected by the GRAV-D mission (NGS, 2021b). In addition, terrain databases are used to estimate higher frequency gravity terms due to irregular topography (Hardy, 2020).
Since the DOV generally decreases in magnitude with altitude (i.e., the Earth becomes more like a point-mass), the reference model errors (computed by differencing the reference model with the higher-accuracy database) determined near the surface are a conservative bound on the errors that would be experienced at a given operational altitude. In addition, with an increase in altitude, the higher frequency components of the DOV are attenuated more than the lower frequency components (Vanderwerf, 1996), thus, further increasing the conservatism.
Previous work by the authors (an example plot is reproduced here in Figure 10) showed that EGM2008 performs well in capturing lower spatial frequencies in the DOV, but exhibits non-trivial residuals in the higher spatial frequency components (Needham & Braasch, 2020). The reference model residuals, such as those depicted in Figure 10 for EGM2008, are used in the simulation architecture described herein as the basis for stochastic model determination.
DOV profiles as determined from EGM2008 and DEFLEC data sets for a simulated approach into Ted Stevens International Airport in Anchorage, Alaska
Simulated trajectories can be generated over the regional area and gravity model errors can then be formed by differencing the regional database and reference gravity model estimates. Thresholds may be set to extract only the data sets containing significant differences. Following this, segments of data are chosen with approximately constant mean and variance. Ideally, segments would be chosen that are statistically stationary, but this is rarely possible in practice and, thus, we chose segments that were second-order weakly stationary.
The formation of these residual data sets is the first step of the model development process as outlined in Figure 11. The residuals are then processed using time series system identification to generate model parameters (with the obvious under-standing that these data sets are actually spatial series rather than time series; Cryer & Chan, 2008; Ljung, 1999). System identification is used to fit each of the individual runs in the data set with candidate models of various orders.
Block diagram outlining the process for using system identification to analyze measured data, identify, and validate representative ARMA models
A generalized model that can be used for system identification is the autoregressive integrated moving average (ARIMA) model. The ARIMA model allows for time-differencing non-stationary data to create stationary time series. In contrast, the ARMA model is sufficient when the original data set is stationary.
The ARMA(p,q) model is a combination of a AR model of order p and a moving average (MA) model of order q, and can be represented mathematically as Equation (4), where the AR coefficients are represented by φ and the MA coefficients are represented by θ. The input to the model is a white noise term with unity variance zk, whose amplitude is scaled by the white noise standard deviation σwn:
4
After creating the measured data set, the next step is selecting the minimum order model that is sufficient to represent the data. Time series system identification is utilized to fit ARMA models for various orders up to ARMA(15,15) to each of the runs. The outputs of system identification include model order selection criteria and model coefficients for each of the time series analyzed.
A variety of model order selection criteria are described in the literature. These include the Akaike Information Criterion (AIC; Akaike, 1970), the Final Prediction Error (FPE; Akaike, 1974), and the Bayesian Information Criteria (BIC; Schwarz, 1978). An advantage of the BIC is that it imposes a greater penalty for increasing model complexity relative to the other criteria. Since an aim of this work is to identify the most efficient (i.e., simplest) models, the BIC was selected for the present work.
The coefficients resulting from system identification represent the structure of the model that best fits a particular data sample for a given (p,q) order combination. Once a model order is selected via the model order selection criteria, a deeper analysis of the model coefficients can be conducted for the selected order. The coefficients can be visualized by viewing the corresponding roots in a pole-zero plot. As mentioned previously, gravity cannot be treated as a single stationary stochastic process. Individual stationary data segments can be processed by system identification, but each segment will produce slightly different roots. To account for the varying roots, a region of interest can be defined around the identified roots.
Lastly, any model must be validated against the original measured data. One validation check is to compare the power spectral densities of the measured data to that of the candidate models. The time series data generated by the model should have similar spectral properties to that of the data. To demonstrate this process, two examples are now given.
The first example is geared towards a high-end user (i.e., DO-384 Category A) such as a commercial aircraft that can coast for extended amounts of time, whereas the second example is intended for applications such as GNSS-enhanced attitude and heading reference system or AHRS (i.e., DO-384 Category C) that utilize lower cost inertial sensors that can coast only for brief periods of time.
6 EXAMPLE HIGH-FIDELITY STOCHASTIC MODEL
An example is now shown to demonstrate the development of cross-track and along-track ARMA models representative of high-order reference model errors. The reference model and regional truth data are EGM2008 and xDEFLEC19 (NGS, 2021c), respectively, and the measured data to be analyzed is created by differencing the two data sets.
In this example, 24 cross-track and 28 along-track DOV time series data were analyzed, each 100 nautical miles long and running north-to-south with a sampling interval of one nautical mile (NM; note that data obtained from east-to-west trajectories yield similar results). The 100-NM distance was chosen to be representative of what an aircraft could potentially coast through during a GNSS outage on approach to an airport.
These series were down-selected from a larger data set of 6,000 random trajectories by selecting those that contained at least one data point larger than 20 arc-seconds. This threshold of 20 arc-seconds was chosen to identify the worst-case scenarios that are of interest in this work. As a result, the benign conditions that made up the majority of the data set were not considered.
The data set was then screened so that only second-order weakly stationary runs were used for system identification. This screening helped to ensure a consistent set of data from which to perform system identification. Example cross-track and along-track reference model errors are shown in Figure 12 and Figure 13, respectively.
Example differences between EGM2008 and xDEFLEC19 cross-track components of DOV for select trajectories
Example differences between EGM2008 and XDEFLEC19 along-track components of DOV for select trajectories
6.1 Example Cross-Track Model Structure Development
Once the aforementioned data set was created, system identification was performed on each of the individual data sets for ARMA models up to ARMA(15,15). The first items of interest were the model order selection criteria values. The results of each of the individual runs varied slightly in terms of BIC and model coefficients.
To provide a measure of overall performance, the BICs for the select runs were averaged for a particular order combination such that a single value was calculated for each (p,q) combination as shown in Figure 14. An inspection of Figure 14 indicates that the minimum BIC value occurred with ARMA(2,0; i.e., AR[2]), thus making it the primary candidate model.
Mean of the BIC for each of the selected cross-track runs
Next, the AR(2) model coefficients resulting from system identification were analyzed. The coefficients of the AR(2) model were used to calculate the corresponding roots and view them using a pole-zero plot as shown in Figure 15. The clustering apparent from a visual inspection of the figure demonstrates an under-lying stochastic structure of the residuals and adds weight to the choice of AR(2) models for the cross-track component. The combination of low BIC and consistent poles indicates that an AR(2) model is sufficient to represent the underlying structure of the cross-track data for the given application.
Pole-zero plot showing the cross-track and along-track pole results from system identification of each of the trajectories
The desire is to create a single, efficient model for each track component. A means to do so is to take the centroid of the cluster of roots and use the corresponding coefficients to create a single model. This approach allows for a single model to be used for each realization.
A more robust approach is to identify the region of roots and draw uniformly from this region for each Monte Carlo realization, thus ensuring that, over the course of a large-scale Monte Carlo simulation, the worst-case scenarios would be exercised. In the present example, the centroids were used to develop a sample model for further analysis.
Once a second-order AR model had been selected and model coefficients were identified via the centroids, the model could be mathematically represented as Equation (5), in which coefficients ϕ1 and ϕ2 control the structure of the signal and σwn is the standard deviation of the driving unity-variance white noise source zk.
5
For AR(2) models, the standard deviation of the driving white noise input (σwn) is related to the standard deviation of the AR model output (σx) as follows (Cryer & Chan, 2008):
6
The final step is validating the model against the original measured data by comparing the power spectral density (PSD) of the measured data to that of the model specified in Equation (5). A comparison of the PSDs is provided in Figure 16 and demonstrates an agreement between the measured and simulated data as the spectrum of the simulated data falls well within the 95% confidence bounds. Validating the AR(2) model completes the final step in developing the cross-track ARMA model. A similar process for the along-track data is conducted and presented next.
Comparison of the measured data PSD to the simulated data to show validation of the AR(2) model specified in Equation (5)
6.2 Example Along-Track Model Structure Development
The along-track model development follows the same process as the cross-track case. The model order selection criterion was first calculated for each (p,q) combination up to ARMA(15,15) as shown in Figure 17. For the along-track models, the ARMA(2,2) had the minimum BIC value. To investigate this model order choice, the pole-zero plot was analyzed and sample roots from the distribution were selected and used to calculate the coefficients in the model shown in Equation (7). Similar to the cross-track example, this model is validated against the measured data by comparing the PSDs of each.
7
Mean of the model order selection criteria (i.e., BIC) over a range of (p,q) combinations up to ARMA(15,15)
Pole-zero plot showing the ARMA(2,2) poles and zeros for model fits in addition to the roots selected for the ARMA(2,2) model used for validation
The measured data PSD with confidence bounds compared to the PSD produced by the ARMA(2,2) model in Equation (7)
6.3 Model Magnitude Determination
Having identified the model coefficients for both the along-track and cross-track cases, attention is now turned to the magnitude of the model output. The AR coefficients dictate the stochastic structure of the data (i.e., serial correlation), but the standard deviation of the model output might be scaled arbitrarily by the standard deviation of the white noise input as shown, for example, in Equation (5) for the AR(2) case.
The magnitude of the aforementioned model errors (i.e., the difference between the EGM2008 model and the DEFLEC data) can be assessed. Since a large portion had small magnitudes, specific focus was given to the tails. It is important to bound the tails as these are the worst-case situations that need to be protected against. In cases where the tails are not well-bounded by a known probability distribution, bounding techniques can be utilized to ensure proper bounding.
Examples of this technique are Gaussian bounding (Blanch et al., 2018) and Gaussian-Pareto bounding (Larson et al., 2019). Although alternative techniques exist, Gaussian (or Gaussian-related) bounding has gained broad acceptance in the civil aviation community for its utility in demonstrating that a given algorithm meets a given set of integrity performance requirements.
An example of Gaussian bounding is shown here. The Gaussian bounding methods and tools used are those described in Blanch et al. (2018). The tools provide a one-sided bound, thus both sides of the distribution need to be bounded separately. To be conservative, the larger of the two standard deviations was used to define the bounding two-sided Gaussian distribution.
In this example, a Gaussian distribution with a standard deviation of 16.3 arc-seconds bounds the residuals for both the along-track and cross-track components (note that the bounding Gaussian distribution is zero-mean). The complementary cumulative distribution function (CDF) is shown in Figure 20 and demonstrates the bounding of the tails of the residual distribution.
Complementary CDF showing that Gaussian distribution bounds the tails of the model errors
Once the model parameters were identified, the gravity model errors could be simulated as shown in Figure 21. Simulated reference model error/residuals were generated using the along-track and cross-track models of Equation (5) and Equation (7) and σx = 16.3 arc-seconds and is shown in Figure 22.
Block diagram demonstrating the use of ARMA models to generate simulated reference model
Example simulated reference gravity model errors using the models given in Equations (5) and (7) and a model standard deviation of σx = 16.3 arc-seconds
The standard deviations for the white noise input were 9.7 arc-seconds for the cross-track model and 7.1 arc-seconds for the along-track model. As shown in Figure 21, the final simulated gravity model residuals were obtained by summing this result with the deterministic gravity model error (such as the difference curve shown in Figure 2).
7 LOW-FIDELITY STOCHASTIC MODEL DEVELOPMENT APPROACH AND EXAMPLE
Although the high-fidelity modeling approach described provides both deterministic and stochastic characterization of the DOV, it might prove computationally burdensome for some applications. If an inertial system utilizes accelerometers with bias instabilities that are larger than the worst-case total DOV, a high-fidelity gravity model would not be needed for testing.
For example, inertial systems containing accelerometers with bias instabilities on the order of one milli-g or more contain inherent system errors that render the total DOV variations minor in consideration. It is further noted that in such systems the navigator gravity models are typically low-order models, which account for the DOV minimally. Thus, the total DOV is essentially the gravity model error.
Even though these systems do not have the stringent accuracy and integrity requirements of low-approach operations such as RNP-AR, gravity compensation errors cannot be ignored completely in safety certification testing. There is, thus, a need for a computationally efficient low-fidelity model that captures the structure and magnitude of the worst-case total DOV variations.
A key feature of the low-fidelity approach is that the actual DOV does not need to be calculated for every location point used in the Monte Carlo simulation, thus reducing the computational burden of the testing process. The remainder of this section discusses low-fidelity approaches followed by an example.
An assumption for the low-fidelity case is that high-order gravity model errors are negligible compared to accelerometer errors. The maximum DOV on Earth is on the order of 100 arc-seconds, or the equivalent of 500 micro-g’s, and the errors in high-order reference models are less than this value (by at least a factor of two). Thus, when testing systems containing accelerometers with relatively large bias instabilities (e.g., 1 milli-g or larger), the DOV from the high-order reference models can be used by itself in simulations.
The goal is to create time series data that is sufficient to represent the residuals remaining after gravity compensation so as to thoroughly stress the navigation system under testing. Although the high-order reference models are computationally burdensome, a simpler approach is to use a worldwide database, such as the EGM2008 database that can be downloaded from the National Geospatial-Intelligence Agency (2021), and to sample the DOV from randomly selected trajectories across the globe and, thus, from gravitationally diverse regions.
Since these databases are only valid at the surface of the reference ellipsoid, they cannot be used for high fidelity applications where the effects of altitude variation must be taken into account. Again, the sample data can be directly used as input to the Monte Carlo testing. In this approach, the sample DOV series is not tied to the actual positions of the Monte Carlo trajectories. Instead, a given trajectory utilized in Monte Carlo testing was subjected to a wide variety of DOV profiles. Since the DOV data is available in a database, the computational burden is small.
Another approach is to create a stochastic model or models of the entire DOV so as to not require a deterministic component (as was the case in the high-fidelity approach). The aim is to create a computationally simple ARMA model that is sufficient to represent the underlying statistics of the DOV. The trajectories in the high-fidelity example were concentrated in areas near mountain ranges so as to focus on the worst-case areas.
The lower frequency terms were ignored by the high-fidelity stochastic model as they were already captured by the deterministic part of the simulator. However, for total DOV modeling, worst-case lower frequency terms like those due to the continental shelf (e.g., left side of Figure 4) or large, more isolated mountains such as Hawaii (Figure 3) need to be included to ensure robustness in the low-fidelity approach.
A dilemma that arises with the low-fidelity approach is that a single model cannot reproduce residuals representative of both the low and high frequencies mentioned above. A robust low-fidelity simulator architecture needs to account for these different worst-case scenarios in some manner. One method is to determine separate models for each case and implement some mechanism to switch between the models throughout the simulation as demonstrated in the block diagram of Figure 23.
Block diagram demonstrating an example of low-fidelity architecture
The choice of switching logic is dependent upon its application. Some users may want to generate an initial quiet segment with little or no DOV, followed by a low-frequency DOV segment to simulate a trajectory that traverses a single large geographic feature. Other users may wish to generate high-frequency DOV segments only to simulate a trajectory over a mountainous region.
An example is now given that shows the development of two models for each of the track components. The first models were based on the same trajectories that were used in the high-fidelity example to capture the higher-frequency errors such as those in mountain ranges. The second set of models was based on analyzing trajectories passing near Hawaii and the Puerto Rico trench to represent areas with geographies that produce large, low-frequency gravity terms. The DEFLEC data was chosen as the source of truth data for DOV.
7.1 Geographic Area #1: Large Mountain Ranges
First, the DOV near large mountain ranges was considered as these areas produced large DOV deflections that could change rapidly as shown in Figure 24 and Figure 25 for the cross-track and along-track directions, respectively.
Sample cross-track total DOV for a select trajectory
Sample along-track total DOV for a select trajectory
The same trajectories identified in the high-fidelity example were examined. System identification was performed in the same way as described earlier and it was found that an ARMA(1,1) yielded the minimum model order selection criteria for both track components (see Figure 26).
Mean BIC values for each ARMA(p,q) combination for the cross-track (left) and along-track (right)
The roots of the ARMA(1,1) model fits shown in Figure 27 were analyzed to gain an understanding of the distribution. The poles and zeros in Figure 27 were not clustered especially tightly and this was due to the diversity of the actual terrain. As mentioned earlier, simulation fidelity can be increased by randomly sampling poles from the range of identified locations to generate simulation models.
Pole-zero plots showing roots of the ARMA(1,1) fits to the low-order cross-track and along-track
To illustrate validation, the centroids were calculated to form an example model and the corresponding ARMA(1,1) coefficients were found as shown in Equation (8) and Equation (9):
8
9
Simulated time-series data was generated using these models and the simulated PSD was compared to that of the original measured data as shown in Figure 28. A visual inspection of the PSDs demonstrated the similarities in spectrum and served as a validation step.
PSDs of the measured data compared to the spectrum of the simulated data generated by Equation X and Equation Y
DOV for a trajectory over Puerto Rico showing the large low-frequency effects and gravitational effects of the trench
For an ARMA(1,1) model, the standard deviation of the white noise input (σwn) is related to the standard deviation of the model output (σx) as follows (Brockwell & Davis, 1996):
10
7.2 Geographic Area #2: Isolated Large Geographic Features
The second area to be studied had large but slowly changing DOV values such as those demonstrated in Figure 3 for an approach into Hawaii or a flight over the Puerto Rico trench, where an 82 arc-second DOV change occurred over an 80-NM span (Needham & Braasch, 2017).
To better capture the longer wavelengths, segments 300-NM long were studied. The cross-track DOV profiles were processed via system identification and resulted in the BIC as shown in Figure 30. An inspection of the figure demonstrated that an increase in MA order provided minimal improvement, while an increase along the AR order axis showed a significant drop at first, then little improvement was gained after order three.
Mean BIC for the cross-track low-fidelity residuals from areas near Hawaii and Puerto Rico
The lowest BIC occurred with a model order combination of ARMA(3,3), but many of the combinations nearby had similar values. Given the low-fidelity nature of the model being developed, the slightly simpler AR(2) and AR(3) models were first selected and the calculated roots were examined further.
The AR(2) roots shown in Figure 31 were tightly clustered, whereas the AR(3) roots in Figure 32 did not exhibit the same tight clustering. Thus, the AR(2) model was selected for the cross-track component.
Pole-zero plot demonstrating the results of the system identification on cross-track DOV of trajectories near Hawaii and the Puerto Rico trench using AR(2)
Pole-zero plot demonstrating the results of the system identification on DOV of trajectories near Hawaii and the Puerto Rico trench using AR(3)
Low-fidelity Geographic Area #2 along-track
Next the along-track direction of the second geographic area was considered. The model order selection criterion yielded similar results as the cross-track case and demonstrated similar properties along the MA and AR order axes. The lowest BIC occurred with ARMA(4,2), however the values of nearby combinations were similar. Similar to the cross-track example, the AR(2) and AR(3) models were selected as candidates for further evaluation due to their relatively low BIC values.
The AR(2) roots in Figure 34 clustered more tightly than the AR(3) roots shown in Figure 35. Given the relatively low value of BIC and the tight clusters, the AR(2) model was deemed sufficient to represent the structure of the along-track model.
Pole-zero plot demonstrating the results of the system identification on along-track DOV of trajectories near Hawaii and the Puerto Rico trench using AR(2)
Low-fidelity Geographic Area #2 along-track
To validate the selected AR(2) models, the centroids of the roots were used to produce coefficients of the example models of Equation (11) and Equation (12).
Example cross-track AR(2) model:
11Example along-track AR(2) model:
12
The models were used to generate data so the PSDs of the measured and simulated data could be compared. The PSDs are shown in Figure 36 for the cross-track direction and Figure 37 for the along-track direction and demonstrate agreement.
The measured data PSD with confidence bounds compared to the PSD produced by AR(2) cross-track model in Equation (10)
The measured data PSD with confidence bounds compared to the PSD produced by AR(2) along-track model in Equation (11)
7.3 Low Fidelity Magnitude Determination
Again, as was done with the high-fidelity models, the amplitude of the models must be addressed. In order to fully characterize the magnitude variations, all of the DOV data were utilized, not just the segments used for coefficient determination. The data distribution and bounding Gaussian distribution are shown in Figure 38 using the complementary CDF.
Complementary CDF showing that Gaussian distribution bounds the tails of the low-fidelity model errors
The standard deviation of the bounding Gaussian was 25.23 arc-seconds, which was used to set the input white noise standard deviation accordingly. Example data generated with the models of Equation (10) and Equation (11) are plotted in Figure 39.
Example cross-track and along-track data generated using the model as specified in Equation (10) and Equation (11)
8 CONCLUSION AND RECOMMENDATIONS
This article presented a method to determine efficient stochastic models that could be used to generate time/spatial series samples representative of gravity model errors. The models provided a means for manufacturers to thoroughly test gravity compensation schemes during certification testing.
This article also described the physical phenomena that induces differences in the statistical properties of the along-track and cross-track errors. Separate ARMA models were developed for both the along-track and cross-track directions. A test framework was outlined to demonstrate how the ARMA models could be incorporated for testing of a navigation system, with examples given for high-fidelity and low-fidelity applications.
HOW TO CITE THIS ARTICLE
Needham, T., & Braasch, M. (2022). Gravity modeling in GNSS-aided inertial navigation system safety certification. NAVIGATION, 69(2). https://doi.org/10.33012/navi.520
ACKNOWLEDGEMENTS
The authors would like to thank the members of the RTCA SC-159/WG-2C working group for their helpful comments and suggestions. The authors would like to thank Honeywell for funding portions of this work. The authors also extend their gratitude to the reviewers of this paper for their many helpful comments and constructive criticisms.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
- Article
- Abstract
- 1 INTRODUCTION
- 2 GRAVITY COMPENSATION AND MIS-MODELING
- 3 EFFECTS OF A STATIONARY MASS ON A MOVING OBJECT
- 4 HIGH-FIDELITY GRAVITY SIMULATION ARCHITECTURE
- 5 HIGH-FIDELITY STOCHASTIC MODEL DEVELOPMENT APPROACH
- 6 EXAMPLE HIGH-FIDELITY STOCHASTIC MODEL
- 7 LOW-FIDELITY STOCHASTIC MODEL DEVELOPMENT APPROACH AND EXAMPLE
- 8 CONCLUSION AND RECOMMENDATIONS
- HOW TO CITE THIS ARTICLE
- ACKNOWLEDGEMENTS
- REFERENCES
- Figures & Data
- Supplemental
- References
- Info & Metrics
Related Articles
Cited By...
- No citing articles found.