The MFAS algorithm

The MFAS algorithm for continuous-wave lidar Doppler sensing

Wind speed and direction can be found from a set of Doppler lidar measurements, even if each single measurement provides only single-axis (line-of-sight) information. One “wind vector algorithm”, popular from the early days of both CW and pulsed scanning lidars, fits a sine wave to a series of Doppler frequency estimates obtained at different view angles and (usually) different times [1]. For pulsed lidars another method, using a “maximum function of the accumulated spectrum” (MFAS for short), was examined by Banakh and Smalikho [2], [3], and a variant was recommended by Augère et al. [4]. Weighted and/or iterative fitting routines and probabilistic algorithms need increased computing power, and are often trial-and-error or proprietary, but there are recent relatively explicit examples for sine-fitting [5]. In this note we illustrate MFAS for both wind estimation and target Doppler signal estimation, for CW lidars.

Wind vector estimation

We assume for the moment that the wind is horizontal and homogeneous and described by a vector V. Its component VLOS along the lidar line-of-sight (LOS) is the scalar

VLOS = |V| sin(φ) cos(θi-θ)

where the lidar cone angle is φ, the lidar azimuth scan angle is θi, and the wind angle is θ.

Suppose that lidar “spectra” (accumulations of many consecutive short-term |DFT|2 outputs) are obtained at different scan angles θi (i = 1, 2, 3,…) with fixed cone angle φ [6]. In each spectrum we expect a peak according to the Doppler relation fi = 2VLOS/λ, where the peak’s frequency fi depends on scan angle:

|V| = λfi/ [2 |sin(φ) cos(θi-θ)|]

With V unknown in practice, the MFAS method starts with a guess or hypothetical wind vector Vguess. This represents a point in a 2D grid (Vx, Vy), or equivalently a magnitude (Vx2 + Vy2)1/2 and a direction θguess = tan-1(Vy, Vx).

Let us assume that the lidar is not direction-sensitive (ZephIR is an example): although VLOS may be positive or negative, the measured frequencies are non-negative.

If θguess = θ (i.e. the guess is correct in angle), then a spectral peak is expected at a frequency

(2/λ)|V||sin(φ) cos(θi-θ)|

and the spectral peaks for different scan angles should appear at a common frequency if each spectrum is rescaled by dividing all the bin frequencies by |cos(θi-θ)|. The spectral peaks will not all coincide if the angle is wrongly guessed. Thus, if for each different θguess in our grid we “stack” or accumulate the rescaled spectra by summing over all θi, we expect the peak bin of the stack to be highest for a narrow range of θguess near the true value θ.

A peak is expected at some frequency of the stack whenever θ is correctly guessed, but not at the true Doppler frequency unless |Vguess|is also correct. So we do not need to seek for each guess (Vx, Vy) the maximum of the entire accumulated stack (that is, all frequency bins, each containing the sum of the powers from all the rescaled spectra). Instead we find for each guess (Vx, Vy) the sum of powers in only the narrow region near the expected peak. In this sense the name “maximum function of accumulated spectra” needs qualification.

DFTs have discrete frequency points (bins), but |V| and the scan angle are in effect continuous variables, so the exact expected frequency (2/λ)|V||sin(φ)cos(θi-θ)| will generally fall between two DFT points, and some approximation or interpolation will be involved. We may choose, according to judgement and experience, only one DFT point, or a value interpolated (as in Banakh/Smalikho) between the nearest two, or some other smoothed estimate. We may iterate the MFAS algorithm, finding an improved estimate by searching through a smaller more densely sampled grid around the current estimate.

Interpolation of the entire frequency scale is convenient for showing complete spectrograms. But, given the assumptions of the algorithm, including the existence of a reasonably homogeneous atmosphere in steady horizontal flow, we need in practice to store, scale, and accumulate the spectra in only a small region (at most a few DFT points or bins) at or near (2/λ)|Vguess||sin(φ)|, as said above.

Note the implications if we think of MFAS as compensating for φ and θ to recover an estimate of what the Doppler spectrum would look like with the LOS parallel with the wind vector. For example, is the wind flow “homogeneous”, and to what extent is the observed “Doppler” spectrum determined by non-homogeneous movement (e.g. turbulence) and other factors such as modulation? The expected or hypothetical shape of the Doppler spectrum depends on our answers and may or may not be given by a simple rescaling.

If we “stretch” the frequency scale by the factor 1/|cos(θi-θ)|, the apparent power densities are altered, as they would be if we observed a moving mass of air (with significant Doppler spread) from varying angles. We can then conserve the total power in the measurement bandwidth (another point that is sometimes not explicit in the literature) by dividing the bin heights by the same factor.

