Category: Research & Development

Innovation: Where Are We?

Positioning in Challenging Environments Using Ultra-Wideband Sensor Networks

By Zoltan Koppanyi, Charles K. Toth and Dorota A. Grejner-Brzezinska

INNOVATION INSIGHTS by Richard Langley — **INNOVATION INSIGHTS** by Richard Langley

QUICK. WHO WAS THE FIRST TO PREDICT THE EXISTENCE OF RADIO WAVES? If you answered James Clerk Maxwell, you are right. (If you didn’t and have an electrical engineering or physics degree, it’s back to school for you.) In the mid-1800s, Maxwell developed the theory of electric and magnetic forces, which is embodied in the group of four equations named after him. This year marks the 150th anniversary of the publication of Maxwell’s paper “A Dynamical Theory of the Electromagnetic Field” in the Philosophical Transactions of the Royal Society of London.

Interestingly, Maxwell used 20 equations to describe his theory but Oliver Heaviside managed to boil them down to the four we are familiar with today. Maxwell’s theory predicted the existence of radiating electromagnetic waves and that these waves could exist at any wavelength. Maxwell had speculated that light must be a form of electromagnetic radiation. In his 1865 paper, he said “This velocity [of the waves] is so nearly that of light, that it seems we have strong reason to conclude that light itself (including radiant heat, and other radiations if any) is an electromagnetic disturbance in the form of waves propagated through the electromagnetic field according to electromagnetic laws.”

That electromagnetic waves with much longer wavelengths than those of light must be possible was conclusively demonstrated by Heinrich Hertz who, between 1886 and 1889, built various apparatuses for transmitting and receiving electromagnetic waves with wavelengths of around 5 meters (60 MHz). These waves were, in fact, radio waves. Hertz’s experiments conclusively proved the existence of electromagnetic waves traveling at the speed of light. He also famously said “I do not think that the wireless waves I have discovered will have any practical application.” How quickly he was proven wrong.

Beginning in 1894, Guglielmo Marconi demonstrated wireless communication over increasingly longer distances, culminating in his bridging the Atlantic Ocean in 1901 or 1902. And, as they say, the rest is history. Radio waves are used for data, voice and image one-way (broadcasting) and two-way communications; for remote control of systems and devices; for radar (including imaging); and for positioning, navigation and time transfer. And signals can be produced over a wide range of frequencies from below 10 kHz to above 100 GHz.

Conventional radio transmissions use a variety of modulation techniques but most involve varying the amplitude, frequency and/or phase of a sinusoidal carrier wave. But in the late 1960s, it was shown that one could generate a signal as a sequence of very short pulses, which results in the signal energy being spread over a large part of the radio spectrum. Initially called pulse radio, the technique has become known as impulse radio ultra-wideband or just ultra-wideband (UWB) for short and by the 1990s a variety of practical transmission and reception technologies had been developed.

The use of large transmission bandwidths offers a number of benefits, including accurate ranging and that application in particular is being actively developed for positioning and navigation in environments that are challenging to GNSS such as indoors and built-up areas. In this month’s column, we take a look at the work being carried out in this area by a team of researchers at The Ohio State University.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas. Email him at lang @ unb.ca.

GNSS technology provides position, navigation and timing (PNT) information with high accuracy and global coverage where line-of-sight between the satellites and receivers is assured. This condition, however, is typically not satisfied indoors or in confined environments. Emerging safety, military, location-based and personal navigation applications increasingly require consistent accuracy and availability, comparable to that of GNSS but in indoor environments.

Most of the existing indoor positioning systems use narrowband radio frequency signals for location estimation, such as Wi-Fi, or telecommunication-based positioning (including GSM and UMTS mobile telephone networks). All these technologies require dedicated infrastructure, and the narrowband RF systems are subject to jamming and multipath, as well as loss of signal strength while propagating through walls. In contrast, using ultra-wideband (UWB) signals can, to some extent, remediate those problems by offering better resistance against interference and multipath, and they feature better signal penetration capability. Due to these properties, the use of UWB has the potential to support a broad range of applications, such as radar, through-wall imagery, robust communication with high frequency, and resistance to jamming. Furthermore the impulse radio UWB (IR-UWB), the subject of this article, can be an efficient standalone technology or a component of positioning systems designed for multipath-challenged, confined or indoor environments, where GNSS signals are compromised.

IR-UWB positioning can be useful in typical emergency response applications such as fires in large buildings, dismounted soldiers in combat situations, and emergency evacuations. In such circumstances, the positioning/navigation systems must determine not only the exact position of any individual firefighter or soldier to facilitate their team-based mission, but also navigate them back to safety. Under these scenarios, a temporary ad hoc network has to be quickly deployed, as the existing infrastructure is usually non-functional, damaged or destroyed at that point. The UWB-based systems may easily satisfy these criteria: (1) nodes placed in the target area can rapidly establish the network geometry even if line-of-sight between nodes is not available, (2) the communication capability allows for sharing measurements, and (3) the node positions may be calculated based on these measured ranges in a centralized or distributed way. Once the node coordinates have been determined, the tracking of the moving units can start. Obviously, the resistance against jamming makes this solution attractive for military applications.

Ad Hoc Network Formation for Emergency Response

Quick deployment
Sufficient positioning accuracy
Robustness against interference (jamming)
Signal penetration through solid structures

Generally, positioning systems, both local and global, require an infrastructure, which defines the implementation of a coordinate frame. For example, the national reference frames and their realizations support conventional land surveying, or the satellite and the GPS tracking subsystems, as well as the beacons in Wi-Fi systems. UWB positioning also follows the same logic; the network infrastructure defines a local coordinate system and allows for range measurements between the network nodes and the tracked unit(s).

Ad Hoc Sensor Network: Ad hoc networks are temporary, and thus, the node coordinates are not expected to be known or measured a priori; consequently, they are calculated based on measuring the ranges between the units in the initial phase, and can be updated subsequently if the network configuration changes.

Anchored Networks: The network nodes’ coordinates are known. If only local coordinates are known, then to connect to a global coordinate frame, at least one node’s global coordinates and a direction vector must be known to anchor and orient the network.

Anchor-Free Networks: No node coordinates are known, thus the localization problem is underdetermined. Nevertheless, the problem is still solvable, if it is extended with additional constraints.

Tracking: Once a network is established, static/moving objects can be positioned in the network coordinate system.

Ultra-Wideband Ranging

At the beginning of the 21st century, the Federal Communications Commission (FCC) introduced new regulations that enabled several commercial applications and initiated research on UWB application to PNT. The current FCC rules for pulse-based positioning or localization implementations require the applied bandwidth be between 3.1 and 10.6 GHz and the bandwidth to be higher than 500 MHz or the fractional bandwidth to be more than 0.2.

The typical IR-UWB ranging system consists of multiple transceiver units, including the transmitter and the receiver components. The transmitter emits a very short pulse (high bandwidth) with low energy, and the receiver detects the signal after it travels through the air, interacting with the environment. After reaching objects, the emitted pulse is backscattered as several signals, which likely reach the receiver at different times. In contrast, conventional RF signals are longer in duration, thus the backscattered waves overlap each other at the receiver, forming a complex waveform, and may not be distinguishable individually. Due to the shortness of the UWB signals, measurable peaks are nicely separated, representing different signal paths.

The wave shape of the impulse response of the transmission medium highly depends on the environment complexity due to multipath. Detections in the received wave are determined by a peak-detecting algorithm. Note that the travel time is generally determined from the first detection, as it is assumed to be from the shortest path, although other peak detection algorithms also exist.

In the experiments discussed in this article, a commercial UWB radio system was used. This sensor’s bandwidth is between 3.1 and 5.3 GHz, with a 4.3-GHz center frequency. Three methods are available to obtain ranges: (1) coarse range estimation, based on the received signal strength with dynamic recalibration; (2) precision range measurement (PRM), which uses the two-way time-of-flight technique; and (3) the filtered range estimates (FRE) method that refines the PRM solution using Kalman filtering. In our investigations, PRM data were used in static situations, when both the unit to be positioned and the reference units were static (such as when determining network node coordinates), and FRE was logged in kinematic scenarios.

Localization in a UWB Network

Commercial UWB products usually provide capabilities for all three applications: communication, ranging and radar imaging. In positioning applications, identical units are used for both the rovers — that is, the units to be localized — and the static nodes of the network. The general terminology, however, is that the rover unit with unknown position is called the receiver, and units deployed at known locations are called transmitters. We will also use the terms rover and stations. The positions are typically defined in a local coordinate system. The usual ranging methods used in RF technologies, including signal strength and fingerprinting, time of arrival, angle of arrival, and time difference of arrival, are also applicable to UWB systems. TABLE 1 lists the ranging methods and typical performance levels; the achievable accuracies are based on external references. Note that the accuracy depends on the sensor hardware and network configuration, applied bandwidth, signal-to-noise ratio, peak detection algorithm, experiment circumstances, formation and the environment complexity.

TABLE 1. Typical accuracy of the different UWB localization techniques. Note that the results depend on the hardware, antenna, applied bandwidth, experiment circumstances and geometric configuration; * denotes indoor environment with area coverage of a few times 10 × 10 meters, with line-of-sight conditions, and ** refers to the maximum error in the outdoor test area of about 100 × 100 meters).

Signal Strength. The received signal strength (RSS) requires modeling of the signal loss, which is a challenging problem since signals at different frequencies interact with the environment in different ways, and thus the resulting accuracy is generally inadequate for most applications. The fingerprinting approach is also applied to UWB positioning; the signal-strength vector received from the transmitters identifies a location by the best match, where the vector-location pairs are measured in a calibration/training phase and stored in a database.

Time of Flight. The time-of-flight method requires the synchronization of the clocks of the UWB units, which is difficult, in particular, in the low-cost systems. Therefore, most UWB systems are based on the two-way time-of-flight method, which eliminates the unknown clock delay between the sensors, although it also has its own challenges. The range between two units is obtained by measuring the time difference of the transmitted and received pulses plus knowing the fixed response time of the responding unit.

Computing Position in a Network. Once the ranges are known in a network environment, the position is determined by circular lateration. The principle for the 2D case with three stations is shown in FIGURE 1. Note that each range determines a circle around the known stations (stations 1, 2 and 3 in the figure), thus, if the stations’ coordinates are known, the unknown position can be calculated as the intersection of these circles. The problem is treated as a system of non-linear equations; note that the lateration requires at least three or four nodes in an adequate spatial distribution for 2D and 3D positioning, respectively. The measured ranges, characterized by the error terms usually modeled with a normal distribution, are depicted by the dotted parallel circles around the solid “perfect” range in Figure 1. Note that this is an optimization problem, which can be solved with direct numerical approximation, such as gradient methods, or by solving the respective linear system after linearizing the problem with close initial position values.

Time Difference and Angle of Arrival. The time difference of arrival (TDoA) approach is useful when the time synchronization is not established. The unknown time delays are eliminated by subtracting the travel times between the rover and the stations, and the response time of the responding unit must be known. The location estimation is similar to the time of arrival case, but rather than the intersection of the circles, hyperbolic function curves representing constant TDoA values are used to determine the rover position. Also, if errors are present in the measurements, the position calculation becomes an optimization problem instead of finding the root of an equation. The TDoA can be combined with the angle of arrival (AoA). This method assumes that the set of UWB antennas are arranged in an array, and the angle can be calculated as the time difference of the first and the last detection from different antennas of the array.

Calibration

The ranges obtained by UWB sensors could be further improved by calibration — for example, by estimating antenna and hardware delays. In our outdoor tests, the joint calibration model (see Two Calibration Models box) was used, and coefficients of various model functions were estimated. During these tests, the UWB units were placed at the corners of a 15 × 15 meter area (see FIGURE 2).

At two diagonal corners, two UWB units with a 1.5-meter vertical separation were installed on poles, while at the two other corners only one unit was used. These six units formed the nodes or the stations of the network. In all cases, a GPS antenna was fixed to the top of the poles to provide reference data. A pushcart with two UWB units, a logging laptop computer, a GPS antenna and a receiver formed the rover system. The reference solution was obtained by using the GPS measurements, with the accuracy around 1 centimeter after kinematic post-processing using precise satellite orbit and clock data. During calibration, the pushcart was collecting stationary data at points 1 to 12, marked on a 5 × 5 meter grid, as shown in Figure 2.

Two Calibration Models

Individual sensor calibration is the approach where the sensor delays are determined separately, for example, , where is the measured range between stations A and B, and are the calibration functions, and is the corrected range.
Joint calibration model is the approach where the calibration function does not provide the offset per station, but rather gives the relative offset between the two stations, where .

The calibration model as a function of the measured distance can be constant, linear or a higher-order polynomial.

After acquiring range data between the rover and network stations, three types of joint calibration functions were investigated: constant, linear and polynomial models. The coefficients of these functions were estimated from the measured ranges and GPS-provided reference positions at all grid points. The estimated functions with respect to the six network nodes are shown in FIGURE 3. Our hypothesis was that the accuracy is assumed to depend on the rover-station distance, and thus, the detected discrepancies between the rover and reference points are expected to be higher if the distance is larger. The results indicate that a constant correction (that is, an antenna delay) is generally sufficient, indicating that the calibration may be applicable to similar installations. In some cases, a linear trend (a distance dependency) may be recognized due to slight data changes, but the observed regression lines are either increasing or decreasing, which clearly rejects the distance-dependency hypothesis. The linear and second-order polynomial functions likely model only local effects. The corrections provided by these functions depend on the environment, and consequently, are valid only in that configuration and where they were observed.

Error surfaces, derived as the approximation of a second-order surface from the residuals at the grid points between the receiver and the six station units, show that the discrepancies can be as large as 0.5 meter. Calibrated results using the constant model show that all the discrepancies are less than 10 centimeters with an empirical standard deviation of 3.6 centimeters. This suggests that, at least, the constant-model-based calibration is needed.

Tracking Outdoors and Indoors

If the coordinates of the network nodes and the calibration parameters are known, the location of the moving rover can be calculated with circular lateration. The experiment described in this section is based on the same field test as presented earlier. For assessing the outdoor tracking performance, a random trajectory of the pushcart inside and outside of the rectangle defined by nodes was acquired (see FIGURE 4). The reference trajectory was obtained by GPS and the UWB trajectory was calculated with circular lateration.

TABLE 2 presents a statistical comparison of the coordinate component differences between the GPS reference and the UWB trajectory based on calibrated ranges. The mean of the X and Y coordinate differences are around 0 centimeters, and their standard deviations are 9.7 and 13.2 centimeters, respectively, with the largest differences being less than half a meter in both coordinate components. Note that the vertical coordinates have large errors due to the small vertical angle, which translates to weak geometric conditions for error propagation.

TABLE 2. Statistical results for the coordinate components.

Indoor UWB positioning is more challenging than outdoor, as propagation through walls modifies the RF signals resulting in attenuations and delays. Furthermore, the geometric error propagation conditions (that is, the shape of the network) may also reduce the quality of positioning. In the indoor tests, a personal navigation system demonstration prototype built in our lab (shown in FIGURE 5) was used as a rover. During the tests, the person was moving at a normal pace, and the rover unit recorded the ranges from the reference stations. Concerning the network, two point types are defined: (1) network nodes depicted by a double circle in the figure, which are used in the tracking phase; and (2) reference points marked by a single circle, which support the validation of the positioning results.

Since no reference solution was available during the indoor testing, the calibration method’s consistency was evaluated based on the relative or internal accuracy metric, which is the a posteriori reference standard deviation error:

where v is the vector of residual errors and r=dim(A^TA) – rank(A^TA) is the degrees of freedom of the network with A being the design matrix describing the geometry of the network. The m₀ values are shown in FIGURE 6. This parameter describes the statistical difference of the measurements from the assumed model (circular lateration). The average m₀ is 7.6 centimeters without calibration, and higher if any of the outdoor calibration models are used.

FIGURE 6. The indoor test results showing values of m0 at the epochs.

To estimate the absolute or external accuracy without a reference trajectory, points 1002 and 1004 were used as checkpoints with known coordinates. Obviously, these points were not part of the network. The UWB rover unit was placed at these points, and data were acquired in a static mode. The coordinates were continuously calculated after measuring at least three ranges. TABLE 3 presents the statistical results. Note that the average is not 0, thus the result is biased, indicating that the signal penetration and/or multipath effects are present in this complex indoor environment. Also, note that no calibration was performed, as no indoor calibration results were available, and using the outdoor calibration models only decreased the positioning accuracy. In addition, the standard deviations indicate the average m₀ is consistent with the external error for point 1002, while this hypothesis is rejected for point 1004.

TABLE 3. Differences between the UWB position estimations and the correct coordinates at points 1002 and 1004.

Taking a closer look at the results of point 1004, the ambiguity problem of the circular lateration can be observed. The random measurement error can be large enough to cover two possible intersections in circular lateration, thus the estimator may oscillate between two solutions. Two main causes for this ambiguity are a weak network configuration and the large ranging errors (see FIGURE 7).

Ad Hoc UWB Sensor Network

We have also carried out tests on an indoor ad hoc sensor network using different coordinate estimation methods. Indoor distance measurements typically do not follow a normal or Gaussian error distribution but rather a Gaussian mixture distribution, which demands the use of a robust estimation method. Our results showed that the maximum likelihood estimation technique performs better than conventional least squares for this type of network.

Conclusion

Ultra-wideband technology is an effective positioning method for short-range applications with decimeter-level accuracy. The coverage area can be extended with increasing network size. The technology can be used independently or as a component of an integrated positioning/navigation system. GPS-compromised outdoor situations and indoor applications can be supported by UWB in permanent and ad hoc network configurations. While UWB technology is relatively less affected by environmental conditions, signal propagation through objects or other non-line-of-sight conditions can reduce the reliability and accuracy.

Acknowledgments

This article is based, in part, on the paper “Performance Analysis of UWB Technology for Indoor Positioning,” presented at the 2014 International Technical Meeting of The Institute of Navigation, held in San Diego, Calif., Jan. 27–29, 2014.

Manufacturer

The experiments discussed in the article used a Time Domain Corp. PulsON 300 UWB radio system.

ZOLTAN KOPPANYI received his B.Sc. degree in civil engineering in 2010 and his M.Sc. in land surveying and GIS in 2012, both from Budapest University of Technology and Economics (BME), Hungary. He also received a B.Sc. in computer science from the Eötvös Loránd University, Budapest, in 2011. He is a Ph.D. student at BME and was a visiting scholar at the Ohio State University (OSU), Columbus, in 2013. His research area is human mobility pattern analysis and indoor navigation.

CHARLES K. TOTH is a research professor in the Department of Civil, Environmental and Geodetic Engineering at OSU. He received an M.Sc. in electrical engineering and a Ph.D. in electrical engineering and geo-information sciences from the Technical University of Budapest, Hungary. His research expertise covers broad areas of 2D/3D signal processing; spatial information systems; high-resolution imaging; surface extraction, modeling, integrating and calibrating of multi-sensor systems; multi-sensor geospatial data acquisition systems, and mobile mapping technology.

DOROTA A. GREJNER-BRZEZINSKA is a professor in geodetic science, and director of the Satellite Positioning and Inertial Navigation (SPIN) Laboratory at OSU. Her research interests cover GPS/GNSS algorithms, GPS/inertial and other sensor integration for navigation in GPS-challenged environments, sensors and algorithms for indoor and personal navigation, and Kalman and non-linear filtering.

“Ultra-Wideband Range Estimation: Theoretical Limits and Practical Algorithms” by I. Guvenc, S. Gezici, and Z. Sahinoglu in Proceedings of ICUWB2008, the 2008 Institute of Electrical and Electronics Engineers International Conference on Ultra-Wideband, Hannover, Germany, September 10–12, 2008, Vol. 3, pp. 93–96, doi: 10.1109/ICUWB.2008.4653424.

• Received Signal Strength Fingerprinting

“Increased Ranging Capacity Using Ultrawideband Direct-Path Pulse Signal Strength with Dynamic Recalibration” by B. Dewberry and W. Beeler in Proceedings of PLANS 2012, the Institute of Electrical and Electronics Engineers / Institute of Navigation 2012 Position, Location and Navigation Symposium, Myrtle Beach, S.C., April 23–26, 2010, pp. 1013–1017, doi: 10.1109/PLANS.2012.6236843.

“Indoor Ultra-Wideband Location Fingerprinting” by H. Kröll and C. Steiner in Proceedings of IPIN 2010, the 2010 International Conference on Indoor Positioning and Indoor Navigation, Zurich, September 15–17, 2010, pp. 1–5, doi: 10.1109/IPIN.2010.5648087.

• Ultra-Wideband Time-of-Arrival and Angle-of-Arrival“Ultra-Wideband Time-of-Arrival and Angle-of-Arrival Estimation Using Transformation Between Frequency and Time Domain Signals” by N. Iwakiri and T. Kobayashi in Journal of Communications, Vol. 3, No. 1, January 2008, pp. 12–19, 10.4304/jcm.3.1.12-19.

• Maxwell’s Equations

“The Long Road to Maxwell’s Equations” by J.C. Rautio in IEEE Spectrum, Vol. 51, No. 12, December 2014, North American edition, pp. 36–40 and 54–56, doi: 10.1109/mspec.2014.6964925.

A Student’s Guide to Maxwell’s Equations by D. Fleisch, Cambridge University Press, Cambridge, U.K., 2008.

March 5, 2015

Innovation: Python GNSS Receiver

An Object-Oriented Software Platform Suitable for Multiple Receivers

By Eliot Wycoff, Yuting Ng, and Grace Xingxin Gao

INNOVATION INSIGHTS by Richard Langley

AND NOW FOR SOMETHING COMPLETELY DIFFERENT. My first introduction to computer programming was during a visit to the Faculty of Mathematics at the University of Waterloo when I was still a high school student. We got to keypunch a simple program onto cards using the FORTRAN programming language and submit the “job” to the university’s IBM 7040 mainframe computer. That visit helped seal the choice of Waterloo for my undergraduate education — but in applied physics, not math.

Once I became an undergraduate, I learned how to properly program in FORTRAN (actually FORTRAN IV with the WATFOR compiler developed at Waterloo) and in assembly language on the SPECTRE virtual computer (written in FORTRAN), both on Waterloo’s new IBM 360 mainframe. Knowing how to program was instrumental in my graduate work on the geodetic application of very long baseline interferometry (VLBI) at York University. Being humble Canadians (and despite the fact that VLBI was invented in Canada), we called it just LBI. My LBI data analysis FORTRAN program was initially on a box full of punched cards that I would have to carry back and forth to the computer center being careful not to drop the box and get the cards out of order.

While I was a graduate student, I also got to use the Spiras-65 minicomputer that controlled the playback of the LBI recorded tapes at the National Research Council in Ottawa. It was programmed using punched paper tape.

I saw the progression from punched tape and cards to the use of terminals to enter programs and magnetic tapes for storing them and the data to be analyzed. The University of New Brunswick, where I came to work in 1981, was one of the first universities in Canada to introduce an interactive terminal- (or work-station-) based time-sharing system for programmers to develop and run their jobs on the central computer. The last card reader at UNB was retired in 1987.

By the time I came to work at UNB, the era of the personal computer had already dawned. Although the Department of Surveying Engineering (as it was then called) acquired an HP 1000 minicomputer for various research tasks, personal computers began to show up on faculty members’ desks and in their labs. Some of us started out with Apple II computers (we used them, for example, for recording data from Transit–U.S. Navy Navigation Satellite System–receivers) and progressed through various Macintosh models.

Once I became a professor, I did less and less programming myself–leaving it up to my graduate students to do the heavy lifting in that area. These days, my personal programming efforts are limited to short scripts mostly using the Python language. Python, which gets its name from the Monty Python’s Flying Circus television series, was first introduced back in 1991 but it is only relatively recently that its popularity has taken off. Python can be run on a wide variety of platforms under many operating systems.

One of the key features of Python is that it supports multiple programming paradigms, including object-oriented programming (OOP). OOP is a programming methodology based on the use of data structures, known as objects, rather than just functions and procedures. The objects, organized into classes, exchange information in a standardized way and their use helps ensures good code modularity.

In this month’s column, we take a look at how Python has been used to develop a software-defined GNSS receiver — one well-suited to processing data from a network of receiver front ends.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas. Email him at lang @ unb.ca.

With billions of GNSS-enabled devices in use today, the potential gains from harnessing data collected over a network of GNSS receivers has never been greater, yet the necessary architectures to handle and extract useful data collected over such networks are not well explored. Traditional uses of GNSS in cooperative positioning treat individual GNSS receivers as “black boxes” that merely output navigation solutions. As such, the wealth of information contained in each receiver’s raw signals is largely discarded.

Of particular interest are ideas such as inter-receiver aiding, in which networked receivers might share acquisition, tracking, and navigation information (possibly in real time) to improve receiver performance. In addition, a network of receivers might also be used as a sensing tool: it is expected that atmospheric parameters, for instance, could be recovered by analyzing the raw signal data arriving at an appropriately sized network.

In light of these interesting research areas, it would be expedient to develop a set of tools that can process and handle the raw data being produced at every receiver in a GNSS receiver network. Existing software-defined receivers (SDRs) have gone a long way towards making the fast prototyping of new receiver architectures possible. An SDR attempts to shift as many receiver functions, such as mixing and tracking, from being implemented in hardware to being implemented in software. This allows for fast prototyping as receiver components can be more quickly modified in software than in hardware. The hardware components that a GNSS SDR still requires are an antenna and a front end including an analog-to-digital converter (ADC). An analog GNSS signal is received at the antenna. It is then mixed to an intermediate frequency and digitized by the ADC. The digital stream is then processed by the SDR’s software component.

But with regard to processing data from a receiver network, existing SDRs have a number of notable flaws. In brief, existing software receivers are designed to process the data arriving at one real-world receiver. Thus a procedural coding design is typically used. While procedural code is a good solution for the linear processes that occur in a single receiver (acquisition, tracking, demodulation of the navigation data, position calculations, and so on), this software design style does not adapt well to the task of performing all of these actions on multiple receivers with the additional goal that each receiver shares tracking data with every other one. In such scenarios, not only is there data being produced for every receiver in the network, but there is also data being produced about the relationships between the receivers in the network. Thus, an SDR that was originally designed to process data from only one receiver will prove difficult to adapt to the task of processing many.

Luckily, object-oriented programming, a well-known and widely used software design philosophy, is well suited to the receiver network problem. Therefore, for this work, we designed and implemented an object-oriented software platform for many receivers. Python was chosen as the programming language because of its support for object-oriented programming, its portability, its free cost, its numerical abilities (using open-source libraries such as NumPy and SciPy), and its ease of use. And as a reference, an existing Matlab software receiver was used as a basis for developing many of the core algorithms in this work. We call our development simply the Python Receiver.

Design

Many of the core functions in the Python Receiver are modeled after those found in the Matlab development. Thus, this particular implementation is suited for the raw GPS L1 signal data mixed to an intermediate frequency by the SDR front end. In addition, the basic algorithms for acquisition, scalar tracking, and navigation are similar to the Matlab ones, with the exception that acquisition is made more robust by using multiple noncoherent integrations. The primary innovation of this software, however, is in the way in which the code is organized. For tracking multiple receivers, the Python Receiver was designed under an object-oriented approach.

FIGURE 1 illustrates the main objects that a user would be expected to use in the Python Receiver. Each object is defined as a class, and as such each object is capable of storing object-specific data as well as performing certain object-specific functions. The hierarchy of Figure 1 roughly illustrates which objects are defined as members of other classes for typical usage. Thus, inside any instance of the network class may exist any number of receiver objects. Likewise, an instance of the constellation class may be home to any number of satellite objects.

FIGURE 1. Typical object (class) hierarchy.

For data coming from a single real-world receiver, use of the Python Receiver would typically be as follows. First, a user would initialize an instance of the receiver class using a dictionary of predefined settings, such as the file location of the data source. Second, the user would initialize a constellation object of satellites by passing the pseudorandom noise (PRN) code values of each satellite to be included in the constellation. At this point, the user could then use built-in functionality in the receiver object to perform acquisition of all of the satellites in the constellation. Results of this acquisition attempt would be stored in the receiver object, where they could then be used to run the receiver’s built-in scalar tracking functionality. Likewise, scalar tracking data would be stored in the receiver object, and again the user could use the receiver’s built-in navigation functionality to decode the navigation bits produced during scalar tracking and perform navigation computations. Satellite-specific ephemerides would be stored in the relevant satellite objects.

Navigation solutions are stored as a part of the receiver’s state object. The state object, which is also used in the satellite class, is a container for holding state information in the Earth-centered Earth-fixed (ECEF) coordinate system (such as position and velocity) and clock terms, and it also provides the ability to return position coordinates in other systems, such as the GPS geodetic system (frame) of WGS 84. While it is not a key feature of the Python Receiver, the state object is designed as an object so that it can be readily used elsewhere should an algorithm need to store state information and have coordinate transformations readily available.

Tracking channels need not be restricted to the hierarchy shown in Figure 1. During operation for just one data source, the scalar tracking function defined at the receiver level will initialize a sufficient number of tracking channels to track all of its observed satellites. However, when operating on multiple sources of data and with the intent to share tracking outputs between channels, it is helpful to place tracking channels into groups, as shown in FIGURE 2. In the example that will be discussed in following sections, two real-world receivers observed a similar set of satellites. It was therefore helpful to define channel groups for each commonly observed satellite, with one channel in the group corresponding to the satellite as tracked by the first receiver, and the other channel corresponding to the satellite as tracked by the second. Tracking groups as a class, however, may be easily modified for other experimental purposes.

FIGURE 2. Left: an independent tracking channel (corresponding to one tracking channel object). Right: a channel group. Note that in the channel group, updates to the code and carrier phase of each channel may be performed cooperatively.

Independent tracking channels have an update function that processes the next segment of raw data in three main steps: computing correlations (early, late, and prompt), producing discriminator outputs, and generating code and carrier-frequency updates. For a group of channels, this sequence of steps is interrupted after discriminator outputs have been computed. At this point, the channel group may instruct the tracking channels to update their code and carrier frequencies independently or through some other cooperative means that considers data across all of the channels.

As for the last few classes: correlators and filters are defined as objects so that they can be easily changed depending on the experimental circumstances. And satellites, in addition to holding satellite-specific ephemerides, have built-in functionality to return their locations given a particular epoch of GPS Time.

Naturally, core functions such as these would be found in traditional software receivers, but by repackaging them into the object-oriented framework, both code reusability and modifiability increase. And in addition, by defining classes for networks of receivers and groups of tracking channels, simulations and experiments involving cooperative positioning of receivers become easier to conduct.

Experiment

To help illustrate how the Python Receiver lends itself to the task of cooperatively tracking multiple receivers, concurrent data from two SDR front ends was collected on a boat in Lake Titicaca just offshore from Puno, Peru. The boat was a small motorized ferry capable of transporting approximately twenty passengers. One antenna and front end, hereafter referred to as “Receiver X” was placed on the port side of the boat, while the other, “Receiver Y” was placed on the starboard side. Maintaining a fixed baseline, both receivers captured raw GPS L1 signals from separate portions of the sky and mixed them to an intermediate frequency of 5.456 MHz. Raw data collection was performed concurrently at both receivers for 15 minutes as the boat returned from the floating islands of the Uros people to the dock at Puno. Finally, while Lake Titicaca is at a high elevation in the Altiplano (the Andean Plateau), the surrounding mountains do not rise far above the horizon, and thus visibility was quite good in most directions.

Some challenges, however, present themselves in this data set. While Receiver X was able to acquire eight satellites, and Receiver Y was able to acquire 10, the signal quality at Receiver Y was generally poor. In Figure 3, in-phase prompt correlator outputs from traditional scalar tracking are shown for both Receivers X and Y and satellites with PRN codes 27 and 29. For satellite 27, Receiver Y loses lock of the signal between code periods 100,000 and 200,000, and for satellite 29, it completely loses track of the signal after only a few thousand code periods. (Recall that the C/A-code period is one millisecond.)

FIGURE 3. The in-phase prompt correlator outputs for both receivers and satellites PRN 27 and 29. The cyan dots are correlator outputs, the red line is the locking metric, and the dashed green and blue lines are the thresholds set for determining good and poor lock, respectively. Locking metric values above the dashed green line represent a good lock, and values below the dashed blue line represent loss-of-lock. Note that y-axis values differ from graph to graph.

To better characterize the tracking performance of each receiver-satellite pair, a locking metric was designed and implemented, the values of which are shown as the red lines in the graphs of Figure 3. Inspired by the earlier use of the square-law detector, we have expressed the metric as:

(1)

where N is the number of most recent correlator samples, I_i and Q_i are the ith in-phase and quadrature-phase prompt correlator outputs, and the square-root operator returns the negative square root of the absolute value of the expression under the radical if that expression is negative.

After visually examining the relationship of this locking metric with the quality of the in-phase prompt correlator outputs, two thresholds were determined in order to better characterize the quality of the tracking loop lock. The first threshold, represented as the dashed green lines in the graphs of Figure 3, is the threshold above which the tracking loops were considered locked well. Its value was set to 250. The second threshold, whose value was set to 150 and is represented by the dashed blue lines, is the threshold below which the tracking loops were considered to be in a complete loss-of-lock situation. Locking metric values between 150 and 250 were considered as representing a situation in which the tracking loops were weakly locked to the incoming signals.

Despite the poor performance of Receiver Y in tracking many of its signals, navigation functionality in the Python Receiver was still able to recover sufficient ephemerides from the tracking data to perform position calculations. FIGURE 4 shows the navigation solutions for Receiver Y over a 13-minute interval, roughly capturing the route that the ferry took westward back to Puno. Note that the moustache-shaped region in the right-hand side of the map is the collection of floating islands of the Uros. Just as the ferry left these islands, the navigation solutions for Receiver Y become much nosier. Possible reasons for this are the slight change in heading that the ferry made, or the thicket of reeds that surrounded the boat during this portion of the journey. Navigation results for Receiver X were much less noisy.

FIGURE 4. The trip back to Puno on the left (west) from the floating islands of the Uros on the right (east) as determined by traditional scalar tracking and navigation at Receiver Y. Image courtesy of Google Earth and the GPS Visualizer.

Cooperative Scalar Tracking

While all of these traditional results were obtained using the Python Receiver, they could have just as easily been obtained using procedurally coded receivers. Assuming, however, that one is interested in performing experiments that involve data sharing between multiple receivers, the Python Receiver lends itself handily to the task.

An experiment was devised in which scalar tracking performed at both Receivers X and Y would be done cooperatively. In particular, it was observed that often when one of the two receivers momentarily lost track of its signal for a particular satellite, the other receiver would be tracking well. In addition, it was noted that because the two receivers maintained a fixed baseline during tracking, their tracking channels should have maintained a steady difference in code phases that changed slowly provided that the receiver-satellite geometry did not change quickly. As shown in FIGURE 5, the only violation of this scenario would occur when one of the two receivers lost lock and thus allowed for drift in its code-tracking loop. It should be noted that unlike the situation in Figure 5, the reported code difference between the two receivers suffered from a bias that grew linearly in time. This bias, which was likely due to clock errors in one or both of the receiver front ends, was eliminated through a linear regression before the plotting of the figure.

FIGURE 5. The code-phase difference between Receivers X and Y for PRN 27 from 300,000 to 500,000 milliseconds. Note the large variance around 400,000 milliseconds corresponding to a loss-of-lock for Receiver Y.

All of these observations motivated the following cooperative scalar tracking design. First, any satellite that was observed by only one receiver would be independently tracked by that receiver in the traditional manner. A single tracking loop object would be allocated in Python for this particular receiver-satellite pair. Second, any satellite that was observed by both receivers would have a channel group object allocated in Python. This channel group would contain two tracking channel objects, one for each receiver.

As shown in Figure 2, this channel group required specific code to be written to handle the cooperative updates of both receivers’ code and carrier frequencies. The algorithm was designed as follows. For each update epoch (generated by a call of the channel group’s update function), if both of the tracking channels were locked to their incoming signals, the channel group would save their code-phase difference for that code period. And since both channels were locked, both would update their code and carrier frequencies in the traditional manner, relying on discriminator outputs only.

If, on the other hand, one of the tracking channels was in a loss-of-lock situation, the channel group would search the previous 5,000 milliseconds of data for code periods during which, presumably, both tracking channels were mutually locked. This data would contain information about the expected code-phase difference between the two tracking channels at the current code period. At this point, a linear regression on the data from the mutually locked code periods was used to determine this expected code-phase difference. Finally, we note again that this expected code-phase difference would only remain valid under the assumption that the receiver-satellite geometry was not changing rapidly, as was the case for this data. But acknowledging that some changes in the geometry might occur (such as a change in heading of the boat) is the reason why the search interval for mutually locked data was limited to five seconds.

Assuming that one of the receivers was in a loss-of-lock situation and that sufficient data from the past five seconds existed to generate an estimate of the current expected code-phase difference, the channel group could then make a cooperative update of the lockless tracking channel. For this channel, the channel group would replace the traditional code-tracking discriminator outputs with the offset of the expected code-phase difference d_exp from the currently observed code-phase difference d_cur. In the following equation, the new discriminator output is denoted as c:

. (2)

Expressing d_cur=y_cur−x_cur and d_exp=y_exp−x_exp, where x_cur/exp and y_cur/exp represent current and expected code phases at two receivers, we can rewrite Equation 2 as

(3)

or

(4)

since we expect the x receiver to be locked, and therefore .

Some finer points to mention include that the “loss-of-lock” and “tracking well” designations were determined by way of the locking metric defined in the previous section. In addition, if a receiver was “tracking weakly,” it would update its code and carrier frequencies by relying solely on its own discriminator outputs. Also, because in traditional scalar tracking loss-of-lock might occur for an extended interval greater than five seconds at one receiver (such as Receiver Y’s tracking of satellite 27 seen in Figure 3 between 300,000 and 400,000 milliseconds), whenever the channel group was called to cooperatively update a lockless tracking channel’s code frequency, it would record the current code-phase difference between both receivers. Under all scenarios, the carrier-frequency update would be done independently at each channel using discriminator outputs alone. And finally, in order for both receivers to share relevant data with each other during tracking, clock bias terms found after traditional scalar tracking were used to align in time the raw data files for each receiver appropriately.

Results and Discussion

Using cooperative scalar tracking, drifting of the code-phase difference during code periods when one of the receivers is experiencing loss-of-lock is expected to be suppressed. And indeed, results such as those shown in FIGURE 6 verify this expectation. Since cooperative scalar tracking does not attempt to modify the way either receiver tracks during periods of good lock, this type of modified scalar tracking is not expected to produce less noisy tracking results. It is expected, however, to help lockless tracking channels to regain track after short signal outages, similar to the benefits of vector tracking.

FIGURE 6. The code-phase difference between Receivers X and Y for PRN 27 from 300,000 to 500,000 milliseconds, this time using cooperative scalar tracking. Presence of the red line indicates code periods during which cooperative code-phase updates were made for Receiver Y. Note that noisy drifting of the code-phase difference is suppressed.

Strikingly, this form of cooperative tracking allowed for Receiver Y to continually track the signal from satellite 29 (albeit with occasional outages) for the full thirteen minutes of data shown in FIGURE 7. Whereas in Figure 3, Receiver Y very quickly loses track of satellite 29, Figure 7 shows that Receiver Y, under cooperative scalar tracking, can maintain a good enough lock on the signal that by roughly 750,000 code periods, it is able to pick up the signal again quite strongly. This change in signal strength may have been due to a slight change in heading that the ferry made near Isla Taquile towards the end of this data set (see Figure 4 and FIGURE 8).

FIGURE 7. The in-phase prompt outputs for Receiver Y and PRN 29 using cooperative scalar tracking. Compare this to the bottom-right graph in Figure 3. Inter-receiver aiding allowed Receiver Y to track this signal for a majority of the code periods.

