A Simulation Experiment to Retrieve the Atmospheric Density and Three-Dimensional Wind Field by Double Falling Spheres

: Falling-sphere sounding remains an important method for in situ determination in the middle atmosphere and is the only determination method within the altitude range of 60–100 km. Traditional single-falling-sphere sounding indicates only the atmospheric density and horizontal wind but not the vertical wind; the fundamental reason is that the equation set for retrieving atmospheric parameters is underdetermined. For tractability, previous studies assumed the vertical wind, which is much smaller than the horizontal wind, to be small or zero. Obtaining vertical wind profiles necessitates making the equations positive definite or overdetermined. An overdetermined equation set consisting of six equations, by which the optimal solution of density and three-dimensional wind can be obtained, can be established by the double-falling-sphere method. Hence, a simulation experiment is designed to retrieve the atmospheric density and three-dimensional wind field by double falling spheres. In the inversion results of the simulation experiment, the retrieved density is consistent with the constructed atmospheric density in magnitude; the density deviation rate does not generally exceed 20% (less than 5% below 60 km). The atmospheric density retrieved by the double-falling-sphere method is more accurate at low altitudes than the single-falling-sphere method. The vertical wind below 50 km and horizontal wind retrieved by double-falling-sphere method is highly consistent with the constructed average wind field. Additionally, the wind field deviation formula is deduced. These results establish the fact that the double-falling-sphere method is effective in detecting atmospheric density and three-dimensional wind.


Introduction
The middle atmosphere, which refers to the region extending from more than 10 km to approximately 100 km above the ground, includes the stratosphere, the mesosphere and the lower thermosphere [1,2]. In recent years, numerous studies of this thin and neutral atmosphere have been carried out [3]. He et al. [4] explored the spectral characteristics of temperature fluctuations and threedimensional wind field fluctuations by a set of near-space high-resolution balloon data and increased the height range of spectral analysis to 38 km. He et al. [5] analyzed the scale interactions between the small-scale gravity wave and turbulence in the middle stratosphere. Sheng et al. [6] applied Thorpe analysis to COSMIC (Constellation Observing System for Meteorology, Ionosphere, and Climate) data to retrieve the strongest mixed layer in the troposphere (SMLT) altitude.
Wind is one of the key parameters of atmospheric dynamics [7]. The neutral wind field not only plays an important role in the energy transmission and atmospheric dynamics but also has a great impact on the safety and orbits of spacecraft along with atmospheric density and temperature [8]. The understanding of the middle atmospheric wind field is of great significance to the study of dynamics and behaviors of the middle atmosphere and the development of forecasting capabilities for weather and space environments [9]. However, there are still few measurement methods for the wind field of the middle atmosphere at present. During the past decade, horizontal winds have been measured from space for the first time from 30 to 90 km using submillimetric limb sounding, and the Swedish Space Agency will launch the Stratospheric Inferred Winds (SIW) satellite in 2023 [10]. Baron et al. [10] conducted a simulation study to assess the measurement performance of SIW. Additionally, the main passive optical instruments used for neutral wind measurement include the Fabry-Perot interferometer [11], Michelson interferometer [12] and Doppler asymmetric spatial heterodyne spectrometer interferograms [13]. Using meteor radar [14], lidar [15,16], medium frequency radar [17], etc. can also measure atmospheric parameters such as the wind field in the middle atmosphere. The measurement methods mentioned above are all remote sensing methods, but in situ measurement methods are still necessary for small-scale or short-term middle atmosphere study.
The falling sphere sounding method is an important measurement method for the middle atmospheric parameters and the only in situ measurement method at altitudes ranging from 60 to 100 km [18,19]. Bartman et al. [20] used the trajectory data of an inflatable sphere with a small antenna and a transponder to retrieve the atmospheric density and temperature profiles successfully. Otterman et al. [21] described the theory of using the falling sphere to retrieve the upper atmosphere density and horizontal wind field. Faucher et al. [22] obtained atmospheric density, pressure and temperature profiles over 88-128 km using an inflatable sphere carrying a tri-axial accelerometer. Yuan et al. [23] obtained atmospheric density, temperature and horizontal wind field profiles by GPS data from rigid falling spheres. Ge et al. [24] used the radar tracking data of the first passive falling sphere experiment in China's northwest region to retrieve the horizontal wind field and extracted wind shear and gravity waves from the horizontal wind field.
In previous studies using falling spheres to retrieve atmospheric parameters, vertical wind is usually ignored [20,25] or assumed to be a small value [23,26] since it is very small compared with horizontal wind and difficult to measure [27]. Although this can ensure the completeness of the equations for retrieving atmospheric parameters, it also introduces new errors at the same time.
The vertical wind field plays a significant role in the vertical transport of momentum, kinetic energy and gravity potential energy [27]. However, few previous studies have discussed the specific methods of obtaining vertical wind using the falling sphere. Jones and Peterson [28] mentioned that vertical wind might be measured by the trajectories of two spheres with different mass-to-area ratios, but to date, no study or experiment on the inversion of vertical wind with double falling spheres has been conducted. Therefore, in this paper, the theory of the double-falling-sphere method is described in detail, a simulation experiment using double falling spheres to retrieve the three-dimensional wind field and atmospheric density is designed to verify the feasibility of this method, and the error of the retrieved three-dimensional wind field is discussed.
In this paper, Section 2 introduces the data used in the simulation experiment of retrieving the atmospheric density and three-dimensional wind field by double falling spheres. Section 3 describes the theory and experimental design of calculating the three-dimensional wind field using double falling spheres, and the inversion results are given in Section 4. The causes and formula of wind field deviation are discussed in detail in Section 5. Finally, our conclusions are given in Section 6.

