Snow Layer Detection by Pattern Matching in a Multipath Radar Interference Scenario

ABSTRACT Snowpack in sloped mountain terrain may present stratification in which layers of snow accumulate over a certain period. The loosely adhering layers in snowpack may exhibit a severe avalanche risk, and detecting such stratification phenomena is key to avalanche detection. The danger of manual, in-situ sampling calls for investigation of alternative approaches, including the use of in-situ radars. Considering a buried radar sensor, this manuscript analyses the impact of multipath phenomena on the results of the sensing operation. Multiple bounces in snow layers may degrade the received signal and corrupt the estimation results. Multipath phenomena are well known and analysed in communication systems, but snow layering measurement is a different problem, with little previous research done on related multipath phenomena. This paper analyses the problem with a multi-receiver radar architecture, investigating the impact of multi-path effects from layered snow on measurements carried out by a multi-static radar of Frequency-Modulated Continuous Wave (FMCW) type. A few simplifying assumptions are formulated, in order to make the problem manageable, and a possible approach for solving the related issues is proposed, using a robust algorithm based on pattern matching to detect snow layering in such a challenging environment. Results obtained from various simulations carried out with realistic environmental and geometrical parameters appear to support the effectiveness of the proposed method. Results show that with 3 metres of snow thickness, the probability of successfully detecting the correct numbers of layers exceeds 80%.


