Kinematic landslide monitoring with Kalman filtering

. Landslides are serious geologic disasters that threat human life and property in every country. In addition, landslides are one of the most important natural phe-nomena, which directly or indirectly affect countries’ econ-omy. Turkey is also the country that is under the threat of landslides. Landslides frequently occur in all of the Black Sea region as well as in many parts of Marmara, East Anatolia, and Mediterranean regions. Since these landslides re-sulted in destruction, they are ranked as the second important natural phenomenon that comes after earthquake in Turkey. In recent years several landslides happened after heavy rains and the resulting ﬂoods. This makes the landslide monitoring and mitigation techniques an important study subject for the related professional disciplines in Turkey. The investigations on surface deformations are conducted to deﬁne the bound-aries of the landslide, size, level of activity and direction(s) of the movement, and to determine individual moving blocks of the main slide. This study focuses on the use of a kinematic deformation analysis based on Kalman Filtering at a landslide area near Istanbul. Kinematic deformation analysis has been applied in a landslide area, which is located to the north of Istanbul city. Positional data were collected using GPS technique. As part of the study, conventional static deformation analysis methodology has also been applied on the same data. The results and comparisons are discussed in this paper.


Introduction
Landslide is a major type of natural hazards killing or injuring a large number of individuals and creating very high costs every year. So, deformation measurements on landslides are Correspondence to: M. T. Ozludemir (tozlu@itu.edu.tr) today a very important task in engineering geodesy (Haberler, 2004;Haberler-Weber, 2005).
The understanding of the behaviour of landslides and the identification of their possible triggering effects (seismic, hydrological) usually requires a good knowledge of the surface and subsurface kinematics of the sliding landmasses. Furthermore, the mitigation of the landslide risks requires the establishment of monitoring systems which can detect early indications of rapid, catastrophic failures, and enable effective stabilization measures. The study of the kinematics of landslides and their monitoring are usually based on geotechnical and geodetic methods (Stiros et al., 2004).
Turkey is a country under the risk of landslides. In recent years several landslides happened after heavy rains and the resulting floods. This makes the landslide monitoring and mitigation techniques an important study subject for the related professional disciplines in Turkey.
In the periphery of Istanbul, which is a metropolis with a population of more than 10 million, there are also some landslide regions. The area considered in this study is one of these landslide regions located nearby Gürpinar village that is to the northwest of Istanbul. In the area, without doing proper geotechnical investigations some buildings, mostly weekend houses have been built. But after the construction work had completed, many damages on the constructions took place as a result of the landslides. In order to investigate the effects of landslides in and around the settlement area, a multidisciplinary project has been realised. Throughout the project geotechnical investigations and geodetic deformation measurements and analysis were done (Acar et al., 2003).
According to the geotechnical investigations, excluding the artificial disturbance of natural equilibrium, the reason of soil movements depends upon the changes in conditions related to underground water, seismic forces arising after earthquakes and the decrease in sliding strength in fissured (capilar fissures) and heavily consolidated clays. The area where the study was carried out is an old landslide region where  (Altan et al., 1994).
slope equilibrium was formed gradually. As a result of movements in the form of small-sized flours due to surface water and the settlement in the area, the equilibrium was disturbed, causing mass movements ( Fig. 1). In order to prevent soil movements, drainage was performed (Acar et al., 2004;Altan et al., 1994).
As seen from the pictures in Fig. 2, following the landslides there were severe damages on the buildings.
The landslide monitoring project was carried out to identify the characteristics of the landslides and their potential risks in the study area. The geodetic measurement campaigns consist of two phases. In the first part of the monitoring, terrestrial geodetic observations were done in 6 epochs between 1990 and 1991. By the time, some network points were lost. In the second phase, the remaining network points were measured using the GPS technique in 4 epochs between 1996 and 1998. The results of the terrestrial works had been published before (Altan et al., 1994). In this study, only the campaigns with GPS observations are evaluated.
As part of landslide monitoring, different deformation analysis algorithms have been applied so far. In this study, a deformation analysis procedure that has been carried out through a kinematic deformation model based on Kalman filtering technique is discussed. The kinematical motion model, in which GPS derived station coordinates are taken as input values has been evaluated by Kalman filtering technique and motion parameters of network stations have been determined.
Then movement determination was made with static and kinematical models. Only moved points and movement quantities were determined with static model. In addition to the point positions, velocities and accelerations of point positions were also determined with kinematical model being a time dependent function using Kalman-Filter technique (Yalçinkaya, 2003).
In the following sections, deformation analysis procedure and numerical applications are given.

Deformation analysis with S-transformation
In geodetic adjustment, a constrained network adjustment procedure is applied when the datum parameters are known in advance. In order to determine these parameters in advance, some coordinate and measurement values are considered as errorless. Since these values contain errors, it is clear that the results of the adjustment would also contain errors. In case of the constrained network adjustment, the measurements would also be affected with these constraints. Therefore the post-adjustment cofactor matrix does not contain real information regarding the inner accuracy of the network. In addition point errors get higher when getting further from the known points. Because of these reasons, in the networks established for deformation monitoring, free network adjustment procedure is chosen.
Free network adjustment reflects the inner accuracy of the network more realistically and also its external parameters do not depend on certain assumptions. The disadvantageous character of constrained adjustment is removed by the application of free network adjustment, in which the external parameters of the network are computed through the adjustment computations (Tanir, 2000).
For the application of deformation analysis in geodetic networks and for inter-comparison between their accuracies, the regarding parameters to be compared should be determined in the same datum. The datum consistency of geodetic networks measured and computed in different periods can be achieved by S-transformation without performing another adjustment procedure. Moreover the displaced points in the network can also be determined through S-transformation (Demirel, 1987;Denli, 1998;Doǧanalp et al., 2005;Erol and Ayan, 2003;Erol et al., 2005;Inal and Ceylan, 2002). S-Transformation from datum k to datum j can be done by the following formulae.
where x j and x k are the coordinate vectors j and k, S j transformation matrix, Q j xx and Q k xx are cofactor matrices.
where E j is the datum determiner matrix of which diagonal elements are 1 for the datum determiner points and 0 for the other points.; I is the identity matrix; and G is the matrix containing the conditions for a free network to allow for the computation of the coordinates (Illner, 1983 and1985;Niemeier, 1985;Welsch, 1993;Denli, 2004;Erol et al., 2005;Inal and Yiǧit, 2006).

Global test through S-transformation
After the evaluation of the measurements made in different periods, the global test includes only the network sections covered by the identical points. The coordinates of network points based on the measurements collected at two different periods, i.e. t j and t k , can be given as follows: where x k c and x k a are coordinates of datum and non-datum points at period k, respectively. Q k xx cofactor matrix of parameter vector x k can be written as follows: Transformation from datum k to j with the help of E j matrix and positioning of the network with respect to the datum points are accomplished by the following formulae: These computations are done separately for each period. Global test for the coordinates of datum points and their cofactor matrix at datum j is performed as follows: If the degree of freedom of R c is shown as h c =u c −d, then, If the test value T is greater then the critical value, i.e. if T >F h c ,f,1−α , then there are significant deformations in the common points. The next step is then localisation of these deformations Erol et al., 2005;Denli, 1998).

Investigation of significant point displacements by Stransformation
The next step after the failing of global test is to determine which points have deformations. Considering that each datum point has moved, sub-vector x k c containing the coordinates of common points in the parameter vector x k , which are determined through free adjustment at datum k, is partitioned into sub-vectors of x k h , which contains point coordinates of a possibly moving point, and x k s containing the other common points that are assumed to be fixed. Since vector x k a contain the other parameters regarding the non-common points and other unknowns, x k vector and Q k xx cofactor matrix can then be arranged as follows (Denli, 1998;Inal and Ceylan, 2002): The network measured in a time t should be positioned with respect to the points, of which coordinates are in x k s and assumed to be fixed. Denoting this datum by j , the Stransformation given in Eqs. (1) and (2) is performed.
The null hypothesis for the points assumed to be fixed is given below: For testing the null hypothesis, with respect to the Eqs. (9), (10), (11) and (12), coordinate differences d s , their cofactor matrix (Q dd ) s and the increase value for residuals in quadratic form R s can be given as follows: The procedure given in Eqs. (15) to (20) is repeated for each point in x k c vector. The point that has the minimum R s value is tested.
Test value for the point with minimum R s value is computed as follows.
If the test value is greater than the critical value, F h s ,f,1−α , the displacement of this point is considered as statistically significant. This point is then excluded from the datum points and included in vector x k a , which contains non-common points. The process, given in Sect. 2.2, is repeated for the remaining datum points. Investigation on other moving points is kept going on until no significant movement is detected (Demirel, 1987;Denli, 1998;Inal and Ceylan, 2002).

Kinematic deformation model
Kalman filtering is an optimal estimation method to analyse a dynamic system and was developed in early 1960s (Celik et al., 2006). In addition, Kalman filtering is an important tool for deformation analysis combining information on object behaviour and measurement quantities (Kuhlmann, 2003). The intention of kinematic models is to find a suitable description of point movements by time functions without regarding the potential relationship to causative forces. Polynomial approaches, especially velocities and accelerations, and harmonic functions are commonly applied (Welsch and Heunecke, 2001).
A time-dependent 3-D kinematic model that contains position, velocity and acceleration can be expressed by the following formula: j : coordinate of point j at time (t k ) period v xj , v yj , v zj : velocities of X, Y, Z coordinates of point j a xj , a yj , a zj : accelerations of X, Y, Z coordinates of point j k=1, 2, . . . , i (i: measurement period number) j =1, 2, . . . , n (n: number of points) Kalman filtering technique is employed for the prediction of present state vector using state vector information of known motion parameters at period t k and the measurements collected at period t k+1 . The state vector of motion parameters consists of position, motion and acceleration variables. The motion and acceleration parameters are the first and the second derivations of the position with respect to time. The matrix form of the motion model used for the prediction of motion parameters by Kalman filtering technique in 3-D networks can be given as follows:  (Acar et al., 2003;Yalçinkaya and Bayrak, 2003;Yalçinkaya and Bayrak, 2005). The system noise is considered as the noise matrix S that consists of the terms of the last column given in Eq. (24).
where α k : the random noise vector between periods t k+1 and t k Q Y − Y − ,k : the cofactor matrix of state vector at time t k Q ww,k : the cofactor matrix of system noises at time t k The random noise vector α is uncertain and as a rule it cannot be measured. Therefore, for α, pseudo observation vector can be used as α=0. The effect of noise on positions can be determined from former experiences. Its effect on the motion and acceleration, however, can be hardly predicted (Acar et al., 2004;Yalçinkaya and Bayrak;2002).
The residual vector for the observations at period k+1 is formed as follows:  (25) and (27) as follows (Acar et al., 2003;Acar et al., 2004;Yalçinkaya, 2003;Yalçinkaya and Bayrak, 2005): By this model, motion parameters and cofactor matrix are computed. Kalman gain matrix is given as follows: Using the equations above, innovation vector d k+1 , state vector filtered at time t k+1 ; Y − k+1 , predicted state vector; vȲ ,k+1 , residual vector of observations at time t k+1 can be computed by the following equation: Actually, the filtering phase is based on classical least squares adjustment. The most important difference from the classical adjustment procedure is that, contrary to the classical approach, in the filtering the number of observations can be less than the number of unknowns. Through the filtering, adjusted values of state unknowns are computed using weighted combination of measurements and a priori estimations (Acar et al, 2004;Yalçinkaya and Bayrak, 2002).

Numerical application
In order to determine the point displacements in the landslide area, a deformation network consisting of 13 points was set up. The control points were established in stable areas out of the landslide region. The locations of the deformation points were determined according to the geotechnical investigations in the landslide region. The geodetic measurements used in the project were GPS measurements that were carried out in four periods between July 1996 and April 1998. In this study, the measurements of periods March 1997, October 1997 and April 1998 have been evaluated. The GPS data was collected in rapid static mode using 6 Leica SR399 and 4 Trimble SSI receivers. In all periods 2 sessions of GPS observations (10 min at each point) were realised. The data has been processed using commercial Leica SKI-Pro software. The measurements in each period have been adjusted through free network adjustment procedure and their adjusted values and variance-covariance matrix have been computed (Acar et al., 2003). As a first step, a static deformation analysis was carried out through the evaluation of adjusted coordinates and their variance-covariance information. In the analysis of all periods, Codeka3D deformation analysis software was used. The results are given in Table 1 As seen from the results, significant displacements in the positions of 5 points have been detected. Out of these points, point 2396 has the largest displacements in its coordinate components.
In the second part of the study, the kinematic deformation analysis procedure based on Kalman filtering technique as described in the previous section has been applied. By the application of this model, not only the point displacements but also the motion parameters of network points have been computed. Afterwards, it has been tested whether or not the obtained results are statistically significant. If parameters have significantly changed in kinematical model, a " √ " sign is given in Table 2 and 3, otherwise a "−"sign. Table 2 shows the results of kinematic deformation analysis for the period between March 1997 and October 1997. The results are almost identical with the static analysis results given in Table 1. In Table 3, the results of kinematic analysis for the period between March 1997 and March 1998 are given. In these tables, point displacements, velocities and accelerations are given. Figures 3 and 4 show horizontal and vertical displacements in the period between March 1997 andMarch 1998, respectively. According to the results given above, kinematic model yields almost identical results with the static deformation analysis. However, additional parameters, i.e. time depen- dent motion parameters, have also been computed by the application of kinematic deformation analysis.

Conclusions
In this study, a kinematic deformation analysis procedure based on Kalman filtering technique has been applied on a data set collected in a landslide area by GPS. In addition to this technique, the data has also been analyzed by static deformation analysis. Two different approaches yielded identical results. However, the kinematic model has some clear advantages. For example, in kinematic model time dependent motion parameters of each point can be determined.
Stepwise computation of motion parameters eases the control of the computations and the interpretation of the results. It is obvious that, for the computation of motion parameters or in other words for modelling the motion, more measurements are required. This is actually the main drawback of kinematic deformation model approach. In this study, in order to overcome this problem, Kalman filtering technique has been conducted for the computation of motion parameters. The main advantage of Kalman filtering technique is that it requires less measurement period. However, since the Kalman filtering technique employs prediction, the kinematic behaviours should not be extended unlimitedly by extrapolation.
The study area discussed in this paper is a landslide area where a multi-disciplinary project had been conducted. The project partners are geodesists, civil engineers and earth scientists. However, this study focused only on the geodetic deformation monitoring process. It is clear that, through the combination of different data sets, a more realistic deformation model for the landslides would be produced.
Following the landslide monitoring project, many buildings in the landslide area were destroyed. In addition, some preventive measures were taken in and around the landslide area. The authorities do no longer give building permission in the study area.