Data
The neutral atmospheric density and temperature profile data provided by the NRLMSISE-00 (NRL: US Naval Research Laboratory, MSIS: Mass Spectrometer and Incoherent Scatter radar, E: Exosphere) model are used as the simulated reference atmospheric density and temperature. The three-dimensional average wind field of the atmosphere required for the simulation experiment is constructed based on the atmospheric horizontal wind profile data provided by the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2).

NRLMSISE-00 Model
The NRLMSISE-00 atmosphere model, which is currently one of the most widely used atmospheric models, can provide physical quantities such as neutral atmospheric density and temperature for the region from 0 to 1000 km [29]. This paper uses the neutral atmospheric density and temperature profile data provided by the NRLMSISE-00 model as the simulated reference atmospheric density and temperature. Additionally, the upper atmosphere temperature of the model is used as the initial atmosphere temperature in the inversion process.

MERRA-2 Atmospheric Reanalysis Data Set
MERRA-2 is a global atmospheric reanalysis dataset that includes various meteorological parameters, such as net radiation, temperature, relative humidity, and wind speed. In this paper, the wind field data of MERRA-2 are used as the simulated horizontal wind field from 20 to 60 km.

Theory
A single sphere is mainly affected by the combined action of gravity, buoyancy, the Coriolis force, and atmospheric drag during the falling process [16,22,24], which can be expressed as where ⃑ is the total acceleration of the falling sphere, ⃑ is the acceleration of gravity, ⃑ is the acceleration of buoyancy, ⃑ is the acceleration of Coriolis force, and ⃑ is the acceleration of atmospheric drag. The most important item in Equation (1) for retrieving the atmospheric wind field is the drag acceleration term, which can be expressed as where is the atmospheric density, is the mass of the sphere, is the cross-sectional area of the sphere, ⃑ is the velocity of the sphere relative to air (airspeed), which equals the ground speed ( ⃑) minus the wind speed ( ⃑), and is the drag coefficient. The drag coefficient is related to the Reynolds number and Mach number [20], that is, = ( , ). In this paper, the empirical formula for the drag coefficient summarized by Henderson [30] is used, and the relationship between and ( , ) is shown in Figure 1. The Reynolds number and Mach number can be expressed as [23,24] where is the dynamic viscosity of the flow medium, is the characteristic length of the falling sphere (the sphere diameter), is the ratio of specific heat, is the gas constant for dry air, and is the atmospheric temperature. Hence, the drag coefficient is related to atmospheric density and temperature in the final analysis. To obtain the drag coefficient, it is necessary to know the temperature, and the inversion method of temperature is discussed in detail by Yuan et al. [23]. It is worth noting that the temperature data used to calculate the drag coefficient is not entirely derived from the NRLMSISE-00 atmospheric model. Only at the first inversion height, the temperature value from the NRLMSISE-00 model is needed. The temperature inversion method is detailed in Appendix A. The vector equation Equation (1) contains 4 unknowns: (zonal wind speed), (meridional wind speed), (vertical wind speed), and . Gravity, the Coriolis force, and the acceleration of the sphere can all be calculated from the trajectory of the falling sphere, and then the equations to solve these 4 unknowns are: where the subscripts x, y and z represent the latitudinal, meridional and vertical components of the vector respectively, and is the sphere volume.
where the subscripts (1) and (2) represent the physical quantities of the different falling spheres.
To ensure that Equation set (6) is overdetermined, the mass-to-area ratios of the two falling spheres should be different. The zonal wind speed and meridional wind speed can be easily solved by proportional relationships from the first three equations of the Equation set (6) if the atmospheric density and vertical wind speed can be determined. In this paper, the brainstorming optimization algorithm (BSO, described in Appendix B) is used to find optimal and values, and then the proportional relationship is used to obtain and , so that the target function ( , ) can be minimized. The target function is defined as where is the weight coefficient and the value range of is 0~1.

