Four-dimensional mesospheric and lower thermospheric wind fields using Gaussian process regression on multistatic specular meteor radar observations

Mesoscale dynamics in the mesosphere and lower thermosphere (MLT) region have been difficult to study from either ground- or satellite-based observations. For understanding of atmospheric coupling processes, important spatial scales at these altitudes range between tens and hundreds of kilometers in the horizontal plane. To date, this scale size is challenging observationally, so structures are usually parameterized in global circulation models. The advent of multistatic specular meteor radar networks allows exploration of MLT mesoscale dynamics on these scales using an increased number of detections and a diversity of viewing angles inherent to multistatic networks. In this work, we introduce a four-dimensional wind field inversion method that makes use of Gaussian process regression (GPR), which is a nonparametric and Bayesian approach. The method takes measured projected wind velocities and prior distributions of the wind velocity as a function of space and time, specified by the user or estimated from the data, and produces posterior distributions for the wind velocity. Computation of the predictive posterior distribution is performed on sampled points of interest and is not necessarily regularly sampled. The main benefits of the GPR method include this non-gridded sampling, the built-in statistical uncertainty estimates, and the ability to horizontally resolve winds on relatively small scales. The performance of the GPR implementation has been evaluated on Monte Carlo simulations with known distributions using the same spatial and temporal sampling as 1 d of real meteor measurements. Based on the simulation results we find that the GPR implementation is robust, providing wind fields that are statistically unbiased with statistical variances that depend on the geometry and are proportional to the prior velocity variances. A conservative and fast approach can be straightforwardly implemented by employing overestimated prior variances and distances, while a more robust but computationally intensive approach can be implemented by employing training and fitting of model hyperparameters. The latter GPR approach has been applied to a 24 h dataset and shown to compare well to previously used homogeneous and gradient methods. Small-scale features have reasonably low statistical uncertainties, implying geophysical wind field horizontal structures as low as 20–50 km. We suggest that this GPR approach forms a suitable method for MLT regional and weather studies.

1 Introduction

The mesoscale neutral dynamics of the mesosphere and lower thermosphere (MLT) region are challenging to study, despite their importance in global circulation models. Due to the lack of observations, these scales are usually parameterized in models (e.g., Liu, 2019). MLT large-scale dynamics have been studied with monostatic specular meteor radars (SMRs) by providing mean horizontal winds over areas with an approximately 200–300 km radius at MLT altitudes and 1–2 h and 2–3 km temporal and vertical resolutions, respectively (e.g., Hocking et al., 2001; Holdsworth et al., 2004). These measurements have made significant contributions to community understanding of the climatological behavior of mean winds, planetary waves, and total tides over a variety of SMR monostatic sites (e.g., Mitchell et al., 1999, 2002; Pancheva et al., 2002; Sandford et al., 2006; Hoffmann et al., 2010). Moreover, when the winds from more than one SMR widely separated in longitude at a similar latitude are combined, spatiotemporal ambiguities of tides and planetary waves have been successfully resolved (e.g., Murphy, 2003; Murphy et al., 2006; He et al., 2018; He and Chau, 2019). Monostatic SMRs have also been used to study MLT gravity wave momentum flux with wide and narrow beam observing configurations, with the caveat that spatial and temporal contributions are combined (e.g., Hocking, 2005; Fritts et al., 2012; Andrioli et al., 2013; Placke et al., 2015).

Recently, multistatic configurations have been proposed to complement these previous studies and to allow the investigation of MLT mesoscale dynamics. These configurations include the MMARIA (Multistatic Multi-frequency Agile Radar Investigations of the Atmosphere) concept (Stober and Chau, 2015; Chau et al., 2017). This concept has been further augmented by the SIMONe (Spread Spectrum Interferometric Multistatic meteor radar Observing Network) approach (Chau et al., 2019). By using recent technological developments in atmospheric radars, such as spread-spectrum, MIMO (multi-input, multiple-output), and compressed sensing approaches (Vierinen et al., 2016; Urco et al., 2018, 2019b), SIMONe allows the implementation of MMARIA with several attractive qualities: it is easier, cheaper, and inherently expandable compared to original proposed configurations using traditional pulsed systems. Examples of SIMONe implementations in Germany, Peru, and Argentina can be found in several studies (Vierinen et al., 2019; Charuvil Asokan et al., 2020; Vargas et al., 2021; Chau et al., 2021; Conte et al., 2021).

Multistatic observing approaches allow a large increase in scattering detections per unit time along with observation of the same volume from different viewing angles. These two features unlock the possibility of estimating the spatial features of the wind within the observed volume. Depending on the resolutions and spatial scales covered, different aspects of MLT mesoscale dynamics and coupling can be studied with the technique. For example, at scales between a few tens of kilometers to a few hundreds of kilometers, the contributions of gravity waves and strongly stratified turbulence can be studied with multistatic approaches (e.g., Roberts and Larsen, 2014; Marino et al., 2015).

The spatial structure of horizontal winds has also been pursued using a variety of other techniques including meteorological radars in the lower atmosphere, coherent scatter radars in the mesosphere, and Fabry–Pérot interferometers in the thermosphere (e.g., Browning and Wexler, 1968; Chau et al., 2017; Meriwether et al., 2008). As in the case of the initial MMARIA analysis, these techniques typically approximate wind fields as analytic, differentiable polynomials in order to obtain gradients of the horizontal winds. Although they provide additional spatial information beyond direct single-point information, these methods can aggressively smooth real spatial structure and, in some cases, can introduce artificial structure, particularly in regions with sparse or noisy measurements. In recent years, a variety of analysis approaches using statistical inverse theory have been applied to these and similar problems. These studies have the goal of estimating the spatial structure of multi-point projected wind velocities and electric fields (e.g., Nicolls et al., 2014; Hysell et al., 2014; Harding et al., 2015; Stober et al., 2018). For example, a Tikhonov regularization originally developed for a optical network of Fabry–Pérot interferometers (Harding et al., 2015) has been adapted to yield MLT wind fields over Peru (Chau et al., 2021).

As in any statistical inverse theory problem, more independent samples are desirable to reduce the impact of regularization constraints and to improve the quality of the estimates. In November 2018, a short observing campaign was conducted in northern Germany, herein denoted SIMONe2018, in which six existing MMARIA links were complemented with eight additional SIMONe links. During this campaign, we obtained on average 200 000 meteor scatter observations per day (e.g., Vierinen et al., 2019; Charuvil Asokan et al., 2020). For reference, a monostatic SMR obtains on average 10 000 meteors per day at a comparable latitude and seasonal time.

Some previous analysis methods have been published on multistatic observations of MLT mesoscale dynamics, such as the gradient method and variants of Tikhonov regularization (Chau et al., 2017; Stober et al., 2018; Chau et al., 2021). However, given the novelty of multistatic measurements and the lack of a reliable ground-truth observation, different wind field approaches still need to be explored, particularly in the properties of resulting statistical measures of bias and variance. In this work, we introduce a multistatic analysis technique based on Gaussian process regression (GPR) (Rasmussen and Williams, 2006). Some of the main benefits of GPR are that analysis predictions essentially interpolate the measurements (within error bounds) and that final output products inherently include quantitative uncertainties.

GPR is a Bayesian and nonparametric approach currently being used in many different machine-learning applications (e.g., Wahlström et al., 2013; Foreman-Mackey et al., 2017). As a Bayesian technique, a key user input is the specification of a prior distribution for the values to be estimated, including hyperparameters of the distribution. Despite needing these hyperparameters, GPR is nonparametric in the sense that it does not compress the training data into a finite-dimensional parameter vector, in contrast to parametric methods like linear regression (Rasmussen and Williams, 2006). GPR is also known in other fields as kriging, and it has a long history of use in geostatistics under that name (Matheron, 1973; Journel and Huijbregts, 1978; Daley, 1991). Deep connections can be found between GPR and interpolation techniques using reproducing kernel Hilbert spaces (Scheuerer et al., 2013), including those that employ regularization. This ties GPR mathematically to the previously mentioned wind field estimation techniques, but the Bayesian viewpoint afforded by GPR can be more natural for expressing prior information and analyzing uncertainty. We direct the reader to Rasmussen and Williams (2006) for a general discussion of GPR and its place in the wider estimation landscape.

In this article, we start by introducing the wind estimation problem and geometrical considerations. Next, we present the wind field estimation method using GPR, including the necessary mathematical expressions. The proposed estimation is subsequently applied to both Monte Carlo simulations and to measurements from the SIMONe2018 campaign in Sects. 5 and 6, respectively. In the latter section, estimated wind fields are compared to the winds obtained with the homogeneous and gradient methods, i.e., to the zero- and first-order Taylor expansion approximations. Finally, we discuss the benefits and challenges of the proposed estimation approach for MLT wind field studies.

2 Specular meteor radar measurements and geometry

