Shield tunneling machine deviation rectification track high-precision autonomous planning method and system

Abstract

The invention belongs to the field of shield tunneling machine deviation rectification track planning, and discloses a shield tunneling machine deviation rectification track high-precision autonomous planning method and system, and the method comprises the steps: determining the minimum turning radius of a shield tunneling machine; carrying out coordinate calculation on the DTA to obtain the coordinate, the direction and the curvature of each point; projecting the starting point state and the end point state of the shield tunneling machine under the Cartesian coordinate system into the established Frenet coordinate system, establishing a correction curve parameterization equation based on a quintic polynomial line type, and converting the state of the shield tunneling machine into the Cartesian coordinate system; through segment pre-typesetting constraint, determining a deviation correction curve end point feasible region; the deviation degree of the deviation correction curve and the DTA and the length and smoothness of the deviation correction curve serve as optimization parameters, normalization processing is carried out according to the threshold value of each target, a comprehensive cost function is established, an improved particle swarm optimization algorithm is used for solving, and a global optimal deviation correction track is output. According to the method, real-time autonomous high-precision deviation rectification track planning of the shield tunneling machine can be realized, the track feasibility is improved, the pose deviation degree is reduced, and the tunnel quality and the construction efficiency are improved.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of trajectory deviation correction planning of a shield tunneling machine, and in particular relates to a high-precision autonomous trajectory deviation correction planning method and system for a shield tunneling machine. BACKGROUND

[0002] A shield tunneling machine has become the core equipment for tunnel construction of subways, highways, water conservancy projects and energy pipelines due to its high efficiency, safety, low cost and small environmental impact. In an ideal construction state, the tunneling trajectory of the shield should be accurately coincident with the designed tunnel axis to ensure the quality of the formed tunnel and the efficiency of the construction. The "Code for Construction and Acceptance of Shield Tunneling" clearly stipulates that the deviation between the shield and the DTA should not exceed 50 mm during the shield construction process. Due to factors such as complex geological conditions, equipment response lag and human operation errors in actual construction, the shield machine often exhibits pose deviation phenomena such as head knocking, upward movement, horizontal deviation and snake-shaped trajectory, which may even lead to tunnel deviation exceeding the limit and even require re-excavation. To achieve high-precision pose deviation control of the shield machine, the deviation between the measured pose of the shield machine and the DTA needs to be planned to correct the trajectory, and the shield machine needs to be controlled to follow the corrected trajectory to quickly and smoothly return to the DTA. Although the pose measurement and propulsion control technologies are relatively mature, there is still a lack of systematic research on high-precision trajectory correction planning. In the face of the trend of shield intelligence and unmannedization, high-precision autonomous trajectory correction planning technology for the shield machine has become the basis for autonomous excavation of the shield.

[0003] Although there have been some related researches on trajectory correction planning of the shield machine, there are still significant technical defects: first, the calculation of the minimum turning radius of the shield machine lacks systematic consideration, affecting the actual feasibility; second, the widely used cubic polynomial or circular curve line type cannot guarantee the continuity of the curvature at the start and end points, resulting in large vibrations during the propulsion process of the shield machine and affecting the quality of the tunnel construction; in addition, since the shield machine needs to assemble segments after completing excavation to form the final tunnel lining, the existing researches do not consider this particularity of the trajectory planning of the shield machine, resulting in the inability to guarantee the smooth assembly of the segments or even causing damage to the segments.

[0004] In view of the above analysis, the existing technology has the following technical problems that need to be solved urgently:

[0005] The existing trajectory correction planning of the shield machine has the problems of poor trajectory feasibility, poor construction quality and the like caused by imperfect calculation of the minimum curvature radius, discontinuity of the curvature of the corrected trajectory and lack of constraints for segment assembly, which cannot meet the precision and efficiency requirements of autonomous excavation of the shield, and therefore there is an urgent need to construct a high-precision autonomous trajectory correction planning method for the shield machine that integrates accurate calculation of the minimum turning radius, continuous curvature of the corrected trajectory and satisfies the constraints of segment pre-arrangement. SUMMARY

[0006] In view of the problems of the prior art, the present application provides a high-precision autonomous planning method and system for a deviation trajectory of a shield tunneling machine.

[0007] The present application is implemented in a high-precision autonomous planning method for a deviation trajectory of a shield tunneling machine, characterized in that the high-precision autonomous planning method for a deviation trajectory of a shield tunneling machine specifically comprises the following steps:

[0008] S1: determining the minimum turning radius of the shield tunneling machine through geometric constraints and mechanical constraints;

[0009] S2: performing coordinate calculation on the DTA according to the start point and end point coordinates, the mileage and the characteristics of the curve of each curve segment on the DTA, to obtain the coordinates, direction and curvature of each point on the DTA;

[0010] S3: establishing a Frenet coordinate system, projecting the start point and end point states of the shield tunneling machine in the Cartesian coordinate system into the Frenet coordinate system, and establishing a parameterized equation of the deviation curve corresponding to the end point of the deviation curve in the Frenet coordinate system based on a quintic polynomial linearization, and converting the state of the shield tunneling machine in the Frenet coordinate system into the Cartesian coordinate system;

[0011] S4: determining the feasible region of the end point of the deviation curve through segment pre-layout constraints;

[0012] S5: taking the deviation degree of the deviation curve from the DTA, the length and smoothness of the deviation curve as optimization parameters, performing normalization processing according to the threshold values of the targets, and establishing a comprehensive cost function in a weighted summation manner;

