Estimation of turbulence from time-lapse imagery

Jack E. McCrae; Santasri R. Bose-Pillai; Steven T. Fiorino

doi:10.1117/1.OE.56.7.071504

7 February 2017 Estimation of turbulence from time-lapse imagery

Jack E. McCrae, Santasri R. Bose-Pillai, Steven T. Fiorino

Author Affiliations +

Optical Engineering, Vol. 56, Issue 7, 071504 (February 2017). https://doi.org/10.1117/1.OE.56.7.071504

Abstract

Atmospheric turbulence parameters are estimated for an imaging path based on time-lapse imaging results. Atmospheric turbulence causes frame-to-frame shifts of the entire image as well as parts of the image. The statistics of these shifts encode information about the turbulence strength (as characterized by Cn2, the refractive index structure function constant) along the optical path. The shift variance observed is simply proportional to the variance of the tilt of the optical field averaged over the area being tracked and averaged over the camera aperture. By presuming this turbulence follows the Kolmogorov spectrum, weighting functions, which relate the turbulence strength along the path to the shifts measured, are derived. These weighting functions peak at the camera and fall to zero at the object. The larger the area observed, the more quickly the weighting function decays. One parameter we would like to estimate is r0 (the Fried parameter or atmospheric coherence diameter.) The weighting functions derived for pixel sized or larger parts of the image all fall faster than the weighting function appropriate for estimating the spherical wave r0. If we were to presume that Cn2 is constant along the path, then an estimate for r0 could be obtained for each area tracked, but since the weighting function for r0 differs substantially from that for every realizable tracked area, it can be expected that this approach would yield a poor estimate. Instead, the weighting functions for a number of different patch sizes can be combined through the Moore–Penrose pseudoinverse to create a weighting function that yields the least-squares optimal linear combination of measurements for the estimation of r0. This approach is carried out for one example and is shown to give noisy results. A modified version of this approach that creates larger patches by averaging several smaller patches together solves this noise issue. This approach can also work to estimate other atmospheric para

1. Introduction

A wide variety of techniques exist to characterize atmospheric turbulence and its effects on optical systems. These techniques include scintillometry; scintillation, detection, and ranging;¹^,² differential image motion monitoring;³ and difference in differential tilt variance.⁴ Recently reported techniques have used weather radar⁵^,⁶ and satellite data.⁷ This wide range of approaches indicates strong interest in this area and high demand for the information provided. This information is critical to predicting the performance of optical systems in real-world environments, better understanding the turbulent atmosphere, and designing systems to work through turbulence.

This paper discusses the use of time-lapse imagery to make estimations of turbulence; this approach offers advantages and disadvantages with respect to other methods. Some of the advantages of using time-lapse imaging include its simplicity, potential low cost, passive nature, need to access only one end of the optical path, and ability to leverage large amounts of data, which can be collected quickly and easily with digital cameras. Disadvantages include the need for a target object at the far side of the path, the potentially large amount of data processing required, and the heavy influence of turbulence near the camera on the observed image motion. Techniques are developed in this paper to mitigate this final point. Important results from the effort reported here include an expression for the expected aperture averaged wavefront tilt variance due to turbulence from an extended source, a technique to allow parameters based on arbitrary path weightings to be estimated from time-lapse imagery, and an approach to combat noise by composing the wavefront tilt variance for large patches on an extended object from smaller ones. In particular, for this last item, this paper further extends the results of our prior work.⁸^,⁹

The remainder of this paper consists of Secs. 2–4. Section 2 describes the experimental configuration, explains the correlation process used to track image motion, derives the weighting functions that related the turbulence along the path to the image shift variances expected, and details how a linear combination of multiple weighting functions will be used to achieve some desired weighting function. Section 3 starts with an example of how the weighting function for the spherical wave Fried parameter, $r_{0}$ , can be approximated with a combination of the weighting functions for a number of tracked patch sizes. This technique is then used to estimate $r_{0}$ from the tilt variances measured using just two different patches. This result is shown to be noise prone, and an improved technique that builds up the tilt variance of large patches from smaller patches is demonstrated. Section 4 of this paper is titled “Conclusions,” where the results are reviewed, there is further analysis of these results, and a discussion about further work in this area is given.

2. Methodology

2.1.

Time-Lapse Imaging Experiment

