Carrier Phase-Based Precise Heading and Pitch Estimation Using a Low-Cost GNSS Receiver

Abstract

Attitude and heading estimation methods using the global navigation satellite system (GNSS) are generally based on multi-antenna deployment, where the installation space and system cost increase with the increase in the number of antennas. Since the single-antenna receiver is still the major choice of the mass market, we focus on precise and reliable heading and pitch estimation using a low-cost GNSS receiver. Carrier phase observations are precise but have an ambiguity problem. A single difference between consecutive epochs can eliminate ambiguity and reduce the measurement errors. In this work, a measurement model based on the time-differenced carrier phases (TDCPs) is utilized to estimate the precise delta position of the antenna between two consecutive epochs. Then, considering the motion constraint, the heading and pitch angles of a moving land vehicle can be determined by the components of the estimated receiver delta position. A threshold on the length of the delta position is selected to avoid large errors in static periods. To improve the reliability of the algorithm, the Doppler-aided cycle slip detection method is applied to exclude carrier phases with possible cycle slips. A real vehicular dynamic experiment using a low-cost, single-frequency GNSS receiver is conducted to evaluate the proposed algorithm. The experimental results show that the proposed algorithm is capable of providing precise vehicular heading and pitch estimates, with both the root mean square errors being better than 1.5°. This also indicates that the cycle slip exclusion is indispensable to avoid unexpected large errors.

1. Introduction

Precise heading and pitch angles are essential building blocks for many vehicular applications. Therefore, in previous works, numerous efforts were devoted to partial or fully three-dimensional (3D) attitude estimation using different sensors and techniques. However, compared with the attitude estimation methods of inertial sensors [1,2,3,4], cameras [5,6], and multi-sensor integrated systems [7,8], the advantages of the attitude estimation methods based on the global navigation satellite system (GNSS) include being simple, cost effective, easy to use, and immune to accumulated errors [9,10].

Gade summarized seven ways to find the heading under different application situations, where four of them depend on GNSS technology [6]. Gross et al. utilized the signal strength of the global positioning system (GPS) and the derived velocity vector for attitude estimation [10]. Goh and Low conducted a survey of the GPS-based attitude estimation methods [11], and it was shown that the estimation-based methods generally showed better accuracy than the analytical methods. The GNSS compass model relies on a short moving baseline to determine the vehicular heading and pitch angles [12,13,14], where the integer ambiguity resolution (AR) is necessary to improve the estimation accuracy. Multi-antenna GNSS attitude estimation also attracted considerable research attention to providing full 3D attitude information [15,16]. Raskaliyev et al. later contributed a literature review on this subject covering various estimation techniques [17]. The system model with baseline constraints is highly recommended since it can estimate the integer ambiguity and attitude simultaneously and improve the AR rate. However, the constrained model induces the attitude determination process at a greater computational cost [18,19,20]. Aside from that, the deployment of a multi-antenna system is generally less favorable than that of a single antenna due to the extra space requirement and the increased system cost.

The time-differenced carrier phases (TDCPs) of a single-antenna GNSS receiver are in essence the between-epoch single difference (SD) of its carrier phase measurements. With the TDCP technique, the time-correlated errors can be eliminated or largely reduced if the sampling interval is short [21,22]. More importantly, the biggest impediment of the carrier phases, the integer ambiguities, can also be eliminated as long as there is no cycle slip between the two consecutive epochs [23]. Therefore, the TDCP measurements have been widely applied for precise velocity estimation, either for the averaged velocity from the phase-derived position increment (receiver delta position) [24,25,26] or the instantaneous velocity directly estimated from GNSS observations [27]. Aside from that, some recent investigations also applied TDCP measurements for vehicular navigation [28,29,30,31].

Graas and Soloviev presented a meticulous discussion on the adjustment and compensation of the time-differenced carrier phase and the satellite–receiver geometry change [24]. Serrano et al. evaluated the achievable accuracy of TDCP-based velocity estimation [32]. Ding and Wang improved the computation efficiency by avoiding the dependency on satellite velocities and reducing the required epochs of the receiver positions [26]. Freda et al. compared the velocity estimation accuracy of three methods: the direct position difference method, the raw Doppler-based method, and the TDCP-based method [25]. To reduce the influence of the positioning error on the accuracy of the line-of-sight (LOS) unit vector, Duan et al. implemented the precise point positioning (PPP) technique instead of the conventional single-point positioning (SPP) method for Doppler-based velocity estimation [33].

However, there are several intrinsic limitations in the TDCP technique: its vulnerability to cycle slips (CSs); the reduced data availability of the carrier phases compared with that of the code measurements; and the estimated solutions from the TDCP technique being discontinuous since they are calculated on an epoch-by-epoch basis. In this work, the possible cycle slips are detected and excluded using high-rate Doppler measurements, whereas the other two shortcomings can be further addressed by introducing aiding sensors in future work.

