## Abstract

A new wavefront sensorless adaptive optics method is presented that can accurately correct for time-varying aberrations using a single focal plane image at each sample instance. The linear relation between the mean square of the aberration gradient and the change in second moment of the image forms the basis of the presented method. The new algorithm results in significant improvements when an accurate model of the aberration’s temporal dynamics is known, by applying a Kalman filter and optimal control. Moreover, where existing wavefront sensorless adaptive optics methods update all modes sequentially, the information of the Kalman filter is used to select and update the modes that are expected to give the greatest improvement in performance. The performance is analyzed in a simulation of an adaptive optics system for atmospheric turbulence. The results show that the new method is able to correct for the aberration more accurately for higher wind speeds and higher noise levels than existing algorithms.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

Wavefront sensorless (WFSless) adaptive optics (AO) systems are systems in which the aberrations of the wavefront have to be corrected without using a dedicated wavefront sensor (WFS). Instead, only the images of a focal plane camera are used. The correction is applied to a deformable mirror (DM) in order to minimize the effect of the aberration on image quality. Finding an accurate correction without a WFS is challenging because of the nonlinearity of the underlying optimization problem. Various WFSless AO algorithms have been developed [1–4]. The common features of these methods are that they are all iterative and require many measurements to converge. Recently, WFSless AO also has gained attention in free-space optical communication leading to new developments [5,6].

Alternatively, an accurate correction can be found by solving the phase retrieval problem [7]. However, this either requires additional constraints, such as knowledge of the field’s amplitude in the pupil plane, or requires multiple simultaneous measurements at different positions along the optical axis. Furthermore, the phase retrieval problem for AO systems is computationally very demanding and would limit the sampling frequency of the control loop significantly [8].

A recent development in WFSless AO is an approach that requires only $m+1$ measurements, where $m$ is the number corrected modes [9,10] and is often referred to as model-based or second-moment (SM)-based WFSless AO. A modal basis is used that is spanned by the influence functions of the DM. This type of method has been shown to converge faster than other optimization algorithms [11]. The key to this approach is the linear relation between the mean square of the phase aberration gradient and the change in SM of the point-spread function (PSF). By exploiting this linear relation, a closed-form expression of the correction can be computed. In contrast with iterative algorithms, this method is, due to its fast convergence, promising for real-time AO applications in which the aberration is time varying.

However, this method is still useful only when the aberrations are static or change very slowly over time. Since the method relies on the assumption that the aberration does not change upon taking all $m+1$ measurement images, the performance decreases rapidly when the aberrations move faster. In [11], a method is used that requires the aberration not to change over only two measurement images. However, this method does not aim to predict the evolution of the aberration over time, nor discusses the effect of measurement noise on the methods’ performance. Modeling the aberration’s temporal behavior and the application of Kalman filter theory have been proven to be successful to deal with time-varying aberrations and measurement noise in the field of AO for astronomy [12,13]. Therefore, the effects of including a temporal aberration model and Kalman filter for WFSless AO are studied in this paper.

This paper presents an extension of the methods described in [9,11,14] for aberrations that continuously change over time and for which an accurate temporal model is available or can be identified. An example of an application that deals with these types of aberrations is AO for astronomy [12]. From now on, these aberrations will be referred to as dynamic or time-varying aberrations. When a dynamic model of the aberration is available, Kalman filter theory and optimal control are applied to close the loop and compute the optimal DM commands. The Kalman filter is used to predict the aberration in the future and to select and update the measurement image(s) that are expected to give the greatest improvement in performance. The performance will be analyzed in a simulation study of an AO system that corrects for aberrations caused by atmospheric turbulence and is compared to the method in [14].

The remainder of this paper is structured as follows. First, the existing framework of SM-based WFSless AO is explained, and its limitation for dynamic aberration correction is discussed. Section 3 presents the new SM-based WFSless AO method for time-varying aberrations. The simulation environment is discussed in Section 4, which considers a case study of dynamic aberrations caused by atmospheric turbulence. The performance of the new method is presented in Section 5, where it is compared to the existing SM-based WFSless AO method [14]. The main conclusions of this paper are summarized in Section 6.

## 2. SECOND-MOMENT-BASED WAVEFRONT SENSORLESS ADAPTIVE OPTICS

