Features of navigational radiosonde coordinate data processing for high resolution windfinding

Upper-air sounding systems (UASS) are the main means of obtaining aerological information to provide short-term and long-term weather forecasts, prevention of natural and anthropogenic catastrophes, etc. The UASS uses the principle of tracking a radiosonde launched into free flight, usually using a gas-filled balloon. Radiosondes are equipped with various sensors for measuring atmospheric parameters. The balloon, in addition to the means of the radiosonde delivery, serves as a horizontal wind sensor. At the end of the 20th century, theodolites and radars were used to track the radiosonde. In the 90s, such companies as Air, Vaisala, Graw, etc. developed UASSs using signals from the American Global Positioning System (GPS). These systems have fundamental advantages in terms of efficiency, size, mobility, radiosonde positioning accuracy and operation on mobile platforms. The use of Global Navigation Satellite System (GNSS) in upper-air sounding systems offer new perspectives of fine wind field measuring and detecting turbulent layers. The navigation radiosonde provides more detailed and high-quality information. The processing of such data has its own specificity in comparison with radar UASS. This paper discusses the coordinate data processing of a navigational radiosonde. To analyze the errors in estimating the wind speed, computer simulation of a balloon-sonde system is used. The coordinate data filter parameters and filtering efficiency was refined during the simulation.


Introduction
Currently, exactly as in the foreseeable future of mankind, the relevance of obtaining reliable, accurate and detailed information about the atmosphere state is beyond any doubt, especially in light of such issues as climate change, environmental modeling and the formation of accurate weather forecasts, warning of natural disasters, air and space exploration. As man develops the environment, develops technology and technology, requirements for the upper-air sounding systems (UASS) characteristics increase. New measurement methods are being developed and existing UASS are being upgraded.
In world practice, two basic radio engineering methods for measuring the coordinates and velocity of a radiosonde have been formed. In the first case, a radar or radio direction finder is used. The second method uses signals from radio navigational systems [1,2].  Foreign UASSs of the first type use radio direction finder for sonde tracking. These stations measure atmospheric pressure instead of distance as the equivalent of altitude. A distinctive feature of the Russian civil UASS is the implementation of the principle of a combined radio channel based on the use of a super-regenerative transceiver in the radiosonde. This allows us to use secondary radar for sounding. But in 2010-2012, with financial support from the Ministry of Education and Science of the Russian Federation, Ural Federal University (UrFU) and Radiy JSC developed the Polyus navigation UASS based on the GLONASS/GPS [3].
The radiosonde is launched on a gas-filled balloon, which is also used as a wind sensor. UASS should determine the basic parameters of the balloon motion: horizontal velocity, direction and height. Additionally, the vertical ascent velocity and spatial coordinates of the balloon, the average velocity and motion direction in the layer may be of interest. The radiosonde is tied to balloon with a cord from 10 to 30 meters long to exclude its influence on temperature and humidity sensors. The balloon friction against the air during ascent leads to the appearance of turbulent vortices in the balloon wake and fluctuations of the radiosonde reference point. These fluctuations and the balloon velocity unevenness lead to radiosonde swinging on the cord. The dynamics of these vibrations depends on the cord length.
Radiosonde oscillations cord, on the one hand, improve the sensor air blow-off and, on the other hand, make it difficult to estimate the fine structure of the wind field, since the coordinates of the radiosonde, and not the balloon, are actually measured. The accuracy and resolution of radar UASSs is worse than navigation ones. For example, the Vector-M aerologic radar computing system distinguishes between radiosonde oscillations at a distance of no more than 2 km. As the radiosonde moves away from the radar, the coordinate measurement error increases. Therefore, radar UASS averages the wind over relatively large atmospheric layers. The accuracy of navigation UASSs, even in the general mode using GNSS signals of standard accuracy, is higher than that of radar. The radiosonde navigation receiver can provide an estimate of the current velocity and coordinates every second. As a result, we obtain a conditionally instantaneous velocity, and we can also calculate the average layer velocity using coordinates. But the navigation receiver data contains fluctuations due to the radiosonde swinging. These fluctuations should be filtered out. In existing UASSs, coordinate data are smoothed using piecewise linear approximation, approximation by small order polynomials or splines, and also a linear lowpass filter (LPF) [2,4]. This increases the averaging layer thickness and adds the dynamic error in the turbulence zones.
Unfortunately, there is no reference measuring system with which it would be possible to determine the measurement errors of atmospheric parameters, including wind. Currently, navigation UASSs are the most accurate means of measuring atmospheric parameters in a given altitude range, so an estimate of their errors can be obtained by indirect methods and using modeling. For example, the World Meteorological Organization (WMO) uses comparative tests of different UASSs to evaluate sounding accuracy [5]. This paper presents the results of radiosonde motion modeling to assess the smoothing efficiency of the sonde pendulum oscillations. As a result, it is possible to estimate the error in the measurement of velocity introduced by radiosonde swinging and subsequent filtering. In addition, the model allows configuring filter parameters more optimally.

