- Email: [email protected]
Development of an observational error model
Icarus 212 (2011) 438–447
Contents lists available at ScienceDirect
Icarus journal homepage: www.elsevier.com/locate/icarus
Development of an observational error model James Baer a,⇑, Steven R. Chesley b, Andrea Milani c a
James Cook University, School of Engineering and Physical Sciences, Townsville QLD 4811, Australia Jet Propulsion Laboratory/Caltech, 4800 Oak Grove Drive, Pasadena, CA 91109, USA c Dipartimento di Matematica, Universit di Pisa, Largo Pontecorvo 5, 56127 Pisa, Italy b
a r t i c l e
i n f o
Article history: Received 4 September 2010 Revised 18 November 2010 Accepted 18 November 2010 Available online 26 November 2010 Keywords: Asteroids Orbit determination Celestial mechanics
a b s t r a c t In calculating the orbit of a minor planet with a least-squares algorithm, current practice is to assume that all observations of a given era have the same uncertainty, and that the errors in these observations are uncorrelated. These assumptions are unrealistic; and they lead to sub-optimal orbits. Our objective is to develop and validate an observational error model that provides realistic estimates of the uncertainties and correlations in asteroid observations. When used to populate the covariance matrix of the least-squares algorithm, the resulting orbits are shown to more accurately and precisely represent asteroid trajectories. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction
2. The least-squares orbit determination algorithm
Our research interests include asteroid mass determination, Yarkovsky-effect modeling, and near-Earth asteroid collision analysis. These applications require the most accurate and precise possible orbits of subject bodies, in order to detect the smallest perturbations, or to provide reliable estimates of impact probability. Unfortunately, the current convention in least-squares orbit determination is to assume that all observations from a given era are of equivalent precision, regardless of source; and that any errors in these observations are uncorrelated. As we will see, these assumptions make it difficult to simultaneously utilize both heritage observations and more precise sightings from the new surveys; and while such simplifications may have been necessary when computational resources were more limited, they actually penalize the accuracy and precision of current work. In this paper, we will describe the development of a statistical error model of asteroid observations, providing realistic, site-specific estimates of observational uncertainties and error correlations in given eras. We will demonstrate that, by properly weighting observations based on each observatory’s demonstrated performance, the resulting orbits more accurately reflect an asteroid’s true trajectory. We begin with a short summary of current techniques.
The conventional least-squares orbit determination algorithm relies upon finding the minimum of the cost function
⇑ Corresponding author. E-mail addresses: [email protected] (J. Baer), [email protected] (S.R. Chesley). 0019-1035/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2010.11.031
T
QðxÞ ¼ b K1 b In this equation, x is the state vector, b = b(x) is a vector containing the ‘‘observed – computed’’ residuals for each observation, and K is the observational covariance matrix, with
Kij ¼ r ij ri rj where ri is the square root of the covariance (i.e., the uncertainty) of the ith observation, and rij is the correlation between errors in the ith and jth observations. We would find the minimum by seeking stationary points of Q(x):
0¼
dQ db T ¼ 2b K1 dx dx
If we define A ¼ db , then the equation to be solved reduces to dx
0 ¼ AT K1 b If we assume a total of N data points, where N = 2(number of optical observations) + (number of radar delay observations) + (number of radar doppler observations), then b is an N 1 vector, and A is an N 6 matrix. The solution of this equation is
dx ¼ ðAT K1 AÞ1 AT K1 b Conventionally, it is assumed that ri is 3 arcsec for optical observations prior to 1890, 2 arcsec for observations from 1890 to 1950,
439
J. Baer et al. / Icarus 212 (2011) 438–447
Second RA Error Correlation Model − Observatory 568
Second Dec Error Correlation Model − Observatory 568
0.1
0.08 0.07
0.08
0.06 0.06
Correlation
Correlation
0.05 0.04 0.02
0.04 0.03 0.02
0 0.01 −0.02 −0.04
0 0
1
2
3
4
5
−0.01
0
1
Separation Time (days)
2
3
4
5
Separation Time (days) Second Dec Error Correlation Model − Observatory 608
Second RA Error Correlation Model − Observatory 608
0.7
0.45 0.4
0.6
0.35 0.5
Correlation
Correlation
0.3 0.25 0.2 0.15 0.1
0.4 0.3 0.2
0.05 0.1 0 −0.05
0 0
1
2
3
4
5
Separation Time (days)
0
1
2
3
4
5
Separation Time (days)
Fig. 1. RA and Dec correlation models for observatories 568 and 608.
and 1 arcsec thereafter. The additional assumption that observations are uncorrelated, i.e.,
rij ¼
1 when i ¼ j; 0
when i – j:
leads to K being a diagonal matrix, thus significantly simplifying the solution. In an era prior to the widespread availability of modern PCs, these assumptions were not merely defensible, but necessary. However, we must now recognize that these assumptions are both flawed and unnecessary. Factors such as aperture size, detector technology, astrometric reduction methodology, and observer experience are clearly relevant in determining the level of observational precision and accuracy. Combining the astrometric products of diverse instruments also becomes problematic. For instance, Pan-STARRS is expected to yield observational uncertainties of 0.1 arcsec (Jedicke, 2004), almost an order of magnitude more precise than current surveys; but if all observations are weighted with a uniform uncertainty of 1 arcsec, the least-squares algorithm will produce a solution whose errors are equally distributed. It is therefore highly likely that the apparent residuals in Pan-STARRS observations will greatly exceed the design specifications, casting doubt on the system’s performance, and nullifying the usefulness of these highly-precise observations.
Finally, Carpino et al. (2003) has demonstrated that closelyspaced observations of the same asteroid made by the same observatory are significantly correlated; ignoring this correlation has the effect of inaccurately weighting such observations in the leastsquares covariance matrix. In short, subject to these assumptions, the least-squares algorithm will converge to a sub-optimal orbit, with an unrealistic error ellipsoid. Addressing these deficiencies requires the development of time-based error models for all contributing observatories, including realistic values for the bias and RMS errors of each observatory, and a model for the correlation of observations. Carpino et al. (2003) outlined a design for an observational error model. Our intent was to build upon this work, extend it to all available optical observations of the numbered asteroids, validate its performance against conventional methods, and apply it to derive the most accurate possible asteroid orbits. 3. Definition of the error model The observational error model is essentially a statistical analysis of the demonstrated performance of each contributing observatory over a specific time interval. It is explicitly assumed that ‘‘observed – computed’’ residuals are normally distributed (after outlier removal).
440
J. Baer et al. / Icarus 212 (2011) 438–447
Second Dec Error Correlation Model − Observatory 683 0.6
0.25
0.5
0.2
0.4
Correlation
Correlation
Second RA Error Correlation Model − Observatory 683 0.3
0.15 0.1
0.3 0.2
0.05
0.1
0
0
−0.05
0
1
2
3
4
−0.1
5
0
1
2
3
4
5
Separation Time (days)
Separation Time (days) Second RA Error Correlation Model − Observatory 689
Second Dec Error Correlation Model − Observatory 689
0.2
0.25
0.18 0.16
0.2
Correlation
Correlation
0.14 0.12 0.1 0.08
0.15
0.1
0.06 0.04
0.05
0.02 0
0
1
2
3
4
5
Separation Time (days)
0
0
1
2
3
4
5
Separation Time (days)
Fig. 2. RA and Dec correlation models for observatories 683 and 689.
We began by calculating the ‘‘observed – computed’’ residuals for all available observations. To ensure that the orbits were sufficiently well-established as to ensure that these residuals could be computed accurately, consideration was restricted to the 52 million optical observations of numbered asteroids in the Minor Planet Center’s database as of December 12, 2008. The post-1995 observations reduced using the USNO A1.0, A2.0, and B1.0 star catalogs were debiased using the procedure described in Chesley et al. (2010); this helped ensure the resulting population of residuals was unbiased and normally distributed. The debiased observations were then separated into statistically homogeneous ‘‘bins’’, designed to isolate any factors that might reasonably be expected to impact astrometric accuracy or precision; thus, all observations assigned to a particular bin must share the following characteristics: they must be from the same observatory; they must occur within the same 30-day period; they must have a measured or predicted apparent magnitude falling within a one-magnitude range; they must either lie within 10° of the galactic equator, lie above galactic latitude +10°, or lie below galactic latitude 10°; the same detector technology (e.g., CCD vs. photographic plate) must be used to collect the observations.
To satisfy the Law of Large Numbers, we required a minimum bin size of 1000 observations: If a bin was smaller than the minimum size, it was combined with the nearest bin (in the sense of time) that shared all other characteristics. If a bin was still smaller than the minimum size, it was combined with the nearest bin (in the sense of apparent magnitude) that shared all other characteristics. If a bin was still smaller than the minimum size, it was combined with the nearest bin (in the sense of galactic latitude) that shared all other characteristics. If a bin was still smaller than the minimum size, it was combined with a bin that matched in all characteristics except detector technology. If any bins smaller than 1000 observations remained, those observations were placed into a ‘‘MIX’’ bin. Initial attempts indicated that a single MIX bin would be dominated by post-1995 observations; as a result, the uncertainties assigned to observations from earlier eras would be significantly too low. Therefore, we elected to create four MIX bins: MIX1 consists of observations from 1800 to 1890 (the visual era).
441
J. Baer et al. / Icarus 212 (2011) 438–447
Second Dec Error Correlation Model − Observatory 691
0.14
0.14
0.12
0.12
0.1
0.1
0.08
Correlation
Correlation
Second RA Error Correlation Model − Observatory 691
0.06 0.04
0.08 0.06 0.04
0.02 0.02
0 −0.02
0
1
2
3
4
0
5
0
1
3
4
5
Second Dec Error Correlation Model − Observatory 699
Second RA Error Correlation Model − Observatory 699 0.14
0.08
0.12
0.07
0.1
0.06
0.08
Correlation
Correlation
2
Separation Time (days)
Separation Time (days)
0.06 0.04
0.05 0.04 0.03
0.02 0.02
0 −0.02
0
1
2
3
4
5
0.01
0
1
2
3
4
5
Separation Time (days)
Separation Time (days)
Fig. 3. RA and Dec correlation models for observatories 691 and 699.
MIX2 consists of observations from 1890 to 1950 (the early photographic era). MIX3 consists of observations from 1950 to 1995 (the late photographic era). MIX4 consists of observations from 1995 to the present (the CCD era). Within each bin, the sum of the ‘‘observed – computed’’ RA residuals yielded the mean RA bias; the sum of the ‘‘observed – computed’’ Dec residuals yielded the mean Dec bias; the standard deviation of the ‘‘observed – computed’’ RA residuals yielded the observational uncertainty in RA; the standard deviation of the ‘‘observed – computed’’ Dec residuals yielded the observational uncertainty in Dec. Modeling the error correlations was more challenging. As noted previously, Carpino et al. (2003) demonstrated that closely-spaced observations of the same asteroid made by the same observatory had significant non-zero error correlations; errors from observations of different asteroids made by the same observatory were only very weakly correlated. Moreover, our intent was to
model error correlation values for use in calculating the orbit of a single asteroid. Therefore, we restricted our attention to closelyspaced observations of the same asteroid made by the same observatory. Carpino et al. (2003) further noted that such error correlations decayed quickly during the first day, then decayed more slowly to near zero over a period of weeks. We therefore required that observations be separated by less than 50 days. Subject to these requirements, a single bin could not hope to contain enough such pairs of observations as to provide a statistically significant sampling. Therefore, the correlation models were only observatory-specific. Since the first three MIX bins were too sparse to allow the derivation of correlation models, the four MIX bins were combined for this purpose, resulting in a single MIX error correlation model. However, separate biases and RMS errors were calculated for each MIX bin. To accurately model each observatory’s performance, observation pairs were placed in ‘‘baskets’’, depending on their separation in time; baskets were spaced at one hour intervals for the first day, and at daily intervals thereafter. The empirical error correlations within each basket B were calculated as
Corr ¼
1 X bi ri bj rj NB i;j2B
442
J. Baer et al. / Icarus 212 (2011) 438–447
Second RA Error Correlation Model − Observatory 704
Second Dec Error Correlation Model − Observatory 704 0.3
0.2
0.25 0.15
Correlation
Correlation
0.2 0.1
0.05
0.15 0.1 0.05
0 0 −0.05
0
1
2
3
4
−0.05
5
0
1
Separation Time (days)
2
3
4
5
Separation Time (days)
Second RA Error Correlation Model − Observatory MIX
Second Dec Error Correlation Model − Observatory MIX 0.2
0.18
0.18
0.16
0.16 0.14
Correlation
0.1 0.08
0.12 0.1 0.08 0.06
0.06 0.04 0.04 0.02
0.02 0
1
2
3
4
5
0
0
1
2
3
Separation Time (days)
Separation Time (days)
Fig. 4. RA and Dec correlation models for observatories 704 and MIX.
1.3 outside galactic plane inside galactic plane
1.2 1.1
Bin RA RMS Error (Arc Sec)
Correlation
0.14 0.12
1 0.9 0.8 0.7 0.6 0.5 0.4 1995
2000
2005
2010
Bin Start Date Fig. 5. RA RMS errors for observations of magnitude 18–19 made by observatory 704 as a function of time.
4
5
443
J. Baer et al. / Icarus 212 (2011) 438–447
where NB is the number of observation pairs in basket B, bi is the ‘‘observed – computed’’ residual in the ith observation, and ri is the uncertainty in the ith observation. As a purely practical matter, however, the resulting 148 empirical data points for each observatory were too unwieldy to be used in an observational error model; it would be necessary to find a set of mathematical functions, each with a small number of coefficients and scaling factors, that best approximated (in a leastsquares sense) the empirical data curves for each observatory. Moreover, the fact that these mathematical functions would be
used to populate the covariance matrix imposed additional constraints. To guarantee that the covariance matrix can be inverted in the least-squares algorithm, it must be positive definite; and, as explained in Carpino et al. (2003), the individual correlations rij must decay to zero in order to guarantee that the covariance matrix will be positive definite. Therefore, we are limited to using exponential decay functions, or the products of quadratic functions times exponential decay functions, to model the error correlations. In practice, we used functions of the form
1.3 outside galactic plane inside galactic plane
Bin Dec RMS Error (Arc Sec)
1.2
1.1
1
0.9
0.8
0.7
0.6
0.5 1995
2000
2005
2010
Bin Start Date Fig. 6. Dec RMS errors for observations of magnitude 18–19 made by observatory 704 as a function of time.
1.1 app mag 15−16 app mag 19−20
1
RA RMS Error (Arc Sec)
0.9
0.8
0.7
0.6
0.5
0.4
1995
2000
2005
Bin Start Date Fig. 7. RA RMS errors for observatory 704 as a function of bin threshold magnitude.
2010
444
J. Baer et al. / Icarus 212 (2011) 438–447
1 app mag 15−16 app mag 19−20
0.9
Dec RMS Error (Arc Sec)
0.8
0.7
0.6
0.5
0.4
1995
2000
2005
2010
Bin Start Date Fig. 8. Dec RMS errors for observatory 704 as a function of bin threshold magnitude.
0.8 0.7
frequency
0.6 0.5 0.4 0.3 0.2 0.1 0 −3
−2
−1
0
1
2
3
standard deviations
the quadratic-exponential product term was used in a handful of cases where the correlation curves showed additional, non-exponential features. The constants A, B, C, D, E, F, and G were determined in a two-stage process. First, a systematic survey over the range of possible values of these constants revealed a combination that minimized the RMS difference between the observed and modeled correlations; and second, a least-squares algorithm was used to refine the solution. Each bin inherits the RA and Dec correlation models of its observatory; so the error model for each bin consists of the seven correlation model constants in both RA and Dec, combined with the mean biases and RMS errors in RA and Dec. Therefore, one need only determine in which bin an individual observation lies in order to determine the corresponding error model RA and Dec biases, uncertainties, and correlations. The entries of the covariance matrix corresponding to observation pairs made by the same observatory, separated by less than 50 days, are populated using Eq. (1); all other entries are assumed zero.
Fig. 9. Relative differences in error model and debiased epoch state vectors for asteroids 1–200,000.
3.1. Numerical considerations
r ij ¼ AeBt þ CeDt þ Eð1 Ft 2 ÞeGt
In our orbit and mass determinations, we do not directly solve the least squares normal equations
ð1Þ
where t is the time separating observations i and j. The first exponential term was used to represent the rapid decay of error correlations during the first day, while the second exponential term was used to represent the slower subsequent decay through t = 50 days;
ADx ¼ Db Instead, we apply a Square Root Information Filtering (SRIF) algorithm (Bierman, 1974) that requires the covariance matrix to be the identity.
Table 1 Characteristics of the MIX bins. Bin
RA RMS (arcsec)
Dec RMS (arcsec)
RA bias (arcsec)
DEC bias (arcsec)
MIX1 MIX2 MIX3 MIX4
2.976 1.837 0.941 0.545
2.440 1.605 0.919 0.534
0.436 0.120 0.054 0.001
0.037 0.018 0.041 0.029
445
J. Baer et al. / Icarus 212 (2011) 438–447
Cumulative Distribution Function of Nominal RMS Residuals 1 DEM debiased
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.2
0.4
0.6
0.8
1
1.2
Fig. 10. Cumulative Distribution Function: RMS errors for the error model and the conventional model with debiased observations.
Cumulative Distribution Function of Nominal Normalized RMS Residuals 1 DEM debiased
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
Fig. 11. Cumulative Distribution Function: Normalized RMS errors for the error model and the conventional model with debiased observations.
However, the observational error model assumes that the errors between closely-spaced observations made by the same observatory have non-zero correlations; so the covariance matrix is no longer diagonal. We therefore apply the Choleski Square Root algorithm (Burden and Faires, 1985) to find a lower triangular matrix L by which we may pre-multiply A and Db such that the covariance matrix becomes the identity, and the SRIF algorithm may be applied. In some real-world cases, we have discovered that the Choleski algorithm fails. Since the Choleski algorithm is valid if the covariance matrix is positive definite, we can only conclude that the correlation model sometimes results in a covariance matrix that is not positive definite. In theory, this is not possible, since the correlation model is based on functions like AeBt and E(1 Ft2)eGt that decay to zero for large t. In practice, however, estimated pairwise correlations, also known as polychoric correlations, can indeed lead to non-positive definite
covariance matrices, because the correlations are not obtained en masse, and therefore may not be statistically consistent. A straightforward (if numerically brutal) solution is to multiply all of the off-diagonal elements of the covariance matrix by 0.9; this procedure can be repeated until the covariance matrix is positive definite, and the Choleski algorithm converges. After introducing this simple algorithm, the problem has never reappeared.
4. Construction and analysis The observational error model is available to other researchers via the web.1 1
home.earthlink.net/jimbaer1/bins.txt.
446
J. Baer et al. / Icarus 212 (2011) 438–447
All MPC-listed observatories (and thus all observations in the MPC database) were included in the error model. We attempted to create correlation models for as many observatories as possible. As described in Section 3, however, very large data sets are required to find enough observation pairs to populate the 148 correlation baskets. In the end, correlation models were created for the following observatories:
106 – Crni Vrh 291 – LPL/Spacewatch II 568 – Mauna Kea 608 – Haleakala-AMOS 644 – Palomar Mountain/NEAT 673 – Table Mountain Observatory, Wrightwood 683 – Goodricke-Pigott Observatory, Tucson 689 – US Naval Observatory, Flagstaff 691 – Steward Observatory, Kitt Peak-Spacewatch 699 – Lowell Observatory-LONEOS 703 – Catalina Sky Survey 704 – Lincoln Laboratory ETS, New Mexico E12 – Siding Spring Survey G96 – Mt. Lemmon Survey J75 – OAM Observatory, La Sagra
those resulting from the conventional default historical uncertainties with only catalog debiasing applied. While the majority of differences are less than 1r, there is a non-trivial proportion of cases where the differences are significantly greater. Finally, as Table 1 illustrates, the RA and Dec RMS errors for the MIX1, MIX2, and MIX3 bins closely match the corresponding default historical values. While not endorsing the default assumptions, we find this result reassuring. Indeed, the default historical values were not arbitrary, but rather based upon extensive experience in orbit determination; any error model that yielded terribly different mean values for these eras would have been suspect. Moreover, this result further illustrates the advantages of our bin design, as new technology leads to smaller residuals over time.
5. Error model validation To validate performance, the orbits of the numbered asteroids were calculated both using the error model, and using conven-
Cumulative Distribution Function of Opposition RMS errors 1 DEM
Several of these correlation models are illustrated in Figs. 1–4. Notice that the error correlations are quite small, and quickly decay to zero. This is precisely as expected. And since the correlations were significant only for approximately the first five days, particular care was taken in fitting the correlation models to the data in this period. Analysis of the bin-specific residuals in the error model reveal several significant features. First, as noted in Chesley et al. (2010, Icarus, in press), more than 60% of observations in the MPC database were collected by the modern survey observatories 699, 703, 704, and G96. Perusal of their entries in the error model reveals that their RA and Dec RMS errors are significantly less than 1 arcsec, typically ranging from 0.2 to 0.8 arcsec. Absent this statistical treatment, observations from these observatories would have been assigned unrealistically high uncertainties, thus deweighting these observations in the least squares solutions. Not only would the resulting uncertainty ellipsoids have been too large, but the solutions themselves would have been sub-optimal, favoring less precise observations. Second, as illustrated by Figs. 5 and 6, the RA and Dec RMS errors for observations of apparent magnitude 18–19 made by observatory 704 decrease significantly with time. This is not unexpected; a survey’s initial test observations are likely to be less precise, as both equipment and software are calibrated and refined for operational use. Moreover, these figures also demonstrate that observations more than 10° outside the galactic plane have significantly higher RMS errors than those within the plane, likely due to the relative absence of background stars from which the asteroid’s position is deduced. Absent this statistical treatment, however, these observations would have been assigned a uniform uncertainty, likely overweighting the less-precise test observations (and those outside the galactic plane), and likely underweighting the more-precise operational observations (and those near the galactic plane). Third, as illustrated by Figs. 7 and 8, the RA and Dec RMS errors for observatory 704 increase as the threshold magnitude of the bin increases. Again, this is not unexpected; as an asteroid becomes fainter, the signal-to-noise ratio drops, making it more difficult to precisely locate the image centroid. However, these plots clearly illustrate that failing to account for this effect again results in improper observational weighting, sub-optimal orbits, and inaccurate uncertainty ellipsoids. Fourth, Fig. 9 illustrates the relative differences for asteroids 1– 200,000 between error model epoch state vector components, and
debiased 0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
Fig. 12. Cumulative Distribution Function: RMS errors for the error model and the conventional model with debiased observations for opposition n 3.
Cumulative Distribution Function of Opposition RMS errors 0.9 debiased error model
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
RMS Error (arc sec) Fig. 13. Cumulative Distribution Function: RMS errors for the error model and the conventional model with debiased observations for opposition n 4.
J. Baer et al. / Icarus 212 (2011) 438–447
tional default historical uncertainties with only catalog debiasing applied. The residuals and RMS errors from these two sets of orbits were compared. In 99.1% of cases, the error model RMS errors were smaller, with a mean improvement of 0.118 arcsec. Additionally, the normalized error model RMS errors (compiled by dividing observational residuals by their corresponding assumed uncertainties) were closer to 1 in 98.2% of cases. These results are illustrated in Figs. 10 and 11. Obviously, these are very desirable characteristics, as they imply that the error model orbits agree with observations significantly better, and that the uncertainty ellipsoids are more realistic. Nevertheless, these results are insufficient to validate the error model. By selectively excluding observations, it is possible for the least-squares algorithm to converge to a minimum residual solution that is nonetheless physically inaccurate. In order to prove that the orbits calculated using the error model are superior, we must demonstrate that they better predict an asteroid’s position at times outside of the observational baseline. Assuming that an asteroid was observed at n oppositions, we examined two scenarios: We selected 4530 numbered asteroids that had 50 or more observations spread over at least 50 days in both oppositions n 1 and n 2. We calculated the error model and debiased orbits based on observations from those two oppositions, and used these orbits to calculate the ‘‘observed – predicted’’ residuals for the observations of opposition n 3. The results are illustrated in Fig. 12, which shows that the error model orbits predict the n 3 observations better, with a mean improvement of 0.091 arcsec. We selected 7250 numbered asteroids that had 30 or more observations spread over at least 30 days in each of oppositions n 1, n 2, and n 3. We calculated the error model and debiased orbits based on observations from those three oppositions, and used these orbits to calculate the ‘‘observed – predicted’’ residuals for the observations of opposition n 4. The results are illustrated in Fig. 13, which shows that the error model orbits predict the n 4 observations better, with a mean improvement of 0.080 arcsec. We therefore conclude that the observational error model is statistically valid, in that it more accurately models the physical trajectory of asteroids, and their corresponding region of uncertainty. 6. Interpretation It is important to remember that this is a statistical error model, based upon Gaussian probability distributions. The resulting least-
447
squares standard deviations should therefore not be interpreted as absolute uncertainties, but rather as 1r uncertainties. This change in perception is crucial to placing such high-precision applications as asteroid impact analysis, Yarkovsky-effect modeling, and mass determination onto a firm statistical foundation. Assume, for instance, that we are attempting to recover a newly-discovered asteroid. Then there is a 32% probability that the asteroid will lie outside the 1r uncertainty ellipsoid, a 5% probability that the asteroid will lie outside the 2r uncertainty ellipsoid, and a 0.3% probability that the asteroid will lie outside the 3r uncertainty ellipsoid. The observer must calibrate the size of the search region against the risk they are willing to accept of missing the target. On the other hand, the user can rest assured that, on average, the asteroid will indeed lie within the 1r uncertainty ellipsoid 68% of the time, etc. 7. Future plans With the construction and performance of the observational error model validated, we look forward to its operational implementation. However, it will require regular maintenance, both to accommodate the introduction of new observatories, and to account for changes in performance in existing observatories, as new technology, software, and procedures are adopted. Additionally, there is a large volume of photometric asteroid observations. The development of a photometric error model, similar in principle to the astrometric model developed herein, could potentially enable the development of more precise asteroid shape models, and thus more precise volume, density, and porosity estimates. References Bierman, G., 1974. Sequential square root filtering and smoothing of discrete linear systems. Automatica 10, 147–158. Burden, R.L., Faires, J.D., 1985. Numerical Analysis, third ed. Weber and Schmidt, Prindle. Carpino, M., Milani, A., Chesley, S.R., 2003. Error statistics of asteroid optical astrometric observations. Icarus 166, 248–270. Chesley, S., Baer, J., Monet, D., 2010. Treatment of star catalog biases in asteroid astrometric observations. Icarus 210, 158–181. Jedicke, R., and The Pan-Starrs Collaboration, 2004. The Pan-STARRS Solar System survey. In 35th COSPAR Scientific Assembly, COSPAR, Plenary Meeting, vol. 35, p. 999.
Contents lists available at ScienceDirect
Icarus journal homepage: www.elsevier.com/locate/icarus
Development of an observational error model James Baer a,⇑, Steven R. Chesley b, Andrea Milani c a
James Cook University, School of Engineering and Physical Sciences, Townsville QLD 4811, Australia Jet Propulsion Laboratory/Caltech, 4800 Oak Grove Drive, Pasadena, CA 91109, USA c Dipartimento di Matematica, Universit di Pisa, Largo Pontecorvo 5, 56127 Pisa, Italy b
a r t i c l e
i n f o
Article history: Received 4 September 2010 Revised 18 November 2010 Accepted 18 November 2010 Available online 26 November 2010 Keywords: Asteroids Orbit determination Celestial mechanics
a b s t r a c t In calculating the orbit of a minor planet with a least-squares algorithm, current practice is to assume that all observations of a given era have the same uncertainty, and that the errors in these observations are uncorrelated. These assumptions are unrealistic; and they lead to sub-optimal orbits. Our objective is to develop and validate an observational error model that provides realistic estimates of the uncertainties and correlations in asteroid observations. When used to populate the covariance matrix of the least-squares algorithm, the resulting orbits are shown to more accurately and precisely represent asteroid trajectories. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction
2. The least-squares orbit determination algorithm
Our research interests include asteroid mass determination, Yarkovsky-effect modeling, and near-Earth asteroid collision analysis. These applications require the most accurate and precise possible orbits of subject bodies, in order to detect the smallest perturbations, or to provide reliable estimates of impact probability. Unfortunately, the current convention in least-squares orbit determination is to assume that all observations from a given era are of equivalent precision, regardless of source; and that any errors in these observations are uncorrelated. As we will see, these assumptions make it difficult to simultaneously utilize both heritage observations and more precise sightings from the new surveys; and while such simplifications may have been necessary when computational resources were more limited, they actually penalize the accuracy and precision of current work. In this paper, we will describe the development of a statistical error model of asteroid observations, providing realistic, site-specific estimates of observational uncertainties and error correlations in given eras. We will demonstrate that, by properly weighting observations based on each observatory’s demonstrated performance, the resulting orbits more accurately reflect an asteroid’s true trajectory. We begin with a short summary of current techniques.
The conventional least-squares orbit determination algorithm relies upon finding the minimum of the cost function
⇑ Corresponding author. E-mail addresses: [email protected] (J. Baer), [email protected] (S.R. Chesley). 0019-1035/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2010.11.031
T
QðxÞ ¼ b K1 b In this equation, x is the state vector, b = b(x) is a vector containing the ‘‘observed – computed’’ residuals for each observation, and K is the observational covariance matrix, with
Kij ¼ r ij ri rj where ri is the square root of the covariance (i.e., the uncertainty) of the ith observation, and rij is the correlation between errors in the ith and jth observations. We would find the minimum by seeking stationary points of Q(x):
0¼
dQ db T ¼ 2b K1 dx dx
If we define A ¼ db , then the equation to be solved reduces to dx
0 ¼ AT K1 b If we assume a total of N data points, where N = 2(number of optical observations) + (number of radar delay observations) + (number of radar doppler observations), then b is an N 1 vector, and A is an N 6 matrix. The solution of this equation is
dx ¼ ðAT K1 AÞ1 AT K1 b Conventionally, it is assumed that ri is 3 arcsec for optical observations prior to 1890, 2 arcsec for observations from 1890 to 1950,
439
J. Baer et al. / Icarus 212 (2011) 438–447
Second RA Error Correlation Model − Observatory 568
Second Dec Error Correlation Model − Observatory 568
0.1
0.08 0.07
0.08
0.06 0.06
Correlation
Correlation
0.05 0.04 0.02
0.04 0.03 0.02
0 0.01 −0.02 −0.04
0 0
1
2
3
4
5
−0.01
0
1
Separation Time (days)
2
3
4
5
Separation Time (days) Second Dec Error Correlation Model − Observatory 608
Second RA Error Correlation Model − Observatory 608
0.7
0.45 0.4
0.6
0.35 0.5
Correlation
Correlation
0.3 0.25 0.2 0.15 0.1
0.4 0.3 0.2
0.05 0.1 0 −0.05
0 0
1
2
3
4
5
Separation Time (days)
0
1
2
3
4
5
Separation Time (days)
Fig. 1. RA and Dec correlation models for observatories 568 and 608.
and 1 arcsec thereafter. The additional assumption that observations are uncorrelated, i.e.,
rij ¼
1 when i ¼ j; 0
when i – j:
leads to K being a diagonal matrix, thus significantly simplifying the solution. In an era prior to the widespread availability of modern PCs, these assumptions were not merely defensible, but necessary. However, we must now recognize that these assumptions are both flawed and unnecessary. Factors such as aperture size, detector technology, astrometric reduction methodology, and observer experience are clearly relevant in determining the level of observational precision and accuracy. Combining the astrometric products of diverse instruments also becomes problematic. For instance, Pan-STARRS is expected to yield observational uncertainties of 0.1 arcsec (Jedicke, 2004), almost an order of magnitude more precise than current surveys; but if all observations are weighted with a uniform uncertainty of 1 arcsec, the least-squares algorithm will produce a solution whose errors are equally distributed. It is therefore highly likely that the apparent residuals in Pan-STARRS observations will greatly exceed the design specifications, casting doubt on the system’s performance, and nullifying the usefulness of these highly-precise observations.
Finally, Carpino et al. (2003) has demonstrated that closelyspaced observations of the same asteroid made by the same observatory are significantly correlated; ignoring this correlation has the effect of inaccurately weighting such observations in the leastsquares covariance matrix. In short, subject to these assumptions, the least-squares algorithm will converge to a sub-optimal orbit, with an unrealistic error ellipsoid. Addressing these deficiencies requires the development of time-based error models for all contributing observatories, including realistic values for the bias and RMS errors of each observatory, and a model for the correlation of observations. Carpino et al. (2003) outlined a design for an observational error model. Our intent was to build upon this work, extend it to all available optical observations of the numbered asteroids, validate its performance against conventional methods, and apply it to derive the most accurate possible asteroid orbits. 3. Definition of the error model The observational error model is essentially a statistical analysis of the demonstrated performance of each contributing observatory over a specific time interval. It is explicitly assumed that ‘‘observed – computed’’ residuals are normally distributed (after outlier removal).
440
J. Baer et al. / Icarus 212 (2011) 438–447
Second Dec Error Correlation Model − Observatory 683 0.6
0.25
0.5
0.2
0.4
Correlation
Correlation
Second RA Error Correlation Model − Observatory 683 0.3
0.15 0.1
0.3 0.2
0.05
0.1
0
0
−0.05
0
1
2
3
4
−0.1
5
0
1
2
3
4
5
Separation Time (days)
Separation Time (days) Second RA Error Correlation Model − Observatory 689
Second Dec Error Correlation Model − Observatory 689
0.2
0.25
0.18 0.16
0.2
Correlation
Correlation
0.14 0.12 0.1 0.08
0.15
0.1
0.06 0.04
0.05
0.02 0
0
1
2
3
4
5
Separation Time (days)
0
0
1
2
3
4
5
Separation Time (days)
Fig. 2. RA and Dec correlation models for observatories 683 and 689.
We began by calculating the ‘‘observed – computed’’ residuals for all available observations. To ensure that the orbits were sufficiently well-established as to ensure that these residuals could be computed accurately, consideration was restricted to the 52 million optical observations of numbered asteroids in the Minor Planet Center’s database as of December 12, 2008. The post-1995 observations reduced using the USNO A1.0, A2.0, and B1.0 star catalogs were debiased using the procedure described in Chesley et al. (2010); this helped ensure the resulting population of residuals was unbiased and normally distributed. The debiased observations were then separated into statistically homogeneous ‘‘bins’’, designed to isolate any factors that might reasonably be expected to impact astrometric accuracy or precision; thus, all observations assigned to a particular bin must share the following characteristics: they must be from the same observatory; they must occur within the same 30-day period; they must have a measured or predicted apparent magnitude falling within a one-magnitude range; they must either lie within 10° of the galactic equator, lie above galactic latitude +10°, or lie below galactic latitude 10°; the same detector technology (e.g., CCD vs. photographic plate) must be used to collect the observations.
To satisfy the Law of Large Numbers, we required a minimum bin size of 1000 observations: If a bin was smaller than the minimum size, it was combined with the nearest bin (in the sense of time) that shared all other characteristics. If a bin was still smaller than the minimum size, it was combined with the nearest bin (in the sense of apparent magnitude) that shared all other characteristics. If a bin was still smaller than the minimum size, it was combined with the nearest bin (in the sense of galactic latitude) that shared all other characteristics. If a bin was still smaller than the minimum size, it was combined with a bin that matched in all characteristics except detector technology. If any bins smaller than 1000 observations remained, those observations were placed into a ‘‘MIX’’ bin. Initial attempts indicated that a single MIX bin would be dominated by post-1995 observations; as a result, the uncertainties assigned to observations from earlier eras would be significantly too low. Therefore, we elected to create four MIX bins: MIX1 consists of observations from 1800 to 1890 (the visual era).
441
J. Baer et al. / Icarus 212 (2011) 438–447
Second Dec Error Correlation Model − Observatory 691
0.14
0.14
0.12
0.12
0.1
0.1
0.08
Correlation
Correlation
Second RA Error Correlation Model − Observatory 691
0.06 0.04
0.08 0.06 0.04
0.02 0.02
0 −0.02
0
1
2
3
4
0
5
0
1
3
4
5
Second Dec Error Correlation Model − Observatory 699
Second RA Error Correlation Model − Observatory 699 0.14
0.08
0.12
0.07
0.1
0.06
0.08
Correlation
Correlation
2
Separation Time (days)
Separation Time (days)
0.06 0.04
0.05 0.04 0.03
0.02 0.02
0 −0.02
0
1
2
3
4
5
0.01
0
1
2
3
4
5
Separation Time (days)
Separation Time (days)
Fig. 3. RA and Dec correlation models for observatories 691 and 699.
MIX2 consists of observations from 1890 to 1950 (the early photographic era). MIX3 consists of observations from 1950 to 1995 (the late photographic era). MIX4 consists of observations from 1995 to the present (the CCD era). Within each bin, the sum of the ‘‘observed – computed’’ RA residuals yielded the mean RA bias; the sum of the ‘‘observed – computed’’ Dec residuals yielded the mean Dec bias; the standard deviation of the ‘‘observed – computed’’ RA residuals yielded the observational uncertainty in RA; the standard deviation of the ‘‘observed – computed’’ Dec residuals yielded the observational uncertainty in Dec. Modeling the error correlations was more challenging. As noted previously, Carpino et al. (2003) demonstrated that closely-spaced observations of the same asteroid made by the same observatory had significant non-zero error correlations; errors from observations of different asteroids made by the same observatory were only very weakly correlated. Moreover, our intent was to
model error correlation values for use in calculating the orbit of a single asteroid. Therefore, we restricted our attention to closelyspaced observations of the same asteroid made by the same observatory. Carpino et al. (2003) further noted that such error correlations decayed quickly during the first day, then decayed more slowly to near zero over a period of weeks. We therefore required that observations be separated by less than 50 days. Subject to these requirements, a single bin could not hope to contain enough such pairs of observations as to provide a statistically significant sampling. Therefore, the correlation models were only observatory-specific. Since the first three MIX bins were too sparse to allow the derivation of correlation models, the four MIX bins were combined for this purpose, resulting in a single MIX error correlation model. However, separate biases and RMS errors were calculated for each MIX bin. To accurately model each observatory’s performance, observation pairs were placed in ‘‘baskets’’, depending on their separation in time; baskets were spaced at one hour intervals for the first day, and at daily intervals thereafter. The empirical error correlations within each basket B were calculated as
Corr ¼
1 X bi ri bj rj NB i;j2B
442
J. Baer et al. / Icarus 212 (2011) 438–447
Second RA Error Correlation Model − Observatory 704
Second Dec Error Correlation Model − Observatory 704 0.3
0.2
0.25 0.15
Correlation
Correlation
0.2 0.1
0.05
0.15 0.1 0.05
0 0 −0.05
0
1
2
3
4
−0.05
5
0
1
Separation Time (days)
2
3
4
5
Separation Time (days)
Second RA Error Correlation Model − Observatory MIX
Second Dec Error Correlation Model − Observatory MIX 0.2
0.18
0.18
0.16
0.16 0.14
Correlation
0.1 0.08
0.12 0.1 0.08 0.06
0.06 0.04 0.04 0.02
0.02 0
1
2
3
4
5
0
0
1
2
3
Separation Time (days)
Separation Time (days)
Fig. 4. RA and Dec correlation models for observatories 704 and MIX.
1.3 outside galactic plane inside galactic plane
1.2 1.1
Bin RA RMS Error (Arc Sec)
Correlation
0.14 0.12
1 0.9 0.8 0.7 0.6 0.5 0.4 1995
2000
2005
2010
Bin Start Date Fig. 5. RA RMS errors for observations of magnitude 18–19 made by observatory 704 as a function of time.
4
5
443
J. Baer et al. / Icarus 212 (2011) 438–447
where NB is the number of observation pairs in basket B, bi is the ‘‘observed – computed’’ residual in the ith observation, and ri is the uncertainty in the ith observation. As a purely practical matter, however, the resulting 148 empirical data points for each observatory were too unwieldy to be used in an observational error model; it would be necessary to find a set of mathematical functions, each with a small number of coefficients and scaling factors, that best approximated (in a leastsquares sense) the empirical data curves for each observatory. Moreover, the fact that these mathematical functions would be
used to populate the covariance matrix imposed additional constraints. To guarantee that the covariance matrix can be inverted in the least-squares algorithm, it must be positive definite; and, as explained in Carpino et al. (2003), the individual correlations rij must decay to zero in order to guarantee that the covariance matrix will be positive definite. Therefore, we are limited to using exponential decay functions, or the products of quadratic functions times exponential decay functions, to model the error correlations. In practice, we used functions of the form
1.3 outside galactic plane inside galactic plane
Bin Dec RMS Error (Arc Sec)
1.2
1.1
1
0.9
0.8
0.7
0.6
0.5 1995
2000
2005
2010
Bin Start Date Fig. 6. Dec RMS errors for observations of magnitude 18–19 made by observatory 704 as a function of time.
1.1 app mag 15−16 app mag 19−20
1
RA RMS Error (Arc Sec)
0.9
0.8
0.7
0.6
0.5
0.4
1995
2000
2005
Bin Start Date Fig. 7. RA RMS errors for observatory 704 as a function of bin threshold magnitude.
2010
444
J. Baer et al. / Icarus 212 (2011) 438–447
1 app mag 15−16 app mag 19−20
0.9
Dec RMS Error (Arc Sec)
0.8
0.7
0.6
0.5
0.4
1995
2000
2005
2010
Bin Start Date Fig. 8. Dec RMS errors for observatory 704 as a function of bin threshold magnitude.
0.8 0.7
frequency
0.6 0.5 0.4 0.3 0.2 0.1 0 −3
−2
−1
0
1
2
3
standard deviations
the quadratic-exponential product term was used in a handful of cases where the correlation curves showed additional, non-exponential features. The constants A, B, C, D, E, F, and G were determined in a two-stage process. First, a systematic survey over the range of possible values of these constants revealed a combination that minimized the RMS difference between the observed and modeled correlations; and second, a least-squares algorithm was used to refine the solution. Each bin inherits the RA and Dec correlation models of its observatory; so the error model for each bin consists of the seven correlation model constants in both RA and Dec, combined with the mean biases and RMS errors in RA and Dec. Therefore, one need only determine in which bin an individual observation lies in order to determine the corresponding error model RA and Dec biases, uncertainties, and correlations. The entries of the covariance matrix corresponding to observation pairs made by the same observatory, separated by less than 50 days, are populated using Eq. (1); all other entries are assumed zero.
Fig. 9. Relative differences in error model and debiased epoch state vectors for asteroids 1–200,000.
3.1. Numerical considerations
r ij ¼ AeBt þ CeDt þ Eð1 Ft 2 ÞeGt
In our orbit and mass determinations, we do not directly solve the least squares normal equations
ð1Þ
where t is the time separating observations i and j. The first exponential term was used to represent the rapid decay of error correlations during the first day, while the second exponential term was used to represent the slower subsequent decay through t = 50 days;
ADx ¼ Db Instead, we apply a Square Root Information Filtering (SRIF) algorithm (Bierman, 1974) that requires the covariance matrix to be the identity.
Table 1 Characteristics of the MIX bins. Bin
RA RMS (arcsec)
Dec RMS (arcsec)
RA bias (arcsec)
DEC bias (arcsec)
MIX1 MIX2 MIX3 MIX4
2.976 1.837 0.941 0.545
2.440 1.605 0.919 0.534
0.436 0.120 0.054 0.001
0.037 0.018 0.041 0.029
445
J. Baer et al. / Icarus 212 (2011) 438–447
Cumulative Distribution Function of Nominal RMS Residuals 1 DEM debiased
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.2
0.4
0.6
0.8
1
1.2
Fig. 10. Cumulative Distribution Function: RMS errors for the error model and the conventional model with debiased observations.
Cumulative Distribution Function of Nominal Normalized RMS Residuals 1 DEM debiased
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
Fig. 11. Cumulative Distribution Function: Normalized RMS errors for the error model and the conventional model with debiased observations.
However, the observational error model assumes that the errors between closely-spaced observations made by the same observatory have non-zero correlations; so the covariance matrix is no longer diagonal. We therefore apply the Choleski Square Root algorithm (Burden and Faires, 1985) to find a lower triangular matrix L by which we may pre-multiply A and Db such that the covariance matrix becomes the identity, and the SRIF algorithm may be applied. In some real-world cases, we have discovered that the Choleski algorithm fails. Since the Choleski algorithm is valid if the covariance matrix is positive definite, we can only conclude that the correlation model sometimes results in a covariance matrix that is not positive definite. In theory, this is not possible, since the correlation model is based on functions like AeBt and E(1 Ft2)eGt that decay to zero for large t. In practice, however, estimated pairwise correlations, also known as polychoric correlations, can indeed lead to non-positive definite
covariance matrices, because the correlations are not obtained en masse, and therefore may not be statistically consistent. A straightforward (if numerically brutal) solution is to multiply all of the off-diagonal elements of the covariance matrix by 0.9; this procedure can be repeated until the covariance matrix is positive definite, and the Choleski algorithm converges. After introducing this simple algorithm, the problem has never reappeared.
4. Construction and analysis The observational error model is available to other researchers via the web.1 1
home.earthlink.net/jimbaer1/bins.txt.
446
J. Baer et al. / Icarus 212 (2011) 438–447
All MPC-listed observatories (and thus all observations in the MPC database) were included in the error model. We attempted to create correlation models for as many observatories as possible. As described in Section 3, however, very large data sets are required to find enough observation pairs to populate the 148 correlation baskets. In the end, correlation models were created for the following observatories:
106 – Crni Vrh 291 – LPL/Spacewatch II 568 – Mauna Kea 608 – Haleakala-AMOS 644 – Palomar Mountain/NEAT 673 – Table Mountain Observatory, Wrightwood 683 – Goodricke-Pigott Observatory, Tucson 689 – US Naval Observatory, Flagstaff 691 – Steward Observatory, Kitt Peak-Spacewatch 699 – Lowell Observatory-LONEOS 703 – Catalina Sky Survey 704 – Lincoln Laboratory ETS, New Mexico E12 – Siding Spring Survey G96 – Mt. Lemmon Survey J75 – OAM Observatory, La Sagra
those resulting from the conventional default historical uncertainties with only catalog debiasing applied. While the majority of differences are less than 1r, there is a non-trivial proportion of cases where the differences are significantly greater. Finally, as Table 1 illustrates, the RA and Dec RMS errors for the MIX1, MIX2, and MIX3 bins closely match the corresponding default historical values. While not endorsing the default assumptions, we find this result reassuring. Indeed, the default historical values were not arbitrary, but rather based upon extensive experience in orbit determination; any error model that yielded terribly different mean values for these eras would have been suspect. Moreover, this result further illustrates the advantages of our bin design, as new technology leads to smaller residuals over time.
5. Error model validation To validate performance, the orbits of the numbered asteroids were calculated both using the error model, and using conven-
Cumulative Distribution Function of Opposition RMS errors 1 DEM
Several of these correlation models are illustrated in Figs. 1–4. Notice that the error correlations are quite small, and quickly decay to zero. This is precisely as expected. And since the correlations were significant only for approximately the first five days, particular care was taken in fitting the correlation models to the data in this period. Analysis of the bin-specific residuals in the error model reveal several significant features. First, as noted in Chesley et al. (2010, Icarus, in press), more than 60% of observations in the MPC database were collected by the modern survey observatories 699, 703, 704, and G96. Perusal of their entries in the error model reveals that their RA and Dec RMS errors are significantly less than 1 arcsec, typically ranging from 0.2 to 0.8 arcsec. Absent this statistical treatment, observations from these observatories would have been assigned unrealistically high uncertainties, thus deweighting these observations in the least squares solutions. Not only would the resulting uncertainty ellipsoids have been too large, but the solutions themselves would have been sub-optimal, favoring less precise observations. Second, as illustrated by Figs. 5 and 6, the RA and Dec RMS errors for observations of apparent magnitude 18–19 made by observatory 704 decrease significantly with time. This is not unexpected; a survey’s initial test observations are likely to be less precise, as both equipment and software are calibrated and refined for operational use. Moreover, these figures also demonstrate that observations more than 10° outside the galactic plane have significantly higher RMS errors than those within the plane, likely due to the relative absence of background stars from which the asteroid’s position is deduced. Absent this statistical treatment, however, these observations would have been assigned a uniform uncertainty, likely overweighting the less-precise test observations (and those outside the galactic plane), and likely underweighting the more-precise operational observations (and those near the galactic plane). Third, as illustrated by Figs. 7 and 8, the RA and Dec RMS errors for observatory 704 increase as the threshold magnitude of the bin increases. Again, this is not unexpected; as an asteroid becomes fainter, the signal-to-noise ratio drops, making it more difficult to precisely locate the image centroid. However, these plots clearly illustrate that failing to account for this effect again results in improper observational weighting, sub-optimal orbits, and inaccurate uncertainty ellipsoids. Fourth, Fig. 9 illustrates the relative differences for asteroids 1– 200,000 between error model epoch state vector components, and
debiased 0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
Fig. 12. Cumulative Distribution Function: RMS errors for the error model and the conventional model with debiased observations for opposition n 3.
Cumulative Distribution Function of Opposition RMS errors 0.9 debiased error model
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.5
1
1.5
2
2.5
3
RMS Error (arc sec) Fig. 13. Cumulative Distribution Function: RMS errors for the error model and the conventional model with debiased observations for opposition n 4.
J. Baer et al. / Icarus 212 (2011) 438–447
tional default historical uncertainties with only catalog debiasing applied. The residuals and RMS errors from these two sets of orbits were compared. In 99.1% of cases, the error model RMS errors were smaller, with a mean improvement of 0.118 arcsec. Additionally, the normalized error model RMS errors (compiled by dividing observational residuals by their corresponding assumed uncertainties) were closer to 1 in 98.2% of cases. These results are illustrated in Figs. 10 and 11. Obviously, these are very desirable characteristics, as they imply that the error model orbits agree with observations significantly better, and that the uncertainty ellipsoids are more realistic. Nevertheless, these results are insufficient to validate the error model. By selectively excluding observations, it is possible for the least-squares algorithm to converge to a minimum residual solution that is nonetheless physically inaccurate. In order to prove that the orbits calculated using the error model are superior, we must demonstrate that they better predict an asteroid’s position at times outside of the observational baseline. Assuming that an asteroid was observed at n oppositions, we examined two scenarios: We selected 4530 numbered asteroids that had 50 or more observations spread over at least 50 days in both oppositions n 1 and n 2. We calculated the error model and debiased orbits based on observations from those two oppositions, and used these orbits to calculate the ‘‘observed – predicted’’ residuals for the observations of opposition n 3. The results are illustrated in Fig. 12, which shows that the error model orbits predict the n 3 observations better, with a mean improvement of 0.091 arcsec. We selected 7250 numbered asteroids that had 30 or more observations spread over at least 30 days in each of oppositions n 1, n 2, and n 3. We calculated the error model and debiased orbits based on observations from those three oppositions, and used these orbits to calculate the ‘‘observed – predicted’’ residuals for the observations of opposition n 4. The results are illustrated in Fig. 13, which shows that the error model orbits predict the n 4 observations better, with a mean improvement of 0.080 arcsec. We therefore conclude that the observational error model is statistically valid, in that it more accurately models the physical trajectory of asteroids, and their corresponding region of uncertainty. 6. Interpretation It is important to remember that this is a statistical error model, based upon Gaussian probability distributions. The resulting least-
447
squares standard deviations should therefore not be interpreted as absolute uncertainties, but rather as 1r uncertainties. This change in perception is crucial to placing such high-precision applications as asteroid impact analysis, Yarkovsky-effect modeling, and mass determination onto a firm statistical foundation. Assume, for instance, that we are attempting to recover a newly-discovered asteroid. Then there is a 32% probability that the asteroid will lie outside the 1r uncertainty ellipsoid, a 5% probability that the asteroid will lie outside the 2r uncertainty ellipsoid, and a 0.3% probability that the asteroid will lie outside the 3r uncertainty ellipsoid. The observer must calibrate the size of the search region against the risk they are willing to accept of missing the target. On the other hand, the user can rest assured that, on average, the asteroid will indeed lie within the 1r uncertainty ellipsoid 68% of the time, etc. 7. Future plans With the construction and performance of the observational error model validated, we look forward to its operational implementation. However, it will require regular maintenance, both to accommodate the introduction of new observatories, and to account for changes in performance in existing observatories, as new technology, software, and procedures are adopted. Additionally, there is a large volume of photometric asteroid observations. The development of a photometric error model, similar in principle to the astrometric model developed herein, could potentially enable the development of more precise asteroid shape models, and thus more precise volume, density, and porosity estimates. References Bierman, G., 1974. Sequential square root filtering and smoothing of discrete linear systems. Automatica 10, 147–158. Burden, R.L., Faires, J.D., 1985. Numerical Analysis, third ed. Weber and Schmidt, Prindle. Carpino, M., Milani, A., Chesley, S.R., 2003. Error statistics of asteroid optical astrometric observations. Icarus 166, 248–270. Chesley, S., Baer, J., Monet, D., 2010. Treatment of star catalog biases in asteroid astrometric observations. Icarus 210, 158–181. Jedicke, R., and The Pan-Starrs Collaboration, 2004. The Pan-STARRS Solar System survey. In 35th COSPAR Scientific Assembly, COSPAR, Plenary Meeting, vol. 35, p. 999.