Shield data batch processing and cleaning method based on physical mechanics constraints

Abstract

The application discloses a shield data batch processing and cleaning method based on physical mechanics constraints, comprising: obtaining shield data fragments, extracting parameter actual characteristic value domain and matching with a pre-constructed physical typical value domain library, and actively inferring and unifying dimensions; based on extreme value characteristics of data, self-adaptive calculation of working condition threshold, division of physical tunneling working conditions and extraction of pure steady-state data in the advancing stage; performing multi-dimensional physical mechanics constraint verification including array sensor spatial distribution area anomaly detection, multi-parameter physical coupling constraint verification and long-time sequence sensor gradual drift detection on the pure steady-state data, and generating system-level noise labels; based on the labels, calculating data quality scores including physical credibility dimensions, and persistently storing cleaning data carrying labels and scores. The application can effectively identify and eliminate hidden system-level noise that violates engineering mechanics constraints.

Description

Technical Field

[0001] This invention relates to the field of industrial data processing and intelligent construction technology, and in particular to a batch processing and cleaning method for industrial data based on the physical and mechanical constraints of tunnel boring machines. Background Technology

[0002] The construction of large-diameter tunnel boring machines (TBMs) generates massive amounts of multi-dimensional time-series sensor data, including thrust, soil pressure, and torque. High-quality industrial data is the foundation for optimizing TBM tunneling parameters, sensing complex geological formations, and creating digital twin simulations of construction risks. Accurate denoising and standardization of this heterogeneous raw data, characterized by strong temporal sequence and multi-physics coupling, directly determines the convergence stability and physical boundary fidelity of the upper-level artificial intelligence prediction model.

[0003] Current methods for cleaning data from the tunnel boring machine (TBM) industry primarily rely on general time-series preprocessing frameworks. Regarding unit unification, static data dictionaries entered manually are typically used, with fixed coefficient conversions applied to the original values ​​through pre-defined rule tables. For anomaly removal, the Laida criterion (3σ criterion) or median absolute deviation algorithm based on a global data window is widely employed. This involves calculating the statistical mean and variance of the entire dataset to establish upper and lower limit pipelines, thereby directly eliminating outliers exceeding the statistical distribution probability. Simultaneously, operations such as missing value interpolation and timestamp alignment are used to maintain the format integrity of the data matrix.

[0004] However, existing general data cleaning solutions are detached from the physical operation of tunnel boring machines (TBMs), and mainly face problems such as passive reliance on quantities, failure of steady-state statistics, and failure to detect hidden physical false signals. Specifically, during long-term construction, upgrades to the acquisition system are prone to cause abrupt changes in sensor quantities, and simply relying on static dictionaries to perform passive conversions can lead to catastrophic system-level errors. In addition, TBM tunneling naturally involves cyclical alternations of propulsion, assembly, and shutdown. Using global or single-window calculations of statistical thresholds will greatly inflate the data variance due to non-stationary changes during the start-up and shutdown transitions, causing true anomalies within the steady-state during the propulsion period to be masked by excessive tolerance.

[0005] Even more challenging is that conventional univariate probabilistic models struggle to effectively identify the mechanical relationships between cylinder pressures distributed along a circular array or coupled parameters with tight energy conservation relationships. This results in their inability to eliminate hidden hardware noise that conforms to probability distributions but violates physical laws in spatial and mechanical combinations. In fact, they are prone to misjudging and removing multi-zone systematic compliance fluctuations generated when the tunnel boring machine performs serpentine correction. Summary of the Invention

[0006] Purpose of the invention: To provide a batch processing and cleaning method for shield tunneling data based on physical and mechanical constraints, in order to solve the above-mentioned problems existing in the prior art.

[0007] Technical solution: A batch processing and cleaning method for industrial data based on the physical and mechanical constraints of tunnel boring machines, comprising:

[0008] Obtain the shield tunneling industry data fragments, extract the actual characteristic value range of each parameter, and match and standardize them with a pre-built physical typical value range library to obtain dimensionally unified data;

[0009] The working condition threshold is calculated based on the extreme value characteristics of the dimensionless data. The dimensionless data is divided into different physical tunneling working conditions according to the working condition threshold, and pure steady-state data is extracted from the advancement stage within the physical tunneling working condition.

[0010] Multidimensional physical and mechanical constraint verification is performed on pure steady-state data to obtain system-level noise labels. Among them, multidimensional physical and mechanical constraint verification includes anomaly detection for spatial distribution regions of array sensors, multi-parameter physical coupling constraint verification for heterogeneous parameters, and sensor gradual drift detection for long time series.

[0011] Calculate a data quality score that includes the physical credibility dimension based on system-level noise labeling;

[0012] The dimensionally unified data carrying system-level noise labels and data quality scores is persistently stored as cleaned data.

[0013] Beneficial effects: This invention can effectively identify and eliminate hidden system-level noise that violates engineering mechanics constraints. Attached Figure Description

[0014] Figure 1 This is a flowchart of the shield tunneling data batch processing and cleaning method based on physical and mechanical constraints of the present invention.

[0015] Figure 2 This is a flowchart of the present invention for constructing dimensionally unified data.

[0016] Figure 3 This invention is a flowchart of how dimensionally unified data is divided into different physical tunneling conditions based on working condition thresholds.

[0017] Figure 4 This is a simplified flowchart of a processing procedure in one embodiment of the present invention. Detailed Implementation

[0018] Before describing the specific steps, it is necessary to explain the execution environment of this invention. This method is preferably executed in a distributed stream-batch integrated computing engine. For example, a batch processing execution environment can be built based on the data stream application programming interface of a stream processing framework. In this environment, by setting the runtime execution mode to batch processing mode, phased scheduling and resource optimization of bounded datasets are achieved. In this case, the fault tolerance mechanism adopted by the framework is intermediate result materialization and phase-level task re-execution, rather than the periodic snapshot mechanism in streaming mode. Furthermore, incremental batch processing jobs can be triggered by an external scheduler at a preset period to achieve near real-time data access and processing.

[0019] Example 1: A batch processing and cleaning method for shield tunneling data based on multidimensional physical and mechanical constraints, mainly including the following steps:

[0020] Step 101: Obtain shield tunneling industrial data fragments;

[0021] The aforementioned shield tunneling industrial data sharding refers to the formation of self-contained data blocks by cutting massive amounts of raw shield tunneling construction data according to specific rules. In practice, the system can read raw historical data from a time-series database by time range or equipment measurement dimension, and adaptively slice the data based on a set threshold for the amount of data in each shard. To ensure the effectiveness of subsequent spatial and temporal analysis, a hash mapping mechanism is used during the sharding process to ensure that data with the same tunneling ring number belongs to the same shard. In some embodiments, each shard not only contains the target ring data to be cleaned but also pre-loads contextual state information for cross-ring analysis. The data structure can be represented as a multi-dimensional time-series data table containing timestamps and measurement parameters such as thrust, soil chamber pressure, and cutterhead torque.

[0022] Furthermore, after acquiring data shards, optional preprocessing operations on basic fields can be performed. Regular expressions and pre-trained semantic models are used to map various non-standard raw sensor field names to standard unified names, and automatic data type conversion is completed based on confidence level.

[0023] Step 102: Extract the actual feature value range of each parameter in the shield tunneling industrial data segment, match it with the pre-built physical typical value range library, infer and unify the dimensions of each parameter, and obtain dimension-unified data.

[0024] Because sensor labels are often missing or incorrect in practical engineering, it is necessary to infer the physical dimensions from the numerical range of the data itself. The actual characteristic value range characterizes the main numerical fluctuation range of the parameter within the current segment, which can be determined by extracting the high quantile interval of the parameter sequence. The physical typical value range library pre-stores the upper and lower limits of physical empirical extreme values ​​for various types of sensors under different commonly used engineering units. By calculating the overlap between the actual extracted range and the standard range in the library, the system can proactively determine the true acquisition unit of the current data.

[0025] In some optional implementations, if the identified actual dimension is found to be inconsistent with the expected international standard unit, the system will automatically call the corresponding conversion factor to multiply and convert the original value, thereby generating dimension-unified data with clear physical meaning and consistent scale, avoiding the failure of statistical thresholds due to unit confusion in the future.

[0026] Step 103: Adaptively calculate the operating condition threshold based on the extreme value characteristics of the dimensional unified data;

[0027] The tunneling process of a tunnel boring machine (TBM) exhibits strong alternating phases, with completely different data distributions across different phases. To customize a benchmark for each working condition, a dynamic scale, i.e., a working condition threshold, needs to be calculated for the current loop of data. The extreme value characteristic represents a statistical quantity in the data sequence that reflects the true physical upper limit. Considering that directly using the maximum value is easily affected by spike pulses generated by electromagnetic interference from individual sensors, extreme outliers are first removed during the calculation, and the filtered high percentile is selected as the reference benchmark. Then, the threshold is obtained by scaling it according to a certain ratio.

[0028] As an example, this step can adaptively generate a thrust determination threshold to distinguish whether the machine is under thrust by calculating the percentile characteristics of the thrust parameter sequence in the current segment, thereby adapting to the differences in absolute thrust values ​​caused by different geological conditions and machine models.

[0029] Step 104: Divide the dimensionless data into different physical tunneling conditions according to the working condition threshold, and extract pure steady-state data from the propulsion stage within the physical tunneling conditions.

[0030] By performing a joint logical comparison between the parameter values ​​at each moment and the calculated working condition thresholds, the continuous time series can be segmented into physical tunneling working condition segments with clear operational semantics, such as advancement, assembly, and shutdown. Among these, the advancement phase constitutes the majority of the data and contains the most important equipment status information. However, the beginning and end of the advancement phase are usually accompanied by sharp increases or decreases in parameters, and this non-stationary process can pollute the statistical variance. Therefore, by calculating the derivative or rate of change characteristics of the parameters within this phase, the system accurately locates the start and end points of the interval where the parameter changes tend to level off, removes the transitional segments with sharp fluctuations at both ends, and the remaining intermediate stable operating interval is the pure steady-state data.

