The characteristics of atmospheric media influence the propagation of laser beams, leading to a variety of optical turbulence effects^[1–3]. The strength of optical turbulence is measured by the refractive index structure parameter $C_n^2$, and the optical turbulence effects are mainly governed by the distribution of $C_n^2$ on the propagation path. Observations have revealed that $C_n^2$ varies with height, characterized by bursty extremum and distinct layers^[4–9], which, experimentally, are symbols of intermittency and stratification structures, respectively. As the shape of $C_n^2$ profiles bring influence to optical turbulence effects, it is suggested that an appropriate optical turbulence model would be one where stratification structures are superposed on a background optical turbulence. However, hindered by scarcity of theoretical models and valid statistical methods, the current assessment of optical turbulence effects usually rely on smoothed models generally ignoring the fine structures^[9–12].

In this Letter, we introduce a novel method named complete ensemble empirical mode decomposition (CEEMD) to process high-resolution optical turbulence profiles. Then, we establish a novel optical turbulence model named the intrinsical optical turbulence model (IOTM), and further evaluate the effects of intermittency and stratification structures on optical turbulence effects.

Among all of the techniques used to measure the $C_n^2$ profile (DIMM, GSM, SCIDAR, MASS and SLODAR)^[13], a radiosonde balloon equipped with microthermal and meteorological sensors has become the standard and most delicate facility. The microthermals can measure the $C_n^2$ profiles with an exceptionally high vertical resolution, which give a precise idea of the shape of the turbulent layers. The microthermal sensors used in our experiments were developed at the Key Laboratory of Atmospheric Composition and Optical Radiation, Chinese Academy of Sciences. During a flight, the microthermal sensors probe and calculate the temperature structure constant $C_T^2$$({\mathrm{K}}^2{\mathrm{m}}^{-2/3})$, using $$C_T^2=\langle {[T(x)-T(x+r)]}^2\rangle ·r^{-2/3},l_0\ll r\ll L_0,$$(1)where the triangle brackets denote an ensemble average, $r$ is the separation between two temperature probes in the position of $x$ and $x+r$, $l_0$ is the inner scale, and $L_0$ is the outer scale that define the lower and upper limit of the inertial range. The noise level of the probes in microthermal sensors corresponds to a temperature difference of 0.002 K, implying a minimum detectable $C_T^2$ of $4.0\times {10}^{-6}{\mathrm{K}}^2{\mathrm{m}}^{-2/3}$. The meteorological radiosonde measures the ambient temperature, pressure, and humidity, and GPS is used to obtain position, wind velocity, and wind direction. The detectable precision of the temperature is 0.2 K, and the precision of the pressure is 1.5 hPa. Figure 1 gives the simultaneously measured temperature, pressure, and $C_T^2$ at noon on 4 November 2015.

As a function of height $h$, optical turbulence strength $C_n^2$$({\mathrm{m}}^{-2/3})$ is an expression of absolute temperature $T$ (K), pressure $P$ (hPa), and temperature structure constant $C_T^2$^[1]: $$C_n^2(h)={[79\times {10}^{-6}\frac{P(h)}{T^2(h)}]}^2{C}_T^2(h).$$(2)

During November 2015, we launched three balloons on each sunny day within 1–2 h of local time 8, 14, and 20 o’clock in the Northwest Plateau of China. Figure 2 shows the $C_n^2$ profile at noon on 4 November 2015. One can see that $C_n^2$ fluctuates intermittently around the mean, which is averaged over a scale of 1000 m. The original data represented by gray dots illustrate the frequent fluctuations due to intermittency, and the variation of the red mean line shows the stratification structures.

The calculation of $C_n^2$, based on measured $C_T^2$, is subjected to uncertainties in pressure, temperature, and $C_T^2$. Assuming that these errors are uncorrelated with each other, the uncertainty analysis of Eq. (2) results in the estimate of the error in $C_n^2$, $${\left(\frac{\sigma }{C_n^2}\right)}^2=4{\left(\frac{{\sigma }_P}{P}\right)}^2+16{\left(\frac{{\sigma }_T}{T}\right)}^2+{\left(\frac{{\sigma }_{C_T^2}}{C_T^2}\right)}^2,$$(3)where $\sigma$ represents the error of $C_n^2$. ${\sigma }_P$, ${\sigma }_T$ and ${\sigma }_{C_T^2}$ are the uncertainties in pressure, temperature, and $C_T^2$, respectively. As shown in Fig. 1(c), $\log (C_T^2)$ varies around $-5$ at low altitudes and increases with height, staying around $-2$ in the upper air.

