Multi-robot collaborative operation simulation control method and system based on digital twinning

Abstract

The invention relates to the technical field of digital twinning, in particular to a multi-robot collaborative operation simulation control method and system based on digital twinning. The method comprises the following steps: acquiring robot structure data; constructing a robot three-dimensional twinborn model based on robot structure data; according to the robot three-dimensional twinborn model, warehousing collaborative operation transportation simulation is carried out, and collaborative operation transportation data is obtained; performing path offset anomaly analysis according to the collaborative operation transportation data to obtain path offset data; performing terrain pit detection according to the path offset data to obtain terrain pit data; robot chassis abrasion detection is conducted according to the terrain pothole data, and chassis abrasion data are obtained; performing ultrasonic crack detection according to the chassis wear data to obtain chassis crack data; and collaborative assembly is carried out according to the chassis crack data to obtain collaborative assembly data. The cooperative assembly efficiency and the operation response rate of multiple robots in a complex environment are improved based on the digital twinning technology.

Description

Technical Field

[0001] The present invention relates to the field of digital twin technology, and in particular to a multi-robot collaborative operation simulation control method and system based on digital twin. Background Art

[0002] Traditional multi-robot collaborative operation simulation control typically focuses only on the geometric modeling of the robot's external structure during the model construction phase, ignoring the assembly relationship between components and the internal stress-bearing state, resulting in inaccurate feedback from the simulation model on physical behavior. During the collaborative operation simulation process, it is often based on static path planning algorithms, lacking the ability to dynamically respond to path deviations, inertial interference, and sudden terrain changes during actual operation, making it difficult to effectively identify deviation anomalies and terrain risks during robot operation. In terms of equipment wear detection, it relies only on time periods or number of operations for prediction, failing to dynamically evaluate the actual terrain impact, chassis stress, and structural fatigue, lacking real-time and accuracy. During assembly task execution, traditional systems often use a fixed assembly sequence, lacking quantitative identification and optimization analysis of assembly interference risks, and making it difficult to flexibly adjust the assembly process sequence to improve collaborative efficiency. Furthermore, in multi-robot collaborative simulation, there is a lack of comprehensive analysis and scheduling control of the assembly task priorities of each robot, making it difficult to achieve an efficient and reliable collaborative operation process. The overall system has a low level of intelligence and cannot meet the requirements for high-precision dynamic response to operation paths, structural states, assembly interference, and task scheduling in complex operation scenarios. Summary of the Invention

[0003] Based on this, it is necessary for the present invention to provide a multi-robot collaborative operation simulation control method and system based on digital twins to solve at least one of the above technical problems.

[0004] To achieve the above objectives, a multi-robot collaborative operation simulation control method based on digital twins includes the following steps:

[0005] Step S1: Acquire robot structural data; construct a three-dimensional twin model of the robot based on the robot structural data; perform warehouse collaborative operation and transportation simulation based on the three-dimensional twin model of the robot to obtain collaborative operation and transportation data;

[0006] Step S2: performing path deviation anomaly analysis based on the collaborative operation transport data to obtain path deviation data; performing terrain pothole detection based on the path deviation data to obtain terrain pothole data; performing robot chassis wear detection based on the terrain pothole data to obtain chassis wear data;

[0007] Step S3: performing ultrasonic crack detection based on the chassis wear data to obtain chassis crack data; performing collaborative assembly based on the chassis crack data to obtain collaborative assembly data; performing assembly interference detection based on the collaborative assembly data to obtain assembly interference data;

[0008] Step S4: Optimize the assembly process sequence based on the assembly interference data to obtain the assembly process optimization sequence; perform task priority analysis based on the assembly process optimization sequence to obtain task priority data; simulate the multi-robot collaborative assembly operation based on the task priority data to obtain collaborative assembly operation data.

[0009] This invention, by integrating digital twin technology with multidimensional data fusion and sensing, achieves high-precision modeling and dynamic feedback control of robot structure, path, terrain, component status, assembly interference, and task scheduling throughout the multi-robot collaborative operation process. This addresses the problems of coarse static modeling, delayed state assessment, and inefficient task collaboration inherent in traditional simulation methods. In robot structural modeling, the system not only captures the robot's external geometry but also integrates joint drive data, material properties, and transmission path structure information to construct a three-dimensional twin model with physical response characteristics, improving the model's accuracy in responding to dynamic behavior during actual operation. During the collaborative operation simulation phase, trajectory point cloud data and inertial measurement unit (IMU) data are collected in real time during robot operation. Path deviation detection is achieved by comparing these data with the theoretical path, enhancing the simulation platform's dynamic recognition of path anomalies and effectively detecting path deviations caused by turning inertia, collision disturbances, or navigation errors. In terrain assessment, the system integrates laser radar scanning data from the robot's underbody with track force distribution data to identify subtle ground undulations and the depth and width of potholes. This allows for centimeter-level accuracy in detecting sudden terrain changes, providing fundamental support for chassis damage risk assessment. In chassis wear detection, the cumulative impact fatigue of the wheel support structure and steering knuckle load-bearing areas is dynamically calculated by combining terrain impact intensity, the robot's actual operating speed, and the frequency of continuous collisions. Using metal fatigue formulas and the driving load-frequency accumulation coefficient, the chassis wear state is accurately determined, avoiding inaccurate judgments based solely on operating time or shift cycles. For chassis crack detection, multi-point ultrasonic phased array technology is used to perform high-resolution crack imaging at stress hotspots. By analyzing the time delay, amplitude variation, and velocity distortion of the acoustic echo, the initial crack propagation path is determined at the millimeter level, providing complete structural data for subsequent collaborative assembly. In the assembly phase, a 3D assembly constraint pose library and an overlap analysis method for assembly interference geometry are used to quantitatively identify potential interference risks in the assembly sequence. Assembly logic prioritizing the minimum interference path is reconstructed, improving the flexibility and accuracy of assembly action sequence adjustments. At the task scheduling level, a task priority weight matrix is ​​constructed based on the reachable area distribution, structural redundancy, load capacity, and assembly time window requirements of each robot. This is combined with a multi-objective optimization algorithm to allocate task resources and reconstruct collaborative scheduling paths, improving the convergence speed and work balance of the overall collaborative simulation. Overall, it has achieved high-precision modeling, high-dynamic response analysis, and high-intelligence scheduling control for the entire process of multi-robot assembly tasks in complex structural scenarios, effectively improving collaborative work efficiency and system stability.

[0010] Preferably, this specification also provides a multi-robot collaborative operation simulation control system based on digital twins, which is used to execute the multi-robot collaborative operation simulation control method based on digital twins as described above. The multi-robot collaborative operation simulation control system based on digital twins includes:

[0011] The warehouse collaborative operation and transportation simulation module is used to obtain robot structural data; build a 3D twin model of the robot based on the robot structural data; and perform warehouse collaborative operation and transportation simulation based on the 3D twin model of the robot to obtain collaborative operation and transportation data.

[0012] The robot chassis wear detection module is used to perform path deviation anomaly analysis based on collaborative operation transportation data to obtain path deviation data; perform terrain pothole detection based on path deviation data to obtain terrain pothole data; and perform robot chassis wear detection based on terrain pothole data to obtain chassis wear data;

[0013] The assembly interference detection module is used to perform ultrasonic crack detection based on chassis wear data to obtain chassis crack data; perform collaborative assembly based on chassis crack data to obtain collaborative assembly data; and perform assembly interference detection based on the collaborative assembly data to obtain assembly interference data.

[0014] The collaborative assembly operation simulation module is used to optimize the assembly process sequence based on assembly interference data to obtain the assembly process optimization sequence; perform task priority analysis based on the assembly process optimization sequence to obtain task priority data; and simulate multi-robot collaborative assembly operations based on the task priority data to obtain collaborative assembly operation data.

[0015] The multi-robot collaborative operation simulation control system based on digital twin of the present invention can realize any one of the multi-robot collaborative operation simulation control methods based on digital twin of the present invention, and is used to combine the operation and signal transmission medium between each module to complete the multi-robot collaborative operation simulation control based on digital twin. The internal modules of the system cooperate with each other to improve the collaborative assembly efficiency and operation response rate of multiple robots in complex environments. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Other features, objects and advantages of the present invention will become more apparent upon reading the detailed description of non-limiting embodiments thereof made with reference to the following drawings:

[0017] Figure 1 This is a schematic flow chart of the steps of a multi-robot collaborative operation simulation control method based on digital twins of the present invention;

[0018] Figure 2 Detailed step flow diagram of step S3 in the present invention;

[0019] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings and in conjunction with the embodiments. DETAILED DESCRIPTION

[0020] The following is a clear and complete description of the technical method of the present invention in conjunction with the accompanying drawings. Obviously, the embodiments described are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without making any creative work are within the scope of protection of the present invention.

[0021] In addition, the accompanying drawings are merely schematic illustrations of the present invention and are not necessarily drawn to scale. Identical reference numerals in the figures denote identical or similar parts, and thus repetitive descriptions thereof will be omitted. Some of the block diagrams shown in the accompanying drawings are functional entities that do not necessarily correspond to physically or logically separate entities. These functional entities may be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor and / or microcontroller approaches.

[0022] It should be understood that although the terms "first," "second," and the like may be used herein to describe various elements, these elements should not be limited by these terms. These terms are used solely to distinguish one element from another. For example, a first element may be referred to as a second element, and similarly, a second element may be referred to as a first element, without departing from the scope of the exemplary embodiments. The term "and / or" as used herein includes any and all combinations of one or more of the listed associated items.

[0023] To achieve this, please refer to Figures 1 to 2 The present invention provides a multi-robot collaborative operation simulation control method based on digital twins, the method comprising the following steps:

[0024] Step S1: Acquire robot structural data; construct a three-dimensional twin model of the robot based on the robot structural data; perform warehouse collaborative operation and transportation simulation based on the three-dimensional twin model of the robot to obtain collaborative operation and transportation data;

[0025] In this embodiment, in the process of obtaining the robot structure data, a high-precision three-dimensional laser scanning device is used to perform a structural scan on each robot in the multi-robot system, and the scanning resolution must be controlled within 0.1mm to ensure that each structural detail can be fully extracted. After the scanning is completed, the laser point cloud data is imported into the three-dimensional reconstruction tool MeshLab for noise removal and edge repair, and then the solid geometry conversion is performed using a CAD reconstruction plug-in (such as Geomagic Design X) to obtain complete structural data of structural components such as the chassis, transmission system, assembly arm, wheel set and motor drive unit. The structural data must be exported in accordance with the ISO 10303 standard (STEP format) to ensure that the structural data is parseable in subsequent platforms. Based on the above structural data, an industrial-grade digital twin platform (such as Siemens Tecnomatix or Dassault 3D EXPERIENCE) is used to build a three-dimensional twin model of the robot. During the construction process, the kinematic constraint relationship between each structural component needs to be set. For example, the RP pair (rotation-translation pair) is used between the assembly arm and the chassis, and the R pair is used between the drive wheel and the chassis to ensure that the virtual model has motion response characteristics consistent with the physical device. Dynamic parameters must be set according to the equipment manufacturer's instructions, including the maximum torque of each joint, the joint limit angle, and the acceleration threshold. The torque setting is based on 10Nm, the joint limit angle is set to ±180 degrees, and the maximum joint linear velocity does not exceed 1.5m / s. After construction is completed, the three-dimensional data of the warehouse scene is introduced based on the above-mentioned digital twin model. This scene data is acquired in the form of a point cloud and then constructed into a virtual work warehouse environment with elements such as obstacles, slopes, and platforms using terrain modeling tools such as ReCap and AutoCAD Civil 3D. Then, a path control engine (such as ROS Navigation Stack) is deployed on the platform. By setting fixed transportation start and end points (for example, the starting point coordinates X1=0, Y1=0, Z1=0, and the end point X2=15, Y2=5, Z2=0) and an obstacle avoidance strategy (using the Dijkstra algorithm for path calculation), a multi-robot collaborative operation and transportation simulation is performed. The simulation process records the robot's position changes (recorded every 100ms), waypoint sequence, acceleration, wheel speed and other transport-related data, and finally forms the collaborative operation transport data and exports it in CSV format for use in the next stage.

