Innovation: How Deep Is That White Stuff?

September 1, 2014  - By

Using GPS Multipath for Snow-Depth Estimation

By Felipe G. Nievinski and Kristine M. Larson



FRINGES. No, I’m not talking about the latest celebrity hairstyles nor the canopy of an American doorless, four-wheeled carriage from yesteryear (think Oklahoma!). I’m talking about interference fringes. But there is a connection to these other uses of the word fringe as we’ll see. You’ve all seen interference fringes at your local gas station, typically after it has just rained. They are the alternating bands of color we perceive when looking at a gasoline or oil slick in a puddle of water. They are caused by the white light from the Sun or artificial lighting reflected from the top surface of the slick and that from the bottom surface at the slick-water interface combining or interfering with each other at our eyeballs. The two sets of light waves arrive slightly out of phase with each other, and depending on the wavelengths of the reflected light and our angle of view, produce the colorful fringes. If the incident light was monochromatic, consisting of a single frequency or wavelength, then we would perceive just alternating bright and dark bands. The bright bands result from constructive interference when the phase difference is a near a multiple of 2π whereas the dark bands result from destructive interference when the difference is near an odd multiple of π.

Interference fringes had been seen long before the invention of the automobile. They are clearly seen on soap bubbles and the iridescent colors of peacock feathers, Morpho butterflies, and jewel beetles are also due to the interference phenomenon rather than pigmentation. Sir Isaac Newton did experiments on interference fringes (amongst other things) and tried to explain their existence — wrongly, it turned out. But he did coin the term fringes since they resembled the decorative fringe sometimes used on clothing, drapery, and, yes, surrey canopies.

It was the English polymath, Thomas Young, who, in 1801, first demonstrated interference as a consequence of the wave-nature of light with his famous double-slit experiment. You may have replicated his experiment in a high-school physics class. I did and I think I did it again as an undergraduate student taking a course in optics. Already by that point I was aiming for a career in physics or space science but I didn’t know that as a graduate student I would do research involving interference fringes. But not using light waves.

My research involved the application of very long baseline interferometry or VLBI to geodesy. VLBI had been developed by radio astronomers to better understand the structure of quasars and other esoteric celestial objects. At either ends of a baseline connecting large radio telescopes, perhaps stretching between continents, the quasar signals were recorded on magnetic tape and precisely registered using atomic clocks. When the tapes were played back and the signals aligned, one obtained interference fringes as peaks and troughs in an analog or digital waveform. Computer analysis of these fringes not only provided information on the structure of the observed radio source but also on the distance between the radio telescopes — eventually accurate enough to measure continental drift. 

But what has all of this got to do with GPS? In this month’s column, we look at a technique that uses fringes generated by signals arriving at an antenna directly from GPS satellites and those reflected by snow surrounding the antenna to measure its depth and how it varies over time. GPS for measuring snow depth; who would have thought?

“Innovation” is a regular feature that discusses advances in GPS technology and its applications as well as the fundamentals of GPS positioning. The column is coordinated by Richard Langley of the Department of Geodesy and Geomatics Engineering, University of New Brunswick. He welcomes comments and topic ideas.

Snowpacks are a vital resource for human existence on our planet. They provide reservoirs of fresh water, storing solid precipitation and delaying runoff. One sixth of the world population depends on this resource. Both scientists and water-supply managers need to know how much fresh water is stored in snowpack and how fast it is being released as a result of melting.

Snow monitoring from space is currently under investigation by both NASA and ESA. Greatly complementary to such spaceborne sensors are automated ground-based methods; the latter not only serve as essential independent validation and calibration for the former, but are also valuable for climate studies and flood/drought monitoring on their own. It is desirable for such estimates to be provided at an intermediary scale, between point-like in situ samples and wider area pixels.

In the last decade, GPS multipath reflectometry (GPS-MR), also known as GPS interferometric reflectometry and GPS interference-pattern technique, has been proposed for monitoring snow. This method tracks direct GPS signals, those that travel directly to an antenna, that have interfered with a coherently reflected signal, turning the GPS unit into an interferometer (see FIGURE 1). Its main variant is based on signal-to-noise ratio (SNR) measurements, although GPS-MR is also possible with carrier-phase and pseudorange observables. Data are collected at existing GPS base stations that employ commercial-off-the-shelf receivers and antennas in a conventional, antenna-upright setup. Other researchers have used a custom antenna and/or a dedicated setup, with the antenna tipped for enhanced multipath reception.