Images of the Good Samaritan Hospital were captured with a Canon 40D camera ( $5.7 - μ m$ pixel pitch) using a 300-mm telephoto lens ( $f / #$ 5.6) from a ground floor laboratory at the Air Force Institute of Technology (AFIT). The distance from the AFIT lab to the hospital is 12.8 km. The camera was mounted on a sturdy tripod, pointed out through a closed window, and located a few feet behind the window glass. Window blinds and a cardboard shade were configured to minimize any effect from solar heating on the camera and tripod. Images were captured at a rate of one per minute. It is presumed that wavefront tilts due to atmospheric turbulence are uncorrelated at this timescale. The path from AFIT to the hospital is close to horizontal and the bulk of the path is roughly 60 m above the floor of the Miami Valley. This is shown in Fig. 1; the imaging path is a little closer to horizontal than depicted since the camera is pointed at the upper floors of the hospital, and the lowest couple of floors are obscured by the ground. Frame-to-frame shifts of selected portions of the image are taken as the change in the wavefront tilt due to (i.e., averaged over) that portion of the image and averaged over the camera aperture. Figure 2 shows a $256 \times 256$ image of the target hospital cut from the full image. This image portion is about 60 m on a side at the hospital, and each pixel thus spans about 0.24 m at the hospital.

Fig. 1

Experimental path and elevation profile.

Fig. 2

Image of Good Samaritan Hospital collected at 6 a.m. on July 25, 2014.

2.2.

Correlation Technique

Image shifts were evaluated using a subpixel precision correlation technique. First, a desired patch on the hospital was selected, and its coordinates were identified in a representative image from the set. A Gaussian intensity taper centered on this patch was applied to all images in the set. Since the patches used were large, the patch position was not corrected over the course of the day to account for image motion; several pixels of motion attributed mostly to changes in the refractive index gradient were present. The cross correlation was then computed between adjacent pairs of images in the time series using the fast Fourier transform. The location of the peak value of the cross correlation was identified, and an $11 \times 11$ sample region centered on the peak was excised from the cross-correlation matrix. A least-squares fit was used to match a paraboloid to this subregion, and the position of the peak of this paraboloid was taken as the image shift. More precisely, the image shift was computed as the sum of the peak location shift in the correlation and the paraboloid peak location position. The peak location shift in the correlation was an integral number of pixels, usually zero, whereas the paraboloid peak location gave a fractional correction. Across the 1000 or so images collected in a day, this technique gave results that showed little drift when the shifts were summed to give a net motion. Other techniques were explored, including those using phase correlation, but these provided results that appeared noisier.

2.3.

Path Weighting Functions

The expected variance of the wavefront tilt was evaluated using an approach outlined by Fried¹⁰ and Winick and Marquis¹¹ but modified and extended for the case under study here. More precisely, the referenced technique to evaluate the expected tilt difference variance between two sources averaged across a receiver aperture was modified for spherical rather than plane waves, then turned around and applied at the source plane to generate the expected tilt variance as averaged over both the aperture and a Gaussian source patch. This derivation is given as follows.

The $z$ -tilt over an aperture of diameter $D$ , when viewing a source in the direction $θ$ can be expressed as Ref. 12

Eq. (1)

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

where

λ

is the wavelength,

φ (r, θ)

is the turbulence-induced wavefront distortion at aperture coordinate

r

, and

Eq. (2)

W (r) = {\begin{cases} 1 & | r | \leq 0.5 D \\ 0 & | r | > 0.5 D \end{cases} .

The mean correlation between tilts observed at the aperture due to two sources at viewing directions $θ_{1}$ and $θ_{2}$ can be written as follows:

Eq. (3)

⟨ α (θ_{1}) \cdot α (θ_{2}) ⟩ = {(\frac{32 λ}{π^{2} D^{4}})}^{2} ⟨ \iint d r d r^{'} W (r) W (r^{'}) r \cdot r^{'} φ (r, θ_{1}) φ (r^{'}, θ_{2}) ⟩,

where the angular brackets indicate ensemble averaging. Interchanging the order of integration and ensemble averaging yields

Eq. (4)

⟨ α (θ_{1}) \cdot α (θ_{2}) ⟩ = {(\frac{32 λ}{π^{2} D^{4}})}^{2} \iint d r d r^{'} W (r) W (r^{'}) r \cdot r^{'} ⟨ φ (r, θ_{1}) φ (r^{'}, θ_{2}) ⟩ .

Because the integrand is an odd function of both

r

and

r^{'}

,

\iint d r d r^{'} W (r) W (r^{'}) r \cdot r^{'} = 0

, terms that are functions of either

r

or

r^{'}

, but not both, can be added without changing the result of the integration.

Hence, Eq. (4) can be expressed as

Eq. (5)

⟨ α (θ_{1}) \cdot α (θ_{2}) ⟩ = (- \frac{1}{2}) {(\frac{32 λ}{π^{2} D^{4}})}^{2} \iint d r d r^{'} W (r) W (r^{'}) r \cdot r^{'} D (r - r^{'}, θ_{1} - θ_{2}),

