Time-Varying GPS Displacement Network Modeling by Sequential Monte Carlo

Abstract

Geodetic observations through high-rate GPS time-series data allow the precise modeling of slow ground deformation at the millimeter level. However, significant attention has been devoted to utilizing these data for various earth science applications, including to determine crustal velocity fields and to detect significant displacement from earthquakes. The relationships inherent in these GPS displacement observations have not been fully explored. This study employs the sequential Monte Carlo method, specifically particle filtering (PF), to develop a time-varying analysis of the relationships among GPS displacement time-series within a network, with the aim of uncovering network dynamics. Additionally, we introduce a proposed graph representation to enhance the understanding of these relationships. Using the 1-Hz GEONET GNSS network data of the Tohoku-Oki Mw9.0 2011 as a demonstration, the results demonstrate successful parameter tracking that clarifies the observations’ underlying dynamics. These findings have potential applications in detecting anomalous displacements in the future.

1. Introduction

The Global Positioning System (GPS) provides permanent static displacement information which is useful in the verification of physically based models in the study of tectonic and volcanic systems [1] and complements seismological data in earthquake-related studies. It can aid in determining earthquake rupture geometry [2], estimating the time-varying distribution of fault slip [3], assessing earthquake magnitude for early warning systems [4], and detecting ground motion caused by earthquakes [5,6]. Additionally, GPS data are used to calculate velocity fields, contributing to the description of crustal deformation in various regions [7,8,9].

As the significance of GPS observations in earthquake-related studies and crustal deformation analyses becomes evident, the spatiotemporal dynamics of surface displacement from GPS network data become crucial for effective geological hazard assessment and mitigation. Previous studies have explored correlations between seismic activity and surface deformations [10], used machine learning frameworks to integrate the spatiotemporal dependencies of GPS displacements for landslide displacement prediction [11], and employed the spatiotemporal fields of GPS time-series for earthquake prediction [12].

Despite these efforts, there has been minimal attention paid to modeling the relationships among GPS measurements in a network and tracking network dynamics. This study builds upon the potential of high-rate GPS time-series data, employing the sequential Monte Carlo method, specifically particle filtering (PF), to develop a time-varying analysis of the relationships among GPS displacement time-series within a network. The aim is to uncover the network dynamics and enhance the understanding of these relationships through a proposed graph representation. Our focus is on utilizing the 1-Hz GEONET GNSS network data of the Tohoku-Oki Mw9.0 2011 earthquake as a demonstration, with results highlighting the potential of our approach for anomalous displacement detection and geological hazard assessment in the future.

2. GPS Data

The post-processed GPS records used in this study obtained from [13], originally composed of 1-Hz GPS displacement data from 847 GEONET [14] stations in Japan in the north–south, east–west, and up–down components. The dataset covers a period prior, during, and after the Tohoku-Oki Mw9.0 earthquake in Japan on 11 March 2011, and it was originally used in a study to detect significant ground motion from a GPS data network [6], from which the data were also shown to reliably express the ground motion caused by the propagation of seismic motion. For an illustrative purpose, a subset of these data was selected for the experiment, as discussed in Section 5.1. Additionally, the location coordinates of the GPS stations were separately obtained from [15].

3. Model

In this study, a network consists of a subset of N GPS stations. The displacement observed at any station at time epoch t is assumed to be related to the displacement observed at the previous time epoch $t-1$ of all stations including itself. A simple linear relationship for an observation at the ith GPS station is assumed as follows:

(1)$\begin{array}{cc}\hfill x_{i,t}& =a_{i1,t}{x}_{1,t-1}+a_{i2,t}{x}_{2,t-1}+\cdots +a_{iN,t}{x}_{N,t-1}+{\eta }_{i,t}\hfill \\ \hfill & ={\mathit{a}}_{i,t}·{\mathit{x}}_{t-1}+{\eta }_{i,t},\hfill \end{array}$

where $x_{i,t}$ denotes an observation of the ith GPS station at time epoch t. The vector ${\mathit{x}}_{t-1}$ denotes the observations of all GPS stations in a network at time epoch $t-1$. $a_{ij,t}$ are coefficients of the linear equation, which we want to recover. These coefficients reflect the influence of the previous observation at station j on the current observation at station i. Moreover, they are time-varying and can be of different values at a different time epoch t. Lastly, ${\eta }_{i,t}$ are noise terms.

Observation Equation (1) can be written in a vector form as

(2)$\begin{array}{cc}\hfill {\mathit{x}}_t& ={\mathit{A}}_t{\mathit{x}}_{t-1}+{\mathit{\eta }}_t\hfill \\ \hfill \left[\begin{array}{c}{x}_{1,t}\\ x_{2,t}\\ ⋮\\ x_{N,t}\end{array}\right]& =\left[\begin{array}{cccc}{a}_{11,t}& a_{12,t}& \cdots & a_{1N,t}\\ a_{21,t}& a_{22,t}& \cdots & a_{2N,t}\\ ⋮& ⋮& \ddots & ⋮\\ a_{N1,t}& a_{N2,t}& \cdots & a_{NN,t}\end{array}\right]\left[\begin{array}{c}{x}_{1,t-1}\\ x_{2,t-1}\\ ⋮\\ x_{N,t-1}\end{array}\right]+\left[\begin{array}{c}{\eta }_{1,t}\\ {\eta }_{2,t}\\ ⋮\\ {\eta }_{N,t}\end{array}\right],\hfill \end{array}$