[0026] Step S2: performing path deviation anomaly analysis based on the collaborative operation transport data to obtain path deviation data; performing terrain pothole detection based on the path deviation data to obtain terrain pothole data; performing robot chassis wear detection based on the terrain pothole data to obtain chassis wear data;

[0027] In this embodiment, after obtaining the collaborative operation transport data, it is necessary to perform path deviation anomaly analysis on the simulated operation trajectory of each robot. The Euclidean distance of each robot's simulated position sequence is compared with the theoretical trajectory point by point, and the deviation judgment threshold is set to 30mm. When the deviation value exceeds this threshold at any moment, it is determined to be a path deviation anomaly point. This analysis process is processed using the Scipy space vector library in Python. The path deviation rate and deviation peak of each robot are recorded in a tabular form to form path deviation data. Based on the above path deviation data, the geographic coordinates of the deviation anomaly points are analyzed and compared with the original terrain data to detect terrain potholes. A terrain height grid map (with 0.2mx 0.2m as the unit grid) is used to analyze the regional height. When the maximum height difference in the same grid exceeds 50mm, it is determined to be a terrain pothole area. Matplotlib is used to visualize and form a terrain pothole distribution map, which is used as the terrain pothole data input into the subsequent stage. Based on the terrain pothole data, the chassis wear detection is performed in combination with the wheel force data of the robot passing through the corresponding area in the simulation. The impact force F during each pothole is calculated as: F = m × a, where m is the local mass of the wheel, set to 5 kg, and a is the acceleration recorded in the simulation, with a maximum value of 0.6 g. If the number of consecutive impacts n exceeds 30 and the impact force F exceeds 50 N, the chassis structure is considered to be at risk of wear. Test results are output as wheel number, position, peak force, and number of impacts, forming chassis wear data.

[0028] Step S3: performing ultrasonic crack detection based on the chassis wear data to obtain chassis crack data; performing collaborative assembly based on the chassis crack data to obtain collaborative assembly data; performing assembly interference detection based on the collaborative assembly data to obtain assembly interference data;

[0029] In this embodiment, according to the chassis wear data, an industrial ultrasonic flaw detection module is used to perform crack detection. The coordinates of the wear area marked in the simulation are imported into the three-dimensional twin platform and located correspondingly to the virtual robot chassis part. The CIVA ultrasonic simulation tool is used, the probe frequency is set to 5MHz, the pulse width is 0.2μs, the scanning mode is scanning + angle reflection method, and the incident angle is set to 45°. In the simulation environment, the wear area is scanned and simulated by adjusting the simulated probe path to obtain the echo signal. When the peak amplitude of the reflected signal exceeds the background noise + 12dB and there is a double peak phenomenon (first reflection and crack reflection), it is identified as a crack area. The final output crack length, depth and position constitute the chassis crack data. Combined with the crack location and structural model, a multi-robot collaborative assembly operation simulation is performed. The collaborative goal is set to dismantle the crack area components and replace them with new structural parts. Tecnomatix Process Simulate is used to perform collaborative disassembly and assembly tasks. By specifying two robots to perform disassembly and conveying operations respectively, the path is generated by setting a virtual collision boundary. During the operation, the action sequence, execution time, collision data, etc. are output to form collaborative assembly data. Assembly interference detection is performed based on collaborative assembly data. The 3D interference analysis module is called within Process Simulate to detect spatial overlap between assembly actions frame by frame. The decision rule is: interference is flagged when the gap between any two components is less than 1.0 mm or the 3D bounding box overlap exceeds 10%. Interference data is statistically output by time series, part number, and interference type (linear / planar) to generate assembly interference data.

[0030] Step S4: Optimize the assembly process sequence based on the assembly interference data to obtain the assembly process optimization sequence; perform task priority analysis based on the assembly process optimization sequence to obtain task priority data; simulate the multi-robot collaborative assembly operation based on the task priority data to obtain collaborative assembly operation data.

[0031] In this embodiment, assembly process sequence optimization is performed based on assembly interference data. The original assembly sequence is represented as a directed graph structure, with nodes representing assembly actions and edges representing dependencies. Non-critical paths are rearranged through topological sorting, and the order of conflicting nodes is adjusted based on the interference data. A genetic algorithm-based sorting optimization method is used, with the fitness function set as the total number of interferences + the total length of path moves, the population size set to 50, the maximum number of iterations set to 200, the crossover rate set to 0.8, and the mutation rate set to 0.05. Ultimately, an assembly process optimization sequence with no interference and the shortest total path is obtained. Based on the assembly process optimization sequence, all tasks are prioritized. Each assembly task is quantitatively scored based on the critical path dependency depth, the robot idle time matching degree, and the degree of workspace conflict. The scoring formula is: P = 0.5 × D + 0.3 × T + 0.2 × S, where D is the dependency depth (higher layers have higher priority), T is the idle matching degree (more idle robots have higher scores), and S is the spatial conflict (fewer conflicts have higher scores). Task priority data is imported into the simulation platform in JSON format for simulation scheduling. Based on this priority data, a multi-robot collaborative assembly operation simulation was conducted. The optimized sequence and priority scheduling strategy were loaded into Process Simulate. Using an event-driven simulation engine, each robot was simulated for 30 minutes, with a sampling interval of 100ms. The initial positions of all robots were uniformly set to the work area surrounding the assembly platform. During the simulation, information such as each robot's task sequence, conflict avoidance actions, and execution efficiency was recorded. Ultimately, collaborative assembly operation data was generated and output as a visual process diagram and structure change log for analysis.

[0032] Preferably, step S1 is specifically as follows:

[0033] Step S11: Obtain robot structure data;

[0034] In this embodiment, based on the robot manufacturing standards and design drawings, detailed structural parameters related to the robot body structure are collected. Specifically, the position parameters of each motion joint of the robot (such as DH parameters: connecting rod length a, connecting rod offset d, connecting rod torsion angle α and joint variable θ), the installation position of the drive unit, the range of motion, the maximum torque (unit: Nm), the load capacity (unit: kg), etc.; the chassis structure parameters need to be collected, including wheeled or tracked drive configuration, chassis dimensions (length, width and height, unit: mm), minimum turning radius (unit: m), ground pressure (unit: kPa), etc. The structural data needs to be measured by a laser measurement system (such as Leica TS60) or a high-precision three-dimensional coordinate measuring machine (such as ZEISS PRISMO). The measurement error must be controlled within ±0.01mm to ensure data accuracy. Material information (such as aluminum alloy 6061-T6, stainless steel 304, carbon fiber composite material, etc.) strength, elastic modulus (unit: GPa), density (unit: g / cm 3 ) should also be uniformly coded and obtained for use in subsequent analysis processes for mass distribution and deformation simulation. The above data is uniformly stored in the structural data table, which uses fields to mark component numbers, position coordinate systems, dimension values, connection relationship identifiers, etc.

[0035] Step S12: performing component space modeling based on the robot structure data to obtain the robot's three-dimensional component structure;

[0036] In this embodiment, the modeling is performed in a CAD modeling environment, such as SolidWorks 2023 or Siemens NX. First, according to the geometric parameters obtained in step S11, the geometric modeling tools such as "extrude", "rotate", "sweep" and "loft" are used to build each segment of the robot arm one by one, including the base, rotating arm, connecting rod, joint housing, end effector bracket, etc. The modeling of each component should be based on the original dimensions (length, diameter, wall thickness, etc.) calibrated in the structural drawing, for example, the connecting rod length L = 520mm, the joint housing inner diameter D = 140mm, the thickness t = 5mm, etc. The modeling of each structural component requires defining a local coordinate system, calibrating the origin position and orientation, and embedding connection interface data, such as the joint interface male and female key value, flange diameter, screw hole spacing and other parameters. After the modeling is completed, the built-in geometric consistency verification tool of the CAD tool needs to be used to verify the continuity of the boundary surface and the integrity of the solid volume. The constructed three-dimensional component structure needs to be exported as a standard format file (STEP or IGES) and the structure number and version number must be bound in the PDM (product data management) system.

[0037] Step S13: reconstructing the component assembly relationship according to the robot's three-dimensional component structure to obtain component assembly relationship data;

[0038] In this embodiment, the assembly is carried out in an assembly environment, supported by the SolidWorks Assembly module or the assembly design module in PTC Creo. According to the three-dimensional component structure constructed in step S12, the assembly operation is performed in the order of numbers in the design drawing (such as starting from the chassis, to the arm body, joint, and end). The assembly of each component depends on its connection method. If it is a bolt connection, the input parameters include the bolt specification (such as M6×30), the tightening torque (unit: N·m), if it is a plug-in connection, the insertion depth (unit: mm) and the fit tolerance (such as H7 / k6) are defined; if it is a welded connection, the weld type (such as V-groove, butt weld) and weld size (unit: mm) need to be specified. The assembly relationship data is represented by a "constraint matrix" between parts, which includes position constraints (such as coaxial, facing), degree of freedom restrictions (such as rotational freedom constraint is θ=±180°) and interference check marks. After assembly is completed, the interference detection tool is used to verify the physical conflict of all interfaces, requiring the interference area to be <0.01mm 2 The final assembly relationship data is exported in XML format, where the node information includes fields such as component ID, assembly constraint type, coordinate transformation matrix, and assembly sequence index.

[0039] Step S14: constructing a three-dimensional twin model of the robot based on the component assembly relationship data and the three-dimensional component structure of the robot to obtain a three-dimensional twin model of the robot;

[0040] In this embodiment, a three-dimensional modeling engine (such as Unity 2022.3LTS or Blender 3.6) is used to reconstruct the robot twin model. The component structure (STEP file) is used as the geometric basis and imported into the modeling engine. Combined with the assembly relationship XML file in step S13, the relative position and assembly posture of each component are parsed. The assembly constraint data is loaded through the scene script and the entity hierarchy structure (such as base_link→link1→joint1→link2) is constructed. In this process, each level component needs to define its transformation matrix (4×4SE(3) matrix) in the local coordinate system, including translation (unit: mm) and rotation (unit: °) parameters. All movable joints should be embedded in the real-time drive interface. By adding a joint controller script (such as a custom JointController class), the joint angle range (such as θ∈[-90°,+90°]) is defined and the control signal input channel (such as CAN bus ID address) is bound. Each component in the twin model needs to specify material properties (such as diffuse reflectance, glossiness, color RGB value) to achieve visual simulation. The final twin model needs to be verified for consistency in a simulation environment to ensure that the deviation between the actual measured size and the virtual model is controlled within ±0.5mm, and the time error between the simulated motion and the actual behavior does not exceed 50ms.

