Real-time monitoring and early warning system for foundation pit slope stability

Abstract

The application discloses a kind of real-time monitoring and early warning method for foundation pit slope stability, which belongs to geotechnical engineering monitoring technical field.For the problem of multi-source data fusion difficulty, coupling analysis deficiency and single early warning mode, displacement, deep deformation, groundwater level, soil moisture content, rainfall intensity and construction vibration data are collected, aligned, converted and fused by Kalman filter without trace through edge gateway, and the coupling strength matrix is extracted to select the dominant space-time sub-block.When the coupling strength increases above the threshold value, the mechanism-data dual-driven hybrid early warning model is triggered, the extended finite element seepage-stress coupling equation is embedded into the LSTM loss function as a physical constraint term, and the conditional Wasserstein generative adversarial network is used to generate virtual samples on the Riemannian manifold to expand the training set.The instability probability and early warning level are outputted, and the closed-loop response is automatically triggered in the digital twin platform dynamically mapped to BIM model.The application is used for all-weather real-time monitoring and grading early warning of foundation pit slope, and improves the prediction accuracy and response efficiency.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering monitoring technology. More specifically, this invention relates to a method for real-time monitoring and early warning of the stability of foundation pit slopes. Background Technology

[0002] In monitoring the stability of foundation pit slopes, existing technologies typically deploy sensors for surface displacement, deep deformation, groundwater level, soil moisture content, rainfall intensity, and construction vibration, with each sensor collecting data independently. In actual engineering projects, differences in sampling frequency, start-up time, and internal clock drift among different sensors make it difficult to accurately align data message timestamps. Furthermore, each sensor uses its own local coordinate system, making manual registration cumbersome and prone to introducing errors. This results in a lack of accurate temporal and spatial correlation between displacement and data such as water level and moisture content.

[0003] In terms of data fusion, traditional methods often rely on simple statistical averaging or threshold judgment, which are difficult to effectively handle the nonlinear coupling relationships between multiple physical quantities. In particular, construction vibration signals are transient and have high amplitude characteristics, which can easily mask the changing trends of slowly varying physical quantities such as displacement. Conventional filtering methods are unable to accurately remove vibration interference without losing effective information.

[0004] For early warning mechanisms, existing technologies generally employ single or multiple independent thresholds. Displacement changes may be caused by rising groundwater levels or by instantaneous impacts such as construction vibrations, and the corresponding instability risks differ. Existing methods struggle to distinguish the causes, easily leading to false alarms or missed alarms. More critically, the coupling strength between displacement and other physical quantities changes dynamically, and fixed-threshold early warning methods cannot track this evolutionary trend, making it difficult to issue effective warnings in the early stages of instability.

[0005] In building early warning models, data-driven machine learning methods face the problem of insufficient samples: slope instability is a low-probability event, complete cases are limited, and the instability process includes multiple evolution rate stages, with insufficient training samples for transitional forms. Furthermore, purely data-driven models lack constraints from physical and mechanical mechanisms, and prediction results may violate the seepage-stress coupling law, reducing reliability. Attempts to introduce physical equations face difficulties in balancing computational complexity and data loss.

[0006] The existing system operates in a relatively rigid mode, analyzing at a fixed frequency during the daily stable phase, resulting in wasted resources. When entering a potentially unstable phase, it cannot automatically switch to a high-precision mode, leading to a delayed response. Early warning results are typically expressed numerically or in simple charts, lacking effective integration with the engineering information model, resulting in low closed-loop efficiency from early warning issuance to on-site verification, handling, and feedback. Summary of the Invention

[0007] One object of the present invention is to solve at least the above-mentioned problems and to provide at least the advantages that will be described later.

[0008] Another objective of this invention is to provide a method for real-time monitoring and early warning of foundation pit slope stability, which aims to solve the technical problems in the prior art, such as the spatiotemporal asynchrony and difficulty in fusion of multi-source heterogeneous sensor data, the difficulty in effectively eliminating transient construction vibration interference, the inability of fixed threshold early warning to track the dynamic evolution of the coupling strength of multiple physical quantities, the lack of physical mechanism constraints in data-driven models under small sample conditions, and the inability of early warning modes to adaptively switch and the low efficiency of closed-loop response.

[0009] To achieve these objectives and other advantages of the present invention, a method for real-time monitoring and early warning of foundation pit slope stability is provided, comprising the following steps: Step 1: Collect surface displacement, deep deformation, groundwater level, soil moisture content, rainfall intensity, and construction vibration data of the foundation pit slope; each sensor adds type, accuracy level, and GNSS timestamp to the data message header and connects to the edge computing gateway through a unified interface; Step 2: The edge computing gateway aligns the timestamps with GNSS integer second pulses, transforms all coordinates to the BIM coordinate system, and fuses multiple physical quantities using unscented Kalman filtering. The data is constructed into a third-order tensor according to physical quantities, time windows, and spatial measurement points. The core tensor and factor matrix are obtained through Tucker decomposition. The coupling strength matrix is ​​extracted from the core tensor. Spatiotemporal sub-blocks dominated by displacement and simultaneously satisfying the coupling strength conditions with groundwater level and soil moisture content are selected. The corresponding original signals are extracted, and transient impact sub-blocks dominated by construction vibration are removed. A standardized dataset is formed and uploaded to the cloud. Step 3: Track the coupling strength matrix in the cloud. When the coupling strength between displacement and other physical quantities increases and exceeds the threshold in multiple consecutive time windows, the mechanism-data dual-driven hybrid early warning model is triggered to switch from low-power inference mode to high-precision early warning mode. The extended finite element seepage-stress coupling equation is used as a physical constraint term, and LSTM loss function is embedded for joint training. A physical residual loss term is added in addition to the mean square error loss. At the same time, the measured instability cases are divided into multiple types according to the evolution speed. Between adjacent instability stages, a conditional Wasserstein generative adversarial network is used to learn the Riemannian manifold and perform geodesic interpolation to generate virtual samples to expand the training set. The model outputs the slope instability probability and early warning level. Step 4: Dynamically map the warning level and instability probability to the BIM model on the digital twin platform, and automatically trigger a closed-loop response process of verification, handling and feedback based on the warning level.

[0010] Preferably, step two, which involves obtaining the core tensor and factor matrix using Tucker decomposition, includes: a. Expand the third-order tensor into a two-dimensional matrix along the monitoring physical quantity mode, perform singular value decomposition on the two-dimensional matrix, and take the first few columns of the left singular vector matrix as the monitoring physical quantity factor matrix. b. Expand the same third-order tensor into a two-dimensional matrix along the time window dimension, perform singular value decomposition on the two-dimensional matrix, and take the first few columns of the left singular vector matrix as the time window factor matrix. c. Expand the same third-order tensor into a two-dimensional matrix along the spatial measurement point mode, perform singular value decomposition on the two-dimensional matrix, and take the first few columns of the left singular vector matrix as the spatial measurement point factor matrix. d. Perform a product operation on the first dimension of the original third-order tensor and the transpose of the monitoring physical quantity factor matrix. The operation method is as follows: for each element in the tensor, multiply the index on the monitoring physical quantity dimension with the corresponding value in the transpose of the factor matrix and sum them to obtain a new third-order tensor. e. Perform a product operation in the second dimension on the third-order tensor obtained in step four and the transpose of the time window factor matrix. Multiply the index in the time window dimension with the corresponding value in the transpose of the factor matrix and sum them to obtain another third-order tensor. f. Perform a product operation on the third dimension of the third-order tensor obtained in step 5 and the transpose of the spatial measurement point factor matrix. Multiply the indices on the spatial measurement point dimension with the corresponding values ​​in the transpose of the factor matrix and sum them to obtain the core tensor. g. Extract the absolute value of each element in the core tensor, and arrange these absolute values ​​according to the index positions of the corresponding physical quantities, time windows, and spatial measurement points in the original third-order tensor to form a three-dimensional array with the same structure as the original third-order tensor. This array is the coupling strength matrix.

[0011] Preferably, step two, which involves selecting spatiotemporal sub-blocks whose dominant mode is displacement and which simultaneously satisfy the coupling strength conditions with groundwater level and soil moisture content, includes: In the monitoring physical quantity factor matrix A, the row vectors corresponding to displacement, groundwater level, and soil moisture content are located. The cosine similarity between the displacement row vector and the groundwater level row vector, and the cosine similarity between the displacement row vector and the soil moisture content row vector are calculated. In the core tensor G, the slice matrix corresponding to the displacement is extracted along the dimension of the monitoring physical quantity. The elements in the slice matrix are thresholded, and the index positions with the absolute value of the element greater than the preset threshold are marked as strong coupling candidate positions. Among the strong coupling candidate positions, the index positions that simultaneously satisfy the displacement-groundwater level cosine similarity exceeding the first threshold and the displacement-soil moisture content cosine similarity exceeding the second threshold are selected to form the target spatiotemporal sub-block. The original signal components corresponding to the third-order tensor at the target spatiotemporal sub-block positions are retained, and the signal components corresponding to the other positions are set to zero to obtain the signal after removing the transient impact sub-block.

[0012] Preferably, the calculation steps for the physical residual loss term in step three include: On the extended finite element discrete nodes, the displacement field u of the i-th node predicted by the Long Short-Term Memory network is extracted. i and pore water pressure field p i Substituting into the discrete form of the seepage-stress coupling control equation: K u u i +K c p i -f i =R u，i ; K c T (du) i / dt)+K p p i -q i =R p，i ; Among them, K u K is the stiffness matrix. c Let K be the coupling matrix. p Let f be the penetration matrix. i Let q be the nodal force vector of the i-th node. i Let R be the node flow vector of the i-th node. u，i Let R be the displacement field residual of the i-th node. p，i Let be the residual of the pore water pressure field at the i-th node; The physical residual loss term is: ; Where N is the total number of discrete nodes, and λ is the weighting coefficient for the difference in dimensions between the residual of the equilibrium displacement field and the residual of the pore water pressure field. The overall loss function for joint training is: ; Where L MSE α is the mean squared error loss term, and α is the physical residual penalty coefficient.

[0013] Preferably, step three, which involves generating a Riemannian manifold of the latent space using conditional Wasserstein and performing geodesic interpolation on the Riemannian manifold, includes: The measured samples of various instability evolution paths were divided into initial creep stage, constant creep stage, accelerated creep stage, and critical instability stage according to time series. A conditional Wasserstein generative adversarial network was constructed. The generator input was random noise z and stage condition label c, and the output was the virtual sample of the corresponding stage. The discriminator output was the Wasserstein distance estimate. After the conditional Wasserstein generative adversarial network training converged, a set of random noise z1 and z2 was sampled for adjacent instability stages, and mapped to the latent space by the generator to obtain latent vectors w1 and w2. The geodesic path between the latent vectors was calculated. ; Where exp and log are the exponential and logarithmic mappings on the Riemannian manifold, respectively; Interpolation progress parameters from the starting point to the ending point ( When =0, it corresponds to w1. When =1 corresponds to w2); in Uniform sampling within the value range M Each interpolation point is used to decode the latent vector corresponding to each interpolation point through a generator, resulting in a virtual instability evolution sample that smoothly transitions between adjacent instability stages.

[0014] Preferably, the switching steps between low-power inference mode and high-precision warning mode in step three include: In low-power inference mode, the mechanism-data dual-driven hybrid early warning model loads a lightweight branch of the Long Short-Term Memory (LSTM) network. This lightweight branch is obtained from the original LSTM network after channel pruning and weight quantization, and the physical residual loss term is not calculated during inference. When the cloud server detects that the displacement-groundwater level coupling strength and the displacement-soil moisture content coupling strength are monotonically increasing within k consecutive time windows, and the coupling strength of the current time window exceeds their respective preset thresholds, the mode switch is triggered. The switching process includes: loading the complete LSTM network weights from the model storage area to replace the lightweight branch weights, and activating the computation node corresponding to the physical residual loss term in the computation graph, so that the model enters the high-precision early warning mode.

[0015] Preferably, the deep deformation monitoring in step one is achieved by deploying an array-type microelectromechanical inclinometer, the packaging structure of which includes: MEMS accelerometer chips are soldered onto a flexible printed circuit board, which is fixed to the inner wall of a glass fiber reinforced epoxy resin composite tube by rigid partitions arranged at equal intervals. The outer wall of the composite tube is formed by spirally winding epoxy resin-impregnated glass fiber cloth along the tube axis at an angle of 45° to 55°. Both ends of the composite tube are sealed with stainless steel end caps and O-rings. The temperature self-compensation circuit includes a thermistor from the same batch as the MEMS accelerometer chip. The thermistor and the MEMS accelerometer chip are arranged adjacent to each other on the same flexible printed circuit board. The temperature self-compensation circuit corrects the accelerometer output in real time through a pre-calibrated temperature-drift curve. The redundant MEMS sensing unit consists of two sets of MEMS accelerometer chips arranged orthogonally on the same rigid partition.

