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

Abstract

The application discloses a modeling method based on shield tunneling data feature analysis and parameter correlation, and relates to the field of tunnel engineering data processing. The method comprises the following steps: acquiring shield tunneling time sequence parameters, resampling non-uniform time sequences into spatially aligned standardized footage domain sequences through state cleaning and coordinate domain transformation; using mixed variable elimination and lag correlation analysis to strip environmental common cause interference and identify the physical response delay between parameters, and constructing a time-lag directed correlation graph model; inputting the footage domain sequence and the graph model into a graph neural network, using a time-lag compensation aggregation mechanism to learn features, and outputting a signed key parameter influence degree set based on a prediction gradient. The application solves the problems of data space-time dislocation caused by fluctuation of the advancing speed and misjudgment of parameter correlation caused by physical response lag, and realizes accurate identification and interpretation of key parameters of shield tunneling.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of tunnel engineering data processing, and particularly relates to a modeling method based on shield tunneling data feature analysis and parameter correlation. BACKGROUND

[0002] As the mainstream method for developing urban underground space, the shield method involves the coupling of geotechnical mechanics and machinery in the tunneling process. With the improvement of the digital level of shield machines, feature analysis can be performed on multi-source data such as real-time collected thrust, torque, speed and geological parameters, key control parameters can be identified, and the tunneling strategy can be optimized, so as to ensure construction safety and improve tunneling efficiency.

[0003] At present, shield data analysis research mainly relies on time domain statistical methods and basic deep learning models. The existing technology usually directly uses time series data, uses statistical quantities such as mean and standard deviation to describe the working condition change, or uses long short-term memory network (LSTM) to perform end-to-end regression prediction on time series data. In terms of parameter correlation analysis, the existing scheme usually calculates a global correlation matrix using Pearson correlation coefficient or Spearman rank correlation coefficient, and uses the feature variables as inputs of a neural network model to evaluate the tunneling performance.

[0004] However, in a complex non-steady-state tunneling scenario, the above method has problems of time-space benchmark misalignment and physical causal distortion. Therefore, further research and innovation are needed to solve the above problems existing in the prior art. SUMMARY

[0005] The application provides a modeling method based on shield tunneling data feature analysis and parameter correlation.

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

[0007] Obtaining multi-source shield tunneling time series parameters in the shield tunneling process, performing state cleaning and coordinate domain transformation on the shield tunneling time series parameters, and constructing a standardized footage domain sequence;

[0008] Based on the standardized footage domain sequence, using mixed variable elimination and lag correlation analysis, identifying the causal time lag and coupling strength between parameters, and constructing a time lag directed association graph model;

[0009] Inputting the standardized footage domain sequence and the time lag directed association graph model into a graph neural network model, using the time lag compensation aggregation mechanism in the graph neural network model to learn features, and obtaining a time lag compensation feature vector corresponding to each parameter;

[0010] The time delay compensation characteristic vector is used for regression prediction of the tunneling performance index, and a key parameter influence degree set is obtained according to gradient information and edge weight contribution in a prediction process.

[0011] Beneficial effects: The application solves the problems of data space-time dislocation caused by fluctuation of advancing speed and physical causal distortion caused by physical response lag, and realizes accurate identification and interpretation of key parameters of shield tunneling. Related technical effects will be described in detail below in combination with specific embodiments. BRIEF DESCRIPTION OF DRAWINGS

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

[0013] Figure 2 A flowchart of state cleaning and coordinate domain transformation of shield tunneling time sequence parameters and construction of a standardized footage domain sequence provided by the embodiment of the application.

[0014] Figure 3 A flowchart of feature learning by using a time delay compensation aggregation mechanism in a graph neural network model and obtaining time delay compensation characteristic vectors corresponding to each parameter provided by the embodiment of the application.

[0015] Figure 4 A flowchart of obtaining a key parameter influence degree set according to gradient information and edge weight contribution in a prediction process provided by the embodiment of the application.

[0016] Figure 5 A training flowchart of a graph neural network model provided by the embodiment of the application.

[0017] Figure 6 A flowchart of feature processing and optimization based on the ST-GCN-LSTM model provided by the embodiment of the application. DETAILED DESCRIPTION

[0018] In order for those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be described clearly and completely below in combination with the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor should belong to the scope of protection of the present application.

[0019] It is to be understood that the terms first, second, etc. used in the description and the claims of the present application are used to differentiate similar objects, and are not necessarily used to describe a particular sequential or chronological order. It should be understood that the data used in this way can be interchanged under appropriate circumstances, so that the embodiments of the application described herein can be implemented in an order other than that illustrated or described herein. In addition, the terms include and have and any variations thereof are intended to cover non-exclusive inclusion, for example, a process, method, system, product or device including a series of steps or units does not necessarily limit to the clearly listed steps or units, but can include other steps or units not clearly listed or inherent to these processes, methods, products or devices.

[0020] To solve the above problems, the applicant has conducted in-depth search and analysis, and found that:

[0021] Correspondingly, the nonlinear fluctuation of the shield propulsion speed causes the spatial distortion of the time domain data. Since the propulsion speed of the shield machine changes at any time, the corresponding excavation distance in the physical space of the same time window is not constant. Direct feature extraction based on the time axis causes the parameter waveforms under the same geological section to be squeezed or stretched, which destroys the geometric comparability of the data.

[0022] Further, the static correlation analysis ignores the hysteresis and common cause interference of the physical response. There is an inherent physical delay between the mechanical transmission and the soil response, such as the change of the thrust which needs to be reflected on the torque after a certain footage, and the mixed factors such as stratum change will cause multiple parameters to fluctuate at the same time (pseudo correlation). The existing technology cannot strip these interferences and align the causal time lag, so that the parameter correlation model constructed often deviates from the real physical coupling mechanism.

[0023] To solve these problems, in combination with Figures 1 to 6 The present application is specifically illustrated by the following embodiments.

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

[0025] Step 101, acquiring multi-source shield tunneling time sequence parameters in the shield tunneling process, performing state cleaning and coordinate domain transformation on the shield tunneling time sequence parameters, and constructing a standardized footage domain sequence.

[0026] In this embodiment, the multi-source shield tunneling time sequence parameters refer to the original running data collected in real time from the PLC control system, the guide system and various additional sensors of the shield machine. Specifically, the parameters include but are not limited to the thrust of the propulsion system, the propulsion speed, the rotating speed of the cutterhead system, the torque, and the flow and pressure of the slurry circulation system. The above data in the original collection state is usually a time sequence indexed by time.

[0027] Among them, the state cleaning is mainly used to eliminate invalid non-excavation data, such as data during the period of shutdown maintenance, segment assembly, etc., so that the subsequent analysis is only for the effective excavation process. The coordinate domain transformation is used to solve the time-space misalignment problem caused by the change of the advancing speed. Specifically, this step maps the original sequence with time t as the independent variable to the standardized sequence with the advancing footage u as the independent variable.

[0028] Through this transformation, the data points originally unevenly distributed on the time axis due to the fast and slow advancing are resampled to the uniformly distributed footage grid, realizing the alignment in the physical space.

[0029] In step 102, based on the standardized footage domain sequence, the mixed variable elimination and the lag correlation analysis are used to identify the causal time lag and the coupling strength between parameters, and a time-lag directed association graph model is constructed.

[0030] In this embodiment, the time-lag directed association graph model is a graph structure G=(V, E, A, T) that can express the complex dynamic relationship between shield parameters. Among them, V represents the parameter node set; E represents the directed edge set; A is the adjacency matrix, and the element a _ij represents the coupling strength of parameter j to parameter i; T is the time-lag matrix, and the element τ _ij represents the physical lag distance required for the change of parameter j to be transmitted to parameter i, in footage.

