Multimodal data dynamic sensing and intelligent precision control method for sludge solidification construction process

Abstract

This invention discloses a method for dynamic perception and intelligent precise control of multimodal data in the sludge solidification construction process, belonging to the field of intelligent control of engineering machinery. This method collects the global pose of the machine body, driving power parameters, and visual image sequences containing high-position markers. Utilizing a vision-inertial fusion algorithm for engineering machinery with kinematic constraints, it calculates the real-time three-dimensional trajectory and attitude of the mixing head underground while resisting sludge obstruction and vibration interference. Furthermore, it discretizes the mixing trajectory to construct the effective sweep volume of micro-element trajectory segments and, combined with adaptive soil resistance indicators, establishes a physical closed-loop control model based on the target unit volume grouting density. This invention abandons traditional linear empirical control logic and solves the problem of uneven distribution of solidifying agent caused by the tilting and variable speed motion of the mixing head by real-time matching of sweep volume and grouting energy, achieving microscopic calculation control and full-process traceability of the construction quality of concealed works.

Description

Technical Field

[0001] This invention belongs to the field of intelligent control of engineering machinery, and in particular to a method for dynamic perception and intelligent precise control of multimodal data in the sludge solidification construction process. Background Technology

[0002] In-situ solidification and improvement technology, as an efficient method for soft soil foundation treatment, is widely used in concealed engineering fields such as river dredging, silt solidification, and soft soil subgrade reinforcement. It is of great significance for improving land utilization and engineering safety. This technology uses a specially designed mixing head to forcibly mix the solidifying agent with the in-situ soil, causing a physicochemical reaction in the soil to form a composite foundation with a certain strength. It is one of the key construction methods in modern geotechnical engineering.

[0003] Current construction monitoring primarily employs an indirect positioning mode using GNSS fuselage positioning and boom tilt sensors, which calculates the coordinates of the mixing head based on the fuselage position. For grouting control, open-loop or semi-closed-loop adjustments are often made using linear empirical formulas based on depth or lowering speed. Existing technologies generally assume the mixing head is in an ideal vertical position and set fixed grouting flow rates or mixing speeds based on preset soil parameters.

[0004] However, existing technologies suffer from two major drawbacks under complex working conditions: a disconnect between perception and physical reality, and a lack of closed-loop control logic. The indirect geometric positioning method ignores the cumulative errors introduced by machine vibration and hydraulic arm deformation, and once the mixing head penetrates the ground, it is in a blind spot of posture perception, lacking the ability to continuously correct the actual underground trajectory, leading to deviations from the design axis. Traditional linear grouting control only considers one-dimensional depth changes and fails to calculate the nonlinear changes in the effective swept volume when the mixing head is tilted or moving at a non-uniform speed. When the mixing head is tilted, the cross-sectional area of ​​the actual disturbed soil increases in an elliptical shape. If grouting is still controlled by vertical advance, it will directly lead to insufficient grout density per unit volume of soil (leakage) or excessive energy input (waste), severely affecting the uniformity and integrity of the mixing quality. Summary of the Invention

[0005] The purpose of this invention is to provide a method for dynamic perception and intelligent precise control of multimodal data in the sludge solidification construction process, in order to solve one of the problems mentioned above in the existing technology.

[0006] Technical solution: A method for dynamic sensing and intelligent precise control of multimodal data in the sludge solidification construction process, comprising:

[0007] Real-time acquisition of global pose data of the excavator body, visual image sequence including the marking structure arranged on the mixing arm, and driving power parameters of the mixing drive shaft; wherein, the marking structure is arranged along the length of the mixing arm and the arrangement position is higher than the bottom of the mixing head by a set distance, so that the marking structure is within the visual range when the mixing head is submerged underground for operation;

[0008] Based on global pose data and visual image sequences, the real-time three-dimensional motion trajectory and attitude of the stirring head are calculated using visual-inertial fusion.

[0009] Based on the driving force parameters and the vertical motion velocity derived from the real-time three-dimensional motion trajectory, the soil resistance characteristic index that characterizes the current soil properties is calculated.

[0010] Based on the effective swept volume generated by the real-time three-dimensional motion trajectory, combined with the soil resistance characteristic index, the target grouting control volume at the current moment is calculated.

[0011] Adjust the execution parameters of the grouting equipment according to the target grouting control volume.

[0012] Beneficial effects: This invention solves the problem of uneven distribution of curing agent caused by tilting and variable speed of the mixing head, and realizes microscopic calculation control and full-process traceability of the construction quality of concealed works. Attached Figure Description

[0013] Figure 1 This is a flowchart illustrating the overall hardware architecture and processing flow of the multimodal data dynamic perception and intelligent precise control method for the sludge solidification construction process in this application embodiment.

[0014] Figure 2 This is a diagram illustrating the detailed calculation process of soil resistance characteristic indicators in the embodiments of this application.

[0015] Figure 3 This is a flowchart illustrating the steps of the physical control model in the embodiments of this application.

[0016] Figure 4 This is a three-dimensional structural diagram of the tracked hydraulic excavator in the embodiments of this application.

[0017] Figure 5 This is a side view of the tracked hydraulic excavator in an embodiment of this application.

[0018] Figure 6 This is a flowchart illustrating the steps of the central control unit in this application for evaluating construction quality and recording data.

[0019] In the picture: 1. Body; 2. Walking mechanism; 3. Upper arm; 4. Lower arm; 5. Mixing arm; 6. Mixing head; 7. Front-view wide-angle camera. Detailed Implementation

[0020] Example 1, such as Figure 1 As shown, this paper presents the overall hardware architecture and processing flow of a method for dynamic sensing and intelligent precision control of multimodal data during sludge solidification construction. A physical hardware environment was constructed, providing a foundation for the subsequent implementation of high-precision positioning and closed-loop control algorithms.

[0021] Step 101: Real-time acquisition of global pose data of the excavator body, including visual image sequences of the marking structure arranged on the mixing arm, and driving power parameters of the mixing drive shaft; wherein, the marking structure is arranged along the length of the mixing arm and the arrangement position is higher than the bottom of the mixing head by a set distance, so that the marking structure is within the visible range when the mixing head is submerged underground.

[0022] In this embodiment, the acquisition of global pose data primarily relies on an RTK (Real-Time Kinematic) positioning module and an IMU (Inertial Measurement Unit) module mounted on the top of the excavator. The RTK positioning module provides centimeter-level three-dimensional position coordinates (X, Y, F, G, I) of the excavator in an engineering coordinate system (typically ENU, i.e., the North-South coordinate system). RTK ,Y RTK Z RTK The IMU module is used to measure the fuselage's attitude angles (including roll, pitch, and yaw) and its three-axis angular velocities and three-axis accelerations in real time. The global pose data includes not only static position information but also high-frequency dynamic information reflecting fuselage vibration and tilt.

[0023] The acquisition of visual image sequences is achieved through a forward-looking wide-angle camera mounted on an external bracket of the cab. Specifically, this camera is configured to have a field of view covering the excavator's boom, arm, and mixing head. To address the technical challenge of not being able to directly observe the mixing head after it has penetrated underground, this embodiment employs a special marking strategy. Marking structures (such as high-contrast QR codes or reflective strips) are not placed on the mixing head body, but rather spaced along the length of the mixing arm. The "distance above the bottom of the mixing head" specifically refers to the vertical height difference h between the lowest marking structure and the bottom of the mixing head. offset The distance is set to a fixed value, such as 0.42 meters or longer (e.g., 2 to 5 meters), depending on the geometry of the maximum designed mixing depth and the camera's field of view. This setup ensures that when the mixing head is submerged deep underground, the marker structure located above the mixing arm remains exposed above the surface and within the camera's field of view, guaranteeing continuous observation and tracking of the marker structure.

[0024] The driving power parameters mainly include the real-time torque T and real-time speed n of the mixing drive shaft. This data can be directly acquired by torque and speed sensors installed on the hydraulic motor output shaft, or indirectly calculated by reading the hydraulic system pressure and flow values ​​from the excavator controller's local area network bus. All sensor data is aggregated to the central processor via the vehicle gateway for time synchronization and fusion processing.

[0025] Step 102: Based on global pose data and visual image sequences, use visual-inertial fusion to calculate the real-time three-dimensional motion trajectory and attitude of the stirring head.

[0026] This embodiment addresses the problem that relying solely on GPS or vision alone cannot accurately determine the location of the underground mixing head in complex construction environments. The system establishes a rigid body transformation relationship from the machine's global pose to the engineering coordinate system. Simultaneously, it uses the marker structure in the visual image sequence to calculate the local relative pose of the mixing arm relative to the machine. By cascading the global and local transformations and incorporating the geometric kinematic constraints of the excavator boom, the real-time three-dimensional spatial position P(t)=[X(t),Y(t),Z(t)] of the grout nozzle at the bottom of the mixing head in the engineering coordinate system is calculated. T And the spatial orientation R(t) of the mixing head. This process specifically employs a multi-frame optimized engineering machinery vision-inertial fusion algorithm to resist strong vibrations and mud obstruction interference during construction.

[0027] Step 103: Based on the driving force parameters and the vertical motion velocity derived from the real-time three-dimensional motion trajectory, calculate the soil resistance characteristic index that characterizes the current soil properties.