Figure 3 shows the relative error $\sigma /C_n^2$ in $C_n^2$ of the profile in Fig. 2, plotted against the height and $\log (C_T^2)$. The value of $\sigma /C_n^2$ decreases with the increase of $C_T^2$, and increases with the height as the temperature and pressure fall off. The uncertainty in $C_n^2$ is about 40% at low altitudes, about 10% at modest altitudes, and increases to about 30% above 25 km.

Analyzing the nonlinear and nonstationary series like $C_n^2$ profiles, traditional data-analysis methods all suffer one or more flaws, such as predetermined assumptions and the questionable approximation of ensemble mean using temporal or spatial mean. In 1998, Huang et al.^[14] proposed empirical mode decomposition (EMD) to eschew the defects of traditional approaches. Latterly, Wu and Huang^[15] developed a noise-assisted technique called ensemble EMD (EEMD) to solve the mode mixing problem of EMD, but it created new ones. Torres et al.^[16] recently upgraded EEMD to CEEMD, which includes a negligible error and less computational cost. EMD and its upgraded versions have become an established tool for time-frequency analysis of nonlinear and nonstationary signals in recent years^[17,18]. EMD decomposes a signal into several oscillating intrinsic mode functions (IMFs) by a process called the sifting algorithm. Based on EEMD, CEEMD adds white noise at each stage of the decomposition and computes a unique residue as each mode. At the end, a given signal $x$ can be represented as a sum of IMFs and a residue $R$, $$x=\sum _{i=1}^K{c}_i+R,$$(4)where $K$ is the total number of IMFs, $c_i$ the $i$th IMF, $R$ is also named as a trend term. As $\log (C_n^2)$ rather than $C_n^2$ is generally used in applications, optical or $C_n^2$ profiles refer to the spatial series of $\log (C_n^2)$ specifically hereinafter.

Figure 4 shows the CEEMD decomposition result of the $C_n^2$ profile in Fig. 2, which is decomposed into nine IMFs and a residue. The low-order IMFs contain the highest and strongest fluctuations, while the subsequent ones contain smaller and slower fluctuations.

Indeed, we have examined totally 25 $C_n^2$ flights under various meteorological conditions, which are valid after washing out the bad data, and have demonstrated the suitability of CEEMD. For all of the decompositions, 4 profiles produce 8 IMFs, 15 profiles produce 9 IMFs, and the rest produce 10 IMFs. We have found the physical meaning of every IMF, which is consistent with other studies^[19]. Though the number and shape of IMFs are different between different profiles, we find that the fourth to 8 IMFs in all of the profiles are associated with temperature gradients and wind shears. Our experiments indicate that the results of CEEMD are insensitive to the amplitude of the added noise. Besides, a large realization is suggested to get stable numbers of the modes. We set the level of noise as 0.2 multiples to standard deviation of the input signal, the number of realizations as 2000, and the maximum sifting iterations as 2000.

Wu and Huang established a statistical significance test method for IMFs against white noise^[20]. Assuming the first IMF as a referenced noise, the averaged period of each IMF and its energy density are plotted in Fig. 5. Only those lying above the 80% lines contain physical meaning. It is clear that the second and third IMFs likely contain little useful statistically significant information. Most of the fourth and all of the rest of the IMFs locate above the 80% confidence, indicating that they contain significantly quasi-cyclical variations with height.

The contribution of IMFs at different scales on the whole profile can be evaluated quantitatively by the variance contribution rate (VCR). The VCR of the $i$th IMF component ${\mathrm{VCR}}_i$ can be calculated via the formula $${\mathrm{VCR}}_i=\frac{\mathrm{var}(c_i)}{\sum _{j=1}^K\mathrm{var}(c_j)+\mathrm{var}(R)}\times 100\%.$$(5)

Table 1 lists the mean and standard deviation of quasi-cyclical scales and the VCR of IMFs. We can see that the first three IMFs maintain frequent variations with scales smaller than 113 m, contributing about 50% to whole profiles. The scales of $\mathrm{IMF}4-\mathrm{IMF}8$ range from 232 to 4244 m with little VCR. IMF9 and IMF10 are large-scale variations whose VCR is 31%.

Thus, we get the scales and statistical significance of each IMF. Judging from the statistical significance test and quasi-cyclical scale, we believe that the first three IMFs are intermittent oscillations, and the quasi-cyclical IMFs maintain variations due to stratification structures.

Based on CEEMD, we obtain an IOTM, which includes stratification structures by superimposing specific scale of IMFs. Note the $\log (C_n^2)$ series as $x$ in Eq. (4), $P_k$ as the $k$-order profile at the chosen scale of the IOTM, $$P_k=c_k+c_{k+1}+\cdots +R.$$(6)

