Simulation and Analysis of Mie-Scattering Lidar-Measuring Atmospheric Turbulence Profile

Based on the residual turbulent scintillation theory, the Mie-scattering lidar can measure the intensity of atmospheric turbulence by detecting the light intensity scintillation index of the laser return signal. In order to evaluate and optimize the reliability of the Mie-scattering lidar system for detecting atmospheric turbulence, the appropriate parameters of the Mie-scattering lidar system are selected and optimized using the residual turbulent scintillation theory. Then, the Fourier transform method is employed to perform the numerical simulation of the phase screen of the laser light intensity transformation on the vertical transmission path of atmospheric turbulence. The phase screen simulation, low-frequency optimization, and scintillation index calculation methods are provided in detail, respectively. Based on the phase distribution of the laser beam, the scintillation index is obtained. Through the relationship between the scintillation index and the atmospheric turbulent refractive index structure constant, the atmospheric turbulence profile is inverted. The simulation results show that the atmospheric refractive index structure constant profile obtained by the iterative method is consistent with the input HV5/7 model below 6500 m, which has great guiding significance to carry out actual experiments to measure atmospheric turbulence using the Mie lidar.


Introduction
Atmospheric turbulence is a kind of irregular vortex motion and is generated mainly due to three reasons: wind shear caused by the drag of the air on the Earth's surface, thermal convection caused by surface heat radiation, and fluctuation in the temperature and velocity fields due to heat-releasing phase transition processes (deposition and crystallization) [1]. The physical properties of turbulence, such as pressure, speed, and temperature, are random. The atmospheric turbulence effect will cause a series of optical phenomena, including light intensity scintillation, beam drift, beam expansion, and arrival angle fluctuation [2]. Among them, the light intensity scintillation will cause the loss of received power, thereby reducing the signal-to-noise ratio of the system. In fact, atmospheric turbulence has a certain impact on the study of adaptive optics [3,4], atmospheric laser communication [5], astronomical site selection, and astronomical observation image analysis [6]. However, there is no complete model to fully describe the impact of atmospheric turbulence [7].
The Navier-Stokes (NS) equation is the basis for many atmospheric turbulence calculations. According to the different processing scales of atmospheric turbulence, it is divided into direct numerical simulation (DNS), Reynolds averaged Navier-Stokes (RANS), and large eddy simulation (LES). Among them, the direct numerical simulation method does not need any turbulence model, artificial assumptions, or empirical constants, and in principle can solve all atmospheric turbulence problems [8]. Therefore, numerical modeling is usually used as one of the main research methods to study atmospheric turbulence. Numerical modeling usually includes the Fourier transform method (FFT), Zernike polynomial higher-order method (Zernike), fractal method (Fractal), and some hybrid methods, among which the Fourier transform method has been applied widely in the turbulent phase screen simulation due to its fast calculation advantage, especially in the occasions with high real-time requirements, such as space optical communication and atmospheric laser transmission. However, the Fourier transform method has serious defects in the low-frequency part, which makes the phase screen have a large error in characterizing large-scale turbulent eddies [9]. Therefore, some research has been performed to overcome the defect. In 1992, Lane et al. proposed a multi-order harmonic method to improve the low-frequency compensation, simulation accuracy, and calculation speed [10]. In 2008 and 2009, Qian et al. studied the distribution of the phase screen and the selection of calculation parameters [11,12]. In 2013, Herman and Strungala proposed a method of adding sub-harmonics to perform low-frequency compensation on the phase screen generated by the FFT method [13]. In 2013, Charnotskii et al. proposed a sparse spectral model, which superimposed the dense sampling of the low-frequency part and the sparse sampling of the high-frequency part, and achieved good results [14]. In 2017, Feng et al. proposed a phase screen simulation of atmospheric turbulence based on wavelet analysis and found that it was in good agreement with the von Karman spectrum and the computational complexity was reduced [15]. In 2019 and 2020, Zhang et al. proposed an optimization algorithm for sparsely resampling the low-frequency region with the help of the gravitational search algorithm to form a low-frequency compensation screen [16,17].
Moreover, in using the scintillation index to study atmospheric turbulence, some research has been carried out. In 2006, Rao et al. conducted a simulation study on the scintillation index of Gaussian beams in atmospheric turbulence and verified that the scintillation index increases with the distance from the optical axis [18]. In 2017, Han et al. used the step Fourier method to simulate the two-way transmission of the laser and analyzed the backward enhancement effect of the plane mirror and the angle of reflection on the scintillation index by calculating the scintillation index of the retroreflection spot [19]. In 2018, Li et al. simulated the scintillation index of plane waves propagating in atmospheric turbulence with the non-Kolmogorov spectrum and compared the variation of the scintillation index of the Kolmogorov spectrum and the non-Kolmogorov spectrum alone on the turbulent path [20].
An active optical remote sensing detection method is important for atmospheric turbulence measurement. The differential image movement measurement method (DIM) uses artificial light sources, which makes the measurement susceptible to light sources and has low resolution [21]. The Doppler wind profiler (MST) measurement method has high sensitivity and versatility but is easily affected by weather and has a limited measurement range [22]. As an important active remote sensing method, lidar has the characteristics of high resolution and high precision and is a reliable and effective means to detect atmospheric parameters. The backscattering amplification effect of the lidar (BSA) measurement method has high spatial, temporal, and resolution features as well as good real-time performance, and is not affected by time integration [23], while the scintillation lidar (DCIM) has accurate detection and low measurement cost [24]. In this paper, considering the advantages and disadvantages of various measurement methods, a Mie-scattering lidar system is selected for atmospheric turbulence measurement. According to the residual turbulent scintillation (RTS) theory [25], it was improved by adding small aperture diaphragms with appropriate apertures [26]. The atmospheric turbulence profile is obtained by inversion using the light intensity information at different heights of the return signal of the Mie-scattering lidar system.
In this paper, based on the RTS theory, the Mie-scattering lidar is selected to measure the intensity of atmospheric turbulence and the numerical simulation is carried out by using the phase screen simulation. In the simulation, the improved fast Fourier transform method with sub-harmonics is used to simulate the transmission of laser beams in atmospheric turbulence. After atmospheric extinction and backscattering, telescope imaging spots at different heights from the Mie-scattering lidar return signal are obtained. By analyzing the imaging spot, the scintillation change of its light intensity can be obtained. Moreover, by understanding the influence of specific measurement parameters on atmospheric turbulence detection, the reliability of its detection performance can be judged and follow-up guidance for system evaluation and optimization can be provided.

