## Abstract

The performance of a proposed high-frequency (HF) navigation concept is analyzed using simulated data. The method relies on pseudorange and beat carrier-phase measurements of signals that propagate in the ionosphere along curved trajectories, where signals are refracted back downwards from the ionosphere. It has been demonstrated that the location of a receiver can be determined if several signals, broadcast from beacons at different locations, are received and processed at a user receiver. A challenge of determining exact signal paths is the uncertainty in the ionosphere’s electron density distribution. This is addressed by a batch filter that simultaneously estimates the receiver position along with corrections to a parametric model of the ionosphere. A previous paper developed the theory and batch filter for this concept. The present study examines its potential performance. Total horizontal position errors on the order of tens to hundreds of meters are achieved, depending on the case’s characteristics.

## 1 INTRODUCTION

Alternatives to satellite-based navigation have been widely explored in the past two decades. This search is motivated by concerns about the rise of threats from jamming and spoofing. High-frequency (HF) navigation is a concept that has been proposed by Baumgarten and Psiaki (2017) and Baumgarten et al. (2021) as an alternative radio-based positioning method. It relies on the determination of the exact paths traveled by electromagnetic signals with carrier frequencies in the range 2-10 MHz. These signals, which are broadcast from static ground stations, are refracted from the ionosphere to arrive at a receiver whose location is to be determined. The present study explores the potential performance of the proposed method.

Radio-navigation methods are generally based on the fusion of multiple measurements of radio-frequency (RF) signals and require mathematical models for the measured observables. In the context of HF navigation, these observables include range-equivalent group delays and range-equivalent beat carrier phases. The measurement models for these two observables require refractive ray-tracing calculations through a three-dimensional model of the ionosphere’s free-electron density profile and models of signal reflection off the Earth’s surface. The needed models have been introduced and developed in Baumgarten et al. (2021). This reference also discusses and demonstrates how uncertainties in the electron density profile can induce significant errors in the computation of signal trajectories. Consequently, errors in the computed range-equivalent group delays and carrier phases will occur, degrading positioning accuracy markedly.

A common approach for the reduction of environmental uncertainties relies on the use of data that are obtained from external auxiliary systems. Huang and Reinisch (2006) and Zaalov et al. (2017) demonstrated the use of ionosonde data (Juras, 1985) to generate an enhanced local model for electron density. In Fridman et al. (2006), GPS data inversion computations were used for generating a refined local electron density model. Nickisch et al. (2016) demonstrated the assimilation of HF near-vertical propagation time delay and additional observables in electron density characterization.

A different strategy is considered within the current study—one, however, that does not depend on external aiding for corrections to the ionosphere electron density profile. Instead, this approach assumes that, if sufficiently many group delay and beat carrier-phase measurements are available, then a solution can be obtained not only for the unknown location and clock offset of the receiver, but also for a parameterized model of corrections to the electron density distribution that applies in the vicinity of the propagating signals. A solution for this coupled positioning/ionosphere characterization problem is obtained using a method that combines an a-priori ionosphere model, ray-tracing calculations for propagating signals, and model inversion calculations to compute an enhanced estimate for the unknown quantities. The details of this technique are provided in Baumgarten et al. (2021), which also demonstrates the feasibility of the proposed concept.

The present paper focuses on the potential performance of the batch filter of Baumgarten et al. (2021) and on its sensitivity to problem parameters. In particular, a key question concerns the way that positioning accuracy varies with the set of parameters that characterize a given estimation scenario. This set includes the placement of the broadcasting ground stations, the number of signals that are received by the user receiver, and the size of the differences between the truth and the assumed (a-priori) parameterized ionosphere models. The system’s potential performance is evaluated using both Monte Carlo simulation calculations and a linearized Cramér-Rao covariance lower-bound analysis.

A comprehensive description of the fundamental models that lie in the core of this study is provided by Baumgarten et al. (2021). This includes physical models of propagating RF signals in the ionosphere that are implemented using refractive ray-tracing techniques of the Earth’s geometry and magnetic field as well as of the three-dimensional electron density distribution in the ionosphere. Baumgarten et al. (2021) also derived mathematical models for single-hop and multi-hop ray paths. These models combine to yield the batch filter’s models of the measured observables. The reader is encouraged to review this reference for a broader understanding of these topics. Alternatively, Section 2 of the present paper provides summarized descriptions of both physical and mathematical models. The material covered in this section should suffice for a basic understanding of concepts that are considered later in this paper.

The remainder of this paper is divided into five sections. As already stated, Section 2 briefly reviews the physical models that have been developed by Baumgarten et al. (2021). Section 3 reviews the batch filter of Baumgarten et al. (2021) that computes estimates for the unknown receiver location and ionosphere parameters. Section 4 briefly reviews previous results, defines this paper’s performance analysis methodology, presents detailed results for three representative test cases, and briefly discusses 18 additional cases that have been considered. These results are further discussed and analyzed in Section 5. The discussion focuses on positioning accuracy as well as the fidelity of the a-posteriori estimate of the ionosphere model. This section also proposes follow-up studies. Section 6 summarizes the study’s developments and draws conclusions about the proposed new system.

## 2 AN OVERVIEW OF MODELS

This section provides a succinct overview of concepts that constitute the foundation of this work: HF signal propagation, Earth and ionosphere modeling, and measurement modeling. Comprehensive discussions about these topics are provided in Baumgarten et al. (2021) and Baumgarten (2018).

### 2.1 Signal Propagation in the Ionosphere and Ray-Tracing Computations

This study considers transmitted RF signals with carriers in the range 2 MHz through 10 MHz. The signals are presumed modulated by a binary phase shift keying (BPSK) pseudorandom code, so that pseudorange measurement is facilitated. Pseudorange measurement accuracy is assumed to be 1 kilometer in 1-sigma, with a signal bandwidth of 100 KHz. Beat carrier-phase measurements are utilized as well, facilitated through the use of a stable internal oscillator and a phase-locked loop, so that the expected accuracy for a range-equivalent beat carrier-phase measurement is 1 meter.

