OSA's Digital Library

Optics Express

Optics Express

  • Editor: Michael Duncan
  • Vol. 14, Iss. 17 — Aug. 21, 2006
  • pp: 7487–7498
« Show journal navigation

Fast optimal wavefront reconstruction for multi-conjugate adaptive optics using the Fourier domain preconditioned conjugate gradient algorithm

Curtis R. Vogel and Qiang Yang  »View Author Affiliations


Optics Express, Vol. 14, Issue 17, pp. 7487-7498 (2006)
http://dx.doi.org/10.1364/OE.14.007487


View Full Text Article

Acrobat PDF (178 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

We present two different implementations of the Fourier domain preconditioned conjugate gradient algorithm (FD-PCG) to efficiently solve the large structured linear systems that arise in optimal volume turbulence estimation, or tomography, for multi-conjugate adaptive optics (MCAO). We describe how to deal with several critical technical issues, including the cone coordinate transformation problem and sensor subaperture grid spacing. We also extend the FD-PCG approach to handle the deformable mirror fitting problem for MCAO.

© 2006 Optical Society of America

1. Introduction

Multi-conjugate adaptive optics (MCAO) [1

1. J. M. Beckers, “Increasing the size of the isoplanatic patch with multi-conjugate adaptive optics,” in Proceedings of European Southern Observatory Conference and Workshop on Very Large Telescopes and Their Instrumentation, M. H. Ulrich, ed., Vol. 30 of ESO Conference and Workshop Proceedings (European Southern Observatory, Garching, Germany, 1988), pp. 693–703.

, 2

2. D. C. Johnston and B. M. Welsh, “Analysis of multi-conjugate adaptive optics,” J. Opt. Soc. Am. A 11, 394–408 (1994). [CrossRef]

] refers to the use of several deformable mirrors (DMs), each conjugate to a different altitude, to correct for volume atmospheric turbulence and thereby increase the corrected field of view in an astronomical adaptive optics system. By wavefront reconstruction we mean the determination of DM actuator commands, given wavefront sensor signals. In this paper we assume a linear model for sensor signals as a function of the turbulence profile, we assume DM displacements depend linearly on actuator commands, and we employ an optimal, or minimum variance, approach to wavefront reconstruction [3

3. T. Fusco, J. M. Conan, G. Rousset, L. M. Mugnier, and V. Michau, “Optimal wave-front reconstruction strategies for multi-conjugate adaptive optics,” J. Opt. Soc. Am. A 18, 2527–2538 (2001). [CrossRef]

]. This gives rise to a pair of linear systems—one to estimate the volume turbulence profile and second to determine the actuator commands. Solution of the first system is referred to as the estimation step, or tomography, and solution of the second is called the fitting step.

With planning for a future 30-meter class telescope (TMT) that will rely on MCAO underway [4

4. R. G. Dekany, M. C. Britton, D. T. Gavel, B. L. Ellerbroek, G. Herriot, C. E. Max, and J-P. Veran, “Adaptive optics requirements definition for TMT,” Advancements in Adaptive Optics, edited by D. B. Calia, B. L. Ellerbroek, and R. Ragazzoni, Proc. SPIE 5490, 879–890 (2004). [CrossRef]

], it is important to be able to solve both the estimation step and the fitting step very quickly. Since TMT will require tens of thousands of sensors and actuators, conventional wavefront reconstructors, which are based on direct matrix-vector multiplication, may be prohibitively expensive. This has motivated the development of alternative approaches that rely on sparse direct matrix equation solvers [5

5. B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A , 19, 1803–1816 (2002). [CrossRef]

] and on iterative methods like the preconditioned conjugate gradient algorithm (PCG) [6

6. G. Golub and C. VanLoan, Matrix Computations, 2nd Edition, Johns Hopkins University Press, 1989.

] with multigrid preconditioning (MG-PCG) [7

7. 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). [CrossRef]

, 8

8. L. Gilles, B. L. Ellerbroek, and C. R. Vogel, “Preconditioned conjugate gradient wave-front reconstructors for multi-conjugate adaptive optics,” Appl. Opt. 42, 5233–5250 (2003). [CrossRef] [PubMed]

, 9

9. B. L. Ellerbroek, L. Gilles, and C. R. Vogel, “Numerical simulations of multi-conjugate adaptive optics wavefront reconstruction on giant telescopes,” Appl. Opt. 42 (2003), pp. 4811–4818. [CrossRef] [PubMed]

] and Fourier domain preconditioning (FD-PCG) [10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

].

