1. Introduction

PROGENY (PROvision of Galileo Expertise, Networking and support for International Initiatives) is a research and technological development project launched by the European GNSS Supervisory Authority (GSA), in the frame of the 6th Framework Programme. PROGENY consists of a series of activities supporting the innovation and international initiatives in relation to the Galileo programme. In particular, the project has established a platform for scientific and technological cooperation with different regions worldwide, and has run a set of targeted studies in cooperation with international partners.

This paper presents the results of the study performed by the University of Calgary Department of Geomatics Engineering related to the definition of a method for LEO satellite attitude determination, using Multiple Carrier Ambiguity Resolution (MCAR).

GNSS-based attitude determination is accomplished by kinematic carrier-phase GNSS techniques. Namely, a number of short baselines are established on the vehicle with known coordinates in the vehicle body frame. Carrier-phase GNSS is used to determine local level frame (east, north, vertical) components of the same baselines and the knowledge of the baselines in both frames is then used to establish the rotation angles between the two frames.

In order for attitude to be determined precisely, the GNSS baseline components must be estimated using fixed carrier-phase ambiguities. There are presently two main approaches to dual-frequency kinematic double-difference ambiguity resolution, and these two methods can be generalized to the case where modernized GPS and Galileo provide additional observations on additional frequencies.

The first approach is to estimate a position based on pseudorange measurements and also form the wide lane observable by differencing the L1 and L2 phase measurements. Either the pseudorange or the pseudorange-derived position can be used to provide an initial estimate of the widelane ambiguity. The widelane can generally be resolved quickly over short baselines typically associated with attitude determination (i.e., order of a few metres). Once resolved, the next step is either to use the fixed widelane phase range, or alternately the fixed widelane position estimate as a starting point to estimate a float solution for the L1 ambiguity. Because a fixed widelane phase range is more precise than a pseudorange measurement, it becomes possible to estimate the L1 float ambiguity with sufficient confidence to allow it to be resolved quickly. The sequence of steps between pseudorange, widelane, and L1 has led to this method being called “cascading.” Triple-frequency variations of this algorithm have been proposed for both modernized GPS and Galileo where an additional step is added to make use of the third frequency to form an even longer wavelength widelane observable before cascading down to the widelane observable. These methods are generally referred to as either three carrier ambiguity resolution (TCAR), multiple carrier ambiguity resolution (MCAR), or simply cascading integer resolution (CIR) methods [1].

The second approach to multiple frequency kinematic ambiguity resolution involves using Teunissen’s Least squares AMBiguity Decorrelation Adjustment (LAMBDA) method to determine an optimal linear combination of a set of float ambiguities for the purposes of ambiguity resolution [2, 3]. In this method, any set of float ambiguities may be estimated and then a linear combination is found that minimizes the correlation between the set of ambiguity states.

In both these methods, either geometry-free or geometry-based processing may be used. In geometry-free processing, double-difference observations from each satellite are treated independently until ambiguities are resolved at which point the fixed phase range is used to compute a position solution. In Geometry-based processing, the float ambiguities and the baseline vector are estimated together and the baseline vector is then improved once the ambiguities are resolved [4, 5].

In the geometry-based case, the inclusion of observations from more satellites, or satellites from more than one system (e.g., GPS and Galileo) has been previously shown to result in improved ambiguity resolution performance. Likewise, the addition of a third frequency has been shown to improve geometry-free ambiguity resolution [6].

The major contribution of this paper is that it presents, to our knowledge, the first large-scale comparison through simulation of geometry-based LAMBDA to geometry-based cascading for any application (terrestrial or orbiting, positioning or attitude determination). Most previous work comparing the two methods was theoretical only and addressed only the geometry-free case [7] though the geometry-based case is addressed in [8].

This paper aims to demonstrate the effectiveness of various multicarrier ambiguity resolution methods for the purpose of attitude determination of a Low Earth Orbiting (LEO) satellite. The remainder of the paper is divided into three sections. Following a brief review of carrier phase GNSS in Section 2, a set of simulation scenarios are defined in Section 3 in order to assess several ambiguity resolution methods for two typical LEO orbiting users tracking GNSS satellites. Performance is assessed by simulating the use of Galileo alone, and Galileo and GPS together with one, two and three frequencies. The results of this series of simulations are presented in Section 4 with an emphasis on the performance of each method in each scenario in terms of probably of correct ambiguity resolution.

2. Background

GNSS data processing for attitude determination can be divided into four steps listed in what follows. The first three steps are identical to the process of carrier-phase GNSS positioning without attitude determination. The difference is that for attitude determination the baseline(s) being determined connects two or more points on a vehicle with known coordinates in the body frame of the vehicle. Attitude determination is then implemented in the fourth step.

(1) Float Ambiguity Solution
This is the process of using the available observations to estimate a real-valued (float) estimate of the carrier phase ambiguities. These ambiguity estimates may, if necessary, be filtered over time in order to reduce their uncertainty.

(2) Ambiguity Resolution
This is the process of resolving the float ambiguities to integer values. The output from the ambiguity resolution process is a set of integer carrier phase ambiguities. It is noted, however, that the ambiguities are not necessarily guaranteed to be correct.

