Loading operation motion trail optimization method for skid steer loader

Abstract

The invention relates to the technical field of engineering machinery control, in particular to a skid loader loading operation motion trail optimization method. The method comprises the following steps: carrying out loading attribute feature identification processing of environment perception based on loading operation demand data and multi-source environment perception data, and generating loading attribute environment perception monitoring data; acquiring design parameters of the skid steer loader; based on the design parameters of the skid steer loader, optimizing the multi-body dynamic model of the loader; establishing a panoramic sensing coordinate system of the loading attribute based on the loading attribute environment sensing monitoring data to obtain a loading attribute panoramic sensing coordinate system; performing global path analysis of slip loading according to the loading attribute panoramic sensing coordinate system, and generating slip loading global path data; and carrying out loading motion trail optimization processing on the slip loading global path data to obtain slip loading global optimization motion trail data. According to the invention, intelligent optimization control of the loading operation movement track of the skid steer loader is realized.

Description

Technical Field

[0001] This invention relates to the field of engineering machinery control technology, and in particular to a method for optimizing the motion trajectory of a skid steer loader during loading operations. Background Technology

[0002] Skid steer loaders, as flexible and efficient construction machinery, are widely used in construction, mining, and agricultural production. The efficiency and precision of their loading operations directly impact the overall project progress and operating costs. Skid steer loader operations largely rely on operator experience and judgment, manually controlling actions such as shoveling, lifting, and unloading. This operating mode is not only labor-intensive but also significantly affected by human factors, easily leading to low loading efficiency, increased energy consumption, and even accelerated mechanical wear and increased operational safety risks due to insufficient operational precision. With the deepening application of intelligent technology in the construction machinery field, higher demands are placed on the automation and precision of skid steer loaders, especially in complex operating environments. Optimizing the loading trajectory has become a key technical bottleneck for improving equipment performance. However, existing methods for optimizing the loading trajectory of skid steer loaders are designed for single working conditions or static environments and lack adaptability to dynamic operating scenarios. They are based on simple planning of preset paths and do not fully consider the real-time changes in the working environment and the dynamic characteristics of the loaded materials. This results in a large deviation between the optimized trajectory and the actual operating requirements. In the modeling process, the multibody dynamics of the loader are often simplified and the mechanical response during the contact process between the bucket and the material is ignored, making it difficult to guarantee the accuracy and feasibility of trajectory optimization. Summary of the Invention

[0003] Based on this, the present invention provides a method for optimizing the motion trajectory of a skid steer loader during loading operations, in order to solve at least one of the above-mentioned technical problems.

[0004] To achieve the above objectives, a method for optimizing the loading trajectory of a skid steer loader includes the following steps: Step S1: Obtain loading operation demand data; use the multi-source image monitoring equipment built into the skid steer loader to perform multi-source environmental perception monitoring processing on the loading operation area to obtain multi-source environmental perception data; based on the loading operation demand data and the multi-source environmental perception data, perform loading attribute feature identification processing for environmental perception to generate loading attribute environmental perception monitoring data. Step S2: Obtain the design parameters of the skid steer loader; based on the design parameters of the skid steer loader, perform multibody dynamics modeling and operational characteristic parameter identification and optimization processing of the loader to obtain an optimized multibody dynamics model of the loader; Step S3: Establish a panoramic perception coordinate system for loading attributes based on the environmental perception monitoring data of loading attributes to obtain the panoramic perception coordinate system for loading attributes; perform global path analysis of sliding loading based on the panoramic perception coordinate system for loading attributes to generate global path data for sliding loading. Step S4: Extract target load images at each local loading moment based on loading attribute environmental perception monitoring data and skid loader global path data to obtain local target load image data; perform bucket loading contact optimization trajectory analysis based on local target load image data to obtain group local optimized bucket contact trajectory data; optimize the loading motion trajectory of the skid loader global path data based on the group local optimized bucket contact trajectory data to obtain skid loader global optimized motion trajectory data; execute skid loader loading intelligent control operation based on skid loader global optimized motion trajectory data.

[0005] The beneficial effects of this application are as follows: By acquiring loading operation requirements and using multi-source image monitoring equipment on the airborne platform to collect data on the operation area, perform extrinsic parameter correction and temporal synchronization, and identify loading attribute features based on binarization / segmentation, the present invention can significantly improve the perception accuracy and spatiotemporal consistency of the operation scene. Multi-source fusion and geometric extrinsic parameter / temporal correction reduce projection errors and time deviations between different sensors, making depth information, image contours, and point cloud data comparable in the same panoramic coordinate system, thereby improving the recognition rate and confidence of loading attributes such as material boundaries, stack height, slope, and retrievable locations. Accurate loading attribute environmental perception and monitoring data not only provide reliable input for subsequent path feasibility domain analysis and local contact decision-making, but also reduce ineffective shoveling, repetitive operations, or safety risks caused by perception misjudgments, thereby improving operational efficiency, reducing energy consumption, and improving operational stability and safety. Based on the design parameters of the skid steer loader, component connection relationship analysis, actuator capability profile, lifting characterization, and multibody dynamics modeling are performed. Combined with adaptive parameter identification and order reduction processing based on operational characteristics, an optimized multibody dynamics model that reflects the nonlinear characteristics of the machine body and is suitable for online calculation is obtained. This model, by explicitly considering linkage kinematic constraints, actuator capability boundaries, and lifting height-capacity mapping, ensures that path planning and contact trajectory generation are strictly feasible in terms of dynamics, avoiding commands exceeding hydraulic or structural capabilities. Through lifting dynamic order reduction and linearization approximation, both model accuracy and real-time performance are balanced, facilitating online identification and controller integration. The adaptively identified and optimized model can perform parameter corrections as actuators age, load changes, and site conditions change, improving trajectory prediction accuracy, energy consumption estimation reliability, and lifting stability assessment capabilities, contributing to the generation of safer, more energy-efficient, and more effective loading trajectories. A panoramic perception coordinate system is established based on loading attribute environmental perception monitoring data, and skid path feasibility domain analysis, skid characteristic line supplementation, and cost feature evaluation are conducted on it, ensuring that the generated global path is highly consistent in spatial representation and dynamic constraints. The panoramic coordinate system unifies multi-source sensing data into a single reference frame, reducing projection and position errors. Feasibility domain analysis explicitly considers the differential / skid steering characteristics, turning radius, ground friction, and slope limitations of the skid steer loader, thus eliminating dynamically infeasible or high-risk path candidates. Line completion and skid characteristic fitting smooth the path and match the machine's actual response, reducing path tracking errors and control difficulties caused by lateral skid. Furthermore, by constructing multi-dimensional cost features including energy consumption, slope, short-term obstacle probability, lifting stability, and operational efficiency, and adaptively adjusting weights based on operational requirements, the resulting global path balances energy minimization, operational efficiency, and safety, improving the decision-making quality and applicability of global planning. This provides robust and interpretable initial trajectory data for subsequent local optimization and closed-loop execution.Based on the global path, image extraction of the target load at each local loading moment and group correlation and force-potential energy field analysis of local bucket contact can significantly improve the success rate of single scooping and overall loading efficiency. Local image extraction and shape feature analysis make the evaluation of candidate contact points and maximum loading depth more accurate, while force-potential energy field analysis provides a physically quantifiable criterion for the superiority or inferiority of contact strategies, thus prioritizing the effective loading contact trajectory while ensuring stability. Group prediction of multiple local contact trajectories and the establishment of the state transition and observation matrix of the load can predict and filter the short-term evolution of material morphology, reducing erroneous decisions caused by sensor noise or single observation errors. Mapping the group-optimized local contact trajectory back to the global path and performing local corrections makes the final output globally optimized motion trajectory more closely match the actual material changes and actuator capabilities while maintaining dynamic feasibility, reducing heavy scooping, empty loading, and energy waste. Combined with intelligent control execution, a closed-loop feedback of perception-planning-execution is realized, enhancing the response capability to dynamic obstacles and sudden scenarios, thereby improving operational robustness, safety, and overall operational efficiency.

[0006] Therefore, the skid steer loader loading trajectory optimization method of this application can significantly improve the accuracy and feasibility of trajectory optimization. By introducing a multi-source environmental perception and temporal calibration mechanism, real-time dynamic monitoring and fusion of the working area environment and material status are achieved, ensuring high precision and spatiotemporal consistency of the perceived data. This allows trajectory planning to be adjusted in real time according to changes in the environment and materials, avoiding deviations caused by traditional static path planning. Simultaneously, by combining multibody dynamics modeling based on design parameters and adaptive parameter identification, the nonlinear kinematics and dynamic characteristics of the loader are fully preserved, ensuring that trajectory generation strictly meets actuator capabilities and structural constraints, guaranteeing the trajectory's physical feasibility. By introducing mechanical response analysis and potential energy field modeling of the bucket-material contact process, the force characteristics of the scooping action and material deformation characteristics are explicitly considered in trajectory optimization, effectively improving the stability of the scooping action and loading efficiency. Attached Figure Description

[0007] Figure 1 This is a flowchart illustrating the steps of a method for optimizing the motion trajectory of a skid steer loader during loading operations, as described in this invention. Figure 2 for Figure 1 A detailed flowchart illustrating the implementation steps of step S3. Figure 3 for Figure 1 A detailed flowchart illustrating the implementation steps of step S4. The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0008] The technical method of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0009] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. Functional entities may be implemented in software, or in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods. The term "and / or" as used herein includes any and all combinations of one or more of the associated items listed.

[0010] To achieve the above objectives, please refer to Figures 1 to 3 This invention provides a method for optimizing the loading trajectory of a skid steer loader. In an embodiment of this invention, please refer to... Figure 1 The diagram shown is a flowchart illustrating the steps of a method for optimizing the loading trajectory of a skid steer loader according to the present invention. The method includes the following steps: Step S1: Obtain loading operation demand data; use the multi-source image monitoring equipment built into the skid steer loader to perform multi-source environmental perception monitoring processing on the loading operation area to obtain multi-source environmental perception data; based on the loading operation demand data and the multi-source environmental perception data, perform loading attribute feature identification processing for environmental perception to generate loading attribute environmental perception monitoring data. In this embodiment of the invention, the loading operation requirement information is read by the task description unit and stored in structured fields. These fields include material type (bulk or block), target loading quantity, target unloading location coordinates, operation time window, maximum allowable roll angle, and target weight. The airborne sensing hardware consists of a monocular or binocular camera array, a temporal depth camera, a short-range lidar, and an inertial measurement unit. All sensors are sampled and synchronized using hardware triggering or a precise time protocol. Sampled data is buffered according to a timestamp sequence and used for subsequent registration. Geometric extrinsic parameter calibration uses a calibration board and calibration target for static calibration, combined with hand-eye calibration methods to obtain the transformation matrix from each sensor to the aircraft reference. Temporal calibration uses hardware trigger signal comparison or timestamp interpolation to correct sampling delay. Based on the extrinsic parameters, the laser point cloud is projected into the camera pixel coordinates. Depth and color information are fused through voxel grid downsampling and truncated band sign distance field fusion to generate a color voxel model and elevation raster. Image preprocessing first performs radial distortion correction on the projected image, then performs joint threshold binarization of grayscale and depth, uses morphological opening and closing operations to remove noise, and extracts contours using edge detection and connected component analysis, determining the stockpile boundary using convex hull and minimum bounding rectangle. The stockpile height is calculated based on the voxel height field, and the slope and aspect are estimated based on point cloud normals and RANSAC fitting of the normal surface. A regression mapping is established using visual texture roughness and historical friction observations to estimate the friction coefficient. A reachability test is performed in the panoramic coordinate system using a robotic arm / bucket attitude reachability model to generate a détente polygon. Finally, the location information, stockpile height, slope, détente polygon, friction estimation, and various confidence fields are structured to generate loading attribute environmental perception monitoring data.

