Multi-dimensional risk supervision method for bridge and underground pipe gallery project

Abstract

The invention discloses a multi-dimensional risk supervision method for bridge and underground pipe gallery engineering, and particularly relates to the field of underground pipe gallery engineering supervision, and the method comprises the steps: an underground structure slippage risk identification and prediction process based on pressure data driving: firstly, collecting an original pressure signal and screening effective data; establishing a standardized two-dimensional pressure field matrix by using the variation function; local gradient mutation features are recognized through Sobe l convolution, evaluation units are divided to calculate entropy distribution, entropy change abnormity is extracted, and a multi-dimensional phase space trajectory is constructed; the maximum Lyapunov index is further applied to estimate and recognize the stability state, gradient and entropy features are combined to construct a parameter inversion observation set, and a constrained optimization algorithm is adopted to iteratively solve an equivalent friction angle; constructing the time sequence into a time sequence, inputting the time sequence into an ARI MA model, executing slip rolling prediction, finally generating a risk thermodynamic diagram with confidence, and achieving the dynamic evaluation and early warning of the slip trend of the underground pipe gallery.

Description

Technical Field

[0001] The present invention relates to the technical field of underground pipeline corridor engineering supervision, and more specifically, to a multi-dimensional risk supervision method for bridge and underground pipeline corridor engineering. Background Art

[0002] Existing monitoring systems for bridge and underground utility tunnel foundations often focus on vertical settlement as their primary assessment metric, overlooking the complex mechanical interactions at the lateral contact interface between the structure and the surrounding soil. In real-world engineering environments, particularly in areas with weak foundations or collapsible soils, factors such as groundwater level fluctuations, seasonal saturation fluctuations, earthquakes, and load disturbances can lead to dynamic evolution of interfacial friction characteristics, resulting in uneven contact pressure distribution and, consequently, interfacial slip. Current monitoring systems lack real-time perception and evolutionary modeling of the interfacial state between the utility tunnel and the soil, often failing to identify early signs of structural misalignment or deformation caused by the accumulation of interfacial slip. This interfacial mechanical rebalancing process typically manifests as increased pressure concentration, sudden changes in stress at boundary nodes, or gradual separation of localized contact areas. However, these manifestations are often dismissed as error fluctuations or instrument drift in traditional early warning systems focused on vertical settlement, potentially evolving into irreversible horizontal slip.

[0003] Therefore, there is an urgent need to construct a structure-foundation slip risk perception method based on the identification of the evolution law of contact pressure distribution, so as to improve the forward identification capability of non-explicit risk patterns such as progressive lateral slip, and at the same time serve as a reference for safety supervision of bridge and underground pipeline corridor projects.

[0004] In order to solve the above problems, a technical solution is now provided. Summary of the Invention

[0005] In order to overcome the above-mentioned defects of the prior art, an embodiment of the present invention provides a multi-dimensional risk management method for bridge and underground pipeline corridor projects to solve the problems raised in the above-mentioned background technology.

[0006] To achieve the above object, the present invention provides the following technical solutions:

[0007] S1: Collect the original pressure signal sequence and filter the valid signal set to generate discrete pressure data with time and space identification;

[0008] S2: The spatial autocorrelation of the pressure data is calculated using the variogram. When the spatial distribution of the data conforms to the characteristics of the spherical model, a standardized two-dimensional pressure field matrix is ​​generated.

[0009] S3: Perform convolution operations on the two-dimensional pressure field matrix in the horizontal and vertical directions respectively to generate horizontal and vertical gradient component matrices, identify mutation points and generate a two-dimensional gradient mutation feature map;

[0010] S4: Divide the two-dimensional pressure field matrix into several evaluation units, calculate the entropy distribution of each evaluation unit, identify abnormal units with entropy changes, and generate entropy change feature vectors;

[0011] S5: Map the entropy change eigenvector to a multidimensional phase space trajectory, estimate the maximum Lyapunov exponent of the phase space trajectory by fitting the divergence rate of adjacent trajectories through least squares, compare the maximum Lyapunov exponent with a preset stability threshold, and generate an interface stability state identifier based on the comparison result;

[0012] S6: Predefine the friction angle parameter space, generate the observation data set for parameter inversion based on the gradient mutation characteristics and entropy change distribution pattern of the current pressure field matrix, and combine the interface stability state identification to iteratively solve the equivalent friction angle;

[0013] S7: Establish the time variation sequence of the equivalent friction angle to construct the slip prediction input vector. Use the ARIMA model to perform rolling prediction on the input vector to obtain the slip prediction value and generate a slip risk heat map with confidence.

[0014] In a preferred embodiment, in S1, collecting the original pressure signal sequence and screening the valid signal set to generate discrete pressure data with time and space identifiers specifically includes:

[0015] An array of piezoresistive sensors is placed on the outer wall of the underground tunnel to collect the original pressure signal sequence;