[0028] In this embodiment, the vertical motion velocity v down It is obtained by performing time differentiation or difference calculations on the Z-axis coordinates in the real-time three-dimensional motion trajectory. The soil resistance characteristic index R is a dimensionless physical quantity that comprehensively reflects the soil hardness, density, and construction difficulty at the current mixing depth. In specific calculations, the system considers not only the contributions of torque and rotational speed but also the vertical motion velocity. For example, with the same torque, a faster drilling speed indicates lower soil resistance; conversely, a slower speed indicates greater resistance. By constructing a functional relationship that includes torque, rotational speed, and velocity terms, soil properties can be calculated more accurately, providing a basis for subsequent grouting control.

[0029] Step 104: Based on the effective sweep volume generated by the real-time three-dimensional motion trajectory, and combined with the soil resistance characteristic index, calculate the target grouting control volume at the current moment.

[0030] This embodiment focuses on achieving refined construction. Unlike the extensive approach of existing technologies that only linearly adjusts the flow rate based on depth, this embodiment introduces the physical concept of effective swept volume. The system considers the trajectory of the mixing head per unit time as the swept volume and calculates the actual volume of soil ΔV occupied or disturbed by this swept volume in three-dimensional space. Then, based on the soil resistance characteristic indicators, the target grouting density ρ per unit volume required for this type of soil is determined. target (For example, how many liters of grout need to be injected per cubic meter of soil). This is determined using the formula Q. target =ΔV*ρ targetThe theoretically required grout flow rate for the current moment is calculated, i.e., the target grouting control volume. This method enables on-demand grouting, ensuring uniform distribution of the curing agent within the soil.

[0031] Step 105: Adjust the execution parameters of the grouting equipment according to the target grouting control volume.

[0032] In this embodiment, the central processing unit sends the calculated flow command to the grouting flow controller. The grouting equipment adjusts the frequency converter frequency of the grouting pump or the opening of the proportional flow valve in the pipeline according to the command, so that the actual output grout flow rate approximates the target grouting control quantity. Simultaneously, the system collects the actual grouting volume feedback from the flow meter in real time, forming a closed-loop control to ensure that the construction quality meets design requirements.

[0033] Example 2 describes a robust trajectory calculation based on vision-inertial fusion. It utilizes a multibody kinematics chain model combined with a prediction-correction mechanism and an outlier removal mechanism to achieve high-precision and robust three-dimensional trajectory tracking under harsh working conditions such as mud obstruction and machine shaking.

[0034] Step 201: Construct a multibody kinematics model of the excavator, including the fuselage, boom, arm and mixing arm.

[0035] In this embodiment, the multibody kinematics chain model is a mathematical model describing the geometric connections and motion transmission relationships between the rigid body components of the excavator. Specifically, the system defines the fuselage coordinate system O. body Large arm coordinate system O boom Forearm coordinate system O stick and the coordinate system O of the stirring arm stir The transformation relationships between the coordinate systems can be represented by rotation matrices and translation vectors. Although the hinge angles between the boom, forearm, and fuselage are usually unknown (no angle sensors are installed), this embodiment treats them as unknown state variables to be solved by the image acquisition and processing system, or directly uses the end pose calculated based on image data to replace the intermediate joint angle calculation.

[0036] During coordinate transfer, the influence of fuselage attitude on the final coordinates must be considered. Global pose data includes the fuselage rotation matrix R relative to the engineering coordinate system. body_global For the Z-axis coordinate (depth) of the stirring head in the engineering coordinate system, its calculation formula cannot simply be the sum of the ship's altitude and relative depth, but should employ a full matrix transformation: P global =P RTK +R body_global *P visual_body , where P global The coordinates of the mixing head in the engineering coordinate system; P RTK These are the coordinates of the positioning module in the engineering coordinate system; P visual_bodyIt is the coordinate vector of the mixing head relative to the machine body, calculated based on the pose of image data. This can effectively eliminate the coupling error of roll angle or pitch angle on depth measurement caused by the machine body operating on rough terrain.

[0037] Based on the current frame's visible marker set, and according to the boom position where the visible marker is located, the rigid body components (i.e., the observable parts of the boom, forearm, and stirring arm) participating in the kinematic chain calculation are determined. Then, using the pose information of the visible markers, combined with the geometric connections between the rigid body components, the stirring arm pose and the initial pose estimate P of the stirring head are obtained. 0 (t), R 0 (t).

[0038] Using the estimated values ​​P(t-Δt) and R(t-Δt) from the previous moment, and the motion ΔP calculated by the IMU, IMU ,ΔR IMU Construct kinematic prediction values:

[0039] P'(t)=P(t-Δt)+ΔP IMU ;

[0040] R'(t)=R(t-Δt)·ΔR IMU ;

[0041] P is calculated from the pose of a single frame based on image data. 0 (t),R 0 (t) is used as the observed value, and together with the predicted value, it constitutes a Kalman filter nonlinear least squares optimization problem: min P(t),R(t) |P(t)-P'(t)|W p 2 +∣∣R(t)-R'(t)∣∣W r 2 +λ·φ indicates the reprojection error, where W p 2 W is the location prediction covariance weight matrix. r 2 Let be the attitude prediction covariance weight matrix; λ be the image observation weight coefficients; and φ be the loss function robust to gross residuals.

[0042] The reprojection error of each visible marker is evaluated, and suspicious markers with errors exceeding the threshold are automatically removed (they may be contaminated by mud or partially obscured).

[0043] When the number of available markers falls below a preset lower limit, the system switches to a short-term inertial trajectory mode based on IMU prediction and reconverges when the markers become visible again.

[0044] Step 202: Identify the currently visible signage structure from the visual image sequence, and calculate the local observation pose of the signage structure relative to the forward-looking wide-angle camera based on perspective geometry.

[0045] In this embodiment, the system processes each frame of image to extract high-contrast feature points of the marker structure. Using the Perspective-n-Point (PnP) algorithm or similar computer-based image data pose calculation methods, combined with a pre-calibrated camera intrinsic parameter matrix and the geometric dimensions of the marker structure, the rotation and translation vectors of the marker structure coordinate system relative to the camera coordinate system are calculated. Since the marker structure is positioned on the stirring arm, and the relative positional relationship between the stirring head and the marker structure (including the height offset h) is considered... offset Since the position of the nozzle of the mixing head relative to the camera is known and fixed, the local observation pose of the nozzle relative to the camera can be further calculated.

[0046] Step 203: Based on the multibody kinematics chain model, kinematic transfer and coordinate transformation are performed on the local observed pose and global pose data to obtain the initial state value of the stirring head in the engineering coordinate system.

[0047] After obtaining the local observation pose, the system transforms it to the fuselage coordinate system to obtain P. visual_body Combined with the fuselage position P measured by RTK RTK and the fuselage attitude R measured by the IMU body_global The absolute position and attitude of the stirring head in the engineering coordinate system are calculated using the rigid body transformation formula, and used as the initial state value X at the current moment. init Although the initial value contains absolute position information, pose calculation based on image data in a single frame is easily affected by factors such as changes in illumination and motion blur, resulting in high-frequency noise or jumps. Therefore, subsequent steps are needed for smoothing optimization.

[0048] Step 204: The initial state value is smoothed and optimized using the motion continuity constraint in the time dimension to obtain the real-time three-dimensional motion trajectory and attitude of the stirring head.

[0049] This embodiment utilizes the principle of continuity in physical motion, meaning that the position and velocity of the stirring head will not change abruptly within a very short time, to filter noisy initial state values. This is typically achieved through Kalman filtering (KF), extended Kalman filtering (EKF), or a graph-optimized sliding window smoothing algorithm.

[0050] Step 205: Obtain the estimated state of the stirring head at the previous moment and extract the inertial measurement data contained in the global pose data.

[0051] In this embodiment, the system maintains the state vector X. kIt includes the three-dimensional position, three-dimensional velocity, and attitude quaternions of the stirring head. The state estimate X from the previous moment... k-1 Read from historical records. Inertial measurement data includes acceleration *a* and angular velocity *ω* sensed by the fuselage and boom systems. This inertial data features high-frequency sampling (e.g., 100Hz-200Hz) and high short-time accuracy, making it suitable for short-term state prediction.

[0052] Step 206: Based on inertial measurement data and multibody kinematics chain model, perform differential deduction on the estimated state of the stirring head at the previous moment to generate the prior predicted state at the current moment.

[0053] In this embodiment, the system uses kinematic equations to update the state over time. Specifically, the position prediction follows P... k_pred =P k-1 +v k-1 *Δt+0.5*a*Δt 2 The physical laws, in which P k-1 The position optimized from the previous moment; v k-1 'a' represents the estimated velocity at the previous moment; 'a' represents the estimated acceleration at the current moment; and 't' represents the time interval between two image observations. Attitude prediction is updated by integrating the attitude quaternion from the previous moment using the angular velocity ω. Through this differential derivation, the system provides a high-confidence prediction of the stirring head's state at the current moment, i.e., the prior prediction state, even before the visual image has been fully processed.

[0054] Step 207: Calculate the single-point reprojection error of each identified feature point in the visual image sequence.