[0011] Step S2: Obtain the design parameters of the skid steer loader; based on the design parameters of the skid steer loader, perform multibody dynamics modeling and operational characteristic parameter identification and optimization processing of the loader to obtain an optimized multibody dynamics model of the loader; In this embodiment of the invention, the design parameters of the skid steer loader include link length, three-dimensional coordinates of each hinge point, description of the bucket's geometric surface, hydraulic cylinder diameter and stroke, maximum valve flow rate and rated pressure, vehicle center of gravity and mass distribution, wheelbase and axle load distribution, installation offset and assembly tolerance. A topology diagram of the loading components is constructed based on the geometric parameters, treating components as nodes and hinges and slide rails as connecting edges. A homogeneous transformation matrix is ​​used to represent the relative pose of each joint, and closed-chain constraint equations are derived to generate linkage kinematic constraint data. The constraint Jacobian matrix is ​​used for numerical solution and handling of redundant degrees of freedom. The actuator analysis establishes a force-pressure mapping and a flow-velocity mapping based on hydraulic principles, and transforms the actuator capability profile into a force-velocity grid representation. Each point in the grid represents continuous output capability, response delay, and saturation boundary. Simultaneously, a nonlinear boundary model is established for dead zone and hysteresis behavior to form actuator capability data. The lifting characterization analysis calculates the lever arm and torque balance through discrete lifting angle step calculations. At each angle point, the maximum load capacity that can be supported is derived using static equilibrium, and the lifting stability limit is calculated by combining the centroid projection and support polygon criterion, outputting a lifting height-capacity sample curve. Multibody dynamics modeling establishes the Lagrange equation using generalized coordinate vectors, deriving the mass matrix, Coriolis and centrifugal terms, gravity term and contact Jacobian. The first-order inertia of the actuator, valve flow dynamics, and hydraulic saturation model are incorporated into the dynamic equations. The Lagrange multiplier method or penalty function method is used to handle closed-chain and contact constraints, establishing a multibody dynamics framework model with nonlinear boundary constraints. The adaptive parameter identification for operational characteristics employs a phased process: offline, joint angles, speeds, hydraulic pressures, and flow rates are recorded using a preset excitation sequence; inertial parameters and friction terms are estimated using least squares or weighted least squares; online, recursive least squares with a forgetting factor or extended Kalman filtering is applied to continuously correct key parameters; residual monitoring and convergence criteria are used to trigger parameter updates; and the final output is an optimized multibody dynamics model of the loader that meets real-time calculation requirements and reflects mechanical nonlinearity. To reduce the online computational load, the lifting dynamics are reduced in order, and the reduced-order results are mapped back to the architecture model to maintain the accuracy of the lifting characteristic representation and computational efficiency.

[0012] Step S3: Establish a panoramic perception coordinate system for loading attributes based on the environmental perception monitoring data of loading attributes to obtain the panoramic perception coordinate system for loading attributes; perform global path analysis of sliding loading based on the panoramic perception coordinate system for loading attributes to generate global path data for sliding loading. In this embodiment of the invention, a three-dimensional panoramic perception coordinate system is constructed based on the environmental perception monitoring data of the loading attributes. This coordinate system takes the loader body reference point as the origin and defines the forward axis, lateral axis, and vertical axis. Through extrinsic parameter transformation, all sensor data are uniformly mapped to the body coordinates. Elevation grids and occupancy grids are generated on the panoramic coordinate system. The terrain reference surface is extracted using ground plane fitting and the RANSAC method, and the occupancy representation and elevation map are constructed using the voxelization method. The feasible region analysis of the slip path is performed based on the occupancy and elevation maps: the loader's outline is projected and a buffer radius convolution is performed on the outline to obtain the mechanical space occupancy area. The slope threshold, lateral friction constraint, center of mass stability margin during lifting, and drive torque requirement are evaluated grid by grid. Grids that meet all constraints are marked as passable grids. For key intermediate nodes in the feasible region, skeleton extraction or key point extraction is performed to form a node graph. Spline curves or Bézier curves are used to continuously complete the paths between nodes. Then, based on the optimized multibody dynamics model of the loader, the lateral slip angle, lateral force, and wheel-ground interaction force are estimated on the curve, and penalty terms are established for curvature, acceleration, and slip response. The curve parameters are adjusted through constraint optimization to obtain a slip-optimized feasible path that conforms to the slip response. Path cost feature analysis calculates the value terms of each generation using the dynamic model: path length is obtained by integrating the arc length; energy consumption is estimated by integrating pressure multiplied by flow rate over time; slope cost is mapped as an increment based on the additional energy consumption and stability risk caused by the slope; short-term dynamic obstacle probability penalty is obtained by integrating the obstacle probability heatmap generated by the short-term target tracker on the path; dynamic environmental risk penalty is obtained by comparing the overturning matrix and recovery matrix to obtain the overturning risk score; lifting stability cost is calculated based on the stability margin calculated based on the distance between the center of mass and the support polygon during lifting; short-term operational efficiency cost is estimated based on the expected number of shoveling cycles and time. A linear cost mapping is constructed based on cost characteristics, and the weights of each cost item are adjusted according to the loading operation demand weight rules. A two-stage global path search strategy is adopted: at the grid level, a heuristic path search algorithm is used to quickly screen candidate paths, and then at the continuous level, nonlinear programming with dynamic constraints is used to constrain and optimize the candidate paths, outputting global path data for skid loading that satisfies both dynamic feasibility and demand trade-off.

[0013] Step S4: Extract target load images at each local loading moment based on loading attribute environmental perception monitoring data and skid loader global path data to obtain local target load image data; perform bucket loading contact optimization trajectory analysis based on local target load image data to obtain group local optimized bucket contact trajectory data; optimize the loading motion trajectory of the skid loader global path data based on the group local optimized bucket contact trajectory data to obtain skid loader global optimized motion trajectory data; execute skid loader loading intelligent control operation based on skid loader global optimized motion trajectory data.

[0014] In this embodiment of the invention, an extraction window is constructed on the global path of sliding loading, using each local loading moment as an anchor point. A spatial constraint box is generated based on the bucket's arrival posture and the feeding direction. The point cloud and image within the constraint box are cropped for local high-precision segmentation. The segmentation uses depth-guided connected component clustering combined with a surface segmentation algorithm to extract the target loading surface point set. Subsequently, convex hull calculation and principal component analysis are performed to obtain the surface principal direction, and the shape feature vector is calculated through local curvature and normal distribution. Bucket contact point difference analysis is performed by generating a set of candidate contact nodes on the target surface and using point-by-point collision detection with the bucket geometry and the loader's inverse kinematics to evaluate the achievable maximum loading depth of each contact point. The collision detection simulates the bucket feeding trajectory at discrete depth steps and checks joint limits, machine interference, and ground contact constraints. The maximum depth value that satisfies all constraints is recorded, and this depth, along with indicators such as the incident angle and contact normal deviation, forms the contact point difference feature. Based on the differences in contact points and the material mechanics model, a potential energy field is established. The mechanics model uses friction-cohesion and cutting resistance terms. The resistance is integrated along the proposed contact trajectory to estimate the required work, and the expected loading capacity is calculated based on the bucket's geometric volume. The potential energy ratio of the amount of material removed per unit of work is constructed as an evaluation index to rank candidate trajectories. Group local time correlation is achieved by treating several consecutive scoops as a state sequence. The state vector is defined to include the remaining volume, surface height distribution, and centroid position. Based on the contact trajectory prediction, a state transition matrix and observation matrix are constructed. Kalman filtering or particle filtering is used to fuse real-time observations and predictions to obtain a robust material state estimate, reducing the impact of sensor noise on single decisions. The contact trajectories after group correlation are optimized using multi-objective constraints to find the local optimal trajectory. The objectives are to minimize the unit work, maximize the stability margin, and satisfy the hydraulic and structural moment boundaries. The solver uses sequential quadratic programming and starts with the previous solution as a warm start. The solution results are inserted as local trajectory segments into the global path. After insertion, continuity and dynamic feasibility are verified. If necessary, the neighborhood path is reprogrammed locally to maintain overall feasibility. During the execution phase, the trajectory generator outputs attitude references and valve flow references. The position loop and force loop are controlled in parallel to achieve closed-loop tracking. The actuator feedback is monitored in real time, and when the deviation exceeds the threshold, it triggers a second local optimization or replanning.

[0015] Furthermore, step S1 includes the following steps: Step S11: Obtain loading operation requirement data; In this embodiment of the invention, the loading operation requirement data is generated into a standardized work order by the operation management unit before the task is issued. The work order fields adopt fixed enumeration and numerical range constraints, including information such as material category code, target loading mass, expected filling rate per bucket, unloading posture, operation time window, and safety thresholds (maximum tilt angle, minimum safe distance from obstacles, maximum slope). After the work order enters the parsing process, the dimensions are first unified, and the coordinate system is mapped from the construction coordinate system to the machine base coordinate system through a fixed transformation matrix. Then, consistency verification is carried out. If there is a conflict between the loading mass and the upper limit of the unloading point, the target quantity is compressed according to the safety threshold rules. When the time window and the path passage period do not overlap, the time sequence reordering strategy is triggered. After the verification is completed, an "operation requirement vector" is constructed in vector form. The vector contains weights, thresholds, geometric regions, and target indicators, and is bound to the unique task identifier, serving as a unified input reference for subsequent perception, modeling, and planning stages.

[0016] Step S12: Use the multi-source image monitoring equipment built into the skid steer loader to perform multi-source environmental perception monitoring and processing on the loading operation area, generate multi-source environmental perception data, and perform geometric external parameter correction and time sequence synchronization processing on the multi-source environmental perception data through the internal and external participation calibration records of the multi-source image monitoring equipment to generate fused environmental perception monitoring data. In this embodiment of the invention, the multi-source image monitoring device consists of a forward-facing binocular camera, a top-mounted structured light depth camera, a short-range lidar, and an IMU. The sampling timing is synchronized via hardware pulse triggering and a precise time protocol. Intrinsic parameter calibration is performed on a stable support, using a calibration board to obtain distortion and focal length. Extrinsic parameter calibration employs a hand-eye calibration process, aligning the coordinate systems of each sensor to the body's base coordinates. The calibration results generate a rigid body transformation matrix and record the temperature drift compensation coefficient. After sampling, the image undergoes radial and tangential distortion correction, depth is filled with holes and filtered to preserve edges, and the lidar point cloud is subjected to statistical outlier removal and voxel downsampling. The IMU attitude is constructed using median integral fusion to create a high-frequency attitude trajectory for motion compensation, and different sensors are interpolated and aligned on a common time grid. Subsequently, spatial registration is achieved under the action of extrinsic parameter transformation: the lidar point cloud is projected onto the image plane to form a dense depth-guided image, while simultaneously establishing a ground elevation grid and an occupancy probability grid. To suppress metallic reflections and dust noise, a secondary screening process is introduced using a reflection intensity threshold and spatiotemporal consistency constraints. The final output is fused environmental perception monitoring data containing 3D point sets, normals, intensity, timestamps, and image texture indexes, providing a spatiotemporally consistent and geometrically accurate environmental description for subsequent segmentation and attribute labeling.

