Why are uncertainties of orbital state vectors provided as covariance matrixes?

Question

I have recently become interested in testing a numerical orbital propagator with ISS ephemerides, following @RyanC's suggestion in this answer.

Ephemerides for the ISS seem to be publicly available in a format known as OEM (Orbital Ephemeris Message). The format specifications are defined here. After reading it, this seems like an ideal format to distribute ephemerides!

However, something caught my attention. Uncertainties for the provided state vectors can be optionally provided in the form of covariance matrixes. From this answer, as well as from the definition in Wikipedia, I learnt that covariance matrixes give the set of all possible covariances between a number of variables, with the values in the diagonal of the matrix being the squared standard deviations of each of the variables. In the case of the covariance matrixes in OEM files, these are 6x6 symmetric matrixes, since we have 6 variables: the X, Y and Z components of position and velocities.

I understand that the diagonal values of the matrix therefore provide the standard deviations for each of the components of position and velocity (at a given time, which is in accordance with the fact that the covariance matrixes in an OEM file are always accompanied by a specified epoch).

However, how should the other values of the matrix (i.e., covariances between different variables, some being even covariances between a position component and a velocity component) be interpreted? Since the file specification was decided to include the full covariance matrix, and not just the standard deviation for each variable, I believe these covariances must be providing valuable information. But under what circumstances or for what purpose could such information be used?

Answer 1

But under what circumstances or for what purpose could such information be used?

If you transform the 6x6 covariance matrix to the nominal perifocal (PQW) frame, the information contained therein becomes much clearer. For example, the along-track position uncertainty will inevitably be very much higher than the cross-track and radial uncertainties in the position vector. With regard to predicting possible collisions, you are looking for the possible intersection between two cigar-shaped objects rather than two blob-shaped objects. The cross-correlations will also make a lot more sense in the perifocal frame. For example, radial position uncertainty is typically highly correlated with along-track velocity uncertainty, and vice versa.

Answer 2

