Measurement of Snow Water Equivalent Using Drone-Mounted Ultra-Wide-Band Radar

: The use of unmanned aerial vehicle (UAV)-mounted radar for obtaining snowpack parameters has seen considerable advances over recent years. However, a robust method of snow density estimation still needs further development. The objective of this work is to develop a method to reliably and remotely estimate snow water equivalent (SWE) using UAV-mounted radar and to perform initial ﬁeld experiments. In this paper, we present an improved scheme for measuring SWE using ultra-wide-band (UWB) (0.7 to 4.5GHz) pseudo-noise radar on a moving UAV, which is based on airborne snow depth and density measurements from the same platform. The scheme involves autofocusing procedures with the frequency–wavenumber (F–K) migration algorithm combined with the Dix equation for layered media in addition to altitude correction of the ﬂying platform. Initial results from ﬁeld experiments show high repeatability (R > 0.92) for depth measurements up to 5.5m, and good agreement with Monte Carlo simulations for the statistical spread of snow density estimates with standard deviation of 0.108g/cm 3 . This paper also outlines needed system improvements to increase the accuracy of a snow density estimator based on an F–K migration technique.


Introduction
There is a vast variety of applications for drone-mounted radar including archaeological investigations; detection of buried mines [1,2]; soil moisture mapping [3], snow, ice, and glacier measurements [4,5]; and mapping of civil infrastructure [6,7]. With regards to snow, one application is snow water equivalent (SWE) measurements, which is of great interest for the hydropower industry and other disciplines in need of meteorological data.
Today, ground-based SWE surveys (manual or with ground penetrating radar (GPR)) are usually conducted by means of snowmobile, where avalanche safety and accessibility might reduce the survey area [8][9][10].
The ground below the snowpack in mountainous and marshland areas often contains sparse scattering objects, potentially producing diffraction hyperbolas in a radar B-scan. These objects are usually rocks with a relative permittivity (∼4-7) that is different from snow. Migration methods applied on radar imagery at the correct propagation velocity of the intermediate medium cause the hyperbolas to collapse at their focal point. Previous studies using commercial GPR, mounted on a snowmobile, show that SWE can be estimated using frequency-wavenumber (F-K) migration and manual velocity picking [8]. A similar method also demonstrates autofocusing using the varimax norm to automatically pick the velocity [11], and similar results can be produced from offset antenna arrays [12]. Furthermore, Kirchhoff's time migration with a two-layered variable-depth velocity model was used to focus radar image GPR data from a helicopter platform [5]. Recent work shows the use of a commercially available, lightweight, and small form-factor 24-GHz radar for This section covers the main method used for SWE estimation and the theory behind the different steps. We assume that the radar image is already preprocessed with both matched filtering and frequency-domain noise filtering. These basic preprocessing steps for the UWB radar data are described in more detail in [25,26].

Altitude Correction
As the UAV performs an airborne survey of an area, the aircraft's altitude will inevitably vary somewhat, at least on a medium to large spatial scale. The altitude variations depend on what sensors the UAV computer has available to feed into the autopilot. In our setup, we have a laser rangefinder [27] that measures relative altitude on a centimeter scale and provides real-time feedback to the autopilot. Nevertheless, minor deviations in the altitude regulation will distort the radar image. Hence, the radar data algorithms need some correction to "level out" these variations.
The method used in this work is to circularly shift the radar data in the fast-time direction according to the relative variations in altitude. The rectification is performed by combining the laser range data and the first surface pulse reflection in the radar return. The combination of altitude information in these two signals reduces spurious deviations significantly and minimizes the influence of signal drop-out in the laser rangefinder.
Prior to migration, the sectioned data segment is expanded in fast-time according to the distance measured from the UAV to the snow surface. This method is explained in more detail in [26].

Methods to Estimate Propagation Velocity from Diffraction Hyperbolas
There are several different methods to estimate the propagation velocity of a radar signal return from an object [28]. Nonetheless, this paper will focus on estimating the propagation velocity by analyzing diffraction hyperbolas in B-scan radar images. In a B-scan GPR image, diffraction hyperbolas are manifested as "south-opening" branches [18]. If diffraction hyperbolas are present in the radar image, there are generally two different procedures to extract the bulk propagation velocity information: • Curve fitting, which attempts to draw a hyperbolic function in the image by optimizing fit with underlying data [18]; • Autofocusing techniques, which use a migration algorithm and generate performance metrics to find the optimum value of the velocity [29].
As previously stated, hyperbolas contain information about the mean propagation velocity in the medium. Hence, if we know the parameters that mathematically describe the hyperbola, we can estimate the propagation velocity [30] by the simple relation between the hyperbola asymptotic constants a and b. Thus, the velocity v rms can be calculated from [31]: where a and b need to be in their correct units, namely, seconds and meters, respectively. Autofocusing is performed by testing different propagation velocities in a migration process where we assess how well the image is focused in order to determine the mean propagation velocity (see example in Section 2.3).

Autofocusing Metrics
Autofocusing techniques are widely used in, for example, synthetic aperture radar (SAR) applications for phase error correction [32][33][34], and for GPR, usually to estimate the dielectric constant [29,[35][36][37]. Autofocusing generally works by testing different propagation velocities v t with a migration algorithm outputting a migrated image s(x, y) v t and choosing the best fitting velocity based on some performance metric.
Since the apex of the hyperbolas should have a maximum at the correct propagation velocity, we can look at the average image intensity AI to evaluate the focusing [29]. For a migrated radar image s(x, y) v t at test velocity v t , this parameter can be stated as where k ∈ [2,4]. The size of the image is represented by m and n as the slow-time and fast-time samples, respectively. However, this metric is known to have poor performance for increasing signal to noise ratio (SNR) [38]. Thus, higher-order techniques that involve the variance of the migrated image [38] can be expressed by where k ≥ 1, andμ andσ are the mean and variance of the migrated data, respectively. As in [38], a k value of 10 was found to be optimal. AH will have a maximum at the best fitting propagation velocity representing the mean propagation velocity v rms along the antenna-to-target trajectory.

Dix's Equation
We emphasize that the best fit velocity v rms represents the average velocity from the antennas to the hyperbola. Hence, to remove the influence of the air section to calculate the average propagation velocity v s in the snow itself, Dix's equation for layered media can be used [39]: where t tot is the total TWT in the medium, t air is the TWT of the air layer, and v air is the approximate propagation velocity in the atmosphere (i.e., approximately 0.2997 m/ns).

Estimation of Snow Parameters
Snow depth is usually measured by evaluating the TWT from the snow surface to the ground and calculating the distance traveled based on an estimate of the propagation velocity.
For lossless, homogeneous, isotropic materials, the relative propagation velocity in the snow v s is related to the relative dielectric constant r by [30]: where c is the propagation velocity in free space. The depth d s is then simply written as where t snow is the TWT from the snow surface to the ground. Equation (5) states that the relative dielectric constant is nonlinearly related to the propagation velocity. From this relation, there are several different models for estimating the density of snow. The relation between snow permittivity and density has been modeled both theoretically [40] and empirically [41]. Considering dry snow, a linear relationship is an acceptable approximation for snow densities below 0.5 g cm −3 , where the relation between the real part of the relative dielectric constant r and relative snow density ρ s can be modeled as [41] where ρ s is relative density compared with the density of water (1 g/cm 3 ). SWE is the depth of water resulting from the mass of melted snow and typically expressed in millimeters of equivalent water [42]. The SWE is the product of the snow depth and the vertically integrated snow density, and can be expressed as Furthermore, SWE is often used to estimate the total precipitation in a given location and the resulting water volume available for the melt season.

Radar System
The ultra-wideband snow sounder (UWiBaSS) is a custom-developed radar system for drone-mounted snow measurements. Papers [25,26,43] detail recent advances of the radar system. New developments include retrofitting the radio frequency (RF) operation band as well as digital modules with 3D printed casings coated in conductive paint to reduce weight while still offering electromagnetic interference (EMI) protection. Additionally, 500-MHz high-pass filters are added to the receiving (RX) channels to reduce low-frequency antenna cross-talk. Table 1 summarizes main features of the UWB-system while Figure 1 shows the UAV "Cryocopter FOX" during flight.

Processing Flow
The basic preprocessing steps for the sampled radar data involve match filtering to correlate the received signal with the transmitted signal, high-pass filtering to remove low-frequency cross-talk between receiving (RX) and transmitting (TX) antennas, reference subtraction, and altitude correction involving rectification of the snow surface return. A detailed description can be found in [26].
The laser altimeter measures the distance to the snow surface. By locating the first interface reflection, the radar is capable of similar distance measurements. We have found that the best practice is to use a combination of both sensors in case of laser signal drop out, which possibly occurs when the laser altimeter beam strikes large snow crystals at an angle.
A distance-vector representing traveled distance over the snow surface is calculated using the position data collected from the UAV autopilot. Using the Haversine formula, we calculate the distance between each coordinate point to generate the distance vector. Due to the non-equidistant sampling of this vector, the radar data goes through a nonlinear interpolation in the slow-time direction. This procedure interpolates the radar data into a domain with equidistant position samples in space, allowing us to fix the sampling interval dx instead of calculating a mean value for each segment.
Hyperbolas in the radar image are segmented manually using a simple "draw rectangle" function, creating a slow-time section of the image to focus within. Furthermore, fast-time data outside the rectangular window is set to zero. Although not implemented here, this step can be automated [18].
The method averages altimeter data within each short section (5 to 20 m). Further, we add zeros to the top of the radar image corresponding to the actual scanning platform altitude above the snow surface, effectively recreating the "real" radar scenario. Figure 2 shows an example image of the processed radar data before the autofocusing procedure.

Autofocusing
The autofocusing procedure involves performing F-K migration for a set of test velocities, where autofocusing metrics are calculated for each test velocity. In order to reduce the numerical calculation time, the search for the optimal velocity can be done in several different ways, such as first performing a low-resolution scan over a broad spread of velocities before performing a high-resolution scan over a much more narrow interval. For the present tests, we used a brute force linear vector in two steps. First a coarse search was performed with step value of 0.01 m/ns to scan through an interval from 0.1 to 0.4 m/ns. Thereafter, a fine search was used with step value of 0.0005 m/ns to scan through an interval from 0.25 to 0.3 m/ns. F-K migration was performed with the CREWES MatLab toolbox written by G. F. Margrave for the CREWES project (University of Calgary) [44,45]. The F-K domain interpolation routine was modified with a convolution, omitting a for-loop in order to reduce the computational time. This modification improved the speed of the function by approximately 200%.
After migration over all test velocities, we selected the autofocusing metric's (see Equation (3)) maximal value, corresponding to a "best fit" propagation velocity. This velocity represents the average velocity from the radar antennas to the ground; therefore, we need to remove the influence of the air section with Dix's equation. This estimation is performed for each manually selected section. The section length is typically in the range of 5 to 20 m containing one or more hyperbolas. Figure 3 shows an example of the autofocusing result for a segment gathered as in Figure 2. Figure 4 shows the combined autofocusing result of the coarse and fine search.

Monte Carlo Simulation of FK-Migration Method
Monte Carlo simulation is a numerical experimentation technique to obtain the statistics of the output variables of a computational model [46]. The proposed method discussed above involves several steps using error-prone parameters. The parameters include altitude and position data from the UAV dataflash-log applied to calculate the distance and the altitude vectors. These vectors position the data in 2D-space and are essential inputs to the migration algorithm.
In this section, we evaluate the uncertainty of the method by creating synthetic data sets of a single hyperbola in a homogeneous medium. The synthetic data is F-K migrated and combined with Equation (3), resulting in an estimate for v rms . Further, using Equations (4) (giving v s ), (5), and (7), we calculate the bulk dielectric constant r and density ρ s of the snowpack.
The altitude and distance vectors are given random errors set to adequate values for our imaging system to make the simulations relevant for real data results. Specifically, the synthetic data parameters are given values similar to the field experiment data in Section 3.1, with a fixed snow-depth of 2 m and an initial altitude (above the snow surface) of 7 m over a horizontal distance of 15 m. The mean propagation velocity is set to 0.29 m/ns, resulting in a snow velocity of 0.258 m/ns after Dix equation analysis.

Laser Altimeter Error Sources
The SF-11 rangefinder data-sheet [27] reports an accuracy of ±10 cm. Additionally, the laser is not always pointing in nadir due to the attitude of the UAV, which is corrected through the attitude data (roll, pitch, yaw). However, given that the UAV often flies over uneven terrain, additional errors are expected.
The rangefinder error is assumed to be normally distributed in the Monte Carlo analysis, with a standard deviation of 15 cm deemed reasonable. The altitude error contributes twice to the erroneous estimations since it is part of two parameters in the method: firstly, prior to the focusing when the length of the fast-time vector is adjusted according to the altitude above the surface; secondly, in the time parameter when Dix's equation is applied.

Distance Error Sources
The distance vector is calculated using the Haversine formula [47] on the filtered position data produced by the UAV autopilot. The position estimate is the output of the autopilot extended Kalman filter (EKF). The EKF estimates vehicle position, velocity, and angular orientation based on rate gyroscopes, accelerometer, compass (magnetometer), global positioning system (GPS), airspeed, and barometric pressure measurements. Concerning the method proposed in this paper, the error can be defined as the deviation from linear movement (constant horizontal velocity), as the absolute position is not relevant for the focusing procedure. In real data, the deviation from a perfectly linear equidistant distance vector (i.e., constant horizontal velocity) was approximately 0.037 m, when looking at real data sets from autonomous flights. Therefore, we choose a standard deviation slightly larger than this (0.045 m) as a conservative assumption for the error analysis.
The distance error is applied to the data by interpolating the synthetic data set along the distance axis according to the linear movement deviation, producing non-equidistant sampling. As expected, the interpolation procedure causes some skewness in the hyperbola that influences the statistical results (see the following subsection).

Monte Carlo Simulation Results
An essential prerequisite of the autofocusing method under study is that it is unbiased to produce a close-to-correct estimate of the sought parameter (snow density or SWE) mean value.
The Anderson-Darling, Jarque-Bera, and Chi-square tests [48][49][50] all indicate that the results from the Monte Carlo simulations are normally distributed both in the bivariate case and with the errors combined, as in Figure 5. As seen in Figure 5, the velocity estimates are unbiased, with a standard deviation of 0.0031 m/ns. Introducing Dix's equation in Figure 5b causes further uncertainties in the estimate since we use the erroneous altitude once more. However, the method is still unbiased, with an increased standard deviation of 0.0147 m/ns. Since the estimation of r in Figure 5c involves squaring of the velocity estimates, the distribution (if assumed normal) is transformed to a Chi-square distribution of order 1. However, if the standard deviation is significantly smaller than the mean, the Chi-square distribution can be approximated as normally distributed. Finally, the approximate linear relationship in Equation (7) is used to obtain the distribution in Figure 5d.
The Monte Carlo simulations emphasize the need for accurate positioning for the presented method as the relative spread of the density parameter is relatively large compared to the mean value (σ ρ /µ ρ = 0.33). Nonetheless, using the average of several estimates to reduce the estimator spread could be a useful approach as the snow density spatial variability is expected to be much smaller than local variations in snow depth [51].

Field Experiments
Due to Covid-19, two major field campaigns were canceled in the winter of 2020. However, a 1-day campaign 10 km from UiT The Arctic University of Norway was carried out. The main goal of the campaign was to test the latest iteration of the radar system and preliminarily investigate methods to measure snow density.
A transect of approximately 250 m was flown autonomously for four passes to test the reproducibility of the method. A section of the transect showed several overlapping hyperbolas at the ground level (shown in Figure 2). These hyperbolas were assumed to be caused by reflections from rocks, as is well known from traditional GPR [30], but this is not necessarily the case given that the UWiBaSS would be able to detect hyperbolas from targets buried beneath several meters of snow during flight. As shown in Figures 6 and 7, the snow depth in this transect varied from 1 m to almost 6 m.  After the spring snowmelt, radar data images containing hyperbolas were confirmed to be from a rocky area, as shown in Figure 8.
Four passes over this rocky area of the transect were segmented into 40 smaller segments of varying size (5 to 20 m) used to estimate the density.
For this data set, the returning signal from the snowpack was taken at the second period of the transmitted signal owing to the radar system having a 5.9-m unambiguous range, and the data collection was taken from approximately 7-m relative altitude, excluding snow depth. For more information on performing measurements outside the unambiguous range from UAV radar, see [26].
The in situ density for each noticeable layer was collected in a snow pit close to the transect. We found the weighted density average (weighted mean based on the thickness of each layer) to be 0.327 g cm −3 . Figure 9 shows the snow pit density profile along with hardness and grain type. Grain size was 1 to 2 mm.   Figure 10 shows the distribution of density estimations from the autofocusing procedure. Figure 6 shows measured snow depth along the 250-m transect for four passes. Identification of the snow-ground interface is described in [26], where we obtained the TWT used in Equation (6). Snow depth is then calculated using the TWT and the estimated propagation velocity of the snow.

Discussion
The accuracy of the method seems to depend heavily on the accuracy of the recorded UAV position, since the spatial sampling rate and correct position and altitude above snow are important for the F-K migration algorithm. Hence, real-time kinematic (RTK) GPS should be used for better distance calculations and data interpolation [52]. Taking the average of the altitude for each section will, in most cases, not introduce significant errors, as the typical variation in altitude across the sections is on the order of a few cm. However, for robustness, modifying the code to create a variable altitude above the snow should be implemented in future work [5].
Ideally, several snow pits should have been dug along the transect. However, the time constraints of the project limited the in situ collection to one full snow pit.
Comparing Figure 5d with Figure 10, we see that the statistical spread (standard deviation) is of the same order, which leads us to believe that the error assessments made to the scanning system prior to the Monte Carlo simulation were reasonable and that other contributing error sources are of minor importance in the campaign.
The standard deviation (0.1 g/cm 3 from Figure 10) of the estimations are comparable to what is observed in previous papers [8], reporting a standard deviation of 0.07 g/cm 3 from a snowmobile platform with a fixed antenna-to-snow air gap of approximately 50 cm.
The degree to which the snow depth measurements can be trusted is of vital importance when determining the SWE. A comparison of depth measurements for several passes along the same transect gives a good indication of the repeatability of depth measurements. This can be seen for four passes in Figure 6 and for two passes overlaid on the radar B-scan in Overall, the correlation is very high, with the lowest value of 0.92 between passes 2 and 4. The slight skewness between the depth measurements in Figure 6 appears to come from drift in the UAV positioning. Hence, the repeatability of snow depth measurements will likely also benefit from more accurate positioning such as RTK GPS. The snow depth is measured beyond 5 m, which is deeper than the 3-m depth probes we had available.
However, in situ snow depth validation of the same radar system has been presented in [26].
Notice in Figure 2 that the reconstructed radar pulse delay becomes almost 90 ns, resulting in a total distance from the antennas to the ground of approximately 10 to 14 m. With an unambiguous range of 5.9 m, this is accomplished by exploiting the lack of significant reflections in the air-section except for the cross-talk [26] and, thus, using the previously transmitted pulse for detection.
The area coverage of this UAV-mounted radar system is small compared to satellite remote sensing products [24]. However, the high spatial resolution generated by the UAVmounted radar has the opportunity to detect local variability, resulting in more accurate SWE measurements.

Conclusions
In this paper, we present a noninvasive method to estimate snow density from a UAV that could be further used to estimate SWE. Density estimations from a limited field trial show good agreement with in situ snow density, and snow depth measurements show high repeatability. However, the sample size for both in situ and radar data is somewhat limited.
Future work will include field campaigns with RTK-ready UAV for improved position and altitude information, as well as high-resolution in situ density measurements using the SnowMicroPen (SMP) or similar devices. Additionally, automatic image segmentation should be implemented to reduce the amount of manual labor to analyze the data sets. Multithreaded programming should be implemented on the interpolation part of the F-K migration routine to improve speed, as this section of the code occupies approximately 95% of the run time.  Acknowledgments: The authors would like to thank M. Eckerstorfer and A. Kjellstrup for collecting in situ data and operating the UAV during the field experiments. The authors acknowledge the useful comments by four anonymous reviewers that helped improve the readability of this study.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Mapping of Squared Inverse Gaussian Variable
Consider a Gaussian variable X with standard deviation σ and mean value µ. Then, defining χ 2 is centrally chi-square-distributed with one degree of freedom.
Knowing that E[χ 2 ] = 1 and E[χ 4 ] = 2 in a central chi-squared distribution, the first and second moments of X can be stated: where the approximations are valid if σ µ-that is, a narrow probability distribution of X.
Mapping a Gaussian variable X with the relation w = RX 2 , we obtain the probability density function (pdf) [53]: Notice that, for situations where σ µ, the last term in Equation (A4) can be neglected unless w and µ are close to zero.
Then, Equation (A4) can be stated: Furthermore, Taylor series expanding the argument in the exponential function to second order around the central value w 0 = µ 2 R of the distribution, and similarly approximating the denominator √ Rw ≈ µR, we obtain p(w) ≈ 1 2µRσ We recognize Equation (A6) as a Gaussian distribution with standard deviation σ = 2µRσ and mean value µ = µ 2 R, in accordance with the results obtained in Equations (A2) and (A3), if R = 1.
The permittivity of snow is related to propagation velocity v through = (c/v) 2 , where c is the speed of light in air. In order to relate this mapping to Equation (A6), define the reciprocal variable u = 1/w. Now, it can readily be shown that Once again, Taylor expanding u 2 to zeroth order and Taylor expanding the argument of the exponential function to second order, both around u 0 = 1/µ , Equation (A7) reduces to p(u) = 1 √ 2π(σ /µ 2 ) exp − (u − 1/µ ) 2 2(σ /µ 2 ) 2 . (A8) From Equation (A8), we see that u is also normally distributed with a mean value µ u = 1/µ and standard deviation σ u = σ /µ 2 .
Furthermore, the constant R can be identified as R = 1/c 2 and X = v, leading to the final result that is close to normally distributed with first and second moments: provided that σ µ. According to [41], there exists a linear relation between the dielectric constant and density ρ in dry snow: This leads to the first and second moments of ρ: Hence, we obtain a constant change in mean value and spread in the probability distribution of 1/2 as expected from a linear transformation.
As an example, assume a statistical draw of v from a Gaussian distribution. Figure A1 shows the histogram from 100,000 statistical realizations mapped by using = (c/v) 2 with first and second moments of snow velocity of 0.244 m/nsecs and 0.0127 m/nsecs, respectively. Observe that the fit between the histogram and the theoretical approximate Gaussian distribution is not perfect, as the assumption that σ µ (0.16 1.51) is not entirely fulfilled.