Estimating Turbulence Distribution over a Heterogeneous Path Using Time-Lapse Imagery from Dual Cameras

: Knowledge of turbulence distribution along an experimental path can help in effective turbulence compensation and mitigation. Although scintillometers are traditionally used to measure the strength of turbulence, they provide a path-integrated measurement and have limited operational ranges. A technique to proﬁle turbulence using time-lapse imagery of a distant target from spatially separated cameras is presented here. The method uses the turbulence induced differential motion between pairs of point features on a target, sensed at a single camera and between cameras to extract turbulence distribution along the path. The method is successfully demonstrated on a 511 m almost horizontal path going over half concrete and half grass. An array of Light-Emitting Diodes (LEDs) of non-uniform separation is imaged by a pair of cameras, and the extracted turbulence proﬁles are validated against measurements from 3D sonic anemometers placed along the path. A short-range experiment with a heat source to create local turbulence spike gives good results as well. Because the method is phase-based, it does not suffer from saturation issues and can potentially be applied over long ranges. Although in the present work, a cooperative target has been used, the technique can be used with non-cooperative targets. Application of the technique to images collected over slant paths with elevated targets can aid in understanding the altitude dependence of turbulence in the surface layer.


Introduction
For compensation and mitigation schemes to work effectively, it helps to know how turbulence is distributed along a path.Scintillometers have been the gold standard for measuring turbulence [1], but scintillometers provide an integrated turbulence measurement, and have limited range capability.Several irradiance-based as well as phase-based techniques to profile turbulence has been discussed in the literature.The Scintillation Detection and Ranging (SCIDAR) technique has been used by astronomers to vertically profile turbulence [2].The technique uses correlation of scintillation due to a binary star pair at the pupil plane of the telescope.A phase-based technique known as Slope Detection and Ranging (SLODAR) profiles turbulence by measuring the cross-correlation of wavefront tilts due to a binary star pair using a Shack-Hartmann sensor [3].Whiteley et al. introduced another phase-based method known as Difference in Differential Tilt Variance (DDTV) that relied on differential tilt measurements between two receiver apertures due to three beacons [4].Gladysz et al. have also used differential tilt measurements at a single receiver due to an array of LEDs to characterize turbulence along a horizontal path [5].
In an earlier work, a time-lapse imaging technique was introduced to measure integrated turbulence parameters remotely from a single site [6].The technique was applied to images collected over a 7 km slant path from Fitz Hall at University of Dayton to Dayton VA Medical Center.By tracking the motion of window corners on the hospital building, path-weighted estimates of the refractive index structure constant, 2 n C were evaluated.
The estimates were in good agreement with co-located scintillometer measurements.Time-lapse estimates from a second imaging experiment conducted over a 1 km horizontal path close to the ground with cooperative and non-cooperative targets at the Laser Experimental Range at Wright Patterson Air Force Base agreed with scintillometer measurements as well.
In all of the above-mentioned experiments in time-lapse, a single camera was used to capture images of the scene which enabled estimation of integrated turbulence parameters.In the present work, two spatially separated cameras were used to image an array of LEDs over a 511 m heterogeneous path, half being concrete, the rest grass.The differential tilt variances due to a pair of LEDs sensed by a single camera and between two cameras form a rich set of measurements from which turbulence distribution along the path can be extracted.The profiling results were compared to 3D sonic anemometers placed every 100 m along the path.A second experiment was conducted over a shorter grassy path with a gas heater placed in the middle of the path to generate a localized spike in turbulence.A poster board with white circles against a black background was imaged by the two cameras.Both the experiments were done over an approximately horizontal path very close to the ground.It is worth mentioning here that commercially available systems such as DELTA, developed by MZA Associates Corporation profiles turbulence along a path from time-lapse imagery of a single camera [7].A richer set of weighting functions, and hence better prospect at profiling can be obtained by use of multiple cameras, as described in this work.The imaging approach is also a low cost, portable approach to profiling when compared to more sophisticated systems such as MZA's PROPS [8].Methodology for the present work is described in Section 2. Theoretical expressions for the path weighting functions, which are an important component of this work have been derived here.Section 2 also provides a description of the imaging experiments.Profiles derived from timelapse imagery measurements are compared to anemometer measurements in Section 3. A summary of the work and future research directions is discussed in Section 4.

