What you’ll build in this tutorial

A complete numerical baseline solution from one real Qelaro execution: every observation, matrix, and diagnostic remains inspectable.

One epoch, fully traced

Start with two receivers observing the same GPS satellites. Subtract within each receiver, then between receivers, to isolate the relative geometry. Build the covariance alongside the observations, solve float unknowns, search integer ambiguities, and inspect the constrained result.

Numbers before abstraction: first G04 relative to G09, then the full system.

Common satellites5

ReferenceG09

DD pairs4

DD observations8

Float unknowns7

Float DOF1

Lesson convention

For this fully traced one-epoch example, the rover RINEX approximate ECEF position is held fixed and the rover-to-base baseline is estimated. This follows the traced Qelaro implementation. Operational PPK often uses a known base and estimates the rover instead.

\mathbf b=\mathbf r_{base}-\mathbf r_{rover},\qquad \hat{\mathbf r}_{base}=\mathbf r_{rover,fixed}+\hat{\mathbf b}

Double-difference processing pipeline

Inputs and sign convention

2024-11-14 09:00:00 GPS. GPS L1, C1C metres, L1C cycles. λ = 0.19029367279836487 m/cycle.

DD=(B_j-B_m)-(R_j-R_m)

Where this lesson ends

Candidate baseline length: 810.359239 m. Integer candidate: [13, 7, 27, 1]. Ratio diagnostic: 1.445276.

Warning · Candidate-fixed; not validated

The candidate and centimetre-level formal sigmas do not establish external accuracy or operational fixing reliability.

1. Problem and model

Code measures a travel-time-based distance. Carrier phase measures a fractional wave plus an unknown whole-cycle count; multiplying by wavelength expresses it in metres.

Code pseudorange

P_r^s=\rho_r^s+c(\delta t_r-\delta t^s)+T_r^s+I_r^s+\mathrm{other\ terms}+\epsilon_P

Code gives the distance scale without a whole-cycle ambiguity, but is noisier than phase.

Carrier phase in metres

\lambda\Phi_r^s=\rho_r^s+c(\delta t_r-\delta t^s)+T_r^s-I_r^s+\lambda N_r^s+\mathrm{other\ terms}+\epsilon_\Phi

Phase is precise, but the integer N must be estimated. The ionosphere enters code and phase with opposite signs.

Read the observation equation

Symbols and physical units
Symbol	Meaning	Unit
r / s	receiver / satellite	—
ρ	exact geometric range ‖satellite − receiver‖	m
c	speed of light, 299792458	m/s
δtᵣ / δtˢ	receiver / satellite clock offset	s
T / I	tropospheric / ionospheric delay	m
λ / Φ	wavelength / measured phase	m/cycle / cycle
N	whole-cycle ambiguity	cycle
ε	measurement noise and remaining error	m

First subtraction: same receiver

SDᵣ(j,m) = yᵣ,j − yᵣ,m. The same receiver-clock term cancels. The satellite-clock difference between j and m remains.

Second subtraction: two receivers

DD(j,m) = SDbase − SDrover. Simultaneous observations remove the common satellite-clock terms. Spatially correlated atmosphere and orbit errors are reduced, not guaranteed to vanish.

Warning · Structural equations, limited correction scope

These equations explain cancellation. This traced DD_linear lesson uses broadcast NAV states and the core DD numerical path; it does not apply a complete operational zero-difference correction ledger. Full PPK may explicitly model Sagnac, clocks, antenna effects, troposphere and other terms.

2. Input data (RINEX)

Read the real observation files, extract the shared epoch, and propagate broadcast navigation parameters to satellite states.

Base observations · 36.24O

RINEX 3.02. The first GPS observation fields are C1C (m), then L1C (cycles); trailing digits encode signal quality.

     3.02           OBSERVATION DATA    M (MIXED)           RINEX VERSION / TYPE
  3246705.6964  4052996.6836  3692516.1666                  APPROX POSITION XYZ 
G   16 C1C L1C D1C S1C C2W L2W D2W S2W C2S L2S D2S S2S C5Q  SYS / # / OBS TYPES 
> 2024 11 14 09 00  0.0000000  0 39       0.000000000000
G04  21325621.891 8 112066897.717 8     -1691.258 8        48.383 8  21325630.250 9  87324885.628 9     -1317.838 9        56.226 9  21325630.336 8  87324885.635 8     -1317.857 8        53.559 8  21325634.332 8  83686366.378 8     -1262.959 8        51.124 8
G09  20149741.844 8 105887607.793 8       -95.203 8        52.968 8  20149748.410 9  82509850.730 9       -74.205 9        56.226 9  20149748.785 9  82509851.734 9       -74.146 9        54.030 9  20149749.148 9  79071943.958 9       -70.984 9        54.060 9

Rover observations · 15.24O

RINEX 3.02. The first GPS observation fields are C1C (m), then L1C (cycles); trailing digits encode signal quality.

     3.02           OBSERVATION DATA    M (MIXED)           RINEX VERSION / TYPE
  3246374.5864  4052665.3939  3693176.0609                  APPROX POSITION XYZ 
G   16 C1C L1C D1C S1C C2W L2W D2W S2W C2S L2S D2S S2S C5Q  SYS / # / OBS TYPES 
> 2024 11 14 09 00  0.0000000  0 33       0.000000000000
G04  21325674.301 8 112067164.390 8     -1691.398 8        51.422 8  21325680.582 9  87325076.873 9     -1318.094 9        56.226 9  21325680.453 9  87325081.886 9     -1318.113 9        54.067 9  21325684.051 8  83686546.923 8     -1262.975 8        52.228 8
G09  20149603.512 8 105886884.690 8       -95.852 8        52.498 8  20149608.793 9  82509277.969 9       -74.344 9        56.226 9  20149609.379 8  82509273.972 8       -74.344 8        53.030 8  20149608.891 8  79071387.938 8       -71.197 8        53.381 8

Convert cycles into metres