A SM-based approach to WFSless AO has already been shown to outperform the iterative algorithms in terms of convergence speeds when applying correction for a static aberration. Recently, this was also applied to the case of dynamic aberrations [14]. For the completeness of this paper, this section summarizes the theory of this SM-based WFSless AO approach. A scheme of a WFSless AO setup is shown in Fig. 1.

The main goal of WFSless AO is to correct for a residual phase aberration in the pupil plane $\phi (\boldsymbol{\chi})$, with pupil plane coordinates $\boldsymbol{\chi}=[{\chi_1}\;{\chi_2}{]^T}\in {{\mathbb R}^2}$. Define the PSF $I(\boldsymbol{\zeta},\phi)$ with image plane coordinates $\boldsymbol{\zeta}\in {{\mathbb R}^2}$ as

#### A. Static Aberration Correction

Define the DM influence functions as ${E_j}(\boldsymbol{\chi})$ for each actuator $j=1,\ldots,m$. Each actuator is poked independently, such that this phase will be added to the existing residual aberrated wavefront $\phi (\boldsymbol{\chi})$. It is important to notice that, since a residual wavefront is observed by the focal plane camera, computed control signals will always be incremental and have to be added on top of the current control signal.

An inherent property of the SM-based WFSless AO approach for control is that it will always be unable to compensate for the part of the aberration that is orthogonal to the DM basis. However, the true aberration, $\tilde\phi (\boldsymbol{\chi})$, actually is

where such that $\phi (\boldsymbol{\chi})$ represents the part of the aberration within the span of the DM influence functions, and ${\phi_\bot}(\boldsymbol{\chi})$ is the part orthogonal to it. Generally, in the derivation of the control law, it is assumed that $\tilde\phi (\boldsymbol{\chi})=\phi (\boldsymbol{\chi})$, such that only $\phi (\boldsymbol{\chi})$ will be estimated. However, it has to be noted that this is an approximation and that ${\phi_\bot}(\boldsymbol{\chi})$ will have an influence on the measurements in practice [15]. Also in the simulations of this paper, ${\phi_\bot}(\boldsymbol{\chi})$ will affect the measurements.SM-based WFSless AO aims to find the coefficients ${u_j}$ describing the aberration most accurately. In the general SM-based WFSless AO method [9], each actuator is poked sequentially. When poking actuator $j$ with an amplitude $\beta$ on top of the current control signals, the total phase of the field measured by the camera will be

Define the matrix $S\in {{\mathbb R}^{m\times m}}$ and vector $\boldsymbol{s}\in {{\mathbb R}^m}$ as#### B. Challenges for Time-Varying Aberrations

When the method described in Section 2.A is applied to time-varying aberrations, some difficulties arise. Denote the dynamic aberrated wavefront by $\phi (\boldsymbol{\chi},t)$, with $t$ describing the current discrete time instance. Furthermore, assume it takes ${T_i}$ seconds to collect an image and compute its SM. A total of $m+1$ images have to be taken. First, a reference image is taken, followed by $m$ images, each corresponding to poking a different actuator. When sequentially poking all actuators and including the fact that $\phi (\boldsymbol{\chi},t)$ is time varying, the total phase aberration corresponding to the image with the $j$th actuator poked, previously Eq. (7), becomes time varying:

One way to decrease the effect of this time delay is to update the actuators one by one, taking a new reference image in between. However, the linear system of Eq. (10) is in general not a decoupled system. A new basis can be formulated in order to have a diagonal matrix $C$. Such a diagonalization is proposed in [11] and used for dynamic aberrations in [14]. For the completeness of this paper, this is shortly summarized in the next subsection. It should be noted that, although the maximum time difference between the images can be reduced from $(m+1){T_i}$ to just $2{T_i}$, not all problems are resolved. First of all, there is still a (small) delay between the two images that is not taken into account in [11]. Second, there is a time of $2m{T_i}$ between updating the same mode again. These problems will be treated by the new method proposed in Section 3.

#### C. Diagonalizing the Linear System

The closed-form solution of Eq. (11) can be separated along the elements of $\boldsymbol{u}$ when ${C_m}$ is diagonal. The singular value decomposition (SVD) of $S$,

