Reply to Reviewer comments on Manuscript doi:10.5194/amt-2016-365-RC1,2017 Mean wind vector estimation using the Velocity-Azimuth-Display (VAD) method: An explicit algebraic solution

The paper deals with a well-known practical issue of determining a vertical wind profile from a set of known radar measurements. The authors are able by using the measurement directions as a frame to provide a new derivation and an explicit solution (17). The outcome is similar to a Fourier transform of the equidistant radial velocities. The results that the vertical and horizontal winds are the first Fourier components of the azimuthal wind field is wellknown although the Fourier expansion is not usually presented explicitly (but see Browning Wexler, page 107). Hence the analytical result is not really new. We thank the reviewer for pointing out that an example for the explicit solution for the wind vector retrieval in our equation (17) is given in [Browning and Wexler(1968)], namely between their equations (7) and (8). However, these authors have not mentioned the rather general significance of these formulae since they appear only for a special case and are neither discussed nor numbered. We can only speculate whether or not Browning and Wexler fully appreciated the importance of their example and would therefore respectfully disagree with the referee that equation (17) in our paper is well-known. The derivation of an explicit formula has, to the best of


Introduction
The wind vector is a fundamental variable for describing the state of the atmosphere (Dutton, 1986); it can be measured with a variety of in situ and remote sensors.The latter are typically ground-based active systems emitting artificially generated electromagnetic (lidar, radar) or acoustic waves (sodar).In the so-called Doppler systems, the frequency shift of the scattered waves is used to measure the motion of the scattering medium directly.While the details of the measurement process differ between the various instruments, the common feature of all Doppler instruments is that the velocity can only be determined along a single direction, called the line of sight or radial direction.This provides merely one component of the full 3-D wind vector.Of course it is possible to use three Doppler instruments to sample the same volume from different directions, but this approach is impractical for operational meteorology (Stephens, 1994) and only used for special applications (Fuertes et al., 2014).Vertical wind pro-filing attempts to estimate the wind vector as a function of height using data from a single instrument.
A number of different retrieval methods have been proposed for uniform or linear wind fields (Browning and Wexler, 1968;Waldteufel and Corbin, 1979;Koscielny et al., 1984;Caya and Zawadzki, 1992;Stephens, 1994).These methods are known as velocity-azimuth display (VAD), volume velocity processing (VVP) or Doppler beam swinging (DBS), and its different variants are successfully used in operational meteorology.
It is obvious that such simplifying assumptions can in general not be made for the instantaneous wind field in a turbulent flow.However, it is customary to decompose a turbulent wind field into a mean and a fluctuating component describing the deviations from the mean (Salby, 1996;Davidson, 2004;Vallis, 2006).This is justified by the claim of an existing spectral gap between the mean flow and (microscale) turbulence (Stull, 1988).Indeed, evidence for such a spectral gap in the boundary layer was recently reported by Larsén et al. (2016).The wind measurement can thus be split up into two tasks, namely first the determination of the mean value and, second, the estimation of the Reynolds stress tensor or other statistical parameters describing the turbulent part of the flow (Sathe and Mann, 2013;Sathe et al., 2015;Newman et al., 2016).
For operational applications, like data assimilation for numerical weather prediction models, the interest is clearly in the mean wind.This is due to the fact that such models are unable to resolve the small (turbulent) scales directly.Processes on such scales must instead be parame-terized (Warner, 2011).For operational Doppler profilers it is therefore enough to aim at the determination of the mean wind vector profile, with typical averaging times of O(10 min).This restriction makes it more likely that simplifying assumptions like homogeneity or linearity of the wind field hold at least on average without incurring large errors.In fact the assumption of statistical (local) homogeneity and quasi-steadiness can often be applied in boundary layer meteorology (Wyngaard, 2010) even though it is clear that there are limits to these assumptions (Maurer et al., 2016).
Practical experiences from comparisons of various wind retrieval methods with Doppler radars suggest that the simplest methods for the retrieval of the horizontal winds give the best accuracy in comparison with independent wind sensors (Holleman, 2005) -a seemingly counterintuitive result, given the large area scanned in comparison with special wind profiling instruments like radar wind profilers or Doppler lidars (Cifelli et al., 1996).A possible explanation is that the retrievals using more complex wind field models are ill-conditioned (Shenghui et al., 2014).Given these results and the importance for wind profile measurements for operational meteorology, it seems therefore appropriate to further investigate the wind retrieval methods for the rather simple assumption of horizontal homogeneity and quasi-steadiness (or stationarity).
The paper is organized as follows: the first section is concerned with an algebraic description of the wind retrieval problem for a Doppler profiler.This includes an extension/recasting to a frame-based sensing concept.In the second section, an explicit solution for the case of symmetric VAD sampling with constant elevation is provided.This allows a direct calculation of the retrieval error and provides a guideline for an optimal sampling configuration.