Radiosonde oscillation modeling
Simplified, the balloon-sonde system can be represented as a mathematical pendulum with a movable suspension point (see Fig. 1). We represent the sonde as a material point A, having mass m, forced to move around a circle (in this plane) of radius l (cord length) centered at point O. The equation for this system can be written as [6]: , (1) where g is the gravitational acceleration, x(t), y(t) are functions that determine the horizontal and vertical movement of the suspension point,  is the deviation angle of the cord from the vertical ( > 0, if the cord deviates to the right), t is the time.
To simplify the equation, we can assume that vertical acceleration is close to zero and can be neglected. Then the equation of the balloon-sonde system takes the following form: . (2) Although Eq. (2) does not take into account the damping of the radiosonde oscillations, the fluctuation of the suspension point, and not the idealness of the pendulum, it allows us to study the system behavior and to study the possibilities of filtering these oscillations. The solution of this equation for a given function x(t) allows us to simulate the radiosonde movement, estimate the error of navigation data filtering by comparing the given function x(t) or with its estimate at the filter output. Simulation also allows you to select the appropriate filter parameters. Adaptive filtering can be applied to reduce dynamic error [7].
To obtain an analytical solution of the equation (2), it is necessary to linearize it near the pendulum's resting point. As a result, we can get a solution with an acceptable error for  in the range of ± 30 degrees. But the solutions presented in this paper were obtained by numerical methods. It is necessary to determine the function x(t) or for modeling. Let's consider several options. Let the horizontal balloon coordinate change as follows , where S is a function of the average motion of the radiosonde balloon, A is the amplitude of the oscillating movements component, ω is the angular frequency of the oscillating movements.
Using the sine function, a smooth change in the position and speed of the ball is modeled. Figure 2 shows fragments of wind profiles obtained by sounding of the Polyus UASS. Obviously, the wind field is not uniform. Wind velocity varies with altitude, sometimes increasing, sometimes decreasing relative to some average profile. The thickness of the layer in which the velocity changes first to one side, then to the other, may vary from 100 meters to several kilometers. Typically, the sonde velocity is in the range of 4 to 9 m s -1 . Therefore, for modeling a dynamic air environment, ω can be selected in the range from 2π⁄25 to 2π⁄11 rad s -1 . In calmer sections, ω is between 2π⁄400 to 2π⁄60 rad s -1 .  The balloon velocity along one of the horizontal coordinates is obtained by differentiating the expression (3) .
(4) At altitudes of up to 40 km, wind velocity usually does not exceed 50 m s -1 . The velocity limit is considered to be 100 m s -1 . Analysis of the Polyus upper-air sounding data shows that the modulus of the average vertical wind gradient can reach 0.04 m s -1 per 1 m of height. Vertical wind shear (i.e. increased velocity gradient) of high-altitude jet stream can reach 10 m s -1 per 1 km [8]. The International Civil Aviation Organization (ICAO) recommends the following gradation of vertical wind shear for safety analysis of aircraft take-off and landing, i.e. for the lower troposphere [9]: weak (0..2 m s -1 per 30 m), moderate (2..4 m s -1 per 30 m), strong (4..6 m s -1 per 30 m) and very strong (more than 6 m s -1 per 30 m). Table 1 shows the data on the wind distribution in the surface layer of the atmosphere with a height of 300 m in the vicinity of the Baikonur Cosmodrome [10] with abnormally strong winds at its upper boundary. The table shows that for this territory, the most probable gradients with strong winds at an altitude of 300 m lie in the range from 4 to 6 m s -1 per 100 m, which ICAO gradations correspond to a weak wind shear, and the maximum observed vertical gradient of 10.7 belongs to the category of moderate wind shear. Thus, the following norms of vertical gradients can be adopted for modeling: -the average gradient of the extended sections of the profile is in the range from 0.5 to 5 m s -1 per 100 m, -gradients of cyclically or short-term wind changing is in the range from 5 to 20 m s -1 per 100 m or more. It must be taken into account that in equation (4) the environment dynamics is determined by both the frequency ω and the amplitude Aω. The most difficult for measurements will be the environment with the maximum possible wind velocity gradient (modulo). In order to fix the vertical velocity gradient changing according to a sinusoidal law, it is necessary to set the amplitude Aω inversely proportional to the frequency ω. To do this, replace A in equation (4) by and obtain the following .
Let the average velocity of the radiosonde ascent is 5 m s -1 . Then, based on the norms adopted above for vertical velocity gradients, parameter B can be associated with a wind velocity gradient at 100 m of height G 100 as follows .
Thus, parameter B can be selected in the range from 0.4 to 2. Expression (5) allows us to fix the velocity gradient between adjacent sine extrema points with changing frequency ω. In a special case, you can increase the parameter B to 2.5 (wind shear of more than 30 m s -1 per 100 m).
S(t) can be chosen to increase linearly, then the velocity , and . Or quadratic S(t)=kt 2 , so that the average balloon velocity increases linearly as it ascents, S(t)=kt , and . For k = 0, the motion model with a quadratic dependence S(t) becomes equivalent to a linearly increasing one. If we limit the average increase in speed over the rise time to 60 m s -1 , then 2k = 60/8000 = 0.0075 m s -2 . This corresponds to a gradient G 100 = 0.15. In a real atmosphere, the average wind velocity can increase, decrease or remain without significant changes at different heights. Therefore, the gradients are usually greater. Above, we restricted our consideration of the average velocity gradient G 100 in the range from 0.5 to 5. Based on this, the coefficient k can be set in the range from 0.0125 up to 0.125. It should be taken into account that under normal conditions, it rarely exceeds 0.03.
The effect of the abrupt wind velocity change on the balloon-sonde system was also evaluated in the study. The abrupt change was set using the Gauss function: where a is the amplitude of the peak, b is the position of the peak center, and c determines the width of the bell-shaped curve. As a result, for the abrupt change in wind velocity, we obtain the following . (8) As noted above, piecewise-linear approximation, approximation by small-order polynomials or splines, as well as linear LPF are used to filter coordinates and velocity of radiosonde. Piecewise linear approximation is qualitatively inferior to other methods. The knees at the junction points of the linear approximation segments do not correspond to the physical properties of the air. Modern UASSs use computers with good computing power to filter coordinate-telemetric information. The receipt speed of the sonde data is relatively small. Therefore, the use of linear approximation is not justified, and we will not consider it.

