AU2020101276A4 - Systems and methods for processing GNSS data streams for determination of hardware and atmosphere-delays - Google Patents

Abstract

Methods and systems for precise determination of hardware, ionosphere and troposphere-delays in the Global Navigational Satellite System (GNSS) code and phase signals from a single reference receiver computing to multiple reference computing modes. The single receiver computing mode resolves single differenced (SD) and undifferenced (UD) integers and computes SD and UD ionosphere-delay and troposphere-delay terms, with two initially given phase uncalibrated hardware delays (UHD) in each set of code and phase signals. Substituting the obtained SD and UD results and initial UHD values into the combined equations, one obtaining the remaining UHD parameters epoch by epoch. Multiple-receiver computing mode resolving the double-differenced integers from the UD integers and UHD delays for consistent integer and UHD solutions between all reference stations and deriving the consistent slant ionosphere- and troposphere-delays. The satellite and receiver specific UHDs being estimated through least-square adjustment among all the combined UHD raw measurements. Fitting or filtering processing for each UHD delay time series is performed individually. The system creating raw or fitted state space representation (SSR) corrections from these precise delay estimates from a single reference receiver or multiple receivers, these SSR corrections supporting both real time kinematic positioning and precise point positioning with ambiguity resolutions at the user-terminals with single, dual and multiple frequency receivers. 9/10 550 SSR corrections for SD Receiver8data 05 ionosphere & streams and precise troposphere delays and GNSS orbits and UHDs from a single clocks reference receiver -800 812 810A tApplysingle-epochSR 822 Apply fitted/predicted SSR corrections 820 c ffor SD UHDs, SD ionosphere & SD ionosphere & troposphere delays troposphere delays to correct user SD observables to calibrate user SD observables RTK method PPP-AR method Code-offset determination & 824. Code-offset determination! calibration with SD code and calibration with SD code and phase measurements phase measurements 816 User-receiver RTK processing User-end PPP-AR/PPP-RTK processing with DD single-, dual- ortriple- 826- with SD single, dual- or triple frequency measurements frequency data FIG. 9

Description

9/10

550 SSR corrections for SD Receiver8data 05 ionosphere & streams and precise troposphere delays and GNSS orbits and UHDs from a single clocks reference receiver -800

812 810A tApplysingle-epochSR 822 Apply fitted/predicted SSR corrections 820 c ffor SD UHDs, SD ionosphere & SD ionosphere & troposphere delays troposphere delays to correct user SD observables to calibrate user SD observables RTK method PPP-AR method

Code-offset determination & 824. Code-offset determination! calibration with SD code and calibration with SD code and phase measurements phase measurements

816

User-receiver RTK processing User-end PPP-AR/PPP-RTK processing with DD single-, dual- ortriple- 826- with SD single, dual- or triple frequency measurements frequency data

FIG. 9

Systems and methods for processing GNSS data streams for determination of hardware and atmosphere-delays

FIELD OF INVENTION

[0001] The present invention relates to generally to processing of Global Navigation Satellite

Systems (GNSS) data streams collected from every reference station and more specifically

relates to determination of state-space representation corrections for satellite specific hardware delays, ionosphere-delays and troposphere-delays in GNSS code and phase signals.

BACKGROUND

[0002] The Global Navigation Satellite Systems (GNSS) generally include the US Global Positioning System (GPS), Russia's GLONASS, Europe's Galileo and China's Beidou system. Each of these

systems transmits code and phase signals at three or more frequencies. In the observation equations for a single receiver tracking S satellites with three code and phase signals, there are a

total of 15S+10 parameters to be determined, including 4 satellite position and time states, 1 troposphere-delay, 1 ionosphere-delay, 6 satellite-specific hardware delays, 3 code offsets and 3

phase integer ambiguity terms for each satellite, 4 receiver position and time states and 6 receiver-specific hardware delays for the receiver. There are also 6S receiver noise and multipath

terms. Each type of state parameters and delays has different characteristics and time variations.

To determine these parameters in real time, there are basically two GNSS computing modes in processing GNSS observation data: network-based computing mode and single-receiver computing

mode. Network-based computing makes use of data streams or data files from a large number of regionally or globally distributed receivers to precisely redetermine GNSS satellite orbits and

receiver coordinates, satellite and receiver clock offsets, differential code biases alongside the global ionosphere-delay model, and uncalibrated hardware delays (UHD) in transmitters and

receivers. To support single-receiver based processing, the redetermined biases of satellite orbits and clocks relative to broadcast ephemeris, UHDs for respective code and phase measurements,

ionosphere-delay models for each line-of-sight and troposphere-delays for each station are grouped into State-Space Representation (SSR) corrections (Wubbena et al. , 2005). The integer ambiguities in phase measurements are determined as intermediary parameters alongside the state and delay parameters. The network-based computing problems are complex and have drawn serious research and development efforts in both academia and industry communities worldwide since early 1990s. Single-receiver based computing can refer to a reference-receiver based computing and a user-receiver based computing. A reference receiver may simply output raw observations to users who may perform baseline processing, real time kinematic (RTK) or differential positioning. Alternatively, as shown in Feng et al., (2013), single-epoch data streams from each GNSS reference receiver station are processed to generate station-based solutions, or reference receiver-specific parameters. These may include precise receiver clock offset, zenith tropospheric delay (ZTD), Differential Code Biases (DCB), ambiguity parameters, vertical ionospheric delays. At the user end, a single-receiver computing may use the single-epoch data streams from a reference station or virtual reference station to determines the states of the receiver and integer ambiguities, such as RTK positioning mentioned above. If the single-receiver computing treats the satellites orbit, clock and hardware delays for a network or reference station as given values and determine the receiver states and integers, the processing is known as precise point positioning (PPP) with ambiguity resolutions (PPP-AR). If all SSR corrections from the reference network are provided to the user receiver as known values, the user positioning processing is known as PPP-RTK.

[0003] As part of the research and development outcomes, various network-based real time GNSS

processing platforms and products have been developed for research and commercial PPP-AR services in the past decade. These may refer to the Centre National d'Etudes Spatiales (CNES)'s

CLK93, (Laurichesse & Privat, 2015), Deutsches Zentrum fOr Luft- und Raumfahrt (DLR)'s CLK22, Wuhan University phasebias (Geng et al.,2019) and Geoscience Australia's ACS services. These

products are openly available for research and purposes and often used by hardware integrators and end users to generate PPP-AR solutions in emerging markets. In parallel with public offerings,

there are commercial PPP-AR services with global coverage, such as Trimble's OmniSTAR, Fugro's Starfix and NanCom's StarFire (NavCom Technology, 2016) focusing on marine and agricultural

applications. Trimble's RTX products provide PPP-AR and PPP-RTK services using dense regional

networks (Chen et al., 2015; Trimble Navigation Limited, 2016). Veripos provides six additional PPP-AR and PPP-RTK products for marine applications (VERIPOS, 2017).

[00004] In general, a GNSS computing problem is to solve a subset set of parameters in the

original code and phase signals, depending on users' applications, with the remaining parameters being constrained, given, or cancelled. The constrained or given parameters in a GNSS equation system play the role of datum, although the datum settings could be arbitrary, to certain extent, and can introduce errors or biases for the derived parameters. Cancellation of biases and errors in

GNSS equations is a very effective way in GNSS computing and can be effectively realised by combination and differencing treatment among the original code and phase measurements. One

example is the double-difference operation to remove satellite and receiver specific clock offsets and hardware delays common to the same code and phase signals in both reference and user

receivers. As a result, RTK processing can determine the coordinates to the precision of centimetres to millimetres through kinematic or filtering processing. Theoretically if all the

cancelled bias and delay terms are redetermined from the reference receiver or network epoch by epoch and applied on the same bias and delays terms in the user equations, the user PPP-RTK

processing can achieve the same outcome as RTK processing. In other words, with the same

reference and user receiver data sets or data streams, RTK and PPP-RTK should give equivalent results. The fact is that the existing PPP-based processing takes longer time to resolve the

ambiguity parameter than the RTK-based processing with the same reference stations. The hypothesis is that theoretically, there must be disputes in the estimation or treatment of the

biases and delays in the network or single receiver-based processing.

[0005] Many recent research and development efforts and commercial systems turn to use UD

and uncombined single-epoch code and phase measurements in their network-based computing problems (Gu et al., 2013; Lou et al., 2016; Zhao et al., 2018). One claimed benefit is to allow full

