PaperlayerView exact revision on arXiv

Source-linked AI summary

The astrometric core solution for the Gaia mission. Overview of models, algorithms and software implementation

Lennart Lindegren, Uwe Lammers, David Hobbs, William O'Mullane, Ulrich Bastian, José Hernández

arXiv:1112.4139v1astro-ph.IM

TL;DR

Gaia requires a global astrometric solution that can estimate source parameters together with satellite attitude and instrument calibration at unprecedented scale. The paper develops the mathematical models, iterative algorithms, and AGIS software framework for this task, reporting successful simulated-data processing and a feasible path to an accurate mission solution. Its scope depends on selected primary sources and models that may need refinement with flight data.

  • Problem

    Gaia’s astrometric core solution must handle roughly half a billion coupled unknowns and enormous, entangled observation datasets while preserving the mission’s astrometric accuracy.

  • Method

    The paper formulates a global weighted least-squares problem and solves it iteratively through source, attitude, calibration, and global-processing components implemented in AGIS.

  • Results

    A first AGIS run using simulated five-year data with 50 million primary sources was successfully completed, and test runs demonstrated the global iterative approach’s theoretical validity and practical feasibility.

  • Takeaways & Limitations

    The astrometric core solution supplies primary-source astrometry, the instrument attitude reference frame, and geometric calibration used by other Gaia scientific-reduction processes.

Abstract

from arXiv · show

The Gaia satellite will observe about one billion stars and other point-like sources. The astrometric core solution will determine the astrometric parameters (position, parallax, and proper motion) for a subset of these sources, using a global solution approach which must also include a large number of parameters for the satellite attitude and optical instrument. The accurate and efficient implementation of this solution is an extremely demanding task, but crucial for the outcome of the mission. We provide a comprehensive overview of the mathematical and physical models applicable to this solution, as well as its numerical and algorithmic framework. The astrometric core solution is a simultaneous least-squares estimation of about half a billion parameters, including the astrometric parameters for some 100 million well-behaved so-called primary sources. The global nature of the solution requires an iterative approach, which can be broken down into a small number of distinct processing blocks (source, attitude, calibration and global updating) and auxiliary processes (including the frame rotator and selection of primary sources). We describe each of these processes in some detail, formulate the underlying models, from which the observation equations are derived, and outline the adopted numerical solution methods with due consideration of robustness and the structure of the resulting system of equations. Appendices provide brief introductions to some important mathematical tools (quaternions and B-splines for the attitude representation, and a modified Cholesky algorithm for positive semidefinite problems) and discuss some complications expected in the real mission data.

1. Introduction

Gaia’s astrometric core solution addresses the difficult problem of extracting a precise reference frame from enormous, entangled datasets. It estimates source, attitude, and instrument-related parameters for a carefully selected subset of well-behaved sources.

  • Mission scope: Gaia is designed to measure about one billion objects, with typical astrometric accuracies of 8–25 µas for simple stars down to 15th magnitude.The astrometric measurements are complemented by dedicated photometric and spectroscopic observations.
  • Computational challenge: The data-processing task is difficult because it combines enormous data volumes with complex relationships across data types and observation epochs.The paper situates the astrometric core solution within DPAC’s broader effort to construct the Gaia Catalogue.
  • Core unknowns: The astrometric core solution simultaneously determines source astrometry, satellite attitude, and geometric instrument calibration, with optional global parameters.These unknowns respectively represent the reference frame, instrument pointing over time, and instrument geometry; global parameters may describe effects such as deviations from General Relativity.
  • Primary sources: About 10^8 primary sources are assumed for the core solution, selected as stable, point-like, astrometrically well-behaved objects from slightly more than one billion observed sources.The selection emphasizes effectively single stars and extragalactic sources such as quasars and AGNs.
  • Scale of the solution: The problem involves roughly 5 × 10^8 unknowns, 10^11 observations, and about 70 TB of raw satellite data, ruling out a direct solution.The paper therefore presents a feasible mathematical formulation, practical solution method, and software implementation.

2. Outline of the approach