The loss of lock indicator (LLI) provided by the receiver is generally referred to for revealing possible cycle slips, but Zhao et al. showed that the LLI is not always reliable [34]. A survey of cycle slip detection methods for single-frequency GNSS receivers was conducted by Farooq et al. [35]. Li et al. applied the positional polynomial constraint for single-frequency GNSS cycle slip estimation [36]. The TDCP technique has also been implemented for the geometry-based cycle slip detection and repair method [37]. This process is akin to ambiguity determination in relative positioning that requires integer estimation and validation. Du and Gao applied the inertial aided cycle slip detection and identification method [38]; however, this method mainly relies on the performance of the inertial sensors. The hypothesis test method was also applied for cycle slip detection by Zangeneh-Nejad et al. [39], but the calculations are more complicated. Theoretically, the variation of the carrier phase between consecutive epochs can be predicted by the Doppler measurements, which are immune to cycle slips. Thus, Doppler-aided cycle slip detection is straightforward for a stand-alone single-frequency receiver [40]. It has been demonstrated that Doppler measurements with a sampling rate higher than 1 Hz are capable of cycle slip detection [41].

Different from the popularity of TDCP-based precise velocity estimation, little research has been carried out on TDCP-based attitude estimation. Although the single antenna GNSS receiver has been reported to determine the attitude, either specially designed receivers are utilized or aiding sensors are integrated [42,43,44]. Based on the precise estimation of the receiver’s delta position, Sun et al. investigated both the TDCP and time difference pseudorange (TDPR) techniques for precise vehicle dynamic heading and pitch estimation [9]. Nevertheless, a high-end GNSS receiver was utilized, and the tests were conducted in good conditions. More contributions could be added to this method using a low-cost GNSS receiver, since it has taken over the majority of mass market applications. Strategies to reduce errors and improve system reliability are indispensable regarding the low-cost receiver applied to a practical environment.

Therefore, in this work, we mainly contributed toward precise heading and pitch angle estimation with a low-cost GNSS receiver for land vehicle applications that benefit from the motion constraint. First, the precise delta position of the receiver between two consecutive epochs is estimated based on the time-differenced carrier phases of a single-antenna GNSS receiver. The vehicular heading and pitch angles are subsequently computed from the components of the delta position. To eliminate the impact of cycle slips, the Doppler-aided cycle slip detection method is applied. Data continuity of the carrier phases is essential for the TDCP method; thus, the exclusion of incomplete data is also analyzed. A real data set of a vehicular dynamic test with a low-cost, single-frequency GNSS receiver is utilized to evaluate the proposed algorithm, and the experimental results show that the root mean square (RMS) of the estimated heading and pitch errors are 1.07° and 1.25°, respectively, with the solution availability being about 80%. Aside from that, the TDCP-based heading and pitch angles can also be utilized as an additional aiding source for the inertial navigation system (INS).

2. Methodology

We first present the overview regarding the proposed algorithm, followed by the detailed derivation of the heading and pitch estimation method based on the TDCP technique. The carrier phases are utilized to estimate the precise delta position of the receiver, and the heading and pitch angles are derived from the estimated delta position afterward. This is followed by the description of cycle slip detection using high-rate Doppler measurements, and carrier phases with possible cycle slips are excluded. Finally, the data continuity of the carrier phase measurements is also discussed, and incomplete data between two consecutive epochs are excluded from attitude estimation.

2.1. System Overview

We structured the data processing flow of the proposed vehicular heading and pitch angles estimation algorithm as shown in Figure 1. The satellite positions computed from the broadcast ephemeris are used for the positioning process and to form the satellite delta position and LOS unit vector. As for the GNSS observations, the code measurements are used for positioning to obtain the approximate receiver positions and to form the LOS unit vector. The carrier phases first go through a check to ensure data continuity, followed by CS detection using the Doppler measurements. Then, the formed TDCP, LOS unit vector, and satellite delta position are utilized to estimate the receiver delta position. Finally, the vehicular heading and pitch angles are computed if the length of the estimated delta position is greater than a certain threshold.

2.2. Carrier Phase-Based Heading and Pitch Estimation

2.2.1. Receiver Delta Position Estimation

The GNSS carrier phase measurement obtained by a receiver can be generally modeled as [45]

$$\Phi =\rho +d\rho +c\left(dT-dt\right)+\lambda N-d_{ion}+d_{trp}+{\epsilon }_{\Phi }$$

(1)

where $\Phi$ is the carrier phase-based range in the unit of length (i.e., meters); $\rho$ is the geometric distance between the receiver and the satellite; $d\rho$ denotes the satellite orbital error; $c$ is the speed of light in a vacuum; $dt$ and $dT$ are the satellite clock bias and the receiver clock offset, respectively, both with respect to the GPS time; $\lambda$ is the carrier wavelength of a specified frequency; $N$ is the related integer ambiguity; $d_{ion}$ and $d_{trp}$ are the ionospheric error and tropospheric error, respectively; and ${\epsilon }_{\Phi }$ includes the noise, multipath error, and other unmodeled errors.

As shown in Figure 2, the geometric distance between the receiver and the satellite at epoch $k$, according to the green triangle, can be expressed as

$${\rho }_k=\Vert {\mathit{r}}_k^s-{\mathit{r}}_{r,k}\Vert ={\mathit{u}}_k·\left({\mathit{r}}_k^s-{\mathit{r}}_{r,k}\right)$$

