Real-time updating method of digital twin parameter field based on ensemble Kalman filter

Abstract

The present application relates to a kind of digital twin parameter field real-time updating method based on ensemble Kalman filter.It is driven by spatial discrete sampling data to establish initial twin scene with implicit surface interpolation and triangulation, and based on Fisher information gain criterion optimization sensor arrangement to obtain high information content observation signal.Subsequently, ensemble Kalman filter is used in combination with adaptive noise adjustment to carry out robust Bayesian assimilation, solve the problem that static initial model cannot track parameter field dynamic evolution.Further, introduce spatial compatibility weight constraint, mechanical coupling operator joint gate and Laplace smoothing processing, so that the update result considers spatial continuity and mechanical compatibility.Finally, the updated parameter field is used for forward stratum response deduction and risk assessment, and the fusion visualization of stratum model, risk cloud chart and construction suggestion is realized through three-dimensional rendering engine, forming complete digital twin closed loop from real-time perception to prediction decision.

Description

Technical Field

[0001] This invention belongs to the field of digital twin and computer simulation technology, specifically relating to a method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering. Background Technology

[0002] In the field of engineering construction, whether it is shield tunneling, foundation pit excavation, slope treatment, or underground utility tunnel construction, the operation process faces core challenges such as complex and variable geological conditions, difficulty in accurately obtaining stratum parameters, and real-time evolution of construction disturbances. During the operation, the mechanical properties of the strata ahead of the work site often deviate significantly from the spatial discrete sampling data in the exploration phase. Moreover, as engineering equipment continues to advance or excavate, construction disturbances will continuously change the stress state and physical properties of the surrounding strata. If construction decisions still rely on static initial geological exploration models, they cannot reflect the dynamic changes in stratum parameters in a timely manner, which can easily lead to safety accidents such as excessive surface settlement, structural deformation, or even collapse.

[0003] Digital twin technology provides a technical path to solve the above problems. However, existing methods have significant technical shortcomings in the dynamic updating of parameter fields. Specifically, when using posterior parameters obtained from Bayesian inversion to progressively update a 3D stratigraphic model, a local algebraic correction strategy centered on parameter components is usually adopted. That is, the difference between the posterior mean and the previous parameter field is directly used as the evolution increment and written into the parameter field element by element. This approach ignores the inherent spatial correlation of the strata as a continuous medium in engineering operation scenarios and does not consider that the influence of the current step observation information on different spatial locations should be transmitted in a decaying manner along the distance from the operation surface. This can easily lead to geologically unreasonable step-like parameter jumps between adjacent grid nodes, which is more obvious at the interface of soft and hard heterogeneous strata and in groundwater sensitive areas. Meanwhile, mechanical parameters such as elastic modulus, Poisson's ratio, cohesion, and internal friction angle in the formation are coupled with each other, controlled by common geological origins and stress release processes. If updates are performed independently dimension by dimension without applying cross-parameter linkage constraints, a significant jump in one parameter may occur while its cooperating parameters are not adjusted synchronously. This leads to a loss of mechanical coordination in the local parameter combination, resulting in numerical anomalies and prediction instability in subsequent response predictions. Furthermore, residual sampling errors and local noise disturbances in the observation data themselves will be directly incorporated into the parameter field without spatial smoothing constraints, forming local spikes with grid sensitivity. This weakens the twin model's ability to express the gradual characteristics of the real formation and affects the reliability of risk assessment and construction parameter recommendations.

[0004] Therefore, an optimized scheme for real-time updating of digital twin parameters in engineering construction is desired. Summary of the Invention

[0005] To address the aforementioned technical problems, this application is proposed. Embodiments of this application provide a real-time update method for digital twin parameter fields based on ensemble Kalman filtering, comprising: S1, based on spatial discrete sampling data and engineering equipment model data, performs spatial interpolation and triangulation reconstruction on the stratigraphic interface control points, and performs spatial registration and scene integration on the stratigraphic entity mesh and engineering equipment model to obtain the initial twin scene; S2, based on the initial twin scenario and preset operating parameters, performs numerical simulation and sensitivity analysis on the position response of candidate sensors, and optimizes sensor layout, heterogeneous signal acquisition and synchronous preprocessing based on the analysis results to obtain the observation signal vector; S3, based on the observed signal vector, prior parameter set and response surrogate model, performs robust assimilation and Bayesian inversion on the predicted signal and the measured signal, and combines adaptive noise adjustment to perform set recursive update to obtain the posterior parameter set; S4, based on the posterior parameter set, historical operation sequence, previous loop parameter field and engineering equipment pose data, performs joint discrimination and mask update on parameter information gain and historical evolution law, and performs selective correction of parameter field according to the discrimination result to obtain updated parameter field; S5, based on the updated parameter field, engineering equipment pose data and candidate construction parameters, performs evolution prediction and risk assessment of the ground response ahead, and combines the updated parameter field for 3D rendering and parameter feedback to obtain a visualized twin scene and construction parameter suggestions.

[0006] Compared with existing technologies, this invention proposes a real-time update method for digital twin parameter fields based on ensemble Kalman filtering. It establishes an initial 3D twin scene through implicit surface interpolation and triangulation reconstruction driven by spatially discrete sampling data. Based on finite element multi-condition simulation and Fisher's information gain criterion, it achieves optimal spatial arrangement of sensors and synchronous acquisition of heterogeneous signals, thus providing high-information observation input for subsequent parameter inversion. Furthermore, it employs an ensemble Kalman filtering framework combined with an adaptive noise covariance adjustment mechanism to robustly assimilate and Bayesianly invert multi-source observation signals and surrogate model prediction signals, solving the problem that traditional static geological exploration models cannot track the dynamic evolution of stratigraphic parameters. To address the spatial discontinuities and cross-parameter mechanical imbalances that easily occur when writing inverted parameters into the 3D stratigraphic model, it introduces a joint gating correction strategy based on compatibility weight constraints and mechanical coupling operators according to the spatial distance of the working face, and superimposes Laplace spatial smoothing processing, ensuring that the parameter field update simultaneously satisfies the requirements of spatial continuity and mechanical coordination, suppressing spike anomalies driven by local noise. Ultimately, the updated parameter field is used for proxy simulation and risk quantification assessment of the forward stratum response. The 3D rendering engine is used to achieve the fusion visualization of the stratum model, risk cloud map and construction parameter suggestions, forming a complete digital twin closed loop from real-time perception, parameter inversion, scene update to predictive decision-making. Attached Figure Description

[0007] Figure 1 This is a flowchart of a method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering, according to an embodiment of this application. Figure 2 This is a schematic diagram of the data flow of the real-time update method for digital twin parameter fields based on ensemble Kalman filtering according to an embodiment of this application; Figure 3 The flowchart describes a real-time update method for digital twin parameter fields based on ensemble Kalman filtering according to an embodiment of this application. It involves robust assimilation and Bayesian inversion of the predicted signal and the measured signal, combined with adaptive noise adjustment for ensemble recursive update, to obtain the posterior parameter set. Figure 4 To provide a flowchart for the real-time update method of digital twin parameter field based on ensemble Kalman filtering according to the embodiments of this application, the method performs joint discrimination and mask update on parameter information gain and historical evolution law based on posterior parameter set, historical operation sequence, previous loop parameter field and engineering equipment pose data, and performs selective correction of parameter field based on discrimination results, so as to obtain the update parameter field flowchart. Figure 5This is a flowchart illustrating the dynamic update mask-based method for real-time updating of digital twin parameter fields using ensemble Kalman filtering, which involves gating and spatially mapping the posterior mean of the parameters and the parameter field of the previous loop to obtain the updated parameter field. Figure 6 To provide a flowchart for the real-time updating method of digital twin parameter field based on ensemble Kalman filtering according to the embodiments of this application, the method uses the updated parameter field, engineering equipment pose data and candidate construction parameters to predict the evolution and risk assessment of the ground response ahead, and combines the updated parameter field for three-dimensional rendering and parameter feedback to obtain a visualized twin scene and construction parameter suggestions. Detailed Implementation

[0008] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0009] Figure 1 This is a flowchart of a method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering, according to an embodiment of this application. Figure 2 This is a schematic diagram of the data flow in the real-time update method for digital twin parameter fields based on ensemble Kalman filtering according to an embodiment of this application. Figure 1 and Figure 2As shown, the real-time update method for digital twin parameter fields based on ensemble Kalman filtering according to an embodiment of this application includes the following steps: S1, based on spatial discrete sampling data and engineering equipment model data, spatial interpolation and triangulation reconstruction are performed on the control points of the stratigraphic interface, and spatial registration and scene integration are performed on the stratigraphic entity mesh and the engineering equipment model to obtain an initial twin scene; S2, based on the initial twin scene and preset operating parameters, numerical simulation and sensitivity analysis are performed on the position response of candidate sensors, and sensor optimization, heterogeneous signal acquisition and synchronous preprocessing are performed based on the analysis results to obtain the observation signal vector; S3, based on the observation signal vector, prior parameter set and response data... The model robustly assimilates and performs Bayesian inversion on the predicted and measured signals, and combines adaptive noise adjustment for ensemble recursive updates to obtain the posterior parameter set; S4, based on the posterior parameter set, historical operation sequences, previous loop parameter field, and engineering equipment pose data, jointly discriminates and updates the parameter information gain and historical evolution patterns, and selectively corrects the parameter field based on the discrimination results to obtain the updated parameter field; S5, based on the updated parameter field, engineering equipment pose data, and candidate construction parameters, performs evolution prediction and risk assessment of the forward stratum response, and combines the updated parameter field for 3D rendering and parameter feedback to obtain a visualized twin scene and construction parameter suggestions.

