Fundamental performance of transverse wind estimator from Shack-Hartmann wave-front sensor measurements

Zhenghan Li; Xinyang Li

doi:10.1364/OE.26.011859

1. Introduction

Adaptive optics (AO) system compensate for wave-front distortions due to atmospheric turbulence by using deformable mirror (DM) [1]. The Shack-Hartmann wave-front sensor (SHWS) detects the distorted wave-front due to atmospheric turbulence and outputs the distorted wave-front information by slope measurements. The control system calculates voltage commands to DM using slope measurements. DM generates corresponding shape that is conjugate to the distorted wave-front based on the voltage commands, partially compensating the distortions.

With the assumption that wave-front state is temporally static, traditional adaptive optics control systems suffer from 2 to 3 sampling periods servo lag [2]. The delay time comes from reading out the wave-front sensor and computing control commands. A great amount of effort has been invested into developing predictive correction methods to compensate the time delay error [3–5]. To solve the problems of heavy calculation cost in traditional predictive correction methods, predictive correction methods using transverse wind information are becoming a hot research area [6,7].

Under frozen-flow assumption, the atmospheric turbulence is located in several thin horizontal layers. The turbulence spatial pattern does not change dramatically within the encountered time delays, but slides according to the transverse wind. Estimating the frozen layers and individual transverse wind will remarkably reduce the computational tasks of predictive correction of AO. We proposed an open loop predictive correction method, see Fig. 1. And the key component of the method is the transverse wind estimation.

Fig. 1 Principle of AO predictive correction.

Download Full Size | PDF

Application of SHWS to wind profile measurement is not new. In addition to benefit predictive correction, transverse wind measurement is conducive to investigate atmospheric properties and statistics. Considerable work has been reported on methods to the measurement [8–12]. Unfortunately, the standards of performance evaluation of these estimators are not uniform. The main purpose of this work is to find the performance limits in transverse wind estimation using Cramér–Rao lower bound (CRLB). Developing such performance bounds provides a fundamental benchmark to compare the performance of these estimators. In addition, the analysis of the details of the bound provides an insight into the nature of the problem itself, thereby suggesting how to improve the performance of estimator. A preliminary work on performance analysis on wind estimator is given in [12], but did not consider the relationship between performance and atmospheric parameters such as $r_{0}$ and wind velocity. Another work on performance analysis on wind estimator is given in [6] and has the following defects: first of all, the analysis did not consider the bias in the wind estimator; second, the derived limits cannot predict the performance of estimator accurately in high signal-to-noise regimes. Although the estimation bias is often difficult to express, we will derive the bias expression of gradient-based transverse wind estimator and analyze the source of the bias for the first time. Finally, we will derive the biased CRLB for the gradient-based transverse wind estimator which can predict the performance of estimator more accurately.

2. Models

2.1 Shack-Hartmann wave-front sensor

Shack-Hartmann wave-front sensor (SHWS) is commonly used in adaptive optics systems. It consists of an array of lenses of the same focal length, see Fig. 2.

Fig. 2 SHWS subapertures layout ( $14 \times 14$ ).

Download Full Size | PDF

Each lens is focused onto a CCD camera. The local tilt of the wave-front across each lens can then be calculated from the position of the focal spot on the sensor, see Fig. 3.

Fig. 3 Principle of SHWS.

Download Full Size | PDF

Any phase aberration can be approximated by a set of discrete tilts. By sampling an array of lens, all of these tilts can be measured and the whole wave-front can be approximated.

\hat{s} (z, t) = H ϕ (z, t) + s_{n} (z, t)

where

ϕ (z, t)

is the disturbed wave-front phase,

\hat{s} (z, t)

is the slope measurements, and

s_{n} (z, t)

is the noise on the slope measurements. And

z

represents the spatial coordinates, with

x

and

y

two dimensions. In the open loop predictive correction scheme, see Fig. 1, the SHWS is located before the corrector and detects the aberration in real time. We assume the slope noise to be Gaussian noise that is spatially and temporally uncorrelated [6,13–15]. The covariance of the noise is

σ^{2} I

, and

I

is the identity matrix. And

H

is the operator of slope observation that transform the distorted wave-front into slope measurements.

2.2 Frozen-flow turbulence assumption

This work is mainly focusing on an atmospheric model that is driven by frozen-flow. Actual temporal dynamics of the phase aberrations involve other processes such as the static effects due to telescope inherent geometric abnormalities. However, frozen-flow is expected to be the most important component of the aberration dynamics on a time scale of milliseconds. Previous works describe how to identify and separate static phase aberrations from frozen-flow turbulence using a recursive mean estimator [6]. In this article, it is assumed that all phase aberrations not due to frozen-flow turbulence have been removed.

Frozen-flow turbulence assumption is a commonly used modeling technique that has been verified through astronomical observations [16–19]. These studies show that the assumption is generally valid on a time scale of 10-20 milliseconds and is occasionally on a time scale of 50-100 milliseconds. In this paper, the sampling period of the AO system is 1 milliseconds and the time delay is three sampling periods, so the frozen-flow turbulence assumption is reasonable. Under frozen-flow turbulence assumption, the atmospheric turbulence is located in several thin horizontal layers. For a single layer of frozen-flow turbulence, the temporal dynamics of the phase aberrations are translational shift. The turbulence spatial pattern do not change dramatically within the encountered time delays, but slide according to the transverse wind.

