Estimation of turbulence from time-lapse imagery

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 parameters.


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; 1,2 differential image motion monitoring; 3 and difference in differential tilt variance. 4 Recently reported techniques have used weather radar 5,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 timelapse 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, 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.

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 × 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.

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 × 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.

Path Weighting Functions
The expected variance of the wavefront tilt was evaluated using an approach outlined by Fried 10 and Winick and Marquis 11 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 1 7 7 αðθÞ ¼ 32λ where λ is the wavelength, ϕðr; θÞ is the turbulence-induced wavefront distortion at aperture coordinate r, and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; 1 0 7 WðrÞ ¼ 1 jrj ≤ 0.5D 0 jrj > 0.5D : The mean correlation between tilts observed at the aperture due to two sources at viewing directions θ 1 and θ 2 can be written as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 7 1 9 2 ZZ dr dr 0 WðrÞWðr 0 Þr · r 0 ϕðr; θ 1 Þϕðr 0 ; θ 2 Þ ; (3) where the angular brackets indicate ensemble averaging.
Interchanging the order of integration and ensemble averaging yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 6 1 6 2 ZZ dr dr 0 WðrÞWðr 0 Þr · r 0 hϕðr; θ 1 Þϕðr 0 ; θ 2 Þi: Because the integrand is an odd function of both r and r 0 , R R dr dr 0 WðrÞWðr 0 Þr · r 0 ¼ 0, terms that are functions of either r or r 0 , but not both, can be added without changing the result of the integration.
Hence, Eq. (4) can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 4 9 5 where Dðr − r 0 ; θ 1 − θ 2 Þ ¼ h½ϕðr; θ 1 Þ − ϕðr 0 ; θ 2 Þ 2 i 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 3 6 9 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 ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 2 6 8 The integrations over r and r 0 can be done using techniques outlined by Fried 10 and Winick and Marquis. 11 Applying those techniques, Eq. (7) reduces to Each pixel in the time-lapse imagery corresponds to a patch of ∼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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 6 3 ; 5 0 with d as the 1∕e patch diameter and θ ¼ jθj.
Hence, the patch-averaged tilt variance is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 6 3 ; 4 2 9 hα 2 p i¼ 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), E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; 3 3 7 ; t e m p : i n t r a l i n k -; e 0 1 2 ; 6 3 ; Changing the variables of integration, Eq. (11) reduces to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 3 2 6 ; 7 4 1 , where x ¼ jxj and y ¼ jyj. Substituting the above result in Eq. (13), the expression for patch-averaged tilt variance becomes E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 3 2 6 ; 5 3 0 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 3 2 6 ; 4 7 6 is the path weighting function. The weighting function for a point source is recovered by letting d → 0 in Eq. (15). The weighting function is seen to depend only on geometry, with the 5∕6th 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.
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.

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 2 n . 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 2 n 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 6 3 ; 3 7 4 r 0;sw ¼ 0.423k 2 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 2 n 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 6 3 ; 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 MW, 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 3 2 6 ; 6 8 6 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.

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. 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.

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.
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.
From these shifts, C 2 n 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 2 n 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. 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 (α 2 p10 and α 2 p20 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.
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.
There is no distribution of C 2 n 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.

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 20m 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.
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 α 2 p10 and α 2 p20 , the desired variances. The small patch size variance, α 2 p10 , 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, α 2 p20 , 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 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.

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.   Optical Engineering 071504-9 July 2017 • Vol. 56 (7)