Multi-sensor observation fusion scheme based on 3D variational assimilation for landslide monitoring

Abstract Multi-sensor observation is very important for monitoring landslide disasters. Since various surveying techniques are currently available for detecting variational slope activities from different perspectives, studies have focused on integration of multi-source information for the analysing landslide displacements. In this study, a general multi-source data fusion scheme for landslide monitoring based on three-dimensional variation (3DVar) data assimilation was developed. The scheme was used to fuse different observations of Xishancun Landslide in Li County, Sichuan Province, China. First, the displacement observations obtained by a Global Positioning System (GPS) and Borehole Inclinometers (BIs) were assimilated for accurate evaluation of slope activities. Then, slope Stability Index (SI) was introduced to validate the assimilation results within a time interval. SIAssi values calculated using the integration model developed in the present study were compared with SIFS simulated by a physically based landslide model. The correlation coefficient between them ss 0.75, which is larger than those with SIGPS (0.45) or SIBIs (0.41) values determined by the GPS and BIs respectively. The assimilation results are thus confirmed to be more credible for slope stability simulation.


Introduction
A landslide is a hazardous geological phenomenon that frequently occurs globally (Schuster 1996;Pradhan and Buchroithner 2012). Landslides pose a threat to millions of lives and diverse structures and infrastructure such as buildings, roads, and reservoirs (Pourghasemi et al. 2012;Lu et al. 2015). They sometimes occur on a very large scale, resulting in heavy casualties and economic losses. As a landslide-prone country, China currently contends with the threat of great landslides, especially in the south-western mountainous region (Li et al. 2016a;Li et al. 2017). It has been reported that there are presently more than 90,000 potential landslides threating about 70 cities and villages of China with at least 10% of the residents in direct danger (Zhou et al. 2005;Huang and Li 2011;Li et al. 2016b). It is very important to construct the systematic platform for monitoring landslide processes and identifying landslide characteristics.
A variety of monitoring techniques for slope failure survey have been developed and used for landslide monitoring in recent years (Gili et al. 2000;Travelletti et al. 2012;Zhao et al. 2018). Among them, remote sensing techniques such as those that utilize a synthetic aperture radar (SAR) ) and very high resolution (VHR) satellite images (Debella-Gilo and K€ a€ ab 2011) are used to monitor regional slope deformation. However, their low temporal resolution, which is due to the long revisit cycle of the utilized satellite systems, may result in the inability to capture slope failure variations (Greif and Vlcko 2012). Nevertheless, some in situ devices are particularly known for their high-temporal and continuous sequential observations. For example, the Global Positioning System (GPS) has been increasingly employed in the study of landslide surface movements (Gili et al. 2000;Wang 2011;Huang et al. 2017), but its accuracy is affected by some inherent errors such as clock bias, ionosphere error, and troposphere error (Parkinson 1996). Underground inclinometers are used to determine internal slope changes (Chang et al. 2005;Segalini and Carini 2013) but the used inclinometers to calculate the surface displacement may also cause some uncertainties. These devices are thus limitedly effective for slope failure survey through direct monitoring. Accordingly, multi-sensor networks have been developed for gathering diverse information from underground layers and determining the topographic characteristics of the slope surface and the meteorological conditions (Qiao et al. 2013). Multi-source data can be used to enhance the effectiveness of sensors and achieve better research analyses from different perspectives (Lee and Pradhan 2007;Dewitte et al. 2008;Corsini et al. 2013;Bardi et al. 2016). Landslide monitoring using multi-source data particularly facilitates the illustration of the underlying cause of the slope activities (Du et al. 2013;Xia et al. 2013;Jiang et al. 2016). However, multi-source data may also contain significant discrepancies regarded to data scale, resolution, format, and other factors, and it may complicate the data processing. Moreover, observation errors and uncertainties negatively affect the accuracies of the results, leading to a biased estimation of the slope stability (Liao et al. 2012).
In this study, we developed a multi-sensor observation fusion scheme based on three-dimensional variational (3DVar) assimilation to fuse two types of data, namely, GPS and Borehole Inclinometers (BIs) data. First, 3DVar assimilation method is mathematically premised on optimal state estimation, in which historical data is used to reduce analysis error under constraints of dynamic consistency to infer present state through best solutions (Dee and Da Silva 1998;Renzullo et al. 2008). This method was first applied in the fields of meteorology and oceanography (Charney et al. 1969;Thompson 1969). So far, variational assimilation has been implemented and investigated in numerous studies such as those on weather prediction (Barker et al. 2004), wind field modelling (Ulazia et al. 2016), land surface temperature inversion (Carrera et al. 2016), and remote sensing of soil moisture (Mechri et al. 2016). On the other hand, there are kinds of errors and inherited weaknesses in landslide observations, which may bring disadvantages to slope stability estimation. For example, random errors in BIs could be smaller than 1mm but the their measurement magnitude is as much as a few centimetres per 30 m deep (Stark and Choi 2008); on the contrary, GPS affords millimetre-accurate landslide movement observations (Wang 2013), but the rough results contain more random noises in related to ionospheric and tropospheric influences (Huang et al. 2016).
The rest of this paper is organized as follows. Section 2 firstly presents the general framework to fuse multi-source observations based on 3DVar assimilation method. Section 3 introduces study area and relevant data including GPS and BIs data. Results are discussed and summarized in Section 4 and Section 5, respectively.