(2)

$${\mathit{u}}_k=\frac{{\mathit{r}}_k^s-{\mathit{r}}_{r,k}}{\Vert {\mathit{r}}_k^s-{\mathit{r}}_{r,k}\Vert }$$

(3)

where ${\mathit{r}}_k^s$ and ${\mathit{r}}_{r,k}$ are the satellite position vector and the receiver position vector at epoch k, respectively, both expressed in the Earth-centered-Earth-fixed (ECEF) frame, hereafter referred to as the e-frame for brevity; ‖·‖ denotes the norm of a vector; ${\mathit{u}}_k$ is the line-of-sight unit vector (or equivalently the directional cosine vector) at epoch $k$, pointing from the receiver to the satellite; and the dot operator denotes the inner product of two vectors (i.e., $\mathit{u}·\mathit{r}={\mathit{u}}^{\mathrm{T}}\mathit{r}$), with the superscript T representing the transposes of the matrices. Note that in this article, vectors and matrices are denoted in bold italic font, whereas scalars are denoted in italic font for distinction.

Similar to Equation (2), the geometric distance of the same receiver–satellite pair at epoch $k-1$, according to the blue triangle in Figure 2, can be expressed as

$${\rho }_{k-1}={\mathit{u}}_{k-1}·\left({\mathit{r}}_{k-1}^s-{\mathit{r}}_{r,k-1}\right)$$

(4)

where the mathematical notations are defined the same as Equation (2), except for the subscript to distinguish epochs.

The satellite orbital error, satellite clock bias, and atmospheric error presented in Equation (1) are slow varying errors [46], and thus the changes of these errors between two consecutive epochs are negligible if the sampling interval is short (e.g., $\le 1 \mathrm{s}$). Aside from that, the integer ambiguity is a constant if there is no cycle slip at two adjacent epochs. Therefore, the TDCP is calculated by the difference between the carrier phase observations at epoch $k$ and epoch $k-1$:

$$\begin{array}{cc}\delta \Phi & ={\Phi }_k-{\Phi }_{k-1}\\ & \approx {\rho }_k+cd{T}_k-\left({\rho }_{k-1}+cd{T}_{k-1}\right)+\delta {\epsilon }_{\Phi }\\ & =\left({\rho }_k-{\rho }_{k-1}\right)+\left(cd{T}_k-cd{T}_{k-1}\right)+\delta {\epsilon }_{\Phi }\\ & =\delta \rho +c\delta dT+\delta {\epsilon }_{\Phi }\end{array}$$

(5)

where the $\delta$ operator denotes that the difference is made between the two consecutive epochs and $\delta {\epsilon }_{\Phi }$ consists of the residuals caused by neglecting the common errors in addition to the noise and multipath-related error.

According to Equations (2) and (4), the change of the receiver–satellite geometric distance is

$$\begin{array}{cc}\delta \rho & ={\mathit{u}}_k·\left({\mathit{r}}_k^s-{\mathit{r}}_{r,k}\right)-{\mathit{u}}_{k-1}·\left({\mathit{r}}_{k-1}^s-{\mathit{r}}_{r,k-1}\right)\\ & \approx \mathit{u}\cdot \left({\mathit{r}}_k^s-{\mathit{r}}_{k-1}^s\right)-\mathit{u}·\left({\mathit{r}}_{r,k}-{\mathit{r}}_{r,k-1}\right)\\ & =\mathit{u}·\left(\delta {\mathit{r}}^s-\delta {\mathit{r}}_r\right)\end{array}$$

(6)

where in the second line, the approximation $\mathit{u}\approx {\mathit{u}}_k\approx {\mathit{u}}_{k-1 }$ is presumed based on the fact that the two line-of-sight unit vectors ${\mathit{u}}_k$ and ${\mathit{u}}_{k-1}$ are merely the same for land vehicle applications. Although both the satellite and the receiver are moving, the position change of the satellite (about 3 km/s) and that of the receiver (much smaller than that of the satellite) between the two adjacent epochs are negligible when compared with the geometric distance between the satellite and the receiver (about 20,000 km).

By substituting Equation (6) into Equation (5) and transposing the satellite-related term to the left-hand side, we can write the formula as

$$\delta \Phi -\mathit{u}·\delta {\mathit{r}}^s=-\mathit{u}·\delta {\mathit{r}}_r+c\delta dT+\delta {\epsilon }_{\Phi }$$

(7)

where the unknowns are only shown on the right-hand side of the equation, whereas the left-hand side is the compound measurement.

Supposing there are $m$ satellites simultaneously observed by the receiver at two consecutive epochs, then the $m$ equations, like Equation (7), can be organized in matrix-vector form as

$$\mathit{y}=\mathit{A}\mathit{x}+\mathit{v}$$

(8)

where the unknown vector $\mathit{x}={[\begin{array}{cc} {(\delta {\mathit{r}}_r)}^{\mathrm{T}}& c\delta dT\end{array}]}^{\mathrm{T}}$ consists of two parts: the delta position vector of the receiver expressed in the e-frame and the change of the receiver clock offset. The compound measurement vector $\mathit{y}$ and the design matrix A thus can be derived accordingly:

$$\mathit{y}=\left[\begin{array}{c}\delta {\Phi }^1-{\mathit{u}}^1·\delta {\mathit{r}}^1\\ \delta {\Phi }^2-{\mathit{u}}^2·\delta {\mathit{r}}^2\\ ⋮\\ \delta {\Phi }^m-{\mathit{u}}^m·\delta {\mathit{r}}^m\end{array}\right], \mathit{A}=\left[\begin{array}{cc}-{({\mathit{u}}^1)}^{\mathrm{T}}& 1\\ -{({\mathit{u}}^2)}^{\mathrm{T}}& 1\\ ⋮& ⋮\\ -{({\mathit{u}}^m)}^{\mathrm{T}}& 1\end{array}\right]$$

(9)

where $\mathit{v}={[\begin{array}{cccc}\delta {\epsilon }_{\Phi }^1& \delta {\epsilon }_{\Phi }^2& \cdots & \delta {\epsilon }_{\Phi }^m\end{array}]}^{\mathrm{T}}$ includes the measurement noise, multipath error, and residuals of neglecting common errors. The superscript in Equation (9) denotes the element related to a specific satellite $s$ ($s=1,2,\cdots ,m$).

As is commonly known, the classical least squares estimation (LSE) method can be applied to solve the unknown vector:

$$\mathit{x}={({\mathit{A}}^{\mathrm{T}}\mathit{A})}^{-1}{\mathit{A}}^{\mathrm{T}}\mathit{y}$$

(10)

This holds true as long as at least four satellites (i.e., $m \ge 4$) are simultaneously observed at two consecutive epochs.

Equation (9) indicates that the calculation of the precise delta position depends on the satellite position in relation to the ephemeris as well as the receiver position in relation to the positioning solution. To determine the satellite delta positions $\delta {\mathit{r}}^s$ and the LOS unit vectors ${\mathit{u}}^s$ ($s=1,2,\cdots ,m$), the satellite positions at two consecutive epochs and the approximate receiver position at the current epoch are required.

2.2.2. Heading and Pitch Derivation

Distinct from the multi-propeller drones whose attitude may suddenly change, for general land vehicle applications, the vehicular attitude is constrained by the ground motion, assuming no side slip occurs. Therefore, the heading and pitch angles of a land vehicle can be derived from its delta positions. Note that the roll angle cannot be determined in this method, since the delta position vector of the receiver is located within the longitudinal section of the vehicle. As illustrated in Figure 3, the vehicular heading and pitch angles (i.e., $\psi$ and $\theta$), are respectively computed as follows:

$$\psi ={\mathrm{atan}}_2\left(\delta r_{r,E},\delta r_{r,N}\right)$$

(11)

$$\theta =\mathrm{atan}\left(\frac{\delta r_{r,U}}{\sqrt{\delta r_{r,E}^2+\delta r_{r,N}^2}}\right)$$

(12)

where $\psi \in \left(-\pi , \pi \right]$, $\theta \in \left(-\pi /2, \pi /2\right)$; $\mathrm{atan}()$ and ${\mathrm{atan}}_2()$ are the inverse tangent function and the four-quadrant inverse tangent function, respectively; and $\delta r_{r,E}$, $\delta r_{r,N}$, and $\delta r_{r,U}$ are the components of the receiver delta position expressed in the local level east-north-up (ENU) frame (also known as the navigation frame or simply the n-frame), which is denoted as $\delta {\mathit{r}}_r^n$ for clarity.

$\delta r_r^n$ is transformed from the previously estimated receiver delta position vector expressed in the e-frame, denoted as $\delta r_r^e$ herein to distinguish it. A rotation matrix $R_e^n$ is resorted to to realize the coordinate transformation [47,48,49]. The relationship between the e-frame and the n-frame is shown in Figure 4, and the origins and axes of different frames are distinguished by subscripts:

$$\delta {\mathit{r}}_r^n={\mathit{R}}_e^n\delta {\mathit{r}}_r^e$$

(13)

$${\mathit{R}}_e^n={\mathit{R}}_{\mathit{x}}\left(90°-\phi \right){\mathit{R}}_{\mathit{z}}\left(90°+\lambda \right)=\left[\begin{array}{ccc}-\sin \lambda & \cos \lambda & 0\\ -\sin \phi \cos \lambda & -\sin \phi \sin \lambda & \cos \phi \\ \cos \phi \cos \lambda & \cos \phi \sin \lambda & \sin \phi \end{array}\right]$$

(14)

In Equation (14), $\phi$ and $\lambda$ are the geodetic latitude and longitude of the receiver in degrees, respectively; ${\mathit{R}}_{\mathit{x}}\left(i\right)$ denotes the rotation matrix that rotates along the x-axis at $i$ degrees; and ${\mathit{R}}_{\mathit{z}}\left(j\right)$ is defined accordingly.

