Modeling method based on shield tunneling data feature analysis and parameter relevance

Abstract

The invention discloses a modeling method based on shield tunneling data feature analysis and parameter relevance, and relates to the field of tunnel engineering data processing. The method comprises the steps that shield tunneling time sequence parameters are obtained, and a non-uniform time sequence is resampled into a space-aligned standardized footage domain sequence through state cleaning and coordinate domain transformation; by means of mixed variable rejection and lagging correlation analysis, environment common cause interference is stripped, physical response delay among parameters is recognized, and a time-delay directed correlation graph model is constructed; and inputting the footage domain sequence and the graph model into a graph neural network, performing feature learning by using a time delay compensation aggregation mechanism, and outputting a key parameter influence degree set with symbols based on a prediction gradient. According to the method, the problem of data space-time dislocation caused by propelling speed fluctuation and the problem of parameter relevance misjudgment caused by physical response lag are solved, and accurate identification and explanation of shield tunneling key parameters are achieved.

Description

Technical Field

[0001] This invention belongs to the field of tunnel engineering data processing, and in particular, it is a modeling method based on the feature analysis and parameter correlation of shield tunneling data. Background Technology

[0002] As the mainstream construction method for urban underground space development, the tunneling process involves the coupling of geotechnical mechanics and mechanical processes. With the improvement of the digitalization level of tunnel boring machines, feature analysis can be performed on multi-source data such as thrust, torque, rotational speed, and geological parameters collected in real time to identify key control parameters and optimize tunneling strategies, thereby ensuring construction safety and improving tunneling efficiency.

[0003] Currently, research on tunnel boring machine (TBM) data analysis primarily relies on time-domain statistical methods and basic deep learning models. Existing techniques typically directly use time-series data, employing statistics such as mean and standard deviation to describe changes in operating conditions, or using Long Short-Term Memory (LSTM) networks for end-to-end regression prediction of time-series data. Regarding parameter correlation analysis, existing methods often use Pearson correlation coefficients or Spearman rank correlation coefficients to calculate a global correlation matrix, thereby selecting characteristic variables and using them as input to neural network models for tunneling performance evaluation.

[0004] However, in complex, unsteady tunneling scenarios, the aforementioned methods face problems of spatiotemporal reference misalignment and physical causality distortion. Therefore, further research and innovation are needed to address these issues in existing technologies. Summary of the Invention

[0005] Purpose of the invention: In view of the above-mentioned problems in the prior art, this application provides a modeling method based on shield tunneling data feature analysis and parameter correlation.

[0006] Technical solution: According to one aspect of this application, a modeling method based on shield tunneling data feature analysis and parameter correlation includes:

[0007] Obtain multi-source shield tunneling time sequence parameters during the shield tunneling process, perform state cleaning and coordinate domain transformation on these parameters, and construct a standardized advance domain sequence.

[0008] Based on the standardized advance domain sequence, using confounding variable elimination and lag correlation analysis, the causal time lag and coupling strength between parameters are identified, and a time lag directed correlation graph model is constructed.

[0009] The standardized advance domain sequence and the time-delay directed correlation graph model are input into the graph neural network model. The time-delay compensation aggregation mechanism in the graph neural network model is used to learn features and obtain the time-delay compensation feature vectors corresponding to each parameter.

[0010] Based on the time-delay compensation feature vector, regression prediction of tunneling performance indicators is performed, and the set of key parameter influence degrees is obtained according to the gradient information and edge weight contribution in the prediction process.

[0011] Beneficial effects: This invention solves the problem of data spatiotemporal misalignment caused by fluctuations in propulsion speed, and the problem of physical causal distortion caused by physical response lag, achieving accurate identification and interpretation of key parameters of tunnel boring machine (TBM) excavation. The related technical effects will be described in detail below with reference to specific embodiments. Attached Figure Description

[0012] Figure 1 A flowchart illustrating a modeling method based on shield tunneling data feature analysis and parameter correlation provided in this application embodiment.

[0013] Figure 2 This application provides a flowchart for performing state cleaning and coordinate domain transformation on shield tunneling time sequence parameters and constructing a standardized advance domain sequence.

[0014] Figure 3 This is a flowchart illustrating how a time delay compensation feature vector is obtained by using a time delay compensation aggregation mechanism in a graph neural network model to perform feature learning and obtain the time delay compensation feature vector corresponding to each parameter, as provided in an embodiment of this application.

[0015] Figure 4 This application provides a flowchart for obtaining a set of key parameter influence degrees based on gradient information and edge weight contributions during the prediction process.

[0016] Figure 5 This is a flowchart illustrating the training process of a graph neural network model provided in an embodiment of this application.

[0017] Figure 6 This is a flowchart of feature processing and optimization based on the ST-GCN-LSTM model provided in the embodiments of this application. Detailed Implementation

[0018] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0019] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0020] To address the aforementioned issues, the applicant conducted in-depth searches and analyses, and discovered:

[0021] Correspondingly, the nonlinear fluctuations in the tunnel boring machine's (TBM) advance speed lead to spatial distortion of the time-domain data. Because the TBM's advance speed changes constantly, the tunneling distance corresponding to the same time window in physical space is not constant. Feature extraction based directly on the time axis causes the parameter waveforms under the same geological section to be compressed or stretched, destroying the geometric comparability of the data.

[0022] Furthermore, static correlation analysis neglects the lag and common-cause interference of physical response. Mechanical transmission and soil response have inherent physical delays; for example, changes in thrust require a certain advance to be reflected in torque, and confounding factors such as soil strata changes can cause multiple parameters to fluctuate simultaneously (pseudo-correlation). Existing techniques cannot isolate these interferences and align causal delays, leading to parameter correlation models that often deviate from the true physical coupling mechanism.

[0023] To solve these problems, combined with Figures 1 to 6 The present invention will be specifically described through the following embodiments.

[0024] On the one hand, an exemplary scheme is provided for a modeling method based on the feature analysis and parameter correlation of shield tunneling data, specifically including:

[0025] Step 101: Obtain the multi-source shield tunneling time sequence parameters during the shield tunneling process, perform state cleaning and coordinate domain transformation on the shield tunneling time sequence parameters, and construct a standardized advance domain sequence.

[0026] In this embodiment, the multi-source shield tunneling timing parameters refer to the raw operational data collected in real time from the shield machine's PLC control system, guidance system, and various additional sensors. Specifically, the parameters include, but are not limited to, the thrust and propulsion speed of the propulsion system, the rotational speed and torque of the cutterhead system, and physical quantities such as the flow rate and pressure of the slurry circulation system. In the raw acquisition state, this data is typically a time series indexed by time.

[0027] State cleaning is primarily used to remove invalid non-tunneling data, such as data from periods of shutdown for maintenance or segment assembly, ensuring that subsequent analysis focuses only on the valid tunneling process. Coordinate domain transformation addresses the time-space misalignment caused by variations in advance speed. Specifically, this step maps the original sequence with time t as the independent variable to a standardized sequence with advance depth u as the independent variable.

[0028] Through this transformation, the data points that were originally unevenly spaced on the time axis due to the different speeds of advancement are resampled onto a uniformly distributed advance grid, achieving physical alignment.

[0029] Step 102: Based on the standardized advance domain sequence, using confounding variable elimination and lag correlation analysis, identify the causal time delay and coupling strength between parameters, and construct a time-delayed directed correlation graph model.

[0030] In this embodiment, the time-delayed directed graph model is a graph structure G=(V, E, A, T) that can express the complex dynamic relationships between shield tunneling parameters. Here, V represents the set of parameter nodes; E represents the set of directed edges; and A is an adjacency matrix whose element a _ij The coupling strength of parameter j to parameter i is represented by τ; T is the time delay matrix, whose element τ _ij This represents the physical hysteresis distance required for a change in parameter j to propagate to parameter i, expressed in advance.

[0031] Specifically, confounding variable removal is used to eliminate spurious correlations caused by common external factors such as formation changes and tool wear. For example, harder formations lead to increased thrust and decreased tunneling speed. Without removing formation factors, directly calculating the correlation between thrust and speed may result in a false strong correlation. Introducing and removing the influence of confounding variables allows for the extraction of purer causal relationships between parameters.

[0032] Building upon this, hysteresis correlation analysis calculates the maximum cross-correlation between parameter pairs in the advance domain, determining the optimal response delay and the corresponding coupling strength. This process upgrades static correlation analysis to causal inference incorporating spatiotemporal dynamic information, providing a physically-compliant graphical prior for subsequent model learning.

[0033] Step 103: Input the standardized advance domain sequence and the time-delay directed correlation graph model into the graph neural network model, and use the time-delay compensation aggregation mechanism in the graph neural network model to perform feature learning, and obtain the time-delay compensation feature vectors corresponding to each parameter.

[0034] In this embodiment, the graph neural network model is a deep learning model specifically designed for processing spatiotemporal graph data. For example, it is a variant of the Spatiotemporal Graph Convolutional Network (ST-GCN), which is also known as a spatiotemporal graph convolutional neural network and is pre-built. A time-delay compensation aggregation mechanism is one of the components of this model. In standard graph convolution, node aggregation typically assumes that the information of neighboring nodes is synchronized.

[0035] However, in tunnel boring machines (TBMs), there is a physical delay between actions (such as increasing thrust) and responses (such as velocity changes). Therefore, when aggregating the features of neighboring nodes, the time-delay compensation aggregation mechanism performs a translation operation on the feature sequence of neighboring nodes in the advance dimension based on the identified optimal response time delay, i.e., it reads u-τ. _ij The features at the current position u are used instead of those at the current position u. This mechanism allows the model to capture the true causal chain, and the generated time-delay compensated feature vector can more accurately reflect the actual contribution of the parameters to the system state.