\Phi_{m}=\Phi_{cycle}\,\lambda_{L1}

Base G04 phase

112066897.717 cycles × 0.19029367279836487 m/cycle

= 21325621.56568662 m

The phase is now in the same length unit as code. This does not remove its integer ambiguity.

Base navigation · 36.24N

Actual G04 broadcast ephemeris excerpt. D denotes a decimal exponent. These parameters feed the Kepler propagation, not the final DD solver directly.

     3.02           N: GNSS NAV DATA    G: GPS              RINEX VERSION / TYPE
G04 2024 11 14 10 00 00 0.479244161397D-03 0.648014975013D-11 0.000000000000D+00
     0.115000000000D+03 0.178125000000D+02 0.498413618071D-08-0.136268651039D+00
     0.728294253349D-06 0.314607680775D-02 0.479444861412D-05 0.515367835236D+04
     0.381600000000D+06 0.745058059692D-08 0.114317989652D+01 0.167638063431D-07
     0.965789022041D+00 0.288125000000D+03-0.300042464309D+01-0.830748889740D-08
     0.106790162525D-09 0.100000000000D+01 0.234000000000D+04 0.000000000000D+00
     0.200000000000D+01 0.000000000000D+00-0.465661287308D-08 0.115000000000D+03
     0.374406000000D+06 0.400000000000D+01 0.000000000000D+00 0.000000000000D+00

Rover navigation · 15.24N

Actual G04 broadcast ephemeris excerpt. D denotes a decimal exponent. These parameters feed the Kepler propagation, not the final DD solver directly.

     3.02           N: GNSS NAV DATA    G: GPS              RINEX VERSION / TYPE
G04 2024 11 14 10 00 00 0.479244161397D-03 0.648014975013D-11 0.000000000000D+00
     0.115000000000D+03 0.178125000000D+02 0.498413618071D-08-0.136268651039D+00
     0.728294253349D-06 0.314607680775D-02 0.479444861412D-05 0.515367835236D+04
     0.381600000000D+06 0.745058059692D-08 0.114317989652D+01 0.167638063431D-07
     0.965789022041D+00 0.288125000000D+03-0.300042464309D+01-0.830748889740D-08
     0.106790162525D-09 0.100000000000D+01 0.234000000000D+04 0.000000000000D+00
     0.200000000000D+01 0.000000000000D+00-0.465661287308D-08 0.115000000000D+03
     0.376896000000D+06 0.400000000000D+01 0.000000000000D+00 0.000000000000D+00

NAV → transmit time → ECEF → LOS

The Qelaro wrench path estimates signal transmit time from code and broadcast clock information, selects ephemeris, and propagates satellite ECEF coordinates. For base G04:

Epoch satellite state

Epoch receive time = 378000 s GPS week

Transmit time = 377999.92838615703 s

Signal time difference = 0.07161384297069162 s

Broadcast satellite clock = 0.00047922526818618074 s

ECEF = [846359.717513945, 24131563.233674087, 10905032.09222802] m

The signal time difference above is from the traced path; it is not simply P/c. Its clock handling belongs to satellite-state propagation.

NAV orbit tracks

Oblique ECEF projection, not a skyplot. NAV propagated every 120 s over ±30 min; large dots show epoch GPS time. DD states use signal transmit time instead.

A numeric line of sight

The design model uses base-frame satellite states and the fixed rover coordinates, exactly as geometry_row does. Normalize the ECEF difference to obtain u.

\mathbf u_j={\mathbf s_j-\mathbf r_{rover}\over\|\mathbf s_j-\mathbf r_{rover}\|}

G04: three subtractions, a norm, three divisions

X: 846359.717513945 − 3246374.5864 = -2400014.868886055 m

Y: 24131563.233674087 − 4052665.3939 = 20078897.839774087 m

Z: 10905032.09222802 − 3693176.0609 = 7211856.03132802 m

ρ = √((-2400014.868886055)² + (20078897.839774087)² + (7211856.03132802)²) = 21469352.045359004 m

uX = -2400014.868886055 / 21469352.045359004 = -0.11178795074091966

uY = 20078897.839774087 / 21469352.045359004 = 0.9352353903067376

uZ = 7211856.03132802 / 21469352.045359004 = 0.3359140050473482

A unit direction is dimensionless. Rotate a receiver-specific LOS to ENU; azimuth = atan2(E,N), elevation = asin(U). The receiver-frame Az/El determines weights and reference choice.

All satellite states at transmit time

Receiver satellite states; ECEF metres, GPS week seconds, degrees
Receiver	PRN	t transmit	X	Y	Z	Az	El
Base	G04	377999.928386	846359.717514	24131563.233674	10905032.092228	99.187117	47.102246
Base	G07	377999.929523	16535451.776103	9559024.856856	19168552.314801	307.453856	65.599553
Base	G08	377999.921132	-5883706.490332	25640105.566921	582780.305009	113.553399	17.923853
Base	G09	377999.932345	8585785.851513	17860591.721371	17557100.591124	55.056628	74.419671
Base	G16	377999.921858	-10291804.145127	11037476.839404	21489206.413478	42.581565	19.597189
Base	G20	377999.920802	14572581.830868	-7889378.980343	20804966.630072	313.072698	19.924584
Base	G30	377999.925168	24270697.178789	1612147.408554	11130157.986078	268.557928	36.899820
Rover	G04	377999.928386	846359.717622	24131563.233453	10905032.092712	99.198397	47.101383
Rover	G07	377999.929523	16535451.774717	9559024.857325	19168552.315762	307.436441	65.604710
Rover	G09	377999.932345	8585785.851206	17860591.722367	17557100.590263	55.084021	74.425578
Rover	G16	377999.921860	-10291804.148830	11037476.835244	21489206.413922	42.584842	19.603357
Rover	G30	377999.925168	24270697.179019	1612147.408434	11130157.985592	268.550240	36.899059

3. Satellite selection

A double difference needs both receivers to observe the target and the reference at the same epoch.