(3) Fixed Ambiguity Baseline Solution
The integer ambiguities are used, along with the corresponding carrier phase measurements, to generate an estimate of the relative position vector between the two receivers involved in the double-difference. This relative position vector is usually called the baseline vector, or simply “baseline.”

(4) Convert Baseline Solution to Attitude Solution
This is the process of using the known baseline vectors in the body frame of the vehicle (spacecraft) and baseline vectors obtained from the GNSS solution in the local-level frame to estimate the attitude of the vehicle.

Each of these steps is discussed in more detail in the following sections.

2.1. Float Solution

The primary objective of the float solution is to obtain an initial, real-valued, estimate of the carrier phase ambiguities. The actual implementation of the float solution will depend on the data processing strategy adopted, but these strategies can be broken down into two categories: geometry-free and geometry-based [9, 10].

2.1.1. Geometry-Free Approach

The geometry-free approach is not concerned with the position of the receiver and instead aims to estimate the double-difference range and ambiguity to each satellite, along with any significant systematic errors. The pseudorange and carrier phase measurements made on one or more frequencies are the inputs into the system. The state vector usually consists of the range to the satellite, the ambiguities to be estimated and optionally an ionospheric error term. The latter is usually only included when the residual double-difference ionospheric error is nonnegligible. For the case at hand however, because the receivers are all located within a few metres of each other, the ionosphere term can be safely neglected.

The fact that each satellite is treated separately is both an advantage and a disadvantage. It is an advantage because it provides a relatively simple implementation and it does not depend on the number of satellites in view, nor their distribution in the sky. It is a disadvantage because it does not exploit the fact that the measurements to all of the different satellites are related via the position of the receiver. In other words, no information is shared between filters estimating each double-difference ambiguity, which generally degrades performance. A further disadvantage is that the pseudorange errors—particularly multipath—can significantly degrade reliability.

2.1.2. Geometry-Based Approach

In contrast to the geometry-free approach, the geometry-based approach explicitly estimates the baseline vector between the two receivers along with the ambiguities and any other systematic errors. Again, for the short baselines involved in this application, these systematic errors need not be considered. The state vector is usually divided into two components, a vector of ambiguities to the various satellites and the remaining states such as position, velocity, and so forth. In this way, all of the observations are linked together via the position information which provides geometric strength to the solution. This also implies that the ambiguities for all the satellites are estimated together, instead of on a satellite-by-satellite basis, as with the geometry-free approach.

The main disadvantage of the geometry-based approach is that it is dependent on the number and distribution of the satellites being tracked. As such, if the number of satellites tracked decreases below four, or if the distribution of satellites in the sky is unsatisfactory, performance will suffer. That said, for the application at hand, and for the planned number of GNSS satellites in orbit (Galileo with/without the addition of GPS), this is not expected to play a major role.

2.2. Carrier Phase Ambiguity Resolution

In practice, there are several strategies of ambiguity resolution, and in the present section, the implementation details of each will be described. As discussed earlier, ambiguity estimation techniques can be broadly classified as either geometry-based or geometry-free, depending on whether or not baseline vector components are estimated simultaneously with the ambiguity states. Ambiguity resolution methods can be further classified as instantaneous or filtered depending on whether or not more than a single epoch of observations is used in the estimation process. Finally, multifrequency methods can be further divided between those that use specific linear combinations of the various frequencies (referred to as cascading methods in this paper) and those that attempt to estimate the optimal combination on the fly (LAMBDA).

2.2.1. Widelaning, Cascading Methods and Lambda

In the simplest sense, ambiguity resolution can be accomplished by comparing an absolute measurement (a code pseudorange) with a relative measurement (a phase measurement) to determine the bias between the two (the ambiguity). Conceptually, this requires that the code pseudorange be accurate enough that one can confidently determine the phase ambiguity. In practice, this means that the uncertainty on the code measurement must be significantly less than the carrier wavelength so that the code measurement can place you definitively in a particular carrier phase cycle. In low frequency radio-navigation systems this is relatively easy but in GNSS, the wavelengths are short and the code measurements are relatively noisy. For example, the GPS L1 wavelength is approximately 19 cm but a typical double-differenced code measurement may have errors on the order of 50 cm or more.

Widelaning
If measurements on more than one frequency are available, it is possible to form linear combinations of these measurements that have larger effective wavelengths. Specifically, a widelane (WL) observation is formed when two phase observations are subtracted from each other. The resulting linear combination has a frequency equal to the difference between the frequencies of the two original observations, With current dual-frequency GPS, a widelane observation can be formed. The resulting combination has an effective wavelength of approximately 86 cm, making it more reasonable to estimate the widelane ambiguity from a between receiver single-differenced code pseudorange measurement. Once the widelane ambiguity has been determined, a fixed widelane phase range can be formed that can then be used to estimate the L1 ambiguity.
With the addition of a third frequency on GPS and the deployment of Galileo, new widelane phase combinations are possible. From (2) it should be noted that a widelane wavelength is inversely proportional to the difference in the frequencies of the two signals involved. For this reason L2 and L5 in GPS and E5a and E5b in Galileo can be used to form very long wavelength widelanes that are often referred to as the extra-widelane (EWL) combination. For GPS L2-L5, the EWL wavelength is 5.861 m. For Galileo E5a-E5b it is 9.768 m. The existence of these very long wavelengths forms the basis for a number of integer cascading ambiguity resolution algorithms. Table 1 lists the modernized GPS and Galileo frequencies and wavelengths while possible widelane combinations are listed in Table 2.