[0036] Step 104: Based on the time delay compensation feature vector, regression prediction of tunneling performance indicators is performed. According to the gradient information and edge weight contribution in the prediction process, the set of key parameter influence degrees is obtained.

[0037] In this embodiment, the tunneling performance index can be preset, typically selected as a physical quantity that comprehensively reflects tunneling efficiency and safety, such as single-loop penetration. The model learns the nonlinear mapping relationship between the parameters and this index through a regression prediction task. During the prediction process, gradient analysis methods (such as calculating the partial derivative of the output with respect to the input) combined with the edge weights in the graph model can quantify the magnitude and direction of each parameter's contribution to the prediction result.

[0038] The set of key parameter influences is the final output, which not only includes a list of parameters sorted by importance but also clarifies whether each parameter has a positive or negative inhibitory effect. This provides decision-making support for on-site operators, such as identifying whether the key factor currently limiting tunneling efficiency is insufficient torque or excessive soil chamber pressure.

[0039] On the other hand, an alternative implementation of the data standardization method based on ring-level advance domain transformation is described, which solves the data quality problem caused by working condition fluctuations and speed changes during shield tunneling, especially by achieving spatiotemporal alignment of data through state machine and advance domain transformation.

[0040] Step 201: Based on the running status codes recorded by the PLC system, remove the downtime data and assembly time data from the shield tunneling timing parameters, and retain the continuous tunneling section data.

[0041] In this embodiment, the PLC system of the tunnel boring machine records the equipment's operating mode in real time. Typically, status codes clearly distinguish between tunneling mode, assembly mode, and standby mode. The elimination operation specifically involves traversing the original time series and retaining only data segments where the status code corresponds to the tunneling mode and the advance speed is greater than a preset minimum threshold (e.g., 5 mm / min). This step removes invalid zero values ​​or noise data generated by equipment stoppage, sensor drift, or human error, allowing subsequent analysis to focus on the actual tunneling physical process.

[0042] Step 202: Calculate the mean and standard deviation of the continuous tunneling section data using Bessel's formula, and detect and remove outlier data points whose absolute residuals exceed three times the standard deviation based on the 3σ criterion.

[0043] In this embodiment, statistical filtering is employed to further eliminate potential instantaneous spikes or outliers during sensor acquisition. Specifically, for each parameter sequence x, its mean μ and standard deviation σ are calculated over the continuous tunneling section. Bessel's formula is the unbiased estimate of the sample standard deviation: σ = sqrt[(1 / N-1) × ∑ i=1 N (x _i -μ) 2 In other words, Bessel's formula can be expressed as:

[0044] σ=sqrt((1 / (N-1))×∑((x _i -μ) 2 ));

[0045] Where σ is the sample standard deviation; sqrt is the square root operation; N is the total number of samples; x _i is the data value of the i-th sampling point; μ is the arithmetic mean of the data in this segment; ∑ is the summation operation, which accumulates i from 1 to N.

[0046] Furthermore, for each data point x in the sequence _i Calculate its absolute residual |x _i -μ|. If the residual is greater than 3σ, then x is determined to be... _i If an outlier is identified, it is removed and set to null.

[0047] Step 203: Linear interpolation is used to fill and repair the missing data points after removing anomalies, and a smoothing filtering algorithm is used to eliminate high-frequency noise to obtain a standardized advance domain sequence.

[0048] In this embodiment, missing values ​​are filled using linear interpolation with adjacent valid data points to ensure sequence continuity. Furthermore, to eliminate interference from high-frequency random noise in differential calculations (such as velocity and acceleration calculations), algorithms such as moving average filtering or Savitzky-Golay filtering can be used to smooth the data. Here, the standardized advance domain sequence temporarily refers to the cleaned time series, ensuring data quality for subsequent coordinate transformations.

[0049] Step 204: Use multi-parameter joint threshold criteria or PLC status codes to identify the status of shield tunneling timing parameters and extract the ring tunneling data segment under effective propulsion status.

[0050] Accordingly, in addition to relying on PLC status codes, a joint criterion based on physical parameters can also be introduced. Specifically, the system can set the following logical rule: a valid propulsion state is determined when and only when the propulsion speed > speed threshold, total thrust > thrust threshold, and cutterhead rotation speed > rotation speed threshold are simultaneously satisfied. For example, the speed threshold can be set to 10 mm / min, and the thrust threshold to 5000 kN. Through multi-parameter cross-validation, false judgments caused by PLC signal delay or single-point sensor failure can be eliminated. Each segment of data that continuously satisfies this criterion is marked as a ring tunneling data segment, typically corresponding to the tunneling process of one ring of tunnel segments.

[0051] Step 205: Divide the ring tunneling data segment into the starting phase, the stable phase, and the ending phase according to the physical characteristics of the advancement process.

[0052] In this embodiment, considering that the stress state and operation mode of the tunnel boring machine have obvious stage characteristics during the excavation of each ring, directly analyzing the data of the entire ring may obscure local patterns. Therefore, this embodiment introduces a phase segmentation mechanism. Specifically, the advance of each ring can be normalized to the interval [0, 1].

[0053] Based on experience or historical data, the first 20% of the footage, i.e., u∈[0, 0.2], can be defined as the starting phase, which usually involves the gradual loading of thrust and attitude adjustment; the middle 60% of the footage, i.e., u∈(0.2, 0.8], can be defined as the stable phase, where the parameters are relatively stable and best reflect the formation characteristics; and the last 20% of the footage, i.e., u∈(0.8, 1.0], can be defined as the closing phase, which involves unloading and attitude fine-tuning.

[0054] Step 206: For each loop tunneling data segment within a phase, establish a monotonic mapping relationship between the time dimension and the advance footage dimension. Based on the normalized advance footage grid, resample the shield tunneling time series parameters to generate a standardized advance footage domain sequence that eliminates the influence of advance speed fluctuations. This is used to achieve spatiotemporal alignment.

[0055] In the original data, since the advance speed v(t) is variable, the advance increment Δu≈v(t) corresponding to the same time interval Δt is not fixed. To eliminate this non-uniformity, this embodiment constructs a monotonic mapping function u=f(t) from time t to cumulative advance u.

[0056] Furthermore, a fixed advance step size Δu is set. _grid For example, generating a uniform sequence of feed grid points U every 10 millimeters or every 1 centimeter. _grid ={0,Δu _grid ,2Δu _grid , ..., nΔu _grid}, where n is a positive integer, and nΔugrid does not exceed the maximum excavation depth. For each target grid point u _q ∈U _grid Find the two adjacent sampling points (t) in the original data. _a ,u(t _a )) and (t _b ,u(t _b )), making u(t) _a )≤u _q ≤u(t _b The parameter value corresponding to this position is calculated using the linear interpolation formula, specifically:

[0057] x _new (u _q )=x(t _a )+(x(t _b )-x(t _a )) / (u(t _b )-u(t _a ))×(u _q -u(t _a ));

[0058] In other words, the calculation of linear interpolation resampling can be performed using the following formula:

[0059] x _new (u _q )=x(t _a )+((x(t _b )-x(t _a )) / (u(t _b )-u(t _a )))×(u _q -u(t _a ));

[0060] Where, x _new (u _q ) for the target advance grid point u_q The interpolation parameter results at t; _a and t _b The immediate neighbor of u in the original time series _q The two time points before and after; x(t) _a ) and x(t _b () represent time points t _a and t _b The corresponding original parameter value; u(t) _a ) and u(t) _b () represent time points t _a and t _b The corresponding cumulative advance value; u _q This refers to the preset standardized advance grid position.

[0061] By performing the above resampling operation on all parameters, all data can be unified under the same advance coordinate system. The resulting standardized advance domain sequence, indexed by the physical advance, reflects the spatial behavior of the tunnel boring machine and eliminates the problem of data waveform scaling and distortion caused by different operator advance speeds. This lays the foundation for subsequent excavation based on physical location parameter correlations, such as the influence of strata on the cutterhead.

[0062] On the other hand, this paper describes an optional implementation process of a time-delayed directed graph construction method based on confounding elimination and stability constraints, used to solve the spurious correlation problem caused by common cause interference, and the causal misalignment problem caused by time lag. Accordingly, this embodiment includes:

[0063] Step 301: Based on the standardized advance domain sequence, calculate the maximum cross-correlation coefficient between each parameter pair within the advance hysteresis window, and determine the initial response time delay of each parameter pair.

[0064] In other words, the advance lag window is preset and refers to the maximum range of physical response delays searched within the advance domain. Considering the physical characteristics of the tunnel boring machine's mechanical transmission and soil response, this window can be set, for example, to the advance length of two rings at the front and rear, such as ±3 meters. For each pair of tunneling parameters x _i and x _j The system calculates its cross-correlation function R. _ij (τ). This function quantitatively describes the change in parameter x. _j After translation τ on the advance axis, and with parameter x _i The degree of waveform similarity. The specific calculation formula can be expressed as:

[0065] R _ij (τ)=(E[(x _i (u)-μ _i (x) _j (u+τ)-μ _j )]) / (σ_i ×σ _j );

[0066] Where u represents the advance index, τ represents the hysteresis, and μ _i and σ _i The parameters are x respectively _i The mean and standard deviation within the calculation window, E[…] represents the expected value.

[0067] In other words, the cross-correlation function can be described as:

[0068] R _ij (τ)=E[(x _i (u)-μ _i )×(x _j (u+τ)-μ _j )] / (σ _i ×σ _j );

[0069] Among them, R _ij (τ) is the cross-correlation coefficient between parameters i and j with a lag of τ; E represents the expected value; x _i (u) is the value of parameter i at the advance u; x _j (u+τ) is the value of parameter j at the advance u+τ; μ _i and μ _j σ represents the mean values ​​of parameters i and j, respectively; _i and σ _j Let i and j be the standard deviations of parameters i and j, respectively.