Skyward-propagating HF signals that travel in the ionosphere are refracted towards the Earth. Their trajectories are characterized by a curved path and frequency-and-path-dependent propagation speeds of both modulated code and carrier phase. This propagation mechanism, generally modeled through the Appleton-Hartree equation (Baumgarten et al., 2021), is considered for ray-tracing calculations that lie at the core of this study. The fundamental set of ray-tracing equations that are utilized is provided by Jones and Stephenson (1975) in the form of a coupled system of six non-linear ordinary differential equations. Jones and Stephenson (1975), combining the work of Lighthill (1970) and Haselgrove (1955), proposed different Hamiltonians that could be used in the implementation of these equations. Alternative Hamiltonian formulations are presented by Psiaki (2019).

### 2.2 Earth and Ionosphere Models

The physical environment that surrounds propagating signals is considered in the form of parametric models for the Earth’s shape, the Earth’s magnetic flux, and for the free-electron density distribution in the ionosphere.

The Earth’s shape is modeled using the WGS-84 ellipsoid model (National Geospatial-Intelligence Agency, 2011) whose equation in Cartesian Earth-centered, Earth-fixed (ECEF) coordinates is given in Baumgarten et al. (2021). This model has been chosen for its relative simplicity, although a more realistic method for modeling the shape of the Earth may be required when dealing with real data. Such methods could utilize two-dimensional spline surfaces to approximate existing digital representations of the Earth.

The International Geomagnetic Reference Field’s (IGRF) 11th generation is used to model Earth’s magnetic flux vector field in the vicinity of propagating HF signals. This model is embedded in the ray-tracing computations that are utilized in this study. Additional information about the IGRF model can be found in IAGA Division V-MOD (2019). It should be noted that, in the scope of this simulation-based analysis, unlike when working with real data, results are not expected to be sensitive to the particular IGRF and electron density models used.

Ionosphere modeling plays a significant role in this study due to the strong dependence of the propagating signals’ directions and speeds on ionosphere characteristics, namely electron density and magnetic field. The electron density distribution is modeled starting from a three-parameter Chapman beta vertical profile. It was expanded into a three-dimensional spatial model that contains electron density data for the entire near-Earth upper atmosphere using latitude and longitude maps of the vertical profile’s three parameters. The extent to which the Chapman profile can be regarded as the exact solution for an ionosphere density profile is discussed in Baumgarten et al. (2021) and Stankov et al. (2003).

A function *N _{e}*(

*r*) is defined that maps the three-dimensional ECEF Cartesian location vector

*r*to its corresponding electron density value in units of electrons/m

^{3}. This mapping procedure relies on a preceding computation of auxiliary quantities: altitude of maximum electron density,

*h*

_{max}[

*r*]; vertical total electron content,

*VTEC*[

*r*]; and altitude scale factor,

*h*[

_{sf}*r*]. These quantities constitute the Chapman profile’s set of three parameters (Baumgarten, 2018). Note that this computation implicitly accounts for transforming the location’s Cartesian representation,

*r*, to geodesic latitude, longitude, and altitude coordinates. The latitude and longitude dependencies of the three Chapman profile parameters were modeled using a special bi-quintic spline which works with data that are defined at spline nodes. Spline nodes were placed at predefined latitudes and longitudes along circles of constant latitude. Additional information on spline computations, spline stored data, and nodes placement is provided in Baumgarten et al. (2021), Baumgarten (2018), and Psiaki et al. (2019).

Note that distinct D and E layers cannot be represented using a single Chapman profile. Degraded results using real daytime data are likely unless the paper’s models are modified to include multiple Chapman profiles that could model the D and E layers in addition to the main F layer modeling that is assumed here. However, the use of a Chapman profile appears to be a reasonable choice for this stage of simulation-based analysis in which this method’s potential accuracy is evaluated.

### 2.3 Ray Paths

Signals are assumed to be transmitted from ground stations with known locations. The signals follow curved, refracted paths in the ionosphere and are bent back towards the Earth. They may be reflected from the Earth’s surface and refracted back from the ionosphere multiple times before finally arriving at the location of a receiver. It is assumed that the signals are reflected from the Earth’s surface in a lossless specular manner at *bounce points*. These bounce points’ locations are determined as part of a ray-path solution process, as discussed in depth in Baumgarten et al. (2021) and Baumgarten (2018).

Bounce points are connected through curved signal trajectories that are termed *ray-hops*. Ray-hops additionally connect the location of a transmitter with the first bounce point, and the last bounce point with the location of the receiver. Calculation of individual ray-hop trajectories, given a ray-hop’s known start and end locations and given a set of applicable ionosphere parameters that characterize the electron density profile in the vicinity of that ray-hop, was performed using an iterative numerical solution method that has been developed for a two-point boundary value problem (TPBVP). This method is discussed in Psiaki (2019).

Sequenced ray-hops that connect the transmitters, intermediate bounce points, and the receiver constitute *ray paths*. Their characterizing attributes, including total group delays, range-equivalent beat carrier phases, and these quantities’ sensitivities to inputs, are required for the computation of a receiver’s location. Determination of these quantities involves finding the solution of coupled, non-linear equations that define the physical characteristics of the signals’ trajectories. As part of the solution, the location of a ray-path’s bounce points is determined using an algorithm that is called a *ray-path solver*. The ray-path solver is input with locations for the signal’s start and end points that are assumed to be fixed and known. It is additionally given fixed ionosphere parameters. Additional details about the ray-path solver, along with a discussion about ray paths’ feasibility and solution uniqueness, are provided in Baumgarten (2018) and Baumgarten et al. (2021).

### 2.4 Measured Observables

In the context of this study, range measurements are based on signal propagation time and on wave phase measurements. A single transmitter may generate multiple ray paths that differ in their transmission frequencies or in the number of hops of each ray path.

The first type of observable that is utilized in this work is the ray path’s total range-equivalent group delay, which equals the measured signal propagation time multiplied by the speed of light, *c*. Propagation time is computed as the difference between the measured reception time according to the erroneous receiver clock and the true transmission time at the beacon transmitter, whose clock is assumed to be pre-calibrated.

Beat carrier phase is the second type of measurement that is used in this study. It is defined as the difference between measured changes in the received signal’s phase during a fixed period of time and the phase change that occurs in the transmitted signal over the corresponding period of time. This measured quantity is equivalent to the negative of the time integral of the received carrier Doppler shift (Bennett, 1968). Baumgarten (2018) shows how the use of a frequency-hopped continuous-phase signal and multiple beat carrier-phase measurements at multiple times associated with different frequencies can be used to enable estimation of the unknown bias term that is typical for this type of measurement.

