## 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.^{1} 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.^{1}

The simulation of imaging under isoplanatic atmospheric turbulence has been well studied.^{1} Commercial and free open-source software^{2}^{,}^{3} 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^{4} and Praus et al.^{5} Further wave propagation-based anisoplanatic simulation works include that of Bos and Roggemann,^{6}^{,}^{7} and Monnier et al.^{8}^{,}^{9} Anisoplanatic simulations that do not involve wave propagation have also been developed.^{10}11.^{–}^{12} 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.^{4}^{,}^{13}14.^{–}^{15} 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.^{7} 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,^{7} 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.^{7} Our approach is based on the constrained least squares optimization presented by Schmidt,^{2} 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,^{16} 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.^{17} 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,^{1} this is given by

## (1)

$${D}_{n}(\mathbf{r})=\u27e8{[n(\mathbf{x}+\mathbf{r})-n(\mathbf{x})]}^{2}\u27e9={C}_{n}^{2}{r}^{2/3},$$^{1}The units of ${C}_{n}^{2}$ are ${\mathrm{m}}^{-2/3}$ and typical values tend to range from $1\times {10}^{-13}$ to $1\times {10}^{-17}$.

^{2}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)^{2} can be expressed as a weighted integral of ${C}_{n}^{2}(z)$ yielding

## (2)

$${r}_{0}={[0.423{k}^{2}{\int}_{z=0}^{z=L}{C}_{n}^{2}(z){\left(\frac{z}{L}\right)}^{5/3}\mathrm{d}z]}^{-3/5},$$Another important statistic is the isoplanatic angle.^{1} 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.^{2}^{,}^{18} 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^{2}^{,}^{18}

## (3)

$${\theta}_{0}={[2.91{k}^{2}{L}^{5/3}{\int}_{z=0}^{z=L}{C}_{n}^{2}(z){(1-\frac{z}{L})}^{5/3}\mathrm{d}z]}^{-3/5}.$$^{2}given by

## (4)

$${\sigma}_{\chi}^{2}=0.563{k}^{7/6}{L}^{5/6}{\int}_{z=0}^{z=L}{C}_{n}^{2}(z){\left(\frac{z}{L}\right)}^{5/6}{(1-\frac{z}{L})}^{5/6}\mathrm{d}z.$$## 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

where $\rho =\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^{1}

## (6)

$${H}_{\mathrm{atm}}(\rho )=\mathrm{exp}\{-3.44{\left(\frac{\lambda l\rho}{{r}_{0}}\right)}^{5/3}[1-\alpha {\left(\frac{\lambda l\rho}{D}\right)}^{1/3}\left]\right\},$$^{19}

## (7)