exploitation of the information contained in each individual observation type and preserves the original measurement accuracy. However, the implications of UD and uncombined models must

be clarified. First, mathematically, differenced observations can be completely equivalent to UD observations, while uncombined and linearly combined observations can be completely equivalent

as well. All information and precision can be preserved unless a treatment is chosen not to. GNSS processing with uncombined and UD models must determine all or a subset of state and delay

parameters altogether. The approach depends on precise modelling of every state and delay term

and precisely given constraints for precise solutions. This could be possible for the satellite states based on well-understood satellite dynamics knowledges and static receiver states. However, for

the ionosphere-delays subject to spatial and temporal variations, even complicated treatments with both deterministic and stochastic models do not necessarily lead to centimetre-level resolution needed for fast ambiguity resolution. In general, mismodelling of one type of delays may affect the estimation of all or other state or delay parameters in the whole network. The challenge is how to model different types of delays individually in time and geographic domains without affecting each other.

[0006] This invention aims to address the above problems through both single-receiver computing and network-based computing procedures. The single-receiver computing approach regenerate

the SSR corrections for satellite-specific uncalibrated hardware delays (UHD) in GNSS code and

phase signals, as well as ionosphere- and troposphere-delays. These SSR components are derived with respect to a set of given precise GNSS orbits, clocks and widelane and narrowlaneUHD delays.

With the SD integer solutions from PPP-AR, the SD ionosphere-delays are directly determined with the integer-fixed phase combination and the SD troposphere-delays with the geometry-based

ionosphere-free phase combination. The remaining UHD delays are derived from geometry-free and ionosphere-free combinations epoch by epoch. The single-receiver kinematic computing

approach is then extended to the network-based computing mode. The SD and Unidifferenced (UD) integers and UHD delays from a network of multiple stations are processed to derive DD

integers, DD code offsets, DD ionosphere-delays and DD troposphere-delays for all the essential baselines. Appending the UD quantities for each receiver-reference satellite, SD from one receiver

and DD quantities from baselines, are mapped to all line-of-sight (LOS) directions. Thus, the single

epoch network-adjusted SSR corrections are obtained to retain the consistence of corrections for network-based RTK services. Fitted or filtered SSR corrections then support network-based PPP-AR

and PPP-RTK services.

SUMMARY OF THE INVENTION

[0007] In one aspect, the invention comprises the system of using combined GNSS measurements from a single receiver to determine the hardware delays and atmosphere delays in code and

phase signals at dual or three frequencies. The hardware refers to satellite signal transmitters and antenna. The atmosphere-delays comprise ionosphere and troposphere-delays. The resulting

delays can form three types of state space representation (SSR) corrections from dual- or triple frequency data streams from a single receiver and supporting both PPP and RTK positioning

services. The method comprises the steps of a. reformatting the original code and phase observables into three groups of linear combinations, including geometry-free ionosphere-free (GFIF) combinations, geometry-free ionosphere-present (GFIP) combinations, and geometry-based ionosphere-free (GBIF) combinations, the combinations comprising only linearly independent combinations for the processing in two cases: o For a triple frequency receiver and in each Line of Sight (LOS), the six independent combinations include 4 GFIF observables, 1 GFIP observable, and 1 GBIF observable.

They are mathematically equivalent to the original 3 code and 3 phase signals. Of the six combinations, there two widelane code- combinations from GFIF group and

three phase-only combinations from GFIF, GFIP and GBIF groups, thus preserving the phase precision of the states and delays to be determined.

o For a dual-frequency receiver in each LOS, the combinations include 2 GFIF, 1 GFIP

and 1 GBIF combinations, there two code-dominated combinations and two phase only combinations.

b. Making use of three GFIF models to determine Undifferenced (UD) integers for the receiver to-satellite directions and SD integers between satellites. The initial combined UHD

solutions from the previously determined or externally provided sources are applied to reduce the biases in of the float UD and SD integer ambiguities. For each receiver, using a

least-square estimation, the UD float ambiguities can be adjusted with SD integers as constraints and result in UD integers in LOS directions.

c. Giving the known precise GNSS orbits, clock and station coordinates and two initial phase hardware delay for each line of sight from the receiver, determining slant and SD

ionosphere and troposphere-delays with phase-only measurements, and computing UD uncalibrated hardware delays (UHD) in code and phase signals epoch by epoch in all

receiver-satellite directions.

[0008] In one embodiment, the obtained single-epoch time series for UHDs, SD ionosphere-delays,

SD troposphere-delays from a single receiver form a set of single-epoch SSR corrections, which can

be accessed a user receiver to correct the same bias terms in the user measurements. In another embodiment, the obtained single-epoch SSR corrections are fitted as functions of time, which can

then be accessed by a user receiver to correct the same bias terms in the user measurements a. wherein a set of single-epoch SSR corrections are applied to calibrate the SD code and phase measurements, users obtain DD linear equations for original phase measurements for single, dual or triple frequency, so the users can perform single-reference RTK positioning with single frequency, dual-frequency or triple frequency measurements. b. wherein the fitted SSR corrections are applied to correct the same delay terms in the SD code and phase measurements, the users obtain SD linear equations for original code and phase measurements for single, dual or triple frequency signals, so the users can perform single reference PPP-RTK positioning with single-frequency, dual-frequency or triple frequency measurements c. In both cases, the code bias offsets that cannot be cancelled in the treatment can be estimated with long-term or historical data before the receiver is entered the RTK or PPP-RTK services with respect to the reference receiver, thus could be removed from the code signals.

[0009] In another aspect, the invention comprises an additional system of using combined GNSS

measurements from a plurality of receivers to adjust the UHD in LOS code and phase signals at dual or three frequencies. As a result, the resulting SSR corrections over a network of receivers

maintain consistence between stations and can support both PPP and RTK positioning services within a network coverage. The additional system comprising,

a. Performing DD operation on the vector of UHDs and integers for all LOS directions to obtain DD integers for the necessary baselines to enable consistence between all LOS

integers b. Mapping the baseline DD integers and the first station SD integers to obtain consistent LOS

integers by appending the UD integers of each station to the reference satellite c. Performing DD operation onto code UHD measurements and obtaining DD code UHD

offset samples. Taking average or fitting over historical data period and gives precise

estimation of the DD code offsets. d. Removing the obtained LOS integers, LOS ionosphere-delay and LOS troposphere-delays

from the selected UD combined measurements, forming linear equations for all satellite and receiver specific UHDs and performing least-square adjustment to derive the satellite

and receiver specific UHD time sequences of the whole network. e. Converting combined UHDs to original code and phase UHDs and performing fitting

process to each combined or UHD time series for a function of time.

[0010] In one embodiment wherein the obtained single-epoch and adjusted time series for

satellite specific UHDs, slant ionosphere-delays, slant troposphere-delays from a network of receivers are grouped as the single-epoch and network adjusted SSR corrections. The single-epoch

ionosphere-delays and troposphere-delays for all stations and satellites can be fitted by appropriate interpolation models, respectively, such as linear combination models and low order

surface models (Dai, 2001). In another embodiment, the obtained single-epoch SSR corrections are fitted by an appropriate interpolation model, such as polynomial models individually as

functions of time. For user receivers within the network coverage, wherein there are two options:

a. the users perform the ionosphere and troposphere interpolations using the single-epoch

SSR corrections from the surrounding stations and satellite specific UHDs to calibrate the same bias terms in the measurements and proceed with network-based RTK positioning

with single-frequency, dual-frequency or triple frequency measurements. b. the users perform the ionosphere and troposphere interpolations using the single-epoch

SSR corrections from the surrounding stations and fitted satellite-specific UHDs to calibrate the same bias terms in the measurements and proceed with network-based PPP-RTK

positioning with single-frequency, dual-frequency or triple frequency measurements.

[0011] The SSR components have different formations in various code and phase or combinations.

The uncertainty of fitted lumped SSR correction to SD phase observables is in the order of millimetres to a few centimetres, thus theoretically providing the basis for SSR corrections in

supporting instant or fast ambiguity resolution or fast position convergence. The uncertainty of fitted lumped SSR correction to SD code observables can be in the order of a few to several

centimetres, thus theoretically providing the basis for the SSR corrections in supporting decimetre to submetre positioning accuracy with single or combined code measurements.

BRIEF DESCRIPTIONS OF THE DRAWINGS

[0012] In the drawings, like elements are assigned like reference numerals. The drawings are not necessarily to scale, with the emphasis instead placed upon the principles of the present invention.