FIGURE 1. Standard geodetic receiver installation. The antenna is protected by a hemispherical radome. The monument (tripod structure) is ~ 2 meters above the ground. GPS satellites rise and set in ascending and descending sky tracks, multiple times per day. The specular reflection point migrates radially away from the receiver for decreasing satellite elevation angle. The total reflector height is made up of an a priori value and an unknown bias driven by the thickness of the snow layer.

FIGURE 1. Standard geodetic receiver installation. The antenna is protected by a hemispherical radome. The monument (tripod structure) is ~ 2 meters above the ground. GPS satellites rise and set in ascending and descending sky tracks, multiple times per day. The specular reflection point migrates radially away from the receiver for decreasing satellite elevation angle. The total reflector height is made up of an a priori value and an unknown bias driven by the thickness of the snow layer.

In this article, we summarize the SNR-based GPS-MR technique as applied to snow sensing using geodetic instruments. This forward/inverse approach for GPS-MR is new in that it capitalizes on known information about the antenna response and the physics of surface scattering to aid in retrieving the unknown snow conditions in the site surroundings. It is a statistically rigorous retrieval algorithm, agreeing to first order with the simpler original methodology, which is retained here for the inversion bootstrapping. The first part of the article describes the retrieval algorithm, while the second part provides validation at a representative site over an extended period of time. 

Physical Forward Model

SNR observations are formulated as SNR = Ps/Pn. In the denominator, we have the noise power, Pn, here taken as a constant, based on nominal values for the noise power spectral density and the noise bandwidth. The numerator is composite signal power:

Eq-1.   (1)

Its incoherent component is the sum of the respective direct and reflected powers (although direct incoherent power is negligible). In contrast, the coherent composite signal power follows from the complex sum of direct and reflection average voltages (not to be confused with the electromagnetic propagating fields, which neglect the receiving antenna response and also the receiver tracking process):


It is expressed in terms of the coherent direct and reflected powers, as well as the interferometric phase,

Eq-3 , (3)

which amounts to the reflection excess phase with respect to the direct signal.

We decompose observations, SNR = tSNR + dSNR, into a trend

Eq-4  (4)

over which interference fringes are superimposed:

Eq-5. (5) 

From now on, we neglect the incoherent power, which only impacts tSNR, not dSNR, and drop the coherent power superscript, for brevity.

The direct or line-of-sight power is formulated as

Eq-6  (6)

where  Eq-6-a  is the direction-dependent right-hand circularly polarized (RHCP) power component incident on an isotropic antenna; the left-handed circularly polarized (LHCP) component is negligible. The direct antenna gain, Eq-6-b, is obtained evaluating the antenna pattern in the satellite direction and with RHCP polarization.

The reflection power,

Eq-7, (7)

is defined starting with the same incident isotropic power, Eq-6-a, as in the direct power. It ends with a coherent power attenuation factor, 

Eq-8  (8)

where  θ  is the angle of incidence (with respect to the surface normal), k = 2π/λ, is the wave number, and λ = 24.4 centimeters is the carrier wavelength for the civilian GPS signal on the L2 frequency (L2C). This polarization-independent factor accounts only for small-scale residual height above and below a large-scale trend surface. The former/latter results from high-/low-pass filtering the actual surface heights using the first Fresnel zone as a convolution kernel, roughly speaking. Small-scale roughness is parameterized in terms of an effective surface standard deviation s (in meters); its scattering response is modeled based on the theories of random surfaces, except that the theoretical ensemble average is replaced by a sensing spatial average. Large-scale deterministic undulations could be modeled, but their impact on snow depth is canceled to first-order by removing bare-ground reflector heights.

At the core of Eq-Pr, we have coupled surface/antenna reflection coefficients,  Eq-X=, producing respectively RHCP and LHCP fields (under the assumption of a RHCP incident field). These terms include antenna response power gain and phase patterns, evaluated in the reflection direction, and separately for each polarization. The surface response is represented by complex-valued Fresnel coefficients for cross- and same-sense circular polarization, respectively. The medium is assumed to be homogeneous (that is, a semi-infinite half-space). Material models provide the complex permittivity, which drives the Fresnel coefficients.

The interferometric phase reads:


The first term accounts for the surface and antenna properties of the reflection, as above. The last one is the direct phase contribution, which amounts to only the RHCP antenna phase-center variation evaluated in the satellite direction. The majority of the components present in the direct RHCP phase (such as receiver and satellite clock states, the bulk of atmospheric propagation delays, and so on) are also present in the reflection phase, so they cancel out in forming the difference.

At the core of the interferometric phase, we have the geometric component, φI = i, the product of the wave number and the interferometric propagation delay. Assuming a locally horizontal surface, the latter is simply:

Eq-10  (10)

in terms of the satellite elevation angle, e, and an a priori reflector height, HA. Snow depth will be measured in terms of changes in reflector height.

The physical forward model, based only on a priori information, can then be summarized as:

Eq-11a  (11)

where interferometric power and phase are, respectively:

 Eq-12 (12)

Eq-13. (13)

In all of these terms the pseudorandom-noise-code modulation impressed on the carrier wave can be safely neglected, given the small interferometric delay and Doppler shift at grazing incidence, stationary surface/receiver conditions, and short antenna installations.

Parameterization of Unknowns

There are errors in the nominal values assumed for the physical parameters of the model (permittivity, surface roughness, reflector height, and so on). Ideally we would estimate separate corrections for each one, but unfortunately many are linearly dependent or nearly so. Because of this dependency, we have kept physical parameters fixed to their optimal a priori values, and have estimated a few biases. Each bias is an amalgamation of corrections for different physical effects. In a later stage, we rely on multiple independent bias estimates (such as for successive days) to try and separate the physical sources.

Each satellite track is inverted independently. A track is defined by partitioning the data by individual satellite and then into ascending and descending portions, splitting the period between the satellite’s rise and set at the near-zenith culmination. Each satellite track has a duration of ~1–2 hours. This configuration normally offers a sufficient range of elevation angles, unless the satellite reaches culmination too low in the sky (less than about 20°), in which case the track is discarded. In seeking a balance between under- and over-fitting, between an insufficient and an excessive number of parameters, we estimate the following vector of unknown parameters:

Eq-14. (14)

FIGURE 2 shows the effect of the constant and linear biases on the SNR observations. Reflector height bias, HB , changes the number of oscillations; phase shift, φB , displaces the oscillations along the horizontal axis; reflection power, Eq-14-a   , affects the depth of fades; zeroth-order noise power,   Eq-14-b  , shifts the observations up or down as a whole; and first-order noise power,  Eq-14-c  , tilts the SNR curve. A good parameterization yields observation sensitivity curves as unique as possible for each parameter.

FIGURE 2. Effect of each parameter on SNR observations; curves are displaced vertically (6 dB) for clarity.

FIGURE 2. Effect of each parameter on SNR observations; curves are displaced vertically (6 dB) for clarity.

The forward model, now including the biases, can be summarized as follows:

Eq-15 (15)

where the modified interferometric power and phase are given by:

Eq-16, (16)

Eq-17. (17)

The total reflector height, H = HAHB (a priori value minus unknown bias), is to be interpreted as an effective value that best fits measurements, which includes snow and other components.

Bootstrapping Parameter Priors. Biases and SNR observations are involved non-linearly through the forward model. Therefore, there is the need for a preliminary global optimization, without which the subsequent final local optimization will not necessarily converge to the optimal solution.

SNR observations would trace out a perfect sinusoid curve in the case of an antenna with isotropic gain and spherical phase pattern, surrounded by a smooth, horizontal, and infinite surface (free of small-scale roughness, large-scale undulations, and edges), made of perfectly electrically conducting material, and illuminated by constant incident power. Thus, in such an idealized case, SNR could be described exactly by constant reflector height, phase shift, amplitude, and mean values.

As the measurement conditions become more complicated, the SNR data start to deviate from a pure sinusoid. Yet a polynomial/spectral decomposition is often adequate for bootstrapping purposes. 

Statistical Inverse Model Formulation

Based on the preliminary values for the unknown parameters vector and other known (or assumed) values, we run the forward model to obtain simulated observations. We form pre-fit residuals comparing the model values to SNR measurements collected at varying satellite elevation angles (separately for each track). Residuals serve to retrieve parameter corrections, such that the sum of squared post-fit residuals is minimized. This non-linear least squares problem is solved iteratively using both a functional model and a stochastic model. The functional modeling includes a Jacobian matrix of partial derivatives, which represents the sensitivity of observations to parameter changes where the partial derivatives are defined element-wise. Instead of deriving analytical expressions, we evaluate them numerically, via finite differencing. The stochastic model specifies the uncertainty and correlation expected in the residuals. Their a priori covariance matrix modifies the objective function being minimized. 

