Characterization of Localized Atmospheric Turbulence Layer Using Laser Light Backscattered o ﬀ Moving Target

: A concept of atmospheric turbulence characterization using laser light backscattered o ﬀ a moving unresolved target or a moving target with a glint is considered and analyzed through wave-optics numerical simulations. The technique is based on analysis of the autocorrelation function and variance of the power signal measured by the target-in-the-loop atmospheric sensing (TILAS) system composed of a single-mode-ﬁber-based optical transceiver and the moving target. It is shown that the TILAS received power signal autocorrelation function strongly depends on the turbulence distribution and is weakly sensitive to the turbulence strength, while the signal variance equally depends on these parameters. Assuming the atmospheric turbulence model can be represented by a single spatially localized turbulence layer and the target position and speed are known independently, consecutive analysis of the autocorrelation function and variance of the TILAS signal allows evaluation of both the turbulence layer strength and position along the optical propagation path. It is also demonstrated that the autocorrelation function can potentially be used for the atmospheric turbulence outer scale estimation. Author Contributions: Conceptualization, V.A.K., M.A.V., and V.S.R.G.; Methodology, V.A.K. and S.L.L.; Software, V.A.K.; Validation, V.A.K.; Formal analysis, V.A.K., S.L.L., and M.A.V.; Investigation, V.A.K., S.L.L., M.A.V., and V.S.R.G.; Writing—original draft preparation, V.A.K. and S.L.L.; Writing—review and editing, V.A.K., S.L.L., M.A.V., V.S.R.G. the


Introduction
Atmospheric turbulence characterization based on remote sensing of the target-plane laser beam characteristics is typically performed along the optical propagation path between optical transceiver and a target. Assuming the turbulence is uniform (homogeneous) and described by the Kolmogorov power spectrum [1], these characteristics can then be used to estimate the refractive index structure parameter C 2 n averaged over the target line of sight. However, in many cases, atmospheric turbulence is non-uniform [2] and composed of relatively thin, horizontally extended layers of high values of the refractive index caused, e.g., by evening transition of the boundary layer [3] or jet streams and gravity waves in the upper troposphere and lower stratosphere [4]. These layers have frequently been observed with balloons [5] and radars [6]. The knowledge of atmospheric turbulence layer characteristics, such as their strength and position along the optical propagation path, is important for improving the efficiency of laser systems in beam projection and free-space optical communication applications using, e.g., multi-conjugate adaptive optics technique [7].
Appl. Sci. 2020, 10, 6887; doi:10.3390/app10196887 www.mdpi.com/journal/applsci In this paper, we present a concept which allows characterization of atmospheric turbulence layers along the trajectory of a moving target as well as real-time evaluation of major laser beam characteristics such as intensity scintillations directly on the remote target. The technique is based on the previously developed target-in-the-loop atmospheric sensing (TILAS) concept [8], which was experimentally proven over a 7-km atmospheric propagation path to a stationary retro-target [9]. Here, we consider a more general case when the TILAS system, a single-mode-fiber-based optical transceiver, is operating with a moving target. As shown in the paper, the use of the moving target offers several advantages for atmospheric characterization such as estimation of strength and position of a spatially localized turbulence layer, as well as evaluation of the turbulence outer scale L 0 , based on analysis of both the autocorrelation function and variance of the power signal measured by the TILAS system.
The principal assumptions required for the TILAS technique presented here are as follows. First, the target is supposed to either be unresolved (smaller than the diffraction-limited beam footprint at the target plane), have a glint caused by reflection from the target's surface, or have a small-size retro-reflector attached. It is also presumed that the distance to the target as well as the target's speed are known. Besides, we utilize the Taylors' frozen turbulence hypothesis and assume that the target is moving sufficiently fast, so that we may neglect the impact of the wind along the propagation path.
The paper is organized as follows. A concept of the TILAS system operating with a moving target with a glint is described in Section 2. Numerical simulation approach for modeling and simulation of the TILAS system is detailed in Section 3. Numerical simulation results are discussed in Section 4. The paper concludes with a summary and suggestions for future work.

