## 1.

## 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. Irradiance-based 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}_{n}^{2}$ 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}_{n}^{2}$ that are derived from such intermediate quantities as the velocity structure constant ${C}_{v}^{2}$ and the temperature structure constant ${C}_{T}^{2}$ 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 single-point 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}_{n}^{2}$ 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 path-weighted ${C}_{n}^{2}$ 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}_{n}^{2}$ 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}_{n}^{2}$ 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.

## 2.

## Methodology

## 2.1.

### 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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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 $\sim 1.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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. 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.

## 2.2.

### 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 $\mathbf{\theta}$, can be expressed as^{13}

## (1)

$$\mathbf{\alpha}(\mathbf{\theta})=\frac{32\lambda}{{\pi}^{2}{D}^{4}}\int \mathrm{d}\mathbf{r}W(\mathbf{r})\mathbf{r}\varphi (\mathbf{r},\mathbf{\theta}),$$## (3)

$$\u27e8\mathbf{\alpha}({\mathbf{\theta}}_{1})\xb7\mathbf{\alpha}({\mathbf{\theta}}_{2})\u27e9={\left(\frac{32\lambda}{{\pi}^{2}{D}^{4}}\right)}^{2}\u27e8\iint \mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime}W(\mathbf{r})W({\mathbf{r}}^{\prime})\mathbf{r}\xb7{\mathbf{r}}^{\prime}\varphi (\mathbf{r},{\mathbf{\theta}}_{1})\varphi ({\mathbf{r}}^{\prime},{\mathbf{\theta}}_{2})\u27e9,$$Interchanging the order of integration and ensemble averaging results in

## (4)

$$\u27e8\mathbf{\alpha}({\mathbf{\theta}}_{1})\xb7\mathbf{\alpha}({\mathbf{\theta}}_{2})\u27e9={\left(\frac{32\lambda}{{\pi}^{2}{D}^{4}}\right)}^{2}\iint \mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime}W(\mathbf{r})W({\mathbf{r}}^{\prime})\mathbf{r}\xb7{\mathbf{r}}^{\prime}\u27e8\varphi (\mathbf{r},{\mathbf{\theta}}_{1})\varphi ({\mathbf{r}}^{\prime},{\mathbf{\theta}}_{2})\u27e9.$$Hence, Eq. (4) becomes:

## (5)

$$\u27e8\mathbf{\alpha}({\mathbf{\theta}}_{1})\xb7\mathbf{\alpha}({\mathbf{\theta}}_{2})\u27e9=(-\frac{1}{2}){\left(\frac{32\lambda}{{\pi}^{2}{D}^{4}}\right)}^{2}\iint \mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime}W(\mathbf{r})W({\mathbf{r}}^{\prime})\mathbf{r}\xb7{\mathbf{r}}^{\prime}{\mathfrak{D}}_{\varphi}(\mathbf{r}-{\mathbf{r}}^{\prime},{\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2}),$$For a spherical wave propagating through turbulence characterized by the Kolmogorov power spectrum, the wave structure function can be written as^{14}

## (6)

$$\mathfrak{D}(\mathbf{r}-{\mathbf{r}}^{\prime},{\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2})=2.91{k}^{2}{\int}_{0}^{L}\mathrm{d}z{C}_{n}^{2}(z){|(\mathbf{r}-{\mathbf{r}}^{\prime})(1-\frac{z}{L})+z({\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2})|}^{5/3},$$## (7)

$$\u27e8\mathbf{\alpha}({\mathbf{\theta}}_{1})\xb7\mathbf{\alpha}({\mathbf{\theta}}_{2})\u27e9=(-\frac{2.91}{2}){\left(\frac{64}{\pi {D}^{4}}\right)}^{2}{\int}_{0}^{L}\mathrm{d}z{C}_{n}^{2}(z)\iint \mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime}W(\mathbf{r})W({\mathbf{r}}^{\prime})\mathbf{r}\xb7{\mathbf{r}}^{\prime}\phantom{\rule{0ex}{0ex}}\times {|(\mathbf{r}-{\mathbf{r}}^{\prime})(1-\frac{z}{L})+z({\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2})|}^{5/3}.$$^{15}and Winick and Marquis.

