A Computational Model of Cn2 Profile Inversion for Atmospheric Laser Communication in the Vertical Path

In this paper, an atmospheric structure constant Cn2 model is proposed for evaluating the channel turbulence degree of atmospheric laser communication. First, we derive a mathematical model for the correlation between the atmospheric coherence length r0, the isoplanatic angle θ0 and Cn2 using the Hufnagel–Valley (HV) turbulence model. Then, we calculate the seven parameters of the HV model with the actual measured r0 and θ0 data as input quantities, so as to draw the Cn2 profile and the θ0 profile. The experimental results show that the fitted average Cn2 contours and single-day Cn2 contours have superior fitting performance compared with our historical data, and the daily correlation coefficient between the single-day computed θ0 contours and the measured θ0 contours is up to 87%. This result verifies the feasibility of the proposed method. The results validate the feasibility of the proposed method and provide a new technical tool for the inversion of turbulence Cn2 profiles.


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 C 2 n profile through inverting the upper wind speed parameters of the turbulence parameters and the surface atmospheric refractive index structure constant using real-time measured r 0 and θ 0 data. From the above-mentioned research, it could be seen that when inverting the C 2 n 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 C 2 n 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 C 2 n profile. This research method simplifies the inversion method of the C 2 n 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 C 2 n in satellite-to-ground laser communication, in order to better evaluate communication error rate and design laser emission systems.

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. propagating through atmospheric turbulence. The atmospheric isoplana cates the angular correlation of the wavefront after the beacon light propa atmospheric turbulence, and the measurement diagram is shown in Figure  Both of them contain the path integral term of 2 n C . The generalized H the theoretical formulas of the whole atmosphere coherence length and the phere isoplanatic angle [15,16] are as follows: , k is the wave number, λ represents the wavelength an transmission path. Both Equations (1) and (2) contain the path integral term can be represented using the generalized HV [17,18] model as follows: ( )  Both of them contain the path integral term of C 2 n . 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 k = 2π/λ, k is the wave number, λ represents the wavelength and z denotes transmission path. Both Equations (1) and (2) contain the path integral term of C 2 n , which can be represented using the generalized HV [17,18] model as follows: where h = z cos Θ indicates the height and Θ is zenith angle. a 1 , b 1 and c jointly characterize the variation of turbulence intensity in the region at and above the top of the troposphere; a 2 and b 2 together characterize the variation of turbulence intensity in the tropospheric range; a 3 and b 3 are combined to characterize the turbulence intensity change in the boundary layer. a 2 indicates the intensity of turbulence at the beginning of the troposphere, while a 3 represents the variation of near-surface turbulence. b 1 , b 2 and b 3 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 Γ(z) function is the Gamma function [19]. Through considering Equations (4) and (5) in the case of vertical channels, we have Sensors 2023, 23, 5874 4 of 13 From the above-mentioned Equation (6), it can be known that there is a certain relationship between the coherence length r 0 and the isoplanatic angle θ 0 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 (a 1 , c, b 1 , a 2 , b 2 , a 3 , b 3 ) are determined, θ 0 can be obtained via solving r 0 . With such consideration, we propose a method to solve for the seven parameters, which can be expressed as where Represents the average value of the measured data r 0 . When R = r 0 , the variance value of the measured r 0 data is the smallest. R i is the i-th R value obtained via simulation calculation. R and R i should be as close as possible, and the degree of closeness is controlled by the accuracy G. θ i is the calculated representative value of the i-th θ 0 data at the same period as the measured r 0 data, and the relationship between this value and R i satisfies Equation (6), which can be used as a boundary condition. The scale factor M is determined according to the ratio of the average value of the measured area r 0 and θ 0 data, i.e., M = θ 0 /r 0 . 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 M within a reasonable range for filtering the calculated results. Therefore, the upper and lower limits M max and M min of the scale factor M should be selected according to the actual situation and fluctuate around the scale factor M. Once the value range of the scale factor M and the accuracy G are determined, and the range of seven parameters (a 1 , c, b 1 , a 2 , b 2 , a 3 , b 3 ) 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.