Observation fusion method
2.1. General framework Figure 1 illustrates the multi-sensor observation fusion scheme developed in the present study. Considering the errors of historical multi-sensor data, the cost functions derived from the parameters are implemented to generate the assimilated data, which contains corrections for the current epoch based on the modelled error information and real-time information.

Assumption and error modelling
It is assumed that the error (e) values of the multi-sensor observation at a location (x, y) in a long time series are normally distributed with the density function f (e | l, r), where l is the expectation, and r is the square root of the variance. As expressed by Equation (1), e is mutually independent. These parameters are modelled using a monotonic (e.g. linear/non-linear or smooth) link function, and the error distribution can be subsequently determined.

Cost function obtained by 3DVar assimilation
3DVar is a classic data assimilation method that involves the minimization of a prescribed cost function (Ide et al. 1997). It provides an optimal solution by considering the fields of the background and observation data. Based on 3DVar assimilation and the modelled parameters (l, r), the cost function can be formulated as follows: where X denotes the location of the desired data, and Z i is the ith kind of observation.
In the case of two input data series, the 3DVar cost function can be formulated as follows: where J b is the cost function value for the analysis data (X a ) and background data (X b ), and J o is the value for the analysis data and observation data (Y o ) in Figure 2. The matrixes B and R respectively represent the error variances of X b and Y o , and H is the linear projection function of Y o into the dimensional space of X b . H is an identity matrix when X b and Y o are in the same dimensional space.
To seek for the solution (analysis values X a ) of the assimilation problem in Equation (3), the differential Equation (4) is used to minimize the cost function. Equation (5) expresses the iterative method for solving the problem.

Description of experimental area
The study area of the present work is located in South-western China, east of the Tibetan Plateau (Figure 3). It is 22 km from Wenchuan City, and is one of the most prone to severe geo-hazard areas, after being exposed to the 5.12 Wenchuan Earthquake in 2008 (Cui et al. 2011;Dai et al. 2011). The specifically considered 'Xishancun Landslide' located north of the Zagunao River is a large soil accumulation with a relative elevation of about 1790 m, slope length of about 3800 m, and average width of approximately 700 m (Xu et al. 2016).