Table 1: Characteristics of modernized GPS and GALILEO open service frequencies and signals.

Table 2: GPS and Galileo widelane wavelengths in metres.

Cascading Algorithms
The use of the L2-L5 widelane as the first step in a cascaded ambiguity resolution for GPS has been proposed by many researchers, particularly for the geometry-free case [11, 12], and similar algorithms have been proposed for Galileo [13, 14]. When the EWL ambiguity is resolved, the resulting phase-range is then used to estimate a shorter wavelength widelane ambiguity and this process is continued until the L1 ambiguity is resolved at which point a carrier-phase position solution is computed from the fixed L1 phase measurements. GPS methods have traditionally been called CIR (cascaded integer resolution) methods while Galileo methods have been referred to as TCAR or MCAR (triple- or multicarrier ambiguity resolution) methods.
Most additional previous work focuses on finding better linear combinations of the three phase measurements (other than EWL, WL, L1) to increase the likelihood of successful ambiguity resolution in the presence of large differential atmospheric errors for long baseline surveying applications [15–19]. There have been some attempts to use triple-frequency GPS and triple-frequency Galileo simultaneously in geometry-based ambiguity resolution techniques [4, 5] including the use of the two common frequencies (L1 and L5) to obtain an additional double-differenced measurement between the two systems [6].

The Lambda Approach to Multifrequency Ambiguity Resolution
The Least squares AMBiguity Decorrelation Adjustment (LAMBDA) method is a generic method for ambiguity resolution [20]. The method can be applied to any set of float ambiguities that have been jointly estimated (meaning that the float ambiguities share a covariance matrix). The method can either be applied to the multifrequency geometry-free case or to single- or multifrequency geometry-based cases. Its excellent performance on short baselines is well known and is described in [21]. In the geometry-free case two or three ambiguities are being estimated for a single double-difference satellite pair. In the geometry-based case, several ambiguities and baseline vectors are being estimated. In both cases, the key to the use of the LAMBDA method is in the fact that double-differenced ambiguities that are estimated together are usually highly correlated. Earlier search-based techniques would estimate float ambiguities for either of the two above-mentioned cases and then search a large volume of possible ambiguity sets around the float solution for the best fitting fixed ambiguity solution. Unfortunately, due to the high correlation of the ambiguity states, the search volume is typically very large making these methods very time consuming [22]. The LAMBDA method also employs a search, but prior to the search, it attempts to find a linear transformation of the ambiguity set being estimated that decorrelates the ambiguities from each other while maintaining their integer nature. Unlike the cascading methods mentioned earlier, which rely on specific linear combinations of the phase measurements to facilitate ambiguity resolution, the LAMBDA method uses the covariance matrix of the float ambiguity solution to find the optimal linear combination to facilitate ambiguity resolution. It has been previously shown that the various cascading schemes (where the linear combinations allowed are restricted the EWL, WL, and L1) are theoretically suboptimal compared to the LAMBDA- derived linear combinations in the geometry free case [7]. The LAMBDA method’s application to the multifrequency ambiguity resolution problem has been studied extensively in [6, 23–28].

2.2.2. Geometry-Based AR and Assessing Probability of Correct Fix

The general mathematical model for geometry-based carrier-phase GNSS positioning can be described as follows:

where is the observation vector, is the state vector to be estimated in the filter, comprised of the ambiguities sub-vector and the other parameters to be estimated, for example, position, ionospheric and/or tropospheric estimates, is the design matrix for the full state vector , and and are the corresponding design matrices for states and , is the measurement noise, generally assumed to be white and to follow a Gaussian normal distribution. There are three steps to resolve the ambiguities to their integer values, which can be described as follows.

(1) Estimate the ambiguities as real values, , with the other state parameters, . The integer nature of the ambiguities is ignored and the estimates are referred to as float estimates. The corresponding covariance matrix of the errors of the state vector can be the partitioned into the covariance of the ambiguity errors, the covariance matrix of the other states, and a cross-covariance term

(2) Determine the integer value of ambiguities based on the float estimates . Various methods have been proposed in the process of integer fixing, as discussed earlier.

(3) Compute the fixed estimates of parameter based on the earlier fixed integer ambiguities. The fixed estimates and their variance can be formulated from the solution in the aforementioned Step 1 as

where the covariance of the fixed estimates is valid only under the assumption that the ambiguities have been fixed to their correct values. In the aforementioned Step 2, the process of obtaining the integer values of the ambiguities can be defined as a mapping of the real space to the integer space. Then the probability that a given integer vector is equal to a particular integer vector , and can be assessed as

where is the probability density function (PDF) of the float ambiguities [29]. This defines the probability of correctly resolving the ambiguities, or the probability of correct fix (PCF).

However, it is difficult to quantify (8) numerically because of the complexity of the so-called pull-in region and the computational load of the integration process. Thus a simplification or approximation is required. There have been various bounds proposed to approximate the PCF [29]. A method of ambiguity bootstrapping is widely used and adopted to determining a lower bound of the PCF [29–31]. The evaluation is based on the following expressions:

In (9), is the bootstrapped integer ambiguity vector, is the conditional standard deviation of ambiguity conditioned on the previous ambiguities, and describes the area under the normal distribution.