FIGURE 8. The trip back to Puno as determined by Receiver Y after cooperative scalar tracking and navigation computations. Compared to Figure 4, the navigation solutions are less noisy. Image courtesy of Google Earth and the GPS Visualizer.

Given the locking metric defined in the section “Experiment,” quantitative measures of how often each channel spent locked or in loss-of-lock can be made. In total, both receivers tracked six common satellites (with each receiver also tracking other satellites independently). TABLE 1 shows the locking frequencies for each commonly tracked satellite.

TABLE 1. Percent of time each tracking channel spent locked. Lock was designated if the locking metric was above 150. The best values for Receiver Y are highlighted in green, with the most notable improvement occurring for satellite 29.

Granted that the drift in the code phase for lockless tracking channels is curtailed in cooperative scalar tracking, an improvement in navigation solutions is also expected. This expectation is verified by comparing the qualitative level of noise in the solutions of Figure 8 to the solutions in Figure 4. Notably, the noise in the reed thicket (the section of the route immediately after leaving the moustache-shaped floating islands region) is suppressed. Not shown are the navigation solutions for the port side receiver, Receiver X, which by comparison to Receiver Y were relatively good in both forms of scalar tracking.

Conclusion

The experiment we carried out highlighted the abilities of the Python Receiver. Data from two SDR front ends and associated antennas placed on either side of a small transport ferry was used to track both receivers by using groups of tracking channels that could cooperatively modify their individual channels’ code and carrier frequencies. In this way, loss-of-lock in many of the tracking channels was avoided leading to improved navigation precision. More importantly, it is expected that future experiments like these can be easily implemented within the framework of the Python Receiver, and thus topics like cooperative vector tracking might be more easily investigated.

Acknowledgments

This article is based, in part, on the paper “A Python Software Platform for Cooperatively Tracking Multiple GPS Receivers” presented at ION GNSS+ 2014, the 27th International Technical Meeting of the Satellite Division of The Institute of Navigation, held in Tampa, Florida, September 8–12, 2014.

Manufacturers

The Python Receiver uses SiGe GN3S v3 Samplers, developed by the University of Colorado and SiGe Semiconductor (acquired by Skyworks Solutions Inc., Woburn, Massachusetts) and marketed by SparkFun Electronics, Niwot, Colorado.

ELIOT WYCOFF received his B.S. in applied mathematics from Columbia University, New York, in 2011. While working on the Python Receiver, he was a graduate student in the Department of Aerospace Engineering at the University of Illinois at Urbana-Champaign (UIUC).

YUTING NG obtained a B.S. in electrical and computer engineering from UIUC in 2014. She is currently a graduate student in the Department of Aerospace Engineering, UIUC.

GRACE XINGXIN GAO is an assistant professor in the Department of Aerospace Engineering, UIUC. She received her B.S. in mechanical engineering in 2001 and her M.S. in electrical engineering in 2003, both from Tsinghua University, China. She obtained her Ph.D. in electrical engineering at Stanford University in 2008. Before joining UIUC in 2012, Gao was a research associate at Stanford University.

FURTHER READING

• Authors’ Conference Paper

“A Python Software Platform for Cooperatively Tracking Multiple GPS Receivers” by E. Wycoff and G.X. Gao in Proceedings of ION GNSS+ 2014, the 27th International Technical Meeting of the Satellite Division of The Institute of Navigation, Tampa, Florida, September 8–12, 2014, pp. 1417–1425.

• Software-Defined GNSS Receivers

Digital Satellite Navigation and Geophysics: A Practical Guide with GNSS Signal Simulator and Receiver Laboratory by I.G. Petrovski and T. Tsujii with foreword by R.B. Langley, published by Cambridge University Press, Cambridge, U.K., 2012.

“Software GNSS Receiver: An Answer for Precise Positioning Research” by T. Pany, N. Falk, B. Riedl, T. Hartmann, G. Stangl, and C. Stöber in GPS World, Vol. 23, No. 9, September 2012, pp. 60–66.

A Software-Defined GPS and Galileo Receiver: A Single-Frequency Approach by K. Borre, D.M. Akos, N. Bertelsen, P. Rinder, and S.H. Jensen, published by Birkhäuser Engineering, Springer-Verlag GmbH, Heidelberg, 2007.

“GNSS Software Defined Radio: Real Receiver or Just a Tool for Experts?” by J.-H. Won, T. Pany, and G. Hein in Inside GNSS, Vol. 1, No. 5, July–August 2006, pp. 48–56.

“Satellite Navigation Evolution: The Software GNSS Receiver” by G. MacCougan, P.L. Normark, and C. Ståhlberg in GPS World, Vol. 16, No. 1, January 2005, pp. 48–55.

• Python

Learn Python in One Hour by V.R. Volkman, published by Modern Software Press, L.H. Press Inc., Ann Arbor, Michigan, 2014.

A Primer on Scientific Programming with Python by H.P. Langtangen, published by Springer-Verlag GmbH, Heidelberg, 2009.

“Python for Scientific Computing” by T.E. Oliphant in Computing in Science & Engineering, Vol. 9, No. 3, May–June 2007, pp. 10–20, doi: 10.1109/MCSE.2007.58.

• Noncoherent Integration

“GNSS Radio: A System Analysis and Algorithm Development Research Tool for PCs” by J.K. Ray, S.M. Deshpande, R.A. Nayak, and M.E. Cannon in GPS World, Vol. 17, No. 5, May 2006, pp. 51–56.

Fundamentals of Global Positioning System Receivers: A Software Approach, 2nd edition, by J. B.-Y. Tsui, published by Wiley-Interscience, John Wiley & Sons, Inc., Hoboken, New Jersey, 2005.

“An Assisted GPS Acquisition Method Using L2 Civil Signal in Weak Signal Environment” by D.J. Cho, C. Park, and S.J. Lee in Journal of Global Positioning Systems, Vol. 3 No. 1-2, December 2004, pp. 25–31.

• GPS Position Display

“GPS Visualizer: Do-It-Yourself Mapping” website by A. Schneider.

• Square Law Detector

“Lock Detection in Costas Loops” by A. Mileant and S. Hinedi in IEEE Transactions on Communications, Vol. 40, No. 3, March 1992, pp. 480–483, doi: 10.1109/26.135716.

February 1, 2015
Innovation: A Bright Idea

Testing the Feasibility of Positioning Using Ambient Light

By Jingbin Liu, Ruizhi Chen, Yuwei Chen, Jian Tang, and Juha Hyyppä

INNOVATION INSIGHTS by Richard Langley

AND THEN THERE WAS LIGHT. Well, the whole electromagnetic (EM) spectrum, actually. Visible light occupies only a small portion of the spectrum, which extends from below the extremely low frequency (ELF) 3 to 30 hertz band with equivalent wavelengths of 100,000 to 10,000 kilometers through infrared, visible, and ultraviolet light and x-rays to gamma rays in the 30 to 300 exahertz band (an exahertz is 10¹⁸hertz) with wavelengths of 10 to 1 picometers and beyond. The radio part of the spectrum extends to frequencies of about 300 gigahertz or so, but the distinction between millimeter radio waves and long infrared light waves is a little blurry.

Natural processes can generate electromagnetic radiation in virtually every part of the spectrum. For example, lightning produces ELF radio waves, and the black hole at the center of our Milky Way Galaxy produces gamma rays. And various mechanical processes can be used to generate and detect EM radiation for different purposes from ELF waves for communication tests with submerged submarines to gamma rays for diagnostic imaging in nuclear medicine.

Various parts of the EM spectrum have been used for navigation systems over the years. For example, the Omega system used eight powerful terrestrial beacons transmitting signals in the range of 10 to 14 kilohertz permitting global navigation on land, in the air, and at sea. At the other end of the spectrum, researchers have explored the feasibility of determining spacecraft time and position using x-rays generated by pulsars — rapidly rotating neutron stars that generate pulses of EM radiation.

But the oldest navigation aids, lighthouses, used the visible part of the EM spectrum. The first lighthouses were likely constructed by the ancient Greeks sometime before the third century B.C. The famous Pharos of Alexandria dates from that era. And before the construction of lighthouses, mariners used fires built on hilltops to help them navigate. The Greeks also navigated using the light from stars, or celestial navigation. Records go back to Homer’s Odyssey where we read “Calypso, the lovely goddess had told him to keep that constellation [the Great Bear] to port as he crossed the waters.” By around 1500 A.D., the astrolabe and the cross-staff had been developed sufficiently that they could be used to measure the altitudes of the sun or stars to determine latitude at sea. Celestial navigation was further advanced with the introduction of the quadrant and then the sextant. And determining longitude was possible by observing the moons of Jupiter (but not easily done at sea), measuring distances between the moon and other celestial bodies and, once it was developed, using a chronometer to time altitude observations.

How else is light used for positioning and navigation? Early in the space age, satellites were launched with flashing beacons or with large surface areas to reflect sunlight so that they could be photographed from the ground against background stars with known positions to determine the location of the camera. We also have laser ranging to satellites and the moon and the related terrestrial LiDAR technology, as well as the total stations used by surveyors. And in this month’s column, we take a look at the simple, innovative method of light fingerprinting: the use of observations of the artificial light emitted by unmodified light fixtures as well as the natural light that passes through windows and doorways in a technique for position determination inside buildings.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas.

Over the years, various localization technologies have been used to determine locations of people and devices in an absolute or relative sense. Relative positioning methods determine a location relative to another one in a local coordinate framework, while absolute positioning techniques fix an absolute location in a specific coordinate framework.

In the past, people observed the positions (orientation angles) of a celestial body (such as the sun, the moon, or a star) to determine their locations on the Earth, which is known as celestial navigation (see FIGURE 1). The locations are resolved by relating a measured angle between the celestial body and the visible horizon to the Nautical Almanac, which is a knowledge base containing the coordinates of navigational celestial bodies and other relevant data. Other than an observation device, celestial navigation does not rely on any infrastructure, and hence it can be used virtually anywhere on the globe at anytime, weather permitting. Nowadays, an increasing number of applications, location-based services, and ambient intelligence largely require positioning functions across various environments due to increasing mobility of people and devices. In particular, the development of robotics for a number of purposes requires the support of localization capability in various conditions where positioning infrastructure may be missing.

Various positioning technologies share an intrinsic characteristic that a positioning solution is resolved by using the dependency between spatial locations and a set of physical observables. The dependency may be expressed in the form of either a deterministic function model or a probabilistic model. A deterministic model expresses the dependency between locations and observables in a closed-form function, while a probabilistic model defines the dependency between locations and observables in the Bayesian sense. Depending on the form of dependency, different mathematical models have been used for position resolution.

For example, satellite-based GNSS positioning derives the location of a user’s receiver based on radio frequency (RF) signals transmitted by the satellite systems. GNSS positioning is grounded in accurate time determination: the time differences between the transmitted and the received radio signals denote signal travel times (observables), which are then converted into distance measurements between the satellite and the user antenna. Using the distance measurements between the user antenna and four different satellites, the receiver can obtain three-dimensional receiver coordinates in a global reference frame and the time difference between the receiver and satellite clocks. The dependency between user location and a set of distance observables can be expressed in a simplified equation:

(1)

where ρⁱ is an observed range between the ith satellite and the receiver, (x,y,z)ⁱ is the position of the ith satellite, (x,y,z) is the position of the receiver to be estimated, γ denotes errors in the range observable, δ_t and c are receiver clock error and the speed of light, respectively (the sign of the clock term is arbitrary, but must be used consistently).

It is obvious that GNSS positioning relies strongly on the visibility of the GNSS constellation — the space infrastructure — as it requires line-of-sight visibility of four or more satellites. The positioning capability is degraded or totally unavailable in signal-blocked environments, such as indoors and in urban canyons.

An example of Bayesian positioning is to use various signals of opportunity (SOOP) — signals not originally intended for positioning and navigation. They include RF signals, such as those of cellular telephone networks, digital television, frequency modulation broadcasting, wireless local area networks, and Bluetooth, as well as naturally occurring signals such as the Earth’s magnetic field and the polarized light from the sun. Indicators of these signals, such as signal strengths and signal quality, are dependent on locations in the Bayesian sense. The dependency between signal indicators and locations is expressed in a probabilistic model:

(2)

where signifies a dependency between a set of physical signals and locations, I denotes indicators of SOOP signals, L denotes location, and P(i|l) is the probability that signal indicators (i) are observed at location (l).

Positioning resolution involves finding a location that yields the maximum a posteriori probability given a specific set of observables. Bayes’ Rule for computing conditional probabilities is applicable in the positioning estimation, and a family of Bayesian inference methods has been developed (see Further Reading).

An inertial navigation system (INS) is a typical relative positioning technology, and it provides the estimation of moved distance, direction, and/or direction change. A commonly used INS consists of accelerometers, gyroscopes, and a compass. It is self-contained and needs no infrastructure in principle to operate. However, the sensors yield accumulated positioning errors, and they need extra information for calibration. For example, in a GNSS/INS combined system, the INS needs to be calibrated using GNSS positioning results. To achieve an enhanced positioning performance in terms of availability, accuracy, and reliability, different positioning technologies are commonly integrated to overcome the limitations of individual technologies in applicability and performance.

This article discusses the feasibility of ambient light (ambilight) positioning, and we believe it is the first time that ambilight has been proposed as a positioning signal source. We propose the use of two types of observables of ambient light, and correspondingly two different positioning principles are applied in the positioning resolution. Our solution does not require any modifications to commonly used sources of illumination, and it is therefore different from other indoor lighting positioning systems that have been proposed, which use a modulated lighting source.

Ambilight positioning does not require extra infrastructure because illumination infrastructure, including lamps and their power supply and windows, are always necessary for our normal functioning within spaces. Ambilight exists anywhere (indoor and outdoor), anytime, if we consider darkness as a special status of ambient light. Ambilight sensors have been sufficiently miniaturized and are commonly used. For example, an ambilight sensor is used in a modern smartphone to detect the light brightness of the environment and to adaptively adjust the backlight, which improves the user vision experience and conserves power. Additionally, ambilight sensors are also widely used in automotive systems to detect the light intensity of environments for safety reasons. Therefore, ambilight positioning can use existing sensors in mobile platforms. This article presents the possibilities and methods of ambilight positioning to resolve both absolute and relative positioning solutions, and which can be integrated as a component in a hybrid positioning system.

Absolute Positioning Using Ambilight Spectral Measurements

The essence of localization problems is to resolve the intrinsic dependency of location on a set of physical observables. Therefore, a straightforward idea is that the type of observables applicable to positioning can be determined once the location-observables dependency is established. The feasibility is validated when the location-observables dependency is confirmed in the sense of necessary and sufficient conditions.

Ambient light is a synthesis of artificial light sources and natural light. The light spectrum is defined by the distribution of lighting intensity over a particular wavelength range. Researchers have reported development of sensor technology that has a spectral response from 300 to 1450 nanometers (from ultraviolet through infrared light). The spectrum of ambient light is mainly determined by colors of reflective surfaces in the circumstance, in addition to that of artificial and natural light sources. Therefore, intensity spectrum measurements are strongly correlated with surrounding environments of different locations. The traditional fingerprinting method can be used to resolve the positioning solution.

The fingerprinting approach makes use of the physical dependency between observables and geo-locations to infer positions where signals are observed. This approach requires the knowledge of observable-location dependency, which comprises a knowledge database. The fingerprinting approach resolves the most likely position estimate by correlating observed SOOP measurements with the knowledge database. The related fingerprinting algorithms include K-nearest neighbors, maximum likelihood estimation, probabilistic inference, and pattern-recognition techniques. These algorithms commonly consider moving positions as a series of isolated points, and they are therefore related to the single-point positioning approach. In addition, a “hidden Markov” model method has been developed to fuse SOOP measurements and microelectromechanical systems (MEMS) sensors-derived motion-dynamics information to improve positioning accuracy and robustness.

In the case of ambilight positioning, prior knowledge is related to structure layout information, including the layout of a specific space, spatial distribution of lighting sources (lamps), types of lighting sources, and windows and doors where natural light can come through. Spatial distribution of lighting sources is normally set up together with power supplies when the structure is constructed, and their layout and locations are not usually changed thereafter. For example, illumination lamps are usually installed on a ceiling or a wall in fixed positions, and the locations of doors and windows, through which light comes, are also typically fixed throughout the life of a building. Therefore, the knowledge database of lighting conditions can be built up and maintained easily through the whole life cycle of a structure.

In practice, a specific working region is divided into discrete grids, and intensity spectrum measurements are collected at grid points to construct a knowledge database. The grid size is determined based on the required spatial resolution and spatial correlation of spectrum measurements. The spatial correlation defines the degree of cross-correlation of two sets of spectrum measurements observed at two separated locations.

We measured the spectrum of ambient light with a two-meter grid size in our library. The measurements were conducted using a handheld spectrometer. FIGURE 2 shows a set of samples of ambilight spectrum measurements, and the corresponding photos show the circumstances under which each spectrum plot was collected. These spectral measurements show strong geo-location dependency. Spectrum differences of different locations are sufficiently identifiable. TABLE 1 shows the cross-correlation coefficients of spectral measurements of different locations. The auto-correlation coefficients of spectral measurements of a specific location are very close to the theoretical peak value of unity, and the cross-correlation coefficients of spectra at different locations are significantly low. Therefore, the correlation coefficient is an efficient measure to match a spectrum observable with a geo-referred database of ambilight spectra.

FIGURE 2. Ambilight spectral measurements of nine locations in the library of the Finnish Geodetic Institute (arbitrary units). The photos below the spectrum plots show the circumstances under which the corresponding spectral measurements were collected.

TABLE 1. Correlation coefficient matrix of spectral measurements of different locations.

Relative Positioning Using Ambilight Intensity Measurements

Total ambilight intensity is an integrated measure of the light spectrum, and it represents the total irradiance of ambient light. In general, a lamp produces a certain amount of light, measured in lumens. This light falls on surfaces with a density that is measured in foot-candles or lux. A person looking at the scene sees different areas of his or her visual field in terms of levels of brightness, or luminance, measured in candelas per square meter. The ambilight intensity can be measured by a light detector resistor (LDR), and it is the output of an onboard 10-bit analog-to-digital converter (ADC) on an iRobot platform, which is the platform for a low-cost home-cleaning robot as shown in FIGURE 3.

FIGURE 3. The iRobot-based multi-sensor positioning platform, which is equipped with a light sensor and other versatile positioning sensors as marked in the figure.

We designed a simple current-to-voltage circuit based on an LDR and a 10-kilohm resistor, and the integrated analog voltage is input into the iRobot’s ADC with a 25-pin D-type socket, which is called the Cargo Bay Connector. FIGURES 4 and 6 show that the LDR sensor was not saturated during the test whenever we turned the corridor lamps on or off. Since the output of the light sensor was not calibrated with any standard light source, the raw ADC output rather than real values of physical light intensity was used in this study. During the test, the iRobot platform ran at a roughly constant speed of 25 centimeters per second, and the response time of the LDR was 50 milliseconds according to the sensor datasheet. The sampling rate of light intensity measurements was 5 Hz. Thus, the ADC could digitalize the input voltage in a timely fashion.

FIGURE 4. Total irradiance intensity measurements of ambient light in a closed space. The estimated lamp positions (magenta points) can be compared to the true lamp positions (green points).

FIGURE 6. Total irradiance intensity measurements of ambient light in the open corridor of the third floor.

We conducted the experiments with the iRobot platform in two corridors in the Finnish Geodetic Institute building. The robot was controlled to move along the corridors, and it collected measurements as it traveled. The two corridors represent two types of environment. The corridor of the first floor is a closed space where there is no natural light, and the corridor of the third floor has both natural light and artificial illuminating light. The illuminating fluorescent lamps are installed in the ceiling. In a specific environment, fluorescent lamps are usually installed at fixed locations, and their locations are not normally changed after installation. Therefore, the knowledge of lamp locations can be used for positioning.

Ambilight positioning is relatively simple in the first case where there is no natural light in the environment and all ambilight intensity comes from artificial light. Because the fluorescent lamps are separated by certain distances, the intensity measurements have a sine-like pattern with respect to the horizontal distance along the corridor. The sine-like pattern is a key indicator to be used for detecting the proximity of a lamp. As shown in Figures 4 and 6, raw measurements of ambilight intensity and smoothed intensity have a sine-like pattern. Because raw intensity measurements have low noise, either raw measurements or smoothed intensity can be used to detect the proximity of a lamp. Figure 4 also shows the results of detection and the comparison to the true lamp positions. There are four fluorescent lamps in this corridor test. The first three were detected successfully, and the estimated positions are close to true positions with a root-mean-square (RMS) error of 0.23 meters. The fourth lamp could not be detected because its light is blocked by a shelf placed in the corridor just below the lamp as shown in FIGURE 5. Figure 4 shows the sine-like intensity pattern of the fourth lamp did not occur due to the blockage.

FIGURE 5. The light of the fourth lamp in the corridor is blocked by shelves, and the corresponding sine-like light pattern does not appear.

On the third floor, the situation is more complicated because there is both natural light and incandescent lamps in the corridor. Natural light may come in from windows, which are located at multiple locations on the floor. In addition, the light spectrum in the corridor may be interfered with by light from office rooms around the floor. To recover the sine-like intensity pattern of the lamps, the intensity of the background light was measured when the incandescent lamps were turned off. Therefore, the calibrated intensity measurements of illuminating lamps can be calculated as follows:

(3)

where I_a is the intensity measurements of composite ambient light, I_b is the intensity measurements of background light, and I_c is the intensity measurements of the calibrated ambient light of the illuminating lamps.

Figure 6 shows the intensity measurements of composite ambient light, background light, and calibrated lamp light. In addition, the intensity measurements of calibrated lamp light are smoothed by an adaptive low-pass filter to mitigate noise and interference. The intensity measurements of smoothed lamp light were used to estimate the positions of the lamps according to the sine-like pattern. The estimated lamp positions were compared to the true lamp positions, and the errors are shown in FIGURE 7. The estimated lamp positions have a mean error of 0.03 meters and an RMS error of 0.79 meters. In addition, for the total of 15 lamps in the corridor, only one lamp failed to be detected (omission error rate = 1/15) and one lamp was detected twice (commission error rate = 1/15).

Discussion and Conclusion

Ambilight positioning needs no particular infrastructure, and therefore it does not have the problem of infrastructure availability, which many other positioning technologies have, limiting their applicability. For example, indoor positioning systems using Wi-Fi or Bluetooth could not work in emergency cases when the power supply of these devices is cut off. What ambilight positioning needs is just the knowledge of indoor structure and ambilight observables. The lighting conditions of an indoor structure can be reconstructed based on the knowledge of the layout structure whenever illuminating lamps are on or off. Thus, ambilight observables can be related to the layout structure to resolve positioning estimates as we showed in this article.

Besides indoor environments, the methods we have presented are also applicable in many other GNSS-denied environments, such as underground spaces and long tunnels. For example, the Channel Tunnel between England and France has a length of 50.5 kilometers, and position determination is still needed in this kind of environment. In such environments, there is usually no natural light, and the intensity of illuminating lamps has a clear sine-like pattern.

In particular, ambient light positioning is promising for robot applications when a robot is operated for tasks in a dangerous environment where there is no infrastructure for other technical systems such as Wi-Fi networks. Given the knowledge of the lighting infrastructure acquired from the construction layout design, the method of ambilight positioning can be used for robot localization and navigation. Our tests have shown also that the proposed ambilight positioning methods work well with both fluorescent lamps and incandescent lamps, as long as the light intensity sensor is not saturated.

A clear advantage of the technique is that the illuminating infrastructure and the structure layout of these environments are kept mostly unchanged during their life cycle, and the lighting knowledge can be constructed from the structure design. Hence, it is easy to acquire and maintain these knowledge bases. The hardware of ambient light sensors is low-cost and miniature in size, and the sensors can be easily integrated with other sensors and systems.

Although a spectrometer sensor is not currently able to be equipped with a mobile-phone device, the proposed ambilight positioning techniques can still be implemented with a modern mobile phone in several ways. For example, an economical way would be to form a multispectral camera using a selection of optical filters of selected bands or a miniature adjustable gradual optical filter. The spectral resolution then is defined by the bandwidth of the band-pass optical filters and the optical characteristics of the gradual optical filter. Other sensors, such as an acousto-optic tunable filter spectrometer and a MEMS-based Fabry-Pérot spectrometer, could also be used to measure the spectrum of ambilight in the near future. With such techniques, ambilight spectral measurements can be observed in an automated way and with higher temporal resolution.

Acknowledgments

The work described in this article was supported, in part, by the Finnish Centre of Excellence in Laser Scanning Research (CoE-LaSR), which is designated by the Academy of Finland as project 272195. This article is based on the authors’ paper “The Uses of Ambient Light for Ubiquitous Positioning” presented at PLANS 2014, the Institute of Electrical and Electronics Engineers / Institute of Navigation Position, Location and Navigation Symposium held in Monterey, California, May 5–8, 2014.

JINGBIN LIU is a senior fellow in the Department of Remote Sensing and Photogrammetry of the Finnish Geodetic Institute (FGI) in Helsinki. He is also a staff member of the Centre of Excellence in Laser Scanning Research of the Academy of Finland. Liu received his bachelor’s (2001), master’s (2004), and doctoral (2008) degrees in geodesy from Wuhan University, China. Liu has investigated positioning and geo-reference science and technology for more than ten years in both industrial and academic organizations.

RUIZHI CHEN holds an endowed chair and is a professor at the Conrad Blucher Institute for Surveying and Science, Texas A&M University in Corpus Christie. He was awarded a Ph.D. degree in geophysics, an M.Sc. degree in computer science, and a B.Sc. degree in surveying engineering. His research results, in the area of 3D smartphone navigation and location-based services, have been published twice as cover stories in GPS World. He was formerly an FGI staff member.

YUWEI CHEN is a research manager in the Department of Remote Sensing and Photogrammetry at FGI. His research interests include laser scanning, ubiquitous LiDAR mapping, hyperspectral LiDAR, seamless indoor/outdoor positioning, intelligent location algorithms for fusing multiple/emerging sensors, and satellite navigation.

JIAN TANG is an assistant professor at the GNSS Research Center, Wuhan University, China, and also a senior research scientist at FGI. He received his Ph.D. degree in remote sensing from Wuhan University in 2008 and focuses his research interests on indoor positioning and mapping.

JUHA HYYPPA is a professor and the head of the Department of Remote Sensing and Photogrammetry at FGI and also the director of the Centre of Excellence in Laser Scanning Research. His research is focused on laser scanning systems, their performance, and new applications, especially those related to mobile laser scanning and point-cloud processing.

FURTHER READING

• Authors’ Conference Paper

“The Uses of Ambient Light for Ubiquitous Positioning” by J. Liu, Y. Chen, A. Jaakkola, T. Hakala, J. Hyyppä, L. Chen, R. Chen, J. Tang, and H. Hyyppä in Proceedings of PLANS 2014, the Institute of Electrical and Electronics Engineers / Institute of Navigation Position, Location and Navigation Symposium, Monterey, California, May 5–8, 2014, pp. 102–108, doi: 10.1109/PLANS.2014. 6851363.

• Light Sensor Technology

“High-Detectivity Polymer Photodetectors with Spectral Response from 300 nm to 1450 nm” by X. Gong, M. Tong, Y. Xia, W. Cai, J.S. Moon, Y. Cao, G. Yu, C.-L. Shieh, B. Nilsson, and A.J. Heeger in Science, Vol. 325, No. 5948, September 25, 2009, pp. 1665–1667, doi: 10.1126/science.1176706.

• Light Measurement

“Light Intensity Measurement” by T. Kranjc in Proceedings of SPIE—The International Society for Optical Engineering (formerly Society of Photo-Optical Instrumentation Engineers), Vol. 6307, Unconventional Imaging II, 63070Q, September 7, 2006, doi:10.1117/12.681721.

• Modulated Light Positioning

“Towards a Practical Indoor Lighting Positioning System” by A. Arafa, R. Klukas, J.F. Holzman, and X. Jin in Proceedings of ION GNSS 2012, the 25th International Technical Meeting of the Satellite Division of The Institute of Navigation, Nashville, Tennessee, September 17–21, 2012, pp. 2450–2453.

• Application of Hidden Markov Model Method

“iParking: An Intelligent Indoor Location-Based Smartphone Parking Service” by J. Liu, R. Chen, Y. Chen, L. Pei, and L. Chen in Sensors, Vol. 12, No. 11, 2012, pp. 14612-14629, doi: 10.3390/s121114612.

• Application of Bayesian Inference

“A Hybrid Smartphone Indoor Positioning Solution for Mobile LBS” by J. Liu, R. Chen, L. Pei, R. Guinness, and H. Kuusniemi in Sensors, Vol. 12, No. 12, pp. 17208–17233, 2012, doi:10.3390/s121217208.

• Ubiquitous Positioning

“Getting Closer to Everywhere: Accurately Tracking Smartphones Indoors” by R. Faragher and R. Harle in GPS World, Vol. 24, No. 10, October 2013, pp. 43–49.

“Hybrid Positioning with Smartphones” by J. Liu in Ubiquitous Positioning and Mobile Location-Based Services in Smart Phones, edited by R. Chen, published by IGI Global, Hershey, Pennsylvania, 2012, pp. 159–194.

“Non-GPS Navigation for Security Personnel and First Responders” by L. Ojeda and J. Borenstein in Journal of Navigation, Vol. 60, No. 3, September 2007, pp. 391–407, doi: 10.1017/S0373463307004286.

November 4, 2014
Innovation: Scintillating Statistics
A Look at High-Latitude and Equatorial Ionospheric Disturbances of GPS Signals

By Yu Jiao, Yu (Jade) Morton, Steve Taylor, and Wouter Pelgrum

INNOVATION INSIGHTS by Richard Langley

THE EARTH’S IONOSPHERE. It’s both a blessing and a curse. Together with the magnetosphere, it helps to protect life on our planet from the damaging outpour of particle and electromagnetic radiation from the sun. In particular, it absorbs a lot of the extreme-ultraviolet (EUV) radiation arriving at the Earth. In fact, that is primarily how the ionosphere is formed. The EUV energy strips off the outer electrons of atmospheric gases producing a plasma of free electrons and ions.

The ionosphere has another beneficial role in that it permits long distance radio communication using high-frequency (HF) or shortwave signals. Although its use is in decline since the advent of the Internet, HF is still in use by some broadcasters and military organizations and is indispensible during natural disasters when electricity grids and network links go down.

But the ionosphere can be a pain, too, particularly for GNSS users. The signals from GNSS satellites must travel though the ionosphere on their way to receivers on or near the Earth’s surface. The signals are perturbed by the presence of the free electrons causing an advance in the phase of a signal’s carrier and a delay in the arrival of the pseudorandom noise code modulation (due to the refractive index being frequency dependent or dispersive) and so there is a contribution to carrier-phase and pseudorange (code) measurements, which must be accounted for when determining positions, velocities, and time (PVT) from the measurements.

Again, since the ionosphere is a dispersive medium, by linearly combining simultaneous measurements (either pseudoranges or carrier phases) on two frequencies such as the GPS L1 and L2 frequencies, an observable virtually free of ionospheric effects can be constructed and used for PVT determinations. This approach does require, however, a dual- or multi-frequency receiver. Single-frequency receivers (or the post-processing of single-frequency data) require the use of a model to account for the ionospheric biases as much as possible. The GPS navigation message, for example, includes values of the parameters of a simple ionospheric model. But, on average, its accuracy is only around 50%. More accurate ionospheric corrections can be acquired from elsewhere, even in real time, such as those from satellite-based augmentation systems.

But there is another ionospheric effect that can play havoc with GNSS signals: scintillations. These are rapid fluctuations in the amplitude and phase of the signals caused by small-scale irregularities in the ionosphere. When sufficiently strong, scintillations can result in the strength of a received signal dropping below the threshold required for acquisition and tracking or in causing problems for the receiver’s phase lock loop resulting in many cycle slips.

The occurrence of scintillations depends on many factors including solar and geomagnetic activity, time of year, time of day, and geographical location. In particular, scintillations are most prevalent in equatorial and polar (Arctic and Antarctic) regions. And the processes involved are not fully understood, hindering our ability to model and predict scintillations.

In an effort to help improve the monitoring, mapping, and modeling of scintillations, a team of researchers led by Prof. Jade Morton is monitoring high-latitude and equatorial scintillations and they discuss some of their preliminary results in this month’s column.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas. Write to him at lang @ unb.ca.

Among other effects of the Earth’s ionosphere on GPS and other GNSS signals, scintillation is potentially the most problematic. Ionospheric scintillation refers to the random amplitude and phase fluctuations of radio signals after propagating through plasma irregularities. These irregularities occur more frequently in high-latitude and equatorial regions, especially during solar maxima. Occurrence of scintillation is difficult to predict and model because of the complexity of the ionosphere’s internal mechanisms and solar activities that are the driving forces of space weather phenomena. GNSS signals are particularly vulnerable to scintillation, as strong scintillation can severely impact the acquisition and tracking processes in GNSS receivers, causing degradation in positioning accuracy and even loss-of-lock. With the increasing reliance on GNSS applications, understanding the characteristics of ionospheric scintillation and its effects on GNSS signals and receivers has become an important topic and has gained worldwide attention from both ionospheric scientists and GNSS engineers.

Since 2009, our research group has established several ionospheric scintillation monitoring and data collection systems located in high-latitude and equatorial regions. The results presented here are based on data collected from a specialized commercial dual-frequency GPS ionospheric monitoring receiver at Gakona, Alaska (62.4°N, 145.2°W), and a commercial multi-system, multi-frequency GNSS ionospheric monitoring receiver located at Jicamarca, Peru (11.9°S, 76.9°W).

Measurements are filtered to remove slowly varying trends caused by satellite-receiver dynamics, receiver oscillator errors, the background ionosphere and troposphere gradient, and other potential contributions from multipath and man-made interferences. Scintillation events above preset threshold levels from the filter outputs are extracted for analysis. The threshold levels are set based on two commonly used scintillation indices, the S₄ index and σ_φ, which are defined as the standard deviations of the detrended signal amplitude and carrier phase to represent the magnitude of signal intensity and phase fluctuation, respectively. In the study discussed in this article, the thresholds for S₄ and σ_φ are 0.15 and 15°, respectively for high-latitude measurements. For low-latitude data, the threshold for S₄ is raised to 0.2 to accommodate stronger amplitude scintillation, while the threshold for σ_φ remains 15°.

From data collected at Gakona, between August 2010 and March 2013, we extracted 655 amplitude and 2,355 phase-scintillation events from 657 equivalent days of data, while from data collected at Jicamarca, we extracted about 830 amplitude and 1,100 phase-scintillation events from 190 days of data collected from November 2012 to June 2013. Based on these events, we established a number of amplitude and phase scintillation distributions, which include scintillation-index-magnitude distributions, event-duration distributions, and event-occurrence frequency distributions. These results show very different characteristics of scintillation observed at low latitudes and high latitudes, indicating that there must be different mechanisms contributing to the formation and evolution of ionosphere plasma irregularities in the two regions. These characteristics are useful for scintillation-event prediction and modeling in the future.

Data Collection System and Event Thresholds

FIGURE 1 illustrates the general architecture of the event-driven GNSS data collection system. The system hardware consists of a multi-band GNSS antenna, a commercial ionospheric scintillation monitor (ISM) receiver, an array of reconfigurable software-defined radio (SDR) radio-frequency (RF) front-end devices capable of sampling intermediate-frequency (IF) signals, one or multiple data collection servers, a data storage array, timing signal distribution hardware to ensure both time and frequency consistency across all RF front ends and receivers, and network/communication devices that allow remote access of the receivers and servers to monitor the status of the hardware, to query recorded data, and reset and reconfigure the data collection system.

FIGURE 1. General architecture of the event-driven GNSS data collection system deployed at several high-latitude and equatorial sites since 2009.

Custom-designed space weather event monitoring and trigger software resides on the data collection and control server. The ISM receiver operates continuously to produce and record routine measurements such as I and Q channel accumulator outputs, pseudorange, carrier phase, Doppler frequency, C/N₀, and scintillation indices, while the SDR RF front ends only temporarily store the latest one-minute worth of IF samples in each device’s circular buffer. Scintillation event thresholds are pre-determined based on analysis of baseline data collected at the same local site using the same hardware. The real-time event trigger software compares ISM receiver measurements with the pre-set event threshold. If the measurements exceed the thresholds, the contents of the circular buffers will be written to the data storage array until after the event subsides. These raw IF samples are then further post-processed using a wide range of receiver processing algorithms for analysis of scintillation features and robust receiver algorithm development.

The high-latitude GNSS receiver array at Gakona, was initially established in 2009 and has been continuously evolving into a four-antenna array capable of collecting GPS L1, L2C, and L5 and GLONASS L1 and L2 signal data until its recent relocation to and upgrade at Poker Flat Research Range, north of Fairbanks. Several publications have discussed the system setup, receiver signal processing of data collected by the system, and characterization of high-latitude scintillations based on analysis of the array outputs (see Further Reading). In this article, only the data collected using the commercial ISM receiver are discussed because this is the longest operating receiver at this site. The receiver outputs L1C/A signal intensity and carrier-phase measurements at a rate of 50 Hz and semi-codeless tracking results of L2P(Y) at 1 Hz.

Since 2011, several GNSS data collection systems have been deployed at low-latitude locations, including Hong Kong, Singapore, Peru, Ascension Island, and Puerto Rico. In this article, we use results from the ISM receiver at Jicamarca, Peru, close to the geomagnetic equator. FIGURE 2 shows the data-collection-system-setup block diagram at Jicamarca. The ISM receiver used in this location generates 100-Hz carrier-phase measurements and I/Q channel correlator outputs; the latter are further processed to generate 50-Hz signal-intensity measurements for GPS L1C/A, L2C, and L5 signals and GLONASS, Galileo, and BeiDou open signals. Seven SDR front ends driven by the same oven-controlled crystal oscillator (OCXO) signal from the ISM receiver sample GPS, GLONASS, Galileo, and BeiDou open signals. Preliminary results obtained from these and other low-latitude SDR data have been presented in several papers in the archived literature (see Further Reading).

FIGURE 2. Current multi-GNSS data collection system configuration at Jicamarca Radio Observatory in Peru. (GLO = GLONASS, BDS = BeiDou System, VPN = virtual private network, ISMET = ionospheric scintillation monitoring event triggering, RAID = redundant array of independent disks)

The raw carrier-phase and signal-intensity measurements obtained from the two ISM receivers at Gakona and Jicamarca were detrended, from which the two scintillation indices S₄ and σ_φ were computed using Equations (1) and (2). In the two equations, I and φ stand for detrended signal intensity and carrier phase, respectively, and <·> represents the expected value that is essentially the average value over the interval of interest. In this study, the interval of interest was set to 10 seconds to most effectively highlight scintillation features based on evaluations of several different time intervals between 10 and 60 seconds.

(1)

(2)

As we mentioned earlier, the characterization of scintillation was carried out on the basis of scintillation events extracted from the raw data. After the evaluation of non-scintillation events and baseline indicators, a set of criteria has been established to extract interesting events through a semi-automated process from a large amount of data while keeping the number of selected events caused by non-scintillation factors (such as multipath and interference) low. A brief summary and explanations of the criteria are listed as follows:
- The elevation angle mask is 30° to reduce multipath effects.
- The thresholds for S₄ and σ_φ are 0.15 and 15° respectively for data collected at Gakona.
- For Jicamarca data, the thresholds are 0.2 and 15° respectively.
- To exclude interference cases, the index value has to remain above the threshold value for a minimum of 30 seconds to qualify as a scintillation event.
- An event detected within 5 minutes of the end of another event is combined as one event with the previous one.
- Scintillations experienced by multiple satellite signals simultaneously are treated separately, and events experienced simultaneously for all visible satellites are further analyzed to ensure that they are not caused by interferences.
- Carrier cycle slip/loss-of-lock detection and repair procedures are implemented to determine whether these cases are caused by scintillation or other factors.
It is important to note that the above criteria and procedures contain some degrees of arbitration, especially the last two, as they were applied based on visual inspections. These artificially imposed rules nevertheless are necessary for statistical analysis and comparison of scintillation observations.