Considering the denominators of the inverse tangent functions in Equations (11) and (12), it is apparent that the heading and pitch angles cannot be determined if $\delta {\mathit{r}}_r=0$, and large errors may arise when $\delta {\mathit{r}}_r$ is small, which correspond to situations where the land vehicle is stationary and moving slowly, respectively. Therefore, a threshold $T_r$ can be set on the length of the receiver’s delta position to avoid the singularity problem and to resist possible large errors. The GNSS heading and pitch angles are calculated only if $\delta {\mathit{r}}_r$ is in accordance with the condition described by

$$\Vert \delta {\mathit{r}}_r\Vert \ge T_r$$

(15)

where the value of $T_r$ can be determined by analyzing real test data.

2.3. Doppler-Aided Cycle Slip Detection

The GNSS Doppler measurement reflects the instantaneous frequency shift of the receiver-obtained signal relative to the satellite-transmitted signal [34,49]; thus, it is immune to the influence of cycle slips. Supposing there is no cycle slip between two adjacent epochs, then the expected carrier phase at the current epoch can be predicted using the Doppler measurements and the carrier phase measurement at the last epoch:

$${\widehat{\varphi }}_k={\tilde{\varphi }}_{k-1}-\frac{{\tilde{D}}_k+{\tilde{D}}_{k-1}}{2}\delta t$$

(16)

where $\varphi$ is the carrier phase in cycles; $D$ denotes the Doppler measurement as a frequency (i.e., hertz); the subscripts $k$ and $k-1$ indicate the two adjacent epochs; the overhead scripts ~ and ^ on symbols represent the observed and computed quantities, respectively; and $\delta t$ denotes the sampling interval. As previously defined, let $\delta$ denote the difference operator between two consecutive epochs. Thus, the second term on the right-hand side of Equation (16) represents the estimated carrier phase variation using Doppler measurements:

$$\delta \widehat{\varphi }=-\frac{{\tilde{D}}_k+{\tilde{D}}_{k-1}}{2}\delta t$$

(17)

where the negative sign is due to the Doppler measurement being defined as positive when the satellite and the receiver are approaching each other. As a result, the carrier phase is decreased.

Based on Equations (16) and (17), the observed-minus-computed (OMC) carrier phase residual at epoch $k$ can be written as

$${\tilde{\varphi }}_k-{\widehat{\varphi }}_k={\tilde{\varphi }}_k-\left({\tilde{\varphi }}_{k-1}+\delta \widehat{\varphi }\right)=\delta \tilde{\varphi }-\delta \widehat{\varphi }$$

(18)

Theoretically, supposing no cycle slip occurs at epoch $k$ and epoch $k-1$, the result of Equation (18) should be a small value if not exactly zero due to noise. $\delta \tilde{\varphi }={\tilde{\varphi }}_k-{\tilde{\varphi }}_{k-1}$ is the variation of the measured carrier phase. Therefore, a threshold $T_{\varphi }$ on the difference of $\delta \tilde{\varphi }$ and $\delta \widehat{\varphi }$ can be selected for cycle slip detection. A possible cycle slip is claimed if it corresponds to the condition of Equation (19):

$$\left|\delta \tilde{\varphi }-\delta \widehat{\varphi }\right|>T_{\varphi }$$

(19)

where the value of $T_{\varphi }$ is considered the sum of the measurement noise and modeling errors.

It is worth mentioning that the variation of the carrier phase would decrease as the sampling interval decreased, and the related OMC residual would decrease accordingly. Therefore, a higher sampling rate is helpful for the accuracy improvement of the Doppler-aided cycle slip detection method.

To illustrate the Doppler-aided cycle slip detection method, an example based on the GPS satellite with the pseudorandom number (PRN) G01 is depicted in Figure 5. The carrier phases and Doppler measurements of the L1 signal are utilized to evaluate the algorithm. This shows that the difference of the measured and estimated carrier phase variations varies from time to time, and it is mostly limited to the range of ±0.2 cycles. However, there are five exceptions, with differences of about 0.5 cycles (epochs 767, 887, 1427, 1847, and 4433). These exceptions are considered possible cycle slips that occurred. Therefore, the carrier phase observations at the related epochs are excluded for attitude estimation.

2.4. Incomplete Data Exclusion

Another problem that needs to be considered is the data continuity of the carrier phase observations, which is important to guarantee smooth attitude estimation solutions. Therefore, satellites that drop in and out of view between two consecutive epochs should be excluded; otherwise, some unexpected large spikes in the estimated attitude may appear frequently due to the incomplete data. For instance, if a satellite is visible at epoch $k$ but is lost track of at epoch $k+1$, then the carrier phase of that satellite should be temporarily excluded for attitude estimation.

The carrier phase and code measurements are obtained in different loops [50]. The acquisition and tracking of the carrier phase are generally more difficult compared with those of the code, especially when the receiver is operated under a challenging environment. The visibility of the satellite observations of a vehicular dynamic test is depicted in Figure 6, and the number of epochs with available code and carrier phase measurements are counted in

Table 1. Number of epochs of the code and carrier phase observations, counted per satellite, and the corresponding percentages relative to the total epochs.

Satellite PRN	Code		Carrier Phase
	Epochs	Percentage	Epochs	Percentage