ϕ (z, t) = ϕ (z - w, t - 1)

where

ϕ (z, t)

is the disturbed wave-front phase at spatial position

z

, and

t

is the discrete time index. Transverse wind vector

w

has

u

and

v

two dimensions,

w = {(u, v)}^{T}

. Without the consideration of noise, we can derive the similar temporal dynamics for the slope measurements using Eq. (1):

s (z, t) = s (z - w, t - 1)

2.3 Gradient-based transverse wind estimator

We introduced the non-iterative gradient-based transverse wind estimator. We assume that the measurement data is in this form:

\hat{s} (z, t) = s (z, t) + s_{n} (z, t)

\hat{s} (z, t + 1) = s (z + w, t) + s_{n} (z, t + 1)

We calculate the time-domain difference of two measurement data as:

s_{t} (z) = s (z + w, t) - s (z, t) + ε

where

s_{t}

represents time-domain difference, and

ε

is a Gaussian white noise. We linearize

s (z + w, t)

about a point

w = {(0, 0)}^{T}

in a Taylor series as:

s (z + w, t) - s (z, t) = s_{z} w + R (z, w)

where

s_{z}

is spatial gradient

[s_{x}, s_{y}]

. And

R (z, w)

represents the remainder term in the Taylor expansion. With the remainder term

R (z, w)

ignored, Eq. (7) becomes:

s_{t} = s_{z} w + ε

Using the approximated gradients, we can calculate the transverse wind using least squares:

\hat{w} = \sum_{x, y} \frac{s_{z}^{T} s_{t}}{s_{z} s_{z}^{T}}

where the sum is taken to be over all slope measurements.

2.4 Performance bounds in transverse wind estimation with SHWS

We proposal to use mean square error (MSE) to measure the performance of the transverse wind estimator. The $MSE (w)$ is the diagonal terms of $E [(\hat{w} - w) {(\hat{w} - w)}^{T}]$ . We explore the fundamental performance bound for transverse wind estimator using CRLB [20]. CRLB provides a lower limit for the MSE of any methods used to estimate a deterministic parameter vector from a given set of data. The CRLB on transverse wind $w$ for any estimator is given by:

MSE (w) \geq \frac{\partial E [\hat{w}]}{\partial w} J^{- 1} (w) \frac{\partial E {[\hat{w}]}^{T}}{\partial w} + (E [\hat{w}] - w) {(E [\hat{w}] - w)}^{T}

where the matrix

J (w)

is the Fisher information matrix (FIM).

\hat{w}

is the transverse wind estimation, and

E [\hat{w}]

represents the expectations of the estimation [21]. And the bias of the estimation is given by:

b (w) = E [\hat{w}] - w

For the unbiased estimator, Eq. (10) simplifies to the form as:

MSE (w) \geq J^{- 1} (w)

The bound expresses that the variance of any unbiased estimator is at least as high as the inverse of the Fisher information.

Since Fisher information plays an important role in constraining the variance of estimators, we will first explore the details of Fisher information in transverse wind estimation in section 3. Then we will show the existence of systematic bias in the gradient-based transverse wind estimator and analyze the source of the bias in section 4. And we will explore the details of the systematic bias in section 5. In the end, we will derive the biased CRLB which can predict the performance of estimator more accurately in section 6.

3. Fisher information for transverse wind estimation with SHWS

We use $s (z, t)$ represents ideal slope measurements without noise. The joint likelihood function for the $\hat{s} (z, t)$ and $\hat{s} (z, t - 1)$ is represented by $P (\hat{s}; w)$ . And the natural logarithm of the likelihood function is given by:

\ln P (\hat{s}; w) = \frac{- 1}{2 σ^{2}} \sum_{x, y} {[\hat{s} (z, t) - s (z, t)]}^{2} + {[\hat{s} (z, t) - s (z - w, t - 1)]}^{2} + const

Fisher information is defined as the expected curvature of the log likelihood:

J (w) = - E [\frac{\partial^{2} \ln P (\hat{s}; w)}{\partial w \partial w^{T}}]

The FIM in transverse wind estimation is a 2 by 2 matrix, and we define the matrix as:

J (w) = [\begin{matrix} a_{1} & a_{2} \\ a_{2} & a_{3} \end{matrix}]

The elements of the $J (w)$ is calculated as:

a_{1} = \sum_{x, y} - E [\frac{\partial^{2} \log P (\hat{s}; w)}{\partial u^{2}}] = \sum_{x, y} \frac{1}{σ^{2}} {(\frac{\partial s (z - w, t - 1)}{\partial u})}^{2}

a_{2} = \sum_{x, y} - E [\frac{\partial^{2} \log P (\hat{s}; w)}{\partial u \partial v}] = \sum_{x, y} \frac{1}{σ^{2}} (\frac{\partial s (z - w, t - 1)}{\partial u}) (\frac{\partial s (z - w, t - 1)}{\partial v})