[0016] Perform moving average filtering on the original signal, calculate the signal-to-noise ratio index of each node in real time, set the signal-to-noise ratio standard value to filter the valid signal set;

[0017] The effective signal set is bound to the spatial coordinates of the tunnel sensor array, and the timestamp synchronization technology is used to align the pressure data stream to generate discrete pressure data with time and space identification.

[0018] In a preferred embodiment, in S2, the spatial autocorrelation of the pressure data is calculated using a variogram. When the spatial distribution of the data conforms to the characteristics of a spherical model, generating a standardized two-dimensional pressure field matrix specifically includes:

[0019] Map the discrete pressure data into a gridded initial matrix, with missing points marked as null values;

[0020] The spatial autocorrelation of pressure data is calculated using the variogram to verify whether the data distribution conforms to the characteristics of the spherical model. If so, the interpolation process is started. If not, the valid signal set is re-screened.

[0021] Applying ordinary Kriging interpolation algorithm to the null values ​​in the initial matrix, a complete pressure distribution matrix is ​​generated based on the spatial correlation of adjacent sensor nodes;

[0022] The absolute pressure values ​​of the interpolated matrix are converted into relative pressure coefficients in [0, 1], and the linear scaling method is used to eliminate the dimension difference to generate a standardized two-dimensional pressure field matrix.

[0023] In a preferred embodiment, in S3, performing convolution operations on the two-dimensional pressure field matrix in the horizontal and vertical directions respectively to generate horizontal and vertical gradient component matrices, identifying mutation points and generating a two-dimensional gradient mutation feature map specifically includes:

[0024] The Sobel operator is used to perform convolution operations on the two-dimensional pressure field matrix in the horizontal and vertical directions, and the absolute value of the pressure gradient at each node is calculated to generate the horizontal and vertical gradient component matrices.

[0025] Based on the statistical distribution characteristics of the gradient component matrix, the local neighborhood gradient mean and standard deviation are calculated according to the sliding window, and the gradient threshold is set to the mean plus the set w times the standard deviation;

[0026] Compare the gradient component value with the dynamic threshold node by node. When the gradient value in any direction exceeds the threshold, it is marked as a mutation point, and a set of abnormal points with spatial coordinates is generated.

[0027] The set of abnormal points is mapped back to the original pressure field coordinate system, and the gradient intensity information is superimposed to generate a two-dimensional gradient mutation feature map.

[0028] In a preferred embodiment, in S4, dividing the two-dimensional pressure field matrix into a plurality of evaluation units, calculating the entropy distribution of each evaluation unit, identifying abnormal units with entropy change, and generating an entropy change feature vector specifically include:

[0029] Divide the standardized two-dimensional pressure field matrix into evaluation units of n×n node size and establish an index coordinate system of the evaluation units;

[0030] Calculate the probability distribution of the pressure value of each node in each evaluation unit to the total value of the unit, and generate the probability density matrix of the unit;

[0031] Apply the Shannon entropy formula to calculate the discrete entropy value of the pressure distribution unit by unit, output the entropy distribution of the evaluation unit, and complete the initial calculation of the unit entropy value;

[0032] The difference between the current unit entropy value and the baseline entropy value established based on the sliding average of the previous m days is calculated. When the difference between z consecutive samplings exceeds the set threshold, it is marked as an entropy change abnormal unit. The spatial distribution information of the abnormal unit and the entropy change amplitude data are integrated to generate the entropy change feature vector.

[0033] In a preferred embodiment, in S5, the entropy change characteristic vector is mapped to a multidimensional phase space trajectory, the divergence rate of adjacent trajectories is fitted by least squares, the maximum Lyapunov exponent of the phase space trajectory is estimated, the maximum Lyapunov exponent is compared with a preset stability threshold, and based on the comparison result, the interface stability state identifier is generated, which specifically includes:

[0034] The entropy change feature vectors are arranged in time sequence to construct a unit entropy value time series dataset, which is then converted into a multi-dimensional phase space trajectory using the phase space reconstruction method.

[0035] The Rosenstein algorithm is used to estimate the maximum Lyapunov exponent of the phase space trajectory, and the divergence rate of adjacent trajectories is fitted by least squares to output the exponential change curve;

[0036] The current Lyapunov exponent is compared with the preset stability threshold. When the exponent value exceeds the threshold continuously and the duration meets the periodic condition, the instability judgment result and the spatial position information are integrated to generate an interface stability state identification matrix with confidence.

[0037] In a preferred embodiment, in S6, a friction angle parameter space is predefined, and based on the gradient mutation characteristics and entropy change distribution pattern of the current pressure field matrix, an observation data set for parameter inversion is generated in combination with the interface stability state identifier. The iterative solution of the equivalent friction angle specifically includes:

[0038] Based on the interface stability state identification matrix and the pressure field matrix, a finite element model of pressure distribution is established to define the friction angle parameter space;