where

D (r - r^{'}, θ_{1} - θ_{2}) = ⟨ {[φ (r, θ_{1}) - φ (r^{'}, θ_{2})]}^{2} ⟩

is the phase structure function.

For a spherical wave propagating through turbulence characterized by the Kolmogorov power spectrum, the phase structure function can be written as in Ref. 13

Eq. (6)

D (r - r^{'}, θ_{1} - θ_{2}) = 2.91 k^{2} \int_{0}^{L} d z C_{n}^{2} (z) {| (r - r^{'}) (1 - \frac{z}{L}) + z (θ_{1} - θ_{2}) |}^{5 / 3},

where

k = 2 π / λ

is the wave number and

L

is the path length. The receiver is located at

z = 0

, and the source is located at

z = L

. Substituting Eq. (6) in Eq. (5),

Eq. (7)

⟨ α (θ_{1}) \cdot α (θ_{2}) ⟩ = (- \frac{2.91}{2}) {(\frac{64}{π D^{4}})}^{2} \int_{0}^{L} d z C_{n}^{2} (z) \times \iint d r d r^{'} W (r) W (r^{'}) r \cdot r^{'} {| (r - r^{'}) (1 - \frac{z}{L}) + z (θ_{1} - θ_{2}) |}^{5 / 3} .

The integrations over

r

and

