GNSS-IR Snow Depth Retrieval from Multi-GNSS and Multi-Frequency Data

: Global navigation satellite system interferometric reﬂectometry (GNSS-IR) represents an extra method to detect snow depth for climate research and water cycle managing. However, using a single frequency of GNSS-IR for snow depth retrieval is often found to be challenging when attempting to achieve a high spatial and temporal sensitivity. To evaluate both the capability of the GNSS-IR snow depth retrieved by the multi-GNSS system and multi-frequency from signal-to-noise ratio (SNR) data, the accuracy of snow depth retrieval by different frequency signals from the multi-GNSS system is analyzed, and a joint retrieval is carried out by combining the multi-GNSS system retrieval results. The SNR data of the global positioning system (GPS), global orbit navigation satellite system (GLONASS), Galileo satellite navigation system (Galileo), and BeiDou navigation satellite system (BDS) from the P387 station of the U.S. Plate Boundary Observatory (PBO) are analyzed. A Lomb–Scargle periodogram (LSP) spectrum analysis is used to compare the difference in reﬂector height between the snow-free and snow surfaces in order to retrieve the snow depth, which is compared with the PBO snow depth. First, the different frequency retrieval results of the multi-GNSS system are analyzed. Then, the retrieval accuracy of the different GNSS systems is analyzed through multi-frequency mean fusion. Finally, the joint retrieval accuracy of the multi-GNSS system is analyzed through mean fusion. The experimental shows that the retrieval results of different frequencies of the multi-GNSS system have a strong correlation with the PBO snow depth, and that the accuracy is better than 10 cm. The multi-frequency mean fusion of different GNSS systems can effectively improve the retrieval accuracy, which is better than 7 cm. The joint retrieval accuracy of the multi-GNSS system is further improved, with a correlation coefﬁcient (R) between the retrieval snow depth and the PBO snow depth of 0.99, and the accuracy is better than 3 cm. Therefore, using multi-GNSS and multi-frequency data to retrieve the snow depth has a good accuracy and feasibility.


Introduction
Snow is an important part of the land hydrological cycle and global climate system. Accurate real-time snow depth data are an important reference indicator for water resource management and climate disaster warning [1,2]. Among the existing snow depth detection methods, in situ snow sensor measurement lacks time resolution, and global navigation satellite system interferometric reflectometry (GNSS-IR), as a new microwave sensing technology, has proven to be able to realize snow depth detection [3]. GNSS-IR is a kind of satellite remote sensing technology that uses GNSS signals as the transmitting source to realize the retrieval of the physical parameters of surface targets by receiving and processing the interference effect of GNSS signals formed by direct and surface reflection [4][5][6]. At present, this technology is mainly employed to retrieve the soil moisture content (SMC), snow depth, and vegetation parameters [7][8][9][10][11][12].
In recent years, researchers have made remarkable achievements in GNSS-IR snow depth detection. Larson et al. first proposed to extend the application of signal-to-noise ratio (SNR) data observed by the traditional geodetic global positioning system (GPS) receiver in order to detect snow depth [3]. Larson et al. conducted snow depth retrieval using SNR data of a GPS L2C signal at multiple stations of the Plate Boundary Observatory (PBO) and further verified the feasibility of this technology [13]. Later, Larson et al. developed a snow depth retrieval algorithm based on GPS L1 SNR data, compared it with the snow depth results of L2C signal retrieval, and found that the accuracy was improved [14]. Tabibi et al. evaluated the value of GPS L5 SNR data in snow depth detection, compared the results with the L2C signal, and found that there was no detectable deviation in the L5 retrieval results [15]. Tabibi et al. proposed the global orbit navigation satellite system multipath reflectometry (GLONASS-MR) SNR retrieval model, which extended GPS-MR to multiple GNSS, and showed a strong correlation by comparing the retrieval results of GPS L2C and GLONASS R2-coarse acquisition (C/A) signals [16]. Later, Tabibi et al. used simulation and field measurements to evaluate the accuracy of GPS and GLONASS combined multi-GNSS-MR snow depth retrieval. At the same time, the variance factor was used to form the optimal multi-GNSS combined snow depth daily sequence retrieval. Compared with the single signal snow depth retrieval results, the accuracy was significantly improved [17]. Jin et al. used GPS L2P SNR data to retrieve the snow depth, and compared it with GPS L1 C/A code results, which showed a high correlation, indicating that using L2P SNR data can better estimate the snow depth [18]. Zhou et al. used GLONASS L1 SNR data for snow depth detection, and the accuracy reached centimeter-level and showed a strong correlation with the measured data [19]. Zhou et al. also studied the retrieval of different signal combinations and proposed using GPS L1, L2, and L5 signal multipath reflection and SNR combination for the retrieval of snow depth. The results showed that this method can be effectively used for snow depth detection [20].
The above scholars' research on snow depth retrieval is based on GPS and GLONASS observation data and has not been extended to other GNSS systems. To further analyze the potential of other GNSS systems in snow depth detection, Wang et al. used the SNR data of GPS L1 and BeiDou navigation satellite system (BDS) B1I signals to retrieve the snow depth, finally reaching an accuracy of 5 cm in a day [21]. Wang et al. used multi-GNSS system data to retrieve the snow depth and found that the trend of single signal retrieval results of multi-GNSS system constellations was in good agreement, except for the GPS precise code (P-code) signal. Then, the multi-GNSS system combination method based on robust regression was used to combine the signal retrieval between constellations. The results showed that the accuracy, availability, and time sampling of multi-GNSS system combination retrieval were improved [22].
From the current research status, snow depth retrieval is mainly concentrated in single or dual GNSS systems and single frequency SNR data, which is often found to be challenging when attempting to meet the accuracy and time resolution requirements of snow depth detection. Therefore, based on previous studies, this article conducts snow depth retrieval using multi-GNSS and multi-frequency SNR data. For the case that the antenna height of the GNSS receiver is unknown, the mean value of the multi-day Lomb-Scargle periodogram (LSP) spectrum analysis results in the snow-free surface is used as the initial reference reflector height of the multi-GNSS and multi-frequency GNSS-IR in the article. Then, the snow depth retrieval capability of the Galileo satellite navigation system (Galileo) and BDS multi-frequency signal, which are rarely used to retrieve the snow depth, except the GPS and GLONASS signal, are evaluated. A mean fusion of multi-frequency retrieval results of different GNSS systems is proposed to improve the accuracy compared with different frequencies of multi-GNSS system retrieval results. Finally, the multi-GNSS system retrieval results are fused to further evaluate the accuracy of the GNSS combination retrieval of snow depth. In this article, through the above process, the feasibility and accuracy of multi-GNSS and multi-frequency GNSS-IR snow depth retrieval are evaluated.