Radiosonde velocity smoothing
Matlab and LabView approximate (smooth) a set of observations (X, Y) by a cubic spline by minimizing the following function (9) In our case, the set of observational data (X, Y) consists of the values of the measurement time X and the values of the velocity measurement or coordinate Y at the corresponding time instants. Next, p is the balance coefficient, is the weight coefficient of the i-th element, is the i-th element of Y, is the i-th element of X. is the second derivative of the cubic spline function f(x). λ(x) is a piecewise linear constant smoothness function , для i = 0, 1, …, n-2 (10) where is the i-th element of smoothness.
If you do not want to vary the degree of smoothness in different parts of the observations, then all the samples of the function can be equal to 1. Weighting factors by default are also equal to 1.
The balance coefficient p can be set in the range from 0 to 1 setting the desired balance between smoothness and a minimum standard deviation from the original data. If p = 0, then the cubic spline approximation will be equivalent to linear. If p = 1, then there is no smoothing and the approximation coincides with the original data. The closer p is to zero, the greater the smoothing. And vice versa, the closer p is to 1, the closer the approximation to the original function. Fig. 3 shows examples of smoothing by a cubic spline with different coefficients p. Consider smoothing by a cubic spline, when the horizontal velocity of the radiosonde balloon changes according to the expression (5) with linearly increasing and quadratic S(t). The maximum distinguishable frequency ω is limited by the frequency of the radiosonde oscillations on the cord. The frequency of free radiosonde oscillations should be at least 2 times higher than the ω. The minimum cord length is 10 meters, as it is believed that at a smaller distance, the balloon will affect the radiosonde sensors. With a 10-meter cord, the period of free radiosonde oscillations is approximately 6.35 s. Therefore, the maximum frequency ω, which we can distinguish, is approximately 2π⁄13 rad s -1 , which corresponds to disturbance with a thickness of 100 m at a radiosonde ascent velocity about 7.5 m s -1 , or a thickness of 65 m at an ascent velocity about 5 m s -1 . Figures 4-5 show the results of modeling the radiosonde motion and smoothing its velocity with a cubic spline for , and . The step of the input sequence (radiosonde speed) is 1 s, which corresponds to the standard GNSS receiver maximum output rate. Fig. 4.a shows the velocity profiles of the balloon and the radiosonde obtained as a result of modeling, and the smoothing result as an estimate of the wind (balloon) velocity. In fig. 4.b. shows a graph of the error in estimating wind velocity, calculated as the difference between a given balloon velocity and a smoothed radiosonde velocity. Fig. 5.a shows the dependence of the phase of radiosonde oscillations on the flight time. The phase diagram of the radiosonde velocity is shown in Fig. 5  The graphs show that for the limiting frequency w, it is possible to obtain an estimate of the fine wind field structure with an error of at least 1.5 m s -1 without taking into account the error of the navigation receiver, modeling, and other errors introduced by factors not taken into account in the model. At the boundaries of the smoothing interval, typical With a further increase in the frequency ω, it becomes impossible to restore the wind profile by smoothing; the sonde motion phase diagram becomes more complicated (see Fig. 6-7).