r^{'}

can be done using techniques outlined by Fried¹⁰ and Winick and Marquis.¹¹ Applying those techniques, Eq. (7) reduces to

Eq. (8)

⟨ α (θ_{1}) \cdot α (θ_{2}) ⟩ = (- \frac{2.91}{2}) {(\frac{16}{π})}^{2} D^{- 1 / 3} \int_{0}^{L} d z C_{n}^{2} (z) \int_{0}^{2 π} d ϑ \int_{0}^{1} d u [(u \cos^{- 1} u) - u^{2} (3 - 2 u^{2}) \sqrt{1 - u^{2}}] \times {[u^{2} {(1 - \frac{z}{L})}^{2} + {(\frac{z}{D} | θ_{1} - θ_{2} |)}^{2} + 2 u (1 - \frac{z}{L}) (\frac{z}{D} | θ_{1} - θ_{2} |) \cos ϑ]}^{5 / 6} .

Each pixel in the time-lapse imagery corresponds to a patch of $\sim 0.24 m$ , not just a point on the target. Hence, the shift (or the tilt) measured from the whole image or even a pixel on the image is not a tilt due to a single point source, but rather an average tilt due to several incoherent point sources over a patch. Since a Gaussian window is applied to the images before the cross correlation is computed, to make the analysis consistent, the source patch is modeled with the same Gaussian. The patch-averaged tilt $α_{p}$ is defined as

Eq. (9)

α_{p} = \frac{\int d θ α (θ) P_{G} (θ)}{\int d θ P_{G} (θ)},

where

P_{G} (θ) = e^{- \frac{4 θ^{2} L^{2}}{d^{2}}}

, with

d

as the

1 / e

patch diameter and

θ = | θ |

.

Hence, the patch-averaged tilt variance is

Eq. (10)

⟨ α_{p}^{2} ⟩ = \frac{\iint d θ_{1} d θ_{2} ⟨ α (θ_{1}) \cdot α (θ_{2}) ⟩ P_{G} (θ_{1}) P_{G} (θ_{2})}{\iint d θ_{1} d θ_{2} P_{G} (θ_{1}) P_{G} (θ_{2})} .

The order of integration and ensemble averaging is interchanged to obtain the above expression. The term with angle brackets in Eq. (10) is just the covariance of tilts due to two point sources. Substituting Eq. (8) in Eq. (10),

Eq. (11)

⟨ α_{p}^{2} ⟩ = (- \frac{2.91}{2}) {(\frac{16}{π A_{p}})}^{2} D^{- 1 / 3} \int_{0}^{L} d z C_{n}^{2} (z) \iint d θ_{1} d θ_{2} P_{G} (θ_{1}) P_{G} (θ_{2}) \int_{0}^{2 π} d ϑ \int_{0}^{1} d u [(u \cos^{- 1} u) - u^{2} (3 - 2 u^{2}) \sqrt{1 - u^{2}}] \times {[u^{2} {(1 - \frac{z}{L})}^{2} + {(\frac{z}{D} | θ_{1} - θ_{2} |)}^{2} + 2 u (1 - \frac{z}{L}) (\frac{z}{D} | θ_{1} - θ_{2} |) \cos ϑ]}^{5 / 6},

where

A_{p} = \int d θ P_{G} (θ) = \frac{π d^{2}}{4 L^{2}}

.

Let

Eq. (12)

θ_{1} - θ_{2} = \frac{d}{L} x θ_{1} + θ_{2} = \frac{2 d}{L} y .

Changing the variables of integration, Eq. (11) reduces to

Eq. (13)

⟨ α_{p}^{2} ⟩ = (- \frac{2.91}{2 π}) {(\frac{16}{π})}^{3} D^{- 1 / 3} \int_{0}^{L} d z C_{n}^{2} (z) \iint d x d y [e^{- {(x + 2 y)}^{2}} e^{- {(2 y - x)}^{2}}] \times \int_{0}^{2 π} d ϑ \int_{0}^{1} d u [(u \cos^{- 1} u) - u^{2} (3 - 2 u^{2}) \sqrt{1 - u^{2}}] \times {[u^{2} {(1 - \frac{z}{L})}^{2} + {(\frac{z d}{D L} | x |)}^{2} + 2 u (1 - \frac{z}{L}) (\frac{z d}{D L} | x |) \cos ϑ]}^{5 / 6} .

Now, $\int d y e^{- {(x + 2 y)}^{2}} e^{- {(2 y - x)}^{2}} = \int d y e^{- 8 (y^{2} + \frac{x^{2}}{4})} = \int_{0}^{\infty} d y \int_{0}^{2 π} d β [y e^{- 8 (y^{2} + \frac{x^{2}}{4})}] = \frac{π}{8} e^{- 2 x^{2}}$ , where $x = | x |$ and $y = | y |$ . Substituting the above result in Eq. (13), the expression for patch-averaged tilt variance becomes

Eq. (14)

⟨ α_{p}^{2} ⟩ = \int_{0}^{L} d z C_{n}^{2} (z) f_{α} (z),

where

Eq. (15)

f_{α} (z) = - 5.82 {(\frac{16}{π})}^{2} D^{- 1 / 3} \int_{0}^{\infty} d x (x e^{- 2 x^{2}}) \int_{0}^{2 π} d ϑ \int_{0}^{1} d u [(u \cos^{- 1} u) - u^{2} (3 - 2 u^{2}) \sqrt{1 - u^{2}}] \times {[u^{2} {(1 - \frac{z}{L})}^{2} + {(\frac{z d}{D L} x)}^{2} + 2 u (1 - \frac{z}{L}) (\frac{z d}{D L} x) \cos ϑ]}^{5 / 6},

is the path weighting function. The weighting function for a point source is recovered by letting

d \to 0

in Eq. (15). The weighting function is seen to depend only on geometry, with the

5 / 6

th power in the integrand hinting at the presumption that the turbulence obeys the Kolmogorov power spectrum locally.

No simpler form for the weighting function could be obtained; hence in the present work, Eq. (15) was evaluated numerically for cases of interest. Figure 3 shows a plot of $f_{α} (z)$ for several different patch sizes on the target.

Fig. 3

Weighting functions for computing the expected tilt variance for differently sized patches of an image from the turbulence along the path.

The aperture size and the path length used in the evaluation were the same as those used in the experiment, i.e., $D = 53.6 mm$ and $L = 12.8 km$ . The weight is maximum at the receiver (camera) ( $z = 0$ ) and drops to zero at the target ( $z = L$ ). This implies that the turbulence near the target hardly affects the image shifts seen by the camera. It is also apparent from Fig. 3 that the larger the size of the patch being tracked, the smaller the variance of the motion. This is unsurprising as each subarea within a patch can have larger, random motion, but these random motions often average to a considerably smaller net motion for the patch. The tilt variance expected for even a single pixel ( $d = 0.24 m$ ) in the image is significantly less than that for a point source.

2.4.

Creating an Arbitrary Weighting Function

If we were to presume that the turbulence is constant along the path, then the tilt variance associated with each image patch would provide an estimate for $C_{n}^{2}$ . If we do not make this presumption, then a set of tilt variances from differently sized patches can be seen as encoding differences in $C_{n}^{2}$ along the path. Members of a set of path weighting functions, from a variety of differently sized image patches, each characterize the turbulence along (almost) the same path, but each weights that path differently. We wish to determine what linear combination of weighting functions from this set will best reproduce the path weighting function associated with some atmospheric parameter we desire to estimate. For example, the spherical wave atmospheric coherence diameter or $r_{0, sw}$ is given by Refs. 13 and 14

Eq. (16)

r_{0, sw} = {[0.423 k^{2} \int_{0}^{L} C_{n}^{2} {(1 - \frac{z}{L})}^{5 / 3}]}^{- 3 / 5} .

This equation has been written so the integration takes place from the camera to the object, so the camera is at zero and the hospital is at $L$ as above. In this expression, $C_{n}^{2}$ is weighted by $z^{5 / 3}$ along the path. This is the same weighting as Eq. (15) for $f_{α} (z)$ in the case where the patch size, $d$ , goes to zero. As there is no beacon or point source at the target, no part of the image corresponds to a weighting function of this form. In fact, the weighting function for a single pixel sized patch has nearly 25% less area than this desired weighting function. To determine the appropriate linear combination of weighting functions, the Moore–Penrose pseudoinverse can be used to find the least-squares optimal way to achieve the desired weighting function in terms of the basis set available. By sampling the weighting functions along the path and generating a matrix, $M$ , where the rows are formed from the individual sampled weighting functions, and the columns then correspond to the range where these functions are sampled, the desired weights can be computed as

Eq. (17)

W = M^{+} F,

where

W

is the weight to be applied to each weighting function, and

F

is the desired weighting function sampled the same way as

M

, and

M^{+}

is the pseudoinverse. The success of this operation can then be revealed by examining

M W

, which is the actual weighting function generated by attempting to match

F

with a linear combination of functions drawn from

M

. When the rows of

M

are linearly independent, the pseudoinverse may be computed as

Eq. (18)

M^{+} = {(M^{T} M)}^{- 1} M^{T} .

The Moore–Penrose pseudoinverse instead uses a singular value decomposition to avoid inverting small singular values and still gives useful results when

M^{T} M

is singular or nearly so. For the results presented next, both pseudoinversion techniques gave nearly the same results. With further experimentation, the Moore–Penrose pseudoinverse gave better results for cases where very similar weighting functions were included in

M

, but these cases are not included here.

3. Results

3.1.

Approximated Weighting Function Example for Atmospheric Coherence Diameter

Figure 4 shows the path weighting functions for patch diameters of 1, 2, 3, 4, and 5 m along with the desired weighting function with a $z^{5 / 3}$ dependence plus the weighting function achieved with just these five patch sizes. These functions are sampled every 100 m along the 12.8 km path.

Fig. 4

Weighting function achieved when trying to match a desired function with a linear combination of weighting functions.

Although the match of the derived weighting function to the desired is far from perfect, it is reasonable given the small number of degrees of freedom available, especially considering that the basis set used is neither very variegated nor particularly close to what is desired. A bigger concern is the large values of the coefficients applied to the weighting functions to achieve this result. The values found for this example are 22.64, $- 90.66$ , 124.42, $- 58.01$ , and 2.58 for patch sizes of 1 to 5 m, respectively. The weighting functions are being multiplied by large numbers and subtracted. This will likely amplify the noise present substantially when this approach is used to estimate parameters numerically.

3.2.

Example Estimation of $r_{0}$

To carry through an example of this process with actual data, a somewhat simpler case will be used. Just two patch sizes will be employed, 10 and 20 m. Figure 5 shows the weighting functions for these patch sizes as well as the desired function and the least squares match, as in Fig. 4.

Fig. 5

Weighting function achieved for the best match of the desired function using just two patch sizes (10 and 20 m).

If only the 10-m patch were used, the best match would be with a multiplication by 3.8, and the overshoot for small distances would be about twice as high. The weights determined were 10.13 and $- 9.19$ for the 10- and 20-m patch sizes, respectively. Figure 6 shows the frame-by-frame shifts in pixels extracted from these patch sizes for patches centered on the hospital.

Fig. 6

Shifts computed from image data collected on July 25, 2014. The vertical scale is marked in pixels; the horizontal units are minutes since 5:50 a.m. (a) and (b) For the 10-m patch size and (c) and (d) for the 20-m patch size. (a) and (c) For the vertical shifts and (b) and (d) for the horizontal shifts.

From these shifts, $C_{n}^{2}$ can be computed under the presumption that turbulence is constant along the path. Figure 7 shows this result as computed from the shifts in Fig. 6. Likewise, Eq. (16) can be used to compute $r_{0}$ under this same presumption of constant turbulence along the path. This result, for just the 20-m patch size, is shown in Fig. 8. Note that it looks a lot like the curve for $C_{n}^{2}$ but flipped upside down. Since turbulence can be expected to be higher near the ground, only near the camera is the path near the ground (excepting near the hospital where the weights are all very small), and this approach overweights turbulence near the camera, it can be expected that this will underestimate $r_{0}$ . The technique we wish to use is capable, at least in principle, of overcoming this limitation.

Fig. 7

$C_{n}^{2}$ from image shifts presuming constant turbulence along the path.

Fig. 8

Estimated $r_{0}$ from image shifts presuming constant turbulence along the path.

Figure 9 shows the result of using a linear combination of the measured tilt variance with the weights determined above to estimate $r_{0}$ . This estimate of $r_{0}$ was computed using Eq. (19), which takes the measured tilt variances ( $α_{p 10}^{2}$ and $α_{p 20}^{2}$ in a somewhat shorthand notation) and the weights to generate the result. As before, $D$ is the telescope diameter (53 mm in this case). The wavelength of light, $λ$ , is taken as 500 nm for this experiment.

Eq. (19)

r_{0, sw} = {[\frac{D^{1 / 3}}{0.36 λ^{2}} (10.13 ⟨ α_{p 10}^{2} ⟩ - 9.19 ⟨ α_{p 20}^{2} ⟩)]}^{- 3 / 5} .

When the expression inside the square brackets is negative, this approach gives no estimate for

r_{0}

. This happens quite a few times in Fig. 9. In particular, in the morning, the tilt variance measured for the 20-m patch was larger than that of the 10-m patch, and this results in no estimate. This issue can be seen by comparing Figs. 6(a) with 6(c) and 6(b) with 6(d). Equation (19) was used on the horizontal and vertical components separately, and the variances were computed from the shifts between frames. These details contribute competing factors of 2, which cancel here. The variances were computed from a 1-h long (i.e., 60 frame) boxcar moving window.

Fig. 9

Estimated $r_{0}$ from a linear combination of weighting functions using a single 10-m patch and a single 20-m patch.

There is no distribution of $C_{n}^{2}$ along a path that can account for observing a larger tilt variance on larger patches than smaller ones. However, this issue can be explained by a number of things. For one, the paths for the large patch and the small patch are not identical, since they necessarily sample different air along the way because they have different sizes. It could be that a large tilt near the edge of the large patch is not captured by the smaller patch, which was centered on the large patch in this case. Something like this must certainly happen occasionally. Computing the variance over a larger number of frames might help this issue, but the 60 point moving average used already means an hour passes between independent estimates of $r_{0}$ . The substantial amount of noise present, plus the similarity of the path weighting functions for the two patch sizes used, contributes to the high variability seen in the results in Fig. 9. Some of the coarse features appear to be reasonable and in general agreement with the results in Fig. 8 (which should be somewhat incorrect themselves, but not so noisy). Notably, there is floor in $r_{0}$ around 1 cm over much of the day, with a climb at the end of the day as the turbulence settles down. The horizontal and vertical data were kept separate as a gauge on the noise present plus as an indicator of experimental problems. The horizontal and vertical data sets were expected to give the same results.

The patch sizes were chosen somewhat arbitrarily; large patches are easier to work with since the image shifts will not move a significant proportion of the image out from under the patch window. This resulted in the two path weighting functions being perilously similar to each other and quite dissimilar from the desired function. It was observed that the lowest weights were achieved using only two patches: the smallest available patch and a large one.

3.3.

Modified Approach for Estimation of $r_{0}$

To better deal with the noise seen in the previous results, a different approach was taken. Using a number of small patches that sum together to give a large patch, the patches use what is more nearly an identical path through the air; plus, it is mathematically guaranteed that the small patches will have an average tilt variance greater than or equal to that of the larger patch. This will be guaranteed since the small patch tilt variance will be computed as the average of the tilt variances over the small patches, while the large patch tilt variance will be computed as the variance of the average tilt over the small patches.

The same 10- and 20-m patch sizes will be used as above, but now, the 20-m patch will be composed of seven 10-m patches arranged in a hexagonal pattern. A nonlinear minimization routine was used to determine the heights and positions of seven 10-m Gaussian patches that best fit a single 20-m Gaussian patch. There were only three parameters to be fit since a centered hexagonal pattern was presumed: the height of the central 10-m Gaussian plus the height and position of the six Gaussians surrounding the center. The optimal solution was found to be with a central Gaussian with 94% of the height of the target Gaussian, surrounded by six Gaussians with 45% of the height of the target Gaussian and located 84% of the way from the center to the target Gaussian’s $1 / e$ point. Figure 10 shows a surface plot of the 20-m Gaussian function as well as the sum of the seven smaller Gaussians. The fit appears to be quite good.

Fig. 10

A 20-m Gaussian patch along with seven 10-m Gaussian patches fit to this shape.

Equation 19 will be used again to estimate $r_{0}$ for this path; this time, however, the image shifts from the seven different 10-m patches will be used in two different ways to compute $α_{p 10}^{2}$ and $α_{p 20}^{2}$ , the desired variances. The small patch size variance, $α_{p 10}^{2}$ , will be computed as the weighted average of the tilt variances computed from these seven patches projected onto the hospital images using the previously computed heights as the weights. The tilt variance for the larger patch, $α_{p 20}^{2}$ , will be computed as the variance of the identically weighted average of the image shifts. This procedure gave substantially better results, at least in that the value within the square brackets of Eq. (19) was always positive. The results from this procedure are shown in Fig. 11.

Fig. 11

Estimated $r_{0}$ from a linear combination of weighting functions, using seven 10-m image patches and a composition of these seven patches to form a 20-m patch.

The results are now seen to be very comparable to those from Fig. 8, perhaps even a little bit quieter; this seems more likely to be due to the larger number of patches used rather than anything else. Although the result might also show somewhat higher values of $r_{0}$ than those seen in Fig. 8, they are in very good agreement. This implies that the overly heavyweight applied to the first part of the path when only a single patch size is used did not make a great deal of difference for this data set. The disagreement between the vertical and the horizontal results at the start of the data, and again somewhat at the end, is presently unexplained. Although there is some image motion due to refractive effects, especially in the morning, it does not appear to be to blame. To rule this out, the tilt variances were recomputed using second nearest neighbor images, and no significant changes were observed.

4. Conclusions

In conclusion, time-lapse imagery has been used to estimate turbulence along the imaging path. Weighting functions were derived relating the turbulence along the path to the measured image shifts as a function of the experimental parameters and the size of the image patch tracked. Approaches to estimate atmospheric parameters, such as $r_{0}$ , from the imagery were presented and discussed.

Some consideration was given to attempting to profile turbulence along the path using the techniques here. In the case presented here, the quick drop off of the weighting functions with motion away from the camera and their low values deep in the path make this seems to be difficult. Instead, it seemed better to attempt to estimate desired atmospheric parameters more directly. The technique presented here could easily be used for other parameters to allow time-lapse measurements to mimic the results of a scintillometer or other instrument. Presently, the authors are working to deploy this time-lapse system along a path in conjunction with a scintillometer to allow these techniques to be directly compared.

There is substantial room to take this approach further still. A larger and richer set of weighting functions would allow desired functions to be approximated more closely, and hopefully with smaller coefficients; differential motion between image patches appears to be one such set. A larger aperture would clearly help this technique work better, but the idea here was to extract maximum value from readily available data. There is a presumption underpinning this approach that the Gaussian window is well supported by the image data (i.e., there are features sufficiently densely placed within each Gaussian window for the correlation to give results representative of the window). Although this seems reasonable for large window sizes, it might not be a good presumption for smaller ones. This idea, and others, will be investigated with synthetic image data.

Acknowledgments

This work was supported by the Air Force Office of Scientific Research (AFOSR) through a Multidisciplinary Research Program of the University Research Initiative (MURI) Grant and by the High Energy Laser Joint Technology Office (HEL-JTO). The views expressed in this paper are those of the authors and do not necessarily reflect the official policy or position of the Air Force, the Department of Defense, or the U.S. Government.

References

1.

R. Avila, J. Vernin and E. Masciadri, “Whole atmospheric-turbulence profiling with generalized SCIDAR,” Appl. Opt., 36 (30), 7898 –7905 (1997). http://dx.doi.org/10.1364/AO.36.007898 APOPAI 0003-6935 Google Scholar

2.

V. Kluckers et al., “Profiling of atmospheric turbulence strength and velocity using a generalized SCIDAR technique,” A&A, 130 (1), 141 –155 (1998). http://dx.doi.org/10.1051/aas:1998217 AAESB9 0365-0138 Google Scholar

3.

M. Sarazin and F. Roddier, “The ESO differential image motion monitor,” A&A, 227 (1), 294 –300 (1990). Google Scholar

4.

M. Whiteley, D. Washburn and L. Wright, “Differential- tilt technique for saturation- resistant profiling of atmospheric turbulence,” Proc. SPIE, 4494 221 –232 (2002). http://dx.doi.org/10.1117/12.454795 PSISDG 0277-786X Google Scholar

5.

S. Fiorino, “Satellite and radar measurement of CT2, Cn2, and Cv2,” in Imaging and Applied Optics 2014, OSA Technical Digest (online), (2014). Google Scholar

6.

L. Burchett and S. Fiorino, “Measurement of Cn2 profiles from weather radar data and other microwave signals and conversion to visible and NIR Cn2 profiles,” in Imaging and Applied Optics 2013, OSA Technical Digest (online), (2013). Google Scholar

7.

D. Meier and S. Fiorino, “Correlated satellite-derived turbulence, clouds and aerosol data,” in Imaging and Applied Optics 2014, OSA Technical Digest (online), (2014). Google Scholar

8.

S. Basu, J. E. McCrae and S. T. Fiorino, “Estimation of the path-averaged atmospheric refractive index structure constant from time-lapse imagery,” Proc. SPIE, 9465 94650T (2015). http://dx.doi.org/10.1117/12.2177330 PSISDG 0277-786X Google Scholar

9.

J. E. McCrae, S. Basu and S. T. Fiorino, “Estimation of atmospheric parameters from time-lapse imagery,” Proc. SPIE, 9833 983303 (2016). http://dx.doi.org/10.1117/12.2223986 PSISDG 0277-786X Google Scholar

10.

D. Fried, “Varieties of isoplanatism,” Proc. SPIE, 75 20 –29 (1976). http://dx.doi.org/10.1117/12.954733 PSISDG 0277-786X Google Scholar

11.

K. Winick and D. Marquis, “Stellar scintillation technique for the measurement of tilt anisoplanatism,” J. Opt. Soc. Am. A, 5 (11), 1929 –1936 (1988). http://dx.doi.org/10.1364/JOSAA.5.001929 JOAOD6 0740-3232 Google Scholar

12.

D. Fried, “Differential angle of arrival: theory, evaluation, and measurement feasibility,” Radio Sci., 10 (1), 71 –76 (1975). http://dx.doi.org/10.1029/RS010i001p00071 Google Scholar

13.

M. C. Roggemann, B. M. Welsh and B.R. M. Hunt, Imaging Through Turbulence, CRC Press, Boca Raton, Florida (1996). Google Scholar

14.

J. D. Schmidt, Numerical Simulation of Optical Wave Propagation, SPIE Press, Bellingham, Washington (2010). Google Scholar

Biography

Jack E. McCrae received his BS degree in physics from Massachusetts Institute of Technology in 1984, his MS degree in physics (Optics) from Air Force Institute of Technology in 1993, and his PhD in physics from Air Force Institute of Technology in 1997. He is a retired Air Force Colonel with 27 years of service and currently with the Center for Directed Energy as a research assistant professor in the Engineering Physics Department at AFIT.

Santasri R. Bose-Pillai received her PhD in electrical engineering with a focus in optics in 2008 and her MSEE degree in 2005, both from New Mexico State University. She received her BSEE degree with honors from Jadavpur University, India, in 2000. Currently, she is a senior research associate at AFIT's Center for Directed Energy within the Engineering Physics Department. She is a member of OSA, SPIE, and DEPS.

Steven T. Fiorino received his BS degree in geography and meteorology from Ohio State in 1987 and Florida State University in 1989. He additionally holds MS degree in atmospheric dynamics from Ohio State in 1993, and his PhD in physical meteorology from Florida State in 2002. He is a retired Air Force Lieutenant Colonel and currently an associate professor of atmospheric physics within the Engineering Physics Department at AFIT and he is the director of the Center for Directed Energy.

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.

Citation Download Citation

Jack E. McCrae, Santasri R. Bose-Pillai, and Steven T. Fiorino "Estimation of turbulence from time-lapse imagery," Optical Engineering 56(7), 071504 (7 February 2017). https://doi.org/10.1117/1.OE.56.7.071504

Received: 3 November 2016; Accepted: 13 January 2017; Published: 7 February 2017

Access the abstract

JOURNAL ARTICLE
9 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 19 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Turbulence

Cameras

Optical engineering

Atmospheric turbulence

Image analysis

Spherical lenses

Wavefronts

1.

Introduction

2.

Methodology

2.1.

Time-Lapse Imaging Experiment

Fig. 1

Fig. 2

2.2.

Correlation Technique

2.3.

Path Weighting Functions

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

Eq. (15)

Fig. 3

2.4.

Creating an Arbitrary Weighting Function

Eq. (16)

Eq. (17)

Eq. (18)

3.

Results

3.1.

Approximated Weighting Function Example for Atmospheric Coherence Diameter

Fig. 4

3.2.

Example Estimation of r0

Fig. 5

Fig. 6

Fig. 7

Fig. 8

Eq. (19)

Fig. 9

3.3.

Modified Approach for Estimation of r0

Fig. 10

Fig. 11

4.

Conclusions

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years

Example Estimation of $r_{0}$

Modified Approach for Estimation of $r_{0}$