GNSS-IR Snow Depth Retrieval Principle
The snow depth retrieval method used in the article is based on SNR data processing of traditional geodetic GNSS receivers. The SNR is an indicator used to measure the signal strength of global navigation satellites, which is mainly affected by antenna gain, satellite transmission power, and multipath [23][24][25]. The multipath effect of the SNR decreases with increasing satellite elevation angle. The direct and surface reflected signals will have obvious interference effects at the receiver antenna when the satellite is at low elevation angle. At the same time, the frequency of the reflected signal will also change with the change in antenna height. The snow depth parameter can be obtained by comparing the vertical distance difference between the reflection surface and the receiver antenna phase center under snow-free and snow conditions. Figure 1 shows that direct signal and reflected signal generate corresponding interference effects at the receiver to form a composite interference signal, which can be expressed as [24]: where A d and A r are the amplitudes of the direct signal and the reflected signal, and ϕ is the difference between the phases of the direct signal and the phases of the reflected signal (the unit is rad), which can be expressed as [26]: where λ is the wavelength of the signal; θ is the elevation angle of the satellite; h is the vertical distance from the reflector to the antenna phase center; and the full text is uniformly called the reflector height. Because the change in snow depth parameters is only related to the reflected signal in the composite SNR data, it is necessary to eliminate the direct signal in the composite SNR to obtain the reflected signal part. In the article, the composite SNR data are fitted by a cubic low-order polynomial, and the composite SNR is linearized before fitting [27]: After the trend of the direct signal sequence is fitted, the SNR sequence of the reflected signal that removes the influence of the direct signal can be obtained, which can be expressed as SNR r [28]: f is the signal frequency of the multipath effect part in the SNR. After simplification, the relationship reflector height and the satellite signal wavelength can be obtained as follows: Figure 1. Schematic diagram of GNSS-IR snow depth retrieval. After the satellite sends the signal, the right-handed circular polarized (RHCP) antenna receives the direct signal and the surface reflected signal, and produces interference effect at the receiver. The snow surface reflector height (H s ) and snow-free surface reflector height (H sf ) are calculated, respectively, by analyzing the oscillation effect, and the snow depth (h sd ) is calculated by comparing the differences between them. θ is the elevation angle of the satellite.
In the article, the LSP method is used to analyze SNRr , obtain f, and extract h [3,29,30]. As shown in Figure 1, the snow depth parameter is calculated by using the above method by comparing reflector height in the case of snow-free and snow surfaces: where Hsf and Hs are the snow-free and snow surface reflector height, and hsd is the snow depth. Snow depth retrieval is mainly determined by snow and snow-free surface reflector height. Wang et al. mentioned that non-planar reflecting surface and atmospheric refraction in GNSS-IR techniques can lead to errors in reflector height. These two errors are mainly considered in the sea level height retrieval, but rarely in the snow depth, as the snow depth changes slowly, and the reflector height to snow depth is usually smaller than that to the sea level [22]. The purpose of the article is to evaluate the ability of multi-GNSS and multi-frequency GNSS-IR snow depth retrieval, without focusing on the actual error impact caused by these two factors.

Station Location and Surrounding Environment
This article selects GNSS observation data and snow depth data collected at the P387 station ( Figure 2) of the PBO, which is located in Sisters, Oregon, the U.S., with an altitude of 963.041 m. The specific location of the station is 44.29675°N and 121.57446°W. The data of the P387 station is collected by SEPTENTRIO (SEPT) POLARX5 receivers and TRM59800.80 antennas pointing toward the zenith with a sampling frequency of 15 s. The terrain around this station is flat without trees, and the signal acquisition conditions are good.  Figure 1. Schematic diagram of GNSS-IR snow depth retrieval. After the satellite sends the signal, the right-handed circular polarized (RHCP) antenna receives the direct signal and the surface reflected signal, and produces interference effect at the receiver. The snow surface reflector height (H s ) and snow-free surface reflector height (H s f ) are calculated, respectively, by analyzing the oscillation effect, and the snow depth (h sd ) is calculated by comparing the differences between them. θ is the elevation angle of the satellite.
In the article, the LSP method is used to analyze SNR r , obtain f , and extract h [3,29,30]. As shown in Figure 1, the snow depth parameter is calculated by using the above method by comparing reflector height in the case of snow-free and snow surfaces: where H s f and H s are the snow-free and snow surface reflector height, and h sd is the snow depth. Snow depth retrieval is mainly determined by snow and snow-free surface reflector height. Wang et al. mentioned that non-planar reflecting surface and atmospheric refraction in GNSS-IR techniques can lead to errors in reflector height. These two errors are mainly considered in the sea level height retrieval, but rarely in the snow depth, as the snow depth changes slowly, and the reflector height to snow depth is usually smaller than that to the sea level [22]. The purpose of the article is to evaluate the ability of multi-GNSS and multi-frequency GNSS-IR snow depth retrieval, without focusing on the actual error impact caused by these two factors.

