Statistical Interpolation Method of Turbulent Phase Screen

Han-Ling Wu; Hai-Xing Yan; Xin-Yang Li; Shu-Shan Li

doi:10.1364/OE.17.014649

1. Introduction

Multiple phase screen model has been widely used in the simulation of light propagation through the turbulent atmosphere, in the design of adaptive optics (AO) systems and free space optical communication links [1–5]. A turbulent phase screen is a two-dimensional array of random phase values on a grid of sample points and represents phase perturbations of a propagating wave-front through the atmosphere in accordance with Kolmogorov theory [6,7]. Thanks to the insertion of a number of phase screens along the propagation path of an optical wave, the phase perturbations on a horizontal plane are introduced into the propagating optical field by the phase screen [8]. To be specific, the effects of atmospheric turbulence on optical fields are generally realized by adding random phase perturbations at the grid points of optical fields corresponding to those of phase screens. The relationship between the grid points of optical fields and those of phase screens is generally a one-to-one correspondence at spatial positions.

But a relative displacement between the grid points of optical fields and those of phase screens may occur in the simulation of many practical scenarios, especially in time-dependent simulation studies and the simulation of an astronomical telescopes with a laser guide star (LGS) adaptive optics system. The former is due to the movement of phase screens and the latter is caused by the mismatch between the grid points of LGS optical fields and those of phase screens. The relative displacement makes their grid points to be staggered when an optical field arrives at the positions of the phase screens. As a result, the phase perturbations of the optical field caused by atmosphere turbulence cannot directly be obtained from the grid points of phase screens. How to obtain the phase perturbations is an important problem in the simulation of these scenarios.

In this paper, we address the problem of the relative displacement between the grid points of optical fields and those of phase screens. A statistical interpolation approach to obtain phase perturbations is presented, based on the spatial phase structure function, which allows us to take into account the spatial statistics of the turbulent atmosphere. A comparison between the statistical interpolation method and the well-known linear interpolation method is made and analyzed. In order to test the validity of the statistical interpolation method, we also compare the phase structure function predicted by the base phase screens with the interpolated results. In addition, we carry out numerical experiments of light propagation through atmospheric turbulence without/with an adaptive optics system, i.e. open-loop/closed-loop propagation. The long-exposure Strehl ratio is used as the evaluation parameter in the numerical experiments. Emphasis is given to the effects of different interpolation methods on the evaluation parameter. These numerical results can be used to check the accuracy of the statistical interpolation method.

2. Problem statement

In this section, we describe how a relative displacement between the grid points of optical fields and those of phase screens occurs in the simulation of practical scenarios. Actually, these scenarios can be divided into two main categories: one-dimensional relative displacement and two-dimensional relative displacement.

2.1 One-dimensional relative displacement

In many time-dependent simulation studies, such as long exposure imaging in the ground-based astronomy [9] and the estimation of the dynamic control performance in an AO system [10], it is very important how to effectively generate a time series of atmospherically distorted wave-fronts.

According to the Taylor frozen turbulence hypothesis, a common approach for simulating temporally evolving wave-fronts is to produce a phase screen larger than the aperture area of interest and then to shift the screen to simulate the effect of wind blowing the turbulent media over the system pupil [9–12]. Although moving phase screens is a simple, fast and computer effective method for generating time series of atmospherically distorted wave-fronts, a disadvantage of this approach is that the screen may only be shifted by integer multiples of the grid spacing [11]. In general, the movement distances of phase screens due to the wind blowing should be determined by the wind speed and the response time of an optical system and they are not always integer multiples of the grid spacing. For example, in the simulation of a practical astronomical AO system [5], the movement velocity of each phase screen along the propagation path is different because the wind speed usually varies with height. Thus, the movement distance of each phase screen is also different at each time step. This will cause some of phase screens to shift by non-integer multiples of the grid spacing. Generally speaking, the phase screen almost always shifts by non-integer multiples of the grid spacing in many practical scenarios. A similar problem can be met with as well under the condition that the light source and/or receiver are in motion. Only in some particular cases, the screen may shift by integer multiples of the grid spacing [4].

In order to avoid this shortcoming, other methods have been proposed for time series of wave-front. One of them is based on the space-time covariance matrix of the phases over the pupil [13,14]. Another method is realized by a set of Fourier-series-transform-based modal expansions of the wave-front over the pupil [11]. However, they are limited in that the major drawback of them is the high demand of computing time or memory space. It is noteworthy that a different method is recently presented by performing extrapolations of successive columns from an initially phase screen [15,16]. The main advantage is that it can produce an indefinitely long screen with optimized memory storage. Nevertheless, the method for generation of time series of the wave-fronts by moving phase screens is still attractive in practical applications. When the phase screen shifts by non-integer multiples of the grid spacing, a relative displacement between the grid points of optical fields and those of phase screens will occur. In these cases, the shift of phase screens by non-integer multiples of the grid spacing generally makes the grid points of the optical field and those of the phase screen staggered.

To be specific, the grid points of optical fields are in one-to-one correspondence with those of phase screens at the initial time. But, after the screen moves, the grid points of the optical field and those of the screen may be staggered. That is to say, a relative displacement occurs in this case. It is shown in Fig. 1 . At the initial time t = t₀, the phase perturbations of an optical field can be obtained from the corresponding grid points of a phase screen due to one-to-one correspondence between their grid points. See Fig. 1 (a). If the phase screen shifts by one grid spacing or integer multiples of the grid spacing, the grid points of the optical field are still in one-to-one correspondence with those of the phase screen at the time t = t₁. See Fig. 1 (b). But, if the screen shifts by non-integer multiples of the grid spacing, the grid points of the optical field and those of the screen are staggered at the time t = t₂. See Fig. 1 (c). All grid points of the optical field are among those of the screen. Their grid points are staggered in the direction of the wind. Obviously, the phase perturbations of the optical field induced by atmospheric turbulence cannot be directly obtained from the grid points of the phase screen in this case. This is a one-dimensional relative displacement problem. How to obtain the phase perturbations of the grid points of the optical field is an important problem in the simulation of these scenarios.