[0041] Step S15: Perform warehouse collaborative operation and transportation simulation based on the robot's three-dimensional twin model to obtain collaborative operation and transportation data.

[0042] In this embodiment, it is carried out in an industrial-grade virtual simulation platform (such as FlexSim 2023 or Gazebo 11). First, the robot 3D twin model is imported into the simulation environment, and the storage environment information is set, including the site grid coordinates (each grid size is 1000mm×1000mm), shelf distribution coordinate points, ground friction coefficient (μ=0.6), cargo size (for example, 600mm×400mm×300mm) and weight (such as 12kg). Set the multi-robot task distribution strategy, for example, based on the Round-Robin strategy for task allocation, each robot is assigned a starting position, a target shelf, and a path node coordinate sequence. The path planning uses the A* algorithm. During the planning process, the obstacle area cost is defined as 1000, the pass area cost is defined as 1, and the Euclidean distance is used as the heuristic function. The robot movement speed is simulated during transportation (the maximum linear speed is 1.2m / s and the acceleration is 0.5m / s 2 ), turning radius (minimum 1.2m), and load-to-speed coefficient (weight increases by 1kg, speed decreases by 0.02m / s). During the simulation, the robot path point coordinates (x, y, t), number of obstacle avoidances, total path length, task completion time, offset value, etc. are collected and output as a structured data table. The collaborative operation transportation data is exported as a CSV file with fields such as robot number, task number, number of path points, total path length (unit: m), cumulative offset (unit: mm), running time (unit: s), etc., for subsequent anomaly detection and analysis.

[0043] Preferably, step S15 is specifically as follows:

[0044] Step S151: uploading the robot 3D twin model to the robot collaborative operation simulation platform;

[0045] In this embodiment, the constructed three-dimensional twin model of the robot is uploaded to the multi-robot collaborative operation simulation platform, and a structured interface protocol is required to complete data docking. The three-dimensional twin model is exported in glTF (GL Transmission Format) format, which supports embedded mesh data, texture data, and joint level data. The upload operation calls the model registration interface of the simulation platform through the POST method in the HTTP protocol. The interface URL is / robot / uploadModel, and the required fields include: model ID, file path, model resolution (unit mm, value 0.1-1.0), number of motion joints (integer, range 1-12), and model coordinate system definition (right-handed system, Z-axis upward). The upload process needs to verify the integrity of the file, and use the MD5 digest algorithm to generate a 32-bit hash value for the model file before and after uploading. If it is inconsistent, the upload is terminated and prompted to re-export. After the upload is completed, the simulation platform will return a unique model reference ID (modelRefID), which must be referenced in the subsequent simulation parameter configuration.

[0046] Step S152: Set the robot's maximum speed to 0.5m / s-2.0m / s and the acceleration range to 0.1m / s 2 -1.0m / s 2 , minimum turning radius 0.3m-1.0m;

[0047] In this embodiment, when setting the physical motion parameters of the robot, the physical boundary values ​​are directly input through the platform parameter configuration module. The maximum driving speed is set to 0.5m / s to 2.0m / s, with a step of 0.1m / s, forming a set of discrete speed levels, a total of 16 speed values, to facilitate the evaluation of path selection at different speeds during the simulation process. The acceleration range is set to 0.1m / s 2 Up to 1.0m / s 2 , step size is 0.1m / s 2 , a total of 10 levels of acceleration parameters. Each acceleration parameter needs to be calculated in combination with the robot's load mass and the rated power of the drive unit to see if it is achievable. The power formula P = F × v is used to derive the upper limit of acceleration. The driving force used is F = m × a. The mass m is the robot plus the current load mass, and the total mass range is set to 30kg-50kg. The minimum turning radius is determined by the robot's chassis length and wheelbase. The setting range is 0.3m to 1.0m, with a step of 0.1m. The corresponding vehicle body physical characteristics such as wheelbase length of 0.45m and wheel spacing of 0.55m need to match the actual data. All parameters are entered through the parameter setting module in the platform configuration file, and the parameter file is saved in YAML format.

[0048] Step S153: Setting the terrain slope to 0°-5°, the ground friction coefficient to 0.2-0.8, and the load variation to 5kg-20kg;

[0049] In this embodiment, terrain slope parameters are input through the scenario building module in the simulation platform. The slope setting range is 0° to 5°, in 1° increments. Five different simulation scenarios are named slope_0, slope_1, ..., slope_5 and loaded as different test maps. The ground friction coefficient input range is 0.2 to 0.8, in 0.1 increments. The static friction factor between the ground and the robot tires is set according to the Coulomb friction model μ = Ff / N. The friction coefficient is set through the simulation platform's material editing tool. The tire material is set to hard rubber (μ reference value 0.6), and the ground is set to cement (μ reference value 0.5) or steel plate (μ reference value 0.3). Material combinations are set by specifying them in the interface. The load variation is set to 5kg to 20kg, in 5kg increments. The load is represented by a standard cubic load model attached to the back of the robot. The model dimensions are 300mm × 300mm × 300mm. The density is set to 18.5kg / m according to the mass formula ρ = m / V. 3 、37kg / m 3 、55.6kg / m 3 , 74kg / m 3 All data is written into the platform scenario parameter configuration file: terrain:{slope:3,friction:0.4,loadMass:15} and bound to the corresponding scenario number.

[0050] Step S154: Set the scheduling instruction delay to 10ms-200ms and the path planning refresh period to 1s-10s;

[0051] In this embodiment, the scheduling instruction delay parameter is set via the message queue delay injection module in the simulation control module. The delay range is set to 10ms to 200ms, in 10ms increments. The message queue controller introduces a random delay value before each robot receives a control instruction. Specifically, this is achieved by storing pending instructions in a FIFO buffer queue and then pushing the instructions to the execution module according to the set delay time. The path planning refresh cycle is set to 1s to 10s, in 1s increments. The corresponding control logic is: after the refresh cycle time expires, the A path planning algorithm module is re-invoked to update the path node sequence. Path planning uses the A algorithm, with a map grid resolution of 0.1m×0.1m. Obstacle information is extracted from 3D environment modeling data and converted into a Boolean grid with a grid size of 128×128. Each round of planning requires completing node traversal and cost function calculation within the refresh cycle. The cost function f(n) = g(n) + h(n), where g(n) is the current path cost and h(n) is the heuristic function (set as the Euclidean distance). All scheduling and path parameters are written to the scheduling control module configuration file. The JSON structure example is as follows: {delayMs:50,pathRefreshSec:5}.

[0052] Step S155: Run the robot behavior control module and the path interaction module in the collaborative operation simulation platform to obtain collaborative operation transportation data.

[0053] In this embodiment, when running the robot behavior control module and path interaction module in the simulation platform, the following two independent threads are launched in sequence: Thread 1, robot_behavior_controller, is responsible for reading the control instructions output by the scheduling module and performing real-time position updates based on the set acceleration, maximum speed, and path node sequence; Thread 2, path_interaction_engine, is responsible for parsing the path planning results, converting them into executable action instruction sequences (such as "move forward 0.5m, turn left 30°"), and sending them to the behavior control thread. The behavior control thread has a control cycle of 100ms and uses a PID controller to correct the speed and direction of the machine position feedback for each round, with a proportional coefficient Kp = 1.2, an integral coefficient Ki = 0.01, and a differential coefficient Kd = 0.05. The control target is the center position coordinate of the path node, and the control error is less than 0.05m. The collaborative operation transportation data is stored in real time by the data recording module, including each robot's spatial position (x, y, θ) at each moment, control instructions, path node sequence number, transportation task status (loading, transporting, unloading), and interaction event log. The data is saved as a CSV file, with each line representing a 1-second snapshot of the state. For example, it contains: timestamp, robotID, x, y, theta, cmd, pathIndex, and taskStatus. All data is used for subsequent behavior evaluation and path performance analysis.

[0054] Preferably, the path deviation anomaly analysis in step S2 is specifically as follows:

[0055] Extract robot operation trajectory based on collaborative operation transportation data;

[0056] In this embodiment, a complete simulation environment including multiple operating robots (such as AGV transport robots, crawler loading robots) and virtual working environments (such as mining road networks or port loading and unloading sites) is constructed in the digital twin simulation platform. By accessing the TF transformation broadcaster and Odometry module in the ROS system, the coordinate position data, azimuth (Yaw), timestamp information and speed information of each robot are collected. In the actual data extraction process, the position coordinates (x, y, z) and timestamps collected by the robot every 0.1 seconds are structured and stored to form time series trajectory data. If the robot has an IMU component, its inertial posture data (acceleration, angular velocity) can be further extracted as trajectory auxiliary data, and the recording period is also 0.1 seconds. The running trajectory of each robot is output in the form of a "timestamp-position-azimuth" triplet for use in subsequent trajectory analysis steps.

[0057] Perform time alignment processing on the trajectory point sequence based on the robot's running trajectory to obtain time series standard trajectory data;

[0058] In this embodiment, after obtaining the trajectory data of multiple robots, in order to ensure the consistency of trajectory curvature and offset calculation, all trajectories need to be time synchronized. The time alignment process uses a linear interpolation method, sets a unified standard time interval Δt to 0.1 seconds, and resamples the trajectory points according to a unified start time. The specific operation process is: using the time axis of the first robot trajectory as a reference, interpolate the missing time points in the trajectories of other robots in turn, and the interpolation content is the (x, y) coordinates and Yaw angles corresponding to the missing time points; use the linear interpolation method to fill in the missing data points. During the processing, the trajectory of each robot is uniformly resampled into an equally spaced time series to obtain a standardized trajectory sequence. The data structure is: [t0, t1, ..., tn] corresponds to the (x, y, θ) data at each moment, ensuring that the trajectories of all robots have the same time reference.

[0059] Calculate trajectory curvature based on time-series standard trajectory data;

[0060] In this embodiment, the sliding window method is used to analyze the spatial steering changes of adjacent trajectory points for the standard trajectory that has completed time alignment. The sliding window width is set to 3 time points, and the trajectory segment composed of three consecutive points is used as the unit for analysis. In each sliding window, the three position points A (x1, y1), B (x2, y2), and C (x3, y3) are recorded at the front, middle, and back, and the trajectory segment steering structure is constructed based on these three points. The local steering feature of the trajectory is extracted by calculating the angle change between point B and points A and C. To ensure continuity and robustness, a sliding window advancement operation is performed every 0.1 seconds, and a series of curvature change values ​​corresponding to time are output in units of 1 / m. During the calculation process, the minimum length of the sliding window is limited to 0.3m and the maximum length is limited to 0.6m to adapt to the trajectory feature density in different operating scenarios.

[0061] Identify track curves according to track curvature and obtain track curve data;

[0062] In this embodiment, curvature data is used to identify trajectory segments whose curvature continuously exceeds a set threshold as "curve" areas. The curve identification threshold is set to 0.15 / m, that is, when the curvature of three or more consecutive trajectory points is greater than 0.15 / m and the direction is continuous (the curvature sign is consistent), the segment is marked as the starting point of the curve, and the start time is recorded; when the curvature is continuously lower than 0.1 / m, it is marked as the end point of the curve, and the end time is recorded. By screening and merging these trajectory segments that meet the conditions, multiple curve segment indexes are formed, and the corresponding trajectory coordinate information is extracted to form a trajectory curve data set. Each curve data is saved in the format of "curve number, start and end time, start and end coordinates, maximum curvature, average curvature" for subsequent offset analysis.

