Batch Orbit Determination

4.3.13 Batch Orbit Determination

Batch orbit determination estimates a spacecraft state, and optionally additional dynamic or measurement parameters, by processing an entire set of observations over a finite observation arc together.

Instead of updating the state one measurement at a time, batch OD seeks the single reference-epoch state that best explains all available measurements simultaneously.

The central idea is

\[ \boxed{ \text{Reference-Epoch State Guess} + \text{Complete Observation Arc} \rightarrow \text{Best-Fit State and Parameters} } \]

A batch estimator therefore combines:

• orbit propagation,
• measurement prediction,
• STM-based sensitivity mapping,
• residual stacking,
• weighting,
• differential correction,
• iteration,
• convergence testing,
• covariance calculation,
• post-fit residual analysis.

Batch OD is widely used for:

• precise orbit reconstruction,
• post-processing,
• trajectory refinement,
• calibration,
• maneuver estimation,
• parameter estimation,
• generating accurate initial conditions for sequential filters.

4.3.13.1 Basic Batch Estimation Problem

Suppose measurements are available at times

\[ t_1,t_2,\ldots,t_N. \]

Let the unknown reference state be

\[ \boxed{ \mathbf x_0 = \mathbf x(t_0) } \]

at reference epoch t_0.

Each measurement is modeled as

\[ \boxed{ \mathbf z_k = h_k[ \phi( t_k,t_0,\mathbf x_0 )] + \mathbf v_k } \]

where

• \phi(\cdot) propagates the state from t_0 to t_k,
• h_k(\cdot) predicts the measurement,
• \mathbf v_k is measurement noise.

The objective is to estimate the value of \mathbf x_0 that makes the predicted measurements best match the complete measurement set.

4.3.13.2 Why Batch OD Uses a Reference Epoch

The measurements occur at many different times, but a single common state is estimated.

For each observation time,

\[ \boxed{ \mathbf x_k = \phi( t_k,t_0,\mathbf x_0 ) } \]

Therefore all measurements are dynamically linked to the same reference state.

This is one of the defining features of orbit determination.

It ensures that the estimated trajectory is consistent with the equations of motion instead of fitting an unrelated position independently at every time.

4.3.13.3 Initial Reference-State Guess

The batch process begins with an initial guess

\[ \boxed{ \hat{\mathbf x}_0^{(0)} } \]

where the superscript indicates the iteration number.

This initial guess may come from

• an initial orbit determination solution,
• GNSS,
• previous OD,
• propagated ephemeris,
• another navigation source.

The closer the initial guess is to the true state, the more reliable the linearized iterative solution generally becomes.

4.3.13.4 Propagation to Measurement Epochs

At iteration j, propagate

\[ \hat{\mathbf x}_0^{(j)} \]

to every observation epoch:

\[ \boxed{ \hat{\mathbf x}_k^{(j)} = \phi( t_k,t_0, \hat{\mathbf x}_0^{(j)} ) } \]

This propagation uses the selected force model.

For precise OD, the STM and any required parameter sensitivities are usually propagated at the same time.

4.3.13.5 Measurement Prediction

At every observation time, predict the measurement:

\[ \boxed{ \hat{\mathbf z}_k^{(j)} = h_k( \hat{\mathbf x}_k^{(j)} ) } \]

The predicted quantity must use the same physical definition as the actual measurement, including, when required,

• correct coordinate frame,
• sensor geometry,
• station location,
• timing,
• bias corrections,
• light-time,
• calibration.

4.3.13.6 Pre-Fit Residual

The difference between observed and predicted measurement is the pre-fit residual:

\[ \boxed{ \mathbf r_k^{(j)} = \mathbf z_k - \hat{\mathbf z}_k^{(j)} } \]

If the reference-state guess were exact and the models perfect, the residual would contain only measurement noise.

In practice it also contains the effects of

• state error,
• parameter error,
• force-model mismatch,
• sensor bias,
• timing error.

4.3.13.7 Linearization of the Batch Problem

The nonlinear measurement equation is linearized around the current reference-state estimate.

Let

\[ \delta\mathbf x_0 = \mathbf x_0 - \hat{\mathbf x}_0. \]

Then approximately,

\[ \boxed{ \mathbf r_k \approx \mathbf G_k \delta\mathbf x_0 + \mathbf v_k } \]

where

\[ \mathbf G_k \]

is the measurement sensitivity with respect to the reference-epoch state.

4.3.13.8 Measurement Jacobian at the Observation Epoch

At measurement time t_k,

\[ \boxed{ \mathbf H_k = \frac{\partial h_k} {\partial\mathbf x_k} } \]

describes how the measurement changes with the state at that same epoch.

However, batch OD usually estimates \mathbf x_0, not \mathbf x_k.

Therefore the STM is required.

4.3.13.9 STM-Based Design Matrix

The state perturbation at measurement epoch t_k is

\[ \boxed{ \delta\mathbf x_k = \boldsymbol\Phi( t_k,t_0 ) \delta\mathbf x_0 } \]

Therefore,

\[ \delta\mathbf z_k \approx \mathbf H_k \boldsymbol\Phi( t_k,t_0 ) \delta\mathbf x_0. \]

Define

\[ \boxed{ \mathbf G_k = \mathbf H_k \boldsymbol\Phi( t_k,t_0 ) } \]

Then

\[ \boxed{ \mathbf r_k \approx \mathbf G_k \delta\mathbf x_0 + \mathbf v_k. } \]

This is the central STM-based relationship used in batch orbit determination.

4.3.13.10 Physical Meaning of `\mathbf G_k`

The matrix

\[ \mathbf G_k \]

answers:

If I perturb the reference-epoch state slightly, how will the measurement at time t_k change?