There are other possible approaches, such as a uniform shift of all frequencies by an amount equal to the shift expected for the wind Doppler peak. Note that, with “circular” DFT conventions, frequencies shifted beyond one edge of the band set by the sampling frequency reappear at the other edge.

A “best” approach cannot be decided without further information, and several aspects must be balanced:

- The scaling becomes unrealistic at angles θi nearly normal to the wind (θ ± π/2), although in practice this “blow up” may not matter because the default processing excludes a small range of near-DC Doppler where instrumental and detector noise can dominate any desired signal.

- We may expect the spectrum of the lidar atmospheric return to depend on angle, but we do not usually expect the spectrum of detector shot noise to depend on angle; yet both are likely to be significant at the medium or low signal-to-noise for which MFAS is recommended.

- The apparent spectral width of the lidar atmospheric return, even if it can be well characterised by a single number (such as the width of a Gaussian “spectral bump”), is a result of several filtering and convolution processes. These include the isolation of short (~5 μs) segments of time series for the |DFT|2 routine – which creates the bin width (DFT point spacing) of ~ 200 kHz, for any scan angle and no matter how wide or narrow the range of frequencies in the detector output [6].

Thus some factors or features in the presented “spectra” may depend strongly on scan angle, and others may not, and a single rescaling choice will not treat them all equally.

In our MATLAB examples each “spectrum” Si, meaning the set of frequency bin contents obtained by summing many consecutive short-time |DFT)|2 at a fixed (or nearly fixed) scan angle θi, is “rescaled” by dividing all the bin frequencies by |cos(θi-θ)|. The rescaled spectrum is interpolated by MATLAB’s pchip function, yielding new frequency points from DC up to some reasonable maximum wind speed; typically we choose 100 equally spaced points from DC to 20 ms-1. Several of these rescaled spectra, typically 10-40 covering a substantial fraction of one complete azimuthal scan (50 scans/second), are stacked.

Below is an example of MFAS applied to CW lidar data from airport trials, similar to those in [7] and chosen to show a fairly clear Doppler track with moderately high carrier-to-noise and without severe distortions caused by e.g. wake turbulence.

Figure 1: Spectrogram of CW scanning lidar data.

Fifteen consecutive spectral estimates (20 msec each) are plotted together in Figure 1 in a “spectrogram” time-frequency format. The track is indeed “fairly clear” to the eye, and we note again that these estimates emerge from considerable preprocessing [6], including hardware filtering, digitising, whitening, and discarding of information in the lowest bin(s) near DC. It seems likely that the LOS component peaks around estimates 9/10/11, with a magnitude of about 25 bins (about one fifth of the DFT range of 128 bins covering 0-15 ms-1). Given a cone angle near 30 degrees [6], [7] we have sin(φ) ~ ½ and expect that the true wind magnitude for this example was about 6 ms-1.

A broad search through possible wind vectors shows:

Figure 2: Power in MFAS accumulated |DFT)|2.

This has the two peaks expected for a lidar that (as we assumed) does not detect the sign of Doppler shifts. The best-fit wind vector cannot be distinguished, without more information, from the vector with the same magnitude but opposite sense.

A fine search around one peak suggests a “best” vector near (5.3 ms-1, 2.6 ms-1) with magnitude near 6.0 ms-1:

Figure 3: As Figure 2, with a finer grid around one of the two peaks.

For this chosen wind vector, the same time-frequency format shows the scaled spectral estimates aligned so that, as intended, most of the Doppler peaks fall within a narrow range of DFT bins:

Figure 4: Spectrogram of Figure 1’s CW lidar data after rescaling according to the MFAS estimate of wind vector.

Smalikho’s two-bin interpolation [2] omits both the stretching and the pchip interpolation, and his plots of accumulated power are less smooth. The precise fine detail will depend on our pre-processing choices:

Figure 5: Two-bin MFAS interpolation for the same Figure 1 data.

Augère et al. [4] recommend “a rough MFAS estimate in order to initiate a WVML [wake vortex maximum likelihood] algorithm dedicated to wind accuracy”.

Note also that the underlying timescales can, in principle, be adjusted. The usual ZephIR processing [6] accumulates about 4000 DFTs during the 20 ms for each spectral estimate; each estimate yields a mean Doppler estimate, or (in MFAS) a contribution to the “stack”, labelled with a single instant of time. Thus ZephIR’s continuous (not stop-stare) scan introduces a small timing mismatch; its effect on the accuracy of ordinary ZephIR wind estimation is very small; nonetheless it can be reduced in MFAS, if we wish, since smaller groups or even the individual small DFTs (200000 per second) can be individually rescaled with their own accurate θi, so that each contributes to the “stack” and is labelled with its own time instant (with a negligible error of a few μs).