[0039] Extract the gradient mutation characteristics and entropy change distribution pattern of the current pressure field matrix to generate an observation data set for parameter inversion. The observation data set contains the spatial distribution characteristics of pressure gradient and entropy value.

[0040] The LM optimization algorithm with boundary constraints is applied to minimize the residual sum of squares between the model output and the observed data, and the optimal equivalent friction angle is solved iteratively.

[0041] In a preferred embodiment, in S7, a time variation sequence of the equivalent friction angle is established to construct a slip prediction input vector, an ARIMA model is used to perform a rolling prediction on the input vector to obtain a slip prediction value, and a slip risk heat map with confidence is generated, specifically including:

[0042] Continuously solve and establish the time variation series of equivalent friction angle, calculate the equivalent friction angle attenuation rate series, generate non-stationary time series data sets by timestamp alignment, and construct the slip prediction input vector;

[0043] Perform first-order difference processing on the equivalent friction angle series of the input vector and then input it into the prediction model;

[0044] The autoregressive order of the ARIMA model is determined by applying autocorrelation analysis, and the weight coefficient of the moving average term is dynamically adjusted in combination with the time-varying characteristics of the equivalent friction angle.

[0045] Execute the ARIMA model rolling forecast algorithm to output the slippage forecast value and confidence interval within the future set time, and generate a slippage risk heat map with confidence.

[0046] The technical effects and advantages of the multi-dimensional risk supervision method for bridge and underground pipeline corridor projects of the present invention are as follows:

[0047] A high-precision pressure field is constructed through variogram screening. Convolution is combined with Shannon entropy analysis to extract gradient mutations and identify abnormal areas. The Lyapunov exponent is introduced to discriminate slip evolution trends. A high-dimensional observation dataset is constructed to invert the equivalent friction angle, and its time series is used as ARIMA input to achieve rolling slip prediction, ultimately generating a confidence risk heat map. To address the problem that traditional bridge and underground utility corridor monitoring methods focus solely on vertical settlement and ignore the lateral mechanical evolution of the interface, a slip risk identification and prediction method is proposed that integrates gradient mutation characteristics, entropy change distribution patterns, and nonlinear stability analysis. This method has enhanced capabilities for early identification of interface instability and forward-looking assessment of slip trends. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 This is a schematic diagram of a multi-dimensional risk management method for bridge and underground pipeline corridor projects according to the present invention. DETAILED DESCRIPTION

[0049] The following will provide a clear and complete description of the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0050] Example 1

[0051] Figure 1 The present invention provides a multi-dimensional risk management method for bridge and underground pipeline corridor projects, which includes the following steps:

[0052] S1: Collect the original pressure signal sequence and filter the valid signal set to generate discrete pressure data with time and space identification;

[0053] S2: The spatial autocorrelation of the pressure data is calculated using the variogram. When the spatial distribution of the data conforms to the characteristics of the spherical model, a standardized two-dimensional pressure field matrix is ​​generated.

[0054] S3: Perform convolution operations on the two-dimensional pressure field matrix in the horizontal and vertical directions respectively to generate horizontal and vertical gradient component matrices, identify mutation points and generate a two-dimensional gradient mutation feature map;

[0055] S4: Divide the two-dimensional pressure field matrix into several evaluation units, calculate the entropy distribution of each evaluation unit, identify abnormal units with entropy changes, and generate entropy change feature vectors;

[0056] S5: Map the entropy change eigenvector to a multidimensional phase space trajectory, estimate the maximum Lyapunov exponent of the phase space trajectory by fitting the divergence rate of adjacent trajectories through least squares, compare the maximum Lyapunov exponent with a preset stability threshold, and generate an interface stability state identifier based on the comparison result;

[0057] S6: Predefine the friction angle parameter space, generate the observation data set for parameter inversion based on the gradient mutation characteristics and entropy change distribution pattern of the current pressure field matrix, and combine the interface stability state identification to iteratively solve the equivalent friction angle;

[0058] S7: Establish the time variation sequence of the equivalent friction angle to construct the slip prediction input vector. Use the ARIMA model to perform rolling prediction on the input vector to obtain the slip prediction value and generate a slip risk heat map with confidence.

[0059] In S1, the original pressure signal sequence is collected and the valid signal set is screened to generate discrete pressure data with time and space identification.

[0060] An array of piezoresistive sensors is evenly distributed across the exterior of the underground tunnel. Each sensor has a resolution of 0.1 kPa and a data sampling frequency of 1 Hz. They are spaced 0.5 meters apart, covering the tunnel floor, sidewalls, and roof, ensuring a full pressure sensing range with no blind spots. The raw pressure signals collected by the sensors are continuous time series, with each record accompanied by the original timestamp and corresponding sensor number.