[0017] Step S13: Perform contour segmentation processing on the environmental perception image based on the image binarization parameters of the fused environmental perception monitoring data to generate environmental perception image segmentation data; In this embodiment of the invention, based on the fused image, depth, and point cloud information, denoising and edge-preserving filtering are first performed on the depth and brightness channels respectively. A bilateral filter is used to suppress high-frequency noise while preserving geometric edges. Hole filling and depth confidence mapping are performed on the depth map, with confidence calculated jointly by depth validity, reflection intensity, and temporal consistency. Binarization adopts a joint judgment strategy of global threshold and local adaptive threshold: first, the global threshold is determined by histogram segmentation, and then the threshold is corrected within a local window using local mean or median statistics. Differences between the two are resolved by using depth gradient and normal abrupt change as weights to form an initial binary mask. Morphological opening operations are sequentially performed on this mask to remove isolated noise points, and closing operations are performed to fill small holes. Subsequently, the area and shape of connected components are verified, and connected component analysis is used to remove small-area false targets and label multi-scale candidate regions. Edge detection obtains an edge confidence map by using a multi-channel gradient operator that combines brightness and depth gradients. Based on high-confidence edges, a chain code contour tracking algorithm is applied to generate an ordered vertex sequence, and then spline fitting is used to smooth the contour while preserving the curvature information at the corners. To address occlusion and short-term viewpoint changes, contour extrapolation and repair are performed on the time axis using pose compensation and point cloud registration methods between adjacent frames to maintain segmentation coherence. The final output environment-aware image segmentation data includes the 2D contour of each segmented region, the polygon vertex sequence, the corresponding depth boundary point set, the normal distribution description, the confidence index of each region, and the timestamp and sensor extrinsic parameter index, satisfying the topological and geometric constraints required for subsequent 3D reconstruction and attribute calculation.

[0018] Step S14: Based on the loading operation requirement data, perform loading attribute feature identification processing on the environmental perception image segmentation data to generate loading attribute environmental perception monitoring data.

[0019] In this embodiment of the invention, the loading attribute labeling process is performed on the environmental perception image segmentation data region by region using the task requirement vector. First, a 3D reconstruction is performed on each segmented region: boundary point clouds are generated by mapping the segmentation contour and depth, and the normal and local curvature are obtained in the local neighborhood through principal component analysis. Then, a height grid is constructed in the segmented region using a voxelization method, and the volume is estimated by integral calculation of the height grid. Loading accessibility assessment is performed by performing discrete attitude checks on candidate contact directions using the inverse kinematics model of the loading mechanism. The checks include joint angle limits, mechanical interference judgment, and stability margin constraints of the lifting attitude. Spatial locations that do not meet any of the constraints are marked as unreachable. Cutting resistance is estimated on the material surface: the cutting resistance per unit depth is calculated by calling the friction and cohesion regression mapping obtained from the experiment based on the surface normal distribution, curvature, and brightness texture features. The resistance is then integrated along the feed trajectory in combination with the bucket geometry to obtain the expected work required and the maximum shovelable depth. Mass estimation is calculated using volume and material density benchmarks combined with humidity or reflectivity correction factors. The center of gravity position is calculated by combining the volumetric center of gravity with the remaining volume after shoveling for stability assessment. Stability assessment uses the minimum distance projected from the center of gravity to the vehicle support polygon as a metric, calculates safety margin values, and marks high-risk areas with thresholds. Dynamic risk assessment spatially convolves the obstacle probability field of short-term target tracking across segmented regions to generate a risk heatmap. Comprehensive operational indicators construct attribute scoring functions based on efficiency, energy consumption, and stability weights in task requirements. For each segmented region, a standardized attribute set is output, including accessibility mask, estimated volume and mass, surface normal and slope distribution, maximum shoveling depth and expected work, stability margin value, dynamic risk value, and region confidence label, recorded with timestamps and coordinate indices as direct constraints and scoring inputs for global path planning and local contact trajectory optimization.

[0020] Furthermore, step S2 includes the following steps: Step S21: Obtain the design parameters of the skid steer loader; In this embodiment of the invention, the design parameters of the skid steer loader are obtained through a standardized parameter extraction process. The parameter set includes three-dimensional geometric descriptions of components, coordinates of connection points, articulation types, center of mass position and mass distribution, material properties, hydraulic cylinder diameter and stroke, valve rated flow and pressure, wheelbase, vehicle moment of inertia, installation offset, and manufacturing tolerances. First, the dimensions of each parameter are standardized and units are converted. A parameter matrix is ​​established according to the machine's base coordinate system, and the coordinate system transformation relationships are recorded. Three-dimensional geometry is expressed using geometric primitives (patch, surface equations, and boundary representations). The volume, surface area, and center of mass coordinates of each component are calculated analytically, and the local mass distribution is obtained using volume segmentation. The parameters of hydraulic and actuator components are annotated, defining the rated operating point, short-term maximum output, temperature-related performance degradation coefficient, and fatigue life index. Finally, a structured parameter package is generated as a parameter document table for subsequent kinematic and dynamic model and capability boundary derivation.

[0021] Step S22: Analyze the connection relationship of loading components based on the design parameters of the skid steer loader to obtain the connection relationship data of the loading components, and perform kinematic constraint analysis of the loading components through the connection relationship data to generate the kinematic constraint data of the loading components. In this embodiment of the invention, the connection relationship of the loading components is analyzed according to the design parameters. First, a geometric primitive model is established for all components, represented by patches and boundaries, and the three-dimensional coordinates and joint axis are recorded for each hinge point. A homogeneous transformation matrix method is used to establish a relative pose expression for each joint, and a parameterized joint table lists the degrees of freedom of rotation or displacement and the joint limit values. A mechanism topology diagram is constructed using the component's centroid and hinge points as graph nodes, and hinges and slide rails as graph edges. This topology diagram is used to automatically generate a set of kinematic constraint equations. The Jacobian matrix of the constraint equations is derived for the closed-chain structure, and the forward kinematics, inverse kinematics, and constraint satisfaction solutions are solved using a numerical iteration method. The iterative method uses Newton-Raphson iteration combined with a linear solution subroutine to ensure convergence. Joint limits and singular poses are identified, and the rank loss of the constraint Jacobian is analyzed using singular value decomposition, and the singular domain boundary is recorded. To support collision detection and interference verification, a hierarchical bounding box (OBB) tree is constructed, and swept-body collision detection is performed on key pose sets. The swept bodies are generated from joint trajectory samples and their spatial intersection with the environment model is detected. All analysis results form connection relationship data for the loading components, including joint types and parameters, constraint equations, Jacobian matrix templates, joint limits and singular domain annotations, interference thresholds, and collision rejection vectors. Based on this data, linkage kinematic constraint analysis is performed. The constraint consistency of representative poses is verified using a numerical solver, and the inequality constraint set is recorded. Finally, the linkage kinematic constraint data of the loading components is output for subsequent dynamic modeling and path feasibility verification.

[0022] Step S23: Analyze the loader actuators according to the skid steer loader design parameters to obtain loader actuator data, and perform actuator capability profile abstraction processing on the loader actuator data to obtain loader actuator capability data; In this embodiment of the invention, the physical parameters of the actuators, from the hydraulic cylinder, pump, and valve to the drive motor, are analyzed. First, a force-pressure and velocity-flow mapping model is established for the hydraulic cylinder. This model is mapped onto standard operating points and solved using static equilibrium formulas and fluid continuity equations. Valve flow coefficients, pipeline pressure drop, and viscous losses are considered to obtain the actual output relationship. A first- or second-order transfer function model is constructed for the transient response of the electro-hydraulic components, and the time constant and damping ratio are recorded. For hysteresis behavior, a nonlinear hysteresis description is used, and dead zone boundaries are marked with piecewise functions. The continuous output capability, short-term peak capability, and thermal capacity constraints are represented in a multi-dimensional discrete manner to generate a capability grid. The horizontal axis of the grid represents the execution speed or displacement rate, and the vertical axis represents the output force or torque. The sustainable power, maximum allowable duration, and thermal accumulation coefficient are recorded at each grid point. Fatigue life is mapped based on stress cycle fatigue accumulation. The Miner's rule is used to derive a life loss estimate, and a life weight is added to the capability grid. The measurement resolution, noise spectral density, and delay characteristics of the sensor interface are labeled for observation function and filter design. The final result is a set of actuator capability data, including force-velocity capability surfaces, transient response templates, thermal and lifetime boundaries, inequality constraints, and observation accuracy specifications. This set is used as a rigid constraint in dynamic simulation and trajectory feasibility assessment to ensure that the trajectory does not exceed the actuator capability boundaries.

[0023] Step S24: Perform loading and lifting characterization analysis based on the skid steer loader design parameters to obtain loading and lifting characterization data; In this embodiment of the invention, the loading and lifting characterization analysis employs a discrete attitude scanning method. Mesh sampling is performed on the lifting angle and bucket attitude within a defined range, and static analysis and stability assessment are conducted at each sampling point. The static analysis calculates the required output force of the hydraulic cylinder to maintain the attitude based on the lever arm principle, and simultaneously calculates the projection position of the vehicle's center of gravity onto the wheel support polygon. The minimum distance between the projection and the support boundary is used as a stability margin index. The maximum mass that the bucket can bear at each height sample is obtained through inverse kinematics of torque balance, and this value is used to construct a lifting height-lifting capacity mapping. To reflect short-term dynamic effects, an inertial correction term is added to the static curve. The inertial correction term estimates the transient additional load using the acceleration spectrum and the equivalent mass of the bucket, and a dynamic margin surface is formed using a correction coefficient. The generated lifting characterization data includes a height-capacity discrete point set, stability margin contour lines, a dynamic margin correction matrix, and key attitude descriptions. It also records the center of gravity displacement contribution and symmetry imbalance factor under each attitude, used for rigorous determination of the safety margin during the lifting stage during path planning. This characterization data serves as the physical boundary of the lifting action during implementation, and is used for trajectory verification and local contact strategy constraints.

[0024] Step S25: Based on the kinematic constraint data of the loading components, the capability data of the loader actuators, and the loading and lifting characterization data, perform multibody dynamics modeling of the loader to obtain the multibody dynamics model of the loader; In this embodiment of the invention, a multibody dynamics model is constructed based on linkage kinematic constraint data, actuator capability data, and lifting characterization data. The modeling framework adopts the Lagrange equation, and the generalized coordinates cover the vehicle body degrees of freedom and the loading mechanism joint degrees of freedom. The mass matrix M(q) is derived, and the Coriolis-eccentric term is calculated. The equation includes the gravity term G(q), and the actuator input τ and contact force term are added to the right-hand side of the equation. The actuator dynamics are incorporated into the equations using a nonlinear first-order dynamic model with hysteresis and saturation. Valve flow dynamics are coupled to the cylinder velocity equations through the flow-pressure difference relationship. Closed-loop constraints are applied using the Lagrange multiplier method, with the multipliers solved along with the coordinates during numerical integration. An implicit numerical integration scheme (such as a variation of the midpoint) is used to ensure the numerical stability of the rigid system. The contact model employs an elastic-damped normal model combined with a Coulomb friction tangential model, and contact determination is based on nearest neighbor search and hierarchical spatial indexing. To meet real-time requirements, the model is reduced in order: modal truncation is performed on the high-frequency elastic modes, and a singular perturbation method is used to approximate the hydraulic transmission delay. The lifting subsystem retains key low-order dynamics to maintain mechanical accuracy. The output loader multibody dynamics model includes continuous-time dynamic equations, constraints, actuator sub-models, contact and friction models, and reduced-order mappings, serving as the mathematical basis for trajectory optimization, energy consumption estimation, and stability analysis.

