Vallado & Crawford (2008): SGP4 Orbit Determination

Source

Vallado, D.A. and Crawford, P., “SGP4 Orbit Determination,” AIAA 2008-6770, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, Hawaii, August 18-21, 2008. Vallado: Analytical Graphics Inc. / Center for Space Standards and Innovation. Crawford: University of Dundee.

Download PDF

Why This Matters for SGP4

Every other document in this archive addresses the forward problem: given a TLE, compute the satellite’s position and velocity. This paper addresses the inverse problem: given observations or an externally generated ephemeris, produce a TLE that best fits the data.

The relationship is not symmetric — the inverse is harder. SGP4’s analytical perturbation model includes secular, long-period, and short-period corrections derived from Brouwer’s 1959 theory, atmospheric drag from STR#2, and deep-space resonance from STR#1. Deriving analytical partial derivatives through all of this is impractical. Vallado and Crawford sidestep the problem with finite differencing, treating SGP4 as a black box.

The compatibility constraint is absolute: STR#3 itself stated that “the NORAD element sets must be used with one of the models described in this report.” A TLE generated by inverting SGP4 is only valid when propagated forward with the same SGP4. This paper completes the ecosystem — forward and inverse, using the identical model.

Brouwer (1959) --> Lane & Hoots (STR#2) --> Hoots & Roehrich (STR#3, 1980)
                                                  |
                                                  v
                                        Vallado et al. (2006/2008, Rev-1)
                                                  |    Forward: TLE --> state
                                                  |
                                                  v
                                        Vallado & Crawford (2008) <-- THIS PAPER
                                                       Inverse: state --> TLE
                                                       (differential correction
                                                        using SGP4 as the
                                                        propagation model)

The Differential Correction Approach

The method is batch weighted least-squares, the standard orbit determination technique adapted for TLE fitting. Observations are accumulated across a fit span, residuals are computed against the SGP4-predicted positions, and the orbital elements are iteratively corrected until the residuals converge.

The core update equation:

$\delta x = (A^T W A)^{-1} A^T W b$

where:

Symbol	Meaning
$\delta x$	Correction to the state vector (orbital elements)
$A$	Jacobian matrix — partials of observations with respect to state
$W$	Weighting matrix — diagonal, $1/\sigma^2$ per observation
$b$	Residual vector — observed minus computed
$(A^T W A)^{-1}$	Covariance matrix — uncertainty estimate as a byproduct

The iteration continues until three simultaneous criteria are met: the relative change in RMS drops below $\epsilon = 0.0002$, the absolute RMS drops below $\epsilon$, and the iteration count stays under 25. Physical bounds are enforced at each step — the semimajor axis cannot drop below 1 Earth radius, eccentricity cannot exceed 1.0, and mean motion cannot go negative.

Two matrix solution methods are implemented: LU decomposition (efficient but fragile near singularities) and SVD (robust, recommended for operational use). The matrix sizes are small (6x6 or 7x7), so the SVD overhead is negligible.

State Representations

The choice of state variables determines where the singularities fall and how fast the iteration converges. The code supports three options:

Keplerian elements $(a, e, i, \Omega, \omega, M)$ — the most intuitive representation since these are the TLE variables directly. But circular orbits ($e = 0$) make $\omega$ undefined, and equatorial orbits ($i = 0$) make $\Omega$ undefined. These singularities cause the Jacobian matrix to become ill-conditioned.

Equinoctial elements $(a, k_e, h_e, p_e, q_e, \lambda_e)$ — the preferred representation. The singularities are pushed to retrograde equatorial orbits ($i = 180\degree$), which are not physically encountered in the satellite catalog. The definitions:

$k_e = e \cos(\Omega + \omega), \quad h_e = e \sin(\Omega + \omega)$

$p_e = \tan(i/2) \sin \Omega, \quad q_e = \tan(i/2) \cos \Omega$

$\lambda_e = M + \omega + \Omega$

The inverse transformation recovers the Keplerian elements needed for TLE output:

$e = \sqrt{k_e^2 + h_e^2}, \quad i = 2\arctan\sqrt{p_e^2 + q_e^2}, \quad \Omega = \arctan\frac{p_e}{q_e}$

Cartesian state vectors $(x, y, z, \dot{x}, \dot{y}, \dot{z})$ — no geometric singularities, but all six elements change rapidly over an orbit, requiring more iterations to converge.

Finite Differencing for the Jacobian

SGP4’s perturbation model is a layered stack of Brouwer corrections, Lyddane reformulations, atmospheric drag terms, and deep-space resonance integration. Deriving analytical partial derivatives through this entire model would be a substantial mathematical undertaking — and would need to be re-derived every time the propagator changes.

Instead, the code treats SGP4 as a black box. For each element of the state vector, it:

1. Perturbs the element by 0.1% of its current value (floored at $10^{-7}$)
2. Re-initializes SGP4 with the perturbed state
3. Propagates to the observation time
4. Forms the numerical derivative from the difference

This costs 6 extra SGP4 propagations per observation per iteration (7 if solving for $B^*$). For a typical fit with 144 observations and 10 iterations, that is roughly 10,000 SGP4 calls — trivial on modern hardware, but it explains why the original AFSPC system used a simpler propagator (SGP, not SGP4) for operational orbit determination when computational cost mattered.

Why Not Analytical Partials?

Analytical partials would be faster per evaluation but brittle to maintain. Every code correction in Rev-1’s 17 changes or mode switch between AFSPC and improved operation would require re-deriving the partials. Finite differencing works with any SGP4 variant unchanged.

The Kozai-Brouwer Round Trip

TLEs publish mean motion in Kozai’s formulation, but SGP4 internally converts to Brouwer’s formulation during initialization. The forward direction (Kozai $n_0$ to Brouwer $a_0''$) is handled by the existing sgp4init code from STR#3.

The inverse direction — recovering the Kozai mean motion from a Brouwer semimajor axis after the DC process has converged — is not provided by the forward code. This paper adds a Newton-Raphson root solver:

$g(n_0) = f(n_0) - a_0'' = 0$

where $f(n_0)$ is the Kozai-to-Brouwer conversion from Hoots & Roehrich (1980, p. 10). This typically converges in 2-4 iterations to machine precision. Without this step, the TLE output would encode incorrect mean motion values, and forward propagation with the standard SGP4 would produce the wrong trajectory.

Step Limiting for Robustness

Unconstrained Newton-Raphson corrections can diverge, especially in early iterations when the initial guess is poor. The paper implements a tiered correction limiter based on the ratio of correction to nominal value:

Correction/Nominal Ratio	Action
> 1000	Limit correction to 10% of nominal
> 200	Limit correction to 30% of nominal
> 100	Limit correction to 70% of nominal
> 10	Limit correction to 90% of nominal

The authors acknowledge these thresholds are “somewhat arbitrary, and not fully researched.” Two different limiting strategies each fail on a different set of ~20-30 satellites out of the full catalog, with very little overlap between failure sets. This suggests an adaptive strategy could solve nearly the entire catalog.

Catalog-Scale Results

The acid test: processing the entire unclassified satellite catalog (16,235 objects) through a TLE-to-ephemeris-to-TLE round trip. The baseline configuration recovered 16,198 of 16,235 objects (99.6%) to sub-meter accuracy — the regenerated TLE produces positions within 1 meter of the original across the fit span.

Performance across orbit regimes using precision ephemerides:

Orbit	Example	Fit Accuracy	Prediction (2 days)
LEO	ICESat	~500-600 m	~4-5 km
MEO	GPS	~150 m	~10 km (30 days)
GEO	Intelsat	Variable	Maneuver-dependent

The ephemeris-derived TLEs outperformed the AFSPC catalog TLEs, confirming Kelso’s (2007) independent finding that higher-quality observations yield better TLEs even through SGP4’s simplified force model.

GEO Low-Inclination Objects

GEO satellites with inclination < 0.1 degrees stress the singularity handling in both Keplerian and equinoctial element formulations and may need special treatment: 6-state solution without $B^*$, zeroed lunar-solar perturbation terms at initialization, and careful monitoring of the Jacobian condition number. Nearly all of the ~30 problem cases in the full catalog fall in this category.

Covariance as a Byproduct

The $(A^T W A)^{-1}$ matrix from the least-squares solution is the formal covariance of the fitted state. This is significant because AFSPC does not distribute covariance data with TLEs — the standard two-line format has no covariance fields.

Representative semimajor axis uncertainties from ephemeris fits:

Object	$\sigma_a$ (meters)
ICESat (LEO)	0.6
GPS (MEO)	3.4
GEO satellite	~6.9

Covariance Interpretation

The authors caution that these values “seem to be more a measure of the frequency of the data than the accuracy.” Batch least-squares covariance reflects how well the data constrains the fit, not the absolute accuracy of the result. For conjunction screening, the covariance requires careful calibration against independent truth data.

Observation Processing

The code supports three observation types:

Type	Input	Status at Publication
TLE values	Direct TLE element extraction (test mode)	Operational
Ephemeris	Position/velocity in TEME, ECEF, or ECI	Operational
Raw observations	Range, azimuth, elevation	Incomplete

Raw Observation Processing

The raw observation pipeline (radar/optical data directly to TLE) was not fully functional at publication. The initial orbit determination front-end did not converge reliably. The authors noted this would be updated in the web distribution or a follow-on paper. The ephemeris input mode is the primary use case — organizations with independent tracking data would typically process observations to an ephemeris first, then feed that into the SGP4 DC.

Practical Significance

This paper fills two gaps in the SGP4 ecosystem:

1. Independent TLE generation. Any organization with tracking data (optical telescopes, radars, or access to precision ephemerides) can now produce TLEs that are guaranteed compatible with the standard SGP4 forward propagator.

2. Post-maneuver recovery. When a satellite maneuvers, the pre-maneuver TLE becomes stale. Organizations with independent observations can generate a post-maneuver TLE faster than waiting for the next catalog update.

3. Covariance for conjunction assessment. The covariance matrix output enables uncertainty quantification for collision probability calculations — something the standard TLE format cannot provide.

The same two authors (Vallado and Crawford) who produced the definitive forward propagator in Rev-1 now provide its inverse. The shared authorship guarantees that the forward and inverse pipelines are built on identical code, eliminating the compatibility risk that plagued SGP4 for 25 years.

---

Rev-1: The Forward Propagator The forward propagation code being inverted by this paper

Vallado Rev-3 (2013) The latest SGP4 code with dual AFSPC/improved operation modes

Crawford (1995): Kepler Solver Crawford's bounded Newton-Raphson technique, from the same co-author

Brouwer (1959): Mean Element Theory The mean element theory that the differential correction process fits