## 3 BATCH ESTIMATION OF RECEIVER POSITION, RECEIVER CLOCK OFFSET, AND IONOSPHERE PARAMETER CORRECTIONS

The solution for the unknown receiver location and corrections to a local model for electron density was obtained through the minimization of a cost function that considers the differences between measured pseudoranges and beat carrier phases and their computed counterparts. In addition, it includes a term for the differences between a-priori and current estimates of parameters that characterize the electron density distribution in the ionosphere.

### 3.1 Batch Filter Problem Definition and the Solution Algorithm

The batch filter is designed to minimize the cost function: 1

where is a vector containing three coordinates for the unknown location of the receiver, a term for the unknown receiver clock offset, and potentially, unknown beat carrier-phase measurement bias terms. The vector contains the unknown ionosphere parameters—these are the parameters that characterize electron density distribution in the vicinity of the true ray paths, as described in Baumgarten (2018). The vector characterizes an assumed model for the ionosphere. It holds a-priori values for . The vector contains the measured pseudoranges and range-equivalent beat carrier phases for the given ray paths. The matrix *R* is the measurement noise covariance matrix and *M* is the covariance matrix that models the uncertainty in the a-priori ionosphere parameter vector . The scaling parameter *ζ* is used for dynamic relative weighting of the two terms whose sum constitutes the cost function. It is important to note that the right-hand side term of Equation (1) does not include a-priori values of the elements of so that no prior information about the receiver position is assumed. Additional details about this weighted nonlinear least-squares cost function are provided in Baumgarten et al. (2021) and in more detail in Baumgarten (2018).

This batch estimation problem is solved using the Gauss-Newton method. In its basic form, this method is described in Gill et al. (1981) and Nocedal and Wright (2006). The actual implementation of this method, however, includes several adaptations that address certain characteristics of the present estimation problem. One requirement for the modified method was an ability to handle situations in which intermediate computations of a signal’s group delay and beat carrier phase were unsuccessful. This might be the case if the current guess for one or more of the unknown parameters was poor, rendering a problem setup that includes physically infeasible ray paths. This requirement has been addressed through the development of a dynamic comparison mechanism that evaluates the cost reduction for varying sets of measurements. A second adaptation is the inclusion of a measurement rejection mechanism designed to detect and reject bad measurements using likelihood tests. The prime motivation for such a mechanism originated from the need to identify measurements that were contaminated by significant, yet un-modeled, measurement errors that affect sensor readings. A third adaptation, which has been designed to ensure a high rate of solution convergence, is the use of a two-stage method. In this method, the algorithm only allows limited corrections to the set of unknowns in its early iterations. Additional details on these modifications to the Gauss-Newton method are provided in Baumgarten (2018).

### 3.2 Estimation Errors and the A-Posteriori Covariance Matrix

The estimation errors for the vectors and are defined as: 2

where and are the filter’s estimates that minimize the cost function of Equation (1) and and are the corresponding unknown true values. The computed set of range-equivalent group delays and beat carrier phase, , can be related to through a first-order approximation that is computed by expanding in a Taylor series about and : 3

This approximation can be used to determine an approximate formula for the dependence of the estimation errors on the measurement error vector and the a-priori ionosphere parameters’ error vector . The resulting predicted estimation error takes the form: 4

where *N _{x}* is the number of elements in , and

*N*is the number of elements in .

_{p}Estimation error covariance analysis considers two different cases that will be discussed later in the performance analysis section. In the first case, for which the ionosphere model error vector is sampled from a Gaussian distribution that satisfies and , the estimation error covariance matrix is given by the standard form: 5

This covariance matrix constitutes the Cramér-Rao lower bound for the batch filter’s mean-square estimation error. In a second case that assumes a constant , the covariance matrix is given by: 6

and the estimation error mean value is non-zero. The corresponding mean square error matrix takes the form: 7

## 4 EVALUATION METHODS AND RESULTS FOR BATCH-FILTER TEST CASES

### 4.1 Review of Previous Results

Baumgarten and Psiaki (2017) considered a simplified, straight-segment, ionosphere-thin-sheet-reflections ray-path model. The nine test cases that were studied in that work were characterized by differing numbers of ray paths, differing discrepancies between the true and a-priori ionosphere models, as well as differing ground station locations. The study considered a limited number of batch filter solutions. Preliminary results suggested that the paper’s simplified-model problem was sufficiently observable to make such a system a candidate for navigation. Position errors, ranging from tens to thousands of meters, appeared to be consistent with the corresponding computed Cramér-Rao bounds. At the same time, the filtered estimates of the ionosphere electron density profile parameters tended to have significantly reduced errors in comparison to the a-priori estimates. A principal goal of the present study has been to explore whether such results hold up with much more complex and realistic models of the ionosphere and signals’ refractive ray paths.

### 4.2 Ionosphere Parameter Variability Matrix and the Ionosphere Error Index

The ionosphere parameter variability matrix is an empirical covariance matrix that models the manner at which ionosphere parameters vary over time as well as the correlations between them. The method used for generating this matrix, and its role with the batch estimator, are discussed in Baumgarten et al. (2021).

In the context of performance evaluation, the ionosphere parameter variability matrix takes two roles. First, it is used in generating random a-priori ionosphere models in the form of random parameter perturbation vectors. These models are used with the type of test case analysis that has random a-priori ionosphere variations, demonstrated later with Test Case H0.

The ionosphere parameter variability matrix is additionally used in evaluating the magnitude of errors in the estimates of a-priori and a-posteriori ionosphere model parameters. Error vectors’ squares are normalized by this matrix’s inverse, followed by taking the 10-based logarithm of the resulting scalar quantity. The resultant non-dimensional valued *Ionosphere Error Index* (IEI) can be regarded as small or big, indicating small or significant differences between evaluated ionosphere parameters’ values and their true counterparts, respectively. See Baumgarten (2018) for additional statistical characteristics of the IEI.

### 4.3 Performance Evaluation Methods

