4.1

Turbulence and wind database

The profiles database has been provided by Durham University. Wind and C_n² profiles were obtained through a global turbulence model capable of converting global meteorological re-analysis models (in this case ERA5⁵) into 3-Dimensional optical turbulence maps. This model is similar to the one described in Osborn and Sarazin⁶ but was improved to include a separate boundary layer and enable stronger turbulence strength near the ground to be modelled⁷.

The ECMWF model from which the turbulence is being calculated has a spatial resolution of 0.3° along latitude and longitude and a temporal resolution of 1 hour. We chose a grid of 11 by 11 points around Tenerife’s island, which is a site of interest for a potential future ground station, and focused on the first 19 days of march 2018. It leaves us with 121 simultaneous measurements for each hour, with some missing values. On the overall 19 days considered it leaves us with 37059 profiles on 113 pressure levels each.

In the following we consider that these profiles are representative of field data and could be obtained with instruments directly measuring turbulence.

4.2

Power injected time series database and limitation

Correlated time series of received optical power in a single mode fiber are obtained from the C_n² and wind profile using SAOST (Simplified Adaptive Optics Simulation Tool)^8,9,10 a pseudo-analytic adaptive optics simulation tools.

This tool assumes decorrelation between phase and scintillation perturbations. Scintillation influence is simulated assuming the small perturbation approximation in the Rytov regime. In practice this constrains the validity of the approach to limited Rytov variance (for horizontal propagation, typically ). Correlated time series of ρ_φ, the coupling efficiency neglecting scintillation, and ρ_I, the term of scintillation are generated independently. The received optical power (ROP) attenuation is obtained through the product of these two quantities.

The residual phase error is computed in SAOST using a Monte-Carlo approach. Following the algorithm described in Roddier¹¹ random occurrences of Zernike coefficients are sampled to describe the corrected phase. The temporal correlation of ρ_φ is obtained by filtering the raw Zernike coefficients by their temporal power spectral model¹² in the Fourier domain. Some hypothesis have to be made on the turbulence condition of which we can cite the use of von Kármán wind turbulence model and a fixed outer scale. More details and comparisons with an end-to-end simulation can be found in Canuet¹³.

The temporal impact of the scintillation contribution ρ_I to the overall attenuation is approximated by:

in the weak fluctuation regime of the scintillation index χ(r, t) where is the variance of the punctual log amplitude and χ_AP(t) is the log-amplitude averaged on the receiver aperture.

Under weak irradiance fluctuation hypothesis, ¹⁴’¹⁵. Series are randomly generated following this normal law and the temporal correlation is obtained by filtering with the spectrum given in Robert¹⁶.

In the following all simulations have been made for a state of the art mid-range adaptive optics system. The number of corrected radial order is fixed to 10 and the AO frequency to 2kHz. The diameter of the ground station telescope is 60cm.

4.3

Choice of relevant integrated parameters

We have based our choice of integrated parameters on the integrated parameters regularly used to describe the error budget of an adaptive optics.

The first chosen parameter is Fried’s parameter¹⁷, a key parameter when one is interested in the effects of turbulence. Fried’s parameter r₀ is defined as the typical diameter of a telescope that would be limited by atmospheric turbulence. In the case of a plane wave considering a Kolmogorov spectrum, we have along the line of sight:

Where λ is the wavelength (1.55μm). Experimentally, Fried parameter can be estimated either from the amplitude jitter of a star at the focal plane of an imager or more robustly thanks to differential imaging such as performed with a DIMM¹⁸.

The second parameter, denoted , is a measure of the height dispersion of atmospheric layers, homogeneous to an altitude:

It provides an assessment of the physical origin for the angular decorrelation of the phase perturbations (the influence of distant turbulence layers) while being independent from the turbulence strength. It is related to the Fried parameter and the isoplanatic patch θ₀²⁰ by:

The estimation of can therefore be performed thanks to a measurement of θ₀ which is derived for instance from limited aperture averaged scintillation by night-time²¹ or thanks to a Shabar measurement by day-time²². The last parameter, denoted describes an average wind speed over the turbulent layers²³:

It is related to the turbulence coherence time τ according to:

The turbulence coherence time can be extracted for instance from the temporal analysis of the jitter of a bright enough point source image. This evaluation can then be exploited to provide an estimation of .

In the following the point source log-amplitude variance is supposed to be known (measured for instance from the scintillation of a bright point source with a small diameter instrument).

We report on figure 5 the distribution and cumulative distribution of the Fried parameter for all the profiles of the database. Same quantities for and are plotted figure 6 and figure 7 respectively. Point source log-amplitude variance distribution and cumulative distribution are plotted figure 8. Evaluations are calculated at 1.55 μm for a 33.2° elevation. Point source log-amplitude variance is very rarely above 0.3 which tends to consolidate the Rytov regime hypothesis.

Figure 5:

Distribution and cumulative distribution of r₀ (m).