Fig. 1 Schematic illustration of the one-dimensional relative displacement. (a) At the initial time t = t₀, the grid points of an optical field are in one-to-one correspondence with those of a phase screen. (b) At the time t = t₁, the grid points of the optical field are still in one-to-one correspondence with those of the phase screen because the phase screen shifts by one grid spacing. (c) At the time t = t₂, a relative displacement between the grid points of the optical field and those of the phase screen occurs because the phase screen shifts by non-integer multiples of the grid spacing.

Download Full Size | PDF

2.2 Two-dimensional relative displacement

The two-dimensional relative displacement mainly occurs in the simulation of an astronomical telescope with a laser guide star (LGS) adaptive optics system.

The light rays from an astronomical object, which is at infinite or at a very great distance, generally form a parallel bundle and pass through the atmosphere as a cylinder, while an LGS emits a spherical wave-front and the rays from the LGS trace out a cone because of the finite height of LGS above the telescope [5]. The difference in the propagation paths (i.e. between the cylinder and the cone) leads to an incorrect sampling of the turbulence below the LGS. That is to say, the region of the light rays from the LGS is smaller than that of the astronomical object (see Fig. 2 (a) ). The grid points of the object optical field are in one-to-one correspondence with those of the phase screen. But, the grid points of the LGS optical field is among those of the phase screen and its grid spacing is smaller than that of the phase screen or the object optical field. A two-dimensional relative displacement between the grid points of the LGS optical field and those of the phase screen occur in the simulation of an astronomical telescope with a laser guide star (LGS) adaptive optics systems (see Fig. 2 (b)). Obviously, the phase perturbations of the object optical field can be obtained from the grid points of the phase screen, while the phase perturbations of the LGS optical field cannot be directly obtained from the phase screen. This is a two-dimensional relative displacement problem. When the temporal response of the adaptive optics system must be considered, there also exists a one-dimensional relative displacement problem. How to obtain these phase perturbations of LGS optical field is an important problem in the simulation of an astronomical telescope with a laser guide star (LGS) adaptive optics system [5,17].

Fig. 2 Schematic illustration of the two-dimensional relative displacement. (a) The propagation geometry in the simulation of an astronomical telescope with a laser guide star adaptive optics system. (b) The grid points of a LGS optical field is among those of the phase screen and the grid spacing of the LGS optical field is smaller than that of the phase screen or the object optical field due to the finite altitude of LGS.

Download Full Size | PDF

2.3 Solution scheme of the relative displacement problem

How to solve the above-mentioned questions? A fast and practical method is interpolation. By interpolation, we can obtain new sample points within the range of a discrete set of known sample points. Actually, these new sample points construct a new phase screen which is different from that old phase screen which consists of those known sample points. In practical numerical simulations, it is not the old screen but the new screen that provides a propagating wave-front with phase perturbations because the spatial positions of new sample points is one-to-one correspondence with those of the grid points of the propagating wave-front. On the one hand, the phase value of the new screen and that of the old screen is highly correlated in that the nearest neighbors of a new sample point belong to the old screen [18]. On the other hand, the phase value of the new screen is different from that of the old screen due to the difference of their spatial positions. This is why Lane’s method relies on the subset of prior data represented by the m nearest neighbors when rendering a new point [18]. It is noteworthy that Lane’s method is mainly used to increase the resolution of a phase screen by consecutive subdivision [19,20] while the interpolation method in this paper is mainly used to obtain the phase value of any intermediate point in a line segment or square. For convenience, the value of a new sample point is obtained from two or four nearest points in the old screen in this paper. To be specific, the interpolation geometry is shown in Fig. 3 for the above-mentioned scenarios. The first kind of scenario is a one-dimensional interpolation problem. We have known the phase values of two endpoints: P₁ and P₂. Then, the phase value of the unknown sample point P can be obtained by interpolation. The second kind of scenario is a two-dimensional interpolation problem. We have known the phase values of the four corner points: P₁, P₂, P₃, and P₄. Then, the phase value of the unknown sample point P can be also obtained by interpolation.

Fig. 3 Graphical depiction of interpolation. The grid spacing is r. (a) One-dimensional interpolation of the first kind of scenario, r₁ and r₂ are the distances between P₁ and P, P₂ and P, respectively. (b) Two-dimensional interpolation of the second kind of scenario. uand v are the distances from a point P to the line P₁P₄, P₃P₄, respectively.

Download Full Size | PDF

In other words, a new phase screen, which represents the phase perturbations or the effects of turbulence on the optical field of LGS at a later time step, can be generated from an already existing phase screen by interpolation. For convenience, the existing phase screen will be called as the base phase screen while the new phase screen will be called as the interpolated phase screen in the following text. The ensemble statistics of the interpolated screen should match that of the base screen. The most commonly used ensemble statistics is the phase structure function. An established method for evaluating the accuracy of the interpolated phase screen is to compare the structure function derived from the interpolated screen with that from the base screen. The phase structure function is given by [21]:

D_{φ} (r) = E {{(φ (R + r) - φ (R))}^{2}} = 6.88 {(r / r_{0})}^{5 / 3}