Station Location and Surrounding Environment
This article selects GNSS observation data and snow depth data collected at the P387 station ( Figure 2) of the PBO, which is located in Sisters, Oregon, the U.S., with an altitude of 963.041 m. The specific location of the station is 44.29675 • N and 121.57446 • W. The data of the P387 station is collected by SEPTENTRIO (SEPT) POLARX5 receivers and TRM59800.80 antennas pointing toward the zenith with a sampling frequency of 15 s. The terrain around this station is flat without trees, and the signal acquisition conditions are good. Figure 2 shows that, regarding the P387 site location and surrounding environment, the terrain around P387 is flat. The ability of multi-GNSS and multi-frequency GNSS-IR snow depth retrieval is key to the article. In order to reduce the influence of terrain on the reflector height, the experimental region with small surface fluctuation is selected.  Figure 2 shows that, regarding the P387 site location and surrounding environment, the terrain around P387 is flat. The ability of multi-GNSS and multi-frequency GNSS-IR snow depth retrieval is key to the article. In order to reduce the influence of terrain on the reflector height, the experimental region with small surface fluctuation is selected.
At the same time, it can be seen that the vegetation around the P387 site is rare, and that the vegetation type is lawn. The surface can be defined as bare soil in winter snow stage. Therefore, with the melting of snow, the surface is gradually exposed, and the signal reflected by the surface is less affected by vegetation attenuation. The surface should be selected to be close to the bare soil in the subsequent analysis of snow-free reflection, so as to reduce the influence of vegetation error caused by the reflector height of snow-free and snow surface.
The roughness of snow surface will also lead to a retrieval error. In the article, the snow surface is regarded as a plane when extracting the reflector height, and it is not corrected temporarily.
For the above description, the calculation results of the reflector height in the snowfree and snow surfaces will not cause significant errors due to the change in the position of the mirror point. Therefore, the rise and fall stages of the satellite can be used as a signal source when selecting the experimental data. Nevertheless, it is necessary to determine whether it is a continuous observation period in advance in order to select the available arc segment. Figure 3 shows the PBO snow depth data between days of year (DOYs) 024 and 065 of 2017. The article also selects GNSS observation data at this time, and the observation period selected in the article is when the snow has reached the deepest state, followed by the process of ablation. The feasibility and accuracy of snow depth retrieval using multi-GNSS and multi-frequency SNR data are verified by the change in the snow depth. At the same time, it can be seen that the vegetation around the P387 site is rare, and that the vegetation type is lawn. The surface can be defined as bare soil in winter snow stage. Therefore, with the melting of snow, the surface is gradually exposed, and the signal reflected by the surface is less affected by vegetation attenuation. The surface should be selected to be close to the bare soil in the subsequent analysis of snow-free reflection, so as to reduce the influence of vegetation error caused by the reflector height of snow-free and snow surface.

Selection and Analysis of Experimental Data
The roughness of snow surface will also lead to a retrieval error. In the article, the snow surface is regarded as a plane when extracting the reflector height, and it is not corrected temporarily.
For the above description, the calculation results of the reflector height in the snowfree and snow surfaces will not cause significant errors due to the change in the position of the mirror point. Therefore, the rise and fall stages of the satellite can be used as a signal source when selecting the experimental data. Nevertheless, it is necessary to determine whether it is a continuous observation period in advance in order to select the available arc segment. Figure 3 shows the PBO snow depth data between days of year (DOYs) 024 and 065 of 2017. The article also selects GNSS observation data at this time, and the observation period selected in the article is when the snow has reached the deepest state, followed by the process of ablation. The feasibility and accuracy of snow depth retrieval using multi-GNSS and multi-frequency SNR data are verified by the change in the snow depth.

Selection and Analysis of Experimental Data
During the experiment period, when the GNSS signal is transmitted to the surface receiver, the signal will pass through the atmosphere, and a signal refraction effect will occur when passing through the troposphere. Williams et al. considered that tropospheric delay will cause height error in the obtained vertical reflection distance [31]. Aiming at this problem, this article gives the tropospheric delay information during the experiment, as shown in Figure 4. During the experiment period, when the GNSS signal is transmitted to the surface receiver, the signal will pass through the atmosphere, and a signal refraction effect will occur when passing through the troposphere. Williams et al. considered that tropospheric delay will cause height error in the obtained vertical reflection distance [31]. Aiming at this problem, this article gives the tropospheric delay information during the experiment, as shown in Figure 4.   During the experiment period, when the GNSS signal is transmitted to the surface receiver, the signal will pass through the atmosphere, and a signal refraction effect will occur when passing through the troposphere. Williams et al. considered that tropospheric delay will cause height error in the obtained vertical reflection distance [31]. Aiming at this problem, this article gives the tropospheric delay information during the experiment, as shown in Figure 4.  As can be seen from Figure 4, the troposphere delay slowly changes during the 42 days of experiment, and is basically the same. He et al. corrected the tropospheric delay error in the process of retrieving coastal typhoon storm surge by using GNSS-IR signal, and the final accuracy was only improved by approximately 0.5 cm, which can ignore its influence [32]. Therefore, the error caused by tropospheric delay is not especially corrected in the process of snow depth retrieval.