Additionally, each of the embodiments depicted are but one of several possible arrangements

utilising the fundamental concepts of the present invention. Several aspects of the present invention are illustrated by way of example, and not by way of limitation, in detail in the figures,

wherein

[0013] FIG.1 is a representative illustration of a Global Navigation Satellite System (GNSS) with 4

satellites in view. Each satellite transmits code and phase signals at dual or multiple frequencies.

[0014] FIG.2 is a schematic flowchart showing the overall structure of the sequential GNSS computing system of implementing the present invention.

[0015] FIG.3 is a schematic flowchart showing the major steps of forming combined UD and SD observation equations that implementing the present invention.

[0016] FIG. 4 is a schematic flowchart showing the major steps of the computing UD integers and

SD integers of implementing the present invention

[0017] FIG. 5 is a schematic flowchart showing the major steps of the computing slant ionosphere

and troposphere-delays of implementing the present invention.

[0018] FIG. 6 is a flowchart of a method for determination of UHD in multiple frequency code and phase signals and fitting them as function of time of implementing the present invention.

[0019] FIG. 7 is a flowchart of a method for determination of LOS integers for multiple stations, then recomputing slant ionosphere and troposphere-delays in implementing the present invention.

[0020] FIG. 8 is a flowchart of a network-based computing method for redetermination of satellite

and receiver specific UHD in multiple frequency code and phase signals and fitting them as function of time of implementing the present invention.

[0021] FIG. 9 is a flowchart of a system for estimation of user states and integers with single-epoch and fitted SSR correction messages from a single reference receiver and user-observation data

streams that implementing the present invention.

[0022] FIG. 10 is a flowchart of a system for estimation of user states and integers with SSR

correction messages from a network of reference receivers and user-observation data streams that implementing the present invention.

BRIEF DESCRIPTIONS OF THE TABLES

[0023] In the tables, like equations are assigned like reference numerals. Each of the Tables lists some of possible linear combinations in each type or group utilising the fundamental concepts of the present invention. They are given by way of example, and not by way of limitation, in detail in the tables, after wherein

[0024] Table 1 lists the original code and phase observation equations in three frequencies represented by fi, f2 and f3, and the description of each notation or each term in the equations.

[0025] Table 2 gives the definitions of all combined code and phase signals, observation equations, and descriptions of the notations to be used in implementation of this invention.

DETAILED DESCRIPTION OF VARIOUS EMBODIMENTS

[0026] The invention relates to a system and method for processing GNSS data files and streams

from a single receiver or multiple receivers. When describing the present invention, all terms not defined herein have their common art-recognised meanings. To the extent that the following

description is of a specific embodiment or a particular use of the invention, it is intended to be illustrative only, and not limiting of the claimed invention. The following description is intended to

cover all alternatives, modifications and equivalents that are included in the spirit and scope of the invention, as defined in the appended claims.

[0027] FIG. 1 illustrates a representative view of Global Navigation Satellite System (GNSS) 10. In one embodiment, the GNSS 10 comprises the Global Positioning System (GPS). The GNSS 10

includes a plurality of satellites 11 orbiting the earth with each satellite 11 transmitting signals 12. The GNSS 10 may broadcast the signals 12 on multiple frequencies. For example, if the GNSS 10 is

a GPS system, each satellite 11 can broadcast the signals 12 using two or three frequency (i.e. the

fi, f2 and f3 frequencies used by satellites in the GPS constellation). A GNSS receiver 15 is provided that is operative to receive the signals 12 transmitted by the plurality of satellites 11. A receiver

tracks dual or triple frequency GNSS signal 12 continuously from one more GNSS constellations. The constellations can include any combinations of GPS, Galilleo, Beidou, Glonass, and QZSS.

[0028] FIG. 2 is a schematic illustration of the overall embodiment of the GNSS computing system

of the present invention that processes the data files or data streams 110 for a single GNSS receiver 15 as shown in FIG. 1. The single receiver computing embodiment 100 can include the

preparation method 200 for physical observation equations for selected combined signals as

shown in FIG. 3, Integer ambiguity resolution subsystem 300 as shown in FIG 4, the computing subsystem 400 for determination of slant ionosphere-delays and troposphere-delay as shown in

FIG 5, the computing subsystem 500 for determination of single-epoch satellite-specific UHDs and

their fitting or filtering processing as shown in FIG. 6, the computing subsystem 600 for network based treatment for DD integers and consistent LOS integers as shown in FIG.7; the computing

subsystem 700 for determination of network-adjusted single-epoch UHD observables and fitted/filtered UHD results, the user-end RTK and PPP-RTK processing module 800 with SSR

corrections from a single reference receiver and the user-end RTK and PPP-RTK processing module 900 with SSR corrections from a network of multiple reference receivers. The methods and

computing subsystems are described in detail in the rest of the description.

[0029] FIG.3 is a schematic flowchart of Method 200 comprising a process to establish six linearly

independent combinations over three code and phase observations in each line of sight direction and SD equations with respect to a reference satellite.

[0030] Table 2 lists the key combinations to be used in this invention, including geometry-free and

ionosphere-free (GFIF) combinations 21, 22, 23, 24, geometry-based ionosphere-free (GBIF) combination 25, and geometry-free and ionosphere-preserved (GFIP) combination 27. Such a set

of the six types of combinations is completely equivalent to the set of six original code and phase

observables for each receiver-satellite path.

For the selected combinations, the observation equation system for all types of delays and noise terms is expressed as follows:

P -0#13 T 1 ifl5 0] b1 T 0 0 0 42 0 eT P2 -32 P P2 0 b2 0 0 N 2 q F_ F 0 b3 0 0 + ND #A -P 21 3 8 232 13 #10 0 ' -4B. 0 0 A34 )p2 32 432 -6 AIF 032 02 1 2 2B 2 I A + f-( ) 0 -2 # -#3 j, # ,i #A ]j j, 0) -0, -A 43 0 -4

for s=1,2,...,S (Equation 1)

where the subscript "r" represents a single receiver until otherwise specified.

p = X -X"+CLK,-CLK'+Ts, where T is the slant troposphere-delays to all the satellites,

computed with the same empirical model. For instance, the model adopted by the Space-Based Augmentation Systems (SBAS) standards [RTCA-MOPS, 2006] is as follows:

T - (T + T,, (Equation 2)

where Tz,ry and Tz,wet are calculated from the receiver's height and estimates of five

meteorological parameters: pressure [P(mbar)], temperature [T(K)], water vapour pressure [e (mbar)], temperature "lapse" rate [P(K/m)] and water vapour "lapse rate" [A (dimensionless)]. The obliquity factor is the Black and Eisner mapping function (Black and Eisner,1984):

1.001 g(z,) = - 101CS(: 0.002001+ cos 2 Equation 3) ((z)

The term dT represents a residual troposphere-delay with respect to the computed value from

the above empirical model.

[0031] Step 210 computes the matrix A as follows:

1_ 0 0 0 0 1 0 0 -1 0 0 0 f0 0 010 0 -1 0 0 A 0 0 0

0 0 0 1 , f2 f3 2, = A fl+ A f±+A ffl 0 2 0 0 0 -, 0 0 0 _182t -)622 1 A3_f A = 2

0 0 f 0 0 0 0 A 0 0 f,+A f.+f A AA f fA3 0 0 0 -1 0 1j 0 0 0 0 A _f

0 0 0 1 0 0 j (Equation 4)

where A 1, A 2 and A are invertible. Therefore, once the six UHD delays in the combined signals are determined, the UHDs in the original code and phase signals can be uniquely obtained, while the noise terms in the original signals can be fully propagated to the combined signals

[0032] Comparing to use of the original signals, the combined models in (Equation 1) allow various types of parameters to be estimated sequentially or separately, while preserving all parameters and precision. In the dual frequency case, the second and fourth combinations in Equation 1 do not exist or can be removed.

[0033] Step 220 applies single-difference (SD) between satellites "s" and the reference satellite "1" and linearization, six SD linear observation equations for the position states, SD integer

ambiguities, SD hardware delays, SD ionosphere and troposphere-delays will be obtained.

[0034] The computed SD linear questions for each satellite are expressed as follows:

013 0 0] iPT 0] 0 ~0] 0~| 0 232 0 Vb T eT e P2 P2 0 0 0 0 3) ( 0 VN 1 Vb2 + 2 e2 A( - )- p = dX+ Vd?" + 0 + VN +A -B3 +A( - 2 13NVB

) 0, q5, 0 0 0 0 fl 21-213 6l 22'ql 2 ~~ #2 02 1 H 1 Af3 j0 -VB 2 62 62 031 031 -o -( ) 01 f f -2' 3VB 3j1, -'31, -'31, -2 13 0