Target-in-the-Loop Atmospheric Sensing Concept
The theoretical foundation of the TILAS concept is based on the optical reciprocity principle [10] and the associated conservation law for the so-called interference metric J int coupling the complex amplitudes of the counter-propagating optical waves [11][12][13]. As was recently shown [8,9,14], the square modulus of the interference metric can be measured using the TILAS system depicted in Figure 1a. This system is composed of a laser, fiber-optics circulator, fiber-coupled transceiver, and fiber-coupled photodetector. The circulator separates the outgoing and incoming laser beams by sending the outgoing, linearly-polarized Gaussian beam from the laser to the transceiver, while the entering wave is directed to the photodetector. The power signal measured at the TILAS photodetector is given by the expression [8] where α is a coefficient dependent on the transmitted power and the photodetector's sensitivity, A(r, z = 0) is the stationary complex amplitude of the outgoing wave transmitted along the optical z-axis, ψ(r, z = 0, t) is the amplitude of the incoming wave backscattered off a target located at the plane z = L back to the receiver (transceiver) plane z = 0, r = {x, y} is the spatial vector in the plane orthogonal to the propagation direction, t is the time, and the integration is performed over the entire transverse plane. In this paper, we consider the stationary TILAS transceiver and the target with a glint, moving perpendicularly to the sensing path direction with velocity T v ( Figure 1b). For this type of target, the following quantities are computed based on the TILAS received power signal measurements: 1. Autocorrelation function of the signal, where τ is the time delay between the power signal values measured at different moments of time t and angle brackets denote the averaging over the integration time interval T. Here we assume that at T τ , where τat is the characteristic correlation time of the atmospheric turbulence pattern change inside the transceiver aperture [8].
2. Interference metric variance, also referred to here as the TILAS received power signal variance, given by the expression [8] In this paper, we focus on the atmospheric turbulence model that can be represented by a single spatially localized turbulence layer, referred to as here as the single layer model. The homogeneously distributed turbulence model is also briefly considered. We show that parameters B(τ) and 2 J σ given by the equations above can be used for estimation of both the strength and position of the turbulence layer along the optical propagation path, and, potentially, can lead to the turbulence profiling. We also demonstrate that consecutive analysis of the correlation properties of the TILAS received power signal PTILAS (t) results in the estimation of the atmospheric turbulence outer scale L0 of the layer.
Numerical simulation approach for the analysis is explained in the following section.

Numerical Simulation Approach and Settings
Our approach is based on direct numerical integration of the parabolic wave equation (the socalled wave-optics technique). This integration is performed using the split-step operator method [15][16][17], where the turbulence-induced refractive index perturbations along the propagation path of length L are represented by a set of thin, two-dimensional (2D), statistically independent phase screens.
In our simulations, the number of phase screens M is varied from one for the single layer model to five for the distributed turbulence model. For the single layer model, the screens are positioned at five selected distances from the TILAS transceiver, while for the distributed turbulence model, they In this paper, we consider the stationary TILAS transceiver and the target with a glint, moving perpendicularly to the sensing path direction with velocity v T (Figure 1b). For this type of target, the following quantities are computed based on the TILAS received power signal measurements:

1.
Autocorrelation function of the signal, where τ is the time delay between the power signal values measured at different moments of time t and angle brackets denote the averaging over the integration time interval T. Here we assume that T τ at , where τ at is the characteristic correlation time of the atmospheric turbulence pattern change inside the transceiver aperture [8].

2.
Interference metric variance, also referred to here as the TILAS received power signal variance, given by the expression [8] In this paper, we focus on the atmospheric turbulence model that can be represented by a single spatially localized turbulence layer, referred to as here as the single layer model. The homogeneously distributed turbulence model is also briefly considered. We show that parameters B(τ) and σ 2 J given by the equations above can be used for estimation of both the strength and position of the turbulence layer along the optical propagation path, and, potentially, can lead to the turbulence profiling. We also demonstrate that consecutive analysis of the correlation properties of the TILAS received power signal P TILAS (t) results in the estimation of the atmospheric turbulence outer scale L 0 of the layer.
Numerical simulation approach for the analysis is explained in the following section.