Directional Dependence

It is important to know at which elevation angles the parameter estimates are best determined. Here, we focus on the phase parameters instead of reflection power or noise power parameters. 

We can utilize the estimated reflector height and phase shift to evaluate the full phase bias function over varying elevation angles. Similarly, we can extract the corresponding 2-by-2 portion of the parameters’ a posteriori covariance matrix, containing the uncertainty for reflector height and for phase shift, as well as their correlation, which is then propagated to obtain the full phase uncertainty (see FIGURE 3).

FIGURE 3. Uncertainty of full phase function, propagated from the uncertainty of reflector height and of phase shift, as well as their correlation.

FIGURE 3. Uncertainty of full phase function, propagated from the uncertainty of reflector height and of phase shift, as well as their correlation.

The uncertainty attains a clear minimum versus elevation angle. The least-uncertainty elevation angle pinpoints the observation direction where reflector height and phase shift are best determined (in combined form, not individually). The azimuth and epoch coinciding with the peak elevation angle act as track tags, later used for clustering similar tracks and analyzing their time series of retrievals.

If we normalize phase uncertainty by its value at the peak elevation angle, then plot such sensing weights (between 0 and 1) versus the radial or horizontal distance to the center of the first Fresnel zone at each elevation angle, we obtain FIGURE 4. It can be interpreted as the reflection footprint, indicating the importance of varying distances, with a longer far tail and a shorter near tail (respectively regions beyond and closer than the peak distance). The implications for in situ data collection are clear: one should sample more intensely near the peak distance (about 15 meters) and less so in the immediate vicinity of the GPS antenna, tapering it off gradually away from the antenna. As a caveat, these conclusions are not necessarily valid for antenna setups other than the one considered here.

FIGURE 4. Reflection footprint in terms of a sensing weight (between 0 and 1) defined as the normalized reciprocal of full phase uncertainty, plotted versus the radial or horizontal distance from the receiving antenna to the center of the first Fresnel zone at each elevation angle; valid for an upright 2-meter-tall antenna; the receiving antenna is at zero radial distance.

FIGURE 4. Reflection footprint in terms of a sensing weight (between 0 and 1) defined as the normalized reciprocal of full phase uncertainty, plotted versus the radial or horizontal distance from the receiving antenna to the center of the first Fresnel zone at each elevation angle; valid for an upright 2-meter-tall antenna; the receiving antenna is at zero radial distance.


We now examine the snow-depth retrievals from the GPS multipath retrieval algorithm and assess both the precision and accuracy of the method. Multiple metrics have been developed to assess the quality of the results. The accuracy of the method has been evaluated by comparing with in situ data over a multi-year period. Three field sites were chosen to highlight different limitations in the method, both in terms of terrain and forest cover: grassland, alpine, and forested. We will look at the forested site in some detail.

Satellite Coverage and Track Clustering. All GPS-MR retrievals reported here are based on the newer GPS L2C signal. Of the approximately 30 GPS satellites in service, 8-10 L2C satellites were available between 2009 and 2012 (8, 9, and 10 satellites at the end of 2009, 2010, and 2011, respectively). Satellite observations were partitioned into ascending and descending portions, yielding approximately twenty unique tracks per day at a site with good sky visibility. GPS orbits are highly repeatable in azimuth, with deviations at the few-degree range over a year, translating into ~50-100-centimeter azimuthal displacement of the reflecting area (corresponding to the first Fresnel zone at 10°-15° elevation angle for a 2-meter high antenna). This repeatability permits clustering daily retrievals by azimuth. It also allows the simplification that estimated snow-free reflector heights are fairly consistent from day to day, facilitating the isolation of the varying snow depth during the snow-covered period.

For a given track, its revisit time is also repeatable, amounting to practically one sidereal day. The deficit in time relative to a calendar day results in the track time of the day receding ~4 minutes and 6 seconds every day. This slow but steady accumulation eventually makes the time of day return to its starting value after about one year. As all GPS satellites drift approximately at the same rate, the time between successive tracks remains nearly repeatable. Its reciprocal, the sampling rate, has a median equal to approximately one track per hour, with a low value of one track within two hours and a high of one track within 15 minutes; both extremes occur every day, with low-rate idle periods interspersed with high-rate bursts. The time of the day reduced to a fixed day (such as January 1, 2000) could also be used to cluster tracks. Neighboring clusters, which are close in azimuth and/or in reduced time of the day, are expected to be more comparable, as they sample similar conditions and are subject to similar errors.