The full covariance matrices are useful because you need the whole matrix if you want to apply a change of coordinates (as David Hammen's answer alludes to). For a state vector $x$ covariance matrix $C$ and a linear (e.g. coordinate) transformation $H$ applied to $x$, we have:

$$ \begin{align} x' &= H \cdot x \\ C' &= H \cdot C \cdot H^T \end{align} $$

After the coordinate transformation, you can then extract the standard deviations from the diagonal of $C'$. This would not be possible if you only had the diagonal elements of $C$.

More specifically, since the covariance matrix is a positive definite, full-rank matrix, you can decompose $C$ as $C = V \cdot L \cdot V^{-1}$, where $L$ is a diagonal matrix of with the eigenvalues of $C$ and the columns of $V$ are the eigenvectors of $C$. This allows you to see the covariance matrix $C$ as an uncertainty ellipsoid around a state $x$, of which the principle axes are aligned with the eigenvectors and the size along those axes is given by the eigenvalues:

By Nicoguaro - Own work, CC BY 4.0

This is particularly useful in visualizing and analyzing uncertainties in tracking (full state, not only position!), as David Hammen illustrates. You can for example assess the probability of collision and time of collision.

Answer 3

The full covariance matrix is useful in itself for much more than just drawing ellipsoids, and it can be a lot bigger than 6x6. Covariance is essential to the basic operations of orbit determination, in particular to how data from measurements are fused with each other and with the estimated state. I'm going to wave my hands a lot, and generally avoid (or defer to footnotes) the rigorous mathematics, but I hope I can still give some insight about the motivation.

This is the core subject matter of entire textbooks, like Tapley, Schutz, and Born, Statistical Orbit Determination (2004), which weighs in at a bit over 500 pages. There are much shorter review papers (20ish pages), like Stacey & D'Amico (2021) or Hinks & Psiaki (2012), but they are very dense and hard to follow if you haven't already read one of the 500-page versions.

The ellipsoid representation of a covariance matrix says how uncertain the estimated position and velocity are along each direction. Directions near eigenvectors with large eigenvalues are relatively less certain, and directions near eigenvectors with small eigenvalues are relatively more certain. Descriptions of the shapes characteristic of individual radar and optical measurements can be found in this answer.

Covariance matrices are real and symmetric by definition, so unlike most random matrices, all covariance matrices have inverses<1>. The inverse of a covariance matrix is called an information matrix. It can also be drawn as an ellipsoid, but for it, directions near eigenvectors with large eigenvalues are relatively more certain, and directions near eigenvectors with small eigenvalues are relatively less certain. One reason information matrices are interesting is that they are additive. That is, if you make a bunch of different measurements of your spacecraft, and then you want to combine them all into one best estimate of its position and velocity, you weight them by adding up their information.

In the scalar version of the method of least squares (invented by Gauss in 1801 to determine the orbit of the asteroid Ceres), you weight each observation by the reciprocal of its standard deviation squared --- which is another name for the inverse of its variance. Then the estimated variance of the answer is the reciprocal of the sum of the inverses of the individual variances, just like the total resistance of a group of resistors in parallel. The identical procedure is used for combining observations of vector quantities, whose uncertainties combine in exactly the same way, by summing individual information matrices and taking the inverse (matrix equivalent of scalar reciprocal) of the total, forming a covariance matrix which is smaller (less uncertain) than all of the individual uncertainties you started with.

One good way to think about the size of a covariance matrix is by considering the product of all its eigenvalues, which is a constant times the ellipsoid volume in N dimensions; it is also the same number as the determinant of the matrix<2>. A way to turn that size into a distance metric is given by Förstner and Moonen (1999), but that's more complicated than I want to get into just yet<3>. In a typical least squares solution, the total uncertainty is spread equally over time, so the total volume of the error ellipsoid at a given confidence stays constant, but it is allowed to rotate over time. Some representations may even allow it to slosh volume between axes, becoming more or less spherical over time. In a sequential solution, like a Kalman filter, you can see observations happen by watching the uncertainty volume, which grows smoothly in continuous time (in part since nonzero velocity uncertainties directly cause growing position uncertainties), and then becomes suddenly and discretely smaller every time a new measurement is put into the filter.

The off-diagonal elements of the covariance matrix are to me the most interesting part, because those express the correlations among the dimensions. Each measurement is less sensitive to certain dimensions than to others, but the mapping between the dimensions in which the measurement is made and the dimensions in which it is applied can be complicated<4>, and vary with time. This is an important part of planning observations to support orbit determination, especially if you have a widespread observing network. Your predicted orbit comes with an estimated uncertainty (the covariance matrix), and the optimal choice of which measurements to make from which locations at which times is determined by calculating which of those possible measurements would produce the greatest reduction in that covariance. Equivalently, by finding which measurement would give the greatest increase in information, which depends not only on the measurements, but also on what you already know, as described by the information matrix.

Note well that covariance matrices, like all other estimates of uncertainty, are only theoretical calculations based on rough guesses about how good we usually are at making guesses. They never describe the amount by which a measurement is actually off --- if we knew that, we wouldn't bother drawing an ellipsoid, we'd just move the dot! They describe our best guess at about how wrong we are often likely to be, based (hopefully) on past tests which compared our predictions with reality to give a statistical characterization of our typical accuracy (but, as every investment prospectus is careful to say, past performance is not a guarantee of future performance).

I have a few more things to add (like dimensions higher than 6) later, but this is a good start for now.

Footnotes, mostly mathematical

<1> : Covariance matrices are also positive definite, which together with their symmetry implies not only that all their eigenvalues are nonzero (required for them to have inverses), but also that every eigenvalue they have is individually real and positive, and the same is true of their inverses (whose eigenvalues are the reciprocals of those of the original matrix). When you are computing covariance numerically, it can happen that multiple small errors in your computation combine to produce a matrix which is not invertible, or has negative or imaginary eigenvalues, and thus breaks the nice properties that covariance is supposed to have. There are a variety of ways to defeat this, most of which are based on the idea of taking the square root of the covariance matrix, and operating only on that, so that when it comes time to output the computed covariance, they can square their internal representation to generate the output, hopefully guaranteeing that it will have the right properties. There is more than one way to do this, and indeed more than one way to define square root of a matrix; the typical ones used in orbit determination are the Cholesky decomposition and the Square-Root Information Filter (SRIF), which are among the examples treated in Chapter 5 of Tapley, Schutz, and Born, as part of a hundred pages devoted to the topic.

<2> : If you are computing the determinant of a matrix, keep in mind that the way usually taught in schools (at least, in American ones) is only applicable to 2x2 and 3x3 matrices. It is generalizable (expansion by minors), but for an NxN matrix it requires a sum over N factorial terms, which is generally a bad idea. Since you will probably want the eigenvalues themselves, finding them and multiplying them together is not a bad way to proceed, but there are even better ones, some of which appear in this discussion on Math Stack Exchange.

<3> : For covariances $A$ and $B$, the distance $d$ between them may be calculated as $$d(A,B)=\sqrt{\sum_{i=1}^n\ln^2 \lambda_i(A, B)}$$ the square root of the sum of the squares of the natural logarithms of the generalized eigenvalues of $A$ and $B$. In Python, this is numpy.sqrt( numpy.sum( numpy.square (numpy.log( scipy.linalg.eigvalsh(A, B))))); Matlab, or anything else that has a built-in solver for the generalized eigenvalue problem, has an equivalent one-liner. The generalized eigenvalue problem is, given matrices $A$ and $B$, find all scalars $\lambda$ such that $\det(A - \lambda B) = 0$. The usual eigenvalue problem is the case $B=I$, the identity matrix. When $B$ is invertible, this equals finding the eigenvalues of the matrix ($A$ times the inverse of $B$); or at least it should, but in numerical computation frequently the answers are not quite the same, and one obtains $d(B,B)=10^{-14}$, rather than the exactly zero it ought to be. The proof that this is indeed a distance (obeys the triangle inequality $d(A,C)\le d(A,B)+d(B,C)$, etc.) and is the best such formula, in the sense that it equals the Riemannian metric of the "natural torsion-free $GL(N,\mathbb{R})$-invariant connection" on the space of real, symmetric, positive-definite matrices (the set of all possible covariances) is partly in Förstner W., Moonen B. A Metric for Covariance Matrices. In: Grafarend E.W., Krumm F.W., Schwarze V.S. (eds) Geodesy-The Challenge of the 3rd Millennium. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-662-05296-9_31, and partly by reference to Kobayashi & Nomizu, Foundations of Differential Geometry (1963). This paper makes me very happy. :)