^{16}Applying those techniques, Eq. (7) reduces to

## (8)

$$\u27e8\mathbf{\alpha}({\mathbf{\theta}}_{1})\xb7\mathbf{\alpha}({\mathbf{\theta}}_{2})\u27e9=(-\frac{2.91}{2}){\left(\frac{16}{\pi}\right)}^{2}{D}^{-1/3}{\int}_{0}^{L}\mathrm{d}z{C}_{n}^{2}(z)\phantom{\rule{0ex}{0ex}}\times {\int}_{0}^{2\pi}\mathrm{d}\vartheta {\int}_{0}^{1}\mathrm{d}u[(u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]\phantom{\rule{0ex}{0ex}}\times {[{u}^{2}{(1-\frac{z}{L})}^{2}+{(\frac{z}{D}|{\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2}|)}^{2}+2u(1-\frac{z}{L}\left)\right(\frac{z}{D}|{\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2}|\left)\mathrm{cos}\text{\hspace{0.17em}}\vartheta \right]}^{5/6}.$$## (9)

$${\mathbf{\alpha}}_{p}=\frac{\int \mathrm{d}\mathbf{\theta}\mathbf{\alpha}(\mathbf{\theta}){P}_{G}(\mathbf{\theta})}{\int \mathrm{d}\mathbf{\theta}{P}_{G}(\mathbf{\theta})},$$The correlation between two patch-averaged tilts ${\mathbf{\alpha}}_{p1}$ and ${\mathbf{\alpha}}_{p2}$ is

## (10)

$$\u27e8{\mathbf{\alpha}}_{p1}.{\mathbf{\alpha}}_{p2}\u27e9=\frac{\iint \mathrm{d}{\mathbf{\theta}}_{1}\mathrm{d}{\mathbf{\theta}}_{2}\u27e8\mathbf{\alpha}({\mathbf{\theta}}_{1})\xb7\mathbf{\alpha}({\mathbf{\theta}}_{2})\u27e9{P}_{G}({\mathbf{\theta}}_{1}){P}_{G}({\mathbf{\theta}}_{2}-\mathrm{\Delta}\mathbf{\theta})}{\iint \mathrm{d}{\mathbf{\theta}}_{1}\mathrm{d}{\mathbf{\theta}}_{2}{P}_{G}({\mathbf{\theta}}_{1}){P}_{G}({\mathbf{\theta}}_{2}-\mathrm{\Delta}\mathbf{\theta})},$$## (11)

$$\u27e8{\mathbf{\alpha}}_{p1}.{\mathbf{\alpha}}_{p2}\u27e9=(-\frac{2.91}{2}){\left(\frac{16}{\pi {A}_{p}}\right)}^{2}{D}^{-1/3}{\int}_{0}^{L}\mathrm{d}z{C}_{n}^{2}(z)\iint \mathrm{d}{\mathbf{\theta}}_{1}\mathrm{d}{\mathbf{\theta}}_{2}{P}_{G}({\mathbf{\theta}}_{1}){P}_{G}({\mathbf{\theta}}_{2}-\mathrm{\Delta}\mathbf{\theta})\phantom{\rule{0ex}{0ex}}\times {\int}_{0}^{2\pi}\mathrm{d}\vartheta {\int}_{0}^{1}\mathrm{d}u[(u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]\phantom{\rule{0ex}{0ex}}\times {[{u}^{2}{(1-\frac{z}{L})}^{2}+{(\frac{z}{D}|{\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2}|)}^{2}+2u(1-\frac{z}{L}\left)\right(\frac{z}{D}|{\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2}|\left)\mathrm{cos}\text{\hspace{0.17em}}\vartheta \right]}^{5/6},$$Equation (11) can be equivalently expressed as

## (12)

$$\u27e8{\mathbf{\alpha}}_{p1}.{\mathbf{\alpha}}_{p2}\u27e9=(-\frac{2.91}{2}){\left(\frac{16}{\pi {A}_{p}}\right)}^{2}{D}^{-1/3}{\int}_{0}^{L}\mathrm{d}z{C}_{n}^{2}(z)\iint \mathrm{d}{\mathbf{\theta}}_{1}\mathrm{d}{\mathbf{\theta}}_{2}^{\prime}{P}_{G}({\mathbf{\theta}}_{1}){P}_{G}({\mathbf{\theta}}_{2}^{\prime})\phantom{\rule{0ex}{0ex}}\times {\int}_{0}^{2\pi}\mathrm{d}\vartheta {\int}_{0}^{1}\mathrm{d}u[(u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]\phantom{\rule{0ex}{0ex}}\times {[{u}^{2}{(1-\frac{z}{L})}^{2}+{(\frac{z}{D}|{\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2}^{\prime}-\mathrm{\Delta}\mathbf{\theta}|)}^{2}+2u(1-\frac{z}{L}\left)\right(\frac{z}{D}|{\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2}^{\prime}-\mathrm{\Delta}\mathbf{\theta}|\left)\mathrm{cos}\text{\hspace{0.17em}}\vartheta \right]}^{5/6}.$$Let

## (13)

$${\mathbf{\theta}}_{1}-{\mathbf{\theta}}_{2}^{\prime}=\frac{d}{L}\mathbf{x},\phantom{\rule{0ex}{0ex}}{\mathbf{\theta}}_{1}+{\mathbf{\theta}}_{2}^{\prime}=\frac{2d}{L}\mathbf{y}.$$## (14)

$$\u27e8{\mathbf{\alpha}}_{p1}.{\mathbf{\alpha}}_{p2}\u27e9=(-\frac{2.91}{2\pi}){\left(\frac{16}{\pi}\right)}^{3}{D}^{-1/3}{\int}_{0}^{L}\mathrm{d}z{C}_{n}^{2}(z)\iint \mathrm{d}\mathbf{x}\mathrm{d}\mathbf{y}[{e}^{-{(\mathbf{x}+2\mathbf{y})}^{2}}{e}^{-{(2\mathbf{y}-\mathbf{x})}^{2}}]\phantom{\rule{0ex}{0ex}}\times {\int}_{0}^{2\pi}\mathrm{d}\vartheta {\int}_{0}^{1}\mathrm{d}u[(u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]\phantom{\rule{0ex}{0ex}}\times {[{u}^{2}{(1-\frac{z}{L})}^{2}+{\left(\frac{zd}{DL}\right|\mathbf{x}-\frac{L}{d}\text{\hspace{0.17em}}\mathrm{\Delta}\mathbf{\theta}\left|\right)}^{2}+2u(1-\frac{z}{L}\left)\right(\frac{zd}{DL}|\mathbf{x}-\frac{L}{d}\text{\hspace{0.17em}}\mathrm{\Delta}\mathbf{\theta}|\left)\mathrm{cos}\text{\hspace{0.17em}}\vartheta \right]}^{5/6}.$$Substituting the above result in Eq. (14), the expression for patch-averaged tilt correlation becomes

## (15)

$$\u27e8{\mathbf{\alpha}}_{p1}.{\mathbf{\alpha}}_{p2}\u27e9=-\left(\frac{2.91}{\pi}\right){\left(\frac{16}{\pi}\right)}^{2}{D}^{-1/3}{\int}_{0}^{L}\mathrm{d}z{C}_{n}^{2}(z){\int}_{0}^{2\pi}\mathrm{d}\delta {\int}_{0}^{\infty}\mathrm{d}x(x{e}^{-2{x}^{2}})\phantom{\rule{0ex}{0ex}}\times {\int}_{0}^{2\pi}\mathrm{d}\vartheta {\int}_{0}^{1}\mathrm{d}u[(u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]\phantom{\rule{0ex}{0ex}}\times {\{{u}^{2}{(1-\frac{z}{L})}^{2}+{\left(\frac{zd}{DL}\right)}^{2}[{x}^{2}+{\left(\frac{L}{d}\mathrm{\Delta}\theta \right)}^{2}-2x\frac{L}{d}\mathrm{\Delta}\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\delta ]\phantom{\rule{0ex}{0ex}}+2u(1-\frac{z}{L}\left)\right(\frac{zd}{DL}\left)\mathrm{cos}\text{\hspace{0.17em}}\vartheta \sqrt{{x}^{2}+{\left(\frac{L}{d}\mathrm{\Delta}\theta \right)}^{2}-2x\frac{L}{d}\mathrm{\Delta}\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\delta}\right\}}^{5/6}.$$## (16)

$$\u27e8{({\mathbf{\alpha}}_{p1}-{\mathbf{\alpha}}_{p2})}^{2}\u27e9=2[\u27e8{\mathbf{\alpha}}_{p}^{2}\u27e9-\u27e8{\mathbf{\alpha}}_{p1}.{\mathbf{\alpha}}_{p2}\u27e9]\phantom{\rule{0ex}{0ex}}={\int}_{0}^{L}\mathrm{d}z{C}_{n}^{2}(z){f}_{d}(z),$$## (17)

$${f}_{d}(z)=-11.64{\left(\frac{16}{\pi}\right)}^{2}{D}^{-1/3}{\int}_{0}^{\infty}\mathrm{d}x(x{e}^{-2{x}^{2}})\phantom{\rule{0ex}{0ex}}\times {\int}_{0}^{2\pi}\mathrm{d}\vartheta {\int}_{0}^{1}\mathrm{d}u[(u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]\phantom{\rule{0ex}{0ex}}\times {[{u}^{2}{(1-\frac{z}{L})}^{2}+{\left(\frac{zd}{DL}x\right)}^{2}+2u(1-\frac{z}{L}\left)\right(\frac{zd}{DL}x\left)\mathrm{cos}\text{\hspace{0.17em}}\vartheta \right]}^{5/6}\phantom{\rule{0ex}{0ex}}+\left(\frac{5.82}{\pi}\right){\left(\frac{16}{\pi}\right)}^{2}{D}^{-1/3}{\int}_{0}^{2\pi}\mathrm{d}\delta {\int}_{0}^{\infty}\mathrm{d}x(x{e}^{-2{x}^{2}})\phantom{\rule{0ex}{0ex}}\times {\int}_{0}^{2\pi}\mathrm{d}\vartheta {\int}_{0}^{1}\mathrm{d}u[(u\text{\hspace{0.17em}}{\mathrm{cos}}^{-1}\text{\hspace{0.17em}}u)-{u}^{2}(3-2{u}^{2})\sqrt{1-{u}^{2}}]\phantom{\rule{0ex}{0ex}}\times {\{{u}^{2}{(1-\frac{z}{L})}^{2}+{\left(\frac{zd}{DL}\right)}^{2}[{x}^{2}+{\left(\frac{L}{d}\mathrm{\Delta}\theta \right)}^{2}-2x\frac{L}{d}\mathrm{\Delta}\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\delta ]\phantom{\rule{0ex}{0ex}}+2u(1-\frac{z}{L}\left)\right(\frac{zd}{DL}\left)\mathrm{cos}\text{\hspace{0.17em}}\vartheta \sqrt{{x}^{2}+{\left(\frac{L}{d}\mathrm{\Delta}\theta \right)}^{2}-2x\frac{L}{d}\mathrm{\Delta}\theta \text{\hspace{0.17em}}\mathrm{cos}\text{\hspace{0.17em}}\delta}\right\}}^{5/6},$$## 2.3.

### 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}_{n}^{2}$. 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}_{n}^{2}$ 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}

## (18)

$${r}_{0,sw}={[0.423{k}^{2}{\int}_{0}^{L}\mathrm{d}z{(1-\frac{z}{L})}^{5/3}{C}_{n}^{2}(z)]}^{-3/5}.$$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}_{n}^{2}$ estimates that compare directly to scintillometer measurements.

## 3.

## Results

## 3.1.

### 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}_{n}^{2}$. In essence, the path-weighting function for the ${C}_{n}^{2}$ estimate was the weighting function shown in Fig. 7(b). The ${C}_{n}^{2}$ estimates were compared to the BLS 2000 scintillometer measurements. Figures 8Fig. 9Fig. 10–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}_{n}^{2}$, 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}_{n}^{2}$ 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}_{n}^{2}$ 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.

## 3.2.

### Runway Imaging Experiment: Comparison of Time-Lapse Estimates to Scintillometer Measurements

Figure 12 shows a set of path weighting functions for the runway experiment for different patch separations. Two independent estimates of ${C}_{n}^{2}$ were obtained from the time-lapse images of July 18 and July 20. For July 18, four patches of $1/e$ diameter of 46 cm (40 pixels) centered at the four corners of the trailer window were used for the estimation. Hence there were two pairs of patches separated by 1.17 m and two pairs of patches separated by 2.44 m in each image. A 36-point moving window (corresponding to 3 min of imagery) was used to compute the variance. Since the path was nearly horizontal over the same surface type, it was presumed that ${C}_{n}^{2}$ was constant along the path. The path-averaged ${C}_{n}^{2}$ was obtained by dividing the total variance by the area of the corresponding weighting function 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}_{n}^{2}$ 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}_{n}^{2}$ 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.

## 4.

## 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 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}_{n}^{2}$ 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}_{n}^{2}$ 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}_{n}^{2}$.

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 time-lapse 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}_{n}^{2}$ 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.

## Acknowledgments

This work was supported in part by grants from the Directed Energy Joint Transition Office (DE-JTO) and the Air Force Office of Scientific Research (AFOSR). The authors thank Dr. Mikhail Vorontsov and his research group from UD for allowing them to set up the imaging experiment inside the Intelligent Optics Laboratory and for providing the BLS 2000 scintillometer data. The views expressed in this paper are those of the authors and do not necessarily reflect the official policy or position of the Air Force, the Department of Defense or the U.S. Government.

## References

## Biography

**Santasri R. Bose-Pillai** received her BSEE degree with honors from Jadavpur University, India, in 2000, her MSEE degree in 2005, and her PhD in electrical engineering with a focus in optics in 2008, both from New Mexico State University. Currently, she is a research assistant professor at AFIT’s Center for Directed Energy within Engineering Physics Department. Her research interests include propagation and imaging through atmosphere, telescope pointing and tracking, and synthesis of partially coherent sources. She is a senior member of SPIE and a regular member of OSA and DEPS.

**Jack E. McCrae** received his BS degree in physics from Massachusetts Institute of Technology in 1984, his MS degree in physics (optics), and his PhD in physics from the Air Force Institute of Technology in 1993 and 1997, respectively. 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. His research interests include optics, lasers, quantum and nonlinear optics, laser radar, atmospheric propagation, and imaging. He is a member of SPIE, OSA, and DEPS.

**Christopher A. Rice** received his BS degree in electrical engineering from Cedarville University, his MS degree in electrical engineering from the Air Force Institute of Technology (AFIT), and his PhD in optical sciences and engineering also from AFIT. Currently, he serves as a research assistant professor in the Department of Engineering Physics at the AFIT. He is interested in topic areas related to high-energy lasers, remote sensing, laser lethality, and optical diagnostics. His work on specific research topics currently include atmospheric propagation of high-energy lasers; diode pumped alkali and rare gas laser gain cell construction; innovative laser demonstration; aerosol field measurements; modeling, simulation, and validation of directed energy simulations; pulsed laser ablation of titanium and carbon; and new techniques of turbulence characterization over open paths.

**Ryan A. Wood** is a third-year undergraduate studying statistics in the Faculty of Arts and Sciences at Harvard University. He worked at AFIT’s Center for Directed Energy over the summer of 2017. His interests include computer vision, statistical learning, and data mining.

**Connor E. Murphy** graduated from high school in 2015 and has completed his degree in physics and mathematics double major at Grove City College. He has worked as an intern at the Air Force Institute of Technology’s Center for Directed Energy in the summer of 2017 and at the Air Force Research Laboratory’s Directed Energy Department in the summer of 2018. His research interests have included characterization and modeling of atmospheric turbulence, nanotechnology research at Grove City College, and the development and design of aluminum-air fuel cells. He is a member of the Society of Physics Students (SPS) and Sigma Pi Sigma, the physics honor society.

**Steven T. Fiorino** received his BS degrees in geography and meteorology from Ohio State and Florida State universities in 1987 and 1989, respectively. He additionally holds his 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 with 21 years of service and currently an associate professor of atmospheric physics within the Engineering Physics Department at AFIT and is the director of the Center for Directed Energy. His research interests include microwave remote sensing, development of weather signal processing algorithms, and atmospheric effects on military systems such as high-energy lasers and weapons of mass destruction. He is a member of SPIE, OSA, AMS, AIAA, and DEPS.