It combines two sensitivities:

\[ \boxed{ \text{Reference-State Change} \xrightarrow{\boldsymbol\Phi} \text{State Change at }t_k \xrightarrow{\mathbf H_k} \text{Measurement Change} } \]

Thus,

\[ \boxed{ \mathbf G_k = \mathbf H_k\boldsymbol\Phi_k } \]

connects orbital dynamics directly to measurement information.

4.3.13.11 Residual Stacking

Instead of treating the residuals individually, batch OD stacks them into one large vector:

\[ \boxed{ \mathbf r = \begin{bmatrix} \mathbf r_1\\ \mathbf r_2\\ \vdots\\ \mathbf r_N \end{bmatrix} } \]

Similarly, stack the design matrices:

\[ \boxed{ \mathbf G = \begin{bmatrix} \mathbf G_1\\ \mathbf G_2\\ \vdots\\ \mathbf G_N \end{bmatrix} } \]

Then the entire linearized batch problem becomes

\[ \boxed{ \mathbf r \approx \mathbf G \delta\mathbf x_0 + \mathbf v. } \]

4.3.13.12 Dimensions of the Batch System

Suppose the unknown state contains n elements and there are m total scalar measurements.

Then

\[ \mathbf G \in \mathbb R^{m\times n} \]

and

\[ \delta\mathbf x_0 \in \mathbb R^n. \]

Usually,

\[ m>n, \]

so the system is overdetermined.

The estimator therefore seeks the correction that best fits all measurements in a least-squares sense.

4.3.13.13 Unweighted Least Squares

The simplest objective is

\[ \boxed{ J = \mathbf r^T\mathbf r } \]

or equivalently

\[ \boxed{ J = \sum_{k=1}^{N} \mathbf r_k^T \mathbf r_k. } \]

The estimate minimizes the sum of squared residuals.

The corresponding normal equations are

\[ \boxed{ \mathbf G^T \mathbf G \, \delta\mathbf x_0 = \mathbf G^T \mathbf r. } \]

4.3.13.14 Unweighted Least-Squares Solution

If the required inverse exists,

\[ \boxed{ \delta\mathbf x_0 = ( \mathbf G^T\mathbf G )^{-1} \mathbf G^T \mathbf r } \]

This is the classical least-squares correction.

However, orbit-determination measurements rarely all have identical uncertainty.

Therefore weighted least squares is usually more appropriate.

4.3.13.15 Weighted Least Squares

If the stacked measurement-noise covariance is

\[ \mathbf R, \]

define the weight matrix

\[ \boxed{ \mathbf W = \mathbf R^{-1}. } \]

Then the weighted cost function is

\[ \boxed{ J = \mathbf r^T \mathbf W \mathbf r. } \]

Measurements with lower uncertainty receive greater statistical influence.

4.3.13.16 Weighted Normal Equations

Minimizing the weighted cost gives

\[ \boxed{ \mathbf G^T \mathbf W \mathbf G \, \delta\mathbf x_0 = \mathbf G^T \mathbf W \mathbf r. } \]

The formal solution is

\[ \boxed{ \delta\mathbf x_0 = ( \mathbf G^T \mathbf W \mathbf G )^{-1} \mathbf G^T \mathbf W \mathbf r. } \]

This is the standard weighted batch least-squares correction.

4.3.13.17 Normal Matrix

Define the normal matrix

\[ \boxed{ \mathbf N = \mathbf G^T \mathbf W \mathbf G. } \]

Then

\[ \boxed{ \mathbf N \delta\mathbf x_0 = \mathbf b } \]

with

\[ \boxed{ \mathbf b = \mathbf G^T \mathbf W \mathbf r. } \]

The matrix \mathbf N contains accumulated information from the complete observation arc.

4.3.13.18 Information Interpretation

Because

\[ \mathbf N = \mathbf G^T \mathbf R^{-1} \mathbf G, \]

it has the form of an information matrix.

Thus,

\[ \boxed{ \text{More Informative Measurements} \rightarrow \text{Larger Information} \rightarrow \text{Smaller State Uncertainty} } \]

in an appropriate matrix sense.

Weakly observed directions appear as poorly conditioned or nearly singular directions of \mathbf N.

4.3.13.19 Why Weighting Matters

Suppose one range sensor has

\[ \sigma_\rho=1\ \text{m} \]

while another has

\[ \sigma_\rho=100\ \text{m}. \]

Treating both equally would ignore a large difference in measurement quality.

Using

\[ W_i=\frac{1}{\sigma_i^2} \]

naturally gives greater influence to the more precise measurement.

4.3.13.20 Measurement Covariance in Batch OD

If measurements are independent,

\[ \boxed{ \mathbf R = \operatorname{blockdiag} ( \mathbf R_1, \mathbf R_2, \ldots, \mathbf R_N ) } \]

If measurement errors are correlated,

\[ \mathbf R \]

must contain the appropriate off-diagonal terms.

Ignoring correlation can cause the batch solution to become overconfident.

4.3.13.21 A Priori State Information

Batch OD may also include prior knowledge of the reference state.

Let

\[ \boxed{ \mathbf x_0 \sim \mathcal N( \bar{\mathbf x}_0, \mathbf P_0 ) } \]

Then the batch objective can include both measurement residuals and a prior penalty.

4.3.13.22 Batch Cost Function with Prior

A common objective is

\[ \boxed{ J = ( \mathbf x_0-\bar{\mathbf x}_0 )^T \mathbf P_0^{-1} ( \mathbf x_0-\bar{\mathbf x}_0 ) + \mathbf r^T \mathbf R^{-1} \mathbf r. } \]

Thus the solution balances

• prior orbit knowledge,
• new measurement information.

