Turbulence profiling using pupil plane wavefront data derived Fried parameter values for a dynamically ranged Rayleigh beacon

Steven M. Zuraski; Elizabeth Beecher; Jack E. McCrae; Steven T. Fiorino

doi:10.1117/1.OE.59.8.081807

14 May 2020 Turbulence profiling using pupil plane wavefront data derived Fried parameter values for a dynamically ranged Rayleigh beacon

Steven M. Zuraski, Elizabeth Beecher, Jack E. McCrae, Steven T. Fiorino

Author Affiliations +

Optical Engineering, Vol. 59, Issue 8, 081807 (May 2020). https://doi.org/10.1117/1.OE.59.8.081807

Abstract

Long-range optical imaging applications are typically hindered by atmospheric turbulence. The effect of turbulence on an imaging system can manifest itself as an image blur effect usually quantified by the phase distortions present in the system. The blurring effect can be understood on the basis of the measured strength of atmospheric optical turbulence along the propagation path and its impacts on phase perturbation statistics within the imaging system. One method for obtaining these measurements is by the use of a dynamically ranged Rayleigh beacon system that exploits strategically varied beacon ranges along the propagation path, effectively obtaining estimates of the aberrations affecting an optical imaging system. We developed a method for extracting tomographic turbulence strength estimations from a dynamically ranged Rayleigh beacon system that uses a Shack–Hartmann sensor as the phase measurement device. The foundation for extracting tomographic information from strategically range-varied beacon measurements obtained in rapid sequence is presented along with modeled example cases for typical turbulence scenarios. Additionally, the processing algorithm was used to simulate identification of isolated strong turbulence layers. We present the chosen processing algorithm’s foundation and provide discussion of the utility of this algorithm as an atmospheric turbulence profiling methodology.

1. Introduction

Long-range optical imaging applications are typically hindered by atmospheric turbulence. The pseudorandom variations in the index of refraction are the physical observables resulting from atmospheric turbulence. This effect leads to unwanted optical blurring of the image. As a result, there is a need to understand how these pseudorandom variations in the index of refraction behave and evolve along the viewing path. Furthermore, knowledge of the profile strength of the index of refraction structure parameter can lead to mitigation strategies for combating the undesired effects. A previously presented conceptual strategy¹ has been developed for measuring the strength of the refractive index structure parameter profile utilizing a dynamically ranged Rayleigh beacon system. A dynamically ranged Rayleigh beacon system is a modification of a traditional Rayleigh beacon system design to allow for the originating location of the backscattered field to change in the range from the collecting aperture of the telescope. When this dynamically ranged methodology is employed, a new set of data processing algorithms can be used to extract information about localized turbulence strength contributions as opposed to just a single volumetric estimation of turbulence strength, which is the output from traditional laser beacon systems. This paper presents an approach for extracting localized turbulence strength information based on the methodology employed by a dynamically ranged Rayleigh beacon; the associated data exploitation algorithm that has been developed and investigated can accurately produce tomographic refractive index structure parameter strength profiles. Two example cases are used in this paper to highlight this capability.

Other similar methods for sensing the tomographic strength of turbulence over a path or overcoming turbulence effects in a discretized step manner using beacon-based measurements have previously been proposed, and those have provided an inspirational basis for the methods presented in this paper. One such method involves three-dimensional (3-D) tomographic projection where wavefront sensor measurements are used to directly measure the phase in a two-dimensional (2-D) plane orthogonal to the direction of propagation. Then measurements taken at various angles through the turbulent flow were used to reconstruct the 3-D structure of the index of refraction.² A second method is a byproduct of a proposed multiconjugate adaptive optics (MCAO) system in which multiple Rayleigh beacons are used to measure needed phase corrections over a wide field of regard. By intentionally choosing beacons at varied ranges, turbulent layer range misregistration in the MCAO system can be sensed and avoided. It was also proposed that intentionally varied ranges from multiple beacons could be used to find strong contributing turbulent layers, and thus reduce the effect of focal anisoplanatism and improve stability of the point spread function across larger field angles.³ Additionally, wavefront sensor data from Rayleigh beacons have been used to produce propagated laser beams with improved performance using phase estimates to correct scintillation effects. This type of measurement and response has been proposed to increase the effective range of a laser weapon system through sequential compensation iterations within a slowly changing or nearly static atmosphere. This proposes a laser beam projection system with a near-field phase retrieval method for delivering more power to a target at an extended range.⁴ These works have provided inspirational aspects influencing the design concept of a dynamically ranged Rayleigh beacon system¹; however, each individual method presents differences that influence that way the collected data is analyzed to produce estimates of the strength of the refractive index structure parameter, associated turbulence-related metrics, or near-real-time phase correction. The concept of a dynamically ranged Rayleigh beacon system is unique compared to these because of its ability to control the backscattered field originating range on a pulse-by-pulse basis employed by a single on-axis beacon. This system design introduces a new set of possibilities for extracting localized turbulence strength information and consequently presents an opportunity for innovative algorithm development. Specific to this paper, the data exploitation algorithm that was developed is a new means for extraction of tomographic estimations of the refractive index structure parameter specifically tied to the dynamically ranged Rayleigh beacon-based data collection system.