Experimental Analysis and Discussion
In order to verify the feasibility of the above theoretical method, in December 20 whole-layer atmospheric coherence length differential image motion monitor (DI and isoplanatic angle meter were tested in the Nanshan area of Xinjiang. Both the D and the isoplanatic angle measuring instrument use the stars in the air as beacons and the differential image motion method and the starlight scintillation method to mea respectively. That is, the plane wave emitted by a star propagates through the turb atmosphere and its wavefront distorts, and the wavefront distortion changes the pr gation direction and energy of the light wave.

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 r 0 and θ 0 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 r 0 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, D = 100 mm represents the pupil diameter of DIMM and d = 200 mm denotes the distance between the two pupils. σ 2 l and σ 2 t indicate the vertical and horizontal position variance, respectively. and light intensity of the star image change with the influence of a and the values of   When measuring θ 0 with an isoplanatic angle meter, the key technology is to fit the weighting function W(z) = Cz 5/3 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. ( ) where λ = 500 nm is the detection wavelength, D = 100 mm represents the pupil diameter of DIMM and d = 200 mm denotes the distance between the two pupils.   In (a,b), the three rings connected by solid blue dots are opaque (the black ring corresponds to the bright silver ring) and the three rings connected by the orange hollow circle are transparent parts (the white ring corresponds to the transparent ring). In (b), the inner and outer radii of the transparent rings, from inside to outside, are as follows: for the innermost first bright ring, the inner ring radius is 37.389 mm and the outer ring radius is 43.840 mm; for the middle second bright ring, the inner ring radius is 62.890 mm and the outer ring radius is 69.240 mm; for the outermost third bright ring, the inner ring radius is 81.940 mm and the outer ring radius is 101.600 mm; finally, the outermost black ring is the reserved installation allowance, and the size can be set according to the actual situation.
The weighting function of the three-ring apodizing mirror is shown below [22][23][24]: where ρ is the radius, 0 J represents the zero-order Bessel function and z denotes  In (a,b), the three rings connected by solid blue dots are opaque (the black ring corresponds to the bright silver ring) and the three rings connected by the orange hollow circle are transparent parts (the white ring corresponds to the transparent ring). In (b), the inner and outer radii of the transparent rings, from inside to outside, are as follows: for the innermost first bright ring, the inner ring radius is 37.389 mm and the outer ring radius is 43.840 mm; for the middle second bright ring, the inner ring radius is 62.890 mm and the outer ring radius is 69.240 mm; for the outermost third bright ring, the inner ring radius is 81.940 mm and the outer ring radius is 101.600 mm; finally, the outermost black ring is the reserved installation allowance, and the size can be set according to the actual situation. The weighting function of the three-ring apodizing mirror is shown below [22][23][24]: where ρ is the radius, J 0 represents the zero-order Bessel function and z denotes transmission path. k is the wave number, k = 2π/λ, λ is the wavelength, κ represents the space wave number, κ max = 2π/l 0 and κ min = 2π/L 0 . l 0 and L 0 indicate the inner and outer scales of atmospheric turbulence, respectively. P(κρ) is the transmittance function, as shown in Equation (13): where R 1 , R 2 , R 3 , R 4 , R 5 and R 6 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, l 0 = 0.005 m and L 0 = 10 m, the result of fitting calculation C is C = 8.847 × 10 −17 m 4 . Combining the weighting function W(z) = Cz 5/3 obtained through fitting with the normalized variance σ 2 s (0) of the light intensity fluctuation, we have where A represents the light transmission area of the three-ring apodizing mirror, A = 0.0156 m 2 , and C denotes the fitting coefficient of the three-ring apodizing mirror. Considering Equations (14) and (3), we can obtain θ 0 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 σ 2 s (0), the value of θ 0 can be obtained, as shown in Equation (15). The optical receiving system is shown in Figure 5.
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 ( ) 2 s 0 σ , the value of 0 θ can be obtained, as shown in Equation (15). The optical receiving system is shown in Figure 5.

Model Analysis and Experimental Discussion
The average C 2 n profile and single-day C 2 n 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 r 0 and θ 0 .
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 C 2 n profile in the Xianghe model. As the latitudes of Nanshan area and Xianghe area in Xinjiang are relatively close, the value ranges of a 1 , c and b 1 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 a 2 is expanded to [1 × 10 −18 , 1 × 10 −15 ], the range of b 2 is also expanded to [1500, 3000] and the range of a 3 is changed to [1 × 10 −17 , 1 × 10 −14 ] while the range of b 3 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 C 2 n profile and the single-day C 2 n profile, the expression of the C 2 n profile is as follows.
To further analyze the turbulence change, we plotted the daily C 2 n profile of the Nanshan area, the average C 2 n profile and the turbulence profile of the Beijing Xianghe Model. The expression of the Xianghe model is  (20) 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 a 3 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.

1000
2300 52 Observing Figure 6, we can know that the overall trend of the turbulence s model in the Nanshan area is similar to that in the Xianghe area. At heights of 10 km, the "trough" and "peak" of the turbulence intensity with height are which is consistent with the variation pattern of the HV turbulence model. In ge average turbulence profile in the Nanshan area is more in line with the Xiang because the data is averaged after multiple measurements, and the simulated t profile from the averaged data is more consistent with the long-term turbulence variation in the Nanshan region. The Nanshan area is closer to the Xianghe area sion, but Nanshan is at a high altitude (around 2000 m), resulting in both simila differences in the details of turbulence intensity variation between the two areas seen from Figure 6 that the intensity of near-surface turbulence in the Nansh greater than that in the Xianghe area, and the variation of near-surface turbule sity mainly depends on the 3 a parameter, which is greatly influenced by the ne wind speed, ground temperature, humidity and other climatic conditions. The intensity of near-surface turbulence varies with changing climatic conditions in regions.  Moreover, the single-day C 2 n 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 a 2 when other parameters are constant. The decrease in a 2 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 a 1 when other parameters are unchanged, and an increase in the value of a 1 causes the profile to shift to the right above 10 km. However, the trend is not highlighted in the figure because the parameters a 1 , c and b 1 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 θ 0 was derived from the measured r 0 data. The calculated θ 0 values on a single day were compared with the actual measured θ 0 values, as shown in Figure 7. Moreover, the single-day 2 n C 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 2 a when other parameters are constant. The decrease in 2 a 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 1 a when other parameters are unchanged, and an increase in the value of 1 a causes the profile to shift to the right above 10 km. However, the trend is not highlighted in the figure because the parameters 1 a , c and 1 b 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)   As shown in Figure 7, the trend of the theoretically calculated values of the wholeatmosphere 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 R xy was obtained from the following Equation (21). The variation trend is shown in Figure 8.
where x i and y i represent the actual measurements and calculated value of θ 0 , respectively. n represents the quantity value of the data. X and Y denote the average of actual measurements and calculated θ 0 value, respectively. As shown in Figure 7, the trend of the theoretically calculated values of the wholeatmosphere 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 xy R was obtained from the following Equation (21). The variation trend is shown in Figure 8.
where i x and i y represent the actual measurements and calculated value of 0 θ , respectively. n represents the quantity value of the data. X and Y denote the average of actual measurements and calculated 0 θ 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. 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.

Conclusions
In this paper, based on the generalized HV model, a theoretical relationship equation between r 0 and θ 0 is derived, which establishes a certain connection between them numerically and provides a reference for related studies involving r 0 and θ 0 . First, a new method is proposed to solve the seven parameters of the generalized HV model using the whole-atmosphere coherence length and isoplanatic angle to invert the C 2 n profile. The average C 2 n profile and single-day C 2 n profile of the Nanshan area are obtained using the proposed method's inversion with the measured r 0 and θ 0 data as inputs, and the trend is in good agreement with that of the Xianghe model, which is in accordance with the turbulence variation law. Moreover, there is a high correlation between the calculated daily variation profile of the whole-atmosphere isoplanatic angle and the measured θ 0 profile, and the correlation coefficient's average value of 16 sets of data has reached 87%. The analytical results of the inverse C 2 n profile and the calculated θ 0 profile better support the feasibility and correctness of the proposed inversion method, which could provide a new reference for the better study of the C 2 n profile inversion method.