Highway tunnel construction ground surface settlement prediction system based on digital twinning

Abstract

The application discloses a highway tunnel construction ground surface settlement prediction system based on digital twinning, which comprises a twinning modeling module, a tunnel face digital twinning grid model and a ground surface settlement digital twinning grid model are constructed, and a dynamic correlation is established through a space-time coordinate mapping relationship; a data acquisition module, which acquires real-time tunneling parameter sequences and ground surface settlement monitoring sequences in the tunneling process; a load mapping module, which maps the real-time tunneling parameter sequences to corresponding nodes of the tunnel face digital twinning grid model to generate a tunnel face dynamic load distribution field; a time series extrapolation module, which performs time series extrapolation processing on a ground surface settlement prediction cloud image, combines the spatial distribution characteristics of key settlement area identification, and generates a settlement rate evolution curve and a settlement threshold overrun time prediction value of each monitoring section. The system realizes bidirectional dynamic iterative prediction of tunneling load and settlement, and gives an overrun time in advance, thereby providing support for construction parameter optimization.

Description

Technical Field

[0001] This invention relates to the field of surface settlement prediction technology in tunnel construction, specifically a surface settlement prediction system for highway tunnel construction based on digital twins. Background Technology

[0002] During highway tunnel construction, the excavation operation disturbs the strata and triggers surface settlement. Accurate prediction of this settlement is crucial for ensuring construction safety. Traditional methods, based on empirical formulas or two-dimensional single-step finite element analysis, treat the strata as a static loading object, failing to reflect the dynamic interaction between real-time changes in excavation parameters and the continuous evolution of surface settlement. With the introduction of the digital twin concept, some schemes collect parameters such as thrust and torque during the excavation process and map them unidirectionally onto the tunnel's three-dimensional mesh model, driving displacement field updates and generating predicted surface settlement results. This unidirectional mapping mode does not feed back the surface settlement that has already occurred to the stress balance system of the tunnel surrounding rock. This results in a lack of real-time coupling between excavation loading, strata deformation, settlement feedback, and load redistribution in the simulation process. Especially when encountering weak interlayers or strata with large displacements, surface settlement has a significant reaction effect on the spatial distribution of subsequent loads at the excavation face. Conventional methods cannot capture this closed-loop effect, leading to significant prediction bias. Current solutions only output settlement distribution cloud maps or linear extrapolation curves at a certain moment, lacking the ability to analyze and extrapolate future development trends from the perspective of settlement differentiation rates. They cannot predict when settlement at each monitoring section will exceed the preset safety threshold, nor can they provide quantitative basis directly for tunneling parameter control. Construction teams can only take passive remedial measures after settlement exceeds the limit alarm, delaying the optimal adjustment window. The lack of a two-way coupled iterative and accurate prediction mechanism between the dynamic load at the tunnel face and the surface settlement field, and the inability to extrapolate time-series rates based on the coupling results, calibrate the time of exceeding limits, and generate adjustment amounts for construction parameters in reverse, are problems that urgently need to be solved by current technology. Summary of the Invention

[0003] This paper proposes a digital twin-based surface settlement prediction system for highway tunnel construction to address the shortcomings of existing technologies, such as the lack of a two-way real-time iterative coupling mechanism between tunnel face load and surface settlement leading to inaccurate predictions, and the inability to provide early warnings of settlement rate evolution and over-limit moments at monitoring sections.

[0004] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides a surface settlement prediction system for highway tunnel construction based on digital twins. The system includes a twin modeling module, a data acquisition module, a load mapping module, a coupled calculation module, and a time series extrapolation module.

[0005] The twin modeling module is used to construct digital twin mesh models of the tunnel face and surface settlement, and to establish a dynamic relationship between the two models through a spatiotemporal coordinate mapping relationship. Preferably, the construction process is as follows: a three-dimensional hexahedral mesh structure of the tunnel face is generated with the tunnel design axis as the baseline, and each mesh node stores the node coordinates and the lithological parameters of the strata to which it belongs; a triangular mesh surface model of the surface is constructed based on the original topographic point cloud data of the surface, and it is resampled into a regular rectangular mesh structure, with each mesh point storing the point coordinates and initial elevation value; then, a spatial interpolation mapping table is established between the three-dimensional hexahedral mesh nodes and the regular rectangular mesh points, and the table contains the weight coefficient of each three-dimensional hexahedral mesh node to the neighboring regular rectangular mesh points.

[0006] The data acquisition module is used to collect real-time tunneling parameter sequences and surface settlement monitoring sequences during the tunnel excavation process. Specifically, using the center point of the tunnel boring machine cutterhead as the reference origin, it collects cutterhead thrust, cutterhead torque, and tail grouting pressure parameters along the excavation direction at a fixed excavation ring spacing, generating a real-time tunneling parameter sequence. Simultaneously, a settlement monitoring sensor array deployed at various points in a regular rectangular grid on the ground surface collects real-time settlement displacement values ​​at each grid point at equal time intervals, generating a surface settlement monitoring sequence. The parameter recording time is time-aligned with the nearest monitoring time to generate synchronized tunneling settlement parameter pairs.

[0007] The load mapping module maps real-time tunneling parameter sequences to corresponding node positions in the digital twin mesh model of the tunnel face, generating a dynamic load distribution field for the tunnel face. Based on the tunneling ring spacing corresponding to the parameter recording time, the set of node layers representing the current tunnel face in the model is determined. The cutterhead thrust parameters are weighted according to a Gaussian function, with the central region nodes receiving a higher proportion of thrust than the edge region nodes. The cutterhead torque parameters are converted into tangential friction load components and applied to the outer edge nodes of the tunnel face. The tail grouting pressure parameters are converted into normal expansion load components and applied to nodes within a fixed ring spacing range behind the node layer set.

[0008] The coupled calculation module is used to perform bidirectional coupled iterative calculations based on the boundary condition constraints between the dynamic load distribution field of the tunnel face and the digital twin mesh model of surface settlement, generating a predicted surface settlement cloud map and key settlement area markers. As a technical solution of the present invention, the calculation process includes: inputting the load values ​​of each node in the dynamic load distribution field of the tunnel face as displacement boundary conditions into the equilibrium equation of the digital twin grid model of the tunnel face, and calculating the initial displacement response field of the rock mass surrounding the tunnel; extracting the displacement response components of the surface location from the field, and transferring them to each grid point of the digital twin grid model of surface settlement according to the weight coefficients in the spatial interpolation mapping table, generating the surface displacement prediction increment; superimposing the prediction increment onto the current surface elevation value, updating the spatial morphology of the digital twin grid model of surface settlement, and obtaining the updated surface settlement distribution field; comparing the residual vectors of the surface settlement distribution fields before and after the update, if the residual modulus is greater than the preset convergence threshold, then feeding the updated distribution field back to the tunnel face model as a new boundary condition to recalculate the displacement response field; iterating until the residual modulus is less than or equal to the convergence threshold, outputting the converged surface settlement distribution field as a prediction cloud map, and extracting the continuous grid point area with settlement greater than the preset risk threshold as the key settlement area identifier. Preferably, when the number of iterations exceeds the preset maximum number of iterations, the calculation is forcibly terminated and the current surface subsidence distribution field is output as the suboptimal prediction result.