[0055] In actual construction, mud splashes often obscure some markings, leading to false detections by image recognition algorithms. To address this issue, the system introduces a reprojection error verification mechanism. The system utilizes prior predicted states and the camera's projection model to virtually project all known marking feature points on the mixing arm onto the current image plane, obtaining a set of virtual pixel coordinates u. proj Simultaneously, the image processing algorithm extracts the pixel coordinates u of the feature points actually observed in the current frame. obs Single-point reprojection error e i Defined as the Euclidean distance between the actual observation point and the corresponding virtual projection point, i.e., e i =||u obs_i -u proj_i ||.

[0056] Step 208: Compare the single-point reprojection error with a preset confidence threshold.

[0057] The confidence threshold is a preset pixel distance value, such as 3.0 pixels or 5.0 pixels. This threshold reflects the system's tolerance to image observation noise.

[0058] Step 209: When the single-point reprojection error of a certain identification feature point exceeds the confidence threshold, the identification feature point is determined to be an outlier point that is contaminated or occluded. The outlier point is removed when constructing the observation residual equation, and the remaining effective feature points are used to minimize the solution.

[0059] If the reprojection error of a feature point is greater than that of other points or exceeds a threshold, it indicates that the point is highly likely to be a mismatched point caused by mud obscuring it, or a false corner point caused by light spots. The system will mark such points as outliers and remove them from the observation set. The outlier removal mechanism based on dynamic priors improves the robustness of the system under harsh operating conditions.

[0060] Step 210: Construct the observation residual equation, which characterizes the reprojection error between the feature points projected onto the image plane from the prior predicted state and the actual identified feature points in the visual image sequence.

[0061] After removing outliers, the system constructs the global observation residual equation using the remaining valid feature points. This equation is essentially a nonlinear least squares problem, and its goal is to find the optimal state estimate that minimizes the sum of squared projection errors in that state.

[0062] Step 211: Correct the prior predicted state by minimizing the observation residual equation to obtain the real-time three-dimensional motion trajectory and attitude at the current moment.

[0063] The system employs the Gauss-Newton method or the Levenberg-Marquardt algorithm to iteratively solve the aforementioned residual equations, or performs a measurement update step within the extended Kalman filter framework. Through this process, the prior predicted state (primarily determined by inertial data) is corrected by the image observation data (the effective portion), and the fused result is a high-precision, smooth, and robust real-time three-dimensional motion trajectory and attitude. This result is not only used to record the construction trajectory but also serves as input for subsequent calculations of the swept volume and control of the grouting volume.

[0064] Example 3, such as Figure 2 As shown, this paper describes the refined calculation process of soil resistance characteristic indicators. It includes a basic resistance model incorporating a dip angle penalty term, introduces an adaptive weight optimization method based on the coefficient of variation, and a weight matching strategy based on an offline library, addressing the problem that a single fixed formula is insufficient to adapt to complex and variable geological conditions.

[0065] Step 301: Extract the stirring torque and stirring speed from the driving power parameters, and analyze the real-time tilt angle of the stirring head relative to the vertical direction from the attitude data of the real-time three-dimensional motion trajectory.

[0066] In this embodiment, the real-time tilt angle γ of the stirring head relative to the vertical direction is determined by comparing the stirring head axis vector k. axis It is calculated using the gravity direction vector in the global coordinate system (usually the negative Z-axis). The specific calculation formula can be the inverse cosine of the dot product of two unit vectors. This angle directly reflects the perpendicularity deviation of the stirring.

[0067] Step 302: Normalize the stirring torque, stirring speed and vertical motion speed to construct torque component, speed component and velocity component respectively.

