Mathematical Modeling of Landslide Slope Dynamics and Sustainable Development of Territories of Asian Russia

The possibilities of applying the results of modeling the landslide stress-strain state to predict landslide hazard and make management decisions under intensive and spatial development of the territories of Asian Russia are considered. Economic and spatial development of the territories concerning reorganization or formation of economically viable components of various production chains require careful consideration of their safe and continuous operation. In this respect, great attention is paid to the studies on exogenous geological processes under the complex geomorphic conditions that can have a substantial impact on the operation of production chain links. The authors present a technique for modeling a landslide slope movement caused by anthropogenic impact. The revealed patterns in which influencing factors affect vertical displacements of landslide points created the conditions for their description at a new qualitative level. Mathematical modeling of landslide slope movement is based on the dynamic-type prediction model. This made it possible to predict the displacements of observed landslide points with sufficient accuracy. Based on a geo-cognitive safety monitoring of production chain links primarily important for the social and economic development of territories, the possibility of using the prediction models of the landslide slope stress-strain behavior as a "solver" in expert systems for forecasting landslide danger when making management decisions is shown.


Introduction
Cutting-edge management models for regional and inter-regional economic and spatial development call for new scientifically valid methods and technologies.
The best practices in economic management prove that competitive strength can be achieved through the processes of horizontal and vertical regional integration where the development is triggered by the groups of companies and production chains rather than individual enterprises. Such an approach makes it possible to improve production and territorial management systems on various levels (Diao et al. 2009;Bell and Morse 2018;Han et al. 2021) but requires taking into consideration numerous factors that can slow down or even knock out the whole production chain for a long while. Considering that in Asian Russia the number of duplicate links of production chains is rather limited, control, monitoring, and safety of their operation become particularly important.
Presently, when numerous programs for the development of territories in Asian Russia are being elaborated, structured spatial data allowing automation and significant improvement of the quality of management decisions are of great demand. Such decisions, based on spatial expert systems (ES) and knowledge bases, make it possible to reveal the attractiveness of a certain territory, its economic capacity, possible weaknesses and environmental factors affecting the object under study, as well as establish their cause-effect relations (Diao et al. 2009;Atif et al. 2021).
One of the dangerous factors that may highly negatively impact on a production chain link is geological phenomena -the most common natural and anthropogenic processes able to cease the operation of not only a single link but the whole production chain, causing major damage to regional or inter-regional industrial cluster. Hungr et al. (2014) believe that formation, identification and monitoring mechanisms, including timely identification of dangerous areas and mitigation of landslide consequences are the most important components for carrying out a qualitative assessment of the dynamics of natural and anthropogenic processes. Schuster and Krizek (1978), and Van Den Eeckhaut et al. (2012) admit that the nature of landslide dynamics is exceptional and may be systematized and classified only within the limits of the monitoring zone having similar physical and geological characteristics.
The significance of complex research and mathematical modeling of slope processes is immense however, the engagement of geodetic observation methods is still underestimated. The latter would encompass not only observational data on downslope mass transfer, but some topographic information about land relief, slope orientation, other landscape properties, which have a significant impact on the downslope mass transfer. Simonyan (2011) summarizes and substantiates geodetic methods for studying landslide trends that allow carrying out high-degree control of different movements within general landslide areas, as well as their peculiarities. At the same time, geodetic data on movement compilation of monitored landslide points constitute an extraordinary valuable integral and quantitative characteristic is a result of a complex number of factor interaction, which defines the landslide behavior and dynamics. Such data are the basis for identifying certain patterns and features of a landslide process development. Moreover, the topography of land relief reflects the degree to which the gravitational forces influence landslide movement, while the orientation and other properties of the landslide define insolation level and seasonal pattern of the process dynamics.
Along with geodetic observations over landslide-prone slopes, there are other approaches and methods to register their movement and dynamics. At the same time, the type of a landslide process and the mechanism of landslide mass displacement define the choice of calculation equations, which is consequential to its geologic structure and soil elastic properties. Baborykin et al. (2015) describe some of the currently used methods for predicting landslide dynamics based on mathematical analysis under a number of assessment criteria of a digital terrain model (DTM) taken with certain weights.
An interdisciplinary approach that uses geological, geomorphological and geotechnical data together with GIS capabilities (Federici et al. 2007;Ali et al. 2019) and includes data on a slope geometry and its direction, as well as landslide movement type and degree of its activity, along with the statistical analysis of a number of factors associated with the landslide and various efficient algorithms for its study, enables spatial modeling of certain layers of a landslide slope, identifying the degree of their instability, possible deadlines of their occurrence, movement, etc. Wen et al. (2017) used a least square support vector machine model (GA-LSSVM) together with an algorithm for streamlining model parameters to predict dynamics of landslide hazards in the vicinity of the Three Gorges Dam Reservoir; hence, trustworthy results of early prediction were obtained based on the case of Shuping landslide.
Analysis of literature shows that methods of monitoring and predicting means for landslide hazards are being developed in two areas: development of methods to predict the proneness of the territories under study to landslide hazard, and building efficient models for predicting dynamics of certain landslides. In a number of cases, the development of secondtype forecast models provides better results because they can predict future landslide behavior.
The issues related to mapping landslide-prone territories in the regions of Asian Russia are most relevant considering their geographical location, amounts of precipitation, significant temperature drops, rugged terrain, possible earthquakes, etc. That is why numerous research is devoted to efficient risk mitigation: implementation of machine learning methods and creation of algorithms for automatic detection of geologically hazardous areas. In this respect, the further types of machine learning models have proved themselves quite worthy. Spatial modeling of landslide hazards based on statistical analysis of the exposure factors used for solving multidimensional and nonlinear hazard prediction problems by building models based on decision trees is offered by Hosseinalizadeh et al. (2019). We consider support vector machine model (SVM) (Yousefi et al. 2015(Yousefi et al. , 2020 together with the maximum entropy method (Boogar et al. 2019) as high-quality machine learning algorithms used for mapping landslide-prone areas with a quite high level of prediction accuracy. Greco et al. (2007), Yousefi et al. (2020), and Pourghasemi et al. (2020) analyzed the capabilities of a generalized linear model of multiple regression (GLM) and logistic regression model (LR). A model of boosted regression tree was developed by Yousefi et al. (2020), while Pourghasemi and Rossi (2017) used multivariate adaptive regression splines (MARS). To predict landslide hazards Yousefi et al. (2020) and Huang et al. (2018) also used the model of modified analytical hierarchy process (M-AHP), support vector machine model (SVM), frequency ratio model (FR) and the model of artificial neural network ); a model based on the algorithm of fuzzy logic theory was offered by Akgun et al. (2012). It is important to note that there is a certain trend in applying an integrated approach to detecting landslide-prone areas that combines empirical and physical research methods Hsu et al. 2018). Bearing that in mind, we can see a connection between empirical landslide parameters (e.g. landslide body volume, terrain topography, geological and morphological distinctive feature) and physical models integrated in GIS (e.g. AHP, WoE, LR models) for further landslide susceptibility analysis (LSA).
Presently some papers provided profound research on identification of efficient methods for detecting landslideprone areas. Thus, after a comprehensive analysis of GLM, MARS GAM, and M-AHP machine learning models Pourghasemi and Rossi (2017), and Gayen et al. (2019) chose the M-AHP model as the most efficient. While studying the efficiency of SVM, MLP, RBF Neural Nets, KLR, and LMT models (Bui et al. 2016) for detecting landslide-prone areas, the MLP model was most efficient, while relatively small landslides were better detected with the KLR model and LMT logistic model tree. Huang et al. (2017) analyzed a number of models to predict landslide dynamics; they showed that an ELM multivariate model has some advantages over a linear and Chaotic PSO-SVM models, since predictions consider such triggers as amounts of rainfall and water level in the reservoir, which contribute to the displacement of a landslide body. Hong et al. (2015) compared two models used for mapping landslide hazards in Yihuang area, China: 2-class Kernel model of logistics regression (KLR) and a model of alternative decision tree; it is emphasized that when modeling landslide hazard, the ADT model ensures better high-quality results. Dou et al. (2020) discuss special aspects of SVM machine learning algorithms and their hybrid modifications for mapping the areas prone to landslide hazards; while comparing the four models, SVM-Boosting hybrid model proved to be the most efficient one. Application features of the two types of neuron network for mapping landslide-prone areas with the help of Levenberg-Marquardt learning algorithm and Bayesian regularization were discussed by Bui et al (2012). The calculated at the learning stage weights of the factors under analysis were later used to determine indices of landslide proneness; here the network model of Bayesian regularization turned out to be more efficient. Pradhan (2013), Kuzin and Sannikova (2016), and Huang et al. (2018) used the advantages of GIS for efficient landslide susceptibility analysis, the latter enables storing and processing topographic data, satellite images, field surveys, as well as features of geology, hydrology, climatology, vegetation, land cover, etc. Data analysis and processing with using various models based on GIS and machine learning algorithms allowed creating thematic layers of object-oriented data, calculating various indices and identifying weights of parameters for landslide susceptibility assessment.
Traditional approaches that use remote sensing technologies for building digital terrain models (DTM), based on joint use of photogrammetric, airborne (ALS) and terrestrial (TLS) laser scanning data, as well as instrumental data from GPS reference point network (Mantovani et al. 2013), provide researchers with tools for high-quality identification of various-type landslides and majority of the surface morphology features (Baborykin and Zhidilyaeva 2014; Shan and Toth 2018). At the same time, coprocessing of ALS data (before and after the event) together with long-period seismic recordings provides detailed landslide dynamics (Yamada et al. 2013). Godone et al. (2018) showed the capabilities of aerial laser scanning (LiDAR) for identifying unsteady slopes, their morphometric filtration, GIS processing and further mapping. Jaboyedoff et al. (2012) consider building precise high-resolution digital elevation models (HRDEM) or 3D models as highly attractive. Carter et al. (2001) showed that their combination with geological data allows predicting landslide areas that are susceptible to displacement, improving the quality of landslide mapping, revealing their morphological features (e.g., benches, slope profile, etc.). As for TLS technology, a number of specialized applications for interpretation of landslide formation mechanism, geotechnical properties of slope stability, their border detection and volume assessment were developed (Rowlands et al. 2003;Dunning et al. 2009;Mazzanti et al. 2018). Lately unmanned aerial vehicles (UAV) became ubiquitous, they allow getting real-time high-resolution aerial images that help identify and describe natural / anthropogenic hazards and terrain topography (Hu et al. 2019;Guo et al. 2020). UAVs equipped with a GPS-positioning system (RTK mode) ensure the capability for building a DTM, carrying out its automated analysis and following identification of terrain specific features in images (Ovsiuchenko and Akopov 2012;Godone et al. 2018). Westoby et al. (2012) and Cignetti et al. (2019) recommend using this technology for small-area landslides, which in combination with a MVS method allows building a high-resolution DTM capable to identify and map many geomorphological features generated by landslide evolution. For instance, Zeybek and Sanlioglu (2020) present their research related to building time-dependent DTMs based on post-processed UAV survey data with the use of radial basis function method to detect landslide displacements, which enables assessing the quality of the built DTMs by comparing them with the results of GPS measurements.
Numerous researchers apply various methods of mathematical modeling for predicting dynamic changes in ground water level to assess landslide stability. Chen et al. (2010) described the use of a SOM-RBFN multinodal model, combining self-organizing map theory (SOM) and radial basis function network (RBFN). Application of artificial neuron network (ANN) taught with the help of Bayesian regularization (BR) was studied by Mohanty et al. (2010). Gundogdu and Guney (2007) carried out spatial analysis of variation in ground water level with the help of multipurpose kriging and further empirical selection of models aligned with experimental models. To assess landslide slope stability Iverson (2000) offered to use a transitional model that has certain allowances related to variation in ground water pressure caused by rainfall; the One of the key methods in studying landslide slope dynamics is numerical mathematical modeling. As source data, many researchers use a number of indices characterizing rock physical and mechanical properties. Thus, Liu et al. (2015) made a detailed analysis of several calculation options for studying slope stability, where they considered limit equilibrium method, elevated limit strength method, and strength reduction method (SRM).The methods differ in their input data and analytical models. As the result of the comparative analysis, the factors of safety, sliding surface shapes and positions were calculated. Based on the information about geometry and slope feature input parameters Kozhogulov et al. (2017) (2013) and Moretti et al. (2015) used signal conversion results together with field survey data to predict the time of landslide occurrence. In a number of cases they managed to associate high-frequency signals with topography changes due to increased stirring of landslide forming materials during its movement. Lin (2015) and Fuchs et al. (2018) applied direct modeling of signals for assessing landslide mass and acceleration, as well as interpretation of sequence of events. The second approach is based on studying high-frequency signals that are more typical for a wide range of different size landslides, provided they are in the close proximity to the scene. Due to a rather complicated mechanism of a signal source, topography, and sideline diverse landslide energy scattering effects, such as friction, cracking, and plastic deformation, this approach is less applicable for further interpretation (Manconi et al. The use of machine learning decision algorithms pre-trained on well-known datasets and data received from regional seismic networks is quite efficient and will allow identifying seismic waves taken from already existing database that covers all possible types of events (Hilbert et al. 2014) and implementing automatic detection of an event type with high accuracy. For instance, Dammeier et al. (2011) presented results of their work where they analyzed signals received from the network of continuously operating seismic stations. Within their research, they used multivariate linear regression methods to calculate characteristics of seismic signals and key parameters of a landslide (drop height, volume, beating, potential energy), which further were verified on known datasets. The obtained results showed that this technique may be used for automatic classification when detecting such events in the future.
Over the past decades, a 3D analysis has been used when calculating a slope stability. Meanwhile a number of 3D slope stability evaluation methods based on numerical simulation and other approaches were offered. Dong, Hu and Song (2018) presented some predicted results for landslide slope stability based on 3D geological and geotechnical technologies and models of numerical simulation, while Sari et al. (2020) evaluated slope stability analysis by using the concept of the ultimate rock equilibrium and calculating the displacement resistance coefficient. Zhang et al. (2016) examined some approaches to detecting landslide process activation and identified key factors that cause landslide mass displacement, which later were used for building physical and mathematical models and carrying out necessary calculations. They also showed strengths and weaknesses of different calculation methods (finite element limit equilibrium method, slices method and finite element strength reduction method) when evaluating the three cases of stability analysis, i.e. natural slope, anchored slope with seepage, and excavation anchored slope. After the calculation of sliding surface depth, safety factor and anchoring effect parameters, the SRM method proved to be most efficient.
Nowadays, in order to detect the beginning of landslide occurrence and movement, special emphasis is given to the development of real-time landslide monitoring systems. An example would be an AI-based monitoring system containing different monitoring elements to obtain the data needed for analysis and DeepAR mathematical probabilistic-type prediction model (Dong et al., 2021)  Special interest is paid to insufficiently studied landslide slope dynamics prone to anthropogenic influence. As a rule, there are only nonsystematic geodetic observations of a landslide movement and a lack of quantitative and descriptive characteristics of spatiotemporal parameters presenting anthropogenic influence. It is obvious that the lack of source data makes the identification of patterns for landslide development particularly difficult. The opportunity to use mathematical models to describe landslide movement taking into consideration such influencing anthropogenic factors as blasting and changes caused by excavated soil creates unique conditions for mathematical modeling when predicting landslide processes. Thus, we believe that any stage of landslide development should be simulated with mathematical models and data of geodetic and monitoring observations over the changes in anthropogenic impacts.
The main idea of the proposed approach is that mathematical modeling of landslide dynamics influenced by natural and anthropogenic factors in the context of limited and insufficient source data could become an efficient "solver" used in spatial expert systems to oversee the state and manage safety of industrial facilities. The research is pursuing the following objectives: (i) develop a technique to single out positive and negative landslide movement components caused by anthropogenic factors; (ii) justify the choice of key influencing factors causing landslide movement, i.e. blasting operations and removal of great amount of soil; (iii) build a dynamic mathematical prediction model represented by two first conditional moment functions of the movement process to predict displacements of the observed landslide points; (iiii) apply the developed technique using it as a "solver" in a spatial ES generating management decisions in order to prevent emergencies at the key links of various production chains which may have a negative impact on the economic sustainability of regional and inter-regional interaction.

Investigated processes: general information
The analyzed landslides developed on the right-side slope of the Angara River, in the construction section of the rockfill dam of the Boguchanskaya hydroelectric power plant (HPP) -a key link of several production chains that have strategic significance. The left riverbank is composed of basaltic rocks, the displacements of fixed geodetic points occurred there only during quarry blasting operations.
To a certain extent, landslides are related to the geology of the right bank where sedimentary rocks of binominal structure slop down to the river at angles varying from 5 to 15º. The upper part of the slope is composed of Ordovician hard rocks while the underlying Cambrian sediments are soft clayed argillite and marlstones prone to plastic deformations caused by changing groundwaters. Thus, the underlying clayed rocks form a sliding surface. Overall, the analyzed landslide zone can be considered homogenous in terms of its geological structure. Natural heterogeneity of landslide displacements in the analyzed area is caused by differences in slope steepness, seasonal climatic changes and tectonic disturbance. According to the classification offered by Hungr et al. (2014) landslides on the right bank slope of the Angara River in the dam section of the Boguchanskaya HPP can be classified as plastic.
The study object was the only landslide slope where geodetic observations (III class leveling) of vertical displacements of two groups of intact landslide points (8021,8022,8023,8024,8025) and (8061, 8062, 8063) ( Fig. 1) were carried out for nine years. The observations covered a 2.5 (along) by 1.5 (across) km strip of the riverbank. Initially, there were reasons to assume that major factors affecting the movement of large soil masses to the earth deposits of the dam and groundwater-level changes were caused by blasting works carried out within the period of geodetic observations. Particular attention should be paid to assessing the impact of a landslide slope unloading due to removal of soil from rock quarries since it has an ambiguous effect on landslide development. On the one hand, removal of soil decreases the weight of bank masses and reduces the slope tangent (active force component), which enhances the slope sustainability. On the other hand, it affects gravitational force and decreases friction making the landslide move faster.
Water-level fluctuations in the Angara and groundwater levels, changes of solar radiation and ambient temperature, precipitation, number and power of blasts, volume of removed soil can also be considered as important factors, which are often interconnected and affect landslide point displacements indirectly. Further studies in this field will give better understanding of a land-sliding mechanism and help substantiate reliable forecasts based on the analysis and assessment of the latter.

Hydrological instability assessment of the right bank of the Angara in the dam landfall area
The general picture and possibility to assess landslides based on accumulated observations can be represented as:  rocks fall towards the Angara at an angle up to 15º, series of tectonic disturbance in the work area create preconditions for the occurrence of landslides;  having low inner friction and soil adhesion coefficients weak rocks are prone to plastic deformations, which determines scale, place and time of landslide occurrence depending on natural and anthropogenic impacts;  due to significant heterogeneity of the lithological composition of the rocks, their jointing, degree and time of water intrusion, as well as the influence of adverse stochastic factors, deterministic forecast of the location, scale and duration of a landslide event applying the known methods of rock and soil mechanics is impossible.
To provide reliable assessment and prediction of landslide behavior based on quantitative data of geodetic and other observations it is necessary to carry out a qualitative assessment of the impact of a number of natural and anthropogenic factors upon general regularities of landslide development. In this context, we will assess the impact of hydrological conditions, atmospheric precipitation and blasting on the landslide processes. Subterranean waters combined with the outlined above conditions determine a wide range of landslide events to a larger extent:  the position of soft clayed semi-rocks below the subterranean water level, weakening effect of hydrostatic weighing due to confined groundwaters and regional tectonic disturbance create conditions for the occurrence of big regional landslides with sliding step-like or curvilinear surface below the Angara water line. These landslide movements, combined with smaller surface shifting, are the greatest threat for the right-bank dam abutment of the Boguchanskaya HPP, which is practically the only energy hub for all regional production chains. In this regard, we will point out that for nine years of geodetic observations we registered slight displacements of the existing geodetic reference network points, which confirms the existence of regional blocs.
Atmospheric precipitation and moisture of the weakened layers have little impact upon landslide behavior in the aeration area, but they form a perched water table, which influences the development of local landslides.
Blasting operations on the right bank of the Angara influence landslide development through subterranean waters.
Direct impact is exercised through primary and cross waves in rocks, the degree of their impact depends on the explosion capacity and the distance from the weakened areas. Blasting impact occurs through fluctuations of ground or pressure waters of the piezometric level when primary waves influence upon water-bearing rocks. At the same time, the impact of an-thropogenic groundwater fluctuations upon landslide behavior is insignificant if compared with natural ones, since their impact time and amplitude are relatively small.
Changing the piezometric level, which characterizes the varying hydrostatic weighing of rocks (interchange influence of overlying rocks upon a weakened layer), aggravates the impact of confined groundwaters towards the development of regional landslides. To a certain extent, this circumstance can be regarded as the occurrence of induced seismicity.
The above geotechnical assessment of the impact of prevailing conditions and main influencing factors upon landslide development reveals some obstacles preventing univocal interpretation of the observed landslide processes using geodetic methods of observation.

Dynamic mathematical model
Dynamic mathematical models describe deformations taking into consideration joint effect of time and main influencing factors. Predicted results of mathematical models can become one of the key tools used as a "solver" in spatial expert systems when shaping knowledge about the state of a link of a particular regional or interregional production chain. In this regard, dynamic models are more advantageous since they have flexible structure that matches the physical essence of the developed process and takes into account inertial interaction both landslide process and influencing factors and their temporal changes (Khoroshilov 2018). For instance, a second-order lag input-output model that describes the movement of a randomly selected landslide point caused by two main influencing factors is defined by the following recurrent expression: where xk is an output variable, i.e., the value of a landslide point shift at the k th quantization interval; Zk, Vk -input expressed as values of two main influencing factors (changes in the blasting charge mass and the volume of   2 , , , 1 2 1 2 /1 2 , for which conditional mathematical expectation is presented as  ˆ/ , At the same time, parameters ˆ, , , 1 2 1 2     are estimated by solving the relevant system of normal equations. At the final stage, the functionality expressed through the conditional correlation function is minimized. To define the order of an auto-regression model and execute the final estimation stage, asymptotically unbiased estimates of the correlation function are calculated on residual modeling errors with using the expression (Box et al. 2015): where time shift m = 0, 1, 2… M < N.
Residual errors εk can be found as difference between actual displacements and their relevant approximated values calculated on the models: The choice of an auto-regression order is based on the correlation function graphs. In our case, the following model can be used (Box et al., 2015): where: µ, η are estimated parameters.

Analysis of landslide geodetic observation
In most cases, data were processed by applying methods of mathematical statistics (correlation analysis). Statistically homogenous groups of displacements were removed from the data array. To assess the homogeneity of landslide dis- Analysis of geodynamic processes monitored with multiple cycle observations over the displacement of local geodetic network points must be accompanied by assuring the stability of particular network points. However, it is extremely difficult and often impossible to provide absolute point stability within local geodetic networks; at the same time, mutual point stability can be provided. In this case, definition of displacement values and their directions lose their initial meaning making the vector representation of displacements useless for analysis and interpretation, since transition from one conditional network to another entails the reorientation of vectors.
At the same time, observation inter-cycle time intervals differed significantly and did not follow the extremums of anthropogenic impacts. Even in conditions of such incomplete information, negative and positive vertical displacements of reference points were found out (Khoroshilov et al. 2013). It is obvious that natural downward sliding, considerably increased by blasting, caused negative displacements. In its turn, soil softening after the removing and following discharge of so-called overburden pressure caused positive displacements.
At the initial stage, the positive and negative trends of temporal developments upon the impact of anthropogenic factors were approximated with 4 th and 5 th degree polynoms having sufficiently tight dependence. Such approximation was only a mathematical description of the observed process that did not synthesize overall trends of its development sufficiently. Therefore, the trends can be used only for a rough forecast of a landslide slope dynamics, since it is impossible to identify place, time and impact level with using the graph of the observed response.
To confirm the discovered specifics of the response of the landslide slope to blasting and unloading, we adjusted the heights of all reference points of the analyzed landslide network using the method offered by Fedoseev (Klyushin et al., 1993). Under the conditions of mutual instability of shifting points, this enabled to define relative displacements of reference points and their actual heights. By using adjusted heights, we calculated the speeds of reference point inter-cycle displacements and estimated their accuracy. In such a way, all monitored parameters of the landslide process were normalized, which enabled to carry out the comparative analysis of its development regardless the duration of the inter-cycle intervals.
The regularities of cinematic development of the isolated centered process were defined by using the correlation theory of random functions (random parameter distribution law). First, we estimated and ensured the process normality for each section. Zero mathematical expectation of the centered process trend determined its linearity at all modeled sections.
As a result, the temporal development of standard and normalized auto-correlation function was approximated. Later, the identified trends allowed assessing the expected landslide dynamics in similar conditions.

Spectral assessment of the effect of influencing factors
To estimate the impact of anthropogenic and some meteorological factors (having harmonic temporal development) upon displacement of landslide reference points, mathematical tools for spectral analysis of time series were used. The es-sence of spectral analysis was to represent a set of operations enabling to present periodic function Х(t) as a sum of harmonic components and define their parameters. According to Jenkins and Watts (1969) the time sequence was represented as series: where fk=K/T -frequency; Тharmonic process period; К=1, 2, 3…N; ak, bk -actual and false spectrum components that can be calculated as: here N is the number of samples.
The amplitude characteristic was considered as Amplitude-frequency dependence can be represented graphically, where harmonic frequencies are plotted along one axis, and their amplitude characteristics are plotted along the other axis. Analysis of a harmonic spectrogram, by breaking it into several characteristic zones, makes it possible to find out the frequency of each harmonic, its period and amplitude. By restoring the analyzed function with inverse Fourier transform in the low-frequency vibration spectra, it is possible to identify the main trend of the analyzed process (i.e., function inclination to time axis), and from calculated characteristics to carry out the analysis of landslide slope deformation process described with a dynamic model. At this stage, the problems related to finding out the impact level of the influencing factors upon the value of displacements of landslide points are solved.
The following analysis procedure was used in the study:  calculation of Fourier transform of the analyzed function;  graphical representation of the amplitude spectrum (spectrogram);  breaking down the spectrogram into characteristic zones;  calculation of the period of each harmonics (based on its frequency);  performing inverse Fourier transform in the low-frequency vibration zone to get the trends of the studied process;  identification of correlating processes by comparing their spectrograms. enous groups of landslide points (8021, 8022, 8023, 8024, 8025 -group one; 8061, 8062, 8063 -group two) as well as relevant influencing factors (water level in the Angara, monthly mean temperature, average daily atmospheric precipitation, amount of removed soil, amount of blasting charge). We also considered the spectra of other meteorological factors (changes in the level of groundwaters, changes in atmospheric pressure) and spectra of horizontal displacements of landslide reference points. The obtained spectrograms are shown in Fig. 2-4. a b Figure 2. Vertical displacement spectrograms of (a) group 1 landslide points; (b) group 2 landslide points a b Figure 3. Vertical displacements triggered by the (a) volume of removed soil; (b) amount of blasting charge a b Figure 4. Vertical displacements caused by changes of (a) air temperature; (b) the Angara water level

Landslide slope movement: separate modeling
The technique for the assessment of landslide slope dynamics consists of two parts (Khoroshilov et al. 2013). The first was aimed at finding invariants of the most stable reference points within the group of landslide points; the second, based on the results of geodetic observations, was focused on landslide susceptibility assessment. Finding the invariants of the most stable landslide points increased the credibility of an actual picture of a landslide development because of the lack of geodetic information.
To study mutual stability of the reference points we applied the method developed by Yu. Fedoseev (Klyushin et al., 1993) where using inter-cycle fluctuations of heights of the observed group of landslide points in the system of their average height enables to not only identify unstable reference points but also determine the size of their displacement. The The next stage was the transition from the values of displacements observed at different times in inter-cycle intervals to the displacement speed, relying on the adjusted heights of the landslide points. Thus, the monitored parameters of the landslide process were normalized regardless the inter-cycle time intervals. A comparative analysis of the obtained results showed that displacement speed values of only 24 of 120 cases were less than their calculation error. The proximity of the displacement speed values and their errors were caused by low speeds at short inter-cycle time intervals and insufficient accuracy of the third class leveling used for geodetic observations. At the same time, the fact that the values of positive displacement speed decreased due to the natural downward landslide movement during unloading of the slope, while negative speed decreased under the inertial impact of the preceding unloading of the slope was also taken into account. This proved the opposite-sign response of the landslide slope to blasting and unloading impacts, which formed the grounds for separate consideration of negative and positive elements of vertical displacements at the next stage of the study in accordance with the determination and probabilistic nature of the landslide process. Fig. 6 shows the development of negative and positive trends of temporal landslide displacements. Figure 6. Results of separate modeling of negative and positive components of the vertical displacement process Then, separate centering of the process components was carried out; this allowed revealing the nature of two types of anthropogenic impacts as temporal changes of the average displacement values fixed by observations (Fig. 7, 8).   The follow-up procedure was the modeling of five centered dynamic implementations of the landslide process represented as time approximated mathematical expectation.

Modeling anthropogenic displacements (prediction mathematical modeling)
Responses of the landslide slope to anthropogenic impacts were analyzed with taking into account nine- Within the study, the following input impacts were used to build up prediction mathematical models: Zk -change in charge mass when carrying out blasting operations throughout the entire cycle of geodetic observations; Vk-change in the volume of removed soil after blasting operations throughout all cycles of geodetic observations.
As output factor xk, the displacement values of landslide points at k th sampling step were used (adjustment by Yu. Fedoseev).
The overall period for prediction models was 9 years; the observation sampling interval was equal to 0.2 years.
Below is a group of prediction mathematical models with two input impacts x(Zk, Vk) observed every 0.
Tab. 1 gives the values of root-mean-square errors m, calculated upon discrepancies between adjusted and approximated displacements of landslide points.

Use of the results of mathematical modeling in a spatial knowledge-based system
Presently, the development of economic and financial situation depends on the use of geospatial data. One of the main goals of infrastructural organization of geospatial data is to provide users to be capable of acquiring complete, exact and updated dataset at the right time. With the advance of information technology, the capabilities of data storage, data mining, and data productivity have triggered a rise in the use of spatial information, especially in the area of regional social and economic development, management and safety (Ahsan et al. 2016;Chan et al., 2018). Efficient decision making at the level of regional or interregional spatial development and management relies on the collected information (format, source, attributes etc.) and the knowledge of decision makers. To assist the decision makers in using, selecting, and processing in- formation, spatial knowledge base and knowledge-based systems are developed (Harmon and King 1985;Cockcroft 2004). Fig. 10 represents an example of key elements of a knowledge-based system and their interdependency that enable to optimize workflows and provide safety and sustainability of an economically important link of a production chain (LPC), which may be unique, especially in such regions as Siberia and the Far East, supporting the operation of a number of production chains and industries. Frank et al. (2007) note that the keystone elements of the system are to be regularly updated with spatial data (observations), which are used in decision-making and management. Mathematical models describing, simulating, and predicting the LPC state, its significant facilities and their environment, and "solvers" that automatically trigger alerts, lists and protocols of needed actions and warnings when something may get or is getting out of control, either endangers operation of the LPC, allow specialists and decision makers to undertake preventive actions to keep it in accident-free operation.  Nguyen et al. 2020). In our case, the Rela-Ops model offered by Nguyen et al. (2020) seems to be most appropriate for combining all kinds of spatial data into knowledge. The model is based on not only object-oriented approaches but also ontology, and includes foundation components consisting of concepts, relations, operators, and inference rules. Besides the structure, each concept of this model is a class of objects, which solve problems on their own. Its processing of algorithms for solving problems combines the knowledge of relations and operators in the reasoning. What is more a knowledge model for multiple knowledge domains, in which each sub-domain has the form as the Rela-Ops model, may be built increasing ability and efficiency of the automated monitoring and maintaining of a complicated network of dependent LPCs.
During the study, modified Rela-Ops model and its algorithms were tested on the experimental data, which proved their validity for building reliable intelligent spatial problem solvers providing from 95.2 to 96.4% of correct management decisions.

Discussion
Analyzing the spectrograms (Fig. 2 a, b) of vertical displacements of the two groups of landslide points, two main harmonics with frequencies 40 and 75 were identified that match the periods of 12 and 3 months. The spectrograms of influencing factors enabled to find out bursts at frequency 40 that is associated with natural impacts and matches the period of 12 months. It is obvious, that this harmonics is caused by seasonal fluctuations. The spectrogram of vertical displacements of the first group of landslide points demonstrates a sharp peak in the vicinity of zero frequency (Fig. 2 a) and intensive low-frequency fluctuations. When restoring this part of the spectrum, a straight line with a slight time axis inclination was obtained, which can be interpreted as the development of vertical displacements of landslide points.
The spectrogram of the second group of landslide points (Fig. 2 b) has a gentler peak in the area of zero frequency and more intensive low-frequency fluctuations. When restoring this part of the spectrum, a straight line with a more significant time axis inclination was obtained, which means that the speed of vertical displacements of the second group of landslide points is higher than the first one.
The spectrograms on changing the amount of soil taken out during slope unloading and changing the charge value in the course of blasting operations (Fig. 3 b) show a peak in the area of frequency 75 that matches the period of 8 months, and that both spectrograms are practically identical. The same frequency was observed in the spectrograms of vertical displacements of landslide points, which makes it possible to conclude that blasting operations and slope unloading are the major influencing factors responsible for vertical displacements of landslide points. It should be noted that apart from the main seasonal frequency 40, the spectrogram of changes in the Angara water level also has a peak at frequency 75, which means that blasting operations influence the changes in the Angara water level (Fig. 4 b). This is additionally confirmed by the fact that spectrograms of other meteorological factors do not have peaks at this frequency.
Then the spectrograms of horizontal displacements of both groups of landslide points were analyzed. The obtained results demonstrated that they correlate with vertical displacements and are influenced by the same factors as in case of vertical displacements. The described periodicity of the influencing factors upon displacements of the analyzed landslide points was revealed across the whole area of the right-bank landslide slope.
To specify the periodicity of these impacts, we analyzed the nature of blasting operations and slope unloading impacts. To this purpose we identified time dependencies on changing the charge in the course of blasting operations; chang-ing the amount of removed soil in the course of slope unloading in the landslide area; changing the amount of bulk soil in the area of rocky earth deposits. These areas were chosen because landslide slope unloading and filling in rocky earth deposits took place in the close vicinity to the landslide points used for geodetic observations (8061-8063) and the fact that these observations were carried out for 9 years. The analysis of the above processes showed that earth works on the slope unloading were carried out every two-three years with following decrease of their intensity, which gave a broader and fuller picture of the processes in the site.
The heights of all landslide points in the analyzed group were adjusted according to the method offered by Yu. Fedoseev, i.e., estimates of their actual heights were obtained by separate modeling of negative and positive displacement trends reflecting the response of the landslide slope to blasts and unloading (Fig. 5). It was found that in 21 of 120 cases, the speed values exceed their calculation errors significantly. The closeness of the speed values and their errors is due to low speeds that take place at short inter-cycle time intervals and insufficient accuracy of geodetic observations (III class leveling).
It should be taken into account that estimated values of positive speeds were slightly reduced because of landslide slope unloading during its natural downward sliding, while negative speeds decreased under the inertial impact of the preceding unloading. Thus, different response of landslide slope to blasting and unloading impacts was proven. Later this formed the grounds for considering negative and positive components of vertical displacement process separately according to their deterministic and probabilistic nature. The above process components were centered separately, which revealed the nature of the two types of anthropogenic impacts as temporal changes of average displacement values, registered with geodetic observations. The centered values of the negative component of the process had, mainly, a positive development trend, while the positive component showed a negative development trend, confirming the above-described compensatory and inertial impact of the influencing factors.
The centered values of both components of the process mainly reflect natural character of its development and are slightly prone to anthropogenic impacts. Therefore, they were close to each other, which allowed us to combine them by averaging on closely adjacent sections, including centered values of some remote sections. Five integrated implementations of the centered process were obtained as well as graphs of the development of average negative and positive vertical displacements characterizing the degree of influence of each of the considered anthropogenic impacts.
As a result of the studies carried out, analytical validation and experimental assessment of the technique describing the process of landslide slope movement represented by two separate components related to blasting operations and slope unloading due to removal of soil from rock quarries was done. The technique of transition to centered normally distributed process of vertical landslide displacements with zero mathematical expectation, independent from blasting and unloading impacts, and the procedure of landslide process dynamic modeling with inverse verification of reference and model values of centered displacements were also developed.
At the final stage, dynamic mathematical modeling was carried out according to recurrent expression (1). The modeling matches the physical essence of the process development by taking into account the inertial nature of displacements monitored by geodetic observations as well as influencing factors and their temporal changes. To describe the dynamics of the analyzed process, we used an "input-output" second-order inertia model, representing the displacement of a randomly selected landslide point under the influence of two main anthropogenic factors.
Dynamic mathematical modeling of a landslide processes (prediction models) was used in algorithms forming a "solver" (intelligent problems solver) of a spatial expert system and later tested with a spatial knowledge base of a modified Rela-Ops model. Managing decisions made by the spatial knowledge base of the modified Rela-Ops model over 251 situations simulated for 5 spatial classes proved to be highly efficient showing their reliability in 239 cases.

Conclusions
Summarizing the results of the studies on dynamic mathematical modeling of landslide processes and the use of the built prediction models as an intelligent problem solver of a spatial expert system, we shall underline the following:  detailed representation and analysis of spatial information depends on the location and density of landslide reference points and time frame of geodetic observations over their displacements;  to carry out credible and high-quality mathematical modeling it is necessary to have complete quantitative data on the development of major factors influencing landslide activity, and zones of their influence;  partial loss and sparseness of landslide reference points hampered reconstruction of an unbiased spatiotemporal picture of landslide movements, because of that we cannot clime the completeness and detail analysis of the studied processes;  correlation of vertical and horizontal displacements confirms the overall landslide behavior revealed by geodetic observations; the analysis showed regularities of the impact of influencing factors upon the displacements, which enabled to describe them at a new level;  correct selection of the input influencing factors and the type of a mathematical model allowed using built-up mathematical models to forecast the future displacements of observed landslide points when studying the landslide slope dynamics; the forecast RMSE for the displacement of different landslide points varied from 0.37 to 0.68 mm, which indicates the high accuracy of the forecast;  presented technique for identifying features and regularities of landslide behavior was used as an intelligent problems solver of a spatial expert system, based on the modified Rela-Op model under conditions of incomplete information. Verification of the methodology on simulated situations showed that it can be efficiently used as a "solver" of a spatial expert system when dealing with similar and other types of slope impacts supported by sufficient geodetic observations.