[0009] As a preferred embodiment of the present invention, a material parameter feedback correction step is included before performing the above-mentioned coupled calculation. The initial values ​​of the elastic modulus and Poisson's ratio corresponding to the lithological parameters of each stratum in the digital twin mesh model of the tunnel face are obtained to construct an initial material parameter field. Based on the difference between the actual settlement and the predicted settlement at the corresponding time in the surface settlement monitoring sequence, the settlement prediction error field at each monitoring time is calculated. The settlement prediction error field is back-projected to the node positions of the digital twin mesh model of the tunnel face according to the inverse mapping relationship of the spatial interpolation mapping table to generate the elastic modulus correction coefficient for each node. The initial value of the elastic modulus of each node is multiplied by the corresponding correction coefficient to generate an updated material parameter field and replace the initial field. The calculation method for the settlement prediction error field is as follows: extract the actual settlement displacement values ​​of the latest preset number of monitoring times from the surface settlement monitoring sequence to form a measured settlement vector, and extract the predicted settlement displacement values ​​of the corresponding times to form a predicted settlement vector. Calculate the element-wise difference between the two to generate a difference vector, and then perform Kriging space interpolation on the difference vector to generate a continuous error distribution surface covering the entire regular rectangular grid area, which serves as the settlement prediction error field for the current monitoring time.

[0010] The temporal extrapolation module is used to perform temporal extrapolation processing on the surface subsidence prediction cloud map. Combining the spatial distribution characteristics of key subsidence area markers, it generates subsidence rate evolution curves and predicted times when subsidence thresholds exceed limits for each monitoring section. Subsidence trough curves at each monitoring section location are extracted from the surface subsidence prediction cloud map; these curves follow the Peck formula. Based on historical surface subsidence monitoring sequences from multiple monitoring times, the least squares method is used to fit the subsidence trough depth growth coefficient and subsidence trough width expansion coefficient for each monitoring section, establishing a function for the subsidence trough curve's variation with increasing tunneling ring spacing. Substituting the tunneling ring spacing values ​​for the future prediction period into the function, the predicted subsidence trough curves for each section at future times are calculated. The sequence of maximum subsidence variation over time is extracted, and after differentiation, the rate of change of maximum subsidence at each time point is generated, i.e., the subsidence rate evolution curve.

[0011] Simultaneously, a continuous set of grid points covered by the key settlement area markers is obtained, and the average value of the settlement rate evolution curves of each grid point in this set is calculated to generate the average settlement rate curve for the key area. The starting point of the acceleration segment of the rate is determined based on the difference between the rate values ​​at the current time and the previous time, and the curve is extrapolated and extended from this point. Set a settlement accumulation threshold and a settlement rate threshold, compare the average settlement rate curve of the key area with the settlement rate threshold, and record the crossing time of the curve across the rate threshold; accumulate the predicted settlement values ​​of each grid point in the key settlement area markers to the accumulation threshold, and record the accumulation time when the accumulation threshold is reached; compare the crossing time and the accumulation time, and take the smaller one as the predicted value when the settlement threshold exceeds the limit.

[0012] Furthermore, after generating the predicted time when the settlement threshold of the monitoring section exceeds the limit, the system compares the difference between this predicted value and the planned passage time of the tunnel excavation plan section to calculate the safe construction margin time. If the safe construction margin time is less than the preset safe margin threshold, the system back-calculates the required adjustment amount of the tunneling parameters based on the settlement rate evolution curve, including the cutterhead thrust reduction ratio and the shield tail grouting pressure increase ratio, and transmits this adjustment amount to the tunnel boring machine control system to guide the setting of subsequent tunneling ring spacing construction parameters, thereby achieving active control of surface settlement.

[0013] The technical effects and advantages provided by the present invention in the above technical solution are as follows: The real-time acquired cutterhead thrust is applied to the node layers of the corresponding annular surface of the tunnel face with radial weights distributed according to a Gaussian function, with higher weights at the center and lower weights at the edges. The cutterhead torque is converted into a tangential friction load and applied to the outer edge nodes, while the tail grouting pressure is converted into a normal expansion load and applied to the rear annular nodes, generating a dynamic load distribution field at the tunnel face. During coupled calculation, this load field is used as the displacement boundary condition. The tunnel surrounding rock equilibrium equation is input to solve for the initial displacement response. The surface component is extracted based on the weight coefficients of the spatial interpolation mapping table and superimposed onto the surface settlement grid to update the surface elevation and generate an updated settlement field. The norm of the residual vector of the two settlement fields is compared with a preset convergence threshold. If convergence is not met, the updated settlement field is fed back to the tunnel grid as a new boundary constraint, and the displacement calculation is performed again. The process is iterated until convergence or the maximum number of iterations is reached, and a suboptimal result is output. This two-way iterative coupling mechanism enables surface subsidence changes to constrain the stress rebalancing and load transfer path of the tunnel surrounding rock, solving the problem of neglecting the subsidence feedback effect in traditional one-way prediction. It has a stronger ability to capture anomalous changes in subsidence distribution when traversing heterogeneous rock strata, and the consistency between predicted cloud maps and actual measurements is significantly improved.

[0014] In the temporal extrapolation process, settlement troughs for each section are extracted from the predicted cloud map. The Peck formula is used to perform least-squares fitting on the historical measured settlement trough sequence to determine the depth growth coefficient and width expansion coefficient of the settlement troughs. For the grid points covered by the key settlement area markers, the average settlement rate curve is calculated, the starting point of the acceleration segment is identified and extrapolated, and the curve is compared with the cumulative settlement threshold and the rate threshold respectively. The smaller value between the crossing time and the cumulative time is taken as the predicted value of the over-limit time. The difference between this over-limit time and the planned tunnel passage time is compared to obtain the safe construction margin time. When the margin is lower than the preset safety threshold, the required cutterhead thrust reduction ratio and shield tail grouting pressure increase ratio are calculated back based on the rate curve, and the adjustment amount is automatically generated and sent to the tunnel boring machine control system. This achieves accurate prediction of the entire evolution of the future settlement rate and the risk moments for each section, transforming settlement control from empirical trial-and-error to precise intervention through quantitative calculation, and preventing over-limit accidents caused by parameter adjustment lags. Attached Figure Description

[0015] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0016] Figure 1 This is a schematic diagram of a digital twin-based surface settlement prediction system for highway tunnel construction. Figure 2This is a flowchart of the construction process for a digital twin mesh model of the tunnel face and surface settlement; Figure 3 This is a flowchart of the generation process for the dynamic load distribution field at the tunnel face based on tunneling parameters and settlement monitoring sequences. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] See Figure 1 This invention provides a digital twin-based system for predicting surface settlement during highway tunnel construction. The system includes a twin modeling module, a data acquisition module, a load mapping module, a coupled calculation module, and a time-series extrapolation module. The twin modeling module constructs a digital twin mesh model of the tunnel face and a digital twin mesh model of surface settlement, establishing a dynamic relationship between them through a spatiotemporal coordinate mapping. The data acquisition module collects real-time tunneling parameter sequences and surface settlement monitoring sequences during tunnel excavation. The load mapping module maps the real-time tunneling parameter sequences to corresponding node positions in the digital twin mesh model of the tunnel face, generating a dynamic load distribution field for the tunnel face. The coupled calculation module performs bidirectional coupled iterative calculations based on the boundary condition constraints between the dynamic load distribution field of the tunnel face and the digital twin mesh model of surface settlement, generating a predicted surface settlement cloud map and identifying key settlement areas. The temporal extrapolation module performs temporal extrapolation processing on the predicted cloud map of land subsidence, and combines the spatial distribution characteristics of key subsidence area markers to generate subsidence rate evolution curves and predicted values ​​of subsidence threshold exceedance times for each monitoring section.

[0019] In specific implementation, please refer to Figure 2 The steps for constructing digital twin mesh models of the tunnel face and surface settlement are implemented in the following way.

[0020] Acquire tunnel design axis data and original surface topographic point cloud data. Tunnel design axis data is a three-dimensional coordinate sequence of the tunnel centerline in the construction coordinate system, with each coordinate point accompanied by a mileage value, representing the arc length distance from the starting point along the tunnel excavation direction. Original surface topographic point cloud data is acquired by a 3D laser scanner or aerial photogrammetry system. The point cloud data contains the three-dimensional coordinates of discrete points on the surface, with a point cloud density of no less than four points per square meter.

