Modeling and simulation of multispectral imaging through anisoplanatic atmospheric optical turbulence

Russell C. Hardie; Michael A. Rucci; Santasri R. Bose-Pillai; Richard Van Hook; Barry K. Karch

doi:10.1117/1.OE.61.9.093102

10 September 2022 Modeling and simulation of multispectral imaging through anisoplanatic atmospheric optical turbulence

Russell C. Hardie, Michael A. Rucci, Santasri R. Bose-Pillai, Richard Van Hook, Barry K. Karch

Author Affiliations +

Optical Engineering, Vol. 61, Issue 9, 093102 (September 2022). https://doi.org/10.1117/1.OE.61.9.093102

Abstract

We explore the modeling and simulation of multispectral imaging through anisoplanatic atmospheric optical turbulence. We analyze the impact of wavelength on a number of key atmospheric optical turbulence statistics. This includes the impact of wavelength on tilt and tilt variance. The modeling analysis also includes the impact of wavelength on the atmospheric optical transfer function. Here, we investigate the balance between diffraction and turbulence degradation as a function of wavelength. We also present a method for simulating atmospheric degradation for multispectral imagery using numerical wave propagation. Our approach uses a phase screen resampling method to allow for modeling the same atmospheric realization but with sampling parameters tailored to each wavelength. A number of multispectral simulation results, along with a validation study that compares the empirical statistics from the simulation to their theoretical counterparts, are presented. Real image data are also studied to validate theoretical multispectral tilt statistics.

1. Introduction

Atmospheric optical turbulence is a major source of image degradation in long range imaging.¹^,² Temperature and pressure variations along the optical path lead to changes in the refractive index of air. This leads to amplitude and phase disturbances in the wavefront at the pupil plane of the optics that ultimately produce blurring and warping artifacts in acquired imagery. Wind and convection make this a time varying phenomenon. In very narrow field-of-view (FOV) acquisition scenarios, the degradation in each frame is isoplanatic and can be effectively modeled with a single point spread function (PSF) for each frame. When the FOV is larger, we get the more complex anisoplanatic scenario with spatially varying PSF degradation.

The simulation of realistic atmospheric turbulence degradation has been studied extensively. This is done to better understand the phenomenology and how it impacts image acquisition. It is also done to help develop and evaluate turbulence mitigation (TM) methods. Isoplanatic turbulence simulation has been widely used for astronomical imaging applications.² More recently, anisoplanatic turbulence simulation has been developed to model terrestrial incoherent imaging over long horizontal paths. Numerical wave propagation methods are widely used for both isoplanatic³ and anisoplanatic turbulence simulations.⁴^,⁵ A comparison of numerical simulation techniques is presented in Ref. 6. A fast warping-only anisoplanatic simulator is presented in Ref. 7. This method generates random tilt fields with the proper turbulence tilt statistics by filtering white noise arrays. A fast simulation method based on sampling spatially correlated Zernike coefficients was recently presented in Ref. 8.

Addressing TM in the anisoplanatic case is a particularly challenging problem due to the spatially and temporally varying nature of the degradation. Modeling and simulation of anisoplanatic optical turbulence is an important part of developing and evaluating TM algorithms. For example, anisoplanatic turbulence simulators have been used recently to provide training data for machine learning-based TM.⁹ Such simulators have been used in many other works¹⁰^–¹⁴ for TM algorithm development, tuning, and quantitative performance evaluation. Most of the focus of anisoplantic turbulence modeling, simulation, and mitigation has been for a single wavelength.

In this paper, we explore the modeling and simulation of multispectral imaging through anisoplanatic atmospheric optical turbulence. We analyze the impact of wavelength on a number of key atmospheric optical turbulence statistics. This includes the impact of wavelength on tilt and tilt variance. The modeling analysis also includes the impact of wavelength on the atmospheric optical transfer function (OTF). We investigate the balance between diffraction and turbulence degradation as a function of wavelength. We also present a method for simulating atmospheric degradation for multispectral imagery using numerical wave propagation. Our approach uses a phase screen resampling method to allow for modeling the same atmospheric realization but with sampling parameters tailored to each wavelength. We argue that our approach produces PSFs with appropriate spatial and spectral correlation.

The organization of the remainder of this paper is as follows. In Sec. 2, we focus on multispectral atmospheric turbulence statistics. In particular, we explore the impact of wavelength on a number of well-known statistics for atmospheric characterization. Multispectral anisoplanatic turbulence simulation is addressed in Sec. 3. Experimental results with the proposed simulator and real data are presented in Sec. 4. Finally, we offer conclusions in Sec. 5.

2. Atmospheric Turbulence Characterization

In this section, we examine the impact of wavelength on a number of key atmospheric optical turbulence statistics, including the OTF.

2.1.

Refractive Index Structure Function

The Kolmogorov atmospheric turbulence model gives rise to the refractive index structure function.² This provides a statistical description of the variation in the index of refraction as a function of the distance between points for a locally isotropic atmosphere. Here, we express this function with wavelength dependence as follows:

Eq. (1)

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

where

x

and

r

are three-dimensional (3D) spatial coordinate vectors and

r = | r |

. The wavelength-dependent atmospheric index of refraction is given by

n_{λ} (\cdot)

, and the notation

⟨ \cdot ⟩

represents an ensemble mean operator. The turbulence strength in Eq. (1) is captured by the refractive index structure parameter,

C_{n, λ}^{2}

, in units of

m^{- 2 / 3}

.² To represent variability along the optical path, this parameter is expressed with the function

C_{n, λ}^{2} (z)

, where

z

is the distance along the optical path.

Note that, in the visible to short-wave infrared spectral range, the spectral ratio of refractive index structure parameters can be fairly accurately estimated with the spectral ratio of squared atmospheric refractivities.¹⁵^,¹⁶ This is expressed as

Eq. (2)

C_{n, λ_{2}}^{2} (z) = \frac{N^{2} (λ_{2})}{N^{2} (λ_{1})} C_{n, λ_{1}}^{2} (z),

where the atmospheric refractivity is given by

N (λ) = [n_{λ} - 1] \times 10^{6}

and

n_{λ}

is the index of refraction of the atmosphere for wavelength

λ

at a specified pressure, temperature, and composition. The refractivity for standard dry air using Edlén’s relation¹⁷ is given as

Eq. (3)

N (λ) = (n_{λ} - 1) \times 10^{6} = 83.4213 + \frac{24060.3}{130 - 1 / λ^{2}} + \frac{159.97}{38.9 - 1 / λ^{2}},

for wavelength in microns. This corresponds to a temperature of 288.15 K (15°C), pressure of 101325 Pa (760 Torr), and

{CO}_{2}

abundance of 0.0003 by volume.¹⁷ Edlén provides a multiplicative scaling factor to account for other pressures and temperatures for dry air. This is expressed as

Eq. (4)

N (λ, p, t) = N (λ) N (p, t) = N (λ) \times \frac{0.00138823 p}{1 + 0.003671 t},

where

p

is the pressure in Torr and

t

is the temperature in °C. Note that this pressure and temperature scaling factor,

N (p, t)

, is not a function of wavelength. Thus, although atmospheric turbulence produces disturbances in the pressure and temperature that change the index of refraction, the refractivity ratio for two wavelengths remains constant using Edlén’s dry air model. We use this fact in Appendix A to explicitly derive the relation in Eq. (2).

Note that other models for atmospheric refractivity may be used, such as Ciddor’s model.¹⁸ Although the refractivity values may vary for different models and atmopsheric conditions, the wavelenth ratio values tend to be very similar to that found with the simpler model in Eq. (3). For this reason, we use Eq. (3) to model refractivity ratios for the purposes of this paper.

The atmospheric refractivity as given by Edlén’s relation in Eq. (3) is shown in Fig. 1 for wavelengths ranging from 0.4 to $1.5 μ m$ . The wavelength scaling parameter for the refractive index structure parameter from Eq. (2), $N^{2} (λ_{2}) / N^{2} (λ_{1})$ , is shown in Fig. 2 using Edlén’s relation. Here, the reference wavelength is $λ_{1} = 0.88 μ m$ and $λ_{2}$ is shown on the horizontal axis. Compared with the reference wavelength, the refractive index structure parameter is $\sim 3 %$ higher at 0.5 μm and 1% lower at $1.5 μ m$ . This has implications for all of the parameters that depend on the refractive index structure parameter, as shown in the following sections.

Fig. 1

Atmospheric refractivity using Edlén’s relation in Eq. (3) for standard dry air.

Fig. 2

