Geological disaster pre-disaster risk early warning method based on multi-source heterogeneous data fusion

Abstract

This invention discloses a method for early warning of geological disaster risks based on multi-source heterogeneous data fusion, belonging to the field of geological disaster monitoring and early warning technology. A three-dimensional geomechanical numerical model of the target geological disaster body is established, a multi-source heterogeneous data stream is constructed, and a high-confidence observation sequence is generated. Using a physical information neural network as a surrogate model, the residuals of the geomechanical control equations are embedded into a loss function to invert unsteady mechanical parameters in real time, constructing a digital twin that matches the actual state of the geological body. Based on the digital twin and short-term weather forecasts, multi-time-step forward dynamics extrapolation is carried out to calculate the probability of the cumulative plastic zone and potential sliding surface connection. By reconstructing the phase space, the Lyapunov exponent spectrum is calculated, and the chaotic instability critical state is determined based on the maximum exponent turning from negative to positive, generating an early warning signal. This invention achieves deep integration of data-driven and physical mechanism constraints, improving the accuracy, timeliness, and reliability of early warning, and is suitable for high-precision early warning of geological disasters such as slope erosion.

Description

Technical Field

[0001] This invention relates to the field of geological disaster monitoring and early warning technology, specifically to a method for early warning of impending geological disaster risks based on the fusion of multi-source heterogeneous data. Background Technology

[0002] my country's geological environment is complex and fragile, with frequent geological disasters such as landslides, collapses, and debris flows. These disasters are characterized by their numerous locations, wide distribution, suddenness, and high risk of causing damage, posing a continuous threat to engineering construction, ecological security, and people's lives and property. Traditional geological disaster early warning systems often rely on empirical indicators such as rainfall thresholds and displacement rates, or use simplified limit equilibrium models for stability calculations. These methods struggle to depict the entire nonlinear evolution of soil and rock masses under the influence of rainfall infiltration, stress redistribution, and seepage-stress coupling, often resulting in delayed warnings, frequent false alarms and missed alarms, and ambiguous critical criteria. This makes it difficult to meet the high-precision and timeliness requirements for early warning of disasters on steep slopes in mountainous areas and complex geological formations along reservoir banks. With the rapid development of integrated space-air-ground monitoring technology, multi-source methods such as space-based InSAR, UAV LiDAR, ground-based GNSS, and deep underground stress sensing are becoming increasingly widespread, enriching the data dimensions. However, multi-source heterogeneous data suffers from bottlenecks such as inconsistent spatiotemporal benchmarks, strong noise interference, and insufficient utilization of information redundancy and complementarity. Simply stacking data cannot be transformed into reliable understanding of disaster evolution. There is an urgent need to establish an efficient fusion and dynamic assimilation mechanism under a unified spatiotemporal framework.

[0003] Multi-source heterogeneous data fusion has become a core technology in the field of geological disaster monitoring and early warning. Existing research mainly focuses on data layer or feature layer fusion, using methods such as Kalman filtering, Bayesian estimation, and deep learning to improve observation consistency and robustness, but generally lacks rigid constraints from geomechanical mechanisms. Xue Hongwen et al. (2026) constructed a multi-source data fusion early warning system, using convolutional neural networks to extract features and combining them with deep kernel support vector machines to achieve classification prediction, solving the problem of multi-source data consistency, but did not embed physical equations such as the constitutive relationship of soil and rock and fluid-structure interaction into the model training process. Du Yan et al. (2026) conducted an early warning research situation analysis based on a large language model, pointing out the common defects of data-driven models, such as weak generalization ability and poor interpretability under small samples and extreme conditions. Huang Runqiu and Xu Qiang (2000) systematically proposed to use phase space reconstruction and Lyapunov exponent to quantitatively determine the chaotic characteristics of geological disaster systems, providing a nonlinear dynamic basis for identifying instability critical states, but this method is mostly used for post-event analysis and has failed to couple with real-time numerical models to achieve closed-loop early warning. Existing academic research either focuses too much on data-driven approaches and lacks a clear mechanism, or it adheres to physical simulations but lacks sufficient data assimilation capabilities, making it difficult to support the accurate extrapolation of time-varying parameters and unsteady geological bodies.

[0004] Physical Information Neural Networks (PINN) and digital twins offer new pathways to solving the challenges of mechanism and data fusion. Related patented technologies attempt to integrate physical equation constraints and dynamic parameter inversion into the early warning process. Chinese patent CN121505787A discloses a geological early warning method for rockfall slopes based on digital twins. Through multi-dimensional physical perception and dynamic updates of the mechanical twin model, a digital twin is constructed to accurately replicate the slope state, which to some extent solves the problems of disconnect between traditional static modeling and actual working conditions, as well as data fragmentation. This patent uses digital twins and real-time monitoring data to drive updates, improving the matching degree between the model and the real geological body. However, there are still significant shortcomings: the residuals of partial differential equations are not embedded into the loss function to achieve physical soft constraints; the mechanical parameter inversion does not consider the time-varying characteristics of cohesion, internal friction angle, and permeability coefficient; multi-time-step coupled forward extrapolation is not carried out; and chaotic criteria such as Lyapunov exponent spectrum are not introduced to achieve quantitative identification of the critical point of disaster. The early warning threshold still relies on empirical settings, making it difficult to capture the sudden change behavior from stability to chaotic instability.

[0005] Therefore, existing technologies and methods have significant shortcomings in areas such as spatiotemporal unified fusion of multi-source heterogeneous data, deep coupling of physical mechanisms and data-driven approaches, real-time inversion of time-varying mechanical parameters, and quantitative determination of chaotic instability critical states. Traditional numerical models suffer from fixed parameters and computational time consumption; purely data-driven models lack physical constraints and have poor generalization ability; existing digital twin and intelligent algorithm fusion schemes have not achieved the organic unity of PDE constraints, adaptive parameter updates, and nonlinear dynamic criteria. Summary of the Invention

[0006] To address the aforementioned technical issues, this application discloses a method for early warning of impending geological disaster risks based on multi-source heterogeneous data fusion, specifically including:

[0007] A three-dimensional geomechanical numerical model of the target geological hazard body is established. The numerical model defines the physical framework of hazard evolution based on the constitutive relationship of rock and soil.

[0008] Construct a multi-source heterogeneous data stream, which includes space-based InSAR deformation data, airborne UAV lidar point cloud data, ground-based GNSS displacement data, and deep underground stress monitoring data.

[0009] Using the Physical Information Neural Network (PINN) as a surrogate model, the residuals of the partial differential equations of the numerical model are added as constraints to the loss function. By minimizing the difference between the monitoring data and the model predictions, the unsteady mechanical parameters in the numerical model are inverted and updated in real time. The parameters include at least time-varying cohesion, internal friction angle and permeability coefficient, thus obtaining a digital twin that matches the real state of the current geological body in real time.

[0010] Based on the updated digital twin, combined with short-term weather forecast data for future periods, forward dynamic extrapolation at multiple time steps is performed to calculate the probability of the cumulative plastic zone and potential sliding surface connecting at future moments.

[0011] Based on the simulation results, the Lyapunov index spectrum of the state is calculated. When the maximum Lyapunov index turns from negative to positive, it is determined that the geological body has entered a chaotic and unstable state, and a disaster risk warning signal is generated accordingly.

[0012] Preferably, the construction of the multi-source heterogeneous data stream further includes:

[0013] Spatiotemporal references were unified for space-based InSAR deformation data, airborne UAV lidar point cloud data, ground-based GNSS displacement data, and deep underground stress monitoring data. A Kalman filter fusion algorithm was then used to generate high-confidence observation sequences. Its state update equation is: ,in, For the true state vector, For the observation matrix, To observe the noise, Let be the noise covariance matrix.

[0014] Preferably, the step of using the Physical Information Neural Network (PINN) as a surrogate model and incorporating the partial differential equation residuals of the numerical model as constraint terms into the loss function includes:

[0015] Construct the total loss function of the physical information neural network The total loss function Data mismatch Physical equation residuals and boundary condition residuals The weighted composition is calculated using the following formula: ,in, For neural networks at position and time Output displacement prediction value, For observations in a multi-source heterogeneous data stream, For the residual operator of the governing equation, For boundary condition operators, , , These are the weighting coefficients for each item. , , These represent the number of observation data points, the number of physical equation configuration points, and the number of boundary condition points, respectively.

[0016] Preferably, the residual term of the physical equation The construction includes:

[0017] Based on Biot's consolidation theory, the residual between the mass conservation equation and the momentum conservation equation is defined as: ,in, Let ∇ be the effective stress tensor and ∇ be the divergence operator. For Biot coefficient, For the pore water pressure gradient, The vector of gravitational acceleration. displacement vector The second derivative with respect to time, Porosity For fluid density, For fluid velocity, For source and sink items.