Equipment installation and data calibration
To detect the real-time slope displacement, three types of monitoring devices were installed on the slope surface (Figure 4), namely, a Global Positioning System, Borehole Inclinometers (BIs), and two rain gauges (RG). Some local instabilities such as positive and negative road subsidence and surface debris could occur in the neighbour of these sensors. All the devices are powered by solar energy so that they continuously work in 24-h 1 day. The data obtained by all the sensors is transmitted to a data centre through a General Packet Radio Service (GPRS). Table 1 gives the details of the three types of employed monitoring devices. The BIs and RG have the same data collection interval of 10 min, while the corresponding interval for the GPS is 2 h. Due to minimal data loss during the GPRS data transmission, the numbers of available data records were 2,560 for the GPS, 24,883 for the BIs, and 25,104 for the each RG. The duration time of the data acquisition is from January 1st, 2016 to November 17th, 2016.
After difference calibration, the horizontal and vertical accuracies of the GPS data are respectively 2.5 mm ±1 ppm and 5 mm ±1ppm. BIs comprised six inclinometers connected by an aluminium alloy tube, with the one located deepest underground being at a depth of 42.78 m ( Figure 5). Although BIs may not directly measure the horizontal movements of the casing (Stark and Choi, 2008), the tilt angle could be converted into horizontal displacement, as illustrated in Figure 5. The following equation is used. where L i is the distance between two neighbouring BIs. Because the tilt angle h i is rather small for each movement, L i is considered to be close to the initially set pitch (K in Figure 5). Equation (6) could be used to determine the incremental displacement at each inclinometer and the cumulative horizontal displacement. Figure 6 shows the coordinate calibrations of the BIs and GPS. The northward and eastward coordinates of the GPS are respectively denoted by X and Y, while the coordinates of the BIs in the horizontal plane are denoted by A and B (Figure 6(a)). There is a rotation angle h between the two coordinate systems in Figure 6(b). The resultant horizontal displacement is calculated using Equation (7).

Possible displacement based on GPS and BIs error estimation
The variations of the slope surface displacements determined by the GPS and BIs after space and time calibration are shown in Figure 7. GPS has been proven enable millimetre-accurate recording of slope failure, but more random errors, introduced during signal transmission and by interferences, which still need to be carefully mitigated to achieve millimetre accuracy (Wang 2013). The process of it could be complicated sometimes in related to factors of ionospheric and tropospheric influences (Huang et al. 2016). Conversely, while BIs only afford measurement magnitudes, being as much as a few centimetres per 30 m (Stark and Choi 2008), their associated random errors, as in situ instruments, are smaller. Therefore, considering insufficient mitigation of random errors, the secondary errors from GPS are larger than from BIs.  To verify the assumptions mentioned in Section 2.2, that the random errors of the GPS and BIs observations are normally distributed, the moving average (MA) method was used to fit the 'true' data for a 24 h duration. The obtained errors distributions for the GPS and BIs are shown in Figure 8. As can be observed, the errors for the two types of equipment range within À5 to 5 mm and À2 to 2 mm, respectively.

Assimilation results for GPS and BIs
Based on the above analysis of the probable error distribution, two monitoring data were input to a 3DVar assimilation model, with the background (observation) field input being the data obtained by the GPS (or BIs). Figure 9 shows the cumulative slope surface displacements obtained by the GPS, BIs, and assimilation results. As can be observed, the sequence of the assimilation results is much closer to that of the GPS sequence at the 20th iteration step. With continuous assimilation, it subsequently becomes closer to the BIs sequence. In addition, the noise in the assimilation results gradually attenuates with increasing number of iterations.
With increasing number of iterations, the derived cost function J finally coverages to 1410.2147 from the initial value of 4906.4473 ( Figure 10). Considering that a threshold of 0.001 selected as the terminal iteration signal, the final iteration step (the 92nd iteration) degrades to 0.0009 from 111.4385 at the 20th iteration. emphasizes on the rainfall-induced micro-deformations of the slope in the simulation of the slope stability as follows.
where DL þ (or DL À ) is the load (or unload) increment, DR þ (or DR À ) is the load (or unload) response increment. In our test, the rainfall intensity is regarded as the load (L) acting on the slope system, and the induced displacement as the response (R) to the load.
In addition, c þ (or c À ) is defined as the rainfall-induced deformation rate due to the load (or unload). In the elastic-plastic slope deformation state, the rate of the deformation response to the load (c þ ), is mostly equal to that of the unload (c À ). However, the deference between c þ and c À increases with increasing load. If the slope system is continuously loaded to destabilization, the load deformation ratio would tend to infinity (c þ ! 1) while c À would remain constant. Under this condition, a very small load could induce a very large displacement response, and g would be close to 0. The slope is thus unstable when g is close to 0. In the present study, g is normalized as the stability index (SI) with a value range of 0-1. Figure 11 shows the variations of the daily rainfall and SI between October 4th and 14th, 2016. As can be observed, the daily rainfall is maximum at 7.8 mm/day on October 11th, and minimum at 1.0 mm/day on October 6th. In addition, more unstable processes may occur after heavier rainfall, such as on October 6th and 12th. Variations of SI extracted by three kinds of data indicate different changes in the test time window. BIs show that the highest unstable process occurs at the 12th. Moreover, SI BIs exhibits more days (5 days) to be lower than SI GPS and SI Assi .