[0021] A three-dimensional hexahedral mesh structure for the tunnel excavation face is generated using the tunnel design axis as a baseline. Tunnel cross-sections are extracted along the tunnel design axis at fixed mileage intervals. These fixed mileage intervals are set based on the tunnel boring machine's (TBM) excavation ring distance, typically half of the ring distance. On each cross-section, multiple layers of quadrilateral rings are generated outwards from the projection point of the tunnel design axis onto the cross-section center, forming a two-dimensional quadrilateral mesh. The radial dimension of the quadrilateral mesh is determined by the tunnel excavation radius and the expected disturbance zone range. Corresponding quadrilateral nodes on adjacent mileage cross-sections are connected along the mileage direction to form three-dimensional hexahedral units. All three-dimensional hexahedral units are combined to form a three-dimensional hexahedral mesh structure. Each grid node in the three-dimensional hexahedral mesh structure stores node coordinates and the lithological parameters of the strata to which the node belongs. The node coordinates are determined by the cross-sectional position, radial distance, and circumferential angle. The lithological parameters of the strata to which the node belongs are calculated and determined based on the three-dimensional spatial equations of the stratigraphic interfaces in the geological exploration report and the node coordinates. The node coordinates are substituted into the stratigraphic interface equations, and the stratigraphic position of the node is determined according to the positive or negative sign of the calculated values. Then, the corresponding stratigraphic lithological parameters are assigned, including cohesion, internal friction angle, elastic modulus, and Poisson's ratio.

[0022] A triangular mesh surface model of the land surface is constructed based on the original topographic point cloud data. After denoising and thinning preprocessing of the original topographic point cloud data, the Delaunay triangulation algorithm is used to connect discrete points into a set of triangular facets. Each triangular facet is determined by the coordinates of three vertices and a normal vector. Triangular facets share edges and do not overlap, forming a continuous triangular mesh surface model covering the entire surface monitoring area. The triangular mesh surface model is then resampled into a regular rectangular mesh structure at a preset resolution. The preset resolution is determined based on the deployment spacing of the surface subsidence monitoring sensor array, ranging from one-tenth to one-fifth of the deployment spacing to ensure the ability to characterize the spatial features of subsidence deformation. The regular rectangular grid structure is formed by the intersection of longitudinal grid lines along the tunnel axis and transverse grid lines perpendicular to the tunnel axis. Each grid point stores the point coordinates and the initial elevation value. The point coordinates are the planar coordinates of the grid line intersections. The initial elevation value is calculated by interpolation of the vertical projection of the triangular mesh surface model at the point coordinates. The interpolation algorithm uses centroid coordinate interpolation, that is, determining the triangle in which the point coordinates are located in the triangular mesh surface model, and using the coordinates and elevation values ​​of the three vertices of the triangle to calculate the point elevation value.

[0023] A spatial interpolation mapping table is established between the node numbers of a 3D hexahedral mesh structure and the mesh point numbers of a regular rectangular mesh structure. For each node in the 3D hexahedral mesh structure, the node coordinates are extracted. A search radius is set with the horizontal projection point of the node coordinates as the center. The search radius is the minimum side length of the regular rectangular mesh structure multiplied by a preset multiple, ranging from two to five. All regular rectangular mesh points are selected within this search radius. For each found regular rectangular mesh point, the spatial distance between the 3D hexahedral mesh node and the regular rectangular mesh point is calculated using the Euclidean distance formula. Weighting coefficients are calculated based on the spatial distance using an inverse distance weighted form. The weighting coefficient is equal to the negative power of the spatial distance, with an exponent of two. All weighting coefficients are normalized so that the sum of the weighting coefficients of all regular rectangular mesh points corresponding to a 3D hexahedral mesh node is one. The spatial interpolation mapping table contains three records: the 3D hexahedral mesh node number, the regular rectangular mesh point number, and the weighting coefficient. The node numbers of the 3D hexahedral mesh and the regular rectangular mesh are encoded using consecutive integers starting from zero. Each 3D hexahedral mesh node corresponds to multiple entries, and each entry records a matching regular rectangular mesh node number and its corresponding weight coefficient. The spatial interpolation mapping table is stored in system memory for quick lookup during subsequent coupled calculations.

[0024] In specific implementation, please refer to Figure 3 The real-time tunneling parameter sequence and surface settlement monitoring sequence during the tunnel excavation process are collected through the following methods.

[0025] Using the center point of the tunnel boring machine's (TBM) cutterhead as the reference origin, located at the geometric center of the cutterhead's front face, its spatial coordinates are updated in real time with the tunneling mileage. Cutterhead thrust, cutterhead torque, and tail grouting pressure parameters are collected along the tunneling direction at a fixed tunneling ring spacing. The fixed tunneling ring spacing is the distance advanced by the TBM in one complete tunneling cycle, determined according to the TBM model and construction process, and pre-set in the construction parameter configuration file. During each data collection, the current values ​​of the cutterhead thrust sensor, cutterhead torque sensor, and tail grouting pressure sensor are read from the TBM's programmable logic controller (PLC). The cutterhead thrust sensor range is no less than 120% of the TBM's maximum design thrust, the cutterhead torque sensor range is no less than 120% of the maximum design torque, and the tail grouting pressure sensor range is no less than 200% of the maximum grouting pressure. The three read values ​​are combined with the tunneling mileage value corresponding to the current cutterhead center point to form a record, generating a real-time tunneling parameter sequence according to the collection order. Each tunneling ring distance corresponds to a parameter recording time. The parameter recording time is the time when the tunneling process ends and the segment assembly begins. This time point is triggered by the travel counter of the tunnel boring machine's control system.

[0026] Settlement monitoring sensor arrays are deployed at each grid point of a regular rectangular grid structure on the ground surface. The settlement monitoring sensors utilize hydrostatic levels or total stations with automatic monitoring modules. One sensor is installed at each regular rectangular grid point, with its base buried at a fixed depth below the ground surface. This fixed depth is determined based on the local frost line depth, and is no less than 0.3 meters below the frost line. The sensors collect real-time settlement displacement values ​​for each grid point at equal time intervals. These time intervals are set according to the construction stage and the sensitivity of the soil layer, with intervals no greater than five minutes in critical crossing sections and no greater than fifteen minutes in general sections. After each data collection, the sensors transmit the data wirelessly to a field data aggregation base station. The aggregation base station constructs a vector from the settlement displacement values ​​of all grid points collected at the same time. The length of the vector is equal to the total number of regular rectangular grid points, and the position of each element in the vector corresponds one-to-one with the grid point number, generating a surface settlement monitoring sequence. Each monitoring moment corresponds to a set of settlement displacement vectors for all grid points, and the monitoring moment is marked as the average value of the time when the aggregation base station receives all the data collected by the sensors for this period.

[0027] The time of each parameter recording time is time-aligned with the nearest monitoring time. For each parameter recording time in the real-time tunneling parameter sequence, all monitoring times in the surface settlement monitoring sequence are traversed, the absolute time difference between the parameter recording time and the monitoring time is calculated, and the monitoring time with the smallest absolute time difference is selected as the matching monitoring time. The cutterhead thrust, cutterhead torque, and tail grouting pressure parameters of the matched parameter recording time are combined with the settlement displacement vector of the full grid points at that monitoring time to form a synchronized tunneling settlement parameter pair. The time alignment process is completed in the system data acquisition module, and the synchronized tunneling settlement parameter pair is used for subsequent load mapping and coupling calculations.

[0028] In practice, the steps for mapping the real-time tunneling parameter sequence to the corresponding node positions of the digital twin mesh model of the tunnel face to generate the dynamic load distribution field of the tunnel face are as follows.

