A wide variety of techniques exist to characterize atmospheric turbulence and its effects on optical systems. These techniques include scintillometry; scintillation, detection, and ranging;1,2 differential image motion monitoring;3 and difference in differential tilt variance.4 Recently reported techniques have used weather radar5,6 and satellite data.7 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.8,9
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, , can be approximated with a combination of the weighting functions for a number of tracked patch sizes. This technique is then used to estimate 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.
Time-Lapse Imaging Experiment
Images of the Good Samaritan Hospital were captured with a Canon 40D camera ( pixel pitch) using a 300-mm telephoto lens ( 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 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.
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 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.
Path Weighting Functions
The expected variance of the wavefront tilt was evaluated using an approach outlined by Fried10 and Winick and Marquis11 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 -tilt over an aperture of diameter , when viewing a source in the direction can be expressed as Ref. 12
The mean correlation between tilts observed at the aperture due to two sources at viewing directions and can be written as follows:
Hence, Eq. (4) can be expressed as
For a spherical wave propagating through turbulence characterized by the Kolmogorov power spectrum, the phase structure function can be written as in Ref. 1310 and Winick and Marquis.11 Applying those techniques, Eq. (7) reduces to
Each pixel in the time-lapse imagery corresponds to a patch of , 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 is defined as
Hence, the patch-averaged tilt variance is
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),
Changing the variables of integration, Eq. (11) reduces to
Now, , where and . Substituting the above result in Eq. (13), the expression for patch-averaged tilt variance becomes
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 for several different patch sizes on the target.
The aperture size and the path length used in the evaluation were the same as those used in the experiment, i.e., and . The weight is maximum at the receiver (camera) () and drops to zero at the target (). 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 () in the image is significantly less than that for a point source.
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 . If we do not make this presumption, then a set of tilt variances from differently sized patches can be seen as encoding differences in 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 is given by Refs. 13 and 14
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 as above. In this expression, is weighted by along the path. This is the same weighting as Eq. (15) for in the case where the patch size, , 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, , 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
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 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.
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, , 124.42, , 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.
Example Estimation of
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.
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 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.
From these shifts, 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 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 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 . The technique we wish to use is capable, at least in principle, of overcoming this limitation.
Figure 9 shows the result of using a linear combination of the measured tilt variance with the weights determined above to estimate . This estimate of was computed using Eq. (19), which takes the measured tilt variances ( and in a somewhat shorthand notation) and the weights to generate the result. As before, is the telescope diameter (53 mm in this case). The wavelength of light, , is taken as 500 nm for this experiment.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.
There is no distribution of 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 . 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 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.
Modified Approach for Estimation of
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 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.
Equation 19 will be used again to estimate 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 and , the desired variances. The small patch size variance, , 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, , 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.
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 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.
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 , 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.
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.
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.