[0063] Performing trajectory deviation detection on the trajectory curve data according to preset standard trajectory curve data to obtain trajectory deviation data;

[0064] In this embodiment, a standard trajectory template for the curve section is established. The standard template is a trajectory path generated in advance in the digital twin system to simulate ideal operating conditions. The standard trajectory curve is compared point by point with the current curve trajectory of the robot to be tested, and the coordinate error corresponding to the same time point is calculated. The reference direction of the lateral offset detection is set to the normal direction of the current trajectory point (perpendicular to the tangent direction of the trajectory), and the coordinate projection method is used to calculate the offset distance between the current point and the standard trajectory in the normal direction. The offset threshold is set to 0.25m, and the points exceeding this value are recorded as "offset points". Each offset point constitutes an offset area data, and its starting point, end point, maximum offset value and average offset value are extracted to form a complete trajectory offset data set.

[0065] Calculate the lateral offset amplitude based on the trajectory offset data; calculate the longitudinal offset amplitude based on the trajectory offset data;

[0066] In this embodiment, each section of offset area data is traversed and its lateral offset distance is statistically analyzed. The calculation of the lateral offset amplitude is based on the lateral projection distance of each time point in the offset segment. The lateral error values ​​of all offset points are extracted (positive values ​​indicate leftward offset, negative values ​​indicate rightward offset), the maximum value, minimum value and average value are recorded, and output in the format of "offset segment number, maximum lateral offset, average lateral offset". The maximum lateral offset limit value is set to 0.6m to mark abnormal offset segments. The amplitude data is used as a lateral stability evaluation indicator for the next step. In the same offset segment, the longitudinal offset along the tangent direction of the trajectory is calculated. The tangential distance error between the standard trajectory point and the actual trajectory point is used as the longitudinal offset value. The longitudinal offset calculation window length is set to 3 time points. For each time point, the tangential difference between it and the same time point of the standard trajectory is calculated to obtain a continuous longitudinal offset sequence. The maximum value and average value in the sequence are counted to represent the longitudinal maximum error and stability error, respectively. If the longitudinal offset amplitude exceeds 0.8m, the segment is marked as a potential inertial interference segment.

[0067] Calculate the lateral deviation volatility based on the lateral deviation amplitude;

[0068] In this embodiment, the lateral offset value sequence of all offset segments is extracted and a second-order difference operation is performed on the sequence to obtain the change amplitude sequence of the offset value between each consecutive time point. The standard deviation of the sequence is calculated and defined as the lateral offset volatility of the segment. The volatility evaluation period is set to 1 second, and the change amplitude of the offset value within each second is evaluated separately. The volatility exceeds 0.3m / s. 2 The segments are recorded as "high-frequency horizontal fluctuation segments". Each segment of volatility is classified and managed based on the time index.

[0069] Perform inertial offset detection based on the longitudinal offset amplitude to obtain inertial offset abnormal data;

[0070] In this embodiment, inertial analysis is performed by combining longitudinal offset amplitude and velocity information. The robot's inertial offset assessment criteria are set: when the longitudinal offset exceeds 0.8m and the speed is simultaneously greater than 1.2m / s, it is considered inertial-driven offset behavior. During the detection process, combined with the acceleration data from the IMU, if the acceleration direction is consistent with the offset direction, it is further confirmed as an inertial offset anomaly. The time periods that meet the above conditions are numbered, and their start and end times, maximum longitudinal offset, and corresponding speed values ​​are extracted and recorded as inertial offset anomaly segments to form an inertial offset anomaly data table.

[0071] The path deviation data are obtained by integrating the lateral deviation fluctuation rate and the inertial deviation anomaly data.

[0072] In this embodiment, the time axis is used as the main index, and the "horizontal volatility value" and "whether there is an inertia abnormality mark" fields are horizontally connected. 2 Time periods with abnormal inertial offset are recorded as high-risk segments for path offset. The start and end times, corresponding robot IDs, start and end coordinates, and maximum lateral and longitudinal offset values ​​for these high-risk segments are structured and recorded. This outputs a path offset dataset to support analysis and use in subsequent path planning or control feedback modules.

[0073] Preferably, the terrain pothole detection in step S2 is specifically as follows:

[0074] Identify the deviation mutation point based on the path deviation data;

[0075] In this embodiment, in the digital twin environment, a standard path data set for the work site is pre-constructed and represented in the form of a two-dimensional coordinate sequence. The path offset data is collected in real time by the working robot at a frequency of 1 Hz during the operation process, and the acquisition frequency is determined by the periodic parameters set in the robot control system. The path offset value is calculated as follows: the Euclidean distance between the robot's current position coordinates (x_actual, y_actual) and the corresponding standard path point coordinates (x_ref, y_ref) is used as the offset value to form a continuous offset time series. In order to identify the mutation point, a sliding window difference algorithm is used, the window length is set to 5 sampling points, the difference between the maximum offset value and the minimum offset value in each window is calculated, and the mutation threshold Δthresh = 0.35m is set. When the offset difference of a window exceeds the threshold, the timestamp corresponding to the center point of the window is marked as the offset mutation point. To ensure the accuracy of the mutation point, data with a distance between adjacent mutation points less than 2 seconds need to be eliminated to avoid misjudgment due to short-term jitter.

[0076] Calculate the terrain slope based on the offset mutation point to obtain slope data;

[0077] In this embodiment, a distance of 1.0m is extended before and after each identified offset mutation point, and the elevation data (z coordinate) collected by the robot within this range is extracted. An elevation point is collected every 0.1m to form a data window with a length of 21. Based on the elevation points in the window, the slope value is calculated using the first-order difference, that is, Δz / Δx calculation is performed on continuous points, and Δx is a fixed value of 0.1m. The average slope value of the terrain section is obtained by averaging the absolute values ​​of all slopes in the window. An average slope value is generated for each offset mutation point to form a discrete slope data set. If the standard deviation of the slope value is greater than 0.05, it means that the slope of the section changes dramatically and is recorded as a high slope area. Throughout the process, all coordinate data must be standardized using the coordinates provided by the RTK-GPS system, the unified coordinate reference is the WGS-84 system, and the projection is converted into the UTM coordinate system for Euclidean calculations.

[0078] Calculate terrain height difference based on slope data;

[0079] In this embodiment, based on the average slope value obtained in the previous stage, each mutation point is taken as the center, and the elevation extreme values ​​within the range of 1.0m are taken forward and backward, which are recorded as z_max and z_min respectively. The terrain height difference value is defined as z_max-z_min. In order to accurately extract extreme values, the number of sampling points is required to be no less than 15 (sampling at a spacing of 0.1m is 21 points), and abnormal peak points are excluded by setting upper and lower thresholds. The threshold is set to the slope average value ±3×slope standard deviation. The calculated height difference values ​​are stored in the terrain feature data set for subsequent area screening. If the difference value is greater than 0.2m, the terrain changes in the area are large, and the next step of suspected area identification is entered.

[0080] Identify suspected pothole areas based on terrain height differences;

[0081] In this embodiment, the terrain height difference values ​​generated above are screened, and the pothole recognition threshold is set to Δh_thres = 0.3m. Any area that satisfies the condition that the difference value is greater than the threshold and the length of the difference position interval is less than 1.5m is defined as a suspected pothole area. The length of the pothole area is calculated by the lateral distance between the elevation extremes. If it exceeds 1.5m, it is classified as a slope or step rather than a pothole. To further prevent misidentification, the situation where there are more than two adjacent areas with continuous differences in height is excluded, that is, if the distance between the centers of two potholes is less than 1m, only the one with a larger depth is retained. The output of this step is a plurality of segment coordinate intervals (start and end UTM coordinates), and the lower left corner is uniformly used as the reference point and saved in the list of suspected pothole areas.

[0082] Collect robot posture angle data based on suspected pothole areas;

[0083] In this embodiment, the working robot posture acquisition module is started in the suspected pothole area, the acquisition frequency is set to 50Hz, and the acquisition period is the time required for the robot to pass through the area. The acquisition parameters are the pitch angle (pitch), roll angle (roll) and yaw angle (yaw) output by the robot's IMU module. In order to ensure the accuracy of the posture angle data, Kalman filter fusion processing of the three-axis gyroscope and accelerometer data is required, and the filter window is set to 0.2 seconds. The data recording time range is 0.5 seconds before entering the pothole area to 0.5 seconds after leaving, ensuring that the complete dynamic response is included. Each suspected area corresponds to a set of posture angle time series, which are used for subsequent impact assessment.

[0084] Evaluate terrain impact intensity based on robot posture angle data;

[0085] In this embodiment, the rate of change of the posture angle when the robot passes through a suspected pothole area is used as the basis for evaluating the impact intensity. The time derivatives of the pitch angle and roll angle are calculated separately (differential form: Δθ / Δt, Δt=0.02s), and the maximum absolute value is taken as the impact index θ_rate_max of the area. If θ_rate_max is greater than 10° / s, it is marked as a "strong impact area". The impact intensity is divided into three levels: θ_rate_max<5° / s is low impact, 5° / s≤θ_rate_max<10° / s is medium impact, and θ_rate_max≥10° / s is high impact. This classification is used by the control system to generate dynamic response control instructions and adjust the robot's collaborative path and speed according to the impact intensity level.

[0086] The suspected pothole areas are screened according to the terrain impact intensity to obtain the terrain pothole data.

[0087] In this embodiment, in all suspected pothole areas, areas with impact intensity levels of medium impact and above are screened out as the final terrain pothole data. At the same time, it is required that the terrain drop value Δh>0.3m and the average slope value in the area are greater than 5° and not more than 30° (to prevent misjudgment as a steep slope). If any of these items is not met, the area is removed from the pothole data. The final screening results are stored in the form of regional center coordinates and start and end UTM coordinates, and the save format is a GeoJSON file. The fields include: area number, center coordinates, attitude impact intensity value, slope value, drop value, and impact level. This data will be used by the subsequent multi-robot simulation control module to adjust the path planning and task allocation strategy.

[0088] Preferably, the robot chassis wear detection in step S2 is specifically as follows:

[0089] Identify high-impact areas based on topographic pothole data;

[0090] In this embodiment, the terrain pothole data obtained in the previous step is time-series processed, and the interaction path between the robot and the pothole area is reconstructed in chronological order. By extracting the impact acceleration value of each robot at a specific terrain pothole unit (each unit has a side length of 1 meter), the instantaneous acceleration value recorded by the three-axis acceleration sensor is used as the evaluation parameter. The acceleration threshold is set to ±3g, which comes from the correlation statistics of impact damage and acceleration values ​​in multiple field operation experiments. Acceleration events greater than the threshold are screened, and the corresponding pothole coordinates and impact time are recorded. Combined with the correspondence between the robot path and the terrain, the area where two or more impact events occur within a unit time (1s) is identified as a high impact area. To ensure accuracy, the sensor data needs to be low-pass filtered at 10Hz to remove high-frequency noise. The Butterworth filter design is adopted, the filter order is 2nd order, and the cutoff frequency is set to 10Hz.

[0091] Determine the chassis stress condition based on the high impact area and obtain chassis stress data;

