III.1

2D (azimuth, range) cross-correlations.

This approach computes 2-D cross-correlations in the (azimuth, range) domain as previously done in Ref. 1,. The kernel corresponds to near-instantaneous⁴ sections of features by the scan cone, and the goal is to find this section again at a later (or previous) time. In practice, there are two cases:

Figure 4.

Two cases for 2D cross-correlations in the (azimuth, range) domain: “across” mode (left) and “tangent” mode (right). An idealized aerosol feature (ellipse) is depicted crossing the lidar scan domain (dotted circle) at two different times t (purple) and t + Δt (green). Its section is imaged by the lidar at each time (bold lines), these are what is being correlated.

Figure 5.

Kernel and search region for the 2D (azimuth, range) cross-correlation (left). The 3D (scan, azimuth, range) data volume is in practice seen as a 2D (time, range) plane to simplify computations (right).

From the 3D data point of view, the kernel is defined as a subset of a 2D (azimuth, range) slice, and the search region extends around it in all three dimensions (Figure5, left). However, with full 360° scans the data is continuous in the azimuth dimension from one scan to the next. To reflect this without having to resort to complicated boundary condition, the 3D data volume is “unfolded” as a 2D (time, range) plane (Figure 5, right).

III.2

“Beam-to-beam” 1D (scan) and 2D (scan, range) cross-correlations

The beam-to-beam approach considers the time series observed at each azimuth position – that is to say, in the scan dimension of the 3-D volume.

For a given beam, the kernel is selected at a given azimuth angle in, either the scan dimension at fixed range (1D variant), or in the (scan, range) dimensions (2D variant). The 1D (scan) or 2D (scan, range) cross-correlations are computed between this kernel and a search region for every other beam — Figure 6. The overall peak corresponds to the best match.

Figure 6.

Illustration of beam-to-beam correlations for the 1D case (scan dimension). The lidar senses at fixed azimuth positions (blue dots). For a given azimuth (purple dot), the kernel (gray window) is extracted from the time-series of backscatter values (purple curve) recorded at this position. It is correlated with time-series from every other azimuth. The best match is here represented in green.

IV.1

Input data

This is a collection of backscatter profiles along with their metadata (headers). The metadata should either directly provide, or contain all necessary information to compute, each data point space-time coordinates:

IV.2

Outputs Data

Each estimation provides the following outputs:

Optional outputs can be useful for visualizations.

V.1

Processor Functions Implementation

Fig. 9 below shows the realised Electronic Hardware (Fig. 9a) with embedded Processor Functions (Fig. 9b)

Fig. 9a:

Electronic Board

Fig. 9b:

Processor Functions on the Raspberry Pi

V.2

Pre-processing

Raw data from the GS are processed before running the motion estimation algorithms. The steps are standard in lidar data processing; they are applied to each profile independently:

Coherent aerosol features are required to run the motion estimation. The usual signal-to-noise ratio (SNR) is not the best indicator, a good SNR is not necessarily related to the presence of features (in particular in the near range). Authors in [3] proposed the iSNR which estimates a local “coherent texture”-to-noise ratio; regions of the scan where the iSNR is below a given threshold (set to 3 in the article) are not considered for the motion estimation process.

Figure 9 presents examples of 2D (time, range) data recorded by the GS at azimuth angles 0° and 180° along with the corresponding binary masks of iSNR > 3. The iSNR is here computed as in [3], in practice it will need to be tuned according to the data recorded by the GS after the upgrade modifications.

Figure 10.

The upgraded CWL system with scanner & scanning geometry algorithm-optimised

VI.1

Wind measurements with scanner-optimized Ground-based System

The tests done on synthetic data sets have enabled to validate the “scan-range 2-D cross-correlation” method as the most robust algorithm for motion estimation by using the new processing-optimized scanning system & sampling strategy, Figs. 11 & 12

Figure 11.

The figure on the left shows a Photek 24-Channel array PMT-MCP array detector and its electronics

Figure 12

shows the 16 Channel Eventech single-photon-counting board & pre-processing electronics

The choice of the altitudes of interest z_min and z_max, of the desired observable velocity bounds U_min, U_max and W_max, as well as hardware constraints will lead to setting the best parameters for both the GS (elevation angle, scan rate …) and the motion estimation algorithm (in particular, sizes of kernel and search regions).

VI.2