Common-satellite filtering

Base · 7 satellites

G04G07G08G09G16G20G30

∩ Rover · 5 satellites

G04G07G09G16G30

= Common · 5 satellites

G04G07G09G16G30

G08 and G20 disappear because no usable matching rover observation exists in this epoch. They cannot form DD pairs.

Choose the highest common base elevation

Real elevations at both receivers
Satellite	Base elevation (°)	Rover elevation (°)	Selection
G04	47.102246	47.101383	Selected target
G07	65.599553	65.604710	Selected target
G09	74.419671	74.425578	★ Reference
G16	19.597189	19.603357	Selected target
G30	36.899820	36.899059	Selected target

Maximum over the intersection

max(base elevations) = 74.4196711687952°

= G09 reference

A high-elevation reference has smaller assumed measurement noise. G09 is reused in every DD, so its noise will also create correlations.

Skyplot / reference selection

Receiver

Zenith is at the center. Radius = 90° − elevation. Dashed ring: 15° mask. G09 ★ is the selected reference.

Data summary

5 common satellites − 1 reference = 4 DD pairs. Each pair supplies code and phase: 4 × 2 = 8 observations. Next, subtract reference G09 from each target within each receiver.

4. Single differences from scratch

At one receiver, subtract reference G09 from target G04. This is satellite-to-satellite differencing, not a between-receiver single difference.

Same receiver, two satellites

SD_r(j,m)=y_{r,j}-y_{r,m}

Both measurements contain the same receiver-clock term cδtᵣ, so it cancels. The satellite-clock difference remains until the second subtraction.

Applied to code and carrier phase

Code is already in metres. Convert phase cycles to metres using λ before combining the observables with the geometric model.

Double-difference geometryG04 − G09

Solid: BaseDashed: RoverG09 referenceG04 target

Target

Rover → Base baseline. Five common satellites; selected target and reference highlighted. Conceptual, not to scale.

Code worked example: G04 − G09

Base code SD

21325621.891 − 20149741.844

= 1175880.0469999984 m

Rover code SD

21325674.301 − 20149603.512

= 1176070.7890000008 m

The similar million-metre SDs mostly describe the satellite geometry. Their difference will be only hundreds of metres.

Phase worked example: cycles first

Base phase SD

(112066897.717 − 105887607.793) cycles

= 6179289.923999995 cycles

× 0.19029367279836487 m/cycle

= 1175879.774923887 m

Rover phase SD

(112067164.39 − 105886884.69) cycles

= 6180279.700000003 cycles

× 0.19029367279836487 m/cycle

= 1176068.1230341755 m

The trace converts each phase observation before differencing. Subtracting cycles first is algebraically equivalent but can differ by a few nanometres in floating-point arithmetic.

All single differences in the trace

All target/reference SDs; G09 is the reference, values in metres
Target	Base code SD	Rover code SD	Base phase SD	Rover phase SD
G04	1175880.047000	1176070.789000	1175879.774924	1176068.123034
G07	984870.074000	984818.992000	984870.257671	984817.818253
G08	3370389.133000	No matching rover	3370393.154134	No matching rover
G16	3315272.558000	3314817.668000	3315272.828546	3314813.173228
G20	3482779.679000	No matching rover	3482778.495685	No matching rover
G30	2368905.371000	2369098.687000	2368905.139357	2369098.361649

Single-difference observation view

The numerical values are available in the accompanying tables. The interactive chart loads as you approach it.

What one subtraction achieves

Receiver clock cancels

The shared receiver clock is removed from both code and phase SDs.

Warning · What remains

Geometric differences, the satellite-clock difference, atmospheric differences, ambiguities and measurement noise still remain. The next step combines matching SDs from both receivers.

5. Double differences and differencing matrix

Subtract rover SD from base SD. The same reference G09 appears in each pair; retain the exact sign convention throughout the observations and covariance.

Double-difference definition

DD(j,m)=SD_B(j,m)-SD_R(j,m)

DD=+y_{B,j}-y_{R,j}-y_{B,m}+y_{R,m}

Under simultaneous observations, the satellite-clock terms remaining in the SD cancel between receivers. Atmospheric and orbit effects are only reduced according to how correlated they are.

Worked G04 double differences

Code DD

1175880.0469999984 − 1176070.7890000008

= -190.742000002414 m

Phase DD

1175879.774923887 − 1176068.1230341755

= -188.3481102883816 m

Phase − code = 2.393889714 m. This contains ambiguity plus code/phase noise and remaining model differences; it is not ambiguity alone.

All real DD observations

Four target/reference pairs, G09 reference; metres
Target	Code DD	Phase DD	Phase − code
G04	-190.742000	-188.348110	2.393890
G07	51.082000	52.439418	1.357418
G16	454.890000	459.655318	4.765318
G30	-193.316000	-193.222292	0.093708

DD observation plots

The numerical values are available in the accompanying tables. The interactive chart loads as you approach it.

Differencing matrix D · exact 8 × 20 structure

Stack the 20 independent undifferenced observations in the exported order. A row of D selects four measurements; every other coefficient is zero. Phase and code rows alternate for G04, G07, G16, G30.

\underbrace{\mathbf y_{DD}}_{8\times1}=\underbrace{\mathbf D}_{8\times20}\underbrace{\mathbf y_{UD}}_{20\times1}

D: blue −1, teal +1; negative is a valid coefficient8 × 20

Row / column	B G04 Φ	R G04 Φ	B G09 Φ	R G09 Φ	B G04 P	R G04 P	B G09 P	R G09 P	B G07 Φ	R G07 Φ	B G07 P	R G07 P	B G16 Φ	R G16 Φ	B G16 P	R G16 P	B G30 Φ	R G30 Φ	B G30 P	R G30 P
G04 phase
G04 code
G07 phase
G07 code
G16 phase
G16 code
G30 phase
G30 code

G04 phase × B G04 Φ = 1

Show full mathematical matrix

Explain a row

G04 phase: row-by-column product