The $C_{n, λ}^{2}$ scaling from Eq. (2) as a function of wavelength for a reference wavelength of $0.88 μ m$ using Edlén’s relation.

2.2.

Fried Parameter

The wavelength-specific atmospheric coherence diameter (or Fried parameter)¹^,³ is expressed as a weighted integral of $C_{n, λ}^{2} (z)$ yielding

Eq. (5)

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

where

L

is the path length. This expression is for spherical wave propagation with

z = 0

being at the source. Note that a smaller Fried parameter corresponds to higher turbulence strength. This parameter plays a key role in atmospheric OTF models and atmospheric tilt statistics.

Based on Eq. (5) and the approximation in Eq. (2), the Fried parameter at one wavelength is related to that at another as

Eq. (6)

r_{0} (λ_{2}) = (\frac{N (λ_{1}) λ_{2}}{N (λ_{2}) λ_{1}})^{6 / 5} r_{0} (λ_{1}) .

The derivation of Eq. (6) is provided in Appendix B. The

r_{0} (λ)

scaling parameter in Eq. (6) using Edlén’s relation in Eq. (3) is shown in Fig. 3 for a reference wavelength of

λ_{1} = 0.88 μ m

with

λ_{2}

shown on the horizontal axis. This plot clearly illustrates the wavelength-dependent nature of the Fried parameter. This is true even if one assumes a constant atmospheric refractivity spectral ratio, as seen in Eq. (6).

Fig. 3

Fried parameter scaling from Eq. (6) using Edlén’s relation and a reference wavelength of $0.88 μ m$ .

2.3.

Tilt and Tilt Variance

The next statistic that we consider is tilt variance. This statistic provides a method to characterize the warping effects of turbulence. Let an observed point-source direction angle relative to the optical axis be defined as $θ = [θ_{x}, θ_{y}]^{T}$ . The two-axis tilt vector is then defined as $α_{λ} (θ) = [α_{λ, x} (θ), α_{λ, y} (θ)]^{T}$ . The two-axis $Z$ -tilt variance is defined as $T_{Z}^{2} (λ) = ⟨ α_{λ} (θ) \cdot α_{λ} (θ) ⟩$ . Tyson¹⁹ provides an expression for the two-axis $Z$ -tilt variance in terms of the Fried parameter as

Eq. (7)

T_{Z} (λ)^{2} = 0.3641 (\frac{D}{r_{0} (λ)})^{5 / 3} (\frac{λ}{D})^{2} = 0.3641 r_{0} (λ)^{- 5 / 3} D^{- 1 / 3} λ^{2},

where

D

is the aperture diameter. Equivalently, this is expressed in terms of the refractive index structure parameter as

Eq. (8)

T_{Z} (λ)^{2} = 6.0802 D^{- 1 / 3} [\int_{z = 0}^{z = L} C_{n, λ}^{2} (z) (\frac{z}{L})^{5 / 3} d z] .

Thus, the tilt variance is proportional to any constant scaling of the refractive index structure parameter, and the root mean squared (RMS)

Z

-tilt is proportional to the square root of the scaling constant. Also, based on Eqs. (8) and (2), the wavelength dependence of the tilt variance is expressed as

Eq. (9)

T_{Z}^{2} (λ_{2}) = \frac{N^{2} (λ_{2})}{N^{2} (λ_{1})} T_{Z}^{2} (λ_{1}) .

Thus, the tilt variance wavelength scaling is the same as that of the refractive index structure parameter in Eq. (2) and shown in Fig. 2. Note that, unlike the Fried parameter, tilt variance only shows a wavelength dependence resulting from a change in

C_{n, λ}^{2}

. Because the dependence on

C_{n, λ}^{2}

is the same for tilt correlation and differential tilt variance,⁵^,²⁰ their wavelength scaling is also the same as that in Eq. (2) and Fig. 2. The same is true for the two-dimensional (2D) tilt correlation, patch tilt variance, and residual tilt variance derived in Ref. 14.

We may also directly relate a tilt realization at one wavelength to that at another for the same atmosphere as

Eq. (10)

α_{λ_{2}} (θ) = \frac{N (λ_{2})}{N (λ_{1})} α_{λ_{1}} (θ) .

This is a consequence of the refractivity model in Eq. (4) in which the dispersive factor,

N (λ)

, is independent of pressure and temperature. This allows fluctuations in refractivity as a result of pressure and temperature along the optical path to cancel for wavelength ratios. We derive Eq. (10) in Appendix C. Note that Eq. (10) is consistent with Eq. (9) because the tilt mean is zero.

We believe that Eq. (10) is an important relationship for several reasons. One is that it may be used to help validate multispectral simulations. Furthermore, the warping at different wavelengths scales according to a constant refractivity ratio can be used to aid in image registration for multispectral imagery. Finally, one can use Eq. (10) in conjunction with a fast turbulence warping simulator such as that proposed by Schwartzman et al.⁷ In particular, a realization of statistically accurate warping tilts can be generated for one wavelength and then simply scaled according to Eq. (10) for a different wavelength in a multispectral simulation. Note that a single-band fast warping simulator based on Ref. 7 is used in Ref. 9 for training a machine learning TM algorithm. This method employs the tilt statistics reported in Ref. 14 with added blurring from the average short-exposure OTF. This provides a fast approximate simulation with anisoplanatic warping and isoplanatic blurring.

2.4.

Isoplanatic Angle

The isoplanatic angle is another important statistic that gives a measure of how spatially invariant the effects of turbulence are. A large isoplanatic angle indicates a higher degree of spatial invariance, and vice versa. Two point sources separated by less than the isoplanatic angle will have an average wave function phase difference at the aperture of <1 rad.²^,³ The wavelength-dependent version of the isoplanatic angle is expressed as

Eq. (11)

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

Using steps similar to those in Appendix B, we express the isoplanatic angle wavelength scaling as

Eq. (12)

θ_{0} (λ_{2}) = (\frac{N (λ_{1}) λ_{2}}{N (λ_{2}) λ_{1}})^{6 / 5} θ_{0} (λ_{1}) .

Note that the scaling in Eq. (12) is the same as that of the Fried parameter in Eq. (6), which is plotted in Fig. 3.

2.5.

Log Amplitude Variance

The last statistic that we consider in this section is the log-amplitude variance.³ This statistic reflects fluctuations in the amplitude of the wave function in the pupil plane and is given as

Eq. (13)