[0092] In this embodiment, within the identified high-impact area, six strain gauge arrays (each array contains four strain gauges, arranged in a full-bridge measurement format) are deployed on the robot chassis to collect strain data in real time. The sampling frequency is set to 1000 Hz. The collected strain data is converted into actual force values ​​using the material mechanics formula: σ = E·ε, where the unit of stress σ is MPa and E is the elastic modulus of the material (taking Q345 steel as an example, the value is 2.1×10 11 Pa), where ε is the dimensionless strain value. The force direction is determined by the orientation of the strain gauges at each measuring point (0°, 45°, and 90°), ultimately forming a three-dimensional vector-based chassis force data matrix. This data matrix is ​​matched with the pothole impact timestamp to obtain the actual force distribution during the high-impact event.

[0093] Identify chassis stress concentration areas based on chassis force data;

[0094] In this embodiment, the maximum stress value of each measuring point in the force data matrix per unit time is used as input to construct a two-dimensional grid model of the chassis structure (unit size is 10cm×10cm), and the stress value is mapped to the corresponding unit of the model. The stress concentration identification threshold is set to 200MPa. This threshold is derived from 58% of the yield strength of chassis steel (Q345 yield strength is 345MPa). When the stress values ​​of three or more consecutive adjacent grid units exceed the threshold, they are marked as stress concentration areas. After the marking is completed, the boundary range is expanded by 5cm through the expansion algorithm to include the edge force-affected area, and finally the stress concentration area position information and unit number list are output.

[0095] Calculate surface roughness based on chassis stress concentration areas;

[0096] In this embodiment, the central unit (with the greatest force) in the stress concentration area is selected as the measurement reference point, and a contact surface profiler (model Mitutoyo SJ-410) is used to measure the roughness in the area. The spacing between each measuring point is 1 mm, the measuring length is set to 5 cm, and the roughness calculation parameters are Ra (arithmetic mean deviation) and Rz (maximum height). The Ra calculation formula is: Ra = (1 / L) ∫ 0^L | Z (x) | dx, where Z (x) is the distance of the profile curve from the reference line, and L is the sampling length. The results of all measuring points are reconstructed into a two-dimensional roughness distribution map according to the spatial position. The area with Ra value greater than 1.6 μm or Rz value greater than 6.3 μm is marked as a rough surface area. After the data is linearly interpolated and reconstructed, the roughness image and key numerical data are output.

[0097] Evaluate the fatigue state of the main load-bearing base plate based on the surface roughness to obtain fatigue data of the main load-bearing base plate;

[0098] In this embodiment, the Ra value and Rz value of the rough surface area are substituted into the empirical formula of steel fatigue life: N = A (1 / Ra)^b, where N is the expected fatigue life (in terms of the number of impacts), A and b are empirical constants, and the values ​​for Q345 steel are 1.5×10 7 Similar to 1.2. Combined with the impact count data for the same location in the historical operation path, the fatigue remaining life of each location is estimated. If the remaining impact count N is less than 5000, it is recorded as a potential fatigue area. The location coordinates and remaining fatigue life value of this area are output to generate the fatigue data table of the main load-bearing base plate.

[0099] Collecting the main load-bearing base plate image according to the main load-bearing base plate fatigue data;

[0100] In this embodiment, after locating the coordinates of the fatigue zone, the visual platform is controlled to move to the specified position by the robotic arm control system, and an industrial camera (model Basler acA1920-40gc, resolution 1920×1200) is used for image acquisition. The light source uses a ring LED fill light, and the illumination angle is set to 45° to avoid reflections. The camera focal length is set to 12mm, the acquisition distance is 30cm, the exposure time is set to 0.5ms, the image acquisition frame rate is 5fps, and 5 frames are taken continuously to ensure image integrity. The image is saved in TIFF format, and grayscale equalization processing is performed after acquisition to enhance the surface texture features and provide a clear image basis for subsequent erosion detection.

[0101] Perform erosion detection based on the main bearing base plate image to obtain erosion data;

[0102] In this embodiment, image preprocessing is performed on the main support baseplate image. This image is acquired using an industrial-grade camera mounted on the bottom of the work robot. The image acquisition resolution should be no less than 2048×2048 pixels, and the image acquisition frequency is set to once every 30 seconds. During the image preprocessing step, a 5×5 median filter is first used to denoise the image to remove high-frequency interference caused by on-site dust, water stains, or light reflections. A contrast enhancement method based on histogram equalization is then used to redistribute the image's grayscale range, improving the clarity of grayscale boundaries in eroded areas and thereby enhancing subsequent edge recognition accuracy. After preprocessing, the eroded area identification process begins. Specifically, a region segmentation method based on image grayscale value distribution differences is used. By sliding an analysis window across the image, the local grayscale mean and standard deviation within each window are extracted and compared with the grayscale statistical characteristics of known normal areas. When the local grayscale mean is at least 30% lower than the normal area mean and the standard deviation is greater than a set threshold (e.g., 50), the area is identified as a suspected eroded area. In order to further verify the regional characteristics, the three-dimensional morphology of the suspected erosion area was collected by structured light projection. The depth of the tiny surface depressions was reconstructed by a structured light scanning system composed of a camera and a projector. When the illumination angle was controlled to 30 degrees and the viewing angle was kept orthogonal, the acquired image was converted into a depth map through a grayscale mapping relationship, and the average depth and maximum depth of the suspected erosion area were extracted in micrometers. Subsequently, the boundary contour, center position, area (in square millimeters), average depth, and maximum depth information of each identified erosion area were sorted to form a unified data structure and stored as an erosion data record. The data record includes fields such as the erosion area number, two-dimensional image coordinates, corresponding three-dimensional spatial position, area, depth parameters, etc., and is uploaded to the digital twin control platform for further structural wear assessment.

[0103] The wear degree of the robot chassis is determined based on the erosion data to obtain chassis wear data.

[0104] In this embodiment, the erosion area data is spatially aligned with the CAD model coordinates of the chassis. In the specific operation, a one-to-one mapping relationship between the image coordinate system and the model coordinate system is established through multiple fixed identification points on the main load-bearing base plate, and the coordinate unification is completed using the standard rigid body alignment method. After the alignment is completed, the area, depth and distribution range of each erosion area are statistically analyzed, and classified according to their positional relationship with the main load-bearing area of ​​the chassis. If the erosion area falls on the central support beam of the chassis, the load transfer contact area or the suspension connection part, and the area exceeds 50 square millimeters and the depth exceeds 80 microns, it is determined to be a key wear area; if the erosion area is located in the edge package plate or non-structural connection area, even if the area reaches 100 square millimeters and the depth exceeds 100 microns, it is not temporarily classified as severe wear, but is regarded as non-critical wear. After counting all the erosion areas, the cumulative erosion area and total erosion volume of the entire platform base plate are further calculated, and the distribution density of the erosion area on the chassis plane is analyzed. When the distribution of erosion areas covers more than 10% of the chassis surface area, or there are more than two areas with erosion depth greater than 120 microns and an area greater than 150 square millimeters in the high-intensity area, the chassis wear level is judged to be a severe wear state. Finally, the relative position, size, depth and cumulative data of various erosion areas in the chassis structure are combined to form a complete chassis wear assessment result. The result includes the wear level (divided into mild, moderate and severe), the coordinate distribution of the wear area, the total erosion area, the erosion volume, the wear distribution of key areas, and is output in a structured manner for the synchronous update of the virtual model in the digital twin platform and the input basis for the subsequent adjustment of the multi-robot operation path.

[0105] Preferably, step S3 is specifically as follows:

[0106] Step S31: planning ultrasonic sensor placement based on chassis wear data to obtain ultrasonic placement location information;

[0107] In this embodiment, all severely worn areas in the chassis wear data are located and numbered, and the wear areas are mapped using the chassis 3D structure data embedded in the digital twin platform. The principle adopted for point planning is to cover the geometric center of the severely worn area and its edge extension area, where the wear area is greater than 100mm 2And areas where the erosion depth exceeds 120μm need to be equipped with no less than 3 ultrasonic sensor monitoring points. All layout points must meet the lower limit of the layout density of no less than 30mm between sensors, and be laid outward in a spiral expansion manner to identify potential crack propagation paths. The sensor type is a piezoelectric ultrasonic probe with a center frequency of 5MHz, a beam width of 10°, a minimum detection depth of 0.5mm, and a maximum detection depth of 20mm. The sensor and the chassis surface are in contact with a coupling agent for acoustic resistance reduction. The final layout location information includes three-dimensional coordinates (in millimeters), sensor number, action area number, installation direction angle and other fields to form the ultrasonic layout location information.

[0108] Step S32: transmitting ultrasonic signals according to the ultrasonic point location information to obtain original ultrasonic echo signal data;

[0109] In this embodiment, when transmitting ultrasonic signals based on ultrasonic point location information, the ultrasonic transmitter sequentially drives the piezoelectric ultrasonic sensors at each location through a multi-channel switching matrix, triggering pulses at a frequency of 1kHz at each sensor. The single transmitted signal is a pulse wave with a center frequency of 5MHz, a pulse width of 2μs, and a voltage amplitude of 80V. After the signal enters the chassis material through the coupling agent, the reflected wave is recovered in the same sensor channel. The sampling frequency is set to 100MHz, the sampling length is 2048 points, and the acquisition time window for each channel is set to 20μs, corresponding to the detection of echo characteristics within a depth of 20mm. The collected raw ultrasonic echo signal data is stored in the time domain. Each location corresponds to a set of one-dimensional time series signals, including a complete curve of the signal amplitude changing over time. All signal data is recorded in a JSON structure and packaged and uploaded to the analysis module as the basis for subsequent echo feature recognition.

[0110] Step S33: identifying acoustic wave dual-path reflections based on the original ultrasonic echo signal data to obtain acoustic wave dual-path reflection data;

[0111] In this embodiment, when identifying dual-path acoustic reflections based on raw ultrasonic echo signal data, waveform normalization is performed on each set of ultrasonic echo signals, normalizing the maximum signal amplitude to 1, and performing zero-mean correction for baseline drift. Subsequently, envelope detection is used to extract the peak positions of each primary and secondary echo. The presence of multipath reflection is determined by calculating the time delay intervals between different echoes. The judgment criteria are: the presence of a secondary peak with an amplitude exceeding 30% of the primary peak after the primary reflection peak, and the time delay between this secondary peak and the primary peak is less than 10 μs, indicating the presence of dual-path acoustic reflections. To further verify path overlap, short-time Fourier transform (SFT) is used to perform time-frequency analysis on the signals. The frequency domain variation consistency index is set to 90%, meaning that a true dual-path signal is identified when the frequency energy distribution overlaps between the primary and secondary peaks by 90%. For each identified dual-path echo pair, information such as the first echo time, second echo time, amplitude ratio, and frequency matching is recorded to form a dual-path acoustic reflection data table.

[0112] Step S34: Calculating the crack size based on the acoustic wave dual-path reflection data; determining the degree of ultrasonic crack damage based on the crack size to obtain chassis crack data;