It should be noted that these iterative methods yield approximations to the optimal wavefront reconstruction. In simulations [9

9. B. L. Ellerbroek, L. Gilles, and C. R. Vogel, “Numerical simulations of multi-conjugate adaptive optics wavefront reconstruction on giant telescopes,” Appl. Opt. 42 (2003), pp. 4811–4818. [CrossRef] [PubMed]

, 10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

] both MG-PCG and FD-PCG have been shown to be very rapidly convergent, so few iterations are needed. In addition, the cost per iteration is low, so both algorithms provide essentially optimal reconstructions at relatively low cost.

This paper is organized as follows. In section 2 we review basic MCAO concepts and introduce notation. Section 3 contains a review of optimal volume turbulence estimation. In subsection 3.1 we present the cone-coordinate transformation problem. In subsections 4.1 and 4.2 we review the basic ideas underlying the FD-PCG algorithm and we present 2 different implementations of the algorithm. In subsection 4.3 we describe how to deal with computational grids that are finer than the sensor subaperture grid. Section 5 contains a review of the fitting step and an outline of how to efficiently apply FD-PCG to fitting. In subsection 5.1 we demonstrate FD-PCG fitting for a simulated MCAO system for a 30-meter class telescope. We close with a summary in section 6.

2. Basic Concepts and Notation

We assume a layered atmosphere with a refractive index profile which we denote by ψ(x,zl ), l = 1,…,n layer, where the discrete layer heights zl are known, and x = (x,y) denotes lateral position. Since rapid refractive index variations are concentrated in regions of high turbulence, ψ is often referred to as the atmospheric volume turbulence profile. (Note that one can obtain the Cn2 profile, a standard measure of turbulence strength as a function of height z, by integrating the square of ψ.) We further assume that at any instant in time, ψ is a realization of a wide-sense stationary stochastic process characterized by a Kolmogorov power spectral density function. See [11

11. J. W. Hardy, Adaptive Optics for Astronomical Telescopes, Oxford University Press, 1998.

] for details.

The phase, or wavefront aberration, is obtained by propagating light from an idealized point light source, either at infinity in the case of a natural guidestar (NGS) or at a known finite altitude in the case of a laser guidestar (LGS), through the atmosphere to the pupil plane of the telescope. We denote dependence of phase on volume turbulence by

ϕ(x,θ)=[]xθ,

where θ denotes direction of the guidestar and x denotes location in the pupil plane. Using a geometric optics approximation, we model the propagator P as

[]xθ==1nlayerψ(x+zθ,z).
(1)

In our model, both ψ and ϕ represent optical path differences and have units of length.

We assume that wavefront sensor measurements can be modelled as noisy, discrete approximations to the gradient of the phase. For example, suppressing the dependence on direction θ, the output of the (i,j)th element of a Shack-Hartmann wavefront sensor can be modelled [12

12. D. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. , Vol. 67 No. 3, (1977), pp.370–375 [CrossRef]

] as

si,jx=1h[(ϕ(xi+h,yj)ϕ(xi,yj))2+(ϕ(xi+h,yj+h)ϕ(xi,yj+h))2]+ηi,jx,
si,jy=1h[(ϕ(xi,yj+h)ϕ(xi,yj))2+(ϕ(xi+h,yj+h)ϕ(xi+h,yj))2]+ηi,jy,
(2)

where h represents the sensor subaperture width and ηi,jx and ηi,jy represent sensor noise and grid discretization effects. Thus si,jxϕx(xi,yj) and si,jyϕy(xi,yj). We use the symbol Γ to represent the (discrete gradient) mapping from phase ϕ to sensor measurements s = (sx , sy ).

Discretization in the direction variable ψ arises naturally because of the finite number of guidestars in a practical MCAO system. We impose a nodal discretization on ψ in the lateral variable x at each layer zl and nodal discretization of ϕ in each (discrete) direction θ. Hence the operators P and Γ have discrete matrix (approximate) representations. In addition, the Kolmogorov spectrum for (discretized) volume turbulence has a block diagonal matrix representation in the Fourier domain, with one diagonal block for each layer, due to independence of the layers.

In order to efficiently implement computational techniques based on the discrete Fourier transform, it is necessary to impose equispaced rectangular grids on the lateral variable x = (x,y) at each altitude zl . However, telescope apertures are typically circular or annular. Hence it is necessary to introduce an additional masking operator M. We take [Ms](xj ,yj ) to equal s(xi ,yi ) if (xi ,yj ) lies within the aperture and we take it to equal zero for (xi ,yi ) outside the telescope aperture. In an abuse of notation, we will also refer to the masking function whose value M(xj ,yj ) = 1 for (xj ,yj ) inside the aperture and M(xj ,yj ) = 0 for (xj ,yj ) outside the aperture.

The composition of discretized propagation, discrete gradient (corresponding to sensing), and masking will be represented by the volume turbulence-to-sensor measurement operator

