Equations 3, 4 & 5Working Papers explore the technical and scientific themes that underpin GNSS programs and applications. This regular column is coordinated by Prof. Dr.-Ing. Günter Hein, head of Europe’s Galileo Operations and Evolution.
Working Papers explore the technical and scientific themes that underpin GNSS programs and applications. This regular column is coordinated by Prof. Dr.-Ing. Günter Hein, head of Europe’s Galileo Operations and Evolution.
GNSS reflectometry (GNSS-R) is nowadays a well-established method to remotely sense various properties of the Earth’s surface via reflected GNSS signals. Ground-based, airborne, and satellite-based GNSS-R systems have been used in numerous applications such as estimating wind speed (and eventually wind direction), soil moisture, ice properties, or altimetry data.
These applications have employed various types of GNSS receivers, modified commercial receivers, FPGA-based receivers, and software receivers, which particularly benefit post-processing techniques. To date, most of these have employed the GPS C/A signal for obvious reasons. However, the next generation of GNSS-R systems may exploit signals from multiple GNSSes and on multiple frequencies. Furthermore, software receiver technology and hard disk storage space is now mature enough to process even wideband signals in real-time on a conventional PC.
GNSS-R Demonstration Application
In order to demonstrate the utility of using other GNSS signals for reflectometry techniques, we undertook a demonstration of an altimetry application that processed GPS L1 and Galileo E1/E5a/E5b signals.
We used a multi-frequency and multi-sensor software receiver with appropriate GPS and Galileo dual RF frontends that operated in real-time. The receiver supports GNSS-R by means of the following elements:
- Dual-antenna input. Ideally two RF frontends can be synchronized to allow reception of up to four RF frequency bands from each of the two antennas. (See the accompanying photo of our equipment setup for the demonstration test in Austria.)
- Acquisition and tracking. The receiver supports all GNSS and millions of correlators to acquire signals quickly. Important tracking parameters such as integration time, correlator positions and loop parameters need to be freely configurable.
- Channel slaving (open loop tracking). It must be possible to link two tracking channels such that the slave (reflected signal) channel uses the numerically controlled oscillator (NCO) values from a master (line-of-sight) channel. Offsets need to be applied and slaving can be used for carrier-only or for carrier and code processing.
- Multi-correlator (Doppler/delay map). The receiver computes the correlation function on a configurable grid of Doppler and code delay values. A coherent integration time of several seconds is possible from the receiver side.
- Parameter estimation. Based on the multi-correlator values, the receiver shall perform an estimation of code delay, carrier phase, Doppler, and signal power of the reflected signal.
Altimetry
Altimetry, which estimates the height of the antennas over reflection points, is one possible application of a multi-constellation reflectometry system. The estimated height depends on the geometry and the path delays determined by measuring the time lapse between the correlation peaks of the direct and reflected signals (Figure 1) or by analyzing carrier phases of the direct and reflected signals.
The geometry changes in time only with the elevation angles if the antenna position and the reflection surface are fixed and the antenna setup is symmetric. The geometric situation is described in Figure 2.
Altimetry Observation Equations
The height can be estimated either directly or indirectly by using a code or carrier phase altimetry model. Those models are given by Equations (2) and (3), respectively. These equations include various hardware delays, which must be determined in advance to perform precise measurements.
All hardware elements in the RF cascade introduce a hardware delay. Therefore, the system needs to be calibrated at the beginning and end of the measurement period to determine the hardware delay and the behavior over time in order to compensate for it during analysis. This delay is indicated as ΔdhwDelay in the equations.
At the reflection point, a phase shift due to the reflection can occur. This phase shift is considered in equations with the parameter ΔdreflDelay. The value of this phase shift is treated as constant (for similar elevation angles) while performing the altimetry analysis.
Antennas can introduce further phases. This happens when the direct signal is received through a right-handed-circular-polarized (RHCP) antenna and the reflected signal is received with a left-handed-circular-polarized (LHCP) antenna. So, equipment operators must determine in advance if the antennas used in the setup will introduce a phase shift. This phase shift appears as ΔdantDelay in the equations and is also taken as a constant value because the behavior does not change over time.
In our case, we have a phase change of 0 degrees introduced by the antenna when using the RHCP port and 90 degrees when using the LHCP port. The reason for this is that such dual linear antennas often use a 90-degree hybrid that introduces this phase shift.
A geometric delay can be introduced by a nonsymmetrical setup of the antennas. This delay is represented by the parameter ΔdgeomDelay. In our case, we have chosen a symmetric setup (one antenna pointing towards zenith, the other towards nadir). In general, this delay varies over time with the geometric situation.
Equation (1) is a summation of all occurring delays which influence code and carrier phase measurements.
Δdcode/carr = ΔdhwDelay + ΔdreflDelay + ΔdantDelay + ΔdgeomDelay (1)
Equation (2) models code phase altimetry. This model can be solved on an epoch-per-epoch basis due to the fact that all necessary parameters for calculating the receiver height are given.
ΔS = 2 h sin(e) + Δdcode (2)
h Antenna height above ground (meters)
e Satellite elevation (radians)
ΔS (= ds in Figure 2) Additional path length of the surface reflected signal (meters)
Δdcode Representation for all occurring delays (meters)
Equation (3) considers the carrier phase altimetry model. This model can be solved indirectly by minimizing a cost function (4) in order to obtain the receiver height h and the carrier ambiguity difference ΔN.
Equation (3) (see inset photo, above right)
Δdcarr Introduced hardware delay difference of cascade elements and geometric placement of the antennas (radians)
ΔN Ambiguity difference (radians)
ΔΦ Unwrapped carrier phase difference (radians)
λ Carrier wavelength (meters)
The cost function in (4) is minimized in order to find the antenna height h and the ambiguity difference ΔN.
Equation (4) (see inset photo, above right)
Since we are working on a multi-frequency basis, we also want to introduce a method to work for that scenario. The multiple frequency approach gives a gain in height accuracy and more reliable integer ambiguity determination. The minimized cost function is similar to Equation (4) and is shown as Equation (5).
This introduces a dual-frequency cost function using E1 and E5b, which can easily be adapted for more frequencies. The only difference from the single-frequency approach is the need to work with the same height when minimizing the cost function. The elevation can be the same if the signals from one satellite are used.
Equation (5) (see inset photo, above right)
Demonstration Test RF Configuration
As mentioned earlier, the reflectometry system used in our demonstration is based on a software receiver with a dual frontend configuration, described in more detail in the publication by T. Pany et alia listed in the Additional Resources section near the end of this article.)
Figure 3 shows the reflectometry configuration that we used in our demonstration test. The master frontend samples the GNSS signals (in our case L1/E1/E5a/E5b) from the up-looking antenna, while the slave frontend samples from the down-looking antenna. Standard tracking is employed for the master channels, and they control via channel slaving the tracking channels of the down-looking antenna. The slaved channels output multi-correlator values, which are analyzed to derive the code and phase delay of the reflected signal with respect to the line-of-sight signal.
Generally, the receiver allows for real-time operation of multi-frequency reflectometry, but for our tests, data analysis in post-processing was more convenient. The ADC sampling is synchronized between the two frontends to better than ±14 centimeters. The frontend measures the remaining offset and outputs it at an accuracy of better than ±2 centimeters.
The RHCP + LHCP antennas used in our reflectometry experiment are generally L1/L2 antennas, but allow for reception of E5a/L5 and E5b signals through the L2 path. The signal degradation on L5/E5a is six decibels, and much less on E5b. The antenna manufacturer also offers L1/L2/E5a/E5b RHCP+LHCP antennas, but those were not available at the time of our experiments.
Signal Processing Methodology
We will now present the signal processing scheme and main parameters of the software receiver, the output of which was used for further analysis.
Master Slave Tracking. The receiver has a master-slave tracking algorithm for tracking direct and reflected signals either from file or directly from the USB interface of the RF frontend. The slave channel uses code and carrier NCO parameters from the direct tracking loop. These parameters assist the Slave channel to keep track of very weak and fluctuating reflected signals.
Figure 4 shows the general structure of the software receiver’s tracking channel, which uses intermediate frequency (IF) samples as input. The green NCO tracking loop keeps track of the signal by using delay/frequency/phase lock loop (DLL/FLL/PLL) values from its own discriminators, by employing vector tracking position/velocity/time (PVT) parameters, or by incorporating the NCO values from a master channel. Outputs of a tracking channel are the pseudoranges, Doppler, phases, and correlator values.
Multi-Correlator. The parameters for the slave multi-correlators can be set individually for each signal type. The list of multi-correlator parameter values used in this and other demonstrations can be found in Table 1.
Post-Processing Dump Log Files. The outputs from the tracking channels are saved in scientific ASCII log files. These dump log files are used for post analyzing purposes. In the case of altimetry, these dump log files are used for calculating code and carrier phase difference in order to calculate the antenna height by applying different models. Further intermediate data, for example, satellite and receiver position, are written into raw log files, which can also be used for post-processing analysis.
Altimetry Analysis
In this section we will explain how altimetry analysis can be performed using the output of the software receiver and the methods needed to achieve acceptable results. The altimetry analysis is currently done with MATLAB scripts (possible in real-time with the application programming interface). The overall schematic process for a single frequency altimetry approach can be seen in Figure 5. A dual-frequency approach with E1 and E5b is also shown later in this article.
Analyzing Multi-Correlator Log Files. First the *.dumplog and *.rawlog files, where the multi-correlator values and satellite/position information of all satellites are stored, are read in. Valid correlator values are filtered and synchronized with the rawlog data of the current processed satellite. After synchronization, the code correlation peaks will be determined in order to calculate code and carrier phase differences.
Polyfit for AltBOC Code Measurements. The software receiver implements a post-correlation FFT-based multi-correlator, which is described further in the article by C. Stöber et alia (Additional Resources). In the context of reflectometry, this is known as a Doppler/delay map and is computed with respect to the direct signal. Figure 6 shows part of the reflected AltBOC correlation function. In this static scenario, only correlator values at zero Doppler offset are considered.
The determination of the peak position for calculating the code phase difference is accomplished by applying a second-order polynomial in the code phase direction at zero Doppler if narrow-correlator spacing is used (e.g., for the AltBOC signal). The curve is fitted around the largest calculated correlator value of the reflected signal, and the peak position of the polynomial is determined by searching the maximum. If the antenna height over the reflection point is very low, only AltBOC code measurements make sense because of improved accuracy compared to C/A code.
Carrier Phase Retrieval. In order to determine carrier phase differences, another peak search method is used. In this case, the peak is found by a triangle interpolation method. The maximum peak position is searched first. Then the gradients on the left and right sides of the next correlator position are determined. The higher gradient is used to interpolate from both sides to the center in order to achieve a more accurate peak position. The interpolation is done for the real (I) and imaginary (Q) part separately.
At this peak position the energy is highest in order to achieve the most accurate phase angle. The phase angles of direct and reflected signals are subtracted in order to obtain the carrier phase difference. This carrier phase difference is used to perform altimetry analysis by applying the carrier phase altimetry model.
System Calibration
For system calibration we take a signal from one antenna, split it up, and forward it to both front ends (RF switch position RF cal. in Figure 3). This means that the direct signal is received at exactly the same time with both frontends. Therefore, the carrier phase difference and code phase difference should also both be zero. If these values are not zero, a delay is present in the system that must be estimated.
The system calibration needs to be performed for code and carrier phase separately, and the delay is frequency-dependent. This calibration delivers information about internal system delays, which have an influence on the applied models. These delays might be time-dependent because the delays occurring in the electrical elements may be temperature- dependent. Therefore, a calibration is done before and after the measurement sequence in order to obtain information about the behavior over time.
Figure 7 shows a carrier phase calibration plot. In the upper figure, the two flat parts at the beginning and end are the calibration periods (in red). In between, signals reflected from the lake surface are shown. The lower figure shows only the calibration parts in black and a fourth-order polynomial is fitted to these to model the hardware delay as a function of time (red line).
Various methods can be chosen that best fit the delay. This can be a polynomial of order zero up to sixth order. Higher orders are also possible but in our experience no gain from higher orders can be achieved. As mentioned previously, these delays are then taken into account when applying the altimetry models.
Proofing System Calibration. We can check to see if we have estimated our delays correctly. We know that the height of the antenna is zero for the calibration phase because no phase difference exists. If the internal delays are estimated correctly, the altimetry model with delays should deliver a height of nearly zero for the calibration sequence. If not, we can assume that the altimetry estimation for the measurement sequence has an offset due to internal delays.
Our goal here is to model the delay as well as possible to minimize offsets in height estimation. Such height estimations are applied on the calibration sequences as shown later in Figure 14 where they are indicated with blue and green lines.
Cycle Slip Detection and Compensation
When carrier phase unwrapping is performed using a simple MATLAB function, cycle slips of 2π radians can often be seen in the carrier phase difference. Such cycle slips can occur if the reflection is not ideal, which could be caused by a rough reflecting surface or a weak reflection.
If cycle slips occur, we cannot apply altimetry models because the cycle slips influence the results drastically. Therefore, cycle slips must be detected and eliminated before applying the model. We can do this by searching the lowest number of samples where a certain threshold is exceeded. This threshold is near a cycle slip of 2π radians and can be adjusted. In Figure 8 in red a simple cycle slip of 2π radians can be seen.
Sometimes double cycle slips occur. These are two cycle slips within a sequence. Attempting to detect the cycle slips in two iterations is not ideal because doing so doubles the threshold, which is lower than 2π radians. Therefore, we try to detect rarely occurring double cycle slips with a defined threshold slightly below 4π radians before detecting the individual cycle slips. Figure 9 shows a detected double cycle slip in green.
Once a cycle slip has been detected, we know the number of samples in which the cycle slip occurs. If we take equally distributed steps, with the number of samples from 0 to 2π radians, the possibility exists of introducing very sharp and high peaks when adding to the existing signal.
We attempt to smooth the cycle slips as much as possible by weighting the added parts, which depend on the distance between steps in a cycle slip (the greater the step, the greater the compensation). At the end of the cycle slip, the compensation is always 2π radians.
Situations may also arise in which the phase makes a 2π radians jump without introducing a cycle slip. This phase change will also be detected as a cycle slip; so, we then have to determine if it really is. This is done by averaging a multiple of the number of samples where the suspected cycle slip occurs. This averaged level before and after the cycle slip is compared to a reduced cycle slip detection threshold because the average can be lower due to adverse phase point distribution. This method prevents the introduction of cycle slips by the method itself, as seen in blue in Figure 10.
A cycle slip may also occur over numerous samples, which makes them more difficult to detect. In this case we can increase only the parameter that defines the number of samples within which we look for a cycle slip. When using this method, we recommend starting with a low number and increasing it iteratively.
In Table 2 we have listed common values for the method to detect and remove all cycle slips in most cases. If all cycle slips are not detected, the parameters must be adjusted. We describe the parameters in Table 2 thus:
- Points are the number of samples looked for in a cycle slip.
- Threshold is the level at which a cycle slip is detected, value is lower than 2π radians.
- AvPtsFac is the factor for averaging when detecting wrong cycle slips.
- Threshold % is the reduction of the threshold when trying to detect incorrect cycle slips.
A complete example of this “non–cycle slip” detection method can be seen in Figure 11. The upper plot shows the unwrapped original signal and the lower plot shows the non-slipped carrier phase difference.
Applying the Carrier Phase Altimetry Model
Once a data sequence is prepared, the carrier altimetry model can be applied for the single-frequency case.
When applying the model we must distinguish between a float solution and a fixed solution. To do this, we first estimate the antenna height and the integer ambiguity by minimizing the cost function (4) with MATLAB’s “fminsearch” function. After this step, we have the float solution because the carrier ambiguity is a real-valued number (float) at this moment.
Our experience indicates that if not all delays are known, the height of the float solution is more accurate than the height of the fixed solution. If all delays are known, however, the carrier ambiguity must be an integer multiple of 2π radians. So, the integer ambiguity is rounded in terms of cycles to the next integer and fixed. Then the model is applied a second time with fixed integer ambiguity in order to get the fixed height. If all delays are modeled correctly, this latter height solution should be the more accurate of the two.
Dual-Frequency Approach
With the dual-frequency approach, we try to find a more reliable and accurate height. The way to achieve this is partly similar to the single-frequency approach. We take the cost function described in Equation (4) and calculate a float solution for each frequency. Then the float ambiguities are rounded to integer ambiguities.
Now there are two fixed integer ambiguities for the two frequencies. Next the cost function described in Equation (5) is minimized with varying integer ambiguities. The cost function delivers a height with a certain cost. The height at minimized cost — or in other words, where the model best fits with the calculated integer ambiguities — is considered as the best solution. Later in the article, we show an example of this method.
Lake Experiment and Results
We conducted three experiments (AltBOC, Measurement 1, Measurement 2) on lakes near Graz, Austria, to demonstrate our achievements with these methods. Table 3 presents the experiment parameters. In this experiment we concentrate our measurements on the Galileo PRN11 satellite, focusing on the carrier phases of E1 and E5b as well as the E5 AltBOC code.
The signal strength and location are ideal for receiving a reflected signal over a water surface. Figure 12 shows the skyplot of the surrounding environment for one of the experiments.
Galileo AltBOC Code Results
Generally, code phase measurements are expected to be not as precise as the carrier phase measurements under smooth water surface conditions. However, for the wideband Galileo AltBOC signal on E5 we would expect a centimeter-level noise for altimetry measurements and decimeter-level degradation due to secondary multipath. The software receiver implements a unique method to track the AltBOC by a post-correlation method based on E5a and E5b data.
The achieved results for AltBOC code altimetry based on an epoch-per-epoch solution of Equation (2) can be seen in Figure 13. For this analysis, the data was processed with the software receiver and a coherently slaved multi-correlator. The coherent integration time was 2.56 seconds; such long integration times are possible for flat water surfaces. Each star in Figure 14 corresponds to a different set of multi-correlator values.
The noise is at the centimeter level. The overall accuracy for this first attempt is at the level of about half a meter.
Single-Frequency Carrier Phase Altimetry Results
The carrier phase altimetry was performed on E1 and E5b for the data and pilot channels. Figure 14 shows four plots for the pilot channel on E1. The first plot gives an overview of the whole signal carrier phase difference. The next three plots show the sequences that we selected to analyze. The blue and green lines in the second and fourth plots indicate the calibration sequences at beginning and end of the signal sample.
These calibration sequences should yield an estimated antenna height of zero if calibration has been done correctly. The antenna height for the E1 pilot channel for the blue sequence indicates 0.0006682 meters and for the green, 0.000633 meters. This is a good indication that calibration has been performed correctly.
For hardware delay calibration, we used a polynomial fit of fourth order, shown previously in Figure 8. The black parts in Figure 14 correspond to changing the RF switch manually from Calib. to Meas. as shown in Figure 3. In the final product, this will be done with a USB-based RF switch within a few milliseconds.
The red sequences in Figure 14 comprise our measurement sequence. This sequence delivers information about antenna height and reads 2.7103 meters (fixed solution). Here we must note that the estimated float height matches the true height very well. We also should note that the measured true height can vary a bit because it is difficult to determine the true height of water with small waves.
The values and also the plots for the E1 data channel are nearly identical. The blue calibration reads zero meters, the green calibration reads 0.000398 meters, and the red estimated height reads 2.7103 meters.
The same single frequency altimetry analysis is done with the E5b signal. For calibration, a first order fit was chosen. Results for the E5b pilot signal can be seen in Figure 15. The blue and green calibration sequences read -1.14*10-5 meters and -0.0773 meters. The red estimated height (fixed solution) reads 2.7262 meters.
The values and plots for the E5b data channel are nearly identical as for the pilot. The blue and green are not exactly zero. The blue calibration reads -1.1*10-5 meters and the green reads -0.075 meters. The estimated height is 2.7262 meters (fixed solution).
The single frequency approach delivers an estimated height for each frequency. With only single-frequency approaches, we cannot determine the best antenna height estimation. Therefore, we need a dual-frequency approach that combines these single-frequency approaches.
Dual-Frequency Results
This approach uses the cost function described in (5). In Figure 16, the results of the dual frequency approach of the E1 and E5b pilot channels are shown. The blue line indicates the carrier phase difference. The black lines are the varied integer ambiguities and resolved models. The green line is the model which fits best to both frequencies.
The estimated height derived by the dual-frequency approach is 2.7183 meters, which is taken as the most reliable and accurate solution within this analysis. With this result, we have an offset of 8.3 millimeters from the true height (again, it is difficult to measure the true height to millimeter accuracy on a waving water surface).
At this point we also want to show the results of a second dual-frequency analysis to get more confidence in this method. Figure 17 shows a second measurement with the true height of 2.635 meters, which was obtained by mounting the antenna right beside the lake and using a measuring tape. This measurement was more difficult due to the presence of higher waves of about three to six centimeters during the observation period.
The red lines in Figure 17 show the phase model based on the integer ambiguities from the respective single-frequency solutions. They are incorrectly resolved and only the dual-frequency approach (green lines) gives the correct ambiguities and the correct height. This is related to the fact that, in the single-frequency approach, the height parameter can compensate any incorrect ambiguity.
This is not the case for two frequencies and only certain ambiguity combinations are near the measurements. This limits the ambiguity search space drastically, in the same way as dual-frequency RTK positioning is superior to single-frequency RTK positioning.
The estimated heights of the second run are:
- Estimated height of E1 (fixed): 2.6352 meters
- Estimated height of E5b (fixed): 2.7915 meters
- Estimated height by dual-frequency model: 2.6418 meters
- True height: 2.635 meters
The offset of the estimated height using the dual-frequency approach is about 6.8 millimeters from the measured true height in this case.
Water Wave Analysis
The carrier phase analysis provided additional information about water waves. In the second carrier phase experiment we had waves about three to six centimeters in height (peak-peak). We estimated that the waves had a frequency of one to two waves per second.
This estimation matches the results in Figure 18. This plot shows the impact of waves on the water surface on the carrier phase difference, with a wave height of about three to four centimeters and a period of about 0.75 seconds.
From these results, some additional information can be derived. If waves from spatially different reflections points on the water surface can be correlated, the wave’s speed and direction — which is related to the surface wind — can be estimated.
Outlook and Applications
Given a high-performance receiver, the new Galileo signals offer the possibility of measuring distances from reflecting surfaces very precisely. Decimeter accuracy might be achievable with Galileo AltBOC code measurements, while millimeter accuracy for a solution with fixed carrier phase ambiguities may be possible.
The technique should be applicable when the reflecting surface is flat (compared to the carrier wavelength) on an area corresponding to the first Fresnel zone (which in our experiment is several meters). In principle, one could think of monitoring large man-made structures such as water dams and skyscrapers. Monitoring of landslides or ice thickness might also be possible. GNSS signals are also reflected from airplanes or other vehicles. If these are strong enough, a precision passive radar could be built using the methods described in this article.
Acknowledgment
The authors would like to acknowledge Daniel Llorens for helpful discussions of the hardware delay inside a dual port RHCP/LHCP antenna. The support of the University of Applied Sciences is also acknowledged.
Additional Resources
1. Borre, K., and D. M. Akos, N. Bertelsen, P. Rinder, and S. H. Jensen, “A Software-Defined GPS and Galileo Receiver: A Single-Frequency Approach,” Springer, ISBN: 0-8176-4390-7, 2007
2. Irsigler, M., “Characterization of Multipath Phase Rates in Different Multipath Environments,” GPS Solutions, Vol. 14, pp. 305-317, 2010
3. Pany, T., and N. Falk, B. Riedl, T, Hartmann, G. Stangl, and C, Stöber, “Software GNSS Receiver, An Answer for Precise Positioning Research,” GPS World, Vol. 3, No. 29, September 2012
4. Soulat, F., and O. Germain, G. Ruffini, E. Farrés , T. Sephton, I. Raper, and S. Kemble, STERNA – A Feasibility Study of PARIS Tsunami Detection – Final Report, STARLAB, ESA/ESTEC Contract 19016/05/NL/JA, 2005
5. Stöber, C., and F. Kneissl, B. Eissfeller, and T. Pany, T., “Analysis and Verification of Synthetic Multicorrelators,” Proceedings of the 24th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2011), Portland, OR, pp. 2060-2069, September 2011