for s=2,3,...S (Equation 5)

where the symbol "V" stands for the between-satellite single difference operator, VdT"is the SD

VqS' troposphere-delay of the satellite "s" with respect to satellite "1";. - is the SD ionosphere f2

delay of the satellite s with respect to the satellite 1 on frequency fi; H is a 3-dimensional row

vector for three partial derivatives of the SD range with respect to 3D receiver coordinates; dX is the 3-dimensionsal column vector for the derivations of the receiver coordinates with respect to

the given states; Vp' is the SD computed range with precise satellite orbits and clocks between

sth and 1' satellites from the rth receiver

Vp'" = X,- X" -X,-X 1 -(CLKs -CLK 1 )+(T -T1 ) for s=2, 3,...S (Equation 6)

[0035] In addition, according to the structure of bL and BL in Table 1, the following results are obtained:

1 VbL = -VbL' +Vbj.;, VBL = VBja for L=1,2,3, s=2,3,...,S (Equation 7)

[0036] Step 230 establishes the covariance matrix for the combined measurement error term in

Equation 5. To bring together all the SD elements in one observational equation system for a single station, a single-difference matrix D is set for all satellites in view, where the first satellite is

set to the reference satellite.

-1 1 0 -- 0]

S= 1 0 1 ... ,D=SiIi ,=S 1 ®A (Equation 8)

L-1 0 0 -.. 1j

where 16 is the 6-dimensional identity matrix; the symbol "@ " is the Kronecker product. For a

receiver tracking S satellites, there are a total of 6S original or 6S combined measurements, respectively. The combined noise vectors and original noise vector 6 are related as follows: el 61 1l e2

.2 E2 ,=A e= 3 s=1,...,S, (Equation 9) _ S = _ ,=A c S 2, La Es]C _63i

[0037] Therefore, the covariance matrix of the SD combined models (Equation 4) can be obtained

by variance propagation from the six original code and phase measurements.

Cov(DS)= DCov(6)D =DCov(E)D (Equation 10)

[0038] The equivalence between SD combined equation 4 and their original SD models in terms of receiver states and integer ambiguity parameters is then completely preserved due to use of the

observation vector (Equation 5) and covariance matrix (Equation 10) altogether.

[0039] FIG. 4 describes the system 300 making use of precise orbits and clocks, and initial phase UHD corrections for integer ambiguity resolution.

[0040] In another embodiment, Step 310 may implement the geometry-free approach with (Equation 1) to fix the UD and SD integers. With GFIF models 21 and 22 and 24, the UD extra-wide

lane, wide-lane and narrow-lane floating ambiguities can be obtained the average over their samples of size M over a data window,