G01	4738	98.7%	4400	91.6%
G10	4594	95.7%	4495	93.6%
G11	4718	98.3%	4588	95.6%
G21	4761	99.2%	3961	82.5%
G22	4699	97.9%	4097	85.3%
G31	4755	99.0%	3815	79.5%
G32	4761	99.2%	4745	98.8%

, together with the related percentages relative to the total epochs, which were also calculated. It should be noted that only satellites with elevation angles greater than 15° were considered.

Based on Figure 6 and Table 1, it can be found that there were some epochs where signals from all satellites showed a loss of lock (e.g., from epoch 3210 to epoch 3225), as both the code and carrier phase observations at those epochs were interrupted. Among all satellites, the availability of code observations was greater than that of the carrier phase observations. However, a satellite with higher availability for its code observations did not necessarily guarantee a higher availability of its carrier phase observations, which indicates that they were not correlated. Satellite G32 showed the best data availabilities for both kinds of observations, whereas satellite G31 showed the lowest carrier phase availability, although its code availability was pretty good.

The properties of the carrier phase availability could be explained by considering the elevation angles of the selected satellites. As shown in Figure 7, there were seven selected satellites in total, which could be straightforwardly divided into three groups: satellite G32 alone as the first group, which owned the highest elevation; satellites G01, G10, and G11 as the second group, with medium relative elevation angles; and the remaining satellites as the third group with relatively low elevation angles, in which the satellite G31 had the lowest elevation. This demonstrates that the availability of carrier phase observations was affected by the satellite elevation, whereas that of the code observations was not. Aside from that, the code observations of particular satellites contained similar epoch numbers; however, this was not the case with the carrier phase observations.

3. Experimental Results and Analysis

The experimental set-up is presented first, including descriptions related to the utilized devices, testing trajectory, and collected data set, followed by the analysis of the approximate receiver positions that were obtained from the single-point positioning. Then, the vehicular heading and pitch angels are estimated based on the TDCP technique, and two additional methods are implemented for comparison. Finally, the influence of cycle slips on the heading and pitch estimation results is discussed.

3.1. Experimental Set-Up and Data Collection

As shown in Figure 8, an experimental platform based on a land vehicle was built to evaluate the performance of the proposed algorithm. A low-cost, single-frequency GNSS receiver (i.e., the ublox NEO-M8T) was used as the primary receiver to collect the experimental data. A high-performance synchronized position attitude navigation (SPAN) system from Novatel was utilized to provide the reference solutions, which consisted of an OEM6 receiver and a KVH 1750. Note that the OEM6 is a high-performance GNSS receiver and the KVH 1750 is a tactical-grade inertial measurement unit (IMU) with fiber optic gyros (FOGs) inside. A Trimble R9 receiver installed on the roof of the engineering building at the University of Calgary (UC) was also recorded as the base station. The SPAN system combined with the base station was post-processed in the GNSS/INS tightly integrated mode, which was able to provide high-performance position, velocity, and attitude reference solutions.

The vehicular experiment was carried out on 2 August 2020 in Calgary, Canada. The vehicle started and ended at two different parking lots on campus at UC, and the test trajectory included driving along straight lines, turns, and stops at traffic lights. During the test, the ublox NEO-M8T and OEM6 were connected to the same GNSS antenna fixed on the car roof via a signal splitter. The data set of the ublox was collected at a 10-Hz sampling rate, and only the measurements of the GPS L1 signal were utilized. The traveled distance of the car was approximately 4.4 km with a maximum speed of 52 km/h, and 8 min of data were extracted from the original data set by partially discarding the stationary sections at the start point and the end point, since the proposed algorithm could not estimate the heading and pitch angles of the vehicle with stationary data. In addition, the reference solutions were extracted and aligned with the corresponding test data set according to the GPS time.

3.2. Positioning Results

As derived in the methodology section, the estimation of the precise delta position of the receiver relied on the satellite position and the approximate receiver position. Therefore, the positioning process was implemented first based on the code measurements to obtain the approximate receiver position. In the particular processing flow, the satellite clock bias and the relativistic effect were corrected using the broadcast ephemeris. The ionospheric error and the tropospheric error were compensated by the Klobuchar model and the Saastamoinen model, respectively. Satellites with elevation angles lower than 15° were temporarily excluded as long as the situation remained.

The signal-to-noise ratio (SNR) variations of the selected satellites are presented in Figure 9. This shows that the SNR values of the satellites with lower elevation angles, as depicted in Figure 7, were generally lower and fluctuated more acutely compared with the SNR values of the satellites with higher elevation angles. Aside from that, by considering the vehicular dynamics and the test environment, it was found that the changes in the SNR values in the stationary stages and open-sky environment were much better than those in the dynamic period or across high buildings.

We found that setting an SNR mask was helpful to improve the positioning accuracy. Consequently, based on trial and error, satellites with SNRs lower than 35 dBHz at a specific epoch were temporarily rejected. Note that the selection of this threshold was a trade-off between the availability and accuracy of the solutions. A higher SNR threshold implied more satellites being excluded, which caused more epochs without enough satellites (less than four) for positioning, and therefore, the solution availability was decreased. On the contrary, a lower SNR threshold took more satellites with bad observations into consideration, which would reduce the positioning accuracy as a result.

