Slope digital twin modeling method based on multi-source heterogeneous data fusion

Abstract

The invention provides a multi-source heterogeneous data fusion side slope digital twin modeling method, which comprises the following steps of: acquiring side slope multi-dimensional monitoring data by arranging a GNSS (Global Navigation Satellite System) sensor, a multi-point displacement meter, a distributed optical fiber strain sensor, an accelerometer, an osmometer, a monocular camera and satellite remote sensing image equipment; the collected data is converted into a unified format through time alignment, space registration and standardization processing and serves as modeling input; the method comprises the following steps: constructing an initial digital twinborn model reflecting the real form and physical characteristics of a slope by utilizing a three-dimensional modeling and finite element simulation technology; in combination with real-time sensing data, model evolution is dynamically driven based on a space-time fusion algorithm, boundary conditions and material parameters are automatically corrected through actual measurement deviation feedback, and continuous twin iteration updating of the model is achieved; and finally, extracting a landslide risk index to realize real-time early warning of the side slope. According to the invention, multi-source sensing and digital twinborn fusion is realized, and the accuracy, real-time performance and intelligent level of slope monitoring are improved.

Description

Technical Field

[0001] This invention relates to the field of geological disaster monitoring and modeling technology, specifically to a method for modeling digital twins of slopes by fusing multi-source heterogeneous data. Background Technology

[0002] Landslides are a widespread, sudden, and highly destructive geological hazard. Their occurrence not only damages the landform and disrupts the ecological environment, but also poses a significant threat to the safety of engineering facilities and the lives and property of people. In recent years, with the advancement of infrastructure construction in mountainous areas, the frequency of landslide disasters has been continuously increasing, bringing severe challenges to infrastructure operation, mountain transportation, and geological safety management. To reduce the risk losses caused by landslides and improve early warning and response capabilities, research on precise and efficient dynamic monitoring and modeling technologies for landslide bodies is particularly crucial.

[0003] Currently, landslide monitoring technologies mainly include various sensors such as satellite remote sensing, unmanned aerial vehicles (UAVs), GNSS (Global Navigation Satellite System), slope radar, inclinometers, strain gauges, multi-point displacement gauges, and pore water pressure gauges. Jiang et al. (2022) used UAV ground laser scanning technology to acquire three-dimensional topographic data of slopes, which can overcome terrain limitations to obtain surface deformation data of slopes in blind areas. Ouellet et al. (2024) used distributed optical fiber equipment to acquire continuous strain data of slow creep slopes, realizing long-term continuous strain monitoring along the failure direction. Chen Mohan (2025) et al. used displacement sensors to capture the rate of slope displacement change and set a sliding early warning threshold based on the displacement change data to realize early warning of slope instability. However, single-source monitoring sensors have limitations in sampling frequency, spatial resolution, or observed physical quantities, making it difficult to comprehensively reflect the complex evolution law of landslide bodies. Meanwhile, different types of sensors have different sampling frequencies, monitoring ranges, unit dimensions, and signal-to-noise characteristics, lacking a unified expression, which directly affects the accuracy of multi-field information fusion modeling of slopes.

[0004] For modeling and simulating the evolution of landslide disasters, numerical simulation methods such as the finite element method, discrete element method, and particle method are currently available. These simulation methods can reproduce the landslide instability mechanism relatively well under ideal initial conditions and known material parameters. However, in actual engineering, the complexity of the internal structure of the landslide body, the time-varying nature of material parameters, and the dynamic changes in boundary conditions are difficult to accurately obtain. This results in traditional numerical models lacking dynamic coupling with the actual situation, leading to large prediction errors. In addition, most existing numerical simulation models are "open-loop" calculation processes, where slope model parameters and boundary conditions are set once at the beginning of modeling and remain fixed during the analysis process, making it impossible to dynamically correct the parameter model based on monitoring data. This method is difficult to reflect the real dynamic changes of the landslide body, limiting the model's adaptability and timeliness, and making it difficult to form a slope early warning system with adaptive and real-time evolution capabilities.

[0005] Due to the limitations of traditional technologies in real-time perception, dynamic simulation, and closed-loop early warning, digital twin technology, as a new concept in the industrial field, provides a new approach to problem-solving. A digital twin is a technology that enables real-time interaction and collaboration between the physical and virtual worlds by constructing a virtual copy of a physical entity. Introducing this concept into landslide monitoring aims to build a spatiotemporal digital mapping system for slopes, achieving continuous state perception, dynamic deduction, risk assessment, and early warning. Currently, landslide digital twin modeling technology is still in its initial exploratory stage. Xue et al. (2025) constructed an AI-assisted digital twin framework for highway slope stability, using lidar technology to generate a real-time digital elevation model of the slope and using neural networks to deduce pore pressure and slope strength variables to identify landslide risks. Ju et al. (2025) constructed a digital twin simulation model based on real-time rainfall data to study rainfall-induced landslides, and assessed landslide risk coefficients in real-time within a digital twin environment. Yang et al. (2024) constructed a multi-temporal digital twin of the Baige landslide on the Jinsha River using 10 UAV aerial surveys over 5 months, analyzed the spatiotemporal evolution of the landslide, and achieved quantitative analysis of landslide deformation area and collapse volume. However, the digital twin constructed based on single-source data differs significantly from the actual situation on-site, and there is still a lack of a multi-state variable coupling evolution mechanism and a comprehensive indicator system for landslide risk identification.

[0006] In summary, current research on landslide digital twins is still in its early exploratory stage, with limited research on technologies such as multi-source data fusion, model parameter self-evolution, and multi-variable coupled early warning. Therefore, there is an urgent need to propose a novel digital twin method for slopes that integrates multi-source heterogeneous monitoring data, achieves closed-loop dynamic correction of the three-dimensional finite element model, and constructs multi-physical variable coupled early warning indicators, thereby enabling high-precision real-time modeling and advanced early warning of landslide disasters. Summary of the Invention

[0007] To address the technical problems existing in the prior art, this invention provides a method for modeling digital twins of slopes using multi-source heterogeneous data fusion for slope safety monitoring and early warning. Specifically, it provides a method for multi-source heterogeneous monitoring data fusion and landslide risk identification of slope digital twins. First, a multi-source sensor network, including GNSS, multi-point displacement gauges, accelerometers, piezometers, and satellite remote sensing, is deployed on the target slope to collect data such as three-dimensional displacement, deep displacement, structural vibration, pore pressure, and three-dimensional point cloud data of the slope. The sensor data is then spatiotemporally aligned by unifying the time axis and spatial coordinate system. Next, the Z-score algorithm is used to remove abrupt data changes, ultimately outputting structured time-series data. Finally, based on the cleaned and aligned data, acceleration responses such as slope displacement, strain, and pore water pressure are extracted. Based on the state characteristic values, a nonlinear evolution model considering displacement-strain geometry and pore pressure diffusion is established, and the state vector characteristic values ​​are updated using a Kalman filter method combined with a multi-field spatiotemporal evolution model. Furthermore, a three-dimensional finite element model is constructed using image point cloud data, and ABAQUS finite element numerical simulation is performed to obtain the modeling calculation results. Subsequently, by defining displacement, strain, and pore pressure residual thresholds to trigger a model correction mechanism, an inversion function with residual minimization as the objective is constructed to dynamically correct the modeling material parameters of shear modulus, friction angle, and cohesion, achieving twin iterative updates of the finite element model. Finally, the dynamic state vector is extracted, and a weighted landslide hazard index function is constructed to achieve landslide risk level identification and early warning. Periodic closed-loop monitoring and early warning of slopes are completed through new data input. This technology integrates multi-source heterogeneous monitoring information, extracts state characteristic vectors from multiple fields such as displacement, strain, and seepage pressure, realizes iterative updates of twin models based on residual triggering and parameter inversion, and establishes a comprehensive landslide risk identification index. It features high modeling accuracy, strong model self-updating capability, good early warning timeliness, and applicability to complex slope engineering.

[0008] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0009] A method for modeling digital twins of slopes using multi-source heterogeneous data fusion includes the following steps:

[0010] S1. Multi-source sensor deployment and data acquisition: GNSS sensors, multi-point displacement gauges, distributed fiber optic strain sensors, accelerometers, piezometers, monocular cameras and satellite remote sensing imagery equipment are deployed in the target slope area to acquire three-dimensional displacement, deep displacement, strain change, vibration response, pore pressure state, image information and slope three-dimensional topography data, respectively.

[0011] S2. Data Synchronization and Cleaning: The raw sensor data obtained in step S1 is standardized and preprocessed, specifically including the following steps:

[0012] S21. Time Alignment: By constructing a unified time axis, data with different sampling frequencies are time aligned. The sliding window method is used to extract feature values ​​from high-frequency data, and low-frequency data is aligned through linear interpolation.

[0013] S22. Spatial Registration: Register the spatial position of the sensor to the slope engineering coordinate system according to the sensor type;