Numerical Simulation Approach and Settings
Our approach is based on direct numerical integration of the parabolic wave equation (the so-called wave-optics technique). This integration is performed using the split-step operator method [15][16][17], where the turbulence-induced refractive index perturbations along the propagation path of length L are represented by a set of thin, two-dimensional (2D), statistically independent phase screens.
In our simulations, the number of phase screens M is varied from one for the single layer model to five for the distributed turbulence model. For the single layer model, the screens are positioned at five selected distances from the TILAS transceiver, while for the distributed turbulence model, they are distributed equidistantly along the propagation path at distances z j = ∆z(2 j − 1)/2, where ∆z = L/M is the layer width and j = 1, . . . , M. The latter case is illustrated in Figure 2 for M = 5.
are distributed equidistantly along the propagation path at distances (2 1) 2 is the layer width and j = 1, …, M. The latter case is illustrated in Figure 2 for M = 5. For the distributed turbulence model, we assume that the turbulence is homogeneous, which is a commonly used assumption for relatively short or nearly horizontal propagation paths [18]. In this case, 2 ( ) const n C z = along the entire path and To model engagements involving moving targets, the technique for generation of infinitely-long phase screens [20] is used. The phase screens are assumed to obey the von Karman power spectrum law; it allows us to vary the value of the turbulence outer scale L0. Note that the use of infinitely-long phase screens permits extending the range of L0 values beyond the numerical simulation area size (up to four times the size). The turbulence strength is defined by the refractive index structure parameter C 2 n (z); the corresponding plane-wave Fried parameter computed over the distance L is given by [18], where k = 2π/λ is the wavenumber. For the phase screen located at z = z j , the Fried parameter is correspondingly defined as r For the single layer model, we suppose that C 2 n (z) = const withing the layer of width ∆z and zero otherwise, i.e., r 0 = 1.68(k 2 C 2 n ∆z) For the distributed turbulence model, we assume that the turbulence is homogeneous, which is a commonly used assumption for relatively short or nearly horizontal propagation paths [18]. In this case, C 2 n (z) = const along the entire path and For fair comparison of the results obtained for the single layer model with those for the distributed turbulence model, which corresponds to the same value of r 0 for both models, the refractive index structure parameter values for the single layer model are M-fold higher than those for the distributed turbulence model.
In our simulations of the single layer model, the layer width is fixed at ∆z = L/5 and the turbulence strength C 2 n is varied in the interval 10 −16 m −2/3 ≤ C 2 n ≤ 10 −14 m −2/3 . The Fried parameter is respectively ranging from 65 cm to 4.1 cm. For the distributed turbulence model, the corresponding C 2 n values are decreased five times.
To simulate a target moving with the speed v T = |v T | relative to a stationary laser platform with the TILAS sensor, the position of the target is fixed but the jth phase screen is shifted in the opposite direction with the speed v j = v T z j /L. Over the time interval ∆t = t − t , this corresponds to the shift distance ∆s j = v j ∆t, as illustrated in Figure 2. Thus, during the same time, screens located at different positions along the path are shifted different distances.
To model engagements involving moving targets, the technique for generation of infinitely-long phase screens [20] is used. The phase screens are assumed to obey the von Karman power spectrum law; it allows us to vary the value of the turbulence outer scale L 0 . Note that the use of infinitely-long phase screens permits extending the range of L 0 values beyond the numerical simulation area size (up to four times the size).
As the TILAS outgoing beam, a collimated Gaussian beam with fixed wavelength λ = 1.064 µm and radius W 0 = 1.5 cm, truncated by the circular aperture of diameters d TILAS = 3.3 cm and 7 cm is considered. The propagation distance to the target is fixed at L = 7 km. The target's glint is modeled by the circular retro-reflector of diameter d T = 1 cm. Note that for the set of the parameters selected, d T is smaller than the diffraction-limited target-plane beam footprint radius 16 cm, i.e., the retro-reflector is unresolved. The target speed relative to the TILAS system platform with the TILAS sensor is set to v T = 200 m/s. All computations are performed using numerical grids with resolution of 2048 × 2048 pixels; the simulation area width and sampling size are fixed at 4.096 m and 2 mm, respectively, to meet sampling requirements [21]. The sampling time, i.e., the time between the consecutive phase screen shifts, is fixed at 0.1 ms. The integration interval T is chosen to be equal to 5 s. During this time, the target moves over a distance of 1 km; this results in about 70-m propagation distance change, which is insignificant compared to L = 7 km and can be neglected. At the same time, as mentioned earlier, this interval should exceed the characteristic correlation time τ at , which can be estimated for the jth screen as τ 0 is the jth screen's phase fluctuations correlation length, also known as the Tatarskii parameter for spherical wave [22]. Assuming the turbulence layer of width ∆z = 1.4 km is located at z = 0.1 L = 0.7 km (closest to the TILAS transceiver aperture) and moving with the speed v = 0.1v T = 20 m/s, τ at ≈ 28 ms for C 2 n = 1 × 10 −16 m −2/3 and 1.7 ms for C 2 n = 1 × 10 −14 m −2/3 , which means that the condition T τ at is satisfied. The results of our numerical simulations are presented and discussed in the following section.

