Evaluation of filtering methods for use on high frequency measurements of 1 landslide displacements 2

Sohrab Sharifi, Michael T. Hendry, Renato Macciotta, Trevor Evans 3 Department of Civil and Environmental Engineering, University of Alberta, Edmonton, AB, 4 Canada 5 Canadian National Railway, Kamloops, British Columbia, Canada 6 7 Abstract 8 Displacement monitoring is a critical control for risks associated with potentially sudden slope 9 failures. Instrument measurements are, however, obscured by the presence of scatter. Data 10 filtering methods aim to reduce the scatter and therefore enhance the performance of early 11 warning systems (EWSs). The effectiveness of EWSs depends on the lag time between the onset 12 of acceleration and its detection by the monitoring system, such that a timely warning is issued 13 for implementation of consequence mitigation strategies. This paper evaluates the performance 14 of three filtering methods (simple moving average, Gaussian-weighted moving average, and 15 Savitzky-Golay), and considers their comparative advantages and disadvantages. The evaluation 16 utilized six levels of randomly generated scatter on synthetic data as well as high-frequency global 17 navigation satellite system (GNSS) displacement measurements at the Ten-mile landslide in 18 British Columbia, Canada. The simple moving average method exhibited significant 19 disadvantages compared to the Gaussian-weighted moving average and Savitzky-Golay 20 approaches. A framework is presented that can be followed to evaluate the adequacy of different 21 algorithms for minimizing monitoring data scatter. 22