Simulation Experimental Design
The simulation experiment mainly consists of three parts: constructing the atmospheric parameter profile, forward-modeling the trajectories of double falling spheres, and retrieving the atmospheric three-dimensional wind field and density profile based on the trajectory data of double spheres. Figure 2 shows the flow chart of the simulation experiment. Since it is almost impossible for the two spheres to pass through the same spatial position at the same time during the falling process and since the time for the two falling spheres to reach the same height is also different, there must be certain variation in the wind field in different spaces and times. Therefore, it is necessary to add random deviation to the wind field in the process of the forwardmodeling of trajectories of the double spheres.
Since the trajectories of the double spheres are inevitably inconsistent, the two spheres actually reflect the wind field and atmospheric density at different spatial positions. To reflect the average condition of the three-dimensional wind field and atmospheric density over small scales of time and space, the mass-to-area ratios of the two falling spheres should be different but not too different. If the two mass-to-area ratios are too different, the trajectories of the two spheres are far apart, and the retrieved three-dimensional wind field and atmospheric density cannot represent the average condition of a small temporal and spatial scale. The horizontal wind field from 20 to 60 km provided by the MERRA-2 data set was used as the simulated corresponding height horizontal wind field, while the horizontal wind field from 60 to 100 km and the vertical wind field were artificially set. Some researchers have pointed out that the vertical wind speed can reach several meters per second at heights of 30-50 km [31] and even reach approximately 25 m/s at heights of approximately 80 km [28]. Based on these conclusions, the simulated three-dimensional wind field ( Figure 4) was used as the background three-dimensional wind field within a small temporal and spatial scale.

Forward-Modeling Trajectories of Double Falling Spheres
This paper simulates the trajectories of two falling spheres with different mass-to-area ratios in the simulated atmospheric environment mentioned above. The sphere with a smaller mass-to-area ratio was used as the datum sphere, while the other sphere with a larger mass-to-area ratio was used as the auxiliary sphere. The trajectories of the two falling spheres in the simulated atmosphere were calculated with a time step of 0.5 s; that is, the position, velocity and acceleration of the falling spheres were calculated every half second, and the movement of the falling spheres was regarded as uniform acceleration motion in each calculation step.
The masses of the falling spheres used in the simulation experiment were both 0.28 kg, and the radii of the datum and auxiliary sphere were 0.6 m and 0.5 m. Before simulating the descending processes of the two falling spheres separately, the highest positions of the two falling spheres and the horizontal motion speeds at the highest positions need to be given as the initial condition of the simulation. When calculating the trajectories of the two spheres, the random deviation should be added to the three-dimensional wind field at each height layer. The amplitude of the random deviations is 10% of the true value of the wind speed and satisfies the random distribution rather than the normal distribution. Figure 5 shows the trajectories of the double falling spheres in which the values of x, y and h are the components of the zonal, meridional, and vertical displacements of the falling spheres, respectively, relative to the observation station.