[0013] S6: solving the comprehensive cost function by using an improved particle swarm optimization algorithm, and outputting a globally optimal deviation trajectory.

[0014] Further, the step S1 comprises the following steps of determining the minimum deviation curvature radius of the shield tunneling machine through geometric constraints and mechanical constraints:

[0015] Step S11: calculating the minimum deviation curvature radius under the shield tail gap constraint:

[0016]

[0017] where D g is the segment diameter, T n is the shield tail gap in the turning direction of the shield tunneling machine, l g is the length of the segment inside the shield tail, which is taken as 1.5B g according to engineering experience, B g is the segment width, and D g is the segment diameter.

[0018] Calculate the minimum turning radius under the maximum thrust cylinder difference constraint:

[0019]

[0020] Where D F is the thrust cylinder installation diameter, l s is the distance of the shield machine advancing along the correction curve when turning, generally taken as B g , U nmax is the maximum thrust cylinder stroke difference allowed in construction.

[0021] Calculate the minimum turning radius under the shield machine hinge device constraint:

[0022]

[0023] Where, is the limit hinge angle of the shield machine, D s is the front shield cutterhead diameter, L1 and L2 are the lengths of the front shield and rear shield of the shield machine respectively.

[0024] The minimum turning radius under geometric constraint is:

[0025] R minG = max(R mT ,R mU ,R mc )

[0026] Step S12, according to the front soil resistance f1, the friction force f2 between the shield shell and the surrounding soil, the turning soil resistance f3 and the thrust force F t , a statics model of the turning process of the shield machine is established, and the minimum turning radius R tmax under the mechanical constraint is obtained by the maximum thrust force F minF ;

[0027] Step S13, introduce the safety factor, and the corrected minimum turning radius of the shield machine is:

[0028] R min =1.2max(R minG ,R minF )

[0029] Furthermore, in step S3, a Frenet coordinate system is established, and the starting and ending states (x, y, θ, κ) of the shield machine in the Cartesian coordinate system are projected into the Frenet coordinate system (s, l, l', l'), where x and y are the coordinates of the cutterhead center of the shield machine in the Cartesian coordinate system, θ is the heading of the shield machine, κ is the curvature of the current tunneling trajectory of the shield machine, s is the tunneling mileage of the cutterhead center of the shield machine in the Frenet coordinate system, l is the deviation between the cutterhead center of the shield machine and the DTA, l' and l' are the first and second derivatives of l with respect to s, respectively. In the Frenet coordinate system, based on the quintic polynomial line type, a parameterized equation of the correction curve corresponding to the end point of the correction curve is established. The specific steps are as follows:

[0030] Step S31, set the boundary conditions: the deviation at the starting mileage s0 is l0, the slope is l'0, the curvature is l'0', and the end mileage s i The deviation, slope and curvature are all 0.

[0031] Step S32: Substitute the boundary conditions into the polynomial equation and describe it in the form of a matrix as follows:

[0032]

[0033] in, is the coefficient vector of the polynomial, is the boundary condition vector,

[0034] Step S33: Solve the polynomial coefficients using the following equation to generate a quintic polynomial curve equation uniquely determined by the end point of the correction curve.

[0035]

[0036] Furthermore, in step S4, the feasible region of the end point of the correction curve is determined by the pre-layout constraints of the segments. The specific steps are as follows:

[0037] Step S41, initializing the range s of the correction end point i ∈[s min ,s max ], where s min Take s0, s max Determined based on engineering experience, generally the larger value is taken.

[0038] Step S42, calculate the median value Solve the end point of the correction curve as s mid The corresponding correction curve equation is obtained, and the feasibility of the correction curve is verified through segment pre-layout;

[0039] Step S43: Update the search interval boundary value. If the segment pre-layout verification is passed, then smax = s mid ; if the pre-assembly of the segment cannot be verified, s min = s mid .

[0040] Steps S42 and S43 are continued until the search interval width is less than the preset precision threshold ∈, and the minimum value of the output end point of the correction curve is The maximum value is the initial s max .

[0041] Further, the feasibility of the correction curve is verified by the pre-assembly of the segment in step S42, specifically:

[0042] (1) The relative pose (x s , y s , z s , α s , β s , γ s ) of the last ring of assembled segments and the shield machine is calculated by the pose (x u,n , y d,n , z l,n , α r,n , β u,n , γ d,n ) of the shield machine, the front, rear, left and right tail gaps and the push cylinder strokes (T l,n , T r,n , T g,n , T g,n ) and (U g,n , U g,n , U g,n , U g,n ):

[0043]

[0044] Wherein, L F is the distance from the installation plane of the push cylinder to the cutter head, c y is the length when the push cylinder has zero stroke, i n is the point position of the last ring of assembled segments, and the number of point positions of the segment is N g .

[0045] (2) The homogeneous transformation matrix of the last ring of assembled segments relative to the shield machine and the homogeneous transformation matrix of the shield machine relative to the world coordinate system are calculated to obtain the homogeneous transformation matrix of the last ring of assembled segments relative to the world coordinate system

[0046]

[0047] (3) The pose transfer matrix between two adjacent ring segments related to the segment points Obtain the pose transfer matrix of the next ring of segments to be assembled

[0048]

[0049] (4) By matrix The inverse solution is used to obtain the position of the next ring of segments to be assembled (x n+1 ,y n+1 ,z n+1 ,α n+1 ,β n+1 ,γ n+1 ).

