Modeling and simulation of multispectral imaging through anisoplanatic atmospheric optical turbulence

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.


Introduction
Atmospheric optical turbulence is a major source of image degradation in long range imaging. 1,2 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. 2 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 3 and anisoplanatic turbulence simulations. 4,5 A comparison of numerical simulation techniques is presented in Ref. 6. A fast warping-only anisoplanatic simulator is presented in *Address all correspondence to Russell C. Hardie, rhardie1@udayton.edu 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. 9 Such simulators have been used in many other works [10][11][12][13][14] 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.

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

Refractive Index Structure Function
The Kolmogorov atmospheric turbulence model gives rise to the refractive index structure function. 2 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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 2 8 8 D n;λ ðrÞ ¼ h½n λ ðx þ rÞ − n λ ðxÞ 2 i ¼ C 2 n;λ r 2∕3 ; where x and r are three-dimensional (3D) spatial coordinate vectors and r ¼ jrj. The wavelength-dependent atmospheric index of refraction is given by n λ ð·Þ, and the notation h·i represents an ensemble mean operator. The turbulence strength in Eq. (1) is captured by the refractive index structure parameter, C 2 n;λ , in units of m −2∕3 . 2 To represent variability along the optical path, this parameter is expressed with the function C 2 n;λ ð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. 15,16 This is expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 1 4 2 C 2 n;λ 2 ðzÞ ¼ where the atmospheric refractivity is given by NðλÞ ¼ ½n λ − 1 × 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 17 is given as (3) 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. 17 Edlén provides a multiplicative scaling factor to account for other pressures and temperatures for dry air. This is expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 6 5 3 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. 18 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 ∼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.

Fried Parameter
The wavelength-specific atmospheric coherence diameter (or Fried parameter) 1,3 is expressed as a weighted integral of C 2 n;λ ðzÞ yielding E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 3 5 7 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 4 7 7 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).

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 2 Z ðλÞ ¼ hα λ ðθÞ · α λ ðθÞi. Tyson 19 provides an expression for the two-axis Z-tilt variance in terms of the Fried parameter as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 7 1 0 where D is the aperture diameter. Equivalently, this is expressed in terms of the refractive index structure parameter as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 4 1 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 5 4 8 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 2 n;λ . Because the dependence on C 2 n;λ is the same for tilt correlation and differential tilt variance, 5,20 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 4 0 3 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. 7 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 singleband 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.

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. 2,3 The wavelength-dependent version of the isoplanatic angle is expressed as Using steps similar to those in Appendix B, we express the isoplanatic angle wavelength scaling as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 6 7 6 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.

Log Amplitude Variance
The last statistic that we consider in this section is the log-amplitude variance. 3 This statistic reflects fluctuations in the amplitude of the wave function in the pupil plane and is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 5 3 4 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 6 ; 4 6 5 The scaling term in Eq. (14) is plotted in Fig. 4 for a fixed reference wavelength of λ 1 ¼ 0.88 μm.

Optical Transfer Functions
Consider an OTF model that includes the effects of both diffraction and the atmosphere. Such an OTF is expressed as where ρ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and u and v are the spatial frequencies in units of cycles per unit distance. The atmospheric OTF model developed by Fried 1,2 is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 6 ; 7 0 9 H atm ðρ; λÞ ¼ exp −3.44 λlρ 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. 12 This parameter has been used to model the level tilt variance correction achieved when image registration and fusion are employed. [12][13][14] 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 wellknown model 21 to explicitly show the wavelength dependence yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 6 ; 5 6 9 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 6 ; 4 3 9 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 wellvalidated simulation results in prior papers. 5,22 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. 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 2 n;λ 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. 23 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 diffractionlimited 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). Table 1 Optical parameters used for the example and simulation results.  Table 1 for two turbulence levels: (a) C 2 n;0.88 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.

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  Table 1. The PSF peak is divided by the peak of the diffraction-limited PSF for (a) each wavelength and (b) the shortest wavelength.
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. 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. 5 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.

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. 5 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. 22 In that work, a sampling rule is proposed that sets the phase screen sample spacing as follows: 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. 3,22 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. 3,22 We found that we achieve a high level of agreement with the theoretical multispectral statistics studied here using a g ¼ 5 and c ¼ 2. Fig. 7 Phase screen geometry for the anisoplanatic numerical wave propagation simulation method. 5 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.
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 3 and also used in Ref. 22 is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 6 ; 6 9 9 The number of samples, N ðλÞ, must be such that N ðλÞ ≥ 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.
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).

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. 5 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. 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 2 n;0.88 μm ¼ 1 × 10 −16 m −2∕3 , 5 × 10 −16 m −2∕3 , and 10 × 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. 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 2 n;0.88 μm values of 1 × 10 −16 m −2∕3 , 5 × 10 −16 m −2∕3 , and 10 × 10 −16 m −2∕3 . Quantitative validation results for these PSFs are summarized in Tables 3-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 oneaxis 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.

