3.1 Implementing the Earth geometry in the wind analysis

Figure 2Illustration on how local coordinates (zonal, meridional) change with geographic position with respect to a radar location.

Download

A new aspect of the MMARIA concept is the necessity to consider the geometry of the Earth. At present our domain area in Germany has a horizontal extension of approximately 600 km×600 km. These distances are too large to assume a plane geometry, which is very often the case for classical monostatic MRs, where all observations are referred to the location of the radar itself. However, it is straightforward to at least consider that the altitude or height above the surface needs to be corrected for the Earth curvature using a mean Earth radius (R_E=6 378 137.0 m). Further, it is possible to obtain a local elevation angle for each meteor, instead of the one observed relative to the receiver location. However, this simple correction turned out to still be insufficient for dealing with large domain areas. Therefore, we outline a more detailed procedure, taking into account the Earth shape with the WGS84 geoid model (National Imagery and Mapping Agency, 2000).

In the following we outline the procedure how to compute new local coordinates (ENU: east–north–up) for each meteor to reduce potential errors in the wind field estimation due to projection issues. Considering that the Earth is not a perfect sphere, we have to deal with two different coordinates, the geodetic coordinates (longitude and latitude) and the Earth-Centered, Earth-Fixed (ECEF) coordinates, also called geocentric coordinates. The geodetic coordinates are determined by the normal to the ellipsoid, whereas the ECEF coordinates are defined by the Earth center using a (X, Y, Z)-coordinate system. Thus, the geodetic and geocentric latitude can be different.

We need to transform the observed meteor positions relative to the radar into a local coordinate frame (ENU) by determining the geodetic longitude, latitude and height of the meteor itself. The corresponding transformations are listed in the appendix. The procedure contains four steps: firstly the geodetic coordinates of the radar (longitude-ϕ,latitude-λ, height-h) are converted into ECEF, which means we receive a vector ${\mathit{x}}_R=(X_R,Y_R,Z_R)$, with respect to the Earth center. From the interferometric solution we also know the position vector ${\mathit{x}}_P=(x,y,z)$ with respect to the receiver location for each meteor. This vector x_P is then transformed into ECEF coordinates and we obtain a vector pointing from the Earth center towards the meteor position with coordinates ${\mathit{x}}_M=(X_M,Y_M,Z_M)$. Further, we convert the vector x_M given in ECEF coordinates back into a geodetic position given by the latitude, longitude, and height above the Earth's surface for each meteor. Finally, we use the ECEF vector x_M and the geodetic reference to compute a local coordinate set ENU at the position of the meteor itself. A detailed description of all applied coordinate transformations is summarized in the appendix.

In summary we perform the following steps:

1. conversion of the geodetic coordinates of the radar $(\mathit{\varphi },\mathit{\lambda },h)\to (X_R,Y_R,Z_R)$ by using the transformation geodetic to ECEF;

2. transformation of meteor coordinates into ECEF $(x,y,z)\to (X_M,Y_M,Z_M)$ by using the transformation of ENU to ECEF;

3. conversion of ECEF frame meteor position into geodetic coordinates $(X_R,Y_R,Z_R)\to (\mathit{\varphi },\mathit{\lambda },h)$ ECEF to geodetic;

4. determination of local ENU using the geodetic position of each meteor in ECEF coordinates $(\mathit{\varphi },\mathit{\lambda },h)\to (x_M,y_M,z_M)$.

Figure 2 shows an example of the difference between the ENU coordinates (black cross) of the radar location and the ENU coordinates (red cross) of a meteor observed at a horizontal distance of 300 km. The blue circle marks the 300 km range around the radar, which is assumed to be located at Juliusruh. Although the difference appears to be small, it introduces an error of a few meters per second in the derived zonal and meridional and vertical wind speeds. Depending on the range and geographic latitude of the measurements the local azimuth and zenith shows differences up to 4^∘ compared to the radar site. As our wind measurements are supposed to be aligned along the zonal and meridional direction it is beneficial and straightforward to remove this bias. In the case of the standard SMR wind analysis technique (Hocking et al., 2001), where a homogenous wind is assumed within the observation volume, the error is almost compensated due to the large number of meteors used for the wind estimate. However, it turns out that the random distribution of meteors within the observation volume is sometimes not favorable to compensate for this bias, thus, it also has an effect on the standard analysis. In particular, altitudes where only a small number of meteors are used for the wind estimation procedure are more prone to this type of error. In particular, it appears to be very critical or in fact almost impossible to obtain a reliable vertical wind velocity or momentum flux, if the full Earth geometry is not taken properly into account.