[0016] Preferably, the implementation steps of the unified data interface in step one include: A predefined data dictionary covering displacement, pressure, moisture content, and vibration parameters is used. Each entry in the data dictionary includes a data code, physical quantity name, unit, data type, and value range. The edge computing gateway allocates independent communication ports and message queues for each type of sensor. Each sensor encapsulates data packets according to the data dictionary format. The header of the data packet contains the sensor type code, accuracy level code, and GNSS timing timestamp, while the body of the packet contains monitoring data in key-value pair format. The data parsing module of the edge computing gateway reads the sensor type code from the packet header, queries the data dictionary to obtain the key-value pair parsing template corresponding to that type of sensor, extracts the monitoring data according to the parsing template, and writes it into the time-series database.

[0017] Preferably, step two, which involves using unscented Kalman filtering to perform optimal fusion estimation of multiple physical quantities, includes: Constructing the state vector: ; Where d x d y d z For three-dimensional displacement components, v x v y v z The displacement rate components are three-dimensional; the state transition equation adopts a constant velocity model; in the observation equation, the observation vector of the GNSS receiver is [d]. x d y d z ] T The observation vector of the ground-based synthetic aperture radar is transformed into the BIM coordinate system by the line-of-sight projection matrix as d. x d y d zA linear combination; by querying the preset accuracy level-variance mapping table using the accuracy level labels attached to each sensor message, the observation noise covariance matrix R of each sensor is obtained. k The diagonal element values; calculating the Sigma point set and its weights during the time update phase, and using R during the measurement update phase. k To measure the noise covariance, calculate the Kalman gain and perform state updates.

[0018] Preferably, in step one, surface displacement monitoring simultaneously deploys a BeiDou / GNSS receiver and a ground-based synthetic aperture radar (SAR). In areas where satellite signals are obscured, the SAR interferometry is used to obtain a line-of-sight deformation sequence. Pixels with amplitude deviation indices less than a preset threshold in the time-series interferogram are selected as candidate permanent scatterers. After time-dimensional phase unwrapping of the candidate points, a differential phase model between adjacent points is established. ; Where Δv is the deformation rate difference, T is the time baseline, B⊥ is the vertical baseline, R is the slant distance, θ is the incident angle, Δh is the elevation error difference, and Δϕ is the vertical distance. atm For atmospheric delay phase difference, Δϕ noise This is the noise phase difference; The atmospheric delay phase component of each permanent scatterer candidate point is separated and estimated by spatiotemporal filtering. The atmospheric delay phase component is then subtracted from the total phase to obtain the corrected deformation phase sequence.

[0019] The present invention has at least the following beneficial effects: 1. The real-time monitoring and early warning method for foundation pit slope stability provided by this invention solves the problem of time and space synchronization difficulties for multi-source heterogeneous sensor data through a unified data interface and edge computing gateway. Employing unscented Kalman filtering and tensor decomposition techniques, it achieves deep fusion and coupling strength analysis of multiple physical quantities such as displacement, water level, and rainfall, intelligently filtering out spatiotemporal sub-blocks closely related to instability and effectively eliminating transient interference signals such as construction vibrations. Based on this, a mechanism-data dual-driven hybrid early warning model is constructed, embedding physical equations into an LSTM neural network, overcoming the lack of physical constraints in purely data-driven models and improving prediction reliability. Simultaneously, generative adversarial networks are used to generate virtual samples on Riemannian manifolds, solving the small-sample learning problem. Finally, visualization and closed-loop response are achieved on a digital twin platform, significantly improving the early warning accuracy and emergency response efficiency of all-weather monitoring.

[0020] 2. The Tucker decomposition method employed in this invention provides a systematic and repeatable mathematical path for extracting core features from complex third-order monitoring tensors (physical quantities, time, and space). By sequentially performing singular value decomposition along each mode and obtaining the factor matrix, and then obtaining the core tensor through dimensional product operations, this process can compress the original high-dimensional redundant monitoring data into a low-dimensional "coupling strength matrix." This matrix intuitively quantifies the strength of the interaction between different physical quantities at different times and spatial locations, thus providing a quantitative and structured basis for subsequent identification of dominant instability factors and strongly coupled regions. This method solves the problem of effectively handling nonlinear coupling relationships of multiple physical quantities in existing technologies, has high computational efficiency, and retains the main information of the original tensor, making it a key technical link in realizing multi-field coupling analysis.

[0021] 3. This invention further refines the method of screening key spatiotemporal sub-blocks based on the coupling strength matrix, achieving intelligent and precise removal of transient high-amplitude interference such as construction vibration. Traditional filtering methods struggle to distinguish between displacement abrupt changes caused by vibration and the true accelerated deformation that is a precursor to instability. This method, however, calculates the cosine similarity between displacement and water level / water content vectors in the factor matrix, and combines this with threshold segmentation in the core tensor to accurately locate the spatiotemporal positions where "displacement-water level-water content" are strongly coupled. These positions are precisely the potential instability regions truly related to seepage-stress coupling; conversely, signals that do not meet the strong coupling conditions are eliminated by setting them to zero. This significantly improves the signal-to-noise ratio of the source data, effectively avoiding false alarms caused by interference signals, and providing a clean and reliable data foundation for subsequent high-precision early warning.

[0022] 4. This invention addresses the problem of scarce slope instability samples, particularly the lack of transitional samples between different evolution stages, by proposing a virtual sample generation method based on Riemannian manifold geodesic interpolation. This method first divides a limited number of measured instability cases into creep stages, then uses Conditional Wasserstein to generate an adversarial network to learn the Riemannian manifold structure in the latent space. By calculating the geodesic paths between latent vectors of adjacent instability stages on the manifold and performing uniform interpolation, physically smooth, continuous, and in line with the instability evolution law intermediate transitional samples can be generated. This differs from simple linear interpolation; the geodesic ensures that the generated samples lie on the shortest path of the data manifold, resulting in higher realism and diversity. This technique effectively expands the training set, especially enriching the key transitional data from slow deformation to accelerated instability, significantly alleviating the small sample learning problem and improving the early warning model's ability to capture the entire instability evolution process.

[0023] 5. The adaptive switching mechanism between low-power inference mode and high-precision early warning mode proposed in this invention achieves a dynamic balance between system energy efficiency and early warning performance. When the foundation pit slope is in a stable or weakly coupled state, the system runs in low-power mode, loading only a lightweight LSTM network that has undergone pruning and quantization, and does not calculate complex physical residual terms, significantly reducing the computational load and energy consumption in the cloud. When the cloud continuously tracks and finds that the displacement-water level and displacement-water content coupling strengths show a monotonically increasing trend and exceed the threshold in multiple time windows, the system automatically judges that the slope is entering a potential instability trend, thereby triggering mode switching, loading the complete high-precision model, and activating physical constraints. This mechanism ensures the economy of daily monitoring and can provide the highest accuracy early warning analysis in a timely manner when the risk increases, avoiding resource waste caused by long-term high-load operation or safety risks caused by delayed early warning.

[0024] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Detailed Implementation

[0025] The present invention will be further described in detail below with reference to embodiments, so that those skilled in the art can implement it based on the description.

[0026] It should be understood that terms such as “having,” “comprising,” and “including” as used herein do not exclude the presence or addition of one or more other elements or combinations thereof.

[0027] It should be noted that, unless otherwise specified, the experimental methods described in the following implementation plan are all conventional methods, and the reagents and materials described are all commercially available unless otherwise specified.

[0028] This invention provides a method for real-time monitoring and early warning of the stability of foundation pit slopes, which includes the following steps: Step 1: Collect surface displacement, deep deformation, groundwater level, soil moisture content, rainfall intensity, and construction vibration data of the foundation pit slope; each sensor adds type, accuracy level, and GNSS timestamp to the data message header and connects to the edge computing gateway through a unified interface; Step 2: The edge computing gateway aligns the timestamps with GNSS integer second pulses, transforms all coordinates to the BIM coordinate system, and fuses multiple physical quantities using unscented Kalman filtering. The data is constructed into a third-order tensor according to physical quantities, time windows, and spatial measurement points. The core tensor and factor matrix are obtained through Tucker decomposition. The coupling strength matrix is ​​extracted from the core tensor. Spatiotemporal sub-blocks dominated by displacement and simultaneously satisfying the coupling strength conditions with groundwater level and soil moisture content are selected. The corresponding original signals are extracted, and transient impact sub-blocks dominated by construction vibration are removed. A standardized dataset is formed and uploaded to the cloud. Step 3: Track the coupling strength matrix in the cloud. When the coupling strength between displacement and other physical quantities increases and exceeds the threshold in multiple consecutive time windows, the mechanism-data dual-driven hybrid early warning model is triggered to switch from low-power inference mode to high-precision early warning mode. The extended finite element seepage-stress coupling equation is used as a physical constraint term, and LSTM loss function is embedded for joint training. A physical residual loss term is added in addition to the mean square error loss. At the same time, the measured instability cases are divided into multiple types according to the evolution speed. Between adjacent instability stages, a conditional Wasserstein generative adversarial network is used to learn the Riemannian manifold and perform geodesic interpolation to generate virtual samples to expand the training set. The model outputs the slope instability probability and early warning level. Step 4: Dynamically map the warning level and instability probability to the BIM model on the digital twin platform, and automatically trigger a closed-loop response process of verification, handling and feedback based on the warning level.

[0029] This technical solution specifically includes: Step 1: Various sensors are deployed on the slope of the foundation pit to collect monitoring data. Surface displacement monitoring simultaneously utilizes a BeiDou GNSS receiver and a ground-based synthetic aperture radar (SAR). In areas where satellite signals are blocked, line-of-sight deformation sequences are obtained using ground-based SAR interferometry. Pixels with amplitude deviation indices less than a preset threshold in the time-series interferogram are selected as candidate permanent scatterers. After time-dimensional phase unwrapping of the candidate points, a differential phase model between adjacent points is established. This model includes parameters such as deformation rate difference, time baseline, vertical baseline, slant range, incident angle, elevation error difference, atmospheric delay phase difference, and noise phase difference. The atmospheric delay phase component of each candidate permanent scatterer is separated and estimated using spatiotemporal filtering. This component is subtracted from the total phase to obtain the corrected deformation phase sequence. Deep deformation monitoring is achieved by deploying an array-type MEMS inclinometer. The inclinometer's packaging structure is as follows: a MEMS accelerometer chip is soldered onto a flexible printed circuit board, which is fixed to the inner wall of a glass fiber reinforced epoxy resin composite tube by rigid partitions arranged at equal intervals. The outer wall of the composite tube is formed by spirally winding epoxy resin-impregnated fiberglass cloth along the tube's axial direction, with the winding angle controlled between 45 and 55 degrees. Both ends of the composite tube are sealed with stainless steel end caps and O-rings. The temperature self-compensation circuit includes a thermistor from the same batch as the MEMS accelerometer chip. This thermistor and the MEMS accelerometer chip are located adjacent to each other on the same flexible printed circuit board. The temperature self-compensation circuit corrects the accelerometer output in real time using a pre-calibrated temperature-drift curve. The redundant MEMS sensing unit consists of two sets of MEMS accelerometer chips arranged orthogonally on the same rigid partition. Each sensor appends its type, accuracy class, and GNSS timestamp to the data packet header and connects to the edge computing gateway through a unified interface. The unified data interface is implemented by pre-defining a data dictionary covering displacement, pressure, moisture content, and vibration parameters. Each entry in the data dictionary includes a data code, physical quantity name, unit, data type, and value range. The edge computing gateway allocates independent communication ports and message queues for each type of sensor. Each sensor encapsulates data packets according to the data dictionary format. The header of the data packet contains, in sequence, the sensor type code, the accuracy level code, and the GNSS timing timestamp. The body of the packet contains monitoring data in key-value pair format. The data parsing module of the edge computing gateway reads the sensor type code from the packet header, queries the data dictionary to obtain the key-value pair parsing template corresponding to that sensor type, extracts the monitoring data according to the parsing template, and writes it to the time-series database.

