Basic orbit determination, application to TLE estimation
The objective here is to solve the following problem:
Find X0 such that:
F(X0, t) = M(t) for t = t1, t2, ... tn
Where:
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 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 | ![]() | ![]() |
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:
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.