A Design of Differential-Low Earth Orbit Opportunistically Enhanced GNSS (D-LoeGNSS) Navigation Framework

Abstract

Considering the problem of GNSS service interruption caused by the insufficient number of available satellites in complex environments, Low Earth Orbit (LEO) satellites can supplement GNSS effectively. To eliminate the unknown satellite clock error and the atmospheric delay error with spatial correlation in LEO observations, a Differential-Low Earth Orbit opportunistically enhancing GNSS (D-LoeGNSS) navigation framework is proposed. Firstly, because of the uncertainty of the LEO orbit, we derive the effect of the LEO orbit error on the differential measurement model. Secondly, aiming at the noise amplification and correlation in double-difference (DD), we propose a Householder-Based D-LoeGNSS (HB-DLG) algorithm, which suppresses noise by introducing an orthogonal matrix. Thirdly, in D-LoeGNSS, the typical measurement of LEO is Doppler, which is heterogeneous with the GNSS pseudorange, rendering the Dilution of Precision (DOP) evaluation method unsuitable. Given the unbiasedness of differential measurements, the Cramer Rao Lower Bound (CRLB) is derived as a metric to characterize the positioning accuracy and satellite spatial distribution. Finally, a field experiment using Orbcomm (ORB) and GPS is conducted. The experimental results show that the performance of the HB-DLG algorithm is superior to DD. Especially when the number of satellites is insufficient or the measurement redundancy is poor; the D-LoeGNSS framework has advantages of rapid convergence and high accuracy compared with a single constellation.

1. Introduction

Over the past decade, the Global Navigation Satellite System (GNSS), represented by GPS, GLONASS, BeiDou and Galileo, has been the primary means of providing navigation and positioning services. However, GNSS is inherently defected by low received signal power and is easily blocked resulting in frequent service interruptions or even unavailability in challenging environments. Therefore, improving GNSS performance in complex environments has been the main driving force for the development of satellite navigation technology.

Nowadays, various Satellite-Based Augmentation Systems (SBAS) have been successively proposed which can be divided into information enhancement and signal enhancement [1,2]. Information enhancement aims to improve positioning accuracy and reliability by broadcasting error correction in GNSS observations and integrity information to users. It does not provide additional measurements. Typical systems include the Wide Area Augmentation System (WAAS), and the BeiDou regional correction broadcasted by the Geostationary Earth Observation satellites [3]. Nevertheless, since they generally use medium or high-orbit satellites, the availability of location services in complex environments has not been ensured.

In contrast, signal enhancement can provide additional measurement to GNSS signals and improve positioning availability by increasing measurement redundancy. Low Earth Orbit (LEO) satellites are a typical representation of signal enhancement, possessing ideal properties for aiding GNSS positioning in complex environments. Specifically, LEO is about 20 times closer to the earth, making the received signal power 30 dB higher than GNSS signals and less susceptible to interference. Furthermore, companies such as OneWeb, SpaceX, and Telesat plan to deploy tens of thousands of broadband internet satellites into space, increasing the number of LEO satellites significantly. Besides, each broadband provider deploys broadband internet satellites into specific constellations and broadcasts signals on different frequency bands, making LEO signals diverse in frequency and difficult to be targeted. Therefore, using LEO satellites for GNSS navigation enhancement, known as LEO enhanced GNSS (LeGNSS), has become a current research hotspot [4,5,6,7].

The working mode of LEO satellites in the LeGNSS can also be categorized into two types. The first is LEO navigation satellites, which broadcast dual-frequency pseudorange and carrier-phase observations. Leveraging the significant spatial geometry changes of LEO satellites in a short period, the correlation between equations across epochs is weakened, solving the problem of long initialization times for GNSS carrier-phase ambiguity in PPP. By designing different constellation configurations and the number of satellites, LeGNSS can shorten the PPP convergence time to seconds [8,9,10]. However, according to [11] existing broadband internet LEO constellations may not support pseudorange-based navigation for two reasons: First, due to cost constraints they will not carry atomic clocks. Second, they will not broadcast downlink signals that support pseudorange measurements through a calibrated connection between transmission time tags and transmitter clock time. Further, various PPP error models for LEO still adopt empirical models from GPS, so the research under this mode is still at the simulation stage.

Another mode of LEO is opportunistically to assist GNSS. Opportunistic fashion is a method of extracting useful observational information from public radio signals intended for non-navigation purposes, i.e., Signal-of-Opportunity (SOPs). Compared with the first mode, the opportunistic approach has additional advantages: (i) It protects the privacy of the user, as only the downlink LEO satellite signal is used without communicating back with the LEO; (ii) it has better universality, which can be extended to any number of LEO constellations without being limited to a single constellation [12]. There are few studies using LEO satellites to assist GNSS opportunistically. The authors of [13] combined Orbcomm (ORB) with GPSs, and analyzed the impact of different observation errors on positioning accuracy. Since the LEO clock error was not considered, the positioning accuracy was poor. The authors of [14] proposed to use the Doppler measurement extracted from the Iridium single-tone signals combined with GPS for positioning, and a simple experiment was carried out. In particular, the following challenges should be considered when considering LEO as SOP beacons:

• Limited information about LEO satellite signals is available. As [15,16] mentioned, broadband internet constellations such as Starlink generally do not disclose signal structure specifications to protect intellectual property. Therefore, it is difficult to construct pseudoranges by decoding pseudo-random sequences, and the typical measurement is Doppler;
• The specific parameters for the LEO error model are unknown. Although the Doppler error model can be regarded as the time differential of the pseudorange correction model, its accurate equipment parameters and environmental parameters cannot be obtained. In addition, due to cost reduction, LEO onboard clocks use pace-qualified OCXO and TXCO, which are far less stable than GNSS atomic clocks [17,18]. Therefore, the LEO satellite clock drift cannot be ignored;
• The uncertainty of LEO satellite orbits. As the LEO constellations do not have the massive network of ground stations that form the control segment of GNSS. The precise ephemeris of the LEO satellites is not broadcast or known. Although the North American Aerospace Defense Command (NORAD) makes publicly available estimates of Keplerian elements for LEO satellite orbits that are updated daily in two-line element (TLE) files. However, TLE is an inaccurate ephemeris. The satellite state calculated by TLE files has been proven to accumulate an error of one kilometer or more within a day [19].

The differential framework was introduced to remove unmodeled common errors; performance in LEO constellations has been partially studied recently. The authors of [20] proposed an altimeter-assisted Orbcomm carrier-phase double-difference (DD) framework. Using two Orbcomms, the 2D position error of the user station reaches 11.93 m after waiting for 114 s. Since the Iridium signals use TDMA/FDMA modulation, the authors of [21] used Iridium as the SOP beacon and measured TDOA and FDOA as the observation. However, the Iridium downlink code rate is only 25 kHz, and the measurement accuracy of TDOA is limited. Finally, the 2D positioning accuracy reaches 29.1 m within 50 min. The authors of [22] proposed to construct the line-of-sight vector in the Jacobian matrix by measuring the AOA of the LEO signal, thereby circumventing orbital error on positioning accuracy. However, AOA requires complex antenna arrays to obtain accurate relative azimuth angles. The authors of [23] presented the preliminary results of differential Doppler positioning using Starlink satellites, and designed a Doppler tracking algorithm based on a generalized likelihood ratio (GLR) detector. Under the 1 km and 9 m baselines, the user’s 2D position error is 5.6 m and 2.6 m, respectively. The authors of [24] developed a navigation framework with carrier-phase differential measurements from the mega-constellation LEO satellites. The simulation results showed that using 44 Starlink satellites, an unmanned aerial vehicle (UAV) achieved a 3D position RMSE of 2.2 m, and a 2D RMSE of 32.4 cm. However, this precision comes at the cost of added ambiguities that need to be resolved. Moreover, in a complex environment the carrier-phase is prone to cycle slip and loss of lock, resulting in frequent initialization of ambiguities. In contrast, LEO Doppler information has better availability.

This paper proposes a Differential-LEO opportunistically enhanced GNSS (D-LoeGNSS) navigation framework when there is an insufficient number of GNSS satellites and evaluates the performance of the framework. A reference station whose state is known a priori and an unknown user station observes GNSS and LEO simultaneously. Compared with the single system which has to observe for a long time until the satellite geometry changes sufficiently, LoeGNSS can effectively improve the availability and real-time performance of the system. The main contributions of this paper are as follows:

• The influence of atmospheric delay on LEO Doppler measurement is analyzed;
• A D-LoeGNSS navigation framework is proposed to remove unknown satellite clock errors and atmospheric delays with spatial correlation, and the measurement residuals caused by LEO ephemeris errors are derived;
• Based on the orthogonal transformation [24], we proposed a householder-based D-LoeGNSS algorithm (HB-DLG) for the problems of high measurement noise correlation and prominent noise in DD. The noise amplification caused by the difference is suppressed by introducing an orthogonal matrix;
• Since the GNSS pseudorange and the LEO Doppler are heterogeneous measurements, DOP is unsuitable for characterizing positioning accuracy. Although [11,25,26] proposed the concept of generalized GDOP, which unifies the units of Doppler and pseudorange by normalizing the scale factor, whose computational complexity is high. Considering that the measurements error after difference only contains white noise, which satisfies the assumption of unbiased estimation, the Cramer Rao Lower Bound (CRLB) is introduced into the difference-system. We derived the CRLB of HB-DLG as a metric reflecting the spatial distribution of satellites in the LoeGNSS constellation and the positioning accuracy.

