Estimation of Atmospheric Turbulence Parameters using Differential Motion of Extended Features in Time-lapse Imagery

Abstract

A system and method provide improved remote turbulence measurement. The system includes an image capturing device that captures time-lapse images of a distant target anywhere from a km to more than a 100 km away. A processor of the system tracks relative motion due to atmospheric turbulence of some number of patches of definite size on each of these images using a subpixel accurate correlation technique. The processor computes differential tilt variances between every pair of patches from the image collection and evaluates the theoretical weighting functions between turbulence along the path and differential tilt variances. The processor determines weights to linearly combine the weighting functions such that the combined weighting function closely resembles the weighting function corresponding to a turbulence parameter of interest. The processor then combines the differential tilt variances using the determined weights to obtain the desired turbulence parameter.

Description

BACKGROUND

1. Technical Field

The present disclosure generally relates to optical or near infra-red air turbulence measurements system.

2. Description of the Related Art

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. This analysis is relevant to laser communications, surveillance or tactical applications. Irradiance based methods such as the Scintillation, Detection and Ranging (SCIDAR) 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² 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² that are derived from such intermediate quantities as the velocity structure constant, C_v² and the temperature structure constant, C_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.

SUMMARY

The present innovation overcomes the foregoing problems and other shortcomings, drawbacks, and challenges of measuring turbulence for surveillance systems that are affected by turbulence or a tactical weapon engagements scenario where accuracy is dependent on turbulence monitoring. While the present innovation will be described in connection with certain embodiments, it will be understood that the invention is not limited to these embodiments. To the contrary, this invention includes all alternatives, modifications, and equivalents as may be included within the spirit and scope of the present invention.

Scintillometers are gold standards for measuring turbulence. However, they have limited operational ranges since they suffer from saturation issues. They also require access to both ends of the path. Other turbulence measurement techniques such as the Differential Image Motion Monitor (DIMM), Scintillation, Detection and Ranging (SCIDAR), Slope, Detection and Ranging (SLODAR) and MZA Associates' Path Resolved Optical Profiler System (PROPS) require active sources at the target location. Some of these instruments require sophisticated instrumentation as well. The proposed technique is a low cost approach to remotely characterize turbulence from a single site without deployment of sensors or sources at the target location. MZA's Delayed Tilt Anisoplanatism (DELTA) is an imaging based technique that uses the turbulence induced motion of point features on a distant target to remotely profile turbulence. However, since point features are tracked, the operational range is limited. Our imaging technique is phase-based (hence no saturation) and it tracks motion of extended patches rather than point features which allow the technique to be applied over longer ranges. The measured tilt variances can be linearly combined in such a way that the output can be made to mimic the output of any traditional turbulence measuring system as well. [5-8]