There are multiple technologies for producing tomographic estimations of the strength of turbulence. These technologies did not provide inspiration for the system or processing techniques in this paper, but are the alternative technologies that aim at the same goal of estimating the profiled strength of turbulence along a viewing path. These alternative technologies utilize techniques involving slope detection and ranging (SLODAR), sound detection and ranging (SODAR), scintillation detection and ranging (SCIDAR), differential image motion LiDAR (DIM LiDAR), and multiaperture scintillation system (MASS). Each technique has distinctive advantages and disadvantages as compared to the dynamically ranged Rayleigh beacon technique. SLODAR is an optical triangulation technique utilizing two or more beacon sources. Systems using laser beacon sources are still quite new, but more established techniques utilizing paired natural stars have been implemented at select observatories for a number of years. The geometric regional overlap in the propagated waves from the two separated sources allows for estimation of $C_{n}^{2}$ relatively near the aperture in resolution bins defined by the spatial location of subapertures in the pupil plane relative to the geometric separation of the beacon sources. SLODAR processing is dependent on the statistical correlation between wavefront sensor slopes.⁵^–⁷ A dynamically ranged beacon system is similar in sensor design to a SLODAR system; however, the way in which the data are processed to produce a turbulence profile estimate is vastly different. A SLODAR system analyzes slope correlations from angularly separated beacons, whereas a dynamically ranged beacon system utilizes spatial variances and differentiation between sequential beacons from separated ranges along the optical axis of the sensing system. Next, a SODAR system utilizes sound waves in a similar way to radar and measures the structure constant of temperature $C_{T}^{2}$ . The structure constant of temperature is then related to the structure constant of the index of refraction.⁸ This is an inferred method for producing a profile estimation of atmospheric turbulence and does not directly measure any perturbed wavefront. SCIDAR is another optical triangulation technique similar to SLODAR, except the measurement taken is a scintillated intensity pattern, which is a second-order effect as it is dependent on the second derivative of the optical wavefront. From the intensity pattern, a scintillation index can be calculated through correlation of the multiple intensity patterns, and a resultant profile of turbulence strength can be estimated.⁸^–¹⁰ An adaptation of a SCIDAR system is the MASS. The MASS was designed to overcome the deficiencies of a SCIDAR system, which include the requirement for a bright double star and a large aperture telescope. The MASS utilizes a single star and multiple apertures as spatial filters to sense the photon counts associated with scintillation. The mapping of the photon statistics in each aperture can be traced to turbulence effects originating from specific ranges.¹¹^,¹² Finally, a DIM LiDAR is the incorporation of a DIM monitor and a LiDAR system. The DIM part utilizes a technique that measures the variance of the differential wavefront tilt by two or more spatially separated apertures. The variance statistics are gathered over a relatively long timeframe to produce an averaged value of the turbulence strength between the collecting apertures and guide star. The LiDAR part comes from changing the range of the beacon. This allows for the discretization of the averaged turbulence strength estimates to build up a profile estimate.¹³

These alternative methods described have underlying processing techniques and inherent phenomenology characteristics that are distinctively different from those described by the methods proposed for a dynamically ranged beacon. Multiple methods rely on statistical spatial correlations to build up an estimate of turbulence strength. A dynamically ranged beacon system could exploit similar correlations under subaperture crossed-path geometries between two different range beacons; however, the methodology is not reliant on that to produce a turbulence strength profile estimate. Instead, a metric of wavefront variance across an aperture is used for volumetric estimates of turbulence strength and then a differentiation algorithm is employed based on the dynamically ranged beacon location to produce localized turbulence strength estimates. The dynamically ranged beacon system utilizes a series of direct measurements of the wavefront present in the system’s pupil as the input to producing a profile estimate of the turbulence strength. This is inherently different from methods where associated measurements are made, such as $C_{T}^{2}$ , and inferred relations are used to get back to an index of refraction structure parameter estimate. It is also inherently different from the case of building up variance statistics from a few apertures measuring modal tilts and applying DIM algorithms to assess turbulence strength. The fact that a wavefront is captured at a finer resolution so that the zonal tilts are the dominating aberrations in each section means a wavefront spatial variance statistic can be used to produce a near-instantaneous estimate of the integrated turbulence strength. This coupled with a rapidly changed beacon location can build up a turbulence strength profile.

The exploitation of measurements produced by a dynamically ranged Rayleigh beacon is based on the formulation of the Fried parameter. The estimation of the Fried parameter from each wavefront produced provides a metric for assessing the strength of the turbulence along the integrated viewing path. On an individual laser pulse basis and corresponding single measurement, this is nominal treatment for data of this type.¹⁴ However, tuning the range gate timing for a turbulence profiling purpose or investigation of a concentrated layer’s structure can be done to achieve further characterization than that of an integrated volume measurement. Moreover, under the assumption that the atmosphere is frozen between a measurement and the immediate next measurements in sequence within the specified fraction of a second time period, a discretized range resolved estimation of the strength of the optical turbulence can be produced. The data processing algorithm for constructing a turbulence strength profile, the basis for the algorithm’s formulation, and the results from modeled scenarios are presented in this paper.

2. Key Concepts to Support Algorithm Development and Evaluation

2.1.

Optical Turbulence Metrics

Introduced by Fried (1966), the Fried parameter $r_{0}$ defines the diameter of a circular pupil that would produce an equivalent diffraction-limited full-width at half-maximum of a point source image as the atmospheric turbulence would with an infinite in an extent mirror. This metric is correlated with the strength of the atmospheric refractive index fluctuations. The root-mean-square (rms) phase variation over a circular aperture of diameter $r_{0}$ is $⟨ σ^{2} ⟩ = 1.03 {rad}^{2}$ . Spatial wavefront patches of size $r_{0}$ can be regarded as planar phase regions within the circular pupil.³ The pupil used for measurements should be divided into many subapertures that are smaller than $r_{0}$ in extent resulting in many individual planar regions that can be accurately measured by a conventional Shack–Hartmann wavefront sensor (SHWFS). The individual spot measurements of a SHWFS can then be used to reconstruct the nonplanar phase of the pupil commonly known as the wavefront. This nonplanar phase is used for estimation of the atmospheric $r_{0}$ for the sensed volume under the Rayleigh beacon. For a spherical wave, $r_{0}$ can be expressed in terms of a weighted integral:

Eq. (1)

r_{0} = {[0.423 k^{2} \int_{0}^{L} C_{n}^{2} (z) {(\frac{z}{L})}^{5 / 3} d z]}^{- 3 / 5},

where

k = 2 π / λ

is the wavenumber,

λ

is the wavelength,

L

is the distance to the source, and

C_{n}^{2}

is the refractive index structure parameter.

C_{n}^{2}

is a metric used to quantify the strength of turbulence at specified distances

z

along the viewing path. The units are in

m^{- 2 / 3}

and are typically in the range of

1 \times 10^{- 13}

to

1 \times 10^{- 17}

in the lower atmosphere, but turbulence can be stronger or weaker depending on the location and viewing scenario. For this treatment, we will assume that