a_{3} = \sum_{x, y} - E [\frac{\partial^{2} \log P (\hat{s}; w)}{\partial v^{2}}] = \sum_{x, y} \frac{1}{σ^{2}} {(\frac{\partial s (z - w, t - 1)}{\partial v})}^{2}

Noting that by way of chain rule, we can write Eq. (16)-(18) as:

a_{1} = \sum_{x, y} \frac{1}{σ^{2}} {(\frac{\partial s (z, t)}{\partial x})}^{2} = \sum_{x, y} \frac{1}{σ^{2}} s_{x}^{2} (z, t)

a_{2} = \sum_{x, y} \frac{1}{σ^{2}} (\frac{\partial s (z, t)}{\partial x}) (\frac{\partial s (z, t)}{\partial y}) = \sum_{x, y} \frac{1}{σ^{2}} s_{x} (z, t) s_{y} (z, t)

a_{3} = \sum_{x, y} \frac{1}{σ^{2}} {(\frac{\partial s (z, t)}{\partial y})}^{2} = \sum_{x, y} \frac{1}{σ^{2}} s_{y}^{2} (z, t)

where

s_{x}

and

s_{y}

represent the partial derivative in the

x

and

y

directions respectively. In practice, the partial derivatives of the ideal slope are not available, since it is hard to obtain the formal slope distributions. In most applications, only the discrete slope measurements are available. And we need to approximate the partial derivatives using the spatial gradient operator.

3.1 Fisher information and atmospheric coherence length

Although it is hard to obtain the formal slope distribution, we can analysis the problem theoretically with the approximation that the slope distribution is similar to wave-front distribution. In this article, we discuss the Kolmogorov turbulence. For further understanding the FIM, we consider the problem in the Fourier domain. We assume the discrete Fourier transform (DFT) of the slope $s (z)$ is given by:

s (z) = \sum_{f_{x}} \sum_{f_{y}} S (f) exp (j f^{T}) Δ f_{x} Δ f_{y}

where

S (f)

are the coefficients of the DFT of the slope

s (z)

. And

f

represents spatial frequency

f = {(f_{x}, f_{y})}^{T}

, where

f_{x}

and

f_{y}

are the spatial frequencies in the

x

and

y

directions. Therefore, the partial derivative

s_{x} (z, t)

can be written as

j f_{x} S (f) exp (j f^{T})

and similarly for the

y

partial derivative. Then, we can write Eq. (14) as:

J (w) = \frac{1}{σ^{2}} \sum_{f_{x}} \sum_{f_{y}} {| S (f) |}^{2} f f^{T} Δ f_{x} Δ f_{y}

When considering Kolmogorov turbulence, we can approximate $S (f)$ as

S (f) = h (f) \sqrt{0.023 r_{0}^{- 5 / 3} {‖ f ‖}^{- 11 / 3}}

where

h (f)

is Hermitian matrix with zero mean and unit variance, and

r_{0}

is the atmospheric coherence length.

Using Eq. (24) we can write Eq. (23) as:

J (w) = \frac{0.023}{σ^{2} r_{0}^{5 / 3}} \sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} f f^{T} Δ f_{x} Δ f_{y}

We can see that FIM is a function of $σ$ and $r_{0}$ . Considering the $MSE (w)$ is the diagonal terms of $E [(\hat{w} - w) {(\hat{w} - w)}^{T}]$ , we use the square root of the trace of the inverse FIM as the scalar bounds of performance. We can obtain the theoretical conclusions that the square root of the trace of the inverse FIM will be proportional to $r_{0}^{5 / 6}$ based on Eq. (25).

In order to verify the analysis results, it is necessary to obtain the atmospheric turbulence phase screen. The simulation conditions for generating the phase screen are shown in Table 1.

Table 1. Summary of parameters of simulation

View Table

The Kolmogorov phase screens are generated based on power spectrum method. To simulate the frozen flow slides according to the transverse wind, the Fourier domain method is used. The phase screen is translated by the specified horizontal wind ranging from 0 grid/frame to 22 grid/frame. Transverse wind velocity for simulation is 1 grid/frame, while the actual transverse wind velocity is 6 m/s based on the simulation conditions. The SHWS with 14 subapertures across the 1.2-meter primary aperture is used to measure the average wave-front slope over each of its subapertures. The subaperture length is about 8.6 cm. The resolution of the SHWS is 200 × 200, and there are 14 pixels across each subaperture.

We explore $r_{0}$ ranging from 10 cm to 30 cm for the subaperture size of 8.6 cm. Transverse wind velocity for simulation is 5 grid/frame (30m/s). The square root of the trace of the inverse FIM is computed for each of the $r_{0}$ using Eqs. (19)-(21). The partial derivatives are approximated by central difference operator, and the gradient operator will be discussed in detail later in this article.

As shown in Fig. 4, the red circles represent the numerical results of the square root of the trace of the inverse FIM in each value of $r_{0}$ . The dashed line represents the theoretical result that the unbiased CRLB is proportional to $r_{0}^{5 / 6}$ . We find that numerical results agree well with the theoretical analysis based on Eq. (25). The unbiased CRLB suggests improvement of performance as the turbulence increases. However, we should point out that the FIM depends on the inverse 5/3 power of $r_{0}$ is not general, but limited to the cases where $r_{0}$ not less than about 3/4 of the subaperture size. In practice, the subaperture size would affect the FIM. For $r_{0}$ less than about 3/4 subaperture length, decreasing $r_{0}$ would damage the performance.

