Open Access
25 April 2014 Extended high-angular-frequency analysis of turbulence effects on short-exposure imaging
Author Affiliations +
Abstract
An improved analysis of optical turbulence effects on short-exposure passive (SE) imaging is described, resulting in a new analytic expression for the SE modulation transfer function (MTF). This analysis expands on a 2011 study that examined characteristics of a tilt-phase component discovered in the standard theory of SE turbulence effects characterization. The analysis introduces an improved integration technique and a reformulated phase structure function, facilitating computation of a 38,007 element database of MTF results at low- to high-angular frequencies covering a wide range of diffraction and turbulence conditions. Analysis of this database is described, yielding a new analytic SE MTF, accurate to a root-mean-square error of 0.000218 versus the database. Comparisons show that the new expression is well correlated to an alternative computationally intensive method, and it is a factor 29 to 64 improvement over prior analytic expressions. Limits of applicability of the approach for incoherent imaging are also discussed. The low-computational cost of the new method is suitable for systems performance modeling of turbulence impacts, including path-varying turbulence scenarios.

1.

Introduction

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 Paralleling this inverse problem is the direct simulation of turbulence on clear images to create artificially degraded imagery at known turbulence levels912 for benchmarking purposes. A third research area involves system intercomparisons for cost-effectiveness evaluation.1315 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 modelers1317 had been Fried’s SE modulation transfer function (MTF).1820 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.2123 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 Charnotskii2528 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.

2.

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 z=L. Photons pass through the system aperture’s lens at z=0 to reach the imaging plane at z=R. The system’s thin lens optics are considered focused on the object plane. The main propagation axis is denoted by z; transverse two dimensional (2-D) vectors are denoted by x.

Fig. 1

Imaging scenario geometry. 2-D vectors x oriented transverse to main z optical axis.

OE_53_4_044112_f001.png

In general, light from a point rO=(xO,L) is evaluated for scatter off the ensemble of turbulent fluctuations at (x,z), passage through the aperture at (xE,0) and its effect on the light at image point (xI,R). Each object point (xo,L) will have a reciprocal image point at (xi,+R) 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 xI, given the current optical turbulence pattern (considered frozen) in the region L<z<0.

2.1.

Propagation Model

The Rytov approximation considers a propagating scalar wavefunction

Eq. (1)

UU0exp(+iϕ),
consisting of a zeroth-order free-space solution, U0, plus a complex turbulence multiplier. U0 contains the primary longitudinal phasor exp[ik(z+L)], where i2=1 and k=2π/λ is the wavenumber, plus the phase information necessary to determine the zero-turbulence focal point of the source point in the image plane, and spherical wave phase information that will be focused by the system lens effect.

Turbulent amplitude () and phase (ϕ) effects are imposed on the propagating wave via weak scattering from position-varying refractive index fluctuations

Eq. (2)

n1(x,z)=n(x,z)/n1,
where n1 is the expectation (mean) of the refractive index n(x,z). Turbulent fluctuations will be modeled (exceptions noted) assuming Kolmogorov turbulence, characterized by the power spectrum26

Eq. (3)

Φn(Cn2,κ)=0.033Cn2κ11/3,
where κ is the turbulent spatial frequency m1. For natural turbulence, n10.001, such that the refractive index structure parameter, Cn2, is typically less than 1011m2/3, leading to forward scattering (in the +z direction), and use of the parabolic wave equation

Eq. (4)

2ikzU+2U+2k2n1U=0,
with z the on-axis derivative operator, and 2 the transverse Laplacian.

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 rO=(xO,L) at system entrance pupil points, rE=(xE,0),

Eq. (5)

[(rO;xE)ϕ(rO;xE)]=[ReIm]{k22πTVG(rO;r)n1(r)G(r;rE)G(rO;rE)dr}.

This equation reflects scattering from turbulence at points r=(x,z), L<z<0, using free-space Green’s function propagators16

Eq. (6)

G(r1;r2)=1(z2z1)exp{ik[(z2z1)+|x2x1|22(z2z1)]},
with ri=(xi,zi) and z2>z1, to propagate the field from the source, through the scattering point, to the observation point. The Green’s function in the denominator normalizes the results according to the free space effect. Note, particularly, that and ϕ are source and observation point dependent quantities; therefore they are anisoplanatic.

2.2.

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 rO

Eq. (7)