Reflection Region Analysis
The effective reflection region of GNSS signal to the surface can be described by the first Fresnel reflection region, which is a group of ellipses related to the receiver antenna height, satellite azimuth, and satellite elevation angle. Figure 5a shows the Fresnel reflection region of GPS G10 satellite with DOY of 024 in 2017 at P387 station. Assuming that receiver antenna height is 2 m, different color lines represent the reflection regions with different satellite elevation angles. With the increase in the elevation angle, the Fresnel reflection region will be smaller. The figure shows the Fresnel reflection region map of the satellite elevation angle of 5-25 degrees.
It can be seen that the effective reflection region is related to the satellite elevation angle. When the satellite elevation angle gradually increases, the effective reflection region will be decrease in size, and will gradually approach the receiver antenna. Large interference will occur when the satellite is at a low elevation angle. The larger the satellite elevation angle, the less affected the satellite is by the multipath of surrounding signals. Therefore, the signal data with low satellite elevation angle should be selected in the process of snow depth retrieval.  Figure 5b shows that the trajectory of reflection points changes with the change in the relative position between GNSS satellite and receiver, showing both different directions and distributions at different arcs and the ground reflection point trajectory formed by some satellites of the four GNSS systems. The combined signal sensing range of the four GNSS systems is significantly expanded, which can provide more data sources and wider sensing region, which is conducive to improving the time resolution of snow depth retrieval.

SNR Types
PBO provides multi-GNSS and multi-frequency SNR data in the observation data of receiver independent exchange format (RINEX) version 3.03. The SNR types of P387 station mainly include S1C, S1W, S2L, S2W, and S5Q for GPS; S1C and S2C for GLONASS; S1C, S5Q, S6C, S7Q, and S8Q for Galileo; S2I, S6I, and S7I for BDS; S1C, S2L, and S5Q for quasi-zenith satellite system (QZSS); and S1C and S5I for satellite-based augmentation It can be seen that the effective reflection region is related to the satellite elevation angle. When the satellite elevation angle gradually increases, the effective reflection region will be decrease in size, and will gradually approach the receiver antenna. Large interference will occur when the satellite is at a low elevation angle. The larger the satellite elevation angle, the less affected the satellite is by the multipath of surrounding signals. Therefore, the signal data with low satellite elevation angle should be selected in the process of snow depth retrieval. Figure 5b shows that the trajectory of reflection points changes with the change in the relative position between GNSS satellite and receiver, showing both different directions and distributions at different arcs and the ground reflection point trajectory formed by some satellites of the four GNSS systems. The combined signal sensing range of the four GNSS systems is significantly expanded, which can provide more data sources and wider sensing region, which is conducive to improving the time resolution of snow depth retrieval.

SNR Types
PBO provides multi-GNSS and multi-frequency SNR data in the observation data of receiver independent exchange format (RINEX) version 3.03. The SNR types of P387 station mainly include S1C, S1W, S2L, S2W, and S5Q for GPS; S1C and S2C for GLONASS; S1C, S5Q, S6C, S7Q, and S8Q for Galileo; S2I, S6I, and S7I for BDS; S1C, S2L, and S5Q for quasi-zenith satellite system (QZSS); and S1C and S5I for satellite-based augmentation system (SBAS). The data sampling interval is 15 s, and the specific information is shown in Table 1. Table 1 shows that the P387 station provides SNR data of different signal types of six satellite systems. After reading the data, it is found that the QZSS and SBAS have fewer satellites and no available arc segments. Therefore, this article does not consider using QZSS and SBAS data for snow depth retrieval; instead, the SNR data of GPS, GLONASS, Galileo, and BDS are used. In the article, multi-GNSS and multi-frequency SNR data are used to retrieve the snow depth. In addition to GPS, the observation satellites will be extended to other systems, which is of great help in improving the retrieval accuracy of snow depth and expanding the observation range and time resolution. Note: When reading RINEX 3.03, BDS 1I/Q/X and 2I/Q/X can be regarded as the same as 2I/Q/X in the current RINEX standard, and the AS in Table 1 is anti-spoofing. Figure 6 shows the experimental technical scheme of multi-GNSS and multi-frequency GNSS-IR snow depth retrieval. It can be seen that the technical route of the article can be divided into three parts: (1) GNSS-IR data preprocessing is carried out, where the SNR, pseudo-random noise (PRN), satellite elevation angle, azimuth angle, and other data parameters are extracted from the observation (OBS) file and navigation (NAV) file collected by GNSS receivers; (2) the LSP method is used to analyze both the reflector height of snow-free and snow surfaces and the difference between them in order to retrieve the snow depth; (3) the multi-GNSS and multi-frequency GNSS-IR snow depth retrieval results and PBO snow depth data are compared, and then the mean fusion analysis of the multi-GNSS and multi-frequency GNSS-IR snow depth retrieval results is carried out. Figure 6 shows that RTKLIB software is used for data processing in the article. By reading the OBS file and NAV file in RINEX 3.03 format, the corresponding elevation angle, azimuth angle, SNR, and other data can be extracted. In order to select the satellites in the four GNSS systems that have available observation arc data in 42 days of the experimental stage, this article finally sets the G10 satellite of GPS, R17 satellite of GLONASS, E12 satellite of Galileo, and C14 satellite of BDS as the experimental data source satellites. In the article, the SNR data in the rising stage are mainly selected for processing. When the retrieval results show abnormal values, the SNR data in the falling stage can be selected as a supplement for retrieval. At the same time, the minimum elevation angle threshold of 5 degrees can also be increased accordingly in order to achieve a more practical retrieval value.

