Simulation of anisoplanatic imaging through optical turbulence using numerical wave propagation with new validation analysis

Russell C. Hardie; Jonathan D. Power; Daniel A. LeMaster; Douglas R. Droege; Szymon Gladysz; Santasri Bose-Pillai

doi:10.1117/1.OE.56.7.071502

6 February 2017 Simulation of anisoplanatic imaging through optical turbulence using numerical wave propagation with new validation analysis

Russell C. Hardie, Jonathan D. Power, Daniel A. LeMaster, Douglas R. Droege, Szymon Gladysz, Santasri Bose-Pillai

Author Affiliations +

Optical Engineering, Vol. 56, Issue 7, 071502 (February 2017). https://doi.org/10.1117/1.OE.56.7.071502

Abstract

We present a numerical wave propagation method for simulating imaging of an extended scene under anisoplanatic conditions. While isoplanatic simulation is relatively common, few tools are specifically designed for simulating the imaging of extended scenes under anisoplanatic conditions. We provide a complete description of the proposed simulation tool, including the wave propagation method used. Our approach computes an array of point spread functions (PSFs) for a two-dimensional grid on the object plane. The PSFs are then used in a spatially varying weighted sum operation, with an ideal image, to produce a simulated image with realistic optical turbulence degradation. The degradation includes spatially varying warping and blurring. To produce the PSF array, we generate a series of extended phase screens. Simulated point sources are numerically propagated from an array of positions on the object plane, through the phase screens, and ultimately to the focal plane of the simulated camera. Note that the optical path for each PSF will be different, and thus, pass through a different portion of the extended phase screens. These different paths give rise to a spatially varying PSF to produce anisoplanatic effects. We use a method for defining the individual phase screen statistics that we have not seen used in previous anisoplanatic simulations. We also present a validation analysis. In particular, we compare simulated outputs with the theoretical anisoplanatic tilt correlation and a derived differential tilt variance statistic. This is in addition to comparing the long- and short-exposure PSFs and isoplanatic angle. We believe this analysis represents the most thorough validation of an anisoplanatic simulation to date. The current work is also unique that we simulate and validate both constant and varying Cn2(z) profiles. Furthermore, we simulate sequences with both temporally independent and temporally correlated turbulence effects. Temporal correlation is introduced by generating even larger extended phase screens and translating this block of screens in front of the propagation area. Our validation analysis shows an excellent match between the simulation statistics and the theoretical predictions. Thus, we think this tool can be used effectively to study optical anisoplanatic turbulence and to aid in the development of image restoration methods.

1. Introduction

Long range imaging in the atmosphere is impacted, and often limited, by optical turbulence. Random variations in the index of refraction along the optical path are caused by temperature variations and convection.¹ This leads to warping and blurring degradations. When the warping and blurring are spatially varying over the field of view of an imaging sensor, this is referred to as anisoplanatic imaging. Anisoplanatic conditions typically prevail when imaging extended scenes. With very small fields of view, the optical turbulence may be reasonably modeled as isoplanatic, using a spatially invariant point spread function (PSF). This is often the case with astronomical imaging.¹

The simulation of imaging under isoplanatic atmospheric turbulence has been well studied.¹ Commercial and free open-source software²^,³ is available for this purpose. However, anisoplanatic imaging simulations using numerical wave propagation have only recently been described in the literature. Some of the first such papers include those of Carrano⁴ and Praus et al.⁵ Further wave propagation-based anisoplanatic simulation works include that of Bos and Roggemann,⁶^,⁷ and Monnier et al.⁸^,⁹ Anisoplanatic simulations that do not involve wave propagation have also been developed.¹⁰^–¹² The ability to accurately simulate the degradation effects of optical turbulence is important for several reasons. First, such simulations allow us to study the impact of atmospheric turbulence under a wide variety of imaging scenarios. Second, we are able to use simulated data to quantitatively evaluate turbulence mitigation methods.⁴^,¹³^–¹⁵ Quantitative analysis is possible because we have an objective truth image from which the degraded images are generated. Without such a truth image, assessment of restoration algorithms is subjective, and optimization cannot be automated. We believe that quantitative performance analysis is critical to the advancement of image restoration methods.

The simulation methodology presented here is based on that of Bos and Roggemann.⁷ We expand on, and extend, that work in several ways. We also provide a complete description of the simulation, including the details of the wave propagation method used. Like Bos and Roggemann,⁷ our approach computes an array of PSFs for a two-dimensional (2-D) grid on the object plane. The PSFs are then used in a spatially varying weighted sum operation with an ideal image to produce a simulated image with realistic optical turbulence degradation. The degradation includes spatially varying warping and spatially varying blurring.

To produce the PSF array, we generate a series of extended phase screen realizations. Simulated point sources are numerically propagated from different positions on the object plane, through the phase screens, and ultimately to the focal plane of the simulated camera. Note that the optical path for each PSF will be different, and thus, pass through a different portion of the extended phase screens. These different paths give rise to spatially varying, but spatially correlated PSFs. As we shall show, these PSFs may be used to generate accurate anisoplanatic effects. The method we use to define the individual phase screen statistics is distinct from that of Bos and Roggemann.⁷ Our approach is based on the constrained least squares optimization presented by Schmidt,² but is extended to include isoplanatic angle in the cost function. Another unique feature of our method is that we exclude the phase screen at the pupil plane. This aides in generating the appropriate level of anisoplanatic PSF correlations.

We also present a validation analysis here. In particular, we compare the simulated outputs with the theoretical anisoplanatic tilt correlation,¹⁶ and a derived differential tilt variance (DTV) statistic. This is in addition to comparing the long- and average short-exposure PSFs and isoplanatic angle. We believe this analysis represents the most thorough validation of an anisoplanatic simulation to date. The current work is also unique that we simulate and validate both constant and varying $C_{n}^{2} (z)$ profiles. Furthermore, we simulate sequences with both temporally independent and temporally correlated turbulence effects. Temporal correlation is introduced by generating larger extended phase screens and translating this block of screens in front of the propagation area. This approach is similar to that described by Dios et al.¹⁷ Our validation analysis shows an excellent match between the simulation statistics and the theoretical predictions. Thus, we believe this approach can be used effectively to study optical anisoplanatic turbulence and to aid in the development of image restoration methods.

The rest of this paper is organized as follows. Key anisoplanatic turbulence statistics are presented and discussed in Sec. 2. These statistics are used in setting up the simulation and in the validation analysis. The details of the simulation are presented in Sec. 3. The experimental results are presented in Sec. 4. Finally, we offer conclusions in Sec. 5.

2. Anisoplanatic Optical Turbulence

2.1.

Optical Turbulence Statistics

Variations in the index of refraction in the atmosphere can be modeled with a refractive index structure function. Using a Kolmogorov model,¹ this is given by

Eq. (1)

D_{n} (r) = ⟨ {[n (x + r) - n (x)]}^{2} ⟩ = C_{n}^{2} r^{2 / 3},

where

x

and

r

are three-dimensional spatial coordinate vectors

r = | r |

,

n (\cdot)

is the index of refraction, and

⟨ \cdot ⟩

represents an ensemble mean operator. The parameter

C_{n}^{2}

is the refractive index structure parameter.¹ The units of

C_{n}^{2}

are

m^{- 2 / 3}

and typical values tend to range from

1 \times 10^{- 13}

to

1 \times 10^{- 17}

.² When this quantity varies along the optical path, it may be expressed as a function

C_{n}^{2} (z)

, where

z

is the distance from the source.

The atmospheric coherence diameter (or Fried parameter)² can be expressed as a weighted integral of $C_{n}^{2} (z)$ yielding

Eq. (2)

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

where

k = 2 π / λ

is the wavenumber and

λ

is the wavelength. Note that this is for spherical wave propagation where

z = 0

is at the source, and

z = L

is at the camera. The Fried parameter is directly linked to the atmospheric PSF and the tilt variance for a point source. A small

r_{0}

, relative to the camera aperture, indicates a high level of atmospheric blurring and a high tilt variance.

Another important statistic is the isoplanatic angle.¹ Two point sources that are separated by less than the isoplanatic angle will have a mean wave function phase difference at the aperture of $< 1$ radian.²^,¹⁸ Another way to think of this is that points separated by less than this angle will have approximately the same PSF. The isoplanatic angle can also be expressed as a weighted integral of $C_{n}^{2} (z)$ , yielding²^,¹⁸