Two analysis methods have been used in this study. In the first method, simulation-based analysis is conducted through the examination of individual batch-filter solutions that operate on simulated input data. Error statistics are generated by processing the outcome of sets of Monte Carlo runs that consider different random measurement errors for a given set of input parameters and, possibly, different yet constant random differences between the a-priori and true ionosphere model parameters. This approach is useful in assessing the batch-filtering algorithm’s performance given particular IEI values. As such, it provides in-depth yet somewhat narrow-in-scope information.

The second type of analysis assumes that the difference between the true and a-priori ionosphere parameter vectors is a random vector sampled from a zero-mean Gaussian distribution. This addresses the need to account for ionosphere errors’ randomness at the cost of certain required simplifications that are applied to the assumed ionosphere’s error model. In this approach, error statistics are generated using the linearized covariance analysis that was presented earlier.

The first of the two analyses is similar to a traditional Monte Carlo analysis that uses a truth-model simulation of a filter. The resulting error statistics that characterize the filter’s performance are, therefore, regarded as having high fidelity. The large number of ionosphere parameters and the large cost in filter computation time per simulation, however, makes it impractical to run a sufficient number of simulations to fully explore all the effects of randomness in the measurement errors and randomness in the differences between the a-priori ionosphere model and the “truth” ionosphere model. Therefore, the second covariance-based analysis was ultimately used to further characterize the possible effects of ionosphere model uncertainty using, in effect, a traditional Cramer-Rao-lower-bound-type analysis of the possible filter accuracy.

Our assessment of the proposed navigation system’s effectiveness evaluated the batch filter’s performance in terms of positioning accuracy and its ability to estimate corrections to an erroneous a-priori model of the ionosphere. Assessment was performed through a comprehensive analysis of 21 test cases that differed in both analysis method and various input parameters.

The 21 test cases that were studied are divided into four classes. These classes differ by the type of measurements that were processed in each and by the analysis method that was used for assessing the filter’s performance. Test cases of Class 1 and Class 2 considered a fixed ionosphere model, whereas test cases of Class 3 and Class 4 considered randomly varied ionosphere models consistent with the two types of analysis that have been described. Both group delay and beat carrier-phase measurements were processed with the test cases of Class 2 and Class 4, while only group delay measurements were processed with the test cases of Class 1 and Class 3.

Each class of test cases was divided into two-to-three groups that differed by the number and placement of ground transmitters and by the number of signal measurements that were assumed to be available for the receiver. Test case groups of low signal availability included scenarios where the receiver was assumed to process 21 different signal measurements. Test case groups of medium signal availability considered scenario setups wherein the receiver was assumed to process up to 68 different observables. In the test cases that were included in groups of high signal availability, it was assumed that either 128 or 256 observables, including group delay and beat carrier-phase measurements, were processed by the receiver. In the latter two cases, transmitters were assumed to be broadcasting signals with varying frequencies. Finally, each group of test cases contained test cases that differed by setup parameters that included ray-path geometry, the number of hops for each ray path, signal frequencies, true and a-priori ionospheric models, measurement noise, receiver clock error, and the true location of the receiver. Additional details about all 21 cases, including a description of the particular setup used with each test case, are contained in Baumgarten (2018).

The resulting ensemble of test cases allowed for a comprehensive study of the system’s performance and performance sensitivity to setup parameters. Of the studied test cases, three representative examples are considered in depth in the following three subsections. These consist of a representative test case (Test Case E0), a test case wherein performance was evaluated for a more challenging setup (Test Case D2), and a test case in which the impact of a-priori ionosphere variations on performance were considered and assessed (Test Case H0). The rest are discussed in brief summary at the end of this section for the sake of brevity.

### 4.4 Truth-Model Simulation

Algorithm validation and performance assessment were performed using a truth-model simulation. The simulation was designed to enable testing of any desired combination of ground station placement, ray-path characteristics, measurement error models, ionosphere error models, and other parameters. It has been shown by Baumgarten (2018) that not all such combinations yield physically feasible ray paths. The feasibility of the ray paths of a given configuration is tested during the first stage of a simulation study.

Figure 1 illustrates a typical scenario. It shows the different signals’ curved ray paths, starting from ground stations and eventually arriving at user equipment (UE). Different ray paths transmitted from the same ground station are shown in different shades of green with gray circles denoting their ground bounce-points. The blue circles in the figure denote ground broadcast stations, with the corresponding broadcast signals’ identifying indices shown next to the circles. The magenta diamonds with adjacent three-digit numbers denote ionosphere model bi-quintic spline nodes with their identifying numbers next to them. The North American coastline is shown as a thin blue line. The receiver’s true location is marked with a thin black X at the convergence point of the various ray paths.

Computation for an electron density distribution truth model utilizes a Chapman profile that is fit to an International Reference Ionosphere (IRI) model for a particular time. The model used in the simulation utilized the 2012 release IRI Fortran code available on the official NASA website (Bilitza, n.d.). A Chapman vertical electron density profile is fit to the entire IRI vertical profile, including the F layer, at a set of latitude and longitude points, wherein each such Chapman profile was assumed to be the truth profile at the given latitude and longitude. The fit is performed using a nonlinear least-squares technique (Baumgarten, 2018). This fitting procedure is carried out at each bi-quintic spline node to determine the three parameters. The three-dimensional ionosphere model is defined by the three Chapman vertical profile parameter values’ natural logarithms at each of the latitude/longitude nodes along with all of their latitude and longitude partial derivatives of up to second order in latitude, longitude, or both. This parameterization and its calculation from the IRI model are described in Psiaki et al. (2015, 2019). It should be noted that the Chapman model ignores the possibility of distinct D and E layers, including a sporadic E layer. This level of simplification would likely produce unsatisfactory results if working with real daytime data, but it is reasonable to use a Chapman profile at this stage of simulation-based study of the proposed system’s potential accuracy.

The simulation uses truth values of the receiver location, receiver clock offset, and ionosphere model parameters in the pseudorange measurement model for each ray path. It also uses these truth values, along with the true beat carrier-phase measurement bias terms, in the beat carrier-phase measurement model. These measurement values are input directly into the main batch filtering algorithm with or without random measurement errors added, depending on the type of simulation-based evaluation that is being performed. The simulation also generates an a-priori estimate of the ionosphere parameter vector for use in the cost function of Equation (1). This a-priori vector is typically different from the truth vector. The method of generating appropriate differences, perhaps differences that are even a bit larger than one would expect in a real situation, is to use the IRI model to generate via Chapman-profile fitting at a different date than the date used to generate the truth using the same fitting technique. Such a choice ensures that the truth-model simulation is not using an unreasonably optimistic model of how well the filter’s known would approximate the truth ionosphere.