the C_{n}^{2}

power spectrum follows the Kolmogorov model form.¹⁵

The Earth’s atmosphere refractive index structure evolves over time and space in a random fashion under the Kolmogorov statistical model. This causes light to be distorted as it propagates through the atmosphere. In most theories of atmospheric turbulence, turbulent flow kinetic energy is transferred from large eddies to small eddies in a cascading fashion until the energy is dissipated. The average size of the large eddies is the outer scale $L_{0}$ , and the average size of the smallest eddies is the inner scale $l_{0}$ . The range of eddy sizes between the inner and outer scales is called the inertial subrange. The Kolmogorov model only pertains to the inertial subrange and ignores inner and outer scale effects. Electromagnetic propagation takes place within Earth’s atmosphere at the speed of light. Therefore, within relatively short periods of time, the refractive index changes present within regions of the atmosphere can be considered stationary in time and fixed in position. This is because the speed of light is comparatively fast, and light can travel much farther than even the span of the largest of turbulent eddy cells. Consequently, the temporal properties of atmospheric turbulence can be regarded as static for these same short periods of time.¹⁶

In atmospheric propagation scenarios that are modeled using wave optics techniques, turbulence can be treated as a finite number of discrete layers. Each layer is represented by a phase screen that is a flattened projection representation of phase variations in a much larger volume. Each phase screen is a singular realization of the atmospheric turbulence experienced by the propagating wave for a designated volume. This effect can be summarized by dividing the total volume into discrete segments. According to Andrews and Phillips,¹⁷ the discrete sum version of the Fried parameter for a spherical wave is

Eq. (2)

r_{0, SW} = {[0.423 k^{2} \sum_{i = 1}^{n} C_{n_{i}}^{2} {(\frac{z_{i}}{Δ z})}^{5 / 3} Δ z_{i}]}^{- 3 / 5},

where

Δ z_{i}

is the thickness of the turbulence volume and

z_{i}

is the location along the propagation path. A segmented form of an effective coherence diameter is

Eq. (3)

r_{0_{i}} = [0.423 k^{2} C_{n_{i}}^{2} Δ z_{i}]^{- 3 / 5} .

This expression for the effective coherence diameter is the plane wave-based solution, which differs from the spherical wave solution. This removes the effect of weighting the metric toward the receiver.¹⁷ Since the employment of the described algorithms is based on differentiation of subsequent measurements, the effects near the receiver are canceled out and the simpler form of a plane wave can be used. However, the physical nature of the wave propagation from a beacon is that of a spherical wave, so it is important to utilize spherical wave-based equations for total integrated volumes. This can consequently be substituted into the discrete sum version of the Fried parameter for a spherical wave to yield

Eq. (4)

r_{0, SW} = [\sum_{i = 1}^{n} r_{0_{i}}^{- 5 / 3} {(\frac{z_{i}}{Δ z})}^{5 / 3}]^{- 3 / 5} .

Since it is shown that the total path Fried parameter can be thought of as a sum of discrete segment measurements or estimates, it is consequently justifiable to utilize subsequent measurements of the Fried parameter where the atmosphere is considered stationary to derive the Fried parameter segment strength. This is represented mathematically as $r_{0_{i}} = r_{0_{j + 1}} - r_{0_{j}}$ , where $j + 1$ and $j$ represent two measurements of $r_{0}$ from different beacon ranges. Under these conditions, the Fried parameter segment strength can be used to calculate the refractive index structure parameter in discretized layers. This metric is calculated as

Eq. (5)

C_{n_{seg}}^{2} = \frac{r_{0_{i}}^{- 5 / 3}}{0.423 k^{2} Δ z_{i}},

where

r_{0_{i}}

is the Fried parameter segment for two subsequent path measurements. The refractive index structure parameter segments can then be recombined to build up the total profile based on the choices of ranges used as part of the concept of operations (CONOPS) employed by the dynamically ranged Rayleigh beacon system. This treatment lays the foundation for the simulation and evaluation of the profiled nature of atmospheric turbulence strength as measured by a dynamically ranged Rayleigh beacon system that operates utilizing

n

phase measurements configured such that structure effects are measurable.

2.2.

Simulated Wave Propagation

Simulations were carried out in a multistep propagation method described by Schmidt¹⁸ utilizing the metrics presented previously. For a dynamically ranged Rayleigh beacon system scenario, a point source is projected at distance $L$ , and a multistep propagation is performed where a sequence of $n$ phase screens is utilized. The $n$ phase screens used for simulation are chosen in number to at minimum match the resolution desired by the methodology described in Eq. (5). The number of phase screens chosen in a simulation to match practice would be limited by the repetition rate of the laser system used and the assumed time during which the atmosphere can be considered statistically frozen. Standard Fourier optic wave propagation methods were used to simulate the beacon projection and returned light field. This method accurately captured the modeled effects of the layered atmosphere and how a Shack–Hartmann sensor system captured the phase information present in the telescope’s pupil plane from the returned optical energy originating from varied ranges. The data from multiple measurements were aggregated within an assumed time window constraint and were used to feed to the methodology associated with Eq. (5) to form a profile estimate. Many Shack–Hartmann-based measurements were collected to build up an ensemble of realizations for randomly different atmospheres all having the same turbulence strength injections. This collection of data was used to build up statistically significant metrics for utilizing atmosphere propagated wavefronts originating at strategically varied Rayleigh beacon ranges. This was done to ensure a specific SHWFS-based phase measurement of a single turbulence scenario was not an anomaly and was representative of how the methodology presented could perform.¹⁹ Within the framework of a Fourier optic simulation of the wave propagation of a dynamically ranged Rayleigh beacon, accumulated measurements from individualized propagations can be used to provide measurement-based tomographic turbulence strength estimations.

2.3.

Shack–Hartmann Wavefront Sensor Measurements