can be used to formulate a new basis that results in a decoupled linear equation. $\Sigma$ is a diagonal matrix and is shown to correspond to the correlation matrix [formerly $S$ for the basis of Eq. (6)] for the basis described by basis functions ${\tilde E_j}(\boldsymbol{\chi})=E(\boldsymbol{\chi}){U_j},\;j=1,\ldots,m$ [11], where $E=[{E_1}(\boldsymbol{\chi}),{E_2}(\boldsymbol{\chi}),\ldots,{E_m}(\boldsymbol{\chi})]$ and ${U_j}$ the $j$th column of $U$. Consequently, the update in terms of the new basis becomes where ${\boldsymbol{y}_d}$ contains the changes of the SM of the PSF corresponding to actuating the new set of modes ${\tilde E_j}(\boldsymbol{\chi})$, ${C_d}=2\beta {c_0}\Sigma$, and ${\boldsymbol{y}_{d,0}}$ is a vector containing the diagonal elements of ${\beta^2}{c_0}\Sigma$. Since $C_d^{-1}$ is diagonal, it is no longer necessary to wait until $m$ measurement images are taken, but the elements in ${\boldsymbol{u}_d}$ can be updated after a reference image and a single measurement image.## 3. PREDICTIVE SECOND-MOMENT-BASED WAVEFRONT SENSORLESS AO

In order to have a more systematic approach of dealing with time-varying aberrations, a dynamic model that exploits the spatio-temporal relations in the wavefront will be used. Combining the aberration dynamics with Eq. (10), a linear state-space model is obtained. The pupil plane coordinates $\boldsymbol{\chi}$ are discretized and sampled on an $n$-by-$n$ grid, such that the wavefront at time instance $t$ can be represented by the vector $\boldsymbol{\phi}(t)\in {{\mathbb R}^{{n^2}}}$. Similarly, the DM influence functions ${E_j}(\boldsymbol{\chi})$ are sampled on the same square grid, and each function is represented by the vector ${\boldsymbol{e}_j}\in {{\mathbb R}^{{n^2}}}$. The influence matrix is defined as $E=[{{\boldsymbol{e}_1},\;{\boldsymbol{e}_2},\ldots,\;{\boldsymbol{e}_m}}]\in {{\mathbb R}^{{n^2}\times m}}$. For reasons discussed in the previous section, it is assumed that the wavefront can be written as Eq. (6), such that only the part of the aberration within the span of the DM influence functions will be modeled.

First, without considering the time index yet, a general modal basis with coefficient vector $\boldsymbol{x}\in {{\mathbb R}^m}$ is defined to represent any wavefront $\phi (\boldsymbol{\chi})$ in Eq. (6). The relation between the coefficient vector $\boldsymbol{u}$ and $\boldsymbol{x}$ of either basis is given by an invertible matrix $B\in {{\mathbb R}^{m\times m}}$, such that

Using this modal basis, the observable part of the phase, $\boldsymbol{\phi}$, can per definition be modeled as In the simulations presented in this paper, both $\boldsymbol{x}$ and $\boldsymbol{u}$ will correspond to the same modal basis, i.e., $B=I$. It should be noted that although both $\boldsymbol{x}$ and $\boldsymbol{u}$ represent coefficients belonging to a basis spanning the same space, they will be used to represent different processes $\boldsymbol{x}(t)$ and $\boldsymbol{u}(t)$ in the remainder of this paper. This will be clarified in Sections 3.A and 3.B, in which a dynamic model of the wavefront will be expressed in terms of $\boldsymbol{x}(t)$, and a Kalman filter is derived.The output $\boldsymbol{y}(t)\in {{\mathbb R}^p}$ will be similar to Eq. (10). Due to the use of a Kalman filter, it is possible to update the DM with a smaller number of images by poking a selection of actuators, i.e., any $p\le m$ can be taken without diagonalizing the system as in Section 2.C. The exact definition of $\boldsymbol{y}(t)$ and the reason that $p$ can be different from $m$ will be explained further in the following subsections. Since the collection of one image is done once every ${T_i}$ seconds, the total time it takes to collect the data for the measurement vector $\boldsymbol{y}(t)$ is $(p+1){T_i}$ seconds. The output is updated every $T=(p+1){T_i}$ seconds, while the input will still be updated every ${T_i}$ seconds. As a result, the model becomes a so-called multi-rate linear time-invariant system, where the input and output are obtained over different sample periods. Although the input and output sampling rates are different, they are uniformly sampled, and the sample times coincide every $p+1$ samples. A schematic representation of one output sample time is shown in Fig. 2.

