The problem of imaging in the presence of degrading optical turbulence conditions has been an ongoing area of active research for many years. Much recent research has focused on correcting imagery through image processing (dewarping and/or deconvolution) techniques.18.104.22.168.6.7.–8 Paralleling this inverse problem is the direct simulation of turbulence on clear images to create artificially degraded imagery at known turbulence levels910.11.–12 for benchmarking purposes. A third research area involves system intercomparisons for cost-effectiveness evaluation.1314.–15 In each case, undergirding the research is the core issue of establishing the baseline effect of turbulence on image quality, particularly for short-exposure (SE) imaging.
The standard method used by systems performance modelers1314.15.16.–17 had been Fried’s SE modulation transfer function (MTF).1819.–20 For example, Holst15 devotes an entire chapter to Fried’s SE MTF. This model was the common choice due to its simplicity, despite the existence of more accurate (and numerically intensive) techniques (e.g., Ref. 21).
Critics of Fried’s approach pointed to its high-frequency response failure when modeling high turbulence degradation influences.2122.–23 Nonetheless, the simplicity of Fried’s result suggested a re-examination of Fried’s approach might yield new findings. The resulting study24 examined the properties of an overlooked tilt-phase correlation term and developed improved expressions for diffraction influences on the turbulent phase structure function (PSF), angle-of-arrival variance, and therefore on the SE MTF. However, the numerical integration technique used was not optimal for evaluating the tilt-phase influence at combined high-angular frequency and high-turbulence conditions. Consequently, a modeling compromise was adopted based on Fried’s quadratic correction to the long-exposure (LE) MTF. This approach was subsequently critiqued by Charnotskii2526.27.–28 who argued its quadratic term would produce inaccurate results at high frequency.
This issue is addressed here through a new analysis of high-frequency turbulence effects and development of a new SE MTF analytic model. The new model exhibits diffraction-limited behavior at all frequencies while applying to a wider range of optical scenarios. The resulting model is compared with published calculations from a path-integral technique21 and a Markov method.26 The model’s Rytov approximation (RA)29 is also described, and criteria developed by Tatarski29 and Dashen30 are examined to determine the extent of validity of modeled conditions.
The remainder of the paper is designed as follows: Section 2 describes the propagation model used, its justification, and numerical processing techniques developed to handle high-turbulence conditions. Section 3 describes the analysis process used to translate the computed database into the new analytical SE MTF expression. Section 4 presents comparisons with alternative approaches, followed by conclusions in Sec. 5.
Propagation and Imaging Models
The propagation model used is based on the Rytov approximation,16,29 involving a multiplicative turbulent scattering effect, using the imaging geometry of Fig. 1. Monochromatic incoherent light of wavelength emerges from the transverse object plane at . Photons pass through the system aperture’s lens at to reach the imaging plane at . The system’s thin lens optics are considered focused on the object plane. The main propagation axis is denoted by ; transverse two dimensional (2-D) vectors are denoted by .
In general, light from a point is evaluated for scatter off the ensemble of turbulent fluctuations at , passage through the aperture at and its effect on the light at image point . Each object point will have a reciprocal image point at following a chief ray through the center of the system entrance pupil. Due to the incoherent nature of the light, each photon travels independently, emitted as a quantum event, and arriving at the image plane based on conditional passage through the system entrance pupil. The instantaneous point spread function represents the electro-magnetic equivalent of the photon’s quantum probability of arrival at image point , given the current optical turbulence pattern (considered frozen) in the region .
The Rytov approximation considers a propagating scalar wavefunction
Turbulent amplitude () and phase () effects are imposed on the propagating wave via weak scattering from position-varying refractive index fluctuations26
Using Eq. (1)’s model propagating wave, the RA method’s single-scattering approximation is used to evaluate and instantaneous amplitude and phase perturbations for radiation emerging from a source at at system entrance pupil points, ,
This equation reflects scattering from turbulence at points , , using free-space Green’s function propagators16
Modulation Transfer Function
Following the methodology of Ref. 24, we analyze the wave perturbations at the system’s entrance pupil associated with the instantaneous SE MTF for the source point
Here, the field variables reflect the situation just after truncation by the aperture
According to Fried’s SE theory, the instantaneous turbulent phase tilt is removed, based on variable denoting the linear moment of the instantaneous phase perturbation
Vector variable is a normalized angular frequency, where denotes the diffraction-limited maximum angular frequency16 passed by a circular aperture, such that . is the entrance pupil diameter, dimensionalizing temporarily, before is used to transform integral (10) to dimensionless form.
The mean SE MTF sought is then evaluated by taking the expectation of Eq. (7), and16 analysis (pp. 395–407).
For brevity of presention, only results are shown here, based on more complete descriptions in Refs. 20 and 24. First, amplitude and phase perturbations are considered independent, permitting the expectation of the amplitude term to be evaluated as31 The remaining three terms are phase related, involving phase variance, tilt variance, and tilt-phase terms. The second term is the phase variance term, appearing as Ref. 24. This combined term quantifies the turbulence LE effect. While the instantaneous perturbation effects were source specific (), therefore anisoplanatic, satisfying an energy conservation constraint,32 the second-order WSF is source independent, evaluated by weighted path integral of ,31 therefore isoplanatic.
Assuming constant Kolmogorov turbulence throughout the following20,24
The remaining two terms reflect SE corrections to the LE MTF (due to ). The tilt variance term that can be written24Ref. 24).24 24
When turbulence is low (), the turbulent exponential argument approaches zero, and the SE MTF approaches the behavior of the system MTF,
When , the Eq. (11) integral was evaluated numerically using the form24 evaluation of was handled using a relatively inefficient integration technique, but here, based on the PSF of Eq. (18), was removed from the exponential integral, allowing the remaining term
Third, the evaluations were modified to handle the impact of large effects inside the exponential through use of an offset factor
By adjusting , math overflow errors in evaluating the integral could always be eliminated. Thereafter, calculations could be interpreted in the form
Combining this term with the tilt variance effect of Eq. (16) yields the function
Calculations of the function were performed for a wide range of conditions: from to by factors of 2 (9 steps), from to in increments of tenths of decades (41 steps), and from 0.01 to 0.99 in steps of 0.01, and four additional steps at (103 steps), for a database of 38,007 calculations. The analysis of this database is described in Sec. 3, using these results to formulate a new analytic SE MTF.
Sections 5 through 7 of Ref. 24 provides further details of the analysis of Fried’s approach.
Range of Rytov Approximation Validity
Before examining the database, some discussion of the extent of validity of the Rytov method seems appropriate. Clearly, the Rytov approximation is applicable to coherent radiation studies under weak turbulence conditions, but is typically considered inappropriate for stronger turbulence scenarios, when the Rytov variance21,26 in Sec. 4, where comparable results are achieved at Rytov variance values seemingly beyond the capability of the RA, one may question the validity of such restrictions for incoherent radiation, since temporal and spatial correlations do not exist between propagating photons.
Tatarski’s analysis29 only required that , and that perturbation field gradients be small relative to wave number , conditions easily met. Similarly, Dashen30 developed full and partial saturation conditions for judging RA validity, formed as ratios of a turbulence parameter, , to a diffraction parameter, .
For the full saturation case, Dashen used,33 based on Kaimal’s 1968 Kansas experiment analysis,34 . For the turbulence parameter, Dashen used a phase variance metric also based on an outer-scale influenced covariance
More stringent conditions were associated with a partial saturation case involving multiple scattering scales. Dashen considered the effect of imposing a scale size cutoff that would affect both turbulent and diffraction effects. For the SE imaging problem, the most natural cutoff is the finite aperture. Dashen’s , corresponding to an RMS phase variance, is equivalent to one of Fried’s18 mean aperture phase variance terms, where the order indicates the number of phase perturbation terms corrected during the imaging process. In SE imaging, piston, tip, and tilt phase-expansion terms do not alter image quality (corresponding to Fried’s order-L case), leading to the condition,
Setting , for a given , a maximum value results
This condition translates to a maximum Rytov variance function of
Typical terrestrial imaging scenarios involve , implying the RA remains valid up to high-turbulence conditions.
Analytic Model of Short-Exposure Atmospheric Modulation Transfer Function
Using the previous section’s data set, this section describes development of an analytic expression for . This is possible due to significant correlations that exist between results at different values. This behavior is illustrated in Fig. 2 for two sets of plots at values of and 16, for a series of increasing values of .
The key feature of these curves is that for moderate values often encountered for typical terrestrial imaging scenarios, exceeds unity at low to moderate frequencies, resulting in enhanced behavior versus the LE response.
In developing a model of this behavior, it was first observed that , permitting to be expressed
The maximum value of (as ) must also be a function of , . and are plotted in Fig. 3.
Figure 3 curves both suggest asymptotic behaviors at both extremes of . These behaviors can be modeled using a sigmoidal function (a variant of the hyperbolic tangent)
To facilitate the modeling process, a rescaled and shifted sigmoid was also introduced
The center of this function occurs at with an overall shift of between asymptotes, and a central derivative of .
The function can be expressed using as
However, exhibits different asymptotic behaviors at large and small . Therefore, a third sigmoidal form was developed, featuring a central splice (piecewise sigmoidal):39) and (41) are plotted along with database derived points in Fig. 3. Variables and will frequently appear hereafter as surrogates for variables and in arguments to sigmoidal functions, as the majority of the modeled behaviors exhibit logarithmic dependence.
Figure 4 illustrates plots of for the same families of curves ( and 16) as in Fig. 2, highlighting the similarities in behaviors at different values. These plots also reveal how develops increasingly complex behaviors as increases.
To model these behaviors, it was observed that function curves at high maintained simpler forms with increasing than at low , suggesting the functions could be divided at the peak and studied separately at high and low frequencies. The next step was, therefore, to model the form of the peak curves (as illustrated in Fig. 4) as functions of and .
Note that the peak curves at different values exhibit similar dependencies. Thus, the behaviors of the different peak curves could be modeled as shifted versions of the peak curve for , described by (abscissa) and (ordinate) functions. The abscissa curve is modeled as
Several curves at different levels illustrate this shifting behavior in Fig. 5. A pair of analytic expressions that approximate this behavior are given by
Several ordinate function curves, designated (using the non-normalized form), are illustrated for varying in Fig. 6, for values ranging from to 16. This component exhibits both shifting and scaling properties, but again was quantified relative to its behavior at , using
The behavior at was then rescaled using , and shifted using , in the form,
The function supports a further parameterization through introduction of47) and (39), where . The curve shapes appear to vary based on height , while dependence may be modeled using a Legendre-polynomial-related series expansion.
The Legendre polynomials begin as , , and use recurrence relation35
Using this basis set, each curve could be separated into two functions at the peak position , followed by affine transformation to stretch each monotonic function into its own half-range, and reverse copied to fill the opposing half-range. Thereby, two functions, each even about , could be analyzed using Legendre-derived basis functions. Due to their symmetry, only even-order series terms were needed to characterize each function . The model could, thus, be splined using the coefficient sets: for , and for , as
Sets of coefficient behaviors are illustrated in Fig. 7. behaviors are plotted in Fig. 8. Coefficients through , and , exhibit approximately hyperbolic dependence. The remaining coefficients exhibit approximately linear dependence.
Hyperbolic behaviors of and coefficients were expressed using constants through in the models
SPGD hyperbolic Cm and Dm coefficients.
SPGD linear Cm and Dm coefficients.
With the development of the and functional forms, the new analytic SE MTF model form was completed. However, given the number of free constants (80) involved, an open question remained whether a slightly modified set of coefficients could produce a better fit. To explore this possibility, a stochastic parallel gradient descent (SPGD) procedure36 was applied to improve the fit to the database.
The coefficient set derived in the initial analysis (described above) was used to “seed” the SPGD optimizer. The SPGD consisted of randomly perturbing the full 80-parameter coefficient set, recomputing the RMS error for the full database, and comparing that error with the previous best fit error. Any new set that outperformed the prior best fit was adopted as the new best fit. Use of the SPGD reduced the overall RMS error of from 0.006795 to 0.002728. For simplicity of presentation, coefficients listed in Tables 1 and 2 represent the SPGD best-fit coefficients. For the remaining expressions, the original results (used in the plotted figures) were retained in Eqs. (39) through (49), so curves in the various figures can be replicated. Optimized versions of these equations are given below.
To summarize, the new SE MTF model is formulated starting with the general expression in Eq. (27), then using from (20), from (36), from (42), expressions in Eqs. (44), (47), and (50) through (55), and results of the SPGD analysis
Comparisons with Alternative Approaches
The analytic model developed in the previous section is compared here with results obtained by Charnotskii21,26 and several prior analytic methods. Charnotskii21 used a Feynmann-path integral technique.
Figures 9 through 11 intercompare the new model with Charnotskii’s21 Figs. 8 through 10. Charnotskii’s parameter is related to via , such that values 0, , , , and correspond to values 0.000, 1.241, 2.851, 4.638, and 6.550, respectively. Charnotskii’s parameter is related to as . Reported values of the three plots were 10.0, 1.0, and 0.1, respectively, corresponding to values 2.523, 0.798, and 0.252. The curves (the uppermost line in each graph) correspond to the system MTF, and match identically. Remaining lines in Figs. 9 and 11 were plotted using values that produced the closest fit to Charnotskii’s results, in Fig. 9 and in Fig. 11. Lines in Fig. 10 were constructed based on the reported (). Systematic deviations employed in these plots appear consistent with Charnotskii comments that and constituted limiting cases, but might be due to approximations used when Charnotskii developed Eq. (31) of Ref. 21.
The comparisons of Figs. 9Fig. 10–11 provide partial verifications of both the present method and Charnotskii’s calculations. A similar comparison is made in Fig. 12 between data from Fig. 2 of Ref. 26 (symbols) and simulations (lines) based on,27), transformed into distances in the entrance pupil. Results similar to Charnotskii’s calculations were achieved using values 0.02, 0.21, 0.78, and 12.0, and (corresponding to ). Similar stretching as in Figs. 9 and 11 was again observed. Curiously, the lowest line (second lowest ) exhibited the worst fit to Charnotskii’s results, suggesting discrepancies between the methods are not explicable simply due to RA failure.
A key element of these results is the range of Rytov variance where similar results were obtained using RA versus the alternative approaches of Refs. 21 and 26. While Figs. 9 and 11 involved stretching, a single adjustment appeared to correct all cases, while no adjustment was used in Fig. 10 where the two strongest turbulence cases corresponded to values of 12.7 and 22.8. This ability of the RA to perform at high levels is consistent with the Sec. 2.3 analysis.
Lastly, four prior SE MTF models, (1) (Fried’s far-field case), (2) (Fried’s near-field case), (3) Holst’s15 method where either Fried’s near-field or far-field approximation was chosen based on best fit, and (4) the method of Ref. 24, were compared with the new analytic SE MTF model (5). An RMS metric was used to compare each estimate against database-derived MTF values, weighted by frequency . The RMS results computed were , , , , and . The prequel model (4) outperformed methods (1) through (3), but the new model offered dramatic improvement (by a factor of 29) over the prequel, and greater than 100 better than the far-field approach.
This paper’s extended high angular frequency analysis, using the methodology described in Sec. 2, has lead to a new analytic expression for the atmospheric SE MTF. This model’s RMS error of 0.000218 versus numerically integrated results represents factors of 29 improvement over the prequel,24 49 over Holst,15 and 64 over Fried’s near-field case.20 The new model exhibits diffraction-limited behavior at all angular frequencies.
The primary focus of this effort was to provide systems performance engineers with an easily computed SE MTF that accounted for a thorough range of system and observational scenario conditions. Specifically, it improves evaluation of SE MTF degradation at high frequencies. This will aid trade studies in evaluating image reconstruction technique effectiveness where high frequency edge information losses are of critical concern.
Augmenting the present model, a previous study37 showed that the new model, based on structure functions, can treat path-varying turbulence effects, including slant-path and valley geometries. A planned future paper will expand on these results, considering the annular aperture case (Newtonian telescopes, etc.), where the circular aperture results represent the zero central obscuration limit of the more general solution.
The main caveats of the model regard the Kolmogorov spectrum used and the validity of the RA method. Regards the spectrum, we impose the requirement that , ensuring occurs in the inertial subrange. Regarding the RA method, conditions derived based on Dashen’s 30 analysis [Eqs. (31) and (35)] suggest a relatively wide range of applicability. This suitability was further tested through intercomparison with path-integral21 and Markov26 techniques.
Results of the current study might also be used in an improved inverse filtering procedure [e.g., Ref. 5], though such a study is beyond the scope of the present work. In such a study, multiple range effects could be simulated, but special handling would be required to subsegment the images to apply different MTF’s in different portions of the image field.
David H. Tofsted is a physicist with the U.S. Army Research Laboratory at White Sands Missile Range, New Mexico, with over 33 years in service. He received his BS degree in physics from Pennsylvania State University, Pennsylvania, in 1979, and MS and PhD degrees in electrical engineering from New Mexico State University, New Mexico, in 1992 and 2013. His research has focused on optical propagation problems, including numerical simulations of optical turbulence in the atmospheric boundary layer and its impacts on optical imaging systems. He has authored over 60 publications.