[0031] In some alternative implementations, if the raw data lacks mechanical signal fields that can directly indicate the assembly process, the system can classify the intervals that do not meet the thrust and speed threshold conditions into non-propulsion stages. This does not affect the accurate extraction of the core propulsion stage and pure steady-state data.

[0032] Step 105: Perform multidimensional physical and mechanical constraint verification on the pure steady-state data to obtain system-level noise labels; wherein, the multidimensional physical and mechanical constraint verification includes anomaly detection for the spatial distribution region of the array sensor, multi-parameter physical coupling constraint verification for heterogeneous parameters, and sensor gradual drift detection for long time series; wherein, the multidimensional physical and mechanical constraint verification is performed based on pre-constructed statistical baseline model parameters;

[0033] This step constitutes a deep data cleansing process. Conventional isolated point detection based on probability distribution cannot detect systemic faults in engineering. Therefore, this step comprehensively utilizes the tunnel boring machine mechanism to construct a three-dimensional physical and mechanical constraint verification system. Specifically, for similar sensors arranged in a circumferential array, anomalies in the spatial distribution area are detected using their inherent spatial topological relationships to distinguish between individual hardware failures and overall mechanical attitude deflection. For heterogeneous parameters with different properties but fixed relationships in fluid mechanics or energy conservation, multi-parameter physical coupling constraint verification is used to verify the fit of their mathematical equations. For macroscopic periods of data sequences exceeding tens of rings, statistical control charts are used to detect gradual sensor drift and capture slowly aging zero-point offsets. The system-level noise marker is a structured set of the output results of the above verification modules, used to identify which physical law each data point violates.

[0034] Optionally, before performing physical verification, the system can first apply the Gaussian distribution rule of statistics to the pure steady-state data to perform a single-parameter anomaly screening. Then, the results of multi-dimensional physical and mechanical constraint verification are used to perform a second verification of these initially screened suspicious points. If the physical model shows that it is reasonable, the anomaly mark is removed, thereby significantly reducing the false alarm rate of the algorithm caused by geological changes.

[0035] Step 106: Calculate a data quality score including the physical credibility dimension based on the system-level noise label;

[0036] To provide quantifiable trust metrics for downstream data analysis models, the underlying validation results need to be transformed into intuitive scores. The physical trust dimension is a metric used to penalize data that violates common engineering principles. During calculation, the system parses system-level noise markers, counts how many spatial, coupling, or temporal physical rules the data point actually violates, and performs a proportional division with the total number of rules it can participate in. This process ensures mathematical fairness in the scoring mechanism between critical parameters with rich constraints and ordinary parameters lacking constraints.

[0037] Furthermore, the score of this physical credibility dimension can be weighted and integrated with common quality indicators such as non-empty rate and format compliance according to the set weights, and finally output a comprehensive data quality score with a standardized range.

[0038] Step 107: Persistently store the dimension-unified data carrying the system-level noise label and the data quality score as cleaned data.

[0039] After end-to-end processing, the data has undergone a transformation from raw, dirty data to high-quality assets. The system updates the original measurement values ​​to values ​​with corrected dimensions and attaches the generated anomaly classification tags and comprehensive scores as metadata to the time series. The persistent storage refers to writing this structured and tagged content into a high-concurrency database, forming a standardized wide table that can be directly accessed by upper-layer business logic.

[0040] In specific system deployment, it is preferable to adopt a storage architecture in which one engineering project corresponds to one time-series database instance and each tunneling ring corresponds to one data table, and adopt a batch asynchronous submission writing mode to support the efficient storage and traceability query of massive measurement point data at the shield tunneling construction site.

[0041] In some embodiments, the clean steady-state data carrying the system-level noise tag and the data quality score, together with its corresponding operating condition tag, are persistently stored as cleaned data along with the remaining data in the dimensionless data.

[0042] Example 2: In shield tunneling engineering data cleaning, the missing units are actively deduced by leveraging the extreme value distribution characteristics of the data itself, and a process of dynamically detecting and correcting sudden changes in units during construction is described in detail below:

[0043] Step 201: Calculate the 90th and 99th quartiles of the target parameter in the shield tunneling industrial data segment, and use the interval formed by the 90th and 99th quartiles as the actual feature value range.

[0044] In shield tunneling data acquisition systems, due to equipment aging or protocol loss, target parameters often only record pure numerical values ​​without unit metadata. To infer their units, it is necessary to extract a numerical range that represents the actual physical working state of the parameter. The target parameter can specifically be a force parameter such as total thrust, soil chamber pressure, or cutterhead torque. During construction, the shield machine spends a significant amount of time assembling segments or undergoing maintenance, during which time the readings of these parameters are numerous and close to zero. If the upper and lower limits of the global data or the conventional fifth to ninety-fifth quantiles are used directly for statistical analysis, the actual high-level working values ​​will be severely diluted by the massive number of zero values ​​from downtime, resulting in a significantly undervalued feature range.

[0045] Specifically, this invention employs a decoupling mechanism, directly arranging the target parameter sequence from the current data slice in ascending order, and extracting data between the 90th and 99th percentiles. This high-percentile interval stably falls within the physical stress range of the tunnel boring machine's full-load propulsion, while the 99th percentile perfectly avoids occasional single-point extreme spike failures of the sensor. In this way, without the need for precise pre-identification of the tunneling conditions, the true magnitude range of the parameter during propulsion work, i.e., the actual characteristic value range, can be independently and purely extracted.

[0046] Step 202: Calculate the interval matching degree between the actual feature value range and the typical value range of each candidate unit in the physical typical value range library;

[0047] The physical typical value range library pre-defines the conventional physical limit ranges of various core parameters of the tunnel boring machine (TBM) under different geological conditions and machine types, in the correct units. For example, the typical value range of total thrust when using kilonewtons as the unit of measurement can be set to 5,000 to 80,000, while the typical value range when using tons as the unit of measurement is set to 500 to 8,000. The candidate unit is the original recording unit that may appear, such as kilonewtons or tons. The system needs to evaluate what proportion of the currently extracted actual feature value range falls within the typical range of a certain candidate unit.

[0048] In the specific calculation, the geometric relationship of interval overlap is used to quantify this possibility, and the calculation formula is as follows:

[0049] m_k = len(R_actual ∩ R_k) / len(R_actual);

[0050] Where m_k is the interval matching degree between the actual feature range and the typical range of the k-th candidate unit, R_actual is the extracted actual feature range, R_k is the typical range of the k-th candidate unit, len represents the length range of the corresponding interval, and ∩ represents taking the intersection of the two intervals on the number axis.

[0051] This formula yields a ratio between 0 and 1. The larger the value, the better the distribution scale of the current data matches the physical common sense corresponding to the candidate unit.

[0052] The physical typical value range library can be pre-constructed based on the model and specifications of the tunnel boring machine (TBM) and engineering geological conditions. It covers commonly used parameter types in TBM tunneling (including but not limited to total thrust, thrust of each zone, cutterhead torque, cutterhead speed, soil chamber pressure, grouting pressure, and grease injection pressure), and pre-stores the reasonable upper and lower limits of each parameter under various commonly used engineering units. The specific values ​​in the physical typical value range library can be determined based on equipment technical specifications, industry engineering standards, and historical construction experience data. Those skilled in the art can expand and adjust it according to actual machine parameters and geological conditions. When a target parameter is not covered by the physical typical value range library, the system can mark the parameter as having an unidentifiable unit, retain the original value, and trigger a manual confirmation prompt in the log.

[0053] Step 203: Select the candidate unit corresponding to the highest interval matching degree as the initial identification unit;

[0054] After completing the traversal calculation, the system selects the candidate unit with the highest score by comparing the interval matching degree of each candidate unit. This unit is then identified as the most likely original collection unit of the current data segment and recorded as the initial identification unit.

[0055] In some optional implementations, a confidence protection mechanism can be set for dimensional identification. When the highest interval matching degree is less than a preset lower limit ratio, or the difference between the highest interval matching degree and the second highest interval matching degree is too small to be clearly distinguishable, the system does not forcibly assign dimensions, but marks the parameter as an unidentifiable unit and triggers an alarm in the system log to prompt manual intervention for confirmation.

[0056] Step 204: Calculate the ring-level average jump ratio of the target parameter between adjacent tunneling rings;

[0057] Long-cycle tunnel boring machine (TBM) projects, after months or years of excavation, are prone to requiring the replacement of all sensors or upgrades to the host computer's data acquisition program. This can lead to sudden changes in the recording unit for the same target parameter during construction. To capture such abrupt temporal jumps, the system calculates the evolution trend of the parameter sequence along the time axis. Specifically, using the tunneling ring as the statistical granularity, the system calculates the average value of the target parameter within each tunneling ring. Then, it divides the average value of the current ring by the average value of the previous ring to obtain the ring-level mean jump ratio, which reflects the order-of-magnitude leap.

[0058] The specific calculation formula is: ρ(k) = x_mean(k) / x_mean(k - 1);

[0059] Where ρ(k) is the ring-level mean jump ratio of the target parameter in the k-th tunneling ring, x_mean(k) is the numerical mean of the target parameter in the k-th tunneling ring, and x_mean(k - 1) is the numerical mean of the target parameter in the immediately preceding tunneling ring.

[0060] Step 205: When the difference between the ring-level mean jump ratio and the preset conversion coefficient ratio between two candidate units is less than the preset tolerance, it is determined that the target parameter has undergone a dimensional change.

[0061] When units change abruptly, the physical forces on adjacent rings usually remain continuous and stable, but the numerical readings will jump by a factor strictly corresponding to the unit conversion factor. The system pre-stores fixed conversion factors for converting various candidate units to SI units. The core logic for determining whether there is a sudden change in dimensions lies in verifying whether the actual ring-level mean jump ratio is consistent with the ratio of the theoretical conversion factor.

[0062] The logical inequality formula for this determination is: |ρ(k) - (c_Ub / c_Ua)| < 0.1 * (c_Ub / c_Ua);