A SHWFS is utilized to measure the incident phase on the collecting aperture of the sensing system. The SHWFS consists of a lenslet array mapped to the pupil plane of the collecting aperture and a camera that is placed at the focal plane of the lenslets. The lenslet focal length and spacing determine the systems sensitivity and dynamic range for sensing wavefront tilt. The number of lenslets in the array determines the resolution of the wavefront measuring system as projected to the pupil plane of the collecting aperture. Under the simulation scenario utilizing Kolmogorov turbulence statistics, the incoming wave has zonal phase tilt aberrations that will produce a nonuniform grid of spots imaged by the SHWFS system. The spot locations, one mapped to each lenslet, are used to evaluate the local gradient of the aberrated wave and can be reconstructed to produce a 2-D segmented version of the incoming wave. Measurement of local gradients and production of a resultant measured wavefront are formed utilizing tilt sensing methods and gridded optical path difference reconstruction matrices. The reconstruction matrices relate the measured local gradients to the wavefront nodes to produce an estimated wavefront fit to the uniform grid.²⁰

Within the simulation, lenslet mapping sizes are chosen such that they are smaller than the smallest Fried parameter experienced. This is a key design parameter because if the lenslet mapping to the telescope’s pupil plane was significantly larger than $r_{0}$ , then the atmospheric induced tilt would no longer be the dominate aberration in each lenslet and the measurements would not be representative of the true distorted wavefront present in the telescope’s pupil. The Fried parameter and associated requirement for lenslet array size are precalculated after the chosen profile is generated and the atmospheric path propagation has taken place for the longest beacon range. In practice, the range of typical turbulence strengths for a measurement site is assumed and from this the lenslet mapping size can be determined. The lenslet mapping is typically a fixed value based on SNR constraints combined with the desired spatial resolution required for accurately sensing the wavefront present in the pupil plane. For simulation purposes, a $10 \times 10$ or larger grid of lenslets was used. This ensured enough measurements were made to reconstruct a meaningful wavefront and also matched well with the LiDAR equation-based SNR model associated with the proposed dynamically ranged Rayleigh beacon system.¹⁹^,²⁰

2.4.

Profiled Turbulence Metrics

In simulation, the turbulence strength modeled along the optical path and the phase screens used at each multistep junction have known optical properties. These optical properties stem from an $r_{0}$ data metric, which is calculated from an input $C_{n}^{2}$ strength profile applied to a controlled random phase screen generation.¹⁵ These input metrics are used as truth data for comparison against the modeled data and outlined algorithm output. How well the algorithm’s segmented $C_{n}^{2}$ outputs matched the input $C_{n}^{2}$ strength profile served as a merit function for evaluation. Two $C_{n}^{2}$ profiles used for simulation presented in this paper are shown in Fig. 1.

Fig. 1

The two turbulence strength profiles used for simulation in this paper.

From the simulated data, a method for exploitation of a dynamically ranged Rayleigh beacon system was evaluated. We utilized the SHWFS zonal measurements to produce an $r_{0}$ metric based on the flatness of the wave, consistent with Fried’s requirements for the measurement of $r_{0}$ . This method is further explained by Eqs. (6) and (7) and the associated descriptions. We then utilized those metrics for estimation of $C_{n}^{2}$ based on subsequent layer measurements and Eq. (5). This segmented estimate of $C_{n}^{2}$ was used as the final output and compared against the simulation inputs. Due to the random nature of simulating phase screen representations of pupil distortions that result from a modeled turbulent volume, one would expect small deviations between the estimated $C_{n}^{2}$ profile and the input data.

3. Algorithm Implementation and Evaluation

The data collection scenario consists of a dynamically ranged Rayleigh beacon where pupil plane wavefront estimation data are taken on a laser pulse by laser pulse basis in rapid succession such that the turbulent volume can reasonably be assumed to be stationary. Each individual pupil plane wavefront estimation data collect is done by a SHWFS model and corresponds to a turbulent volume integrated along the viewing path. A projected phase screen is reconstructed from the SHWFS data. Using the mean of the variance of the phase across the aperture $r_{0}$ can be estimated for the designated volume corresponding to the conic volume between the Rayleigh beacon and collecting aperture. Utilizing the collection of rapid succession measurements and Eq. (5), $C_{n}^{2}$ segments can be estimated and combined to form an estimate of the $C_{n}^{2}$ strength profile.

3.1.

Overview of Simulation Implemented

The simulation consisted of an ensemble of 200 realizations for each input profile. Two of the profiles used are presented and are shown in Fig. 1. The first is a standard Hufnagel-Valley 5/7 profile and the second is a Hufnagel-Valley 5/7 profile with enhanced turbulence strength localized in two altitude regions. The second profile is used to probe the algorithm’s effectiveness for locating strong changes in turbulence strength in neighboring stratified layers.

These refractive index structure parameters were used as the seed strengths for generation of pseudorandom Kolmogorov phase screens. The same random number generator seed was used for a single realization, which consisted of multiple propagations from the beacon at the varied ranges to the collecting aperture. This assumption built into the simulation is intended for the use of probing the algorithm’s effectiveness. In practice, the number of sequential measurements taken from varied path length ranges will depend on how long the atmosphere can be considered frozen, the speed of the optical shutter, the repetition rate of the laser system used, and the SNR of the returned laser light from a beacon making a phase estimation measurement possible.

The phase present at the collecting aperture was analyzed using a SHWFS model. For the simulations presented in this paper, a $10 \times 10$ lenslet grid was masked to a 60.96-cm circular aperture. This resulted in each subaperture’s size extent being 6.096 cm, which is smaller than any of the $r_{0}$ input values for a single beacon measurement. This ensures that wavefront tilt on a subaperture is the dominant aberration present. In the simulation, each lenslet region was $\sim 100 pixels$ across and each resultant spot produced from a lenslet covered a $3 \times 3 pixel$ box. The intensity distribution with the $3 \times 3$ box was used to estimate subpixel centroid location. The subpixel centroid, the location within the lenslet region of interest, and the lenslet focal length were used together to deduce the slope of the optical field present within a lenslet region. An example of the resultant slope vectors is shown in Fig. 3. For this analysis, focal shifts within the sensing system due to changing the range between the beacon and the aperture were neglected. In practice, changing the range of the beacon will shift the location of the focus in a telescope system. The design of the optical relay and location of optical components such as the shutter used need to be chosen with careful consideration specific to the telescope used. It is possible under constraints to keep the Shack–Hartmann sensor in collimated space with mapped lenslet sizes remaining smaller than $r_{0}$ . Therefore, for the simulated treatment, this effect was ignored.