[0029] Based on the tunneling ring distance corresponding to the recorded time of each parameter in the real-time tunneling parameter sequence, the set of node layers corresponding to the current tunneling face position in the digital twin mesh model of the tunnel face is determined. The tunneling ring distance is obtained from the tunnel boring machine control system and converted into a tunneling mileage value, which represents the curved distance of the cutterhead center point along the tunnel design axis. In the digital twin mesh model of the tunnel face, all three-dimensional hexahedral mesh nodes store their respective mileage coordinate components. Centered on the tunneling mileage value, a longitudinal search tolerance is set, which is half the unit size of the three-dimensional hexahedral mesh structure in the mileage direction. All nodes whose mileage coordinates fall within the range of tunneling mileage value minus the longitudinal search tolerance and tunneling mileage value plus the longitudinal search tolerance are selected. These nodes constitute the node layer set. The node layer set contains all three-dimensional hexahedral mesh nodes located within the same tunneling mileage section, including excavation face boundary nodes and surrounding rock nodes within a certain radial range from the excavation face.

[0030] The cutterhead thrust parameters are distributed to each node in the node layer set according to the radial distance from the center to the edge of the tunnel face. The radial distance of each node in the node layer set is extracted; the radial distance is the planar projection distance of the node coordinates from the cutterhead center point. A Gaussian function is used for the allocation weights, and the specific formula for calculating the allocation weights is as follows: in, This represents the radial distance between nodes, in meters. Indicates distance from the center as The weight values ​​assigned to the nodes are dimensionless. The standard deviation of the Gaussian function. The value is taken as one-third of the excavation radius of the tunnel face, which is determined based on the outer diameter of the tunnel's designed cross-section. This is the base of the natural logarithm. During the calculation, the radial distance of each node in the node layer set is first calculated. Substituting into the above formula yields the weight value. Then, the weights of all nodes are summed to obtain the total weight. The weight of each node is then normalized by dividing the total weight, resulting in a normalized allocation weight. The cutterhead thrust parameter value is multiplied by the normalized allocation weight to obtain the thrust load value applied to that node. The radial distance to the center node of the tunnel face is zero, resulting in the largest calculated allocation weight; the radial distance to the edge nodes is close to the excavation radius, resulting in the smallest allocation weight. The proportion of allocated thrust values ​​to nodes in the central region is higher than that to nodes in the edge region, which conforms to the actual contact pressure distribution between the cutterhead and the tunnel face.

[0031] The cutterhead torque parameter is converted into a tangential friction load component. The cutterhead torque reflects the cutting and torsional action of the cutterhead and scraper on the rock and soil at the tunnel face. The cutterhead torque parameter value is divided by the number of nodes at the outer edge of the tunnel face, and then by the excavation radius of the tunnel face, to obtain the tangential force amplitude acting on each outer edge node. For nodes located at the outer edge of the tunnel face in the node layer set, which are nodes with a radial distance equal to the excavation radius, a tangential direction vector is defined. This tangential direction vector is the tangential direction of the node's rotation around the cutterhead center within the excavation section, and its direction is obtained by rotating the node's position vector relative to the cutterhead center by 90 degrees around the tunnel mileage direction axis. The tangential force amplitude is multiplied by the tangential direction vector to generate the tangential friction load component, which is then applied to the corresponding outer edge node.

[0032] The tail grouting pressure parameters are converted into normal expansion load components. The tail grouting pressure is the radial pressure formed by the grout within the tail void during synchronous tail grouting. In the digital twin mesh model of the tunnel face, nodes within a fixed ring spacing range behind the node layer set are determined. The fixed ring spacing is determined based on the geometric distance from the tail to the cutterhead, which is the distance from the tail grouting port to the center of the cutterhead, obtained from the tunnel boring machine design drawings. Within the nodes within this fixed ring spacing range, the normal direction vector at the corresponding position of each node is extracted. The normal direction vector is the perpendicular direction from the node position to the tunnel design axis and points towards the interior of the surrounding rock. The tail grouting pressure parameter is applied in a uniform or non-uniform circumferential distribution manner. When uniformly distributed, the tail grouting pressure value is used as the normal load amplitude, multiplied by the normal direction vector of each node to generate the normal expansion load component, which is then applied to the corresponding nodes. When non-uniformly distributed, the influence of grout gravity is considered, and the tail grouting pressure value is linearly corrected at the crown and crown according to the hydrostatic pressure gradient, and then multiplied by the normal direction vector of each node to generate the load component.

[0033] After completing the cutterhead thrust distribution, cutterhead torque conversion and application, and shield tail grouting pressure conversion and application, the three types of load components experienced by all nodes at each parameter recording time are vector-superimposed to obtain the resultant load value of each node in the digital twin mesh model of the tunnel face at that time. The resultant load values ​​of all nodes constitute the dynamic load distribution field of the tunnel face. The dynamic load distribution field of the tunnel face is dynamically updated according to the recording time of each parameter in the real-time tunneling parameter sequence.

[0034] In practice, based on the boundary condition constraints between the dynamic load distribution field of the tunnel face and the digital twin grid model of surface settlement, a two-way coupled iterative calculation process is performed to generate a surface settlement prediction cloud map and key settlement area markers through the following steps.

[0035] The nodal load values ​​in the dynamic load distribution field of the tunnel face are input as displacement boundary conditions into the equilibrium equations of the digital twin mesh model of the tunnel face. The equilibrium equations of the digital twin mesh model of the tunnel face are established based on the finite element method of elasticity, and the discretized form of the equilibrium equations is as follows: in, The overall stiffness matrix of the digital twin mesh model of the tunnel face is represented by force divided by length. The overall stiffness matrix is ​​obtained by assembling the element stiffness matrices of each element in the three-dimensional hexahedral mesh structure according to the node number. The element stiffness matrix is ​​calculated by integrating the shape function of the three-dimensional hexahedral isoparametric element based on the elastic modulus and Poisson's ratio in the lithological parameters of the strata to which each element node belongs. This represents the nodal displacement vector to be solved, with the dimension of length. The nodal displacement vector contains the three orthogonal displacement components of all nodes in the three-dimensional hexahedral mesh structure. The equivalent nodal load vector, with the dimension of force, is composed of the resultant load values ​​of each node in the dynamic load distribution field of the tunnel face, arranged according to node number. For nodes without applied loads, the corresponding element in the equivalent nodal load vector is set to zero. When displacement boundary conditions are applied, the outer boundary nodes of the three-dimensional hexahedral mesh structure are constrained according to the actual constraints. The displacements in all three directions of the bottom boundary nodes are constrained to zero, the horizontal normal displacement of the side boundary nodes is constrained to zero, and the top boundary nodes are treated as free boundaries before being covered by the surface settlement digital twin mesh model. Solving the equilibrium equations yields the initial displacement response field of the rock mass surrounding the tunnel, which is the set of displacement components of all nodes in three directions.

[0036] Displacement response components of the surface location are extracted from the initial displacement response field. In the digital twin mesh model of the tunnel face, the surface location corresponds to the top node of a three-dimensional hexahedral mesh structure. The vertical displacement component of the top node is extracted, with the vertical direction representing the gravity direction and downward displacement being positive. The extracted vertical displacement component of the top node is used as the source data for the surface displacement response. The surface displacement response components are then transferred to each mesh point of the surface settlement digital twin mesh model according to the weighting coefficients in the spatial interpolation mapping table. The spatial interpolation mapping table is constructed using the above embodiment and stored in system memory. For each regular rectangular mesh point in the surface settlement digital twin mesh model, all entries containing the number of that regular rectangular mesh point are searched in the spatial interpolation mapping table. The corresponding three-dimensional hexahedral mesh node number and weighting coefficient are read from each entry. The surface displacement response components of each three-dimensional hexahedral mesh node are multiplied by their corresponding weighting coefficients and then summed to obtain the predicted surface displacement increment of that regular rectangular mesh point. Perform the above weighted summation operation on all regular rectangular grid points to generate the surface displacement prediction increment field.