SMRs receive echoes from meteor trails when the radar Bragg vector (k_B) points perpendicular to them. The Doppler shift (f) of the received signal of a meteor echo at time t and location given by longitude, latitude, and altitude $(\mathrm{\Lambda },\mathrm{\Phi },z)$ results from the projection of the atmospheric wind vector (u) in the Bragg vector k_B (e.g., Hocking et al., 2001; Holdsworth et al., 2004), i.e.,

$$\begin{array}{}\text{(1)}& f(\mathrm{\Lambda },\mathrm{\Phi },z,t)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\left[\begin{array}{ccc}{k}_u& k_v& k_z\end{array}\right]\left[\begin{array}{c}u(\mathrm{\Lambda },\mathrm{\Phi },z,t)\\ v(\mathrm{\Lambda },\mathrm{\Phi },z,t)\\ w(\mathrm{\Lambda },\mathrm{\Phi },z,t)\end{array}\right],\end{array}$$

where k_u, k_v, and k_w are the Bragg vector components of k_B, and u, v, and w are wind vector components of u in the zonal (east), meridional (north), and vertical (up) directions, respectively. The Bragg vector is given by the difference of the scattered and incident wave vectors, i.e., ${\mathit{k}}_{\mathrm{B}}={\mathit{k}}_{\mathrm{s}}-{\mathit{k}}_{\mathrm{i}}$. Using interferometry on reception, the angle of arrival (AOA) is obtained. In the case of MIMO systems, interferometry is also implemented on transmission, allowing measurement of the angle of departure (AOD) (e.g., Chau et al., 2019). By combining these angles along with range information, the meteor location $(\mathrm{\Lambda },\mathrm{\Phi },z)$ and Bragg wave vectors are obtained. In the reductive case of monostatic systems, ${\mathit{k}}_{\mathrm{B}}=-\mathrm{2}{\mathit{k}}_{\mathrm{i}}$ and its magnitude is equal to $\mathrm{4}\mathit{\pi }/\mathit{\lambda }$, where λ is the radar wavelength.

As mentioned above, traditionally a mean horizontal wind has been obtained from analysis that simultaneously solves N equations of the form of Eq. (1), with the assumption that the wind is constant in the observed volume (zero-order Taylor approximation or homogeneous method). The data for the N equation set were obtained by binning desired observations with regular altitude and temporal resolutions. In general, with a sufficient number of meteors and viewing angles, the method yields spatial information on the wind inside the observed volume. For example, Chau et al. (2017) implemented a gradient method, whereby the wind field estimation includes the first-order Taylor expansion terms.

In multistatic geometries, both the observed volumes and separations of the multistatic links are relatively large. For this reason, it is necessary to take the Earth's geoid shape into account. Moreover, the GPR model described in the next section is directly dependent on calculating coordinate distances accurately. This implies that altitudes and horizontal distances that account for the Earth's curvature, which is the measurement goal, must also try to minimize mapping distortions, particularly in distance scaling. Use of a naive geometric projection such as the equirectangular projection, in which latitude and longitude are simply scaled to yield x–y coordinates in meters, does not satisfy these requirements. Therefore, in this work, we use a local azimuthal equidistant projection centered in the observing region, with Earth shape based on the well-known WGS84 geoid model. This projection is used to transform longitude and latitude into local x and y coordinates, where horizontal distance in x and y reasonably approximates the true geodesic distance. Subsequently, we use these (x,y) projected coordinates in place of (Λ,Φ) geodetic coordinates from Eq. (1). Note that this does not change the definitions of $(u,v,w)$ and k_B, which remain aligned with a local east–north–up coordinate system and not, in general, with the projected x and y coordinates.

To represent a set of Doppler wind measurements, we use the following notation for the measurement equation. Let ${\mathit{x}}_m=(t_m,z_m,y_m,x_m)$ denote the coordinates for a measurement m of M. Then the ensemble of coordinates is given by the matrix X as

$$\begin{array}{}\text{(2)}& \mathbf{X}=\left[\begin{array}{c}{\mathit{x}}_{\mathrm{1}}^{\mathit{⊺}}\\ \mathrm{⋮}\\ {\mathit{x}}_M^{\mathit{⊺}}\end{array}\right]=\left[\begin{array}{cccc}{t}_{\mathrm{1}}& z_{\mathrm{1}}& y_{\mathrm{1}}& x_{\mathrm{1}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}& \mathrm{⋮}\\ t_M& z_M& y_M& x_M\end{array}\right],\end{array}$$

and the corresponding wind vectors are given by

$$\begin{array}{}\text{(3)}& \mathit{u}=\left[\begin{array}{c}u\left({\mathit{x}}_{\mathrm{1}}\right)\\ \mathrm{⋮}\\ u\left({\mathit{x}}_M\right)\end{array}\right]\mathit{v}=\left[\begin{array}{c}v\left({\mathit{x}}_{\mathrm{1}}\right)\\ \mathrm{⋮}\\ v\left({\mathit{x}}_M\right)\end{array}\right]\mathit{w}=\left[\begin{array}{c}w\left({\mathit{x}}_{\mathrm{1}}\right)\\ \mathrm{⋮}\\ w\left({\mathit{x}}_M\right)\end{array}\right].\end{array}$$

We group the Bragg vectors of a set of measurements by component and combine with the $\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}$ scaling to give u, v, and w measurement vectors:

$$\begin{array}{}\text{(4)}& {\mathit{a}}_u={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\left[\begin{array}{c}{k}_{u_{\mathrm{1}}}\\ \mathrm{⋮}\\ k_{u_M}\end{array}\right]{\mathit{a}}_v={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\left[\begin{array}{c}{k}_{v_{\mathrm{1}}}\\ \mathrm{⋮}\\ k_{v_M}\end{array}\right]{\mathit{a}}_w={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\left[\begin{array}{c}{k}_{w_{\mathrm{1}}}\\ \mathrm{⋮}\\ k_{w_M}\end{array}\right].\end{array}$$

Finally, using ⊙ to denote the element-wise (Hadamard) vector product, our measurement equation following from Eq. (1) for the ensemble of Doppler measurements f is

$$\begin{array}{}\text{(5)}& \mathit{f}\left(\mathbf{X}\right)=\left[\begin{array}{c}f\left({\mathit{x}}_{\mathrm{1}}\right)\\ \mathrm{⋮}\\ f\left({\mathit{x}}_M\right)\end{array}\right]={\mathit{a}}_u\odot \mathit{u}+{\mathit{a}}_v\odot \mathit{v}+{\mathit{a}}_w\odot \mathit{w}+\mathit{ϵ},\end{array}$$

where $\mathit{ϵ}\sim \mathcal{N}(\mathrm{0},{\mathbf{\Sigma }}_n)$ is zero-mean Gaussian measurement uncertainty with covariance Σ_n.

3 Estimation problem

The estimation task is to take a set of Doppler measurements f and infer wind values u(x^′), v(x^′), w(x^′) at a chosen location x^′ using the measurement model from Eq. (5). We employ Gaussian process regression (GPR) to model the winds and hence Doppler measurements as a stochastic process. This approach allows estimation at arbitrary coordinates (convenient for random meteor locations and non-gridded prediction) and produces statistical uncertainty as an output product.

Our GPR method is implemented as a three-stage process. First, one defines the form for the model, which includes mean and covariance functions and their hyperparameters. Then, one fully specifies the model by setting hyperparameter values, either through prior knowledge or a separate fitting process. Finally, one applies the specified model to a set of measurements to calculate the posterior predictive distribution and make an estimate at points of interest. Figure 1 summarizes our implementation in a block diagram. In the following paragraphs, we describe the method in detail.

Figure 1 Block diagram of processing flow. The blocks in orange indicate input from the user; blocks in green belong to the GPR model, and the estimates are obtained in the red block (see text for details).

3.1 Gaussian process definitions

For a function f(x) drawn from a Gaussian process, we write

$$\begin{array}{}\text{(6)}& f\left(\mathit{x}\right)\sim \mathcal{G}\mathcal{P}\left(m\right(\mathit{x}),\mathit{\kappa }(\mathit{x},{\mathit{x}}^{\mathbf{\prime }}\left)\right).\end{array}$$

This representation is fully defined by mean and covariance functions, which describe the first- and second-order statistics:

$$\begin{array}{}\text{(7)}& {\displaystyle }m\left(\mathit{x}\right)& {\displaystyle }=\mathbf{E}\left[f\right(\mathit{x}\left)\right],\text{(8)}& {\displaystyle }\mathit{\kappa }(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})& {\displaystyle }=\mathbf{E}\left[\right(f\left(\mathit{x}\right)-m\left(\mathit{x}\right)\left)\right(f\left({\mathit{x}}^{\mathbf{\prime }}\right)-m\left({\mathit{x}}^{\mathbf{\prime }}\right)\left)\right],\end{array}$$

where E denotes the expected value. Gaussian processes are convenient because evaluating them at a set of points leads to a Gaussian random vector:

$$\begin{array}{}\text{(9)}& \left[\begin{array}{c}f\left({\mathit{x}}_{\mathrm{1}}\right)\\ \mathrm{⋮}\\ f\left({\mathit{x}}_N\right)\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}m\left({\mathit{x}}_{\mathrm{1}}\right)\\ \mathrm{⋮}\\ m\left({\mathit{x}}_N\right)\end{array}\right],\left[\begin{array}{ccc}\mathit{\kappa }({\mathit{x}}_{\mathrm{1}},{\mathit{x}}_{\mathrm{1}})& \mathrm{\cdots }& \mathit{\kappa }({\mathit{x}}_{\mathrm{1}},{\mathit{x}}_N)\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ \mathit{\kappa }({\mathit{x}}_N,{\mathit{x}}_{\mathrm{1}})& \mathrm{\cdots }& \mathit{\kappa }({\mathit{x}}_N,{\mathit{x}}_N)\end{array}\right]\right),\end{array}$$