The next subsection will discuss the temporal model of the open-loop aberration. Afterwards, the resulting closed-loop state-space system and a Kalman filter implementation are discussed followed by an optimal controller using the predictions from the Kalman filter. In the last subsection of this section, it is explained how the information given by the Kalman filter can be used to select which actuator to poke for the next measurement. Table 1 gives an overview of important notations that are used throughout this section.

#### A. Dynamic Aberration Model

The temporal dynamics of the aberration caused by the turbulence, denoted by ${\boldsymbol{x}_t}(k)\in {{\mathbb R}^m}$, will be described by a vector auto regressive (VAR) model of order 1. Two different models are defined. One has the output sample time $T$, and the other has the input sample time ${T_i}$. The two models are the following:

where ${\boldsymbol{w}_f}(k)$ and $\boldsymbol{w}(k)$ are zero-mean Gaussian processes with covariance matrices ${Q_f}$ and $Q$ respectively, i.e., ${{\boldsymbol{w}}_f}(k)\sim {\cal N}(0,{Q_f})$ and $\boldsymbol{w}(k)\sim {\cal N}(0,Q)$. The simulations in this paper will focus on the example of aberrations introduced by atmospheric turbulence. However, it should be noted that the proposed method works for any type of dynamic aberration that can be accurately represented by this type of model.The matrices $A,\;{A_f},\;Q$, and ${Q_f}$ can be derived from first principles, requiring knowledge of turbulence statistics, wind direction, and wind speed. When this knowledge is not available, the system matrices can be identified when a sufficiently large dataset of open-loop aberration data $\{{\boldsymbol{x}_t}(k{T_i}),\;k=1,2,\ldots,{N_{id}}\}$ is available. With this data, identification of the matrices $A,\;{A_f},\;Q$, and ${Q_f}$ follows from a linear least-squares problem [16]. The simulations in this paper will use a simulated identification dataset to identify the system dynamics. The exact method to obtain this identification dataset in practice is beyond the scope of this paper. Since the model identification is done offline, a WFS can be temporally included in the AO system, or computational complex methods can be used to collect this dataset. When a WFS is included in the system during the identification data collection, the WFS measurements can be used to reconstruct a wavefront that is then mapped onto the desired modes to form the identification dataset. Alternatively, when additional constraints are available, such as a sparsity constraint or knowledge of the amplitude in the pupil plan, or when multiple images along the optical axis can be taken, there are existing methods that can obtain the identification dataset from solving the phase-retrieval problem on a time series of focal plane images (see [7] for an overview).

#### B. Kalman Filter Implementation

The loop is closed by the DM. The influence of the DM on the wavefront is defined as

where the delay represents the fact that the DM cannot respond instantaneously, and $B$ is the transformation matrix defined in Eq. (15). In closed loop, the residual aberration, denoted by $\boldsymbol{x}(t)\in {{\mathbb R}^m}$, is defined as For the output $\boldsymbol{y}(kT)\in {{\mathbb R}^p}$, only a selection of images is taken and processed. Within each output sample time, a selection of $p$ actuators is poked and the corresponding images are taken. Defining ${\cal I}(kT)=\{{{\cal I}_1},\ldots,{{\cal I}_p}\}\subseteq\{1,2,\ldots,m\}$ as the set of $p$ distinct integer elements corresponding to the actuators that will be poked for the next measurement, the output equation $\boldsymbol{y}(kT)$, based on Eq. (10), isCombining Eq. (21) with Eqs. (18)–(20), the following single-rate state-space model can be derived, which is sampled at the output sampling rate:

To arrive at an optimal prediction of the state vector, the Kalman filter essentially performs two steps: a measurement update, in which a newly obtained measurement is used to improve the current estimate of the state vector $\boldsymbol{x}(t)$, and a time update, where the model is used to predict that state vector. The measurement update is computed every output sample time and is given by

#### C. Optimal Control

The derived optimal state prediction in Eqs. (25) and (26) is used to create an optimal controller. The control law aims to minimize the norms of the predicted residual wavefront coefficients, i.e.,

Solving Eq. (27) and using $\boldsymbol{\hat x}(kT+j{T_i}|kT)=0$ for $\,j=1,2,\ldots,p+1$ gives the following simplified optimal control actions:

#### D. Actuator Selection Algorithm