[0061] The raw pressure signal was smoothed using a moving average filter with a sliding window length of 10 seconds to eliminate short-term fluctuations. The signal-to-noise ratio (SNR) within each window was calculated for each node. This ratio is defined as the ratio of the signal mean to the standard deviation. The SNR threshold was set at 2.5. If the SNR exceeded the threshold for three consecutive windows, the node data was considered valid. Invalid signals were marked and excluded from subsequent modeling.

[0062] The resulting valid signal set is bound to the three-dimensional spatial coordinates of the corresponding sensor to generate a space-time-value triplet data structure. Timestamps are aligned using a unified time base, with an error controlled within ±0.2 seconds, ensuring temporal synchronization of data streams across multiple nodes. The resulting discrete pressure dataset is organized in a matrix format, with each row representing the complete spatial pressure distribution at a point in time. This serves as the basic input for subsequent spatial field reconstruction and dynamic evolution analysis.

[0063] In S2, the spatial autocorrelation of pressure data is calculated using the variogram, and a standardized two-dimensional pressure field matrix is ​​generated when the spatial distribution of the data conforms to the characteristics of the spherical model.

[0064] The discrete pressure data is mapped to a two-dimensional grid matrix according to the spatial distribution of the sensors. The grid spacing is consistent with the sensor layout, and missing data points are marked as null values ​​to form the initial pressure field matrix. The spatial pressure difference is modeled using a variogram, which is expressed as:

[0065]

[0066] Where, γ(h) is the variation value when the spacing is h, Z(x i ) and Z(x i +h) are the positions x i The number of data pairs calculated near the distance h and the pressure data value N(h) at the distance h.

[0067] The spherical model is a common variogram fitting model used to describe the spatial autocorrelation characteristics of variables as they change with distance. The variograms in each direction are calculated and fitted to determine whether they conform to the characteristics of the spherical model, that is, γ(h) tends to increase first and then tend to saturation as the distance h increases.

[0068] If the variogram distribution satisfies a spherical structure, the interpolation process is initiated. Otherwise, the sensor node density and effective signal set coverage are reassessed, and nodes are removed or compensated before resampling. When interpolation is initiated, Ordinary Kriging is used to interpolate the empty points in the initial matrix. Based on the spatial autocovariance matrix and the weighted minimum variance estimation method, data from adjacent effective nodes are used for inference.

[0069] After taking the absolute value of the pressure value in the two-dimensional matrix completed by Kriging interpolation, linear normalization processing is performed to standardize the value and map all pressure data to the interval [0,1] to eliminate the original pressure unit difference and sensor calibration deviation. Finally, a standardized two-dimensional pressure field matrix is ​​generated as the input for subsequent spatial gradient feature extraction and dynamic stability analysis.

[0070] In S3, convolution operations are performed on the two-dimensional pressure field matrix in the horizontal and vertical directions respectively to generate horizontal and vertical gradient component matrices, identify mutation points and generate a two-dimensional gradient mutation feature map.

[0071] The input is the normalized two-dimensional pressure field matrix P'(x,y) generated in the previous step. The Sobel operator is used to perform convolution operations in the horizontal (x direction) and vertical (y direction) directions to extract local gradient changes. The following Sobel kernel is used:

[0072] Gx=[[-1,0,1],[-2,0,2],[-1,0,1]],

[0073] Gy=[[-1,-2,-1],[0,0,0],[1,2,1]].

[0074] Perform two-dimensional convolution operations on the matrix P' to obtain the horizontal gradient matrix Mx(x,y) and the vertical gradient matrix My(x,y), both of which are in absolute value form to represent the gradient strength.

[0075] By default, a 5×5 pixel (where the pixel is the specific pressure data corresponding to the x, y coordinates) gradient intensity sliding window is used to perform local statistics on Mx and My, calculate the neighborhood mean and standard deviation point by point, and construct a dynamic threshold model.

[0076] The gradient anomaly recognition thresholds Tx = μ(x) + wσ(x) and Ty = μ(y) + wσ(y) are set at each pixel position to achieve adaptive mutation point recognition, where Tx and Ty are the gradient anomaly recognition thresholds in the horizontal and vertical directions, μ(x) and μ(y) are the means of the window in the horizontal and vertical directions, σ(x) and σ(y) are the standard deviations of the window in the horizontal and vertical directions, and w is the set standard deviation multiple. The specific value is set according to the sliding window size and the actual gradient change, and the default value is 2.5.

[0077] Traverse each coordinate point in the entire matrix. If Mx(x,y)>Tx(x,y) or My(x,y)>Ty(x,y), mark the point as an abnormal mutation point. Record its two-dimensional spatial coordinates and the maximum gradient strength value G(x,y)=max(Mx,My), and construct the abnormal point set E={(x,y,G)}.

[0078] Initialize the mutation feature map with the same size as the original matrix, map the coordinates of all abnormal points and corresponding gradient values ​​in the set E to the feature map, and assign the value of non-abnormal points to 0. The final generated image reflects the local mutation intensity in spatial grayscale, which serves as one of the input features for subsequent stability analysis, entropy change distribution identification, and risk heat map construction.