[0037] The predicted surface displacement increment is superimposed onto the current surface elevation value of the digital twin grid model for surface subsidence. Each regular rectangular grid point in the digital twin grid model stores a current surface elevation value, initially set to the initial elevation value of the point in one of the embodiments described above. The predicted surface displacement increment corresponding to each regular rectangular grid point is algebraically added to the current surface elevation value according to the rule of positive sign for downward movement, resulting in the updated point elevation value. The updated point elevation values ​​of all regular rectangular grid points constitute the spatial morphology of the updated digital twin grid model for surface subsidence, generating the updated surface subsidence distribution field.

[0038] The residual vector is compared between the updated surface subsidence distribution field and the surface subsidence distribution field generated in the previous iteration. The surface subsidence distribution field generated in the previous iteration is defined as containing the elevation values ​​of all points in the previous iteration's regular rectangular grid. For the first iteration, the elevation values ​​of the points in the previous iteration are taken as the initial elevation values. The residual vector is defined as the vector formed by the differences between the elevation values ​​of each point in the updated surface subsidence distribution field and the corresponding elevation values ​​in the surface subsidence distribution field generated in the previous iteration. The dimension of the residual vector is equal to the total number of regular rectangular grid points. The magnitude of the residual vector is calculated by taking the square root of the sum of the squares of its components, and the unit of magnitude is length. The preset convergence threshold is set according to the required accuracy of subsidence prediction, and is set to one-hundredth of the preset resolution of the regular rectangular grid structure.

[0039] If the magnitude of the residual vector exceeds a preset convergence threshold, the updated surface settlement distribution field is fed back as a new boundary condition to the digital twin mesh model of the tunnel face. Specifically, the feedback method involves comparing the elevation values ​​of each point in the updated surface settlement distribution field with the predicted surface displacement increment transmitted to that point, calculating the stiffness attenuation coefficient of the surrounding rock, and using this coefficient to correct the material parameters of the corresponding elements at the top node in the digital twin mesh model of the tunnel face. This correction includes reducing the elastic modulus. The elastic modulus of the elements in the neighborhood of the top node is multiplied by the stiffness attenuation coefficient, and the overall stiffness matrix is ​​reassembled. The displacement response field is recalculated under the new global stiffness matrix. The above steps are executed iteratively, completing a full cycle from load input, displacement response calculation, surface displacement transfer, surface elevation update to residual comparison and stiffness correction in each iteration. During the iteration process, the current iteration number is recorded; this number is the cumulative value from the start of the bidirectional coupled iterative calculation to the completion of the current iteration. If the number of iterations exceeds the preset maximum number of iterations, the calculation is forcibly terminated, and the current surface settlement distribution field generated in the last iteration is output as the suboptimal prediction result. The preset maximum number of iterations is set based on computational resources and convergence efficiency, and is typically between twenty and fifty.

[0040] The calculation terminates iteratively when the magnitude of the residual vector is less than or equal to a preset convergence threshold. The converged surface subsidence distribution field is output as a surface subsidence prediction cloud map. The surface subsidence prediction cloud map uses a regular rectangular grid structure as its spatial basis. Each grid point is accompanied by its converged elevation value. The subsidence is converted into color information for visualization through chromatographic mapping. A negative subsidence value indicates that the grid point has shifted downwards relative to its initial elevation, while a positive subsidence value indicates uplift. Continuous grid point areas with subsidence exceeding a preset risk threshold in the surface subsidence prediction cloud map are extracted as key subsidence areas. The preset risk threshold is set according to the engineering safety level and the protection standards for surrounding buildings and structures. The extraction process employs a region growth algorithm: grid points with settlement exceeding a preset risk threshold are marked as seed points. The settlement of adjacent grid points is checked using the seed point as the center. If the settlement of an adjacent grid point also exceeds the preset risk threshold, that adjacent grid point is added to the region, and growth continues outwards until the settlement of all adjacent grid points is no greater than the preset risk threshold. This process is repeated until all grid points meeting the criteria are classified into a specific region. The boundary lines of each continuous grid point region constitute the envelope boundary of the critical settlement region identifier. The region within the boundary is the critical settlement region requiring focused attention. The region number, the list of grid points covering the region, and the maximum and average settlement within the region are appended as attribute information to the critical settlement region identifier.

[0041] In practice, the steps for extrapolating the time series of the predicted surface subsidence cloud map are implemented in the following way.

[0042] Settlement trough curves at each monitoring section location are extracted from the surface settlement prediction cloud map. The monitoring section is a vertical plane perpendicular to the tunnel axis. The location of the monitoring section is determined by the mileage value along the tunnel's design axis. The spacing between adjacent monitoring sections is set according to the grid line interval of the regular rectangular grid structure in the surface settlement digital twin grid model along the tunnel axis. For each monitoring section, a row of regular rectangular grid points perpendicular to the tunnel axis is extracted from the surface settlement prediction cloud map corresponding to the section's mileage value. The settlement value of this row of grid points is extracted; the settlement value is the difference between the converged elevation value and the initial elevation value. Downward settlement results in a negative settlement value. A surface settlement distribution curve along the direction perpendicular to the tunnel axis is formed by using the vertical distance from the grid point in this row to the projection line of the tunnel axis on the horizontal plane as the abscissa and the settlement value as the ordinate. This curve is the settlement trough curve.

[0043] The settling trough curve was fitted using Peck's formula, which is expressed as follows: in, This represents the vertical distance from the calculation point to the projection line of the tunnel axis on the horizontal plane, in meters. This indicates the vertical distance from the tunnel axis. The amount of surface subsidence at a given location is expressed in millimeters, with negative values ​​for downward subsidence. This represents the maximum surface settlement above the tunnel axis, in millimeters, and is taken from the settlement trough curve. Settlement at the location; This parameter represents the width of the settling channel, in meters. It characterizes the width characteristic of the settling channel curve, which decreases from the center to both sides. The value is taken as the amount of settlement on the settling channel curve. times The vertical distance corresponding to the location; is the base of the natural logarithm.

[0044] In practice, based on historical surface settlement monitoring sequences from multiple monitoring times, a function is fitted to represent the change of settlement trough curves for each monitoring section as the tunneling ring spacing increases. The historical surface settlement monitoring sequences from multiple monitoring times are obtained from the historical records of the data acquisition module, covering all monitoring data from the start of tunnel excavation to the current time. For each monitoring section, at each historical monitoring time, the real-time settlement displacement values ​​of each regular rectangular grid point at that monitoring section location are extracted from the surface settlement monitoring sequence. The settlement trough curve for that historical monitoring time is constructed in the same manner as described above, and the parameters are fitted using the Peck formula to obtain the settlement trough curve for that historical monitoring time. Value and Value. Record the tunneling ring distance corresponding to this historical monitoring moment. The tunneling ring distance is the cumulative tunneling distance from the starting mileage of tunneling to the mileage of this monitoring section. For the same monitoring section, as the tunneling ring distance increases at different historical monitoring moments, a series of values ​​are obtained. Discrete data pairs of values ​​and tunneling ring spacing, and a series of Discrete data pairs of values ​​and tunneling ring spacing.

[0045] Will Discrete data related to the tunneling ring spacing are fitted to a linear or power function using the least squares method to obtain the settlement channel depth growth coefficient. The settlement channel depth growth coefficient is expressed as a symbol... This means that if a linear function is used... Fit, where, Indicates the tunneling ring spacing. If the intercept is a linear function, then The dimension is the reciprocal of the length dimension, and its value is obtained by solving a system of linear equations using the least squares method. Discrete data related to the tunneling ring spacing were fitted to a linear function using the least squares method to obtain the settlement trough width expansion coefficient. The settlement trough width expansion coefficient is denoted by the symbol... This indicates that a linear function is used. Fit, where, The intercept of the linear function. This is a dimensionless slope coefficient, the value of which is obtained by solving a system of linear equations using the least squares method. The variation function includes the settling tank depth growth coefficient. and the expansion coefficient of the settling trough width The specific form of the change function is the linear expression obtained from the above fitting. After the fitting is completed, the change function parameters of each monitoring section are... , , , Stored in system memory for use in timing extrapolation.

