Model-based wavefront sensorless adaptive optics system for large aberrations and extended objects

Huizhen Yang; Oleg Soloviev; Michel Verhaegen

doi:10.1364/OE.23.024587

1. Introduction

Wavefront sensorless (WFSless) adaptive optics (AO) systems have been widely studied in recent years. The most important key to successful correction is an efficient control method for WFSless AO system. It is difficult to find a fast correction method because there is no direct solution between the correction elements and the metric function. Various optimization algorithms have been developed over the past decades, such as stochastic parallel gradient descent algorithm (SPGD-Stochastic Parallel Gradient Descent) [1,2], Alopex algorithm [3], simulated annealing [4], genetic algorithm [5] and recent method based on the trust region [6]. Above algorithms are iterative optimization essentially and need many intensity measurements or evaluations of the metric function, which limits the convergence rate of AO system. Even though the rapid control algorithm up to know still needs several hundred intensity measurements [6] and the number of intensity measurements increases as the turbulence strength magnifies [6,7].

It has been shown that the model-based control algorithms for WFSless AO system feature faster convergence [8–12], with only N+1 photo-detector measurements, where N is the number of corrected aberration terms. Therefore, this kind of control algorithms have a great potential in real time wavefront correction fields. Reference [9] introduces a generalized model-based control method for WFSless AO systems. Unlike the model-based method proposed by M. J. Booth, the general approach is insensitive to the selection of sets of functions as well as the bias values and is suitable not only for small, but also for large aberrations. But the method of Reference [9] assumes the imaging of a point source. In this paper, we extend the general model-based WFSless AO system to the imaging of extended sources and investigate the restoration capability and the convergence rate of the AO system for wavefront aberrations corresponding to various turbulence strength.

The paper is organized as follows. In Section 2, we prove more accurate results than those of Reference [9] and generalize it for imaging of an extended object. A model-based WFSless AO system is built with a CCD and a 61-element DM in Section 3. Closed loop simulation results for extended objects are given in Section 4.

2. Theoretical basis

The linear relation between the second moment of the intensity distribution in the focal plane and the mean-square gradient magnitude of the wavefront aberration has been shown and applied for a wavefront sensorless AO system in the paper of Linhai and Rao [9]. The original reasoning uses some general geometrical optical considerations. In this section, we prove the similar and more accurate results using physical optics, and obtain generalisation for an extended object. We then check the theoretical results on the computer simulated images of extended objects.

2.1. Second moment of an aberrated PSF

Let’s consider an incoherent imaging system, with a pupil aberration φ(x), which images the object intensity distribution o(ξ) as image intensity distribution i(u), where x = (x, y) are the pupil coordinates expressed in wavelength and ξ = (ξ, η), u = (u, v) are the angular coordinates in the object and image planes respectively. The PSF of this system is an image of a point source located on the optical axis with amplitude a, i.e. o(ξ) = aδ(ξ), where δ(ξ) is the Dirac delta function. The PSF is given by the squared absolute value of the Fourier transform of the aberrated pupil function:

I (u) = {| ℱ (a P (x) e^{- i φ (x)}) |}^{2} \equiv {| ℱ (A (x)) |}^{2},

where ℱ is the Fourier transform, P(x) is the unaberrated pupil function, and A(x) = aP(x)e⁻ⁱ^φ⁽^x⁾. Then, using the following Fourier transform properties for autocorrelation and derivative

ℱ (f ⋆ g) = {(ℱ (f))}^{*} \cdot ℱ (g),

where ⋆ denotes cross-correlation and ·* denotes complex conjugate, and

ℱ ({f^{'}}_{x} (x)) = 2 π i u ℱ (f (x)),

where

f_{x}^{'} (x)

is the partial derivative of f(x) with respect to x, one can easily get

\begin{array}{l} I (u) u^{2} = {| ℱ (A (x)) |}^{2} u^{2} = {(ℱ (A (x)) u)}^{*} \cdot (ℱ (A (x)) u) \\ = \frac{1}{4 π^{2}} {(ℱ (A_{x}^{'} (x)))}^{*} \cdot ℱ (A_{x}^{'} (x)) = \frac{1}{4 π^{2}} ℱ (A_{x}^{'} (x) ⋆ A_{x}^{'} (x)) . \end{array}

Using the identity

\int_{ℝ^{2}} ℱ (f) (u) d u = f (0),

one can finally get:

\begin{array}{l} \int_{ℝ^{2}} I (u) u^{2} d u = \frac{1}{4 π^{2}} \int_{ℝ^{2}} ℱ ({A^{'}}_{x} (x) ⋆ {A^{'}}_{x} (x)) d u \\ = \frac{1}{4 π^{2}} ({A^{'}}_{x} ⋆ {A^{'}}_{x}) (0) = \frac{1}{4 π^{2}} \int_{ℝ^{2}} {A^{'}}^{*}_{x} (x) \cdot {A^{'}}_{x} (x + 0) d x = \frac{1}{4 π^{2}} \int_{ℝ^{2}} {| {A^{'}}_{x} (x) |}^{2} d x, \end{array}

provided all the integrals exist.

The integrand in the last expression can be further transformed as

\begin{array}{l} {| {A^{'}}_{x} (x) |}^{2} = {| a {P^{'}}_{x} (x) e^{- i φ (x)} - i {φ^{'}}_{x} (x) \cdot a P (x) e^{- i φ (x)} |}^{2} \\ = {| A (x) (\frac{{P^{'}}_{x} (x)}{P (x)} - i {φ^{'}}_{x} (x)) |}^{2} = {| A (x) |}^{2} \cdot {| \frac{{P^{'}}_{x} (x)}{P (x)} - i {φ^{'}}_{x} (x) |}^{2} \\ = a^{2} P^{2} (x) ({(\frac{{P^{'}}_{x} (x)}{P (x)})}^{2} + {({φ^{'}}_{x} (x))}^{2}) = a^{2} ({({P^{'}}_{x} (x))}^{2} + {(P (x) {φ^{'}}_{x} (x))}^{2}) \end{array}

and thus

\int_{ℝ^{2}} I (u) u^{2} d u = \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} ({({P^{'}}_{x} (x))}^{2} + {(P (x) {φ^{'}}_{x} (x))}^{2}) d x .

This formula cannot be directly used for the apertures with sharp boundaries, as in this case the integral $\int_{ℝ^{2}} {({P^{'}}_{x} (x))}^{2} d x$ does not converge, and thus the second moment of the intensity distribution is not limited.

Denoting an unaberrated PSF (for which φ(x) = const and thus the second term in the integrand of Eq. (8) vanishes) as I₀, the last equation can be rewritten as

\int_{ℝ^{2}} I (u) u^{2} d u = \int_{ℝ^{2}} I_{0} (u) u^{2} d u + \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} {(P (x) {φ^{'}}_{x} (x))}^{2} d x,

or, rearranging the terms to get rid of possibly divergent integrals as

\int_{ℝ^{2}} (I (u) - I_{0} (u)) u^{2} d u = \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} {(P (x) {φ^{'}}_{x} (x))}^{2} d x,

we have proved that the second moment of the intensity distribution change in the focal plane is proportional to the integral of the square of the phase derivative multiplied by pupil function.

Repeating the same for the second image plane coordinate and adding the result to Eq. (10), one can get a relation similar to Eq. (5) of the original paper of Linhai and Rao [9], but more accurate:

\begin{array}{l} \int_{ℝ^{2}} (I (u) - I_{0} (u)) {| u |}^{2} d u = \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} P^{2} (x) ({({φ^{'}}_{x} (x))}^{2} + {({φ^{'}}_{y} (x))}^{2}) d x, \\ \equiv \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} P^{2} (x) {| \nabla φ |}^{2} d x . \end{array}

The equation is also valid for the apertures with sharp edges, for which the right part reduces to an averaged over the aperture value of the square of the phase gradient (erroneously named in the article of Linhai and Rao as its second moment).

Taking into account that total intensity does not depend on the phase if integrated over the whole image plane

\int_{ℝ^{2}} I (u) d u = \int_{ℝ^{2}} I_{0} (u) d u = a^{2} \int_{ℝ^{2}} P^{2} (x) d x,

one can get rid of unknown squared amplitude a² in the right-hand side of (11) by dividing it by the total power:

\frac{\int_{ℝ^{2}} (I (u) - I_{0} (u)) {| u |}^{2} d u}{\int_{ℝ^{2}} I (u) d u} = \frac{1}{4 π^{2}} \frac{\int_{ℝ^{2}} P^{2} (x) {| \nabla φ |}^{2} d x}{\int_{ℝ^{2}} P^{2} (x) d x} .

In a similar manner, one can prove the well-known identity widely used for restoring the phase from the Shack-Hartmann tests:

\frac{\int_{ℝ^{2}} (I (u) - I_{0} (u)) u d u}{\int_{ℝ^{2}} I (u) d u} = \frac{1}{2 π} \frac{\int_{ℝ^{2}} P^{2} (x) \nabla φ d x}{\int_{ℝ^{2}} P^{2} (x) d x},

which states that the average of the wavefront gradient over a subaperture is proportional to the change of the focal spot centroid (expressed through the first moment of the intensity).

From practical considerations, the integration domain ℝ² in the left-hand side of Eq. (11) should be limited to at least the detector size, and the equation changes to an approximation. Also, for small aberrations, the main energy of I(x) is mainly contained near the origin, and for the large values of u, the measured difference between the unaberrated and aberrated PSFs I(u) − I₀(u) is dominated by noise and should be discarded. This justifies a change of the integration kernel proposed by Booth [8] and Linhai and Rao [9] from |u|² to m(u) defined as.

m (u) = {\begin{array}{l} (1 - {| u |}^{2} / R^{2}) & if | u | \leq R \\ 0 & if | u | > R, \end{array}

which has bounded support and thus also addresses the convergence problem of the integrals. Then, if the integration domain is big enough, one can use the approximate version of identity (12)

\int_{| u | \leq R} I (u) d u \approx \int_{| u | \leq R} I_{0} (u) d u \approx a^{2} \int_{ℝ^{2}} P^{2} (x) d x,

to express the “masked detector signal” (MDS) from the paper [9] defined as

MDS = \frac{\int_{| u | \leq R} I (u) (1 - \frac{{| u |}^{2}}{R^{2}}) d u}{\int_{| u | \leq R} I (u) d u}

in terms of the mean-square wavefront gradient:

\begin{array}{l} MDS = 1 - \frac{1}{R^{2} a^{2} \int_{ℝ^{2}} P^{2} (x) d x} \int_{| u | \leq R} I (u) {| u |}^{2} d u \\ \approx 1 - \frac{1}{R^{2} a^{2} \int_{ℝ^{2}} P^{2} (x) d x} (\int_{| u | \leq R} I_{0} (u) {| u |}^{2} d u + \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} {(P (x) {φ^{'}}_{x} (x))}^{2} d x) \\ \approx 1 - \frac{\int_{| u | \leq R} I_{0} (u) {| u |}^{2} d u}{R^{2} a^{2} \int_{ℝ^{2}} P^{2} (x) d x} - \frac{1}{4 π^{2} R^{2} \int_{ℝ^{2}} P^{2} (x) d x} \int_{ℝ^{2}} {(P (x) {φ^{'}}_{x} (x))}^{2} d x \\ = {MDS}_{0} - c_{0} (R) \int_{ℝ^{2}} {(P (x) {φ^{'}}_{x} (x))}^{2} d x, \end{array}

where c₀, MDS₀ depend only on R and do not depend on the aberration. From this equation it’s clear that MDS is always less than 1, and reaches its maximum MDS₀ when there is no aberration: φ(x) = const.

A note on the approximate equality sign in Eq. (18): similarly as the aperture size defines the resolution of the imaging system, the size R of the integration domain defines the “resolution” of the phase profile φ(x) — all equations above can be made precise if one substitutes the original phase multiplied by the pupil function by the result of its convolution with the inverse Fourier image of the integration domain. From this point of view, one should always use an as large integration domain as possible, especially when large aberration gradients are expected.

2.2. Second moment for the image of a point source with arbitrary coordinates

Now consider an unaberrated image I_0,_ξ(u) of a point source with the angular coordinates ξ = (ξ, η) of the amplitude a = a(ξ) located not on the optical axis :

I_{0, ξ} (u) = {| ℱ (a P (x) e^{- i k ξ \cdot x}) |}^{2},

which is a shifted copy of I₀ :

I_{0, ξ} (u) = I_{0} (u - ξ) .

For the shifted unaberrated PSF I_0,_ξ, its second moment is related to the second moment of the central PSF I₀ by the same Eq. (11):

\int_{ℝ^{2}} I_{0, ξ} (u) {| u |}^{2} d u = \int_{ℝ^{2}} I_{0} (u) {| u |}^{2} d u + \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} {(P (x) k | ξ |))}^{2} d x

Any tip/tilt component in the pupil aberration φ(x) would be undistinguishable from the position-related phase term e⁻ⁱ^k^ξ^·^x, so it’s natural to restrict the analysis only to the aberrations with zero average tip and tilt terms (see Eq. (14)):

\int_{ℝ^{2}} {φ^{'}}_{x} (x) P^{2} (x) d x = \int_{ℝ^{2}} {φ^{'}}_{y} (x) P^{2} (x) d x = 0.

For these aberrations, one can be sure that any linear component in the phase is caused only by the position of the point source, and that the aberrated PSF I_ξ is the shifted copy of the central aberrated PSF I:

I_{ξ} (u) = I (u - ξ) \equiv {| ℱ (a P (x) e^{- i (φ (x) + k ξ \cdot x)}) |}^{2} .

Then the second moment of the aberrated PSF I_ξ is also related to the unaberrated PSF I_0,_ξ similar to Eq. (11) as:

\begin{array}{l} \int_{ℝ^{2}} I_{ξ} (u) {| u |}^{2} d u = \int_{ℝ^{2}} I_{0} (u) {| u |}^{2} d u + \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} {(P (x) (\nabla φ (x) + k | ξ |))}^{2} d x \\ = \int_{ℝ^{2}} I_{0} (u) {| u |}^{2} d u + \frac{a^{2} k^{2}}{4 π^{2}} \int_{ℝ^{2}} P^{2} (x) \cdot {| ξ |}^{2} d x \\ + \frac{2 a^{2} k}{4 π^{2}} \int_{ℝ^{2}} P^{2} (x) \cdot | ξ | \nabla φ (x) d x + \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} {(P (x) \nabla_{φ} (x))}^{2} d x \\ = \int_{ℝ^{2}} I_{0, ξ} (u) {| u |}^{2} d u + 0 + \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} {(P (x) \nabla_{φ} (x))}^{2} d x, \end{array}

i.e for the shifted PSF we have the same relation as Eq. (11):

\int_{ℝ^{2}} (I_{ξ} (u) - I_{0, ξ} (u)) {| u |}^{2} d u = \frac{a^{2}}{4 π^{2}} \int_{ℝ^{2}} P^{2} (x) {| \nabla φ |}^{2} d x .

Again, the integration domain should be large enough to include the informative part of I_ξ(u). In practice, the integration domain is limited by the image detector size, and for the point sources located close to its edge, Eq. (25) is not valid for larger aberrations (see Fig. 1).

Fig. 1 Dependence of the second moment of the point source on the averaged square of the phase gradient, plotted for the different angular coordinates of the point source. In the simulation, the CCD size was 1024×1024 pixels, the 500 turbulent phase realisation has been used as described in Section 2.4. The position of the point source was chosen to be in the centre of the image (bottom plot), and sequentially shifted 5,10,…,50 % of the image size towards the image edge (bottom to top). The plots are artificially separated for better visibility. The green lines have slope of $\frac{1}{4 π^{2}}$ . Please note the different behaviour for the two upper plots (corresponding to the point sources located 51 pixel from the image edge and on the image edge).

Download Full Size | PDF

2.3. Extended object

The incoherent image of an extended source is given by the convolution of the object intensity function and the intensity PSF [13], which is equivalent to the superposition of the images of all point sources forming the object. Therefore one can expect from the results of the previous subsection the same linear relation between the second moment of the image of an extended object and the averaged square of the phase gradient — the sum of linear plots on Fig. 1 will result in a linear plot.

More formally, denoting an incoherent (unaberrated) image of an object with intensity o(ξ) as i₍₀₎(u), and writing it as a combination of all (unaberrated) point source images with amplitudes a(ξ), where o(ξ) = a²(ξ),

i_{(0)} (u) = \int_{ℝ^{2}} I_{(0,) ξ} d ξ,

one gets

\begin{array}{l} \int_{ℝ^{2}} (i (u) - i_{0} (u)) {| u |}^{2} d u = \int_{ℝ^{2}} \int_{ℝ^{2}} (I_{ξ} (u) - I_{0, ξ} (u)) {| u |}^{2} d u d ξ \\ = \frac{1}{4 π^{2}} \int_{ℝ^{2}} o (ξ) d ξ \int_{ℝ^{2}} {(P (x) | \nabla_{φ} |)}^{2} d x, \end{array}

and after dividing the both parts by the total power

\int_{ℝ^{2}} i (u) du = \int_{ℝ^{2}} i_{0} (u) d u = \int_{ℝ^{2}} o {(ξ)}^{2} d ξ \int_{ℝ^{2}} P^{2} (x) dx,

one has

\frac{\int_{ℝ^{2}} (i (u) - i_{0} (u)) {| u |}^{2} d u}{\int_{ℝ^{2}} i (u) d u} = \frac{1}{4 π^{2}} \frac{\int_{ℝ^{2}} {(P (x) {φ^{'}}_{x} (x))}^{2} d x}{\int_{ℝ^{2}} P^{2} (x) d x},

that is the normalized change of the second moment of the image of the extended object is proportional to the mean-square gradient (denoted as MSG further) of wavefront aberration.

To keep consistence with the original terms of article [9], and to discriminate the point source and extended object imaging, we define the masked detector signal for extended sources (MDSE) as

MDSE = \frac{\int_{ℝ^{2}} i (u) {| u |}^{2} d u}{\int_{ℝ^{2}} i (u) d u},

and Eq. (29) can be written now as

MDSE - {MDSE}_{0} = c_{0} MSG,

where MDSE0 corresponds to unaberrated image and

c_{0} = \frac{1}{4 π^{2}}

. From this definition, we can see that the MDSE is always bigger than the MDSE₀.

Unlike the situation for a centred point source, the extended object contains terms I_ξ(u), for which the phase gradient magnitude is determined not only by the aberration, but also by the point position. Thus the integration domain for the extended source should be always kept as large as possible. To minimize the influence of the object points located near the image edge, one can also use the apodisation technique.

2.4. The linear dependence of MDSE on MSG of the wavefront gradient: simulation results

To verify the linear relation between MDSE and MSG of the wavefront for extended sources, we have set up an incoherent imaging model under atmospheric turbulence. The PSF is obtained by Fast Fourier Transform of the wavefront in simulation, with the dimensions chosen in such a way that the diffraction limited PSF has width of 5 pixels. Two 8-bit grey-scale images, 1024 × 1024 pixels each, referred further as “the Satellite” and “the Whirlpool”, were chosen as extended objects, where the Satellite has a clear background and the Whirlpool scatters in the whole image. The objects and their diffraction limited images are shown in Fig. 2. One 1024 × 1024 blurred image can be obtained by convolving the object image with the PSF. Considering the size of CCD in real application may vary, we set the detector size as 896 × 896, 768 × 768, 640 × 640 and 512 × 512 respectively, obtaining the corresponding sub-images by taking the central part of the biggest distorted image.

Fig. 2 Extended objects to be imaged and their diffraction-limited images

Download Full Size | PDF

For the turbulence simulation, we used the method proposed by N. Roddier [14], based on Zernike expansion of randomly weighted Karhunen-Loeve functions. The wavefront aberration is composed by Zernike modes number 3–104 (starting from the defocus term; piston and tip-tilt excluded to satisfy conditions (22)), and its size is 128 × 128 under the condition of uniform illumination. 500 phase screens with different RMS values (corresponding to the different MSG values) in the range from 0 to 1.81λ were used to test the relation between MDSE and MSG. Figures 3(a) and Figures 3(b) show two different phase screens, their PSF profiles and corresponding images of two objects respectively.

Fig. 3 Two examples of the wavefront aberrations, PSF and corresponding images. The wavefront scale is in λ; the PSF are normalized according to the ideal PSF of the imaging system.

Download Full Size | PDF

To reduce the effect of the edge pixels, we used apodisation with the Tukey window. The apodisation width α can be set to a small value for an object on a clear background, and should be increased for an object with details near the image edges. Relation between the MSG and the MDSE for the Satellite and the Whirlpool calculated for 500 random atmospheric aberrations with various MSG are given in Figs. 4(a) and Figs. 4(b) respectively. The apodisation width α was set to 0.1 for the Satellite and 0.75 for the Whirlpool. MDSE was calculated according to Eq. (30). Please note that the apodisation has affected the slope of the relation (c₀ from Eq. (31)), but hasn’t changed its linearity. For the same image, the bigger the parameter α is, the larger the deviation of the fit slopes from the theoretical value is.

Fig. 4 MDSE plotted as a function of MSG for Satellite and Whirlpool. The green lines have slope of $\frac{1}{4 π^{2}}$ .

Download Full Size | PDF

From Figs. 4(a) and Figs. 4(b), we can see that the mean square of the wavefront gradients MSG and the masked detector signals MDSE have an obviously approximate linear relation for two different extended sources and 4 different CCD sizes, and simulation values are consistent with theoretical results. The bigger the CCD size is, the better the linear relation is. At the same time, the area acquired can not be too small. There are some errors when the CCD size is 512 × 512, such as small discrepancy between the theoretical and simulation values begins to appear in Figs. 4(a) and Figs. 4(b). From the above analysis, we can get that there is a good linear dependence between MDSE and MSG when the CCD size is greater than or equal to 25 percent of the initial object area. One should acquire the area as big as possible for a better linearity in applications. However, a relatively small area provides also satisfactory results, e.g. 640 × 640 of object in the Figs. 4(a) and Figs. 4(b). We have set the CCD size is 640 × 640 to investigate wavefront correction speed and capability in the following simulation.

3. Model-based WFSless AO system

As shown in Fig. 5, the wavefront-sensorless AO system consists of a 61-element DM, an imaging lens, a CCD, and the control algorithm. The wavefront to be corrected is reflected to the imaging lens, and the lens focuses it into the CCD. The intensity is detected by the CCD and provided to the control algorithm. The control algorithm generates a correction voltage vector. The voltage vector is applied to the actuators of 61-element deformable mirror.

Fig. 5 Block diagram of WFSless AO system for extended sources

Download Full Size | PDF

We used the mirror influence functions E_k(x, y), (k = 1,…, 61) as basis functions in this paper. From the experimental measurements [15], we know the actuator influence function of a 61-element DM actuators is approximately Gaussian:

E_{k} (r) = E_{k} (x, y) = e^{\ln ω} {(\sqrt{{(x - x_{k})}^{2} + {(y - y_{k})}^{2}} / d)}^{α},

where (x_k, y_k) is the location of the k-th actuator, ω is the coupling value between actuators and is set to 0.08, and α is the Gaussian index and is set to 2. The distance between actuators is d which is about 0.1367 normalized in a unit circle. The layout of all actuators is hexagonal. We removed the tip/tilt component of the influence functions and supposed the stroke of the DM to be large enough to correct the wave-front aberrations in the simulation.

The control algorithm is divided into two parts: the preprocessing step and the iteration step. Here, the MSG of the aberrated wavefronts is shorten as S. S and its inverse matrix S ⁻¹ are calculated according to the measured influence functions E_k(x, y) of DM in the preprocessing step. The component S (i, j) of the matrix is given by

S (i, j) = D^{- 1} \iint_{D} (\frac{\partial}{\partial x} E_{i} (x, y) \frac{\partial}{\partial x} E_{j} (x, y) + \frac{\partial}{\partial y} E_{i} (x, y) \frac{\partial}{\partial y} E_{j} (x, y)) d x d y .

Matrices S, S⁻¹ and vector S_m can be calculated in advance from the influence functions; they do not depend on the wavefront to be corrected. The S_m is the diagonal vector of matrix S.

One iteration step is as follows: obtain the MDSE corresponding to the initial wavefront and record it as MDSE_init; get MDSE₁, MDSE₂,…, MDSE_k sequentially through superimposing the basis function shape introduced by DM into the wavefront to be corrected and write the coefficient vector of basis functions as β; compute the control signal ν of this iteration as

ν = \frac{S^{- 1} M}{2 c_{0} β} - \frac{β S^{- 1} S_{m}}{2},

where vector M is given by the MDSE change for each of the actuators:

M = (\begin{matrix} {MDSE}_{1} - {MDSE}_{i n i t} \\ {MDSE}_{2} - {MDSE}_{i n i t} \\ ⋮ \\ {MDSE}_{k} - {MDSE}_{i n i t} \end{matrix}),

Note that Eq. (34) is slightly different from Eq. (16) from Reference [9] because we use the relation (31) to derive the control signal ν of DM.

4. Closed loop simulation results

Based on the analysis in Section 2.4, a CCD with big size is more favourable to the wavefront correction than a CCD with a small size. In this section, we have set the CCD size to a relatively small value (640 × 640) to investigate wavefront correction speed and capability. To zero out the average tip-tilt component of the aberration (condition (22)), the 640 × 640 images were cropped relatively to the computed centroid coordinates of a bigger distorted image. The apodisation width for the Satellite was set to 0.1 and 0.5 for the Whirlpool respectively. We performed the simulation with the random phase screens calculated for five different turbulence strengths, namely D/r₀ = 5, 10, 15, 20, 25, with the corresponding initial RMS equal to 0.21λ, 0.37λ, 0.54λ, 0.68λ, and 0.81λ respectively. The averaged evolution curves over 200 frames per turbulence strength were recorded as the simulation results. The averaged RMS of the residual wavefronts and the averaged MR (Mean Radius) of the far-field intensity distribution were used to measure the correction capability of the AO system respectively. The smaller MR value is, the more focused the far field energy is and the smaller wavefront aberration is.

The adaptive curves of the control algorithm for the Satellite and the Whirlpool are given in Fig. 6 and Fig. 7 respectively.

Fig. 6 Averaged adaptation process over 200 different phase screens under different turbulence strengths for Satellite.

Download Full Size | PDF

Fig. 7 Averaged adaptation process over 200 different phase screens under different turbulence strengths for the Whirlpool.

Download Full Size | PDF

From the averaged curves of different metrics in Figs. 6(a) and Figs. 6(b) and Figs. 7(a) and Figs. 7(b), we can conclude that the AO system needs only one iretation to converge. One iteration of control algorithm means that the model-based WFSless AO system just needs 62 measurements of the CCD data in this paper. This conclusion can be extended to other DMs which have different number of actuators. The convergence speed of the system only depends on the number of actuators of DM and is independent of the turbulence strength. Therefore, this kind of control algorithm can be not only applied to static or slowly changing wavelength correction, but also to the real-time wavefront correction, provided the AO system is equipped with high-speed optoelectronic detection devices and data processor.

For investigating the correction capability of the wavefront sensorless AO system, we also calculated the averaged RMS of the residual wavefront by the ideal correction 61-element deformable mirror over 200 different phase screens under various turbulent strength. The ideal correction ability can be obtained by the minimum variance solution. The corresponding data are given in Figs. 8(a) and Figs. 8(b). RMS values of the residual wavefront aberrations after correction are very close to those of the ideal correction of 61-element DM. The maximum difference of RMS is 0.01λ, the minimum difference is 0.002λ, and the averaged relative error (∥RMS_corr − RMS_fit∥₂/∥RMS_fit∥₂, where ∥.∥₂ is the Frobenius norm) is 6%. The results above show the correction capability of DM has been used almost completely.

Fig. 8 Comparisons of initial RMS, fitted RMS and corrected RMS under different turbulences.

Download Full Size | PDF

Figures 9 – 11 give the correction effect for the wavefronts and imaging of extended source by single phase screen under different turbulence strength. Considering the length of the paper, we show correction effects for D/r₀ equal to 5,15,25 respectively. The unit of colorbar for wavefront is wavelength λ. From Fig. 11, we can see good correction results even for the strong turbulence.

Fig. 9 Correction effect of wavefront and imaging for D/r₀ = 5

Download Full Size | PDF

Fig. 10 Correction effect of wavefront and imaging for D/r₀ = 15

Download Full Size | PDF

Fig. 11 Correction effect of wavefront and imaging for D/r₀ = 25

Download Full Size | PDF

5. Conclusions

The masked intensity distribution change in the focal plane is proportional to the integral of the square of the phase derivative multiplied by pupil function for the point source. We proved this relation based on optical and mathematical theory and generalized it to the imaging of extended sources, where the mean square of the wavefront gradients MSG is approximately linear to the masked imaging of extended sources MDSE. Simulations results justified the linear dependence of MDSE on MSG of the wavefront for extended sources and the simulation values are consistent with the theoretical values when the CCD size is greater than or equal to 25% of initial object area.

To make full use of the correction capability of DM and obtain rapid convergence, we used influence functions of the deformable mirror as the basis functions of the system control algorithm. An adaptive optics system simulation platform is set up with a 61-element deformable mirror and a CCD. The convergence speed and the correction capability are investigated through correcting wavefront aberrations under different turbulence strength. Research results show that the general model-based AO system can be extended to the imaging correction of extended objects and has strong correction ability of the wavefront aberrations under different turbulence. Additionally, the convergence speed of the system depends only on the number of actuators of DM and is independent of the turbulence strength when influence functions of deformable mirror are used as basis functions.

Acknowledgments

The authors are grateful to Dr. Gleb Vdovin for the fruitful discussions. The work is partly sponsored by the European Research Council (Advanced Grant Agreement No. 339681), National Natural Science Foundation of China (Grant No. 11573011) and the Key Lab on Adaptive Optics of Chinese Academy of Science (Grant No. LAOF 201302).

References and links

1. M. A. Vorontsov, G. W. Carhart, M. Cohen, and G. Cauwenberghs, “Adaptive optics based on analog parallel stochastic optimization: analysis and experimental demonstration,” J. Opt. Soc. Am. A 17, 1440–1453 (2000). [CrossRef]

2. C. Geng, W. Luo, Y. Tan, H. Liu, J. Mu, and X. Li, “Experimental demonstration of using divergence cost-function in spgd algorithm for coherent beam combining with tip/tilt control,” Opt. Express 21, 25045–25055 (2013). [CrossRef] [PubMed]

3. M. Zakynthinaki and Y. Saridakis, “Stochastic optimization for a tip-tilt adaptive correcting system,” Comput. Phys. Commun. 150, 274–292 (2003). [CrossRef]

4. S. Zommer, E. N. Ribak, S. G. Lipson, and J. Adler, “Simulated annealing in ocular adaptive optics,” Opt. Lett. 31, 939–941 (2006). [CrossRef] [PubMed]

5. P. Yang, M. Ao, Y. Liu, B. Xu, and W. Jiang, “Intracavity transverse modes controlled by a genetic algorithm based on zernike mode coefficients,” Opt. Express 15, 17051–17062 (2007). [CrossRef] [PubMed]

6. Q. Yang, J. Zhao, M. Wang, and J. Jia, “Wavefront sensorless adaptive optics based on the trust region method,” Opt. Lett. 40, 1235–1237 (2015). [CrossRef] [PubMed]

7. H. Yang, Z. Zhang, and J. Wu, “Performance comparison of wavefront-sensorless adaptive optics systems by using of the focal plane,” International Journal of Optics, 1–8 (2015). [CrossRef]

8. M. J. Booth, “Wavefront sensorless adaptive optics for large aberrations,” Opt. Lett. 32, 5–7 (2007). [CrossRef]

9. H. Linhai and C. Rao, “Wavefront sensorless adaptive optics: a general model-based approach,” Opt. Express 19, 371–379 (2011). [CrossRef] [PubMed]

10. L.-H. Huang, “Coherent beam combination using a general model-based method,” Chinese Phys. Lett. 31, 094205 (2014). [CrossRef]

11. J. Antonello, T. van Werkhoven, M. Verhaegen, H. H. Truong, C. U. Keller, and H. C. Gerritsen, “Optimization-based wavefront sensorless adaptive optics for multiphoton microscopy,” J. Opt. Soc. Am. A 31, 1337–1347 (2014). [CrossRef]

12. J. Antonello and M. Verhaegen, “Modal-based phase retrieval for adaptive optics,” J. Opt. Soc. Am. A 32, 1160–1170 (2015). [CrossRef]

13. J. Goodman, Fourier Optics (Roberts & Company, 2005).

14. N. A. Roddier, “Atmospheric wavefront simulation using zernike polynomials,” Opt. Eng. 29, 1174–1180 (1990). [CrossRef]

15. W. Jiang, N. Ling, X. Rao, and F. Shi, “Fitting capability of deformable mirror,” Proc. SPIE 1542, 130–137 (1990). [CrossRef]

Model-based wavefront sensorless adaptive optics system for large aberrations and extended objects

Abstract

1. Introduction

2. Theoretical basis

2.1. Second moment of an aberrated PSF

2.2. Second moment for the image of a point source with arbitrary coordinates

2.3. Extended object

2.4. The linear dependence of MDSE on MSG of the wavefront gradient: simulation results

3. Model-based WFSless AO system

4. Closed loop simulation results

5. Conclusions

Acknowledgments

References and links

Cited By

Figures (11)

Equations (35)

Optics Express