3.2.

Zonal Measurements to Estimate $r_{0}$

Each propagation produced one wavefront present at the collecting aperture. This was measured by the SHWFS. An example of the SHWFS system model used is shown in Fig. 2 along with an example of the input wavefront seen by the sensing system, which is the turbulence-induced phase distortions present in the system’s pupil used for analysis. Figure 3 shows an example of the zonal gradients produced from the sensed wavefront and the resultant segmented reconstructed wavefront. The wave is sensed using a SHWFS, which converts focused spots to local tilt measurements using a spot centroiding algorithm on the imaged focal plane. The resultant slopes from each centroid are representative of zonal tilts in the pupil. These tilts are used to reconstruct the wave using a matrix inversion reconstruction method.²⁰ From each reconstructed wave, an estimate of the Fried parameter was produced based on the wavefront variance statistics. In simulation, this was done by tiling the simulated phase present at the collection aperture with circular pupils representing the lenslets of the SHWFS. Each tile produced a tilt measurement that was converted into a zonal curvature through reconstruction. The mean of the variance of the phase across the aperture was then calculated to derive the Fried parameter. This is shown in Eq. (6). The mean of the variance of the phase across the aperture is

Eq. (6)

⟨ σ^{2} ⟩ = 0.134 {(\frac{D}{r_{0}})}^{5 / 3},

where

D

is the aperture diameter.¹⁴ Manipulation of Eq. (6) to solve for

r_{0}

results in

Eq. (7)

r_{0} = \frac{0.299 D}{{⟨ σ^{2} ⟩}^{3 / 5}},

where this

r_{0}

is a representative metric of the turbulence strength within the volume of a single beacon measurement,

D

is the aperture diameter, and

⟨ σ^{2} ⟩

is the mean of the variance of the estimated phase across the aperture that is built from the zonal tilt tiles reconstructed from the SHWFS measured gradients. These single beacon measurements are accumulated and referenced as

1, 2, \dots, j, j + 1

, where each measurement comes from a different beacon range. These are then used as inputs to Eq. (5) to produce a segmented profile estimate of the atmospheric turbulence strength.

Fig. 2

(a) SHWFS model used for simulations and (b) example of turbulence-induced phase distortions present in system’s pupil.

Fig. 3

Example of a wavefront (a) sensed and (b) reconstructed.

This methodology of using a SHWFS to estimate $r_{0}$ from a pupil image whose phase is corrupted by turbulence is susceptible to small errors that can propagate forward into the dynamically ranged Rayleigh beacon profiling algorithm. These errors could stem from the gridded lenslet array that is mapped to the pupil. Each lenslet provides a point estimate of the average tilt contour within the mapped region, which is the dominant aberration. However, high-order aberrations do exist that are not captured by this methodology. In addition, the mapping of the lenslets to the pupil inherently has discrete regions with sharp edges that may fall into nonideal locations depending on the turbulence-induced phase distortions present in the pupil. This has the potential to average a large tilt between two subapertures, thus reducing contributing terms to the variance. The opposite is also true, and a large tilt may be entirely captured within one subaperture resulting in a large contributing term to the variance calculation. Since the turbulence-induced phase distortions present in the pupil are random in nature, these effects should average out through many measurements. Additionally, the wavefront reconstruction process naturally reduces the influence of a single large spike in a local tilt measurement. However, the dynamically ranged Rayleigh beacon profiling algorithm relies on a singular series of adjacent measurements to create a profile, and therefore, is susceptible to small errors in singular turbulence strength profile estimations.

The phase incident on the collecting aperture associated with each realization is discretely sampled at a high rate as compared to the Fried parameter or the number of lenslets in the sensing system. This is evident in comparing the example phase shown in Fig. 2 to the reconstructed phase in Fig. 3. The discretization level is chosen such that the width of a subaperture is smaller than the Fried parameter. This chosen region width mapped to the collecting aperture is equal to the physical diameter of a lenslet projected onto the primary mirror of the collecting telescope system. Due to the discretized down sampling nature of this physical process, small deviations from the truth metric could be induced. However, due to the constraints imposed, these small deviations should be minimized and will have minimal overall effect on the system.

From the measured zonal gradients and resultant reconstructed wavefronts, the Fried parameter was estimated. Fried parameter estimates compared to true Fried parameters for the individualized beacon propagations are shown in Fig. 4. Each data point in Fig. 4 consisted of a propagation from a single altitude beacon to the sensing system. The sensing system then estimates an $r_{0}$ value from a reconstructed wavefront like that shown in Fig. 3 to produce a single $r_{0}$ value. The $r_{0}$ value is estimated by manipulating Eq. (6) to solve for $r_{0}$ as shown in Eq. (7). Each black data point is from a series of $r_{0}$ estimates over the altitude range where the random number generation seed was controlled to be the same. The whole ensemble consists of 200 independent realizations each drawn as black data points. The true $r_{0}$ is plotted as a solid line and was used as the input metric for the simulation. The true $r_{0}$ is the approximate average of the estimated $r_{0}$ values comprising the whole ensemble.

Fig. 4

Ensemble of $r_{0}$ versus altitude estimates (data points) compared to truth data (solid line). Note: $r_{0}$ decreases with altitude because the path length increases with each altitude step. (a) The standard Hufnagel-Valley 5/7 profile and (b) the Hufnagel-Valley 5/7 profile with enhanced turbulence strength localized in two altitude regions.

3.3.

Influence of Focal Anisoplanatism

When choosing the locations of the dynamically ranged beacons, careful consideration needs to be taken to avoid negative effects associated with focal anisoplanatism. Focal anisoplanatism arises from the beacon geometry. The effects of focal anisoplanatism are typically minimized by placing the beacon as far away as possible.²⁰ However, for a dynamically ranged Rayleigh beacon system, the location of the beacon is intentionally varied and not at a maximized distance from the telescope. Additionally, the treatment of focal anisoplanatism is applied in a differential manner as opposed to the traditional treatment against a plane wave. This has to be done because of the unique dynamically ranged operation of the beacon system. Hardy²⁰ admittedly states that calculation of the focal anisoplanatism error has proved to be a difficult task, but provides a method developed by Belsher and Fried in 1994,²¹ which has matched numerical evaluations within 1%. This metric for focal anisoplanatism error from a single beacon measurement is