The paper formulates Gaia’s astrometric analysis as a global weighted least-squares problem linking source parameters with nuisance parameters through an observation model. Because the system is too large and entangled for direct solution, the implementation uses structured iterative processing.

  • Global formulation: The global analysis minimizes discrepancies between observed data and model-predicted data over source and nuisance parameters.The metric is defined by data statistics, yielding a weighted least-squares solution with robustness considerations.
  • Parameterization: The unknown source vector describes barycentric source motions, while nuisance parameters represent instrument and other incidental factors required to model observations.The observations may include measured detector coordinates at specific times.
  • Implementation constraints: The implementation depends on Gaia’s instrument and mission characteristics and on practical constraints imposed by data processing.Thus the abstract global formulation must be adapted to the specific observation system.
  • Model fitting: The model computes expected CCD output from source and nuisance parameters, while processing adjusts those parameters to fit observed CCD data.The paper primarily addresses the geometrical analysis and gives only a brief account of CCD-level data modelling and processing.
  • Observation variables: The minimization is implemented in field angles η and ζ rather than directly measured pixel coordinates κ and µ.This choice is part of the practical observation-equation formulation.
  • Paper roadmap: The paper organizes the treatment into mathematical modelling and least-squares equations, iterative solution methods, and detailed source, attitude, calibration, and global-processing components.Complications omitted from the main description are discussed separately in an appendix.

3. Mathematical formulation of the basic observation model

Gaia’s basic observation model combines relativistic reference systems, source-motion parameters, satellite attitude, and geometric instrument calibration. It translates source and instrument states into field-angle observations while accounting for mission-specific time, viewing, and calibration conventions.

  • 3.1. Reference systems: Gaia’s high target accuracy requires a consistent relativistic model for the observer, source, light propagation, and transformations between reference systems.The adopted formulation uses the parametrized post-Newtonian framework described by Klioner.
  • 3.1. Reference systems: The Gaia orbit and light propagation are modelled in the BCRS with TCB as the time coordinate, while satellite attitude and observed directions use the co-moving CoMRS.Ground-based ranging and on-board-clock calibration require the GCRS, but those processing aspects are outside this paper’s scope.
  • 3.1. Reference systems: In the CoMRS, Gaia’s three-dimensional attitude is represented as a spatial rotation, with quaternions providing one possible representation and the instrument-aligned frame called the SRS.The SRS is defined from the preceding and following viewing directions through the relations x = ⟨f_F + f_P⟩, z = ⟨f_F × f_P⟩, and y = z × x.
  • 3.2. Astrometric model: The astrometric model calculates each source’s proper direction at an observation time from astrometric parameters and auxiliary data including Gaia’s barycentric ephemeris and relevant solar-system ephemerides.The model uses the standard Hipparcos-derived formulation, with barycentric time corrected for the Römer delay.
  • 3.2. Astrometric model: Each source is assumed to have uniform space velocity relative to the solar-system barycentre, represented through six astrometric parameters related to position, parallax, and proper motion.The parameters include right ascension, declination, annual parallax, two transverse proper motions, and radial proper motion.
  • 3.2. Astrometric model: The reference epoch is preferably near the mission midpoint to reduce statistical correlations between position and proper-motion parameters.Transforming between kinematic and astrometric parameters is non-trivial because light-propagation timing is not fully modelled.
  • 3.3. Attitude model: Attitude deviations can reach about 1 arcmin and vary over seconds to minutes, so the exposure-smoothed effective attitude must be estimated at sub-mas accuracy.The finite CCD integration time prevents direct observation of the instantaneous physical attitude.
  • 3.4. Geometric instrument model: Geometric calibration separates large-scale AL, small-scale AL, and large-scale AC effects, using models suited to their spatial structure and expected time variation.The formulation may require modification after flight data, although a more flexible generic calibration model has been implemented for testing alternatives and systematics.

4. Solving the global minimization problem