[0079] In S4, the two-dimensional pressure field matrix is ​​divided into several evaluation units, the entropy distribution of each evaluation unit is calculated, the abnormal units with entropy change are identified, and the entropy change feature vector is generated.

[0080] The standardized two-dimensional pressure field matrix is ​​divided into evaluation units of size n×n (n is greater than 0 and is consistent with the size of the gradient intensity sliding window). During the division process, an evaluation unit index coordinate system is established, and the numbering method adopts row-first ascending index to ensure spatial reversible mapping.

[0081] The sum of the node pressure value sets in each unit is calculated, and the probability density vector is calculated according to the node pressure ratio to generate the probability density matrix in the unit.

[0082] Apply the Shannon entropy formula to calculate the discrete entropy value of the pressure distribution unit by unit, output the entropy distribution of the evaluation unit, and complete the initial calculation of the unit entropy value. The calculation expression is:

[0083]

[0084] Where H(k) is the unit entropy value of the kth unit, p k (j) is the probability density of the jth node of the kth unit. k (j) Items that are 0 are excluded or minimized to avoid numerical errors. Each unit gets an entropy value to form a two-dimensional entropy distribution matrix.

[0085] A time series entropy map is constructed, recording the historical entropy values ​​of each unit over the continuous monitoring period. A baseline entropy mean is calculated over a sliding window over the past m days (specified by the data update frequency and the rate of entropy change, with an initial default of 30). The real-time entropy value is then differentiated from the baseline entropy value. If the entropy exceeds a set threshold (specified by the desired monitoring sensitivity, with a default of 25% of the baseline) for z consecutive monitoring times (z ≥ 3), the unit is marked as an entropy change anomaly. The location, entropy difference amplitude, and current normalized pressure center value of all units marked as anomaly are integrated to construct an entropy change feature vector of a spatial-numerical composite structure, which is used as input for subsequent phase space reconstruction and nonlinear stability discrimination models.

[0086] In S5, the entropy change eigenvector is mapped to the multidimensional phase space trajectory. The divergence rate of adjacent trajectories is fitted by least squares to estimate the maximum Lyapunov exponent of the phase space trajectory. The maximum Lyapunov exponent is compared with the preset stability threshold. Based on the comparison results, the interface stability state identifier is generated.

[0087] The extracted entropy change feature vectors are classified and arranged in time series according to the evaluation unit index to form an independent entropy value time series for each unit. An entropy value time series set is established, and the multidimensional trajectory vector set {V(t)} is reconstructed using a delayed embedding method.

[0088] The maximum Lyapunov exponent (λ) describes how adjacent trajectories diverge over time in phase space:

[0089] λ>0: The system is very sensitive to initial conditions and exhibits chaotic characteristics, such as earthquakes, fluid turbulence, and structural vibrations;

[0090] λ≈0: The system trajectory remains essentially unchanged, indicating quasi-periodic motion, such as a pendulum;

[0091] λ<0: All trajectories eventually converge to a stable state, such as a damped oscillating system.

[0092] In the structural health monitoring of bridge-underground pipeline corridors, we use the Lyapunov exponent to determine whether the system vibration enters a nonlinear chaotic state, so as to monitor possible structural damage or abnormal load effects.

[0093] The Rosenstein algorithm is used to estimate the maximum Lyapunov exponent. The specific process is as follows:

[0094] Find the nearest neighbor trajectory for each trajectory point in the phase space, track the distance change between the two over time, and calculate the distance change of each pair of nearest neighbor points over time t:

[0095] D(t)=∥V(t+Δt)-V(t near +Δt)‖;

[0096] Where D(t) is the distance expansion change between the two, V(t) is any vector in the multidimensional trajectory vector set {V(t)}, t is the time index of the current main trajectory, V(t near ) is the nearest trajectory point vector of vector V(t), t near The time index of the nearest trajectory point, Δt is the minimum change time step, which is set based on the data update frequency.

[0097] The natural logarithm of the distance over time is sampled and the least squares straight line fitting is performed. The fitting formula is:

[0098] lnD(t)=λ·t+C;

[0099] Where λ is the maximum Lyapunov exponent (slope) we estimate, C is the intercept, which represents the logarithm of the initial distance, and the final output is the exponential change curve λ(t).

[0100] A stability threshold is set (specifically set based on the slope change rate, with a default setting of 0.15). If the Lyapunov exponent λ(t) exceeds this threshold for a set number of consecutive time periods, and the duration continuously exceeds the minimum period constraint (specifically set based on the data update frequency, with a default setting of 10 minutes), the corresponding unit is considered to be in a dynamic instability trend. The spatial index of the unstable unit, the duration of the threshold exceeding the limit, and the rate of change of the Lyapunov curve are integrated into the interface stability state data, and the output is a two-dimensional state identification matrix with confidence scores. The confidence levels are represented by 0 and 1, with 0 indicating no instability and 1 indicating instability.