+1 × B G04 Φ (21325621.56568662 m)

−1 × R G04 Φ (21325672.311871227 m)

−1 × B G09 Φ (20149741.790762734 m)

+1 × R G09 Φ (20149604.18883705 m)

= (1) × (21325621.56568662) + (-1) × (21325672.311871227) + (-1) × (20149741.790762734) + (1) × (20149604.18883705)

Sum of signed contributions = 21325621.56568662 + -21325672.311871227 + -20149741.790762734 + 20149604.18883705

= -188.3481102883816 m

This is the same DD as the SD subtraction. Round-off may affect the last few nanometres if four large phase values are combined in a different arithmetic order.

Next: one geometric row per observation

D constructs observed differences. H constructs the baseline and ambiguity model for those differences. D also propagates undifferenced variances into the correlated DD covariance.

6. Linear design matrix

Linear baseline-domain DD model with LOS vectors fixed at the lesson epoch. The exact range norm is nonlinear in receiver coordinates; double differencing does not make it exactly linear.

Freeze the LOS, then construct H

Evaluate satellite states and LOS at the approximate/fixed coordinates. Under this frozen, first-order baseline representation, y = Hx + ε is linear in the baseline components and ambiguity states.

\mathbf x=[dX,dY,dZ,N_{04},N_{07},N_{16},N_{30}]^T,\qquad\mathbf h_{xyz}=-\mathbf u_j+\mathbf u_m

Observations8

Unknowns7

Float DOF = 8 − 71

Row-by-row geometry walkthrough

Target satellite

Target unit LOS = [-0.111787951, 0.935235390, 0.335914005]. Reference G09 = [0.263253170, 0.680782993, 0.683543916]. Both use the same fixed rover origin.

−uG04 + uG09

X: −(-0.11178795074091966) + (0.26325317028501133) = 0.375041121025931

Y: −(0.9352353903067376) + (0.6807829927840786) = -0.25445239752265914

Z: −(0.3359140050473482) + (0.6835439159782231) = 0.34762991093087486

Each coefficient tells how a metre of baseline along that ECEF axis changes the modelled DD.

Phase row versus code row

H_{\Phi,j}=[h_X,h_Y,h_Z,\ldots,\lambda_j,\ldots]

\begin{bmatrix}0.3750411 & -0.2544524 & 0.3476299 & 0.1902937 & 0.0000000 & 0.0000000 & 0.0000000\\0.3750411 & -0.2544524 & 0.3476299 & 0.0000000 & 0.0000000 & 0.0000000 & 0.0000000\end{bmatrix}

Rows above: G04 phase, then code. Place λ = 0.19029367279836487 m/cycle only in this target’s ambiguity column. Code has no carrier integer ambiguity, so all four ambiguity coefficients are zero.

Full design matrix and observation vector

H: geometry columns followed by four cycle ambiguities8 × 7

Row / column	dX	dY	dZ	N04	N07	N16	N30
G04 phase
G04 code
G07 phase
G07 code
G16 phase
G16 code
G30 phase
G30 code

G04 phase × dX = 0.375041121025931

Show full mathematical matrix

Show y and explicit ordering

Observation vector y; same row order as H, in metres
Row	Value
G04 phase	-188.3481102884
G04 code	-190.7420000024
G07 phase	52.4394177087
G07 code	51.0819999985
G16 phase	459.6553181112
G16 code	454.8899999969
G30 phase	-193.2222924232
G30 code	-193.3160000034

\begin{bmatrix}-188.348110288\\-190.742000002\\52.439417709\\51.081999999\\459.655318111\\454.889999997\\-193.222292423\\-193.316000003\end{bmatrix}

One complete H row × x example

G04 phase prediction using the float state

dX: (0.375041121025931) × (329.96449120681626) = 123.75025268095533

dY: (-0.25445239752265914) × (331.82978936717865) = -84.43488547391759

dZ: (0.34762991093087486) × (-661.7168607071129) = -230.03257334907178

N04: (0.19029367279836487) × (12.449682739387208) = 2.3690958536524005

N07: (0) × (6.983324467405504) = 0

N16: (0) × (25.00173284312381) = 0

N30: (0) × (0.48040378538633366) = 0

= -188.34811028838166 m

Sum all seven products. Predicted DD plus the post-fit residual equals the observed −188.34811028838158 m. H is now ready, but observations are correlated: solve using R_DD, not equal independent weights.

7. Stochastic model and DD covariance

Measurements have different precision. Low elevation increases the assumed noise, and a shared reference makes double differences correlated.

Elevation-dependent observation variance

Nominal phase σ₀ = 0.006 m; code σ₀ = 0.7 m. Use the receiver-specific satellite elevation E in degrees.

\sigma_i={\sigma_0\over\max(\sin E_i,0.25)},\qquad \operatorname{Var}(y_i)=\sigma_i^2

Base G04 phase weighting

E = 47.10224609536898°

sin(E) = 0.7325695836802872

max(0.7325695836802872, 0.25) = 0.7325695836802872

σ = 0.006 / 0.7325695836802872 = 0.008190348239490323 m

Variance = (0.008190348239490323)² = 0.00006708180428412224 m²

Independent UD weight = 1 / 0.00006708180428412224 = 14907.172081486371 m⁻²

Lower elevation increases σ and decreases weight. The floor 0.25 avoids unbounded noise near the horizon.

From independent UD to correlated DD

Assume the 20 undifferenced code/phase measurements are independent in this lesson. R_UD is diagonal in m². Differencing propagates their noise using the same D as the observations.

\underbrace{R_{DD}}_{8\times8}=\underbrace D_{8\times20}\underbrace{R_{UD}}_{20\times20}\underbrace{D^T}_{20\times8}

For one entry, take row i of D R_UD and column j of Dᵀ: sum Dᵢₖ · σₖ² · Dⱼₖ across the 20 columns. Most products are zero.

Keep the off-diagonal entries

Using a diagonal approximation to R_DD would discard the shared-reference correlation and change the weighted problem.