Gaia’s astrometric core solution addresses an infeasible direct least-squares problem by iteratively updating source, attitude, calibration, and global parameters. The analysis explains convergence despite a singular normal matrix and motivates practical choices among iterative schemes.

  • About 5 × 10^8 unknowns and 8 × 10^10 observations make a direct solution infeasible, motivating iterative methods.The data are highly entangled and include calibration connectivity.
  • AGIS cycles through source, attitude, calibration, and optional global updates, each using estimates from the other blocks.The blocks are iterated because each process depends on data produced by the others.
  • The normal matrix is singular because small reference-frame rotations and inertial spins leave the differential observations unchanged.The resulting solution is aligned with the ICRS afterward using a frame rotator.
  • The simple iteration converges when the spectral radius of the iteration matrix, excluding null-space eigenvalues, is below 1.The null-space component remains unchanged, while the other error components decay; frame rotation restores the required null-space components after convergence.
  • The implemented scheme lies between block Jacobi and Gauss–Seidel, while practical AGIS experience favors S[ACG] because source–attitude coupling dominates convergence.Calibration generally converges faster, so including its newest update has limited effect on the overall convergence rate.
  • Krylov methods offer more efficient approximations, but changing observation weights can destabilize conjugate gradients, whereas simple iteration is extremely stable.The conjugate-gradient algorithm assumes a constant normal matrix, which is violated when weights change with residuals.

5. Updating processes

The AGIS iteration loop updates source, attitude, calibration, and global parameters through distinct processing blocks, with robust inner iterations handling source outliers and excess noise. Attitude updates exploit local spline structure, overlapping segments, and regularization to remain computationally feasible and numerically stable.

  • The AGIS kernel process consists of source, attitude, calibration, and global updating blocks forming the iteration loop.
  • Source updating: Each source update solves for the five astrometric-parameter corrections one source at a time using weighted normal equations.The source parameters are updated with attitude, calibration, and global parameters held fixed.
  • Source updating: The source inner iteration jointly estimates excess source noise and robust observation weights, repeating the process until the updated residuals and weights stabilize.The excess-noise estimate is constrained by the expected chi-square behavior of the objective function, while downweighting factors reduce the influence of outliers.
  • Attitude updating: Overlapping attitude segments duplicate some observation processing but reconcile updates in their overlap; with one-day segments and L = 12, the duplicated fraction is 0.4%.Segment overlap allows consecutive attitude solutions to produce essentially matching updates in the middle of the overlap region.
  • Attitude updating: A small regularization parameter λ = 10^-3 to 10^-2 yields an attitude solution that is numerically stable and insensitive to the precise value of λ.Regularization is needed because long attitude intervals can approach an ill-posed problem, while quaternion scaling creates a rank defect.

6. Auxiliary processes

The auxiliary processes align AGIS results with the ICRS, select stable primary sources, and quantify astrometric uncertainties and correlations. These tasks address the solution’s arbitrary reference frame, source suitability, and singular normal matrix.

  • 6.1. Frame rotator: The frame rotator transforms source positions, proper motions, and attitude parameters from AGIS’s arbitrary frame into a frame closely representing the ICRS.AGIS converges with source and attitude parameters in the same but largely arbitrary frame; the rotator applies the required common transformation.
  • 6.1. Frame rotator: The AGIS and ICRS frames differ by a time-dependent rotation modeled as a uniform spin, parameterized by orientation and spin vectors ε and ω.The six frame-rotator coordinates are evaluated at an adopted epoch tfr and support rigorous transformations between the frames.
  • 6.1.4. Determination of the frame rotator parameters: 4 µas yr−1 is the expected amplitude of the proper-motion pattern from solar-system acceleration, which Gaia should detect using sufficient primary quasars.The frame-rotator solution includes the acceleration vector components in the ICRS.
  • 6.2. Selection of primary sources: Primary-source selection uses a relegation factor U, where larger values indicate less suitable sources and threshold crossings can demote or promote sources.Potential primary sources are processed through the source update before their primary or secondary status is decided.
  • 6.3. Computation of standard uncertainties and correlations: The catalogue must provide standard uncertainties and correlations for astrometric parameters, including covariance matrices for pairs of sources derived from the pseudo-inverse of the singular normal matrix.Correlations arise from attitude-modelling errors and finite astrometric weight in attitude determination.

7. Software implementation and demonstration solutions