[0101] In S6, the friction angle parameter space is predefined, and based on the gradient mutation characteristics and entropy change distribution pattern of the current pressure field matrix and the interface stability state identification, an observation data set for parameter inversion is generated, and the equivalent friction angle is solved iteratively.

[0102] The generated interface stability state identification matrix and the standardized two-dimensional pressure field matrix were input, and a two-dimensional elastic mechanics finite element model was established based on this to simulate the contact behavior between the underground pipeline corridor and the foundation soil. The model was divided using quadrilateral plane stress elements, and the boundary conditions set the bottom of the pipeline corridor to be fixed and the side walls to allow horizontal elastic slip. Coulomb friction constraints were introduced on the contact surface to define the equivalent friction angle parameter space. The step size is 0.5°, which is used to construct the friction angle-pressure field response mapping.

[0103] From the current normalized pressure matrix, we extract the gradient mutation signature and entropy change distribution maps for the regions corresponding to the anomaly stability markers, forming the observation dataset O(x,y) = {G(x,y),H(x,y)}. G is the local gradient magnitude, and H is the local entropy. The dataset coverage is limited to regions with a confidence score of 1, and the first-order neighboring domain is expanded to enhance inversion stability.

[0104] Initialize the friction angle search starting value Run finite element forward simulation based on the current model parameters and calculate the simulated pressure field matrix Compute the residuals between the simulated and observed data.

[0105] Applying the Levenberg-Marquardt optimization algorithm with boundary constraints, The friction angle estimate is iteratively updated within the range. The iteration process is terminated when the rate of decrease of the residual sum of squares is less than the set step size ε (specifically set according to the inversion residual sum of squares, with a default setting of 0.01) or when the maximum number of iterations exceeds the maximum set number of iterations (set according to the time constraint of the iterative calculation), and the current optimal friction angle estimate is output.

[0106] In S7, the time variation sequence of the equivalent friction angle is established to construct the slip prediction input vector. The ARIMA model is used to perform rolling prediction on the input vector to obtain the slip prediction value and generate a slip risk heat map with confidence.

[0107] The time variation series of the equivalent friction angle is continuously solved and established, the equivalent friction angle attenuation rate series is calculated, the non-stationary time series data set is generated by timestamp alignment, and the slip prediction input vector is constructed.

[0108] The ADF stationarity test is performed on the original series of equivalent friction angle. After confirming that it is non-stationary, the first-order difference processing is performed as the premise of ARIMA modeling.

[0109] The autoregressive order p and moving average order q of the ARIMA model, as well as the difference order d, are determined by the autocorrelation function ACF and partial autocorrelation function PACF graphs as follows:

[0110] If the original sequence is stationary, d = 0; if the first-order difference is stationary, d = 1 (1 is selected by default based on the characteristics of the slip data);

[0111] The ACF function is used to select q. If the ACF plot truncates quickly after the qth order (clearly drops to zero), the truncation position is used as the recommended value of q.

[0112] If the PACF plot is truncated quickly after the pth order, the position of the truncation is used as the recommended value of p.

[0113] Combined with the fluctuation characteristics of the friction angle series, the weight coefficient of the moving average term is dynamically adjusted to adapt to the time-varying response of the system. Finally, p, d, and q are brought into the determined ARIMA (p, d, q) model structure (when the data update frequency is low and the slip trend is relatively stable, the simplified model order ARIMA (1, 1, 1) can be used).

[0114] Based on the final model, the 95% confidence interval of the slip prediction value is calculated, and a spatial slip risk distribution map is constructed. Combined with the prediction credibility annotation, a risk heat map is generated for the risk supervision system to call and geographical visualization display.

[0115] The above formulas are all dimensionless and numerical calculations. The formulas are obtained by collecting a large amount of data and performing software simulation to obtain the most recent real situation. The preset parameters and thresholds in the formulas are set by technicians in this field according to actual conditions.

[0116] The above embodiments can be implemented in whole or in part by software, hardware, firmware or any other combination. When implemented using software, the above embodiments can be implemented in whole or in part in the form of a computer program product. The computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer program are loaded or executed on a computer, the process or function described in the embodiment of the present application is generated in whole or in part. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another computer-readable storage medium. For example, the computer instructions can be transmitted from one website, computer, server or data center to another website, computer, server or data center via a wired (e.g., infrared, wireless, microwave, etc.) method. The computer-readable storage medium can be any available medium that can be accessed by a computer or a data storage device such as a server or data center that contains one or more available media sets. The available medium can be a magnetic medium (e.g., a floppy disk, a hard disk, a tape), an optical medium (e.g., a DVD), or a semiconductor medium. The semiconductor medium can be a solid-state drive.