The obtained trajectory of the positioning solution is depicted in Figure 10. The real vehicular dynamic test shows that there were several outages throughout the trajectory, which indicates that the test trajectory was under a typical suburban environment, consisting of both the open-sky and partially sheltered situations. The outages included those caused by full signal blockages and those due to the exclusion of satellites with bad SNR values.

The positioning results and the related position errors are presented in Figure 11, with both expressed in the n-frame. Note that the position errors were computed by comparing the SPP solutions with the reference solutions. The statistical properties of the SPP solutions including the RMS, absolute maximum error, and solution availability are listed in

Table 2. Statistical properties of the positioning results for the root mean square (RMS) error, absolute maximum (Max) error, and solution availability.

RMS (m)			Max (m)			Solution Availability
E	N	U	E	N	U
0.54	0.73	1.43	18.43	14.16	48.45	97.8%

. This shows that although the maximum position errors, especially for the vertical direction, were large, the RMS of the horizontal and vertical errors was less than 1 m and 1.5 m, respectively. The positioning solution was not always available since there were some epochs where the number of visible satellites was less than four, due to the combined effect of the signal blockage and satellite exclusion. Taking the vehicular dynamics into consideration, it was found that the position errors at static periods were generally smaller than those at dynamic periods. The long gap of outages and the abrupt large position errors between epochs 2500 and 2800 were caused by poor observations, as the SNR values fluctuated dramatically during the corresponding period.

3.3. Heading and Pitch Estimation Results

With the obtained approximate receiver positions from the positioning process, the TDCP-based GNSS heading and pitch angles could be estimated afterward. In the processing flow, the incomplete data were rejected first as described in the methodology section. Then, the Doppler-aided cycle slip detection method was conducted. The threshold for cycle slip detection $T_{\phi }$ was set to 0.4 cycles based on the analysis of the test data, and thus the related carrier phase with an absolute OMC residual greater than 0.4 cycles was excluded. After the precise delta position of the receiver was estimated, the threshold for the delta position $T_r$ was set to 0.1 m, which was equivalent to a velocity constraint of 1 m/s or 3.6 km/h. Delta positions with a length smaller than $T_r$ were refused for the heading and pitch estimation.

The three data exclusion steps reduced the solution availability but improved the accuracy and smoothness of the heading and pitch estimates. The time series of the number of valid satellites is presented in Figure 12, where a valid satellite means that its carrier phase observations pass all three data exclusion steps. It is shown that the number of valid satellites changed frequently. The maximum valid satellite number was seven, since akin to the positioning process, only satellites with elevation angles greater than 15° were considered. Aside from that, there were several epochs with less than four valid satellites, which would result in the inability to estimate the precise delta position of the receiver and, consequently, the inability to provide the heading and pitch estimates at these epochs.

In addition to the proposed TDCP-based vehicular heading and pitch angle estimation method, two extra algorithms were also used for comparison: the direct position difference of the single-point positioning (DSPP) solutions and the TDPR-based heading and pitch angle estimation method. Note that the TDPR-based method is much like the TDCP-based method; the only difference is that the code-based pseudorange measurements instead of the carrier phases are utilized for the receiver delta position estimation. The estimated vehicular heading and pitch angles of the three algorithms are compared in Figure 13, with the related attitude errors presented in Figure 14. The RMS, absolute maximum error, and solution availability of the attitude errors are listed in

Table 3. RMS, absolute maximum error, and solution availability of the heading and pitch estimates from the DSPP, TDPR, and TDCP methods.

Algorithm	RMS (°)		Max (°)		Solution Availability
	$\mathit{\psi }$	$\mathit{\theta }$	$\mathit{\psi }$	$\mathit{\theta }$
DSPP	18.58	18.06	171.44	85.36	84.1%
TDPR	17.56	18.54	179.01	87.30	84.3%
TDCP	1.07	1.25	5.99	12.87	79.7%

. Note that the heading and pitch angles from the high-performance SPAN system were utilized as the reference.

Among the three algorithms used for comparison, the differences included two parts: the delta positions of the receiver between consecutive epochs were computed from different data, and the data availabilities were different accordingly. The results show that the heading and pitch estimates from the DSPP method and those from the TDPR-based method embodied similar unacceptably large errors, considering that these two methods are based on the pseudorange measurements. The RMS values of the heading and pitch angles were around 18° for both methods, while the absolute maximum errors for the heading and pitch were larger than 170° and 80°, respectively. Although the estimated angles generally followed the references, the tremendous fluctuations were presented at multiple epochs. More epochs with fluctuations were experienced by the DSPP method compared with the TDPR method, especially for the heading estimates. This indicates that both the DSPP and TDPR methods were unsuitable for vehicular heading and pitch angle estimation.