[0025] Step S26: Perform adaptive parameter identification and optimization processing on the multibody dynamics model of the loader to obtain an optimized multibody dynamics model of the loader.

[0026] In this embodiment of the invention, adaptive parameter identification and optimization of the operating characteristics of the established multibody dynamics model are implemented. The identification process is divided into two stages: offline identification and online adaptation. Offline identification involves designing an excitation sequence to cover the working domain of each joint, recording joint angles, velocities, accelerations, hydraulic pressures, and flow rates. A weighted least squares solver is used to estimate inertia, friction terms, flow coefficients, and valve time constants. During the solution process, a robust weighting matrix is ​​used to weight measurement noise to improve the robustness of the estimation. Online adaptation uses recursive least squares with a forgetting factor or extended Kalman filtering to continuously update the core parameters. A projection operator is used to constrain parameter estimation within the physically feasible region to prevent divergence, thereby obtaining an optimized multibody dynamics model for the loader.

[0027] Furthermore, step S25 includes the following steps: Step S251: Perform lifting attitude dimension analysis based on the loading lifting characterization data to obtain lifting attitude dimension data; In this embodiment of the invention, based on the lifting characterization data, grid sampling is first performed in the attitude domain according to the joint variable space of the lifting mechanism. The sampling dimensions include boom extension, boom segment angle, bucket pitch angle, and overall vehicle tilt angle. For each sampled attitude, the homogeneous transformation relationship of the aforementioned geometric primitives is calculated, an attitude vector sequence is constructed, and singular value decomposition and principal component analysis are performed to identify the dominant degrees of freedom and their energy distribution in the attitude set, thereby forming a basis vector set for the lifting attitude dimension. For each basis vector, the corresponding stability margin gradient, actuator torque sensitivity, and contact geometry change rate are calculated. The attitude sensitivity matrix is ​​obtained by numerical difference method, and the norm is used to measure the sensitivity of key directions. Accessibility and interference scanning are also performed in the attitude dimension. Collision detection is performed on discrete attitudes using bounding box hierarchical partitioning, and the joint limit approach rate and interference margin are recorded for each attitude. The output results include: the attitude basis vector set, the stability margin curve corresponding to the basis vector, the attitude sensitivity matrix, the boundary of the infeasible attitude set, and the dimensionality reduction mapping relationship of the attitude representation, which are used by the subsequent decision module to limit the lifting action space and provide a metric for attitude selection.

[0028] Step S252: Perform linearization characteristic analysis of lifting height and lifting capacity based on the loading lifting characterization data to obtain lifting height-capacity characteristic data; In this embodiment of the invention, based on the height-capacity curve of lifting characterization data, the original nonlinear relationship is linearized using piecewise linear regression. First, the height interval is divided into several sub-segments. The linear relationship of each sub-segment is fitted using the least squares method, and the slope, intercept, and fitting residual of each sub-segment are calculated. The residual statistics are used to determine whether the segment boundaries need to be redefined. Then, continuity constraints are applied to smooth the endpoints of adjacent sub-segments, ensuring that the piecewise linear mapping has no jumps at the nodes. To account for dynamic effects, a dynamic correction term is introduced on top of the static linearization. By comparing typical lifting acceleration curves with the transient response of the hydraulic cylinder, a height-related acceleration sensitivity coefficient is derived and added to the linear mapping in the form of a first-order hysteresis filter, forming a family of linear approximate functions for height-capacity. For each approximate linear segment, the confidence interval and maximum error bound are calculated. The output includes linear segment parameters, a dynamic correction template, error bounds, and node smoothing coefficients. The subsequent planning module uses these linear expressions as constraints on the bucket lifting capacity during path generation and uses the error bounds for safety margin constraints.

[0029] Step S253: Based on the kinematic constraint data of the loading components and the capability data of the loader actuators, perform nonlinear boundary constraints and architecture modeling of the loader's multibody dynamics to obtain the loader's multibody dynamics architecture model; In this embodiment of the invention, a nonlinear boundary constraint set and system architecture model are established based on the kinematic constraint data and actuator capability data. Joint limits, geometric interference thresholds, actuator saturation surfaces, thermal accumulation limits, and lifetime constraints are compiled into a constraint set in the form of inequalities. Complementary conditions are added for contact scenarios to reflect the bidirectional constraints of contact state and contact force. The Lagrange multiplier method or complementary relaxation method is used to simultaneously model the equality and inequality. After discretization, a constraint projection step based on quadratic programming is used in the numerical solver to ensure that each integration step satisfies the boundary constraints. The system architecture adopts a hierarchical multibody representation: the base layer consists of the vehicle body rigid body and drive structure; the intermediate layer consists of the multi-joint chain of the loading mechanism; and the top layer consists of the actuator model and thermal / lifetime submodules. Each layer is connected through a formal interface matrix, which describes the points of action, transmission coupling, and constraint mapping relationships. The system performs parameter consistency checks on the architecture model, evaluates the collaborative response of different sub-modules under extreme conditions, and generates a boundary constraint mapping table. The output includes constraint equation templates, complementary condition statements, actuator saturation mappings, and architecture interface matrices, which serve as structural constraint inputs for dynamic solution and trajectory optimization.

[0030] Step S254: Map the lifting posture dimension data and lifting height-capacity characteristic data to the loader multibody dynamics architecture model to perform loading and lifting dynamics characteristic mapping processing to obtain the loader multibody dynamics model.

[0031] In this embodiment of the invention, lifting attitude dimension data and lifting height-capacity characteristic data are mapped to a multibody dynamics architecture model. First, a parameterized representation is specified for lifting-related variables in generalized coordinates. The high-dimensional attitude is expressed by a linear combination of the dimensionality-reduced attitude basis vectors, thereby replacing the high-dimensional joint description with low-dimensional parameters in the dynamic equations to reduce computational load. Simultaneously, the piecewise linearized height-capacity expression is added as an inequality constraint to the drive input constraint set. Specifically, in each height segment, the upper bound of the execution torque or hydraulic pressure is written as a linear function and used as a linear constraint in the optimizer. Dynamic correction terms are incorporated into the state equations as additional state variables, forming a linear parameter time-varying or affine parameterized model. Parameter updates based on the loaded mass estimation are performed on the mass matrix and inertia tensor: during the calculation process, the overall center of mass position and inertial elements are corrected by translating the center of mass term for the loaded mass of the bucket. Subsequently, numerical integration is used to verify the stability margin of the corrected system under extreme attitudes. A multibody dynamics model with lift mapping is formed, which includes dimensionality-reduced attitude expression, piecewise linear force constraints, dynamic correction sub-model, corrected mass and inertia matrices, and constraint matrices.

[0032] Furthermore, step S252 includes the following steps: Based on the loading and lifting characterization data, dynamic lifting reduction processing is performed to obtain dynamic lifting reduction data of the loader. Then, the lifting height and lifting capacity are linearized and approximated using the dynamic lifting reduction data of the loader to obtain lifting height-capacity characteristic data.

[0033] In this embodiment of the invention, dynamic order reduction is performed on the dynamic equations of the full-order lifting subsystem based on the Lagrange framework. First, frequency and time domain response analysis is conducted on the full-order model under typical operating conditions to identify the dominant modes. Impulse response or pseudo-random binary sequence excitation is used to sample the response of joint angles, hydraulic cylinder displacement, and cylinder pressure. A finite-time Hankel matrix is ​​constructed and singular value decomposition is performed to obtain the Hankel singular value spectrum. The number of retained modes is determined based on the breakpoints of the singular value amplitudes. The retained modes are projected using orthogonal projection to obtain the projection matrix. After projection, a reduced-order state equation is formed, and the output consists of a reduced-order matrix set, mode vectors, and corresponding time constants. To ensure numerical stability and energy consistency, equilibrium truncation or energy preservation orthogonalization is used during the projection process. Singular perturbation decomposition is applied to the hydraulic nonlinear saturation to approximate fast dynamics as algebraic constraints, while slow dynamics are retained as state equations, thus maintaining the passivity and asymptotic stability of the reduced-order model. The reduced-order results include modal shape, natural frequency, damping ratio, modal contribution rate, upper bound of the estimated approximation error, and modal mapping matrix. Additionally, parameter update rules for online correction are generated. These rules map the current bucket load and center of mass position to the mass matrix of the reduced-order model via a linear correction term, thus updating the dynamic response prediction with extremely low overhead during operation. The reduced-order dataset serves as the computational benchmark for subsequent trajectory optimization and real-time simulation, replacing the full-order model, to achieve a high-fidelity approximation of lifting dynamics within a limited computational budget. A family of linear approximations of lifting height and lifting capacity is constructed using the steady-state mapping relationship of the reduced-order model. First, the quasi-static equilibrium of the reduced-order model is solved at several discrete points within the height range. At each discrete height, the state derivative is set to zero, and considering actuator saturation and thermal / lifetime boundaries, the corresponding required driving force or hydraulic pressure is solved, forming a height-driving force mapping point set. Piecewise linear regression is performed on the mapped point set: segments are divided using the fitting error and Hankel residual as the segmentation criteria. For each segment, linear coefficients (slope and intercept) are calculated, along with the confidence interval and maximum deviation of the linear approximation. To reflect dynamic transition characteristics, the primary time constant related to lifting in the reduced-order model is extracted as a first-order dynamic correction term and appended to the driving response of each linear segment using a first-order hysteresis model, forming a linear approximation with dynamic correction. The generated lifting height-lifting capacity characteristic data includes a linear segment parameter table, dynamic time constants, error limits, and safety margins for application. A rapid update mechanism is also included: when the loaded mass estimate changes, the slope of the linear segment is adjusted analytically using an inertia correction factor introduced by the mass. The adjustment formula is based on the lever arm's linear response and the center of mass translation, ensuring that the corrected linear parameters can be obtained within a constant time under online conditions. The final output is a piecewise linear constraint set used by the trajectory optimizer. The optimization problem is written in the form of linear constraints, and dynamic correction terms are used for compensation during the control execution phase, thus balancing computational efficiency and mechanical accuracy.

[0034] Furthermore, as an embodiment of the present invention, reference is made to... Figure 2 As shown, Figure 1 A detailed flowchart of step S3 is shown below. In this embodiment, step S3 includes the following steps: Step S31: Establish a panoramic perception coordinate system for loading attributes based on the environmental perception monitoring data of loading attributes, so as to obtain the panoramic perception coordinate system for loading attributes; In this embodiment of the invention, the process of establishing a panoramic perception coordinate system based on loading attribute environmental perception monitoring data includes coordinate reference definition, sensor frame fusion, and terrain datum calibration. First, the loader body reference point is taken as the origin, the forward axis along the loader's forward direction is defined as the X-axis, the lateral axis as the Y-axis, and the vertical axis as the Z-axis. The attitude angles provided by the inertial measurement unit are used to correct the body reference direction to ensure initial alignment between the vehicle's heading and the panoramic coordinate system. Intrinsic and intrinsic parameter mapping is performed on the multi-source sensors, and the rigid body transformation matrices of each sensor are used for coordinate transformation. A point cloud registration algorithm is used to register the short-range lidar point cloud and the depth camera point cloud within the body coordinate system. The registration process primarily uses the iterative nearest-point method, supplemented by initial pose estimation based on key features to accelerate convergence. The ground datum is extracted from the point cloud through RANSAC plane fitting, and a local zero-elevation surface is defined using the fitted plane. During datum extraction, non-ground points determined by the normal angle threshold are removed to eliminate the influence of material piles and obstacles. After registration, the multi-source data is represented in a panoramic coordinate system as a voxel grid or elevation raster. The voxel resolution is set proportionally according to the scale of the working area and the planning accuracy requirements. The voxels record the occupancy probability, average elevation, normal statistics, and texture index. Temporal registration and timestamp interpolation are used to eliminate motion blur for continuous temporal frames, generating a panoramic perception field with time stamps. A short-term probability decay layer is maintained for dynamic obstacles to reflect temporal correlation. The final output is a panoramic perception coordinate system file containing rigid body transformation matrices, terrain datum equations, occupancy rasters, and elevation rasters. This coordinate system serves as a unified spatial reference for subsequent path planning and local contact assessment.