Eq. (8)

σ_{FA}^{2} = {(\frac{D}{d_{0}})}^{5 / 3},

where

D

is the aperture diameter and

d_{0}

is the diameter over the aperture in which the wavefront error due to focal anisoplanatism is

< 1 {rad}^{2}

. For comparing subsequent dynamically ranged beacon measurements against each other while considering the effects of focal anisoplanatism,

d_{0}

needs to be derived. A single beacon

d_{0}

can be described as

Eq. (9)

d_{0} = {k^{2} {[0.057 μ_{0}^{+} (H) + 0.500 \frac{μ_{5 / 3}^{-} (H)}{H^{5 / 3}} - 0.452 \frac{μ_{2}^{-} (H)}{H^{2}}]}}^{- 3 / 5},

where

H

is the range to the beacon from the telescope and

μ_{m}

is the turbulence moment described in terms of the upper and lower moments for

m = 0

, 5/3, and 2. The upper turbulence moment is

Eq. (10)

μ_{m}^{+} (H) = \int_{H}^{\infty} C_{n}^{2} (z) z^{m} d z,

and the lower turbulence moment is

Eq. (11)

μ_{m}^{-} (H) = \int_{0}^{H} C_{n}^{2} (z) z^{m} d z .

Relating this to the dynamically ranged beacon scenario, a differential metric can be developed which results in

Eq. (12)

Δ d_{0} = (k^{2} {0.057 [\int_{H_{k}}^{\infty} C_{n}^{2} (z) d z - \int_{H_{k + 1}}^{\infty} C_{n}^{2} (z) d z] + 0.500 [\frac{\int_{0}^{H_{k}} C_{n}^{2} (z) z^{5 / 3} d z}{{H_{k}}^{5 / 3}} - \frac{\int_{0}^{H_{k + 1}} C_{n}^{2} (z) z^{5 / 3} d z}{{H_{k + 1}}^{5 / 3}}] - 0.452 [\frac{\int_{0}^{H_{k}} C_{n}^{2} (z) z^{2} d z}{{H_{k}}^{2}} - \frac{\int_{0}^{H_{k + 1}} C_{n}^{2} (z) z^{2} d z}{{H_{k + 1}}^{2}}]})^{- 3 / 5},

where

k

and

k + 1

represent the two adjacent measurements from the dynamically ranged Rayleigh beacon system. This differential metric accounts for errors in focal anisoplanatism associated between subsequent beacon measurements and can flow back into the treatment shown in Eq. (8) to yield

Eq. (13)

Δ σ_{FA}^{2} = {(\frac{D}{Δ d_{0}})}^{5 / 3} .

It is shown in Eq. (12) that the focal anisoplanatism error has a dependence on the integrated turbulence strength and beacon range. Without a priori knowledge of the turbulence strength profile, it is difficult to calculate the effect that the focal anisoplanatism error will have on the estimation accuracy of the segmented $C_{n}^{2}$ metric proposed by Eq. (5). However, trends can be elucidated through ratio analysis. For subsequent dynamically ranged beacon locations that are in close spatial proximity, the overall effect of focal anisoplanatism error is small. The backscattered field from the two beacons would experience nearly the same turbulence. If the distance between subsequent beacon measurement locations is large, the effect of focal anisoplanatism error could be relatively large. This could result in very different turbulence experienced by the backscattered fields from the two beacons and consequently influence the accuracy of estimating the segmented turbulence strength along the path as described by Eq. (5). To mitigate the negative effects of focal anisoplanatism, it is recommended that in the CONOPS planning for a data collect utilizing a dynamically ranged Rayleigh beacon the location choices for each individual beacon measurement be spaced in relatively close proximity to each other.

3.4.

Determination of Turbulence Profile ( $C_{n}^{2}$ )

The refractive index structure parameter is estimated using the formulation shown in Eq. (5) from data generated by the individual realization sequences that comprised a single random number generation seed. The $r_{0}$ data used are represented by the data shown in Fig. 4. This resulted in values of the refractive index structure parameter that were compared to the initial data input in the simulation environment. The results are shown in Fig. 5. Figure 5(a) shows a Hufnagel-Valley 5/7 profile and Fig. 5(b) shows a Hufnagel-Valley 5/7 with two enhanced turbulence strength peaks. For the Hufnagel-Valley 5/7 profile, the estimated $C_{n}^{2}$ profile data values have a mean absolute error of $4.00 \times 10^{- 19}$ from the input data values. The mean of the estimated values is 0.11 standard deviations away from the truth input corresponding to a mean percent difference of 0.58%. These metrics show that the data processing algorithm agrees well with the input parameters.

Fig. 5

Ensemble of $C_{n}^{2}$ estimates (data points) compared to truth input data (solid line). (a) The standard Hufnagel-Valley 5/7 profile and (b) the Hufnagel-Valley 5/7 profile with enhanced turbulence strength localized in two altitude regions.

4. Discussion

Utilization of individualized Fried parameter measurements in sequence to produce a discretized estimation of the turbulence strength profile is shown to produce quantitatively similar results in simulation to the inputted truth data. Analyzing the total ensemble of realizations, the $C_{n}^{2}$ estimates were only a fraction of a mean percent difference from the truth. This was expected based on the mathematical foundation of the profiling algorithm coupled with the small error that can propagate into the system from discretely modeling the wave propagations. These margins may grow for an operational system as larger errors in estimation may propagate through the system. Two examples of error creation sources are SHWFS centroiding imperfections in the presence of higher noise and unequal total SNR seen in the SHWFS spots for varied altitude range measurements.