Worked diagonal: G04 phase variance

R_{00}=\sum_{k=0}^{19}D_{0k}\sigma_k^2D_{0k}

Four nonzero products

(1) × 0.00006708180428412224 × (1) = 0.00006708180428412224

(-1) × 0.00006708368307282375 × (-1) = 0.00006708368307282375

(-1) × 0.000038798961684839866 × (-1) = 0.000038798961684839866

(1) × 0.00003879673170860536 × (1) = 0.00003879673170860536

0.00006708180428412224 + 0.00006708368307282375 + 0.000038798961684839866 + 0.00003879673170860536

= 0.0002117611807503912 m²

Convert variance to standard deviation

σ_DD = √(0.0002117611807503912)

= 0.014552016380914062 m

Squaring the signs makes all four independent variance contributions positive.

Off-diagonal covariance: shared G09

R_{02}=\operatorname{Cov}(DD_{04,\Phi},DD_{07,\Phi})=\sum_kD_{0k}\sigma_k^2D_{2k}

Only the two shared reference columns overlap

Base G09: (−1) × 0.000038798961684839866 × (−1) = 0.000038798961684839866

Rover G09: (+1) × 0.00003879673170860536 × (+1) = 0.00003879673170860536

0.000038798961684839866 + 0.00003879673170860536

= 0.00007759569339344522 m²

Correlation coefficient

0.00007759569339344522 / √(0.0002117611807503912 × 0.00016440829112487126)

= 0.41586520580049174

Reusing the reference observations couples DD errors even when their target measurements are independent.

DD covariance heatmap

R_DD (m²); color intensity encodes correlation8 × 8

Row / column	G04 phase	G04 code	G07 phase	G07 code	G16 phase	G16 code	G30 phase	G30 code
G04 phase
G04 code
G07 phase
G07 code
G16 phase
G16 code
G30 phase
G30 code

G04 phase × G04 phase = 0.0002117611807503912 m² · correlation = 1.000000

Show full mathematical matrix

Phase and code alternate. Nonzero off-diagonals form same-observable blocks; phase/code cross-covariances are zero under this model. Select a cell for exact covariance and correlation.

Undifferenced covariance diagonal

R_UD diagonal: all 20 exported variances in m²
UD observation	Variance
phase:base:G04	0.0000670818
phase:rover:G04	0.0000670837
phase:base:G09	0.0000387990
phase:rover:G09	0.0000387967
code:base:G04	0.9130578916
code:rover:G04	0.9130834640
code:base:G09	0.5280969785
code:rover:G09	0.5280666260
phase:base:G07	0.0000434081
phase:rover:G07	0.0000434045
code:base:G07	0.5908320790
code:rover:G07	0.5907838345
phase:base:G16	0.0003200092
phase:rover:G16	0.0003198158
code:base:G16	4.3556813318
code:rover:G16	4.3530484264
phase:base:G30	0.0000998609
phase:rover:G30	0.0000998644
code:base:G30	1.3592179401
code:rover:G30	1.3592659944

Show full 20 × 20 R_UD

Ready for weighted least squares

R_DD has positive eigenvalues (minimum 0.00010396339754041666 m²). Cholesky can now turn this correlated system into independent unit-variance observations.

8. Float solution (weighted least squares)

Solve the 8-observation, 7-unknown system without forming a matrix inverse or squaring its condition number through normal equations.

Cholesky → whitening → QR → solution

y=Hx+\epsilon,\quad R_{DD}=LL^T,\quad LH_w=H,\quad Ly_w=y

H_w=Q R,\quad R x=Q^Ty_w,\quad \min_x\|y_w-H_wx\|^2

Cholesky handles measurement correlation; triangular solves whiten H and y. Reduced QR supplies an upper-triangular R and the transformed right-hand side. Back substitution yields the baseline and float ambiguities.

1. Verify the Cholesky factors

L_{ii}=\sqrt{R_{ii}-\sum_{k<i}L_{ik}^2},\quad L_{ji}={R_{ji}-\sum_{k<i}L_{jk}L_{ik}\over L_{ii}}

First diagonal and a shared-reference entry

L00 = √(0.0002117611807503912) = 0.014552016380914063

L20 = 0.00007759569339344522 / 0.014552016380914063 = 0.005332298381358141

L22 = √(0.00016440829112487126 − (0.005332298381358141)²) = 0.011660826947392565

Verify R02 = L00 × L20 = 0.014552016380914063 × 0.005332298381358141 = 0.00007759569339344522

The positive off-diagonal covariance creates a nonzero lower-triangular coefficient.

Show full L (8 × 8)

2. Whiten by triangular solve

Solve L · y_w = y row by row, and L · H_w = H for each column. The inverse notation is explanatory; no inverse is computed.

First phase observation and first H coefficient

y_w0 = -188.3481102883816 / 0.014552016380914063 = -12943.09361384534

H_w00 = 0.375041121025931 / 0.014552016380914063 = 25.772450443212964

y_w2 = (52.43941770866513 − 0.005332298381358141 × (-12943.09361384534)) / 0.011660826947392565

= 10415.715401959364

A later phase row subtracts the shared-reference contribution already carried by earlier rows. Whitening is more than dividing every DD by its own sigma.

Original and whitened system

Original H and y

Show full whitened H_w (8 × 7) and y_w (8 × 1)

3. QR and back substitution

Reduced Q is 8 × 7, R is 7 × 7, and b = Qᵀy_w is 7 × 1. The unused orthogonal residual direction accounts for one degree of freedom.

One Qᵀ row-by-column product

Q(1,1) × y_w1 = (-0.342730699188212) × (-12943.09361384534) = 4435.995523931694

Q(2,1) × y_w2 = (-0.002937691707327533) × (-112.35085130445847) = 0.3300521641882964

Q(3,1) × y_w3 = (0.5737978481463599) × (10415.715401959364) = 5976.515084549181

