Model-based aberration correction in a closed-loop wavefront-sensor-less adaptive optics system

H. Song; R. Fraanje; G. Schitter; H. Kroese; G. Vdovin; M. Verhaegen

doi:10.1364/OE.18.024070

1. Introduction

In recent years, wavefront-sensor-less (WFSless) adaptive optics (AO) systems have been used in many scientific and medical applications, such as laser systems [1–10] and microscopes [11–18], to improve the laser beam quality or the image resolution, by correcting the air-turbulence-, heat- or specimen-induced wavefront aberrations in the optical path. Unlike the AO systems in astronomy applications [19, 20] where the wavefront aberration can be measured directly with dedicated wavefront sensors (e.g., the Shack-Hartmann WFS), there is no direct wavefront measurement in WFSless AO systems and the sensor signal (e.g., the intensity within a pin hole) is usually nonlinearly related to the wavefront aberration. Aberration correction is performed by adapting the shape of the deformable mirror (DM) such that certain performance metric (e.g., the light intensity measurement or the sharpness of the image) reaches its maximum.

Different optimization algorithms, such as gradient descent optimization algorithm, simplex optimization algorithm, genetic algorithm, simulated annealing algorithm, etc., have been used for aberration correction in WFSless AO systems and the improvements in the performance metric have been demonstrated in [1–12,14–16]. By exploring the structure of the performance metric function, model-based approaches have been proposed to speed up the correction [13, 17,18,21]. In a recent work by Débarre [18], the performance metric is locally represented as a separable quadratic function of the aberration modal coefficients by sophisticated choice of the aberration modal basis, such that N aberration modes can be corrected after 2N + 1 images.

In this paper, we further improve the correction speed of the WFSless AO system by wavefront aberration estimation and correction in three steps. First, with the external aberration absent (e.g., the aberration induced by air turbulence, heat or specimen), the WFSless AO system is calibrated such that the system aberration (e.g., initial aberration in the DM, misalignment of the optical components) is removed. Second, still with the external aberration absent, a nonlinear static model of the calibrated WFSless AO system is identified from the measurement data, which describes the transfer from the DM control signal to the intensity measurement. This step is analogue to determining the influence matrix of the DM in WFS-based AO systems; however, in WFSless AO systems, because the transfer from the DM control signal to the intensity measurement is nonlinear, a nonlinear model identification approach is required. Third, when the external aberration is present, the DM is initially excited by N + 2 predefined control signals and the corresponding N + 2 intensity measurements are collected. Aberration is estimated and corrected based on these N + 2 pairs of input-output data and the model of the WFSless AO system, by solving a nonlinear least squares (NLLS) optimization problem online. With new input-output data available, the aberration estimation and correction are refined iteratively. This approach is validated in a WFSless AO experimental setup and the performance of the resulting closed-loop system is evaluated.

The contribution of our work is that a new model-based approach has been proposed and validated for aberration estimation and correction in WFSless AO systems. The paper is organized as follows. Section 2 analyzes the WFSless AO system. Section 3 explains our approach on wavefront aberration estimation and correction. Section 4 describes the experimental setup. Section 5 reports and evaluates the experimental results. Section 6 concludes the work.

2. System analysis

The schematic of a common closed-loop WFSless AO system [2, 3, 21] under investigation is depicted in Fig. 1. The incident light beam is disturbed in front of the entrance pupil. The entrance pupil is conjugated to the DM by two lenses L1 and L2. After the beam is reflected by the DM, it is focused by the lens L3. A pin hole is placed at the focal point of L3. After the pin hole, a photodiode measures the intensity within the pin hole and feeds the intensity measurement to the control system. The control objective is to maximize the intensity measurement y(k) ∈ ℝ at time k by adapting the control signal u(k) ∈ ℝ^N to the DM, i.e,

max_{u (k)} y (k),

where u(k) can be the zonal or modal representation of the control signal, with dimension N.

Fig. 1 Schematic of a common closed-loop WFSless AO system. The incident light beam is disturbed in front of the entrance pupil. The control system adapts the control signal u(k) to maximize the intensity measurement y(k).

Download Full Size | PDF

By physical modeling [22], the intensity measurement y(k) is related to the incident wavefront aberration and the DM deformation as:

\begin{array}{l} y (k) = \iint_{Σ_{2}} | \iint_{Σ_{1}} a_{i} (ξ, η, k) exp [- j \frac{2 π}{λ} (ϕ_{x} (ξ, η, k) + ϕ_{m} (ξ, η, k))] \\ \cdot exp [- j \frac{2 π}{λ d} (α ξ + β η)] d ξ d η |^{2} d α d β + w (k) . \end{array}

Here (ξ, η) and (α, β) represent the coordinates in the input plane and the focal plane of the lens L3, respectively, and with

j = \sqrt{- 1}

. a_i(ξ, η, k) is the amplitude of the complex optical field at time instant k. ϕ_x(ξ, η, k) and ϕ_m(ξ, η, k) represent the incident wavefront aberration and wavefront manipulation by the DM at time k, respectively. The phase of the complex optical field is given by ϕ_x(ξ, η, k) + ϕ_m(ξ, η, k). λ is the wavelength of the light and d the focal distance of L3. Σ₁ represents the illuminated area of L3 and Σ₂ the area of the pin hole. w(k) is the measurement noise.

Because in many cases wavefront aberration is the main factor for intensity measurement reduction at given incident light power [12, 13, 17], the amplitude variation in the optical field is omitted such that