The AGIS software uses iterative, modular processing to solve simulated Gaia astrometry, with convergence monitored through source, attitude, calibration, and global-parameter updates. A large demonstration run shows that extended iteration reduces errors and handles severe spatial weight contrasts, though convergence can slow in dense regions.

  • AGIS software overview: The RunManager supports multiple iteration schemes, combining provisional updates from source, attitude, calibration, and global computations.The kernel uses a Gauss-Seidel-type preconditioner, while final updates can be formed through different algebraic schemes.
  • AGIS software overview: 50 million primary sources were successfully processed for a simulated 5 yr mission, marking a milestone toward the planned 100 million-source solution.This run was about a factor of two below the envisaged final primary-source count.
  • Convergence and source results: 135 conjugate-gradient iterations reduced example parallax and along-scan attitude updates to at or below 5 × 10^-4 µas.The conjugate-gradient scheme was reinitialised after 40 iterations to avoid a previously observed slow-convergence phase.
  • Convergence and source results: 146 µas was the settled typical parallax-error RSE, reached around iteration 25, while updates decreased to 10^-3 µas around iteration 120.The overall ∼150 µas error is representative of the median magnitude, G ≃19.
  • Convergence and source results: Several hundred µas of spatially correlated parallax error remained near iteration 20, but regional errors virtually disappeared by iteration 135.Early errors reached a few mas in high-density galactic-plane regions, whereas sufficiently long iteration allowed the system to cope with large weight imbalances when source density supported attitude determination.
  • Convergence and source results: Parallax-error patterns vary mainly with ecliptic latitude because the scanning law controls observation counts and scan geometry.The over-density of observations near ±45° ecliptic latitude does little to improve parallaxes because those observations mainly shift the AC direction.
  • Convergence and source results: Median update maps provide a convergence diagnostic for real mission data, where true error maps are unavailable, and iteration-135 truncation errors should be well below 1 µas.The amplitude and spatial distribution of updates indicate systematic truncation errors in the non-converged solution.
  • Convergence and source results: High source-density contrasts slow convergence because attitude errors are corrected mainly by the denser field of view, weakening counterbalancing information from the other field.Other datasets with lower density contrasts converged more rapidly.

8. Conclusions

The astrometric core solution is central to Gaia because it delivers primary-source astrometry, the instrument attitude reference frame, and geometric calibration. AGIS’s fundamental mathematical, algorithmic, and software components were implemented and validated in test runs, although further development was needed before flight data.

  • The core solution will eventually encompass at least 100 million primary sources while producing astrometric results, an attitude-based reference frame, and instrument calibration.
  • AGIS was built within DPAC’s Coordination Unit 3 to implement the astrometric core solution.
  • Test runs demonstrated the theoretical validity and practical feasibility of the global iterative approach, including data management and computations.
  • The fundamental parts of AGIS were already in place, but additional features and complications remained before the software was ready for flight data.
  • The authors expected AGIS to compute an accurate astrometric solution consistent with Gaia’s goals.

A.2. Spatial rotations

Unit quaternions provide a compact representation of spatial orientations and rotations. The appendix distinguishes vector and frame rotations, explains their composition rules, and relates time-dependent quaternions to angular velocity.

  • Unit quaternions represent three-dimensional orientations and spins compactly, with numerical and computational advantages over rotation matrices and Euler angles.
  • A rotation is represented by a unit quaternion, and successive rotations are composed through quaternion multiplication.
  • Quaternion sign ambiguity requires continuity handling for continuously changing orientations such as Gaia’s attitude.
  • Vector rotation: Vector rotations require quaternion multiplication and its inverse in a chosen reference system, while the resulting physical vector is reference-system independent.
  • Frame rotation: Frame rotations transform coordinates of the same physical vector when the reference system itself is rotated.
  • For a time-dependent attitude quaternion, the angular velocity can be calculated in either celestial or instrument coordinates from the quaternion and its time derivative.

A.5. The attitude matrix

The attitude matrix and quaternion representations encode the rotation from Gaia’s celestial reference system to its instrument system. A stable conversion algorithm is provided, with sign handling needed to maintain temporal continuity.

  • The attitude quaternion represents the rotation from the celestial reference system CoMRS to the instrument system SRS.
  • The attitude matrix is expressed through quaternion components, providing the coordinate transformation associated with the attitude.
  • Converting an attitude matrix to a quaternion is numerically less straightforward, so the appendix gives Klumpp’s stable algorithm in pseudocode.
  • The conversion algorithm returns a quaternion with qw ≥ 0, but sign reversal may be required to preserve temporal continuity.