0 = [ (Equation 10a)

1[P,(t)- (t)-(b + B(equationl10b) MA7

] = [#,(t)-(p 82(t)+p f,(t))+(Bf A(B +N )-p2(B.+N 2))] (Equation 10c) MA 1 2

where b' b' B' B'2B or the combined terms are initial UHD values derived from previous or historical data sets through Equation 61 or external products for PPP-AR processing such as Centre CNES's CLK93 (Laurichesse & Privat, 2015) and Wuhan University phasebias (Geng et al.,2019). Removing the effects of these initial combined UHD values can reduce biases in the floating ambiguity solutions so that the UD integers can more reliably and stably fixed. N ,N3 2 are the UD

integers by rounding the float UD ambiguity solutions N,3N32. For Nthe rounding integer is

denoted by N.The data window for averaging may be set to be about tens of minutes, so the effective independent sample size M could be over 100. As a result, the noise level of the average values for the narrow lane integer in Equation 10c can be reduced to about 0.1 cycles.

[0041] Single differencing (SD) between satellites can delete the receiver-specific hardware delays

and clocks-bias. Step 320 determines the SD integers by applying SD operation on Equations 10a, b and 10c. That is, SD integers are similarly obtained from SD code and phase measurements.

The SD extra-wide-lane, wide-lane and narrow-lane floating ambiguities can be obtained by the average over their samples of size M over a data window,

VA," [rv(t)-v (t)-(vb +1 vB°)], (Equationh1a) 13MA"7 VR2 - 3 ,=1IN [vP,(t)-V,(t)-(vb2++2 3VB° 2)], (Equation (equation 11b) 11b) -M vN° - -[v 3 (t)+p,,vq1 (t)-(py (t))+(VB 1 0 -p,A(VB'°+vN 0+vo)](Equationl11c) 12A)-p A (VB 2 +vN 2 ))( MAI =

where Vb°,Vb32VB4,VB,VB10 are initial UHD values which may be obtained from the products

of this invention or a well-designed PPP-AR processing platform. The consistence of SD integers

can be verified by stable observed-computed residuals or kinematic position solutions with the SD GBIF model 21 or 26 phase measurements along with initial SD UHD values.

[42] Step 330 performs the least square adjustment to obtain the consistence between each set of UD integers, through the following equation system for each frequency with the wavelength

N0 - N° = dN+ c (Equation 12) SIdN =VN - SIN°

where N0 is a S-dimensional float ambiguity vector, and its initial UD integer vector by rounding is

N0. VN is the known (S-1)-dimensional SD integer vector and Si is their (S-1)-by-S matrix. The UD ambiguity correction term is represented by dN, its least square with constraint conditions is

solved as follows:

1 J(N° - k NO)+S1(S S)(VN -SN°). (Equation 13) dN= k[Is-SIS1 S 1S

[0043] Therefore the LOS integer vector can be obtained

N= No+[dN]rounding. (Equation 14)

It is noted that the second term of Equation 13 in the right-hand side maps a SD vector term onto

the UD vector term, regardless of integer or real-valued nature. The results 340 include 3(S-1) SD integers and S consistent LOS integers for all directions.

[0044] FIG. 5 illustrates the method 400 comprising necessary steps of computing slant

ionosphere and troposphere-delays.

[0045] After the SD integers are correctly determined, Step 410 computes the slant ionosphere delay for the reference satellite path with GFIP 25

- [ -0,+(AB,°-A1 3 BI°3)+(A 1K1 -2 1 31 3)] , (Equation 15a)

and (S-1) SD ionosphere-delays

= f- [V -V , +(1 ,VB, - 13 A VB°3)+( VK- A 13 VK 3)]', (Equation 15b) f3+ f

where s=2,3,...,5,the terms (B- 3B ) and (VB°- 3VB ) represent the initial UHD values,

which may vary with time, or remain a constant. However, they play the roles of the datum in the

determination of the slant and SD ionosphere-delays in Equations 15a and 15b. Meanwhile, these

terms serve as boundary conditions for UHD determination,

(2B1 - 1 3B3 )r =(2B° - 4 3 B°)I (Equation 16a)

(2VBI - 1 3VBB)' =(VL-VB- YB°I),' s=2,3,...,S (Equation 16b)

which remain unchanged for both reference-station-based and user-based computing. While the

external SD UHD products can be used as initial values to set (Equation 16b) and the datum settings for the reference satellite path (Equation 16a) can have more flexible options. One

straightforward option is to set [( 0) B ,-0 B )] =0. An alternativeistoset

[(i,o 0 0) (1,1,0)- B1,)]'with respect to an external ionosphere model. Different settings can

affect the LOS UHD values and variations, but not the positioning results.

[0046] Step 420 computes the slant troposphere-delay for the reference satellite path with the narrow-lane GBIF model 25 in Table 2,

dT = [V#IF Vp+ XJ )V13 +1 X+1 ( "1 3 )B,+2(B"]., (Equation 17a) f1 -,f3 f; -,f 3

and the SD residual troposphere-delays

Vdl7 )=[VIFVp'+ )V 13 + 1V 1+( )VB, +VB,, s=2,3,...,S(Equation 17b) where the term vpa as shown in Equation 2 contains the effects of satellite orbit and clocks errors, and receiver position coordinate errors as well the troposphere model errors. The terms

(" Bf +x1 B") and "I vB4 +kevB° play the roles of the datum in the slant and SD f-fA f -fA troposphere-delays. Meanwhile, these terms serve as boundary conditions for UHD determination as well and the values can depend on the satellite clock offset definition. In the estimation of six

UHDs, it also serves as a boundary condition for UHD estimation, i.e.,

[(J )B1,+k1 B],=[( " )B03+kB/] (Equation18a) f-f3 fl-f3

)=[( "f )VB/3+ VBO]/ (Equation 18b)

where s=2,3,...,S. The datum settings for the reference satellite path (18a) has more options. One

straightforward option is to set [( (1,1,0) )B(1 _ + 1 iB( 0 )] = 0. Another option is to set the il-fz

[( "") )B,_I o+ ()Bo,) ] value in the reference satellite path with respect to a new troposphere

model and receiver "r" and satellite "1" clock offset.

[0047] The condition Equations 16b and 18b allow the SD widelane and narrowlane UHDs for dual frequency signals to remain unchanged. The SD ionospheric delay and tropospheric delay terms

have the precision of serval millimetres.

[0048] Step 430 compute the slant ionosphere-delays by appending the slant delay for reference

satellite path to all SD ionosphere-delays gives the slant ionosphere-delays in all LOS directions, as follows:

1] 0 0 - 0 2] 1 1 0 qr qi+ 1 V0 1 o0 1 ' , 2 (Equation19) -2 f 2 1 ff M ~ - - 11- s,1

[0049] Step 440 computes the slant troposphere-delay for the receiver-reference satellite direction and appends the UD delays to all SD delays to give the slant troposphere-delays in all LOS

directions, with the following equation

1 0 0]VdT2,1 0 1 0 ... 0 1 dj 0 1VdT ',1 d,= 1dT 1 +0 1 - -0 .(Equation 20)

Vd s,1 1 0 0 --- 1]

[0050] FIG. 6 outlines the procedures of Method 500 for determining the UHDs in code and phase measurements from the outputs of Method 300 and Method 400: UD and SD integers and slant ionosphere and slant troposphere-delays.

[0051] Step 510 moves all the knowns or determined UD quantities to the left-hand-side of Equation 1 for the reference satellite path and the SD quantities to the left-hand side of Equation 5, the measurement equations for code and phase UHD parameters for the slant path is obtained as follows:

yoP33 -# 01 0 0 0 b7 e, '~~2~~32 13 0~ 0 0 ~ -b ~ el Y 2n-03 0 f3 0 K-7 0 b2 e2 y3 _ P1F IF -di- N + 0 =A b +A e y4 1- 21 13 - 22 32 0 ' # A3 3 Jt22232 -0 4B, Ys -IF 3 , 3)B2 2

YJr _ 01-03 Jr 0] f-3 - A -ABj ej

(Equation 21a)

For each SD with respect to a reference satellite "1", the measurement equations for code and phase UHD parameters are:

T Vy - V - 01 10 0 Vb eT el Vy2 VP2 -V#0 2 0 _f 0 V T" 0 Vb2 e2 e2 Vy3 _ VPIF VIF 0d- - 0 ' =A + e3 -e Vy4 VNV2V, 0d' p p2 BA +A VY V#IF rV, I 13 .r -VB 62 -2

Vy J, V#I -VA3 j, Jfl -f - - ,VB3 U, _'I J, _'3 1, 0 -A j - AV~ 6] fi

For s=2,3,...,S (Equation 21b)

[0052] Step 520 inversely determines the UHD values in the original SD code and phase signals through the inverse of the matrix A, b1 1 V s1I V (-I )T 1 b, Y, yb Vy, T'(2,) 1 ]S'1 b2 Y2 V 2 Vy 3 + 2r

ZI. b3 =A1 3s,l Vb3_ =A1 3 ,Z, (Z I+VZ 31) AB Y4 -2VB Vy 4 -22B2 Ys -A2V 2 VY 5 --' 3B3jUr _Y 6Ajr B -2,5 A ~ 3 VB3 jr VYy+VZs,1)T 6 jr PIvZ)

(Equation 22)

[0053] Due to the structure of A matrix, three phase UHD terms can be determined with the three

phase-only combinations epoch-by-epoch. But V5 1(t)=VBj° and VB3(t)=VB3 remain

unchanged. Ei(t)=B0 =0 and 63(t)=B3 =0. The precision of single-epoch phase UHD5 2 (t)

and VB2(t)solutions from (Equation 23) is in the same order as the precision of phase

measurements.

[0054] Step 530 performs fitting or filtering processes gives smoothed estimates for SD code UHDs

Vb(t), V 2(t) 3 and phase UHD V52(t), as functions of time. Fitting V52(t) does not V (t)

necessarily result in much higher precision, but allow the UHD values to be predicted or interpolated for real time applications.

[0055] Fitting or filtering processes over a moving time window will give smoothed estimates as

functions of time, Vb(t), Vi(t) Vb(t), which results in much higher precision and allows the

code UHD values to be predicted or interpolated for real time applications.

[0056] Step 540 collects all SD single-epoch UHDs and fitted UHDs for each reference receiver to

form single-epoch SSR corrections (VZ', Vdt."' and ,for s=2,3,...,S) and fitted SSR corrections f2 (VZ", Vdl,' and , for s=2,3,...,S).

[0057] Step 550 collects all LOS single-epoch UHDs and all slant ionosphere and slant troposphere delays and prepare for the system 600. For every epoch, the combined UHDs in the first to sixth

equation can be obtained. For convenience, for each receiver "r", we specify the S-by-1 vector Ir for slant ionosphere-delays, the vector dTr for slant troposphere-delays and the S-by-6 matrix Z, for the LOS UHDs. For all receivers, we also specify the RS-by-6 matrix Z for LOS UHDs and the RS by-1 vectors I and dT for all slant ionosphere and troposphere-delays in the network, respectively: q 1 dT dTl (Z)T Z q 12 dT, 2 dT 2 (Z)T Z I, q ,I= I ;dT= dT3 ,dT= dT ;Zr= (Z)T ,Z (23) T _qS _RI _dTsj dr] (ZS) ] _ZR]

[0058] The precision of all single-epoch delay samples for LOS ionosphere, troposphere and phase

UHDs is all in the order of several millimetres to 1 centimetre, while the single-epoch code UHD

samples contain the effects of code noise and multipath errors, thus having the uncertainty in the order of a few to several decimetres. The precision of the fitted or filtered code UHD over a data

window of tens of minutes will be in the order of a few to several centimetres. The SD formations of the above delays from a single-receiver can be used as part of SSR corrections for RTK and PPP

RTK services.

[0059] FIG. 7 is a schematic illustration of the embodiment of the GNSS computing system 600 of the present invention that processes the data files or data streams 110 for a plurality of GNSS

receivers 15 shown in FIG 1 and the processing outcomes from the single reference receiver

embodiment 100. The system 600 is a network-based computing embodiment comprising double differenced treatments to maintain the consistence between all LOS integers over a network of

multiple stations, thus the consistence between LOS ionosphere and LOS troposphere-delays.

[0060] Step 610 restructures LOS integers for all the LOS directions for a network of R stations and

S satellites, by introducing the following notations:

P,ADi 1] N,] P p © # N2 P= P2 A2 ,0 = .' N .2

PR] _P, U L1)R I~ s NRI S I _NR for i=1,..,R (Equation 24) NN B1 l BJ] b1 ] b]1 N B2 B b2 N = . ,B= . ,B= . ,b = ,

Ns] _BR] Bs] _bR bs where P, (D are the RS-dimensional vectors for code and phase measurements, N standards for the RS-dimensional integer vector for all LOSs at each frequency; b for the RS-dimensional code UHD vector and B for phase UHDs for all LOS at each frequency, from the column of the matrix Z in

(23). [0061] Defining the UD, SD and DD operator matrices s as follows

1 0 -- 00 0 - 0 00 -- 0] 0 0 --- 0 1 0 -- 0 00 -- 0 Do=. 000 0 10. . . (Equation 25)

0 0 --- 0 0 0 --- 0 --- 1 0 --- 0_]

Di =[S 0 0 --- 0] (Equation 26)

-S1 S 0 -- 03 -1 1 0 . 0]

-So0 S 0.10 .-. (Equation 27)

-Sl 0 0 --- Si) -- 1 0 0 --- 1)

where Do is a R-by-RS matrix; Di is a (S-1)-by-RS matrix and D 2 is a (R-1)(S-1)-by-RS matrix.

[0062] Applying Do and D 1operators to the UD vector N, the following results are obtained:

N1 N 12- N1 VN 1, D N DN.2AD N=1. = (Equation 28)

N } NL1 VN

[0063] Applying the D 2 operation to the UD vector N+B, the following results are obtained:

AVN] NN 2,1 2N21(I -(Ns N) 1)

D 2 (N+B)=D 2 N+D 2B= : = (Equation 29) AVNA2, N 2- NI-(N 12- N1)

AVNs _Ns - NR'-(Ns - N1)]

where DON gives the R-dimensional UD integer vector, DN is the (S-1)-dimensional SD integer

vector, D 2 (B+N) is the (R-1)(S-1)-dimensional DD integer vector. It must be noted D 2N may not be

the same as D 2(B+N) as D2 B may not be zeros in all elements.

[0064] Step 620 maps UD, SD, and DD integers to all LOS directions for LOS integer consistence. Inverting the matrix D and portioning the inverse matrix into three components according to the

column sizes of No, VN and AVN,

D-1=[U0 U1 U 2 ], (Equation 30)

the RS-dimensional LOS integers can be computed with the following equation

N = UODON+ UID 1 N+ U 2 D 2 (N+B). (Equation 31)

[0065] Step 630 adjusts the LOS UHDs with adjusted LOS integers from consistence. The adjusted B is

given by

B = B± N - N, (Equation 32)

where the over-line of N and B indicates the network adjusted LOS integers and phase delays.

[0066] FIG. 8 illustrates the computing system 700 that comprises steps to compute the satellite

and receiver specific hardware delays in combined or original signals.

[0067] Step 710 of FIG. 8 obtains the code UHD offsets in DD code signals. Substituting the LOS integers in the code dominated combinations, performing DD operations, the code offsets are obtained as follows:

AVb 13 =D 2 (P 13- ( 13) - 13 +N 13 )] (Equation 33)

AVb 32= D 2[(P32 -(D32)-)2 32+N 32 )] (Equation 34)

AVbF =D2 PIF IF )13 13 1 1 1 (Equation 35)

[0068] These code offsets remain constants or slowly varying, thus being fitted as function of time over a long period of observations or updated daily and represented by AV6 3 (t), AV632(t) and AVbIF Q)

[0069] Step 720 forms the observation vectors for six types of combined UHD state vectors by removing the code-offsets from the UHD equations,

Y 1(t (P13 - D 13 - N 13 )-U 2AV6 13(t) (Equation 36)

Y2 (t) (P32 -( D- 32 N232) - U 2AV6 32 (t) (Equation 37)

Y 3 W=I -( EIF- IF)- 4 N13 )]- U2AVbIF(t) 1 (Equation 38) fl-f3

Y4 (t)=I fi 21 F13 f 22 I 32 +(2NI- 2 A3 13 22 2 N 32 ) (Equation 39)

4A >1f- (Equation 40) Y 5 (t)= DIF-P-AWt+( N 13 1 1)

1f3 Y6 - =(I - (13 +GNl - 43 13+ -I(t) (Equation 41)

[0070] The observation equations and statistic models for satellite and receiver specific UHDs.

YJ (t)=HZUHDJ +E (Equation 42)

where

-1 0 0 ... 07

UHD1 0 --- 1 0 --- 0 101 0 UHDs

H ,Z UHDJ UHD 2 (Equation 43) R

-1 --- 0 0 --- 1 UHDR]J

0 --- -1 0 -.. 1]

[0071] The term {ZUHDJ} in Equation 51 are unrelated, but their error terms are generally cross

correlated. For Equation 51, it is assumed that error vectors are contemporaneously correlated, the covariance matrix between two USD observable vectors are

E(s,= ;C ss)= J =1, 2,..., 6 I s, JK,J=1,2,..., 6,K =1, 2,..., 6 (Equation 44) Cov(CJ,,)= aJ,KIRS

where I RS is the RS-dimensional unit matrix.

[0072] Step 730 performs least square estimation to obtain time series for each satellite and receiver specific UHD states.

[0073] Because the noise vectors in Equation 42 could be cross correlated, theoretically run the weighted least squares estimation should be applied to solve all the UHD vectors in Equation 42 together as above. However, due to the same matrix H for the USD parameters in all the models, the weighted least square solutions turn out to be theoretically identical to the least squares estimate without considering the correlation between models (Davidson, R., MacKinnon J. G., 1993). That is, we have

ZuDj =(H TH)1 HTY , (Equation 45)

and the residual vector for the Jth equation

r = r - H(HT -1HT]Yj. (Equation 46)

[0074] The residuals contain the multipath errors and computed geometric range errors. In an

alternative embodiment, the historical data of a moving window of hours may be used to train

the error patterns and generate corrections for these errors for these residual terms in the next

few hours. As a result, the systematic residual errors in Equation 51 can be removed before the least-square estimation with Equation 45.

[0075] There are six residual vectors for all UHD solutions. The estimators of the covariance components are given by

62- _ jj -J - i K (Equation 47) J (S -1)(R -1) 1 JK= (S -1)(R -1)

[0076] The covariance matrix for the estimateZUHDJ for each GF equation is

Cov(ZUHDJ) = (H H)~1 (Equation 48)

[0077] The network-adjustment solution (Equation 54) gives the time series for each element of

ZUHD(t). The covariance matrix Equation 48 shows the performance factor of the network

adjustment results from a single epoch. The Network Dilution of Precision (NDOP) defined as the

square root of the diagonal elements of the matrix (HTH)~l can describe this performance

factor. For a satellite-specific UHD parameter, NDOP depends on both the number of stations and the number of satellites commonly observed. This means that network-adjustment process can

also improve the accuracy or precision of single-epoch UHD samples at every epoch. It must be

noted that in the receiver specific UHDs from Y(t), [UHD 2 ,UHD 3,...UHDR 5 , are defined with

respect to their receiver clock offsets.

[078] Step 740 converts the UHD estimates in six combined signals into the UHD estimates in the

original code and phase signals

UHD,P TUHD,1

UHDP2 UHD,2 HDP _ -1 U(Equation 49) UHD , UHD,4

UHD,# 2 UHD,5

HD, 3 HD,6

where XUHD,J(t) is the (S+ R-1)-by-1 UHD vector, include S satellite-specific UHDs and (R-1)

receiver specific UHDs.

[0079] Step 750 collects the single-epoch satellite specific UHDs and fitted UHDs for each reference receiver to form single-epoch SSR corrections for 6S satellite specific UHDs out of

XUHDJ(t), for J=1 ,P2,3,01,02,03 and a RS-dimensional slant ionosphere-delay 1(t) and RS

dimensional troposphere vector dT(t).

[0080] Step 760 estimates of UHDs parameters in the time domain with the network-adjusted UHD samples at each epoch. Regardless of post-processing or real-time processing, the UHDs of

six original signals may be modelled over a continuous observation period to achieve filtered results. How to model the UHDs in time domain is a less studied issue. In this treatment, in order

to model the possible time variations of UHDs, each element of the state vector XUHDJ is more

generally augmented by a p, -degree polynomial function of time with (pJ+1) parameters,

J=PP2, P3,01,2,03. We introduce the following matrix and vector notations:

mHD,Jl

1 is, 2 ' ' I,md] , (Equation 50) L HD,J (tK)

X O' X 0 ,2 X Oma

X X --- X _mX1 Xi, 2 -'- X1,m , (Equation 51)

Lx,, XP1 XPJMd]