[0113] In this embodiment, when calculating the crack size based on the dual-path reflection data of the acoustic wave, the known sound velocity parameters (the longitudinal wave sound velocity in steel is 5900m / s) are combined with the time difference between the main echo and the secondary echo to calculate the depth dimension of the crack in the direction of ultrasonic propagation. The crack length is trigonometrically converted by the beam angle and the sound path length, and the maximum energy area range within the approximate near-field divergence angle is taken as the lateral dimension of the crack. If the time difference is 2μs, the path difference is 11.8mm, and the corresponding crack depth dimension is 5.9mm. According to the crack size classification standard: a depth greater than 5mm and a lateral dimension exceeding 8mm is defined as a medium crack, and a depth exceeding 10mm and a lateral dimension exceeding 12mm is defined as a severe crack. The crack data includes crack number, position coordinates, depth dimension, lateral dimension, and crack grade fields. All data are summarized according to the structure number to generate a chassis crack data table, which is imported into the digital twin system for synchronous annotation of three-dimensional structures.

[0114] Step S35: performing collaborative assembly according to the chassis crack data to obtain collaborative assembly data;

[0115] In this embodiment, when performing collaborative assembly operations based on chassis crack data, the multi-robot collaborative assembly control module in the digital twin platform is first called to locate the crack area determined on the model using the crack data, and generate the assembly path for the chassis maintenance parts or reinforcements. The assembly task is assigned different robot tool head types based on the crack level. For example, for severe crack areas with a depth exceeding 10mm, a robot equipped with an electric torque gun and a riveting tool head is assigned. The assembly holes are preset within a 10mm area extending from both ends of the crack. Three fixed points are required on each reinforcement plate, with an assembly spacing of 20mm. The assembly direction must meet the requirement of an angle of no more than 10° with the structural surface. Path planning uses five-axis motion instructions. The robot number, task ID, path coordinates, motion instructions, assembly tool ID, and other contents form collaborative assembly data in the form of a task package and are uploaded to the collaborative control system scheduling module.

[0116] It is particularly important that step S35 includes the following steps:

[0117] Step S351: Identify crack-affected components based on chassis crack data to obtain crack-affected component data;

[0118] In this embodiment, crack identification uses a method based on grayscale image segmentation and region growing algorithm to process the chassis three-dimensional point cloud crack data. First, the input crack data comes from the laser point cloud data collected by the chassis three-dimensional scanning equipment. The point cloud density is set to 100 sampling points per square millimeter, and the data accuracy is 0.01mm. The edge information of the crack area is enhanced by morphological expansion and corrosion operations, and the continuous crack boundary is extracted. Then, the crack boundary is mapped to the assembly BOM structure diagram, and the geometric intersection area of ​​the crack area and the structural component is calculated by the degree of spatial coordinate overlap. If the intersection area exceeds 10% of the total area of ​​the component, it is marked as a component affected by the crack. The intersection area threshold is set here to 10%, which is set based on the assembly mechanics reliability experiment. The final crack-affected component data includes the component number, the three-dimensional coordinates of the crack interference area, the corresponding structural hierarchy path and the interference area ratio.

[0119] Step S352: performing assembly path adjustment analysis based on the crack-affected component data to obtain assembly path optimization data;

[0120] In this embodiment, the assembly path analysis adopts the discrete assembly path grid analysis method. The assembly path where the affected components are located is spatially discretized in units of 0.1mm to generate a three-dimensional path point set, and the posture change matrix of each path segment, the distance to the surrounding components, and the path pass width are recorded. The assembly path width uses 5mm as the minimum pass threshold, and the posture change is calculated using Euler angles. The path risk points are judged by comparing the spatial conflict probability of the path segments in the original assembly path with the contact rate of the cracked components. All path segments are marked as "safe", "adjustable" or "high risk", and the assembly path sequence is reconstructed according to the risk labels. The path adjustment rule takes the minimum change path length and the minimum change posture as the optimization goals, and adopts a three-dimensional space search algorithm to output a new path segment sequence with non-overlapping and high pass rate, thereby forming assembly path optimization data.

[0121] Step S353: Perform multi-robot task scheduling planning based on the assembly path optimization data to obtain robot task scheduling data;

[0122] In this embodiment, a multi-robot parallel assembly task partitioning algorithm is adopted, and each optimized assembly path segment is divided into task modules according to the operation unit. The task module is divided into 20 action instructions as a basic unit, and each action instruction includes posture data, fixture control instructions, movement path and process waiting time. The execution time of each module is estimated, and the execution time error threshold is set to 0.1s. A DAG (directed acyclic graph) is used to construct the dependency relationship diagram between modules, and each edge in the diagram represents the sequence between modules. Then, a task scheduler based on genetic algorithm is used to match robot resources. The scheduling objectives are the minimum total task execution time and the maximum resource balance rate. The number of genetic algorithm iterations is set to 200 times, the initial population size is 40, the crossover rate is 0.7, and the mutation rate is 0.05. The task scheduling data finally generated includes each robot number, the corresponding task module number, the start time, the end time and the execution path.

[0123] Step S354: Perform dynamic collaborative execution simulation control according to the robot task scheduling data to obtain collaborative assembly data.

[0124] In this embodiment, collaborative simulation control is implemented through a digital twin simulation platform. The platform used is based on Unity3D and ROS (Robot Operating System) for integrated control. A robot structure model that is completely consistent with the physical system is loaded in the simulation system, and task scheduling data is input as a drive instruction. Each robot control unit receives an action instruction and decomposes the execution instruction through a motion controller, which is converted into a sequence of joint angle changes for servo control. The angle, position, speed and load status of each joint of the robot are updated in real time in the simulation platform. In order to verify the conflicts and waiting events in collaborative control, an event detector is set to monitor the spatial overlap between task modules. If the spatial overlap is lower than the set threshold of 1.5mm, the task pause signal is triggered and the path priority is automatically adjusted. The collaborative assembly data finally output includes a robot synchronous execution schedule, an action completion mark, a spatial interference count and an assembly action execution trajectory record.

[0125] Step S36: Perform assembly interference detection based on the collaborative assembly data to obtain assembly interference data.

[0126] In this embodiment, collaborative assembly data is called as input, which includes the collaborative robot number, the sequence of motion path points of each robot (three-dimensional coordinates), the geometric model information of the assembly parts, the type and shape and size parameters of the assembly tools, and the operation sequence of each robot at different time steps. Secondly, relying on the three-dimensional structural model of the chassis built on the digital twin platform, the corresponding static structural data and the CAD geometric models of the structural reinforcement parts, screw holes, brackets, welds, etc. in the assembly area are called in, and all of them are loaded into the virtual simulation environment. Assembly interference detection adopts a fine three-dimensional collision detection process. The detection tool uses a method based on spatial bounding volumes. First, an accurate directional bounding box (OBB) is constructed for each assembly tool, assembly part and robot end effector. The size parameters of each bounding box are strictly defined according to the physical structure. For example, the bounding box size of the riveting gun is 120mm long, 40mm wide and 60mm high. A separate bounding volume needs to be constructed for each key frame on the motion path of the end effector. Secondly, the motion trajectory of each step of the robot operation is interpolated with a step size of 10ms to ensure that assembly interaction detection is performed for each tiny motion step. The interference judgment criteria are as follows: if the robot tool bounding box and the chassis structure bounding box at any time have a volume intersection, and the minimum distance of the intersection volume is less than 0.5mm, it is identified as a potential assembly interference point; if the volume of the intersection area exceeds 5mm 3Or the angle between the bounding boxes is less than 30°, it is further classified as severe assembly interference. In addition, for multi-robot collaborative tasks, it is also necessary to detect whether there is a spatial conflict between different robots due to path intersection or end effector proximity. The detection method is: if the distance between the bounding boxes of two robots in the same time frame is less than 15mm, it is marked as an "assembly path interference warning point". During the interference detection process, a three-dimensional interference body model is generated for each assembly contact area, and the intersection coordinates and interference volume values ​​are extracted. All data are recorded in a structured table, including the interference number, occurrence time step, interference object name, interference type (tool-structure / tool-tool / tool-part), interference volume (unit: mm 3 ), the three-dimensional coordinates of the interference position (in mm), the normal angle of the interference contact surface (in degrees), the corresponding task number and robot ID. The final assembly interference data is used for subsequent path optimization or assembly adjustment control module calls.

[0127] Preferably, step S36 is specifically as follows:

[0128] Step S361: Identify component poses based on collaborative assembly data to obtain component pose data;

[0129] In this embodiment, the implementation process for identifying component poses based on collaborative assembly data first calls upon the assembly component number, assembly time point, robot end-effector coordinate information, and assembly target area coordinate system corresponding to each assembly task recorded in the collaborative assembly data. Using a six-degree-of-freedom coordinate transformation, the robot TCP end-effector pose data in the world coordinate system (including position (x, y, z) and posture (rotation angles θx, θy, θz about the X, Y, and Z axes)) is transformed into the assembly reference coordinate system of the assembly component. The coordinate transformation method employed is a combination of Euler angles and translation matrices. The specific formula is not detailed here. Instead, the position of each assembly node is adjusted by a known tool center point offset (determined by the assembly tool geometry; for example, the end welding gun has an axial offset of 120 mm and a radial offset of 35 mm) to ensure that the pose of the assembled component is consistent with the assembly interface orientation. For each time frame, a timestamp-containing pose sequence data is generated by recording the component's three-dimensional pose (composed of three rotation angles and three translation amounts) and stored as a component pose data file.

[0130] Step S362: constructing a component boundary model based on the component pose data;

[0131] In this embodiment, the six-degree-of-freedom pose data of the parts obtained in the previous step are read and the CAD geometric information of the parts is matched. The information includes the surface boundary point cloud (STL format or OBJ format) and the assembly direction vector of each assembly. The AABB (axis-aligned bounding box) method is used to construct the bounding box. The specific operation is to perform coordinate projection on all boundary points, extract the maximum and minimum boundary points in the three axes of X, Y, and Z, and form a hexahedral boundary body parallel to the coordinate axis. The size of the boundary body is ΔX=xmax-xmin, ΔY=ymax-ymin, ΔZ=zmax-zmin, and the three-axis length data is recorded in millimeters. In order to improve the calculation accuracy of the spatial boundary model, each boundary body is further subdivided into 50×50×50 voxel units, and the side length of each voxel is 2mm to achieve more refined boundary envelope calculation and lay the foundation for subsequent space occupancy analysis.

[0132] Step S363: performing component space occupancy analysis based on the component boundary model to obtain spatial envelope data;

[0133] In this embodiment, voxel scanning is performed on the boundary body, and each voxel is marked as "occupied" or "free" in the virtual space. Three-dimensional space Boolean volume operations are used to perform overlapping voxel statistics on multiple assemblies. If a spatial coordinate point is covered by two or more boundary bodies in the same time period, the point is recorded as a "co-occupied space point". All co-occupied points are integrated to form the spatial envelope data set at that moment, including space occupancy rate (number of occupied voxels / total number of voxels), number of voxels in the overlapping area, three-dimensional coordinate position, spatial coordinate direction distribution density and other information. The spatial envelope data is organized in a time series manner, and the file structure contains timestamp, assembly number, spatial direction (X / Y / Z), volume distribution data and mapping information relative to the reference coordinate system.

[0134] Step S364: restoring the assembly path based on the spatial envelope data;