Eq. (3)

θ_{0} = {[2.91 k^{2} L^{5 / 3} \int_{z = 0}^{z = L} C_{n}^{2} (z) {(1 - \frac{z}{L})}^{5 / 3} d z]}^{- 3 / 5} .

A statistic that relates to fluctuations in the wave function amplitude is the log-amplitude variance,² given by

Eq. (4)

σ_{χ}^{2} = 0.563 k^{7 / 6} L^{5 / 6} \int_{z = 0}^{z = L} C_{n}^{2} (z) {(\frac{z}{L})}^{5 / 6} {(1 - \frac{z}{L})}^{5 / 6} d z .

It is interesting to visualize the three

C_{n}^{2} (z)

weighting functions corresponding to

r_{0}

,

θ_{0}

, and

σ_{χ}^{2}

. This is done in Fig. 1. Note that the isoplanatic angle is most impacted by turbulence at the source, while

r_{0}

is most impacted by turbulence near the camera. The log-amplitude variance is impacted most by the center of the optical path. Since the weighting for these three statistics covers the optical path in this balanced manner, we use them to determine the phase screen Fried parameters in our simulation, as shown in Sec. 3.5.

Fig. 1

Weighting functions for the $C_{n}^{2} (z)$ corresponding to three of the key optical turbulence parameters for $L = 7 km$ . Isoplanatic angle is heavily weighted toward the source end of path. The Fried parameter is weighted toward the camera end of the path, and the log amplitude variance is weighted primarily in the middle.

2.2.

Optical Transfer Functions

The overall optical transfer function (OTF) for an imaging system in optical turbulence may be modeled to include the atmospheric OTF and the diffraction OTF. This is given by

Eq. (5)

H (ρ) = H_{atm} (ρ) H_{dif} (ρ),

where

ρ = \sqrt{u^{2} + v^{2}}

and

u

and

v

are the spatial frequencies in units of cycles per unit distance. The atmospheric OTF is given by¹

Eq. (6)

H_{atm} (ρ) = \exp {- 3.44 {(\frac{λ l ρ}{r_{0}})}^{5 / 3} [1 - α {(\frac{λ l ρ}{D})}^{1 / 3}]},

where

l

is the camera focal length and

D

is the aperture diameter. The parameter

α

here serves as a switch. In particular, the long exposure atmospheric OTF is generated by Eq. (6) when

α = 0

, and the average short exposure OTF is given by Eq. (6) when

α = 1

. The diffraction limited OTF for a circular aperture is given by¹⁹

Eq. (7)