[0018] Preferably, the real-time inversion and updating of the unsteady mechanical parameters in the numerical model includes:

[0019] The time-varying cohesion internal friction angle and permeability coefficient The parameter is parameterized as an implicit function that evolves over time. The sensitivity of the stress tensor to the strain tensor is calculated using the automatic differentiation technique of the physical information neural network. The parameter gradient is solved using the adjoint state method, and the update formula is: ,in, The set of unsteady mechanical parameters to be inverted. The cohesion of the soil and rock mass, The internal friction angle of the rock and soil mass. The permeability coefficient of the rock and soil mass. In the first The unsteady mechanical parameter vector updated in the next iteration. In the first The unsteady mechanical parameter vector of the next iteration Total loss function Regarding parameters In the The gradient vector at the next iteration For learning rate, The momentum coefficient, The regularization coefficient is . To constrain the parameters within the physically feasible region Projection operator within.

[0020] Preferably, the forward dynamics extrapolation with multiple time steps includes:

[0021] The updated dynamic equilibrium equations of the digital twin are solved using explicit-implicit hybrid integration. These equations consider the coupling effect of pore water pressure and skeleton stress, and their discretized form is as follows: ,in, , , These are the mass matrix, damping matrix, and stiffness matrix, respectively. Let be the nodal displacement vector. For external load vectors, The pore water pressure field is derived from short-term weather forecast data. It is a matrix of shape functions.

[0022] Preferably, the explicit-implicit hybrid integral includes:

[0023] The dynamic term is explicitly integrated using the central difference method, and the seepage term is implicitly integrated using the backward difference method. The coupled solution scheme is as follows:

[0024]

[0025]

[0026] in, In order to be in The nodal displacement vector at time t. In order to be in The nodal displacement vector at time t. For time step, In order to be in The nodal velocity vector at time t, It is the inverse of the mass matrix. In order to be in The external nodal load vector at time t. In order to be in The internal node resistance vector at time t; In order to be in The nodal pore water pressure vector at time t. In order to be in The nodal pore water pressure vector at time t. It is the inverse of the penetration matrix. In order to be in The external fluid flow vector at time t. This is the fluid-structure interaction matrix.

[0027] Preferably, the calculation of the probability of the cumulative plastic zone and the potential sliding surface connecting at future times includes:

[0028] Define local damage variables The local damage variable is based on cumulative plastic strain. and stress triaxiality The construction, its evolution equation is: ,in, and Let be the material damage constant. For at any time The equivalent plastic strain rate;

[0029] Based on local damage variables Calculate the probability of full connectivity. The global penetration probability is determined by calculating the seepage threshold of the damage field gradient: ,in, As a potential connecting route, This is the critical damage threshold.

[0030] Preferably, the calculation of the Lyapunov exponent spectrum of the state based on the deduction results includes:

[0031] Constructing the reconstructed phase space vector of the state space ,in It is in a displacement state. In terms of speed state, The state of time-varying mechanical parameters;

[0032] Apply small perturbations to the reconstructed phase space trajectory Solve the variational equation ,in In the state Jacobian matrix at the location;

[0033] The Gram-Schmidt orthogonalization method is used to periodically reorthogonalize the disturbance vectors of the discretized time series, and the maximum Lyapunov exponent is calculated. The formula is: Calculated by sliding time window The time series curve is used to determine the bifurcation of dynamic behavior from steady state to chaotic state when the curve crosses the preset zero value criterion threshold and the second derivative is greater than zero, thus generating a disaster risk warning signal.

[0034] Preferably, the determination that the geological body has entered a chaotic and unstable state, and the generation of a pre-disaster risk warning signal accordingly, includes:

[0035] Set tiered early warning thresholds ,when And duration At that time, a Level 1 emergency disaster warning was triggered. This is a preset critical value;

[0036] Simultaneously calculate the fractal dimension. ,like and If the characteristics of non-integer dimensions are exhibited, it is determined that the system has entered a chaotic evolution stage dominated by strange attractors.

[0037] Compared with the prior art, the technical solution of this application has the following technical effects:

[0038] This invention constructs a three-dimensional geomechanical numerical model and integrates multi-source heterogeneous monitoring data, which can fully characterize the real mechanical state and spatial distribution characteristics of geological hazard bodies. It achieves efficient integration and unified expression of space-based, air-based, ground-based, and deep monitoring information, providing a stable and reliable data and model foundation for accurate identification of geological body conditions. It ensures that subsequent inversion and early warning judgments are carried out in an orderly manner under a unified physical framework, and improves the consistency and integrity of the entire early warning process.

[0039] This invention uses a physical information neural network as a proxy model and incorporates the residuals of geomechanical partial differential equations into the loss function constraint. It can strictly follow the physical evolution law of soil and rock while being data-driven, effectively realize the real-time inversion and dynamic updating of unsteady mechanical parameters, and quickly generate a digital twin that is synchronously matched with the real state of the geological body. This significantly improves the model's adaptability and computational efficiency to complex working conditions and ensures that the model state is highly synchronized with the actual geological body.

[0040] This invention uses an updated digital twin combined with short-term meteorological data to conduct forward dynamic extrapolation over multiple time steps. It can accurately calculate the distribution of the cumulative plastic zone and the probability of potential sliding surface connection in the future time period of a geological body, clearly restore the entire process of disaster evolution from stability to local damage and then to connection and instability, provide intuitive and reliable quantitative basis for risk assessment, and fully support the early identification and trend prediction of the impending disaster state.

[0041] This invention achieves quantitative determination of chaotic instability by calculating the Lyapunov exponent spectrum, which can accurately capture the critical abrupt change characteristics of geological bodies from stability to instability. It automatically generates early warning signals of impending disaster risks based on exponent changes, improves the scientific nature and accuracy of early warning judgment, effectively realizes timely early warning of geological disasters in the impending stage, provides stable and efficient technical support for geological disaster emergency prevention and control, and ensures reliable early warning results and timely response.

[0042] The above description is only an overview of the technical solution of this application. In order to better understand the technical means of this application and implement it in accordance with the contents of the specification, and to make the above and other objects, features and advantages of this application more obvious and understandable, the preferred embodiments of this application are described in detail below with reference to the accompanying drawings.

[0043] The above and other objects, advantages and features of this application will become more apparent to those skilled in the art from the following detailed description of specific embodiments in conjunction with the accompanying drawings. Attached Figure Description

[0044] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort. In all drawings, similar elements or parts are generally identified by similar reference numerals. In the drawings, the elements or parts are not necessarily drawn to scale.

[0045] Based on the description of the figures and their corresponding technical content in the document, the titles of the figures are as follows:

[0046] Figure 1 Flowchart of a geological disaster risk early warning method based on multi-source heterogeneous data fusion;

[0047] Figure 2 A schematic diagram of the structure and loss function of the PINN (Physical Information Neural Network) model.

[0048] Figure 3 A schematic diagram of the process for inverting unsteady mechanical parameters and updating digital twins;

[0049] Figure 4 Comparison curves of fitting degree between different method models and measured mechanical states;

[0050] Figure 5 A comparison curve of the time consumption of multi-time-step dynamics extrapolation;

[0051] Figure 6 A bar chart comparing the accuracy and lead time of various early warning methods;

[0052] Figure 7 Spatiotemporal evolution cloud map of damage variables across the entire geological body;

[0053] Figure 8 The temporal variation of the maximum Lyapunov exponent and the determination of chaotic instability are shown in the figure. Detailed Implementation

[0054] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. In the following description, specific details such as specific configurations and components are provided merely to help fully understand the embodiments of this application. Therefore, those skilled in the art should understand that various changes and modifications can be made to the embodiments described herein without departing from the scope and spirit of this application. In addition, for clarity and brevity, descriptions of known functions and structures are omitted in the embodiments.

[0055] It should be understood that the phrase "an embodiment" or "this embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of this application. Therefore, "an embodiment" or "this embodiment" appearing throughout the specification does not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments.

[0056] Furthermore, reference numerals and / or letters may be repeated in different examples within this application. Such repetition is for the purpose of simplification and clarity and does not in itself indicate a relationship between the various embodiments and / or settings discussed.

[0057] In this article, the term "and / or" is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can mean: A exists alone, B exists alone, and A and B exist simultaneously. The term " / and" in this article describes another type of relationship between related objects, indicating that two relationships can exist. For example, A / and B can mean: A exists alone, and A and B exist alone. In addition, the character " / " in this article generally indicates that the related objects before and after it are in an "or" relationship.