Fig. 4 RMSE as it relates to $r_{0}$ .

Download Full Size | PDF

We find that CRLB is inversely proportional to the total number of subapertures in the system. Because the high-order SHWS can accurately detect more high frequency information which increase the performance of the transverse wind estimator. The conclusion is consistent with the observation that estimator performance will improve as the turbulence increases for $r_{0}$ not less than about 3/4 of the subaperture size. It suggests that high-order SHWS would get better transverse wind estimation performance.

3.2 Fisher information and slope measurement noise

We find that slope measurement noise plays an important role in the Fisher information Eq. (25). We can obtain the theoretical conclusions that the square root of the trace of the inverse FIM will be proportional to the standard deviation of the slope measurement noise based on Eq. (25). We explore the standard deviation of the slope measurement noise $σ$ ranging from 0.01 pixel to 1 pixel. The simulation conditions are the same as 3.1, but we fixed $r_{0}$ of 15 cm. The square root of the trace of the inverse FIM is computed for each of the noise condition using Eqs. (19)-(21). The partial derivatives are approximated by central difference operator.

As shown in Fig. 5, the red circles represent the numerical results of the square root of the trace of the inverse FIM in each value of the slope measurement noise $σ$ . The dashed line represents the theoretical result that the unbiased CRLB is proportional to $σ$ . We find that numerical results agree well with the theoretical analysis drawn from Eq. (25). The unbiased CRLB suggests improvement of performance in better noise condition.

Fig. 5 RMSE as it relates to slope noise.

Download Full Size | PDF

We find that SHWS would get better transverse wind estimation performance with brighter SHWS reference. Because the brighter SHWS reference decrease the slope measurement noise. The conclusion is consistent with the observation that estimator performance will improve in better noise condition.

3.3 Fisher information and wind direction

Now we explore the relationship between Fisher information and wind direction. Assuming 0 degree is parallel to the $x$ axis of the SHWS grid and wind direction is $φ$ , we express the transverse wind in the following linear combination form:

w_{φ} = d^{T} w = u \cos (φ) + v \sin (φ)

where

d

is

{[\cos (φ), \sin (φ)]}^{T}

. And we can derive the unbiased CRLB for

w_{φ}

as:

MSE (w_{φ}) \geq \frac{σ^{2}}{a_{1} a_{3} - a_{2}^{2}} [a_{3} \cos^{2} (φ) + a_{1} \sin^{2} (φ) - 2 a_{2} \sin (φ) \cos (φ)]

We can obtain the theoretical conclusions that the unbiased CRLB will be as a function of wind direction $φ$ based on Eq. (27). We explore the wind direction $φ$ ranging from 0 degree to 360 degree. The simulation conditions are the same as 3.1, but we fixed $σ$ of 0.1 pixel. The $a_{1}$ , $a_{2}$ and $a_{3}$ are computed using Eqs. (19)-(21). The partial derivatives are approximated by central difference operator. As shown in Fig. 6, the unbiased CRLB will be slightly different in different wind direction, which agrees well with the theoretical analysis based on Eq. (27). We can use derivative to get the extreme value for RMSE and find that the extreme value will be obtained at the wind direction calculated by $φ = \arcsin (2 a_{2} / (a_{1} - a_{3}))$ . It suggests that the performances of the transverse wind estimator will come to an extreme value in these directions.

Fig. 6 RMSE as it relates to wind direction: (a) rectangular coordinate; (b) polar coordinate.

Download Full Size | PDF

4. Bias in transverse wind estimation with SHWS

In this section, we will show the existence of the deterministic bias in the estimator and explore the source of the deterministic bias.

4.1 Bias in gradient-based transverse wind estimator

By comparing the actual MSE performance with unbiased CRLB, we show the existence of bias in the gradient-based transverse wind estimator using slope measurements of SHWS. White Gaussian noise was added to the slope measurements prior to estimation and the entire process was repeated 500 times at each noise condition. We explore the standard deviation of the slope measurement noise $σ$ ranging from 0.01 pixel to 1 pixel. We obtain the actual performance by computing the square root of the trace of the MSE matrix. The spatial gradients are approximated by central difference operator. Figure 7 shows the actual estimator performance as a function of slope measurement noise. The dashed line indicates the performance limit for the unbiased estimator. Although this limit indicates a constant improvement as the noise decreases, but the performance of the estimator flattens below certain $σ$ value. The ﬂattening of the actual performance curve indicates that there is a bias in the estimator. It means that the bound given by Eq. (12) is imprecise and that the complete bound Eq. (10) must be used to predict the performance of the estimator more accurately.

Fig. 7 RMSE versus slope noise.

Download Full Size | PDF

4.2 Source of deterministic bias

As mentioned in section 2.3, one source of the deterministic bias is from the ignorance of the remainder term $R (z, w)$ . The linear approximation error of the gradient-based transverse wind estimator will increase as the wind velocity increase. Another source of the deterministic bias is from the gradient approximation. For further understanding the estimator bias, we consider the estimation process in the Fourier domain. The DFT of $s (z + w, t)$ is $S (f) [\exp (j f^{T} w) - 1]$ , and Eq. (6) can be written as:

S_{t} (f) = S (f) [\exp (j f^{T} w) - 1] + Z (f)

where

Z (f)

is the Gaussian white noise in the Fourier domain. In practice, we calculate the transverse wind using the approximated gradients. The DFT of the approximated gradients is

| S (f) | G (f)

, where

G (f)

is the gradient operator in the Fourier domain. We obtain the estimator as

\hat{w} = Q^{- 1} \sum_{f_{x}} \sum_{f_{y}} {| S (f) |}^{2} j G (f) S_{t}^{*} (f) Δ f_{x} Δ f_{y}

where

Q = \sum_{f x} \sum_{f x} {| S (f) |}^{2} G (f) G {(f)}^{T} Δ f_{x} Δ f_{y}

.

The expected value of the estimation result is

E [\hat{w}] = Q^{- 1} \sum_{f_{x}} \sum_{f_{y}} {| S (f) |}^{2} G (f) \sin (f^{T} w) Δ f_{x} Δ f_{y}

We can obtain the estimator bias as:

b (w) = Q^{- 1} \sum_{f_{x}} \sum_{f_{y}} {| S (f) |}^{2} [G (f) \sin (f^{T} w) - G (f) G {(f)}^{T} w] Δ f_{x} Δ f_{y}

Using Eq. (24) we can write Eq. (31) as:

b (w) = \frac{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} [G (f) \sin (f^{T} w) - G (f) G {(f)}^{T} w] Δ f_{x} Δ f_{y}}{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} G (f) G {(f)}^{T} Δ f_{x} Δ f_{y}}

We can obtain the theoretical conclusions that the estimator bias is a function of $w$ and the bias is dependent on gradient operator and independent on $r_{0}$ drawn from Eq. (32).

In practice, we can calculate the estimator bias as:

b (w) = \sum_{x, y} \frac{s_{z}^{T} s_{t}}{s_{z} s_{z}^{T}} - w

The estimator bias is calculated using noiseless slope measurement. We explore the transverse wind velocity ranging from 0 grid/frame to 8 grid/frame (0 subaperture/frame to 0.56 subaperture/frame). We explore $r_{0}$ ranging from 10 cm to 25 cm (1.2 subaperture to 2.9 subaperture). There is no temporal aliasing of the gradient sampling in the range of $r_{0}$ and wind velocity. The estimator bias is calculated by Eq. (33). The spatial gradients are approximated by central difference operator. As shown in Fig. 8, the estimator bias is a function of wind velocity, which agrees well with the theoretical analysis drawn from Eq. (32). Consistent with the theoretical analysis drawn from Eq. (32), the simulation results demonstrate that the estimator bias is independent on $r_{0}$ . However, we should point out that the finding of the relationship between bias and $r_{0}$ is not general, but limited to the cases where $r_{0}$ not less than about 3/4 of the subaperture size and there is no temporal aliasing of the gradient sampling.

Fig. 8 Bias versus wind velocity in different $r_{0}$ .

Download Full Size | PDF

5. Analysis of gradient-based estimator bias

In this section, we further explore the deterministic bias approximation Eq. (31). We will explore how the transverse wind velocity and gradient operator affect the bias of the gradient-based transverse wind estimator.

In order to gain a deeper understanding of the bias, we Taylor expand the $\sin$ function in Eq. (31). Considering the actual wind velocity is smaller than 1 subaperture/frame (84 m/s) under most scenarios, the high order terms will add with little effect on the bias function. The bias Eq. (31) can be simplified by truncating the series to that of a cubic function of:

b (w) = Q^{- 1} \sum_{f_{x}} \sum_{f_{y}} {| S (f) |}^{2} G (f) [f^{T} w - \frac{1}{6} {(f^{T} w)}^{3}] Δ f_{x} Δ f_{y} - w

Noting that $f^{T} w = ‖ w ‖ f^{T} n$ , where $n$ is the unit vector ${[\cos ψ, \sin ψ]}^{T}$ , we can write Eq. (34) as:

\begin{array}{l} b (w) = ‖ w ‖ Q^{- 1} \sum_{f_{x}} \sum_{f_{y}} {| S (f) |}^{2} G (f) [f^{T} - G (f)] n Δ f_{x} Δ f_{y} \\ - \frac{1}{6} {‖ w ‖}^{3} Q^{- 1} \sum_{f_{x}} \sum_{f_{y}} {| S (f) |}^{2} G (f) (^{f^{T} n) 3} Δ f_{x} Δ f_{y} \end{array}

Using Eq. (24) we can write Eq. (35) as

\begin{array}{l} b (w) = ‖ w ‖ \frac{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} G (f) [f^{T} - G (f)] n Δ f_{x} Δ f_{y}}{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} G (f) G {(f)}^{T} Δ f_{x} Δ f_{y}} \\ - \frac{1}{6} {‖ w ‖}^{3} \frac{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- \frac{11}{3}} G (f) {(f^{T} n)}^{3} Δ f_{x} Δ f_{y}}{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- \frac{11}{3}} G (f) G {(f)}^{T} Δ f_{x} Δ f_{y}} \end{array}