where it can be noted that each ith row of ${\mathit{A}}_t$, the denoted ${\mathit{a}}_{i,t}=\left(\begin{array}{cccc}{a}_{i1,t}& a_{i2,t}& \cdots & a_{iN,t}\end{array}\right)$ is a hidden state vector for an observation of the ith GPS station at time epoch t, or $x_{i,t}$. Later, in Section 5.4, a graph representation based on the recovered ${\mathit{A}}_t$ for each time epoch is introduced.

The coefficients $a_{ij,t}$ in Equation (1) are assumed to have linear transitions from time epoch $t-1$ to time epoch t as follows:

(3)$a_{ij,t}=a_{ij,t-1}+v_{ij,t},$

where $v_{ij,t}$ are noise terms.

For simplicity, the state noise terms $v_{ij,t}$ in Equation (3) are assumed to be i.i.d. (independently and identically distributed) Gaussians. Conversely, the GPS observation noises ${\eta }_{ij,t}$ in Equation (1) have been shown not to necessarily follow Gaussian distributions [16,17,18]. These observation noises might exhibit heavier tails or other non-Gaussian characteristics. In this study, we model these observation noises ${\eta }_{ij,t}$, using three different distributions, i.i.d. Gaussian, i.i.d. Laplace, and i.i.d. Cauchy, and we present results for each. However, it can be noted that these noise terms can be modeled using other kinds of noise such as Gaussian mixtures [19], and alpha-stable [20] distributions, depending on the specific application requirements. This flexibility in choosing the noise terms allows for more accurate and tailored representations in various application scenarios.

Equations (1) and (3) define a model [21] widely used in time-series models, which has been applied in various fields such as computational biology [22], biogeophysics [23], brain connectivity [24], and petroleum science [25].

More precisely, our model is defined by two stochastic processes in the forms

(4)${\mathit{a}}_t=f_t\left({\mathit{a}}_{t-1},{\mathit{v}}_t\right),$

(5)${\mathit{x}}_t=h_t\left({\mathit{a}}_t,{\mathit{\eta }}_t\right),$

where a state Equation (4) represents a process [26], in which a hidden parameter vector at time t depends on that of the previous time instant $t-1$. Equation (5) is an observation equation which is related to the hidden parameter vector ${\mathit{a}}_t$ of the state equation. ${\mathit{v}}_t$ and ${\mathit{\eta }}_t$ are noise terms. The intuition of the model is to capture the behaviors of an observation vector ${\mathit{x}}_t$ in terms of an unobserved state vector ${\mathit{a}}_t$.

In the case of a linear Gaussian model, where $f_t$ and $h_t$ are linear functions and the noise terms ${\mathit{v}}_t$ and ${\mathit{\eta }}_t$ are normally distributed, one of the classical methods to solve this problem is the Kalman Filter [27]. However, the aim of this study is to incorporate any existing information and model hidden parameters across any underlying distribution. Furthermore, our approach can be generalized to nonlinear observation or state equations, thereby offering enhanced flexibility and applicability across a broader spectrum of scenarios.

4. Sequential Monte Carlo

We propose to apply the sequential Monte Carlo method or particle filtering, a Bayesian method based on an importance sampling and resampling technique. This method is used to compute the posterior distributions of the hidden parameters, while it also allows the utilization of prior information. Importantly, this method allows nonlinearities and non-Gaussian noises in the state and observation equations. This offers flexibility to the modeling of geophysical phenomena, which may not always follow a Gaussian distribution, and deviations from the normal distribution can influence actual dynamics [28,29].

More precisely, in this study, a sequential Monte Carlo method or particle filtering (PF) is used to sequentially find the following posterior of the hidden parameter vector at each time epoch t, according to Bayes’ rule:

(6)$p\left({\mathit{a}}_{i,t}|{\mathit{x}}_{1:t}\right)=\frac{p\left({\mathit{x}}_t|{\mathit{a}}_{i,t}\right)p\left({\mathit{a}}_{i,t}|{\mathit{x}}_{1:t-1}\right)}{p\left({\mathit{x}}_t|{\mathit{x}}_{1:t-1}\right)},$

where ${\mathit{x}}_{1:t}$ denoted observations at all GPS stations from time epochs 1 to t, while ${\mathit{x}}_t$ denoted the observations at time epoch t. Recall that ${\mathit{a}}_{i,t}=\left(\begin{array}{cccc}{a}_{i1,t}& a_{i2,t}& \cdots & a_{iN,t}\end{array}\right)$ is a hidden parameter vector, which we want to recover, for an observation at the ith GPS station at time t or $x_{i,t}$.

For a Gaussian observation noise assumption, an observation has the following likelihood:

(7)$p\left(x_{i,t}\right|{\mathit{a}}_{i,t})=\frac{1}{{\left(2\pi {\sigma }_{\eta }^2\right)}^{1/2}}exp\left(-\frac{{\left(x_{i,t}-{\widehat{x}}_{i,t}\right)}^2}{2{\sigma }_{\eta }^2}\right),$

where ${\widehat{x}}_{i,t}$ is derived from ${\mathit{a}}_{i,t}$ and ${\mathit{x}}_{t-1}$ through Equation (1), and ${\sigma }_{\eta }$ is a standard deviation of the observation noise.

For a Laplace observation noise assumption, an observation has the following likelihood:

(8)$p\left(x_{i,t}\right|{\mathit{a}}_{i,t})=\frac{1}{2\beta }exp\left(-\frac{\left|x_{i,t}-{\widehat{x}}_{i,t}\right|}{\beta }\right),$

where ${\widehat{x}}_{i,t}$ is derived from ${\mathit{a}}_{i,t}$ and ${\mathit{x}}_{t-1}$ through Equation (1). $\beta >0$ is a scale parameter, and $\sqrt{2}\beta$ is a standard deviation of the observation noise.

For a Cauchy observation noise assumption, an observation has the following likelihood:

(9)$p\left(x_{i,t}\right|{\mathit{a}}_{i,t})=\frac{1}{\pi \gamma \left[1+{\left(\frac{x_{i,t}-{\widehat{x}}_{i,t}}{\gamma }\right)}^2\right]},$

where ${\widehat{x}}_{i,t}$ is derived from ${\mathit{a}}_{i,t}$ and ${\mathit{x}}_{t-1}$ through Equation (1), and $\gamma >0$ is a scale parameter that determines the distribution’s spread.

Equation (6) provides the optimal Bayesian solution for the hidden parameters for Equations (4) and (5). However, the denominator in Equation (6) is intractable, and the solution often cannot be determined [30]. Particle filtering solves for the solution of the model in Equations (4) and (5) via a sampling scheme. It provides a Monte Carlo approximation for the posterior in Equation (6), using a finite number M of weighted samples or particles:

(10)$p\left({\mathit{a}}_{i,t}|{\mathit{x}}_{1:t}\right)\approx \sum _{m=1}^M{w}_{i,t}^{\left(m\right)}\delta \left({\mathit{a}}_{i,t}-{\mathit{a}}_{i,t}^{\left(m\right)}\right),$

where ${\mathit{a}}_{i,t}^{\left(m\right)}$ are particles, $w_{i,t}^{\left(m\right)}$ are their weights, and $\delta$ denotes the delta-Dirac function, which concentrates probability density at the particles. As the number of particles, M, grows and tends toward infinity, the accuracy of the approximation improves and converges towards the true distribution.

More precisely, at any time epoch t, the algorithm has a set of filtering particles $\{{\mathit{a}}_{i,t-1}^{\left(m\right)},w_{i,t-1}^{\left(m\right)}\}m=1\dots M$, which represent samples from the previously estimated posterior distribution $p\left({\mathit{a}}_{i,t-1}\right|{\mathit{x}}_{1:t-1})$. To estimate the posterior $p\left({\mathit{a}}_{i,t}\right|{\mathit{x}}_{1:t})$ in a current iteration, we choose to sample from a proposal distribution q, which is perhaps convenient to sample from and approximates the target posterior distribution in some sense:

(11)${\mathit{a}}_{i,t}^{\left(m\right)}\sim q\left({\mathit{a}}_{i,t}\right|{\mathit{a}}_{i,t-1}^{\left(m\right)},{\mathit{x}}_t).$

To ensure that particles approximate samples from the target distribution, the algorithm utilizes the sequential importance sampling method [30], where weights assigned to particles are determined by a correction factor: $p/q$. This is to adjust more weights to particles from critical regions, effectively reducing the overall sampling variance of the estimator. Furthermore, this particular sampling method requires fewer samples compared to alternative methods such as rejection sampling. More precisely, the importance weight [26] of a particle ${\mathit{a}}_{i,t}^{\left(m\right)}$ is assigned as

(12)$w_{i,t}^{\left(m\right)}\propto \frac{p\left({\mathit{a}}_{i,t}^{\left(m\right)}\right|{\mathit{x}}_{1:t})}{q\left({\mathit{a}}_{i,t}^{\left(m\right)}\right|{\mathit{a}}_{i,t-1}^{\left(m\right)},{\mathit{x}}_t)},$

which, to avoid recalculation when new data arrives, is equivalent to the following sequential update [26]:

(13)$w_{i,t}^{\left(m\right)}\propto w_{i,t-1}^{\left(m\right)}\frac{p\left(x_{i,t}\right|{\mathit{a}}_{i,t}^{\left(m\right)})p\left({\mathit{a}}_{i,t}^{\left(m\right)}\right|{\mathit{a}}_{i,t-1}^{\left(m\right)})}{q\left({\mathit{a}}_{i,t}^{\left(m\right)}\right|{\mathit{a}}_{i,t-1}^{\left(m\right)},{\mathit{x}}_t)}.$

The proposal distribution, q, should be selected based on the characteristics of the problem and the target distribution. The popular choice is a bootstrap filter [31], which uses the state transition density as the proposal distribution, namely to let $q\left({\mathit{a}}_{i,t}^{\left(m\right)}\right|{\mathit{a}}_{i,t-1}^{\left(m\right)},{\mathit{x}}_t)=p\left({\mathit{a}}_{i,t}^{\left(m\right)}\right|{\mathit{a}}_{i,t-1}^{\left(m\right)})$. This results in a simplified weight update, requiring only the likelihoods as follows:

(14)$w_{i,t}^{\left(m\right)}\propto w_{i,t-1}^{\left(m\right)}·p\left(x_{i,t}\right|{\mathit{a}}_{i,t}^{\left(m\right)}).$

The particle weights in Equation (14) are then normalized so that ${\sum }_{m=1}^M{w}_{i,t}^{\left(m\right)}=1$ to ensure that the weighted samples represent a valid probability distribution for the estimation of the posterior in Equation (6). The normalization [30] is as follows:

(15)$w_{i,t}^{\prime \left(m\right)}=\frac{w_{i,t}^{\left(m\right)}}{{\sum }_{m=1}^M{w}_{i,t}^{\left(m\right)}}.$

The final weighted samples ${\{{\mathit{a}}_{i,t}^{\left(m\right)},w_{i,t}^{\prime \left(m\right)}\}}_{m=1...M}$ represent samples which estimate the posterior distribution in Equation (6).

It is important to note that in high-dimensional state spaces, it can be difficult to sample particles that adequately cover the state space. This limited number of particles may struggle to represent the target distribution accurately, leading to particle weights becoming concentrated on a few particles. This problem, known as degeneracy, can be resolved by resampling [32], which involves replicating particles with higher weights, and removing particles with lower weights. This prevents the algorithm from being dominated by a few particles. Typically, the resampling step is triggered when $N_{eff}=\frac{1}{{\sum }_{m=1}^M{\left(w_t^{\left(m\right)}\right)}^2}$ is below a user-set threshold [26].

The particle filtering method employed in this study is summarized in Algorithm 1. It can be noted that this algorithm is applicable to real-time data. Additionally, the first set of particles are generated from a prior distribution which represents an initial belief or knowledge about the possible states of a system.

Algorithm 1 GPS Displacement Network Learning

Input: $X=[{\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_T]\in {\mathbb{R}}^{NxT}$   Output: M samples from $p\left({\mathit{a}}_{i,t}\right|{\mathit{x}}_{1:t})$ for the ith GPS station at $t=1\cdots T$ for all $i=1\cdots N$ 1:. $M\leftarrow$ number of particles 2:. for$i=1$ to N do 3:.     Sample ${\mathit{a}}_{i,0}^{\left(m\right)}\sim Prior\left(i\right)$ for $m=1\cdots M$ 4:.     Set weight $w_{i,0}^{\left(m\right)}\leftarrow 1/M$ for $m=1\cdots M$ 5:. end for 6:. for $i=1$ to N do 7:.     for $t=1$ to T do 8:.         Sample ${\mathit{a}}_{i,t}^{\left(m\right)}\sim q\left({\mathit{a}}_{i,t}\right|{\mathit{a}}_{i,t-1}^{\left(m\right)},{\mathit{x}}_t)$ for $m=1\cdots M$ (Prediction step) 9:.         ${\widehat{x}}_{i,t}^{\left(m\right)}\leftarrow$ Equation (1) using ${\mathit{a}}_{i,t}^{\left(m\right)}$ and ${\mathit{x}}_{t-1}$ for $m=1\cdots M$ (Prediction step) 10:.         $w_{i,t}^{\left(m\right)}\leftarrow$ Equation (14) using ${\widehat{x}}_{i,t}^{\left(m\right)}$ and ${\mathit{x}}_t$ for $m=1\cdots M$ (Update step) 11:.         ${\{{\mathit{a}}_{i,t}^{\left(m\right)},w_{i,t}^{\left(m\right)}\}}_{m=1\cdots M}\leftarrow$ Equation (15) with resampling if needed. 12:.     end for 13:. end for

5. Results and Discussion

5.1. Selected Network Data

Two networks were selected for modeling and discussion. The first network, Network 1, is a clustered network of 10 GPS stations near the earthquake epicenter. The second network, Network 2, is a sparse network of 10 GPS stations. Locations of GPS stations in both networks are shown in Figure 1.

Figure 2 shows snapshots of post-processed measurements of the north displacement in meters, retrieved from [13], from GPS stations in the two selected networks. Vertical lines mark the earthquake’s origin time, the time where the earthquake originates at its source, at 14:46:18 on 11 March 2011 (Japan local time) [33]. It can be noted that stations whose locations are near the earthquake epicenter experienced the shaking first; hence, significant displacements were observed at an earlier time.

Time-series were selected from GPS measurements of the north component at each station in both networks from 09:00:05 to 15:25:25 (23,131 data points) from the original 1-Hz data retrieved [13]. It can be noted that vertical displacements were not used in this study because their accuracy is usually less than that of the horizontal ones [34,35].

The number of particles used is $M=$ 10,000. The state noise terms $v_{ij,t}$ in the state update Equation (3) are modeled as i.i.d. Gaussian distributions with zero mean, and a standard deviation of ${10}^{-2}$ meter (1 cm).

For observation noise, we utilize three different distributions:

First, for the Gaussian observation noise, the observation noise terms ${\eta }_{i,t}$ in Equation (1) are assumed to have a zero mean, and a standard deviation of ${\sigma }_{\eta }={10}^{-2}$ m (1 cm). This value was chosen since it was reported that large coseismic ground displacement could be detected by a real-time GPS network (RTK mode) once the displacement exceeds approximately this threshold (1 cm), which represents the GPS data noise level [6].

Second, for the Laplace observation noise, the mean of the observation noise terms ${\eta }_{i,t}$ in Equation (1) is similarly set to zero, and the standard deviation is ${10}^{-2}$ m (1 cm). Consequently, this setting results in a $\beta$ value of $\frac{{\sigma }_{\eta }}{\sqrt{2}}$ for the likelihood Equation (8).

Third, for the Cauchy observation noise, the scale parameter $\gamma$ in the likelihood Equation (9) is chosen to be ${10}^{-2}$ m (1 cm).

The particles were initialized for the first iteration, which represents the prior information of the matrix ${\mathit{A}}_t$ which is a diagonal matrix with unit values along its diagonal, with zeros elsewhere, added with a normal perturbation with a zero mean and a standard deviation of ${10}^{-2}$ m (1 cm) to enhance the variability in the initial state estimates.

5.2. Modeling Results

At each time epoch t, the estimated parameters are present in the form of particles. These particles serve as the basis for deriving valuable statistics, including the mean and standard deviation. Additionally, the distribution can be visually examined through representations such as histograms, which provide richer representations of the entire probability distribution. Figure 3 shows histograms of particles for the hidden parameters in the north components of Network 1 and Network 2 at one instant, under all three distributions for the observation noise. The mean value of particles for a given $a_{ij,t}$ in the matrix ${\mathit{A}}_t$ is considered to be the estimation for the hidden parameter $a_{ij,t}$.

The prediction residuals are calculated as the difference between the observed and estimated values. From the start of the estimation results at 9:00:05 until the earthquake’s origin time at 14:46:18, the mean values of the residuals of both networks are at zero, and the standard deviations of the residuals are 0.22 cm (Gaussian assumption) and 0.21 cm (Laplace and Cauchy assumptions) for Network 1, and 0.19 cm (Gaussian assumption) and 0.18 cm (Laplace and Cauchy assumptions) for Network 2, respectively. These figures indicate a good performance of the PF estimates, with a marginally improved accuracy observed under Laplace and Cauchy noise assumptions.

Figure 4 shows snapshots of the prediction residuals, in centimeters, of GPS stations in Network 1 (top) and Network 2 (bottom), each under three different distribution assumptions for the observation noise. The first vertical lines in all residual graphs mark the earthquake’s origin time. Notably, following the earthquake, there is a marked fluctuation in the estimation performance, indicative of strong disturbances. This sudden fluctuation in the prediction residuals is interpreted as an anomalous event, signifying a deviation from an expected estimation. For example, in the north component of Network 1, an anomaly is detected at 14:47:10 (52 seconds after the origin time) for both Gaussian and Laplace assumptions, and at 14:47:09 for Cauchy assumption. In Network 2, an anomaly is detected at 14:47:33 (1 min and 15 seconds after the origin time) for both Gaussian and Laplace assumptions, and at 14:47:42 for Cauchy assumption. The criteria for the anomaly detection are discussed in Section 5.3.

After the anomalies, up to the 6 min mark after the origin time (as indicated by the third vertical line of each residual graph in Figure 4), Network 1 exhibited a wider range of prediction residuals (−29.1 to 12.8 cm under Gaussian assumption, −11.23 to 15.04 cm under Laplace assumption, and −24.93 to 25.81 cm under Cauchy assumption) compared to Network 2’s range (−11.9 to 4.72 cm under Gaussian assumption, −14.50 to 7.67 cm under Laplace assumption, and −18.44 to 11.58 cm under Cauchy assumption). Note that the largest observed displacements during this period were 1.93 and 0.99 m for Network 1 and Network 2, respectively. The predictions improved around the 6 min post-origin time, marked by the third vertical line in each residual graph, as expected since the disturbance caused by the earthquake started to lessen.

It can be noted that, as the number of GPS stations, N, in a network increases, the number of hidden parameters in the matrix ${\mathit{A}}_t$ will increase quadratically. Consequently, the selection of a subset of GPS stations becomes crucial and poses a challenging task.

5.3. Parameter Choices for Anomaly Detection

In Figure 4, anomalies are identified using specific criteria: The analysis begins by defining a leading period (l) of 260 seconds immediately preceding the earthquake’s origin time. Anomalies are those instances where the the absolute value of the prediction residual surpasses a specified anomaly threshold (z), multiplied by the standard deviation from the mean residual value calculated within the leading period. This condition must be met for a predetermined number of consecutive counts (n). The anomaly threshold z essentially represents the number of standard deviations by which a residual at a given time point deviates from the mean residual value of the leading period. The use of consecutive counts n helps to account for site-specific GPS errors. Specifically, in Figure 4, the anomalies are marked using the following parameters for both networks: l = 260 s window before the origin time, $z=3.0$, and $n=5$.

With these parameters, an anomaly in the east component of Network 1 is detected slightly earlier, at 14:47:05. This is consistent with expectations, as the seismic displacement from the Tohoku-Oki earthquake was first detected in the east component [6].

Our approach for identifying anomalies, which utilizes prediction residuals derived from the estimation of hidden parameters, can offer more sensitivity to the anomaly detection than if solely relying on raw displacement data. For example, when comparing anomaly detection based on prediction residuals with a straightforward approach that applies a short-term moving variance directly to the displacement observations, using the parameters previously discussed for the north component of Network 2, our method detects the anomaly the earliest at 14:47:33 (as shown in Figure 4). In contrast, with the same parameter settings, direct thresholding of the displacement observations identifies the anomaly later, at 14:47:43.

It is essential to note that the selected anomaly threshold z, the leading period l, and the consecutive counts n are arbitrary and can be adjusted for each network to accommodate the varying behaviors in estimation noise. For instance, in the north component of both Network 1 and 2, changing the leading period to either a 360 s or 460 s window before the earthquake’s origin time still results in detecting an anomaly within 1 s of the initial detection time, under all distribution assumptions, for the same anomaly threshold.

5.4. Graph Representation

The hidden parameters $a_{ij,t}$ represent relationships between GPS displacements in a network at different time epochs t. Instead of comparing raw values, we assess z-scores, denoted $z_{ij,t}$. A z-score is a statistical measure that indicates the number of standard deviations a data point is from the mean of its distribution. Specifically, at each time epoch t, a z-score of an estimated hidden parameter $a_{ij,t}$ denoted $z_{ij,t}$ is calculated as $z_{ij,t}=\frac{(a_{ij,t}-{\mu }_t)}{{\sigma }_t}$, where ${\mu }_t$ and ${\sigma }_t$ are the mean and the standard deviation of all particles $a_{ij,t}^{\left(m\right)}$, in a matrix ${\mathit{A}}_t$ of Equation (2).

A graphical representation of a network at time epoch t consists of N nodes, each representing a GPS station, and $NxN$ directed edges. Directed edges from j to i represent $z_{ij,t}$. A high positive $z_{ij,t}$ reflects a positive influence of measurement of GPS station j on i, while near-zero $z_{ij,t}$ indicates proximity to the mean. Negative $z_{ij,t}$ signifies a negative influence of the measurement of GPS station j on i.

Figure 5 depicts graph representations of the hidden parameters for the north displacements of Network 1 and Network 2 under the three different assumptions for the observation noise, each at three distinct time instants. The graphs for the first two time instants (the origin time and the anomalies) of each network exhibit similar behaviors, noticeable by the similar colors of their edges. However, at the third time instant shown, both networks display graphs distinct from those in the previous two time instants. This divergence is understandable since, at the marked anomalies of both networks (center column), the networks began to undergo changes, reflected by high estimate residuals (noticeable in Figure 4), caused by the disturbance from the earthquake that the model did not capture well initially. Subsequently, at the third time instant (the rightmost column), the model learns and presents different behaviors, as indicated by graphs with distinctly different edge colors. It can be noted that in Network 2, at the third time instant shown, each GPS station’s measurement heavily relies solely on its own previous measurements. This may be due to the highly sparse locations of the stations and hence less dependency after the disturbance from the earthquake. Additionally, it can be noted that the Cauchy observation noise assumptions lead to sparser networks with less significant branches, as shown by the lighter colors of the graph edges.

The proposed graph representation of a network enables the tracking of hidden parameters in a compact manner, with edges defined using the z-scores of the estimated parameters. This representation should also aid the understanding of networks following weaker earthquakes, whose hidden parameters are expected to fluctuate more subtly. Other potential graph representations include a graph whose edges are defined based on the amount of change in the parameters relative to those in the previous time step.

6. Conclusions

We showcased the capabilities of a sequential Monte Carlo method, specifically particle filtering, in tracking changes in the relationships among GPS measurements in a network before, during, and after the disturbance caused by an earthquake. Our proposed model effectively captures the time-varying behaviors of the network, which can be useful for anomalous displacement detection. The proposed graph representations aid in understanding and facilitate the tracking of network dynamics. The versatility of the method allows for extensions to model other time-varying geodetic data networks, and it facilitates the adoption of different model equations and assumptions. In the future, regarding this work, we will also consider other distribution models [36] and nonlinear state and observation models.

Author Contributions

Conceptualization, E.E.K.; Methodology, E.E.K.; Formal analysis, S.P.; Writing—original draft, S.P.; Writing—review & editing, S.P. and E.E.K.; Visualization, S.P.; Supervision, E.E.K. and M.S.O.; Funding acquisition, E.E.K. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original code presented in the study is openly available at https://github.com/SuchanunP/pf_gps_dynamics/ (accessed on 2 April 2024). The animated graph representations of Network 1 under three distribution assumptions for the observation noise are included in the Supplementary Materials, and further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/e26040342/s1.

References

1. Segall, P. Earthquake and Volcano Deformation; Princeton University Press: Princeton, NJ, USA, 2010.

2. Segall, P.; Davis, J.L. GPS Applications for Geodynamics and Earthquakes Studies. Annu. Rev. Earth Planet. Sci.; 1997; 25, pp. 301-336. [DOI: https://dx.doi.org/10.1146/annurev.earth.25.1.301]

3. Segall, P.; Matthews, M. Time dependent inversion of geodetic data. J. Geophys. Res. Solid Earth; 1997; 102, pp. 22391-22409. [DOI: https://dx.doi.org/10.1029/97JB01795]

4. Allen, R.M.; Ziv, A. Application of real-time GPS to earthquake early warning. Geophys. Res. Lett.; 2011; 38, [DOI: https://dx.doi.org/10.1029/2011GL047947]

5. Bock, Y.; Prawirodirdjo, L.; Melbourne, T.I. Detection of arbitrarily large dynamic ground motions with a dense high-rate GPS network. Geophys. Res. Lett.; 2004; 31, [DOI: https://dx.doi.org/10.1029/2003GL019150]

6. Psimoulis, P.A.; Houlié, N.; Habboub, M.; Michel, C.; Rothacher, M. Detection of ground motions using high-rate GPS time-series. Geophys. J. Int.; 2018; 214, pp. 1237-1251. [DOI: https://dx.doi.org/10.1093/gji/ggy198]

7. Tsai, M.C.; Yu, S.B.; Shin, T.C.; Kuo, K.W.; Leu, P.L.; Chang, C.H.; Ho, M.Y. Velocity Field Derived from Taiwan Continuous GPS Array (2007–2013). Terr. Atmos. Ocean. Sci.; 2015; 26, 527. [DOI: https://dx.doi.org/10.3319/TAO.2015.05.21.01(T)]

8. Zubovich, A.V.; Wang, X.Q.; Scherba, Y.G.; Schelochkov, G.G.; Reilinger, R.; Reigber, C.; Mosienko, O.I.; Molnar, P.; Michajljow, W.; Makarov, V.I. et al. GPS velocity field for the Tien Shan and surrounding regions. Tectonics; 2010; 29, [DOI: https://dx.doi.org/10.1029/2010TC002772]

9. Bitharis, S.; Pikridas, C.; Fotiou, A.; Rossikopoulos, D. GPS data analysis and geodetic velocity field investigation in Greece, 2001–2016. GPS Solut.; 2023; 28, 16. [DOI: https://dx.doi.org/10.1007/s10291-023-01549-8]

10. Sobolev, G.A.; Zakrzhevskaya, N.A.; Akatova, K.N.; Gitis, V.G.; Derendyaev, A.B.; Bragin, V.D.; Sycheva, N.A.; Kuzikov, S.I. Dynamics of interaction between fields of seismicity and surface deformations (Bishkek geodynamic test area). Izv. Phys. Solid Earth; 2010; 46, pp. 817-838. [DOI: https://dx.doi.org/10.1134/S1069351310100034]

11. Jiang, Y.; Luo, H.; Xu, Q.; Lu, Z.; Liao, L.; Li, H.; Hao, L. A Graph Convolutional Incorporating GRU Network for Landslide Displacement Forecasting Based on Spatiotemporal Analysis of GNSS Observations. Remote Sens.; 2022; 14, 1016. [DOI: https://dx.doi.org/10.3390/rs14041016]

12. Gitis, V.; Derendyaev, A. A Technology for Seismogenic Process Monitoring and Systematic Earthquake Forecasting. Remote Sens.; 2023; 15, 2171. [DOI: https://dx.doi.org/10.3390/rs15082171]

13. Psimoulis, P. 1 Hz GPS Japan Network (GEONET) Data for Tohoku-Oki Mw9.0 2011 Earthquake; The University of Nottingham: Nottingham, UK, 2019; [DOI: https://dx.doi.org/10.17639/nott.6999]

14. Sagiya, T. A decade of GEONET: 1994–2003. Earth Planets Space; 2004; 56, pp. xxix-xli. [DOI: https://dx.doi.org/10.1186/BF03353077]

15. The Geospatial Information Authority of Japan (GSI). 2011 Tohoku Earthquake: Dataset; GSI: Tsukuba, Japan, 2011.

16. Mao, A.; Harrison, C.; Dixon, T. Noise in GPS coordinate time series. J. Geophys. Res.; 1999; 104, pp. 2797-2816. [DOI: https://dx.doi.org/10.1029/1998JB900033]

17. Langbein, J. Noise in GPS displacement measurements from Southern California and Southern Nevada. J. Geophys. Res.; 2008; 113, [DOI: https://dx.doi.org/10.1029/2007JB005247]

18. Klos, A.; Bogusz, J.; Figurski, M.; Kosek, W. Irregular variations in GPS time series by probability and noise analysis. Surv. Rev.; 2015; 47, pp. 163-173. [DOI: https://dx.doi.org/10.1179/1752270614Y.0000000133]

19. Kotecha, J.; Djuric, P. Gaussian sum particle filtering. IEEE Trans. Signal Process.; 2003; 51, pp. 2602-2612. [DOI: https://dx.doi.org/10.1109/TSP.2003.816754]

20. Misra, N.; Kuruoglu, E.E. Stable Graphical Models. J. Mach. Learn. Res.; 2016; 17, pp. 1-36.

21. Kantas, N.; Doucet, A.; Singh, S.S.; Maciejowski, J.; Chopin, N. On Particle Methods for Parameter Estimation in State-Space Models. Stat. Sci.; 2015; 30, pp. 328-351. [DOI: https://dx.doi.org/10.1214/14-STS511]

22. Ancherbak, S.; Kuruoglu, E.E.; Vingron, M. Time-Dependent Gene Network Modelling by Sequential Monte Carlo. IEEE/ACM Trans. Comput. Biol. Bioinform.; 2016; 13, pp. 1183-1193. [DOI: https://dx.doi.org/10.1109/TCBB.2015.2496301]

23. Chen, J.; Hubbard, S.S.; Williams, K.H.; Pride, S.; Li, L.; Steefel, C.; Slater, L. A state-space Bayesian framework for estimating biogeochemical transformations using time-lapse geophysical data. Water Resour. Res.; 2009; 45, [DOI: https://dx.doi.org/10.1029/2008WR007698]

24. Ambrosi, P.; Costagli, M.; Kuruoğlu, E.E.; Biagi, L.; Buonincontri, G.; Tosetti, M. Modeling brain connectivity dynamics in functional magnetic resonance imaging via particle filtering. Brain Inform.; 2021; 8, 19. [DOI: https://dx.doi.org/10.1186/s40708-021-00140-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34586519]

25. Miao, Y.; Kowal, D.R.; Panchal, N.; Vila, J.; Vannucci, M. Nonlinear state-space modeling approaches to real-time autonomous geosteering. J. Pet. Sci. Eng.; 2020; 189, 107025. [DOI: https://dx.doi.org/10.1016/j.petrol.2020.107025]

26. Arulampalam, M.; Maskell, S.; Gordon, N.; Clapp, T. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Trans. Signal Process.; 2002; 50, pp. 174-188. [DOI: https://dx.doi.org/10.1109/78.978374]

27. Kalman, R.E. A New Approach to Linear Filtering and Prediction Problems. J. Basic Eng.; 1960; 82, pp. 35-45. [DOI: https://dx.doi.org/10.1115/1.3662552]

28. Sura, P.; Gille, S. Interpreting wind-driven Southern Ocean variability in a stochastic framework. J. Mar. Res.; 2003; 61, pp. 313-334. [DOI: https://dx.doi.org/10.1357/002224003322201214]

29. Peinke, J.; Böttcher, F.; Barth, S. Anomalous statistics in turbulence, financial markets and other complex systems. Ann. Der Phys.; 2004; 13, pp. 450-460. [DOI: https://dx.doi.org/10.1002/andp.200451607-809]

30. Doucet, A.; Godsill, S.; Andrieu, C. On sequential Monte Carlo sampling methods for Bayesian filtering. Stat. Comput.; 2000; 10, pp. 197-208. [DOI: https://dx.doi.org/10.1023/A:1008935410038]

31. Candy, J. Bootstrap Particle Filtering. Signal Process. Mag. IEEE; 2007; 24, pp. 73-85. [DOI: https://dx.doi.org/10.1109/MSP.2007.4286566]

32. Liu, J.S.; Chen, R. Sequential Monte Carlo Methods for Dynamic Systems. J. Am. Stat. Assoc.; 1998; 93, pp. 1032-1044. [DOI: https://dx.doi.org/10.1080/01621459.1998.10473765]

33. Fujinawa, Y.; Noda, Y. Japan’s Earthquake Early Warning System on 11 March 2011: Performance, Shortcomings, and Changes. Earthq. Spectra; 2013; 29, pp. 341-368. [DOI: https://dx.doi.org/10.1193/1.4000127]

34. Murray-Moraleda, J. GPS: Applications in Crustal Deformation Monitoring. Encyclopedia of Complexity and Systems Science; Meyers, R.A. Springer: New York, NY, USA, 2009; pp. 4249-4283.

35. Yuan, L.; Chao, B.F.; Ding, X.; Zhong, P. The tidal displacement field at Earth’s surface determined using global GPS observations. J. Geophys. Res. (Solid Earth); 2013; 118, pp. 2618-2632. [DOI: https://dx.doi.org/10.1002/jgrb.50159]

36. Kuruoglu, E. Analytical representation for positive alpha-stable densities. Proceedings of the 2003 IEEE International Conference on Acoustics, Speech, and Signal Processing; Hong Kong, China, 6–10 April 2003; Volume 6, VI-729. [DOI: https://dx.doi.org/10.1109/ICASSP.2003.1201785]