[0058] In this article, the term "at least one" is merely a description of the relationship between related objects, indicating that there can be three relationships. For example, "at least one of A and B" can mean: A exists alone, A and B exist simultaneously, or B exists alone.

[0059] It should also be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion.

[0060] Example 1

[0061] This embodiment mainly describes a geological disaster risk early warning method based on multi-source heterogeneous data fusion, such as... Figure 1 As shown, it specifically includes:

[0062] A three-dimensional geomechanical numerical model of the target geological hazard body is established. The numerical model defines the physical framework of hazard evolution based on the constitutive relationship of rock and soil.

[0063] Construct a multi-source heterogeneous data stream, which includes space-based InSAR deformation data, airborne UAV lidar point cloud data, ground-based GNSS displacement data, and deep underground stress monitoring data.

[0064] Using the Physical Information Neural Network (PINN) as a surrogate model, the residuals of the partial differential equations of the numerical model are added as constraints to the loss function. By minimizing the difference between the monitoring data and the model predictions, the unsteady mechanical parameters in the numerical model are inverted and updated in real time. The parameters include at least time-varying cohesion, internal friction angle and permeability coefficient, thus obtaining a digital twin that matches the real state of the current geological body in real time.

[0065] Based on the updated digital twin, combined with short-term weather forecast data for future periods, forward dynamic extrapolation at multiple time steps is performed to calculate the probability of the cumulative plastic zone and potential sliding surface connecting at future moments.

[0066] Based on the simulation results, the Lyapunov index spectrum of the state is calculated. When the maximum Lyapunov index turns from negative to positive, it is determined that the geological body has entered a chaotic and unstable state, and a disaster risk warning signal is generated accordingly.

[0067] Furthermore, the three-dimensional geomechanical numerical model takes the entire spatial domain of the target geological hazard body as the analysis object. Based on high-density geological survey profile data, indoor soil and rock mechanics test results, structural plane occurrence and distribution data, and aquifer spatial distribution data, it completes high-precision spatial discretization. The entire model domain is divided into spatial units and spatial nodes, and independent calculation partitions are divided according to soil and rock type, structural plane characteristics, weak interlayer distribution, and pore structure characteristics. Each partition is independently configured with constitutive relations, stress-strain correlation rules, seepage conduction rules, and energy transfer rules. The model is supported by the coupled theories of continuum mechanics, soil mechanics, rock mechanics, and seepage mechanics. It unifies the deformation, damage, seepage, stress redistribution, and energy dissipation of the geological body within the same physical framework. Through the constitutive relations of soil and rock, it fully defines the full-cycle evolution path from elastic deformation, plastic development, local damage, structural expansion to overall instability, forming a closed and self-consistent physical calculation system.

[0068] All field quantities within the model, including stress, strain, displacement, pore water pressure, velocity, and damage, are transmitted and iterated within a unified topological grid. Rigid mechanical and flexible seepage correlations are established between grid nodes. Rigid mechanical correlations ensure that stress, strain, and displacement satisfy equilibrium and coordination conditions between nodes, while flexible seepage correlations ensure that pore water pressure, velocity, and flow rate satisfy continuity and conservation conditions between nodes. External loads, fluid seepage, and parameter changes can be uniformly transmitted and accurately responded to throughout the entire space. The state update of each node synchronously affects the state calculation of surrounding neighboring units, ensuring that changes in field quantities throughout the entire space are lag-free, distortion-free, and misaligned.

[0069] Furthermore, the multi-source heterogeneous data stream consists of four types of observation data: space-based, airborne, ground-based, and deep underground. Space-based InSAR deformation data covers the entire spatial range, outputting high-resolution planar deformation field time series data. Each time series contains a large number of deformation observation points, reflecting the large-scale deformation trend and local abnormal deformation areas of the geological body. Airborne UAV lidar point cloud data outputs three-dimensional geometric information of the surface with high density, acquiring massive point cloud data per sortie, accurately depicting topographic relief, slope curvature, micro-topographic changes, and spatial distribution of cracks. Ground-based GNSS displacement data outputs three-dimensional displacement information of key points at high frequency, deploying multiple monitoring stations across the entire region, acquiring a large amount of displacement time series data per day, capturing the dynamic deformation characteristics and deformation rate of the geological body. Deep underground stress monitoring data outputs stress time series at multiple depths and orientations within the geological body, acquiring a large amount of stress data per measurement point per day, reflecting the internal stress state and load adjustment process.

[0070] All data undergoes timestamp alignment, spatial coordinate transformation, and observation dimension normalization to eliminate differences in time sampling, spatial representation, and numerical scale among heterogeneous data, achieving a unique representation of the same physical quantity at the same time and location. Time alignment employs linear interpolation and resampling to unify all data to the same time interval. Spatial coordinate transformation achieves seamless multi-coordinate system conversion. Observation dimension normalization maps displacement, deformation, stress, and point cloud geometric features to the same numerical range, eliminating computational interference caused by dimensional differences. Based on this, a Kalman filter fusion algorithm is used to generate a high-confidence observation sequence. The state update equation satisfies the optimal estimation theory, with the specific formula as follows: , This indicates that after fusion processing, at time... The obtained high-confidence observation sequence is generated by combining all multi-source heterogeneous data after normalization. The data dimension is completely matched with the dimension of the monitored physical quantity. The sequence length covers the entire monitoring period, and the total number of data points meets the input requirements of subsequent models. The vector represents the true state of the geological body, containing 12 dimensions of state information, including displacement, stress, deformation gradient, and seepage gradient components. Each dimension corresponds to the actual physical state of the geological body, and the vector dimension is consistent with the degree of freedom of the model nodes. The observation matrix uses a fixed dimension of 12 rows and 6 columns to complete the linear mapping from the high-dimensional real state vector to the low-dimensional observation vector. The matrix elements are determined by the deployment location of the observation equipment, the monitoring principle, and the spatial response function. The matrix values ​​are calibrated and kept fixed after the spatiotemporal reference is unified. The noise is a 6-dimensional observation noise that strictly follows a zero-mean Gaussian distribution. The noise components include equipment noise, environmental interference noise, and data conversion error. The distribution characteristics are determined by statistical analysis of measured data. The noise covariance matrix is ​​6×6. The diagonal elements represent the dispersion of noise in each dimension, and the off-diagonal elements represent the statistical correlation between noise in different dimensions. The matrix value is obtained by long-term statistical calculation of the original monitoring data, which can truly reflect the actual distribution characteristics of the observed noise. The fusion process is recursively updated through this formula to continuously correct observation bias, suppress random noise and outlier interference, and transform the original multi-source data into a standardized observation set with high consistency, high smoothness and high reliability.

[0071] Furthermore, such as Figure 2 As shown, the Physical Information Neural Network (PINN) is based on a multi-layer fully connected network architecture. The network input layer contains 18 input dimensions, including spatial three-dimensional coordinates, time variables, boundary condition vectors, and external load vectors. The network hidden layer contains 12 layers of neurons, with the number of neurons in each layer being 512, 256, 256, 128, 128, 64, 64, 32, 32, 16, 16, and 8, respectively. A layer-by-layer shrinking structure is used to improve the accuracy and convergence stability of high-dimensional mechanical field mapping. The network output layer simultaneously outputs 24-dimensional physical fields, including displacement field, stress field, pore water pressure field, and mechanical parameter field, realizing high-dimensional nonlinear mapping and fast approximation of geomechanical systems. The network activation function adopts a hybrid mechanism of Swish and Tanh, with the first 6 layers using Swish. Activation functions enhance nonlinear expressiveness; the last six layers use the Tanh activation function to ensure that the output physical quantities are bounded and conform to mechanical constraints. Xavier uniform distribution is used for network initialization to avoid gradient vanishing and gradient exploding. During network training, the partial differential equations, boundary conditions, and initial conditions of the numerical model are embedded as rigid soft constraints into the optimization objective, no longer relying solely on data fitting, ensuring that the model output simultaneously satisfies data matching accuracy and consistency with physical laws. The training batch size is set to 256, the optimizer is AdamW, the initial learning rate is 1e-4, decays by a factor of 0.8 every 200 training epochs, the weight decay coefficient is 1e-5, and the total number of training epochs is 10,000. Training is stopped early when the total loss function drops below 1e-6. The total loss function is constructed as follows: Due to data mismatch items Physical equation residuals Boundary condition residuals The weighted combination, with the three components constraining the network's convergence direction according to preset weights, is defined by the following formula: The formula is the core loss function for training the PINN model, which consists of three weighted parts, each of which corresponds to a specific physical constraint and data matching target; The total loss function is the sole optimization objective for network training. The smaller the value, the higher the degree to which the model output simultaneously satisfies data matching and physical constraints. , , These are weighting coefficients used to balance the contribution of data fitting and the contribution of physical constraints. The weight values ​​were determined through multiple experiments to ensure that physical constraints are not covered by data fitting during training. For neural networks in spatial coordinates With time The output displacement prediction value contains three-dimensional displacement components, and the value is directly calculated by the network forward propagation; The observed displacement values ​​are obtained after fusing multi-source heterogeneous data, and their dimensions are perfectly matched with the predicted values. The values ​​are directly taken from the observation sequence after Kalman filtering fusion. The total number of observation data points is evenly distributed across the entire model space and the entire time series to ensure data coverage integrity. For the data mismatch term, the deviation between the network-predicted displacement and the measured displacement is calculated using the mean square error form, which characterizes the degree of fit of the model to the real monitoring data.