A.6. Differential rotation

Differential rotation expresses the small orientation difference between two nearly co-aligned reference systems through three rotation angles. The appendix derives the first-order quaternion transformation used for attitude error angles and summarizes spline representation for time-dependent attitude modelling.

  • Differential rotation: The differential-rotation approximation applies only when the rotation-angle norm is much smaller than one.
  • Differential rotation: Three small angles φx, φy, and φz describe the spatial rotation bringing one nearly co-aligned reference system into coincidence with another.
  • Differential rotation: If q0 and q1 represent the two systems, the relative frame rotation d satisfies q1 = q0d.
  • Differential rotation: Equations A.21–A.22 transform the quaternion pair into the three differential-rotation angles, with the sign of dw ensuring correct angle signs.
  • Spline representation: The attitude representation uses splines as piecewise polynomials joined at knots, typically with cubic order M = 4.

B.1. B-splines

B-splines represent spline functions as linear combinations of minimally supported basis functions defined by a knot sequence. Their local support yields sparse, banded least-squares systems and efficient evaluation.

  • B-splines are basis functions with minimal support, uniquely defined by the adopted knot sequence for a given order and smoothness.
  • For cubic splines, the non-zero parts extend across four knot intervals, while the spline remains continuous through its second derivative at ordinary knots.
  • For order M, each B-spline spans at most M consecutive knot intervals, making the least-squares normal matrix symmetric and banded with half-bandwidth M −1.Cholesky decomposition preserves this band structure and is therefore efficient for solving the normals.
  • In attitude determination, each of the four quaternion components uses a spline, producing 4N parameters and a 4N × 4N normal matrix represented in 4 × 4 blocks.
  • At any time t, at most M B-splines are non-zero, and de Boor’s algorithm computes them simultaneously in a numerically stable way.The same algorithm can also compute derivatives with little additional effort.

B.3. Use of multiple knots

Multiple knots control endpoint representation and derivative continuity in B-splines. Their placement can preserve the fitted spline while materially affecting the numerical conditioning of the least-squares system.

  • A knot of multiplicity m removes continuity for derivatives of order M −m and higher; ordinary multiplicity-one knots preserve continuity through derivative order M −2.
  • The spline interval uses τM−1 = tbeg and τN = tend so that exactly M B-splines are non-zero throughout the interval.
  • The spline interval [tbeg, tend] is the interval over which the spline is fitted to the data, while interior knots determine its flexibility, including through multiplicities where required.
  • The first and last M−1 knots may be placed in different ways, including M-fold endpoint knots, without changing the spline between tbeg and tend.
  • Collapsing anterior and posterior knots onto the interval endpoints produces a much smaller condition number than a regular arrangement.

C.1. The use of normal equations

The paper solves least-squares systems through normal equations and Cholesky factorization, exploiting the sparse band structure produced by spline fitting. This choice prioritizes computational efficiency for very large, well-posed problems.

  • The normal-equation system Nx = b is solved by Cholesky decomposition of the symmetric normal matrix N.The factorization and triangular solves can be performed in place.
  • Cholesky is well suited to the sparse band matrix arising in B-spline fitting because it preserves the band structure.
  • Normal equations are generally less stable than methods operating directly on Ax ≃ h, but the authors use them for very large problems regarded as inherently well-posed.
  • The three computational steps are N = U′U, followed by solving U′y = b and Ux = y.
  • Selected elements or submatrices of N−1 can be obtained without computing the full inverse.

C.3. Application to semidefinite systems