[0035] Step S32: Perform a sliding path feasible region analysis based on the panoramic perception coordinate system of the loading attributes to obtain sliding path feasible region data; In this embodiment of the invention, the process of conducting skid path feasibility domain analysis in the panoramic perception coordinate system of loading attributes includes occupancy assessment, vehicle kinematic envelope calculation, terrain constraint determination, and actuator capability comparison. First, the panoramic occupancy grid undergoes spatial dilation: after morphological dilation using the vehicle's outer rectangle and safety edge, occupancy convolution is performed on the grid to generate occupancy masks for the vehicle in different orientations. Terrain constraints are calculated using the elevation grid to determine the local slope and aspect of each grid point. Lateral and longitudinal tilt are quantified using normal distribution and elevation neighborhood fitting. Grids are marked as slope infeasible or slope feasible using slope thresholds and stability margin thresholds. Skid characteristics are calculated using the loader's dynamics model to determine the required lateral traction and wheel-end torque. For each candidate path, the skid angle and lateral friction requirement are estimated for a micro-segment. This friction requirement is compared with the ground friction coefficient field estimated from perception. When the friction requirement exceeds the ground friction capacity or the torque of the actuator in that attitude exceeds the capacity grid, the corresponding grid is marked as dynamically infeasible. Temporal convolution is performed on the dynamic obstacle probability field on the grid to quantify short-term obstacle risk. Grids with high risk have their infeasibility weight increased based on a threshold. The final output of the slip path feasible domain data includes a set of feasible grids, infeasibility cause labels (occupancy, slope, friction, actuator capability, dynamic risk), and a feasibility cost map oriented towards path search.

[0036] Step S33: Perform path completion and slip characteristic optimization processing on the slip path feasible region data to obtain slip optimized feasible path data; In this embodiment of the invention, the process of performing path completion and slip characteristic optimization on the feasible region data of the slip path includes key node extraction, continuous curve fitting, dynamic response fitting, and curvature-constrained smoothing. First, skeletonization or centerline extraction is performed in the feasible region to obtain a topology map of the free channel. Several key intermediate nodes are extracted from the channel topology based on the shortest connectivity and obstacle gaps. Parameterized line completion is performed between nodes using cubic splines or B-splines. During the line completion process, positional and tangential continuity constraints are applied to the nodes to ensure path smoothness. The completed curves are used for forward dynamic simulation with discrete arc length steps. The simulation uses the optimized multibody dynamics model to calculate the desired lateral slip, lateral force, and actuator torque at discrete points. For segments where lateral slip or torque exceeds limits, a constraint optimizer is used to adjust the spline control point positions. The optimization objective is to minimize the rate of curvature change and the L2 distance from the initial path, while simultaneously satisfying curvature and curvature derivative constraints. The curvature constraint is jointly determined by vehicle stability and the slip critical value. To ensure the robustness of the path in actual tracking, a path redundancy band is introduced to describe the adjustable buffer bands on both sides of the path, and a robustness analysis is performed on the line completion results. Path segments with insufficient robustness trigger local node redistribution or shape correction. Finally, the slip optimization feasible path data that meets the dynamic response characteristics is output. This path data includes continuous parameterized curves, control point sets, curvature constraint annotations, and local robustness scores.

[0037] Step S34: Based on the optimized multibody dynamics model of the loader and the environmental perception monitoring data of loading attributes, perform global path cost feature analysis of skid loading to obtain global path cost feature data of skid loading. The global path cost feature data of skid loading includes path length cost data, skid loading energy consumption cost data, path slope cost data, short-term dynamic obstacle probability penalty cost data, dynamic environmental risk penalty cost data, loading lifting stability cost data, and short-term operation efficiency cost data. In this embodiment of the invention, the process of cost feature analysis of the global path of skid loading based on the optimized multibody dynamics model of the loader and environmental perception monitoring data of loading attributes adopts a method combining discrete integration and local simulation to calculate multidimensional cost terms. The path is discretized into equidistant micro-segments, and the instantaneous traction force, bucket attitude energy consumption, and actuator flow demand are estimated using a dynamic model on each micro-segment. The path length cost is obtained by summing the arc lengths of the micro-segments. The skid loading energy consumption cost is obtained by estimating the instantaneous power (pressure multiplied by flow rate) of the hydraulic system in each micro-segment and integrating it over time to obtain the total energy consumption. The path slope cost is a composite index of slope energy consumption and stability risk formed by multiplying the micro-segment slope and slope sensitivity coefficient by the cumulative sum. The short-term dynamic obstacle probability penalty is obtained by integrating the probability grid on the path, with high-probability areas incurring significant penalties. The dynamic environmental risk penalty is calculated based on the slope change rate, surface roughness, and humidity estimates to determine the probability of skid instability and convert it into a penalty value. The loading and lifting stability cost is obtained by calling the lifting characterization data at the expected lifting position on the path, calculating the minimum distance between the centroid projection and the support polygon during lifting, and accumulating it according to the inverse safety margin function. The short-term operation efficiency cost is obtained by estimating the operation cycle and time cost per unit path segment using the expected number of scoops, average loading capacity, and round-trip time. To ensure consistency in the dimensions of cost terms, the cost value is normalized on a global scale, and local maximum and mean values ​​are recorded. The output of global path cost feature data for slip loading includes the micro-segment distribution curves of each cost term, cumulative cost value, and normalization index, providing accurate numerical input for subsequent cost model construction.

[0038] Step S35: Establish the mapping relationship between the global path and cost of skid loading based on the global path cost feature data of skid loading, generate a preliminary global path cost model of skid loading, and perform adaptive weight parameter adjustment on the preliminary global path cost model of skid loading through loading operation demand data to obtain the global path cost model of skid loading. In this embodiment of the invention, during the process of establishing the mapping relationship between the global path and cost of skid loading, various cost items are normalized based on path length cost data, skid loading energy consumption cost data, path slope cost data, short-term dynamic obstacle probability penalty cost data, dynamic environmental risk penalty cost data, loading and lifting stability cost data, and short-term operational efficiency cost data. The normalization process employs a piecewise scaling method, using the maximum and minimum values ​​of each cost item in historical or simulation data as boundaries to ensure that cost items of different dimensions are compared uniformly under the same scale. Subsequently, a linear weighting function is constructed to realize the mapping relationship between the path and the cost. This function uses the normalized value of each cost item as input to calculate the comprehensive cost value and maps it one by one onto the path data set to generate a preliminary global path cost model for skid loading. The preliminary model at this stage reflects the basic contribution relationship of different cost items but does not dynamically adapt to actual task requirements. After obtaining the preliminary cost model, adaptive weight parameters are adjusted in conjunction with loading operation requirement data. The loading operation requirement data includes multiple operation-oriented parameters such as efficiency priority, energy consumption limitations, operational stability constraints, and environmental risk control. The system parses demand data into a priority matrix and uses a weighting algorithm to map the priority matrix to the weight parameters of each cost item. For example, when the operational demand emphasizes improving loading efficiency in a short period of time, the system increases the weight of the short-term operational efficiency cost item and decreases the weight of the energy consumption cost item. When the operational demand involves complex slopes or high-risk scenarios, the system increases the weight of the path slope cost and the dynamic environmental risk penalty cost to ensure the stability and safety of the operation process. The weight adjustment process adopts an exponential smoothing and real-time feedback correction mechanism: on the one hand, it achieves a smooth transition by weighted balancing of historical task preferences and current operational instructions; on the other hand, it adjusts the weight offset through real-time path execution feedback, enabling the model to respond instantly to environmental changes and operational demand switching. The resulting global path cost model for skid loading not only includes the mapping function between path and cost but also records the weighting factors, normalized boundary parameters, and adjustment strategy rules of each cost item. This model can be used as an evaluation function in subsequent path search to achieve quantitative comparison and optimization of different candidate paths, thereby outputting the skid loading global path that best meets the task objectives in a dynamic environment.

[0039] Step S36: Transfer the slip optimization feasible path data to the slip loading global path cost model for global path search processing of slip loading, and generate slip loading global path data.

[0040] In this embodiment of the invention, after inputting the slip-optimized feasible path data into the slip-loading global path cost model, a global path search is performed. During the path search process, the system calculates the global cost value of different candidate paths through a cost function and sorts and filters them. A two-stage retrieval and refinement strategy is adopted to balance real-time performance and optimality. In the first stage, a heuristic graph search algorithm is applied at the grid level to obtain candidate discrete paths. The heuristic function is composed of Euclidean distance and dynamic obstacle probability weighting to improve the sensitivity to dynamic risks. The candidate paths are quickly sorted at the discrete layer using the cost model, and several low-cost paths are retained as refinement candidates. In the second stage, the candidate paths are parameterized as spline curves and refined in continuous space using a constrained numerical optimization method. The optimization objective is the continuous integral value of the cost model. Constraints include curvature and curvature derivative restrictions, velocity acceleration boundaries, and actuator torque boundaries. The numerical solution adopts sequential quadratic programming or interior point method and a hot-start strategy is used to improve the convergence speed. For each optimized path, a forward simulation based on a multibody dynamics model is performed for verification. The simulation output includes trajectory tracking error, energy consumption estimation, and stability margin. If any indicator exceeds a threshold, local replanning is triggered, and the process returns to the second-stage iteration. The path with the minimum cost and that has passed verification from the refined path set is used as the global path data output for slip loading. The path is then time-parameterized to generate a velocity-attitude reference curve for the controller to execute. When the perception layer detects dynamic obstacles or a significant update to the cost model input, incremental replanning in the first stage is triggered to maintain online response capability.

[0041] Furthermore, step S33 includes the following steps: The loader's slip characteristics were analyzed based on the optimized multibody dynamics model to obtain slip characteristic data. The intermediate nodes of the feasible region data of the skid path are selected for skid path line supplementation processing to obtain skid path line supplementation data. Then, the skid path line supplementation data is subjected to loader skid characteristic line fitting optimization processing based on the loader skid characteristic data to obtain skid optimized feasible path data.