[0068] To eliminate the influence of different physical dimensions, each parameter needs to be normalized. For example, the torque component C T This can be expressed as the ratio of the real-time torque T to the soft soil calibration reference torque T0, i.e., C T =T / T0. Similarly, the speed component C n =n0 / max(n,n min It is usually inversely proportional to the resistance, so it is taken as the reciprocal. The velocity component C v =v0 / max(|v down |,v min This treatment makes the components numerically comparable, and all show a positive correlation with soil resistance. Among them, The recommended value for the minimum effective drilling speed threshold is [value to be filled in]. (i.e., approximately 0.6 m / min). If the speed is lower than this value, the stirring head is considered to be in a paused or quasi-stationary state, and the normal control mode should be exited. The recommended value is the minimum effective speed threshold. .when or When the system enters a static / low-speed protection mode, it locks the current grouting flow rate to the calculated value from the previous normal time point and displays a message on the interface indicating that the mixing head speed is too low and the grouting flow rate is locked. Normal calculation resumes once the speed or rotation speed returns to above the threshold.

[0069] Step 303: Obtain the preset maximum allowable tilt angle threshold.

[0070] Maximum permissible tilt angle threshold γ max This is a quality control red line set according to construction specifications, such as 5° or 8°. When the actual inclination angle approaches or exceeds this value, it usually indicates serious potential construction quality problems.

[0071] Step 304: Calculate the ratio of the real-time tilt angle to the maximum allowable tilt angle threshold.

[0072] The system calculates the dimensionless ratio. γ =γ / γmax .

[0073] Step 305: The square or higher power of the ratio is used as the tilt angle penalty term to achieve non-linear growth of the soil resistance characteristic index when the real-time tilt angle approaches the maximum allowable tilt angle threshold.

[0074] Step 306: Construct a tilt angle penalty term based on the real-time tilt angle. The tilt angle penalty term is configured to increase non-linearly as the real-time tilt angle increases.

[0075] Describe the tilt angle penalty term P γ The construction logic is as follows. Preferably, a quadratic relation is used for construction: P γ =(γ / γ max ) 2 The square term is introduced because when the inclination angle is small, its impact on the resistance coefficient is negligible; however, as the inclination angle increases, the lateral compression resistance and equipment wear risk caused by the inclination increase exponentially. Through nonlinear penalty, the calculated resistance value under inclination conditions can be amplified, forcing the subsequent control system to respond (such as increasing the grouting volume to compensate for uneven mixing, or prompting the operator to correct the deviation), truly reflecting the comprehensive construction difficulty including attitude factors.

[0076] Step 307: Establish a sliding time window that includes historical data of a set length prior to the current time.

[0077] To capture the dynamic characteristics of the soil layers, the system maintains a first-in, first-out data queue with a time length Δt. window The sampling frequency can be determined based on the actual sampling frequency and the rate of soil layer change; in this embodiment, it is set to 200 milliseconds. This window contains the most recent torque, rotational speed, and velocity sequence.

[0078] Step 308: Calculate the time series variation coefficients of stirring torque, stirring speed and vertical motion speed within the sliding time window.

[0079] The coefficient of variation (CV) is defined as the ratio of the standard deviation to the mean, i.e., CV = σ / μ. The system calculates the coefficient of variation CV for the torque series within the window. T Coefficient of variation (CV) of the rotational speed sequence n and the coefficient of variation (CV) of the velocity sequence v These statistics reflect the degree of fluctuation or sensitivity of various physical quantities in the current soil layer.

[0080] Step 309: Pre-construct a soil layer feature model library, which contains a variety of typical soil types.

[0081] To support dynamic weight adjustment, this embodiment pre-constructs an offline database. This database contains characteristic descriptions of typical geological types such as silt, soft soil, clay, and stiff plastic clay.

[0082] Step 310: For each typical soil type, pre-calibrate its corresponding torque sensitivity, rotational speed sensitivity, and velocity sensitivity, and establish a preset benchmark combination of weighting coefficients for that soil type.

[0083] For example, in silt layers with high water content, torque changes may not be significant, but drilling speed is extremely sensitive. Therefore, the preset weighting combination is (torque weight 0.3, rotation speed weight 0.2, speed weight 0.4, inclination angle weight 0.1). For stiff plastic clay, torque is extremely sensitive, and the preset weighting combination is (torque weight 0.65, rotation speed weight 0.2, speed weight 0.1, inclination angle weight 0.05). These combinations are obtained in advance through analysis of a large amount of historical construction data.

[0084] Step 311: Identify the sensitivity of the current soil layer to each parameter based on the time series coefficient of variation, and dynamically adjust the weight coefficients used for weighted summation accordingly, so that parameters with larger coefficients of variation within the sliding time window have larger weight coefficients.

[0085] Step 312: When dynamically adjusting the weight coefficients used for weighted summation, the current soil layer is matched to the typical soil type in the soil feature model library based on the current time series variation coefficient, and the preset benchmark combination corresponding to the typical soil type is called as the adjustment benchmark.

[0086] Describes the real-time update mechanism for weights. Based on the currently calculated coefficient of variation eigenvector, the system searches the model library for the best-matching soil category (e.g., by calculating Euclidean distance) and fine-tunes the weights based on the baseline weight combination for that category, taking into account the magnitude of the coefficient of variation. For example, if the currently detected coefficient of variation for torque suddenly increases, the system automatically increases the weight of the torque component in the total resistance formula.

[0087] Step 313: Weighted summation of torque component, rotational speed component, velocity component and tilt angle penalty term to obtain soil resistance characteristic index.

[0088] The system mainly provides three calculation schemes for soil resistance characteristic indicators: adaptive weight method, fixed weight method, and soil layer model switching method.

[0089] Adaptive weighting method: The system uses dynamically adjusted weight coefficients w T ,w n ,w v ,w γ Calculate the soil resistance characteristic index R=w T *CT +w n *C n +w v *C v +w γ *P γ This index R can adaptively reflect the actual construction impedance under different geological conditions.

[0090] Fixed-weight method: As an alternative or simplified implementation of the above scheme, in scenarios with limited computing resources or simple geological conditions, a linear formula with fixed weights can also be used for calculation. For example, directly setting... Although its adaptability is not as good as that of adaptive algorithms, it can still provide an effective resistance reference in conventional soft soil construction.

[0091] Soil layer model switching method: As an alternative to the adaptive weighting method, a method of soil layer segmentation identification and resistance coefficient model switching can be adopted.

[0092] Based on engineering geological survey data, multiple soil layer resistance coefficient models are pre-established. Each model corresponds to a typical soil layer type and includes the characteristic weight combination of that soil layer.

[0093] Extremely soft sludge model: w T =0.3,w n =0.2,w v =0.4,w γ =0.1;

[0094] Normal soft soil model: w T =0.5,w n =0.3,w v =0.2,w γ =0.15;

[0095] Harder clay model: w T =0.6,w n =0.25,w v =0.15,w γ =0.15;

[0096] Hard clay model: w T =0.65,w n =0.2,w v =0.1,w γ =0.2;

[0097] During implementation, the average torque μ within the sliding window is used as the basis. T Average rotational speed μ n Average drop speed μ v Construct the feature vector F=[μ T / T0,n0 / μn ,v0 / μ v ].

[0098] Calculate the Euclidean distance between the feature vector F and the feature center of each soil layer model, and select the soil layer model with the smallest distance as the current identification result.

[0099] To avoid abrupt changes in the resistance coefficient due to frequent jumps in soil layer identification results, an exponential smoothing switch is introduced:

[0100] Let w be the weight used in the previous time step. prev The weight corresponding to the currently identified soil layer is w. new The actual weight used is: w actual =β×w new +(1-β)×w prev , where β is the smoothing coefficient (ranging from 0.1 to 0.3, preferably 0.2).

[0101] The soil resistance coefficient is calculated using smoothed weights:

[0102] R'=w T,actual ×T / T0+w n,actual ×n0 / n+w v,actual ×v0 / ∣v down ∣+w γ,actual ×(γ / γ max ) 2 .

[0103] Example 4, such as Figure 3 The diagram illustrates the physical control model. The construction control logic is upgraded from traditional experience-based linear parameter adjustment to a dual closed-loop control system based on effective swept volume and physical density. By constructing a dynamic geometric model of the mixing head's movement underground, the actual disturbed soil volume per unit time is accurately calculated, and the grout flow rate and mixing energy are matched accordingly, solving the problem of uneven grouting caused by variations in mixing speed, path inclination, or rotational operation.

[0104] Step 401: Discretize the real-time three-dimensional motion trajectory in the time dimension to obtain the micro-element trajectory segment at the current moment.

[0105] In this embodiment, the central processing unit samples the continuous trajectory at a fixed control period Δt (e.g., 50 milliseconds or 100 milliseconds). The infinitesimal trajectory segment starts from the bottom position P of the stirring head at the previous control moment. k-1 and the current position P k The spatial vector ΔP formed kThis vector contains not only the vertical advance component Δz, but also the horizontal displacement component Δxy. For example, when an excavator performs point mixing, Δxy is close to zero; however, when performing lateral movement mixing similar to milling, Δxy cannot be ignored.

[0106] Step 402: The real-time tilt angle of the stirring head is analyzed from the spatial pose of the micro-element trajectory segment.

[0107] Step 403: Construct a virtual cylinder based on the designed stirring radius of the stirring head, and perform spatial sweep modeling of the virtual cylinder along the micro-element trajectory segment to obtain a three-dimensional spatial geometry.

[0108] Describe the construction process of the swept geometry. A virtual cylinder refers to a cylinder with a base radius equal to the designed stirring radius r. design A geometric unit with an infinitesimally small height. When the stirring head moves along the infinitesimal trajectory segment, the trajectory swept by this virtual cylinder in three-dimensional space forms an oblique cylinder or a more complex swept body.

[0109] Step 404: Correct the projected area of ​​the three-dimensional spatial geometry on the horizontal plane using the real-time tilt angle, so that the projected area increases with the increase of the real-time tilt angle, and calculate the corrected effective sweep volume accordingly to compensate for the expansion of the soil disturbance range caused by the tilt of the mixing head.

[0110] Step 405: Based on the geometric properties of the stirring head and the spatial pose of the micro-element trajectory segment, construct the three-dimensional spatial geometry swept by the stirring head at the current moment, and calculate the volume of the three-dimensional spatial geometry as the effective sweep volume; wherein, the geometric properties include at least the design stirring radius of the stirring head.

[0111] The calculation logic for the effective swept volume ΔV is described. Under ideal vertical conditions, the swept volume is simply the base area multiplied by the vertical advance. However, when the mixing head has an inclination angle γ, or when there is horizontal displacement, the actual disturbed soil cross-section is no longer a standard circle, but rather an ellipse or a composite shape with a larger projected area. This embodiment uses the projected area correction method to calculate the equivalent cross-sectional area S. eq Specifically, the correction formula can be expressed as follows: ;

[0112] in, The projection displacement components of the current infinitesimal trajectory segment onto the horizontal plane are given by the position vector. The horizontal components are obtained as follows:

[0113]

[0114] The physical meaning of the second term is: when the stirring head is displaced in the horizontal direction, the area (width) of the additional rectangular strip swept by the cylinder along the horizontal direction. ,length ).

[0115] The corresponding effective sweep volume correction is:

[0116]

[0117] When the stirring head only performs pure vertical movement ( When it degenerates into:

[0118]

[0119] when and It degenerates into the volume of a standard cylinder. It aligns with intuition.

[0120] Step 406: Based on the soil resistance characteristic index, determine the target grouting density per unit volume required for the current soil.

[0121] In this embodiment, the system determines the grouting requirements corresponding to the current soil type based on the soil resistance characteristic index R, using a lookup table or mapping function. The target grouting density per unit volume is ρ. target (Unit: liters / cubic meter) reflects the volume of solidifying agent added to each cubic meter of in-situ soil as required by the design. Generally, the greater the resistance R, the harder or more cohesive the soil, and the higher the required solidifying agent admixture ratio, i.e., ρ. target It should be increased accordingly.

[0122] The lookup table or mapping function can be pre-established by those skilled in the art based on construction design specifications and engineering geological survey data. Specifically, based on the design solidifier admixture ratio for different soil types, the soil resistance characteristic index R and the target unit volume grouting density ρ can be established. target The correspondence between them. Those skilled in the art can determine this correspondence based on the design requirements of the actual project.

[0123] Step 407: Calculate the target grouting control amount based on the effective swept volume and the target unit volume grouting density, so that the ratio of the target grouting control amount to the effective swept volume is close to the target unit volume grouting density.

[0124] The system is based on the physical equation Q target =(ΔV / Δt)*ρ target Calculate the instantaneous grout flow rate required at the current moment, i.e., the target grouting control quantity. The physical essence of this control logic is: regardless of whether the mixing head is moving quickly in and out, slowly in and out, or inclined during drilling, as long as the volume ΔV it sweeps is determined, the system ensures that the total amount of grout injected into that volume remains constant (i.e., constant density) by adjusting the flow rate Q. This eliminates the common quality problems of excessive grouting due to slow speed or insufficient grouting due to fast speed.

[0125] Step 408: Calculate the energy input value of the stirring head within the current micro-element trajectory segment based on the stirring torque, stirring speed, and micro-element time interval.

[0126] In addition to controlling the grouting volume, this embodiment also introduces a parallel energy control loop. The system calculates the work done by the mixing drive system on the soil within a time interval Δt, i.e., the energy input value E. in =T*2π*(n / 60)*Δt.

[0127] Step 409: Calculate the real-time energy density per unit volume based on the ratio of the energy input value to the effective swept volume.

[0128] The system calculates the energy density index E. density =E in / ΔV. This index reflects how much mechanical mixing energy a unit volume of soil receives, and is a physical quantity used to evaluate the uniformity of mixing and the effect of cutting and crushing.

[0129] Step 410: Compare the real-time energy density per unit volume with the preset target energy density, and adjust the rotational speed or torque output of the mixing drive shaft according to the comparison result to ensure that the mixing energy obtained per unit volume of soil meets the preset standard.

[0130] The target energy density can be predetermined based on the requirements for mixing uniformity in the construction design specifications, combined with the soil type and the geometric parameters of the mixing head. Those skilled in the art can determine a suitable target energy density value based on the actual engineering design requirements, either through on-site mixing calibration or by conversion from the Binder Reaction Number (BRN) index. If the real-time energy density E... density If the drilling speed is lower than the design threshold (e.g., when the soil is too hard and the drilling speed is too fast), the system will issue a command to reduce the lowering speed of the excavator's hydraulic arm or prompt the operator to reduce the drilling speed, thereby increasing the mixing time and energy input per unit volume of soil to prevent the soil from being partially or unevenly mixed.

[0131] As an alternative to the aforementioned advanced control method based on volume models, in embodiments lacking high-precision trajectory calculation capabilities or with limited sensor configuration, the linear empirical formula from the original technical solution can also be used for control. Specifically, the target flow rate Q actual It can be calculated as Q design *(0.6+0.008*R)*(z cur / z design )*(1+k tilt *sin(γ)). Although this formula lacks a rigorous volumetric physical meaning, it can still provide an empirical control method with tilt compensation in conventional engineering.

[0132] Step 411, 3D trajectory discretization and local disturbance volume estimation: Based on the original 3D positioning, for each sampling interval Δt:

[0133] Obtain the current spatial position P(t) of the stirring head and the previous time P(t-Δt);

[0134] Calculate the displacement Δs on the horizontal plane. h (t) and vertical advance Δz(t);

[0135] By combining the stirring radius and tilt angle, the r value is improved. actual Formulas can be used to establish more reasonable local volume increments ΔV(t), for example: ΔV(t)≈πr eq 2 (t)·∣Δz(t)∣, where r eq (t) is no longer simply represented by r design cosγ+L stir sinγ, but r eq The improvement approximates the spatial path of the stirring head within Δt as a line segment; a swept volume is constructed around this line segment, equivalent to a local cylinder or elliptical cylinder; r is defined. eq (t) is the equivalent radius of the swept section on the horizontal plane: r eq (t)=sqrt A sweep (t) / π, where A sweep (t) Calculated based on the rolling path length and the design blade radius.

[0136] r eq It is no longer a simple linear superposition of cosγ and sinγ, but is linked to the actual swept area in space.

[0137] Local energy density / slurry density index E d (t),Q d (t), for each Δt:

[0138] Calculate the actual shotcrete volume ΔV during this time period. grout (t) (integrated by the flow meter);

[0139] The work done by torque or the equivalent energy consumption is calculated as E(t)≈T(t)·ω(t)·Δt;

[0140] Calculate the density and energy density Q of the slurry per unit volume in a localized area. d (t)=ΔV grout (t) / ΔV(t),E d (t)=E(t) / ΔV(t).

[0141] Target density curve and closed-loop control:

[0142] Based on the geological survey and design requirements, the target slurry density per unit volume, Q, is given in advance. d , design (z) (varies with depth);

[0143] Or target energy density E d , design (z);

[0144] In real-time control: the current Q... d (t) and Q d , design (z cur )Compare;

[0145] If Q d If (t) is lower than the target, increase the shotcrete flow rate in the next Δt.

[0146] If Q d If (t) is significantly higher than the target, the amount of sprayed grout can be appropriately reduced to avoid waste and local over-spraying.

[0147] Q at this time actual The calculation is no longer just an empirical linear correction, but can be written as:

[0148] Q actual (t+Δt)=Q actual (t)+K q ·(Q d,design (z cur )-Q d (t))+K r ·(R(t)-R target (z cur ))+Kγ·f(γ(t)); where:

[0149] K q ·(Q d,design (z cur )-Q d (t) is the closed-loop density of slurry per unit volume;

[0150] K r ·(R(t)-R target (z cur This is a correction for the deviation of the drag coefficient (the original R concept is used as a secondary factor).

[0151] Kγ·f(γ(t)) uses f(γ) (such as sinγ or γ²) to compensate for the increase in volume due to the tilt.

[0152] Mass score based on volume-energy density:

[0153] After the mixing construction at a certain point is completed, the analysis is no longer based solely on the average soil resistance characteristic index R.avg Stirring time t stir Instead, it calculates Q within the stirring range. d (t), E d Distribution of (t) (mean, variance, minimum);

[0154] Define a new rating:

[0155] Score * =100-α·σ(Q d ) / Q d,design -β·σ(Ed) / E d,design -λ1·∣1-Q' d / Q d,design |-λ2⋅(maximum tilt factor).

[0156] The score reflects whether the distribution of slurry and energy within the overall volume is uniform and sufficient, rather than simply a few linear deviations.

[0157] Example 5 describes a predictive compensation mechanism for the large lag characteristics of engineering systems, and a post-construction quality evaluation method based on density distribution.

[0158] Step 501: Based on the soil resistance characteristic index calculated at multiple consecutive time points, calculate the rate of change of the soil resistance characteristic index with time.

[0159] Due to the physical response delay in hydraulic systems and grouting pipelines, directly controlling grouting based on the current resistance often results in a lag. This embodiment calculates the first derivative or differential rate of change of the soil resistance characteristic index dR / dt=(R k -R k-N ) / (N*Δt). This rate of change reflects the trend of soil hardness changes; for example, dR / dt>0 indicates that the soil is entering a hard soil layer.

[0160] Step 502: Based on the rate of change and the preset system response delay time of the grouting equipment, predict the estimated soil resistance characteristic index for the next execution time after the current time.

[0161] The system uses a linear extrapolation model to predict the resistance value at future times: R pred =R curr +(dR / dt)*τ d Among them, the preset system response delay time τ d It is a constant calibrated based on the valve response time of the grouting pump and the transmission time of the grout in the pipeline, with a typical value of 0.5 seconds to 1.5 seconds.

[0162] Step 503: Use the estimated soil resistance characteristic index to perform advance compensation correction on the target grouting control volume to obtain the compensated execution command, and send the compensated execution command to the grouting equipment.

[0163] The system uses the predicted resistance value R pred Replace the current resistance value R curr Substitute the density into the formula to calculate the flow command after advance compensation. For example, when it is predicted that the system will enter a hard soil layer (requiring a large amount of slurry), it will advance the flow command by τ. d Increased flow rate over time offsets pipeline transmission delays, ensuring that the increased slurry reaches the nozzle precisely when the mixing head actually contacts the hard soil.

[0164] Step 504: Calculate the actual unit volume grouting density corresponding to each micro-element trajectory segment during the entire construction process. The actual unit volume grouting density is calculated by the ratio of the actual grouting volume to the effective sweeping volume.

[0165] Step 505: Calculate the distribution dispersion of the actual unit volume grouting density relative to the target unit volume grouting density.

[0166] After the mixing construction at a certain point is completed, the system no longer only focuses on whether the total grouting volume meets the standard, but also retrospectively examines the microscopic quality of the mixing point depth. The system extracts the actual grouting density data ρ at each recording moment. actual (t)=Q measured (t) / ΔV(t), and calculate its value relative to the design target value ρ. target The dispersion index. The dispersion of a distribution can be characterized by the standard deviation σ or the coefficient of variation CV.

[0167] Step 506: Generate the construction quality evaluation result of the mixing construction at this point based on the distribution dispersion, reflecting the uniformity of the distribution of the curing agent inside the solidified whole during the mixing construction.

[0168] If the dispersion is less than the preset threshold (e.g., 0.1), the point is judged to be uniformly stirred and the quality grade is excellent; if the dispersion is large, even if the total grouting volume meets the standard, it indicates that there is local over-mixing and local under-mixing (the "candied hawthorn" phenomenon), and the system will evaluate it as unqualified or require re-mixing.

[0169] As a simplified alternative evaluation method, a linear scoring formula based on weighted deductions can be used when data dimensions are insufficient. For example, Score = 100 - 30 * |R| avg -R target | / R target -20*(γ max / 5°). This method is simple to calculate and provides an intuitive percentage score, making it easy for construction workers to make quick judgments.

[0170] Example 6 describes the interaction logic of the feedback terminal. In this example, the feedback terminal in the driver's cab provides a visual monitoring interface.

[0171] The interface displays the calculated 3D trajectory of the mixing head in real time and shows the current cutting depth in the form of a side profile. To visually represent soil layer information, the system uses soil type to render the trajectory lines in different colors: for example, it displays blue when drilling through a silt layer with R < 30, and red when drilling through a hard soil layer with R > 70.

[0172] In addition, the interface integrates an alarm mechanism. When the system detects the real-time energy density E... density If the value remains below the threshold for an extended period, or if the real-time inclination angle γ exceeds 8 degrees, the corresponding indicator on the interface will flash a red warning light, accompanied by a buzzer alarm, prompting the operator to adjust the machine's attitude or reduce the drilling speed. This intuitive, WYSIWYG feedback reduces the operator's cognitive load and ensures the effective implementation of complex algorithms in actual construction.

[0173] Example 7 describes the detailed process of improving the robustness and smoothing of the algorithm.

[0174] Step 701: Construct feature vectors for soil layer model matching.

[0175] The system constructs the feature vector F in the following way: Within the sliding time window, the mean torque μ is calculated. T Average rotational speed μ n and the average drop speed μ v Subsequently, normalization is performed using the calibration values ​​of each parameter to construct a three-dimensional feature vector F=[μ T / T0,n0 / μ n ,v0 / μ v The vector F represents the centroid location of the current soil layer in the feature space. The system calculates the Euclidean distance between this vector F and the centroid vectors of each typical soil layer in the model library; the soil layer type with the smallest distance is identified.

[0176] Step 702: Perform exponential smoothing on the dynamically adjusted weighting coefficients.

[0177] To prevent abrupt and drastic fluctuations in the drag coefficient R due to sudden changes in soil layer identification results (e.g., a sudden change from silt to clay), the system introduces a smooth switching mechanism. Specifically, let w be the weight combination used in the previous control cycle. prev The new weight combination calculated based on the coefficient of variation or model matching is w. new The weight w actually used in resistance calculation actual Update using the following formula: w actual =β*w new +(1-β)*wprev Wherein, β is the smoothing coefficient, preferably ranging from 0.1 to 0.3 (e.g., 0.2). It acts as a low-pass filter, ensuring a smooth transition for the control system. The selection of the smoothing coefficient β requires a balance between system response speed and stability. When β is large (e.g., 0.3), the system responds quickly to changes in soil layers, but may produce large control fluctuations; when β is small (e.g., 0.1), the system has good stability, but its adaptability to rapidly changing soil layers decreases. Under typical soft soil construction conditions, a β value of 0.15~0.25 is recommended. At this value, the system can complete the weight switching of new soil layers within approximately 0.5 seconds while maintaining stability.

[0178] Step 703: Limit the rate of change of soil resistance characteristic indicators.

[0179] After calculating the rate of change dR / dt, the system must perform a limiting check. Because sensor noise can cause the instantaneously calculated derivative value to be abnormally large, directly using it for prediction would introduce overshoot. Therefore, the system defines a maximum permissible rate of change dR. max (For example, 20 / second). If the calculated |dR / dt|>dR max Then the rate of change is forcibly truncated to dR. max or -dR max The rate of change after amplitude limiting is used for subsequent prediction and compensation R. pred =R curr +(dR / dt) limited *τ d This can effectively prevent malfunctions in grouting volume caused by signal glitch.

[0180] Step 704, emergency handling when image tracking is lost.

[0181] When no valid identifying feature points can be detected in multiple consecutive frames (e.g., 10 frames) of visual images, the system determines that image tracking has been lost. At this time, the system switches to a pure inertial estimation mode, relying solely on IMU data for short-term state prediction, while simultaneously issuing an alarm to the operator. In pure inertial mode, the system locks the grouting flow rate to the calculated value at the most recent valid image observation time, and automatically reduces the stirring speed until image tracking is restored after the loss continues for more than a preset time threshold (e.g., 3 seconds).

[0182] Step 705, protection handling when tilt angle exceeds limit.

[0183] When the real-time tilt angle γ exceeds the maximum permissible tilt angle threshold γ maxWhen the tilt angle exceeds the limit, the system enters the tilt protection mode. In this mode, the system issues an alarm to the operator through the interface and a buzzer, prompting them to correct the attitude of the mixing head. At the same time, the system automatically increases the grouting flow rate to 1.2 to 1.5 times the normal value to compensate for the uneven mixing caused by excessive tilt. If the tilt angle continues to exceed the limit for more than a preset time (e.g., 5 seconds), the system will automatically pause the drilling operation until the operator confirms and restores the normal attitude.

[0184] Example 8, as Figure 4 As shown in the figure, this embodiment provides a three-dimensional positioning and trajectory tracking system for an excavator, which is installed on a tracked hydraulic excavator.

[0185] An excavator consists of basic structures such as the body, traveling mechanism, boom, arm, mixing arm, and mixing head.

[0186] An RTK+IMU positioning module is installed on the machine body. A forward-looking wide-angle camera is installed on the upper left front side of the cab via an extension bracket. Multiple sets of visual image recognition marking structures are arranged on the outer surfaces of the boom, forearm, mixing arm, and mixing head. A torque and speed sensing module is installed on the mixing shaft. A slurry flow controller is connected in series in the slurry delivery pipeline leading to the mixing head. A central control unit and a display and alarm terminal are arranged in the cab.

[0187] Each module communicates with the central control unit via a wired bus or wireless link, forming an overall data acquisition, computing, and control system.

[0188] Before construction begins, the operator imports the electronic plan and grid division information of the construction area into the central control unit via a display terminal, and inputs the design plan location, design depth, and design total grouting volume V for each mixing point. design Design stirring time t design Design stirring radius r design and the target drag coefficient R target Parameters such as these.

[0189] The central control unit calls upon the pre-calibrated internal and external parameters of the camera and the geometric positions of each marker structure in the fuselage coordinate system to establish the transformation relationship between the engineering coordinate system, the fuselage coordinate system, and the camera coordinate system. Simultaneously, the height offset of the mixing head's spray nozzle relative to the geometric center of the visual image marker is set to h. offset This lays a solid foundation for the subsequent three-dimensional positioning and attitude calculation of the stirring head.

[0190] like Figure 5 As shown, the field of view of the forward-looking wide-angle camera covers the movement areas of the upper arm, forearm, mixing arm, and mixing head.

[0191] Visual image markings on the boom, arm, and mixing arm are arranged along the length of the components. The markings on the mixing arm are positioned higher than the bottom of the mixing head, so that the markings remain visible on the ground and within the camera's field of view even after the mixing head has penetrated the soil.

[0192] The outer surface of the stirring head is marked with high-contrast markers to indicate the axial direction. A camera continuously captures video images and transmits them to the central control unit. The central control unit runs an image recognition algorithm to extract the pixel coordinates of each marker point from the image. Combining the camera's orientation and the marker's geometric position, it calculates the three-dimensional coordinates of the marker point in the machine's coordinate system and the unit direction vector k of the stirring head's axis. axis .

[0193] The fuselage RTK+IMU positioning module outputs the three-dimensional position (X) of the fuselage reference point in the engineering coordinate system in real time. RTK ,Y RTK Z RTK The central control unit constructs the global attitude rotation matrix R of the fuselage based on the attitude angle information. body_global The pose sensing system simultaneously calculates the relative position (x, y) of the mixing head's spray nozzle in the machine's coordinate system. v ,y v ,z v And define the relative position vector: P visual_body =[x v y v z v +h offset ].

[0194] Then, the central control unit calculates the three-dimensional coordinates (X, Y, Z) of the mixing head injection nozzle in the engineering coordinate system according to the following relationship:

[0195] X=X RTK +(R body global ·P visual body ) x ;

[0196] Y=Y RTK +(R body global ·P visual body ) y ;

[0197] Z=Z RTK +(R body global ·P visual body ) z .

[0198] in(·) x ,(·) y ,(·) zThese represent the components of the vector in the X, Y, and Z directions of the engineering coordinate system, respectively. By differentiating the Z values ​​at consecutive time points, the central control unit can obtain the vertical downward velocity v of the stirring head. down , as the drilling speed.

[0199] In terms of attitude calculation, the central control unit represents the attitude matrix of the stirring head in the global coordinate system as follows:

[0200] R visual_global =R body_global ·R visual_body ,

[0201] Where R visual_body Let e ​​be the attitude matrix of the stirring head relative to the fuselage. Let the vertically downward unit vector in the global coordinate system be e. z =[0 0 -1],

[0202] The real-time tilt angle (verticality deviation angle) γ of the stirring head is defined as:

[0203] γ=arccos(∣(R visual global ·e z )·k axis |);

[0204] Where e z =[0,0,-1] T Let k be a vertically downward unit vector. axis This is the unit vector pointing from the nozzle to the top of the mixing head. When γ = 0°, the axis of the mixing head coincides with the vertical direction. Through proper arrangement of markings and precise calibration, the tilt angle measurement accuracy in this embodiment can be better than 0.4°. Based on this, the central control unit combines the current mixing depth L... stir Calculate the equivalent disturbance radius of the stirring head in the horizontal plane for use in sweep volume calculation: This value represents the semi-major axis of the horizontal cross-section ellipse, characterizing the instantaneous perturbation range of the stirring head at the current depth. When hour, The increase is not significant; when hour, .

[0205] Calculate the horizontal offset of the mixing point, which is used for quality evaluation and overlap inspection of the mixing point:

[0206] At the bottom of the stirring point ( The maximum offset at () is .

[0207] Calculate the maximum horizontal coverage area of ​​the mixing points for checking the overlap between adjacent mixing points: the horizontal distance from the top of the mixing point to the farthest disturbed edge at the bottom of the mixing point.

[0208]

[0209] The determination of a decrease in the effective overlap width between adjacent mixing points should be based on The comparison is made with the spacing between stirring points, rather than a single point. .

[0210] As γ increases, the projection of the stirring section on the horizontal plane changes from an approximately circular shape to an elliptical shape, and the effective overlap width between adjacent stirring points decreases. Based on this, the system can subsequently determine areas where there may be insufficient overlap or missed stirring.

[0211] The torque and speed sensing module measures the average torque T and speed n of the mixing motor in real time and sends the data to the central control unit. The central control unit sets the soft soil calibration torque T0, rated idling speed n0, design lowering speed v0, and normalized maximum tilt angle γ. max Based on this, the soil resistance coefficient is constructed:

[0212] R=0.5×T / T0+0.3×n0 / n+0.2×v0 / ∣v down |+0.15×(γ / γ) max ) 2 .

[0213] The first three parameters comprehensively reflect the impact of torque, rotational speed, and lowering speed on soil hardness and mixing difficulty. The last parameter is the square of the inclination angle penalty, which amplifies the R value when the mixing head is tilted significantly, allowing this indicator to more accurately represent the construction difficulty. The central control unit can also classify the soil into extremely soft silt, normal soft soil, relatively hard clay, and stiff plastic clay grades based on the value of R, and mark them with different colors or labels on the interface.

[0214] Before calculating the soil resistance coefficient, the sensitivity of the current soil layer to each sensor parameter is identified based on the characteristics of sensor data changes within the sliding time window, and the weighting coefficients are dynamically adjusted to form an adaptive soil resistance coefficient R*.

[0215] In this embodiment, the sliding window length Δt window If the time interval is set to 200 milliseconds, the number of sampling points within the window is N = Δt / δt.

[0216] Let the torque, rotational speed, lowering speed, and tilt angle at the i-th sampling moment within the window be denoted as T(i), n(i), v(i), and γ(i), respectively, where i = 1, 2, ..., N.

[0217] The coefficient of variation is calculated for each sensor parameter sequence within the window to characterize the sensitivity of the parameter to fluctuations in the current soil layer.

[0218] The torque coefficient of variation is defined as: CVT =σ T / μ T ;

[0219] Where μ T The mean torque within the window, μ T =(1 / N)×Σ i=1 N T(i); σ T σ represents the standard deviation of torque within the window. T =sqrt((1 / N)×Σ i=1 N (T(i)-μ T ) 2 ).

[0220] Similarly, calculate the coefficient of variation of rotational speed CV. n Coefficient of variation (CV) of the decommissioning speed v Inclination coefficient of variation (CV) γ .

[0221] Define the original weights of each parameter as w T 0 =0.5、w n 0 =0.3、w v 0 =0.2、w γ 0 =0.15.

[0222] Calculate the weighting factor for the coefficient of variation:

[0223] f T =1+α×(CV T -CV ref );

[0224] Among them, CV ref The reference coefficient of variation (range 0.05~2.0, preferably 1.0); α is the sensitivity adjustment coefficient, which can be determined by those skilled in the art through experiments on the verification data based on the actual sensor noise level and the rate of soil layer change.

[0225] Similarly, calculate f n f v f γ .

[0226] Lower bound truncation of the weighting factors ensures f T f n f v f γ All values ​​should be no less than 0.5 to avoid excessive weight decay.

[0227] Calculate adaptive weights:

[0228] w T * =w T 0 ×f T / (w T 0 ×f T +w n 0 ×f n +w v 0 ×f v +w γ 0 ×f γ );

[0229] Similarly, calculate w n * w v * w γ * This normalizes the sum of the four weights to 1.

[0230] R * =w T * ×T / T0+w n * ×n0 / n+w v * ×v0 / ∣v down ∣+w γ * ×(γ / γ max ) 2 ;

[0231] The adaptive resistance coefficient R* automatically adjusts the contribution ratio of each component in different soil layers: when torque fluctuates, the torque weight is increased; when the ground velocity fluctuates, the velocity weight is increased, thus more accurately reflecting the actual construction difficulty of the current soil layer.

[0232] In terms of shotcrete control, the central control unit calculates the design shotcrete intensity based on the total design grouting volume and the design mixing time.

[0233] Q design =V design / t design (L / s).

[0234] Let the current stirring depth be z. cur The designed stirring depth is z design Inclination compensation coefficient k tilt Real-time shotcrete setting flow rate Q actual Given by the following formula:

[0235] Qactual =Q design ×(0.6+0.008R)×z cur / z design ×(1+k tilt sinγ).

[0236] Among them, (0.6+0.008R) makes the shotcrete volume vary with soil resistance. The coefficients 0.6 and 0.008 are empirical coefficients, which can be determined by those skilled in the art through regression analysis based on historical construction data of similar projects, or can be adaptively adjusted according to actual engineering conditions. When R is within the normal range, the corresponding shotcrete volume is approximately 0.6 to 1.4 times the design value; (1+k tilt sinγ) is used to compensate for the increased stirring path and increased disturbance volume caused by the tilt of the stirring head, k tilt This is the tilt angle compensation coefficient, and its specific value can be adaptively adjusted according to the actual geometry of the mixing head and the construction conditions. The shotcrete flow controller is based on Q. actual By adjusting the valve opening or pumping frequency, its built-in flow meter feeds back the actual grouting flow rate and cumulative grouting volume to the central control unit in real time, forming a closed-loop control of the grouting process.

[0237] The spatial trajectory of the mixing head and the movement of the excavator are controlled by the operator; this system only automatically adjusts the amount of shotcrete.

[0238] When calculating the target shotcrete flow rate, not only are the soil resistance coefficient R* and inclination angle γ at the current moment considered, but the rate of change of resistance coefficient is also introduced for prediction compensation, so as to achieve early response to sudden changes in soil layers.

[0239] Based on the resistance coefficient sequence R*(i) (i=1,2,...,N) within the sliding window, linear regression is used to calculate the rate of change of the resistance coefficient over time:

[0240] dR / dt=(N×Σ i=1 N (i×R * (i))-Σ i=1 N i×Σ i=1 N R * (i)) / (N×Σ i=1 N i 2 -(Σ i=1 N i) 2 ) / δt;

[0241] Where δt is the sampling period.

[0242] Let the system execution delay be τ. d(Including valve response time and pipeline pressure build-up time, ranging from 0.3 to 1.0 seconds, preferably 0.5 seconds), the estimated resistance coefficient at the predicted time is:

[0243] R pred * =R * +(dR / dt)×τ d ;

[0244] For R pred * Cut off the upper and lower limits to ensure that it is within a reasonable range (0~120).

[0245] Q actual * =Q design ×(0.6+0.008×R * pred )×z cur / z design ×(1+k tilt ×sinγ)

[0246] When dR / dt>0 (resistance is increasing, indicating an approach to a harder soil layer), Q actual * Increase in advance to compensate for execution delays; when dR / dt < 0 (resistance decreasing, about to enter a softer soil layer), Q actual * Reduce the volume in advance to avoid excessive slurry.

[0247] To prevent misjudgments caused by sensor noise, a limit is set on the rate of change of the drag coefficient:

[0248] ∣dR / dt∣ limited =min(∣dR / dt∣,dR max );

[0249] Where dR max The maximum allowable rate of change is set (range 10~30 / second, preferably 20 / second). When this limit is exceeded, the prediction compensation is calculated based on the limit.

[0250] The display terminal in the driver's cab serves as a human-machine interface and alarm unit, used to display the construction status and system calculation results.

[0251] The interface displays the construction site grid in the form of an electronic plan. Green grids indicate areas that have been mixed, while red grids indicate areas that have not yet been mixed. The current position of the mixing head is marked with a prominent icon.

[0252] The interface simultaneously displays the real-time depth curve and 3D trajectory projection of the stirring head, as well as the numerical value and bar indicator of the tilt angle γ. When γ exceeds the preset limit, the interface will issue an alarm by changing color or flashing.

[0253] Torque, resistance coefficient R, target grouting flow rate and actual grouting flow rate are compared using curves or bar charts. The completion rate of grouting volume per unit volume and the current mixing quality score are also continuously updated on the interface.

[0254] Operators can adjust the excavator's travel path and mixing rhythm based on this information to avoid under-mixing or over-mixing.

[0255] like Figure 6 As shown, the central control unit performs real-time control as well as construction quality evaluation and data recording.

[0256] After each mixing point is completed, the system calculates the time-averaged resistance coefficient R at that mixing point. avg Actual mixing time t stir Root mean square fluctuation of stirring speed Δn rms Actual cumulative grouting volume V actual And the maximum tilt angle γ during the mixing construction process. max_in_pile And calculate the mixing quality score based on this:

[0257] .

[0258] When γ > 8° or the maximum tilt angle of the mixing arm γ is at any moment during construction, max_in_pile When the angle is greater than 6°, the system determines that the mixing point is a serious tilt defect, forcibly limits the score to 60 points or below, and marks it in red on the construction plan, indicating that a re-inspection or rework is required.

[0259] The central control unit records data such as mixing location, start and end time, spatial coordinates of the mixing head, depth, tilt angle, spraying volume, and quality score according to the mixing point coordinates, forming a complete construction log. The construction log can be stored locally or uploaded to the on-site office computer or remote server via the network, allowing project managers to monitor construction progress and quality in real time, and providing a basis for later acceptance and dispute resolution.

[0260] This embodiment illustrates how to achieve three-dimensional positioning, trajectory tracking, and intelligent control of shotcrete in an excavator mixing head with a single camera and a chassis RTK+IMU configuration.

[0261] Those skilled in the art can adjust the camera model, identification structure, weighting coefficient, or threshold parameters according to engineering needs. Such equivalent substitutions or improvements do not change the basic principles of the present invention and should be considered to fall within the protection scope of the present invention.

[0262] To address the blind spots and positioning errors in underground pose perception, a vision-inertial fusion algorithm for engineering machinery was adopted, combined with the placement of high-level markers on the mixing arm. Utilizing multi-body kinematics chain constraints and the short-term prediction capabilities of the IMU, not only was high-frequency noise caused by machine vibration compensated, but the underground trajectory could still be calculated from the visible markers above the mixing head even when it was submerged underground. This achieved continuous, high-precision 3D positioning throughout the entire process, solving the problems of large indirect positioning errors and the inability to track blind spots in pose perception.

[0263] To address the issues of uneven grouting and control model distortion caused by tilting, a physical control model based on effective swept volume and energy density is introduced. By discretizing the trajectory and introducing a nonlinear correction for the projected area based on the tilt angle, the system can calculate the actual soil volume disturbed by the mixing head in real time and dynamically adjust the flow rate based on the target grouting density per unit volume. This closed-loop control based on physical mechanisms replaces traditional linear empirical formulas, ensuring precise matching between the injected solidifying agent and the actual amount of soil cut under tilting or variable-speed movement of the mixing head. This solves the problem of localized under-mixing or grout waste caused by the failure to calculate changes in swept volume.

Claims

1. A method for dynamic sensing and intelligent precise control of multimodal data in the sludge solidification construction process, characterized in that, include: Real-time acquisition of global pose data of the excavator body, visual image sequence including the marking structure arranged on the mixing arm, and driving power parameters of the mixing drive shaft; wherein, the marking structure is arranged along the length of the mixing arm and the arrangement position is higher than the bottom of the mixing head by a set distance, so that the marking structure is within the visual range when the mixing head is submerged underground for operation; Based on global pose data and visual image sequences, the real-time three-dimensional motion trajectory and attitude of the stirring head are calculated using visual-inertial fusion. Based on the driving force parameters and the vertical motion velocity derived from the real-time three-dimensional motion trajectory, the soil resistance characteristic index that characterizes the current soil properties is calculated. Based on the effective swept volume generated by the real-time three-dimensional motion trajectory, combined with the soil resistance characteristic index, the target grouting control volume at the current moment is calculated. Adjust the execution parameters of the grouting equipment according to the target grouting control volume; Based on the effective swept volume generated from the real-time three-dimensional motion trajectory, and combined with soil resistance characteristic indicators, the target grouting control volume at the current moment is calculated, including: Discretize the real-time three-dimensional motion trajectory in the time dimension to obtain the micro-element trajectory segment at the current moment; Based on the geometric properties of the stirring head and the spatial pose of the infinitesimal trajectory segment, a three-dimensional spatial geometry swept by the stirring head at the current moment is constructed, and the volume of the three-dimensional spatial geometry is calculated as the effective sweep volume; wherein, the geometric properties include at least the design stirring radius of the stirring head; Based on the soil resistance characteristics, determine the target grouting density per unit volume required for the current soil. Based on the effective swept volume and the target unit volume grouting density, calculate the target grouting control amount so that the ratio of the target grouting control amount to the effective swept volume approximates the target unit volume grouting density. Based on the geometric properties of the stirring head and the spatial pose of the infinitesimal trajectory segments, a three-dimensional spatial geometry swept by the stirring head at the current moment is constructed, and the volume of the three-dimensional spatial geometry is calculated as the effective sweep volume, including: The real-time tilt angle of the stirring head is analyzed from the spatial pose of the infinitesimal trajectory segment; A virtual cylinder is constructed based on the designed stirring radius of the stirring head, and the virtual cylinder is modeled by spatial sweeping along the micro-element trajectory segment to obtain a three-dimensional spatial geometry; The projected area of ​​a 3D spatial geometry on a horizontal plane is corrected by using the real-time tilt angle, and the corrected effective sweep volume is calculated accordingly.

2. The method according to claim 1, characterized in that, The real-time 3D motion trajectory and attitude of the stirring head are calculated using vision-inertial fusion, including: Construct a multibody kinematics chain model of an excavator that includes the fuselage, boom, arm and mixing arm; Identify the currently visible signage structure from a sequence of visual images and calculate the local observation pose of the signage structure relative to the forward-looking wide-angle camera based on perspective geometry. Based on the multibody kinematics chain model, the kinematics transfer and coordinate transformation of the local observed pose and the global pose data are performed to obtain the initial state value of the stirring head in the engineering coordinate system. By utilizing the motion continuity constraint in the time dimension to smooth and optimize the initial state value, the real-time three-dimensional motion trajectory and attitude of the stirring head are obtained.

3. The method according to claim 2, characterized in that, The initial state is smoothly optimized by utilizing the motion continuity constraint in the time dimension, including: Obtain the estimated state of the stirring head at the previous moment and extract the inertial measurement data contained in the global pose data; Based on inertial measurement data and a multibody kinematics chain model, the estimated state of the stirring head at the previous moment is differentially extrapolated to generate the prior predicted state at the current moment. An observation residual equation is constructed, which characterizes the reprojection error between the feature points projected onto the image plane from the prior predicted state and the actual identified feature points in the visual image sequence. By minimizing the observation residual equation to correct the prior predicted state, the real-time three-dimensional motion trajectory and attitude at the current moment are obtained.

4. The method according to claim 1, characterized in that, Based on the driving dynamic parameters and the vertical motion velocity derived from the real-time three-dimensional motion trajectory, soil resistance characteristic indicators characterizing the current soil properties are calculated, including: The stirring torque and stirring speed are extracted from the driving power parameters, and the real-time tilt angle of the stirring head relative to the vertical direction is analyzed from the attitude data of the real-time three-dimensional motion trajectory. The stirring torque, stirring speed, and vertical motion speed are normalized to construct torque components, speed components, and velocity components, respectively. A tilt penalty term is constructed based on the real-time tilt angle, and the tilt penalty term is configured to increase non-linearly with the increase of the real-time tilt angle; The soil resistance characteristic index is obtained by weighted summation of the torque component, rotational speed component, velocity component, and tilt angle penalty term.

5. The method according to claim 4, characterized in that, Before performing a weighted summation of the torque component, speed component, velocity component, and tilt angle penalty term, the following is also included: Establish a sliding time window that includes historical data of a set length prior to the current moment; Calculate the time series variation coefficients of stirring torque, stirring speed, and vertical motion speed within the sliding time window; Based on the time series coefficient of variation, the sensitivity of the current soil layer to various parameters is identified, and the weight coefficients used for weighted summation are dynamically adjusted accordingly, so that parameters with larger coefficients of variation within the sliding time window have larger weight coefficients.

6. The method according to claim 1, characterized in that, Before adjusting the execution parameters of the grouting equipment according to the target grouting control volume, the following steps are also included: Based on the soil resistance characteristic index calculated at multiple consecutive time points, the rate of change of the soil resistance characteristic index with time is calculated. Based on the rate of change and the preset system response delay time of the grouting equipment, the predicted soil resistance characteristic index is predicted for the next execution time after the current time. The target grouting control volume is corrected by using the estimated soil resistance characteristic index to obtain the compensated execution command, and then the compensated execution command is sent to the grouting equipment.

7. The method according to claim 1, characterized in that, Based on the current moment's infinitesimal trajectory segment, after completing the mixing construction at the corresponding mixing point, the process also includes: The actual unit volume grouting density corresponding to each micro-element trajectory segment during the entire construction process is statistically analyzed. The actual unit volume grouting density is calculated by the ratio of the actual grouting volume to the effective sweeping volume. Calculate the dispersion of the actual unit volume grouting density relative to the target unit volume grouting density; The evaluation result of the mixing construction quality at this mixing point is generated based on the distribution dispersion, reflecting the uniformity of the internal curing agent distribution.

8. The method according to claim 5, characterized in that, Before identifying the sensitivity of the current soil layer to various parameters, the following steps are also included: A soil layer feature model library is pre-built, which contains a variety of typical soil types; For each typical soil type, its corresponding torque sensitivity, rotational speed sensitivity, and velocity sensitivity are pre-calibrated, and a preset benchmark combination of weighting coefficients for that soil type is established accordingly. When dynamically adjusting the weighting coefficients used for weighted summation, the current soil layer is matched to the typical soil type in the soil feature model library based on the current time series variation coefficient, and the preset benchmark combination corresponding to the typical soil type is called as the adjustment benchmark.

9. The method according to claim 1, characterized in that, The method also includes closed-loop control of stirring energy: Based on the stirring torque, stirring speed, and infinitesimal time interval, the energy input value of the stirring head within the current infinitesimal trajectory segment is calculated; The real-time energy density per unit volume is calculated based on the ratio of energy input value to effective swept volume. The real-time energy density per unit volume is compared with the preset target energy density, and the speed or torque output of the stirring drive shaft is adjusted according to the comparison result.

10. The method according to claim 3, characterized in that, Before correcting the prior predicted state by minimizing the observation residual equation, an outlier removal step is also included: Calculate the single-point reprojection error of each identified feature point in the visual image sequence; The single-point reprojection error is compared with a preset confidence threshold. When the single-point reprojection error of a certain identifier feature point exceeds the confidence threshold, the identifier feature point is determined to be an outlier point that is contaminated or occluded. The outlier point is then removed when constructing the observation residual equation, and the remaining valid feature points are used to minimize the solution.

11. The method according to claim 4, characterized in that, The specific calculation method for the tilt penalty is as follows: Obtain the preset maximum allowable tilt angle threshold; Calculate the ratio of the real-time tilt angle to the maximum permissible tilt angle threshold; By using the square or higher power of the ratio as the tilt angle penalty term, the soil resistance characteristic index increases non-linearly when the real-time tilt angle approaches the maximum allowable tilt angle threshold.

Images