Obviously, $P_1$ is the original profile, while $P_k$ holds intermittency in the $k$th IMF and subsequent ones. Figure 6(a) gives four intrinsical optical turbulence profiles of the original profile in Fig. 2. We can add stratification structures into the model referring to the scales in Table 1. Figure 6(b) plots a comparison between $P_7$, $P_9$, and running mean profiles. Obviously, the variation of $P_7$ is similar to that of the mean profile over 1 km, while $P_7$ corresponds to a mean profile over 10 km.

To describe wave propagation in atmosphere, one has to use a comprehensive description of the media that is as realistic as possible. Fine structures are a natural feature of the $C_n^2$ signal, which lead to a discrepancy between the averaged model and measured signal^[21]. In adaptive optics and site-testing engineering, the coherent length and isoplanatic angle are pivotal parameters to the design and evaluation of systems. For plane waves translating vertically from height $h=0$ to $h=L$, coherence length $r_0$ and isoplanatic angle ${\theta }_0$ are expressed as^[1]$$r_0={(0.423k^2{\int }_0^L{C}_n^2(h)\mathrm{d}h)}^{-3/5},$$(7)$${\theta }_0={(2.914k^2{\int }_0^L{C}_n^2(h)h^{5/3}\mathrm{d}h)}^{-3/5}.$$(8)

Here, wavenumber $k=2\pi /\lambda$, and wavelength $\lambda =0.5\times {10}^{-6}\mathrm{m}$.

Figure 7 shows the coherence length and isoplanatic angle of the $C_n^2$ profile in Fig. 2, calculating from the IOTM. The value of $r_0$ raises from 2.30 at $P_1$ to 5.48 at $P_4$, finally reaching 7.72 at $P_{10}$. The value of ${\theta }_0$ raises from 0.74 at $P_1$ to 1.31 at $P_4$, finally reaching 1.97 at $P_{10}$. Note that $r_0$ and ${\theta }_0$ increase drastically with the deletion of IMF9 to the model, emphasizing the importance of the commonly detected bulge at the tropopause for coherence length and isoplanatic angle. It is interesting that though coherence length is sensitive to low height turbulence, and the isoplanatic angle is sensitive to upper air turbulence, they both fall off when stratification structures are added into the $C_n^2$ model.

The calculations of optical turbulence effects are always related to the integrals of $C_n^2$ weighed by height, which are called turbulence moments^[11]. Here, we consider the turbulence moment $I_1$, $I_2$, and $I_3$, corresponding, respectively, to coherence length, isoplanatic angle, and scintillation^[1]: $$I_1={\int }_0^L{C}_n^2(h)\mathrm{d}h,$$(9)$$I_2={\int }_0^L{C}_n^2(h)h^{5/3}\mathrm{d}h,$$(10)$$I_3={\int }_0^L{C}_n^2(h){(L-h)}^{5/6}\mathrm{d}h.$$(11)

The contributions of different intrinsic optical turbulence profiles to $I_1$, $I_2$, and $I_3$ can be quantified by the ratio $$m_1(k)={\int }_0^L{10}^{P_k(h)}\mathrm{d}h/I_1\times 100\%,$$(12)$$m_2(k)={\int }_0^L{10}^{P_k(h)}{h}^{5/3}\mathrm{d}h/I_2\times 100\%,$$(13)$$m_3(k)={\int }_0^L{10}^{P_k(h)}{(L-h)}^{5/6}\mathrm{d}h/I_3\times 100\%.$$(14)

The blue dots in Fig. 8 show a scatter plot of $m_1$, $m_2$, and $m_3$ versus the order of intrinsic optical turbulence profiles. For convenience, the 10th IMF and the residue are amalgamated into the ninth one and can be regarded as the large-scale background profile. Generally, at a given order, $m_1$, $m_2$, and $m_3$ of the 25 different profiles are highly variable due to the diverse distribution of different profiles. However, in most of the cases, they vary in the same decreasing trend. Judging from the mean level that is represented as the red line, $m_1$ drops down to 25.6% at $P_4$, highlighting the impact of intermittent parts on propagation. $m_1$ decreases from 21.1% at $P_5$ to 15.9% at $P_8$, showing the contribution of stratification structures. $P_9$, regarded as the background profile, provides a nearly 13.4% contribution to $I_1$. $I_2$ and $I_3$ share a similar trend with a minor difference in value.

In conclusion, we affirm the feasibility of CEEMD for analyzing optical turbulence profiles. The established IOTM successfully adds fine structures into the optical turbulence model. The calculations of turbulence moments using the IOTM show the effects of intermittency and stratification. Future investigations are needed to verify the results with other measuring methods such as the differential image motion monitor.