Only the TDCP method resulted in reasonable heading and pitch estimates, with RMS values better than 1.5° (i.e., more than ten times better than those from the other two methods). The absolute maximum errors of the TDCP method were also much smaller compared with those of the DSPP and TDPR methods. The vast improvement of the attitude estimation accuracy was due to the fact that the carrier phases were much more accurate and more robust for multipath errors than the pseudorange measurements. However, the TDCP method held the smallest solution availability, since fewer carrier phases were observed than for the pseudorange measurements, and carrier phases with cycle slips are excluded. The solution availabilities of the other two methods were similar, as the availability of the TDPR method was slightly larger than that of the DSPP method, since the pseudorange measurements were directly utilized in the TDPR method to estimate the delta position of the receiver.

3.4. Influence of Cycle Slips

To show the influence of cycle slips on the vehicular heading and pitch angle estimation, the TDCP-based algorithm was evaluated again but without cycle slip exclusion. The comparison of the heading and pitch errors of the TDCP method with and without cycle slips are presented in Figure 15, and the related statistical properties are listed in

Table 4. RMS, absolute maximum error, and solution availability of heading and pitch estimates from TDCP with and without cycle slips.

Cycle Slip Status	RMS (°)		Max (°)		Solution Availability
	$\mathit{\psi }$	$\mathit{\theta }$	$\mathit{\psi }$	$\mathit{\theta }$
CS included	1.31	2.45	22.90	58.93	81.3%
CS excluded	1.07	1.25	5.99	12.87	79.7%

. Compared with the solutions of the TDCP method where the possible cycle slips were detected and excluded, there were several large spikes in the TDCP method without cycle slip exclusion both for the heading and pitch estimates. Thus, the absolute maximum errors also increased a lot. The RMS of the heading estimates increased approximately 20%, while that of the pitch estimates almost doubled.

However, more available solutions were obtained, since more epochs (i.e., those with possible cycle slips specifically) were involved in the estimation of the receiver delta position. This shows that cycle slip detection and exclusion is an indispensable part of the proposed TDCP-based vehicular heading and pitch angle estimation algorithm. Otherwise, some unexpected large estimation errors would occur.

4. Discussion

GNSS observations under a challenging environment suffer from signal blockage, reflection, and attenuation, especially when a low-cost receiver is operated in dynamic applications. The satellite elevation angle and SNR value of the received signal are two indicators that reflect the data quality, but the availability of code observations is not necessarily consistent with the availability of carrier phase observations. Regarding the positioning test, it was found that setting elevation mask and SNR mask was necessary to reject errors and guarantee the positioning accuracy. Aside from that, a higher data sampling frequency was good for not only the Doppler-aided cycle slip detection but also the TDCP-based heading and pitch estimation. However, the carrier phases with possible cycle slips were excluded in the current research. The introduction of a cycle slip repair technique would increase the solution availability. Aside from that, precise point positioning could be implemented instead of the standard point positioning to improve the positioning accuracy, which would be of benefit to improve the accuracy of attitude estimation.

In terms of the vehicular heading and pitch angle estimation with a low-cost GNSS receiver, the test results show that only the TDCP-based algorithm could provide an acceptable estimation accuracy, whereas the results of the DSPP and TDPR methods showed large variations, although the solution availabilities of these two methods were higher than that of the TDCP method. It was also found that the exclusion of incomplete data and carrier phases with possible cycle slips was essential to avoid gross errors in the TDCP method. However, since the heading and pitch estimations were derived from the receiver’s delta position, the TDCP technique could not work properly if the vehicle was static or moving slowly, and the solutions were not continuous. These two limitations can be further mitigated by integration with aiding sensors (e.g., the inertial sensors) which apply the continuous attitude estimation technique.

5. Conclusions

A precise heading and pitch estimation algorithm based on the time-differenced carrier phases (TDCPs) of a single-antenna GNSS receiver is proposed. First, the precise delta positions of the receiver between two consecutive epochs are estimated by the TDCP technique. The heading and pitch estimates are derived afterward according to the components of the receiver delta position, where a threshold on the length of the delta position is set to avoid large errors. Since the cycle slip problem is a threat to the TDCP technique, the Doppler-aided cycle slip detection method is used to temporarily exclude satellites with possible cycle slips. Aside from that, the continuity of the carrier phases is also considered, and satellites with incomplete data are temporarily refused.

To assess the performance of the proposed algorithm, we carried out a car-borne dynamic experiment in a typical suburban environment. The experimental data were collected by a low-cost single-frequency GNSS receiver (i.e., the ublox NEO-M8T). In terms of the recorded observations, there were epochs with all satellite signals showing loss of lock, but the availability of the code observations was better than the availability of the carrier phases. With the reference solutions provided by a high-performance SPAN system, the root mean square error (RMS) of the positioning results was 0.54 m, 0.73 m and 1.43 m in the east, north, and vertical directions, respectively. The heading and pitch estimation results demonstrated the feasibility of the proposed TDCP-based algorithm for land vehicular applications, as the RMS of the attitude errors was less than 1.5°, whereas massive fluctuations existed in the solutions of the direct position difference of the single-point positioning (DSPP) and time-differenced pseudorange (TDPR) methods. In addition, cycle slip detection and exclusion was an indispensable part of the TPCP-based algorithm to avoid some unexpected large errors.