5.1 Bias and gradient operator

As Eq. (36) indicates, gradient operator plays an important role in the bias function. Gradient operator controls the coefficients of the bias function. We calculate the bias of the gradient-based transverse wind estimator using three different gradient operators: central difference operator, Prewitt operator and Sobel operator. The kernels of these gradient operators are shown in Fig. 9.

Fig. 9 The gradient kernel: (a)-(b) central difference operator; (c)-(d) Prewitt operator; (e)-(f) Sobel operator.

Download Full Size | PDF

The estimator bias is calculated using noiseless slope measurement. We explore the transverse wind velocity ranging from 0 subaperture/frame to 0.8 subaperture/frame. Because there is no temporal aliasing of the gradient sampling for the fixed $r_{0}$ of 15 cm (1.75 subaperture). Figure 10 shows the bias of the gradient-based transverse wind estimator using three different gradient operators.

Fig. 10 Bias versus wind velocity using different gradient operator.

Download Full Size | PDF

When the wind velocity is small, the linear term in the bias Eq. (36) will greatly affect the bias function. The performance of Prewitt operator and Sobel operator are poor because of their larger linear coefficient. The central difference operator minimizes the coefficient of the linear term and seems to be the optimal gradient operator in the velocity range.

5.2 Bias and wind velocity

The bias of the gradient-based transverse wind estimator will be divided into two regions, positive and negative, by examining the bias function Eq. (36).

‖ w^{*} ‖ = \pm {(\frac{6 \sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} G (f) [f^{T} - G (f)] n Δ f_{x} Δ f_{y}}{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} G (f) {(f^{T} n)}^{3} Δ f_{x} Δ f_{y}})}^{\frac{1}{2}}

where

‖ w^{*} ‖

is the transverse wind velocity that makes the bias is 0. We can find that the relative bias varies linearly with square of transverse wind velocity by examining the relative bias:

\begin{array}{l} η = \frac{b (w)}{w} = \frac{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} G (f) [f^{T} - G (f)] n Δ f_{x} Δ f_{y}}{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} G (f) G {(f)}^{T} Δ f_{x} Δ f_{y}} \\ - \frac{1}{6} {‖ w ‖}^{2} \frac{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} G (f) {(f^{T} n)}^{3} Δ f_{x} Δ f_{y}}{\sum_{f_{x}} \sum_{f_{y}} {| h (f) |}^{2} {‖ f ‖}^{- 11 / 3} G (f) G {(f)}^{T} Δ f_{x} Δ f_{y}} \end{array}

We explore the transverse wind velocity ranging from 0 subaperture/frame to 1.6 subaperture/frame. We should point out that there will be temporal aliasing of the gradient sampling when the transverse wind velocity larger than 0.875 subaperture/frame for the fixed $r_{0}$ of 15 cm (1.75 subaperture). Examination of the absolute bias in Fig. 11, we find there is a critical velocity. When the transverse wind velocity is faster than the critical velocity, the bias tends to worsen dramatically.

Fig. 11 Absolute bias versus wind velocity using different gradient operator.

Download Full Size | PDF

As shown in Fig. 11, the critical transverse wind velocities are around 0.9 subaperture/frame for the central difference operator, 1.3 subaperture/frame for the Prewitt operator and 1.5 subaperture/frame for Sobel operator. The finding of rapidly-increasing bias for wind velocities in the range 0.9-1.5 subaperture/frame is consistent with the onset of temporal aliasing of the gradient sampling. This finding is understandable and is a comprehensive result of the linear approximation, gradient approximation, and temporal aliasing of the gradient sampling.

Therefore, the gradient operator should be designed such that the transverse wind to be estimated is less than the critical velocity. Considering that the transverse wind velocity is less than 1 subaperture/frame (84 m/s) in most realistic scenarios, the central difference is the optimal gradient operator.

6. MSE performance of the gradient-based estimator

In this section, we will derive the biased CRLB Eq. (10) for the gradient-based transverse wind estimator. We will explore how the transverse wind velocity, slope noise, $r_{0}$ ,the number of subaperures and SHWS reference brightness affect the MSE performance.

First, we calculate the derivative of the expectation of transverse wind

\frac{\partial E [\hat{w}]}{\partial w^{T}} = \sum_{x, y} \frac{s_{z} s_{z} {(z + w)}^{T}}{s_{z} s_{z}^{T}}

Using Eq. (39), we obtain the biased CRLB as

\begin{array}{l} MSE (w) \geq (\sum_{x, y} \frac{s_{z} s_{z} {(z + w)}^{T}}{s_{z} s_{z}^{T}}) J^{- 1} (w) {(\sum_{x, y} \frac{s_{z} s_{z} {(z + w)}^{T}}{s_{z} s_{z}^{T}})}^{T} \\ + (\sum_{x, y} \frac{s_{z}^{T} s_{t}}{s_{z} s_{z}^{T}} - w) {(\sum_{x, y} \frac{s_{z}^{T} s_{t}}{s_{z} s_{z}^{T}} - w)}^{T} \end{array}