For
, which corresponds to a layer thickness of 100 m at an ascent velocity of 5 m s -1 , choosing the optimal value of the balance parameter p, we can achieve a filtering error of less than 1 m s -1 . Unfortunately, the balance coefficient chosen for is not optimal for laminar streams with frequencies in the range from 2π⁄400 to 2π⁄60 rad s -1 . For example, if , then smoothing with the same coefficient p gives an error about ± 0.5 m s -1 . By reducing the balance coefficient by 10 times, us get an error within ± 0.05 m s -1 , not considering edge effects. Of course, with such p, the error for a layer 100 m thick ( ) increases up to ± 4 m s -1 at B = 2. Decreasing B to typical values in the range from 0.2 to 0.5 does not affect measurement error for slow processes, but proportionally reduces the error for high frequencies ω at the optimal balance coefficient p for slower processes. For example, if in the case considered above with and p, selected for slow processes, we reduce B to 0.5, then the error decreases to 1 m s -1 .  In general, the value of the optimal coefficient p depends on the input data frequency characteristic, the ratio of this characteristic to the natural radiosonde frequency and the relative sampling rate.
Modeling with a quadratic dependence S(t) expectedly showed that an additional linear increase in velocity at k in the range from 0 to 0.125 does not affect the error of the velocity estimate.
Smoothing the sonde velocity with a low-pass filter with a finite impulse response (FIR) for balloon motion according to expression (5) gives similar results to the cubic spline. Cubic spline smoothing gives a slightly smaller error for small ω. The following FIR filter parameters were selected: the cutoff frequency is 0.078 Hz, the sampling frequency is 1 Hz, and the filter length is 60.
Consider the case when the wind velocity does not change harmoniously, but abruptly (see equation 7). Parameters a and c should be set based on the norms adopted above for possible values of the vertical wind velocity gradient. For example, if we consider a jump 100 m wide in height, then with a sonde ascent speed of 5 m s -1 , you can set c = 20/6 = 3.33. The amplitude a in this case can be selected in the range from 10 to 40 m s -1 .
Smoothing with a cubic spline gives an error of more than 3 m s -1 at a = 10, p = 0.005 and c = 3.33. And at a = 40, the error modulus increases up to 15 m s -1 . An increase in the balance coefficient p allows us to reduce errors at the jump location to 4 m s -1 at a = 40, but at the same time the filtering error of the radiosonde pendulum oscillations in "quiet" sections increases. By increasing the weight coefficients w i (see formula 9) in the zone of the velocity jump, the error can be slightly reduced. The dynamic error of the FIR filter at a velocity jump location is less than that of a cubic spline. For a = 40 and c = 3.33, the error was ± 3 m s -1 , and for a = 10 less than 1 m s -1 .

Conclusions
Using computer simulation, we obtained estimates of wind velocity measurement errors introduced by radiosonde swinging on the balloon cord. The largest errors are expected in small thickness turbulence zones.
Smoothing with a cubic spline is preferable for laminar layers of the atmosphere and for obtaining average estimates. FIR filter smoothing is more suitable for small turbulent zones. Dynamic filtering errors can be reduced by adaptively changing the smoothing coefficients.
The maximum resolution of the wind field for navigational UASSs is limited by the natural frequency of the sonde oscillations, i.e. cord length. The shorter the cord, the better the resolution.