[0050] (5) Calculate the deviation between the next ring of segments to be assembled and the DTA at different assembly points, and select the assembly point that minimizes the deviation to pre-assemble the next ring of segments.

[0051] Repeat steps (3) to (5) until the segment pre-layout for the entire correction curve is completed.

[0052] In the entire pre-layout plan, if the maximum axis fitting deviation between the segment and DTA is less than 20mm, the correction curve passes the segment pre-layout verification, otherwise it fails.

[0053] Furthermore, in step S5, the deviation degree between the correction curve and the DTA, the length and smoothness of the correction curve are used as optimization parameters, normalization is performed according to the physical threshold of each target, and a comprehensive cost function is established in a weighted summation manner, specifically including:

[0054] (1) Deviation cost:

[0055]

[0056] Among them, max|l(s)| is the maximum deviation between the correction curve and DTA, l max is the maximum allowable deviation between the shield machine and DTA, s0 is the starting point of the correction curve, and s i The end point of the correction curve.

[0057] (2) Length cost:

[0058]

[0059] L max The operator can set the length threshold of the correction curve according to the tendency.

[0060] (3) Smoothness cost:

[0061]

[0062] r min (s i ) is the minimum curvature radius of the correction curve, R min is the minimum turning radius of the shield machine.

[0063] (4) The penalty term:

[0064] P s = μ s I(r min (s i )< R min )

[0065] P d = μ d I(max |l(s)|>l max )

[0066] Wherein, P s and P d are the smoothness penalty term and the deviation penalty term, respectively, μ s = 200, μ d = 10 are the smoothness penalty factor and the deviation penalty factor, respectively, and I(·) is an indicative function, which is 1 when the condition is met, and 0 when the condition is not met.

[0067] (5) The comprehensive objective function:

[0068] C total (s) = λ d C d + λ l C l + λ s C s + P s + P d

[0069] Wherein, λ d , λ l and λ s are the weight coefficients of the deviation cost, the length cost and the smoothness cost, respectively, all belonging to [0, 1], and satisfying λ d + λ l + λ s = 1.

[0070] Further, the step S6 uses the improved particle swarm optimization algorithm to solve the comprehensive cost function, and outputs the globally optimal correction trajectory, and the specific steps are as follows:

[0071] Step S61, parameter setting: including particle swarm size N, maximum iteration number G m , acceleration coefficients c1 and c2, extreme values of inertia weight ω min and ω max , extreme values of particle position and s max , the maximum value of particle velocity V max , the fractal parameter μ and the disturbance intensity coefficient σ.

[0072] Step S62, population initialization: the population is initialized using Tent chaotic mapping, the dimension of the individual is 1, the particle position represents the mileage of the end point of the correction curve, and the speed of the individual is where 0 < k ≤ N, G is the current iteration number of the algorithm, 0 ≤ G ≤ G m .

[0073] Step S63, fitness value calculation: the corresponding correction curve equation of each individual is calculated, which is brought into the comprehensive objective function, and the fitness value f k of each particle is calculated.

[0074] Step S64, individual optimal solution update: for each particle, the fitness value is compared with the fitness value of the individual historical optimal position p kb , if better, update p kb , otherwise maintain the individual optimal position unchanged.

[0075] Step S65, global optimal solution update: for each particle, the fitness value is compared with the fitness value of the global historical optimal position g b , if better, update g b , otherwise maintain the global optimal position unchanged.

[0076] Step S66, update the particle position and speed of the next generation, the update formula is as follows:

[0077]

[0078] where c1 and c2 are acceleration coefficients, respectively controlling cognitive component and social component, r1 and r2 are uniform random numbers in the range of [0, 1]. ω is the inertia weight, N(0, 1) represents a standard Gaussian distribution random variable, and σ is the disturbance intensity coefficient. The inertia weight ω is adaptively adjusted, specifically:

[0079] The inertia weight is gradually reduced with iteration, as shown in the following formula:

[0080]

[0081] Step S67, boundary condition processing: when the position or speed of the particle exceeds the set value, the boundary condition processing strategy is used to limit the particle in the feasible search space, as shown in the following formula:

[0082]

[0083] When the speed exceeds the boundary, set it to the maximum value V max When the position exceeds the boundary, reverse the speed and update the position.

[0084] Repeat steps S63-S67 until the maximum number of iterations is reached, and output the global optimal solution g b , and obtain the optimal curve correction point and the correction curve equation.

[0085] Further, the step S62 uses Tent chaotic mapping to initialize the population, specifically:

[0086] Generate an initial chaotic variable a1 randomly in the interval [0, 1], and generate a chaotic sequence {a k}

[0087]

[0088] Based on the obtained chaotic sequence, map the chaotic variable to the search space of the optimization problem through linear transformation to obtain the initial population position of the particle swarm, as shown in the following formula:

[0089]

[0090] In the formula, represents the initial position of the kth particle.

[0091] Another object of the present application is to provide a high-precision autonomous planning system for a shield tunneling machine, which specifically comprises:

[0092] A turning radius determination module is configured to determine the minimum turning radius of the shield tunneling machine through geometric constraints and mechanical constraints.

[0093] A coordinate calculation module is configured to perform coordinate calculation on the DTA to obtain the coordinates, direction and curvature of each point on the DTA.

[0094] A coordinate system conversion module is configured to convert the state of the shield tunneling machine in the Frenet coordinate system and the Cartesian coordinate system.

[0095] A correction curve endpoint feasible region determination module is configured to determine the correction curve endpoint feasible region through segment pre-layout constraints.