$${H}_{\mathrm{dif}}(\rho )=\{\begin{array}{ll}\frac{2}{\pi}\left[{\mathrm{cos}}^{-1}\right(\frac{\rho}{2{\rho}_{c}})-\frac{\rho}{2{\rho}_{c}}\sqrt{1-{\left(\frac{\rho}{2{\rho}_{c}}\right)}^{2}}],& \rho \le {\rho}_{c}\\ 0,& \text{otherwise}\end{array},$$## 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 $\mathbf{\alpha}(\mathbf{\theta})={[{\alpha}_{x}(\mathbf{\theta}),{\alpha}_{y}(\mathbf{\theta})]}^{T}$ for a source viewed from the direction angle vector $\mathbf{\theta}={[{\theta}_{x},{\theta}_{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.,^{16} following techniques outlined by Fried^{20} and Winick.^{21} The tilt correlation can be expressed as a weighted integral of ${C}_{n}^{2}(z)$ as follows:

## (8)

$$\u27e8\mathbf{\alpha}({\mathbf{\theta}}_{1})\xb7\mathbf{\alpha}({\mathbf{\theta}}_{2})\u27e9={\int}_{z=0}^{L}{C}_{n}^{2}(z){f}_{c}(z)\mathrm{d}z.$$## (9)

$${f}_{c}(z)=(-\frac{2.91}{2}){\left(\frac{16}{\pi}\right)}^{2}{D}^{-1/3}{\int}_{\phi =0}^{2\pi}{\int}_{u=0}^{1}[u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}(u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]\times {[{\left(\frac{uz}{L}\right)}^{2}+{\left(\frac{L-z}{D}\mathrm{\Delta}\theta \right)}^{2}+2(\frac{uz}{L}\left)\right(\frac{L-z}{D}\mathrm{\Delta}\theta )\mathrm{cos}(\phi )]}^{5/6}\mathrm{d}u\text{\hspace{0.17em}}\mathrm{d}\phi ,$$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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$) and drops down to zero at the source ($z=0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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.

## Table 1

Optical parameters.

Parameter | Value |
---|---|

Aperture | $D=0.2034\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Focal length | $l=1.2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

$F$-number | $f/\#=5.9$ |

Wavelength | $\lambda =0.525\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ |

Object distance | $L=7\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{km}$ |

Nyquist pixel spacing (focal plane) | ${\delta}_{f}=1.5488\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ |

Nyquist pixel spacing (object plane) | ${\delta}_{o}=9.0344\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ |

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

## (10)

$$\u27e8{[\mathbf{\alpha}({\mathbf{\theta}}_{1})-\mathbf{\alpha}({\mathbf{\theta}}_{2})]}^{2}\u27e9=\u27e8{[{\mathbf{\alpha}}_{x}({\mathbf{\theta}}_{1})-{\alpha}_{x}({\mathbf{\theta}}_{2})]}^{2}+{[{\alpha}_{y}({\mathbf{\theta}}_{1})-{\alpha}_{y}({\mathbf{\theta}}_{2})]}^{2}\u27e9=2[\u27e8{\mathbf{\alpha}}^{2}({\mathbf{\theta}}_{1})\u27e9-\u27e8\mathbf{\alpha}({\mathbf{\theta}}_{1})\xb7\mathbf{\alpha}({\mathbf{\theta}}_{2})\u27e9].$$## (11)

$$\u27e8{\mathbf{\alpha}}^{2}({\mathbf{\theta}}_{1})\u27e9=\u27e8{\mathbf{\alpha}}^{2}({\mathbf{\theta}}_{2})\u27e9={\int}_{z=0}^{L}{C}_{n}^{2}(z){f}_{s}(z)\mathrm{d}z,$$## (12)

$${f}_{s}(z)=(-\frac{2.91}{2}){\left(\frac{16}{\pi}\right)}^{2}{D}^{-1/3}(2\pi ){\int}_{u=0}^{1}[(u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]{\left(\frac{uz}{L}\right)}^{5/3}\mathrm{d}u.$$^{22}expressed in terms of ${r}_{0}$ as

## (13)

$${T}_{Z}^{2}=.3641{\left(\frac{D}{{r}_{0}}\right)}^{5/3}{\left(\frac{\lambda}{D}\right)}^{2}.$$## (14)

$$\u27e8{[\mathbf{\alpha}({\mathbf{\theta}}_{1})-\mathbf{\alpha}({\mathbf{\theta}}_{2})]}^{2}\u27e9={\int}_{z=0}^{L}{C}_{n}^{2}(z){f}_{d}(z)\mathrm{d}z,$$## (15)

$${f}_{d}(z)=2.91{\left(\frac{16}{\pi}\right)}^{2}{D}^{-1/3}{\int}_{\phi =0}^{2\pi}{\int}_{u=0}^{1}[u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}(u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]\times \{{[{\left(\frac{uz}{L}\right)}^{2}+{\left(\frac{L-z}{D}\mathrm{\Delta}\theta \right)}^{2}+2(\frac{uz}{L}\left)\right(\frac{L-z}{D}\mathrm{\Delta}\theta )\mathrm{cos}(\phi )]}^{5/6}-{\left(\frac{uz}{L}\right)}^{5/3}\}\mathrm{d}u\text{\hspace{0.17em}}\mathrm{d}\phi .$$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).

## 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.^{7} 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.

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 ${\delta}_{f}=1/(2{\rho}_{c})$, and the Nyquist spacing at the object is ${\delta}_{o}=\lambda L/(2D)$. 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.^{7}

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.

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

## (16)

$${w}_{i}(x,y,t)=\mathrm{Re}\{{u}_{i}(x,y){e}^{j2\pi \nu t}\}=|{u}_{i}(x,y)|\mathrm{cos}[2\pi \nu t+\angle {u}_{i}(x,y)],$$^{2}This is given by

## (18)

$${u}_{0}(x,y)=\lambda L{\alpha}^{2}{e}^{-\frac{jk}{2\text{\hspace{0.17em}\hspace{0.17em}}L}({x}^{2}+{y}^{2})}\text{\hspace{0.17em}}\mathrm{sinc}(\alpha x,\alpha y){e}^{-\frac{{\alpha}^{2}}{16}({x}^{2}+{y}^{2})},$$^{2}

We can express the free-space propagation, plus phase screen perturbation,^{2} as

^{2}for Fresnel diffraction. This is given by

## (20)

$${h}_{\mathrm{\Delta}{z}_{i}}(x,y)=\frac{{e}^{jk\mathrm{\Delta}{z}_{i}}}{j\lambda \mathrm{\Delta}{z}_{i}}\text{\hspace{0.17em}}\mathrm{exp}[\frac{jk}{2\mathrm{\Delta}{z}_{i}}({x}^{2}+{y}^{2})].$$^{2}In particular, the method described by Schmidt as angular-spectrum propagation is used.

^{2}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.

^{2}

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).^{23} This gives a sample spacing for the phase screens of