a_{i} (ξ, η, k) = a_{i},

where a_i is a constant. Apart from that, if the aberration is corrected within a short time, it is reasonable to consider the wavefront aberration as constant (e.g., when any single point in the specimen is imaged in scanning-type microscopes under normal operational conditions). This simplifies ϕ_x(ξ, η, k) as

ϕ_{x} (ξ, η, k) = ϕ_{x} (ξ, η),

such that ϕ_x(ξ, η) is time-independent.

The speed of aberration correction generally depends on the correction algorithm and the sampling rate of the WFSless AO system. As the sampling rate increases, the dynamics in the DM becomes more significant. Since the static nonlinearity in the intensity measurement is a common bottleneck for efficient aberration correction in WFSless AO systems while the DM dynamics is device- and sampling-rate- dependent, in this paper we focus on the static nonlinearity in the intensity measurement. Dynamics in the DM at high sampling rate is left for future research. In this case, the DM wavefront manipulation ϕ_m(ξ, η, k) can be written as

ϕ_{m} (ξ, η, k) = D (ξ, η) u (k)

where D(ξ, η) represents the static linear response of the DM. This linear representation of the DM response is valid for most commonly-used DMs when they are appropriately linearized, e.g, taking the square root of the voltage as the control signal for electrostatic-actuated DM [2], or by hysteresis compensation in piezo-driven DM [23]. Each column of D(ξ, η) can be considered as a mode of the DM deformation and u(k) contains all the modal coefficients. The column space of D(ξ, η) forms a basis for ϕ_m(ξ, η, k). Different basis can be used (e.g., DM actuator basis, Zernike basis, Lukosz basis), depending on how the control signal u(k) is defined. For instance, if u(k) is same as the voltage applied to each actuator of the DM (i.e, zonal control), then D(ξ, η) is the influence function of the DM; otherwise, Zernike modal control or Lukosz modal control can also be applied.

Combining Eq. (2)–(5), the behavior of the WFSless AO system can be represented as

\begin{array}{l} y (k) = \iint_{Σ_{2}} | \iint_{Σ_{1}} a_{i} exp [- j \frac{2 π}{λ} (ϕ_{x} (ξ, η) + D (ξ, η) u (k))] \\ \cdot exp [- j \frac{2 π}{λ d} (α ξ + β η)] d ξ d η |^{2} d α d β + w (k) . \end{array}

Static nonlinearity is visible in Eq. (6) in the wavefront-intensity mapping. Because this mapping is surjective (i.e., different wavefronts can give the same intensity measurement) and not invertible, the wavefront can not be obtained from single intensity measurement. However, with the model of the WFSless AO system describing the transfer from u(k) to y(k) with the aberration ϕ_x(ξ, η) absent, and at least N + 2 pairs of u(k) and y(k) collected with ϕ_x(ξ, η) present, the aberration ϕ_x(ξ, η) can be estimated in the basis defined by D(ξ, η), as will be explained in Section 3.

3. Model-based aberration estimation and correction

3.1. Modeling of the WFSless AO system

Because the DM deformation ϕ_m(ξ, η) can not be measured in the WFSless AO system and D(ξ, η) can not be obtained with high accuracy, it is difficult to get an accurate model of the real system from Eq. (6). The artifacts in the optical components may also degrade the accuracy of Eq. (6). As will be shown later on, since hundreds of times of intensity calculations are needed by our proposed algorithm to estimate the aberration, the computational complexity in Eq. (6) (e.g., two double integrals for each intensity calculation) will slow down the aberration correction speed. Therefore in our work the AO model is identified directly from u(k) and y(k) by black-box identification [24, 25].

In this sense, the system description in Eq. (6) is represented by

y (k) = g (ϕ_{x} (ξ, η) + D (ξ, η) u (k)) + w (k),

where g represents the static nonlinear wavefront-intensity mapping, including the double integral over the coordinate (ξ, η). The wavefront aberration ϕ_x(ξ, η) can be split into two parts as

ϕ_{x} (ξ, η) = \underset{ϕ_{1} (ξ, η)}{\underset{︸}{D (ξ, η) x}} + Δ ϕ_{x} (ξ, η) .

Here ϕ₁(ξ, η) represents the part of ϕ_x(ξ, η) lying within the range of D(ξ, η) and Δϕ_x(ξ, η) represents the part of ϕ_x(ξ, η) which is orthogonal to the range of D(ξ, η). It is assumed that the wavefront aberration can be represented by a finite low-order Zernike aberrations [13, 26, 27], then it is possible that the DM can generate these low-order Zernike modes efficiently and Δϕ_x(ξ, η) can be neglected. As a result, Eq. (8) can be approximated by

ϕ_{x} (ξ, η) \approx D (ξ, η) x .

Substitute Eq. (9) into (7), we have

y (k) \approx g (D (ξ, η) (x + u (k))) + w (k) .

Merging D(ξ, η) and g into one static nonlinear mapping f, we can further simplify the system description as

y (k) \approx f (x + u (k)) + w (k) .

Equation (11) considers the aberration as a disturbance directly applied on the input u(k), which allows to identify the model of the WFSless AO system only based on u(k) and y(k) but meanwhile accounting for the influence of the aberration.