[0031] Specifically, the mixed variable elimination is used to eliminate the pseudo-correlation caused by common external factors such as stratum change and tool wear. For example, the hardening of the stratum will cause the increase of the thrust and the decrease of the advancing speed. If the stratum factor is not eliminated, the correlation between the thrust and the speed may be falsely strong. By introducing and eliminating the influence of the mixed variable, the more pure causal relationship between parameters can be extracted.

[0032] On this basis, the lag correlation analysis calculates the maximum cross-correlation between parameter pairs in the footage domain, determines the best response time lag and the corresponding coupling strength. This process upgrades the static correlation analysis to causal inference containing time-space dynamic information, providing a graph structure prior that conforms to the physical law for subsequent model learning.

[0033] In step 103, the standardized footage domain sequence and the time-lag directed association graph model are input into the graph neural network model, and the time-lag compensation aggregation mechanism in the graph neural network model is used for feature learning to obtain the time-lag compensation feature vector corresponding to each parameter pair.

[0034] In this embodiment, the graph neural network model is a deep learning model specially designed for processing spatio-temporal graph data, for example, a variant of the spatio-temporal graph convolutional network ST-GCN, which can also be referred to as a spatio-temporal graph convolutional neural network, which is pre-constructed. The time delay compensation aggregation mechanism is one of the components of the model. In the standard graph convolution, the node aggregation usually assumes that the information of the neighbor nodes is synchronous.

[0035] However, in the shield system, there is a physical delay between the action (such as increasing the thrust) and the response (such as the change in speed). Therefore, when aggregating the features of the neighbor nodes, the time delay compensation aggregation mechanism performs a shift operation on the feature sequence of the neighbor nodes in the footage dimension according to the identified optimal response time delay, that is, reads the features at u-τ _ij instead of the features at the current position u. This mechanism enables the model to capture the real causal chain, and the generated time delay compensation feature vector can more accurately reflect the actual contribution of the parameters to the system state.

[0036] Step 104, regression prediction of the tunneling performance index based on the time delay compensation feature vector, and obtaining a key parameter influence degree set according to the gradient information and edge weight contribution in the prediction process.

[0037] In this embodiment, the tunneling performance index can be pre-set, and is usually selected to be a physical quantity that can comprehensively reflect the tunneling efficiency and safety, for example, the single-circle penetration. The model learns the nonlinear mapping relationship between the parameters and the index through the regression prediction task. In the prediction process, the gradient analysis method (such as calculating the partial derivative of the output with respect to the input) is used in combination with the edge weight in the graph model to quantify the contribution size and direction of each parameter to the prediction result.

[0038] The key parameter influence degree set is the final output result, which not only contains a parameter list sorted by importance, but also explicitly shows whether each parameter plays a positive promoting role or a negative inhibiting role. It provides a decision basis for the on-site operators, for example, identifying whether the key factor restricting the tunneling efficiency is insufficient torque or excessive soil pressure.

[0039] On the other hand, an optional implementation of a data standardization method based on the ring-level footage domain transformation is described, which solves the data quality problem caused by the fluctuation of working conditions and the change of speed in the shield advancing process, especially through the state machine and the footage domain transformation to realize the spatio-temporal alignment of the data.

[0040] Step 201, based on the running state code recorded by the PLC system, the shutdown period data and the assembly period data in the shield tunneling time series parameters are removed, and the continuous tunneling section data is retained.

[0041] In the embodiment, the PLC system of the shield tunneling machine records the working mode of the equipment in real time. Generally, the state code can clearly distinguish the tunneling mode, the assembly mode, the standby mode and the like. The elimination operation is specifically traversing the original time sequence, and only keeping the data segment corresponding to the tunneling mode and the pushing speed greater than the preset minimum threshold (for example, 5 mm / min). This step can remove invalid zero values or noise data generated due to equipment stop, sensor drift or human error, so that the subsequent analysis focuses on the real tunneling physical process.

[0042] In step 202, the mean and standard deviation of the continuous tunneling segment data are calculated by using the Bessel formula, and the abnormal data points with absolute residual exceeding three times the standard deviation are detected and eliminated based on the 3σ criterion.

[0043] In the embodiment, in order to further clean the transient spikes or outliers that may occur in the sensor acquisition process, a statistical method is used for filtering. Specifically, for each parameter sequence x, the mean μ and the standard deviation σ of the continuous tunneling segment are calculated. The Bessel formula is an unbiased estimation formula of the sample standard deviation σ = sqrt[(1 / N-1) × ∑ i=1 N (x _i -μ) 2 ]. Alternatively, the Bessel formula can be expressed as:

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

[0045] Wherein, σ 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 segment data; ∑ is the summation operation, which is accumulated from i to N.

[0046] Further, for each data point x _i in the sequence, the absolute residual |x _i -μ| is calculated. If the residual is greater than 3σ, x _i is determined as an abnormal point and is removed and set to null.

[0047] In step 203, the linear interpolation method is used to fill and repair the missing data points after removing the abnormal points, the high-frequency noise is eliminated by using the smoothing filter algorithm, and the standardized footage domain sequence is obtained.

[0048] In this embodiment, for the generated missing values, linear interpolation is used to fill in the adjacent valid data points before and after, to ensure the continuity of the sequence. Further, in order to eliminate the interference of high-frequency random noise on the differential calculation (such as speed and acceleration calculation), a sliding average filter or Savitzky-Golay filter algorithm can be used for smoothing processing. The standardized footage domain sequence at this stage is temporarily referred to as the cleaned time sequence, which prepares the data quality for subsequent coordinate transformation.

[0049] Step 204, using multi-parameter joint threshold criterion or PLC state code to identify the state of shield tunneling time sequence parameters, and extracting the ring tunneling data segment in the effective pushing state.

[0050] Correspondingly, in addition to relying on the PLC state code, a joint criterion based on physical parameters can also be introduced. Specifically, the system can set the following logical rule: only when the pushing speed > speed threshold, the total thrust > thrust threshold, and the cutterhead speed > speed threshold are met simultaneously, is it determined as an effective pushing state. For example, set the speed threshold to 10 mm / min and the thrust threshold to 5000 kN. Through multi-parameter cross-validation, false positives caused by PLC signal delay or single-point sensor failure can be excluded. Each continuous data segment that meets the criterion is marked as a ring tunneling data segment, which usually corresponds to the tunneling process of a ring segment.

[0051] Step 205, dividing the ring tunneling data segment into start-up phase, stable phase and end phase according to the physical characteristics of the pushing process.

[0052] In this embodiment, considering that the stress state and operation mode of the shield machine have obvious stage characteristics during the tunneling process of each ring, directly analyzing the whole ring data may mask the local rules. Therefore, this embodiment introduces a phase segmentation mechanism. Specifically, the pushing footage of each ring can be normalized to the [0, 1] interval.

[0053] According to experience or historical data statistics, the first 20% of the footage, i.e. u∈[0, 0.2], can be defined as the start-up phase, which usually involves 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, which is relatively stable in parameters and best reflects the stratum characteristics; the last 20% of the footage, i.e. u∈(0.8, 1.0], can be defined as the end phase, which involves unloading and attitude fine-tuning.

