Estimation of atmospheric turbulence using differential motion of extended features in time-lapse imagery

Abstract. Accurate characterization of atmospheric turbulence is useful for performance assessment of optical systems operating in real environments and for designing systems to mitigate turbulence effects. Irradiance-based techniques such as scintillometry, suffer from saturation, and hence commercial scintillometers have limited operational ranges. A method to estimate turbulence parameters, such as path weighted Cn2 and Fried’s coherence diameter r0 from turbulence-induced random, differential motion of extended features in the time-lapse imagery of a distant target is presented. Since the method is phase-based, it can be applied to longer paths. It has an added advantage of remotely sensing turbulence without the need for deployment of sensors at the target location. The approach uses a derived set of path weighting functions that drop to zero at both ends of the imaging path, the peak location depending on the size of the imaging aperture, and the relative sizes and separations of the features whose motions are being tracked. Using different sized features separated by different amounts, a rich set of weighting functions can be obtained. These weighting functions can be linearly combined to approximate a desired weighting function such as that of a scintillometer or that of r0 in inhomogeneous turbulence. The time-lapse measurements can thus mimic the measurements of a scintillometer or any other instrument. The method is applied to images captured along two different paths, and the estimates are compared to co-located scintillometer measurements.


Introduction
Accurate characterization of atmospheric turbulence and its effects is extremely important for performance evaluation of optical systems operating in real environment and for designing of systems to mitigate turbulence effects. Irradiancebased methods such as the scintillation, detection, and ranging technique have been used in the past by astronomers to obtain a low-resolution vertical profile of turbulence. 1,2 The path-weighted refractive index structure constant C 2 n is traditionally measured using scintillometers. 3 However, irradiance-based techniques are of limited use in the saturation regime and hence are not suitable for measurements over long paths through turbulence. Also direct point estimates of C 2 n that are derived from such intermediate quantities as the velocity structure constant C 2 v and the temperature structure constant C 2 T require measurements of wind speed and temperature at high-temporal resolution. 4 Use of these methods requires physical sensors to be deployed in and around the region of interest. Alternate techniques that can provide reliable turbulence information over paths through strong turbulence are being investigated. 5 Meier and Fiorino 6 have demonstrated methods to extract turbulence information from polar orbiting satellite data. A phase-based technique built on a process known as the difference in differential tilt variance was introduced by Whiteley et al. 7 The technique requires two receiver apertures and an accompanying focal plane camera at one end of the path and three point source beacons at the other end. The apertures were used to observe either a singlepoint source or two individual point sources. Since differential tilt measurements are phase related, conventional propagation theory can be used even for strong turbulence paths where the variance of irradiance saturates. Gladysz et al. 8 measured turbulence strength over a horizontal path from the differential angle of arrival fluctuations of an array of light emitting diodes.
An imaging experiment was done at the Air Force Institute of Technology (AFIT) in the summer of 2014 to monitor the effects of atmosphere over a period of time. A digital single-lens reflex camera fitted with a telephoto lens was used to capture images of a hospital building 12.8 km away. The imaging path was almost horizontal, with most of the path 60 m above the ground. Path-averaged C 2 n was estimated from turbulence-induced random motion of different parts of the building. 9,10 In this work, a scheme to directly estimate turbulence parameters such as pathweighted C 2 n and Fried's coherence diameter r 0 , from the differential motion of target features is described. The differential mode of measurement cancels the effect of any common mode disturbance such as platform vibrations. In addition to being phase-based, this technique has the added advantage of being able to remotely measure turbulence from a single location. Some commercially available systems estimate turbulence by tracking differential motion of feature points on the target. 11 Tracking the motion of extended features on the target rather than point features allows measurement of C 2 n over a longer range. In this paper, the differential technique has been applied to images captured along two different propagation paths.
The remainder of this paper consists of three sections: methodology, results, and conclusions. Section 2 describes the experimental configuration, derives the path weighting functions that relate the turbulence to the differential tilt variances, and demonstrates how a linear combination of multiple weighting functions can be used to approximate a desired weighting function, such as that of a scintillometer or that of r 0 in inhomogeneous turbulence. Section 3 discusses how the ideas developed in Sec. 2 can be applied to actual experimental imagery in order to obtain C 2 n estimates. These estimates are compared to scintillometer measurements taken along the same path. A brief summary of the findings and future research directions are given in Sec. 4.