Figure 6:

Distribution and cumulative distribution of (m) for all profiles.

Figure 7:

Distribution and cumulative distribution of (m.s^-1) for all profiles.

Figure 8:

Distribution and cumulative distribution of .

5.1

Method

Using the same notation as in 4.2, we want a machine-learning algorithm that would give us the probability density function (PDF) of ρ_φ and of ρ_I. In a second step, we will also want to estimate the half correlation time of ρ_φ and of ρ_I. A usual and effective way to describe a distribution using machine learning is to parameterize the PDF and then use a model to estimate the parameters. In 4.2, we gave the parametric form of ρ_I that depends only on one unknown parameter assuming the point source log-amplitude variance is measured. For ρ_φ, it is a bit more complex as there is no known parametric form so we conducted a statistical study to determine which density function best fits the data.

We were led to study , the loss in power induced by the phase fluctuation in dB as it gives a bigger weight to the low values of ρ_φ(t) that are the critical values for our application.

Studying the distribution of this loss on our 37k profiles, it appears that the tail of the distribution shows an exponential decay, which is consistent with the closed form of the distribution proposed by Canuet¹⁰. Among all light-tailed distributions, the one that statistically best fits our data is a Gumbel distribution of the following form:

where . Comparing real and fitted distribution quantiles, we show that for the 1% quantile the median relative error is 0.48% and never exceeds 2.01%, and the bigger the quantile the lower the error.

This result is particularly interesting as, with the Gumbel distribution being a good approximation of the density probability of phase related power attenuation we can describe the latest with only two parameters, μ and β and then describe power attenuation distribution with only 3 parameters :μ, β and σ_χAP(4 parameters if we include the point source log-amplitude variance which is supposed to be directly measured).

5.2

Result for probability density

Different machine learning metamodels have been tested of which we can cite gradient boosting, multilayer perceptron, support vector machine regression and Gaussian process²⁴. Best results were obtained for Gaussian process²⁵.

15% of our database was used for a training set and the 85% remaining for a test set. Using r₀, and as inputs of the model we were able to assess μ with a coefficient of determination of 0.9994, β with a coefficient of determination of 0.9992 and σ_χAP with 0.9969.

In order to characterize the goodness of our estimation on the received optical power’s statistic, the absolute error made on the mean and standard deviation of the reconstructed PDF of the received optical power are computed. We also look at the absolute error on the 1% quantile of the received optical power (see red dashed line on figure 9) as we want a faithful reproduction of the tail of the distribution. Statistics on the absolute error made on these values on our 37k profiles can be seen in figure 10. On all profiles, the prediction error on the value of the 1% quantile is less than 0.7 dB.

Figure 9:

Example of PDF computed from SAOST (blue) time series and PDF computed from our model (yellow) with, from left to right: L_φ, ρ_I and the received optical power.

Figure 10:

Histogram and CDF of the absolute error on prediction of the mean, the standard deviation and the 1% quantile.

5.3

Result for temporal content

Temporal correlation of the power attenuation will affect the duration on which interleaving shall be performed to mitigate fading losses. We investigate here the possibility to estimate the half correlation of L_φ and ρ_I using only r₀, and . Once again, best results were obtained with the use of a Gaussian process metamodel. On figure 11 where the estimated value is compared to the expected one, we can see that the coefficient of determination is 0.88 for the phase demi-correlation and 0.84 for the scintillation²⁵.

Figure 11:

Prediction of the half-correlation time using Gaussian process on inputs r₀, and , x axis shows the true value and y the predicted one. On the left the half-correlation neglecting scintillation effects and on the right neglecting phase effects.

This low score is due to the lack of temporal information found in r₀, and . One way to recover the missing temporal information is to add a measurement of the temporal power spectral density (PSD)of the scintillation to our inputs. As inputs of our model we added 5 moments of the PSD of scintillation measured by a 5 cm pupil instrument, these results were simulated using the formulas of Shen²⁶. As seen on figure 12, adding information on the PSD to our model enabled much better results with coefficient of determination of 0.99 for the phase and 0.98 for the scintillation²⁵.

Figure 12:

Prediction of the half-correlation time using Gaussian process on inputs r₀, , and information on the PSD of scintillation for a 5cm pupil. Same other parameters as for figure 11.

6.1

Inputs and limitation

The main challenge in the estimation of the power margin is to describe the impact of using an interleaver. Indeed, the interleaver acts as a convolutional filter on the mutual information time series seen by the receiver and this non-linear operation cannot be translated in simple terms by knowing only the half correlation and the probability density. In addition, we are looking for a valid method regardless of the size of the interleaver.

The best way to achieve our goal is to generate a stochastic process equivalent to any time series of ROP obtained with the same exact atmospheric and instrumental conditions. By being stochastically equivalent to the ground truth, the generated series is sharing the same finite-dimensional law and thus the same power margin.