[0014] S23. Anomaly Removal: Anomaly detection and cleaning are performed on outliers, drift interference, and duplicate records in the original sensor data. Specifically, for continuous signals of displacement, strain, and tilt angle, the Z-score algorithm combined with a sliding time window is used to remove abrupt values ​​based on three times the standard deviation above and below the mean. For pore water pressure and rainfall monitoring data, physically reasonable interval thresholds are set for interval filtering. For high-frequency acceleration signals, a low-pass filtering process is preset to suppress spike noise, ultimately retaining physically reasonable and statistically stable data samples.

[0015] S24. Unit conversion and format unification: Standardize the units of physical quantities and output them in a structured manner as a time series data table with a unified field format;

[0016] S3. State Variable Extraction: Based on the cleaned and aligned structured multi-source monitoring data from step S2, key physical quantities of slope state evolution are extracted to construct a set of state variables for model identification, stability assessment, and evolution simulation. This includes the following steps:

[0017] S31. Three-dimensional displacement field reconstruction: Based on discrete displacement observations from GNSS monitoring points, multi-point displacement gauges, and image feature points, a spatial correlation function between each observation point is first constructed. Then, the displacement of grid points without sensors in the slope area is estimated using the Kriging interpolation method to form a spatially continuous initial displacement field. Subsequently, a displacement evolution state transition equation is established on a unified time axis. The spatial interpolation results of the Kriging method are used as the observation input, and the displacement estimate of each new node is recursively calculated using the Kalman filtering method. Finally, the slope mechanical equilibrium equation is introduced as a soft constraint to obtain three-dimensional displacement field data that satisfies physical consistency and spatiotemporal continuity.

[0018] S32. Strain Tensor Field Inversion: Based on distributed fiber strain data, the parameters are determined by using the Tikhonov regularization method and combining it with the L-curve criterion to invert the strain tensor field inside the slope, obtain a continuous and smooth strain distribution, and evaluate its uncertainty.

[0019] S33. Pore pressure state calculation: Based on the monitoring data of the piezometer, the pore water pressure is converted into groundwater head and effective stress to form groundwater state variables, and its spatial distribution is determined in combination with the deployment depth.

[0020] S34. Acceleration Response Extraction: The acceleration signal is filtered to extract the effective response components, and the peak value, energy, and frequency index are calculated to characterize the structural disturbance and vibration response features.

[0021] S4. Multi-source heterogeneous data fusion: Based on the state variables extracted in step S3, state information from different types of sensors is fused to construct a unified state representation for physical modeling and stability prediction, achieving spatiotemporal joint estimation. This includes the following steps:

[0022] S41. Constructing a state vector: Constructing a state vector composed of three-dimensional displacement, principal strain components, pore pressure, and acceleration eigenvalues;

[0023] S42. State-space modeling: The observation equation Z based on the state vector. t =h(x t )+v t and state prediction equation In the observation equation, Z t Let represent the observation data from multiple sensors at time t, h(·) represent the nonlinear mapping relationship between various sensors and the state, and x t The state vector at time t represents the three-dimensional displacement, principal strain components, pore water pressure, and acceleration eigenvalues, v. t For the observation noise term; in the state prediction equation, f(·) represents the nonlinear evolution function considering displacement-strain geometry, pore water pressure diffusion, and quasi-static or dynamic evolution processes. Represents the state vector x based on the previous time step. t-1 The predicted state of the slope, w t This is the state prediction noise term;

[0024] S43, Observation-based update: Where K t Z is the Kalman gain coefficient. t The observed values ​​are used to achieve spatiotemporal joint estimation of multi-source states;

[0025] S5. Construction of the 3D finite element model: based on the update in step S4 Multi-source heterogeneous data are then used to construct a three-dimensional slope geometric model using digital elevation model or point cloud data, and a finite element model containing geological units, physical parameters and support structures is established, with mesh refinement implemented in the potential sliding zone.

[0026] S6. Simulation-driven and initial solution: Apply boundary loads and pore pressure conditions to the established finite element model, carry out nonlinear static or dynamic analysis, and obtain the predicted displacement and stress response distribution.

[0027] S7. Model feedback correction mechanism:

[0028] S71. State error triggering condition: Define the residual threshold ε of multidimensional state variables such as displacement, strain, and pore pressure, and compare the mean square error (MSE) of the simulation results with the measured state monitoring parameters in real time. When MSE > ε, the model correction mechanism is triggered.

[0029] S72. Parameter Inversion and Model Update: Construct a nonlinear least squares objective function with the goal of minimizing the residuals. The optimization function is F = min θ ||Y sim (θ)-Y obs || 2 +λ||θ-θ prior || 2 Where θ is the model parameter and θ={G,φ,c}, representing the numerical simulation parameters of shear modulus, friction angle and cohesion, respectively; Y sim (θ) represents the simulation response, Y obs Let λ be the observed state variable, λ be the regularization factor, and θ be the θ value. prior The parameters are prior estimates, i.e., empirical values ​​selected based on field conditions; the optimization results are used to correct material parameters and pore water distribution field, and then reapplied to the finite element model to achieve dynamic closed-loop iterative updates.

[0030] S8. Construction of Landslide Risk Index:

[0031] S81. Risk Variable Extraction: Extract the characteristic quantities, including the derivative, increment and volatility of multi-source state variables over time, and construct a dynamic indicator set for landslide trend identification and risk assessment.

[0032] S82. Risk Formula Expression: Based on the comprehensive simulation of variables including safety factor, acceleration response, pore pressure growth rate, and strain increment, a weighted landslide hazard index function is constructed to determine the landslide evolution stage and risk level range.

[0033] S9. Landslide Risk Identification and Intelligent Early Warning: When the risk indicators exceed the threshold, an early warning is automatically triggered, and the pre-landslide area and risk level are output.

[0034] S10, Model Closure and Periodic Update: Newly collected data is input into step S5 to complete the state update and risk reassessment, forming a closed-loop iteration.

[0035] Furthermore, in step S1, the GNSS sensor obtains the displacement using the three-dimensional coordinate difference method as follows:

[0036]

[0037] Where Δx, Δy, and Δz are the displacement changes in the x, y, and z directions, respectively;

