key: cord-0668740-ouf9voib
authors: Ma, Yiming; Anastasiou, Andreas; Wang, Ting; Montiel, Fabien
title: Detecting change-points in noisy GPS time series with continuous piecewise structures
date: 2022-02-24
journal: nan
DOI: nan
sha: bad7c56c951c4ce6c9f4ef469ec240b7b1daf881
doc_id: 668740
cord_uid: ouf9voib

Detecting change-points in noisy data sequences with an underlying continuous piecewise structure is a challenging problem, especially when prior knowledge of the exact nature of the structural changes is unknown. One important application is the automatic detection of slow slip events (SSEs), a type of slow earthquakes, in GPS measurements of ground deformation. We propose a new method based on Singular Spectrum Analysis to obscure the deviation from the piecewise-linear structure, allowing us to apply Isolate-Detect to detect change-points in SSE data with piecewise-non-linear structures. We demonstrate its effectiveness in both simulated and real SSE data.

Slow slip events (SSEs) are fault slips occurring at the interface between tectonic plates (the subduction zone), generally recorded by the Global Positioning System (GPS) networks (Mitsui and Hirahara, 2006; Ito et al., 2007) . The aim of identifying SSEs is to arXiv:2202.12414v1 [stat.AP] 24 Feb 2022 monitor the change in their ongoing recurrence pattern, which could potentially provide a feasible way to forecast impending large earthquakes (Obara and Kato, 2016) . In the past decades, over 15, 000 GPS stations have been deployed in many subduction zones worldwide (Bedford and Bevis, 2018) . The drastic proliferation of GPS data makes it impossible to inspect SSEs visually, this being extremely time-consuming and introducing subjective human errors. An automated detection tool for SSEs is thus crucial.

In this study, we refer to the GPS data recording SSEs as SSE data. Fig. 1 shows observed SSE data in the Hikurangi subduction zone, New Zealand. In periods when no SSEs occur, the overall trend of the GPS data is a noisy linear process with some seasonal effects from known sources. However, the occurrence of an SSE redirects the original trend to a different direction. Upon completion of the SSE, the trend reverses back to its previous state. The start and the end times of an SSE are change-points in the SSE data, so we can regard the detection of SSEs in GPS data as a problem of offline change-point detection (CPD) , where, in general, the idea is to test whether pre-recorded data are homogeneous or not in a statistical sense. (Wallace, 2020) , which are discussed in §5.3.

Two important elements of a parametric CPD method are the model characterizing the underlying signal (e.g. piecewise-constant, piecewise-linear, or other) and the associated contrast function (Cho and Kirch, 2020; Yu, 2020; Truong et al., 2020) , which is derived based on the assumed underlying structure and it is employed in order to decide whether there are any change-points. Existing CPD methods are mainly designed for signals with known structures, e.g. piecewise-constant or piecewise-linear (Truong et al., 2020) . However, the exact form of the underlying signal of many real-world processes is often unknown and likely non-linear. The SSE data are one of the representative examples where the structure of the underlying signal is not known . If the assumed model is not consistent with the underlying data structure, many existing detection methods fail to correctly segment the data sequence (Truong et al., 2020) .

As shown in Fig. 1 , SSE data appear to have a piecewise structure. However, the exact form of the underlying signal is unknown, making it impossible to construct a reliable contrast function to detect the change-points. Most geophysical models for SSEs suggest that the signal is piecewise non-linear, although the nature of the non-linearity is unknown. Therefore, the problem of detecting change-points in SSE data can be treated as CPD for a certain class of signals, which have a piecewise structure with an unknown form, but may be approximated as piecewise-linear. We refer to this kind of signals as piecewise-non-linear signals.

In this paper, we propose a new detection method for signals with unknown continuous piecewise structures, called singular spectrum analysis isolate-detect (SSAID), which does not require knowledge of the exact structure of the underlying signal in the noisy data. Instead of developing a specific model for the underlying signal, SSAID aims at obscuring the differences between the piecewise-linear model and the real model by first decomposing the signal into spectral components, then adding noise to these components, and finally applying existing CPD methods for piecewise-linear signals. To be more specific, we use singular spectrum analysis (SSA; Ghil et al., 2002) for decomposition and Isolate-Detect (ID; for CPD. We conduct a range of simulations to evaluate the performance of the new method and compare it with existing methods.

In §2, based on simulated SSE data, we show how the detection methods commonly used for continuous piecewise-linear signals fail to detect change-points in continuous signals with unknown piecewise structures. In §3, the new SSAID method along with some model assumptions is introduced. In §4, we present results of applying SSAID to simulated SSE data, and conduct sensitivity tests to explore the influence of parameters on the detection performance of SSAID. In §5, we demonstrate its detection capability with the real SSE data observed in the Cascadia subduction zone. A sliding window approach is also introduced into SSAID to improve its adaptability to real SSE signals with more complicated structures. Discussions and conclusions are included in §6.

In this section, we aim to quantify performance of existing CPD methods designed for piecewise-linear signals at detecting change-points in a noisy signal with an unknown piecewise structure. As an example, we use simulated SSE data obtained from a geophysical process-based model. We evaluate the accuracy of these methods with respect to both the estimated number of change-points and the estimated change-point locations. The new method introduced in §3 uses the findings of this section to extend the applicability of existing CPD methods to signals with different piecewise structures.

To simulate GPS data that contain SSEs, we adopt a deterministic geophysical model of slip at the interface between two tectonic plates, i.e. the subduction zone. This approach is commonly used to simulate the evolution of SSEs in the field of geoscience, even though the mechanism of generating SSEs still remains unclear . The evolution of the SSEs is simulated by a system of coupled non-linear ordinary differential equations (ODEs), which is solved numerically. Therefore, the exact form of the simulated SSE data is not known. More details on the ODE system and numerical solution methods can be found in and the supplement. Fig. 2 (a) shows a simulated signal with 5 SSEs, in a one-year period. The recurring periodic pattern is consistent with direct SSE observations from GPS data. By changing the model parameters, we simulated other SSE signals with different numbers of SSEs per year (see Fig. S2 in the supplement).

In these simulated signals, we define the start of an SSE when the slip velocity becomes 20% higher than the plate velocity, and the end of an SSE when the slip velocity becomes lower than 1.2 times the plate velocity. The plate velocity refers to the slip velocity of the subducting plate in our model (for details see Fig. S1 in the supplementary file).

We construct the noisy simulated data X t using

where T is the length of the data sequence, and f t is the simulated pure SSE data (generated by the deterministic geophysical model), standardised through the Z-score normalisation for ease of comparison. The second term C wn × t in Eq. (1) denotes the noise model contained in X t . We assume that t are independent, Gaussian random variables with mean zero and variance one. The noise level C wn is the standard deviation of the noise model. Fig. 2 (b) and (c) show two examples of simulated noisy SSE data with different noise levels. 

Both the observational and simulated noisy SSE data, shown in Fig. 1, Fig. 2 

(c), appear to have a piecewise-linear structure, even though the pure signal ( Fig. 2 (a) ) does not. Therefore, the existing CPD methods for continuous piecewise-linear signals might be useful in detecting SSEs in GPS data. We now test the performance of three well-established CPD algorithms for piecewise-linear signals on our simulated data: the Narrowest-Over-Threshold (NOT) algorithm (Baranowski et al., 2019) , the Continuouspiecewise-linear Pruned Optimal Partitioning (CPOP) algorithm (Fearnhead et al., 2019) and the Isolate-Detect (ID) method . We focus on the simulated signal with 5 SSEs (i.e. 10 change-points; see Fig. 2 ), which have a duration of approximately one week each. We carry out 10, 000 simulations for each noise level, with noise levels C wn = 1%, 2%, · · · , 100%. Fig. 3 (a)-(c) show that the performance of NOT is not satisfactory for the data as NOT overestimates or underestimates the number of true change-points depending on the noise level C wn , while CPOP and ID can consistently detect the number of true change-points for some noise levels in a certain range. We compare the three methods with respect to the percentage of successful detections (R sd ) for each noise level, which is calculated by

where ξ is the number of simulations for each noise level (i.e. ξ = 10, 000 here), and α is the number of successful detections. Here, we define a successful detection when two conditions are met: (1) the number of the estimated change-points,N , is exactly the correct number of change-points, N , (i.e. 10 here), and (2) the root mean squared error (RMSE) of the detected change-points in those cases withN = N is less than a predefined threshold value v. The RMSE between the estimated change-point locations, p i , and the true locations, q i , is calculated by RMSE = N i=1 (p i − q i ) 2 /N , only for simulations that satisfyN = N . In general, the larger the threshold values v, the higher the successful percentage R sd , but we need to control the threshold to be small enough so that the detection error is acceptable. We show the effect of different threshold values v on R sd in the supplement (see Fig. S3 (a)). In our simulated noiseless signal, the duration of each SSE is about 7 days, and the corresponding recurrence time is 74 days.

An RMSE of 3 days (i.e. v=3) is an acceptable error for detecting such SSEs (Holtkamp and Brudzinski, 2010; Nishimura et al., 2013) . We also observe that the CPD accuracy both with respect to the estimated number and the estimated change-point locations decreases with the noise level (see Fig. S3 (b) in the supplement). Fig. 3 (d) shows that, for noise levels in a certain range, the ID and CPOP methods are able to detect all SSEs in over 50% of the simulations. We refer to this noise level range as the suitable noise level (SNL) range. For the simulated signal containing five SSEs, the SNL range for ID and CPOP are 25−47% and 48−85%, respectively. There is no SNL range for the NOT method. Fig. 3 (e) shows the SNL range for simulated signals with different numbers of true change-points (based on the simulated signals shown in the supplement; see Fig. S2 ). We observe that the SNL range depends on the number of true change-points. The NOT method has a much broader SNL range than the CPOP and ID methods when the number of change-points is 8 or less, while no SNL range exists for NOT when the number of change-points exceeds 8. This suggests that NOT does not seem to be suitable for detecting SSEs with short durations. In contrast, SNL ranges for the CPOP method only occur when the number of change-points is 6 or more.

Interestingly, the ID method is the only one among the three methods that exhibits an SNL range for all the simulated SSEs. However, its extent varies depending on the signal.

We also observe that the values for SNLs generally decrease as N increases. When an SSE has a longer duration, the difference between the piecewise-non-linear shape of the SSE signal and a piecewise-linear signal becomes larger. It is sensible that more noise is needed in such cases to cover up the difference between the signals' actual structure and that of a continuous piecewise-linear signal.

We have observed that an SNL range can be found for accurate detection of changepoints in complex piecewise signals such as SSEs using existing CPD methods for continuous piecewise-linear signals. Among the algorithms considered, the ID method seems to have the best behaviour overall when a range of different signals is considered with the number of change-points ranging from 2 to 12. However, the SNL is not consistent for different methods and signal types. Since the noise level and the underlying SSE signal in real-world GPS data are not known, the three CPD methods considered here cannot be directly employed to consistently detect SSEs. Note that these three methods have been shown to outperform other CPD methods for continuous, piecewise-linear signals , which was the reason for only using these methods for the comparison carried out in this section. Motivated by widening SNL ranges, we aim to develop a new algorithm based on the ID method to detect change-points in continuous signals with unknown piecewise structures such as those governing the behaviour of SSEs.

We propose a new method, called singular-spectrum-analysis isolate-detect (SSAID), to determine automatically the SNL for input data with an unknown noise level. This method removes and then adds noise systematically to produce in-SNL data from the input data. Here, the input data refers to the original noisy data on which the CPD is conducted, and the in-SNL data refers to noisy data, the noise level of which is within its SNL range. SSAID is roughly divided into three steps: (1) decomposing the input data into different components by singular spectrum analysis (SSA) and then reconstructing data with different noise levels; (2) adding white noise with various noise levels back to each denoised signal to generate new noisy data, some of which are in-SNL data;

(3) identifying in-SNL data from the new noisy data generated in Step 2, and then outputting the locations of estimated change-points for the input data. The pseudocode of SSAID can be found in Section 6 of the supplement. and Ω 3 are defined in §3.3.

Step 1: Decomposition process SSA is a powerful non-parametric tool for separating underlying signals from the noise, without the need of a priori knowledge of the underlying dynamics (Ghil et al., 2002; Chen et al., 2013; Walwer et al., 2016) . However, it is not designed for detecting changepoints. SSA decomposes the noisy data into different components, and then chooses some of these components in order to reconstruct the signal for the underlying true dynamics. We first use SSA to decompose X t into M components,

The components R j about 1/10 of the time needed by NOT and about 1/600 of the time needed by CPOP for correctly detecting 7 change-points in the slope of a data sequence of length 1408, making ID also the best choice for our framework in terms of computational cost.

We carry out two steps in our improvement scheme: (1) determiningN k,s , the number of estimated change-points in Z k,s t , and (2) identifying the locations of the estimated

For ease of presentation, this set G k,s is called a group. Every realisation Z k,s,m t in this group is called a member and has the same noise level as Z k,s t (see Step 3 in Fig. 4 ).

Firstly, we determineN k,s by a majority voting rule based on the following results.

Let F denote the number of realisations in G k,s with successful detections (see the definition of a successful detection in §2.2). Let P s be the probability that a detection is successful for a given noise level a s . As Z k,s,m t (m = 1, · · · , Q) are independent of each other, the probability that at least half of these detections are successful is

If Z k,s t is an in-SNL data, by the definition of the SNL in §2.2, P s is over 0.5. This gives

is the expectation, and hence P (F ≥ Q/2) will converge to 1 if Q is large enough. For example, P (F ≥ Q/2) is 0.9832 if P s = 0.6 and Q = 100.

We estimate the number of change-points for Z k,s t using the mode of theN k,s,m values, denoted byĤ k,s = M o{N k,s,1 , · · · ,N k,s,Q }, whereN k,s,m is the number of estimated change-points for Z k,s,m t and M o{·} denotes the mode. According to Eq. (5), the

Secondly, we identify the locations of estimated change-points for Z k,s t . For ease of discussion, we call a member of G k,s a qualified member if it satisfiesN k,s,m =Ĥ k,s . Let κ denote the number of qualified members in group G k,s . The locations of the estimated change-points for the j-th qualified member are collected in a vector u k,s,j . All these u k,s,j have the same length; this beingĤ k,s . We store these vectors into a matrix

where I j,i is the location of the i-th change-point for the j-th qualified member, i.e.

u k,s,j = {I j,1 , · · · , I j,Ĥ k,s }, j = 1, · · · , κ. We take the mode of the i-th column in D as the

We verify the validity of the scheme to improve the successful percentage R sd , by conducting tests similar to those of Fig. 3(b) , but now using our method and exploring a much broader range of noise levels. For each noisy data X t in Eq. (1), we generate Q realisations by simulating different noise models, i.e.

where m t is the m-th realisation of t , C wn changes from 1% to 200%, with increments of 1%, and the underlying signal f t is kept unchanged. This set of {X 1 t , · · · , X Q t } is a group for X t , the same as that mentioned before (i.e. G k,s ; also see Step 3 in Fig.   4 ). Following our approach in §2.2, 10, 000 groups are randomly generated for each noise level C wn . Each of these groups contains Q realisations of the white noise. We estimate change-points for each group by using the mode as discussed above. We choose Q = 100 here. Fig. 5 (a) and (b) confirm that our majority voting rule can significantly increase the successful percentage R sd to 100%, when the input data has an SNL. We also calculate the percentage R 1 for each noise level by

where β is the number of groups which satisfy that the mode of the number of estimated change-points,Ĥ, is equal to the number of true change-points (i.e.Ĥ = N = 10 here), and Λ is the total number of groups for the same noise level (i.e. Λ = 10, 000 here). We observe that R sd overlaps with R 1 , which means that the performance of the majority voting rule only depends on the noise level C wn . 

After adding noise, we have generated L×M ×Q new noisy data Z k,s,m t (k = 1, · · · , M ; s = 1, · · · , L; m = 1, · · · , Q) (see Step 3 in Fig. 4) to produce in-SNL data from the input data. However, only some of these noisy data are in-SNL data. Based on the tests conducted in §3.2, we impose three conditions observed to identify in-SNL data: (1)

Here, R 2 refers to the percentage of qualified members (see the definition in §3.2, i.e.N k,s,m =Ĥ k,s ) in a given group (i.e. R 2 = κ/Q), H is the mode of the number of estimated change-points for each group (also seeĤ k,s in §3.2), Ω 3 is the third quartile of the RMSE calculated for each group and v is a pre-defined threshold (see §2.2). The aim of the first condition is to locate in-SNL data. However, R 2 ≥ 50% can occur when the noise level is an SNL or when it is much larger than the SNL range, for which the number of estimated change-pointsN = 0. This situation is demonstrated in Fig. 5(c) . The second condition remedies this pathology. Finally, the third condition aims to remove groups with low accuracy. When calculating the RMSE, the locations of true change-points in the real-world data is estimated through the approach shown in Eq. (6). Members of a group for which the three conditions are met are all in-SNL data. Otherwise, none of them is. For cases in which no change-points are present in the input data X t , no groups have in-SNL data sinceĤ = 0, which means that SSAID will not output any change-points (i.e.N = 0). After obtaining all identified in-SNL data, the scheme for improving R sd described in §3.2 is adopted to determine the locations of change-points for the input data (for details, see the pseudocode of SSAID in the supplement). The Type I error of SSAID is the same as that of ID, which has been considered under the null hypothesis of no change-points in .

In this section, we demonstrate the performance of our method for a range of simulated SSE data X t , and discuss factors affecting the detection performance. For verification of the general applicability of our method to detect change-points in other types of signals with known piecewise structures, including piecewise-linear and piecewise-non-linear (i.e.

piecewise-quadratic, piecewise-exponential and sinusoidal), see the additional simulation studies in the supplement.

The SSE data are generated in the form of Eq. (1), in which f t is the pure SSE signal shown in Fig. 2 (a) , and the noise level C wn varies from 0 to 100%, with increments of 1%.

We also considered f t from the other simulated pure SSE signals in §2, and the results are included in the supplement (see Fig. S5 ) as they are consistent with the current tests. We create 100 data sequences of independent standard Gaussian random variables t (t = 1, 2, . . . , T ). In total, we have 100 × 101 noisy time series X t (t = 1, 2, . . . , T ). Fig. 6 (a) shows the number of estimated change-points for each time series. SSAID correctly estimates the number of true change-points in over 80% of all cases analysed. In particular, the number of estimated change-points is correct for all the cases with noise levels lower than 40%. Fig. 6 (b) shows that the percentage of successful detections R sd (see Eq.

(2)) and the percentage R 1 (see Eq. (8)) are different. They are both 100% when C wn ≤ 25%, and decrease with C wn , when C wn > 25%. The percentage of successful detections R sd decreases faster than R 1 when C wn increases, indicating that the accuracy of the detected change-point locations fades with increasing C wn values. Fig. 6 (c) demonstrates the high accuracy of the change-points detected using our method for data with a low noise level. Fig. 6 (d) shows that when the noise level is very high (C wn = 100%), the locations of some detected change-points is not as accurate. The effect of the noise level C wn on the performance of our method comes from a deficiency of SSA, which generally fails to distinguish the underlying signal from the noise when the signal-to-noise ratio (SNR) of the input data is too low.

We now compare the detection performance of SSAID with a commonly-used method in the field of geosciences for detecting SSEs in GPS data, i.e. linear regression combined with Akaike's Information Criterion (AIC) (Nishimura et al., 2013) . This method (1) uses a sliding window of a fixed width and moves the window at each time point; (2) fits a linear model to the data in the window; (3) divides the data in the window into equal halves and fits a linear model to each half, and (4) calculates the AIC difference (i.e. ∆AIC) between the single linear model and the two linear models at each middle point of the window. If this midpoint of a window is a change-point, e.g. the start-or end-point of an SSE, the two linear models fit the observational data better than a single linear model, thus resulting in a negative ∆AIC. As a negative ∆AIC does not always correspond to change-points of SSE signals, we must specify an appropriate threshold, denoted by ζ, in order to detect change-points of SSEs. If ∆AIC is lower than ζ, its corresponding time is regarded as a change-point. We take a sliding time window of 15 days to calculate ∆AIC for each data point (except the first and last 7 observations) in the simulated SSE data in Fig. 6 (c) and (d). Fig. 6 (e) and (f) show the changepoints detected using three thresholds (high, medium and low). We observe that none of them succeeds in finding all the true change-points accurately. When ζ is too high, only the most significant SSEs can be detected, while for smaller ζ, the detection generally overestimates the number of change-points. The selection of the threshold value depends on the signal itself, making it impossible to detect all the change-points in multiple time series or even within a single time series by using a single threshold. Also, the selection of the sliding window width is subjective. The test confirms that our approach outperforms the linear regression method.

Three key parameters may affect the performance of SSAID: the number of decomposed components M in SSA, the number of realisations Q, and the highest level L of added Gaussian noise (see Eq. (4)). The first parameter M comes from the well-developed SSA algorithm, and has been widely discussed in the literature (e.g. Ghil et al., 2002; Walwer et al., 2016) . Based on these studies, M = 100 is a reasonable choice. We will mainly focus on the selection of L and Q. We seek to guarantee the existence of in-SNL data among the Z k,s,m t (see Step 3 in Fig. 4 ) while mitigating computing time.

We can see from Eq. (5) that the percentage of successful detections R sd increases with Q if P s is fixed. We reproduce the tests shown in Fig. 5 (a) , but using different Q values. Fig. 7 (a) and (b) show that the SNL range varies little for Q ≥ 30 and R sd reaches 100% for in-SNL data. To ensure convergence, we take Q = 50 in our subsequent tests. As discussed in §3.2, L should be sufficiently large so that we can obtain in-SNL data (see Fig. 4 ); meanwhile, L should not be too large to limit the computation time. We reproduce the tests from Fig. 6 (a) and (b), with L = 10, 20, · · · , 100. Fig. 7 (c) shows that the dependence of R sd on C wn converges rapidly with L for L ≥ 30. We choose L = 80 in our tests of current simulated SSE data, which is large enough to guarantee the existence of in-SNL data.

The SSE data were obtained from the Scripps Orbit and Permanent Array Center 

Based on numerical studies in §4, we choose M = 100, Q = 50 and L = 80 for the current applications. The detected change-points are indicated by blue vertical lines in Fig.8 . As opposed to simulated SSE data, we do not have the information on when an SSE starts and ends for real data sets. However, we can utilize the tremor catalogue to check if the detected change-points coincide with possible SSEs or not, since tremors often accompany SSEs. An increasing number of tremors generally indicates that an SSE is probably occurring (Ito et al., 2007) . We include the total number of tectonic tremors per day in Cascadia in the bottom row of Fig. 8 (data from Idehara et al., 2014) .

Four obvious peaks in the number of tremors imply four possible SSEs (labelled as SSE1, SSE2, SSE3 and SSE4, shaded areas in Fig. 8 ). Associated SSEs detected using SSAID in different SSE data are highlighted by the same color. Our current detections succeed in detecting change-points in a subset of the GPS data for all 4 events. Note that there may be a time lag of up to 0.05 years (i.e. about 18 days) between SSEs detected using GPS data and SSEs indicated by tremor events (Crowell et al., 2016; Wang et al., 2018) .

In addition, it is not expected that all the SSEs can be recorded in a single GPS station.

The best detection was for the fourth SSE. That event lasted for more than 1 month and caused significant ground movements, which was detected in a number of other studies as well (Crowell et al., 2016) . However, past studies didn't detect the other SSEs that our method could detect. These SSEs have a much shorter duration of the order of a few days and cause small ground displacements, explaining why the change-points can only be detected at some GPS stations.

In the aforementioned numerical studies, the underlying signals have a regular periodicity, and the amplitude jump in each cycle remains nearly constant. In real-world data, however, the amplitude jumps at different change-points could differ. For example, in Fig. 1 , the amplitude jump caused by the SSE indicated by the red arrow is much smaller than that indicated by the yellow arrow. We also notice such an obvious amplitude jump difference in our current SSE data (e.g. the shaded-yellow and shaded-blue areas for the GPS station ALBH in Fig. 8 ).

We now introduce a sliding window into SSAID to improve its adaptability to real SSE signals with more complicated structures. First, we reproduce the tests from Fig. 6 (a), with 2 newly generated pure SSE data (both contain two SSEs, see Fig. 9 (a)), to investigate the performance of SSAID in such situations. The amplitude jumps caused by the two SSEs in the first simulated SSE data are very different with each other, but those in the second SSE data are similar to each other. Dotted lines in Fig. 9 (b) show that SSAID can consistently detect change-points in the second SSE data with a high R sd > 80% for all noise levels C wn , while it fails to detect true change-points in the first SSE data in most cases. This observation can be explained by the SNL range for the two SSE data (see Section 7 in the supplement).

We introduce a sliding window into SSAID to overcome this issue. We first divide the input data into several non-overlapping segments of equal length denoted as L s , and then apply SSAID to detect change-points in three consecutive segments iteratively. We refer to the three consecutive segments as a sliding window, with the length L w equal to 3L s .

For example, if we choose the segment length L s = 80 for the SSE data in Fig. 9 (a), we divide the input data into 6 segments, separated by the vertical dotted black lines.

We then have 4 sliding windows, which consist of segment 1 − 3, segment 2 − 4, segment 3 − 5 and segment 4 − 6, respectively (see the top of Fig. 9 (a) ). SSAID conducts the change-point detection on the data in these sliding windows in turn. Fig. 9 (b) illustrates that the introduction of a sliding window drastically increases the successful detection ratio R sd from 0% to over 80% when C wn ≤ 15%, but fails to improve the successful detection ratio R sd when C wn ≥ 25%. When the noise level of the first SSE data is too large, the amplitude jump in the red subset will be fully covered up by the noise, so that SSAID cannot identify change-points, as expected. Fig. 9 (b) also shows that, the introduction of a sliding window has limited influence on the performance of SSAID for the second SSE data, as in-SNL data exist without the implementation of the sliding window approach. The current tests confirm that the introduction of a sliding window in SSAID can help us deal with more complicated cases when the amplitude jumps at different change-points greatly vary.

We now investigate the effect of the segment length L s on the detection performance, by reproducing the tests shown in Fig. 9 (b), but with a range of L s . The results demonstrate that the L s value should be in the range 80 ≤ L s ≤ 100 for the first SSE data, as small L s < 80 (large L s > 100) overestimates (underestimates) the number of change-points. For the second SSE data, L s should be larger than 80 (see Fig. S8 in the supplement).

We then apply SSAID with a sliding window to detect change-points in the observed SSE data from the previous section. We choose L s = 80 for these observed data. Fig. 8 shows that the new results (see red dotted lines) mostly overlap with those detected by SSAID without a sliding window. However, the new approach detects more accurately the start and end times of SSE2 and SSE3 at KTBW station and that of SSE4 at UCLU station than the method without a sliding window. It is well known that there are longterm and short-term SSEs which cause different amplitude jumps in GPS data, which explains why a sliding window could be preferable. 

We have presented a method called SSAID for detecting change-points in noisy data sequences with an unknown continuous piecewise structure. This method creates in-SNL data through a decomposition process and a follow-up process systematically adding white noises, so that we can then apply Isolate-Detect, which is an existing method for continuous piecewise-linear signals, to detect change-points in the newly constructed noisy data. SSAID exhibits very good performance in both simulated and real data sets.

We verify the validity of SSAID on various simulated SSE data, and apply it to detect SSEs in observed SSE data. We further compare its detection performance with a commonly-used method (i.e. linear regression combined with AIC). Our key findings are: (1) Our method is able to consistently detect the true change-points in all the simulated SSE data; (2) SSAID outperforms the linear regression combined with AIC when detecting change-points in SSE data, particularly in terms of detection accuracy We also investigated the general adaptability of SSAID to detect change-points in other fields in the supplement, including simulated data with various piecewise structures (i.e. piecewise-quadratic, piecewise-exponential, sinusoidal function and piecewiselinear) and real data sets from other fields (i.e. COVID-19 daily confirmed cases in the US and the monthly S&P 500 close price index). We observe that SSAID can obtain as good performance in both simulated and real data sets as that for the SSE data.

SSAID can also be applied to detect change-points in piecewise-linear signals, with a more competent detection accuracy to that of Isolate-Detect.

We conclude that SSAID can be used to detect change-points in data from various disciplines, as it does not need prior information of the underlying dynamics, while the development of an objective method for selecting the length of the sliding window for SSAID is also needed in the subsequent work.

We consider a simplified Cascadia-like model, which consists of a planar frictional interface simulating the thrust fault between a subducting oceanic plate and an overlying continental plate in an elastic half-space, as shown in Fig. S1 . It is composed of a seismogenic zone where megathrust earthquakes spontaneously arise, a transition zone where SSEs spontaneously arise, a creeping zone where instabilities cannot spontaneously arise, and an infinitely long zone which is enforced to slide at a constant speed (i.e. the plate velocity v pl ). We use the following dimensionless system of ODEs to simulate the slip velocity history for the i-th subfault,

where indices i and j refer to different subfaults; N is the number of the subfaults in our model; K ij is the kernel coefficient, which quantifies the interactions between the i-th and j-th subfault; v i and θ i are the slip rate and a state variable for the i-th subfault, respectively; a and b are experimentally determined constitutive constants. This system is solved by using the adaptive fourth order Runge-kutta method ( see more details in ). Except for the SSE data shown in the main context, we simulate more SSE data with different number of SSEs in a one-year period, shown in Fig. S2.   Fig. S1 . Geometry of a 2D simplified subduction fault, which is composed of a seismogenic zone, a transition zone (velocity-weakening, where SSEs should occur), a creeping zone (velocity-strengthening) and an infinitely long zone forced to slide at the constant velocity of v plate .

When defining a successful detection in §2.2, the threshold v should be designated by the users according to their own fields. Fig. S3 shows the dependence of R sd on different v.

We will obtain a higher R sd if adapting a larger v. To investigate the effect of noise levels on the detection accuracy of Isolate-Detect (ID) , we plot the quartile distributions of RMSE for each noise level in Fig. S3(b) . We first identify the successful detections among the corresponding 10, 000 tests for each noise level in the tests shown in Fig. 3(b) of the main context, and then calculate the RMSE for each detection. The quartile values are picked from these calculated RMSE. 

We verify the validity of taking the mode of each column in D shown in Eq. (6) to estimate the locations of underlying change-points in noisy data sequences. We first recall tests for a suitable noise level (i.e. C wn = 35%) shown in Fig. 3 (b) as an example. It is sensible that the mode of the number of detected change-points is equal to the number of true change-points (i.e. 10). Among these 10, 000 tests, we identify the ones withN = 10, and store the locations of their detected change-points to form the matrix D shown in Eq. (6). The final output of estimated change-points for this noise level is (20, 26, 94, 99, 165, 173, 239, 246, 316, 321) , which is consistent with the locations of true change-points (i.e. 19, 26, 93, 100, 166, 173, 240, 246, 313, 320) . The histogram for each estimated change-point is shown in Fig. S4 . The current result confirms that the mode is suitable to estimate the location of true change-points.

We further demonstrate the performance of our method based on the same tests shown in Fig. 6 (a) in the main context, but using different underlying SSE signal f t . The results are shown in Fig. S5 , which are consistent with the observations in the main text (see §4.1.1).

We first verify the general applicability of SSAID to detect change-points in other types of piecewise-non-linear data sequences using simulated signals with known structures. We consider signals with the following structures: (1) piecewise-quadratic;

(2) piecewiseexponential; and (3) sinusoidal (see Fig. S6 (a), (c) and (e), respectively). We have 100 × 101 noisy data sequences for each scenario. Fig. S6 (b) and (f) show that for the piecewise-quadratic and the sinusoidal signals, SSAID consistently detects the true change-points successfully when C wn is sufficiently small (approximately ≤ 20% for piecewise-quadratic; ≤ 40% for sinusoidal funtion). SSAID can guarantee successful detection of true change-points, while the successful detection percentage monotonically decreases with C wn . Our approach performs even better for the piecewise-exponential signal (see Fig. S6 (d) ). It finds the correct number of change-points even when C wn is close to 100%. We also check its adaptability to piecewise-linear signals, and further compare its performance with that of ID. Fig. S6 (h) illustrates SSAID can accurately detect true change-points, and improve the R sd of ID for C wn ≥ 60%. These four cases indicate that our approach can be applied to a wide variety of signals.

Algorithm 1: The pseudocode of SSAID Result: Estimated change-points for the input data X t 1

Step 1 (Decomposition process): Obtain denoised data with different noise levels Y k t (k = 1, · · · , M );

Step 2 (Adding noise in the way of the R sd improvement scheme): Generate a range of new noisy data Z k,s,m t (k = 1, · · · , M ; s = 1, · · · , L; m = 1, · · · , Q) to guarantee the existence of in-SNL data; Condition 3:

27 Ω 3 ← The third quartile (75%) of the RMSE for each group by assuming U k,s as the real change-points for Z k,s t ;

if Ω 3 > v then 29 All the members in the current group G k,s are marked as 'NOT in-SNL data'; Output the number of change-points for Y k t :Ĥ k ← M o{Ĥ k,1 , · · · ,Ĥ k,L }; 35 end We further elaborate possible reasons that could explain the detected change-points, which can often be associated with the issued dates of government interventions (e.g. lockdown, mass testing or lifting restrictions). Detecting these change-points can further help monitor COVID-19 phase changes, e.g. in terms of its growth rate, which can assist the government assess the effectiveness of their policies on the evolution of the pandemic (Jiang et al., 2020; Pei et al., 2020) . The interventions described in the following context mainly come from the timeline of the COVID-19 pandemic in the United States (2020) in Wikipedia (see https://en.wikipedia.org/wiki/Timeline_of_the_ COVID-19_pandemic_in_the_United_States_(2020)#cite_note-121) and the report from Schuchat et al. (2020) . The first change-point is detected on 9 March 2020, which might result from the order issued by the CDC (Centers for Disease Control and Prevention) on 5 March. This order lifted the restrictions for a coronavirus test so that anyone in the US was allowed to get an COVID-19 test, resulting in a rapid increase. This change-point is also consistent with the syndromic surveillance in counties affected early by the COVID-19 in the US, which did not show an increase in visits for COVID-19-like illness (e.g. fever, cough) until 28 February 2020. We find the second change-point on 02

April 2020, which coincides with the emergent orders from New York state government (i.e. the leading confirmed cases in the US at that time. The New York state government closed all nonessential construction sites on 28 March 2020, and the U.S. Navy medical ship arrived in the New York city to assist with non-COVID operations, relieving land hospitals to stop the city's growing COVID-19 pandemic, on 30 March 2020 (see https://en.wikipedia.org/wiki/COVID-19_pandemic_in_New_York_(state)). The third change-point is detected on 14 June 2020, the same day as the date when the House Appropriations Committee required the passengers to wear facial masks on public transportation. On 17 July 2020, the fourth change-point, the highest single-day increase in new confirmed cases worldwide was recorded in the US at that time. More possible explanations about detected change-points can be found in the timeline of COVID-19 in the US (2020) (see https://en.wikipedia.org/wiki/Timeline_of_the_COVID-19_ pandemic_in_the_United_States_(2020)#cite_note-121).

We finally investigate the performance of SSAID on detecting change-points in the monthly S&P 500 close price index. These change-points can be assosiated with the connecting points between the bull and bear stock markets, at which the previous trend shows either a unidirectional increase or decrease, while the later trend reverses to an 

Detecting multiple generalized change-points by isolating single ones

Narrowest-over-threshold detection of multiple change points and change-point-like features

Greedy automatic signal decomposition and its application to daily gps time series

Singular spectrum analysis for modeling seasonal signals from gps time series

Data segmentation algorithms: Univariate mean change and beyond

Single-station automated detection of transient deformation in gps time series with the relative strength index: A case study of cascadian slow slip

Detecting changes in slope with an l 0 penalty

Advanced spectral methods for climatic time series

Determination of slow slip episodes and strain accumulation along the cascadia margin

Regional and global variations in the temporal clustering of tectonic tremor activity

Slow earthquakes coincident with episodic tremors and slow slip events

Slow slip events controlled by the slab dip and its lateral change along a trench

Detection of short-term slow slip events along the nankai trough, southwest japan, using gnss data

Connecting slow earthquakes to huge earthquakes

Episodic slow slip events and rate-and-state friction

Selective review of offline change point detection methods

Slow slip events in new zealand

Data-adaptive detection of transient deformation in geodetic networks

Identifying the recurrence patterns of nonvolcanic tremors using a 2-d hidden markov model with extra zeros

2020) A review on minimax rates in change point detection and localisation

37 Collect all the groups which are not marked as 'NOT in-SNL data

38 Step 4: Output the final estimated locations of change-points for X t

39 Pick up all the members withN k,s,m =N from in-SNL data, and then store their detected change-points into a new matrix D f , which is the same as D in Eq

U ← The mode of each column in D f is used to estimate the location of change-points in X t

The SNLs for the green and red subsets have no common values. This explains why no SNL exists for the first SSE data. The SNLs for the golden and cyan subsets partially overlap in the noise levels 32 − 152%

We now apply SSAID to detect change-points in other real-world data sets with piecewisenon-linear structures

COVID-19 daily new confirmed cases in the US The first case study here is to look for change-points in COVID-19 daily new confirmed cases in the US. We perform the Anscombe transform on the original data, to make its distribution more consistent with the assumption of Gaussianity for the noise (Anscombe, 1948). Fig. S7 (a) shows that the detected change-points by SSAID with and without a sliding window mostly overlap, except two change-points indicated by the shaded green areas

As described in §4.4, our method cannot find an SNL for this kind of data set without a sliding window, thus failing to detect any change-points. The detected change-points by our method with a sliding window of 40 observations are mostly consistent with those detected by

We investigate the effect of the segment length on the detection performance, by reproducing the tests shown in Fig. Fig. 9 (b), but with a range of L s . The results are shown in Fig

Detecting multiple generalized change-points by isolating single ones

The transformation of poisson, binomial and negative-binomial data

An analysis of bull and bear markets in the us and norway. Master's thesis

Rhetoric, risk, and markets: The dot-com bubble

Time series analysis of covid-19 infection curve: A change-point perspective

The global financial crisis: Causes and consequences

Differential effects of intervention timing on covid-19 spread in the united states

Episodic slow slip events and rate-and-state friction

Public health response to the initiation and spread of pandemic covid-19 in the united states

Predicting the turning points of a time series