Results and Discussion

In this section, we discuss the data sets we have collected and analyzed.

Available Dataset from Alaska and Peru. The ISM receiver at Gakona, started recording effective GPS data in August 2010. Environmental issues and human factors lead to a few intermittent data gaps during the more than three and a half years of data recording.

TABLE 1 lists monthly normal operation days and the percentage of time when data were collected. In all, the results presented in this article are based on approximately 3,000 scintillation events extracted from 657 days’ worth of data that was collected in a time span of 32 months.

Similarly, the number and percentage of days of effective data from Jicamarca, are summarized in Table 2. The dataset from this location runs from November 2012 until June 2013. Roughly 2,000 scintillation events have been extracted to enable statistical comparison of characteristics of scintillation observed in high- and low-latitude regions.

Scintillation Indicator Distributions. The magnitudes of the two scintillation indices, S₄ and σ_φ, are often used to indicate the intensity of ionospheric scintillation, as their values directly reflect the disturbance rate of received power and carrier-phase measurements. Although there have been discussions regarding the suitability of σ_φ as a phase scintillation indicator, it is, nevertheless, a measure of the magnitude of carrier variations in a certain spectral range that are related to scintillation activities. In the absence of a commonly accepted new indicator for phase scintillation, we will use σ_φ in this study simply as a means to measure the phase fluctuations. FIGURE 3 compares the intensity distributions of amplitude and phase scintillation observed at the Alaska (square markers) and Peru (triangle markers) sites. MaxS₄/σ_φ in the figures is the peak S₄ or σ_φ value during an amplitude or phase scintillation event, which is a more practical indicator of scintillation impact on GNSS receivers.

FIGURE 3. Maximum S4 and σφ distributions of (a) amplitude and (b) phase scintillation observed at Gakona, Alaska, and Jicamarca, Peru.

Figure 3a shows that amplitude scintillation events observed at Jicamarca are generally more intense than those observed at Gakona. This is consistent with most previous studies, which concluded that scintillation is the most intense in the equatorial region. Figure 3b, on the other hand, shows that the intensity of phase scintillation at Jicamarca is slightly lower than that at Gakona. Nevertheless, this result does not necessarily reflect scintillation intensity observed in other parts of the equatorial region, as Jicamarca is not located close to the equatorial anomaly crest where scintillation activity is the strongest.

The duration of a scintillation event is another indicator of scintillation’s negative impact on the acquisition and tracking processes in receivers. FIGURE 4 plots the amplitude and phase event duration probability distributions, with the mean event durations at each site shown in the plots. The results show that at Gakona (square markers), phase scintillation lasts much longer than amplitude scintillation. At Jicamarca (triangle markers), amplitude scintillation events last slightly longer than the phase ones on average, and both types have much longer durations than those at high latitudes.

FIGURE 4. Duration distributions of (a) amplitude and (b) phase scintillation events observed at Gakona, Alaska, and Jicamarca, Peru.

Ionospheric scintillation of combined high intensity and long duration is usually considered a big threat to signal processing in GNSS receivers. Unfortunately, these two aspects are often correlated, especially at low latitudes. Moderate correlation coefficient values have been observed between scintillation durations and the magnitudes of scintillation indicators at Jicamarca (FIGURE 5b). The correlations, however, are much smaller at Gakona (FIGURE 5a), especially for amplitude scintillation events. These results further confirm that scintillation is a more severe issue in the equatorial region.

FIGURE 5. Scintillation duration vs. intensity at (a) Gakona, Alaska, and (b) Jicamarca, Peru.

Scintillation Occurrence Frequency and Relating Factors. We define the scintillation occurrence frequency as the number of scintillation events recorded during a certain time interval, which can be an hour, a day, a month, a season, and so on. The occurrence frequency is an important indicator in scintillation monitoring and forecasting, as it helps to identify the periods when scintillation events are most likely to occur.

FIGURE 6 illustrates scintillation hourly occurrence probabilities at the two sites with respect to Coordinated Universal Time (UTC) (upper) and hours post sunset (lower). Also consistent with numerous previous research findings, scintillation at high latitudes was more frequent during nighttime than at other times. Scintillation observed at Jicamarca occurred more frequently at night as well, but was greatly concentrated between one and two hours post sunset and midnight. Statistics show that 98% of Jicamarca’s scintillation events were observed from one to six hours after local sunset.

FIGURE 6. Scintillation occurrence frequency with respect to UTC hours and hours after sunset at (a) Gakona, Alaska, and (b) Jicamarca, Peru.

As demonstrated in Figure 6, scintillation occurrence frequency is largely influenced by solar inputs, which are the main driving force in atmospheric ionization and ionospheric irregularity formation. Scintillation occurrence can also be affected by geomagnetic activities. FIGURE 7 shows how scintillation occurrence frequency was affected by solar activity and seasons. The four seasons are defined as: spring (SP) – March to May; summer (SU) — June to August; fall (FA) — September to November; and winter (WI) – December to February. The intensity of solar activity is indicated by the smoothed average sunspot numbers, which are marked as black dots in the plot.

FIGURE 7. Seasonal scintillation occurrence frequency and smoothed sunspot number.

Several phenomena can be observed in Figure 7. At Gakona, scintillation occurrence frequency is clearly influenced by solar activity. The occurrence frequency is also modulated by season, with equinoxes generally more active than adjacent solstices. In contrast to the half-a-year cycle at high latitudes, scintillation occurrence frequency at Jicamarca more closely follows a one-year cycle as described in previous research, and decreases largely in the summer.

Our analysis also shows that the level of geomagnetic field activity also directly impacts scintillation occurrence frequency. FIGURE 8 shows the correlations between scintillation daily occurrence frequencies and Ap index values at the two sites. Ap is a widely used index that linearly reflects the daily average level of global geomagnetic field activity. Ap can be converted to the conventional Kp index using a quasi-logarithmic conversion table. The result in Figure 8a was obtained using data collected during seven months at Gakona: March and November 2011; March, July, October, and November 2012; and March 2013. During these months, scintillation activity was generally high. Figure 8b was generated using all the data listed in Table 2. Clearly shown in the plots, scintillation occurrence frequency at high latitudes is strongly correlated with geomagnetic field activities, while at Jicamarca such correlations do not exist. This result also confirms many previous research findings.

FIGURE 8. Daily scintillation occurrence frequency with respect to Ap index value at (a) Gakona, Alaska, and (b) Jicamarca, Peru.

Summary and Conclusions

This article presented comparative work on ionospheric scintillation characterization using data collected at Gakona, Alaska, and Jicamarca, Peru, during the current solar maximum to investigate the different natures of scintillation at high latitude and in equatorial regions. Scintillation intensity, duration, and occurrence frequency distributions were analyzed to demonstrate the differences at the two locations.

Scintillation in the equatorial region is typically more severe with deeper and faster signal power fadings and longer durations. Also, low-latitude scintillation with stronger intensity usually lasts longer, which further contributes to its negative impact on receivers. At high latitudes, phase fluctuations overwhelmed amplitude scintillation by the number of occurrences and their duration.

Scintillation is more frequent during nighttime, and almost all low-latitude scintillation events occur within six hours after local sunset. The overall occurrence frequency of scintillation not only increases with high solar activity, but also follows certain seasonal patterns. In general, scintillation is more active around the equinoxes. Additionally, high-latitude scintillation is also closely correlated to geomagnetic field activity, while the relationship is not obvious in the equatorial region.

Lastly, we would like to point out that the results presented here are preliminary and may be restricted to local effects, especially at low latitudes. As more data become available from Jicamarca and other equatorial sites where SDR data collection systems ensure quality inputs during strong scintillation events, a more comprehensive analysis and comparison can be made to facilitate global scintillation monitoring, mapping, and modeling.

Acknowledgments

The data collection and analysis project discussed in this article was supported by the U.S. Air Force Office of Scientific Research and Air Force Research Laboratory grants. The authors appreciate the support of High Frequency Active Auroral Research Program (HAARP) staff and the University of Alaska Fairbanks Geophysical Institute for organizing and sponsoring the HAARP campaign and HAARP staff support of the GNSS receiver data collection system setup. The authors would also like to acknowledge Jicamarca Radio Observatory for hosting the GNSS equipment. The Jicamarca Radio Observatory is a facility of the Instituto Geofisico del Peru, operated with support from the U.S. National Science Foundation through Cornell University. This article is based, in part, on the paper “Comparative Studies of High-latitude and Equatorial Ionospheric Scintillation Characteristics of GPS Signals” presented at PLANS 2014, the Institute of Electrical and Electronics Engineers / Institute of Navigation Position, Location and Navigation Symposium held in Monterey, California, May 5–8, 2014.

Manufacturers

The commercial ISM receivers used at Gakona and Jicamarca were a GPS Silicon Valley — now NovAtel Inc. — GSV4004B and a Septentrio N.V. PolaRxS Pro, respectively.

YU JIAO is a Ph.D. candidate at the Colorado State University (CSU), Fort Collins, Colorado. She received her master’s degree in computational science and engineering from Miami University, Oxford, Ohio, in 2013 and her bachelor’s degree in electronic and information engineering from Beihang University (previously known as the Beijing University of Aeronautics and Astronautics), Beijing, China, in 2011. Her research interests are in GNSS signal processing and ionosphere effects on GNSS in both high-latitude and equatorial regions.

YU (JADE) MORTON is an electrical engineering professor at CSU. She received a Ph.D. in electrical engineering from Pennsylvania State University (Penn State), State College, Pennsylvania, and was a post-doctoral research fellow in the Space Physics Research Laboratory of the University of Michigan, Ann Arbor, Michigan. Prior to joining CSU, she was a professor in the Department of Electrical and Computer Engineering at Miami University. Her research interests are advanced GNSS receiver algorithms for accurate and reliable operations in challenging environments, studies of the atmosphere using radar and satellite signals, and development of new applications using satellite navigation technologies.

STEVE TAYLOR is a graduate student in the Department of Electrical and Computer Engineering at Miami University. He received his B.S. in computer science from Miami University in 2011. Taylor developed software systems for ionosphere space weather monitoring and has been involved in deployment of Dr. Morton’s research team’s GNSS data collection system in Alaska, Peru, Hong Kong, Ascension Island, and Puerto Rico.

WOUTER PELGRUM is an assistant professor of electrical engineering at Ohio University, where he conducts research in and teaches about topics in electronic navigation, such as GNSS, Distance Measuring Equipment or DME, and time and frequency transfer. Before joining Ohio University in 2009, he worked in private industry, where he contributed to the development of an integrated GPS-eLoran receiver and antenna. From 2006 until 2008 he operated his own company, specializing in navigation-related research and consulting.

FURTHER READING

• Authors’ Conference Paper

“Comparative Studies of High-latitude and Equatorial Ionospheric Scintillation Characteristics of GPS Signals” by Y. Jiao, Y. Morton, and S. Taylor in Proceedings of PLANS 2014, the Institute of Electrical and Electronics Engineers / Institute of Navigation Position, Location and Navigation Symposium, Monterey, California, May 5–8, 2014, pp. 37–42, doi: 10.1109/PLANS.2014.6851355.

• Introduction to Ionospheric Scintillation and GNSS

“Ionospheric Scintillations: How Irregularities in Electron Density Perturb Satellite Navigation Systems” by the Satellite-Based Augmentation Systems Ionospheric Working Group in GPS World, Vol. 23, No. 4, April 2012, pp. 44–50.

“GNSS and Ionospheric Scintillation: How to Survive the Next Solar Maximum” by P. Kintner, Jr., T. Humphreys, and J. Hinks in Inside GNSS, Vol. 4, No. 4, July/August 2009, pp. 22–30.

“GPS and Ionospheric Scintillations” by P. Kintner, B. Ledvina, and E. de Paula in Space Weather, Vol. 5, S09003, 2007, doi: 10.1029/2006SW000260.

A Beginner’s Guide to Space Weather and GPS by P. Kintner, Jr., unpublished article, October 31, 2006.

“Limitations in GPS Receiver Tracking Performance Under Ionospheric Scintillation Conditions” by S. Skone, K. Knudsen, and M. de Jong in Physics and Chemistry of the Earth, Part A: Solid Earth and Geodesy, Vol. 26, No. 6-8, 2001, pp. 613–621, doi: 10.1016/S1464-1895(01)00110-7.

“Radio Wave Scintillations in the Ionosphere” — a review paper by C.K. Yeh and C.-H. Liu in Proceedings of the IEEE, Vol. 70, No. 4, 1982, pp. 324–360, doi: 10.1109/PROC.1982.12313.

High-Latitude Scintillations

“Characterization of High Latitude Ionospheric Scintillation of GPS Signals” by Y. Jiao, Y. Morton, S. Taylor, and W. Pelgrum in Radio Science, Vol. 48, 2013, pp. 698–708, doi: 10.1002/2013RS005259.

Equatorial Scintillations

“Statistics of GPS Scintillations over South America at Three Levels of Solar Activity” by A.O. Akala, P.H. Doherty, C.E. Valladares, C.S. Carrano, and R. Sheehan in Radio Science, Vol. 46, No. 5, October 2011, doi: 10.1029/2011RS004678.

“Measuring Ionospheric Scintillation in the Equatorial Region over Africa, Including Measurements from SBAS Geostationary Satellite Signals” by A.J. Van Dierendonck and B. Arbesser-Rastburg in Proceedings of ION GNSS 2004, the 17th International Technical Meeting of the Satellite Division of The Institute of Navigation, Long Beach, California, September 21–24, 2004, pp. 316–324.

“Effects of the Equatorial Ionosphere on GPS” by L. Wanninger in GPS World, Vol. 4, No. 7, July 1993, pp. 48–54.

Scintillation-Triggering Data Collection

“An Improved Ionosphere Scintillation Event Detection and Automatic Trigger for GNSS Data Collection Systems” by S. Taylor, Y. Morton, Y. Jiao, J. Triplett, and W. Pelgrum in Proceedings of ION ITM 2012, The Institute of Navigation 2012 International Technical Meeting, Newport Beach, California, January 30 – February 1, 2012, pp. 1563–1569.

Software Defined Radio Processing of GPS Scintillation Data

“Triple Frequency GPS Signal Tracking During Strong Ionospheric Scintillations over Ascension Island” by M. Carroll, Y.J. Morton, and E. Vinande in Proceedings of PLANS 2014, the Institute of Electrical and Electronics Engineers / Institute of Navigation Position, Location and Navigation Symposium, Monterey, California, May 5–8, 2014, pp. 43–49, doi: 10.1109/PLANS.2014.6851356.

Forecasting Scintillations

“A Forecasting Ionospheric Real-time Scintillation Tool (FIRST)” by R.J. Redmon, D. Anderson, R. Caton, and T. Bullett in Space Weather, Vol. 8, No. 12, December 2010, doi: 10.1029/2010SW000582.

“Specification and Forecasting of Scintillations in Communication/Navigation Links: Current Status and Future Plans” by S. Basu, K.M. Groves, Su. Basu, and P.J. Sultan in Journal of Atmospheric and Solar-Terrestrial Physics, Vol. 64, 2002, pp. 1745–1754, doi: 10.1016/S1364-6826(02)00124-4.

Alternative Scintillation Indices

“Improved Amplitude- and Phase-scintillation Indices Derived from Wavelet Detrended High-latitude GPS Data” by S.C. Mushini, P.T. Jayachandran, R.B. Langley, J.W. MacDougall, and D. Pokhotelov in GPS Solutions, Vol. 16, No. 3, July 2012, pp. 363–373, doi: 10.1007/s10291-011-0238-4

“Perils of the GPS Phase Scintillation Index (s_f)” by T.L. Beach in Radio Science, Vol. 41, RS5S31, 2006, doi: 10.1029/2005RS003356.

“Problems in Data Treatment for Ionospheric Scintillation Measurements” by B. Forte and S.M. Radicella in Radio Science, Vol. 37, No. 6, 1096, 2002, pp. 8-1–8.5, doi: 10.1029/2001RS002508.
October 2, 2014
Innovation: How Deep Is That White Stuff?

Using GPS Multipath for Snow-Depth Estimation

By Felipe G. Nievinski and Kristine M. Larson

INNOVATION INSIGHTS by Richard Langley

FRINGES. No, I’m not talking about the latest celebrity hairstyles nor the canopy of an American doorless, four-wheeled carriage from yesteryear (think Oklahoma!). I’m talking about interference fringes. But there is a connection to these other uses of the word fringe as we’ll see. You’ve all seen interference fringes at your local gas station, typically after it has just rained. They are the alternating bands of color we perceive when looking at a gasoline or oil slick in a puddle of water. They are caused by the white light from the Sun or artificial lighting reflected from the top surface of the slick and that from the bottom surface at the slick-water interface combining or interfering with each other at our eyeballs. The two sets of light waves arrive slightly out of phase with each other, and depending on the wavelengths of the reflected light and our angle of view, produce the colorful fringes. If the incident light was monochromatic, consisting of a single frequency or wavelength, then we would perceive just alternating bright and dark bands. The bright bands result from constructive interference when the phase difference is a near a multiple of 2π whereas the dark bands result from destructive interference when the difference is near an odd multiple of π.

Interference fringes had been seen long before the invention of the automobile. They are clearly seen on soap bubbles and the iridescent colors of peacock feathers, Morpho butterflies, and jewel beetles are also due to the interference phenomenon rather than pigmentation. Sir Isaac Newton did experiments on interference fringes (amongst other things) and tried to explain their existence — wrongly, it turned out. But he did coin the term fringes since they resembled the decorative fringe sometimes used on clothing, drapery, and, yes, surrey canopies.

It was the English polymath, Thomas Young, who, in 1801, first demonstrated interference as a consequence of the wave-nature of light with his famous double-slit experiment. You may have replicated his experiment in a high-school physics class. I did and I think I did it again as an undergraduate student taking a course in optics. Already by that point I was aiming for a career in physics or space science but I didn’t know that as a graduate student I would do research involving interference fringes. But not using light waves.

My research involved the application of very long baseline interferometry or VLBI to geodesy. VLBI had been developed by radio astronomers to better understand the structure of quasars and other esoteric celestial objects. At either ends of a baseline connecting large radio telescopes, perhaps stretching between continents, the quasar signals were recorded on magnetic tape and precisely registered using atomic clocks. When the tapes were played back and the signals aligned, one obtained interference fringes as peaks and troughs in an analog or digital waveform. Computer analysis of these fringes not only provided information on the structure of the observed radio source but also on the distance between the radio telescopes — eventually accurate enough to measure continental drift.

But what has all of this got to do with GPS? In this month’s column, we look at a technique that uses fringes generated by signals arriving at an antenna directly from GPS satellites and those reflected by snow surrounding the antenna to measure its depth and how it varies over time. GPS for measuring snow depth; who would have thought?

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas.

Snowpacks are a vital resource for human existence on our planet. They provide reservoirs of fresh water, storing solid precipitation and delaying runoff. One sixth of the world population depends on this resource. Both scientists and water-supply managers need to know how much fresh water is stored in snowpack and how fast it is being released as a result of melting.

Snow monitoring from space is currently under investigation by both NASA and ESA. Greatly complementary to such spaceborne sensors are automated ground-based methods; the latter not only serve as essential independent validation and calibration for the former, but are also valuable for climate studies and flood/drought monitoring on their own. It is desirable for such estimates to be provided at an intermediary scale, between point-like in situ samples and wider area pixels.

In the last decade, GPS multipath reflectometry (GPS-MR), also known as GPS interferometric reflectometry and GPS interference-pattern technique, has been proposed for monitoring snow. This method tracks direct GPS signals, those that travel directly to an antenna, that have interfered with a coherently reflected signal, turning the GPS unit into an interferometer (see FIGURE 1). Its main variant is based on signal-to-noise ratio (SNR) measurements, although GPS-MR is also possible with carrier-phase and pseudorange observables. Data are collected at existing GPS base stations that employ commercial-off-the-shelf receivers and antennas in a conventional, antenna-upright setup. Other researchers have used a custom antenna and/or a dedicated setup, with the antenna tipped for enhanced multipath reception.

FIGURE 1. Standard geodetic receiver installation. The antenna is protected by a hemispherical radome. The monument (tripod structure) is ~ 2 meters above the ground. GPS satellites rise and set in ascending and descending sky tracks, multiple times per day. The specular reflection point migrates radially away from the receiver for decreasing satellite elevation angle. The total reflector height is made up of an a priori value and an unknown bias driven by the thickness of the snow layer.

In this article, we summarize the SNR-based GPS-MR technique as applied to snow sensing using geodetic instruments. This forward/inverse approach for GPS-MR is new in that it capitalizes on known information about the antenna response and the physics of surface scattering to aid in retrieving the unknown snow conditions in the site surroundings. It is a statistically rigorous retrieval algorithm, agreeing to first order with the simpler original methodology, which is retained here for the inversion bootstrapping. The first part of the article describes the retrieval algorithm, while the second part provides validation at a representative site over an extended period of time.

Physical Forward Model

SNR observations are formulated as SNR = P_s/P_n. In the denominator, we have the noise power, P_n, here taken as a constant, based on nominal values for the noise power spectral density and the noise bandwidth. The numerator is composite signal power:

. (1)

Its incoherent component is the sum of the respective direct and reflected powers (although direct incoherent power is negligible). In contrast, the coherent composite signal power follows from the complex sum of direct and reflection average voltages (not to be confused with the electromagnetic propagating fields, which neglect the receiving antenna response and also the receiver tracking process):

(2)

It is expressed in terms of the coherent direct and reflected powers, as well as the interferometric phase,

, (3)

which amounts to the reflection excess phase with respect to the direct signal.

We decompose observations, SNR = tSNR + dSNR, into a trend

(4)

over which interference fringes are superimposed:

. (5)

From now on, we neglect the incoherent power, which only impacts tSNR, not dSNR, and drop the coherent power superscript, for brevity.

The direct or line-of-sight power is formulated as

(6)

where is the direction-dependent right-hand circularly polarized (RHCP) power component incident on an isotropic antenna; the left-handed circularly polarized (LHCP) component is negligible. The direct antenna gain, , is obtained evaluating the antenna pattern in the satellite direction and with RHCP polarization.

The reflection power,

, (7)

is defined starting with the same incident isotropic power, , as in the direct power. It ends with a coherent power attenuation factor,

(8)

where θ is the angle of incidence (with respect to the surface normal), k = 2π/λ, is the wave number, and λ = 24.4 centimeters is the carrier wavelength for the civilian GPS signal on the L2 frequency (L2C). This polarization-independent factor accounts only for small-scale residual height above and below a large-scale trend surface. The former/latter results from high-/low-pass filtering the actual surface heights using the first Fresnel zone as a convolution kernel, roughly speaking. Small-scale roughness is parameterized in terms of an effective surface standard deviation s (in meters); its scattering response is modeled based on the theories of random surfaces, except that the theoretical ensemble average is replaced by a sensing spatial average. Large-scale deterministic undulations could be modeled, but their impact on snow depth is canceled to first-order by removing bare-ground reflector heights.

At the core of , we have coupled surface/antenna reflection coefficients, , producing respectively RHCP and LHCP fields (under the assumption of a RHCP incident field). These terms include antenna response power gain and phase patterns, evaluated in the reflection direction, and separately for each polarization. The surface response is represented by complex-valued Fresnel coefficients for cross- and same-sense circular polarization, respectively. The medium is assumed to be homogeneous (that is, a semi-infinite half-space). Material models provide the complex permittivity, which drives the Fresnel coefficients.

The interferometric phase reads:

.(9)

The first term accounts for the surface and antenna properties of the reflection, as above. The last one is the direct phase contribution, which amounts to only the RHCP antenna phase-center variation evaluated in the satellite direction. The majority of the components present in the direct RHCP phase (such as receiver and satellite clock states, the bulk of atmospheric propagation delays, and so on) are also present in the reflection phase, so they cancel out in forming the difference.

At the core of the interferometric phase, we have the geometric component, φ_I = kτ_i, the product of the wave number and the interferometric propagation delay. Assuming a locally horizontal surface, the latter is simply:

(10)

in terms of the satellite elevation angle, e, and an a priori reflector height, H_A. Snow depth will be measured in terms of changes in reflector height.

The physical forward model, based only on a priori information, can then be summarized as:

(11)

where interferometric power and phase are, respectively:

(12)

. (13)

In all of these terms the pseudorandom-noise-code modulation impressed on the carrier wave can be safely neglected, given the small interferometric delay and Doppler shift at grazing incidence, stationary surface/receiver conditions, and short antenna installations.

Parameterization of Unknowns

There are errors in the nominal values assumed for the physical parameters of the model (permittivity, surface roughness, reflector height, and so on). Ideally we would estimate separate corrections for each one, but unfortunately many are linearly dependent or nearly so. Because of this dependency, we have kept physical parameters fixed to their optimal a priori values, and have estimated a few biases. Each bias is an amalgamation of corrections for different physical effects. In a later stage, we rely on multiple independent bias estimates (such as for successive days) to try and separate the physical sources.

Each satellite track is inverted independently. A track is defined by partitioning the data by individual satellite and then into ascending and descending portions, splitting the period between the satellite’s rise and set at the near-zenith culmination. Each satellite track has a duration of ~1–2 hours. This configuration normally offers a sufficient range of elevation angles, unless the satellite reaches culmination too low in the sky (less than about 20°), in which case the track is discarded. In seeking a balance between under- and over-fitting, between an insufficient and an excessive number of parameters, we estimate the following vector of unknown parameters:

. (14)

FIGURE 2 shows the effect of the constant and linear biases on the SNR observations. Reflector height bias, H_B , changes the number of oscillations; phase shift, φ_B , displaces the oscillations along the horizontal axis; reflection power, , affects the depth of fades; zeroth-order noise power, , shifts the observations up or down as a whole; and first-order noise power, , tilts the SNR curve. A good parameterization yields observation sensitivity curves as unique as possible for each parameter.

FIGURE 2. Effect of each parameter on SNR observations; curves are displaced vertically (6 dB) for clarity.

The forward model, now including the biases, can be summarized as follows:

(15)

where the modified interferometric power and phase are given by:

, (16)

. (17)

The total reflector height, H = H_A – H_B (a priori value minus unknown bias), is to be interpreted as an effective value that best fits measurements, which includes snow and other components.

Bootstrapping Parameter Priors. Biases and SNR observations are involved non-linearly through the forward model. Therefore, there is the need for a preliminary global optimization, without which the subsequent final local optimization will not necessarily converge to the optimal solution.

SNR observations would trace out a perfect sinusoid curve in the case of an antenna with isotropic gain and spherical phase pattern, surrounded by a smooth, horizontal, and infinite surface (free of small-scale roughness, large-scale undulations, and edges), made of perfectly electrically conducting material, and illuminated by constant incident power. Thus, in such an idealized case, SNR could be described exactly by constant reflector height, phase shift, amplitude, and mean values.

As the measurement conditions become more complicated, the SNR data start to deviate from a pure sinusoid. Yet a polynomial/spectral decomposition is often adequate for bootstrapping purposes.

Statistical Inverse Model Formulation

Based on the preliminary values for the unknown parameters vector and other known (or assumed) values, we run the forward model to obtain simulated observations. We form pre-fit residuals comparing the model values to SNR measurements collected at varying satellite elevation angles (separately for each track). Residuals serve to retrieve parameter corrections, such that the sum of squared post-fit residuals is minimized. This non-linear least squares problem is solved iteratively using both a functional model and a stochastic model. The functional modeling includes a Jacobian matrix of partial derivatives, which represents the sensitivity of observations to parameter changes where the partial derivatives are defined element-wise. Instead of deriving analytical expressions, we evaluate them numerically, via finite differencing. The stochastic model specifies the uncertainty and correlation expected in the residuals. Their a priori covariance matrix modifies the objective function being minimized.

Directional Dependence

It is important to know at which elevation angles the parameter estimates are best determined. Here, we focus on the phase parameters instead of reflection power or noise power parameters.

We can utilize the estimated reflector height and phase shift to evaluate the full phase bias function over varying elevation angles. Similarly, we can extract the corresponding 2-by-2 portion of the parameters’ a posteriori covariance matrix, containing the uncertainty for reflector height and for phase shift, as well as their correlation, which is then propagated to obtain the full phase uncertainty (see FIGURE 3).

FIGURE 3. Uncertainty of full phase function, propagated from the uncertainty of reflector height and of phase shift, as well as their correlation.

The uncertainty attains a clear minimum versus elevation angle. The least-uncertainty elevation angle pinpoints the observation direction where reflector height and phase shift are best determined (in combined form, not individually). The azimuth and epoch coinciding with the peak elevation angle act as track tags, later used for clustering similar tracks and analyzing their time series of retrievals.

If we normalize phase uncertainty by its value at the peak elevation angle, then plot such sensing weights (between 0 and 1) versus the radial or horizontal distance to the center of the first Fresnel zone at each elevation angle, we obtain FIGURE 4. It can be interpreted as the reflection footprint, indicating the importance of varying distances, with a longer far tail and a shorter near tail (respectively regions beyond and closer than the peak distance). The implications for in situ data collection are clear: one should sample more intensely near the peak distance (about 15 meters) and less so in the immediate vicinity of the GPS antenna, tapering it off gradually away from the antenna. As a caveat, these conclusions are not necessarily valid for antenna setups other than the one considered here.

FIGURE 4. Reflection footprint in terms of a sensing weight (between 0 and 1) defined as the normalized reciprocal of full phase uncertainty, plotted versus the radial or horizontal distance from the receiving antenna to the center of the first Fresnel zone at each elevation angle; valid for an upright 2-meter-tall antenna; the receiving antenna is at zero radial distance.

Results

We now examine the snow-depth retrievals from the GPS multipath retrieval algorithm and assess both the precision and accuracy of the method. Multiple metrics have been developed to assess the quality of the results. The accuracy of the method has been evaluated by comparing with in situ data over a multi-year period. Three field sites were chosen to highlight different limitations in the method, both in terms of terrain and forest cover: grassland, alpine, and forested. We will look at the forested site in some detail.

Satellite Coverage and Track Clustering. All GPS-MR retrievals reported here are based on the newer GPS L2C signal. Of the approximately 30 GPS satellites in service, 8-10 L2C satellites were available between 2009 and 2012 (8, 9, and 10 satellites at the end of 2009, 2010, and 2011, respectively). Satellite observations were partitioned into ascending and descending portions, yielding approximately twenty unique tracks per day at a site with good sky visibility. GPS orbits are highly repeatable in azimuth, with deviations at the few-degree range over a year, translating into ~50-100-centimeter azimuthal displacement of the reflecting area (corresponding to the first Fresnel zone at 10°-15° elevation angle for a 2-meter high antenna). This repeatability permits clustering daily retrievals by azimuth. It also allows the simplification that estimated snow-free reflector heights are fairly consistent from day to day, facilitating the isolation of the varying snow depth during the snow-covered period.

For a given track, its revisit time is also repeatable, amounting to practically one sidereal day. The deficit in time relative to a calendar day results in the track time of the day receding ~4 minutes and 6 seconds every day. This slow but steady accumulation eventually makes the time of day return to its starting value after about one year. As all GPS satellites drift approximately at the same rate, the time between successive tracks remains nearly repeatable. Its reciprocal, the sampling rate, has a median equal to approximately one track per hour, with a low value of one track within two hours and a high of one track within 15 minutes; both extremes occur every day, with low-rate idle periods interspersed with high-rate bursts. The time of the day reduced to a fixed day (such as January 1, 2000) could also be used to cluster tracks. Neighboring clusters, which are close in azimuth and/or in reduced time of the day, are expected to be more comparable, as they sample similar conditions and are subject to similar errors.

Observations. FIGURE 5 shows several representative examples of SNR observations. A typical good fit between measured and modeled values is shown in Figure 5(a), corresponding to the beginning of the snow season. Generally the model/measurement fit is good when the scattering medium is homogeneous; it deteriorates as the medium becomes more heterogeneous, particularly with mixtures of soil, snow, and vegetation. There are genuine physical effects as well as more mundane spurious instrumental issues that degrade the fit but do not necessarily cause a bias in snow-depth estimates. These include secondary reflections, interferometric power effects, direct power effects, and instrument-related issues.

FIGURE 5. Examples of observations: (a) good fit; (b) presence of secondary reflections; (c) vanishing interference fringes; (d) atypical interference fringes.

Secondary reflections originate from disjoint surface regions. Interference fringes become convoluted with multiple superimposed beats (see Figure 5(b)). As long as there is a unique dominating reflection, the inversion will have no difficulty fitting it, as the extra reflections will remain approximately zero-mean.

Random deviations of the actual surface with respect to its undulated approximation, called roughness or residual surface height, will affect the interferometric power. SNR measurements will exhibit a diminishing number of significant interference fringes, compared to the measurement noise level (see Figure 5(c)). This facilitates the model fit but the reflector height parameter may become ill-determined: its estimates will be more uncertain. Changes in snow density also affect the fringe amplitude.

Snow precipitation attenuates the satellite-to-ground radio link, which affects SNR measurements through the direct power term. First, this shifts the SNR measurements up or down (in decibels); second, it tilts the trend tSNR as attenuation is elevation-angle dependent; third, fringes in dSNR will change in amplitude because of the decrease in the coherent component of the direct power.

Partial obstructions can affect either or both direct and interferometric powers. In this case, SNR measurements, albeit corrupted, are still recorded. This situation is in contrast to complete blockages as caused by topography. The deposition of snow and the formation of a winter rime on the antenna are a particularly insidious type of obstruction, as their presence in the near-field of the antenna element can easily distort the gain pattern in a significant manner. In the far-field, trees are another important nuisance, so much so that their absence is held as a strong requirement for the proper functioning of multipath reflectometry.

Satellite-specific direct power offsets and also long-term power drifts are to be expected as spacecraft age and modernized designs are launched. In addition, noise power depends on the state of conservation of receiver cables and on their physical temperature. Less subtle incidents are sudden ~3-dB SNR steps, hypothesized to originate in the receiver switching between the L2C data and pilot subcodes, CM and CL.

Quality Control. Anomalous conditions may result in measurement spikes, jumps, and short-lived rapidly-varying fluctuations. For snow-depth-sensing purposes, it is necessary and sufficient to either neutralize such measurement outliers through a statistically robust fit or detect unreliable fits and discard the problematic ones that could not otherwise be salvaged.

The key to quality control (QC) is in grouping results into statistically homogeneous units, having measurements collected under comparable conditions. In our case, azimuth-clustered tracks are the natural starting unit. Secondarily, we must account for genuine temporal variations in the tendency of results, from beginning to peak to the end of the snow season. The detection of anomalous results further requires an estimate of the statistical dispersion to be expected. Considering that the sample is contaminated with outliers, robust estimators (running median instead of the running mean, and median absolute deviation over the standard deviation) are called for, if the first- and second-order statistical moments are to be representative. Given estimates of the non-stationary tendency and dispersion, a tolerance interval can then be constructed such that it bounds, say, a 99% proportion of the valid results with 95% confidence level. We also desire QC to be judicious, or else too many valid estimates will be lost. Notice that in the present intra-cluster QC, we compare an individual estimate to the expected performance of the track cluster to which it belongs; later, we complement QC with an inter-cluster comparison of each cluster’s own expected performance.

Based on our practical experience, no single statistic detects all the outliers. We use four particular statistics that we have found to be useful: 1) degrees of freedom, essentially the number of observations per track (modulo a constant number of parameters); 2) using the scaled root-mean-square error (RMSE) to test for goodness-of-fit, that is, how well measurements can be explained adjusting the unknown values for the parameters postulated in the model; 3) reflector height uncertainty; and 4) peak elevation angle, which behaves much like a random variable, as it is determined by a multitude of factors.

Combinations. We combine multiple clusters to average out random noise. Noise mitigation aims at not only coping with measurement errors but also compensating for model deficiencies, to the extent that they are not in common across different clusters. Before we combine different clusters, we have to address their long-term differences. The initial situation is that snow surface heights will be greater downhill and smaller uphill; we take this into account on a cluster-by-cluster basis by subtracting ground heights from their respective snow surface heights, resulting in snow thickness values, which is a completely physically unambiguous quantity. Snow thickness is more comparable than snow heights across varying-azimuth track clusters. Yet snow tends to fill in ground depressions, so thickness exhibits variability caused by the underlying ground surface, even when the overlying snow surface is relatively uniform. Further cluster homogeneity can be achieved by accounting for the temporally permanent though spatially non-uniform component of snow thickness.

The averaging of snow depths collected for different track clusters employs the inversion uncertainties to obtain a preliminary running weighted median, calculated for, say, daily postings, with overlapping windows or not. The preliminary post-fit residuals then go through their own averaging, necessarily employing a wider averaging window (say, monthly), which produces scaling factors for the original uncertainties. The running weighted median is then repeated, producing final averages. The variance factors reflect the fact that some clusters are better than others.

Thus, the final GPS estimates of snow depth follow from an averaging of all available tracks, whose individual snow depth values were previously estimated independently. A new average is produced twice daily utilizing the surrounding 1–2 days of data (depending on the data density), that is, 12-hour posting spacing and 24-hour moving window width. The averaging interval must be an integer number of days, so as to minimize the possibility of snow-depth artifacts caused by variations in the observation geometry, which repeats daily.

Site-Specific Results

We explored GPS-MR snow-depth retrieval at three stations over a long period (up to three years). Throughout, we assessed the performance of the GPS estimates against independent nearly co-located in situ measurements. We also compared the GPS estimates to the nearest SNOTEL station. SNOTEL (from snowpack telemetry) is an automated system for collecting snowpack and related data in the western U.S. operated by the U.S. Department of Agriculture. Although not co-located with GPS, SNOTEL data are important because they provide accurate information on the timing of snowfall events.

The three sites we used were 1) a site in the T.W. Daniel Experimental Forest within the Wasatch Cache National Forest in the Bear River Range of northeastern Utah, with an elevation of 2,600 meters; 2) one of the stations of the EarthScope Plate Boundary Observatory, a grassland site located near Island Park, Idaho; and 3) an alpine site in the Niwot Ridge Long-term Ecological Research Site near Boulder, Colorado. While we have fully documented the results from each site, due to space limitations we will only discuss the results from the forested site (known as RN86) in this article. This is a more challenging site than the other two, due to the presence of nearby trees. Furthermore, it was subject to denser in situ sampling of 20-150 measurements spatially replicated around the GPS antenna, and repeated approximately every other week for about one year.

We show results for the 2012 water-year, the period starting October 1 through September 30 of the following year. Where GPS site RN86 was installed, topographical slopes range from 2.5° to 6.5° (at the 2-meter spatial scale), with average of ~5° within a 50-meter radius around the GPS antenna. RN86 was specifically built to study the impact of trees on GPS snow depth retrievals (see FIGURE 6). Ground crews manually collected in situ measurements around the GPS antenna approximately every other week starting in November 2011. Measurements were made every 1–2 meters from the antenna up to a distance of 25-30 meters. In the second half of the year, the sampling protocol was changed to azimuths of 0° (N), 45° (NE), 135° (SE), 180° (S), 225° (SW), and 315° (NW). With these data it is possible to obtain in situ average estimates, with their own uncertainties (based on the number of measurements), which allows a more meaningful comparison.