To identify an accurate nonlinear model of the WFSless AO system from u(k) and y(k), the nonlinearity in the system should be excited persistently by the input u(k). Random signals can then be used to excite the system for data collection. Since f is identified only based on u(k) and y(k), y(k) should be collected with x = 0. If x = x₀ ≠ 0 (x₀ is an unknown nonzero constant vector) during the data collection, then there is an offset of x₀ in the estimated aberration, as will be seen in the next section. In practice, this aberration-free condition may be achieved after the calibration of the WFSless AO system, when the aberration of the WFSless AO system itself (system aberration, e.g., initial aberration in the DM, misalignment of the optical components) has been corrected and the aberration induced by external sources (e.g., air turbulence, high power heating or specimen) is still absent. The system aberration can be corrected by optimization algorithms like simplex algorithm, genetic algorithm, etc. Although optimization algorithm is used here for system aberration correction, the system aberration only needs to be corrected once during the operation of the WFSless AO system. Significant time can still be saved in correcting the external aberrations.

With the input-output data u(k) and y(k), the model structure needs to be selected for the nonlinear black-box model. There is a very rich spectrum of possible descriptions for nonlinear black-box models, e.g., neural network [28, 29], fuzzy models [30], etc. Because a 2-layer neural network is able to model a broad range nonlinearities and, from practical point of view, it can be implemented and trained with the MATLAB Neural Network Toolbox [31] very conveniently, a 2-layer neural network is built in our work, which has N_Q neurons in the first layer and one in the second. The output ŷ(k) of the neural network is determined as

\hat{y} (k) = \hat{f} (u (k)) = W_{1} tanh (W_{2} u (k) + s_{1}) + s_{2} .

W₂ ∈ ℝ^N_Q×N and W₁ ∈ ℝ^1×N_Q contain the input and output weights of the neural network, respectively; s₁ ∈ ℝ^N_Q×1 and s₂ ∈ ℝ are biases on the input and output neurons, respectively. tanh is the hyperbolic tangent function.

The number of neurons N_Q should be defined by the user when constructing the neural network. Parameters W₁, W₂, s₁ and s₂ are then optimized by training the neural network with sufficient data points u(k) and y(k). Details on training and validating the neural network can be found, for instance, in [28, 29].

3.2. Aberration estimation and correction

With the unknown aberration x present, if the WFSless AO system is excited by a certain number of inputs u(k), k = 1, ⋯ , K (K is the number of data points) and the intensity y(k), k = 1, ⋯ , K, are collected, then x can be estimated by solving a set of nonlinear equations as