[0063] Where ρ(k) is the ring-level mean jump ratio calculated above, c_Ub is the theoretical conversion factor of the old candidate unit to the international standard unit before the mutation, c_Ua is the theoretical conversion factor of the new candidate unit to the international standard unit after the mutation, |...| represents taking the absolute value of the internal calculation result, and 0.1 represents the preset tolerance ratio.

[0064] For example, suppose the total thrust parameter was originally recorded in tons, with a conversion factor of 9.81 for kilonewtons (kN). Subsequently, the data acquisition system changes, and the target parameter is recorded directly in kN, with a conversion factor of 1.0. At the boundary of the abrupt change, because the recorded value is magnified approximately 9.81 times from representing tons to representing kN, the actual calculated ring-level mean jump ratio is approximately 9.81. Substituting this into the above formula for verification, the old conversion factor 9.81 divided by the new conversion factor 1.0 equals 9.81. The actual jump ratio closely matches the theoretical ratio, with the difference falling within a 10% tolerance range. Therefore, the system definitively determines that the target parameter underwent a dimensional abrupt change from tons to kN at this point.

[0065] Step 206: Based on the initial identification dimension and the judgment result, call the corresponding unit conversion factor to convert the target parameter into international standard units. Specifically, when no dimensional change is determined, the overall conversion is performed directly using the unit conversion factor corresponding to the initial identification dimension; when a dimensional change is determined, the conversion is performed segment by segment using the change ring number as the boundary, calling the corresponding unit conversion factors before and after the change for separate conversions.

[0066] After clarifying the initial identity of the data and potential temporal shifts, the system proceeds to the numerical correction execution phase. If the above determination results indicate that no abrupt change has occurred within the current segmented time period, the system directly uses the unique conversion factor corresponding to the initial identification units to multiply all values ​​of the target parameter sequence by this factor. If the determination results detect a change point, the system segments the data based on the ring number where the change occurred. The old unit conversion factor before the change point is applied to the data segment before the change point, and the new unit conversion factor after the change point is applied to the data segment after the change point.

[0067] Through the above multi-dimensional verification and dynamic compensation, the chaotic engineering readings can be completely transformed into a unified numerical sequence based on international standard units, thus clearing away the obstacle of dimensional confusion for subsequent high-precision noise reduction based on physical constraints.

[0068] Example 3: This example describes how to use the robust extreme value characteristics and derivative change rate of data to sequentially extract the macroscopic physical tunneling conditions from the mixed global shield tunneling data, and further segment out the pure steady-state data at the microscopic level.

[0069] Step 301: Use the median absolute deviation algorithm to filter extreme values ​​in the thrust sequence and speed sequence of the dimensionless data;

[0070] When tunnel boring machines operate in harsh underground environments, the data acquisition system may occasionally record spike pulse values, such as 100,000 kilonewtons, which violate physical limits, due to loose sensor wiring or transient electromagnetic interference. If directly used in statistical analysis, these rare outliers can severely contaminate the calculation of the baseline threshold. This step employs the median absolute deviation algorithm as a coarse-grained pre-filter. This algorithm uses the median of the data sequence as a central reference point, penalizing data that deviates excessively from the median, identifying outlier intervals without assuming a normal distribution. Through this filtering operation, erroneous readings are temporarily excluded from the sample pool for subsequent threshold calculations, ensuring that the subsequent thrust and rotational speed sequences accurately reflect the true nature of the engineering operation.

[0071] Step 302: Extract the 95th quantile of the filtered thrust sequence and speed sequence respectively as the extreme reference benchmarks for thrust and speed.

[0072] After removing destructive spikes, the system needs to find a benchmark value that can represent the full-load working state of the tunneling ring. Traditional techniques usually take the maximum value of the parameters within the ring as the judgment boundary, but this is easily inflated by some local abnormal high values ​​whose amplitude is not large enough to be intercepted by the median absolute deviation algorithm, but still do not represent normal stable working conditions.

[0073] In this embodiment, the system sorts the filtered sequence in ascending order of numerical value and directly extracts the value at the 95th percentile. This extreme value reference represents the upper limit of operation achievable under most normal and stable tunneling conditions. This high percentile value not only faithfully reflects the natural fluctuations in working thrust and cutterhead rotation speed caused by differences in strata hardness, but also avoids the occasional surge interference that may exist in the tail 5%.

[0074] Step 303: Multiply the extreme value reference benchmark by a preset scaling factor to obtain the thrust threshold and speed threshold, which are used as the operating condition threshold.

[0075] To define the physical boundaries of different operational stages, the system generates specific decision thresholds based on the aforementioned extreme value reference benchmark. When operators initiate shield tunneling or cutterhead cutting, parameters are typically gradually increased from zero to the target working value. Multiplying the extreme value reference benchmark by a small, pre-defined scaling factor aims to identify the critical states when the machine is just beginning to be stressed or is about to stop and unload stress.

[0076] For example, the preset proportional coefficient can be set to 0.1. The system calculates the thrust threshold by multiplying the extreme thrust reference base by 0.1, and calculates the speed threshold by multiplying the extreme speed reference base by 0.1. This setting of the working condition threshold means that the system only recognizes that the equipment has truly entered the heavy-load physical operation range when the current force or motion state measured by the sensor exceeds 10% of the maximum normal working capacity of the tunneling ring.

[0077] Step 304: When the data satisfies that the thrust is greater than the thrust threshold and the rotational speed is greater than the rotational speed threshold, it is determined to be in the propulsion stage;

[0078] This step precisely identifies the core construction section of the tunnel boring machine by jointly judging the force and motion dimensions. During normal tunneling, the main hydraulic cylinders must simultaneously exert forward thrust to overcome soil resistance, and the cutterhead must rotate at a sufficient speed to cut the soil. Only when the thrust reading at the same time point crosses the thrust threshold (force threshold) and the rotation speed reading also crosses the rotation speed threshold (motion threshold) does it indicate that the equipment has both forward driving force and cutting action. Based on this, the system classifies the data sequence within that time period as the aforementioned propulsion stage.

[0079] Step 305: When the data satisfies that the thrust is not greater than the thrust threshold, the rotation speed is not greater than the rotation speed threshold, and there is a valid assembly machine start signal in the shield tunneling industrial data segment, it is determined to be in the assembly stage.

[0080] After completing one stage of advancement, the tunnel boring machine (TBM) needs to pause, and the segment assembler installs concrete segments at the tail of the shield. During this period, the main thrust is typically significantly reduced to zero or an extremely low pressure state, and the cutterhead stops rotating. Once the system detects that both the thrust and rotational speed have fallen below their respective extremely low thresholds, and receives a digital switch signal from the equipment's programmable logic controller (PLC) indicating the assembler system's operational status, the system reliably determines this time period as the assembly stage due to the mutual verification of multiple action logics.

[0081] Step 306: Determine the remaining data as part of the shutdown phase;

[0082] Besides continuous tunneling and segment installation, tunnel boring machines also experience long periods of downtime, such as equipment maintenance, standby, or geological exploration. After excluding time segments that clearly conform to propulsion and assembly characteristics, the remaining fragmented time periods that lack heavy thrust and speed support and have not received assembly operation instructions are collected and archived by the system as the shutdown phase.

[0083] In some optional implementations, data acquisition systems for older models or those using specific data protocols may have missing fields, making it impossible to obtain valid assembly machine start signals. As an alternative, when encountering such signal loss, the system can employ a coarse-grained binary search logic. This involves abandoning the distinction between assembly and shutdown actions, and merging all data where the thrust is not greater than a thrust threshold and the rotational speed is not greater than a rotational speed threshold, classifying them as non-propulsion conditions. While this alternative approach sacrifices some precision in condition analysis, it does not affect the lossless extraction and subsequent noise reduction of core propulsion stage data.

[0084] Step 307, the advancing stage, the assembly stage and the shutdown stage constitute the physical tunneling condition.

[0085] At this point, the original time series, which was continuous and lacked semantic labels, was completely deconstructed. The system used the dynamically generated thresholds mentioned above to assign an identity label to each data point, reflecting its physical mechanical action at that time. This completed the data dimensionality reduction at the macro level and formed a set of physical tunneling conditions with a clear structure.

[0086] Step 308: Calculate the difference quotient of the thrust parameters at both ends of the set sliding window during the propulsion phase to obtain the smooth change rate;

[0087] While the above steps isolate the propulsion phase, the propulsion phase itself is not a uniform process from beginning to end. Within the first few tens of seconds of each propulsion, the thrust rapidly increases from zero to tens of thousands of kilonewtons; before the propulsion ends, the thrust is rapidly unloaded. If this transient transition process, involving drastic monotonous rises and falls, is not eliminated, it will greatly inflate the overall data variance. Therefore, the system undergoes a secondary segmentation at the microscopic level, using derivative characteristics to find the true steady state.

[0088] To avoid interference from high-frequency noise in derivative calculations, the system abandons point-to-point direct differencing and instead employs a method of calculating the difference quotient at both ends of a set time window. This method is equivalent to calculating the slope of a secant line over a macroscopic time span.

[0089] The specific calculation formula is: F_dot(t) = (F(t) - F(t - W_s * Δt)) / (W_s * Δt);

[0090] Where F_dot(t) is the smoothed rate of change at time t, F(t) is the thrust amplitude record at the current time t, F(t - W_s * Δt) is the thrust amplitude record at a set sliding window position along the time axis backward, W_s is the set sliding window consisting of the number of consecutive sampling points, and Δt is the fixed time sampling interval between two adjacent data points.

[0091] By calculating the absolute difference between the two endpoints across a window length and dividing it by the span time, the system obtains an extremely smooth and physically meaningful derivative sequence. Positive values ​​in this sequence represent the equipment being loaded and climbing, negative values ​​represent the equipment being unloaded and shutting down, and values ​​approaching zero represent the equipment entering a state of force equilibrium.