G=MΓP.
(3)

By atmospheric turbulence tomography, or volume turbulence estimation, we mean the task of estimating the (discretized) volume turbulence profile ψ from noisy sensor measurements

s=+η.

In the implementation of the fitting step, one first propagates the light to the pupil plane from guidestars at infinity in several specified sample directions to obtain an estimated phase, ϕ est. The sampling must be dense enough to well-represent the phase, but not so dense as to make the computations overly expensive. The estimation plus this first stage of the fitting step can together be viewed as a glorified interpolation scheme to extend the phase to directions that have not been sampled by the guidestars.

ϕDM=Pm.
(4)

Here P is again a geometric optics propagation operator with a representation as in Eq. (1); the DM displacements are analogous to turbulence layers placed at the DM conjugate altitudes. We assume the DM displacements can be represented as

m(x,zk)=[Hkak](x),k=1,,nDM.
(5)

The vector a k represents the actuator commands for the kth DM, and Hk is the actuator-to-DM influence operator. We represent the combination of Eqs. (4) and (5) as

ϕDM=PHa,
(6)

where a = (a 1,…,a nDM) is the concatenation of the actuator commands. As in the estimation step, nodal discretization of the DM displacements gives rise to a discrete matrix to represent H. The second stage of fitting can then be formulated as follows: Find the actuator command vector a for which H a best matches (in a sense that we will precisely define in section 5) the estimated phase ϕ est.

3. Turbulence Estimation for an LGS-NGS MCAO System

To obtain sufficiently bright light sources to perform accurate volume turbulence profile estimation, it is necessary to make use of LGS beacons. Unfortunately, these beacons provide limited information about certain low-order components of the profile like tip-tilt, so additional NGS beacons are needed [5

5. B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A , 19, 1803–1816 (2002). [CrossRef]

]. In this section we address technical issues that arise from a mixture of LGS and NGS beacons.

s=[shst]=[GhψGtψ]+[ηhηt].
(7)

Each of Gh , Gt has form similar to G in Eq. (3). From [3

3. T. Fusco, J. M. Conan, G. Rousset, L. M. Mugnier, and V. Michau, “Optimal wave-front reconstruction strategies for multi-conjugate adaptive optics,” J. Opt. Soc. Am. A 18, 2527–2538 (2001). [CrossRef]

], the optimal volume turbulence estimate is then given by

ψest=argminψ{GhψshCh12+GtψstCt12+ψCψ12}
(8)

where xA2=defxTAx,Ch is the covariance matrix for high order sensor noise,Ct is the covariance matrix for low order sensor noise, and Cψ is the covariance matrix for volume turbulence. The first two terms in (8) correspond to the pair of entries in (7), while the third is a prior, or regularization term, whose role is to incorporate prior statistical information about the volume turbulence.

By setting the gradient of the bracketed quantity in Eq. (8) equal to zero, we obtain

ψest=(GhTCh1Gh+GtTCt1Gt+Cψ1)1(GhTCh1sh+GtTCt1st)
(9)
=[CψGhTCψGtT][GhCψGhT+ChGhCψGtTGtCψGhTGtCψGhT+Ct]1[shst]
(10)

Following [5

5. B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A , 19, 1803–1816 (2002). [CrossRef]

, 10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

], we express the covariance matrix for high-order noise in the presence of tip-tilt uncertainty as

Ch1=Nh1Nh1T(TTNh1T)1TTNh1
(11)

where Nh describes the noise within the high-order sensors, and T is a low-rank block matrix with 2n LGS columns (n LGS denotes number of LGSs). The product Ch1 sh can be interpreted as a noise-weighted, tip-tilt-removed version of the high-order sensor signals sh. We then rewrite Eq. (29) as

ψest=(AhAlr)1(bh+bt)
(12)

where

Ah=GhTNh1Gh+Cψ1,
(13)
Alr=GhTNh1T(TTNh1T)1TTNh1GhGtTCt1Gt,
(14)
bh=GhTCh1sh,bt=GtTCt1st.
(15)

The two matrix products on the right hand side of Eq. (14) have rank 2n LGS and nt , respectively, where nt is the total number of low-order sensor measurements. Hence we can decompose

Alr=U1U1TU2U2T,
(16)

where U 1 has 2n LGS columns and U 2 has nt columns. On the other hand, the matrix Ah in Eq. (13) has full rank. The Sherman-Morrison formula [6

6. G. Golub and C. VanLoan, Matrix Computations, 2nd Edition, Johns Hopkins University Press, 1989.

] allows to write

(AhAlr)1=[Ah(U1U2)(U1TU2T)]1
=Ah1+[W1W2][IW1TU1W1TU2W2TU1I+W2TU2]1[W1TW2T]L
(17)