FIGURE 6. Aerial view of the forested site (RN86) around the GPS antenna (marked with a circle).

There is reduced visibility at the current site, compared to other sites. Track clusters are concentrated due south, with only two clusters located within ±90° of north. Therefore, the GPS average snow depth is not necessarily representative of the azimuthally symmetric component of the snow depth. In the presence of an azimuthal asymmetry in the snow distribution around the antenna, the GPS average would be expected to be biased towards the environmental conditions prevalent in the southern quadrant. To rule out the possibility of an azimuthal artifact in the comparisons, we have utilized only the in situ data collected along the SE/S/SW quadrant.

The comparison shows generally excellent agreement between GPS and in situ data (see FIGURE 7). The first four and the last one in situ data points were collected with coarser spacing and/or smaller azimuthal coverage, which may be partially responsible for different performance in the first and second halves of the snow season. The correlation between GPS and in situ snow depth at RN86 amounts to 0.990, indicating a very strong linear relationship. Carrying out a regression between in situ and GPS values, the RMS of snow-depth residuals improves from 9.6 to 3.4 centimeters. The regression intercept and slope (with corresponding 95% uncertainties) amount to 15.4 ± 9.11 centimeters and 0.858 ± 0.09 meters per meter, respectively. According to these statistics, the null hypotheses of zero intercept and unity slope are rejected at the 95% confidence level. This implies that at this location GPS snow-depth estimates exhibit both additive and multiplicative biases. The latter is proportional to snow depth itself, meaning that, compared to an ideal one-to-one relationship, GPS is found to under-estimate in situ snow depth at this site by 14 ± 9%, although the uncertainty is somewhat large.

FIGURE 7. Snow-depth measurement at the forested site (RN86) for the water-year 2012

The SNOTEL sensors are exceptionally close to the GPS antenna at this site, about 350 meters horizontally distant with negligible vertical separation. Yet the former is located within trees, while the latter is located at the periphery of the forest and senses the reflections scattered from an open field. Therefore, only the timing of snowfall events agrees well, not the amount of snow. Although forest density is generally negatively correlated with snow depth, exceptions are not uncommon, especially in localized clearings exposed to intense solar radiation, where shading of the snow by the trees reduces ablation.

Conclusions

In this article, we have discussed a physically based forward model and a statistical inverse model for estimating snow depth based on GPS multipath observed in SNR measurements. We assessed model performance against independent in situ measurements and found they validated the GPS estimates to within the limitations of both GPS and in situ measurement errors after the characterization of systematic errors. The assessment yielded a correlation of 0.98 and an RMS error of 6–8 centimeters for observed snow depths of up to 2.5 meters at three sites, with the GPS underestimating in situ snow depth by ~5–15%. This latter finding highlights the necessity to assess effects currently neglected or requiring more precise modeling.

Acknowledgments

The research reported in this article was supported by grants from the U.S. National Science Foundation, NASA, and the University of Colorado. Nievinski has been supported by a Capes/Fulbright Graduate Student Fellowship and a NASA Earth System Science Research Fellowship. The article is based, in part, on two papers published in the IEEE Transactions on Geoscience and Remote Sensing: “Inverse Modeling of GPS Multipath for Snow Depth Estimation – Part I: Formulation and Simulations” and “Inverse Modeling of GPS Multipath for Snow Depth Estimation – Part II: Application and Validation.”

Manufacturers

For the forested site (RN86), a Trimble NetR9 receiver was used with a Trimble TRM57971.00 (Zephyr Geodetic II) antenna with no external radome.

FELIPE G. NIEVINSKI is a faculty member at the Federal University of Santa Catarina, Florianópolis, Brazil. He has also been a post-doctoral researcher at São Paulo State University, Presidente Prudente, Brazil. He earned a B.E. in geomatics from the Federal University of Rio Grande do Sul, Porto Alegre, Brazil, in 2005; an M.Sc.E. in geodesy from the University of New Brunswick, Fredericton, Canada, in 2009; and a Ph.D. in aerospace engineering sciences from the University of Colorado, Boulder, in 2013. His Ph.D. dissertation was awarded The Institute of Navigation Bradford W. Parkinson Award in 2013.

KRISTINE M. LARSON received a B.A. degree in engineering sciences from Harvard University and a Ph.D. degree in geophysics from the Scripps Institution of Oceanography, University of California at San Diego. She was a member of the technical staff at the Jet Propulsion Lab from 1988 to 1990. Since 1990, she has been a professor in the Department of Aerospace Engineering Sciences, University of Colorado, Boulder.

FURTHER READING

• Authors’ Journal Papers

“Inverse Modeling of GPS Multipath for Snow Depth Estimation—Part I: Formulation and Simulations” by F.G. Nievinski and K.M. Larson in IEEE Transactions on Geoscience and Remote Sensing, Vol. 52, No. 10, 2014, pp. 6555–6563, doi: 10.1109/TGRS.2013.2297681.

“Inverse Modeling of GPS Multipath for Snow Depth Estimation—Part II: Application and Validation” by F.G. Nievinski and K.M. Larson in IEEE Transactions on Geoscience and Remote Sensing, Vol. 52, No. 10, 2014, pp. 6564–6573, doi: 10.1109/TGRS.2013.2297688.

• More on the Use of GPS for Snow Depth Assessment

“Snow Depth, Density, and SWE Estimates Derived from GPS Reflection Data: Validation in the Western U.S.” by J.L. McCreight, E.E. Small, and K.M. Larson in Water Resources Research, published first on line, August 25, 2014, doi: 10.1002/2014WR015561.

“Environmental Sensing: A Revolution in GNSS Applications” by K.M. Larson, E.E. Small, J.J. Braun, and V.U. Zavorotny in Inside GNSS, Vol. 9, No. 4, July/August 2014, pp. 36–46.

“Snow Depth Sensing Using the GPS L2C Signal with a Dipole Antenna” by Q. Chen, D. Won, and D.M. Akos in EURASIP Journal on Advances in Signal Processing, Special Issue on GNSS Remote Sensing, Vol. 2014, Article No. 106, 2014, doi: 10.1186/1687-6180-2014-106.

“GPS Snow Sensing: Results from the EarthScope Plate Boundary Observatory” by K.M. Larson and F.G. Nievinski in GPS Solutions, Vol. 17, No. 1, 2013, pp. 41–52, doi: 10.1007/s10291-012-0259-7.

• GPS Multipath Modeling and Simulation

“Forward Modeling of GPS Multipath for Near-Surface Reflectometry and Positioning Applications” by F.G. Nievinski and K.M. Larson in GPS Solutions, Vol. 18, No. 2, 2014, pp. 309–322, doi: 10.1007/s10291-013-0331-y.

“An Open Source GPS Multipath Simulator in Matlab/Octave” by F.G. Nievinski and K.M. Larson in GPS Solutions, Vol. 18, No. 3, 2014, pp. 473–481, doi: 10.1007/s10291-014-0370-z.

“Multipath Minimization Method: Mitigation Through Adaptive Filtering for Machine Automation Applications” by L. Serrano, D. Kim, and R.B. Langley in GPS World, Vol. 22, No. 7, July 2011, pp. 42–48.

“It’s Not All Bad: Understanding and Using GNSS Multipath” by A. Bilich and K.M. Larson in GPS World, Vol. 20, No. 10, October 2009, pp. 31–39.

“GPS Signal Multipath: A Software Simulator” by S.H. Byun, G.A. Hajj, and L.W. Young in GPS World, Vol. 13, No. 7, July 2002, pp. 40–49.

September 1, 2014
Innovation: Not Just a Fairy Tale

A Hansel and Gretel Approach to Cooperative Vehicle Positioning

By Scott Stephenson, Xiaolin Meng, Terry Moore, Anthony Baxendale, and Tim Edwards

MEET GEORGE JETSON.Those of us of a certain age will remember the animated TV sitcom The Jetsons, which featured George Jetson, “his boy Elroy, daughter Judy, and Jane, his wife.” It portrayed life in 2062, 100 years after the series debuted in 1962. George and his family used many futuristic gadgets including robot maids, talking alarm clocks, flat-screen TVs, and flying automated cars. Many of those devices are already available, well ahead of schedule. But flying cars are not quite with us yet. However, asphalt-hugging automated vehicles are already here, albeit still in limited numbers. Google created a buzz recently with tests of its self-driving car. Google’s cars were developed as an outcome of the Defense Advanced Research Projects Agency’s 2005 Grand Challenge in which teams created autonomous vehicles and raced them through a challenging road course.

Self-driving cars use a host of sensors to determine their position with respect to their surroundings and to navigate a chosen route legally and safely. Although wide-spread ownership of self-driving cars might still be a ways off, drivers of conventional vehicles will soon benefit from the research being conducted to provide them with positional awareness of other vehicles in their vicinity. This work may be characterized as part of the larger effort in developing intelligent transportation systems or ITS.

What is ITS? In the words of ITS Canada, it’s “the application of advanced and emerging technologies (computers, sensors, control, communications, and electronic devices) in transportation to save lives, time, money, energy and the environment.” This definition applies to all modes of transportation, including ground transportation such as private automobiles, commercial vehicles, and public transit, as well as rail, marine, and air modalities. The term ITS includes consideration not only of the vehicle, but also the infrastructure, and the driver or user, interacting together dynamically.

Just looking at ground transportation, there are many ITS developments underway, some of which are already implemented to some degree including systems for vehicle navigation, traffic-signal-control, automatic license-plate recognition, parking guidance, and road lighting to name but a few.