[0092] In one optional implementation, the set sliding window can cover approximately 30 to 60 consecutive sampling points (the specific value can be determined based on the sensor sampling frequency and the thrust loading rate of the tunnel boring machine); the preset rate of change threshold can be estimated based on the rated thrust of the tunnel boring machine and the typical propulsion loading time (for example, for a model with a rated thrust of approximately 30,000 kN and a loading time of approximately 60 seconds, the rate of change in the steady-state phase is typically less than 50 kN / s, and this value can be set as a reference benchmark); the preset ratio of the thrust extreme value reference benchmark can be set to a value between 0.6 and 0.8 to ensure that steady-state is only confirmed when the thrust reaches the normal tunneling working level. Those skilled in the art can adaptively adjust the above parameters based on specific machine parameters and geological conditions using engineering experience or data-driven methods, without requiring creative effort.

[0093] Step 309: Search for the first moment within the propulsion phase where the absolute value of the smooth change rate is not greater than a preset change rate threshold and the thrust amplitude is greater than a preset ratio of the thrust extreme value reference benchmark, and record it as the steady-state start moment.

[0094] A thrust amplitude greater than a preset ratio of the thrust extreme value reference benchmark can, in some scenarios, be understood as: the thrust amplitude greater than the product of the thrust extreme value reference benchmark and a preset high-level ratio coefficient. Alternatively, it can be understood as the thrust amplitude greater than the extreme value reference benchmark multiplied by a preset ratio.

[0095] The system searches forward point by point along the time axis from the starting point of the propulsion phase. To confirm that the equipment has truly entered a stable physical working state, two extremely stringent conditions must be met simultaneously. First, the absolute value of the smoothed rate of change must shrink to an extremely low set range, meaning that the transient process of a sharp thrust increase has ended, and the parameter curve exhibits horizontal oscillation characteristics. Second, the absolute value of the thrust at this point must have climbed to a high level sufficient for normal cutting operations. The system compares the thrust amplitude with the extreme thrust reference value extracted in the previous steps. Only when the amplitude exceeds a specific high proportion of this extreme reference value, coupled with the condition that the slope approaches zero, will the system mark that time point, establishing it as the start of steady state. Without this amplitude constraint, the system is highly susceptible to misjudging the brief pause in the initial propulsion phase, where the thrust is maintained at a low level to adjust attitude, as a normal steady state.

[0096] Step 310: Search for the last time in the propulsion phase that satisfies the same judgment condition, and record it as the steady-state end time;

[0097] Similar to the forward search, the system performs a backward search along the timeline from the end of the advancement phase. Since the tunnel boring machine typically experiences a relatively long period of stable tunneling before the shutdown command is issued, the system searches backward for the critical turning point where the thrust begins to experience a substantial monotonic decline. This point, found during the backward search, is the last point where the thrust amplitude remains within a high range while the rate of change slope remains within a slightly fluctuating range. This point is identified as the end point of the current stable operation and recorded as the steady-state termination moment.

[0098] Step 311: Extract the data segment between the steady-state start time and the steady-state end time as the pure steady-state data to remove the start-up transition segment and shutdown transition segment within the propulsion phase.

[0099] The data located between these two anchor points, having eliminated the disturbances of sudden increases during loading and sudden drops during unloading, exhibits an excellent unimodal Gaussian clustering characteristic and is extracted and named pure steady-state data. The data segments remaining before the steady-state start time and after the steady-state end time are isolated and defined as the start-up transition segment and the shutdown transition segment, respectively. This dual segmentation method, based on differential derivatives and physical boundaries, reshapes the data environment for the statistical denoising algorithm in subsequent embodiments, preventing method failure due to mixed distributions.

[0100] Example 4: After extracting clean steady-state data, cross-validation is performed by combining single-parameter statistical features with multi-parameter physical models to achieve the process of identifying abnormal data and preventing false alarms and downgrading. This mainly includes the following steps:

[0101] Step 401: Calculate the steady-state mean and steady-state standard deviation of the pure steady-state data under the current physical tunneling conditions;

[0102] After excluding transient processes during startup and shutdown, the fluctuations of a single parameter in the pure steady-state data better conform to the characteristics of a normal distribution. The system calculates the arithmetic mean of all discrete sampling points for this data segment, denoted as the steady-state mean; simultaneously, it calculates the standard deviation of the dispersion of each sampling point relative to the steady-state mean, denoted as the steady-state standard deviation. This step provides the basic calculation parameters for subsequently constructing a filtering boundary based on statistical probability.

[0103] Step 402: Construct a piecewise steady-state filtering threshold by adding or subtracting three times the steady-state standard deviation from the steady-state mean;

[0104] Based on the normal distribution theory, the system uses the three-standard-deviation rule to determine the reasonable distribution range of normal data. Specifically, the system uses the aforementioned steady-state mean as a benchmark, shifting it upwards and downwards by three steady-state standard deviations respectively to generate the upper and lower limits of the interval. These upper and lower limits together constitute the segmented steady-state filtering threshold. Since this threshold is dynamically generated based on the current independent steady-state segment data, it can adapt to the overall shift of the data baseline values ​​under different geological conditions.

[0105] Step 403: Initially mark the data points in the pure steady-state data that exceed the segmented steady-state filtering threshold as data with suspected statistical noise.

[0106] The system compares each value in the clean steady-state data sequence with the aforementioned segmented steady-state filtering thresholds point by point. When a value at a certain timestamp is greater than the upper limit or less than the lower limit, it indicates that the value is statistically an anomaly, usually caused by random sensor noise or signal interference. The system extracts these data points that exceed the threshold range and adds an anomaly status label, recording them as data with suspected statistical noise. It should be noted that, due to the limitations of the univariate statistical model, actual stress mutations caused by the tunnel boring machine encountering localized hard geology can also lead to values ​​exceeding the limits. Therefore, this labeling is a preliminary screening result and needs to proceed to the subsequent physical verification process.

[0107] Step 404: For the data of the start-up transition segment and the shutdown transition segment that are stripped out during the propulsion phase, the segmented steady-state filtering threshold is not checked.

[0108] During the overall implementation phase, the data from the start-up transition and shutdown transition phases exhibit a monotonically rapid upward or downward trend, failing to meet the statistical distribution assumption of steady fluctuations. Applying the aforementioned steady-state filtering threshold to these two data segments would result in the incorrect labeling of normal loading and unloading data points. Therefore, the system employs isolation rules in its logical control to retain these two data segments while skipping the statistical filtering step, ensuring the effectiveness and accuracy of statistical anomaly detection.

[0109] Step 405: Substitute the data carrying the suspected statistical noise label into the pre-fitted linear coupling baseline model and calculate its standardized coupling residuals on each physical parameter coupling pair.

[0110] For suspected abnormal data identified in the initial screening, the system calls all physical correlation models containing that parameter for secondary verification. The thrust and soil pressure, torque and rotational speed of the tunnel boring machine (TBM) system are constrained by rigid mechanical laws. The system extracts the tagged data, inputs the actual recorded values ​​of its correlated parameters into the corresponding linearly coupled baseline model, and calculates the theoretical predicted value of the target parameter. Subsequently, the difference between the actual recorded value and the theoretical predicted value of the target parameter is calculated and divided by the standard deviation of the residuals of the historical fitted model to obtain a quantitative index, namely the standardized coupling residual. This residual value is used to measure whether the current combination state of this set of parameters deviates from the normal physical correlation law.

[0111] Step 406: If the absolute value of the standardized coupling residuals of the data in all corresponding linear coupled baseline models is not greater than the preset residual threshold, it indicates that the data fluctuations of the data conform to the physical and mechanical coupling law. Then, its suspected statistical noise label is removed and it is downgraded to normal observation deviation data.

[0112] In the physical verification logic, if the absolute values ​​of the standardized coupling residuals calculated for a target data point in multiple independent but interconnected coupled models are all less than or equal to a preset residual threshold, it indicates that the data exhibits synergistic variation within the multidimensional physical space. The physical essence of this phenomenon is that the numerical abrupt change does not originate from independent random noise in the sensor hardware, but rather from the overall mechanical response generated when the tunnel boring machine's cutting system encounters a sudden change in geological resistance. Based on this objective judgment, the system performs logical correction, removing any suspected statistical noise markers attached to the data point and updating its status label to normal observation deviation data. This cross-validation degradation mechanism effectively prevents extreme data containing important geological stress characteristics from being incorrectly rejected by statistical algorithms.

[0113] Step 407: If the absolute value of the standardized coupling residual in any corresponding linear coupled baseline model is greater than the preset residual threshold, then the original determination is maintained, it is confirmed as real noise and incorporated into the system-level noise label, or in other words, it is incorporated into the system-level noise label.

[0114] Conversely, if a data point marked with suspected statistical noise exhibits a sudden numerical change in a single time series, but its associated parameters do not show a coordinated physical response, resulting in the calculated absolute value of the standardized coupling residual exceeding the preset residual threshold after being substituted into any linearly coupled baseline model, this isolated abrupt change, violating the law of physical conservation, is typically caused by local hardware failures or signal crosstalk in the sensor. Therefore, the system determines that the statistical anomaly of this data point is valid, maintains the original anomaly determination, classifies it as real noise, and merges this anomaly record into the global system-level noise label for subsequent data isolation and removal during the data entry process.

[0115] Example 5: A sensor array anomaly detection method based on the spatial isolation criterion is provided. The main process is as follows:

[0116] Step 501: For array sensor parameters with spatial partitioning structure, calculate the spatial deviation of each partition parameter in the pure steady-state data relative to the mean of all partitions at the same time.

[0117] The front-end cutting and propulsion systems of large-diameter tunnel boring machines (TBMs) are not composed of a single actuator, but rather employ a highly modular array design. For example, to propel the heavy shield and control the tunneling direction, the main propulsion cylinders are typically divided into multiple independent sections along the circumference of the shield shell, arranged evenly in the upper, lower, left, and right directions. Similarly, to sense the slurry resistance at the excavation face, soil chamber pressure sensors are deployed in a multi-point array at the top, waist, and bottom of the chamber wall behind the cutterhead. The values ​​collected by these ring-arranged measuring elements constitute the parameters of the array sensor with the spatial partitioning structure.