4.3.13.23 Prior-Weighted Normal Equations

If the correction is defined around current estimate \hat{\mathbf x}_0, prior information contributes an additional term.

Conceptually,

\[ \boxed{ \mathbf N_{\rm total} = \mathbf P_0^{-1} + \mathbf G^T \mathbf R^{-1} \mathbf G. } \]

The corresponding posterior covariance is approximately

\[ \boxed{ \mathbf P_{\rm post} = \mathbf N_{\rm total}^{-1}. } \]

4.3.13.24 Differential Correction

Because spacecraft dynamics and measurement functions are nonlinear, one least-squares solution is generally not sufficient.

The procedure is iterative.

At iteration j:

\[ \boxed{ \hat{\mathbf x}_0^{(j+1)} = \hat{\mathbf x}_0^{(j)} + \delta\mathbf x_0^{(j)}. } \]

This process is called differential correction.

4.3.13.25 Why Differential Correction Is Required

The equation

\[ \mathbf r \approx \mathbf G\delta\mathbf x_0 \]

is only a local linear approximation.

If the initial state error is large, the first linearization may be inaccurate.

After correcting the state, the trajectory and measurement sensitivities must therefore be recomputed.

Thus,

\[ \boxed{ \text{Propagate} \rightarrow \text{Linearize} \rightarrow \text{Solve} \rightarrow \text{Correct} \rightarrow \text{Repeat} } \]

4.3.13.26 Iterative Batch OD Cycle

A typical iteration follows:

\[ \boxed{ \begin{array}{c} \text{Current Reference State}\\ \hat{\mathbf x}_0^{(j)} \\ \downarrow\\ \text{Propagate State + STM}\\ \downarrow\\ \text{Predict Measurements}\\ \downarrow\\ \text{Compute Residuals}\\ \downarrow\\ \text{Form }\mathbf G\\ \downarrow\\ \text{Solve Weighted Least Squares}\\ \downarrow\\ \delta\mathbf x_0^{(j)} \\ \downarrow\\ \hat{\mathbf x}_0^{(j+1)} = \hat{\mathbf x}_0^{(j)} + \delta\mathbf x_0^{(j)} \\ \downarrow\\ \text{Repeat Until Converged} \end{array} } \]

4.3.13.27 Gauss-Newton Interpretation

Differential correction is closely related to the Gauss-Newton method for nonlinear least squares.

Given

\[ J(\mathbf x) = \mathbf r(\mathbf x)^T \mathbf W \mathbf r(\mathbf x), \]

the local quadratic approximation leads to

\[ \boxed{ \delta\mathbf x = ( \mathbf G^T \mathbf W \mathbf G )^{-1} \mathbf G^T \mathbf W \mathbf r. } \]

Thus batch OD can be understood as a nonlinear optimization problem.

4.3.13.28 Sign Convention in Differential Correction

Care must be taken with residual definition.

If

\[ \mathbf r = \mathbf z_{\rm obs} - \mathbf z_{\rm calc}, \]

the correction sign differs from the convention obtained if

\[ \mathbf r = \mathbf z_{\rm calc} - \mathbf z_{\rm obs}. \]

Therefore the residual convention, Jacobian sign, and update equation must be internally consistent.

4.3.13.29 Parameter Augmentation in Batch OD

The unknown vector can include both the reference state and model parameters:

\[ \boxed{ \mathbf y = \begin{bmatrix} \mathbf x_0\\ \boldsymbol\theta \end{bmatrix}. } \]

Then

\[ \boxed{ \delta\mathbf y = \begin{bmatrix} \delta\mathbf x_0\\ \delta\boldsymbol\theta \end{bmatrix} } \]

is solved simultaneously.

4.3.13.30 Augmented Batch Design Matrix

For parameters \boldsymbol\theta, the measurement sensitivity becomes

\[ \boxed{ \mathbf G_k^{\rm aug} = \begin{bmatrix} \frac{\partial\mathbf z_k} {\partial\mathbf x_0} & \frac{\partial\mathbf z_k} {\partial\boldsymbol\theta} \end{bmatrix}. } \]

The state term is

\[ \boxed{ \frac{\partial\mathbf z_k} {\partial\mathbf x_0} = \mathbf H_k\boldsymbol\Phi_k. } \]

Parameter sensitivities may contain both direct measurement dependence and dynamic dependence.

4.3.13.31 Dynamic Parameter Sensitivity

If parameter \theta affects the dynamics,

\[ \boxed{ \frac{\partial\mathbf z_k} {\partial\theta} = \mathbf H_k \mathbf S_{\theta,k} + \mathbf H_{\theta,k}. } \]

Here,

\[ \mathbf S_{\theta,k} = \frac{\partial\mathbf x_k}{\partial\theta}. \]

This is how parameters such as drag coefficients or force-model terms enter the batch design matrix.

4.3.13.32 Batch Estimation of Biases

Suppose a range bias b_\rho is estimated.

The model is

\[ z_\rho = h_\rho(\mathbf x) + b_\rho + v. \]

Then

\[ \boxed{ \frac{\partial z_\rho} {\partial b_\rho} = 1. } \]

A column of ones may therefore appear in the design matrix for that measurement bias.

4.3.13.33 Batch Estimation of Ballistic or Dynamic Parameters

A drag-related parameter may influence the trajectory over the full arc.

Its sensitivity can accumulate with time.

Therefore long observation arcs may improve parameter estimation.

However, long arcs also increase sensitivity to unmodeled dynamics.

The same arc-length tradeoff discussed previously remains important.

4.3.13.34 Maneuver Estimation in Batch OD

If a maneuver occurs during the observation arc, its correction may be estimated.

For an impulsive event,