The presence of localized turbulence strength spikes was realizable utilizing this framework as shown by the example case used in this paper. This example case was specifically designed with two high-strength localized turbulence regions to test the proposed algorithm’s ability to both isolate the strength spikes and estimate their strength increases. This was achievable through simulation since prior knowledge of the location of the spikes was known and the dynamically ranged beacon altitude locations for subsequent wave propagations could be chosen such that the sampling and resultant profile resolution was adequate to capture the turbulence strength spike. In practice, there could be practical limitations such as not knowing the location, thickness, or strength of a spike to capture the locational cut-off edges in range from the collecting telescope. Additionally, having a finite number of subsequent measurements achievable due to repetition rate limitations of the laser system coupled with the changing turbulence atmosphere time frames could bring the system outside of the frozen flow assumption. For example, if there are only 10 laser pulses available before the atmosphere changes to create a beacon source, then the measurements used will be few in number and will only be able to interrogate a fraction of the whole path in high resolution. To still have utility, the system could be configured so that long benign paths could be set as a single integrated strength measurement since this methodology relies on the $r_{0}$ metric, and the part of the path where the profile structure exists could be configured to obtain closely spaced beacon measurements. This is a CONOPS choice an operator of a dynamically ranged Rayleigh beacon system would have to make. With proper CONOPS choices and small prior knowledge assumptions, it is possible to effectively utilize a low number of measurements in an efficient way to make measurement-based estimates of the turbulence strength profile utilizing the algorithm presented in this paper.

Furthermore, an algorithm like the one presented that is reliant on assessing the Fried parameter for the metric of merit may not be the most optimal solution, although it produces agreeable results. A dynamically ranged Rayleigh beacon system that utilizes wavefronts measured across the entire pupil of the imaging system at a resolution smaller than $r_{0}$ that also has an aperture that is many $r_{0}$ s across inherently has spatial diversity in the tilt measurements. One could theorize that a hybrid sensing algorithm could be employed that utilizes the differential-range measurements such as those presented here enhanced by the use of the geometries of the spatial correlations of close and far neighboring tilt measurements on the SHWFS. The enhancement could come from treatment that is similar to that of a DIM LiDAR,¹³ but acting on the Shack–Hartmann slope measurements from differential range-stacked beacons. Alternatively, SLODAR-like processing could be employed, but the separation used to build up a measurement-based weighting function would come from subaperture separation and correlation from different range beacons. To accomplish this, a new derivation of cross-path weighting functions would need to be developed and the consequences of beacon location choices would need to be considered. These methods could have an advantage of localized high-resolution estimates of the $C_{n}^{2}$ parameter while also sweeping through a broad range of gates. The details of employing these types of algorithms are under development and will be presented in a future publication.

5. Conclusions

Atmospheric turbulence is known to be the limiting cause of image blur in many optical imaging systems. Consequently, a body of research has evolved over the years with ways to measure turbulence and mitigate its effects. As the imaging systems grow in resolution capability, overcoming turbulence effects will become all the more important. A crucial step to that will be understanding the profiled evolution of the turbulence along the viewing path. With this knowledge, advanced mitigation techniques and prediction techniques could be employed. To support future system capability trends, this paper presented an algorithm that could be applied to dynamically ranged Rayleigh beacon data to extract information about the strength changes along the propagation path. The modeling framework was created to investigate vignette scenarios that are represented of possible turbulence environments. Results were shown that supported correct estimation of turbulence strength profiles and identification of localized strength spikes under properly configured range gated set ups. In conclusion, the algorithm presented is able to adequately estimate the strength profile of optical turbulence along the viewing path utilizing wavefront data obtained from a dynamically ranged Rayleigh beacon-based system.

Acknowledgments

This work was supported by the Air Force Research Labs Sensors Directorate.

References

1.

S. M. Zuraski et al., “Electro-optic testbed utilizing a dynamic range gated Rayleigh beacon for atmospheric turbulence profiling,” Proc. SPIE, 10002 1000207 (2016). https://doi.org/10.1117/12.2240980 PSISDG 0277-786X Google Scholar

2.

M. C. Roggemann et al., “Sensing three-dimensional index-of-refraction variations by means of optical wavefront sensor measurements and tomographic reconstruction,” Opt. Eng., 34 (5), 1374 –1385 (1995). https://doi.org/10.1117/12.199870 Google Scholar

3.

M. Lloyd-Hart et al., “Dynamically refocused Rayleigh laser beacons for atmospheric tomography,” Proc. SPIE, 4494 259 –270 (2002). https://doi.org/10.1117/12.454800 PSISDG 0277-786X Google Scholar

4.

Y. Zhang et al., “Compensation of laser beam projection through turbulence in multimirror adaptive optics system with phase retrieval method,” Proc. SPIE, 5160 98 –106 (2004). https://doi.org/10.1117/12.502892 PSISDG 0277-786X Google Scholar

5.

A. Cortés et al., “Atmospheric turbulence profiling using multiple laser star wavefront sensors,” Mon. Not. R. Astron. Soc., 427 (3), 2089 –2099 (2012). https://doi.org/10.1111/j.1365-2966.2012.22076.x MNRAA4 0035-8711 Google Scholar

6.

D. J. Laidlaw et al., “Turbulence profiling for adaptive optics tomographic reconstructors,” Proc. SPIE, 9909 99093I (2016). https://doi.org/10.1117/12.2233323 PSISDG 0277-786X Google Scholar

7.

A. Lambert et al., “Turbulence profiling using extended objects for slope detection and ranging (SLODAR),” Proc. SPIE, 6316 631601 (2006). https://doi.org/10.1117/12.682428 PSISDG 0277-786X Google Scholar

8.

A. de la Nuez Cruz et al., “Atmospheric turbulence profiling at the Teide Observatory: comparison and calibration of SODAR and SCIDAR measurements,” Proc. SPIE, 6747 67470T (2007). https://doi.org/10.1117/12.737749 PSISDG 0277-786X Google Scholar

9.

J. Osborn, “Scintillation correction for astronomical photometry on large and extremely large telescopes with tomographic atmospheric reconstruction,” Mon. Not. R. Astron. Soc., 446 (2), 1305 –1311 (2015). https://doi.org/10.1093/mnras/stu2175 MNRAA4 0035-8711 Google Scholar

10.