An important aspect of ITS is cooperative vehicle communication, which includes transmission of data vehicle–to–vehicle or vehicle–to–infrastructure (and vice versa — known by the abbreviation V2X. Data from vehicles can be acquired and transmitted to other vehicles or to a server for central fusion and processing. These data can include accurate real-time vehicle coordinates, which can be used to improve driver situational awareness and to monitor traffic flow for example. This use of V2X is known as cooperative vehicle positioning.

Several technologies are being developed for accurate cooperative vehicle positioning including lidar, radar, image-based cameras, ultra-wideband, and signals of opportunity. But GNSS also has a role to play. In this month’s column, team of British researchers turn to a children’s fairy tale for inspiration in their development of a cooperative vehicle positioning approach using carrier-phase observations — another innovative application of real-time kinematic or RTK GNSS technology.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas.

There is little doubt in the benefit gained from cooperative modes of road transport, as agents working together generally perform better. In simple terms, this is the holistic idea that the whole is greater than the sum of its parts, commonly known as synergy. On top of this clear advantage, the complex systems theory of emergence suggests that novel strategies will develop from the as-yet-undefined patterns and structures. It is clear, however, that to facilitate this development certain technological advances need to be achieved. In this case, individual road agents need to accurately identify their location, and communicate easily and safely with other agents. This is a shift away from protective and passive systems toward preventative and active transport safety.

Cooperative driving, or vehicle-to-vehicle or vehicle-to-infrastructure driving (V2X), is proposed as the next major safety breakthrough in road transport. An example of the concept is shown in FIGURE 1. It involves agents in the road transport environment communicating on local and national levels in real time, to maximize the efficiency of movement, dramatically reduce the number of accidents and fatalities, and make transportation more environmentally friendly.

Figure 1. Vehicle-to-vehicle communications as envisioned by the United States Department of Transportation.

In the U.S., the National Highway Traffic and Safety Administration has commented that connected vehicle technology “can transform the nation’s surface transportation safety, mobility and environmental performance,” with industry experts predicting the widespread uptake of the technology within five to six years. This provides an opportunity for road vehicles to share GNSS information.

To an extent, this is possible with current technology. Communication is fairly pervasive and pretty robust, with the explosion in personal handheld mobile devices, using the GSM/GPRS, 3G, and 4G cellular communications networks. Positioning systems exist now that will provide a reasonably accurate and reliable location most of the time. However, the type of applications included in cooperative driving demand much higher performance from these positioning systems. For instance, as shown in the example in FIGURE 2, two vehicles approaching an intersection at relatively high speeds require accurate and reliable high output position information, and an ability to communicate with one another, in order to assess the likelihood of collision.

Figure 2. Vehicles approaching a road intersection would benefit from V2X communication.

These requirements are partly inter-linked, and can be mutually beneficial. For instance, communications methods can be used to share information to aid positioning, and some existing positioning systems can also be utilized to share information.

Many recent solutions in vehicle tracking research have shifted the GNSS receiver to a supplemental role in the positioning system, favoring an inertial device as the core of the integrated solution. The clear advantage is that an inertial device operates continuously, although other sensors are required to achieve the required navigation performance. The GNSS receiver is demoted because of its inherent limitations, namely the requirement of a clear view of the satellites and the availability of correctional information.

Most vehicle positioning research over the past two decades has focused attention on GNSS-centered systems, as evidenced by the abundant use of satnav devices used to assist in-car navigation. Despite its apparent monopoly over vehicle positioning in the commercial sector, the most
successful systems developed to guide autonomous vehicles either relegate GNSS to one of a suite of sensors, or almost disregard it altogether. This is often due to its apparent lack of positioning accuracy or availability. Popular terrestrial positioning sensors include lidar, radar, image-based cameras, ultra-wideband (UWB), and signals of opportunity. Clearly, the combination of different complementary sensors is important, but it would be a mistake to discount the more advanced GNSS positioning techniques that are available, especially with the expansion of the four global GNSS services.

Cooperative Positioning

The positioning of GNSS receivers relative to one another is a common application in transportation, such as during the aerial refueling of an airborne fighter jet by a tanker. In this case, it is important to know accurately the relative position of the two airplanes, but not necessarily their absolute position.

Relative positioning of road vehicles is more complex. By their nature, road vehicles are almost always close to other vehicles or road infrastructure, and there are many separate agents in each scenario. Vehicles can also travel large distances, and in terms of GNSS positioning, this may mean vastly different atmospheric conditions. Hence, relative positioning in road transport is useful if all GNSS receivers relate to the same datum, which in most cases is effectively absolute positioning.

Some previous work carried out by others concentrated on using GNSS code (pseudorange) and Doppler measurements for the relative positioning of vehicles, because it offers a simpler implementation method and is not susceptible to the cycle slips attributed to carrier-phase measurements. However, this means sacrificing the higher accuracy solution available from carrier-phase measurements. A major obstacle to GNSS positioning for V2X applications is the likely scenario of mixed receiver and antenna technology between vehicles. This has a major influence on the performance of relative positioning. By comparing various V2X relative positioning solutions, researchers found that an increase in positioning accuracy was typically accompanied by a decrease in availability and an increased demand for transmission bandwidth between the vehicles.

RTK GNSS Positioning. Real-time kinematic (RTK) GNSS positioning can be used to provide a solution at an accuracy of better than 5 centimeters (horizontal). This relies on the static reference receiver being located within 20 kilometers of the roving receiver, observing a good selection of common satellites with dual-frequency receivers.

When RTK positioning is used, the distance to the reference station has a bearing on the successfulness of the integer ambiguity resolution. A short baseline will benefit from a closer correlation of errors, due to the GNSS signals traveling through very similar parts of the atmosphere. Assuming each receiver is observing common satellites, this similarity will typically result in a higher success rate in the ratio test using the common Least Squares Ambiguity Decorrelation Adjustment, or LAMBDA, technique. This is particularly important following a GNSS outage.

GNSS positioning of road vehicles using RTK or network RTK (where a network of reference stations replaces a single RTK reference station) can provide highly accurate (< 5 centimeters), high integrity, real-time tracking information with little delay and at a high output rate. The proliferation of network RTK GNSS positioning systems has increased dramatically over the last decade. Network RTK GNSS positioning can minimize the spatial decorrelation of errors that is a characteristic of single-reference RTK positioning as distance increases between reference and rover receivers. This allows the wide mobility range demanded from automotive applications.

The transmission protocol of network RTK corrections is typically RTCM v3.0 or higher, and the composition of the correction information varies depending on the commercial service provider. The most common type of correction message format is that for a virtual reference station (VRS), although the most comprehensive and versatile method is the master-auxiliary concept (MAC). See references in Further Reading for details.

In V2X and other intelligent transportation systems (ITS) applications, the position must be accurate, reliable, available, and continuous. Previous research has shown that network RTK GNSS positioning can deliver a highly accurate and precise solution in an ideal observation environment. In one test, more than 99 percent of the observations lay within 2 centimeters of the truth solution, with a very small number of anomalous results of up to 20 centimeters.

The availability of a network RTK solution is determined by the availability of GNSS signals and the network RTK corrections. As network RTK positioning uses carrier-phase observations, GNSS outages and cycle slips significantly affect the performance of a receiver. However, the re-initialization of the fixed integer ambiguity resolution following a GNSS outage (such as caused by an overhead bridge) can be relatively fast. But from a cold start, the ambiguity resolution can take up to two minutes. This limits the widespread adoption of the technology for vehicle positioning.

NGI Road Vehicle and Electric Locomotive Testbeds. We have carried out research at the Nottingham Geospatial Institute (NGI) using state-of-the-art testing facilities. These bespoke in-house facilities allow repeated controlled experiments, and are a useful tool in the development of ITS and V2X technology.

To test the positioning performance thoroughly and under real-world conditions, we carried out experiments using the NGI’s road vehicle, which is equipped with a collection of on-board ground-truth systems.

Also, the roof of the Nottingham Geospatial Building (home of NGI) is the location of a remotely operated electric locomotive running on a 200-millimeter-gauge railway track. A photograph of the locomotive and plan of the track are shown in FIGURE 3. The locomotive can carry a selection of various positioning instruments, such as GNSS receivers, inertial navigation system (INS) devices, and tracking prisms, and can travel at a speed of over three meters per second. The position of the track is accurately known, and has previously been scanned at a resolution of 2 millimeters.

Figure 3. The NGB2 reference base station and electric locomotive track on the roof of the Nottingham Geospatial Building.

Three control solutions are used to assess the performance of the cooperative positioning techniques in real-world tests: An RTK GNSS control solution provided by a local static continuously operating reference station (CORS); a network RTK GNSS solution based on the MAC standard; and a
dual-frequency GPS/INS system. Each vehicle also can be independently tracked using survey-grade total stations or a proprietary UWB positioning system.

Sharing Network RTK Corrections

If vehicles could communicate with one another on the road, this would help overcome the communication system limitation in network RTK positioning of road vehicles. For instance, if vehicle A has an external connection to a network RTK service provider (such as a mobile Internet connection) and a local connection to a second vehicle (B), then it could share its network RTK correction messages directly. Effectively, vehicle A would re-broadcast the correction information it has received from the corrections provider to the receiver on vehicle B. However, this would rely on the functional capability of the receiver of vehicle B, as network RTK real-time processing can be computationally intensive.

Not all network RTK correction messages can be shared in this way, and the range over which the correction messages are still valid needs to be determined. As vehicles communicating with V2X devices are likely to be relatively close (a few hundred meters at most), the feasibility of sharing network RTK information is good.

However, the network RTK VRS technique may offer more advantages. It is the most common form of network RTK used around the world, and requires significantly less bandwidth (approximately 10 kilobits per second at 10 Hz). The rover receiver is also less burdened by processing requirements. A VRS system operating on buses in Minnesota restricts the baseline to 2 miles, by updating the VRS location every 2 minutes.

Correction messages typically have a lifespan of 10 seconds. After this time, the receiver determines the messages to be too old and does not compute a fixed-integer position. It can, however, use the information to calculate a differential GNSS (DGNSS) position. Therefore, the relayed message must arrive at the receiver on vehicle B well within 10 seconds. Previous trials at NGI found that the typical message latency of the original correction message reaching vehicle A via a GSM/GPRS connection is 0.85 seconds. The additional V2X communication to transfer the message to vehicle B should not add a significant delay.

Capturing Network RTK Messages. To demonstrate the potential benefit of sharing network RTK messages between vehicles, network RTK messages were captured on board a vehicle and shared with a second vehicle. Vehicle A is the NGI van, and vehicle B is the NGI electric train. Most off-the-shelf network-RTK-enabled GNSS receivers are designed to communicate directly with the network RTK server using a connected communication device (GSM modem, UHF/VHF radio, cell phone, and so on), which typically provides a stable connection to minimize data loss.

To intercept the network RTK correction message, the GNSS receiver was set up to simply accept the correction message from a smartphone via Bluetooth. In this case, the connection to the network RTK service provider is established between the smartphone and the network RTK server. An application running on the smartphone (as shown in FIGURE 4) requests information from the network RTK server, logs the data, and passes the message directly to the Bluetooth-connected GNSS receiver on vehicle A. By intercepting the correction message, it can also be forwarded on to a second receiver, in this case on vehicle B.

Figure 4. Flowchart showing the capturing and sharing of network RTK correction messages (left), and the NTRIP client program running on an Android smartphone (right).

Sharing Messages with Second Receiver. FIGURE 5 shows the positioning solutions generated by a shared-network-RTK correction message. The original message was captured by the smartphone application operating on board vehicle A (the NGI van), and applied to GNSS observations made by a receiver on vehicle B (the NGI train). The baseline between the two vehicles was less than 100 meters, and the location of the VRS requested from the network RTK server was the NGI building (in geodetic coordinates to three decimal places). As Figure 5 clearly shows, the shared VRS corrections are equally valid for any receiver operating in the vicinity of the VRS. The thick red line is the fixed position of the train track, and the thin blue line represents the positions generated by the GNSS receiver using the shared network RTK corrections.

Figure 5. Sharing the network RTK message from vehicle A to vehicle B.

The VRS message type was chosen because it requires much less bandwidth, takes less processing capacity, and is prevalent among legacy receivers. Network RTK users typically require download speeds of 1.8 kilobits per second (VRS) and 5.6 kilobits per second (MAC). This is well within the typical speeds available from cellular wireless communications, which offer 80 kilobits per second downlink speeds from 2.5G systems to beyond 40 megabits per second for recent 4G systems.

The GNSS receiver on vehicle B is operating in an ideal location, with a clear view of the sky and a high number of visible satellites, which improves the probability of successful RTK ambiguity resolution.

Generating Pseudo-VRS Corrections

The potential benefit to GNSS positioning of using V2X communication between various road vehicles and infrastructure can be expanded by the implementation of pseudo-VRS positioning. This system resembles the children’s fairy tale Hansel and Gretel, where in order to help remember the route through a forest that guides them back to their home, Hansel drops markers along the path (in separate cases small white pebbles, and then breadcrumbs). By using the markers, the children can navigate their way through the forest, but without them they are left lost and disoriented.

The pseudo-VRS system uses a similar principle, where vehicle A marks its path by leaving behind small packets of information that can be used by other nearby vehicles. The small packets of information are VRS-like, and are broadcast using V2X communication devices and technology. Like the breadcrumbs in the fairy tale that are eaten by birds shortly after being dropped by Hansel, these VRS-like packets of information have a short lifespan.

VRS Requirements. It has been long established that a short baseline between reference and rover receivers leads to more accurate and successful relative GNSS positioning. A short baseline can effectively deal with satellite orbit and atmospheric errors, which become difficult to deal with as the baseline length grows, and is the reason why RTK GNSS positioning is typically limited to baselines shorter than 20 kilometers. A typical RTK baseline may be between 1 and 10 kilometers long, but it is still beneficial to reduce the baseline further, particularly if there is a large difference in elevation. This is enabled by the VRS network RTK technique. By using the observation data from several permanent reference stations that surround the rover location, a virtual reference station is created close to the location of the rover, including virtual observation measurements and position. This VRS information is transmitted to the rover, and the rover receiver treats the information like that of a real reference station. This technique can deliver better than 5-centimeter accuracy up to 35 kilometers.

The principle builds on the transfer of measurements made at the real reference stations to the VRS. The carrier-phase measurement at the real reference station (  ), shown in Equation 1, is made up of the geometric distance between the receiver and satellite (   ), the integer ambiguity (   ), and the receiver and satellite clock bias ( ). The key to the VRS technique is that the integer ambiguity and the receiver and satellite clock bias are not location dependent, so they can be transferred directly to the virtual reference station from the real reference station.

   (1)

By differencing the carrier-phase equation of the real and virtual reference stations ( and , respectively), the ambiguity and clock errors are canceled. The result is shown in Equation 2.

(2)

By combining the carrier-phase measurement equations at the real and virtual reference stations, only two unknown terms remain. The first includes the position of the VRS (   ), which is, in principle, arbitrary and is typically the approximate location of the rover receiver. The second is the observable of the VRS ( ), which can now be obtained without actually measuring it. (In practice, the technique is a little more complex, as satellite orbit and atmospheric errors and biases need to be modeled for the VRS position). The VRS information can then be packaged using the RTCM standards and delivered to the rover receiver to enable network RTK VRS positioning.

Pseudo-VRS. Using the established VRS techniques and standards described above, we propose to use the GNSS observations and subsequent position information to simulate the existence of a VRS (see FIGURE 6). Imagine vehicle A carries a GNSS receiver together with the means to calculate its position accurately (for instance, it is also receiving differential corrections or has other positioning devices on board). So long as the receiver can successfully resolve the integer ambiguity, it can also produce each component required to describe a VRS. Clearly in this case, the receiver on vehicle A is a “real” reference station, but the existing VRS standards can be exploited to transfer this information to other local GNSS receivers. For instance, a receiver operating on vehicle B can use the information as a local real-time differential correction service.

Figure 6. The flow of data during the generation and sharing of pseudo-VRS data.

Because the VRS technique is well established (the most popular form of network RTK positioning), legacy receivers are able to take advantage of this pseudo-VRS information. RTCM standards are also well defined for the transfer of GNSS information in this form.

The pseudo-VRS information is valid for several seconds, so the delays introduced in transferring the information from one vehicle to a second can easily be accommodated. Like any communication device based on radio waves, V2X communication devices are likely to be subject to a level of delay and message loss that requires redundancy in the system. It is important that during one epoch the whole pseudo-VRS message is delivered, as there is little similarity between one epoch and the next. The original reference receiver is likely to be on a moving vehicle.

Effectively, the pseudo-VRS imitates the VRS in Equation 2 by providing the virtual reference station coordinates and carrier-phase observable. The information is also delivered to the second receiver in the same format RTCM message. A slight difference here is that only one-way communication is needed — the original coordinates of the VRS do not need to be supplied by the second receiver.

The pseudo-VRS processing is carried out using the RTKLIB open source software. RTKLIB has limited options to vary the position of the base station during RTK positioning, so the program is seeded with customized configuration files and run independently for each epoch. This creates an additional feature: The processing of each epoch has no effect on any other.

Vehicle-to-Vehicle Communication. As we just consider the exploitation of V2X devices in this article, the nature of the communication medium is not under test. For this reason, off-the-shelf wireless routers (2.4 GHz) were used to communicate between vehicles, using fixed local IP addresses. However, the performance of the routers under cooperative driving tests is limited by range, multipath, and signal obstruction.

Real-World Tests

To generate significant test results, some of the following tests use recorded and replayed data.

Test Setup. To test the performance of a pseudo-VRS positioning system, and the success of different configurations, real-world tests were carried out at the Nottingham Geospatial Institute. Two vehicles were used. Vehicle A was the NGI’s road vehicle, and vehicle B was the NGI’s electric locomotive. As the position of the locomotive test track is very accurately known, this can be used to measure the performance of the pseudo-VRS system.

Vehicle A was equipped with six GNSS receivers, a tactical-grade INS system, and a wheel odometer, and tracked using a total station and 360º prism. This provided multiple position solutions to ensure significant results.

Vehicle B was equipped with a GNSS receiver, and tracked using a proprietary UWB system for related V2X tests.

Also, on the roof of the NGB, and lying inside the track perimeter, is the NGB continuously operating reference
station. This hyper-local reference station allows local RTK solutions, and acts as a barometer of GNSS activity when tests are episodically carried out.

FIGURE 7 shows an aerial image of the test scenario. The Google background shows the NGB to the west, and surrounding roads to the south and west (still under construction during the image acquisition). The thin yellow line is a ground distance of 100 meters. The red dots signify the position of vehicle A (in the east), and the purple dots show the position of vehicle B (on the roof of the NGB building). The accuracy of the Google image is unknown, and is used here purely for illustrative purposes.

Figure 7. Aerial image of the test.

Test Results. These tests are designed to show the performance of a pseudo-VRS system using a V2X communication system. However, the results shown here were created using recorded raw data. The test results will help to design the correct RTCM message to share between vehicles in future tests.

To simulate the operation of a pseudo-VRS system, vehicle A must share its known absolute position and some raw RINEX information for each epoch with vehicle B. Vehicle B can then use this information, together with its own observed RINEX data, for the same epoch to calculate its known absolute position. In practice, there will be a slight delay in the delivery of the information from vehicle A (much like in a traditional RTK system), so that information from concurrent epochs are unlikely to be used.

The RTKLIB software cannot directly handle the variation of a base station’s coordinates (and output an absolute solution), so a small separate script was designed to utilize the processing capability of the software in a pseudo-VRS system.

FIGURE 8 shows the results of pseudo-VRS positioning. During dual-frequency tests, 99.67 percent of observations achieved fixed ambiguity (1197/1201). During single-frequency (broadcast ionosphere) RTK, 61.45 percent (738/1201) observations achieved fixed ambiguity. The ratio test threshold was 2.0. Around the area of 454930E 339708N, the number of common visible satellites dropped from eight to seven, and then again from seven to six three seconds later. This caused each of the three solutions to degrade slightly. The dual-frequency RTK solution briefly lost its fixed ambiguity solution (for two epochs, or 0.1 seconds), before regaining the fixed solution. The single-frequency RTK solution could not achieve a fixed ambiguity solution again until the number of common visible satellites returned to seven (five seconds after the initial satellite was lost). The DGNSS solution saw a similar degradation in its solution during this period.

Figure 8. Results from pseudo-VRS positioning.

The mean coordinate errors for the three solutions are 0.054, 0.707, and 0.323 meters (1 standard deviation, 3D), as shown in Table 1. This is compared to a solution calculated using the local CORS base station. The error in horizontal and vertical follows the typical ratio of 1:2. Test results were also completed using a lower pseudo-VRS update rate. At 1 Hz, the results prove even better. Although the latency of the correction is up to 1 second (positioning is calculated epoch by epoch), the results were better than updates at 20 Hz. The dual-frequency RTK solution achieved a fixed ambiguity at every epoch (100 percent), and when compared to the known track position appeared correctly fixed. The single-frequency RTK solution achieved a fixed ambiguity for 70.02 percent (897/1201) of the observations; a slight improvement over the 20-Hz results.

Table 1. Results from pseudo-VRS positioning.

Table 2 shows the performance of the pseudo-VRS system under different latency scenarios. This is important because a message transmitted by vehicle A may be delayed or newer messages may be disrupted. Once the latency of the correction message reaches 8 seconds, the performance of the positioning solution begins to drop. The number of fixed ambiguity solutions falls, and the resulting positioning accuracy also decreases. However, the solution can still deliver 20- to 30-centimeter accuracy with a message latency of up to 30 seconds.

Table 2. Effect of message latency on positioning quality.

Conclusions

This article has outlined the potential benefit of V2X technology to cooperative vehicle positioning. A vehicle that knows its absolute position accurately can assist a second vehicle to position itself using established GNSS techniques.

The pseudo-VRS base-station location must have reasonably accurate coordinates. Without this, the correct integer ambiguity cannot be resolved, and there is the risk of an incorrect resolution giving false success. This requires good reliability and integrity of the position of vehicle A, a characteristic that can be provided by network RTK positioning but likely needs further support from alternative positioning solutions.

Acknowledgments

The authors acknowledge Leica Geosystems for the provision of an academic license for the SmartNet network RTK service. We thank Yang Gao and Qiuzhao Zhang of the University of Nottingham for their assistance and detailed discussion during the experimental tests. The work was supported by the U.K.’s Engineering and Physical Sciences Research Council. This article is based on the paper “A Fairy Tale Approach to Cooperative Vehicle Positioning” presented at the 2014 International Technical Meeting of The Institute of Navigation held in San Diego, California, January 27–29, 2014.

Manufacturers

For our tests, vehicle A (NGI’s road vehicle) was equipped with six Leica Geosystems AG GS10 GNSS receivers with individual AS10 antennas, an Applanix Corp. POS RS with Honeywell International Inc. CIMU tactical grade INS system, and was tracked using a Leica Nova TS50 total station. Vehicle B (NGI’s electric locomotive) was equipped with a Leica GS10 GNSS receiver and AS10 antenna.

SCOTT STEPHENSON is a postgraduate student at the Nottingham Geospatial Institute (NGI) within the University of Nottingham, Nottingham, U.K.

XIAOLIN MENG is an associate professor, theme leader for positioning and navigation technologies, and an M.Sc. course director at NGI.

TERRY MOORE is the director of NGI at UoN, where he is the professor of satellite navigation and an associate dean within the Faculty of Engineering.

ANTHONY BAXENDALE is head of Advanced Technologies & Research at MIRA Ltd. (formerly the Motor Industry Research Association), an automotive consultancy company headquartered near Nuneaton in Warwickshire, U.K.

TIM EDWARDS is a principal engineer responsible for intelligent mobility research activities within the Future Transport Technologies Group at MIRA Ltd.

FURTHER READING

• Authors’ Conference Paper

“A Fairy Tale Approach to Cooperative Vehicle Positioning” by S. Stephenson, X. Meng, T. Moore, A. Baxendale, and T. Edwards in Proceedings of ION ITM 2014, the 2014 International Technical Meeting of The Institute of Navigation, San Diego, California, January 27–29, 2014, pp. 431–440.

• Intelligent Transportation Systems

Proceedings of IEEE ITSC 2013, the 16th International IEEE Conference on Intelligent Transportation Systems, “Intelligent Transportation Systems for All Modes,” The Hague, The Netherlands, October 6–9, 2013.

Overview of Intelligent Transport Systems (ITS) Developments in and Across Transport Modes by G.A. Giannopoulos, E. Mitsakis, and J.M. Salanoca, Joint Research Centre Scientific and Policy Report EUR 25223 EN, Institute for Energy and Transport, Joint Research Centre, European Commission, 2012, doi: 10.2788/12881.

“How Google’s Self-Driving Car Works” by E. Guizzo in IEEE Spectrum Blog, October 18, 2011.

“Elbow Room on the Shoulder: DGPS-Based Lane-Keeping Enlists Laser Scanners for Safety and Efficiency” by C. Shankwitz in GPS World, Vol. 21, No. 7, July 2010, pp. 30–37.

“Driverless Cars” by R. Murray in Computing and Control Engineering, Vol. 18, No. 3, June-July 2007, pp. 14–17.

• GNSS and Inertial Navigation Systems

“GPS and Inertial Systems for High Precision Positioning on Motorways” by J.E. Naranjo, F. Jiménez, F. Aparicio, and J. Zato in Journal of Navigation, Vol. 62, No. 2, April 2009, pp. 351–363, doi: 10.1017/S0373463308005249.

• Vehicle-to-Vehicle and Vehicle-to-Infrastructure Technologies

“Implementation of V2X with the Integration of Network RTK: Challenges and Solutions” inProceedings of ION GNSS 2012, the 25th International Technical Meeting of The Satellite Division of the Institute of Navigation, Nashville, Tennessee, September 17–21, 2012, pp. 1556–1567.

DOT Launches Largest-Ever Road Test of Connected Vehicle Crash Avoidance Technology, National Highway Traffic Safety Administration press release, August 21, 2012.

“Relative Positioning for Vehicle-to-Vehicle Communication-enabled Vehicle Safety Applications” by C. Basnayake, G. Lachapelle, and J. Bancroft in Proceedings of the 18th ITS World Congress, Orlando, October 16–20, 2011.

“Can GNSS Drive V2X” by P. Alves, T. Williams, C. Basnayake, and G. Lachapelle in GPS World, Vol. 21, No. 10, October 2010, pp. 35–43.

• Network RTK

“Network RTK for Intelligent Vehicles” by S. Stephenson, X. Meng, T. Moore, A. Baxendale, and T. Edwards in GPS World, Vol. 24, No. 2, February 2013, pp. 61–67.

“A Comparison of the VRS and MAC Principles for Network RTK” by V. Janssen in Proceedings of IGNSS2009, the 2009 Symposium of the International Global Navigation Satellite Systems Society, Gold Coast, Queensland, Australia, December 1–3, 2009.

“Introduction to Network RTK” by L. Wanninger, IAG Working Group 4.1: Network RTK (2003–2007). Online article. Last modified June 16, 2008.

RTCM Standard 10403.1 for Differential GNSS (Global Navigation Satellite Systems) Services – Version 3, developed by RTCM Special Committee No. 104, Radio Technical Commission for Maritime Services, Arlington, Virginia, October 27, 2006.

“Accuracy Performance of Virtual Reference Station (VRS) Networks” by G. Retscher in Journal of Global Positioning Systems, Vol. 1, No. 1, 2002, pp. 40–47.

“An Overview of Multi-Reference Station Methods for cm-Level Positioning” by G. Fotopoulos and M.E. Cannon in GPS Solutions, Vol. 4, No. 3, January 2001, pp. 1–10, doi: 10.1007/PL00012849.

July 1, 2014
Innovation: The European Way
Performance of the Galileo Single-Frequency Ionospheric Correction During In-Orbit Validation

By Roberto Prieto-Cerdeira, Raül Orús-Pérez, Edward Breeuwer, Rafael Lucas-Rodriguez, and Marco Falcone

OFF TO A GOOD START. That’s how we might characterize the European Galileo satellite navigation system. The official beginning of the Galileo program occurred on May 26, 2003, when the European Union and the European Space Agency officially agreed on the first stage of the program (although some work on system concepts took place earlier). The first two prototype or development satellites, Galileo In-Orbit Validation Element-A (GIOVE-A) and GIOVE-B, were launched on December 28, 2005, and April 26, 2008, respectively. The satellites successfully validated key technologies for the full Galileo constellation and secured the system’s frequency allocations.

The first two In-Orbit Validation (IOV) satellites were launched by a single rocket on October 21, 2011, and the third and fourth IOV satellites were similarly launched on October 12, 2012. The two GIOVE satellites and first two IOV satellites provided an opportunity to use Galileo-only receiver measurements and after-the-fact precise satellite orbit and clock data to compute the position of a receiver’s antenna. Joined by two colleagues, I was pleased to report our successful attempt using dual-frequency carrier-phase and pseudorange data collected on May 17, 2012, in an article in the September 2012 issue of this magazine. The two GIOVE satellites were subsequently retired.

The four IOV satellites began transmitting navigation messages with valid ephemerides in March, 2013, and this paved the way for the first real-time single-frequency pseudorange Galileo position fix using the broadcast messages on the morning of March 12 at the Navigation Laboratory of the European Space Research and Technology Centre in Noordwijk, the Netherlands. The position fix included compensation for the effect of the ionosphere on the Galileo signals.

The signals from GNSS satellites travel through the ionosphere on their way to receivers on or near the Earth’s surface. The free electrons populating this region of the atmosphere affect the propagation of the signals, changing their speed and direction of travel. This results in a delay in the arrival of the modulated components of the signals (from which pseudorange measurements are made) and an advance in the phases of the signals’ carrier waves (affecting carrier-phase measurements). The ionosphere is a dispersive medium for radio signals, so by making measurements simultaneously on two frequencies transmitted by a satellite, most of the effect of the ionosphere can be removed. However, single-frequency devices such as most vehicle navigation and handheld receivers don’t have the luxury of dual-frequency correction. These devices must rely on a single-frequency correction model. The coefficients for such a model are included in the navigation messages transmitted by all GPS satellites. Known as the Ionospheric Correction Algorithm or Klobuchar Algorithm, it removes at least 50 percent of the ionosphere’s effect.

The Galileo satellites also include the parameters of an ionospheric algorithm, called NeQuick G, in their navigation messages. In this month’s column, the Galileo system design team describes the novel European way for modeling the ionosphere for single-frequency users and compares its performance to the current GPS approach.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas. Write to him at lang @ unb.ca.

Radiowave propagation of GNSS signals is affected by the Earth’s atmosphere and the characteristics of the local environment surrounding the receiver. GNSS systems are based on the broadcasting of radiowave ranging signals in the microwave domain (mainly in the so-called L-band, although some new systems like the Indian Regional Navigation Satellite System are also expected to broadcast in the S-band). These electromagnetic signals may suffer from a number of impairments as they propagate from a satellite to a receiver. In considering these effects, we can divide the Earth’s atmosphere into two parts: the electrically neutral atmosphere (primarily the lowest part, the troposphere), whose main effect is a group delay on the navigation signal due to water vapor and the gas components of the dry air, which, for microwave frequencies, is non-dispersive (independent of frequency); and the ionosphere, the ionized part of the atmosphere. The local environment may affect the navigation signal in various ways, too, such as signal fading or complete signal blockage by vegetation or obstacles such as buildings, and multipath, where the signal is broadened in the time and frequency domains due to reflections and diffraction by surrounding objects. In this article, we will discuss the effect of the ionosphere on GNSS signals and how it is being dealt with by the Galileo satellite navigation system.

The ionosphere owes its existence to solar radiation, primarily extreme ultraviolet light. The radiation ionizes the atoms and molecules in the upper atmosphere at heights of less than a hundred kilometers to a few kilometers above the Earth’s surface, producing a sea of ions and free electrons (collectively known as a plasma). This region is responsible for a number of dispersive (frequency-dependent) effects on navigation signals. Chief among these is a persistent delay of the pseudorandom noise (PRN) ranging codes (and the advance of the phase of the underlying carrier waves), thereby introducing positioning and timing errors if not compensated for. Signals are also susceptible to scintillations — rapid variations of amplitude and/or phase of the signals due to diffraction and refraction caused by plasma irregularities. Furthermore, the ionosphere can bend the signal path, resulting in a slightly longer path than the straight path, and rotate the polarization of the signal.

The ionospheric refractive index (the ratio of the speed of propagation of electromagnetic waves in a vacuum to the speed of their propagation in a medium) is related to the number of free electrons along the propagation path. For this purpose, the total electron content (TEC) is defined as the electron density in a cross-section of 1 square meter, integrated along a slant (or vertical) path between two points (such as a satellite and a receiver). It is often expressed in TEC units (TECU) where 1 TECU = 10¹⁶ electrons per meter squared = 0.1624 meters of delay at the GPS L1 frequency. According to the electron density, the mechanisms responsible for such ionization, and the dynamics, the ionosphere is usually sub-classified in layers of different characteristics: D, E, F1, and F2, with the latter largely responsible for the ionospheric effects on GNSS.

All of the propagation effects due to the ionosphere depend on a number of factors such as time of the day, location, season, and solar activity. There is also an interaction between solar activity, the ionosphere, and the Earth’s magnetic field, which, at times, can result in a significant disturbance of the ionosphere, as happens during geomagnetic storms. On a long timescale, solar activity follows a periodic, approximately 11-year, cycle. And spatially, the behavior of the ionosphere can be broadly classified into four main regions: the equatorial anomaly regions, located at around ±15-20º on either side of the magnetic equator, usually presenting the largest TEC values; mid-latitude regions, where the daytime TEC values are usually less than half the values found in the equatorial anomaly regions; and the auroral and polar regions, which present moderate TEC values but with larger variability than at mid-latitudes due to the characteristics of the geomagnetic field.

If we ignore some smaller, higher-order terms, the ionospheric group delay (the delay of the “group” of waves making up the PRN ranging code modulations) may be expressed in meters as 40.3 sTEC / f², where sTEC is slant TEC in electrons per meter squared, calculated along the straight propagation path between receiver and satellite, and f is the carrier frequency in hertz. This effect introduces ranging errors of several meters if not corrected. The higher order terms usually account for differences at the millimeter level (rising to centimeter level during extreme ionospheric disturbances) and may be safely neglected for code ranging. The effect on the carrier phase has the same magnitude as the code delay, but of opposite sign, meaning that the carrier phase is advanced while propagating through the ionosphere. Since the group delay is dispersive, its effect can be mitigated using linear combinations of signals at two separate frequencies.

For single-frequency receivers, GNSSes often rely on correction models driven by broadcast data. For example, with GPS, the Ionospheric Correction Algorithm (ICA, also known as the Klobuchar algorithm) uses eight broadcast coefficients to describe the ionosphere, which is represented as a two-dimensional thin-shell model (the vTEC is assumed to be concentrated in a two-dimensional shell at a given height, relying on an analytical mapping or obliquity function to convert between vTEC and sTEC depending on the elevation angle of the received signal). This model is very efficient in terms of computational complexity, and it usually removes more than 50 percent of the ionospheric error, particularly at mid-latitudes.

Galileo and NeQuick G

Galileo provides dual-frequency services able to mitigate the effects of the ionosphere, but also services to single-frequency users. For a Galileo single-frequency receiver, an algorithm has been developed based on an adaptation of the NeQuick electron density model.

With the launch of the Galileo In-Orbit Validation (IOV) satellites and the initial navigation message broadcast, for the first time the end-to-end performance of the single-frequency correction algorithm for Galileo could be analyzed. The objective of the IOV phase was to launch the first four operational Galileo satellites and to deploy the first version of a completely new ground segment. During this phase, the European Space Agency (ESA) needed to validate — in the operational environment — all space, ground, and user components and their interfaces, prior to full system deployment, including the single-frequency correction algorithm performance starting from April 2013. Results were obtained for the period up to March 2014, coinciding with the maximum of solar cycle 24 and including three equinoxes with increased solar activity. In this article, we present performance results showing that the algorithm is capable of correcting more than 70 percent of the ionospheric group delay error under nominal ionospheric conditions, using only the reduced Galileo infrastructure during IOV (four satellites and a partial set of the Galileo sensor or monitoring stations).

The Algorithm. The Galileo single-frequency correction algorithm is based on an adaptation of the three-dimensional NeQuick electron density model, driven by an effective ionization level calculated with three broadcast ionospheric coefficients.

The original NeQuick model is a three-dimensional and time-dependent ionospheric electron density model based on an empirical climatological representation of the ionosphere, which predicts monthly mean electron density from analytical profiles, depending on solar-activity-related input values: sunspot number or solar flux, month, geographic latitude and longitude, height and UT. It allows us to calculate the TEC through numerical integration of electron density along a path between a beginning and an end point crossing the ionosphere. As an example, a global vTEC map obtained with NeQuick is illustrated in FIGURE 1. The first version of this model (NeQuick1) was incorporated into a previous version of the International Telecommunication Union (ITU) recommendation ITU-R P.532 for TEC estimation in radiowave propagation predictions. Researchers have continued development of the model with updated formulations, and version NeQuick2 is the one currently recommended by the ITU.

FIGURE 1. Global vTEC map obtained with the NeQuick electron density model for a sunspot number of 150 at 13h UT in the month of April (grid resolution 2.5 degrees × 2.5 degrees).

The NeQuick model has been adapted for Galileo single-frequency ionospheric corrections (for convenience, the Galileo version is known as NeQuick G) in order to derive real-time predictions based a single input parameter, Az, which is determined using three coefficients broadcast in the navigation message. The three coefficients are used in a second-degree polynomial as a function of the modified dip latitude (MODIP) of the receiver, to determine Az, which replaces the solar flux input parameter of the parent NeQuick model, with the following equation:

(1)

where a_i0-2 are the three broadcast coefficients. MODIP is expressed in degrees. A grid table of MODIP values versus geographical location is provided together with the algorithm. A map showing five different MODIP regions is presented in FIGURE 2, each region usually presenting different behavior.

FIGURE 2. MODIP regions. Contours are modified dip latitudes.

The performance of the Galileo single-frequency ionospheric algorithm, designed to reach a correction capability of at least 70 percent of the ionospheric code delay, had been assessed in the past using GPS data only and using GPS plus Galileo In-Orbit Validation Element satellite data for an offline estimation of the broadcast parameters.

Since the first successful autonomous real-time Galileo-based position fix on March 12, 2013, the Galileo navigation messages have been broadcast by the four IOV spacecraft to the external user community, including the ionospheric broadcast parameters determined with IOV-only observations.

Experiment Period and Performance Indicators

To analyze the performance of the single-frequency ionospheric correction, a number of performance indicators were used:
- The root-mean-square (RMS) error of the ionospheric model in meters of L1 code delay, for one station and one day.
- The relative correction capability, expressed as an RMS percentage, defined as:
(2)
where STEC_ref is the reference STEC and STEC_NeQuickG is the STEC obtained with the Galileo correction model. The factor 66 is used to avoid the fact that small absolute errors, which are relatively large due to small reference values, inflate the correction capability; it is linked to a target correction of 70 percent with a minimum absolute threshold of 20 TECU (30 percent of 66 TECU is about 20 TECU).

Performance verification has been assessed for the period from April 2013 to March 2014, which includes the secondary peak of the current solar maximum. The Galileo broadcast data used for this test are the Az coefficients broadcast by the four Galileo IOV satellites. It is important to remember that during the period of this assessment, the IOV infrastructure was reduced with respect to the target full operational capability, including the generation of the ionospheric parameters: four IOV satellites (no other GNSS satellites were used in the estimation) and a reduced number of monitoring stations.

Since the ionospheric correction performance assessment can be done independently of the Galileo signals and analysis of performance is preferred over independent data and locations, reference STEC estimated using dual-frequency observables from GPS at stations from the International GNSS Service (IGS), distributed around the world, were selected for the correction capability performance assessment. This resulted in observations of six to nine satellites for any epoch and with more than 120 stations per day, which assured good global coverage for the test. Performance has been computed individually for each set of broadcast parameters. For this aspect of ionospheric correction assessment, the differences between GPS and full constellation Galileo geometries are considered to be negligible.

As a reference for comparative purposes, for some cases the results have been compared to those obtained with the GPS ICA correction model using the broadcast parameters from GPS satellites.

The reference ionosphere STEC values were computed using dual-frequency carrier-phase GPS observables from IGS stations at a sampling rate of 300 seconds, and using IGS final global ionospheric maps (GIMs) to level the geometry-free combination of carrier phases. In this context, the IGS GIMs are employed to align the geometry-free or ionospheric combination, LI, to compute the ambiguity term (B_I) for each satellite-to-receiver arc:
(3)

where LI represents the linear combination between signals at frequencies f₁ and f₂; is the ionospheric delay in meters of LI; and B_I is composed of several terms: station and satellite phase inter-frequency biases ( and respectively), LI phase ambiguity (λ₁N₁^j – λ₂N₂^j), phase wind-up, multipath, and noise. And i corresponds to the station and j to the satellite.

Then, in order to compute the corresponding B_I term for each satellite-receiver continuous arc, the sTEC prediction of the GIM (sTEC_{GIM_map}) is computed for each satellite ionospheric pierce point, and then the average is computed as follows:
(4)

where the indices i, j, and α correspond to the receiver, satellite, and arc indicator respectively, and the average is performed over the corresponding continuous (no cycle slips) arc (α) of data. is estimated following the mapping function and the procedures to interpolate in space and time recommended by IGS for GIM maps represented in ionosphere-exchange (IONEX) format.

With this estimation, the aligned STEC can be obtained as:
(5)

which is the STEC used as an accurate sTEC estimation or “truth” reference value.

Results

The first analysis that we performed was the daily RMS error and correction capability for all stations. Most days have shown very promising performance. To see different levels of performance, results for one “bad” day and one typical “good” day, in the period of experimentation, are presented in FIGURE 3. It is observed that even for the “bad” day, the correction capability is above 70 percent, except for some stations in the equatorial regions. This performance is exceeded significantly for the “good” day, with RMS residual ionospheric errors below 1.5 meters for L1 even at low latitudes.

FIGURE 3a. Performance of the Galileo single-frequency ionospheric correction when using the E11 satellite broadcast, “bad day” RMS error in meters of L1.

FIGURE 3b. Performance of the Galileo single-frequency ionospheric correction when using the E11 satellite broadcast, “good day” RMS error in meters of L1.

FIGURE 3c. Performance of the Galileo single-frequency ionospheric correction when using the E11 satellite broadcast, “good day” correction capability in percent.

FIGURE 3d. Performance of the Galileo single-frequency ionospheric correction when using the E11 satellite broadcast, “good day” correction capability in percent.

The evolution of the RMS residual error both for Galileo NeQuick G and GPS ICA from April 2013 to March 2014 are presented in FIGURE 4. In this figure, ionospheric activity at the equinoxes is clearly observed in the degradation of performance, and the influence of increased solar activity from October 2013 to March 2014 is also evident.

FIGURE 4. Global daily RMS ionospheric residual error in meters of L1 after correction with Galileo NeQuick G (red) and GPS ICA (blue) from April 2013 to March 2014.

The residual error of the Galileo correction model is already at the level of the expected capability for the full constellation. It also shows better performance as compared to the GPS ICA model, especially at equatorial latitudes.

The level of correction capability for each station for the Galileo NeQuick G model and the GPS ICA model are presented in FIGURE 5 for a quiet day in May 2013 and an active day during the spring equinox in 2014.

FIGURE 5a. RMS correction capability (percent, with a lower bound of 20 TECU) of Galileo NeQuick G correction models for day 127, 2013.

FIGURE 5b. RMS correction capability (percent, with a lower bound of 20 TECU) of GPS ICA correction models for day 127, 2013.

FIGURE 5c. RMS correction capability (percent, with a lower bound of 20 TECU) of Galileo NeQuick G correction models for day 80, 2014.

FIGURE 5d. RMS correction capability (percent, with a lower bound of 20 TECU) of GPS ICA (right) correction models for day 80, 2014.

Effect in the Positioning Domain. We have performed two analyses to assess the correction performance in the positioning domain: one using GPS observables and one with Galileo-only observables. In both cases, we used three ionospheric delay mitigation methods: the dual-frequency ionosphere-free combination, the single-frequency GPS ICA correction algorithm, and the single-frequency Galileo NeQuick G correction algorithm.

The performance of the correction algorithm in the positioning domain using GPS observables was performed with data from two stations: Noordwijk in The Netherlands (a mid- to high-latitude station) and Malindi in Kenya (a low-latitude station) for the day of year (doy) 172 of 2013. Results are presented in FIGURES 6 and 7 showing good performance of the NeQuick G correction, in particular at low latitude. The results do not include code smoothing neither for single-frequency nor dual-frequency positioning. In the results, it may be observed that, as expected, the noise level for single-frequency positioning is much lower than that of ionosphere-free, but a higher bias may be present (the residual mean ionospheric error).

FIGURE 6a. Horizontal GPS positioning error on L1 using single-frequency NeQuick G correction (blue), L1 and GPS ICA (red) and dual-frequency ionosphere-free (green) for mid-latitude station in Noordwijk (doy 172, 2013).

FIGURE 6b. Vertical GPS positioning error on L1 using single-frequency NeQuick G correction (blue), L1 and GPS ICA (red) and dual-frequency ionosphere-free (green) for mid-latitude station in Noordwijk (doy 172, 2013).

FIGURE 7a. Horizontal GPS positioning error on L1 and single-frequency NeQuick G correction (blue), L1 and GPS ICA (red) and dual-frequency ionosphere-free (green) for low-latitude station in Malindi (doy 172, 2013).

FIGURE 7b. Vertical GPS positioning error on L1 and single-frequency NeQuick G correction (blue), L1 and GPS ICA (red) and dual-frequency ionosphere-free (green) for low-latitude station in Malindi (doy 172, 2013).

Positioning domain analysis with Galileo-only observations using the four Galileo IOV satellites, and applying the NeQuick G correction, was evaluated for a station in Washington, D.C., for doy 245, 2013, including E1-only, E5a-only, and dual-frequency E1-E5a ionosphere-free observations. (E1 is centered at the GPS L1 frequency, while E5a is centered at the GPS L5 frequency.) These results are presented in FIGURE 8. The single-frequency positioning performance is considered promising considering the limited number of satellites.

FIGURE 8a. Horizontal Galileo IOV positioning error on E1 and single-frequency NeQuick G correction (blue), E5a and single-frequency NeQuick G correction (red) and dual-frequency E1-E5a ionosphere-free (green) for mid-latitude station in Washington (doy 245, 2013).

FIGURE 8b. Vertical Galileo IOV positioning error on E1 and single-frequency NeQuick G correction (blue), E5a and single-frequency NeQuick G correction (red) and dual-frequency E1-E5a ionosphere-free (green) for mid-latitude station in Washington (doy 245, 2013).

Conclusions

The performance of the Galileo single-frequency ionospheric correction algorithm, based on NeQuick G, was evaluated using the broadcast navigation messages from the four Galileo IOV satellites, both in correction capability and in the positioning domain for the period April 2013 to March 2014. Despite the reduced infrastructure (broadcast ionospheric parameters estimated using only the IOV satellites at a limited number of monitoring stations), the performance shows promising results, in particular for low-latitude regions where the ionosphere is more problematic and, as expected, it has been confirmed that the correction performance is correlated with solar activity.

Acknowledgments

The NeQuick electron density model was developed by the Abdus Salam International Center of Theoretical Physics in Trieste, Italy, and the University of Graz in Austria. The adaptation of NeQuick for the Galileo single-frequency ionospheric correction algorithm (NeQuick G) was performed by ESA and involved the original developers of NeQuick and other European ionospheric scientists under various ESA projects.

Note to Manufacturers

The publication of the NeQuick G model and the Galileo single-frequency correction algorithm is under preparation for public release by the European Commission.

ROBERTO PRIETO-CERDEIRA is a propagation engineer in the European Space Agency (ESA) at the European Space Research and Technology Centre (ESTEC) in Noordwijk, The Netherlands, responsible for the activities related to radiowave propagation for GNSS and satellite mobile communications.

RAUL ORUS-PEREZ is a propagation engineer at ESTEC, working on activities related to radiowave propagation in the troposphere and ionosphere for GNSS and other ESA projects.

EDWARD BREEUWER is the system integration and verification manager in the Galileo Project Office at ESTEC, responsible for the organization and coordination of all testing activities at the system level. He had overall responsibility for the IOV test campaign.

RAFAEL LUCAS-RODRIGUEZ is the Galileo Services Engineering Manager for the Galileo project at ESTEC.

MARCO FALCONE is the System Manager in the Galileo Project Office at ESTEC.

FURTHER READING

• Development of NeQuick Ionospheric Model

“A New Version of the NeQuick Ionosphere Electron Density Model” by B. Nava, P. Coïsson, and S.M. Radicella in Journal of Atmospheric and Solar-Terrestrial Physics, Vol. 70, No. 15, December 2008, pp. 1856–1862, doi: 10.1016/j.jastp.2008.01.015.

“A Family of Ionospheric Models for Different Uses” by G. Hochegger, B. Nava, S.M. Radicella, and R. Leitinger in Physics and Chemistry of the Earth, Part C: Solar, Terrestrial & Planetary Science, Vol. 25, No. 4, 2000, pp. 307–310, doi : 10.1016/S1464-1917(00)00022-2.

“An Analytical Model of the Electron Density Profile in the Ionosphere” by G. Di Giovanni and S.M. Radicella in Advances in Space Research, Vol. 10, No. 11, 1990, pp. 27–30, doi: 10.1016/0273-1177(90)90301-F.

• Evaluation of the Galileo Single-Frequency Ionospheric Model

“Assessment of NeQuick Ionospheric Model for Galileo Single-Frequency Users” by A. Angrisano, S. Gaglione, C. Giola, M. Massaro, and U. Robustelli in Acta Geophysica, Vol. 61, No. 6, December 2013, pp. 1457–1476, doi: 10.2478/s11600-013-0116-2.

Ionosphere Modelling for Galileo Single Frequency Users by B. Bidaine, Ph.D. thesis, Université de Liège, Liège, Belgium, October 2012.

“GIOVE-A Experimentation Campaign: Ionospheric Related Data Analysis” by R. Orus and R. Prieto-Cerdeira in Proceedings of NAVITEC 2008, the 4th ESA Workshop on Satellite Navigation User Equipment Technologies: GNSS User Technologies in the Sensor Fusion Era, Noordwijk, The Netherlands, December 10–12, 2008.

“Assessment of the Ionospheric Correction Algorithm for GALILEO Single Frequency Receivers” by R. Prieto-Cerdeira, R. Orus, and B. Arbesser-Rastburg in Proceedings of NAVITEC 2006, the 3rd ESA Workshop on Satellite Navigation User Equipment Technologies, Noordwijk, The Netherlands, December 11–13, 2006.

“Advanced Ionospheric Modelling for GNSS Single Frequency Users” by M.A Aragón Ángel and F. Amarillo Fernández in the Proceedings of PLANS 2006, the Institute of Electrical and Electronics Engineers / Institute of Navigation Position, Location and Navigation Symposium, San Diego, California, April 24–27, 2006, pp. 110–120, doi: 10.1109/PLANS.2006.1650594.

• GPS Ionospheric Model

“Ionospheric Time-delay Algorithm for Single-frequency GPS Users” by J.A. Klobuchar in IEEE Transactions on Aerospace and Electronic Systems, Vol. AES-23, No. 3, May 1987, pp. 325–331, doi: 10.1109/TAES.1987.310829

“Ionospheric Effects on GPS” by J.A. Klobuchar in GPS World, Vol. 2, No. 4, April 1991, pp. 48–51.

• Ionospheric Effects on GNSS

“GPS, the Ionosphere, and the Solar Maximum” by R.B. Langley in GPS World, Vol. 11, No. 7, July 2000, pp. 44–49.

• International GNSS Service Ionosphere Map Exchange Format

IONEX: The IONosphere Map EXchange Format Version 1 by S. Schaer, W. Gurtner, and J. Feltens, February 25, 1998.
June 1, 2014
Innovation: Reducing the Jitters
Chip-scale atomic clock.

How a Chip-Scale Atomic Clock Can Help Mitigate Broadband Interference

Small low-power atomic clocks can enhance the performance of GPS receivers in a number of ways, including enhanced code-acquisition capability that precise long-term timing allows. And, it turns out, such clocks can effectively mitigate wideband radio frequency interference coming from GPS jammers. We learn how in this month’s column.

By Fang-Cheng Chan, Mathieu Joerger, Samer Khanafseh, Boris Pervan, and Ondrej Jakubov

INNOVATION INSIGHTS by Richard Langley

THE GLOBAL POSITIONING SYSTEM is a marvel of science and engineering. It has become so ubiquitous that we are starting to take it for granted. Receivers are everywhere. In our vehicle satnav units, in our smart phones, even in some of our cameras. They are used to monitor the movement of the Earth’s crust, to measure water vapor in the troposphere, and to study the effects of space weather. They allow surveyors to work more efficiently and prevent us from getting lost in the woods. They navigate aircraft and ships, and they help synchronize mobile phone and electricity networks, and precisely time financial transactions.

GPS can do all of this, in large part, because the signals emitted by each satellite are derived from an onboard atomic clock (or, more technically correct, an atomic frequency standard). The signals from all of the satellites in the GPS constellation need to be synchronized to within a certain tolerance so that accurate (conservatively stated as better than 9 meters horizontally and 15 meters vertically, 95% of the time), real-time positioning can be achieved by a receiver using only a crystal oscillator. This requires satellite clocks with excellent long-term stability so that their offsets from the GPS system timescale can be predicted to better than about 24 nanoseconds, 95% of the time. Such a performance level can only be matched by atomic clocks.

The very first atomic clock was built in 1949. It was based on an energy transition of the ammonia molecule. However, it wasn’t very accurate. So scientists turned to a particular energy transition of the cesium atom and by the mid-1950s had built the first cesium clocks. Subsequently, clocks based on energy transitions of the rubidium and hydrogen atoms were also developed. These initial efforts were rather bulky affairs but in the 1960s, commercial rack-mountable cesium and rubidium devices became available. Further development led to both cesium and rubidium clocks being compact and rugged enough that they could be considered for use in GPS satellites. Following successful tests in the precursor Navigation Technology Satellites, the prototype or Block I GPS satellites were launched with two cesium and two rubidium clocks each. Subsequent versions of the GPS satellites have continued to feature a combination of the two kinds of clocks or just rubidium clocks in the case of the Block IIR satellites.

While it is not necessary to use an atomic clock with a GPS receiver for standard positioning and navigation applications, some demanding tasks such as geodetic reference frame monitoring use atomic frequency standards to control the operation of the receivers. These standards are external devices, often rack mounted, connected to the receiver by a coaxial cable—too large to be embedded inside receivers.

But in 2004, scientists demonstrated a chip-scale atomic clock, and by 2011, they had become commercially available. Such small low-power atomic clocks can enhance the performance of GPS receivers in a number of ways, including enhanced code-acquisition capability that precise long-term timing allows. And, it turns out, such clocks can effectively mitigate wideband radio frequency interference coming from GPS jammers. We learn how in this month’s column.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas. Write to him at lang @ unb.ca.

Currently installed Local Area Augmentation System (LAAS) ground receivers have experienced a number of disruptions in GPS signal tracking due to radio frequency interference (RFI). The main sources of RFI were coming from the illegal use of jammers (also known as personal privacy devices [PPD]) inside vehicles driving by the ground installations. Recently, a number of researchers have studied typical properties of popular PPDs found in the market and have concluded that the effect of PPD interference on the GPS signal is nearly equivalent to that of a wideband signal jammer, to which the current GPS signal is most vulnerable. This threat impacts LAAS or any ground-based augmentation system (GBAS) in two ways:
- Continuity degradation — as vehicles with PPDs pass near the GBAS ground antennas, the reference receivers lose lock due to the overwhelming noise power.
- Integrity degradation — the code tracking error will increase when the noise level in the tracking loop increases.
Numerous interference mitigation techniques have been studied for broadband interference. The interference mitigation methods can be separated according to the two fundamental stages of GPS signal tracking: the front-end stage, in which automatic gain control and antenna nulling/beam forming techniques are relevant, and the baseband stage, where code and carrier-tracking loop algorithms and aiding methods are applicable.

In our current work, the baseband strategy and resources that are practically implementable at GBAS ground stations are considered. Among those resources, we focus on using atomic clocks to mitigate broadband GNSS signal interference. For GPS receivers in general, wide tracking loop bandwidths are used to accommodate the change in signal frequencies and phases caused by user dynamics. Unfortunately, wide bandwidths also allow more noise to enter into the tracking loop, which will be problematic when wideband inference exists. The general approach to mitigate wideband interference is to reduce the tracking loop bandwidth. However, a reference receiver employing a temperature-compensated crystal oscillator (TCXO) needs to maintain a minimum loop bandwidth to track the dynamics of the clock itself, even when all other Doppler effects are removed. The poor stability of TCXOs fundamentally limits the potential to reduce the tracking loop bandwidth. This limitation becomes much less constraining when using an atomic clock at the receiver, especially in the static, vibration-free environment of a GBAS ground station.

Integrating atomic clocks with GPS/GNSS receivers is not a new idea. Nevertheless, the practical feasibility of such integration remained difficult until recent advancements in atomic clock technology, such as commercially available compact-size rubidium frequency standards or, more recently, chip-scale atomic clocks (CSACs). Most of the research using atomic clock integrated GPS receivers aims to improve positioning and timing accuracy, enhance navigation system integrity, or coast through short periods of satellite outages. In these applications, the main function of the atomic clock is to improve the degraded system performance caused by bad satellite geometries. As for using narrower tracking loop bandwidths to obtain better noise/jamming-resistant performance, the majority of work in this area has focused on high-dynamic user environments with extra sensor aiding, such as inertial navigation systems, pseudolites, or other external frequency-stable radio signals. These aids alone do not permit reaching the limitation of tracking loop bandwidth reduction since the remaining Doppler shift from user dynamics still needs to be tracked by the tracking loop itself. Our research intends to explore the lower end of the minimum tracking loop bandwidth for static GPS/GNSS receivers using atomic clocks.

High-frequency-stability atomic clocks naturally reduce the minimum required bandwidth for tracking clock errors (since clock phase random variations are much smaller). We have conducted analyses to obtain the theoretical minimum tracking loop bandwidths using clocks of varying quality. Carrier-phase tracking loop performance under deteriorated C/N₀ conditions (that is, during interference) was investigated because it is the most vulnerable to wideband RFI. The limitations on the quality of atomic clocks and on the receiver tracking algorithms (second- or third-order tracking loop bandwidths) to achieve varying degrees of interference suppression at the GBAS reference receivers are explored. The tracking loop bandwidth reductions and interference attenuations that are achievable using different qualities of atomic clocks, including CSACs and commercially available rubidium receiver clocks, are also discussed in this article.

In addition to the theoretical analyses, actual GPS intermediate frequency (IF) signals have been sampled using a GPS radio frequency (RF) frond-end kit, which is capable of utilizing external clock inputs, connected to a commercially available atomic clock. The sampled IF data are fed into a software receiver together with and without simulated wideband interference to evaluate the performance of interference mitigation using atomic clocks. The wideband interference is numerically simulated based on deteriorated C/N₀. The actual tracking errors generated from real IF data are used to validate the system performance predicted by the preceding broadband interference mitigation analyses.

Signal Tracking Loop and Tracking Error

The carrier-phase tracking phase lock loop (PLL) is introduced first to understand the theoretical connection between the carrier-phase tracking errors and the signal noise plus receiver clock phase errors. A simplified PLL is shown in FIGURE 1 with incoming signals set to zero. In the figure, n(s), c(s), and δθ(s) are receiver white noise, clock phase error or clock disturbance, and tracking loop phase error respectively, with s being the Laplace transform parameter. G(s) is the product of the loop filter F(s) and the receiver clock model 1/s.

FIGURE 1. Simplified tracking loop diagram.

From Figure 1, the transfer functions relating the white noise and clock disturbance to the output can be derived as:
(1)

The frequency response of H(s) is complementary to 1-H(s). Therefore, the PLL tracking performance is a trade-off between the noise rejection performance and the clock disturbance tracking performance.

Total PLL errors resulting from different error sources are presented as phase jitter, which is the root-mean-square (RMS) of resulting phase errors. Equation (2) shows the definition of the standard deviation of phase jitter resulting from the error sources considered in this work:
(2)

where , and are standard deviations of receiver white noise, receiver clock errors, and satellite clock error, respectively, for static receivers.

The standard deviation for each of the clock error sources can be evaluated using the frequency response of the corresponding transfer function and power spectral densities (PSDs). The equations to evaluate the phase error from each error source are:
(3)

where S_rx and S_sv are one-sided PSDs for receiver clock and satellite clock, respectively. B_w is the bandwidth of the tracking loop and T_c is the coherent integration time.

Receiver and Satellite Clock Models

In general, the receiver noise can be reasonably assumed to be white noise with constant PSD with magnitude (noise density) of N₀. However, it is not the case for clock errors. The clock frequency error PSD is usually formulated in the form of a power-law equation and has been used to describe the time and frequency behaviors of the random clock errors in a free running clock:

(4)

where s_y(f) represents the PSD of clock frequency errors and is a function of frequency powers.

The clock phase error PSD can be analytically derived from the frequency PSD equation because the phase error is the time integral of the frequency error:
(5)

where f₀ is the nominal clock frequency. The h coefficients of the clock phase error PSD are the product of the h coefficients from the clock frequency error PSD and the nominal frequency.

We have adopted the PSD clock error models in our work to perform tracking loop performance analysis. The PSD of the CSAC is derived from an Allan deviation figure published by the manufacturer and is shown in FIGURE 2. We took three piecewise Allan deviation straight lines, which are slightly conservative, and converted them to a PSD.

FIGURE 2. Allan deviations for chip-scale atomic clock.

Three PSDs of clock error models are listed in TABLE 1, which represent spectrums of the well known TCXO, the CSAC, and a rubidium standard. Phase noise related h₀ and h₁ coefficients in the CSAC model are assumed to be the same as the TCXO because they can’t be obtained from the Allan deviation figure. The rubidium clock phase noises resulting from h₀ and h₁ coefficients are assumed to be two times smaller than those of the TCXO, and the same model is also used as the satellite clock error model in our tracking loop analysis.

TABLE 1. Coefficients of power-law model.

Theoretical Carrier Tracking Loop Performance

Second- and third-order PLLs are used to study the tracking loop performance. The loop filters for each PLL are given by:
(6)

where F₂(s) and F₃(s) are second- and third-order loop filters respectively. Typical coefficients for the second- and third-order loop filters are a₂ = 1.414; w_o,2 = 4×B_w,2× a₂/[(a₂)²+1]; a₃ = 1.1; b₃ = 2.4; w_o,3 = B_w,3/0.7845. B_w,2 and B_w,3 are the second- and third-order tracking loop bandwidths accordingly.

As stated earlier, three error sources are considered for static receivers. Using the clock error models described earlier, the contribution of different error sources to phase jitter is a function of PLL tracking bandwidth. The resulting phase tracking errors from different error sources are evaluated based on Equation (3) and shown in FIGURE 3.

FIGURE 3. Phase error contribution from different error sources.

The third-order PLL performance using 2-, 1-, 0.5- and 0.1-Hz tracking loop bandwidths were analyzed as a function of C/N₀ and are shown in FIGURES 4 and 5. For each selected bandwidth, three different qualities of receiver clocks were analyzed, and a conventional 15-degree performance threshold was adopted. The second-order PLL performs similarly to the third-order PLL. However, the phase jitter tends to be more biased when the tracking loop bandwidth becomes smaller. This phenomenon will be observed later on using signal data for performance validation. Therefore, only the third-order loop performance analysis is shown in Figures 4 and 5. It is obvious from these two figures that the minimum tracking loop bandwidth for a TCXO receiver PLL is about 2 Hz, and the PLL can work properly only while C/N₀ is above 24 dB-Hz.

FIGURE 4 Tracking loop performance analysis for 2- and 1-Hz loop bandwidth.

FIGURE 5. Tracking loop performance analysis for 0.5- and 0.1-Hz loop bandwidth.

As for the receiver using atomic clocks, CSAC and a rubidium frequency standard in our analysis, the PLL bandwidth can be reduced down to at least 0.1 Hz while C/N₀ is above 15 dB-Hz.

Experimental Tracking Loop Performance

Experimental data were collected at Nottingham Scientiﬁc Limited. The experiment was conducted using a GPS/GNSS RF front end with a built-in TCXO clock. The RF front end also has the capability of accepting atomic clock signals through an external clock input connector to which the CSAC (see Photo) was connected during data collection. All data (using the built-in TCXO clock or the CSAC) were sampled at a 26-MHz sampling rate and at a 6.5-MHz IF with 2-MHz front-end bandwidth and four quantization levels.

A MatLab-coded software defined receiver (SDR) was used to process collected IF samples for tracking loop performance validation. TCXO phase jitters resulting from different tracking loop bandwidths are shown in FIGURE 6 for a typical second-order PLL under a nominal C/N₀, which is about 45 dB-Hz. A 45-degree loss-of-lock threshold was adopted (three times larger than the standard deviation threshold used in an earlier performance analysis). In our work, all code tracking delay lock loops (DLLs) are implemented using a second-order loop filter with 20-millisecond coherent integration time and 0.5-Hz loop bandwidth without any aiding. The resulting phase jitters in the figure become biased when the tracking loop bandwidth is reduced. This observed phenomenon implies that a second-order PLL time response cannot track the clock dynamics when the loop bandwidth approaches the minimum loop bandwidth (where loss of lock occurs).

FIGURE 6. Second-order PLL phase jitter using TCXO.

The same IF data was re-processed by the SDR using the third-order PLL with the same range of tracking loop bandwidths. The resulting phase jitters are shown in FIGURES 7 and 8. There is no observable phase jitter bias before the PLLs lose lock in the figures. These results demonstrate that a third-order PLL performs better in terms of capturing the clock dynamics when the tracking loop bandwidth is reduced close to the limitation. Therefore, only the third-order PLL will be considered further.

FIGURE 7. Third-order PLL phase jitter using TCXO.

FIGURE 8. Third-order PLL phase jitter using CSAC.

The performance of the TCXO PLL can be evaluated from the results in Figure 7. It demonstrates that the minimum loop bandwidth is 2 Hz, which is consistent with the previous analysis shown in figure 4. However, the minimum bandwidth using the CSAC is shown to be 0.5 Hz in Figure 8. This result does not meet the performance predicted by the analysis, which shows that the working bandwidth can be reduced to 0.1 Hz.

Analysis and Tracking Performance under PPD Interference

The motivation of our work, as described earlier, is to improve the receiver signal tracking performance under PPD interference, or equivalently, wideband interference. We carried out a simple analysis first to understand how much signal deterioration a GBAS ground receiver could expect. A 13-dBm/MHz PPD currently available on the market was used to analyze the signal deterioration based on the distance between the PPD and the GBAS ground receiver. A simple analysis using a direct-path model shows that noise power roughly 30 dB higher than the nominal noise level (about -202 dBW/Hz) could be experienced by the GBAS ground receiver if the nearest distance is assumed to be 0.5 kilometers. In this case, any wideband interference mitigation method to address PPD interference has to handle C/N₀ as low as 10 to 15 dB-Hz.

Gaussian distributed white noises were simulated and added on top of the original IF samples, then re-quantized to the original four quantization levels to mimic the PPD interference signal condition. A 20-dB higher noise level was simulated to demonstrate the effectiveness of this signal deterioration technique.

The tracking loop performance using the third-order PLL under low C/N₀ conditions was evaluated using the IF sampling and PPD interference simulation technique just described. The evaluation results show that the minimum PLL bandwidth using the TCXO is still 2 Hz. This result is roughly consistent with a previous analysis showing a 24-dB-Hz C/N₀ limitation using 2-Hz tracking bandwidth. The PLL using the CSAC performs better than that using the TCXO, which is expected.

After raising the noise level 5 dB higher to achieve an average of C/N₀ of 18 dB-Hz, phase jitters using the TCXO exceed the threshold at all bandwidths as shown in FIGURE 9. The same magnitude of noise was also added to the CSAC IF samples. The resulting phase jitters are shown in FIGURE 10, which demonstrates that the minimum bandwidth is 1 Hz for this deteriorated signal condition. Any further increase in noise level will result in loss of lock for PLLs using a CSAC at all tracking bandwidths.

FIGURE 9. Phase jitter using TCXO under 18 dB-Hz C/N₀.

FIGURE 10. Phase jitter using CSAC under 18 dB-Hz C/N₀.

Summary and Future Work

We explored a baseband approach for an effective wideband interference mitigation method in this article. We have presented the theoretical analysis and actual data validation to study the possible improvement of the PLL tracking performance under PPD interference, which has been experienced by LAAS ground receivers.

The limitations of reducing PLL tracking loop bandwidths using different qualities of receiver clocks have been analyzed and compared with the experimental results generated by processing IF samples using an SDR. We conclude that the PLL tracking performance using a TCXO is consistent between theoretical prediction and data validation under both nominal and low C/N₀ conditions. However, the PLL tracking performance using the CSAC was not as good as the analysis prediction under both conditions.

In our future work, to understand the reason for the tracking performance inconsistency using the CSAC, we will carefully examine and evaluate the hardware components in line between the external clock input and the IF sampling chip. In this way, we will exclude the clock performance degradation due to any hardware incompatibility.

Other types of high quality clocks, such as extra-low-phase-noise oven-controlled crystal oscillators and low-phase-noise rubidium oscillators, will also be tested to explore the limitation of PLL tracking bandwidth reduction. If the results using other clocks exhibit good consistency between performance analysis and data validation, it is highly possible that the CSAC clock error model mis-represents the available commercial products.

In our future work, we will also consider simulating PPD interference more closely to the real scenario, by adding analog interference signals on top of GPS/GNSS analog signals before taking digital IF samples.

Acknowledgments

The authors would like to thank the Federal Aviation Administration for supporting the work described in this article. Also, the authors would like to extend their thanks to all members of the Illinois Institute of Technology NavLab and to the collaborators from Nottingham Scientific Limited for their insightful advice. This article is based on the paper “Using a Chip-scale Atomic Clock-Aided GPS Receiver for Broadband Interference Mitigation” presented at ION GNSS+ 2013, the 26th International Technical Meeting of the Satellite Division of The Institute of Navigation held in Nashville, Tennessee, September 16–20, 2013.

Manufacturers

The CSAC used in our tests is a Symmetricom Inc., now part of Microsemi Corp. (www.microsemi.com), model SA.45s. We used a Nottingham Scientific Ltd. (www.nsl.eu.com) Stereo GPS/GNSS RF front end with the MatLab-based SoftGNSS 3.0 software from the Danish GPS Center at Aalborg University (gps.aau.dk).

FANG-CHENG CHAN is a senior research associate in the Navigation Laboratory of the Department of Mechanical and Aerospace Engineering at the Illinois Institute of Technology (IIT) in Chicago. He received his Ph.D in mechanical and aerospace engineering from IIT in 2008. He is currently working on GPS receiver integrity for Local Area Augmentation System (LAAS) ground receivers, researching GPS receiver interference detection and mitigation to prevent unintentional jamming using both baseband and antenna array techniques, and developing navigation and fault detection algorithms with a focus on receiver autonomous integrity monitoring or RAIM.

MATHIEU JOERGER obtained a master’s in mechatronics from the National Institute of Applied Sciences in Strasbourg, France, in 2002, and M.S. and Ph.D. degrees in mechanical and aerospace engineering from IIT in 2002 and 2009 respectively. He is the 2009 recipient of the Institute of Navigation Bradford Parkinson award, which honors outstanding graduate students in the field of GNSS. He is a research assistant professor at IIT, working on multi-sensor integration, on sequential fault-detection for multi-constellation navigation systems, and on relative and differential RAIM for shipboard landing of military aircraft.

SAMER KHANAFSEH is a research assistant professor at IIT. He received his M.S. and Ph.D. degrees in aerospace engineering at IIT in 2003 and 2008, respectively. He has been involved in several aviation applications such as autonomous airborne refueling of unmanned air vehicles, autonomous shipboard landing, and ground-based augmentation systems. He was the recipient of the 2011 Institute of Navigation Early Achievement Award for his contributions to the integrity of carrier-phase navigation systems.

BORIS PERVAN is a professor of mechanical and aerospace engineering at IIT, where he conducts research focused on high-integrity satellite navigation systems. Prof. Pervan received his B.S. from the University of Notre Dame, M.S. from the California Institute of Technology, and Ph.D. from Stanford University.

ONDREJ JAKUBOV received his M.Sc. in electrical engineering from the Czech Technical University (CTU) in Prague in 2010. He is a postgraduate student in the CTU Department of Radio Engineering and he also works as a navigation engineer for Nottingham Scientiﬁc Limited in Nottingham, U.K. His research interests include GNSS signal processing algorithms and receiver architectures.

FURTHER READING

• Authors’ Conference Paper

“Performance Analysis and Experimental Validation of Broadband Interference Mitigation Using an Atomic Clock-Aided GPS Receiver” by F.-C. Chan, S. Khanafseh, M. Joerger, B. Pervan and O. Jakubov in the Proceedings of ION GNSS+ 2013, the 26th International Technical Meeting of the Satellite Division of The Institute of Navigation, Nashville, Tennessee, September 16–20, 2013, pp. 1371–1379.

• Chip-Scale Atomic Clocks

“The SA.45s Chip-Scale Atomic Clock–Early Production Statistics” by R. Lutwak in the Proceedings of the 43rd Annual Precise Time and Time Interval (PTTI) Systems and Applications Meeting, Long Beach, California, November 14–17, 2011, pp. 207–219.

“Time for a Better Receiver: Chip-Scale Atomic Frequency References” by J. Kitching in GPS World, Vol. 18, No. 11, November 2007, pp. 52–57.

“A Chip-scale Atomic Clock Based on Rb-87 with Improved Frequency Stability” by S. Knappe, P.D.D. Schwindt, V. Shah, L. Hollberg, J. Kitching, L. Liew, and J. Moreland in Optics Express, Vol. 13, No. 4, 2005, pp. 1249–1253, doi: 10.1364/OPEX.13.001249.

• Atomic Clocks and GNSS Receivers

“Three Satellite Navigation in an Urban Canyon Using a Chip-scale Atomic Clock” by R. Ramlall, J. Streter, and J.F. Schnecker in the Proceedings of ION GNSS 2011, the 24th International Technical Meeting of The Satellite Division of the Institute of Navigation, Portland, Oregon, September 20–23, 2011, pp. 2937–2945.

“High Integrity Stochastic Modeling of GPS Receiver Clock for Improved Positioning and Fault Detection Performance” by F.-C. Chan, M. Joerger, and B. Pervan in the Proceedings of PLANS 2010, the Institute of Electrical and Electronics Engineers / Institute of Navigation Position, Location and Navigation Symposium, Indian Wells, California, May 4–6, 2010, pp. 1245–1257, doi: 10.1109/PLANS.2010.5507340.

“Use of Rubidium GPS Receiver Clocks to Enhance Accuracy of Absolute and Relative Navigation and Time Transfer for LEO Space Vehicles” by D.B. Cox in the Proceedings of ION GNSS 2007, the 20th International Technical Meeting of the Satellite Division of The Institute of Navigation, Fort Worth, Texas, September 25–28, 2007, pp. 2442–2447.

• Clock Stability

“Signal Tracking,” Chapter 12 in Global Positioning System: Signals, Measurements, and Performance, Revised Second Edition by P. Misra and P. Enge. Published by Ganga-Jamuna Press, Lincoln, Massachusetts, 2011.

“Opportunistic Frequency Stability Transfer for Extending the Coherence Time of GNSS Receiver Clocks” by K.D Wesson, K.M. Pesyna, Jr., J.A. Bhatti, and T.E. Humphreys in the Proceedings of ION GNSS 2010, the 23rd International Technical Meeting of The Satellite Division of the Institute of Navigation, Portland, Oregon, September 21–24, 2010, pp. 2937–2945.

“Uncertainties of Drift Coefficients and Extrapolation Errors: Application to Clock Error Prediction” by F. Vernotte, J. Delporte, M. Brunet, and T. Tournier in Metrologia, Vol. 38, No. 4, 2001, pp. 325–342, doi: 10.1088/0026-1394/38/4/6.

• Tracking Loop Filters and Inertial Navigation System Integration

“Kalman Filter Design Strategies for Code Tracking Loop in Ultra-Tight GPS/INS/PL Integration” by D. Li and J. Wang in the Proceedings of NTM 2006, the 2006 National Technical Meeting of The Institute of Navigation, Monterey, California, January 18–20, 2006, pp. 984–992.

“Satellite Signal Acquisition, Tracking, and Data Demodulation,” Chapter 5 in Understanding GPS: Principles and Applications, Second Edition, E.D. Kaplan and C.J. Hegarty, Editors. Published by Artech House, Norwood, Massachusetts, 2006.

“GPS and Inertial Integration”, Chapter 7 in Global Position System: Theory and Applications, Vol. 2, by R.L. Greenspan. Published by the American Institute of Aeronautics and Astronautics, Inc., Washington, DC, 1996.

• GNSS Jamming

“Know Your Enemy: Signal Characteristics of Civil GPS Jammers” by R.H. Mitch, R.C. Dougherty, M.L. Psiaki, S.P. Powell, B.W. O’Hanlon, J.A. Bhatti, and T.E. Humphreys in GPS World, Vol. 23, No. 1, January 2012, pp. 64–72.

“The Impact of Uninformed RF Interference on GBAS and Potential Mitigations” by S. Pullen, G. Gao, C. Tedeschi, and J. Warburton in the Proceedings of ION GNSS 2012, the 25th International Technical Meeting of the Satellite Division of The Institute of Navigation, Nashville, Tennessee, September 17–21, 2012, pp. 780–789.

“Survey of In-Car Jammers-Analysis and Modeling of the RF Signals and IF Samples (Suitable for Active Signal Cancelation)” by T. Kraus, R. Bauernfeind, and B. Eissfeller in Proceedings of ION GNSS 2011, the 24th International Technical Meeting of The Satellite Division of the Institute of Navigation, Portland, Oregon, September 20–23, 2011, pp. 430–435.
May 1, 2014
Innovation: Ground-Based Augmentation
Combining Galileo with GPS and GLONASS

By Mirko Stanisak, Mark Bitter, and Thomas Feuerle

INNOVATION INSIGHTS by Richard Langley

GPS = SAFER FLIGHT. While reviewing material for an article celebrating the 25th anniversary of the launch in February 1989 of the first Block II or operational GPS satellite, I was yet again annoyed by many articles on the Web stating that GPS only became available for civil use after the launch of this satellite. Some sources get closer to the truth when they say that GPS was opened for civil use in 1983, following the shoot-down of the Korean Airlines Flight 007. In fact, GPS was designed to serve the needs of both the military and civil communities from the outset. A government memo from April 1973 clearly states: “Civil user needs should be considered in the design of the spaceborne equipment.”

One of the first demonstrations of the use of GPS for aircraft navigation occurred in July 1983, when a Sabreliner business jet was flown in stages from Cedar Rapids, Iowa, to the Paris Air Show, flying only when a sufficient number of the experimental or Block I satellites were in view. The first standalone GPS receivers certified for aviation use (with Receiver Autonomous Integrity Monitoring or RAIM) became available by the mid-1990s. But already the Federal Aviation Administration had been looking into the development of a system to provide higher accuracies and better integrity than that afforded by standalone receivers. In 1994, the FAA announced the development of the Wide Area Augmentation System, its brand of a system generically known as satellite-based augmentation. Geostationary satellites transmit corrections and integrity information to GPS receivers, permitting GPS use for en route navigation all the way down to traditional Category I approach and landing. CAT I approaches can be flown down to a decision height of 61 meters (200 feet). WAAS was declared operational on July 10, 2003, but enhancements to the system continue. Japan, Europe, and India also have operational SBAS based on GPS.

Ground-based GPS augmentation was first developed for maritime applications with the U.S. Coast Guard’s low-frequency system coming on line in the mid-1990s. Also in the mid-1990s, the FAA began the development of the Local Area Augmentation System, generically known as a ground-based augmentation system (GBAS), to provide aircraft with approach and landing capabilities from CAT I down through CAT II (30-meter or 100-foot decision height) and CAT III (no decision height but certain visual range minima) using a VHF datalink. Initial CAT I systems are being operated at Bremen, Germany, and at Newark Liberty International Airport and Houston George Bush Intercontinental Airport.

While a GPS-based GBAS will definitely offer improved navigation services for aircraft, might these services be even better if the systems were to use satellites from other constellations besides GPS? In this month’s column, we look at a straw-man concept for modifying the GBAS protocols to accommodate multiple constellations and the results of preliminary tests using GPS, GLONASS, and Galileo simultaneously.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas. Write to him at lang @ unb.ca.

Ever since the declaration of Full Operational Capability (FOC) of the U.S. Global Positioning System in April 1995, GPS has dominated satellite navigation, especially in aviation applications. By contrast, the Russian GLONASS system cannot be used in western aviation because no approval guidelines exist for GLONASS equipment. Thus GPS has been the de-facto standard in aviation for years.

However, within the last few years, major changes have evolved in the field of GNSS, providing a wide variety of useable satellite navigation systems. The European Union launched its Galileo project, which will provide global multi-frequency services in the near future. China is upgrading its BeiDou system (formerly called Compass) to provide global coverage with more medium-Earth-orbit (MEO) satellites. The operators of GPS and GLONASS have started modernization programs that will enable multi-frequency operations in the future, too. Therefore, a large number of usable satellites and signals from multiple systems will soon be available.

In aviation, almost all phases of flight can be assisted by satellite navigation systems nowadays. The most challenging phase of flight with respect to accuracy, continuity, availability, and integrity is the approach and landing phase. The Ground Based Augmentation System (see FIGURE 1; courtesy of the European Organization for Civil Aviation Equipment) allows precision approaches to be performed using satellite navigation. It uses a VHF data link to broadcast differential GNSS corrections, integrity information, and approach definitions to approaching aircraft. These aircraft combine the differential corrections with their own GNSS measurements, calculate a GBAS-corrected position solution, and determine path deviations based on the selected approach.

FIGURE 1. GBAS principle. (Source: EUROCAE WG 28, ED-114)

From a technical perspective, GBAS can use either GPS or GLONASS for differential corrections. For this, the International Civil Aviation Organization (ICAO) Standards and Recommended Practices (SARPs) include GPS and GLONASS side by side. On the other hand, some standardization documents (for example, those from RTCA) are limited to GPS only, effectively excluding GLONASS from being used in the western world. Nevertheless, Russian GBAS systems provide differential corrections for GPS and GLONASS, and are expected to be certified in Russia in the near future. Additional GNSS such as Galileo or BeiDou are not yet included within these documents, as these systems are not approved for aviation use themselves. This article will focus on how a multi-constellation GBAS with GPS, GLONASS, and Galileo could work.

GBAS installations can provide multiple services for different kinds of operation, based on GNSS L1 corrections only. On the one hand, the differentially corrected positioning service (DCPS) is intended to be a generic service for high accuracy positioning. On the other hand, two different GBAS approach services have been defined. GBAS Approach Service Type C (GAST-C) allows Category I (CAT I) procedures and is already in operation. GAST-D is still under development and will enable precision approaches and landings down to CAT II/III minima once certified. To mitigate all possible hazards, GAST-D will require some additional broadcast messages.

VHF Data Broadcast

The VHF Data Broadcast (VDB) is used to communicate binary GBAS messages to approaching aircraft. It operates in the VHF band (108.025 – 117.975 MHz) and uses time-division multiple access (TDMA) to allow the operation of multiple GBAS ground stations on a single frequency. As shown in FIGURE 2, VDB uses UTC time to have a common time frame. Two frames are transmitted each second, lasting 0.5 seconds each. Within each frame, eight slots with durations of 62.5 milliseconds can be used for transmission. Binary application data is encoded using a differentially encoded eight-phase-shift-keying modulation (D8PSK) and a symbol rate of 10,500 symbols per second. With three bits transmitted per symbol, up to 31,500 bits per second can be transmitted. Each slot can contain up to 222 bytes of binary application data. Usually, only a subset of slots is allocated to a particular ground facility. This way, multiple GBAS ground facilities can share a common VDB frequency.

FIGURE 2. VDB timing structure. (Source: RTCA SC-159, DO-246D)

Within each slot, multiple VDB messages can be transmitted as application data. The coding of information in VDB messages is defined in the RTCA’s GNSS-Based Precision Approach Local Area Augmentation System (LAAS) Signal-in-Space Interface Control Document (ICD) and depends on the VDB message type. (LAAS is the U.S. GBAS.) Currently, message types (MT) 1, 2, 3, 4 and 11 are defined. Figure 2 is derived from this document.

Message Type 1 – MT1. Within VDB Message Type 1, differential corrections based on 100-second smoothing are transmitted. These corrections are required by all GBAS approach services (GAST-C and GAST-D). Aside from the differential corrections, additional information for the first broadcast satellite is transmitted. This includes an ephemeris cyclic redundancy check (CRC), mitigating the effects of wrongly received GNSS navigation data, and the Issue of Data (IOD) flag, indicating the time of applicability for the ephemeris data to be used. To transmit this information for all satellites, the satellite for which differential corrections are transmitted first has to be alternated continuously.

Each MT1 message can contain up to 18 pseudorange- and range-rate corrections for individual satellites. Nevertheless, it is possible to link two consecutive MT1 messages using the Additional Message Flag (AMF). The value of this parameter indicates whether this is a single message (0), or the first (1) or second (3) part of a linked MT1 message. Up to 36 differential corrections can be transmitted using two consecutive VDB time slots with 18 corrections each.

All MT1 measurement blocks must be transmitted at least once per frame. The maximum transmission rate is once per slot for all measurement blocks.

Message Type 2 – MT2. VDB Message Type 2 contains station and integrity parameters such as the coordinates of the reference point to which all differential corrections refer. MT2 messages can include (next to a “core” MT2 message) multiple Additional Data Blocks (ADBs) to transmit information required for different GBAS services. At the moment, the Additional Data Blocks 1, 3, and 4 are defined.

ADB1 contains the maximum distance to the reference point at which the corrections may be used (Dmax) as well as parameters to calculate the remaining risk of incorrect GNSS ephemeris data (Kmd,e). Within ADB3, additional information required for GAST-D is transmitted. ADB4 implements the VDB authentication feature. If this ADB is broadcast by a ground facility, MT2 messages must be transmitted first and contain additional indications about which VDB slots are allocated to the ground facility.

MT2 messages must be transmitted at least each 20th frame, but may be repeated up to once per frame.

Message Type 3 – MT3. The VDB Message Type 3 is a fill message, which is only used in conjunction with the GBAS authentication feature (MT2, ADB4). Among other things, this feature requires a minimum slot occupancy of at least 95 percent. Thus, MT3 messages are broadcast only by ground facilities that support the authentication feature and are completely ignored by airborne GBAS receivers.

Message Type 4 – MT4. With VDB Message Type 4, approach information can be broadcast to approaching aircraft. A pilot can select a specific approach by simply tuning to a given channel number.

Currently, GBAS only uses Instrument Landing System look-alike straight-in approaches called Final Approach Segments (FAS). Each FAS represents one approach. This way, a single GBAS ground facility can provide multiple approaches for all runways of an airport. All approaches must be broadcast at least once per 20 consecutive frames.

Message Type 11 – MT11. The VDB Message Type 11 provides differential corrections in a way very similar to MT1 messages. The main difference is that MT11 corrections are based on 30-second smoothing, which is required for GAST-D service. As for MT1, all MT11 measurement blocks must be transmitted at least once per frame.

Enhancements for GBAS with Galileo

At the moment, the GBAS standardization documents include information on GPS, GLONASS, and SBAS ranging sources. No information on Galileo or other constellations has been added yet. Thus, to include Galileo for GBAS, some Galileo-specific experimental additions to the standards are necessary. These proposed modifications have been made in such a way as to keep as close to the other system standards as possible to preserve consistency. This way, hardly any new functionality is added, but additional satellites can be used. The additional Galileo signals (E5a, E5b, E6) are not used at the moment; however, they might be highly beneficial for multi-frequency applications in the future.

All modifications presented here are purely experimental and will most probably not be exactly the same as those in future standards documents. Nevertheless, they provide a way to test Galileo together with GPS and GLONASS for GBAS on an experimental basis.

Ranging Source ID. The Ranging Source ID uniquely addresses a single satellite. It is used in MT1 and MT11 to transmit the differential corrections and other information for each ranging source. In ICAO Annex 10, Standards and Recommended Practices, the Ranging Source ID is defined for GPS, GLONASS, and SBAS only. To provide Galileo corrections as well, an experimental mapping for Galileo satellites was added; see TABLE 1.

TABLE 1. GBAS Ranging Source IDs.

In this way, up to 36 Galileo satellites can be addressed.

Navigation Data. Galileo provides two different sets of navigation data. The I/NAV data corresponds to the Safety-of-Life (SoL) service and is broadcast on E1 and E5b. The F/NAV data corresponds to the Open Service (OS) and is broadcast on E5a. In order to remain as close as possible to the legacy navigation systems, we selected the I/NAV navigation data for use, as it is broadcast on the E1 frequency and can thus be received with an L1-only GNSS receiver.

The navigation data is primarily used in VDB MT1. For the first transmitted correction in this message, the ephemeris set that shall be used in the aircraft is identified via the Issue of Data (IOD) field. To be consistent with the GPS ephemeris, we used Galileo’s IODnav parameter.

Together with the identification of the navigation data, a CRC parameter is transmitted in MT1 for the first satellite within the differential corrections. This parameter ensures that the receiver as well as the ground facility use identical navigation data for all calculations. The CRC algorithm uses the raw navigation data to generate a distinct CRC value.

For GPS and GLONASS, two ephemeris masks are defined. These masks ensure that only information relevant for GBAS processing are covered by the CRC. For Galileo, a similar mask had to be designed.

Additional Data Blocks in MT2. Within VDB MT2, station parameters and integrity information are transmitted. Some parameters for the over-bounding of possible ephemeris errors are specific to each satellite navigation system.

To extend MT2 to Galileo, parameters for the DCPS, GAST-C, and GAST-D must be added for Galileo. For downward compatibility, these parameters cannot be included in the existing Additional Data Blocks beside the existing parameters. Thus, a new Additional Data Block (ADB5) was defined on an experimental basis. This Additional Data Block is dedicated to Galileo and is structured as shown in TABLE 2. The coding of all values corresponds to the coding of the parameters for the existing systems.

TABLE 2. Additional Data Block 5 in Message Type 2 for Galileo parameters.

Optimized VDB Transmission Scheme

Having available a large number of ranging sources for differential corrections, the VHF VDB is a bottleneck for the transmission of this data. To demonstrate this, we first consider the number of visible satellites that there will be in the future. This leads to construction rules for an optimal VDB transmission scheme, which allows transmitting the maximum number of differential corrections.

Number of Satellites Available. To demonstrate the number of differential corrections enabled by the different systems in the future, we computed the number of visible satellites over a day for a stationary GNSS receiver in Braunschweig, Germany. Even though only four Galileo satellites were in orbit at that time, up to 26 different satellites (GPS, GLONASS, and Galileo) were in view simultaneously. Keeping in mind the preliminary Galileo constellation, it is obvious that more than 30 satellites will be available simultaneously in the future — considering only GPS, GLONASS, and Galileo. Adding BeiDou satellites for GBAS would further boost these numbers.

The broadcast of such a large number of differential corrections is limited by the capacity of the VDB and thus by the number of slots assigned to a GBAS ground facility. The number of assigned slots for a facility should be limited as far as possible to be able to use the same frequency for other GBAS ground facilities. Thus, the available capacity must be used as effectively as possible.

Number of Bytes Required. Each VDB message is framed by a message block header (6 bytes) and the message block CRC (4 bytes).

The length of each message depends on the message type and the amount of information to be transmitted. The resulting length for a message of each type is given in TABLE 3.

TABLE 3. Size of different VDB message types (including message block header and CRC). Variable length message types are dependent on the number of corrections, N.

VDB Constraints. A GBAS ground facility must transmit the VDB data following some constraints. These are:
- MT2 messages (including all Additional Data Blocks required) must be transmitted at least each 20th frame (that is, every 10 seconds).
- If authentication is required, each MT2 message must be transmitted in the first slot assigned to the GBAS ground facility.
- All differential corrections (both MT1 and MT11) must be transmitted at least once in each frame. However, it is possible to split the differential corrections into two adjacent slots using the Additional Message Flags in MT1 and MT11 messages.
- Within each MT1 message, the ephemeris decorrelation parameter (P_eph), the Issue of Data (IOD), and the ephemeris CRC is transmitted for the first satellite in the message. Thus, the first satellite must be alternated in order to broadcast the ephemeris information for all satellites.
- Approach definitions are transmitted in MT4 messages. All MT4 messages must be transmitted within at least each 20th slot.
Based on these constraints, a VDB encoding scheme has been developed, which allows us to fulfill all the requirements listed above while optimizing the number of differential corrections that can be transmitted. Even though it is optimized for GAST-D-like services (including authentication parameters, MT11 messages, and experimental Galileo extensions), it can be used for legacy GAST-C systems, too.

Rules for Optimal VDB Transmission. To fulfill the requirement for the MT2 message to be transmitted first, a complete MT2 message must be transmitted each 20th frame at the beginning of the first slot assigned. If no MT2 message has to be transmitted, an MT4 message is transmitted instead. Thus, all messages are arranged in proper order by three simple rules:
1. MT2 (each 20th frame) or MT4 (otherwise)
2. MT11 (all corrections; can be split into two messages)
3. MT1 (all corrections; can be split into two messages).
Additionally, two more rules must be fulfilled. On the one hand, if supporting the authentication feature, each slot in which the ground facility may transmit VDB data must be filled to at least 95 percent. For this, MT3 null messages may be used to ensure that each slot is filled sufficiently. On the other hand, an additional rule for MT1 messages is necessary if more than three slots are assigned to the GBAS ground facility. In this case, to maximize the number of differential corrections the MT1 messages may be transmitted in the last two assigned slots only. This rule is necessary because the Additional Message Flag is limited to two slots for differential corrections.

Using this transmission scheme, the number of differential corrections is maximized while fulfilling the minimum requirements on the VDB data. Even in case of the maximum number of differential corrections, MT4 approach definitions can still be broadcast. However, in this case, the number of transmittable FAS segments is limited to 19. If more approaches (or different approach types such as Terminal Area Paths (TAPs)) have to be transmitted, the VDB generation scheme must be adapted.

Number of Transmittable Corrections. Using the optimized transmission scheme explained earlier, the number of transmittable corrections can be calculated easily for different numbers of assigned slots for GAST-C as well as for GAST-D services (see TABLE 4).

TABLE 4. Number of differential corrections that can be broadcast.

The exact distribution of VDB messages for the maximum number of differential corrections (18) is shown in FIGURE 3 for an MT1/MT11 configuration and two assigned slots.

FIGURE 3. VDB messages for two slots and 18 satellites (MT1 and MT11).

Experimental Realization of Multi-Constellation GBAS

The experimental GBAS multi-constellation extensions described earlier have been implemented in software for further testing. As these enhancements are purely experimental and might change in the future, we have ensured that these definitions can be changed easily.

Navigation Software. The Institute of Flight Guidance at Technische Universität Braunschweig has been developing an experimental navigation framework for many years. This software, called TriPos, can handle and combine different navigation technologies. TriPos can be used for simulations, post-processing of recorded data, and even for live (online) processing. It is written in C++ and supports various platforms.

The navigation framework can be extended easily. Originally, only GPS was supported within the software, but support for GLONASS and Galileo as well as augmentation systems like SBAS and GBAS were added over the past few years. Additionally, the software handles GNSS data of multiple frequencies internally and can thus be used for multi-constellation and multi-frequency applications. TriPos includes decoders for the binary protocols of most GNSS receivers currently available.

For GBAS research, two components can be simulated using the software. On the one hand, the Ground Facility simulation calculates the differential corrections and provides simulated VDB data. On the other hand, the GBAS receiver simulation emulates the behavior of an airborne GBAS receiver and uses VDB data and GNSS measurements to calculate a GBAS solution. Both simulations can use either recorded data in post-processing or live data for online-processing. This allows complete simulation of GBAS.

Multi-Constellation GBAS Ground Facility Simulation. The GBAS ground facility simulation uses raw binary data from multiple stationary GNSS receivers to calculate binary VDB data. The simulation can be freely configured to process either live or pre-recorded GNSS data. Even though it features all algorithms required by the standards, it does not contain additional monitor algorithms at the moment.

Nevertheless, it can provide a valid VDB signal-in-space (SIS), which can be used by GBAS receivers and simulation tools (such as Eurocontrol’s PEGASUS tool). The ground facility simulation supports legacy GBAS CAT-I (GAST-C) as well as GAST-D (including all additional VDB information required) using GPS and GLONASS. Support for Galileo has been added according to the experimental definitions described earlier. In addition to FAS data blocks, the ground facility simulation is also capable of providing curved approaches using TAP data blocks.

Multi-Constellation Airborne GBAS Receiver Simulation. The GBAS receiver simulation has been used for various GBAS-related projects. It supports GAST-C as well as GAST-D and can be configured flexibly to use GPS, GLONASS, and/or Galileo (using the experimental enhancements as described earlier). For GAST-D, all airborne monitoring algorithms required are present. Thus, the aircraft-specific parameters (for example for the airborne geometry screening) can be configured together with the other parameters.

Flight Trials

The practicability of the multi-constellation GBAS approach has been tested in flight trials. To ensure that all four Galileo satellites were in view and capable of providing valid data during our trials, an orbit prediction tool and the Notice Advisory to Galileo Users (NAGU) service of the European GNSS Service Center (GSC) were used prior to the flight.

The data processing configuration is shown in FIGURE 4 and includes the GBAS simulation components explained earlier. All processing is done in real time while recording all data for later post processing.

FIGURE 4. Schematic data processing for the flight experiments (ground components in orange, airborne components in blue).

Ground Processing. On the ground, two Septentrio AsteRx3 GNSS receivers connected to two roof-top antennas were used. The GNSS receivers were connected to the GBAS ground facility simulation via a network and provided binary GPS, GLONASS, and Galileo raw measurements with an update rate of 2 Hz as well as navigation data. Using this data, the ground facility simulation generated binary VDB data. The GBAS ground facility simulation was configured to generate multi-constellation GAST-D VDB data for a three-slot configuration. All required messages (MT1, MT2 including all required ADBs, MT3, MT4 and MT11) were generated and sent to the telemetry facility via the network.

Telemetry. Official VHF data broadcasts operate in a frequency band between 108 and 118 MHz, which is reserved for authorized aviation applications. However, for our experimental system, an alternative data link was used. The Institute of Flight Guidance operates a full-duplex telemetry system to share data between ground and aircraft. Even though the operating frequencies are different, the telemetry system allows the generated binary VDB data to be transmitted to research aircraft. The airborne telemetry receiver outputs data as if it were a VDB receiver to allow us to switch between a real VDB receiver and the telemetry receiver easily.

Research Aircraft. The Institute of Flight Guidance operates the research aircraft of the Technische Universität Braunschweig. The Dornier Do 128-6 with the call sign D-IBUF (see FIGURE 5) is a twin-engine turboprop aircraft without a pressurized cabin and has been used multiple times for GBAS-related research over the years.

FIGURE 5. Research aircraft D-IBUF (Dornier Do 128-6).

The research aircraft allows us to flexibly integrate experimental equipment for specific flight trials. For the multi-constellation GBAS flights, a JAVAD Delta GNSS receiver (capable of multiple constellations and frequencies), a telemetry receiver, and an experimental cockpit display were installed temporarily.

Airborne Processing. The online GBAS receiver simulator uses GNSS data from the JAVAD Delta GNSS receiver together with the VDB data received via telemetry. The receiver was configured to output raw GPS, GLONASS, and Galileo measurements with an update rate of 10 Hz. The simulator was configured to use this data to calculate a multi-constellation GAST-D solution. Based on the selected approach definition, the resulting information (deviations, distance to threshold, and so on) was displayed in the cockpit using an experimental cockpit display.

Results. The flight test was conducted in the evening of November 6, 2013 (16:52 – 17:58 UTC), at Research Airport Braunschweig (EDVE). We performed five approaches with a 10 nautical mile final segment. The flight path as calculated by the GBAS receiver subsystem is shown in FIGURE 6.

FIGURE 6. Flight trial trajectory. (Map data © OpenStreetMap contributors)

FIGURE 7 shows the number of satellites used for the GBAS receiver simulation, and distinguishes between the different satellite navigation systems used. Up to 22 satellites have been used simultaneously for GBAS processing, including up to 10 GPS satellites, eight GLONASS satellites, and four Galileo satellites.

FIGURE 7. Number of satellites used by the multi-constellation GBAS receiver simulation.

Even though no certified GBAS equipment was used for the flight trials, FIGURE 8 shows the resulting vertical and lateral protection levels (VPL and LPL) of the online multi-constellation GBAS receiver simulation. Both values fluctuate due to the differences between 100- and 30-second smoothing position solutions, which have to be added to the protection levels for GAST-D. Nevertheless, both sets of values remain clearly below the corresponding Alert Limits (FAS Lateral Alarm Limit (FASLAL): 40 meters, FAS Vertical Alarm Limit (FASVAL): 10 meters). A valid GAST-D service was achieved continuously.

FIGURE 8. Vertical and lateral protection levels (VPL and LPL).

FIGURE 9 shows a vertical integrity diagram, commonly known as a Stanford plot, for the integrity of the multi-constellation GBAS simulation. This plot shows the Vertical Protection Level (VPL) as determined by the GBAS receiver simulation against the actual Vertical Position Error (VPE). The Vertical Position Error is a direct measure for the Vertical Navigation System Error (V-NSE). This has been determined using a precise point positioning reference trajectory. Both values are normalized by the current VAL as these values change during the approaches. During the flight, the GBAS online processing ran at a rate of 10 Hz, resulting in 43,670 GAST-D epochs and an availability of 100 percent.

FIGURE 9. Normalized vertical Stanford plot of flight trials (GAST-D using GPS, GLONASS, and Galileo). Color scale indicates number of occurrences.

Of course, these results must not be misinterpreted as a multi-constellation GBAS performance assessment. The ground facility simulation was highly experimental and lacked any kind of long-term analysis. Even the GNSS antennas used do not meet formal requirements. However, aside from a quantitative judgment, these results show the practicability of this multi-constellation GBAS approach on an experimental basis.

Conclusion and Outlook

In this article, experimental extensions to GBAS have been developed to support GPS, GLONASS, and Galileo simultaneously. Based on these extensions, an optimized VDB transmission scheme has been created. In this way, the number of transmittable differential corrections could be maximized. Using flight trials, the multi-constellation GBAS concept has successfully been verified. The experimental airborne GBAS subsystem was able to calculate a valid GBAS solution including GPS, GLONASS, and Galileo satellites continuously.

It has been shown that multi-constellation GBAS is possible from a purely technical perspective. On the other hand, neither operational nor approval aspects for satellite navigation systems other than GPS have been addressed yet. Additionally, further testing would be necessary to ensure the compatibility with legacy GPS-only GBAS equipment. However, in theory, all modifications for Galileo are backward compatible. Nevertheless, it has to be assured that certified GBAS multi-mode receivers only use the GPS part of the VDB data and are not disturbed by additional VDB messages or additional ranging sources, for example. The required tests are planned for the future.

The operational benefit of multi-constellation GBAS systems cannot be foreseen yet. A certification for this will take several years and could only be addressed by the GBAS community after the completion of the GAST-D certification. Most probably, the use of GNSS signals on multiple frequencies could provide a highly improved GBAS service and will allow much more operational benefit. Many of the satellite navigation systems have already introduced additional frequencies, including signals in the protected L5 aviation band. The use of multiple frequencies for satellite navigation in aviation can remove most ionospheric errors effectively and mitigate a major source of uncertainty. Thus, multi-constellation GBAS can just be seen as a preliminary step on the way towards multi-frequency GBAS. The concepts and infrastructure described in this article will serve as a basis for more research in this area.

Acknowledgments

Most of our work on multi-constellation GBAS was done within the research project “Bürgernahes Flugzeug,” which was established in 2009 and is partly funded by the German federal state of Lower Saxony. This is gratefully acknowledged by the authors. Additionally, the authors would like to thank all colleagues involved for constructive discussions and their support. This article is based on the paper “Mulitple Satellite Navigation for the Ground Based Augmentation System” presented at ITM 2014, The Institute of Navigation 2014 International Technical Meeting, held in San Diego, California, January 27-29, 2014.

MIRKO STANISAK is a research assistant at the Institute of Flight Guidance (IFF) at the Technische Universität (TU) Braunschweig in Germany. He received his diploma in mechanical engineering (Dipl.-Ing.) in 2009 from TU Braunschweig.

MARK BITTER holds a Dipl.-Ing. in mechanical engineering from TU Braunschweig and has been employed as a research engineer at TU Braunschweig IFF since 2003.

THOMAS FEUERLE received his Dipl.-Ing. in mechanical engineering in 1997 from TU Braunschweig. He joined the TU Braunschweig IFF in May 1997. Since 2005, he has been the leader of the Air Traffic Management Team at the IFF. In April 2010, he completed his Ph.D. dissertation at TU Braunschweig.

FURTHER READING

• Authors’ Conference Paper

“Multiple Satellite Navigation Systems for the Ground Based Augmentation System,” by M. Stanisak, M. Bitter, and T. Feuerle in Proceedings of ITM 2014, the 2014 International Technical Meeting of The Institute of Navigation, San Diego, California, January 27–29, 2014, pp. 254–264.

• Standards Documents

Aeronautical Communications, Vol. 1, Radio Navigation Aids, Annex 10 to the Convention on International Civil Aviation, International Standards and Recommended Practices, International Civil Aviation Organization, Montreal, Draft Version, May 2010.

GNSS-Based Precision Approach Local Area Augmentation System (LAAS) Signal-In Space Interface Control Document (ICD), DO-246D, RTCA Special Committee 159, Global Positioning Systems, RTCA Inc. Washington, D.C., December 2008.

Minimum Operational Performance Standards for GPS Local Area Augmentation System Airborne Equipment, DO-253C, RTCA Special Committee 159, Global Positioning Systems, RTCA Inc. Washington, D.C., December 2008.

Minimum Operational Performance Specification for Global Navigation Satellite Ground Based Augmentation System Ground Equipment to Support Category I Operations, ED-114, EUROCAE Working Group 28 on Global Navigation Satellite System, European Organisation for Civil Aviation Equipment, Malakoff, France, September 2003.

• GBAS Research and Development

“Conception, Implementation and Validation of a GAST-D Capable Airborne Receiver Simulation” by M. Stanisak, R. Schork, M. Kujawska, T. Feuerle, and P. Hecker in Proceedings of ION GNSS 2012, the 25th International Technical Meeting of the Satellite Division of The Institute of Navigation, Nashville, Tennessee, September 17–21, 2012, pp. 250–257.

“Making the Case for GBAS: Experimental Aircraft Approaches in Germany,” by U. Bestmann, P.M. Schachtebeck, T. Feuerle, and P. Hecker in Inside GNSS, Vol. 1, No. 7, October 2006, pp. 42–45.

“Initial GBAS Experiences in Europe” by A. Lipp, A. Quiles, M. Reche, W. Dunkel, and S. Grand-Perret in Proceedings of ION GNSS 2005, the 18th International Technical Meeting of the Satellite Division of The Institute of Navigation, Long Beach, California, September 13–16, 2005, pp. 2911–2922.

• GPS Use in Aviation

“Aircraft Landings: The GPS Approach,” by G. Dewar in GPS World, Vol. 10, No. 6, June 1999, pp. 68–74.

“GPS in Civil Aviation” by K.D. McDonald in GPS World, Vol. 2, No. 8, September 1991, pp. 52–59.
April 2, 2014
Innovation: A PET Project from Finland
Automating GNSS Receiver Testing

By Sarang Thombre, Jussi Raasakka, Tommi Paakki, Francescantonio Della Rosa, Mikko Valkama, Laura Ruotsalainen, Heidi Kuusniemi, and Jari Nurmi

INNOVATION INSIGHTS by Richard Langley

WE HAVE A CAT. My wife and I do, that is. One with a voracious appetite. She likes to be fed on demand, even at the most inopportune times. Like three o’clock in the morning. No, it doesn’t help to close the bedroom door. Her squeaking (yes, some cats squeak) still wakes us up. I was designated as the one to get up in the night to feed her. Sometimes twice. Each night, every night. That got tiresome (literally) very quickly. Automation came to the rescue. We now have a microprocessor-controlled cat feeder, which rotates food compartments into the feeding position at pre-programmed times. Just fill up one or two of the compartments with “crunchies” before retiring, set the activation time to 3:00 a.m., say, and no more middle-of-the-night squeaking interrupting blissful sleep.

This is just one example of how automation — machines replacing (or supplementing) human activity to perform repetitive, difficult, undesirable, or even humanly-impossible tasks — can affect (and benefit) our everyday lives.

As noted on Wikipedia, two common types of automation are ones that involve feedback control, which is usually continuous and involves making measurements using one or more sensors and computing adjustments to keep the measured variables within a set range, and those that involve sequence control, in which a programmed sequence of discrete operations is performed, often based on system logic. An aircraft autopilot is an example of the former while our cat-feeding machine is an example of the latter. Some systems, such as Earth-orbiting satellites, can involve both types.

Automation applications range from the (now) mundane (such as point-and-shoot cameras, smart phones, home control, and factory assembly lines) to the (now) exotic (such as robots to assist the elderly and the infirm and robots to explore space). Laboratories have also benefited from increasing automation, making rapid clinical and analytical testing, for example, possible.

The testing of GNSS receivers can also benefit from automation. This work typically requires the active participation of humans to initiate, control, monitor, and terminate test cases. These manual operations are often inefficient and inaccurate, rendering the test results unreliable. Furthermore, accessing the internal signals of a receiver at different stages of processing is necessary to pinpoint the exact location of any anomalies. Using traditional black-box testing techniques, it is only possible to test the final outputs of a receiver. In this month’s column, we take a look at an automated test bench for analyzing the overall performance of multi-frequency, multi-constellation GNSS receivers. The system includes a data-capture tool to extract internal process information and controlling software, called the Automated Performance Evaluation Tool or AutoPET, which is able to communicate between all modules of the system for hands-free, one-button-click testing of GNSS receivers. Would my cat appreciate the benefit? Likely not, but GNSS engineers and scientists certainly will.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas.

The prototype GNSS receiver developed at the Department of Electronics and Communications Engineering of Tampere University of Technology (TUT), called TUTGNSS, is now in the performance-testing phase. TUTGNSS is a GPS L1/L5 + Galileo E1/E5a dual-frequency dual-constellation receiver jointly developed by TUT and its international partners under two European Union Framework Programme research grants.

During the manual testing of the receiver, it was noticed that the results were often contaminated with errors due to imprecise time-keeping and inconsistent test environments.

It was also strenuous and time consuming to perform repetitive tests over multiple iterations, with decreasing personnel efficiency as the number of iterations increased. The aforementioned problems led to the results being deemed unreliable and unrepeatable. There was thus a need to innovate and automate the testing process and environment. In addition, there was also the need to study the signals as they flowed through the internal signal processing chain, so that the exact location of anomalies could be detected.

Currently, few solutions are available in the commercial and academic domains, which can perform end-to-end fully automated, yet customizable testing of GNSS receivers. A couple of commercial testing tools were recently unveiled, which claim to perform similar automated testing of GNSS receivers. However, these are not fully customizable by the end-user, having the limitation that they can be used only with their parent company’s proprietary signal simulators. Other commercial automated testing tools are available nowadays. However, they are targeted towards electronic systems other than GNSS receivers. It was due to these reasons that we decided to implement an in-house solution. Consequently, we devised the Automated Performance Evaluation Tool (AutoPET), along with a data capture tool.

AutoPET is implemented completely in software (Qt, with C++) and communicates with the receiver under test (RUT) via RS-232 and a National Marine Electronics Association (NMEA) protocol and with a commercial GNSS signal simulator via an RS-232 link. It handles the GNSS test cases with user-defined iterations and other system settings. AutoPET has already been used for making test runs on the TUTGNSS receiver with positive results. It is possible to initiate the overall testing of the receiver with a single button-click and the results are stored in the computer without any human intervention. Test scenarios currently included in the tool’s library are: time-to-first-fix (TTFF), position accuracy, acquisition sensitivity, tracking sensitivity, and reacquisition time. By changing the scenarios in this library, the tool can be used with different simulator models. Another innovative aspect of AutoPET is that it uses multi-threading to perform the receiver testing. Multiple software processing threads are necessary to keep track of the receiver operations and simulator feeds simultaneously, so that an appropriate interrupt can be generated when the receiver has performed the desired operation. This feature is explained in further detail later on.

Data Capture Tool (dCAP) is a hybrid (software-controlled hardware) entity capable of extracting the user-defined internal process data from the different modules (acquisition, tracking, bit decoding, and so on) of the GNSS RUT and stores it in a computer via a 100-Mbps Ethernet link. The dCAP hardware is independent of the receiver module (although implemented on the same softcore) and operates through minimal interference with the receiver operation. This data can then be post-processed to monitor and record the behavior of the receiver and to investigate any anomalies in its intermediate stages. An experimental version of dCAP has already been used to monitor the carrier-to-noise-density ratio (C/N₀), carrier Doppler, and code delay from the internal tracking channels, and the raw GNSS signals in I/Q format entering the baseband processing unit (BPU) of the TUTGNSS receiver from its radio front end.

The benefits of AutoPET over state-of-art approaches are that it is portable (software platform independent), easy to use, suitable for testing most receivers using a variety of simulators (provided each of them can communicate with the outside world using some form of communication protocol), and its operational parameters are easy to modify through an external configuration file. dCAP is designed specifically for the TUTGNSS receiver; however, it can be easily replicated for most experimental embedded system receivers. Once implemented, dCAP offers a clear view of the internal operation of the receiver by accessing intermediate signals between the input and output terminals. The speed and size of data capture are limited only by the type of Ethernet connection and the size of the internal and external memories. Additional details of AutoPET and dCAP are provided in the next two sections of this article, while the third section describes the application of these tools in testing the GPS L1 operation of the TUTGNSS receiver.

Automated Performance Evaluation Tool

AutoPET is a software program developed using the Qt platform and the C++ language, which communicates between the GNSS receiver, signal simulator, and its associated computer through a remote PC that houses AutoPET. The set-up is shown in FIGURE 1. This figure also denotes the different communication protocols used between the different modules.

FIGURE 1. Block schematic of the AutoPET assembly.

At the center is the GNSS receiver, which accepts RF signals from the GNSS signal simulator. These signals represent signals from the sky in accordance with the scenario loaded in the simulator, and therefore represent unidirectional communication. On the other hand, the receiver communicates with the remote PC housing AutoPET using the NMEA-0183 protocol. This is bidirectional communication, as the receiver continuously updates its status via NMEA messages to AutoPET and, in turn, AutoPET sends a response/control command to the receiver. The receiver sends the $GPGGA NMEA message every second, and through reading this message, AutoPET can determine the current status (acquisition, tracking, position fix, and so on) of the receiver.

The TUTGNSS receiver has the capability to perform a cold start to initiate the next test iteration when commanded by AutoPET. For this purpose, we have designed a simple custom message string, which can be identified by the TUTGNSS receiver as a cold-start command. In response, the receiver sends a custom NMEA message, $GPTXT, which identifies that it has successfully performed a cold start. Performing a cold start involves erasing all a priori navigation-related information from the receiver memory. This includes erasing the ephemeris, almanac, and timing information, and ensuring that all satellite tracking is lost.

AutoPET communicates with the GNSS signal simulator through its controlling computer, called the Sim-PC (which runs the control software for the simulator). This communication is bidirectional using a 100-Mbps Ethernet link. The AutoPET library holds the scenario files, through which it remotely controls the simulator. In turn, the Sim-PC returns responses or error messages in the form of Extensible Markup Language (XML) strings to the AutoPET. The communication between the Sim-PC and the simulator is through its proprietary protocols.

AutoPET makes extensive use of multi-threading. The receiver, AutoPET, and the simulator function autonomously of each other and hence are independently controlled using their own processing threads running in parallel. Examples of some processing threads are:
- Thread 1 – monitors the receiver operation through the received NMEA messages. This thread is responsible for identifying, for example, if the receiver achieves a position fix or if it performs a successful cold start.
- Thread 2 – monitors the simulator through the received XML error messages and response messages from the Sim-PC. It is responsible for identifying, for example, if the simulator scenario is successfully set up or if the satellite signals are turned on and off when demanded by the test case.
- Thread 3 – monitors the internal operation of AutoPET itself to check, for example, if a timer has expired or if the user performs any operation on the GUI during the progress of a test.
Each thread generates an internal software interrupt within AutoPET based on which future course of action has to be dynamically determined.

Later in the article, the application of AutoPET for single-frequency, single-constellation operation and testing of the TUTGNSS receiver is described. However, it can just as easily be applied for more complex, multi-frequency, multi-constellation testing. The scenarios are stored in the library of AutoPET, and they can be easily updated without requiring any changes in the tool itself. On the other hand, the receiver operation needs to be updated to perform position fixes with multiple signals and constellations. If the receiver allows updating of its operation mode using software commands, as is the case in TUTGNSS, these commands can also be included within AutoPET.

In the case of TUTGNSS, two configuration settings control the mode of operation and the manner in which it has to be turned on (cold, warm, or hot start) via a 32-bit control word. Table 1 describes the various options and the digital control word bits corresponding to each option. There are eight possible modes of operation that would require three bits to be uniquely represented. However, we have assigned five bits, to accommodate any planned future increase in operating modes. Similarly, there are three ways to turn on the TUTGNSS receiver, and they can be uniquely represented by two bits. Therefore, out of the 32 available bits, only seven bits are currently utilized. The rest of the bits are in reserve for future use. The mode selection bits are in least significant bit positions of the control word. For example, if the receiver should perform a position fix after a warm start using GPS L1 and Galileo E1 signals, the 32 bit control word would be 00000000_0000000_00000000_00100010. Using this control word at the beginning of every test, AutoPET can be used for a simple single constellation or more advanced multi-constellation testing of the receiver.

TABLE 1. Control words for multi-frequency, multi-constellation testing of TUTGNSS.

Data Capture Tool

The overall set up of dCAP is shown in FIGURE 2. The TUTGNSS receiver consists of the radio front end and the BPU implemented on an Altera Stratix-II development board. This board consists of the NIOS-II softcore embedded processor controlled by the MicroC operating system within a field-programmable gate array (FPGA) board. The hardware is programmed using VHSIC Hardware Description Language (VHDL) and consists of the system entity and a few peripheral entities, such as a phase-locked loop (PLL), which are not shown in the figure for sake of simplicity. The system entity consists of (among others) two software-controlled hardware entities, one for the TUTGNSS receiver BPU and the other for the dCAP server, called CPU-0 and CPU-1 respectively. The Control-PC is responsible for the overall programming of the FPGA board through a USB link. It also holds a Qt-based user interface acting as the dCAP client implementation.

FIGURE 2. Overall block schematic of the dCAP assembly.

The dCAP client (in the Control-PC) establishes an Ethernet connection with the dCAP server (on the FPGA) and requests a user-specified internal data sample. As an example, let us assume the user requests raw I/Q samples input to the TUTGNSS BPU from the radio front end. The dCAP server software communicates with the TUTGNSS software, which in turn allows the dCAP server hardware access to the requested data from the appropriate region of the TUTGNSS hardware, similar to how a signal across a resistor on a dense printed circuit board is viewed by placing oscilloscope probes across it. The only limitation with dCAP is that the user has to predict, in advance, which internal data parameters are of interest and create access points within the correct hardware entities. The dCAP server hardware will connect to the respective access point when demanded by the client.

This data snapshot is first buffered in the local shared memory entity on the FPGA board due to the requirements of speed, size, and time synchronization. The dCAP server software is responsible for transferring this data from the internal memory to the Control-PC through the Ethernet link. The data is stored on the Control-PC hard drive in the form of a *.bin file. Therefore, the size of each data-packet that can be accessed at a time is limited by the size of the FPGA memory entity, while the total data size is limited only by the size of the hard drive of the Control-PC. The speed of data capture is restricted by the maximum speed of the Ethernet link between the dCAP client and server.

In FIGURE 3, the internal operation of the dCAP server is illustrated, assuming that we would like to access the raw samples from the radio front end. The first block that the samples enter inside the TUTGNSS BPU is the baseband converter unit (BCU). This is where the dCAP hardware probes listen in on the signal samples. Through these probes, the signals are diverted to the first-in-first-out (FIFO) data collector on the dCAP server (CPU-1) in addition to their usual route through the further baseband processing blocks of the receiver. After the FIFO collector, the data undergoes clock arbitration, time synchronization, and master-slave synchronization, before being buffered into the on-chip Synchronous Dynamic Random Access Memory (SDRAM), where it waits until the dCAP server transfers it through the Ethernet-based local network to the requesting dCAP client within the Control-PC. In the case where different internal data has to be monitored, the probes simply reorient to the correct access point within the correct hardware entity (for example, to monitor the signal C/N₀, the probes access the tracking loops).

FIGURE 3. Block schematic of an example of the dCAP internal operation.

TUTGNSS Receiver Performance Testing

During the GPS L1 performance testing of the TUTGNSS receiver, the reference receiver position in the simulator was set randomly. Ionosphere and troposphere errors were turned off in the simulator. On average, 100 iterations were performed for each test, and the total duration to complete all tests was two weeks. dCAP was used in monitoring the tracking channels and extracting information such as the C/N₀, carrier Doppler, and code-delay estimates for the satellites being tracked. Access to these parameters enabled testing the acquisition and tracking sensitivity of the TUTGNSS receiver, thus confirming the results of the tests performed using AutoPET.

Acquisition Sensitivity. Acquisition sensitivity for the TUTGNSS receiver was measured to be -141.5 dBm via AutoPET and -141 dBm via dCAP. Each coherent integration interval was 4 milliseconds, and 256 such intervals were integrated non-coherently. Using AutoPET, 100 acquisition iterations were performed at every power level, and the average number of satellites acquired was recorded. It was observed that no satellites were acquired at -142 dBm. The acquisition sensitivity test using dCAP involved extracting the carrier Doppler and code-delay estimates. A successful acquisition was assumed only if the code-delay estimate error was less than ±1 chip (300 meters) and the carrier Doppler estimate error was less than ±150 Hz. Based on these criteria, 96.72% of acquisitions were found to be successful when the satellite power was maintained at -141 dBm in the simulator as shown in the histograms in FIGURES 4 and 5.

FIGURE 4. Code-delay estimate within ±1 chip (300 meters).

FIGURE 5. Carrier Doppler estimate within ±150 Hz.

Tracking Sensitivity. Tracking sensitivity for the TUTGNSS receiver was measured to be -151 dBm via both tools, assuming a coherent integration interval of 20 milliseconds. Using AutoPET, 100 tracking iterations were performed at every power level and the average number of satellites tracked was recorded. Using dCAP, this test was performed by selecting one satellite and observing how the receiver C/N₀ tracked this satellite during high and low signal power conditions. Twenty tracking iterations of 90 seconds each were performed for a particular satellite. In each iteration, the satellite power in the simulator was maintained at the nominal condition of -130 dBm (equivalent to 38 dB C/N₀ in the receiver) for the first 30 seconds. Subsequently, the power of the satellite was dropped to -151 dBm (equivalent to 17 dB C/N₀ in the receiver).

As visible in Figure 6, the receiver was able to continue tracking the satellite at -51 dBm in 19 out of the 20 iterations. In the case where tracking was lost, the C/N₀ can be seen to diverge rapidly to 0. To make sure that in the rest of the 19 cases the receiver was really tracking the satellite at low power, the power of the satellite was increased again after an additional 30 seconds. In each of the 19 cases, the receiver successfully continued to track the satellite.

FIGURE 6. Tracking C/N₀ in one tracking channel using dCAP.

3D Position Accuracy and TTFF. Computation of the position fix was performed using a least-squares algorithm without any filtering. Using only AutoPET, 100 position fix iterations were performed and the average 3D error in meters was computed. Within the same test case, the time for achieving a position fix was also recorded. The initial (0–30 seconds) position fix estimates are not very accurate. This is because only the first four acquired satellites are used for the position computation. As more satellites are acquired and tracked, their inclusion into the computation gradually improves the position accuracy to within 1 meter. The average TTFF was computed to 60.59 seconds.

Validity of C/N₀ Estimator. FIGURE 7 presents a comparison of C/N₀ measurements between the TUTGNSS receiver (extracted using dCAP) and a commercial receiver. The input power from the simulator was varied between -130 dBm and -151 dBm with steps of around 2 dB for 10 seconds each. The C/N₀ readings from the two receivers were measured at each power level and plotted on the same scale. The reference power level represents the C/N₀ readings of a hypothetical (ideal) receiver with zero radio front-end losses. As the figure shows, on average there is close conformance between the estimated values of C/N₀ in the two receivers. The difference between the two receivers and the reference is approximately 5 dB, which includes radio front-end noise and other losses. The TUTGNSS receiver displays lower C/N₀ estimation peak-to-peak inconsistency than the commercial receiver.

FIGURE 7. C/N₀ measurement using dCAP: Comparison between TUTGNSS, a commercial, and a hypothetical receiver.

Other Uses of dCAP. During initial prototype validation, we noticed that satellite tracking was inconsistent even under high C/N₀ conditions. dCAP was used to extract detailed baseband tracking information that helped to identify the source of the problem: signal anomalies due to insufficient clock buffering on an experimental RF front end, as shown in FIGURE 8. Such anomalies would have been impossible to detect with traditional black-box testing practices. Once the problem was rectified, dCAP was used once again to monitor the RF front-end signals and performance of the baseband tracking loops, where FIGURES 9 and 10 show a marked improvement in the receiver signal processing and satellite tracking performance.

FIGURE 8. Signal anomaly in the Q-branch signal due to insufficient clock buffering in the experimental RF front end: detected using dCAP.

FIGURE 9. Code Doppler extracted from one tracking loop.

FIGURE 10. Carrier Doppler extracted from one tracking loop using dCAP.

Conclusion

In this article, we have demonstrated the results of the TUTGNSS prototype receiver testing using AutoPET and dCAP. Results were presented, analyzed, and conclusions drawn for the GPS L1 performance of the receiver. Furthermore, the procedures can be easily replicated through software modifications for testing more advanced multi-frequency, multi-constellation modes of the receiver.

Added to the benefits of automation in terms of improved accuracy and personnel efficiency, the proposed AutoPET is a cost-effective solution to anyone working on GNSS receiver technology to understand its most important performance parameters. This tool is portable (software platform-independent), easy to install, and easy to execute on any computer with the basic scientific software. From an academic point of view, dCAP is useful for teaching the spectral characteristics of GNSS signals at every stage from deep inside the receiver to researchers or university students in laboratory exercises. Together, these tools have assisted in the complete characterization of the TUTGNSS receiver at TUT, and can be easily adapted, enhanced, and applied to other research-based receivers as well. In other words, the proposed research has an academic as well as practical appeal.

Acknowledgments

This research work received support from the Tampere Doctoral Programme in Information Science and Engineering (TISE), Nokia Foundation, and the Ulla Tuominen Foundation. It has also been partially supported by the Academy of Finland (under the projects: 251138 “Digitally-Enhanced RF for Cognitive Radio Devices”, and 256175 “Cognitive Approaches for Location in Mobile Environments”). We wish to gratefully acknowledge each of these institutions. This article is based on the paper “Automated Test-bench Infrastructure for GNSS Receivers – Case Study of the TUTGNSS Receiver” presented at the 26th International Technical Meeting of the Satellite Division of The Institute of Navigation held in Nashville, Tennessee, September 16–20, 2013.

Manufacturers

The tests described in this article used a Spirent Federal Systems STR4500 multi-channel GPS/SBAS simulator and a u-blox AG EVK-5P GNSS receiver evaluation kit with a LEA-5P receiver module.

SARANG THOMBRE is a GNSS research scientist in the Department of Navigation and Positioning at the Finnish Geodetic Institute (FGI), Helsinki.

JUSSI RAASAKKA is a GNSS R&D scientist at Honeywell International s.r.o. in the Czech Republic.

TOMMI PAAKKI is a teaching assistant and a doctoral student at the Department of Electronics and Communications Engineering, Tampere University of Technology (TUT).

FRANCESCANTONIO DELLA ROSA is the project manager of the Multitechnology Positioning Professionals (MULTI-POS) Marie Curie Initial Training Network and a research scientist at TUT.

MIKKO VALKAMA is a full professor and the head of the Department of Communications Engineering at TUT.

LAURA RUOTSALAINEN is the deputy head of the Department of Navigation and Positioning and aspecialist research scientist at FGI.

HEIDI KUUSNIEMI is a professor and the acting head of the Department of Navigation and Positioning at FGI.

JARI NURMI is a professor in the Department of Electronics and Communications Engineering at TUT.

FURTHER READING

• Authors’ Conference Paper

“Automated Test-bench Infrastructure for GNSS Receivers – Case Study of the TUTGNSS Receiver” by S. Thombre, J. Raasakka, T. Paakki, F. Della Rosa, M. Valkama, and J. Nurmi in Proceedings of ION GNSS+ 2013, the 26th International Technical Meeting of the Satellite Division of The Institute of Navigation, Nashville, Tennessee, September 16–20, 2013, pp. 1919–1930.

• TUTGNSS

“TUTGNSS – University Based Hardware/Software GNSS Receiver for Research Purposes” by T. Paakki, J. Raasakka, F. Della Rosa, H. Hurskainen, and J. Nurmi, in Proceedings of Ubiquitous Positioning Indoor Navigation and Location Based Service (UPINLBS), 2010, Helsinki, Finland, October 14–15, 2010, doi: 10.1109/UPINLBS.2010.5654337.

• Automated GNSS Receiver Testing

“GPS Interference Testing: Lab, Live, and LightSquared” by P. Boulton, R. Borsato, B. Butler, and K. Judge in InsideGNSS, Vol. 6, No. 4, July/August 2011, pp. 32–45.

“Software-based GNSS Signal Simulators: Past, Present and Possible Future” by S. Thombre, E.S. Lohan, J. Raquet, H. Hurskainen, and J. Nurmi, in Proceedings of ENC GNSS 2010, the European Navigation Conference 2010, Braunschweig, Germany, October 19–21, 2010.

• GNSS Receiver Testing in General

“Simulating GPS Signals: It Doesn’t Have to Be Expensive” by A. Brown, J. Redd, and M.-A. Hutton in GPS World, Vol. 23, No. 5, May 2012, pp. 44–50.

“Realistic Randomization: A New Way to Test GNSS Receivers” by A. Mitelman in GPS World, Vol. 22, No. 3, March 2011, pp. 43–48.

“Record, Replay, Rewind: Testing GNSS Receivers with Record and Playback Techniques” by D.A. Hall in GPS World, Vol. 21, No. 10, October 2010, pp.28–34.

• NMEA 0183

NMEA 0183, The Standard for Interfacing Marine Electronic Devices, Ver. 4.10, published by the National Marine Electronics Association, Severna Park, Maryland, June 2012.

“NMEA 0183: A GPS Receiver Interface Standard” by R.B. Langley in GPS World, Vol. 6, No. 7, July 1995, pp. 54–57.

Unofficial online NMEA 0183 descriptions: “NMEA data”; “NMEA Revealed” by E.S. Raymond, Ver. 2.13, November 2013.
March 1, 2014
Innovation: Ionospheric Modeling Using GPS
Greater Fidelity Using a 3D Approach

By Wei Zhang, Attila Komjathy, Simon Banville, and Richard B. Langley

INNOVATION INSIGHTS by Richard Langley

MAY YOU LIVE IN INTERESTING TIMES. So goes the purported Chinese proverb and curse. When it comes to the ionosphere, an interesting time might indeed be a curse for most users of GPS. The ionosphere – that region of the upper atmosphere where free electrons exist in sufficient numbers to affect the propagation of radio waves – owes its existence primarily to the extreme ultraviolet (EUV) and x-ray photons emitted by the sun. They ionize atoms and molecules in the upper atmosphere, freeing the outer electrons. Mostly the ionosphere is well behaved but it can get quite interesting when it is disturbed by space weather events such as solar flares or coronal mass ejections.

The signals from the GPS satellites are perturbed as they transit the ionosphere. Pseudorange measurements are increased in value (an additional delay) and carrier-phase measurements are decreased (a phase advance). If not fully modeled or otherwise accounted for, the perturbations can decrease the accuracy of GPS positioning, navigation, and timing (PNT). For highest PNT accuracies, observations are made at the two frequencies transmitted by all GPS satellites and because the ionosphere’s effect on radio signals is dispersive, a linear combination of the measurements removes almost all of the ionospheric perturbations. On the other hand, the ionosphere’s effect on single frequency observations must be corrected using a model. Most commonly, the model assumes that all of the electrons in the ionosphere can be compressed into a thin shell at a certain height above the receiver. This permits the computation of an estimate of the vertical ionospheric delay. Then, a mapping function is used to predict the slant delay, the delay contributing to a GPS measurement. The approach works reasonably well, particularly if near-real-time values of vertical delay can be provided to users as is done by the Wide Area Augmentation System and other satellite-based augmentation systems. However, this two-dimensional approach ignores the fact that the electron content of the ionosphere is actually spread out in the vertical direction and so has certain inaccuracies, which can increase when the ionosphere is disturbed.

In an effort to improve ionosphere modeling with potential application to single-frequency GNSS users, a couple of my current graduate students together with a former student, have investigated a three-dimensional approach to ionospheric modeling using empirical orthogonal functions or EOFs to describe the vertical structure of the ionosphere. EOFs reduce the dimensionality of a data set or an empirical model consisting of a large number of interrelated variables, while retaining as much of the variance present in the data set as possible. This is achieved by transforming to a new set of variables, the orthogonal functions, which are uncorrelated (orthogonal), and which are ordered so that the first few retain most of the variation present in all of the original variables. Only three functions are required to account for more than 99 percent of the variability in the International Reference Ionosphere – 2007, for example.

In this month’s column we look at the performance of this 3D approach to modeling the ionosphere including times when the ionosphere is particularly interesting.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas.

Ionospheric modeling plays an important role in improving the accuracies of positioning and navigation, especially for current civil aircraft navigation and mass-market single-frequency users. Measurement-driven models are considered to be among the best candidates for real-time single-frequency positioning owing to their real-time applicability and relatively higher accuracy compared to empirical models, such as the GPS broadcast (also known as Klobuchar) and NeQuick models. A good example of a real-time positioning application is satellite-based augmentation systems (SBAS), including the Wide Area Augmentation System (WAAS), the European Geostationary Navigation Overlay Service (EGNOS), the Japanese MSTAT Satellite-based Augmentation System (MSAS), and the Indian GPS Aided Geo Augmented Navigation system (GAGAN). Because the ionosphere can be the largest error source in single-frequency positioning, the accuracy of ionospheric modeling is critical for single-frequency applications.

Several organizations have been routinely providing ionospheric products to correct errors caused by the ionosphere in the form of ionospheric maps — that is, vertical total electron content (vTEC) at grid points (including regional and global products), such as those from WAAS and the International GNSS Service (IGS), with various processing time delays ranging from near real time to a couple of weeks. Among the earliest works of ionosphere modeling, the University of New Brunswick-Ionospheric Modeling Technique (UNB-IMT) was developed in the mid-1990s. This technique was demonstrated to effectively derive both regional and global total electron content (TEC) maps. However, most of the models, including the current version of UNB-IMT, approximate the ionosphere using a single thin-shell approach with an altitude set at, for example 350 km, which may introduce additional modeling errors up to several TEC units (1 TECU = 10¹⁶ electrons/m²), corresponding to meter-level errors of measurement delay or advance at the GPS L1 frequency.

To overcome any downside of such models, three-dimensional (3D) ionospheric tomographic modeling methods have been proposed and implemented by several groups since the late 1990s. Different from the two-dimensional (2D) single thin-shell ionospheric models, where the parameters to be estimated are associated with TEC, the modeled variables in the tomographic model are related to electron density functions. Therefore, we may expect more complex structures of electron densities (such as those observed during ionospheric storms or in the highly variable equatorial anomaly) to be revealed by the models. A commonly accepted modeling approach is to describe the ionospheric horizontal (longitudinal and latitudinal) variability by a spherical harmonic (SH) expansion up to a specific degree and its vertical dimension modeled by empirical orthogonal functions (EOFs).

However, SH models are not ideal for capturing local variability in the ionosphere because each basis function of spherical harmonics exists over the entire geographic region of interest, such as the entire globe in the case of global modeling. In other words, localized measurements will have influence on the estimated state across the whole globe. As alternative approaches, wavelet, finite element (meshes/pixels), and local-basis-function models have been proposed and implemented to capture the localized information content in the measurements and pass this information on to the end user. On the other hand, the inversion process can occasionally become singular as many of the parameters to be estimated tend to be ineffective and less meaningful. This is especially the case when our goal is to obtain better accuracies with higher order wavelet bases or smaller meshes/pixels. Due to the potential computing and transmitting burden, the two modeling techniques may have more difficulties associated with real-time applications, such as real-time single-frequency positioning, although they have advantages for capturing localized structures in the ionosphere.

Aiming for potential real-time applications of 3D tomographic models, we have extended the UNB-IMT from 2D to 3D by modeling the vertical dimension of the ionosphere using EOFs. In this article, we discuss our approach and report on some initial tests including comparing its performance with the 3D SH approach.

The 2D UNB-IMT was demonstrated to work with various network sizes: regional, baseline-by-baseline, and even single standalone stations. Therefore, it is expected that this technique will help in capturing localized ionospheric structures above small regional networks or above a single standalone station. Additional benefits may be expected for disturbed ionospheric conditions. For assessing the two modeling techniques, a small regional network was chosen to perform station-by-station and batch processes. The performance of both methods with the two processing scenarios has been compared by analyzing the post-fit residuals and vTECs of the state estimation process, as well as the repeatability of estimates of differential code biases (DCBs) for both quiet and disturbed ionospheric conditions.

3D UNB-IMT

Because of the limited number of ionospheric parameters to be estimated, the 2D UNB-IMT was considered suitable for real-time applications, such as real-time single-frequency precise point positioning (PPP) and SBASs. In fact, it can be proven that the modeling method of the current 2D UNB-IMT is identical to the original planar fit of WAAS in nature if the locations of reference stations tend to collocate with WAAS ionospheric grid points (IGPs). Although additional parameters are involved, we believe the 3D UNB-IMT approach with its potential for improved modeling accuracy is still suitable for real-time applications. In this section, we will introduce the 3D UNB-IMT modeling strategy and demonstrate its applicability with a regional network and single standalone stations.

Model Description. In order to clearly present the technique demonstrated in our recent work, we first briefly review the 2D UNB-IMT. Linear polynomial functions were initially proposed for describing the spatial variability of the ionosphere. We model the observed slant TEC (sTEC) between a satellite and a receiver from carrier-phase and pseudorange (code) observations at some epoch as the product of a bilinear polynomial representing the vTEC at the thin-shell ionospheric pierce point (IPP) of the signal raypath and a mapping function that projects the vTEC to sTEC plus receiver and satellite instrumental biases (DCBs). The input variables are the geographic longitude of the IPP referenced to the solar-geomagnetic coordinate system (in other words, the difference between the longitude of the IPP and the longitude of the mean sun) and the difference between the geomagnetic latitude of the IPP and the geomagnetic latitude of the station. We consequently have three polynomial coefficients to estimate for each station: a constant term, one to describe the longitude variations, and one for the latitude variations.

The mapping function used in the model is the standard geometric mapping function, which computes the secant of the zenith angle of the signal geometric ray path at the IPP at a specified shell height. Because of the dependence of the ionosphere on solar radiation and the geomagnetic field, the solar-geomagnetic reference frame is used to compute the TEC over each station in this technique. Since the ionosphere changes more slowly in the sun-fixed reference frame than in the Earth-fixed one, such a reference frame is ideal for producing more accurate TEC estimates.

The initial version of UNB-IMT ignored the non-linear spatial variation of the ionosphere. Non-linear terms are expected to be able to absorb more complex variability of the ionosphere and thus more properly describe the ionosphere in disturbed conditions. Regarding this issue, the drawbacks of some modeling methods based on linear models have been reported: for example, the highly variable ionosphere might be absorbed by the estimated DCBs, making the repeatability of the estimated DCBs (day-to-day variability) correlated with the variability of the ionosphere. To enhance the performance of UNB-IMT, especially under disturbed ionospheric conditions, UNB researchers extended the linear version of UNB-IMT to a quadratic one and assessed it by using a wide-area regional network in North America. This modified approach reduced the post-fit residuals significantly by better modeling the ionospheric variations with the help of the additional second order (non-linear) terms.

To better use a priori information in the development of 3D UNB-IMT, we separate the TEC into a background reference or “known” part and a perturbation or to-be-modeled part. The background reference part of TEC could be calculated from an a priori source of electron density, such as any kind of ionospheric model, including empirical and theoretical ionospheric models. The density, as a function of latitude, longitude, height, and time, is integrated along the raypath between the receiver and a satellite.

Then, the perturbation part of the electron density is modeled by the inner product of EOFs and polynomial functions with associated estimated coefficients to depict the variability of the ionosphere in the vertical and horizontal directions respectively. And this part is similarly integrated along the raypath and added to the reference part along with the DCBs.

Empirical Orthogonal Functions. The EOF method is a method of choice for analyzing the variability of a single field (with only one scalar variable). Variability of the ionosphere with respect to height is needed for the 3D models. The method finds the spatial patterns of variability based on historical data sets (as reflected in empirical or theoretical models). In other words, the modes of variability decomposed by the method are primarily “data modes,” and not necessarily physical or actual real-time models. Due to its noted ability in describing the background ionosphere, the data sets output from the empirical Ionospheric Reference Ionosphere 2007, were utilized to form the EOFs in our technique.

Thus, the data sets of electron densities are realized by uniform sampling at the following variant time scale intervals and specific geographic locations:
- Solar cycle: [1998:1:2008] (year)
- Season of year: [Dec., Mar., Jun., Sep.] (month)
- Time of day: [1:1:24] (hour)
- Day of month: [1:9:28] (day of month)
- Geographic latitude: [30:5:60] (degree)
- Geographic longitude: [280:5:300] (degree),
where the numbers separated by colons correspond to minimum:increment:maximum. The data sets cover the whole spatial area of interest. The data sets of a whole solar cycle in typical equinox and solstice months are used to ensure that the EOFs span the range of profile variations that include the variation in solar EUV and x-ray output. Each electron density profile with respect to height at these locations and time points is sampled in the vertical dimension at [100:2:2000] (km). Figure 1 shows the first three third-order normalized EOFs based on the data sets. The first three eigenvalues account for 92.22, 6.69, and 0.78 percent of the total respectively. Provided the solution is nonsingular, the choice of the highest order of EOFs is a trade off between processing time and modeling accuracy as to the specific network and capability of computer(s). In our current work, the highest order of three was chosen. In this case, the neglected vertical variation of the ionosphere corresponding to higher order EOFs is 0.31 percent.

FIGURE 1. The normalized first three dominant EOFs extracted from the IRI-2007 empirical model.

Once the modeling approach has been constructed, the following task is to estimate the coefficients. Considering the potential real-time applications, a Kalman filter is employed to solve the TEC observation equation. To be specific, the following settings are used. The correlation time is set to five minutes, which correspond to the WAAS update interval for ionospheric grid points. The uncertainty of the dynamic model, 0.008 TECU2/second, is chosen to characterize the potential rapid change of the ionosphere.

Finally, the estimated coefficients provided by the Kalman filter are then used to reconstruct the electron density field.

Testing the Approach

In this section, we report on tests of the 3D UNB-IMT and compare its performance with that of the 3D SH approach. Because of the advantages of sensitivity of 2D UNB-IMT, especially with the single-station processing strategy, it is expected that this technique will help in better capturing localized ionospheric structures above small regional networks or above a single standalone station compared to the 3D SH approach. Additional benefits may be expected for disturbed ionospheric conditions.

For assessing the two modeling techniques, a small regional network of four IGS reference stations located from geographic latitude 39.0° N to 48.1° N and longitude 66.7° W to 77.6° W was chosen to perform single-station and multi-station (network) processing. The stations are GODZ in Greenbelt, Maryland; UNBJ and FRDN in Fredericton, New Brunswick; and VALD in Val d’Or, Quebec. Figure 2 shows the locations of the reference stations chosen for the modeling. The dual-frequency GPS data used for the tests was obtained from October 13–25 (day of year (doy) 286–298) in 2011 with the sampling time interval of 30 seconds. The corresponding values of the interplanetary magnetic field Bz component; the planetary geomagnetic index, Kp; the auroral electrojet index, AE; and the disturbance storm-time index, Dst on these days are shown in FIGURE 3. It is seen that a severe ionospheric storm triggered by a coronal mass ejection from the sun occurred late on October 24 (doy 297), 2011, and continued through the entire day of October 25 (doy 298), 2011. The other days seem relatively quiet. Thus, we chose October 16, 2011, as a typical day with quiet ionospheric conditions and October 25, 2011, as a typical day with disturbed ionospheric conditions in the following tests. The performance of both methods (3D UNB-IMT and SH model) with the two processing scenarios will be compared by analyzing the post-fit residuals and TEC of the state estimation process for both quiet and disturbed ionospheric conditions.

FIGURE 2. The network of the four stations used in the evaluation procedures.

All four reference stations in the small network have the ability to provide both C/A- and P-code pseudorange measurements. In our tests, the P-code observable is used to extract TEC through leveling the corresponding carrier-phase measurements. We used a 15°-elevation-angle cut-off in our study.

Single Station Experiment. The estimated parameters of 2D and 3D UNB-IMT have different physical meanings due to the different modeling strategies. In theory, the 3D UNB-IMT can reproduce the electron densities for any location (horizontal and vertical) at any epoch. Figure 4 shows an example of the electron density profile produced by the linear 3D UNB-IMT in the zenith direction of station FRDN at 12:00 UT on October 16 (doy 289), 2011. Therefore, we will have to integrate electron densities into TEC for the 3D UNB-IMT modeling results if we want to compare how the two approaches have modeled the ionosphere side by side. For the purpose of sensitivity comparison, the results from 2D and 3D UNB-IMT are compared in terms of post-fit residuals as well as time series of estimated vTEC in the single-station processing scenario. As discussed above, we use the GPS data from station FRDN only for October 16 and 25, 2011, in this subsection. The post-fit residuals are calculated as the difference between the measured and estimated biased sTEC.

FIGURE 4. The electron density profile produced by linear 3D UNB-IMT overhead FRDN at 12:00 (UT) on October 16 (doy 289) in 2011.

From the top to bottom panels, Figure 5 shows the estimated vTEC in the zenith direction over the station, post-fit residuals, estimated satellite and receiver DCBs, and unbiased sTEC with respect to local mean solar time series obtained with linear 2D (left-hand panels) and 3D (right-hand panels) UNB-IMT approaches respectively. We use a different color for each satellite to see individual improvement of satellites in terms of post-fit residuals, estimated DCB, and unbiased sTEC. As for the potential improvement of 3D UNB-IMT, we supposed, if the 2D model with single-shell assumption does not depict the variability of the ionosphere quite well (especially the vertical variability of the ionosphere), we should expect to see an improvement from the 3D model in terms of post-fit residuals. As seen in this figure, the 3D UNB-IMT improves the results in terms of post-fit residuals. The means and standard deviations of the residuals with the 2D and 3D UNB-IMT are shown in Table 1.

FIGURE 5. Sensitivity test (the panels from the top to the bottom correspond to: estimated vertical TEC, post-fit residuals, satellite and receiver DCB, slant TEC with respect to local time) between linear 2D (the left-hand panels) and 3D (the right-hand panels) models at FRDN on October 16 (doy 289) in 2011.

TABLE 1. The means and standard deviations of the residuals under the quiet (Q, October 16, 2011) and disturbed or storm (S, October 25, 2011) ionospheric conditions with linear (L) and quadratic (Q) modeling approaches. Units = TECU.

The 3D UNB-IMT with three times as many parameters is allowed to “accommodate” more (vertical) variations of the ionosphere. The benefits are also manifest in the improvement of the estimated vTEC and estimated satellite and receiver DCBs. In terms of estimated vTEC, the smooth variation of TEC may be expected at mid-latitudes during quiet ionospheric conditions without any ionospheric anomaly. The unmodeled variation of TEC in 2D UNB-IMT seen in the post-fit residuals is also manifest as “artificial small jumps” in the vTEC panel. In other words, the 3D UNB-IMT is able to better represent the measurements from low-elevation-angle satellites owing to the EOFs replacing the mapping function. It is the typical case when a satellite comes into or goes out of view of the receiver. The estimated DCBs are relatively constant over the entire day. But it is also found from the estimated DCBs that the results from 2D UNB-IMT have slightly more variability. Both effects seem to be related to the unmodeled errors. The post-fit residuals in the 3D UNB-IMT are closer to the zero mean Gaussian distribution.

Then, we further evaluated the performance of 2D and 3D UNB-IMT under significantly disturbed conditions. Figure 6 shows the results with the same modeling strategies as demonstrated in Figure 5 but on October 25, 2011. Similar conclusions can be drawn from Figure 6, where better results in terms of post-fit residuals are obtained with 3D UNB-IMT (Table 1). In terms of estimated vTEC, the results from both strategies under the disturbed conditions look more irregular than those under the quiet conditions and deviate a little from the sine-wave-like daily variation. Some actual variation of the ionosphere during disturbed conditions may be captured and correctly illustrated as the bumps for both approaches. Furthermore, the unmodeled errors may also be explained as artificial small jumps/bumps in vTEC curves (revealed by the magnitude of post-fit residuals). It is seen that 3D linear UNB-IMT explains more variation of the ionosphere than 2D linear UNB-IMT. However, some residual unmodeled errors may still exist with the 3D model.

FIGURE 6. Sensitivity test (the panels from the top to the bottom correspond to: estimated vertical TEC, residuals, satellite and receiver DCB, slant TEC with respect to local time) between linear 2D (the left-hand panels) and 3D (the right-hand panels) models at FRDN on October 25 (doy 298) in 2011.

As concluded by other investigators, a higher order model could explain more spatial (non-linear) variations of the ionosphere, especially for geomagnetic storm conditions. The results with 2D and 3D quadratic UNB-IMT approaches are shown in Figure 7. In the post-fit residual panels, it can be seen that the residuals with 3D quadratic UNB-IMT are mostly within ±2 TECU except for several small spikes that happened between 0:00 and 4:00 local mean solar time and reflect that not all the electron density variations had been correctly represented by the model used. But it is clear that the 3D quadratic UNB-IMT can significantly improve the modeling precision compared to the 2D quadratic/linear UNB-IMT and 3D linear UNB-IMT. The magnitude of the post-fit residuals shown in this panel is even comparable with the results for the quiet condition shown in Figure 5. In terms of vTEC, a few spurious spikes are occasionally found when processing the data from the four stations with the 3D quadratic model and single-station processing strategy. Other data sources, such as data from incoherent backscatter measurements, may be needed to confirm if the spikes are caused by the instability of the model or actual ionospheric structures. Still, the vTEC curves with 3D quadratic UNB-IMT look smoother than 2D UNB-IMT. In terms of estimated DCBs, it is found that the results with 3D quadratic UNB-IMT approach exhibit relatively fewer perturbations than the other three approaches tested.

FIGURE 7. Sensitivity test (the panels from the top to the bottom correspond to: estimated vertical TEC, residuals, satellite and receiver DCB, slant TEC with respect to local time) between quadratic 2D (the left-hand panels) and 3D (the right-hand panels) models at FRDN on October 25 (doy 298) in 2011.

As we found for the 2D modeling approaches, the single thin-shell assumption with a fixed ionospheric shell height may introduce additional modeling errors. That is mainly because the layer with highest electron density (F2 layer) is not always located at a fixed height. Especially in disturbed ionospheric conditions, such as the case shown in Figures 6 and 7, the layer height would change significantly. Some methods have been proposed and tested with the help of more reliable “true” heights from other resources, such as ionosondes. However, due to the limited number of the instruments deployed and limited information provided (only information from overhead), the applications with these methods would have to be limited to the specific area covered by stations or networks equipped with the instruments. In addition, as to real-time application, the data processing time delay of ionosondes might be another technical issue these methods have to face. Compared with these methods, one benefit of the 3D UNB-IMT is its potential for real-time application for any size of network. Another benefit is its vertical modeling capability to depict vertical variation of electron density so the improved results would also be expected for disturbed ionospheric conditions. It is clearly seen from Figures 6 and 7 that the lowest vTECs around 4:00 LT reach down to 0 TECU with the 2D linear/quatratic UNB-IMT, which are considered as unphysical results. It is confirmed that small biases still exist in the results with the 2D model likely due to the improper shell height chosen (fixed at 350 km for the results shown in this article).

Multi-Station Experiment. When using the modeling scheme for a network solution, we will generally have two possible processing scenarios. One is processing the data of all the stations as a batch, and the other is processing station by station (or baseline by baseline).

The advantages and disadvantages of the batch process can be summarized as follows. It has more redundancies in the Kalman filter to estimate a more stable and reliable set of satellite and receiver DCBs. Due to more measurements as an input (state) of the Kalman filter, the convergence time would be shorter in terms of the estimated DCBs. It would be of benefit for real-time applications if we have limited a priori information about the estimated ionospheric parameters and/or DCBs. However, the batch solution seems to be less sensitive to localized information content than the station-by-station solution. The overall effect of the batch solution is smoothing over the network, reducing the size of some small perturbations. Theoretically, localized measurements should not have significant influence on the estimated state across an extended area or even the entire globe. In other words, the batch solution may be beneficial for relatively small local-area networks, but may not be ideally suited for networks as large as wide-area ones. Another straightforward disadvantage of the batch process is its relatively longer processing time, which might be a downside if it is used for real-time applications.

In the multi-station experiment, we tested the 3D UNB- IMT with a small regional network of four IGS reference stations (Figure 2) to investigate its performance with localized ionospheric variations. We performed tests with two scenarios: batch and station-by-station. Due to space restrictions, we cannot thoroughly report the results we obtained here. Please see the conference paper listed in Further Reading for the full details. Overall, the results we obtained in terms of post-processing residuals were similar to those in the single station experiment. We also found that the 3D UNB-IMT with EOFs seems to be able to better model the measurements with low elevation angles than the 2D UNB-IMT with a mapping function.

Comparing 3D UNB-IMT with SH Model. We have compared the results using the batch processing strategy with those from the SH model. The reason for this approach is that we intended to compare the results of the two processing strategies (UNB-IMT and SH) with identical conditions. That is, both methods processed the data using a batch scheme and estimated both ionospheric parameters and DCBs simultaneously, instead of using some other source or processed results. Therefore, in this case, we can compare the results side by side and evaluate the effectiveness of the estimated ionospheric parameters.

Based on the data from the network of the four stations, the sensitivity of the SH models is lower than that of 3D UNB-IMT, although the number of ionospheric parameters of the SH models is comparable or even larger than that of 3D UNB-IMT. In other words, the ionospheric parameters in 3D UNB-IMT to describe the variability of the ionosphere are more effective and meaningful to such a network scale than those in the 3D SH model.

Given the nature of its basis functions, the SH model is an excellent tool for global modeling, but it has some shortcomings for localized variability modeling. As to larger regional networks with longer baselines, such as those used for WAAS, which covers North America, the difference of the sensitivities between the batch and the station-by-station solutions should be larger than the results we have obtained. However, we cannot conclude that the sensitivity of 3D UNB-IMT is better than that of the 3D SH model with the batch processing strategy for such large regional networks before more tests are conducted. Still, it is clearly seen in our tests that the 3D SH model is not always ideal for regional networks in terms of sensitivity.

We reached similar conclusions for October 25, 2011, where the residuals spread more widely compared with quiet-condition residuals. In the storm conditions, the residuals of the quadratic 3D UNB-IMT spread relatively less than those of other modeling strategies. This is especially the case for the several hours at the beginning of the day, which corresponds to the peak of the Dst and Kp indices shown in Figure 3. The quadratic 3D UNB-IMT seems to have the capacity to handle the ionospheric spatial and temporal variation even during severe storm conditions.

FIGURE 3. Interplanetary magnetic field Bz component, Kp index, AE index, and Dst index during October 13–25 (doy 286–298) in 2011; nT = nanoteslas (Data from World Data Center for Geomagnetism, Kyoto and Goddard Space Flight Center Space Physics Data Facility).

Repeatability of Estimated DCBs. The DCBs not only have influence on the quality (accuracy) of the vTEC estimation, but their repeatability can also provide information to evaluate ionospheric models. This implies that the ionospheric models that have the capability to estimate/eliminate more accurate DCBs, independent of ionospheric variability, are preferable. We carried out a number of tests to evaluate the repeatability of estimated DCB values using the 2D and 3D UNB-IMT approaches as well as the 3D SH technique under both quiet and disturbed ionospheric conditions. For quiet ionospheric conditions, the performance of all the tested models looks comparable, although the quadratic 3D UNB-IMT performs slightly better than the others. As to the disturbed conditions, the quadratic 2D/3D UNB-IMT seems be able to provide more stable DCBs than the other models. However, the improvement of the extension from 2D to 3D is slight for the quadratic models, although it is significant for the linear models. The performance of the 3D SH model looks fairly poor compared to 3D UNB-IMT for regional modeling. Consult the conference paper for further details.

Conclusions and Future Research

In the work described in this article, we extended the UNB-IMT from 2D to 3D and compared the performance between them in station-by-station and batch processing scenarios for both quiet and storm ionospheric conditions. We used the data from a small regional network of dual-frequency GPS receivers. The DCBs and ionospheric delays were estimated at the same time by a Kalman filter. The newly developed approach was evaluated by analyzing the post-fit residuals, TEC of the state estimation process, and the repeatability of estimates of DCBs.

In the single-station processing, the improvement of 3D UNB-IMT has been demonstrated in both quiet and disturbed ionospheric conditions in terms of post-fit residuals. The 3D UNB-IMT with more parameters allows the depiction of more complex (vertical) variability of the ionosphere. The 3D UNB-IMT is able to better deal with the measurements from low-elevation-angle satellites owing to EOFs replacing the mapping function. The artificial jumps with 2D UNB-IMT when satellites come into or go out of view of the receiver have been properly handled by the 3D UNB-IMT. In addition, the time series of estimated DCBs with 3D UNB-IMT exhibit less perturbation than the results with 2D UNB-IMT.

As to the multi-station (network) processing, it is confirmed that the station-by-station solution is more sensitive to localized information than the batch solution. Based on the results from our research, station-by-station processing with 3D UNB-IMT is suggested to increase chances to catch localized ionospheric structures.

The repeatability of estimated DCBs was investigated as another indicator to evaluate the viability ofionospheric models.

Before the 3D UNB-IMT is tested in the positioning domain for single-frequency positioning, it is worth validating the model with other data sources. In addition, the potential benefits of 3D UNB-IMT during extremely disturbed ionospheric conditions is worth investigating further.

Acknowledgments

We would like to thank the IGS and the Crustal Dynamics Data Information System for providing the GPS data, and we acknowledge the financial contribution of the Natural Sciences and Engineering Research Council of Canada for supporting the first and last authors. This article is based on the paper “Eliminating Potential Errors Caused by the Thin Shell Assumption: An Extended 3D UNB Ionospheric Modelling Technique” presented at the 26th International Technical Meeting of the Satellite Division of The Institute of Navigation, Nashville, Tennessee, September 16–20, 2013.

WEI ZHANG received his M.Sc. degree in space science in 2009 from the School of Earth and Space Science of Peking University, China. He is currently an M.Sc.E. student in the Department of Geodesy and Geomatics Engineering at University of New Brunswick (UNB) under the supervision of Dr. Richard B. Langley.

ATTILA KOMJATHY is a principal investigator at the California Institute of Technology Jet Propulsion Laboratory and an adjunct professor at UNB, specializing in remote sensing techniques using GPS. He received his Ph.D. from the Department of Geodesy and Geomatics Engineering of UNB in 1997.

SIMON BANVILLE works for the Geodetic Survey Division of Natural Resources Canada on real-time precise point positioning (PPP) using global navigation satellite systems. He is also in the process of completing his Ph.D. degree at UNB under the supervision of Dr. Langley.

FURTHER READING

• Authors’ Conference Paper

“Eliminating Potential Errors Caused by the Thin Shell Approximation: An Extended 3D UNB Ionospheric Modelling Technique” by W. Zhang, R.B. Langley, A. Komjathy, and S. Banville in Proceedings of ION GNSS+ 2013, the 26th International Technical Meeting of the Satellite Division of The Institute of Navigation, Nashville, Tennessee, September 16–20, 2013, pp. 2447–2462.

• 2D Ionosphere Modeling

“SBAS Ionospheric Modeling with the Quadratic Approach: Reducing the Risks” by H. Rho, R. Langley, and A. Komjathy in Proceedings of ION GNSS 2005, the 18th International Technical Meeting of the Satellite Division of The Institute of Navigation, Long Beach, California, September 13–16, 2005, pp. 723–734.

Global Ionospheric Total Electron Content Mapping Using the Global Positioning System by A. Komjathy, Ph.D. dissertation, Technical Report No. 188, Department of Geodesy and Geomatics Engineering, University of New Brunswick, Fredericton, New Brunswick, Canada, 1997.

“Improvement of a Global Ionospheric Model to Provide Ionospheric Range Error Corrections for Single-frequency GPS Users” by A. Komjathy and R. Langley in Proceedings of the 52nd Annual Meeting of The Institute of Navigation, Cambridge, Massachusetts, January 22–24, 1996, pp. 557–566.

• 3D (4D) Ionosphere Modeling

“Comparison of 4D Tomographic Mapping Versus Thin-shell Approximation for Ionospheric Delay Corrections for Single-frequency GPS Receivers over North America” by D.J. Allain and C.N. Mitchell in GPS Solutions, Vol. 14, No. 3, 2009, pp. 279–291, doi: 10.1007/s10291-009-0153-0.

“Regional 4-D modeling of the Ionospheric Electron Density” by M. Schmidt, D. Bilitza, C. Shum, and C. Zeilhofer in Advances in Space Research, Vol. 42, No. 4, 2008, pp. 782–790, doi: 10.1016/j.asr.2007.02.050.

“History, Current State, and Future Directions of Ionospheric Imaging” by G.S. Bust and C.N. Mitchell in Reviews of Geophysics, Vol. 46, No. 1, RG1003, March 2008, doi: 10.1029/2006RG000212.

“Development of the Global Assimilative Ionospheric Model” by C. Wang, G. Hajj, X. Pi, I.G. Rosen, and B. Wilson in Radio Science, Vol. 39, No. 1, RS1S06, February 2004, doi: 10.1029/2002RS002854.

Contributions to the 3D Ionospheric Sounding with GPS Data by M. García-Fernández, Ph.D. dissertation, Research Group of Astronomy and Geomatics, Universitat Politècnica de Catalunya, Barcelona, Spain, 2004. Available online in three parts:

http://www.tesisenred.net/bitstream/handle/10803/7015/01Mgf01de03.pdf?sequence=1

http://www.tesisenred.net/bitstream/handle/10803/7015/01Mgf01de03.pdf?sequence=2

http://www.tesisenred.net/bitstream/handle/10803/7015/01Mgf01de03.pdf?sequence=3.

• Ionospheric Reference Models

“The NeQuick Model Genesis, Uses and Evolution” by S.M. Radicella in Annals of Geophysics, Vol. 52, No. 3/4, June/August 2009, pp. 417–422, doi: 10.4401/ag-4597.

“International Reference Ionosphere 2007: Improvements and New Parameters” by D. Bilitza and B. Reinisch in Advances in Space Research, Vol. 42, No. 4, 2008, pp. 599–609, doi: 10.1016/j.asr.2007.07.048.

• Space Weather and the Ionosphere

“GNSS and the Ionosphere: What’s in Store for the Next Solar Maximum” by A.B.O. Jensen and C. Mitchell in GPS World, Vol. 22, No. 2, February 2011, pp. 40–48.

“Space Weather: Monitoring the Ionosphere with GPS” by A. Coster, J. Foster, and P. Erickson in GPS World, Vol. 14, No. 5, May 2003, pp. 42–49.

“GPS, the Ionosphere, and the Solar Maximum” by R.B. Langley in GPS World, Vol. 11, No. 7, July 2000, pp. 44–49.

• Empirical Orthogonal Functions

“Empirical Orthogonal Functions and Related Techniques in Atmospheric Science: A Review” by A. Hannachi, I.T. Jolliffe, and D.B. Stephenson in International Journal of Climatology, Vol. 27, No. 9, July 2007, pp. 1119–1152, doi: 10.1002/joc.1499.

“Empirical Orthogonal Functions: The Medium is the Message” by A.H. Monahan, J.C. Fyfe, M.H.P. Ambaum, D.B. Stephenson, and G.R. North in Journal of Climate, Vol. 22, No. 24, December 2009, pp. 6501–6514, doi: 10.1175/2009JCLI3062.1.

A Manual for EOF and SVD Analyses of Climatic Data by H. Bjornsson and S. Venegas, Report No. 97-1, Department of Atmospheric and Oceanic Sciences and Centre for Climate and Global Change Research, McGill University, Montreal, February 1997.
February 1, 2014
Innovation: Cycle Slips

Detection and Correction Using Inertial Aiding

By Malek O. Karaim, Tashfeen B. Karamat, Aboelmagd Noureldin, Mohamed Tamazin, and Mohamed M. Atia

A team of university researchers has developed a technique combining GPS receivers with an inexpensive inertial measuring unit to detect and repair cycle slips with the potential to operate in real time.

INNOVATION INSIGHTS by Richard Langley

DRUM ROLL, PLEASE. The “Innovation” column and GPS World are celebrating a birthday. With this issue, we have started the 25th year of publication of the magazine and the column, which appeared in the very first issue and has been a regular feature ever since. Over the years, we have seen many developments in GPS positioning, navigation, and timing with a fair number documented in the pages of this column.

In January 1990, GPS and GLONASS receivers were still in their infancy. Or perhaps their toddler years. But significant advances in receiver design had already been made since the introduction around 1980 of the first commercially available GPS receiver, the STI-5010, built by Stanford Telecommunications, Inc. It was a dual-frequency, C/A- and P-code, slow-sequencing receiver. Cycling through four satellites took about five minutes, and the receiver unit alone required about 30 centimeters of rack space. By 1990, a number of manufacturers were offering single or dual frequency receivers for positioning, navigation, and timing applications. Already, the first handheld receiver was on the market, the Magellan NAV 1000. Its single sequencing channel could track four satellites. Receiver development has advanced significantly over the intervening 25 years with high-grade multiple frequency, multiple signal, multiple constellation GNSS receivers available from a number of manufacturers, which can record or stream measurements at data rates up to 100 Hz. Consumer-grade receivers have proliferated thanks, in part, to miniaturization of receiver chips and modules. With virtually every cell phone now equipped with GPS, there are over a billion GPS users worldwide. And the chips keep getting smaller. Complete receivers on a chip with an area of less than one centimeter squared are common place. Will the “GPS dot” be in our near future?

The algorithms and methods used to obtain GPS-based positions have evolved over the years, too. By 1990, we already had double-difference carrier-phase processing for precise positioning. But the technique was typically applied in post-processing of collected data. It is still often done that way today. But now, we also have the real-time kinematic (or RTK) technique to achieve similar positioning accuracies in real time and the non-differenced precise point positioning technique, which does not need base stations and which is also being developed for real-time operation. But in all this time, we have always had a “fly in the ointment” when using carrier-phase observations: cycle slips. These are discontinuities in the time series of carrier-phase measurements due to the receiver temporarily losing lock on the carrier of a GPS signal caused by signal blockage, for example. Unless cycle slips are repaired or otherwise dealt with, reduction in positioning accuracy ensues. Scientists and engineers have developed several ways of handling cycle slips not all of which are capable of working in real time. But now, a team of university researchers has developed a technique combining GPS receivers with an inexpensive inertial measuring unit to detect and repair cycle slips with the potential to operate in real time. They describe their system in this month’s column.

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas.

GPS carrier-phase measurements can be used to achieve very precise positioning solutions. Carrier-phase measurements are much more precise than pseudorange measurements, but they are ambiguous by an integer number of cycles. When these ambiguities are resolved, sub-centimeter levels of positioning can be achieved.

However, in real-time kinematic applications, GPS signals could be lost temporarily because of various disturbing factors such as blockage by trees, buildings, and bridges and by vehicle dynamics. Such signal loss causes a discontinuity of the integer number of cycles in the measured carrier phase, known as a cycle slip. Consequently, the integer counter is reinitialized, meaning that the integer ambiguities become unknown again. In this event, ambiguities need to be resolved once more to resume the precise positioning and navigation process. This is a computation-intensive and time-consuming task. Typically, it takes at least a few minutes to resolve the ambiguities.

The ambiguity resolution is even more challenging in real-time navigation due to receiver dynamics and the time-sensitive nature of the required kinematic solution. Therefore, it would save effort and time if we could detect and estimate the size of these cycle slips and correct the measurements accordingly instead of resorting to a new ambiguity resolution. In this article, we will briefly review the cause of cycle slips and present a procedure for detecting and correcting cycle slips using a tightly coupled GPS/inertial system, which could be used in real time. We will also discuss practical tests of the procedure.

Cycle Slips and Their Management

A cycle slip causes a jump in carrier-phase measurements when the receiver phase tracking loops experience a temporary loss of lock due to signal blockage or some other disturbing factor. On the other hand, pseudoranges remain unaffected. This is graphically depicted in FIGURE 1. When a cycle slip happens, the Doppler (cycle) counter in the receiver restarts, causing a jump in the instantaneous accumulated phase by an integer number of cycles. Thus, the integer counter is reinitialized, meaning that ambiguities are unknown again, producing a sudden change in the carrier-phase observations.

FIGURE 1. A cycle slip affecting phase measurements but not the pseudoranges.

Once a cycle slip is detected, it can be handled in two ways. One way is to repair the slip. The other way is to reinitialize the unknown ambiguity parameter in the phase measurements. The former technique requires an exact estimation of the size of the slip but could be done instantaneously. The latter solution is more secure, but it is time-consuming and computationally intensive. In our work, we follow the first approach, providing a real-time cycle-slip detection and correction algorithm based on a GPS/inertial integration scheme.

GPS/INS Integration

An inertial navigation system (INS) can provide a smoother and more continuous navigation solution at higher data rates than a GPS-only system, since it is autonomous and immune to the kinds of interference that can deteriorate GPS positioning quality. However, INS errors grow with time due to the inherent mathematical double integration in the mechanization process. Thus, both GPS and INS systems exhibit mutually complementary characteristics, and their integration provides a more accurate and robust navigation solution than either stand-alone system. GPS/INS integration is often implemented using a filtering technique. A Kalman filter is typically selected for its estimation optimality and time-recursion properties.

The two major approaches of GPS/INS integration are loosely coupled and tightly coupled. The former strategy is simpler and easier to implement because the inertial and GPS navigation solutions are generated independently before being weighted together by the Kalman filter. There are two main drawbacks with this approach: 1) signals from at least four satellites are needed for a navigation solution, which cannot always be guaranteed; and 2) the outputs of the GPS Kalman filter are time correlated, which has a negative impact upon the system performance. The latter strategy performs the INS/GPS integration in a single centralized Kalman filter. This architecture eliminates the problem of correlated measurements, which arises due to the cascaded Kalman filtering in the loosely coupled approach. Moreover, the restriction of visibility of at least four satellites is removed. We specifically use a tightly coupled GPS/reduced inertial sensor system approach.