{\begin{matrix} y (1) & = & a \hat{f} (x + u (1)) \\ y (2) & = & a \hat{f} (x + u (2)) \\ ⋮ & ⋮ & ⋮ \\ y (K) & = & a \hat{f} (x + u (K)) \end{matrix}

Here f̂ is the model of the WFSless AO system identified in previous step. a is a scaling factor, accounting for the possible variation in the incident light power between the modeling and aberration estimation. For instance, in microscopes, the light power emitted or reflected by the specimen may vary from point to point. The obstructing layers of the specimen may also scatter, reflect or absorb the light passing through. Although we are only interested in x for aberration correction, a should also be estimated because it is unknown in Eq. (13).

To obtain an analytical solution of Eq. (13) may be infeasible in practice, for instance, if the nonlinearity f̂ has some high-degree components. Alternatively, a numerical solution can be obtained by solving a nonlinear least squares (NLLS) problem as

(\hat{a}, \hat{x}) = \underset{\hat{a}, \hat{x}}{argmin} \underset{J (\hat{a}, \hat{x})}{\underset{︸}{{|| Y_{[1, K]} - {\hat{Y}}_{[1, K]} ||}_{2}^{2}}},

with Y_[1,K] and Ŷ_[1,K] constructed as

Y_{[1, K]} = [\begin{matrix} y (1) \\ y (2) \\ ⋮ \\ y (K) \end{matrix}], {\hat{Y}}_{[1, K]} = [\begin{matrix} \hat{y} (1) \\ \hat{y} (2) \\ ⋮ \\ \hat{y} (K) \end{matrix}] = [\begin{matrix} \hat{a} \hat{f} (\hat{x} + u (1)) \\ \hat{a} \hat{f} (\hat{x} + u (2)) \\ ⋮ \\ \hat{a} \hat{f} (\hat{x} + u (K)) \end{matrix}] .

Here â and x̂ are the estimates of a and x, respectively. For given â and x̂, the intensity is estimated by ŷ(k) = âf̂(x̂ + u(k)).

To have an efficient aberration correction, a compromise should be made in K concerning the accuracy of the aberration estimation and the correction speed. From one hand, inadequate data points can not give an accurate aberration estimation, for instance, more than one solutions may exist in Eq. (13) and the cost function J(â, x̂) in Eq. (14) does not have a unique global minimum (see Fig. 2 for an illustration). From the other hand, if more data points are collected than necessary, then the correction speed will be slowed down. A theoretical analysis on this is difficult because several factors should be considered, e.g., the nonlinearity f, the model uncertainty in f̂, the measurement noise in y(k) and the values of the K inputs. However, as a practical solution, aberration estimation and correction can be implemented in an iterative manner and the model-based aberration correction (MBAC) algorithm is described below.

Before the aberration estimation, the WFSless AO system is initially excited by N + 2 control signals u(k) and the corresponding intensity measurements y(k) are collected. Here N + 2 data points are collected for initialization concerning that N + 1 unknowns need at least N + 1 equations in Eq. (13) to have a unique solution if f were a linear function, and that nonlinear functions may need more equations in general. Since the aberration estimation and correction will be refined iteratively later on, these N + 2 data points serves as an initial trial for the MBAC algorithm. A natural option for the first control signal is u(1) = 0, i.e., no correction by the DM. The other N + 1 control signals should excite the aberrated system in such a way that rich information can be collected on the aberration x. Selection of such N + 1 inputs has been investigated in [21]. The optimum distribution of the N + 1 inputs can be geometrically interpreted as the N + 1 vertices of a regular simplex in the N-dimensional space (see Appendix B of [21]).
From time k = N + 2 on, the aberration estimation (denoted as x̂(k – 1)) is given by Eq. (14), based on previous K = k – 1 control inputs and intensity measurements. The control input is then set as u(k) = –x̂(k – 1) to counter-react on the aberration and the corresponding intensity y(k) is measured. The newly-collected y(k) and u(k) are added into Y_[1,K] and Ŷ_[1,K] respectively in Eq. (15) and the aberration estimation can be refined by solving Eq. (14) with the latest Y_[1,K] and Ŷ_[1,K]. This estimation-correction-collection procedure can be repeated iteratively. The algorithm can be stopped when a certain criterion is met, for instance, when the improvement over the previous intensity measurement is less than a certain threshold value, or when the maximum number of intensity measurements is exceeded.

Fig. 2 Cost function J(x̂) depends on the number of data points for solving the NLLS problem in Eq. (14). For clarity of explanation, the intensity variation a is not considered. In (a), the nonlinearity is represented as $y = f (x + u) = {(\frac{2 J_{1} (x + u)}{x + u})}^{2}$ to simulate the intensity distribution in the Airy disk [32], with J₁ the Bessel function of the first kind. The aberration shift the original system y = f (u) horizontally by x = −1.25. The model uncertainty is neglected, i.e., f̂ = f. With single data point P1, the cost function J(x̂) has two minima at x̂ = −1.25 and x̂ = 3.25 as plotted in (b). With points P1 and P2, J(x̂) has one unique global minimum at x̂ = −1.25 but there is a local minimum at x̂ = 2.5. This local minimum vanishes when P3 is added and the domain of convex is increased.

Download Full Size | PDF

Due to the modeling uncertainty in f̂ and the measurement noise in y(k), the accuracy of the aberration estimation may be limited and the intensity may not reach its maximum by the MBAC algorithm. In this situation, other optimization algorithms like simplex algorithm, genetic algorithm, etc., can be used to continue searching for the optimum. Under the assumption that f̂ is a close approximation of f, the MBAC algorithm will steer the DM to a point close to its optimum. This point can then be used as a new initial condition for desired nonlinear optimization method, like the simplex algorithm described in [33]. The initial simplex of the simplex algorithm is constructed around the control signal which gives the maximum intensity measurement in the MBAC algorithm. The hybrid algorithm (MBAC+Simplex) is described in pseudo code below. The MBAC algorithm stops after a fixed number of intensity measurements P (P is a user-defined number), to distinguish the intensity improvements due to the MBAC algorithm and due to the simplex algorithm. The simplex algorithm stops at time P̂ (P̂ is a user-defined number).

MBAC+Simplex algorithm (general description and pseudo code implementation):

Initialization of MBAC, i.e., collecting N + 2 data points
- Set u(1) = 0.
- Set u(k) as in Appendix B of [21], with k = 2, ⋯ , N + 2.
- Set â(k) = 1, with k = 1, ⋯ , N + 2.
- for k = 1 : N + 2
  Excite the WFSless AO system with u(k) and collect y(k).
- end
Aberration estimation and correction by MBAC
- for k = N + 3 : P
  p = argmax_p y(p) ;
  â_init = â(p), x̂_init = −u(p);
  [â(k – 1), x̂(k – 1)] = argmin_â,x̂ J(â, x̂) as in Eq.(14), with initial conditions â_init and x̂_init.
  Set u(k) = –x̂(k – 1), excite the system with u(k) and collect y(k).
- end
Aberration correction by the simplex algorithm
- p = argmax_p y(p);
- u_c = u(p);
- Construct simplex around u_c as u(k) = u(k – P + 1) + u_c with k = P + 1, ⋯ , P + N + 1.
- for k = P + 1 : P̂
  Run simplex algorithm as in [33].
- end

4. Experimental setup

The closed-loop WFSless AO experimental setup is the same as in Fig. 1. The collimated laser beam is generated by a He-Ne laser with a wavelength of 632 nm. Aberration is generated by a circular glass plate. One side of the glass plate is polished in such a way that the resulting wavefront aberration has a spatial Kolmogorov distribution [20]. The intensity transmission of the disturbance generator is about 78% as measured by a power meter (PM100, Thorlabs, Germany). During the modeling of the WFSless AO system, this aberration generator is removed. The entrance pupil has a diameter of 6 mm. It is conjugated to the PDM by lenses L1 and L2. The focal distances of L1 and L2 are 6 cm and 20 cm, respectively. The PDM (37-actuator, OKOTech, The Netherlands) has a clear aperture of 30 mm and only the central area with a radius of 20 mm is illuminated to generate Zernike modes efficiently [27]. Lens L3 has a focal distance of 400 mm. The pin hole (NT56-282, Edmunds Optics, with a diameter of 50 μm) is placed at the focal point of L3, followed by a photodiode (TSL250R-LF, TAOS, Korea) measuring the light intensity inside the pin hole. The high voltage amplifier (HVA, OKOTech, The Netherlands) has 40 channels, each with an output range of 0∼300 V, a voltage amplification of 80 at low frequencies and a −3dB bandwidth of 1 kHz. The control algorithm is implemented in MATLAB (Version 7.5.0.342, The MathWorks). Signal generation and data acquisition is accomplished by a dSPACE system (DS1006, dSPACE, Germany) with the digital-to-analog card (DS2103) output range of ±10 V, 14-bit and analog-to-digital card (DS2004) input range of ±10 V, 16-bit. Interfacing between MATLAB and the dSPACE system is done via MLIB (dSPACE, Germany).

Figure 3 depicts the block diagram of the closed-loop WFSless AO system. The physical input of the WFSless AO system is the voltage V(k) ∈ ℝ³⁷, which is applied to 37 actuators of the PDM. The output of the WFSless AO system is the light intensity measurement y(k) ∈ ℝ from the photodiode. To reduce the uncertainty in the AO setup, a hysteresis compensator Ĥ⁻¹ is implemented to compensate for the hysteresis in the PDM as described in [23]. To reduce the dimension of the control signal u(k) ∈ ℝ^N, the PDM is controlled in Zernike basis by N = 9 modes. This is accomplished by the matrix L ∈ ℝ^37×N which transforms the modal control signal u(k) to the pseudo voltage V̂(k). L is derived according to the Zernike polynomials description in [26] and the theoretical model of the PDM in [34]. The indexing of Zernike modes is the same as in [26]. Only Zernike-Mode 2 to 10 are controlled (i.e., piston is neglected). With the hysteresis compensator and the modal transformer, the WFSless AO system is conceptually considered to have the modal control signal u(k) as input and intensity measurement y(k) as output. The intensity measurement is fed into the controller and the control signal u(k) is calculated.

Fig. 3 Block diagram of the closed-loop WFSless AO system. The physical WFSless AO system has voltage V(k) as input, but conceptually u(k) can be considered as its input because of the hysteresis compensator and the modal transform.

Download Full Size | PDF

5. Experiments and results

Experiments have been carried out in the setup described in Section 4 to validate the proposed approach for aberration correction, which mainly consist of three steps as follows:

With the aberration generator absent, the WFSless AO system is calibrated using a simplex optimization algorithm. The system aberration is corrected by adapting the shape of the PDM such that the intensity measurement is maximized.
The WFSless AO system is excited by random control signals u(k) and the intensity measurements y(k) are collected. Based on u(k) and y(k), the WFSless AO system is modeled by a neural network as described in Section 3.1.
Aberration is introduced in the WFSless AO system by the aberration generator and corrected by the proposed MBAC+Simplex algorithm as described in Section 3.3. For a comparison, the simplex algorithm alone is also used to correct the aberration. Intensity improvements by these two algorithms are evaluated and compared.

5.1. System calibration

To allow for bi-directional operation of the PDM in later experiments, all the actuators in the PDM are biased by 150 V initially. A simplex optimization algorithm is then used to correct the system aberration, which maximizes the intensity measurement y(k) by adapting the control signal u(k) as in Eq. (1). The sampling rate of the system during the calibration is f_s = 50 Hz, which is much less than the resonance frequency of the PDM (about 1 kHz), so that the AO system is considered static. The maximum intensity measurement is denoted as y_max, which is used to normalize intensity measurement in Section 5.3. The control signal which results in the maximal intensity measurement, denoted as u₀, is used as a bias in all the following experiments.

5.2. Modeling of the AO system

To collect enough input-output data for modeling the WFSless AO system, the system is excited by 10000 control signals u(k) in open-loop with the aberration generator absent and the intensity measurements y(k) are collected. The control signals u(k) distribute randomly within the operational range of the PDM, to give a persistent excitation. The sampling rate of the system is also 50 Hz.

Among the 10000 collected data points, 6000 are randomly selected for identification of the AO model and the rest 4000 are for validation. The AO system is modeled as a 2-layer feedforward neural network with N_Q neurons in its first layer and one neuron in its second layer as in Eq. (12). The neural network is implemented and trained by MATLAB Neural Network Toolbox [31]. Parameters W₁, W₂, s₁ and s₂ in Eq. (12) are optimized by minimizing the mean square of the fitting error, using Levenberg-Marquardt (LM) backpropagation algorithm, i.e.,

(W_{1}, W_{2}, s_{1}, s_{2}) = arg min_{(W_{1}^{*}, W_{2}^{*}, s_{1}^{*}, s_{2}^{*})} \frac{1}{N_{t}} Σ_{k = 1}^{N_{t}} {(y (k) - \hat{y} (k))}^{2} .

N_t is the number of data points for identification, in our case, N_t = 6000.

The accuracy of the model is evaluated by calculating the variance accounted for (VAF) of the model, which is defined as

VAF (\hat{y}, y) = (1 - \frac{var (\hat{y} - y)}{var (y)}) \times 100 % .

Here var(y) is the variance of y. Figure 4 shows the VAFs of the AO model with different number of neurons in the first layer. From this plot, it can be seen that VAF already reaches as high as 98.2% at N_Q = 20 for the identification set and 97.8% for the validation set, indicating that the neural network can model the AO system very accurately. The difference in VAF is negligible for N_Q > 20. Therefore 20 neurons are used in the first layer, to have a good balance between the model accuracy and the model complexity. Experiments show that the the number of neurons N_Q needed to accurately model the system is about twice the number of modes in the system, i.e., N_Q ≈ 2N.

Fig. 4 Accuracy of the neural network model for different number of neurons N_Q. VAF increases with N_Q in both identification and validation sets for N_Q ≤ 20. The difference in VAF is negligible for N_Q > 20. Hence 20 neurons are used.

Download Full Size | PDF

5.3. Aberration correction

The aberration generator is inserted in the optical path as in Fig. 1. The MBAC+Simplex algorithm is used to correct the aberration. To have a statistics of the performance, experiments have been carried out for 20 static aberrations, which are generated by rotating the circular glass plate such that the beam is disturbed by different regions of the glass plate.

Figure 5 shows the time line of the WFSless AO system. In each experiment, during the initialization, the aberrated system is excited by N + 2 = 11 control signals u(k), k = 1, ⋯ , N + 2, at a rate of 50 Hz. Inputs u(k) are initialized as in Section 3.3. The amplitude of the simplex is selected as half of the operational range of the PDM. After the intensity y(k), k = 1, ⋯ , N + 2, are collected, the aberration is estimated by solving a NLLS optimization problem as in Eq. (14), using the function fmincon in MATLAB Optimization Toolbox. fmincon is used in our work because: (1) it is computationally very efficient and can be called in MATLAB very conveniently; (2) the convexity of J(â, x̂) improves with more data points so that a local optimization algorithm like fmincon may already be enough to get an accurate estimation â and x̂. â is constrained to be within [0, 1] during the estimation. As time keeps going, more data points are available and the aberration is estimated and corrected iteratively as in the MBAC+Simplex algorithm. After P = 19 data points, the simplex algorithm (named as Simplex 1) is switched on. For a comparison, the intensity is also maximized by the simplex algorithm alone (Simplex 2). Simplex 1 and Simplex 2 are the same except that the initial guess for Simplex 1 comes from the MBAC algorithm, but the initial guess for Simplex 2 is zero. Both simplex algorithms stop after P̂ = 200 intensity measurements, when they have converged. The sampling intervals between the 11th and the 19th samples vary because of the computational time of the NLLS algorithm, as will be discussed later. After Simplex 1 is switched on, the sampling rate returns to 50 Hz.

Fig. 5 Time line of the WFSless AO system with the MBAC algorithm, including initialization and aberration correction. The initial sampling interval is t_s = 20 ms. The computational time t_c_,1 for the first aberration estimation takes about 40 ms, while the estimation time t_c,2 afterwards takes about 20 ms because a better initial guess is provided for the solving the NLLS problem.

Download Full Size | PDF

Figure 6 shows the convergence curve for one static aberration which gives the lowest initial intensity. The intensity has been normalized as ỹ(k) = y(k)/(y_max * 0.78), where ỹ(k) is the normalized intensity and the intensity transmission ratio (78%) of the disturbance generator is accounted for. The initial intensity without correction is 0.17. After N + 2 = 11 samples are collected, the aberration is estimated and corrected by the MBAC algorithm. The intensity increases to 0.38 (about 2.2 times of the initial value) at the 12th time sample. With one more data sample acquired, the intensity jumps to 0.83 at the 13th time sample, which is almost 5 times of the initial value. At the 14th time sample, the intensity already converges to 0.86 and the intensity keeps at about 0.86 from the 15th and 19th samples.

Fig. 6 Aberration correction with the MBAC+Simplex algorithm and with the simplex algorithm alone, for one static aberration. The MBAC algorithm consists of the initialization and the aberration correction. The initial intensity is 0.17. With the MBAC algorithm, the intensity converges to 0.86 at the 14th time sample, which it takes 30 time samples for the simplex algorithm alone to reach 0.8. The simplex algorithm after MBAC also shows faster convergence than the simplex algorithm alone.

Download Full Size | PDF

The MBAC algorithm stops after 19 time samples and Simplex 1 is switched on thereafter. Simplex 1 is initialized from the 20th to the 29th time samples. The initial simplex of Simplex 1 is constructed around the input point which gave the highest intensity in the past 19 samples, as described at the end of Section 3.2. Since the initialization of the simplex algorithm is only for data collection, intensity fluctuation is observed from the 20th to the 29th time samples as expected. However, after the initialization of Simplex 1 is completed, the intensity is further improved by Simplex 1 as can be seen from the small plot in Fig. 6. This plot shows that Simplex 1 converges faster than Simplex 2 because the MBAC algorithm provides a better initial value for Simplex 1.

Figure 7 shows the convergence curve averaged over 20 experiments and the standard deviation of ỹ(k) for k ≥ 12. The initial intensity is 0.49 in average. With the MBAC algorithm, the intensity increases to 0.82 (an improvement of 67%) and 0.87 (an improvement of 78%) at the 12th and 13th time sample, respectively. The intensity converges to 0.89 at the 15th time sample, while it takes Simplex 2 about 45 time samples to reach the same level. Because Simplex 1 starts at a better initial condition provide by MBAC, the intensity reaches 0.95 at the 60th time sample, while Simplex 2 takes 90 time samples to reach the same level. A significant improvement has been achieved in correction speed. The standard deviation of ỹ(k) with MBAC is also smaller than with the simplex algorithm. For instance, at the 15th time sample, the standard deviation of ỹ(k) with the MBAC algorithm is about 0.02 while that with Simplex 2 is 0.08, about 3 times larger. This indicates that the MBAC algorithm can improve the intensity in a more deterministic manner than simplex.

Fig. 7 Correction of 20 static aberrations. The initial intensity is 0.49 in average. With the MBAC algorithm, the intensity increases to 0.82 at the 12th time sample and to 0.87 at the 13th time sample. The intensity converges to 0.89 by the MBAC at the 15th time sample, while Simplex 2 needs 45 time samples to reach the same level. The standard deviation of ỹ(k) is also reduced with the MBAC algorithm, indicating that MBAC can give a more deterministic intensity improvement than simplex.

Download Full Size | PDF

5.4. Computational complexity

Referring to Fig. 5, the computational time varies from each time when the aberration is estimated. In the first aberration estimation after 11 data samples, the cost function J(â, x̂) is evaluated for about 578 times by the function fmincon and t_c,1 is about 40 ms in average. The sampling interval between the 11th and the 12th time sample is then equal to T_s,1 = t_c,1 + t_s = 40 + 20 = 60 ms. In the aberration estimations afterwards, because a better initial guess is provided for â and x̂, the number of cost function evaluations is reduced to 251 in average and the computational time t_c_,2 reduces to about 20 ms. The sampling interval becomes T_s,2 = 20 + 20 = 40 ms.

In applications where the correction speed is the most important, the MBAC algorithm alone can be used and the correction may stop, e.g., after 15 time samples in our experiments where the intensity reaches 0.89. This leads to a total correction time of t_s × 11 + T_s,1 + T_s,2 × 3 = 400 ms, while the simplex algorithm alone needs 45 time samples (i.e., t_s × 45 = 900 ms) to reach the same intensity level. A reduction of 56% has been achieved in the correction time. If a higher intensity end value is desired, e.g., 0.95, simplex alone needs 90 time samples in average (i.e., t_s × 90 = 1.80 s). The hybrid MBAC+Simplex algorithm needs 60 time samples (19 time samples by MBAC and 41 by Simplex 1), which takes t_s × 11+ T_s,1 + T_s,2 × 7+ t_s × 41 = 1.24 s in average. The time needed by the MBAC+Simplex algorithm is only 70% of that by the simplex algorithm alone.

6. Conclusion

A new approach has been proposed for aberration estimation and correction in WFSless AO systems. The wavefront aberration is estimated by solving a NLLS problem online, based on the model of the WFSless AO system and a minimum number of N + 2 intensity measurements. Experimental results show that in average 82% of the maximum intensity can be achieved at the N + 3 = 12th time sample by the MBAC algorithm and intensity converges to 89% at the 15th time sample. With the better initial condition provided by the MBAC algorithm, the simplex algorithm also shows faster convergence than used alone.

Future work will further improve the correction speed by increasing the sampling rate of the control system and considering the dynamics of the DM.

Acknowledgments

This work is supported by Delft Center for Mechatronics and Microsystems (DCMM). We would like to thank Mr. Arjan van Dijke from TUDelft for his contribution in the implementation of the experimental setup and Dr. Niek Doelman from TNO for the Komolgorov aberration generator.

References and links

1. M. A. Vorontsov, G. W. Carhart, D. V. Pruidze, J. C. Ricklin, and D. G. Voelz, “Adaptive imaging system for phase-distorted extended source and multiple-distance objects,” Appl. Opt. 36(15), 3319–3328 (1997). [CrossRef] [PubMed]

2. G. Vdovin, “Optimization-based operation of micromachined deformable mirrors,” Proc. SPIE 3353, 902–909 (1998). [CrossRef]

3. 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(8), 1440–1453 (2000). [CrossRef]

4. W. Lubeigt, G. Valentine, J. M. Girkin, E. Bente, and D. Burns, “Active transverse mode control and optimization of an all-solid-state laser using an intracavity adaptive-optic mirror,” Opt. Express 10(13), 550–555 (2002). [PubMed]

5. U. Wittrock, I. Buske, and H. M. Heuck, “Adaptive aberration control in laser amplifiers and laser resonators,” Proc. SPIE 4969, 122–136 (2003). [CrossRef]

6. M. de Boer, K. Hinnen, M. Verhaegen, R. Fraanje, G. Vdovin, and N. Doelman, “Control of a thermal deformable mirror: correction of a static disturbance with limited sensor information,” in Proceedings of the 4th International Workshop on Adaptive Optics for Industry and Medicine, pages 61–71, Münster, Germany, 2003.

7. R. El-Agmy, H. Bulte, A. H. Greenaway, and D. Reid, “Adaptive beam profile control using a simulated annealing algorithm,” Opt. Express 13(16), 6085–6091 (2005). [CrossRef] [PubMed]

8. A. A. Aleksandrov, A. V. Kudryashov, A. L. Rukosuev, T. Yu. Cherezova, and Yu. V. Sheldakova, “An adaptive optical system for controlling laser radiation,” J. Opt. Technol. 74(8), 550–554 (2007). [CrossRef]

9. P. Yang, Y. Liu, W. Yang, M. W. Ao, S. J. Hu, B. Xu, and W. H. Jiang, “Adaptive mode optimization of a continuous-wave solid-state laser using an intracavity piezoelectric deformable mirror,” Opt. Commun. 278(2), 377–381 (2007). [CrossRef]

10. W. Lubeigt, S. P. Poland, G. J. Valentine, A. J. Wright, J. M. Girkin, and D. Burns, “Search-based active optic systems for aberration correction in time-independent applications,” Appl. Opt. 49(3), 307–314 (2010). [CrossRef] [PubMed]

11. O. Albert, L. Sherman, G. Mourou, T. B. Norris, and G. Vdovin, “Smart microscope: an adaptive optics learning system for aberration correction in multiphoton confocal microscopy. Opt. Lett. , 25(1):52–54, 2000. [CrossRef]

12. L. Sherman, J. Y. Ye, O. Albert, and T. B. Norris, “Adaptive correction of depth-induced aberrations in multiphoton scanning microscopy using a deformable mirror,” J. Microsc. 206(1), 65–71 (2002). [CrossRef] [PubMed]

13. M. J. Booth, M. A. A. Neil, R. Juskaitis, and T. Wilson, “Adaptive aberration correction in a confocal microscope,” Proc. Nat. Acad. Sci. U.S.A. 99(9), 5788–5792 (2002). [CrossRef]

14. P. Marsh, D. Burns, and J. M. Girkin, “Practical implementation of adaptive optics in multiphoton microscopy,” Opt. Express 11(10), 1123–1130 (2003). [CrossRef] [PubMed]

15. A. J. Wright, D. Burns, B. A. Patterson, S. P. Poland, G. J. Valentine, and J. M. Girkin, “Exploration of the optimisation algorithms used in the implementation of adaptive optics in confocal and multiphoton microscopy,” Microsc. Res. Tech. 67(1), 36–44 (2005). [CrossRef] [PubMed]

16. S. P. Poland, A. J. Wright, and J. M. Girkin, “Evaluation of fitness parameters used in an iterative approach to aberration correction in optical sectioning microscopy,” Appl. Opt. 47(6), 731–736 (2008). [CrossRef] [PubMed]

17. D. Débarre, E. J. Botcherby, M. J. Booth, and T. Wilson, “Adaptive optics for structured illumination microscopy,” Opt. Express 16(13), 9290–9305 (2008). [CrossRef] [PubMed]

18. D. Débarre, E. J. Botcherby, T. Watanabe, S. Srinivas, M. J. Booth, and T. Wilson, “Image-based adaptive optics for two-photon microscopy,” Opt. Lett. 34(16), 2495–2497 (2009). [CrossRef] [PubMed]

19. F. Roddier, Adaptive Optics in Astronomy, (Cambridge University Press, Cambridge, UK, 1999). [CrossRef]

20. J. W. Hardy, Adaptive Optics for Astronomical Telescopes,(Oxford University Press, New York, USA, 1998).

21. M. J. Booth, “Wave front sensor-less adaptive optics: a model-based approach using sphere packings,” Opt. Express 14(4), 1339–1352 (2006). [CrossRef] [PubMed]

22. J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, USA, 1996).