J. Osborn et al., “Characterising atmospheric optical turbulence using stereo-SCIDAR,” J. Phys. Conf. Ser., 595 (1), 012022 (2015). https://doi.org/10.1088/1742-6596/595/1/012022 JPCSDZ 1742-6588 Google Scholar

11.

A. Tokovinin, “Turbulence profiles from the scintillation of stars, planets and moon,” in Workshop Astronomical Site Evaluation, (2007). Google Scholar

12.

A. Tokovinin et al., “Accurate seeing measurements with MASS and DIMM,” Mon. Not. R. Astron. Soc., 381 (3), 1179 –1189 (2007). https://doi.org/10.1111/j.1365-2966.2007.12307.x MNRAA4 0035-8711 Google Scholar

13.

G. G. Gimmestad et al., “LIDAR system for monitoring turbulence profiles,” Proc. SPIE, 6551 65510H (2006). https://doi.org/10.1117/12.720805 PSISDG 0277-786X Google Scholar

14.

T. J. Brennan, “Anatomy of the slope discrepancy structure function: characterization of turbulence,” Proc. SPIE, 5087 103 –114 (2003). https://doi.org/10.1117/12.487749 PSISDG 0277-786X Google Scholar

15.

S. K. Saha, Diffraction-Limited Imaging with Large and Moderate Telescopes, World Scientific Publishing Co. Pvt. Ltd., Hackensack, New Jersey (2007). Google Scholar

16.

C. W. Higgins et al., “The effect of scale on the applicability of Taylor’s frozen turbulence hypothesis in the atmospheric boundary layer,” Boundary Layer Meteorol., 143 (2), 379 –391 (2012). https://doi.org/10.1007/s10546-012-9701-1 Google Scholar

17.

L. C. Andrews and R. L. Phillips, Laser Beam Propagation through Random Media, SPIE Press, Bellingham, Washington (2005). Google Scholar

18.

J. D. Schmidt, Numerical Simulation of Optical Wave Propagation with Examples in MATLAB, SPIE Press, Bellingham, Washington (2010). Google Scholar

19.

S. M. Zuraski et al., “Turbulence and aerosol research dynamic interrogation system testing,” in Adv. Maui Opt. and Space Surveillance Technol. Conf., (2018). Google Scholar

20.

J. W. Hardy, Adaptive Optics for Astronomical Telescopes, Oxford University Press, Inc., New York (1998). Google Scholar

21.

D. L. Fried and J. F. Belsher, “Analysis of fundamental limits to artificial-guide-star adaptive-optics-system performance for astronomical imaging,” J. Opt. Soc. Am. A, 11 277 –287 (1994). https://doi.org/10.1364/JOSAA.11.000277 JOAOD6 0740-3232 Google Scholar

Biography

Steven M. Zuraski received his BS degree from Michigan State University in 2010, his MS degree in electro-optics from the University of Dayton in 2014, and is currently a PhD candidate in engineering physics at the Air Force Institute of Technology. He is a researcher at the Air Force Research Labs Sensors Directorate. His research interests include atmospheric turbulence sensing and electro-optic systems for remote sensing.

Elizabeth Beecher received her doctorate degree in electrical engineering from The Ohio State University. She is a research engineer at the Air Force Research Laboratory. She is a research team lead, manages multiple contracts, and mentors junior personnel. Her technical areas of expertise/interest include space domain awareness, small satellites and cube satellites, atmospheric turbulence, remote sensing, electro-optics, signal and image processing, space sensing, electromagnetics, communications, quantum physics, and autonomous optimization.

Jack E. McCrae received his BS degree in physics from Massachusetts Institute of Technology in 1984, his MS degree in physics (optics) from Air Force Institute of Technology (AFIT) in 1993, and his PhD in physics from AFIT in 1997. He is a retired Air Force Colonel currently at the Center for Directed Energy and a research assistant professor in the Engineering Physics Department at AFIT. His research interests include optics, lasers, quantum and nonlinear optics, laser radar, and atmospheric propagation.

Steven T. Fiorino received his BS degrees from The Ohio State in 1987 and Florida State in 1989, his MS degree in atmospheric dynamics from The Ohio State in 1993, and his PhD in physical meteorology from Florida State in 2002. He is a professor of atmospheric physics at AFIT and a director at AFIT Center for Directed Energy. His research interests include microwave remote sensing, weather signal processing, and atmospheric effects on weapon systems. He is a member of SPIE, OSA, AMS, AIAA, and DEPS.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Steven M. Zuraski, Elizabeth Beecher, Jack E. McCrae, and Steven T. Fiorino "Turbulence profiling using pupil plane wavefront data derived Fried parameter values for a dynamically ranged Rayleigh beacon," Optical Engineering 59(8), 081807 (14 May 2020). https://doi.org/10.1117/1.OE.59.8.081807

Received: 18 November 2019; Accepted: 29 April 2020; Published: 14 May 2020

Access the abstract

JOURNAL ARTICLE
14 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 7 scholarly publications.

Explore citations on Lens.org

KEYWORDS

Turbulence

Wavefronts

Rayleigh guide stars

Atmospheric propagation

Profiling

Wave propagation

Atmospheric turbulence

1.

Introduction

2.

Key Concepts to Support Algorithm Development and Evaluation

2.1.

Optical Turbulence Metrics

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

Eq. (5)

2.2.

Simulated Wave Propagation

2.3.

Shack–Hartmann Wavefront Sensor Measurements

2.4.

Profiled Turbulence Metrics

Fig. 1

3.

Algorithm Implementation and Evaluation

3.1.

Overview of Simulation Implemented

3.2.

Zonal Measurements to Estimate r0

Eq. (6)

Eq. (7)

Fig. 2

Fig. 3

Fig. 4

3.3.

Influence of Focal Anisoplanatism

Eq. (8)

Eq. (9)

Eq. (10)

Eq. (11)

Eq. (12)

Eq. (13)

3.4.

Determination of Turbulence Profile (Cn2)

Fig. 5

4.

Discussion

5.

Conclusions

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years

Zonal Measurements to Estimate $r_{0}$

Determination of Turbulence Profile ( $C_{n}^{2}$ )