This finer resolution in scan angle and time requires more memory and processing, and increases the bin size (of the DFTs and their stack), but incurs no extra penalty from nonlinearities in an estimation algorithm. In contrast, the usual mean-Doppler algorithms applied to estimated spectra become strongly nonlinear as the carrier-to-noise falls below some threshold, so we often accumulate over substantial times before attempting the estimation.

Alternatively, we can search for the wind vector by estimating first the best-fit angle and then the best-fit magnitude. In Figure 6 the same 15 scans are used and the bin contents are summed for various hypothetical angles. The peak occurs near 0.46 rad. Then a search through possible magnitudes shows that the best fit to the observed data of Figure 1 occurs for |V| ~ 6 ms-1; thus the same result as above is obtained.

Figure 6: Result of MFAS search for best-fit wind angle.

Our initial assumptions were very simple, and in practice we expect inhomogeneities and vertical wind components.

Accumulated spectra in target Doppler sensing

In lidar Doppler sensing, a two-step method is common [8]: the first step involves a time series of short-time frequency estimates, and the second a power spectral estimate (or spectrogram etc.) of that series. Sharp spectral features (tonals, near-sinusoids) are then identified with “signals” or vibrations of the target. As with MFAS for scanning lidars, it is possible to estimate the parameters of a sinusoidal frequency variation without performing the intermediate step of Doppler frequency estimation: that is, we form a series of short-time spectral estimates – with no explicit judgement of their peak or mean Doppler values – and rescale/realign them in a search for the “best” value of a vector. The vector now represents the amplitude and phase of a periodic movement rather than the wind magnitude and direction.

In a way this is merely curve-fitting or best-fit peak-searching, without the above issues of spectral densities and LOS dependence; but note that individual nonlinear estimates of peak Doppler are unnecessary, because the “best” or most plausible signal can be determined in a single stage (resembling MFAS or maximum likelihood etc.) rather than in two stages.

The simulated example below has an FM signal Acos(2πft + ϕ) with modulation frequency f = 1 Hz, frequency deviation A = 5000 Hz, and initial phase ϕ = π. The sampling frequency is 50,000 s-1 and multiplicative Gaussian speckle of bandwidth 5000 Hz is present [8]. A three-parameter search (for A, f, and ϕ) correctly identifies the signal; in contrast with the scanning-lidar example above we do not assume in advance that the sinusoid frequency f is 1 Hz.

Figure 7: (left) Spectrogram of simulated carrier; (right) Frequency modulation signal estimated from MFAS-type algorithm.

References

[1] K A Browning and R Wexler, “The determination of kinematic properties of a wind field using Doppler radar”, J. Appl. Meteor., vol. 7, pp. 105–113, 1968.

[2] I Smalikho, “Techniques of wind vector estimation from data measured with a scanning coherent Doppler lidar”, J. Atmos. Oceanic Technology, vol. 20, pp. 276-291, 2003.

[3] V Banakh and I Smalikho, Coherent Doppler wind lidars in a turbulent atmosphere, Artech House, Boston / London, 2013.

[4] B Augère, M Valla, A Durécu, A Dolfi-Bouteyre, D Goular, F Gustave, C Planchat, D Fleury, T Huet and C Besson, “Three-dimensional measurements with the fibered airborne coherent Doppler wind lidar LIVE”, Atmosphere, vol. 10, 2019.

[5] Xunbao Rui, Pan Guo, He Chen, Siying Chen, and Yinchao Zhang, “Adaptive iteratively reweighted sine wave fitting method for rapid wind vector estimation of pulsed coherent Doppler lidar”, Optics Express, vol. 27, pp. 21319-21334, 2019.

[6] C A Hill and M Harris, “Lidar measurement report” in Technical Report for UpWIND Work Package 6; Technical University of Denmark: Roskilde, Denmark, 2010.

[7] C A Hill, J Bennett and D Smith, “Airport trials with the Aviation ZephIR lidar”, Proc. 15th Coherent Laser Radar Conference, Toulouse, France, 22-26 June 2009.

[8] C A Hill, M Harris, K D Ridley, E Jakeman and P Lutzmann, “Lidar frequency modulation vibrometry in the presence of speckle”, Applied Optics, vol. 42, pp. 1091-1100, 2003.

(September 2020)