Besides a more accurate prediction, having a model of the aberration dynamics results in another important advantage of this method. In [14], all modes are actuated sequentially, and this is repeated after the last mode is actuated. In this section, a method is proposed that uses the information from the Kalman filter, rather than sequentially poking all modes. In other words, it uses the information in the Kalman filter to decide which set ${\cal I}(kT)$ in Eq. (21) will give the most informative measurements.

The selection method is based on the realization that the state error covariance matrix of the Kalman filter,

The basis of the actuator selection method in this paper lies in the fact that both the time and measurement update of the state error covariance do not require an actual measurement, but can be computed from the so-called Riccati equation [16]. Consequently, at time $t=kT$, the values of, e.g., $P(kT+T|kT+T)$, can be computed before actually observing the measurement at $t=kT+T$, assuming the selection of actuators ${\cal I}(kT+T)$, i.e., $C(kT+T)$ and ${\boldsymbol{y}_\textbf{0}}(kT+T)$ in Eq. (21), is known.

Of course, $I(kT+T)$ is not known, but is the unknown still to be determined. However, it is known that there is only a finite number of sets possible. Theoretically, it is possible to compute $P(kT+T|kT+T)$ for all possible sets ${\cal I}(kT+T)$, but for larger values of $p$, there are too many possible combinations of actuators, and this is therefore not practical. However, it is possible to compute the state error covariance matrix corresponding to poking a single actuator, i.e., $P(kT+2{T_i}|kT+2{T_i})$. This leads to the following simple algorithm.

The first step is to compute, at $t=kT$, all possible values of $P(kT+2{T_i}|kT+2{T_i})$. This requires solving the Riccati equation $m$ times, such that $m$ different matrices are obtained. The second step is to compare the value of ${\rm trace}[{P(kT+2{T_i}|kT+2{T_i})}]$ for each matrix. The $p$ actuators that have led to the $p$ smallest values of ${\rm trace}[{P(kT+2{T_i}|kT+2{T_i})}]$ are defined as the next set of actuators to be updated, ${\cal I}(kT+T)$.

As a result of this actuator selection method, it is expected that the actuators located in an area where the dynamic model of the aberration is less accurate are updated more frequently. The Kalman filter including the optimal control law and actuator selection method are implemented in a simulation study that is discussed in the following section.

## 4. SIMULATION OF AO FOR ATMOSPHERIC TURBULENCE

The performance of the method is shown in a simulation of an AO system for atmospheric turbulence, where the aberrations in the wavefront shown in Fig. 1 are caused by atmospheric turbulence and are compared to another SM-based WFSless AO algorithm for dynamic aberrations [14]. This section will discuss the simulation conditions and the implementation details of the algorithm. Table 2 summarizes the most important simulation parameters. The results of the simulations are presented in Section 5.

#### A. Adaptive Optics Simulation Design

The phase aberrations caused by atmospheric turbulence is simulated using Oriented Matlab Adaptive Optics (OOMAO) [17]. A single turbulence layer is considered with Fried parameter ${r_0}$, outer scale ${L_0}$, and wind speed $v$. The telescope diameter will be fixed at $D=1\,\,\rm{m}$ and sample frequency at ${f_s}=1000\,\,\rm{Hz}$, i.e., ${T_i}=\frac{1}{1000}\,\,{\rm s}$. In order to have a fair comparison between different wind speeds, a collection of wavefronts on an $n\times n$ grid of pixels has been generated at a speed of $1\,\,{\rm pixel}/T_i$ over a period of $\tilde N$ time samples. The simulated sequence of wavefronts is stored in a three-dimensional tensor of dimension $n\times n\times\tilde N$. Afterwards, linear interpolation along the third dimension of this tensor is performed to obtain wavefronts at slower or faster wind speeds. A set of 20 turbulent wavefronts, each containing $N=2000$ time samples, is created for each combination of parameters, and the performance of the algorithm is tested in a Monte Carlo simulation.

The DM consists of a square grid of ${m_1}\times {m_1}$ actuators with Gaussian influence functions:

where $d$ is the distance between actuators in the pupil plane, and $\lambda\gt 0$ is the coupling parameter, defining the width of the functions. A schematic representation is shown in Fig. 3. During the simulations in this paper, ${m_1}=7$, $\lambda=0.1$, and $d=D/({m_1}+1)=0.125\,\,\rm{m}$. The actuators at the corners of the square are removed, since they have little influence inside the circular aperture, such that a total of $m=37$ active actuators are used. The influence of the actuator on the closed-loop aberration is simulated as in Eq. (20).The PSF images are simulated according to the definition in Eq. (1) and sampled on an equally spaced square grid. As a result, the SM of this discretized image becomes a weighted sum of the PSF pixel values.

#### B. Measurement Noise

The camera noise in each pixel can be modeled as a combination of Gaussian (read-out noise) and Poisson noises (shot noise). For pixels with large mean values, Poisson distributions can be accurately approximated by a Gaussian distribution. For small pixel values, the influence of the shot noise becomes less significant with respect to the Gaussian read-out noise. Furthermore, low-valued pixels have a very low SNR and are therefore often truncated when processing the images. Since the SM is a weighted sum of the pixel values, the noise contribution to the SM is expected to be approximately Gaussian and can be estimated based on the intensity measurement and camera properties.

In practice, methods such as truncation of low-valued pixels and filtering have to be used to decrease the effect of the measurement noise on the SM. However, it has been decided that this will not be included in the analysis for this paper. Instead, in order to have a clear analysis of the noise sensitivity of the algorithm itself, the noise signal $\boldsymbol{v}(kT)$ will be simulated directly as a Gaussian noise, $\boldsymbol{v}(kT)\sim {\cal N}(0,\sigma_y^2I)$.

## 5. SIMULATION RESULTS

The results of the simulations will be discussed in this section. A number of parameters will be varied in order to study the performance of the methods under different circumstances. First, the wind velocity, $v$, is varied to see whether the improvements of the new method are indeed more clear for faster moving turbulence. Second, the Fried parameter, ${r_0}$, is changed to see which method is better to deal with more severe aberrations. Third, the influences of increasing the number of actuators that is updated per measurement, $p$, are discussed. Finally, the noise sensitivity is analyzed by varying the parameter ${\sigma_y}$ in the measurement noise covariance $R=\sigma_y^2I$. The rest of the parameters will be kept constant in order to have a fair comparison. The standard values of the parameters can be found in Table 2.

In the legends of the figures, the new method will be referred to as *Dynamic SM*, as it includes a temporal dynamic model of the aberrations plus the SM model Eq. (10). It is compared to an existing method [14], which is referred to as *Static SM*, as it is based on the assumption that the change in aberration between two input sampling times is negligible, i.e., the aberration is approximately static over this time period. The parameter $\beta$ is seen as an important tuning parameter in both methods. When $\beta$ is taken too small, the SNR will be very low. Too large values of $\beta$ will create too many speckles in the PSF for the output to still be informative. Therefore, it has an optimal value that depends on the current simulation conditions and measurement noise level. During the following simulations, $\beta$ has been tuned for each different set of parameters using a grid search in order to improve the performance for both the static and dynamic SM-based method.

The main goal is to decrease the effect of the aberration on the image. Therefore, the Strehl ratio is used as a measure of performance. The Strehl ratio is computed as the fraction of the maximum pixel value of the aberrated PSF over the maximum pixel value of the unaberrated PSF. A higher Strehl ratio indicates a better performance. Also, since the method is based on minimizing the SM of the PSF, this will be occasionally used as a measure. Since each simulation takes $N$ time samples, the mean Strehl value and SM for each simulation will be used, and the Monte Carlo simulation is visualized in boxplots.

#### A. Increasing Wind Speed

It is expected that the new method, due to its predictive capabilities and smart actuator selection algorithm, can handle much higher wind speeds than the static method for the same AO system. It should be noted that for both the static and dynamic SM-based methods, besides the Fried parameter ${r_0}$ and sample frequency ${f_s}$, the performance at a certain wind speed is influenced by the number of actuators $m$ and the spatial sample distance of the actuators $d$. This is due to the fact that only one actuator can be poked at the same time. So when $m$ is increased, more modes have to be updated within the same time span, and when $d$ is decreased, the turbulence moves the distance between actuators in a shorter time. Therefore, in order to improve the performance under higher wind speeds or when increasing the number of DM actuators, it might be necessary to increase the sample frequency ${f_s}$. For the simulations in this paragraph, $m$ and ${f_s}$ are fixed at the values in Table 2.

