1. Introduction
Free space optical communication (FSOC) is one of the main communication technologies for future 6G, with high speed, no electromagnetic interference, high bandwidth, etc. However, turbulence effects in the atmosphere cause wavefront distortion of optical waves, which increases the communication BER. Since the atmosphere is changing in real time, how to better detect changes in atmospheric turbulence is the focus of our research and the basis for our adaptive correction of communication using the nature of atmospheric turbulence [
1,
2].
The current measurement methods of atmospheric turbulence intensity are mainly divided into two categories: direct measurement and indirect measurement. The most common method of direct measurement is to release a sounding balloon at the measurement site and measure the atmospheric turbulence profile through the sensors mounted on it. However, due to the influence of wind and other factors, its movement direction is uncontrollable and the real-time performance is poor. In addition, the use of optical methods to directly measure turbulence profiles mainly includes radar [
3], SCIDAR [
4] (SCIntillation Detection and Ranging), MASS [
5] (Multi-Aperture Scintillation Sensor), etc. Lucie Rottner et al. [
6] proposed a wind reconstruction method applied to a five-beam wind Doppler lidar (i.e., Leosphere’s WindCube model) that relies on particle filtering to associate lidar data with numerical particles to obtain turbulence estimates available for each new observation. Indirect measurement is made mainly through measuring meteorological parameters and turbulence parameters such as atmospheric coherence length, isoplanatic angle, etc., to invert, predict [
7,
8] and estimate the atmospheric turbulence profile. Rafalimanana A et al. [
9] proposed a forecasting study based on the Weather Research and Forecasting (WRF) model. It can predict and describe a set of useful meteorological parameters related to atmospheric physics (pressure, temperature, relative humidity, wind speed, direction, etc.), and then inject the predicted parameters into the optical turbulence model to calculate the refractive index structure constant. Santasri R. Bose-Pillai et al. [
10,
11] proposed a method for estimating turbulent parameters through deriving a weighting function that relates the turbulence intensity along the path to the shifts measured and then estimates the turbulent parameters. Wang Yao [
12] et al. took five conventional meteorological parameters as input and used an artificial neural network to predict the profile of the sea surface near Mauna Loa, Hawaii, for one month. Zhang et al. [
13] based their study on the artificial neural network algorithm and established an artificial neural network model based on the data to predict the upper atmospheric turbulence profile. The predicted value simulated using the neural network algorithm is in good agreement with the actual turbulence profile in the Maoming area, which proves the feasibility and reliability of using a neural network to simulate the atmospheric turbulence profile. Based on the Hufnagel–Valley (HV) model, Robert K. Tyson et al. [
14] obtained the
profile through inverting the upper wind speed parameters of the turbulence parameters and the surface atmospheric refractive index structure constant using real-time measured
and
data. From the above-mentioned research, it could be seen that when inverting the
profile, the input includes not only meteorological parameters such as pressure and temperature, but also turbulence parameters such as atmospheric coherence length and isoplanatic angle. There are many types of initial input parameters, and the overall amount of data is huge, resulting in a large amount of calculation in the inversion process, which is complicated and cannot obtain useful information for the inversion
profile from a single type of parameter.
In this paper, we propose a new method to invert the atmospheric turbulence profile based on the generalized HV mode, taking the atmospheric coherence length and atmospheric isoplanatic angle as inputs. Based on the generalized Hufnagel–Valley turbulence model, the method deduces the theoretical relationship between the atmospheric coherence length and the isoplanatic angle, solves the seven parameters of the generalized Hufnagel–Valley turbulence model through the inversion algorithm, and then obtains the profile. This research method simplifies the inversion method of the profile. It not only provides a new idea for the inversion of the turbulence profile, but also develops a method to determine the parameters of the HV model turbulence profile mode using the coherence length and isoplanatic angle of the whole atmosphere layer, which can ensure high accuracy and require fewer input data. This work can provide a theoretical reference for evaluating the profile performance of atmospheric turbulence structure parameter in satellite-to-ground laser communication, in order to better evaluate communication error rate and design laser emission systems.
2. Correlation Model Establishment and Inversion
Both the atmospheric isoplanatic angle and the atmospheric coherence length can represent the variation of the intensity of atmospheric turbulence in the transmission path, and the atmospheric coherence length represents the diffraction limit of light waves propagating through atmospheric turbulence. The atmospheric isoplanatic angle indicates the angular correlation of the wavefront after the beacon light propagates through atmospheric turbulence, and the measurement diagram is shown in
Figure 1.
Both of them contain the path integral term of
. The generalized HV model and the theoretical formulas of the whole atmosphere coherence length and the whole atmosphere isoplanatic angle [
15,
16] are as follows:
where
,
is the wave number,
represents the wavelength and
denotes transmission path. Both Equations (1) and (2) contain the path integral term of
, which can be represented using the generalized HV [
17,
18] model as follows:
where
indicates the height and
is zenith angle.
,
and
jointly characterize the variation of turbulence intensity in the region at and above the top of the troposphere;
and
together characterize the variation of turbulence intensity in the tropospheric range;
and
are combined to characterize the turbulence intensity change in the boundary layer.
indicates the intensity of turbulence at the beginning of the troposphere, while
represents the variation of near-surface turbulence.
,
and
represent the attenuation speed of each turbulent layer with the increase of height.
Substituting Equation (3) into Equations (1) and (2), we can get
where the
function is the Gamma function [
19]. Through considering Equations (4) and (5) in the case of vertical channels, we have
From the above-mentioned Equation (6), it can be known that there is a certain relationship between the coherence length
and the isoplanatic angle
of the whole atmosphere layer. Therefore, Equation (6) can be written as
From Equations (7) and (8), we can know that as long as the values of the seven parameters (
,
,
,
,
,
,
) are determined,
can be obtained via solving
. With such consideration, we propose a method to solve for the seven parameters, which can be expressed as
where
,
,
R Represents the average value of the measured data
. When
, the variance value of the measured
data is the smallest.
is the
i-th
R value obtained via simulation calculation.
R and
should be as close as possible, and the degree of closeness is controlled by the accuracy
.
is the calculated representative value of the
i-th
data at the same period as the measured
data, and the relationship between this value and
satisfies Equation (6), which can be used as a boundary condition. The scale factor
is determined according to the ratio of the average value of the measured area
and
data, i.e.,
. However, as for the simulation calculation, sometimes the calculated value and the result are close but not equal, which will lead to errors in the calculation results. The problem can be avoided through controlling the scaling factor
within a reasonable range for filtering the calculated results. Therefore, the upper and lower limits
and
of the scale factor
should be selected according to the actual situation and fluctuate around the scale factor
. Once the value range of the scale factor
M and the accuracy
G are determined, and the range of seven parameters (
,
,
,
,
,
,
) is input, the seven parameter values can be simulated using Equations (9) and (10), which greatly reduces the number of initial data required for inversion of turbulence profile parameters. Therefore, large amounts of meteorological data input are no longer necessary, avoiding the tediousness of data collection and processing, and having wider applicability in practical engineering applications. The entire process is shown in
Figure 2.
3. Experimental Analysis and Discussion
In order to verify the feasibility of the above theoretical method, in December 2020, a whole-layer atmospheric coherence length differential image motion monitor (DIMM) and isoplanatic angle meter were tested in the Nanshan area of Xinjiang. Both the DIMM and the isoplanatic angle measuring instrument use the stars in the air as beacons and use the differential image motion method and the starlight scintillation method to measure, respectively. That is, the plane wave emitted by a star propagates through the turbulent atmosphere and its wavefront distorts, and the wavefront distortion changes the propagation direction and energy of the light wave. On the imaging target surface, the position and light intensity of the star image change with the influence of atmospheric turbulence, and the values of
and
are obtained through measuring the statistics of the change of the position and light intensity. The specific measurement method is shown in
Figure 3:
The value of
is obtained through calculating the horizontal and vertical position variance of the stellar image imaged on the target surface of the CCD camera and substituting it into Equation (11), utilizing DIMM [
20,
21].
where
= 500 nm is the detection wavelength,
= 100 mm represents the pupil diameter of DIMM and
= 200 mm denotes the distance between the two pupils.
and
indicate the vertical and horizontal position variance, respectively.
When measuring
with an isoplanatic angle meter, the key technology is to fit the weighting function
using a three-ring apodizing mirror. The aperture of the three-ring apodizing mirror is circularly symmetrical, and its physical diagram and structure diagram are shown in
Figure 4.
The weighting function of the three-ring apodizing mirror is shown below [
22,
23,
24]:
where
is the radius,
represents the zero-order Bessel function and
denotes transmission path.
is the wave number,
,
is the wavelength,
represents the space wave number,
and
.
and
indicate the inner and outer scales of atmospheric turbulence, respectively.
is the transmittance function, as shown in Equation (13):
where
,
,
,
,
and
are the ring radii of the three-ring apodizing mirror from inside to outside. When the zenith angle is set to 0°, with
= 500 nm,
= 0.005 m and
= 10 m, the result of fitting calculation
is
. Combining the weighting function
obtained through fitting with the normalized variance
of the light intensity fluctuation, we have
where
represents the light transmission area of the three-ring apodizing mirror,
, and
denotes the fitting coefficient of the three-ring apodizing mirror. Considering Equations (14) and (3), we can obtain
as follows:
Equation (15) indicates that the solution of the isoplanatic angle has nothing to do with the wavelength. When the stellar light wave is transmitted through the turbulent atmosphere, the distorted wavefront is modulated by the three-ring apodizing mirror, received by the optical receiving system and finally converged on the target surface of the charge coupled device (CCD) camera to form a star point image. Through measuring the light intensity of the star point image and calculating its normalized light intensity fluctuation variance
, the value of
can be obtained, as shown in Equation (15). The optical receiving system is shown in
Figure 5.
4. Model Analysis and Experimental Discussion
The average profile and single-day profile (e.g., data 1: 11 December 2020, data 2: 13 December 2020, data 3: 16 December 2020) in the Nanshan area during the measurement period were calculated based on the proposed inversion profile method, after sorting out the measurement data of and .
From the simulation, the selected parameter ranges are shown in
Table 1. The selection of the seven parameter ranges comprehensively considered the changes in the profile in Xinjiang [
25] and the
profile in the Xianghe model. As the latitudes of Nanshan area and Xianghe area in Xinjiang are relatively close, the value ranges of
,
and
remain consistent. According to Ref. [
25], the near-surface turbulence intensity in the Altay and Korla regions of Xinjiang is about 10
−16 m
−2/3, and it declines rapidly with height. The turbulence intensity in the range of 5–30 km changes within [1 × 10
−18, 1 × 10
−16], and the degree of turbulence intensity decline cannot be accurately estimated. In order to make the simulation calculation range closer to the real situation, the range of
is expanded to [1 × 10
−18, 1 × 10
−15], the range of
is also expanded to [1500, 3000] and the range of
is changed to [1 × 10
−17, 1 × 10
−14] while the range of
is reduced to [200, 800]. Based on the above-mentioned method,
R and
M are determined according to the measured atmospheric coherence length and isoplanatic angle on the day of measurement. The specific parameter ranges are shown in
Table 1.
Obtaining the seven parameter values of the average
profile and the single-day
profile, the expression of the
profile is as follows.
To further analyze the turbulence change, we plotted the daily
profile of the Nanshan area, the average
profile and the turbulence profile of the Beijing Xianghe Model. The expression of the Xianghe model is
Observing
Figure 6, we can know that the overall trend of the turbulence simulation model in the Nanshan area is similar to that in the Xianghe area. At heights of 5 km and 10 km, the “trough” and “peak” of the turbulence intensity with height are reflected, which is consistent with the variation pattern of the HV turbulence model. In general, the average turbulence profile in the Nanshan area is more in line with the Xianghe model, because the data is averaged after multiple measurements, and the simulated turbulence profile from the averaged data is more consistent with the long-term turbulence intensity variation in the Nanshan region. The Nanshan area is closer to the Xianghe area in dimension, but Nanshan is at a high altitude (around 2000 m), resulting in both similarities and differences in the details of turbulence intensity variation between the two areas. It can be seen from
Figure 6 that the intensity of near-surface turbulence in the Nanshan area is greater than that in the Xianghe area, and the variation of near-surface turbulence intensity mainly depends on the
parameter, which is greatly influenced by the near-surface wind speed, ground temperature, humidity and other climatic conditions. Therefore, the intensity of near-surface turbulence varies with changing climatic conditions in different regions.
Moreover, the single-day
profiles in
Figure 6a–c can well reflect the turbulence variation in the Nanshan area. As shown in the figure, the overall trend of the single-day turbulence profile is in accordance with the law of turbulence change, the variability is mainly reflected at 5 km and 10 km and there is no significant change in the altitude region above 15 km, which indicates that the intensity of atmospheric turbulence changes drastically in the range of 15 km.
It is worth noting that the cause of the “trough” at 5 km is mainly caused by parameter when other parameters are constant. The decrease in makes the “trough of the wave” sink even further, i.e., the intensity of turbulence at the beginning of the troposphere changes. Similarly, the “crest” at 10 km is the result of the action of when other parameters are unchanged, and an increase in the value of causes the profile to shift to the right above 10 km. However, the trend is not highlighted in the figure because the parameters , and have different degrees of variation and their combined effect causes this phenomenon.
To further explore the rationality of this theoretical formula, the single daily profile parameter values in
Table 1 were substituted into Equation (10), and then
was derived from the measured
data. The calculated
values on a single day were compared with the actual measured
values, as shown in
Figure 7.
As shown in
Figure 7, the trend of the theoretically calculated values of the whole-atmosphere isoplanatic angle throughout the day is essentially consistent with the actual measured values, showing alternating up and down in the local time range, (i.e., the theoretical value of the isoplanatic angle in
Figure 7a is slightly larger than the measured value in the period of 11:15∼13:30, and the theoretical values of the isoplanatic angle are smaller than the measured values in the period of 15:00–18:00 in
Figure 7b,c).
The above-mentioned situation may be caused by the difference in the measurement principle between DIMM and the isoplanatic angle meter. DIMM inverts the atmospheric coherence length value based on the position variance caused by the jitter of the measured star image, and the isoplanatic angle meter inverts the isoplanatic angle value on the basis of the variance of the measured star image’s light intensity fluctuation, which causes the difference in the details of the result. However, on the other hand, the consistency of the overall trend also reflects the accuracy of the overall measurement of the two approaches.
The theoretical data of the whole-atmosphere isoplanatic angle were calculated based on 16 sets of measurements, and the correlation coefficient between the measured isoplanatic angle
was obtained from the following Equation (21). The variation trend is shown in
Figure 8.
where
and
represent the actual measurements and calculated value of
, respectively.
represents the quantity value of the data.
and
denote the average of actual measurements and calculated
value, respectively.
Observing
Figure 8, we can know that the correlation coefficients between the calculated and measured values of the atmospheric isoplanatic angle are above 80%, with an average value of 0.8195 and the maximum value reaching 0.8708. Consequently, the isoplanatic angle data obtained from the theoretical equation have a fine correlation between the overall trend and the measured values, which further proves the correctness of the theoretical equation and the inversion method.