[0118] During normal, stable linear tunneling, since the same rigid shield is being pushed or the same face is being cut into the soil, theoretically, the readings of these ring-shaped zone sensors should remain at a similar level at any given moment. If the value of one zone deviates significantly from the others, this abrupt distribution indicates a strong anomaly. Therefore, the system extracts the clean, steady-state data stream and, at each specific time snapshot slice, calculates the average of all zone readings, representing the overall macroscopic center of gravity of the tunnel boring machine at that moment. Then, the system subtracts this global center of gravity value from the reading of each individual zone, thus identifying the degree to which each zone deviates from the overall value; this is the spatial deviation.

[0119] The specific calculation formula is: d_i(t) = F_i(t) - F_mean(t);

[0120] Where d_i(t) is the spatial deviation of the i-th partition calculated at time t, F_i(t) is the actual parameter reading of the i-th partition at time t, and F_mean(t) is the arithmetic mean of the readings of all partitions of the array sensor system at time t.

[0121] Through this subtraction operation, the system reduces the absolute value of the thrust, which was originally tens of thousands of kilonewtons, into a relatively clean error fluctuation band, greatly highlighting local anomalies.

[0122] Step 502: Using a pre-constructed historical spatial deviation baseline, the spatial deviation is converted into a standardized spatial residual;

[0123] Obtaining only an absolute deviation value is insufficient to support anomaly detection. Even under ideal tunneling conditions, due to minute manufacturing tolerances in the internal frictional resistance of each cylinder, there is a natural fixed static error of several hundred kilonewtons between different sections. These biases cannot be treated as faults and intercepted. Therefore, the system introduces the pre-built historical spatial deviation baseline as a reference.

[0124] This baseline is calculated by the system in the background using clean data confirmed to be healthy from dozens or even hundreds of past tunneling loops, and is statistically analyzed for each specific partition. It includes the inherent mean and allowable variance of each partition. The system first subtracts the historical mean of the partition from the real-time spatial deviation calculated at the current moment, and then divides it by the historical allowable standard deviation of the partition, thereby completing the dimensionless normalization process and outputting the standardized spatial residual.

[0125] The specific standardized formula is: z_sp_i(t) = (d_i(t) - μ_di) / σ_di;

[0126] Where z_sp_i(t) is the standardized spatial residual of the i-th partition at time t, d_i(t) is the current spatial deviation calculated above, μ_di is the inherent mean habit of the partition statistically derived from the historical spatial deviation baseline, and σ_di is the allowable standard deviation of the partition statistically derived from the historical spatial deviation baseline.

[0127] This transformation process not only eliminates the comparison barriers caused by different dimensions and absolute magnitudes, but also eliminates the systematic fixed bias caused by the inherent mechanical tolerances of the equipment.

[0128] Step 503: When the absolute value of the standardized spatial residual of a single partition parameter is greater than the preset residual threshold, and the total number of partitions with the absolute value of the standardized spatial residual greater than the preset residual threshold at the same time does not exceed one-third of the total number of partitions, it is determined that the spatial isolation criterion is met.

[0129] The spatial isolation criterion aims to address the false alarm problem that easily occurs in the detection of spatial dimensional anomalies. During actual tunneling operations, when performing axis correction or attitude control operations, the system typically increases the thrust of the cylinder on one side of the circumference and correspondingly decreases the thrust of the cylinder on the opposite side. This compliant engineering attitude intervention can cause the spatial residuals of multiple zones in the array sensor to simultaneously exceed the normal fluctuation range. If only the residual exceeding the standard of a single zone is used as the absolute criterion, it is highly likely that valid data reflecting the systematic attitude changes of the equipment will be misjudged as sensor noise and discarded.

[0130] To effectively distinguish between normal attitude control and genuine sensor malfunctions, this method constructs the aforementioned verification logic by combining local deviations and global spatial topology relationships. Specifically, the system not only detects whether the absolute value of the standardized spatial residual of a single partition is greater than the preset residual threshold, but also counts the total number of partitions exceeding the preset residual threshold at the same timestamp. The spatial isolation criterion is satisfied only if the total number of such partitions is less than or equal to one-third of the total number of partitions in the array sensor. The physical meaning of this criterion is that the out-of-bounds behavior of a few partitions indicates a local sensor hardware anomaly, while the coordinated out-of-bounds behavior of more than one-third of the partitions represents a macroscopic adjustment of the overall rigid body attitude of the tunnel boring machine. This achieves precise false alarm protection against global attitude compliance fluctuations at the algorithm level.

[0131] Assume the current monitored machine is a main tunnel boring machine equipped with six thrust zones, meaning the total number of zones M equals six. In this case, the red line that triggers the judgment logic is a rounded-down integer of two.

[0132] Scenario 1: In a single second, only one of the six zones, located on the upper right side of the shield, experiences an abnormal drop in sensor readings, with a residual error reaching 5.8, severely exceeding the acceptable range. Meanwhile, the residual errors of the other five zones remain well near zero. Since the total number of zones exceeding the acceptable range is only one, less than or equal to the two specified by the red line, this perfectly conforms to the spatial isolation criterion. In the physical world, this translates to only one sensor hardware failure.

[0133] Scenario Two: In another second, the operator sharply turns the steering wheel to climb an incline. At this moment, the two thrust zones at the bottom begin to rapidly increase pressure to provide lift, while the two thrust zones at the top begin to desperately depressurize. Instantly, the residual values ​​of all four zones exceed the threshold alarm. Since the total number of simultaneous exceedances is as high as four, far exceeding the red line limit of two, this scenario is problematic.

[0134] The above phenomena indicate that multiple stress zones of the tunnel boring machine (TBM) are undergoing coordinated, systemic stress changes. The physical essence of this is a macroscopic adjustment of the overall posture of the equipment, rather than an isolated malfunction of a single sensor. Therefore, the system determines that this situation does not meet the spatial isolation criterion and does not mark it as an anomaly. This logic effectively avoids misjudging the systemic compliance fluctuations generated by the TBM during posture correction as sensor noise, ensuring the complete preservation of valid engineering data.

[0135] Step 504: Data points that satisfy the spatial isolation criterion are marked as spatial distribution anomalies caused by local sensor faults and incorporated into the system-level noise label. Alternatively, this label is merged into the system-level noise label.

[0136] For data points that meet the spatial isolation criterion, the deviation exceeds the limit and lacks the spatial coordination that the array sensor system should have, indicating that the numerical mutation does not originate from the actual physical force change of the tunnel boring machine, but is usually caused by sensor circuit noise, communication interference or local mechanical detection end failure.

[0137] The system identifies and marks these isolated outlier data points as spatial distribution anomalies caused by local sensor malfunctions, and incorporates these anomaly labels into the system-level noise labeling. This processing effectively transforms spatial topological relationships into computable discrimination rules, thereby providing an objective judgment basis based on both physical and spatial constraints for downstream data quality physical reliability scoring calculations, as well as for the final isolation and removal of outlier data.

[0138] Example 6: In a multi-dimensional physical and mechanical constraint verification system, the specific process for detecting hidden anomalies such as long-term slow aging or minute zero-point shift of sensors is as follows:

[0139] Step 601: Using the shield tunneling ring as the aggregation granularity, calculate the ring-level mean of each parameter in the pure steady-state data, and construct a ring-level mean sequence;

[0140] The high temperature, high humidity, and continuous vibration in the tunnel boring machine (TBM) construction environment can easily cause gradual zero-point drift in sensors. Because the magnitude of this type of drift is extremely small, typically far below the statistical threshold built from single-point data, conventional high-frequency transient anomaly detection methods struggle to identify it. Therefore, the system employs a frequency reduction and aggregation strategy, using a complete TBM tunneling ring as a time statistical unit. For the extracted clean steady-state data, the system calculates the arithmetic mean of the values ​​at each discrete sampling point within the same tunneling ring, obtaining a representative characteristic value reflecting the overall parameter level of that ring. By concatenating the means of multiple consecutive tunneling rings in chronological order, high-frequency fluctuation interference can be eliminated, forming the ring-level mean sequence used for macro-trend analysis.

[0141] Step 602: Extract the sliding median of the ring-level mean sequence within the preset historical ring number window as the drift reference benchmark;

[0142] As the tunnel boring machine (TBM) traverses different geological layers, the overall average values ​​of parameters such as thrust or torque will change gradually and systematically with variations in soil hardness. To accurately distinguish between gradual changes in soil physics and sensor element drift, the system extracts a continuous historical loop window of a set length on the time axis, specifically the fifty tunneling loops prior to the current moment. The system extracts the average values ​​of all loop levels within this window and sorts them by value, selecting the moving median at the middle position as the drift reference benchmark. Compared to the arithmetic mean, the median is more robust to the data from a few extreme fault loops that may be mixed in within the window, and it can adaptively track and absorb normal baseline fluctuations caused by changes in soil physics as the time window progresses.

[0143] Step 603: Calculate the upper and lower cumulative sum control chart statistics of the ring-level mean sequence using the drift reference benchmark;

[0144] To address the cumulative amplification effect of small deviations, the system employs a cumulative sum control chart algorithm for verification. The system establishes an upper-side cumulative sum to monitor a continuous upward trend in values ​​and a lower-side cumulative sum to monitor a continuous downward trend. For the latest input ring-level mean, the system calculates its difference from the drift reference baseline and subtracts a preset tolerance bias representing the allowable fluctuation of the device. If the difference after subtraction is positive, it is added to the upper-side statistic; if it is negative, it is converted to a positive value and added to the lower-side statistic. If the parameter only exhibits random fluctuations around the baseline, positive and negative differences will cancel each other out, and the cumulative sum will remain near zero; if there is a unidirectional gradual drift, small differences will be continuously accumulated and amplified.

[0145] Its calculation formula can be expressed as:

[0146] C_plus(k) = max(0, C_plus(k - 1) + x_mean(k) - x_ref - δ);

[0147] C_minus(k) = max(0, C_minus(k - 1) - x_mean(k) + x_ref - δ);

[0148] Where C_plus(k) is the upper cumulative sum control chart statistic of the k-th ring, C_minus(k) is the lower cumulative sum control chart statistic of the k-th ring, x_mean(k) is the ring-level mean of the k-th ring, x_ref is the drift reference benchmark, δ is the set tolerance bias, and max indicates taking the maximum value among the values ​​inside the brackets.