[0117] Those skilled in the art will appreciate that the modules and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0118] Those skilled in the art will clearly understand that, for the convenience and brevity of description, the specific working processes of the systems, devices and modules described above can refer to the corresponding processes in the aforementioned method embodiments and will not be repeated here.

[0119] In the several embodiments provided in this application, it should be understood that the disclosed systems, devices and methods can be implemented in other ways. For example, the device embodiments described above are merely schematic. For example, the division of the modules is only a logical function division. In actual implementation, there may be other division methods, such as multiple modules or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of devices or modules, which can be electrical, mechanical or other forms.

[0120] The modules described as separate components may or may not be physically separate, and the components shown as modules may or may not be physical modules, and may be located in one place or distributed across multiple network modules. Some or all of the modules may be selected to achieve the purpose of this embodiment according to actual needs.

[0121] In addition, each functional module in each embodiment of the present application may be integrated into one processing module, or each module may exist physically separately, or two or more modules may be integrated into one module.

[0122] If the functions are implemented in the form of software function modules and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present application, or the part that contributes to the prior art or the part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present application. The aforementioned storage medium includes various media that can store program codes, such as a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk.

[0123] The above description is merely a specific embodiment of the present application, but the scope of protection of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in this application should be included in the scope of protection of this application. Therefore, the scope of protection of this application should be based on the scope of protection of the claims.

[0124] Finally: The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A multi-dimensional risk supervision method for bridge and underground pipeline corridor projects, characterized in that: The steps include: S1: Collect the original pressure signal sequence and filter the valid signal set to generate discrete pressure data with time and space identification; S2: The spatial autocorrelation of the pressure data is calculated using the variogram. When the spatial distribution of the data conforms to the characteristics of the spherical model, a standardized two-dimensional pressure field matrix is ​​generated. S3: Perform convolution operations on the two-dimensional pressure field matrix in the horizontal and vertical directions respectively to generate horizontal and vertical gradient component matrices, identify mutation points and generate a two-dimensional gradient mutation feature map; S4: Divide the two-dimensional pressure field matrix into several evaluation units, calculate the entropy distribution of each evaluation unit, identify abnormal units with entropy changes, and generate entropy change feature vectors; S5: Map the entropy change eigenvector to a multidimensional phase space trajectory, estimate the maximum Lyapunov exponent of the phase space trajectory by fitting the divergence rate of adjacent trajectories through least squares, compare the maximum Lyapunov exponent with a preset stability threshold, and generate an interface stability state identifier based on the comparison result; S6: Predefine the friction angle parameter space, generate the observation data set for parameter inversion based on the gradient mutation characteristics and entropy change distribution pattern of the current pressure field matrix, and combine the interface stability state identification to iteratively solve the equivalent friction angle; S7: Establish the time variation sequence of the equivalent friction angle to construct the slip prediction input vector. Use the ARIMA model to perform rolling prediction on the input vector to obtain the slip prediction value and generate a slip risk heat map with confidence.

2. A multi-dimensional risk management method for bridge and underground pipeline corridor projects according to claim 1, characterized in that: In S1, the original pressure signal sequence is collected and the valid signal set is screened to generate discrete pressure data with time and space identification, which specifically includes: An array of piezoresistive sensors is placed on the outer wall of the underground tunnel to collect the original pressure signal sequence; Perform moving average filtering on the original signal, calculate the signal-to-noise ratio index of each node in real time, set the signal-to-noise ratio standard value to filter the valid signal set; The effective signal set is bound to the spatial coordinates of the tunnel sensor array, and the timestamp synchronization technology is used to align the pressure data stream to generate discrete pressure data with time and space identification.

3. A multi-dimensional risk management method for bridge and underground pipeline corridor projects according to claim 2, characterized in that: In S2, the spatial autocorrelation of pressure data is calculated using the variogram. When the spatial distribution of the data conforms to the characteristics of the spherical model, a standardized two-dimensional pressure field matrix is ​​generated, specifically including: Map the discrete pressure data into a gridded initial matrix, with missing points marked as null values; The spatial autocorrelation of pressure data is calculated using the variogram to verify whether the data distribution conforms to the characteristics of the spherical model. If so, the interpolation process is started. If not, the valid signal set is re-screened. Applying ordinary Kriging interpolation algorithm to the null values ​​in the initial matrix, a complete pressure distribution matrix is ​​generated based on the spatial correlation of adjacent sensor nodes; The absolute pressure values ​​of the interpolated matrix are converted into relative pressure coefficients in [0, 1], and the linear scaling method is used to eliminate the dimension difference to generate a standardized two-dimensional pressure field matrix.