[0072] The residual operator for the geomechanical governing equations is calculated by combining network output with physical partial differential equations. A total of 384,000 points are configured for the physical equations, which are randomly and uniformly distributed throughout the model's global space to ensure that physical constraints cover the entire domain. The residual term of the physical equation represents the degree to which the model output deviates from the geomechanical control equation; The boundary condition operator includes multi-dimensional constraints such as displacement boundary, stress boundary, and seepage boundary. The values ​​are obtained by comparing the network output with the preset boundary conditions. The total number of boundary condition points is 24,000, which are densely distributed along the entire boundary of the model to ensure the effectiveness of boundary constraints.

[0073] The boundary condition residual term represents the degree to which the model output deviates from the boundary conditions; each term of the loss function represents the degree of deviation in the form of mean square error. Through iterative optimization, the total loss value is continuously reduced, so that the network output can simultaneously approximate the observed data and the physical equations.

[0074] Physical equation residual terms Based on Biot's consolidation theory, the momentum and mass conservation equations in the seepage-stress coupling process of soil and rock are uniformly transformed into residual forms, ensuring that the network output strictly satisfies the coupling laws of soil and rock mechanics and seepage. The specific formula for the residual operator is as follows:

[0075]

[0076] The formula is the core residual equation for the physical constraints, consisting of two lines, corresponding to the two major physical laws of conservation of momentum and conservation of mass, respectively; the first line corresponds to the residual of the momentum conservation equation for the soil and rock skeleton. The divergence of the effective stress tensor characterizes the rate of change of stress in the rock and soil skeleton in space. The effective stress tensor contains six independent components, including three-dimensional normal stress and three-dimensional shear stress. Each component is calculated by combining the network output with the constitutive relation. The value of is the Biot coefficient, and its range is constrained within the physically feasible region. Its value is determined by the pore structure and compressibility characteristics of the soil and rock mass. The pore water pressure gradient represents the rate of change of pore water pressure in space and is the core driving force for the movement of pore fluids. It is a gravitational volume force. Density of the rock and soil mass This is the vector of gravitational acceleration, directed vertically downwards; For inertial force, The first line represents the second derivative of the displacement vector with respect to time, characterizing the inertial effect of the geological body's motion, which cannot be ignored during dynamic deformation; the second line corresponds to the residuals of the pore fluid mass conservation equation. The rate of change of pore fluid mass over time. Porosity reflects the proportion of pore space in soil and rock masses. Pore ​​fluid density; The divergence of fluid mass flux characterizes the balance between fluid inflow and outflow in space. This is the pore fluid velocity vector, which includes three-dimensional velocity components; The fluid source and sink term characterizes the intensity of external fluid injection and outflow, and its value is determined by meteorological supply and boundary seepage conditions. This residual term transforms the laws of continuum mechanics and seepage mechanics into computable soft constraints, enabling the surrogate model output to have strict physical consistency and interpretability, with the residual mean controlled within 1e-5. The network learns the data distribution and physical laws simultaneously during training, without relying on a large number of labeled samples, and still maintains stable generalization ability under small sample conditions.

[0077] Furthermore, such as Figure 3 As shown, the unsteady mechanical parameters include time-varying cohesion. Time-varying internal friction angle Time-varying permeability coefficient The three types of parameters exhibit continuous spatiotemporal evolution characteristics as geological bodies deform, accumulate damage, and change seepage states. These parameters are parameterized as high-dimensional implicit functions that vary with time, spatial coordinates, and mechanical state, embedded between layers 6 and 7 of a physical information neural network, and updated synchronously with forward and backward propagation. The parameterized function is constructed using a multilayer perceptron, containing three hidden layers with 32, 16, and 8 neurons per layer, respectively. Inputs include spatial coordinates, time, equivalent plastic strain, and stress triaxiality; outputs are time-varying mechanical parameters at the corresponding locations, achieving a strong coupling relationship between parameters and mechanical state. Automatic differentiation technology using neural networks is employed to calculate higher-order partial derivatives of the stress tensor with respect to the strain tensor, mechanical parameters, and spatial coordinates, enabling automatic differentiation. A combination of forward and backward propagation is employed. The first derivative is used for gradient calculation, the second derivative for calculating the residuals of the physical equations, and the third derivative for sensitivity analysis. The computational accuracy reaches machine-level precision with no numerical discretization errors. The sensitivity matrix of the mechanical field output to parameter changes is obtained, perfectly matching the number of nodes and parameter dimensions, accurately reflecting the impact of parameter perturbations on the model output. The sensitivity values ​​remain stable within a controllable range. The adjoint state method is used to solve the gradient vector of the total loss function with respect to the unsteady mechanical parameters. This gradient vector contains the gradient components of all spatially partitioned parameters. The complexity of solving the adjoint equation is consistent with forward propagation and does not significantly increase with the number of parameters. This guides the iterative update of parameters along the loss descent direction. The parameter update formula is:

[0078]

[0079] This formula is the core equation for the iterative update of unsteady mechanical parameters, enabling adaptive optimization of parameters and constraints on the physical feasible region. The set of unsteady mechanical parameters to be inverted covers 18 partitions with a total of 54 independent parameters. Each parameter corresponds to the time-varying mechanical properties of a computational partition in the model. For the first The parameter vector updated in the next iteration is the output of this iteration and will be directly used to update the mechanical parameters of the numerical model. For the first The parameter vector of the next iteration is the basis of the input for this iteration, and retains the parameter state of the previous iteration; The learning rate is dynamically adjusted within the range of 1e-5 to 1e-3 to control the single-step update amplitude of the parameters. If the learning rate is too large, it will easily lead to parameter oscillation, while if it is too small, the convergence will be slow. The total loss function at the th At the next iteration, the parameters are... The gradient vector is solved by the adjoint state method and automatic differentiation technique. The gradient direction points to the direction in which the total loss function decreases the fastest, ensuring that the parameter update direction is optimal. The momentum coefficient is fixed at 0.9 and is used to accumulate the update direction of the previous iteration, improve the stability of the iteration convergence, suppress oscillations, and accelerate the model convergence. For the momentum term, the trend of parameter updates in the previous step is retained to avoid the training process getting stuck in local optima; The regularization coefficient is fixed at 1e-4 and is used to suppress overfitting and abnormal fluctuations in parameter updates, thereby improving the generalization ability and physical rationality of the parameters. As a projection operator, the updated parameters are constrained within a feasible region that conforms to the geophysical properties. Internally, the projection rules are determined based on geotechnical mechanical test data, ensuring that the parameter values ​​do not exceed the physically reasonable range and have clear engineering significance. The parameter update process continuously reduces the multidimensional deviation between model predictions and measured data. The average deviation of displacement prediction is less than 1e-3, the average deviation of stress prediction is less than 1e-2, and the average deviation of pore water pressure prediction is less than 1e-2. This ensures that the mechanical state, seepage state, and damage state of the numerical model are kept in real-time synchronized with the actual state of the geological body, with a synchronization delay of less than 1 second. This forms a digital twin that can dynamically evolve, be mapped in real time, and be adaptively adjusted. The digital twin and the physical entity are completely matched in terms of spatial structure, mechanical properties, and evolutionary trends. The parameter inversion process does not require manual intervention and can run independently on edge computing terminals, meeting the low power consumption and high real-time requirements of field monitoring scenarios.

[0080] Furthermore, using the updated digital twin as the computational foundation, short-term weather forecast data for the next 120 time periods are input. Meteorological elements are transformed into mechanical excitations such as pore water pressure fields, external load fields, and fluid recharge fields. The excitation field data covers the entire simulation cycle and 47,260 nodes across the entire domain. The excitation intensity changes continuously over time, and its spatial distribution highly matches the permeability characteristics and topographic features of the geological body. Multi-time-step forward dynamic simulations are conducted, employing an adaptive mechanism: the time step is increased during stable phases and decreased during deformation acceleration phases. The future mechanical and seepage responses of the geological body are solved time-by-time, reconstructing the entire evolution process of the geological body's future state. The simulation process solves the dynamic equilibrium equations considering the coupling effect of pore water pressure and rock-soil skeleton stress. The equations are numerically solved using a combination of spatial finite element discretization and time-step discretization, resulting in the following discretized form:

[0081]

[0082] in, The system quality matrix is ​​determined by the density of the soil and rock mass and the spatial discrete grid. The matrix adopts a sparse storage format to reduce memory usage and improve solution efficiency. for The acceleration vector at each time point is obtained from the second derivative of displacement with respect to time, and it characterizes the dynamic acceleration properties of the geological body. The system damping matrix is ​​represented by Rayleigh damping, which characterizes the energy dissipation of the system. The matrix value is determined by the damping of the soil and rock material and the structural damping. for The velocity vector at each time point is obtained from the first derivative of displacement with respect to time and represents the speed of geological body movement. The system stiffness matrix has dimensions of 47260×47260 and varies with unsteady mechanical parameters. Dynamically updated, reflecting the deterioration process of the mechanical properties of geological bodies in real time, is the core matrix that embodies the dynamic characteristics of digital twins; for The displacement vector at each time point, which contains three-dimensional displacement components, is the core output physical quantity of the dynamic deduction. for The external load vector at any time includes self-weight load, meteorological load, and boundary load. The load values ​​change dynamically with time and space. The equivalent volumetric load generated by pore water pressure. The pore water pressure field is derived from meteorological data. For a shape function matrix, an eight-node isoparametric element form is used to achieve interpolation mapping from nodal variables to variables within the element. The pore water pressure gradient is a complete representation of the coupling effect of the seepage field on the mechanical field.

[0083] To balance computational stability, accuracy, and efficiency, a hybrid explicit-implicit integration scheme is employed to solve the dynamic and seepage terms in a coupled manner. The dynamic response term is explicitly integrated using the central difference method to quickly capture dynamic deformation and impact characteristics; the seepage diffusion term is implicitly integrated using the backward difference method to ensure computational stability and convergence over large time steps. The coupled solution scheme comprises two formulas; the first is the formula for solving the dynamic term:

[0084]

[0085] in, For those in demand The displacement vector at each time point is the output of the integration in this step; For known The displacement vector at each time point is provided by the result of the previous iteration. For the time step, an adaptive adjustment strategy is adopted: the step size is increased to improve efficiency during the stable phase, and the step size is decreased to ensure accuracy during the deformation acceleration phase. For known The velocity vector at each time point is provided by the result of the previous iteration; The coefficients related to the time step are determined by the discretization scheme of the central difference method; The inverse of the mass matrix is ​​calculated quickly using sparse matrix solving techniques, avoiding the computational overhead of direct inversion. For known The external node load vector at any given time includes all external excitations; For known The resistance vector of internal nodes at any given time is determined by the stress and deformation state of the soil and rock skeleton, reflecting the soil and rock's ability to resist deformation; The unbalanced force at the nodes is the core driving force that causes displacement of geological bodies;

[0086] The second is the formula for solving the seepage term:

[0087]

[0088] in, For those in demand The pore water pressure vector at each time point is the core output for solving the seepage field. For known The pore water pressure vector at each time point is provided by the result of the previous iteration. The time step is kept consistent with the dynamic integral step to ensure spatiotemporal coupling consistency. It is the inverse of the permeability matrix, representing the seepage and conduction characteristics of the soil and rock mass. The matrix value is dynamically updated with the permeability coefficient. For known The external fluid flow vector at any given time is determined by meteorological recharge and boundary conditions; The fluid-structure interaction matrix characterizes the interaction strength between the displacement field and the pore water pressure field. The matrix value is determined by the relationship between the deformation of the rock and soil skeleton and the change of pore structure. This term represents the coupling effect of the displacement field and the seepage field, reflecting the squeezing effect of skeletal deformation on pore fluids. By solving the coupled formulas, the dynamic field and the seepage field are synchronously iterated, yielding full-dimensional field information on stress, displacement, seepage, and damage of the geological body at future moments. The field data covers the entire spatial and temporal simulation process, with a single-moment output of 1,134,240 values. The simulation results can be visualized, intuitively displaying the deformation trend of the geological body, the damage propagation path, and the seepage field distribution.

[0089] Furthermore, based on the plastic strain time series data and stress state time series data obtained from forward dynamics deduction, the global spatial local damage variables... This variable comprehensively reflects the degree of plastic damage in the soil and rock mass, the influence of stress state, and the rate of deformation development, and is determined by cumulative plastic strain. With stress triaxiality Co-driven evolution, the evolutionary equation is: , Spatial location ,time The local damage variable ranges from 0 to 1. The larger the value, the more severe the damage to the soil and rock mass and the more significant the structural destruction. It is a natural exponential function, which ensures that the damage variable is monotonically increasing and bounded, which is consistent with the physical property of irreversible damage in soil and rock masses. The damage accumulation integral term covers the entire evolution process from the initial time to the current time. and This is the material damage constant, determined by geotechnical mechanics tests. Its value is fixed and physically unique. Controlling the overall rate of damage development, The extent to which the triaxiality of the control stress affects the damage; for The equivalent plastic strain rate at any given time characterizes the rate of plastic deformation development, and its value is directly calculated from the strain results derived from dynamics. for The stress triaxiality at any time reflects the three-dimensional stress state of the rock and soil mass. The more unfavorable the stress state, the faster the damage develops. The stress triaxiality influence term amplifies the damage development rate under unfavorable stress conditions; the local damage variable forms a continuous and smooth damage field in space, and the damage value increases monotonically with deformation and stress development;

[0090] Based on the distribution of local damage variables across the entire geological body, the probability of the potential sliding surface connecting across the entire geological body is calculated. The possibility of penetration is determined by comparing the path integral of the damage field with the critical damage threshold. The formula is as follows: , This represents the probability of potential slip surface connection across the entire geological body, with a value ranging from 0 to 1. The closer the value is to 1, the higher the risk of connection instability. It is a probabilistic computation operator, determined based on the statistical results of potential interconnected paths across the entire domain; Indicates that in the computation of the entire domain There is at least one potential through path in memory. A total of 1280 candidate paths were generated across the entire domain, covering all sliding directions and structural surface orientations; Potential connecting routes The damage variable space integral represents the average damage degree of the entire path, and the integration interval is the spatial range of the entire through path; The critical damage threshold is determined by the mechanical properties of soil and rock and engineering experience. When the path damage integral exceeds this threshold, the path is judged to have failed through. When any potential path to fail meets the condition of exceeding the damage integral limit, the risk of failure through is included. The probability of failure through the entire domain is derived by combining the risks of all paths. It can quantitatively identify the damage concentration area, the direction of the sliding surface expansion, the level of failure through risk, and the overall instability probability, providing spatial distribution basis and quantitative indicators for the determination of the state of imminent disaster.

[0091] Furthermore, based on the displacement time series, velocity time series, and mechanical parameter time series obtained from dynamic derivation, a high-dimensional state space and reconstructed phase space vector of the geological body system are constructed. This transforms the continuous evolution process of the geological body into trajectory motion in a high-dimensional phase space. The formula for the phase space vector is: , For a moment The reconstructed phase space vector, with a total dimension of 9, is the fundamental object of nonlinear dynamics analysis. for The displacement state vector at any time contains three-dimensional displacement components, representing the spatial deformation state of the geological body. for The velocity state vector at any given time contains three-dimensional velocity components, representing the motion state of the geological body; for The state vector of unsteady mechanical parameters at any given time includes time-varying cohesion, internal friction angle, and permeability coefficient components, characterizing the material degradation state of the geological body; This is a vector transpose operation, which converts a row vector into a column vector to meet the format requirements of subsequent matrix operations;

[0092] Apply a small initial perturbation to the phase space trajectory Solve the system variational equation to obtain the disturbance propagation law and amplification rate. The variational equation is:

[0093]

[0094] in, perturbation vector The derivative with respect to time characterizes the rate of change of the disturbance with time; To provide a small initial disturbance applied to the phase space trajectory, the disturbance amplitude is set to 1e-8 to ensure that the disturbance is small enough and does not change the essential characteristics of the system. For the system in state The Jacobian matrix at the position has a dimension of 9×9. The matrix elements are precisely calculated by automatic differentiation technology. Each element represents the sensitivity of one state variable to another, reflecting the coupling strength and disturbance amplification characteristics between state variables. The product of the Jacobian matrix and the perturbation vector represents the driving effect of the system's own dynamic characteristics on the perturbation. By solving this equation, the evolution law of the perturbation vector over time can be obtained, providing a basis for the calculation of the Lyapunov exponent.