[0054] Step 206, for each ring tunneling data segment in each phase, a monotonic mapping relationship between time dimension and pushing footage dimension is established, and the shield tunneling time sequence parameters are resampled based on the normalized footage grid to generate a standardized footage domain sequence that eliminates the influence of pushing speed fluctuations. Used to realize space-time 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 interpolation parameter results at time u _a and t _b are the two time points immediately before and after u _q ; x(t _a ) and x(t _b ) are the original parameter values at time points t _a and t _b ; u(t _a ) and u(t _b ) are the cumulative footage values at time points t _a and t _b ; u _q is the preset normalized footage grid position.

[0061] By performing the above resampling operation on all parameters, all data can be unified under the same footage coordinate system. The normalized footage domain sequence generated thereby has an index of physical footage, reflects the behavior of the shield machine in space, eliminates the problem of stretching and deformation of data waveform caused by different pushing speeds of the operation hand, and lays a foundation for subsequent excavation of parameter correlation based on physical position, such as the influence of strata on the cutter head.

[0062] In yet another aspect, an optional implementation process of a time-lag directed graph construction method based on hybrid rejection and stability constraint is described, which is used to solve the problem of false correlation caused by common cause interference and the problem of causal misplacement caused by time lag. Accordingly, the present embodiment comprises:

[0063] Step 301, based on the normalized footage domain sequence, calculate the maximum cross-correlation coefficient between each pair of parameters within the footage lag window to determine the initial response time lag of each pair of parameters.

[0064] In other words, the footage lag window is preset, which refers to searching for the maximum range of physical response delay on the footage domain. Considering the physical characteristics of mechanical transmission of the shield machine and soil response, the window can be set to, for example, a footage length of 2 rings before and after, 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 degree of similarity between the waveform of parameter x _j and parameter x _i when parameter x _j is shifted by τ on the footage axis. The specific calculation formula can be expressed as:

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

[0066] wherein u denotes the footage index, τ denotes the lag, μ _i and σ _i are the mean and standard deviation of x _i respectively, and E[...] denotes the mathematical expectation.

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

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

[0069] wherein R _ij (τ) is the cross-correlation coefficient of parameter i and parameter j at lag τ; E denotes the mathematical expectation; x _i (u) is the value of parameter i at footage u; x _j (u + τ) is the value of parameter j at footage u + τ; μ _i and μ _j are the mean of parameters i and j respectively; σ _i and σ _j are the standard deviation of parameters i and j respectively.

[0070] The system traverses all τ values within a preset window range, for example, τ ∈ [-300, 300], unit: cm, to find the lag that makes the absolute value |R _ij (τ)| reach the maximum value, and record it as the initial response lag. For example, if it is found that the fluctuation pattern of the thrust pressure at the current footage u is most consistent with the fluctuation pattern of the cutterhead torque at footage u + 50 cm, then record the time lag τ = 50 cm between them. That is, the change of the thrust force is transmitted and reflected in the torque after about 50 cm of tunneling process.

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

[0072] That is, the set of confounding variables can be pre-set to solve the problem of pseudo-correlation caused by third-party common factors. The set C specifically includes the following three types of physical quantities:

[0073] Stratum label data C _geoFor example, the hardness coefficient of the stratum (such as the uniaxial compressive strength UCS), the water content, or the rock type code;

[0074] The tool wear index C reflecting the cutting performance degradation of the device _wear The index can be estimated by the cumulative revolutions of the cutter head, the cumulative tunneling distance, or the specific energy consumption model;

[0075] The tunneling control mode code C reflecting the change of the human operation strategy _mode For example, the state identification of switching from the earth pressure balance mode to the semi-open mode, which can be represented by 0 / 1 dummy variables.

[0076] Further, the common influence of the mixed variable set on the standardized footage domain sequence is removed, specifically, the linear regression of each tunneling parameter is performed using the mixed variable set as the independent variable, and the residual part after the regression is extracted as the pure residual sequence.

[0077] Specifically, the operation of removing uses a multiple linear regression method. For each tunneling parameter x _k in the standardized footage domain sequence, 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, and ε _k (u) is the regression residual. Alternatively, β _0 , β _1 , β _2 , and β _3 are regression coefficients, respectively representing the influence degree and direction of the corresponding variable on the target output value, and u is the footage.

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

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

[0082] wherein, # represents a hat symbol.

[0083] In other words, the mixed variable regression residual can be expressed as:

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

[0085] wherein, x _tilde_k (u) is the pure residual value of the parameter k at the footage u; x _k (u) is the original observation value; β _0 is the regression intercept; β _m is the regression coefficient corresponding to the mth mixed variable; C _m (u) is the value of the mth mixed variable at the footage u; and ∑ represents the summation of all mixed variables m. β # _0 , β # _m are the estimated values of β _0 , β _m , and M is the number of mixed variables.

[0086] wherein, x _tilde_k represents the fluctuation component of the parameter x _k itself after removing the influence of the environment and the device baseline. For example, when the shield machine enters a hard rock section, the thrust and torque will naturally increase due to the hardening of the stratum, which belongs to the covariation caused by the stratum; through the above removal step, the system strips off this covariation component caused by the stratum and focuses on analyzing how the active adjustment of the thrust causes the change of the torque under the same stratum condition, i.e., the direct causality between the excavation parameters.

[0087] In step 303, 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 satisfy the causality significance test are mapped as directed edges and edge attributes of the graph model, and a time-lag directed association graph model is obtained; wherein the maximum cross-correlation coefficient and the lag correlation correspond to the coupling strength between each parameter.

[0088] In other words, the Granger causality significance test is performed on the lag correlation of the pure residual sequence, and the parameter pairs that do not pass the test are removed by setting a significance level threshold, such as P value < 0.05.

[0089] Optionally, the normalized footage domain sequence is divided into a plurality of overlapping footage sliding windows, and in each footage sliding window, a local lag correlation coefficient between the pure residual sequences corresponding to each parameter is calculated respectively.

[0090] In this embodiment, in order to capture the dynamic changes of the parameter relationship and verify its stability, window analysis is needed. The system sets the window length as L _w , for example, 10 ring footage, about 15 meters, and the sliding step as S _w , for example, 1 ring footage. In the wth analysis window, based on the initial response time lag, the correlation coefficient between the pure residual sequences x _tilde_i and x _tilde_j of the parameters i and j is calculated, denoted as p ij w Here, the pure data is used, so the correlation coefficient can better reflect the direct physical coupling. Further, the calculation of the local correlation coefficient can also reuse the cross-correlation formula, but only on the data subset of the local window.

[0091] Optionally, the proportion of the same sign direction of the local lag correlation coefficient in all footage sliding windows is calculated to obtain the parameter correlation stability score.

[0092] In this embodiment, the parameter correlation stability score S _ij is used to measure whether the correlation between two parameters is stable, which is a key indicator of the present application to distinguish between accidental noise and physical law. If there is a real physical cause between the parameters i and j, the correlation sign (positive or negative) of the two should remain consistent in most windows. In specific calculation, the correlation sign Sign _global = sign (p ij total ) in the global range is calculated; wherein Sign _global is the global correlation sign, sign (…) is the sign function, and p ij is the global total correlation value of the ith object and the jth object.

[0093] Further, the number N _same of windows in which the sign of the local correlation coefficient p w is the same as the global sign is counted in all W windows. The calculation formula of the parameter correlation stability score is:

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

[0095] where I(...) is an indicator function that takes 1 when the condition is met and 0 otherwise, and == denotes the equality test. The score S _ij ranges from 0 to 1. For example, if a pair of parameters always shows positive correlation in 10 analysis windows, then S _ij = 1, indicating a relatively stable relationship; if 5 windows show positive correlation and 5 windows show negative correlation, then S _ij = 0.5, indicating that the correlation is very unstable, which might be random noise or calculation error.