Test of assimilation results based on slope mechanism
The slope stability determined by different methods may differ, resulting in errors in the monitoring of the slope failure. In further investigation, an attempt is made to incorporate the mechanism of the slope stability into the analysis to verify the present assimilation results.
A physical model referred as the Slope-Infiltration-Distributed Equilibrium (SLIDE) model is used to model the slope stability. SLIDE is a shallow rainfallinduced landslide model. Based on the rainfall infiltrating process, it simplifies the saturated-unsaturated seepage equation (Richards 1931) and integrates the contribution of the apparent cohesion to the shear strength of the soil and water level when the surface of the soil attains saturation (Liao et al. 2012;He et al. 2016). The detailed theory is available in the work of Liao et al. (2012). The model simulates the variations of the slope factor of safety (FS) in Equation (10). The value of the SI FS can be calculated using Equation (11).
The SLIDE model has been previously used to simulate rainfall-induced shallow landslides in Li County, where the Xishancun Landslide is located (Li et al. 2017).

Figure 11
Rainfall and slope stability based on the ULDRR model.
Therefore, soil parameters in the previous work were adopted in the present study to determine the variables in Equation (10). As can be observed from Figure 12, the variation of SI FS over the entire test duration is closer to SI Assi than those of SI GPS and SI BIs even when those are divergent for example on 7th and 12th. Further, the assimilation results regarding to slope failure are more accurate in SI Assi than in SI GPS and SI BIs . For example, SI GPS is 0.03 on 8th and SI BIs is 0.01 on 12th, this kind of low values generally indicated a slop failure but it did not happen actually. Therefore, based on the analysis above, the assimilation results are thus confirmed to be more credible for slope stability simulation, by integrating the advantages afforded by the use of both GPS and BIs.
Furthermore, we have compared the assimilated results with the traditional fusion method-polynomial method (PM). In the fitted PM, the data from GPS and BIs is considered as sample data during 4th-10th, and other data during 11th-14th as tested data. Meanwhile, SI PM from each tested data have been estimated by PM. As shown in Table. 2, the assimilated results (SI Assi ¼ 0.28, 0.23, 0.30, and 0.74) perform better than the results from PM in each estimated four days (SI PM ¼ 0.49, 0.36, 0.39, and 0.39), and are closer to the results from SLIDE model (SI FS ).

Summary and conclusion
Multi-source observation is an advanced scientific method for monitoring landslides (Herrera et al. 2017). However, it faces challenge of the need to fuse multiple heterogeneous observation data for processing (Gravina et al. 2017). In the present study, a multi-sensor observation fusion scheme based on 3DVar assimilation is developed for the fusion of GPS and BIs data. Some major conclusions of this study can be summarized as follows: 1. Observations of the two considered sensors may differ for serial slope displacements, due to internal and external factors such as systematic data errors, noise, and data frequency. After estimating the observation errors, the cost function is constructed initially and iteration steps of it converge to less than 0.001 mm through the proposed multi-sensor observation fusion scheme 2. The results show that the proposed scheme based on 3DVar stands out in fusing multiple heterogeneous data. The utilized assimilation process particularly enables more accurate monitoring of slope failure. During a typical period with continuous rainfall, the slope stability extracted from the assimilation results (SI Assi ) was observed to be closer to that determined by simulation using a physical model than those observed by GPS (SI GPS ) and BIs (SI BIs ), respectively.
In further study, we will attempt to use the proposed scheme to fuse additional remotely sensed data such as those obtained by Interferometry Synthetic Radar (In-SAR) and Light Detection and Ranging (LiDAR). This is expected to enrich the acquired data to achieve spatial-temporal assimilation (Hunt et al. 2007).

Disclosure statement
No potential conflict of interest was reported by the authors.