Weighting Functions
Given the size of the LEDs was only 7% of the imaging aperture diameter, they were treated as point sources in the present framework [9].The source-aperture geometry permits 3 different scenarios: (i) sensing by a pair of apertures where the sensing paths cross each other somewhere along the path, (ii) sensing by a pair of apertures such that the sensing paths tend to cross behind the apertures, and (iii) sensing by a single aperture where the sensing paths converge at the aperture.These three geometries are shown in Figure 1.Consider the crossed-path geometry of Figure 1a.The two apertures, each of diameter D are located at position coordinates 1 r and 2 r , respectively, with their centerto-center separation being 2 1 s r -r  .The z-tilt over the 1st aperture when viewing the 2nd source in the direction 2 θ can be expressed as [10].

32
, , where  is the wavelength,    Similarly, the z-tilt over the 2nd aperture when viewing the 1st source in the direction 1 θ can be expressed as: where the angled brackets indicate ensemble averaging.By interchanging the order of integration and ensemble averaging, Equation (4) can be rewritten as, Equation ( 5) can be equivalently expressed as, Since , terms which are functions of either r or r and not both, can be added without changing the result of the integration.Hence, Equation ( 6) becomes, where are the phase structure functions.Considering a spherical wave propagating from each source through turbulence characterized by the Kolmogorov power spectrum, and assuming that the wave structure function is approximately equal to the phase structure function [11], where Δθ θ -θ  is the angular separation between the two sources, wave number and L is the path length.The apertures are located at 0 z = and the sources are located at z L = .To make the analysis simple, it is assumed that the source separation vector, Δθ and the aperture separation vector, s point in the same direction as shown in Figure 1a.
The integrations over r and  r can be done using techniques outlined by Fried [12] and Winick et al. [13].Applying those techniques, Equation ( 9) reduces to The differential tilt variance is thus a path-weighted integral of C , the path weighting function being.
The path weighting function for the differential tilt variance due to a pair of sources in the sensing configuration shown in Figure 1b is obtained by replacing s -Δθz in Equation (11) with s Δθz  : Here, Δθ is the angular separation between the two sensing paths shown in Figure 1b.
The differential tilt variance due to the two sources sensed by a single aperture is also a path-weighted integral of 2 n C , the path weighting function being [14], where Δθ is the angular separation between the two sources at a single aperture.Note that Equation ( 13) can be obtained from Equation (11) by setting the aperture separation to be zero.15 different source separations result out of the LED array configuration used in the longer range, mixed path experiment.Figure 2 shows the path weighting functions of differential tilt variance corresponding to the 15 different source separations for the 3 sensing scenarios described in Figure 1.The weighting functions for all three cases drop to zero at the source end of the path.This implies that the apertures cannot sense the turbulence at the source end.The weighting functions for case (i), i.e., the crossed-path geometry, all have a characteristic notch, where they dip to zero.The notch appears where the two sensing paths from the LEDs to the apertures cross along the path.The tilt contribution due to turbulence at this location is exactly the same for both paths, and hence has zero weight on the differential signal.The locations of the notches change with LED separation.As the separation between the LEDs increase, the notch locations move closer to the aperture end of the path.The single aperture weighting functions have characteristic peaks whose locations depend on the separation between the pair of LEDs.Larger the separation, the closer is the peak to the aperture.Since in this geometry, the two sensing paths converge at the aperture, the effect of turbulence at the aperture is canceled.This results in zero weighting at the aperture end in Figure 2c.The weighting functions for case (ii) show less diversity, as can be seen in Figure 2b.All three sets of weighting functions have very similar values beyond 370 m.This clearly suggests that the present technique will fail to profile turbulence for z > 370 m.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, 2 n C along the path can be estimated: where M  is the pseudo-inverse of M and V is a set of measured differential tilt vari- ances.Equation ( 14) can be rewritten as: Equation ( 15) shows a relationship between the estimated and actual values of along the path.The matrix M M  is in essence the influence function matrix, or the impulse response matrix.Each column of the matrix describes the response of the estimation method to a turbulence impulse at that range.It also describes how turbulence estimated at a particular location is influenced by turbulence elsewhere along the path.Figure 3 shows these influence functions at 5 different locations for the mixed path experiment.Near the apertures (Figure 3a), the profiling resolution is satisfactory, and improves between 120 m-220 m as the influence function starts getting narrower (Figure 3b).The resolution again starts degrading as we move closer and closer to the source end (Figure 3c).The nature of the influence functions at and beyond 370 m (Figure 3d,e) clearly suggest that the present technique fails to profile turbulence here.