[0095] The Gram-Schmidt orthogonalization method is used to periodically reorthogonalize the disturbance vectors of the discrete time series. This orthogonalization operation is performed every 5 time steps to suppress computational errors caused by disturbance divergence and ensure the accuracy of the exponent calculation. The maximum Lyapunov exponent is calculated based on the divergence rate of the disturbance trajectory, using the following formula: The formula is the core calculation equation for the maximum Lyapunov exponent and serves as the quantitative basis for determining chaotic instability. The maximum Lyapunov exponent is used, and the sign of the value directly determines the stability: a negative number indicates stable convergence, zero indicates critical stability, and a positive number indicates chaotic divergence. For the limit operation when time approaches infinity, it represents the long-term evolutionary characteristics; The reciprocal of time is used to normalize the average divergence rate of the perturbation; The natural logarithm operation transforms the exponential growth of the perturbation into linear growth, which facilitates numerical calculation. For a moment The L2 norm of the perturbation vector characterizes the overall magnitude of the perturbation; The L2 norm of the initial perturbation vector represents the magnitude of the initial perturbation. This represents the amplification factor of the disturbance amplitude, reflecting the degree of divergence of the disturbance over time; continuous calculation is performed using a sliding time window. The time series curve, with a step size of 1 time, outputs the maximum Lyapunov exponent value point by point; when the maximum Lyapunov exponent turns from negative to positive and crosses the zero value criterion, and the second derivative of the curve is greater than zero, it is determined that the geological body has entered a chaotic unstable state from a stable convergent state.

[0096] Set tiered early warning thresholds With duration threshold ,when And duration At that time, a Level 1 emergency disaster warning is triggered; the fractal dimension is calculated simultaneously. The box counting method is used for calculation, the calculation interval covers the entire phase space, and the number of boxes increases exponentially. Furthermore, when the fractal dimension exhibits non-integer dimension characteristics, it is determined that the geological body has entered a chaotic evolution stage dominated by strange attractors, and the instability is irreversible. Based on this, a standardized disaster risk warning signal is generated to complete the entire process of warning judgment. The entire judgment process takes less than 3 seconds, meeting the real-time requirements of disaster warning. This judgment method breaks away from the dependence on traditional empirical thresholds and is based on nonlinear dynamics theory, making the warning results more scientific and reliable.

[0097] This embodiment combines multi-source heterogeneous data fusion, three-dimensional geomechanical modeling, and physical information neural networks to achieve real-time inversion of unsteady mechanical parameters and dynamic matching of digital twins. It relies on explicit and implicit hybrid integrals to complete multi-time-step fluid-structure coupling inference, accurately calculate damage distribution and sliding surface penetration probability, and quantitatively identify chaotic instability critical states based on Lyapunov exponent spectrum. This can significantly improve the accuracy, timeliness, and reliability of geological disaster early warning, and achieve efficient early warning that is physically reliable, dynamically adaptive, and quantitatively predictable throughout the entire process.

[0098] Based on Embodiment 1, Embodiment 2 details the technical verification of this application. Currently, multiple technical routes have been formed in the field of geological disaster early warning. Different methods have significant differences in data utilization, physical mechanism expression, parameter update mechanism, and early warning criterion design. In order to objectively illustrate the comprehensive performance advantages of this application, six advanced early warning methods that have been disclosed and widely used in the industry are selected for comparative analysis, namely: the geological early warning method based on physical information neural network (PINN-EW), the slope early warning method based on digital twin (DT-EW), the numerical simulation method based on fluid-structure interaction (FSC-NS), the statistical learning method based on multi-source data fusion (MDF-SL), the chaotic identification method based on Lyapunov exponent (LE-CI), and the parameter identification method based on finite element inversion (FEI-PI). The three-dimensional geological model used includes spatial units and spatial nodes, and is divided into 18 geotechnical mechanics zones. Multi-source monitoring data are input, including InSAR deformation data, UAV LiDAR point cloud, GNSS displacement time series, and deep underground stress data.

[0099] To quantify the overall performance of different methods, six key indicators were statistically analyzed: goodness of fit, average error of parameter inversion, time consumption of single-step inference, accuracy of early warning, early warning lead time, and physical consistency. As shown in Table 1, the core indicators of the seven methods are compared quantitatively.

[0100] Table 1. Quantitative Comparison of Core Indicators for 7 Methods

[0101] This application 99.2 1.5 2.1 96.7 82 98.1 PINN-EW 94.6 3.8 3.1 89.3 51 91.2 DT-EW 93.1 4.2 3.5 87.1 47 89.7 FSC-NS 89.7 5.6 4.8 83.6 39 85.3 MDF-SL 85.2 6.9 2.0 81.2 33 72.4 LE-CI 88.4 6.1 2.2 85.9 45 78.6 FEI-PI 90.5 5.3 4.2 84.8 41 86.9

[0102] Note: FD fit; PIE parameter inversion average error; ST single-step extrapolation time; WA early warning accuracy; WLE early warning lead time; PC physical consistency.

[0103] All values ​​in Table 1 are the average values ​​of five independent experiments under 120 consecutive simulations, demonstrating statistical stability. The results show that this application outperforms the six control methods in all six indicators, achieving a goodness of fit of 99.2%, highly reproducing the true mechanical state of the geological body. The average error of parameter inversion is only 1.5%, far lower than other methods. The single-step simulation time remains at 2.1 seconds, balancing accuracy and real-time performance. The early warning accuracy rate is 96.7%, and the early warning lead time is 82 minutes, providing ample time for emergency response. The physical consistency reaches 98.1%, ensuring that the calculation results conform to the basic laws of geotechnical mechanics. Overall, the application significantly outperforms existing technologies.

[0104] The fusion quality of multi-source heterogeneous data directly determines the accuracy of subsequent model inversion and inference. To evaluate the processing capabilities of different methods for space-based, air-based, ground-based, and deep data, four indicators, namely spatiotemporal uniformity rate, noise suppression rate, effective data rate, and observation consistency, are compared. Table 2 shows the comparison of multi-source data fusion effects.

[0105] Table 2 Comparison of Multi-Source Data Fusion Effects

[0106] This application 100 94.3 98.6 97.8 PINN-EW 97.2 86.5 95.1 92.3 DT-EW 96.8 85.7 94.2 91.5 FSC-NS 95.3 83.1 92.6 89.7 MDF-SL 94.1 80.2 91.3 87.2 LE-CI 95.7 82.6 92.9 88.4 FEI-PI 96.1 83.9 93.5 89.1

[0107] Note: STR spatiotemporal consistency rate; NSR noise suppression rate; EDR effective data rate; OC observation consistency.

[0108] Table 2 shows that, based on the spatiotemporal benchmark unification algorithm and Kalman filter fusion strategy, this application achieves a spatiotemporal unification rate of 100%, a noise suppression rate of 94.3%, an effective data rate of 98.6%, and an observation consistency of 97.8%. All indicators are significantly higher than the control method, indicating that this application can effectively eliminate the differences in time, space, and numerical scale of multi-source data, suppress observation noise, improve data reliability, and provide high-quality input for high-precision modeling.

[0109] In continuous time-series simulations, the degree of matching between the model output and the actual mechanical state of the geological body can intuitively reflect the stability of the method. Based on the comprehensive calculation results of global nodal displacement, stress, and pore water pressure, such as... Figure 4 As shown, the model-measured mechanical state fitting degree comparison was obtained. Throughout the entire time series, the curve corresponding to this application consistently remained above 99.2%, with a smooth curve and no significant fluctuations. It remained highly stable even in the later stages of the simulation. In contrast, PINN-EW and DT-EW began to show a slight decrease after 60 steps. FSC-NS and FEI-PI, due to their fixed parameters, could not adapt to changes in the geological body state, resulting in a more significant decrease. MDF-SL, relying entirely on data-driven analysis and lacking physical constraints, exhibited the largest curve fluctuations. LE-CI, without incorporating a mechanical model and only performing time series analysis, showed an overall low matching degree. This result fully demonstrates that this application, by embedding physical partial differential equations into loss function constraints, real-time inversion of time-varying mechanical parameters, and dynamic updating of the digital twin, can continuously maintain a high degree of synchronization between the model and the real geological body state. This avoids the fitting degree decay problem caused by parameter solidification, missing mechanisms, and data noise in traditional methods, and maintains high-precision matching capability even in long-term, strongly coupled, and unsteady geological evolution processes.

[0110] The strength of physical constraints directly determines the reliability of the early warning results. To quantify the differences in the degree of compliance with mechanical laws among different methods, four indicators are used for comparison: mean PDE residual, mean boundary condition residual, damage calculation deviation, and penetration probability error. Table 3 shows the comparison of the effectiveness of physical constraints.