There are many methods to fix the float ambiguities to integer values in Step 2. In the case of ambiguity rounding, the conditional nature of the conditional standard deviations are ignored and an estimate of the PCF can be obtained from the variances of the float ambiguities. It has been previously observed that (9) is not invariant to ambiguity parameterizations [29–32]. In the remainder of this paper the PCF is quantified in terms of Probability of Incorrect Fix (PIF), which is equal to one minus PCF. Likewise, the lower bound of probability of correct fix becomes an upper bound on the probability of incorrect fix. In this case “good” performance is represented as a small upper bound (significantly less than 1) while “poor” performance is represented by an upper bound on PIF that is close to one (suggesting that incorrect ambiguity resolution is likely).

2.3. Fixed Ambiguity Baseline Solution and Transformation to Attitude Solution

Generally speaking, if the carrier phase ambiguities are resolved to their correct integer values, the accuracy of the baseline estimate obtained from (6) is on the order 1–3 cm in each coordinate direction. However, if one or more ambiguities are resolved to an incorrect integer value, then the baseline solution will contain errors on the order of the carrier phase wavelength (e.g., approximately 19 cm at L1). This is why the ambiguity resolution process is so critical. Fixed position accuracy can be obtained from the covariance matrix of the fixed solution obtained with (7).

The baseline solution obtained from (6) can then be used to compute the attitude of the vehicle. We denote this solution in the local level frame as . Assuming the baseline vector between any two antennas is also known in the body frame, , then the following relationship will hold:

where is the rotation matrix that transforms a vector in the body frame into the local level frame. The assumed convention for the rotation matrix is

where , , and are respectively the roll, pitch and azimuth of the vehicle, and , , and are rotation matrices about the primary , secondary , and tertiary axes respectively. It is noted that other orders of rotations can also be used without loss of generality.

Assuming two noncolinear (nonparallel) baseline vectors are available, (11), (12) allow for the estimation of the attitude parameters using the known baseline vectors in the body frame and the measured baseline vectors in the local level frame. Mathematically,

where are the errors in the baseline vector estimated from the GNSS data. Estimation of the attitude parameters can be performed using least-squares or Kalman filtering. The covariance matrix of the attitude solution can be obtained from the covariance matrices of baselines though covariance propagation [33].

3. Simulation Test

In this section, the accuracy of a short baseline is simulated for two different LEO satellites. The accuracy is highly dependent on the ability to correctly resolve integer ambiguities, and as such the first step is to assess the ambiguity resolution performance.

3.1. Simulation Parameters

The signal characteristics of modernized GPS and Galileo have been investigated in great detail by many researchers [34]. Only the GPS civilian signals and Galileo open service (OS) signals will be considered in this study, their characteristics are listed in Table 1.

Note that the open services of modernized GPS and Galileo have two frequencies in common, but the third frequency on each is unique. Because of this, it could be possible to form double-difference observations between the two systems using L1 and E1/E2 phase measurements and also L5 and E5a [6]. However, this is not possible for L2 and E5b. Double-differencing between systems was not conducted in this study and for all three frequencies a base satellite was assigned for each system.

No official documents on the range accuracy of future GPS and Galileo signals have been released but several studies quantify the signal performances of future GNSS systems (e.g., [35]). The code and carrier observations are affected by systematic errors and random noise. However, for attitude determination applications, where the baseline is very short, the only error sources that will not be completely cancelled by differencing are the code and phase multipath and the receiver noise. Code multipath and noise depend on the structure of the code and the design or the receiver. Phase noise can be assumed to be a few percent of a carrier cycle for all of the signals and similarly phase multipath can be shown to have a maximum amplitude of one quarter of the wavelength and typical values of less than a few centimetres [36].

3.1.1. Frequency Combinations

In all of the cascading methods described earlier only a single code measurement is used for each satellite. In principle, this should be the code measurement with the smallest measurement noise and the smallest multipath. The same can be said for the LAMBDA-based methods. In principle, it is possible to use multiple code observations in a geometry-based least-squares or filtered LAMBDA solution, but in practice, little is gained from this, particularly in short baselines where there are no residual ionospheric errors.

In the geometry-based LAMBDA and cascading simulations presented in what follows, a single code observation from each satellite is used in conjunction with one, two or three phase measurements from each satellite. The following three frequency combinations are assessed.

Single-Frequency
In this case, only L1 and E1/E2 code and phase measurements are used. LAMBDA is applied to decorrelate the estimated ambiguities and obtain PCF estimates. This is contrasted with what could be called “single-frequency cascading” or simply single-frequency ambiguity resolution without the application of LAMBDA. In both cases, the total number of double-difference ambiguities estimated is where and are the number of GPS and Galileo satellites tracked, respectively.

Dual-Frequency
In the dual-frequency case, L2 and E5b phase observations are added to the single-frequency case described earlier. In both the LAMBDA and cascading scenarios this doubles the number of ambiguities to be estimated to . With LAMBDA, the ambiguities consist of some linear combination of all the ambiguities on each frequency (as determined by the algorithm) while with cascading the ambiguities are specifically the Widelane ambiguities (which are resolved first) followed by the L1 and E1/E2 ambiguities.