[0030] Step 2: The edge computing gateway aligns the timestamps of each sensor with GNSS integer-second pulses and transforms all coordinates to the BIM coordinate system. Unscented Kalman filtering is used for optimal fusion estimation of multiple physical quantities: a state vector containing three-dimensional displacement components and three-dimensional displacement rate components is constructed, and the state transition equation adopts a constant velocity model. In the observation equation, the observation vector of the GNSS receiver is the three-dimensional displacement component, and the observation vector of the ground-based synthetic aperture radar is transformed to a linear combination of the three-dimensional displacement components in the BIM coordinate system via a line-of-sight projection matrix. The pre-defined accuracy level-variance mapping table is queried using the accuracy level label attached to each sensor message to obtain the diagonal element values ​​of the observation noise covariance matrix of each sensor. During the time update phase, the Sigma point set and its weights are calculated. During the measurement update phase, the Kalman gain is calculated using this covariance matrix as the measurement noise covariance, and the state is updated. After filtering and fusion, the data is constructed into a third-order tensor according to three dimensions: physical quantity type, time window sequence number, and spatial measurement point number. Tucker decomposition is performed on this third-order tensor to obtain the core tensor and factor matrix. The decomposition process is as follows: First, expand the third-order tensor into a two-dimensional matrix along the monitoring physical quantity mode. Perform singular value decomposition (SVD) on this two-dimensional matrix, and take the first few columns of the left singular vector matrix as the monitoring physical quantity factor matrix. Second, expand the same third-order tensor into a two-dimensional matrix along the time window dimension. Perform SVD on this two-dimensional matrix, and take the first few columns of the left singular vector matrix as the time window factor matrix. Third, expand the same third-order tensor into a two-dimensional matrix along the spatial measurement point mode. Perform SVD on this two-dimensional matrix, and take the first few columns of the left singular vector matrix as the spatial measurement point factor matrix. Fourth, multiply the original third-order tensor with the transpose of the monitoring physical quantity factor matrix in the first dimension. This is done by multiplying each element of the tensor by its index in the monitoring physical quantity dimension with the corresponding value in the transpose of the factor matrix and summing the results, resulting in a new third-order tensor. Fifth, multiply this new tensor with the transpose of the time window factor matrix in the second dimension. This is done by multiplying the index in the time window dimension with the corresponding value in the transpose of the factor matrix and summing the results, resulting in yet another third-order tensor. Next, multiply the tensor with the transpose of the spatial measurement point factor matrix in the third dimension. Multiply the indices of the spatial measurement points with the corresponding values ​​in the transpose of the factor matrix and sum them to obtain the core tensor. Extract the absolute value of each element in the core tensor and arrange these absolute values ​​according to the corresponding physical quantity, time window, and spatial measurement point index position in the original third-order tensor to form a three-dimensional array with the same structure as the original third-order tensor. This array is the coupling strength matrix. From the coupling strength matrix, select spatiotemporal sub-blocks that are displacement-dominant and simultaneously satisfy the coupling strength conditions with groundwater level and soil moisture content. The specific selection method is as follows: locate the row vectors corresponding to displacement, groundwater level, and soil moisture content in the monitoring physical quantity factor matrix, and calculate the cosine similarity between the displacement row vector and the groundwater level row vector, as well as the cosine similarity between the displacement row vector and the soil moisture content row vector.In the core tensor, a slice matrix corresponding to the displacement is extracted along the dimension of the monitored physical quantity. Elements in this slice matrix are thresholded, and index positions with absolute values ​​greater than a preset threshold are marked as strong coupling candidate positions. Among these strong coupling candidate positions, index positions that simultaneously satisfy a displacement-groundwater level cosine similarity exceeding a first threshold and a displacement-soil moisture content cosine similarity exceeding a second threshold are selected to form target spatiotemporal sub-blocks. The original signal components corresponding to the third-order tensor positions in the target spatiotemporal sub-blocks are retained, while the signal components corresponding to other positions are set to zero, thereby eliminating transient impact sub-blocks dominated by construction vibration. A standardized dataset is then uploaded to the cloud.

[0031] Step 3: The cloud server continuously tracks the received coupling strength matrix. The system operates normally in low-power inference mode. In this mode, the mechanism-data dual-driven hybrid early warning model loads a lightweight branch of the LSTM. This lightweight branch is obtained from the original LSTM network after channel pruning and weight quantization. During inference, the physical residual loss term is not calculated to reduce computational resource consumption. When the cloud server detects that the displacement-groundwater level coupling strength and the displacement-soil moisture content coupling strength both exhibit monotonically increasing trends over multiple consecutive time windows, and the coupling strength in the current time window exceeds their respective preset thresholds, a mode switch is triggered. The switch process includes loading the complete LSTM network weights from the model storage area to replace the lightweight branch weights, and activating the computation nodes corresponding to the physical residual loss term in the computation graph, causing the model to enter a high-precision early warning mode. In high-precision early warning mode, the model uses the extended finite element seepage-stress coupling equation as a physical constraint term embedded in the LSTM loss function for joint training. The physical residual loss term is calculated as follows: On the discrete nodes of the extended finite element method, the displacement field and pore water pressure field of the i-th node predicted by the LSTM are extracted and substituted into the discrete form of the seepage-stress coupling control equation to obtain the displacement field residual and the pore water pressure field residual, respectively. The two control equations are: the displacement field residual is obtained by multiplying the stiffness matrix by the displacement field, adding the coupling matrix by the pore water pressure field, and subtracting the node force vector; the pore water pressure field residual is obtained by multiplying the transpose of the coupling matrix by the time derivative of the displacement field, adding the seepage matrix by the pore water pressure field, and subtracting the node flow rate vector. The physical residual loss term is equal to the sum of the squares of the displacement field residuals and the squares of the pore water pressure field residuals on all discrete nodes, multiplied by the dimensional balance weight coefficient, divided by the total number of nodes. The overall loss function for joint training is the mean square error loss term plus the physical residual penalty coefficient multiplied by the physical residual loss term. To address the challenge of small-sample learning, the measured instability cases were categorized into multiple types based on their evolution rate. Furthermore, the measured samples of each type's instability evolution path were divided into initial creep stage, constant-rate creep stage, accelerated creep stage, and critical instability stage based on time series. A Conditional Wasserstein Generative Adversarial Network (GAN) was constructed. The generator's input consisted of random noise and stage conditional labels, and its output was virtual samples corresponding to the stage. The discriminator's output was the Wasserstein distance estimate. After network training converged, a set of random noise was sampled for each adjacent instability stage and mapped to the latent space by the generator to obtain latent vectors. The geodesic path between latent vectors was calculated. This path was defined on the Riemannian manifold by exponential and logarithmic mappings, where the interpolation progress parameter corresponds to the starting latent vector at time zero and the ending latent vector at time one. Multiple interpolation points were uniformly sampled within the parameter's range. The latent vectors corresponding to each interpolation point were decoded by the generator to obtain virtual instability evolution samples showing a smooth transition between adjacent instability stages, used to expand the training set. The model ultimately outputs the slope instability probability and warning level.

[0032] Step four: Dynamically map the warning level and instability probability to the BIM model on the digital twin platform for visualization. An automatic closed-loop response process is triggered based on the warning level: when the warning level reaches the verification level, the system sends a verification instruction to on-site management personnel; when the warning level reaches the handling level, the system automatically generates a handling plan suggestion and pushes it to the relevant responsible persons; after on-site verification or handling is completed, a feedback report is submitted in the system, the system records the feedback information and updates the warning status, forming a complete closed loop of verification, handling, and feedback.

[0033] This technical solution extracts the coupling strength matrix through tensor decomposition and filters spatiotemporal sub-blocks related to seepage-stress coupling, effectively eliminating transient interference; it triggers a high-precision early warning mode by tracking the continuous increasing trend of coupling strength, achieving adaptive switching of the early warning mode; it embeds physical equations into the neural network loss function, making the prediction results conform to the laws of mechanics; and it uses geodesic interpolation on the Riemannian manifold to generate virtual samples, solving the small sample learning problem, thus outperforming existing methods in both prediction accuracy and response efficiency.

[0034] In another technical solution, step two, which involves obtaining the core tensor and factor matrix using Tucker decomposition, includes: a. Expand the third-order tensor into a two-dimensional matrix along the monitoring physical quantity mode, perform singular value decomposition on the two-dimensional matrix, and take the first few columns of the left singular vector matrix as the monitoring physical quantity factor matrix. b. Expand the same third-order tensor into a two-dimensional matrix along the time window dimension, perform singular value decomposition on the two-dimensional matrix, and take the first few columns of the left singular vector matrix as the time window factor matrix. c. Expand the same third-order tensor into a two-dimensional matrix along the spatial measurement point mode, perform singular value decomposition on the two-dimensional matrix, and take the first few columns of the left singular vector matrix as the spatial measurement point factor matrix. d. Perform a product operation on the first dimension of the original third-order tensor and the transpose of the monitoring physical quantity factor matrix. The operation method is as follows: for each element in the tensor, multiply the index on the monitoring physical quantity dimension with the corresponding value in the transpose of the factor matrix and sum them to obtain a new third-order tensor. e. Perform a product operation in the second dimension on the third-order tensor obtained in step four and the transpose of the time window factor matrix. Multiply the index in the time window dimension with the corresponding value in the transpose of the factor matrix and sum them to obtain another third-order tensor. f. Perform a product operation on the third dimension of the third-order tensor obtained in step 5 and the transpose of the spatial measurement point factor matrix. Multiply the indices on the spatial measurement point dimension with the corresponding values ​​in the transpose of the factor matrix and sum them to obtain the core tensor. g. Extract the absolute value of each element in the core tensor, and arrange these absolute values ​​according to the index positions of the corresponding physical quantities, time windows, and spatial measurement points in the original third-order tensor to form a three-dimensional array with the same structure as the original third-order tensor. This array is the coupling strength matrix.

[0035] In step a of this technical solution, the constructed third-order tensor is expanded into a two-dimensional matrix along the dimension of the monitored physical quantity. The third-order tensor originally has three dimensions: the type of monitored physical quantity, the time window number, and the spatial measurement point number. Expanding along the monitored physical quantity dimension means arranging all data in this dimension into rows of the matrix, and arranging all combinations of the other two dimensions—time window and spatial measurement point—into columns. The resulting two-dimensional matrix has the number of rows equal to the number of monitored physical quantity types, and the number of columns equal to the product of the number of time windows and the number of spatial measurement points. Then, singular value decomposition is performed on this two-dimensional matrix to obtain a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix. The first few columns of the left singular vector matrix are taken as the monitored physical quantity factor matrix. The number of columns selected here is determined by the cumulative singular value contribution rate; typically, the first columns, where the sum of singular values ​​accounts for more than 95% of the total sum of all singular values, are retained. This reduces the data dimensionality while retaining the main information.

[0036] In step b, the same third-order tensor is expanded into a two-dimensional matrix along the time window dimension. At this point, the number of rows in the matrix equals the number of time windows, and the number of columns equals the product of the number of monitored physical quantities and the number of spatial measurement points. After performing singular value decomposition on this matrix, the first few columns of the left singular vector matrix are taken as the time window factor matrix, with the selection rule also based on the singular value contribution rate.

[0037] In step c, the same third-order tensor is expanded into a two-dimensional matrix along the spatial measurement point dimension. At this point, the number of rows in the matrix equals the number of spatial measurement points, and the number of columns equals the product of the number of monitored physical quantities and the number of time windows. After performing singular value decomposition on this matrix, the first few columns of the left singular vector matrix are taken as the spatial measurement point factor matrix.

[0038] In step d, the original third-order tensor is multiplied in its first dimension by the transpose of the monitoring physical quantity factor matrix. This operation is called modular multiplication in tensor analysis. Specifically, it iterates through all elements of the original tensor. For each index position in the monitoring physical quantity dimension, the index is multiplied by the corresponding value in the transpose of the factor matrix. Then, the product results of all monitoring physical quantity indices are summed. After the operation, the first dimension of the original tensor is replaced by the number of columns in the monitoring physical quantity factor matrix, thus obtaining a new third-order tensor. Its dimension in the monitoring physical quantity direction is reduced from the original number of categories to the number of columns in the factor matrix, while the dimensions of the time window and spatial measurement points remain unchanged.

[0039] In step e, the third-order tensor obtained in step four is multiplied along its second dimension by the transpose of the time window factor matrix. The operation is similar to step d, but this time it is performed along the time window dimension. For each index in the time window dimension of the new tensor, the index is multiplied by the corresponding value in the transpose of the time window factor matrix, and the results are summed. After the operation, the second dimension of the new tensor is reduced to the number of columns in the time window factor matrix, resulting in another third-order tensor.

[0040] In step f, the third-order tensor obtained in step five is multiplied along its third dimension by the transpose of the spatial measurement factor matrix. Along the spatial measurement dimension, each index is multiplied by the corresponding value in the transpose of the spatial measurement factor matrix and the results are summed. After the operation, the third dimension is reduced to the number of columns in the spatial measurement factor matrix, and the resulting tensor is the core tensor. The core tensor has a much smaller dimension than the original third-order tensor, and each element represents the interaction strength under the corresponding factor combination.

[0041] In step g, the absolute value of each element in the core tensor is extracted. Since the core tensor has a different dimension than the original third-order tensor, they cannot be directly mapped by position. Therefore, it is necessary to reconstruct the core tensor and each factor matrix back to the original dimensional space through tensor multiplication. This involves performing a modular multiplication of the core tensor with the monitoring physical quantity factor matrix, the time window factor matrix, and the spatial measurement point factor matrix to obtain a three-dimensional reconstructed tensor with the same structure as the original third-order tensor. The absolute value of each element in this reconstructed tensor is then taken and arranged according to the index positions of the monitoring physical quantity, time window, and spatial measurement point, thus obtaining the coupling strength matrix. Each value in this matrix reflects the strength of the coupling effect between multi-source data at a specific monitoring physical quantity, a specific time window, and a specific spatial measurement point location.