### 4.5 Test Case E0: A Representative Case

Test Case E0 considered an array of 11 ground station transmitters at various locations across the continental United States (CONUS). This set of ground transmitter stations was a subset of an array of ground station transmitters that is based on a grid of small circles of constant latitude, spaced 5 degrees apart (Baumgarten et al., 2021). The nominal longitudinal difference between two neighboring stations that lie on the same small circle is 10 degrees, though not all nominal station slots on a given small circle of latitude are occupied by transmitters. The nominal longitude grid of stations that lie on two neighboring small circles are offset by 5 degrees. The user receiver was located at latitude/longitude/altitude (LLA) [40.1°, −95.1°, 10,000 m]. This setup is similar to that shown in Figure 1. The total number of nominal signal paths was 32. Each signal path was assumed to be received at four closely spaced times, each with a different carrier frequency, so the total number of effective received signals was 128. A common beat carrier-phase bias was assumed to apply for each set of four signals from the same transmitter with the same number of bounces. Thus, there were 32 unknown beat carrier-phase biases. Both pseudorange and beat carrier-phase measurements were processed. Therefore, the total number of observables was 256. The HF signals for this test case had frequencies in the range 3.0–6.0 MHz. The number of hops for each ray path was 1–4, with a mixture of signals arriving from above the UE and from below the UE.

This test case considered a fixed ionosphere model for which the truth ionosphere electron density profile was based on the IRI model computed for October 23, 2009, at UTC 14:22. The a-priori model was based on the model computed for September 23, 2009, at the same UTC, such that the total seasonal discrepancy was one month, and the corresponding IEI (*ξ*) took a medium value of −0.2276, as with many other test cases that have been studied. It should be noted that these dates in the year 2009 represent a time of solar minimum.

It is reasonable to assume that both the truth and a-priori ionosphere electron density models would have taken significantly different values had this setup considered a year of solar maximum, such as the year 2003. The use of models that had been generated for a different level of solar activity would not necessarily have had an impact on the basic results of this type of case-based study, except that they would likely change the actual geometry of the various simulated ray paths. There are, however, open questions about the fidelity with which a batch filter’s parameterized ionosphere model could represent the true electron density spatial distribution and about whether the possible fidelity depends on the time point within a solar cycle. The investigation of such questions is left for future study.

Figure 2 plots position errors obtained with a 100-run Monte Carlo analysis that used 100 different random measurement error instantiations but a single constant difference between the a-priori and truth ionosphere parameter vectors. The mean position error for Test Case E0 was about 5 meters in the north direction and less than 0.3 meters in both the east and vertical directions. The lengths of the horizontal 90% error ellipse’s semi-major and semi-minor axes were 24 and 13 meters, respectively. The standard deviation for the vertical error was 0.8 meters. The biases in this figure result from the constant error between the a-priori and true ionosphere parameter vectors. The scatter and the corresponding error ellipses were the result of the random measurement errors.

The a-posteriori ionosphere model estimation errors could be assessed by considering the summary statistics that are shown in Figure 3. Each of the three panels plots histograms for the 100 Monte Carlo cases with a separate panel for each of the three Chapman profile parameters. Each panel plots two histograms for the corresponding Chapman profile parameter—a blue one labeled the 80% histogram and a gray one labeled the 95% histogram. Note how each individual histogram’s bar heights cumulatively add up to 100, the total number of Monte Carlo cases.

Each histogram point was calculated based on another statistical analysis as follows. For the corresponding Monte Carlo sample, the errors in the given Chapman profile parameter were calculated on a dense set of latitude and longitude points over the CONUS. The cumulative distribution of these errors is plotted for the Monte Carlo sample, and two points are read off of this cumulative distribution—its 80% point and its 95% point. The corresponding errors indicated that 80% or 95% of the CONUS had errors no greater than the respective value in it’s a-posteriori ionosphere estimate. For example, by considering the 80th percentile histogram in the middle panel of Figure 3, one can infer from the heights of the left-most two bars (34 and 50 Monte Carlo counts, respectively, with the 50-count bar at 0.32 km) that in 34 + 50 = 84 of 100 Monte Carlo runs, a-posteriori h_{sf} errors had values of 0.32 km or less over at least 80% of the area above CONUS. The blue dashed line and the gray dotted line mark the a-priori ionosphere model error values for the 80th and 95th percentiles, respectively. Hence, residual a-posteriori estimation errors for the three Chapman parameters were remarkably small compared to their a-priori values for all 100 runs of this Monte Carlo analysis. That is, the solid blue and solid grey histograms in all three panels of Figure 3 lie well to the left of each panel’s corresponding dashed blue and grey vertical lines.

Another method of assessing ionosphere estimation accuracy considers maps that plot estimation error residuals for the three Chapman parameters. These plots are generated for nominal test case runs (i.e., runs for which no errors were added to the raw measurements). The only source of error was the difference between the truth and a-priori ionosphere parameter estimates. Figure 4 plots errors from truth for the a-priori (top) and the a-posteriori (bottom) estimates of the vertical total electron content parameter, VTEC. These errors have been computed and plotted for a region that contained all active grid nodes for this test case. In other words, this is the region where the ionosphere has been probed by propagating signals. Other regions of the plot have been left blank/white. The red square at latitude/longitude (40.1°, −95.1°) indicates the true position of the receiver. Blue circles with white edges denote the locations of the ground stations. Magenta diamonds denote the locations of the bi-quintic spline grid nodes. The small green squares mark computed truth Earth bounce-points. North America’s coastline is shown in white with the borders of the states shown in gray.

It is evident that the a-posteriori errors in the VTEC were dramatically reduced relative to their a-priori values above the vast majority of CONUS. For the a-priori data, 80% of the errors above CONUS were below 1.25 TECU and 95% were below 1.93 TECU. For the a-posteriori (estimated) model, 80% of the errors were below 0.04 TECU and 95% were below 0.09 TECU.

### 4.6 Test Case D2: One With Fewer Available Signals and a Worse A-Priori Ionosphere Model