Mie-Scattering Lidar
The schematic diagram of the Mie-scattering lidar is shown in Figure 1. The collimated laser beam emitted by the transmitting system is expanded by the beam expander and then transmitted to the atmosphere [27]. Through the combined action of atmospheric Fresnel diffraction, atmospheric extinction and atmospheric turbulence [28], the backscattered return signal is received by the receiving system and then passes through the spectroscopic system. The return signals of each wavelength are separated and the Mie-scattering signal with a wavelength of 532 nm is selected for our goal. Finally, the return signal is detected, amplified, and displayed by the signal acquisition system. Table 1 shows the system parameters of the Mie-scattering lidar [29].
In this paper, based on the RTS theory, the Mie-scattering lidar is selected to me the intensity of atmospheric turbulence and the numerical simulation is carried o using the phase screen simulation. In the simulation, the improved fast Fourier tran method with sub-harmonics is used to simulate the transmission of laser bea atmospheric turbulence. After atmospheric extinction and backscattering, tele imaging spots at different heights from the Mie-scattering lidar return signal are obt By analyzing the imaging spot, the scintillation change of its light intensity c obtained. Moreover, by understanding the influence of specific measurement param on atmospheric turbulence detection, the reliability of its detection performance c judged and follow-up guidance for system evaluation and optimization can be pro