[0042] The steps in the Tucker decomposition process described above have a clear logical progression. First, principal component directions in the three dimensions are extracted through three independent singular value decompositions, resulting in factor matrices. Then, the original tensors are modularly multiplied sequentially using the transposes of the factor matrices. The physical meaning of this step is to peel away the main change patterns in each dimension from the original data one by one; the remaining core tensors represent the pure interaction effects between the principal components. Finally, the coupling strength matrix, consistent with the spatial structure of the original data, is obtained by reconstructing and taking the absolute value. Throughout the process, the column vectors of the factor matrix maintain orthogonality and normalization, ensuring the uniqueness and numerical stability of the decomposition.

[0043] The Tucker decomposition method defined in this technology provides a systematic and repeatable mathematical path for the coupled analysis of multi-source monitoring data. By unifying the processing of physical quantities, time, and space through tensor decomposition, the resulting coupling strength matrix not only quantifies the interaction strength between displacement and physical quantities such as water level and moisture content, but also indicates which time windows and spatial measurement points this strong coupling occurs in. This three-dimensional joint analysis capability is not available in existing pairwise analysis methods. Furthermore, the factor matrix extraction based on singular value decomposition possesses optimal low-rank approximation properties, capable of compressing data dimensions while preserving the main energy of the original tensor to the greatest extent, thus effectively filtering out random noise interference. Subsequent screening of spatiotemporal sub-blocks dominated by displacement and strongly coupled with water level and moisture content relies on the coupling strength matrix and factor matrix output in this step. Therefore, this step is the key technical foundation for achieving deep fusion of multiple physical quantities and intelligent removal of interference signals in the entire early warning method. The standardized dataset processed by this step exhibits a significantly improved signal-to-noise ratio, providing high-quality input for high-precision early warning models in the cloud.

[0044] In another technical solution, step two, which involves selecting spatiotemporal sub-blocks that exhibit displacement as the dominant mode and simultaneously satisfy the coupling strength conditions with groundwater level and soil moisture content, includes: In the monitoring physical quantity factor matrix A, the row vectors corresponding to displacement, groundwater level, and soil moisture content are located, and the cosine similarity between the displacement row vector and the groundwater level row vector, and the cosine similarity between the displacement row vector and the soil moisture content row vector are calculated. In the core tensor G, the slice matrix corresponding to the displacement is extracted along the dimension of the monitoring physical quantity. The elements in the slice matrix are thresholded, and the index positions with the absolute value of the element greater than the preset threshold are marked as strong coupling candidate positions. Among the strong coupling candidate positions, the index positions that simultaneously satisfy the displacement-groundwater level cosine similarity exceeding the first threshold and the displacement-soil moisture content cosine similarity exceeding the second threshold are selected to form the target spatiotemporal sub-block. The original signal components corresponding to the third-order tensor at the target spatiotemporal sub-block positions are retained, and the signal components corresponding to the other positions are set to zero to obtain the signal after removing the transient impact sub-block.

[0045] In this technical solution, the row vectors corresponding to the three physical quantities—displacement, groundwater level, and soil moisture content—need to be located in the monitoring physical quantity factor matrix obtained in step two. The monitoring physical quantity factor matrix is ​​a two-dimensional matrix with the same number of rows as the number of monitored physical quantities and the same number of columns as the number of principal components selected in step two. Each row vector represents the load distribution of that physical quantity in each principal component direction. To quantify the directional similarity between displacement and groundwater level, the cosine similarity between the displacement row vector and the groundwater level row vector is calculated. The cosine similarity originates from the cosine value of the angle between vectors in geometry. Its calculation formula is the dot product of two vectors divided by the product of their respective moduli. The result ranges from -1 to +1; the closer the value is to one, the more consistent the directions of the two vectors are, meaning the two physical quantities exhibit highly similar change patterns in the principal component space. Similarly, the cosine similarity between the displacement row vector and the soil moisture content row vector is calculated. These two cosine similarity values ​​reflect the overall correlation strength between displacement and water level, and between displacement and moisture content, serving as global constraints for subsequent screening.

[0046] Next, a slice matrix corresponding to the displacement is extracted from the core tensor along the dimension of the monitored physical quantities. The core tensor is a multidimensional array, whose first dimension corresponds to the monitored physical quantity. However, after Tucker decomposition, this dimension is usually much smaller than the original number of physical quantity types because only the principal components are retained in step two. The coupling strength matrix is ​​constructed by reconstructing a three-dimensional array with the same structure as the original third-order tensor. Therefore, in this coupling strength matrix, a two-dimensional slice matrix corresponding to the displacement also exists along the dimension of the monitored physical quantities, with rows corresponding to time windows and columns corresponding to spatial measurement points. This slice matrix is ​​then thresholded, and the index positions where the absolute value of an element is greater than a preset threshold are marked as strong coupling candidate positions. The preset threshold can be set based on engineering experience, usually taking a certain percentile of all element values ​​in the coupling strength matrix, such as 80% or 90%, to ensure that only positions with sufficiently strong coupling are included as candidates. The physical significance of thresholding is that only those time windows and spatial measurement points exhibiting sufficiently strong coupling effects in the displacement direction can truly reflect the precursory information of slope instability.

[0047] After obtaining the strong coupling candidate locations, further global constraints based on cosine similarity need to be applied. Specifically, among the strong coupling candidate locations, index locations that simultaneously satisfy two conditions are selected. The first condition is that the cosine similarity between the displacement row vector and the groundwater level row vector exceeds a first threshold; the second condition is that the cosine similarity between the displacement row vector and the soil moisture content row vector exceeds a second threshold. The first and second thresholds can be determined based on statistical analysis of actual monitoring data, for example, both can be set to 0.7 or 0.8. These two conditions ensure that the selected locations are not only sufficiently strong in local coupling, but also that, from a global modal perspective, the trend of displacement changes is highly consistent with the trends of water level and moisture content changes. Index locations that simultaneously satisfy these three constraints constitute the target spatiotemporal sub-block.

[0048] Finally, the signal components in the original third-order tensor are processed according to the index position of the target spatiotemporal sub-block. Specifically, the original signal components corresponding to the target spatiotemporal sub-block position of the original third-order tensor are retained, while the signal components corresponding to all other positions in the original third-order tensor are set to zero. After this selective retention operation, only the data from time windows and spatial measuring points that simultaneously satisfy the strong coupling condition of "displacement-water level-water content" are retained from the original monitoring data. Transient impact signals such as construction vibrations, which do not satisfy this strong coupling condition, have their corresponding spatiotemporal sub-blocks removed. The resulting new data tensor is the signal after removing the transient impact sub-blocks, which can be further used for subsequent standardized dataset construction.

[0049] This technical solution, through joint analysis of the coupling strength matrix and factor matrix, identifies signal components truly related to seepage-stress coupling from the perspective of the intrinsic correlation of multiple physical quantities, achieving intelligent and precise elimination of transient impact interference. The technical effects of this step include: significantly improving the signal-to-noise ratio of monitoring data and avoiding false alarms caused by external interference such as construction vibration; simultaneously, by retaining only spatiotemporal sub-blocks closely related to the instability mechanism, the computational load of subsequent early warning models is correspondingly reduced, improving overall response efficiency. This step provides a crucial signal filtering method for forming a standardized dataset for uploading to the cloud, and is an important intermediate link connecting data fusion and the dual-drive early warning model.

[0050] In another technical solution, the calculation steps for the physical residual loss term in step three include: On the extended finite element discrete nodes, the displacement field u of the i-th node predicted by the Long Short-Term Memory network is extracted. i and pore water pressure field p i Substituting into the discrete form of the seepage-stress coupling control equation: K u u i +K c p i-f i =R u，i ; K c T (du) i / dt)+K p p i -q i =R p，i ; Among them, K u K is the stiffness matrix. c Let K be the coupling matrix. p Let f be the penetration matrix. i Let q be the nodal force vector of the i-th node. i Let R be the node flow vector of the i-th node. u，i Let R be the displacement field residual of the i-th node. p，i Let be the residual of the pore water pressure field at the i-th node; The physical residual loss term is: ; Where N is the total number of discrete nodes, and λ is the weighting coefficient for the difference in dimensions between the residual of the equilibrium displacement field and the residual of the pore water pressure field. The overall loss function for joint training is: ; Where L MSE α is the mean squared error loss term, and α is the physical residual penalty coefficient.

[0051] This technical solution further defines the calculation steps for the physical residual loss term in step three. This step uses the discrete form of the extended finite element seepage-stress coupling control equation as the loss function of the physical constraint embedded in the long short-term memory network, so that the displacement field and pore water pressure field predicted by the model simultaneously satisfy the mechanical laws and monitoring data.

[0052] The two discrete equations underlying the physical residual loss term are derived from Biot's consolidation theory, the foundation of porous elasticity, which describes the bidirectional coupling between fluid seepage and solid skeleton deformation. After spatial discretization, for the first... i Each node has a displacement field u i and pore water pressure field p i The following two discrete forms of the governing equations should be satisfied: The first equation is the force balance equation: K u ⋅u i +K c ⋅ p i -f i =R u,i Where: Ku The stiffness matrix, derived from the principle of virtual work in elasticity, has elements determined by material parameters such as Young's modulus and Poisson's ratio, as well as the element shape function. Its dimension is force per length; K c The coupling matrix is ​​derived from the Biot coupling coefficient and the strain-pressure coupling term, and its dimensions are the same as K. u The dimensionless quantity obtained by multiplying by the characteristic length is also force per length; f i R is a nodal force vector, comprising the soil's self-weight, support load, and external surcharge, with the dimension of force; u,i For displacement field residuals, the dimensions are the same as f. i The same. This equation represents the margin of the sum of the elastic force at the node, the contribution of coupled pressure, and the external load, when the predicted u... i and p i The residual is zero when the force balance is precisely satisfied.

[0053] The second equation is the flow continuity equation: Where: K c T The transpose of the coupling matrix converts the displacement rate of change into seepage volume change; is the derivative of the displacement field with respect to time, i.e., the deformation rate; K p The permeation matrix is ​​derived from Darcy's law and the mass conservation equation, and its dimension is volumetric flow rate per pressure gradient. q i The nodal flow vector includes boundary flows caused by rainfall infiltration, changes in groundwater level, etc., and is expressed in volumetric flow rate. R p,i The residual of the pore water pressure field, with dimensions of... q i The same. The equation means the rate of change of seepage volume at the node, the margin of the sum of seepage flow and external flow, and the residual is zero when the predicted value accurately satisfies the requirement of continuous flow.

[0054] The logical relationship between these two equations is as follows: the force balance equation controls the mechanical response of the solid skeleton, and the flow continuity equation controls the transport of the pore fluid. The two are coupled through the coupling matrix K. c These equations and their transposes are interconnected. Changes in displacement lead to alterations in pore volume, resulting in pore pressure. Changes in pore pressure, in turn, affect the effective stress, causing additional deformation. This is the essence of seepage-stress coupling. Satisfying both equations simultaneously at discrete nodes means that the prediction results conform to the physical mechanism.

[0055] Based on this, the physical residual loss term is defined. L physics The weighted average of the sums of squared residuals at all discrete nodes is expressed as: in: NTo expand the total number of discrete nodes in the finite element method; ∥⋅∥ 2 The square of the L2 norm is the sum of the squares of the vector components. l This is the dimensional balance weighting coefficient. Because R... u,i The dimension of force is... R p,i The dimension of is volumetric flow rate. Since the two have different dimensions, they cannot be directly added together. l Its function is to adjust the two residual terms to the same numerical order of magnitude. l The determination method is as follows: First, calculate the characteristic forces based on the material parameters and mesh scale of typical working conditions. F 0 and characteristic flow Q 0, take l = F 0 2 / Q 0 2 The characteristic force can be taken as the product of the elastic modulus and the characteristic area, and the characteristic flow rate can be taken as the permeability coefficient multiplied by the characteristic head difference and then multiplied by the characteristic area. l The specific values ​​depend on the material parameters and geometric dimensions of the actual project and should be obtained through preliminary finite element calibration for each specific foundation pit project.

[0056] The overall loss function for joint training is defined as: L total = L MSE + αL physics ; in: L MSE The mean squared error loss term is the mean squared error between the predicted value of the long short-term memory network and the actual value of the monitoring, which reflects the model's fitting accuracy to the time series data. α This is the physical residual penalty coefficient, used to adjust the relative importance of data-driven terms and physical constraint terms. α The value of should be such that, in the initial stage of training L physics and L MSE Being of the same order of magnitude, they are typically determined through a hyperparameter scan on the validation set, generally ranging from 0.001 to 0.1. When α When the size is too large, the model will become overly biased towards physical constraints and ignore monitoring data; when α When the size is too small, physical constraints have almost no effect. Reasonable... α It enables the model to fit the data without violating basic mechanical laws.