[0096] Alternatively, the calculation of the stability score can be:

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

[0098] where S _ij is the stability score of the correlation between parameters i and j; W is the total number of sliding windows; ∑ denotes the summation over all windows w; I(...) is an indicator function that takes 1 when the condition is met and 0 otherwise; sign(...) is the sign function; ρ _ij_w is the local correlation coefficient in the w-th window; Sign _global is the sign of the global correlation coefficient.

[0099] Optionally, only the parameter pairs with a correlation stability score exceeding a confidence threshold are retained, a directed edge from the cause parameter to the result parameter is established, the parameter correlation stability score and the initial response time lag are jointly used as the edge attributes of the directed edge, and a time-lag directed correlation graph model is generated.

[0100] Specifically, the system sets a strict confidence threshold S _th , for example, 0.8. Only when S _ij > S _th , it is considered that there is an effective edge between parameters i and j. For the direction of the edge, the sign of the time lag τ _ij can be used to determine: if τ _ij > 0, i.e., the change of j lags behind i, then the physical causal direction is i→j, and a directed edge from i to j is established.

[0101] On this basis, the generated time-lag directed correlation graph model G not only contains the connection relationship, but also contains the attribute tuple of each edge. Specifically, the attributes of edge e _i→j can be represented as:

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

[0103] where Attr _ij is the attribute of the directed edge e _i→j , Weight _ij is the absolute value of the global correlation coefficient, τ _ij is the physical time delay, and S _ij is the stability score. The above edge attributes will guide the subsequent graph neural network to perform feature aggregation with spatio-temporal physical constraints.

[0104] An example provides an exemplary scheme of a graph neural network comprising a time delay compensation mechanism and a training method thereof, in particular, an improved ST-GCN architecture, which solves the defect that the traditional graph convolutional network GCN cannot handle physical delay. Specifically, the scheme comprises:

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

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

[0107] In this embodiment, the cascaded structure is used to capture the spatial dependence (i.e. who influences whom) and the temporal dependence (i.e. the historical evolution of itself) between parameters at the same time. Specifically, the input data passes through the ST-GCN layer. The layer takes the time delay directed correlation graph model as the skeleton, and performs graph convolution operation on the data on each frame of cutting slice.

[0108] Unlike CNNs that process images, ST-GCN performs aggregation on non-Euclidean graph structures. After graph convolution processing, the original scalar parameter value is mapped to a high-dimensional hidden layer feature vector, which contains the interaction information between the parameter at this moment and its neighbor parameters. Further, the spatial feature sequence will be sent to the subsequent LSTM layer for time dimension modeling.

[0109] Step 402, for each target parameter node in the time delay directed correlation graph model, retrieve all neighbor nodes and the initial response time delay in the corresponding directed edge attribute.

[0110] For any target node in the graph, the algorithm finds all source nodes (neighbor nodes) pointing to it. Meanwhile, the corresponding initial response time lag is read from the edge attribute, which has been converted into a discretized index offset expressed in the number of footage grid. For example, if the footage grid step is 10mm and the physical time lag is 500mm, then the discrete time lag = 500 / 10 = 50.

[0111] At step 403, for each neighbor node corresponding to the 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.

[0112] In conventional graph convolution, the feature of node i at footage u is usually the aggregation of the features of neighbor j at the same footage u. In this embodiment, the system performs an index translation operation, i.e., reads the feature vector h _ij (u-τ _j ) of neighbor j at footage u-τ _ij . If u-τ _ij exceeds the sequence boundary, zero padding or boundary value replication strategy can be used. If the change of thrust j needs to travel a footage distance of τ _ij to cause the change of torque i, when analyzing the current state of torque i, the state of thrust j before τ _ij should be referred to, rather than the current thrust.

[0113] At step 404, the translated and aligned neighbor node sequence and the target parameter node sequence are weighted and aggregated using a graph convolution kernel to generate a time lag compensated 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] where h _i l (u) represents the feature vector of node i at the lth layer and footage u; W _self and W _neigh are weight matrices that need to be learned by the model during training; σ(…) is a nonlinear activation function, such as a rectified linear unit ReLU or a hyperbolic tangent function Tanh; A _ij is the normalized weight of the edge, and h_i_l+1 (u) the feature vector of the corresponding node i at the l+1th layer, at the footage u, h _j_l (u-τ _ij ) the feature vector of the corresponding u-τ _ij the node embedding feature of the lower source node j at the lth layer, τ _ij is the physical time lag of node j to node i, ∑ j∈N(i) denotes the summation operation on all source nodes j pointing to the target node i, i.e., the in-neighbor set N(i). The time-lag compensated feature vector generated in this way fuses the effective information on the physical causal chain, avoids feature ambiguity caused by signal misplacement, and realizes the aggregation of correct information at the correct time (footage).

[0117] In another aspect, a calculation method of time-lag compensated graph convolution aggregation is provided, specifically:

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

[0119] wherein h _i_new (u) is the updated feature vector of node i at the footage u; Activation is a nonlinear activation function, such as ReLU; W _self is the transformation weight matrix of the node's own feature; h _i (u) is the current input feature of node i; ∑ denotes the summation on all neighbor nodes j; A _ij is the normalized weight of edge i-j; W _neigh is the transformation weight matrix of the neighbor feature; h _j (u-τ _ij ) is the feature vector of neighbor node j at the footage position (u-τ _ij ), which reflects the physical time lag translation; τ _ij is the discretized footage lag obtained from the edge attribute.

[0120] Step 405: input the high-dimensional embedding vector sequence into a long short-term memory (LSTM) layer, use the gating mechanism to capture the long-distance evolution dependence of the excavation parameters in the footage dimension, and output the time-lag compensated feature vector.

[0121] In this embodiment, the spatial feature sequence extracted by the ST-GCN is sent to the LSTM layer. The LSTM unit can remember the long-distance historical information, such as the influence of the stratum trend several tens of cycles ago on the current, through the forget gate, the input gate and the output gate. The ST-GCN focuses on the lateral coupling between parameters, while the LSTM focuses on the longitudinal evolution of the parameters, and the combination of the two constitutes the spatio-temporal modeling capability. Alternatively, a gated recurrent unit (GRU) or a temporal convolutional network (TCN) can also be used to replace the LSTM layer.

[0122] Step 406: A training sample set containing a historical standardized footage domain sequence is called to calculate the parameter correlation stability score of each parameter pair in the training sample set. This step is usually performed in the offline training phase.

[0123] Correspondingly, using the accumulated historical project data, the full amount of stability score matrix is calculated. This matrix is fixed as prior knowledge in this model training process and is used to guide the weight update of the model. Specifically, the stability score is a real number between 0 and 1, and the closer the value is to 1, the more stable and reliable the edge is.

[0124] Step 407: 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 correlation graph model and the instability coefficient, and the instability coefficient is determined based on the complementary value of the parameter correlation stability score.

[0125] Correspondingly, in order to integrate the physical stability constraint into the optimization process of deep learning, a special loss function L is designed. The function is composed of two parts, a prediction error term L _pred and a stability regularization term L _reg , and the specific formula can be expressed as:

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