Mie-Scattering Lidar
The schematic diagram of the Mie-scattering lidar is shown in Figure 1 collimated laser beam emitted by the transmitting system is expanded by the expander and then transmitted to the atmosphere [27]. Through the combined act atmospheric Fresnel diffraction, atmospheric extinction and atmospheric turbulenc the backscattered return signal is received by the receiving system and then through the spectroscopic system. The return signals of each wavelength are sep and the Mie-scattering signal with a wavelength of 532 nm is selected for our goal. F the return signal is detected, amplified, and displayed by the signal acquisition sy Table 1 shows the system parameters of the Mie-scattering lidar [29].   The improved Mie-scattering lidar system is mainly based on the RTS theory, and a suitable aperture diaphragm is added to the original Mie-scattering lidar system to meet the needs of atmospheric turbulence detection. In RTS theory [25], when light propagates in the atmosphere, if the radius of the scatterer is smaller than the spatial correlation scale l I (l I = √ λL, where λ is the wavelength of the incident laser light, L is the laser detection distance), for any large receiving aperture D, the backscattered light intensity signal can reduce the aperture averaging effect, and there will be some light intensity fluctuations.
The following four conditions are required to be satisfied: (1) the single scattering is approximate and the aerosol particle size is uniform; then, the laser detection distance L, the average radius of aerosol particles a, and the incident light wavelength λ, satisfy L a 2 /λ; (2) atmospheric optical thickness τ 1; (3) the radial light intensity correlation scale l I (l I ≈ L) is much larger than the scatterer scale cτ/2 (where c is the speed of light, τ is the pulse width), and the scatterer scale is much larger than the wavelength, that is, L cτ/2 λ; and (4) under weak fluctuations (σ I 2 0.6), the telescope receiving aperture D, the turbulent coherence length of plane waves ρ 0 , l I and ρ cv satisfy D ρ 0 l I ρ cv (where σ I 2 is the logarithmic light intensity fluctuation variance; where ρ cv is the light intensity coherence length derived from the Van-Cittert-Zernike theory, a ef is the effective size of the beam from the light source L, and k is the wave number, k = 2π/λ). The lidar can be distinguished by two characteristic scales, namely, aerosol spot scale, l a = λF/D, and turbulence-related scale, l t = Fl I /L, where F is the focal length of the telescope. Note that l a /l t = D/ √ λL, so it generally satisfied l t l a for most large-aperture receiving systems (D √ λL). Because the spatial scales l a and l t of the scintillation spot are completely different, three different expressions for the fluctuation of the backscattered light intensity can be obtained by changing the diameter of the telescope's field of view diaphragm d 0 : where, σ p 2 is the variance of the intensity fluctuation of the backscattered light passing through the aperture; σ I 2 is the variance of the intensity fluctuation of the probe beam. Equation (1) shows that the variance of light intensity fluctuations includes aerosol speckle fluctuations and turbulence fluctuations. Equation (2) shows that the diaphragm smooths the aerosol speckle fluctuations but includes turbulence fluctuations. Equation (3) shows that aerosol speckle fluctuations and turbulence fluctuations are smoothed by the diaphragm, and the total received energy does not fluctuate.
In this paper, the selected telescope aperture D is 10 inches (254 mm), the repetition frequency f is 10 Hz, and the adjustable focal length F is 2500 mm. Considering the experimental operability and the actual laboratory conditions, the selected field of view of the telescope is d 0 ≥ 0.2 mm. According to the RTS theory, the spatial correlation scale l I (l I = √ λL) of the light intensity fluctuation is calculated as shown in Figure 2. The spatial correlation scale l I increases with the increase in the detection distance and incident wavelength. When the detection distance is the maximum value of 3 km and the incident laser wavelength takes the maximum value of 1064 mm, l I gets the maximum value of 57.5 mm, which obviously satisfies D l I . In addition, according to our laboratory's long-term measurement results of aerosol particle spectral distribution in the Yinchuan area, China, it can be known that the aerosol particle size a = 0.9 µm. Yinchuan area is located in the junction of the four deserts, namely Badain Jaran, Ulan Buhe, Tengger, and Mu Us deserts, combined with the typical temperate continental semi-arid climate, the aerosol particle size in spring is slightly larger than the average; obviously, all meet the condition a r l I . Therefore, the improved Mie-scattering lidar system can satisfy the RTS theory. By changing the diameter of the field of view diaphragm, the information characterizing atmospheric turbulence can be obtained from the return signal.
ors 2022, 22, x FOR PEER REVIEW China, it can be known that the aerosol particle size = 0.9 μm. Yinchuan are in the junction of the four deserts, namely Badain Jaran, Ulan Buhe, Tengger, deserts, combined with the typical temperate continental semi-arid climate, particle size in spring is slightly larger than the average; obviously, all meet th ≪ . Therefore, the improved Mie-scattering lidar system can satisfy the R By changing the diameter of the field of view diaphragm, the information cha atmospheric turbulence can be obtained from the return signal. Gaussian beam is one of the most commonly used models in lasers. In fa intensity of a Gaussian beam in a plane is perpendicular to the direction of p A Gaussian beam in a homogeneous medium is expressed as: where, ( ) is the beam radius (i.e., ( ) = 1 + ⁄ ), is the minim of ω(z) at z = 0, which is called the light waist radius, r is the distance between space (x, y, z) and the light source; and k is the wave number (i.e., = 2 ⁄ shows the light intensity distribution of the Gaussian beam with a beam wai 20 mm, wavelength = 532 nm, and = 0. It can be seen that the light spot circle, and the light spot is at the center. Gaussian beam is one of the most commonly used models in lasers. In fact, the light intensity of a Gaussian beam in a plane is perpendicular to the direction of propagation. A Gaussian beam in a homogeneous medium is expressed as: where, ω(z) is the beam radius (i.e., ω(z) = ω 0 1 + z 2 /z 2 0 ), ω 0 is the minimum value of ω(z) at z = 0, which is called the light waist radius, r is the distance between a point in space (x, y, z) and the light source; and k is the wave number (i.e., k = 2π/λ). Figure 3 shows the light intensity distribution of the Gaussian beam with a beam waist radius of 20 mm, wavelength λ = 532 nm, and z = 0. It can be seen that the light spot is a regular circle, and the light spot is at the center.

Numerical Simulation of Gaussian Beam Transmission in Atmosphere
When the laser is transmitted in the atmosphere, it is subjected to the dual action of the turbid atmosphere and continuous turbulence and then diffracted in the atmospheric space, resulting in atmospheric scattering and absorption and a series of turbulent effects caused by turbulent refraction. In the phase screen model, the light propagation in atmospheric turbulence mainly starts from the light propagation equation. Assuming that the transmission path of the light beam is composed of several thin phase screens distributed in a vacuum environment, and when the phase change caused by the refractive index of the medium is particularly small, the process of light propagation in a vacuum environment and the phase modulation of the medium can be regarded as completely independent of each other and completed simultaneously [1]. That is, the multi-layer phase screen model of light propagation in a continuous random medium is shown in Figure 4, in which the circles or ellipses represent random media in the atmosphere, such as atmospheric molecules, aerosols, etc. The phase disturbance is firstly generated by a thin phase screen with negligible thickness. After propagating through the vacuum of distance Δ , the phase disturbance of the previous phase screen is superimposed on this phase screen to repeat the phase modulation, with the above process repeating until the end of the transmission.

Numerical Simulation of Gaussian Beam Transmission in Atmosphere
When the laser is transmitted in the atmosphere, it is subjected to the dual action of the turbid atmosphere and continuous turbulence and then diffracted in the atmospheric space, resulting in atmospheric scattering and absorption and a series of turbulent effects caused by turbulent refraction. In the phase screen model, the light propagation in atmospheric turbulence mainly starts from the light propagation equation. Assuming that the transmission path of the light beam is composed of several thin phase screens distributed in a vacuum environment, and when the phase change caused by the refractive index of the medium is particularly small, the process of light propagation in a vacuum environment and the phase modulation of the medium can be regarded as completely independent of each other and completed simultaneously [1]. That is, the multi-layer phase screen model of light propagation in a continuous random medium is shown in Figure 4, in which the circles or ellipses represent random media in the atmosphere, such as atmospheric molecules, aerosols, etc. The phase disturbance is firstly generated by a thin phase screen with negligible thickness. After propagating through the vacuum of distance ∆z, the phase disturbance of the previous phase screen is superimposed on this phase screen to repeat the phase modulation, with the above process repeating until the end of the transmission.
Assuming that the light travels along the z direction, the light field can be expressed as E = ue ikz . When the laser passes through the atmospheric turbulence with a refractive index of n = 1 + n 1 , under the condition of paraxial approximation, the light propagation satisfies the parabolic equation: where, k = 2π/λ, λ is the wavelength of the light wave, k is the wave number, and n 1 is the refractive index fluctuation. Let the transmission distance L of the light wave be divided into N z segments according to ∆z equal intervals, so the head and tail coordinates of the i-th segment are z i−1 and z i , respectively, and the i-th phase screen is set at the endpoint z i of each segment. For each phase screen, the grid spacing of ∆x is equally divided into N × N grids, and the above formula is solved by the step-by-step method. It can be known that the field distribution of the I + 1-th phase screen is: where, F −1 and F are the inverse Fourier transform and Fourier transform, u i is the field distribution of the i-th phase screen, ∆z i = z i+1 − z i is the distance between the phase screen, and κ x and κ y are the phase space wave number with the unit of rad/m, which are related to the scale of turbulence. Assuming that the light travels along the direction, the light field can be expressed as = . When the laser passes through the atmospheric turbulence with a refractiv index of = 1 + , under the condition of paraxial approximation, the light propagation satisfies the parabolic equation: is the wavelength of the light wave, is the wave number, and is the refractive index fluctuation.
Let the transmission distance L of the light wave be divided into segments ac cording to Δ equal intervals, so the head and tail coordinates of the i-th segment ar and , respectively, and the i-th phase screen is set at the endpoint of each seg ment. For each phase screen, the grid spacing of ∆ is equally divided into × grids and the above formula is solved by the step-by-step method. It can be known that the field distribution of the I + 1-th phase screen is: where, ℱ and ℱ are the inverse Fourier transform and Fourier transform, is th field distribution of the i-th phase screen, ∆ = − is the distance between th phase screen, and and are the phase space wave number with the unit of / which are related to the scale of turbulence.
The most important problem in the numerical simulation method of light propaga tion is to construct the appropriate phase screen, so that it can accurately reflect the chang in the atmospheric turbulent refractive index as much as possible [30], that is, to choos the appropriate refractive index fluctuation . We use the spatial spectrum model of at mospheric turbulence to obtain a random field of the phase space, and then use Fourie transform to get the spatial distribution of the two-dimensional phase. The specific pro cess is as follows: generate a complex Gaussian random number matrix , select an at Δ Phase screens The most important problem in the numerical simulation method of light propagation is to construct the appropriate phase screen, so that it can accurately reflect the change in the atmospheric turbulent refractive index as much as possible [30], that is, to choose the appropriate refractive index fluctuation n 1 . We use the spatial spectrum model of atmospheric turbulence to obtain a random field of the phase space, and then use Fourier transform to get the spatial distribution of the two-dimensional phase. The specific process is as follows: generate a complex Gaussian random number matrix a R , select an atmospheric turbulence power density spectral function φ n to filter it to obtain a random field in the complex space, and perform inverse Fourier transform on the discrete complex random field to obtain the real phase distribution in space.

Direction of propagation
Atmospheric turbulence refractive index spectral function φ n and atmospheric turbulence phase spectral function φ φ using Kolmogorov spectrum [2] are written by: where, C n 2 is the atmospheric refractive index structural constant when the laser travels along the horizontal path, and the unit is m −2/3 ; κ is the spatial frequency; L 0 is the turbulent outer scale, and l 0 is the turbulent flow inner scale, which determines the largest and small eddies in the turbulent eddy.
The complex space random field obtained after filtering the complex Gaussian random number matrix a R is given by: where, A R and B R are Gaussian random numbers whose real and imaginary parts have a mean of 0 and a variance of 1. The real and imaginary parts of the phase field obtained by Equation (10) are independent of each other, and both satisfy the spatialspectral distribution.
From the above, the field distribution in the continuous space can be obtained. However, the numerical simulation needs to generate a discrete random field, that is, it needs to perform discrete Fourier transform on it. Taking a phase screen as an example, as shown in Figure 5, the square phase screen with side length L is divided into N equal parts to form N × N square grids with a spacing of ∆x. Then, the discrete Fourier transform of Equation (10) can be obtained: where, ∆k is the wavenumber interval of the phase space, according to the sampling theorem, ∆k = 2π/(N∆x).
a mean of 0 and a variance of 1. The real and imaginary parts of the phase field obta by Equation (10) are independent of each other, and both satisfy the spatial-spectra tribution. From the above, the field distribution in the continuous space can be obtained. H ever, the numerical simulation needs to generate a discrete random field, that is, it n to perform discrete Fourier transform on it. Taking a phase screen as an exampl shown in Figure 5, the square phase screen with side length L is divided into N equal to form × square grids with a spacing of ∆ . Then, the discrete Fourier transfor Equation (10) can be obtained: where, ∆ is the wavenumber interval of the phase space, according to the sampling orem, ∆ = 2 ( ∆ ) ⁄ . In order to obtain the phase distribution in real space, the inverse Fourier trans of Equation (12) can be given by: In order to obtain the phase distribution in real space, the inverse Fourier transform of Equation (12) can be given by:

Low-Frequency Sampling Optimization
For the simulated turbulent phase screen generated by the Fourier transform method, the range of the minimum wave number interval ∆k of the phase screen does not include (−∆κ x /2, ∆κ x /2) and −∆κ y /2, ∆κ y /2 , so it has the disadvantage of insufficient low frequency [31]. We use the method of adding sub-harmonics to perform low-frequency compensation on the phase screen by the Fourier transform method. In essence, the high-frequency phase screen simulated by the FFT method is interpolated and fitted to improve the statistical characteristics of low frequency. The specific process is as follows: firstly, divide the square area surrounded by the first sampling point of the low-frequency part limited to the high-frequency part into nine areas with equal area; then, set the sampling points contained in the area (including the edge part) to zero; and finally, take eight small areas outside the center as new sampling points to form a p-order (p ≥ 1) harmonic network. As shown in Figure 6, the sampling interval of the p-th harmonic network becomes ∆ f p = ∆k/3 p , that is, the original high-frequency sampling area is replaced by 3 −2p lowfrequency sampling small areas. According to the Kolmogorov spectrum, the subharmonic low-frequency compensation screen can be expressed as: S low (p ∆x, q ∆x) = a R φ φ (p ∆x, q ∆x)/∆k (15) points contained in the area (including the edge part) to zero; and finally, take eight s areas outside the center as new sampling points to form a -order ( ≥ 1) harmonic work. As shown in Figure 6, the sampling interval of the p-th harmonic network beco ∆ = ∆ 3 ⁄ , that is, the original high-frequency sampling area is replaced by 3 l frequency sampling small areas. According to the Kolmogorov spectrum, the sub monic low-frequency compensation screen can be expressed as: Therefore, the phase screen expression after sub-harmonic compensation is written The phase structure function is a quantity that describes the statistical properties o atmospheric turbulence phase screen, and the accuracy of the simulation is generally jud by a large number of phase screens. The phase structure function is defined by [32]: When the laser is transmitted in the atmospheric turbulence whose spectral m takes the Kolmogorov spectrum, the theoretical expression of the structure functio given by [32]: grid points in high frequency grid points in low frequency Therefore, the phase screen expression after sub-harmonic compensation is written by: The phase structure function is a quantity that describes the statistical properties of the atmospheric turbulence phase screen, and the accuracy of the simulation is generally judged by a large number of phase screens. The phase structure function is defined by [32]: When the laser is transmitted in the atmospheric turbulence whose spectral model takes the Kolmogorov spectrum, the theoretical expression of the structure function is given by [32]: For the inner scale, taking infinity ( κ 0 → ∞ ) and infinitely small, there are: where K 5/6 is the modified Bessel function of the third kind, and Γ(11/6) is the Gamma function.
From a statistical point of view, when giving a large number of random phase screen samples, its structure function will gradually converge to its expected structure function. In order to reduce the calculation amount, the expected structure function can be calculated directly [17]. The specific process is as follows: perform inverse Fourier transform on the phase spectral function φ φ → κ r to obtain the auto correlation function κ r , and then substituting it into Equation (17) can get its structure function. The auto correlation functions of the high-frequency phase screen and low-frequency compensation can be written by: where, f (p, q) and f (p , q ) are the spatial filter functions, f (p, q) = 2π , and r 0 is the atmospheric coherence length, which is taken as 0.185.

Scintillation Index of Laser Return Signal
For the RTS theory, the light intensity change of the lidar return spot is the research focus. Under weak turbulence conditions, the scintillation index β 2 I is one of the commonly used physical quantities to describe atmospheric turbulence. β 2 I represents the normalized light intensity fluctuation difference, which is defined as: where, < > represents the ensemble mean of the light intensity, and ρ represents the radial distance from the optical axis of the Gaussian beam. It is known that the Gaussian beam will produce spot drift after transmission in atmospheric turbulence. In order to obtain its scintillation index β 2 I (ρ i , L i ) through several light spots with a transmission distance of L i , one need to judge the drift distance ρ i of the center of mass of the light spot [33]. The spot centroid is defined as: Usually, its first moment is used to represent the coordinates of the centroid of the spot, namely [1]: x c = xI(x, y)dxdy/ I(x, y)dxdy (25) y c = yI(x, y)dxdy/ I(x, y)dxdy In the simulation, ten light spots with a transmission distance of ρ i are taken. For each light spot, if ρ i ≤ 15∆x, the center of mass of the drifted light spot is selected as the center of the circle. In addition to this, take 15∆x as the radius to make a circle, and count the light intensity of points in the circle. If ρ i > 15∆x, the center of mass of the spot after the drift is also taken as the circle center. Then, ρ i − 1.5∆x and ρ i + 1.5∆x are used as the radii to make a circle, respectively. Through counting the light intensity of N nodes in the ring, the scintillation index can be calculated by the definition formula β 2 I (ρ i , L i ) = Finally, ten statistical averages are performed to obtain the final scintillation index at a certain distance. For a Gaussian beam, according to the Rytov theory, the spot scintillation index of the beam is defined as: (27) where, k = 2π/λ, γ = (1 + iaz)/(1 + iaL) = γ r − iγ i , a = λ/πW 2 0 + i/R 0 .

Inversion of Turbulent Profile
The refractive index structure constant C n 2 can directly describe the fluctuation of the refractive index and is one of the most intuitive physical quantities to measure the intensity of atmospheric optical turbulence. The atmospheric turbulence intensity profile is the variation of atmospheric refractive index structure constant C n 2 with height [33]. In this paper, the refractive index structure constant C n 2 can be calculated from the correlation relationship between the known spherical wave scintillation index β I,S 2 and the atmospheric refractive structure constant C n 2 in the horizontal and vertical directions, so as to obtain the atmospheric turbulence intensity profile.
Kolmogorov believes that atmospheric turbulence is composed of turbulent eddies with great difference and different scales. Under the large Reynolds number Re ( Re = UL/ν , where U is the characteristic velocity of the fluid, L is the overall characteristic scale of the fluid, and ν is the molecular kinematic viscosity coefficient), the turbulent eddies of different scales coexist. After the cascade process, the small-scale turbulence finally reaches a statistical equilibrium process, that is, local isotropic turbulence [2]. The threedimensional power spectrum function of the locally isotropic turbulent refractive index fluctuation is expressed as: where, φ n (κ) is the atmospheric turbulence power spectrum function, κ is the spatial frequency, l 0 is the turbulent inner scale, that is, the scale that determines the smallest vortex in the turbulent vortex, and f (κl 0 ) is the inner scale correction model factor. The Kolmogorov atmospheric turbulence power spectrum is an ideal power spectrum model, which does not consider the inner scale effect of turbulence. It considers that the outer scale L 0 is infinite and the inner scale l o is zero [34]. At this time, f (κl 0 ) = 1. So the three-dimensional power spectrum function is expressed as: Under the condition of weak fluctuation (β I 2 < 1), the scintillation index satisfies: where, σ 2 χ is the logarithmic amplitude fluctuation variance, k = 2π/λ is the wave number, and κ is the spatial frequency.
In the vertical direction, according to the Kolmogorov spectrum, the relationship between the spherical wave scintillation index β I,S 2 and the refractive index structure constant C n 2 is given by: The atmospheric turbulence intensity profile can be inverted by the iterative method. The entire detection range 0~L km is evenly divided into several parts, and the corresponding length interval of each part is ∆L. Assuming that the refractive index structure constant C n 2 in each segment is a constant, then the C n 2 in each segment can be directly extracted from the integral equation. Here, taking Equation (32) as an example, firstly, integral the range from zero to L 1 , the corresponding spherical wave scintillation index β I,S 2 , C n1 2 within zero to L 1 can be obtained. Then, by iterating in turn, the curve of C n 2 changing with the detection height under the vertical detection path can be obtained. Therefore, using the hierarchical iteration method, Equation (32) can be converted to Equation (33).

Comparison of Low-Frequency Sampling Optimization Results
In the simulation model, the transmission of laser light in atmospheric turbulence is regarded as a series of phase screens placed in the Gaussian beam, including highfrequency phase screens and low-frequency compensated phase screens. Assuming that the turbulent spectrum structure is the Kolomogorov spectrum, the corresponding phase screen is constructed by the Fourier transform method, here, C n 2 = 2 × 10 −15 , L = 0.5 m, N × N = 256 × 256, ∆z = 200 m, and the transmission distance L = 2000 m. Figure 7 shows the phase distribution of the high-frequency phase screen numerically simulated in the turbulent atmosphere.

Comparison of Low-Frequency Sampling Optimization Results
In the simulation model, the transmission of laser light in atmospheric turb regarded as a series of phase screens placed in the Gaussian beam, including quency phase screens and low-frequency compensated phase screens. Assumin turbulent spectrum structure is the Kolomogorov spectrum, the correspondi screen is constructed by the Fourier transform method, here, = 2 × 10 ， ， × = 256 × 256，∆ = 200 m, and the transmission distance = 2000 m 7 shows the phase distribution of the high-frequency phase screen numerically s in the turbulent atmosphere.  In view of the shortcomings of the lack of low-frequency information in the phase screen constructed by the Fourier transform method, the low-order harmonic compensation is performed. According to the steps mentioned above, the schematic diagram of lowfrequency compensation for numerical simulation of the phase distribution in a turbulent atmosphere is shown in Figure 8. It can be seen that the transition of the phase distribution is smoother. Figure 9 shows the comparison of the Kolmogorov phase screen structure function. It can be seen that the phase distribution obtained by the spectral inversion method is consistent in the high-frequency part, but obviously insufficient in the low-frequency part. However, the phase distribution after the third harmonic compensation is significantly improved in the low-frequency part. tion is performed. According to the steps mentioned above, the schematic diagram o frequency compensation for numerical simulation of the phase distribution in a turb atmosphere is shown in Figure 8. It can be seen that the transition of the phase distrib is smoother. Figure 8. The three-dimensional (left) and two-dimensional (right) schematic diagrams of t merically simulated phase distribution in the turbulent atmosphere after the third harmoni pensation under the same conditions. Figure 9 shows the comparison of the Kolmogorov phase screen structure fun It can be seen that the phase distribution obtained by the spectral inversion meth consistent in the high-frequency part, but obviously insufficient in the low-frequency However, the phase distribution after the third harmonic compensation is signifi improved in the low-frequency part.   Figure 9 shows the comparison of the Kolmogorov phase screen struc It can be seen that the phase distribution obtained by the spectral inversi consistent in the high-frequency part, but obviously insufficient in the low-fr However, the phase distribution after the third harmonic compensation is improved in the low-frequency part.  Dr r/m harmonic compensation with three harmonic compensation with zero theory Figure 9. Kolmogorov phase screen structure function comparison diagram.

Result Verification of Turbulent Phase Screen Simulation Model
In order to test the reliability of the results of the turbulent phase screen simulation model, the scintillation index is chosen as the physical quantity to test this method. In order to simulate the transmission of the laser beam in atmospheric turbulence, we select the Gaussian distribution of the laser beam to pass through the phase screen shown in Figure 8. Figure 10 shows the light intensity distribution of the laser beam passing through the atmospheric turbulence on the vertical path. It can be seen that the phase disturbance generated by the Gaussian beam after passing through the turbulent phase screen makes the light intensity distribution no longer satisfy the Gaussian distribution. When the distance between phase screens (∆z) remains the same, the discrete degree of the spot increases with the increase in the transmission distance (L); when the transmission distance (L) remains the same, the discrete degree of the spot becomes larger with the increase in the distance between the phase screen (∆z). order to simulate the transmission of the laser beam in atmospheric turbulence, we the Gaussian distribution of the laser beam to pass through the phase screen show Figure 8. Figure 10 shows the light intensity distribution of the laser beam passing thr the atmospheric turbulence on the vertical path. It can be seen that the phase distur generated by the Gaussian beam after passing through the turbulent phase screen m the light intensity distribution no longer satisfy the Gaussian distribution. When th tance between phase screens (Δz) remains the same, the discrete degree of the sp creases with the increase in the transmission distance (L); when the transmission dis (L) remains the same, the discrete degree of the spot becomes larger with the incre the distance between the phase screen (Δz).
When the conditions are 2 For the light intensity distribution of the laser transmission on the vertical pat spot drift phenomenon occurs. Figure 11 shows the spot drift phenomenon. In Figu the hexagon indicates the original spot centroid located in the middle, while the tri shows the centroid of the drifted spot. It is clear that the larger the value of the p screen spacing Δz is, the greater the discrete degree of the spot and the more obviou phenomenon of light spot drifts. For the light intensity distribution of the laser transmission on the vertical path, the spot drift phenomenon occurs. Figure 11 shows the spot drift phenomenon. In Figure 11, the hexagon indicates the original spot centroid located in the middle, while the triangle shows the centroid of the drifted spot. It is clear that the larger the value of the phase screen spacing ∆z is, the greater the discrete degree of the spot and the more obvious the phenomenon of light spot drifts. In this paper, the scintillation index of different heights is simulated and 10 sp ages of the simulation results for each height are employed to obtain the correspo scintillation index from the definition formula for each spot image. The above resu statistically averaged to obtain the scintillation index simulation value representin height, and then compared with the Rytov theoretical results as shown in Figure 12 be seen that the change trend of the theoretical value and the simulation value is ba consistent. In this paper, the scintillation index of different heights is simulated and 10 spot images of the simulation results for each height are employed to obtain the corresponding scintillation index from the definition formula for each spot image. The above results are statistically averaged to obtain the scintillation index simulation value representing each height, and then compared with the Rytov theoretical results as shown in Figure 12. It can be seen that the change trend of the theoretical value and the simulation value is basically consistent. , the distance between the phase screens (Δz) and the transmission distance (L) are changed, respectively, and the spot drift of the laser beam passing through the phase screen are shown in (a) and (b).
In this paper, the scintillation index of different heights is simulated and 10 spot images of the simulation results for each height are employed to obtain the corresponding scintillation index from the definition formula for each spot image. The above results are statistically averaged to obtain the scintillation index simulation value representing each height, and then compared with the Rytov theoretical results as shown in Figure 12. It can be seen that the change trend of the theoretical value and the simulation value is basically consistent.

Atmospheric Turbulence Profile
The Hufnagel-Valley 5/7(HV 5/7 ) model is a common statistical model representing the turbulence profile in the mid-latitude inland region, and the expression is as follows [35]: (34) Figure 13 shows the profile of the atmospheric refractive index structure constant obtained by the iterative inversion algorithm. Compared with the input HV 5/7 model, it can be seen that below the detection height of 4500 m, the C n 2 obtained by inversion is slightly smaller than that of the HV 5/7 model. When the detection height is 4500~6500 m, the C n 2 inverted is slightly larger than that of the HV 5/7 model. Therefore, the atmospheric refractive index structure constant profile obtained by the inversion algorithm is consistent with the input model below 6500 m, which preliminarily verified that the Mie-scattering lidar improved by the RTS theory has certain reliability in detecting atmospheric turbulence.

Atmospheric Turbulence Profile
The Hufnagel-Valley 5/7(HV5/7) model is a common statistical model representing th turbulence profile in the mid-latitude inland region, and the expression is as follows [ Figure 13 shows the profile of the atmospheric refractive index structure constan obtained by the iterative inversion algorithm. Compared with the input HV5/7 model, can be seen that below the detection height of 4500 m, the obtained by inversion slightly smaller than that of the HV5/7 model. When the detection height is 4500~6500 m the inverted is slightly larger than that of the HV5/7 model. Therefore, the atmo pheric refractive index structure constant profile obtained by the inversion algorithm consistent with the input model below 6500 m, which preliminarily verified that the Mi scattering lidar improved by the RTS theory has certain reliability in detecting atmo pheric turbulence.

Conclusions
Numerical simulation is the most commonly used method to study atmospheric tu bulence, and it can be used to verify and evaluate the atmospheric turbulence detectio capability of lidar systems. The Fourier transform method is one of the methods in nu merical simulation. In this paper, using RTS theory combined with a Mie-scattering lida system in the laboratory, a phase screen simulation experiment for detecting atmospher turbulence was designed. The Gaussian laser beam transmission in atmospheric turbu lence was simulated by the Fourier transform method with the addition of the third ha monic, and a series of phase screens were used to simulate the phase disturbance and ligh intensity scintillation caused by atmospheric turbulence. The scintillation index can b calculated from the light intensity distribution of the phase screen. From the relationshi between the scintillation index and the atmospheric refractive index structure constant i the vertical direction, the simulated atmospheric turbulence profiles are obtained by ite ative inversion, and the accuracy of the phase screen construction is tested by the Ryto theory. The simulation results are shown as follows:

Conclusions
Numerical simulation is the most commonly used method to study atmospheric turbulence, and it can be used to verify and evaluate the atmospheric turbulence detection capability of lidar systems. The Fourier transform method is one of the methods in numerical simulation. In this paper, using RTS theory combined with a Mie-scattering lidar system in the laboratory, a phase screen simulation experiment for detecting atmospheric turbulence was designed. The Gaussian laser beam transmission in atmospheric turbulence was simulated by the Fourier transform method with the addition of the third harmonic, and a series of phase screens were used to simulate the phase disturbance and light intensity scintillation caused by atmospheric turbulence. The scintillation index can be calculated from the light intensity distribution of the phase screen. From the relationship between the scintillation index and the atmospheric refractive index structure constant in the vertical direction, the simulated atmospheric turbulence profiles are obtained by iterative inversion, and the accuracy of the phase screen construction is tested by the Rytov theory. The simulation results are shown as follows: (1) Compared with the traditional Fourier transform method, the Fourier transform method with the addition of the third harmonic makes the phase distribution of the phase screen smoother, and it makes up for the shortcomings of the traditional method of insufficient low-frequency distribution; (2) Comparing the spot scintillation index of Gaussian beam obtained by the Rytov theory, the trend of the scintillation index obtained by simulation is basically consistent with the theoretical value. The scintillation index gradually increases with the increase in the measurement distance, and it can be seen that the constructed phase screen is accurate; (3) Comparing the atmospheric turbulence profile of the HV5/7 model, the atmospheric refractive index structure constant profile obtained by iterative inversion is in good agreement with the HV5/7 model below 6500 m; it fluctuates with the increase in the detection distance, generally increases and tends to be stable.
This study preliminarily verifies the reliability of the Mie-scattering lidar system for detecting atmospheric turbulence from a theoretical point of view, and provides a lowcost and easy-to-operate measurement method for detecting atmospheric turbulence in the future. In fact, the simulation of atmospheric turbulence has practical applications, such as the determination of the height of the top of the atmospheric boundary layer and the study of the wind load characteristics of the atmospheric boundary layer turbulence on low-rise buildings. At the same time, it provides a feasible research method for indepth understanding of the direct and indirect effects of aerosols. Further research on the inversion algorithm is required in the later stage.  Data Availability Statement: The aerosol particles size distribution and concentration data used to support the findings of this study are available from the corresponding author upon request.