Introduction
When avalanche risk is concerned, the most important factor impacting the stability of the snowpack is the stratification of snow. Stratification refers to the presence of snow or ice layers with different characteristics, possibly loosely adhering to each other. This situation is due to thermodynamic changes within the snowpack, and frequently connected with freshly fallen snow (see Quno et al. (2018)). Snowpack is a non-homogeneous and thermodynamically very active material; therefore, it reacts promptly to any environmental or temporal variation. Density layering emerges as a result of overburden pressure and heat transport within the snowpack. Between successive snow deposits, the specific temperature and humidity conditions uniquely transform the deposited snow, which will result in different features with respect to the next layer. As a result, in the long term, a snowpack ends up being made of multiple layers, whose mutual interfaces may fail to offer a steady grip, especially when the two adjacent layers differ strongly. Avalanches, for example, occur when the underlying layer can not withstand the load of the overlying layer. Consequently, the underlying layer collapses under the pressure; this results in snow mass slide, also known as avalanche. Other mechanisms may activate, but in many cases, the presence of layers with different thickness and density is an important factor. Snow thickness estimation has already been investigated in scientific literature (Botteron et al. (2013); Eriksrd et al. (2019); Heilig et al. (2009a); Tsai et al. (2019)), considering different sensing devices. While estimating the snow thickness with microwave sensing devices like an FMCW (Frequency Modulated Continuous Wave) radar, it is important to determine the snowpack stratigraphy, as it is definitely relevant to avalanche forecasting (Myhre 2018). In addition to that, for accurate estimation, it is essential that practical scenarios like multipaths in snow medium be considered. These multipaths are basically present because of the non-homogeneity of the medium.
Microwave systems are a natural choice for monitoring the snowpack because they can provide non-destructive and remote monitoring almost instantaneously. Therefore, they represent a valid solution to overcome at least a part of the inconveniences related to manual analysis. Systems based on airborne or satellite platforms, both active and passive radar, have been exploited by Ulaby (1982). A major advantage is their ability to cover large areas, but in most cases, they do not represent the best option for continuous snow monitoring in real time. Moreover, the achievable resolution of the internal stratigraphy of snow is far behind the potential of ground-based systems, and the related computational and system costs could be high. Ground-based radar systems with quasi-mono-static or bi-static configuration are suited for remotely controlled applications with continuous monitoring of the snowpack. They can provide adequate resolution, and when buried they represent the most convenient solution to provide adequate protection in the avalanche release zones. Although installation and maintenance of such systems could be complicated on the steep sides of a mountain, the performance gain in terms of accuracy and computational complexity still makes ground-based radar the best choice for snowpack monitoring. Upward-looking radars have received significant attention due to the possibility of buried installation, as previously discussed. Buried radars   Eriksrd et al. (2019); Okorn et al. (2014)) are used to determine some properties of snow, i.e. height or density; however, either they require some additional external aids or devices, i.e. ultrasound and laser metres, or they rely on a priori assumptions. The need for additional external aid has been considered a critical limiting factor for operational installations. Besides that, external aid devices like ultrasonic and laser metres can not be employed in the buried installation, i.e. upward-looking radar. In addition to that, because of additional devices, the overall system to estimate the snowcover becomes more complicated. An innovative radar measurement system based on multi-receiver architecture was presented by Espín-López et al. (2018); Kulsoom, Dell'Acqua, and Pasian (2018); Pasian et al. (2019) using FMCW upward-looking radar. The system was based on FMCW radars with two independent receivers, which can simultaneously calculate the total snow thickness and wave speed into the medium. This multi-receiver architecture is theoretically capable of determining snow pack thickness and wave speed into medium simultaneously, without any external aid or apriori assumption. The location of the receiver, however, may in principle heavily impact on the measurement accuracy. The present work is based on the same system architecture as that of our previous work (Espín-López et al. (2018); Kulsoom, Dell'Acqua, and Pasian (2018); Pasian et al. (2019)), which discussed the problem of snowpack estimation with buried FMCW bi-static radar made up of multi-receiver design. In this paper, however, we will focus on the problem of snow layer detection in the presence of multipath interference. Multipath interference issue is observed also in radio communication systems; we have attempted and proposed enhancements leveraging on the solutions that were put forward in a communication scenario. Because of successive snowfalls and temperature changes, the snowpack becomes non-homogeneous and consequently gets stratified, as already mentioned. Unavoidably, reflections and multiple paths are present in the system, which generate disturbances and deteriorate the estimation process. Although in some cases multipath propagation can be useful to collect additional information as pointed out by Espín-López and Pasian (2020), in general, it constitutes a source of disturbances. The presence of multipath interference is common among various technologies, including wireless communications, radar tracking, and positioning. The multi-path interference problem is recognized and well investigated in flight tracking radars. It has been proven that because of multipath interference the system struggles in tracking targets at low elevation angles (Mrstik and Smith 1978). A significant research effort has been spent to investigate and mitigate or suppress the effects of multi-path interference in radars and communication systems. The research work reported by Colone et al. (2009b), Ducheng Wu et al. (2012), Cardinali et al. (2007), Colone et al. (2009a), Wang and Shao (2014) have employed least mean square (LMS), constant modulus algorithm (CMA), or some other adaptive filters and attempted to cancel or suppress multipath interference in radar considering either a Lineof-Sight (LOS) reference antenna or reference signal. Besides their remarkable complexity, the above-mentioned algorithms cannot be employed in cases where no LOS antenna is available. In our case, multipaths are originated inside the snowpack due to the different layer reflections, and the LOS antenna is not a feasible solution in snowpack sensing. In the work by Schonken, de Swardt, and van der Merwe (2016), the deterioration in FMCW signal amplitude and phase due to multipath interference is discussed. Although the authors have set up a very good model for the system and deeply investigated the problem, they do not propose a specific approach to mitigating it. The research work by Roehr, Gulden, and Vossiek (2007) suggests to increase the bandwidth to tackle the problem, which in turn improves the resolution. This is expected to enable discrimination of the actual target in the received spectrum, even with the superimposition of timeshifted versions of the same signal, due to multipath reflection. This solution, however, is only suitable for those cases where a substantial increase in bandwidth is viable. To reduce the effect of multipaths interference, in the work by Wagner et al. (2011), the FMCW radar used a detector and positioning error was analysed for different amplitude ratios between LOS and No-Line-of-Sight (NLOS) signal. A trade-off was found between positioning error and windowed functions. Despite the reduction in the power level of side lobes due to the window functions, the resolution becomes somehow poor because of the bandwidth decrease.
The problem of multipath interference is also discussed by Lucas et al. (2017) for ground-based, upward-looking radar. The authors avoid multipath interference either by setting up an appropriate measurement geometry or through the shielding of the refractive surface of antennas. This approach reportedly mitigates interference to a significant extent. It is however not designed to deal multipath reflection originated from inside a snowpack, which is a different problem that cannot be mitigated by just shielding or using some particular antenna geometry. Then, in the work by Gurgel, Barbin, and Schlick (2007), an algorithm is proposed to suppress interference in FMCW radar, but starting from an assumption that the peaks in the spectrum must be separated from each other and there should be no complete overlap in peaks. However, this assumption is unrealistic in the case of multiple layer detection due to the short distances involved in snow layers. Based on the findings of all the previous works referred above, in this paper, we propose a novel method to detect multiple snow layers, in the presence of multipath interference and also to estimate their thickness. To the best of our knowledge, we are the first to investigate this issue in snowpack monitoring based on an FMCW, upward-looking radar.