To verify the biased CRLB Eq. (40), White Gaussian noise was added to the slope measurements prior to estimation and the estimation process was repeated 500 times at each noise condition. We explore the standard deviation of the slope measurement noise $σ$ ranging from 0.01 pixel to 1 pixel. We contain the actual performance by computing the square root of the trace of the MSE matrix. The spatial gradients are approximated by central difference operator. As shown in Fig. 12, we compared the actual performance with both the unbiased CRLB and the biased CRLB. When the slope noise is small, the actual performance is almost the same as the biased CRLB and the RMSE is largely decided by estimation bias. The biased CRLB can predict MSE performance of the gradient-based transverse wind estimator more accurately. As suggested in Fig. 11, the bias of estimator using central difference operator will be less for 3 grid/frame (0.21 subaperture/frame) than for 6 grid/frame (0.42 subaperture/frame), therefore the MSE performance will be better for smaller transverse wind velocity. The numerical result shown in Fig. 12 is consistent with the conclusion.

Fig. 12 RMSE versus slope noise as it relates to wind velocity.

Download Full Size | PDF

As illustrated in the previous section, the bias is independent on $r_{0}$ when $r_{0}$ not less than about 3/4 of the subaperture size and there is no temporal aliasing for the encountered wind velocity. For $r_{0}$ not less than about 3/4 of the subaperture size, the unbiased CRLB is proportional to $r_{0}^{5 / 6}$ . We can obtain the theoretical conclusions that the biased CRLB will be a function of $r_{0}$ based on Eq. (40). To verify the conclusion, we explore $r_{0}$ ranging from 10 cm to 25 cm (1.2 subaperture to 2.9 subaperture) and the wind is fixed on 5 grid/frame (0.353 subaperture/frame). There is no temporal aliasing in the simulation condition. The standard deviation of the slope measurement noise ranging from 0.01 pixel to 1 pixel. The spatial gradients are approximated by central difference operator. As shown in Fig. 13, the MSE performance will be better for smaller $r_{0}$ when $r_{0}$ not less than about 3/4 of the subaperture size, which is consistent with our theoretical analysis drawn from Eq. (40). It suggests that we can get better estimation result when turbulence increase when $r_{0}$ not less than about 3/4 of the subaperture size and there is no temporal aliasing.

Fig. 13 RMSE versus slope noise as it relates to $r_{0}$ .

Download Full Size | PDF

As mentioned in previous, SHWS would get better transverse wind estimation performance with brighter SHWS reference, larger total number of subapertures in the system.

7. Conclusions and future work

This paper derived the fundamental performance limits for the transverse wind estimator from slope measurements of SHWS using the Cramér–Rao bound. We derived and analyzed the unbiased CRLB. The first finding is that the unbiased CRLB is proportional to $r_{0}^{5 / 6}$ for $r_{0}$ not less than about 3/4 of the subaperture size. It suggests that the MSE performance of the estimator will improve as the turbulence increase when $r_{0}$ not less than about 3/4 of the subaperture size. As we expected, the unbiased CRLB is proportional to the standard deviation of the slope measurement noise. The transverse wind estimator would get better performance with brighter SHWS reference, larger total number of subapertures in the system.

Then we introduced the non-iterative gradient-based transverse wind estimator. The source of the deterministic bias of the gradient-based transverse wind estimators is analyzed for the first time. One source of the deterministic bias is from the ignorance of the remainder term. The linear approximation error of the gradient-based transverse wind estimator will increase as the wind velocity increases. Another source of the deterministic bias is from the gradient approximation. The bias is independent on atmospheric coherence length $r_{0}$ when $r_{0}$ not less than about 3/4 of the subaperture size and there is no temporal aliasing. As a comprehensive result of the linear approximation, gradient approximation, and temporal aliasing of the gradient sampling, there is a critical velocity. When the transverse wind velocity is larger than the critical velocity, the bias tends to worsen dramatically.

Finally, we derived biased CRLB for the gradient-based transverse wind estimators from Shack-Hartmann wave-front sensor measurements and the bound can predict the performance of estimator more accurately. The biased CRLB is a function of wind velocity, $r_{0}$ , subaperture size, SHWS read noise and SHWS reference brightness. For $r_{0}$ less than about 3/4 of the subaperture size, decreasing $r_{0}$ would no longer increase the performance of the estimator. Further, because decreasing $r_{0}$ requires smaller subapertures, a brighter SHWS reference is required. As power increases in higher spatial frequencies, the coherence time of the atmosphere decreases. Faster sampling generally requires a brighter SHWS reference.

As future work, we would like to introduce the multilayer transverse wind estimator and analyze the MSE performance for the multilayer estimator.

Funding

National Natural Science Foundation of China (61675205);

Acknowledgments

We gratefully acknowledge the comments and suggestions of Linhai Huang, Institute of Optics and Electronics, Chinese Academy of Sciences.

References and links

1. J. Hardy, Adaptive Optics for Astronomical Telescopes (Oxford University, 1998)

2. W. Jiang and H. Li, “Hartmann-Shack wave-front sensing and control algorithm,” Proc. SPIE 1271, 82–93 (1990). [CrossRef]