[0127] Wherein, L _pred can be calculated by mean square error (MSE), that is, L _pred =(1 / N)∑(y _true -y _pred ) 2 ; L _reg is a penalty term based on the edge weight, and its calculation formula is: L _reg =∑ i，j (1-S _ij )×|A _ij |; in the formula, (1-S _ij ) is the instability coefficient, and for the stability score S _ijVery low, i.e. the coefficient of instability is high edge, if the model in the training is given a greater weight |A _ij , then it will lead to L _reg large increase, punished. Conversely, for stability is high (S _ij ≈1) edge, (1-S _ij ) close to 0, the penalty is small, allowing the model to learn its weight freely, λ is the balance coefficient, used to control the strength of regularization.

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

[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 sample batch size; y _true is the real label value; y _pred is the model prediction value; λ is the balance coefficient of regularization term; ∑ represents the sum of all edges i, j; S _ij is the stability score of edge i-j; abs(A _ij ) is the absolute value of the edge weight A _ij learned by the model.

[0131] Step 408, using the back propagation algorithm to iteratively update the weight parameters of the graph neural network model, so that the model can reduce the prediction error while suppressing the numerical size of the low stability edge weight.

[0132] In this embodiment, through this optimization process with physical constraints, the model will be forced to forget the accidental and unreliable correlation, i.e. S _ij low), and focus on statistically robust physical laws. For example, if a pair of parameters only shows strong correlation in individual abnormal samples, but is irrelevant in most normal samples (resulting in S _ij low), the conventional training may retain the edge due to overfitting abnormal samples; but after introducing the regularization term of this embodiment, the model will tend to compress the weight A _ij of this edge to 0. This mechanism improves the generalization ability of the model under unknown working conditions, and also enhances the interpretability of the model.

[0133] Another example, provides a specific implementation process of key parameter influence degree evaluation and index definition, accordingly, includes:

[0134] Step 501, the preset tunneling performance index is single-circle penetration; the calculation method of single-circle penetration is: based on the standardized footage domain sequence, the ratio of the pushing speed to the cutter head rotating speed at the current pushing footage is calculated.

[0135] In the embodiment, single-circle penetration PPR is selected as the tunneling performance index. It is the distance of the forward pushing of the shield machine in one rotation of the cutter head, and the unit is usually mm / r. The calculation formula is PPR(u)=v(u) / n(u), wherein v is the pushing speed, and n is the cutter head rotating speed. Compared with the pushing speed alone, the single-circle penetration can better reflect the objective hardness and softness of the stratum and the cutting efficiency, and eliminates the influence of the human adjustment of the rotating speed by the operator.

[0136] In other words, the calculation of the single-circle penetration KPI can be expressed as:

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

[0138] Wherein, PPR(u) is the single-circle penetration at the footage u; v(u) is the pushing speed; and n(u) is the cutter head rotating speed.

[0139] Step 502, wherein when the single-circle penetration is calculated, a rotating speed effectiveness threshold is set, data points with the cutter head rotating speed lower than the rotating speed effectiveness threshold are detected and removed, and only the calculation results of the ratio in the stable rotating state are reserved as the target values for the regression prediction.

[0140] In the embodiment, in order to prevent numerical calculation errors, the rotating speed effectiveness threshold is set to 0.5 rpm. When n(u)<rotating speed effectiveness threshold, it is determined that the cutter head is in the state of stopping or starting instant, at this time, the calculation of v / n will produce infinite or meaningless values, therefore the target value of the point is set to invalid or the value of the last time is used. The model only calculates the loss and updates the gradient on the effective data segment.

[0141] Step 503, the partial derivative of the tunneling performance index with respect to the time lag compensation feature vector corresponding to each parameter is calculated, and a gradient vector containing positive and negative direction information is obtained.

[0142] In the embodiment, the automatic differentiation function of the deep learning framework is used to calculate the partial derivative g _pred of the predicted PPR y _k with respect to the input feature x _k , that is, g _pred = Яy _k / Яx _k , wherein Я corresponds to the partial derivative. The gradient value g _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, according to the absolute value of the signed contribution degree, the parameters are sorted, and the parameters are classified as positive promotion factors or negative inhibition factors according to the sign attribute, forming a set of key parameter influence degrees.

[0150] In this embodiment, the system finally outputs an ordered list. Each item in the list contains: parameter name, contribution degree absolute value (representing importance), sign (representing direction of action). For example: total thrust, contribution degree 0.85, direction (+), determined as a major promotion factor; cutter torque, contribution degree 0.62, direction (-), determined as a major obstacle factor / warning factor; foam injection amount, contribution degree 0.45, direction (+), determined as a secondary promotion factor. It can directly guide the field operation, for example, when the operator sees that the torque is a negative inhibition factor, he will understand that the current torque is not a power source, but a manifestation of formation resistance, and will take measures to reduce the thrust or increase the modifier, rather than blindly increase the torque.

[0151] The above three importance evaluation indexes all measure the influence strength of the feature on the model output, and their values are always non-negative. To further identify the contribution direction of the feature to the prediction target, that is, to judge whether the increase of the feature promotes or inhibits the prediction output, the invention introduces a signed contribution degree index. The calculation formula of the signed contribution degree is:

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

[0153] Wherein, C i sign represents the signed contribution degree of the i-th feature; N is the total number of training samples; y #n / y _i n represents the partial derivative of the model prediction output y # to the i-th feature x _i at the n-th sample, reflecting the instantaneous influence direction of the feature on the output at the current sample point; x _i n is the value of the i-th feature of the n-th sample; μ _i and σ _i are the mean and standard deviation of the i-th feature in all samples, respectively.

[0154] When C i sign> 0 indicates that the feature overall contributes positively to the predicted output; when C i sign < 0 indicates that the feature overall contributes negatively to the predicted output; when C i sign ≈ 0 indicates that the positive and negative effects of the feature cancel each other out.

[0155] In another example, an example scheme of a dual-channel composition method based on physical priori and data-driven fusion is described, and the embodiment provides a hybrid modeling idea of white-box physical mechanism + black-box data mining, which is suitable for engineering scenes with high requirements for model interpretability.

[0156] Unlike the human skeleton action recognition scene where the graph structure is naturally determined by anatomical physical connection, shield tunneling parameters such as cutter torque, pushing speed, excavation bin pressure, total thrust, etc. are distributed in different physical subsystems and do not have a predefined connection topology between them. The correlation between them is implicit and dynamic.

[0157] To enable the spatio-temporal graph convolutional network to capture the spatio-temporal correlation between parameters, the present application constructs a shield tunneling parameter graph G=(V, E, A), wherein V is a node set, each node corresponding to a tunneling parameter; E is an edge set, representing the correlation between parameters; A is an adjacency matrix, A _ij represents the edge weight between parameter i and parameter j. The present application adopts a dual-channel strategy combining physical priori and data-driven to construct the graph structure.

[0158] The first channel is physical priori graph construction, which defines the physical correlation edges between parameters based on the mechanical transmission relationship and energy flow path of the shield machine system. Let the shield machine tunneling parameter set contain M parameters, which constitute the node set V={v _1 , v _2 ,..., v _M}. According to the structure of the shield machine system, the parameters are divided into the following subsystems: cutterhead subsystem (cutterhead torque, cutterhead speed, penetration), propulsion subsystem (total thrust, pushing speed, cylinder pressure in each zone), soil bin subsystem (excavation bin pressure, soil bin pressure, screw speed).

[0159] The element definition rule of the physical priori adjacency matrix A phys is as follows:

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

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

[0162] A phys _ij =0, if parameter i is not physically related to parameter j;

[0163] wherein, the strong correlation edge reflects the direct mechanical transmission relationship within the subsystem, such as the power transmission relationship between the cutter torque and the cutter rotating speed; the weak correlation edge reflects the coupling effect between the subsystems, such as the influence of the cutter cutting resistance on the propulsion resistance, and the influence of the excavation chamber pressure on the propulsion speed.