where

W1=Ah1U1,W2=Ah1U2.
(18)

Then setting b = bh + b lr, we obtain from (12) and (17)

ψest=Ah1b+Lb.
(19)

From (17) we see that the product Lb can be computed by (i) taking a few dot products to obtain W1T b and -W 2 T b; (ii) applying the inverse of a small (2n LGS + nt × 2n LGS + nt ) matrix to the result of (i); and (iii) adding together scalar multiples of the columns of W 1, W 2, where the scalars come from (ii). The product Ah1 b is computed iteratively using FD-PCG.

3.1. The Cone-Coordinate Transformation

With geometric optics propagation, the contribution to the pupil-plane phase at position x due to light that has passed through a turbulence layer at altitude zl from a guidestar at height H with orientation θ is proportional to

ψ(cx+s,z)=defψ˜(x˜+s˜).
(20)

Here c = 1 - zl /H and s = zlθ are the cone compression factor and the shift, respectively, and the mapping

xcx=defx˜.
(21)

is referred to as the cone coordinate transformation.

From Eq. (20) we see that the cone coordinate transformation allows us to replace the combination of a rescaling and a shift with a simple shift. This simplifies the representation of the high-order geometric optics propagator Ph and facilitates the FD-PCG method. Figyre 1 illustrates the discrete grids that conform to cone coordinates and to standard coordinates.

Fig. 1. Illustration of cone-coordinate grid vs conventional equispaced grid. Cone-coordinate grid points in figure on the left have grid spacings that decrease with layer height. The conventional grid on the right has the same number of grid points at each layer height, but the grid spacing does not vary with layer height.

The difficulty that arises with a (nodal) cone coordinate representation of the volume turbulence profile is that the low-order propagator Pt then has a more complicated representation. To overcome this problem, we interpolate back to standard coordinates before applying Pt . In cone coordinates, Eq. (8) takes the form

ψ˜est=argminψ˜{G˜hψ˜shCh12+GtIcψ˜stCt12+ψ˜Cψ˜12}
(22)

The h in the first term on the right hand side denotes the cone-coordinate representation of Gh in Eqn. (7). The Ic in the second term denotes the transformation from cone coordinates back to standard coordinates. In practice all operators have discrete matrix representations and Ic is a 2-D interpolation matrix, which is very sparse. In cone coordinates, Eqs. (12)-(15) must be slightly modified, with h taking the place of Gh and GtIc taking the place of Gt . Given ψ˜est, one computes

ψest=Icψ˜est
(23)

to get an estimate for the volume turbulence profile in standard coordinates.

4. FD-PCG Implementation of the Estimation Step

The effectiveness of the FD-PCG algorithm is a consequence of the fact that the components of the matrix Ah in Eq. (13) nearly all have nice Fourier domain representations. In particular, the discrete gradient operator Γ corresponding to the Fried model (2) has a diagonal Fourier representation due to the Fourier shift theorem. The Fourier representer for the inverse Kolmogorov turbulence covariance matrix Cψ1 is block diagonal with diagonal blocks. In cone coordinates, the high-order propagator Ph has a Fourier representer which has a block decomposition with diagonal blocks, again as a consequence of the shift theorem. The mask M unfortunately does not have a compact representation, but it can be well-approximated by a scalar multiple of a diagonal matrix, as we demonstrated in [10].

We decompose the propagator Ph into blocks Pkj , where index k corresponds to LGS direction and index j corresponds to turbulence layer. Then from Eqs. (13) and (3), Ah also has a block decomposition,

[Ah]ij=k=1nLGSPkiTSkPkj+δijBj,i,j=1,,nL,
(24)

with

Sk=ΓxTMNh1MΓx+ΓyTMNh1MΓy,k=1,,nLGS,
(25)

the δij is the Kronacker delta, the regularization term Bj corresponds to the inverse covariance matrix for the jth layer of turbulence, and Γx, Γy are the x- and y-components of the discrete gradient Γ.

In cone coordinates, each component propagator Pkj is a simple shift operator. Hence it has a Fourier representation

Pkj=F1P̂kjF,
(26)

where F represents the 2-D discrete Fourier transform on an nx × nx grid and kj is a simple filter (component-wise scalar multiplication) operator. We assume the components Γx, Γy of the gradient operator have analogous Fourier-domain filter representations Γ^ x, Γ̂y. In addition, we assume the tip-tilt removed noise covariance is scalar, Nh1 = σh2 I. Then

[Ah]ij=FÂijF1withÂij=defkP̂ki*ŜkP̂kj+δijB̂j,
(27)

where j = FBjF -1, ∗ denotes complex conjugate transpose, and