[0135] In this embodiment, the center coordinate point of the component at each moment in the spatial envelope data is read (as the assembly position reference point), and the assembly path is restored using the third-order spline interpolation method in combination with the robot path point sequence in the collaborative assembly data. The path interpolation step size is 5mm. During the interpolation process, the path smoothness must be guaranteed (continuous second-order derivatives do not change suddenly), and a gap of at least 3mm must be maintained between the path and the non-traversable area in the spatial envelope. For assembly paths with curved surface placement tasks, the path must also be projected onto the reference surface, and the assembly direction must be corrected in real time based on the surface normal. The final restored path is output in the form of a three-dimensional path point set (including path point number, spatial coordinates, and path direction vector) for subsequent interference analysis.

[0136] Step S365: performing assembly space overlap calculation based on the assembly path and the spatial envelope data to obtain assembly interference data.

[0137] In this embodiment, the path point set is voxelized, and each path point is converted into the index number of the corresponding voxel unit in the three-dimensional voxel coordinate system of the spatial envelope. Subsequently, a Boolean intersection operation is performed on the path voxel set and the spatial envelope voxel set constructed in step S363. If the voxel number where the path point is located has been marked as "occupied", it is considered that a spatial overlap has occurred. For each spatial overlapping point, its path point number, the assembly path segment number to which it belongs, the time step, the overlapping voxel number, the coordinate position, and the Euclidean distance to the nearest boundary point are recorded. If two consecutive points on the path overlap, and the corresponding spatial overlapping volume is greater than 10mm 3 , it is marked as an assembly interference point, and interference angle data is recorded based on the direction vector of the overlapping area (calculated from the angle between the path direction and the normal of the envelope boundary surface). All overlapping information is compiled into an assembly interference data file, which includes fields such as path number, assembly ID, interference volume, number of overlapping voxels, minimum spacing, interference location coordinates, and interference direction angle. This data will serve as the basic input for assembly task scheduling and path reconstruction.

[0138] Preferably, step S4 is specifically as follows:

[0139] Step S41: calculating the interference volume based on the assembly interference data;

[0140] In this embodiment, the assembly interference data includes the position coordinates, voxel numbers, path point numbers, interference direction angles, and boundary point cloud data of the overlapping spatial area where multiple assembly paths overlap with the spatial envelope. The calculation of the interference volume is based on the three-dimensional voxel unit, and the volume synthesis is performed using the voxel accumulation method. First, the spatial coordinates of each number marked as an interference voxel are decoded, and the unit volume (8mm) of each voxel is calculated according to the side length of the voxel divided in space (for example, the side length of each voxel is 2mm). 3 ). Then, the number of interference voxels corresponding to each set of assembly paths is summarized and multiplied by the volume of a single voxel to obtain the total interference volume of the path. If the interference area is a common voxel at the intersection of multiple paths, the voxel unique counting principle is used to prevent repeated counting. In order to enhance the volume accuracy, the boundary mesh is generated for adjacent interference voxels using the octahedron merging rule, and the boundary mesh points are subdivided and reconstructed using the MarchingCubes algorithm, and the boundary volume is calculated according to the triangular mesh volume formula. Finally, the interference volume of each pair of assembly part interference pairs is output in cubic millimeters. The fields include the assembly path number, the number of interference voxels, the rough volume estimate (voxel counting method), the fine volume calculation value (mesh reconstruction method), and the corresponding timestamp.

[0141] Step S42: sorting the parts interference risks according to the interference volume to obtain the parts interference risk sequence number;

[0142] In this embodiment, by analyzing the interference volume data of each set of parts obtained in the previous step, an interference risk measurement index of the assembled parts is constructed. The risk measurement index is calculated using the following formula: R = V / D × (1 + θ / 90), where R is the interference risk coefficient and V is the interference volume (unit: mm 3 ), D is the distance from the minimum interference point to the boundary (unit: mm), θ is the angle between the path direction and the boundary normal (unit: degree), and the range of θ is limited to 0 to 90 degrees. The value of D is calculated by the minimum Euclidean distance between the assembly path and the boundary point. The smaller the distance, the closer the path is to the boundary wall, and the higher the interference risk. θ is the angle obtained by calculating the inner product of the path vector and the normal vector of the interference envelope surface. All interference risk coefficients R are sorted from large to small, and the number of each corresponding assembly part is marked to form an interference risk sequence table. The fields of the sequence table include: part number, interference path number, interference volume V, minimum spacing D, angle θ, risk coefficient R, and sort number (serial number). After sorting is completed, it is output in CSV format for subsequent assembly sequence optimization.

[0143] Step S43: Optimizing the assembly process sequence based on the part interference risk sequence number to obtain an optimized assembly process sequence;

[0144] In this embodiment, the interference risk sequence number is read, and the parts corresponding to the assembly paths with interference risk coefficients greater than 100 are defined as "high-risk assembly parts". Based on the original assembly process sequence table (including workstation number, part ID, assembly sequence, and process dependency), the high-risk parts and their dependency paths are reordered. If a part has a dependency relationship with multiple high-risk parts, its assembly sequence is adjusted to the back. The reordering uses a directed acyclic graph (DAG) to construct a process dependency graph, where nodes represent part assembly steps and edges represent dependency relationships. By performing topological sorting on the graph and combining it with the reverse priority order of interference risk, a new assembly process optimization sequence is generated. During the sorting process, the following restrictions are enforced: ① Only one part can be assembled at the same time at the same workstation; ② Dependencies cannot be violated; ③ The assembly of high-risk parts must be placed after low-risk parts. The final output optimization sequence fields include: part ID, original assembly sequence number, optimized assembly sequence number, interference risk coefficient, dependency path, and execution time point.

[0145] It is particularly important that step S43 includes the following steps:

[0146] Step S431: extracting the sequence of high-risk interference parts based on the part interference risk serial number to obtain high-risk part sequence data;

[0147] In this example, we obtain data on the serial numbers of parts with interference risks and sort them by interference volume from high to low. Based on this sequence, we use a Python script to call structured query instructions to match and extract part code information from the database. Using an ORDER BY statement in SQL combined with a JOIN query, we associate each part's unique number with its spatial location, assembly stage, and interference characteristics, and filter them based on an interference volume threshold. This threshold is set to 2.5 cm. 3 The threshold, derived from preliminary experimental statistics, represents the boundary value at which spatial overlap significantly affects the assembly process. All parts whose interference volume exceeds this threshold are classified as high-risk parts and are output as high-risk part sequence data in the order of their original serial numbers. This output data is formatted in JSON format and includes fields such as part number, interference volume value, spatial block number, and assembly area code. This data is used in subsequent assembly logic relationship screening processes.

[0148] Step S432: Screening the assembly order logic relationship based on the high-risk parts sequence data to obtain assembly logic constraint data;

[0149] In this embodiment, based on the sequence data of high-risk parts, the assembly design BOM structure tree is called for parsing to determine the superior-subordinate relationship of each part in the assembly process. The BOM structure defines the master-slave assembly connection relationship, nesting level and connection surface constraint parameters of the parts. The BOM tree node depth information (Depth) and connection direction parameter (ConnectionAxis) are used to determine the assembly order relationship between the parts. If a high-risk part A is the parent node of another high-risk part B, and its connection direction is consistent with the interference direction (i.e., the spatial vector angle is less than 30°), it is inferred that A must be assembled before B. At the same time, combined with the predefined sequence rules in the assembly process specifications, such as tightening parts after locating parts, covering parts after basic parts, etc., the screening results are supplemented. The final output assembly logic constraint data is recorded in a table form, and the fields contain part number pairs (predecessor ID, subsequent part ID), connection direction vector, sequence relationship label and impact level label, which are used for subsequent assembly sequence judgment.

[0150] Step S433: performing assembly sequence feasibility judgment based on the assembly logic constraint data to obtain assembly sequence feasibility data;

[0151] In this embodiment, assembly logic constraint data is utilized to call a topological sorting algorithm to determine the feasibility of the assembly sequence. The specific algorithm adopts a topological directed graph loop detection method (a variant of the Kahn algorithm), constructing a directed graph structure with parts as graph nodes and constraint relationships as edges, and determining whether there are loops. If a loop is detected (i.e., a conflicting assembly sequence exists), the part pairs in the loop are returned and marked as "conflict relationships" for subsequent path reconstruction to avoid them. For paths without loops, the node arrangement order is extracted as the currently feasible assembly path candidate sequence. Simultaneously, by calling the assembly reachability calculation engine, based on the three-dimensional assembly space data (obtained from the previous step), it is calculated whether the space required for part assembly in each path is blocked by parts in subsequent steps. If the blockage degree is greater than 85%, it is considered unreachable and marked as an infeasible path. Ultimately, all sequences that meet the logical relationship and are spatially reachable are marked as "feasible" and output to a feasibility data set in the format of an XML document containing path number, sequence array, feasibility label, and conflict identification fields.

[0152] Step S434: reconstructing the interference avoidance priority path according to the assembly sequence feasibility data to obtain assembly path reconstruction data;

[0153] In this embodiment, based on the set of assembly sequence paths marked as "feasible", a path reconstruction algorithm based on the interference volume penalty function is called to select the priority path. The path length, interference penalty factor and process sequence priority are defined as weight parameters in the heuristic function. The interference penalty factor is obtained by reading the interference volume data of each part and assigning a weight coefficient (set to 1.0 / cm 3 ), used for path cost evaluation. The assembly path length is calculated by the physical distance between the assembly action nodes, in mm, and is derived from the coordinate difference of the workstation positions marked in the CAD assembly layout file. The priority level is given by the assembly step priority coefficient set in the process manual. The coefficient range is [1,10], and the lower the coefficient, the higher the priority. The system calculates the total path cost value based on the three factors and selects the path with the lowest total cost from all feasible paths as the reconstructed path. The output assembly path reconstruction data is in CSV format, and the fields contain information such as the assembly order array, the total cost value of the path, the penalty distribution of each step, and the priority distribution.

[0154] Step S435: Assign assembly action sequence numbers based on the assembly path reconstruction data to obtain an assembly process optimization sequence.

[0155] In this embodiment, the assembly path reconstruction data output in CSV format is read, and an action sequence number is assigned to each assembly action in the order of the path. The number assignment is based on the assembly station coding and execution time priority strategy, and the assembly line logic is adopted, that is, the action numbers in the same assembly area must be incremented and cannot jump. The action sequence number format is "A-XXX", where A represents the station number and XXX is a three-digit sequence number. For example, B-012 represents the 12th step action of station B. The system is numbered according to the assembly action definition of each part in the station in the path. Actions such as positioning, insertion, and tightening are subdivided into different steps for numbering. A forward scanning algorithm is used in the numbering process to ensure that there are no conflicts and duplications. After all the numbering is completed, it is summarized into an assembly process optimization sequence with a JSON structure. The fields include action number, part number, station number, action type, preceding action ID, expected execution time, etc., for use by the subsequent robot scheduling module.

[0156] Step S44: performing task priority analysis according to the assembly process optimization sequence to obtain task priority data;

[0157] In this embodiment, each assembly task is defined as (part ID, workstation ID, assembly time period), and a task scheduling matrix is ​​constructed. Each cell in the matrix records the time binding status of the part and the workstation. The priority scoring rule is: P = α × R + β × T + γ × W, where P is the task priority score, R is the interference risk coefficient (derived from step S42), T is the assembly sequence of the task in the sequence (the smaller the sequence number, the higher the score), and W is the resource conflict index (if the task shares a workstation with other tasks, W = 1, otherwise W = 0). α, β, and γ are taken as 5, 2, and 3 respectively, which are used to strengthen the weight of the interference factor in the priority. After the calculation is completed, the priority score of each task is sorted from high to low and numbered as the task priority number. The task priority data contains fields: task ID, part ID, workstation ID, priority score P, sorting number, starting time point, duration, dependent task list, etc., which serve as input parameters for multi-robot task scheduling.