[0042] In this embodiment of the invention, based on the optimized multibody dynamics model, a parameter scanning domain for slip analysis is first constructed, including the vehicle speed range, left and right drive differential speed range, longitudinal load state, ground friction coefficient and terrain slope range, and bucket lifting height and attitude combination. Time-domain numerical simulation is performed on each parameter vector. The numerical integration uses an implicit integration scheme to ensure the stability of the rigidity term. The contact force is solved using a complementary constraint method and nearest-neighbor contact determination to obtain the positive pressure and tangential shear force between the tire or drive surface and the ground. During the simulation, the longitudinal slip ratio and lateral slip angle are calculated for each wheel / side set. The longitudinal slip ratio is expressed as the ratio of the wheel circumferential speed to the ground forward speed, and the lateral slip angle is calculated as the arctangent of the wheel lateral speed and longitudinal speed. Based on the ground-wheel / track contact, the instantaneous traction force, lateral force, traction efficiency loss, and ground settlement energy are evaluated to form a micro-segment energy consumption spectrum and cumulative slip energy consumption. The simulation results are post-processed in both the frequency and time domains: the dominant frequency and damping ratio are extracted in the frequency domain, and the slip angle distribution, lateral force utilization rate, and longitudinal traction coefficient are statistically analyzed in the time domain. The final output slip characteristic data includes: slip angle surface mapped by velocity and differential speed, lateral and longitudinal force ratio matrix, slip energy loss curve per unit distance, critical differential speed threshold (when the lateral force utilization rate is continuously greater than the set upper limit), cumulative ground settlement rate, and corresponding traction attenuation factor. This data is used for dynamic feasibility assessment, energy consumption estimation, and slip warning threshold setting for real-time control during the path phase, providing a rapid calculation method for the online planning module in the form of numerical lookup tables or polynomial fitting. A series of intermediate nodes are obtained in the feasible region of the slip path through skeleton extraction or topology analysis. The nodes are mainly key points of the channel, projection points of obstacle boundaries, and points close to the material pile. An initial broken path is constructed using the set of intermediate nodes, and a line completion technique is used to obtain a continuous curve. The curve after completion is discretized into micro-segments based on arc length. Forward simulation is performed on each micro-segment using an optimized loader multibody dynamics model. The simulation output includes instantaneous slip angle, lateral force demand, actuator torque demand, and instantaneous energy consumption. The simulation output is compared with slip characteristic data to calculate the slip overload factor (lateral force demand / friction limit) and energy consumption penalty for each micro-segment. For curve segments with overload or high energy consumption concentration areas, local control point adjustments are implemented. These adjustments minimize the objective function, which consists of the curvature L2 norm, the slip overload penalty, and a weighted sum of path length. Constraints include spatial boundary limits, curvature derivative limits, and inequality constraints ensuring the actuator torque does not exceed the capacity grid. A constrained sequential quadratic programming strategy is used to solve the problem, employing previous solutions as a warm start to accelerate convergence. If any micro-segment still exceeds the slip or torque limits after adjustment, the local optimization radius is expanded or intermediate nodes are rearranged to change the path topology.The final output of the slip optimization feasible path data is a parameterized continuous curve, with piecewise velocity profiles, predicted slip angle spectra, expected energy consumption and stability margin indices for each segment, and a buffer zone describing the redundancy space on the left and right sides of each segment, which is used for subsequent time parameterization and safety margin assessment of closed-loop tracking control.

[0043] Furthermore, as an embodiment of the present invention, reference is made to... Figure 3 As shown, Figure 1 A detailed flowchart of step S4 is shown below. In this embodiment, step S4 includes the following steps: Step S41: Extract the target loading image at each local loading moment based on the loading attribute environmental perception monitoring data and the skid loading global path data to obtain the target loading image data at the local moment. In this embodiment of the invention, the desired bucket pose for each local loading moment is pre-calculated along the global sliding loading path. The bucket pose is obtained from path time parameterization and the inverse kinematics solution of the loading mechanism. Based on the desired bucket pose and the extrinsic parameters of each onboard sensor, a frustum and perception window for each sensor are constructed. Motion compensation and viewpoint correction are performed on multi-source images using time-synchronized attitude information. Image processing first performs spatial filtering and hole filling on pixel-level distortion and depth noise, and then projects the depth point cloud onto each camera plane using projection mapping to obtain depth-guided region candidate boxes. The candidate boxes are clipped with bucket geometric bounding boxes and reachability boundaries to generate a subset of target loading images at local moments. Geometric extrinsic parameter fine calibration is performed on subsets from different sensors at the same local moment. A multi-view dense registration method is used to run iterative nearest-point registration on the 3D point set, supplemented by normal consistency constraints, to eliminate slight perspective errors and temporal deviations. For occluded pixels, inter-frame extrapolation and multi-view fusion strategies are employed to recover missing regions. The fusion output includes calibrated color images, corresponding depth maps, a subset of 3D point clouds, a normal vector field, and a confidence index per pixel. The final output of local time target load image data is presented by time index, listing the image set, depth set, point cloud set, camera calibration index, and confidence sequence, providing a geometrically complete and spatiotemporally consistent local observation foundation for subsequent bucket contact trajectory analysis.

[0044] Step S42: Analyze the bucket loading contact trajectory based on the target load image data at a local time to obtain the bucket loading contact trajectory data at a local time; In this embodiment of the invention, when performing bucket loading contact trajectory analysis on target load image data at a local time, the image and depth point cloud are reconstructed with high precision in a local window region. The reconstruction output includes a closed triangular mesh, a normal field, and a local curvature field. Based on the triangular mesh, several candidate contact nodes are selected on the local surface. For each node, a parameterized intrusion curve is generated along several intrusion directions. The parameterized variables include the intrusion depth profile, the intrusion angle, and the intrusion velocity profile. For each intrusion curve, collision detection and contact criterion calculation are performed iteratively with small step sizes. Collision detection relies on bounding box hierarchical structure acceleration and nearest neighbor triangulation retrieval. During the determination process, the contact surface patch, normal distribution, and tangential friction distribution are calculated simultaneously. The above geometric information is mapped to optimize the multibody dynamics model of the loader. The dynamic simulation calculates the contact reaction force, the torque distribution along the bucket structure, the vehicle center of gravity offset caused by the contact force, and the instantaneous pressure / torque value required by the actuator at each discrete depth step. For each candidate trajectory, an objective function term is calculated: unit work, peak actuator load percentage, minimum stability margin, and trajectory time cost. These physical quantities are used to sort the candidate trajectories and eliminate those that violate actuator capabilities or the lower limit of the stability margin. The output local-moment bucket loading contact trajectory data includes a sequence of trajectory points, bucket pose at each point, corresponding contact force distribution, estimated load increment, required actuator pressure curve, and trajectory confidence index, providing physically consistent and verifiable input for swarm prediction and subsequent optimization.

[0045] Step S43: Perform group local time prediction processing of bucket loading contact trajectory based on local time bucket loading contact trajectory data and local time target load image data to obtain bucket loading contact trajectory prediction data; In this embodiment of the invention, when performing group local time prediction processing on local moment bucket loading contact trajectory data and local images, a trajectory set is first established using sequence parameterized trajectories as the basic unit. The trajectory parameter vector includes intrusion angle, intrusion depth, intrusion rate, estimated loading mass, and confidence level. A parallel physical-statistical hybrid prediction framework is adopted: physical prediction is based on a reduced-order dynamic model and the principle of volume conservation to derive the impact of short-time trajectories on local volume, while statistical prediction uses a state-space autoregressive model to perform time-related modeling of historical observation residuals. Multi-model parallel forward prediction is implemented on the heterogeneous trajectory set, using Monte Carlo sampling or particle representation to characterize uncertainties. The particles carry local volume distribution, centroid displacement, and surface roughness estimates during time progression, and maintain distribution representativeness through resampling at each time step. The output is a time-series trajectory prediction set, with each prediction containing time-series trajectory pose, expected removal volume, estimated contact force time series, and confidence level distribution. A cross-consistency test is further performed on the predicted trajectories: the local volume reduction caused by each predicted trajectory is spatially accumulated and conservatively compared with the maximum shovelable volume inferred from the local imagery to mark the over-allocation area and incorporate conflict resolution constraints in subsequent optimization. This predicted data provides temporal priors and uncertainty measures for establishing a dynamic model of cargo changes and optimizing swarm decisions.

[0046] Step S44: Establish a dynamic change state transition matrix and a dynamic change observation matrix of the load based on the bucket loading contact trajectory prediction data. The dynamic change state transition matrix of the load is based on the target load change relationship at each local loading moment derived from the bucket loading contact trajectory prediction data. The dynamic change observation matrix of the load reflects the natural deformation observation value of the target load at each local loading moment. In this embodiment of the invention, based on the predicted contact trajectory data of the bucket loading, a local state vector of the load is defined. The state vector contains several components: local remaining volume, local height field coefficients (displayed using orthogonal basis function expansion), local centroid position components, and local surface roughness indices. For each predicted trajectory, the increment of influence on the state vector is calculated within a short time step according to the principle of mass conservation. The increment calculation is based on the voxel integral of the trajectory sweep volume, and the amount of mass removal is allocated to each voxel. The changes at the voxel level are then mapped back to the changes in the basis function coefficients of the state vector. The linearized mapping under single-step action is approximated by a Jacobian matrix. The Jacobian matrix terms consist of the partial derivatives of the trajectory parameters with respect to the state components, and the partial derivative values ​​are obtained through finite difference or analytical derivation. The linear increments of each trajectory at the same time step are combined according to the superposition principle to obtain the overall state transition matrix A_k for that time step. The overall state transition matrix represents the change in the target load due to each bucket loading step, such as the change in the load at the moment after the bucket has finished loading the target load. During matrix construction, volume conservation and non-negativity constraints are explicitly verified. If the superposition result shows a negative volume or exceeds the upper bound of the local volume, the contribution of the corresponding trajectory is conservatively truncated and a correction term is added to the matrix to maintain physical consistency. The observation matrix H_k is derived based on the sensor observation model. The observation matrix represents the observations of the target load due to natural dynamic changes after loading, such as the natural loss error of sand and gravel after loading, and the natural fall of the load. The observation model linearly maps the depth histogram, local point cloud density, and normal statistics to the observation components of the state vector. The mapping Jacobian is obtained based on the projection of the observation operator onto the basis function coefficients. To suppress numerical defects, A_k and H_k are analyzed using singular value decomposition after construction, and Tikhonov regularization is applied to the critical directions. The regularization parameters are determined by both observation confidence and physical constraints. The process noise covariance Q_k and observation noise covariance R_k are estimated statistically based on sensor calibration errors and historical residuals to fully describe randomness and uncertainty. The final outputs are time series A_k, H_k, Q_k, and R_k, which are used for filtering and subsequent optimization.

[0047] Step S45: Based on the dynamic change state transition matrix of the load and the dynamic change observation matrix of the load, perform group local time correlation optimization processing on the bucket loading contact trajectory data to obtain the group local time optimized bucket contact trajectory data. In this embodiment of the invention, the dynamic change state transition matrix and the dynamic change observation matrix of the load are used as priors to perform group-level joint optimization on the set of bucket contact trajectories at local moments. First, in the prediction-update loop, extended Kalman filtering or unscented Kalman filtering is used to recursively estimate the state: the prediction step uses A_k and Q_k to predict the state and covariance of the next moment, and the update step uses H_k and R_k to correct the current observations into the state estimate and output the filtered mean and covariance. Based on the filtering results, a finite-time domain joint optimization problem is constructed, with the optimization window covering several future local moments. The decision variables include the intrusion depth time profile, intrusion rate, and slight adjustment of the intrusion angle for each candidate trajectory. The objective function is in the form of a weighted sum, with weights including maximizing the cumulative removed volume, minimizing the work per unit, minimizing the L2 norm of the overall centroid shift, and a collision penalty term between trajectories; the constraint set includes the execution element capability inequality, the lower bound of the lifting stability margin, the non-negativity of local volume, and the space occupancy mutual exclusion constraint. The solution employs a sequential quadratic programming approach within a recursive finite-time optimization framework. The quadratic subproblems are solved in each iteration by a quadratic programming solver with inequality constraints. The solver uses the previous solution as a warm-start to accelerate convergence. To enhance the robustness of the solution, a soft constraint on observation consistency is introduced into the optimization process. The filter covariance is used as a weight to incorporate the observation consistency bias as a quadratic penalty term into the objective function, balancing the discrepancies between model predictions and sensor observations. The solution results include the jointly optimized trajectory sequence, the state estimate updates at corresponding time steps, covariance convergence information, and the proposed corrections to the global path. A conflict allocation matrix is ​​also output for subsequent global coordination.