Ŝk=σh2(Γ̂x*M̂Γ̂x+Γ̂y*M̂Γ̂y),M̂=defFMF1
(28)
σh2(Γ̂x2+Γ̂y2)
=defS˜k.
(29)

As in [10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

], we have approximated the Fourier transformed pupil mask by a scalar multiple of the identity to obtain (29) from (28). For simplicity of presentation we have taken the scalar to be 1 in Eq. (29). We will revisit the issue of masking in Section 4.3 below.

4.1. Direct Implementation of FD-PCG

PCG iteration to solve Ahx = b requires storage of 4 vectors, each the size of the right-hand-side b; see [6

6. G. Golub and C. VanLoan, Matrix Computations, 2nd Edition, Johns Hopkins University Press, 1989.

] for details. The dominant costs at each iteration are typically from the operator multiplications h = Ahd and the preconditioner applications z = C -1 r.

PCG iteration can be applied directly to compute x = Ah1b in Eq. (19). In this case, motivated by Eq. (27) and approximation (29), we take the preconditioner to be the matrix C having block components

[C]ij=F1ĈijFwithĈij=k=1nLGSP̂ki*S~kP̂kj+δijB̂i,
(30)

Let C denote the matrix comprised of blocks Ĉij . We demonstrated in [10] that there exists a permutation, or reordering of rows and columns, for which Ĉ is block diagonal and the diagonal block size is quite small. (Diagonal block size will be addressed in detail in Section 4.3 below.) Hence Ĉ can be inverted directly and efficiently stored. The preconditioning step z = C -1 r is then implemented as follows:

  1. 2-D Fourier transforms are applied to the blocks of r (note that blocks correspond to turbulence layers; there are n layer blocks), yielding a vector f.

    The entries of are permuted, yielding a vector .

  2. is multiplied by the block diagonal inverse of the reordered Ĉ, yielding .

    The inverse permutation is applied to , yielding z.

  3. Inverse Fourier transforms are applied to the blocks of to obtain z.

In the direct implementation of PCG, computations h = Ahd are carried out as sparse matrix-vector multiplications using the block decomposition (24). The propagators Pkj correspond to simple shift-and-adds and need not be stored, and the matrices Sk are very sparse. We employ Ellerbroek’s biharmonic approximation to the inverse turbulence covariance [5

5. B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A , 19, 1803–1816 (2002). [CrossRef]

], so the regularization matrices Bj are also very sparse. Computations z = C -1 r are dominated by the layer-wise 2-D Fourier transforms in steps (i) and (v).

4.2. Transformed Implementation of FD-PCG and Comparison with Direct Implementation

Alternatively, one can apply PCG to the transformed system Âx̂ = b, where  has blocks given in (27) and b has blocks j = F[b]j. In this case the preconditioning step = C -1 is carried out as in the direct PCG approach, but the Fourier transforms in steps (i) and (v) are omitted.

To carry out the step ĥ = Âd̂, we see from (27) that Fourier domain propagations can be carried out as component-wise product, or filter operations, involving (complex-valued) scalar multiplications. Similarly, the regularization terms i give rise to filtering operations. However, we see from Eq. (28) that since the transformed pupil mask has a full matrix representation, a Fourier transform / inverse Fourier transform pair is required for multiplication by each Ŝk .

If n layer = n LGS, then the number of Fourier transforms applied during each PCG iteration is the same for the direct implementation and the transformed implementation. The cost of propagation and regularization is very similar using both approaches. However, the transformed implementation requires n layer additional Fourier transforms to obtain from b and n layer additional inverse Fourier transforms to obtain x from . Moreover, the transformed approach requires 4 complex-valued PCG vectors, while the direct approach requires 4 real-valued vectors, which need half as much storage.

4.3. Grid Masking for High Order Sensor Subapertures

In our model for high-order sensors we assume sensor subapertures are square and have vertices (corners) which lie on the computational grid. To facilitate very accurate turbulence profile estimation, we take a computational grid that contains additional points that are not subaperture vertices. In order for the preconditioner C to closely approximate the operator Ah , we must incorporate information about this grid structure in the Fourier domain. To this end we extend the subaperture vertex grid to the entire computational domain, as shown in Fig. 2, and we define the subaperture mask

MS(x)={1,ifxisasubaperturevertex0,otherwise.
(31)

We then replace the k in Eq. (29) with

S̃k=σh2(Γ̂x*M̂SΓ̂x+Γ̂y*M̂SΓ̂y),M̂S=F1MSF,

and we construct the preconditioner C from Eq. (30). As in [10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

], we reorder the rows and columns of the Fourier representer Ĉ according to spatial frequency, so the permuted Ĉ is block diagonal with diagonal blocks that are nb × nb , where nb = n layersx)2, with Δs equal to the sensor subaperture spacing and Δx equal to the spacing of the computational grid. For example, for a 6-layer atmospheric model with the grid shown in Fig. 2, Δsx = 2 and the diagonal blocks of the permuted Ĉ are 24 × 24.