[0038] Displacement of discrete points measured by multi-point displacement gauge {u i The continuous displacement field is calculated using an interpolation function as follows:

[0039]

[0040] Where, φ i (z) is the interpolation shape function corresponding to the i-th measuring point, representing the displacement at position z caused by the displacement u at the i-th measuring point. i The contributed weighting factors are determined by the selected interpolation algorithm, including linear, spline, or finite element shape functions, and satisfy the interpolation condition φ. i (z j )=δ ij And the constraint that the sum of the weights is 1;

[0041] A monocular camera continuously acquires a sequence of images of the slope surface under fixed installation conditions. Pixel displacement is calculated using the optical flow or feature point matching results of images from adjacent time points. Let I... t (x,y) and I t-1 (x, y) represent the grayscale images acquired at times t and t-1, respectively. Assuming that the grayscale remains constant over time, the optical flow constraint equation is as follows:

[0042] I x ·u+I y ·v+I t =0

[0043] Among them, I x I y I t These are the gradients of the image in the x, y and time directions, respectively, and (u, v) are the planar displacement components of the pixel.

[0044] By solving for (u,v) through feature matching or dense optical flow, and combining the camera imaging model with ground calibration, pixel displacement is converted into physical displacement:

[0045]

[0046] Among them, s x s y Given the pixel size and Z as the depth distance of the corresponding point, the surface displacement vector field Δr(x,y)=(ΔX,ΔY) of the monitoring area can be obtained.

[0047] Furthermore, the strain ε acquired by the distributed fiber optic strain sensor in step S1 t Calculated by the following formula:

[0048]

[0049] Where ΔL(s,t) represents the displacement increment of the monitoring segment, L0 is the initial length, and s is the position variable along the fiber optic axis.

[0050] Furthermore, the pore water pressure p measured by the piezometer in step S1 t The following relationship must be satisfied:

[0051] p t =ρ w ·g·h t

[0052] Where, ρ w Let g be the density of water, g be the acceleration due to gravity, and h be the acceleration due to gravity. t This refers to the groundwater head height.

[0053] Furthermore, the acceleration a obtained by the accelerometer in step S1 t The result is obtained by the following formula after filtering:

[0054] a t =LPF(a raw (t))

[0055] Where LPF represents the low-pass filter function, a raw (t) represents the original sampled signal.

[0056] Furthermore, the state variable extraction method in step S3 includes the following modeling and reconstruction mechanism:

[0057] (1) Construct a structured multi-source monitoring dataset and integrate GNSS displacement data u GNSS (t), satellite imagery data u SID (t), Image recognition feature point displacement u IMG (t), fiber strain data ε DFOS The monitoring state vector is constructed by aligning the pore pressure data p(t) from the piezometer and the accelerometer data a(t) with the spatial location (x, y, z) according to a unified time series t:

[0058] x(x,y,z,t)=[u(x,y,z,t),ε(x,y,z,t),p(x,y,z,t),φ(t)]

[0059] Wherein, φ(t) represents the external excitation acting on the slope system at time t. Its value can include rainfall time series, ground motion acceleration, construction load changes, etc., to reflect the influence of the external environment on the evolution of the slope state.

[0060] (2) Based on a three-dimensional spatial grid, the continuous three-dimensional displacement field u is reconstructed using the Kriging interpolation method. 3D(x,y,z,t) is used, and the time series is dynamically corrected using the Kalman filter algorithm to improve its continuity and consistency in the spatiotemporal dimensions, satisfying the following optimization objectives:

[0061]

[0062] Among them, w i ·||u i -u(x p ,y p ,z p ,t i ) 2 The objective is to minimize the residuals of the observed data. To optimize the obtained three-dimensional displacement field estimation vector, which contains three components: x, y, and z, u(x p ,y p ,z p ,t i ) represents the spatial location (x) of the observation point. p ,y p ,z p 0 and time t i The displacement vector is obtained by interpolation or solving for the continuous displacement field to be determined; U obs ={u GNSS ,u SID ,u IMG} represents the fused multi-source displacement observation data set, w i The weighting coefficient for the i-th observation sensor reflects the data confidence level of different sensors; u i This represents the displacement vector measured by the i-th observation sensor; λ1 is the soft-wire term in the slope mechanics equilibrium equation, and λ1 is the regularization coefficient used to weight the equilibrium observation fitting term and the physical constraint term. ρ represents the divergence of the stress tensor σ corresponding to the displacement field u, and ρ represents the bulk density of the slope material.

[0063] (3) Based on fiber strain observations, the continuous spatial strain tensor field ε(x,y,z,t) is obtained by inversion using the Tikhonov regularization method. The optimized model is as follows:

[0064]

[0065] Where A is the measurement projection matrix, d is the observed strain data, λ2 is the regularization parameter, and L is the regularization operator;

[0066] (4) Convert the pore water pressure data p(t) into head variable. And calculate the effective stress variable σ'=σ-p t These constitute groundwater state variables;

[0067] (5) After the acceleration signal is processed by bandpass or lowpass filtering, the typical characteristic parameters of the disturbance response are extracted as follows:

[0068] Peak acceleration: A peak =max|a(t)|

[0069] Effective value (RMS):

[0070] Energy index: E=∫|a(t)| 2 dt

[0071] Clock speed specifications:

[0072] Where T is the signal analysis duration, i.e., the length of the time interval covered by integration or statistical calculation, and f is the frequency variable used to describe the distribution of the acceleration signal in the frequency domain. This represents the Fourier transform operator, which transforms a time-domain signal to the frequency domain. This represents the complex amplitude of the acceleration signal a(t) at frequency f;

[0073] Ultimately, a multidimensional spatiotemporal state variable set x(x,y,z,t) is formed, which is used to drive the modeling and evolution simulation of the digital twin of the slope.

[0074] Furthermore, the data fusion algorithm in step S4 employs extended Kalman filtering, including:

[0075] Constructing the state vector:

[0076]

[0077] Status information:

[0078]

[0079] P t =F t P t-1 F t +Q t

[0080] Where, x t-1 This represents the system's state vector at the previous time t-1, which serves as the input for predicting the current state; w t For state prediction noise; P t To estimate the state at time t The magnitude of uncertainty and the correlation between data in each dimension; F t Q is the Jacobian matrix for the state transition; t The measured noise covariance of the sensor;

[0081] Observation Update:

[0082]

[0083] Among them, K t H is the Jacobian matrix of the observation model. t Let be the Jacobian matrix of the observation model, representing the sensitivity of the state vector at time t to the linearization of the observation equation. R is the transpose of the Jacobian matrix of the observation model, with dimensions n×m, used to map the observation space error back to the state space. t To observe the noise covariance;

[0084]

[0085] P t =(IK t H t )P t

[0086] Among them, K t Z is the Kalman gain coefficient. t The observed values, where I is the identity matrix.

[0087] Furthermore, in step S5, a three-dimensional slope geometric model is constructed from the digital elevation model or point cloud data. This three-dimensional slope geometric model is reconstructed from a digital elevation model formed by fusing multi-source surveying and mapping data. The data includes: reflectance images acquired from satellite remote sensing images. sat (x,y), image sequence acquired by a monocular slope camera The elevation difference data obtained by the structured optical flow SfM method, and the high-precision GNSS three-dimensional coordinate points p deployed at the slope control points. i =(x i ,y i ,z i The fused elevation field is modeled using the following interpolation expression:

[0088]

[0089] Where N represents the number of high-precision control points provided by the GNSS data source, M represents the number of elevation difference data points obtained by optical flow SfM method inversion from a monocular camera, and L represents the number of elevation points provided by satellite remote sensing image data. Let be the elevation value of the i-th GNSS control point. Let j be the elevation value of the SfM inversion point of the j-th monocular camera. Let ω be the elevation value of the k-th satellite remote sensing elevation point. i η j μ kThese are the normalized weighting coefficients allocated based on the accuracy and spatial coverage density of each data source, satisfying the following normalization condition:

[0090]

[0091] Based on the soil and rock stratification data provided by the engineering survey, the slope area Ω is divided into k geological units Ω. k Each element is assigned the following material parameter vector:

[0092] θ k =(E k ,v k ,c k ,φ k ,ρ k ,k k )

[0093] Among them, E k For the elastic modulus, v k For Poisson's ratio, c k For cohesion, φ k ρ is the internal friction angle. k For density, k k Permeability coefficient;

[0094] If a potential slip surface exists in the geological structure Γ slip Therefore, a zero-thickness contact element is used for modeling. This contact interface follows the Mohr-Coulomb shear failure criterion, expressed as:

[0095] τ≤c slip +σ·tanφ slip

[0096] Where τ is the shear stress, σ is the normal stress, and c slip φ slip These are the cohesion of the slip surface and the internal friction angle, respectively.

[0097] If the slope is supported by anchor bolts, frame beams, and anti-slide piles, and finite element modeling is performed using either rod elements or shell elements respectively, then the overall stiffness matrix can be uniformly represented as:

[0098] K global =K soil +K support +K interface

[0099] Among them, K soil K is the global stiffness matrix of the soil element. support For the overall stiffness matrix of the support structure unit, K interface For interface element stiffness matrix, represent the independent contribution terms of the stiffness of the anchor-soil coupling interface;

[0100] The model region Ω is discretized using tetrahedral elements, and the slip surface Γ slip The region undergoes mesh refinement processing, and the cell scale within the refined region satisfies:

[0101]

[0102] Among them, h global This represents the typical element size in the unrefined region. The refinement process aims to improve the analytical accuracy of the displacement and shear stress field distributions, and to meet the numerical stability and resolution requirements near the slip surface.

[0103] Regarding boundary conditions, a fully fixed constraint is set at the bottom of the slope model, namely μ. x =μ y =μ z =0, μ x μ y μ z These represent the displacement constraint coefficients in the three directions; the sliding surface is configured with contact-friction boundary conditions.

[0104] The final three-dimensional finite element model includes: the fused surface geometry Ω surf Multi-layered geological unit structure Ω k and its parameter distribution θ k Slip surface Γ slip The support structure system and its coupling relationship with the slope soil and rock mass, the non-uniform mesh division scheme and the corresponding boundary conditions provide a physical basis for subsequent numerical simulation evolution and model feedback iteration.

[0105] Furthermore, the model feedback mechanism in step S7 includes the following judgment conditions and parameter update steps:

[0106] If the following error conditions are met:

[0107]

[0108] The inversion process is then initiated, and the following optimization objective is fitted using a nonlinear minimization square method:

[0109]

[0110] in, This represents the simulated displacement of the i-th monitoring point. δ represents the measured displacement of the i-th monitoring point. u δ ε This represents the relative error between displacement and strain. The simulated and measured strains are represented by N, which represents the number of observation data points; θ = {G, φ, c} represent the numerical simulation parameters of shear modulus, friction angle, and cohesion, respectively.

[0111] Furthermore, in step S8, the landslide hazard index K risk Satisfy the following expression:

[0112]

[0113] Where w1, w2, w3, and w4 are the weighting coefficients of each indicator. For displacement rate, For strain rate, The rate of change of pore water pressure. This represents the rate of change of the safety factor.

[0114] Compared with existing technologies, the slope digital twin modeling method based on multi-source heterogeneous data fusion provided by this invention integrates various heterogeneous monitoring data, including those from GNSS sensors, multi-point displacement gauges, distributed fiber optic strain sensors, accelerometers, piezometers, monocular cameras, and satellite remote sensing imagery. Through time alignment, spatial registration, and standardization, a unified state vector is constructed, and key physical field quantities such as three-dimensional displacement field, strain tensor field, and pore pressure response are extracted to establish a digital twin modeling system coupled with multi-field evolution. On this basis, a dynamic correction mechanism for model parameters based on "monitoring-simulation" residual triggering is proposed to achieve adaptive updates of shear modulus, friction angle, and cohesion parameters, thereby improving the model's predictive ability. Furthermore, a weighted landslide hazard index function is constructed, integrating the derivatives and fluctuation characteristics of multi-source state variables to divide risk level intervals, achieving advanced identification and intelligent early warning of landslide risks. This method has the advantages of multiple information fusion types, high model building accuracy, fast twin simulation response, and strong risk identification capability, which significantly improves the real-time performance, accuracy, and engineering adaptability of slope disaster monitoring and early warning, and can be well applied to slope risk management scenarios. Attached Figure Description

[0115] Figure 1 This is a flowchart of the landslide monitoring and early warning system driven by multi-source sensors provided by the present invention.

[0116] Figure 2 This is a time series variation curve of multi-source monitoring data provided by the present invention.

[0117] Figure 3 This is a visualization diagram of the three-dimensional modeling and numerical analysis process of landslide areas provided by the present invention.

[0118] Figure 4 This is a flowchart of a virtual-real coupled landslide risk assessment driven by multi-source information fusion provided by the present invention. Detailed Implementation

[0119] To make the technical means, creative features, objectives and effects of this invention easier to understand, the invention will be further described below with reference to specific illustrations.

[0120] like Figure 1 As shown, this invention provides a method for modeling a digital twin of a slope using multi-source heterogeneous data fusion, comprising the following steps:

[0121] S1. Multi-source sensor deployment and data acquisition: GNSS sensors, multi-point displacement gauges, distributed fiber optic strain sensors, accelerometers, piezometers, monocular cameras and satellite remote sensing imagery equipment are deployed in the target slope area to acquire three-dimensional displacement, deep displacement, strain change, vibration response, pore pressure state, image information and slope three-dimensional topography data, respectively.

[0122] S2. Data Synchronization and Cleaning: The raw sensor data obtained in step S1 is standardized and preprocessed, specifically including the following steps:

[0123] S21. Time Alignment: By constructing a unified time axis, data with different sampling frequencies are time aligned. The sliding window method is used to extract feature values ​​from high-frequency data, and low-frequency data is aligned through linear interpolation.

[0124] S22. Spatial Registration: Register the spatial position of the sensor to the slope engineering coordinate system according to the sensor type;

[0125] S23. Anomaly Removal: Anomaly detection and cleaning are performed on outliers, drift interference, and duplicate records in the original sensor data. Specifically, for continuous signals of displacement, strain, and tilt angle, the Z-score algorithm combined with a sliding time window is used to remove abrupt values ​​based on three times the standard deviation above and below the mean. For pore water pressure and rainfall monitoring data, physically reasonable interval thresholds are set for interval filtering. For high-frequency acceleration signals, a low-pass filtering process is preset to suppress spike noise, ultimately retaining physically reasonable and statistically stable data samples.

[0126] S24. Unit conversion and format unification: Standardize the units of physical quantities and output them in a structured manner as a time series data table with a unified field format;

[0127] S3. State Variable Extraction: Based on the cleaned and aligned structured multi-source monitoring data from step S2, key physical quantities of slope state evolution are extracted to construct a set of state variables for model identification, stability assessment, and evolution simulation. This includes the following steps:

[0128] S31. Three-dimensional displacement field reconstruction: Based on discrete displacement observations from GNSS monitoring points, multi-point displacement gauges, and image feature points, a spatial correlation function between each observation point is first constructed. Then, the displacement of grid points without sensors in the slope area is estimated using the Kriging interpolation method to form a spatially continuous initial displacement field. Subsequently, a displacement evolution state transition equation is established on a unified time axis. The spatial interpolation results of the Kriging method are used as the observation input, and the displacement estimate of each new node is recursively calculated using the Kalman filtering method. Finally, the slope mechanical equilibrium equation is introduced as a soft constraint to obtain three-dimensional displacement field data that satisfies physical consistency and spatiotemporal continuity.

[0129] S32. Strain Tensor Field Inversion: Based on distributed fiber strain data, the parameters are determined by using the Tikhonov regularization method and combining it with the L-curve criterion to invert the strain tensor field inside the slope, obtain a continuous and smooth strain distribution, and evaluate its uncertainty.

[0130] S33. Pore pressure state calculation: Based on the monitoring data of the piezometer, the pore water pressure is converted into groundwater head and effective stress to form groundwater state variables, and its spatial distribution is determined in combination with the deployment depth.

[0131] S34. Acceleration response extraction: The acceleration signal is filtered to extract the effective response components, and the peak value, energy, frequency and other indicators are calculated to characterize the structural disturbance and vibration response features.

[0132] S4. Multi-source heterogeneous data fusion: Based on the state variables extracted in step S3, state information from different types of sensors is fused to construct a unified state representation for physical modeling and stability prediction, achieving spatiotemporal joint estimation. This includes the following steps:

[0133] S41. Constructing a state vector: Constructing a state vector composed of three-dimensional displacement, principal strain components, pore pressure, and acceleration eigenvalues;

[0134] S42. State-space modeling: The observation equation Z based on the state vector. t =h(x t )+v t and state prediction equation In the observation equation, Z t Let represent the observation data from multiple sensors at time t, h(·) represent the nonlinear mapping relationship between various sensors and the state, and x t The state vector at time t represents the three-dimensional displacement, principal strain components, pore water pressure, and acceleration eigenvalues, v. t For the observation noise term; in the state prediction equation, f(·) represents the nonlinear evolution function considering displacement-strain geometry, pore water pressure diffusion, and quasi-static or dynamic evolution processes. Represents the state vector x based on the previous time step. t-1 The predicted state of the slope, w t This is the state prediction noise term;

[0135] S43, Observation-based update: Where K t Z is the Kalman gain coefficient. t The observed values ​​are used to achieve spatiotemporal joint estimation of multi-source states;

[0136] S5. Construction of three-dimensional finite element model: Based on the multi-source heterogeneous data updated in step S4, a three-dimensional slope geometric model is then constructed using digital elevation model or point cloud data. A finite element model containing geological units, physical parameters and support structure is established, and mesh refinement is implemented in the potential sliding zone.

[0137] S6. Simulation-driven and initial solution: Apply boundary loads and pore pressure conditions to the established finite element model, carry out nonlinear static or dynamic analysis, and obtain the predicted displacement and stress response distribution.

[0138] S7. Model feedback correction mechanism:

[0139] S71. State error triggering condition: Define the residual threshold ε of multidimensional state variables such as displacement, strain, and pore pressure, and compare the mean square error (MSE) of the simulation results with the measured state monitoring parameters in real time. When MSE>ε, the model correction mechanism is triggered.

[0140] S72. Parameter Inversion and Model Update: Construct a nonlinear least squares objective function with the goal of minimizing the residuals. The optimization function is F = min θ ||Y sim (θ)-Y obs || 2 +λ||θ-θ prior || 2 Where θ is the model parameter and θ={G,φ,c}, representing the numerical simulation parameters of shear modulus, friction angle and cohesion, respectively; Y sim (θ) represents the simulation response, Y obs Let λ be the observed state variable, λ be the regularization factor, and θ be the θ value. prior The parameters are prior estimates, i.e., empirical values ​​selected based on field conditions; the optimization results are used to correct material parameters and pore water distribution field, and then reapplied to the finite element model to achieve dynamic closed-loop iterative updates.

[0141] S8. Construction of Landslide Risk Index:

[0142] S81. Risk Variable Extraction: Extract the characteristic quantities, including the derivative, increment and volatility of multi-source state variables over time, and construct a dynamic indicator set for landslide trend identification and risk assessment.

[0143] S82. Risk Formula Expression: Based on the comprehensive simulation of variables including safety factor, acceleration response, pore pressure growth rate, and strain increment, a weighted landslide hazard index function is constructed to determine the landslide evolution stage and risk level range.

[0144] S9. Landslide Risk Identification and Intelligent Early Warning: When the risk indicators exceed the threshold, an early warning is automatically triggered, and the pre-landslide area and risk level are output.

[0145] S10, Model Closure and Periodic Update: Newly collected data is input into step S5 to complete the state update and risk reassessment, forming a closed-loop iteration.

[0146] In a specific embodiment, the GNSS sensor obtains the displacement in step S1 using the three-dimensional coordinate difference method as follows:

[0147]

[0148] Where Δx, Δy, and Δz are the displacement changes in the x, y, and z directions, respectively;

[0149] Displacement of discrete points measured by multi-point displacement gauge {u i The continuous displacement field is calculated using an interpolation function as follows:

[0150]

[0151] Where, φ i (z) is the interpolation shape function corresponding to the i-th measuring point, representing the displacement at position z caused by the displacement u at the i-th measuring point. i The contributed weighting factors are determined by the selected interpolation algorithm, including linear, spline, or finite element shape functions, and satisfy the interpolation condition φ. i (z j )=δ ij And the constraint that the sum of the weights is 1;

[0152] A monocular camera continuously acquires a sequence of images of the slope surface under fixed installation conditions. Pixel displacement is calculated using the optical flow or feature point matching results of images from adjacent time points. Let I... t (x,y) and I t-1 (x, y) represent the grayscale images acquired at times t and t-1, respectively. Assuming that the grayscale remains constant over time, the optical flow constraint equation is as follows:

[0153] I x ·u+I y ·v+I t =0

[0154] Among them, I x I y It These are the gradients of the image in the x, y, and time directions, respectively, where u and v0 are the planar displacement components of the pixels.

[0155] By solving for (u,v) through feature matching or dense optical flow, and combining the camera imaging model with ground calibration, pixel displacement is converted into physical displacement:

[0156]

[0157] Among them, s x s y Given the pixel size and Z as the depth distance of the corresponding point, the surface displacement vector field Δr(x,y)=(ΔX,ΔY) of the monitoring area can be obtained.

[0158] As a specific embodiment, the strain ε acquired by the distributed fiber optic strain sensor in step S1 t Calculated by the following formula:

[0159]

[0160] Where ΔL(s,t) represents the displacement increment of the monitoring segment, L0 is the initial length, and s is the position variable along the fiber optic axis.

[0161] As a specific embodiment, the pore water pressure p measured by the piezometer in step S1 t The following relationship must be satisfied:

[0162] p t =ρ w ·g·h t

[0163] Where, ρ w Let g be the density of water, g be the acceleration due to gravity, and h be the acceleration due to gravity. t This refers to the groundwater head height.

[0164] As a specific embodiment, the acceleration a obtained by the accelerometer in step S1 t The result is obtained by the following formula after filtering:

[0165] a t =LPF(a raw (t))

[0166] Where LPF represents the low-pass filter function, a raw (t) represents the original sampled signal.

[0167] As a specific embodiment, the state variable extraction method in step S3 includes the following modeling and reconstruction mechanism:

[0168] (1) Construct a structured multi-source monitoring dataset and integrate GNSS displacement data uGNSS (t), satellite imagery data u SID (t), Image recognition feature point displacement u IMG (t), fiber strain data ε DFOS The monitoring state vector is constructed by aligning the pore pressure data p(t) from the piezometer and the accelerometer data a(t) with the spatial location (x, y, z) according to a unified time series t:

[0169] x(x,y,z,t)=[u(x,y,z,t),ε(x,y,z,t),p(x,y,z,t),φ(t)]

[0170] Wherein, φ(t) represents the external excitation acting on the slope system at time t. Its value can include rainfall time series, ground motion acceleration, construction load changes, etc., to reflect the influence of the external environment on the evolution of the slope state.

[0171] (2) Based on a three-dimensional spatial grid, the continuous three-dimensional displacement field u is reconstructed using the Kriging interpolation method. 3D (x,y,z,t) is used, and the time series is dynamically corrected using the Kalman filter algorithm to improve its continuity and consistency in the spatiotemporal dimensions, satisfying the following optimization objectives:

[0172]

[0173] Among them, w i ·||u i -u(x p ,y p ,z p ,t i )|| 2 The objective is to minimize the residuals of the observed data. To optimize the obtained three-dimensional displacement field estimation vector, which contains three components: x, y, and z, u(x p ,y p ,z p ,t i ) represents the spatial location (x) of the observation point. p ,y p ,z p ) and time t i The displacement vector is obtained by interpolation or solving for the continuous displacement field to be determined; U obs ={u GNSS ,u SID ,u IMG} represents the fused multi-source displacement observation data set, w i The weighting coefficient for the i-th observation sensor reflects the data confidence level of different sensors; u i This represents the displacement vector measured by the i-th observation sensor; λ1 is the soft-wire term in the slope mechanics equilibrium equation, and λ1 is the regularization coefficient used to weight the equilibrium observation fitting term and the physical constraint term. ρ represents the divergence of the stress tensor σ corresponding to the displacement field u, and ρ represents the bulk density of the slope material.

[0174] (3) Based on fiber strain observations, the continuous spatial strain tensor field ε(x,y,z,t) is obtained by inversion using the Tikhonov regularization method. The optimized model is as follows:

[0175]

[0176] Where A is the measurement projection matrix, d is the observed strain data, λ2 is the regularization parameter, and L is the regularization operator;

[0177] (4) Convert the pore water pressure data p(t) into head variable. And calculate the effective stress variable σ'=σ-p t These constitute groundwater state variables;

[0178] (5) After the acceleration signal is processed by bandpass or lowpass filtering, the typical characteristic parameters of the disturbance response are extracted as follows:

[0179] Peak acceleration: A peak =max|a(t)|

[0180]

[0181] Energy index: E=∫|a(t)| 2 dt

[0182]

[0183] Where T is the signal analysis duration, i.e., the length of the time interval covered by integration or statistical calculation, and f is the frequency variable used to describe the distribution of the acceleration signal in the frequency domain. This represents the Fourier transform operator, which transforms a time-domain signal to the frequency domain. This represents the complex amplitude of the acceleration signal a(t) at frequency f;

[0184] Ultimately, a multidimensional spatiotemporal state variable set x(x,y,z,t) is formed, which is used to drive the modeling and evolution simulation of the digital twin of the slope.

[0185] As a specific embodiment, the data fusion algorithm in step S4 employs extended Kalman filtering, including:

[0186] Constructing the state vector:

[0187]

[0188] Status information:

[0189]

[0190] P t =F t P t-1 F t +Q t

[0191] Where, x t-1 This represents the system's state vector at the previous time t-1, which serves as the input for predicting the current state; w t For state prediction noise; P t To estimate the state at time t The magnitude of uncertainty and the correlation between data in each dimension; F t Q is the Jacobian matrix for the state transition; t The measured noise covariance of the sensor;

[0192] Observation Update:

[0193]

[0194] Among them, K t H is the Jacobian matrix of the observation model. t Let be the Jacobian matrix of the observation model, representing the sensitivity of the state vector at time t to the linearization of the observation equation. R is the transpose of the Jacobian matrix of the observation model, with dimensions n×m, used to map the observation space error back to the state space. t To observe the noise covariance;

[0195]

[0196] P t =(IK t H t )P t

[0197] Among them, K t Z is the Kalman gain coefficient. t The observed values, where I is the identity matrix.

[0198] In a specific embodiment, the digital elevation model or point cloud data in step S5 constructs a three-dimensional slope geometric model. The three-dimensional slope geometric model is reconstructed from a digital elevation model formed by fusing multi-source surveying and mapping data. The data includes: reflectance images I acquired from satellite remote sensing images. sat (x,y), image sequence acquired by a monocular slope camera The elevation difference data obtained by the structured optical flow SfM method, and the high-precision GNSS three-dimensional coordinate points p deployed at the slope control points. i =(x i ,y i ,z i The fused elevation field is modeled using the following interpolation expression:

[0199]

[0200] Where N represents the number of high-precision control points provided by the GNSS data source, M represents the number of elevation difference data points obtained by optical flow SfM method inversion from a monocular camera, and L represents the number of elevation points provided by satellite remote sensing image data. Let be the elevation value of the i-th GNSS control point. Let j be the elevation value of the SfM inversion point of the j-th monocular camera. Let ω be the elevation value of the k-th satellite remote sensing elevation point. i η j μ k These are the normalized weighting coefficients allocated based on the accuracy and spatial coverage density of each data source, satisfying the following normalization condition:

[0201]

[0202] Based on the soil and rock stratification data provided by the engineering survey, the slope area Ω is divided into k geological units Ω. k Each element is assigned the following material parameter vector:

[0203] θ k =(E k ,v k ,c k ,φ k ,ρ k ,k k )

[0204] Among them, E k For the elastic modulus, v k For Poisson's ratio, c k For cohesion, φ k ρ is the internal friction angle. k For density, k k Permeability coefficient;

[0205] If a potential slip surface exists in the geological structure Γ slip Therefore, a zero-thickness contact element is used for modeling. This contact interface follows the Mohr-Coulomb shear failure criterion, expressed as:

[0206] τ≤c slip +σ·tanφ slip

[0207] Where τ is the shear stress, σ is the normal stress, and c slip φ slip These are the cohesion of the slip surface and the internal friction angle, respectively.

[0208] If the slope is supported by anchor bolts, frame beams, anti-slide piles, etc., and finite element modeling is performed using rod elements or shell elements respectively, then the overall stiffness matrix is ​​uniformly represented as:

[0209] K global =K soil +K support +K interface

[0210] Among them, K soil K is the global stiffness matrix of the soil element. support For the overall stiffness matrix of the support structure unit, K interface For interface element stiffness matrix, represent the independent contribution terms of the stiffness of the anchor-soil coupling interface;

[0211] The model region Ω is discretized using tetrahedral elements, and the slip surface Γ slip The region undergoes mesh refinement processing, and the cell scale within the refined region satisfies:

[0212]

[0213] Among them, h global This represents the typical element size in the unrefined region. The refinement process aims to improve the analytical accuracy of the displacement and shear stress field distributions, and to meet the numerical stability and resolution requirements near the slip surface.

[0214] Regarding boundary conditions, a fully fixed constraint is set at the bottom of the slope model, namely μ. x =μ y =μ z =0, μ x μ y μ z These represent the displacement constraint coefficients in the three directions; the sliding surface is configured with contact-friction boundary conditions.

[0215] The final three-dimensional finite element model includes: the fused surface geometry Ω surf Multi-layered geological unit structure Ω k and its parameter distribution θ k Slip surface Γ slip The support structure system and its coupling relationship with the slope soil and rock mass, the non-uniform mesh division scheme and the corresponding boundary conditions provide a physical basis for subsequent numerical simulation evolution and model feedback iteration.

[0216] As a specific embodiment, the model feedback mechanism in step S7 includes the following judgment conditions and parameter update steps:

[0217] If the following error conditions are met:

[0218]

[0219] The inversion process is then initiated, and the following optimization objective is fitted using a nonlinear minimization square method:

[0220]

[0221] in, This represents the simulated displacement of the i-th monitoring point. δ represents the measured displacement of the i-th monitoring point. u δ represents the relative error between displacement and strain. The simulated and measured strains are represented by N, which represents the number of observation data points; θ = {G, φ, c} represent the numerical simulation parameters of shear modulus, friction angle, and cohesion, respectively.

[0222] As a specific embodiment, the landslide hazard index K in step S8 risk Satisfy the following expression:

[0223]

[0224] Where w1, w2, w3, and w4 are the weighting coefficients of each indicator. For displacement rate, For strain rate, The rate of change of pore water pressure. This represents the rate of change of the safety factor.

[0225] To better understand the multi-source heterogeneous data fusion method for slope digital twin modeling provided by this invention, a detailed explanation will be given below with reference to a specific implementation case. This implementation case involves field testing of the slope of Longmenqiao Reservoir in Liangjiang New Area–Changshou District. The specific steps are as follows:

[0226] S1. On-site multi-source sensor deployment and geometric modeling:

[0227] In the slope monitoring area, GNSS-V2 reference points and moving points, multi-point displacement gauges in vertical shafts, distributed fiber optic strain cables, accelerometers (RF), piezometers, high-support wireless inclinometers, crack gauges, rain gauges, monocular cameras, and satellite remote sensing systems were deployed to simultaneously collect information on the three-dimensional displacement, strain changes, pore water pressure, and image changes of the slope surface and deep layers. DEM elevation data was obtained using satellite remote sensing, and the slope geometry was parameterized by combining it with geological survey data to complete the initial three-dimensional slope modeling.

[0228] Figure 2 As shown, the GNSS displacement curves, tilt angle changes, pore water pressure and rainfall information collected by the multi-source sensors on site reflect the spatiotemporal evolution characteristics of slope deformation, strain and seepage state.

[0229] Figure 3 As shown, remote sensing imagery, on-site scene data acquisition, and GNSS sensor data acquisition were conducted, followed by the generation of DEM data for 3D modeling. A mesh was then generated using software, and a Python script for finite element preprocessing was written to refine the mesh. Finally, a nonlinear static-dynamic coupled solution was performed based on the ABAQUS platform, outputting multi-field response results including stress, displacement, and shear stress.

[0230] S2. Data Cleaning and Standardization:

[0231] A unified timestamp system is used to align asynchronous sampled data. High-frequency data is extracted for peak features through a sliding window, while low-frequency data is synchronized by linear interpolation. The spatial positions of various sensors are uniformly registered to the engineering three-dimensional coordinate system. The Z-score algorithm and interval discrimination method are used to remove abrupt changes, drifts, and outliers. All data have unified dimensions and standardized encoding, and the structured output is time-series data in a unified format. Figure 3 Steps 4-5 demonstrate the DEM generation and orthophoto matching, laying the foundation for subsequent 3D modeling.

[0232] S3. State variable extraction and physics field construction:

[0233] Based on discrete observation data from GNSS sensors, displacement gauges, and image feature points, an initial displacement field distribution is first established using Kriging spatial interpolation. Then, a Kalman filter is introduced to dynamically correct the time series, and slope equilibrium constraints are fused to reconstruct a three-dimensional displacement field with spatiotemporal continuity and physical consistency. Subsequently, an observation matrix is ​​constructed based on strain data acquired by BOTDA (a distributed fiber optic strain sensor based on Brillouin scattering). Tikhonov regularization is introduced to invert the strain tensor field, addressing local measurement gaps and ill-posed field inversion issues. For pore pressure monitoring values, groundwater state variables are calculated using head conversion and effective stress theory, and the pore pressure distribution surface is reconstructed by combining burial depth and interlayer structure. Acceleration signals are processed by bandpass filtering to extract dynamic response indices such as dominant frequency, RMS, and short-time energy, used to characterize disturbance intensity and propagation direction, ultimately forming a dynamic response model including x=[u t ,ε t ,p t ,a t The set of multiple state variables, including [ ], serves as the core input for subsequent twin model-driven and evolutionary feedback.

[0234] S4. Constructing a state space and multi-source data fusion mechanism:

[0235] Based on the state vector, a nonlinear evolution model and an observation model are constructed. The extended Kalman filter (EKF) method is used for state estimation and updating to achieve a dynamic expression of the actual slope state. This step provides real-time driving data for the digital twin.

[0236] S5. Construction and Simulation Solution of Three-Dimensional Finite Element Model:

[0237] A finite element model was constructed based on the on-site DEM and survey data, dividing the model into different soil and rock layers. Support structures (frame beams, anchors) and interface contact models were introduced. Boundary conditions included reservoir water pressure and pore water pressure loading, and the mesh was refined in the subsidence zone. The ABAQUS platform was used for nonlinear static-dynamic coupled solution to obtain multi-field response results for slope stress, displacement, and shear stress.

[0238] S6. Model feedback correction and parameter inversion mechanism:

[0239] After each simulation, the simulation results are compared with the measured state data. The mean squared residual is defined as the model error index. When it exceeds a set threshold, the parameter inversion mechanism is triggered, and the following least squares objective function is constructed for correction:

[0240] F = min θ ||Y sim (θ)-Y obs || 2 +λ||θ-θ prior || 2

[0241] The inversion results are fed back into the finite element model to achieve iterative updates of the twin.

[0242] S7. Landslide Risk Index Construction and Early Warning Mechanism:

[0243] Extracting the rate of change of key variables Construct the following weighted risk indicator function:

[0244]

[0245] The weighting coefficients for each item are selected based on on-site parameter adjustments. If K... risk A value >0.8 triggers an early warning, and the system outputs the pre-slip zone range, high-risk level, and recommended response measures.

[0246] S8. Model closure and continuous update verification:

[0247] The entire monitoring period lasted 3 months, with the twin system performing 2 data updates and 1 feedback correction daily. This embodiment verifies that the present invention possesses excellent multi-source fusion modeling capabilities, dynamic twin response capabilities, and practical engineering applicability.

[0248] like Figure 4 As shown, the virtual-real coupling closed loop of the present invention is illustrated: the left side represents the physical space, and after cleaning and spatiotemporal registration of multi-source data such as GNSS, pore permeation pressure, and acceleration, the state variable x = [u] is extracted. t ,ε t ,p t ,a t ] T The observation field was then fused using EKF (Extended Kalman Filter) to obtain a continuous and consistent observation field. The right side is the virtual simulation space, where 3D geometric modeling and mesh refinement of the potential slip zone are completed based on DEM and geological survey data. Material parameters and boundary / loads are assigned, and then the model parameters θ are corrected according to the least squares objective function constructed in step S6, where θ is the model parameter and θ = {G,φ,c]. When MSE > ε, the parameter θ is updated and written back to the simulation model. The bottom is the risk assessment and early warning, based on {du / dt,dε / dt,dp / dt,dF}. s / dt]Calculate K risk Exceeding the threshold K th This triggers a tiered warning system.

[0249] The above process is executed at a fixed cycle, forming a closed loop of "observation-driven - simulation - calibration - write-back - evaluation", which enables the digital twin to adapt to new data and maintain physical consistency.

[0250] Compared with existing technologies, the slope digital twin modeling method based on multi-source heterogeneous data fusion provided by this invention integrates various heterogeneous monitoring data, including those from GNSS sensors, multi-point displacement gauges, distributed fiber optic strain sensors, accelerometers, piezometers, monocular cameras, and satellite remote sensing imagery. Through time alignment, spatial registration, and standardization, a unified state vector is constructed, and key physical field quantities such as three-dimensional displacement field, strain tensor field, and pore pressure response are extracted to establish a digital twin modeling system coupled with multi-field evolution. On this basis, a dynamic correction mechanism for model parameters based on "monitoring-simulation" residual triggering is proposed to achieve adaptive updates of shear modulus, friction angle, and cohesion parameters, thereby improving the model's predictive ability. Furthermore, a weighted landslide hazard index function is constructed, integrating the derivatives and fluctuation characteristics of multi-source state variables to divide risk level intervals, achieving advanced identification and intelligent early warning of landslide risks. This method has the advantages of multiple information fusion types, high model building accuracy, fast twin simulation response, and strong risk identification capability, which significantly improves the real-time performance, accuracy, and engineering adaptability of slope disaster monitoring and early warning, and can be well applied to slope risk management scenarios.

[0251] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for modeling digital twins of slopes using multi-source heterogeneous data fusion, characterized in that, Includes the following steps: S1. Multi-source sensor deployment and data acquisition: GNSS sensors, multi-point displacement gauges, distributed fiber optic strain sensors, accelerometers, piezometers, monocular cameras and satellite remote sensing imagery equipment are deployed in the target slope area to acquire three-dimensional displacement, deep displacement, strain change, vibration response, pore pressure state, image information and slope three-dimensional topography data, respectively. S2. Data Synchronization and Cleaning: The raw sensor data obtained in step S1 is standardized and preprocessed, specifically including the following steps: S21. Time Alignment: By constructing a unified time axis, data with different sampling frequencies are time aligned. The sliding window method is used to extract feature values ​​from high-frequency data, and low-frequency data is aligned through linear interpolation. S22. Spatial Registration: Register the spatial position of the sensor to the slope engineering coordinate system according to the sensor type; S23. Anomaly Removal: Anomaly detection and cleaning are performed on outliers, drift interference, and duplicate records in the original sensor data. Specifically, for continuous signals of displacement, strain, and tilt angle, the Z-score algorithm combined with a sliding time window is used to remove abrupt values ​​based on three times the standard deviation above and below the mean. For pore water pressure and rainfall monitoring data, physically reasonable interval thresholds are set for interval filtering. For high-frequency acceleration signals, a low-pass filtering process is preset to suppress spike noise, ultimately retaining physically reasonable and statistically stable data samples. S24. Unit conversion and format unification: Standardize the units of physical quantities and output them in a structured manner as a time series data table with a unified field format; S3. State Variable Extraction: Based on the cleaned and aligned structured multi-source monitoring data from step S2, key physical quantities of slope state evolution are extracted to construct a set of state variables for model identification, stability assessment, and evolution simulation. This includes the following steps: S31. Three-dimensional displacement field reconstruction: Based on the discrete displacement observations of GNSS monitoring points, multi-point displacement gauges and image feature points, the spatial correlation function between each observation point is first constructed, and the displacement of grid points in the slope area without sensor deployment is estimated using the Kriging interpolation method to form a spatially continuous initial displacement field. Subsequently, a displacement evolution state transition equation is established on a unified time axis. The spatial interpolation results of the Kriging method are used as the observation input. The displacement estimates of each new node are calculated recursively using the Kalman filter method. Then, the slope mechanical equilibrium equation is introduced as a soft constraint. Finally, three-dimensional displacement field data that meets physical consistency and spatiotemporal continuity are obtained. S32. Strain Tensor Field Inversion: Based on distributed fiber strain data, the parameters are determined by using the Tikhonov regularization method and combining it with the L-curve criterion to invert the strain tensor field inside the slope, obtain a continuous and smooth strain distribution, and evaluate its uncertainty. S33. Pore pressure state calculation: Based on the monitoring data of the piezometer, the pore water pressure is converted into groundwater head and effective stress to form groundwater state variables, and its spatial distribution is determined in combination with the deployment depth. S34. Acceleration Response Extraction: The acceleration signal is filtered to extract the effective response components, and the peak value, energy, and frequency index are calculated to characterize the structural disturbance and vibration response features. S4. Multi-source heterogeneous data fusion: Based on the state variables extracted in step S3, state information from different types of sensors is fused to construct a unified state representation for physical modeling and stability prediction, achieving spatiotemporal joint estimation. This includes the following steps: S41. Constructing the state vector: Constructing a state vector composed of three-dimensional displacement, principal strain components, pore pressure, and acceleration eigenvalues; S42. State-space modeling: The observation equation Z based on the state vector. t =h(x t )+v t and state prediction equation In the observation equation, Z t Let represent the observation data from multiple sensors at time t, h(·) represent the nonlinear mapping relationship between various sensors and the state, and x t The state vector at time t represents the three-dimensional displacement, principal strain components, pore water pressure, and acceleration eigenvalues, v. t For the observation noise term; in the state prediction equation, f(·) represents the nonlinear evolution function considering displacement-strain geometry, pore water pressure diffusion, and quasi-static or dynamic evolution processes. Represents the state vector x based on the previous time step. t-1 The predicted state of the slope, w t This is the state prediction noise term; S43, Observation-based update: Where K t Z is the Kalman gain coefficient. t The observed values ​​are used to achieve spatiotemporal joint estimation of multi-source states; S5. Construction of the 3D finite element model: based on the update in step S4 Multi-source heterogeneous data are then used to construct a three-dimensional slope geometric model using digital elevation model or point cloud data, and a finite element model containing geological units, physical parameters and support structures is established, with mesh refinement implemented in the potential sliding zone. S6. Simulation-driven and initial solution: Apply boundary loads and pore pressure conditions to the established finite element model, carry out nonlinear static or dynamic analysis, and obtain the predicted displacement and stress response distribution. S7. Model feedback correction mechanism: S71. State error triggering condition: Define the residual threshold ε of multidimensional state variables such as displacement, strain, and pore pressure, and compare the mean square error (MSE) of the simulation results with the measured state monitoring parameters in real time. When MSE > ε, the model correction mechanism is triggered. S72. Parameter Inversion and Model Update: Construct a nonlinear least squares objective function with the goal of minimizing the residuals. The optimization function is F = min θ ||Y sim (θ)-Y obs || 2 +λ||θ-θ prior || 2 Where θ is the model parameter and θ={G,φ,c}, representing the numerical simulation parameters of shear modulus, friction angle and cohesion, respectively; Y sim (θ) represents the simulation response, Y obs Let λ be the observed state variable, λ be the regularization factor, and θ be the θ value. prior The parameters are prior estimates, i.e., empirical values ​​selected based on field conditions; the optimization results are used to correct material parameters and pore water distribution field, and then reapplied to the finite element model to achieve dynamic closed-loop iterative updates. S8. Construction of Landslide Risk Index: S81. Risk Variable Extraction: Extract the characteristic quantities, including the derivative, increment and volatility of multi-source state variables over time, and construct a dynamic indicator set for landslide trend identification and risk assessment. S82. Risk Formula Expression: Based on the comprehensive simulation of variables including safety factor, acceleration response, pore pressure growth rate, and strain increment, a weighted landslide hazard index function is constructed to determine the landslide evolution stage and risk level range. S9. Landslide Risk Identification and Intelligent Early Warning: When the risk indicators exceed the threshold, an early warning is automatically triggered, and the pre-landslide area and risk level are output. S10, Model Closure and Periodic Update: Newly collected data is input into step S5 to complete the state update and risk reassessment, forming a closed-loop iteration.

2. The method for modeling a digital twin of a slope using multi-source heterogeneous data fusion according to claim 1, characterized in that, In step S1, the GNSS sensor obtains the displacement using the three-dimensional coordinate difference method as follows: Where Δx, Δy, and Δz are the displacement changes in the x, y, and z directions, respectively; Displacement of discrete points measured by multi-point displacement gauge {u i The continuous displacement field is calculated using an interpolation function as follows: Where, φ i (z0 is the interpolation shape function corresponding to the i-th measuring point, indicating that the displacement at position z is caused by the displacement u of the i-th measuring point.) i The contributed weighting factors are determined by the selected interpolation algorithm, including linear, spline, or finite element shape functions, and satisfy the interpolation condition φ. i (z j )=δ ij And the constraint that the sum of the weights is 1; A monocular camera continuously acquires a sequence of images of the slope surface under fixed installation conditions. Pixel displacement is calculated using the optical flow or feature point matching results of images from adjacent time points. Let I... t (x,y) and I t-1 (x, y) represent the grayscale images acquired at times t and t-1, respectively. Assuming that the grayscale remains constant over time, the optical flow constraint equation is as follows: I x ·u+I y ·v+I t =0 Among them, I x I y I t These are the gradients of the image in the x, y, and time directions, respectively, where u and v0 are the planar displacement components of the pixels. By solving for (u,v) through feature matching or dense optical flow, and combining the camera imaging model with ground calibration, pixel displacement is converted into physical displacement: Among them, s x s y Given the pixel size and Z as the depth distance of the corresponding point, the surface displacement vector field Δr(x,y)=(ΔX,ΔY) of the monitoring area can be obtained.

3. The method for modeling a digital twin of a slope using multi-source heterogeneous data fusion according to claim 1, characterized in that, In step S1, the strain ε acquired by the distributed fiber optic strain sensor t Calculated by the following formula: Where ΔL(s,t) represents the displacement increment of the monitoring segment, L0 is the initial length, and s is the position variable along the fiber optic axis.

4. The method for modeling a digital twin of a slope using multi-source heterogeneous data fusion according to claim 1, characterized in that, The pore water pressure p measured by the piezometer in step S1 t The following relationship must be satisfied: p t =ρ w ·g·h t Where, ρ w Let g be the density of water, g be the acceleration due to gravity, and h be the acceleration due to gravity. t This refers to the groundwater head height.

5. The method for modeling a digital twin of a slope using multi-source heterogeneous data fusion according to claim 1, characterized in that, In step S1, the acceleration a obtained by the accelerometer t The result is obtained by the following formula after filtering: a t =LPF(a raw (t)) Where LPF represents the low-pass filter function, a raw (t) represents the original sampled signal.

6. The method for modeling a digital twin of a slope using multi-source heterogeneous data fusion according to claim 1, characterized in that, The state variable extraction method in step S3 includes the following modeling and reconstruction mechanism: (1) Construct a structured multi-source monitoring dataset and integrate GNSS displacement data u GNSS (t), satellite imagery data u SID (t), Image recognition feature point displacement u IMG (t), fiber strain data ε DFOS The monitoring state vector is constructed by aligning the pore pressure data p(t) from the piezometer and the accelerometer data a(t) with the spatial location (x, y, z) according to a unified time series t: x(x,y,z,t)=[u(x,y,z,t),ε(x,y,z,t),p(x,y,z,t),φ(t)] Wherein, φ(t) represents the external excitation acting on the slope system at time t. Its value can include rainfall time series, ground motion acceleration, construction load changes, etc., to reflect the influence of the external environment on the evolution of the slope state. (2) Based on a three-dimensional spatial grid, the continuous three-dimensional displacement field u is reconstructed using the Kriging interpolation method. 3D (x,y,z,t) is used, and the time series is dynamically corrected using the Kalman filter algorithm to improve its continuity and consistency in the spatiotemporal dimensions, satisfying the following optimization objectives: Among them, w i ·||u i -u(x p ,y p ,z p ,t i )|| 2 The objective is to minimize the residuals of the observed data. To optimize the obtained three-dimensional displacement field estimation vector, which contains three components: x, y, and z, u(x p ,y p ,z p ,t i ) represents the spatial location (x) of the observation point. p ,y p ,z p ) and time t i The displacement vector is obtained by interpolation or solving for the continuous displacement field to be determined; U obs ={u GNSS ,u SID ,u IMG } represents the fused multi-source displacement observation data set, w i The weighting coefficient for the i-th observation sensor reflects the data confidence level of different sensors; u i This represents the displacement vector measured by the i-th observation sensor; λ1 is the soft-wire term in the slope mechanics equilibrium equation, and λ1 is the regularization coefficient used to weight the equilibrium observation fitting term and the physical constraint term. ρ represents the divergence of the stress tensor σ corresponding to the displacement field u, and ρ represents the bulk density of the slope material. (3) Based on fiber strain observations, the continuous spatial strain tensor field ε(x,y,z,t) is obtained by inversion using the Tikhonov regularization method. The optimized model is as follows: Where A is the measurement projection matrix, d is the observed strain data, λ2 is the regularization parameter, and L is the regularization operator; (4) Convert the pore water pressure data p(t) into head variable. And calculate the effective stress variable σ'=σ-p t These constitute groundwater state variables; (5) After the acceleration signal is processed by bandpass or lowpass filtering, the typical characteristic parameters of the disturbance response are extracted as follows: Peak acceleration: A peak =max|a(t)| Effective value (RMS): Energy index: E=∫|a(t)| 2 dt Clock speed specifications: Where T is the signal analysis duration, i.e., the length of the time interval covered by integration or statistical calculation, and f is the frequency variable used to describe the distribution of the acceleration signal in the frequency domain. This represents the Fourier transform operator, which transforms a time-domain signal to the frequency domain. This represents the complex amplitude of the acceleration signal a(t) at frequency f; Ultimately, a multidimensional spatiotemporal state variable set x(x,y,z,t) is formed, which is used to drive the modeling and evolution simulation of the digital twin of the slope.

7. The method for modeling a digital twin of a slope using multi-source heterogeneous data fusion according to claim 6, characterized in that, The data fusion algorithm in step S4 employs an extended Kalman filter, including: Constructing the state vector: Status information: P t =F t P t-1 F t +Q t Where, x t-1 This represents the system's state vector at the previous time t-1, which serves as the input for predicting the current state; w t For state prediction noise; P t To estimate the state at time t The magnitude of uncertainty and the correlation between data in each dimension; F t Q is the Jacobian matrix for the state transition; t The measured noise covariance of the sensor; Observation Update: Among them, K t H is the Jacobian matrix of the observation model. t Let be the Jacobian matrix of the observation model, representing the sensitivity of the state vector at time t to the linearization of the observation equation. R is the transpose of the Jacobian matrix of the observation model, with dimensions n×m, used to map the observation space error back to the state space. t To observe the noise covariance; P t =(I-K t H t )P t Among them, K t Z is the Kalman gain coefficient. t The observed values, where I is the identity matrix.

8. The method for modeling a digital twin of a slope using multi-source heterogeneous data fusion according to claim 1, characterized in that, In step S5, a three-dimensional slope geometric model is constructed from digital elevation model or point cloud data. This three-dimensional slope geometric model is reconstructed from a digital elevation model formed by fusing multi-source surveying data. The data includes: reflectance image I acquired from satellite remote sensing imagery. sat (x,y), image sequence acquired by a monocular slope camera The elevation difference data obtained by the structured optical flow SfM method, and the high-precision GNSS three-dimensional coordinate points p deployed at the slope control points. i =(x i ,y i ,z i The fused elevation field is modeled using the following interpolation expression: Where N represents the number of high-precision control points provided by the GNSS data source, M represents the number of elevation difference data points obtained by optical flow SfM method inversion from a monocular camera, and L represents the number of elevation points provided by satellite remote sensing image data. Let be the elevation value of the i-th GNSS control point. Let j be the elevation value of the SfM inversion point of the j-th monocular camera. Let ω be the elevation value of the k-th satellite remote sensing elevation point. i η j μ k These are the normalized weighting coefficients allocated based on the accuracy and spatial coverage density of each data source, satisfying the following normalization condition: Based on the soil and rock stratification data provided by the engineering survey, the slope area Ω is divided into k geological units Ω. k Each element is assigned the following material parameter vector: i k =(E k ,v k ,c k ,f k ,r k ,k k ) Among them, E k For the elastic modulus, v k For Poisson's ratio, c k For cohesion, φ k ρ is the internal friction angle. k For density, k k Permeability coefficient; If a potential slip surface exists in the geological structure Γ slip Therefore, a zero-thickness contact element is used for modeling. This contact interface follows the Mohr-Coulomb shear failure criterion, expressed as: τ≤c slip +σ·tanφ slip Where τ is the shear stress, σ is the normal stress, and c slip φ slip These are the cohesion of the slip surface and the internal friction angle, respectively. If the slope is supported by anchor bolts, frame beams, and anti-slide piles, and finite element modeling is performed using either rod elements or shell elements respectively, then the overall stiffness matrix can be uniformly represented as: K global =K soil +K support +K interface Among them, K soil K is the global stiffness matrix of the soil element. support For the overall stiffness matrix of the support structure unit, K interface For interface element stiffness matrix, represent the independent contribution terms of the stiffness of the anchor-soil coupling interface; The model region Ω is discretized using tetrahedral elements, and the slip surface Γ slip The region undergoes mesh refinement processing, and the cell scale within the refined region satisfies: Among them, h global This represents the typical element size in the unrefined region. The refinement process aims to improve the analytical accuracy of the displacement and shear stress field distributions, and to meet the numerical stability and resolution requirements near the slip surface. Regarding boundary conditions, a fully fixed constraint is set at the bottom of the slope model, namely μ. x =μ y =μ z =0, μ x μ y μ z These represent the displacement constraint coefficients in the three directions; the sliding surface is configured with contact-friction boundary conditions. The final three-dimensional finite element model includes: the fused surface geometry Ω surf Multi-layered geological unit structure Ω k and its parameter distribution θ k Slip surface Γ slip The support structure system and its coupling relationship with the slope soil and rock mass, the non-uniform mesh division scheme and the corresponding boundary conditions provide a physical basis for subsequent numerical simulation evolution and model feedback iteration.

9. The method for modeling a digital twin of a slope using multi-source heterogeneous data fusion according to claim 1, characterized in that, The model feedback mechanism in step S7 includes the following judgment conditions and parameter update steps: If the following error conditions are met: The inversion process is then initiated, and the following optimization objective is fitted using a nonlinear minimization square method: in, This represents the simulated displacement of the i-th monitoring point. δ represents the measured displacement of the i-th monitoring point. u δ ε This represents the relative error between displacement and strain. The simulated and measured strains are represented by N, which represents the number of observation data points; θ = {G, φ, c} represent the numerical simulation parameters of shear modulus, friction angle, and cohesion, respectively.

10. The method for modeling a digital twin of a slope using multi-source heterogeneous data fusion according to claim 1, characterized in that, In step S8, the landslide hazard index K risk Satisfy the following expression: Where w1, w2, w3, and w4 are the weighting coefficients of each indicator. For displacement rate, For strain rate, The rate of change of pore water pressure. This represents the rate of change of the safety factor.

Images