Dual-Camera Time-Lapse Imaging Experiment
Two experiments were conducted to demonstrate the dual-camera profiling approach.The first experiment was done at the Laser Experimental Range at Wright Patterson Air Force Base in Dayton, Ohio over 2 days in June and 1 day in July of 2019.The 511 m experimental path was approximately horizontal and about 1.5 m above the ground.An LED target board constructed using several 100 W, 5 mm diameter LEDs of wavelength 850 nm and varying separations was placed on one end of the path.The LED array configuration resulted in 15 different source separations of 4, 6, 8, 10, 14, 16, 20, 24, 30, 36, 40, 44, 50, 60 and 80 cm.At the other end, two science cameras (FLIR GS3-U3-89S6M-C), fitted with Nikon lenses of focal length 300 mm (f/#4), were mounted on a tripod within a trailer.The cameras were spatially separated by 10 cm.The first 250 m of the imaging path from the cameras was concrete, and the remainder was grass.The experimental path and the setup are shown in Figure 4.A set of 10 short exposure images of the LED board was captured by both cameras every 10 s.Due to changes in brightness conditions throughout the day, the exposure times of the cameras had to be varied, but care was taken to keep the exposures short such that the turbulence was frozen during each frame capture.Four Sonic Anemometers from Applied Technologies, Inc. were set up every 100 m along the path.Two anemometers were thus on grass and two were on concrete.The anemometers measured sonic temperature and wind speed in three directions, and a 2 n C value could be derived from the temperature structure constant obtained from these measurements at each location.

Mapdata: Google
The images were cropped around the LED board and the cropped images were 250 pixels × 250 pixels.Figure 5 shows a cropped image captured by one camera at 14:50 local time (UTC-4) on 27 June 2019.Due to camera setup imperfections, the images from the two cameras appeared to be rotated by different amounts.To estimate the angle of rotation and to correct for it, a mean frame was generated by averaging 300 frames captured over a 5-min period.Centroids were found for 12-pixel × 12-pixel blocks surrounding each of the five spots in the top row of the pattern.A linear fit was done on the centroid locations, and the angle of rotation was determined from the slope of the line fit.This exercise was done for the images captured by each camera.Next, x-and y-centroids (in pixels) were computed for every 12 pixels × 12 pixels blocks surrounding each spot on every frame.The centroids were corrected for the rotation.Angular tilts were obtained by multiplying the x-and y-centroids (in pixels) by the ratio of the camera pixel pitch to the focal length.In weak turbulence conditions, each spot in the image could be clearly identified.However, when the turbulence grew strong, the spots with the least separation merged, making it difficult to track these spots.Hence, only 7 out of the 15 separations could be exploited during periods of strong turbulence.
Difference between tilts obtained from a single camera or between the two cameras due to every pair of sources with a specific separation were computed.This was repeated for all 300 frames.The differential tilt variances were computed across frames for each source pair.The x-and y-variances were added and averaged across all source pairs that had the correct separation, resulting in 15 different variances, corresponding to the 15 separations, for each of the 3 sensing geometries of Figure 1.The 45 different weighting functions shown in Figure 2 were sampled at 0.5 m intervals, so a 45 × 1023 matrix M was generated.By using the Moore-Penrose technique, a pseudo-inverse of M was computed with a threshold.The thresholding was done to suppress noise effects, since it set all the singular values of M less than the threshold to zero.The pseudo-inverse matrix along with the computed tilt variances were used to profile turbulence along the path.Through this inversion process several unrealistic negative values of  The second experiment was over a shorter range.The 48 m experimental path was once again approximately horizontal, and over grass, about 1.7 m above the ground.The experiment was conducted during the evening hours of 15 November 2019 in the front yard of a house in Beavercreek, OH.The target board in this experiment consisted of 5 rows of 0.8 mm diameter white dots, with varying separations, printed against a black background.The imaging system was the same as that used in the previous experiment.At 18 m from the cameras, a propane heater was placed.The purpose of the short-range experiment was to further test the dual-camera profiling approach by checking if the extracted profile would show the local spike in turbulence, caused by the propane heater.An LED source illuminated the target board and the cameras captured 10 frames every 15 s, with an exposure time of 30 micro-seconds.The experimental setup is shown in Figure 6.The arrangement of dots on the target board resulted in 81 different separations, and hence 243 weighting functions, corresponding to the 3 sensing geometries could be obtained.The weighting functions for the short-range experiment are shown in Figure 7.The differential tilt variances were computed from the centroid motion of the dots, in the same way as in the long-range experiment.Only the top row of dots was considered, owing to concerns that turbulence might not be fully developed closer to the ground, where the heater was.Turbulence was very strong, too for the lower rows, resulting in neighboring dots merging together.The straight up matrix inversion method again gave some negative values of 2 n C , and hence the constrained optimization technique was used with initial guess values from the matrix inversion method.They follow the same trend, and the values are reasonably close.It is important to note here that due to a power issue, one of the 4 anemometers stopped recording during some data collections.It is apparent from Figure 8 that the present technique fails to profile turbulence beyond 370 m from the cameras.As mentioned before, the weighting functions start looking very similar beyond this point.In strong turbulence conditions, the spots with the smallest separation merged, making it impossible to track their motion.This resulted in only 7 usable separations, and hence 21 weighting functions that could be used for extraction of profile in the strong turbulence cases.The reduced set of weighting functions for the strong turbulence case makes the technique fail even earlier.Thus, it is not possible to compare the 2 n C value obtained from the anemometer at 400 m with a corresponding estimate derived from timelapse measurements.The profiles from time-lapse measurements (red trace) in Figure 8b,c have 2 n C values lower than those from the anemometers.Both these profiles correspond to times when the turbulence was strong and hence the assumption that the phase structure function is approximately equal to the wave structure function does not hold good.This assumption used in the derivation of the weighting functions probably caused the underestimation of 2 n C in these cases.However, the differences are not significant.Although effects due to any common mode disturbances were eliminated by using a differential tilt signal, any difference in the motion of the two cameras could have affected the results.During strong turbulence situations, it was difficult to get the cameras focused.This defocus could have been an additional source of noise.

Short-Range Imaging Experiment
Figure 9 shows the profiles from the short-range experiment.The red and the blue traces correspond to the matrix inversion method and the constrained optimizer, respectively.Again, the straight up matrix inversion method gives noisy results, with several negative values of 2 n C .In Figure 9, the negative values have been replaced by a floor of 5 × 10 -13 m -2/3 .The optimizer performs better, and clearly shows a spike at 18 m from the cameras, where the gas heater was placed.Different data sets were processed and they all showed this spike in turbulence at 18 m, thus validating the dual camera profiling technique.It is evident from Figure 9 that the technique fails to profile turbulence beyond 30 m from the cameras.

Discussion
A method for obtaining turbulence distribution along a path using time-lapse imagery from a pair of spatially separated cameras was described.Given that the method is phase-based, it has the potential to be applied over longer paths, where irradiance-based methods fail due to saturation issues.Another advantage is it is a low-cost, portable approach, capable of remote estimation from a single site without deployment of sensors or sources at the target location.Although cooperative targets have been used in this work, the method can work with targets of opportunity as long as there are enough trackable features on the target.
Two experiments were conducted to verify this approach.One experiment was over a heterogeneous path of part grass and part concrete.The time-lapse estimates were compared to sonic anemometer derived 2 n C values and there was good agreement.The change in 2 n C from concrete to grass was however, not dramatic.A second shorter range experiment was conducted over grass with a gas heater placed in the middle of the path.The idea was to verify if the profiling technique would show the local spike in turbulence due to the gas heater.The profile extracted from the time-lapse imagery showed the spike in turbulence at the location of the heater.The technique, however, fails to profile turbulence close to the source end of the path.
The straight up matrix inversion technique to extract the profiles gives noisy results, with several negative values of 2 n C .A constrained nonlinear optimization technique was used to improve upon the profiles obtained from the matrix inversion method.The results can be improved by refining the present optimization technique or more sophisticated algorithms can be used that are more robust to noise.It was difficult to get the cameras focused during periods of strong turbulence and this introduced additional errors.The way the two cameras were mounted in the long-range experiment, differential motion between them was possible.This could have negatively affected the results as well.The assumption in the mathematical model that the wave structure function is approximately equal to the phase structure function probably resulted in underestimation of 2 n C during strong turbulence conditions.Results can be greatly improved by focusing the cameras when turbulence levels are low and by using a stable imaging platform such that the differential motion between cameras is minimized.Multiple spatially separated cameras, instead of a pair of cameras can help in improving the profiling resolution as well as the fraction of the path over which meaningful 2 n C estimates can be obtained.
The dual-camera profiling technique can be useful in studying how turbulence changes with altitude in the surface layer.There is very limited information about the surface layer and these type of measurement campaigns can help improve existing atmospheric models.Future work will include imaging elevated targets over slant paths to understand the altitude dependence of turbulence.Multiple spatially separated cameras will also be used to profile over longer ranges using targets of opportunity.


is the turbulence induced wavefront distortion at aperture coordinate r due to the 2nd source and,

Figure 1 .
Figure 1.Imaging path geometries.(a) Sensing paths cross; (b) Sensing paths do not cross; (c) Sensing at a single aperture.

Figure 2 .
Figure 2. Path weighting functions of differential tilt variance for 15 different source separations in the long-range experiment.The weighting functions are for the imaging path geometries shown in Figure 1.(a) Sensing paths cross; (b) Sensing paths do not cross; (c) Sensing at a single aperture.

Figure 4 .
Figure 4. Long-range experimental setup.Stars denote the location of the anemometers.

2 nC 2 nC
were generated.A constrained nonlinear optimization technique was then applied to drive the profile towards having only positive values of The objective function to be minimized was the squared difference between the measured differential tilt variances and those arising from a given choice of 2 n C values along the path.The optimizer's initial guess values for 2 n C were the same as those obtained through the matrix inversion technique and the negative values were replaced by a floor.The choice of floor was based on the range of 2 n C values obtained from the matrix inversion technique and it varied depending on weather conditions and time of day.An upper and a lower bound were set as well, based on the estimates from the matrix inversion technique.

Figure 5 .
Figure 5.A sample cropped image of the LED board captured by one of the cameras at 14:50 local time (UTC-4) on 27 June 2019.

Figure 6 .
Figure 6.Experimental setup for the short-range experiment.(a) Experimental path; (b) The gas heater used to create a local spike in turbulence; (c) Target board.

Figure 7 .
Figure 7. Weighting functions corresponding to 3 sensing geometries and 81 separations on the target board used in the short-range imaging experiment.

Figure 8
Figure 8 shows 4 example profiles from the long-range experiment.The blue trace corresponds to the profile extracted using the matrix inversion technique and the red trace corresponds to the profile from the optimizer.The black diamonds correspond to the 2 n C outputs from the 4 anemometers.Gaps in the blue trace are due to negative values of 2 n C , which sometimes appear in the matrix inversion method.The oscillations in the blue trace are due to noise from the measurement and from the profile extraction method.The constrained optimization technique performs better in presence of noise.Overall, there is good agreement between the optimizer profiles and the anemometer derived 2 n C values.

Figure 9 .
Figure 9. Turbulence profiles from the short-range imaging experiment.Profiles were derived from images captured at 20:05 local (UTC-5) on 15 November 2019.