4. A multi-dimensional risk management method for bridge and underground pipeline corridor projects according to claim 3, characterized in that: In S3, convolution operations are performed on the two-dimensional pressure field matrix in the horizontal and vertical directions to generate horizontal and vertical gradient component matrices, identify mutation points, and generate a two-dimensional gradient mutation feature map. Specifically, the following steps are performed: The Sobel operator is used to perform convolution operations on the two-dimensional pressure field matrix in the horizontal and vertical directions, and the absolute value of the pressure gradient at each node is calculated to generate the horizontal and vertical gradient component matrices. Based on the statistical distribution characteristics of the gradient component matrix, the local neighborhood gradient mean and standard deviation are calculated according to the sliding window, and the gradient threshold is set to the mean plus the set w times the standard deviation; Compare the gradient component value with the dynamic threshold node by node. When the gradient value in any direction exceeds the threshold, it is marked as a mutation point, and a set of abnormal points with spatial coordinates is generated. The set of abnormal points is mapped back to the original pressure field coordinate system, and the gradient intensity information is superimposed to generate a two-dimensional gradient mutation feature map.

5. A multi-dimensional risk management method for bridge and underground pipeline corridor projects according to claim 4, characterized in that: In S4, the two-dimensional pressure field matrix is ​​divided into several evaluation units, the entropy distribution of each evaluation unit is calculated, the abnormal units with entropy change are identified, and the entropy change feature vector is generated. Specifically, the following steps are performed: Divide the standardized two-dimensional pressure field matrix into evaluation units of n×n node size and establish an index coordinate system of the evaluation units; Calculate the probability distribution of the pressure value of each node in each evaluation unit to the total value of the unit, and generate the probability density matrix of the unit; Apply the Shannon entropy formula to calculate the discrete entropy value of the pressure distribution unit by unit, output the entropy distribution of the evaluation unit, and complete the initial calculation of the unit entropy value; The difference between the current unit entropy value and the baseline entropy value established based on the sliding average of the previous m days is calculated. When the difference between z consecutive samplings exceeds the set threshold, it is marked as an entropy change abnormal unit. The spatial distribution information of the abnormal unit and the entropy change amplitude data are integrated to generate the entropy change feature vector.

6. A multi-dimensional risk management method for bridge and underground pipeline corridor projects according to claim 5, characterized in that: In S5, the entropy change eigenvector is mapped to the multidimensional phase space trajectory. The divergence rate of adjacent trajectories is fitted by least squares to estimate the maximum Lyapunov exponent of the phase space trajectory. The maximum Lyapunov exponent is compared with the preset stability threshold. Based on the comparison results, the interface stability state identification is generated, which specifically includes: The entropy change feature vectors are arranged in time sequence to construct a unit entropy value time series dataset, which is then converted into a multi-dimensional phase space trajectory using the phase space reconstruction method. The Rosenstein algorithm is used to estimate the maximum Lyapunov exponent of the phase space trajectory, and the divergence rate of adjacent trajectories is fitted by least squares to output the exponential change curve; The current Lyapunov exponent is compared with the preset stability threshold. When the exponent value exceeds the threshold continuously and the duration meets the periodic condition, the instability judgment result and the spatial position information are integrated to generate an interface stability state identification matrix with confidence.

7. A multi-dimensional risk management method for bridge and underground pipeline corridor projects according to claim 6, characterized in that: In S6, the friction angle parameter space is predefined. Based on the gradient mutation characteristics and entropy change distribution pattern of the current pressure field matrix and the interface stability state identification, the observation data set for parameter inversion is generated. The iterative solution of the equivalent friction angle specifically includes: Based on the interface stability state identification matrix and the pressure field matrix, a finite element model of pressure distribution is established to define the friction angle parameter space; Extract the gradient mutation characteristics and entropy change distribution pattern of the current pressure field matrix to generate an observation data set for parameter inversion. The observation data set contains the spatial distribution characteristics of pressure gradient and entropy value. The LM optimization algorithm with boundary constraints is applied to minimize the residual sum of squares between the model output and the observed data, and the optimal equivalent friction angle is solved iteratively.

8. A multi-dimensional risk management method for bridge and underground pipeline corridor projects according to claim 7, characterized in that: In S7, the time variation sequence of the equivalent friction angle is established to construct the slip prediction input vector. The ARIMA model is used to perform rolling prediction on the input vector to obtain the slip prediction value and generate a slip risk heat map with confidence. Specifically, the following steps are performed: Continuously solve and establish the time variation series of equivalent friction angle, calculate the equivalent friction angle attenuation rate series, generate non-stationary time series data sets by timestamp alignment, and construct the slip prediction input vector; Perform first-order difference processing on the equivalent friction angle series of the input vector and then input it into the prediction model; The autoregressive order of the ARIMA model is determined by applying autocorrelation analysis, and the weight coefficient of the moving average term is dynamically adjusted in combination with the time-varying characteristics of the equivalent friction angle. Execute the ARIMA model rolling forecast algorithm to output the slippage forecast value and confidence interval within the future set time, and generate a slippage risk heat map with confidence.

Images