We know that a Gaussian process is defined entirely by its auto covariance function²⁷. However, this is not necessarily true for all types of stochastic processes and there is no reason it should be the case for our process, product of two non-Gaussian stochastic processes (phase and scintillation).

6.2

Downlink budget excluding turbulence

The optical downlink signal between the on-board terminal described in Berceau et al.²⁸ and the ESA OGS described section 3 is considered. The communication chain equipment (Laser Communication Electronic, Modems and Laser Power Electronics) considered are the one of Poulenard et al.¹.

6.3

Method for link margin estimation

Even though it is very easy to generate a Gaussian process, there are very few methods for generation of stochastic process with any given marginal distribution. An intuitive method was developed by Tsoukalas²⁹ and comes with a R library anySim³⁰. It is based on Nataf’s joint distribution model where the distribution of a correlated random variable is mapped to a Gaussian distribution via its inverse cumulative distribution function.

Using this library and the previously estimated distribution and half-correlation, we generated time series of coupled flux with phase effects only and with scintillation effects only. Combining the two series with the downlink budget given in section 6.2, we obtained time series of ROP using no prior knowledge on the profiles but r₀, , and .

AnySim requires an autocorrelation function; the simplest interpolation was made using an exponential decay matching the predicted half correlation time.

6.4

Result for link margin estimation

The generation of random time series of 5 seconds with a 4 kHz sampling, i.e. 20.000 occurrences to ensure good statistical convergence, is very memory intensive. In the following, we limited our simulation to 5.000 pairs of profiles randomly selected in our database.

Generation of random time series was run two times, once using the half-correlation time computed using only r₀, , and and once with the moments of the temporal PSD added into the prediction of the half correlation time. Following results of power margin are given for an interleaver of 100ms and the MODCOD 9/10 with a required mutual information of 0.94 and MODCOD 3/10 with a required mutual information of 0.12.

On figure 13 and 14 we compared the probability density function of the “true” power margin (after SAOST) and the estimated power margin using only r₀, , and for MODCOD 3/10 and 9/10. When computing the error made on the power margin we see that it never exceeds 0.7 dB for MODCOD 3/10 and 0.6 dB for MODCOD 9/10 (in blue figure 13 and 14).

Figure 13:

PDF true power margin vs predicted using only r₀, and for MODCOD 3/10-interleaver 100ms.

Figure 14:

PDF true power margin vs predicted using only r₀, and for MODCOD 9/10-interleaver 100ms.

When the temporal PSD of scintillation is added to the input (yellow on figure 15 and 16), no significant improvement is observed, we even have a maximal error a tiny bit bigger with 0.8dB in both MODCOD case that might be attributed to the convergence error of the estimation due to the limited length of the considered time series.

Figure 15:

Distribution of power margin error [dB] for MODCOD 3/10 - interleaver 100ms. In blue using the model with r₀, and in yellow when the PSD of scintillation is added.

Figure 16:

Distribution and cumulative distribution of power margin error [dB] for MODCOD 9/10 - interleaver 100ms. In blue using the model with r₀, and in yellow when the PSD of scintillation is added.

In this specific case of 100ms interleaver, significantly bigger than the half correlation time (ranging in the 10^th of milliseconds), better prediction of the half correlation does not translate into better prediction of the power margin which confirms the relevance of the interleaving process.

6.5

Model summary and sensitivity analysis

We proposed a model that takes as input 4 moments (r₀, , and ) and gives as output a power margin estimation with an error bellow 1dB. This model includes intermediate steps of estimating PDF and autocorrelation time as well as random time series generation.

An interesting thing to look at is which one of the inputs has the biggest impact on the power margin estimation. Knowledge on the impact of each of the moments would enable us to improve the precision on the measurement of a given moment or on the contrary to reduce the required resolution. To do so we use classical tools of sensitivity analysis that are first and total order Sobol indices as well as Shapley values.

The Sobol indices³¹ are also known as variance-based sensitivity analysis. The variance of the output of our model is decomposed into fractions, which can be attributed to each of the moments used as inputs. The first order Sobol indices characterize the contribution of a given input to the output variance and the total Sobol indices measure the contribution to the output variance of the studied input, including all variance caused by its interactions, with other input variables.

Sobol indices are good to determine the sensitivity of the model caused by a specific input but works under the assumption of independence between inputs, their interpretation becomes hazardous in the case of correlated inputs. In our case, r₀, and are strongly correlated being all moments of the same profiles. To deal with correlated inputs, methods have been developed around Shapley indices^32,33.

On figure 17 we can see that r₀ has the biggest impact on the predicted value of the power margin, accounting for more than half of the variance according to Shapley values. This result was expected, as in the downlink case, r₀ intervenes in every contributors to the adaptive optics error budget considered here. We expect that , that has little impact here, will have a bigger impact on the uplink case as the anisoplanatic error (related to ) will become an important term in the error budget.

Figure 17

Sensitivity analysis to describe the impact of r₀, and on the prediction of the power margin.