Reduced Inertial Sensor System. Recently, microelectromechanical system or MEMS-grade inertial sensors have been introduced for low-cost navigation applications. However, these inexpensive sensors have complex error characteristics.

Therefore, current research is directed towards the utilization of fewer numbers of inertial sensors inside the inertial measurement unit (IMU) to obtain the navigation solution.

The advantage of this trend is twofold. The first is avoidance of the effect of inertial sensor errors. The second is reduction of the cost of the IMU in general. One such minimization approach, and the one used in our work, is known as the reduced inertial sensor system (RISS). The RISS configuration uses one gyroscope, two accelerometers, and a vehicle wheel-rotation sensor. The gyroscope is used to observe the changes in the vehicle’s orientation in the horizontal plane. The two accelerometers are used to obtain the pitch and roll angles. The wheel-rotation sensor readings provide the vehicle’s speed in the forward direction. FIGURE 2 shows a general view of the RISS configuration.

FIGURE 2. A general view of the RISS configuration.

A block diagram of the tightly coupled GPS/RISS used in our work is shown in FIGURE 3. At this stage, the system uses GPS pseudoranges together with the RISS observables to compute an integrated navigation solution. In this three-dimensional (3D) version of RISS, the system has a total of nine states. These states are the latitude, longitude, and altitude errors ( ; the east, north, and up velocity errors ; the azimuth error ; the error associated with odometer-driven acceleration ; and the gyroscope error .

The nine-state error vector x_k at time t_k is expressed as:
    (1)

FIGURE 3. Tightly coupled integration of GPS/RISS using differential pseudorange measurements.

Cycle Slip Detection and Correction

Cycle slip handling usually happens in two discrete steps: detection and fixing or correction. In the first step, using some testing quantity, the location (or time) of the slip is found. During the second step, the size of the slip is determined, which is needed along with its location to fix the cycle slip. Various techniques have been introduced by researchers to address the problem of cycle-slip detection and correction. Different measurements and their combinations are used including carrier phase minus code (using L1 or L2 measurements), carrier phase on L1 minus carrier phase on L2, Doppler (on L1 or L2), and time-differenced phases (using L1 or L2). In GPS/INS integration systems, the INS is used to predict the required variable to test for a cycle slip, which is usually the true receiver-to-satellite range in double-difference (DD) mode, differencing measurements between a reference receiver and the roving receiver and between satellites. In this article, we introduce a tightly coupled GPS/RISS approach for cycle-slip detection and correction, principally for land vehicle navigation using a relative-positioning technique.

Principle of the Algorithm. The proposed algorithm compares DD L1 carrier-phase measurements with estimated values derived from the output of the GPS/RISS system. In the case of a cycle slip, the measurements are corrected with the calculated difference. A general overview of the system is given in FIGURE 4.

FIGURE 4. The general flow diagram of the proposed algorithm.

The number of slipped cycles is given by
   (2)
where
is the DD carrier-phase measurement (in cycles)
is DD estimated carrier phase value (in cycles).
is compared to a pre-defined threshold μ . If the threshold is exceeded, it indicates that there is a cycle slip in the DD carrier-phase measurements.

Theoretically, would be an integer but because of the errors in the measured carrier phase as well as errors in the estimations coming from the INS system, will be a real or floating-point number.

The estimated carrier-phase term in Equation (2) is obtained as follows:
    (3)
where
λ is the wavelength of the signal carrier (in meters)
are the estimated ranges from the rover to satellites i and j respectively (in meters)
are known ranges from the base to satellites i and j respectively (in meters).
What we need to get from the integrated GPS/RISS system is the estimated range vector from the receiver to each available satellite ( ). Knowing our best position estimate, we can calculate ranges from the receiver to all available satellites through:
(4)
where
is the calculated range from the receiver to the mth satellite
x^KF is the receiver position obtained from GPS/RISS Kalman filter solution
x^m is the position of the mth satellite
M is the number of available satellites.
Then, the estimated DD carrier-phase term in Equation (3) can be calculated and the following test quantity in Equation (2) can be applied:
   (5)
If a cycle slip occurred in the ith DD carrier-phase set, the corresponding set is instantly corrected for that slip by:
   (6)
where s is the DD carrier-phase-set number in which the cycle slip has occurred.

Experimental Work

The performance of the proposed algorithm was examined on the data collected from several real land-vehicle trajectories. A high-end tactical grade IMU was integrated with a survey-grade GPS receiver to provide the reference solution. This IMU uses three ring-laser gyroscopes and three accelerometers mounted orthogonally to measure angular rate and linear acceleration. The GPS receiver and the IMU were integrated in a commercial package. For the GPS/RISS solution, the same GPS receiver and a MEMS-grade IMU were used. This IMU is a six-degree of freedom inertial system, but data from only the vertical gyroscope, the forward accelerometer, and the transversal accelerometer was used. TABLE 1 gives the main characteristics of both IMUs. The odometer data was collected using a commercial data logger through an On-Board Diagnostics version II (OBD-II) interface. Another GPS receiver of the same type was used for the base station measurements. The GPS data was logged at 1 Hz.

Table 1. Characteristics of the MEMS and tactical grade IMUs.

Several road trajectories were driven using the above-described configuration. We have selected one of the trajectories, which covers several real-life scenarios encountered in a typical road journey, to show the performance of the proposed algorithm. The test was carried out in the city of Kingston, Ontario, Canada. The starting and end point of the trajectory was near a well-surveyed point at Fort Henry National Historic Site where the base station receiver was located. The length of the trajectory was about 30 minutes, and the total distance traveled was about 33 kilometers with a maximum baseline length of about 15 kilometers. The trajectory incorporated a portion of Highway 401 with a maximum speed limit of 100 kilometers per hour and suburban areas with a maximum speed limit of 80 kilometers per hour. It also included different scenarios including sharp turns, high speeds, and slopes.

FIGURE 5 shows measured carrier phases at the rover for the different satellites. Some satellites show very poor presence whereas some others are consistently available. Satellites elevation angles can be seen in FIGURE 6.

FIGURE 5. Measured carrier phase at the rover.

FIGURE 6. Satellite elevation angles.

Results

We start by showing some results of carrier-phase estimation errors. Processing is done on what is considered to be a cycle-slip-free portion of the data set for some persistent satellites (usually with moderate to high elevation angles). Then we show results for the cycle-slip-detection process by artificially introducing cycle slips in different scenarios. In the ensuing discussion (including tables and figures), we show results indicating satellite numbers without any mention of reference satellites, which should be implicit as we are dealing with DD data.

FIGURE 7 shows DD carrier-phase estimation errors whereas FIGURE 8 shows DD measured carrier phases versus DD estimated carrier phases for sample satellite PRN 22.

FIGURE 7. DD-carrier-phase estimation error, reference satellite with PRN 22.

FIGURE 8. Measured versus estimated DD carrier phase, reference satellite with PRN 22.

As can be seen in TABLE 2, the root-mean-square (RMS) error varies from 0.93 to 3.58 cycles with standard deviations from 0.85 to 2.47 cycles. Estimated phases are approximately identical to the measured ones. Nevertheless, most of the DD carrier-phase estimates have bias and general drift trends, which need some elaboration. In fact, the bias error can be the result of more than one cause. The low-cost inertial sensors always have bias in their characteristics, which plays a major role in this. The drift is further affecting relatively lower elevation angle satellites which can also be attributed to more than one reason. Indeed, one reason for choosing this specific trajectory, which was conducted in 2011, was to test the algorithm with severe ionospheric conditions as the year 2011 was close to a solar maximum: a period of peak solar activity in the approximately 11-year sunspot cycle.

Table 2. Estimation error for DD carrier phases (in cycles).

Moreover, the time of the test was in the afternoon, which has the maximum ionospheric effects during the day. Thus, most part of the drift trend must be coming from ionospheric effects as the rover is moving away from the base receiver during this portion of the trajectory. Furthermore, satellite geometry could contribute to this error component. Most of the sudden jumps coincide with, or follow, sharp vehicle turns and rapid tilts. Table 2 shows the averaged RMS and standard deviation (std) DD carrier-phase estimation error for the sample satellite-pairs. We introduced cycle slips at different rates or intensities and different sizes to simulate real-life scenarios. Fortunately, cycle slips are usually big as mentioned earlier and this was corroborated by our observations from real trajectory data. Therefore, it is more important to detect and correct for bigger slips in general.

Introducing and Detecting Cycle Slips. To test the robustness of the algorithm, we started with an adequate cycle slip size. Cycle slips of size 10–1000 cycles were introduced with different intensities. These intensities are categorized as few (1 slip per 100 epochs), moderate (10 slips per 100 epochs), and severe (100 slips per 100 epochs). This was applied for all DD carrier-phase measurement sets simultaneously. The threshold was set to 1.9267 (average of RMS error for all satellite-pairs) cycles. Four metrics were used to describe the results. Mean square error (MSE); accuracy, the detected cycle slip size with respect to the introduced size; True detection (TD) ratio; and Mis-detection (MD) ratio. Due to space constraints and the similarity between results for different satellites, we only show results for the reference satellite with PRN 22. FIGURES 9–12 show introduced versus calculated cycle slips along with the corresponding detection error for sample satellites in the different scenarios. TABLES 3–5 summarize these results.

FIGURE 9. Introduced and calculated cycle slips (upper plot) and detection error (lower plot). Few cycle slips case, reference satellite with PRN 22.

FIGURE 10. Introduced and calculated cycle slips (upper plot) and detection error (lower plot). Moderate cycle slips case, reference satellite with PRN 22.

FIGURE 11. Introduced and calculated cycle slips (upper plot) and detection error (lower plot). Intensive cycle slips case, reference satellite with PRN 22.

FIGURE 12. Introduced and calculated cycle slips (upper plot) and detection error (lower plot). Small cycle slips case, reference satellite with PRN 22.

Table 3. Few slips (1 slip per 100 epochs).

Table 4. Moderate slips (10 slips per 100 epochs).

Table 5. Intensive slips (100 slips per 100 epochs).

All introduced cycle slips were successfully detected in all of the few, moderate, and severe cases with very high accuracy. A slight change in the accuracy (increasing with higher intensity) among the different scenarios shows that detection accuracy is not affected by cycle-slip intensity. Higher mis-detection ratios for smaller cycle-slip intensity comes from bigger error margins than the threshold for several satellite pairs. However, this is not affecting the overall accuracy strongly as all mis-detected slips are of comparably very small sizes. MD ratio is zero in the intensive cycle-slip case as all epochs contain slips is an indicator of performance compromise with slip intensity.

It is less likely to have very small cycle slips (such as 1 to 2 cycles) in the data and usually it will be hidden with the higher noise levels in kinematic navigation with low-cost equipment. However, we wanted to show the accuracy of detection in this case. We chose the moderate cycle slip intensity for this test. TABLE 6 summarizes results for all satellites.

Table 6. Small slips (1–2 cycles) at moderate intensity (10 slips per 100 epochs).

We get a moderate detection ratio and modest accuracy as the slips are of sizes close to the threshold. The MSE values are not far away from the case of big cycle slips but with higher mis-detection ratio.

Conclusions

The performance of the proposed algorithm was examined on several real-life land vehicle trajectories, which included various driving scenarios including high and slow speeds, sudden accelerations, sharp turns and steep slopes. The road testing was designed to demonstrate the effectiveness of the proposed algorithm in different scenarios such as intensive and variable-sized cycle slips.

Results of testing the proposed method showed competitive detection rates and accuracies comparable to existing algorithms that use full MEMS IMUs. Thus with a lower cost GPS/RISS integrated system, we were able to obtain a reliable phase-measurement-based navigation solution. Although the testing discussed in this article involved post-processing of the actual collected data at the reference station and the rover, the procedure has been designed to work in real time where the measurements made at the reference station are transmitted to the rover via a radio link. This research has a direct influence on navigation in real-time applications where frequent cycle slips occur and resolving integer ambiguities is not affordable because of time and computational reasons and where system cost is an important factor.

Acknowledgments

This article is based on the paper “Real-time Cycle-slip Detection and Correction for Land Vehicle Navigation using Inertial Aiding” presented at ION GNSS+ 2013, the 26th International Technical Meeting of the Satellite Division of The Institute of Navigation held in Nashville, Tennessee, September 16–20, 2013.

Manufacturers

The research reported in this article used a Honeywell Aerospace HG1700 AG11 tactical-grade IMU and a NovAtel OEM4 GPS receiver integrated in a NovAtel G2 Pro-Pack SPAN unit, a Crossbow Technology (now Moog Crossbow) IMU300CC MEMS-grade IMU, an additional NovAtel OEM4 receiver at the base station, a pair of NovAtel GPS-702L antennas, and a Davis Instruments CarChip E/X 8225 OBD-II data logger.

Malek Karaim is a Ph.D. student in the Department of Electrical and Computer Engineering of Queen’s University, Kingston, Ontario, Canada.

Tashfeen Karamat is a doctoral candidate in the Department of Electrical and Computer Engineering at Queen’s University.

Aboelmagd Noureldin is a cross-appointment professor in the Departments of Electrical and Computer Engineering at both Queen’s University and the Royal Military College (RMC) of Canada, also in Kingston.

Mohamed Tamazin is a Ph.D. student in the Department of Electrical and Computer Engineering at Queen’s University and a member of the Queen’s/RMC NavINST Laboratory.

Mohamed M. Atia is a research associate and deputy director of the Queen’s/RMC NavINST Laboratory.

FURTHER READING

• Cycle Slips

“Instantaneous Cycle-Slip Correction for Real-Time PPP Applications” by S. Banville and R.B. Langley in Navigation, Vol. 57, No. 4, Winter 2010–2011, pp. 325–334.

“GPS Cycle Slip Detection and Correction Based on High Order Difference and Lagrange Interpolation” by H. Hu and L. Fang in Proceedings of PEITS 2009, the 2nd International Conference on Power Electronics and Intelligent Transportation System, Shenzhen, China, December 19–20, 2009, Vol. 1, pp. 384–387, doi: 10.1109/PEITS.2009.5406991.

“Cycle Slip Detection and Fixing by MEMS-IMU/GPS Integration for Mobile Environment RTK-GPS” by T. Takasu and A. Yasuda in Proceedings of ION GNSS 2008, the 21st International Technical Meeting of the Satellite Division of The Institute of Navigation, Savannah, Georgia, September 16–19, 2008, pp. 64–71.

“Instantaneous Real-time Cycle-slip Correction of Dual-frequency GPS Data” by D. Kim and R. Langley in Proceedings of KIS 2001, the International Symposium on Kinematic Systems in Geodesy, Geomatics and Navigation, Banff, Alberta, June 5–8, 2001, pp. 255–264.

“Carrier-Phase Cycle Slips: A New Approach to an Old Problem” by S.B. Bisnath, D. Kim, and R.B. Langley in GPS World, Vol. 12, No. 5, May 2001, pp. 46-51.

“Cycle-Slip Detection and Repair in Integrated Navigation Systems” by A. Lipp and X. Gu in Proceedings of PLANS 1994, the IEEE Position Location and Navigation Symposium, Las Vegas, Nevada, April 11–15, 1994, pp. 681–688, doi: 10.1109/PLANS.1994.303377.

Short-Arc Orbit Improvement for GPS Satellites by D. Parrot, M.Sc.E. thesis, Department of Geodesy and Geomatics Engineering Technical Report No. 143, University of New Brunswick, Fredericton, New Brunswick, Canada, June 1989.

• Reduced Inertial Sensor Systems

“A Tightly-Coupled Reduced Multi-Sensor System for Urban Navigation” by T. Karamat, J. Georgy, U. Iqbal, and N. Aboelmagd in Proceedings of ION GNSS 2009, the 22nd International Technical Meeting of the Satellite Division of The Institute of Navigation, Savannah, Georgia, September 22–25, 2009, pp. 582–592.

“An Integrated Reduced Inertial Sensor System – RISS / GPS for Land Vehicle” by U. Iqbal, A. Okou, and N. Aboelmagd in Proceedings of PLANS 2008, the IEEE/ION Position Location and Navigation Symposium, Monterey, California, May 5–8, 2008, pp. 1014–1021, doi: 10.1109/PLANS.2008.4570075.

• Integrating GPS and Inertial Systems

Fundamentals of Inertial Navigation, Satellite-based Positioning and their Integration by N. Aboelmagd, T. B. Karmat, and J. Georgy. Published by Springer-Verlag, New York, New York, 2013.

Aided Navigation: GPS with High Rate Sensors by J. A. Farrell. Published by McGraw-Hill, New York, New York, 2008.

Global Positioning Systems, Inertial Navigation, and Integration, 2nd edition, by M.S. Grewal, L.R. Weill, and A.P. Andrews. Published by John Wiley & Sons, Inc., Hoboken, New Jersey, 2007.

January 1, 2014