Excavation face stability judgment method and device based on soil bin pressure gradient field distribution

Abstract

The invention discloses an excavation face stability judgment method and device based on soil bin pressure gradient field distribution, and the method relates to the field of constructional engineering.The method comprises the steps that an initial soil bin pressure plane gradient field distribution model is built through a non-uniform cubic B-spline basis function; taking a prediction error of the initial soil bin pressure plane gradient field distribution model as a target function, and determining an optimal control point; obtaining a soil bin pressure plane gradient field distribution model according to the optimal control point; obtaining an excavation face stable region of the pressure gradient field in the soil bin; according to the soil bin pressure plane gradient field distribution model, the maximum value and the minimum value of the gradient change of the soil bin pressure field are determined; if the maximum value and the minimum value are both in the excavation face stability domain, it is determined that the excavation face of the soil bin is stable, and if the maximum value or the minimum value is not in the excavation face stability domain, it is determined that the excavation face is unstable. According to the method, the stability of the excavation face in the shield construction process can be judged in real time.

Description

Technical Field

[0001] The present invention relates to the field of construction engineering, and particularly to the determination of the stability of the excavation face based on the distribution of the soil bin pressure gradient field. Background Art

[0002] One of the key indicators for evaluating the performance of shield tunneling is the settlement control index above the tunnel. Especially during the construction of urban subways, the settlement control index is highly concerned during the shield tunneling process, and the stability of the shield excavation face is the key factor determining the settlement control index.

[0003] Currently, the stability of the shield excavation face focuses on the determination of the excavation face failure mode and the support pressure of the excavation face. The key to controlling the stability of the shield excavation face is to balance the water and soil pressure in front of the shield and the pressure inside the soil bin. According to the instability failure mechanism of the excavation face, microscopic and mesoscopic analysis models, plastic limit theory analysis methods, wedge mechanics analysis models, etc. have emerged. The above theories have all achieved certain research results, but they cannot be well applied to guide the actual on-site construction, and cannot judge the stability of the excavation face in real time to guide the shield operators at the construction site to operate in a timely manner.

[0004] Therefore, there is an urgent need for a method that can determine the stability of the excavation face during the shield tunneling process in real time. Summary of the Invention

[0005] Based on this, in view of the above technical problems, it is necessary to provide a method and device for determining the stability of the excavation face based on the distribution of the soil bin pressure gradient field, which can determine the stability of the excavation face during the shield tunneling process in real time.

[0006] The present invention adopts the following technical solutions:

[0007] The present invention provides a method for determining the stability of the excavation face based on the distribution of the soil bin pressure gradient field, including:

[0008] Constructing an initial soil bin pressure plane gradient field distribution model with non-uniform cubic B-spline basis functions; the initial soil bin pressure plane gradient field distribution model includes a function of the soil pressure value in the soil bin changing with position;

[0009] Taking the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function, optimizing the control points of the initial soil bin pressure plane gradient field distribution model according to the soil pressure values at the soil bin monitoring points to determine the optimal control points; the control points are the coefficients of the non-uniform cubic B-spline basis functions;

[0010] Taking the optimal control points as the parameters in the initial soil bin pressure plane gradient field distribution model to obtain the soil bin pressure plane gradient field distribution model;

[0011] According to the variation range of the pressure gradient in the vertical and horizontal directions inside the soil bin, the stable region of the excavation face of the pressure gradient field inside the soil bin is obtained;

[0012] According to the distribution model of the plane pressure gradient field of the soil bin, determine the maximum and minimum values of the gradient change of the soil bin pressure field;

[0013] If both the maximum and minimum values are within the stable region of the excavation face, it is determined that the excavation face of the soil bin is stable. If there is a maximum or minimum value not within the stable region of the excavation face, it is determined that the excavation face is unstable.

[0014] Preferably, the initial distribution model of the plane pressure gradient field of the soil bin is:

[0015]

[0016] Among them, p(x, y, t) is a non-uniform cubic B-spline basis function, B i,α (x) is the B-spline basis function defined by the knot vector x, B j,β (y) is the B-spline basis function defined by the knot vector y, d ij is the control point, m is the number of control vertex boundaries in the x direction, n is the number of control vertex boundaries in the y direction, i is the random value of the control vertex in the x direction, and j is the random value of the control vertex in the y direction.

[0017] Preferably, the soil bin monitoring points include the existing soil bin monitoring points and the newly added soil bin monitoring points, specifically including:

[0018] Obtain the soil pressure value and angle of the existing soil bin monitoring points;

[0019] Select the angle at equal intervals between any two existing soil bin monitoring points as the angle of the newly added soil bin monitoring points;

[0020] Create a periodic cubic spline interpolation model according to the soil pressure values of the existing soil bin monitoring points;

[0021] According to the periodic cubic spline interpolation model, calculate the soil pressure value corresponding to the angle of the newly added soil bin monitoring points;

[0022] Take the existing soil bin monitoring points and the newly added soil bin monitoring points as the soil bin monitoring points.

[0023] Preferably, the method further includes:

[0024] For any newly added soil bin monitoring point, take the newly added soil bin monitoring point as the test set, and the remaining soil bin monitoring points as the training set; the test set includes the soil pressure test value of the newly added soil bin monitoring point; the training set includes the soil pressure test values of the remaining soil bin monitoring points;

[0025] Construct a cubic spline interpolation model through the training set;

[0026] According to the cubic spline interpolation model and the test set, the predicted values of the soil bin pressure for the test set are obtained;

[0027] Calculate the mean square error between the predicted values of the soil bin pressure and the measured values of the soil bin pressure;

[0028] If the mean square error is within the preset error threshold range, it is determined that the setting of the new pressure monitoring point in the soil bin is reasonable.

[0029] Preferably, taking the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function, and optimizing the control points of the initial soil bin pressure plane gradient field distribution model according to the soil pressure values of the soil bin monitoring points to determine the optimal control points, specifically including:

[0030] Taking the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function and the control points as variables, optimize the control points through the chaotic adaptive particle swarm - sequential quadratic programming algorithm, and take the control points corresponding to the minimum objective function value in the iterative process of the chaotic adaptive particle swarm - sequential quadratic programming algorithm as the optimal control points.

[0031] Preferably, taking the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function, and optimizing the control points of the initial soil bin pressure plane gradient field distribution model according to the soil pressure values of the soil bin monitoring points to determine the optimal control points, specifically further including:

[0032] The objective function is:

[0033]

[0034] where, min x,y (p(x k , y k , t) - p k (t)) is the prediction error of the initial soil bin pressure plane gradient field distribution model, p k (t) is the true value of the soil bin pressure, and p(x, y, t) is the predicted value of the soil bin pressure;

[0035] Convert the optimal solution of the objective function into solving the optimal control points, and define the first function as:

[0036]

[0037] where, B i,α (x) is the B - spline basis function defined by the knot vector x, B j,β (y) is the B - spline basis function defined by the knot vector y, d ij is the control point, p k (t) is the true value of the soil bin pressure, x kLet \(x\) be the random optimization value in the \(x\)-direction, \(i\) be the random value in the \(x\)-direction, \(j\) be the random value in the \(y\)-direction, \(m\) be the number of control vertex bounds in the \(x\)-direction, \(n\) be the number of control vertex bounds in the \(y\)-direction, \(k\) be the number of control vertex optimization times, and \(r\) be the total number of control vertex optimization times.

[0038] When the first function is minimized, we can obtain

[0039] Substitute into the first function, and the first equation is obtained as:

[0040]

[0041] Simplify the first equation to obtain the first matrix expression as:

[0042] B T BD(t) = B T P(t);

[0043] where B T is the transpose of the B-spline basis function matrix, B is the B-spline basis function matrix, D(t) is the control vertex value at time t, and P(t) is the soil bin pressure value at time t;

[0044] Use the singular value decomposition method to solve the first matrix expression, and the control points are obtained as:

[0045]

[0046] where D(t) is the control vertex value at time t, U is an orthogonal matrix, V T is the transpose of the orthogonal matrix, P(t) is the soil bin pressure value at time t, is the matrix singular value.

[0047] Calculate the singular value of matrix B and substitute the singular value into the control points to obtain the optimal control points.

[0048] Preferably, according to the pressure gradient change ranges in the vertical and horizontal directions inside the soil bin, the excavation face stability region of the soil bin internal pressure gradient field is obtained, specifically including:

[0049] Determine the minimum value of the internal pressure change in the soil bin according to the minimum value of the pressure gradient change range in the vertical direction inside the soil bin and the minimum value of the pressure gradient change range in the horizontal direction inside the soil bin;

[0050] Determine the maximum value of the internal pressure change in the soil bin according to the maximum value of the pressure gradient change range in the vertical direction inside the soil bin and the maximum value of the pressure gradient change range in the horizontal direction inside the soil bin;

[0051] Determine the stable area of the excavation face of the soil pressure gradient field in the soil bin according to the minimum and maximum values of the internal pressure change in the soil bin.

[0052] Preferably, the range of the pressure gradient change in the vertical direction inside the soil bin is:

[0053]

[0054] Wherein, is the range of the pressure gradient change in the vertical direction inside the soil bin, ρ m is the density of the muck inside the soil bin, τ a is the cohesion of the muck inside the soil bin, L is the vertical distance of the soil bin, and g is the acceleration due to gravity;

[0055] The range of the pressure gradient change in the horizontal direction inside the soil bin is:

[0056]

[0057] Wherein, is the range of the pressure gradient change in the horizontal direction inside the soil bin;

[0058] The calculation formulas for the minimum and maximum values of the internal pressure change in the soil bin are:

[0059]

[0060] Wherein, is the value of the internal pressure change in the soil bin, is the value of the pressure gradient change in the vertical direction inside the soil bin, is the value of the pressure gradient change in the horizontal direction inside the soil bin.

[0061] Preferably, the soil bin pressure plane gradient field distribution model constructs multiple B-spline curves based on the same non-uniform cubic B-spline basis function, and the soil bin pressure values on each B-spline curve are the same.

[0062] The present invention provides a device for determining the stability of the excavation face based on the distribution of the soil bin pressure gradient field, which is characterized by including:

[0063] An initial model construction module, which is used to construct an initial soil bin pressure plane gradient field distribution model with a non-uniform cubic B-spline basis function; the initial soil bin pressure plane gradient field distribution model includes a function of the soil pressure value in the soil bin changing with position;

[0064] A first determination module, which is used to take the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function, optimize the control points of the initial soil bin pressure plane gradient field distribution model according to the soil pressure values of the soil bin monitoring points, and determine the optimal control points; the control points are the coefficients of the non-uniform cubic B-spline basis function;

[0065] A substitution module, configured to use the optimal control points as parameters in the initial soil bin pressure plane gradient field distribution model to obtain the soil bin pressure plane gradient field distribution model;

[0066] A second determination module, configured to obtain the excavation face stability region of the soil bin internal pressure gradient field according to the pressure gradient change ranges in the vertical and horizontal directions inside the soil bin;

[0067] A third determination module, configured to determine the maximum and minimum values of the soil bin pressure field gradient change according to the soil bin pressure plane gradient field distribution model;

[0068] A determination module, configured to determine that the excavation face of the soil bin is stable if both the maximum value and the minimum value are within the excavation face stability region, and determine that the excavation face is unstable if there is a maximum value or a minimum value not within the excavation face stability region.

[0069] The present invention provides a computer-readable storage medium storing a computer program, and when the computer program is executed by a processor, the above-described method for determining the stability of an excavation face based on the distribution of a soil bin pressure gradient field is implemented.

[0070] The present invention provides a computer device, including a memory, a processor, and a computer program stored on the memory and executable on the processor, and when the processor executes the program, the above-described method for determining the stability of an excavation face based on the distribution of a soil bin pressure gradient field is implemented.

[0071] At least one of the above technical solutions adopted by the present invention can achieve the following beneficial effects:

[0072] The present invention constructs an initial soil bin pressure plane gradient field distribution model based on non-uniform cubic B-spline basis functions, providing a mathematical basis for constructing and representing complex soil bin pressure environments. It has higher approximation ability and better set characteristics, and can more accurately depict the changes in the soil bin pressure gradient field. The prediction error of the initial soil bin pressure plane gradient field distribution model is used as the objective function, and the control points of the initial soil bin pressure plane gradient field distribution model are optimized according to the soil pressure values at the soil bin monitoring points to determine the optimal control points. The optimal control points are used as parameters in the initial soil bin pressure plane gradient field distribution model to obtain the soil bin pressure plane gradient field distribution model, making the prediction accuracy of the soil bin pressure plane gradient field distribution model higher by determining the optimal control points. According to the pressure gradient change ranges in the vertical and horizontal directions inside the soil bin, the excavation face stability region of the soil bin internal pressure gradient field is obtained. According to the soil bin pressure plane gradient field distribution model, the maximum and minimum values of the soil bin pressure field gradient change are determined. If both the maximum and minimum values are within the excavation face stability region, it is determined that the excavation face of the soil bin is stable; if there is a maximum or minimum value not within the excavation face stability region, it is determined that the excavation face is unstable. This method can determine the stability of the excavation face during shield tunneling construction in real time. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] The drawings described herein are used to provide a further understanding of the present invention and form a part of the present invention. The schematic embodiments and descriptions thereof are used to explain the present invention and do not constitute an improper limitation to the present invention. In the drawings:

[0074] Figure 1 B-spline curves of different degrees obtained from the same control polygon provided by the present invention;

[0075] Figure 2 Schematic flow chart of the excavation face stability determination method based on the soil bin pressure gradient field distribution provided by the present invention;

[0076] Figure 3 Distribution map of shield soil bin monitoring points provided by the present invention;

[0077] Figure 4 Right-line shield soil bin pressure parameters provided by the present invention;

[0078] Figure 5 Plan view of the shield tunnel section provided by the present invention;

[0079] Figure 6 Schematic longitudinal section of the shield tunnel geology provided by the present invention;

[0080] Figure 7 Shield soil bin pressure numerical distribution map provided by the present invention;

[0081] Figure 8The fitting result of the soil bin pressure gradient field provided by the present invention;

[0082] Figure 9 The ideal and actual soil bin pressure distribution diagram provided by the present invention;

[0083] Figure 10 A schematic diagram of an excavation surface stability determination device based on soil bin pressure gradient field distribution provided by the present invention;

[0084] Figure 11 A schematic diagram of a computer device for implementing a method for determining excavation surface stability based on soil bin pressure gradient field distribution provided by the present invention. DETAILED DESCRIPTION

[0085] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention will be clearly and completely described below in conjunction with the specific embodiments of the present invention and the corresponding drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0086] B-spline basis functions constitute a set of special polynomial functions defined on the knot vectors. These knot vectors consist of a non-decreasing sequence of knot coordinates in the parameter space, providing a mathematical basis for constructing and representing complex shapes. The knot vector definition is shown in formula (1):

[0087] E=[ξ1,ξ2,…,ξ n+k+1 ] (1)

[0088] Among them, ξ i ∈R is the i-th node, i is the node vector sequence number, k is the number of basis functions defined on the node vector E, n is the number of basis functions defined on the node vector E, [ξ i ,ξ i+1 ] is the node interval.

[0089] If the lengths of the intervals between all adjacent nodes are equal, it is called a uniform node, otherwise it is a non-uniform node.

[0090] The method of defining B-spline basis function usually uses recursive definition method or non-recursive definition method. The present invention uses "truncated power function difference quotient" to define B-spline basis function, which is a non-recursive definition method. The value of the basis function can be directly calculated without step-by-step recursive calculation. Therefore, it has more advantages in computing efficiency, flexibility and generalization ability.

[0091] First, we introduce the truncated power function The definition of , the truncated power function is shown in formula (2):

[0092]

[0093] Among them, \(t\) is the position of the calculation node, \(x\) is the truncation position, and \(k\) is the power of the function.

[0094] As can be seen from formula (2), only the part where \(t\geq x\) in \((t - x)\) k is intercepted, and the value of the part where \(t\lt x\) is zero.

[0095] The B-spline basis function defined by the difference quotient of the truncated power function requires the use of the difference quotient operation formula, which is shown in formula (3), and the truncated power function identity is shown in formula (4):