\[ \boxed{ \mathbf v^+ = \mathbf v^- + \Delta\mathbf v. } \]

An uncertain

\[ \Delta\mathbf v \]

may be included as an estimated parameter.

Measurements before and after the maneuver then constrain the maneuver correction.

4.3.13.35 State Discontinuity at a Maneuver

A maneuver creates a state discontinuity in velocity.

Therefore the STM and sensitivity propagation must account for the event.

Conceptually,

\[ \boxed{ \delta\mathbf x^+ = \mathbf M \delta\mathbf x^- + \mathbf B_{\Delta v} \delta\Delta\mathbf v } \]

where \mathbf M is the event transition mapping.

The details depend on the maneuver model.

4.3.13.36 Multiple-Shooting or Multi-Arc Batch Estimation

For long arcs or maneuver-rich trajectories, a single reference state may become inconvenient.

The arc can be divided into segments with separate local states.

Continuity constraints then connect the segments.

This leads to

• multiple-shooting methods,
• multi-arc estimation.

These approaches can improve numerical robustness for complex trajectories.

4.3.13.37 Multi-Arc Orbit Determination

Suppose separate arcs have states

\[ \mathbf x_0^{(1)}, \mathbf x_0^{(2)},\ldots. \]

Each arc may have its own measurements while selected global parameters are shared.

This allows simultaneous estimation of

• local states,
• common physical parameters,
• measurement biases.

Multi-arc estimation is useful when continuity is imperfect or model behavior changes with time.

4.3.13.38 Residual Weighting by Measurement Type

Different measurement types may have different uncertainties.

For example,

\[ \sigma_\rho,\quad \sigma_{\dot\rho},\quad \sigma_A,\quad \sigma_E \]

may differ substantially.

Each residual should therefore be weighted according to its covariance.

Otherwise one measurement type may dominate the solution simply because of its numerical units.

4.3.13.39 Unit Scaling in Batch OD

Suppose the measurement vector contains

• meters,
• radians,
• hertz.

Directly minimizing unweighted numerical residuals would be meaningless because the units differ.

Weighting with \mathbf R^{-1} converts each residual into a statistically appropriate scale.

Thus weighted least squares also provides a natural solution to mixed measurement units.

4.3.13.40 Whitening the Batch Problem

If

\[ \mathbf R = \mathbf L\mathbf L^T, \]

define

\[ \boxed{ \tilde{\mathbf r} = \mathbf L^{-1}\mathbf r } \]

and

\[ \boxed{ \tilde{\mathbf G} = \mathbf L^{-1}\mathbf G. } \]

Then the weighted problem becomes an ordinary least-squares problem:

\[ \boxed{ \min \| \tilde{\mathbf r} - \tilde{\mathbf G} \delta\mathbf x \|^2. } \]

This is often numerically useful.

4.3.13.41 Solving Normal Equations Directly

Although

\[ \delta\mathbf x = ( \mathbf G^T\mathbf W\mathbf G )^{-1} \mathbf G^T\mathbf W\mathbf r \]

is useful mathematically, explicitly computing the matrix inverse is generally discouraged.

Instead solve

\[ \boxed{ \mathbf N\delta\mathbf x = \mathbf b } \]

using an appropriate linear solver.

4.3.13.42 Why Explicit Matrix Inversion Is Usually Avoided

Explicit inversion can

• increase numerical error,
• waste computation,
• hide conditioning problems.

Prefer

• Cholesky decomposition,
• QR decomposition,
• SVD,
• other stable factorization methods.

Thus,

\[ \boxed{ \text{Solve Linear System} \quad\text{rather than}\quad \text{Explicitly Form } \mathbf N^{-1}. } \]

4.3.13.43 QR Solution

For the whitened problem,

\[ \tilde{\mathbf G} \delta\mathbf x \approx \tilde{\mathbf r}, \]

perform

\[ \boxed{ \tilde{\mathbf G} = \mathbf Q\mathbf R. } \]

Then solve the triangular system

\[ \boxed{ \mathbf R \delta\mathbf x = \mathbf Q^T \tilde{\mathbf r}. } \]

QR avoids explicitly forming \mathbf G^T\mathbf G, which can improve numerical conditioning.

4.3.13.44 SVD Solution

For a poorly conditioned design matrix,

\[ \boxed{ \tilde{\mathbf G} = \mathbf U \boldsymbol\Sigma \mathbf V^T. } \]

Then the least-squares correction can be computed using the singular values.

SVD reveals directly

• rank deficiency,
• weak estimation directions,
• parameter correlations.

It is especially useful for diagnostic analysis.

4.3.13.45 Condition Number

The condition number provides a measure of numerical sensitivity:

\[ \boxed{ \kappa(\mathbf G) = \frac{\sigma_{\max}} {\sigma_{\min}} } \]

for a full-rank design matrix.

A very large condition number indicates that small changes in measurements may produce large changes in some estimated state combinations.

4.3.13.46 Rank Deficiency

If

\[ \operatorname{rank}(\mathbf G)<n, \]

then some state or parameter combination cannot be uniquely determined from the available measurements.

In that case,

\[ \boxed{ \mathbf G^T \mathbf W \mathbf G } \]

is singular or nearly singular.

This is an observability or identifiability problem, not merely a numerical inconvenience.

4.3.13.47 Weak Estimation Directions

Even if the design matrix is technically full rank, one singular value may be very small.

Then the corresponding direction in parameter space is weakly constrained.

This may result in

• large covariance,
• strong parameter correlation,
• unstable corrections.

Such directions should be identified rather than hidden by aggressive numerical regularization.

4.3.13.48 State Scaling

Batch variables may have very different numerical magnitudes.

For example,

\[ r\sim10^6\,{\rm m} \]

and a calibration parameter might be

