WindSeer: Real-time volumetric wind prediction over complex terrain aboard a small UAV

WindSeer is a convolutional neural network (CNN) with four-channel input. The input is composed of a binary mask indicating cells containing input measurements, the terrain model stored as a distance field, and the sparse horizontal wind speed measurements (two channels). Note that vertical wind speed measurements are not an input to the model, since weather stations typically measure only the horizontal wind and the vertical wind is not observable with a standard fixed-wing sUAV sensor set. In Appendix C we show empirically that adding vertical wind as an input, if it were available, has a limited impact on prediction quality. The percentage of observed cells in the input data varies across experiments ranging from \qty3.5e-6 to \qty0.19, thus WindSeer always operates on highly sparse observations.

The four-channel output of WindSeer has the same spatial dimension and resolution as the input. The first three channels contain the three-dimensional wind prediction ($W_{x},W_{y},W_{z}$) and the fourth channel contains the turbulence kinetic energy (TKE) prediction — a metric for the strength of the turbulent velocity fluctuations in the wind field that is proportional to the sum of the variances in each dimension. Thus, the prediction contains properties ($W_{z}$, TKE) that are not available as input measurements.

We evaluated WindSeer on held-out CFD-simulated flows. The dense label data allowed for evaluation over the full domain, thus evaluating the influence of different measurement locations as well as qualitatively characterizing the prediction quality.

We evaluated WindSeer on real wind data available from three published measurement campaigns. The long measurement periods allow filtering out short-term effects such as wind gusts and reduces the overall measurement noise, enabling an evaluation of WindSeer with clean input wind measurements over larger-scale domains.

We evaluated WindSeer on noisy, local wind measurements collected by multiple fixed-wing sUAV. This mirrors the targeted use case of WindSeer and challenges it with high input noise levels due to real sensor noise and short term wind effects such as gusts or turbulence.

We evaluated prediction times of WindSeer on an NVIDIA Jetson Orin AGX, a light-weight and low-power single-board computer, to show real-time performance on sUAV flight-grade hardware. The average inference times over 100 runs on a $64^{3}$ and $384\times 384\times 192$ prediction domain were \qty[separate-uncertainty-units = single, separate-uncertainty=true]0.021(2) and \qty[separate-uncertainty-units = single, separate-uncertainty=true]3.577(15) respectively. Mixed-precision inference reduced the inference times to \qty[separate-uncertainty-units = single, separate-uncertainty=true]0.021(5) and \qty[separate-uncertainty-units = single, separate-uncertainty=true]1.700(5). These inference times show that WindSeer is capable of low-latency wind predictions over large domains with limited compute to quickly recalculate predictions in response to new measurements.

We developed a pipeline (Fig. 1) to train and deploy \@iaciCNN CNN, WindSeer, that predicts the dense wind and turbulence around complex terrain. The network training consists of two steps: First we generated a dataset of dense flows over terrain patches from Switzerland using \@iaciRANS RANS CFD solver (Fig. 1 A)). We then trained WindSeer using the label flows to simulate local wind measurements along randomly generated piecewise-linear trajectories, robustifying the predictions by adding noise to the measurements along the trajectories (Fig. 1 B)). The trained WindSeer was evaluated on (i) held back CFD-simulated flows on previously unobserved terrains, (ii) real wind data gathered in measurement campaigns [9, 10, 12, 13, 14], and (iii) real wind data measured by multiple sUAV around mountainous terrain.