Fig. 2. Sensor subaperture grid for a simulated 4-meter telescope. Green circles represent points in the computational grid. Blue stars represent vertices of 1/2 meter × 1/2 meter square high-order sensor subapertures. The red circle represents the outer edge of the clear aperture, or pupil.

5. FD-PCG for the Fitting Step

Given a volume turbulence profile ψ est, one propagates to the pupil plane to obtain an estimated phase, ϕ est. Slightly modifying the notation introduced at the end of section 2, we model the masked DM corrections to the phase as

ϕDM=MPDMHa,
(32)

where M is a pupil mask, a is the actuator command vector, H a represents PDM figures at the conjugate altitudes, and P DM represents propagation from the nsample virtual guidestars, through “layers” that correspond to DM displacements at the conjugate altitudes, to the pupil plane. P DM has an n sample × n DM block decomposition, the actuator influence operator H is block diagonal with n DM blocks, and a can be decomposed into n DM blocks.

The optimal actuator command vector is given by

apot=argmina{MPDMHaϕestW2+aR2}
=Afit1(HTPDMTMWϕest).
(33)

where

Afit=HTPDMTMWMPDMH+R.
(34)

The purpose of the matrix W is to allow selective weighting in certain sampling directions, and the role of R is to stabilize the actuator commands. One can also accommodate fast tip-tilt mirrors or other low-order corrective elements. These give rise to low-rank terms and can be handled as in the estimation step using the Sherman-Morrison formula.

The matrix A fit takes the form of Ah in Eq. (13) provided we identify MP DM H, W, and R in (34) with Gh , Nh1, and Cψ1 in (13). In order to efficiently implement FD-PCG, components of A fit must be well-approximated by matrices with sparse Fourier domain representations. For this reason, we select W and R to have blocks with translation invariant structure (the matrix representers are then block circulant with circulant blocks, or BCCB). The propagator P DM should automatically be translation invariant, and many DM models assume that the blocks of H are translation invariant. As was the case in the estimation step, only the pupil mask does not have a sparse Fourier representer, but it can be approximated as before. We can again employ either of the two FD-PCG solutions strategies described in Sections 4.1 and 4.2.

5.1. An Illustrative Example

We present results for the fitting step applied to the simulated MCAO system for a 30-meter telescope described in [10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

]. This incorporates a 6-layer Cerro Pachon model for the atmosphere, 5 laser guidestars at 90 Km, 4 natural guidestars to deal with tip-tilt uncertainty, two DMs (one conjugate to ground and the second to 12 Km) with a corrected field of view of 30 arcseconds, and 1/2 meter sensor subaperture and actuator spacings. For the estimation step we now use 256 × 256 computational grids at each of the 6 layers, with twice the grid resolution in [10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

] (1/4 meter grid spacing on a square domain with a 64 meter width). The performance of the FD-PCG algorithm for estimation is much the same as in [10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

], so we present no additional details here.

We propagated the estimated volume turbulence profile to the pupil plane and then carried out the rest of the fitting step using the transformed implementation of the FD-PCG algorithm; see Section 4.2. The actuator influence functions, which determine the matrix H in Eq. (4), are taken to be bilinear splines with control points that coincide with the high-order sensor subaperture vertices (see Fig. 2). The matrix H has translation invariant structure and is quite sparse. FD-PCG performance is summarized in Fig. 3.

Let a k denote the approximation to a opt in Eq. (33) after k FD-PCG iterations, and let ϕkDM denote the corresponding DM correction to the phase, obtained by substituting a k for a in Eq. (32). By the residual phase error we mean ϕ true - ϕkDM, where ftrue is the true phase, obtained by propagating through the (simulated) true volume turbulence profile to the pupil plane.

Note that at most 3 iterations are required to reduce the residual phase error to its asymptotic level, where the difference between a k and a opt is negligible and produces negligible change in ϕkDM. Any further reduction in phase error would require adding more DMs or more DM actuators to reduce the fitting error or decreasing sensor noise to reduce the estimation error.

Fig. 3. FD-PCG performance for the fitting step. The blue line shows the on-axis residual phase error, and the red line shows phase error at the edge of the 30-arcsecond corrected field of view of the telescope.

In [10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

] we observed similar behavior of the RMS phase error as we varied the number of FD-PCG iterations in the estimation step. The main difference was that 6 to 8 iterations were required to achieve asymptotic error levels for the estimation step, while only 2 or 3 iterations are required to reach asymptotic error levels in the fitting step.

5.2. Cost of Estimation and Fitting for the Example