which represents the average squared difference of the phase of the screen for pairs of points of varying location and separation r.

E {\cdot}

is the expectation operator.

r_{0}

is the atmospheric coherence length and is given by

r_{0} = {(0.423 k^{2} \int_{0}^{L} C_{n}^{2} (z) d z)}^{- 3 / 5}

. It characterizes the strength of the atmospheric turbulence. k is the wave number

2 π / λ

. λ is the wavelength.

C_{n}^{2}

is the refractive index structure constant dependent on location, time and altitude. L denotes the thickness of turbulent layer.

3. Description of interpolation methods

Here, we will describe two interpolation methods which can obtain a new data point from a discrete set of known data points. The first method is the well-known linear interpolation method and the second method is the statistical interpolation method based on the phase structure function, which is a statistical description of phase fluctuations caused by atmospheric turbulence. The statistical interpolation method is proposed in this paper for the first time within our knowledge.

3.1 Linear interpolation

Linear interpolation is a simple form of interpolation. If we have known the phase values of four points: P₁, P₂, P₃, and P₄. See Fig. 3 (b). Then, the phase value of the point P can be obtained by linear interpolation. It is commonly expressed in the form:

P = \frac{v}{r} (1 - \frac{u}{r}) P_{1} + \frac{u v}{r^{2}} P_{2} + \frac{u}{r} (1 - \frac{v}{r}) P_{3} + (1 - \frac{u}{r}) (1 - \frac{v}{r}) P_{4}

where P₁, P₂, P₃, and P₄ are the four corner points of a rectangle or a square and denote the values of these points at the same time (the same below). The interpolation geometry is similar to Fig. 3(a) when

u = 0

or

v = 0

in Eq. (2).

3.2 Statistical interpolation

The turbulence-induced fluctuations of optical phase are caused almost exclusively by small random fluctuations of atmospheric refractive-index [21]. The statistical nature of phase fluctuations can be well characterized by the phase structure function. Being different from the linear interpolator, the statistical interpolation is based on the linear interpolator between two or four existing points on the phase screen followed by a random displacement of the interpolated point. The coefficients and the magnitude of random displacement are determined by the phase structure function.

3.2.1 One-dimensional statistical interpolation

The geometry of one-dimensional interpolation is shown in Fig. 3 (a). Because the distance r, r₁ and r₂ are known, the following phase structure functions can be also known from Eq. (1):

D_{0} = E {{(P_{1} - P_{2})}^{2}} = 6.88 {(r / r_{0})}^{5 / 3}

D_{1} = E {{(P - P_{1})}^{2}} = 6.88 {(r_{1} / r_{0})}^{5 / 3}

D_{2} = E {{(P - P_{2})}^{2}} = 6.88 {(r_{2} / r_{0})}^{5 / 3}

We need to obtain the phase value of the sample point P from phase values of the known sample points P₁, P₂. We assume:

P = a P_{1} + b P_{2} + R

where a and b are weight coefficients and satisfy a + b = 1, R is a zero mean Gaussian variable with a variance

σ^{2}

, From Eq. (1), two equations can be derived:

{\begin{array}{c} E {{(P - P_{1})}^{2}} = {(1 - a)}^{2} D_{0} + σ^{2} = D_{1} \\ E {{(P - P_{2})}^{2}} = a^{2} D_{0} + σ^{2} = D_{2} \end{array}

Because D₀, D₁ and D₂ are known, we can solve Eq. (3) and obtain:

{\begin{array}{c} a = (D_{0} - D_{1} + D_{2}) / 2 D_{0} \\ σ^{2} = D_{2} - {(D_{0} - D_{1} + D_{2})}^{2} / 4 D_{0} \end{array}

Having known the coefficient a, another coefficient b can be obtained by the relation a + b = 1. Thus, by using Eqs. (3)-(5) the phase value of the sample point P can be obtained from the phase values of the known sample points P₁, P₂.

During the above derivation, we use the easily proven fact that for any three random variables, sayα,β, andχ, having zero mean and equal variances [16]:

〈 (α - χ) (β - χ) 〉 = \frac{1}{2} [〈 {(α - χ)}^{2} 〉 + 〈 {(β - χ)}^{2} 〉 - 〈 {(α - β)}^{2} 〉]

The terms on the right-hand side of Eq. (6) correspond to the structure function of relating points.

Two special cases can be considered. The first case is the point P coincides with the endpoint P₁ or P₂. The coefficient (a, b) is equal to (1, 0) or (0, 1) and the variance $σ^{2}$ is equal to 0. The second case is the point P is the midpoint of P₁P₂. The coefficient (a, b) is equal to (0.5, 0.5) and D₁ is equal to D₂. The variance $σ^{2}$ is equal to D₁-0.25D₀. Obviously, the results of the second case are consistent with those of the known random mid-point displacement method [19,20].

3.2.2 Two-dimensional statistical interpolation

The geometry of two-dimensional interpolation is shown in Fig. 3 (b). We need obtain the value of the sample point P from the values of the known sample points: P₁, P₂, P₃, and P₄. Similarly to section 3.2.1, we also assume:

P = a P_{1} + b P_{2} + c P_{3} + d P_{4} + R

where a, b, c and d are weight coefficients and a + b + c + d = 1, R is zero mean Gaussian variables with a variance

σ^{2}

. From Eq. (1), there are four equations:

{\begin{array}{c} E {{(P - P_{1})}^{2}} = D_{1} \\ E {{(P - P_{2})}^{2}} = D_{2} \\ E {{(P - P_{3})}^{2}} = D_{3} \\ E {{(P - P_{4})}^{2}} = D_{4} \end{array}

And assume:

E {{(P_{1} - P_{2})}^{2}} = E {{(P_{2} - P_{3})}^{2}} = E {{(P_{3} - P_{4})}^{2}} = E {{(P_{4} - P_{1})}^{2}} = D_{0}

E {{(P_{1} - P_{3})}^{2}} = E {{(P_{2} - P_{4})}^{2}} = D_{5}

D₀, D₁, D₂, D₃, D₄ and D₅ are all known.

We can obtain the following equations by using the Eqs. (6), (7) and (8):

{\begin{array}{c} {(a - 1)}^{2} D_{0} + b^{2} D_{5} + c^{2} D_{0} + b (a - 1) D_{5} + c (a - 1) (2 D_{0} - D_{5}) + b c D_{5} + σ^{2} = D_{1} (9.1) \\ a^{2} D_{0} + {(b - 1)}^{2} D_{5} + c^{2} D_{0} + a (b - 1) D_{5} + a c (2 D_{0} - D_{5}) + c (b - 1) D_{5} + σ^{2} = D_{2} (9.2) \\ a^{2} D_{0} + b^{2} D_{5} + {(c - 1)}^{2} D_{0} + a b D_{5} + a (c - 1) (2 D_{0} - D_{5}) + b (c - 1) D_{5} + σ^{2} = D_{3} (9.3) \\ a^{2} D_{0} + b^{2} D_{5} + c^{2} D_{0} + a b D_{5} + a c (2 D_{0} - D_{5}) + b c D_{5} + σ^{2} = D_{4} (9.4) \end{array}

Solving Eq. (9).1)~(9.4) and obtaining:

{\begin{array}{c} a = \frac{2 D_{0} D_{5} + 4 D_{0} D_{3} - 4 D_{0} D_{1} - D_{5}^{2} + D_{4} D_{5} - 3 D_{3} D_{5} + D_{2} D_{5} + D_{1} D_{5}}{4 D_{5} (2 D_{0} - D_{5})} \\ b = \frac{2 D_{0} D_{5} + 4 D_{0} D_{4} - 4 D_{0} D_{2} - D_{5}^{2} + D_{3} D_{5} - 3 D_{4} D_{5} + D_{2} D_{5} + D_{1} D_{5}}{4 D_{5} (2 D_{0} - D_{5})} \\ c = \frac{2 D_{0} D_{5} + 4 D_{0} D_{1} - 4 D_{0} D_{3} - D_{5}^{2} + D_{3} D_{5} - 3 D_{1} D_{5} + D_{2} D_{5} + D_{4} D_{5}}{4 D_{5} (2 D_{0} - D_{5})} \\ d = \frac{2 D_{0} D_{5} + 4 D_{0} D_{2} - 4 D_{0} D_{4} - D_{5}^{2} + D_{4} D_{5} - 3 D_{2} D_{5} + D_{3} D_{5} + D_{1} D_{5}}{4 D_{5} (2 D_{0} - D_{5})} \end{array}

After having known the coefficients a, b, c and d, the variance

σ^{2}

can be obtained from Eq. (9).4).

Similarly to the one-dimensional statistical interpolation, two special cases can be also considered. The first case is the point P coincides with the corner point P₁, P₂, P₃, or P₄. The coefficient (a, b, c, d) is equal to (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0) or (0, 0, 0, 1) and the variance $σ^{2}$ is equal to 0. The second case is the point P is the center point of the square P₁P₂P₃P₄. The coefficient (a, b) is equal to (0.25, 0.25, 0.25, 0.25) and D₁ = D₂ = D₃ = D₄. The variance $σ^{2}$ is equal to $σ^{2} = D_{1} - \frac{1}{4} D_{0} - \frac{1}{8} D_{5}$ . Obviously, the results of the second case are consistent with those of the known random mid-point displacement method [19,20].

In comparison with the linear interpolator, the statistical interpolator introduces an additional term according to the phase structure function. Compared to the random mid-point displacement method, the statistical interpolator is a generalization of that method. Initially, the mid-point displacement method is used to generate a Kolmogorov phase screen by a process of successive subdivision in order to increase the resolution of the phase screen [19]. In order to achieve a more accurate statistics, an improved method (for convenience, we call it Harding’ algorithm hereafter) was proposed that firstly generates an exact low resolution phase screen by factorizing its covariance matrix and then increases the resolution of the phase screen by the random mid-point displacement method [20]. But the interpolator must be done in several consecutive steps in order to produce a higher resolution and only be used to obtain the phase value of the mid-point in a line segment or the center point in a square. In this paper, our statistical interpolator can obtain the value of any point in a line or a square and not limited to the midpoint in a line or the center point in a square. It can be used to overcome the disadvantage of the so-called moving phase screens method [11]. That is to say, it can make the screen be shifted by non-integer multiples of the grid spacing in the simulations.

4. Results and discussions

In order to provide a complete evaluation of the statistical interpolation method, we will implement some comparison analysis. At first, we will compare the phase structure function derived from the two interpolation methods with that of the base screen. Then, some numerical experiments are also designed to test the accuracy of the two interpolation methods. These numerical experiments are mainly for the light propagation through atmospheric turbulence, including open-loop (without an AO system) and closed-loop (with an AO system).

4.1 Comparison of the phase structure function

4.1.1 One-dimensional interpolation

As mentioned above, the phase perturbations caused by atmospheric turbulence are actually a new phase screen which is obtained from an existing phase screen by interpolation. Its ensemble statistics, namely the phase structure function, should match very closely the phase structure function of the base screen. Following that, we can evaluate the accuracy of the two interpolation methods by comparing the average phase structure function derived from an ensemble of the interpolated phase screens with that derived from an ensemble of the base phase screen. First, we generate a base phase screen. Then, an interpolated phase screen is generated by means of interpolating between the points of the base screen. In order to obtain the ensemble statistics, two sets of samples of 10000 base phase screens and interpolated phase screens are generated to compute the phase structure function. We define a relative error of the structure functions as an evaluation parameter. It is given by:

E r r o r (r) = 1 - D_{ϕ}^{s i m u l a t i o n} (r) / D_{ϕ}^{base} (r)

where

D_{ϕ}^{base} (r)

and

D_{ϕ}^{s i m u l a t i o n} (r)

denotes the average phase structure function computed from the base phase screens and the average phase structure function computed from the interpolated phase screens, respectively.

Covariance method [20] can obtain more accurate statistical property of a distorted wave-front and it is used to generate the base phase screens here. We can use Eqs. (2) and (3) to generate a new phase screen by interpolating a base phase screen. For the one-dimensional interpolation method, the grid spacing of the interpolated screen should be equal to that of the base screen. Some parameters which are used to generate the base phase screens are as follows. The atmospheric coherence length r₀ = 0.20 m. The grid spacing is 0.01 m. The effects of the position of the interpolated point P on the accuracy of two interpolation methods are also investigated. Assuming $r_{1} = h * r$ ( $0 \leq h \leq 1$ ) in Fig. 3(a), three different values for h are chosen: 0.3, 0.5, 0.7.The grid points of the base phase screen are 33 × 34 and the grid points of the interpolated phase screen are 33 × 33. A comparison of the relative error is shown in Fig. 4 .

Fig. 4 The relative error of the phase structure functions for different interpolation methods. (a) h = 0.3, (b) h = 0.5, (c) h = 0.7.

Download Full Size | PDF

It is shown in Fig. 4 that there exists substantial difference in the high frequency region. The relative error brought by the linear interpolator is bigger than that brought by the statistical interpolator in the high frequency region. Compared with the structure function of the base screen, the maximum error brought by statistical interpolation is about 3%, while the maximum error brought by linear interpolation is about 16% when the parameter h is equal to 0.5. In the low frequency region, there is no substantial difference between the linear interpolator and the statistical interpolator. In addition, a maximum error occurs when the parameter h is equal to 0.5 in the high frequency, especially for the linear interpolator.

4.1.2 Two-dimensional interpolation

For the two-dimensional interpolation method, we can use Eqs. (2) and (7) to generate a new phase screen by interpolating a base phase screen. That is to say, each grid point of the new screen is obtained by interpolating from the grid points of the base screen. In comparison with the one-dimensional interpolation method, the two-dimensional interpolation method is more complicated. On the one hand, the positions of the interpolated points may affect the accuracy of the two interpolation methods. It is shown in Fig. 3 (b). The position of the interpolated point P is determined by the parameters u and v. Similarly to the case of one-dimensional interpolation, we can assume $u = m * r$ and $v = n * r$ ( $0 \leq m, n \leq 1$ ). So the position of the interpolated point P can be changed in the square grid P₁P₂P₃P₄ by adjusting the parameter mand n. On the other hand, the ratiolof the grid spacing of the base screen to that of the interpolated screen also affect the accuracy of the interpolation methods. Being different from the one-dimensional interpolator, the grid spacing of the interpolated screen is generally less than or equal to that of the base screen in the case of two-dimensional interpolation. That is $l \geq 1$ .This means that the interpolated screen has more grid points than the base screen in the same area. See Fig. 2 (b). In order to obtain an interpolated screen from a base screen, one or more grid points are obtained by interpolating between the four corner points of a square in the base screen. The grid spacing of the interpolated screen is equal to that of the base screen when only one point can be interpolated from the four corner points of a square in the base screen. Otherwise, the grid spacing of the interpolated screen is less than that of the base screen when more points are interpolated from the four corner points of a square in the base screen. Thus, we will first investigate how the position of the interpolated point and the ratio l of the grid spacing affect the accuracy of the interpolation methods.

Some parameters which are used to generate the base screen are as follows. The atmospheric coherence length r₀ = 0.14 m. The grid spacing of the base screen is 0.01 m. The grid points of the base screen are 34 × 34 and those of the interpolated screen are 33 × 33.

In order to investigate the effect of the ratio l on the accuracy of the interpolation method, we change the value of l under the fixed value m and n. The ratio l is chosen as 1.0, 1.2, 1.6 and 2.3 while $(m, n) = (0.3, 0.8)$ .The comparison of the relative error of the phase structure function is shown in Fig. 5 . The results show that the relative error increases with the ratio l in the high frequency region while remains comparatively small in the low frequency region. The ratio l must be equal to 1.0 in order to minimize the relative error in the high frequency region. But, it is generally greater than 1.0 in the simulation of actual scenarios, i.e. the grid spacing of the base screen is greater than that of the interpolated screen. Fortunately, the random mid-point displacement method can be used to solve this problem. By using this method, the local region of the base screen, which is corresponding to the interpolated screen, can be subdivided into a phase screen with higher resolution until the grid spacing of this screen is equal to that of the interpolated screen, i.e. $l = 1.0$ . Then the interpolated screen can be obtained by interpolating the screen with higher resolution. Thus, we only discuss the case of $l = 1.0$ in the following text.

Fig. 5 Relative error of the phase structure function for different interpolation methods at (a) l = 1.0, (b) l = 1.2, (c) l = 1.6, (d) l = 2.3.

Download Full Size | PDF

Similarly to the one-dimensional interpolator, the relative error brought by the linear interpolator is also bigger than that brought by the statistical interpolator in the high frequency region.

In order to investigate the effect of the position of the interpolated point on the accuracy of the interpolation method, we change the value of m and n under a fixed value of $l = 1.0$ . The comparison of the relative error of the phase structure function is shown in Fig. 6 . Similarly to that of the one-dimensional interpolator, the results show that the relative error arrives at a maximum value when $(m, n) = (0.5, 0.5)$ . The maximum error brought by the statistical interpolator is about 7% while the maximum error brought by the linear interpolator is about 25% when the interpolated point lies in the center of the square grid. The closer to the endpoint the interpolated point is, the smaller the relative error is. The relative error in the high frequency region is bigger than that in the low frequency region. In addition, the relative error brought by the linear interpolator is also bigger than that brought by the statistical interpolator in the high frequency region.

Fig. 6 The relative error of the phase structure function for different interpolation methods at (a) (m, n) = (0.1, 0.1), (b) (m, n) = (0.5, 0.5), (c) (m, n) = (0.8, 0.2), (d) (m, n) = (0.9, 0.8).

Download Full Size | PDF

4.2 Comparison of simulation results of light propagation through turbulent atmosphere

In addition to the comparison of the phase structure function, numerical experiments are carried out for a plane wave propagating in the atmospheric turbulence with/without AO for some scenarios. By reasonably designing a numerical experimental scheme, the effect of different interpolation methods on the light propagation through the turbulent atmosphere with/without AO is investigated. The results in Section 4.1 have shown that the relative error reaches a maximum in the high frequency region when the interpolated point is located at the mid-point in a line segment or at the center point in a square grid. Thus, we only discuss this specific case.

This numerical experimental scheme is designed as follows. Firstly, a high resolution phase screen is created by Harding’s algorithm [20]. Secondly, we extract a point every a point from the high resolution screen. These extracted points constitute a low resolution phase screen. Thus, the grid spacing of the low resolution screen is two times that of the high resolution screen. It is the low resolution screen that is used in the simulation computation of the light propagation through atmospheric turbulence with/without AO. The long-exposure Strehl ratio is used as the evaluation parameter in the numerical experiments.

The scheme is illustrated in Fig. 7 . For convenience, only 5 × 5 grid points are displayed in Fig. 7 (a). The dark points denote the grid points of the high resolution screen while the hollow points denote the grid points of the low resolution screen. It is shown in Fig. 7 (b). The grid point P is a grid point of the low resolution screen. The corresponding point in the high resolution screen (i.e. the base screen) is the grid point P₁₀. The phase value of each grid point in the low resolution screen can be obtained by three approaches. The first approach is that the value of the grid point P₁₀ in the high resolution screen is directly taken as the phase value of the grid point P in the low resolution screen. The second approach is that the phase value of the grid point P is obtained by interpolating between two nearest points P₂ and P₆ or P₄ and P₉. The grid point P is the interpolated point and it is also the mid-point of the line segment P₂P₆ and P₄P₉ at the same time. This is a case of one-dimensional interpolation. The third approach is that the phase value of the grid point P is obtained by interpolating between four nearest corner points P₁, P₃, P₅ and P₇. The grid point P is the interpolated point and it is also the center point of the square grid P₁P₃P₅P₇. This is a case of two-dimensional interpolation. Two interpolators (i.e. the linear interpolator and the statistical interpolator) are used to realize the interpolation process. Obviously, the first approach can provide an accurate result because this approach doesn’t use any interpolation. Thus, this accurate result provided by the first approach will be used as a standard value for comparing to evaluate two interpolation methods. The result provided by the linear interpolator and that provided by the statistical interpolator can be compared to the standard value. A relative error of the long-exposure Strehl ratio is computed. In this way, we can quantitatively estimate the accuracy and reliability of two interpolation methods.

Fig. 7 Generation of the low-resolution screen (hollow points) by extracting a point every pixel from the high-resolution screen (dark points).

Download Full Size | PDF

The simulation algorithm of laser propagation in the atmosphere with/without AO is described by Yan et al [2] in details. Here we will use this simulation algorithm to investigate the effects of interpolation methods on the simulation results of the light propagation. The simulated beam is chosen as a plane wave beam. It propagates along a 3000 m horizontal path, which was divided into 15 equivalent phase screens, each covering a distance z = 100 m. The wavelength is $λ = 1.55 μ m$ . The computational grid is a 256 × 256 mesh. The Rytov variance $σ_{R}^{2}$ can be used as a measure of the strength of the atmospheric turbulence [21]. For plane waves, it is defined as $σ_{R}^{2} = 1.23 C_{n}^{2} k^{7 / 6} L^{11 / 6}$ . Here, different Rytov variances can be obtained by changing the turbulence strength $C_{n}^{2}$ . The Strehl ratio is defined as the ratio of the intensity of the center point on the target after propagation through the turbulent atmosphere to that through a vacuum.

Considering the inevitable time delay of an adaptive optics system [22], random phase screen must move in a certain direction at a definite speed and a rectangular phase screen must be used in the simulation. The effects of the time delay are treated as follows. The phase screens are assumed to move laterally in the delay time period because of lateral wind and/or the lateral movements of the laser and/or target. Thus, the laser beam with compensated phase will propagate through the laterally moved phase screens to reach the target. Firstly, an initial phase screen with 1025 × 1025 grid points is generated by Harding’ algorithm. Secondly, a high resolution screen with 1025 × 513 grid points is formed by choosing one half of the initial screen. Finally, a low resolution phase screen with 512 × 256 can be created by extracting a point every pixel from the high resolution phase screen. The phase value of each point in the low resolution screen is obtained by three above-mentioned approaches.

Firstly, we investigate the effects of the interpolation methods on the light propagation results under the condition of one-dimensional interpolation. The one-dimensional interpolation method includes the linear interpolator and the one-dimensional statistical interpolator. In Fig. 8 , the relative error of the long-exposure Strehl ratio is shown. It can be seen that the result provided by the linear interpolator is greater than that provided by the one-dimensional statistical interpolator under the condition of open-loop (without AO) and closed-loop (with AO). For two interpolators, the closed-loop result is greater than the open-loop result. Under the condition of open-loop and closed-loop, the relative error both increases with the Rytov variance.

Fig. 8 Relative error of the long-exposure Strehl ratio from the first approach (the standard value) and that from the second approach which uses the linear interpolation given by Eq. (2) or the one-dimensional statistical interpolation given by Eq. (3). (a) Open-loop; (b) Closed-loop.

Download Full Size | PDF

Secondly, we also investigate the effects of the interpolation methods on the light propagation results under the condition of two-dimensional interpolation. The two-dimensional interpolation method includes the linear interpolator and the two-dimensional statistical interpolator. In Fig. 9 , the relative error of the long-exposure Strehl ratio is shown. Similarly to the case of one-dimensional interpolation, it is also shown that the result provided by the linear interpolator is greater than that provided by the two-dimensional statistical interpolator under the condition of open-loop and closed-loop. For the linear interpolator, the closed-loop result is greater than the open-loop result. For the two-dimensional statistical interpolator, the relative error seems to smoothly fluctuate with the Rytov variance increasing and both remain at a rather small value under the condition of open-loop and closed-loop.

Fig. 9 Relative error of the long-exposure Strehl ratio from the first approach (the standard value) and that from the second approach which uses the linear interpolator given by Eq. (2) or the two-dimensional statistical interpolator given by Eq. (6). (a) Open-loop; (b) Closed-loop.

Download Full Size | PDF

From the above analysis, it is shown that the closed-loop result of the relative error is greater than the open-loop result for the linear interpolator in two interpolation cases, i.e., the one-dimensional interpolation and the two-dimensional interpolation. The reason for this difference lies mainly in that only one interpolation is performed under the condition of open-loop while two interpolations are performed under the condition of closed-loop. The latter makes the result of the relative error larger than the former does. For the statistical interpolator, the result of one-dimensional interpolation is larger than that of two-dimensional interpolation under the condition of open-loop and closed-loop. This is because the latter (4-nearest-neighbor statistical interpolator) can provide a more accurate statistical property than the former (2-nearest-neighbor statistical interpolator) by introducing a correlation between more adjacent points in the high resolution screen. Due to this fact, the 4-nearest-neighbor statistical interpolator produces a more accuracy under the condition of open-loop and closed-loop (See Fig. (9)). Although in theory a statistical interpolator utilizing more than 4 nearest neighbors would provide better results, it increases the complexity of the interpolator greatly.

Concerning the comparison between the results of the linear interpolator and those of the statistical interpolator employed in the light propagation through the turbulent atmosphere with/without AO, we observed that the statistical interpolator can provide a more accurate result than the linear interpolator. It is well known that the linear interpolator is a kind of weighted average interpolator and smoothes phase fluctuations. Thus, it artificially alleviates the effect of atmospheric turbulence. That is, the interpolated screen is smoother than the base screen, especially in the high frequency region. This will impose an important effect on the numerical simulation of light propagation through the turbulent atmosphere with/without AO. Considering the random fluctuation property of atmospheric turbulence, the statistical interpolator is based on a linear interpolation between two/four nearest existing points on the phase screen plus a random displacement of the interpolated point. The magnitude of the random displacement is determined by the turbulence strength. From the comparison of phase structure function and light propagation through atmospheric turbulence, it concludes that the statistical interpolator is more accurate and reliable than the linear interpolator.

5. Summary

Turbulent phase screen has been widely used to represent the phase perturbations induced by a turbulent media in numerical simulation of light propagation through the atmosphere. But these perturbations cannot be directly obtained from the phase screen in lots of actual scenarios, such as scenarios arising from temporal evolution conditions or astronomical telescope with a laser guide star (LGS) adaptive optics system, due to the movement of phase screen(s) or the mismatch between the optical field of astronomical object and that of the LGS. A statistical interpolation method is proposed to address this problem in this paper. This method is actually an extension of the known random midpoint displacement algorithm. It can be used to obtain the value of any point in a line or a square but nor limited to the midpoint. Taking the results provided by the base phase screen(s) as the standard for comparing, the phase structure function and the numerical results of light propagation through the turbulent atmosphere with/without AO are used to evaluate the statistical interpolator. The statistical interpolator is also compared with the well-known linear interpolator. For the comparison of the relative error of the phase structure function, it is shown that there exists substantial difference in the high frequency region. In addition, a maximum error occurs when the interpolated point is the mid-point of the line segment or the center point of the square grid in the high frequency, especially for the linear interpolator. For the comparison of light propagation through the turbulent atmosphere, it is also shown that the statistical interpolator can obtain better results than the linear interpolator. The difference between the open-loop results of the statistical interpolator and those of the linear interpolator is larger than that between the close-loop results of the statistical interpolator and those of the linear interpolator. It can conclude that the statistical interpolator is more accurate and reliable than the linear interpolator.

References and links

1. J. A. Fleck Jr, J. R. Morris, and M. D. Feit, “Time-dependent propagation of high energy laser beams through the atmosphere,” Appl. Phys. (Berl.) 10(2), 129–160 (1976). [CrossRef]

2. H.-X. Yan, S.-S. Li, D.-L. Zhang, and S. Chen, “Numerical simulation of an adaptive optics system with laser propagation in the atmosphere,” Appl. Opt. 39(18), 3023–3031 (2000). [CrossRef]

3. L. C. Andrews, R. L. Philips, R. J. Sasiela, and R. R. Parenti, “Strehl ratio and scintillation theory for uplink Gaussian-beam waves: beam wander effects,” Opt. Eng. 45(7), 076001–1 (2006). [CrossRef]

4. F. Dios, J. Recolons, A. Rodríguez, and O. Batet, “Temporal analysis of laser beam propagation in the atmosphere using computer-generated long phase screens,” Opt. Express 16(3), 2206–2220 (2008). [CrossRef] [PubMed]

5. H.-X Yan, Han-Ling Wu, Shu-Shan Li and She Chen, “Cone effect in astronomical adaptive optics system investigated by a pure numerical simulation,” Proc. SPIE 5903, 5903OU1–12 (2005).

6. E. M. Johansson and D. T. Gavel, “Simulation of stellar speckle imaging,” Proc. SPIE 2220, 372–383 (1994). [CrossRef]

7. V. Sriram and D. Kearney, “An ultra fast Kolmogorov phase screen generator suitable for parallel implementation,” Opt. Express 15(21), 13709–13714 (2007). [CrossRef] [PubMed]

8. A. Beghi, A. Cenedese, and A. Masiero, “Stochastic realization approach to the efficient simulation of phase screens,” J. Opt. Soc. Am. A 25(2), 515–525 (2008). [CrossRef]

9. H. Jakobsson, “Simulations of time series of atmospherically distorted wave fronts,” Appl. Opt. 35(9), 1561–1565 (1996). [CrossRef] [PubMed]

10. H.-X. Yan, S.-S. Li, and S. Chen, “Numerical simulation investigations of the dynamic control process and frequency response characteristics in an adaptive optics system,” Proc. SPIE 4494, 156–166 (2002). [CrossRef]

11. B. M. Welsh, “Fourier-series-based atmospheric phase screen generator for simulating anisoplanatic geometries and temporal evolution,” Proc. SPIE 3125, 327–338 (1997). [CrossRef]

12. G. Sedmak, “Performance analysis of and compensation for aspect-ratio effects of fast-fourier-transform-based simulations of large atmospheric wave fronts,” Appl. Opt. 37(21), 4605–4613 (1998). [CrossRef]

13. M. C. Roggemann, B. M. Welsh, D. Montera, and T. A. Rhoadamer, “Method for simulating atmospheric turbulence phase effects for multiple time slices and anisoplanatic conditions,” Appl. Opt. 34(20), 4037–4051 (1995). [CrossRef] [PubMed]

14. H.-X. Yan, S. Chen, and S.-S. Li, “Turbulent phase screens generated by covariance approach and their application in numerical simulation of atmospheric propagation of laser beam,” Proc. SPIE 6346, 634628 (2006). [CrossRef]

15. F. Assémat, R. W. Wilson, and E. Gendron, “Method for simulating infinitely long and non stationary phase screens with optimized memory storage,” Opt. Express 14(3), 988–999 (2006). [CrossRef] [PubMed]

16. D. L. Fried and T. Clark, “Extruding Kolmogorov-type phase screen ribbons,” J. Opt. Soc. Am. A 25(2), 463–468 (2008). [CrossRef]

17. C. S. Gardner, B. M. Welsh, and L. A. Thompson, “Design and performance analysis of adaptive optical telescopes using laser guide stars,” Proc. IEEE 78(11), 1721–1743 (1990). [CrossRef]

18. B. Formwalt and S. Cain, “Optimized phase screen modeling for optical turbulence,” Appl. Opt. 45(22), 5657–5668 (2006). [CrossRef] [PubMed]

19. R. G. Lane, A. Glindemann, and J. C. Dainty, “Simulation of Kolmogorov phase screen,” Waves Random Media 2(3), 209–224 (1992). [CrossRef]

20. C. M. Harding, R. A. Johnston, and R. G. Lane, “Fast simulation of a kolmogorov phase screen,” Appl. Opt. 38(11), 2161–2170 (1999). [CrossRef]

21. L. C. Andrews, and R. L. Phillips, Laser Beam Propagation through Random Media, 2nd ed. (SPIE Press, Washington, 2005).

22. R. K. Tyson, Principles of Adaptive Optics, 2nd ed. (Academic Press, Boston, 1997).

Statistical Interpolation Method of Turbulent Phase Screen

Abstract

1. Introduction

2. Problem statement

2.1 One-dimensional relative displacement

2.2 Two-dimensional relative displacement

2.3 Solution scheme of the relative displacement problem

3. Description of interpolation methods

3.1 Linear interpolation

3.2 Statistical interpolation

3.2.1 One-dimensional statistical interpolation

3.2.2 Two-dimensional statistical interpolation

4. Results and discussions

4.1 Comparison of the phase structure function

4.1.1 One-dimensional interpolation

4.1.2 Two-dimensional interpolation

4.2 Comparison of simulation results of light propagation through turbulent atmosphere

5. Summary

References and links

Cited By

Figures (9)

Equations (16)

Optics Express