[0149] Step 604: When the cumulative sum of the upper and lower sides and the control chart statistics exceed the preset control limit, it is determined that the corresponding parameter has undergone gradual drift.

[0150] The system sets a circuit breaker threshold for the aforementioned accumulation process, namely the preset control limit. In specific implementation, this control limit can be set to five times the historical standard deviation within the window containing the aforementioned drift reference benchmark. The system continuously monitors and updates the statistical values. When the statistics of the upper or lower cumulative control chart exceed the preset control limit, it indicates that the cumulative effect of the unidirectional shift of the parameter has exceeded the reasonable boundary of statistical random fluctuation. Based on this, the system determines that the parameter corresponding to the sensor has undergone a systematic gradual drift.

[0151] Step 605: Determine the drift start loop number by backtracking the cumulative starting point of the control chart, mark the data points in the corresponding drift interval as gradual drift anomalies, and incorporate them into the system-level noise label.

[0152] After confirming the occurrence of drift, to achieve source tracing and cleaning of contaminated historical data, the system traces backward along the time axis in the historical records of the cumulative data and control charts. The system identifies the first tunneling ring whose statistical value continuously accumulates from zero without returning to zero as the drift initiation ring number. Subsequently, the system extracts the complete time range from this drift initiation ring number to the ring number where the current alarm is triggered, classifying all detailed data points falling within this drift interval as anomalies and assigning them a gradual drift anomaly label. This label is ultimately incorporated into the system-level noise tag, thereby seamlessly integrating the system deviation detection results from the long-term time series dimension into the global multi-dimensional physical constraint evaluation system.

[0153] Example 7: This describes how to construct a scoring calculation mechanism with adaptive attributes after obtaining multi-dimensional anomaly labels, and how to transform multi-dimensional physical verification results into a comprehensive data quality score that includes a physical credibility dimension.

[0154] Step 701: Analyze the system-level noise markers, count the number of array sensor spatial partitions, physical parameter coupling pairs, and gradual drift detection items actually triggered by the current data point, and sum them up to obtain the total number of dynamic verification items;

[0155] When dealing with shield tunneling parameters exhibiting highly heterogeneous characteristics, the dimensions and number of physical constraints imposed on different sensor parameters vary. For example, the total thrust parameter is configured with multiple spatial partitions and multiple physical-mechanical coupling pairs, while the temperature parameter may lack a spatial array structure and have fewer associated coupling pairs. To eliminate the impact of these parameter attribute differences on the quality evaluation benchmark, the system employs a calculation mechanism that dynamically adjusts the evaluation benchmark.

[0156] In practice, the system parses the metadata of the system-level noise marker corresponding to the data point and extracts the physical verification dimension that the specific parameter actually participates in the calculation at the current moment. The system arithmetically sums the number of array sensor spatial partitions, the number of physical parameter coupling pairs participating in the comparison, and the number of gradual drift detection items executed for the parameter to generate the total number of dynamic verification items. This total number represents the total amount of effective physical rule constraints that the specific parameter is subject to in the system.

[0157] Step 702: Count the number of verification items confirmed as abnormal in the system-level noise markers, calculate the ratio of the number of abnormal verification items to the total number of dynamic verification items, and obtain the dimensionless penalty ratio.

[0158] After determining the effective verification dimension cardinality of the current parameter, the system quantifies the degree of data deviation based on the judgment result. The system traverses the system-level noise markers, extracts the entries marked as spatial distribution anomalies, hidden combination anomalies, and gradual drift anomalies, and counts the absolute number of independent verification items confirmed as anomalies.

[0159] Subsequently, the system performs a division operation, using the number of abnormal verification items as the numerator and the total number of dynamic verification items as the denominator. The calculated quotient, between zero and one, is the dimensionless penalty ratio. This division calculation uniformly maps deviation events with different physical dimensions and different fault types to a mathematical ratio, ensuring the consistency of the evaluation scale between core parameters with more physical constraints and ordinary parameters with fewer constraints.

[0160] Step 703: Calculate the physical credibility score based on the dimensionless penalty ratio;

[0161] After obtaining the dimensionless penalty ratio, the system converts it into a percentage score, which serves as a quantitative indicator representing the degree to which the data conforms to the laws of physics and mechanics.

[0162] The specific formula for physical credibility scoring is: S_phys = 100 * (1 - N_err / N_total);

[0163] Where S_phys is the calculated physical reliability score, N_err is the number of verification items confirmed as abnormal, and N_total is the total number of dynamic verification items.

[0164] Taking the total thrust parameter, which includes six spatial partitions, two coupling pairs, and one time-series drift detection, as an example, the total number of dynamic verification items is nine. If this parameter experiences only one residual exceeding the standard anomaly during the verification process, its physical reliability score is approximately 88.9 points when substituted into the formula. This calculation logic reflects that if a parameter exhibits only a few local deviations in multiple physical cross-validations, it still possesses high physical reliability on a macroscopic level. By dynamically adjusting the denominator, the system effectively avoids excessive penalties for core parameters with many configured physical rules, ensuring the fairness of the scoring.

[0165] Step 704: Calculate the integrity score, consistency score, and validity score based on the unified dimensional data; combine the physical credibility score with the pre-calculated integrity score, consistency score, and validity score according to preset weights and sum them to obtain the data quality score.

[0166] When outputting the final data quality assessment, physical reliability, as one of the feature dimensions, needs to be integrated with conventional data cleaning evaluation indicators. Specifically, the system obtains in parallel the integrity score representing the proportion of non-empty fields, the consistency score representing the compliance of field format and units, and the validity score representing the fact that the values ​​are within the reasonable fluctuation range of univariate statistics.

[0167] The system performs a weighted sum of the scores from the four dimensions based on pre-defined weights. For example, the pre-defined weight for the integrity score can be set to 0.2, the pre-defined weight for the consistency score to 0.25, the pre-defined weight for the validity score to 0.25, and the pre-defined weight for the physical reliability score to 0.3. The weighted sum outputs a data quality score between zero and one hundred. This score can serve as a quantitative benchmark for downstream data-driven models, such as those for shield tunneling construction risk warning, to determine whether to accept the time series data.

[0168] The integrity score, consistency score, and validity score are all normalized to a percentage range of 0 to 100. Their calculation methods can adopt the data quality assessment methods commonly used in the field (such as the integrity score being the proportion of non-empty fields multiplied by 100). Those skilled in the art can determine the specific calculation method for each sub-dimension according to the actual data quality management specifications.

[0169] In some embodiments, the data quality score can be obtained by weighting and summing the physical credibility score with the integrity score pre-calculated based on the non-empty ratio of the field in the unified dimensional data, the consistency score pre-calculated based on the field format and dimensional compliance of the unified dimensional data, and the validity score pre-calculated based on the reasonable fluctuation range of single-parameter statistics according to a preset weight combination.

[0170] Example 8 describes a distributed state broadcast optimization and baseline cold start hierarchical strategy. This strategy is used in a distributed batch processing computing environment to optimize the memory overhead of cross-loop historical data through a state broadcast mechanism, and to address the baseline model cold start degradation when the sample size of historical data is insufficient in the early stages of the project.

[0171] Step 801: Extract the statistical baseline model parameters required for cross-loop analysis; wherein, the statistical baseline model parameters include the historical spatial deviation baseline, linear coupling baseline model, and drift reference benchmark required for the multidimensional physical and mechanical constraint verification;

[0172] In the data processing of large-scale tunnel boring machine (TBM) projects, high-frequency sampling by sensors generates massive amounts of time-series records. Verification of multidimensional physical and mechanical constraints (such as spatial distribution consistency, parameter coupling, and gradual drift verification) heavily relies on cross-ring analysis of data from dozens of consecutive historical tunneling rings. To avoid directly reading and processing this vast amount of detailed data on distributed computing nodes, the system pre-executes lightweight data feature extraction before the main cleaning process. Specifically, the system extracts pure mathematical statistics from the cleaned historical data, including the mean and standard deviation of each partition, the slope and intercept of the multivariate least squares fit, and the moving median of the sequence. This set of numerical values, which highly summarizes the historical physical laws, constitutes the parameters of the statistical baseline model.

[0173] Step 802: Store the statistical baseline model parameters into the global state cache of the distributed batch processing execution environment;

[0174] After obtaining the above parameters, the system uses them as shared basic rule configurations and stores them in a high-speed access area pre-allocated by the underlying computing framework. When selecting a distributed computing engine with integrated stream and batch processing capabilities, this area can be specifically configured as a global state cache for key-value pair data structures to facilitate efficient subsequent distribution and retrieval.

[0175] Step 803: The updated statistical baseline model parameters are pushed to each parallel processing node using a state broadcast mechanism so that the multidimensional physical and mechanical constraint verification can be directly called for calculation, avoiding the loading and splicing of original historical data across nodes in the shield tunneling industrial data sharding.

[0176] During concurrent data cleaning, for each shard of shield tunneling industrial data to be processed, the system employs a state broadcast mechanism to completely copy and push the statistical baseline model parameters from the global state cache to the local memory of each parallel processing node. When a computing node executes verification logic such as multi-parameter physical coupling, it can directly call the slope and intercept in its local memory to calculate the residual formula. This architecture design uses lightweight model parameters to replace cross-node data transfer, completely avoiding the massive network input / output overhead caused by physically splicing dozens of rings of original historical data with the current data shard, and effectively preventing memory overflow failures in parallel nodes.

[0177] However, the construction of the aforementioned baseline model relies on sufficient historical sample accumulation. For newly constructed tunnel boring machine (TBM) projects in the initial tunneling stage, the system is configured with the following cold-start layered degradation strategy to ensure the smooth operation of the algorithm.

[0178] Step 804: Obtain the cumulative number of completed cycles for the current project;

[0179] Before processing each batch, the system queries the project master database to obtain the cumulative number of completed cycles, reflecting the absolute progress of the project. This value serves as the sole criterion for evaluating the degree of historical data accumulation and is used to trigger the corresponding model running state.