Time-Lapse Imaging Experiment
Two time-lapse imaging experiments were done over two different ranges at different times of the year. In the first experiment, conducted in February 2017, images of the Dayton VA Medical Center were captured from a window at University of Dayton's (UD) Intelligent Optics Laboratory. Figure 1 shows the layout and a ground elevation profile of the 7-km path from the VA Medical Center to UD. Example images are shown in Fig. 2. The images were captured and saved every 40 s using a Manta G609 Allied Vision Technologies camera with a pixel pitch of 4.54 μm, mounted at the back of a 14-in. (35 cm) telescope with focal length of 3.91 m. The exposure time used was 10 ms. Turbulence was estimated by measuring the motion of four patches, marked as white circles in Fig. 2(a), each centered at one of the corners of the two windows in the right wing of the hospital. A Scintec BLS 2000 scintillometer measured turbulence along nearly the same path. The BLS transmitter was a floor above the two tracked windows at the VA Medical Center and can be seen as the bright source at the center of Fig. 2(a).
Although the goal of the first experiment was to estimate turbulence by looking at targets of opportunity, the second experiment used specially designed targets over a much shorter path close to the ground. The experiment was conducted over 3 days in July 2017 at the Laser Experimental Range at Wright Patterson Air Force Base, Ohio. 12 The 1-km path was nearly horizontal over an asphalt runway ∼1.5 m above the ground. Targets were four posters with white circles spray painted randomly on a black background. The imaging system comprised a Lumenera camera (pixel pitch 4.65 μm) and a Nikon lens with focal length 400 mm (system ƒ# 2.8). The targets were affixed to the side of a trailer and images were captured every 10 s. A BLS 900 scintillometer was used to monitor the same path. Figures 3 and 4 show the experimental path and sample images from July 18 and July 20, respectively. The exposure time on the camera was set to 500 μs on July 18 and 5 ms on July 20. The 5-ms exposure time caused the images taken later in the day on July 20 to be badly saturated and unusable for reliable turbulence estimation. Strong winds on July 18 blew away the target posters during the course of the experiment, and turbulence was then estimated by tracking the differential motion of four patches, each centered at one of the corners of the trailer window.
The images were first cropped to isolate the region of interest. In the VA Medical Center imaging experiment, the region of interest was the right wing of the hospital building. The trailer was the region of interest for the runway experiment. A reference image was chosen from the images collected and a cross-correlation algorithm was used to estimate the image shift between each image and the reference image. A Gaussian window was applied to the images before the cross-correlation algorithm was run. This reduced the effects of the frame edges on the correlation result.    Bose-Pillai et al.: Estimation of atmospheric turbulence using differential motion of extended features in time-lapse imagery A parabolic fit was applied to the correlation peak to provide a subpixel estimate of the shift between the images. The correlation with a reference image provided information about slow image drift during the course of the day due to refractive bending. This long-term drift information was used to adjust the tracking window keeping it locked on to a feature through the collection. The cross-correlation algorithm was then applied to each image and its neighboring image to estimate the random motion of the feature due to turbulence. The shifts of two features were subtracted to get the differential signal.

Path Weighting Functions for Differential Patch-Averaged Tilt Variance
The z-tilt over an aperture of diameter D, when viewing a source in the direction θ, can be expressed as 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 1 ; 6 3 ; 5 7 1 αðθÞ ¼ 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 ; 5 The mean correlation between tilts observed at the aperture due to two sources at viewing directions θ 1 and θ 2 can be written 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 3 ; 6 3 ; 4 2 8 where the angled brackets indicate ensemble averaging.
Interchanging the order of integration and ensemble averaging results in 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 ; 6 3 ; 3 2 2 Since RR drdr 0 WðrÞWðr 0 Þr · r 0 ¼ 0, terms that are functions of either r or r 0 , and not both, can be added without changing the result of the integration.
Hence, Eq. (4) becomes: 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 ; 6 3 ; 2 1 6 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 wave structure function can be written as 14 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 ; 7 5 2 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. Assuming that the phase structure function is approximately equal to the wave structure function and substituting Eq. (6) into Eq. (5): ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 6 2 7 The integrations over r and r 0 can be done using techniques outlined by Fried 15 and Winick and Marquis. 16 Applying those techniques, Eq. (7) 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 0 8 ; 3 2 6 ; 5 0 4 Each pixel in the time-lapse imagery corresponds to a patch of finite size (∼8 mm for the VA imaging experiment and ∼1 cm for the runway experiment) 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 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 as 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 ; 3 2 6 ; 2 5 3 where P G ðθÞ ¼ e − 4θ 2 L 2 d 2 , d being the 1∕e patch diameter and θ ¼ jθj.
The correlation between two patch-averaged tilts α p1 and α p2 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 ; 3 2 6 ; 1 6 3 where Δθ is the angular separation between the two patches. The order of integration and ensemble averaging is interchanged to obtain the above expression. In the following analysis, the patches are assumed to be equal in size. The term with the angle brackets in Eq. (10) is nothing but the correlation of tilts due to two point sources. Substituting Eq. (8) into 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 ; 7 1 9 (11) can be equivalently 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 1 2 ; 6 3 ; 5 4 0 Let 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 ; 6 3 ; Changing the variables of integration, Eq. (12) 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 4 ; 6 3 ; 3 1 8 Now, , where x ¼ jxj and y ¼ jyj.
Substituting the above result in Eq. (14), the expression for patch-averaged tilt correlation 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 5 ; 3 2 6 ; 7 5 2 The expression for patch-averaged tilt variance hα 2 p i derived in Refs. 9 and 10 can be obtained from Eq. (15) by setting Δθ ¼ 0. The differential patch-averaged tilt variance can hence 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 1 6 ; 3 2 6 ; 5 3 7 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 7 ; 3 2 6 ; 4 6 6 is the path weighting function. No simpler form for the weighting function could be obtained, and hence in the present work, Eq. (17) was evaluated numerically for cases of interest. Figure 5 shows a plot of f d ðzÞ for different patch sizes and separations. The aperture size and the path length used in the evaluation were the same as those used in the VA Medical Center imaging experiment, i.e., D ¼ 35 cm and L ¼ 7 km. The weighting function drops to zero at both ends of the path. This implies that the turbulence near the target hardly affects the differential motion seen by the camera. The tilt due to turbulence at the camera is the same for both patches since the two sensing paths converge at the camera. This explains the zero weighting of turbulence at the camera end of the path. The time-lapse method, thus, is not sensitive to turbulence variations either at the source or at the camera. The peak locations of the weighting functions depend on the patch size and separation relative to the imaging aperture. Subaperture sized patches and separations have weighting functions that peak toward the target end of the path. Larger patch sizes and separations have weighting functions that peak toward the camera end of the path.

Creating Arbitrary Weighting Function from Linear Combination of Path Weighting Functions
If the turbulence is presumed to be constant along the path, then the differential tilt variance associated with each pair of image patches would provide an estimate for C 2 n . If this presumption is not made, a set of tilt variances from pairs of patches of different sizes and separations 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 and separated image patches, each characterizes the turbulence along (almost) the same path, but each weights that path differently. However, this rich set of weighting functions can be used as a basis set such that the weighting functions from this set can be linearly combined to reproduce the path weighting function associated with some atmospheric parameter of interest. To determine the appropriate linear combination of weighting functions, the Moore-Penrose pseudoinverse technique can be used to find the least-squares optimal way to achieve the desired weighting function in terms of the basis set available. This technique was used in Ref. 10 to get an unbiased estimate of r 0 from the motion of image patches of two different sizes.
Here it is shown that the weighting functions for differential patch-averaged tilt variances can be used to better reproduce the weighting functions for the spherical wave r 0 and a scintillometer.
The spherical wave r 0 is given by 17 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 ; 6 3 ; 1 3 9 r 0;sw ¼ 0.423k 2 This equation has been written so the integration takes place from the camera to the target, i.e., the camera is at 0 and the target 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 that of tilt variance due to a point source. Since there is no beacon or point source at the target, no part of the image corresponds to a weighting function of this form. By sampling the weighting functions along the path and generating a matrix M, where the rows are formed from the individual weighting functions and the columns correspond to the ranges where these functions are sampled, the desired weights can be computed as where W is the weight to be applied to each weighting function, 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. The Moore-Penrose technique is used here to obtain the pseudoinverse. Figure 6 shows how three different weighting functions corresponding to pairs of image patches of different sizes and separations are linearly combined to reproduce the weighting function for r 0;sw . The aperture size and the path length used in the evaluation are the same as those used in the VA Medical Center imaging experiment. The weighting function obtained from linear combination of weighting functions, matches the desired weighting function well, except for a small portion of the path close to the camera. So unless there is significantly strong turbulence near the camera, it is possible to obtain a fair estimate of r 0 from the differential motion of these pairs of patches.
The same technique can be used to obtain a weighting function that closely matches that of a scintillometer. The approximate functional form of the weighting function for Scintec BLS scintillometers can be found in Ref. 18. Figure 7 shows how weighting functions corresponding to pairs of patches of two different separations from the VA Medical Center images can be linearly combined to approximate the scintillometer weighting function. Except for the region near the camera, the two weighting functions agree reasonably well. This weighting function has been used with the VA Medical Center imagery in the following section to obtain C 2 n estimates that compare directly to scintillometer measurements. Fig. 5 Path weighting functions of differential patch-averaged tilt variances for different patch sizes and separations. The camera is at z ¼ 0 and the target is at z ¼ 7000 m.

VA Medical Center Imaging Experiment: Comparison of Time-Lapse Estimates to Scintillometer Measurements
Four patches of 1∕e diameter 16 cm (20 pixels), each centered at one of the corners of the two windows in the right wing of the building (shown in Fig. 2), were tracked. Hence there were two pairs of patches separated by 48 cm, and two pairs of patches separated by 1.8 m at the target for each image. A 30-frame moving window, corresponding to 20 min of imagery, was used to compute the variance from the differential motion. The horizontal and vertical variances were added to obtain the total variance. The tilt variances for the two different patch separations were multiplied by their corresponding weights (computed using the pseudoinverse technique mentioned in Sec. 2) and added together to obtain a path-weighted estimate of C 2 n . In essence, the path-weighting function for the C 2 n estimate was the weighting function shown in Fig. 7(b). The C 2 n estimates were  Optical Engineering 104108-7 October 2018 • Vol. 57 (10) compared to the BLS 2000 scintillometer measurements. Figures 8-11 show the results of these comparisons for four days over different weather conditions. According to meteorological observations, 19 February 14, 2017 was a clear day. As evident from Fig. 8, the time-lapse estimates agree very well with the scintillometer measurements. Throughout the experiment, the two windows being tracked were in the shadow for most of the day. Poor contrast during sunrise and sunset times and shadows sweeping across the windows during the early mornings [as seen in Fig. 2(a)] contributed to estimation error during these times. February 15 was clear earlier during the day, but it became cloudy by midmorning. Overcast conditions prevailed from around noon to 1:30 PM (local standard time, UTC-5 h), causing the dip in C 2 n , seen in Fig. 9. Again, there is excellent agreement between the time-lapse estimates and the scintillometer measurements. February 18 was marked by alternating cloudy and clear conditions throughout the day. This resulted in multiple peaks and troughs in the C 2 n profile, as seen in Fig. 10. The time-lapse estimates and scintillometer show agreeable trends and values. However, due to software issues, the camera failed to capture images intermittently, causing gaps in the time-lapse C 2 n profile. Figure 11 shows a comparison between time-lapse estimates and scintillometer measurements for February 21. This was a cloudy day, with clear conditions developing very late in the afternoon. The time-lapse estimates agree reasonably well with the scintillometer, the gaps in the profile suggesting  again lack of data due to software issues during collection. Poor contrast during cloudy conditions could have contributed to the slight differences in the two profiles. Since the path was nearly horizontal over the same surface type, it was presumed that C 2 n was constant along the path. The path-averaged C 2 n was obtained by dividing the total variance by the area of the corresponding weighting function Fig. 10 Comparison of C 2 n estimates from time-lapse imagery with C 2 n from scintillometer: February 18, 2017. Differences in the two profiles are possibly due to estimation errors resulting from moving shadows, contrast change in neighboring images due to changing cloud cover, and poor image contrast during early morning and late afternoon. Fig. 11 Comparison of C 2 n estimates from time-lapse imagery with C 2 n from scintillometer: February 21, 2017. Differences in the two profiles are possibly due to estimation errors resulting from moving shadows, contrast change in neighboring images due to changing cloud cover, and poor image contrast during early morning and late afternoon.

Runway Imaging Experiment: Comparison of Time-Lapse Estimates to Scintillometer Measurements
Optical Engineering 104108-9 October 2018 • Vol. 57 (10) for that separation. For July 20, 12 patches of 1∕e diameter 16 cm (14 pixels) from different parts of the 4 target posters were used. This resulted in four pairs of patches of separation 1.16 m (100 pixels) and four pairs of patches with separation 4.79 m (412 pixels) from each image. A 60-point moving window was used to compute the variance. Figures 13 and 14 show comparisons of the time-lapse estimates to the scintillometer measurements. Just as in the VA Medical Center experiment, there is good agreement here too. If C 2 n were constant along the path, all three estimates in Figs. 13 and 14 would be the same. The difference in estimates can be attributed to slight variations in C 2 n along the path, in addition to estimation noise. The estimation noise can be due to blurring of features, moving shadow effects [as evident in Fig. 4(a) with the overhang casting a shadow], and changes in contrast between adjacent images as cloud shadows pass over the target. The over-exposure in the images of July 20 led to some estimation error as well.

Conclusions
A method for obtaining turbulence information from differential motion of extended features in the time-lapse imagery of a distant target was described. A potentially significant advantage of this method, although not tested in the Fig. 12 Path weighting functions of differential patch-averaged tilt variances for three different patch separations. The camera is at z ¼ 0 and the target is at z ¼ 1000 m. Fig. 13 Comparison of C 2 n estimates from time-lapse imagery with C 2 n from scintillometer: July 18, 2017. Estimation errors are possibly due to moving shadows, blurring of features, and changes in contrast between neighboring images due to changing cloud cover around noon. experiments described here is that it is phase-based, and hence can be applied to long paths through turbulence where irradiance-based techniques suffer from saturation problems. Additionally, it allows turbulence to be measured remotely from a single site. The variance of the differential image motion was shown to be a weighted integral of C 2 n over the propagation path. An analytical expression for the path weighting function was developed and was shown to depend on the aperture size, path length, and the feature sizes and separations. These weighting functions form a rich set of a family of functions that tapers to zero at both ends of the measurement path. Weighting functions for different patch sizes and separations can be linearly combined to approximate any desired weighting function, such as that of Fried's coherence diameter r 0 or that of a scintillometer. The time-lapse measurements can thus mimic the measurements from a scintillometer or any other instrument. The technique was applied to time-lapse images captured along two different paths. The path-averaged estimates of C 2 n from the experimental data agreed very well with scintillometer measurements, even when noncooperative targets were imaged. The slight differences were mainly due to turbulence-induced blurring of features, moving shadows, poor contrast in some images, and contrast changes due to changing cloud cover. With a noncooperative target such as the VA Medical Center, there were few features to track in a single frame and hence several frames were required to obtain tilt variances. This introduced some temporal averaging effects in the estimates.
In the derivation of the analytical expression for differential patch-averaged tilt variance, it was assumed that the wave structure function is equal to the phase structure function. This might influence the results at high scintillation levels. The Gaussian window used on the images was chosen such that it dropped to zero at the edges of the tracking window. Often the Gaussian window might not have sufficient trackable features under it, e.g., a patch centered at a window corner. In such cases, the effective patch size would be much smaller than the size of the Gaussian window used for the analysis, resulting in overestimation of C 2 n . The time-lapse technique works well with noncooperative targets as long as there are enough trackable features on the target. If the technique is used at night for turbulence characterization, either the target needs to be actively illuminated or some light sources should be present in the scene. In the future, a method to better estimate the effective patch sizes will be developed. Simulations will be done with different turbulence strengths to create synthetic imagery. The timelapse technique will be applied on these imageries to get a better understanding of the nature of the estimation errors. Experiments will be done with elevated targets to get an idea of how C 2 n varies with altitude. Meteorological measurements taken simultaneously will provide a better understanding of how turbulence depends on different weather parameters. Additionally, a scheme to use the time-lapse measurements to profile turbulence along the path will be investigated. n estimates from time-lapse imagery with C 2 n from scintillometer: July 20, 2017. Estimation errors are possibly due to blurring of features, saturation effects, and changes in contrast between neighboring images due to changing cloud cover around noon.