Test Case D2 has characteristics that are significantly different from those of Test Case E0. Unlike the previous test case that had 32×4 ray paths, limited signal availability was assumed in this case so that the total number of available ray paths was 17×4. Moreover, a more significant difference between the a-priori and true ionosphere models was assumed, resulting in an IEI value of 0.0848. Baumgarten (2018) showed that this value indicates an excessive error in the a-priori ionosphere model, suggesting that this test case represents a scenario of extremely poor a-priori knowledge of the ionosphere.

The position error scatter plot of Figure 5 is very different from that of Test Case E0 in Figure 2. It is characterized by a significantly larger mean horizontal position error. In addition, its 90% poison error ellipse is about three times as large as for Test Case E0 as the lengths of the 90% error ellipse’s semi-major and semi-minor axes were 75 and 40 meters, respectively. The standard deviation for the vertical error was 2 meters.

The accuracies of the ionosphere corrections are similarly inferior to those presented for the previous test case. For the a-priori data, 80% of the VTEC errors above CONUS were below 1.75 TECU and 95% were below 3.02 TECU. For the a-posteriori (estimated) model, 80% of the VTEC errors were below 0.10 TECU and 95% were below 0.31 TECU (Baumgarten, 2018). This still represents a significant improvement of the estimated ionosphere model relative to the a-priori model, but not as much improvement as in the previous case.

### 4.7 Test Case H0: Consideration of Random A-Priori Ionosphere Variations

A different perspective on performance was obtained by expanding the scope of the test cases so that their a-priori ionospheric model errors were taken as a random vector rather than a constant one. Test Case H0 is the random-ionosphere equivalent of Test Case E0. It considered the same array of ground stations, the same nominal ray paths, and the same receiver location. Group delay and beat carrier-phase measurement noise had the same statistics as those of Test Case E0.

The a-priori ionosphere parameter vector error model utilizes a covariance matrix of the form *γM*_{0}. The scaling factor *γ* takes five different values: 1, 0.5, 0.1, 0.001, and 10^{−9}. The first case, for which *γ* = 1, is a worst-case scenario in which the uncertainty of the ionosphere model is assumed to equal the original ionospheric variability matrix, *M*_{0}. The last case could be regarded as a case of extremely good knowledge of the ionosphere.

Figure 6 plots the horizontal and vertical 90% error ellipses for Test Case H0. These were generated using the Cramér-Rao covariance calculation of Equation (5), as described in Baumgarten (2018). Each panel’s five ellipses correspond to the five different values of *γ*. For the case of *γ* = 10^{−9}, the resulting position errors were solely due to random measurement errors. The upper panel indicates that horizontal position accuracy was sensitive to *γ*, while the lower panel indicates that vertical accuracy was somewhat insensitive to *γ*.

Position error results can be compared to the results that were obtained for the fixed-ionosphere-model Test Case E0. The position error pattern in that case was characterized by a fixed-ionosphere-induced mean error that was relatively small compared to measurement-noise-induced errors (see again Figure 2). The resulting pattern is similar to that of the present Test Case H0 for *γ* = 0.1.

Figure 7 presents ten 80th percentile value maps above CONUS for the estimation errors in the ionospheric VTEC parameter for Test Case H0. Note that the information that is presented in these plots is different from the information that was presented in Figure 4. The error value plotted at each latitude and longitude point indicates that there was an 80% probability of having a VTEC error below the plotted value. Figure 4, on the other hand, plots the actual errors that occurred for a given particular error between the a-priori ionosphere parameter vector and the corresponding truth vector.

The left-hand column of Figure 7 shows 80th percentile error maps for the a-priori estimates of VTEC, while the right-hand column contains maps for its a-posteriori 80th percentile errors. Each row is associated with a different value of *γ*—so, the top row is associated with *γ* = 1 and the bottom row is associated with *γ* = 10^{−9}. Note the different color-coded scales for the ten maps of Figure 7. This approach has been favored over using the same color-coded scales for neighboring panels in an effort to keep the map plots as informative as possible. The reader should note that, in each row, the right-hand plot of the a-posteriori error map is characterized by values that are significantly smaller than those of the left-hand a-priori error map.

### 4.8 Results for Additional Test Cases

A thorough study of the performance of the batch filter is reported in Baumgarten (2018). It includes an analysis of a series of test cases that differ in the sets of parameters that define them. These parameters include the following: type of available measurements, number of ground stations and their placement, number of ray paths, ray-path geometry, number of hops for each ray path, signal frequencies, true and a-priori ionospheric models, receiver clock error, and the true location of the receiver. Analyses of results for these test cases explore positioning sensitivity to scenario parameters. They also explore the extent to which errors in an a-priori parameterized ionosphere model can be reduced. Additional analyses study other aspects of the filter’s performance. These include filter convergence characteristics, scenario setup feasibility, and algorithm robustness.

For the test cases of Class 1, for which only group delay measurements were processed, positioning accuracy was degraded compared to the equivalent test cases of Class 2, in which both group delay and beat carrier-phase measurements were processed. The mean horizontal error for Class 1’s test cases ranged from hundreds to thousands of meters and the length of their 90% error ellipse semi-major axis was typically in an order of kilometers. In contrast, Class 2’s test cases had both typical mean errors and 90% error ellipse semi-major axis lengths in an order of tens of meters. Similar observations have been made for a comparison between results obtained for test cases of Class 3 and Class 4 that consider random ionosphere model variations. 90% error ellipse semi-major axis lengths ranged from tens to hundreds of meters for the test cases of Class 3, while ranging from a few meters to tens of meters for their Class 4 equivalents. The improved accuracy that was observed for test cases of Class 2 and Class 4 was attributed to the use of beat carrier-phase observables in addition to the pseudorange observables that were used in all classes.

Looking at different groups of test cases within a given class, the impact of having a larger number of available measurements was clearly and consistently evident. This was the case, for instance, when results for test cases of Group D (medium signal availability) and Group E (high signal availability) were compared. The accuracies for Group E’s test cases, which have 32 ray paths and 128 observables, were significantly improved compared to their Group D equivalents, for which only 17 ray paths and 68 observables were available. The lengths of the 90% error ellipse semi-major axis were 75 meters and 24 meters for test cases D2 and E2, respectively. The corresponding horizontal mean errors were 137 and 16 meters. This was also the observation for a comparison between test cases of Group G (medium signal availability) and Group H (high signal availability) of Class 4, for which a random ionosphere model was considered and for which the number of available ray paths differed. The dimensions of the 90% error ellipses that characterized the horizontal positioning accuracy for Test Case H0, shown in Figure 6, roughly tripled for Test Case G0 that had 17 rather than 32 ray paths.