The wind speed is varied between 0.2 m/s and 16 m/s, and the other parameters are kept constant at their standard values according to Table 2. The results of this simulation are shown in Fig. 4. From observing the Strehl ratio in Fig. 4(a), the new method is clearly better, most notably for higher wind speeds. Even for lower wind speeds, the new method outperforms the existing approach. This can be explained by the fact that the Kalman filter is more suited in dealing with measurement noise. In Fig. 4(b), the SM of the PSF is shown. Even when the static method improves the image quality in terms of the Strehl ratio at lower wind speeds, the SM of the PSF is larger than in the case without AO. This is caused by the fact that under the noisy conditions in this simulation, the static method introduces high-frequency aberrations erroneously in the compensation. The oscillating trend of the static method’s line is due to the different step lengths of the interpolation when generating the wavefront data as discussed in Section 4.A. Interpolation of a time-varying aberration has in general a slightly smoothing effect on the SM of the aberrated PSF. This effect is usually very small, as is seen in the line corresponding to no control, but seems to be amplified by the static method. Since this smoothing effect can be seen as a small amount of noise on the data, this issue indicates again the high noise sensitivity of the static method compared to the dynamic method.

It should be noted that due to the small number of DM actuators, the fitting error is relatively large, and the maximum achievable Strehl ratio for this AO configuration is approximately 0.83. So at very low wind speeds, the new method approaches the theoretical optimum. In order to improve the theoretical maximum closed-loop performance, the fitting error should be decreased by using a DM with more actuators. However, as discussed in the beginning of this paragraph, having more actuators within the same size aperture corresponds to more modes to be updated in the same time. This would require that also the sample frequency ${f_s}$ is increased in order to maintain a good performance for high wind speeds.

#### B. Number of Images in Output

In the previous simulation, it was assumed that $p=3$, such that three output channels are observed within one output sample time. The method allows for any $p$ that satisfies $1\le p\le m$. Varying $p$ can possibly lead to a better performance and needs to be investigated. Decreasing $p$ corresponds to a smaller output sampling time $T$, having the advantage that the measurement updates are more frequent. However, the drawback is that, since per $(p+1){T_i}$ seconds $p$ output channels are created, more time is spent per acquired output channel when $p$ is small. On the other hand, increasing $p$ increases the output sampling time $T$, but means that less time is spent per acquired output channel. For the sake of a complete comparison, a multivariate output is also considered for the existing method of [14]. It should be mentioned that in the original method, varying $p$ was not discussed, and $p$ was chosen to equal 1.

The results for varying $p$ while keeping the other parameters equal to their standard value in Table 2 are shown in Fig. 5(a). It clearly shows how $p=1$ is not the optimal value in this case, but larger values are more optimal. The effect of $p$ for other wind speeds is shown in Fig. 5(b). Only the medians over the Monte Carlo runs are plotted. Although all lines have similar trends, it is visible that $p=1$ does not always lead to the best results. Especially in the region where the wind speed is low, $p=5$ is better than $p=1$ as is also clear in Fig. 5(a).

The effect of the actuator selection method described in Section 3.D is very clearly visible in Fig. 6. For the same actuator configuration as in Fig. 3 and a single turbulence layer moving from left to right over this aperture, this figure shows the amount of times a certain actuator was deemed to be in the set of most informative actuators. The actuators located at the edge where the new turbulence comes in are much more regularly chosen than the other actuators. This is expected, since the aberration towards the center and right side of the aperture is merely a shift of the wavefront at a previous time instance. The aberration at the left side of the aperture was unknown at this previous time instance and more difficult to model.

#### C. Stronger Aberrations

As the aberrations become more severe, the PSF images become more distorted. Therefore, a decrease in performance is expected for any algorithm. If the same performance for a lower Fried parameter is desired, the number of actuators has to be increased. Figure 7 shows the performance of the method for ${r_0}$ varying between 0.1 m and 0.35 m, while all other parameters are fixed at their standard value given in Table 2. The top line indicated by “DM optimal” corresponds to a controller that assumes perfect knowledge of the residual wavefront and maps it onto the DM. So it represents the theoretical optimal performance when using this DM. When the turbulence strength increases, the performance of both methods decreases as expected. Especially for turbulence strengths ${r_0}\ge 0.2$, the new method clearly outperforms the existing method.

#### D. Noise Sensitivity