Reconstruction of constant wind vector
Under the assumption of a stationary, horizontally homogeneous and vertically piecewise constant wind field, the wind retrieval can be described algebraically.The new aspect is an interpretation of the sensing setup based on the mathematical concept of frames (Christensen, 2008) which allows, for specific sensing scenarios, an explicit computation of involved matrices and therewith an explicit derivation of the associated eigenvalues.

Sensing model and retrieval
For a given azimuth α and zenith angle φ, the beam direction can be described by a unit vector given as The goal is to retrieve an unknown wind vector v ∈ R 3 from projections of v on a set of different beam vectors {e k } N k=1 .This set of beam vectors defines the spatial sampling.The assumption is that within the sampling volume and sampling time the wind vector to be determined v ∈ R 3 is constant, i.e. within the sampling volume we assume a constant wind.
This assumption appears to be overly restrictive, but the goal is not to determine the instantaneous wind vector in an arbitrary turbulent wind field, but rather the mean (horizontal) wind vector over an averaging time of O(10-30 min).For the average wind field, horizontal homogeneity has to be assumed over the area spanned by the beam directions, which is for Doppler profilers typically O(0.1-10 km), and stationarity has to be assumed over the averaging time.In the vertical, the wind field is assumed to be piecewise constant over layers with a thickness of the order of the radial resolution of the Doppler profiler, namely O(10-100 m).
The sampling process can be described through the application of projection matrices.Geometrically the projection of a vector v onto a vector e can be described by the 3×3 matrix P = ee T , which easily follows since the projection of v onto e can be expressed as e, v e = e e, v = ee T v = P v, where •, • denotes the inner product; i.e. for two given vectors By construction, the projection P is a rank one matrix, and, moreover, P is idempotent and symmetric; i.e.P P = P and P T = P .Assume now that the spatial sampling consists of N beam vectors, e 1 , . .., e N , for which we can associate N projections, P 1 , . .., P N .Each beam direction provides us with one radial velocity vector, denoted by p k , k = 1, . .., N ; hence for each k we can write a 3 × 3 linear system, p k = P k v.Note that the magnitude of p k is equal to the (radial) component of the wind field in the beam direction.Combining all N linear systems into one single system results in where p ∈ R 3N and P ∈ R 3N,3 .As each P k is of rank one and as we are usually faced with noisy measurements, directly solving Eq. ( 1) is impossible.A stable retrieval of v can be achieved through a minimization of p − P v 2 with respect to v. The optimal v is given through the solution of the normal equation: A unique solution requires invertibility of P T P in Eq. ( 2), which can be achieved if the rank of P T P equals three.Hence, at least three linearly independent beam directions are (obviously) required to obtain a unique solution.To obtain feasible numerical approximations of v, one has to ensure numerical stability of the inversion process especially in the case of noisy data; i.e. we have to ensure reasonable approximation quality also for the case p = p + with ≤ δ.As we have to solve the normal equation, we first express the symmetric map P T P by its eigensystem, P T P = U DU T , where U is the orthogonal matrix of eigenvectors of P T P and D = diag(λ 1 , λ 2 , λ 3 ) is the diagonal matrix of eigenvalues of P T P .Then, it follows that v = (P T P ) −1 P T p = U D −1 U T P T p. ( With v = (P T P ) −1 P T p , we obtain where λ min denotes the smallest eigenvalue.Therefore, the recovery error can be minimized by maximizing the smallest eigenvalue of P T P .This can be achieved by a proper choice of the corresponding beam vectors e 1 , . .., e N .Hence, the main question to answer is how to set up the beam vectors determining the spatial sampling.

Frame-based recast of the sampling design
In order to answer this question, we consider the set of beam vectors as a frame.Without mathematical rigour, a frame can be seen as a collection of vectors that span the full vector space and which are not necessarily linearly independent.Such a system of vectors is called overcomplete or redundant and allows given vectors to be represented in different ways (non-uniqueness).The redundancy has useful errorsuppressing effects.Using this approach, the goal is to find a simple description of reconstruction stability and reconstruction error dependent on the sampling design.A set of vectors {e k } N k=1 forms a frame for R 3 if there exist constants 0 < A ≤ B < ∞, the so-called frame bounds, such that for all v ∈ R 3 where S is the frame operator introduced below.This frame condition ensures first that all radial components have finite energy and second that the set of beam directions is complete; i.e. there exists no (wind) vector in R 3 that is orthogonal to all beam directions.Let us first reformulate the reconstruction problem.Let e 1 , . .., e N ∈ R 3 denote the individual unit vectors of beam directions, and consider the so-called pre-frame operator T : Then the frame operator defined as S = T T * : R 3 → R 3 is given by which is self-adjoint and symmetric.The frame operator (Eq.6) relates to the above-mentioned projections as follows: = P T 1 P 1 + . . .+ P T N P N = P T P , and thus the invertibility of S is ensured by selecting three linearly independent projections (as already mentioned).In what follows we aim to elaborate how the number of beam directions might change the frame bounds of S, which coincide with the smallest and largest eigenvalues of P T P ; i.e. for the bounds in Eq. ( 5) we have A = λ min and B = λ max .
In order to provide an explicit computation of the solution and therewith an explicit stability analysis, we recast the optimization problem by means of the pre-frame operator T ; i.e. we aim to find an equivalent formulation for p − P v 2 R 3N .First, we have T * = (e 1 , . .., e N ) T : R 3 → R N , and by p k = e k V k = e k (e k ) T v = P k v the normal equation reads as where V ∈ R N is comprised of the radial wind components for the given beam configuration; see e.g.Päschke et al. (2015).This holds true due to G. Teschke and V. Lehmann: Sampling strategies for wind vector retrieval The equivalence of the optimization problems is immediate, and we have a reduction of dimension by a factor of 3.

Explicit solution and error analysis
In practice, the linear system (Eq.3) can be solved numerically through the singular value decomposition (see e.g.Päschke et al., 2015), to minimize errors from finite computational accuracy.This method provides numerical solutions in the general case, and it can therefore be implemented in operational Doppler systems.Nevertheless, an explicit solution of Eq. ( 3) would provide more insight into error propagation and thus allow a further investigation of optimal sampling conditions.Such an explicit solution can indeed be given for a VAD-like sampling scenario.In the following section it is shown that all the involved quantities and error bounds can be explicitly calculated.As these error bounds depend directly on the sensing parameters, the sampling design can be optimized towards a minimal error in the retrieval.

Equispaced circular VAD-like sampling
With pre-assigned equispaced azimuth angles α k = 2π k/N, k = 0, . .., N − 1 and constant zenith angle φ we have The minimization of V − T * v 2 results in Sv = T V or, equivalently, in P T P v = P T p. Hence, in order to provide an explicit expression for the solution of this linear system, we have to derive P T P , which is given by The key for evaluating this matrix is interpreting each of the entries as finite geometric series.Therefore, with the help of the following summations, the matrix P T P simplifies to This means that the frame operator S is diagonal for each φ and N ≥ 3 with frame bounds (see Fig. 2) In the case of A = B, i.e. sin 2 φ = 2cos 2 φ, the frame is called tight.The corresponding φ then satisfies sin 2 φ = 2/3 and hence φ tight = arcsin √ 2/3 ≈ 54.7356.The wind retrieval vector is now easily computed as where with bounds for φ < φ tight : Nsin 2 φ and A −1 = 1 Ncos 2 φ , and for φ = φ tight : A = B = N 3 .With the help of Eqs. ( 15) and ( 16), the explicit algebraic solution is obtained as The matrix multiplications in Eq. ( 17) of the inverse frame operator S −1 with the pre-frame operator T , whose columns are comprised of the unit vectors describing the beams, yield the explicit solution for the wind vector: (18)

Estimation of the retrieval error
Since the wind retrieval for the equispaced VAD sampling case can be explicitly expressed as v = S −1 T V , it is possible to investigate the propagation of measurement errors in the radial wind components to the final wind vector directly.In what follows, the deterministic as well the stochastic error model will be discussed.Assume, as before, the deterministic error model V δ = V + V , where V ≤ δ.For the reconstruction error we then obtain Therefore, with the help of the Cauchy-Schwarz inequality, Consequently, from Eq. ( 20) we deduce The essential observation in Eq. ( 21) is that an increase of the number of beams leads to a smaller reconstruction error and that the smallest error (for any N) is achieved for φ = φ tight ; see Fig. 3. Now assume that the measured radial wind components follow the simple stochastic model, V δ = V + V , with V ∼ N (β, ), where N (β, ) is the N-dimensional normal distribution with expectation vector β and variance matrix .If we assume that the components of β are constant, β i = β for i = 0, . .., N − 1, and = diag(σ 2 , . .., σ 2 ).By computing the expectation of the bias, E( v), one obtains It can clearly be seen that a constant bias in the radial wind estimates affects only the estimation of the vertical wind component, whereas the horizontal wind vector components remain bias-free.This is due to the symmetry of the sampling which leads to a cancellation of any existing bias in the radial winds.
To compute the mean square error (MSE), standard arguments lead to which is obvious as by independency it holds for j = k that E V j V k = E V j • E V k = β 2 , and therefore In the stochastic regime, the deterministic error estimate can be reproduced.Indeed, it can be observed that This estimate verifies the deterministic recovery error, and for growing N this error component can also be made arbitrarily small.The second summand, however, cannot be compensated as it is independent of N.
The MSE can be explicitly calculated as follows: For β = 0 and fixed N, the choice φ = φ tight yields the smallest value for the MSE.The case β = 0 changes the situation.Let , where c = 1 + N β 2 σ 2 .For extremal values, F (φ) = 0 must be evaluated, which is equivalent to evaluating Equation ( 24) provides us for each given N, β and σ with an optimal (MSE-minimizing) zenith distance angle φ.
Finally, from the computation of E v 2 in Eq. ( 23) it follows that supporting and explaining results obtained by Cheong et al. (2008), who have experimentally shown that the MSE or likewise the RMSE of the wind retrieval is significantly reduced by increasing the number of off-vertical beams in the Doppler beam-swinging technique in the presence of wind field inhomogeneities.Note, however, that for the vertical wind component only the random error can be reduced by an increase of N.

Conclusions
In this note, the mathematical concept of frames is applied to the analysis of the spatial (beam configuration) sampling setup for Doppler profilers for the case of a horizontally homogeneous and stationary wind field.It could be shown that it is possible to derive a compact explicit least-squares wind retrieval solution for a typical symmetric VAD scanning scheme.Such an explicit formula had hitherto not been published yet.Besides its simplicity, it allows for a straightforward stability analysis in the practically relevant case of noisy data.The explicit solution exhibits the known fact that the VAD-based estimate for the horizontal wind components is unbiased even if the radial wind components have a constant (direction-independent) bias.Furthermore it was shown that the MSE retrieval error is ∝ 1/N up to a constant offset due to the bias, which means that a larger number of offzenith beam directions is beneficial to reduce the variance of the wind vector components.The total retrieval error is dependent upon the zenith angle φ.For the most relevant case β = 0 it is minimal if the beam vectors form a tight frame.The optimal beam zenith angle for this case was calculated as φ = arcsin √ 2/3 = 54.7356• .It must be noted that the corresponding elevation angle of 35.264 • is a much lower value than what is used in practical configurations of most Doppler systems.However, it is important to appreciate that for the practically relevant case of an estimate of only the horizontal components of the mean wind (since the mean vertical wind component will mainly be close to zero, except for special meteorological conditions) no such optimum can be derived from purely geometrical arguments.Another reason to deviate from this optimal elevation angle is due to the requirement to keep the sampled volume small, in an attempt to minimize deviations from a constant wind field.Technical constraints like the usable Nyquist velocity range and limited scanning capabilities of phased array radar antennas are further reasons why the theoretically optimal elevation is not used in practice.