[0096] [x i ,x i+1 ,…,x i+k+1 (t - x) k =0 (3)

[0097] [x i ,x i+1 ,…,x i+k+1 (t - x) k =0 (4)

[0098] Thus, the \(k\)-th B-spline basis function can be defined by the difference quotient of the truncated power function, as shown in formula (5):

[0099]

[0100] If the basis function is defined on a uniform knot vector, it is called a normalized B-spline basis function. Therefore, the definition of the \(k\)-th normalized B-spline basis function is shown in formula (6):

[0101] B i,j (x)=(x i+k+1 -x i )M i,j (x) (6)

[0102] For the convenience of understanding and subsequent calculations, the 0-th B-spline basis function and the 1-st B-spline basis function will be introduced below.

[0103] In the calculation process of the 0-th B-spline basis function, first substitute \(k = 0\) into formula (5), and the result is shown in formula (7):

[0104]

[0105] Thus, the expression of the 0-th B-spline basis function is shown in formula (8):

[0106]

[0107] Similarly, combining formula (6), the expression of the 0th-order normalized B-spline basis function can be obtained as shown in formula (9):

[0108]

[0109] The calculation process of the 1st-order B-spline basis function is the same as that of the 0th-order. Substituting k = 1 into formula (5), we can get as shown in formula (10):

[0110]

[0111] Expanding and calculating formula (10), we can get as shown in formula (11):

[0112]

[0113] Combining the definition of the truncated power function in formula (2), formula (11) can be further simplified to obtain as shown in formulas (12) and (13):

[0114] M i,1 (x) = 0, x > x i+2 (12)

[0115]

[0116] Substituting the divided-difference operation formula of formula (3) and the truncated power function identity of formula (4) into the calculation expression of the 1st-order B-spline basis function in formula (10) respectively, we can obtain the 1st-order B-spline basis function under the condition of the divided-difference definition of the truncated power function, as shown in formula (14):

[0117]

[0118] Further inductive arrangement can obtain the piecewise expression of the 1st-order B-spline basis function, as shown in formula (15):

[0119]

[0120] Combining formula (6), the expression of the 1st-order normalized B-spline basis function can be obtained again, as shown in formula (16):

[0121]

[0122] Similar to the above derivation process, other high-order B-spline basis functions can be obtained. Since high-order B-spline basis functions have higher approximation ability and better geometric properties, they can more accurately depict the change of the earth pressure gradient field. However, higher-order B-spline basis functions will cause a relatively large amount of computation, a complex calculation process, and are prone to overfitting. Therefore, combined with the research results of relevant references, this paper will select cubic B-spline basis functions for calculation, and the calculation results are expressed as shown in formula (17):

[0123]

[0124] The following gives the basic properties and applicable conditions of B-spline basis functions:

[0125] (1) Normativity: It refers to the standardization and consistency in definition and use, specifically including the consistency of mathematical definitions, the standardization of data structures and representations, the unity of calculation algorithms, and the interoperability of applications and tools. The general expression of the normativity of B-spline basis functions is shown in formula (18).

[0126]

[0127] (2) Positivity and local support: Positivity means that the B-spline basis function is non-negative within its function domain, that is, the value of each function is non-negative, ensuring that the linear combination of functions is always non-negative, so that the generated curve is reasonable in mathematical logic and physical meaning; local support means that the B-spline basis function has only a finite non-zero interval within the domain, that is, each basis function has a value only in the local interval of its domain, and takes zero value in other intervals. This can improve the editing efficiency and flexibility of spline curves. The calculation expressions of positivity and local support are as follows: exactly has s different zeros within (x j , x j+n ) and is always zero outside [x j , x j+n . In particular, there is a relationship expressed as formula (19):

[0128]

[0129] (3) Recursiveness: High-order B-spline basis functions are defined recursively, which not only simplifies the calculation of high-order functions but also makes it easier to change the properties of B-spline curves during the design and editing process, such as adjusting the smoothness or shape of the curve. Recursiveness is an important property for B-spline basis functions to flexibly and efficiently generate curves of different orders. The recursive relationship can be expressed by formula (20) and formula (21):

[0130]

[0131] (4) Differentiability: It characterizes the property that the B-spline basis functions are continuously differentiable of all orders within their domain, which is beneficial for optimizing and controlling the appearance and properties of curves. Its computational expression can be described by formulas (22) and (23):

[0132] B′ j,1 (x) = 0, x ≠ x j , x ≠ x j+1 (22)

[0133]

[0134] It should be noted that only when n = 2, formula (23) can hold at non-knot points.

[0135] On the knot vector in a specific parameter direction, the spline basis functions of corresponding orders can be calculated through recursive formulas. The linear combination of these spline basis functions forms a B-spline curve, where the coefficients of each basis function are called control points. When these control points are set as a series of coordinate points, a B-spline curve in space can be constructed.

[0136] The figure formed by connecting each control point in sequence by straight lines is called a control polygon. By introducing the knot vector and n k-th order B-spline basis functions, combined with the corresponding control vertices, the accurate expression of the B-spline curve is realized. As shown in formula (24):

[0137]

[0138] Among them, P(ξ) is the B-spline curve, B i,k (ξ) is the n k-th order B-spline basis functions, P i is the control vertex.

[0139] Under the framework of formula (24), each basis function corresponds to a control vertex respectively, ensuring that the geometric shape of the B-spline curve can vary flexibly through the adjustment of the control vertices.

[0140] Since the parameter coordinate ξ takes values at the knots ξ = ξ i , and its value is not equal to 1, the B-spline curve does not interpolate through the control vertices. Figure 1 shows the fitting effects of B-spline curves of different orders based on the same knot vector and defined by the same control polygon.

[0141] Such as Figure 1As shown, B-spline curves of all orders interpolate through control points at the start and end. At the same time, the curve is tangent to the control polygon at these endpoints, and this property stems from the definition of B-spline curves based on open knot vectors. It should be noted that as the order of the curve increases, its degree of freedom increases, resulting in a greater deviation between the curve and the control points. In Figure 1 The quadratic and cubic B-spline curves shown in

[0142] exhibit high continuity in the regions outside the endpoints, while at the endpoints, the curve only ensures basic continuity. ij (i = 0, 1, ..., m; j = 0, 1, ..., n) of the control points and the control grid (m + 1) × (n + 1), combined with the parameters u and v and the degrees k and l, and using the knot vectors U = [u0, u1, ..., u m+k+1 and V = [ν0, ν1, ..., ν n+l+1 , the shape and smoothness of the curve can be precisely controlled, thus achieving highly flexible geometric modeling. The control point matrix plays a crucial role in defining the overall structure of the curve, and the selection of parameters and the configuration of knot vectors directly affect the detailed performance of the curve. By finely adjusting these elements, the B-spline curve equation of degree k × l can be obtained as shown in formula (25):

[0143]

[0144] where B i,k (u) is the B-spline basis function defined by the knot vector U, B j,l (v) is the B-spline basis function defined by the knot vector V, u e ≤ u < u e+1 , ν f ≤ ν < ν f+1 is the sub-rectangular domain.

[0145] The local properties of B-spline curves provide a theoretical basis for surface generalization, enabling the construction of B-spline surfaces with high flexibility and accuracy through the control points d e ≤ u < u e+1 , ν f ≤ ν < ν f+1 defined within the sub-rectangular domain. By modifying some vertices, only the local area of the surface is affected, without causing extensive changes to the entire surface. Then the above surface equation can be rewritten as shown in formula (26): ij

[0146]

[0147] ​The non-uniform B-spline basis function technology demonstrates excellent performance in the control of the overall surface shape and the ability of local regulation, providing an efficient method for data fitting. Through the application of the least squares method, this technology can find a balance between global adaptability and accurate capture of local details, optimizing the matching degree between the surface and data points.

[0148] The following will, in conjunction with the accompanying drawings, elaborate on the technical solutions provided by the embodiments of the present invention.

[0149] Figure 2 The following is a schematic flow chart of a method for determining the stability of an excavation face based on the distribution of the soil bin pressure gradient field in the present invention, specifically including the following steps:

[0150] S201: Construct an initial soil bin pressure plane gradient field distribution model using non-uniform cubic B-spline basis functions; the initial soil bin pressure plane gradient field distribution model includes a function of the soil pressure value in the soil bin changing with position.

[0151] Collect the real-time monitoring data set of pressure sensors at each position in the soil bin as L = {(x k , y k , p k ) | k = 1, 2,... r}. Let a = min k∈{1,2,...,r} x k ; b = max k∈{1.2....,r} x k ; c = min k∈{1,2,...,r} y k ; d = max k∈{1,2,...,r} y k . According to the sensor position distribution (x k , y k ), the soil bin pressure plane can be divided by the straight lines x = u i , y = v j , where, u i = a + ih x ; ν j = c + jh y ; If the B-spline basis function S(x, y) in the defined region R is α, β times respectively along the x, y directions, then its parametric knot vectors are shown in formulas (27), (28):

[0152] U = [u0, u1,..., u m+α+1 (27)

[0153] V = [v0, v1,..., v n+β+1 (28)

[0154] where, u0 = a, u m+α+1 = b, ν0 = c, ν n+β+1 = d.

[0155] Therefore, the B-spline functions B i.α (u) and B i,β (ν) are calculated and expressed as shown in Formulas (29) and (30) respectively:

[0156]

[0157] In an exemplary embodiment, under the condition of the domain R, the initial earth pressure plane gradient field distribution model is as shown in Formula (31):

[0158]

[0159] where p(x, y, t) is a non-uniform cubic B-spline basis function, B iα (x) is the B-spline basis function defined by the knot vector x, and B j,β (y) is the B-spline basis function defined by the knot vector y, is, m is the number of control vertex boundaries in the x direction, n is the number of control vertex boundaries in the y direction, i is the random value of the control vertex in the x direction, and j is the random value of the control vertex in the y direction.

[0160] S202: Use the prediction error of the initial earth pressure plane gradient field distribution model as the objective function, optimize the control points of the initial earth pressure plane gradient field distribution model according to the earth pressure values of the earth pressure monitoring points in the earth bin, and determine the optimal control points; the control points are the coefficients of the non-uniform cubic B-spline basis function.

[0161] In an exemplary embodiment, the earth pressure monitoring points in the earth bin include the existing earth pressure monitoring points and the newly added earth pressure monitoring points in the earth bin, specifically including: obtaining the earth pressure values and angles of the existing earth pressure monitoring points in the earth bin; equally spacing and selecting angles between any two existing earth pressure monitoring points in the earth bin as the angles of the newly added earth pressure monitoring points in the earth bin; creating a periodic cubic spline interpolation model according to the earth pressure values of the existing earth pressure monitoring points in the earth bin; calculating the earth pressure values corresponding to the angles of the newly added earth pressure monitoring points in the earth bin according to the periodic cubic spline interpolation model; using the existing earth pressure monitoring points and the newly added earth pressure monitoring points in the earth bin as the earth pressure monitoring points in the earth bin.

[0162] Specifically, adding new monitoring points in the earth pressure balance chamber can be achieved through code. Based on the CubicSpline function in the SciPy library, periodic cubic spline interpolation is implemented. The input is the angles of the existing monitoring points in the earth pressure balance chamber and the corresponding earth pressure values. By setting the boundary condition as 'periodic', it is ensured that the interpolation results are smooth and continuous at the starting and ending points, meeting the periodic boundary conditions. New angles are evenly selected between the existing monitoring points in the earth pressure balance chamber, and the earth pressure values corresponding to these angles are calculated using the interpolation model, which are the predicted earth pressure values of the newly added monitoring points. Based on the code operation results, the specific earth pressure values of the newly added monitoring points E7 - E18 in the earth pressure balance chamber can be obtained, as shown in Table 1, which is the statistical result of the pressure values of the newly added monitoring points in the earth pressure balance chamber.

[0163] Table 1

[0164]

[0165] Specifically, when there are relatively few existing monitoring points in the shield earth pressure balance chamber, the cubic spline interpolation method is used to fit the earth pressure field, effectively increasing the number of monitoring points in the earth pressure balance chamber. By inserting new earth pressure monitoring points between the existing monitoring points in the earth pressure balance chamber, the understanding of the earth pressure distribution in the earth pressure balance chamber is enhanced, providing support for accurately monitoring and controlling the shield earth pressure balance chamber pressure. The distribution of the earth pressure balance chamber pressure monitoring points is as Figure 3 shown, where E1 - E6 are the existing monitoring points in the earth pressure balance chamber, and E7 - E18 are the newly added monitoring points in the earth pressure balance chamber.

[0166] The process of adding new monitoring points in the earth pressure balance chamber using the cubic spline interpolation method is introduced in detail as follows:

[0167] The cubic spline interpolation function depends on piecewise-defined cubic polynomials to approximate the path between data points. Each polynomial is valid within its corresponding data point interval, ensuring not only the continuity of the curve at the data points but also the continuity of the first and second derivatives, thus generating a smooth and natural transition curve. The expression relationship of the cubic spline interpolation polynomial is shown in Equation (32):

[0168] S(x i )=y i (i = 1, 2, …, n - 1) (32)

[0169] where, S(x i ) is the cubic spline interpolation polynomial in the i-th interval, and y i is the polynomial coefficient.

[0170] The piecewise cubic interpolation polynomial expression of S(x) is shown in Equation (33):

[0171]

[0172] Among them, x is an internal node, and S(x) is an interpolation polynomial.

[0173] Based on the natural boundary, fixed boundary, and periodic boundary, the cubic spline interpolation function effectively realizes the smooth transition between data points through piecewise two-point cubic Hermite interpolation polynomials. The interpolation function within each segment depends not only on the endpoint values and first-order derivatives but also involves the calculation expressions of the second-order derivatives, as shown in Formulas (34) and (35):

[0174]

[0175] Among them, h k is the difference between the (k + 1)-th and k-th nodes, and m k - the first derivative value at the k-th node.

[0176] In Formula (35), let k = k - 1, and divide Formula (36) by Formula (35). Formula (36) is as follows:

[0177]

[0178] After simplification, Formula (37) can be obtained:

[0179] λ k m k-1 + 2m k + μ k m k+1 = g k (k = 1, 2, …, n - 1) (37)

[0180] Among them, each physical quantity satisfies the equal relationship in Formulas (38) to (40):

[0181]

[0182] In the construction process of the cubic spline interpolation function, the traversal of the k value involves solving a system of n - 1 equations to determine n + 1 unknowns, where the unknowns represent the first derivative values of the interpolation function at three adjacent nodes. The establishment of this basic system of equations depends on the continuity and smoothness requirements of the interpolation function at each node, ensuring that the curve has second-order continuous smooth characteristics within each segment and at the segment joints, thus realizing high-quality and smooth transition between data points. The following is a detailed explanation of the three types of boundaries.

[0183] The natural boundary condition requires that the second derivative is equal to zero at both ends of the curve. This means that the curve is free at both ends and is not restricted by additional slopes or shapes, so that the curve can transition as flat as possible near the endpoints. It is necessary to satisfy m0 = f′0, m n = f′ n .

[0184] Under fixed boundary conditions, the first derivative of the curve at both ends is set to a specific value. This allows more control over the start and end of the curve. It is necessary to satisfy S”(x0)=f0″,S”(x n ) = f n ″, in the subintervals [x0,x1] and [x n-1 ,x n ] on m0, m1, m n-1 With m n In formula (34), let k = 0, x = x0, and we can get the following equation:

[0185]

[0186] Similarly, in formula (34), let k = n-1, x = x n , the calculation can be obtained as shown in formulas (42) and (43):

[0187]

[0188] Where f' is the first-order derivative of f(x), f n is the nth-order derivative of f(x).

[0189] Combining equations (41) to (43) with the basic equation (37), we can obtain equation (44):

[0190]

[0191] The periodic boundary condition requires that the starting point and the end point of the curve have not only the same value, but also the same first-order and second-order derivatives, thus ensuring the smooth periodic continuation of the curve, which needs to be satisfied as shown in formula (45):

[0192] μ n m1+λ n m m-1 +2m n =g n (45)

[0193] Among them, m n By replacing m0 and combining formula (37) and formula (45), we can get the cubic matrix expression, as shown in formula (46):

[0194]

[0195] Based on the above theoretical analysis, taking the actual data of the right-line shield tunnel as an example, Figure 4As shown, a total of 394 rings of data values of the earth pressure monitoring points E1 - E6 are collected. It should be noted that for each ring of earth pressure monitoring values, the box method is used to process the data, and the average value of the processed data is obtained. As shown in Table 2, Table 2 shows the earth pressure values of E1 - E6.

[0196] Table 2

[0197]

[0198] To simplify the interpolation calculation, based on Figure 3 the distribution of the earth pressure monitoring points, it can be approximately considered that the newly added earth pressure monitoring points are annularly distributed on the circumference of the shield cross-section, that is, they meet the periodic boundary requirements of formula (45). According to the given earth pressure values and periodic boundary conditions, a linear equation system is constructed and solved to obtain the coefficients of each cubic polynomial. Using the obtained cubic polynomial function, the earth pressure values of the newly added earth pressure monitoring points are calculated. The earth pressure values of the newly added earth pressure monitoring points are shown in Table 3.

[0199] Table 3

[0200]

[0201] In an exemplary embodiment, the method further includes:

[0202] For any newly added earth pressure monitoring point, the newly added earth pressure monitoring point is used as a test set, and the remaining earth pressure monitoring points are used as a training set; the test set includes the earth pressure test values of the newly added earth pressure monitoring points; the training set includes the earth pressure test values of the remaining earth pressure monitoring points; a cubic spline interpolation model is constructed through the training set; according to the cubic spline interpolation model and the test set, the earth pressure prediction value of the test set is obtained; the mean square error between the earth pressure prediction value and the earth pressure test value is calculated; if the mean square error is within the preset error threshold range, it is determined that the setting of the newly added earth pressure monitoring point is reasonable.

[0203] Specifically, to evaluate the rationality of the values of the newly added earth pressure monitoring points, based on the characteristic of fewer earth pressure monitoring points, the Leave One Out Cross Validation (LOOCV) is selected for evaluation. The basic idea of the Leave One Out Cross Validation is to leave out one of all the monitoring points as a test set each time, and the rest as a training set, which is used to construct a cubic spline interpolation model, then predict the value of the left-out point, and finally evaluate the accuracy of all predictions. The Leave One Out Cross Validation is implemented through code, and the Mean Square Error (MSE) is used as the performance evaluation index.

[0204] In an exemplary embodiment, the prediction error of the initial soil bin pressure plane gradient field distribution model is used as the objective function, and the control points of the initial soil bin pressure plane gradient field distribution model are optimized according to the soil pressure values at the soil bin monitoring points to determine the optimal control points. Specifically, it includes: using the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function and the control points as variables, optimizing the control points through the chaotic adaptive particle swarm - sequential quadratic programming algorithm, and taking the control points corresponding to the minimum objective function value during the iterative process of the chaotic adaptive particle swarm - sequential quadratic programming algorithm as the optimal control points.

[0205] Specifically, the difference between the soil pressure value predicted by the soil bin pressure plane gradient field distribution model and the true soil pressure value is used as the objective function for optimization. With the minimum objective function value as the optimization goal, the chaotic adaptive particle swarm - sequential quadratic programming algorithm is used to find the control point values, so that the objective function value is minimized. The control points corresponding to the minimum objective function value are the optimal control points.

[0206] In an exemplary embodiment, the prediction error of the initial soil bin pressure plane gradient field distribution model is used as the objective function, and the control points of the initial soil bin pressure plane gradient field distribution model are optimized according to the soil pressure values at the soil bin monitoring points to determine the optimal control points. Specifically, it further includes:

[0207] The objective function is shown in formula (47):

[0208]

[0209] where min x,y (p(x k , y k , t) - p k (t)) is the prediction error of the initial soil bin pressure plane gradient field distribution model, p k (t) is the true value of the soil bin pressure, p(x, y, t) is the predicted value of the soil bin pressure, i is the random value in the x - direction, j is the random value in the y - direction, k is the number of optimization times of the control vertices, and d(t) is the optimization value of the control vertices at time t.

[0210] Converting the optimal solution of the objective function to solve the optimal control points, defining the first function as shown in formula (48):

[0211]

[0212] where B i,α (x) is the B - spline basis function defined by the knot vector x, B j,β (y) is the B - spline basis function defined by the knot vector y, d ij is the control point, p k (t) is the true value of the soil bin pressure, xk is the random optimization value in the x direction, i is the random value in the x direction, j is the random value in the y direction, m is the number of control vertex boundaries in the x direction, n is the number of control vertex boundaries in the y direction, k is the number of optimization times of the control vertex, and r is the total number of optimization times of the control vertex.

[0213] Based on formula (47), the objective of the solution is to find a set of control vertices d ij (t) such that the value of the function f(d ij (t)) is minimized. For this purpose, equalities are introduced as shown in formula (49) and formula (52):

[0214]

[0215] p(t) = [p1(t), p2(t),..., p r (t)] (52)

[0216] When the first function is minimized, the following can be obtained as shown in formula (53):

[0217]

[0218] Substitute formula (53) into formula (48) to obtain the following as shown in formula (54):

[0219]

[0220] Simplify formula (54) to obtain the following as shown in formula (55):

[0221] B T BD(t) = B T P(t) (55)

[0222] where B T is the transpose of the B-spline basis function matrix, B is the B-spline basis function matrix, D(t) is the control vertex value at time t, and P(t) is the soil bin pressure value at time t.

[0223] Since the shield soil bin pressure data points are scattered, the real-time monitoring data are all discrete data, which cannot fully ensure that each divided node area contains at least one data, and may cause the matrix values in a certain row of formula (55) to be all zero, resulting in a singular situation, and further causing the equation to be unsolvable. In view of this, the present invention introduces the method of singular value decomposition to solve formula (55).

[0224] Perform singular value decomposition on the matrix B. There exist orthogonal matrices U and V of order m, and the following equality can be established as shown in formula (56):

[0225] B = USV T (56)

[0226] Among them, S r = diag(σ1, σ2,..., σ r ), λ i is the eigenvalue of B T B.

[0227] Then the expression solution of the optimal control vertex d ij (t) is shown in formula (57) as follows:

[0228]

[0229] Among them, D(t) is the numerical value of the control vertex at time t, U is an orthogonal matrix, V T is the transpose of the orthogonal matrix, P(t) is the numerical value of the earth pressure in the hopper at time t, is the singular value of the matrix.

[0230] S203: Use the optimal control point as a parameter in the initial earth pressure plane gradient field distribution model to obtain the earth pressure plane gradient field distribution model.

[0231] Specifically, as can be seen from formula (57), by calculating the singular value σ i of matrix B, it can be substituted into formula (57) for solution, and the unique solution of the control vertex can be obtained. Substituting the control vertex solution into formula (58) can obtain the characteristic distribution model of the earth pressure field. Formula (58) is as follows:

[0232]

[0233] In an exemplary embodiment, the earth pressure plane gradient field distribution model constructs multiple B-spline curves based on the same non-uniform cubic B-spline basis function, and the earth pressure values on each B-spline curve are the same.

[0234] Specifically, the earth pressure plane gradient field distribution model includes multiple B-spline curves, and the pressing force on each B-spline curve is the same, which can clearly mark the pressure distribution state inside the hopper.

[0235] S204: According to the pressure gradient change ranges in the vertical and horizontal directions inside the hopper, obtain the excavation surface stability domain of the earth pressure gradient field inside the hopper.

[0236] In an exemplary embodiment, according to the variation ranges of the pressure gradients in the vertical and horizontal directions inside the soil bin, the excavation face stability domain of the pressure gradient field inside the soil bin is obtained, specifically including: determining the minimum value of the pressure change inside the soil bin according to the minimum value of the variation range of the pressure gradient in the vertical direction inside the soil bin and the minimum value of the variation range of the pressure gradient in the horizontal direction inside the soil bin; determining the maximum value of the pressure change inside the soil bin according to the maximum value of the variation range of the pressure gradient in the vertical direction inside the soil bin and the maximum value of the variation range of the pressure gradient in the horizontal direction inside the soil bin; and determining the excavation face stability domain of the pressure gradient field inside the soil bin according to the minimum and maximum values of the pressure change inside the soil bin.

[0237] Specifically, obtain the variation ranges of the pressure gradients in the numerical and horizontal directions inside the soil bin, determine the minimum value of the pressure change inside the soil bin according to the minimum values of the variation ranges of the pressure gradients in the vertical and horizontal directions, determine the maximum value of the pressure change inside the soil bin according to the maximum values of the variation ranges of the pressure gradients in the vertical and horizontal directions, and determine the excavation face stability domain of the pressure gradient field inside the soil bin according to the minimum and maximum values of the pressure change inside the soil bin.

[0238] In an exemplary embodiment, the variation range of the pressure gradient in the vertical direction inside the soil bin is as shown in formula (59):

[0239]

[0240] Where, is the variation range of the pressure gradient in the vertical direction inside the soil bin, ρ m is the density of the muck inside the soil bin, τ a is the cohesion of the muck inside the soil bin, L is the vertical distance of the soil bin, and g is the acceleration due to gravity;

[0241] The variation range of the pressure gradient in the horizontal direction inside the soil bin is as shown in formula (60):

[0242]

[0243] Where, is the variation range of the pressure gradient in the horizontal direction inside the soil bin;

[0244] The calculation formulas for the minimum and maximum values of the pressure change inside the soil bin are as shown in formula (61):

[0245]

[0246] Where, is the value of the pressure change inside the soil bin, is the value of the pressure gradient change in the vertical direction inside the soil bin, is the value of the pressure gradient change in the horizontal direction inside the soil bin.

[0247] Substitute Equation (59) and Equation (60) into Equation (61), and the expression of the stable region of the excavation face of the internal pressure gradient field in the shield soil bin can be obtained, as shown in Equation (62):

[0248]

[0249] S205: Determine the maximum and minimum values of the gradient change of the soil bin pressure field according to the distribution model of the soil bin pressure plane gradient field.

[0250] The calculation expressions for the maximum and minimum values of the gradient change of the soil bin pressure plane field are shown in Equations (63) and (64) respectively:

[0251]

[0252] Among them,

[0253] Solve Equations (63) and (64) using the chaotic adaptive particle swarm - sequential quadratic programming algorithm.

[0254] Since the objective functions in Equations (63) and (64) have non - smooth characteristics and belong to constrained optimization problems, it is difficult to solve them using general numerical calculation methods. Therefore, a chaotic adaptive particle swarm algorithm is introduced to solve the problem of finding the optimal value of spatial optimization. At the same time, the sequential quadratic programming algorithm is used to transform the constrained optimization problem into an unconstrained optimization problem.

[0255] The particle swarm optimization algorithm (PSO) is an optimization technique based on swarm intelligence. It solves optimization problems by simulating the social behavior of bird flocks or fish schools. Each "particle" represents a potential solution in the problem space. The search process is guided by tracking and updating the individual and global optimal solutions. Assume that the dimension of the particle search space is D, the number of particles is N, and x i is the current position of particle i, and v i is the current speed of particle i. Then the evolution equations of the particle swarm algorithm are shown in Equations (65) - (68):

[0256]

[0257] Among them, t is the number of iterations, ω is the inertia weight, which is used to control the influence of the current speed of the particle on the future speed, is the best position experienced by particle i, is the best position experienced by all particles, c1 and c2 are learning factors, which are used to control the tendency of the particle to move towards the individual best position and the global best position. rand1 and rand2 are independent random numbers in the range [0, 1], in order to increase the randomness of the search.

[0258] The sequential quadratic programming method demonstrates its unique efficiency and accuracy in dealing with nonlinear optimization problems, especially constrained optimization problems. By transforming the nonlinear original problem into a series of quadratic programming sub-problems, this method realizes the gradual approximation to the optimal solution. Each sub-problem is a quadratic approximation of the original problem at the current iteration point, and this transformation makes the solution of the nonlinear problem more feasible and efficient. The descriptions of the original problem and the sub-problem are shown in equations (69) and (70) respectively:

[0259]

[0260] where f(x) is the objective function, g(x) is the gradient vector of the original function, W(x, λ, μ) is the Hessian matrix, J E (x) is the Jacobi matrix of the equality constraint C E (x), and J I (x) is the Jacobi matrix of the inequality constraint C I (x).

[0261] The logical chaos search algorithm is an optimization method based on chaos theory, which uses the nonlinear dynamic characteristics of the chaos system to enhance the diversity of the search process and avoid falling into local optimal solutions. The chaos system has characteristics such as determinism, randomness, and periodicity, and can search globally. Therefore, the search algorithm guided by chaos variables can effectively explore the solution space and accelerate the convergence to the global optimal solution. Its logical mapping formula is shown in equation (71):

[0262] x t+1 = μx t (1 - x t ) (71)

[0263] where x t is the state value at iteration step t, usually in the range of [0, 1], x t+1 is the state value of the next iteration t + 1, and μ is the system parameter, usually in the range of [0, 4], which controls the dynamic behavior of the logical mapping.

[0264] S206: If both the maximum value and the minimum value are within the stable domain of the excavation face, it is determined that the excavation face of the soil bin is stable; if there is a maximum value or a minimum value not within the stable domain of the excavation face, it is determined that the excavation face is unstable.

[0265] The maximum and minimum values of the gradient change of the soil bin pressure plane field are calculated by using the chaotic adaptive particle swarm - sequential quadratic programming algorithm. If both the maximum value and the minimum value are within the stable domain of the excavation face, it is determined that the excavation face of the soil bin is stable; if there is a maximum value or a minimum value not within the stable domain of the excavation face, it is determined that the excavation face is unstable.

[0266] In an exemplary embodiment, the engineering background of the present invention is the civil construction section of a certain rail transit project. The shield tunnel line is arranged in the east-west direction, with a total length of 620 m for the section. The longitudinal section of the line is set with a "V" - shaped slope, and the slopes on both sides are 25.401‰ and 27‰. The line plane is a circular curve with R = 2400 m, and the overburden of the section is about 9.5 - 16.8 m. The shield tunneling method is adopted for the section tunnel. The plane layout of the line between the starting station and the terminal station is as shown in Figure 5 shown. The shield tunnel of the section mainly passes through fine sand and fine-medium sand. The top of the tunnel is fine sand and the bottom is fine-medium sand. The geological profile of the shield tunnel is as shown in Figure 6 shown. According to the geological and hydrogeological conditions of this project, the cross-sectional dimensions of the tunnel, the tunnel burial depth, the line design and other conditions, 1 pressure balance shield machine is invested in this project. The main technical parameters of the balance shield machine are shown in Table 4.

[0267] Table 4

[0268]

[0269] As shown in Figure 3 shown, based on the existing monitoring points E1 - E6 in the soil bin, 12 new soil bin monitoring points E7 - E18 are added by using the cubic spline interpolation method. The soil bin pressure values of E7 - E18 are shown in Table 5.

[0270] After calculation, the mean square error MSE = 0.28 bar. Compared with the actual soil bin pressure numerical range of 1.27 - 3.14 bar for E1 - E6, this error result indicates that the prediction accuracy of the model is within an acceptable range, especially considering the variation range of the soil pressure value. All in all, the above analysis results show that the cubic spline interpolation model exhibits good interpolation and prediction capabilities, and can relatively accurately capture the main change trends of the newly added soil pressure monitoring points E7 - E18 based on the original soil bin pressure data of E1 - E6. Although there are certain prediction errors, the above errors are within an acceptable range.

[0271] Based on the data values in Tables 2 and 3, on the basis of Figure 3 , the distribution of soil bin monitoring points with pressure values is redrawn, and the drawing result is as shown in Figure 7 shown. A B - spline surface is constructed by using the non - uniform B - spline least - squares algorithm shown in formula (58), where the values are as follows: α = β = 3, m = 3, n = 2, r = 18. The input data list of the parameter model is constructed through the coordinate positions and corresponding values of the soil bin pressure monitoring points, and then the fitting is realized through a Python computer program. The input data of the parameter model is shown in Table 4. Based on the data parameters in Table 5, the LSQBivariateSpline function from the scipy.interpolate library is used to fit a B - spline surface defined by control points to the given data points. The fitting result of the pressure gradient field is as shown in Figure 8as shown

[0272] Table 5

[0273]

[0274] Figure 8 shows the contour distribution characteristics of the pressure gradient field in the soil bin of the shield machine under standard working conditions. In the figure, each smooth pressure curve is formed by connecting points with the same soil bin pressure value, and the numbers on the curve mark the corresponding pressure magnitudes, with the unit of bar. It can be intuitively seen from Figure 8 that the pressure field is mainly divided into three major regions, namely the upper, left, and right regions. Among them, the soil bin pressure value in the upper pressure region is significantly less than that in the left and right pressure regions, but the value in the upper pressure region is often the most concerned part by shield construction technicians during the actual construction process and is the main factor affecting the surface settlement above the tunnel. The left and right pressure regions are basically symmetrically distributed, but the area of the left pressure region is slightly larger than that of the right pressure region. The counterclockwise rotation of the cutter head is an important reason for this phenomenon. Analyzing the distribution of the gradient contour line of the entire soil bin pressure field, the pressure gradient is basically symmetric on both sides of the y-axis and can be divided into upper and lower parts along the x-axis direction. The dividing line is approximately at y = -1.0 m. In the upper half of the soil bin, the pressure basically shows a gradually increasing trend and reaches a peak at a position slightly lower than the dividing line in the middle. In addition, in the regions of x ∈ [-2.5, -1.0], y ∈ [-1.5, 0.2] and x ∈ [0.5, 2.5], y ∈ [-1.7, 0.3], two high-pressure regions are formed, and the soil bin pressure values are all above 3 bar. It should be noted that there is a pressure annular region in the upper soil bin pressure region. Based on engineering experience and existing literature, it is speculated that this is because the soil in front of the shield during excavation is not evenly excavated. After the original soil mass in the upper part of the soil bin is stabilized by the cutter, due to the influence of the self-weight of the soil mass, the upper soil mass will fall, resulting in partial collapse. The fallen soil mass accumulates at the lower part of the soil bin. At the same time, due to the influence of the slag discharge of the screw conveyor, the soil bin pressure value at the position near the slag discharge port of the screw conveyor at the lower part of the soil bin decreases instead, and shows a rapid downward trend, thus leading to the above pressure field distribution characteristics. The schematic diagram of the ideal soil bin pressure distribution and the actual soil bin pressure distribution is as Figure 9 shown

[0275] Calculate the stable area of the excavation face under normal working conditions based on formulas (59) - (62). Among them, considering the formation conditions through which the shield passes, the density ρ m of the muck in the soil bin is 1750 kg / m 3 , the cohesion τ a of the soil mass is 2.20 kPa, and the width L of the soil bin is 1 m. Substituting these values into the calculation, the range of the stable area of the excavation face is obtained as

[0276] The chaotic adaptive particle swarm - sequential quadratic programming algorithm is used to evaluate the gradient change. Table 6 shows the algorithm parameter settings. The chaotic adaptive particle swarm - sequential quadratic programming algorithm is implemented by code to determine the stability of the shield excavation face.

[0277] Table 6

[0278]

[0279] Take the case where the soil pressure value above the screw conveyor increases by 30% compared with the normal working condition as an example for the study of the stability determination of the excavation face, that is, the monitoring values of the soil pressure gauges E11 - E14 Figure 7 increase by 30%. The input data of the example parameters are shown in Table 7.

[0280] Table 7

[0281]

[0282] Substitute the data in Table 7 into the code implementation of the stability determination of the excavation face, and the result can be calculated. Compared with the stable area of the normal excavation face of the case project calculated based on the theories of equations (59) - (62), it has exceeded its threshold range, indicating that there is an abnormal situation on the excavation face. The reason is most likely that the slag discharge of the screw conveyor is not timely, resulting in a large soil pressure value in the soil bin above the screw conveyor, and then causing a sharp increase in the local pressure gradient, resulting in an abnormal situation on the excavation face.

[0283] The gradient distribution model of the soil bin pressure field is constructed by the non - uniform cubic higher - order B - spline least - squares method. On this basis, the chaotic adaptive particle swarm - sequential quadratic programming algorithm is introduced to construct a method for determining the stability of the excavation face. The research conclusions are as follows:

[0284] (1) Considering that the number of original soil pressure monitoring points is small and the distribution is relatively scattered, new points are inserted between the existing monitoring points by the cubic spline interpolation method. The leave - one - out cross - validation is used to evaluate the rationality of the numerical values of the newly added soil pressure monitoring points. The calculated mean square error MSE = 0.28 bar. Compared with the actual soil pressure value range of 1.27 - 3.14 bar, the error is well controlled.

[0285] (2) The soil bin pressure field model is constructed by the non - uniform cubic higher - order B - spline least - squares method. The calculated average error is 7.11%. The soil bin pressure gradient field is mainly divided into three major regions. In the upper annular pressure region, the pressure value gradually decreases from the outside to the inside. The left and right pressure regions are basically symmetrically distributed, forming two high - pressure regions with pressure values above 3.0 bar. Due to the counter - clockwise rotation of the cutter head, the area of the left pressure region is slightly larger than that of the right pressure region.

[0286] (3) Simulate the abnormal situation of the soil bin pressure field. Set the soil bin pressure value above the screw conveyor to increase by 30% compared to the normal working condition. Using the chaotic adaptive particle swarm - sequential quadratic programming algorithm, it is calculated that the maximum pressure field gradient change has exceeded the stable domain value of the normal working condition. Therefore, the stability state of the excavation face can be judged by the change of the entire soil bin pressure gradient field.

[0287] When applying the method for judging the stability of the excavation face based on the distribution of the soil bin pressure gradient field provided by the present invention, it is not necessary to execute according to Figure 2 the order of the steps shown. The specific execution order of each step can be determined according to needs, and the present invention does not limit this.

[0288] The above is the method for judging the stability of the excavation face based on the distribution of the soil bin pressure gradient field provided by one or more embodiments of the present invention. Based on the same idea, the present invention also provides a corresponding device for judging the stability of the excavation face based on the distribution of the soil bin pressure gradient field, as Figure 9 shown.

[0289] Figure 10 is a schematic diagram of the device for judging the stability of the excavation face based on the distribution of the soil bin pressure gradient field provided by the present invention, including:

[0290] An initial model construction module 1001, used to construct an initial soil bin pressure plane gradient field distribution model with non - uniform cubic B - spline basis functions; the initial soil bin pressure plane gradient field distribution model includes a function of the soil pressure value in the soil bin changing with position;

[0291] A first determination module 1002, used to take the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function, optimize the control points of the initial soil bin pressure plane gradient field distribution model according to the soil pressure values of the soil bin monitoring points, and determine the optimal control points; the control points are the coefficients of the non - uniform cubic B - spline basis functions;

[0292] A substitution module 1003, used to take the optimal control points as the parameters in the initial soil bin pressure plane gradient field distribution model to obtain the soil bin pressure plane gradient field distribution model;

[0293] A second determination module 1004, used to obtain the stable domain of the excavation face of the soil bin internal pressure gradient field according to the pressure gradient change ranges in the vertical and horizontal directions inside the soil bin;

[0294] A third determination module 1005, used to determine the maximum and minimum values of the soil bin pressure field gradient change according to the soil bin pressure plane gradient field distribution model;

[0295] A determination module 1006 is configured to determine that the excavation face of the soil bin is stable if both the maximum value and the minimum value are within the stable region of the excavation face, and determine that the excavation face is unstable if either the maximum value or the minimum value is not within the stable region of the excavation face.

[0296] For the specific limitations of the device for determining the stability of the excavation face based on the distribution of the soil bin pressure gradient field, reference can be made to the limitations of the method for determining the stability of the excavation face based on the distribution of the soil bin pressure gradient field in the above text, which will not be elaborated here. Each module in the above device for determining the stability of the excavation face based on the distribution of the soil bin pressure gradient field can be implemented in whole or in part by software, hardware, and their combination. The above modules can be embedded in the processor of the computer device in hardware form or be independent of it, or be stored in the memory of the computer device in software form, so as to facilitate the processor to call and execute the operations corresponding to the above modules.

[0297] The present invention also provides a computer-readable storage medium storing a computer program, which can be used to execute the Figure 2 method for determining the stability of the excavation face based on the distribution of the soil bin pressure gradient field provided above.

[0298] The present invention also provides Figure 10 a schematic structural diagram of the computer device as shown in Figure 10 shown. At the hardware level, the computer device includes a processor, an internal bus, a network interface, a memory, and a non-volatile memory. Of course, it may also include other hardware required for other services. The processor reads the corresponding computer program from the non-volatile memory into the memory and then runs it to implement the Figure 2 method for determining the stability of the excavation face based on the distribution of the soil bin pressure gradient field provided above.

[0299] Those of ordinary skill in the art can understand that all or part of the processes in the methods of the above embodiments can be completed by instructing relevant hardware through a computer program. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above methods. Among them, any reference to a memory, storage, database, or other medium used in the various embodiments provided by the present invention can include at least one of non-volatile and volatile memories. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical memory, etc. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0300] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered to be within the scope recorded by the present invention.

Claims

1. A method for determining the stability of an excavation surface based on the distribution of soil bin pressure gradient field, characterized in that: include: The initial soil bin pressure plane gradient field distribution model is constructed using non-uniform cubic B-spline basis functions; The initial soil bin pressure plane gradient field distribution model includes a function of the soil pressure value of the soil bin changing with position; Taking the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function, optimizing the control points of the initial soil bin pressure plane gradient field distribution model according to the soil pressure value of the soil bin monitoring point, and determining the optimal control point; the control point is the coefficient of the non-uniform cubic B-spline basis function; The optimal control point is used as a parameter in the initial soil bin pressure plane gradient field distribution model to obtain the soil bin pressure plane gradient field distribution model; According to the variation range of pressure gradient in the vertical and horizontal directions inside the soil bin, the excavation surface stability domain of the pressure gradient field inside the soil bin is obtained; Determine the maximum and minimum values ​​of the soil bin pressure field gradient change according to the soil bin pressure plane gradient field distribution model; If the maximum value and the minimum value are both within the stable domain of the excavation surface, the excavation surface of the soil bin is determined to be stable; if the maximum value exists or the minimum value is not within the stable domain of the excavation surface, the excavation surface is determined to be unstable.

2. The method according to claim 1, characterized in that The initial soil bin pressure plane gradient field distribution model is: Among them, p(x,y,t) is the non-uniform cubic B-spline basis function, B i,α (x) is the B-spline basis function defined by the knot vector x, B j,β (y) is the B-spline basis function defined by the knot vector y, d ij is the control point, m is the number of vertex limits controlled in the x direction, n is the number of vertex limits controlled in the y direction, i is the random value of the vertex controlled in the x direction, and j is the random value of the vertex controlled in the y direction.

3. The method according to claim 1, characterized in that The soil bin monitoring points include existing soil bin monitoring points and newly added soil bin monitoring points, specifically including: Obtaining the soil pressure value and angle of the existing monitoring point of the soil bin; An angle is selected at equal intervals between any two existing monitoring points of the soil bin as the angle of the newly added monitoring point of the soil bin; Creating a periodic cubic spline interpolation model based on the soil pressure values ​​of the existing monitoring points of the soil bin; According to the periodic cubic spline interpolation model, the soil pressure value corresponding to the angle of the newly added monitoring point of the soil bin is calculated; The existing monitoring points of the soil bin and the newly added monitoring points of the soil bin are used as the soil bin monitoring points.

4. The method according to claim 3, characterized in that The method further comprises: For any newly added monitoring point of the soil bin, the newly added monitoring point of the soil bin is used as a test set, and the remaining monitoring points of the soil bin are used as a training set; the test set includes the soil bin pressure test value of the newly added monitoring point of the soil bin; the training set includes the soil bin pressure test values ​​of the remaining monitoring points of the soil bin; Constructing a cubic spline interpolation model through the training set; According to the cubic spline interpolation model and the test set, a predicted value of the soil bin pressure of the test set is obtained; Calculating the mean square error between the predicted value of the soil bin pressure and the tested value of the soil bin pressure; If the mean square error is within the preset error threshold, it is determined that the setting of the additional pressure monitoring point of the soil bin is reasonable.

5. The method according to claim 1, characterized in that The method uses the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function, optimizes the control points of the initial soil bin pressure plane gradient field distribution model according to the soil pressure values ​​of the soil bin monitoring points, and determines the optimal control points, specifically including: The prediction error of the initial soil bin pressure plane gradient field distribution model is used as the objective function, and the control point is used as a variable. The control point is optimized through the chaotic adaptive particle swarm-sequential quadratic programming algorithm, and the control point corresponding to the minimum objective function value during the iteration process of the chaotic adaptive particle swarm-sequential quadratic programming algorithm is used as the optimal control point.

6. The method according to claim 1, characterized in that The prediction error of the initial soil bin pressure plane gradient field distribution model is used as the objective function, and the control point of the initial soil bin pressure plane gradient field distribution model is optimized according to the soil pressure value of the soil bin monitoring point to determine the optimal control point, which specifically includes: The objective function is: Among them, min x,y (p(x k ,y k ,t)-p k (t)) is the prediction error of the initial soil bin pressure plane gradient field distribution model, p k (t) is the actual value of soil bin pressure, and p(x,y,t) is the predicted value of soil bin pressure; The optimal solution of the objective function is converted into the optimal control point, and the first function is defined as: Among them, B i,α (x) is the B-spline basis function defined by the knot vector x, B j,β (y) is the B-spline basis function defined by the knot vector y, d ij is the control point, p k (t) is the true value of soil bin pressure, x k is the random optimization value in the x direction, i is the random value in the x direction, j is the random value in the y direction, m is the number of control vertex boundaries in the x direction, n is the number of control vertex boundaries in the y direction, k is the number of control vertex optimization times, and r is the total number of control vertex optimization times. When the first function is minimum, we can get Will Substituting the first function into the equation, we get the first equation: Simplifying the first equation, the first matrix expression is: B T BD(t)=B T P(t); Among them, B T is the transpose of the B-spline basis function matrix, B is the B-spline basis function matrix, D(t) is the value of the control vertex at time t, and P(t) is the value of the soil bin pressure at time t; The singular value decomposition method is used to solve the first matrix expression, and the control points are obtained as follows: Among them, D(t) is the value of the control vertex at time t, U is the orthogonal matrix, V T is the orthogonal matrix transpose, P(t) is the soil bin pressure value at time t, are the matrix singular values. The singular values ​​of the matrix B are calculated, and the singular values ​​are substituted into the control points to obtain the optimal control points.

7. The method according to claim 1, characterized in that The excavation surface stability domain of the pressure gradient field inside the soil bin is obtained according to the pressure gradient variation range in the vertical direction and the horizontal direction inside the soil bin, specifically including: Determine the minimum value of the pressure change inside the soil bin according to the minimum value of the pressure gradient change range in the vertical direction inside the soil bin and the minimum value of the pressure gradient change range in the horizontal direction inside the soil bin; The maximum value of the pressure change inside the soil bin is determined according to the maximum value of the pressure gradient change range in the vertical direction inside the soil bin and the maximum value of the pressure gradient change range in the horizontal direction inside the soil bin; According to the minimum and maximum values ​​of the pressure change inside the soil bin, the excavation surface stability domain of the pressure gradient field inside the soil bin is determined.

8. The method according to claim 7, characterized in that The pressure gradient variation range in the vertical direction inside the soil bin is: in, is the range of vertical pressure gradient inside the soil bin, ρ m is the density of the soil inside the soil bin, τ a is the cohesion of the soil inside the soil bin, L is the vertical distance of the soil bin, and g is the acceleration of gravity; The pressure gradient variation range in the horizontal direction inside the soil bin is: in, is the range of pressure gradient variation in the horizontal direction inside the soil bin; The calculation formula for the minimum and maximum values ​​of the internal pressure change of the soil bin is: Among them, ‖▽P‖ is the value of the pressure change inside the soil bin, is the change value of the vertical pressure gradient inside the soil bin, It is the change value of the horizontal pressure gradient inside the soil bin.

9. The method according to claim 1, characterized in that The soil bin pressure plane gradient field distribution model constructs multiple B-spline curves based on the same non-uniform cubic B-spline basis function, and the soil bin pressure value on each B-spline curve is the same.

10. An excavation surface stability determination device based on soil bin pressure gradient field distribution, characterized in that: include: The initial model building module is used to build the initial soil bin pressure plane gradient field distribution model using non-uniform cubic B-spline basis functions; The initial soil bin pressure plane gradient field distribution model includes a function of the soil pressure value of the soil bin changing with position; The first determination module is used to use the prediction error of the initial soil bin pressure plane gradient field distribution model as the objective function, optimize the control points of the initial soil bin pressure plane gradient field distribution model according to the soil pressure value of the soil bin monitoring point, and determine the optimal control point; the control point is the coefficient of the non-uniform cubic B-spline basis function; A substitution module is used to use the optimal control point as a parameter in the initial soil bin pressure plane gradient field distribution model to obtain the soil bin pressure plane gradient field distribution model; The second determination module is used to obtain the excavation surface stability domain of the pressure gradient field inside the soil bin according to the pressure gradient variation range in the vertical direction and the horizontal direction inside the soil bin; A third determination module is used to determine the maximum and minimum values ​​of the soil bin pressure field gradient change according to the soil bin pressure plane gradient field distribution model; A judgment module is used to determine that the excavation surface of the soil bin is stable if the maximum value and the minimum value are both within the stable domain of the excavation surface, and to determine that the excavation surface is unstable if the maximum value exists or the minimum value is not within the stable domain of the excavation surface.

Images