the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Data assimilation with an improved particle filter and its application in the TRIGRS landslide model
Changhu Xue
Guigen Nie
Haiyang Li
Jing Wang
Particle filters have become a popular algorithm in data assimilation for their ability to handle nonlinear or nonGaussian statespace models, but they have significant disadvantages. In this work, an improved particle filter algorithm is proposed. To overcome the particle degeneration and improve particles' efficiency, the processes of particle resampling and particle transfer are updated. In this improved algorithm, particle propagation and the resampling method are ameliorated. The new particle filter is applied to the Lorenz63 model, and its feasibility and effectiveness are verified using only 20 particles. The rootmeansquare difference (RMSD) of estimations converges to stable when there are more than 20 particles. Finally, we choose a peristaltic landslide model and carry out an assimilation experiment of 20 days. Results show that the estimations of states can effectively correct the running offset of the model and the RMSD is convergent after 3 days of assimilation.
 Article
(996 KB) 
Supplement
(305 KB)  BibTeX
 EndNote
Mountainous areas all over the world suffer frequent landslide disasters. Works of landslide monitoring, analysis and forecasting are crucial. Many numerical modeling methods of slope evolution, such as discontinuous deformation analysis (DDA) (Shi, 1992; Jing et al., 2001; Ma et al., 2011) and the distinct element methods (DEMs) (Lorig and Hobbs, 1990; Marcato et al., 2007; Li et al., 2012), have been proposed and developed recently. Iverson proposed a mathematical model that uses Richards' equation to evaluate effects of landslides in response to rainfall infiltration (Iverson, 2000). The Transient Rainfall Infiltration and GridBased Regional SlopeStability (TRIGRS) model is a rasterbased model and depends on time of transient rainfall infiltration (Baum et al., 2008). Jiang adopted the ensemble Kalman filter landslide movement model in relation to hydrological factors, which introduced data assimilation (DA) to landslide evaluation (Jiang et al., 2016).
DA is a common approach to estimating optimal states in dynamic systems. With DA algorithms and operators, DA merges different scales of observations into dynamic models to take advantage of all the information. Many DA algorithms have been developed and improved in recent years, and particle filters (PFs) are a popular algorithm for their ability to handle nonlinear and nonGaussian distributed models (Arulampalam et al., 2002; Moradkhani et al., 2005). The application and improvement of PFs has been researched recently in DA and other fields.
Salamon and Feyen (2009) applied the residual resampling particle filter (RRPF) to assess parameter, precipitation and predictive uncertainty in a rainfall–runoff model. Thirel et al. (2013) assimilated snowcovered areas in physical distributed hydrological models and MODIS satellite data to improve panEuropean flood forecasts. Mattern et al. (2013) carried out assimilation experiments for a threedimensional biological ocean model and satellite observations and verified the feasibility of biological state estimation with sequential importance resampling (SIR) for realistic models.
However, large computational complexity and particle degradation or collapse are still disadvantages of PFs. To solve these problems, some resampling algorithms have been proposed. One improvement is adding an item related to observations to make the proposal density dependent on future observations; accordingly most particles can situate into the range of observation error (van Leeuwen, 2010). This method can achieve good results using only 10–20 particles in highdimensional assimilation experiments. But the number of key particles is reduced when the system variance is larger than the observed variance, and the values of added items are uncertain. Another improvement is to replace the duplicating process by generating a Halton sequence in residual resampling (Zhang et al., 2013). The disordered particle sets are turned into ordered sets and too few particles can hardly describe the posterior probability density function (PDF) better.
In Sect. 2, a new resampling approach is proposed to improve the above method, maintaining both particle diversity and efficiency. The new algorithm formula and implementation process are listed in the text. To predict the safety factor of peristaltic landslide, a simulation experiment, applied to the Lorenz63 model using different numbers of particles, ranging between 10 and 200, is explained in Sect. 3, which demonstrates that the new method shows efficiency and sensitivity to the number of particles. Finally, a rainfall infiltration landslide model case is analyzed. We choose an experimental landslide model with a 10×10 size grid, of which the side length of each grid cell is 10 m. The improved assimilation algorithm is applied to the TRIGRS program to evaluate the change of factor of safety (FS) in the experimental model.
In sequential importance sampling, the state vector is represented by a set of particles:
where x is the state vector with initial PDF p(x_{0}), k is the subscript of time steps, ε_{k−1} is system noise with zero mean at step k−1 and f(⋅) is the model operator. Initial N particles are sampled from p(x_{0}). The observation equation is
where z_{k} is the observation vector at step k, and h(⋅) is the observation operator. Weights of particles are calculated by Eq. (3) and normalized to obtain ${w}_{k}^{i}$ by Eq. (4)
where i is the index of particle number, $p\left({\mathit{z}}_{k}\mathrm{}{\mathit{x}}_{k}^{i}\right)$ is the likelihood of observation and $q({\mathit{x}}_{k}^{i}\mathrm{}{\mathit{x}}_{k\mathrm{1}}^{i},{\mathit{z}}_{k})$ is the proposal function.
Residual resampling is a way to solve particle degeneracy, which is an unavoidable problem in PFs. To keep most particles effective, lowweight particles are removed and highweight particles are duplicated. But with recursive progress the particle sets can hardly represent the prior PDF due to declining particle diversity.
Some improvements to the residual resampling algorithm are proposed in this paper. Firstly, in the process of particle transferring, we choose
where J_{k} is a coefficient like the “gain” in an extended Kalman filter:
in which A_{k}, B_{k} are the linearization parameters of f(⋅) and h(⋅), respectively:
D_{k∕k} is the estimation variance of state x_{k} at step k. This process is equal to translating particles close to observations. But the value of J_{k} is difficult to determine because the variance of state estimation ${\mathbf{D}}_{k\mathrm{1}/k\mathrm{1}}$ in PFs is difficult to compute. To simplify the calculation, suppose that the translated particles are a series of virtual observations about the state at step k. Write the particle set as
and replace ${\mathbf{D}}_{k\mathrm{1}/k\mathrm{1}}$ with the variance of particles. To keep the value of ${\mathbf{D}}_{k\mathrm{1}/k\mathrm{1}}$ unchanged before and after translation, we choose the posterior particles at step k−1:
Secondly, using the method of Zhang et al. (2013) to compute accumulative copy times (ACTs), each parent particle with high weight regenerates a set of new particles. Differently, instead of duplicating or generating a Halton sequence, it generates a series of normally distributed particles:
where ACT_{i} is the ACT of the ith particle, and the mean and variance are related to the value of the parent. Accordingly, the resampled particle set is composed of some different particle sets that obey normal distribution. Assuming that the jth particle of ${\mathit{x}}_{k}^{i}$ is written as ${\mathit{x}}_{k}^{ij}$, the formula (3) can be written as
Briefly, the improved RRPF in this section can be implemented by the following steps:

Step 1. Draw initial particles $\mathit{\left\{}{\mathit{x}}_{\mathrm{0}}^{i}\mathit{\right\}}$ from the prior PDF p(x_{0}).

Step 2. Compute the mean and variance of posterior particles at step k1.
$$\begin{array}{}\text{(11)}& {\stackrel{\mathrm{\u203e}}{\mathit{x}}}_{k\mathrm{1}/k\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^{N}{\mathit{x}}_{k\mathrm{1}/k\mathrm{1}}^{i}\end{array}$$$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{\mathbf{D}}_{k\mathrm{1}/k\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{N\mathrm{1}}}\sum _{i=\mathrm{1}}^{N}\left({\mathit{x}}_{k\mathrm{1}/k\mathrm{1}}^{i}{\stackrel{\mathrm{\u203e}}{\mathit{x}}}_{k\mathrm{1}/k\mathrm{1}}\right)\\ \text{(12)}& {\displaystyle}& {\displaystyle}{\left({\mathit{x}}_{k\mathrm{1}/k\mathrm{1}}^{i}{\stackrel{\mathrm{\u203e}}{\mathit{x}}}_{k/k\mathrm{1}}\right)}^{T}\end{array}$$ 
Step 3. Using the new method in this section, compute the gains of particles.
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{\mathbf{D}}_{k/k\mathrm{1}}=\left[{\displaystyle \frac{\partial {f}_{k}}{\partial {\mathit{x}}_{k}}}\left({\widehat{\mathit{x}}}_{k}\right)\right]{\mathbf{D}}_{k\mathrm{1}/k\mathrm{1}}{\left[{\displaystyle \frac{\partial {f}_{k}}{\partial {\mathit{x}}_{k}}}\left({\widehat{\mathit{x}}}_{k}\right)\right]}^{T}\\ \text{(13)}& {\displaystyle}& {\displaystyle}+{\mathbf{G}}_{k\mathrm{1}}\left({\widehat{\mathit{x}}}_{k\mathrm{1}}\right){\mathbf{Q}}_{k\mathrm{1}}{\mathbf{G}}_{k\mathrm{1}}^{T}\left({\widehat{\mathit{x}}}_{k\mathrm{1}}\right)\end{array}$$$$\begin{array}{ll}{\displaystyle}{\mathbf{J}}_{k}=& {\displaystyle}{\mathbf{D}}_{k/k\mathrm{1}}\left[{\displaystyle \frac{\partial {f}_{k}}{\partial {\mathit{x}}_{k}}}\left({\widehat{\mathit{x}}}_{k/k\mathrm{1}}\right)\right]\left(\right)open="\{">\left[{\displaystyle \frac{\partial {f}_{k}}{\partial {\mathit{x}}_{k}}}\left({\widehat{\mathit{x}}}_{k/k\mathrm{1}}\right)\right]\end{array}\text{(14)}& {\displaystyle}& {\displaystyle}{\left({\mathbf{D}}_{k/k\mathrm{1}}{\left[{\displaystyle \frac{\partial {f}_{k}}{\partial {\mathit{x}}_{k}}}\left({\widehat{\mathit{x}}}_{k/k\mathrm{1}}\right)\right]}^{T}+{\mathbf{R}}_{k}\right\}}^{\mathrm{1}}$$ 
Step 4. Transfer the particles close to the observation:
$$\begin{array}{}\text{(15)}& {\mathit{x}}_{k}^{i}=f\left({\mathit{x}}_{k\mathrm{1}}^{i}\right)+{\widehat{\mathit{\epsilon}}}_{k\mathrm{1}}+{\mathbf{J}}_{k}[{\mathit{z}}_{k}h({\widehat{\mathit{x}}}_{k\mathrm{1}}\left)\right].\end{array}$$ 
Step 5. In residual resampling, each particle generates a set of normally distributed progeny particles, and all progeny sets make
up the resampled particle set:
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}\left\{{\mathit{x}}_{k}^{i\mathrm{1}},{\mathit{x}}_{k}^{i\mathrm{2}},\mathrm{\dots},{\mathit{x}}_{k}^{i{\text{ACT}}_{i}}\right\}={\mathit{X}}_{k}^{i{\text{ACT}}_{i}}\\ \text{(16)}& {\displaystyle}& {\displaystyle}\sim N\left({\mathit{x}}_{k}^{i},{\mathbf{G}}_{k}\left({\mathit{x}}_{k}^{i}\right)\right),\end{array}$$$$\begin{array}{}\text{(17)}& {\displaystyle}\left\{{\mathit{X}}_{k}^{\mathrm{1}{\text{ACT}}_{\mathrm{1}}},{\mathit{X}}_{k}^{\mathrm{2}{\text{ACT}}_{\mathrm{2}}},\mathrm{\dots},{\mathit{X}}_{k}^{N\phantom{\rule{0.25em}{0ex}}{\text{ACT}}_{N}}\right\}={\left\{{\mathit{x}}_{k}^{\ast i}\right\}}_{i=\mathrm{1},\mathrm{2},\mathrm{\dots},N},\end{array}$$when ACT_{i}=0, ${\mathit{X}}_{k}^{i{\text{ACT}}_{i}}$ is an empty set.

Step 6. Compute and normalize weights.
$$\begin{array}{}\text{(18)}& {\stackrel{\mathrm{\u0303}}{w}}_{k}^{i}={w}_{k\mathrm{1}}^{i}\cdot p\left({\mathit{z}}_{k}\mathrm{}{\mathit{x}}_{k}^{i}\right)\end{array}$$$$\begin{array}{}\text{(19)}& {w}_{k}^{i}={\displaystyle \frac{{\stackrel{\mathrm{\u0303}}{w}}_{k}^{i}}{\sum _{j=\mathrm{1}}^{N}{\stackrel{\mathrm{\u0303}}{w}}_{k}^{j}}}\end{array}$$ 
Step 7. Compute the state estimation.
$$\begin{array}{}\text{(20)}& {\widehat{\mathit{x}}}_{k/k}=\sum _{i=\mathrm{1}}^{N}{\mathit{x}}_{k}^{\ast i}\cdot {w}_{k}^{i},\end{array}$$a measure to assess the accuracy of calculation is the root mean square difference (RMSD), which is defined as
$$\begin{array}{}\text{(21)}& \text{RMSD}=\sqrt{{\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^{T}({\widehat{\mathit{X}}}_{t}{\mathit{X}}_{t}^{\text{obs}}{)}^{\mathrm{2}}},\end{array}$$where T is the period of assimilation, and ${\widehat{\mathit{X}}}_{t}$ and ${\mathit{X}}_{t}^{\text{obs}}$ are the assimilated value and the observation of state at time t, respectively.
We choose the Lorenz63 model as an example to test the improved algorithm (Baines, 2008).
where the constants σ, ρ and β are system parameters proportional to the Prandtl number, Rayleigh number and certain physical dimensions of the layer itself. Parameters are given by dt=0.01, σ=10, ρ=28, $\mathit{\beta}=\mathrm{8}/\mathrm{3}$, the observation error ${\mathit{\sigma}}_{\text{obs}}=\sqrt{\mathrm{2}}$ and the model transmission error based on time interval ${\mathit{\sigma}}_{\text{mod}}=\mathrm{2}\sqrt{\mathrm{\Delta}t}$. Initialize the filter with the starting point, which is set to (x_{0}, y_{0}, ${z}_{\mathrm{0}})=(\mathrm{1.50887},\mathrm{1.531271},\mathrm{25.46091})$. The truth is obtained by the formula of the model recursively. Observations are generated from the truth by adding a disturbance every 40 steps, with 1000 recurring steps, and assimilating the observation with the model when observation exists at the current step and moving to the next step when there is no observation.
Figure 1 shows the results of the x component using the new PF with 20 particles. Note that the new PF procedure is close to the truth with much fewer particles, which is more efficient than the standard PF procedure with hundreds of particles. Compute the confidence interval with the 95 % level using the posterior particles every step. Figure 2 shows that the intervals contain observations at almost all the steps at which observations exist. That means particle sets after translation are closer to observations and true states. The evolution of all particles is displayed in Fig. 3, in which most particles are very close to observations except for several ones at moments when the state changed obviously. The RMSD sequence is shown in Fig. 4; it tends to be stable when the number of particles is more than 20. This means the improved algorithm only needs no fewer than 20 particles.
TRIGRS is a program modeling rainfall infiltration, using analytical solutions for partial differential equations that represent onedimensional vertical flow in isotropic, homogeneous materials for simply saturated or unsaturated conditions. It computes changes of rainfall pore pressure and FS with rainfall infiltration. The FS is computed using a simple infiniteslope model, cell by cell.
In this experiment, the FS is applied to assimilation. It is calculated as follows:
in which c is soil cohesion, α is slope angle, ϕ is soil friction angle, φ is the groundwater pressure head depending on depth Z and time t, γ_{W} is groundwater unit weight, and γ_{S} is soil unit weight at saturation.
An example of the 10×10 grid TRIGRS model is set to be the background, and each grid is a square with a length of 10 m. The simulated observations are generated from the FS by adding a disturbance with normal distribution N(0.2, 0.3). Due to the difficulties of determining the parameter φ, the soil friction angle and its high sensitivity to results, we now generate a set of particles $\left\{{\mathit{\phi}}_{k}^{i}\right\}$ to form φ, in which k and i are indices of step and particle number, respectively. The input model variance of φ is 2 and observation variance of FS is 0.3. At each step, φ and FS will be updated, and the updated parameters continue to participate in the next step of operation as initial parameters. The number of particles is set to 20 in the PF program. Figure 5 shows the model running results and the assimilated results of the FS running for 5 days, 10 days, 15 days and 20 days. In the model running results, the value of FS is smaller and decreases rapidly, while in the assimilated results the change is relatively gentle.
To evaluate the distribution variation in φ, we propose that the estimation of φ is calculated as formula
in which ${w}_{k}^{i}$ is calculated using Eqs. (18) and (19). Figure 6 shows the distribution variation in φ running for 5 days, 10 days, 15 days and 20 days. Actually, the estimation of φ uses the same method and particles of the estimation of FS. Figure 6 shows the distribution variation in φ running for 5 days, 10 days, 15 days and 20 days. The change of φ estimation in a single cell is illustrated in Fig. 7, considering the middle unit, grid cell (5, 5).
To assess estimations of all grid cells, the RMSD of the whole grid of points to measure the estimated error is modified to
where N_{p} is the total number of grid points, and i and j are the indices of the row and column number, respectively. The RMSD curve with assimilating days is shown in Fig. 8, which suggests the value is large on the first 2 days of initialization, fluctuates in the next days and is steady when there are no observations.
The problems of particle degeneration and efficient expression of posterior PDF are longterm difficulties that affect the performance of particle filters. Many resampling methods can improve effectiveness of particles, but they still need a large number of samples resulting in a large amount of computation.
In this study, we propose two approaches to improve the particle filter process. Firstly, for the problem of particle degeneration, new Gaussiandistributed offspring particles are generated for each mother particle. This avoids particle duplication and maintains particle diversity. Secondly, in order to improve the propagating efficiency of a priori particles into a posteriori particles, an additional item is added that is similar to the Kalman gain at the step of particle propagation, which greatly reduces the number of particles required. It uses only dozens of particles to achieve good results. A simulation experiment of the Lorenz63 model is carried out to validate the feasibility of these methods. The TRIGRS landslide model is first proposed for application to the assimilation system. Results show that the assimilation process can make the estimation close to observations, which proves the feasibility of applying the improved particle filter to the landslide model.
However, some disadvantages are still present. Grid cells are independent of each other in TRIGRS, and this leads to the FS estimations possibly being greater than the actual values. Therefore, the FS estimations only provide a reference for the actual values. The experiment needs improvement.
The data in this experiment are mainly collected by observation or generated by simulation. Data sets and figures have been uploaded to the Supplement. The model information and program of TRGIRS are obtained from https://pubs.usgs.gov/of/2008/1159/ (US Geological Survey, 2017).
The supplement related to this article is available online at: https://doi.org/10.5194/nhess1828012018supplement.
CX conceived the idea of this article and completed the paper. GN proposed some important suggestions of algorithm applications in the TRIGRS assimilation. Validation of the algorithm, chart sorting and editing were completed by HL and JW.
The authors declare that they have no conflict of interest.
This work is financially supported by the National Key
Basic Research Program of China (grant no. 2013CB733205).
Edited by: Thomas Glade
Reviewed by: two anonymous referees
Arulampalam, M. S., Maskell, S., Gordon, N., and Clapp, T.: A tutorial on particle filters for online nonlinear/nonGaussian Bayesian tracking, IEEE Trans. Signal Process., 50, 174–188, https://doi.org/10.1109/78.978374, 2002.
Baines, P. G.: Lorenz, EN 1963: Deterministic nonperiodic flow. Journal of the Atmospheric Sciences 20, 130–41, Prog. Phys. Geogr., 32, 475–480, https://doi.org/10.1177/0309133308091948, 2008.
Baum, R. L., Savage, W. Z., and Godt, J. W.: TRIGRS – A Fortran Program for Transient Rainfall Infiltration and GridBased Regional SlopeStability Analysis, version 2.0, U.S. Geological Survey OpenFile Report, 2008–1159, 75, 2008.
Iverson, R. M.: Landslide triggering by rain infiltration, Water Resour. Res., 36, 1897–1910, https://doi.org/10.1029/2000wr900090, 2000.
Jiang, Y. A., Liao, M. S., Zhou, Z. W., Shi, X. G., Zhang, L., and Balz, T.: Landslide Deformation Analysis by Coupling Deformation Time Series from SAR Data with Hydrological Factors through Data Assimilation, Remote Sens., 8, 179, https://doi.org/10.3390/rs8030179, 2016.
Jing, L. R., Ma, Y., and Fang, Z. L.: Modeling of fluid flow and solid deformation for fractured rocks with discontinuous deformation analysis (DDA) method, Int. J. Rock Mech. Min. Sci., 38, 343–355, https://doi.org/10.1016/S13651609(01)000053, 2001.
Li, X. P., He, S. M., Luo, Y., and Wu, Y.: Simulation of the sliding process of Donghekou landslide triggered by the Wenchuan earthquake using a distinct element method, Environ. Earth Sci., 65, 1049–1054, https://doi.org/10.1007/s1266501109538, 2012.
Lorig, L. J. and Hobbs, B. E.: Numerical Modeling of Slip Instability Using the Distinct Element Method with State Variable Friction Laws, Int. J. Rock Mech. Min. Sci. Geomech. Abstr., 27, 525–534, https://doi.org/10.1016/01489062(90)91003P, 1990.
Ma, G. C., Kaneko, F., Hori, S., and Nemoto, M.: Use of Discontinuous Deformation Analysis to Evaluate the Dynamic Behavior of Submarine TsunamiGenerating Landslides in the Marmara Sea, Int. J. Comput. Methods, 8, 151–170, https://doi.org/10.1142/S0219876211002526, 2011.
Marcato, G., Fujisawa, K., Mantovani, M., Pasuto, A., Silvano, S., Tagliavini, F., and Zabuski, L.: Evaluation of seismic effects on the landslide deposits of Monte Salta (Eastern Italian Alps) using distinct element method, Nat. Hazards Earth Syst. Sci., 7, 695–701, https://doi.org/10.5194/nhess76952007, 2007.
Mattern, J. P., Dowd, M., and Fennel, K.: Particle filterbased data assimilation for a threedimensional biological ocean model and satellite observations, J. Geophys. Res.Oceans, 118, 2746–2760, https://doi.org/10.1002/jgrc.20213, 2013.
Moradkhani, H., Hsu, K. L., Gupt, H. A., and Sorooshian, S.: Uncertainty assessment of hydrologic model states and parameters: Sequential data assimilation using the particle filter, Water Resour. Res., 41, 17, https://doi.org/10.1029/2004wr003604, 2005.
Salamon, P. and Feyen, L.: Assessing parameter, precipitation, and predictive uncertainty in a distributed hydrological model using sequential data assimilation with the particle filter, J. Hydrol., 376, 428–442, https://doi.org/10.1016/j.jhydrol.2009.07.051, 2009.
Shi, G. H.: Discontinuous Deformation Analysis: A New Numerical Model for the Statics and Dynamics of Deformable Block Structures, Eng. Comput., 9, 157–168, https://doi.org/10.1108/eb023855, 1992.
Thirel, G., Salamon, P., Burek, P., and Kalas, M.: Assimilation of MODIS Snow Cover Area Data in a Distributed Hydrological Model Using the Particle Filter, Remote Sens., 5, 5825–5850, https://doi.org/10.3390/rs5115825, 2013.
US Geological Survey: TRIGRS – A Fortran Program for Transient Rainfall Infiltration and GridBased Regional SlopeStability Analysis, Version 2.0, available at: https://pubs.usgs.gov/of/2008/1159/, last access: 6 December 2017.
van Leeuwen, P. J.: Nonlinear data assimilation in geosciences: an extremely efficient particle filter, Q. J. Roy. Meteorol. Soc., 136, 1991–1999, https://doi.org/10.1002/qj.699, 2010.
Zhang, H. J., Qin, S. X., Ma, J. W., and You, H. J.: Using Residual Resampling and Sensitivity Analysis to Improve Particle Filter Data Assimilation Accuracy, IEEE Geosci. Remote Sens. Lett., 10, 1404–1408, https://doi.org/10.1109/Lgrs.2013.2258888, 2013.