Experimental Technical Scheme
Remote Sens. 2021, 13, 4311 9 of 21 snow depth; (3) the multi-GNSS and multi-frequency GNSS-IR snow depth retrieval results and PBO snow depth data are compared, and then the mean fusion analysis of the multi-GNSS and multi-frequency GNSS-IR snow depth retrieval results is carried out. Figure 6. The technical process of multi-GNSS and multi-frequency GNSS-IR snow depth retrieval. Among the ele is the satellite elevation angle. Figure 6 shows that RTKLIB software is used for data processing in the article. By reading the OBS file and NAV file in RINEX 3.03 format, the corresponding elevation angle, azimuth angle, SNR, and other data can be extracted. In order to select the satellites in the four GNSS systems that have available observation arc data in 42 days of the experimental stage, this article finally sets the G10 satellite of GPS, R17 satellite of GLONASS, E12 satellite of Galileo, and C14 satellite of BDS as the experimental data source satellites. In the article, the SNR data in the rising stage are mainly selected for processing. When the retrieval results show abnormal values, the SNR data in the falling stage can be selected as a supplement for retrieval. At the same time, the minimum elevation angle Retrieval accuracy analysis Figure 6. The technical process of multi-GNSS and multi-frequency GNSS-IR snow depth retrieval. Among the ele is the satellite elevation angle.

Extraction h 3.2.1. Multi-GNSS and Multi-Frequency SNR Sequence Extraction
In the article, multi-GNSS and multi-frequency SNR sequences are extracted by RTK-LIB software. Figure 7 shows the rising SNR sequences of G10 (S1C, S2L, and S5Q), R17 (S1C and S2C), E12 (S1C, S5Q, S6C, S7Q, and S8Q), and C14 (S2I, S6I, and S7I) satellites in the four GNSS systems on DOY 024 in 2017. Through data preprocessing, it is found that the SNR values of G10 (S1W and S2W) are the same, and that the SNR difference with other signal frequencies is significant. Therefore, these two types of SNR data are not used for the experimental in the article. Figure 7 shows the SNR sequence change in the elevation angle from low to high. R17 (S1C and S2C), E12 (S1C, S5Q, S6C, S7Q, and S8Q), and C14 (S2I, S6I, and S7I) satellites in the four GNSS systems on DOY 024 in 2017. Through data preprocessing, it is found that the SNR values of G10 (S1W and S2W) are the same, and that the SNR difference with other signal frequencies is significant. Therefore, these two types of SNR data are not used for the experimental in the article. Figure 7 shows the SNR sequence change in the elevation angle from low to high.  Figure 7 shows that the SNR sequence has strong oscillation at a low elevation angle. It is greatly affected by multipath at a low elevation angle, and the interference caused by the direct signal and the reflected signal is obvious. When increasing the elevation angle, the interference degree gradually decreases. Therefore, after removing the direct signal from the SNR sequence, this article mainly selects the SNR sequence of the reflected signal in the range of elevation angles of 5-30 degrees for snow depth retrieval. At the same time, it can be seen that the SNR sequences of the four GNSS systems at different frequencies have specific differences.  Figure 7 shows that the SNR sequence has strong oscillation at a low elevation angle. It is greatly affected by multipath at a low elevation angle, and the interference caused by the direct signal and the reflected signal is obvious. When increasing the elevation angle, the interference degree gradually decreases. Therefore, after removing the direct signal from the SNR sequence, this article mainly selects the SNR sequence of the reflected signal in the range of elevation angles of 5-30 degrees for snow depth retrieval. At the same time, it can be seen that the SNR sequences of the four GNSS systems at different frequencies have specific differences.

SNR Sequence Data Processing
In the article, the composite SNR sequence is linearized first, and then the linearized SNR sequence is fitted by a cubic low-order polynomial to obtain the direct signal part. The SNR sequence of the reflected signal is obtained by subtracting the composite SNR sequence from the direct signal part. Figure 8 shows the processing of GPS S1C SNR data.

LSP Analysis Results of the Snow Surface
Based on the LSP method analysis of the SNR sequence of the multi-GNSS and multi-frequency reflection signal of DOY 024 in 2017, the results are shown in Figure 9. Figure 9 shows that there are some differences in the results of the multi-GNSS and multi-frequency LSP analysis. Specifically, the S2L and S5Q results in GPS are 1.335 m and 1.330 m, which are consistent and close to 0.1 m compared with the S1C results. The results of S1C and S2C in GLONASS are 1. deviation of S5Q, S6C, S7Q, and S8Q is slight. The results of S6I and S7I in BDS are 1.275 m and 1.310 m, respectively, and the difference is slight. The difference between the results of S6I and S2I is 0.095 m. From the above data, it can be seen that the LSP analysis results of different GNSS systems are different, and that the results of different frequencies in each GNSS system are also different, but the overall consistency is good. In the article, through the LSP method analysis, the reflector height in the multi-GNSS and multi-frequency can be obtained, and the LSP of the snow surface can be obtained.

SNR Sequence Data Processing
In the article, the composite SNR sequence is linearized first, and then the linearized SNR sequence is fitted by a cubic low-order polynomial to obtain the direct signal part. The SNR sequence of the reflected signal is obtained by subtracting the composite SNR sequence from the direct signal part. Figure 8 shows the processing of GPS S1C SNR data.

LSP Analysis Results of the Snow Surface
Based on the LSP method analysis of the SNR sequence of the multi-GNSS and multifrequency reflection signal of DOY 024 in 2017, the results are shown in Figure 9.

SNR Sequence Data Processing
In the article, the composite SNR sequence is linearized first, and then the linearized SNR sequence is fitted by a cubic low-order polynomial to obtain the direct signal part. The SNR sequence of the reflected signal is obtained by subtracting the composite SNR sequence from the direct signal part. Figure 8 shows the processing of GPS S1C SNR data.

LSP Analysis Results of the Snow Surface
Based on the LSP method analysis of the SNR sequence of the multi-GNSS and multifrequency reflection signal of DOY 024 in 2017, the results are shown in Figure 9.