[0111] Table 3 Comparison of the effectiveness of physical constraints

[0112] This application 8.7 1.2 0.021 0.024 PINN-EW 15.2 2.4 0.034 0.037 DT-EW 16.8 2.7 0.039 0.041 FSC-NS 21.3 3.6 0.045 0.048 MDF-SL 42.7 6.9 0.072 0.075 LE-CI 28.6 4.2 0.051 0.053 FEI-PI 20.5 3.3 0.042 0.045

[0113] Note: PDE-RM residual mean; BC-RM boundary residual mean; DC-D damage calculation bias; CP-E penetration probability error.

[0114] In Table 3, smaller values ​​indicate stricter physical constraints and higher computational accuracy. The average PDE residual of this application is only 8.7 × 10⁻⁶. -6 The mean boundary residual is 1.2 × 10⁻⁶. -5 The damage calculation deviation was 0.021 and the penetration probability error was 0.024, both significantly lower than other comparative methods. This indicates that the soft constraints constructed in this application based on Biot consolidation theory and the equations of conservation of momentum and mass can strictly control the model output and ensure that the evolution of displacement field, stress field, seepage field and damage field conforms to the real physical laws. Other methods either lack strong physical constraints, fail to achieve fluid-structure interaction, or fail to embed partial differential equations, resulting in low physical consistency and large errors in the calculation results.

[0115] Computational efficiency is a key indicator for the successful implementation of disaster early warning systems, such as... Figure 5 The curves showing the time consumption comparison of multi-time-step dynamics derivation cover the entire cycle from 1 to 120 steps. This application adopts a hybrid solution format of explicit integration of dynamic terms and implicit integration of seepage terms. The single-step time is stable at about 2.1 seconds throughout the entire derivation process, with a total time of 252 seconds. PINN-EW, due to the lack of a digital twin real-time synchronization mechanism, has a single-step time of 3.1 seconds and a total time of 368 seconds. DT-EW, without using a surrogate model for acceleration, has a single-step time of 3.5 seconds and a total time of 415 seconds. FSC-NS and FEI-PI use traditional finite element solutions, which have large computational loads and long time consumption, with single-step times reaching 4.8 seconds and 4.2 seconds, respectively. Although MDF-SL has a slightly lower single-step time, it completely lacks physical constraints, and the results do not have engineering credibility. LE-CI only performs time series analysis without conducting mechanical derivation, and cannot provide complete early warning basis. The curve clearly shows that this application achieves a significant improvement in computational efficiency while ensuring high accuracy and strong physical constraints by replacing the traditional numerical model with a surrogate model, improving solution efficiency through explicit and implicit hybrid integration, and optimizing the computational process with adaptive step size. This fully meets the real-time requirements for geological disaster early warning.

[0116] Early warning accuracy and early warning lead time are the core performance indicators of early warning methods. Using a preset simulated instability moment as a benchmark, the early warning results of different methods are statistically analyzed, such as... Figure 6 The comparison of the accuracy and lead time of the disaster early warning shown in the figure reveals that the accuracy of this application reaches 96.7%, with an average early warning lead time of 82 minutes, far exceeding the other six existing technologies. PINN-EW has an accuracy of 89.3% and a lead time of 51 minutes, DT-EW has an accuracy of 87.1% and a lead time of 47 minutes, while FSC-NS, MDF-SL, LE-CI, and FEI-PI all have lower accuracy and lead times. This result stems from the fact that this application combines multi-field coupled extrapolation, damage penetration probability calculation, Lyapunov exponent spectrum, and fractal dimension joint criteria, enabling it to accurately capture the critical point of geological bodies transitioning from stable evolution to chaotic instability. Traditional methods either lack mechanical extrapolation, damage analysis, or nonlinear dynamic criteria, making it difficult to accurately identify critical instability states, leading to delayed early warnings or frequent false alarms. This fully demonstrates the significant advantages of this application in critical identification and early warning reliability.

[0117] This embodiment, through multi-dimensional comparison and verification with mainstream advanced early warning methods, fully demonstrates the comprehensive advantages of this application in data fusion, physical constraints, parameter inversion, computational efficiency, and early warning performance. Its model-measured state fitting degree, parameter inversion accuracy, single-step inference time, early warning accuracy, and lead time are all significantly better than existing technologies. It verifies the synergistic effect of multi-source data processing, PINN physical constraints, and digital twin update mechanism, and exhibits higher computational stability and physical consistency.

[0118] Example 3: This verification was conducted based on typical slope geological conditions. The same three-dimensional geological model structure and mesh generation method as in Example 2 were used. The model includes complete stratigraphic structure, weak interlayers, aquifer distribution and boundary load conditions. The input data comes from multi-source heterogeneous information collected by the actual monitoring system. After spatiotemporal alignment, noise suppression and standardization, it is input into the calculation process. All parameter settings, formula calls and iteration processes are executed according to the method of this application. No manual correction or empirical adjustment is introduced in the calculation process, and the original running state of the method is completely maintained.

[0119] The physical information neural network was constructed according to the design structure of this application. The input layer contains 18 dimensions, including three-dimensional spatial coordinates, time, boundary conditions, and external loads. There are 12 hidden layers with the following neuron counts: 512, 256, 256, 128, 128, 64, 64, 32, 32, 16, 16, and 8. The output layer synchronously outputs 24 physical quantities, including displacement, stress, pore water pressure, and mechanical parameters. The activation function uses a hybrid approach of Swish in the front layer and Tanh in the back layer. The network is initialized using a uniform distribution of Xavier signals. The training batch size is 256. The optimizer is AdamW with an initial learning rate of 1e-4, decaying by a factor of 0.8 every 200 epochs, for a total of 10,000 training epochs. The final total loss function converges stably to 9.4 × 10⁻⁻⁻⁶. 7 This meets the requirements for high-precision calculations.

[0120] To comprehensively reflect the performance of the method in actual operation, the indicators of key links such as model construction, data fusion, network training, parameter inversion, twin synchronization, dynamics deduction, and early warning output are uniformly summarized, forming Table 4, which shows the full-process verification results of the method in this application.

[0121] Table 4. Validation results of the entire process of the method in this application.

[0122] MCA 98.7% TSR 99.3% DFQ 97.8% EC 97.6% NTL <![CDATA[9.4×10 -7 ]]> WA 96.7% PIA 98.5% WLE 82min

[0123] Note: MCA model construction accuracy; DFQ data fusion quality; NTL network training loss; PIA parameter inversion accuracy; TSR twin synchronization rate; EC inference consistency; WA early warning accuracy; WLE early warning lead time.

[0124] All values ​​in Table 4 are from actual calculation output. The model construction accuracy is 98.7%, the data fusion quality is 97.8%, the network training loss reaches a high-precision convergence level, the parameter inversion accuracy is 98.5%, the twin synchronization rate is 99.3%, the inference consistency is 97.6%, the early warning accuracy is 96.7%, and the early warning lead time is 82 minutes. All indicators have reached the expected design goals, proving that the method of this application can operate stably and reliably under real working conditions.

[0125] In the multi-timestep forward dynamic extrapolation process, the internal damage of the soil and rock mass continuously develops with plastic strain and stress state. Based on the local damage variable evolution equation proposed in this application, a spatiotemporal evolution cloud map of the global damage variable can be obtained, such as... Figure 7 As shown, the data includes three key states: the initial stage, the middle stage of evolution, and the pre-disaster stage. Damage values ​​are displayed using a three-color grading system: blue for low damage, green for medium damage, and red for high damage. At the initial stage, damage is mainly concentrated at the slope toe and weak interlayers, with overall low values ​​consistent with the initial stress state of the geological body. During the middle stage, damage gradually expands along the potential sliding direction, with the area of ​​high-damage regions continuously increasing, consistent with the development path of the plastic zone. At the pre-disaster stage, the high-damage area forms a continuous, interconnected zone, completely coinciding with the potential sliding surface. The damage distribution, expansion direction, and development rate are highly consistent with geotechnical mechanics theory and the actual landslide evolution law, proving that the damage variable definition and evolution equation proposed in this application can accurately reflect the deterioration process of the internal structure of the geological body, providing a reliable basis for calculating the probability of interconnection.