[0046] The tunneling ring spacing values ​​within the future prediction period are sequentially substituted into the variation function. The future prediction period is the time interval from the current moment to the completion of tunnel excavation. The tunneling ring spacing values ​​within the future prediction period are generated according to the daily tunneling progress schedule in the tunnel excavation plan. Specifically, the generation method is to start from the current excavation mileage value at the cutterhead center point and generate the daily predicted tunneling ring spacing value by increasing the daily tunneling ring spacing value until the tunnel's final mileage. For each predicted tunneling ring spacing value corresponding to a future prediction moment, it is substituted into the variation function of each monitoring section to calculate the predicted value for each monitoring section at that prediction moment. Values ​​and Predictions Value, will predict Values ​​and Predictions Substituting the values ​​into the Peck formula, the predicted settlement trough curves for each monitoring section at each future time are calculated. That is, for each monitoring section, a predicted settlement trough curve is generated at each predicted time.

[0047] Extract the time sequence of maximum settlement from each predicted settlement trough curve. The maximum settlement is the predicted value. For each monitoring section, the predicted values ​​for future times are... The values ​​are arranged in chronological order to form a time series, with the horizontal axis representing future times and the vertical axis representing the predicted maximum settlement. The series of maximum settlement changes over time is differentiated to calculate the rate of change of maximum settlement at each time point. The differentiation process uses the central difference method. The rate of change of the maximum settlement at a future moment The calculation method is as follows ,in Indicates the first The maximum settlement at a future time. Indicates the first The time value at a future moment. The unit is millimeters per day. After calculations are completed for all future times, a curve is plotted with the future time as the x-axis and the maximum settlement rate as the y-axis, generating the settlement rate evolution curve for each monitoring section. The settlement rate evolution curve for each monitoring section is output in both graphical and data table formats. The data table contains three columns of data: future time value, predicted maximum settlement value, and maximum settlement rate of change value.

[0048] In practice, the following steps are taken to generate the settlement rate evolution curves and predicted values ​​of settlement threshold exceedance times for each monitoring section, taking into account the spatial distribution characteristics of key settlement area markers.

[0049] Obtain the set of continuous grid points covered by the key settlement area identifier. The key settlement area identifier is generated by the region growing algorithm in the above embodiment, and the output is the envelope boundary of multiple continuous grid point regions and a list of regular rectangular grid point numbers covered by each region. For each key settlement area identifier, read all the regular rectangular grid point numbers covered by the region and store the grid point numbers in the continuous grid point set. In the settlement rate evolution curve data, each regular rectangular grid point corresponds to a settlement rate evolution curve, which is generated by the above embodiment and contains the settlement rate values ​​at future times. Extract the settlement rate evolution curve corresponding to each grid point from the continuous grid point set, and calculate the arithmetic mean of the settlement rate values ​​of all grid points at the same future time using the following formula: in, Indicates the key area in the first place. The average settlement rate at a future time, expressed in millimeters per day; This represents the total number of grid points in a continuous set of grid points; it is dimensionless. Represents the i-th node in the set of continuous grid points The grid point at the th The settlement rate values ​​at future times are calculated in millimeters per day. After performing the above averaging calculation on all future times, a key area average settlement rate curve is generated with the future times as the x-axis and the key area average settlement rate as the y-axis.

[0050] The starting point of the acceleration phase is determined by the difference between the current and previous velocity values ​​on the average settlement rate curve of the key area. The current time is the latest completed monitoring time, and the previous time is the monitoring time immediately preceding the current time. The average settlement rate value corresponding to the current time is extracted from the average settlement rate curve of the key area. The average settlement rate value corresponding to the previous moment , Indicates the current moment. This represents the previous moment. Calculate the difference between the two rate values. If the difference If the value is greater than zero and the difference between three consecutive monitoring times is positive, the settlement rate is determined to have entered the acceleration phase, and the earliest of the three consecutive monitoring times is determined as the starting point of the acceleration phase. Starting from the starting point of the acceleration phase, the average settlement rate curve of the key area is extrapolated and extended. The extrapolation and extension adopts the linear extrapolation method, using the settlement rate values ​​of a fixed number of monitoring times before the extrapolation starting point as fitting samples. The fixed number of monitoring times is no less than five. The least squares method is used to fit a linear regression equation of the settlement rate with respect to time. Based on the slope and intercept of the linear regression equation, the settlement rate values ​​of each extrapolated time in the future are calculated, generating the extended average settlement rate curve segment of the key area. The extended segment is then spliced ​​with the original curve segment to form a complete extrapolated average settlement rate curve of the key area.

[0051] Set a cumulative settlement threshold and a settlement rate threshold. The cumulative settlement threshold is the maximum cumulative settlement limit determined based on the allowable settlement deformation values ​​of surrounding buildings and underground pipelines. The specific value of the cumulative settlement threshold is obtained from the settlement control standards in the engineering design documents, and the unit is millimeters. The settlement rate threshold is the maximum allowable settlement development rate determined according to the construction safety regulations. The specific value of the settlement rate threshold is obtained from the warning rate value in the industry construction specifications, and the unit is millimeters per day. Compare the average settlement rate curve of the key area with the settlement rate threshold. Traverse the average settlement rate curve of the key area along the time axis in a forward direction to extrapolate the settlement rate value at each moment. When the settlement rate value first exceeds the settlement rate threshold, record that moment as the crossing moment of the curve's rate threshold. If the settlement rate value never exceeds the settlement rate threshold after traversing all moments, mark the crossing moment as infinity, indicating that there is no situation where the settlement rate threshold is crossed.

[0052] The predicted settlement values ​​of each grid point in the key settlement area marker are accumulated to the settlement accumulation threshold. In the continuous grid point set covered by the key settlement area marker, a settlement prediction value sequence is established for each grid point. This sequence is extracted from the predicted settlement trough curves generated in the previous embodiment, specifically by extracting the settlement value at each future time point at the corresponding grid point location. For each future time point, the average of the predicted settlement values ​​of all grid points in the continuous grid point set at that time is taken to obtain a sequence of the average settlement of the key area changing over time. Starting from the current average settlement of the key area, the incremental values ​​of the average settlement of the key area at each time point are accumulated sequentially. After each accumulation, the value is compared with the settlement accumulation threshold. When the accumulated value first reaches or exceeds the settlement accumulation threshold, that time is recorded as the accumulation time. The smaller of the crossing time and the cumulative time is used as the predicted value of the settlement threshold exceeding the limit. That is, the earlier of the crossing time and the cumulative time is taken. The predicted value of the settlement threshold exceeding the limit indicates the expected time when the settlement index will first exceed the limit during future construction.