[0057] In practice, the high-precision early warning mode in step three activates the calculation of the physical residual loss term. During each forward propagation, the Long Short-Term Memory network outputs the predicted displacement field and pore water pressure field at each discrete node, while simultaneously reading the pre-calculated K from the extended finite element model. u K c K p f i , q i And the time step. Since these matrices and vectors depend only on material parameters, boundary conditions, and mesh generation, they can be determined and remain invariant before model training. Then, R is calculated according to the two discrete equations mentioned above. u,i and R p,i And thus obtain L physics During backpropagation, the gradient of the physical residual loss term is passed to the various weight parameters of the Long Short-Term Memory network through the chain rule, thereby guiding the network to learn a mapping relationship that conforms to both the data and physical laws.

[0058] This technical solution incorporates the discrete residuals of the Biot consolidation equation as an additional regularization term into the loss function, effectively imposing a physically constrained subspace on the data manifold, forcing the model output to fall within this subspace. This mechanism-data dual-driven approach maintains the nonlinear expressive power of deep learning while introducing prior knowledge from geotechnical mechanics, significantly improving the generalization ability and physical reliability of the early warning model beyond the training sample coverage. Furthermore, since the physical residual loss term is only enabled in high-precision early warning mode and not calculated in routine low-power inference mode, the allocation of computing resources is optimized while ensuring the reliability of early warnings.

[0059] In another technical solution, step three, which involves using conditional Wasserstein to generate an adversarial network that learns a Riemannian manifold of the latent space and performing geodesic interpolation on the Riemannian manifold, includes: The measured samples of various instability evolution paths were divided into initial creep stage, constant creep stage, accelerated creep stage, and critical instability stage according to time series. A conditional Wasserstein generative adversarial network was constructed. The generator input was random noise z and stage condition label c, and the output was the virtual sample of the corresponding stage. The discriminator output was the Wasserstein distance estimate. After the conditional Wasserstein generative adversarial network training converged, a set of random noise z1 and z2 was sampled for adjacent instability stages, and mapped to the latent space by the generator to obtain latent vectors w1 and w2. The geodesic path between the latent vectors was calculated. ; Where exp and log are the exponential and logarithmic mappings on the Riemannian manifold, respectively; Interpolation progress parameters from the starting point to the ending point ( When =0, it corresponds to w1. When =1 corresponds to w2); in Uniform sampling within the value range M Each interpolation point is used to decode the latent vector corresponding to each interpolation point through a generator, resulting in a virtual instability evolution sample that smoothly transitions between adjacent instability stages.

[0060] In this technical solution, during implementation, the measured samples of various types of instability evolution paths are first divided into four stages according to time series. This division is based on the typical three-stage or four-stage curve in creep theory, namely, the initial creep stage is characterized by a gradually decreasing deformation rate, the isochronous creep stage by a constant deformation rate, the accelerated creep stage by a significantly increased deformation rate, and the critical instability stage by a sharp increase in deformation rate until failure. The division method can be automatic identification using a threshold of deformation rate or the zero-crossing point of the second derivative, or it can be manually labeled by geotechnical engineers based on the monitoring curves. After the division is completed, each measured sample is assigned a stage condition label.

[0061] Construct a conditional Wasserstein generative adversarial network, whose generator accepts a random noise vector. z and stage condition labels c As input, the generator outputs virtual samples corresponding to the specified stage. The discriminator accepts samples and conditional labels, outputting a real number as an estimate of the Wasserstein distance. During training, the generator aims to minimize the Wasserstein distance, while the discriminator aims to maximize it; the two learn adversarially. To ensure the Lipschitz continuity of the Wasserstein distance, gradient penalty or weight clipping is typically used. After network training converges, the generator can map any point in the latent space to a virtual sample that conforms to the characteristics of the specified stage, and the generated samples corresponding to similar points in the latent space are also similar in the data space.

[0062] After the Conditional Wasserstein Generative Adversarial Network (CGN) training converges, a set of random noise is sampled for each adjacent instability stage. Let the initial creep stage sampling noise be... z 1. Sampling noise during constant-rate creep stage z 2. Map the generator to the latent space to obtain the latent vector. w 1 and w 2. The latent space refers to the implicit representation space of the generator's input space after encoding. It is generally assumed that the data distribution within this space follows a simple prior, such as a Gaussian distribution. However, since the generator is a nonlinear mapping, the image of the data manifold in the latent space is not a flat Euclidean space, but rather has a Riemannian manifold structure. To compute the connectivity... w 1 and w The geodesic path of 2 is explained by introducing exponential and logarithmic mappings. The exponential mapping exp... w (v ) indicates from point w Starting along the tangent vector v A point on a manifold reachable by traveling a unit distance; logarithmic mapping log w ( w ′) is the inverse operation, returning from w arrive w The tangent vector of ′. The geodesic path is represented as: w ( t )=exp w1 ( t ⋅log w1 ( w 2)); among which t The interpolation progress parameter ranges from 0 to 1. t =0 corresponds to w 1, t =1 corresponds to w 2. In practical calculations, since the analytical form of the manifold is unknown, exponential and logarithmic mappings cannot be directly calculated. A common approximation method is to perform linear interpolation in the latent space and then map it back to the data space through a generator. However, this approximation is only effective when the manifold is approximately flat. To obtain more accurate geodesics, a metric tensor can be defined in the latent space, and numerical methods for solving the geodesic differential equations can be used, such as Euler-Lagrange discretization or a neural network-based geodesic solver. Another simplified approach is to perform linear interpolation directly in the latent space and then use the output of the discriminator as a correction signal, iteratively optimizing the process to keep the points on the interpolation path near the real data manifold. Within the scope of this invention, those skilled in the art can choose any of the above-mentioned numerical methods for geodesic calculation.