[0096] A comprehensive cost function solving module is configured to establish a comprehensive cost function in a weighted summation manner, and solve the comprehensive cost function using an improved particle swarm optimization algorithm to output a globally optimal correction trajectory.

[0097] In combination with the above technical solutions and the technical problems solved, the technical solution to be protected by the present application has the following advantages and positive effects:

[0098] Firstly, the present application accurately calculates the minimum turning radius of the shield tunneling machine by fusing geometric constraints (shield tail gap, maximum cylinder stroke difference, hinged device) and mechanical constraints, improves the accuracy of deviation correction capability evaluation; adopts a quintic polynomial to construct a parameterized deviation correction curve, ensures the continuity and derivability of the curvature, and prevents vibration of the machine body caused by sudden pressure changes in the propulsion system; establishes segment pre-layout constraints based on the real-time state of the shield tunneling machine, limits the feasible region of the deviation correction curve endpoint, and ensures the compatibility of the planned trajectory and the lining structure; and adopts an improved particle swarm optimization algorithm to optimize the deviation correction curve, ensuring the global optimality of the deviation correction curve.

[0099] Through the present application, real-time autonomous high-precision deviation correction trajectory planning of the shield tunneling machine can be achieved, trajectory feasibility can be significantly improved, the degree of pose deviation can be reduced, and the tunnel quality and construction efficiency can be improved.

[0100] Secondly, the high-precision autonomous planning method for the deviation correction trajectory of the shield tunneling machine provided by the present application can significantly suppress pose abnormal phenomena such as shield tunneling machine tapping, upward movement and snake-shaped trajectory by generating a smooth and continuous deviation correction curve, reduce construction rework rate and equipment vibration, and effectively prolong the service life of key components. At the same time, as the basis for intelligent and unmanned tunneling of the shield tunneling machine, the technology can improve construction efficiency, reduce manual intervention, and create significant engineering economic benefits.

[0101] There is a lack of systematic research in the field of deviation correction trajectory planning of the shield tunneling machine at present: current researches generally focus on optimization of the propulsion control system, and simplified models are mostly used for trajectory generation algorithms. The present application fuses geometric constraints such as shield tail gap, cylinder stroke difference and hinged device with mechanical constraints such as soil resistance and friction, establishes an accurate calculation model for the minimum turning radius of the shield tunneling machine, ensures the continuity of the curvature of the deviation correction curve everywhere through a quintic polynomial parameterization equation, eliminates sudden pressure changes in the propulsion system, and corrects the deviation curve based on the real-time shield tunneling machine pose, shield tail gap and segment pre-layout constraints of the propulsion cylinder stroke, narrows the feasible region of the deviation correction curve endpoint by bisection method, and searches for the global optimal solution of the comprehensive cost function in the feasible region by using an improved particle swarm optimization algorithm.

[0102] The present application first realizes comprehensive autonomous planning of "minimum turning radius calculation, high-order curve generation, segment pre-layout restriction, and global optimization", filling the technical gap in high-precision systematic trajectory generation for intelligent tunneling of the shield tunneling machine.

[0103] Existing methods regard trajectory planning and segment assembly as independent links, the present application introduces segment pre-layout based on real-time shield tunneling machine pose, shield tail gap and cylinder stroke to restrict the feasible region of the deviation correction curve endpoint, so that the deviation correction endpoint meets the needs of lining assembly, ensures the feasibility of the deviation correction curve, and avoids accidents such as segment damage. BRIEF DESCRIPTION OF DRAWINGS

[0104] Figure 1 is a flowchart of the high-precision autonomous planning method for the deviation correction trajectory of the shield tunneling machine provided by the present application.

[0105] Figure 2 is a feasible region flowchart for determining a correction curve endpoint provided by an embodiment of the present application based on segment pre-layout constraints;

[0106] Figure 3 is a schematic diagram for calculating the relative pose of the last ring assembled segment and the shield machine provided by an embodiment of the present application;

[0107] Figure 4 is a shield machine correction trajectory algorithm flowchart based on an improved particle swarm optimization provided by an embodiment of the present application;

[0108] Figure 5 is a shield machine correction trajectory high-precision autonomous planning system module diagram provided by an embodiment of the present application. DETAILED DESCRIPTION

[0109] In order to make the purpose, technical solutions and advantages of the present application clearer and more apparent, the present application will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and do not limit the present application.

[0110] When the shield machine advances in a complex stratum, it is often subject to the constraints of the excavation torque and the stratum reaction force, resulting in difficulty in accurately controlling the deviation trajectory. The traditional linear interpolation or simple curve fitting method cannot take into account both the smoothness and the feasibility. The present method first analyzes the minimum turning radius of the shield machine under the synergistic action of geometric and mechanical constraints, so as to determine the limit motion capability of the advancing device and ensure that the trajectory planning meets the requirements of the section assembly pre-layout while not exceeding the safety red line of the mechanical performance.

[0111] According to the characteristics of each curve segment of the DTA (design tunnel axis), the present method uses a refined coordinate solving process to obtain the spatial coordinates, direction and curvature of the continuous point cloud through joint analysis of the mileage and the curvature. This process lays a foundation for high-precision data support for subsequent parameterized curve construction.

[0112] In the Frenet coordinate system, the correction segment is parameterized by using a quintic polynomial curve, so that the curve is second-order smooth in terms of the curvature and the tangential continuity at the start and end. Through bidirectional mapping between coordinate systems, the shield machine reference frame and the global Cartesian frame can be dynamically switched, realizing the decoupling of curve generation and mechanical motion, while taking into account the dual demands of construction reference and measurement feedback.