Algorithm validation on synthetic data

Synthetic datasets enable to quantitatively evaluate the various motion estimation methods. The objective here is not to physically simulate the remote lidar sensing of aerosols, but rather to study the results of the algorithms assuming a perfect signal (coherent features everywhere, no noise …) and for a given sampling configuration. The aerosol distribution is simulated by a 3D texture (x, t) transported by a constant velocity (u(x, t) = u):

Physical and lidar acquisition processes such as turbulent diffusion, integration time, gate width, range attenuation, instrument noise … are ignored. The texture is here based on a multiscale Perlin noise. It is a coherent procedural texture of fractural nature mathematically defined at every point of space as shown in Figure 13.

Figure 13:

example of synthetic scalar field from multiscale Perlin noise.

A first series of numerical experiments on synthetic data, not reported here, led to select the “scan-range 2-D crosscorrelation” method (section III.2) as the most robust algorithm for motion estimation with CWL-like scanning configuration.

Additional numerical experiments were conducted to further validate the wind estimation algorithm prior to their use with the “real” CWL observational data. The objective was to check the motion estimation module with controlled “back and forth” scan data (reciprocating scan pattern), using a scan configuration similar to that used for the 2022-07-06 CWL data acquisition:

Scan configuration:

Wind Configuration:

Synthetic data Configuration:

Examples of generated synthetic scans and associated estimated motions are presented in Figures 14 and 15. Estimated wind vectors are coherent with the ground truth, except for the South-West quadrant showing spurious estimates. This was expected, as this quadrant is a “blind spot” for the given scan settings (0° to 270°) and wind direction (180°): aerosol features are seen entering, but never leaving the scan cone. In operational context, such spurious values would be discarded by Quality Control methods.

Figure 14:

Synthetic data for the “reciprocating” scan pattern – time-series at various azimuth positions.

Figure 15:

Motions from synthetic reciprocating scans of Figure 13. Left: estimated wind vectors at each azimuth position of the scan at range 1500 m, axes units are in meter in the horizontal plane. Spurious estimates are visible in the lower-left (South-West) sector. Right: time-series of estimated wind speeds at different ranges, horizontal axis is the time (s) and vertical axis the wind speed (m/s). The ground truth is a 5 m/s West wind.

These experiments confirmed that the proposed motion estimation algorithm is able to handle this CWL scanning configuration.

VI.3

CWL Wind Data Verification with Real wind data

The algorithm was then applied to CWL data collected during the night of 2022-07-06 to 2022-07-07. This dataset features a dense cloud layer at range about 1500 to 2000 m (Figures 16 and 17).

Figure 16:

CWL data of 2022-07-07 showing a dense cloud layer at about 2000 m range (“20220706” dataset, file indices above 14000). Dashed white lines denote scan bounds.

Figure 17:

CWL data of 2022-07-07 (same as Figure 16), here presented as time series at various azimuths positions, focusing on the cloud layer.

The proposed cross-correlation used the following dimensions for the correlation kernel: 21 points in the “scan” (also time) domain (about 22 minutes); 31 points in the range domain (about 310 m).

Results gathered over about 10 min of operations are presented Figure 18; only the vectors for which the correlation peak is above 0.7 are retained.

Figure 18:

Raw wind analysis for CWL data of 2022-07-07 (“20220706” dataset, file indices above 15000 – see also Figures 16 and 17) at about 1500 metres range. Left: wind vectors estimated at various azimuth positions. Right: time-series of wind speeds (top) and direction (bottom). Bottom: visualization on the underlying scan data, dashed lines correspond to the beginning of a back-and-forth scan, black markers denote the (time, range) locations of estimates shown in the time-series.

Wind directions appear to gather around 150°, while most speed estimates lye between 8 to 12 m/s. Values may seem scattered, but the following points must be kept in mind:

The winds calculated on 07.07.22, @ 1:40 am at the Park Cottage Observatory, Salehurst location, East Sussex (South London) by the newly developed processing algorithm with the scanner optimized CWL have been verified with winds reported at the same time & location from the US Nat. Met. forecast data (https://earth.nullschool.net) and from the Lydd airport station from Wunderground (https://www.wunderground.com/history/daily/gb/romney-marsh/EGMD/date/2022-7-7). Both sources fit very well with our calculated intensity (8-14 m/sec) and direction (150 degree).