Observations. FIGURE 5 shows several representative examples of SNR observations. A typical good fit between measured and modeled values is shown in Figure 5(a), corresponding to the beginning of the snow season. Generally the model/measurement fit is good when the scattering medium is homogeneous; it deteriorates as the medium becomes more heterogeneous, particularly with mixtures of soil, snow, and vegetation. There are genuine physical effects as well as more mundane spurious instrumental issues that degrade the fit but do not necessarily cause a bias in snow-depth estimates. These include secondary reflections, interferometric power effects, direct power effects, and instrument-related issues.

FIGURE 5. Examples of observations: (a) good fit; (b) presence of secondary reflections; (c) vanishing interference fringes; (d) atypical interference fringes.

FIGURE 5. Examples of observations: (a) good fit; (b) presence of secondary reflections; (c) vanishing interference fringes; (d) atypical interference fringes.

Secondary reflections originate from disjoint surface regions. Interference fringes become convoluted with multiple superimposed beats (see Figure 5(b)). As long as there is a unique dominating reflection, the inversion will have no difficulty fitting it, as the extra reflections will remain approximately zero-mean.

Random deviations of the actual surface with respect to its undulated approximation, called roughness or residual surface height, will affect the interferometric power. SNR measurements will exhibit a diminishing number of significant interference fringes, compared to the measurement noise level (see Figure 5(c)). This facilitates the model fit but the reflector height parameter may become ill-determined: its estimates will be more uncertain. Changes in snow density also affect the fringe amplitude.

Snow precipitation attenuates the satellite-to-ground radio link, which affects SNR measurements through the direct power term. First, this shifts the SNR measurements up or down (in decibels); second, it tilts the trend tSNR as attenuation is elevation-angle dependent; third, fringes in dSNR will change in amplitude because of the decrease in the coherent component of the direct power.

Partial obstructions can affect either or both direct and interferometric powers. In this case, SNR measurements, albeit corrupted, are still recorded. This situation is in contrast to complete blockages as caused by topography. The deposition of snow and the formation of a winter rime on the antenna are a particularly insidious type of obstruction, as their presence in the near-field of the antenna element can easily distort the gain pattern in a significant manner. In the far-field, trees are another important nuisance, so much so that their absence is held as a strong requirement for the proper functioning of multipath reflectometry.

Satellite-specific direct power offsets and also long-term power drifts are to be expected as spacecraft age and modernized designs are launched. In addition, noise power depends on the state of conservation of receiver cables and on their physical temperature. Less subtle incidents are sudden ~3-dB SNR steps, hypothesized to originate in the receiver switching between the L2C data and pilot subcodes, CM and CL.

Quality Control. Anomalous conditions may result in measurement spikes, jumps, and short-lived rapidly-varying fluctuations. For snow-depth-sensing purposes, it is necessary and sufficient to either neutralize such measurement outliers through a statistically robust fit or detect unreliable fits and discard the problematic ones that could not otherwise be salvaged.

The key to quality control (QC) is in grouping results into statistically homogeneous units, having measurements collected under comparable conditions. In our case, azimuth-clustered tracks are the natural starting unit. Secondarily, we must account for genuine temporal variations in the tendency of results, from beginning to peak to the end of the snow season. The detection of anomalous results further requires an estimate of the statistical dispersion to be expected. Considering that the sample is contaminated with outliers, robust estimators (running median instead of the running mean, and median absolute deviation over the standard deviation) are called for, if the first- and second-order statistical moments are to be representative. Given estimates of the non-stationary tendency and dispersion, a tolerance interval can then be constructed such that it bounds, say, a 99% proportion of the valid results with 95% confidence level. We also desire QC to be judicious, or else too many valid estimates will be lost. Notice that in the present intra-cluster QC, we compare an individual estimate to the expected performance of the track cluster to which it belongs; later, we complement QC with an inter-cluster comparison of each cluster’s own expected performance.