Triple-Frequency
Finally, for the triple-frequency case, the number of ambiguities increases again to . Again, for LAMBDA, the ambiguities being estimated are some linear combination of the ambiguities of the three frequencies while for triple-frequency cascading the extra-widelane (EWL) ambiguities are first resolved, followed by the widelane (WL) ambiguities, and finally the L1 and E1/E2 ambiguities.
In both the cascaded and LAMBDA cases, the same input covariance matrix of the ambiguities is used, however, in the LAMBDA cases the LAMBDA algorithm is allowed to determine the optimal decorrelating transformation while in the cascading cases a fixed set of linear combinations (EWL, WL, and L1) are used.

3.1.2. Constellations

According to the Galileo Mission High Level Definition document [37], the space segment comprises a constellation of a total of 30 MEO satellites in 3 orbital planes inclined at 56 degrees at 23616 km altitude. For the simulation, the satellites were assumed to be equally spaced in each plane. The satellites in the second and third planes where advanced by 12 and 24 degrees in mean anomaly with respect to the first plane. The relative orientation of the GALILEO constellation with respect to the GPS constellation has not yet been determined and due to the different inclination and orbital radius of the two systems, will not be constant over time either [38]. As a result, the ascending nodes of the three orbital planes of the simulated GALILEO constellation were arbitrarily assigned right ascensions of 0, 120, and 240 degrees respectively. For fair comparison, the real GPS constellation consisting of 30 satellites as it existed at the start of GPS week 1430 was used. GPS is officially a 24 satellite constellation but has had on the order of 30 satellites for the past several years. At the time of writing, there were 32 GPS satellites but a 30 GPS satellite constellation was chosen because previous studies have shown that small differences in the number of satellites (e.g., between 24 and 32) does not greatly affect the results in terms of positioning accuracy or ambiguity resolution. However, differing numbers of satellites can lead to conclusions about one constellation providing better performance than another and in a vast majority of cases the constellation with more satellites provides better performance [23, 39] . To distinguish the two constellations, the Galileo satellites have been arbitrarily assigned the numbers 36 through 65. A combined constellation is shown in Figure 1.

Figure 1: Combined GPS and GALILEO constellation on September 30, 2007 (GPS PRN: 1~32 (noninclusive), GALILEO PRN: 36~65).

3.1.3. Masking Environment

It is assumed that the satellite-borne antennas are facing in the radial direction, have hemispheric gain patterns and an unobstructed view of the sky. An isotropic mask angle of 5 degrees is assumed to limit use of low elevation signals that are likely to be corrupted by large multipath (from other surfaces on the satellite, e.g.).

3.1.4. User Motion and Orbit Descriptions

The simulated users are two low earth orbiting (LEO) satellites in either highly inclined (near polar) or moderately inclined orbits. Two typical LEO users were modelled after two existing LEO satellites. ENVISAT, a European space agency earth observation satellite was selected as a model for a polar orbiting LEO. The satellite is in a near-polar sun-synchronous orbit with an orbital period of 101-minutes and an altitude of approximately 790 km. The International Space Station (ISS) was selected as an example of a LEO in a moderately inclined orbit. The Keplerian elements used for each satellite and some other information are given in Table 3. It should be noted that these orbital elements are approximate values taken from public sources and are neither synchronized with each other nor with the simulated GNSS constellations. They are meant to be representative of each type of user on an arbitrary day and time. Typical ground tracks for the two satellites as well as the ground tracks for the GPS and Galileo satellites are shown in Figures 2 and 3. Note that the GNSS satellites are in medium earth orbit, far above the user LEOs. In the two figures, the blue portions of the GNSS satellite ground tracks indicate periods when they are in the field of view of the LEO and are being used in the solution.

Table 3: Keplerian elements and other orbital information used to simulate LEO users.

Figure 2: Typical ground track for ENVISAT shown in magenta with the Galileo constellation. The circle containing the black cross indicates the initial ENVISAT location (S, W). The GNSS satellites are labelled by PRN number. The blue portion of their ground tracks indicates the GNSS satellite is being used in the solution.

Figure 3: Typical ground track for the ISS shown in magenta with the combined GPS + Galileo constellation. The circle containing the black cross indicates the initial ISS location (N, E). The GNSS satellites are labelled by PRN number. The blue portion of their ground tracks indicates the GNSS satellite is being used in the solution.

3.1.5. Simulation Time and Data Rates

The trajectories of the two user LEOs and the two GNSS constellations were modelled for an arbitrary 24-hour period. In this 24-hour period, filtered geometry-based LAMBDA and cascading ambiguity resolution schemes were implemented with a reset interval of one-minute such that 1440 intervals were evaluated. This was done to create many samples of the initial convergence phase of the filter, which is of interest in terms of ambiguity resolution performance. A data rate of one observation per ten seconds (0.1 Hz) was assumed for the purposes of evaluating ambiguity resolution methods.

At the beginning of each one-minute interval, the Kalman filters used in each method were reset and re-dimensioned to handle the number of satellites that were above the elevation mask at the start of the interval and were above the elevation mask at the end of the interval. This was done to avoid having to introduce or remove ambiguity states part way though the one-minute filtering interval.

4. Data Analysis and Results

Results from using Galileo alone are now presented and discussed followed by a section where the same simulations are repeated but with GPS and Galileo being used together. Each section begins with satellite availability and dilution of precision results, followed by 1, 2, and 3 frequency geometry-based cascading probability of correct fix results, followed by 1, 2, and 3 frequency geometry-based LAMBDA probability of correct fix results. The first section concludes with a discussion of fixed ambiguity baseline solution estimates and corresponding attitude angle error estimates and these results can be generalized to the dual GNSS case.

4.1. Galileo Results

4.1.1. Satellite Availability and Dilution of Precision

Figure 4 shows a plot of the number of available satellites and the corresponding HDOP and VDOP values as a function of time for ENVISAT. Results for the ISS were very similar, and have been omitted to save space. The main difference between the two is that the VDOP of the polar orbiting ENVISAT tests tends to be worse as the satellite travels over the poles due to the lack of overhead GNSS satellites in the polar regions.

Figure 4: Availability, HDOP, and VDOP of Galileo as viewed from a polar orbiting LEO (ENVISAT) over a 24 hour period.

As expected, VDOP values are larger than HDOP. This is due the geometrical distribution of satellites in GNSS and the fact that observability of the vertical direction is limited by the correlation between the vertical position state and the receiver clock offset. This effect remains despite the elimination of the clock offset term in the double-differencing process. A major result of this is that generally the vertical baseline component will be poorer than the horizontal components, making pitch and roll more difficult to estimate than azimuth when receiver antennas are configured in the horizontal plane.

The DOP values are also periodic with the orbital period of the LEO as can be seen by comparing the DOP time series to the latitude of the LEO shown in Figure 5. This is especially true for polar satellites since when they are over the poles, they will only see GNSS satellites on the horizon around them (since GNSS satellites are generally in 54 to 56 degree inclined orbits). As a result, baseline solutions and corresponding attitude estimation quality will vary periodically with the LEO orbit period.

Figure 5: Latitude of ENVISAT as a function of time.

4.1.2. Cascading with Galileo

Results of the geometry-based cascading schemes are now presented. A single filtered 1-minute interval (of the 1440 1-minute intervals simulated) consisting of 6 simulated observations occurring at time = 0, 10, 20, 30, 40, and 50 seconds is presented in detail first. The results for all 1440 trials are then shown together. An arbitrary 1-minute segment for the ENVISAT was selected. For this, and every 1-minute segment, the following procedure is used. First, the assumption is made that the system is warm started meaning that the receiver has approximate coordinates for itself and has acquired and is tracking all visible satellite before the start of the 1-minute segment. At the beginning of the 1-minute segment, the “rover” receiver in the two receiver pair estimates its position with respect to the “base” receiver to an accuracy of 1 m (1σ) in each dimension in the satellite body frame. This initial value is based on typical differential code positioning accuracy. A Kalman filter is then initialized with this initial estimate for the antenna position, and an initial ambiguity estimate variance that corresponds also to 1 m (but is expressed in terms of cycles).

The combined position and float ambiguity filter is then updated once every 10 seconds with a single code observation from each satellite, and a phase observation for each frequency being used from each satellite for the particular scenario. E1/E2 phase observations are used to estimate both an E1/E2 ambiguity and WL, E5b is used to estimate the WL and the EWL, and E5a contributes only the EWL. Due to the short baseline length and fixed nature of the antennas, the position (baseline component) states are modelled as a random walk processes with a process noise of 0.01 /s (about a mean value determined from an existing orbital and attitude model for the satellite) and ambiguity states are modelled as random constant processes. After each Kalman update, the covariance matrix of the ambiguities is analysed to determine the lower bound on the probability of correct fix of the full set of ambiguities using the technique described in Section 2.2.

The Kalman covariance update equation is given by

where is the design matrix, is the measurement error covariance matrix, and and are the covariance matrices of the errors of the states before and after the update respectively. Note that for simplicity the Kalman gain matrix is included in, but not shown as a separate quantity in (14). The prediction step assumes that the baseline components, and thus also the direction cosines in are expressed in terms of the local level frame of the satellite, thus the transition matrix of the system is an identity matrix and is not shown. The Kalman covariance prediction equation is simply

where contains process noise for the baseline component states and is zeros in all rows and columns corresponding the ambiguity states.

For the cascading technique, the probability of partially fixing each step (the EWL alone, the EWL and WL, and all three) is also evaluated. As a final step, the estimated covariance of the fixed position states is computed using (7). The results of this can then be used to compute the estimated accuracy of the derived attitude angles.

Figure 6 shows the overall probability of incorrect fix upper bounds for single-, dual-, and triple-frequency cascading (over one 1-minute interval). To interpret this figure, consider the red line near the top. This is the probability of incorrectly resolving the E1/E2 ambiguity as a function of time in the case that only a single-frequency receiver is used. The line decreases as a function of time indicating increasing odds of correct ambiguity resolution. However, the dual- and triple-frequency cascading options, shown by the green and blue lines, respectively, show that dual and triple-frequency cascading methods offer several orders of magnitude improvement. The -axis is logarithmically scaled and for example indicates a probability of correct fix of . Note that these are the probabilities of fixing all of the ambiguities, and in dual- and triple-frequency cases there are two and three times as many ambiguities to fix. But even though there are more ambiguities to fix, the ability to cascade substantially increases the chances of correctly resolving the ambiguities.

Figure 6: Sample Probability of Incorrect Fix upper bounds for single-frequency, dual and triple-frequency cascading for ENVISAT.

Figure 7 shows the probability of incorrect fix for each stage of the cascading process for the triple-frequency case. The blue line represents PIF for the EWL ambiguities as a function of time. Note that resolving the EWL ambiguities is more or less assured. Though it is not visible, there is a green line plotted almost directly under the red line. This green line represents the probability of correctly resolving the WL after resolving the EWL. The reason it is not visible is that resolving the WL more or less ensures resolving the E1/E2 ambiguity (in the short baseline for attitude determination case) meaning that the probability of resolving the E1/E2 ambiguity after resolving with WL is close to unity. The result is that the probability of resolving all three is only very slightly less than the probability of resolving the EWL and the WL alone. Note that the red line, representing the probability of resolving the E1/E2 ambiguities after resolving the EWL and WL is identical to the probability of resolving all of the ambiguities in triple-frequency cascading which is plotted as the blue line in Figure 6.

Figure 7: Sample Probability of Incorrect Fix upper bounds for the EWL, WL, and E1/E2 ambiguities in triple-frequency cascading for ENVISAT.

Figure 8 shows the corresponding float solution baseline component accuracies as a function of time for the three scenarios shown in Figure 6. The three baseline component labels, E, N, and h, correspond to the longitude, latitude, and vertical direction in the local level frame of the spacecraft. Note that the scale of the vertical subplot is double those of the two horizontal components corresponding the decreased accuracy provided by GNSS in the vertical direction. Note how float solution is relatively poor, on the order of 10 to 20 cm. In this regime, attitude determination on a short baseline would not be feasible.

Figure 8: Sample estimated float solution accuracy for single, dual and triple-frequency cascading. E, N, h indicate longitude, latitude and vertical components in the local level frame of the spacecraft.

Figure 9 shows the corresponding baseline accuracy components for the fixed solution. Note that the accuracy is more or less constant as a function of time. These estimated accuracies are the accuracies assuming the ambiguities have been correctly resolved and this figure does not take into account the fact that at the beginning of the filtering interval the likelihood, for example, of correctly resolving the E1/E2 ambiguity with a single-frequency receiver is quite low. However, this accuracy level does represent the steady state assuming the ambiguities have been correctly resolved and is useful as the attitude accuracy can be derived directly from this baseline accuracy. For small errors, such as those shown in Figure 9, and a simple antenna configuration with two or three antennas in the horizontal plane, the corresponding attitude errors expressed in radians are well approximated by the ratio of the horizontal positioning accuracy to baseline length for azimuth, and by the ratio of the vertical error to baseline length for the roll and pitch errors.

Figure 9: Sample estimated fixed solution accuracy for single, dual, and triple-frequency cascading. , , indicate longitude, latitude, and vertical components in the local level frame of the spacecraft.

Figure 10 shows the PIF as a function of time for single, dual and triple-frequency cascading for ENVISAT. All 1440 simulations are plotted on one figure to demonstrate the range of PIF for varying LEO to Galileo satellite geometry. As can be seen, the sample result shown in Figure 6 was in fact one of the best cases and even with three frequencies, the PIF can still be relatively high (e.g., only ) after a minute of filtering. Also, the larger performance gain is from one frequency to two, with all of the PIF samples for the single-frequency case being worse than any of the dual or triple-frequency samples whereas there is some overlap between good dual-frequency results and relatively poor triple-frequency results. Figure 11 shows all of the samples of EWL, WL and E1/E2 PIF for the triple-frequency case as discussed earlier in the context of Figure 7. Again, the sample result shown earlier is one of the better samples. From Figure 11 it can be seen that the probability of resolving the EWL in a triple-frequency cascading scheme is very high and very consistent between samples while there is more variability between the probabilities of resolving the WL in the second step in the cascading scheme. It can also be seen that there is little difference in the probability of resolving the WL and the probability of resolving all three, indicating again that once the widelanes are resolved in triple-frequency cascading, resolving the E1/E2 ambiguities is relatively easy.

Figure 10: Probability of Incorrect fix upper bounds for single-frequency, dual and triple-frequency cascading for ENVISAT.

Figure 11: Probability of Incorrect Fix upper bounds for the EWL, WL, and E1/E2 ambiguities in triple-frequency cascading for ENVISAT.

4.1.3. Lambda with Galileo

With the LAMBDA method, as opposed to cascading, all of the ambiguity states are estimated together and ambiguity resolution is facilitated by the LAMBDA algorithm deciding what linear combination of the original ambiguities most decorrelates the ambiguities thus making them easiest to fix. In terms of the covariance simulation, the same Kalman filter equations are used, only in this case the design matrix reflects the fact that the original ambiguities (and not widelanes) are being estimated. The result of the float solution is then sent to the LAMBDA algorithm for decorrelation.

Figure 12 shows the PIF results for the LAMBDA method and can be compared directly the corresponding cascading results shown in Figure 10. Note that there is improvement in all three cases (single-, dual-, and triple-frequency) but most notably the triple-frequency ambiguity resolution PIF (the blue lines) is significantly reduced compared to cascading demonstrating the potential of the LAMBDA method to find optimal linear combinations of the ambiguities that provide more ambiguity resolution potential compared the EWL, WL, and E1/E2 used in the cascading method.

Figure 12: Probability of Incorrect Fix upper bounds for single-, dual- and triple-frequency LAMBDA for ENVISAT.

One important result to note is that less is gained by adding a third frequency in terms of ambiguity resolution for very short baselines. In other words, there is a larger change in PIF from the single-frequency case to the dual-frequency case than there is from the dual- to triple-frequency case. The addition of the second frequency allows the LAMBDA algorithm to automatically form a widelane combination, which for a short baseline such as this, is more than sufficient for ambiguity resolution. The addition of the third frequency provides only a marginal improvement as the code error in this case is smaller that the extra-wide lane wavelength and there is no differential ionosphere error affecting either the code or phase measurements in this application. This is similar to the results of using cascading (found in the TCAR/CIR literature) where the largest gains from using three frequencies are found in longer baseline applications. A final note about the LAMBDA results is that close inspection shows a small number of dual and triple-frequency trials where the PIF appears to increase counter-intuitively from one epoch to the next. This is not an error. An unfortunate feature using the bootstrapped lower bound for the PCF as an estimate of the integer least squares PCF is that the lower bound is not invariant to the decorrelating linear combination found using LAMBDA. Occasionally, the LAMBDA algorithm changes its linear combination from one filter epoch to the next and this can result in discontinuities, both decreasing and increasing, in the PCF lower bound. It should be noted that the actual PCF does not change, only value of the bound.

4.2. GPS/Galileo Results

In this section, the results presented for Galileo are now repeated for the case that both GPS and Galileo are being used.

4.2.1. Availability and Dilution of Precision

With two systems, the number of satellites observed roughly doubles, as shown in Figure 13. The corresponding DOP values also decrease. As in the Galileo only case, the largest (poorest) DOP values occur when the fewest satellites are tracked and the results are periodic with the orbital period of the LEO satellite.

Figure 13: Availability, HDOP, and VDOP of GPS/Galileo as viewed from a polar orbiting LEO (ENVISAT) over a 24-hour period.

4.2.2. Cascading with GPS/Galileo

Figure 14 shows the dual GNSS probabilities of incorrect fix bounds for cascading for ENVISAT analogous to those presented in Figure 10 for the Galileo only case. The same patterns are observed, though with roughly double the number of observations and double the number of ambiguities to estimate, the performance improvement is only marginal. Likewise, the probability of successfully resolving each stage in the three-frequency case, shown in Figure 15, is only marginally better than the results for Galileo only shown in Figure 11.

Figure 14: Probability of Incorrect Fix upper bounds for single-frequency, dual and triple-frequency cascading for ENVISAT using two GNSS.

Figure 15: Probability of Incorrect Fix upper bounds for the EWL, WL, and E1/E2 ambiguities in triple-frequency cascading for ENVISAT using two GNSS.

4.2.3. Lambda with GPS/Galileo

Use of the LAMBDA method with two systems results in a major improvement in the probability of correct fix than can be obtained with a single-frequency receiver. Though the dual- and triple-frequency cases are also improved with LAMBDA, the significant improvement can be seen in the red lines in Figure 16 when compared to the Galileo only case shown in Figure 12. Most striking is the absence of cases where the probability of incorrect fix hovers near the top of the figures (in the 0.1 to 1.0 range). This can be explained by the fact that the LAMBDA method is able to exploit the geometric diversity of the single-frequency signals from the two systems by forming linear combinations of all of the L1 and E1/E2 ambiguities while in the Galileo only case, approximately half as many ambiguities are available. There is less improvement with the dual and triple-frequency cases because these cases already have a diversity of ambiguities on different frequencies for the LAMBDA algorithm work with even in the Galileo only case.

Figure 16: Probability of Incorrect Fix upper bounds for single-frequency, dual and triple-frequency LAMBDA for ENVISAT using two GNSS.

4.3. Galileo versus GPS/Galileo: Recommendations

Based on these results it can be concluded that the LAMBDA approach is superior to the cascading approach in all cases. However the LAMBDA method is particularly useful in the case of a single-frequency dual-GNSS application. However to obtain the most reliable ambiguity resolution, multiple frequencies are required when using both one or two GNSS constellations. If power and weight of the GNSS receiver payload is no object, a dual-system triple-frequency receiver using LAMBDA will provide the best performance. However, if power and weight are an issue, a dual-frequency single-system receiver offers similar performance to a single-frequency dual-system one. Similar results were obtained with the ISS as the simulated user and are therefore not shown.

While these results consider the case of attitude determination for an orbiting user, they can be generalized to ground based users. Of course while the absolute performance of each algorithm will depend on measurement accuracies, satellite geometry, and user dynamics, it is reasonable to expect the relative performance of the algorithms to remain unchanged.

5. Conclusions

In this paper, the effectiveness of two approaches to multiple-frequency carrier phase ambiguity resolution was evaluated for case of attitude determination onboard a low-earth orbiting satellite. The evaluation is based on simulations of a polar orbiting LEO and an inclined orbiting LEO. The overall conclusion is that the geometry-based LAMBDA method is superior to the geometry-based cascading method in terms of probability of correct ambiguity resolution, and time required to achieve a particular probability of correct fix. Further conclusions may also be made about the various frequencies and systems available. The use of two GNSS provides a significant increase in ambiguity resolution performance while the addition of the third frequency on each provides only marginal improvement compared to the dual-frequency case.