Q(4,1) × y_w4 = (0.004918267269825945) × (88.92456611407272) = 0.43735478300231717

Q(5,1) × y_w5 = (-0.454603687891473) × (18710.946577367285) = -8506.065318011502

Q(6,1) × y_w6 = (-0.003896603039069769) × (159.15062306634192) = -0.6201468015101552

Q(7,1) × y_w7 = (0.5886881349760987) × (-12589.85232486183) = -7411.496684747411

Q(8,1) × y_w8 = (0.0050458982997951315) × (-107.27341045280696) = -0.5412907194170439

= -5505.445424851773

Show R_QR, Q, and actual Qᵀy_w

Numerical back substitution: every unknown, last row first

Float solution: baseline + ambiguities

Real float state and unscaled formal 1σ from the trace
Unknown	Float value	Formal σ	Unit
dX	329.964491207	2.554462048	m
dY	331.829789367	2.903071546	m
dZ	-661.716860707	3.380692046	m
N04	12.449682739	6.668804648	cycles
N07	6.983324467	3.907815928	cycles
N16	25.001732843	16.319926290	cycles
N30	0.480403785	10.195412531	cycles

Exact solution vector

Conditioning and degrees of freedom

κ(H_w)1636.8

8 − 7 observations/unknowns1 DOF

Warning · Weak conditioning

1636.8 indicates substantially weaker conditioning than the candidate-fixed geometry-only problem. With only one residual degree of freedom, single-epoch fit diagnostics are limited.

Float means the four ambiguity states are real numbers, not yet constrained to integers.

4. Formal state covariance

Q_{xx}=R_{QR}^{-1}R_{QR}^{-T},\qquad \sigma_{x_i}=\sqrt{Q_{xx,ii}}

Compute the factors with triangular solves. This is covariance under the declared measurement model, without rescaling by the tiny post-fit variance factor.

Construct Q_xx from triangular covariance factors

First baseline sigma

Q_xx00 = 6.525276353430393 m²

σX = √(6.525276353430393)

= 2.554462047756904 m

Full Q_xx (7 × 7), including baseline–ambiguity correlations

The lower-right 4 × 4 block Q_NN is the ambiguity covariance passed into LAMBDA.

5. Post-fit residuals

The numerical values are available in the accompanying tables. The interactive chart loads as you approach it.

Why phase residuals are nearly zero

v=y-H\hat x,\qquad v_w=L^{-1}v

G04 code residual

-190.742000002414 − ((0.375041121025931) × (329.96449120681626) + (-0.25445239752265914) × (331.82978936717865) + (0.34762991093087486) × (-661.7168607071129) + (0) × (12.449682739387208) + (0) × (6.983324467405504) + (0) × (25.00173284312381) + (0) × (0.48040378538633366))

= -0.02479386037992981 m

Warning · Not perfect carrier-phase data

Each of the four phase DD observations has its own ambiguity parameter. Those parameters can absorb phase differences at a single epoch, leaving phase residuals near zero. This does not prove flawless phase measurements.

9. LAMBDA ambiguity fixing

Find the closest integer vector under the full correlated ambiguity covariance, rather than rounding each component independently.

Float ambiguity vector

\begin{bmatrix}12.4496827394\\6.9833244674\\25.0017328431\\0.4804037854\end{bmatrix}

Ordering: N04, N07, N16, N30, all relative to G09. These real-valued estimates include substantial covariance and cannot be judged one component at a time.

Q_NN: real ambiguity covariance (cycles²)4 × 4

Row / column	N04	N07	N16	N30
N04
N07
N16
N30

N04 × N04 = 44.472955436810025

Show full mathematical matrix

Integer least squares versus naive rounding

\hat n=\arg\min_{n\in\mathbb Z^4}(a-n)^TQ_{NN}^{-1}(a-n)

Full four-dimensional candidate comparison
Candidate	Integer vector	Squared norm
Best LAMBDA	[13, 7, 27, 1]	0.067325666147
Second	[11, 7, 22, -1]	0.097304202353
Naive rounding	[12, 7, 25, 0]	0.163160530559

Comparison · Rounding is feasible, but inferior

[12, 7, 25, 0] is an integer vector, but it does not minimize this correlated metric. The best vector is [13, 7, 27, 1].

Naive squared norm: derived from exported Q_NN

This value is derived from the package inputs, not an additional raw exported result. Solve Q_NN · w = (a − n_round), then dot the error with w.

Weighted four-dimensional distance

e = a − n_round = [0.4496827393872085, -0.016675532594495834, 0.0017328431238112785, 0.48040378538633366]

Solve Q_NN w = e → w = [2.3866232595980623, 4.946024743825871, -0.292327924687441, -1.7216321448728586]

e1 × w1 = (0.4496827393872085) × (2.3866232595980623) = 1.0732232852612855

e2 × w2 = (-0.016675532594495834) × (4.946024743825871) = -0.08247759682885121

e3 × w3 = (0.0017328431238112785) × (-0.292327924687441) = -0.0005065584341926535

e4 × w4 = (0.48040378538633366) × (-1.7216321448728586) = -0.8270785994397141

Sum = 1.0732232852612855 + -0.08247759682885121 + -0.0005065584341926535 + -0.8270785994397141

= 0.1631605305585275

The full correlated metric makes this larger than both exported LAMBDA candidate norms. No explicit covariance inverse is needed.

LAMBDA integer transformation

Z: integer unimodular matrix4 × 4

Row / column	1	2	3	4
N04
N07
N16
N30

N04 × column 1 = 17

Show full mathematical matrix

det(Z) = +1. An integer unimodular transformation preserves the integer lattice while reducing covariance correlation for search.

z=Z^Ta,\qquad Q_z=Z^TQ_{NN}Z

First transformed ambiguity: column 1 of Z × a

(17) × (12.449682739387208) = 211.64460656958255

(35) × (6.983324467405504) = 244.41635635919263

(-2) × (25.00173284312381) = -50.00346568624762

(-12) × (0.48040378538633366) = -5.764845424636004