Based on our practical experience, no single statistic detects all the outliers. We use four particular statistics that we have found to be useful: 1) degrees of freedom, essentially the number of observations per track (modulo a constant number of parameters); 2) using the scaled root-mean-square error (RMSE) to test for goodness-of-fit, that is, how well measurements can be explained adjusting the unknown values for the parameters postulated in the model; 3) reflector height uncertainty; and 4) peak elevation angle, which behaves much like a random variable, as it is determined by a multitude of factors. 

Combinations. We combine multiple clusters to average out random noise. Noise mitigation aims at not only coping with measurement errors but also compensating for model deficiencies, to the extent that they are not in common across different clusters. Before we combine different clusters, we have to address their long-term differences. The initial situation is that snow surface heights will be greater downhill and smaller uphill; we take this into account on a cluster-by-cluster basis by subtracting ground heights from their respective snow surface heights, resulting in snow thickness values, which is a completely physically unambiguous quantity. Snow thickness is more comparable than snow heights across varying-azimuth track clusters. Yet snow tends to fill in ground depressions, so thickness exhibits variability caused by the underlying ground surface, even when the overlying snow surface is relatively uniform. Further cluster homogeneity can be achieved by accounting for the temporally permanent though spatially non-uniform component of snow thickness. 

The averaging of snow depths collected for different track clusters employs the inversion uncertainties to obtain a preliminary running weighted median, calculated for, say, daily postings, with overlapping windows or not. The preliminary post-fit residuals then go through their own averaging, necessarily employing a wider averaging window (say, monthly), which produces scaling factors for the original uncertainties. The running weighted median is then repeated, producing final averages. The variance factors reflect the fact that some clusters are better than others.

Thus, the final GPS estimates of snow depth follow from an averaging of all available tracks, whose individual snow depth values were previously estimated independently. A new average is produced twice daily utilizing the surrounding 1–2 days of data (depending on the data density), that is, 12-hour posting spacing and 24-hour moving window width. The averaging interval must be an integer number of days, so as to minimize the possibility of snow-depth artifacts caused by variations in the observation geometry, which repeats daily.

Site-Specific Results

We explored GPS-MR snow-depth retrieval at three stations over a long period (up to three years). Throughout, we assessed the performance of the GPS estimates against independent nearly co-located in situ measurements. We also compared the GPS estimates to the nearest SNOTEL station. SNOTEL (from snowpack telemetry) is an automated system for collecting snowpack and related data in the western U.S. operated by the U.S. Department of Agriculture. Although not co-located with GPS, SNOTEL data are important because they provide accurate information on the timing of snowfall events.

The three sites we used were 1) a site in the T.W. Daniel Experimental Forest within the Wasatch Cache National Forest in the Bear River Range of northeastern Utah, with an elevation of 2,600 meters; 2) one of the stations of the EarthScope Plate Boundary Observatory, a grassland site located near Island Park, Idaho; and 3) an alpine site in the Niwot Ridge Long-term Ecological Research Site near Boulder, Colorado. While we have fully documented the results from each site, due to space limitations we will only discuss the results from the forested site (known as RN86) in this article. This is a more challenging site than the other two, due to the presence of nearby trees. Furthermore, it was subject to denser in situ sampling of 20-150 measurements spatially replicated around the GPS antenna, and repeated approximately every other week for about one year.

We show results for the 2012 water-year, the period starting October 1 through September 30 of the following year. Where GPS site RN86 was installed, topographical slopes range from 2.5° to 6.5° (at the 2-meter spatial scale), with average of ~5° within a 50-meter radius around the GPS antenna. RN86 was specifically built to study the impact of trees on GPS snow depth retrievals (see FIGURE 6). Ground crews manually collected in situ measurements around the GPS antenna approximately every other week starting in November 2011. Measurements were made every 1–2 meters from the antenna up to a distance of 25-30 meters. In the second half of the year, the sampling protocol was changed to azimuths of 0° (N), 45° (NE), 135° (SE), 180° (S), 225° (SW), and 315° (NW). With these data it is possible to obtain in situ average estimates, with their own uncertainties (based on the number of measurements), which allows a more meaningful comparison.

FIGURE 6. Aerial view of the forested site (RN86) around the GPS antenna (marked with a circle).

FIGURE 6. Aerial view of the forested site (RN86) around the GPS antenna (marked with a circle).