System description
Typical dry snow densities for mountain environments range from around 90 kg m −3 to around 450 kg m −3 (Hallikainen, Ulaby, and Abdelrazik 1986), with a negligible Liquid Water Content (LWC). On the other hand, typical wet snow densities' range is 350 kg m −3 to 550 kg m −3 , with values changing during the day with changes in climatic variables (Domine (2011);Seibert et al. (2015)). LWCs for wet snow are often evaluated up to around 10−12% of the volume, since values in excess of 12% are usually no more relevant to the instability of the snowpack (Heilig, Schneebeli, and Eisen 2009b). A notable exception to the above dry and wet cases is that of ice-like layers, with a density approaching 800 kg m −3 to 900 kg m −3 and a negligible LWC. Table 1 summarizes the physical and dielectric properties for a snowpack.
For efficient snowpack stratigraphy estimation, all practical scenarios including multiple paths at the receiving radar must be considered. In this context, when a radar transmits an electromagnetic wave, the echo signal undergoes multiple-bounce reflections because of uneven snow surface and stratigraphy. The focus of this work is to investigate the effects of multipath interference because of indirect reflections which add to direct specular reflection in the snowpack. Although we have examined the snowpack up to 10 layers, for detailed analysis and clear visibility of the figure, only three layers with Table 1. Dielectric properties of snow at 1 GHz calculated according to Hallikainen, Ulaby, and Abdelrazik (1986 different dielectric properties and densities have been demonstrated. The purpose of this paper is to propose a method rather than to include all possible cases.
When dealing with such multi-path scenarios, the spectrum of the echoed mixed signal will no longer contain a single peak but rather one peak per significant reflection. Because of the aforementioned facts, each receiver will collect multiple copies of the signal, each one being a delayed version of the actual signal. Delays cannot be determined apriori and these indirect multiple reflections act defacto as added noise. As such, they affect the estimation accuracy, and may ultimately result in false detection if the algorithm is not robust enough. The concept is illustrated on the specific example of an FMCW radar with a central frequency f of 2.75 GHz and a bandwidth B of 1.10 GHz.

FMCW basics
In an FMCW radar, several types of modulation schemes can be used, e.g. saw-tooth, triangular, and sinusoidal. In this work, we focus on the saw-tooth modulation, as it offers the advantage of better detecting smaller and adjacent targets due to the finer resolution achievable with respect to other techniques (Winarko et al. (2017) Stove (1992). With this type of radar, the transmitted frequency increases linearly as a function of time during a period called sweep time to then restart from the initial frequency again. Typically, the signal is continuously transmitted, with linearly varying frequency. Time-frequency analysis of FMCW-based transmitted signal with sawtooth modulation is shown in Figure 1. The frequency of the signal during one sweep cycle can be expressed using the following signal parameters: sweep frequency F sweep or sweep bandwidth B sweep and sweep time T sweep , together with the starting frequency f 0 T sweep is termed chirp rate and can be considered as the speed of the frequency change. The variations in frequency will also effect the phase as shown in the formula below: Therefore, the transmitted signal in the first sweep can be written as By integrating the above equation, we have where A T x and ϕ 0 are the amplitude and phase of the transmitted signal, respectively. After reflection from each layer, an echoed signal will be delayed by τ L : where τ L is the round trip time; if the targeted ith layer is at a distance L i then τ L i ' 2L i c ffi ffi ffi ffi ffi ffi ffi ε i;m p , where c is the speed of light and ffi ffi ffi ffi ffi ffi ffi ε i;m p is the mean dielectric permittivity of layers up to the ith layer. Because of possible multiple layers in snowpack, more than one signal reflections are expected; each reflected signal is assumed to undergo changes in phase and amplitude because of channel effects. According to the FMCW radar principle, the returned signal is demodulated by mixing it with the transmitted signal. This mixing will produce some higher frequency components, which will be filtered out.
In the above equation, A m is the amplitude of the mixed signal. 1 From Equation 7 it can be seen that the signal frequency changes as a consequence of filtering and mixing operation. The resulting frequency is usually named beat frequency, which is the difference in frequency between the transmitted and the received signals. This difference is denoted f b , and the following relationship holds: After determining the beat frequency from the mixed spectrum, the round trip time τ L for each layer can be computed as

Signal in case of multi-layered snow
In this section, we will discuss the formulation of the received signal for a multi-layered snowpack. A simplified system architecture is shown in Figure 2, for upward-looking bistatic radars. The radar system consists of one transmitter T x and two receivers, R x1 and R x2 . The receivers are placed at distances s 1 and s 2 from the transmitter, respectively. We hypothesize the snowpack is composed of multiple layers, each with a different distance from the receiving radar. Therefore, each receiver will receive multiple echoes of the transmitted signal, due to multiple reflections among layers. To simplify the problem and make it manageable in this first treatment, we will assume that: • changes in the propagation direction that are observed for waves crossing an interface produce negligible geometric effects over the short distances separating the transmitting and receiving antennas; from our previous experience, this assumption is legitimate and sensible.
• that the only significant contributions to the received signal come from multiple bounces on each (single) interface, i.e. we neglect multiple bounces involving different layer interfaces. The echoed signal at each receiver contains more than one beat frequency in its mixed spectrum. For simplicity, if assume a snowpack made up of three layers and each layer with one direct and indirect reflection as well, then we can model the system behaviour as follows:  Layer thickness also unchanged. The small black circles are marked to indicate direct reflections, calculated by later procedure of pattern matching (Section 4).
In Equation 10, the subscript r is used to mark indirect reflection. Subsequently, each receiver R will receive a signal containing all direct X RL i and indirect reflections X RL ir from ith layer: The realistic scenarios foresee multiple paths for each layer, each with its specific delay and path loss. If the difference in delay between the directly and indirectly reflected signals is small, the peaks in the spectrum related to different layers start to overlap, making it impossible to distinguish the correct target. From Equation 8 it can be seen that beat frequencies are proportional to the sweep bandwidth B. In order to overcome the superposition problem, a higher bandwidth is theoretically required Roehr, Gulden, and Vossiek (2007). As already mentioned, a snowpack is generally composed of multiple layers, with each layer interface potentially generating direct and indirect reflections. Because of the different thickness of layers, the signal delay τ d will be different for each layer; the delay or round-trip time from upper layers will be longer than from lower layers. Because of the larger distance of upper layers to the receiver, the directly reflected signal's delay can be sorted as follows: Moreover, this delay is proportional to layer thickness: the thicker the layer, the more the expected delay. A possible solution to the above-mentioned problem is the subject of the remaining part of this paper.

Snow layers modelling
While modelling the system, we have made a few assumptions to make the problem manageable and its treatment concise. For example, although the snowpack could be composed of many layers, only three layers are considered and shown in Figure 2. Moreover, in a realistic scenario, there would be theoretically infinite indirect reflections from the ground and between layers, whereas only one indirect reflection per layer is considered here. Another important point to mention is we assume refraction phenomena are negligible. For sure, changes in refraction indexes between layers will have some impact on the system performance, but for the sake of clarity and problem management, we are ignoring it in this first investigation. Moreover, to make the theoretical discussion of the problem more manageable, we focus on measuring snow layers in the presence of multipath interference, using paired receivers. Other problems and non-ideal features will be tackled in future publications. The geometric propagation distance d contributed by each layer one on the top of other, and having different dielectric properties as given in supplementary material.
The system undergoes different propagation delays and losses and these are proportional to the propagation distance d as in the following equations: In the above equation, the attenuation constant a is a real number. Suppose the first layer has dielectric constant ε 1 , for the second layer ε 2 , third layer ε 3 , and so on. Then, the propagation delay for each layer having varying densities and in consequences dielectric properties ε i can be calculated as where P delay is the propagation delay, and the wavelength λ in the above equation can be written as λ ¼ c 0 ffi ffiffi ε i p f 0 . It can be observed from Equation 12 and Equation 13 that these losses will have a direct impact on propagation distance. These medium losses and propagation delays are calculated for each fraction of d and for all layers. The problem becomes more interesting when we add all the reflected signals for each intermediate layer.
The round-trip time can be calculated as follows from beat frequency of mixed spectrum: Then the snow layer thickness 2 can be calculated based on these round trip times: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi s 2 2 T 2 1 À s 2 1 T 2 2 T 2 2 À T 2 1 r 2 (16) In the above formula, the distance L i is meant as a weighted rather than distance, which incorporates dielectric parameters as weights for each layer. The system model has been built to understand how the signal would be in the real case. This model will be used to develop a strategy to separate the various layers despite multiple reflections.

Snow layering detection
In Section 2, how the features of the received signal are linked to the presence of reflective layer interfaces within the snowpack is described. The location of the different targeted layers can in principle be reconstructed if beat spectra peaks are detected by differentiating them from the system noise. To do so, in our case we used a local-maximum detector based on the zero-derivative method, which searches for peaks in the spectrum. To distinguish the side lobe peaks from the main lobe of the reflected signal, we have defined an amplitude threshold, based on the expected level of sidelobes. The strategy for multipath effect suppression is obviously critical and deserves in-depth analysis. In this section, our simulation results in case of multipath propagation are presented, and their effects in snow-pack stratigraphy are discussed. If the snow has single-layer structure then the effect of multipath reflection is basically harmless, as the greater propagation times of the indirect reflections are sufficient to discriminate them. The problem becomes easier as the distance between receivers increases. The situation becomes more complicated in the presence of multiple layers. Let us consider a realistic scenario, where reflected signals can follow multiple possible paths before being detected. The snowpack will observe multiple reflections. These reflections will include direct and indirect (bouncy) reflection from each layer of the snowpack. The signal from each path will reach the receiver after a specific delay. If the difference in delay between directly reflected and indirectly reflected signal is small enough, they start to superimpose with each other in the frequency spectrum, thus making it impossible to detect the direct reflector. Let us consider an example. From Figures 3-8, the layer densities are 350 kg m −3 , 375 kg m −3 , 387 kg m −3 for layers 1, 2, 3, respectively, based on measurements taken on 20 April 2017, North-Western Alpine environment in Italy. The dielectric constant for each layer can be calculated by using the density ρ i (Hallikainen, Ulaby, and Abdelrazik 1986) and following the equation below: The layer thicknesses are 65 cm, 50 cm, and 35 cm, respectively, for layers 1 to 3. In Figures  3-6, the spectrum with beat frequencies from all three layers is shown with direct reflections and only one indirect reflection from each layer at a time. Then, all direct and indirect reflections together are shown in Figure 6. It can be clearly seen from these figures that the frequency spectrum contains more peaks than the existing layer, three peaks being due to direct reflection from each targeted layer and three more from indirect reflections from each layer. Directly reflected peaks are circled to differentiate them from indirect reflections, calculated by later procedure of pattern matching (Section 4). As we consider reflections from higher layers, the disturbance effect of indirect reflections is smaller; reflections from upper layers accumulate increasing propagation distances and hence the reflected signal becomes weaker. This effect can be seen in Equation 12 and Equation 13, where the propagation distance d drives propagation loss and delay. In order to gain a better insight about where densities change in the snowpack, the signal phase can also be taken into consideration. The spectrum phase diagram is presented in Figures 7 and Figures 8 to show that density layering in snowpack produces visible changes in the phase of the spectral components in the received signal. In other words, a change in snow density is reflected in a change of the local statistics for the phase component of the frequency spectrum. A layer-average line is drawn to highlight the effect. This latter is unfortunately not actionable because phase changes due to layering cannot be clearly discriminated from intra-layer phase changes; still, it gives the sense of the layer information being directly reflected into the spectrum of the back scattered signal.
The effect of multipath interference with the FMCW frequency spectrum can be appreciated in Figures 3-6. Table 2 elaborates the idea further by calculating the thickness of each layer L i measured in metres (m), using the beat frequencies obtained from the frequency spectrum. For this table, it is assumed that the snowpack is composed of three layers; the effect of multipath interference is investigated and shown without applying the proposed pattern matching approach. The thickness estimate for different snowpack layers starting from the bottom layer L1 to top layer L3 is presented in Table 2. Indeed, direct reflection is present on all layers but indirect reflection is added considering one layer at a time to show how the problem changes.    is also worth noticing that, although L est is obtained without considering multipath effects, still it is not a perfect noise-free or loss-free scenario. The L est will contain the impact of propagation delay and medium losses. The purpose of Table 2 is to investigate the combined effect of direct and indirect reflection. This table is also meant to investigate which layers are contributing more error in the estimation process due to their indirect reflections. While investigating the effect of indirect reflection on estimation process, it has also been observed that the top layer computation is more sensitive to interference as it directly overlaps with indirect reflection from the bottom layer. On the contrary, the bottom layer thickness is estimated with a relatively greater accuracy under these circumstances.
In order to see after-effects of multipath indirect reflections in detail, a test was made increasing the number of indirect reflections from two to four bounces. It has been observed that as the number of bounces increase, the propagation distance consequently increases, and thus the effects of interference decrease. Therefore, we can safely assume that more than four bounces do not add any significant degradation to the system performance.
Whenever a complete or partial overlap of peaks in the spectrum occurs, it becomes difficult to determine the actual layer thickness. In Table 2, poor estimates due to this overlap are reported in bold to highlight the error with multipath reflections. In addition to that, as the number of reflections increases, their impact tends to reduce because multiple reflections arrive at the receiver later than the meaningful signal, and with reduced power due to attenuation. For example, with three layers, after five bounces, the effect on the estimation can be considered negligible.
In addition to the above-mentioned problems, some other non-ideal conditions also exist. For example, layers thinner than 15 cm are not visible in principle, as their thickness falls below the range resolution of the radar device. The problem gets even worse if the thin layer is sandwiched between two thicker layers, as its response will deeply mix with the responses from its neighbouring layers.

Pattern matching
So far, the impact of multipath effects on the estimation of snow layering has been analysed in detail in this paper. In this section, instead, a possible new approach is presented to mitigating the multipath interference and to calculating the actual snow density layering in the presence of such interference, as shown in block diagram of Figure  9. There exist different possible approaches to tackle multipath interference. It can be reduced for example by increasing the bandwidth, which in turn increases the resolution, thus making the target peaks more separable. This approach, however, results in increased cost and additional technical requirements placed on the system, which may not be possible to satisfy in all cases. Another possible approach, for a multi-receiver radar architecture as in our case, is to increase the distance between the transmitter and the receivers, as this extends the time and frequency lag between components of the multiple reflections, thus decreasing the chances that reflections interfere with each other. Also, if the distance between receivers increases, then the difference between time of flight of Figure 9. Pattern-similarity-based system block diagram. both receivers also increases, making estimation of snowpack parameters more accurate. However, practical limitations of the hardware and installation will not allow to separate them by an indefinitely long distance, and constraints of a completely different nature will determine the actual location of receivers.
When dealing with multipath scenarios, besides the problem of frequency overlap, it is also difficult to distinguish the actual beat frequencies of layers. The main problem is the mixed spectrum being composed of many peaks from direct and indirect reflections from each layer, as it can be seen in Figure 6. Therefore, the targeted beat frequencies can not be separated with commonly used target detection methods, such as thresholding. Considering all the above-mentioned facts, in order to determine snowcover thickness in the presence of these realistic multipath effects, in this paper a method based on pattern similarity is proposed, distinguishing two possible cases, i.e. the 'simplified case' and 'full case.' These will be defined and explained in the following.
In order to decide which case applies, a comparison is made with a large number of pre-computed peak pattern vectors (PPVs) including most common combinations of layer numbers and thicknesses. These PPVs are the result of simulations carried out in advance on hypothetical, virtual snow layering with different combinations of thickness measures; after each simulation, the resulting PPV is stored with the corresponding layer thickness vector, i.e. the thickness measures of the layers that resulted into the stored peak pattern. The whole set of generated PPVs and their ancillary information is then made available to the thickness estimation algorithm. When a real measurement is carried out, after FMCWbased processing, the prominent peaks of the layers in frequency spectrum and their corresponding beat frequencies are determined. Then, the reference PPVs (PPV ref ) is compared with the PPV obtained in the real case at hand (PPV act ). Their difference is termed dissimilarity index (DSI) and can be calculated as follows for N prominent peaks: The reference case peak pattern with the smallest dissimilarity index is taken as the 'nearest case' and its parameters (e.g. number of layers) are used to specialize further processing in order to reduce the number of unknowns. From this pattern matching process, the number of layers and the first three beat frequencies of the prominent peaks are considered as targeted beat frequencies. The reason for considering the first three peaks is that in FMCW radar direct reflection from a targeted layer always comes earlier than indirect reflections from the same layer, as direct reflections have a smaller distance to cover. In the next step, based on dissimilarity index, a decision will take place to move control to 'full case' or 'simplified case.' The aforementioned decision process is elaborated in detail in Figure 10.

Full case
As visible in Figure 10, if the dissimilarity index is greater than 10%, then the 'full case' is activated. In the 'full case,' for sake of more accurate estimation, the beat frequencies of the actual signal are used to determine the time of flight T 1 for receiver 1 and T 2 receiver 2 for each layer, using Equation 14 and Equation 15. Then, these times of flight are used to calculate the actual layer thickness as in Equation 16.

Simplified case
On the other hand, if the dissimilarity index is less than 10%, then the 'simplified case' is activated. In this case, where patterns match quite well, the previously calculated and stored layer thicknesses of nearest case are directly assigned and considered as the target layer thicknesses, thus reducing the on-site processing load. In the 'simplified case,' the difference with respect to the 'full case' is found only after the decision process, whereas the rest of the procedure will remain the same. Therefore, snow layering estimation block in Figure 9 is superseded in the simplified case.

Overall operation
Both cases are operationally used to determine the number and characteristics of layers based on the direct and indirect reflections pattern. Where usable, the 'simplified case' is faster and more efficient, as in this case the previously stored snow thickness is directly assigned at the expense of precision. An optimal trade-off may be sought by tuning the threshold value for deciding whether the full or the simplified case applies.
An extensive amount of simulations have been carried out to verify this method, but for sake of conciseness, in this paper, only four representative examples are presented. In Figure  11, each example fits into one row, and three database reference layer thickness (DB-RLT) cases are presented on the three columns. In this figure, the pattern matching step of Figure  9 is simulated precisely. The comparison between PPV of the actual signal layer (AL) and the reference case layer is shown for a snowpack with double layer structure; both snow layers are assumed to generate direct and indirect reflections. The difference in PPVs of three reference test cases with the corresponding actual case PPV is calculated and shown in  Figure 11. In all of the four represented examples, the test case 1 was selected as the one with the least difference with actual case and is thus considered as the 'nearest case. ' In examples 1 and 2 of Figure 11, the dissimilarity index is smaller than 10% and thus the 'simplified case' is activated, whereas in the remaining two examples the dissimilarity index is greater than 10%; thus, the 'full case' will be activated, as shown in flow graph Figure 10.
There is one more interesting point about Figure 11. Consider the absolute difference of the 'nearest case' database entry (leftmost column) to the actual layer thickness reported on the rightmost script on each row: |DB-RLT-ALT|. Such difference is the same in the bottom two examples as in the top two examples (i.e. 0.1 (m)). Still, the examples in the two bottom lines fall into the 'simplified case,' whereas those in the two top lines fall into the 'full case.' The reason as described earlier and shown in Table 2 is connected to the fact that the top layer is subjected to greater error because of indirect reflection, hence it is weighted less on the difference computation.

More simulations on real data
In Figures 12 and Figures 13, real Figure 12, the results of applying the pattern similarity method to the snow pack with only one layer are presented. The reported PPVs' difference trend reveals that as the selected hypothesized thickness vector, or DataBase Reference Layer Thickness (DB-RLT) approaches the Actual Layer Thickness (ALT), the dissimilarity index decreases, as expected.
Similarly, in the upper example of Figure 13, the presence of the same layer ALT [0.45 0.6] with the same density can be observed in the database. The aforesaid is considered as a 'simplified case' since the DS Index will be approximately zero.
In the bottom example, for ALT [0.65 0.5], there was no exact match within the available layer thickness combination in the database of the presented example; hence, the method had to look for the nearest pattern combination. The figure shows that PPVs' difference for the first four DB-RLT is minimal, with DS Index [8.9826, 13.4360, 9.9649, 10.1612]; hence, we will consider the first database entry as the nearest one and assign its thickness directly to the actual signal as its DS Index is smaller than 10%. This is also illustrated in the flow graph of Figure 10: when the DS index value is smaller than 10%, the corresponding DB-RLT will be directly assigned, and it will be considered as the targeted thickness.

Building the reference database
As mentioned previously, a database has been formed for all possible layer thickness combinations between a minimum of 20 cm to a maximum of 200 cm, with a step size of 10 cm or 5 cm. From an implementation standpoint, the database will consist of a matrix whose rows will be associated with the layers and the columns with size steps. For example, for a single-layer structure, this means 19 possible columns of snow layer thickness with a step size of 10 cm and 37 columns with a step size of 5 cm. Similarly, for double-layer structures, we will have 361 combinations for a step size of 10 cm and 1369 combinations for a step size of 5 cm. Additional layers will produce even more combinations, with a total number growing quickly with the number of layers. In order to build a generic database for all layers, we make two simplifying assumptions. First, we assume we have a maximum of k layers. Second, we assume the received mixed spectrum will contain 2k prominent peaks, one from direct and one from indirect reflection from each of the k layers. Hence, the total number of columns will be 2k and the number of rows will depend on the numbers of thickness combinations we are considering according to step size and number of layers. As expected, the use of finer step sizes in generating the database will translate into more precise layer thickness estimates, but at the expense of larger database size. A trade-off is unavoidable between accuracy and storage requirements; for example, with an increase in step size from 5 cm to 10 cm, on average the error rate also increased by about 48%. Figure 16 shows instead the normalized error rate as a function of the number of layers and the step size. The normalized error rate is calculated by dividing the total error by the number of layers. It can be noted that for increasing stratification (i.e. number of snow pack layers) the normalized error rate increases. Moreover, as expected, the full case provides more accurate estimation in general, but the gap tends to decrease as the number of layer increases. The analysis included also changing the step size from 10 cm to 5 cm. The performance is always better with the 5 cm step but obviously at the cost of a heavier computational burden.

Error analysis on density change
In order to check the robustness and validity of proposed method with changing densities, the error rate has been analysed for both layer 1 and layer 2 in a doublelayer snowpack in Figures 14 and Figures 15, respectively. The purpose of this analysis is to determine how the system performance changes if the combination of DB reference layer thickness and actual signal layer thickness remain the same, but the densities are changed. Simulations were performed with receivers located at distances RX 1 = 85 cm and RX 2 = 200 cm from the transmitter. For each density value, the thickness was estimated, and the error computed after the estimation. In order to reduce other sources of error, in our experiments, we defined a 5-cm step for building the reference database. The layer thickness values are [0.63 m 0.75 m] for both figures. Figure 14 reports the results of varying density of layer 1 while keeping layer 2 at 300 kg m −3 . Figure 15 is similar except now it refers to varying density of layer 2, with layer 1 density fixed at 300 kg m −3 . It can be observed that, generally, the estimation error increases with snow density. Moreover, changes in the density of layer 2 (top layer) are more likely to impact on error.

Analysis of error and success probability with increasing number of layers
An error analysis has been carried out to assess how the presented system performs for increasing numbers of layers. In general, as visible in Figure 16, more layers in snowpack result in an increase in the normalized error as expected. For a given number of layers, however, finer sampling in the database will lead to smaller errors, as a close match in the database becomes easier to find. Again in Figure 16 one can notice that 5 cm sampling leads to smaller error rates compared with the corresponding 10 cm case. Surely, the system may generate an error during the estimation of the thickness of layers, but the probability of accurately detecting the number of layers is still very high as it can be seen in Figure 17. The foremost limitation in detecting the number of layers is linked to the radar range resolution; layers thinner than 15 cm will generally be invisible. It has been observed that as the number of layers in snowpack increase, the success probability of detection decrease. In any case, above the threshold, thicker snowpack are more likely to be correctly detected, having more possibility of wider layers.

Conclusion
This paper presents an investigation on the effects of multipath interference in radarbased detection of layers in a snowpack using an FMCW bi-static radar architecture. In this context, multipath interference is due to indirect reflections linked to snowpack density layering, with consequent reflection at layer interfaces. For accurate detection of snowpack stratification in such a challenging multipath environment,   a system model has been built and simulation results have been analysed to achieve deeper insight about snow density layering in the frequency spectrum. Numerical and graphical results confirm that multipath reflections from different snow layers deteriorate the snowpack estimation process. A novel algorithm, based on the similarity between actual and sample patterns of spectral power peaks, is presented to distinguish snow density stratification in the presence of multipath interference. This algorithm foresees two different operational modes, i.e. full mode and simplified mode, which helps to reduce the on-site processing load and makes the system faster and more efficient. Our test results appear to support the point that the method can effectively detect layer patterns and help setting an effective starting point for finer and more precise computation of layer thickness and parameters using more refined techniques. Our future research will aim at improving the pattern matching method to accommodate more noisy, real-world signals, and at developing versions of the algorithm suitable to operate in an energy-limited computational environment of the in-situ sensors.

Notes
1. Please see supplementary material for proof of Equation (7). 2. For detailed derivation of thickness L i formula, please see supplementary material.