= 400.2926518178916 cycles

All four z entries and Q_z

Candidates in ambiguity space

The full covariance sets the scale. The close view separates the actual candidates without distorting either axis.

01 / Covariance overviewMarginal covariance · 1σ

02 / Candidate detailInteger lattice · 1 cycle

Best LAMBDA[N04, N07, N16, N30] = [13, 7, 27, 1] cyclesFull 4D squared norm = 0.06732566614671152

Equal axis scales. The dashed ellipse is a 1σ contour of the N04/N16 marginal covariance, not a validation boundary or the full 4D search region. A shorter distance in this projection does not determine the best candidate; ranking uses the complete Q_NN.

Independent rounding can choose a poorer candidate because the objective depends on all four ambiguities and their correlations.

Search and back-transform

The Qelaro lambda_method path decorrelates the float ambiguities and performs integer search in z-space. The covariance metric guides enumeration and pruning; it is not a nearest-point search in a two-dimensional picture.

Exported candidates transformed by Zᵀ; derived integer vectors
Candidate	z candidate	Original n	Norm
Best	[400, -157, -164, -190]	[13, 7, 27, 1]	0.067325666147
Second	[400, -158, -165, -191]	[11, 7, 22, -1]	0.097304202353

Verify the best and second full 4D squared norms numerically

n=Z^{-T}z_{integer}

Show integer back-transform and first row product

For the full algorithmic decorrelation and search walkthrough, continue to the LAMBDA Method tutorial →.

Ratio diagnostic

\mathrm{ratio}={s_2\over s_1}

Second / best squared norm

0.09730420235313753 / 0.06732566614671152

= 1.445276488480616

Warning · Candidate-fixed; not validated

LAMBDA produced a best candidate. The ratio 1.445276488480616 is a diagnostic, not an accepted operational fix. A single epoch does not validate ambiguity reliability, continuity, or multi-epoch post-fix behavior. Next, constrain the geometry to this candidate for teaching.

10. Candidate-fixed solution and baseline

Hold [N04,N07,N16,N30] = [13,7,27,1] and estimate the three baseline components using the same full R_DD.

Remove integer phase contributions

y_{fixed}=y-H_Nn,\qquad y_{fixed}=H_{xyz}b+\epsilon

G04 phase adjusted observation

-188.3481102883816 − (0.19029367279836487 × 13)

0.19029367279836487 × 13 = 2.473817746378743

= -190.82192803476036 m

Code rows do not change. Only the phase ambiguity contributions are subtracted; all eight observations and their covariance remain.

Observations8

Free unknowns: dX,dY,dZ3

Degrees of freedom: 8 − 35

Whiten and QR the three-column system

L H_{xyz,w}=H_{xyz},\quad L y_{fixed,w}=y_{fixed},\quad R_3 b=Q_3^Ty_{fixed,w}

First candidate-fixed whitening and last QR row

y_fixed,w0 = -190.82192803476036 / 0.014552016380914063 = -13113.09189323316

R33 × dZ = b3

(34.51096719520549) × dZ = -22843.033070292375

dZ = -22843.033070292375 / 34.51096719520549

= -661.9064873228442 m

Full H_xyz,w (8 × 3), whitened y, QR R and right-hand side

Numerical back substitution: every unknown, last row first

Baseline vector: rover → base

Candidate-fixed ECEF baseline and formal conditional standard deviations
Axis	Baseline (m)	Formal σ (m)
X	329.63936360562	0.02189458470
Y	331.51136991397	0.02488255636
Z	-661.90648732284	0.02897629598

Baseline length

‖b‖ = √((329.63936360562184)² + (331.5113699139712)² + (-661.9064873228442)²)

dX² = 108662.11003831937 m²

dY² = 109899.78838223784 m²

dZ² = 438120.19796006655 m²

Sum = 656682.0963806238 m²

= 810.3592390912956 m

Base ECEF: fixed rover + baseline

Coordinate addition in metres; base is estimated, rover held fixed
Axis	Fixed rover	+ baseline	= estimated base
X	3246374.5864	329.639363606	3246704.225763605
Y	4052665.3939	331.511369914	4052996.905269914
Z	3693176.0609	-661.906487323	3692514.154412677

Improved conditioning

κ(H_w): 1636.8 → 3.01124. Ratio = 543.567×. Constraining four ambiguity parameters strengthens the geometry problem; it does not validate those integers.

Formal covariance and baseline-length uncertainty

Q_b=R_3^{-1}R_3^{-T},\quad g={b\over\|b\|},\quad \sigma_{\|b\|}^2=g^TQ_bg

Candidate-fixed baseline covariance Q_b (m²)3 × 3

Row / column	X	Y	Z
X
Y
Z

X × X = 0.0004793728393903604

Show full mathematical matrix

Length variance: full covariance propagation

g = [0.4067817675223475, 0.409091861883027, -0.8168062451723012]

Q_b g = [0.0003660057369139119, 0.00038247516328917267, -0.0006597058197951598]

g1 × (Q_b g)1 = 0.4067817675223475 × 0.0003660057369139119 = 0.0001488844605851604

g2 × (Q_b g)2 = 0.409091861883027 × 0.00038247516328917267 = 0.0001564674766739824

g3 × (Q_b g)3 = -0.8168062451723012 × -0.0006597058197951598 = 0.0005388518335851993

Variance = 0.0008442037708443421 m²

σlength = √(0.0008442037708443421)

= 0.029055184921874824 m

Warning · Formal precision ≠ external accuracy

These centimetre-level formal conditional sigmas assume fixed rover approximate coordinates, the selected R_DD, and the chosen exact integer candidate. They exclude rover coordinate uncertainty, absolute base truth error, atmosphere/model error, orbit model error, multipath and wrong-candidate risk. This is not a claim of centimetre external accuracy.

11. Results and statistical diagnostics

A fit must be interpreted against its stochastic assumptions. Neither small residuals nor small formal covariance prove that the position or integer candidate is correct.