There is reduced visibility at the current site, compared to other sites. Track clusters are concentrated due south, with only two clusters located within ±90° of north. Therefore, the GPS average snow depth is not necessarily representative of the azimuthally symmetric component of the snow depth. In the presence of an azimuthal asymmetry in the snow distribution around the antenna, the GPS average would be expected to be biased towards the environmental conditions prevalent in the southern quadrant. To rule out the possibility of an azimuthal artifact in the comparisons, we have utilized only the in situ data collected along the SE/S/SW quadrant.

The comparison shows generally excellent agreement between GPS and in situ data (see FIGURE 7). The first four and the last one in situ data points were collected with coarser spacing and/or smaller azimuthal coverage, which may be partially responsible for different performance in the first and second halves of the snow season. The correlation between GPS and in situ snow depth at RN86 amounts to 0.990, indicating a very strong linear relationship. Carrying out a regression between in situ and GPS values, the RMS of snow-depth residuals improves from 9.6 to 3.4 centimeters. The regression intercept and slope (with corresponding 95% uncertainties) amount to 15.4 ± 9.11 centimeters and 0.858 ± 0.09 meters per meter, respectively. According to these statistics, the null hypotheses of zero intercept and unity slope are rejected at the 95% confidence level. This implies that at this location GPS snow-depth estimates exhibit both additive and multiplicative biases. The latter is proportional to snow depth itself, meaning that, compared to an ideal one-to-one relationship, GPS is found to under-estimate in situ snow depth at this site by 14 ± 9%, although the uncertainty is somewhat large.

FIGURE 7. Snow-depth measurement at the forested site (RN86) for the water-year 2012

FIGURE 7. Snow-depth measurement at the forested site (RN86) for the water-year 2012

The SNOTEL sensors are exceptionally close to the GPS antenna at this site, about 350 meters horizontally distant with negligible vertical separation. Yet the former is located within trees, while the latter is located at the periphery of the forest and senses the reflections scattered from an open field. Therefore, only the timing of snowfall events agrees well, not the amount of snow. Although forest density is generally negatively correlated with snow depth, exceptions are not uncommon, especially in localized clearings exposed to intense solar radiation, where shading of the snow by the trees reduces ablation.


In this article, we have discussed a physically based forward model and a statistical inverse model for estimating snow depth based on GPS multipath observed in SNR measurements. We assessed model performance against independent in situ measurements and found they validated the GPS estimates to within the limitations of both GPS and in situ measurement errors after the characterization of systematic errors. The assessment yielded a correlation of 0.98 and an RMS error of 6–8 centimeters for observed snow depths of up to 2.5 meters at three sites, with the GPS underestimating in situ snow depth by ~5–15%. This latter finding highlights the necessity to assess effects currently neglected or requiring more precise modeling.


The research reported in this article was supported by grants from the U.S. National Science Foundation, NASA, and the University of Colorado. Nievinski has been supported by a Capes/Fulbright Graduate Student Fellowship and a NASA Earth System Science Research Fellowship. The article is based, in part, on two papers published in the IEEE Transactions on Geoscience and Remote Sensing: “Inverse Modeling of GPS Multipath for Snow Depth Estimation – Part I: Formulation and Simulations” and “Inverse Modeling of GPS Multipath for Snow Depth Estimation – Part II: Application and Validation.”


For the forested site (RN86), a Trimble NetR9 receiver was used with a Trimble TRM57971.00 (Zephyr Geodetic II) antenna with no external radome.

FELIPE G. NIEVINSKI is a faculty member at the Federal University of Santa Catarina, Florianópolis, Brazil. He has also been a post-doctoral researcher at São Paulo State University, Presidente Prudente, Brazil. He earned a B.E. in geomatics from the Federal University of Rio Grande do Sul, Porto Alegre, Brazil, in 2005; an M.Sc.E. in geodesy from the University of New Brunswick, Fredericton, Canada, in 2009; and a Ph.D. in aerospace engineering sciences from the University of Colorado, Boulder, in 2013. His Ph.D. dissertation was awarded The Institute of Navigation Bradford W. Parkinson Award in 2013.

KRISTINE M. LARSON received a B.A. degree in engineering sciences from Harvard University and a Ph.D. degree in geophysics from the Scripps Institution of Oceanography, University of California at San Diego. She was a member of the technical staff at the Jet Propulsion Lab from 1988 to 1990. Since 1990, she has been a professor in the Department of Aerospace Engineering Sciences, University of Colorado, Boulder.


• Authors’ Journal Papers