We generated flow data over real terrain patches with a pipeline based on the open source solver OpenFOAM [50, 51] with the steady-state RANS model and the popular $k-\epsilon$ two-equation turbulence closure [52]. The automated pipeline ingests terrain patches and outputs the time-averaged flow solutions for multiple wind speeds [53]. We extracted 563 terrain patches each with an extent of \qtyproduct1.5 x 1.5\kilo from the GeoVite service, which provides access to the swissALTI3D digital elevation map (DEM) for Swiss researchers, with a lateral resolution of \qty0.5m²²2https://geovite.ethz.ch/, recently also available on ArcGIS:
https://elevation.arcgis.com/arcgis/rest/services/WorldElevation/Terrain/ImageServer. The terrain patches exhibit at least one side with near-constant elevation allowing us to simulate a formed boundary layer flow (logarithmic profile) entering into the domain from that face. Some terrains allowed for multiple flow directions leading to 866 terrain/flow direction pairs. The vertical extent of the simulation domain was three times the height difference of the terrain with a lower bound of \qty1100 minimizing the boundary effects on the flow. Each case was simulated with up to 15 different wind speeds if the automatic meshing succeeded, resulting in 7361 executed CFD runs. We initialized the subsequent simulation for the higher wind speed cases with the previous solution to speed up computation. Only solutions that met a required optimization tolerance were accepted as fully converged solutions, which was the case in \qty92.9 of the runs. We enhanced our dataset with one zero-velocity flow for each terrain that had at least one converged CFD simulation, resulting in a total of 7285 flows.

The CFD solutions are computed on an automatically generated irregular mesh with OpenFOAM’s SnappyHexMesh utility. We resampled each case up to a height of \qty1100 to a regular \numproduct91 x 91 x 96 grid resulting in a resolution of \qty16.5 horizontally and \qty11.5 vertically (Fig. 1 A)).

Generating CFD flows is a computationally and labor-intensive task. For reference, our 7361 CFD runs required \qty9168 CPU compute time (\qty782 creating the meshes and \qty8386 solving the flow, average compute time: \qty1.25). Unfortunately, deep networks are notoriously data-hungry and, for a complex modeling problem such as wind prediction, would typically require orders of magnitude more training data to achieve good performance. In computer vision, image augmentation methods are widely used when training deep CNN [54]. These methods aim to improve the quality and size of the datasets when only limited data is available to prevent the networks from over-fitting. In this work, we showed, for the first time, that geometric transformations can be applied to CFD flows to augment the WindSeer training data.

We randomize the locations of terrain features and flow directions by generating $64^{3}$ subdomains sampled from each full \numproduct91 x 91 x 96 grid. The subdomains are constructed by sampling from a range of rotations and origin translation offsets inside the full domain. In a first step the horizontal shift and a rotation around the $z$-axis are sampled from bounded uniform distributions to ensure that the shifted and rotated $64^{2}$ subdomain is fully contained within the full $91^{2}$ domain. Then in a second step the vertical shift is sampled from a triangle distribution with lower limit and mode of 0 and an upper limit of 32. Smaller vertical offsets are favoured to focus on the complex flow regions closer to the terrain. The flow data is linearly interpolated to the coordinates of the subdomain grid, which is the same spatial resolution as the full grid.

The input to WindSeer consists of four volumetric channels, one of which corresponds to the terrain encoding $T$ created as a Euclidean distance transform with zeros inside the terrain. That representation propagates the terrain information over the full domain and even allows us to include terrain features outside of the domain if they are accounted for in the distance field calculation. The remaining three channels include the sparse and noisy horizontal wind measurements ($U_{x,in},U_{y,in}$) and a binary mask $B$ indicating cells containing measurements. Previous work has shown the value of providing binary input masks to CNN handling sparse input data [55, 56].

Measurements from weather stations or realistic flight scenarios only cover a small percentage of the prediction volume along a connected path, e.g. a \qty30 flight segment with our sUAV covers approximately 20 cells at the default grid resolution. Consequently, for a practical onboard wind prediction scenario, we expect the available input wind data to be very sparse and thus construct our network input to reflect this sparsity. We create the input based on the augmented dense flow by creating a mask and then selecting the measurements based on the mask. We emulate the characteristics of \@iacisUAV sUAV flight path by filling the mask along sequential randomly-selected piecewise linear segments with a length of 3 to 500 cells.