[0070] Within a preset window range, for example, τ∈[-300, 300], in cm, the system iterates through all values ​​of τ to find the value that makes the absolute value |R| equal to τ. _ij The lag (τ) that reaches its maximum value is recorded as the initial response time lag. For example, if the fluctuation pattern of the thrust pressure at the current advance u is found to best match the fluctuation pattern of the cutterhead torque at the advance u+50cm, then the time lag between the two is recorded as τ=50cm. That is, the change in thrust is fully transmitted and reflected in the torque after approximately 50cm of tunneling process.

[0071] Step 302: Introduce a set of confounding variables, construct a regression model to eliminate the common influence of the set of confounding variables on the standardized advance domain sequence, and obtain the pure residual sequence corresponding to each parameter.

[0072] In other words, the set of confounding variables can be pre-defined to address spurious correlations caused by common third-party factors. This set C specifically includes the following three categories of physical quantities:

[0073] Stratigraphic label data C reflecting changes in geological conditions _geoFor example, the hardness coefficient of the formation (such as uniaxial compressive strength UCS), water content, or soil type coding;

[0074] Tool wear index C reflects the degradation of equipment cutting performance. _wear This index can be estimated using the cumulative number of cutterhead revolutions, cumulative tunneling distance, or specific energy consumption model.

[0075] Tunneling control mode code C reflecting changes in human operation strategy _mode For example, the status indicator for switching from earth pressure balance mode to semi-open mode can be represented by a 0 / 1 dummy variable.

[0076] Furthermore, to eliminate the common influence of the confounding variable set on the standardized advance range sequence, specifically by using the confounding variable set as the independent variable, performing linear regression on each tunneling parameter, and extracting the residual part after regression as the pure residual sequence.

[0077] Specifically, the elimination operation employs a multiple linear regression method. For each tunneling parameter x in the standardized advance range sequence... _k The following regression equation is established:

[0078] x _k (u)=β _0 +β _1 ×C _geo (u)+β _2 ×C _wear (u)+β _3 ×C _mode (u)+ε _k (u);

[0079] Where β is the regression coefficient, ε _k (u) represents the regression residuals. Or, in other words, β _0 β _1 β _2 β _3 are regression coefficients, representing the degree and direction of the influence of the corresponding variable on the target output value, respectively, and u is the advance.

[0080] After estimating the regression coefficients using the least squares method, the pure residual sequence x is calculated. _tilde_k (u) is ε _k (u). Its calculation formula is:

[0081] x _tilde_k (u)=x _k (u)-(β # _0 +∑ m=1 M β # _m ×C _m (u));

[0082] in, # This represents the _hat symbol.

[0083] In other words, the regression residuals of confounding variables can be expressed as:

[0084] x _tilde_k (u)=x _k (u)-(β _0 +∑(β _m ×C _m (u)));

[0085] Where, x _tilde_k (u) represents the pure residual value of parameter k at the advance u; x _k (u) represents the original observation value; β _0 β is the regression intercept; _m C represents the regression coefficient corresponding to the m-th confounding variable; _m (u) represents the value of the m-th confounding variable at the advance u; ∑ represents the summation over all confounding variables m. β # _0 β # _m For β _0 β _m The estimated value, M is the number of confounding variables.

[0086] Where, x _tilde_k This represents the parameter x after removing the influence of environmental and equipment baselines. _k The system also considers its own fluctuation components. For example, when a tunnel boring machine enters a hard rock section, both thrust and torque will naturally increase due to the hardening of the strata, which is a covariance caused by the strata. Through the above elimination steps, the system removes this covariance component caused by the strata and focuses on analyzing how the active adjustment of thrust causes changes in torque under the same strata conditions, that is, the direct causality between the excavation parameters.

[0087] Step 303: Based on the initial response time delay, calculate the lag correlation between the pure residual sequences corresponding to each parameter, and map the parameter pairs that satisfy the causal significance test and their corresponding initial response time delays to the directed edges and edge attributes of the graph model to obtain the time-delayed directed correlation graph model; where the maximum cross-correlation coefficient and the lag correlation correspond to the coupling strength between each parameter.

[0088] Alternatively, a Granger causality significance test can be performed on the lagged correlation of the pure residual sequence, and a significance level threshold can be set, such as P-value < 0.05, to remove parameter pairs that fail the test.

[0089] Optionally, the standardized advance domain sequence is divided into multiple overlapping advance sliding windows, and the local hysteresis correlation coefficient between the pure residual sequences corresponding to each parameter is calculated within each advance sliding window.

[0090] In this embodiment, windowed analysis is required to capture the dynamic changes in parameter relationships and verify their stability. The system sets the window length to L. _w For example, a 10-ring advance, approximately 15 meters, has a sliding step length of S. _w For example, a 1-ring advance. Within the w-th analysis window, based on the initial response delay, calculate the pure residual sequence x of parameters i and j. _tilde_i and x _tilde_j The correlation coefficient between them is denoted as ρ. ij w The clean data used here provides a more accurate reflection of direct physical coupling. Furthermore, the local correlation coefficient can be calculated using the cross-correlation formula, but only on a subset of data within a local window.

[0091] Optionally, the consistency ratio of the sign direction of the local hysteresis correlation coefficient across all advance sliding windows is statistically analyzed to obtain the parameter correlation stability score.

[0092] In this embodiment, the parameter correlation stability score S _ij The robustness of the correlation between two parameters is a key indicator used in this invention to distinguish between random noise and physical laws. If parameters i and j have a real physical causal relationship, their correlation sign (positive or negative) should remain consistent across most windows. Specifically, the correlation sign (Sign) is calculated globally. _global =sign(ρ ij total ); among which, Sign _global The global correlation symbol is denoted by , and sign(...) is the symbol function, which is the global total correlation value between the i-th object and the j-th object.

[0093] Furthermore, the local correlation coefficient ρ is statistically analyzed across all W windows. ij w The number N windows with the same symbol as the global symbol _same The formula for calculating the parameter correlation stability score is:

[0094] S _ij =(N _same ) / W=(∑ w=1 W I(sign(ρ ij w )==Sign _global )) / W;

[0095] Where I(·) is an indicator function, taking the value 1 when the condition is met and 0 otherwise, and == indicates equality. The score S _ij The value range is from 0 to 1. For example, if a pair of parameters consistently shows a positive correlation across 10 analysis windows, then S... _ij =1 indicates a relatively stable relationship; if 5 windows are positively correlated and 5 windows are negatively correlated, then S _ij =0.5 indicates that the association is very unstable and may be due to random noise or calculation error.

[0096] Alternatively, the stability score can be calculated as follows:

[0097] S _ij =(1 / W)×∑(I(sign(ρ _ij_w )==Sign _global ));

[0098] Among them, S _ij The correlation stability score between parameters i and j is given; W is the total number of sliding windows; ∑ represents the summation over all windows w; I(...) is the indicator function, which takes the value 1 when the condition is met and 0 otherwise; sign(...) is the sign function; ρ _ij_w The local correlation coefficient within the w-th window; Sign _global This is the symbol for the global correlation coefficient.

[0099] Optionally, only parameter pairs whose parameter correlation stability scores exceed the confidence threshold are retained, and directed edges are established from the cause parameter to the result parameter. The parameter correlation stability score and the initial response delay are used together as edge attributes of the directed edges to generate a time-delayed directed correlation graph model.

[0100] Specifically, the system sets a strict confidence threshold S. _th For example, 0.8. Only when S _ij >S _th Only when the time delay τ is reached is a valid edge considered to exist between parameters i and j. The direction of the edge can be determined based on the time delay τ. _ij The sign of τ is determined as follows: if τ _ij If the value of j is greater than 0, meaning that the change of j lags behind that of i, then the physical causal direction is i→j, and a directed edge is established from i to j.

[0101] Based on this, the generated time-delayed directed graph model G not only contains connectivity relationships but also attribute tuples for each edge. Specifically, edge e _i→j The attributes can be represented as:

[0102] Attr _ij =(Weight _ij , τ _ij S_ij );

[0103] Among them, Attr _ij There is a directed edge e _i→j The attribute of Weight _ij The absolute value of the global correlation coefficient, τ, can be taken. _ij For physical time delay, S _ij This is used to score stability. The edge properties mentioned above will guide subsequent feature aggregation with spatiotemporal physical constraints in the graph neural network.

[0104] As an example, an exemplary scheme for a graph neural network including a time delay compensation mechanism and its training method is provided, particularly an improved ST-GCN architecture that addresses the limitation of traditional graph convolutional networks (GCNs) in handling physical latency. Specifically, this scheme includes:

[0105] Step 401: The graph neural network model adopts a cascaded structure of spatial graph convolution and temporal memory network; the model uses spatiotemporal graph convolutional layer (ST-GCN) to extract the spatial coupling features of each tunneling parameter on the parameter association graph topology, and generates a high-dimensional embedding vector sequence containing the interaction information between nodes.

[0106] In other words, the spatiotemporal graph convolutional layer is the basic unit of the spatiotemporal graph convolutional network. Each layer contains sub-modules such as spatial graph convolution, temporal convolution, activation and normalization, which complete a local aggregation and transformation of spatiotemporal features.

[0107] In this embodiment, the cascaded structure is used to simultaneously capture the spatial dependencies (i.e., which parameters influence which) and temporal dependencies (i.e., their own historical evolution) between parameters. Specifically, the input data passes through the ST-GCN layer. This layer uses a time-delayed directed graph model as its skeleton and performs graph convolution operations on the data in each frame's advance slice.

[0108] Unlike CNNs that process images, ST-GCN aggregates data on a non-Euclidean graph structure. After graph convolution, the original scalar parameter values ​​are mapped to high-dimensional hidden feature vectors, containing the interaction information between the parameter and its neighboring parameters at that moment. Furthermore, the spatial feature sequence is fed into subsequent LSTM layers for temporal modeling.

[0109] Step 402: For each target parameter node in the time-delayed directed association graph model, retrieve all its neighbor nodes and the initial response time delay in the corresponding directed edge attributes.

[0110] For any target node in the graph, the algorithm finds all nodes pointing to its source node (neighbor node). Simultaneously, it reads the corresponding initial response delay from the edge attributes; this initial response delay has been converted into a discretized index offset expressed in terms of the mesh advance number. For example, if the mesh advance step is 10mm and the physical delay is 500mm, then the discrete delay = 500 / 10 = 50.

[0111] Step 403: For each neighbor node's corresponding standardized advance domain sequence, perform a translation operation based on the initial response delay in the advance dimension to align the feature sequences of the neighbor nodes with the target parameter node in terms of physical causality.

[0112] In conventional graph convolution, computing the features of node i at advance u usually involves aggregating the features of its neighbor j at the same advance u. However, in this embodiment, the system performs an index shift operation, that is, it reads the features of neighbor node j at advance u-τ. _ij The eigenvector h at that location _j (u-τ _ij If u-τ _ij If the sequence exceeds the boundary, a zero-padding or boundary value replication strategy can be used. If the change in thrust j requires passing through τ... _ij Only a certain advance distance can transmit and cause a change in torque i. Therefore, when analyzing the current state of torque i, τ should be taken into account. _ij The previous state of thrust j, not the current thrust.

[0113] Step 404: Use graph convolution kernels to perform weighted aggregation on the translated and aligned neighbor node sequence and the target parameter node sequence to generate the time delay compensation feature vector corresponding to the target parameter node.

[0114] In this embodiment, the specific aggregation formula can be expressed as:

[0115] h _i_l+1 (u)=σ(W _self ×h _i_l (u)+∑ j∈N(i) A _ij ×W _neigh ×h _j_l (u-τ _ij ));

[0116] Among them, h _i l (u) represents the feature vector of node i at layer l and step u; W _self and W _neigh It is the weight matrix that the model needs to learn during training; σ(...) is the non-linear activation function, such as the rectified linear unit ReLU or the hyperbolic tangent function Tanh; A _ij It is the normalized weight of the edge, h_i_l+1 (u) is the feature vector of node i at layer (l+1) and advance u, h _j_l (u-τ _ij Corresponding to u-τ _ij The node embedding features of the lower source node j in the l-th layer, τ _ij It is the physical time delay from node j to node i, ∑ j∈N(i) This represents the summation operation on the set N(i) of all source nodes j that point to the target node i, i.e., the set of in-neighbors of i. The time-delay compensated feature vector generated in this way integrates effective information from the physical causal chain, avoids feature ambiguity caused by signal misalignment, and achieves the aggregation of correct information at the correct time (footprint).

[0117] On the other hand, a method for calculating time-delay compensated graph convolution aggregation is provided, specifically as follows:

[0118] h _i_new (u)=Activation(W _self ×h _i (u)+∑(A _ij ×W _neigh ×h _j (u-τ _ij )));

[0119] Among them, h _i_new (u) is the updated feature vector of node i at the advance u; Activation is a non-linear activation function, such as ReLU; W _self h is the transformation weight matrix for the node's own features; _i (u) represents the current input feature of node i; ∑ represents the summation over all neighboring nodes j; A _ij W represents the normalized weight of edge ij; _neigh h is the transformation weight matrix for the neighbor features; _j (u-τ _ij ) represents the neighbor node j at the advance position (u-τ) _ij The eigenvector at position τ reflects the physical time delay translation; _ij This is the discretized advance hysteresis obtained from the edge attributes.

[0120] Step 405: Input the high-dimensional embedding vector sequence into the Long Short-Term Memory (LSTM) network layer, use the gating mechanism to capture the long-distance evolution dependence of the tunneling parameters in the advance dimension, and output the time-delay compensated feature vector.

[0121] In this embodiment, the spatial feature sequence extracted by ST-GCN is fed into an LSTM layer. LSTM units, through forget gates, input gates, and output gates, can remember long-term historical information, such as the influence of stratigraphic trends from dozens of rings ago on the present. ST-GCN focuses on the lateral coupling between parameters, while LSTM focuses on the longitudinal evolution of the parameters themselves; the combination of the two constitutes spatiotemporal modeling capabilities. Optionally, gated recurrent units (GRUs) or temporal convolutional networks (TCNs) can be used instead of LSTM layers.

[0122] Step 406 involves calling the training sample set containing historical standardized advance domain sequences and calculating the parameter correlation stability score for each parameter pair in the training sample set. This step is typically performed during the offline training phase.

[0123] Accordingly, a stability score matrix is ​​calculated for the entire dataset using accumulated historical project data. This matrix is ​​fixed as prior knowledge during the model training process to guide the model's weight updates. Specifically, the stability score is a real number between 0 and 1, with values ​​closer to 1 indicating greater stability and reliability of the edge.

[0124] Step 407: Call the loss function that includes a prediction error term and a stability regularization term. The stability regularization term is set as the weighted sum of the absolute values ​​of the edge weights of all directed edges in the time-delayed directed graph model and the instability coefficient. The instability coefficient is determined based on the complementary value of the parameter association stability score.

[0125] Accordingly, to incorporate physical stability constraints into the deep learning optimization process, a dedicated loss function L was designed. This function consists of two parts: a prediction error term L0 and a prediction error term L0. _pred and stability regularization term L _reg The specific formula can be expressed as:

[0126] L=L _pred +λ×L _reg ;

[0127] Among them, L _pred The mean squared error (MSE) can be used for calculation, i.e., L _pred =(1 / N)∑(y _true -y _pred ) 2 L _reg This is a penalty term based on edge weights, and its calculation formula is: L _reg =∑ i，j (1-S _ij )×|A _ij |;In this formula, (1-S _ij This is the instability coefficient, which is used for the stability score S. _ijVery low, meaning edges with a high instability coefficient, if the model assigns them large weights |A during training. _ij |, which will lead to L _reg A significant increase results in penalties. Conversely, for high stability (S... _ij The edge ≈1), (1-S _ij The penalty is small when the value is close to 0, allowing the model to learn its weights freely. λ is a balance coefficient used to control the strength of regularization.

[0128] In other words, the stability regularization loss function can be specifically expressed as the following formula:

[0129] Loss=(1 / N)×∑((y _true -y _pred ) 2 )+λ×∑((1-S _ij )×abs(A _ij ));

[0130] Where Loss is the total loss value; N is the batch size; y _true y is the actual label value; _pred The model predicts the value; λ is the balance coefficient of the regularization term; ∑ represents the summation over all edges i, j; S _ij The stability score for edge ij; abs(A _ij ) represents the edge weights A learned by the model. _ij The absolute value of.

[0131] Step 408: With the goal of minimizing the loss function, the weight parameters of the graph neural network model are iteratively updated using the backpropagation algorithm, so that the model can reduce the prediction error while suppressing the magnitude of the weights of unstable edges.

[0132] In this embodiment, through this physically constrained optimization process, the model is forced to forget accidental and unreliable associations, i.e., S. _ij Instead of focusing on low-level (or high-level) physical laws, we focus on statistically robust physical laws. For example, if a pair of parameters shows a strong correlation only in a few outlier samples, but not in most normal samples (leading to S...), we would consider the following: _ij (Low), conventional training might retain this edge due to overfitting to outliers; however, with the introduction of the regularization term in this embodiment, the model will tend to reduce the weight A of this edge. _ij Compress to 0. This mechanism improves the model's generalization ability under unknown conditions and also enhances the model's interpretability.

[0133] Another example provides a specific implementation process for assessing the impact of key parameters and defining indicators, which includes:

[0134] Step 501: The preset tunneling performance index is single-circle penetration. The single-circle penetration is calculated as follows: based on the standardized advance range sequence, the ratio of the advance speed at the current advance point to the cutterhead rotation speed is calculated.

[0135] In this embodiment, single-circle penetration ratio (PPR) is selected as the tunneling performance index. It represents the distance the tunnel boring machine advances in one rotation of the cutterhead, typically expressed in mm / r. The calculation formula is PPR(u) = v(u) / n(u), where v is the advance speed and n is the cutterhead rotation speed. Compared to simple advance speed, single-circle penetration ratio better reflects the objective hardness of the strata and cutting efficiency, eliminating the influence of operator adjustments to the rotation speed.

[0136] Alternatively, the calculation of the single-loop penetration KPI can be expressed as:

[0137] PPR(u) = v(u) / n(u);

[0138] Where PPR(u) is the single-turn penetration at the feed rate u; v(u) is the feed rate; and n(u) is the cutterhead rotation speed.

[0139] Step 502, wherein when calculating the single-turn penetration, a rotational speed effectiveness threshold is set, data points with a cutter head rotational speed lower than the rotational speed effectiveness threshold are detected and removed, and only the ratio calculation results under stable rotational conditions are retained as the target value for regression prediction.

[0140] In this embodiment, to prevent numerical calculation errors, a rotational speed validity threshold of 0.5 rpm is set. When n(u) < the rotational speed validity threshold, it is determined that the cutter head is at a stop or start-up moment. At this time, calculating v / n will produce an infinite or meaningless value. Therefore, the target value at this point is set to invalid or the value of the previous moment is used. The model only calculates the loss and updates the gradient on the valid data segment.

[0141] Step 503: Calculate the partial derivatives of the tunneling performance index with respect to the time delay compensation eigenvectors corresponding to each parameter, and obtain the gradient vector containing positive and negative directional information.

[0142] In this embodiment, the automatic differentiation function of the deep learning framework is used to calculate the target output y. _pred That is, the predicted PPR, for input feature x _k partial derivative g _k =Яy _pred / Яx _k Where Я corresponds to the partial derivative. The gradient value g _k Reflects parameter x _kSmall changes can lead to an increase or decrease in KPIs. For example, if the thrust gradient is positive, it means that increasing the thrust is beneficial to improving penetration; if the torque gradient is negative, it may mean that excessive torque reflects excessive formation resistance, indicating a decrease in penetration.

[0143] Step 504: Standardize the gradient vector, combine it with the edge weights in the time-delayed directed graph model to perform path aggregation, and calculate the signed comprehensive contribution of each parameter to the tunneling performance index.

[0144] In this embodiment, in order to comprehensively measure the global impact of a parameter, the signed comprehensive contribution C is calculated using the following formula. _k =(1 / N)∑ u ((Яy _pred / Яx _k (u))×((x _k (u)-μ _k ) / σ _k ));

[0145] On the other hand, the formula for calculating the signed comprehensive contribution can also be:

[0146] C _k =(1 / M)×∑((Яy / Яx _k (u))×((x _k (u)-μ _k ) / σ _k ));

[0147] Among them, C _k The signed overall contribution of parameter k; M is the total number of time steps in the evaluation sample, and N corresponds to M; ∑ represents the accumulation of the advance steps u; Яy / Яx _k (u) represents the target output y in relation to the input parameter x. _k The partial derivative (gradient) at the advance u, Я corresponds to the partial derivative; x _k (u) is the input value; μ _k and σ _k These are the mean and standard deviation of parameter k, respectively.

[0148] The multiplication by the standardized input value here is to account for the fluctuation range of the parameters themselves. For features aggregated by graph convolution, a weighted summation can be performed using edge weights to propagate the influence of neighboring nodes back to the source node. In contrast, existing gradient sensitivity methods typically calculate the sum of squared gradients, which loses sign information and cannot distinguish between promotion and inhibition. The method in this embodiment preserves the sign and has stronger physical guidance.

[0149] Step 505: Sort the parameters according to the absolute value of the signed comprehensive contribution, and classify the parameters into positive promoting factors or negative inhibiting factors according to the sign attribute, forming a set of key parameter influence degrees.

[0150] In this embodiment, the system ultimately outputs an ordered list. Each item in the list includes: parameter name, absolute value of contribution (representing importance), and symbol (representing the direction of action). For example: total thrust, contribution 0.85, direction (+), determined as a major promoting factor; cutterhead torque, contribution 0.62, direction (-), determined as a major hindering factor / warning factor; foam injection volume, contribution 0.45, direction (+), determined as a minor promoting factor. This can directly guide on-site operations. For example, when the operator sees that the torque is a negative inhibiting factor, they will understand that the current torque is not a power source, but a manifestation of formation resistance, and will take measures to reduce thrust or increase amendments, rather than blindly increasing torque.

[0151] The three importance assessment indices mentioned above all contribute to the model output by balancing the influence of the feature, and their values ​​are always non-negative. To further identify the direction of the feature's contribution to the prediction target, i.e., to determine whether an increase in the feature promotes or inhibits the prediction output, this invention introduces a signed contribution index. The formula for calculating the signed contribution is:

[0152] C i sign =(1 / N)×∑ n=1 N [(Яy #n / Яx _i n )×(x i_n -μ _i ) / σ _i ];

[0153] Among them, C i sign Represents the signed contribution of the i-th feature; N is the total number of training samples; Яy #n / Яx _i n This represents the model's predicted output y in the nth sample. # For the i-th feature x _i The partial derivative of x reflects the direction of the instantaneous influence of the feature on the output at the current sample point; _i n For the i-th feature of the n-th sample; μ _i and σ _i are the mean and standard deviation of the i-th feature across all samples, respectively.

[0154] When C i signWhen C > 0, it indicates that the feature contributes positively to the overall predicted output; when C i sign When C < 0, it indicates that the feature as a whole has a negative contribution to the predicted output; when C i sign When ≈0, it means that the positive and negative effects of this feature cancel each other out.

[0155] Another example describes an exemplary scheme for a dual-channel mapping method based on the fusion of physical priors and data-driven approaches. This embodiment provides a hybrid modeling approach combining white-box physical mechanisms and black-box data mining, which is suitable for engineering scenarios with high requirements for model interpretability.

[0156] Unlike human skeletal motion recognition, where the graph structure is naturally determined by anatomical and physical connections, shield tunneling parameters, such as cutterhead torque, propulsion speed, excavation chamber pressure, and total thrust, are distributed across different physical subsystems. There is no predefined connection topology between them; their relationships are implicit and dynamic.

[0157] To enable spatiotemporal graph convolutional networks to capture the spatiotemporal relationships between parameters, this invention constructs a shield tunneling parameter graph G=(V, E, A), where V is a set of nodes, each node corresponding to a tunneling parameter; E is a set of edges, representing the relationships between parameters; and A is an adjacency matrix. _ij This represents the edge weight between parameter i and parameter j. This invention employs a dual-channel strategy combining physical priors and data-driven approaches to construct the graph structure.

[0158] The first approach involves constructing a physical priori graph, based on the physical association edges between parameters defining the mechanical transmission relationships and energy flow paths of the tunnel boring machine (TBM) system. Let the TBM tunneling parameter set contain M parameters, forming a node set V = {v...} _1 v _2 , ..., v _M Based on the structure of the tunnel boring machine system, the parameters are divided into the following subsystems: cutterhead system (cutterhead torque, cutterhead speed, penetration depth), propulsion subsystem (total thrust, propulsion speed, hydraulic cylinder pressure in each zone), and soil chamber subsystem (excavation chamber pressure, soil chamber pressure, screw conveyor speed).

[0159] Physical prior adjacency matrix A phys The element definition rules are as follows:

[0160] A phys _ij =1.0, if parameter i and parameter j belong to the same subsystem (strongly related edge);

[0161] A phys _ij =0.5, if parameter i and parameter j have cross-subsystem coupling (weakly related edge);

[0162] A phys _ij =0, if parameter i and parameter j have no physical relationship;

[0163] Strongly correlated edges reflect the direct mechanical transmission relationship within a subsystem, such as the power transmission relationship between cutterhead torque and cutterhead speed; weakly correlated edges reflect the coupling effect between subsystems, such as the influence of cutterhead cutting resistance on propulsion resistance, and the influence of excavation chamber pressure on propulsion speed.

[0164] The second approach involves data-driven graph construction, which analyzes historical tunneling data to mine statistical correlations between parameters. Compared to graph construction methods based on simple correlation coefficients, this method eliminates spurious correlations caused by common causes by removing confounding variables, thus avoiding the introduction of false edges; and it captures the temporal transit characteristics between parameters through lag correlation analysis to determine the directionality and time delay of the correlation.

[0165] Step 601: Based on the mechanical structure and mechanical transmission principle of the shield tunneling system, define the static physical coupling strength between parameter nodes and construct the physical prior adjacency matrix;

[0166] In this embodiment, the physical prior adjacency matrix is ​​not obtained through data calculation, but is predefined by a domain expert knowledge base. The element value 'a' in the matrix... _phy_ij This represents the theoretical physical influence of parameter j on parameter i. For example, according to the physical formula (power = force × velocity), thrust has a direct physical contribution to power, so the corresponding matrix element is set to 1.0. For parameter pairs without a direct physical relationship, such as grouting pressure to cutterhead temperature, their connection weight is set to 0.

[0167] The construction of the physical prior adjacency matrix includes:

[0168] The physical coupling relationships between parameter nodes are divided into internal coupling within the propulsion system, internal coupling within the cutterhead system, and cross-system coupling between propulsion and cutterhead. Preset prior weights are assigned to the three types of coupling relationships.

[0169] Each parameter node is mapped to four physical roles: instruction parameter, execution parameter, status parameter, and indicator parameter. A physical mask matrix is ​​constructed based on physical causal logic.

[0170] The physical mask matrix is ​​used to prohibit the connection of indicator parameters to other parameters, and only retain the one-way connection permission of instruction parameters to execution parameters and execution parameters to status parameters;

[0171] By using the physical mask matrix to perform hard constraint filtering on the prior weights, a valid physical prior adjacency matrix is ​​obtained.

[0172] To standardize the graph structure, physical role constraints are introduced. Specifically, parameters are divided into the following categories: Command parameters (CMD), such as the speed command and thrust command set by the operator; Execution parameters (ACT), such as the actual measured cylinder pressure and motor current; State parameters (STATE), such as soil chamber pressure and attitude deviation; and Key Performance Indicator (KPI), such as tunneling speed and penetration depth.

[0173] The physical mask matrix M is a matrix composed of 0s and 1s. Its generation rules are as follows: if the source node is KPI, then the corresponding row is all 0s; the indicator is the result and cannot be used as a cause. If both the source and target nodes are CMD, then the corresponding positions are 0s; there is no causal relationship between the instructions. Only when connecting causal chains that conform to CMD→ACT→STATE→KPI is the physical mask matrix element M... _ij =1. Based on this, the physical prior adjacency matrix is ​​calculated using the Hadamard product, i.e.:

[0174] A _phy_final =A _phy ×M;

[0175] In other words, physical mask constraints are specifically:

[0176] A _phy_final_ij =A _phy_init_ij ×M _ij ;

[0177] Among them, A _phy_final_ij A _phy_final A represents the elements of the constrained physical prior matrix; _phy_init_ij A _phy The initial physical coupling strength; M, M _ij These are elements of the physical mask matrix, with values ​​of 0 or 1, determined by the role rules.

[0178] Step 602: Calculate the adaptive fusion coefficient that reflects the fluctuation of the current tunneling conditions, and use the adaptive fusion coefficient to perform weighted fusion of the physical prior adjacency matrix and the data-driven adjacency matrix to generate a dual-channel parameter correlation graph model.

[0179] The construction of a data-driven adjacency matrix reflecting dynamic coupling relationships includes:

[0180] Optionally, the Spearman rank correlation coefficients of each parameter pair in the standardized advance domain sequence are calculated within the sliding time window and used as the instantaneous correlation weights at the current moment.

[0181] In this embodiment, the Spearman correlation coefficient is calculated to capture the nonlinear monotonic relationship between parameters. Let the window length be W, and the rank correlation coefficient matrix calculated on the window data at time t is denoted as P._t .

[0182] Optionally, the exponential moving average (EMA) algorithm is used to smoothly update the instantaneous association weights and the weights of the data-driven adjacency matrix at the previous time step to generate the data-driven adjacency matrix at the current time step, thereby eliminating the interference of short-term noise on the stability of the graph structure.

[0183] In this embodiment, a time smoothing mechanism is introduced to prevent drastic changes in the graph structure caused by individual abnormal windows. The specific calculation formula for the exponential moving average is as follows:

[0184] A _data (t)=α×P _t +(1-α)×A _data (t-1);

[0185] Here, α is a smoothing factor, typically ranging from 0.1 to 0.3. For example, if α = 0.2, the data at the current time drives the adjacency matrix A. _data (t) retains 80% of the historical structural information and only incorporates 20% of the current instantaneous correlation. This processing is similar to low-pass filtering, which can obtain a more robust time-varying graph structure.

[0186] Or, in other words, the EMA smoothing formula is as follows:

[0187] A _data_t =α×P _t +(1-α)×A _data_t_minus_1 ;

[0188] Among them, A _data_t The data-driven adjacency matrix at time t; α is the smoothing factor (between 0 and 1); P _t The instantaneous correlation coefficient matrix calculated for time t; A _data_t_minus_1 Let A be the data-driven adjacency matrix of the previous time step t-1. _data (t-1).

[0189] Step 603: Calculate the adaptive fusion coefficient that reflects the fluctuation of the current tunneling conditions, and use the adaptive fusion coefficient to perform weighted fusion of the physical prior adjacency matrix and the data-driven adjacency matrix to generate a dual-channel parameter correlation graph model.

[0190] The calculation of adaptive fusion coefficients and weighted fusion includes:

[0191] Optionally, the coefficient of variation of the standardized advance domain sequence within the current sliding time window is calculated, where the coefficient of variation is the mean of the ratios of the standard deviations to the means of each parameter.

[0192] In this embodiment, the coefficient of variation η(t) is used to quantify the stability of the current tunneling state. For each parameter x within the window... _k Calculate its coefficient of variation cv _k =σ _k / |μ _k The overall operating condition coefficient of variation of the system is the average of the coefficients of variation of all parameters: η(t) = (1 / K)∑cv _k .

[0193] The coefficient of variation for operating conditions can also be calculated using the following formula:

[0194] η(t) = (1 / K) × ∑(σ _k_win / abs(μ _k_win ));

[0195] Where η(t) is the coefficient of variation of the operating conditions at time t; K is the total number of parameters; ∑ represents the summation over all parameters k; σ _k_win σ is the standard deviation of parameter k within the current window. _k ;abs(μ _k_win ) represents the absolute value of the mean of parameter k within the current window, corresponding to |μ _k |

[0196] Optionally, the Sigmoid function is used to map the coefficient of variation of the operating conditions to an adaptive fusion coefficient with a value between 0 and 1; when the coefficient of variation of the operating conditions is less than the preset stability threshold, the adaptive fusion coefficient approaches 0, so that the fusion result is dominated by the physical prior adjacency matrix.

[0197] In this embodiment, the formula for calculating the adaptive fusion coefficient β(t) can be designed as follows:

[0198] β(t)=1 / (1+exp(-k×(η(t)-η _0 )));

[0199] Where, η _0 η is the preset stability threshold (inflection point), and k is the slope parameter that controls the steepness of the transition. When η(t) < η _0 When the operating condition is stable, the exponential term is large, and β(t) approaches 0; when η(t) > η _0 When the operating conditions fluctuate drastically, β(t) approaches 1.

[0200] In other words, the calculation of the adaptive fusion coefficient is as follows:

[0201] β(t)=1 / (1+exp(-k×(η(t)-η _0 )));

[0202] Where β(t) is the adaptive fusion coefficient, taking a value between 0 and 1; exp is the natural exponential function; k is the slope parameter controlling the steepness of the curve; η(t) is the coefficient of variation of the current operating condition; η _0 This is the preset stable threshold inflection point.

[0203] Optionally, the data-driven adjacency matrix is ​​superimposed onto the physical prior adjacency matrix based on the adaptive fusion coefficient to obtain a dual-channel parameter correlation graph model that dynamically adjusts with the stability of the operating conditions.

[0204] In this embodiment, the final graph model adjacency matrix A(t) is calculated as follows:

[0205] A(t) = A _phy +β(t)×A _data (t);

[0206] During steady-state tunneling, when β is close to 0, the model mainly relies on physical priors and ignores small random fluctuations in the data to ensure the robustness of the model. Under non-steady-state or unknown conditions, when β is close to 1, the physical priors may fail, and the model automatically increases the proportion of data-driven weights to capture new coupling relationships using real-time data.

[0207] In other words, performing dual-channel fusion can also be described by the following formula:

[0208] A _final (t)=A _phy_final +β(t)×A _data_t ;

[0209] Among them, A _final (t) is the final merged adjacency matrix; A _phy_final The physical prior adjacency matrix; β(t) is the adaptive fusion coefficient; A _data_t This is a data-driven adjacency matrix.

[0210] In some embodiments, exemplary schemes describing model enhancement, uncertainty analysis, and post-processing provide auxiliary support for the entire modeling method. Specifically, these include:

[0211] Step 701: For any two tunneling parameters in the standardized advance field sequence, calculate the Pearson correlation coefficient, which reflects linear correlation, and the Spearman rank correlation coefficient, which reflects nonlinear monotonic correlation, respectively.

[0212] Construct a linear correlation matrix and a rank correlation matrix containing all parameter pairs, and sort the parameter pairs by correlation strength according to the absolute value of the correlation coefficient;

[0213] Parameter pairs with the highest correlation strength are selected by a predetermined ratio. The selected parameter pairs are mapped to edges of a graph model, and the corresponding correlation coefficients are used as edge weights to construct a parameter correlation graph model.

[0214] Specifically, the system computes two matrices in parallel: the P matrix (Pearson) and the S matrix (Spearman). The Pearson coefficients are calculated using the formula ρ(X,Y)=cov(X,Y) / (σ). _X ×σ _Y The Pearson coefficient is used to capture linear relationships; the Spearman coefficient is the Pearson coefficient calculated for rank sequences. When constructing a graphical model, the maximum absolute value of the two, max(|P_{i=1}^{p}), can be taken. _ij |,|S _ij The edge weight '|' is used to balance linear and nonlinear relationships.

[0215] Above, ρ(X, Y) is the Pearson correlation coefficient between variables X and Y, used to measure the strength and direction of the linear correlation between the two variables; Cov(X, Y) is the covariance between variables X and Y, reflecting the overall error of the two variables and the trend of their co-variance; σ _X σ represents the standard deviation of variable X, characterizing the dispersion of data for variable X; _Y Let P be the standard deviation of variable Y, representing the degree of dispersion of the data for variable Y. _ij |、|S _ij | correspond to the i-th row and j-th column element of the Pearson matrix P and the Spearman matrix S, respectively.

[0216] Step 702: While keeping the Dropout layer of the graph neural network model active, perform multiple forward propagation samplings on the standardized advance domain sequence to obtain multiple sets of feature importance scores.

[0217] In this embodiment, Monte Carlo Dropout (MCDropout) is used to estimate the uncertainty of the model. Typically, Dropout is only enabled during training and disabled during inference. However, in this embodiment, Dropout is forcibly enabled during the inference phase, for example, with a dropout rate p=0.2. The same input data is repeated N times, for example, N=50, in a forward propagation. Because the neurons dropped each time are different, the importance scores of the outputs will also differ each time, forming a score distribution.

[0218] Step 703: Calculate the standard deviation of the importance scores of multiple sets of features as an uncertainty index, and calculate the inverse normalized value of the uncertainty index as the fusion weight.

[0219] In this embodiment, for each parameter, the standard deviation σ of its 50 scores is calculated. _kA larger standard deviation indicates greater uncertainty in the model's assessment of the parameter's importance (i.e., the model hesitates at that point). Define the fusion weight w. _k =1 / (σ _k +ε), where ε is a small quantity that is excluded from zero.

[0220] Alternatively, a formula for calculating the uncertainty fusion weights can be provided, namely:

[0221] w _k =1 / (σ _score_k +ε);

[0222] Among them, w _k The fusion weights for parameter k; σ _score_k ε is the standard deviation of the importance score obtained by parameter k in multiple MCDropout samplings; ε is a small positive number to prevent the denominator from being zero.

[0223] Step 704: Use the fusion weight to perform a weighted average of the feature importance scores of multiple groups to obtain a robust final feature importance score, and generate a set of key parameter influence based on this (final feature importance score).

[0224] In this embodiment, the final score is... _final_k =(∑w _k ×Score _n_k ) / (∑w _k By reducing the weight of samples with high uncertainty, this method can filter out misjudgments caused by random model errors and improve the reliability of the results.

[0225] The robust feature importance score can also be calculated using the following formula:

[0226] Score _final_k =∑(w _k ×Score _n_k ) / ∑(w _k );

[0227] Among them, Score _final_k The final weighted score for parameter k; ∑ represents the summation over multiple samples n; w _k For weight fusion; Score _n_k This represents the original importance score obtained from the nth sampling.

[0228] Step 705: Based on the parameter ranking in the key parameter influence set, select the parameter sequence with high contribution from the standardized advance domain sequence and construct a high-dimensional feature matrix.

[0229] In this embodiment, based on the aforementioned sorting, the first K, for example the first 10, most critical parameters are selected, and sequence data over a period of time are extracted to form a matrix X with dimensions [T, K].

[0230] Step 706: Calculate the covariance matrix of the high-dimensional feature matrix, perform eigenvalue decomposition on the covariance matrix, and determine the variance contribution rate of each principal component.

[0231] In this embodiment, X is subjected to mean removal. The covariance matrix C is calculated as C = (1 / (T-1)) × X. T ×X. Perform eigenvalue decomposition on C to obtain eigenvalues ​​and corresponding eigenvectors. The variance contribution rate of the i-th principal component is r. _i =λ _i / (∑λ _j ), where λ _i , λ _j The i and j eigenvalues ​​are obtained after the eigenvalue decomposition of the covariance matrix C.

[0232] In other words, the PCA covariance matrix can be:

[0233] Cov=(1 / (T-1))×X T ×X;

[0234] Where Cov is the covariance matrix; T is the sample sequence length; X is the high-dimensional feature matrix after removing the mean; X T Let X be the transpose of X.

[0235] Step 707: Select the top K principal components whose cumulative variance contribution rate reaches the preset threshold as projection basis, and map the high-dimensional feature matrix into a low-dimensional comprehensive feature vector as a compressed representation of the shield tunneling state.

[0236] In this embodiment, a cumulative contribution rate threshold is set, for example, 95%. The top k principal components are selected such that ∑ i=1 k r _i ≥0.95. Using k eigenvectors to construct a projection matrix PY, the compressed feature Z = X × PY is calculated. This low-dimensional vector Z retains most of the information from the original key parameter set and can be used for subsequent transmission, storage, or visualization, such as displaying a two-dimensional scatter plot on a large monitoring screen.

[0237] In other embodiments, specific implementation methods for various feature importance evaluation operators are provided, particularly for calculating feature importance scores. In practical applications, one or more of the following basic operators can also be used to generate initial scores, which can then be enhanced by combining them with an uncertainty fusion mechanism.

[0238] Step 801: Calculate feature importance scores based on gradient sensitivity.

[0239] In this embodiment, the gradient sensitivity method focuses on the impact of small changes in the input on the output. For the k-th tunneling parameter x... _k Its importance score I _grad_k The calculation is as follows:

[0240] I _grad_k =(1 / N)∑ u |Яy _pred / Яx _k (u)|;

[0241] It should be understood that this embodiment preferably uses absolute value averaging to avoid the excessive influence of abnormal gradients. This operator mainly reflects the sensitivity of the parameters within a local range.

[0242] Alternatively, the gradient sensitivity can be calculated as follows:

[0243] I _grad_k =(1 / N)×∑(abs(Яy / Яx _k (u)));

[0244] Among them, I _grad_k The gradient sensitivity score for parameter k; N is the number of samples; ∑ represents the summation over all samples u; abs is the absolute value function; Яy / Яx _k (u) is the partial derivative.

[0245] Step 802: Calculate the feature importance score based on the feature perturbation method.

[0246] The feature perturbation method assesses the importance of a feature by disrupting its information. Specifically, keeping other parameters constant, the numerical sequence of the parameters in the test set is randomly shuffled to obtain a perturbed sequence. The perturbed data is then input into the model for prediction, yielding the prediction error E'. If the original prediction error is E, then the importance score I of that parameter is determined. _pert_k =|E'-E| / E; If changing a parameter causes a significant drop in the model's predictive performance, it means that the parameter contains key information; otherwise, it means that the parameter is unimportant.

[0247] In some embodiments, the calculation of the feature perturbation score can also be as follows:

[0248] I _pert_k =abs(Error _perturbed -Error _original ) / Error _original ;

[0249] Among them, I _pert_k Score the importance of feature perturbations for parameter k; Error _perturbedTo shuffle the parameter k, the model's prediction error; Error _original This represents the original prediction error of the model.

[0250] Step 803: Calculate the feature importance score based on the weight contribution.

[0251] In this embodiment, the method directly analyzes the weight parameters inside the graph neural network. For an edge connecting parameter node j to parameter node i, its weight in the first layer of graph convolution is W. _ij_1 Then the global importance score I of parameter j. _weight_j It can be defined as the sum of the magnitudes of the weights of all edges originating from j: I _weight_j =∑ i |W _ij_1 This operator directly reflects the topological connectivity strength learned by the model during training.

[0252] On the other hand, the weighted contribution can also be described as:

[0253] I _weight_j =∑(abs(W _ij_1 ));

[0254] Among them, I _weight_j The weight contribution score of source node j; ∑ represents the summation over all target nodes i; abs(W _ij_1 ) represents the absolute value of the weights connecting j to i in the first layer graph convolution.

[0255] Step 804: Integrate the above operators using an uncertainty fusion mechanism.

[0256] In this embodiment, the above three operators (I) _grad I _pert I _weight This can be considered a source of importance scores for multiple sets of features. The system can calculate three scores separately, use MCDropout to estimate the uncertainty (standard deviation) of each method, and fuse them through inverse variance weighting to obtain a comprehensive importance index that includes both local sensitivity and global robustness.

[0257] In calculating the overall importance score S _j Based on this, further combining the signed contribution C i sign A dual-criteria feature screening method is implemented. The first criterion is importance screening: S features are eliminated. _j Below the preset threshold θ _s The low contribution feature, where θ _s It can be set as the total feature S _j 0.5 times the mean. The second criterion is to filter by contribution direction; for features that pass the first criterion, C is eliminated.i sign <0 and |C i sign |>θ _c The significant negative contribution feature, where θ _c The negative contribution threshold can be set to 0.1.

[0258] Furthermore, define the feature comprehensive evaluation index E. _i for:

[0259] E _i =S _i ×sign(C i sign );

[0260] Among them, S _i The overall importance score for the i-th feature; sign(..) is the sign function, when C i sign When C > 0, the value is increased by 1. i sign When <0, it takes the value -1; when C i sign When =0, it takes the value 0. E _i >0 indicates a positive high contribution characteristic, E _i <0 indicates a negative high contribution characteristic.

[0261] According to one aspect of this application, some methods of the present invention can also be implemented using the following methods:

[0262] Optionally, the shield machine operation data under normal project conditions can be obtained and preprocessed.

[0263] The data collected by the tunnel boring machine (TBM) is acquired through a PLC system or sensors. It contains a large amount of data related to downtime and segment assembly. The downtime and assembly data needs to be removed to obtain operational status data for each segment. This data may contain null values, which must be removed.

[0264] Outliers in tunnel boring machine (TBM) construction data can affect the model's predictive performance. This paper uses the 3σ criterion for outlier removal. For data following a normal distribution, the probability that the data values ​​fall within the range (μ-3σ, μ+3σ) is 0.9974, where σ represents the standard deviation and μ represents the mean. The calculation formula is as follows:

[0265] σ=sqrt[1 / (n-1)×[∑ i=1 n (x _i ) 2 -(1 / n)×(∑ i=1 n x_i ) 2 ]];

[0266] In the formula, n represents the number of samples of shield tunneling construction data, x _i The observation value representing the i-th shield tunneling construction data sample is a specific single construction data record, such as the shield thrust, earth pressure, and other parameter values ​​at a certain moment. i is the sequence number of the sample.

[0267] If the residual error of a key parameter data of a tunnel boring machine meets the condition that the absolute value is greater than 3, then the data is considered abnormal and should be removed.

[0268] Based on this, for the few missing data points caused by rejection or collection failure, linear interpolation or the mean of the previous and next valid data is used to fill in the gaps to ensure the continuity of the data sequence.

[0269] Optionally, the distribution characteristics and dynamic variation characteristics of the main tunneling parameters can be analyzed.

[0270] After data preprocessing, key tunneling parameters such as cutterhead rotation speed, cutterhead torque, penetration depth, total thrust, propulsion speed, and excavation chamber pressure are organized according to ring number or fixed time window, and their mean, standard deviation, interquartile range, and other statistical quantities are calculated to reveal the distribution characteristics of parameters as they change with different propulsion stages and loads.

[0271] Furthermore, the mean, standard deviation, quartiles, and other indicators are calculated for the complete sequence of each major parameter. As needed, this can be extended to higher-order statistical indicators such as skewness and kurtosis. The variation patterns in different stages of advancement can be displayed through line graphs or box plots, which facilitates the identification of the stable range and load level of the parameters.

[0272] Furthermore, based on the fluctuation trend of parameters with ring number, combined with different propulsion stages, load change conditions, tool wear and pressure control changes, potential abnormal patterns or gradual trends are identified, providing an explanatory basis for subsequent correlation and feature importance analysis.

[0273] Furthermore, for abrupt changes, spikes, or abnormal segments discovered in the statistical analysis, the preprocessing results are reviewed and verified to eliminate false anomalies caused by sensor jitter or short-term interference, thus ensuring the reliability of the distribution analysis.

[0274] Optionally, Pearson and Spearman correlation analysis of the tunneling parameters can be performed. To quantify the coupling relationship between the key parameters, the Pearson correlation coefficient and the Spearman rank correlation coefficient are calculated respectively to form a dual-index correlation matrix, which is used to identify the strength of linear and nonlinear correlations between the parameters.

[0275] For each tunneling parameter, Pearson and Spearman coefficients are calculated with all other parameters to form two symmetric matrices that reflect the dependencies between the operational data.

[0276] The correlation coefficients are sorted by absolute value, and the top 30% of highly correlated variables are selected to eliminate weakly correlated and redundant parameters, thereby narrowing the feature space and improving the efficiency of subsequent neural network analysis.

[0277] Optionally, feature importance assessment can be performed based on a neural network model.

[0278] Accordingly, based on the correlation screening, the obtained variable set is input into the constructed neural network model, and the contribution of each feature to the tunneling efficiency or key response indicators is evaluated through the weight changes and gradient distribution during the training process.

[0279] The selected parameters are normalized or standardized to ensure consistent dimensions and avoid interference from physical scale differences in network training. An ST-GCN-LSTM (Spatiotemporal Graph Convolutional Neural Network-Long Short-Term Memory Network) model is constructed, using key tunneling parameters as input and outputs such as tunneling efficiency, penetration depth, or torque stability. The backpropagation algorithm is used to iteratively optimize the network weights.

[0280] Based on this, gradient sensitivity analysis, feature perturbation methods, or methods based on network weight contribution are used to evaluate the strength of the influence of each input parameter on the output index, obtaining normalized feature importance scores. Based on the importance scores, parameters with the highest contribution are selected as the final optimal feature set, providing compact and efficient input variables for PCA dimensionality reduction.

[0281] Optionally, PCA can be used for dimensionality reduction of high-dimensional features. A covariance matrix is ​​constructed for the parameter set filtered by the neural network. Principal component analysis is then used to extract the principal components that best represent the changes in the tunneling state, thus mapping high-dimensional features to low-dimensional comprehensive features. A covariance matrix is ​​constructed for the standardized key parameters, and its eigenvalues ​​and eigenvectors are solved to determine the contribution of each principal component. Based on a cumulative contribution rate threshold (e.g., 98%), the top few principal components are selected to compress the data dimensionality while preserving the main changes in the tunneling state.

[0282] The obtained principal components are used as new low-dimensional input features for subsequent prediction modeling, runtime evaluation, or tunneling strategy optimization, thereby improving model training efficiency and reducing noise interference.

[0283] This application employs a ring-level advance domain transformation technique. By identifying the tunneling section through a state machine and establishing a monotonic mapping between time and advance, the non-uniform time series affected by advance speed fluctuations is resampled into a spatially strictly aligned standardized advance domain sequence, eliminating geometric distortions in the data waveform and enabling subsequent analysis to be based on a unified physical spatial benchmark.

[0284] Furthermore, a time-delay causal relationship graph and compensation model were constructed. On the one hand, the concept of causal inference was introduced, and common interferences from confounding variables such as formation and wear were eliminated through regression, thus removing spurious correlations. On the other hand, the physical response time delays between parameters were calculated, and explicit time-delay index shifting was performed during feature aggregation in the graph neural network. This mechanism enables the model to learn the true physical transmission laws, rather than spurious statistical correlations, achieving accurate quantification and interpretable assessment of the positive and negative impacts of key parameters.

[0285] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A modeling method based on shield tunneling data feature analysis and parameter correlation, characterized in that, The application relates to a shield tunneling performance prediction method based on a time-delay directed association graph model. The method comprises the following steps: acquiring multi-source shield tunneling time sequence parameters in a shield tunneling process, performing state cleaning and coordinate domain transformation on the parameters, and constructing a standardized footage domain sequence; based on the standardized footage domain sequence, using mixed variable elimination and lag correlation analysis to identify the causal time lag and coupling strength between parameters, and constructing a time-delay directed association graph model; inputting the standardized footage domain sequence and the time-delay directed association graph model into a graph neural network model, using the time-delay compensation aggregation mechanism in the graph neural network model to learn features, and obtaining a time-delay compensation feature vector corresponding to each parameter; 2. The method of claim 1, wherein, based on the time-delay compensation feature vector, performing regression prediction on the tunneling performance index, and obtaining a key parameter influence degree set according to gradient information and edge weight contribution in the prediction process. The method comprises the following steps: using a multi-parameter joint threshold criterion or a PLC state code to identify the state of the shield tunneling time sequence parameters, and extracting ring tunneling data segments in an effective advancing state; dividing the ring tunneling data segments into a starting phase, a stable phase and an ending phase according to the physical characteristics of the advancing process; 3. The method according to any one of claims 1 or 2, characterized in that, for the ring tunneling data segments in each phase, a monotonic mapping relationship between the time dimension and the advancing footage dimension is established, the shield tunneling time sequence parameters are resampled based on a normalized footage grid, and a standardized footage domain sequence eliminating the influence of advancing speed fluctuation is generated. using mixed variable elimination and lag correlation analysis to identify the causal time lag and coupling strength between parameters, and constructing a time-delay directed association graph model, comprising: based on the standardized footage domain sequence, the maximum cross-correlation coefficient between each parameter pair is calculated within the footage lag window to determine the initial response time lag of each parameter pair; a set of mixed variables is introduced to construct a regression model to eliminate the common influence of the set of mixed variables on the standardized footage domain sequence, and a pure residual sequence corresponding to each parameter is obtained; based on the initial response time lag, the lag correlation between the pure residual sequences corresponding to each parameter is calculated, the parameter pairs and the corresponding initial response time lag that meet the causal significance test are mapped as the directed edges and edge attributes of the graph model, and the time-delay directed association graph model is obtained; 4. The method of claim 3, wherein, wherein the maximum cross-correlation coefficient and the lag correlation correspond to the coupling strength between parameters. calculating the lag correlation between the pure residual sequences corresponding to each parameter, mapping the parameter pairs and the corresponding initial response time lag that meet the causal significance test as the directed edges and edge attributes of the graph model, comprising: dividing the standardized footage domain sequence into multiple overlapping footage sliding windows, calculating the local lag correlation coefficient between the pure residual sequences corresponding to each parameter in each footage sliding window; statistically analyzing the consistency of the sign direction of the local lag correlation coefficient in all footage sliding windows to obtain a parameter association stability score; only keeping the parameter pairs with the parameter association stability score exceeding a confidence threshold, establishing a directed edge from a cause parameter to a result parameter, and taking the parameter association stability score and the initial response time lag as the edge attributes of the directed edge to generate the time-delay directed association graph model.

5. The method of claim 3, wherein, The graph neural network model comprises at least one time-lag compensation graph convolution layer; feature learning is performed by using a time-lag compensation aggregation mechanism in the graph neural network model to obtain a time-lag compensation feature vector corresponding to each parameter, including: For each target parameter node in the time-lag directed association graph model, all neighbor nodes and initial response time lags in the corresponding directed edge attributes are retrieved; For each neighbor node corresponding to a standardized footage domain sequence, a translation operation based on the initial response time lag is performed in the footage dimension, so that the feature sequence of the neighbor node is physically and causally aligned with the target parameter node; The time-lag compensation feature vector corresponding to the target parameter node is generated by using a graph convolution kernel to perform weighted aggregation on the neighbor node sequence and the target parameter node sequence after translation alignment.

6. The method of claim 1, wherein, According to the gradient information and edge weight contribution in the prediction process, a key parameter influence degree set is obtained, including: The partial derivative of the tunneling performance index with respect to the time-lag compensation feature vector corresponding to each parameter is calculated to obtain a gradient vector containing positive and negative direction information; The gradient vector is standardized and combined with the edge weight in the time-lag directed association graph model to perform path aggregation, and the signed comprehensive contribution degree of each parameter to the tunneling performance index is calculated; According to the absolute value of the signed comprehensive contribution degree, the parameters are sorted, and according to the sign attribute, the parameters are classified as positive promoting factors or negative inhibiting factors to form the key parameter influence degree set.

7. The method of claim 1, wherein, The state cleaning and coordinate domain transformation of the shield tunneling time sequence parameters are performed to construct a standardized footage domain sequence, including: Based on the running state code recorded by the PLC system, the shutdown period data and assembly period data in the shield tunneling time sequence parameters are removed, and the continuous tunneling section data is retained; The mean and standard deviation of the continuous tunneling section data are calculated using the Bezier formula, and the abnormal data points with absolute residual exceeding three times the standard deviation are detected and removed based on the 3σ criterion; The missing data points after removing the abnormal data are filled and repaired using the linear interpolation method, and high-frequency noise is eliminated using the smoothing filter algorithm to obtain the standardized footage domain sequence.

8. The method of claim 4, wherein, The graph neural network model is obtained by training the following steps, specifically including: A training sample set containing historical standardized footage domain sequences is calculated to calculate the parameter correlation stability score of each parameter pair in the training sample set; A loss function containing a prediction error term and a stability regularization term is called, wherein the stability regularization term is set as the weighted sum of the absolute value of the edge weight of all directed edges in the time-lag directed association graph model and the instability coefficient, and the instability coefficient is determined based on the complementary value of the parameter correlation stability score; The weight parameters of the graph neural network model are iteratively updated using the back propagation algorithm to minimize the loss function, so that the model reduces the prediction error while suppressing the numerical size of the low stability edge weight.

9. The method of claim 1, wherein, After obtaining the key parameter influence degree set, further including: According to the key parameter influence degree set, the top-ranked key parameter sequence is selected from the standardized footage domain sequence to construct a high-dimensional feature matrix; The covariance matrix of the high-dimensional feature matrix is calculated, and the eigenvalue decomposition of the covariance matrix is performed to determine the variance contribution rate of each principal component; The first K principal components with a cumulative variance contribution rate reaching a preset threshold are selected as projection bases, and the high-dimensional feature matrix is mapped to a low-dimensional comprehensive feature vector.

10. The method of claim 1, wherein, The generation of the key parameter influence degree set includes: In the case that the Dropout layer of the graph neural network model is activated, the normalized footage domain sequence is sampled multiple times by forward propagation to obtain multiple sets of feature importance scores; The standard deviation of the multiple sets of feature importance scores is calculated as an uncertainty indicator, and the reciprocal normalized value of the uncertainty indicator is calculated as a fusion weight; The multiple sets of feature importance scores are weighted and averaged using the fusion weight to obtain a final feature importance score with robustness, and a key parameter influence degree set is generated accordingly.

Images