Multispectral PSF Analysis
Note in Tables 3-5 that the wavelength-dependent C 2 n;λ 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-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 2 n;0.88 μm ¼ 5 × 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. 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  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). 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.

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. 24 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 2 n;0.88 μm values of 1 × 10 −16 m −2∕3 , 5 × 10 −16 m −2∕3 , and 10 × 10 −16 m −2∕3 . The increased blurring may be observed for the higher values of C 2 n;0.88 μm . 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. 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).
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-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. 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 5,14 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). 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 2 n;0.88 μm ¼ 5 × 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.

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 gradientbased image registration method. 25 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 2 n;λ ðzÞ profile and using an imagesized patch. These theoretical values are compared with the empirical values in Table 9, from which we see a good level of agreement.

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.

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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 1 1 6 ; 3 0 0 h½ðn λ ðx þ rÞ − 1Þ − ðn λ ðxÞ − 1Þ 2 i ¼ C 2 n;λ r 2∕3 : Next, the spectral ratio by which we cancel common terms is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 6 ; 2 4 8 C 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 1 1 6 ; 1 7 0 C 2 n;λ 2 C 2 n;λ 1 ¼ h½Nðλ 2 ÞNðpðx þ rÞ; tðx þ rÞÞ − Nðλ 2 ÞNðpðxÞ; tðxÞÞ 2 i h½Nðλ 1 ÞNðpðx þ rÞ; tðx þ rÞÞ − Nðλ 1 ÞNðpðxÞ; tðxÞÞ 2 i : Factoring out the wavelength-dependent terms of the refractivities yields Finally, canceling the common pressure and temperature terms, we obtain the desired result E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 1 1 6 ; 7 2 3 C 2

Appendix B: Fried Parameter Scaling Derivation
The Fried parameter as a function of wavelength λ 2 is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 1 1 6 ; 6 3 4 Using Eq. (2), we express this in terms of the refractive index structure parameter for wavelength λ 1 defined as C 2 n;λ 1 ðzÞ. This is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 1 1 6 ; 5 6 3 Now we separate the terms that define r 0 ðλ 1 Þ from scaling constants to produce E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 8 ; 1 1 6 ; 5 0 5 Noting that the second term in Eq. (28) equals r 0 ðλ 1 Þ, we write E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 9 ; 1 1 6 ; 4 4 8 Finally, combining the exponents and removing the negative by inverting the argument, we obtain 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 27 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 ds is interpreted as an elementary arc length. 28 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. 28 This is expressed as where xðr; θ; zÞ is the 3D optical path as a function of z and x 0 ðr; θ; zÞ ¼ ∂xðr; θ; zÞ∕∂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 (34) for 0 ≤ z ≤ 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 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 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 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 1 ; 1 1 6 ; 6 0 0 n λ 1 ðxðr; θ; zÞÞ − 1 n λ 2 ðxðr; θ; zÞÞ − 1 ¼ Nðλ 1 ÞNðpðxðr; θ; zÞ; tðxðr; θ; zÞÞ Nðλ 2 ÞNðpðxðr; θ; zÞ; tðxðr; θ; zÞÞ ¼ 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 2 ; 1 1 6 ; 5 1 8 or equivalently E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 4 3 ; 1 1 6 ; 4 6 2 α λ 2 ðθÞ ¼ Nðλ 2 Þ Nðλ 1 Þ α λ 1 ðθÞ; which is the same as Eq. (10).