Basic orbit determination

1 - Introduction

The objective here is to solve the following problem:

Find X0 such that:

F(X0, t) = M(t) for t = t1, t2, ... tn

Where:

X0 = X(t0) is the state vector at (initial) time t0 (size = P x 1). It may contain initial orbital elements plus additional parameters (such as drag coefficient ...).
M(t) is the measurements vector at time t (size = K x 1).
F is a function that propagates the state vector from time t0 to time t and computes the theoretical measurements vector at time t.

In general the above equation has no exact solution and is solved in the "least squares" sense such that:

Σ ∥M(tk) - F(X0, tk)∥² is minimum. Note that in this case, all measurements components have the same weight, which is not necessarily the case.

One common way to proceed is to linearize the equations around a "reference" trajectory, that is for a particular value of X0. After linearization, the equations can be written:

H * dX0 = dZ

Where:

dX0 denotes state vector increments (size = P x 1),
dZ is a one-column vector containing M(tk) - F(X0, tk) for all tk (size = Kn x 1),
H is the sensitivity matrix: matrix of the partial derivatives of the measurements with respect to the initial state (size = Kn x P).

dX0 can then be computed using a least squares method. A few iterations are necessary to reach convergence provided the initial guess is not too far from the solution.

First execute the code snippet below that defines some utility functions.

// ================================================
// Utility functions
// ================================================

// ------------------------------------------------
// Plot final results (position gaps in local frame) 
// cjd: dates/times (TREF time scale) (1xN)
// pos: position vectors to be checked (ECI) (3xN)
// pvref: reference position/velocity vectors (ECI) (6xN)
// ------------------------------------------------
function plot_gaps_qsw(cjd, pos, pvref); 
  posref = pvref(1:3,:); 
  velref = pvref(4:6,:); 
  dpos_qsw = CL_fr_inertial2qsw(posref, velref, pos - posref); 

  scf(); 
  plot(cjd - cjd(1), dpos_qsw(1,:), "r"); 
  plot(cjd - cjd(1), dpos_qsw(2,:), "g"); 
  plot(cjd - cjd(1), dpos_qsw(3,:), "b"); 
  xtitle("Position gaps in qsw local frame (m)", ...
         "Time from initial date (days)");
  CL_g_legend(gca(), ["q", "s", "w"]); 
  CL_g_stdaxes(); 
endfunction

// Compute a reference trajectory using "lydlp" model 
// cjd (1xN): dates (TREF time scale)
// pvref (6xN): reference position/velocity (ECI frame) 
function [cjd, pvref]=build_trajectory()
  cjd0 = CL_dat_cal2cjd(2018,1,1); 
  mean_cir0 = [7.e6; 0; 1.e-3; 1.7; 0; 0]; 

  cjd = cjd0 + (0 : 60 : 2 * 86400) / 86400; 
  cir = CL_ex_propagate("lydlp", "cir", cjd0, mean_cir0, cjd, "o"); 
  pvref = CL_oe_convert("cir", "pv", cir); 
endfunction

2 - Application 1: Eckstein-Hechler versus Lyddane

We want to compare two propagation models: Eckstein-Hechler ("eckhech") and Lyddane ("lydlp"). The model to be adjusted (or rather the associated initial conditions) is Eckstein-Hechler, whereas the reference trajectory is computed using Lyddane.

The measurements are the position vectors at each propagation date.

The initial state vector (X0) consists of the mean (equinoctial) orbital elements (size = 6x1).

The main function is orb_fit that computes the initial mean orbital elements for "eckhech" model best fitting the reference trajectory.

// ================================================
// Application 1
// ================================================

// ------------------------------------------------
// Measurements function
// ------------------------------------------------
function [pos]=compute_pos(cjd0, equin0, cjd)
  // only 5 zonal terms (same as Lyddane)
  J1JN = CL_dataGet("j1jn"); 
  equin = CL_ex_propagate("eckhech", "equin", cjd0, equin0, cjd, "o", ...
                          j1jn = J1JN(1:5)); 
  pv = CL_oe_convert("equin", "pv", equin); 
  pos = pv(1:3,:); 
endfunction
 