[0164] The second channel is data-driven graph construction, which analyzes the statistical correlation between parameters by analyzing historical tunneling data. Compared with the graph construction method based on simple correlation coefficient, this method eliminates false correlation caused by common causes by mixed variable elimination, avoiding the introduction of pseudo edges; through lag correlation analysis, the time sequence transmission characteristics between parameters are captured, and the directionality and time delay of the correlation are determined.

[0165] Step 601, based on the mechanical structure and mechanical transmission principle of the shield tunneling system, the static physical coupling strength between the parameter nodes is defined, and the physical prior adjacency matrix is constructed;

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

[0167] wherein, the physical prior adjacency matrix is constructed, including:

[0168] The physical coupling relationship between the parameter nodes is divided into propulsion system internal coupling, cutter system internal coupling, and propulsion-cutter cross-system coupling, and preset prior weights are assigned to the three types of coupling relationships respectively;

[0169] Map each parameter node to four physical roles of instruction parameter, execution parameter, state parameter and index parameter, and construct a physical mask matrix according to the physical causal logic;

[0170] The physical mask matrix is used to prohibit the connection of index parameters to other parameters, and only retains the one-way connection permission of instruction parameters to execution parameters and execution parameters to state parameters;

[0171] The physical mask matrix is used to prohibit the connection of index parameters to other parameters, and only retains the one-way connection permission of instruction parameters to execution parameters and execution parameters to state parameters;

[0172] To regulate the graph structure, physical role constraints are introduced. Specifically, parameters are divided into the following categories: command parameters (CMD), such as the rotational speed command and the thrust command set by the operator; execution parameters (ACT), such as the actual measured cylinder pressure and motor current; state parameters (STATE), such as the soil bin pressure and attitude deviation; and index parameters (KPI), such as the tunneling speed and penetration depth.

[0173] The physical mask matrix M is a matrix composed of 0 and 1. Its generation rule is: if the source node is KPI, the corresponding row is all 0, and the index is the result, which cannot be used as the cause in reverse; if the source node is CMD and the target node is CMD, the corresponding position is 0, and there is no causality between the commands; only when the connection meets the causal chain of CMD→ACT→STATE→KPI, the physical mask matrix element M _ij =1. On this basis, the physical prior adjacency matrix is calculated by Hadamard product, that is:

[0174] A _phy_final =A _phy ×M;

[0175] In other words, the physical mask constraint is specifically:

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

[0177] Wherein, A _phy_final_ij , A _phy_final are the physical prior matrix elements after constraint; A _phy_init_ij , A _phy are the physical coupling strengths defined initially; M, M _ij are the physical mask matrix elements, which take values of 0 or 1 and are determined by the role rules.

[0178] Step 602, an adaptive fusion coefficient reflecting the fluctuation degree of the current tunneling working condition is calculated, and the adaptive fusion coefficient is used to weight and fuse the physical prior adjacency matrix and the data-driven adjacency matrix to generate a dual-channel parameter correlation graph model.

[0179] Wherein, the data-driven adjacency matrix reflecting the dynamic coupling relationship is constructed, including:

[0180] Optionally, the Spearman rank correlation coefficient of each parameter pair in the standardized footage domain sequence is calculated in the sliding time window as the instantaneous correlation weight at the current time.

[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 be P_t .

[0182] Optionally, the instantaneous correlation weight is updated by an exponential moving average (EMA) algorithm and the data-driven adjacency matrix weight of the previous time is smoothed to generate the data-driven adjacency matrix of the current time, so as to eliminate the interference of short-time noise on the stability of the graph structure.

[0183] In the embodiment, in order to prevent the graph structure from jumping sharply due to individual abnormal windows, a time smoothing mechanism is introduced. The specific calculation formula of the exponential moving average is as follows:

[0184] A _data (t) = a x P _t + (1-a) x A _data (t-1);

[0185] Wherein, a is a smoothing factor, and the value range is usually 0.1 to 0.3. For example, a = 0.2, the data-driven adjacency matrix A _data (t) retains 80% of the historical structure information and only fuses 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, the EMA smoothing formula is as follows:

[0187] A _data_t = a x P _t + (1-a) x A _data_t_minus_1 ;

[0188] Wherein, A _data_t is the data-driven adjacency matrix at time t; a is a smoothing factor (between 0 and 1); P _t is the instantaneous correlation coefficient matrix calculated at time t; A _data_t_minus_1 is the data-driven adjacency matrix at the previous time t-1, corresponding to A _data (t-1).

[0189] Step 603, calculating an adaptive fusion coefficient reflecting the fluctuation degree of the current tunneling working condition, and using the adaptive fusion coefficient to weight and fuse the physical prior adjacency matrix and the data-driven adjacency matrix to generate a dual-channel parameter correlation graph model.

[0190] Wherein, the adaptive fusion coefficient is calculated and weighted and fused, including:

[0191] Optionally, the working condition variation coefficient of the standardized footage domain sequence in the current sliding time window is calculated, and the working condition variation coefficient is the mean value of the ratio of the standard deviation to the mean value of each parameter.

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

[0193] The working condition variation coefficient can also be solved by the following formula:

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

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

[0196] Optionally, the working condition variation coefficient is mapped to an adaptive fusion coefficient with a value between 0 and 1 using a Sigmoid function; when the working condition variation coefficient is less than a preset stable threshold, the adaptive fusion coefficient tends to 0, so that the fusion result is dominated by the physical prior adjacency matrix.

[0197] In the embodiment, the calculation formula of the adaptive fusion coefficient β(t) can be designed as:

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

[0199] where η _0 is a preset stable threshold (inflection point), and k is a slope parameter that controls the steepness of the transition. When η(t)<η _0 , i.e., the working condition is stable, the exponential term is large, and β(t) tends to 0; when η(t)>η _0 , i.e., the working condition is severely fluctuating, β(t) tends to 1.

[0200] Alternatively, the calculation of the adaptive fusion coefficient is specifically:

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

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

[0203] Optionally, the data-driven adjacency matrix is superimposed to the physical prior adjacency matrix based on the adaptive fusion coefficient, to obtain a dual-channel parameter correlation graph model dynamically adjusted according to the operating condition stability.

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

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

[0206] In the steady-state tunneling, that is, β is close to 0, the model mainly relies on the physical prior, and the small random fluctuations of the data are ignored, to ensure the robustness of the model; in the non-steady-state or unknown operating condition, that is, β is close to 1, the physical prior may fail, the model automatically increases the proportion of the data-driven weight, and the new coupling relationship is captured by using the real-time data.

[0207] In other words, the dual-channel fusion is performed, which can also be described as the following formula:

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

[0209] Wherein, A _final (t) is the final fused adjacency matrix; A _phy_final is the physical prior adjacency matrix; β(t) is the adaptive fusion coefficient; and A _data_t is the data-driven adjacency matrix.

[0210] In some embodiments, an exemplary scheme for describing model enhancement, uncertainty analysis and post-processing is provided to assist the entire modeling method. Specifically, it includes:

[0211] Step 701, for any two tunneling parameters in the standardized footage domain sequence, respectively calculating the Pearson correlation coefficient reflecting linear correlation and the Spearman rank correlation coefficient reflecting nonlinear monotonic correlation;

[0212] Constructing a linear correlation matrix and a rank correlation matrix containing all parameter pairs, and sorting the parameter pairs according to the absolute value of the correlation coefficient;

[0213] The parameter pairs with a correlation strength ranking in a preset proportion are screened out, the screened parameter pairs are mapped as edges of the graph model, and the corresponding correlation coefficients are taken as edge weights, so as to construct the parameter correlation graph model.

[0214] Specifically, the system simultaneously and in parallel calculates two matrices: a P matrix (Pearson) and an S matrix (Spearman). The Pearson coefficient calculation formula is ρ(X, Y) = cov(X, Y) / (σ _X Xσ _Y ), which is used to capture linear relationships; the Spearman coefficient is a Pearson coefficient calculated on rank order sequences. When constructing the graph model, the maximum absolute value max(|P _ij |, |S _ij |) of the two can be taken as the edge weight, taking into account both linear and nonlinear relationships.

[0215] Above, ρ(X, Y) is the Pearson correlation coefficient of variable X and variable Y, which is used to measure the strength and direction of the linear correlation between the two variables; Cov(X, Y) is the covariance of variable X and variable Y, which reflects the overall error and the trend of cooperative change of the two variables; σ _X is the standard deviation of variable X, which represents the dispersion degree of the data of variable X itself; σ _Y is the standard deviation of variable Y, which represents the dispersion degree of the data of variable Y itself. |P _ij |, |S _ij | correspond to the i-th row and j-th column elements of the Pearson matrix P and the Spearman matrix S, respectively.

[0216] Step 702, under the condition that the Dropout layer of the graph neural network model is activated, the standardized footage domain sequence is sampled multiple times in forward propagation to obtain multiple groups of feature importance scores.

[0217] In this embodiment, Monte Carlo Dropout (MC Dropout) technology is used to estimate the uncertainty of the model. Generally, Dropout is only enabled during training and disabled during inference. However, in this embodiment, Dropout is forcibly enabled during the inference stage, for example, the dropout rate p = 0.2. For the same input data, N times, for example, N = 50, forward propagation is repeated. Since the neurons discarded each time are different, the importance scores output each time will also differ, forming a score distribution.

[0218] Step 703, the standard deviation of the multiple groups 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.

[0219] In this embodiment, for each parameter, the standard deviation σ _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, according to the aforementioned ranking, the top K, for example, the top 10, most critical parameters are selected, the sequence data in a period of time is extracted to form a matrix X with a dimension of [T, K].

[0230] In step 706, a covariance matrix of the high-dimensional feature matrix is calculated, and eigenvalue decomposition is performed on the covariance matrix to determine the variance contribution rate of each principal component.

[0231] In this embodiment, X is subjected to mean removal processing. The covariance matrix C=(1 / (T-1))×X T ×X is calculated. Eigenvalue decomposition is performed 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 represent the i th and j th eigenvalues obtained after 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] Wherein, Cov is the covariance matrix; T is the sample sequence length; X is the high-dimensional feature matrix after mean removal; X T is the transpose matrix of X.