3.2 Retrieving arbitrary non-homogenous wind fields

The retrieval of arbitrary and non-homogenous wind fields is mathematically more demanding as the number of unknowns exceeds the number of measurements, which does not allow to directly solve the equations applying standard least squares or singular value decomposition algorithms. However, it is possible to constrain the problem by additional assumptions or a priori knowledge. Very often the smoothness is used to regularize the problem so that it can be solved by applying statistical inversion algorithms (Aster et al., 2013).

At first we define a spatial grid and a domain area. In the case of the German MMARIA network we use a 30 km×30 km grid spacing in zonal and meridional direction. The total domain area is about 600 km×600 km. The spatial grid is fixed for each altitude and follows the Earth's surface. A schematic view of the spatial grid is given in Fig. 3. There are many other possibilities to define the spatial grid, e.g., using fixed longitude and latitude bins or arbitrary grids by using each individual meteor position. It turns out that spatial grids with fixed horizontal distance have benefits for the diagnostic of the wind fields as they allow to use discrete Fourier transforms or wavelength based spatial spectral analysis techniques. In this study, we have adopted a regular grid.

Figure 3Schematic of 3-D gridding to compute horizontally resolved wind fields including the Earth surface.

Download

A first step of the wind field inversion procedure is to assign each observed meteor to a grid cell j centered at time t and position x_j. This is equivalent to averaging measurements in time and space. The temporal and vertical weighting due to the random occurrence of the meteors is done similarly to the general case (see above). At present we still use a mean total temporal and vertical shear amplitude for the whole domain area. In order to take into account that each observed meteor does not occur exactly at the position of the grid point x_j, which denotes the center of our grid cells, we assign a weight to each observation (Shepard, 1968), i.e.,

The weight w_i for meteor i at position x_i and at time t_i is inversely proportional to its distance to the center of the grid cell j. The exponent p is used to control how fast the weight is reduced as a function of distance. Assuming a value p=0 results in a box car with equal weight for each meteor independent of its distance from the grid cell center. We use the distance in meters and p=0.2. The value of p is selected in such a way that a meteor at 30 km distance from the grid cell center enters the retrieval with a non negligible weight. The term γ is used to control the slope of the space distance, we use γ_x=1.0. The main reason for the averaging is that a single meteor, which lasts for 20–200 milliseconds, does not necessarily provide a representative mean wind velocity for a 30 min time bin at a grid point or cell. At least two meteors have to occur within one grid cell, otherwise we do not attempt to estimate the wind speed for that grid cell. We will discuss later how this weight is applied in the inversion procedure.

We can relate the measured radial velocity of each meteor to the three dimensional wind velocity by using local ENU coordinates for each measurement i as follows:

where v_r,i is radial velocity for meteor i; $u_j,v_j,$ and w_j are the zonal, meridional, and vertical wind components in grid cell j, corresponding to the meteor location. The angles θ_i and ϕ_i, corresponding to meteor i, are the local ENU coordinates corresponding to the line of sight velocity along the direction vector from the radar. In the case of a forward scatter geometry the measurement of the radial velocity is more complicated. For multi-static geometries it has to be considered to obtain the effective Bragg wavelength and pointing vector (see Stober and Chau, 2015).

The radial wind equation for arbitrary measurements and grid cells can be expressed as a linear matrix equation. The mapping from the zonal, meridional, and vertical components to observed radial velocities is given by a geometry matrix $\mathbf{G}\in {\mathbb{R}}^{n\times \mathrm{3}m}$. All the measurements during an analysis interval are represented as a vector ${\mathit{v}}_r\in {\mathrm{R}}^{n\times \mathrm{1}}$, where n is the number of measurements. The radial velocity vector v_r contains all observed radial velocities, either for each individual meteor weighted by its distance from the grid cell, further referred to as “full wind retrieval”, or an averaged value that is already interpolated to the defined grid cells, further referred to as “packed wind retrieval”. The packed retrieval has the advantage that the inversion matrix only scales with the number of grid cells and not with the number of meteors, as is the case for the “full wind retrieval”. The unknown 3-D wind components for each grid cell are also expressed as a vector