Post-fit double-difference residuals

The numerical values are available in the accompanying tables. The interactive chart loads as you approach it.

All exact float and candidate-fixed residuals

Formal conditional precision

The numerical values are available in the accompanying tables. The interactive chart loads as you approach it.

1σ coordinate precision in metres, conditional on the model
Axis	Float	Candidate-fixed
X	2.554462	0.021895
Y	2.903072	0.024883
Z	3.380692	0.028976

Chi-square: whitened residual energy

\chi^2=v^T R_{DD}^{-1}v=v_w^Tv_w,\qquad \nu=n-p

A two-sided 95% consistency interval uses χ² quantiles at 0.025 and 0.975 for ν degrees of freedom. Too little residual energy is also inconsistent with the declared model.

Numerical sum of squared whitened residuals

Float chi-square · DOF 1

χ² = 0.0004832536841015236 · DOF = 1. Interval [0.0009820691171752555, 5.02388618731489]. Two-sided consistency: below lower bound.

Warning · Below lower bound

0.0004832536841015236 < 0.0009820691171752555. This does not pass the two-sided test.

Candidate-fixed chi-square · DOF 5

χ² = 0.06780891982717298 · DOF = 5. Interval [0.8312116134866626, 12.832501994030023]. Two-sided consistency: below lower bound.

Warning · Below lower bound

0.06780891982717298 < 0.8312116134866626. This also does not pass the two-sided test.

Warning · Why unusually small residual energy matters

Possible explanations include an over-conservative noise model, single-epoch parameter flexibility, or mismatched stochastic assumptions. The float phase ambiguity parameters absorb the four phase observations. An unusually small χ² does not prove superior accuracy, and the candidate remains unvalidated.

12. What this lesson proves — and what it does not

The full numerical path is reproducible. Its conclusions are conditional on a deliberately small, one-epoch teaching model.

This lesson demonstrates

✓Real one-epoch RINEX processing
✓Broadcast NAV satellite propagation
✓Common satellite and reference selection
✓Satellite-to-satellite single differences
✓Between-receiver double differences
✓H construction with fixed LOS vectors
✓D construction and covariance propagation
✓Correlated R_DD and Cholesky whitening
✓QR weighted least squares and float ambiguities
✓LAMBDA candidate search
✓Constrained candidate-fixed baseline
✓Residual and chi-square diagnostics

These claims are not proven

×Validated multi-epoch PPK performance
×Production ambiguity fixing reliability
×Cycle-slip management over time
×Reference-change continuity
×External absolute coordinate accuracy
×Full atmospheric / antenna correction validation
×Equivalence to a complete production RTKLIB solution
×Centimetre external accuracy

Carry the evidence forward

Keep the observation signs, covariance correlations, baseline direction and validation status together. The strongest result of this lesson is that you can inspect every stage from raw code/phase through to the conditional coordinate covariance.

Resources and numerical provenance

Use the immutable exports to check the numbers independently, then continue through the related tutorials.

Download the actual lesson data

Complete teaching package: data, Python reproducer and citationopengnsslab-dd-teaching-v1.zip ↓Reproduction steps, conventions and scopeREADME.md ↓File checksums and original export provenancemanifest.json ↓All normalized lesson data, RINEX snippets and derived arithmeticlesson.json ↓Four real DD rowsdouble_differences.csv ↓H and observation vector ydesign_matrix_and_observations.csv ↓155 NAV orbit samplesorbit_track_samples.csv ↓

Source provenance and implementation notes

SPP from Scratch →Single-receiver modelling and weighted least squares.LAMBDA Method →Detailed integer decorrelation and search.DOP Calculation →Geometry and covariance interpretation.Broadcast Ephemeris →NAV-to-ECEF numerical propagation.

What is a GNSS double difference?

What you’ll build in this tutorial

One epoch, fully traced

Double-difference processing pipeline

Inputs and sign convention

Where this lesson ends

1. Problem and model

Code pseudorange

Carrier phase in metres

Read the observation equation

2. Input data (RINEX)

Base observations · 36.24O

Rover observations · 15.24O

Convert cycles into metres

Base G04 phase

Base navigation · 36.24N

Rover navigation · 15.24N

NAV → transmit time → ECEF → LOS

Epoch satellite state

NAV orbit tracks

A numeric line of sight

G04: three subtractions, a norm, three divisions

3. Satellite selection

Common-satellite filtering

Base · 7 satellites

∩ Rover · 5 satellites

= Common · 5 satellites

Choose the highest common base elevation

Maximum over the intersection

Skyplot / reference selection

4. Single differences from scratch

Same receiver, two satellites

Code worked example: G04 − G09

Base code SD

Rover code SD

Phase worked example: cycles first

Base phase SD

Rover phase SD

All single differences in the trace

Single-difference observation view

What one subtraction achieves

5. Double differences and differencing matrix

Double-difference definition

Worked G04 double differences

Code DD

Phase DD

All real DD observations

DD observation plots

Differencing matrix D · exact 8 × 20 structure

G04 phase: row-by-column product

6. Linear design matrix

Freeze the LOS, then construct H

Row-by-row geometry walkthrough

−uG04 + uG09

Phase row versus code row

Full design matrix and observation vector

One complete H row × x example

G04 phase prediction using the float state

7. Stochastic model and DD covariance

Elevation-dependent observation variance

Base G04 phase weighting

From independent UD to correlated DD

Worked diagonal: G04 phase variance

Four nonzero products

Convert variance to standard deviation

Off-diagonal covariance: shared G09

Only the two shared reference columns overlap

Correlation coefficient

DD covariance heatmap

Undifferenced covariance diagonal

8. Float solution (weighted least squares)

Cholesky → whitening → QR → solution

1. Verify the Cholesky factors

First diagonal and a shared-reference entry

2. Whiten by triangular solve

First phase observation and first H coefficient

Original and whitened system

3. QR and back substitution

One Qᵀ row-by-column product

Float solution: baseline + ambiguities