<4> : By "the dimensions in which the measurement is made", I mean things like azimuth and elevation angles, range and range-rate, etc. By "the dimensions in which it is applied" I mean the parameters in which the state vector is expressed, which can be position & velocity, or Keplerian elements, or equinoctial, or anything else you may desire. The mapping between them can be very complicated, but for many kinds of measurement processing, we want to use the tools of linear algebra because they are so convenient, but this generally means we need to linearize the problem in the neighborhood of an approximate solution, compute a change that will move our answer closer to correct, and iterate until we reach some convergence criterion. Therefore, when you look at this relationship written as a matrix, usually in the form $y=Hx$, $y$ and $x$ really represent not the full measurement ($Y$) and state ($X$) vectors, but rather the residual in the current iteration of the approximate difference between the actual measurement $Y_\mathrm{obs}$ and the measurement predicted from the current state estimate $Y_\mathrm{calc}$ (known), and of the difference between the current estimate of the state $X_\mathrm{est}$ and the actual state $X_\mathrm{true}$ (unknown), in which case $H$ is the matrix of partial derivatives of the measurements with respect to the state. That is, for some generally nonlinear model function $G$, which gives instructions for how to predict what the measurements will be if you know the state, such that $Y_\mathrm{calc} = G(X_\mathrm{est})$, then defining $y=Y_\mathrm{obs}-Y_\mathrm{calc}$ and $x=X_\mathrm{true}-X_\mathrm{est}$, we find $Y_\mathrm{obs}=G(X_\mathrm{true})$ may be approximated as $$Y_\mathrm{obs}=G(X_\mathrm{est})+\frac{\partial G}{\partial X}(X_\mathrm{true}-X_\mathrm{est})+\frac{\partial^2 G}{\partial X^2}(X_\mathrm{true}-X_\mathrm{est})^2+\cdots\\ Y_\mathrm{obs}\approx Y_\mathrm{calc}+\frac{\partial G}{\partial X}(X_\mathrm{true}-X_\mathrm{est})\\y \approx H x$$ The operation of estimating the state given the measurements is trying to find the inverse of this process, which can get very tricky, especially as it often does not have a unique inverse.

Answer 4

If you were to only use the standard deviations of the components (the diagonal elements of the covariance), you wouldn't be able to catch some important correlations.

I won't go into the state matrix transformation math here, but try to illustrate the concept with a simple 1-D example: Imagine an object that now is at

x(t=0) = (0 +/- 1)m
v(t=0) = (1 +/- 0.2)m/s

If no acceleration happens, where is it after 100s?

x(t=100) = (100 +/- 20.02)m   [sqrt(1*1 + 20*20)]
v(t=100) = (1 +/- 0.2)m/s

But this isn't the whole truth: the major part of the x uncertainty comes from the initial v uncertainty.

If you just used the individual standard deviations, and plotted the probability density in an x/v graph, you'd see an axis-parallel ellipse 40.04m wide, and 0.4m/s tall, e.g. indicating that x=120,v=1 has a higher probability than x=120,v=1.2, although it should be the other way round (a speed of v=1.2 most likely results in an x=120 after 100s, while with v=1 we'd expect x=100, not 120).

This means that over time, x and v get more and more correlated, and this correlation is found in the non-diagonal elements of the covariance matrix.

And as orbital state is very much about predicting future state from some initially-known situation, the correlation between position and velocity becomes a very important factor, one that is best captured using a covariance matrix.