[0010] Specifically, in step S1, based on spatially discrete sampling data and engineering equipment model data, spatial interpolation and triangulation reconstruction are performed on the stratigraphic interface control points. Spatial registration and scene integration are then performed between the stratigraphic entity mesh and the engineering equipment model to obtain an initial twin scene. It should be noted that, since only discretely distributed borehole exploration data can be obtained before construction (e.g., in shield tunnel engineering, only geological exploration information from a limited number of boreholes along the route is available before construction), the stratigraphic boundary information between boreholes is spatially discontinuous and cannot directly form a three-dimensional continuous stratigraphic model. Furthermore, the engineering equipment model and the stratigraphic data belong to different coordinate systems, lacking a unified spatial reference. Therefore, the technical solution of this application first performs spatial interpolation and triangulation reconstruction on the stratigraphic interface control points based on spatially discrete sampling data and engineering equipment model data, and then performs spatial registration and scene integration between the stratigraphic entity mesh and the engineering equipment model to obtain an initial twin scene. Through the above processing, discrete borehole information can be transformed into a topologically complete three-dimensional geological entity and integrated with the engineering equipment model in the same coordinate system, providing an operable digital twin base for subsequent sensor deployment simulation and dynamic parameter updates.

[0011] More specifically, in a specific example of this application, the spatially discrete sampling data is first analyzed to extract the three-dimensional spatial coordinates of each borehole and its corresponding stratigraphic elevation information. The top and bottom coordinates at the interfaces of each stratigraphic unit are used as the original interface control points. Since the interface control points alone are insufficient to uniquely determine the normal trend of the surface, they are further offset upwards and downwards by a preset distance along the normal direction of the local fitting plane where each control point is located, generating two sets of auxiliary constraint points. The original control points are assigned a function value of zero, the upwardly offset auxiliary constraint points are assigned a function value of positive one, and the downwardly offset auxiliary constraint points are assigned a function value of negative one. These three types of points are aggregated to form a control point dataset with spatial directional constraints. Based on this, implicit surface spatial interpolation is performed on the control point dataset using radial basis functions. Specifically, with the coordinates of each node and the corresponding function values ​​in the control point dataset known, the following implicit surface interpolation equation is established: in, Represents any evaluation point in three-dimensional space The implicit function scalar value at a point indicates that the point lies on a stratigraphic boundary surface when the value is zero. This represents the total number of nodes in the control point dataset. Indicates the first The radial basis function weight coefficients corresponding to each control point This represents the radial basis function kernel, used to describe how the influence between two points varies with distance. Indicates assessment points With the Control points The Euclidean distance between them The additional low-order polynomial drift term is used to ensure the global translation and rotation invariance of the interpolation. After obtaining all weight coefficients by solving the above equations, a uniform 3D prediction mesh is generated within the spatial bounding box of the target work area. The implicit function scalar value is calculated node by node, and discrete points on the zero isosurface are extracted to obtain the 3D prediction point cloud of each stratum interface. Subsequently, Delaunay triangulation is performed on the 3D prediction point cloud to generate a closed triangular patch mesh that satisfies the empty circle property. Degenerate triangles with aspect ratios exceeding the threshold are removed to ensure mesh quality. Based on the lithology of each stratum, the closed entities are assigned corresponding material rendering labels to form a topologically complete stratum entity mesh. Finally, a homogeneous coordinate transformation matrix is ​​constructed based on the engineering coordinates and initial attitude angles of the engineering equipment's starting position (e.g., the coordinates of the shield machine's starting shaft and initial attitude angles). 3D spatial rotation and translation transformations are performed on the engineering equipment model data to achieve spatial registration with the stratum entity mesh in the same global geographic coordinate system. The two are then integrated into the 3D rendering engine and encapsulated into an interactive scene file to obtain the initial twin scene.

[0012] Specifically, in step S2, based on the initial twin scenario and preset operating parameters, numerical simulation and sensitivity analysis are performed on the candidate sensor position responses. Based on the analysis results, sensor layout optimization, heterogeneous signal acquisition, and synchronous preprocessing are then performed to obtain the observation signal vector. It should be noted that, due to the dual constraints of physical space and economic cost on the number and installation location of sensors during engineering operations, if the sensor layout lacks quantitative optimization basis, the collected observation data may not adequately constrain key stratigraphic parameters, affecting the accuracy and convergence of subsequent Bayesian inversion. Therefore, the technical solution of this application further performs numerical simulation and sensitivity analysis on the candidate sensor position responses based on the initial twin scenario and preset operating parameters, and performs sensor layout optimization, heterogeneous signal acquisition, and synchronous preprocessing based on the analysis results to obtain the observation signal vector. Through the above processing, the information gain of the observation data on the target inversion parameters can be maximized under limited sensor resources, and multi-source heterogeneous signals can be unified into a standardized vector that can be directly input into the assimilation algorithm.

[0013] More specifically, in a specific example of this application, a finite element coupled model is first established based on the geometric topological relationship between the 3D geological entity mesh and the engineering equipment model in the initial twin scene. After setting boundary conditions, multiple combinations of geotechnical mechanics parameters and construction parameters included in the preset working condition parameters are imported to drive multi-working condition forward numerical simulation. After each simulation calculation, physical quantities such as stress, strain, displacement, and pore water pressure at all candidate sensor virtual nodes are extracted. The extraction results of each working condition are aggregated and stored to form a simulation response database. After obtaining the simulation response database, sensitivity analysis and information gain evaluation are further performed on the sensor responses at each candidate location. Specifically, the partial derivatives of the sensor responses at each candidate location with respect to the target inverted geological parameters are calculated using the finite difference method, and the Fisher information matrix is ​​constructed by combining the preset basic noise variance of each candidate sensor. The calculation relationship is as follows: in, Represents the Fisher information matrix. This represents the total number of candidate sensor nodes. Indicates the first The fundamental noise variance of each candidate position sensor Represents the first in the simulation response database Predicted physical response values ​​for each candidate location under specific operating conditions. This represents the column vector of target formation parameters to be inverted, including mechanical parameters such as cohesion, internal friction angle, and elastic modulus. Based on this, using the D-optimal criterion as the objective function, i.e., maximizing the determinant of the Fisher information matrix, a greedy algorithm iteratively selects the sensor locations that bring the maximum information increment from all candidate nodes until the preset upper limit of the number of sensors is reached, thereby obtaining the optimal sensor layout scheme.

[0014] Based on the optimal sensor layout scheme, physical monitoring devices such as earth pressure gauges, strain gauges, and settlement gauges are deployed at corresponding locations on the engineering equipment body, pre-formed structural components, and the ground surface (e.g., on the tunnel boring machine body, assembled segments, and corresponding locations on the ground surface). A distributed data acquisition link is established, and synchronous online acquisition is performed on each channel during the operation of the engineering equipment. An absolute timestamp and the current working mileage are added to each data stream, forming original heterogeneous signals with different sampling frequencies and physical dimensions. Subsequently, Kalman low-pass filtering for noise reduction and extreme value anomaly removal are performed on the original heterogeneous signals. To address the issue of inconsistent sampling frequencies among the sensors, cubic spline interpolation is used to time-align and upsample the low-frequency channel data, with the global master time axis as the reference. The physical quantity data of all channels under the same synchronous timestamp are spliced ​​and fused in a predetermined dimensional order to finally obtain the observation signal vector.

[0015] Specifically, in step S3, based on the observed signal vector, prior parameter set, and response surrogate model, robust assimilation and Bayesian inversion are performed on the predicted signal and the measured signal, and adaptive noise adjustment is combined to perform set recursive updates to obtain the posterior parameter set. It should be noted that, given that the multi-source observation signals collected during engineering operations inevitably contain abnormal readings introduced by sensor drift, electromagnetic interference, and local faults, if the observation data and the surrogate model prediction results are directly assimilated without identifying and suppressing abnormal channels, noise disturbances will be transmitted to the parameter update results via the Kalman gain, causing the inverted formation parameters to deviate from the true values. Furthermore, finite set samples are prone to generating spurious long-range correlations in high-dimensional parameter spaces, further reducing the stability of the inversion. Therefore, the technical solution of this application further performs robust assimilation and Bayesian inversion on the predicted signal and the measured signal based on the observed signal vector, prior parameter set, and response surrogate model, and combines adaptive noise adjustment to perform set recursive updates to obtain the posterior parameter set. Through the above processing, the interference of abnormal sensor channels on the inversion results can be automatically identified and suppressed in each ring assimilation calculation. At the same time, the spurious correlation introduced by finite samples is suppressed by localization truncation and covariance expansion, so that the obtained posterior parameter set is closer to the real formation state in a statistical sense.