^{2}In particular, we use a phase screen width of $X=sD=\mathrm{\Delta}x\times \mathcal{N}$, where $s$ is the multiplier parameter. Using Eq. (21), this means that $\mathcal{N}={(sD)}^{2}/(\lambda L)$. We chose $s$ to be as small as possible such that $s\ge 4$ and $\mathcal{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

## (22)

$$p(x,y)=a(x,y){u}_{N}(x,y)\mathrm{exp}\left[\frac{-j\pi ({x}^{2}+{y}^{2})}{\lambda L}\right].$$^{19}where $P(u,v)=\mathrm{FT}\{p(x,y)\}$ and $\mathrm{FT}\{\xb7\}$ 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.

^{7}The final simulation output is computed with a spatially varying weighted sum as 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}_{\mathrm{sim}}=-l/L$, where the in-focus modeled system would have a magnification of ${M}_{\mathrm{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).^{2} 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

## (25)

$${S}_{{\phi}_{i}}^{mvK}(\rho )=\frac{0.023{\mathrm{e}}^{-{\rho}^{2}/{\rho}_{m}^{2}}}{{r}_{{0}_{i}}^{5/3}{({\rho}^{2}+{\rho}_{0}^{2})}^{11/6}},$$^{2}Note that for ${l}_{0}=0$ and ${L}_{0}=\infty $, the modified von Kármán PSD is equivalent to a Kolmogorov PSD.

^{2}An example of the modified von Kármán PSD is shown in Fig. 6 for ${l}_{0}=.01\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, ${L}_{0}=100\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, and ${r}_{{0}_{i}}=0.01\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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{\mathcal{N}}\times \tilde{\mathcal{N}}$ array of independent and identically distributed Gaussian random samples with standard deviation of 1. Note that $\tilde{\mathcal{N}}$ is the extended phase screen dimension, prior to cropping to a size of $\mathcal{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/\mathrm{\Delta}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^{24} version of a filter with frequency response given by

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.^{2}^{,}^{25} Note that evaluating Eqs. (25) and (26) for use with FFTs is done in frequency increments of $1/(\mathcal{N}{\mathrm{\Delta}}_{x})$. The subharmonic methods seek to produce appropriate spectral content at spatial frequencies below the $1/(\mathcal{N}{\mathrm{\Delta}}_{x})$ limit. We employ a subharmonic generation method based on the technique presented by Schmidt.^{2} 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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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.

## 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,^{2} 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,^{2} the individual screen Fried parameters are related to screen ${C}_{{n}_{i}}^{2}$ parameters using

## (29)

$${\widehat{r}}_{0,sw}={[0.423{k}^{2}\sum _{i=1}^{N}({C}_{{n}_{i}}^{2}){\left(\frac{{z}_{i}}{L}\right)}^{5/3}\mathrm{\Delta}{z}_{i}]}^{-3/5}={\left[\sum _{i=1}^{N}{r}_{{0}_{i}}^{-5/3}{\left(\frac{{z}_{i}}{L}\right)}^{5/3}\right]}^{-3/5}.$$## (30)

$${\widehat{\theta}}_{0}={[2.91{k}^{2}{L}^{5/3}\sum _{i=1}^{N}({C}_{{n}_{i}}^{2}){(1-\frac{{z}_{i}}{L})}^{5/3}\mathrm{\Delta}{z}_{i}]}^{-3/5}={\left[\sum _{i=1}^{N}{L}^{5/3}{(1-\frac{{z}_{i}}{L})}^{5/3}6.8794{r}_{{0}_{i}}^{-5/3}\right]}^{-3/5}.$$## (31)

$${\widehat{\sigma}}_{\chi ,sw}^{2}=0.563{k}^{7/6}{L}^{5/6}\sum _{i=1}^{N}{C}_{{n}_{i}}^{2}{\left(\frac{{z}_{i}}{L}\right)}^{5/6}{(1-\frac{{z}_{i}}{L})}^{5/6}\mathrm{\Delta}{z}_{i}=1.331{k}^{-5/6}{L}^{5/6}\sum _{i=1}^{N}{r}_{{0}_{i}}^{-5/3}{\left(\frac{{z}_{i}}{L}\right)}^{5/6}{(1-\frac{{z}_{i}}{L})}^{5/6}.$$## (32)

$$\left[\begin{array}{c}{({\widehat{r}}_{0,sw})}^{-5/3}\\ \frac{{\widehat{\sigma}}_{\chi ,sw}^{2}}{1.331{k}^{-5/6}{L}^{5/6}}\\ \frac{{({\widehat{\theta}}_{0})}^{-5/3}}{6.8794{L}^{5/3}}\end{array}\right]=\left[\begin{array}{cccc}{\left(\frac{{z}_{1}}{L}\right)}^{5/3}& {\left(\frac{{z}_{2}}{L}\right)}^{5/3}& \cdots & {\left(\frac{{z}_{N}}{L}\right)}^{5/3}\\ {\left(\frac{{z}_{1}}{L}\right)}^{5/6}{(1-\frac{{z}_{1}}{L})}^{5/6}& {\left(\frac{{z}_{2}}{L}\right)}^{5/6}{(1-\frac{{z}_{2}}{L})}^{5/6}& \cdots & {\left(\frac{{z}_{N}}{L}\right)}^{5/6}{(1-\frac{{z}_{N}}{L})}^{5/6}\\ {(1-\frac{{z}_{1}}{L})}^{5/3}& {(1-\frac{{z}_{2}}{L})}^{5/3}& \cdots & {(1-\frac{{z}_{N}}{L})}^{5/3}\end{array}\right]\left[\begin{array}{c}{r}_{{0}_{1}}^{-5/3}\\ {r}_{{0}_{2}}^{-5/3}\\ \vdots \\ {r}_{{0}_{N}}^{-5/3}\end{array}\right].$$We use a constrained least squares optimization, based on minimizing ${\Vert \hat{\mathbf{b}}-\mathbf{b}\Vert}^{2}$, where $\mathbf{b}={[{({r}_{0,sw})}^{-5/3},\frac{{\sigma}_{\chi ,sw}^{2}}{1.331{k}^{-5/6}{L}^{5/6}},\frac{({\theta}_{0}{)}^{-5/3}}{6.8794{L}^{5/3}}]}^{T}$, and $\hat{\mathbf{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 $\mathrm{\Delta}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 $\lambda =0.525\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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\text{\hspace{0.17em}\hspace{0.17em}}\text{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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{km}$ |

Propagation step | $\mathrm{\Delta}z=700\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Cropped screen samples | $\mathcal{N}=256$ |

Propagation screen width | $X=0.9699\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Pupil plane point spread | $\tilde{D}=0.8136\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Propagation sample spacing | $\mathrm{\Delta}x=0.0038\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Outer scale | ${L}_{0}=300\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Image size (pixels) | $257\times 257\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$ |

Image size (object plane) | $2.3218\times 2.3218\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$ to ${C}_{n}^{2}=1.5\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$. Some of the validation results are presented in Table 3. The validation metrics listed are ${r}_{0}$, ${\theta}_{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) | ||||
---|---|---|---|---|---|

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 ${\theta}_{0}$ ($\mu \mathrm{rads}$) | 8.5401 | 4.9283 | 3.2515 | 2.1452 | 1.6819 |

Simulation ${\theta}_{0}$ ($\mu \mathrm{rads}$) | 8.6071 | 5.1919 | 3.3938 | 2.1933 | 1.6928 |

Percent error (%) | 0.78 | 5.35 | 4.38 | 2.24 | 0.65 |

Theoretical ${\theta}_{0}$ (pixels) | 6.6170 | 3.8186 | 2.5193 | 1.6621 | 1.3032 |

Simulation ${\theta}_{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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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.

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.^{15}

Structure functions of phase^{2} are shown in Figs. 12. These curves show the average squared phase difference over the aperture for two points separated by an angle of $\mathrm{\Delta}\theta $. 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.

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.

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}$.

## 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 | |
---|---|---|

(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 ${\theta}_{0}$ ($\mu \mathrm{rads}$) | 1.7114 | 3.0861 |

Simulation ${\theta}_{0}$ ($\mu \mathrm{rads}$) | 1.7385 | 3.1816 |

Percent error (%) | 1.58 | 3.09 |

Theoretical ${\theta}_{0}$ (pixels) | 1.3260 | 2.3912 |

Simulation ${\theta}_{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$ |

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.

## 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}$, ${\theta}_{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,^{16} 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

## 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.