Noise is added to the sampled wind data in order to account for fluctuations in the wind and sensor errors that are not captured by the RANS CFD simulations. Two types of disturbance are added, white Gaussian noise (sampled i.i.d. at each measurement from $\mathcal{N}\left(0,\sigma_{g}^{2}\right)$) and measurement bias (sampled from $\mathcal{U}\left(-0.1,0.1\right)$ and applied to all measurements). The first has the purpose of simulating noise due to sensor measurements [57], while the latter simulates the effects of sensor miscalibration. The standard deviation for the Gaussian noise $\sigma_{g}$ itself is drawn from a uniform distribution: $\sigma_{g}\sim\mathcal{U}\left(0,0.1\right)$, simulating different noise levels. All the noise values are scaled with the mean wind velocity for each sample in the training set to have coherent noise levels from low to high velocity samples. Note that noise is only added to the training inputs and not to the CFD ground truth labels used to compute the network training losses.

The sparse input implies that for most cells in the input wind velocity channels ($U_{x,in},U_{y,in}$) the values are undefined since they do not contain a measurement. We test and evaluate two approaches to filling the missing information. The first naïve approach simply places zeros in all voxels without a measurement. This results in large gradients of the input for high magnitude wind. The second approach uses the per-channel-average of all measurements as the fill value resulting in a smoother input and propagating the information over the whole domain.

The labels are constructed by stacking the four volumetric channels corresponding to the three-dimensional predicted velocity ($U_{x,out},U_{y,out},U_{z,out}$) as well as the TKE at each cell from the CFD ground truth flows.

The wind prediction model is an encoder-decoder CNN with skip connections based on the U-Net architecture [58]. The WindSeer encoder is composed of single 3D convolutions with kernel size 3 and reflection padding to preserve the size. Using skip connections, the information at each depth is relayed to the decoder before utilizing a max-pooling layer with kernel size of 2 to down-sample the feature map. The original domain size is restored by pairing the information from the skip connection with a nearest-neighbor up-sampled feature map followed by two 3D convolutions with kernel size 4. This decoder structure removes checkerboard artifacts sometimes experienced when using a decoder-encoder CNN [59]. Each convolution, except the final one, is followed by the nonlinear ReLu layer with negative slope 0.1 [60].

A majority of the cells in the wind speed input channels contain no measurements. Our first approach was to set the values of all these cells to zero. However, this resulted in a network overfitting to the number of observed cells and did not generalize to larger domains and different resolutions as demonstrated in Appendix C. In the end, we ended up filling the unobserved cells with the average of all measurements per channel which results in a smoother input and helps to propagate the information over the full domain.

A scaled version of the mean squared error (MSE) loss is applied to train the model balancing the loss $L(\cdot)$ between the samples and channels:

where $X$ is the network prediction, $Y$ the label flow, $\hat{Y}_{c}$ the label average per channel of the non-terrain cells, and $N$ the number of non-terrain cells. Normalizing the error by the average label value balances the loss for flows of different magnitudes. Without accounting for the number of terrain cells in the loss, a sample with a high ratio of terrain cells would not contribute much to the overall loss. Thus, scaling according to $N$ prevents these cases from being underrepresented in the training.

The model is trained using the Adam optimizer [61] for 3000 epochs. The initial learning rate of $1.0\times 10^{-5}$ is quartered every 700 epochs.

Each of the three measurement campaign datasets that we used for evaluation are publicly available but require some preprocessing to enable direct comparison with our wind prediction outputs. We convert the data from the different file formats for each measurement campaign to the same gridded format that we use to store the CFD solutions. Each experiment provides terrain data as well as wind measurements collected using static masts equipped with airflow sensors at various heights. The terrain is discretized by querying the raw data using bilinear interpolation in the center of the respective cell. The location of each measurement is converted into the cell coordinates.

WindSeer predicts the wind using the measurements from one mast. The measurements are filled into the corresponding cell and averaged in case of multiple measurements in one cell. The predictions, which are obtained with trilinear interpolation at the sensing locations, are then compared to the measured wind. CNN allow for variable input sizes, a trait we exploit to predict at a higher spatial resolution for domains with smaller length scales (see Bolund Hill below) at an increased domain size of $384\times 384\times 192$ cells. Since the terrain is represented as a Euclidean distance field, this gives WindSeer a sense of the grid resolution and thus the scale of the flow, enabling us to predict the wind at different scales.