The rest of the paper is organized as follows: Section 2 presents the principles and methods. Section 2.1 presents the D-LoeGNSS framework. Section 2.2 describes the D-LoeGNSS framework single difference (SD) measurement model. To remove the receiver clock bias and clock drift, Section 2.3 describes the D-LoeGNSS framework DD measurement model and the problem of noise amplification. Section 2.4 presents the HB-DLG measurement model. Section 2.5 derives the CRLB as the positioning accuracy evaluation metric for the LoeGNSS constellation. In Section 3, we take Orbcomm and GPS as examples to conduct field experiments and present the results of the D-LeoGNSS framework. Finally, the conclusion is given in Section 4.

2. Principles and Methods

This section presents the D-LeoGNSS framework and its measurement model. In particular, the influence of atmospheric delay on LEO Doppler measurement was analyzed. Then, the measurement residuals caused by LEO ephemeris errors were derived. Given the HB-DLG measurement model, we derived the CRLB as the positioning accuracy evaluation metric for the LoeGNSS constellation.

2.1. D-LoeGNSS Framework Description

The framework consists of a reference station with a known position and a user station, operating in space composed of M ($M<4$) GNSS satellites and N ($N\ge 2$) LEO satellites. Each station simultaneously measures LEO and GNSS. The reference broadcasts its position, GNSS/LEO satellite state, and GNSS pseudorange/LEO Doppler measurements to the user. The GNSS satellite identification is determined by the cyclic convolution of the local C/A code sequence and the acquisition signal. For LEO, this can be accomplished by allowing the user and reference to first learn about their operating environment with preloaded information. For example, each LEO corresponds to a specific earth coverage area at a given time, and its orbital information can be predicted by preloading the TLE file combined with the SGP4 model [27]. Therefore, the preliminary screening of satellites is completed by detecting whether the satellite coverage area intersects with the reference station area. Afterwards, the Doppler measurements at the reference are matched with the theoretical Doppler to complete the LEO satellite identification. For the Orbcomm constellation, each satellite broadcasts signals in its own two specific frequency bands so that the satellites can be determined through spectral scanning. Once the satellites are identified, the receiver can calculate satellite orbital information at any given time through the SGP4 model.

In particular, to remove unknown satellite clock errors, atmospheric delays, and orbital errors in LEO measurements, we used a difference framework containing SD and DD. To suppress the noise amplification and correlation in DD, the HB-DLG measurement model was proposed. The experimental scene is shown in Figure 1.

2.2. D-LoeGNSS Framework Single Difference Measurement Model

2.2.1. Effect of Ionospheric Change Rate on LEO Doppler Measurements

The ionosphere is a dispersive medium for electromagnetic waves, ranging from 60 to 2000 km above the earth’s surface [28]. GNSS orbits are typically higher than 2000 km, and GNSS signals tracked by ground receivers penetrate all ionospheres. Ionospheric effects are mainly determined by the total electron content (TEC), whereas most LEO satellites orbit between 500 and 1500 km, i.e., within the ionosphere. Compared to GNSS, which is affected by the entire TEC, LEO satellites are only affected by the bottom TEC [29]. Since the specific bottom TEC parameter is unknown, we analyzed the impact of the ionosphere on LEO Doppler using the GPS empirical model. According to [29,30], the ionospheric delay in the pseudorange in meters is:

$$I^{\left(i\right)}=\frac{\kappa \times VTE{C}^{\left(i\right)}\times {10}^{16}\times {\mathrm{F}}^{\left(i\right)}\left(\theta \right)}{f_i{}^2}$$

(1)

where $f_i$ is the carrier frequency of satellite i, $\kappa$ is the proportionality constant which reads: $\kappa =q^2/\left(8{\pi }^2{m}_e{\epsilon }_0\right)=c^2{r}_e/\left(2\pi \right)$ where $q,m_e,r_e$ are the electron charge, mass, and radius, respectively; c is the vacuum speed of light and ${\epsilon }_0$ is the vacuum permittivity. The value of the constant is approximately $\kappa$ ≈ 40.308193 m³·s⁻² [31]. $VTE{C}^{\left(i\right)}$ is expressed vertical total electron content of satellite i and is assumed to be constant during satellite visibility, $\theta$ indicates the elevation angle of satellite i relative to the user, and ${\mathrm{F}}^{\left(i\right)}\left(\theta \right)$ is the obliquity factor of satellite i, namely:

$${\mathrm{F}}^{\left(i\right)}\left(\theta \right)={\left(1-{\left(\frac{{\mathrm{R}}_{\mathrm{e}}\cdot \cos \left(\theta \right)}{{\mathrm{R}}_{\mathrm{e}}+h}\right)}^2\right)}^{-\frac{1}{2}}$$

(2)

where R_e = 6368 km is the average radius of the earth, h is the mean ionospheric height, generally 350 km. As the Doppler can be regarded as the time derivative of pseudorange pairs, then:

$${\dot{I}}^{\left(i\right)}=-\frac{\kappa \times VTE{C}^{\left(i\right)}\times {10}^{16}}{f_i{}^2\times {\left[1-{\left(\frac{{\mathrm{R}}_{\mathrm{e}}}{{\mathrm{R}}_{\mathrm{e}}+h}\cdot \cos \left(\theta \right)\right)}^2\right]}^{\frac{3}{2}}}\times \frac{{\mathrm{R}}_{\mathrm{e}}\sin \left(2\theta \right)}{2\left({\mathrm{R}}_{\mathrm{e}}+h\right)}\cdot \frac{\partial \theta }{\partial t}$$

(3)

According to (3), the ionospheric change rate is positively correlated with the elevation change rate and inversely proportional to the square of the satellite carrier frequency. We simulated GPS PRN2, Starlink 4272, and Orbcomm FM114 satellites as examples. GPS, Starlink, and Orbcomm use the L1 frequency band, 11.325 GHz and 137.5 MHz, respectively. Since the typically visible period for LEO was around 10 min, we used STK to obtain satellite states; the simulation step was set to 1 s and the cutoff elevation angle was 15 deg. Specific meteorological parameters are provided in [32,33]. The results are shown in Figure 2.

In Figure 2, although the specific values of the ionospheric change rate may differ from the actual situation, due to the high carrier frequency of GPS and the elevation changes slowly, the influence of the ionospheric change rate can be ignored. Starlink has the highest signal frequency, even its elevation angle changes rapidly, the ionospheric change rate value is slightly higher than the GPS. The Orbcomm signal was most affected by the ionosphere change rate since it had the lowest frequency with quick changes in the elevation angle.

2.2.2. Effect of Tropospheric Change Rate on LEO Doppler Measurements

The troposphere is the atmosphere extending from the ground to 50 km, and is a non-dispersible medium. Both LEO and GNSS were affected by the troposphere. The troposphere is divided into two components called dry and wet delay. To analyze the change rate of tropospheric delay, we used the Hopfield model as an example:

$$\Delta {\mathrm{D}}_{tro}^z=\Delta D_{dry}^z+\Delta D_{wet}^z$$

(4)

where $\Delta {\mathrm{D}}_{tro}^z$ is the zenith tropospheric delay; $\Delta D_{dry}^z$ and $\Delta D_{wet}^z$ is the dry delay and wet delay, respectively:

$$\begin{array}{l}\Delta D_{dry}^z=1.552\times {10}^{-5}\times \frac{P}{T}\times \left(h_{dry}-h_0\right)\\ \Delta D_{wet}^z=7.46512\times {10}^{-2}\times \frac{e_0}{T^2}\times \left(h_{wet}-h_0\right)\end{array}$$

(5)

In (5), P is the local atmospheric pressure in mPa, T represents temperature in Kelvin, $h_{dry},h_{wet}$ represent the heights of the dry and wet atmospheric layers, respectively. Additionally, it recommends the following empirical values:

$$\begin{array}{l}{h}_{dry}=40136+148.72\left(T-273.16\right)\\ h_{wet}=11000\\ e_0=RH\times {10}^{\frac{7.5\times \left(T-273.3\right)}{T}}\end{array}$$

(6)

$e_0$ is vapor pressure, and RH represents relative humidity. The zenith tropospheric delay component is multiplied by the corresponding obliquity factors to characterize the tropospheric delay and change rate in the direction of signal transmission:

$$\Delta D_{tro}^{}=F_{dry}\Delta D_{dry}^{}+F_{wet}\Delta D_{wet}^{}$$

(7)

The obliquity factors are:

$$\begin{array}{l}{F}_{dry}=\frac{1}{\sin (\sqrt{{\theta }^2+{\left(\frac{2.5\pi }{180}\right)}^2})}\\ F_{wet}=\frac{1}{\sin (\sqrt{{\theta }^2+{\left(\frac{1.5\pi }{180}\right)}^2})}\end{array}$$

(8)

Finally, the tropospheric delay and its change rate for GPS, Starlink, and Orbcomm are shown in Figure 3. Note that the measured values for the change rate were obtained by differencing the tropospheric delays for adjacent epochs.

In Figure 3, the change rate of the tropospheric delay for Orbcomm and Starlink was greater than that of the GPS, but their magnitudes were still small. In particular, Orbcomm was more affected by the ionosphere than the troposphere.

2.2.3. Single Difference Measurement Model

LEO is equipped with a pace-qualified OCXO or temperature-compensated crystal oscillator (TCXO) as the onboard clock, whereas its specific clock parameters cannot be obtained, therefore, we introduced difference method into the LoeGNSS constellation to remove the satellite clock error. We took the i-th GNSS pseudorange observation and j-th LEO Doppler observation as an example:

$$\begin{array}{l}\left\{\begin{array}{l}{\rho }_{u,G}^{\left(i\right)}=r_{u,G}^{\left(i\right)}+I_{u,G}^{\left(i\right)}+T_{u,G}^{\left(i\right)}+t_u-t_G{}^{\left(i\right)}+{\epsilon }_{u,G}^{\left(i\right)};{\epsilon }_{u,G}^{\left(i\right)}\sim N\left(0,{\sigma }_{\rho }^2\right)\\ {\rho }_{r,G}^{\left(i\right)}=r_{r,G}^{\left(i\right)}+I_{r,G}^{\left(i\right)}+T_{r,G}^{\left(i\right)}+t_r-t_G{}^{\left(i\right)}+{\epsilon }_{r,G}^{\left(i\right)};{\epsilon }_{r,G}^{\left(i\right)}\sim N\left(0,{\sigma }_{\rho }^2\right)\end{array}\right.\\ \left\{\begin{array}{l}{\dot{\rho }}_{u,L}^{\left(j\right)}={\dot{r}}_{u,L}^{\left(j\right)}+{\dot{I}}_{u,L}^{\left(j\right)}+{\dot{T}}_{u,L}^{\left(j\right)}+f_u-f_L{}^{\left(j\right)}+{\epsilon }_{u,L}^{\left(j\right)};{\epsilon }_{u,L}^{\left(j\right)}\sim N\left(0,{\sigma }_{fd}^2\right)\\ {\dot{\rho }}_{r,L}^{\left(j\right)}={\dot{r}}_{r,L}^{\left(j\right)}+{\dot{I}}_{r,L}^{\left(j\right)}+{\dot{T}}_{r,L}^{\left(j\right)}+f_r-f_L{}^{\left(j\right)}+{\epsilon }_{r,L}^{\left(j\right)};{\epsilon }_{r,L}^{\left(j\right)}\sim N\left(0,{\sigma }_{fd}^2\right)\end{array}\right.\end{array}$$

(9)

where the subscript $L$ and $G$ refer to LEO and GNSS, respectively. i and j refer to the given satellites. The subscripts u and r correspond to the user and reference station, respectively. $\rho$ and $\dot{\rho }$ indicate the measured pseudorange and Doppler, respectively. $r$ and $\dot{r}$ represent the geometric range and range rate between the satellite and the receiver, respectively. $I$, $T$, and $\dot{I}$, $\dot{T}$ denote the ionospheric and tropospheric delay and their change rate, respectively. $t_{*}$ and $t_G$ are the clock bias between the receiver and the GNSS satellite. $f_{*}$ and $f_L$ are the clock drifts between the receiver and the LEO satellite. ${\epsilon }_{r,G}^{}$,${\epsilon }_{u,G}^{}$ and ${\epsilon }_{r,L}^{}$,${\epsilon }_{u,L}^{}$ are the corresponding GNSS pseudorange measurement error and LEO Doppler measurement error at the reference and user station, which obey the $N\left(0,{\sigma }_{\rho }^2\right)$ and $N\left(0,{\sigma }_{fd}^2\right)$, respectively. By differencing the corresponding terms, we obtain:

$$\begin{array}{l}{\rho }_{u,G}^{\left(i\right)}-{\rho }_{r,G}^{\left(i\right)}=\left(r_{u,G}^{\left(i\right)}-r_{r,G}^{\left(i\right)}\right)+\left(I_{u,G}^{\left(i\right)}-I_{r,G}^{\left(i\right)}\right)+\left(T_{u,G}^{\left(i\right)}-T_{r,G}^{\left(i\right)}\right)+\left(t_u-t_r\right)+\left({\epsilon }_{u,G}^{\left(i\right)}-{\epsilon }_{r,G}^{\left(i\right)}\right)\\ {\rho }_{ur,G}^{\left(i\right)}=r_{ur,G}^{\left(i\right)}+t_{ur}+{\epsilon }_{ur,G}^{\left(i\right)};{\epsilon }_{ur,G}^{\left(i\right)}\sim N\left(0,2{\sigma }_{\rho }^2\right)\\ {\dot{\rho }}_{u,L}^{\left(j\right)}-{\dot{\rho }}_{r,L}^{\left(j\right)}=\left({\dot{r}}_{u,L}^{\left(j\right)}-r_{r,L}^{\left(j\right)}\right)+\left({\dot{I}}_{u,L}^{\left(j\right)}-{\dot{I}}_{r,L}^{\left(j\right)}\right)+\left({\dot{T}}_{u,L}^{\left(j\right)}-{\dot{T}}_{r,L}^{\left(j\right)}\right)+\left(f_u-f_r\right)+\left({\epsilon }_{u,L}^{\left(j\right)}-{\epsilon }_{r,L}^{\left(j\right)}\right)\\ {\dot{\rho }}_{ur,L}^{\left(j\right)}={\dot{r}}_{ur,L}^{\left(j\right)}+f_{ur}+{\epsilon }_{ur,L}^{\left(j\right)};{\epsilon }_{ur,L}^{\left(j\right)}\sim N\left(0,2{\sigma }_{fd}^2\right)\end{array}$$

(10)

Assuming that the line-of-sight vectors from the user and the reference station to the satellite were parallel to each other. $r_{ur,G}^{\left(i\right)}$ is approximated as the projection of the baseline vector $b_{ur}$ on the normalized line-of-sight vector $E_{r,G}{}^{\left(i\right)}$ at the reference. Finally, the relationship between the bur and the SD measurements of the LoeGNSS constellation is obtained:

$$\begin{array}{l}{\rho }_{ur,G}^{\left(i\right)}=E_{r,G}{}^{\left(i\right)}\cdot b_{ur}+t_{ur}+{\epsilon }_{ur,G}^{\left(i\right)};{\epsilon }_{ur,G}^{\left(i\right)}\sim N\left(0,2{\sigma }_{\rho }^2\right)\\ {\dot{\rho }}_{ur,L}^{\left(j\right)}=\frac{\partial E_{r,L}{}^{\left(j\right)}}{\partial t}\cdot b_{ur}+E_{r,L}{}^{\left(j\right)}\cdot {\dot{b}}_{ur}+f_{ur}+{\epsilon }_{ur,L}^{\left(j\right)};{\epsilon }_{ur,L}^{\left(j\right)}\sim N\left(0,2{\sigma }_{fd}^2\right)\end{array}$$

(11)

Using inter-station difference, the unknown change rate of atmospheric delay and satellite clock drift in LEO measurements were removed.

2.2.4. Effect of LEO Orbit Error on the Measurement Model

The GNSS satellite state can be calculated using orbit interpolation with the real-time ephemeris issued by IGS, and its orbit accuracy is meter level [2,30]. The Kepler parameters of LEO are available through TLE files. Nevertheless, according to [12], satellite orbit errors calculated using TLE can accumulate to the order of kilometers. Therefore, it is necessary to consider the influence of ephemeris error on differential positioning accuracy.

According to (10) we have:

$$r_{ur,G}^{\left(i\right)}=‖x_G^{\left(i\right)}-x_u‖-‖x_G^{\left(i\right)}-x_r‖$$

(12)

$x_G^{\left(i\right)}$ is the position vector of the i-th GNSS satellite. $x_u$ and $x_r$ are the position vectors of the user and the reference, respectively. Assume that there is a perturbation error $\Delta x$ in the satellite position obtained from the extrapolation of IGS ephemeris parameters, namely $x_G^{\left(i\right)}=x_G^{\left(i0\right)}+\Delta x$, where $x_G^{\left(i0\right)}$ is the true satellite position. Let $g_u\left(x_G^{\left(i\right)}\right)=‖x_G^{\left(i\right)}-x_u‖$ be first-order Taylor expansion at point $x_G^{\left(i0\right)}$:

$$\begin{array}{ccc}\hfill g_u\left(x_G^{\left(i\right)}\right)& =& g_u\left(x_G^{\left(i0\right)}\right)+{\left.\frac{\partial g_u\left(x_G^{\left(i\right)}\right)}{\partial x_G^{\left(i\right)}}\right|}_{x_G^{\left(i0\right)}}\cdot \Delta x\hfill \\ & =& ‖x_G^{\left(i0\right)}-x_u‖+E_u^{x^{\left(i0\right)}}\cdot \Delta x\hfill \end{array}$$

(13)

Likewise $g_r\left(x_G^{\left(i\right)}\right)=‖x_G^{\left(i\right)}-x_r‖=‖x_G^{\left(i0\right)}-x_u‖+E_r^{x^{\left(i0\right)}}\cdot \Delta x$, and:

$$\begin{array}{ccc}\hfill r_{ur,G}^{\left(i\right)}& =& ‖x_G^{\left(i\right)}-x_u‖-‖x_G^{\left(i\right)}-x_r‖\hfill \\ & =& ‖x_G^{\left(i0\right)}-x_u‖-‖x_G^{\left(i0\right)}-x_r‖+\left(E_u^{x^{\left(i0\right)}}\cdot \Delta x-E_r^{x^{\left(i0\right)}}\cdot \Delta x\right)\hfill \end{array}$$

(14)

That is, when the GNSS satellite has an orbit error $\Delta x$, it will cause a deviation in the pseudorange differential measurement $\left(E_u^{x^{\left(i0\right)}}\cdot \Delta x-E_r^{x^{\left(i0\right)}}\cdot \Delta x\right)$.

Similarly the first-order Taylor expansion of ${\dot{g}}_u\left(x_L^{\left(j\right)}\right)$ at point $x_L^{\left(j0\right)}$ is:

$${\dot{g}}_u\left(x_L^{\left(j\right)}\right)={\dot{r}}_{u,L}^{\left(j0\right)}+\frac{\partial E_u^{x^{\left(j0\right)}}}{\partial t}\cdot \Delta x+E_u^{x^{\left(j0\right)}}\cdot \Delta \dot{x}$$

(15)

thus

$$\begin{array}{ccc}\hfill {\dot{r}}_{ur,L}^{\left(j\right)}& =& \frac{\partial ‖x_L^{\left(j\right)}-x_u‖}{\partial t}-\frac{\partial ‖x_L^{\left(j\right)}-x_r‖}{\partial t}\hfill \\ & =& {\dot{r}}_{ur,L}^{\left(j0\right)}+\left(\frac{\partial E_u^{x^{\left(j0\right)}}}{\partial t}-\frac{\partial E_r^{x^{\left(j0\right)}}}{\partial t}\right)\cdot \Delta x+\left(E_u^{x^{\left(j0\right)}}-E_r^{x^{\left(j0\right)}}\right)\cdot \Delta \dot{x}\hfill \end{array}$$

(16)

Namely in the LEO Doppler difference model, the satellite orbit error $\Delta x$ will not only cause the measurement deviation of $\left(\frac{\partial E_u^{x^{\left(i\right)}}}{\partial t}-\frac{\partial E_r^{x^{\left(i\right)}}}{\partial t}\right)\cdot \Delta x$, but also the corresponding velocity error $\Delta \dot{x}$ will cause $\left(E_u^{x^{\left(i\right)}}-E_r^{x^{\left(i\right)}}\right)\cdot \Delta \dot{x}$. In particular, in short baseline (for LEO constellation the short baseline range is no more than 1 km), $E_u^{x^{\left(i\right)}}\approx E_r^{x^{\left(i\right)}}$, $\frac{\partial E_u^{x^{\left(i\right)}}}{\partial t}\approx \frac{\partial E_r^{x^{\left(i\right)}}}{\partial t}$, the orbital error is removed.

2.3. D-LoeGNSS Framework Double-Difference Measurement Model

Although SD can eliminate common errors with spatial correlation, receiver clock bias and clock drift introduced by GNSS and LEO still exist. To this end, a common method is DD. Without loss of generality, the first satellite of each system is taken as the reference satellite. Assuming that the user is stationary, the LoeGNSS DD measurements are yielded:

$$\begin{array}{l}{\rho }_{ur,G}^{\left(i1\right)}={\rho }_{ur,G}^{\left(i\right)}-{\rho }_{ur,G}^{\left(1\right)}=\left(E_{r,G}{}^{\left(i\right)}-E_{r,G}{}^{\left(1\right)}\right)\cdot b_{ur}+{\epsilon }_{ur,G}^{\left(i1\right)}\\ {\dot{\rho }}_{ur,L}^{\left(j1\right)}={\dot{\rho }}_{ur,L}^{\left(j\right)}-{\dot{\rho }}_{ur,L}^{\left(1\right)}=\left(\frac{\partial E_{r,L}{}^{\left(j\right)}}{\partial t}-\frac{\partial E_{r,L}{}^{\left(1\right)}}{\partial t}\right)\cdot b_{ur}+{\epsilon }_{ur,L}^{\left(j1\right)}\end{array}$$

(17)

In (17), DD removes the receiver clock term at the cost of sacrificing one measurement equation for each system. We define the DD permutation matrices:

$$H_G=\left[-e_{M-1},I_{M-1}\right];H_L=\left[-e_{N-1},I_{N-1}\right]$$

(18)

$e_{M-1}$ and $e_{N-1}$ are all 1-element column vectors of order $\left(M-1\right)$ and $\left(N-1\right)$, respectively. $I_{M-1}$ and $I_{N-1}$ are $\left(M-1\right)\times \left(M-1\right)$ and $\left(N-1\right)\times \left(N-1\right)$ identity matrices. Then:

$$\begin{array}{l}\left[\begin{array}{l}{\rho }_{ur,G}^{\left(21\right)}\\ {\rho }_{ur,G}^{\left(31\right)}\\ ⋮\\ {\rho }_{ur,G}^{\left(M1\right)}\end{array}\right]={\left[\begin{array}{cc}\begin{array}{c}-1\\ -1\\ ⋮\\ -1\end{array}& \left|\begin{array}{cccc}1& & & \\ & 1& & \\ & & \ddots & \\ & & & 1\end{array}\right.\end{array}\right]}_{\left(M-1\right)\times M}\times \left[\begin{array}{l}{\rho }_{ur,G}^{\left(1\right)}\\ {\rho }_{ur,G}^{\left(2\right)}\\ ⋮\\ {\rho }_{ur,G}^{\left(M\right)}\end{array}\right]\\ \left[\begin{array}{l}{\dot{\rho }}_{ur,L}^{\left(21\right)}\\ {\dot{\rho }}_{ur,L}^{\left(31\right)}\\ ⋮\\ {\dot{\rho }}_{ur,L}^{\left(N1\right)}\end{array}\right]={\left[\begin{array}{cc}\begin{array}{c}-1\\ -1\\ ⋮\\ -1\end{array}& \left|\begin{array}{cccc}1& & & \\ & 1& & \\ & & \ddots & \\ & & & 1\end{array}\right.\end{array}\right]}_{\left(N-1\right)\times N}\times \left[\begin{array}{l}{\dot{\rho }}_{ur,L}^{\left(1\right)}\\ {\dot{\rho }}_{ur,L}^{\left(2\right)}\\ ⋮\\ {\dot{\rho }}_{ur,L}^{\left(N\right)}\end{array}\right]\end{array}$$

(19)

We define the SD noise matrices ${\mathsf{\epsilon }}_{ur,G}^S={\left[{\epsilon }_{ur,G}^{\left(1\right)}\cdots {\epsilon }_{ur,G}^{\left(M\right)}\right]}^{\mathrm{T}}$ and ${\mathsf{\epsilon }}_{ur,L}^S={\left[{\epsilon }_{ur,L}^{\left(1\right)}\cdots {\epsilon }_{ur,L}^{\left(N\right)}\right]}^{\mathrm{T}}$. The DD noise matrices ${\mathsf{\epsilon }}_{ur,G}^D={\left[{\epsilon }_{ur,G}^{\left(21\right)}\cdots {\epsilon }_{ur,G}^{\left(M1\right)}\right]}^{\mathrm{T}}$, ${\mathsf{\epsilon }}_{ur,L}^D={\left[{\epsilon }_{ur,L}^{\left(21\right)}\cdots {\epsilon }_{ur,L}^{\left(N1\right)}\right]}^{\mathrm{T}}$. According to the law of propagation of variance:

$$\begin{array}{l}Q\left({\mathsf{\epsilon }}_{ur,G}^D\right)=H_G\times Q\left({\mathsf{\epsilon }}_{ur,G}^S\right)\times H_G^T=2{\sigma }_{\rho }^2\times D_{\left(M-1\right)\times \left(M-1\right)}\\ Q\left({\mathsf{\epsilon }}_{ur,L}^D\right)=H_L\times Q\left({\mathsf{\epsilon }}_{ur,L}^S\right)\times H_L^T=2{\sigma }_{fd}^2\times D_{\left(N-1\right)\times \left(N-1\right)}\end{array}$$

(20)

where $D_M=H\times H^T={\left[\begin{array}{ccccc}2& 1& 1& \cdots & 1\\ 1& 2& 1& \cdots & 1\\ 1& 1& \ddots & \ddots & ⋮\\ ⋮& ⋮& \ddots & 2& 1\\ 1& 1& \cdots & 1& 2\end{array}\right]}_{M\times M}$ is the covariance coefficient matrix. Compared to SD, DD further amplifies noise while introducing a correlation between noise.

2.4. Householder-Based Differential LoeGNSS Measurement Model (HB-DLG)

To remove the receiver clock term and circumvent the problem of prominent noise. We introduced the householder transformation matrix into the SD, and performed QR decomposition on the coefficient matrix of the receiver clock item. Assuming the user is stationary, according to (11), the matrix form of the LoeGNSS single difference equation is:

$$\begin{array}{l}{\mathsf{\rho }}_{ur,G}^{}=E_{r,G}^{}\cdot b_{ur}+e_M\cdot t_{ur}+{\mathsf{\epsilon }}_{ur,G}^S;{\mathsf{\epsilon }}_{ur,G}^S\sim N\left(0,2{\sigma }_{\rho }^2{I}_M\right)\\ {\dot{\mathsf{\rho }}}_{ur,L}^{}=\frac{\partial E_{r,L}^{}}{\partial t}\cdot b_{ur}+e_N\cdot f_{ur}+{\mathsf{\epsilon }}_{ur,L}^S;{\mathsf{\epsilon }}_{ur,L}^S\sim N\left(0,2{\sigma }_{fd}^2{I}_N\right)\end{array}$$

(21)

Constructing householder unitary matrices $P^G\in R^{M\times M}$ and $P^L\in R^{N\times N}$, respectively. Such that:

$$\begin{array}{l}{\mathit{P}}^G{\mathit{e}}_M=‖{\mathit{e}}_M‖{\mathit{e}}_{1,0}^G;{\mathit{P}}^G={\mathit{I}}_M-\frac{2{\mathit{u}}^G{\mathit{u}}^G{}^T}{{\mathit{u}}^G{}^T{\mathit{u}}^G};{\mathit{u}}^G={\mathit{e}}_M-‖{\mathit{e}}_M‖{\mathit{e}}_{1,0}^G\\ {\mathit{P}}^L{\mathit{e}}_N=‖{\mathit{e}}_N‖{\mathit{e}}_{1,0}^L;{\mathit{P}}^L={\mathit{I}}_N-\frac{2{\mathit{u}}^L{\mathit{u}}^L{}^T}{{\mathit{u}}^L{}^T{\mathit{u}}^L};{\mathit{u}}^L={\mathit{e}}_N-‖{\mathit{e}}_N‖{\mathit{e}}_{1,0}^L\end{array}$$

(22)

let ${\mathit{e}}_{1,0}^G={\left[\begin{array}{cc}1& \underset{M-1}{\underbrace{\begin{array}{ccc}0& \cdots & 0\end{array}}}\end{array}\right]}^T$, ${\mathit{e}}_{1,0}^L={\left[\begin{array}{cc}1& \underset{N-1}{\underbrace{\begin{array}{ccc}0& \cdots & 0\end{array}}}\end{array}\right]}^T$, then

$$\begin{array}{c}{\mathit{P}}^G={\left[\begin{array}{cc}\frac{1}{\sqrt{M}}& \frac{{\mathit{e}}_{M-1}^T}{\sqrt{M}}\\ \frac{{\mathit{e}}_{M-1}}{\sqrt{M}}& {\mathit{I}}_{M-1}-\frac{{\mathit{e}}_{M-1}{\mathit{e}}_{M-1}^T}{M-\sqrt{M}}\end{array}\right]}_{M\times M}=\left[\begin{array}{l}{\mathit{p}}_1^G\\ \overline{{\mathit{P}}^G}\end{array}\right]\\ {\mathit{P}}^L={\left[\begin{array}{cc}\frac{1}{\sqrt{N}}& \frac{{\mathit{e}}_{N-1}^T}{\sqrt{N}}\\ \frac{{\mathit{e}}_{N-1}}{\sqrt{N}}& {\mathit{I}}_{N-1}-\frac{{\mathit{e}}_{N-1}{\mathit{e}}_{N-1}^T}{N-\sqrt{N}}\end{array}\right]}_{N\times N}=\left[\begin{array}{l}{\mathit{p}}_1^L\\ \overline{{\mathit{P}}^L}\end{array}\right]\end{array}$$

(23)

Dividing the P matrix into two parts, where ${\mathit{p}}_1^G={\left[\begin{array}{cc}\frac{1}{\sqrt{M}}& \frac{{\mathit{e}}_{M-1}^T}{\sqrt{M}}\end{array}\right]}_{1\times M}$ is the first row of the matrix P, and $\overline{{\mathit{P}}^G}={\left[\begin{array}{cc}\frac{{\mathit{e}}_{M-1}}{\sqrt{M}}& {\mathit{I}}_{M-1}-\frac{{\mathit{e}}_{M-1}{\mathit{e}}_{M-1}^T}{M-\sqrt{M}}\end{array}\right]}_{\left(M-1\right)\times M}$. Similarly, ${\mathit{p}}_1^L={\left[\begin{array}{cc}\frac{1}{\sqrt{N}}& \frac{{\mathit{e}}_{N-1}^T}{\sqrt{N}}\end{array}\right]}_{1\times N}$, $\overline{{\mathit{P}}^L}={\left[\begin{array}{cc}\frac{{\mathit{e}}_{N-1}}{\sqrt{N}}& {\mathit{I}}_{N-1}-\frac{{\mathit{e}}_{N-1}{\mathit{e}}_{N-1}^T}{N-\sqrt{N}}\end{array}\right]}_{\left(N-1\right)\times N}$. Left by multiplying ${\mathit{P}}^G$ and ${\mathit{P}}^L$ by (21) gives:

$$\begin{array}{l}\left[\begin{array}{l}{p}_1^G{\mathsf{\rho }}_{ur,G}^{}\\ \overline{P^G}{\mathsf{\rho }}_{ur,G}^{}\end{array}\right]=\left[\begin{array}{l}{p}_1^G{E}_{r,G}^{}\\ \overline{P^G}{E}_{r,G}^{}\end{array}\right]\cdot b_{ur}+\left[\begin{array}{c}\sqrt{M}\\ 0_{M-1}\end{array}\right]\cdot t_{ur}+\left[\begin{array}{l}{p}_1^G{\mathsf{\epsilon }}_{ur,G}^S\\ \overline{P^G}{\mathsf{\epsilon }}_{ur,G}^S\end{array}\right]\\ \left[\begin{array}{l}{p}_1^L{\dot{\mathsf{\rho }}}_{ur,L}^{}\\ \overline{{\mathsf{P}}^L}{\dot{\rho }}_{ur,L}^{}\end{array}\right]=\left[\begin{array}{l}{p}_1^L\frac{\partial E_{r,L}^{}}{\partial t}\\ \overline{P^L}\frac{\partial E_{r,L}^{}}{\partial t}\end{array}\right]\cdot b_{ur}+\left[\begin{array}{c}\sqrt{N}\\ 0_{N-1}\end{array}\right]\cdot f_{ur}+\left[\begin{array}{l}{p}_1^L{\mathsf{\epsilon }}_{ur,L}^S\\ \overline{P^L}{\mathsf{\epsilon }}_{ur,L}^S\end{array}\right]\end{array}$$

(24)

In (24), the parts containing $t_{ur}$ and $f_{ur}$ are both restricted to the first row of the respective measurement equations. We extracted the part where the clock coefficient was 0 respectively, then:

$$\begin{array}{l}\left[\overline{P^G}{\mathsf{\rho }}_{ur,G}^{}\right]=\left[\overline{P^G}{E}_{r,G}^{}\right]\cdot b_{ur}+\left[\overline{P^G}{\mathsf{\epsilon }}_{ur,G}^S\right]\\ \left[\overline{P^L}{\dot{\mathsf{\rho }}}_{ur,L}^{}\right]=\left[\overline{P^L}\frac{\partial E_{r,L}^{}}{\partial t}\right]\cdot b_{ur}+\left[\overline{P^L}{\mathsf{\epsilon }}_{ur,L}^S\right]\end{array}$$

(25)

In the model above, only the unknown baseline vector $b_{ur}$ was included in the transformed equation. To obtain a unique least-square solution requires that:

$$\begin{array}{l}\left(M-1\right)+\left(N-1\right)\ge 3\\ \Rightarrow M+N\ge 5\end{array}$$

(26)

In addition, according to the unitary matrix property, its noise covariance matrix is:

$$\begin{array}{c}Q\left(\overline{{\mathit{P}}^G}{\mathit{\epsilon }}_{ur,G}^S\right)=\overline{{\mathit{P}}^G}\times E\left({\mathit{\epsilon }}_{ur,G}^S\times {\mathit{\epsilon }}_{ur,G}^S{}^T\right)\times {\overline{{\mathit{P}}^G}}^T\\ =2{\sigma }_{\rho }^2\overline{{\mathit{P}}^G}\times {\overline{{\mathit{P}}^G}}^T=2{\sigma }_{\rho }^2{\mathit{I}}_{M-1}\\ Q\left(\overline{{\mathit{P}}^L}{\mathit{\epsilon }}_{ur,L}^S\right)=\overline{{\mathit{P}}^L}\times E\left({\mathit{\epsilon }}_{ur,L}^S\times {\mathit{\epsilon }}_{ur,L}^S{}^T\right)\times {\overline{{\mathit{P}}^L}}^T\\ =2{\sigma }_{fd}^2\overline{{\mathit{P}}^L}\times {\overline{{\mathit{P}}^L}}^T=2{\sigma }_{fd}^2{\mathit{I}}_{N-1}\end{array}$$

(27)

Comparing (20) and (27), the HB-DLG algorithm can remove the common clock term without changing the statistical properties of the measurement noise.

Note that in (25), $\frac{\partial E_{r,L}^{\left(j\right)}}{\partial t}=\frac{x_L^{\left(j\right)}-x_r}{‖x_L^{\left(j\right)}-x_r‖}\times {\left(\frac{x_L^{\left(j\right)}-x_r}{‖x_L^{\left(j\right)}-x_r‖}\times \frac{v_L^{\left(j\right)}}{‖x_L^{\left(j\right)}-x_r‖}\right)}^T$, $E_{r,G}^{\left(i\right)}=\frac{x_G^{\left(i\right)}-x_r}{‖x_G^{\left(i\right)}-x_r‖}$, which means the Doppler Jacobian matrix has been normalized twice more than the pseudorange, resulting in the condition number of the combined Jacobian matrix $\left[\begin{array}{l}\overline{P^G}{E}_{r,G}^{}\\ \overline{P^L}\frac{\partial E_{r,L}^{}}{\partial t}\end{array}\right]$ exceeding $e^5$ and exhibiting ill-conditioned. To increase the numerical stability, we can multiply (25) by the scaling factor $\left[\begin{array}{cc}{I}_{M-1}& \\ & e^3{I}_{N-1}\end{array}\right]$, making the elements in the Jacobian matrix in the same order of magnitude, i.e.:

$$\left[\begin{array}{l}\overline{P^G}{\mathsf{\rho }}_{ur,G}^S\\ e^3\overline{P^L}{\dot{\mathsf{\rho }}}_{ur,L}^S\end{array}\right]=\left[\begin{array}{l}\overline{P^G}{E}_{r,G}\\ e^3\overline{P^L}\frac{\partial E_{r,L}}{\partial t}\end{array}\right]b_{ur}+\left[\begin{array}{l}\overline{P^G}{\mathsf{\epsilon }}_{ur,G}^S\\ e^3\overline{P^L}{\mathsf{\epsilon }}_{ur,L}^S\end{array}\right];\begin{array}{c}\overline{P^G}{\mathsf{\epsilon }}_{ur,G}^S\sim N\left(0,2{\sigma }_{\rho }{}^2{I}_{M-1}\right)\\ e^3\overline{P^G}{\mathsf{\epsilon }}_{ur,L}^S\sim N\left(0,2e^6{\sigma }_{fd}{}^2{I}_{N-1}\right)\end{array}$$

(28)

2.5. Performance Evaluation of HB-DLG Measurement Model Based on CRLB

In GNSS, by assuming that measurement noise is independent and identically (i.i.d) distributed Gaussian white noise, the geometric dilution of precision (GDOP) can serve as a metric of satellite geometry that is independent of the measurement precision. However, in the D-LoeGNSS framework, pseudorange and Doppler are heterogeneous measurements. Although it is possible to assume that the noise in the pseudorange and Doppler are unbiased and independent after SD, the identically distributed does not hold. Therefore, the same derivation does not apply to the current environment.

For parameter–estimation problems, CRLB establishes a lower bound for the variance of any unbiased estimator, i.e., it is impossible to obtain an unbiased estimator with a variance lower than this bound, and provides a standard for comparing the performance of unbiased estimators. In this paper, the user’s positioning error distribution performance is the CRLB. By deriving the CRLB, the theoretical best performance of positioning accuracy under the current environmental parameters can be obtained.

As mentioned above, the difference operation can eliminate the unknown common errors in measurement. Therefore, the measurement residuals only contain noise, and we can derive the CRLB for HB-DLG. We rewrite (28) as follows:

$$\left[\overline{P}\right]\left[\begin{array}{l}{\mathsf{\rho }}_{ur,G}^S\\ {\dot{\mathsf{\rho }}}_{ur,L}^S\end{array}\right]=\left[\overline{P}\right]\left[\begin{array}{l}{E}_{r,G}\\ \frac{\partial E_{r,L}}{\partial t}\end{array}\right]b_{ur}+\left[\overline{P}\right]\left[\begin{array}{l}{\mathsf{\epsilon }}_{ur,G}^S\\ {\mathsf{\epsilon }}_{ur,L}^S\end{array}\right];\begin{array}{c}{\mathsf{\epsilon }}_{ur,G}^S\sim N\left(0,2{\sigma }_{\rho }{}^2{I}_M\right)\\ {\mathsf{\epsilon }}_{ur,L}^S\sim N\left(0,2{\sigma }_{fd}{}^2{I}_N\right)\end{array}$$

(29)

where $\left[\overline{P}\right]=\left[\begin{array}{cc}{\overline{P^G}}_{\left(M-1\right)\times M}^{}& \\ & e^3{\overline{P^L}}_{\left(N-1\right)\times N}^{}\end{array}\right]$. First deriving the CRLB that does not contain the $\left[\overline{P}\right]$ part in (29):

$$\left[\begin{array}{l}{\mathsf{\rho }}_{ur,G}^S\\ {\dot{\mathsf{\rho }}}_{ur,L}^S\end{array}\right]=\left[\begin{array}{l}{E}_{r,G}\\ \frac{\partial E_{r,L}}{\partial t}\end{array}\right]b_{ur}+\left[\begin{array}{l}{\mathsf{\epsilon }}_{ur,G}^S\\ {\mathsf{\epsilon }}_{ur,L}^S\end{array}\right]$$

(30)

Set the measurement vector as $\widehat{\mathsf{\psi }}={\left[\begin{array}{cc}{\widehat{\mathsf{\rho }}}_{ur,G}^S& {\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^S\end{array}\right]}^T$. Correspondingly, $\mathsf{\psi }\left(b_{ur}\right)={\left[\begin{array}{cc}{\mathsf{\rho }}_{ur,G}^S\left(b_{ur}\right)& {\dot{\mathsf{\rho }}}_{ur,L}^S\left(b_{ur}\right)\end{array}\right]}^T$. The joint probability density function (pdf) is:

$$\begin{array}{ccc}\hfill f\left(\widehat{\mathsf{\psi }},\mathsf{\psi }\left(b_{ur}\right)\right)& =& f\left({\widehat{\mathsf{\rho }}}_{ur,G}^S,{\mathsf{\rho }}_{ur,G}^S\left(b_{ur}\right)\right)f\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^S,{\dot{\mathsf{\rho }}}_{ur,L}^S\left(b_{ur}\right)\right)\hfill \\ & =& \frac{1}{\sqrt{{\left(2\pi \right)}^{\left(M+N\right)}{\left(2{\sigma }^2{}_{\rho }\right)}^M{\left(2{\sigma }^2{}_{fd}\right)}^N}}\hfill \\ & \times & \exp \left(-\frac{{\displaystyle \sum _{i=1}^M{\left({\widehat{\rho }}_{ur,G}^{\left(i\right)}-{\rho }_{ur,G}^{\left(i\right)}\left(b_{ur}\right)\right)}^2}}{2{\sigma }_{\rho }^2}-\frac{{\displaystyle \sum _{j=1}^N{\left({\widehat{\dot{\rho }}}_{ur,L}^{\left(j\right)}-{\dot{\rho }}_{ur,L}^{\left(j\right)}\left(b_{ur}\right)\right)}^2}}{2{\sigma }_{fd}^2}\right)\hfill \\ & =& \frac{1}{{\left(2\pi \right)}^{\frac{\left(M+N\right)}{2}}\det {\left(Q_{ur,G}\right)}^{\frac{1}{2}}\det {\left(Q_{ur,L}\right)}^{\frac{1}{2}}}\hfill \\ & \times & \exp \left(\begin{array}{l}-\frac{1}{2}{\left({\widehat{\mathsf{\rho }}}_{ur,G}^{}-{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,G}^{-1}\left({\widehat{\mathsf{\rho }}}_{ur,G}^{}-{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)\\ -\frac{1}{2}{\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^{}-{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,L}^{-1}\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^{}-{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)\end{array}\right)\hfill \end{array}$$

(31)

then

$$\begin{array}{ccc}\hfill In\left(f\right)& =& -\frac{1}{2}In\left({\left(2\pi \right)}^{\left(M+N\right)}det\left(Q_{ur,G}\right)\det \left(Q_{ur,L}\right)\right)+\hfill \\ & & -\frac{1}{2}{\left({\widehat{\mathsf{\rho }}}_{ur,G}^{}-{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,G}^{-1}\left({\widehat{\mathsf{\rho }}}_{ur,G}^{}-{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)\hfill \\ & & -\frac{1}{2}{\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^{}-{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,L}^{-1}\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^{}-{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)\hfill \end{array}$$

(32)

$$\frac{\partial }{\partial b_{ur}}In\left(f\right)=\left[\begin{array}{ccc}\frac{\partial }{\partial x}In\left(f\right)& \frac{\partial }{\partial y}In\left(f\right)& \frac{\partial }{\partial z}In\left(f\right)\end{array}\right]$$

(33)

and

$$\begin{array}{ccc}\hfill \frac{\partial }{\partial x}In\left(f\right)& =& \frac{\partial }{\partial {\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)}In\left(f\right)\frac{\partial {\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)}{\partial x}+\frac{\partial }{\partial {\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)}In\left(f\right)\frac{\partial {\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)}{\partial x}\hfill \\ & =& \frac{1}{2}\left[\begin{array}{l}{\left(\frac{\partial }{\partial x}{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,G}^{-1}\left({\widehat{\mathsf{\rho }}}_{ur,G}^{}-{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)\\ +{\left({\widehat{\mathsf{\rho }}}_{ur,G}^{}-{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,G}^{-1}\frac{\partial }{\partial x}{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\end{array}\right]\hfill \\ & +& \frac{1}{2}\left[\begin{array}{l}{\left(\frac{\partial }{\partial x}{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,L}^{-1}\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^{}-{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)\\ +{\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^{}-{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,L}^{-1}\frac{\partial }{\partial x}{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\end{array}\right]\hfill \\ & =& {\left(\frac{\partial }{\partial x}{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,G}^{-1}\left({\widehat{\mathsf{\rho }}}_{ur,G}^{}-{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)\hfill \\ & +& {\left(\frac{\partial }{\partial x}{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,L}^{-1}\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^{}-{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill \frac{{\partial }^2}{\partial x\partial x}In\left(f\right)& =& \frac{\partial }{\partial {\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)}\frac{\partial }{\partial x}In\left(f\right)\frac{\partial {\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)}{\partial x}+\frac{\partial }{\partial {\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)}In\left(f\right)\frac{\partial {\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)}{\partial x}\hfill \\ & =& {\left(\frac{\partial }{\partial x}{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,G}^{-1}\left({\widehat{\mathsf{\rho }}}_{ur,G}^{}-{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)\hfill \\ & +& {\left(\frac{\partial }{\partial x}{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,L}^{-1}\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^{}-{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)\hfill \\ & =& \frac{\partial }{\partial {\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)}{\left(\frac{\partial }{\partial x}{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)}^T\frac{\partial {\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)}{\partial x}{Q}_{ur,G}^{-1}\left({\widehat{\mathsf{\rho }}}_{ur,G}^{}-{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)\hfill \\ & -& {\left(\frac{\partial }{\partial x}{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,G}^{-1}\frac{\partial {\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)}{\partial x}\hfill \\ & +& \frac{\partial }{\partial {\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)}{\left(\frac{\partial }{\partial x}{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)}^T\frac{\partial {\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)}{\partial x}{Q}_{ur,L}^{-1}\left({\widehat{\dot{\mathsf{\rho }}}}_{ur,L}^{}-{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)\hfill \\ & -& {\left(\frac{\partial }{\partial x}{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,L}^{-1}\frac{\partial {\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)}{\partial x}\hfill \end{array}$$

(35)

take the negative expectation operator for (35):

$$\begin{array}{ccc}\hfill -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial x\partial x}\right)& =& {\left(\frac{\partial }{\partial x}{\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,G}^{-1}\frac{\partial {\mathsf{\rho }}_{ur,G}^{}\left(b_{ur}\right)}{\partial x}+{\left(\frac{\partial }{\partial x}{\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)\right)}^T{Q}_{ur,L}^{-1}\frac{\partial {\dot{\mathsf{\rho }}}_{ur,L}^{}\left(b_{ur}\right)}{\partial x}\hfill \\ & =& {\left(E_{r,x}\right)}^T{Q}_{ur,G}^{-1}{E}_{r,x}+{\left(\frac{\partial }{\partial t}{E}_{r,x}\right)}^T{Q}_{ur,L}^{-1}\frac{\partial }{\partial t}{E}_{r,x}\hfill \end{array}$$

(36)

Correspondingly:

$$\begin{array}{l}-E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial y\partial y}\right)={\left(E_{r,y}\right)}^T{Q}_{ur,G}^{-1}{E}_{r,x}+{\left(\frac{\partial }{\partial t}{E}_{r,y}\right)}^T{Q}_{ur,L}^{-1}\frac{\partial }{\partial t}{E}_{r,y}\\ -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial z\partial z}\right)={\left(E_{r,z}\right)}^T{Q}_{ur,G}^{-1}{E}_{r,z}+{\left(\frac{\partial }{\partial t}{E}_{r,z}\right)}^T{Q}_{ur,L}^{-1}\frac{\partial }{\partial t}{E}_{r,z}\\ -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial x\partial y}\right)={\left(E_{r,x}\right)}^T{Q}_{ur,G}^{-1}{E}_{r,y}+{\left(\frac{\partial }{\partial t}{E}_{r,x}\right)}^T{Q}_{ur,L}^{-1}\frac{\partial }{\partial t}{E}_{r,y}\\ -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial x\partial z}\right)={\left(E_{r,x}\right)}^T{Q}_{ur,G}^{-1}{E}_{r,z}+{\left(\frac{\partial }{\partial t}{E}_{r,x}\right)}^T{Q}_{ur,L}^{-1}\frac{\partial }{\partial t}{E}_{r,z}\\ -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial y\partial z}\right)={\left(E_{r,y}\right)}^T{Q}_{ur,G}^{-1}{E}_{r,z}+{\left(\frac{\partial }{\partial t}{E}_{r,y}\right)}^T{Q}_{ur,L}^{-1}\frac{\partial }{\partial t}{E}_{r,z}\end{array}$$

(37)

The Fisher information matrix is:

$$\begin{array}{ccc}\hfill F& =& \left[\begin{array}{ccc}-E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial x\partial x}\right)& -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial x\partial y}\right)& -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial x\partial z}\right)\\ -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial y\partial x}\right)& -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial y\partial y}\right)& -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial y\partial z}\right)\\ -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial z\partial x}\right)& -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial z\partial y}\right)& -E\left(\frac{{\partial }^{2}In\left(f\right)}{\partial z\partial z}\right)\end{array}\right]\hfill \\ & =& {\left[\begin{array}{l}{E}_{r,G}\\ \frac{\partial E_{r,L}}{\partial t}\end{array}\right]}^T\left[\begin{array}{cc}{Q}_{ur,G}^{-1}& \\ & Q_{ur,L}^{-1}\end{array}\right]\left[\begin{array}{l}{E}_{r,G}\\ \frac{\partial E_{r,L}}{\partial t}\end{array}\right]\hfill \end{array}$$

(38)

Referring to the Fisher information matrix (38) derived from (30), the Fisher information matrix of (29) is obtained:

$$\begin{array}{ccc}\hfill F& =& {\left(\left[\overline{P}\right]\left[\begin{array}{l}{E}_{r,G}\\ \frac{\partial E_{r,L}}{\partial t}\end{array}\right]\right)}^T{\left(\left[\overline{P}\right]\left[\begin{array}{cc}{Q}_{ur,G}^{}& \\ & Q_{ur,L}^{}\end{array}\right]{\left[\overline{P}\right]}^T\right)}^{-1}\left(\left[\overline{P}\right]\left[\begin{array}{l}{E}_{r,G}\\ \frac{\partial E_{r,L}}{\partial t}\end{array}\right]\right)\hfill \\ & =& {\left[\begin{array}{l}\overline{P^G}{E}_{r,G}\\ e^3\overline{P^L}\frac{\partial E_{r,L}}{\partial t}\end{array}\right]}^T{\left[\begin{array}{cc}\overline{P^G}{Q}_{ur,G}^{}{\overline{P^G}}^T& \\ & e^3\overline{P^L}{Q}_{ur,L}^{}{e}^3{\overline{P^L}}^T\end{array}\right]}^{-1}\left[\begin{array}{l}\overline{P^G}{E}_{r,G}\\ e^3\overline{P^L}\frac{\partial E_{r,L}}{\partial t}\end{array}\right]\hfill \\ & =& {\left[\begin{array}{l}\overline{P^G}{E}_{r,G}\\ \overline{P^L}\frac{\partial E_{r,L}}{\partial t}\end{array}\right]}^T{\left[\begin{array}{cc}{Q}_{ur,G}^{}& \\ & Q_{ur,L}^{}\end{array}\right]}^{-1}\left[\begin{array}{l}\overline{P^G}{E}_{r,G}\\ \overline{P^L}\frac{\partial E_{r,L}}{\partial t}\end{array}\right]\hfill \end{array}$$

(39)

Finally, the positioning error variance of the D-LoeGNSS framework satisfies:

$$\mathrm{var}\left(b\right)\ge F^{-1}={\left({\left[\begin{array}{l}\overline{P^G}{E}_{r,G}\\ \overline{P^L}\frac{\partial E_{r,L}}{\partial t}\end{array}\right]}^T{\left[\begin{array}{cc}2{\sigma }_{\rho }^2{I}_M& \\ & 2{\sigma }_{fd}^2{I}_N\end{array}\right]}^{-1}\left[\begin{array}{l}\overline{P^G}{E}_{r,G}\\ \overline{P^L}\frac{\partial E_{r,L}}{\partial t}\end{array}\right]\right)}^{-1}$$

(40)

Therefore, (40) can be used as an evaluation metric of D-LoeGNSS positioning error. In addition, for a given theta, CRLB can also reflect the spatial geometric distribution of different satellite combinations.

3. Experimental Results

This section presents the results of a short baseline experiment using Orbcomm combined with GPS under the D-LoeGNSS framework. Firstly, the signal system of Orbcomm is introduced, and then the experimental hardware configuration and related results are discussed.

3.1. Orbcomm Signal System Introduction

The main problem in the differential system was how to ensure clock synchronization between stations. Therefore, we achieved inter-station epoch matching by partially decoding the Orbcomm signal and aligning the synchronization packet. We introduce its signaling system below.

According to [34,35], the Orbcomm satellites have a VHF band Subscriber Transmitter (STX) that provides a continuous 4800 bps stream of packetized signal to the Subscriber Communicators (SC).

Specifically, the Orbcomm downlink frequency band is 137–138 MHz, divided into 13 channels. S1-S12 channels are used to communicate with user terminals, and the bandwidth of each channel is 25 kHz. Channel S13 is used to communicate with the gateway station with a bandwidth of 50 kHz. Each Orbcomm satellite broadcasts on its own two specific frequency bands [13,36]. The downlink modulation is Symmetrical Differential Phase Shift Keying (SDPSK). Orbcomm downlink data adopts non-return-to-zero differential coding. Its message format is as follows: each main frame contains 16 minor frames, counting 16 s. Each minor frame is 4800 bits long, counting 1 s, and consists of 50 data packets: one synchronization packet, one null-filled fill packet, and 48 other packets that are either fill, command, poll, or message packets. Each packet is 12 words consisting of a packet-type word followed by 9 data words and 2 Fletcher Checksum (FCS) words. Each word occupies eight bits. Its frame format is shown in Figure 4.

As shown in Figure 4, the first packet of each minor frame is used for word and frame synchronization. Further, the downlink does not continuously transmit command and message data, so it has to maintain synchronism by sending fill packets. These fill packets contain “random” data that repeats every 49 packets. The specific packet-type is shown in

Table 1. Orbcomm Packet-Type and function.

Packet-Type	Title 2
65	Synchronization packet
1A	Message packet
1B	Uplink information packet format
1C	Downlink information packet format
1D	Network Control Center information packet format
1E	Fill packet

.

3.2. Experimental Setup and Result Analysis

3.2.1. Experimental Setup

Our experiments were carried out at Beihang University using two sets of receivers, Receiver A and Receiver B, which were kept stationary. All hardware devices in a set of receivers are shown in Figure 5, including the RTL2823U Orbcomm antenna, Trimble Zephyr 3 GPS antenna, dual-channel signal collector, and portable power supply. The dual-channel signal collector integrates two low-noise amplifiers (LNAs) corresponding to frequencies of 137.50 MHz and 1575.42 MHz, an Orbcomm radio frequency integrated circuit (RFIC), a GPS RFIC, an atomic clock that provides a frequency output of 10 MHz, and the SPECTRUM Dual-channel data acquisition board.

We placed Receiver A near the basketball court as a reference station, and Receiver B in the playground as the user station. Their position benchmark was calibrated by the BDStar Navigation Company OEM628 receiver in advance; the baseline was 255.72 m. The experiment site is shown in Figure 6.

We used the TLE files downloaded from NOARD to obtain the ephemeris parameters of Orbcomm satellites. The file contained orbital information for a total of 58 Orbcomm satellites, i.e., OG1 and OG2. According to [36], not all Orbcomm satellites are available and are only visible at certain periods. Therefore, we collected a total of 300 s of experimental data. During the experiment, the reference and the user simultaneously observed two Orbcomm satellite signals at frequency channels S6–S10 and S9–S11, corresponding to satellites FM113 and FM114. The GNSS receiver tracked 8 GPS satellites, namely, PRN: {2,5,11,13,15,18,20,29}, where {2,5,11,18,29} are stably tracked, {13,15,20} were visible part of the time. GPS satellite states and telemetry information can be obtained by the OEM628 receiver.

Orbcomm uses an SDPSK modulation scheme which needs to be squared to transform into a special BPSK modulation to remove data code. The signal spectrum has symmetrical peak lines on the center frequency on the left and right sides, with an interval of 4800 Hz. By adjusting the position of sampling points, we intercepted the data length of 1 s each time for square operation and then performed FFT transformation. Then, the received signal frequency can be obtained. Figure 7 shows the FFT frequency multiplication result after the S6 frequency band of Orbcomm FM113 is down-converted to 28.1625 MHz.

By correlating with the high-precision local carrier, the acquisition accuracy of the frequency can be improved. For a detailed description of extracting the Orbcomm Doppler, refer to [13]. Figure 8 shows the Doppler of FM113 and FM114 at the user station. Next, to recover the data code, the local carrier compensated by the Doppler was correlated with the capture data to achieve carrier stripping. Afterwards, demodulation was achieved by a low-pass filter, sampling decision, and phase recovery. The authors of [34] discussed details for the specific process. Finally, we obtained the synchronization codes of each block and the Doppler measurement. Figure 9 shows the single difference of Doppler. The two satellites appeared simultaneously for about 150 s.

For the GPS part, we only considered the stable tracking {2,5,11,18,29}, and the corresponding skyplot of Orbcomm/GPS is shown in Figure 10.

3.2.2. Experimental Results and Analysis

Since Orbcomm has two visible satellites, according to (26), the combination of 3 GPS+2 Orbcomm is the minimum requirement to solve (28) in real-time. We set Orbcomm Doppler noise to ${\sigma }_{fd}=0.025\mathrm{Hz}$ and GPS pseudorange noise to ${\sigma }_{\rho }=0.5\mathrm{m}$ as empirical values [30] so the CLRB of each satellite combination could be calculated by (40), which reflects the satellite geometric distribution at different epochs. In addition, to compare the performance of the HB-DLG and DD algorithms, we calculated the baseline errors which refer to the Euclidean distance between the calculated baseline vector $\widehat{b}$_ur and the actual calibration vector $b_{ur}$ for each epoch, as shown in Figure 11.

Note that we took the logarithmic form for ease of observation due to poor observation redundancy and significant fluctuations in positioning errors.

Table 2. RMSE of the baseline error under single-epoch of HB-DLG and DD algorithms.

Satellite Combinations	HB-DLG RMSE (m)	DD RMSE (m)
GPS 2/5/11+2 Orbcomm	140.23	278.26
GPS 2/5/18+2 Orbcomm	100.25	182.30
GPS 2/5/29+2 Orbcomm	3.3329 × 103	4.2103 × 103
GPS 2/11/18+2 Orbcomm	2.9376 × 103	3.5575 × 103
GPS 2/11/29+2 Orbcomm	1.9612 × 103	2.5111 × 103
GPS 2/18/29+2 Orbcomm	746.35	862.10
GPS 5/11/18+2 Orbcomm	865.05	869.30
GPS 5/11/29+2 Orbcomm	264.35	497.47
GPS 5/18/29+2 Orbcomm	1.9211 × 103	2.0416 × 103
GPS11/18/29+2 Orbcomm	87.61	177.90

analyzes the RMSE of the baseline error.

In Table 2, the RMSE of various satellite combinations obtained by the HB-DLG algorithm is superior to DD, which verifies the performance of the HB-DLG algorithm. Although the positioning accuracy of these satellite combinations is poor, its error variation is consistent with the CRLB. As shown in Figure 11, unlike the relatively stable positioning error in GPS, the satellite geometry distribution changes caused by the high-speed movement of LEO satellites result in significant differences in positioning accuracy among different satellite combinations. Moreover, the positioning accuracy also deteriorates even within the same satellite combination. By introducing more LEO satellites to increase observation redundancy and enhance model strength, better results can be achieved.

Fortunately, after the Starlink constellation was networked, the number of visible satellites in a single constellation can reached about 20 to 40. In this paper, we used epoch accumulation to increase the number of visible satellites, i.e., all epochs from the current to the initial are used to solve the baseline vector. The positioning accuracy and corresponding CRLB were compared using the HB-DLG algorithm with only 3 GPS epoch accumulations, only 2 Orbcomm epoch accumulations, and 3 GPS+2 Orbcomm epoch accumulations, as shown in Figure 12.

By comparing the CRLB, the convergence speed and error variance of 3 GPS+2 Orbcomm epoch accumulation were better than those of other schemes. Furthermore, with increasing measurement redundancy, the epoch accumulation scheme did not exhibit positioning degradation caused by the deterioration of satellite distribution.

Specifically, the positioning accuracy of only 2 Orbcomm epoch accumulations converges and remains stable after about 90 s. In contrast, only 3 GPS epoch accumulations have the poorest positioning accuracy due to the slow geometric changes of the satellites and the highly correlated Jacobian matrix elements. The 3 GPS+2 Orbcomm epoch accumulation can converge in about 20 s and maintain the highest accuracy.

Table 3. RMSE of the baseline error under epoch accumulation of HB-DLG algorithms.

Satellite combinations	RMSE (m)	RMSE after 90 s (m)
only 2 Orbcomm	3.2104 × 103	45.47
GPS 2/5/11	414.14	274.22
GPS 2/5/11+2 Orbcomm	20.56	17.15
GPS 2/5/18	1.1163 × 103	553.86
GPS 2/5/18+2 Orbcomm	11.14	9.72
GPS 2/5/29	451.87	210.04
GPS 2/5/29+2 Orbcomm	31.85	22.63
GPS 2/11/18	267.80	142.58
GPS 2/11/18+2 Orbcomm	41.97	26.60
GPS 2/11/29	1.3156 × 103	692.58
GPS 2/18/29+2 Orbcomm	37.18	23.19
GPS 5/11/18	1.3471 × 103	620.65
GPS 5/11/18+2 Orbcomm	8.67	7.27
GPS 5/11/29	175.58	83.20
GPS 5/11/29+2 Orbcomm	15.52	14.54
GPS 5/18/29	2.5203 × 103	1.1244 × 103
GPS 5/18/29+2 Orbcomm	18.01	15.79
GPS 11/18/29	669.69	299.36
GPS 11/18/29+2 Orbcomm	7.97	7.03

presents the RMSE of different satellite combinations under epoch accumulation schemes. In particular, it also includes the RMSE of each satellite combination after 90 s.

According to Table 3, even after 90 s of epoch accumulation, the positioning accuracy of only 2 Orbcomms was still inferior to that of the 3 GPS+2 Orbcomm satellite combination. Therefore, when there are few visible satellites, the LoeGNSS constellation epoch accumulation can be used to achieve rapid positioning convergence and improve positioning accuracy.

4. Conclusions

For the insufficient number of GNSS satellites available in complex environments, we considered the LoeGNSS scheme. Based on this scheme, this paper proposed the D-LoeGNSS framework to eliminate the unknown LEO clock error and atmospheric delay variation rate, and deduced the influence of orbit error on the differential measurement model. To address the noise amplification and correlation in DD, the HB-DLG algorithm was proposed which can suppress noise and remove the receiver clock error by introducing an orthogonal matrix. Since the LEO Doppler and GNSS pseudorange are heterogeneous observations, the DOP becomes unsuitable. Thus, the CRLB of HB-DLG was derived as a new evaluation metric.

The experiment was carried out by combining Orbcomm with GPS at Beihang University. Experimental results showed that in single-epoch positioning we took satellite combinations of 3 GPS+2 Orbcomms as an example. Compared with DD, using the proposed HB-DLG improved the positioning accuracy of the corresponding satellite combinations by 49.6%, 45.0%, 20.8%, 17.4%, 21.9%, 13.4%, 0.4%, 46.7%, 5.9%, and 50.8%, respectively. Moreover, due to the high-speed movement of LEO satellites, in GPS 2/5/29+2 Orbcomm, GPS 2/11/18+2 Orbcomm, GPS 2/11/29+2 Orbcomm, GPS 2/18/29+2 Orbcomm, GPS 5/11/18+2 Orbcomm, and GPS 5/18/29+2 Orbcomm satellite combinations, CRLB reflected the accuracy degradation due to the deterioration of satellite distribution. To improve observation redundancy and the number of visible satellites, the positioning error and CRLB of only three GPSs, only 2 Orbcomms, and 3 GPS+2 Orbcomms by epoch accumulation were verified, respectively. Due to the correlation between epochs, only three GPS combination effects were insignificant, and its positioning accuracy remained on the order of hundreds of meters. Only two Orbcomm combinations can converge after 90 s, and its RMSE was reduced from 3.2104 × 10³ to 45.47 m. The 3 GPS+2 Orbcomm epoch accumulation schemes can converge after 20 s and had the highest positioning accuracy. In particular, the RMSE was 7.03 m in GPS 11/18/29+2 Orbcomm. Meanwhile, there was no phenomenon of decreased positioning accuracy in the epoch accumulation scheme.

Future work will investigate the impact of orbit error on long-baseline positioning and mega-constellation-assisted GPS carrier-phase differential positioning.