H_{dif} (ρ) = {\begin{cases} \frac{2}{π} [\cos^{- 1} (\frac{ρ}{2 ρ_{c}}) - \frac{ρ}{2 ρ_{c}} \sqrt{1 - {(\frac{ρ}{2 ρ_{c}})}^{2}}], & ρ \leq ρ_{c} \\ 0, & otherwise \end{cases},

where

ρ_{c} = 1 / (λ f_{n})

is the optical cut-off frequency, and the

f

-number is

f_{n} = l / D

. Note that the theoretical atmospheric PSF is found by applying the inverse Fourier transform to Eq. (5).

2.3.

Anisoplanatic Tilt Statistics

While the isoplanatic angle provides information about the level of anisoplanatism on a small scale, it does not provide insight into anisoplanatic behavior from points in the object plane with a large separation angle. Below, we describe two statistics that do capture large-scale anisoplanatic behavior. These are the two-axis $Z$ -tilt correlation and the DTV. We shall use these as key validation metrics for the simulation.

To begin, let us define a two-axis $Z$ -tilt vector as $α (θ) = {[α_{x} (θ), α_{y} (θ)]}^{T}$ for a source viewed from the direction angle vector $θ = {[θ_{x}, θ_{y}]}^{T}$ . For a spherical wave characterized by the Kolmogorov power spectrum, an analytical expression for the $Z$ -tilt correlation has been derived by Basu (now Bose-Pillai) et al.,¹⁶ following techniques outlined by Fried²⁰ and Winick.²¹ The tilt correlation can be expressed as a weighted integral of $C_{n}^{2} (z)$ as follows:

Eq. (8)

⟨ α (θ_{1}) \cdot α (θ_{2}) ⟩ = \int_{z = 0}^{L} C_{n}^{2} (z) f_{c} (z) d z .

The path weighting function is given by

Eq. (9)

f_{c} (z) = (- \frac{2.91}{2}) {(\frac{16}{π})}^{2} D^{- 1 / 3} \int_{φ = 0}^{2 π} \int_{u = 0}^{1} [u \cos^{- 1} (u) - u^{2} (3 - 2 u^{2}) \sqrt{1 - u^{2}}] \times {[{(\frac{u z}{L})}^{2} + {(\frac{L - z}{D} Δ θ)}^{2} + 2 (\frac{u z}{L}) (\frac{L - z}{D} Δ θ) \cos (φ)]}^{5 / 6} d u d φ,

where

Δ θ = | θ_{1} - θ_{2} |

.

Figure 2 shows a plot of $f_{c} (z)$ for different source separations, expressed in terms of pixel spacings. Note that the term “Nyquist pixels” here refers to pixel spacings corresponding to spatial sampling at the Nyquist rate, relative to the diffraction-limited optical cut-off frequency. The aperture size, path length, and Nyquist pixel spacing used in the evaluation are the same as those used in the simulations and are listed in Table 1. Note that the weight is maximum at the camera ( $z = 7000 m$ ) and drops down to zero at the source ( $z = 0 m$ ). This implies that the turbulence near the source does not contribute to tilt correlation seen at the camera. It is also apparent from Fig. 2 that the tilt correlations decrease with increasing angular separation between the sources.

Fig. 2

Weighting function for $C_{n}^{2} (z)$ profiles for anisoplanatic tilt angle correlation¹⁶ for $L = 7 km$ and $D = 0.2034 m$ . The angle difference between the two sources, $Δ θ = | θ_{1} - θ_{2} |$ , is expressed here in terms of Nyquist sampled pixel spacings for a focal length of $l = 1.2 m$ and wavelength of $λ = 0.525 μ m$ .

Table 1

Optical parameters.

Parameter	Value
Aperture	$D = 0.2034 m$
Focal length	$l = 1.2 m$
$F$ -number	$f / # = 5.9$
Wavelength	$λ = 0.525 μ m$
Object distance	$L = 7 km$
Nyquist pixel spacing (focal plane)	$δ_{f} = 1.5488 μ m$
Nyquist pixel spacing (object plane)	$δ_{o} = 9.0344 mm$

Let us now consider the DTV. This statistic is defined as

Eq. (10)

⟨ {[α (θ_{1}) - α (θ_{2})]}^{2} ⟩ = ⟨ {[α_{x} (θ_{1}) - α_{x} (θ_{2})]}^{2} + {[α_{y} (θ_{1}) - α_{y} (θ_{2})]}^{2} ⟩ = 2 [⟨ α^{2} (θ_{1}) ⟩ - ⟨ α (θ_{1}) \cdot α (θ_{2}) ⟩] .

We can derive an expression for the DTV using Eqs. (8) and (9). To do this, note that the tilt variance term in Eq. (10)

⟨ α^{2} (θ_{1}) ⟩

can be found by substituting

Δ θ = 0

in Eq. (9). This yields

Eq. (11)

⟨ α^{2} (θ_{1}) ⟩ = ⟨ α^{2} (θ_{2}) ⟩ = \int_{z = 0}^{L} C_{n}^{2} (z) f_{s} (z) d z,

where

Eq. (12)

f_{s} (z) = (- \frac{2.91}{2}) {(\frac{16}{π})}^{2} D^{- 1 / 3} (2 π) \int_{u = 0}^{1} [(u \cos^{- 1} u) - u^{2} (3 - 2 u^{2}) \sqrt{1 - u^{2}}] {(\frac{u z}{L})}^{5 / 3} d u .

Note that this yields numerical values closely matching those from the

Z

-tilt variance equation from Tyson,²² expressed in terms of

r_{0}

as

Eq. (13)

T_{Z}^{2} = . 3641 {(\frac{D}{r_{0}})}^{5 / 3} {(\frac{λ}{D})}^{2} .

Substituting the tilt variance from Eq. (11), and the tilt correlation from Eq. (8), into Eq. (10), the differential

Z

-tilt variance can be expressed as

Eq. (14)

⟨ {[α (θ_{1}) - α (θ_{2})]}^{2} ⟩ = \int_{z = 0}^{L} C_{n}^{2} (z) f_{d} (z) d z,

where the weighting function is given by

Eq. (15)

f_{d} (z) = 2.91 {(\frac{16}{π})}^{2} D^{- 1 / 3} \int_{φ = 0}^{2 π} \int_{u = 0}^{1} [u \cos^{- 1} (u) - u^{2} (3 - 2 u^{2}) \sqrt{1 - u^{2}}] \times {{[{(\frac{u z}{L})}^{2} + {(\frac{L - z}{D} Δ θ)}^{2} + 2 (\frac{u z}{L}) (\frac{L - z}{D} Δ θ) \cos (φ)]}^{5 / 6} - {(\frac{u z}{L})}^{5 / 3}} d u d φ .

Figure 3 shows a plot of $f_{d} (z)$ for different source separations. The aperture size and path length used in the evaluation are same as those in Fig. 2. Note that the weighting functions in Fig. 3 drop down to zero at both ends of the path. The zero weight at the source occurs because we have a point source with spherical wave emanating from it. Angular variation due to turbulence immediately at the source is effectively like rotating the point source. This just directs a different ray emanating from the point source to the observer. It does not change the angle of arrival of the point source observed by the camera. On the other side of the path, tilts due to turbulence near the camera tend to be very similar across the field of view, owing to the convergence of the optical paths near the camera. This causes the differential signal to drop to zero at the camera end. It is also evident from Fig. 2 that the DTV grows with increasing angular separation between the sources. This is because the tilt correlation drops and the DTV approaches $2 \times$ the tilt variance from Eq. (10). Furthermore, it is interesting to compare Fig. 3 to Fig. 1. Note that for small separations, the weighting function in Fig. 3 is weighted more heavily toward the source. This is also the case for the isoplanatic angle weighting in Fig. 1, which relates to small scale anisoplanatism. For large separations, the weighting in Fig. 3 appears to approach the $r_{0}$ weighting in Fig. 1. This is a satisfying result, given that $r_{0}$ governs tilt variance as seen in Eq. (13).

Fig. 3

Weighting function for $C_{n}^{2} (z)$ profiles for differential $Z$ -tilt variance¹⁶ for $L = 7 km$ and $D = 0.2034 m$ . The angle difference between the two sources $Δ θ = | θ_{1} - θ_{2} |$ is expressed here in terms of Nyquist sampled pixel spacings for a focal length of $l = 1.2 m$ and wavelength of $λ = 0.525 μ m$ .

3. Simulation Description

3.1.

Overview of Simulation

As mentioned in Sec. 1, the proposed method is based on the method of Bos and Roggemann.⁷ Extended phase screens are generated as shown in Fig. 4. Points in the object plane are projected to the center of the camera pupil. Two such examples are shown in Fig. 4. The local phase screens are cropped from the extended phase screens within a specified distance of the optical path for each point. These local phase screen portions are shown with the blue and green squares in Fig. 4. The extended phase screen sizes are determined based on the cropped phase screen size and the object size, as shown in Fig. 4. More will be said about these dimensions shortly.

Fig. 4

Phase screen geometry for anisoplanatic PSF generation, similar to that used by Bos and Roggemann.⁷

A simulated point source is numerically propagated from the source to the pupil plane though the cropped phase screens for up to each point in the object. The grid of object points is spaced according to the Nyquist sample spacing in the object plane. Nyquist spacing in the focal plane is given by $δ_{f} = 1 / (2 ρ_{c})$ , and the Nyquist spacing at the object is $δ_{o} = λ L / (2 D)$ . Our simulation includes a skip parameter, whereby we have the option to skip a specified number of Nyquist samples for the purposes of propagation, and then we interpolate the PSFs that are generated to obtain the complete set. A bilinear interpolation is employed here for each sample of the PSF, based on the samples from the 4 PSFs that surround the PSF being interpolated. Once the PSFs are all generated, the output image is generated by using PSFs as weights in a spatially varying weighted sum of pixels from an ideal image.⁷

The propagation method described in Sec. 3.2 is applied to each of a grid of points in the object plane. The point source wave functions are propagated through the phase screens and to the pupil plane. The incoherent PSF is then computed as described in Sec. 3.3. The phase screens are generated with the appropriate statistics as described in Sec. 3.4. Finally, in Sec. 3.5, we describe how the individual phase screen Fried parameters are determined in our simulation.

3.2.

Numerical Wave Propagation

The split-step propagation method used to form each PSF is illustrated in Fig. 5. It involves a point source and a series of $N$ phase screens that have been cropped based on the geometry shown in Fig. 4. Note that the path from the point in the object plane to the center of the camera aperture forms the centers of the cropped phase screen windows. We use nearest-neighbor interpolation here to speed computation. However, other forms of interpolation may be used.

Fig. 5

Illustration of the split-step propagation used to generate each individual PSF. The local phase screens shown are cropped from larger screens based on the point source location in the object plane, as shown in Fig. 4.

The real valued wave function for the propagating field in the $i$ ’th $(x, y)$ plane along the $z$ -axis is given by

Eq. (16)

w_{i} (x, y, t) = Re {u_{i} (x, y) e^{j 2 π ν t}} = | u_{i} (x, y) | \cos [2 π ν t + ∠ u_{i} (x, y)],

where

ν

is the frequency of the source. The complex amplitude phasor of the wave function in the

i

’th plane is given by

Eq. (17)

u_{i} (x, y) = | u_{i} (x, y) | e^{j ∠ u_{i} (x, y)} .

Our spherical wave propagation begins with a point source

u_{0} (x, y)

. We model this as a 2-D Gaussian windowed sinc function with quadratic phase.² This is given by

Eq. (18)

u_{0} (x, y) = λ L α^{2} e^{- \frac{j k}{2 L} (x^{2} + y^{2})} sinc (α x, α y) e^{- \frac{α^{2}}{16} (x^{2} + y^{2})},

where

α = \tilde{D} / (λ L)

, and

\tilde{D}

is the width of the amplitude field at the pupil plane where we want a uniform amplitude. The parameter

\tilde{D}

should be greater than

D

and less than phase screen width. This point source is designed to produce a flat amplitude over a

\tilde{D} \times \tilde{D}

patch at the pupil plane.²

We can express the free-space propagation, plus phase screen perturbation,² as

Eq. (19)

u_{i} (x, y) = [u_{i - 1} (x, y) * h_{Δ z_{i}} (x, y)] e^{j φ_{i} (x, y)},

for

i = 1,2, \dots, N

. The term

h_{Δ z_{i}} (x, y)

is the impulse response of free-space propagation² for Fresnel diffraction. This is given by

Eq. (20)

h_{Δ z_{i}} (x, y) = \frac{e^{j k Δ z_{i}}}{j λ Δ z_{i}} \exp [\frac{j k}{2 Δ z_{i}} (x^{2} + y^{2})] .

The term

φ_{i} (x, y)

in Eq. (19) is the phase imparted from the

i

’th phase screen realization. Note that we implement Eq. (19) on discrete space data using fast Fourier transforms (FFTs).² In particular, the method described by Schmidt as angular-spectrum propagation is used.² Absorbing borders using a Gaussian window may be applied in simulations in which a significant amount of signal energy reaches the borders of the simulation area.²

The sample spacing used to represent the phase screens and to implement Eq. (19) is based on Voelz’s critical sampling rule (best use of bandwidth and spatial support).²³ This gives a sample spacing for the phase screens of

Eq. (21)

Δ x = {(\frac{λ L}{N})}^{1 / 2},

where

N

is the number of samples used in each dimension. Note that this sample spacing is for the phase screens for the purposes of propagation and should not be confused with the Nyquist sample spacing for the camera. The cropped phase screen width is set to be a multiple of the aperture.² In particular, we use a phase screen width of

X = s D = Δ x \times N

, where

s

is the multiplier parameter. Using Eq. (21), this means that

N = {(s D)}^{2} / (λ L)

. We chose

s

to be as small as possible such that

s \geq 4

and

N

is a power of two (to speed up FFT computations). We also use a constant screen spacing so that all propagations can be achieved with the same impulse response in Eq. (20).

3.3.

Incoherent Point Spread Function

Following the propagations, the complex amplitude at the pupil plane $u_{N} (x, y)$ is obtained. This is multiplied by the camera aperture mask $a (x, y)$ , and a collimation-type phase compensation is used to allow the lens operation to focus the image at the focal length. This yields

Eq. (22)

p (x, y) = a (x, y) u_{N} (x, y) \exp [\frac{- j π (x^{2} + y^{2})}{λ L}] .

The incoherent PSF can be found based on Eq. (22) and is given by¹⁹

Eq. (23)

h (x, y) = {({| P (u, v) |}^{2}) |}_{u = \frac{x}{λ l}, v = \frac{y}{λ l}},

where

P (u, v) = FT {p (x, y)}

and

FT {\cdot}

is the Fourier transform. In practice, the PSF is evaluated discretely using the FFT and then resampled to the focal plane Nyquist spacing for the camera. We also normalize the discrete PSFs to all have a sum of 1.⁷ The final simulation output is computed with a spatially varying weighted sum as

Eq. (24)

z (m, n) = \sum_{k} \sum_{l} o (m - k, n - l) h_{m, n} (k, l),

where

o (m, n)

is the ideal discrete object image, and

h_{m, n} (k, l)

is the Nyquist sampled discrete PSF associated with object sample

(m, n)

.

Note that using the collimation in Eq. (22), the PSF image is focused at the focal length and not the image distance. This creates a magnification of $M_{sim} = - l / L$ , where the in-focus modeled system would have a magnification of $M_{\mod} = l / (l - L)$ . For large ranges, the difference is negligible. For shorter distances, the magnification can be corrected during the resampling step.

3.4.

Generating Phase Screen Realizations

Phase screen realizations are designed to follow a modified von Kármán phase power spectral density (PSD).² This PSD includes the Kolmogorov PSD as a special case, but has additional parametric flexibility. The modified von Kármán PSD is given by

Eq. (25)

S_{φ_{i}}^{m v K} (ρ) = \frac{0.023 e^{- ρ^{2} / ρ_{m}^{2}}}{r_{0_{i}}^{5 / 3} {(ρ^{2} + ρ_{0}^{2})}^{11 / 6}},

where

ρ = \sqrt{u^{2} + v^{2}}

,

ρ_{m} = 5.92 / (2 π l_{0})

, and

ρ_{0} = 1 / L_{0}

. Note that the radial spatial frequency

ρ

is measured in units of cycles/meter, and

r_{0_{i}}

is the

i

’th phase screen Fried parameter. The parameters

l_{0}

and

L_{0}

are the inner and outer turbulence scales, respectively.² Note that for

l_{0} = 0

and

L_{0} = \infty

, the modified von Kármán PSD is equivalent to a Kolmogorov PSD.² An example of the modified von Kármán PSD is shown in Fig. 6 for

l_{0} = . 01 m

,

L_{0} = 100 m

, and

r_{0_{i}} = 0.01 m

.

Fig. 6

Modified von Kármán phase PSD with $l_{0} = . 01 m$ , $L_{0} = 100 m$ , and $r_{0_{i}} = 0.01 m$ .

Generating a realization of a random process with specified PSD can be done by generating white noise with constant PSD of 1 and filtering it with a frequency response that is the square root of the desired PSD. Thus, we begin realizing the modified von Kármán phase screens by generating an $\tilde{N} \times \tilde{N}$ array of independent and identically distributed Gaussian random samples with standard deviation of 1. Note that $\tilde{N}$ is the extended phase screen dimension, prior to cropping to a size of $N$ for propagation. We wish for these samples to correspond to a constant PSD value of 1 over the range $- f_{s} / 2$ to $f_{s} / 2$ , where $f_{s} = 1 / Δ x$ . This means the total power (i.e., the integral of the PSD) should be $f_{s}^{2}$ . For our discrete samples, that means the variance should be $f_{s}^{2}$ . Thus, we multiply the unit standard deviation samples by $f_{s}$ . Next, we filter the array with a discrete-space impulse invariant²⁴ version of a filter with frequency response given by

Eq. (26)

H_{i} (u, v) = \sqrt{S_{φ_{i}}^{m v K} (u, v)},

for

- f_{s} / 2 \leq u

,

v \leq f_{s} / 2

. Our method avoids the need for generating complex valued random samples as used by some other methods. All phase screens are computed to the size of the maximum screen required and then cropped to the needed dimensions. This gives us consistent spectral content for all screens, regardless of the size required by the geometry shown in Fig. 4.

One can see in Fig. 6 that the most dynamic part of the modified von Kármán PSD lies at low spatial frequencies. Faithful generation of low spatial frequency content is essential in matching theoretical statistics such as tilt variance and the long-exposure OTF. Since the sampling process limits the resolution at which the spatial frequencies can be evaluated in Eqs. (25) and (26), special subharmonic methods may be required.²^,²⁵ Note that evaluating Eqs. (25) and (26) for use with FFTs is done in frequency increments of $1 / (N Δ_{x})$ . The subharmonic methods seek to produce appropriate spectral content at spatial frequencies below the $1 / (N Δ_{x})$ limit. We employ a subharmonic generation method based on the technique presented by Schmidt.² However, we use a real-valued 2-D Fourier series representation of the subharmonic spatial phase realizations, made up of a sum of 2-D cosines, to simplify the computations.

An example phase screen realization is shown in Fig. 7. This has been generated from the modified von Kármán PSD shown in Fig. 6. The scale is in units of radians, and an aperture of size $D = 0.2034 m$ is shown for reference. A single PSF generated from this phase screen and aperture is shown in Fig. 8. Note that the PSF is down sampled to the Nyquist pixel spacing for the camera model, prior to applying it as part of the spatially varying weighted sum operation in Eq. (24) to produce the output image.

Fig. 7

Example phase screen generated with the modified von Kármán phase PSD from Fig. 6. The scale is in units of radians, and an aperture of size $D = 0.2034 m$ is shown for reference.

Fig. 8

Example of a single PSF generated from the single modified von Kármán phase screen and aperture shown in Fig. 7. Note that the PSF is down sampled to the Nyquist pixel spacing for the camera model prior to applying it as part of the spatially varying weighted sum operation to produce the output image.

3.5.

Determining Phase Screen Parameters

The individual phase screen Fried parameters are determined such that they are consistent with the global Fried parameter, log-amplitude variance, and isoplanatic angle. This approach is based on the method presented by Schmidt,² but we have extended this to also include the isoplanatic angle. Based on the weightings shown in Fig. 1, we see that these three parameters give us a good balance of weighting along the optical path. This allows us to simulate variable $C_{n}^{2} (z)$ profiles more accurately.

Following Schmidt’s method,² the individual screen Fried parameters are related to screen $C_{n_{i}}^{2}$ parameters using

Eq. (27)

r_{0_{i}} = {[0.423 k^{2} C_{n_{i}}^{2} Δ z_{i}]}^{- 3 / 5},

or equivalently

Eq. (28)

C_{n_{i}}^{2} = \frac{r_{0_{i}}^{- 5 / 3}}{0.423 k^{2} Δ z_{i}} .

Using a discrete version of the Fried parameter from Eq. (2), and substituting using Eq. (28), we obtain

Eq. (29)

{\hat{r}}_{0, s w} = {[0.423 k^{2} \sum_{i = 1}^{N} (C_{n_{i}}^{2}) {(\frac{z_{i}}{L})}^{5 / 3} Δ z_{i}]}^{- 3 / 5} = {[\sum_{i = 1}^{N} r_{0_{i}}^{- 5 / 3} {(\frac{z_{i}}{L})}^{5 / 3}]}^{- 3 / 5} .

A similar approach using the isoplanatic angle from Eq. (3), yields

Eq. (30)

{\hat{θ}}_{0} = {[2.91 k^{2} L^{5 / 3} \sum_{i = 1}^{N} (C_{n_{i}}^{2}) {(1 - \frac{z_{i}}{L})}^{5 / 3} Δ z_{i}]}^{- 3 / 5} = {[\sum_{i = 1}^{N} L^{5 / 3} {(1 - \frac{z_{i}}{L})}^{5 / 3} 6.8794 r_{0_{i}}^{- 5 / 3}]}^{- 3 / 5} .

Finally, repeating the process for the log-amplitude variance in Eq. (4), yields

Eq. (31)

{\hat{σ}}_{χ, s w}^{2} = 0.563 k^{7 / 6} L^{5 / 6} \sum_{i = 1}^{N} C_{n_{i}}^{2} {(\frac{z_{i}}{L})}^{5 / 6} {(1 - \frac{z_{i}}{L})}^{5 / 6} Δ z_{i} = 1.331 k^{- 5 / 6} L^{5 / 6} \sum_{i = 1}^{N} r_{0_{i}}^{- 5 / 3} {(\frac{z_{i}}{L})}^{5 / 6} {(1 - \frac{z_{i}}{L})}^{5 / 6} .

Combining Eqs. (29)–(31), we have a set of under-determined linear equations when

N > 3

. These can be written in matrix notation as

Eq. (32)

[\begin{matrix} {({\hat{r}}_{0, s w})}^{- 5 / 3} \\ \frac{{\hat{σ}}_{χ, s w}^{2}}{1.331 k^{- 5 / 6} L^{5 / 6}} \\ \frac{{({\hat{θ}}_{0})}^{- 5 / 3}}{6.8794 L^{5 / 3}} \end{matrix}] = [\begin{matrix} {(\frac{z_{1}}{L})}^{5 / 3} & {(\frac{z_{2}}{L})}^{5 / 3} & \dots & {(\frac{z_{N}}{L})}^{5 / 3} \\ {(\frac{z_{1}}{L})}^{5 / 6} {(1 - \frac{z_{1}}{L})}^{5 / 6} & {(\frac{z_{2}}{L})}^{5 / 6} {(1 - \frac{z_{2}}{L})}^{5 / 6} & \dots & {(\frac{z_{N}}{L})}^{5 / 6} {(1 - \frac{z_{N}}{L})}^{5 / 6} \\ {(1 - \frac{z_{1}}{L})}^{5 / 3} & {(1 - \frac{z_{2}}{L})}^{5 / 3} & \dots & {(1 - \frac{z_{N}}{L})}^{5 / 3} \end{matrix}] [\begin{matrix} r_{0_{1}}^{- 5 / 3} \\ r_{0_{2}}^{- 5 / 3} \\ ⋮ \\ r_{0_{N}}^{- 5 / 3} \end{matrix}] .

We use a constrained least squares optimization, based on minimizing ${‖ \hat{b} - b ‖}^{2}$ , where $b = {[{(r_{0, s w})}^{- 5 / 3}, \frac{σ_{χ, s w}^{2}}{1.331 k^{- 5 / 6} L^{5 / 6}}, \frac{(θ_{0})^{- 5 / 3}}{6.8794 L^{5 / 3}}]}^{T}$ , and $\hat{b}$ is the left side of Eq. (32). Our implementation uses the MATLAB function “fmincon” to perform the minimization. The constraints are user controlled in our simulation tool. For the results presented here, we use two constraints. One constraint is that the log-amplitude variance generated by any phase screen cannot exceed 20% of the total. The other constraint is that the $r_{0_{N}} = \infty$ (i.e., no screen at the pupil). Since all point sources converge to the center of the pupil plane, they all share the exact same pupil phase screen. This tends to create excess tilt correlation. This might be expected by looking at the weighting in Fig. 2. Note that this is an artifact of the discrete phase screen approach. With more screens, the common pupil plane screen will represent a smaller $Δ z$ and cause less excess correlation. However, by setting the final phase screen phases to zero and using Eq. (32), we find that we are able to achieve a good match with the theoretical predictions, even with a relatively small number of screens.

We initialize the minimization search with Eq. (27), where $C_{n_{i}}^{2}$ are the average $C_{n}^{2} (z)$ for the section of the $C_{n}^{2} (z)$ profile represented by the $i$ ’th screen, at $z = z_{i}$ . We found that we are able to get reasonable results using these initial values for the simulation (especially for a large number of screens). However, we observed a better overall match with the validation metrics using the optimization.

4. Experimental Results and Validation

The experimental results presented have been generated using the optical parameters listed in Table 1 and simulation parameters in Table 2. Note that we are simulating a range of 7 km for a visible wavelength telescope and camera with a wavelength of $λ = 0.525 μ m$ . We have elected to use $N = 10$ phase screens (9 nonzero phase screens). We have chosen a very large outer scale, since we are validating against the Kolmogorov-based statistics in Sec. 2. Thus, for these results, the modified von Kármán PSD is acting essentially as a Kolmogorov PSD. The images are of size $257 \times 257 pixels$ , and spatial sampling at the Nyquist frequency is used (based on the diffraction-limited optical cut-off frequency). Note that other degradations such as downsampling and noise can easily be added after the turbulence simulation. We use a pixel skip parameter of four pixels to speed up the processing by a factor of $16 \times$ , compared with generating one PSF for each pixel. This can be adjusted according to the needs of the simulation and processing resources available. Our first set of results, in Sec. 4.1, is for a constant $C_{n}^{2} (z)$ profile, and then we consider varying profiles in Sec. 4.2. For each simulation, 200 frames are generated, and the theoretical validation metrics are compared with the parameters estimated from the simulation. Implementation issues and run times are discussed in Sec. 4.3.

Table 2

Simulation parameters.

Parameter	Value
Path length	$L = 7 km$
Propagation step	$Δ z = 700 m$
Cropped screen samples	$N = 256$
Propagation screen width	$X = 0.9699 m$
Pupil plane point spread	$\tilde{D} = 0.8136 m$
Propagation sample spacing	$Δ x = 0.0038 m$
Number of phase screens	$N = 10$ (9 nonzero)
Phase screen type	Modified von Kármán with subharmonics
Inner scale	$l_{0} = 0.01 m$
Outer scale	$L_{0} = 300 m$
Image size (pixels)	$257 \times 257 pixels$
Image size (object plane)	$2.3218 \times 2.3218 m$
Pixel skip	4 pixels ( $65 \times 65$ PSF array)

4.1.

Constant $C_{n}^{2} (z)$ Profile

We have simulated a constant $C_{n}^{2} (z)$ profile with five levels of turbulence, ranging from $C_{n}^{2} = 0.1 \times 10^{- 15} m^{- 2 / 3}$ to $C_{n}^{2} = 1.5 \times 10^{- 15} m^{- 2 / 3}$ . Some of the validation results are presented in Table 3. The validation metrics listed are $r_{0}$ , $θ_{0}$ , and the root-mean-squared (RMS) $Z$ -tilt. The Fried parameter is estimated from the simulation long-exposure PSF, obtained by averaging the generated PSFs over frames. The isoplanatic angle is estimated by an analysis of the wave function phases over the aperture for point sources of varying angular separations. Finally, the $Z$ -tilt is estimated by performing a correlation-based registration of the simulated PSFs. We have found that this type of correlation PSF registration corresponds well with $Z$ -tilt, and PSF centroids correspond well to $G$ -tilt. Note that we see a high level of agreement with regard to all of the metrics in Table 3.

Table 3

Constant Cn2 simulation results. Comparison of theoretical statistical parameters and those estimated from the simulation output.

Parameter	Cn2×10−15 (m−2/3)
Parameter	0.10	0.25	0.50	1.00	1.50
Theoretical $r_{0}$ (m)	0.1901	0.1097	0.0724	0.0478	0.0374
Simulation $r_{0}$ (m)	0.1897	0.1111	0.0740	0.0491	0.0387
Percent error (%)	$- 0.21$	1.28	2.21	2.72	3.48
Theoretical $θ_{0}$ ( $μ rads$ )	8.5401	4.9283	3.2515	2.1452	1.6819
Simulation $θ_{0}$ ( $μ rads$ )	8.6071	5.1919	3.3938	2.1933	1.6928
Percent error (%)	0.78	5.35	4.38	2.24	0.65
Theoretical $θ_{0}$ (pixels)	6.6170	3.8186	2.5193	1.6621	1.3032
Simulation $θ_{0}$ (pixels)	6.6689	4.0228	2.6296	1.6994	1.3116
Percent error (%)	0.78	5.35	4.38	2.24	0.64
Theoretical RMS $Z$ -tilt (pixels)	0.9026	1.4271	2.0182	2.8542	3.4957
Simulation RMS $Z$ -tilt (pixels)	0.9044	1.4294	2.0151	2.8398	3.4692
Percent error (%)	0.20	0.16	$- 0.15$	$- 0.50$	$- 0.76$

The phase screen Fried parameters, obtained using the procedure described in Sec. 3.5, are shown in Fig. 9 for two levels of turbulence. No value is plotted for the pupil location at $z = 7 km$ since it would be infinity, as described in Sec. 3.5. Note that even though we are simulating a constant $C_{n}^{2} (z)$ profile, the values in Fig. 9 are not constant. This is a consequence of the optimization procedure. One can see that the phase screen immediately before the pupil plane does have a lower $r_{0_{i}}$ than the others. This is in response to the pupil plane screen constraint, and this final nonzero screen is effectively “responsible” for a larger portion of the optical path than the others.

Fig. 9

Individual phase screen atmospheric coherence diameters $r_{0_{i}}$ used for phase screen generation for a constant refractive index structure parameter path of (a) $C_{n}^{2} = 0.25 \times 10^{- 15} m^{- 2 / 3}$ and (b) $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ .

The theoretical long- and short-exposure PSFs (with diffraction) are shown in Figs. 10 and 11, respectively. In particular, amplitude normalized cross sections of the theoretical and simulated PSFs are shown for two turbulence levels. An excellent agreement can be seen here in all cases. Note that without subharmonics, the simulated long-exposure PSFs tend to be too narrow, and the RMS $Z$ -tilts tend to be too low. The relationship between these two parameters is explored in depth by Hardie et al.¹⁵

Fig. 10

Comparison of theoretical and simulated long exposure PSFs for (a) $C_{n}^{2} = 0.25 \times 10^{- 15} m^{- 2 / 3}$ and (b) $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ .

Fig. 11

Comparison of theoretical and simulated short exposure (tilt corrected) PSFs for (a) $C_{n}^{2} = 0.25 \times 10^{- 15} m^{- 2 / 3}$ and (b) $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ .

Structure functions of phase² are shown in Figs. 12. These curves show the average squared phase difference over the aperture for two points separated by an angle of $Δ θ$ . Note that the isoplanatic angle is defined to be the angle where these structure functions have a value of 1 radian. These plots show fairly good agreement between the theory and simulation. Thus, the simulation appears to be capable of accurately capturing small scale anisoplanatic behavior. Note that the deviation of the simulated and theoretical curves for high source separation angles is due in phase wrapping. Phase wrapping is an unavoidable consequence of FFT processing during propagation. This artificially limits the estimated structure functions of phase in the simulation. It does not, however, compromise the generation of the PSFs in any way.

Fig. 12

Comparison of theoretical and simulated structure functions of phase, revealing the isoplanatic angle where $D_{φ} (Δ θ) = 1$ . (a) $C_{n}^{2} = 0.25 \times 10^{- 15} m^{- 2 / 3}$ and (b) $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ .

A comparison of theoretical and simulated $Z$ -tilt correlations is shown in Fig. 13 for two turbulence levels. A similar comparison of the DTV is shown in Fig. 14. These curves show that the simulation is accurately producing the larger scale anisoplanatic behavior predicted by the theoretical expressions. Note that if a nonzero phase screen is placed at the pupil plane, we tend to see simulated correlation that is higher, and DTV that is lower, than the theoretical values. This is explained by the fact that all PSFs share the exact same phase screen at the pupil, because of the converging optical paths.

Fig. 13

Comparison of theoretical and simulated $Z$ -tilt correlations for (a) $C_{n}^{2} = 0.25 \times 10^{- 15} m^{- 2 / 3}$ and (b) $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ . These are plotted versus $Δ θ$ , which is expressed in terms of Nyquist pixel spacings here for the parameters in Table 1.

Fig 14

Comparison of theoretical and simulated differential $Z$ -tilt variances for (a) $C_{n}^{2} = 0.25 \times 10^{- 15} m^{- 2 / 3}$ and (b) $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ . These are plotted versus $Δ θ$ , which is expressed in terms of Nyquist pixel spacings here for the parameters in Table 1.

Let us now examine some image results. The ideal object image used in these simulations is shown in Fig. 15. Degraded images for two levels of turbulence are shown in Fig. 16. Also shown in Fig. 16 are the corresponding $Z$ -tilt motion vectors, scaled by $2 \times$ . The level of blurring and warping clearly increases with $C_{n}^{2}$ .

Fig. 15

Ideal undegraded input image.

Fig. 16

Simulated output images and tilt maps for constant $C_{n}^{2}$ . (a) Output frame for $C_{n}^{2} = 0.25 \times 10^{- 15} m^{- 2 / 3}$ ; (b) output frame for $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ ; (c) tilt map for $C_{n}^{2} = 0.25 \times 10^{- 15} m^{- 2 / 3}$ ; (d) tilt map for $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ . Tilt vectors are scaled by $2 \times$ . All images are shown versus Nyquist pixel spacings. Video 1 shows a sequence of frames with 1 m/s wind-driven temporal correlation with $C_{n}^{2} = 0.25 \times 10^{- 15} m^{- 2 / 3}$ on the left and $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ on the right. (Video 1, MOV, 726 kb [URL: http://dx.doi.org/10.1117/1.OE.56.7.071502.1]).

4.2.

Variable $C_{n}^{2} (z)$ Profile

Validation analysis for the variable $C_{n}^{2} (z)$ profiles is shown in Table 4. The corresponding $C_{n}^{2} (z)$ profiles are shown in Fig. 17. Path A has heavy turbulence at the source, and Path B has heavy turbulence at the camera. The screen Fried parameters for these two cases are shown in Fig. 18. The long- and short-exposure PSFs are shown in Figs. 19 and 20, respectively. The structure functions of phase are shown in Fig. 21. The tilt correlation and DTV plots for the two paths are shown in Figs. 22 and 23, respectively. As with the constant $C_{n}^{2} (z)$ profile results, there is generally good agreement between theory and simulation. The most deviation is seen with the tilt correlation and DTV for Path A. We shall continue to investigate this in future work.

Table 4

Variable Cn2(z) profile simulation results. Comparison of theoretical statistical parameters and those estimated from the simulation output. Path average Cn2=1.00×10−15 m−2/3.

Parameter	Profile
Parameter	(A) Heavy at source	(B) Heavy at camera
Theoretical $r_{0}$ (m)	0.0687	0.0381
Simulation $r_{0}$ (m)	0.0703	0.0392
Percent error (%)	2.33	2.89
Theoretical $θ_{0}$ ( $μ rads$ )	1.7114	3.0861
Simulation $θ_{0}$ ( $μ rads$ )	1.7385	3.1816
Percent error (%)	1.58	3.09
Theoretical $θ_{0}$ (pixels)	1.3260	2.3912
Simulation $θ_{0}$ (pixels)	1.3470	2.4652
Percent error (%)	1.58	3.09
Theoretical RMS $Z$ -tilt (pixels)	2.1080	3.4455
Simulation RMS $Z$ -tilt (pixels)	2.1054	3.4071
Percent error (%)	$- 0.12$	$- 1.11$

Fig. 17

Variable $C_{n}^{2} (z)$ profiles. (a) Path A, heavy at source and (b) Path B, heavy at camera. The path averaged $C_{n}^{2}$ value for both paths is $C_{n}^{2} = 1.00 \times 10^{- 15} m^{- 2 / 3}$ .

Fig. 18

Individual phase screen atmospheric coherence diameters $r_{0_{i}}$ used for phase screen generation for variable refractive index structure parameter path. (a) Path A, heavy at source and (b) path B, heavy at camera.

Fig. 19

Comparison of theoretical and simulated long exposure PSFs for (a) Path A, heavy at source and (b) Path B, heavy at camera.

Fig. 20

Comparison of theoretical and simulated short exposure (tilt corrected) PSFs for (a) Path A, heavy at source and (b) Path B, heavy at camera.

Fig. 21

Comparison of theoretical and simulated structure functions of phase, revealing the isoplanatic angle (i.e., where $D_{φ} (Δ θ) = 1$ ). (a) Path A, heavy at source and (b) Path B, heavy at camera. These are plotted versus $Δ θ$ , which is expressed in terms of Nyquist pixel spacings here for the parameters in Table 1.

Fig. 22

Comparison of theoretical and simulated $Z$ -tilt correlations for (a) Path A, heavy at source and (b) Path B, heavy at camera. These are plotted versus $Δ θ$ , which is expressed in terms of Nyquist pixel spacings here for the parameters in Table 1.

Fig. 23

Comparison of theoretical and simulated differential $Z$ -tilt variances for (a) Path A, heavy at source and (b) Path B, heavy at camera. These are plotted versus $Δ θ$ , which is expressed in terms of Nyquist pixel spacings here for the parameters in Table 1.

Image results for the variable $C_{n}^{2} (z)$ profiles are shown in Fig. 24. The blurring is much more significant for Path B (heavy turbulence at camera). This is supported by the much smaller $r_{0}$ as shown in Table 4. The tilt vectors are also largest for Path B. This is consistent with the higher RMS $Z$ -tilt reported in Table 4. On the other hand, Path A yields a smaller isoplanatic angle.

Fig. 24

Simulated output images and tilt maps for variable $C_{n}^{2} (z)$ profiles. (a) Output frame for Path A, heavy at source; (b) output frame for Path B, heavy at camera; (c) tilt map for Path A; (d) tilt map for Path B. Tilt vectors are scaled by $2 \times$ . All images are shown versus Nyquist pixel spacings. Video 2 shows sequences of frames with $1 m / s$ wind-driven temporal correlation, with Path A on the left and Path B on the right. (Video 2, MOV, 778 kb [URL: http://dx.doi.org/10.1117/1.OE.56.7.071502.2]).

4.3.

Implementation

The simulation has been implemented in MATLAB 2016a and run on a PC with Intel(R) Xeon(R) CPU E5-2620 v3 at 2.40 GHz, 16 GB RAM, and running Windows 10. We employ a NVIDIA Quadro K-4200 GPU. We use the parallel computing toolbox for the PSF generation component, such that the FFTs are all computed on the GPU. This is done by simply assigning the corresponding numerical arrays to be “gpuArrays.” We have taken care to code the components efficiently by minimizing the use of “for” loops and maximizing the use of array operations. We used MATLAB’s profiler to look for bottle necks and made adjustments for improved processing speed. The algorithm component run times are listed in Table 5 using the parameters listed in Tables 1 and 2. Note that the GPU provides a $5.66 \times$ speed up of the PSF array computation. The algorithm processing times scales with the number of PSFs being generated. Processing larger images or using a smaller pixel skip will give a correspondingly larger run time.

Table 5

Algorithm average run times for processing one simulated 257×257 pixel frame using the parameters listed in Tables 1 and 2.

Algorithm component	Run time (s)
Phase screen generation	1.67
$65 \times 65$ PSF array generation (w/GPU)	20.46
$65 \times 65$ PSF array generation (w/o GPU)	115.73
Spatially varying convolution	2.23
Total (w/GPU)	24.36
Total (w/o GPU)	119.63

5. Conclusions

We have presented a numerical wave propagation method for simulating imaging of an extended scene under anisoplanatic conditions. In the simulation, we compute an array of PSFs for a 2-D grid of points on the object plane. An ideal image is then degraded by applying the PSFs in a spatially varying weighted sum operation. This gives rise to a spatially varying blurring and warping degradation. The PSFs are generated by the propagation of point source through an array of phase screens. The optical path for each point is somewhat different, by virtue of its originating position on the object plane and passes through different portions of the phase screens. This produces distinct but spatially correlated PSFs. We have employed a modified approach for determining the phase screen Fried parameters that incorporates $r_{0}$ , $θ_{0}$ and the log-amplitude variance.

To see if the resulting PSFs exhibit accurate anisoplanatic statistical properties, we have conducted an extensive validation analysis. This analysis shows that this simulation is capable of generating accurate anisoplanatic effects on both a small and large scale. Small scale anisoplanaticism is validated with the isoplanatic angle. Large scale anisoplanaticism is validated for the first time using tilt correlation,¹⁶ as well as with a newly derived DTV statistic for spherical waves. We also find excellent agreement between the simulated and theoretical long- and short-exposure PSFs.

We have demonstrated the simulation’s performance for both constant $C_{n}^{2} (z)$ profiles and varying profiles. We have also generated sequences of independent frames and sequences with temporally correlated frames. The temporal correlation is achieved by generating extended phase screens and translating these according to a specified wind speed. The portion of the screens between the object and camera for a given frame are used for that frame’s PSF generation. We have included video results showing the results of the temporally correlated sequence generation. Based on our analysis, we believe that this simulation approach can be used effectively to study anisoplanatic optical turbulence and to aid in the development of image restoration methods.

Acknowledgments

The authors would like to thank Dr. Craig Olson and Dr. Mike Theisen at L-3 Communications Cincinnati Electronics, who provided helpful input regarding wave propagation for the simulation. We would also like to thank Matthew D. Howard at AFRL for providing helpful suggestions related to sampling of the phase screens. Thanks also to Michael Rucci and Dr. Barry Karch for providing feedback regarding the simulation and the statistical validation. This work has been supported in part by funding from L-3 Communications Cincinnati Electronics and under AFRL Award No. FA8650-10-2-7028 and FA9550-14-1-0244. Approved for public release, Case Number 88ABW-2016-5066.

References

1.

M. C. Roggemann and B. M. Welsh, Imaging through Turbulence, Laser and Optical Science and Technology, CRC Press, Boca Raton, Florida (1996). Google Scholar

2.

J. D. Schmidt, Numerical Simulation of Optical Wave Propagation with Examples in MATLAB, SPIE Press, Bellingham, Washington (2010). Google Scholar

3.

J. Beck, C. Bekins and J. P. Bos, “Wavepy: a python package for wave optics,” Proc. SPIE, 9846 984603 (2016). http://dx.doi.org/10.1117/12.2225769 PSISDG 0277-786X Google Scholar

4.

C. J. Carrano, “Anisoplanatic performance of horizontal-path speckle imaging,” Proc. SPIE, 5162 14 –27 (2003). http://dx.doi.org/10.1117/12.508082 PSISDG 0277-786X Google Scholar

5.

R. W. Praus, K. Beardmore and B. P. Venet, “The ‘light tunneling’ optical propagation method in wavetrain,” MZA Associates Corporation Technical Memorandum, MZA Associates, Dayton, Ohio (2007). Google Scholar

6.

J. P. Bos and M. C. Roggeman, “Simulation of extended scenes imaged through turbulence over horizontal paths,” Proc. SPIE, 8161 816106 (2011). http://dx.doi.org/10.1117/12.893869 PSISDG 0277-786X Google Scholar

7.

J. P. Bos and M. C. Roggemann, “Technique for simulating anisoplanatic image formation over long horizontal paths,” Opt. Eng., 51 101704 (2012). http://dx.doi.org/10.1117/1.OE.51.10.101704 Google Scholar

8.

G. Monnier, F.-R. Duval and S. Amram, “GPU-based simulation of optical propagation through turbulence for active and passive imaging,” Proc. SPIE, 8535 85350D (2012). http://dx.doi.org/10.1117/12.981528 PSISDG 0277-786X Google Scholar

9.

G. Monnier, F.-R. Duval and S. Amram, “GPU-based simulation of optical propagation through turbulence for active and passive imaging,” Proc. SPIE, 9242 92421R (2014). http://dx.doi.org/10.1117/12.2067492 PSISDG 0277-786X Google Scholar

10.

E. Repasi and R. Weiss, “Computer simulation of image degradations by atmospheric turbulence for horizontal views,” Proc. SPIE, 8014 80140U (2011). http://dx.doi.org/10.1117/12.883805 PSISDG 0277-786X Google Scholar

11.

K. R. Leonard, J. Howe and D. E. Oxford, “Simulation of atmospheric turbulence effects and mitigation algorithms on stand-off automatic facial recognition,” Proc. SPIE, 8546 85460O (2012). http://dx.doi.org/10.1117/12.979480 PSISDG 0277-786X Google Scholar

12.

K. R. Leonard and R. L. Espinola, “Validation of atmospheric turbulence simulations of extended scenes,” Proc. SPIE, 9071 907118 (2014). http://dx.doi.org/10.1117/12.2054311 PSISDG 0277-786X Google Scholar

13.

J. P. Bos and M. C. Roggemann, “Mean squared error performance of speckle-imaging using the bispectrum in horizontal imaging applications,” Proc. SPIE, 8056 805603 (2011). http://dx.doi.org/10.1117/12.884093 PSISDG 0277-786X Google Scholar

14.

G. E. Archer, J. P. Bos and M. C. Roggemann, “Comparison of bispectrum, multiframe blind deconvolution and hybrid bispectrum-multiframe blind deconvolution image reconstruction techniques for anisoplanatic, long horizontal-path imaging,” Opt. Eng., 53 (4), 043109 (2014). http://dx.doi.org/10.1117/1.OE.53.4.043109 Google Scholar

15.

R. C. Hardie et al., “Multi-frame atmospheric optical turbulence mitigation with quantitative error analysis using simulated and real image sequences,” Opt. Eng., (2016). Google Scholar

16.

S. Basu, J. E. McCrae and S. T. Fiorino, “Estimation of the path-averaged atmospheric refractive index structure constant from time-lapse imagery,” Proc. SPIE, 9465 94650T (2015). http://dx.doi.org/10.1117/12.2177330 PSISDG 0277-786X Google Scholar

17.

F. Dios et al., “Temporal analysis of laser beam propagation in the atmosphere using computer-generated long phase screens,” Opt. Express, 16 2206 –2220 (2008). http://dx.doi.org/10.1364/OE.16.002206 OPEXFF 1094-4087 Google Scholar

18.

L.-K. Yu et al., “Isoplanatic angle in finite distance: experimental validation,” Opt. Eng., 54 (2), 024105 (2015). http://dx.doi.org/10.1117/1.OE.54.2.024105 Google Scholar

19.

J. W. Goodman, Introduction to Fourier Optics, 3rd ed.Roberts and Company Publishers, Greenwood Village, Colorado (2004). Google Scholar

20.

D. L. Fried, “Varieties of isoplanatism,” Proc. SPIE, 0075 20 –29 (1976). http://dx.doi.org/10.1117/12.954733 PSISDG 0277-786X Google Scholar

21.

K. A. Winick, “Stellar scintillation technique for the measurement of tilt anisoplanatism,” J. Opt. Soc. Am. A, 5 1929 –1936 (1988). http://dx.doi.org/10.1364/JOSAA.5.001929 JOAOD6 0740-3232 Google Scholar

22.

R. K. Tyson, Adaptive Optics Engineering Handbook, CRC Press, New York (1999). Google Scholar

23.

D. G. Voelz, Computational Fourier Optics: A MATLAB Tutorial, IOP Publishing Ltd., United Kingdom (2011). Google Scholar

24.

A. V. Oppenheim and R. W. Schafer, Discrete-time Signal Processing, 3rd ed.Prentice-Hall, Inc., Upper Saddle River, New Jersey (2010). Google Scholar

25.

R. G. Lane, A. Glindemann and J. C. Dainty, “Simulation of a Kolmogorov phase screen,” Waves Random Media, 2 (3), 209 –224 (1992). http://dx.doi.org/10.1088/0959-7174/2/3/003 WRMEEV 0959-7174 Google Scholar

Biography

Russell C. Hardie is a full professor in the Department of Electrical and Computer Engineering at the University of Dayton, with a joint appointment in the Department of Electro-Optics. He received the University of Dayton’s Top University-Wide Teaching Award, the 2006 Alumni Award in teaching, the Rudolf Kingslake Medal and Prize from SPIE in 1998 for super-resolution research, the School of Engineering Award of Excellence in teaching in 1999, and the first annual Professor of the Year Award in 2002 from the student chapter of the IEEE.

Jonathan D. Power received his BSEE degree from Cedarville University and his MSEE degree from the University of Dayton. He is an electrical engineer working at Wright-Patterson Air Force Base. For his master’s, he completed a thesis and developed a simulation for modeling anisoplanatic atmospheric turbulence for imagery along slanted optical paths. His research interests include atmospheric turbulence, polarimetric imaging, and image processing.

Daniel A. LeMaster is the technical advisor for the EO Target Detection and Surveillance Branch in the Sensors Directorate of the Air Force Research Laboratory (AFRL). He has served as a conference chair, committee member, and guest editor for SPIE and has taught as an adjunct professor (Wright State University) and as a professional education instructor (Georgia Tech Research Institute). Prior to AFRL, he worked as a research engineer in the defense intelligence community and served in the US Army.

Douglas R. Droege received his BS degree in electrical engineering, his BS degree in computer science, his MS degree in electrical engineering, and his PhD degree in electrical engineering, all from the University of Dayton. He is a director of Advanced Programs at L-3 Communications Cincinnati Electronics (L-3). He received the L-3 Integrated Sensor Systems Sector Leadership Award in 2011 and the Corporate L-3 Engineer of the Year Award in 2012.

Szymon Gladysz received his PhD degree in physics from the National University of Ireland in Galway. He is head of the Adaptive Optics Group at Fraunhofer Institute of Optronics, System Technologies and Image Exploitation, Ettlingen, Germany. He has authored or coauthored around 60 publications on turbulence and adaptive optics in journals and conference proceedings. Additionally, he represents Germany on NATO SET-226 Research Task Group [turbulence mitigation for electro optics (EO) and laser systems] and is a committee member on various OSA and SPIE conferences.

Santasri Bose-Pillai received her BSEE degree (with Hons.) from Jadavpur University, India, in 2000, her MSEE degree in 2005, and her PhD in electrical engineering with a focus in optics in 2008 from New Mexico State University. Currently, she is a senior research associate at AFIT’s Center for Directed Energy within Engineering Physics Department. Her research interests include propagation and imaging through atmospheric turbulence, telescope pointing and tracking, rough surface scattering, and laser communications. She is a member of OSA, SPIE, and DEPS.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Russell C. Hardie, Jonathan D. Power, Daniel A. LeMaster, Douglas R. Droege, Szymon Gladysz, and Santasri Bose-Pillai "Simulation of anisoplanatic imaging through optical turbulence using numerical wave propagation with new validation analysis," Optical Engineering 56(7), 071502 (6 February 2017). https://doi.org/10.1117/1.OE.56.7.071502

Received: 11 October 2016; Accepted: 7 November 2016; Published: 6 February 2017

Access the abstract

JOURNAL ARTICLE
16 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 57 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Point spread functions

Cameras

Optical turbulence

Computer simulations

Wave propagation

Device simulation

Optical imaging

1.

Introduction

2.

Anisoplanatic Optical Turbulence

2.1.

Optical Turbulence Statistics

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Fig. 1

2.2.

Optical Transfer Functions

Eq. (5)

Eq. (6)

Eq. (7)

2.3.

Anisoplanatic Tilt Statistics

Eq. (8)

Eq. (9)

Fig. 2

Table 1

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Fig. 3

3.

Simulation Description

3.1.

Overview of Simulation

Fig. 4

3.2.

Numerical Wave Propagation

Fig. 5

Eq. (16)

Eq. (17)

Eq. (18)

Eq. (19)

Eq. (20)

Eq. (21)

3.3.

Incoherent Point Spread Function

Eq. (22)

Eq. (23)

Eq. (24)

3.4.

Generating Phase Screen Realizations

Eq. (25)

Fig. 6

Eq. (26)

Fig. 7

Fig. 8

3.5.

Determining Phase Screen Parameters

Eq. (27)

Eq. (28)

Eq. (29)

Eq. (30)

Eq. (31)

Eq. (32)

4.

Experimental Results and Validation

Table 2

4.1.

Constant Cn2(z) Profile

Table 3

Fig. 9

Fig. 10

Fig. 11

Fig. 12

Fig. 13

Fig 14

Fig. 15

Fig. 16

4.2.

Variable Cn2(z) Profile

Table 4

Constant $C_{n}^{2} (z)$ Profile

Variable $C_{n}^{2} (z)$ Profile