Working Papers • January/February 2013
Galileo Altimetry Using AltBOC and RTK TechniquesReflectometry systems analyze reflected GNSS signals to investigate properties of the reflecting surface and to derive various parameters from that surface. The emergence of additional GNSS systems has increased the availability of multiple signals on multiple frequencies, which provides the opportunity for developing new techniques for applying reflectometry. This article describes an E1/E5a/E5b groundbased altimetry application that has been tested with Galileo signals. The innovative aspect of this system resolves carrier phase ambiguities as in RTK positioning to achieve millimeter accurate surface height measurements, including the use of Galileo AltBOC code measurements for instantaneous decimeter accurate height values. New applications enabled by this technique may include monitoring of large buildings, landslides, or precision passive radar.
Share via: Slashdot Technorati Twitter Facebook 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 (GNSSR) is nowadays a wellestablished method to remotely sense various properties of the Earth’s surface via reflected GNSS signals. Groundbased, airborne, and satellitebased GNSSR 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, FPGAbased receivers, and software receivers, which particularly benefit postprocessing techniques. To date, most of these have employed the GPS C/A signal for obvious reasons. However, the next generation of GNSSR 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 realtime on a conventional PC.
GNSSR Demonstration Application We used a multifrequency and multisensor software receiver with appropriate GPS and Galileo dual RF frontends that operated in realtime. The receiver supports GNSSR by means of the following elements:
Altimetry 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 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 Δd_{hwDelay} 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 Δd_{reflDelay}. 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 righthandedcircularpolarized (RHCP) antenna and the reflected signal is received with a lefthandedcircularpolarized (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 Δd_{antDelay} 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 90degree 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 Δd_{geomDelay}. 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. Δd_{code/carr }=_{ }Δd_{hwDelay} +_{ }Δd_{reflDelay} +_{ }Δd_{antDelay} +_{ }Δd_{geomDelay }(1) Equation (2) models code phase altimetry. This model can be solved on an epochperepoch basis due to the fact that all necessary parameters for calculating the receiver height are given. ΔS = 2 h sin(e) + Δd_{code} (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) Δd_{code} 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) Δd_{carr} 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 multifrequency 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 dualfrequency cost function using E1 and E5b, which can easily be adapted for more frequencies. The only difference from the singlefrequency 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 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 uplooking antenna, while the slave frontend samples from the downlooking antenna. Standard tracking is employed for the master channels, and they control via channel slaving the tracking channels of the downlooking antenna. The slaved channels output multicorrelator values, which are analyzed to derive the code and phase delay of the reflected signal with respect to the lineofsight signal. Generally, the receiver allows for realtime operation of multifrequency reflectometry, but for our tests, data analysis in postprocessing 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 Master Slave Tracking. The receiver has a masterslave 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. MultiCorrelator. The parameters for the slave multicorrelators can be set individually for each signal type. The list of multicorrelator parameter values used in this and other demonstrations can be found in Table 1. PostProcessing 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 postprocessing analysis.
Altimetry Analysis Analyzing MultiCorrelator Log Files. First the *.dumplog and *.rawlog files, where the multicorrelator 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 postcorrelation FFTbased multicorrelator, 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 secondorder polynomial in the code phase direction at zero Doppler if narrowcorrelator 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 The system calibration needs to be performed for code and carrier phase separately, and the delay is frequencydependent. This calibration delivers information about internal system delays, which have an influence on the applied models. These delays might be timedependent 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 fourthorder 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 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:
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 nonslipped carrier phase difference.
Applying the Carrier Phase Altimetry Model 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 realvalued 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.
DualFrequency Approach 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 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 The achieved results for AltBOC code altimetry based on an epochperepoch 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 multicorrelator. 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 multicorrelator values. The noise is at the centimeter level. The overall accuracy for this first attempt is at the level of about half a meter.
SingleFrequency Carrier Phase Altimetry Results 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 USBbased 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 singlefrequency approaches, we cannot determine the best antenna height estimation. Therefore, we need a dualfrequency approach that combines these singlefrequency approaches.
DualFrequency Results The estimated height derived by the dualfrequency 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 dualfrequency 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 singlefrequency solutions. They are incorrectly resolved and only the dualfrequency approach (green lines) gives the correct ambiguities and the correct height. This is related to the fact that, in the singlefrequency 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 dualfrequency RTK positioning is superior to singlefrequency RTK positioning. The estimated heights of the second run are:
The offset of the estimated height using the dualfrequency approach is about 6.8 millimeters from the measured true height in this case.
Water Wave Analysis 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 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 manmade 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
Additional Resources ManufacturersThe software receiver used in the reflectometry demonstration trial was the SXNSR from IFEN GmbH, Poing, Germany. The antennas were from AntCom, Torrance, California, a wholly owned subsidiary of NovAtel, Inc.: uplooking antenna, RHCP, 1G1215RLALXS1; downlooking antenna, LHCP, 1G1215RLALXS1.Copyright © 2016 Gibbons Media & Research LLC, all rights reserved. 