\[ 10^{-5}. \]

Define scaled variables:

\[ \boxed{ \delta\mathbf y = \mathbf S \delta\tilde{\mathbf y}. } \]

Proper scaling can substantially improve conditioning.

4.3.13.49 Convergence Criteria

A batch iteration may be stopped when one or more conditions are satisfied.

For example,

\[ \boxed{ \| \delta\mathbf x_0 \| < \epsilon_x } \]

or

\[ \boxed{ |J^{(j+1)}-J^{(j)}| < \epsilon_J. } \]

Other criteria may examine

• residual RMS,
• normalized correction,
• parameter change.

4.3.13.50 Component-Wise Convergence

Because position and velocity have different units, using one raw Euclidean norm can be misleading.

One may use normalized corrections:

\[ \boxed{ \eta_i = \frac{ |\delta x_i| }{ s_i } } \]

where s_i is a characteristic scale or requirement.

Then require all relevant components to fall below prescribed tolerances.

4.3.13.51 Residual RMS Convergence

Another metric is residual root-mean-square:

\[ \boxed{ {\rm RMS} = \sqrt{ \frac{1}{m} \mathbf r^T\mathbf r } } \]

or weighted RMS,

\[ \boxed{ {\rm WRMS} = \sqrt{ \frac{1}{m} \mathbf r^T \mathbf W \mathbf r }. } \]

Convergence may be indicated when WRMS changes negligibly between iterations.

4.3.13.52 Cost Function Reduction

For a well-behaved differential-correction process,

\[ J^{(j+1)} < J^{(j)} \]

is generally expected.

If the cost repeatedly increases, possible causes include

• poor linearization,
• incorrect Jacobian,
• bad initial state,
• wrong sign convention,
• model inconsistency.

4.3.13.53 Damping the Correction

When the full Gauss-Newton correction is too aggressive, apply

\[ \boxed{ \mathbf x^{(j+1)} = \mathbf x^{(j)} + \lambda \delta\mathbf x^{(j)} } \]

where

\[ 0<\lambda\le1. \]

This can improve robustness when the initial guess is not close enough for the local linearization.

4.3.13.54 Levenberg-Marquardt Idea

A damped least-squares approach modifies the normal matrix:

\[ \boxed{ ( \mathbf G^T\mathbf W\mathbf G + \lambda\mathbf I ) \delta\mathbf x = \mathbf G^T\mathbf W\mathbf r. } \]

This blends Gauss-Newton behavior with additional regularization.

It can improve convergence for nonlinear or ill-conditioned problems.

4.3.13.55 Convergence Is Not Correctness

A differential-correction algorithm may converge numerically to an incorrect solution.

Possible causes include

• wrong measurement model,
• wrong force model,
• biased measurements,
• local minimum,
• poor observability.

Thus,

\[ \boxed{ \text{Converged Batch Solution} \neq \text{Verified Physical Solution} } \]

Validation remains essential.

4.3.13.56 Post-Fit Residuals

After convergence, propagate the final estimated state and recompute the measurement residuals:

\[ \boxed{ \mathbf r_k^{\rm post} = \mathbf z_k - h_k( \hat{\mathbf x}_k ). } \]

These are the post-fit residuals.

They are among the most important diagnostics of a batch OD solution.

4.3.13.57 Desired Properties of Post-Fit Residuals

Ideally, post-fit residuals should be

• approximately zero mean,
• consistent with expected measurement noise,
• free of obvious trends,
• free of unmodeled periodic structure,
• approximately white when the measurement model assumes white noise.

This does not mean every residual must be small.

Individual noise realizations can naturally be several standard deviations away from zero.

4.3.13.58 Residual Bias

If

\[ \boxed{ E[\mathbf r^{\rm post}] \neq0, } \]

the solution may contain systematic mismatch.

Possible causes include

• measurement bias,
• incorrect station geometry,
• force-model error,
• timing error,
• calibration error.

A biased residual sequence should be investigated physically.

4.3.13.59 Residual Trends

A residual that steadily increases or decreases with time may indicate

• incorrect velocity,
• unmodeled acceleration,
• incorrect dynamic parameter,
• maneuver mismatch.

Thus residual shape contains diagnostic information beyond residual RMS.

4.3.13.60 Periodic Residual Structure

Periodic signatures may arise from

• orbital force-model error,
• station geometry,
• attitude-dependent calibration,
• periodic environmental disturbances.

Frequency-domain or time-series analysis can help identify these patterns.

4.3.13.61 Residual Whitening

For correlated measurement covariance,

\[ \mathbf R = \mathbf L\mathbf L^T, \]

the whitened residual is

\[ \boxed{ \tilde{\mathbf r} = \mathbf L^{-1} \mathbf r. } \]

A correctly modeled solution should produce whitened residuals with approximately unit covariance.

4.3.13.62 Normalized Residuals

For scalar measurement variance \sigma_i^2,

\[ \boxed{ r_{n,i} = \frac{r_i}{\sigma_i}. } \]

Normalized residuals allow different measurement types to be compared on a common statistical scale.

4.3.13.63 Chi-Square Cost

The weighted residual cost

\[ \boxed{ J = \mathbf r^T \mathbf R^{-1} \mathbf r } \]

has a chi-square interpretation under Gaussian assumptions.

This enables statistical consistency checks on the overall residual set.

Care is required to use the appropriate degrees of freedom.

4.3.13.64 Degrees of Freedom

If there are m independent scalar measurements and n estimated independent parameters, a common approximate number of residual degrees of freedom is

\[ \boxed{ \nu = m-n. } \]

Then the weighted post-fit cost can be compared with expected chi-square behavior.

This relation should be interpreted carefully when priors, constraints, correlations, or rank deficiencies are present.