Retrieving 3D Wind Field and Density Profile
According to the trajectories of double falling spheres obtained through forward-modeling, the speed and acceleration of the two spheres at different heights are available. Additionally, the speed and acceleration of the auxiliary falling sphere need to be interpolated to the same heights as those of the datum falling sphere, i.e., the time step (0.5 s) of the inversions is that of the calculated heights of the datum falling sphere during the forward-modeling process. Then, the optimal atmospheric density and three-dimensional wind field can be calculated by Equation set (6). Finally, the retrieved atmospheric parameter profiles need to be smoothed. Figure 6 compares the three-dimensional wind field profile calculated by the trajectories of the double falling spheres with the constructed three-dimensional simulated wind field profile. It should be noted that the simulated wind profiles in Figure 4 or Figure 6 are just individual cases. For the convenience of description, this paper carries out inversion and analysis of these individual cases. The retrieved horizontal wind field is highly consistent with the simulated horizontal wind field, especially at altitudes below 50 km, and both almost completely overlapped. There was a relatively larger deviation between the retrieved horizontal wind field and the simulated horizontal wind field at altitudes of approximately 80 km and 90 km. The retrieved vertical wind was also more accurate at altitudes below 50 km but cannot reflect the actual condition of the vertical wind field above 60 km. In the results of the wind field inversion conducted in this paper, the deviation of the retrieved zonal wind (Figure 7a) was less than 2 m/s below 60 km; the deviation of the retrieved meridional wind field (Figure 7b) was less than 2 m/s below 50 km and was approximately 3 m/s at a height of approximately 55 km; and the deviation of the retrieved vertical wind (Figure 7c) was less than 2 m/s below 58 km. In general, the deviation of the vertical wind was greater than that of the horizontal wind field, and the deviation of the wind field at high altitude was greater than that at low altitude.  Figure 8 shows the inversion results of the atmospheric density using the double-falling-sphere method, which were in agreement with the simulated atmospheric density results (within an order of magnitude). The deviation of the retrieved atmospheric density increases with decreasing altitude, but the deviation rate of density was generally not more than 20%, especially at altitudes below 60 km, where it was not more than 5%.

Comparison with the Inversion Results of the Single-Falling-Sphere Method
As a comparison, the horizontal wind field and atmospheric density profile are also retrieved by using trajectory data of a single falling sphere (datum falling sphere). The vertical wind speed was assumed to be zero during the inversion process. Figure 9 shows the deviation profile of the atmospheric parameters retrieved by the single-falling-sphere method. In general, the deviations of the horizontal wind field retrieved by the double-falling-sphere method and single-falling-sphere method were roughly the same, and the single-falling-sphere method had even higher accuracy than the double-falling-sphere method at some altitude.
The atmospheric density retrieved by the double-falling-sphere method had obviously higher accuracy than that retrieved by the single-falling-sphere method. Additionally, the atmospheric density retrieved by the single-falling-sphere method showed obvious negative deviation below 60 km and an obvious positive deviation above 80 km in Figure 9c, while the density retrieved by the double-falling-sphere method did not show obvious positive or negative deviation under 60 km in Figure 8c. Moreover, the density deviation rate of the double-falling-sphere method basically oscillates around zero below 60 km, which is much smaller than that of the single-falling-sphere method. At low altitudes, the deviation rate of the atmospheric density retrieved by the doublefalling-sphere method was less than 5%, while the maximum deviation rate of density retrieved by the single-falling-sphere method was more than 10%.

Wind Field Deviation Formula
In the inversion results of the three-dimensional wind field and atmospheric density profile by the double-falling-sphere method, both the deviation of the wind field and the deviation rate of the atmospheric density show the same characteristics: a small oscillation amplitude at low altitude (30-60 km), and a large oscillation amplitude at high altitude (60-100 km). In fact, there is a certain relationship between the deviation of the three-dimensional wind field and the deviation rate of the atmospheric density.
Suppose that the retrieved atmospheric density at a certain altitude is times the real atmospheric density, i.e., = , where is the retrieved atmospheric density and is the actual atmospheric density. The drag coefficient is related to density; therefore, the deviation of retrieved density can lead to the deviation of the drag coefficient, which results in = , where is the calculated drag coefficient in the inversion process and is the actual drag coefficient during the falling process of the sphere.
The drag force of the falling sphere is expressed as where ⃑ is the actual airspeed and ⃑ is the actual wind field. In the case of density deviation and let = , the drag force of the falling sphere can also be expressed as where ⃑ and ⃑ are the calculated airspeed and wind field in the inversion process, respectively.
Through (8) and (10), we can obtain the following Equation: Therefore, Equation (11) can be expressed as Through Equations (9) and (13), the relationship between the retrieved 3D wind field and the real 3D wind field can be expressed as Additionally, the deviation between the retrieved wind field and the actual wind field is Therefore, the deviation ∆ ⃑ is positively correlated with (the product of and ) and ⃑ (the actual airspeed). Figure 10 shows the variation curves of airspeed with altitude during the falling process. It presents the characteristics of higher airspeed at high altitudes and lower airspeed at low altitudes. Through simulation experiments, it is found that the deviation of the drag coefficient also increases with the increase in the atmospheric density deviation rate. Therefore, according to the formula in Equation (15), the wind field deviation should meet the characteristics that the deviation is smaller at lower heights but greater at higher heights. The three-dimensional wind field retrieved by the double-falling-sphere method has higher accuracy at lower altitudes. Additionally, the deviation of the calculated wind field can be effectively reduced by decreasing the mass-to-area ratios of the falling spheres to reduce the airspeed of the spheres.

Error Sources
In this simulation experiment, the retrieved parameters have a larger deviation above 60 km. Three main causes of the deviation are summarized as below: (1) The difference in the area-to-mass ratios. Due to the difference, the two spheres in the experiments pass through different spatial positions and there must be a certain difference in the atmospheric density, temperature, and wind field which affect the two spheres (this is the reason for adding random perturbations to the wind field during the forward-modeling process). Therefore, there must be deviation in the optimal atmospheric density calculated by the two spheres' trajectories. Moreover, the added random wind perturbations are more likely to cause a larger density deviation rate at high altitudes due to the thin atmosphere. According to Equation (15), the density deviation can further cause the wind field deviation. (2) The airspeeds of the spheres are very large above 60 km. It is clear that the larger the airspeed is, the larger the deviation of the retrieved wind field will be according to Equation (15). Moreover, the trajectories data of the falling spheres is sampled at 0.5 s interval, so the height interval between two adjacent sampling points above 60 km are obviously larger than that at low altitude, which leads to the inaccuracy of temperature inversion and further lead to the deviation of the drag coefficient. (3) There is no clear physical formula for calculating the drag coefficient for now. The drag coefficient can only be achieved by experimental data or empirical formulas. Due to the difficulty of obtaining the true value of the drag coefficient, it is hard to access the impacts of the deviation of the drag coefficient to the inversion results.

Conclusions
In this paper, a simulation experiment is designed to verify the feasibility of retrieving the atmospheric density and three-dimensional wind field by the double-falling-sphere method. It is remarkable that in this simulation experiment, random deviation should be added to the wind field during the forward-modeling of the trajectories of the double spheres to simulate the difference between the average wind field and the actual wind field at different time and space positions. The results of the simulation experiment show that a credible vertical wind field and more accurate atmospheric density can be obtained by using the double-falling-sphere method at lower altitudes. The conclusions are as follows: (1) The retrieved atmospheric density is consistent with the simulated atmospheric density by an order of magnitude and has higher accuracy at low altitudes. It shows the characteristics of a larger deviation rate at higher altitudes and a smaller deviation rate at lower altitudes. The deviation rate of the retrieved density is generally not more than 20% and not more than 5% below 60 km. (2) Compared with the single-falling-sphere method, the atmospheric density calculated from the double falling spheres has higher accuracy at lower altitudes. The results show that the atmospheric density retrieved by the single-falling-sphere method has obvious positive deviation or negative deviation, while the deviation rate of density retrieved by the doublefalling-sphere method oscillates around zero below 60 km. (3) The retrieved horizontal wind field is in high agreement with the simulated wind field, and the vertical wind field below 60 km is more accurate than that above 60 km. In the results of the simulation experiment in this paper, the deviation of each component of the three-dimensional wind field is generally less than 3 m/s below 60 km. In general, the retrieved three-dimensional wind field has high accuracy at altitudes below 60 km. (4) The deviation of the retrieved three-dimensional wind field is related to the deviation rate of the retrieved density and the actual airspeed. There are smaller density deviation rates and airspeeds at lower altitudes and larger density deviation rates and airspeeds at higher altitudes, so the wind field should also show the characteristics of larger deviations at higher altitudes and smaller deviations at lower altitudes according to the derived wind field deviation formula. This conclusion is consistent with the experimental results. (5) To reduce the deviation of the retrieved wind field, falling spheres with small mass-to-area ratios should be used to retrieve the three-dimensional wind field. Falling spheres with small mass-toarea ratios are more easily affected by wind, so the airspeed is smaller. According to the wind field deviation formula, the smaller the airspeed is, the smaller the deviation of the retrieved wind field will be.
where is pressure, is altitude, is atmospheric density, and is the acceleration due to gravity.
By integrating both sides of Equation (A1), we can get: where ℎ is reference altitude and is a fixed value. The Ideal gas equation of state can be expressed as = , where is the gas constant for dry air, and is the atmospheric temperature. Through Equations (A3) and (A2), we can obtain the following equation: The value of (ℎ ) is provided by the NRLMSISE-00 model. Generally, the top height of the sphere's trajectory is selected as ℎ .
At the altitude of ℎ , using the atmospheric temperature (ℎ ) provided by the NRLMSISE-00 model and brainstorming optimization algorithm, the optimal atmospheric density (ℎ ) can be solved.
At the altitude of ℎ below ℎ , the atmospheric density values between ℎ and ℎ (not include ℎ) have been calculated. By using Equation (A4) and the brainstorming optimization algorithm, the temperature value (ℎ) and the optimal atmospheric density (ℎ) can be achieved. This step can be repeated until the entire temperature profile is obtained.
Fan [17] explains the reasons for selecting the top height of the sphere's trajectory as ℎ as follows.

(A5)
The error variance is: Since the atmospheric density decays exponentially with altitude, when inversion is performed from top to bottom, the error caused by the temperature at the reference altitude will decrease as the altitude decreases.

Brief Description of Brainstorming Optimization Algorithm
The brainstorming optimization algorithm (BSO) is a swarm intelligence algorithm based on human creative thinking, which is inspired by the brainstorming process and proposed by Shi [32].
The main steps of BSO are described as follows [33,34]: Step 1. Randomly initialize N possible solutions of equations (individuals); Step 2. Classify N individuals into M groups by K-means clustering; Step 3. Calculate the fitness (such as Equation (A7) Step 4. Perform replacement operation with a certain probability. (Randomly select two groups and create a new individual, the specific process is as follows) Define a certain probability P2b (P2b = 0.5); Define a random number rand3; If rand3 < P2b: Randomly select the centers of two groups, and merger them (calculate the mean); Add random values to generate a new individual; Else: Randomly select two individuals of two groups, and merger them (calculate the mean); Step 6. Record the individual with the minimal fitness value as the best individual Bestt in this iteration.
Step 7. If the iteration number t reaches the maximum number of iterations T, terminate the program, otherwise return to step 2.