The error bars for Fig. 5 are calculated using error propagation from the standard deviations of each axis if not reported for the magnitude [62].

We ran the inference time experiments on an Orin AGX⁷⁷7https://www.nvidia.com/en-us/autonomous-machines/embedded-systems/jetson-orin/, a low power, light weight (\qty623 including the carrier board and heatsink) and small scale (\qtyproduct105 x 105 x 60\milli) single-board computer that can be carried by a small scale sUAV. We set up the Orin AGX with the Jetpack 5.1 software kit that includes CUDA 11.4 and cuDNN 8.6.0 and installed PyTorch 2.0. During the evaluation we ran the Orin in the maximum power mode (\qty60) using all 12 CPU cores.

We used three Multiplex EasyGlider4 airframes equipped with the Pixhawk 4 autopilot [63] using the high quality ADIS16448 inertial measurement unit (IMU) and the u-blox M9N GNSS module for autonomous navigation. We configured the main height source of our modified PX4 autopilot [64] to the GPS height and use the barometric pressure as a fallback. An extension to the guidance law adjusting the airspeed ensured safety during strong wind conditions [5]. We used a custom designed pitot tube with the Sensirion SDP31 differential pressure sensor and Hall sensor airflow vanes to enable measuring the 3D wind vector. Refer to Section D for more details about the airflow sensing setup and calibration procedure. We used a ground station computer with QGroundControl to control and navigate the sUAV. While the default PX4 state estimator could be extended to estimate the 3D wind we opted for an offline flight path reconstruction (FPR) pipeline using an iterated extended Kalman filter (see a similar problem definition in [65]). The offline FPR pipeline allowed us to generate high quality estimates for validating our approach and to adjust the estimation pipeline post flight.

We gathered wind data from flights at three test sites in Switzerland. The first test site at Chasseral is one of the most topographically isolated mountains in Switzerland and is located in the Jura mountains (47° 07’ 38” N, 7° 02’ 47” E, \qty1548 above mean sea level (AMSL)). The other test sites are located on the ridges of the Oberalppass (46° 39’ 24” N, 8° 40’ 21” E, \qty2069 AMSL) and Gotthardpass (46° 34’ 17” N, 8° 33’ 33” E, \qty1960 AMSL) in the Central Swiss Alps. These were chosen to evaluate the prediction performance for domains surrounded by complex terrain. The spatial constraints allowed for two sUAV to simultaneously collect wind data at Oberalppass and Gotthardpass and three sUAVs at Chasseral. The sUAV were flown simultaneously in circular loiter patterns with a radius of \qty100m leading to lateral separation between the planes of up to \qty800m and measurements in different flow regimes. We planned the flights based on NWP forecasts ensuring good flight (no precipitation, fog or clouds) and stable wind conditions (wind magnitude below the cruise speed of \qty10\per, direction and magnitude near-constant over multiple hours).

We use two modes to convert the raw wind estimates to the WindSeer input. In the first mode we generated the input by averaging the measurements over one loiter pattern to generate a single wind estimate. This time-averaged data, similar to the averaged data in the measurement campaigns with static masts, helped to reduce the noise and sensitivity to sensor calibration on the input measurements. The observations from the different sUAV and multiple loiters enabled us to compute the average prediction error and the correlation for each flight experiment to see if the flow trends are well predicted. In the second mode the input is composed of the wind data from a \qty120 window of one sUAV to predict the wind along the flight path within the next \qty120 window for both the input (itself) and the validation sUAV. This sequential time-windowed setup allowed us to qualitatively evaluate the WindSeer performance along the flight paths. A high-resolution elevation map provided by SwissTopo [31] was used to construct the terrain for the WindSeer input.

The authors declare that they have no competing interests.

The synthetic CFD training data, as well as the sUAV flight measurements will be made available through our dataset server: https://projects.asl.ethz.ch/datasets/doku.php?id=nmi_23_windseer.

The code is available on github: https://github.com/ethz-asl/WindSeer