[0053] After generating the settlement rate evolution curves and predicted settlement threshold exceedance times for each monitoring section, the predicted settlement threshold exceedance times are compared with the planned passage times for the corresponding sections in the tunnel excavation plan. The tunnel excavation plan includes the planned passage times for each monitoring section location, which is the planned date and time when the center point of the tunnel boring machine cutterhead reaches the mileage of that monitoring section. For each monitoring section, the predicted settlement threshold exceedance time is subtracted from the planned passage time to obtain the time difference. A positive time difference indicates that the settlement exceedance occurred after tunneling, while a negative time difference indicates that the settlement exceedance occurred before tunneling. The absolute value of the time difference is called the safe construction margin time. If the safe construction margin time is less than the preset safe margin threshold, the safe margin threshold is determined based on the minimum preparation time required for construction emergency response. The safe margin threshold is set to a seven-day construction preparation period. Then, the required adjustment amount of tunneling parameters is calculated by back-calculating the settlement rate evolution curve. The reverse calculation process is as follows: The rate of change of the maximum settlement in the settlement rate evolution curve is used as the target control variable, and the cutterhead thrust and tail grouting pressure are used as adjustment variables. A calibration relationship comparison table is established between the settlement rate change and the cutterhead thrust and tail grouting pressure. This table is obtained by statistically analyzing the settlement rate changes under different combinations of cutterhead thrust and tail grouting pressure from historical construction data. The parameter combinations that can reduce the settlement rate to below the settlement rate threshold are found in the calibration relationship comparison table, and the cutterhead thrust reduction ratio and tail grouting pressure increase ratio are extracted. The cutterhead thrust reduction ratio is the current cutterhead thrust value multiplied by a reduction percentage, with the reduction percentage ranging from 0% to 40%. The minimum reduction ratio that meets the rate control requirements in the calibration relationship comparison table is selected. The tail grouting pressure increase ratio is the current tail grouting pressure value multiplied by an increase percentage, with the increase percentage ranging from 0% to 30%. The minimum increase ratio that meets the rate control requirements in the calibration relationship comparison table is selected. The cutterhead thrust reduction ratio and the tail grouting pressure increase ratio are packaged into a data message for adjusting tunneling parameters and transmitted to the tunnel boring machine (TBM) control system via the industrial Ethernet communication protocol. Upon receiving the adjusted tunneling parameters, the TBM control system adjusts the cutterhead thrust setting by the cutterhead thrust reduction ratio and the tail grouting pressure setting by the tail grouting pressure increase ratio in the subsequent tunneling loop settings. The adjusted construction parameters are then implemented at the start of the next tunneling loop.

[0054] In practical implementation, before inputting the load values ​​of each node in the dynamic load distribution field of the tunnel face as displacement boundary conditions into the equilibrium equation of the digital twin mesh model of the tunnel face, the following material parameter correction steps need to be performed.

[0055] The initial values ​​of elastic modulus and Poisson's ratio corresponding to the lithological parameters of each stratum in the digital twin mesh model of the tunnel face are obtained. In the above embodiment, when constructing the digital twin mesh model of the tunnel face, each mesh node stores the lithological parameters of the stratum to which the node belongs, including the elastic modulus and Poisson's ratio. The initial values ​​of elastic modulus and Poisson's ratio of all nodes are read out to construct an initial material parameter field. The initial material parameter field is an array, the length of which is equal to the total number of nodes in the three-dimensional hexahedral mesh structure. Each element stores the initial values ​​of elastic modulus and Poisson's ratio of the corresponding node.

[0056] Based on the difference between the actual settlement at each monitoring time in the surface subsidence monitoring sequence and the predicted settlement at the corresponding location in the surface subsidence prediction cloud map at the same time, the subsidence prediction error field for each monitoring time is calculated. The actual settlement displacement values ​​for the latest completed preset number of monitoring times are extracted from the surface subsidence monitoring sequence (the preset number is five). For each of the five latest completed monitoring times, the corresponding full-grid point settlement displacement vectors are extracted to form a measured settlement vector matrix. Each column of the measured settlement vector matrix corresponds to a monitoring time, and each row corresponds to a regular rectangular grid point. The predicted settlement displacement values ​​corresponding to the preset number of monitoring times are extracted from the surface subsidence prediction cloud map. The surface subsidence prediction cloud map is the surface subsidence distribution field output after the convergence of the bidirectional coupled iterative calculation in the above embodiment. For each corresponding monitoring time, the predicted settlement displacement values ​​are extracted from the surface subsidence prediction cloud map generated at that time according to the regular rectangular grid point number, forming a predicted settlement vector matrix. The row and column structure of the predicted settlement vector matrix is ​​the same as that of the measured settlement vector matrix. The element-wise difference between the measured settlement vector and the predicted settlement vector is calculated. For each monitoring time and each regular rectangular grid point, the measured settlement displacement value is subtracted from the predicted settlement displacement value to generate the difference vector. Kriging space interpolation is then performed on the difference vector. A spherical model is chosen as the variation function for Kriging space interpolation, with the range of the spherical model being three times the transverse grid spacing of the regular rectangular grid structure. The sill value is taken as the variance of the difference vector, and the nugget value is set to zero. The Kriging space interpolation algorithm is used to interpolate over the entire area of ​​the regular rectangular grid structure, generating a continuous error distribution surface covering the entire regular rectangular grid area. This continuous error distribution surface is used as the settlement prediction error field for the current monitoring time.

[0057] The settlement prediction error field is back-projected onto the node positions of the digital twin mesh model of the tunnel face according to the inverse mapping relationship of the spatial interpolation mapping table. The spatial interpolation mapping table is established by the above embodiment and contains a forward mapping relationship from the 3D hexahedral mesh node number to the regular rectangular mesh point number and weight coefficient. The inverse mapping relationship is to look up all entries in the spatial interpolation mapping table containing the mesh point number for each regular rectangular mesh point number, read the corresponding 3D hexahedral mesh node number and weight coefficient, and for each 3D hexahedral mesh node, the settlement prediction error values ​​of each regular rectangular mesh point are weighted and summed according to the weight coefficient to obtain the back-projected settlement prediction error value of the 3D hexahedral mesh node. The settlement prediction error is back-projected and compared with the initial value of the current elastic modulus of the corresponding node. If the back-projected value is positive, it indicates that the actual settlement is greater than the predicted settlement, and the surrounding rock stiffness is overestimated. In this case, the elastic modulus correction coefficient is less than one, equal to one minus the ratio of the absolute value of the settlement prediction error back-projection to the initial value of the current elastic modulus, and the elastic modulus correction coefficient is limited to between 0.5 and 1. If the settlement prediction error is negative, it indicates that the actual settlement is less than the predicted settlement, and the surrounding rock stiffness is underestimated. In this case, the elastic modulus correction coefficient is greater than one, equal to one plus the ratio of the absolute value of the settlement prediction error back-projection to the initial value of the current elastic modulus, and the elastic modulus correction coefficient is limited to between 1 and 1.5. The initial value of the elastic modulus of each node is multiplied by the corresponding node's elastic modulus correction coefficient to generate the updated elastic modulus value, while Poisson's ratio remains unchanged, resulting in the updated material parameter field. The original initial material parameter field is replaced with the updated material parameter field. After the replacement, when the load values ​​of the dynamic load distribution field at the tunnel face are subsequently input into the equilibrium equations as displacement boundary conditions, the global stiffness matrix of the equilibrium equations is... The element stiffness matrix will be recalculated based on the updated material parameter field, thereby achieving adaptive correction of the material parameters.

[0058] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A surface settlement prediction system for highway tunnel construction based on digital twins, characterized in that, The system includes: The twin modeling module constructs a digital twin mesh model of the tunnel face and a digital twin mesh model of the surface settlement. The digital twin mesh model of the tunnel face and the digital twin mesh model of the surface settlement are dynamically linked through a spatiotemporal coordinate mapping relationship. The data acquisition module collects real-time tunneling parameter sequences and surface settlement monitoring sequences during the tunnel excavation process; The load mapping module maps the real-time tunneling parameter sequence to the corresponding node positions of the digital twin mesh model of the tunnel face, generating a dynamic load distribution field of the tunnel face; The coupled calculation module performs bidirectional coupled iterative calculation processing based on the boundary condition constraint relationship between the dynamic load distribution field of the tunnel face and the digital twin grid model of surface settlement, generating a surface settlement prediction cloud map and key settlement area markers. The temporal extrapolation module performs temporal extrapolation processing on the predicted surface subsidence cloud map, and combines the spatial distribution characteristics of the key subsidence area markers to generate subsidence rate evolution curves and predicted values ​​of subsidence threshold exceedance times for each monitoring section.