When considering individual test cases that belong to the same group, the sensitivity of the position accuracy to the different parameters that has been listed in Section 4.3 can be evaluated. It has been shown that significant discrepancies between the true and a-priori ionosphere typically result in a biased solution for the receiver’s location. Error may be reduced, however, when a relatively large number of ray paths is available. Additional observations about this set of 21 test cases are contained in Baumgarten (2018).

## 5 DISCUSSION OF RESULTS AND PROPOSALS FOR FOLLOW-UP STUDIES

System performance and its sensitivity to parameters have been studied based on an analysis of 21 test cases and variants that differ in many characteristics, as presented in Baumgarten (2018). Results for this ensemble of cases, of which three have been discussed in detail in the previous section, suggest that the problem is sufficiently observable to make this system a candidate for navigation in GNSS-denied situations. That is, a position solution can be obtained to a reasonable level of accuracy despite uncertainty about the ionosphere. At the same time, the filtered estimates of the ionosphere electron density profile parameters tend to have significantly reduced errors in comparison to their a-priori estimates. Therefore, this method may also be useful for remote-sensing-based ionosphere characterization in cases that the receiver location is known a priori.

### 5.1 A-Posteriori Position Accuracy

The position accuracy for Test Case E0 is navigation grade—within 30 meters horizontal and 2 meters vertical 90% of the time. These results imply that, with a sufficient number of signals received and dual group-delay/beat-carrier-phase measurement processing, the achieved accuracy is adequate for the purpose of navigation and guidance for many significant applications. When a random error model for the ionosphere was considered in Test Case H0, with the same ground transmitters/ray paths setup of Test Case E0, a similar level of accuracy was observed for a *γ* value of 0.1. Based on the plots of Figure 6, it can be concluded that navigation grade positioning, in this case, can be achieved if the true uncertainty for the ionosphere error model can be reliably modeled as *γM*_{0} where *γ* < 0.1.

Results have been considered for several fixed-ionosphere-model test cases that were characterized by a relatively large IEI value and corresponding large a-priori ionosphere model errors. They confirm that, with a poor a-priori ionosphere model, the mean position error (i.e., the ionosphere-induced position bias) is further from zero than with the base test cases that have smaller initial IEI values. When position errors are close to their means (i.e., when error ellipses are relatively small), ionosphere-induced position bias will determine whether the system is capable of providing the required level of accuracy for navigation. As an example, reconsider Test Case D2 presented in the previous section. The enlarged mean error is most likely too large for unaided navigation for many applications. It is caused by the large difference between the a-priori and truth ionosphere models.

Additional random ionosphere error test cases demonstrated how position accuracy is closely related to the magnitude of the a-priori errors (or IEI values) in the ionosphere model. With an exceptionally good a-priori model of the ionosphere (i.e., the case of *γ* = 10^{−9}), position errors are primarily due to measurement noise. It is evident that, as *γ* increases, the horizontal position error increases until it reaches a maximum at *γ* = 1. The dimensions of the horizontal error ellipse are roughly doubled at *γ* = 1 in some test cases relative to the *γ* = 10^{−9} error ellipse; they increase by a factor as large as 10 in other test cases. An important observation is that the impact of increasing *γ* on errors in the local vertical direction is typically very limited. Therefore, ionosphere modeling uncertainty has the largest impact on horizontal position accuracy. Note that this conclusion applies to cases in which the user receiver is above ground level and receiver HF signals both from above and from below.

It has been observed that the impact on position accuracy of having different numbers of available measurements for fixed-ionosphere test cases is significant. Scenarios with low numbers of available signals exhibit notably inferior accuracy, evident in the dimensions of their horizontal 90% error ellipses, and, to some extent, in their vertical accuracy. Having a limited number of available measurements, therefore, sets a clear bound on positioning accuracy. Based on the results presented, it is clear that navigation grade accuracy cannot be achieved with only 17 ray paths available, sampled at four different frequencies, unless the system is given an exceptionally good model for the ionosphere. At the same time, it should be noted that even in the worst-case scenario of *γ* = 1, the horizontal 90% error ellipse for Test Case H0 had a fairly small semi-major axis of 80 meters. This is an important observation that suggests that the anticipated negative impact of very poor a-priori ionosphere models may be alleviated given a sufficiently large number of available ray paths.

The manner in which measurement error standard deviation affects positioning accuracy has been studied through a comparative investigation of the 90% error ellipses that have been obtained for several test cases that differ in their assumed measurement noise standard deviations. For a representative fixed-ionosphere test case for which carrier-phase ranging errors were assumed to be 10 times larger than those of a reference test case, the horizontal 90% error ellipse had semi-major and semi-minor axes that were about 3.2 and 3.4 times larger than those of the reference test case’s semi-major and semi-minor axes, respectively. A similar result was obtained for a variant of the random-ionosphere Test Case H0, for which an increase by a factor of 10 for the presumed beat carrier-phase measurement noise resulted in an increase by a factor of three for the dimensions of the 90% error ellipse’s semi-major axis and semi-minor axis.

The important contribution of beat carrier-phase measurements to position solution accuracy is clearly evident. When test cases that originally used both types of observables were tested with only pseudorange measurements, significantly inferior performance was observed, as reported in Baumgarten (2018). For a representative test case that relies on group delay measurements only, the computed 90% horizontal error ellipse’s semi-major and semi-minor axes were 100 meters and 70 meters, respectively. When beat carrier-phase measurements were added to the batch filter, the 90% horizontal error ellipse’s semi-major and semi-minor axes decreased to 35 meters and 20 meters, respectively. For other test cases, performance degradation resulted in horizontal 90% error ellipses for pseudorange-measurements-only test cases that were up to 15 times larger than those of their equivalent pseudorange-and-beat-carrier-phase-measurements test cases. This result is attributed to the low precision of pseudorange measurements that is caused by the limited useable bandwidth of the spreading code in this type of radio ranging application.