According to one aspect of the present innovation, a turbulence characterization system includes (i) An image capturing device that captures time-lapse images of a distant target (anywhere from a km to more than a 100 km away), and (ii) A turbulence measurement system, comprising of a processor, that is communicatively connected to the image capturing device to receive images. The processor tracks the relative motion due to atmospheric turbulence between some number of patches of definite size on each of these images using a sub-pixel accurate correlation technique, the processor computes the differential motion between all pairs of patches on each image, this differential motion is simply proportional to the differential wavefront tilt. The processor computes the differential tilt variance for each pair of patches from the collection of images and evaluates the theoretical weighting functions between turbulence along the imaging path and the measured differential tilt variances. These weighting functions depend upon the size and separation of the patches being tracked as well as the path length and the size of the camera aperture. The processor determines a linear combination of these weighting functions which best matches the weighting function for some desired turbulence parameter (such as the isoplanatic patch size, the log-amplitude variance, Fried's coherence diameter or the turbulence strength in a neighborhood around some specific range) and finally, linearly combines the measured tilt variances to estimate the turbulence parameter of interest.

According to another aspect of the present innovation, a method includes capturing frames of images of a distant target (anywhere from a km to more than a 100 km away) using an image capturing device. The method includes tracking the relative motion due to atmospheric turbulence between some number of patches of definite size on each of these images using a sub-pixel accurate correlation technique. The method includes computing differential motion between all pairs of patches on each image, this differential motion is simply proportional to the differential wavefront tilt. The method includes computing differential tilt variances for each pair of patches from the collection of frames. The method includes computing theoretical weighting functions between turbulence along the imaging path and the measured differential tilt variances. The method includes determining weights to linearly combine weighting functions such that the combined weighting function closely resembles a desired weighting function that corresponds to the turbulence parameter of interest (such as the isoplanatic patch size, the log-amplitude variance, Fried's coherence diameter or the turbulence strength in a neighborhood around some specific range). The method includes linearly combining the differential tilt variances using the determined weights to obtain the turbulence parameter of interest.

Additional objects, advantages, and novel features of the invention will be set forth in part in the description which follows, and in part will become apparent to those skilled in the art upon examination of the following or may be learned by practice of the invention. The objects and advantages of the invention may be realized and attained by means of the instrumentalities and combinations particularly pointed out in the appended claims.

BRIEF DESCRIPTION OF THE DRAWINGS

The description of the illustrative embodiments can be read in conjunction with the accompanying figures. It will be appreciated that for simplicity and clarity of illustration, elements illustrated in the figures have not necessarily been drawn to scale. For example, the dimensions of some of the elements are exaggerated relative to other elements. Embodiments incorporating teachings of the present disclosure are shown and described with respect to the figures presented herein, in which:

FIG. 1A depicts a diagram of an example ground-based experimental imaging setup showing the imaging system and the target, according to one or more embodiments;

FIG. 1B is a diagram of an example airborne imaging setup including an airborne imaging system and a ground target, according to one or more embodiments;

FIG. 2 is a graphical plot of path weighting functions of differential patch-averaged tilt variances for different patch sizes and separations, according to one or more embodiments;

FIG. 3A is a graphical plot of three path weighting functions, according to one or more embodiments;

FIG. 3B is a graphical plot of the linear combination of the three path weighting functions of FIG. 3A, according to one or more embodiments;

FIG. 4A is a graphical plot of two path weighting functions, according to one or more embodiments;

FIG. 4B is a graphical plot of linear combination of the two path weighting functions of FIG. 4A, according to one or more embodiments;

FIG. 5A is a map of an example imaging path;

FIG. 5B is an elevation profile of the imaging path;

FIG. 6A is an example image of the Dayton VA Medical Center captured early morning;

FIG. 6B is an example image of the Dayton VA Medical Center captured during the afternoon with patches tracked marked as white circles in the lower right of FIG. 6A;

FIG. 7 is a graphical plot of a comparison of C_n² estimates from time-lapse imagery with C_n² from scintillometer, according to one or more embodiments;

FIG. 8 is a graphical plot of a comparison of C_n² estimates from time-lapse imagery with C_n² from scintillometer, according to one or more embodiments;

FIG. 9 is a graphical plot of a comparison of C_n² estimates from time-lapse imagery with C_n² from scintillometer, according to one or more embodiments;

FIG. 10 is a graphical plot of a comparison of C_n² estimates from time-lapse imagery with C_n² from scintillometer, according to one or more embodiments;

FIG. 12A is a functional block diagram of an example system that uses a turbulence measurement system that enables turbulence compensation in an imaging surveillance system; according to one or more embodiments;

FIG. 12B is a functional block diagram of an example system that uses a turbulence measurement system along with a nonlinear optimizer to provide turbulence profiles, according to one or more embodiments;

FIG. 12C is a functional block diagram of an example system that produces refractivity/refractive index gradient information, according to one or more embodiments; and

FIG. 13 presents a flow diagram of a method for remote measurement of an arbitrary turbulence parameter from time-lapse imagery, according to one or more embodiments.

DETAILED DESCRIPTION

According to aspects of the present disclosure, a system and method provide improved remote turbulence measurement. The system includes an image capturing device that captures time-lapse images of a distant target (anywhere from a km to more than a 100 km away). A processor of the system tracks relative motion due to atmospheric turbulence of some number of patches of definite size on each of these images using a subpixel accurate correlation technique. The processor computes differential tilt variances between every pair of patches from the image collection and evaluates the theoretical weighting functions between turbulence along the path and differential tilt variances. The processor determines weights to linearly combine the weighting functions such that the combined weighting function closely resembles the weighting function corresponding to a turbulence parameter of interest. The processor then combines the differential tilt variances using the determined weights to obtain the desired turbulence parameter.

In one aspect of the present disclosure, a method includes obtaining a sequence of images of a distant target (anywhere from a km to more than a 100 km away) using a digital image capture device, or even digitizing film images. The method includes tracking the relative motion due to atmospheric turbulence between some number of patches of definite size on each of these images using a sub-pixel accurate correlation technique. The method includes computing the differential motion between all pairs of patches on each image, this differential motion is simply proportional to the differential wavefront tilt. The method includes computing the differential tilt variance for each pair of patches from the collection of images. The method includes evaluating the theoretical weighting functions between turbulence along the imaging path and the measured differential tilt variances. These weighting functions depend upon the size and separation of the patches being tracked as well as the path length and the size of the camera aperture. An expression for these theoretical weighting functions is derived in the description section under the presumptions of geometric optics and turbulence which follows the Kolmogorov power spectrum. The method includes determining a linear combination of these weighting functions which best matches the weighting function for some desired turbulence parameter (such as the isoplanatic patch size, Fried's coherence diameter, the log-amplitude variance, or the turbulence strength in a neighborhood around some specific range).

FIG. 1A depicts a diagram of an example ground-based experimental imaging setup 100a showing an imaging system 102a and a target 104a. FIG. 1B is a diagram of an example airborne imaging setup 100b including an airborne imaging system 102b and a ground target 104b. In one or more embodiments, an imaging system, comprising of an electronic camera, fitted with an imaging lens is used to capture several short exposure images of a distant target. The target should have sufficient trackable features on it. The imaging system is mounted on a tripod. The exposure time is selected such that the atmosphere is frozen during a single capture, yet sufficient signal is present in the images. The time interval between frames is selected such that the turbulence has changed between frames, and the tilts in adjacent frames are statistically independent.

The images are first cropped to isolate the region of interest. A reference image is chosen from the images collected and a cross-correlation algorithm is used to estimate the image shift between each image and the reference image. A Gaussian window is applied to the images before the cross-correlation algorithm is run. This reduces the effects of the frame edges on the correlation result. A parabolic fit is applied to the correlation peak to provide a sub-pixel estimate of the shift between the images. The correlation with a reference image provides information about the slow motion during the course of the day due to refractive bending. This long-term drift information is used to adjust the tracking window keeping it locked on to a feature through the collection. The cross-correlation algorithm is now applied to each image and its neighboring image to estimate the random motion of the feature due to turbulence. The shifts of two features are subtracted to get the differential signal. The differential mode of measurement eliminates the effects of common mode disturbances, such as platform vibrations.

B.1 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 [9]:

$\begin{array}{cc}a\left(\theta \right)=\frac{32\lambda }{{\pi }^2{D}^4}\int drW\left(r\right)r\varphi \left(r,\theta \right),& \left(1\right)\end{array}$

where λ is the wavelength, ϕ(r,θ) is the turbulence induced wavefront distortion at aperture coordinate r and

$\begin{array}{cc}W\left(r\right)=\{\begin{array}{cc}1& r\le 0.5D\\ 0& r>0.5D\end{array}.& \left(2\right)\end{array}$

The mean correlation between tilts observed at the aperture due to two sources at viewing directions θ₁ and θ₂ can be written as,

$\begin{array}{cc}\langle \alpha \left({\theta }_1\right)·\alpha \left({\theta }_2\right)\rangle ={\left(\frac{32\lambda }{{\pi }^2{D}^4}\right)}^2\langle \int \int \mathrm{dr}{\mathrm{dr}}^{\prime }W\left(r\right)W\left(r^{\prime }\right)r·r^{\prime }\varphi \left(r,{\theta }_1\right)\varphi \left(r^{\prime },{\theta }_2\right)\rangle ,& \left(3\right)\end{array}$

where the angled brackets indicate ensemble averaging. Interchanging the order of integration and ensemble averaging results in,

$\begin{array}{cc}\langle \alpha \left({\theta }_1\right)·\alpha \left({\theta }_2\right)\rangle ={\left(\frac{32\lambda }{{\pi }^2{D}^4}\right)}^2\int \int \mathrm{dr}{\mathrm{dr}}^{\prime }W\left(r\right)W\left(r^{\prime }\right)r·r^{\prime }\langle \varphi \left(r,{\theta }_1\right)\varphi \left(r\text{'},{\theta }_2\right)\rangle .& \left(4\right)\end{array}$

Since ∫∫drdr′W(r)W(r′)r·r′=0 terms which are functions of either r or r′, and not both, can be added without changing the result of the integration.

Hence, Eq. (4) becomes,

$\begin{array}{cc}\langle \alpha \left({\theta }_1\right)·\alpha \left({\theta }_2\right)\rangle =\left(-\frac{1}{2}\right){\left(\frac{32\lambda }{{\pi }^2{D}^4}\right)}^2\int \int \mathrm{dr}{\mathrm{dr}}^{\prime }W\left(r\right)W\left(r^{\prime }\right)r·r^{\prime }{D}_{\varphi }\left(r-r^{\prime },{\theta }_1-{\theta }_2\right),& \left(5\right)\end{array}$

where D_ϕ(r−r′,θ₁−θ₂)=[ϕ(r,θ₁)−ϕ(r′,θ₂)]² 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, [10]

$\begin{array}{cc}D\left(r-r^{\prime },{\theta }_1-{\theta }_2\right)=2.91k^2\underset{0}{\overset{L}{\int }}dz{C}_n^2\left(z\right){\left(r-r^{\prime }\right)\left(1-\frac{z}{L}\right)+z\left({\theta }_1-{\theta }_2\right)}^{5/3},& \left(6\right)\end{array}$

where k=2π/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) in Eq. (5),

$\begin{array}{cc}\langle \alpha \left({\theta }_1\right)·\alpha \left({\theta }_2\right)\rangle =\left(-\frac{2.91}{2}\right){\left(\frac{64}{\pi D^4}\right)}^2\underset{0}{\overset{L}{\int }}dz{C}_n^2\left(z\right)\int \int \mathrm{dr}{\mathrm{dr}}^{\prime }W\left(r\right)W\left(r^{\prime }\right)r·r^{\prime }\times {\left(r-r^{\prime }\right)\left(1-\frac{z}{L}\right)+z\left({\theta }_1-{\theta }_2\right)}^{5/3}.& \left(7\right)\end{array}$

The integrations over r or r′ can be done using techniques outlined by Fried [11] and Winick et. al [12]. Applying those techniques, Eq. (7) reduces to

$\begin{array}{cc}\langle \alpha \left({\theta }_1\right)·\alpha \left({\theta }_2\right)\rangle =\left(-\frac{2.91}{2}\right){\left(\frac{16}{\pi }\right)}^2{D}^{-1/3}\underset{0}{\overset{L}{\int }}dz{C}_n^2\left(z\right)\times \underset{0}{\overset{2\pi }{\int }}\mathrm{d\vartheta }\underset{0}{\overset{1}{\int }}\mathrm{du}\left[\left(u{\cos }^{-1}u\right)-u^2\left(3-2u^2\right)\sqrt{1-u^2}\right]\times \left[{u^2\left(1-\frac{z}{L}\right)}^2+{\left(\frac{z}{D}\left|{\theta }_1-{\theta }_2\right|\right)}^2+2u\left(1-\frac{z}{L}\right)\left(\frac{z}{D}\left|{\theta }_1-{\theta }_2\right|\right)\cos \vartheta \right]\text{5/6}.& \left(8\right)\end{array}$

Regarding the novelty of the weighting functions, we have derived expressions for patch averaged tilt variances and the corresponding weighting functions while Whiteley et al. have expressions for weighting functions corresponding to a point source case. A patch is essentially made up of a bunch of incoherent point sources. The fact that we use the weighting functions for patch averaged tilt variance rather than tilt variance due to a point source allows us to apply this technique over long ranges, where even a pixel in the image is a patch and not a point. Tracking motion of extended patches rather than points on a target (over long ranges or in strong turbulence conditions, it is not possible to track points on a non-cooperative target) and developing weighting functions for the patch-averaged tilt variance described below are novel.

Since each pixel in the time-lapse imagery corresponds to a patch of finite size, not just a point on the target, 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 u, is defined as,

$\begin{array}{cc}{\alpha }_p=\frac{\int \mathrm{d\theta }\alpha \left(\theta \right)P_G\left(\theta \right)}{\int \mathrm{d\theta }{P}_G\left(\theta \right)},\mathrm{where}{P}_G\left(\theta \right)=e^{-\frac{4{\theta }^2{L}^2}{d^2}},d\mathrm{being}\mathrm{the}1/e\mathrm{patch}\mathrm{diameter}\mathrm{and}\theta =\left|\theta \right|.& \left(9\right)\end{array}$

The correlation between two patch-averaged tilts α_p1 and α_p2 is

$\begin{array}{cc}\langle {\alpha }_{p1}·{\alpha }_{p2}\rangle =\frac{\int \int {\mathrm{d\theta }}_1{\mathrm{d\theta }}_2\langle \alpha \left({\theta }_1\right)·\alpha \left({\theta }_2\right)\rangle P_G\left({\theta }_1\right)P_G\left({\theta }_2-\mathrm{\Delta \theta }\right)}{\int \int {\mathrm{d\theta }}_1{\mathrm{d\theta }}_2{P}_G\left({\theta }_1\right)P_G\left({\theta }_2-\mathrm{\Delta \theta }\right)},& \left(10\right)\end{array}$

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) in Eq. (10),

$\begin{array}{cc}\langle {\alpha }_{p1}·{\alpha }_{p2}\rangle =\left(-\frac{2.91}{2}\right){\left(\frac{16}{A_p}\right)}^2{D}^{-1/3}\stackrel{L}{\underset{0}{\int }}{\mathrm{dzC}}_n^2\left(z\right)\int \int {\mathrm{d\theta }}_1{\mathrm{d\theta }}_2{P}_G\left({\theta }_1\right)P_G\left({\theta }_2-\mathrm{\Delta \theta }\right)\times \stackrel{2}{\underset{0}{\int }}\mathrm{d\vartheta }\stackrel{1}{\underset{0}{\int }}\mathrm{du}\left[\left(u{\cos }^{-1}u\right)-u^2\left(3-2u^2\right)\sqrt{1-u^2}\right]\times {\left[{u^2\left(1-\frac{z}{L}\right)}^2+{\left(\frac{z}{D}{\theta }_1-{\theta }_2\right)}^2+2u\left(1-\frac{z}{L}\right)\left(\frac{z}{D}{\theta }_1-{\theta }_2\right)\cos \vartheta \right]}^{5/6},\mathrm{where}{A}_p=\int \mathrm{d\theta }{P}_G\left(\theta \right)=\frac{\pi d^2}{4L^2}.& \left(11\right)\end{array}$

Eq. (11) can be equivalently expressed as,

$\begin{array}{cc}\langle {\alpha }_{p1}·{\alpha }_{p2}\rangle =\left(-\frac{2.91}{2}\right){\left(\frac{16}{\pi A_p}\right)}^2{D}^{-1/3}\underset{0}{\overset{L}{\int }}{\mathrm{dzC}}_n^2\left(z\right)\int {\mathrm{d\theta }}_1{{\mathrm{d\theta }}_2}^{\prime }{P}_G\left({\theta }_1\right)P_G\left({{\theta }_2}^{\prime }\right)\times \underset{0}{\overset{2\pi }{\int }}\mathrm{d\vartheta }\underset{0}{\overset{1}{\int }}du\left[\left(u{\cos }^{-1}u\right)-u^2\left(3-2u^2\right)\sqrt{1-u^2}\right]\times {\left[{u^2\left(1-\frac{z}{L}\right)}^2+\frac{z}{D}{\theta }_1-{{\theta }_2}^{\prime }-\mathrm{\Delta \theta }\right)}^2+2u\left(1-\frac{z}{L}\right)\left(\frac{z}{D}{\theta }_1-{{\theta }_2}^{\prime }-\mathrm{\Delta \theta }\right)\cos \vartheta ]\text{5/6}.& \left(12\right)\\ \mathrm{Let}& \\ {\theta }_1-{{\theta }_2}^{\prime }=\frac{d}{L}x,& \\ {\theta }_1+{{\theta }_2}^{\prime }=\frac{2d}{L}y.& \left(13\right)\end{array}$

Changing the variables of integration, Eq. (12) reduces to,

$\begin{array}{cc}\langle {\alpha }_{p1}·{\alpha }_{p2}\rangle =\left(-\frac{2.91}{2\pi }\right){\left(\frac{16}{\pi }\right)}^3{D}^{-1/3}\underset{0}{\overset{L}{\int }}dz{C}_n^2\left(z\right)\int \int \mathrm{dx}\mathrm{dy}\left[e^{-{\left(x+2y\right)}^2}{e}^{-{\left(2y-x\right)}^2}\right]\times \underset{0}{\overset{2\pi }{\int }}\mathrm{d\vartheta }\underset{0}{\overset{1}{\int }}du\left[\left(u{\cos }^{-1}u\right)-u^2\left(3-2u^2\right)\sqrt{1-u^2}\right]\times \hspace{1em}[{u^2\left(1-\frac{z}{L}\right)}^2+{\left(\frac{\mathrm{zd}}{\mathrm{DL}}x-\frac{L}{d}\mathrm{\Delta \theta }\right)}^2+2u\left(1-\frac{z}{L}\right){\left(\frac{\mathrm{zd}}{\mathrm{DL}}x-\frac{L}{d}\mathrm{\Delta \theta }\cos \vartheta \right]}^{5/6}.& \left(14\right)\\ \mathrm{Now},\int \mathrm{dy}{e}^{-{\left(x+2y\right)}^2}{e}^{-{\left(2y-x\right)}^2}=\int \mathrm{dy}{e}^{-8\left(y^2+\frac{x^2}{4}\right)}=\underset{0}{\overset{\infty }{\int }}\mathrm{dy}\underset{0}{\overset{2\pi }{\int }}d\beta \left[y{e}^{-8\left(y^2+\frac{x^2}{4}\right)}\right]=\frac{\pi }{8}{e}^{-2x^2},\mathrm{where}x=\left|x\right|\mathrm{and}y=\left|y\right|.& \end{array}$

Substituting the above result in Eq. (14), the expression for patch-averaged tilt correlation becomes,

$\begin{array}{cc}\langle {\alpha }_{p1}·{\alpha }_{p2}\rangle =-\left(\frac{2.91}{\pi }\right){\left(\frac{16}{\pi }\right)}^2{D}^{-1/3}\underset{0}{\overset{L}{\int }}dz{C}_n^2\left(z\right)\underset{0}{\overset{2\pi }{\int }}d\delta \underset{0}{\overset{\infty }{\int }}dx\left(x{e}^{-2x^2}\right)\times \underset{0}{\overset{2\pi }{\int }}d\vartheta \underset{0}{\overset{1}{\int }}du\left[\left(u{\cos }^{-1}u\right)-u^2\left(3-2u^2\right)\sqrt{1-u^2}\right]\times \hspace{1em}\left[{u^2\left(1-\frac{z}{L}\right)}^2+{\left(\frac{zd}{DL}\right)}^2\left\{x^2+{\left(\frac{L}{d}\Delta \theta \right)}^2-2x\frac{L}{d}\mathrm{\Delta \theta }\cos \delta \right\}+2u\left(1-\frac{z}{L}\right)\left(\frac{zd}{DL}\right)\cos \vartheta \sqrt{x^2+{\left(\frac{L}{d}\mathrm{\Delta \theta }\right)}^2-2x\frac{L}{d}\mathrm{\Delta \theta }\cos \delta }\right]\text{5/6}.& \left(15\right)\end{array}$

The expression for patch-averaged tilt variance α_p² can be obtained from Eq. (15) by setting Δθ=0. The differential patch-averaged tilt variance can hence be expressed as:

$\begin{array}{cc}\begin{array}{c}\langle {\left({\alpha }_{p1}-{\alpha }_{p2}\right)}^2\rangle =2\left[\langle {{\alpha }_p}^2\rangle -\langle {\alpha }_{p1}·{\alpha }_{p2}\rangle \right]\\ =\underset{0}{\overset{L}{\int }}dz{C}_n^2\left(z\right)f_d\left(z\right),\end{array}& \left(16\right)\\ \mathrm{where}& \\ f_d\left(z\right)=-11.64{\left(\frac{16}{}\right)}^2{D}^{-1/3}\stackrel{\infty }{\underset{0}{\int }}\mathrm{dx}\left({\mathrm{xe}}^{-2x^2}\right)=\times \underset{0}{\overset{2\pi }{\int }}\mathrm{d\vartheta }\underset{0}{\overset{1}{\int }}du\left[\left(u{\cos }^{-1}u\right)-u^2\left(3-2u^2\right)\sqrt{1-u^2}\right]\times {\left[{u^2\left(1-\frac{z}{L}\right)}^2+{\left(\frac{zd}{DL}x\right)}^2+2u\left(1-\frac{z}{L}\right)\left(\frac{zd}{DL}x\right)\cos \vartheta \right]}^{5/6}+\left(\frac{5.82}{\pi }\right){\left(\frac{16}{\pi }\right)}^2{D}^{-1/3}\underset{0}{\overset{2\pi }{\int }}d\delta \underset{0}{\overset{\infty }{\int }}\mathrm{dx}\left(x{e}^{-2x^2}\right)\times \underset{0}{\overset{2\pi }{\int }}\mathrm{d\vartheta }\underset{0}{\overset{1}{\int }}du\left[\left(u{\cos }^{-1}u\right)-u^2\left(3-2u^2\right)\sqrt{1-u^2}\right]\times \hspace{1em}[{u^2\left(1-\frac{z}{L}\right)}^2+{\left(\frac{zd}{DL}\right)}^2\left\{x^2+{\left(\frac{L}{d}\Delta \theta \right)}^2-2x\frac{L}{d}\mathrm{\Delta \theta }\cos \delta \right\}+2u\left(1-\frac{z}{L}\right)\left(\frac{zd}{DL}\right)\cos \vartheta \sqrt{x^2+{\left(\frac{L}{d}\Delta \theta \right)}^2-2x\frac{L}{d}\mathrm{\Delta \theta }\cos \delta }]\text{5/6}& (\text{17)}\end{array}$

is the path weighting function. No simpler form for the weighting function can be obtained, and hence in the present work, Eq. (17) is evaluated numerically for cases of interest.

FIG. 2 is a graphical plot 200 of 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. FIG. 2 shows the graphical plot 200 of fd(z) for different patch sizes and separations. 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 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 peak locations of these weighting functions depend on the patch size and separation relative to the imaging aperture. Sub-aperture sized patches and separations have weighting functions that peak towards the target end of the path. Larger patch sizes and separations have weighting functions that peak towards the camera end of the path.

B.2 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². 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² along the path. Members of a set of path weighting functions, from a variety of differently sized and separated image patches, each characterize 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 can be used to find the least-squares optimal way to achieve the desired weighting function in terms of the basis set available. Here it is shown that the weighting functions for differential patch-averaged tilt variances can be used to reproduce the weighting functions for the spherical wave r₀ (Fried's coherence diameter) and a scintillometer.

The spherical wave r₀ is given by: [13]

$\begin{array}{cc}{r}_{0,\mathrm{sw}}={\left[0.423k^2\underset{0}{\overset{L}{\int }}d{z\left(1-\frac{z}{L}\right)}^{5/3}{C}_n^2\left(z\right)\right]}^{-3/5}.& \left(18\right)\end{array}$

This equation has been written so the integration takes place from the camera to the object, i.e. the camera is at 0 and the target is at L as above. In this expression C_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 sampled weighting functions, and the columns then correspond to the range where these functions are sampled, the desired weights can be computed as W=M⁺F, 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.

FIG. 3A is a graphical plot 300a of three path weighting functions. FIG. 3B is a graphical plot 300b of the linear combination of the three path weighting functions of FIG. 3A. Three path weighting functions in FIG. 3A linearly combined to obtain a weighting function very close to the r₀ weighting function in FIG. 3B. The camera is at z=0 and the target is at z=7000 m.

The Moore-Penrose technique is used here to obtain the pseudo-inverse. FIG. 3 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 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₀ from the differential motion of these pairs of patches.

FIG. 4A is a graphical plot 400a of two path weighting functions. FIG. 4B is a graphical plot 400b of linear combination of the two path weighting functions of FIG. 4A. Two path weighting functions in FIG. 4A linearly combined to obtain a weighting function very close to the BLS 2000 scintillometer weighting function in FIG. 4B. The camera is at z=0 and the target is at z=7000 m.

FIGS. 4A-4B demonstrate how weighting functions corresponding to pairs of patches of two different separations can be linearly combined to approximate the Scintec BLS 2000 scintillometer weighting function. Though the two weighting functions do not match perfectly, the agreement is fairly good. The measured differential tilt variances are multiplied by the corresponding weights obtained from the pseudo-inverse technique and combined together to obtain an estimate of the desired turbulence parameter.

C. How the Invention can be used: The following example demonstrates how the method developed can be used to obtain an estimate of path-weighted C_n², such as that from a scintillometer. [14] The experiment was 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. FIG. 5A depicts the layout of the 7 km path from the VA Medical Center to UD. FIG. 5B depicts a ground elevation profile of 7 km path of FIG. 5A.

Example images are shown in FIGS. 6A-6B. 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 inch (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. 6A, 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. 6A.

The four patches tracked were each of 1/e diameter 16 cm (20 pixels). Hence there were two pairs of patches separated by 48 cm, and two pairs of patches separated by 1.8m at the target for each image. A 30 frame moving window, corresponding to twenty (20) minutes 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 pseudo-inverse technique mentioned in the previous section) and added together to obtain a path-weighted estimate of C_n². In essence, the path-weighting function for the C_n² estimate was the weighting function shown in FIG. 4B. The estimates were compared to the BLS 2000 scintillometer measurements.

FIGS. 7-10 show the results of these comparisons for four days over different weather conditions. FIG. 7 is a graphical plot of a comparison of C_n² estimates from time-lapse imagery with C_n² from scintillometer: Feb. 14, 2017. Differences in the two profiles are possibly due to estimation errors resulting from moving shadows and poor image contrast during early morning and late afternoon. FIG. 8 is a graphical plot of a comparison of C_n estimates from time-lapse imagery with C_n² from scintillometer: Feb. 15, 2017. Differences in the two profiles are possibly due to estimation errors resulting from moving shadows and poor image contrast during early morning and late afternoon. FIG. 9 is a graphical plot of comparison of C_n² estimates from time-lapse imagery with C_n² from scintillometer: Feb. 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. 10 is a graphical plot of comparison of C_n estimates from time-lapse imagery with C_n² from scintillometer: Feb. 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.

According to meteorological observations [15], Feb. 14, 2017 was a clear day. As evident from FIG. 7, 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. 6A) contributed to estimation error during these times. Feb. 15, 2017 was clear earlier during the day, but it became cloudy by mid-morning. Overcast conditions prevailed from around noon to 1:30 PM (local standard time), causing the dip in C_n², seen in FIG. 8. Again, there is excellent agreement between the time-lapse estimates and the scintillometer measurements. Feb. 18, 2017 was marked by alternating cloudy and clear conditions throughout the day. This resulted in multiple peaks and troughs in the profile, as seen in FIG. 9. 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² profile. Feb. 21, 2017 was a cloudy day, with clear conditions developing very late in the afternoon. The time-lapse estimates agree reasonably well with the scintillometer, as seen in FIG. 10, 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.

Alternatives: The derivation of the weighting functions given in section B contains a presumption that the source patch is well modeled by the Gaussian window applied to the images. This presumption can be eliminated by using a measurement of the light intensity within each source window. In cases where this Gaussian window doesn't correspond well with the actual source patch detail, this approach will improve the accuracy of the method. To accomplish this, the Gaussian patch function introduced in equation 9 is replaced by the known intensity distribution of the patch being tracked. The known intensity distribution could be taken from the image set used, higher resolution images, or model data. The integrations which follow are then computed numerically. An intermediate approach is to compute the width (i.e. standard deviation) of the intensity distribution within each windowed source patch on the x and y axes and use these as the widths of the Gaussian windows used for the integrations. When these widths differ on x and y the derivation which follows needs to be modified to account for this bivariate Gaussian. In the cases of corners, edges, or small high contrast features, the width on one or both axes will be smaller than the Gaussian window applied, and this correction will prevent the turbulence from being overestimated.

FIG. 11 depicts optical system 1100 that includes image capturing device 1102 that captures time-lapse images of distant target 1104 that is anywhere from a km to more than 100 km away. Optical system 1100 includes turbulence measurement system 1106, comprising processor 108 that is communicatively coupled to image capturing device (ICD) 1102 to receive images from ICD 1102. Processor 1108 comprises one or more hardware computing devices that execute program code to configure optical system 1100 to estimate atmospheric turbulence parameters using differential motion of extended features in time-lapse imagery. Processor 1108 tracks the relative motion due to atmospheric turbulence between some number of patches of definite size on each of these images using a sub-pixel accurate correlation technique. Processor 1108 computes the differential motion between all pairs of patches on each image, this differential motion is simply proportional to the differential wavefront tilt. Processor 1108 computes the differential tilt variance for each pair of patches from the collection of images. Processor 1108 evaluates the theoretical weighting functions between turbulence along the imaging path and the measured differential tilt variances. These weighting functions depend upon the size and separation of the patches being tracked as well as the path length and the size of the camera aperture. Processor 1108 determines a linear combination of these weighting functions which best matches the weighting function for some desired turbulence parameter (such as the isoplanatic patch size, the log-amplitude variance, Fried's coherence diameter or the turbulence strength in a neighborhood around some specific range). And finally, processor 1108 linearly combines the measured tilt variances to estimate the turbulence parameter of interest.

In one or more embodiments, ICD 1102 is a digital camera operating in the visible or near infrared. In one or more embodiments, ICD 1102 captures images no faster than an atmospheric coherence time, with individual exposures shorter than the atmospheric coherence time. In one or more embodiments, ICD 1102 turbulence measurement is used to a) compensate for turbulence degraded images in a surveillance system and b) to understand how turbulence is distributed along the path in a tactical engagement scenario.

FIG. 12A is a functional block diagram of an example system 1200 that uses a turbulence measurement system 1202 that enables turbulence compensation in a surveillance system 1204. Imaging system 1206 images target 1208 and provides turbulence-degraded images 1210 to the turbulence measurement system 1202 and provides turbulence-degraded images 1210 to the turbulence-compensated surveillance system 1204. Turbulence measurement system 1202 generates a turbulence parameter 1214 that is used by the turbulence-compensated surveillance system 1204 for compensation.

In one or more embodiments, the turbulence measurement system uses the pseudoinverse technique to obtain the weights required to linearly combine the differential tilt variance weighting functions to obtain a weighting function that approximates the weighting function of a desired turbulence parameter. The measured differential tilt variances from the sample images are linearly combined according to the computed weights to get the turbulence parameter of interest. The images are blurred and warped due to turbulence. The calculated turbulence parameter can be used to generate a model of the Optical Transfer Function which can be used to deconvolve the turbulence degraded images to get pristine images for surveillance.

FIG. 12B is a functional block diagram of an example system 1220 that uses a turbulence measurement system 1222 to provide weighting functions 1223 and tilt variances 1224 along with a nonlinear optimizer 1225 to provide turbulence profiles 1226. In one or more embodiments, the turbulence measurement system uses a pattern recognition algorithm to identify trackable features in the images. A set of weighting functions for differential tilt variances is numerically computed for all pairs of feature separations. The corresponding differential tilt variances are computed. The weighting functions are related to the tilt variances through a linear system of equations where the turbulence strength at every sampled location along the path is the unknown to be solved. The turbulence strengths at different locations along the path can be solved through a suitable optimization routine that is robust to measurement noise. The turbulence distribution along the path can be integrated to generate any desirable integrated turbulence parameter or simply used for turbulence monitoring during a High Energy Laser (HEL) field testing or during tactical engagements. It will provide important information about how turbulence at different locations is affecting the HEL performance. If used from an airborne platform, the profiling information can be useful to understand how turbulence varies with height. Information about this is limited currently and the knowledge can help improve existing atmospheric models.

FIG. 12C is a functional block diagram of an example system 1242 that uses a first processor 1244 to compute feature shifts 1246 and a second processor 1248 to produce refractivity/refractive index gradient information 1250. In one or more embodiments, the system uses a pattern recognition algorithm to compute the slow, deterministic motion of trackable features in the images over the course of several days due to refractive bending. Since curvature of a ray is related to refractive index gradient, this motion information can be used to obtain the temperature structure at a particular location, e.g. determine the lapse rate and how local weather conditions can influence it. Knowledge of these parameters will be used to improve existing atmospheric models. Refractivity information is useful for pointing/tracking in Directed Energy or Electro Optics Sensing applications.

FIG. 13 presents a flow diagram of a method for estimation of atmospheric turbulence parameters using differential motion of extended features in time-lapse imagery. Method 1300 includes obtaining a sequence of short exposure images of a distant target using a digital image capture device or by digitizing film images (block 1302). In particular, method 1300 obtains the images of the distant target that is anywhere from a km to more than a 100 km away. Method 1300 includes identifying some number of patches/features of definite size on these images and track their relative motion between image frames using a sub-pixel accurate correlation technique (block 1304). Method 1300 includes computing differential motion (or differential tilt) between all pairs of patches on each image (block 1306). Method 1300 includes computing differential tilt variance for each pair of patches from the collection of images (block 1308). In particular, method 1300 is applicable over long ranges where even a pixel is a patch and not a point. Method 1300 includes evaluating theoretical weighting functions between turbulence along the imaging path and measured differential tilt variances (block 1310). These weighting functions depend upon the size and separation of the patches being tracked as well as the path length and the size of the camera aperture. Method 1300 includes determining a linear combination of these weighting functions which best matches the weighting function for some desired turbulence parameter block 1312). Method 1300 includes linearly combining the measured differential tilt variances using the weights determined above to estimate the turbulence parameter of interest (block 1314). Then method 1300 ends.

In one or more aspects of the present disclosure, a method includes obtaining a sequence of images of a distant target (anywhere from a km to more than 100 km away) using a digital image capture device, or even digitizing film images. The method includes tracking the relative motion due to atmospheric turbulence between some number of patches of definite size on each of these images using a sub-pixel accurate correlation technique. The method includes computing the differential motion between all pairs of patches on each image, this differential motion is simply proportional to the differential wavefront tilt. The method includes computing the differential tilt variance for each pair of patches from the collection of images. The method evaluating the theoretical weighting functions between turbulence along the imaging path and the measured differential tilt variances. These weighting functions depend upon the size and separation of the patches being tracked as well as the path length and the size of the camera aperture. The method includes determining a linear combination of these weighting functions which best matches the weighting function for some desired turbulence parameter (such as the isoplanatic patch size, the log-amplitude variance, Fried's coherence diameter, or the turbulence strength in a neighborhood around some specific range). The method includes linearly combining the measured tilt variances to estimate the turbulence parameter of interest.

The following references cited above are hereby incorporated by reference in their entirety:

• [1] R. Avila, J. Vernin, and E. Masciadri, “Whole atmospheric-turbulence profiling with generalized scidar,” Appl. Opt. 36 (30), 7898-7905 (1997).
• [2] V. Kluckers, N. Wooder, T. Nicholls, M. Adcock, I. Munro, and J. Dainty, “Profiling of atmospheric turbulence strength and velocity using a generalized SCIDAR technique,” Astron. Astrophys. Suppl. Ser 130(1), 141-155 (1998).
• [3] R. J. Hill, “Comparison of scintillation methods for measuring the inner scale of turbulence,” Appl. Opt. 27, 2187-2193 (1988).
• [4] Z. Junwei, Q. Xiwen, F. Shuanglian, and Z. Zhigang, “Measurement of atmospheric refractive index structure constant near ground based on CAN bus,” in 10th International Conference on Electronic Measurements and Instruments (ICEMI), 2, 270-273 (2011).
• [5] M. Sarazin and F. Roddier, “The ESO Differential Image Motion Monitor,” Astron. Astrophys. 227, 294-300 (1990).
• [6] Butterley, T., Wilson, R. W., and Sarazin, M., “Determination of the profile of atmospheric optical turbulence strength from slodar data,” Mon. Not. R. Astron. Soc. 369, 835-845 (2006).
• [7] MZA Corporation, “DELTA Imaging path turbulence monitor PM02-600,” https://www.mza.com/doc/misc/MZA_DELTA_PM-02-600_Specifications.pdf
• [8] MZA Corporation, “Path-resolved optical profiler system (PROPS) PR-05-600,” https://www.mza.com/doc/misc/MZA_PROPS_PR-05-600_Specifications.pdf
• [9] D. Fried, “Differential angle of arrival: theory, evaluation, and measurement feasibility,” Rad. Sci. 10(1), 71-76 (1975).
• [10] M. Roggemann, and B. Welsh, [Imaging through Turbulence], CRC Press, Boca Raton, Fla. (1996).
• [11] D. Fried, “Varieties of isoplanatism,” in Imaging through the Atmosphere, J. C. Wyant, ed., Proc. SPIE 75, 20-29 (1976).
• [12] K. Winick, and D. Marquis, “Stellar scintillation technique for the measurement of tilt anisoplanatism,” J. Opt. Soc. Am. A 5(11), 1929-1936 (1988).
• [13] J. D. Schmidt, [Numerical Simulation of Optical Wave Propagation], SPIE Press, Bellingham, Wash. (2010).
• [14] S. R. Bose-Pillai, J. E. McCrae, C. A. Rice, R. A. Wood, C. E. Murphy, and S. T. Fiorino, “Estimation of atmospheric turbulence using differential motion of extended features in time-lapse imagery,” Opt. Eng. 57(10), 104108 (2018).
• [15] http://vortex.plymouth.edu/myo/sfc/statlog-a.html.

While the disclosure has been described with reference to exemplary embodiments, it will be understood by those skilled in the art that various changes may be made and equivalents may be substituted for elements thereof without departing from the scope of the disclosure. In addition, many modifications may be made to adapt a particular system, device or component thereof to the teachings of the disclosure without departing from the essential scope thereof. Therefore, it is intended that the disclosure not be limited to the particular embodiments disclosed for carrying out this disclosure, but that the disclosure will include all embodiments falling within the scope of the appended claims. Moreover, the use of the terms first, second, etc. do not denote any order or importance, but rather the terms first, second, etc. are used to distinguish one element from another.

In the preceding detailed description of exemplary embodiments of the disclosure, specific exemplary embodiments in which the disclosure may be practiced are described in sufficient detail to enable those skilled in the art to practice the disclosed embodiments. For example, specific details such as specific method orders, structures, elements, and connections have been presented herein. However, it is to be understood that the specific details presented need not be utilized to practice embodiments of the present disclosure. It is also to be understood that other embodiments may be utilized and that logical, architectural, programmatic, mechanical, electrical and other changes may be made without departing from general scope of the disclosure. The following detailed description is, therefore, not to be taken in a limiting sense, and the scope of the present disclosure is defined by the appended claims and equivalents thereof.

References within the specification to “one embodiment,” “an embodiment,” “embodiments”, or “one or more embodiments” are intended to indicate that a particular feature, structure, or characteristic described in connection with the embodiment is included in at least one embodiment of the present disclosure. The appearance of such phrases in various places within the specification are not necessarily all referring to the same embodiment, nor are separate or alternative embodiments mutually exclusive of other embodiments. Further, various features are described which may be exhibited by some embodiments and not by others. Similarly, various requirements are described which may be requirements for some embodiments but not other embodiments.

It is understood that the use of specific component, device and/or parameter names and/or corresponding acronyms thereof, such as those of the executing utility, logic, and/or firmware described herein, are for example only and not meant to imply any limitations on the described embodiments. The embodiments may thus be described with different nomenclature and/or terminology utilized to describe the components, devices, parameters, methods and/or functions herein, without limitation. References to any specific protocol or proprietary name in describing one or more elements, features or concepts of the embodiments are provided solely as examples of one implementation, and such references do not limit the extension of the claimed embodiments to embodiments in which different element, feature, protocol, or concept names are utilized. Thus, each term utilized herein is to be given its broadest interpretation given the context in which that terms is utilized.

The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the disclosure. As used herein, the singular forms “a”, “an” and “the” are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms “comprises” and/or “comprising,” when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof.

The description of the present disclosure has been presented for purposes of illustration and description, but is not intended to be exhaustive or limited to the disclosure in the form disclosed. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope of the disclosure. The described embodiments were chosen and described in order to best explain the principles of the disclosure and the practical application, and to enable others of ordinary skill in the art to understand the disclosure for various embodiments with various modifications as are suited to the particular use contemplated.

Claims

1. A turbulence measurement system comprising: a processor that is communicatively connected to an image capturing device that captures time-lapse images of a distant target that is anywhere from a km to more than a 100 km away, and which: receives the time-lapse images;tracks relative motion due to atmospheric turbulence of some number of patches of definite size on each of the time-lapse images using a sub-pixel accurate correlation technique;computes differential motion between all pairs of patches on each image, this differential motion is simply proportional to the differential wavefront tilt;computes differential tilt variance for each pair of patches from the time-lapse images;evaluates theoretical weighting functions between turbulence along the imaging path and the measured differential tilt variances, wherein the weighting functions depend upon a size and a separation of the patches being tracked as well as a path length and a size of a camera aperture of the image capturing device;determines a linear combination of the theoretical weighting functions that best match a weighting function for a selected turbulence parameter of interest; andlinearly combines the measured tilt variances to estimate the selected turbulence parameter of interest.

2. The turbulence measurement system of claim 1, wherein the selected turbulence parameter of interest comprises one of a group consisting of an isoplanatic patch size, log-amplitude variance, Fried's coherence diameter, and turbulence strength in a neighborhood around a specified range.

3. The turbulence measurement system of claim 1, wherein the image capturing device is a digital camera operating in the visible or near infrared.

4. The turbulence measurement system of claim 1, wherein the image capturing device captures images no faster than an atmospheric coherence time, with individual exposures shorter than atmospheric coherence time.

5. The turbulence measurement system of claim 1, wherein the processor, using the selected turbulence parameter of interest, compensates for degraded images captured by a surveillance system.

6. The turbulence measurement system of claim 1, wherein the processor, using the measured differential tilt variances and the corresponding computed weighting functions, determines how turbulence is distributed along the optical path in a tactical engagement scenario.

7. An optical system comprising: an image capturing device that captures time-lapse images of a distant target (anywhere from a km to more than a 100 km away) along an imaging path; anda turbulence measurement system, comprising a processor that is communicatively connected to the image capturing device and which: receives the time-lapse images;tracks relative motion due to atmospheric turbulence of some number of patches of definite size on each of the time-lapse images using a sub-pixel accurate correlation technique;computes differential motion between all pairs of patches on each image, this differential motion is simply proportional to the differential wavefront tilt;computes measured differential tilt variance for each pair of patches from the time-lapse images;evaluates theoretical weighting functions between turbulence along the imaging path and the measured differential tilt variances, wherein the weighting functions depend upon a size and a separation of the patches being tracked as well as a path length and a size of a camera aperture of the image capturing device;determines a linear combination of the theoretical weighting functions that best match a weighting function for a selected turbulence parameter of interest; andlinearly combines the measured tilt variances to estimate the selected turbulence parameter of interest.

8. The optical system of claim 7, wherein the selected turbulence parameter of interest comprises one of a group consisting of an isoplanatic patch size, log-amplitude variance, Fried's coherence diameter, and turbulence strength in a neighborhood around a specified range.

9. The optical system of claim 7, wherein the image capturing device is a digital camera operating in the visible or near infrared.

10. The optical system of claim 7, wherein the image capturing device captures images no faster than an atmospheric coherence time, with individual exposures shorter than atmospheric coherence time.

11. The optical system of claim 7, wherein the processor, using the selected turbulence parameter of interest, compensates for degraded images captured by a surveillance system.

12. The optical system of claim 7, wherein the processor, using the measured differential tilt variances and the corresponding computed weighting functions, determines how turbulence is distributed along the optical path in a tactical engagement scenario.

13. A method comprising: receiving time-lapse images from an image capturing device that captures time-lapse images of a distant target (anywhere from a km to more than a 100 km away) along an imaging path;tracking relative motion due to atmospheric turbulence of some number of patches of definite size on each of the time-lapse images using a sub-pixel accurate correlation technique;computing differential motion between all pairs of patches on each image, this differential motion is simply proportional to the differential wavefront tilt;computing measured differential tilt variance for each pair of patches from the time-lapse images;evaluating theoretical weighting functions between turbulence along the imaging path and the measured differential tilt variances, wherein the weighting functions depend upon a size and a separation of the patches being tracked as well as a path length and a size of a camera aperture of the image capturing device;determining a linear combination of the theoretical weighting functions that best match a weighting function for a selected turbulence parameter of interest; andlinearly combining the measured tilt variances to estimate the selected turbulence parameter of interest.

14. The method of claim 13, wherein the selected turbulence parameter of interest comprises one of a group consisting of an isoplanatic patch size, log-amplitude variance, Fried's coherence diameter, and turbulence strength in a neighborhood around a specified range.

15. The method of claim 13, wherein the image capturing device is a digital camera operating in the visible or near infrared.

16. The method of claim 13, wherein the image capturing device captures images no faster than an atmospheric coherence time, with individual exposures shorter than atmospheric coherence time.

17. The method of claim 13, further comprising using the selected turbulence parameter of interest to compensate for degraded images captured by a surveillance system.

18. The method of claim 13, further comprising using the measured differential tilt variances and the corresponding computed weighting functions to determine distribution of turbulence along the optical path in a tactical engagement scenario.