which enables tractable computation. We recast this compactly using matrix notation as

$$\begin{array}{}\text{(10)}& \mathit{f}\left(\mathbf{X}\right)\sim \mathcal{N}\left(\mathit{m}\right(\mathbf{X}),\mathbf{K}(\mathbf{X},\mathbf{X}\left)\right).\end{array}$$

It might seem like this model is too simple to be useful, but Gaussian processes have a lot of flexibility to fit a wide variety of functions because the posterior distribution is constructed nonparametrically and directly incorporates the measurements. Additionally, a modeler has a lot of freedom in applying Gaussian processes by choosing the form of the mean and covariance functions, including specifying hyperparameters.

3.2 Wind component prior distributions

Since we want to estimate the wind components, we model them as independent Gaussian processes.

$$\begin{array}{}\text{(11)}& {\displaystyle }u\left(\mathit{x}\right)& {\displaystyle }\sim \mathcal{G}\mathcal{P}\left(m_u\right(\mathit{x}),{\mathit{\kappa }}_u(\mathit{x},{\mathit{x}}^{\mathbf{\prime }}\left)\right)\text{(12)}& {\displaystyle }v\left(\mathit{x}\right)& {\displaystyle }\sim \mathcal{G}\mathcal{P}\left(m_v\right(\mathit{x}),{\mathit{\kappa }}_v(\mathit{x},{\mathit{x}}^{\mathbf{\prime }}\left)\right)\text{(13)}& {\displaystyle }w\left(\mathit{x}\right)& {\displaystyle }\sim \mathcal{G}\mathcal{P}\left(m_w\right(\mathit{x}),{\mathit{\kappa }}_w(\mathit{x},{\mathit{x}}^{\mathbf{\prime }}\left)\right)\end{array}$$

Assuming Gaussianity of the wind processes is not simply for convenience (although it does enable closed-form computation). Given some mean and covariance, a Gaussian distribution has the maximum entropy (Cover and Thomas, 2006). In other words, assuming normality imposes the minimal prior information about the wind processes within a second-order statistical framework. The winds likely have more structure than this, including cross-covariances between the components, but this assumption ensures conservative estimates without prior knowledge of the true statistical structure of the wind processes.

Many choices for the mean functions are possible, but for simplicity we restrict our attention to means that are fixed without tunable hyperparameters. Even under this restriction, one can use a standard parametric model for the mean functions, and as long as the parameter-fitting is done with linear regression prior to GPR analysis, no additional hyperparameters are added to the GPR model. In general, the mean functions have less impact on the GPR results than the covariance functions, and we will see later how the posterior predictive distribution is more strongly driven by the measurements and the covariance functions. Often a zero mean is sufficient to produce good results (Rasmussen and Williams, 2006), and that holds for this case as well. Nevertheless, the mean can be useful for including well-known effects. In the models for subsequent sections, we have used two cubic splines taken as a tensor product over altitude and time to produce a mean that accounts for large-scale tidal components.

For the covariance functions, we choose a functional form for which each wind component has an independent amplitude multiplying a common distance kernel.

$$\begin{array}{}\text{(14)}& {\displaystyle }{\mathit{\kappa }}_u(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})& {\displaystyle }={\mathit{\sigma }}_u^{\mathrm{2}}{\mathit{\kappa }}_{\mathrm{d}}(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})\text{(15)}& {\displaystyle }{\mathit{\kappa }}_v(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})& {\displaystyle }={\mathit{\sigma }}_v^{\mathrm{2}}{\mathit{\kappa }}_{\mathrm{d}}(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})\text{(16)}& {\displaystyle }{\mathit{\kappa }}_w(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})& {\displaystyle }={\mathit{\sigma }}_w^{\mathrm{2}}{\mathit{\kappa }}_{\mathrm{d}}(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})\end{array}$$

Using a common distance kernel is convenient for simplifying computations, and we expect that relaxing this assumption in the future would allow for increased expressiveness at the cost of computational burden. The distance kernel κ_d is chosen to be the Matérn covariance with $\mathit{\nu }=\frac{\mathrm{5}}{\mathrm{2}}$ using length scales given by δ_t, δ_z, δ_y, and δ_x for the coordinate dimensions:

$$\begin{array}{}\text{(17)}& {\mathit{\kappa }}_{\mathrm{d}}(\mathit{x},{\mathit{x}}^{\mathbf{\prime }})=\left(\mathrm{1}+\sqrt{\mathrm{5}}r+{\displaystyle \frac{\mathrm{5}}{\mathrm{3}}}{r}^{\mathrm{2}}\right)e^{-\sqrt{\mathrm{5}}r},\end{array}$$

with

$$\begin{array}{}\text{(18)}& {\displaystyle }r& {\displaystyle }={∥{\displaystyle \frac{\mathit{x}-{\mathit{x}}^{\mathbf{\prime }}}{\mathit{\delta }}}∥}_{\mathrm{2}},\text{(19)}& {\displaystyle }\mathit{\delta }& {\displaystyle }={\left[\begin{array}{cccc}{\mathit{\delta }}_t& {\mathit{\delta }}_z& {\mathit{\delta }}_y& {\mathit{\delta }}_x\end{array}\right]}^{\mathit{⊺}},\end{array}$$

where ∥.∥₂ represents the Euclidean norm. Altogether, this results in a hyperparameter set θ of

$$\begin{array}{}\text{(20)}& \mathit{\theta }={\left[\begin{array}{ccccccc}{\mathit{\sigma }}_u^{\mathrm{2}}& {\mathit{\sigma }}_v^{\mathrm{2}}& {\mathit{\sigma }}_w^{\mathrm{2}}& {\mathit{\delta }}_t& {\mathit{\delta }}_z& {\mathit{\delta }}_y& {\mathit{\delta }}_x\end{array}\right]}^{\mathit{⊺}}\end{array}$$

for the GPR wind model. We chose the Matérn-$\frac{\mathrm{5}}{\mathrm{2}}$ covariance because it is twice-differentiable but not infinitely differentiable, so it provides relatively smooth functions while still allowing for rapid, geophysically driven changes that might be expected in wind fields. It is a typical choice for physical processes for this reason across a wide series of applications (Rasmussen and Williams, 2006).

Jointly and in matrix notation, we then write the Gaussian random vectors for the winds at a set of points X as

$$\begin{array}{}\text{(21)}& \left[\begin{array}{c}\mathit{u}\\ \mathit{v}\\ \mathit{w}\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}{\mathit{m}}_u\left(\mathbf{X}\right)\\ {\mathit{m}}_v\left(\mathbf{X}\right)\\ {\mathit{m}}_w\left(\mathbf{X}\right)\end{array}\right],\left[\begin{array}{ccc}{\mathbf{K}}_u(\mathbf{X},\mathbf{X})& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{K}}_v(\mathbf{X},\mathbf{X})& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{K}}_w(\mathbf{X},\mathbf{X})\end{array}\right]\right).\end{array}$$

Note that since we have defined the wind component processes independently, the cross terms are zero in the joint covariance matrix. However, this is not to say that we strictly enforce zero cross-covariance between the wind terms with this model. Rather, it is more accurate to say that we do not require prior knowledge of the cross-covariance but also cannot benefit from the improved estimation that such knowledge would provide.

3.3 Doppler measurement prior distribution

Since we are taking the wind components as Gaussian processes and Eq. (5) provides a linear relationship between the wind components and Doppler measurements, the Doppler measurements themselves also take the form of a Gaussian process. For a set of measurements f corresponding to the locations X, this produces a formulation as

