Nouvelle génération d’optique adaptative pour l’astronomie. Problème Inverse à grands nombres de degrés de liberté
New generation adaptive optics for astronomy. An inverse problem with large number of degrees of freedom
Observations of fainter and fainter astrophysical objects in the universe require an increasing flux-collecting area, which means larger and larger telescopes. However, the optical turbulence in the atmosphere damages the resolution of the images provided by the ground-based telescopes. In order to combine high sensitivity and angular resolution large ground-based telescopes are nowadays associated with adaptive optics systems (OA). OA must provide real-time correction of the atmospherical perturbation thanks to a servo-loop system including at least one wavefront sensor (AFO), a controller and one deformable mirror (MD) to apply the compensation, as illustrated by figure 1.
The OA design must evolve for the future telescopes with mirror diameters of several tenths of meters planned for 2017. On the one hand, the sensed part of the wavefront (WF) will not always match the part to be corrected for the astrophysical observations, so that the criterion of performance should be expressed in the WF space. On the other hand, the next generation of OA will have to cope with several thousands of actuators, to be controlled in real-time, which requires new and fast algorithms.
The outline of this paper can be decomposed in three parts. A first section introduces the statistics and the particularities of the studied signals in OA. Next, the optimization of a criterion related to astrophysical images quality leads to a command expression including a maximum a posteriori reconstructor  and an internal model control .
A fast iterative algorithm, FrIM , involving a fractal preconditioning is then used to compute in real-time such a control.
We consider here a classical OA system as illustrated by figure 1. From Kolmogorov’s model of the turbulence in the atmosphere, the kinetic energy is transferred from large scale vortices towards the smallest structures. This leads to stochastic spatial and temporal variations of the refractive index of the air, distorting the incident WF and producing local phase delays on the pupil. The statistics of the WF in the pupil plane can be described by the structure function in equation (1), where <.>ρ’represents the expectation over coordinate vector ρ’, and r0 is Fried’s parameter which characterizes the turbulence strength . The statistics of the WF perturbations w have fractal characteristics since they are invariant by a change of scale.
The WFs are sampled on a regular two-dimensional grid larger than the pupil, and are denoted by vectors w and wc in RN, for the incident WF and its compensation respectively. The AFO provides a measurement vector of local spatial derivatives, in two dimensions, d in Rm and the command vector a to be send to the actuators. The measurement equation (2) involves a linear model of the AFO, S, and a vector of measurement noise and modeling error, e. Indices k reveal the temporal discretization and enhance the characteristic delays of the OA system; here the exposure time affects the data. Errors e are supposed to have zero-mean Gaussian statistics, hence a diagonal covariance matrix Ce for the noise. The MD is also modeled by a linear device M (cf. equation (3)).
The aim of the OA is to control a multi-inputs multi-outputs system with several thousands of degrees of freedom N at a frequency close to 1 kHz. The perturbation w is not directly available via the data, only its local spatial derivatives are provided, therefore the computational time required to compute the command increases with N.
Image quality is assessed for the astrophysical observations thanks to the Strehl ratio, which is optimal for the minimum variance of the residual WF distortions in equation (5). P is a linear operator masking the telescope pupil and substracting the WF component of uniform delay over the pupil. This component, called the piston mode, has no impact on the criterion <ε2>. We consider a control law written as a linear combination of the residual measurements and the previous command (cf. equation (8)). This expression takes into account another delay, which is the time required to read the detector, obtain the measurements and compute the command ak. This new delay is considered to be equal to the first one. The control law is expressed thanks to a recursive formula (8), so that the criterion in (7) can be developed as a function of all the past measurements and the first command. Optimizing such a developed equation with respect to R and Q would provide an optimal, although complex to determine, solution. Another approach is adopted here, which consists in optimizing the convex criterion (7) with respect to both R and Q, assuming the previous command ak-1 is well-known. Choosing Q = R · G and R defined by equation (15) optimizes the criterion, so that the optimal command is written in equation (17). This OA control in discrete-time is represented by the diagram on figure 2. The obtained architecture reveals an internal model control . The computation procedure is composed of 3 steps thanks to operators E, Γ and F in equations (18)-(21). They respectively stand for a maximum a posteriori reconstructor , a predictor based on the spatio-temporal covariance Tw and a projector in the mirror actuators space. For N = 104 degrees of freedom and about two times more measurements, if the global matrix R in equation (17) is precomputed, its application to a vector of pseudo open-loop measurements dk-1 at the frequency of 1 kHz would require = 4 × N2 × 103 ≈ 4 × 1011 operations per second. It can be noticed on figure 3 that, although such multiplication (solid line) is implemented on the existing OA systems, the predictions of processors capabilities for the next decade enhance the requirement of faster algorithms.
The details of the computation are studied for the most frequent design, which makes matrices S, Ce, M and F sparse and their application to a vector requires O(N) operations. The covariance matrices of the turbulent WF are not sparse, neither is the inverse C-1w. However, it is important for the new generation of OA to take into account those priors. This is why several new methods [9, 10, 2] start the computation by iteratively solving the linear system (22)-(23), using a preconditioned conjugate gradients (GCP) algorithm.
In order to apply A at every iteration, the Fourier domain GCP method  approximates the regularization C-1w in the Fourier domain with a diagonal matrix. The global cost of this method, in figure 3, scales as O(N ln N). We consider here the faster Fractal Iterative Method, FrIM , which uses a fractal linear change of variable K in equation (25), such that approximation (24) is valid. The method builds K by a mid-point algorithm, according to Kolmogorov structure function in equation (1), and as a one-to-one linear relation w = K · u, where u in RN are independent Gaussian variables of zero-mean and variance equal to unity. The essential asset is that K, KT, and their inverses can be applied to a vector of RN in about 6 × N operations. FrIM benefits from this low cost computation to precondition the GCP, so that the linear system is effectively solved in u variables, using equations (26)-(27), in O(N) operations. The convergence of the GCP is efficiently accelerated by this change of variable and by an additional diagonal
preconditioner. The estimation of wk from the data involves less than 34 N operations per iteration.
It must also be noticed that the coherence time of the turbulent WF is assumed to be greater than the loop period, so that the approximation Γ = Id, the identity matrix, is applied.
FrIM has already demonstrated to be the fastest WF reconstruction method for large N . The turbulence dynamics leads to even less computations when the method is involved in a closed-loop OA. The number of iterations required to converge depends on both the signal to noise ratio  and the evolution speed of the turbulence. The iterative method in closed-loop OA starts with the best estimate uk-1 of the previous loop. The simulation of 3 extreme cases of temporal evolution of the turbulence illustrates this in figure 4, bounding the values of niter required for the convergence of the estimation. The reconstruction mean square error in equation (28) is plotted along the 15 iterations of 6 successive OA loops. It appears that when the temporal correlation of the WF increases, which also means when Γ gets closer to the identity matrix, then the estimation convergence is reached in less iterations. The solid curve is a realistic situation, obtained simulating a translating WF above the pupil. In such case, initializing the iterative algorithm with the previous final estimate leads to a convergence reached in 2 or 3 iterations, and a global cost for the command eventually smaller than the one of an isolated WF reconstruction.
As a conclusion, the optimization of the closed-loop AO correction in the WF space led to a internal model control law involving a turbulent WF reconstruction. The computation of the command for large N is then possible in real-time thanks to the use of FrIM method for the reconstruction step. The closed-loop configuration accelerates the convergence of this iterative algorithm thanks to the temporal coherence of the turbulence. This method decreases the computational cost of the command of a factor of about 100 in comparison with the classical matrix-vector multiplication approach, on a system with N = 104.
Les images astrophysiques issues des télescopes au sol sont dégradées par la turbulence de l’atmosphère terrestre. Un système d’Optique Adaptative (OA) doit corriger en temps réel ces perturbations. L’étude d’une nouvelle génération de télescopes, de plus de 30 mètres de diamètre, avec de nouveaux concepts d’OA ayant de 104 à 105 degrés de liberté, appelle de nouveaux algorithmes pour la commande du système. L’optimisation d’un critère nous a conduit à une commande en boucle fermée par modèle interne. Cette structure permet d’exploiter les a priori sur la turbulence en utilisant l’algorithme FrIM, particulièrement adapté à ces systèmes à grands nombres de degrés de liberté.
Adaptive optics, inverse problem, internal model control, conjugate gradients, preconditioner, fast algorithm
Optique adaptative, problème inverse, commande par modèle interne, gradients conjugués, préconditionnement, algorithme rapide
 F. RODDIER, Adaptive Optics in Astronomy. Cambridge University Press., 1999.
 É. THIÉBAUT and M. TALLON, Ultra fast maximum a posteriori wave front reconstruction for extremely large telescopes. submitted to J. Opt. Soc. Am. A, 2008.
 A. TARANTOLA, Inverse Problem Theory and Methods for Model Parameter Estimation SIAM, 2005.
 M. MORARI and E. ZAFIRIOU, Robust Process Control. Prentice Hall, Englewood Cliffs, New Jersey, 1989.
 V.I. TATARSKI, Wave propagation in a turbulent medium. Dover Publications, Inc. New York, 1961.
 D.L. FRIED, Statistics of a geometric representation of wavefront distortion. J. Opt. Soc. Am. A, 55 :1427-1435, 1965.
 K. B. PETERSEN and M. S. PEDERSEN, The Matrix Cookbook, Version February 10, 2007. Technical University of Denmark, 2007.
 E. P. WALLNER, Optimal wave-front correction using slope measurements J. Opt. Soc. Am., 73:1771-1776, 1983.
 L. GILLES, C.R. VOGEL and B.L. ELLERBROEK, Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction. J. Opt. Soc. Am. A, 19 :1817-1822, 2002.
 Q. YANG, C.R. VOGEL and B.L. ELLERBROEK, Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography. Applied Optics, 45 :5281-5293, 2006.