1 (t,-t") - (t -t 0 )P'1 G = (Equation 52) -1 (tK -t1) --- (t -tO)"3

[0081] Considering the independent samples between epochs and the cross-correlation of

between the error time series after the network adjustment for each USD component, we have the linear equations:

X' 1 =GJXq,+ej 1 i=1,...,md , (Equation 53)

and the covariance matrices

Cov(X,, )= o- I, i= k=1,...,m - ip J') 'i d(Equation 54) Cov(XiJ,k) ik'n ik,i,k =1,...,md

where md=R+S-1 is the dimension ofXUHDJ; n represents the samples if the whole data period ti

to tn.

[0082] The weighted least square solutions of the problem (Equation 40) are given by

, = (G G )GT iXJ, (Equation 55)

Cov(X F) (G Gj) 1, =J (Equation 56)

[0083] Defining the residual vector from the ordinary least square solution,

j,i= [In - G (GGj)-1 G ]X , (Equation 57)

the estimators of the covariance components over the whole data period are given by

2= 1 ,jir r lrJk (Equation 58) n- p,-1'Jk n- P-1

At all the time instants ti to to, the estimated

Xi = GJXi (Equation 59)

its covariance matrix can be expressed by

Cov[X,]=6JG,[GjGJ]'G, for i=1,2,..., md (Equation 60)

This processing is run for. J = PI, P2, P3,01,02,03-

[0084] Step 780 coverts the fitted XUHD,J to the fitted cUHD, can J be also obtained through the matrix

A. As a result, the total hardware delays in the combined GFIF measurements can be obtained as follows:

13 = HZUHD,1 2 13'

32 = HZUHD,2 + U 2AVb 32 , (Equation 61) YIF = HZUHD, 3 + U2AVbe

Y$ = HZUHD,4

[0085] Their predicted values can be used as the initial values in the determination of LOS and integers with Equations 10a, 10b and 10c.

[0086] FIG. 9 illustrates the user receiver embodiment 800 comprising the methods 810 and 820 for estimation of user states and integers with SSR correction messages 540, and with the user

data streams from a single, dual or triple frequency receiver 805. The subsystems 810 and 820 comprises steps of user receiver precise state estimation by applying the SSR corrections from a

single reference station to calibrate the same error terms in the GNSS code and/or phase

observables at a single, dual and triple frequency with a user receiver.

[0087] Method 810 is a single-reference RTK method with the single-epoch SSR corrections. Step 812 applies the single-epoch SSR corrections to calibrate the same terms in code and phase signals

of a user receiver, which forms DD equations for RTK positioning. Noting the code UHD has two terms and the phase UHD takes the satellite specific only as shown in Equation 7, the linear

observation equations are expressed as follows:

P, -P, -Vp,'-VdY7T''- f-(-V +V * H u)= dX+Vb'+e -e1 fL eC

for L=1,2,3; s=2,3,...,S (Equation 62) -Vp' 1 -V dT' += q, A L 1 H-sdX - ALVN + e, -e-,

for L=1,2,3, s=2,3,...,S (Equation 63 where the Vbi'', and Vbs:1 are the code-bias offsets that still exists in the code UHD and user code measurements. Substituting the terms Vdfsland vg;' as defined by Equation 17b and Equation

b into the above phase measurements, the following equations are obtained:

(#,-#) -(i, -#r,)-(Vp -Vp)= H"ldX -(VN VNLr,)+(e-s 1)-(s -s,") for L=1,2,3;s=2,3,...,S (Equation 64)

[0088] The equation 64 are the exact forms of DD phase linear equations. This shows how the

single-epoch SD SSR corrections can support RTK positioning that is based on DD code and phase measurements. But using the single-epoch SSR corrections instead of single-epoch measurements

from the reference stations can reduce the data size because the single-epoch SSR components vary slowly and the lower data rates can be adopted when being delivered to users.

[0089] Step 814 performs the determination and/or calibration processes for the code-bias offset

with respect to the reference station, that is, (Vb' -Vb). Determination is performed from the

difference between (equation 62) and (equation 63) where the integers are correctly resolved

with the historical data over longer period. Calibration is to apply a known code-bias to remove its effect when the same receiver and same reference stations are used in RTK services.

[0090] Step 816 performs the standard RTK processing with single-frequency, or widelane, dual frequency measurements.

[0091] Method 820 is a PPP-RTK processing approach with the fitted SSR corrections. Step 822

applies the fitted or predicted SSR to calibrate the same terms in code and phase signals of a user receiver, to enable user-end PPP-RTK processing.

[0092] When fitted or filtered UHD and SD ionosphere and troposphere-delays from a reference

receiver are applied to calibrate the same terms in code and phase signals of a user receiver, the

linear observation equations are expressed as follows:

P -P -Vp; l -VdYs P~f -P IMVt~ -, L,, -(-VkL+VL, )=HgdX +Vbi +e -e

for L=1,2,3; s=2,3,..., S (Equation 65)

$ -L, -Vp" -VdY7' + -ALVb" = Hs dX, - ALVN +± fL L L, ,

for L=1,2,3; s=2,3,..., S (Equation 66)

[0093] The equations 65 and 66 are the forms of SD code and phase equations for PPP-AR or PPP RTK processing. The precision of DD and SD phase measurements Equations (64) and (66) are

almost the same after SSR corrections being applied, but the code measurements Equation 65 will have higher precision than that of Equation 62, due to use of fitted (-Vb'±+Vbj). The fitted or filtered SSR corrections are slowly varying. Predicted values can be used in real time PPP-AR processing, thus including the effects of data latency in the transmission from a reference station and a rover.

[0094] Step 824 performs the determination and/or calibration processes for the code-bias offset

with respect to the reference station, (Vb1 -Vb ). Determination simply takes average of the

difference between Equations 66 and 65 with the historical data over which the integers are

correctly fixed as well. Calibration is to apply a known code-bias to remove its effect.

[0095] Step 826 performs the standard PPP-AR and PPP-RTK processing with single-frequency

code and phase measurements, dual-frequency or triple-frequency code and phase measurements. For long user-reference distance where the effects of the distance-dependent errors on the rover

receiver grow beyond several centimetres, there are dual-frequency method and triple-frequency methods for achieving reliable ambiguity fixing and decimetre positioning results with widelane

phase combinations. In the step 826 the dual-frequency method makes use of more precise combined code measurements P 13 and widelane phase measurements. The observation equations

for each set of SD code and phase signals are as follows:

'I3,1-W1,- Vp"- Vds'- r__ (Vs+Vbsl)= HsdX,,+AVbiS3 + s, -1 forr 1323.r S e3 67)

for s=2,3,..., S (Equation 67)

0"-0 ," -Vpl;)- V d Vq- 13 VhB`= Hs~'dX,-2" 3 VN 3 ±+,s,, - 4,

for s=2,3,..., S (Equation 68)

[0096] The widelane integer ambiguity term VN 3 should be reliably fixed almost instantly due to

the long wavelength (0.862 m with GPS signals) and the low noise of code measurement P 13 . The

user position estimation with Equation 67 can achieve to the accuracy of ten to a few tens of

centimetres immediately, for normal the dilution of precision (DOP) and growth of distance dependent errors within 10 centimetres.

[0097] In another embodiment, the triple frequency methods can make use of widelane GBIF 26 in

Table 2 for user AR and position estimation, the observation equations are expressed as follows:

#L,, - ,u, Vp,''' S[for Vdi'''r1 -[g&VEI +p8A2VhR2]'''"+p81-Z1VN s=2,3,..., S =H"dX,- 11A 13 VN 3 +WsL,u -WIL,u

(Equation 69)

[0098] As the precision of DD phase measurements WLare about 10 cm and with the wavelength

of 3.9 m, the widelane ambiguity terms in (Equation 76) can be resolved instantly. After the

ambiguity is resolved, the position states should be estimated to the accuracy of a few tens of centimetres immediately, depending on the dilution of precision (DOP) and growth of distance

dependent errors.

[0099] FIG. 10 illustrates the user receiver embodiment 900 comprising the methods 920 and 930 for estimation of user states and integers by with SSR correction messages, and with the user data

streams from single, dual, and triple frequency receivers. The system 900 comprises subsystems and steps of user receiver precise state estimation by applying the SSR corrections from a network

of reference stations to GNSS code and/or phase observables at a single, dual, and triple frequency with a user receiver.

[0100] Step 910 builds the ionosphere interpolation models with single-epoch slant ionosphere delays and the troposphere interpolation models with single-epoch slant troposphere-delays from

the surrounding reference stations. The interpolated SSR correction terms are expressed as

functions of time and location, that is, dTi(tx), (, f2or f2 df"(t, x), q(t,x)where "x"indicates

the horizontal locations of a receiver with respect to the reference station locations.

[0101] Method 920 is a network-based user RTK procedure. Step 922 applies the adjusted UHD

and interpolated ionosphere and troposphere SSR to calibrate the same terms in code and phase signals of a user receiver, to support network-based user RTK positioning processing. The SD linear

observation equations are formed

P -P-;'-Vd''(t,x)- vq (t, x) - '+V )=H"dX, ±+Vb"+e, - eL,

for L=1,2,3;s=1,2,...,S (Equation 70) ##-Vp,''' -V di "(t'x)+ - AVE"'))= H" dX. - AVN , +sq", ELu frL123 s ,

for L=1,2,3;s=2,3,...,S (Equation 71)

[0102] The equations (70) and (71) are consistent to the network-based RTK positioning models.

But using the interpolated ionosphere and troposphere corrections from multiple stations instead of single-epoch measurements from a single receiver can improve the precision of the ionosphere

and troposphere-delays or increase the inter-receiver distances

[0103] Step 924 performs the determination and/or calibration processes for the code-bias offset

with respect to the reference station. Determination simply takes average of the difference between Equations 71 and 70 over the historical data and updated regularly. Calibration is to

apply the determined or updated code-bias to remove its effect.

[0104] Step 926 perform the standard RTK processing with single-frequency, or widelane, dual frequency measurements. Users obtain the network-based RTK positioning solutions.

[0105] Method 930 is a network-based user PPP-AR or PPP-RTK procedure. Step 932 applies the

r 'f) 22 to calibrate the same terms in code and phase fitted and adjusted UHDs and di<(t,x),

signals of a user receiver to support network-based PPP-AR and PPP-RTK processing.

[0106] The linear observation equations are expressed as follows:

1 - - P - Vp'-Vdi' - (t,x)-v='(t,x) -Vdi"Lt, X)- C-x, - - '(t)+V (t)) = H;'dX,,+AVb_,,+ek,,-e,

for L=1,2,3;s=2,3,...,S (Equation 72) $ -$ -V -df"' 1 (t,x)+V4;"1 (t,x) LVbt)= H'dX, - ALVN , +

for L=1,2,3, s=2,3,..., S (Equation 73)

[0107] Equations 72 and 73 are the forms of SD code and phase equations for network-based PPP

AR and PPP-RTK positioning. The network adjusted and fitted/filtered UHD corrections are slowly varying. Predicted values can be used in real time PPP-AR or PPP-RTK processing, thus including

the effects of data latency in the transmission from a reference station to a user receiver. Both

di"(t, x), '2 are time and location dependence, but ' is faster time varying and f s location sensitive.

[0108] Step 934 performs the determination and/or calibration processes for the code-bias offset with respect to a reference station. Determination simply takes average of the difference between

Equation 72 and 73 over with the historical data. Calibration is to apply a known code-bias to

remove its effect.

[0109] Step 936 performs the PPP-AR processing with single-frequency, or widelane, dual frequency measurements at a user-end. This is to obtain the PPP-AR and PPP-RTK solutions.

Claims (5)

1. A method of processing data files and streams from a single GNSS receiver for generation of

space-state representation correction components for the receiver tracking a plurality of satellites transmitting code and phase signals at two or more frequencies, the method

including the steps of: a. reformatting the code and phase signals into three groups of linear combinations, including

geometry-free and ionosphere-free (GFIF) combinations, geometry-free and ionosphere

present (GFIP) combinations, and geometry-based and ionosphere-free (GBIF) combinations;

b. using the known precise GNSS orbits, clock and station coordinates, initial hardware delay products in sequentially reprocessing the data streams from the receiver, determining

undifferenced (UD) and single-differenced (SD) integer ambiguities, and slant and SD ionosphere and troposphere delays with phase-only GFIP and GBIF measurements,

respectively and recomputing uncalibrated hardware delays (UHD) in code and phase signals epoch by epoch, which are raw or single-epoch state-space representation (SSR)

correction components; c. applying the raw and fitted SSR corrections for code and phase UHDs, slant or SD

ionosphere and troposphere delays to calibrate the same delay terms in the user code and

phase measurements to enable Real Time Kinematic (RTK), Precise Point Positioning (PPP) with ambiguity resolutions (PPP-AR) or PPP-RTK positioning services, based on user's

single-frequency, dual-frequency or triple-frequency data; d. substituting the fitted SSR corrections for hardware and atmosphere- delays adjusted from

a network of reference receivers into the original code and phase observation equations of each reference receiver, to correct the same delay terms, to enable recovering the integer

fixed original phase measurements and the code-offset corrected original code measurements, and obtain the residuals of original code and phase measurements which

include mainly the effects of multipath errors and random observation noises.

2. The method of Claim 1 wherein the single-receiver based computing makes use of the GFIF and GBIF observables to resolve and verify UD and SD integers further comprising: a. using three GFIF models to determine the UD float and integer ambiguities for each receiver-satellite direction, but the models may be corrected with initial UD code and phase biases to improve the accuracy of the floating ambiguities; b. performing a least-squares adjustment with known SD integers as constraints to obtain the consistent Line-of-Sight (LOS) integers for each reference receiver; c. using the integer-fixed GFIP phase measurements and initial UHD products to compute SD ionosphere-delay for a reference receiver; d. computing a slant ionospheric delay of the receiver to the reference satellite path with respect to zero UD UHD value or an external ionospheric model value as the datum of the slant ionosphere-delay for each receiver; e. mapping the UD and SD ionosphere delays into all LOS directions; f. using the integer-fixed GBIF phase measurements to compute SD troposphere-delays for a reference receiver; g. setting the slant residual troposphere-delay in the reference satellite direction with respect to zero UD HUD value or a different troposphere-model as the datum of slant troposphere delay for each receiver and mapping the UD and SD troposphere delays into LOS directions; h. removing the determined SD integers from all selected SD GFIF, GFIP and GBIF observables, removing the SD ionosphere-delays from the GFIP observables, and removing the SD troposphere delays from the GBIF measurements, then forming linearly independent equations for combined UHD equations for all SD observables ; i. converting the combined UHD observables to UHDs in the original code and phase signals; j. performing a fitting or filtering process for each single epoch UHD time series for forming function of time for prediction and interpolation; k. formatting the SD delays and SD UHD function of time into SSR correction formats; and I. mapping all SD integers and delays to LOS directions by appending the UD quantities, then obtaining slant delays and LOS UHDs for SSR corrections.

3. The method of Claim 1 wherein the network-based computing makes use of single-epoch combined UHD samples from each single receiver to derive adjusted double-differenced (DD) integers for the necessary baselines, which are then mapped to all LOS directions, by appending UD and SD integers from the single-receiver based processing the method comprising the steps of: a. performing DD operations over the obtained LOS UHD and LOS integers to derive DD integers for all necessary baselines to enable consistence between all LOS integers; b. mapping the DD integers and the first station SD integers to adjusted LOS integers by appending the LOS integers for all receiver-reference satellite directions; c. applying the adjusted LOS integers to lumped sums of LOS integers and UHDs to adjust the

LOS UHDs; d. performing DD operation onto code LOS UHD time series and obtaining DD code UHD

offsets for the network with long-term or historical data; e. removing LOS integers, slant ionosphere delay and slant troposphere delays from

combined measurements of all LOSs, forming linear equations for satellite-specific and receiver-specific UHDs and performing least-square estimation to determine satellite- and

receiver specific UHDs;

f. computing the residuals of the least square solutions of the above linear equations and obtaining the covariance matrices of satellite-specific and receiver specific UHDs; and

g. converting the obtained combined UHDs to original code and phase UHDs and performing fitting process to each original and combined UHD time series for a function of time.

4. A method of applying SSR corrections for hardware and atmosphere delays in a single

reference receiver or a network of multiple receivers to correct the same delays terms in a

user-receiver and performing precise positioning, the method comprising the steps of: a. interpolating the ionosphere and troposphere delays from multiple nearby-reference

receivers at the user receiver location epoch by epoch; b. directly applying the single-epoch and interpolated SSR corrections from the reference

receiver, the user receiver forms DD linear equations for original or combined code and phase measurements, so the users can perform RTK positioning with single-frequency,

dual-frequency or triple frequency measurements epoch by epoch; c. applying the fitted and interpolated SSR corrections from the reference receiver instead,

the user receiver forms SD linear equations for original code and phase measurements,

then perform PPP-AR or PPP-RTK positioning with single-frequency, dual-frequency or triple frequency measurements;

d. calibrating user-receiver code offsets with respect to a reference receiver from the difference of code and phase signals at single, dual, and triple frequencies. This treatment is done with historical data after the SD integers are fixed and removed. The code offsets are updated regularly for RTK, PPP-AR and PPP-RTK processing.

5. The method of Claim 4 comprising the further benefits of: a. using the determined single-epoch or fitted SSR corrections instead of raw code and phase

measurements from the reference stations can significantly reduce the overall data size of messages over a period of operations, thus reducing the required bandwidth for data

transmission to users;

b. using the determined single-epoch or fitted SSR corrections from a reference receiver or network in the user SD original phase signals, theoretically enabling the instant or fast

ambiguity resolution at the user terminals; and c. using the determined single-epoch or fitted SSR corrections in the reference receivers

allows the multipath errors in the residuals in each code and phase measurements to be modelled.

d. Using the single-epoch and fitted SSR corrections from a reference receiver or network, user's RTK, PPP-AR and PPP-RTK positioning services can all be supported