For the example above, elements are required to store the nodal discretization of ψ in the estimation step. Standard implementations of PCG require 4 copies of ψ. There is some additional overhead to store the components of Ah , but by taking advantage of special structure (sparsity and translation invariance), this can be kept relatively small. However, the W 1 and W 2 matrices in the low-rank matrix L in Eq. (17) require a great deal more storage. Each column of each Wi has N est elements, and there are 2n LGS + nt = 2×5 + 8= 18 columns in W 1, W 2 for the simulation. Hence the storage costs in the estimation step are dominated by the low-rank terms.

Again in the estimation step, the low-rank computation Lb (which is O(N est)) is much cheaper than the cost of computing Ah1 b via FD-PCG iteration. The FD-PCG costs are the cost per iteration multiplied by the total number of iterations. Dominant costs per iteration are multiplications by Ah and by the inverse preconditioner C -1. Whether we use the direct approach in section 4.1 or the transformed approach in section 4.2 we must apply forward and inverse discrete Fourier transforms to block vectors with N est entries (each block is 256 × 256 in our simulation). Using the fast Fourier transform, the asymptotic cost for this is then O(N estlogN est). There are some additional order N est costs, e.g., from the multiplication by the relatively small blocks of Ĉ -1.

While a detailed discussion of parallel implementation is beyond the scope of this paper, it should be noted that the layered structure of ψ and the block structure of the components of Ah and of Ĉ -1 provide ample opportunities to dramatically reduce computational costs in a parallel computing environment.

Storage requirements for the fitting step are significantly lower than for the estimation step. In the above simulation, there are only 2 DMs, so with a 256 × 256 nodal discretization of each DM, only Nfit=def6×2562130,000 storage elements are needed to represent the DMs. PCG again requires 4 copies, but there are no low-rank terms.

Computational costs in the fitting step are also significantly smaller than for estimation. The asymptotic costs of the FFTs at eachfitting iteration are much smaller because N fit = N est/3 (due to 6 layers in the estimation vs 2 DMs in the fitting). In addition, the 2 or 3 FD-PCG iterations required for the fitting step are much fewer than the 6 to 8 iteration needed for estimation.

6. Discussion and Conclusions

In [10

10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

] we introduced the FD-PCG algorithm to efficiently carry out the estimation, or tomography, step in minimum variance wavefront reconstruction. Advances in this paper included

  • A demonstration that the FD-PCG algorithm can be adapted to efficiently solve the fitting step in optimal MCAO wavefront reconstruction. We show that for a simulated MCAO system for a 30-meter telescope, at most 3 FD-PCG iterations are needed for the fitting step.
  • Solution to the cone coordinate transformation problem. This problem arises with mixed LGS-NGS systems and can lead to field distortion if it is not handled correctly. Our solution involves a cone-coordinate-to-standard coordinate interpolation and is efficiently implemented with sparse matrix techniques.
  • Solution of problems related to computational grids that have higher resolution than the high-order sensor subaperture grids. These problems are solved with a subaperture mask. This leads to a somewhat more complex preconditioner, but the additional computational overhead is small relative to the cost of the 2-D Fourier transforms.
  • Two separate implementations of the FD-PCG algorithm. In the direct approach, PCG vectors lie in the spatial domain and are real-valued, while with the transformed approach the PCG vectors lie in the Fourier domain and are complex valued. The direct approach requires slightly fewer Fourier transforms and slightly less storage than does the transformed approach. Perhaps the only advantage of the transformed approach is that allows easy implementation of wave optics propagators. These have a sparse Fourier domain representation but have a convolution representation in the spatial domain.

Acknowledgments

This research was supported in part by the Air Force Office of Scientific Research through grant F49620-02-1-0297 and by a grant from the Optical Sciences Company.

References and links

1.

J. M. Beckers, “Increasing the size of the isoplanatic patch with multi-conjugate adaptive optics,” in Proceedings of European Southern Observatory Conference and Workshop on Very Large Telescopes and Their Instrumentation, M. H. Ulrich, ed., Vol. 30 of ESO Conference and Workshop Proceedings (European Southern Observatory, Garching, Germany, 1988), pp. 693–703.

2.

D. C. Johnston and B. M. Welsh, “Analysis of multi-conjugate adaptive optics,” J. Opt. Soc. Am. A 11, 394–408 (1994). [CrossRef]

3.

T. Fusco, J. M. Conan, G. Rousset, L. M. Mugnier, and V. Michau, “Optimal wave-front reconstruction strategies for multi-conjugate adaptive optics,” J. Opt. Soc. Am. A 18, 2527–2538 (2001). [CrossRef]

4.