A modified Cholesky procedure extends factorization and solution to positive semidefinite normal matrices without pivoting. It supports efficient banded computations but provides only an unreliable rank-defect estimate.

  • Positive semidefinite normal matrices can arise from missing observations or deliberately unconstrained calibration parameters, so standard Cholesky may fail despite existing solutions.
  • The modified algorithm continues computation for rank-deficient systems by setting singular directions to zero and producing a valid, generally non-unique solution.
  • The method avoids pivoting, preserving the matrix envelope and making it especially suitable for banded and envelope-based sparse matrices.
  • Limited experiments found the algorithm useful, simple, and efficient as a substitute for more sophisticated semidefinite-system methods.
  • The algorithm estimates rank defect, but that estimate is described as rather unreliable; in the positive definite case it is equivalent to standard Cholesky.
  • The baseline source, attitude, and instrument models are only first-order approximations, with additional effects requiring future modelling for microarcsecond-level results.

D.2. Charge Transfer Inefficiency of the CCDs

Radiation damage causes charge-transfer inefficiency that makes CCD samples depend on prior illumination, complicating the astrometric signal model. Gaia mitigates these effects through a calibrated charge distortion model and periodic charge injection.

  • Detector effects: Radiation-induced traps capture and delayed-release charges during TDI, producing charge-transfer inefficiency and related detector effects.The effects are especially relevant because the CCD signal model assumes a perfectly linear detector.
  • Detector effects: Each damaged sample depends on the undamaged value of that sample and preceding samples because the CCD state reflects illumination history.The CCD state is represented by Ψ, including the accumulated radiation damage.
  • Mitigation: Periodic electronic charge injection fills harmful traps temporarily and resets pixel illumination history, reducing subsequent charge-transfer inefficiency.The injections occur in a few consecutive TDI lines, with an example frequency of once per second.
  • Calibration: The charge distortion model uses parameters estimated alongside LSF or PSF calibration so corrected image locations can support a geometric instrument model.The intended calibrated locations should be achromatic and free of charge-transfer effects.
  • Calibration: Residual centroid shifts primarily depend on source magnitude, time since charge injection, and accumulated CCD radiation dose.Imperfections in charge-distortion-model calibration are expected to appear as systematic shifts governed mainly by these known quantities.

D.3. Effects of the finite CCD integration time

Finite CCD integration means each astrometric measurement averages the attitude and source position over the exposure rather than sampling an instantaneous crossing. The resulting effective-attitude model must account for speed mismatch, exposure moments, gating, and short-timescale attitude irregularities.

  • Exposure averaging: Finite integration makes each CCD readout depend on the average attitude and source position during the preceding exposure time T.This replaces the idealization of an instantaneous measurement at the observation line.
  • Exposure averaging: The exposure function e(τ) describes the electron-production rate over lookback time, and its normalized moments are referenced to the exposure mid-time T/2.These moments enter the finite-integration treatment of the observed image location.
  • Exposure averaging: When image and charge speeds match, the transported image location remains approximately constant throughout the integration and equals the crossing coordinate κc.The matching condition is s˙η∗ + ˙k = 0.
  • Calibration: Speed mismatch and non-constant scan rate create an astrometric correction represented through exposure moments and calibration parameters such as e1.The parameter e1 may depend on CCD and gate and may evolve with radiation damage or source magnitude.
  • Attitude representation: The effective attitude is the physical attitude convolved with the exposure function for maximum integration time T.It is used because most observations are ungated, while gated observations are corrected to refer to the same effective attitude.
  • Attitude representation: Higher-order corrections are unlikely to be profitably included when short-timescale attitude irregularities complicate computation of the physical attitude derivatives.To first order, the needed derivative can be obtained from the effective attitude.
  • Attitude representation: High-frequency attitude jitter is largely removed by TDI exposure averaging, whereas impacts can produce detectable angular-velocity changes identified through observation residuals.The basic spline attitude model may not represent all short-timescale physical irregularities.

Appendix E: Tables of acronyms and variables

Appendix E provides lookup tables for the paper’s acronyms and important variables. The tables connect notation and abbreviations to their descriptions and locations in the text.

  • Variables: Table E.2 lists important variables with short descriptions and references to where they are introduced or explained.The entries include attitude matrices, quaternion coefficients, B-spline functions, weights, and source directions.
  • Acronyms: Table E.1 lists the acronyms used throughout the paper.The table includes terms such as AGIS, CCD, CDM, LSF, and PSF.
Loading 1112.4139v1…