2. The surface settlement prediction system for highway tunnel construction based on digital twins according to claim 1, characterized in that, The steps for constructing the digital twin mesh model of the tunnel face and the digital twin mesh model of surface settlement include: Acquire the tunnel design axis and the original surface topographic point cloud data, and generate a three-dimensional hexahedral mesh structure of the tunnel excavation face with the tunnel design axis as the baseline. Each mesh node of the three-dimensional hexahedral mesh structure stores the node coordinates and the lithological parameters of the strata to which the node belongs. A triangular mesh surface model of the land surface is constructed based on the original topographic point cloud data of the land surface. The triangular mesh surface model is then resampled into a regular rectangular mesh structure at a preset resolution. Each grid point of the regular rectangular mesh structure stores the point coordinates and the initial elevation value of the point. A spatial interpolation mapping table is established between the node numbers of a three-dimensional hexahedral mesh structure and the mesh point numbers of a regular rectangular mesh structure. The spatial interpolation mapping table contains the weight coefficient of each three-dimensional hexahedral mesh node to its neighboring regular rectangular mesh points.

3. The surface settlement prediction system for highway tunnel construction based on digital twins according to claim 2, characterized in that, The steps for collecting real-time tunneling parameter sequences and surface settlement monitoring sequences during the tunnel excavation process include: Using the center point of the cutterhead of the tunnel boring machine as the reference origin, the parameters of cutterhead thrust, cutterhead torque and shield tail grouting pressure are collected along the tunnel excavation direction at a fixed excavation ring spacing. Real-time excavation parameter sequence is generated according to the collection order, with each excavation ring spacing corresponding to a parameter recording time. Settlement monitoring sensor arrays are deployed at each grid point of a regular rectangular grid structure on the ground surface. Real-time settlement displacement values ​​of each grid point are collected at equal time intervals to generate a ground settlement monitoring sequence. Each monitoring moment corresponds to a set of settlement displacement vectors of all grid points. The recording time of each parameter is time-aligned with the monitoring time closest to that time to generate synchronized tunneling settlement parameter pairs.

4. The surface settlement prediction system for highway tunnel construction based on digital twins according to claim 3, characterized in that, The step of mapping the real-time tunneling parameter sequence to the corresponding node positions of the digital twin mesh model of the tunnel face to generate the dynamic load distribution field of the tunnel face includes: Based on the tunneling ring distance corresponding to the recording time of each parameter in the real-time tunneling parameter sequence, determine the set of node layers in the digital twin mesh model of the tunnel face that corresponds to the current tunnel face position. The set of node layers includes all three-dimensional hexahedral mesh nodes located within the same tunneling mileage section. The cutterhead thrust parameters are distributed to each node in the node layer set according to the radial distance from the center to the edge of the tunnel face. The proportion of thrust values ​​allocated to nodes in the central area is higher than that allocated to nodes in the edge area. The cutterhead torque parameters are converted into tangential friction load components, and the tangential friction load components are applied to the nodes located at the outer edge of the tunnel face in the node layer set; The tail grouting pressure parameters are converted into normal expansion load components, and the normal expansion load components are applied to the nodes within a fixed ring spacing range behind the node layer assembly.

5. The surface settlement prediction system for highway tunnel construction based on digital twins according to claim 4, characterized in that, When the cutterhead thrust parameters are distributed according to radial distance, a Gaussian function is used to assign weights, with the center node of the tunnel face receiving the largest weight and the edge nodes receiving the smallest weight.

6. The surface settlement prediction system for highway tunnel construction based on digital twins according to claim 4, characterized in that, Based on the boundary condition constraints between the dynamic load distribution field at the tunnel face and the digital twin mesh model of surface settlement, the steps for performing bidirectional coupled iterative calculations to generate a surface settlement prediction cloud map and key settlement area markers include: The load values ​​of each node in the dynamic load distribution field of the tunnel face are used as displacement boundary conditions and input into the equilibrium equation of the digital twin mesh model of the tunnel face to calculate the initial displacement response field of the rock mass surrounding the tunnel. The displacement response components of the surface location are extracted from the initial displacement response field, and the displacement response components are transferred to each grid point of the digital twin grid model of surface subsidence according to the weight coefficients in the spatial interpolation mapping table to generate the predicted surface displacement increment. The predicted surface displacement increment is superimposed onto the current surface elevation value of the digital twin grid model of surface subsidence to update the spatial morphology of the digital twin grid model of surface subsidence and generate an updated surface subsidence distribution field. Compare the residual vector between the updated surface subsidence distribution field and the surface subsidence distribution field generated in the previous iteration. If the magnitude of the residual vector is greater than the preset convergence threshold, the updated surface subsidence distribution field is fed back as a new boundary condition to the digital twin mesh model of the tunnel face, and the displacement response field is recalculated. The above steps are iteratively executed until the magnitude of the residual vector is less than or equal to the preset convergence threshold, at which point the calculation is terminated. The converged surface subsidence distribution field is output as a surface subsidence prediction cloud map, and the continuous grid point regions in the surface subsidence prediction cloud map with subsidence greater than the preset risk threshold are extracted as key subsidence area identifiers.

7. The surface settlement prediction system for highway tunnel construction based on digital twins according to claim 6, characterized in that, When the magnitude of the residual vector is greater than the preset convergence threshold, the current iteration number is recorded. If the iteration number exceeds the preset maximum number of iterations, the calculation is forcibly terminated and the current surface subsidence distribution field is output as the suboptimal prediction result.

8. The surface settlement prediction system for highway tunnel construction based on digital twins according to claim 6, characterized in that, The steps for performing time-series extrapolation processing on the predicted land subsidence cloud map include: The settlement trough curves at each monitoring section location are extracted from the surface settlement prediction cloud map. The settlement trough curves are the surface settlement distribution curves along the direction perpendicular to the tunnel axis. Based on the surface subsidence monitoring sequence at multiple historical monitoring times, a function is fitted to show the change of the subsidence trough curve of each monitoring section with the increase of the tunneling ring distance. The function includes the subsidence trough depth growth coefficient and the subsidence trough width expansion coefficient. Substitute the tunneling ring distance values ​​within the future prediction period into the change function in sequence to calculate the predicted settlement trough curves of each monitoring section at each future time, and extract the sequence of maximum settlement amount changing with time from each predicted settlement trough curve. The sequence of maximum settlement change over time is differentiated to calculate the rate of change of maximum settlement at each moment, generating a settlement rate evolution curve.

9. The surface settlement prediction system for highway tunnel construction based on digital twins according to claim 8, characterized in that, The settling trough curve is fitted using the Peck formula, and the settling trough depth growth coefficient and settling trough width expansion coefficient are determined by fitting historical monitoring data using the least squares method.

10. The surface settlement prediction system for highway tunnel construction based on digital twins according to claim 8, characterized in that, The steps for generating settlement rate evolution curves and predicted settlement threshold exceedance times for each monitoring section, based on the spatial distribution characteristics of the key settlement area markers, include: Obtain the set of continuous grid points covered by the key settlement area marker, calculate the average value of the settlement rate evolution curve of each grid point in the set, and generate the average settlement rate curve of the key area. Based on the difference between the current rate value and the previous rate value in the average settlement rate curve of the key area, the starting point of the acceleration segment of the rate is determined, and the average settlement rate curve of the key area is extrapolated and extended from the starting point of the acceleration segment. Set a threshold for cumulative settlement and a threshold for settlement rate. Compare the average settlement rate curve of the key area with the settlement rate threshold and record the time when the curve crosses the rate threshold. The predicted settlement values ​​of each grid point in the key settlement area are accumulated to the settlement accumulation threshold. The accumulation time when the accumulation threshold is reached is recorded. The smaller of the crossing time and the accumulation time is used as the predicted value of the settlement threshold exceeding the limit.

Images