GNSS-IR Snow Depth Retrieval Based on the Fusion of Multi-Satellite SNR Data by the BP Neural Network

Abstract

Compared with previous snow depth monitoring methods, global navigation satellite system-interferometric reflectometry (GNSS-IR) technology has the advantage of obtaining continuous daily observation data, and has great application potential. However, since GNSS satellites are in motion, their position in the sky is constantly varying, and the Fresnel reflection areas about the receiver in different periods alter accordingly. As a result, the retrieving results obtained from different GNSS satellites, and data sets collected in different periods, fluctuate considerably, making the traditional single-satellite-based GNSS-IR retrieving method have limitations in accuracy and reliability. Therefore, this paper proposed a novel GNSS-IR signal-to-noise ratio (SNR) retrieving snow depth method for fusing the available GNSS-IR observations to obtain an accurate and reliable result. We established the retrieval model based on the backpropagation algorithm, which makes full use of the back propagation (BP) neural network’s self-learning and self-adaptive capability to exploit the degree of contribution of different satellites to the final results. Then, the SNR observations of the global positioning system (GPS) L1 carrier from the Plate Boundary Observation (PBO) site P351 were collected to experiment for validation purposes. For all available GPS L1 carrier data, the snow depth values retrieved for each satellite were first obtained by the existing single-satellite-based GNSS-IR retrieval method. Then, four groups of comparison results were acquired, based on the multiple linear regression model, random forest model, mean fusion model, and the proposed BP neural network model, respectively. Taking the snow depth in-situ data provided by snow telemetry (SNOTEL) as a reference, the root mean squared error (RMSE) and mean absolute error (MAE) of the proposed solution are 0.0297 m and 0.0219 m, respectively. Furthermore, the retrieving results are highly consistent with the measured data, and the correlation coefficient is 0.9407.

1. Introduction

As an essential component of the hydrological system, snow plays an important role in balancing freshwater resources, and in the process of climate regulation [1]. Therefore, it is important to continuously and accurately monitor the snow depth. The traditional way of snow monitoring mainly relies on the manual measurement on the ground, or the deployment of snow depth detection machines (using ultrasound or laser) [2], though the labor cost and time-consumption are appreciable. With the development of the GNSS, the GNSS-IR remote sensing technique has gradually become a research hotspot, and has been widely applied in sea surface altimetry, sea surface wind field, soil moisture, and snow depth monitoring [3,4,5,6,7,8,9].

In snow monitoring, Larson et al. [10] first proposed the theory of snow depth retrieval using GPS SNR data, and experimentally demonstrated the consistency between GPS snow depth estimates, snow depth in-situ data, and snow detector measurements. Later, Larson et al. [11] selected several stations in the PBO network to illustrate the conditions for the successful application of the technique. Ozeki et al. [12] proposed a non-geometric distance L4 phase method based on GPS satellites to invert the snow depth, and obtained a similar accuracy as the SNR retrieval method. Tabibi et al. [13] used the SNR data of GPS L5 for snow depth retrieval, and showed that L5 has almost the same performance as L2C. Zhang et al. [14] further improved the accuracy of snow depth retrieval by establishing a grid model for retrieval of snow depth to attenuate the influence of surrounding terrain, compared with the results provided by the PBO H2O team. To improve the retrieval accuracy, Zhang et al. [15] used a dynamic clustering algorithm to filter the power spectrum density (PSD) of Lomb–Scargle spectrum (LSP) results, and then processed the filtered results with Grubbs’ criterion. Wang et al. [16] analyzed the differences of multiple signals from multiple GNSS systems for snow depth retrieval, and argued that the differences in the environment around the station affect the performance of satellite retrieval at different azimuths. Deng et al. [17] introduced the wavelet decomposition method to process the SNR sequence to reduce the influence of noise signals on the results. Comparative analysis with the LSP, fast Fourier transform (FFT), and nonlinear least square fitting (NLSF) algorithms, conducted by Li et al. [18], found that the best corresponding-to-different-snow-depth algorithms are also different, and proposed a combined NLSF+FFT algorithm. Wang et al. [19] conducted a detailed discussion of the observation–geometric conditions and their influence on the snow depth measurements retrieved with SNR data from GPS and BeiDou. The factors include elevation angle range, arc segment length, number of satellites, and azimuth angle.

Previous GNSS-IR snow depth retrieval methods can be broadly classified into two categories, i.e., retrieval approaches based on a single satellite, and retrieval methods based on a partial satellite under a specific azimuth. Both categories may implement relatively accurate snow monitoring, but there are two apparent disadvantages. For the former, many experiments are required to determine the best retrieval satellite, since the retrieval performances vary corresponding to different satellites [20,21]. For the latter, the exclusion of a large amount of valid satellite data makes it difficult to utilize the daily observations fully [11,22,23].

Therefore, this paper proposes a GNSS-IR multi-satellite data fusion snow depth retrieval model based on the BP neural network. The optimal weight of each satellite is obtained through the BP neural network, so that all of the observations of each satellite that meet the requirements can be used. It solves the problem that a large amount of observation data cannot be utilized in existing algorithms. In addition, three other fusion methods are introduced and compared with the measured snow depth data provided by SNOTEL to verify the effectiveness and accuracy of the proposed method.

2. Materials and Methods

2.1. GNSS-IR Snow Depth Retrieval Principle

For GPS positioning measurement, the signal received by the receiver is a vector sum of the direct signal from the satellite and the reflected signal that enters the receiver through the reflection of the features around the station [24]. The strength of the synthesized signal can be expressed in terms of SNR, which is mainly affected by various factors, such as the transmitting power of the satellite signal, antenna gain, and multipath effect [25]. Under low elevation angle conditions, the multipath effect is stronger, and a large number of reflected signals carrying surface feature information enter the receiver, resulting in lower positioning accuracy and lower SNR quality. Conversely, the surface parameters that cause the multipath effect can also be obtained by processing the SNR data. Therefore, the SNR data before and after snow cover can be processed to obtain the distance from the antenna phase center to the reflecting surface under the conditions of no snow and with snow, and the difference between the two is the snow depth. Figure 1 shows the schematic diagram of GNSS-IR retrieval of snow depth.

The relationship between the direct signal, reflected signal, and SNR can be expressed as [11].

$$SNR\propto P_d+P_r+\sqrt{P_d{P}_r}\mathit{cos}\varphi$$

(1)

where $P_d$ is the direct signal power, $P_r$ is the reflected signal power, and $\varphi$ is the reflected phase. Previous studies have shown that the antenna gain mode determines the $P_d\gg P_r$, which means that the direct signal determines the overall trend of the SNR, and the reflected signal causes a slight effect on the SNR.

The direct signal is usually obtained using a low-order polynomial fit [26], which is eliminated to obtain a reflected signal carrying surface information. In this paper, a quadratic polynomial is selected to fit the SNR, and the following formula is used to linearize it before fitting:

$${\mathrm{SNR}}_{\left(\mathrm{volts}/\mathrm{volts}\right)}={10}^{\left(\frac{\mathrm{SNR}\left(\mathrm{dB}-\mathrm{Hz}\right)}{20}\right)}$$

(2)

Detrended SNR signals can be expressed [11] as the following equation:

$$\mathrm{SNR}-P_d-P_r=A\mathit{cos}\left(4\pi h{\lambda }^{-1}\mathit{sin}e+\phi \right)$$

(3)

where $A$ denotes the average value of factor $\sqrt{P_d{P}_r}$ over the arc span, $h$ is the vertical distance from the phase center of the antenna to the ground, $\lambda$ represents the carrier wavelength of GPS, and $e$ is the satellite elevation angle.

The SNR sequence with the trend term removed is a function of the satellite elevation angle, and the residual sequence has a constant frequency after the antenna height is determined. For snow depth retrieval, the distance from the antenna phase center to the snow surface is unknown, so LSP analysis is imported to obtain the frequency at the peak of the amplitude [27]. Then, the height from the antenna phase center to the snow surface is calculated using the following equation:

$$h=\frac{\lambda f}{2}$$

(4)

Obtaining the antenna height $H$ when there is no snow, and the reflection height $h$ after snow accumulation; then, the snow depth $h_{snow}$ can be expressed as:

$$h_{snow}=H-h$$

(5)

It is worth noting that not all of the observations are suitable for estimating snow depth. A sufficient amount of data is required to extract the frequency of multipath modulation accurately. Therefore, the satellite needs to rise to a specific elevation (set to 20° in this paper). However, the limited number of geodetic receiver channels, and being set to preferentially track satellites with high elevation angles, resulted in data below 10° not being recorded sometimes [11]. Figure 2 shows the multipath reflection point of some satellites of the P351 station on DOY 10, 2015; the dashed circles from the inside to the outside show the track of the reflection point for satellite elevation of 20°, 15°, 10° and 5°, respectively. It is clear that some tracks do not extend to 20°, e.g., the track for PRN 3, PRN 12, and PRN15 in the south do not extend to 5°, as well as the track for PRN 3 and PRN 12 in the northeast. All of the observations at these times are indicated by dashed lines, and are excluded in the subsequent processing.

2.2. BP Neural Network

BP neural network is a multilayer feedforward network trained by error back propagation. Such a network consists of a series of nodes, and includes two processes of forward propagation and error back propagation. The basic idea is to adjust the weights among the neurons by iteration so that the error mean square difference between the actual output value and the desired output value of the network is minimized. Therefore, it has good self-learning, adaptive, robustness, and generalization capabilities [28]. The structure of the BP neural network is shown in Figure 3, which consists of the input layer, hidden layer, and output layer. There can be more than one hidden layer in a neural network; the neurons in the same layer are independent of each other; and the neurons between adjacent layers are completely connected [29]. The establishment process of the neural network mainly contains two parts, i.e., signal forward propagation and error back propagation. The principle of the BP neural network will be briefly introduced in the following (Figure 3).

Let the input vector be $X=\left[\begin{array}{cccc}{x}_1& x_2& ...& x_i\end{array}\right]$, the output of each neuron in the hidden layer be $H=\left[\begin{array}{cccc}{h}_1& h_2& ...& h_j\end{array}\right]$, the output vector be $Y=\left[\begin{array}{cccc}{y}_1& y_2& ...& y_k\end{array}\right]$, and $W_1, B_1$ denote the weights and biases from the input layer to the hidden layer; then, the activation of each neuron in the hidden layer is:

$$h_j=f\left({\displaystyle \sum }_i{w}_{ij}{x}_i+b_j\right)$$

(6)

$$W_1=\left[\begin{array}{cccc}{w}_{11}& w_{12}& ...& w_{1j}\\ w_{21}& w_{22}& ...& w_{2j}\\ ...& ...& ...& ...\\ w_{i1}& w_{i2}& ...& w_{ij}\end{array}\right],B_1=\left[\begin{array}{c}{b}_1^1\\ b_2^1\\ ...\\ b_j^1\end{array}\right]$$

(7)

where $f$ is the activation function of the neuron.

The process from the hidden layer to the output layer is the same as the process described above, except for the difference in weights and biases. The forward propagation is the process of the input vector reaching the output layer from the input layer through the hidden layer to get the computational result of the neural network. If m samples exist, the error function can be defined as:

$$E=\frac{1}{2m}{\displaystyle \sum }_{i=1}^m{\displaystyle \sum }_{k=1}^n{\left(dk\left(i\right)-yk\left(i\right)\right)}^2$$

(8)

where $y_k\left(i\right)$ is the expected value of the corresponding $x_i$ output through the hidden layer, and $d_k\left(i\right)$ is the reference value of the corresponding $x_i$. In general, the result obtained by one forward propagation process does not achieve the desired result; at which time, the neural network enters the error back propagation process. Then, the algorithm will iteratively calculate to adjust the weights and bias of each neuron according to the following two equations until the allowable error is reached or the set number of training is reached.

$$w_{ij}^{\left(l\right)}=w_{ij}^{\left(l\right)}-\alpha \frac{\partial E}{\partial w_{ij}^{\left(l\right)}}$$

(9)

$$b_i^{\left(l\right)}=b_i^{\left(l\right)}-\alpha \frac{\partial E}{\partial b_i^{\left(l\right)}}$$

(10)

where $l$ is the layer number, $\alpha$ is the learning rate, and its value range is (0, 1).

Nodes are the core of the hidden layer, which are used to extract and store the intrinsic connections in the samples. If the number of nodes is too small, it is difficult for the BP network to extract the intrinsic patterns of the samples, which, in turn, leads to a decrease in prediction performance and low fault tolerance. If the number of nodes is too large, it will lead to a long learning time of the network, and may overfit and reduce the generalization ability. There is no strict mathematical formula to calculate the number of nodes, and an empirical formula is usually applied in determining the number of nodes [30].

$$L<\sqrt{m+n}+a$$

(11)

where $L$ is the number of nodes in the implicit layer, $m$ is the number of nodes in the input layer, $n$ is the number of nodes in the output layer, and $a$ is a constant within 1 to 10.

2.3. Data Source

The GPS data used in the prediction of snow depth came from the P351 station of the Plate Boundary Observation (PBO), located in the Ketchum area of Idaho, USA, at an average elevation of 2692.6 m. The station is surrounded by open and flat areas, with few obstacles affecting the GPS signal. There are almost no obstacles such as buildings and trees nearby that affect the GPS signal, and the snow period is up to more than 200 days per year, which is ideal for snow depth retrieval. The antenna used at station P351 is TRM29659.00 SCIT, and the receiver is TRIMBLE NETRS. The sampling time for this experimental data is from 1 January 2015 (DOY1) to April 10, 2015 (DOY100); the receiver sampling rate was 15s; and the SNR data type was GPS L1 C/A. Figure 4 shows the location of P351.

Snow depth in-situ data were obtained from the SNOTEL, operated and managed by the USDA Resource Conservation Service (NRCS), with over 900 stations widely distributed throughout the western United States. The SNOTEL station, numbered 490 (43.87497° N, 114.71363° W), is the closest station to P351 at approximately 450 m, and provides snow depth data for the day.

3. Experiment and Results

3.1. Experimental Technical Scheme

Figure 5 shows the flow chart of snow depth retrieval using multi-satellite data in this paper. The whole technical process can be divided into two parts. (1) Satellite observation files and navigation files were processed using the GNSS data pre-processing software TEQC, developed by Dr. Lou Estey to obtain SNR and elevation angle (ELE) data for each satellite, followed by processing according to the steps described in Section 2.1 to obtain the snow depth retrieved by all satellites that meet the conditions independently. (2) The data of DOY 1-DOY 80 were used to establish a BP neural network model, whereas the data of DOY 81-DOY110 were used for model verification and comparative analysis with the multiple linear regression model, random forest model, and mean fusion model, respectively.

3.2. Snow Depth Extraction

The SNR of satellites can be extracted from the observation file using TEQC software, and Figure 6 shows the SNR data of a rising arc for satellite pseudo-random noise (PRN) 7 on day of the year (DOY) 4, 2015. As the satellite elevation angle rises, the antenna gain is enhanced, and the multipath effect is reduced [21], the SNR rises from about 32 dB to near 55 dB.

The SNR data were linearized; then, quadratic polynomial fitting was used to obtain the direct signal part, and then removed from the SNR to acquire the reflected signal. Figure 7 shows the processing of SNR data.

Analysis SNR data with the trend term removed using LSP were used to obtain frequency at amplitude peak. In fact, the snow depth will have a great impact on the LSP calculation results. In order to reveal this phenomenon, this paper calculates LSP results for 3 days with different snow depths. The red curve, green curve, and blue curve in Figure 8 represent the LSP on DOY 40, 65, and 110, 2015, respectively. These curves are distributed over the three months of gradually warming weather when the snow depth gradually decreases and the dominant frequency increases. In addition, the peak of each curve is much larger than the background noise, indicating that the specular reflection is strong at this time, which is favorable for snow depth retrieval.

Under the condition that the antenna height is known, the snow depth is calculated from Equations (4) and (5).

3.3. Multi-Satellite Data Fusion

The BP neural network is a model with high self-learning and self-adaptive ability, so it is especially suitable for solving problems with complex and unclear internal mechanisms. In addition, it has certain fault tolerance. Damage to local nodes does not significantly impact the global results. This advantage makes the BP neural network very suitable for snow depth retrieval based on multi-satellite data.

(1) Model Structure

Cybenko demonstrated that a neural network with a single hidden layer could approximate any function with a sufficient number of hidden nodes. Therefore, only one hidden layer was set up for the experiments. The input layer is the snow depth of 25 satellites satisfying the conditions. The output layer is the multi-satellite fusion result. According to the formula mentioned above and the actual test, the number of neurons in the hidden layer is determined to be six. Besides nodes, the function is another factor to be considered, which contains transfer functions, training functions, and learning functions. The choice of function has an important impact on the results. However, there is a lack of systematic research on how to choose the best combination of functions; so, in this paper, we choose to test different functions to determine, and finally determine the transfer function as tansig, the training function as traingda, and the learning function as learngdm based on the experimental results.

(2) Network training

For training and testing the network, the snow depth retrieval result of each satellite of DOY 1-DOY 80, and the snow depth data provided by SNOTEL, on the corresponding dates were input into the network, and divided into three parts, i.e., training set, validation set, and test set with the ratio of 70%, 15%, and 15%, respectively. Figure 9 shows the performance images of the trained neural network. It can be seen from the figure that the best performance of the neural network is reached when the training reaches the 96th generation.

(3) Network validation

The snow depth retrieval result of each satellite of DOY 81-DOY110 is input into the trained BP neural network to obtain the snow depth values of the multi-satellite retrieval. This paper introduced multiple linear regression, random forest, and mean fusion methods for comparative analysis to prove the feasibility and accuracy of the proposed model and solution.

3.4. Experiment Results

3.4.1. Single-Satellite Snow Depth Retrieval Results

Due to the different environments around the stations, and the varying performance of the satellites themselves, when the snow depth alters, the satellites may respond differently. Figure 10 shows the discrepancy between the snow depth in some satellite retrievals of DOY 81-DOY 110 and the snow depth provided by SNOTEL. The retrieval quality of different satellites varies greatly, and the specific results calculated for each satellite are given in

Table 1. RMSE, MAE, and R between single-satellite snow depth retrieval and SNOTEL snow depth.

PRN	RMSE/m	MAE/m	R
2	0.1484	0.1335	0.4579
15	0.1022	0.0906	0.6385
21	0.0622	0.00500	0.7078
27	0.0514	0.0355	0.8220

. From Figure 10 and Table 1, it can be found that the deviation between the retrieval results of PRN 2 and the measured values is noticeable. The deviation on DOY 83 is up to 0.349 m, the correlation coefficient is low, and the retrieval quality is poor also. The quality of the retrieval results of PRN 27 is greatly improved. The correlation with the measured results increases while the RMSE decreases. Even so, it still deviates from the measured values several times. It can be concluded that when using the SNR for snow depth retrieval, it is hard to guarantee the accuracy and reliability of the results by a single satellite.

3.4.2. Multi-Satellite Snow Depth Retrieval Results

Figure 11 shows the retrieval results of the multi-satellite snow depth retrieval model established by the above methods. Overall, the results of all four methods are consistent with the measured snow depth. However, it is evident that t random forest and multiple linear regression results reflect apparent deviation from the in-situ data, and the results of mean fusion and BP neural network retrieval are closer to the measured values.

4. Discussion

As shown in Figure 12, the snow depth values retrieved by all four methods show a reasonably consistent trend with the in-situ data provided by SNOTEL. At the same time, the excessive deviation between them is considerably reduced. The result indicates that these methods all respond well to the changes in snow depth, and overcome the limitation of the single-satellite retrieval. However, there are still shortcomings in these algorithms. Above all, the retrieval results by random forest have significant deviations from the measured values, and the results in the last days of the experiment deviate significantly from the actual values (maximum deviation up to 0.17 m). Besides, since the least squares algorithm is used in the regression, which is known to treat all data equally, when there is a large deviation between the single-satellite retrieved data and the in-situ data of the day, squaring the residual will lead to a stronger impact on the whole, making the regression line offset, so the results of the multiple linear regression jumped several times. For example, for DOY100-DOY105, the MLR predictions show large fluctuations, although the measured snow depths remain relatively stable. In addition, the retrieval results of mean fusion and the BP neural network have significant improvements, and are closer to the in-situ data, and the maximum deviation is approximately 0.08 m. However, because each satellite has the same weight, the satellite with poor retrieval performance is not conducive to error suppression for the mean fusion method. As a result, the deviation to the in-situ data is more significant than that of the BP neural network.

To further evaluate the performance of different methods, this paper introduces three indexes, i.e., correlation coefficient (R), RMSE, and MAE. Figure 13 shows the correlation between the retrieval results of the different methods and the in-situ data provided by SNOTEL.

Table 2. R, RMSE, and MAE between different multi-satellite snow depth retrieval methods and SNOTEL snow depth.

Method	R	RMSE/m	MAE/m
Random Forest	0.8710	0.0669	0.0542
MLR	0.8479	0.0498	0.0372
Mean Fusion	0.9276	0.0376	0.0329
BP Neural Network	0.9407	0.0297	0.0219

shows the specific values of each accuracy metric for the different methods.

As shown in Figure 13 and Table 2, the results corresponding to the four retrieval methods strongly correlate with the measured snow depth provided by SNOTEL. The BP neural network retrieval acquires the best result, with a correlation coefficient of 0.9407. The mean fusion, random forest, and multiple linear regression models have slightly inferior results, with correlation coefficients of 0.9276, 0.8710, and 0.8479, respectively. The RMSE and the MAE with the measured values are all between 0.02–0.07 m, which significantly improved compared with the single-satellite retrieval model. The BP neural network has the best performance, with the smallest values of RMSE of 0.0297 m and MAE of 0.0219 m, respectively. The improvement of RMSE and MAE are at least 20% and 30%, respectively, compared with the other three methods.

5. Conclusions

Long-term and accurate snow depth measurement is a significant reference to the water cycle and climate regulation research. For acquiring continuous daily observations, this paper proposes a GNSS-IR SNR solution for retrieving snow depth based on the fusion of available satellite data by the BP neural network. Subsequently, the SNR observations of the GPS L1 carrier from the PBO site P351 were collected for snow depth retrieval. We completed four experiments for validation purposes, based on the multiple linear regression model, random forest model, mean fusion model, and the proposed BP neural network model, respectively. Furthermore, the comparative analysis was conducted with the in-situ data provided by SNOTEL. The following conclusions can be drawn through relevant statistic data and analysis:

• The snow depth retrieval results according to data sets collected in different periods during the experiment showed great fluctuations. Therefore, it is hard to guarantee the accuracy of snow depth retrieval results by using data obtained from a single GNSS satellite.
• The retrieval results have reached an accuracy of 7 cm, and the correlation has also been greatly improved, which is better than the traditional single-satellite retrieval method.
• Compared with related models, i.e., multiple linear regression model, random forest model, and mean fusion model, the retrieving accuracy and reliability of the BP neural network are optimal. The statistical analysis indicated that the correlation between the retrieval results and the in-situ data reaches 0.94, and the RMSE is less than 3 cm.

The experiment results show that the proposed GNSS-IR SNR solution has better accuracy and reliability than previous methods for snow depth retrieval. However, the relationship between the retrieval quality of SNR data under different azimuths at the same station and surface environmental factors, such as topography and vegetation cover, is unclear, and their potential influence on the algorithm deserves further study in follow-ups. Considering there are four different global navigation systems nowadays, and the available SNR data have increased significantly, making full use of observations from multiple systems to achieve more accurate and reliable snow depth retrieval will be the focus of future research. It should be pointed out that the proposed method is available for obtaining daily snow depth observation data in some local areas. Furthermore, the proposed method and related data acquired are expected to provide data support and new opportunities for a further breakthrough in related research.