[0063] In determining t After taking the value interval [0,1], uniform sampling is performed within this interval. M There are 10 interpolation points, and each interpolation point corresponds to a progress parameter. t m Each t m Substituting into the geodesic path expression, we obtain the corresponding latent vector. w ( t m Then, each latent vector is mapped back to the data space through the decoding part of the generator to generate virtual samples. These virtual samples represent the unstable evolution pattern with a smooth transition from the initial creep stage to the isochronous creep stage. By repeating the above process and performing geodesic interpolation on all adjacent stages, an expanded training set covering the entire unstable evolution process can be obtained.

[0064] The logical relationship between the above steps is as follows: First, a limited number of measured samples are organized into a labeled dataset through stage division. Then, a latent manifold for learning the data distribution of the adversarial network is generated using Conditional Wasserstein. Finally, transition samples are generated on the manifold using geodesic interpolation. The core advantage of geodesic interpolation is that, since the instability evolution process is physically continuous and nonlinear, the transition form between adjacent stages does not necessarily change linearly in Euclidean space. However, since geodesics move along the curvature of the manifold, the interpolated samples generated are more consistent with the actual physical evolution. For example, from the constant-rate creep stage to the accelerated creep stage, the rate of change of deformation itself changes. Geodesic interpolation can capture this acceleration change, while linear interpolation will lose this feature.

[0065] This technical solution, by introducing Riemannian manifolds and geodesic interpolation, generates virtual samples that maintain the physical continuity and geometric consistency of the instability evolution process, significantly expanding the effective information content of the training set. This method is particularly suitable for low-probability event scenarios such as slope instability, where complete instability cases are extremely limited, and transitional samples at different evolution rate stages are even more difficult to obtain through actual measurements. Training the early warning model together with the generated virtual samples and measured samples can effectively alleviate overfitting and improve the model's ability to identify different instability stages, especially critical instability precursors. Simultaneously, because conditional generative adversarial networks can generate samples according to stage labels, the number of samples at each stage in the training set is balanced, avoiding the problem of model bias due to data imbalance.

[0066] In another technical solution, the switching steps between low-power inference mode and high-precision warning mode in step three include: In low-power inference mode, the mechanism-data dual-driven hybrid early warning model loads a lightweight branch of the Long Short-Term Memory (LSTM) network. This lightweight branch is obtained from the original LSTM network after channel pruning and weight quantization, and the physical residual loss term is not calculated during inference. When the cloud server detects that the displacement-groundwater level coupling strength and the displacement-soil moisture content coupling strength are monotonically increasing within k consecutive time windows, and the coupling strength of the current time window exceeds their respective preset thresholds, the mode switch is triggered. The switching process includes: loading the complete LSTM network weights from the model storage area to replace the lightweight branch weights, and activating the computation node corresponding to the physical residual loss term in the computation graph, so that the model enters the high-precision early warning mode.

[0067] In this technical solution, under low-power inference mode, a lightweight branch of the Long Short-Term Memory (LSTM) network is loaded into the mechanism-data dual-driven hybrid early warning model. This lightweight branch is obtained by performing channel pruning and weight quantization on the complete LTM network. Channel pruning is a model compression technique that evaluates the importance of each neuron or convolutional channel and removes channels and their connection weights that contribute less, thereby reducing the number of model parameters and computational cost. Weight quantization converts floating-point weights into low-bit integer representations, such as from 32-bit floating-point to 8-bit integers, which can significantly reduce memory usage and accelerate inference operations. After pruning and quantization, the number of parameters in the lightweight branch is typically only one-tenth to one-fifth of that in the complete model, resulting in a several-fold increase in inference speed. In this mode, the physical residual loss term is not calculated during model inference because calculating the physical residual requires accessing data such as the finite element mesh and stiffness matrix and performing matrix operations, which incurs significant computational overhead. Simultaneously, the output accuracy in low-power mode meets the needs of daily stability monitoring because when the slope is in a stable or weakly coupled state, physical quantities such as displacement and water level change gradually, and the lightweight branch is sufficient to capture their main trends.

[0068] The trigger condition for mode switching is determined by the coupling strength matrix monitored in real time by the cloud server. Specifically, the cloud server monitors two key coupling strength indicators: displacement-groundwater level coupling strength and displacement-soil moisture content coupling strength. Both indicators are extracted from the coupling strength matrix. Mode switching is triggered when the following three conditions are met: First, within multiple consecutive time windows, the displacement-groundwater level coupling strength shows a monotonically increasing trend, meaning the coupling strength value of each time window is greater than the value of the previous time window; second, within the same consecutive time window sequence, the displacement-soil moisture content coupling strength also shows a monotonically increasing trend; third, both coupling strength values ​​in the current time window exceed their respective preset thresholds. The number of "multiple consecutive time windows" is denoted as k, and k should be greater than or equal to 3 to exclude false triggers caused by occasional fluctuations. The preset thresholds can be determined based on the statistical distribution of coupling strength under historical stable operating conditions, such as two or three times the historical average, or set based on engineering experience. The determination of the monotonically increasing trend can be verified using the sign of the linear regression slope within the sliding window or the sign sequence of the difference between adjacent windows.

[0069] When the above conditions are met, the system determines that the slope is transitioning from a stable state to a potentially unstable state and enters a high-precision early warning mode. The mode switching process includes the following two operations: First, the complete Long Short-Term Memory (LSTM) network weights are loaded from the model storage area to replace the lightweight branch weights in the current inference engine. The model storage area can be persistent storage or a distributed cache on a cloud server, pre-saving fully trained, unpruned, and unquantized complete model parameters. The loading process needs to ensure atomicity, that is, lightweight branch inference is still used before the replacement is completed, and a seamless switch to the complete model is performed after the replacement is completed. Second, the computation nodes corresponding to the physical residual loss term are activated in the computation graph of the deep learning inference framework. In low-power mode, these computation nodes are marked as skipped or directly short-circuited and do not participate in forward computation. During switching, the execution path of the computation graph is modified so that the data flows through the physical residual computation nodes, thereby outputting the physical residual loss term at each forward propagation, and this term can be used as an auxiliary output or for subsequent gradient calculations. After activating the physical residual node, the model enters a complete high-precision early warning mode. Its output includes not only the instability probability and early warning level, but also the physical residual information of the displacement field and pore water pressure field, which can be reviewed by engineers.

[0070] After mode switching, the system will remain in high-precision warning mode for a preset duration until it detects that the coupling strength has decreased for multiple consecutive time windows and fallen below the safety threshold. At this point, it can automatically switch back to low-power mode. The process of switching back to low-power mode is the reverse of the above process: unload the full model weights, reload the lightweight branch weights, and disable the physical residual calculation node.

[0071] The dual-mode switching mechanism proposed in this technical solution utilizes the dynamic evolution trend of the coupling strength matrix as the switching basis, rather than a single absolute threshold. This is because a monotonically increasing trend is a better indicator of the system entering an unstable evolution process than a single threshold crossing. This trend-based judgment reduces erroneous switching caused by instantaneous disturbances. Simultaneously, through lightweight design using model pruning and quantization, daily power consumption and computational load are significantly reduced, making it suitable for long-term continuous operation. When high-precision analysis is truly needed, it can automatically and quickly switch to the complete model and activate physical constraints, providing the most reliable early warning information. This adaptive switching strategy balances the economy of all-weather monitoring with the accuracy of emergency response, and is one of the key innovations of this invention.

[0072] In another technical solution, the deep deformation monitoring in step one is achieved by deploying an array-type microelectromechanical inclinometer. The packaging structure of this inclinometer includes: MEMS accelerometer chips are soldered onto a flexible printed circuit board, which is fixed to the inner wall of a glass fiber reinforced epoxy resin composite tube by rigid partitions arranged at equal intervals. The outer wall of the composite tube is formed by spirally winding epoxy resin-impregnated glass fiber cloth along the tube axis at an angle of 45° to 55°. Both ends of the composite tube are sealed with stainless steel end caps and O-rings. The temperature self-compensation circuit includes a thermistor from the same batch as the MEMS accelerometer chip. The thermistor and the MEMS accelerometer chip are arranged adjacent to each other on the same flexible printed circuit board. The temperature self-compensation circuit corrects the accelerometer output in real time through a pre-calibrated temperature-drift curve. The redundant MEMS sensing unit consists of two sets of MEMS accelerometer chips arranged orthogonally on the same rigid partition.

[0073] This technical solution further defines the specific implementation method of deep deformation monitoring in step one by deploying an array of microelectromechanical inclinometers.

[0074] The inclinometer's overall structure utilizes a flexible printed circuit board (PCB) as the carrier for its sensing core. MEMS accelerometer chips are soldered onto the PCB using surface mount technology. The PCB itself is flexible, thin, and lightweight, allowing it to adapt to bending deformations during installation while maintaining the relative positional accuracy between the accelerometer chips. To prevent wrinkles or displacement of the PCB due to soil pressure or temperature changes during long-term use, it is fixed to the inner wall of a glass fiber reinforced epoxy resin composite tube using equally spaced rigid partitions. The rigid partitions are disc-shaped, matching the inner diameter of the composite tube, with a through-hole in the center for cable passage, and their edges are bonded or secured to the inner wall of the composite tube. The rigid partitions are evenly spaced along the tube's axial direction; the spacing can be determined based on the inclinometer's spatial resolution requirements, for example, one partition every 0.5 meters or 1 meter. The PCB passes through fixing slots or clips on each rigid partition and is constrained to a predetermined position on the inner wall of the tube.

[0075] The composite tube, serving as the outer shell of the inclinometer, functions as a waterproof, pressure-resistant, impact-resistant material, and protects the internal circuitry. The composite tube uses glass fiber reinforced epoxy resin as the matrix material, with epoxy resin-impregnated glass fiber cloth spirally wound along the tube's axial direction on its outer wall. The spiral winding angle is controlled between 45 and 55 degrees. This angle range originates from the optimal winding angle theory in composite material mechanics. For thin-walled tubes subjected to axial compression and circumferential pressure, when the winding angle is near 55 degrees, the axial stiffness and circumferential stiffness of the tube are balanced, resulting in the strongest buckling resistance. Simultaneously, the spiral winding structure allows the tube to convert radial earth pressure into tensile stress in the fibers, preventing brittle fracture. Both ends of the composite tube are sealed with stainless steel end caps and O-rings. The stainless steel end caps are connected to the ends of the composite tube via threads or snaps, and an O-ring is embedded between the end cap and the tube body. After tightening, a reliable static seal is formed, preventing groundwater or sediment from seeping into the tube.

[0076] A temperature self-compensation circuit is a key design feature for addressing the temperature drift problem in MEMS accelerometers. The output zero point and scale factor of a MEMS accelerometer chip drift with temperature changes, primarily due to the thermal expansion of the silicon micromechanical structure and the temperature dependence of transistor parameters in the detection circuit. To correct this drift in real time, this invention places a thermistor manufactured in the same batch as the MEMS accelerometer chip on a flexible printed circuit board. "Same batch" means the thermistor and accelerometer chip underwent the same process conditions and packaging, resulting in highly consistent temperature response characteristics. The thermistor is placed adjacent to the accelerometer chip on the same flexible printed circuit board, with the distance between them as small as possible, typically within a few millimeters, to ensure that they sense approximately the same temperature. The temperature self-compensation circuit works as follows: During the production stage, an inclinometer is placed in a temperature-controlled chamber, and the output values ​​of the accelerometer and the resistance values ​​of the thermistor are collected at multiple preset temperature points. A temperature-drift curve is then fitted, describing the functional relationship between the accelerometer's zero-point output and temperature. In actual monitoring, the thermistor measures the ambient temperature in real time. The compensation circuit queries the temperature-drift curve based on the current temperature value, calculates the drift at the current temperature, and subtracts it from the accelerometer's original output to obtain the corrected acceleration value. This process can be hardware compensation in the form of analog circuits or software compensation after digital sampling.

[0077] The design of redundant MEMS sensing units aims to improve the reliability and measurement dimensionality of inclinometers. Two sets of MEMS accelerometer chips are arranged orthogonally on the same rigid partition. Orthogonal directions typically refer to two sensing axes that are perpendicular to each other; for example, one set of chips measures the tilt angle change in the X direction, and the other measures the tilt angle change in the Y direction. In inclinometers, it is usually necessary to measure the tilt angles in two orthogonal directions to calculate the spatial displacement vector. Therefore, the redundancy design not only provides backup but also enables dual-axis measurement. Each set of chips may contain multiple accelerometers to form differential or averaged outputs. When one set of chips fails, the other set can still maintain unidirectional monitoring capability. Although it cannot fully reconstruct the spatial displacement, it can still provide an early warning signal. The redundancy design also includes independent power supply and signal transmission lines for the two sets of chips to avoid common-cause failures.

[0078] The assembly sequence of the entire inclinometer is as follows: First, the MEMS accelerometer chip, thermistor, and related resistors and capacitors are soldered onto a flexible printed circuit board and preliminarily tested. Then, rigid partitions are fixed to the inner wall of the composite tube at the designed spacing. The flexible printed circuit board is then passed through the fixing slots of each rigid partition and positioned using a small amount of glue. The temperature self-compensation circuit and signal conditioning circuit are similarly fixed in appropriate positions inside the tube. Finally, stainless steel end caps are fitted with O-rings and tightened at both ends of the composite tube, and the cables are led out. After installation, an overall sealing test and temperature calibration are performed, generating a temperature-drift curve and storing it in an external data acquisition instrument or edge gateway.

[0079] This technical solution utilizes a spirally wound glass fiber reinforced epoxy resin composite tube, which reduces weight and improves corrosion resistance while maintaining strength; optimized winding angle enhances the tube's resistance to lateral pressure; close-range arrangement of thermistors from the same batch enables precise temperature self-compensation; and the orthogonal dual-group redundant design improves reliability and enables biaxial measurement. These structural features collectively enhance the stability and measurement accuracy of the inclinometer in long-term monitoring of foundation pit slopes.

[0080] In another technical solution, the implementation steps of the unified data interface in step one include: A predefined data dictionary covering displacement, pressure, moisture content, and vibration parameters is used. Each entry in the data dictionary includes a data code, physical quantity name, unit, data type, and value range. The edge computing gateway allocates independent communication ports and message queues for each type of sensor. Each sensor encapsulates data packets according to the data dictionary format. The header of the data packet contains the sensor type code, accuracy level code, and GNSS timing timestamp, while the body of the packet contains monitoring data in key-value pair format. The data parsing module of the edge computing gateway reads the sensor type code from the packet header, queries the data dictionary to obtain the key-value pair parsing template corresponding to that type of sensor, extracts the monitoring data according to the parsing template, and writes it into the time-series database.

[0081] This technical solution further defines the implementation steps of the unified data interface in Step 1. The data dictionary is the core metadata definition of the unified data interface. A predefined data dictionary covering displacement, pressure, moisture content, and vibration parameters is stored in the configuration file of the edge computing gateway in a structured table format, or it can be stored in the cloud and synchronized to each edge gateway. Each entry in the data dictionary contains the following fields: data code, physical quantity name, unit, data type, and value range. The data code is a unique integer value or short string used to identify a specific physical quantity type. For example, surface displacement is encoded as "SD01", deep displacement as "SD02", groundwater level as "GW01", soil moisture content as "WC01", rainfall intensity as "RI01", and construction vibration acceleration as "VV01". The physical quantity name uses a human-readable string, such as "surface displacement_X direction". The unit field is filled in according to the International System of Units (SI) or commonly used engineering units, such as "millimeters", "meters", "kilopascals", "percentage", "millimeters per hour", and "meters per square second". The data type field specifies the data storage type of the physical quantity, such as "floating-point", "integer", or "boolean". The value range field gives the reasonable range of the physical quantity under normal operating conditions. For example, the value range of surface displacement can be "-100 to 100" millimeters. The edge gateway can use this range to perform preliminary outlier detection and filtering.

[0082] The edge computing gateway assigns independent communication ports and message queues to each type of sensor. Communication ports can be physical serial port numbers, such as COM1 and COM2, or network port numbers, such as TCP ports 5001 to 5010. Each sensor is assigned to a corresponding port based on its sensor type code when it connects to the gateway. Message queues are used to buffer raw data packets received from sensors of that type to handle sudden data spikes. Message queues can be implemented using a first-in, first-out (FIFO) circular buffer, with each queue having a maximum length; when the length is exceeded, the oldest or newest message is discarded according to a strategy. This independent queue design avoids data blocking between different types of sensors; for example, a large amount of data from a high-frequency vibration sensor will not affect the data processing of a low-frequency displacement sensor.

[0083] Before transmitting data, each sensor must encapsulate a data message according to the data dictionary format. The data message structure consists of two parts: a header and a body. The header contains three fields: sensor type code, accuracy level code, and GNSS timing timestamp. The sensor type code directly references the data code defined in the data dictionary and is 1 to 2 bytes long. The accuracy level code is a predefined code used to identify the sensor's measurement accuracy level; for example, "01" represents high accuracy, "02" represents medium accuracy, and "03" represents low accuracy. Sensors with different accuracies will be assigned different observation noise covariances in the subsequent unscented Kalman filtering. The GNSS timing timestamp uses Coordinated Universal Time (UTC) format, typically a whole second plus milliseconds or microseconds, and originates from the GNSS timing module integrated into the sensor or is uniformly timed by an edge computing gateway. The message body uses a key-value pair format, i.e., a sequence of "key name = key value," with multiple key-value pairs separated by delimiters, such as "CH1=2.35;CH2=2.40." The key name can be the measurement point number or the channel number, and the key value is the corresponding monitoring data value. This key-value pair format has good scalability, allowing different sensors to flexibly organize the message body according to their own number of channels. Before sending, the sensor converts the entire message into a byte stream according to a predetermined encoding rule.

[0084] The data parsing module of the edge computing gateway is responsible for receiving data streams from each port and performing parsing operations. The parsing module first reads the sensor type code in the message header, uses it as a query keyword to search the data dictionary, and obtains the key-value pair parsing template corresponding to that sensor type. The parsing template describes the separator rules for key-value pairs in the message body, the meaning of the physical quantity corresponding to each key name, and the data type conversion method. For example, for a displacement sensor, the parsing template might indicate that "CH1 corresponds to displacement in the X direction, CH2 corresponds to displacement in the Y direction, and the value needs to be divided by 1000 to convert to meters." Based on the parsing template, the module reads the key-value pairs in the message body line by line, extracts the monitoring data, and performs necessary unit and data type conversions. The converted data, along with the accuracy level code and GNSS timestamp extracted from the message header, is written to the time-series database. The time-series database organizes data by physical quantity type, spatial measurement point number, and timestamp, facilitating subsequent tensor construction and fusion analysis. If the sensor type code in the message header does not exist in the data dictionary, or if the message body format does not match the parsing template, the parsing module marks the message as invalid and records the error log, while continuing to process subsequent messages to ensure that the failure of a single sensor does not affect the operation of the entire system.

[0085] This technical solution achieves plug-and-play sensor access capabilities by pre-defining a unified data dictionary and key-value pair message format. The data dictionary centrally manages the metadata of all physical quantities, making the parsing process entirely configuration-driven rather than hard-coded. Adding a new sensor type only requires adding the corresponding entry to the data dictionary, without modifying the parsing module code. Simultaneously, the addition of accuracy level encoding to the message header provides a basis for dynamically adjusting the observation noise covariance of subsequent unscented Kalman filtering, enabling the fusion algorithm to automatically adapt to sensors with different accuracy levels. An independent message queue design ensures non-blocking processing of multi-source data streams, improving the real-time performance and throughput of the edge computing gateway. These designs collectively constitute a standardized, scalable, and self-describing unified data interface, laying the data foundation for the spatiotemporal alignment and fusion of multiple physical quantities in step two.

[0086] In another technical solution, step two, which involves using unscented Kalman filtering to perform optimal fusion estimation of multiple physical quantities, includes: Constructing the state vector: ; Where d x d y d z For three-dimensional displacement components, v x v y v z The displacement rate components are three-dimensional; the state transition equation adopts a constant velocity model; in the observation equation, the observation vector of the GNSS receiver is [d]. x d y d z ] T The observation vector of the ground-based synthetic aperture radar is transformed into the BIM coordinate system by the line-of-sight projection matrix as d. x d y d z A linear combination; by querying the preset accuracy level-variance mapping table using the accuracy level labels attached to each sensor message, the observation noise covariance matrix R of each sensor is obtained. k The diagonal element values; calculating the Sigma point set and its weights during the time update phase, and using R during the measurement update phase. k To measure the noise covariance, calculate the Kalman gain and perform state updates.