[0113] In combination with the pre-set segment assembly constraints, the correction endpoint is limited within the feasible region. This feasible region not only covers the geometric position, but also introduces engineering variables such as the shield tail gap, the advancing cylinder stroke and the segment misjoint assembly, so as to ensure that the trajectory landing point has construction implementability and reduces the risk of subsequent lining damage, etc.

[0114] To balance the deviation, length and overall smoothness of the trajectory, a multi-objective comprehensive cost function is constructed, and the weight distribution mechanism of each index is refined, so that the optimization target is highly consistent with the actual construction demand. The deviation term in the cost function measures the overall error, the length term is used to control the construction period and the advancing efficiency, and the smoothness term ensures the continuous stability of the shield machine movement and prevents the mutation from intensifying the stratum disturbance.

[0115] On this basis, an improved particle swarm optimization algorithm is used to perform global search on the cost function. The population is initialized by Tent chaotic mapping, the inertia weight is dynamically adjusted, and Gaussian disturbance is added to the speed update, so as to overcome the defects of the traditional PSO that it is easy to fall into local extremum and the optimization time is relatively long. The step size is adaptively adjusted in the iteration process, and the inertia factor is corrected in real time in combination with the feedback of the construction sensor, and finally a high-precision correction trajectory is output, which takes into account the engineering feasibility and the mechanical stability.

[0116] As shown in Figure 1 , the embodiment of the present application provides a high-precision autonomous planning method for a shield machine correction trajectory, which specifically comprises:

[0117] S1: Determine the minimum turning radius of the shield machine through geometric constraints and mechanical constraints;

[0118] S2: According to the start and end point coordinates, mileage and characteristics of the curve of each curve segment on the DTA, perform coordinate calculation on the DTA to obtain the coordinates, direction and curvature of each point on the DTA;