“Inverse Modeling of GPS Multipath for Snow Depth Estimation—Part I: Formulation and Simulations” by F.G. Nievinski and K.M. Larson in IEEE Transactions on Geoscience and Remote Sensing, Vol. 52, No. 10, 2014, pp. 6555–6563, doi: 10.1109/TGRS.2013.2297681.

“Inverse Modeling of GPS Multipath for Snow Depth Estimation—Part II: Application and Validation” by F.G. Nievinski and K.M. Larson in IEEE Transactions on Geoscience and Remote Sensing, Vol. 52, No. 10, 2014, pp. 6564–6573, doi: 10.1109/TGRS.2013.2297688.

• More on the Use of GPS for Snow Depth Assessment

“Snow Depth, Density, and SWE Estimates Derived from GPS Reflection Data: Validation in the Western U.S.” by J.L. McCreight, E.E. Small, and K.M. Larson in Water Resources Research, published first on line, August 25, 2014, doi: 10.1002/2014WR015561.

Environmental Sensing: A Revolution in GNSS Applications” by K.M. Larson, E.E. Small, J.J. Braun, and V.U. Zavorotny in Inside GNSS, Vol. 9, No. 4, July/August 2014, pp. 36–46.

Snow Depth Sensing Using the GPS L2C Signal with a Dipole Antenna” by Q. Chen, D. Won, and D.M. Akos in EURASIP Journal on Advances in Signal Processing, Special Issue on GNSS Remote Sensing, Vol. 2014, Article No. 106, 2014, doi: 10.1186/1687-6180-2014-106.

“GPS Snow Sensing: Results from the EarthScope Plate Boundary Observatory” by K.M. Larson and F.G. Nievinski in GPS Solutions, Vol. 17, No. 1, 2013, pp. 41–52, doi: 10.1007/s10291-012-0259-7.

• GPS Multipath Modeling and Simulation

“Forward Modeling of GPS Multipath for Near-Surface Reflectometry and Positioning Applications” by F.G. Nievinski and K.M. Larson in GPS Solutions, Vol. 18, No. 2, 2014, pp. 309–322, doi: 10.1007/s10291-013-0331-y.

“An Open Source GPS Multipath Simulator in Matlab/Octave” by F.G. Nievinski and K.M. Larson in GPS Solutions, Vol. 18, No. 3, 2014, pp. 473–481, doi: 10.1007/s10291-014-0370-z.

Multipath Minimization Method: Mitigation Through Adaptive Filtering for Machine Automation Applications” by L. Serrano, D. Kim, and R.B. Langley in GPS World, Vol. 22, No. 7, July 2011, pp. 42–48.

It’s Not All Bad: Understanding and Using GNSS Multipath” by A. Bilich and K.M. Larson in GPS World, Vol. 20, No. 10, October 2009, pp. 31–39.

GPS Signal Multipath: A Software Simulator” by S.H. Byun, G.A. Hajj, and L.W. Young in GPS World, Vol. 13, No. 7, July 2002, pp. 40–49.

This article is tagged with , and posted in Innovation, OEM

About the Author: Richard B. Langley

Richard B. Langley is a professor in the Department of Geodesy and Geomatics Engineering at the University of New Brunswick (UNB) in Fredericton, Canada, where he has been teaching and conducting research since 1981. He has a B.Sc. in applied physics from the University of Waterloo and a Ph.D. in experimental space science from York University, Toronto. He spent two years at MIT as a postdoctoral fellow, researching geodetic applications of lunar laser ranging and VLBI. For work in VLBI, he shared two NASA Group Achievement Awards. Professor Langley has worked extensively with the Global Positioning System. He has been active in the development of GPS error models since the early 1980s and is a co-author of the venerable “Guide to GPS Positioning” and a columnist and contributing editor of GPS World magazine. His research team is currently working on a number of GPS-related projects, including the study of atmospheric effects on wide-area augmentation systems, the adaptation of techniques for spaceborne GPS, and the development of GPS-based systems for machine control and deformation monitoring. Professor Langley is a collaborator in UNB’s Canadian High Arctic Ionospheric Network project and is the principal investigator for the GPS instrument on the Canadian CASSIOPE research satellite now in orbit. Professor Langley is a fellow of The Institute of Navigation (ION), the Royal Institute of Navigation, and the International Association of Geodesy. He shared the ION 2003 Burka Award with Don Kim and received the ION’s Johannes Kepler Award in 2007.