σ_{χ}^{2} (λ) = 0.563 (\frac{2 π}{λ})^{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 .

Using the expressions in Eqs. (13) and (2) and steps similar to those in Appendix B, we express the wavelength scaling of the log-amplitude variance as

Eq. (14)

σ_{χ}^{2} (λ_{2}) = \frac{N^{2} (λ_{2})}{N^{2} (λ_{1})} (\frac{λ_{1}}{λ_{2}})^{7 / 6} σ_{χ}^{2} (λ_{1}) .

The scaling term in Eq. (14) is plotted in Fig. 4 for a fixed reference wavelength of

λ_{1} = 0.88 μ m

.

Fig. 4

Log amplitude variance scaling from Eq. (14).

2.6.

Optical Transfer Functions

Consider an OTF model that includes the effects of both diffraction and the atmosphere. Such an OTF is expressed as

Eq. (15)

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 model developed by Fried¹^,² is given as

Eq. (16)

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. When the parameter

α

is set to zero, we get the long-exposure OTF. When

α = 1

, we get the average near-field short-exposure OTF. It is interesting to note that

α

controls a Gaussian component of the atmospheric OTF that is associated with tilt variance.¹² This parameter has been used to model the level tilt variance correction achieved when image registration and fusion are employed.¹²^–¹⁴ Note that

α = 0

corresponds to no registration prior to image fusion and

α = 1

corresponds to ideal registration.

Next, we consider the diffraction-limited OTF for a circular aperture. Writing this well-known model²¹ to explicitly show the wavelength dependence yields

Eq. (17)

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

where

F = l / D

is the

f

-number and

1 / (λ F)

is the optical cutoff frequency. Note that the theoretical atmospheric PSF is found by applying the inverse Fourier transform to Eq. (15).

It is interesting to note that the diffraction component in Eq. (17) becomes more low-pass and less favorable at longer wavelengths, whereas the atmospheric OTF in Eq. (16) becomes more favorable. Thus, these OTF components have competing influences on the overall OTF in Eq. (15). However, diffraction generally tends to be the dominant factor. One way to visualize the impact of wavelength on the OTF is from the OTF ratio

Eq. (18)

H (ρ, λ_{1}, λ_{2}) = \frac{H (ρ, λ_{2})}{H (ρ, λ_{1})} .

If the ratio is greater than one for a spatial frequency,

ρ

, then the OTF at

λ_{2}

is more favorable than that at

λ_{1}

for that frequency. Figure 5 shows the OTFs from Eqs. (17) and (15) and the ratio in Eq. (18) for two wavelengths and two turbulence levels. The optical parameters used in these plots are listed in Table 1. These parameters were chosen because they correspond to well-validated simulation results in prior papers.⁵^,²² Note in Figs. 5(a) and 5(b) that the diffraction OTF is more low pass and has a lower cutoff frequency at

λ_{2} = 1.50 μ m

than at

λ_{1} = 0.88 μ m

. For the lower turbulence level shown in Fig. 5(a), the overall OTF and OTF ratio clearly favor the shorter wavelength. However, at the higher turbulence level shown in Fig. 5(b), the overall OTF and ratio favor the longer wavelength slightly until near the optical cutoff frequency.

Fig. 5

OTFs from Eqs. (17), (15), and (18) for optical parameters in Table 1 for two turbulence levels: (a) $C_{n, 0.88 μ m}^{2} = 1.0 \times 10^{- 15} m^{- 2 / 3}$ and (b) $C_{n, 0.88 μ m}^{2} = 2.5 \times 10^{- 15} m^{- 2 / 3}$ .

Table 1

Optical parameters used for the example and simulation results.

Parameter	Value
Aperture	$D = 0.2034 m$
Focal length	$l = 1.2 m$
$F$ -number	$F = 5.9$
$C_{n, λ}^{2}$ reference wavelength	$λ = 0.88 μ m$
Object distance	$L = 7 km$
Nyquist pixel spacing (focal plane)	$δ_{f} = 1.4749 μ m$ @ $λ = 0.5 μ m$
Nyquist pixel spacing (object plane)	$δ_{o} = 8.6037 mm$ @ $λ = 0.5 μ m$

Another potential use for the OTF ratio in Eq. (18) is to modify a degradation at one wavelength to that corresponding to another wavelength. This could be used to turn a single-band simulation into a simple multispectral simulation. Applying the system in Eq. (18) to an image degraded at $λ_{1}$ would make it more consistent with a degradation for $λ_{2}$ . The average of a sequence of such frames would have the correct average OTF. However, the individual frames would not necessarily have the correct anisoplanatic blur and the tilts would not be adjusted according to Eq. (10) to account for any wavelength change in $C_{n, λ_{2}}^{2}$ . Notwithstanding these factors, such a simple approach may be useful in some applications.

An interesting way to visualize the trade-off between diffraction and atmospheric degradation is with the Strehl ratio.²³ The Strehl ratio is defined as the peak PSF value of an aberrated system divided by the peak PSF value of the corresponding diffraction-limited system. Note that a larger Strehl ratio is desirable and a value of one corresponds to a diffraction-limited system. Examples of Strehl ratios for the system in Table 1 are shown in Fig. 6. In these plots, we show the Strehl ratios for the average short-exposure PSFs based on Eq. (15) for a range of turbulence strengths and wavelengths. In Fig. 6(a), the Strehl ratios are computed with respect to the diffraction-limited PSF at each wavelength. In Fig. 6(b), we show a modified Strehl ratio in which all PSF peak values are divided by the peak of the diffraction-limited PSF at the shortest wavelength. The mesh plot in Fig. 6(a) reveals how the impact of turbulence is less severe at longer wavelengths, relative to diffraction. In Fig. 6(b), we see how increasing turbulence and wavelength each reduce the PSF peak. It is interesting to note that, as the turbulence level increases, the peak of the modified Strehl ratio for that turbulence level occurs at slightly longer wavelengths. This is consistent with the OTF depicted in Fig. 5(b).

Fig. 6

Strehl ratios as a function of wavelength and turbulence strength for the optical parameters in Table 1. The PSF peak is divided by the peak of the diffraction-limited PSF for (a) each wavelength and (b) the shortest wavelength.

3. Multispectral Numerical Wave Propagation

In this section, we describe our method for performing multispectral turbulence simulation using numerical wave propagation. We begin with an overview and then focus on one of the key aspects of the proposed multispectral numerical wave propagation method, phase screen resampling.

3.1.

Overview

Our multispectral method is based on the single-band simulator presented in Ref. 5, but it has some important modifications. As with the single-band simulator, we use the phase screen geometry shown in Fig. 7. Note that this is similar to that first introduced in Ref. 4. A grid of simulated point sources are numerically propagated from the object plane to the pupil plane through a series of phase screens. All of the details of the point source model, phase screen model, and slit-step propagation are provided in Ref. 5. The result of the propagation is an array of PSFs, with one PSF for each pixel location in the simulated image. The PSFs are applied in a spatially varying 2D convolution operation to a pristine truth image to generate the simulated atmospheric blur and warping. Because neighboring point sources travel though common overlapping phase screen regions, the approach generates PSFs with appropriate spatial correlation. The phase screen overlap for two selected point source locations is shown in Fig. 7 with the green and blue boxes. This approach was shown in Ref. 5 to produce simulations with key statistics that closely match those predicted by the Kolmogorov theory. These include the long exposure PSF, average short exposure PSF, tilt variance, tilt correlation, differential tilt variance, and isoplanatic angle.

Fig. 7

Phase screen geometry for the anisoplanatic numerical wave propagation simulation method.⁵ Note that two point source propagation paths are shown from the object plane to the pupil plane; one is blue and one is green. Spatial correlation is imparted due to overlapping portions of the phase screens.

Our goal here is to extend the numerical wave propagation method for multispectral simulation. In particular, we want to generate correlated PSF arrays for different wavelengths corresponding to a common atmosphere realization. Note that the statistics of the phase screens vary with wavelength, as do the point source model and propagation filter.⁵ Furthermore, note that the point source array in the object plane has a spacing governed by the diffraction-limited Nyquist rate for the imaging sensor. We compute the spacing for the shortest wavelength being simulated and use that same spacing for all wavelengths. We do this to model a multispectral sensor with registered spectral channels. Note that we assume that the spectral bands are acquired simultaneously so as to share the same atmospheric realization.

3.2.

Phase Screen Resampling

The phase screens are created by first generating white noise random fields and then filtering these to produce phase screens with the desired statistics.⁵ To model a common atmosphere, we employ the same white noise realizations at all wavelengths. Note that the filtering of those common input random fields will still vary to produce the desired phase screens at each wavelength according to the wavelength-specific phase screen statistics.

One phase screen parameter that is critical for successful numerical wave propagation is the sample spacing. The selection of this parameter is addressed in detail by Rucci et al.²² In that work, a sampling rule is proposed that sets the phase screen sample spacing as follows:

Eq. (19)

Δ_{x} (λ) = \frac{λ L}{a_{g} D + c \frac{λ L}{r_{0} (λ)}},

where

a_{g}

is referred to as the aperture gain and

c

is a turbulence-sensitivity parameter that controls the impact of

r_{0}

on the sample spacing. The aperture gain is a scalar multiplier that determines the bandwidth of the point source model employed.³^,²² This bandwidth is linked to the camera aperture so that is appropriate for the diffraction-limited optics. The parameter

c

controls how much the sample spacing is impacted by the level of turbulence.³^,²² We found that we achieve a high level of agreement with the theoretical multispectral statistics studied here using

a_{g} = 5

and

c = 2

.

The number of samples cropped from the extended phase screen and used for propagation is another parameter that requires careful attention for successful propagation. The rule proposed by Schmidt³ and also used in Ref. 22 is

Eq. (20)

N_{\min} (λ) = \frac{(a_{g} + 1) D + 2 c \frac{λ L}{r_{0} (λ)}}{2 Δ_{x} (λ)} + \frac{λ L}{2 Δ_{x} (λ)^{2}} .

The number of samples,

N (λ)

, must be such that

N (λ) \geq N_{\min} (λ)

. The sampling relations in Eqs. (19) and (20) are plotted in Fig. 8 as a function of wavelength and the refractive index structure parameter. In Fig. 8(a), we see that longer wavelengths call for a larger sample spacing. On the other hand, increased turbulence calls for a smaller sample spacing for a given wavelength. Figure 8(b) shows that increased turbulence and decreased wavelength call for more samples.

Fig. 8

Phase screen sampling parameters²² as a function of wavelength and the refractive index structure parameter. The sample spacing from Eq. (19) is shown in (a) and the number of samples from Eq. (20) is shown in (b).

A challenge for multispectral simulation is that the propagation parameters in Eqs. (19) and (20) are highly wavelength dependent. Thus, we have competing requirements. To model a common atmosphere realization, we need to use the same underlying random number array for the phase screens at each wavelength and have the phase screens span the same physical size. However, for effective propagation, we need different sample spacing at each wavelength. Our solution is to initially generate all phase screens for all wavelengths using the sample spacing for the minimum wavelength, $Δ_{x} (λ_{\min})$ . This provides the appropriate sampling for the minimum wavelength phase screens and oversampling for the longer wavelength phase screens. Using a common sample spacing ensures consistency in terms of the physical dimensions of all screens across wavelength from common random white noise fields. Then, after the screens are generated, we resample the screens for the longer wavelength simulations according to Eq. (19). An antialiasing filter is used prior to resampling because we are reducing the sampling rate of the phase screens for longer wavelengths. When extracting propagation windows from the full phase screens, as shown in Fig. 7, we take care to ensure that the windows are large enough to satisfy Eq. (20) for each wavelength. We find that the phase screen resampling is critical for all but very small steps in wavelength. When the resampling is not applied, we find that artifacts tend to emerge in the PSFs and they no longer have the proper tilt characteristics to conform to Eqs. (9) and (10).

4. Experimental Results

In this section, we present a number of simulation results and a validation analysis to demonstrate the efficacy of the proposed multispectral simulation method. The key simulation parameters used here are listed in Table 2. These parameters closely follow those in the original simulation paper.⁵ This is done because the single-band simulation has been well validated with these parameters. The main difference here lies in the wavelengths used and how the phase screen generation and resampling are performed based on $Δ_{x} (λ)$ in Eq. (19). In Sec. 4.1, we consider the generation of independent PSFs to validate the multispectral propagation method with phase screen resampling. In Sec. 4.2, we consider the full multispectral anisoplanatic image degradation process. Section 4.3 presents an image tilt analysis using real imagery from a multispectral camera.

Table 2

Multispectral simulation parameters.

Parameter	Value
Wavelengths	$λ = 0.5$ , 1.0, and $1.5 μ m$
Turbulence levels	$C_{n, 0.88 μ m}^{2} = 1.00$ , 5.00, and $10.00 \times 10^{- 16} m^{- 2 / 3}$
Path length	$L = 7 km$
Number of phase screens	$N = 10$ (9 non-zero)
Propagation step	$Δ z = 700 m$
Aperture gain	$a_{g} = 5$
Turbulence sensitivity	$c = 2$
Phase screen type	Modified von Kármán with subharmonics
Inner scale	$l_{0} = 0.01 m$
Outer scale	$L_{0} = 300 m$

4.1.

Multispectral PSF Analysis

Figure 9 shows simulated PSFs for three wavelengths and three turbulence strengths for a common random array realization and using phase screen resampling. The rows from top to bottom represent wavelengths of $λ = 0.5$ , 1.0, and $1.5 μ m$ . The columns from left to right represent constant refractive index structure parameters of $C_{n, 0.88 μ m}^{2} = 1 \times 10^{- 16} m^{- 2 / 3}$ , $5 \times 10^{- 16} m^{- 2 / 3}$ , and $10 \times 10^{- 16} m^{- 2 / 3}$ . Note that, as the wavelength increases, moving down in the figure, the PSFs are more dominated by diffraction and appear more circularly symmetric. As the turbulence strength increases, moving left to right in the figure, we see that the PSFs broaden due to increased turbulence. Because we use the same white noise realizations for each scenario, the PSFs are all correlated. Notice that the centroid moves farther from the origin as the turbulence level increases. This is consistent with Eq. (8). However, the centroids are nearly constant as the wavelength changes, as prescribed by Eq. (10). Also, note that the shorter wavelength is more impacted by increasing turbulence than the longer wavelength. This is consistent with the theoretical Strehl analysis in Fig. 6.

Fig. 9

Simulated PSF for three wavelengths and three turbulence strengths for a common random array realization. The rows represent wavelengths of $λ = 0.5$ , 1.0, and $1.5 μ m$ , respectively. The columns represent $C_{n, 0.88 μ m}^{2} = 1 \times 10^{- 16} m^{- 2 / 3}$ , $5 \times 10^{- 16} m^{- 2 / 3}$ , and $10 \times 10^{- 16} m^{- 2 / 3}$ , respectively.

We generated 300 realizations of PSFs for each of the scenarios depicted in Fig. 9. This comprises wavelengths of 0.5, 1.0, and $1.5 μ m$ for $C_{n, 0.88 μ m}^{2}$ values of $1 \times 10^{- 16} m^{- 2 / 3}$ , $5 \times 10^{- 16} m^{- 2 / 3}$ , and $10 \times 10^{- 16} m^{- 2 / 3}$ . Quantitative validation results for these PSFs are summarized in Tables 3 Table 4–5. In these tables, units of pixels refer to diffraction-limited Nyquist pixel spacings at the minimum wavelength as reported in Table 1. We believe these units are more intuitive when it comes to image simulation and analysis. The RMS $Z$ -tilt reported in the tables is the one-axis $Z$ -tilt. The theoretical one-axis RMS $Z$ -tilt is computed in pixels by first computing Eq. (8) to give units of radians squared. Assuming isotropicity, we divide by two to get the one-axis $Z$ -tilt variance and then take the square root to get units of radians. Finally, we use the small angle approximation and multiply these radian angles by the focal length and then divide by the pixel spacing. The simulation Fried parameter is estimated by fitting the theoretical long exposure PSF to the average simulated PSF.

Table 3

Multispectral simulation PSF quantitative results for Cn,0.88 μm2=1×10−16 m−2/3.

Parameter	Value
Wavelength $λ$ ( $μ m$ )	0.5000	1.0000	1.5000
$C_{n, λ}^{2}$ ( $\times 10^{- 16} m^{- 2 / 3}$ )	1.0319	0.9967	0.9904
$Δ_{x} (λ)$ (m) (after resample)	0.0033	0.0067	0.0100
$N (λ)$ (samples) (after resample)	456	226	150
$N_{\min} (λ)$ (samples)	356	176	117
Isoplanatic angle $θ_{0}$ (pixels)	6.4307	15.0851	24.6318
$D / r_{0}$	1.1558	0.4927	0.3017
Theoretical $r_{0}$ (m)	0.1760	0.4128	0.6741
Simulation $r_{0}$ (m)	0.1766	0.4012	0.6754
Theoretical RMS $Z$ -tilt (pixels)	0.9628	0.9462	0.9432
Simulation RMS $Z$ -tilt (pixels)	0.9690	0.9520	0.9450

Table 4

Multispectral simulation PSF quantitative results for Cn,0.88 μm2=5×10−16 m−2/3.

Parameter	Value
Wavelength $λ$ ( $μ m$ )	0.5000	1.0000	1.5000
$C_{n, λ}^{2}$ ( $\times 10^{- 16} m^{- 2 / 3}$ )	5.1596	4.9835	4.9522
$Δ_{x} (λ)$ (m) (after resample)	0.0031	0.0063	0.0096
$N (λ)$ (samples) (after resample)	456	224	148
$N_{\min} (λ)$ (samples)	409	198	130
Isoplanatic angle $θ_{0}$ (pixels)	2.4484	5.7434	9.3781
$D / r_{0}$	3.0357	1.2941	0.7925
Theoretical $r_{0}$ (m)	0.0670	0.1572	0.2566
Simulation $r_{0}$ (m)	0.0683	0.1566	0.2555
Theoretical RMS $Z$ -tilt (pixels)	2.1528	2.1157	2.1091
Simulation RMS $Z$ -tilt (pixels)	2.1657	2.1303	2.1157

Table 5

Multispectral simulation PSF quantitative results for Cn,0.88 μm2=10×10−16 m−2/3.

Parameter	Value
Wavelength $λ$ ( $μ m$ )	0.5000	1.0000	1.5000
$C_{n, λ}^{2}$ ( $\times 10^{- 16} m^{- 2 / 3}$ )	10.3192	9.9669	9.9044
$Δ_{x} (λ)$ (m) (after resample)	0.0030	0.0061	0.0092
$N (λ)$ (samples) (after resample)	456	223	147
$N_{\min} (λ)$ (samples)	456	218	142
Isoplanatic angle $θ_{0}$ (pixels)	1.6153	3.7892	6.1872
$D / r_{0}$	4.6012	1.9615	1.2013
Theoretical $r_{0}$ (m)	0.0442	0.1037	0.1693
Simulation $r_{0}$ (m)	0.0458	0.1038	0.1679
Theoretical RMS $Z$ -tilt (pixels)	3.0445	2.9921	2.9827
Simulation RMS $Z$ -tilt (pixels)	2.9447	3.0182	3.0005

Note in Tables 3 Table 4–5 that the wavelength-dependent $C_{n, λ}^{2}$ decreases with the increasing wavelength, as does the RMS $Z$ -tilt. These relations follow from Eqs. (2) and (9), respectively. Also, the isoplanatic angle and Fried parameter increase with the increasing wavelength, as expressed by Eqs. (12) and (6), respectively. One can also see in Tables 3 Table 4–5 that there is generally good agreement between the simulated and theoretical one-axis RMS $Z$ -tilts. The same is true for the Fried parameter.

The simulated long-exposure and average short-exposure PSF cross-sections are shown in Fig. 10 compared with the theoretical curves for the $λ = 1.0 μ m$ case with $C_{n, 0.88 μ m}^{2} = 5 \times 10^{- 16} m^{- 2 / 3}$ . Note that the PSFs are plotted versus pixel spacings. Again, we use the diffraction-limited Nyquist pixel spacing for $λ = 0.5 μ m$ for all wavelengths (see Table 1). The results show excellent agreement between simulation and theory. Similar results are observed with the other two wavelengths tested.

Fig. 10

Comparison of theoretical and simulated PSF cross-sections for $C_{n, 0.88 μ m}^{2} = 5 \times 10^{- 16} m^{- 2 / 3}$ and $λ = 1.0 μ m$ using the optical parameters in Table 1 and simulation parameters in Table 2. The long exposure PSF is shown in (a) and the average short exposure PSF is shown in (b).

Figure 11 shows the individual $x$ and $y$ PSF tilts for the first 100 realizations in units of pixels. The tilts in Fig. 11(a) correspond to PSFs generated using the proposed phase screen resampling method, that is, the phase screens are initially generated for all wavelengths using $Δ_{x} (λ_{\min})$ using the same random arrays between wavelengths. Then, for all but the minimum wavelength, the phase screens are resampled using $Δ_{x} (λ)$ . Note that the number of samples also changes such that the physical dimensions of the resampled phase screens remain constant and the geometry illustrated in Fig. 7 is not altered between wavelengths. In Fig. 11(b), the PSFs are generated using a constant phase screen sample spacing of $Δ_{x} (λ_{\min})$ with no resampling. Note that the tilts are in near agreement for all three wavelengths in Fig. 11(a) using resampling. The longer wavelengths do exhibit a reduction scaling as prescribed by Eq. (10). On the other hand, Fig. 11(b) shows tilts that appear spectrally uncorrelated. We attribute this to error in the numerical wave propagation as a result of improper sampling parameters for the phase screens. The results in Fig. 11 clearly indicate the importance of the proposed phase screen resampling process. Note that generating the phase screens directly from the same random fields with the desired sample spacing is also not a desirable option. This is because the spatial distribution of the random array in physical distance is not consistent between wavelengths. This does not correspond to having the same atmosphere for all wavelengths.

Fig. 11

Simulated PSF tilts for the first 100 realizations for wavelengths of $λ = 0.5$ , 1.0, and $1.5 μ m$ , and $C_{n, 0.88 μ m}^{2} = 5 \times 10^{- 16} m^{- 2 / 3}$ . Tilts for PSFs using the proposed phase screen resampling are shown in (a) and tilts for PSFs using constant phase screen sampling parameters are shown in (b). Note that the tilts in (a) are more consistent with Eq. (10).

4.2.

Multispectral Anisoplanatic Image Simulation

In this section, we present the anisoplanatic image simulation results using the same parameters listed in Tables 1 and 2. The input image comes from a publicly available dataset acquired with the NASA/JPL airborne visible/infrared imaging spectrometer (AVIRIS) sensor.²⁴ The orthorectified data are from September 30, 2011, over Fisherman’s Wharf in Monterey, California. Simulated PSFs are generated for wavelengths $λ = 0.5 μ m$ , $λ = 1.0 μ m$ , and $λ = 1.5 μ m$ . The AVIRIS image bands nearest in wavelength to these are used. False color composite images are generated by displaying the three bands with blue, green, and red, respectively.

The original image and single degraded frames are shown in Fig. 12 for $C_{n, 0.88 μ m}^{2}$ values of $1 \times 10^{- 16} m^{- 2 / 3}$ , $5 \times 10^{- 16} m^{- 2 / 3}$ , and $10 \times 10^{- 16} m^{- 2 / 3}$ . The increased blurring may be observed for the higher values of $C_{n, 0.88 μ m}^{2}$ . However, even though there is an order of magnitude difference in the refractive index structure parameter between Figs. 12(b) and 12(d), the increased blurring appears subtle. The reason for this is largely because the longer wavelengths are less affected by this increase in turbulence, as shown in Fig. 9. It is mainly the blue channel degraded with turbulence for $λ = 0.5 μ m$ that is significantly impacted.

Fig. 12

Anisoplanatic multispectral simulation using AVIRIS imagery over Fisherman’s Wharf in Monterey, California, as truth. The false color images are created by displaying the image at $λ = 0.5 μ m$ as blue, $1.0 μ m$ as green, and $1.5 μ m$ as red. (a) Truth image, (b) $C_{n, 0.88 μ m}^{2} = 1 \times 10^{- 16} m^{- 2 / 3}$ , (c) $C_{n, 0.88 μ m}^{2} = 5 \times 10^{- 16} m^{- 2 / 3}$ , and (d) $C_{n, 0.88 μ m}^{2} = 10 \times 10^{- 16} m^{- 2 / 3}$ .

The turbulence-induced warping within the center portions of the images from Fig. 12 is shown in Fig. 13 using a quiver plot. Here, the spatially varying tilt for each band is shown with color-coded arrows, with the arrow color matching the color of the image channel. One can clearly see the increased warping at the higher turbulence levels. It is also clear that the simulated tilts are highly correlated between wavelengths as predicted by Eq. (10) and are consistent with the isoplanatic PSF results in Fig. 11(a).

Fig. 13

Anisoplanatic multispectral simulation image ROIs with quiver plots showing the simulated warping. (a) Truth image, (b) $C_{n, 0.88 μ m}^{2} = 1 \times 10^{- 16} m^{- 2 / 3}$ , (c) $C_{n, 0.88 μ m}^{2} = 5 \times 10^{- 16} m^{- 2 / 3}$ , and (d) $C_{n, 0.88 μ m}^{2} = 10 \times 10^{- 16} m^{- 2 / 3}$ . The arrows are scaled in length by a factor of 2 to improve their visibility.

To explore the anisoplanatic PSF statistics, 300 simulated multispectral frames are generated at each of the three turbulence levels in Table 2. The resulting PSF statistics are shown in Tables 6 Table 7–8. Note that, as before, we see good agreement between the theoretical and simulated Fried parameters and RMS tilt values. This illustrates that the phase screen geometry depicted in Fig. 7 used for anisoplanatic simulations does not negatively impact the overall multispectral PSF statistics.

Table 6

Multispectral anisoplanatic simulation quantitative results for Cn,0.88 μm2=1×10−16 m−2/3.

Parameter	Value
Wavelength $λ$ ( $μ m$ )	0.5000	1.0000	1.5000
$C_{n, λ}^{2}$ ( $\times 10^{- 16} m^{- 2 / 3}$ )	1.0319	0.9967	0.9904
Theoretical $r_{0}$ (m)	0.1760	0.4128	0.6741
Simulation $r_{0}$ (m)	0.1736	0.3955	0.6649
Theoretical RMS $Z$ -tilt (pixels)	0.9628	0.9462	0.9432
Simulation RMS $Z$ -tilt (pixels)	0.9696	0.9461	0.9284

Table 7

Multispectral anisoplanatic simulation quantitative results for Cn,0.88 μm2=5×10−16 m−2/3.

Parameter	Value
Wavelength $λ$ ( $μ m$ )	0.5000	1.0000	1.5000
$C_{n, λ}^{2}$ ( $\times 10^{- 16} m^{- 2 / 3}$ )	5.1596	4.9835	4.9522
Theoretical $r_{0}$ (m)	0.0670	0.1572	0.2566
Simulation $r_{0}$ (m)	0.0676	0.1569	0.2565
Theoretical RMS $Z$ -tilt (pixels)	2.1528	2.1157	2.1091
Simulation RMS $Z$ -tilt (pixels)	2.1311	2.0744	2.0476

Table 8

Multispectral anisoplanatic simulation quantitative results for Cn,0.88 μm2=10×10−16 m−2/3.

Parameter	Value
Wavelength $λ$ ( $μ m$ )	0.5000	1.0000	1.5000
$C_{n, λ}^{2}$ ( $\times 10^{- 16} m^{- 2 / 3}$ )	10.3192	9.9669	9.9044
Theoretical $r_{0}$ (m)	0.0442	0.1037	0.1693
Simulation $r_{0}$ (m)	0.0437	0.1000	0.1616
Theoretical RMS $Z$ -tilt (pixels)	3.0445	2.9921	2.9827
Simulation RMS $Z$ -tilt (pixels)	3.1727	3.0860	3.0496

The simulation approach depicted in Fig. 7 gives rise to the desired spatial correlation between the PSFs. Here, we demonstrate that the desired correlation extends to the multispectral case. In particular, we compare the theoretical and simulation total tilt correlation⁵^,¹⁴ in Fig. 14 for the three simulated wavelengths. Note that the total tilt correlation curve has a similar shape for all wavelengths, but it has a wavelength scaling that is the same as that provided in Eq. (9).

Fig. 14

Comparison of theoretical and simulation total tilt correlation as a function of pixel separation for the anisoplanatic multispectral image simulation with $C_{n, 0.88 μ m}^{2} = 10 \times 10^{- 16} m^{- 2 / 3}$ .

Another interesting statistic to explore for the multispectral anisoplanatic case is the Strehl ratio. Because the simulator produces a PSF for each pixel, we can readily compute the Strehl ratio across the spatial extent of each frame. Such Strehl ratio maps are provided in Fig. 15 for the first simulated frame at each wavelength with $C_{n, 0.88 μ m}^{2} = 5 \times 10^{- 16} m^{- 2 / 3}$ . The left column shows the Strehl ratios relative to diffraction-limited for the corresponding wavelength, as shown in Fig. 6(a). The right column in Fig. 15 is relative to $λ = 0.5 μ m$ , as shown in Fig. 6(b). Notice the strong spectral correlation exhibited in the Strehl ratio maps. Also note that the impact of diffraction becomes more dominant at longer wavelengths, diminishing the Strehl ratios in the right column. To better illustrate the spectral correlation, a scatter plot of the left column data is shown in Fig. 16. Here, each point shows the Strehl ratio for the three wavelengths at one pixel location.

Fig. 15

Anisoplanatic Strehl ratio images from frame 1 of the multispectral image simulation for $C_{n, 0.88 μ m}^{2} = 5 \times 10^{- 16} m^{- 2 / 3}$ . Strehl ratios relative to diffraction-limited PSF for the corresponding wavelength are shown in the left column. The right column shows Strehl ratios relative to the shortest wavelength.

Fig. 16

Scatter plot of anisoplanatic Strehl ratios from frame 1 of the multispectral image simulation for $C_{n, 0.88 μ m}^{2} = 5 \times 10^{- 16} m^{- 2 / 3}$ . The Strehl ratios are relative to the diffraction-limited PSF for each corresponding wavelength.

4.3.

Real Multispectral Image Tilt Analysis

The final experiment uses real multispectral imagery to study the impact of wavelength on image tilt to validate the relevant theoretical analysis in Sec. 2. The details of the imaging sensor and various statistics are provided in Table 9. Frame 1 from the camera is shown in Fig. 17. The image shifts were estimated using normalized cross-correlation followed by a subpixel gradient-based image registration method.²⁵ The image shifts are shown for each wavelength over 600 frames in Fig. 18. Note that the tilts align very closely, as prescribed by Eq. (10). Based on a scintillometer measurement, we are able to obtain values for $r_{0} (λ)$ and then compute theoretical tilt statistics. The theoretical RMS tilt over the full image is computed using the patch tilt variance approach derived in Ref. 14 assuming a constant $C_{n, λ}^{2} (z)$ profile and using an image-sized patch. These theoretical values are compared with the empirical values in Table 9, from which we see a good level of agreement.

Table 9

Optical and atmospheric parameters for the real multispectral image data.

Parameter	Value
Wavelength $λ$ ( $μ m$ )	0.5500	0.6500	0.8000
Aperture (m)	0.1333	0.1333	0.1333
Focal length (m)	3.7320	3.7320	3.7320
$F$ -number	28.0000	28.0000	28.0000
Nyquist pitch ( $μ m$ )	7.7000	9.1000	11.2000
FPA pitch ( $μ m$ )	6.4500	6.4500	6.4500
Distance (km)	1.8701	1.8701	1.8701
$N (λ) / N (0.55 μ m)$	1.0000	0.9946	0.9900
Measured $r_{0} (λ)$ (m)	0.0310	0.0381	0.0492
$D / r_{0}$	4.2995	3.4956	2.7095
Path average $C_{n, λ}^{2}$ ( $\times 10^{- 15} m^{- 2 / 3}$ )	8.4435	8.3518	8.2750
Isoplanatic angle $θ_{0}$ (pixels)	3.0154	3.7089	4.7849
Theoretical RMS $Z$ -tilt (pixels)	3.4349	3.4162	3.4004
Theoretical RMS image tilt (pixels)	2.4455	2.4322	2.4209
Empirical RMS image tilt (pixels)	2.5979	2.5965	2.6384

Fig. 17

Real multispectral image. (a) $λ = 0.55 μ m$ , (b) $λ = 0.65 μ m$ , (c) $λ = 0.80 μ m$ , and (d) false color composite with blue, green, and red made from (a), (b), and (c), respectively.

Fig. 18

Estimated full-frame image tilts for the real multispectral image sequence for each wavelength channel. Notice that the estimated tilts track closely as predicted by Eq. (10).

5. Conclusions

This paper has examined the modeling and simulation of multispectral atmospheric optical turbulence in the visible to short-wave infrared regime. In particular, Sec. 2 shows how key atmospheric turbulence statistics vary as a function of wavelength. We have explored the consequences of the model in Eq. (4) in which pressure and temperature changes in the atmosphere do not impact the dispersive component of atmospheric refractivity. One such consequence is that the tilt at one wavelength may be expressed in terms of the tilt at another scaled by the refractivity ratio, as given by Eq. (10). Also, tilt variance, tilt correlation, differential tilt variance, and patch tilt variance all simply scale according to the squared refractivity ratio in Eq. (9).

The OTF and Strehl ratio analyses in Sec. 2.6 show that there is a trade-off between increased diffraction and decreased turbulence blurring as the wavelength increases. Scene phenomenology aside, shorter wavelengths are generally favorable because of decreased diffraction. However, with increased turbulence levels, the peak Strehl ratio occurs at slightly longer wavelengths as shown in Fig. 6.

Finally, we have demonstrated a numerical wave propagation simulation methodology capable of producing spatially and spectrally correlated PSFs. These may be used to simulate multispectral image acquisition in which the spectral channels are acquired simultaneously through a common atmospheric realization. The key to the proposed method is the phase screen resampling described in Sec. 3. This allows for using sampling parameters that are tuned to each wavelength, while at the same time producing correlated phase screens spanning the same physical dimensions to model a common atmosphere. We believe the modeling and simulation tools developed here set the stage for the development and evaluation of new multispectral TM algorithms.

6. Appendix A: Spectral Ratio of Refractive Index Structure Parameters

Consider the refractive index structure function from Eq. (1) written equivalently in terms of refractivities as follows:

Eq. (21)

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

Next, the spectral ratio by which we cancel common terms is given as

Eq. (22)

\frac{C_{n, λ_{2}}^{2}}{C_{n, λ_{1}}^{2}} = \frac{⟨ [(n_{λ_{2}} (x + r) - 1) - (n_{λ_{2}} (x) - 1)]^{2} ⟩}{⟨ [(n_{λ_{1}} (x + r) - 1) - (n_{λ_{1}} (x) - 1)]^{2} ⟩} .

If we assume only pressure and temperature changes are present along the relevant optical paths, we may express the refractivities using Eq. (4), yielding

Eq. (23)

\frac{C_{n, λ_{2}}^{2}}{C_{n, λ_{1}}^{2}} = \frac{⟨ [N (λ_{2}) N (p (x + r), t (x + r)) - N (λ_{2}) N (p (x), t (x))]^{2} ⟩}{⟨ [N (λ_{1}) N (p (x + r), t (x + r)) - N (λ_{1}) N (p (x), t (x))]^{2} ⟩} .

Factoring out the wavelength-dependent terms of the refractivities yields

Eq. (24)

\frac{C_{n, λ_{2}}^{2}}{C_{n, λ_{1}}^{2}} = \frac{⟨ [N (λ_{2}) [N (p (x + r), t (x + r)) - N (p (x), t (x))]]^{2} ⟩}{⟨ [N (λ_{1}) [N (p (x + r), t (x + r)) - N (p (x), t (x))]]^{2} ⟩} .

Finally, canceling the common pressure and temperature terms, we obtain the desired result

Eq. (25)

\frac{C_{n, λ_{2}}^{2}}{C_{n, λ_{1}}^{2}} = \frac{N^{2} (λ_{2})}{N^{2} (λ_{1})} .

7. Appendix B: Fried Parameter Scaling Derivation

The Fried parameter as a function of wavelength $λ_{2}$ is given as

Eq. (26)

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

Using Eq. (2), we express this in terms of the refractive index structure parameter for wavelength

λ_{1}

defined as

C_{n, λ_{1}}^{2} (z)

. This is given as

Eq. (27)

r_{0} (λ_{2}) = [0.423 (\frac{2 π}{λ_{2}})^{2} \frac{N^{2} (λ_{2})}{N^{2} (λ_{1})} \int_{z = 0}^{z = L} C_{n, λ_{1}}^{2} (z) (\frac{z}{L})^{5 / 3} d z]^{- 3 / 5} .

Now we separate the terms that define

r_{0} (λ_{1})

from scaling constants to produce

Eq. (28)

r_{0} (λ_{2}) = (\frac{N^{2} (λ_{2})}{N^{2} (λ_{1})} (\frac{λ_{1}}{λ_{2}})^{2})^{- 3 / 5} [0.423 (\frac{2 π}{λ_{1}})^{2} \int_{z = 0}^{z = L} C_{n, λ_{1}}^{2} (z) (\frac{z}{L})^{5 / 3} d z]^{- 3 / 5} .

Noting that the second term in Eq. (28) equals

r_{0} (λ_{1})

, we write

Eq. (29)

r_{0} (λ_{2}) = (\frac{N^{2} (λ_{2})}{N^{2} (λ_{1})} (\frac{λ_{1}}{λ_{2}})^{2})^{- 3 / 5} r_{0} (λ_{1}) .

Finally, combining the exponents and removing the negative by inverting the argument, we obtain

Eq. (30)

r_{0} (λ_{2}) = (\frac{N (λ_{1}) λ_{2}}{N (λ_{2}) λ_{1}})^{6 / 5} r_{0} (λ_{1}) .

8. Appendix C: Tilt Scaling Derivation

The tilt vector for wavelength $λ$ and source originating from angle $θ$ is given by Ref. 26 as

Eq. (31)

α_{λ} (θ) = \frac{32 λ}{π^{2} D^{4}} \int_{r} r W (r) ϕ_{λ} (r, θ) d r,

where

r = [r_{x}, r_{y}]^{T}

are the spatial coordinates over the aperture with the origin in the center and

W (r)

is the binary pupil function. The phase disturbance at the aperture due to atmospheric turbulence is given by the term

ϕ_{λ} (r, θ)

. The optical path difference (OPD) through the atmosphere relative to vacuum propagation is given by the line integral²⁷

Eq. (32)

{OPD}_{λ} (r, θ) = \int_{C (r, θ)} (n_{λ} (x) - 1) d s,

where

n_{λ} (x)

is the index of refraction of the atmosphere in 3D spatial coordinates and

C (r, θ)

defines the optical path for a source originating from angle

θ

to a point on the pupil plane defined by

r

. The term

d s

is interpreted as an elementary arc length.²⁸ The line integral in Eq. (32) is expressed as an integral over distance along the principal axis,

z

, provided we include the magnitude of the path derivative with respect to

z

.²⁸ This is expressed as

Eq. (33)

{OPD}_{λ} (r, θ) = \int_{z = 0}^{L} (n_{λ} (x (r, θ, z)) - 1) | x^{'} (r, θ, z) | d z,

where

x (r, θ, z)

is the 3D optical path as a function of

z

and

x^{'} (r, θ, z) = \partial x (r, θ, z) / \partial z

.

The optical paths for both the spherical wave and plane wave cases are depicted in Fig. 19. The geometry for the spherical wave propagation case with the source at $z = 0$ is shown in Fig. 19(a) and the 3D path is given as

Eq. (34)

x (r, θ, z) = [\begin{matrix} \tan (θ_{x}) (z - L) + r_{x} (\frac{z}{L}) \\ \tan (θ_{y}) (z - L) + r_{y} (\frac{z}{L}) \\ z \end{matrix}],

for

0 \leq z \leq L

. Note that counterclockwise angles in

θ = [θ_{x}, θ_{y}]^{T}

relative to the principal axis are treated as positive. The magnitude of the derivative with respect to

z

is given as

Eq. (35)

| x^{'} (r, θ, z) | = A (r, θ) = \sqrt{(\frac{r_{x}}{L} + \tan (θ_{x}))^{2} + (\frac{r_{y}}{L} + \tan (θ_{y}))^{2} + 1} .

The plane wave case shown in Fig. 19(b) shows the start of the atmospheric effects at

z = 0

. The 3D paths to the pupil plane in this case are given as

Eq. (36)

x (r, θ, z) = [\begin{matrix} \tan (θ_{x}) (z - L) + r_{x} \\ \tan (θ_{y}) (z - L) + r_{y} \\ z \end{matrix}],

for

0 \leq z \leq L

. The magnitude of the derivative for the plane wave case is given as

Eq. (37)

| x^{'} (r, θ, z) | = A (r, θ) = \sqrt{\tan^{2} (θ_{x}) + \tan^{2} (θ_{y}) + 1} .

Fig. 19

Geometries for the tilt scaling calculations in the (a) spherical wave case and (b) plane wave case. A circular pupil is shown on the bottom, and a source is propagating from top to bottom.

Note that, in either the spherical wave or plane wave case, the magnitude of the derivative term is not a function of $z$ and can be pulled out of the integral in Eq. (33), yielding

Eq. (38)

{OPD}_{λ} (r, θ) = A (r, θ) \int_{z = 0}^{L} (n_{λ} (x (r, θ, z)) - 1) d z .

The phase disturbance function is then expressed as

Eq. (39)

ϕ_{λ} (r, θ) = \frac{2 π}{λ} {OPD}_{λ} (r, θ) = \frac{2 π}{λ} A (r, θ) \int_{z = 0}^{L} (n_{λ} (x (r, θ, z)) - 1) d z .

Therefore, combining Eqs. (31) and (39), the wavelength-dependent tilt is expressed as

Eq. (40)

α_{λ} (θ) = \frac{64}{π D^{4}} \int_{r} r W (r) A (r, θ) \int_{z = 0}^{L} (n_{λ} (x (r, θ, z)) - 1) d z d r .

If we assume that only pressure and temperature changes occur in the atmosphere along the optical paths, then we call upon the relation in Eq. (4) to give us

Eq. (41)

\frac{n_{λ_{1}} (x (r, θ, z)) - 1}{n_{λ_{2}} (x (r, θ, z)) - 1} = \frac{N (λ_{1}) N (p (x (r, θ, z), t (x (r, θ, z))}{N (λ_{2}) N (p (x (r, θ, z), t (x (r, θ, z))} = \frac{N (λ_{1})}{N (λ_{2})} .

This follows because any changes in pressure and temperature that impact the refractivities are not a function of wavelength and will cancel in the ratio. Therefore, using Eqs. (40) and (41), we state that

Eq. (42)

\frac{α_{λ_{1}} (θ)}{α_{λ_{2}} (θ)} = \frac{N (λ_{1})}{N (λ_{2})},

or equivalently

Eq. (43)

α_{λ_{2}} (θ) = \frac{N (λ_{2})}{N (λ_{1})} α_{λ_{1}} (θ),

which is the same as Eq. (10).

Acknowledgments

The authors would like to thank Dr. Steve Fiorino at the Air Force Institute of Technology (AFIT) for providing helpful information related to wavelength dependence of the refractive index structure constant. The authors also thank Dr. Joseph French at Leidos for technical feedback and project management support. We also thank Bruce Wilcoxen and Julie Tollefson with Leidos for project management support. The authors thank the engineering teams at Leidos and MZA Associates for their roles in acquiring the real data used here. This work was supported in part under Air Force Research Laboratory (AFRL) Award Nos. FA8650-18-F-1710 and FA8650-21-F-1125. This paper is approved for public release under Case No. AFRL-2022-2971. The authors have no relevant financial interests in the manuscript and no other potential conflicts of interest to disclose.

Code, Data, and Materials Availability

The truth images from Sec. 4.2 come from a publicly available dataset acquired with the NASA/JPL AVIRIS sensor. The orthorectified data are from September 30, 2011, over Fisherman’s Wharf in Monterey, California. These data can be obtained at the AVIRIS Data Portal.²⁴ The data with turbulence were simulated using the method in Ref. 5, modified as described here. The real multispectral camera data from Sec. 4.3 are not publicly available at this time.

References

1.

D. L. Fried, “Optical resolution through a randomly inhomogeneous medium for very long and very short exposures,” J. Opt. Soc. Am., 56 1372 –1379 (1966). https://doi.org/10.1364/JOSA.56.001372 JOSAAH 0030-3941 Google Scholar

2.

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

3.

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

4.

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

5.

R. C. Hardie et al., “Simulation of anisoplanatic imaging through optical turbulence using numerical wave propagation with new validation analysis,” Opt. Eng., 56 (7), 071502 (2017). https://doi.org/10.1117/1.OE.56.7.071502 Google Scholar

6.

S. L. Lachinova et al., “Comparative analysis of numerical simulation techniques for incoherent imaging of extended objects through atmospheric turbulence,” Opt. Eng., 56 (7), 071509 (2017). https://doi.org/10.1117/1.OE.56.7.071509 Google Scholar

7.

A. Schwartzman et al., “Turbulence-induced 2D correlated image distortion,” in IEEE Int. Conf. Comput. Photography (ICCP), 1 –13 (2017). https://doi.org/10.1109/ICCPHOT.2017.7951490 Google Scholar

8.

N. Chimitt and S. H. Chan, “Simulating anisoplanatic turbulence by sampling intermodal and spatially correlated Zernike coefficients,” Opt. Eng., 59 (8), 083101 (2020). https://doi.org/10.1117/1.OE.59.8.083101 Google Scholar

9.

M. A. Hoffmire et al., “Deep learning for anisoplanatic optical turbulence mitigation in long-range imaging,” Opt. Eng., 60 (3), 033103 (2021). https://doi.org/10.1117/1.OE.60.3.033103 Google Scholar

10.

C. J. Carrano, “Anisoplanatic performance of horizontal-path speckle imaging,” 14 –27 (2003). https://doi.org/10.1117/12.508082 Google Scholar

11.

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). https://doi.org/10.1117/1.OE.53.4.043109 Google Scholar

12.

R. C. Hardie et al., “Block matching and Wiener filtering approach to optical turbulence mitigation and its application to simulated and real imagery with quantitative error analysis,” Opt. Eng., 56 (7), 071503 (2017). https://doi.org/10.1117/1.OE.56.7.071503 Google Scholar

13.

R. C. Hardie et al., “Fusion of interpolated frames superresolution in the presence of atmospheric optical turbulence,” Opt. Eng., 58 (8), 083103 (2019). https://doi.org/10.1117/1.OE.58.8.083103 Google Scholar

14.

R. C. Hardie et al., “Application of tilt correlation statistics to anisoplanatic optical turbulence modeling and mitigation,” Appl. Opt., 60 G181 –G198 (2021). https://doi.org/10.1364/AO.418458 APOPAI 0003-6935 Google Scholar

15.

L. R. Burchett and S. T. Fiorino, “Wavelength correction of refractivity variation measurements,” Opt. Express, 21 31990 –31997 (2013). https://doi.org/10.1364/OE.21.031990 OPEXFF 1094-4087 Google Scholar

16.

S. Basu et al., “Estimation of temporal variations in path-averaged atmospheric refractive index gradient from time-lapse imagery,” Opt. Eng., 55 (9), 090503 (2016). https://doi.org/10.1117/1.OE.55.9.090503 Google Scholar

17.

B. Edlén, “The refractive index of air,” Metrologia, 2 (2), 71 (1966). https://doi.org/10.1088/0026-1394/2/2/002 MTRGAU 0026-1394 Google Scholar

18.

P. E. Ciddor, “Refractive index of air: new equations for the visible and near infrared,” Appl. Opt., 35 1566 –1573 (1996). https://doi.org/10.1364/AO.35.001566 APOPAI 0003-6935 Google Scholar

19.

R. K. Tyson, Adaptive Optics Engineering Handbook, (1999). Google Scholar

20.

S. Basu, J. E. McCrae and S. T. Fiorino, “Estimation of the path-averaged atmospheric refractive index structure constant from time-lapse imagery,” 94650T (2015). https://doi.org/10.1117/12.2177330 Google Scholar

21.

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

22.

M. A. Rucci, R. C. Hardie and R. K. Martin, “Simulation of anisoplanatic lucky look imaging and statistics through optical turbulence using numerical wave propagation,” Appl. Opt., 60 G19 –G29 (2021). https://doi.org/10.1364/AO.427716 APOPAI 0003-6935 Google Scholar

23.

J. Bentley and C. Olson, Field-Guide-to-Lens-Design: Point Spread Function and Strehl Ratio, SPIE(2012). Google Scholar

24.

, “AVIRIS data portal,” https://aviris.jpl.nasa.gov/dataportal/ Google Scholar

25.

B. D. Lucas and T. Kanade, “An iterative image registration technique with an application to stereo vision,” in Proc. 7th Int. Joint Conf. Artif. Intell. – Vol. 2, IJCAI’81, 674 –679 (1981). Google Scholar

26.

K. A. Winick and D. V. L. Marquis, “Stellar scintillation technique for the measurement of tilt anisoplanatism,” J. Opt. Soc. Am. A, 5 1929 –1936 (1988). https://doi.org/10.1364/JOSAA.5.001929 Google Scholar

27.

J. E. Greivenkamp, Field Guide to Geometrical Optics, FG01 SPIE Press(2004). Google Scholar

28.

, “Line integral — Wikipedia, the free encyclopedia,” https://en.wikipedia.org/wiki/Line_integral Google Scholar

Biography

Russell C. Hardie is a full professor in the Department of Electrical and Computer Engineering at the University of Dayton and holds a joint appointment in the Department of Electro-Optics and Photonics. He received the University of Dayton’s top university-wide teaching award in 2006 (the alumni award in teaching). He also shared the Rudolf Kingslake Medal and Prize from SPIE in 1998 for work on super-resolution imaging. His research interests include a wide range of topics in signal and image processing, image restoration, medical image analysis, and machine learning.

Michael A. Rucci received his BS and MS degrees in electrical engineering from the University of Dayton in 2012 and 2014, respectively. He is a research engineer at the Air Force Research Laboratory, Wright-Patterson AFB Ohio, USA, and his current research interests include day/night passive imaging, turbulence modeling and simulation, and image processing. He and three other researchers received the Rudolf Kingslake Medal from SPIE in 2018 for their work in turbulence modeling.

Santasri Bose-Pillai received her PhD in electrical engineering with a focus in optics in 2008 and her MSEE degree in 2005, both from New Mexico State University. She received her BSEE degree with honors from Jadavpur University, India, in 2000. Currently, she is a research assistant professor in AFIT’s Engineering Physics Department. Her research interests include propagation and imaging through atmospheric turbulence, telescope pointing, and tracking and synthesis of partially coherent sources. She is a senior member of Optica (formerly OSA) and SPIE and a regular member of DEPS.

Richard Van Hook is the deputy chief of the Passive Radio Frequency Sensing Branch at the Air Force Research Laboratory, Wright-Patterson AFB, Ohio, USA. He received his MS and BS degrees in computer engineering from Wright State University in 2014 and 2008, respectively, and his PhD in electrical engineering from the University of Dayton in 2021. His research interests include image processing, atmospheric modeling, and turbulence mitigation.

Barry K. Karch is a principal research electronics engineer in the Multispectral Sensing and Detection Division, Sensors Directorate of the Air Force Research Laboratory, Wright-Patterson AFB Ohio, USA. He received his BS degree in electrical engineering in 1987, his MS degree in electro-optics and electrical engineering in 1992/1994, and his PhD in electrical engineering in 2015 from the University of Dayton, Dayton, Ohio, USA. He has worked in the areas of EO/IR remote sensor system and processing development for more than 30 years.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 International 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, Michael A. Rucci, Santasri R. Bose-Pillai, Richard Van Hook, and Barry K. Karch "Modeling and simulation of multispectral imaging through anisoplanatic atmospheric optical turbulence," Optical Engineering 61(9), 093102 (10 September 2022). https://doi.org/10.1117/1.OE.61.9.093102

Received: 23 June 2022; Accepted: 18 August 2022; Published: 10 September 2022

Access the abstract

JOURNAL ARTICLE
28 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY