Evaluation of middle atmosphere temperature and wind measurements and their disturbance characteristics by meteorological rockets

It is necessary to carry out the in-situ detection based on the meteorological rocket to deepen the cognitive level of the middle atmosphere environment. However, the effective utilization of rocket data has been limited by a lack of systematic research into data accuracy and the physical mechanisms influencing measurements. In this study, temperature and wind profiles between 20 and 60 km were obtained over northwest China using two meteorological rockets equipped with thermistors and BeiDou positioning, supplemented by a temperature correction technique. The detection results are compared with satellite, empirical model and reanalysis data. An error analysis was performed by integrating the characteristics of the drop sounding process and atmospheric disturbances. The results indicate that the rocket-derived data are of satisfactory quality, with altitude-dependent trends in temperature and wind profiles consistent with those from other sources. Observed discrepancies come from the deviation of the matching data in time and space and the excessive measurement error in the initial fall stage. Also, it is found that the instability of the parachute causes poor positioning data quality and fast falling speed, leading to significantly larger measurement errors at corresponding altitudes. Additionally, the profile from the first detection exhibited more pronounced fluctuations, attributable to the breaking of high-altitude gravity waves. The dissipation of these waves reduces atmospheric stability and generates denser small-scale layered structures on the profile, making significant wind field changes at lower altitudes through the momentum deposition.

1 Introduction

The near space is located in the region of 20–100 km, which can cover the stratosphere, the mesosphere and the low thermosphere. The near-space atmosphere is far from the Earth's surface and lacks the weather phenomena common in the troposphere (cyclones, thunderstorms, fronts, etc.). Nevertheless, its unique significance continues to attract considerable research attention. First, the near space is the upper boundary of the troposphere, which can be coupled with the troposphere and exert a top-down influence. Due to its relatively slow evolution compared with the troposphere, the stratosphere provides valuable information for predicting extreme tropospheric weather and climate (Gray et al., 2018; Jin et al., 2023). For example, the weakening of the stratospheric polar vortex is often a precursor to the occurrence of cold waves in the Northern Hemisphere. Second, the near space is the lower boundary atmosphere of space weather, which can act as a “display screen” for solar activity, and the influence of solar activity on Earth's weather and climate can be reflected in it. For example, solar activity can alter the ozone concentrations in the middle atmosphere, with these changes subsequently transmitted to the troposphere via planetary waves (Krivolutsky et al., 2015). In addition, the near space is the combination of aerospace and aviation, and changes in the internal environment will directly affect the flight attitude and effect of aerospace vehicles (Chen et al., 2023; Roney, 2007; Sheng et al., 2025). Atmospheric disturbances, as the superposition of waves at different scales (including turbulence, gravity waves, planetary waves, etc.), are one of the main dynamic processes in the near space (He et al., 2025, 2024, 2023). As atmospheric density decreases exponentially with height, the amplitude of upward-propagating disturbances such as gravity waves increases progressively, leading to more pronounced wave-driven effects (Lindzen, 1981; Alexander et al., 2010).

Given the growing significance of the near-space region, there is an urgent need to enhance our understanding of its internal atmospheric environment. The necessary condition to support this demand is to carry out accurate detection and adequate research. Satellite remote sensing provides atmospheric profile data with global coverage. However,its capability for wind field detection remains limited, and the vertical resolution of data is coarse (Ern et al., 2022; Thies and Bendix, 2011). Lidar and MST (Meso-Stratosphere-Troposphere) radars can retrieve three-dimensional wind fields and temperature profiles. Nevertheless, the global distribution of detection sites is limited, and the data quality is affected by atmospheric environment and retrieve accuracy (She et al., 2003; Lu et al., 2018; Qiao et al., 2020). Zero-pressure or super-pressure balloons (often referred to as constant-level or flat-floating balloons) enable continuous horizontal sampling within the stratosphere. However, the characteristics of its own drift in the wind bring the uncertainty of detection, and require strict trajectory control technology (He et al., 2024; Alexander et al., 2021). Radiosonde balloons can detect meteorological elements with long time series and high precision. However, the highest detection height is generally less than 30 km, and cannot cover higher airspace (He et al., 2022; Yoo et al., 2020). In contrast, the meteorological rocket sounding is the only in-situ detection method that can obtain the atmospheric environment in the altitude range of 20–100 km. The effective evaluation and inspection of rocket detection accuracy is an important prerequisite for the correct use of this means.

The meteorological rocket sounding primarily includes two methods: falling spheres detection and thermistor detection. The falling spheres can obtain the atmospheric density profile between 30 and 100 km altitude, from which wind fields, temperature, and pressure are derived. The thermistor measurement can obtain the atmospheric temperature from 20 to 60 km, and then calculate the density, pressure, and wind field (Eckermann et al., 1995; Wang et al., 2006). Due to the large amplitude of atmospheric waves in this height range, the momentum and energy dissipated by wave fragmentation can cause drastic changes in meteorological elements such as wind field, density and temperature in the surrounding atmosphere. Consequently, analyzing the interaction mechanism between atmospheric wave and background flow constitutes a major research focus utilizing in-situ observational data. By comparing with satellite, balloon and reanalysis data, thermistor rockets launched from Hainan Station and East China Sea have shown good detection results, and the atmospheric disturbance characteristics in near space are also extracted (Guoying et al., 2011; Song et al., 2024). Atmospheric density is measured using GPS data on a rigid falling ball and the measured deviation from the model results was less than 10 % (Yuan et al., 2017). Using passive ball falling experiments in northwest China, in-situ wind field and gravity wave information are analyzed from 30 to 100 km (Ge et al., 2019). A comprehensive evaluation of the detection accuracy of the TK-1 meteorological rocket is performed and the reliability is demonstrated (Fan et al., 2013). It can be seen that the current results of near space rocket detection are still few, encouraging researchers to work in greater depth.

In this paper, two meteorological rockets launched in the northwest of China are used to obtain meteorological detection data from 20 to 60 km, error analysis and accuracy evaluation are carried out, and wave disturbance characteristics are also extracted. The structure of the paper is as follows: in the second section, the used data is introduced; in the third section, the temperature correction and error calculation method are given; in the fourth section, the comparison results of rocket detection profile and reference data are discussed, in the fifth section, the error analysis is performed; in the sixth section, the characteristics of wave perturbations and their effects on the background atmosphere are discussed; in the seventh section, the conclusion and prospect are given.

2 Rockets instrument and detection principle

The rocket-borne radiosonde system integrates several key subsystems and components, including temperature sensors, pressure sensors, satellite navigation and positioning modules, data acquisition circuits, transmitters, wireless remote control modules, batteries, switches, fixed frames, insulation boxes and fiberglass reinforced plastic shells, etc. The temperature sensor adopts a bead thermistor, purchased from the shelf, model MF51MP-D (Blue Crystal Electronics). The pressure sensor adopts a high-precision digital pressure sensor, purchased from the shelf, model ms5607 (Switzerland). The navigation and positioning module adopts the high-precision positioning module of Beidou, and the antenna uses a four-arm helical antenna, which is a customized product. The main MCU of the data acquisition circuit adopts a 32 bit processor with ARM core, featuring low power consumption and mixed signal processing capabilities. It has a 14 bit A/D conversion accuracy, which can meet the measurement accuracy requirements of sensors. The digital transmitter is composed of dedicated RF chips and power amplifier modules to form a frequency point digital transmitter. It has the advantages of small size and adjustable frequency. When used in conjunction with ground receiving equipment, it can achieve data transmission within a diagonal distance range of 200 km. The physical configuration and internal layout of the rocket sounding instrument are shown in Fig. 1, and the main performance indicators of the rocket sounding instrument are shown in the Table A1.

Figure 1 (a) The physical configuration and (b) internal layout of the rocket sounding instrument.

Figure 2 Meteorological rocket detection mechanism.

The rocket detection mechanism is shown in Fig. 2. The meteorological sonde is carried up by the rocket, under the action of thrust, it rises at a high speed according to the established trajectory. After the engine stops working, the rocket uses inertia to continue rising. When the rocket rises near the top of its trajectory, the parachute carries the sonde and separates from the arrow body. The sonde pulls the parachute and begins to fall.

During this process, the atmospheric parameters are measured in situ and the data is transmitted down to the ground-based receiving system. The atmospheric temperature profile from 20 to 60 km is directly measured using a thermistor sensor. Atmospheric pressure is then derived through an iterative calculation initialized with a base pressure measurement (obtained by a pressure sensor at ∼20 km). Then the atmospheric density is calculated through the ideal gas equation. The real-time position coordinates (X, Y, Z) of the sonde are obtained by using the Beidou positioning system, and the first derivative is obtained by linear fitting after the smoothing position coordinates point by point to calculate the northward, eastward and vertical velocity (represented by $\dot{x}$, $\dot{y}$, and $\dot{z}$). The corresponding acceleration is obtained by quadratic fitting (represented by $\ddot{x}$, $\ddot{y}$, and $\ddot{z}$). Based on the velocity and acceleration information, the meridional, zonal, and synthetic wind are calculated (represented by W_x, W_y, and W), and the wind direction (θ) can be further obtained. The specific calculation formula are given in Eqs. (1)–(4)

$$\begin{array}{}\text{(1)}& {\displaystyle }{W}_x& {\displaystyle }=\dot{x}-{\displaystyle \frac{\ddot{x}}{\ddot{z}-g}}\dot{z}\text{(2)}& {\displaystyle }{W}_y& {\displaystyle }=\dot{y}-{\displaystyle \frac{\ddot{y}}{\ddot{z}-g}}\dot{z}\text{(3)}& {\displaystyle }W& {\displaystyle }=\sqrt{W_x^{\mathrm{2}}+W_y^{\mathrm{2}}}\end{array}$$

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathit{\theta }{\displaystyle }=\left\{\begin{array}{ll}\arctan \left|\frac{W_y}{W_x}\right|+\mathrm{180}\mathit{°},& (W_x>\mathrm{0},W_y>\mathrm{0})\\ -\arctan \left|\frac{W_y}{W_x}\right|+\mathrm{180}\mathit{°},& (W_x>\mathrm{0},W_y<\mathrm{0})\\ -\arctan \left|\frac{W_y}{W_x}\right|+\mathrm{360}\mathit{°},& (W_x<\mathrm{0},W_y>\mathrm{0})\\ \arctan \left|\frac{W_y}{W_x}\right|,& (W_x<\mathrm{0},W_y<\mathrm{0})\end{array}\right.\end{array}$$

The air pressure at each height layer is calculated from the measured base point air pressure (20 km) using the pressure height formula:

$$P=P_{\mathrm{d}}\exp {\displaystyle \frac{-g_{\mathrm{0}}(H-H_{\mathrm{d}})}{R{T}_{\mathrm{d}}}}$$

Among them, P represents the air pressure at the calculated height, P_d is the air pressure of the adjacent lower layer, H is the geopotential, H_d is the geopotential of the adjacent lower layer, R is the dry air gas constant, and T_d is the temperature of the adjacent lower layer. Given the temperature and air pressure, the atmospheric density can be calculated through the ideal gas state equation.

The specific calculation process of atmospheric parameters is shown in Fig. 3.

Figure 3 Atmospheric parameters calculation process.

3 Data introduction

This study analyzes data from two meteorological rocket soundings conducted in northwestern China in autumn 2023. The primary dataset comprises vertical profiles (20–60 km, the effective analysis interval) of atmospheric temperature, pressure, density, and wind (including synthetic, zonal, and meridional components, as well as wind direction).

The Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) carried on TIMED satellite, can obtain the vertical profile of atmospheric temperature, pressure, geopotential height, ozone and other trace gas volume mixing ratio by limb scanning. For this study, we utilize SABER Version 2.0 Level 2A temperature data.

The Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) is a second-generation, high-precision global atmospheric reanalysis dataset. The data has a temporal resolution of 6 h and contains 42 pressure layers ranging from 1000 to 0.1 hPa. The data used in this paper are zonal wind, meridional wind and atmospheric temperature data. The spatial resolution of the original data was 0.5°×0.625°.

The NRLMSISE-00 atmospheric empirical model provides reference atmospheric states from the ground to the thermosphere (0–1000 km). The input parameters include the solar and geomagnetic activity index, date, latitude, longitude, altitude and local time, and the output elements are the temperature and density profile of the neutral atmosphere.

To facilitate comparison with the rocket sounding data, the reanalysis, empirical model, and satellite datasets must be spatiotemporally collocated. The verification data close to the time (<5 h) and within a certain deviation range of latitude (<4°) and longitude (<4°) are selected and interpolated to the same vertical grid points as the data processed by the rocket.

4 Temperature correction and error calculation

4.1 Temperature correction

During the process of parachute fall, thermistor and the outside atmosphere has been a heat exchange, in unit time, thermistor internal energy ΔE, self-heating L, convection exchange heat H, radiation exchange heat Q, viscous exchange heat M, lead conduction heat exchange N have the following relationships (Wagner, 1964):

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }E=L+H+Q+M+N,\end{array}$$

According to the modified formula given by the World Meteorological Organization on the temperature detection data of the rocket sonde, Eq. (5) is expanded to (Organization, 2008):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{T}_{\mathrm{\infty }}& =T_f-{\displaystyle \frac{r{v}_r^{\mathrm{2}}}{\mathrm{2}{c}_p}}+{\displaystyle \frac{m_{T}C}{Ah}}{\displaystyle \frac{\mathrm{d}{T}_f}{\mathrm{d}t}}-{\displaystyle \frac{A_m{\mathit{\rho }}_m{\mathit{\alpha }}_{s}J}{Ah}}\\ & -{\displaystyle \frac{{\mathit{\alpha }}_t\mathit{\sigma }(A_a{T}_a^{\mathrm{4}}+A_b{T}_b^{\mathrm{4}}+A_c{T}_c^{\mathrm{4}})}{Ah}}+{\displaystyle \frac{\mathit{\epsilon }\mathit{\sigma }{T}_f^{\mathrm{4}}}{h}}\\ & -{\displaystyle \frac{Q_c}{Ah}}-{\displaystyle \frac{W_f}{Ah}}\end{array}\end{array}$$

Where T_f is the original temperature, T_∞ is the temperature after correction. The heating term $\frac{r{v}_r^{\mathrm{2}}}{\mathrm{2}{c}_p}$ reflects the influence of heat exchange between the thermistor and its boundary layer on the temperature indication value. The temperature hysteresis term $\frac{m_{T}C}{Ah}\frac{\mathrm{d}{T}_f}{\mathrm{d}t}$ represents the influence of the hysteresis of thermistor heat exchange on the temperature indication value. The reflected radiation term $\frac{A_m{\mathit{\rho }}_m{\mathit{\alpha }}_{s}J}{Ah}$ represents the influence of the short-wave solar radiation reflected by the ground and clouds to the sonde on its temperature indication. The long wave radiation term $\frac{{\mathit{\alpha }}_t\mathit{\sigma }(A_a{T}_a^{\mathrm{4}}+A_b{T}_b^{\mathrm{4}}+A_c{T}_c^{\mathrm{4}})}{Ah}$ represents the influence of radio frequency radiation and infrared radiation in the environment of the sonde on the temperature indication. The external radiation term $\frac{\mathit{\epsilon }\mathit{\sigma }{T}_f^{\mathrm{4}}}{h}$ represents the influence of the thermal radiation of the sensor to the sonde on its temperature indication. The structural heat conduction term $\frac{Q_c}{Ah}$ represents the influence on the thermistor indication due to the thermal conduction of the sonde support to the thermistor. Measuring current heating term $\frac{W_f}{Ah}$ indicates the amount by which the temperature indication of the resistance changes due to the heating of the current. The sonde takes shading measures to ignore the direct solar radiation. The meanings of each item in Eq. (6) are shown in Table A1.

4.2 Error calculation

Temperature measurement error is composed of thermistor static calibration error σT₁, temperature error caused by position error σT₂, and temperature correction error ΔT₃ (Wagner, 1964, 1961), the calculation formula is as follows:

$$\begin{array}{}\text{(7)}& \mathit{\delta }T=\sqrt{\mathit{\sigma }{T}_{\mathrm{1}}^{\mathrm{2}}+\mathit{\sigma }{T}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{\Delta }{T}_{\mathrm{3}}^{\mathrm{2}}},\end{array}$$

σT₁ and σT₂ are the systematic errors of the instrument, which are fixed values in calculation, ΔT₃ is the residual error after temperature correction (Eq. 6), and the formula is calculated as:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathrm{\Delta }{T}_{\mathrm{3}}& =\mathrm{\Delta }(T_{\mathrm{\infty }}-T_f)\\ & =\mathrm{\Delta }\left(-{\displaystyle \frac{r{v}_r^{\mathrm{2}}}{\mathrm{2}{c}_p}}\right)+\mathrm{\Delta }\left({\displaystyle \frac{m_{T}C}{Ah}}{\displaystyle \frac{\mathrm{d}{T}_f}{\mathrm{d}t}}\right)\\ & +\mathrm{\Delta }\left(-{\displaystyle \frac{A_s{\mathit{\alpha }}_{s}J}{Ah}}\right)+\mathrm{\Delta }\left(-{\displaystyle \frac{A_m{\mathit{\rho }}_m{\mathit{\alpha }}_{s}J}{Ah}}\right)\\ & +\mathrm{\Delta }\left(-{\displaystyle \frac{{\mathit{\alpha }}_t\mathit{\sigma }(A_a{T}_a^{\mathrm{4}}+A_b{T}_b^{\mathrm{4}}+A_c{T}_c^{\mathrm{4}})}{Ah}}\right)\\ & +\mathrm{\Delta }\left({\displaystyle \frac{\mathit{\epsilon }\mathit{\sigma }{T}_f^{\mathrm{4}}}{h}}\right)+\mathrm{\Delta }\left(-{\displaystyle \frac{Q_c}{Ah}}\right)+\mathrm{\Delta }\left(-{\displaystyle \frac{W_f}{Ah}}\right),\end{array}\end{array}$$

Wind speed error is composed of systematic error and random error. Systematic error is written as:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{l}\mathrm{\Delta }{W}_x=\mathrm{\Delta }\dot{x}-\frac{\dot{z}}{\ddot{z}-g}\mathrm{\Delta }\ddot{x}-\frac{\ddot{x}}{\ddot{z}-g}\mathrm{\Delta }\dot{z}+\frac{\ddot{x}\dot{z}}{(\ddot{z}-g{)}^{\mathrm{2}}}\mathrm{\Delta }\ddot{z}\\ \mathrm{\Delta }{W}_y=\mathrm{\Delta }\dot{y}-\frac{\dot{z}}{\ddot{z}-g}\mathrm{\Delta }\ddot{y}-\frac{\ddot{y}}{\ddot{z}-g}\mathrm{\Delta }\dot{z}+\frac{\ddot{y}\dot{z}}{(\ddot{z}-g{)}^{\mathrm{2}}}\mathrm{\Delta }\ddot{z}\end{array}\right.,\end{array}$$

Random error is written as:

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{l}\begin{array}{rl}{\mathit{\sigma }}_{W_x}^{\mathrm{2}}& ={\mathit{\sigma }}_{\dot{x}}^{\mathrm{2}}+{\left({\displaystyle \frac{\dot{z}}{\ddot{z}-g}}{\mathit{\sigma }}_{\ddot{x}}\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{\ddot{x}}{\ddot{z}-g}}{\mathit{\sigma }}_{\dot{z}}\right)}^{\mathrm{2}}+{\left[{\displaystyle \frac{\ddot{x}\dot{z}}{(\ddot{z}-g{)}^{\mathrm{2}}}}{\mathit{\sigma }}_{\ddot{z}}\right]}^{\mathrm{2}}\end{array}\\ \begin{array}{rl}{\mathit{\sigma }}_{W_y}^{\mathrm{2}}& ={\mathit{\sigma }}_{\dot{y}}^{\mathrm{2}}+{\left({\displaystyle \frac{\dot{z}}{\ddot{z}-g}}{\mathit{\sigma }}_{\ddot{y}}\right)}^{\mathrm{2}}\\ & +{\left({\displaystyle \frac{\ddot{y}}{\ddot{z}-g}}{\mathit{\sigma }}_{\dot{z}}\right)}^{\mathrm{2}}+{\left[{\displaystyle \frac{\ddot{y}\dot{z}}{(\ddot{z}-g{)}^{\mathrm{2}}}}{\mathit{\sigma }}_{\ddot{z}}\right]}^{\mathrm{2}}\end{array}\end{array}\right.,\end{array}$$

g is the gravity acceleration, $\mathrm{\Delta }\dot{x}$, $\mathrm{\Delta }\dot{y}$ and $\mathrm{\Delta }\dot{z}$ are velocity fitting deviations, $\mathrm{\Delta }\ddot{x}$, $\mathrm{\Delta }\ddot{y}$, and $\mathrm{\Delta }\ddot{z}$ are acceleration fitting deviations, ${\mathit{\sigma }}_{\dot{x}}$, ${\mathit{\sigma }}_{\dot{y}}$, and ${\mathit{\sigma }}_{\dot{z}}$ are speed random errors, ${\mathit{\sigma }}_{\ddot{x}}$, ${\mathit{\sigma }}_{\ddot{y}}$, and ${\mathit{\sigma }}_{\ddot{z}}$ are acceleration random errors.

The total error of wind speed and direction is calculated as follows:

$$\begin{array}{}\text{(11)}& \left\{\begin{array}{l}\mathit{\delta }{W}_{\mathit{\epsilon }}=\sqrt{\mathit{\delta }{W}_{x\cdot \mathit{\epsilon }}^{\mathrm{2}}+\mathit{\delta }{W}_{y\cdot \mathit{\epsilon }}^{\mathrm{2}}}\\ \mathit{\delta }G=\frac{\mathrm{180}}{\mathit{\pi }}\sqrt{{\left(\frac{W_x\cdot \mathit{\delta }{W}_{y\cdot \mathit{\epsilon }}}{W_x^{\mathrm{2}}+W_y^{\mathrm{2}}}\right)}^{\mathrm{2}}+{\left(\frac{W_y\cdot \mathit{\delta }{W}_{x\cdot \mathit{\epsilon }}}{W_x^{\mathrm{2}}+W_y^{\mathrm{2}}}\right)}^{\mathrm{2}}}\end{array}\right.,\end{array}$$

Where $\mathit{\delta }{W}_{x\cdot \mathit{\epsilon }}=\sqrt{{\mathit{\sigma }}_{W_x}^{\mathrm{2}}+\mathrm{\Delta }{W}_x^{\mathrm{2}}}$ is the meridional wind synthesis error, and $\mathit{\delta }{W}_{y\cdot \mathit{\epsilon }}=\sqrt{{\mathit{\sigma }}_{W_y}^{\mathrm{2}}+\mathrm{\Delta }{W}_y^{\mathrm{2}}}$ is the zonal wind synthesis error.

Figure 4 Time-altitude curves of (a) HJ-1 and (b) HJ-2, and horizontal motion trajectory of (a) HJ-1 and (b) HJ-2 (red for rocket ascent, blue for sonde/parachute drift).

5 Comparison of rocket detection results with reference data

5.1 Data quality and trajectory analysis

The two rockets are referred to as HJ-1 and HJ-2, respectively. HJ-1 is launched at 09:00 UTC on the first day, and HJ-2 is launched at 05:00 UTC on the next day.

The time-altitude profiles of HJ-1 and HJ-2 are shown in Fig. 4 (top). The actual detection altitude of HJ-1 is about 74 km, the ascent time is about 2 min, and the fall time (from the highest point to an altitude of 20 km) is 25 min. HJ-2 can reach a maximum altitude of 76 km, the ascent time is about 2 min, and the fall time is 31 min. Taking the launch point as the central point, the horizontal motion trajectory of the ascending stage of rocket launch and the sonde/parachute drift stage are plotted as shown in Fig. 4 (below). During ascent, both rockets traveled predominantly eastward. Following apogee, the sondes continued to drift eastward during their parachute-assisted descent. This consistent eastward motion at all altitudes is attributable to the prevailing westerly background wind field in the region, as evident from the trajectories. The sonde remained within 100 km from the launch point during the entire detection process (from the beginning of the launch to the 20 km falling height).

In order to further analyze the trajectory characteristics of the sonde during its fall, the vertical distribution of zonal velocity (v_x), meridional velocity (v_y) and vertical velocity (v_z) are shown in Fig. 5. The zonal velocity of the two rockets is positive, and the meridional velocity gradually changes from positive to negative, which corresponds to the characteristics of the falling trajectory drifting first to the northeast and then to the southeast in Fig. 4. It is worth noting that there is an obvious disturbance characteristic (denser small scale layered structure) of vertical velocity for HJ-1, compared with that of HJ-2. After the same data processing method, the obvious difference of v_z profile roughness may reflect the great difference of disturbance in the vertical direction at high altitudes.

Figure 5 Velocity-altitude curves of (a) HJ-1 and (b) HJ-2.

5.2 Wind and temperature measurements

Figure 6 shows the comparison of the zonal and meridional winds obtained by the two rockets with the MERRA2 data. Before the launch of the rocket, the balloon sounding is also carried out. Within the balloon's altitude range, the rocket-derived wind speed and direction closely match the balloon profile, with consistent disturbance details (Fig. A1). This agreement validates the reliability of the rocket-retrieved wind fields. The meridional winds of the two rockets both reach the maximum value near 50 km, exceeding 40 m s⁻¹. As reflected in Fig. 4, in the initial stage of fall after the rocket body-parachute separation, the trajectory turns from south to north, which proves that the strong meridional winds dominate at high altitudes.

Figure 6 (a) The vertical distribution of zonal wind and meridional wind of HJ-1, (b) the difference of HJ-1 velocity component with MERRA2, (c) the vertical distribution of zonal wind and meridional wind of HJ-2, and (d) the difference of HJ-2 velocity component with MERRA2.

Here, deviation refers to the difference between the rocket sounding data and a reference dataset (e.g., satellite or reanalysis). It quantifies the consistency between the rocket measurements and other data sources. HJ-1 and MERRA2 have basically the same variation trend of wind speed components at the altitude of 20–60 km, and the zonal wind deviation is relatively small in the whole altitude, while the meridional wind deviation has large positive and negative fluctuations between 50–60 km. The mean absolute deviations (MAD) of the zonal wind at the whole altitude is 3.3 m s⁻¹, and that of the meridional wind is 5.4 m s⁻¹. In contrast, HJ-2 wind profiles show greater fluctuation than that of MERRA2, with deviations increasing markedly above 40 km. The corresponding MAD of zonal wind is 7.5 m s⁻¹ and that of meridional wind is 7.6 m s⁻¹. In the altitude range of 20–45 km, the variation trend of wind speed is consistent. At higher altitudes, the measured wind speed of the rocket can show more significant fluctuation characteristics. There are maximum wind speed areas near 30 and 55 km for both the two rockets, and the maximum near 55 km is difficult to reflect in the MERRA2 data. This suggests that the reanalysis data may have insufficient observational constraints for assimilation at higher altitudes, and the difference of wind field in the upper stratosphere is obviously greater than that in the lower stratosphere even in the close spatiotemporal range. Considering that the output from the model tends to reflect the average trend, and the transient results of a single detection are more prominent, it is reasonable to have differences between the rocket detection and the model.

Figure 7 shows the vertical distribution of temperature from rocket, SABER, MSISE, and MERRA2 data and the corresponding deviation from them. Results before and after temperature correction and corresponding sub-term correction amount are shown in Fig. A2, the temperature correction is larger above 50 km, and gradually decreases below 50 km. Among the various correction sub-items for rocket detection temperature, the influence degree of pneumatic heating, current heating, and temperature hysteresis are relatively large, and these influences gradually decrease with decreasing altitude. The net effect of these corrections is a reduction in temperature across the entire profile. According to the maximum temperature, the stratopause height measured by the rocket (the height of the inflection point) is around 47 km. The stratopause height is consistent with other reference data for HJ-1, but shows some differences for HJ-2. The temperature profiles of the four data have a consistent trend from 20 to 50 km, with small deviation. The deviation between the reference data and the rocket detection results increases above 50 km. Within this interval, the temperature deviation between HJ-1 and MISIS is the smallest, while the difference between HJ-2 and SABER is the smallest. It is worth noting that the temperature deviation of HJ-1 increases sharply above 57 km, a feature potentially attributable to measurement error (discussed further in Sect. 6). The difference of data comparison may be due to the following reasons: (1) There are deviations in the position and time of the reference data matching with the rocket; (2) The results of the model reflect the average over time and space, which is indeed different from the single-point profile.

Figure 7 (a) The vertical distribution of temperature for HJ-1, (b) the difference of HJ-1 temperature with MERRA2, MSISE, and SABER, (c) the vertical distribution of temperature for HJ-2, and (d) the difference of HJ-2 temperature with MERRA2, MSISE, and SABER.

The vertical distribution of rocket detection density and their relative deviations with MERRA2, MSISE, and SABER are also shown in Fig. A3. The density relative deviation of HJ-1 shows a significant maximum value between 40 and 50 km (the deviation can reach about 10 % for the above three reference data), while the relative density deviation of HJ-2 is significantly smaller, especially with excellent consistency with the SBAER data (the relative deviation within the entire detection height range is within 5 %). The large density deviation of HJ-1 in the upper atmosphere is, on the one hand, due to the significant reduction in the density itself, which makes the difference more prominent. On the other hand, it is very likely that there are other strong atmospheric disturbances causing drastic changes in density (discussed later), which have not been captured by the model and satellite data.

6 Error analysis

Accurate measurement is the prerequisite for conducting further data analysis and application. To investigate the sources of discrepancy between the rocket detection results and other data, as well as the reliability of the disturbance analysis, the error level of the rocket instrument is discussed here, which is a comparison among different heights and sub-items within the rocket own detection results. Temperature and wind measurement errors of HJ-1 and HJ-2 can be obtained according to Eqs. (8) and (11), as shown in Fig. 8. Systematic and random errors of wind speeds are shown in Figs. A3 and A4, respectively. The regional-mean atmospheric temperature errors for HJ-1 is 0.31, 0.53, and 5.5° for the 20–30, 30–50, and >50 km altitude bins, respectively. For HJ-2, the corresponding errors are 0.24, 0.55, and 1.75°. Similarly, the wind speed errors for HJ-1 is 0.63, 1.12 and 4.95 m s⁻¹ across the same altitude bins; the values for HJ-2 are 0.38, 1.19 and 4.0 m s⁻¹. The wind direction error levels of HJ-1 are 0.81, 1.08, and 3.15° at 20–30, 30–50, and above 50 km, respectively, while that of HJ-2 is 0.54, 1.11, and 4.25°. According to Eqs. (9) and (10), when the vertical acceleration and vertical velocity are too large, the denominator $\ddot{z}-g$ decreases and the numerator $\ddot{z}$ increases, which can obviously affect the results of systematic error and random error. In the whole detection section, the same smooth fitting points are used, so the velocity error is consistent. However, due to the large jump of the positioning data, the acceleration ratio in the inertial velocity will also jump. When the falling velocity is large, the product will also increase, resulting in a significantly larger error margin at the high altitudes.

Figure 8 Error-height curves of (a) temperature, (b) wind speed, and (c) wind direction for HJ-1 and HJ-2.

The original temperature vertical gradient and vertical acceleration of HJ-1 and HJ-2 are shown in Fig. 9. During the initial descent phase (50–60 km) following parachute deployment, the high descent velocity is accompanied by significant acceleration fluctuations. For HJ-1, the two vertical acceleration peaks within 50–60 km correspond directly to maxima in both the systematic and random wind error profiles (Fig. A3). Similarly, for HJ-2, the rapid increase in vertical acceleration above 50 km is mirrored by a concurrent rise in its wind speed errors. According to the error equation, the measurement error of wind speed depends largely on the velocity error and acceleration error. At the same time, the temperature error is also related to the vertical gradient of the measured temperature indication value (in HJ-1, the obvious gradient deviation above 58 km and its ratio to the convective heat exchange coefficient cause the temperature error to increase sharply), which is also the reason why the temperature error and wind field error in Fig. 8 have inconsistent trends. Therefore, the magnitude of vertical acceleration fluctuation is a primary determinant of wind field measurement error. At high altitudes (near to 60 km), the parachute swing is large, and the data reception is not stable, resulting in the relatively low positioning data quality and the large position error, which finally lead to the relatively large wind field error. As the detection height gradually decreases, the positioning data quality increases and the measurement error decreases gradually as the parachute falls steadily.

Figure 9 (a) Temperature gradient-height curve and (b) vertical velocity-altitude curve for HJ-1, (c) temperature gradient-height curve and (d) vertical velocity-altitude curve for HJ-2.

7 Disturbance characteristic analysis

In the accuracy analysis of the two rocket detection results, we find that compared with HJ-2, HJ-1 has a more intense falling velocity disturbance, and the deviation from the reference data is larger. Profile deviation phenomena discovered in the result above are inferred to be closely related to the strong disturbance at this height. Therefore, it is necessary to further verify these phenomena in the atmosphere through disturbance characteristic analysis. Conducting wave disturbance analysis here, on the one hand, verify the previous detection results, and on the other hand, it is also an application study of rocket detection data, enhancing the theoretical nature and completeness of the rocket data analysis results.

7.1 Wave energy and background field analysis

Due to a scarcity of measured wind field data at high altitudes (30–60 km), a detailed understanding of fine-scale wind disturbances at corresponding interval remains limited. Many sharp peaks in the wind profile captured by balloon and rocket detections represent real atmospheric perturbations (Figs. 6 and A1), which are smoothed out in the reanalysis. Consequently,, using rocket data may be more suitable for analyzing wave disturbance characteristics at high altitude, since reanalysis data failed to capture these details. The apparent differences in vertical velocity and acceleration of the sonde during its fall (Fig. 9) also indicate significant differences in upper atmospheric disturbances. By analyzing the atmospheric background state and gravity wave (GW) information, we compare the difference characteristics of atmospheric disturbance in two detection processes.

Figure 10 (a) Wind shear, (b) square buoyancy frequency, (c) Richardson number and (d) kinetic energy of HJ-1, and (e) wind shear, (f) square buoyancy frequency, (g) Richardson number and (h) kinetic energy of HJ-2.

GWs are generated by the excitation source at the lower atmosphere, and their amplitudes increase gradually as the atmospheric density decreases during upward propagation. Wind shear is an important disturbance source of high-altitude GWs, which can cause GWs to be generated or broken (Larsen, 2002; Larsen and Fesen, 2009). Vertical wind shear is calculated as:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{d}U}{\mathrm{d}z}}=\sqrt{{\left({\displaystyle \frac{\mathrm{d}u}{\mathrm{d}z}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{d}v}{\mathrm{d}z}}\right)}^{\mathrm{2}}},\end{array}$$

Buoyancy frequency N indicates atmospheric stability, with N²>0 and N²<0 i denoting statically stable and unstable conditions, respectively. It is calculated from:

$$\begin{array}{}\text{(13)}& N^{\mathrm{2}}={\displaystyle \frac{g}{T}}\left[\left({\displaystyle \frac{\mathrm{d}T}{\mathrm{d}Z}}\right)+{\displaystyle \frac{g}{c_p}}\right],\end{array}$$

The gradient Richardson number R_i represents the ratio of buoyancy to shear production terms, and is defined as

$$\begin{array}{}\text{(14)}& R_i={\displaystyle \frac{N^{\mathrm{2}}}{{\left(\frac{\mathrm{d}u}{\mathrm{d}z}\right)}^{\mathrm{2}}+{\left(\frac{\mathrm{d}v}{\mathrm{d}z}\right)}^{\mathrm{2}}}},\end{array}$$

Atmospheric GWs can be regarded as superimposed disturbances to the background field. First, a 20-point sliding average is performed on the profile interpolated with equal spacing (50 m interval) to eliminate errors caused by random motion and turbulence. Then the smoothing profile is fitted by fifth-order polynomial to get the background profile. After the background profile is removed, high-pass filtering with a cut-off wavelength of 10 km is performed to obtain the disturbance profile caused by GWs. The kinetic energy E_k of GW is calculated by the following formula:

$$\begin{array}{}\text{(15)}& E_k=\mathrm{1}/\mathrm{2}(u^{{\prime }^{\mathrm{2}}}+v^{{\prime }^{\mathrm{2}}}),\end{array}$$

Where u^′ and v^′ are the disturbance components of the zonal and meridional wind field caused by GWs, respectively.

In the error analysis, considering that the error becomes significant above 55 km (Fig. 8), the height interval selected for disturbance analysis here is 20–55 km. The vertical distribution of wind shear, square buoyancy frequency, Richardson number and kinetic energy obtained according to HJ-1 and HJ-2 detection results are shown in Fig. 10. The wind shear of HJ-1 has the first peak (strongest) near 45–55 km and the second peak near 30–40 km, while the wind shear peak of HJ-2 is between 30–40 km. N² is positive throughout the profile, indicating overall static stability. However, it shows a general tendency to decrease with increasing altitude. HJ-1 has a buoyancy frequency minimum (even close to 0) between 45 and 55 km, corresponding to large wind shear, resulting in a relatively concentrated area of R_i<0.25, indicating strong dynamic instability. In contrast, HJ-2 has a smoother profile with smaller wind shear and larger buoyancy frequency, resulting in fewer dynamic instability regions. For HJ-1, the peak kinetic energy of GW is above 50 km, corresponding to the maximum value region of wind shear, and the dynamic instability region is relatively concentrated, indicating that Kelvin–Holtzmann instability has produced strong high-altitude wave disturbance. Below 50 km, the GW energy of HJ-1 is significantly smaller than that of HJ-2, which is mainly due to the attenuation of zonal wave disturbance (Fig. A5).

7.2 Spectral analysis

Lomb–Scargle spectrum analysis is applied to the disturbance profile of the synthesized wind speed, and the vertical wave-number spectrum caused by GWs are obtained, as shown in Fig. 11. The amplitudes of GWs in HJ-1 are significantly weaker than those in HJ-2, and the vertical wavelengths of the dominant GWs (amplitudes greater than 90 % confidence) are more dispersed, with scales ranged from 1.7 to 7.1 km present. In contrast, the dominant GWs in HJ-2 have stronger amplitudes and are concentrated at wavelengths around 4–6 and 2.9 km.

Figure 11 Gravity wave information for (a) HJ-1 and (b) HJ-2 obtained from the disturbance profile of the wind field, dashed lines represent 90 % confidence, and dominant wavelengths with amplitudes above this threshold are labeled.

Combined analysis of atmospheric instability and GW spectra indicates that the atmospheric disturbance for HJ-1 is more complex. The GW breaks, resulting in enhanced turbulent activity (more dynamic unstable regions), which also leads to a significant reduction in stratification stability (reduced buoyancy frequency) with more small-scale stratification (Held et al., 2019; van Haren et al., 2015). The GW kinetic energy can be reduced and the amplitude corresponding to the dominant wavelength decreases. Consequently, compared with HJ-2, the measured temperature and wind field profile of HJ-1 have more obvious fluctuations and a denser small-scale layered structure. Furthermore,, the GW kinetic energy of HJ-1 is significantly lower than that of HJ-2 in the range of 40–50 km (Fig. 10), which is considered to be the main region where wave dissipation occurs. Consistent with this energy depletion, both zonal and meridional winds for HJ-1 are weaker than for HJ-2 between 40–50 km, while the profiles are similar below 40 km (Fig. 6), which further indicates that wave dissipation weakens the local winds. Thus, the breaking and dissipation of GWs in the upper stratosphere can reasonably explain the difference of detection profiles in adjacent two days.

7.3 Wave dissipation revealed from Stokes parameter method and ERA5 results

To obtain further evidence for GW breaking in the 40–50 km layer during HJ-1, Stokes parameter method (Vincent and Fritts, 1987; Eckermann., 1996) is used here to extract the typical characteristic parameters of the GW. The main realization path is as follows: Fourier transform is applied to the zonal wind and meridional wind disturbances, and corresponding real and imaginary parts are obtained respectively. Then four Stokes parameters I, D, P, and Q are calculated, and information such as scale, propagation and frequency of polychromatic gravity waves can be further obtained. The specific method can be referred to the previous paper (He et al., 2022).

Given that the wave breaking primarily occurs below 50 km, the GW parameters are calculated for the two height intervals of 40–50 and 20–50 km, corresponding to disturbance information in the local and entire height range, respectively. The kinetic energy, horizontal wavelength, intrinsic frequency, vertical group velocity and horizontal propagation direction extracted from the two detections are shown in Table 1. For a local wave disturbance (40–50 km), there is a low-frequency GW of HJ-1, with an intrinsic frequency (the ratio of wave frequency to inertial frequency) of 2.53. The order of wavelength, kinetic energy and vertical group velocity is within a reasonable range. In contrast, the intrinsic frequency and vertical group velocity of HJ-1 are abnormally large, while the horizontal wavelength is abnormally small, which should belong to the omitted cases. The outliers of the characteristic parameter also reflect the breaking of GWs in this region from the perspective of abnormal high frequency waves (Fritts and Alexander, 2003), meaning that GWs can no longer maintain their normal state and dissipate. For the entire wave disturbance (20–50 km), HJ-2 has no obvious wave breaking, and the parameters such as wavelength and frequency are close to the local disturbance, which means a consistent wave propagation process throughout the entire height. In contrast, the wavelength and kinetic energy of the entire wave disturbance of HJ-1 are smaller than that of HJ-2 due to local wave breaking. The wave propagation direction of HJ-2 is significantly different in the entire and local ranges, possibly due to significant wind speed changes near 40 km (Fig. 6c).

Table 1 Gravity wave parameters extracted by Stokes parameter method.

Figure 12 Regional distribution of ERA5 vertical velocity (w) at 3 hPa for (a) 03:00 UTC, (b) 06:00 UTC, and (c) 09:00 UTC on the first day, and (d) 00:00 UTC, (e) 03:00 UTC, and (f) 06:00 UTC on the second day, with the five-pointed star representing the rocket detection position. The purple arrow represents the direction in which the wave travels, and the purple rectangular box represents the region where the wave dissipates occurs. The launch point of rocket is (X° E, Y° N).

Although the soundings are separated by one day, the continuous but slow process of GW momentum deposition into the mean flow (Liu et al., 1999) suggests that a comparison of the two wind profiles may still reveal GW drag effects. For the local wave breaking of HJ-1 (40–50 km), the propagation direction is northwest (the degrees in Table 1 represent angle measured anticlockwise from x axis), and the deposited momentum produces negative drag (deceleration) on the zonal wind and positive drag (acceleration) on the meridional wind. Compared with the earlier detection, significantly stronger meridional wind and significantly weaker zonal wind can be seen near 40 km in the later detection (Fig. 6). This suggests that the drag effect of local wave breaking through deposited momentum is captured at an altitude of 40 km by HJ-2, and the acceleration of tens of meters per day is also reasonable (Li et al., 2022).

In order to further support the wave breaking at high altitude during HJ-1 detection, ERA5 data is used to plot the longitude-latitude cross section of vertical velocity at 3 hPa (near 41 km) in the corresponding region (10° longitude × 5° latitude), as shown in Fig. 12. HJ-1 detection is close to 09:00 UTC (launched at 09:50 UTC) on the first day (T), and HJ-2 detection is close to 06:00 UTC (launched at 5.5 UTC) on the second day (T+1). In this pressure layer, the vertical velocity (w) has an obvious alternating positive and negative perturbation, which indicates the GW activity. For the first day of detection, at 03:00 UTC and 06:00 UTC, the northwestward movement of the GW is observed. At 09:00 UTC, there is a distinct wave breaking (purple rectangular box). For the second day of detection, at 03:00 and 06:00 UTC, the southwestward movement of the perturbation peaks can be observed, and no wave dissipation occurs in the corresponding region. Both the time of wave dissipation and the direction of wave propagation are consistent with the results calculated by Stokes parameter method from rocket data (Table 1), which further proves the reliability of the results.

8 Summary

In this study, the detection effects and data quality of two meteorological rocket launched in the northwest of China in the autumn of 2023 are analyzed. First, using the modified temperature correction model and wind field retrieval algorithm, the atmospheric temperature, pressure, density, wind speed and wind direction measured by the rocket are obtained, and compared with the matched reanalysis, satellite and empirical model data. Second, measurement errors and the overall data accuracy are quantified through error propagation and synthesis analysis, and their impacts are assessed. Finally, the characteristics of atmospheric instability and GW activity are analyzed. The main conclusions are as follows:

1. Both soundings achieved high data acquisition rates with normal, smooth trajectories during ascent and descent, constituting successful experiments that yielded good-quality meteorological profiles from 20 to 60 km altitude.

2. Rocket-derived wind fields agree well with MERRA-2 below 40 km, with deviations increasing at higher altitudes. Rocket temperature profiles also show good agreement with MERRA-2, NRLMSISE-00, and SABER below 50 km, beyond which deviations grow. These biases likely originate from spatiotemporal mismatches in the data collocation and the fundamental difference between model/reanalysis averages and instantaneous in-situ measurements.

3. Below 50 km, the wind measurement error and temperature measurement error remain small (<2 m s⁻¹ and <1.8°, respectively). Errors increase above 50 km, attributable to significant parachute oscillation and unstable data reception during initial descent, which degrade positioning accuracy and sensor-environment coupling

4. The difference in the intensity of GWs causes the obvious difference in vertical velocity of the dropsonde. For HJ-1, the amplitude of GWs over this region is reduced, and turbulent activity is enhanced, resulting in reduced stability of atmospheric stratification and a denser small scale hierarchical structure on the profile. For HJ-2, the stratification stability of the upper atmosphere is stronger. GWs are more stable and less likely to break, allowing the amplitude to grow to a larger degree.

5. The local breaking of GWs at 40–50 km can be captured ideally from HJ-1. The GWs deposited momentum and energy to the mean flow, and the effect of the wave drag changed the wind field structure below, making HJ-2 with one day delay can detect significant wind field changes near the altitude of 40 km. This reflects the forcing effect of wave dissipation on the background wind field through the observation results. The results of ERA5 data further support the wave dissipation and propagation characteristics extracted by rocket data.

The analysis demonstrates that rocket drop sounding, with its high vertical resolution and in-situ nature, can capture atmospheric fine structure close to its true state. Accurate and detailed wind field results are very valuable, especially in the region above 30 km. The large measurement error above 55 km also indicates that it is necessary to improve the data reception quality at the beginning of the drop, and optimize the high-altitude parachute opening and stability control technology to improve the detection accuracy. The presence of atmospheric GWs causes the local feature difference of the detection profile, meaning that the high-altitude disturbance characteristics need to be considered in the rocket detection. This study provides observational support for GW dissipation theories in the upper stratosphere – a region where such detailed in-situ evidence has been scarce. More rocket soundings across different regions and seasons are encouraged to further improve our understanding of the near-space atmospheric environment.

Appendix A

Table A1 Main performance indicators of rocket radiosondes.

Table A2 Variable meaning in the Eq. (6).

Figure A1 Comparison of wind speeds measured by rockets and balloons for (a) HJ-1 and (b) HJ-2, and comparison of wind directions measured by rockets and balloons for (c) HJ-1 and (d) HJ-2.

Figure A2 The vertical distribution of (a) original and corrected temperature, and (b) each correction subterm.

Figure A3 The vertical distribution of rocket detection density for (a) HJ-1 and (c) HJ-2. The relative deviation of rocket detection density with MERRA2, MSISE, and SABER for (b) HJ-1 and (d) HJ-2.

Figure A4 Random error of (a) meridional wind and (b) zonal wind, and systematic error of (c) meridional wind and (d) zonal wind for HJ-1.

Figure A5 Random error of (a) meridional wind and (b) zonal wind, and systematic error of (c) meridional wind and (d) zonal wind for HJ-2.

Figure A6 (a) Zonal wind and (b) meridional wind disturbance profile caused by GWs for HJ-1, and (c) zonal wind and (d) meridional wind disturbance profile caused by GWs for HJ-2. The solid and dashed lines represent the mean and standard deviation over the entire height, respectively.