R. G. Dekany, M. C. Britton, D. T. Gavel, B. L. Ellerbroek, G. Herriot, C. E. Max, and J-P. Veran, “Adaptive optics requirements definition for TMT,” Advancements in Adaptive Optics, edited by D. B. Calia, B. L. Ellerbroek, and R. Ragazzoni, Proc. SPIE 5490, 879–890 (2004). [CrossRef]

5.

B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A , 19, 1803–1816 (2002). [CrossRef]

6.

G. Golub and C. VanLoan, Matrix Computations, 2nd Edition, Johns Hopkins University Press, 1989.

7.

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). [CrossRef]

8.

L. Gilles, B. L. Ellerbroek, and C. R. Vogel, “Preconditioned conjugate gradient wave-front reconstructors for multi-conjugate adaptive optics,” Appl. Opt. 42, 5233–5250 (2003). [CrossRef] [PubMed]

9.

B. L. Ellerbroek, L. Gilles, and C. R. Vogel, “Numerical simulations of multi-conjugate adaptive optics wavefront reconstruction on giant telescopes,” Appl. Opt. 42 (2003), pp. 4811–4818. [CrossRef] [PubMed]

10.

Q. Yang, C.R. Vogel, and B.L. Ellerbroek, “Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography,” Appl. Opt. 45, No. 21 (2006). [CrossRef] [PubMed]

11.

J. W. Hardy, Adaptive Optics for Astronomical Telescopes, Oxford University Press, 1998.

12.

D. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. , Vol. 67 No. 3, (1977), pp.370–375 [CrossRef]

OCIS Codes
(010.1080) Atmospheric and oceanic optics : Active or adaptive optics

ToC Category:
Focus Issue: Adaptive Optics

History
Original Manuscript: March 13, 2006
Revised Manuscript: May 10, 2006
Manuscript Accepted: June 2, 2006
Published: August 21, 2006

Citation
Curtis R. Vogel and Qiang Yang, "Fast optimal wavefront reconstruction for multi-conjugate adaptive optics using the Fourier domain preconditioned conjugate gradient algorithm," Opt. Express 14, 7487-7498 (2006)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-17-7487


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. J. M. Beckers, "Increasing the size of the isoplanatic patch with multi-conjugate adaptive optics," in Proceedings of European Southern Observatory Conference and Workshop on Very Large Telescopes and Their Instrumentation, M. H. Ulrich, ed., Vol. 30 of ESO Conference andWorkshop Proceedings (European Southern Observatory, Garching, Germany, 1988), pp. 693-703.
  2. D. C. Johnston and B. M. Welsh, "Analysis of multi-conjugate adaptive optics," J. Opt. Soc. Am. A 11, 394-408 (1994). [CrossRef]
  3. T. Fusco, J. M. Conan, G. Rousset, L. M. Mugnier, and V. Michau, "Optimal wave-front reconstruction strategies for multi-conjugate adaptive optics," J. Opt. Soc. Am. A 18, 2527-2538 (2001). [CrossRef]
  4. R. G. Dekany, M. C. Britton, D. T. Gavel, B. L. Ellerbroek, G. Herriot, C. E. Max, J-P. Veran, "Adaptive optics requirements definition for TMT," Advancements in Adaptive Optics, edited by D. B. Calia, B. L. Ellerbroek, and R. Ragazzoni, Proc. SPIE 5490, 879-890 (2004). [CrossRef]
  5. B. L. Ellerbroek, "Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques," J. Opt. Soc. Am. A,  19, 1803-1816 (2002). [CrossRef]
  6. G. Golub and C. VanLoan, Matrix Computations, 2nd Edition, Johns Hopkins University Press, 1989.
  7. 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). [CrossRef]
  8. L. Gilles, B. L. Ellerbroek, and, C. R. Vogel, "Preconditioned conjugate gradient wave-front reconstructors for multi-conjugate adaptive optics," Appl. Opt. 42, 5233-5250 (2003). [CrossRef] [PubMed]
  9. B. L. Ellerbroek, L. Gilles, and C. R. Vogel, "Numerical simulations of multi-conjugate adaptive optics wavefront reconstruction on giant telescopes," Appl. Opt. 42, 4811-4818 (2003). [CrossRef] [PubMed]
  10. Q. Yang, C.R. Vogel, and B.L. Ellerbroek, "Fourier domain preconditioned conjugate gradient algorithm for atmospheric tomography," Appl. Opt. 45, 5281-5293 (2006). [CrossRef] [PubMed]
  11. J. W. Hardy, Adaptive Optics for Astronomical Telescopes, Oxford University Press, 1998.
  12. D. Fried, "Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements," J. Opt. Soc. Am. 67,370-375 (1977). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Figures

Fig. 1. Fig. 2. Fig. 3.
 

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited