BeiDou Inter-Satellite-Type Bias Evaluation and Calibration for Mixed Receiver Attitude Determination

The Chinese BeiDou system (BDS), having different types of satellites, is an important addition to the ever growing system of Global Navigation Satellite Systems (GNSS). It consists of Geostationary Earth Orbit (GEO) satellites, Inclined Geosynchronous Satellite Orbit (IGSO) satellites and Medium Earth Orbit (MEO) satellites. This paper investigates the receiver-dependent bias between these satellite types, for which we coined the name “inter-satellite-type bias” (ISTB), and its impact on mixed receiver attitude determination. Assuming different receiver types may have different delays/biases for different satellite types, we model the differential ISTBs among three BeiDou satellite types and investigate their existence and their impact on mixed receiver attitude determination. Our analyses using the real data sets from Curtin's GNSS array consisting of different types of BeiDou enabled receivers and series of zero-baseline experiments with BeiDou-enabled receivers reveal the existence of non-zero ISTBs between different BeiDou satellite types. We then analyse the impact of these biases on BeiDou-only attitude determination using the constrained (C-)LAMBDA method, which exploits the knowledge of baseline length. Results demonstrate that these biases could seriously affect the integer ambiguity resolution for attitude determination using mixed receiver types and that a priori correction of these biases will dramatically improve the success rate.

In this paper, we consider mixed receiver attitude determination using single-and dual-frequency BeiDou observables. Assuming different receiver types may have different delays/biases for different satellite types, we model the differential ISTBs among three BeiDou satellite types. We develop an extended GNSS model, taking into account these biases, and describe the estimation of these biases. Our analyses using the data from two real data campaigns (one is with Curtin's permanent GNSS array consisting of different types of BeiDou-enabled receivers, and the other is with zero-baseline using BeiDou enabled receivers) reveal the existence of non-zero ISTBs between different BeiDou satellite types. We observe that these biases are stable and constant. Hence, we use a priori estimated biases to correct BeiDou observations, so that one can use classical double difference processing without loosing redundancy. We then analyse the impact of these biases on BeiDou-only attitude determination using the constrained (C-)LAMBDA method, which exploits the knowledge of baseline length. Results demonstrate that ISTBs could seriously affect the integer ambiguity resolution for attitude determination using mixed receiver types and that a priori correction of these biases will dramatically improve the success rate.
This contribution is organized as follows. Section 2 describes the functional and stochastic model for BeiDou observations, with special attention to ISTBs and associated processing approaches. Section 3 formulates the quadratically-constrained BeiDou model for attitude determination and describes the C-LAMBDA method. Section 4 demonstrates the existence of non-zero ISTBs between different BeiDou satellite types using real data and presents the results of attitude determination, revealing the impact of ISTBs. Finally, Section 5 contains the summary and conclusions of this contribution.

BeiDou Observations
This section presents the BeiDou observation model, distinguishing satellite types, namely GEO, IGSO and MEO, to accommodate receiver-dependent delays (biases) for different satellite types. Since the BeiDou system uses the code division multiple access (CDMA) technique, the code and phase observations of receiver r tracking satellite s β of type β at frequency j, denoted by p s β r,j and φ s β r,j , respectively, are given as [50]: with f the number of frequencies, m β the number of type-β satellites tracked, η the number of satellite types, ρ s β r the topocentric distance between receiver r and satellite s β , dt r the receiver clock error, dt s β the satellite clock error, a s β r,j the (frequency-dependent) code atmospheric delay, α s β r,j the (frequency-dependent) phase atmospheric delay, d ,β r,j the (satellite-type-dependent) hardware code delay in the receiver, d s β ,j the hardware code delay in the satellite, λ j the wavelength, δ ,β r,j the (satellite-type-dependent) hardware phase delay in the receiver, ϕ r,j the initial phase in the receiver, δ s β ,j the hardware phase delay in the satellite, ϕ s β ,j the initial phase in the satellite, N s β r,j the (integer) phase ambiguities, e s β r,j all other code errors, including measurement noise, and s β r,j all other phase errors, including measurement noise. For simplicity of formulation, we assume that satellites are ordered, such that the first m 1 satellites are of type 1, the next m 2 satellites are of type 2, and so on ( η β=1 m β = m, the number of tracked satellites). Note that all variables are expressed in meters; except the satellite-type-dependent biases of phase observations and the ambiguity, which are expressed in cycles.
The between-receiver single difference (SD) pseudo-range and carrier-phase observations of two BeiDou receivers, r and 1, at frequency j from satellite s β of type β, denoted as p s β 1r and φ s β 1r , respectively, are given as: where the satellite-specific biases are eliminated. Attitude determination for a small platform considered in this paper is based on multiple BeiDou receivers forming very short-baselines for which the differential atmospheric delays can be ignored, that is, a s β 1r,j = α s β 1r,j = 0 ∀j, s β . Further differencing the above SD observables between satellites eliminates the receiver-dependent biases: with the first satellite of the satellite-type 1 as the pivot, the double difference (DD) observables are given as: In classical double differencing for single satellite-type systems, such as GPS, the terms, d ,1β 1r,j and δ ,1β 1r,j , do not exist. However, as shown in Section 4.2, there may exist non-zero differential inter-satellite-type biases (ISTBs) for the BeiDou system if one uses mixed receiver types. That is, δ ,1β 1r,j = 0 and d ,1β 1r,j = 0 for β = 1.
The linearized DD observation equations corresponding to Equations (5) and (6) read: E(∆φ where E(·) denote the expectation operator, ∆p 1 1 s β 1r,j and ∆φ 1 1 s β 1r,j are the observed-minus-computed code and phase observations, b is the baseline vector containing relative position components and g 1 1 s β r is the geometry vector given as g r the unit line-of-sight vector from receiver r to satellite s β .

Classical DD Model
The naive approach is to simply ignore the third term in Equation (18), resulting in the classical baseline model with a full-rank system: Hence, the redundancy of this classical model is equal to As demonstrated in Section 4.3, ignoring the ISTBs results in catastrophic failure of ambiguity resolution.

Extended DD Model with Estimable ISTBs
The rank deficiency in Equation (18) can be removed by constraining a set of parameters (combinations) as the S-basis [53][54][55]. The number of S-basis constraints equals the size of the rank deficiency. There are many possibilities to choose S-basis constraints. One such choice corresponds to a reparametrization in Equation (6), such that the DD ambiguities of the pivot satellites of second satellite types are combined with corresponding phase ISTBs: 1r,j , is now biased by the inter-satellite-type ambiguity between the pivot satellites. Hence, to avoid the jumps due to the changes of reference satellites and cycle slips, we report the fractional part of the estimable phase ISTBs, which are sufficient for ISTB correction, as discussed in Section 2.4. Instead of lumping phase ISTB with the ambiguity of the first (pivot) satellite, one can lump the phase ISTB with the DD ambiguity of the satellite, other than the first satellite of second satellite type, and end up with different S-bases and different estimable phase ISTBs.
Another option, i.e., another S-basis choice, is to lump the phase ISTB with the average of two or more DD ambiguities of the second satellite type. For example, if we lump the phase ISTB with the average of all DD ambiguities of the second satellite type, the DD phase observations of the second satellite type read: 1r,j + with estimable integer ambiguities: and estimable phase ISTB:δ ,1β For this choice of S-basis, however, the fractional part of the estimable phase ISTB is not necessarily equal to that of the actual phase ISTB. Having different choices, one should be cautious while interpreting or using the estimated phase ISTBs based on a given S-basis choice.
In this contribution, we consider the parametrization of Equations (20) and (21), which has the simple interpretation that the estimable integer ambiguities correspond to satellite type-specific double differencing (see Section 2.3), and moreover, it enables us to observe the fractional parts of the phase ISTBs (see Section 4.2). For our choice of S-basis, the full-rank (extended) model reads: . . , C mη , C n = [0 I n−1 ] T and blockdiag referring to the block diagonal matrix formed by given arguments,z = z T ,1 , . . . ,z T ,f

Type-Specific DD Model
Since ISTBs are nuisance parameters, one can eliminate 2f (η − 1) ISTBs using the differencing operator, is the differencing matrix. Pre-multiplying Equation (26) byD T nullifies the third term and results in a typespecific DD model, which has a reference satellite (pivot) per satellite type and reads: This model has 2f (η − 1) less observations and 2f (η − 1) less unknown parameters than Equation (26). Hence, both models have the same redundancy and are equivalent.

ISTB-Corrected DD Model
As shown in Section 4.2, ISTBs are stable and can be assumed to be constant for a given receiver type pair. Hence, one can correct DD observations in Equation (18) with a priori estimates of ISTBs (calibration): withỹ = y − Hh andh the vector of a priori ISTBs known through calibration and consisting of code and fractional phase ISTBs (see Section 4.2). Note that it is sufficient to use the fractional part of phase ISTB to correct phase observations, as the integer part of the phase ISTB can be lumped with the corresponding integer ambiguities without affecting integer ambiguity resolution. The redundancy of an ISTB-corrected system is equal to f (m − 1) − 3. Hence, this model is stronger than the extended model in Equation (26). Finally, Table 1 summarizes the full-rank models described in the above. For further analyses in this paper, we only consider three models, namely, the classical DD model, ignoring ISTBs, the extended model and the classical DD model with ISTB correction. The type-specific DD model is equivalent to the extended DD model. Table 1. Redundancies of models considered. DD, double differencing; ISTB, inter-satellitetype bias.

Model Redundancy
Classical DD with ignoring ISTBs Equation (19)

Stochastic Model
We assume elevation-dependent noise characteristics [56] for the undifferenced observables in Equations (1) and (2). That is, the standard deviation of the undifferenced observable, ς, can be written as: where θ is the elevation angle of the corresponding satellite and σ ς 0 , a ς 0 and θ ς 0 are the elevation-dependent model parameters. We further assume that the receivers have similar characteristics and that the observation noise standard deviations can be decomposed as follows: where σ r and σ ,j are the receiver, the frequency and the satellite-type-dependent weightings, respectively, and σ φ 0 and σ p 0 are observation-dependent weightings.
To construct the stochastic model for the observations in Equation (18), consider the undifferenced observations reading:  (1) and (2). Using the noise characteristics of Equation (30) and assuming that the observables are normally distributed and mutually uncorrelated, the dispersion matrix of the observation vector, ζ, can be written as: Using the DD operator, the dispersion matrix of the DD observations is then given as: For type-specific DD observations in Equation (27), the dispersion matrix is given as:

BeiDou Attitude Determination
Since attitude determination uses rigidly mounted antennas, the baseline length is a priori known and can be used to strengthen the underlying model. With the inclusion of the baseline length constraint to the models in Equations (19), (26) or (28) and with the stochastic model in Equation (34), we obtain the GNSS compass model [42,47]: where l is the known baseline length, · denotes the unweighted norm and κ is the number of integer ambiguities. The parameters for different models are defined as follows: Classical model Equation (19): Extended model Equation (26): ISTBcorrected classical model Equation (28): In the above, the baseline is constrained to lie on a sphere with radius l (S l = {b ∈ R 3 | b = l}). Our objective is to solve for b in a least-squares sense, thereby taking the integer constraints on z and the quadratic constraint on vector b into account. Hence, the least-squares minimization problem that will be solved reads: It is a quadratically-constrained (mixed) integer least-squares (QC-ILS) problem [41,47], for which no closed-form solution is available. In the following sections, we describe the method for solving Equation (39).

The Ambiguity Resolved Attitude
We now describe the steps for computing the integer ambiguity resolved attitude angles.

The Real-Valued Float Solution
The float solution is defined as the solution of Equation (39) without the constraints. When we ignore the integer constraints on z and the quadratic constraint on b, the float solutions,ẑ,b andĥ, and their variance-covariance matrices are obtained as follows: The z-constrained solution of b and its variance-covariance matrix can be obtained from the float solutions as follows:b The zand b-constrained solution of h and its variance-covariance matrix can be obtained from the float solutions as follows:ĥ Using the above estimates, the original problem in Equation (39) can be decomposed as: withê = y − Aẑ − Gb − Hĥ being the vector of least-squares residuals. Note that the first term on the right-hand side of Equation (46) does not depend on the unknown parameters, z, b and h, and is therefore constant. Unlike the third term, which is constrained, the last term can be reduced to zero for any given z and b.

The Integer Ambiguity Resolution
Based on the orthogonal decomposition Equation (46), the quadratic-constrained integer minimization can be formulated as:ž with ambiguity objective function: where:b The cost function, C(z), is the sum of two coupled terms: the first weighs the distance from the float ambiguity vector,ẑ, to the nearest integer vector, z, in the metric of Qẑẑ, while the second weighs the distance from the conditional float solution,b(z), to the nearest point on the sphere, S l , in the metric of Qb (z)b(z) . Unlike with the standard LAMBDA method [29], the search space of the above minimization problem is non-ellipsoidal, due to the presence of the second term in the ambiguity objective function. Moreover, its solution requires the computation of a nonlinear-constrained least-squares problem (49) for every integer vector in the search space. In the C-LAMBDA method, this problem is mitigated through the use of easy-to-evaluate bounding functions [47]. Using these bounding functions, two strategies, namely the Expansion and the Search and Shrink strategies, were developed; see, e.g., [41,45]. These techniques avoid the computation of Equation (49) for every integer vector in the search space and compute the integer minimizer,ž, in an efficient manner.

The Ambiguity Resolved Parameter Estimation
For a single baseline, b is related to the Euler-angles, ξ = [φ, θ] T , with φ the heading and θ the elevation, as b(ξ) = lu(ξ), where u(ξ) = [c θ c φ , c θ s φ , −s θ ] T , with s α = sin(α) and c α = cos(α). The sought for attitude angles, ξ (ž), are the reparametrized solution of Equation (49). Using a first order approximation, the formal variance-covariance matrix of the ambiguity resolved, estimated heading and elevation angles are given by: with Jacobian matrix: Finally, in the case of the extended model, ambiguity-corrected ISTB estimates,ĥ(ž,b(ž)), and the associated variance-covariance matrix, Qĥ (z,b)ĥ(z,b) , are obtained using Equations (44) and (45). The computation of the fractional phase ISTBs from these estimable ISTBs is discussed in Section 4.2.

Measurement Campaign
The analyses in this paper are based on data sets from Curtin University's permanent GNSS stations and series of zero-baseline experiments. The first data set is from Curtin's permanent GNSS antennas (CUT00 and CUTA0) mounted on the roof of building 402 at the campus of Curtin University in Perth, Australia (Figure 1a). These antennas are connected to six BeiDou-enabled receivers, as summarised in Figure 1b, consisting of three Trimble NETR9, two Javad TRE G3T DELTA and a Septentrio POLARx4 receivers. We considered the data from these receivers for five days from April 4 to 8, 2013. BeiDou satellite visibility for April 4 is shown in Figure 2. The data of various zero baselines for five days with a 30 s sampling interval is used to estimate ISTBs in Section 4.2, while the data of various short baselines between CUT00 and CUTA0 on April 7 with a 1 s sampling interval is used to analyse the impact of ISTBs on attitude determination in Section 4.3.   (Table 2). As shown in Figure 3b, a single antenna (Figure 3a) was connected to two BeiDou-enabled receivers (Trimble NETR9 and Javad Javad TRE G3T DELTA). Figure 4 shows the visibility of BeiDou satellites on April 19, 2013. These zero-baseline data sets (with a 30 s sampling interval) are also used to estimate and validate ISTBs in Section 4.2. The stochastic model parameters of the elevation-dependent model Equation (29) for the receivers are reported in Table 3. Since the receivers, except Trimble NETR9, track only B1 and B2 signals, single-and dual-frequency analyses are considered in the paper.

BeiDou Inter-Satellite-Type Bias (ISTB)
First, we considered the estimation of differential ISTBs using zero baseline data. Since the geometry term vanishes for a zero baseline problem, the estimation of code and phase ISTBs for each frequency are decoupled. Using the extended model in Equation (26) and the associated stochastic model in Equation (34), the decoupled system for the differential code ISTBs at the jth frequency is given as: The estimates of differential code ISTBs are then given by the least-squares solution of the above system. Similarly, using the extended model in Equation (26) and the associated stochastic model in Equation (34), the decoupled system for the differential estimable phase ISTBs at the jth frequency is given as: withL ,j = λ j blockdiag I m 1 −1 , C m 2 , . . . , C mη . First, the float solution of the above full-rank square system is obtained by ignoring integer constraints: Since the above system is driven by phase measurement noise, simple rounding ofẑ ,j yields integer ambiguities,ž ,j . The estimates for the estimable phase ISTBs are then given by: Since these estimable phase ISTBs are the sum of actual phase ISTBs and corresponding ambiguities of reference satellites of the second satellite type, the estimates are affected by integer jumps, due to the cycle slips and the changes of reference satellites. Hence, we report only fractional phase ISTBs (the residuals of integer rounding) that are sufficient for ISTB correction as discussed in Section 2.4. That is, the fractional phase ISTBs,ĥ φ;j (ž ,j ) =ĥ φ;j (ž ,j ) − round ĥ φ;j (ž ,j ) , where 'round' refers to the closest integer to the given estimate. However, these fractional phase ISTB estimates are ambiguous if they are equal to or close to a half cycle (e.g., for a half cycle, simple rounding will yield either +0.5 − or −0.5 + , depending on the noise, ). For this situation, we use either "floor" (resulting in the residual for the nearest integer that is smaller than the given estimate) or "ceiling" (resulting in the residual for the nearest integer that is larger than the given estimate) with the sign convention in Table 4. Note that one is free to choose any sign convention, as long as it preserves the consistency of the double difference ambiguities when the reference receiver and/or the reference satellite type changes.  It is observed that the estimated ISTBs (monitored for several days) are very stable and can be used to calibrate BeiDou observations. In the following, we compute the ISTB corrections using epoch-by-epoch estimates of several days. Let us consider the time series of the ith ITSB and associated standard deviations, {h i:k , σ h i:k } K k=1 , where K is the number of epochs and h i:k is the estimated code or phase ISTB at time k. Assuming these estimates are uncorrelated over time, we formulate the following least-squares problem for the bias estimation: with h i = [h i:1 , . . . , h i:K ] T the K × 1 vector of epoch-by-epoch ISTB estimates for the ith ISTB. The least-squares estimate (weighted mean) of the bias and its standard deviation are given as:  (f)    (f)  (f)  (f) Table 5 and 6 summarize the ISTB estimates and their standard deviations, clearly indicating the existence of non-zero ISTBs (highlighted using bold text) between dissimilar receiver types. It was observed that the estimated ISTBs are constant for given receiver-antenna connectivity and receiver operating environment. Nevertheless, it was found that the observed code ISTBs do not significantly affect ambiguity resolution and the consequent ambiguity resolved phase-only baseline estimation. However, the phase ISTBs severely affect ambiguity resolution, especially in the case of half-cycle phase ISTBs. As summarized in Table 7, GEO satellites have phase ISTBs of half-cycles with respect to IGSO/MEO satellites in the case of mixed receivers. Note that phase ISTBs for other receiver pairs can be deduced from estimates in Table 7. For example, with the sign convention in Table 4, the Septentrio-Javad pair has phase ISTBs of -0.5 cycle and 0.5 cycle for B1 and B2, respectively. Note that, for attitude determination with the ISTB-corrected model in the following section, we only correct phase observations with half cycles, as other (code) biases were found to be small enough to not affect the ambiguity resolution significantly.  Table 6. Estimated ISTBs using data from zero baseline experiments with Trimble-Javad receiver pairs: two in Curtin University and another two in Kalamunda (Table 2).  Table 7. Differential phase ISTB between GEO and IGSO/MEO satellites with Trimble as the pivot receiver (cyc) and based on the sign convention considered in Table 4 Frequency Trimble Septentrio Javad

Effect of ISTBs on Attitude Determination
Next, we analysed the impact of ISTBs on BeiDou single-and dual-frequency instantaneous attitude determination using the standard LAMBDA and C-LAMBDA methods comparing three processing approaches, namely, the classical DD model with ignoring ISTBs in Equation (19), the extended model in Equation (26) and the classical DD model with ISTB correction in Equation (28). Results for different receiver pairs (CUT0-CUTA, CUT0-CUAA, CUT1-CUTA and CUT1-CUAA), which consist of mixed receivers forming a short baseline of 8.418 m, as shown in Figure 1b, are discussed in the following. Note that these receiver pairs are not used in the computation of ISTBs in Section 4.2. We considered two performance measures for our analyses; the first one is the empirical instantaneous ambiguity success fraction (relative frequency), which is defined as: success fraction = number of correctly fixed epochs total number of epochs (64) where the true ambiguities are computed using known antenna coordinates in WGS84, as the antennas used are part of Curtin's permanent stations. However, only length information is used for C-LAMBDA processing. The second performance measure is the ambiguity fixed angular estimation accuracy, which is given by the formal and empirical standard deviations of attitude angular estimates. Table 8-10 report the instantaneous ambiguity success fraction for single-frequency B1, single-frequency B2 and dual frequency B1-B2 processing, respectively. The first row in each table corresponds to the baseline with the same receiver type for which ISTBs are zero and correction is not needed. The third row in Table 8 and the second row in Table 9 correspond to baselines with dissimilar receiver types, which have zero ISTBs for corresponding frequencies. The benefits of using C-LAMBDA, which makes use of known baseline length, are highlighted using bold text. Furthermore, catastrophic failures of ambiguity resolution by ignoring non-zero ISTBs are highlighted with emphasized text. Hence, it is wise to use the extended model (or the type-specific DD model) if one does not have the knowledge of ISTB between dissimilar receiver pairs. However, the best processing strategy is to use ISTB-corrected classical double differencing. With ISTB calibration, the C-LAMBDA method yields single-frequency instantaneous attitude determination.  Table 9. Instantaneous ambiguity success fractions (relative frequencies) for single-frequency (B2) processing.

Baseline
Classical DD model Equation (19) Extended model Equation Finally, Table 11 reports ambiguity fixed angular accuracy for single-and dual-frequency processing. Since the baselines (8.418 m) considered in these analyses are formed by receivers with similar noise characteristics, the ambiguity fixed angular standard deviations are the same for all cases, except the cases with catastrophic failure of ambiguity resolution. Hence, we report the average angular standard deviation of all other cases. Single-frequency processing with either B1 or B2 yields the same angular accuracy. The improved dual-frequency angular accuracy reflects the increased redundancy. Table 10.

Baseline
Classical DD model Equation (19) Extended model Equation  Table 11. Empirical and formal (given in brackets) angular standard deviation (deg).

Summary and Conclusions
In this contribution, we investigated the existence of BeiDou inter-satellite-type biases (ISTBs) and their impact on standalone BeiDou attitude determination with mixed receiver types. We considered an extended GNSS double difference model incorporating all possible differential ISTBs among the three BeiDou satellite types (GEO, IGSO and MEO), together with three processing approaches, namely, one based on the classical double differenced model, ignoring the ISTBs, another based on the extended double differenced model, incorporating the ISTBs, and a third one based on the ISTB-corrected classical double differenced model. Our analyses using two real data sets with three different receiver types demonstrate the existence of non-zero ISTBs between different satellite types. The estimated ISTBs are stable and can be used to correct mixed receiver BeiDou attitude determination. It was observed that the estimated ISTBs are constant for a given receiver-antenna connectivity and receiver operating environment. Nevertheless, it was found that the observed code ISTBs do not significantly affect ambiguity resolution and the consequent ambiguity resolved phase-only baseline estimation. However, the mixed receiver half-cycle phase ISTBs severely affect ambiguity resolution. This finding is an important warning for mixed receiver type users, including precise point positioning users [54,55,57,58]. Moreover, it may also trigger GNSS receiver manufacturers to develop mutually consistent measurement extractions, as they are in the early stage of BeiDou-enabled receiver developments. Furthermore, it is suggested to use the extended model or, equivalently, the type-specific DD model, if one does not have the knowledge of ISTBs between dissimilar receiver pairs. However, the best processing strategy is to use the ISTB-corrected classical double differencing procedure. With ISTB correction, the C-LAMBDA method enables single-frequency, instantaneous attitude determination capability in the Asia-Pacific region with the current BeiDou constellation.