4.3.13.65 Post-Fit Covariance

Near the converged solution, the estimated parameter covariance is approximately

\[ \boxed{ \mathbf P = ( \mathbf G^T \mathbf R^{-1} \mathbf G )^{-1} } \]

for the measurement-only full-rank linearized case.

With prior information,

\[ \boxed{ \mathbf P = ( \mathbf P_0^{-1} + \mathbf G^T \mathbf R^{-1} \mathbf G )^{-1}. } \]

4.3.13.66 Interpretation of Formal Covariance

The covariance is based on the assumed

• measurement statistics,
• force model,
• measurement model,
• local linearization.

Therefore it is sometimes called a formal covariance.

A mathematically small formal covariance does not guarantee small true error if the model assumptions are incorrect.

4.3.13.67 Correlation Matrix

Parameter correlations may be summarized using

\[ \boxed{ \rho_{ij} = \frac{P_{ij}} {\sqrt{P_{ii}P_{jj}}}. } \]

Large

\[ |\rho_{ij}| \]

indicates strong coupling.

Strong correlation may reveal that two estimated quantities are difficult to separate using the available observation geometry.

4.3.13.68 Error Ellipsoid from Batch Covariance

The position block

\[ \mathbf P_{rr} \]

can be used to form an uncertainty ellipsoid.

Eigenvalues and eigenvectors identify principal uncertainty magnitudes and directions.

Therefore batch OD produces both

\[ \boxed{ \text{Best-Fit State} } \]

and

\[ \boxed{ \text{Best-Fit Uncertainty Geometry}. } \]

4.3.13.69 Pre-Fit versus Post-Fit Residuals

Pre-fit residuals use the state before differential correction.

Post-fit residuals use the final fitted state.

Typically,

\[ \boxed{ \|\mathbf r^{\rm post}\| < \|\mathbf r^{\rm pre}\| } \]

in an overall least-squares sense.

However, every individual residual does not have to decrease.

4.3.13.70 Measurement Editing Before Batch Solution

Measurements may first be filtered using

• quality flags,
• physical bounds,
• known sensor-health information.

Clearly invalid observations should not be allowed to dominate the least-squares fit.

However, aggressive editing can bias the resulting solution.

4.3.13.71 Iterative Outlier Rejection

One approach is:

1. solve initial batch,
2. examine normalized residuals,
3. identify statistically inconsistent observations,
4. investigate or reject them,
5. rerun the batch solution.

This process should use consistent statistical criteria.

4.3.13.72 Robust Least Squares

Instead of completely rejecting large residuals, robust estimation can reduce their influence.

A robust cost may replace

\[ r^2 \]

with a function

\[ \rho(r) \]

that grows more slowly for large residuals.

Examples include

• Huber loss,
• Tukey-type weighting.

This is useful when occasional outliers are expected.

4.3.13.73 Effective Measurement Weights

In robust estimation, the effective weight of a measurement becomes residual dependent.

Large residuals can be downweighted rather than fully discarded.

However, robust estimation should not be used to hide systematic model errors.

4.3.13.74 Observation Arc Selection

The observation arc strongly affects batch OD.

A short arc may provide insufficient dynamic information.

A long arc can improve observability but becomes more sensitive to force-model error.

Thus,

\[ \boxed{ \text{Arc Length} = \text{Information Gain} \leftrightarrow \text{Model Fidelity} } \]

4.3.13.75 Reference Epoch Placement

The reference epoch can be placed

• at the beginning of the arc,
• midpoint,
• end.

Placing it near the middle can sometimes reduce state correlation and numerical propagation distances in batch problems.

The best choice depends on the objective.

4.3.13.76 Epoch Transformation of the Final Solution

If the estimated state is required at another epoch t_a,

\[ \boxed{ \hat{\mathbf x}(t_a) = \phi( t_a,t_0,\hat{\mathbf x}_0 ) } \]

and covariance may be mapped as

\[ \boxed{ \mathbf P(t_a) = \boldsymbol\Phi(t_a,t_0) \mathbf P_0 \boldsymbol\Phi^T(t_a,t_0) } \]

with appropriate process-noise terms when applicable.

4.3.13.77 Measurement Geometry Across the Arc

Measurements obtained from a single unchanging geometry may produce a poorly conditioned design matrix.

Changing geometry can improve the rank and singular-value distribution of

\[ \mathbf G. \]

Therefore diverse observation geometry is often more valuable than simply increasing measurement count.

4.3.13.78 Repeated Similar Measurements

If many observations have nearly identical sensitivity rows,

\[ \mathbf G_i \approx \mathbf G_j, \]

they provide partly redundant information.

Increasing measurement count does not always create proportionally more independent information.

4.3.13.79 Batch OD and Observability

The batch design matrix provides a practical observability view.

If

\[ \mathbf G \]

has full column rank, the locally linearized unknowns are identifiable from the measurement set.

If not, some state combination lies in its null space.

This connects batch OD directly to geometry and estimability analysis.

4.3.13.80 Batch Information Matrix Eigenstructure

The eigenvectors of

\[ \boxed{ \mathbf N = \mathbf G^T \mathbf W \mathbf G } \]

identify combinations of parameters.

Large eigenvalues correspond to strongly constrained directions.

Small eigenvalues correspond to weakly constrained directions.

This is useful for diagnosing parameter-estimation quality.

4.3.13.81 Measurement-Type Combination

Batch OD naturally combines heterogeneous measurements.

For example,

\[ \boxed{ \begin{bmatrix} \text{Range Residuals}\\ \text{Doppler Residuals}\\ \text{Angle Residuals}\\ \text{GNSS Residuals} \end{bmatrix} } \]

may all appear in one stacked residual vector.