[0180] Step 805: When the cumulative number of completed rings is lower than the first preset ring number threshold, the multidimensional physical and mechanical constraint verification is not enabled.

[0181] As an example, the first preset ring count threshold can be set to ten rings. When the cumulative number of completed rings is less than ten, the available effective historical sample size is extremely small, which cannot meet the basic degrees of freedom required for linear regression or statistical distribution. If a baseline model is forcibly fitted at this time, it will produce misleading physical parameter equations. Therefore, the system triggers a blocking mechanism at this stage, suspending the use of advanced verification modules such as spatial distribution, parameter coupling, and gradual drift, and only retaining the basic single-parameter filtering function based on extreme values ​​to avoid systemic false alarms caused by the immaturity of the model itself.

[0182] Step 806: When the cumulative number of completed rings is between the first preset ring number threshold and the second preset ring number threshold, the statistical baseline model parameters are fitted using existing historical data, and the judgment threshold of the multidimensional physical and mechanical constraint verification is increased by a preset relaxation ratio to perform confidence compensation.

[0183] For example, the second preset ring count threshold can be set to fifty rings. When the cumulative number of completed rings is between ten and forty-nine rings, the system has a preliminary statistical sample and begins to calculate and fit the statistical baseline model parameters based on the existing data. Considering that the model variance is large under small sample conditions, the system simultaneously executes a confidence compensation mechanism.

[0184] Specific operations include one or more of the following methods: First, increase the preset relaxation ratio (e.g., from the default two times to three times) of the standardized residual judgment threshold used in multidimensional physical and mechanical constraint verification (i.e., the preset residual threshold for judging whether the standardized spatial residual and standardized coupling residual exceed the standard) based on the standard value, making the standardized residual judgment range more lenient; second, increase the sensor gradual drift detection control limit based on the cumulative sum control chart by 50% based on the default value. The above two confidence compensation measures target the judgment parameters of multidimensional physical and mechanical constraint verification and do not affect the three-standard-deviation filtering threshold dynamically calculated based on the current steady-state segment in single-parameter piecewise steady-state filtering.

[0185] Step 807: When the cumulative number of completed rings reaches the second preset ring threshold, the statistical baseline model parameters are updated using historical data within the sliding window, and the standard judgment threshold is restored.

[0186] When the cumulative number of completed rings reaches fifty or more, the sample size meets the statistical stability requirements. The system exits the cold start protection state and enters the standard operation mode. The system utilizes a sliding window with a set length of fifty rings that moves forward synchronously with the project progress to continuously filter the latest historical data to iteratively update the statistical baseline model parameters, adapting them to the current geological conditions. Simultaneously, the system cancels all previously set tolerance compensations, restoring all judgment thresholds to the most stringent standard judgment thresholds, thereby achieving high-precision, comprehensive physical and mechanical constraint cleaning for massive amounts of shield tunneling data.

[0187] Example 9 describes a multi-parameter physical coupling constraint cross-validation mechanism to identify hidden systematic biases.

[0188] Step 901: Obtain the linear coupled baseline model corresponding to multiple physical parameter coupling pairs that have been prefitted using historical steady-state data;

[0189] During the tunneling process, the various subsystems of a tunnel boring machine (TBM), including hydraulic drive, electrical power, and mechanical cutting, exhibit inherent energy and mechanical transmission relationships. For example, the total thrust required to overcome ground resistance and the soil pressure in front of the cutterhead exhibit a positive synergistic characteristic; the rotational speed of the cutterhead drive motor and the torque of the cutterhead used to overcome friction also have a specific power response mechanism. These parameter combinations, which reflect the real physical relationships, constitute the aforementioned physical parameter coupling pairs.

[0190] To transform this qualitative physical relationship into a computable quantitative scale, the system pre-retrieves a set of historical data from a stable tunneling state. For each selected pair of coupling parameters, the system uses linear regression algorithms such as least squares to fit the data, extracting the slope and intercept of the physical correlation, thereby establishing an equation reflecting normal mechanical laws, i.e., the linear coupling baseline model.

[0191] The linear prediction relationship of this model can be expressed as:

[0192] x_hat_i(t) = a_ij * x_j(t) + b_ij;

[0193] Where x_hat_i(t) is the theoretically reasonable value of the i-th target parameter at time t based on the model prediction, x_j(t) is the actual recorded value of the j-th associated parameter at time t, a_ij is the fixed slope coefficient of the coupling pair obtained by fitting historical data, and b_ij is the intercept constant obtained by fitting.

[0194] Step 902: Substitute the pure steady-state data into each of the linearly coupled baseline models, and calculate the standardized coupling residuals of each of the physical parameter coupling pairs respectively;

[0195] In the real-time cleaning phase, the system substitutes the clean steady-state data extracted from the current segment into the previously constructed equation for verification. The system uses one parameter of the coupling pair as the input independent variable to calculate the theoretical target value that the other related parameter should achieve. Subsequently, the system compares this theoretical target value with the actual recorded value returned by the sensor. Since the magnitudes of different parameters (e.g., kilonewtons and revolutions per minute) differ significantly, the system divides the calculated absolute difference by the residual standard deviation statistically obtained during the historical modeling period of the coupling pair, thereby completing the dimensionless transformation of the difference.

[0196] The corresponding calculation formula is:

[0197] z_cp_ij(t) = (x_i(t) - x_hat_i(t)) / σ_rij;

[0198] Where z_cp_ij(t) is the standardized coupling residual of the calculated output, x_i(t) is the actual recorded value of the i-th target parameter at time t, x_hat_i(t) is the theoretical reasonable value obtained by the above formula, and σ_rij is the standard deviation of the residual of the physical parameter coupling pair in the historical steady-state model.

[0199] The standardized coupling residual objectively quantifies the degree to which the coordination state of this set of physical parameters deviates from the normal mechanical laws at the current moment.

[0200] Step 903: When the absolute value of the standardized coupling residual of the data point at the same time on at least two of the physical parameter coupling pairs is greater than the preset residual threshold, the cross-confirmation mechanism is triggered.

[0201] Exceeding the residual limit for a single coupling pair is insufficient to determine an anomaly, as localized transient geological disturbances may cause temporary inconsistencies in a pair of parameters. However, genuine sensor hardware failures or communication crosstalk will inevitably disrupt all associated mechanical equilibrium relationships simultaneously.

[0202] To establish a rigorous judgment logic, the system performs a global review of data at the same timestamp. If the system finds that the absolute value of the standardized residual of a certain parameter in the first physical parameter coupling pair (such as thrust and earth pressure) is greater than a preset residual threshold, and the absolute value of the standardized residual of the same parameter in a second independently set physical parameter coupling pair (such as total thrust and zone thrust) also exceeds the preset residual threshold, then multiple pieces of evidence overlap. In this case, the system determines that the data exhibits inexplicable contradictions across multiple physical laws and decisively triggers the aforementioned cross-verification mechanism.

[0203] Step 904: The data points that trigger the cross-confirmation mechanism are marked as hidden combination anomalies and incorporated into the system-level noise flags.

[0204] For data points that meet the above cross-validation criteria, although their individual readings may still be within the upper and lower limits allowed by conventional statistics, they disrupt the overall rigid mechanical coordination of the equipment. Based on this, the system determines that the data mutation originates from abnormal physical factors such as local sensor malfunctions. The system formally confirms this type of data as a hidden combined anomaly, classifies and records the corresponding anomaly labels, and incorporates them into the overall system-level noise labeling. During subsequent data import and modeling, the system can accurately isolate these highly concealed system-level deviation data based on these labels.

[0205] Example 10: In a certain scenario, the simplified processing flow is as follows:

[0206] Step S1: Configuration loading and environment initialization.

[0207] Read the core parameters, including Flink batch processing parallelism of 8, batch size of 10,000 records / batch,

[0208] Checkpoint interval of 3 minutes and connection parameters for InfluxDB and MySQL.

[0209] The batch processing environment is built based on Flink BatchExecution Environment, the RocksDBStateBackend storage state is configured, and the Checkpoint mechanism and Exactly-Once semantics are enabled to ensure data consistency after fault recovery.

[0210] High-performance dual-database connections are established using the HikariCP connection pool. The InfluxDB connection pool has a maximum of 20 connections and supports 3 timeout retries (5-second timeout). The MySQL connection pool has a maximum of 15 connections and is configured with a 10-minute query cache.

[0211] Step S2: Define and load cleaning rules.

[0212] In MySQL, a table named t_field_mapping is created, which adopts a three-dimensional structure of "original field name - standard field name - adaptation identifier". Special prefix fields are associated through "prefix matching + semantic similarity" and standardization is achieved in two steps.

[0213] Furthermore, the field mapping function is: ;

[0214] The semantic similarity verification threshold is set to 0.9. For example, G100 and dx100 are eventually unified into "Thrust_100" after mapping.

[0215] Create a table named t_data_type in MySQL and define the "field feature - candidate type - confidence threshold" rule, which supports four core types.

[0216] Furthermore, the data type inference function is:

[0217]

[0218] The first N=800 data points are used as samples. The type is determined when the confidence level is ≥90%; otherwise, manual confirmation is triggered or the data is marked as STRING by default.

[0219] Create a table named t_noise_filter in MySQL, which contains statistical and spatiotemporal association rules. Numeric fields are filtered using a 3σ threshold.

[0220] Furthermore, the formula for calculating the 3σ threshold is:

[0221] ; ;

[0222] Where μ is the field mean, σ is the standard deviation, and data outside the range are marked as suspected noise.

[0223] Create a table named t_dimension_unify in MySQL, using the International System of Units (SI) as the standard, to store physical quantities and conversion factors.

[0224] Furthermore, the dimensional conversion formula is as follows: ;

[0225] Where Qsrc is the original data, kU is the conversion factor, and QSI is the converted SI unit data.

[0226] Step S3: Data source access and data sharding.

[0227] Data can be read by "time range + measurement" through the InfluxDB batch query API. Custom query conditions are supported. When there is no ending condition, all historical data can be read and new data can be listened for.