LSP Analysis Results of the Snow-Free Surface
In order to weaken the difference between the multi-GNSS and multi-frequency LSP results, this article proposes that the snow-free surface reflector height is analyzed by the multi-GNSS and multi-frequency LSP results. Using the above method, the four-day accumulated data (satellites G10, R17, E12, and C14) of DOY 225-228 in 2017 were selected for processing; this period is about October and belongs to the defined bare soil period. This article selects the data of this period to process, which can reduce the influence of surface vegetation on signal propagation. The multi-GNSS and multi-frequency reference values were obtained by averaging the four-day LSP results. Figure 10 shows that there are some differences in the results of multi-GNSS and multi-frequency LSP. There are some differences between the results of GPS S1C, S2L, and S5Q, but the results of a single frequency at 4 days are basically the same, with slight deviation. The Galileo S5Q, S6C, S7Q, and S8Q results are basically the same, and S1C has a specific difference, but the results of a single frequency at 4 days are basically the same. The deviation of the BDS S6I and S7I results is slight, and the deviation of the S2I results is significant at 225 and 228 on the annual accumulation day. The results of the GLONASS S1C analysis showed a significant variation at 227 accumulated days, and the rest showed good consistency. In the article, the snow depth is obtained through the comparison of the reflector height of different GNSS systems at different frequencies under snow and snow-free surfaces, which can better adapt to the use of multi-GNSS and multi-frequency GNSS-IR technology. The reference value of the reflector height in different GNSS systems at different frequencies is calculated, as shown in Table 2.
each GNSS system are also different, but the overall consistency is good. In the article, through the LSP method analysis, the reflector height in the multi-GNSS and multi-frequency can be obtained, and the LSP of the snow surface can be obtained.

LSP Analysis Results of the Snow-Free Surface
In order to weaken the difference between the multi-GNSS and multi-frequency LSP results, this article proposes that the snow-free surface reflector height is analyzed by the multi-GNSS and multi-frequency LSP results. Using the above method, the four-day accumulated data (satellites G10, R17, E12, and C14) of DOY 225-228 in 2017 were selected for processing; this period is about October and belongs to the defined bare soil period. This article selects the data of this period to process, which can reduce the influence of surface vegetation on signal propagation. The multi-GNSS and multi-frequency reference values were obtained by averaging the four-day LSP results. Figure 10 shows that there are some differences in the results of multi-GNSS and multi-frequency LSP. There are some differences between the results of GPS S1C, S2L, and S5Q, but the results of a single frequency at 4 days are basically the same, with slight deviation. The Galileo S5Q, S6C, S7Q, and S8Q results are basically the same, and S1C has a specific difference, but the results of a single frequency at 4 days are basically the same. The deviation of the BDS S6I and S7I results is slight, and the deviation of the S2I results is significant at 225 and 228 on the annual accumulation day. The results of the GLONASS S1C analysis showed a significant variation at 227 accumulated days, and the rest showed good consistency. In the article, the snow depth is obtained through the comparison of the reflector height of different GNSS systems at different frequencies under snow and snowfree surfaces, which can better adapt to the use of multi-GNSS and multi-frequency GNSS-IR technology. The reference value of the reflector height in different GNSS systems at different frequencies is calculated, as shown in Table 2. Reflector height/(m) GPS S1C GPS S2L GPS S5Q GLONASS S1C GLONASS S2C Galileo S1C Galileo S5Q Galileo S6C Galileo S7Q Galileo S8Q BDS S2I BDS S6I BDS S7I Figure 10. The snow-free surface reflector height reference value of multi-GNSS and multi-frequency LSP analysis results.  Table 2 shows that the results of LSP at different frequencies of the four GNSS systems are more than 1800 m, which is more in line with the actual situation and can be used as the initial reference value of the reflector height of the snow-free surface.

Multi-GNSS and Multi-Frequency GNSS-IR Snow Depth Retrieval Results
The snow depth is obtained by comparing and analyzing the difference in reflector height under snow-free and snow surfaces. The results of the different frequency retrievals of GPS, GLONASS, Galileo, and BDS are compared with the PBO snow depth data, and the results are shown in Figure 11.   Table 2 shows that the results of LSP at different frequencies of the four GNSS systems are more than 1800 m, which is more in line with the actual situation and can be used as the initial reference value of the reflector height of the snow-free surface.

Multi-GNSS and Multi-Frequency GNSS-IR Snow Depth Retrieval Results
The snow depth is obtained by comparing and analyzing the difference in reflector height under snow-free and snow surfaces. The results of the different frequency retrievals of GPS, GLONASS, Galileo, and BDS are compared with the PBO snow depth data, and the results are shown in Figure 11.  Figure 11 shows that the snow depth retrieval from the multi-GNSS and multifrequency SNR data has a trend that is similar to that of the PBO snow depth data. Among them, the results of S1C in GPS are more biased than those of S2L and S5Q. The trend of the S1C results in GLONASS is worse than that of S2C. The results of S1C in Galileo are worse than those of S5Q, S6C, S7Q, and S8Q. The trend of the BDS S2I results is worse than that of the S6I and S7I results.

Mean Fusion of Multi-Frequency Retrieval Results in the Four GNSS Systems
The results of the mean fusion of different frequencies in the four GNSS systems are shown in Figure 12. Figure 12 shows the consistency between the retrieval results of the four GNSS systems and the PBO snow depth, where the trend is more consistent than the previous single frequency retrieval results.

Mean Fusion Retrieval Results of Multi-GNSS System
After the mean fusion of multi-frequency in the four GNSS systems, the retrieval results of highly similar trends are obtained. To further improve the accuracy, the mean fusion of the snow depth retrieval results of the multi-GNSS system is carried out, and the final GNSS system retrieval results are obtained, as shown in Figure 13. Figure 11 shows that the snow depth retrieval from the multi-GNSS and multi-frequency SNR data has a trend that is similar to that of the PBO snow depth data. Among them, the results of S1C in GPS are more biased than those of S2L and S5Q. The trend of the S1C results in GLONASS is worse than that of S2C. The results of S1C in Galileo are worse than those of S5Q, S6C, S7Q, and S8Q. The trend of the BDS S2I results is worse than that of the S6I and S7I results.