Introduction 26
Landslides are associated with significant losses in terms of mortality and financial consequences 27 in countries all over the world. In Canada, landslides have cost Canadians approximately $10 28 billion since 1841 (Guthrie, 2013) and more than $200 million annually (Clague and Bobrowsky, 29 2010). Essential infrastructure, such as railways and roads that play vital roles in the Canadian 30 economy, can be exposed to damage as they transverse landslide-prone areas. Attempting to 31 completely prevent landslides is typically not feasible, as stabilizing options and realignment may 32 not be cost-effective nor environmentally friendly. This accentuates the significance of adopting 33 strategies that require constant monitoring to mitigate the consequences of sudden landslide 34 collapses (Vaziri et al., 2010;Macciotta and Hendry, 2021).  Table 1 and shown in Fig. 1, respectively. Scenarios considered decreasing 104 trends of y from a value of 1 to 0, reflecting cumulative negative displacements or inverse-105 velocities; however, it was acknowledged that absolute cumulative displacements and absolute 106 velocities could show increasing trends. In this regard, the evaluation of synthetic data focused 107 on timely identification of changes in trends as those associated with accelerating and 108 decelerating periods, and the results are valid if the scenarios are mirrored to vary from 0 to 1. 109 Nine of the scenarios are referred to as harmonic scenarios, which are characterized by gradual 110 changes in the trend of parameter y. The remaining three scenarios show sudden variations at or 111 near x=0.5, and are referred to as instantaneous scenarios. Considering the discrete nature of 112 instrument measurements, and to account for different ranges in measurement frequencies, each 113 scenario was generated several times, each time with a different number of points ( The next step was adding random scatter to the scenarios to represent unfiltered displacement 123 measurements. Macciotta et al. (2016) show the scatter in displacement monitoring for a GNSS 124 system used in their analyses fitted a Gaussian distribution. This was also validated for the data 125 scatter for the case study in this paper and is presented in a subsequent section. Based on this 126 observation, the scatter was randomly produced from a normal distribution centred at zero, with 127 extreme values truncated between −1 and 1 and a standard deviation of 0.20. Random generation 128 of the scatter followed the techniques outlined in Clifford (1994) known as acceptance-rejection 129 method, which generates scatter values through a series of iterations until the initial normal 130 distribution is generated. The amplitude of the scatter around the trend in parameter y was defined 131 for each scenario based on scaling the randomly generated scatter. This allowed investigation of 132 the effect of different scatter magnitudes on the performance of the filters. Scaling was done by 133 defining the ratio n/t, which is the ratio of scatter amplitude (maximum deviation around the trend, 134 termed n) to the range of values of the trend (t) in each scenario. Six levels of n/t (0.001, 0.005, 135 0.010, 0.050, 0.100, and 0.150) were considered when performing the analysis to cover a range 136 of possible levels of scatter in unfiltered measurements. Fig. 2 where y î is the filtered value, y j is the unfiltered value, and p is the window length. The window 151 length is constant across the dataset except for the regions near the boundaries as fewer points 152 are available. Accordingly, p will be adjusted to the number of available points that are indeed 153 less than the value set by the user. This will cause variation in the effectiveness of the method at 154 the extremes, which need to be considered when evaluating the results of this approach. 155

Gaussian-weighted moving average 156
Varying the weights of the measurements within the calculation window in SMA can be used to 157 develop different filtering methods. The highest weight can be given to the measurement at the 158 time for which the calculation is being done, with weights decreasing for measurements farther 159 away in time. One simple weighting function that can be adopted is the Gaussian (normal) 160 distribution. The filter that assigns weights based on a Gaussian distribution for the averaging 161 process is: 162 where wj is the weight coefficient based on the Gaussian distribution and the other terms follow 164 the same definition as per SMA. 165

Savitkzy-Golay 166
S-G fits a low-degree polynomial equation to the unfiltered measurements within a window and 167 defines the filtered measurements using the fitted curve (Schafer, 2011). Although this procedure 168 seems dissimilar from the weighted averaging discussed above, it can be transformed into a 169 kernel concept using the least-squares method if the data points are evenly spaced. The detailed 170 procedure is presented in Appendix A. Fig. 3 shows the weight kernel over a window of seven 171 points attained by fitting a quadratic polynomial. An immediate observation is that some points 172 are given negative weights.  This is to reflect the reality of displacement monitoring information as applied for EWSs. This was 183 achieved by applying the filters using the time for which the calculation is being made as the 184 central value, but only utilizing the first half of kernels to assign the weights (the weights are 185 multiplied by 2 in comparison to a symmetric window to keep the sum of weights equal to 1). 186 All of these filters require the definition of a bandwidth. A roughness factor was defined to aid in 187 the evaluation of the effect of bandwidth in reducing scatter. This factor is defined as: 188 where J 2 is the roughness factor, ŷ '' is the second derivative of filtered measurements, R a is the The filters are not expected to remove all scatter, and the error attributed to the residual scatter 204 can be calculated using the root mean square error (RMSE). Given that velocity values are usually 205 used as thresholds in an EWS, one concern is whether the filter should be applied to displacement 206 values or to velocity values derived from unfiltered displacements. To address this issue, two 207 different approaches to filtering were investigated: direct and indirect. As a result, two different 208 approaches using the RMSE were also utilized here. 209

Direct scatter filtration 210
Direct filtration means the filter is applied to the diagram of interest. If the filtered displacement 211 values are the goal, and the filter is applied to unfiltered displacement values, then the filtering 212 process is called direct filtration. The same concept applies when velocity values are derived 213 using unfiltered displacements and the filters are then directly applied to the velocity values. In 214 this approach, the RMSE follows Eq. 5: where RMSEd is the measurement of error in direct filtration, y i is the value of the true trend (for 217 the synthetic scenario), y î is the filtered value, and m is the total number of points. This approach

Indirect scatter filtration 221
Some EWSs can apply the filter to the displacements but use velocity trends as the metric for 222 evaluation. In this case, the filtered velocity values will be computed using the filtered 223 displacements. Indirect filtration indicates the diagram of interest is the first derivative of the 224 diagram to which the filter is applied. The RMSE in this case is defined as: 225 where RMSEi is the measurement of error in indirect filtration, y i ' is the first derivative of the true 227 trend, ŷ i ' is the first derivative of filtered data (derived velocity after the filter is applied to the 228 displacements), and m is the total number of points. 229

Lag Quantification 230
Only antecedent measurements are fed into the filters, which is expected to result in a lag between 231 the true trend and when these are identified by the filters. This lag means the calculated value of 232 velocity or displacement occurred sometime in the past. Consequently, reducing this lag means 233 less time is lost with respect to providing an early warning. To quantify the induced lag, the filtered 234 diagrams of all scenarios at all n/t ratios and BR values were shifted backwards a number of 235 points equivalent to 0.001 (0.1 %) to 0.1 (10 %) of all generated points. This is referred to as the 236 shift ratio (SR). This shift of filtered diagrams is expected to increase their similarity with the true trend until the best correlation is achieved. The R 2 test was used to determine how well the shifted 238 and filtered results replicate the underlying trend. 239

Geocubes Differential GNSS System 240
A Geocubes system is a network of differential GNSS units that works with a single frequency 241 2018). The ability of this system to achieve real-time positioning, remote data collection, and 246 processing makes it a suitable candidate for incorporation into an EWS. As a result, Geocube 247 data are used in this study to evaluate the performance of the three mentioned filters. 248

Outlier Detection 249
Outlier detection techniques have been proposed based on the statistical characteristics of 250 datasets. One common example is the Z-score method, which calculates the mean and standard 251 deviation of data within a defined interval and identifies outlier data as those beyond three 252 standard deviations from the mean (Rousseeuw and Hubert, 2011). A limitation of this kind of 253 approach is the sensitivity of the mean and standard deviation to the outlier data points, which 254 has led to the development of other methods that use other indices such as the median (Salgado 255 et al., 2016). One such technique that was adopted in this study is the Hampel filter (Hampel, 256 1971). In this method, the median of the displacement measurements within a running bandwidth 257 is calculated and data outside a defined threshold from the median are identified as outliers. The 258 threshold is defined as a constant (threshold factor) multiplied by the median absolute deviation. 259 An asymmetric window with a bandwidth ratio of 0.004 (0.4%) and a threshold factor of three were 260

Study Site -Ten-mile Landslide 264
The Ten-mile landslide is located in southwestern British Columbia (BC), in the Fraser River 265 Valley north of Lillooet (Fig. 4a). It is a reactivated portion of a post-glacial earthflow (Bovis 1985)    shows that, for the NASD, the error depends linearly on the BR for all of the filters and does not 314 depend on the scenario or n/t ratio. SMA shows the greatest difference from the true trend, 315 followed by GWMA (approximately 60% less difference than SMA). S-G, on the other hand, 316 almost lies on the horizontal axis for all of the BRs, which means the filtered results yield near 317 zero error. Fig. 7 also shows how error increases as BR increases. This can be attributed to the 318 fact that an asymmetric window was utilized, which leads to a lagged response of the filter. As 319 more points are included in the filtering procedure (increasing BR), this lag increases and, 320 consequently, causes higher error. The RMSEd of filters for the instantaneous synthetic scenarios 321 are shown in Fig. 8. In Scenario 10, the same behaviour as for the harmonic scenarios can be 322 seen from SMA and GWMA, whereas S-G is not as accurate. This is more noticeable in Scenarios 323 11 and 12 in which S-G becomes less accurate than GWMA at high BRs. This result shows that 324 S-G cannot handle the instantaneous scenarios as satisfactorily as it does the harmonic ones. 325 The errors related to SMA and GWMA for the instantaneous synthetic scenarios show non-linear 326 behavior, and are greater when compared to the harmonic scenarios. The results show the error considerably reduced as the BR increases to 0.01 for SMA and GWMA 335 and 0.02 for S-G, and has an asymptotic tendency above these BR values. S-G has the highest 336 error at low BR values in comparison to SMA and GWMA, but shows the least error at BRs above 337 0.01. At BR values over 0.03, fluctuations do not vary significantly with BR (Fig. 6). In this range 338 of BR values, the error of GWMA is either equal to or slightly less than the error of SMA, and S-339 G shows the least error. The RMSEi results for the instantaneous scenarios (Fig. 10) are similar 340 to those for the harmonic scenarios for high n/t ratios (0.05, 0.10 and 0.15). For low n/t ratios, the 341 GWMA is superior at BRs above 0.06, and S-G has the worst performance.  Scenarios 11 and 12 were further analyzed to evaluate how the filter performance is affected by 347 the presence of sudden peak(s). Fig. 11a shows the true trend of Scenario 11 along with two bandwidth increases, the peak in measurements is identified at a later time than the true trend (x 350 = 0.5) and the magnitude of the peak is reduced (more than 70% reduction at BR=0.10). 351 Furthermore, as BR increases, the "instantaneous" nature of the peak is lost to a more transitional 352 variation. This highlights the disadvantage of SMA when handling sudden changes in 353 displacement trends. The calculated x value of the peak in Scenario 11 is plotted for different BR 354 and for all three filters in Fig. 11b. This figure shows the time at which the peak is identified lags 355 as the BR increases for all filters; however, GWMA and S-G identify the peak within a much 356 smaller lag, independent of the n/t ratio. As an example, for a year of monitoring data at a 357 frequency of 30 s and BR=0.10, SMA, GWMA, and S-G predict the peak point approximately 17, 358 3.5, and 2.7 days after the real peak, respectively. Fig. 11c shows the variation of the peak 359 magnitude with respect to BR for all three filters. Both SMA and GWMA underestimate the peak 360 value, and the difference between the calculated peak and real peak increases as BR increases. 361 SMA calculations underestimate the peak more than twice as much as GWMA. On the contrary, 362 S-G intensifies the peak up to BR=0.04, with the impact tending to diminish for higher BR values; 363 it predicts the true value at a BR value of almost 0.09. 364 Scenario 12 was used for a detailed evaluation of the performance of these filters to conserve the 368 underlying original trend. Fig. 12 shows Scenario 12 and the filtered results for all three filters and 369 for an n/t ratio of 0.15. This scenario and parameters were selected for illustration purposes as However, the peak is artificially intensified at BR=0.04, and a significant drop occurs well beyond 379 the true trend immediately after the second peak for both BR values (pulsating effect), which was 380 also observed in Scenario 11. Increasing the degree of the polynomial fitted as part of the S-G 381 methodology was not effective at eliminating this effect. The pulsating effect was also observed 382 when a symmetrical window was utilized and is attributed to the negative weights in the S-G 383 kernel. 384 quantified by the SR, is higher when the trend change is more pronounced; therefore, the 414 correlation between SR and BR is different for different scenarios. Fig. 14b shows the mean 415 correlation between the SR and BR, for all scenarios and n/t values, bounded by one standard 416 deviation, for GWMA and SMA. Table 3 shows linear and quadratic regressions of this correlation 417 and the strength of the correlation in terms of R 2 and RMSE. Fig. 14b shows quantitatively that 418 GWMA lags less than SMA with respect to identifying changes in measurement trends. Moreover, 419 the uncertainty associated with lag in SMA is greater than in GWMA because of larger standard 420 deviation. Fig. 14b quantifies how increasing BR values increases the lag with respect to 421 identifying true measurement trends, and although high BR values decrease the scatter in data, 422 the BR should carefully balance minimizing both scatter (J 2 ) and lag (SR). S-G is not included in 423 this analysis as the method provided no significant lag in identifying changes in measurement 424 trends; however, it had the disadvantages previously noted including pulsating effects and 425 overestimating peak values. 426 Table 3 Regression correlations between shift ratio (SR) and bandwidth ratio (BR) with the strength of the  impacts the performance of S-G, as the presence of the outlier generates a peak that follows the 458 outlier measurement and is followed by a sudden decrease that goes well beyond the data trend. 459 SMA tends to widen the range affected by the outlier more than GWMA but, for most part, the 460 filtered results are almost parallel to the underlying trend. All filters appear to be significantly 461 impacted by the outlier value, suggesting a pre-processing filter is required to remove outliers 462 regardless of the use of SMA, GWMA, or S-G to reduce scatter. The outliers were successfully 463 identified and removed after application of the Hampel algorithm, and the above-mentioned 464 effects were no longer observed in the filtered results. The lag between unfiltered and filtered data for Geocube 46 (Fig. 15) is consistent with NASD 471 results. The NASD lag quantification results ( Fig. 14b and Table 3) were used to provide a 472 correction value for the filtered Geocube results. To determine whether the results of lag 473 correction using the mean correlations derived from NASD (Table 3)  inferred from the filtered shifted displacement measurements closely agrees with that considered 489 to be true scatter, and therefore the filtered displacement measurements are corrected for lag. 490 This suggests the correlations in Fig. 14b and Table 3 based on NASD are applicable to minimize 491 the lag for the Geocube system at the Ten-mile landslide. 492