Correlation Properties of TILAS Received Power Signal
Consider first the correlation properties of the TILAS received power signal P TILAS (t) for the single layer model, where the turbulence layer of a fixed strength is placed at different distances from the TILAS transceiver. Figure 3a shows the time series of P TILAS (t) computed over a 10-ms time interval with 0.1-ms sampling time for two different relative positions of the layer between the TILAS transceiver and the moving target, z/L = 0.1 (closer to TILAS) and z/L = 0.9 (closer to the target), characterized by the same Fried radius r 0 = 4.1 cm and turbulence outer scale L 0 = 4 m. As can be seen from Figure 3a, the signal obtained for the case of the turbulence layer located closer to the transceiver (blue line) is changing slower than that for the layer located near the target (sky-blue line). This can be explained by the fact that the turbulence-induced refractive index inhomogeneities inside the propagating beam footprint are updated faster when the layer is located farther from the TILAS transceiver aperture, or closer to the moving target. This happens because during the same time interval, the laser beam following the target would scan sections of different length for layers positioned at different distances from TILAS (see Figure 2).
Typical dependencies of the autocorrelation functions B (Equation (2)) on the time delay τ computed for five different turbulence layer positions ranging from z = 0.1 L to z = 0.9 L with the step of 0.2 L are shown in Figure 3b. Note that the autocorrelation functions B(τ) have Gaussian-like shape, which allows us to directly determine their width by the drop in the function level by a factor exp(−1) ≈ 0.37 compared with B(τ = 0). This width is referred to here as the correlation time, τ cor . As seen, the correlation time depends on the position of the turbulence layer along the propagation path from the TILAS transceiver to the moving target. The autocorrelation function corresponded to the turbulence layer located near the moving target is characterized by the smallest correlation time while B(τ) related to the layer near the transceiver has the largest τ cor . Typical dependencies of the autocorrelation functions B (Equation (2)) on the time delay τ computed for five different turbulence layer positions ranging from z = 0.1 L to z = 0.9 L with the step of 0.2 L are shown in Figure 3b. Note that the autocorrelation functions B(τ) have Gaussian-like shape, which allows us to directly determine their width by the drop in the function level by a factor exp(−1) ≈ 0.37 compared with B(τ = 0). This width is referred to here as the correlation time, τcor. As seen, the correlation time depends on the position of the turbulence layer along the propagation path from the TILAS transceiver to the moving target. The autocorrelation function corresponded to the turbulence layer located near the moving target is characterized by the smallest correlation time while B(τ) related to the layer near the transceiver has the largest τcor.
Consider now the distributed turbulence model characterized by the same Fried radius equal to 4.1 cm ( 2 defined using the single layer model. However, the distributed turbulence autocorrelation function has a unique shape with a specific slow-decline region at the function's tail (a 'pedestal' or a 'base'), which is distinct from the Gaussian-like shape of the autocorrelation function computed for the single layer model. This means that the presence of such a base can serve as an indication of the turbulence being distributed, or at least of the presence of more than one turbulence layer along the propagation path.   Figure 3b by the dash-dotted black line. As can be seen from the figure, the distributed turbulence correlation time defined at B(τ) ≈ 0.37 lies between τ cor (z/L = 0.7) and τ cor (z/L = 0.5) defined using the single layer model. However, the distributed turbulence autocorrelation function has a unique shape with a specific slow-decline region at the function's tail (a 'pedestal' or a 'base'), which is distinct from the Gaussian-like shape of the autocorrelation function computed for the single layer model. This means that the presence of such a base can serve as an indication of the turbulence being distributed, or at least of the presence of more than one turbulence layer along the propagation path.  Figure 4b. Further analysis reveals that the τ cor curves stay qualitatively the same for different parameters of both the TILAS system and the moving target. Estimations of τ cor performed with the increased diameter of the TILAS transceiver aperture, d TILAS = 7 cm, demonstrate comparable monotonically decreasing dependencies τ cor (z/L) as in Figure 4a, although the absolute values of τ cor are generally larger for larger d TILAS . Similar observations are made for the target speed v T ; τ cor is simply inversely proportional to v T . Thus, assuming both the TILAS and the target parameters including the distance to the target and target speed are known, the correlation time τ cor , which is computed using the autocorrelation function B(τ) of the received power signal P TILAS (t) measured by the TILAS system, can be used to estimate the position of the turbulence layer. The accuracy of this estimation is somewhat reduced by the dependence of τ cor on the turbulence strength.  Further analysis reveals that the τcor curves stay qualitatively the same for different parameters of both the TILAS system and the moving target. Estimations of τcor performed with the increased diameter of the TILAS transceiver aperture, dTILAS = 7 cm, demonstrate comparable monotonically decreasing dependencies τcor (z/L) as in Figure 4a, although the absolute values of τcor are generally larger for larger dTILAS. Similar observations are made for the target speed T v ; τcor is simply inversely proportional to T v . Thus, assuming both the TILAS and the target parameters including the distance to the target and target speed are known, the correlation time τcor, which is computed using the autocorrelation function B(τ) of the received power signal PTILAS (t) measured by the TILAS system, can be used to estimate the position of the turbulence layer. The accuracy of this estimation is somewhat reduced by the dependence of τcor on the turbulence strength.

Estimation of Turbulence Layer Strength Based on TILAS Measurements
Consider analogous dependencies of the TILAS received power signal variance 2 J σ (Equation

Estimation of Turbulence Layer Strength Based on TILAS Measurements
Consider analogous dependencies of the TILAS received power signal variance σ 2 J (Equation (3))on the turbulence layer position (Figure 5a), as well as the turbulence layer strength (Figure 5b). As seen, σ 2 J demonstrates comparable level of variability for both z/L and C 2 n . Besides, the function σ 2 J (z/L) in Figure 5a is not monotonic; strongest signal fluctuations correspond to the turbulence layer placed in the middle of the propagation path. Therefore, the power signal variance is not suitable to be used to define the layer position. However, σ 2 J (C 2 n ) in Figure 5b is monotonically increasing with the increase of the turbulence strength (at least inside the examined interval of C 2 n ). At the same time, as was confirmed by our numerical analysis, this variance does not depend on the speed of the target, v T . Therefore, if the turbulence layer position z/L is already defined from the correlation times estimations and the TILAS parameters and the distance to the target are known independently, the set of functions σ 2 J (C 2 n ) computed for different layer positions can be used for retrieval of the turbulence layer strength.  Figure 5a). This suggests that the layers positioned close to the middle of the The value of the TILAS received power signal variance for the distributed turbulence case with C 2 n = 2 × 10 −15 m −2/3 and the same Fried parameter as for the single layer model with C 2 n = 10 −14 m −2/3 is close to those obtained for the turbulence layer placed around 0.3 z/L or 0.7 z/L (see dash-dotted and solid lines in Figure 5a). This suggests that the layers positioned close to the middle of the propagation path contribute the most to the power signal variance.

Evaluation of Turbulence Layer Outer Scale Based on TILAS Measurements
In this work, the possibility of direct evaluation of the atmospheric-turbulence outer scale L 0 using TILAS measurements and autocorrelation function analysis is investigated. This possibility can help to resolve a long-standing problem of sensing of the so-called piston phase (optical path difference) [23], which, within the framework of the Kolmogorov turbulence theory, is associated with the outer scale. This is essential, e.g., for long-range electro-optical systems based on coherent detection techniques [24].  Figure 6. Here, the TILAS transceiver aperture diameter is fixed at d TILAS = 7 cm; for the TILAS aperture size d TILAS = 3.3 cm considered in previous sections, the functions dependence on the outer scale is less evident. Analysis of shapes of these functions reveals that the B(τ) decline rate increases with the decrease of L 0 . Besides, if the outer scales are small enough, the negative correlation may be observed. The impact of L 0 seems to be more pronounced on the functions computed for layers located closer to the TILAS transceiver (compare Figure 6a for z/L = 0.1 (near TILAS) and Figure 6b for z/L = 0.3); for turbulence layers that are closer to the target, the width of the autocorrelation functions (i.e., the correlation time τ cor ) becomes nearly indifferent to the turbulence outer scale. To increase the sensitivity of τ cor to L 0 , the correlation time can be re-defined, e.g., by the drop in the intensity level by a factor exp(−2) ≈ 0.14. Alternatively, the size of the transceiver aperture diameter can be further increased; as mentioned in Section 4.2, the values of τ cor increase with the increase of d TILAS , which helps to differentiate the curves with different values of L 0 . Thus, the correlation time can be used to estimate the turbulence outer scale if the layer position and strength are already known.

Discussion and Conclusions
Dependence of the autocorrelation function B(τ) and the variance of the TILAS received power signal 2 J σ on the atmospheric turbulence distribution, strength, and outer scale can serve for turbulence characterization. Sequential analysis of these quantities provides information about the turbulence properties along the laser beam propagation path. Analysis of the shape of the autocorrelation function allows distinguishing between homogeneous and inhomogeneous turbulence distributions, while estimation of the strength of homogeneous turbulence or localized

Discussion and Conclusions
Dependence of the autocorrelation function B(τ) and the variance of the TILAS received power signal σ 2 J on the atmospheric turbulence distribution, strength, and outer scale can serve for turbulence characterization. Sequential analysis of these quantities provides information about the turbulence properties along the laser beam propagation path. Analysis of the shape of the autocorrelation function allows distinguishing between homogeneous and inhomogeneous turbulence distributions, while estimation of the strength of homogeneous turbulence or localized layer can be performed based on the signal variance when the turbulence distribution is known.
The procedure of finding the position and strength of a single turbulence layer based on the τ cor and σ 2 J quantities obtained using the TILAS received power signal measurements can be summarized by the following algorithm: Step 1. The relative turbulence layer position z/L is estimated using the set of functions τ cor (z/L) computed for given TILAS parameters, known distance to the target and target speed, as well as various turbulence strengths of the layer C 2 n , as in Figure 4a. Since the functions τ cor (z/L) are weakly dependent on the turbulence strength, the turbulence layer position will be estimated with some error due to the uncertainty of C 2 n .
Step 2. The turbulence layer strength is defined using the set of functions σ 2 J (C 2 n ) computed for the estimated interval of z/L values, as in Figure 5b. The error in the layer position will introduce the corresponding error in the turbulence layer strength estimates.
To reduce the estimation errors, the algorithm can be iterated using the obtained turbulence strength interval as the initial C 2 n interval in Step 1 until the desired accuracy is achieved. To speed up the turbulence layer position and strength retrieval procedure, a set of functions τ cor (z/L) and σ 2 J (C 2 n ) can be precomputed.
Three examples of the evaluation of the turbulence layer characteristics are given in Table 1. The first two columns show the input parameters used for simulations. The third and fourth columns contain the quantities computed using the TILAS power signal processing, which are used to estimate z/L and C 2 n of a single turbulence layer. The fifth and sixth columns show the intervals of values for the layer position and strength computed after the first iteration of the algorithm. Note that the identification of the layer position appears to be more accurate for layers located on the ends of the propagation path, while the estimation of the layer strength is more accurate for layers located in the middle of the optical path. The second iteration results in noticeable increase of the accuracy, as can be seen by comparing the estimates (last two columns) to the input parameters (first two columns). More rigorous analysis based on simultaneous evaluation of strength and positions of few layers performed with summarizing autocorrelation function with weight coefficients regarding their position and strength can potentially provide detailed information about complex turbulence distribution. Additional opportunities are related to the use of a multi-aperture TILAS transceiver system. Such a system can transmit several laser beams with different initial diameters, which would serve for estimation of the turbulence outer scale.