As discussed before, the influence of camera noise on the SM will be approximately Gaussian. Therefore, the Gaussian measurement noise is added to the output $\boldsymbol{y}(kT)$ in order to simulate noisy conditions. In these results, the measurement noise variance is supposed to be known. In practice, this must be calibrated based on the exact noise properties of the camera and PSF intensity. The influence of the measurement noise on the results is shown in Fig. 8. For low noise, i.e., ${\sigma_y}\le 0.1$, the existing method has an advantage. When there is a significant measurement noise, which will be the case in any practical system, the new method is clearly better. This is expected, since the Kalman filter is designed to optimally deal with measurement noise, whereas the original method ignored any noise present in the system.

## 6. CONCLUSION

A method has been presented to extend the SM-based WFSless AO to the case of time-varying aberrations. It combines the knowledge of an accurate temporal model of the aberration dynamics with the linear relation between the SM of the PSF and the mean square of the residual phase aberration. The result is that the AO problem can be cast into a Kalman filtering and optimal control problem. Where all previous methods had to update the actuators sequentially, the new algorithm automatically selects the actuator that is expected to lead to the most informative update. Actuators placed at locations where the model is accurate can be updated with only a small number of images over time. The improved performance has been shown in a simulation study of an AO system for atmospheric turbulence. It was shown that the new method leads to an improved performance for both lower and higher wind speeds and for higher noise levels.

## Funding

H2020 European Research Council (ERC) (339681); Seventh Framework Programme (FP7) (2007-2013).

## REFERENCES

**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. **S. Zommer, E. Ribak, S. Lipson, and J. Adler, “Simulated annealing in ocular adaptive optics,” Opt. Lett. **31**, 939–941 (2006). [CrossRef]

**3. **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]

**4. **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]

**5. **C. E. Carrizo, R. M. Calvo, and A. Belmonte, “Intensity-based adaptive optics with sequential optimization for laser communications,” Opt. Express **26**, 16044–16053 (2018). [CrossRef]

**6. **X. He, X. Zhao, S. Cui, and H. Gu, “A rapid hybrid wave front correction algorithm for sensor-less adaptive optics in free space optical communication,” Opt. Commun. **429**, 127–137 (2018). [CrossRef]

**7. **Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, “Phase retrieval with application to optical imaging: a contemporary overview,” IEEE Signal Process. Mag. **32**, 87–109 (2015). [CrossRef]

**8. **R. Doelman, N. H. Thao, and M. Verhaegen, “Solving large-scale general phase retrieval problems via a sequence of convex relaxations,” J. Opt. Soc. Am. A **35**, 1410–1419 (2018). [CrossRef]

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

**10. **H. Yang, O. Soloviev, and M. Verhaegen, “Model-based wavefront sensorless adaptive optics system for large aberrations and extended objects,” Opt. Express **23**, 24587–24601 (2015). [CrossRef]

**11. **W. Lianghua, P. Yang, Y. Kangjian, C. Shanqiu, W. Shuai, L. Wenjing, and B. Xu, “Synchronous model-based approach for wavefront sensorless adaptive optics system,” Opt. Express **25**, 20584–20597 (2017). [CrossRef]

**12. **C. Kulcsár, H.-F. Raynaud, C. Petit, J.-M. Conan, and P. V. De Lesegno, “Optimal control, observers and integrators in adaptive optics,” Opt. Express **14**, 7464–7476 (2006). [CrossRef]

**13. **K. Hinnen, M. Verhaegen, and N. Doelman, “A data-driven H_{2}-optimal control approach for adaptive optics,” IEEE Trans. Control Syst. Technol. **16**, 381–395 (2008). [CrossRef]

**14. **W. Lianghua, P. Yang, W. Shuai, L. Wenjing, C. Shanqiu, and B. Xu, “A high speed model-based approach for wavefront sensorless adaptive optics systems,” Opt. Laser Technol. **99**, 124–132 (2018). [CrossRef]

**15. **O. Soloviev, “Optimal basis for modal sensorless adaptive optics,” arXiv:1707.08489 (2017).

**16. **M. Verhaegen and V. Verdult, *Filtering and System Identification: A Least Squares Approach* (Cambridge University, 2007).

**17. **R. Conan and C. Correia, “Object-oriented MATLAB adaptive optics toolbox,” Proc. SPIE **9148**, 91486C (2014). [CrossRef]