3. C. Kulcsár, H.-F. Raynaud, C. Petit, and J.-M. Conan, “Minimum variance prediction and control for adaptive optics,” Automatica 48(9), 1939–1954 (2012). [CrossRef]

4. L. Poyneer and J.-P. Véran, “Predictive wavefront control for adaptive optics with arbitrary control loop delays,” J. Opt. Soc. Am. A 25(7), 1486–1496 (2008). [CrossRef] [PubMed]

5. R. Fraanje, J. Rice, M. Verhaegen, and N. Doelman, “Fast reconstruction and prediction of frozen flow turbulence based on structured Kalman filtering,” J. Opt. Soc. Am. A 27(11), A235–A245 (2010). [CrossRef] [PubMed]

6. L. C. Johnson, D. T. Gavel, and D. M. Wiberg, “Bulk wind estimation and prediction for adaptive optics control systems,” J. Opt. Soc. Am. A 28(8), 1566–1577 (2011). [CrossRef] [PubMed]

7. R. Juvénal, C. Kulcsár, H.-F. Raynaud, and J.-M. Conan, “LQG adaptive optics control with wind-dependent turbulent models,” Proc. SPIE 9909, 99090M (2016). [CrossRef]

8. M. Schöck and E. J. Spillar, “Measuring wind speeds and turbulence with a wave-front sensor,” Opt. Lett. 23(3), 150–152 (1998). [CrossRef] [PubMed]

9. K. Yuan, “Measurement of Path-Averaged Transverse Wind Speed with a Shack-Hartmann Wave-Front Sensor,” Acta Opt. Sin. 29(2), 303–307 (1998). [CrossRef]

10. M. B. Roopashree, V. Akondi, and P. B. Raghavendra, “Real-time wind speed measurement using wave-front sensor data,” Proc. SPIE 7588, 75880A (2010). [CrossRef]

11. S. Abado, S. Gordeyev, and E. J. Jumper, “Approach for two-dimensional velocity mapping,” Opt. Eng. 52(7), 071402 (2012). [CrossRef]

12. L. C. Johnson, D. T. Gavel, and D. M. Wiberg, “Bulk Wind Estimator Performance for AO Systems,” in Frontiers in Optics 2009/Laser Science XXV/Fall 2009 OSA Optics & Photonics Technical Digest, OSA Technical Digest (CD) (Optical Society of America, 2009), paper AOWB5.

13. C. Plantet, S. Meimon, J.-M. Conan, and T. Fusco, “Revisiting the comparison between the Shack-Hartmann and the pyramid wavefront sensors via the Fisher information matrix,” Opt. Express 23(22), 28619–28633 (2015). [CrossRef] [PubMed]

14. F. Rigaut and E. Gendron, “Laser guide star in adaptive optics - the tilt determination problem,” Astron. Astrophys. 261, 677–684 (1992).

15. G. Rousset, “Wave-front Sensors,” in Adaptive Optics in Astronomy, F. Roddier, ed. (Cambridge University., 1999), pp. 91–130.

16. M. Schöck and E. J. Spillar, “Method for a quantitative investigation of the frozen flow hypothesis,” J. Opt. Soc. Am. A 17(9), 1650–1658 (2000). [CrossRef] [PubMed]

17. L. Poyneer, M. van Dam, and J.-P. Véran, “Experimental verification of the frozen flow atmospheric turbulence assumption with use of astronomical adaptive optics telemetry,” J. Opt. Soc. Am. A 26(4), 833–846 (2009). [CrossRef] [PubMed]

18. E. Gendron and P. Lena, “Single layer atmospheric turbulence demonstrated by adaptive optics observations,” Astrophys. Space Sci. 239, 221–228 (1996). [CrossRef]

19. B. Kern, T. A. Laurence, C. Martin, and P. E. Dimotakis, “Temporal coherence of individual turbulent patterns in atmospheric seeing,” Appl. Opt. 39(27), 4879–4885 (2000). [CrossRef] [PubMed]

20. S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory (Prentice-Hall, 1993).

21. H. L. V. Trees, Detection, Estimation, and Modulation Theory (Wiley, 1968).

Parameter	Value
Telescope diameter	1.2 m
Obstruction ratio	0.2
Layers	1
Wavelength	550 nm
Height	0 km
Sampling period	1 ms

Fundamental performance of transverse wind estimator from Shack-Hartmann wave-front sensor measurements

Abstract

1. Introduction

2. Models

2.1 Shack-Hartmann wave-front sensor

2.2 Frozen-flow turbulence assumption

2.3 Gradient-based transverse wind estimator

2.4 Performance bounds in transverse wind estimation with SHWS

3. Fisher information for transverse wind estimation with SHWS

3.1 Fisher information and atmospheric coherence length

3.2 Fisher information and slope measurement noise

3.3 Fisher information and wind direction

4. Bias in transverse wind estimation with SHWS

4.1 Bias in gradient-based transverse wind estimator

4.2 Source of deterministic bias

5. Analysis of gradient-based estimator bias

5.1 Bias and gradient operator

5.2 Bias and wind velocity

6. MSE performance of the gradient-based estimator

7. Conclusions and future work

Funding

Acknowledgments

References and links

Cited By

Figures (13)

Tables (1)

Equations (40)

Optics Express