[0158] Step S45: Perform multi-robot collaborative assembly operation simulation according to the task priority data to obtain collaborative assembly operation data.

[0159] In this embodiment, each robot is assigned a workstation and operating time period, and the simulation environment is synchronously scheduled and controlled via a digital twin platform (e.g., a three-dimensional process workshop built using the Unity3D physics engine). The task scheduling module reads a list of priority tasks and advances them along a timeline, sending tasks with a start time less than the current simulation time point to the robot simulation controller one by one. The robot motion control sequence includes the movement path, tool activation period, and end-operation posture. All actions are transmitted to the simulation actuator via real-time interpolation. The robot path planning uses the RRT-Connect algorithm for three-dimensional collision avoidance. The dynamic path generation interval is 20ms, and each path segment has no more than 100 control points. During the simulation, each robot status is fed back in real time, including current position, task completion status, current load information, and collaborative conflict status. The resulting collaborative assembly operation data includes task execution records (task ID, start and end time, executing robot), robot path sequence (pose data for each time step), interference conflict alarm records, and task rescheduling records. All data is stored in a structured database in a time series format for subsequent execution optimization and task replanning.

[0160] Preferably, this specification also provides a multi-robot collaborative operation simulation control system based on digital twins, which is used to execute the multi-robot collaborative operation simulation control method based on digital twins as described above. The multi-robot collaborative operation simulation control system based on digital twins includes:

[0161] The warehouse collaborative operation and transportation simulation module is used to obtain robot structural data; build a 3D twin model of the robot based on the robot structural data; and perform warehouse collaborative operation and transportation simulation based on the 3D twin model of the robot to obtain collaborative operation and transportation data.

[0162] The robot chassis wear detection module is used to perform path deviation anomaly analysis based on collaborative operation transportation data to obtain path deviation data; perform terrain pothole detection based on path deviation data to obtain terrain pothole data; and perform robot chassis wear detection based on terrain pothole data to obtain chassis wear data;

[0163] The assembly interference detection module is used to perform ultrasonic crack detection based on chassis wear data to obtain chassis crack data; perform collaborative assembly based on chassis crack data to obtain collaborative assembly data; and perform assembly interference detection based on the collaborative assembly data to obtain assembly interference data.

[0164] The collaborative assembly operation simulation module is used to optimize the assembly process sequence based on assembly interference data to obtain the assembly process optimization sequence; perform task priority analysis based on the assembly process optimization sequence to obtain task priority data; and simulate multi-robot collaborative assembly operations based on the task priority data to obtain collaborative assembly operation data.

[0165] The present invention is therefore intended to be illustrative and non-restrictive in all respects, with the scope of the invention being defined by the appended claims rather than the foregoing description, and all changes that come within the meaning and range of equivalents of the application documents are intended to be embraced therein.

[0166] The foregoing description is intended only to provide specific embodiments of the present invention, which will enable those skilled in the art to understand and implement the present 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 present invention. Therefore, the present invention is not intended to be limited to the embodiments shown herein, but is to be construed in the widest possible manner consistent with the principles and novel features disclosed herein.

Claims

1. A multi-robot collaborative operation simulation control method based on digital twin, characterized in that: The following steps are involved: Step S1: Acquire robot structural data; construct a three-dimensional twin model of the robot based on the robot structural data; perform warehouse collaborative operation and transportation simulation based on the three-dimensional twin model of the robot to obtain collaborative operation and transportation data; Step S2: performing path deviation anomaly analysis based on the collaborative operation transport data to obtain path deviation data; performing terrain pothole detection based on the path deviation data to obtain terrain pothole data; performing robot chassis wear detection based on the terrain pothole data to obtain chassis wear data; Step S3: performing ultrasonic crack detection based on the chassis wear data to obtain chassis crack data; performing collaborative assembly based on the chassis crack data to obtain collaborative assembly data; Perform assembly interference detection based on collaborative assembly data to obtain assembly interference data; Step S4: Optimizing the assembly process sequence based on the assembly interference data to obtain an optimized assembly process sequence; Task priority analysis is performed based on the assembly process optimization sequence to obtain task priority data; multi-robot collaborative assembly operation simulation is performed based on the task priority data to obtain collaborative assembly operation data.

2. The multi-robot collaborative operation simulation control method based on digital twin according to claim 1 is characterized in that: Step S1 is specifically as follows: Step S11: Obtain robot structure data; Step S12: performing component space modeling based on the robot structure data to obtain the robot's three-dimensional component structure; Step S13: reconstructing the component assembly relationship according to the robot's three-dimensional component structure to obtain component assembly relationship data; Step S14: constructing a three-dimensional twin model of the robot based on the component assembly relationship data and the three-dimensional component structure of the robot to obtain a three-dimensional twin model of the robot; Step S15: Perform warehouse collaborative operation and transportation simulation based on the robot's three-dimensional twin model to obtain collaborative operation and transportation data.

3. The multi-robot collaborative operation simulation control method based on digital twin according to claim 2 is characterized in that: Step S15 is specifically as follows: Step S151: uploading the robot 3D twin model to the robot collaborative operation simulation platform; Step S152: Set the robot's maximum speed to 0.5m / s-2.0m / s and the acceleration range to 0.1m / s 2 -1.0m / s 2 , minimum turning radius 0.3m-1.0m; Step S153: Setting the terrain slope to 0°-5°, the ground friction coefficient to 0.2-0.8, and the load variation to 5kg-20kg; Step S154: Set the scheduling instruction delay to 10ms-200ms and the path planning refresh period to 1s-10s; Step S155: Run the robot behavior control module and the path interaction module in the collaborative operation simulation platform to obtain collaborative operation transportation data.

4. The multi-robot collaborative operation simulation control method based on digital twin according to claim 1 is characterized in that: The path deviation anomaly analysis in step S2 is specifically as follows: Extract robot operation trajectory based on collaborative operation transportation data; Perform time alignment processing on the trajectory point sequence based on the robot's running trajectory to obtain time series standard trajectory data; Calculate trajectory curvature based on time-series standard trajectory data; Identify track curves according to track curvature and obtain track curve data; Performing trajectory deviation detection on the trajectory curve data according to preset standard trajectory curve data to obtain trajectory deviation data; Calculating the lateral offset amplitude based on the trajectory offset data; Calculating the longitudinal offset amplitude based on the trajectory offset data; Calculate the lateral deviation volatility based on the lateral deviation amplitude; Perform inertial offset detection based on the longitudinal offset amplitude to obtain inertial offset abnormal data; The path deviation data are obtained by integrating the lateral deviation fluctuation rate and the inertial deviation anomaly data.

5. The multi-robot collaborative operation simulation control method based on digital twin according to claim 1 is characterized in that: The specific steps of terrain pothole detection in step S2 are: Identify the deviation mutation point based on the path deviation data; Calculate the terrain slope based on the offset mutation point to obtain slope data; Calculate terrain height difference based on slope data; Identify suspected pothole areas based on terrain height differences; Collect robot posture angle data based on suspected pothole areas; Evaluate terrain impact intensity based on robot posture angle data; The suspected pothole areas are screened according to the terrain impact intensity to obtain the terrain pothole data.

6. The multi-robot collaborative operation simulation control method based on digital twin according to claim 1 is characterized in that: The robot chassis wear detection in step S2 is specifically as follows: Identify high-impact areas based on topographic pothole data; Determine the chassis stress condition based on the high impact area and obtain chassis stress data; Identify chassis stress concentration areas based on chassis force data; Calculate surface roughness based on chassis stress concentration areas; Evaluate the fatigue state of the main load-bearing base plate based on the surface roughness to obtain fatigue data of the main load-bearing base plate; Collecting the main load-bearing base plate image according to the main load-bearing base plate fatigue data; Perform erosion detection based on the main bearing base plate image to obtain erosion data; The wear degree of the robot chassis is determined based on the erosion data to obtain chassis wear data.

7. The multi-robot collaborative operation simulation control method based on digital twin according to claim 1 is characterized in that: Step S3 is specifically as follows: Step S31: planning ultrasonic sensor placement based on chassis wear data to obtain ultrasonic placement location information; Step S32: transmitting ultrasonic signals according to the ultrasonic point location information to obtain original ultrasonic echo signal data; Step S33: identifying acoustic wave dual-path reflections based on the original ultrasonic echo signal data to obtain acoustic wave dual-path reflection data; Step S34: Calculating the crack size based on the acoustic wave dual-path reflection data; determining the degree of ultrasonic crack damage based on the crack size to obtain chassis crack data; Step S35: performing collaborative assembly according to the chassis crack data to obtain collaborative assembly data; Step S36: Perform assembly interference detection based on the collaborative assembly data to obtain assembly interference data.

8. The multi-robot collaborative operation simulation control method based on digital twin according to claim 7 is characterized in that: Step S36 is specifically as follows: Step S361: Identify component poses based on collaborative assembly data to obtain component pose data; Step S362: constructing a component boundary model based on the component pose data; Step S363: performing component space occupancy analysis based on the component boundary model to obtain spatial envelope data; Step S364: restoring the assembly path based on the spatial envelope data; Step S365: performing assembly space overlap calculation based on the assembly path and the spatial envelope data to obtain assembly interference data.

9. The multi-robot collaborative operation simulation control method based on digital twin according to claim 1 is characterized in that: Step S4 is specifically as follows: Step S41: Calculating the interference volume based on the assembly interference data; Step S42: sorting the parts interference risks according to the interference volume to obtain the parts interference risk sequence number; Step S43: Optimizing the assembly process sequence based on the part interference risk sequence number to obtain an optimized assembly process sequence; Step S44: performing task priority analysis according to the assembly process optimization sequence to obtain task priority data; Step S45: Perform multi-robot collaborative assembly operation simulation according to the task priority data to obtain collaborative assembly operation data.

10. A multi-robot collaborative operation simulation control system based on digital twin, characterized in that: For executing the multi-robot collaborative operation simulation control method based on digital twin according to claim 1, the multi-robot collaborative operation simulation control system based on digital twin comprises: The warehouse collaborative operation and transportation simulation module is used to obtain robot structural data; build a 3D twin model of the robot based on the robot structural data; and perform warehouse collaborative operation and transportation simulation based on the 3D twin model of the robot to obtain collaborative operation and transportation data. The robot chassis wear detection module is used to perform path deviation anomaly analysis based on collaborative operation transportation data to obtain path deviation data; perform terrain pothole detection based on path deviation data to obtain terrain pothole data; and perform robot chassis wear detection based on terrain pothole data to obtain chassis wear data; The assembly interference detection module is used to perform ultrasonic crack detection based on chassis wear data to obtain chassis crack data; perform collaborative assembly based on chassis crack data to obtain collaborative assembly data; and perform assembly interference detection based on the collaborative assembly data to obtain assembly interference data. The collaborative assembly operation simulation module is used to optimize the assembly process sequence based on assembly interference data to obtain the assembly process optimization sequence; perform task priority analysis based on the assembly process optimization sequence to obtain task priority data; and simulate multi-robot collaborative assembly operations based on the task priority data to obtain collaborative assembly operation data.

Images