$$\begin{array}{}\text{(22)}& \mathit{f}\sim \mathcal{N}\left({\mathit{m}}_f\right(\mathbf{X}),{\mathbf{K}}_f(\mathbf{X},\mathbf{X}\left)\right),\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }{\mathit{m}}_f\left(\mathbf{X}\right)& {\displaystyle }={\mathit{a}}_u\odot {\mathit{m}}_u\left(\mathbf{X}\right)+{\mathit{a}}_v\odot {\mathit{m}}_v\left(\mathbf{X}\right)+{\mathit{a}}_w\odot {\mathit{m}}_w\left(\mathbf{X}\right)\\ {\displaystyle }{\mathbf{K}}_f(\mathbf{X},\mathbf{X})& {\displaystyle }=\left({\mathit{a}}_u{{\mathit{a}}_u}^{\mathit{⊺}}\right)\odot {\mathbf{K}}_u(\mathbf{X},\mathbf{X})+\left({\mathit{a}}_v{{\mathit{a}}_v}^{\mathit{⊺}}\right)\odot {\mathbf{K}}_v(\mathbf{X},\mathbf{X})\\ {\displaystyle }& {\displaystyle }+\left({\mathit{a}}_w{{\mathit{a}}_w}^{\mathit{⊺}}\right)\odot {\mathbf{K}}_w(\mathbf{X},\mathbf{X})+{\mathbf{\Sigma }}_n.\end{array}$$

Note that the Gaussian process being measured is a linear composition. This is only a minor concern for our application, but it does make the formulation slightly different from the more typical examples. The following subsections provide the explicit formulas necessary to perform hyperparameter fitting and wind estimation using this model.

3.4 Model hyperparameter fitting

Fitting for the model hyperparameters θ involves maximizing the likelihood function for the marginal distribution pertaining to a set of measurements. Assuming Doppler measurements f coming from the distribution defined in Eq. (22), the negative log-likelihood as a function of the hyperparameters is

$$\begin{array}{}\text{(23)}& -l\left(\mathit{\theta }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{f}-{\mathit{m}}_f{)}^{\mathit{⊺}}{{\mathbf{K}}_f}^{-\mathrm{1}}(\mathit{f}-{\mathit{m}}_f)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\log \det {\mathbf{K}}_f-C,\end{array}$$

where C is a fixed scaling constant. Minimizing this function requires evaluating the gradient of the negative log-likelihood. For each hyperparameter θ_i, we thus have

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial (-l(\mathit{\theta }\left)\right)}{\partial {\mathit{\theta }}_{\mathrm{i}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{Tr}\left((\mathit{\alpha }{\mathit{\alpha }}^{\mathit{⊺}}-{{\mathbf{K}}_f}^{-\mathrm{1}}){\displaystyle \frac{\partial {\mathbf{K}}_f}{\partial {\mathit{\theta }}_{\mathrm{i}}}}\right),\end{array}$$

where

$$\begin{array}{}\text{(25)}& \mathit{\alpha }={{\mathbf{K}}_f}^{-\mathrm{1}}(\mathit{f}-{\mathit{m}}_f).\end{array}$$

Continuing down the derivative chain for each type of hyperparameter produces

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\partial {\mathbf{K}}_f}{\partial {\mathit{\sigma }}_{\mathrm{i}}^{\mathrm{2}}}}& {\displaystyle }=\left({\mathit{a}}_{\mathrm{i}}{{\mathit{a}}_{\mathrm{i}}}^{\mathit{⊺}}\right)\odot {\mathbf{K}}_{\mathrm{d}},\text{(27)}& {\displaystyle \frac{\partial {\mathbf{K}}_f}{\partial {\mathit{\delta }}_{\mathrm{i}}}}& {\displaystyle }=\left({\mathit{\sigma }}_u^{\mathrm{2}}\left({\mathit{a}}_u{{\mathit{a}}_u}^{\mathit{⊺}}\right)+{\mathit{\sigma }}_v^{\mathrm{2}}\left({\mathit{a}}_v{{\mathit{a}}_v}^{\mathit{⊺}}\right)+{\mathit{\sigma }}_w^{\mathrm{2}}\left({\mathit{a}}_w{{\mathit{a}}_w}^{\mathit{⊺}}\right)\right)\odot {\displaystyle \frac{\partial {\mathbf{K}}_{\mathrm{d}}}{\partial {\mathit{\delta }}_{\mathrm{i}}}},\end{array}$$

and

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{\partial {\mathit{\kappa }}_{\mathrm{d}}({\mathit{x}}_j,{\mathit{x}}_k)}{\partial {\mathit{\delta }}_{\mathrm{i}}}}={\displaystyle \frac{\mathrm{5}}{\mathrm{3}}}(\mathrm{1}+\sqrt{\mathrm{5}}r)e^{-\sqrt{\mathrm{5}}r}{\displaystyle \frac{\mathrm{1}}{{\mathit{\delta }}_{\mathrm{i}}}}{\left({\displaystyle \frac{({\mathit{x}}_j{)}_{\mathrm{i}}-({\mathit{x}}_k{)}_{\mathrm{i}}}{{\mathit{\delta }}_{\mathrm{i}}}}\right)}^{\mathrm{2}},\end{array}$$

where

$$\begin{array}{}\text{(29)}& r={∥{\displaystyle \frac{{\mathit{x}}_j-{\mathit{x}}_k}{\mathit{\delta }}}∥}_{\mathrm{2}}.\end{array}$$

With the objective and gradient known, fitting for θ then involves feeding these functions into an appropriate optimization routine. We have observed the most reliable convergence using SciPy's implementation of the L-BFGS-B and SLSQP algorithms (Virtanen et al., 2020).

3.5 Wind estimation

Having defined the model hyperparameters either through fitting or prior specification, estimating the winds at a set of prediction points X_* involves evaluating the posterior probability distribution given the measurements.

We start with the joint distribution between the measurements and the winds at the prediction points, which from previous definitions is given by

$$\begin{array}{}\text{(30)}& \left[\begin{array}{c}\mathit{f}\\ {\mathit{u}}_{*}\\ {\mathit{v}}_{*}\\ {\mathit{w}}_{*}\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}{\mathit{m}}_f\left(\mathbf{X}\right)\\ {\mathit{m}}_u\left({\mathbf{X}}_{*}\right)\\ {\mathit{m}}_v\left({\mathbf{X}}_{*}\right)\\ {\mathit{m}}_w\left({\mathbf{X}}_{*}\right)\end{array}\right],{\mathbf{K}}_{\mathrm{tot}}\right),\end{array}$$