Their covariance and Jacobians determine their relative information contribution.

4.3.13.82 Multi-Sensor Bias Estimation

Separate bias parameters may be added for different sensors.

For example,

\[ \boxed{ \mathbf y = \begin{bmatrix} \mathbf x_0\\ b_1\\ b_2\\ b_3 \end{bmatrix}. } \]

But each added bias must be observable.

Otherwise the design matrix can become poorly conditioned.

4.3.13.83 Correlated Measurement Sources

If measurements from different sensors share a common error source, their covariance should contain cross terms.

For example,

• common clock error,
• common atmospheric model error,
• shared preprocessing.

Ignoring those correlations can make the solution too confident.

4.3.13.84 Parameter Priors

Weak parameters may be stabilized using prior information:

\[ \boxed{ \theta \sim \mathcal N( \bar\theta, \sigma_\theta^2 ). } \]

This adds information to the normal matrix.

However, a prior-constrained parameter should not be described as strongly observed by the measurements if the measurements themselves contain little information.

4.3.13.85 Constraint Equations

Some parameters may satisfy physical constraints.

A linear constraint can be written

\[ \boxed{ \mathbf C \delta\mathbf x = \mathbf d. } \]

Constrained least squares can incorporate such relationships.

Constraints should represent real physical or mathematical knowledge.

4.3.13.86 Equality Constraints

An exact equality constraint can be introduced through

• Lagrange multipliers,
• state reduction,
• constrained optimization.

This is different from a probabilistic prior, which allows uncertainty.

4.3.13.87 Soft Constraints

A soft constraint treats the relationship as uncertain.

For example,

\[ \boxed{ \mathbf c(\mathbf x) \approx \mathbf d } \]

with assigned covariance.

It can be included as a pseudo-measurement.

This preserves a statistical interpretation.

4.3.13.88 Pseudo-Measurements

A known approximate relationship can be represented as

\[ \boxed{ \mathbf z_p = h_p(\mathbf x) + \mathbf v_p. } \]

The pseudo-measurement is then included in the same weighted least-squares framework.

Its covariance determines how strongly it constrains the solution.

4.3.13.89 Batch OD for Relative Navigation

The same mathematics applies to a general relative-navigation state:

\[ \boxed{ \mathbf x_{\rm rel,0} = \begin{bmatrix} \mathbf r_{\rm rel,0}\\ \mathbf v_{\rm rel,0} \end{bmatrix}. } \]

Measurements over an arc may include

• range,
• range rate,
• LOS,
• bearing angles,
• optical features.

The complete observation set can be used to estimate a best-fit reference relative state.

This is a general application to relative orbit determination, including formation flying and rendezvous/proximity operations.

4.3.13.90 Relative-State Design Matrix

For relative measurement

\[ \mathbf z_k = h_k( \mathbf x_{{\rm rel},k} ), \]

the reference-state sensitivity is

\[ \boxed{ \mathbf G_k = \mathbf H_{{\rm rel},k} \boldsymbol\Phi_{\rm rel}(t_k,t_0). } \]

The same residual stacking and weighted least-squares equations apply.

4.3.13.91 Batch OD for Optical Navigation

For optical measurements,

\[ \mathbf z_k \]

may contain

• pixels,
• angles,
• LOS vectors.

Batch estimation can improve range observability because measurements collected at different geometries are processed together.

Thus a weak instantaneous range direction may become well constrained over the full observation arc.

4.3.13.92 Batch Calibration

Because the complete observation arc is available, batch estimation is especially useful for calibration.

Possible calibration states include

• camera alignment,
• range bias,
• clock offset,
• sensor scale factor,
• station bias.

The accumulated geometry may provide better parameter separation than a short real-time filter interval.

4.3.13.93 Batch Reconstruction after Flight

Post-flight analysis can use all collected observations, detailed force models, and relaxed computational constraints.

Batch OD can therefore produce highly accurate trajectory reconstruction.

This may be valuable for

• performance assessment,
• anomaly investigation,
• scientific data processing,
• system validation.

4.3.13.94 Batch OD versus Sequential OD

Batch OD uses the entire arc together:

\[ \boxed{ \{\mathbf z_1,\ldots,\mathbf z_N\} \rightarrow \hat{\mathbf x}_0. } \]

Sequential OD updates recursively:

\[ \boxed{ \hat{\mathbf x}_{k-1} \rightarrow \hat{\mathbf x}_k. } \]

Both use the same

• dynamics,
• measurement equations,
• uncertainty concepts.

The difference is primarily how information is accumulated and processed.

4.3.13.95 Advantages of Batch OD

Batch estimation provides several advantages:

• uses all observations simultaneously,
• naturally supports post-fit residual analysis,
• useful for high-accuracy reconstruction,
• convenient for parameter estimation,
• can exploit future measurements relative to earlier states,
• useful for calibration,
• provides strong diagnostic capability.

4.3.13.96 Limitations of Batch OD

Limitations include:

• measurements must typically be available before final processing,
• repeated propagation may be computationally expensive,
• poor initial guesses can cause convergence difficulty,
• long arcs increase model sensitivity,
• large parameter sets can create conditioning problems.

Thus batch OD is especially suitable for offline or ground-based estimation, although batch-like optimization can also be used onboard when resources permit.

4.3.13.97 Common Mistakes in Batch Orbit Determination

Common errors include:

1. Forming the design matrix without the STM.
2. Using \mathbf H_k directly when the unknown is the reference-epoch state.
3. Mixing residual sign conventions.
4. Using unweighted least squares for measurements with different uncertainties.
5. Explicitly inverting poorly conditioned normal matrices.
6. Ignoring measurement correlations.
7. Failing to scale states and parameters.
8. Adding unobservable parameters.
9. Using an overly long observation arc with inadequate force modeling.
10. Declaring convergence based only on a small correction.
11. Ignoring post-fit residual structure.
12. Treating formal covariance as guaranteed truth error.
13. Failing to verify Jacobians.
14. Using sensor resolution instead of realistic covariance.
15. Rejecting outliers without statistical justification.
16. Ignoring maneuver discontinuities.
17. Estimating too many correlated biases simultaneously.
18. Forgetting to recompute the STM and measurement Jacobians after each nonlinear iteration.
19. Reporting low residual RMS without checking bias or temporal correlation.
20. Interpreting prior-constrained parameters as strongly identified from data.

4.3.13.98 Practical Batch OD Workflow

A disciplined batch process is

\[ \boxed{ \begin{array}{c} \text{Choose Reference Epoch}\\ \downarrow\\ \text{Define Initial State / Parameter Guess}\\ \downarrow\\ \text{Define Prior Covariance if Available}\\ \downarrow\\ \text{Define Dynamics and Measurement Models}\\ \downarrow\\ \text{Propagate State + STM + Sensitivities}\\ \downarrow\\ \text{Predict All Measurements}\\ \downarrow\\ \text{Compute Pre-Fit Residuals}\\ \downarrow\\ \text{Build }\mathbf G_k=\mathbf H_k\boldsymbol\Phi_k\\ \downarrow\\ \text{Stack Residuals and Design Matrix}\\ \downarrow\\ \text{Apply Measurement Weighting}\\ \downarrow\\ \text{Solve for Differential Correction}\\ \downarrow\\ \text{Update State and Parameters}\\ \downarrow\\ \text{Repeat Until Converged}\\ \downarrow\\ \text{Compute Post-Fit Residuals}\\ \downarrow\\ \text{Compute Covariance and Correlations}\\ \downarrow\\ \text{Check Residual Statistics / Conditioning}\\ \downarrow\\ \text{Validate Against Independent Data} \end{array} } \]

4.3.13.99 Core Batch-OD Equations

Measurement model

\[ \boxed{ \mathbf z_k = h_k[ \phi( t_k,t_0,\mathbf x_0 )] + \mathbf v_k } \]

Pre-fit residual

\[ \boxed{ \mathbf r_k = \mathbf z_k - \hat{\mathbf z}_k } \]

STM mapping

\[ \boxed{ \delta\mathbf x_k = \boldsymbol\Phi_k \delta\mathbf x_0 } \]

Reference-state measurement sensitivity

\[ \boxed{ \mathbf G_k = \mathbf H_k \boldsymbol\Phi_k } \]

Stacked batch model

\[ \boxed{ \mathbf r = \mathbf G \delta\mathbf x_0 + \mathbf v } \]

Weighted cost

\[ \boxed{ J = \mathbf r^T \mathbf R^{-1} \mathbf r } \]

Normal equations

\[ \boxed{ \mathbf G^T \mathbf R^{-1} \mathbf G \, \delta\mathbf x_0 = \mathbf G^T \mathbf R^{-1} \mathbf r } \]

Differential correction

\[ \boxed{ \hat{\mathbf x}_0^{(j+1)} = \hat{\mathbf x}_0^{(j)} + \delta\mathbf x_0^{(j)} } \]

Formal covariance

\[ \boxed{ \mathbf P = ( \mathbf G^T \mathbf R^{-1} \mathbf G )^{-1} } \]

for the full-rank, measurement-only linearized case.

4.3.13.100 Key Takeaways

Batch orbit determination estimates a single reference-epoch state that best explains an entire observation arc.

The measurements are linked to the reference state through the spacecraft dynamics:

\[ \boxed{ \mathbf x_k = \phi( t_k,t_0,\mathbf x_0 ). } \]

The predicted measurements are compared with observations to form residuals:

\[ \boxed{ \mathbf r_k = \mathbf z_k - h_k( \hat{\mathbf x}_k ). } \]

Because the unknown is the reference state, the correct measurement sensitivity is not simply \mathbf H_k, but

\[ \boxed{ \mathbf G_k = \mathbf H_k \boldsymbol\Phi(t_k,t_0). } \]

All measurement residuals are stacked:

\[ \boxed{ \mathbf r = \begin{bmatrix} \mathbf r_1\\ \vdots\\ \mathbf r_N \end{bmatrix} } \]

and all sensitivities are stacked into

\[ \boxed{ \mathbf G. } \]

The weighted correction is obtained from

\[ \boxed{ \mathbf G^T \mathbf R^{-1} \mathbf G \, \delta\mathbf x_0 = \mathbf G^T \mathbf R^{-1} \mathbf r. } \]

Because the orbit and measurement models are nonlinear, this correction is applied iteratively:

\[ \boxed{ \text{Propagate} \rightarrow \text{Predict} \rightarrow \text{Residual} \rightarrow \text{Linearize} \rightarrow \text{Correct} \rightarrow \text{Repeat}. } \]

After convergence, the most important outputs are not only the best-fit state but also

• post-fit residuals,
• covariance,
• parameter correlations,
• conditioning,
• prediction performance.

The central engineering principle is:

\[ \boxed{ \text{A Good Batch OD Solution Must} \quad \text{Fit the Data, Be Numerically Well Conditioned,} \quad \text{and Remain Physically Consistent.} } \]

A very small residual norm alone is not sufficient. The residual statistics, estimated covariance, observability, parameter correlations, and independent validation must all support the final solution.

Next: 4.3.14 Sequential Orbit Determination

Batch orbit determination processes an observation arc together to obtain a best-fit reference state and associated parameters. The next section moves to recursive processing, where the orbit state and covariance are propagated and updated as measurements arrive in real time.