[0087] This technical solution further defines the step of using unscented Kalman filtering to perform optimal fusion estimation of multiple physical quantities in step two. Constructing the state vector is the first step of filtering. In this invention, the state vector x k It includes the displacement components and displacement rate components of a certain measuring point on the foundation pit slope in three-dimensional space, specifically in the form of x. k =[ d x, d y , d z , v x , v y , v z ] T ,in d x , d y , d z These represent the displacement components of the measuring point along the X, Y, and Z directions of the BIM coordinate system, respectively, with the dimension being length; v x , v y , v z These represent the displacement rate components in three directions, with dimensions of length per time. The state vector has a dimension of 6. (Subscripts are used for indexing.) k This indicates the sequence number of the discrete time step, corresponding to the synchronization time window of the edge computing gateway.

[0088] The state transition equation is used to describe the state transition from time 1 to 2. k -1 until time... k The system state evolution law is described. This invention adopts a constant velocity model, which assumes that the velocity in each direction remains constant within adjacent time windows, and the displacement changes linearly with the velocity. The constant velocity model originates from uniform linear kinematics and is one of the most basic kinematic models. Its discretized form is: d x,k = d x,k−1 + v x,k−1 ⋅Δ t , d y,k = d y,k−1 + v y,k−1 ⋅Δ t , d z,k = d z,k−1 + v z,k−1 ⋅Δ t And the velocity components in all directions remain constant. Where Δ t Let x be the time interval between two adjacent time windows. This can be represented in matrix form as x. k =F⋅x k−1 +w k−1Where F is the state transition matrix. The specific structure of F is as follows: in the rows corresponding to the first three-dimensional displacements, the position coefficient corresponding to its own displacement is 1, and the position coefficient corresponding to the velocity in the same direction is Δ. t In the row corresponding to the three-dimensional velocity, the position coefficient corresponding to its own velocity is 1, and the other position coefficients are 0. k−1 The process noise represents the velocity fluctuations that the model fails to describe, and its covariance matrix is ​​denoted as Q.

[0089] The observation equation describes the relationship between sensor observations and the state vector. This invention employs at least two types of surface displacement sensors: a BeiDou GNSS receiver and a ground-based synthetic aperture radar. For the GNSS receiver, it directly outputs the three-dimensional displacement observations of the measurement point in the BIM coordinate system; therefore, the observation vector z... k GNSS =[ d x obs , d y obs , d z obs ] T Observation matrix H GNSS It is a 3×6 matrix, with the elements being: the first three columns are 3×3 identity matrices, and the last three columns are all zeros. The observation equation is z. k GNSS =H GNSS ⋅x k +v k GNSS , where v k GNSS This refers to the observation noise of GNSS.

[0090] For ground-based synthetic aperture radar, the measurement is the line-of-sight deformation, i.e., the displacement projection along the radar beam direction. Let the unit vector of the radar line-of-sight direction in the BIM coordinate system be [cos...]. i x cos i y cos i z ], then the radar observation value z k GBSAR =cos i x ⋅ d x +cos i y ⋅ d y +cos i z ⋅ dz This is a scalar observation, with the observation matrix H. GBSAR This is a 1×6 row vector: the first three elements are the direction cosines of the line of sight, and the last three elements are zero. The observation equation is... z k GBSAR =H GBSAR ⋅x k + v k GBSAR When both GNSS and radar observations are present, the two observation vectors can be stacked into a composite observation vector, and the corresponding observation matrices are also stacked vertically.

[0091] The observation noise covariance matrix R of each sensor k The accuracy level is determined using the accuracy level tag attached to the sensor message. In the unified data interface step one, each sensor appends an accuracy level code to the message header. The edge computing gateway pre-stores an accuracy level-variance mapping table, in which each accuracy level corresponds to a set of variance values. For example, the displacement variance corresponding to high-precision GNSS can take a smaller value, while the variance corresponding to medium-precision GNSS is larger. For different sensor types, the variance dimension in the mapping table is different: GNSS uses a three-dimensional diagonal variance matrix, while ground-based synthetic aperture radar uses a scalar variance. When the gateway parses sensor data, it reads the accuracy level code from the message header, queries the mapping table to obtain the observation noise variance of the sensor at the current time, and fills it into R. k The diagonal elements. This process enables adaptive noise estimation, allowing the filtering algorithm to dynamically respond to changes in the accuracy of each sensor.

[0092] The recursive process of unscented Kalman filtering consists of two main steps: time update and measurement update. In the time update phase, the state estimate from the previous time step is used... The covariance matrix P k−1 This generates a set of Sigma points and their weights. The number of Sigma points is typically 2^3. n +1, where n Since the state vector has a dimension of 6, there are a total of 13 Sigma points. Each Sigma point is propagated to the current time step through the state transition equation, resulting in a predicted set of Sigma points. Then, the weighted summation is used to calculate the predicted state value. and predicted covariance P k∣k−1 During the measurement update phase, the predicted observation values ​​are calculated based on the predicted Sigma point set and the observation function. Observation covariance P zz And the cross-covariance P between state and observation xz Kalman gain K k Calculated as K k =P xz ⋅Pzz −1 Finally, the actual observation vector z at the current time is used. k Update state estimates and covariance: The advantage of unscented Kalman filtering over traditional extended Kalman filtering is that it does not require the calculation of the Jacobian matrix. It approximates the distribution after nonlinear transformation through deterministic sampling, and has higher accuracy for nonlinear observation equations such as line-of-sight projection of ground-based synthetic aperture radar.

[0093] This technical solution unifies linear GNSS observations and nonlinear radar observations through unscented Kalman filtering, and adaptively adjusts the noise covariance using an accuracy level-variance mapping table, resulting in superior fusion results. In areas with good satellite signal, GNSS observations are prioritized, while radar observations supplement the filter updates in signal-obstructed areas, ensuring continuous high-precision positioning around the clock. The displacement estimate from the filtered output serves as the foundational data for constructing the third-order tensor in step two, providing accurate and smooth displacement field inputs for subsequent coupling strength analysis and early warning models.

[0094] In another technical solution, in step one, surface displacement monitoring simultaneously deploys a BeiDou / GNSS receiver and a ground-based synthetic aperture radar. In areas where satellite signals are obscured, the deformation sequence along the line of sight is obtained using ground-based synthetic aperture radar interferometry. Pixels with amplitude deviation indices less than a preset threshold in the time-series interferogram are selected as candidate points for permanent scatterers. After time-dimensional phase unwrapping of the candidate points, a differential phase model between adjacent points is established. ; Where Δv is the deformation rate difference, T is the time baseline, B⊥ is the vertical baseline, R is the slant distance, θ is the incident angle, Δh is the elevation error difference, and Δϕ is the vertical distance. atm For atmospheric delay phase difference, Δϕ noise This is the noise phase difference; The atmospheric delay phase component of each permanent scatterer candidate point is separated and estimated by spatiotemporal filtering. The atmospheric delay phase component is then subtracted from the total phase to obtain the corrected deformation phase sequence.

[0095] This technical solution further defines the specific implementation method of surface displacement monitoring in step one. Specifically, in the surface displacement monitoring of the foundation pit slope, a BeiDou GNSS receiver and a ground-based synthetic aperture radar are deployed simultaneously. The BeiDou GNSS receiver is deployed at the top of the slope and on each level of the ramp, requiring a wide and unobstructed field of view, to provide an absolute coordinate reference and high-precision three-dimensional displacement observation. The ground-based synthetic aperture radar is erected in a stable area opposite the slope, with scanning tracks arranged along the slope direction to perform continuous interferometric imaging of the slope surface. The two devices complement each other: in areas with good satellite signals, GNSS observation is the primary method; in areas where satellite signals are obstructed, such as the bottom of a deep foundation pit or the shady side of a high slope, the interferometric measurement results of the ground-based synthetic aperture radar are used to obtain the line-of-sight deformation sequence.

[0096] The basic principle of ground-based synthetic aperture radar interferometry is to extract phase difference changes from complex radar images of the same scene acquired at different times to reflect minute surface deformations. Ground-based radar wavelengths are typically on the centimeter scale, and the interferometric phase is highly sensitive to deformation. To extract reliable deformation information from the original interferogram, stable scattering points need to be selected as measurement benchmarks. This invention uses the amplitude deviation index as a screening criterion for permanent scatterer candidate points. The amplitude deviation index is defined as the ratio of the standard deviation to the mean of the amplitude values ​​of the same pixel in a time series. This index originates from the statistical concept of the coefficient of variation, and its physical meaning is: when the scattering characteristics of a pixel are stable, the amplitude changes little over time, and the deviation index is close to zero; when affected by factors such as vegetation growth and changes in surface cover, the amplitude fluctuates drastically, and the deviation index is large. A preset threshold is usually set, such as 0.2 or 0.25, and only pixels with an amplitude deviation index less than this threshold are retained as permanent scatterer candidate points. These candidate points typically correspond to ground features with stable radar reflection characteristics, such as exposed rock, concrete structures, and metal components.

[0097] Temporal phase unwrapping is performed on the selected permanent scatterer candidate points. Phase unwrapping is a crucial step in interferometry because the phase measured by radar is wrapped within the principal value interval, and it is necessary to recover the true integer phase difference. Temporal phase unwrapping utilizes the phase changes between adjacent images in the time series, and obtains the absolute phase by accumulating and summing them. Commonly used algorithms include minimum cost flow and branching methods. The unwrapped phase includes deformation phase, atmospheric delay phase, residual phase caused by elevation error, and noise phase.

[0098] To separate the deformation phase, a differential phase model between adjacent points is established. This model utilizes the strong correlation between atmospheric delays between spatially adjacent permanent scatterers, and weakens the atmospheric influence by subtracting from adjacent points. Let two adjacent candidate permanent scatterer points be considered. i and j The differential phase between them is expressed as: ; Where: Δ ϕ ij The untangling phase difference between the two points; l For radar wavelength; Δ v The difference in deformation rate between the two points; T The time baseline is the time interval between two observations; B ⊥ The vertical baseline is the vertical distance between two radar observation positions. R Slant range, which is the distance at which the radar reaches the target; i Δ is the angle of incidence; h The difference in elevation between the two points; Δ ϕ atm,ij The atmospheric delay phase difference between two points; Δ ϕ noise,ij The first term in the equation is the phase difference caused by the deformation rate, which is proportional to the time baseline and is the signal to be extracted. The second term is the residual phase caused by the elevation error, which is proportional to the vertical baseline and can be estimated and subtracted. The last two terms are interference terms.

[0099] Atmospheric delay phase components of candidate permanent scatterers are separated and estimated using spatiotemporal filtering. Atmospheric delay primarily arises from temporal and spatial variations in water vapor content, exhibiting the following spatiotemporal characteristics: spatially, it possesses low-pass smoothness, meaning adjacent points have similar water vapor content; temporally, it exhibits high-pass randomness, meaning water vapor changes rapidly. Based on these characteristics, the following methods can be used to separate the atmospheric phase: first, low-pass filtering is applied to the spatial domain of each interferogram to remove high-frequency deformation signals, retaining the large-scale phase components as estimates of atmospheric delay; or, high-pass filtering is applied to the time series of each pixel in the temporal domain to separate the components that change rapidly over time as atmospheric delay. A more commonly used method is joint spatiotemporal filtering, which involves first performing principal component analysis on the time series to identify the dominant spatially correlated components, subtracting them as atmospheric delay components. The specific steps are as follows: All interferograms are stacked into a three-dimensional spatiotemporal dataset. Singular value decomposition is performed on this dataset. The first few principal components typically represent the dominant modes of atmospheric delay. The phase of each pixel is projected onto these principal components to calculate atmospheric estimates. These estimates are then subtracted from the total phase to obtain the corrected deformation phase sequence. The corrected deformation phase is then converted into the displacement of the slope surface through the line-of-sight projection relationship and input along with GNSS observations into the unscented Kalman filter in step two.