where

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{\mathbf{K}}_{\mathrm{tot}}=& \left[\begin{array}{cc}{\mathbf{K}}_f(\mathbf{X},\mathbf{X})& {\mathit{a}}_u\odot {\mathbf{K}}_u(\mathbf{X},{\mathbf{X}}_{*})\\ {\mathbf{K}}_u({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_u& {\mathbf{K}}_u({\mathbf{X}}_{*},{\mathbf{X}}_{*})\\ {\mathbf{K}}_v({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_v& \mathbf{0}\\ {\mathbf{K}}_w({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_w& \mathbf{0}\end{array}\right]\\ & \left[\begin{array}{cc}{\mathit{a}}_v\odot {\mathbf{K}}_v(\mathbf{X},{\mathbf{X}}_{*})& {\mathit{a}}_w\odot {\mathbf{K}}_w(\mathbf{X},{\mathbf{X}}_{*})\\ \mathbf{0}& \mathbf{0}\\ {\mathbf{K}}_v({\mathbf{X}}_{*},{\mathbf{X}}_{*})& \mathbf{0}\\ \mathbf{0}& {\mathbf{K}}_w({\mathbf{X}}_{*},{\mathbf{X}}_{*})\end{array}\right].\end{array}\end{array}$$

The posterior predictive distribution follows from conditioning on the measurements:

$$\begin{array}{}\text{(32)}& \left[\begin{array}{c}{\mathit{u}}_{*}\\ {\mathit{v}}_{*}\\ {\mathit{w}}_{*}\end{array}\right]|\mathit{f}\sim \mathcal{N}({\mathit{m}}_{\mathrm{post}},{\mathbf{K}}_{\mathrm{post}}),\end{array}$$

where

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{\mathit{m}}_{\mathrm{post}}=\left[\begin{array}{c}{\mathit{m}}_u\left({\mathbf{X}}_{*}\right)\\ {\mathit{m}}_v\left({\mathbf{X}}_{*}\right)\\ {\mathit{m}}_w\left({\mathbf{X}}_{*}\right)\end{array}\right]& +\left[\begin{array}{c}{\mathbf{K}}_u({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_u\\ {\mathbf{K}}_v({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_v\\ {\mathbf{K}}_w({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_w\end{array}\right]\\ & {\mathbf{K}}_f(\mathbf{X},\mathbf{X}{)}^{-\mathrm{1}}(\mathit{f}-{\mathit{m}}_f\left(\mathbf{X}\right))\end{array}\end{array}$$

and

$$\begin{array}{}\text{(34)}& \begin{array}{ll}{\mathbf{K}}_{\mathrm{post}}=& \left[\begin{array}{ccc}{\mathbf{K}}_u({\mathbf{X}}_{*},{\mathbf{X}}_{*})& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{K}}_v({\mathbf{X}}_{*},{\mathbf{X}}_{*})& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{K}}_w({\mathbf{X}}_{*},{\mathbf{X}}_{*})\end{array}\right]\\ & -\left[\begin{array}{c}{\mathbf{K}}_u({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_u\\ {\mathbf{K}}_v({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_v\\ {\mathbf{K}}_w({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_w\end{array}\right]\\ & {\mathbf{K}}_f(\mathbf{X},\mathbf{X}{)}^{-\mathrm{1}}\left[\begin{array}{c}{\mathit{a}}_u\odot {\mathbf{K}}_u(\mathbf{X},{\mathbf{X}}_{*})\end{array}\right]\\ & \left[\begin{array}{cc}{\mathit{a}}_v\odot {\mathbf{K}}_v(\mathbf{X},{\mathbf{X}}_{*})& {\mathit{a}}_w\odot {\mathbf{K}}_w(\mathbf{X},{\mathbf{X}}_{*})\end{array}\right].\end{array}\end{array}$$

The mean of the posterior predictive distribution forms our estimate for the winds at the chosen points of interest, and this is given by the following.

$$\begin{array}{ll}\text{(35)}& {\displaystyle }\widehat{\mathit{u}}\left({\mathbf{X}}_{*}\right)& {\displaystyle }=\mathbf{E}\left[{\mathit{u}}_{*}\right|\mathit{f}]={\mathit{m}}_u({\mathbf{X}}_{*})+({\mathbf{K}}_u({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_u){\displaystyle }& {\displaystyle }{\mathbf{K}}_f(\mathbf{X},\mathbf{X}{)}^{-\mathrm{1}}(\mathit{f}-{\mathit{m}}_f\left(\mathbf{X}\right))\\ \text{(36)}& {\displaystyle }\widehat{\mathit{v}}\left({\mathbf{X}}_{*}\right)& {\displaystyle }=\mathbf{E}\left[{\mathit{v}}_{*}\right|\mathit{f}]={\mathit{m}}_v({\mathbf{X}}_{*})+({\mathbf{K}}_v({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_v){\displaystyle }& {\displaystyle }{\mathbf{K}}_f(\mathbf{X},\mathbf{X}{)}^{-\mathrm{1}}(\mathit{f}-{\mathit{m}}_f\left(\mathbf{X}\right))\\ \text{(37)}& {\displaystyle }\widehat{\mathit{w}}\left({\mathbf{X}}_{*}\right)& {\displaystyle }=\mathbf{E}\left[{\mathit{w}}_{*}\right|\mathit{f}]={\mathit{m}}_w({\mathbf{X}}_{*})+({\mathbf{K}}_w({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_w){\displaystyle }& {\displaystyle }{\mathbf{K}}_f(\mathbf{X},\mathbf{X}{)}^{-\mathrm{1}}(\mathit{f}-{\mathit{m}}_f\left(\mathbf{X}\right))\end{array}$$

Here we can see that the estimates near measurement locations, where ${\mathbf{K}}_{u,v,w}({\mathbf{X}}_{*},\mathbf{X})$ is large, are dominated by the prior covariance function specification. This is why the choice of prior covariance function is more important than the choice of prior mean function for making useful estimates and why our subsequent analysis is concentrated on the covariance hyperparameters.

Similarly, we obtain an estimate of the prediction uncertainty by using the posterior variance for each wind component, given by the following.

$$\begin{array}{ll}{\displaystyle }{\mathit{\sigma }}_{\widehat{\mathit{u}}}^{\mathbf{2}}\left({\mathbf{X}}_{*}\right)& {\displaystyle }=\mathbf{Var}\left[{\mathit{u}}_{*}\right|\mathit{f}]={\mathit{\sigma }}_u^{\mathrm{2}}-\mathrm{diag}(\left({\mathbf{K}}_u\right({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_u)\\ \text{(38)}& {\displaystyle }& {\displaystyle }{\mathbf{K}}_f(\mathbf{X},\mathbf{X}{)}^{-\mathrm{1}}({\mathit{a}}_u\odot {\mathbf{K}}_u(\mathbf{X},{\mathbf{X}}_{*})\left)\right){\displaystyle }{\mathit{\sigma }}_{\widehat{\mathit{v}}}^{\mathbf{2}}\left({\mathbf{X}}_{*}\right)& {\displaystyle }=\mathbf{Var}\left[{\mathit{v}}_{*}\right|\mathit{f}]={\mathit{\sigma }}_v^{\mathrm{2}}-\mathrm{diag}(\left({\mathbf{K}}_v\right({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_v)\\ \text{(39)}& {\displaystyle }& {\displaystyle }{\mathbf{K}}_f(\mathbf{X},\mathbf{X}{)}^{-\mathrm{1}}({\mathit{a}}_v\odot {\mathbf{K}}_v(\mathbf{X},{\mathbf{X}}_{*})\left)\right){\displaystyle }{\mathit{\sigma }}_{\widehat{\mathit{w}}}^{\mathbf{2}}\left({\mathbf{X}}_{*}\right)& {\displaystyle }=\mathbf{Var}\left[{\mathit{w}}_{*}\right|\mathit{f}]={\mathit{\sigma }}_w^{\mathrm{2}}-\mathrm{diag}(\left({\mathbf{K}}_w\right({\mathbf{X}}_{*},\mathbf{X})\odot {\mathit{a}}_w)\\ \text{(40)}& {\displaystyle }& {\displaystyle }{\mathbf{K}}_f(\mathbf{X},\mathbf{X}{)}^{-\mathrm{1}}({\mathit{a}}_w\odot {\mathbf{K}}_w(\mathbf{X},{\mathbf{X}}_{*})\left)\right)\end{array}$$

Since the measurement covariance K_f term includes the assumed measurement noise, these equations effectively propagate the Doppler uncertainty through the measurement geometry and meteor density to produce the wind estimate uncertainty. However, we note that this uncertainty estimate ignores the cross-terms in the covariance both between test locations and among the wind components. These factors can also be included to give a more complete picture of how the individual estimates are correlated at an increased computational cost. More detailed estimates could also be backed by a fully Bayesian approach that involves Markov chain Monte Carlo sampling of the posterior predictive distribution and includes full distributions for the hyperparameters θ.

Evaluating the posterior mean and covariance is a straightforward numerical linear algebra problem. However, given the potential sizes of the various covariance matrices, this can be computationally expensive. Mitigation of this implementation burden can be achieved with both matrix-free and approximate methods (e.g., Gardner et al., 2018; Wilson and Nickisch, 2015). Application of these methods are the subject of future work, but we note that their use would make practical fitting and evaluating more tractable.

4 SIMONe2018 campaign

Before describing and presenting the simulation and experimental results, in this section we briefly describe the SIMONe2018 measurement campaign that was conducted in northern Germany between 2 and 9 November 2018. As mentioned in the Introduction, the SIMONe2018 campaign added eight SIMONe links to six existing MMARIA links. The MMARIA links consist of two pulsed transmitters located in Juliusruh (13.37^∘ E, 54.63^∘ N) and Collm (13.00^∘ E, 51.31^∘ N) operating at 32.55 and 36.2 MHz, respectively. The signals of these transmitters were received at four receiving stations located in Juliusruh, Neustrelitz (13.07^∘ E, 53.33^∘ N), Bornim (13.02^∘ E, 52.44^∘ N), and Collm.

For the SIMONe links, a coded continuous wave (CW) transmitter was operated from Kühlungsborn (11.77^∘ E, 54.12^∘ N) at 32.55 MHz. The transmitter array consisted of five two-element single polarization antennas arranged in a Pentagon configuration. Each antenna transmitted a different pseudo-random code sequence, with 1000 bauds and 10 µs baud length. On reception, four single antennas were used, yielding MISO (multi-input, single-output) links. In addition, the same 32.55 MHz antennas and receiving systems located in Neustrelitz and Bornim were used to receive the coded CW signals, forming both MISO and SIMO (single-input, multiple-output) links at both sites.

The meteor signals from the pulsed links were detected and identified using a similar methodology as described in Hocking et al. (2001). In the case of the SIMONe links, the meteor signals were decoded and detected using the compressed sensing approach introduced by Urco et al. (2019a). Once the signals were detected, Doppler shift and interferometric angles were obtained from the autocorrelation and cross-correlation (between channels), respectively, in a similar manner as employed by Holdsworth et al. (2004). The interferometric angles were obtained using a combination of beam-forming and nonlinear complex fitting of the time series data following Clahsen (2018) and Chau and Clahsen (2019), which includes estimating statistical uncertainties for the Doppler measurements. Such uncertainty estimates are used as quality checks or weights in fitting procedures. Location of the meteors and representation of the Bragg vector in the local meteor ENU coordinate system were performed using the WGS84 representation for an ellipsoidal Earth coordinate frame. That procedure has been described previously in Clahsen (2018) and Stober et al. (2018). More details of the SIMONe2018 campaign can be found in Vierinen et al. (2019) and Charuvil Asokan et al. (2020).

5 Monte Carlo simulations

Monte Carlo simulations of the wind field ($\mathit{u}=u,v,w$) are essential to gauge the bias and variance properties of the GPR method. To create realistic random wind fields with which we could simulate meteor measurements and compare the GPR estimate, we again made use of Gaussian processes. Instances of $\mathit{u}(t,z,y,x)$ were drawn from the Gaussian random vector distribution described by Eq. (21) for specified sample locations, mean wind functions, and covariance amplitude and length scale hyperparameters. The hyperparameters used were as follows: ${\mathit{\sigma }}_u^{\mathrm{2}}={\mathit{\sigma }}_v^{\mathrm{2}}=\mathrm{900}$ m² s⁻², ${\mathit{\sigma }}_w^{\mathrm{2}}=\mathrm{90}$ m² s⁻², ${\mathit{\delta }}_x={\mathit{\delta }}_y=\mathrm{50}$ km, δ_z=3 km, and δ_t=1800 s. These velocity fields were used with observing geometries taken from 1 d of the SIMONe 2018 campaign, specifically 5 November 2018. At each real detection, the measured projected velocity was replaced by a new projected simulation velocity taking into account both the measured Bragg vector and the simulated $\mathit{u}(t,z,y,x)$. In this way, we are able to test the proposed GPR method on actual measuring geometries.

Using the simulated measurements, we followed the GPR method from Sect. 3 to estimate the 4D wind field for comparison to the simulated winds. We explored fitting with different cubic spline forms for the mean wind functions, and qualitatively we found that the wind estimates were not sensitive to the details of the fit as long as it was reasonable. Even using a constant mean of zero produced qualitatively similar results. Thus, to remove a confounding variable, all of the estimation results presented in this section use the exact mean functions that were used to simulate the winds, which in turn are the same mean functions fitted to the SIMON2018 data as described in Sect. 6. Likewise, we fit for the covariance hyperparameters from the simulated measurements and found that the results were similar (within 10 %) to the values used for the simulation. This was reassuring and showed that the fitting procedure works, at least when the winds can be described exactly by a Matérn covariance Gaussian process. Similar to the mean, the estimated winds showed little qualitative sensitivity to small changes in the covariance hyperparameters, so for the subsequent estimation results we used (as a baseline case) the same values for the amplitudes and length scales between the simulation and estimation Gaussian processes in order to remove fitting noise as a confounding variable. These comparisons should be viewed as a best-case scenario from the perspective of the model, and therefore they can be used primarily to explore the effects of meteor measurement spatial density and geometry on the quality of the wind estimates.

5.1 Qualitative comparison of horizontal winds

Figure 2 shows an example of results for simulated (left) and estimated (right) wind fields for three selected altitudes: 84, 90, and 96 km.

Figure 2 Simulated wind field (a, c, e) compared to the resulting GPR estimate based on SIMONe-derived measurements (b, d, f). Each panel shows the horizontal wind speed as a function of latitude and longitude overlaid by streamlines showing the wind flow. The estimated wind speed is masked at 50 % transparency in areas where there are few meteor detections, and thus the estimate uncertainty is relatively high (i.e., the improvement in posterior predictive variance over the prior variance is less than 4 dB).

The horizontal wind magnitude is color-coded (blue–green–yellow tones), while the direction is indicated by the over-plotted streamlines. The estimated values are also masked (altering transparency) in regions where the posterior predictive variances are high. Such regions are naturally where there are fewer meteor detections. Note that contrary to traditional methods and despite the presentation here as horizontal slices, the estimates are not confined to a regular horizontal grid since solutions are inherently obtained in 4D. At an overall level, there is very good agreement between the horizontal wind magnitude and direction at all altitudes in regions where the posterior predictive variance is reasonably low (full color areas).

5.2 Bias and error variance

For a more quantitative idea of the performance of the GPR method, we have repeated the Monte Carlo simulations 4700 times using 100 instances at each $(t,z,y,x)$ location for 47 different overlapping time intervals throughout the day. This is equivalent to observing over 100 d with the same measurement statistics at each of the 47 time intervals of a given day. We estimated bias and error variance by calculating the sample mean and variance of the error between the estimated and simulated u,v, and w wind values over the n=4700 time or trial instances. In the case of the horizontal winds, the bias is given as the magnitude of the mean error vector composed of both the zonal u and meridional v wind components, and the error variance is the sum of both the u and v error variances.

Figure 3 shows the bias of the horizontal wind error (left) color-coded with red tones and the error variance of the horizontal wind (right) color-coded with purple–yellow tones, in both cases for the same altitudes shown in Fig. 2.

Figure 3 Statistics of the horizontal wind estimator error relative to the simulated truth. Each panel shows the bias (a, c, e) or error variance (b, d, f) as a function of latitude and longitude averaged over n=100 trials at each of 47 measurement geometries taken throughout 1 d. Contours on the bias plots give the posterior predictive variance (in m² s⁻²), indicating more confidence in the central areas where the bias also tends to be a little lower. Contours on the error variance plots correspond to the sample error variance (matching the coloring).

In the mean error panels, the posterior predictive variance is also indicated with green contours. A bias of less than 2 m s⁻¹ is seen across the plots, and generally smaller biases are seen in the regions of lower predictive variance where there are more meteor detections. Note also that the uncertainty contours (left) roughly match the shape of the actual error variance (right), which gives confidence that the uncertainty estimates are useful.

Similarly, the bias and variance results for the vertical wind are shown in Fig. 4.

Figure 4 Statistics of the vertical wind estimator error relative to the simulated truth. Each panel shows the bias (a, c, e) or error variance (b, d, f) as a function of latitude and longitude averaged over n=100 trials at each of 47 measurement geometries taken throughout 1 d. Contours on the bias plots give the posterior predictive variance (in m² s⁻²), indicating more confidence in the central areas where the bias also tends to be a little lower. Contours on the error variance plots correspond to the sample error variance (matching the coloring).

Again, we see low biases that are uniformly less that 1 m s⁻¹ in magnitude, with the lowest biases in the regions of low predictive variance. However, this region is smaller than in the horizontal wind case. We are certain that this difference is mainly due to the configuration geometry that is needed to get accurate vertical winds, and the low-variance region provides a better observing geometry than the rest. Given the differences in magnitudes and the typically observed Bragg vectors, vertical wind estimates are relatively less constrained and more susceptible to horizontal wind contamination. Again, as in the case of the horizontal wind results, the uncertainty contours (left) roughly match the shape of the actual error variance (right).

Figure 5 Mean estimator error relative to the simulated truth when varying the covariance amplitudes. Each panel shows distributions of the estimator error averaged over n=100 random trials, for which the distribution is taken over estimates at time–space grid coordinates where the estimated uncertainty shows meaningful improvement (defined as 1.5 dB). Relative to the simulated values, the estimator covariance amplitudes were scaled by $\frac{\mathrm{1}}{\mathrm{2}}$, 1, and 2 to test nine different combinations by varying values for both the horizontal (${\mathit{\sigma }}_u^{\mathrm{2}}={\mathit{\sigma }}_v^{\mathrm{2}}=[\mathrm{450},\mathrm{900},\mathrm{1800}]$ m² s⁻²) and vertical (${\mathit{\sigma }}_w^{\mathrm{2}}=[\mathrm{45},\mathrm{90},\mathrm{180}]$ m² s⁻²) wind components.

5.3 Effects of scaling the covariance amplitudes

Until now we have presented results using estimator prior covariance amplitudes equal to the simulated values. In Figures 5 and 6, we show the biases and error variances while varying over different values of the estimator covariance amplitudes: (a) half, (b) equal to, and (c) double the true value of the simulated winds. Specifically, we took the same 47 observation windows as before, simulated 100 random trials of measurements using covariance amplitudes of ${\mathit{\sigma }}_u^{\mathrm{2}}=\mathrm{900}$, ${\mathit{\sigma }}_v^{\mathrm{2}}=\mathrm{900}$, and ${\mathit{\sigma }}_w^{\mathrm{2}}=\mathrm{90}$ m² s⁻², and estimated the winds with nine different covariance amplitude combinations by scaling the horizontal and vertical values separately by $\frac{\mathrm{1}}{\mathrm{2}}$, 1, and 2. Note that the horizontal amplitudes for the zonal and meridional wind components were varied together such that ${\mathit{\sigma }}_u^{\mathrm{2}}={\mathit{\sigma }}_v^{\mathrm{2}}$. Finally, we computed the error between the estimated and simulated winds, calculated the mean and variance of the error over the random 100 trials (to give bias and error variance, respectively), and plotted the resulting distributions taken over time–space grid coordinates.

Figure 5 shows the GPR bias statistics for the zonal (top), meridional (middle), and vertical (bottom) wind components, with columns corresponding to halved (left), equal (center), and doubled (right) covariance amplitudes for the given wind component.

The remaining vertical and horizontal covariance amplitude value is indicated with different colors. The salient features of this figure are the following: (a) the mean error has a tight distribution around zero, indicating little or no bias regardless of covariance amplitude scaling; and (b) the differences from scaling the covariance amplitudes are minor, with a slightly tighter bias distribution for the vertical wind component, a doubled vertical amplitude, and halved horizontal amplitudes.

The posterior predictive uncertainties are plotted against the error variance in Fig. 6 for both the horizontal (left) and vertical (right) wind components.

Figure 6 Estimator posterior uncertainty versus error variance relative to the simulated truth when varying the covariance amplitudes. Each panel plots the mean (lines) and 90 % confidence interval (shading) of the distribution of the posterior predictive variance versus the error variance calculated over n=100 random trials, for which the distribution is taken over individual estimates at time–space grid coordinates. Relative to the simulated values, the estimator covariance amplitudes were scaled by $\frac{\mathrm{1}}{\mathrm{2}}$, 1, and 2 to test nine different combinations by varying values for both the horizontal (${\mathit{\sigma }}_u^{\mathrm{2}}={\mathit{\sigma }}_v^{\mathrm{2}}=[\mathrm{450},\mathrm{900},\mathrm{1800}]$ m² s⁻²) and vertical (${\mathit{\sigma }}_w^{\mathrm{2}}=[\mathrm{45},\mathrm{90},\mathrm{180}]$ m² s⁻²) wind components.

In the horizontal case, we show the results of the total horizontal wind speed, i.e., $\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}$. Lines give the mean of the error variance distribution, while the shaded region indicates the 90 % confidence interval. For the horizontal and vertical wind plot, different line styles and labeling indicate the estimator values for the horizontal and vertical covariance amplitude, while different colors indicate values for the vertical and horizontal covariance amplitude, respectively. The estimator covariance amplitudes match the simulated covariance amplitudes at the middle-orange values shown (${\mathit{\sigma }}_u^{\mathrm{2}}={\mathit{\sigma }}_v^{\mathrm{2}}=\mathrm{900}$ and ${\mathit{\sigma }}_w^{\mathrm{2}}=\mathrm{90}$ m² s⁻²), and those cases show good linear agreement between uncertainty and error variance. Halving and doubling the prior covariance amplitude of a given wind component similarly scales the posterior estimator uncertainty, resulting in either underestimating or overestimating the uncertainty relative to the observed error variance.

Based on these Monte Carlo simulations, we recommend one of two approaches for applying GPR depending on the requirements of precision. First, if computational speed is a constraint and relatively large uncertainties are acceptable, then using conservative overestimates of the wind variances to specify the covariance amplitudes will still yield unbiased wind estimates with uncertainties that can be treated as rough upper bounds on the error variances. Second, if more precision is needed and computational time is not a problem, then fitting on the incoming data to get more accurate estimates of the prior covariance amplitudes will yield unbiased wind estimates with more accurate uncertainties. This choice between specifying the covariance hyperparameters and fitting for them is a critical decision for any user of the GPR method, as already seen in the block diagram of Fig. 1.

5.4 Qualitative role of the covariance length scales

We have not yet conducted a systematic study of the covariance length scales in the same manner as our examination of the covariance amplitude hyperparameters. This is because the degrees of freedom in perturbing the values are greater, making the analysis more complex, but also because the length scales are easier to interpret without detailed analysis. Because the model will enforce high correlation for coordinates that are “close” relative to the length scales, the covariance length scales set the effective resolution of the wind estimates. So intuitively, increasing the length scales will lose resolution and blur the estimates, while decreasing the length scales will gain resolution at the cost of increasing uncertainty (due to fewer measurements having a strong effect at a given estimation location). This intuition matches the informal testing that we have done in perturbing the length scales from the fitted values.

We have found that fitting the length scale hyperparameters generally does a good job of maximizing resolution while maintaining a usefully low posterior predictive variance. Those optimal values are determined by both the true covariance length scales of the wind field and the spatiotemporal density of the meteor measurements. For these simulated data, we know that the measurement density can support smaller length scales because the fitted values for the corresponding real data are roughly half for the x, y, and t dimensions (see Sect. 6) compared to the values for the simulated winds. Nevertheless, fitting the estimation hyperparameters to the simulated data produced length scales close to the simulation values, showing that the fitting is responsive to the “true” wind covariance distances and does not just tune to the meteor measurement density.

As an alternative to fitting, one always has the option of setting the covariance length scales according to a desired estimation resolution. This is useful when one is content with sacrificing potentially better resolution for the sake of computational simplicity. In the case that the measurement density is not high enough to support analysis at those fixed length scales, that fact will be made clear by having few or no regions of low posterior predictive variance for the resulting winds. The estimates will likely not have the overall best uncertainty, but they will still be valid and thus useful.

Figure 7 Latitude–longitude slices of the winds estimated from SIMONe campaign data. Each panel represents a separate altitude and time and shows the horizontal wind speed as a function of latitude and longitude overlaid by streamlines, which show the wind flow. The wind speed is shown with 50 % transparency in areas where the estimate uncertainty is large (<2 dB improvement relative to prior uncertainty, i.e., where there are few meteor detections).

6 Experimental results

In this section we implement the proposed wind field estimator on a dataset of 24 h observations collected on 5 November 2018 during the SIMONe2018 campaign. After initial data quality control, almost 200 000 meteor detections were obtained in 24 h. Using a conservative approach and performing further quality checks yielded 100 000 high-quality detections. The filter criteria used in this second reduction required that detections were (a) within 3 standard deviations of the zero-order residuals and (b) more than 30^∘ above the horizon to ensure that good interferometric angle of arrival (AOA) or angle of departure (AOD) estimates were obtained (e.g., Chau et al., 2019). Filtering by a minimum elevation angle also has the effect of ensuring that the errors in AOA and/or AOD, when projected into the vertical direction, have a limited effect on the estimated altitude. Meteor location errors are not incorporated into the current GPR method, so their effect must be limited by ensuring that any potential coordinate deviations are much smaller than the covariance length scales used.

Subsequently, GPR results were obtained by first determining mean wind functions by fitting a 6 knot (altitude) by 6 knot (time) tensor product cubic spline over the entire 24 h of data. The 12 spline parameters were calculated by solving the standard least squares problem completely independently of the GPR model. Then the covariance fitting procedure was applied to overlapping 90 min windows spaced at 30 min intervals to estimate the covariance amplitudes and length scales as they varied throughout the day. With the current procedure that computes the full covariance matrix, limiting to short time intervals like this is necessary for computational feasibility. The hyperparameters were found to be constant enough throughout the day that approximate overestimates would suffice and allow proceeding with a single set of hyperparameters. The resulting covariance hyperparameters are ${\mathit{\sigma }}_u^{\mathrm{2}}={\mathit{\sigma }}_v^{\mathrm{2}}=\mathrm{900}$ m² s⁻², ${\mathit{\sigma }}_w^{\mathrm{2}}=\mathrm{90}$ m² s⁻², ${\mathit{\delta }}_x={\mathit{\delta }}_y=\mathrm{26}$ km, δ_z=3 km, and δ_t=900 s. Finally, the wind estimates were produced by selecting a fixed time, gathering data from the 90 min window around that time (more than enough given the time length scale of 15 min), and computing the posterior predictive values at chosen spatial points.

To get a sense of the scales resolved with the GPR method, Fig. 7 shows latitude–longitude slices of wind fields at three different altitudes (84, 90, and 96 km) and three different times (05:00, 08:00, and 11:00 UT).

Figure 8 Altitude–time slices of the winds estimated from SIMONe campaign data. Zonal and meridional winds are shown at a selection of four latitude–longitude points. The wind speed is shown with gray shading in areas where the estimate uncertainty is large (<1 dB improvement relative to prior uncertainty, i.e., where there are few meteor detections).

The presentation format is similar to Fig. 2; i.e., horizontal wind speeds are color-coded, and streamlines show the direction of flow. Areas of large velocity variance are shaded with 50 % transparency to white. The wind fields show significant complexity, much more than can be well represented by the single mean vector per plot that would be reported by a monostatic meteor radar. On simple inspection, horizontal wind structures of ∼ 20–50 km are successfully resolved, which is commensurate with the horizontal length scale hyperparameter of 26 km.

In Figure 8, altitude–time slices at selected latitude–longitude points are shown for both zonal (left) and meridional (right) wind components.

Figure 9 Comparison of GPR to gradient and homogeneous methods. Panels (a)–(c) show the horizontal wind field obtained with the gradient method using 4 h and 4 km bins. Panels (d)–(f) show the horizontal wind field obtained with the GPR method using fitted covariance hyperparameters. Panels (g)–(i) show the wind field difference between the values in the second row and the mean horizontal wind indicated in all panels with a black arrow. In all cases, normalized statistical variance is indicated as gray contour lines, while the color contour represents the vertical component from the gradient method (a–c), GPR method (d–f), and GPR minus the mean from the gradient method (g–i). The row 2 color bar corresponds to the background vertical wind coloring, while the other two color bars correspond to their respective arrow colors.

The large-scale tidal features are in good agreement with those obtained with the homogeneous method applied to the same data (see Vierinen et al., 2019, Fig. 6). The winds show significant variation between horizontal locations as expected.

Although we do not have a ground truth in this analysis to validate the horizontal scales we are resolving, we conduct an additional comparison to complement earlier identification of the large-scale features (i.e., tides). In Fig. 9, we compare GPR wind fields with those obtained with the homogeneous method (i.e., independent of latitude and longitude) and those obtained with a gradient method. Specifically, the homogeneous method uses a zero-order Taylor expansion, while the gradient method uses a first-order Taylor expansion. Both estimates have been obtained with altitude and temporal bins of 4 km and 4 h, respectively, in order to produce a good representation of large-scale features. The specifics of the two methods can be found in Chau et al. (2017) and Chau et al. (2021), respectively.

The gradient wind fields are shown in the first row of Fig. 9 for three selected altitudes (84, 89 and 94 km).

The arrows are color-coded with the horizontal wind speed (green tones), while the mean vertical wind from the gradient method is color-coded with red–yellow–blue tones. In the second row the GPR 3D wind fields are displayed in a manner similar to the gradient estimates in the first row. The third row shows the difference between the GPR wind fields and those from the gradient method. Note that the arrow colors and color bar in the third row are different from the first two rows and show the difference of the horizontal winds. In all three rows the horizontal wind from the homogeneous method is shown with a thick black arrow in the center.

The salient features of Fig. 9 are the following.

• In general, there is good agreement in the horizontal wind components between the gradient and GPR methods. Note that the gradient estimates have been obtained with relatively large temporal and vertical averaging in order to produce a good representation of large-scale features.

• By subtracting the mean wind obtained with the gradient method (i.e., large-scale features) from the GPR estimates, in the third row, mesoscale structures are identified. Horizontal structures of the order of 20–50 km are clearly identified in all three altitude cuts.

Similar wind field comparisons for different times of the day can be found in Supplement Movie S1.

7 Discussion

We have introduced a robust method based on Gaussian process regression analysis to estimate MLT wind fields in four dimensions. The method has been evaluated using Monte Carlo simulations and implemented successfully on real data. The fast implementation using specified covariance hyperparameters (per-component amplitudes and per-dimension length scales) provides unbiased estimates with estimated uncertainties proportional to the prior velocity variances. In other words, if the prior variances are underestimated, the posterior variances are also underestimated. Using a more resource-intensive training and fitting approach, covariance amplitudes can be estimated, resulting in posterior variances that are in good agreement with expectations from Monte Carlo simulations. The training approach requires more computation time than using fixed prior variances, and we have not routinely applied it in analysis to date. However, for method testing purposes, we have implemented it on the real data shown in this work.

As expected, we have shown that mean values of GPR wind fields are in good agreement with the mean winds obtained with the homogeneous method. Similarly, to a first-order approximation, GPR wind fields are also in good agreement with the wind fields obtained with the gradient method. Based on the simulation results, we expect the differences (i.e., the 20–50 km scales within their posterior variances) to be of a geophysical nature.

Although the GPR method is robust, its region of validity and resolution depends highly on the geometrical configuration used, which influences the location and density of meteor observations and the observable projected wind component. For example, we found that the region of low-variance vertical winds is smaller than the region of low-variance horizontal winds. This result occurs even though the SIMONe2018 configuration has far superior properties in terms of links and diversity of Bragg angles compared to any other multistatic configuration used to date to study MLT winds (e.g., Chau et al., 2017; Stober et al., 2018; Spargo et al., 2019; Chau et al., 2021; Conte et al., 2021). Fortunately, the posterior predictive variances provided by the GPR method can be used in the future to optimize the meteor radar network geometry to achieve a given prediction goal, e.g., covering a specified region so that the estimate uncertainty for the winds reaches a particular value given typical meteor statistics.

Estimating the vertical wind component is still challenging due to two factors: the horizontal wind variability is larger than the vertical wind variability (leading to large contamination of the vertical wind when there are errors in the estimated Bragg vector or meteor location), and the majority of Bragg vectors have angles that are not close to zenith. The absence of zenith-oriented Bragg vectors is intrinsic to all specular meteor radars, since any Bragg vectors with angles close to zenith would require meteor trajectories parallel to the Earth's surface and are therefore very unlikely to be observed. In the particular case of the gradient method, Chau et al. (2017) have previously shown that the mean vertical velocity obtained with the homogeneous method, i.e., an area of ∼ 200 km radius, was contaminated by the mean horizontal divergence. Similar effects would be expected at smaller scales. Our experimental results do produce a vertical wind prior variance of about 90 m² s⁻², and some of the vertical wind estimates do show nonzero vertical velocities congruent with that variance. However, the posterior error bars are still large enough that a zero or nearly zero vertical wind is a plausible explanation, especially considering the possible role of horizontal contamination. The important points relevant to the technique are that GPR is agnostic to the prior assumptions one wants to employ for the vertical winds, and it also provides the necessary uncertainty information to allow for assessing the quality of the vertical wind estimates.

These results represent just the first step toward applying GPR analysis to estimate wind fields from meteor observations. We envision multiple directions of future work to expand and improve on the technique. There are many degrees of freedom in specifying mean and covariance functions to represent the wind components that can be explored. Known physical processes imply more structure in the joint wind component covariance than expressed in Eq. (21), so it would make sense to experiment with adding cross-covariance terms and allowing independent length scales for each component. The spatiotemporally varying sampling density imposed by the meteors argues for using covariance functions or hyperparameters that also vary in time and/or space. This can already be achieved in a crude form by performing fitting and estimation on overlapping subsets of the data, and we would like to explore that more as well as to develop a more elegant approach. We have used the mean functions to essentially remove large-scale tidal effects, but it remains to be seen how to strike the optimal balance between complexity in the mean versus covariance functions or even the model complexity overall. At some point, adding complexity transforms the GPR method from data-based estimation into assimilative modeling, and we see value in prioritizing simplicity and clarity.

Incorporating the uncertainty in the meteor locations and Bragg vector components into the GPR analysis is another important avenue for improving the technique. We have so far removed any low-quality meteor detections from the analysis to limit the effect of this additional error, and the quality of the wind estimates would be improved by being able to incorporate these discarded data and make even better use of the high-quality detections. We anticipate that such a task would be challenging; it would likely entail leaving the closed-form solutions behind and numerically sampling from the distributions (e.g., Markov chain Monte Carlo methods).

Future work will also concentrate on further validation (including cross-validation within a single dataset), although the fact remains that no alternative MLT wind instrument is currently available for comparison with GPR estimates. Therefore, independent of the good comparisons with Monte Carlo simulations, we are planning to conduct special future observing campaigns under different atmospheric conditions and geometric configurations to intercompare our GPR method with other wind field methods such as those employing Tikhonov regularization (e.g., Stober et al., 2018; Chau et al., 2021). Similarly, we plan to compare these techniques using synthetic data from regional weather models with high resolution covering the MLT altitudes, such as the ICON-UA model (e.g., Borchert et al., 2019). This analysis concept would be similar to the one implemented in this work, but with more realistic atmospheric dynamics for the simulated winds.

Finally, we plan to apply the GPR method to selected additional datasets that use a multistatic configuration in order to further investigate the properties of the resolved 20–50 km horizontal wind structures. These investigations will cover both individual case studies and statistical studies: for the former, we expect to analyze special geophysical conditions and/or measurements that are complemented by other ground- or satellite-based instruments (e.g., Davis et al., 2018; Vargas et al., 2021); for the latter, we expect to compare the Reynolds stress tensor statistics of GPR-estimated wind fields to those obtained from second-order statistics of projected wind velocities (Vierinen et al., 2019).

8 Conclusions

We have introduced an alternative observation method based on Gaussian process regression analysis to resolve MLT wind fields in 4D from multistatic radar observations. Based on Monte Carlo simulations of known wind field distributions, our proposed method provides unbiased mean velocity estimates and posterior velocity variances that are proportional to prior velocity variances. By using an adaptive fitting procedure based on input data, unbiased posterior variances can be achieved. This adaptive approach is currently not practical for real-time applications but is ideal for case studies.

The horizontal regions of good GPR method performance in MLT wind determination are dependent on the meteor scatter geometric configuration. On one hand, optimal configurations should ultimately increase the number of detections. However, on the other hand, these same configurations need to provide sufficient Bragg vector diversity. For the particular SIMONe2018 experiment scattering geometry, these factors meant that vertical velocity estimates with relatively small variances were obtained over a much smaller horizontal area than horizontal wind estimates.

Overall, the GPR method has attractive benefits for MLT regional and weather studies: (1) it enables flexible analysis by allowing grid-free wind estimates; (2) it provides statistical uncertainties for the estimated winds that reflect measurement uncertainty and meteor observation geometry; and (3) it adapts to the horizontal, vertical, and temporal scales of the data, accounting for measurement density, and is thus able to resolve winds at relatively small scales.