[0228] Based on the total data volume and Flink's parallelism of 8, the data is dynamically partitioned into 5000-10000 records per partition, and ring number hash mapping is used to ensure that data with the same ring number belongs to the same partition.

[0229] Parse the sharded data, extract core metadata and standardize its storage, and build a sharded index table for subsequent verification and traceability.

[0230] Step S4: Layered batch cleaning.

[0231] The field mapping rule table is called to achieve field mapping through prefix matching and semantic similarity calculation, special characters are removed, and English camelCase naming convention is uniformly adopted.

[0232] Based on the data type inference rule table, the confidence level is calculated through the type inference function to determine the optimal type and automatically convert it. If the conversion fails, it is marked as "type exception".

[0233] Identify the physical quantity in the field, call the dimension conversion formula to convert non-standard units to SI units, and mark units that cannot be identified as "dimensional anomaly".

[0234] Apply statistical noise rules to filter data that exceeds the 3σ range or has an incorrect string length, and mark it as suspected noise.

[0235] Based on spatiotemporal correlation rules, the continuity of ring numbers and the rationality of timestamps are verified, logically invalid data is removed, and suspected noise is confirmed a second time.

[0236] A three-dimensional scoring model is constructed to quantify the score of each data point.

[0237] Furthermore, the scoring formula is as follows: ;

[0238] Scomp represents the integrity score, Scons represents the consistency score, and Seff represents the validity score. Scores below 60 are marked as low-quality data.

[0239] Step S5: Verification and storage of cleaning results

[0240] Perform batch-level verification, including field integrity, data consistency, and quantity reasonableness (the amount of data after cleaning is ≥ 70% of the original data amount, otherwise an alarm is triggered).

[0241] Using the FlinkInfluxDB Sink component, the data is written to the target database in a "batch write + asynchronous commit" mode, with "one database per project and one table per ring number". The timestamp is the primary index, and the scoring and cleaning status are the Tag fields.

[0242] Non-compliant data is categorized by anomaly type and stored in a MySQL temporary backup database, recording the original data, the cause of the anomaly, and the cleaning log, and supports manual re-cleaning.

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

Claims

1. A batch processing and cleaning method for shield tunneling data based on physical and mechanical constraints, characterized in that, include: Obtain the shield tunneling industry data fragments, extract the actual characteristic value range of each parameter, and match and standardize them with a pre-built physical typical value range library to obtain dimensionally unified data; The working condition threshold is calculated based on the extreme value characteristics of the dimensionless data. The dimensionless data is divided into different physical tunneling working conditions according to the working condition threshold, and pure steady-state data is extracted from the advancement stage within the physical tunneling working condition. Multidimensional physical and mechanical constraint verification is performed on pure steady-state data to obtain system-level noise labels. The multidimensional physical and mechanical constraint verification includes anomaly detection for spatial distribution regions of array sensors, multi-parameter physical coupling constraint verification for heterogeneous parameters, and sensor gradual drift detection for long time series. Calculate a data quality score that includes the physical credibility dimension based on system-level noise labeling; The dimensional unified data carrying system-level noise labels and data quality scores is persistently stored as cleaned data. In multidimensional physical and mechanical constraint verification, anomaly detection in the spatial distribution region of array sensors includes: For array sensor parameters with spatial partitioning structure, calculate the spatial deviation of each partition parameter in the pure steady-state data relative to the mean of all partitions at the same time. By utilizing a pre-constructed historical spatial bias baseline, spatial bias is transformed into standardized spatial residuals; When the absolute value of the standardized spatial residual of a single partition parameter is greater than the preset residual threshold, and the total number of partitions with the absolute value of the standardized spatial residual greater than the preset residual threshold does not exceed one-third of the total number of partitions, it is determined that the spatial isolation criterion is met. Data points that satisfy the spatial isolation criterion are marked as spatial distribution anomalies caused by local sensor faults and incorporated into the system-level noise labeling. In multidimensional physical and mechanical constraint verification, the verification of multi-parameter physical coupling constraints for heterogeneous parameters includes: Obtain the linearly coupled baseline model corresponding to multiple physical parameter coupling pairs prefitted using historical steady-state data; Substitute the pure steady-state data into each linearly coupled baseline model and calculate the standardized coupling residuals for each physical parameter coupling pair. When the absolute value of the standardized coupling residual of the data point at the same time on at least two physical parameter coupling pairs is greater than the preset residual threshold, the cross-confirmation mechanism is triggered. Data points that trigger the cross-confirmation mechanism are marked as hidden combination anomalies and incorporated into the system-level noise flag; In multidimensional physical and mechanical constraint verification, for long-term sensor gradual drift detection, the following are included: Using the tunnel boring machine's excavation ring as the aggregation granularity, the ring-level mean values ​​of each parameter in the pure steady-state data are calculated, and a ring-level mean sequence is constructed. Extract the moving median of the mean sequence of ring levels within a preset historical ring number window as the drift reference benchmark; Calculate the upper and lower cumulative sums and control chart statistics of the ring-level mean series using the drift reference baseline; When the cumulative sum of the upper and lower control charts exceeds the preset control limits, it is determined that the corresponding parameter has undergone gradual drift. The drift initiation loop number is determined by backtracking the cumulative starting point of the control chart. Data points within the corresponding drift interval are marked as gradual drift anomalies and incorporated into the system-level noise label.

2. The method according to claim 1, characterized in that, The actual characteristic value range of each parameter in the shield tunneling industrial data segment is extracted and matched and standardized with a pre-built physical typical value range library to obtain dimensionally unified data, including: Calculate the 90th and 99th percentiles of the target parameters in the shield tunneling industrial data slice, and use the interval formed by the 90th and 99th percentiles as the actual feature value range. Calculate the interval matching degree between the actual feature range and the typical range of each candidate unit in the physical typical range library; Select the candidate unit corresponding to the highest interval matching degree as the initial recognition unit; Calculate the ring-level average jump ratio of the target parameter between adjacent tunneling rings; When the difference between the ring-level mean jump ratio and the preset conversion coefficient ratio between two candidate units is less than the preset tolerance, it is determined that the target parameter has undergone a dimensional change. Based on the initial identification units and judgment results, the corresponding unit conversion coefficients are called to convert the target parameters into international standard units.

3. The method according to claim 1, characterized in that, The working condition threshold is calculated based on the extreme value characteristics of the dimensionless data; according to the working condition threshold, the dimensionless data is divided into different physical tunneling working conditions, including: The median absolute deviation algorithm is used to filter extreme values ​​in the thrust and speed sequences of dimensionless data. The 95th quantile of the filtered thrust sequence and speed sequence were extracted as extreme value reference benchmarks. Multiply the extreme value reference benchmark by a preset scaling factor to obtain the thrust threshold and speed threshold, which are used as the operating condition threshold. When the data satisfies the condition that the thrust is greater than the thrust threshold and the rotational speed is greater than the rotational speed threshold, it is determined to be in the propulsion phase; When the data meets the conditions that the thrust is not greater than the thrust threshold, the rotation speed is not greater than the rotation speed threshold, and there is a valid assembly machine start signal, it is determined to be in the assembly stage; The remaining data is classified as being in a shutdown phase. The advance phase, assembly phase, and shutdown phase constitute the physical tunneling working conditions.

4. The method according to claim 1, characterized in that, Pure steady-state data is extracted from the propulsion phase within the physical tunneling conditions, including: The difference quotient of thrust parameters at both ends of a set sliding window during the propulsion phase is calculated to obtain the smooth rate of change; The moment when the absolute value of the smooth rate of change is not greater than the preset rate of change threshold and the thrust amplitude is greater than the preset ratio of the thrust extreme value reference benchmark during the search and propulsion phase is recorded as the steady state start moment. The moment when the same judgment condition is last met during the search advancement phase is recorded as the steady-state end moment. The data segment between the start and end of steady state is extracted as pure steady-state data to separate the start-up and shutdown transition segments within the propulsion phase.

5. The method according to claim 1, characterized in that, Data quality scores, including a physical credibility dimension, are calculated based on system-level noise labeling, including: The system-level noise markers are analyzed, and the number of array sensor spatial partitions, physical parameter coupling pairs, and gradual drift detection items actually triggered at the current data point are counted. The total number of dynamic verification items is then accumulated. The number of verification items identified as anomalous in the statistical system-level noise markers is counted, and the ratio of the number of anomalous verification items to the total number of dynamic verification items is calculated to obtain the dimensionless penalty ratio. The physical credibility score is calculated based on the dimensionless penalty ratio. The data quality score is obtained by weighting and summing the physical credibility score with the pre-calculated integrity score, consistency score, and validity score according to preset weights.

6. The method according to claim 1, characterized in that, The method runs in a distributed batch processing execution environment and includes the following steps before obtaining the tunnel boring machine industrial data fragments: Extract the statistical baseline model parameters required for cross-loop analysis; among which, the statistical baseline model parameters include the historical spatial deviation baseline, linear coupling baseline model and drift reference benchmark required for multidimensional physical and mechanical constraint verification; The statistical baseline model parameters are stored in the global state cache of the distributed batch processing execution environment; A state broadcasting mechanism is adopted to push the updated statistical baseline model parameters to each parallel processing node so that the multidimensional physical and mechanical constraint verification can be directly called for calculation, avoiding the loading and splicing of original historical data across nodes in the shield tunneling industrial data sharding.

7. The method according to claim 6, characterized in that, The statistical baseline model parameters are pre-built and updated using a cold-start hierarchical degradation strategy based on the number of tunneling loops. The cold-start hierarchical degradation strategy includes: Get the cumulative number of completed rings for the current project; When the cumulative number of completed rings is lower than the first preset ring number threshold, multidimensional physical and mechanical constraint verification is not enabled. When the cumulative number of completed rings is between the first preset ring number threshold and the second preset ring number threshold, the existing historical data is used to fit the statistical baseline model parameters, and the judgment threshold of the multidimensional physical and mechanical constraint verification is increased according to the preset relaxation ratio to perform confidence compensation. Once the cumulative number of completed rings reaches the second preset ring threshold, the statistical baseline model parameters are updated using historical data within the sliding window, and the standard judgment threshold is restored.

Images