The impact of other problem characteristics on position accuracy, including ground station placement and signal directions of arrival, has been investigated and discussed in Baumgarten (2018). Discussion of these issues has been omitted from this paper for the sake of brevity.

### 5.2 A-Posteriori Ionosphere Model Accuracy

The algorithm has proven successful in reducing errors in the a-priori ionosphere model parameters in all 21 test cases. As one might expect, smaller errors in terms of the latitude/longitude-dependent 80th percentile were observed near the areas where ray paths travel through the ionosphere. These regions where electron density is probed can be identified in the a-priori/a-posteriori Chapman parameter estimates’ plots in Figure 4 and Figure 7 by the green points that designate ground reflection points.

Smaller a-posteriori IEI values for all test cases indicate a-posteriori ionosphere model improvements in comparison to their a-priori counterparts. Computed a-posteriori error 80th percentile upper bounds for a representative random-ionosphere model test case were less than 0.75 km for the *h*_{max} parameter, 0.35 km for the *h*_{sf} parameter, and 0.1 TECU for the VTEC parameter, whereas the corresponding a-priori upper bounds were 30 km for *h*_{max}, 10 km for *h*_{sf}, and 5 TECU for VTEC.

Smaller 80th percentile values for the VTEC parameter were computed for smaller values of *γ*, as shown in Figure 7 for Test Case H0. For the first four cases of *γ*, a significant reduction in the 80th percentile values for the a-posteriori estimates is evident. For *γ* = 0.5, values up to 8.8 TECU were reduced to about 0.38 TECU. For the case of *γ* = 10^{−9}, however, the a-priori 80th percentile values were very small and so were the differences between a-priori and a-posteriori values.

### 5.3 Batch Filter Solution Convergence

Two concerns have arisen involving the batch filter solution algorithm’s ability to arrive at the problem’s optimal solution. The first stems from the fact that, while the Gauss-Newton method’s line search is theoretically guaranteed to converge, it is not guaranteed that it will converge to a global minimum. Significant testing experience, however, indicates that convergence to a solution that is far from the simulated truth values has never been observed. It has been concluded, therefore, that the algorithm is insensitive to its input initial guess and it is nearly guaranteed to converge to its global minimum for simulation-generated test cases. The validity of this statement is demonstrated and discussed in Baumgarten (2018).

The second concern derives from the modifications that were applied to the Gauss-Newton method, including the property of a varying set of measurements that are considered in different steps of the Gauss-Newton process. It has been observed that, for simulation-generated test cases with big IEI values (i.e., with a poor initial guess for the ionosphere parameters), occasional failure in ray path solving attempts is common in early iterations of the Gauss-Newton process, resulting in varying measurement sets. However, this did not prevent the algorithm from converging and arriving at accurate a-posteriori estimates for the unknowns.

### 5.4 Proposals for Further Study

An obvious extension of this study would be to substitute tests involving actual data from a network of HF beacons and receivers for its truth-model simulation tests. Such a network is being deployed in South America (Hysell et al., 2016). The estimation cases that could be addressed with actual experimental data may differ somewhat from the cases that have been studied in this work due to the limited availability of received signals. One interesting case may assume that receiver location and receiver clock error are known a priori, so the problem is reduced to a case of ionosphere estimation only. The present study’s truth-model simulations imply that elimination of the uncertainty that is associated with the receiver location should result in enhanced accuracy of the unknown ionosphere parameters.

A second study should consider enhanced ionosphere modeling. As discussed in Baumgarten (2018), the current ionosphere parameterization does not allow modeling of distinct D and E layers. A greater physical fidelity would be achievable by incorporating increased complexity in the ionospheric model parameterization. Therefore, one useful extension for this work would be to employ a more realistic ionosphere model that would enable representation of the D and E layers. Such a model would need more than a single Chapman curve for its vertical profile and would entail more than three latitude/longitude-dependent parameters in the resulting profile.

A useful variant of this work would be to consider a case in which the estimator lacks a-priori knowledge of the number of bounces for any given signal. This is the case when the number of bounces is unknown and must be estimated for each received signal while estimating the receiver position and clock error along with the ionosphere corrections. Such an estimator would solve a mixed real/integer batch filtering problem. The navigation community already has experience with such problems when working with the technique known as Carrier-Phase Differential GPS (Psiaki & Mohiuddin, 2007).

## 6 SUMMARY AND CONCLUSION

This paper has studied the potential accuracy of a navigation system based on high-frequency signals from ground transmitters. High-frequency signals transmitted from stationary ground-based beacons at known locations propagate in the ionosphere along non-line-of-sight paths via ionosphere refraction. The signals are received at a user’s receiver whose location is unknown. Measured pseudoranges and beat carrier phases were processed to solve a coupled positioning and ionosphere-characterization problem. A nonlinear batch filer employed a modified Gauss-Newton method to estimate the unknown receiver position and clock error, as well as corrections to an a-priori parameterized model of electron density in the ionosphere. The batch filter was proven to be successful in achieving a high rate of convergence to the optimal value of the underlying cost function.

System performance has been investigated using a Monte Carlo analysis which is based on a truth-model simulation. Theoretical Cramér-Rao covariance analysis was also performed. The simulation and the corresponding batch-filter used an advanced ray-tracing model of HF signals that propagate in the ionosphere. 21 simulated test cases that considered various combinations of problem characteristics were studied. Position accuracy is influenced by the level of uncertainty of the a-priori ionosphere model by the number and geometry of available measurements and by receiver ranging error standard deviations.

Results indicate feasibility for the combined HF navigation/ionosphere-correction concept. It has been shown that, with sufficient availability of received signals and sufficient a-priori ionosphere model accuracy, navigation grade accuracy for positioning may be achievable. That is, a 90% horizontal error ellipse with semi-major and semi-minor axes both less than 25 meters and a 90% vertical error bound less than 5 meters may be achievable. Furthermore, a-posteriori ionosphere model estimates consistently improved for all cases in comparison to their a-priori counterparts.

## HOW TO CITE THIS ARTICLE

Baumgarten, Y., Psiaki, M. L., & Hysell, D. L. (2022). Navigation and ionosphere characterization using high-frequency signals: A performance analysis. *NAVIGATION, 69*(4). https://doi.org/10.33012/navi.546

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.