[0048] Step S46: Optimize the loading motion trajectory of the skid loader global path data based on the optimized bucket contact trajectory data at local moments to obtain the skid loader global optimized motion trajectory data; execute the skid loader loading intelligent control operation based on the skid loader global optimized motion trajectory data.

[0049] In this embodiment of the invention, the optimized contact trajectory sequence of the group at local moments is mapped back to the global path. The parameterization interval and insertion point of each local trajectory in the global path are determined. A continuity correction is established by calculating the difference in pose and velocity acceleration before and after insertion. A time reparameterization strategy is applied to the global trajectory to ensure at least second-order derivative continuity. The reparameterization objective is to minimize the jerk (change in velocity rate) and satisfy the time window and actuator response capability constraints. After reparameterization, a full-process forward simulation is performed on a reduced-order multibody dynamics model. The simulation output includes the slip angle time series, instantaneous contact force spectrum, actuator pressure and power curves, and lifting stability margin curve. If a local over-limit situation is found in the simulation, local re-optimization is performed on the local segment first. Local re-optimization adopts the same finite-time domain optimization framework but is limited to the affected interval, and actuator capability grid constraints and lifting safety margin lower limits are applied during the optimization process. After verifying that all dynamic and safety constraints are met, the final globally optimized motion trajectory data for skid loading is generated. The trajectory data includes spatial path parameterization, time-velocity-attitude profiles, bucket pose sequence, actuator reference values ​​(valve flow / pressure or torque profiles), and confidence indices for each step. The execution layer uses a parallel force-position closed-loop control strategy to execute this trajectory: in the contact phase, a force-priority control law is used to track the force reference and limit displacement error; in the non-contact phase, a position-priority control law is used to track the spatial trajectory. The control loop has a built-in fault monitor that observes the contact force residual, observation residual, and actuator response delay. When any residual exceeds a preset safety threshold, an immediate stop is triggered, and a local replanning process is initiated. The replanning updates the A_k / H_k sequence with the latest observations and resolves the swarm optimization problem within a limited time window to generate an alternative trajectory. Then, closed-loop tracking is resumed to ensure that the operation maintains efficiency and safety and meets physical constraints in a dynamic environment.

[0050] Furthermore, step S42 includes the following steps: Step S421: Based on the local time target load image data and the optimized loader multibody dynamics model, perform bucket loading characteristic analysis of contact point differences to obtain bucket loading characteristic data of contact point differences; In this embodiment of the invention, after 3D reconstruction and surface feature analysis of the target load image data at a local moment, the generated point cloud or triangular mesh is used as input. Combined with an optimized multibody dynamics model of the loader, contact point position difference analysis is performed on each candidate bucket contact point. The contact point position difference analysis is obtained by comparing and calculating the spatial displacement, intrusion angle, and normal vector between candidate nodes point by point, determining the offset of each contact node from its adjacent nodes in the horizontal, longitudinal, and depth directions. Simultaneously, combined with the multibody dynamics model, the maximum intrusion depth, contact force distribution, and overall machine stability changes caused by the contact force at different contact nodes are calculated by simulating the kinematic constraints and actuator capabilities of the bucket in the current loading posture. Through statistical analysis and classification of the position differences of each contact node, a bucket loading characteristic data set of contact point position differences is formed, including node coordinates, normal vector, maximum loading depth, contact force magnitude and direction, estimated center of gravity offset, and stability margin. This characteristic data reflects the interaction differences between the bucket and the material during local loading, providing a quantitative basis for subsequent trajectory optimization.

[0051] Step S422: Based on the bucket loading characteristic data of the contact point differences, perform a force potential energy field characteristic analysis of bucket loading contact specificity to obtain the force potential energy field characteristic data of bucket loading contact specificity. In this embodiment of the invention, after obtaining bucket loading characteristic data showing differences in contact points, a local bucket force model is established, mapping the contact force vector of each contact node onto the bucket structure, and calculating the potential energy contribution generated by the force. The force potential energy field characteristic analysis first applies normal and tangential forces to each node, converting the node forces into a combined potential energy of bucket structure deformation, friction work, and cutting work. Subsequently, the local potential energies of all nodes are superimposed to form a complete bucket force potential energy field, and the local bucket contact force distribution is displayed through a three-dimensional equipotential surface or potential energy field diagram. The analysis incorporates a multibody dynamics model to calculate the bucket's deformation trend, inertial force influence, and impact on overall machine stability under the force potential energy field. By spatially clustering high and low potential energy regions, potentially high-risk contact areas and low-energy-consumption areas are identified, and bucket loading contact-specific force potential energy field characteristic data is generated, including node positions, potential energy values, force vectors, local stability scores, and friction-cutting energy consumption indicators, providing a refined basis for trajectory selection and optimization.

[0052] Step S423: By using the bucket loading characteristic data of the contact point difference and the corresponding bucket loading contact-specific force potential energy field characteristic data, the target load image data at a local moment is processed to select the bucket loading contact trajectory at a local moment, so as to obtain the bucket loading contact trajectory data at a local moment.

[0053] In this embodiment of the invention, a local trajectory candidate space is constructed based on node topology. Using an adjacency graph of the node set as a basis, each possible local trajectory is expressed as a sequence of nodes, and geometric interpolation is performed between nodes using parametric curves (such as cubic splines). The parametric representation of each curve includes a spatial pose sequence (bucket position and pitch angle), an intrusion depth profile, and temporal rules. Each candidate trajectory is spatially subdivided into several micro-segments, and local geometric indices (curvature, normal deviation, intrusion angle, and depth gradient) are calculated for each micro-segment to form a trajectory feature vector. Subsequently, an energy and mechanical evaluation is performed on each micro-segment in the potential energy field: the contact surface of the bucket in that posture is mapped to the potential energy field, the rate of change of potential energy generated by the contact force is calculated, the cutting / friction work is integrated along the trajectory to obtain a micro-segment work estimate, and the overturning moment contribution to the entire vehicle is calculated from the contact force vector through torque mapping. The work done by the micro-segments is accumulated over the entire trajectory time series to obtain a total work estimate, and the minimum stability margin, maximum contact force, and peak pressure or torque required by the actuator at different positions of the trajectory are summarized. Based on the aforementioned physical quantities, a multi-objective evaluation function is constructed. The evaluation function terms include total work done (as a proxy for energy consumption), the mass ratio loaded per unit work done (as a proxy for efficiency), the reciprocal of the minimum stability margin (as a safety penalty), the excess of peak actuator load relative to the capacity grid (as a feasibility constraint penalty), and the geometric deviation of the trajectory from the global path (as a continuity term). Weights are determined by the loading operation demand vector and fused with historical weights according to an exponential smoothing rule to maintain smooth switching. For each candidate trajectory, a constrained numerical optimization solver is used to locally optimize the trajectory parameters: the optimization variables are the minor adjustments to the curve control points and the time-pacing parameters; the objective is to minimize the weighted evaluation function; and the constraints are: no collisions in space (determined through bounding box-level collision detection), limits between the bucket and mechanism joints, the inequality of not exceeding the capacity grid limits for actuators, and the lower limit of the lifting stability margin. The numerical solution employs a sequential quadratic programming method and uses an initial candidate trajectory for a warm start to improve convergence speed; the constraint projection step ensures the rigid satisfaction of geometric and dynamic constraints after each iteration. The output local moment bucket loading contact trajectory data is presented in time series format, providing the bucket pose, intrusion depth, expected contact force, estimated loading increment, actuator reference pressure / torque, and trajectory confidence for each trajectory point. This allows subsequent group prediction and global path correction processes to make continuous decisions and issue control commands based on this structured information.

[0054] Furthermore, step S421 includes the following steps: Based on the image data of the target load at a local time, the shape characteristics of the target load at that local time are analyzed to obtain the shape characteristic data of the target load at that local time. By optimizing the multibody dynamics model of the loader, the contact characteristic nodes of each bucket and the corresponding maximum loading depth characteristic data of the target load at that local time are analyzed. Based on the contact characteristic nodes of each bucket and the corresponding maximum loading depth characteristic data, the bucket loading characteristic analysis of contact point differences is performed to obtain the bucket loading characteristic data of contact point differences.

[0055] In this embodiment of the invention, high-density point cloud fragments are obtained by performing 3D reconstruction processing based on multi-view images and depth point clouds of the target load at local time points. Point cloud preprocessing includes statistical outlier removal, voxel downsampling, and normal estimation. Normal estimation uses neighborhood least squares plane fitting to obtain stable normal directions. The stockpile surface and surrounding stray points are segmented using a clustering algorithm based on Euclidean distance. Subsequently, surface reconstruction is performed on the stockpile surface. The reconstruction method uses implicit surface reconstruction (such as Poisson surface reconstruction or 3D TSDF) to generate closed triangular meshes. Curvature calculation is performed on the triangular meshes. The principal curvature and Gaussian curvature are calculated to identify surface feature points (such as protrusions, notches, and cutting edges). Based on the normals and local curvatures, multiple levels of candidate contact regions are defined: the set of vertices with normals facing the bucket intrusion direction and curvatures within a preset range are marked as first-level candidate nodes; the set of vertices with normals deviating from the intrusion direction but located in the reachable zone are marked as second-level candidate nodes. For each candidate node, a local tangential coordinate system is derived, and a set of candidate direction vectors is generated within this coordinate system to characterize the contact geometry between the bucket and the material under different intrusion angles. Each candidate node is then submitted to the optimized loader multibody dynamics model for feasibility and maximum loading depth analysis. The analysis process is an incremental intrusion simulation: using the candidate intrusion vector as the direction, the bucket geometry model is iteratively advanced along the intrusion vector into the reconstructed surface at small depth steps. In the simulation at each depth step, geometric collision detection is first performed to determine mechanism interference or collision with non-target objects on the ground; secondly, the contact surface area and contact normal distribution formed by the intrusion are calculated; then, the cutting resistance and frictional resistance acting on the bucket's leading edge and sidewalls are estimated based on the material mechanics model, and this resistance vector is projected onto the actuator's line of action to obtain the required driving torque or hydraulic cylinder pressure. Simultaneously, a multibody dynamics model is used to calculate the vehicle center position, wheel load changes, and overturning moment generated by contact force under the bucket posture and loaded mass. The calculated driving force or torque is compared with the boundary values ​​and lifting stability margin thresholds on the actuator capability profile. When the driving force demand reaches or exceeds the actuator's continuous capability boundary, or the overturning moment caused by contact force causes the stability margin to fall below the safety threshold, or mechanism interference occurs, the current intrusion depth is determined to be infeasible, and the maximum loading depth is taken as the depth value of the previous feasible step. Finally, a complete feature set is formed for each candidate contact node, including the node's three-dimensional coordinates (body base coordinate system), normal vector, optimal intrusion vector, maximum loading depth, corresponding required actuator force / pressure, estimated loading volume and mass (obtained by integrating the volume of the swept body of the intrusion trajectory), overturning moment and stability margin value generated by contact force, and node observation confidence (weighted based on point cloud density, normal consistency, and multi-view confidence). This feature set serves as the basic physical quantity input for subsequent contact point difference analysis and contact trajectory optimization.For the candidate contact node set, a spatial adjacency network is first constructed. Adjacency determination is based on the Euclidean distance and the angle between normals between nodes, and efficient nearest neighbor retrieval is achieved using a KD-tree. By performing Delaunay triangulation or Voronoi partitioning on the node set, the contact area is divided into several coverage units, each representing a local material set that may be covered in a single scoop. For each pair of adjacent nodes, a positional difference vector and its component characteristics are calculated: horizontal offset, longitudinal forward offset, tilt angle difference, and depth difference. The positional difference vectors are statistically summarized, and the mean, variance, skewness, and kurtosis of the global positional difference distribution are calculated to quantify the non-uniformity of the contact point distribution; simultaneously, the statistics of the maximum depth distribution are calculated to assess the asymmetry of local loading capacity. Based on the maximum loading depth and corresponding intrusion direction of each node, the overall force distribution of the bucket under simultaneous multi-point contact is solved: the predicted contact force of each node is mapped along the bucket structure as a force distribution vector field. The point load is mapped to the overall force-moment response of the machine using a simplified equivalent beam model of the finite element method or a rigid body force balance method, and then the load center of gravity offset and the roll / pitch moment on the vehicle body are calculated. This generates key loading characteristic indicators, including the expected load imbalance matrix caused by contact point differences, the estimated maximum lateral overturning risk, the expected load variance, and the unevenness of unit work distribution. Furthermore, quantitative suggestions for operation sequence and contact strategy are given based on contact point differences: areas with high elevation gradients are used as the initial cutting zone or segmented scooping zone to form a local contact sequence so that the force and center of gravity movement at each step remain within a stable margin. The output bucket loading feature data set of contact point position differences includes node pair position difference matrix, node group division, estimated loading mass and stability margin curve of each node / group, force distribution vector field, priority scooping sequence and risk score based on position difference, which can be used for group local time correlation optimization and actual execution strategy.