// ------------------------------------------------
// Determination of initial orbital elements 
// cjd: dates/times (TREF time scale) (1xN)
// pvref: reference position/velocity vectors (ECI) (6xN)
// ------------------------------------------------
function [cjd0, equin0]=orb_fit(cjd, pvref)
  // Definition of initial state (first guess: osculating elements) 
  cjd0 = cjd(1); 
  equin0 = CL_oe_convert("pv", "equin", pvref(:,1));

  // Fitting parameters: default values except for: 
  // - svScale: state vector scale factors
  // - verbose: verbose mode 
  params = struct(); 
  params.svScale = [equin0(1); 1 * ones(5,1)]; 
  params.verbose = 1; 

  // Adjustment 
  equin0 = CL_adjustState(cjd0, equin0, cjd, pvref(1:3,:), compute_pos, params); 
endfunction

 
// ================================================
// MAIN
// ================================================
// Reference trajectory (using "Lyddane" model)
[cjd, pvref] = build_trajectory(); 

// Determination of initial state (date + mean equinoctial elements)   
[cjd0, equin0] = orb_fit(cjd, pvref); 

// Plot results 
pos = compute_pos(cjd0, equin0, cjd); 
plot_gaps_qsw(cjd, pos, pvref);

This example enables us to compare the two models ("eckhech" and "lydlp") in the "LEO" case, when the same physical constants are used, and after adjustment of the mean elements:

There is (almost) no long-term relative drift.
There are short period variations (~60 meters) mostly in the radial direction.
There is a slight offset between the mean semi major-axes: about 7 meters. It could be the cause of discrepancies when propagating if this offset is not taken into account: about 1 km along-track per day.

3 - Application 2: TLE estimation

We now apply the generic estimation functions to TLEs.

The initial state vector (X0) is a vector of 7 elements: 6 orbital (equinoctial) elements + bstar. The measurements consist of the position vectors at each propagation date.

The main function is tle_fit that computes a TLE from a reference trajectory defined by cjd (time stamps relative to TREF time scale) and pvref (reference position and velocity vectors in ECI frame).

// ================================================
// Application 2: TLE estimation
// ================================================

// ------------------------------------------------
// State vector + associated time stamp to TLE
// state: 7xN; cjd: 1xN
// ------------------------------------------------
function [tle]=state2tle(cjd, state)
  tle = CL_tle_new(length(cjd)); 
  tle = CL_tle_setElem(tle, "equin", cjd, state(1:6,:)); 
  tle = CL_tle_set(tle, "bstar", state(7,:)); 
endfunction

// ------------------------------------------------
// Propagate initial state and compute measurements
// ------------------------------------------------
function [pos]=compute_pos(cjd0, state0, cjd)
  tle = state2tle(cjd0, state0); 
  [pos, vel] = CL_tle_genEphem(tle, cjd); 
endfunction

// ------------------------------------------------
// TLE determination - high level interface for clarity 
// cjd: dates/times (TREF time scale) (1xN)
// pvref: reference position/velocity vectors (ECI) (6xN)
// ------------------------------------------------
function [tle]=tle_fit(cjd, pvref)
  // Definition of initial state (first guess: osculating elements, bstar = 0)
  cjd0 = cjd(1); 
  equin0 = CL_oe_convert("pv", "equin", pvref(:,1)); 
  state0 = [equin0; 0];  
  
  // Fitting parameters: default values except for: 
  // - svScale: state vector scale factors
  // - verbose: verbose mode 
  params = struct(); 
  params.svScale = [equin0(1); 1 * ones(5,1); 0.1]; 
  params.verbose = 1; 

  // Adjustment 
  state0 = CL_adjustState(cjd0, state0, cjd, pvref(1:3,:), compute_pos, params); 
   
  // Initial state to TLE
  tle = state2tle(cjd0, state0);
endfunction 

// ================================================
// MAIN
// ================================================
// Reference trajectory
[cjd, pvref] = build_trajectory(); 
 
// TLE determination 
tle = tle_fit(cjd, pvref); 
disp(tle); 

// Plot results 
pos = CL_tle_genEphem(tle, cjd); 
plot_gaps_qsw(cjd, pos, pvref);

Of course the quality of the fitting depends on the model used to compute the reference trajectory. In this case the model is rather consistent with SGP4, hence the good results.

Report an issue
<< Comparison of orbit propagation models	Cookbook	Handling time scales >>