[0126] The instability of a geological body is essentially a sudden transition from a stable state to a chaotic state. Based on the phase space reconstruction and Lyapunov exponent calculation method proposed in this application, the temporal variation of the maximum Lyapunov exponent can be obtained, such as... Figure 8 As shown, in the early stage of the simulation, the exponent value remained negative and the absolute value remained stable, indicating that the system was in a convergent and stable state and the geological body was generally safe. As the simulation progressed, the exponent gradually increased and the negative value decreased, indicating that the stability of the geological body continued to decrease. At step 98, the exponent turned from negative to positive and clearly crossed the zero threshold. At the same time, the second derivative of the curve was simultaneously greater than zero, marking that the system entered a chaotic divergent state and the geological body entered the stage of imminent instability. This judgment time was completely consistent with the preset simulation instability time, without any lag or misjudgment. This proves that the maximum Lyapunov exponent spectrum, phase space reconstruction, and perturbation orthogonalization calculation method adopted in this application can accurately and stably identify the critical point of geological body instability, providing a scientific, quantitative, and reliable judgment basis for imminent disaster early warning.

[0127] This embodiment details the full-process simulation verification under real geological conditions, achieving stable closed-loop operation of damage evolution, sliding surface breakthrough probability calculation, and Lyapunov exponent critical determination. It verifies the engineering applicability of the damage evolution equation, coupled deduction format, and chaos criterion, proving that the method has high accuracy, high real-time performance, and strong feasibility under real conditions.

[0128] The above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. For those skilled in the art, the present invention can have various modifications and variations. Any changes, modifications, substitutions, integrations, and parameter changes made to these embodiments within the spirit and principles of the present invention, without departing from the principles and spirit of the present invention, through conventional substitutions or to achieve the same function, fall within the scope of protection of the present invention.

Claims

1. A method for early warning of geological disaster risks based on multi-source heterogeneous data fusion, characterized in that, include: A three-dimensional geomechanical numerical model of the target geological hazard body is established. The numerical model defines the physical framework of hazard evolution based on the constitutive relationship of rock and soil. Construct a multi-source heterogeneous data stream, which includes space-based InSAR deformation data, airborne UAV lidar point cloud data, ground-based GNSS displacement data, and deep underground stress monitoring data. Using the Physical Information Neural Network (PINN) as a surrogate model, the residuals of the partial differential equations of the numerical model are added as constraints to the loss function. By minimizing the difference between the monitoring data and the model predictions, the unsteady mechanical parameters in the numerical model are inverted and updated in real time. The parameters include at least time-varying cohesion, internal friction angle and permeability coefficient, thus obtaining a digital twin that matches the real state of the current geological body in real time. Based on the updated digital twin, combined with short-term weather forecast data for future periods, forward dynamic extrapolation at multiple time steps is performed to calculate the probability of the cumulative plastic zone and potential sliding surface connecting at future moments. Based on the simulation results, the Lyapunov index spectrum of the state is calculated. When the maximum Lyapunov index turns from negative to positive, it is determined that the geological body has entered a chaotic and unstable state, and a disaster risk warning signal is generated accordingly.

2. The geological disaster risk early warning method based on multi-source heterogeneous data fusion according to claim 1, characterized in that, The construction of the multi-source heterogeneous data stream also includes: Spatiotemporal references were unified for space-based InSAR deformation data, airborne UAV lidar point cloud data, ground-based GNSS displacement data, and deep underground stress monitoring data. A Kalman filter fusion algorithm was then used to generate high-confidence observation sequences. Its state update equation is: ,in, For the true state vector, For the observation matrix, To observe the noise, Let be the noise covariance matrix.

3. The geological disaster risk early warning method based on multi-source heterogeneous data fusion according to claim 1, characterized in that, The method of using the Physical Information Neural Network (PINN) as a surrogate model and incorporating the residuals of the partial differential equations of the numerical model as constraints into the loss function includes: Construct the total loss function of the physical information neural network The total loss function Data mismatch Physical equation residuals and boundary condition residuals The weighted composition is calculated using the following formula: ,in, For neural networks at position and time Output displacement prediction value, For observations in a multi-source heterogeneous data stream, For the residual operator of the governing equation, For boundary condition operators, , , These are the weighting coefficients for each item. , , These represent the number of observation data points, the number of physical equation configuration points, and the number of boundary condition points, respectively.

4. The geological disaster risk early warning method based on multi-source heterogeneous data fusion according to claim 3, characterized in that, The physical equation residual term The construction includes: Based on Biot's consolidation theory, the residual between the mass conservation equation and the momentum conservation equation is defined as: ,in, Let ∇ be the effective stress tensor and ∇ be the divergence operator. For Biot coefficient, For the pore water pressure gradient, The vector of gravitational acceleration. displacement vector The second derivative with respect to time, Porosity For fluid density, For fluid velocity, For source and sink items.

5. The geological disaster risk early warning method based on multi-source heterogeneous data fusion according to claim 1, characterized in that, The real-time inversion and updating of the unsteady mechanical parameters in the numerical model includes: The time-varying cohesion internal friction angle and permeability coefficient The parameter is parameterized as an implicit function that evolves over time. The sensitivity of the stress tensor to the strain tensor is calculated using the automatic differentiation technique of the physical information neural network. The parameter gradient is solved using the adjoint state method, and the update formula is: ,in, The set of unsteady mechanical parameters to be inverted. The cohesion of the soil and rock mass, The internal friction angle of the rock and soil mass. The permeability coefficient of the rock and soil mass. In the first The unsteady mechanical parameter vector updated in the next iteration. In the first The unsteady mechanical parameter vector of the next iteration Total loss function Regarding parameters In the The gradient vector at the next iteration For learning rate, The momentum coefficient, The regularization coefficient is . To constrain the parameters within the physically feasible region Projection operator within.

6. The geological disaster risk early warning method based on multi-source heterogeneous data fusion according to claim 1, characterized in that, The aforementioned forward dynamics extrapolation with multiple time steps includes: The updated dynamic equilibrium equations of the digital twin are solved using explicit-implicit hybrid integration. These equations consider the coupling effect of pore water pressure and skeleton stress, and their discretized form is as follows: ,in, , , These are the mass matrix, damping matrix, and stiffness matrix, respectively. Let be the nodal displacement vector. For external load vectors, The pore water pressure field is derived from short-term weather forecast data. It is a matrix of shape functions.

7. The geological disaster risk early warning method based on multi-source heterogeneous data fusion according to claim 6, characterized in that, The explicit-implicit hybrid integral includes: The dynamic term is explicitly integrated using the central difference method, and the seepage term is implicitly integrated using the backward difference method. The coupled solution scheme is as follows: in, In order to be in The nodal displacement vector at time t. In order to be in The nodal displacement vector at time t. For time step, In order to be in The nodal velocity vector at time t, It is the inverse of the mass matrix. In order to be in The external nodal load vector at time t. In order to be in The internal node resistance vector at time t; In order to be in The nodal pore water pressure vector at time t. In order to be in The nodal pore water pressure vector at time t. It is the inverse of the penetration matrix. In order to be in The external fluid flow vector at time t. This is the fluid-structure interaction matrix.

8. The geological disaster risk early warning method based on multi-source heterogeneous data fusion according to claim 1, characterized in that, The calculation of the probability of the cumulative plastic zone and the potential sliding surface connecting at future time steps includes: Define local damage variables The local damage variable is based on cumulative plastic strain. and stress triaxiality The construction, its evolution equation is: ,in, and Let be the material damage constant. For at any time The equivalent plastic strain rate; Based on local damage variables Calculate the probability of full connectivity. The global penetration probability is determined by calculating the seepage threshold of the damage field gradient: ,in, As a potential connecting route, This is the critical damage threshold.

9. The geological disaster risk early warning method based on multi-source heterogeneous data fusion according to claim 1, characterized in that, The calculation of the Lyapunov exponent spectrum of the state based on the deduction results includes: Constructing the reconstructed phase space vector of the state space ,in It is in a displacement state. In terms of speed state, The state of time-varying mechanical parameters; Apply small perturbations to the reconstructed phase space trajectory Solve the variational equation ,in In the state Jacobian matrix at the location; The Gram-Schmidt orthogonalization method is used to periodically reorthogonalize the disturbance vectors of the discretized time series, and the maximum Lyapunov exponent is calculated. The formula is: Calculated by sliding time window The time series curve is used to determine the bifurcation of dynamic behavior from steady state to chaotic state when the curve crosses the preset zero value criterion threshold and the second derivative is greater than zero, thus generating a disaster risk warning signal.

10. The method for early warning of geological disaster risks based on multi-source heterogeneous data fusion according to claim 9, characterized in that, The determination that a geological body has entered a chaotic and unstable state generates a pre-disaster risk warning signal, including: Set tiered early warning thresholds ,when And duration At that time, a Level 1 emergency disaster warning was triggered. This is a preset critical value; Simultaneously calculate the fractal dimension. ,like and If the characteristics of non-integer dimensions are exhibited, it is determined that the system has entered a chaotic evolution stage dominated by strange attractors.

Images