[0235] In step 707, the top K principal components whose cumulative variance contribution rate reaches a preset threshold are selected as a projection base, and the high-dimensional feature matrix is mapped to a low-dimensional comprehensive feature vector as a compressed representation of the shield tunneling state.

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

[0237] In some other embodiments, specific implementation methods of a plurality of feature importance evaluation operators are provided, in particular, the calculation of feature importance scores. In actual applications, one or more of the following basic operators can be used to generate initial scores, and an uncertainty fusion mechanism can be used for enhancement.

[0238] In step 801, the feature importance score is calculated based on the gradient sensitivity.

[0239] In this embodiment, the gradient sensitivity method focuses on the impact of small changes in input on the output. For the kth tunneling parameter x _k , its importance score I _grad_k is calculated as:

[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 excessive impact of abnormal gradients. This operator mainly reflects the sensitivity of the parameter in the local range.

[0242] Alternatively, the calculation of gradient sensitivity can be:

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

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

[0245] Step 802, based on the feature perturbation method, calculate the feature importance score.

[0246] Where the feature perturbation method assesses the importance of a feature by destroying its information. Specifically, while keeping other parameters unchanged, the numerical sequence of the parameter in the test set is randomly shuffled to obtain the perturbed sequence. The perturbed data is input into the model for prediction to obtain the prediction error E'. If the original prediction error is E, then the importance score I _pert_k of the parameter = |E'-E| / E; if shuffling one parameter leads to a significant decline in model prediction performance, it means that the parameter contains key information; otherwise, it means that the parameter is not important.

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

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

[0249] where I _pert_k is the feature perturbation importance score of parameter k; Error _perturbedTo break the prediction error of the model after the parameter k; Error _original For 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 the edge connecting parameter node j to parameter node i, its weight in the first layer graph convolution is W _ij_1 . Then the global importance score I _weight_j of parameter j can be defined as the sum of the weight modules of all edges from j: I _weight_j =∑ i |W _ij_1 | This operator directly reflects the topological connection strength learned by the model during training.

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

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

[0254] Where I _weight_j is the weight contribution score of source node j; ∑ represents the summation of all target nodes i; abs(W _ij_1 ) is the absolute value of the weight 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 ) can be regarded as the source of multiple sets of feature importance scores. The system can calculate the three scores respectively, estimate the uncertainty (standard deviation) of each method using MCDropout, and fuse them by inverse variance weighting to obtain a comprehensive importance index that contains both local sensitivity and global robustness.

[0257] Based on the calculation of the comprehensive importance score S _j , further combine the signed contribution C i sign Implement double-criteria feature screening. The first criterion is importance screening: remove low-contribution features whose S _j is lower than the preset threshold θ _s , where θ _s can be set to 0.5 times the mean value of all features S _j . The second criterion is contribution direction screening. For features that pass the first criterion, remove Ci 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 construction data, x _i represents the observation value of the i-th shield construction data sample, is a specific single construction data record, such as the parameter value of the shield thrust force, earth pressure, etc. at a certain time, and i is the serial number identifier of the sample.

[0267] If the residual error of a certain shield tunneling key parameter data satisfies the absolute value > 3, the data is considered to be abnormal data and should be removed.

[0268] On this basis, for a small amount of data missing points caused by removal or failure to collect, linear interpolation method or mean filling method of effective data before and after is used to repair, to ensure the continuity of the data sequence.

[0269] Optionally, the distribution characteristics and dynamic change characteristics of the main tunneling parameters are analyzed.

[0270] After completing the data preprocessing, the cutterhead speed, cutterhead torque, penetration, total thrust, advancing speed, excavation chamber pressure and other key tunneling parameters are arranged according to the ring number or fixed time window, and the mean, standard deviation, quartile distance and other statistical quantities are calculated to reveal the distribution characteristics of the parameters with different advancing stages and load changes.

[0271] Further, the mean, standard deviation, quartile degree and other indexes of the complete sequence of each main parameter are calculated, which can be extended to skewness, kurtosis and other high-order statistical indexes as needed, and the change law in different advancing stages is displayed through the broken line chart or box plot, which is convenient for identifying the stable interval of the parameter and the load level.

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

[0273] Further, for the mutation points, peaks or abnormal sections found in statistical analysis, the preprocessing results are verified to remove pseudo-abnormal points caused by sensor jitter or short-term interference, and the reliability of the distribution analysis is ensured.

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

[0275] For each tunneling parameter, the Pearson and Spearman coefficients are calculated with all other parameters respectively, forming two symmetric matrices reflecting the dependence relationship between the operation data.

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

[0277] Optionally, feature importance evaluation is performed based on the neural network model.

[0278] Correspondingly, on the basis of 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 index is evaluated through the weight change and gradient distribution in the training process.