[0100] This technical solution employs amplitude deviation index screening to ensure the stability of selected scattering points over long time series. Simultaneously, the differential phase model between adjacent points effectively suppresses common-mode errors caused by atmospheric delay through spatial differentiation, and combined with spatiotemporal filtering, further separates residual atmospheric components, significantly improving the accuracy of deformation measurement. In areas where satellite signals are obscured, this method can independently provide high-precision displacement measurements, ensuring blind-spot-free displacement monitoring of the entire foundation pit slope surface. Finally, GNSS and radar observation data are fused within the unscented Kalman filtering framework of step two, achieving all-weather, all-directional, and highly reliable displacement monitoring.

[0101] The number of devices and processing scale described herein are for the purpose of simplifying the description of the invention. Applications, modifications, and variations of the invention will be readily apparent to those skilled in the art.

[0102] Although the embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details.

Claims

1. A method for real-time monitoring and early warning of foundation pit slope stability, characterized in that, Includes the following steps: Step 1: Collect data on surface displacement, deep deformation, groundwater level, soil moisture content, rainfall intensity, and construction vibration of the foundation pit slope; Each sensor appends its type, accuracy level, and GNSS timestamp to the header of its data message and connects to the edge computing gateway through a unified interface; Step 2: The edge computing gateway aligns the timestamps with GNSS integer second pulses, transforms all coordinates to the BIM coordinate system, and fuses multiple physical quantities using unscented Kalman filtering. The data is constructed into a third-order tensor according to physical quantities, time windows, and spatial measurement points. The core tensor and factor matrix are obtained through Tucker decomposition. The coupling strength matrix is ​​extracted from the core tensor. Spatiotemporal sub-blocks dominated by displacement and simultaneously satisfying the coupling strength conditions with groundwater level and soil moisture content are selected. The corresponding original signals are extracted, and transient impact sub-blocks dominated by construction vibration are removed. A standardized dataset is formed and uploaded to the cloud. Step 3: Track the coupling strength matrix in the cloud. When the coupling strength between displacement and other physical quantities increases and exceeds the threshold in multiple consecutive time windows, the mechanism-data dual-driven hybrid early warning model is triggered to switch from low-power inference mode to high-precision early warning mode. The extended finite element seepage-stress coupling equation is used as a physical constraint term, and LSTM loss function is embedded for joint training. A physical residual loss term is added in addition to the mean square error loss. At the same time, the measured instability cases are divided into multiple types according to the evolution speed. Between adjacent instability stages, the conditional Wasserstein generative adversarial network is used to learn the Riemannian manifold and perform geodesic interpolation to generate virtual samples to expand the training set. The model outputs the slope instability probability and early warning level; Step 4: Dynamically map the warning level and instability probability to the BIM model on the digital twin platform, and automatically trigger a closed-loop response process of verification, handling and feedback based on the warning level.

2. The method for real-time monitoring and early warning of foundation pit slope stability as described in claim 1, characterized in that, Step two, which uses Tucker decomposition to obtain the core tensor and factor matrix, includes: a. Expand the third-order tensor into a two-dimensional matrix along the monitoring physical quantity mode, perform singular value decomposition on the two-dimensional matrix, and take the first few columns of the left singular vector matrix as the monitoring physical quantity factor matrix. b. Expand the same third-order tensor into a two-dimensional matrix along the time window dimension, perform singular value decomposition on the two-dimensional matrix, and take the first few columns of the left singular vector matrix as the time window factor matrix. c. Expand the same third-order tensor into a two-dimensional matrix along the spatial measurement point mode, perform singular value decomposition on the two-dimensional matrix, and take the first few columns of the left singular vector matrix as the spatial measurement point factor matrix. d. Perform a product operation on the first dimension of the original third-order tensor and the transpose of the monitoring physical quantity factor matrix. The operation method is as follows: for each element in the tensor, multiply the index on the monitoring physical quantity dimension with the corresponding value in the transpose of the factor matrix and sum them to obtain a new third-order tensor. e. Perform a product operation in the second dimension on the third-order tensor obtained in step four and the transpose of the time window factor matrix. Multiply the index in the time window dimension with the corresponding value in the transpose of the factor matrix and sum them to obtain another third-order tensor. f. Perform a product operation on the third dimension of the third-order tensor obtained in step 5 and the transpose of the spatial measurement point factor matrix. Multiply the indices on the spatial measurement point dimension with the corresponding values ​​in the transpose of the factor matrix and sum them to obtain the core tensor. g. Extract the absolute value of each element in the core tensor, and arrange these absolute values ​​according to the index positions of the corresponding physical quantities, time windows, and spatial measurement points in the original third-order tensor to form a three-dimensional array with the same structure as the original third-order tensor. This array is the coupling strength matrix.

3. The method for real-time monitoring and early warning of foundation pit slope stability as described in claim 2, characterized in that, Step two involves selecting spatiotemporal sub-blocks that exhibit displacement as the dominant mode and simultaneously satisfy the coupling strength conditions with groundwater level and soil moisture content. This includes: In the monitoring physical quantity factor matrix A, the row vectors corresponding to displacement, groundwater level, and soil moisture content are located. The cosine similarity between the displacement row vector and the groundwater level row vector, and the cosine similarity between the displacement row vector and the soil moisture content row vector are calculated. In the core tensor G, the slice matrix corresponding to the displacement is extracted along the dimension of the monitoring physical quantity. The elements in the slice matrix are thresholded, and the index positions with the absolute value of the element greater than the preset threshold are marked as strong coupling candidate positions. Among the strong coupling candidate positions, the index positions that simultaneously satisfy the displacement-groundwater level cosine similarity exceeding the first threshold and the displacement-soil moisture content cosine similarity exceeding the second threshold are selected to form the target spatiotemporal sub-block. The original signal components corresponding to the third-order tensor at the target spatiotemporal sub-block positions are retained, and the signal components corresponding to the other positions are set to zero to obtain the signal after removing the transient impact sub-block.

4. The method for real-time monitoring and early warning of foundation pit slope stability as described in claim 1, characterized in that, The calculation steps for the physical residual loss term in step three include: On the extended finite element discrete nodes, the displacement field u of the i-th node predicted by the Long Short-Term Memory network is extracted. i and pore water pressure field p i Substituting into the discrete form of the seepage-stress coupling control equation: K u u i +K c p i −f i =R u，i ; K c T (you i / dt)+K p p i −q i =R p，i ; Among them, K u K is the stiffness matrix. c Let K be the coupling matrix. p Let f be the penetration matrix. i Let q be the nodal force vector of the i-th node. i Let R be the node flow vector of the i-th node. u，i Let R be the displacement field residual of the i-th node. p，i Let be the residual of the pore water pressure field at the i-th node; The physical residual loss term is: ； Where N is the total number of discrete nodes, and λ is the weighting coefficient for the difference in dimensions between the residual of the equilibrium displacement field and the residual of the pore water pressure field. The overall loss function for joint training is: ； Where L MSE α is the mean squared error loss term, and α is the physical residual penalty coefficient.

5. The method for real-time monitoring and early warning of foundation pit slope stability as described in claim 1, characterized in that, Step three, which involves using Conditional Wasserstein to generate an adversarial network that learns the Riemannian manifold of the latent space and performing geodesic interpolation on the Riemannian manifold, includes the following steps: The measured samples of various instability evolution paths were divided into initial creep stage, constant creep stage, accelerated creep stage, and critical instability stage according to time series. A conditional Wasserstein generative adversarial network was constructed. The generator input was random noise z and stage condition label c, and the output was the virtual sample of the corresponding stage. The discriminator output was the Wasserstein distance estimate. After the conditional Wasserstein generative adversarial network training converged, a set of random noise z1 and z2 was sampled for adjacent instability stages, and mapped to the latent space by the generator to obtain latent vectors w1 and w2. The geodesic path between the latent vectors was calculated. ； Where exp and log are the exponential and logarithmic mappings on the Riemannian manifold, respectively; Interpolation progress parameters from the starting point to the ending point ( When =0, it corresponds to w1. When =1 corresponds to w2); in Uniform sampling within the value range M Each interpolation point is used to decode the latent vector corresponding to each interpolation point through a generator, resulting in a virtual instability evolution sample that smoothly transitions between adjacent instability stages.

6. The method for real-time monitoring and early warning of foundation pit slope stability as described in claim 1, characterized in that, The switching steps between low-power inference mode and high-precision warning mode in step three include: In low-power inference mode, the mechanism-data dual-driven hybrid early warning model loads a lightweight branch of the Long Short-Term Memory (LSTM) network. This lightweight branch is obtained from the original LSTM network after channel pruning and weight quantization, and the physical residual loss term is not calculated during inference. When the cloud server detects that the displacement-groundwater level coupling strength and the displacement-soil moisture content coupling strength are monotonically increasing within k consecutive time windows, and the coupling strength of the current time window exceeds their respective preset thresholds, the mode switch is triggered. The switching process includes: loading the complete LSTM network weights from the model storage area to replace the lightweight branch weights, and activating the computation node corresponding to the physical residual loss term in the computation graph, so that the model enters the high-precision early warning mode.

7. The method for real-time monitoring and early warning of foundation pit slope stability as described in claim 1, characterized in that, In step one, deep deformation monitoring is achieved by deploying an array-type microelectromechanical inclinometer. The packaging structure of this inclinometer includes: MEMS accelerometer chips are soldered onto a flexible printed circuit board, which is fixed to the inner wall of a glass fiber reinforced epoxy resin composite tube by rigid partitions arranged at equal intervals. The outer wall of the composite tube is formed by spirally winding epoxy resin-impregnated glass fiber cloth along the tube axis at an angle of 45° to 55°. Both ends of the composite tube are sealed with stainless steel end caps and O-rings. The temperature self-compensation circuit includes a thermistor from the same batch as the MEMS accelerometer chip. The thermistor and the MEMS accelerometer chip are arranged adjacent to each other on the same flexible printed circuit board. The temperature self-compensation circuit corrects the accelerometer output in real time through a pre-calibrated temperature-drift curve. The redundant MEMS sensing unit consists of two sets of MEMS accelerometer chips arranged orthogonally on the same rigid partition.

8. The method for real-time monitoring and early warning of foundation pit slope stability as described in claim 1, characterized in that, The implementation steps for the unified data interface in Step 1 include: A predefined data dictionary covering displacement, pressure, moisture content, and vibration parameters is used. Each entry in the data dictionary includes a data code, physical quantity name, unit, data type, and value range. The edge computing gateway allocates independent communication ports and message queues for each type of sensor. Each sensor encapsulates data packets according to the data dictionary format. The header of the data packet contains the sensor type code, accuracy level code, and GNSS timing timestamp, while the body of the packet contains monitoring data in key-value pair format. The data parsing module of the edge computing gateway reads the sensor type code from the packet header, queries the data dictionary to obtain the key-value pair parsing template corresponding to that type of sensor, extracts the monitoring data according to the parsing template, and writes it into the time-series database.

9. The method for real-time monitoring and early warning of foundation pit slope stability as described in claim 8, characterized in that, Step two, which involves using unscented Kalman filtering to perform optimal fusion estimation of multiple physical quantities, includes: Constructing the state vector: ； Where d x d y d z For three-dimensional displacement components, v x v y v z The displacement rate components are three-dimensional; the state transition equation adopts a constant velocity model; in the observation equation, the observation vector of the GNSS receiver is [d]. x d y d z ] T The observation vector of the ground-based synthetic aperture radar is transformed into the BIM coordinate system by the line-of-sight projection matrix as d. x d y d z A linear combination; by querying the preset accuracy level-variance mapping table using the accuracy level labels attached to each sensor message, the observation noise covariance matrix R of each sensor is obtained. k The diagonal element values; calculating the Sigma point set and its weights during the time update phase, and using R during the measurement update phase. k To measure the noise covariance, calculate the Kalman gain and perform state updates.

10. The method for real-time monitoring and early warning of foundation pit slope stability as described in claim 1, characterized in that, In step one, surface displacement monitoring simultaneously deploys a BeiDou / GNSS receiver and a ground-based synthetic aperture radar (SAR). In areas where satellite signals are obscured, the SAR interferometry is used to obtain the line-of-sight deformation sequence. Pixels with amplitude deviation indices less than a preset threshold in the time-series interferogram are selected as candidate permanent scatterers. After time-dimensional phase unwrapping of the candidate points, a differential phase model between adjacent points is established. ； Where Δv is the deformation rate difference, T is the time baseline, B⊥ is the vertical baseline, R is the slant distance, θ is the incident angle, Δh is the elevation error difference, and Δϕ is the vertical distance. atm For atmospheric delay phase difference, Δϕ noise This is the noise phase difference; The atmospheric delay phase component of each permanent scatterer candidate point is separated and estimated by spatiotemporal filtering. The atmospheric delay phase component is then subtracted from the total phase to obtain the corrected deformation phase sequence.

Images