where m is the number of grid cells. The mapping in Eq. (7) can be compactly expressed using the following matrix equation:

which relates all measured radial velocities to 3-D velocities within a grid. More explicitly, this is

The geometry matrix G combines all measurements from all possible viewing geometries, but it is not directly invertible. Although we have several different viewing geometries, we do not get always three independent measurements for each grid cell. Hence, the number of unknowns is still larger than the number of measurements (rows in matrix G). This is the case in particular for the edges of our domain area.

Ill-posed problems can be solved by adding additional constrains. Very often the smoothness of the unknown provides a reasonable way to regularize an ill-posed problem (Aster et al., 2013). The smoothness in our case corresponds to the spatial derivative for each wind component and grid cell. This is equivalent to the assumption that the wind field is only slightly changing between neighbored grid cells. Hence, we define a smoothness matrix $\mathbf{L}\in {\mathbb{R}}^{\mathrm{3}m\times \mathrm{3}m}$ in such a way that we couple neighbored grid cells for each velocity component separately. For one velocity component, and one grid cell, the elements of matrix L would be

The matrix L contains such differences for all grid cells and all velocity components. The L matrix is constructed that the total difference between a grid cell and its for neighboring grid cells is small and gradients within the x–y plane are not damped. In Fig. 4, a scheme of how the smoothness matrix is constructed is shown. The number of neighbor grid cells defines the number of entries in L for each grid cell.

Figure 4Schematic of L matrix for a center grid cell.

Download

Finally, we have to deal with grid cells in the domain area where no measurement is available for a given time–altitude bin. This issue is solved by introducing a mesoscale wind field solution to these grid cells. We tested three possible mesoscale solutions and checked how much the final solution depends on this mesoscale boundary condition. The most trivial way is zero padding or simply not using an explicit a priori for these grid cells, the second one is estimating a mean wind using all radial velocity measurements with the “all-sky” fit and the third possibility is to derive a mesoscale wind field solution by computing local mean winds for each multi-static geometry and to estimate a distance weighted background wind field for each grid cell. A similar result is achieved by applying volume velocity processing (VVP) (Browning and Wexler, 1968; Waldteufel and Corbin, 1979), which was already successfully applied using horizontally resolved radial wind measurements (Stober et al., 2013) or multi-static SMR observations (Chau et al., 2017; Stober and Chau, 2015). The mesoscale solution is considered in the retrieval as diagonal matrix prescribing a velocity (zonal, meridional, vertical) for each grid cell. The mesoscale regularization uses a rather low regularization strength α_meso, which is one order of magnitude smaller than the L regularization for the cases presented here.

Figure 5Comparison of 2-D wind fields for three different α=0.014, α=10, and α=0.000001. The length of the arrows between the images does not scale.

Download

Combining all the information and the smoothness constraints into a set of equations allows solving the ill-posed problem. We obtain an estimate for the 3-D wind components $\widehat{\mathit{u}}$ at all grid cells solving the equation

which is a standard regularized weighted linear least-squares estimator (Aster et al., 2013). The matrix $\mathbf{\Sigma }=\mathrm{diag}({\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}},\mathrm{\cdots },{\mathit{\sigma }}_n^{\mathrm{2}})$ is a diagonal matrix that contains the variance (i.e., measure of uncertainty) given to each measurement ${\mathit{\sigma }}_i^{\mathrm{2}}$. The regularization parameter α provides a weight to the regularization constraint. It describes the coupling strength between neighboring grid cells. It should be noted that a number of alternatives to regularizing the solution also exist.

In the inversion process the weight Σ⁻¹ consists of contributions due to the radial velocity uncertainty, the spatial and temporal weighting functions as outlined above, the angular uncertainties in the measurements of ϕ and θ and the errors in the zonal, meridional vertical wind components after each iteration step. There are grid cells where a sufficient number of measurements is not available (more than 2 meteors are required) and only the mesoscale wind field is used. In this case, the measurement is weighted by a large uncertainty (${\mathit{\sigma }}_i=\mathrm{200}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) to ensure that this does not strongly bias the inversion. However, there are areas where ,due to the radial nature of the meteor measurements, the solution can differ rather significantly due to the mesoscale regularization.