[0279] The screened parameters are normalized or standardized to make their dimensions consistent and avoid physical scale differences interfering with network training. A ST-GCN-LSTM (spatial-temporal graph convolutional neural network-long short-term memory network) model is constructed, with key tunneling parameters as input and tunneling efficiency, penetration degree, or torque stability as output. The network weights are iteratively optimized using the backpropagation algorithm.

[0280] On this basis, gradient sensitivity analysis, feature perturbation method, or method based on network weight contribution are used to evaluate the influence of each input parameter on the output index, and the normalized feature importance score is obtained. According to the importance score from high to low, select several parameters with the highest contribution as the final optimized feature set to provide compact and efficient input variables for PCA dimension reduction.

[0281] Optionally, PCA is used for high-dimensional feature dimension reduction. A covariance matrix is constructed for the parameter set screened by the neural network, and the principal components that can best represent the tunneling state changes are extracted through principal component analysis, realizing the mapping of high-dimensional features to low-dimensional comprehensive features. A covariance matrix is constructed for the standardized key parameters, and the eigenvalues and eigenvectors are solved to determine the contribution size of each principal component. According to the cumulative contribution rate threshold (e.g., 98%), the first several principal components are selected to compress the data dimension while maintaining the main change information of the tunneling state.

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

[0283] In this application, ring level footage domain transformation technology is adopted. By identifying the driving section through the state machine and establishing a monotonic mapping between time and footage, the non-uniform time series affected by the fluctuation of driving speed is resampled into a standardized footage domain sequence strictly aligned in space, eliminating the geometric distortion of the data waveform, and making the subsequent analysis based on a unified physical space reference.

[0284] Further, a time-lag causal correlation graph and a compensation model are constructed. On the one hand, the idea of causal inference is introduced, and the common interference of mixed variables such as stratum and wear is removed through regression to strip false correlation; on the other hand, the physical response time lag between parameters is calculated, and explicit time lag index translation is performed when the graph neural network features are aggregated. This mechanism enables the model to learn the true physical conduction law rather than false statistical correlation, and realizes accurate quantification and explainable evaluation of the positive and negative effects of key parameters.

[0285] The preferred embodiments of the application are described in detail above, but the application is not limited to the specific details of the above-described embodiments. Within the technical concept of the application, various equivalent transformations of the technical solutions of the application can be made, and these equivalent transformations all belong to the protection scope of the application.

Claims

1. A modeling method based on shield tunneling data feature analysis and parameter correlation, characterized in that, The method comprises the following steps: obtaining multi-source shield tunneling timing parameters in a shield tunneling process, performing state cleaning and coordinate domain transformation on the multi-source shield tunneling timing 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-lag directed association graph model; inputting the standardized footage domain sequence and the time-lag directed association graph model into a graph neural network model, using the time-lag compensation aggregation mechanism in the graph neural network model to learn features, and obtaining a time-lag compensation feature vector corresponding to each parameter; based on the time-lag 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; wherein, using mixed variable elimination and lag correlation analysis to identify the causal time lag and coupling strength between parameters, and constructing a time-lag directed association graph model, comprises: based on the standardized footage domain sequence, calculating the maximum cross-correlation coefficient between each parameter pair in the footage lag window, and determining the initial response time lag of each parameter pair; introducing a mixed variable set, constructing a regression model to eliminate the common influence of the mixed variable set on the standardized footage domain sequence, and obtaining a pure residual sequence corresponding to each parameter; based on the initial response time lag, calculating the lag correlation between the pure residual sequences corresponding to each parameter, and mapping the parameter pairs and the corresponding initial response time lag that meet the causal significance test to the directed edges and edge attributes of the graph model, to obtain a time-lag directed association graph model; wherein, the maximum cross-correlation coefficient and the lag correlation correspond to the coupling strength between parameters; wherein, calculating the lag correlation between the pure residual sequences corresponding to each parameter, and mapping the parameter pairs and the corresponding initial response time lag that meet the causal significance test to the directed edges and edge attributes of the graph model, comprises: 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; calculating the proportion of consistent signs of the local lag correlation coefficients in all footage sliding windows to obtain a parameter association stability score; only keeping the parameter pairs with a 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 edge attributes of the directed edge, to generate a time-lag directed association graph model; The graph neural network model comprises at least one time-lag compensation graph convolution layer; the time-lag compensation aggregation mechanism in the graph neural network model is used for feature learning 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 a normalized footage domain sequence corresponding to each neighbor node, a translation operation based on the initial response time lag is performed on the footage dimension, so that the feature sequence of the neighbor node is physically and causally aligned with the target parameter node; the graph convolution kernel is used to perform weighted aggregation on the neighbor node sequence and the target parameter node sequence after translation alignment, to generate a time-lag compensation feature vector corresponding to the target parameter node.

2. The method of claim 1, wherein, The shield tunneling time sequence parameters are subjected to state cleaning and coordinate domain transformation to construct a normalized footage domain sequence, including: The shield tunneling time sequence parameters are subjected to state recognition by using a multi-parameter joint threshold criterion or a PLC state code to extract a ring tunneling data segment in an effective advancing state; The ring tunneling data segment is divided into a start-up phase, a stable phase and an end phase according to the physical characteristics of the advancing process; For the ring tunneling data segment 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 to generate a normalized footage domain sequence that eliminates the influence of advancing speed fluctuations.

3. The method of claim 1, wherein, According to the gradient information and the edge weight contribution in the prediction process, a key parameter influence degree set is obtained, including: The partial derivatives of the tunneling performance index with respect to the time-lag compensation feature vectors corresponding to the parameters are calculated to obtain gradient vectors containing positive and negative direction information; The gradient vectors are subjected to standardization processing, the edge weights in the time-lag directed association graph model are combined for path aggregation, and the signed comprehensive contribution degrees of the parameters to the tunneling performance index are calculated; The parameters are sorted according to the absolute values of the signed comprehensive contribution degrees, and the parameters are classified into positive promoting factors or negative inhibiting factors according to the sign attributes to form the key parameter influence degree set.

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

5. The method of claim 1, wherein, The graph neural network model is obtained by training through the following steps, specifically including: A training sample set comprising historical normalized footage domain sequences is used to calculate parameter association stability scores of each parameter pair in the training sample set; The 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 values of the edge weights of all directed edges in the time-lag directed association graph model and an 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.

6. The method of claim 1, wherein, After obtaining the set of key parameter influence degrees, further comprising: According to the set of key parameter influence degrees, the top-ranked key parameter sequences are selected from the standardized footage domain sequences to construct a high-dimensional feature matrix; Calculate the covariance matrix of the high-dimensional feature matrix, and perform eigenvalue decomposition on the covariance matrix to determine the variance contribution rate of each principal component; Select the top K principal components with cumulative variance contribution rate reaching a preset threshold as the projection basis, and map the high-dimensional feature matrix to a low-dimensional comprehensive feature vector.

7. The method of claim 1, wherein, The generation of the set of key parameter influence degrees comprises: While keeping the Dropout layer of the graph neural network model in an activated state, the standardized footage domain sequence is subjected to multiple forward propagation sampling to obtain multiple groups of feature importance scores; Calculate the standard deviation of the multiple groups of feature importance scores as an uncertainty indicator, and calculate the reciprocal normalized value of the uncertainty indicator as a fusion weight; The multiple groups of feature importance scores are weighted and averaged using the fusion weight to obtain a final feature importance score with robustness, and the set of key parameter influence degrees is generated accordingly.

Images