Mean Fusion of Multi-Frequency Retrieval Results in the Four GNSS Systems
The results of the mean fusion of different frequencies in the four GNSS systems are shown in Figure 12.  Figure 12 shows the consistency between the retrieval results of the four GNSS systems and the PBO snow depth, where the trend is more consistent than the previous single frequency retrieval results.

Mean Fusion Retrieval Results of Multi-GNSS System
After the mean fusion of multi-frequency in the four GNSS systems, the retrieval results of highly similar trends are obtained. To further improve the accuracy, the mean fusion of the snow depth retrieval results of the multi-GNSS system is carried out, and the final GNSS system retrieval results are obtained, as shown in Figure 13. shown in Figure 12.  Figure 12 shows the consistency between the retrieval results of the four GNSS systems and the PBO snow depth, where the trend is more consistent than the previous single frequency retrieval results.

Mean Fusion Retrieval Results of Multi-GNSS System
After the mean fusion of multi-frequency in the four GNSS systems, the retrieval results of highly similar trends are obtained. To further improve the accuracy, the mean fusion of the snow depth retrieval results of the multi-GNSS system is carried out, and the final GNSS system retrieval results are obtained, as shown in Figure 13.  Figure 13 shows that the mean fusion results between the multi-GNSS system are in good agreement with the PBO snow depth results, and that the trend is basically the same.

Accuracy Analysis between Multi-GNSS and Multi-Frequency GNSS-IR Snow Depth Retrieval Results and PBO Snow Depth
The above retrieval results are further analyzed, and the retrieval results of different frequency signals in multi-GNSS system are compared; the correlation coefficient (R) and root mean square error (RMSE) are shown in Figure 14.

Accuracy Analysis between Multi-GNSS and Multi-Frequency GNSS-IR Snow Depth Retrieval Results and PBO Snow Depth
The above retrieval results are further analyzed, and the retrieval results of differen frequency signals in multi-GNSS system are compared; the correlation coefficient (R) and root mean square error (RMSE) are shown in Figure 14.   Figure 14 shows the correlation of multi-GNSS and multi-frequency GNSS-IR snow depth retrieval results compared with the PBO snow depth. Table 3 shows the specific R and RMSE.  Figure 14 and Table 3 show that the snow depth results of the GNSS-IR retrieval at different frequencies of the four GNSS systems have a strong correlation with the PBO snow depth. The R of the GPS S2L and S5Q results were 0.99 and 0.98, respectively, which are highly correlated, whereas the R of the S1C was 0.90, which is relatively low. The R of GLONASS S1C and S2C were 0.83 and 0.97, respectively. The correlation of S2C was strong and that of S1C was weak. The R of Galileo S5Q, S6C, S7Q, and S8Q results was approximately 0.95, and that of S1C was 0.88, which is rather weak. The R of BDS S6I and S7I results was 0.95 and 0.93, showing a strong correlation, and that of S2I was 0.86, which was rather weak. At the same time, the RMSE of the results of different frequencies of the four GNSS systems and PBO data were basically in the range of 5 cm to 10 cm, and the error was small. The above data show that the feasibility of retrieving the snow depth using multi-GNSS and multi-frequency GNSS-IR is high. At the same time, for GPS S1C, GLONASS S1C, Galileo S1C, and BDS S2I, the results in the four GNSS systems are rather weak. From the above results, in addition to GPS and GLONASS signals, the Galileo and BDS signals also have a good ability in snow depth retrieval.

Accuracy Analysis of Multi-Frequency Mean Fusion Results in the GNSS Systems
Mean fusion retrieval is carried out for multi-frequency retrieval results under the four GNSS systems; the accuracy between the retrieval results and PBO snow depth are further analyzed, as shown in Figure 15. Figure 15 shows that the retrieval results of the four GNSS systems are obtained by averaging the snow depth results of the GNSS-IR retrieval at multi-frequency. The R between the GPS, GLONASS, Galileo, BDS, and PBO results was 0.99, 0.94, 0.97, and 0.95, respectively, showing strong correlations. At the same time, the RMSE of the four GNSS systems have been reduced to a certain extent, to basically in the range of 4 cm to 7 cm, indicating that the mean fusion of different frequency retrieval results in the four GNSS systems has a good effect on the improvement of retrieval accuracy, which can eliminate the weak correlation of GPS S1C, GLONASS S1C, Galileo S1C, and BDS S2I.
Compared with the snow depth of the different frequency signal in the multi-GNSS system retrieval results, the GNSS multi-frequency mean fusion method can effectively improve the retrieval accuracy. The specific improvement accuracy is shown in Table 4.
As can be seen from Table 4, the R of the GNSS multi-frequency mean fusion increases by 11.7%, which is higher than a single frequency in multi-GNSS system retrieval results, and the GLONASS S2C results decrease by 3.2%, but most of the results' correlations are improved. At the same time, the RMSE decreases by 55.6%, and though the retrieval accuracy of Galileo S7Q and BDS S6I is not improved, the other accuracy has been greatly improved. This shows that the mean fusion of the GNSS multi-frequency retrieval results can effectively improve the accuracy of the snow depth retrieval.    As can be seen from Table 4, the R of the GNSS multi-frequency mean fusion increases by 11.7%, which is higher than a single frequency in multi-GNSS system retrieval results, and the GLONASS S2C results decrease by 3.2%, but most of the results' correlations are improved. At the same time, the RMSE decreases by 55.6%, and though the retrieval accuracy of Galileo S7Q and BDS S6I is not improved, the other accuracy has been greatly improved. This shows that the mean fusion of the GNSS multi-frequency retrieval results can effectively improve the accuracy of the snow depth retrieval.

Accuracy Analysis of Mean Fusion Retrieval between Multi-GNSS System
The accuracy of multi-frequency fusion retrieval results has been improved to a certain extent. Furthermore, the fusion accuracy of the retrieval results of the four GNSS systems is analyzed, as shown in Figure 16.

Accuracy Analysis of Mean Fusion Retrieval between Multi-GNSS System
The accuracy of multi-frequency fusion retrieval results has been improved to a certain extent. Furthermore, the fusion accuracy of the retrieval results of the four GNSS systems is analyzed, as shown in Figure 16.  Figure 16 shows that the R between the mean fusion results of the multi-GNSS system and PBO snow depth was 0.99, showing a strong correlation. At the same time, the RMSE was 3 cm, and the error was also significantly reduced, indicating the effectiveness of this method.
In order to further improve the combined retrieval accuracy of the GNSS system, the GNSS multi-frequency fusion results are further mean fused. The fusion results are compared with the four GNSS systems, and the results are shown in Table 5. It can be seen from Table 5 that the accuracy of the results after multi-GNSS system fusion is further improved. Compared with the results of GNSS multi-frequency fusion, the R of the four GNSS systems are improved, except for GPS. GLONASS has the highest increase of 5.1%. In terms of RMSE, GLONASS increased by 57.1%. The GPS is not improved on R, but it is improved by 25% on RMSE. It can be seen from the above results that the further fusion of the multi-frequency retrieval results of the four GNSS systems can effectively improve the retrieval accuracy. Therefore, the multi-GNSS system combined snow depth retrieval method has strong reliability.  Figure 16 shows that the R between the mean fusion results of the multi-GNSS system and PBO snow depth was 0.99, showing a strong correlation. At the same time, the RMSE was 3 cm, and the error was also significantly reduced, indicating the effectiveness of this method.
In order to further improve the combined retrieval accuracy of the GNSS system, the GNSS multi-frequency fusion results are further mean fused. The fusion results are compared with the four GNSS systems, and the results are shown in Table 5. Table 5. Comparison of multi-frequency mean fusion; multi-GNSS system fusion accuracy increases in R and decreases in the RMSE. It can be seen from Table 5 that the accuracy of the results after multi-GNSS system fusion is further improved. Compared with the results of GNSS multi-frequency fusion, the R of the four GNSS systems are improved, except for GPS. GLONASS has the highest increase of 5.1%. In terms of RMSE, GLONASS increased by 57.1%. The GPS is not improved on R, but it is improved by 25% on RMSE. It can be seen from the above results that the further fusion of the multi-frequency retrieval results of the four GNSS systems can effectively improve the retrieval accuracy. Therefore, the multi-GNSS system combined snow depth retrieval method has strong reliability.

Conclusions
Snow is an important indicator to measure the global hydrological cycle and climate. The accurate long-term monitoring of snow depth is helpful for water resource management and climate disaster warning, and has important application prospects. Based on the current GNSS-IR snow depth retrieval method, the temporal resolution is affected by the number of sky satellite arcs, so multi-GNSS and multi-frequency GNSS-IR are introduced as a supplement. The snow depth retrieval is carried out by using the multi-GNSS and multi-frequency SNR data. The snow depth parameters are obtained by comparing the LSP results of snow-free and snow surfaces. At the same time, the results are compared and analyzed in terms of the PBO snow depth, and the correlation and error analysis of the multi-GNSS and multi-frequency GNSS-IR mean fusion results are carried out. The following conclusions are drawn through experimental analysis: (1) QZSS and SBAS systems in multi-GNSS and multi-frequency SNR data provided by PBO are not suitable for use due to the lack of observation arcs. The GPS S1W and S2W data values are the same, and the other frequency SNR data difference is too large and should not be used; (2) The LSP results of the snow-free surface can be effectively used as the initial reflector height reference value. The snow depth results of multi-GNSS and multi-frequency GNSS-IR retrieval have a strong correlation with PBO snow depth data, and the RMSE of different frequency retrieval results in the multi-GNSS system is between 5 cm and 10 cm. The correlation between the retrieval results of the GPS L1, GLONASS G1, Galileo E1, and BDS B1 bands in the snow depth retrieval results is rather weak; (3) The mean fusion of multi-frequency retrieval results in GPS, GLONASS, Galileo, and BDS can effectively improve the accuracy and solve the relatively weak results in some bands. The four GNSS systems retrieval results show a strong correlation, and the RMSE is between 4 cm and 7 cm. Comparing the different frequency signals in the multi-GNSS system retrieval results, the multi-frequency mean fusion increase by 11.7% in R and the RMSE decreases by 55.6%, which is the highest; (4) The mean fusion accuracy of the retrieval results of the GPS, GLONASS, Galileo, and BDS is significantly improved. The R between the retrieval results and the PBO results is 0.99, and the retrieval accuracy is better than 3 cm, which significantly enhances the accuracy. In the comparison of the multi-frequency mean fusion, the multi-GNSS system fusion increases by 5.1% in R and the RMSE decreases by 57.1%, which is the highest.
Compared with a single GNSS-IR signal, multi-GNSS and multi-frequency GNSS-IR improves the accuracy, continuity, and time resolution of snow depth retrieval. A mean fusion of multi-GNSS and multi-frequency GNSS-IR retrieval results can further enhance the accuracy. With the development of global navigation systems, more types of signals and perfect constellation structures will be provided. Multi-GNSS and multi-frequency GNSS-IR will play a more critical role in the field of snow depth detection.

Data Availability Statement:
The GNSS site data were provided under the PBO Observation Program of the United States, and the measured snow depth data were obtained from https: //data.unavco.org/archive/gnss/products (accessed on 21 October 2021).