In the following we are going to demonstrate the robustness of our algorithm, independent of the choice of regularization parameter α and the mesoscale boundary conditions. For simplicity we will only focus on horizontal winds, and leave the discussion of vertical mesoscale wind velocity for future work. In order to obtain reliable and physically meaningful vertical velocities, additional regularization constraints might be required.

3.3 Robustness of wind field solution

Solving Eq. (12) is straightforward using singular value decomposition or matrix inversion algorithms. As the wind field inversion is still a linear problem, we just need to find a proper solution for the regularization parameter α. A large α means that our solution is dominated by the regularization constraint, a too small α results in a too weak coupling between neighboring grid cells, making the solution unstable as erratic points start to dominate. This is usually expressed in the so-called “L-curve” (Aster et al., 2013).

Figure 6Comparison of wind field solution in dependence of the background mesoscale solution. (a, b) Test case with zero padding: first (a) and final iteration (b). (c, d) Test case with mesoscale wind field (obtained from the data): first (c) and final iteration (d).

Download

In Fig. 5, we compare the obtained wind fields for different strengths of the regularization parameters α. The left picture shows what we consider as solution with α=0.014, which seemed to be sufficient in this case. The wind field in the center was computed assuming a much too strong α=10. This obviously leads to a much too smooth wind field, but still keeps some mesoscale wind field structure. A much higher value in the order of 100 or 1000 is going to further reduce the shown variability. The right picture shows an example with a regularization constraint of α=0.000001 that is intentionally much too small. This obviously leads to some erratic structures and outlier solutions begin to dominate the wind field.

We tested different possibilities to define an optimal regularization strength α. At present we optimize our solution with a global estimate that is valid for the whole domain and constant in time, instead of estimating a local regularization strength α for each grid cell. The local approach did not suppress erratic structures or outliers in the same way, but is going to be further considered as it is going to be beneficial for larger domain areas. After comparing thousands of images using different strengths and ways of estimating the optimal α, it turned out that α=0.1 is very often a useful value, which leads to a convergence within eight iterations. However, the choice of α depends on the used statistical weights and uncertainties entering the retrieval.

As already mentioned above, there are grid cells where we have no direct measurements from one of the systems. We suggested using a mesoscale solution for these points. Now there is the question whether our solution depends on this pre-described mesoscale wind field. Therefore, we prepared two test cases. In the first test, all grid cells with no direct measurements are zero padded. For the second test, we use a computed mesoscale wind field estimated from VVP. Figure 6 shows four pictures using the two test cases. Different colors label grid cells with direct radial velocity measurements (blue) and grid cells with the mesoscale solution (red). The left plot displays the first iteration step and the right one shows the finally obtained wind field solution. North of 52^∘ N there are almost no differences of the solution if one just compares the blue arrows. The main reason for this good agreement is that almost all points are linked by multiple observing geometries, whereas south of 52^∘ N we basically have only monostatic observations. As a result the obtained wind field in the southern part of the domain area is more prone to be affected by the boundary conditions. However, as there is almost no visible difference between the wind fields at latitudes north of 52^∘ N, we conclude that there is almost no impact on the determined wind fields by the pre-described mesoscale winds.

Figure 7Comparison of wind retrievals. Panel (a) shows the solutions of the “packed wind retrieval”. Panel (b) shows the same time and altitudes for the “full wind retrieval”.

Download

Further, we investigated the differences between the full wind retrieval and the packed wind retrieval. Therefore, we kept the regularization strength fixed for both retrievals. The packed retrieval makes use of the total wind variances for each grid cell and a mesoscale regularization, whereas the full wind retrieval uses each individual radial velocity uncertainty and no explicit mesoscale regularization. The other weights and the error treatment is the same. Figure 7 shows a sequence of three successive time steps for both retrievals. The difference in coverage is because we reduced the limit for the full wave retrieval to one meteor per cell, which would be not sufficient for the packed wind retrieval. Further, the packed wind retrieval indicates the impact of the mesoscale solution at the edges of the domain area. However, the general shape of the wind field is reproduced by both approaches with a remarkable agreement. This is further underlined by a correlation density plot shown in Fig. 8. The scattering around the center line is mainly due to the different weights and the mesoscale regularization, which was used in the packed wind retrieval compared to the full wind retrieval.

Figure 8Comparison of zonal and meridional winds as derived from the full wind retrieval algorithm and the packed wind retrieval.

Download