M^S(ω,rO)=4πdxEU˜(rO;xE)U˜*(rO;xEDω).

Here, the field variables reflect the situation just after truncation by the aperture

Eq. (8)

U˜(rO;xE)=WD(xE)exp{(rO;xE)+i[ϕ(rO;xE)a(rO)·xE]},
and have had the mean U0 wavefront tilt and curvature removed. The system entrance pupil is modeled using a uniformly transparent window function

Eq. (9)

WD(x)={1,|x|<D/2,0,|x|>D/2.

According to Fried’s SE theory, the instantaneous turbulent phase tilt is removed, based on variable a denoting the linear moment of the instantaneous phase perturbation

Eq. (10)

a(rO)=64πD4dxEWD(xE)ϕ(rO;xE)xE.

Vector variable ω=Ω/Ω0 is a normalized angular frequency, where Ω0=D/λ denotes the diffraction-limited maximum angular frequency16 passed by a circular aperture, such that 0<|ω|<1. D is the entrance pupil diameter, dimensionalizing ω temporarily, before u=xE/D is used to transform integral (10) to dimensionless form.

The mean SE MTF sought is then evaluated by taking the expectation of Eq. (7), and ω=|ω|

Eq. (11)

MS(ω)=M^S(ω,rO)=4πdxWD(x)WD(xDω)exp(P0P1P2+P3),
where the turbulent effects are expressed as the P0P3 terms. The expectation operator was passed inside the integral, as it affects only the turbulence terms. It is worth noting that this averaging process transforms the anisoplanatic instantaneous MTF into an isoplanatic second-order function, consistent with Goodman’s16 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 as

Eq. (12)

exp(P0)=exp{(x)+(xDω)}=exp{12D(Dω)},
using D, the amplitude structure function.31 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

Eq. (13)

exp(P1)=exp[iϕ(x)iϕ(xDω)]=exp{12Dϕ(Dω)},
where Dϕ(x) is the PSF. This term combines with D(x) to form the wave structure function (WSF), D(x), shown in path-integral form as Eq. (11) of Ref. 24. This combined term quantifies the turbulence LE effect. While the instantaneous perturbation effects were source specific (rO), therefore anisoplanatic, satisfying an energy conservation constraint,32 the second-order WSF is source independent, evaluated by weighted path integral of Cn2(z),31 therefore isoplanatic.

Assuming constant Kolmogorov turbulence throughout the following

Eq. (14)

P0+P1=D(Dω)/2=[D(Dω)+Dϕ(Dω)]/2=[Dω/(r0/2.1)]5/3=(2.1Xω)5/3,
with X=D/r0, using the Fried coherence diameter20,24

Eq. (15)

r0=3.018(k2LCn2)3/5.

The remaining two terms reflect SE corrections to the LE MTF (due to P0+P1). The tilt variance term that can be written24

Eq. (16)

P2=1.0433(2.1Xω)5/3G(Q)ω1/3,
introducing parameter Q=D/F=D/(λL)1/2 to characterize diffraction influences, and function G(Q) that varies from 1/2 to 1; G1 for Q>3 (G(Q) plotted in Ref. 24).

Since P0, P1, and P2 are independent of x, they do not affect the aperture integral of Eq. (11). The remaining tilt-phase correlation, written 24

Eq. (17)

P3(x,ω)=32πD4dxaWD(xa)[xa·(Dω)][Dϕ(|xDωxa|)Dϕ(|xxa|)],
must be evaluated numerically. An improved version of the PSF model is introduced to facilitate evaluation of P3

Eq. (18)

Dϕ(Du)/2=D(Du)α(Qu)/2=[(2.1X/Q)5/3][(Qu)5/3α(Qu)]=f1(X/Q)f2(Qu),
where u=|x|/D is a normalized separation variable, X and Q have absorbed the D factors, and α(Qu) parameterizes the diffraction influence24

Eq. (19)

α(Qu)=01dc(3/8)0dγ[1J0(γc)]1.118γ8/3×cos2{γ2[c(1c)]4π(Qu)2},
using dimensionless path (c=z/L) and frequency (γ=Duκ) variables. Function α(Qu) transitions from 1/2 to 1 over eight orders of magnitude in Qu.24

When turbulence is low (X0), the turbulent exponential argument approaches zero, and the SE MTF approaches the behavior of the system MTF,

Eq. (20)

M0(ω)=4πduW1(u)W1(uω);={(2/π)[cos1(ω)ω1ω2],ω<1,0,ω1,
which depends only on ω, due to azimuthal symmetry.

When X0, the Eq. (11) integral was evaluated numerically using the form

Eq. (21)

MT(Q,X,ω)=4/πM0(ω)duW1(uω)W1(u)×exp[+P3(u,ω,Q,X)],
where the system MTF, a function independent of Q and X, could be divided out to focus on the turbulence-varying component. In the prequel,24 evaluation of MT was handled using a relatively inefficient integration technique, but here, based on the PSF of Eq. (18), f1(X/Q) was removed from the exponential integral, allowing the remaining term

Eq. (22)

F3(Q,u,ω)=32πduaW1(ua)[ua·j^ω][f2(Q|uj^ωua|)f2(Q|uua|)],
to be tabulated and used for interpolation. The dimensionality of the tabulated results was also reduced through application of azimuthal symmetry that permitted the orientation of ω to be fixed along the y-axis (unit vector, j^).

Third, the MT evaluations were modified to handle the impact of large f1(X/Q) effects inside the exponential through use of an offset factor A

Eq. (23)

MT(Q,X,ω)=4exp(A)πM0(ω)duW1(uj^ω)W1(u)exp[f1(X/Q)F3(Q,u,ω)A].

By adjusting A, math overflow errors in evaluating the integral could always be eliminated. Thereafter, MT calculations could be interpreted in the form

Eq. (24)

MT(Q,X,ω)=exp[+(2.1X)5/3N(Q,X,ω)ω2],

Eq. (25)

N(Q,X,ω)=(2.1X)5/3ln[MT(Q,X,ω)]ω2.

Combining this term with the tilt variance effect of Eq. (16) yields the function

Eq. (26)

V(Q,X,ω)=N(Q,X,ω)1.0433G(Q),
such that the SE MTF can be written

Eq. (27)

MS(Q,X,ω)=M0(ω)exp{(2.1Xω)5/3[1V(Q,X,ω)ω1/3]}.

Calculations of the V(Q,X,ω) function were performed for a wide range of conditions: Q from 24 to 2+4 by factors of 2 (9 steps), X from 101 to 10+3 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 ω<0.02 (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.

2.3.

Range of Rytov Approximation Validity

Before examining the V(Q,X,ω) 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 variance

Eq. (28)

σR2=0.496Cn2k7/6L11/6=0.676(X/Q)5/3,
exceeds 0.5 (spherical wave coefficient used). However, in light of subsequent comparisons with Charnotskii’s results21,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 n11, and that perturbation field gradients be small relative to wave number k, 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,

Eq. (29)

ΩΔ1=6kΛ2/L,
where Λ models a single scatterer scale. For a turbulent outer scale of Lo, and a spectral model33 based on Kaimal’s 1968 Kansas experiment analysis,34 Λ=4.36Lo. For the turbulence parameter, Dashen used a phase variance metric also based on an outer-scale influenced covariance

Eq. (30)

ΦΔ12=k2L(2.309Cn2Lo5/3).
In this case, the RA is valid when ΦΔ1<ΩΔ1, or,

Eq. (31)

Cn2<5634Lo7/3/L3,
easily met for common outer scales of 1 m or greater.

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 ΦΔ2, corresponding to an RMS phase variance, is equivalent to one of Fried’s18 Δj2 mean aperture phase variance terms, where the order j 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,

Eq. (32)

ΦΔ22=ΔL2=(4×0.00489)(2.1X)5/3,
for an aperture of diameter D. For cutoff D, Dashen’s ΩΔ parameter becomes

Eq. (33)

ΩΔ2=12πQ2.

Setting ΦΔ2<ΩΔ2, for a given Q, a maximum X value results

Eq. (34)

Xmax(Q)=12.1[36π2Q40.00489]3/5.

This condition translates to a maximum Rytov variance function of Q

Eq. (35)

σR,max2(Q)=0.676(XmaxQ)5/3=14261Q7/3.

Typical terrestrial imaging scenarios involve Q>1/4, implying the RA remains valid up to high-turbulence conditions.

3.

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 V(Q,X,ω). This is possible due to significant correlations that exist between V(Q,X,ω) results at different Q values. This behavior is illustrated in Fig. 2 for two sets of plots at Q values of 1/16 and 16, for a series of increasing values of X.

Fig. 2

Plots of V(Q,X,ω) for Q=1/16 and Q=16 families of curves at different X values as functions of ω.

OE_53_4_044112_f002.png

The key feature of these curves is that for moderate Q values often encountered for typical terrestrial imaging scenarios, V 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 V(Q,X,0)=V(Q,X,1)=V0(Q), permitting V to be expressed

Eq. (36)

V(Q,X,ω)=V0(Q)+VΔ(Q,X,ω).

The maximum value of V (as X) must also be a function of Q, VM(Q). V0(Q) and VM(Q) are plotted in Fig. 3.

Fig. 3

Plots of V0(Q) and VM(Q) for Q=292+6.

OE_53_4_044112_f003.png

Figure 3 curves both suggest asymptotic behaviors at both extremes of Q. These behaviors can be modeled using a sigmoidal function (a variant of the hyperbolic tangent)

Eq. (37)

Σ(χ)=[exp(χ)1][exp(χ)+1]=tanh(χ/2).

To facilitate the modeling process, a rescaled and shifted sigmoid was also introduced

Eq. (38)

η=ΣS(χ,A,B,C,D)=A+BΣ[C(χD)].

The center of this function occurs at (χ,η)=(D,A) with an overall shift of 2B between asymptotes, and a central derivative of dη/dχ|χ=D=BC.

The VM(Q) function can be expressed using S as

Eq. (39)

VM(Q)=ΣS[log2(Q),0.213,0.072,1.525,0.150].

However, V0(Q) exhibits different asymptotic behaviors at large and small Q. Therefore, a third sigmoidal form was developed, featuring a central splice (piecewise sigmoidal):

Eq. (40)

ΣP(χ,A,B1,B2,C1,C2,D)={A+B1Σ[C1(χD)],χDA+B2Σ[C2(χD)],χD.
This function is continuous in its first derivative at the splice point if B1C1=B2C2. V0 is parameterized using this function via

Eq. (41)

V0(Q)=ΣP[log2(Q),0.870,0.370,0.085,0.355,1.545,1.00],
where 0.370×0.355=0.085×1.545. Equations (39) and (41) are plotted along with database derived points in Fig. 3. Variables q=log2(Q) and x=log10(X) will frequently appear hereafter as surrogates for variables Q and X in arguments to sigmoidal functions, as the majority of the modeled behaviors exhibit logarithmic dependence.

Using Eqs. (39) and (41), our problem reduces to finding an analytic expression for

Eq. (42)

V^Δ(Q,X,ω)=VΔ(Q,X,ω)/VM(Q).

Figure 4 illustrates plots of V^Δ for the same families of curves (Q=1/16 and 16) as in Fig. 2, highlighting the similarities in behaviors at different Q values. These plots also reveal how V^Δ develops increasingly complex behaviors as X increases.

Fig. 4

Plots of V^Δ=VΔ(Q,X,ω)/VM(Q) for Q=1/16 (a) and 16 (b).

OE_53_4_044112_f004.png

To model these behaviors, it was observed that function curves at high ω maintained simpler forms with increasing X 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 Q and X.

Note that the peak curves at different Q values exhibit similar X dependencies. Thus, the behaviors of the different peak curves could be modeled as shifted versions of the peak curve for Q=1, described by ωP(1,X) (abscissa) and VP(1,X) (ordinate) functions. The Q=1 abscissa curve is modeled as

Eq. (43)

ωP(1,x)=ΣP[x,0.360,0.133,0.322,+3.450,+1.424,+1.320].

Several ωP(Q,X) curves at different Q levels illustrate this shifting behavior in Fig. 5. A pair of analytic expressions that approximate this behavior are given by

Eq. (44)

ωP(Q,X)=ωP{Q=1,log10(X)xω[log2(Q)]},

Eq. (45)

xω(q)=ΣS[q,0.004,+0.091,+1.750,0.050].

Fig. 5

Plots of ωP(Q,X) for a range of Q values.

OE_53_4_044112_f005.png

Several VP ordinate function curves, designated VP(Q,X)=VΔ(Q,X,ωP) (using the non-normalized form), are illustrated for varying X in Fig. 6, for Q values ranging from 1/16 to 16. This component exhibits both shifting and scaling properties, but again was quantified relative to its behavior at Q=1, using

Eq. (46)

V1(x)=VP(1,10x)=ΣP(x,+0.122,+0.044,+0.082,+4.257,+2.300,+1.167).

Fig. 6

Plots of VP(Q,X) for Q ranging from 1/16 (bottom line) to 16 (top line), by factors of 2.

OE_53_4_044112_f006.png

The behavior at Q=1 was then rescaled using AR, and shifted using xP, in the form,

Eq. (47)

VP(Q,X)=AR[log2(Q)]V1{log10(X)xP[log2(Q)]},

Eq. (48)

AR(q)=ΣS(q,+1.051,+0.3565,+1.600,+0.150),

Eq. (49)

xP(q)=ΣS(q,+0.005,0.088,+1.550,0.150).

The function VP(Q,X) supports a further parameterization through introduction of

Eq. (50)

UP(Q,X)=VP(Q,X)/VM(Q),
using Eqs. (47) and (39), where 0<UP1. The V^Δ curve shapes appear to vary based on height UP, while ω dependence may be modeled using a Legendre-polynomial-related series expansion.

The Legendre polynomials begin as P0(z)=1, P1(z)=z, and use recurrence relation35

Eq. (51)

(n+1)Pn+1(z)=(2n+1)zPn(z)nPn1(z),
to characterize the interval 1z+1. Rescaling and normalizing these functions using

Eq. (52)

Ln(ω)=2n+1Pn(2ω1),
where z2ω1, yields an orthonormal set over 0ω1.

Using this basis set, each V^Δ curve could be separated into two functions at the peak position ω=ωP(Q,X), 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 ω=1/2, could be analyzed using Legendre-derived basis functions. Due to their symmetry, only even-order series terms were needed to characterize each function [L2n(ω)=L2n(1ω)]. The V^Δ model could, thus, be splined using the coefficient sets: C2n for ωωP, and D2n for ωωP, as

Eq. (53)

V^Δ(Q,X,ω)n=0n=NC2n[UP(Q,X)]L2n{ω/2ωP(Q,X)},ωωP;

Eq. (54)

V^Δ(Q,X,ω)n=0n=ND2n[UP(Q,X)]L2n{1[1ω]/2[1ωP(Q,X)]},ωωP.

Sets of C2n(UP) coefficient behaviors are illustrated in Fig. 7. D2n(UP) behaviors are plotted in Fig. 8. Coefficients C0 through C8, and D0, exhibit approximately hyperbolic UP dependence. The remaining coefficients exhibit approximately linear dependence.

Fig. 7

Plots of Cm(UP) expansion coefficients and associated curve fits.

OE_53_4_044112_f007.png

Fig. 8

Plots of Dm(UP) expansion coefficients and associated curve fits.

OE_53_4_044112_f008.png

Hyperbolic behaviors of Cm and Dm coefficients were expressed using constants am through em in the models

Eq. (55)

Cm(UP)=am(UPem)+bm(UPem)2+cm+dm.

Linear cases were modeled by setting bm=cm=em=0. The curves plotted in Figs. 7 and 8 represent functional fits to computed values. Fit coefficients are listed in Tables 1 and 2.

Table 1

SPGD hyperbolic Cm and Dm coefficients.

Coefficientambmcmdmem
C020.212419.63240.00310.57511.0686
C22.11182.49280.00420.35150.9808
C40.34130.33550.00510.00510.7996
C60.35290.37690.00000.01310.9041
C80.36960.39820.00000.01530.9504
D06.43725.89190.00050.57881.0341

Table 2

SPGD linear Cm and Dm coefficients.

Coefficientamdm
C100.010090.00452
C120.005500.00257
D20.302730.00057
D40.033350.00413
D60.012710.00485
D80.004250.00154
D100.004520.00186
D120.000550.00019

With the development of the Cm and Dm 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 V(Q,X,ω) 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 V 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 M0 from (20), V from (36), VΔ from (42), expressions in Eqs. (44), (47), and (50) through (55), and results of the SPGD analysis

Eq. (56)

V0(Q)=ΣP[q,0.879,0.373,0.085,0.373,1.570,0.984],

Eq. (57)

VM(Q)=ΣS[q,0.221,0.070,1.539,0.150],

Eq. (58)

ωP(1,X)=ΣP[x,0.351,0.133,0.318,+3.451,+1.407,+1.274],

Eq. (59)

log10[Xω(Q)]=ΣS[q,+0.004,0.088,+1.745,0.050],

Eq. (60)

V1(x)=VP(1,10x)=ΣP(x,+0.122,+0.044,+0.082,+4.296,+2.323,+1.179),

Eq. (61)

AR(q)=ΣS(q,+1.066,+0.352,+1.596,+0.150),

Eq. (62)

xP(q)=ΣS(q,0.005,+0.089,+1.550,0.150).

4.

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 d parameter is related to X via X=(d2/2)3/5/2.1, such that d values 0, π, 2π, 3π, and 4π correspond to X values 0.000, 1.241, 2.851, 4.638, and 6.550, respectively. Charnotskii’s N parameter is related to Q as Q=(2N/π)1/2. Reported N values of the three plots were 10.0, 1.0, and 0.1, respectively, corresponding to Q values 2.523, 0.798, and 0.252. The d=0 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 Q values that produced the closest fit to Charnotskii’s results, Q=12.0 in Fig. 9 and Q=0.02 in Fig. 11. Lines in Fig. 10 were constructed based on the reported N=1 (Q=0.798). Systematic deviations employed in these plots appear consistent with Charnotskii comments that N=10 and N=0.1 constituted limiting cases, but might be due to approximations used when Charnotskii developed Eq. (31) of Ref. 21.

Fig. 9

Comparison of SE (solid lines) and LE (dashed lines) MTF model results using Q=12.0 versus data (symbols) digitized from Ref. 21’s Fig. 8.

OE_53_4_044112_f009.png

Fig. 10

Comparison of SE (solid lines) and LE (dashed lines) MTF model results using Q=0.798 versus data (symbols) digitized from Ref. 21’s Fig. 9.

OE_53_4_044112_f010.png

Fig. 11

Comparison of SE (solid lines) and LE (dashed lines) MTF model results using Q=0.02 versus data (symbols) digitized from Ref. 21’s Fig. 10.

OE_53_4_044112_f011.png

The comparisons of Figs. 9Fig. 1011 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,

Eq. (63)

DSA[Q,X,r/D]=2(2.1X)5/3(r/D)5/3×{1V[Q,X,(r/D)](r/D)1/3},
an analogous SE structure function derived from the exponential kernel of Eq. (27), transformed into distances r in the entrance pupil. Results similar to Charnotskii’s calculations were achieved using Q values 0.02, 0.21, 0.78, and 12.0, and X=0.314 (corresponding to d=1.0). Similar Q stretching as in Figs. 9 and 11 was again observed. Curiously, the lowest line (second lowest σR2) exhibited the worst fit to Charnotskii’s results, suggesting discrepancies between the methods are not explicable simply due to RA failure.

Fig. 12

Comparison of digitized data from Fig. 2 of Ref. 26 (symbols) with computed curves (dashed) using X=0 (top line), and X=0.314 with Q values 0.02, 0.21, 0.78, and 12.00 (bottom line).

OE_53_4_044112_f012.png

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 Q stretching, a single Q adjustment appeared to correct all X cases, while no Q adjustment was used in Fig. 10 where the two strongest turbulence cases corresponded to σR2 values of 12.7 and 22.8. This ability of the RA to perform at high σR2 levels is consistent with the Sec. 2.3 analysis.

Lastly, four prior SE MTF models, (1) V=0.5 (Fried’s far-field case), (2) V=1.0 (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 V(Q,X) 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 E1=0.035179, E2=0.014089, E3=0.010788, E4=0.006431, and E5=0.000218. 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.

5.

Conclusions

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 o<D<Lo, ensuring D 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.

References

1. 

D. Sadotet al., “Restoration of thermal images distorted by atmophere, using predicted atmospheric modulation transfer function,” Infrared Phys Technol., 36 (2), 565 –576 (1995). http://dx.doi.org/10.1016/1350-4495(94)00045-M IPTEEY 1350-4495 Google Scholar

2. 

D. H. FrakesJ. W. MonacoM. J. T. Smith, “Suppression of atmospheric turbulence in video using an adaptive control grid interpolation approach,” in Proc. IEEE Int’l Conf. on Acoustics, Speech, and Sig. Proc., 1881 –1884 (2001). Google Scholar

3. 

M. A. VorontsovG. W. Carhart, “Anisoplanatic imaging through turbulent media: image recovery by local information fusion from a set of short-exposure images,” J. Opt. Soc. Am. A, 18 (6), 1312 –1324 (2001). http://dx.doi.org/10.1364/JOSAA.18.001312 JOAOD6 0740-3232 Google Scholar

4. 

D. D. BaumgartnerB. J. Schachter, “Improving FLIR ATR performance in a turbulent atmosphere with a moving platform,” Proc. SPIE, 8391 839103 (2012). http://dx.doi.org/10.1117/12.917197 PSISDG 0277-786X Google Scholar

5. 

J. GillesS. Osher, “Fried deconvolution,” Proc. SPIE, 8355 83550G (2012). http://dx.doi.org/10.1117/12.917234 PSISDG 0277-786X Google Scholar

6. 

J. Gilles, “Mao-Gilles algorithm for turbulence stabilization,” Img. Proc. On Line, 2013 173 –182 (2013). http://dx.doi.org/10.5201/ipol.2013.46 Google Scholar

7. 

J. Liuet al., “An efficient variational method for image restoration,” Abstr. Appl. Anal., 2013 213536 (2013). http://dx.doi.org/10.1155/2013/213536 Google Scholar

8. 

Y. Louet al., “Video stabilization of atmospheric turbulence distortion,” Inverse Probl Imag., 7 839 –861 (2013). http://dx.doi.org/10.3934/ipi 1930-8337 Google Scholar

9. 

D. Tofsted, “Turbulence simulation: on phase and deflector screen generation,” (2001). Google Scholar

10. 

G. PotvinJ. L. ForandD. Dion, “Some space-time statistics of the turbulent point-spread function,” J. Opt. Soc. Am. A, 24 (3), 753 –763 (2007). http://dx.doi.org/10.1364/JOSAA.24.002932 JOAOD6 0740-3232 Google Scholar

11. 

E. RepasiR. Weiss, “Analysis of image distortions by atmospheric turbulence and computer simulation of turbulence effects,” Proc. SPIE, 6941 89410S (2008). http://dx.doi.org/10.1117/12.775600 PSISDG 0277-786X Google Scholar

12. 

A. L. RobinsonJ. SmithA. Sanders, “An efficient turbulence simulation algorithm,” Proc. SPIE, 8355 83550K (2012). http://dx.doi.org/10.1117/12.919541 PSISDG 0277-786X Google Scholar

13. 

N. S. Kopeika, A System Engineering Approach to Imaging, SPIE Press, Bellingham, Washington (1998). Google Scholar

14. 

K. Krapelset al., “Atmospheric turbulence modulation transfer function for infrared target acquisition modeling,” Opt. Eng., 40 (9), 1906 –1913 (2001). http://dx.doi.org/10.1117/1.1390299 OPEGAR 0091-3286 Google Scholar

15. 

Electro-Optical Imaging System Performance, 3rd ed.SPIE Press, Bellingham, Washington (2002). Google Scholar

16. 

J. W. Goodman, Statistical Optics, J. Wiley & Sons, New York (1985). Google Scholar

17. 

M. C. RoggemannB. Welsh, Imaging Through Turbulence, CRC Press, Boca Raton, FL (1996). Google Scholar

18. 

D. L. Fried, “Statistics of a geometric representation of wavefront distortion,” J. Opt. Soc. Am., 55 (11), 1427 –1435 (1965). http://dx.doi.org/10.1364/JOSA.55.001427 JOSAAH 0030-3941 Google Scholar

19. 

D. L. Fried, “Untitled errata for Ref. [18],” J. Opt. Soc. Am., 56 (3), 410E (1966). JOSAAH 0030-3941 Google Scholar

20. 

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

21. 

M. I. Charnotskii, “Anisoplanatic short-exposure imaging in turbulence,” J. Opt. Soc. Am. A, 10 (3), 492 –501 (1993). http://dx.doi.org/10.1364/JOSAA.10.000492 JOAOD6 0740-3232 Google Scholar

22. 

A. McDonaldS. C. Cain, “Parameterized blind deconvolution of laser radar imagery using an anisoplanatic optical transfer function,” Opt. Eng., 45 (11), 116001 (2006). http://dx.doi.org/10.1117/1.2393021 OPEGAR 0091-3286 Google Scholar

23. 

J. R. DunphyJ. R. Kerr, “Turbulence effects on target illumination by laser sources: phenomenological analysis and experimental results,” Appl. Opt., 16 (5), 1345 –1358 (1977). http://dx.doi.org/10.1364/AO.16.001345 APOPAI 0003-6935 Google Scholar

24. 

D. Tofsted, “Reanalysis of turbulence effects on short-exposure passive imaging,” Opt. Eng., 50 (1), 016001 (2011). http://dx.doi.org/10.1117/1.3532999 OPEGAR 0091-3286 Google Scholar

25. 

M. I. Charnotskii, “Statistics of the point spread function for imaging through turbulence,” Proc. SPIE, 8014 80140W (2011). http://dx.doi.org/10.1117/12.884248 PSISDG 0277-786X Google Scholar

26. 

M. I. Charnotskii, “Statistical modeling of the point spread function for imaging through turbulence,” Opt. Eng., 51 (10), 101706 (2012). http://dx.doi.org/10.1117/1.OE.51.10.101706 OPEGAR 0091-3286 Google Scholar

27. 

M. I. Charnotskii, “Twelve mortal sins of the turbulence propagation science,” Proc. SPIE, 8162 816205 (2011). http://dx.doi.org/10.1117/12.894550 PSISDG 0277-786X Google Scholar

28. 

M. I. Charnotskii, “Common omissions and misconceptions of wave propagation in turbulence: discussion,” J. Opt. Soc. Am. A, 29 (5), 711 –721 (2012). http://dx.doi.org/10.1364/JOSAA.29.000711 JOAOD6 0740-3232 Google Scholar

29. 

V. I. Tatarski, Wave Propagation in a Turbulent Medium, McGraw-Hill, New York (1961). Google Scholar

30. 

D. H. Dashen, “Path integrals for waves in random media,” J. Math. Phys., 20 (5), 139 –155 (1979). http://dx.doi.org/10.1063/1.524138 JMAPAQ 0022-2488 Google Scholar

31. 

R. S. LawrenceJ. W. Strohbehn, “A survey of clear-air propagation effects relevant to optical communications,” Proc. IEEE, 58 (10), 1523 –1545 (1970). http://dx.doi.org/10.1109/PROC.1970.7977 IEEPAD 0018-9219 Google Scholar

32. 

M. I. Charnotskii, “Energy conservation: a third constraint on the turbulent point spread function,” Opt. Eng., 52 (4), 046001 (2013). http://dx.doi.org/10.1117/1.OE.52.4.046001 OPEGAR 0091-3286 Google Scholar

33. 

D. Tofstedet al., “Characterization of atmospheric turbulence during the NATO RTG-40 land field trials,” Proc. SPIE, 6551 65510J (2007). http://dx.doi.org/10.1117/12.720696 PSISDG 0277-786X Google Scholar

34. 

J. C. Kaimalet al., “Spectral characteristics of surface-layer turbulence,” Q. J. Roy. Met. Soc., 98 (417), 563 –589 (1972). http://dx.doi.org/10.1002/(ISSN)1477-870X QJRMAM 0035-9009 Google Scholar

35. 

W. Presset al., Numerical Recipes in C, Cambridge University Press, Cambridge, UK (1992). Google Scholar

36. 

M. VorontsovV. Sivokon, “Stochastic parallel gradient descent technique for high-resolution wave-front phase-distortion correction,” J. Opt. Soc. Am. A, 15 (10), 2745 –2758 (1998). http://dx.doi.org/10.1364/JOSAA.15.002745 JOAOD6 0740-3232 Google Scholar

37. 

D. Tofsted, “Short-exposure passive imaging through path-varying convective boundary layer turbulence,” Proc. SPIE, 8355 83550J (2012). http://dx.doi.org/10.1117/12.920840 PSISDG 0277-786X Google Scholar

Biography

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.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
David H. Tofsted "Extended high-angular-frequency analysis of turbulence effects on short-exposure imaging," Optical Engineering 53(4), 044112 (25 April 2014). https://doi.org/10.1117/1.OE.53.4.044112
Published: 25 April 2014
Lens.org Logo
CITATIONS
Cited by 5 scholarly publications and 3 patents.
Advertisement
Advertisement
KEYWORDS
Modulation transfer functions

Turbulence

Databases

Data modeling

Optical engineering

Diffraction

Imaging systems

RELATED CONTENT


Back to Top