23. H. Song, G. Vdovin, R. Fraanje, G. Schitter, and M. Verhaegen, “Extracting hysteresis from nonlinear measurement of wavefront-sensorless adaptive optics system,” Opt. Lett. 34(1), 61–63 (2009). [CrossRef]

24. M. Verhaegen and V. Verdult, Filtering and System Identification: A Least Squares Approach, (Cambridge University Press, Cambridge, USA, 2007). [CrossRef]

25. J. Sjöberg, Q. Zhang, L. Ljung, A. Benveniste, B. Delyon, P. Glorennec, H. Hjalmarsson, and A. Juditsky, “Non-linear black-box modeling in system identification: a unified overview,” Automatica31(12), 1691–1724 (1995). [CrossRef]

26. M. Schwertner, M. J. Booth, and T. Wilson, “Characterizing specimen induced aberrations for high NA adaptive optical microscopy,” Opt. Express 12(26), 6540–6552 (2004). [CrossRef] [PubMed]

27. G. Vdovin, O. Soloviev, A. Samokhin, and M. Loktev, “Correction of low order aberrations using continuous deformable mirrors,” Opt. Express 16(5), 2859–2866 (2008). [CrossRef] [PubMed]

28. S. Y. Kung, Digital Neural Networks, (Prentice-Hall, Upper Saddle River, NJ, USA, 1993).

29. S. Haykin, Neural Networks: a Comprehensive Foundation, (Macmillan, New York, USA, 1994).

30. M. Brown and C. Harris, Neurofuzzy Adaptive Modeling and Control, (Prentice-Hall, New York, USA, 1994).

31. H. Demuth, M. Beale, and M. Hagan, Neural Network Toolbox 5 User’s Guide, (The MathWorks, Inc., 2007). [PubMed]

32. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, 7th ed. (Cambridge University Press, Cambridge, UK, 1999).

33. W. H. Press, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in C; the Art of Scientific Computing, 2nd ed. (Cambridge University Press, New York, USA, 1992).

34. M. Loktev, D. Monteiroa, and G. Vdovin, “Comparison study of the performance of piston, thin plate and membrane mirrors for correction of turbulence-induced phase distortions,” Opt. Commun. 192, 91–99 (2001). [CrossRef]

Model-based aberration correction in a closed-loop wavefront-sensor-less adaptive optics system

Abstract

1. Introduction

2. System analysis

3. Model-based aberration estimation and correction

3.1. Modeling of the WFSless AO system

3.2. Aberration estimation and correction

4. Experimental setup

5. Experiments and results

5.1. System calibration

5.2. Modeling of the AO system

5.3. Aberration correction

5.4. Computational complexity

6. Conclusion

Acknowledgments

References and links

Cited By

Figures (7)

Equations (17)

Optics Express