[0016] Figure 3This document describes a flowchart illustrating a real-time update method for digital twin parameter fields based on ensemble Kalman filtering, according to embodiments of this application. The method involves robust assimilation and Bayesian inversion of the predicted and measured signals using the observed signal vector, prior parameter set, and response surrogate model, combined with adaptive noise adjustment for ensemble recursive updating, to obtain the posterior parameter set. (See flowchart for details.) Figure 3 As shown, step S3 includes: S31, based on the response surrogate model, performing forward prediction and signal mapping on the prior parameter set to obtain the predicted signal set; S32, based on the observed signal vector and the predicted signal set, performing anomaly detection and noise adjustment on the sensor innovation statistics to obtain the adaptive observation noise covariance matrix; S33, through covariance statistics and localization constraints, solving for the gain of the prior parameter set, the predicted signal set, and the adaptive observation noise covariance matrix to obtain the Kalman gain matrix; S34, based on the Kalman gain matrix, performing recursive update and dilation correction on the observed signal vector, the prior parameter set, and the predicted signal set to obtain the posterior parameter set.

[0017] In step S31, based on the response surrogate model, forward prediction and signal mapping are performed on the prior parameter set to obtain the predicted signal set. It should be noted that since the ensemble Kalman filter framework requires obtaining the sensor prediction readings corresponding to each parameter sample member before each assimilation calculation, directly calling the finite element model for forward calculation is far more time-consuming than the update cycle requirements of real-time operations. Therefore, the technical solution of this application further uses the response surrogate model to perform forward prediction and signal mapping on the prior parameter set to obtain the predicted signal set. Through the above processing, millisecond-level computational cost can replace finite element forward simulation, providing a prediction benchmark of the same dimension as the observed signal vector for subsequent assimilation stages.

[0018] More specifically, in a concrete example of this application, each Monte Carlo parameter sample member in the prior parameter set is extracted. Each member contains the formation mechanical parameters to be inverted in the current loop, such as elastic modulus, Poisson's ratio, cohesion, and internal friction angle. Each parameter member is sequentially input into a deep neural network response surrogate model pre-trained based on finite element simulation data. This surrogate model takes the formation mechanical parameters as input and outputs physical quantities such as stress, displacement, and pore water pressure at the corresponding sensor placement locations. Forward mapping calculations are performed on each parameter member to obtain the predicted response value of that member at each sensor virtual node under the current operating conditions. The prediction results of all members are aggregated and arranged according to the set dimension to form a predicted signal set with the same physical dimensions and channel dimensions as the observed signal vector.

[0019] Specifically, the aforementioned response proxy model requires an offline construction and training phase before being deployed for online forward prediction. Training data originates from the finite element coupled model established in step two. Multiple sets of formation mechanical parameters from the preset working condition parameters are used as input samples, and simulation results of stress, displacement, and pore water pressure at the corresponding sensor locations are used as supervision labels, forming a set of input-output sample pairs covering the target parameter space. The network structure employs a multi-layer fully connected deep neural network. The input layer dimension matches the number of formation parameters to be inverted, and the output layer dimension matches the total number of channels at the sensor locations. Multiple hidden layers are set in the middle, with batch normalization and nonlinear activation functions configured layer by layer to enhance the network's ability to fit the nonlinear mapping relationship between formation parameters and sensor responses. A random dropout mechanism is introduced between hidden layers to prevent overfitting. During training, the mean squared error loss function is used to measure the deviation between the network's predicted response and the finite element simulation response. An adaptive moment estimation optimizer is used for gradient descent iterations, and the loss value is monitored on an independently partitioned validation set. Training is terminated when the validation set loss no longer decreases for several consecutive rounds due to an early stopping mechanism. After training, the network weights are solidified and encapsulated. During the online operation phase of the engineering task, the finite element forward simulation is replaced with millisecond-level computational cost, and fast forward mapping is performed on each member in the prior parameter set.

[0020] In step S32, based on the observed signal vector and the predicted signal set, anomaly detection and noise adjustment are performed on the sensor innovation statistics to obtain an adaptive observation noise covariance matrix. It should be noted that, due to factors such as groundwater erosion, mechanical vibration, and electromagnetic interference affecting sensors at engineering sites, some channels may experience reading drift or abrupt changes during specific periods. If a fixed noise weight is assigned to all sensors in the assimilation calculation, the deviation signal of the abnormal channel will directly contaminate the parameter inversion result via the Kalman gain. Therefore, the technical solution of this application further performs anomaly detection and noise adjustment on the sensor innovation statistics based on the observed signal vector and the predicted signal set to obtain an adaptive observation noise covariance matrix. Through the above processing, abnormal sensor channels can be automatically identified before each assimilation calculation, and their weight contribution in the gain solution can be reduced, making the assimilation process robust to local sensor failures.

[0021] More specifically, in a concrete example of this application, firstly, based on the factory calibration parameters of each sensor, a basic noise variance is set for each sensor channel. The basic noise variances of all channels are then arranged diagonally to construct an initial observation noise covariance matrix. Based on this, the arithmetic mean of the predicted signal set is calculated along the sample dimension to obtain the predicted mean of each channel. The difference between the measured value of each channel in the observed signal vector and the corresponding predicted mean is calculated, and normalized using the basic noise variance of that channel. The normalized squared innovation statistic for each sensor channel is then calculated. The normalized squared innovation statistic of each channel is compared one by one with a 99% confidence level and a chi-square distribution threshold of one degree of freedom. If the statistic of a channel exceeds this threshold, it is marked as an anomalous sensor. For sensor channels marked as anomalous, their basic noise variance is multiplied by a preset expansion coefficient for amplification correction. Unmarked channels retain their basic noise variance unchanged. This yields an adaptive observation noise covariance matrix, whose diagonal elements are selected according to the following rules: in, The first element on the diagonal of the adaptive observation noise covariance matrix represents the... Correction noise variance for each sensor channel This represents the preset variance inflation coefficient. Indicates the first The fundamental noise variance of each sensor Represents the first element in the observed signal vector. Measured values ​​for each channel Represents the first in the set of predicted signals The arithmetic mean of each channel over all sample members Indicates the confidence level as Furthermore, a chi-square distribution test threshold with one degree of freedom is used. After the above processing, the noise variance corresponding to the abnormal channel is amplified, and its weight contribution in the subsequent Kalman gain solution is reduced accordingly, thereby achieving adaptive suppression of abnormal sensor readings.

[0022] In step S33, the Kalman gain matrix is ​​obtained by solving the gain of the prior parameter set, the predicted signal set, and the adaptive observation noise covariance matrix through covariance statistics and localization constraints. It should be noted that, because the covariance matrix estimated by ensemble Kalman filtering under finite sample conditions contains sampling noise, it is prone to generating spurious statistical correlations between spatially uncorrelated parameters and observation channels. If used directly for gain calculation without constraints, parameters far from the sensor location will also be unreasonably corrected. Therefore, the technical solution of this application further uses covariance statistics and localization constraints to solve the gain of the prior parameter set, the predicted signal set, and the adaptive observation noise covariance matrix to obtain the Kalman gain matrix. Through the above processing, the correction weights of each observation channel for each parameter to be inverted can be accurately quantified while suppressing spurious long-range correlations.

[0023] More specifically, in a concrete example of this application, covariance estimation is first performed based on the sample statistics of the prior parameter set and the prediction signal set. The parameter bias matrix is ​​obtained by subtracting each member of the prior parameter set from its set mean, and the prediction bias matrix is ​​obtained by subtracting each member of the prediction signal set from its set mean. A normalized outer product operation is performed on the parameter bias matrix and the prediction bias matrix to obtain the parameter-observation cross-covariance matrix. A normalized outer product operation is then performed on the prediction bias matrix itself to obtain the observation-observation autocovariance matrix. Based on this, to eliminate spurious long-range correlations introduced by a finite set of samples, a localized truncation matrix constructed based on the Gaspari-Cohn distance decay function is introduced. This truncation matrix is ​​then multiplied by the parameter-observation cross-covariance matrix using the Hadamard product operation, causing the covariance values ​​between parameters and observation channels with large spatial distances to decay close to zero. The Kalman gain matrix is ​​obtained by multiplying the locally corrected cross-covariance matrix by the inverse of the sum of the observation-observation autocovariance matrix and the adaptive observation noise covariance matrix. The calculation relationship is as follows: in, Represents the Kalman gain matrix. This represents the localized truncation matrix generated based on the Gaspari-Cohn function. The Hadamard product operation represents the element-wise multiplication of matrices. This represents the parameter-observation cross-covariance matrix. This represents the observation-observation autocovariance matrix. This represents the adaptive observation noise covariance matrix. The Kalman gain matrix characterizes the proportional distribution of information from each observation channel to each formation parameter to be inverted during the current ring assimilation calculation.

[0024] In step S34, based on the Kalman gain matrix, the observed signal vector, prior parameter set, and predicted signal set are recursively updated and expanded to obtain the posterior parameter set. It should be noted that during the recursive update process of ensemble Kalman filtering, nonlinear mapping and finite sample effects cause the set members to cluster, leading to a continuous underestimation of the posterior covariance. Without correction, this will gradually result in a loss of responsiveness to new observation data in subsequent cycles. Therefore, the technical solution of this application further uses the Kalman gain matrix to recursively update and expand the observed signal vector, prior parameter set, and predicted signal set to obtain the posterior parameter set. Through the above processing, reasonable dispersion of the set can be maintained while completing the parameter state update, ensuring the filter remains effective in continuous multi-cycle operations.