[0056] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0057] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for optimizing the motion trajectory of a skid steer loader during loading operations, characterized in that, Includes the following steps: Step S1: Obtain loading operation demand data; use the multi-source image monitoring equipment built into the skid steer loader to perform multi-source environmental perception monitoring processing on the loading operation area to obtain multi-source environmental perception data; based on the loading operation demand data and the multi-source environmental perception data, perform loading attribute feature identification processing for environmental perception to generate loading attribute environmental perception monitoring data. Step S2: Obtain the design parameters of the skid steer loader; based on the design parameters of the skid steer loader, perform multibody dynamics modeling and operational characteristic parameter identification and optimization processing of the loader to obtain an optimized multibody dynamics model of the loader; Step S3: Establish a panoramic perception coordinate system for loading attributes based on the environmental perception monitoring data of loading attributes to obtain the panoramic perception coordinate system for loading attributes; perform global path analysis of sliding loading based on the panoramic perception coordinate system for loading attributes to generate global path data for sliding loading. Step S4: Extract target load images at each local loading moment based on loading attribute environmental perception monitoring data and skid loader global path data to obtain local target load image data; perform bucket loading contact optimization trajectory analysis based on local target load image data to obtain group local optimized bucket contact trajectory data; optimize the loading motion trajectory of the skid loader global path data based on the group local optimized bucket contact trajectory data to obtain skid loader global optimized motion trajectory data; execute skid loader loading intelligent control operation based on skid loader global optimized motion trajectory data.

2. The method for optimizing the loading trajectory of a skid steer loader according to claim 1, characterized in that, Step S1 includes the following steps: Step S11: Obtain loading operation requirement data; Step S12: Use the multi-source image monitoring equipment built into the skid steer loader to perform multi-source environmental perception monitoring and processing on the loading operation area, generate multi-source environmental perception data, and perform geometric external parameter correction and time sequence synchronization processing on the multi-source environmental perception data through the internal and external participation calibration records of the multi-source image monitoring equipment to generate fused environmental perception monitoring data. Step S13: Perform contour segmentation processing on the environmental perception image based on the image binarization parameters of the fused environmental perception monitoring data to generate environmental perception image segmentation data; Step S14: Based on the loading operation requirement data, perform loading attribute feature identification processing on the environmental perception image segmentation data to generate loading attribute environmental perception monitoring data.

3. The method for optimizing the loading trajectory of a skid steer loader according to claim 1, characterized in that, Step S2 includes the following steps: Step S21: Obtain the design parameters of the skid steer loader; Step S22: Analyze the connection relationship of loading components based on the design parameters of the skid steer loader to obtain the connection relationship data of the loading components, and perform kinematic constraint analysis of the loading components through the connection relationship data to generate the kinematic constraint data of the loading components. Step S23: Analyze the loader actuators according to the skid steer loader design parameters to obtain loader actuator data, and perform actuator capability profile abstraction processing on the loader actuator data to obtain loader actuator capability data; Step S24: Perform loading and lifting characterization analysis based on the skid steer loader design parameters to obtain loading and lifting characterization data; Step S25: Based on the kinematic constraint data of the loading components, the capability data of the loader actuators, and the loading and lifting characterization data, perform multibody dynamics modeling of the loader to obtain the multibody dynamics model of the loader; Step S26: Perform adaptive parameter identification and optimization processing on the multibody dynamics model of the loader to obtain an optimized multibody dynamics model of the loader.

4. The method for optimizing the loading trajectory of a skid steer loader according to claim 3, characterized in that, Step S25 includes the following steps: Step S251: Perform lifting attitude dimension analysis based on the loading lifting characterization data to obtain lifting attitude dimension data; Step S252: Perform linearization characteristic analysis of lifting height and lifting capacity based on the loading lifting characterization data to obtain lifting height-capacity characteristic data; Step S253: Based on the kinematic constraint data of the loading components and the capability data of the loader actuators, perform nonlinear boundary constraints and architecture modeling of the loader's multibody dynamics to obtain the loader's multibody dynamics architecture model; Step S254: Map the lifting posture dimension data and lifting height-capacity characteristic data to the loader multibody dynamics architecture model to perform loading and lifting dynamics characteristic mapping processing to obtain the loader multibody dynamics model.

5. The method for optimizing the loading trajectory of a skid steer loader according to claim 4, characterized in that, Step S252 includes the following steps: Based on the loading and lifting characterization data, dynamic lifting reduction processing is performed to obtain dynamic lifting reduction data of the loader. Then, the lifting height and lifting capacity are linearized and approximated using the dynamic lifting reduction data of the loader to obtain lifting height-capacity characteristic data.

6. The method for optimizing the loading trajectory of a skid steer loader according to claim 1, characterized in that, Step S3 includes the following steps: Step S31: Establish a panoramic perception coordinate system for loading attributes based on the environmental perception monitoring data of loading attributes, so as to obtain the panoramic perception coordinate system for loading attributes; Step S32: Perform a sliding path feasible region analysis based on the panoramic perception coordinate system of the loading attributes to obtain sliding path feasible region data; Step S33: Perform path completion and slip characteristic optimization processing on the slip path feasible region data to obtain slip optimized feasible path data; Step S34: Based on the optimized multibody dynamics model of the loader and the environmental perception monitoring data of loading attributes, perform global path cost feature analysis of skid loading to obtain global path cost feature data of skid loading. The global path cost feature data of skid loading includes path length cost data, skid loading energy consumption cost data, path slope cost data, short-term dynamic obstacle probability penalty cost data, dynamic environmental risk penalty cost data, loading lifting stability cost data, and short-term operation efficiency cost data. Step S35: Establish the mapping relationship between the global path and cost of skid loading based on the global path cost feature data of skid loading, generate a preliminary global path cost model of skid loading, and perform adaptive weight parameter adjustment on the preliminary global path cost model of skid loading through loading operation demand data to obtain the global path cost model of skid loading. Step S36: Transfer the slip optimization feasible path data to the slip loading global path cost model for global path search processing of slip loading, and generate slip loading global path data.

7. The method for optimizing the loading trajectory of a skid steer loader according to claim 6, characterized in that, Step S33 includes the following steps: The loader's slip characteristics were analyzed based on the optimized multibody dynamics model to obtain slip characteristic data. The intermediate nodes of the feasible region data of the skid path are selected for skid path line supplementation processing to obtain skid path line supplementation data. Then, the skid path line supplementation data is subjected to loader skid characteristic line fitting optimization processing based on the loader skid characteristic data to obtain skid optimized feasible path data.

8. The method for optimizing the loading trajectory of a skid steer loader according to claim 1, characterized in that, Step S4 includes the following steps: Step S41: Extract the target loading image at each local loading moment based on the loading attribute environmental perception monitoring data and the skid loading global path data to obtain the target loading image data at the local moment. Step S42: Analyze the bucket loading contact trajectory based on the target load image data at a local time to obtain the bucket loading contact trajectory data at a local time; Step S43: Perform group local time prediction processing of bucket loading contact trajectory based on local time bucket loading contact trajectory data and local time target load image data to obtain bucket loading contact trajectory prediction data; Step S44: Establish a dynamic change state transition matrix and a dynamic change observation matrix of the load based on the bucket loading contact trajectory prediction data. The dynamic change state transition matrix of the load is based on the target load change relationship at each local loading moment derived from the bucket loading contact trajectory prediction data. The dynamic change observation matrix of the load reflects the natural deformation observation value of the target load at each local loading moment. Step S45: Based on the dynamic change state transition matrix of the load and the dynamic change observation matrix of the load, perform group local time correlation optimization processing on the bucket loading contact trajectory data to obtain the group local time optimized bucket contact trajectory data. Step S46: Optimize the loading motion trajectory of the skid loader global path data based on the optimized bucket contact trajectory data at local moments to obtain the skid loader global optimized motion trajectory data; execute the skid loader loading intelligent control operation based on the skid loader global optimized motion trajectory data.

9. The method for optimizing the loading trajectory of a skid steer loader according to claim 8, characterized in that, Step S42 includes the following steps: Step S421: Based on the local time target load image data and the optimized loader multibody dynamics model, perform bucket loading characteristic analysis of contact point differences to obtain bucket loading characteristic data of contact point differences; Step S422: Based on the bucket loading characteristic data of the contact point differences, perform a force potential energy field characteristic analysis of bucket loading contact specificity to obtain the force potential energy field characteristic data of bucket loading contact specificity. Step S423: By using the bucket loading characteristic data of the contact point difference and the corresponding bucket loading contact-specific force potential energy field characteristic data, the target load image data at a local moment is processed to select the bucket loading contact trajectory at a local moment, so as to obtain the bucket loading contact trajectory data at a local moment.

10. The method for optimizing the loading trajectory of a skid steer loader according to claim 9, characterized in that, Step S421 includes the following steps: Based on the image data of the target load at a local time, the shape characteristics of the target load at that local time are analyzed to obtain the shape characteristic data of the target load at that local time. By optimizing the multibody dynamics model of the loader, the contact characteristic nodes of each bucket and the corresponding maximum loading depth characteristic data of the target load at that local time are analyzed. Based on the contact characteristic nodes of each bucket and the corresponding maximum loading depth characteristic data, the bucket loading characteristic analysis of contact point differences is performed to obtain the bucket loading characteristic data of contact point differences.

Images