Conclusion 499
This study evaluated the suitability of SMA, GWMA, and S-G filters for scatter reduction of 500 datasets targeted for use in an EWS. A total of different 12 scenarios with harmonic and 501 instantaneous changes were synthetically generated and random variations with Gaussian 502 distribution then added to produce unfiltered results. The three filters considered were then each 503 applied with different bandwidths and the error computed. These filters were also successfully 504 applied to the records from two Geocubes installed on the Ten-mile landslide. The results led to 505 the following conclusions: 506 • When used for direct filtration of harmonic scenarios, the error resulting from the GWMA 507 approach was approximately one-third that of the SMA approach. The S-G approach 508 resulted in near zero error regardless of the BR and n/t. When used for direct filtration of 509 instantaneous scenarios, the superiority of S-G is no longer unconditional and depends 510 on the BR; this reflects the fact that S-G cannot appropriately handle peaks in the velocity 511 diagram. 512 • When used for indirect filtration of harmonic scenarios, S-G again outperforms the other 513 methods. The error associated with GWMA is marginally less than for SMA. These 514 observations are not valid when the filters are applied to instantaneous scenarios, as 515 GWMA results in less errors than S-G at BRs above 0.03. 516 • Detailed investigations with Scenarios 11 and 12 demonstrated that SMA distorts the 517 underlying trend by displacing and sometimes neglecting peak(s), while GWMA and S-G 518 tend to preserve them somewhat similarly. 519 • Due to the presence of negative weights in the S-G kernel, some artificial smaller troughs 520 and peaks are created after major peaks. This phenomenon, referred to as pulsating effect 521 here, results in unfavorable performance of S-G on the velocity and displacement 522 diagrams, especially in the presence of outliers.
• Investigations on the roughness factor reveal the BR should be at least 0.04. Taking this 524 into account, GWMA seems to be the most reasonable option as the related uncertainties 525 are much lower than for S-G and the error is acceptably less than for SMA. 526 • A consequence of using asymmetric windows in the filtering process is a lag in the SMA 527 and GWMA results that increases with increasing BR. Lag quantification suggested a 528 correlation between the needed shift and BR that can be used to eliminate the lag. SMA 529 requires approximately three times the shift of GWMA on average. 530 • Application of these filters to displacement data reported by Geocubes illustrates that SMA 531 and S-G are unable to properly handle data points at the beginning of the dataset (i.e., 532 near the boundary) in indirect filtration of the velocity diagram. Moreover, SMA and S-G 533 are inclined to respectively understate and overstate peaks and fluctuations in the velocity 534 diagram. Overall, GWMA provides the most reliable filtered values for velocity with no 535 distinct difference between direct and indirect filtration. 536

Appendix A 537
Consider a polynomial of degree k that is intended to be fitted over an odd number of points 538 denoted as z. The weighting coefficients of the Savitzky-Golay filter can be extracted from the first 539 where T operator is the transpose of a matrix and J is the Vandermonde matrix, with elements at 542 the ith row and jth column (1≤i≤z and 1≤j≤k+1) that can be achieved as follows: 543 J ij =m i j-1 , (8) 544 where m is the local index of points (-(z+1) 2 ⁄ ≤m≤ (z+1) 2 ⁄ ). As an example, the kernel of an S-G 545 filter that fits a quadratic polynomial (k=2) over seven points (z=7) is attained here. In the first 546 step, J is set up as follows: 547 https://doi.org/10.5194/nhess-2021-212 Preprint. Discussion started: 24 August 2021 c Author(s) 2021. CC BY 4.0 License.
(10) 550 The second and third rows of C are the coefficients to find the filtered values' first and second 551 derivations at the point of interest, respectively. 552

Data availability 553
The synthetic database can be generated through the comprehensive steps provided here. The 554 Geocube measurements of Ten-mile landslide displacement are not to be publicly available. 555

Competing interests 561
Ground Hazard Research Program, which is funded by the Natural Sciences and Engineering 566 Research Council of Canada (NSERC ALLRP 549684-19), Canadian Pacific Railway, CN, and 567 Transport Canada. 568