[0025] More specifically, in a concrete example of this application, a random perturbation vector following a zero-mean Gaussian distribution is first generated for each set member. This perturbation vector is then superimposed on the observed signal vector to form the pseudo-observation vector corresponding to that member, in order to maintain the statistical randomness of the observation. Subsequently, the pseudo-observation vector of each member is differentially analyzed channel-by-channel with its corresponding predicted reading in the predicted signal set to obtain a residual vector reflecting the deviation between the model prediction and the actual observation. Using the Kalman gain matrix, this residual vector is mapped from the observation space back to the formation parameter space to generate the parameter correction increment for each member. This increment is then accumulated onto the corresponding member in the prior parameter set, completing a single recursive update. The calculation relationship is as follows: in, Represents the th element in the posterior parameter set. The updated parameter vector of each member. Represents the first in the prior parameter set A parameter vector of 1 member, Represents the Kalman gain matrix. Represents the observed signal vector. Represented as the first Zero-mean Gaussian random perturbation vectors generated by each member Represents the first in the set of predicted signals The predicted signal vectors of each member. After the recursive update of all members is completed, the updated set is subjected to covariance scalar inflation, that is, each member is uniformly expanded outward with the set mean as the center by a preset inflation factor to compensate for the covariance underestimation effect accumulated during the recursion process, and finally the posterior parameter set is obtained.

[0026] Specifically, in step S4, based on the posterior parameter set, historical operation sequence, previous loop parameter field, and engineering equipment pose data, joint discrimination and mask update are performed on parameter information gain and historical evolution law. Selective correction of the parameter field is then performed based on the discrimination results to obtain the updated parameter field. It should be noted that, given that the posterior parameter set obtained from Bayesian inversion includes residual observation noise and the influence of finite sample bias, if the posterior mean is directly written into the 3D stratigraphic model in its entirety, noise-driven spurious parameter fluctuations will also enter the parameter field. Furthermore, the inversion reliability of different parameter components varies, and some parameters have limited information gain under the constraints of the current loop observation data. Blindly updating these parameters may introduce additional disturbances. In addition, in engineering operation scenarios, the stratigraphy is a continuous medium, and parameter field updates must also meet the requirements of spatial continuity and coupling coordination between mechanical parameters. If element-wise algebraic correction is performed only around parameter components without considering the spatial proximity relationship of the work surface and cross-parameter linkage constraints, step-like jumps can easily occur between adjacent grid nodes, leading to a loss of mechanical consistency in local parameter combinations. Based on this, the technical solution of this application further utilizes the posterior parameter set, historical operation sequence, previous parameter field, and engineering equipment pose data to jointly discriminate and mask the parameter information gain and historical evolution law. Based on the discrimination results, the parameter field is selectively corrected to obtain an updated parameter field. Through the above processing, low-confidence parameter fluctuations can be screened out based on both statistical evidence and historical patterns. Furthermore, through spatial compatibility weight constraints, joint gating of mechanical coupling operators, and Laplace spatial smoothing, the parameter field update simultaneously meets the requirements of local correction sensitivity, spatial distribution continuity, and mechanical parameter coordination.

[0027] Figure 4 To illustrate the real-time update method for digital twin parameter fields based on ensemble Kalman filtering according to embodiments of this application, a flowchart is provided for jointly discriminating and masking parameter information gain and historical evolution patterns based on posterior parameter sets, historical operation sequences, previous loop parameter fields, and engineering equipment pose data. The flowchart also describes selective parameter field correction based on the discriminating results. For example... Figure 4 As shown, step S4 includes: S41, based on statistical analysis and standard deviation comparison, performing mean extraction and information gain calculation on the posterior parameter set to obtain the posterior mean of the parameters and the Bayesian uncertainty reduction rate vector; S42, based on historical job sequences and long short-term memory networks, performing time-series modeling and pattern extraction on the posterior mean of the parameters and the Bayesian uncertainty reduction rate vector to obtain the parameter update priority vector; S43, performing mask fusion on the Bayesian uncertainty reduction rate vector and the parameter update priority vector to obtain the dynamic update mask; S44, based on the dynamic update mask, performing gating correction and spatial mapping on the posterior mean of the parameters and the parameter field of the previous loop to obtain the updated parameter field.

[0028] In step S41, based on statistical analysis and standard deviation comparison, the mean of the posterior parameter set is extracted and information gain is calculated to obtain the posterior mean of the parameters and the Bayesian uncertainty reduction rate vector. It should be noted that since the posterior parameter set exists in the form of a Monte Carlo sample cluster, it needs to be transformed into deterministic point estimates before it can be used for parameter field correction. Furthermore, the information gain of each parameter component differs under the constraints of the current ring observation data, and this needs to be quantified to distinguish which parameters have indeed received effective updates. Based on this, the technical solution of this application further extracts the mean of the posterior parameter set and calculates information gain based on statistical analysis and standard deviation comparison to obtain the posterior mean of the parameters and the Bayesian uncertainty reduction rate vector. Through the above processing, the optimal estimates of each parameter and the quantification of the degree of constraint by the observed data can be obtained simultaneously, providing a statistical evidence dimension for subsequent mask construction.

[0029] More specifically, in a concrete example of this application, the arithmetic mean of each parameter in the posterior parameter set is first calculated along the set dimension to obtain the posterior mean of the parameters, which includes the optimal estimates of various stratum mechanical parameters such as elastic modulus, Poisson's ratio, cohesion, and internal friction angle. Based on this, the standard deviation of each parameter component in the posterior parameter set along the set dimension is extracted as the posterior standard deviation of this loop, and the prior standard deviation saved before the inversion update is retrieved. The two are then compared element-wise to calculate the Bayesian uncertainty reduction rate vector, the calculation relationship of which is as follows: in, This represents the Bayesian uncertainty reduction rate vector. Represents a vector that is all one. This represents the posterior standard deviation vector of the current loop, obtained statistically from the set of posterior parameters. This represents the prior standard deviation vector before inversion. When the values ​​of each element of this vector approach one, it indicates that the uncertainty of the corresponding parameter under the constraints of the loop observation data has been sufficiently compressed, i.e., the information gain is high. When the values ​​approach zero or are negative, it indicates that the corresponding parameter has hardly obtained effective constraints from the loop observations.

[0030] In step S42, based on historical operation sequences and a long short-term memory network, temporal modeling and pattern extraction are performed on the posterior mean of the parameters and the Bayesian uncertainty reduction rate vector to obtain a parameter update priority vector. It should be noted that parameter update judgment based solely on the Bayesian uncertainty reduction rate of the current single loop is easily affected by fluctuations in a single observation. However, the evolution of formation parameters during engineering operations exhibits temporal continuity and regularity along the operation direction, and the parameter change trends of historical multiple loops can provide empirical support for the update decision of the current loop. Based on this, the technical solution of this application further performs temporal modeling and pattern extraction on the posterior mean of the parameters and the Bayesian uncertainty reduction rate vector based on historical operation sequences and a long short-term memory network to obtain a parameter update priority vector. Through the above processing, the temporal patterns of parameter evolution can be extracted from historical operation data, providing a historical empirical dimension for subsequent mask construction.

[0031] More specifically, in a concrete example of this application, construction parameter records and sensor reading records from multiple past cycles in the historical operation sequence are first extracted. These are then concatenated along the time dimension with the posterior mean of the parameters calculated for the current cycle and the Bayesian uncertainty reduction rate vector, constructing a multi-dimensional temporal input tensor covering the historical window up to the current moment. This temporal input tensor is then fed into a pre-trained two-layer long short-term memory network for recursive processing. The forget gate, input gate, and output gate work together to capture features and pass hidden states regarding the evolutionary trends of parameters across cycles, extracting temporal features reflecting whether each parameter is in a continuously changing phase during continuous operations. Finally, the hidden state vector output by the long short-term memory network is mapped to an output space with the same dimension as the parameters through a fully connected layer, generating update priority weights for each layer's parameters based on historical evolution patterns, resulting in a parameter update priority vector.

[0032] Specifically, the aforementioned Long Short-Term Memory (LSTM) network requires an offline training phase before being deployed for online inference. Training data is derived from historical engineering operation records (e.g., historical tunnel boring machine (TBM) construction records) under similar geological conditions or in the same geological region. Sequences of construction parameters, sensor readings, and ground-value sequences of geological parameters for each ring obtained through offline inversion are extracted from these data. Fixed-length time-series segments are extracted using a sliding window approach as input samples, and the label vectors indicating whether parameters at the end of the corresponding ring have undergone effective evolution serve as supervision signals. The network structure employs a stacked two-layer LSM unit architecture. The dimension of the hidden state in each layer is set according to the number of parameters to be inverted and the temporal complexity. A random dropout mechanism is introduced between the two layers to suppress overfitting. A fully connected layer at the end maps the hidden state to the same output space as the parameter dimension, and the output value is constrained to between 0 and 1 using a Sigmoid activation function. The training process uses a binary cross-entropy loss function to measure the deviation between the predicted priority weights and the true evolution labels, and an adaptive moment estimation optimizer performs gradient descent iterations until the loss value on the validation set converges. After training, the network weights are fixed. During the online operation phase of the engineering task, the network is directly called to perform forward inference on the temporal input tensor of the current loop and update the priority vector of the output parameters.