[0119] S3: Establish a Frenet coordinate system, and project the start and end point states (x, y, θ, κ) of the shield machine in the Cartesian coordinate system into the Frenet coordinate system (s, l, l', l''). Wherein, x and y are the coordinates of the center of the cutter head of the shield machine in the Cartesian coordinate system, θ is the heading of the shield machine, κ is the curvature of the current tunneling trajectory of the shield machine, s is the tunneling mileage of the center of the cutter head of the shield machine in the Frenet coordinate system, l is the deviation of the center of the cutter head of the shield machine from the DTA, l' and l'' are the first and second derivatives of l with respect to s respectively. The conversion equation is:

[0120] s=s0

[0121]

[0122] l'=(1-lκ m )tan(θ-θ m )

[0123]

[0124] Wherein, s0 is the current position of the shield tunneling machine corresponding to the DTA mileage, sign(t) is a sign function, when t>0, sign(t)=1; when t<0, sign(t)=-1.

[0125] In the Frenet coordinate system, based on the quintic polynomial linearity, the parameterization equation of the correction curve corresponding to the end point of the correction curve is established.

[0126] The state of the shield tunneling machine in the Frenet coordinate system is converted into the Cartesian coordinate system. The conversion equation is:

[0127] x=x m -lsinθ m

[0128] y=y m +lcosθ m

[0129]

[0130] S4: Determine the feasible region of the end point of the correction curve by the segment pre-layout constraint;

[0131] S5: Take the deviation degree of the correction curve from the DTA, the length and smoothness of the correction curve as the optimization parameters, normalize according to the threshold value of each target, and establish a comprehensive cost function in a weighted summation manner.

[0132] S6: Solve the comprehensive cost function by using the improved particle swarm optimization algorithm, and output the global optimal correction trajectory.

[0133] The step S1 geometric constraint and mechanical constraint determine the minimum correction curvature radius of the shield tunneling machine, comprising the following steps:

[0134] Step S11, calculate the minimum correction curvature radius under the shield tail gap constraint:

[0135]

[0136] Wherein D g is the segment diameter, T n is the shield tail gap in the turning direction of the shield tunneling machine, and l is the length of the segment inside the shield tail, which is taken as 1.5B g according to engineering experience, B g is the segment width, and D g is the segment diameter.

[0137] Calculate the minimum turning radius under the maximum push cylinder difference constraint:

[0138]

[0139] Wherein, DF To promote the cylinder installation diameter, l s For the shield machine rotates along the deviation curve to advance the distance, generally taken as B g , U nmax The maximum allowable push cylinder stroke difference in construction.

[0140] The minimum turning radius of the shield machine under the constraint of the articulated device is calculated:

[0141]

[0142] Wherein, The limit articulation angle of the shield machine, D s The diameter of the front shield cutter head, L1 and L2 are the lengths of the front shield and rear shield of the shield machine respectively.

[0143] The minimum turning radius under geometric constraints is:

[0144] R minG =max(R mT ,R mU ,R mc )

[0145] Step S12, according to the front soil resistance f1, the friction force f2 between the shield shell and the surrounding soil, the turning soil resistance f3 and the pushing force F t The statics model of the turning process of the shield machine is established, and the minimum turning radius R tmax Under the mechanical constraint, the maximum pushing force F minF ;

[0146] Step S13, introduce the safety factor, the corrected minimum turning radius of the shield machine is:

[0147] R min =1.2max(R minG ,R minF )

[0148] The step S3 is based on a quintic polynomial linearization in the Frenet coordinate system, and the parameterization equation of the deviation curve corresponding to the end point of the deviation curve is established, and the specific steps are:

[0149] Step S31, the purpose of deviation trajectory planning is to make the shield machine return smoothly from the current position to DTA, so it is necessary to not only ensure the continuity of the shield machine coordinates (l), but also ensure the continuity of the slope (l') and the curvature (l'') of the path, and the boundary conditions are set as: the lateral deviation is l0, the slope is l'0, and the curvature is l'0' at the starting point mileage s0, and the lateral deviation, the slope and the curvature are all 0 at the end point mileage s i .

[0150] Step S32, substitute the boundary condition into the polynomial equation, described in the form of matrix as:

[0151]

[0152] wherein, is the coefficient vector of the polynomial, is the boundary condition vector,

[0153] Step S33, solve the polynomial coefficient by the following formula, to generate the quintic polynomial curve equation uniquely determined by the end point of the correction curve.

[0154]

[0155] The step S4 determines the feasible region of the end point of the correction curve by the segment pre-layout constraint, as shown in Figure 2 The specific steps are as follows:

[0156] Step S41, initialize the range s i ∈[s min ,s max ], wherein s min takes s0, s0 is the mileage of the DTA corresponding to the current position of the shield machine, and s max is determined according to the engineering experience value, generally taking a larger value, and in this embodiment, s0+100m.

[0157] Step S42, calculate the intermediate value Solve the correction curve equation corresponding to s mid , and verify the feasibility of the correction curve by the segment pre-layout;

[0158] Step S43, update the search interval boundary value, if verified by the segment pre-layout, s max =s mid ; if not verified by the segment pre-layout, s min =s mid .

[0159] Continue steps S42 and S43 until the search interval width is less than the preset precision threshold ∈, output the minimum value of the end point of the correction curve as , and the maximum value as the initial s max .

[0160] The step S42 verifies the feasibility of the correction curve by the segment pre-layout, specifically as follows:

[0161] (1) Verify the feasibility of the correction curve by the segment pre-layout according to the shield machine pose (x s , y s , z s , αs ,β s ,γ s ), the shield tail clearance and the thrust cylinder stroke (T u,n , T d,n , T l,n , T r,n ) and (U u,n , U d,n , U l,n , U r,n ) on the current state, the relative pose (x g,n , y g,n , z g,n , a g,n , b g,n , g g,n ) of the last ring of assembled segments and the shield machine is calculated, as shown in formula (1): Figure 3

[0162]

[0163] wherein, L F is the distance from the installation plane of the thrust cylinder to the cutter head, c y is the length of the thrust cylinder when the zero stroke, i n is the point of the last ring of assembled segments, and the point number of the segment is N g .

[0164] (2) The homogeneous transformation matrix of the last ring of assembled segments relative to the shield machine is calculated by the homogeneous transformation matrix of the last ring of assembled segments relative to the shield machine and the homogeneous transformation matrix of the shield machine relative to the world coordinate system .

[0165]

[0166] wherein, the homogeneous transformation matrix of the last ring of assembled segments relative to the shield machine is:

[0167]

[0168] wherein, the sin and cos functions are represented by c and s.

[0169] The homogeneous transformation matrix of the shield machine relative to the world coordinate system is:

[0170]

[0171] (3) The pose transmission matrix of the next ring of segments to be assembled is obtained by the pose transmission matrix of the adjacent two rings of segments related to the segment point ​​

[0172]

[0173] (4) by the matrix The inverse solution is to get the pose (x n+1 ,y n+1 ,z n+1 ,α n+1 ,β n+1 ,γ n+1 ) of the next ring to be assembled pipe piece.

[0174] (5) Calculate the deviation of the next ring to be assembled pipe piece and DTA at different assembly points, and select the assembly point that makes the deviation minimum to pre-assemble the next ring pipe piece.

[0175] Repeat steps (3) to (5) until the whole deviation correction curve is completed.

[0176] In the whole pre-layout scheme, the maximum axis fitting deviation of the pipe piece and the DTA is less than 20mm, then the deviation correction curve passes the pipe pre-layout verification, otherwise it does not pass.

[0177] The step S5 takes the deviation of the deviation correction curve from the DTA, the length and smoothness of the deviation correction curve as optimization parameters, normalizes according to the physical threshold of each target, and establishes a comprehensive cost function in a weighted summation manner, specifically including:

[0178] (1) Deviation cost:

[0179]

[0180] Where, max | l (s) | is the maximum deviation of the deviation correction curve from the DTA, l max is the maximum allowable deviation of the shield machine from the DTA, which is taken as 50mm, s0 is the starting point of the deviation correction curve, and s i is the end point of the deviation correction curve.

[0181] (2) Length cost:

[0182]

[0183] L max is the length threshold of the deviation correction curve set by the operator according to the inclination, which can be taken as 50m in general case.

[0184] (3) Smoothness cost:

[0185]

[0186] r min (s i ) is the minimum curvature radius of the deviation correction curve, R minThe minimum turning radius of the shield machine.

[0187] (4) The penalty term:

[0188] The fault tolerance of various constraints in shield construction is different. The constraint conditions are classified according to the priority: first, the curvature radius of the correction curve must be strictly greater than the minimum turning radius of the shield machine, which is the highest priority constraint to ensure the feasibility of the tunneling path; second, the deviation of the shield machine tunneling axis from the DTA needs to be controlled within the allowable range, which is the sub-optimal priority constraint to ensure the engineering accuracy; finally, the length constraint of the correction curve is the lowest priority constraint, allowing for moderate compromise when necessary. Based on this priority system, a stepwise penalty factor mechanism is introduced to strengthen the rigidity of high-priority constraints through the penalty term, which is expressed as follows:

[0189] P s = μ s I(r min (s i )< R min )

[0190] P d = μ d I(max | l(s) | > l max )

[0191] Where P s and P d are the smoothness penalty term and the deviation penalty term, μ s = 200, μ d = 10 are the smoothness penalty factor and the deviation penalty factor, and I(·) is an indicator function that is 1 when the condition is met and 0 when the condition is not met.

[0192] (5) The comprehensive objective function:

[0193] C total (s) = λ d C d + λ l C l + λ s C s + P s + P d

[0194] Where λ d , λ l and λ s are the weight coefficients of the deviation cost, the length cost and the smoothness cost, all belong to [0, 1], and satisfy λ d + λ l + λ s = 1. According to the inclination on the project, set λ d = 0.5, λl = 0.3, λ s = 0.2, in the absence of beyond the ability of the shield to correct deviation, deviation, smoothness and length of the second cost, a great radius of curvature penalty factor to ensure the feasibility of the correction path, deviation once the limit is also immediately corrected deviation.

[0195] The step S6, the specific steps are:

[0196] Step S61, parameter setting: including particle swarm size N, the maximum number of iterations G m , acceleration coefficient c1 and c2, the extreme value of inertia weight ω min And ω max , the extreme value of particle position And s max , the maximum value of particle velocity V max , fractal parameter μ and disturbance intensity coefficient σ.

[0197] Step S62, population initialization: using Tent chaotic mapping to initialize the population, the dimension of individual is 1, particle position Indicates the mileage of the end point of the correction curve, the speed of individual is V k,G , where 0≤k≤N, G is the current iteration number of the algorithm, 0≤G≤G m .

[0198] Step S63, fitness value calculation: calculate the corresponding correction curve equation of each individual, and bring it into the comprehensive objective function, calculate the fitness value f k of each particle.

[0199] Step S64, individual optimal solution update: for each particle, compare its fitness value with the fitness of the individual historical optimal position p kb , if better, update p kb , otherwise maintain the individual optimal position unchanged.

[0200] Step S65, global optimal solution update: for each particle, compare its fitness value with the fitness of the global historical optimal position g b , if better, update g b , otherwise maintain the global optimal position unchanged.

[0201] Step S66, update the particle position and velocity of the next generation, the update formula is as follows:

[0202]

[0203] In the formula, c1 and c2 are acceleration coefficients, respectively controlling cognitive component and social component, r1 and r2 are uniform random numbers in the range of [0, 1], random disturbance is added to avoid the algorithm falling into local optimum. ω is an inertia weight, used to ensure the global convergence performance of the algorithm, the greater the value, the stronger the global convergence ability, otherwise the stronger the local convergence ability. N(0, 1) represents a standard Gaussian distribution random variable, and σ is a disturbance intensity coefficient. The inertia weight ω is adaptively adjusted, specifically:

[0204] The inertia weight is gradually reduced as the iteration proceeds, as shown in the following formula:

[0205]

[0206] Generally, ω max = 0.9, ω min = 0.4.

[0207] Step S67, boundary condition processing: when the position or speed of the particle exceeds the set value, the boundary condition processing strategy is used to limit the particle in the feasible search space, improve the search efficiency, as shown in the following formula:

[0208]

[0209] V max is the maximum speed, the value is 1.

[0210] When the speed exceeds the boundary, set it to the maximum value V max , when the position exceeds the boundary, the speed is reversed, and the position is updated, which mainly hopes to search around the boundary value to get a better solution.

[0211] Repeat steps S63-S67 until the maximum iteration number is reached, output the global optimal solution g b , get the optimal curve correction point and the offset curve equation. The flow chart is shown in Figure 4 .

[0212] As described in step S62, the Tent chaotic mapping is used to initialize the population, specifically:

[0213] The initial chaotic variable a1 is randomly generated in the interval [0, 1], and the chaotic sequence {a n}

[0214]

[0215] Where μ∈[0, 1] is a fractal parameter, controlling the symmetry of the mapping function, generally taking 0.5 to ensure uniform distribution of the sequence.

[0216] Based on the obtained chaotic sequence, the chaotic variable is mapped to the search space of the optimization problem through a linear transformation to obtain the initial population position of the particle swarm, as shown in the following formula:

[0217]

[0218] In the formula, represents the initial position of the kth particle.

[0219] As Figure 5 indicated, the shield tunneling machine deviation trajectory high-precision autonomous planning system provided by the embodiment of the present application specifically comprises:

[0220] The turning radius determination module is configured to determine the minimum turning radius of the shield tunneling machine through geometric constraints and mechanical constraints.

[0221] The coordinate calculation module is configured to perform coordinate calculation on the DTA to obtain the coordinates, direction and curvature of each point on the DTA.

[0222] The coordinate system conversion module is configured to convert the state of the shield tunneling machine in the Frenet coordinate system and the Cartesian coordinate system.

[0223] The deviation curve endpoint feasible region determination module is configured to determine the deviation curve endpoint feasible region through segment pre-layout constraints.

[0224] The comprehensive cost function solving module is configured to establish a comprehensive cost function in a weighted summation manner, and solve the comprehensive cost function by using an improved particle swarm optimization algorithm, and output a globally optimal deviation trajectory.

[0225] It should be noted that the embodiments of the present application can be realized by hardware, software or a combination of software and hardware. The hardware part can be realized by using special logic; the software part can be stored in a memory and executed by a suitable instruction execution system, such as a microprocessor or a specially designed hardware. Those skilled in the art can understand that the above-mentioned devices and methods can be realized by using computer executable instructions and / or included in processor control code, such as provided on a carrier medium, such as a magnetic disk, CD or DVD-ROM, a programmable memory, such as a read-only memory (firmware), or a data carrier, such as an optical or electronic signal carrier. The device of the present application and its modules can be realized by a hardware circuit, such as a very large scale integrated circuit or a gate array, a semiconductor, such as a logic chip, transistor, etc., or a programmable hardware device, such as a field programmable gate array, programmable logic device, etc., can also be realized by software executed by various types of processors, and can also be realized by a combination of the above-mentioned hardware circuit and software, such as firmware.

[0226] The above merely illustrates the specific embodiments of the present application, but the protection scope of the present application is not limited thereto, and any modification, equivalent replacement and improvement within the technical range disclosed by the present application and within the spirit and principle of the present application should be covered within the protection scope of the present application.

Claims

1. A high-precision autonomous planning method for shield machine deviation correction trajectory, characterized in that: The following steps are involved: S1, determine the minimum turning radius of the shield machine based on geometric constraints and mechanical constraints; S2, based on the coordinates of the starting and ending points of each curve segment in the DTA, the curve mileage and curvature characteristics, coordinate solution is performed to obtain the coordinates, direction and curvature of each point on the DTA; S3, in the Frenet coordinate system, a parameterized equation of the correction curve is established based on a quintic polynomial, and the shield machine state is converted between the Frenet coordinate system and the Cartesian coordinate system; S4, determine the feasible region of the end point of the correction curve based on the pre-layout constraints of the segment; S5, construct a comprehensive cost function with deviation degree, curve length and smoothness as indicators; S6, using an improved particle swarm optimization algorithm to solve the cost function and output a global optimal correction trajectory.

2. The method according to claim 1, wherein In step S3, the spatial position, heading and curvature of the starting and ending points of the shield machine are mapped from the Cartesian coordinate system to the Frenet coordinate system. Based on the quintic polynomial line type, a parameterized equation of the correction curve corresponding to the end point of the correction curve is established, and then back-projected into the Cartesian space to construct the actual excavation trajectory.

3. The method according to claim 1, wherein In step S4, the posture of the assembled segments of the previous ring is calculated through the real-time shield machine posture, shield tail gap and thrust cylinder stroke. The assembly posture of the segments to be assembled in the world coordinate system is derived using a continuous matrix transformation method, and the assembly deviation of the entire segment is evaluated. When the fitting deviation between all assembled ring segments and the DTA does not exceed the allowable deviation, it is determined that the correction curve can pass the pre-typesetting verification.

4. The method according to claim 1, wherein In step S1, the shield tail gap, the maximum propulsion cylinder stroke difference and the articulated structure limitation are calculated respectively to obtain multiple turning radius constraints. The maximum value is selected as the geometric limit radius, and the minimum turning radius under the computational mechanics of the shield statics model is introduced. Finally, the corrected radius is obtained by fusion and application of the safety factor as the actual feasible turning constraint.

5. The method according to claim 1, wherein The cost function in step S5 includes three target items: deviation cost item, curve length cost item and smoothness cost item; and a penalty item is set to impose penalties on solutions whose minimum curvature radius is lower than the threshold or whose deviation exceeds the upper limit. The comprehensive objective function is weighted by the weight coefficient and the sum of the items is output as the optimization target.

6. The method according to claim 1, wherein In step S6, a particle swarm optimization algorithm with a disturbance term and an adaptive adjustment strategy for the inertia weight is used to set the swarm size, acceleration factor, maximum number of iterations, search range, and speed limit. The particle position represents the mileage coordinates of the correction endpoint, and the particle fitness is calculated by a comprehensive cost function.

7. A system for autonomous planning of shield machine deviation correction trajectory, characterized in that: include: Turning radius calculation module, used to derive the minimum feasible turning radius of the shield machine based on structural and mechanical constraints; Coordinate solution module, used to analyze the spatial parameters of the DTA curve and generate the target reference line; State projection module, used to map states between Frenet and Cartesian coordinate systems; Curve verification module, used to optimize the feasible region of the correction curve endpoint based on the pre-layout deviation of the segment; The path optimization module is used to output the optimal correction trajectory based on the cost function and the improved particle swarm optimization algorithm.

8. The system according to claim 7, wherein: The path optimization module includes: The particle initialization subunit is used to generate the initial search solution set using the chaotic map; The cost calculation subunit is used to construct a cost function for each particle based on deviation, length and curvature; The solution updating subunit is used to iteratively update the particle state according to the particle's historical optimal position and the global optimal position until the optimal trajectory is output.

9. The system according to claim 7, wherein: The curve verification module includes: The posture calculation unit is used to deduce the relative posture of the previous ring of assembled segments and the shield machine based on the shield machine posture, thrust cylinder stroke and shield tail gap; The pose transfer unit is used to construct the pose transformation matrix between adjacent rings and accumulate the pose of each ring of segments to be assembled in the world coordinate system; The fitting evaluation unit is used to fit the full segment assembly with the correction curve and make tolerance judgment on the error.

10. The system according to claim 7, wherein: The final output of the path optimization module includes: the end point position of the correction curve, the correction curve parameter set, the pre-layout sequence of the segments and the deviation verification table. It supports importing the results into the shield construction information system for automated excavation control.

Images