[0033] In step S43, the Bayesian uncertainty reduction rate vector and the parameter update priority vector are masked and fused to obtain a dynamic update mask. It should be noted that the Bayesian uncertainty reduction rate vector reflects the strength of the statistical constraints imposed on the parameters by the current ring observation data, while the parameter update priority vector reflects the empirical support of historical operational patterns for parameter evolution. Since they determine whether a parameter should be updated from different dimensions, making a decision based on only one dimension carries the risk of misjudgment. Therefore, the technical solution of this application further performs masking and fusion of the Bayesian uncertainty reduction rate vector and the parameter update priority vector to obtain a dynamic update mask. Through the above processing, statistical evidence and historical patterns can be cross-validated, and updates are only allowed when a parameter simultaneously meets the update conditions of both dimensions, thereby filtering out noise-driven one-sided spurious fluctuations.

[0034] More specifically, in a concrete example of this application, the Bayesian uncertainty reduction rate vector and the parameter update priority vector are multiplied by a Hadamard product, that is, corresponding elements of the two vectors of the same dimension are multiplied one by one to generate a dynamic update mask. The calculation relationship is as follows: in, This indicates that the mask is updated dynamically. This represents the Bayesian uncertainty reduction rate vector. This represents the Hadamard product operation. This represents the parameter update priority vector. Due to the multiplicative nature of the Hadamard product, only when a parameter component takes a high value in the Bayesian uncertainty reduction rate vector (i.e., the current loop observation data sufficiently constrains the parameter) and also takes a high value in the parameter update priority vector (i.e., historical operational patterns support the parameter's evolution in the current stage), will its corresponding element in the dynamic update mask approach one. This indicates that the parameter is allowed to be updated at near full amplitude. Conversely, if any dimension takes a low value, the product approaches zero, and the update of the parameter will be suppressed.

[0035] In step S44, based on the dynamically updated mask, the posterior mean of the parameters and the parameter field of the previous loop are gated, corrected, and spatially mapped to obtain the updated parameter field. It should be noted that, given that the dynamically updated mask only expresses the statistical confidence and historical support of whether each parameter component is worth updating in the current loop, and has not yet incorporated the relationship between parameter updates and spatial continuity, as well as the relationship between parameter component updates and mechanical coupling consistency into the correction link, if the difference between the posterior mean of the parameters and the parameter field of the previous loop is directly gated element-wise using the mask and superimposed as is, non-geologically reasonable step jumps can easily form between adjacent grid nodes of the strata before and after the operation. Furthermore, mechanical parameters such as elastic modulus, Poisson's ratio, cohesion, and internal friction angle, which are controlled by common geological origins, may experience imbalances where one parameter jumps while the coordinating parameters are not adjusted synchronously. In addition, local noise disturbances remaining in the observation data will be directly written into the parameter field to form local spikes when there is a lack of spatial smoothing constraints. Based on this, the technical solution of this application further uses a dynamically updated mask to perform gating correction and spatial mapping on the posterior mean of the parameters and the parameter field of the previous loop to obtain the updated parameter field. Through the above processing, the correction effect of the current step observation information can be attenuated and distributed according to the distance from the working surface through spatial compatibility weight constraints. Through the mechanical coupling operator, each parameter can satisfy the preset linkage ratio relationship during the update. Furthermore, through Laplace spatial smoothing, the spikes and inflections caused by mask boundary differences and local noise are weakened, so that the updated parameter field spatially approximates the continuous variation characteristics of the real strata and maintains the coordination of the mechanical parameter combination.

[0036] Figure 5 This document describes a flowchart illustrating a method for real-time updating digital twin parameter fields based on ensemble Kalman filtering, according to embodiments of this application. It describes a process involving gating and spatially mapping the posterior mean of the parameters and the parameter field from the previous loop to obtain the updated parameter field. For example... Figure 5As shown, step S44 includes: S441, based on the pose data of the engineering equipment and the mesh coordinates, performing relevant weight estimation on the spatial distance between the updated region and the working surface to obtain a spatial weight matrix; S442, based on the dynamic update mask, the spatial weight matrix, the posterior mean of the parameters, and the parameter field of the previous loop, performing coupling constraints and compatibility correction on the parameter evolution increment to obtain a correction increment; S443, performing continuous fusion and node mapping on the correction increment and the parameter field of the previous loop to obtain an updated parameter field.

[0037] In step S441, based on the pose data of the engineering equipment and the grid coordinates, a spatial weight matrix is ​​obtained by estimating the spatial distance between the update area and the work surface. It should be noted that since the dynamic update mask only expresses the statistical confidence level of whether each parameter component is worth updating, and has not yet incorporated the attenuation relationship of the current ring observation information's influence on different spatial locations into the update chain, if the mask takes a high value, it will indiscriminately activate the update at the entire grid level. Grid nodes far from the work surface will also be subject to the same magnitude of correction as nodes near the work surface, resulting in a hard-switching characteristic at the update boundary. Based on this, the technical solution of this application further estimates the spatial distance between the update area and the work surface based on the pose data of the engineering equipment and the grid coordinates to obtain a spatial weight matrix. Through the above processing, the engineering common sense that the current ring information can more reliably act on the area near the work surface, while only retaining limited corrections for areas far from the work surface, can be encoded into the update process, transforming the update boundary from a hard-switching to a flexible attenuation.

[0038] More specifically, in a concrete example of this application, the three-dimensional spatial coordinates of the current working face center are first extracted from the pose data of the engineering equipment, and the three-dimensional coordinates of all stratigraphic grid nodes within the region are iterated and updated. The Euclidean distance between the three-dimensional coordinates of each grid node and the coordinates of the working face center is calculated, and a corresponding weight is generated using a Gaussian spatial correlation kernel function. The closer the distance, the higher the weight; the farther the distance, the lower the weight, thus obtaining a spatial weight matrix. The calculation relationship is as follows: in, Indicates the first The spatial compatibility weights of each grid node describe the relative strength at which the current ring parameter update should be preserved at that node. Indicates the first The three-dimensional spatial coordinates of each grid node Represents the three-dimensional pose coordinates of the center of the working surface of the engineering equipment. The spatial correlation length controls the rate at which the observation's influence decays spatially. A larger value indicates a wider spread of the update's influence, while a smaller value indicates that the update is more concentrated in the vicinity of the work surface. The spatial compatibility weights of all grid nodes are arranged by node index to form a spatial weight matrix. After this step, the dynamic update mask, which originally only contained the semantics of whether to update, is supplemented with field distribution constraints on the spatial extent to which the update should proceed, transforming the update boundary from a hard switch to a flexible decay.

[0039] Taking the tunnel boring machine (TBM) operating on a certain ring as an example, the coordinates of the center of the working face are known. The grid nodes within five meters of the working face have weights close to one after Gaussian kernel function calculations, indicating that this area is within the direct influence zone of the current ring excavation disturbance, and parameter corrections can be fully preserved. The weights of grid nodes beyond fifteen meters of the working face decay to near zero, indicating that the current ring observation data has limited constraint on this far-field region, and parameter corrections are adaptively suppressed. Meanwhile, the weights of nodes within the five- to fifteen-meter transition zone show a continuous and gradual distribution, avoiding abrupt boundary changes between the corrected and uncorrected areas. The specific value of the spatial correlation length is determined based on the stratum type and the TBM diameter. A larger value is used in soft strata to accommodate the wide range of disturbance propagation, while a smaller value is used in hard rock strata to constrain updates to be concentrated in the vicinity of the working face.

[0040] In step S442, the parameter evolution increment is coupled and compatibility-corrected based on the dynamically updated mask, spatial weight matrix, posterior mean of parameters, and the parameter field of the previous loop to obtain the corrected increment. It should be noted that although the parameter increment after dual gating by the dynamically updated mask and spatial weight matrix simultaneously reflects statistical reliability and spatial propagation constraints, this step alone is insufficient to solve the problem of coordinated updating between different mechanical parameters. In engineering operation scenarios, formation parameters such as elastic modulus, Poisson's ratio, cohesion, and internal friction angle do not evolve completely independently, but are controlled by common geological origins and stress release processes. If gating correction is performed dimension by dimension, a jump in one parameter may occur while the coordinating parameters are not adjusted synchronously. Therefore, the technical solution of this application further couples and corrects the parameter evolution increment based on the dynamically updated mask, spatial weight matrix, posterior mean of parameters, and the parameter field of the previous loop to obtain the corrected increment. Through the above processing, the changes in parameters are no longer isolated, but are subject to the common constraint of the compatibility of mechanical properties within the same formation unit.

[0041] More specifically, in a concrete example of this application, the original parameter evolution increment is first obtained by subtracting the posterior mean of the parameters from the parameter field of the previous loop. This increment reflects the pure change magnitude of the current loop inversion result relative to the parameter state of the previous loop. Then, a dynamic update mask and a spatial weight matrix are applied together to this increment to simultaneously reflect statistical confidence and spatial propagation constraints. Specifically, an element-wise Hadamard product operation is performed on the dynamic update mask, the spatial weight matrix, and the original parameter evolution increment to synchronously suppress low-confidence parameter components and node increments far from the working surface. Based on this, a mechanical coupling operator matrix is ​​further introduced to perform a matrix transformation on the increment after double gating, ensuring that each mechanical parameter satisfies a preset coupling ratio and linkage relationship during updating. The calculation relationship of the increment is then corrected as follows: in, This represents the increment of the compatibility correction parameter, that is, the parameter increment that can be practically adopted under the combined constraints of spatial compatibility and mechanical compatibility. This represents the mechanical coupling operator matrix, used to describe the linkage relationships and proportion constraints between various mechanical parameters. This represents a dynamically updated mask, used to reflect the strength of parameter updates after dual verification by Bayesian evidence and historical patterns. This represents the Hadamard product operation, which multiplies corresponding elements. Represents the spatial compatibility weight matrix. This represents the posterior mean of the parameter. This represents the parameter field of the previous loop. Through this operation, the original parameter evolution increment is first filtered out by a dynamically updated mask to remove low-confidence components, then the spatial compatibility weight matrix suppresses unreasonable diffusion far from the working surface, and finally the cross-parameter collaborative correction is completed through the mechanical coupling operator matrix to obtain a corrected increment that is consistent with current observational evidence and maintains the consistency of formation mechanical logic.

[0042] Taking the tunnel boring machine traversing the boundary between soft and hard strata as an example, the posterior increase in the elastic modulus in the current ring inversion results is quite prominent, while the posterior changes in cohesion and internal friction angle within the same stratum unit are relatively gradual. Without applying mechanical coupling constraints, the elastic modulus would be significantly corrected independently, causing the parameter combination at that node to deviate from the mechanically compatible range of this type of stratum. After matrix transformation using the mechanical coupling operator matrix, the correction magnitude of the elastic modulus is constrained according to a preset linkage ratio, while simultaneously causing a coordinated and moderate adjustment in cohesion and internal friction angle, keeping the updated parameter combination within a mechanically reasonable range. Even if a single parameter has a large increase in posterior statistical significance, as long as its correlation with the surrounding spatial location is not strong, or its coupling relationship with other mechanical parameters does not support such a jump, the magnitude of its final inclusion in the parameter field will be suppressed, thereby avoiding the formation of numerically sharp and physically unbalanced local anomalies.

[0043] In step S443, the correction increment and the previous loop parameter field are continuously fused and mapped to obtain the updated parameter field. It should be noted that after calculating the compatibility correction parameter increment, it is not directly superimposed onto the previous loop parameter field. Given that even after robust filtering and double verification, the observation assimilation results in engineering construction scenarios may still have residual high-frequency disturbances in certain areas, these disturbances can be amplified into local spikes in the formation response if the parameter field is used for subsequent forward response prediction and risk assessment, affecting prediction stability. Therefore, the technical solution of this application further performs continuous fusion and node mapping of the correction increment and the previous loop parameter field to obtain the updated parameter field. Through the above processing, the parameter update not only retains the effective corrections brought by the assimilation results, but also weakens spikes and inflections caused by local noise or mask boundary differences through smoothing terms, making the updated parameter field spatially closer to the continuous variation characteristics of the actual formation.

[0044] More specifically, in a particular example of this application, to mitigate this adverse effect, the Laplace operator is used to spatially smooth the compatibility correction parameter increment, and this smoothed increment is superimposed on the original correction increment into the previous ring parameter field, forming an updated parameter field that balances boundary continuity and update sensitivity. The computational relationship is as follows: in, This represents the final updated parameter field. This represents the parameter field of the previous loop. This represents the smoothing coefficient, used to control the relative weight between the smoothing term and the increment term. This represents the spatial Laplacian operator, used to characterize the second-order trend of parameter increments in space. This represents the compatibility correction parameter increment. In this calculation relationship, the Laplacian operator, after applying the correction increment, extracts its second-order gradient distribution characteristics in the 3D grid space. A diffusion smoothing effect is applied to local peaks and abrupt inflection points in the increment field. The smoothing coefficient controls the strength of this diffusion effect; a larger value results in a stronger smoothing effect, while a smaller value tends to preserve the details of the original correction increment. After simultaneously superimposing the smoothing term and the original correction increment onto the previous loop parameter field, the resulting updated parameter field absorbs the effective parameter corrections driven by the current loop observation data and eliminates discontinuous jumps at the correction boundary through second-order smoothing constraints. The updated parameter field is then mapped to the 3D stratigraphic grid space nodes for subsequent visualization rendering, forward stratigraphic response prediction, and construction parameter optimization decisions, thus completing a stable transition from observation assimilation results to a physically consistent parameter field.

[0045] Taking the continuous operation of a tunnel boring machine traversing the interface between a water-bearing sand layer and a clay layer as an example, at the grid nodes near the stratigraphic interface, the incremental values ​​of the compatibility correction parameters differ significantly between the sand layer side and the clay layer side. If these parameters are directly written into the parameter field without smoothing, the elastic modulus and permeability coefficient of the nodes on both sides of the interface will exhibit a step jump. After processing with the Laplace operator, the gradient of the incremental field near the interface is moderately diffused, causing the parameter values ​​on both sides to gradually transition along the normal direction rather than abruptly truncated. The updated parameter field distribution is consistent with the gradual transition characteristics of the sand-clay interface zone in the actual strata. The smoothing coefficient is calibrated based on the estimated thickness of the stratigraphic gradient zone. A larger value is taken in areas with wider gradient zones to enhance the smoothing effect, while a smaller value is taken in areas with clearer interfaces to preserve the spatial resolution of the inversion results.

[0046] Specifically, in step S5, based on the updated parameter field, engineering equipment pose data, and candidate construction parameters, the evolution prediction and risk assessment of the forward stratum response are performed. This is combined with 3D rendering and parameter feedback using the updated parameter field to obtain a visualized twin scene and construction parameter suggestions. It should be noted that, given that the updated parameter field obtained through the preceding steps has already incorporated the effective constraints of the current environmental observation data and meets the requirements of spatial continuity and mechanical coordination, it needs to be further transformed into a response prediction and quantitative assessment of the unexcavated stratum ahead. Simultaneously, the updated stratum state, the real-time pose of the tunnel boring machine, and risk warning information are visualized in a unified 3D scene to form a complete closed-loop feedback from parameter updates to construction decisions. Based on this, the technical solution of this application further performs evolution prediction and risk assessment of the forward stratum response based on the updated parameter field, engineering equipment pose data, and candidate construction parameters, and combines 3D rendering and parameter feedback using the updated parameter field to obtain a visualized twin scene and construction parameter suggestions. Through the above processing, the risk of settlement and deformation in the working area ahead can be predicted in advance based on the real-time updated geological parameters. The parameter combination with the lowest risk can be selected from multiple candidate construction schemes for construction personnel to refer to. The geological model, risk cloud map and construction suggestion information are spatially integrated and presented through a 3D rendering engine, so that the digital twin scene is updated synchronously after each operation is completed, supporting real-time perception and real-time decision-making in shield tunneling construction.

[0047] Figure 6 To illustrate the real-time updating method for digital twin parameter fields based on ensemble Kalman filtering according to embodiments of this application, this method uses an updated parameter field, engineering equipment pose data, and candidate construction parameters to predict the evolution and risk assessment of the forward geological response. It also combines the updated parameter field with 3D rendering and parameter feedback to obtain a flowchart of a visualized twin scene and construction parameter suggestions. For example... Figure 6As shown, step S5 includes: S51, based on the updated parameter field, engineering equipment pose data, and candidate construction parameters, performing proxy inference and response prediction on the forward stratum state to obtain a forward stratum response prediction set; S52, based on the risk assessment function and threshold constraints, performing risk quantification and scheme selection on the forward stratum response prediction set to obtain construction parameter suggestions and early warning risk cloud map data; S53, assembling and aligning the construction parameter suggestions, updated parameter field, engineering equipment pose data, and early warning risk cloud map data to obtain a fusion data package to be rendered; S54, based on the 3D rendering engine, dynamically rendering and presenting the construction parameter suggestions and the fusion data package to be rendered to obtain a visualized twin scene.

[0048] In step S51, based on the updated parameter field, the pose data of the engineering equipment, and candidate construction parameters, a proxy simulation and response prediction of the underlying strata are performed to obtain a set of predicted responses for the underlying strata. It should be noted that during engineering operations, construction personnel need to obtain predictions of the responses of the unexcavated strata under different combinations of construction parameters before advancing, in order to select the lowest-risk operational plan. The updated parameter field has been corrected for spatial compatibility and mechanical coupling constraints, providing parameter inputs closer to the actual state for predicting the response of the underlying strata. Based on this, the technical solution of this application further performs proxy simulation and response prediction of the underlying strata based on the updated parameter field, the pose data of the engineering equipment, and candidate construction parameters to obtain a set of predicted responses for the underlying strata. Through the above processing, the settlement and deformation prediction results of the underlying strata under multiple candidate construction plans can be quickly obtained after each stage of operation, providing a quantitative basis for subsequent risk assessment and plan selection.

[0049] More specifically, in a concrete example of this application, the current three-dimensional coordinates and operating direction angle are first extracted from the pose data of the engineering equipment to determine the spatial range of the forward prediction interval. The latest geological mechanical parameters of each grid node within this interval are then read from the updated parameter field as the geological constraint input for the proxy simulation. Subsequently, each set of virtual operation schemes in the candidate construction parameters is traversed. Each scheme contains different construction control parameter settings (e.g., cutterhead torque, propulsion force, and soil chamber pressure in shield tunneling). The geological constraint input and each set of virtual operation schemes are sequentially fed into the deep neural network response proxy model to perform forward fast calculation, thereby obtaining the physical response prediction values ​​(e.g., segment deformation in shield tunneling) for each scheme, such as surface settlement, structural component deformation, and soil pressure distribution. The prediction results of all candidate schemes are aggregated and arranged according to the scheme index to form a forward geological response prediction set. The surface settlement prediction values ​​of each candidate scheme are calculated by spatial node difference to determine the building's uneven settlement, along with the maximum deformation of structural components and the excess value of soil pressure, and are included in the forward geological response prediction set.

[0050] In step S52, based on the risk assessment function and threshold constraints, the predicted set of forward strata responses is quantified for risk and schemes are selected to obtain construction parameter suggestions and early warning risk cloud map data. It should be noted that since the predicted set of forward strata responses contains the physical response estimates corresponding to multiple candidate construction schemes, it is necessary to quantitatively rank the risk levels of each scheme according to engineering safety control standards. Simultaneously, the spatial distribution information of the risk needs to be extracted into data suitable for subsequent rendering to support construction personnel's intuitive interpretation of the forward risk areas. Based on this, the technical solution of this application further quantifies the risk and selects schemes based on the risk assessment function and threshold constraints to obtain construction parameter suggestions and early warning risk cloud map data. Through the above processing, the construction parameter combination with the lowest overall risk can be selected from all candidate schemes, and cloud map data reflecting the spatial distribution of the forward strata risk level can be generated simultaneously.

[0051] More specifically, in a concrete example of this application, the predicted values ​​of each candidate scheme in the forward stratum response prediction set are first extracted for each evaluation index dimension. Evaluation indicators include maximum surface settlement, uneven settlement of buildings, maximum deformation of structural components, and earth pressure exceeding limits. A multi-dimensional risk assessment objective function is introduced to comprehensively score each scheme. For each candidate scheme, its predicted values ​​for each index are compared with the safety control thresholds specified in the engineering specifications. When the predicted value does not exceed the threshold, the risk penalty for that index is zero; when the predicted value exceeds the threshold, an asymmetric penalty of the square of the excess is applied. The scores are then weighted and summed using preset weighting coefficients for each index to obtain the comprehensive risk score for that scheme. After traversing all candidate schemes, the index of the scheme with the lowest comprehensive risk score is selected, and its corresponding construction parameter combination is extracted as a construction parameter suggestion. Simultaneously with the scheme selection, the spatial node coordinates of high-uncertainty areas in the updated parameter field and the grid nodes in the forward stratum response prediction set whose predicted values ​​are close to the safety critical threshold are extracted. Based on the risk level of each node, color gradient mapping is used to generate early warning risk cloud map data containing spatial coordinates and corresponding risk color values.

[0052] In step S53, the construction parameter suggestions, updated parameter fields, engineering equipment pose data, and early warning risk cloud map data are assembled and aligned to obtain the fusion data package to be rendered. It should be noted that since the geological parameter field, tunnel boring machine pose, risk cloud map, and construction suggestions belong to different data sources and coordinate systems, if they are directly sent to the rendering engine without unified assembly and coordinate alignment, spatial misalignment and hierarchical chaos will occur between the layers, making it impossible to form a consistent 3D visualization scene. Based on this, the technical solution of this application further assembles and aligns the construction parameter suggestions, updated parameter fields, engineering equipment pose data, and early warning risk cloud map data to obtain the fusion data package to be rendered. Through the above processing, the three types of heterogeneous rendering data—geological, equipment, and decision-making data—can be layered and encapsulated under the same global coordinate system, providing structured input for the subsequent unified rendering of the 3D engine.

[0053] More specifically, in a concrete example of this application, firstly, based on the updated parameter field, the latest mechanical parameter values ​​and corresponding stratum lithology information of each grid node are extracted, encoded and encapsulated according to the vertex attribute format required by the 3D rendering engine, and the mechanical parameters are mapped to node color channel values ​​and attached with material texture identifiers to form stratum rendering data. Subsequently, based on the pose data of the engineering equipment, the 3D spatial coordinates and attitude angles of the current engineering equipment are extracted, and the real-time translation and rotation of the engineering equipment model relative to the starting position in the initial twin scene are calculated. This spatial transformation relationship is bound to the 3D model of the engineering equipment, and historical trajectory spline curves are generated along the working direction to form equipment rendering data. On this basis, based on the early warning risk cloud map data and construction parameter suggestions, the risk level color value of each spatial node is encoded into a semi-transparent layer with transparency attributes, and the recommended operation values ​​in the construction parameter suggestions are encapsulated into a label information layer anchored to the current position of the tunnel boring machine. The two are then correlated to form decision rendering data. Finally, the stratum rendering data, device rendering data, and decision rendering data are aligned and verified in the global geographic coordinate system. After confirming that the spatial reference benchmarks of the three types of data are consistent, they are layered, assembled, and packaged according to the rendering priority of the stratum base map, device model, and decision overlay layer to obtain the rendering fusion data package.

[0054] In step S54, based on a 3D rendering engine, the construction parameter suggestions and the fusion data package to be rendered are dynamically rendered and presented to the interface to obtain a visual twin scene. It should be noted that since the geological parameter distribution, equipment pose trajectory, and risk warning information contained in the fusion data package to be rendered are all in structured numerical form, construction personnel cannot directly and quickly interpret the current working status and the distribution of risks ahead from the numerical data. It is necessary to convert it into an interactive 3D visual scene to support on-site decision-making. Based on this, the technical solution of this application further utilizes a 3D rendering engine to dynamically render and present the construction parameter suggestions and the fusion data package to be rendered to obtain a visual twin scene. Through the above processing, the updated geological state of each link, the real-time pose of the tunnel boring machine, the risk warning cloud map, and the construction parameter suggestions can be synchronously visualized in a unified 3D scene, forming a real-time digital twin interface that supports closed-loop decision-making.

[0055] More specifically, in a concrete example of this application, the data package to be rendered and merged is first pushed to a WebGL-based 3D rendering engine, which then loads the ground rendering data, equipment rendering data, and decision rendering data sequentially according to the layered assembly order. In the ground layer rendering stage, the mesh vertices and material textures in the ground rendering data are converted to 3DTiles format for streaming loading, and the surface shading and lighting calculations of the ground entity are completed based on the color mapping values ​​of the mechanical parameters of each node. In the equipment layer rendering stage, the 3D model of the tunnel boring machine is driven to complete the real-time position and attitude repositioning based on the spatial transformation relationships in the equipment rendering data, and the operation path marker line is drawn along the historical trajectory spline curve. In the decision layer rendering stage, the semi-transparent color value layer of the early warning risk cloud map is overlaid onto the ground surface using transparency blending, while the annotation information of the construction parameter suggestions is anchored and rendered to the spatial area above the current position of the engineering equipment. After completing the view frustum culling and frame buffer rendering of all layers, a visual twin scene supporting multi-view roaming and real-time refresh is output on the client. Construction personnel can intuitively view the spatial distribution of risk levels of the strata ahead and obtain corresponding construction parameter suggestions in this scene. After each step of the operation is completed and a new round of parameter updates is triggered, the visual twin scene is updated synchronously, forming a complete digital twin closed loop from real-time perception, parameter inversion, scene update to predictive decision-making.

[0056] In summary, a real-time update method for digital twin parameter fields based on ensemble Kalman filtering, according to embodiments of this application, is elucidated. This method establishes a stratigraphic entity mesh through spatial interpolation and triangulation reconstruction of spatially discrete sampled data, and spatially registers and integrates it with an engineering equipment model to form an initial twin scene. Based on this, finite element multi-condition simulation and Fisher information sensitivity analysis are used to optimize sensor placement, achieving efficient acquisition and synchronous preprocessing of heterogeneous signals, thus solving the problem of a lack of systematic acquisition methods for multi-source sensing information during construction. Subsequently, the observed signals and prior parameter sets are input into a response surrogate model, and robust assimilation and Bayesian inversion are implemented through adaptive noise adjustment and ensemble Kalman filtering, solving the problem of abnormal sensor readings interfering with inversion accuracy. In the parameter field update stage, a spatial compatibility weight matrix is ​​introduced to characterize the attenuation relationship of the influence of the working face location on each grid node. A mechanical coupling operator matrix is ​​introduced to constrain the linkage ratio between different formation parameters, and Laplace spatial smoothing is applied to solve the problems of parameter field spatial step jumps, mechanical parameter combination imbalances, and local noise peak residues caused by element-wise algebraic correction. This ensures that the update results take into account both the sensitivity of local corrections and the continuity of the overall field distribution. Finally, based on the updated parameter field, proxy inference and risk quantification of the forward formation response are performed, and the formation state, equipment trajectory, and early warning cloud map are dynamically presented through a 3D rendering engine. Visualized twin scenes and construction parameter suggestions are output, forming a complete closed loop from real-time perception, parameter inversion, intelligent updating to advanced prediction.

[0057] Those skilled in the art will understand that the steps, measures, and schemes in the various operations, methods, and processes discussed in this application can be alternated, modified, combined, or deleted; furthermore, other steps, measures, and schemes in the various operations, methods, and processes discussed in this application can also be alternated, modified, rearranged, decomposed, combined, or deleted; furthermore, the steps, measures, and schemes in the prior art that are similar to those disclosed in this application can also be alternated, modified, rearranged, decomposed, combined, or deleted. The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as the combination of these technical features does not contradict each other, it should be considered within the scope of this specification. The above-described embodiments are merely illustrative of several implementation methods of this disclosure, and their descriptions are relatively specific and detailed. However, they should not be construed as limiting the scope of the patent for the embodiments of this disclosure. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the embodiments of this disclosure, and these all fall within the protection scope of the embodiments of this disclosure. Therefore, the protection scope of the embodiments of this disclosure should be determined by the appended claims. As described above, although the present invention has been shown and described with reference to specific preferred embodiments, it should not be construed as limiting the present invention itself. Various changes in form and detail can be made without departing from the spirit and scope of the present invention as defined in the appended claims.

[0058] The present invention and its embodiments have been described above. This description is not restrictive, and the accompanying drawings are only one embodiment of the present invention; the actual structure is not limited thereto. In conclusion, if those skilled in the art are inspired by this description and design similar structures and embodiments without departing from the spirit of the present invention, such designs should fall within the protection scope of the present invention.

Claims

1. A method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering, characterized in that, include: S1, based on spatial discrete sampling data and engineering equipment model data, performs spatial interpolation and triangulation reconstruction on the stratigraphic interface control points, and performs spatial registration and scene integration on the stratigraphic entity mesh and engineering equipment model to obtain the initial twin scene; S2, based on the initial twin scenario and preset operating parameters, performs numerical simulation and sensitivity analysis on the position response of candidate sensors, and optimizes sensor layout, heterogeneous signal acquisition and synchronous preprocessing based on the analysis results to obtain the observation signal vector; S3, based on the observed signal vector, prior parameter set and response surrogate model, performs robust assimilation and Bayesian inversion on the predicted signal and the measured signal, and combines adaptive noise adjustment to perform set recursive update to obtain the posterior parameter set; S4, based on the posterior parameter set, historical operation sequence, previous loop parameter field and engineering equipment pose data, performs joint discrimination and mask update on parameter information gain and historical evolution law, and performs selective correction of parameter field according to the discrimination result to obtain updated parameter field; S5, based on the updated parameter field, engineering equipment pose data and candidate construction parameters, performs evolution prediction and risk assessment of the ground response ahead, and combines the updated parameter field for 3D rendering and parameter feedback to obtain a visualized twin scene and construction parameter suggestions.

2. The method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering according to claim 1, characterized in that, Step S1 includes: Based on spatial discrete sampling data, the sampling spatial location and stratigraphic layer information are extracted and expanded, and auxiliary constraint points are generated along the normal direction of the interface control points to obtain the control point dataset. Implicit surface space interpolation based on radial basis functions is performed on the control point dataset to obtain a 3D predicted point cloud; Solid reconstruction and material assignment are performed on the 3D predicted point cloud to obtain the geological solid mesh; Based on coordinate transformation and spatial registration, scene integration and encapsulation are performed on the geological entity mesh and engineering equipment model data to obtain an initial twin scene.

3. The method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering according to claim 1, characterized in that, Step S2 includes: Based on the initial twin scenario and preset operating parameters, finite element modeling and multi-condition simulation of the equipment-stratum coupling relationship are performed, and candidate position response results are extracted to obtain a simulation response database. Information gain evaluation and location filtering are performed on the simulation response database to obtain the optimal sensor placement scheme; Based on the optimal sensor layout scheme, sensors are deployed and online data is collected at the locations of the engineering equipment body, tunnel segments, and ground surface to obtain raw heterogeneous signals. The original heterogeneous signals are subjected to signal denoising and time alignment synchronization preprocessing to obtain the observed signal vector.

4. The method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering according to claim 3, characterized in that, Information gain evaluation and location filtering are performed on the simulation response database to obtain the optimal sensor placement scheme, including: Based on the simulation response database, the finite difference method is used to calculate the sensitivity of the sensor response at each candidate position to the target inversion parameters to obtain the Fisher information matrix. Based on the D-optimal criterion, the Fisher information matrix is ​​greedily maximized to obtain the optimal sensor placement scheme.

5. The method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering according to claim 1, characterized in that, Step S3 includes: Based on the response surrogate model, forward prediction and signal mapping are performed on the prior parameter set to obtain the predicted signal set; Based on the observed signal vector and the predicted signal set, anomaly detection and noise adjustment are performed on the sensor innovation statistics to obtain an adaptive observation noise covariance matrix. The Kalman gain matrix is ​​obtained by solving the gain of the prior parameter set, the predicted signal set, and the adaptive observation noise covariance matrix through covariance statistics and localization constraints. Based on the Kalman gain matrix, the observed signal vector, the prior parameter set, and the predicted signal set are recursively updated and expanded to obtain the posterior parameter set.

6. The method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering according to claim 5, characterized in that, Anomaly detection and noise conditioning are performed on sensor innovation statistics to obtain an adaptive observation noise covariance matrix, including: Based on sensor accuracy calibration, a basic noise variance is set for each sensor to obtain the initial observation noise covariance matrix; Based on the observed signal vector and the predicted signal set, the normalized squared innovation statistic is calculated for each sensor and compared with the chi-square distribution threshold to obtain the abnormal sensor identification. The basic noise variance corresponding to the abnormal sensor identifier is expanded and corrected to obtain the adaptive observation noise covariance matrix.

7. The method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering according to claim 1, characterized in that, Step S4 includes: Based on statistical analysis and standard deviation comparison, mean extraction and information gain calculation are performed on the posterior parameter set to obtain the posterior mean of the parameters and the Bayesian uncertainty reduction rate vector. Based on historical job sequences and long short-term memory networks, time-series modeling and pattern extraction are performed on the posterior mean of parameters and the Bayesian uncertainty reduction rate vector to obtain the parameter update priority vector. A mask fusion is performed on the Bayesian uncertainty reduction rate vector and the parameter update priority vector to obtain a dynamic update mask; Based on the dynamically updated mask, the posterior mean of the parameters and the parameter field of the previous loop are gated and spatially mapped to obtain the updated parameter field.

8. The method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering according to claim 1, characterized in that, Step S5 includes: Based on the updated parameter field, engineering equipment pose data and candidate construction parameters, proxy simulation and response prediction of the front stratum state are performed to obtain the front stratum response prediction set. Based on the risk assessment function and threshold constraints, risk quantification and scheme selection are performed on the front stratum response prediction set to obtain construction parameter suggestions and early warning risk cloud map data; The construction parameter suggestions, updated parameter fields, engineering equipment pose data, and early warning risk cloud map data are assembled and aligned to obtain the fusion data package to be rendered. Based on a 3D rendering engine, construction parameter suggestions and data packages to be rendered are dynamically rendered and presented on the interface to obtain a visual twin scene.

9. The method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering according to claim 8, characterized in that, The construction parameter suggestions, updated parameter fields, engineering equipment pose data, and early warning risk cloud map data are assembled and aligned to obtain the fusion data package to be rendered, including: Based on the updated parameter field, the parameters of the formation grid nodes and the formation category information are mapped and encapsulated to obtain the formation rendering data; Based on the pose data of engineering equipment, the position, attitude and operation trajectory of the engineering equipment model are bound and updated to obtain equipment rendering data; Based on the early warning risk cloud map data and construction parameter suggestions, the risk level information and parameter suggestion information are associated and encapsulated in layers to obtain decision rendering data; Coordinate alignment and layered assembly are performed on the ground rendering data, device rendering data, and decision rendering data to obtain the fused data package to be rendered.

10. The method for real-time updating of digital twin parameter fields based on ensemble Kalman filtering according to claim 7, characterized in that, Based on a dynamically updated mask, gating correction and spatial mapping are performed on the posterior mean of the parameters and the parameter field of the previous loop to obtain the updated parameter field, including: Based on the pose data of engineering equipment and grid coordinates, the spatial distance between the update area and the work surface is estimated by relevant weights to obtain a spatial weight matrix. Based on the dynamically updated mask, spatial weight matrix, posterior mean of parameters and parameter field of the previous loop, the parameter evolution increment is coupled with constraints and compatibility correction to obtain the correction increment. The corrected increment and the parameter field of the previous loop are continuously fused and mapped to nodes to obtain the updated parameter field.

Images